Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Mac Cormack FDM (for 1D compressible nozzle flow)

Subject: Mac Cormack FDM (for 1D compressible nozzle flow)

From: Kenan

Date: 14 Aug, 2010 17:07:03

Message: 1 of 3

Hello,

I've been struggling with writing a code for solution of 1D nozzle flow since a few days, using Mac Cormack finite difference method.

This is what i've managed so far; as you'll realize, i'm a noob in Matlab.

%Nozzle geometry
dx=0.1;
x=0:dx:3;
A=(1+2.2*(x.^2-3*x+2.25))';
cfl=0.5; %Courant number
gama=1.4; %Specific heat ratio
N=length(x);%Number of grid points
tt=18; %NUMBER OF TIME STEPS

%Initial values
rho=(1-0.3146*x)';
T=(1-0.2314*x)';
for i=1:N
    V(i)=((0.1+x(i)*1.09)*(T(i)^0.5));
    deltat(i)=(cfl*(dx/((T(i)^0.5+V(i)))))';
end

dt=min(deltat); %Time step difference

%%%Mac Cormack method%%%

for j=1:tt
for k=2:N-1
    
%Prediction step%

rho_der(k)=(1/dx)*(-rho(k)*(V(k+1)-V(k))-rho(k)*V(k)*(log(A(k+1))...
    -log(A(k)))-V(k)*(rho(k+1)-rho(k)));

V_der(k)=(1/dx)*(-V(k)*(V(k+1)-V(k))-(1/gama)*((T(k+1)-T(k))+...
    (T(k)/rho(k))*(rho(k+1)-rho(k))));

T_der(k)=(1/dx)*(-V(k)*(T(k+1)-T(k))-(gama-1)*T(k)*((V(k+1)-V(k))...
    +(V(k)*(log(A(k+1))-log(A(k))))));

rho_bar(k)=rho(k)+rho_der(k)*dt;
V_bar(k)=V(k)+V_der(k)*dt;
T_bar(k)=T(k)+T_der(k)*dt;

%Correction step%

rho_der_bar(k)=(1/dx)*(-rho_bar(k)*(V_bar(k)-V_bar(k-1))-rho_bar(k)*V_bar(k)...
    *(log(A(k))-log(A(k-1)))-V_bar(k)*(rho_bar(k)-rho_bar(k-1)));

V_der_bar(k)=(1/dx)*(-V_bar(k)*(V_bar(k)-V_bar(k-1))-(1/gama)*...
    ((T_bar(k)-T_bar(k-1))+(T_bar(k)/rho_bar(k))*(rho_bar(k)-rho_bar(k-1))));

T_der_bar(k)=(1/dx)*(-V_bar(k)*(T_bar(k)-T_bar(k-1))-(gama-1)*T_bar(k)*((V_bar(k)-V_bar(k-1))...
    +(V_bar(k)*(log(A(k))-log(A(k-1))))));

rho(k)=rho_bar(k)+0.5*dt*(rho_der(k)+rho_der_bar(k));
V(k)=V_bar(k)+0.5*dt*(V_der(k)+V_der_bar(k));
T(k)=T_bar(k)+0.5*dt*(T_der(k)+T_der_bar(k));

%Boundary conditions%
V(1)=2*V(2)-V(3);
rho(31)=2*rho(30)-rho(29);
V(31)=2*V(30)-V(29);
T(31)=2*T(30)-T(29);

Q(:,1)=rho;
Q(:,2)=V;
Q(:,3)=T;

end

j
Q

end

******************************************
The problem is, there is something wrong with the for loop, but i dont know what is. It's supposed to feed new so-called rho,V and T values at the end of the loop to the beginning values rho_der,V_der and T_der again and so on. But if you choose 20 or bigger time steps (tt) it gives negative numbers, which shouldnt be. What can the problem be?

I've found one code written for mac-cormack method on this website, but i couldnt manage to decode it and tried a new code on my own, but it doesnt seem to work.

