[CFDdummy]

Greetings everyone!

I'm learning CFD currently and very new to it. As of the moment, I'm working on quasi - one dimensional nozzle flows in Chapter 7 of CFD by J. Anderson.

So I was wondering if anyone has a MATLAB code for this? I would truly appreciate if anyone could help me with this.

Thanks so much!

[Gnak]

Hey Jedidiah

I will give you something just to get it started, this is an incomplete code for the exact problem that you are referring to. I am not sure if this was even runnable to be frank. I refraining from giving you a complete code since it is not fun to do so.

Written in python btw

import numpy as np

L=3
N=30
gamma=1.4
x_max=3.1
t_max=5.1
dx=L/N
dt=0.1
C=0.5

x=np.arange(0,x_max,dx)
t=np.arange(0,t_max,dt)

s=len(x)
r=len(t)

rho=np.ones(s)
T=np.ones(s)

A=1+2.2*(x-1.5)**2

rho=np.piecewise(x,[x <= 0.5, (x > 0.5)&(x <= 1.5),],[1 , lambda x: 1.-0.366*(x-0.5), lambda x: 0.634-0.3879*(x-1.5)])
T=np.piecewise(x,[x <= 0.5, (x > 0.5)&(x <= 1.5),],[1 ,lambda x: 1.-0.167*(x-0.5), lambda x: 0.833-0.3507*(x-1.5)])
V=0.59/rho/A

t=np.ones(s)
U1=np.ones(s)
U2=np.ones(s)
U3=np.ones(s)
F1=np.ones(s)
F2=np.ones(s)
F3=np.ones(s)
J2=np.ones(s)
U1_t=np.ones(s)
U2_t=np.ones(s)
U3_t=np.ones(s)
U1_bar=np.ones(s)
U2_bar=np.ones(s)
U3_bar=np.ones(s)
rho_bar=np.ones(s)
T_bar=np.ones(s)
V_bar=np.ones(s)
F1_bar=np.ones(s)
F2_bar=np.ones(s)
F3_bar=np.ones(s)
U1_t_bar=np.ones(s)
U2_t_bar=np.ones(s)
U3_t_bar=np.ones(s)
U1_t_avg=np.ones(s)
U2_t_avg=np.ones(s)
U3_t_avg=np.ones(s)
U1r=np.ones(s)
U2r=np.ones(s)
U3r=np.ones(s)

for n in range(0,1): # time

for i in range(1,s-1):
U1[i]=rho[i]*A[i]
U2[i]=rho[i]*A[i]*V[i]
U3[i]=rho[i]*(T[i]/(gamma-1)+gamma*V[i]**2/2)*A[i]
F1[i]=U2[i]
F2[i]=U2[i]**2/U1[i]+(gamma-1)/gamma*(U3[i]-gamma*U2[i]**2/2/U1[i])
F3[i]=gamma*U2[i]*U3[i]/U1[i]-gamma*(gamma-1)/2*U2[i]**3*U1[i]**2
J2[i]=rho[i]*T[i]/gamma/dx*(A[i+1]-A[i])
U1_t[i]=(F1[i]-F1[i+1])/dx
U2_t[i]=(F2[i]-F2[i+1])/dx+J2[i]
U3_t[i]=(F3[i]-F3[i+1])/dx
t[i]=C*dx/(np.sqrt(T[i])+V[i])
U1_bar[i]=U1[i]+U1_t[i]*min(t)
U2_bar[i]=U2[i]+U2_t[i]*min(t)
U3_bar[i]=U3[i]+U3_t[i]*min(t)
rho_bar[i]=U1_bar[i]/A[i]
T_bar[i]=(gamma-1)*((U3_bar[i]/U1_bar[i])-gamma/2*(U2_bar[i]/U1_bar[i])**2)
F1_bar[i]=U2_bar[i]
F2_bar[i]=U2_bar[i]**2/U1_bar[i]+(gamma-1)/gamma*(U3_bar[i]-gamma*U2_bar[i]**2/2/U1_bar[i])
F3_bar[i]=gamma*U2_bar[i]*U3_bar[i]/U1_bar[i]-gamma*(gamma-1)/2*U2_bar[i]**3*U1_bar[i]**2
U1_t_bar[i]=(F1_bar[i-1]-F1_bar[i])/dx
U2_t_bar[i]=(F2_bar[i-1]-F2_bar[i])/dx
U3_t_bar[i]=(F3_bar[i-1]-F3_bar[i])/dx
U1_t_avg[i]=0.5*(U1_t[i]+U1_t_bar[i])
U2_t_avg[i]=0.5*(U2_t[i]+U2_t_bar[i])
U3_t_avg[i]=0.5*(U3_t[i]+U3_t_bar[i])
U1r[i]=U1[i]+U1_t_avg[i]*min(t)
U2r[i]=U2[i]+U2_t_avg[i]*min(t)
U3r[i]=U3[i]+U3_t_avg[i]*min(t)
rho[i]=U1r[i]/A[i]
V[i]=U2r[i]/U1r[i]
T[i]=(gamma-1)*(U3r[i]/U1r[i]-gamma*V[i]**2/2)
U1[0]=2*U1[1]-U1[2]
U2[0]=2*U2[1]-U2[2]
U3[0]=2*U3[1]-U3[2]

cheers,

[CFDdummy]

Hi Gnak

Thanks a lot for the response. I truly appreciate it.
I've already managed to figure out the matlab code for that problem.
I was able to show the properties along the nozzle using matlab.
Now, I'm working on the plotting the shape of CD nozzle.

Thanks again!

[pritam.gole]

If length in the code is considered as 3, then x' = x/L can't go from 0-3. It should vary from 0 to 1. But, in all the code it is considered as 0 to 3. Please help me out to understand it.

[FMDenaro]

Originally Posted by pritam.gole

If length in the code is considered as 3, then x' = x/L can't go from 0-3. It should vary from 0 to 1. But, in all the code it is considered as 0 to 3. Please help me out to understand it.

In general you are right, if the characteristic lenght coincide with the total extension of the spatial domain the non-dimensional position is from 0 to 1. But you have to check from some papers, not from a code what is referred as L.

[raviteja gonnabathula]

Originally Posted by CFDdummy

Hi Gnak

Thanks a lot for the response. I truly appreciate it.
I've already managed to figure out the MatLab code for that problem.
I was able to show the properties along the nozzle using Matlab.
Now, I'm working on plotting the shape of the CD nozzle.

Thanks again!

Can you please share the code here..