initial_data.m

这几个算例都是MacCormack格式的,喷管计算用时间步进法,外突钝角绕流用空间步进法.都写成模块形式,便于查错.

function initial_data(t,N,tU2)
%initial the data of nozzle flow
global rou T V U1 U2 U3 F1 F2 F3 J2 x dx A plnA_px pU1_pt pU2_pt pU3_pt gama;
gama=1.4;
dx=3/(N-1);
x=0:dx:3;

A=NF_Shape(N);
plnA_px=zeros(1,N);
for i=2:(N-1)
    plnA_px(i)=(log(A(i+1))-log(A(i)))/dx;
end

rou=zeros(t+1,N);
T=zeros(t+1,N);
V=zeros(t+1,N);
init_condition(N,tU2);
U1=zeros(t+1,N);
U2=zeros(t+1,N);
U3=zeros(t+1,N);
F1=zeros(t+1,N);
F2=zeros(t+1,N);
F3=zeros(t+1,N);
J2=zeros(t+1,N);
U1(1,:)=rou(1,:).*A;
U2(1,:)=U1(1,:).*V(1,:);
U3(1,:)=rou(1,:).*(T(1,:)./(gama-1)+(gama./2).*V(1,:).^2).*A;

pU1_pt=zeros(t,N);
pU2_pt=zeros(t,N);
pU3_pt=zeros(t,N);


function y=NF_Shape(N)
%Nozzle Shape
global x;
  y=1+2.2*(x-1.5).^2;
  
function init_condition(N,tU2)
global rou T V x A;
  N1=floor(N/6);
  N2=floor(N/2);
  
  rou(1,1:N1)=1;
  rou(1,N1+1:N2)=1-0.366.*(x(N1+1:N2)-0.5);
  rou(1,N2+1:N)=0.634-0.3879.*(x(N2+1:N)-1.5);
  

  T(1,1:N1)=1;
  T(1,N1+1:N2)=1-0.167.*(x(N1+1:N2)-0.5);
  T(1,N2+1:N)=0.833-0.3507.*(x(N2+1:N)-1.5);
  
  V(1,:)=tU2./(rou(1,:).*A);