How to check my iteration results

I'm working on an interaction project for my propulsion class. I am to do an analysis on a quasi-one-dimensional air flow in a duct that is 1 meter long.

Length 1 meter
Isentropic flow
no work interaction
$$P_0=101325$$Pa
$$T_0=288$$K
Mach=0.5
cross-section area at inlet=0.1 m^2
cross-section area at exit=0.4345 m^2

I'm to perform 100,000 iterations from the inlet to the exit to calculate the exit conditions

P, T, velocity, Mach, P total, T total.

We are to use differential analysis. I wrote MATLAB code to perform this analysis but I'm also required to check my results against analytical results.

Could someone explain to me how to calculate the analytical results?

The exit conditions as calculated from my program are as follow below. Attached is my code.

Exit pressure = 101325.9571 Pa
Exit temperature = 288.0015 K
Exit velocity = 168.8889 m/s
Exit Mach = 0.49997
Exit total pressure = 191802.8588 Pa
Exit total temperature = 86.3911 K
Exit density = 0.28613 kg/m^3
P/P*= 0.99999
rho/rho*= 4.3448
T/T*= 0.99999
U/U*= 1.0001
Pt/Pt*= 0.62665
Tt/Tt*= 3.5004
Attached Files
 aeromidterm2.m (1.7 KB, 7 views)
 Recognitions: Science Advisor I suppose the first thing you would need to do is check to see if your exit area is smaller than the critical area. If not, then you should be able to just use Mach tables to find your exit conditions given your area ratio. Have you done a sanity check on your results? Your exit area is more than 4x larger, but your exit Mach is only 0.0003 different.
 I was thinking the same thing about the area and Mach. I've worked with these equations about 5 times now. I'll have to post later my steps for coming up with the equations used for the iterations.

How to check my iteration results

So I went through the algebra again and came up with a different equation for my velocity but now my Mach number increased. Something is seriously wrong. Here is my work.

As we move along the nozzle for each iteration, the area increases by dA.
steps=100,001

A variable with subscript i represents the value of the variable at the current step.
A variable with subscript j represents the value of the variable at the previous step, i-1.

$$dA=(A_{exit} - A_{inlet})/(steps-1)$$
So at any iteration step, the area can be calculated by,

$$A_i=A_j+dA$$

Given:
$$\frac{d\rho}{\rho}+\frac{dU}{U}+\frac{dA}{A}=0$$

But $$\frac{d\rho}{\rho}=\frac{\rho_i-\rho_j}{\rho_j}$$, like wise for the other terms.

$$\frac{\rho_i-\rho_j}{\rho_j}+\frac{U_i-U_j}{U_j}+\frac{A_i-A_j}{A_j}=0$$

Simplifying and solving for $$\frac{\rho_i}{\rho_j}$$

(1) $$\frac{\rho_i}{\rho_j}=3-\frac{U_i}{U_j}-\frac{A_i}{A_j}$$

Given:
$$\frac{dP}{\rho}+UdU=\frac{-\tau_w*C*dx}{\rho{A}}+\eta*dw$$

(*) $$\frac{P_i-P_j}{P_j}+U_j(U_i-U_j)=\frac{-\tau_w*C_i*dx}{\rho_j*A_j}+\eta*dw$$

$$dS*\tau_w=1/2*C_f*\rho_j*U_j^2*C_i*dx$$

$$dS=C*dx$$

So then $$\tau_w=1/2*C_f*\rho_j*U_j^2$$

Plugging into (*) and simplifying.

$$\frac{P_i-P_j}{P_j}+U_j(U_i-U_j)=\frac{-1/2*C_f*U_j^2*C_i*dx}{A_j}+\eta*dw$$

Solve for $$\frac{P_i}{P_j}$$

$$\frac{P_i}{P_j}=\frac{-.5*C_f*U_j^2*C_i*dx}{A_j}+\eta*dw-U_j(U_i-U_j)+1$$
I then grouped the non-iterative terms of subscript i (including $$C_i$$) into a variable called $$\beta$$.

$$\beta=\frac{-.5*C_f*U_j^2*C_i*dx}{A_j}+\eta*dw+1$$

So now

(2) $$\frac{P_i}{P_j}=\beta-U_j(U-i-U_j)$$

Given:
$$C_p*dT+UdU=d_q+d_w$$

Solving for $$T_i$$

$$T_i=\frac{d_q+d_w}{C_p}-\frac{U_j(U_i-U_j)}{C_p}+T_j$$

Letting $$\alpha=\frac{d_q+d_w}{C_p}+T_j$$

(3) $$T_i=\alpha-\frac{U_j(U_i-U_j)}{C_p}$$

Given:
$$\frac{dP}{P}=\frac{d\rho}{\rho}+\frac{dT}{T}$$

$$\frac{P_i-P_j}{P_j}=\frac{\rho_i-\rho_j}{\rho_j}+\frac{T_i-T_j}{T_j}$$

Simplifying and solving for $$\frac{P_i}{P_j}$$

(4) $$\frac{P_i}{P_j}=\frac{\rho_i}{\rho_j}+\frac{T_i}{T_j}-1$$

So, eq (2) into (4):

(5) $$\beta-U_j(U_i-U_j)=\frac{\rho_i}{\rho_j}+\frac{T_i}{T_j}-1$$

eq 3) into (5):

(6) $$\beta-U_j(U_i-U_j)=\frac{\rho_i}{\rho_j}-\frac{U_j(U_i-U_j)}{T_j*C_p}-1$$

eq (1) into (6):

$$\beta-U_j(U_i-U_j)=3-\frac{U_i}{U_j} - \frac{A_i}{A_j} + \frac{\alpha}{T_j} - \frac{U_j(U_i-U_j)}{T_j*C_p} -1$$

Gathering all $$U_i]/tex] terms: [tex]-U_j(U_i-U_j)+\frac{U_i}{U_j}+\frac{U_j(U_i-U_j)}{T_j*C_p}=3-\beta-\frac{A_i}{A_j}+\frac{\alpha}{T_j}-1$$

Factoring the $$U_i$$ term:

$$U_i\left(-U_j+\frac{1}{U_j}+\frac{U_j}{T_j*C_p}\right)=2-\beta-\frac{A_i}{A_j}+\frac{\alpha}{T_j}-U_i^2+\frac{U_i^2}{T_j*C_p}$$

$$U_i=\left(2-\beta-\frac{A_i}{A_j}+\frac{\alpha}{T_j}-U_j^2+\frac{U_j^2}{T_j*C_p}\right)\left(-U_j+\frac{1}{U_j}+\frac{U_j}{T_j*C_p}\right)^{-1}$$

To find the other values of the parameters my professor suggested to use these equations:

$$\rho_i=\frac{\rho_j*U_j*A_j}{U_i*A_i}$$

$$T_i=T_j+\frac{.5*(U_j^2-U_i^2)+dq+dshaft}{C_p}$$

$$P_i=\rho_i*R*T_i$$

R=283, and $$C_p$$=1004

$$M_i=\frac{U_i}{\sqrt{\gamma*R*T_i}}$$

It would be very much appreciated if someone could check over my steps. I've solved through these equations like 6 times and I can't figure out what I'm doing wrong. The only thing I can think of is the term $$\frac{-\tau_w*C*dx}{\rho{A}}$$ in equation (2). I thinking that maybe the denominator term should be $$\rho_i*A_i$$. Any suggestions would be greatly appreciated.