Any help would be appreciated. Thanks in advance.
Kenan

Subject: Mac Cormack FDM (for 1D compressible nozzle flow)

From: Kenan

Date: 15 Aug, 2010 10:35:23

Message: 2 of 3

If there is something unclear you want me to explain, please let me know it.

Subject: Mac Cormack FDM (for 1D compressible nozzle flow)

From: Kenan

Date: 18 Aug, 2010 11:53:05

Message: 3 of 3

Hello again,

Finally could make it work. Here is the code:


%Solution of quasi-one-dimensional isentropic flow through a
%diverging-converging nozzle, by the explicit Mac Cormack finite difference method.

%Nozzle geometry
dx=0.1;
x=0:dx:3;
A=(1+2.2*(x.^2-3*x+2.25))';
cfl=0.5; %Courant number
gama=1.4; %Specific heat ratio
N=length(x);%Number of grid points
tt=1400; %Number of time steps

%Initial values
rho=(1-0.3146*x)';
T=(1-0.2314*x)';
for i=1:N
    V(i)=((0.1+x(i)*1.09)*(T(i)^0.5));
    deltat(i)=(cfl*(dx/((T(i)^0.5+V(i)))))';
end
dt=min(deltat); %Time step difference

%%Mac Cormack method%%

for j=1:tt
    
    %Prediction step%
    
    for k=2:N-1
    
rho_der(k)=(1/dx)*(-rho(k)*(V(k+1)-V(k))-rho(k)*V(k)*(log(A(k+1))...
    -log(A(k)))-V(k)*(rho(k+1)-rho(k)));

V_der(k)=(1/dx)*(-V(k)*(V(k+1)-V(k))-(1/gama)*((T(k+1)-T(k))+...
    (T(k)/rho(k))*(rho(k+1)-rho(k))));

T_der(k)=(1/dx)*(-V(k)*(T(k+1)-T(k))-(gama-1)*T(k)*((V(k+1)-V(k))...
    +(V(k)*(log(A(k+1))-log(A(k))))));

rho_bar(k)=rho(k)+rho_der(k)*dt;
V_bar(k)=V(k)+V_der(k)*dt;
T_bar(k)=T(k)+T_der(k)*dt;

    end

%Intermediate boundary conditions
V_bar(1)=2*V_bar(2)-V_bar(3);
rho_bar(1)=rho(1);
T_bar(1)=T(1);

    %Correction step%

    for m=2:N-1

rho_der_bar(m)=(1/dx)*(-rho_bar(m)*(V_bar(m)-V_bar(m-1))-rho_bar(m)*V_bar(m)...
    *(log(A(m))-log(A(m-1)))-V_bar(m)*(rho_bar(m)-rho_bar(m-1)));

V_der_bar(m)=(1/dx)*(-V_bar(m)*(V_bar(m)-V_bar(m-1))-(1/gama)*...
    ((T_bar(m)-T_bar(m-1))+(T_bar(m)/rho_bar(m))*(rho_bar(m)-rho_bar(m-1))));

T_der_bar(m)=(1/dx)*(-V_bar(m)*(T_bar(m)-T_bar(m-1))-(gama-1)*T_bar(m)*((V_bar(m)-V_bar(m-1))...
    +(V_bar(m)*(log(A(m))-log(A(m-1))))));

rho(m)=rho(m)+dt*0.5*(rho_der(m)+rho_der_bar(m));
V(m)=V(m)+dt*0.5*(V_der(m)+V_der_bar(m));
T(m)=T(m)+dt*0.5*(T_der(m)+T_der_bar(m));

    end

%Boundary conditions%
V(1)=2*V(2)-V(3);
V(31)=2*V(30)-V(29);
rho(31)=2*rho(30)-rho(29);
T(31)=2*T(30)-T(29);

Q(:,1)=rho;
Q(:,2)=V;
Q(:,3)=T;

j
Q

end

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us