Why am I getting non real values for M.
iFunction file
function dMdx=s(x,M,gamma_b,a1Vec,b1Vec,a11Vec,b11Vec,xVec,xMin)a1=interp1(xVec,a1Vec,x);b1=interp1(xVec,b1Vec,x);a11=interp1(xVec,a11Vec,x);b11=interp1(xVec,b11Vec,x);if xVec==xMin dMdx=(1/4)*(-gamma_b*a1+sqrt((gamma_b*a1)^2-4*((gamma_b-1)*a1^2-(1+gamma_b)*a11+(((1+gamma_b)^2)/2)*b11)));elsedMdx= M*((1+((gamma_b-1)/2)*M^2)/(1-M^2))*(-a1+((1+(gamma_b*M^2))/2)*b1);endend
iScript
xVec=0:0.01:0.4; %axial combuster length (can be varied according to user input and units)
xi=xVec(1); %station 3 (combuster starting point)
x4=xVec(end);%station 4 (combuster end point)
gamma_b=1.36;%Average value of ratio of specific heat capacities during burning
gamma_c=1.4;%Average value of ratio of specific heat capacities during compression
X=((xVec-xi)/(x4-xi));a3=0.0038;%station 3 cross secion area in m2.
A=a3*(1+X.^2);%Area Profile
a=gradient(A,xVec);%dA/dx in equation (6-90)
a1Vec=a./A; %(dA/dx)/A in equation (6-90)
d2a=gradient(a,xVec);a11Vec=d2a./A;tau_b=1.5; %Tt4/Tt2 in eq. (6-91) overall total temerature rise in burner
theta=2;%emperical constant in eq. (6-91)
Tt2=2200;%Total temperature at station 2 depends on inlet and compression conditions (USER SPECIFIED)
Ttx=Tt2*(1+(tau_b-1)*((theta*X)./(1+(theta-1)*X)));%Temperature profile eq. (6-91)
b=gradient(Ttx,xVec);%dTt/dx in equation (6-90)
b1Vec=b./Ttx;%(dTt/dx)/Tt in equation (6-90)
d2b=gradient(b,xVec);b11Vec=d2b./Ttx;fun = @(xVec)((a./A)-(((1+gamma_b)/2)*(b./Ttx))); % Function of x only
[fMin, idx] = min(abs(fun(xVec)));xMin = xVec(idx);initial_M=1;[x,M]=ode45(@(x,M) s(x,M,gamma_b,a1Vec,b1Vec,a11Vec,b11Vec,xVec,xMin),[xMin xi],initial_M);%Solve ODE eq. (6-90)
Best Answer