(精品word)拟一维喷管流动的数值解法(MATLAB)代码(精品word)拟一维喷管流动的数值解法(MATLAB)代码(精品word)拟一维喷管流动的数值解法(MATLAB)代码拟一维喷管流动的数值解法（MATLAB）代码数值计算代码%拟一维喷管流动的数值解%亚声速-超声速，非守恒形式functionmain()clear；clc;r=1。4;％绝热指数N=1001；%时间步长i=31；％网格数目L=3；％喷管长度C=0.5;%柯朗数dx=L/(i-1）;%空间步长dt（N）=0；%时间步长x=linspace(0,L，i)；％网格点横坐标A=1+2。2*（x—1。5).^2；%喷管面积%赋值M(N,i)=0；T（N,i）=0;V（N，i)=0；%初始条件M(1，：）=1—0.3146＊x；T（1,:)=1—0.2314＊x;V（1，：）=（0.1+1。09*x）。*(1—0。2314＊x)。^0。5；％按时间步长推进fork=1:N—1％预估偏导数M_t(1:i—1)=—V(k，1：i—1）.＊(M（k，2：i）—M(k，1：i—1）)/dx-M(k，1：i-1)。＊（V（k，2:i）-V(k,1：i—1））/dx-M(k，1:i-1）。*V(k，1：i—1).*log（A（2：i）./A（1：i-1)）/dx;V_t(1:i-1）=—V（k，1：i—1）。＊(V（k,2：i）-V（k，1:i—1))/dx-1/r.＊((T(k,2:i)—T（k,1：i-1））/dx+T（k,1:i—1）。/M（k，1:i-1).*（M（k,2：i）-M(k,1：i-1))/dx)；T_t(1：i-1）=—V(k，1：i—1）.＊(T(k，2:i）-T（k,1：i—1））/dx-(r-1).*T（k,1:i—1）。＊（（V(k，2：i）-V（k,1:i-1)）/dx+V(k，1:i-1).*log(A（2:i)。/A(1:i—1））/dx）;％求取内部网格点处最小时间步长t=C*dx。/（V(k，2：i-1)+sqrt（T(k,2：i—1）)）;dt(k）=min(t);%预估值M1（1：i—1）=M(k，1：i—1）+M_t(1:i-1)＊dt（k);V1（1:i—1)=V(k，1：i—1）+V_t(1：i—1）*dt（k）；T1（1:i—1）=T（k,1：i-1）+T_t(1：i—1）＊dt（k)；%校正偏导数M_t_1（2：i-1）=-V1（2：i—1).＊（M1(2：i-1）-M1(1:i-2)）./dx—M1(2:i-1）.*(V1(2：i—1)-V1(1:i-2））./dx—M1（2：i-1)。*V1（2:i-1).＊log（A（2:i—1)。/A(1：i—2）)。/dx；V_t_1(2：i—1)=-V1（2：i—1）.*(V1（2:i—1)—V1（1:i-2）)./dx-1/r.＊((T1（2:i—1)—T1（1：i—2))。/dx+T1(2：i—1）。/M1（2：i—1）。＊（M1（2:i—1）—M1(1：i—2））。/dx）；T_t_1(2:i-1)=-V1(2:i-1)。*(T1（2：i—1)—T1(1:i-2）)./dx—（r-1)。＊T1(2：i-1）。*((V1(2：i—1)-V1（1:i-2））./dx+V1(2：i—1）.＊log（A（2:i-1)。/A(1：i-2)）./dx)；％偏导数平均值M_t_av（2:i-1)=0.5＊（M_t（2：i—1)+M_t_1（2：i—1）);V_t_av(2:i—1)=0.5＊（V_t（2:i-1）+V_t_1(2：i—1）);T_t_av(2:i—1)=0.5*(T_t(2：i-1）+T_t_1（2：i—1））；%内部网格点修正值M（k+1,2:i—1)=M(k,2:i—1)+M_t_av(2:i-1）*dt(k)；V(k+1,2：i-1)=V（k，2:i—1）+V_t_av(2：i-1)＊dt(k);T(k+1,2：i-1)=T(k,2：i-1）+T_t_av(2：i-1)＊dt（k);％出口边界值M（k+1,i）=2*M(k+1,i-1)-M（k+1，i-2)；V(k+1,i)=2*V（k+1,i-1）—V（k+1,i—2）；T(k+1,i)=2*T（k+1,i—1)—T（k+1，i—2）;％入口边界值M（k+1，1）=1;V(k+1,1)=2＊V(k+1,2）—V(k+1，3）;T(k+1,1)=1;endend图形处理代码closeall；a=［151101151201701］；plot（x,M（a（1),:).*A（：)'。*V