2024年3月29日发(作者：步念珍)
拟一维喷管流动的数值解法(MATLAB)代码
拟一维喷管流动的数值解法(MATLAB)代码
编辑整理：
尊敬的读者朋友们：
这里是精品文档编辑中心，本文档内容是由我和我的同事精心编辑整理后发布的，发布之前我们对
文中内容进行仔细校对，但是难免会有疏漏的地方，但是任然希望（拟一维喷管流动的数值解法
(MATLAB)代码）的内容能够给您的工作和学习带来便利。同时也真诚的希望收到您的建议和反馈，
这将是我们进步的源泉，前进的动力。
本文可编辑可修改，如果觉得对您有帮助请收藏以便随时查阅，最后祝您生活愉快 业绩进步，以
下为拟一维喷管流动的数值解法(MATLAB)代码的全部内容。
拟一维喷管流动的数值解法(MATLAB)代码
拟一维喷管流动的数值解法（MATLAB）代码
数值计算代码
%拟一维喷管流动的数值解
%亚声速-超声速，非守恒形式
function main（)
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；
拟一维喷管流动的数值解法(MATLAB)代码
T（1，：）=1-0。2314*x；
V(1，：)=(0.1+1。09＊x).*（1—0.2314＊x).^0。5；
％按时间步长推进
for k=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）；
拟一维喷管流动的数值解法(MATLAB)代码
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；
end
end
图形处理代码
close all;
拟一维喷管流动的数值解法(MATLAB)代码
a=［1 51 101 151 201 701］;
plot(x，M（a（1）,:).*A（：)'.*V(a（1）,:），’r--'）
hold on
plot(x，M（a（2)，：).＊A（：）’.*V(a（2)，：）,’m。-’)
plot(x,M(a(3），：）。*A（：)'。＊V(a(3）,:),’g*-’)
plot（x，M（a(4),:)。*A（:)'。*V（a（4），:),'co-')
plot（x,M(a（5）,:).＊A(：)’。*V（a(5),:),'bh-’)
plot(x，M(a（6),：)。*A（:）'。*V(a（6)，：),'k〈-’)
axis equal
legend('0dt’，'50dt',’100dt’，’150dt',’200dt',’700dt'）
figure
［ax,h1,h2］=plotyy（x，M(a(6）,：),x,sqrt（(M(a(6），:）.^(1—r)-1）＊5）,'plot'）;
set（h1,'linestyle','—’，’marker',’o’，'color'，'b’);
set（h2,'linestyle’，’—’,’marker',’*'，’color’，’r’）；
set(ax，’xtick’，0：0。3:3)
set(ax（1),'ytick'，0:0.1：1）
set（ax（2)，’ytick'，0：0。4:4)
set(h1，’linewidth'，1。5)
set（h2,’linewidth',1.5）
set（get(ax（1），’xlabel’），’string'，'无量纲轴向距离（x)’）
set(get（ax(1）,’ylabel’),’string’,’无量纲密度(M）’）
set（get（ax(2），’ylabel')，'string','马赫数（Ma）')
title（’无量纲密度和马赫数定常分布数值解’）
legend('无量纲密度数值解'，’马赫数数值解’）
2024年3月29日发(作者：步念珍)
拟一维喷管流动的数值解法(MATLAB)代码
拟一维喷管流动的数值解法(MATLAB)代码
编辑整理：
尊敬的读者朋友们：
这里是精品文档编辑中心，本文档内容是由我和我的同事精心编辑整理后发布的，发布之前我们对
文中内容进行仔细校对，但是难免会有疏漏的地方，但是任然希望（拟一维喷管流动的数值解法
(MATLAB)代码）的内容能够给您的工作和学习带来便利。同时也真诚的希望收到您的建议和反馈，
这将是我们进步的源泉，前进的动力。
本文可编辑可修改，如果觉得对您有帮助请收藏以便随时查阅，最后祝您生活愉快 业绩进步，以
下为拟一维喷管流动的数值解法(MATLAB)代码的全部内容。
拟一维喷管流动的数值解法(MATLAB)代码
拟一维喷管流动的数值解法（MATLAB）代码
数值计算代码
%拟一维喷管流动的数值解
%亚声速-超声速，非守恒形式
function main（)
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；
拟一维喷管流动的数值解法(MATLAB)代码
T（1，：）=1-0。2314*x；
V(1，：)=(0.1+1。09＊x).*（1—0.2314＊x).^0。5；
％按时间步长推进
for k=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）；
拟一维喷管流动的数值解法(MATLAB)代码
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；
end
end
图形处理代码
close all;
拟一维喷管流动的数值解法(MATLAB)代码
a=［1 51 101 151 201 701］;
plot(x，M（a（1）,:).*A（：)'.*V(a（1）,:），’r--'）
hold on
plot(x，M（a（2)，：).＊A（：）’.*V(a（2)，：）,’m。-’)
plot(x,M(a(3），：）。*A（：)'。＊V(a(3）,:),’g*-’)
plot（x，M（a(4),:)。*A（:)'。*V（a（4），:),'co-')
plot（x,M(a（5）,:).＊A(：)’。*V（a(5),:),'bh-’)
plot(x，M(a（6),：)。*A（:）'。*V(a（6)，：),'k〈-’)
axis equal
legend('0dt’，'50dt',’100dt’，’150dt',’200dt',’700dt'）
figure
［ax,h1,h2］=plotyy（x，M(a(6）,：),x,sqrt（(M(a(6），:）.^(1—r)-1）＊5）,'plot'）;
set（h1,'linestyle','—’，’marker',’o’，'color'，'b’);
set（h2,'linestyle’，’—’,’marker',’*'，’color’，’r’）；
set(ax，’xtick’，0：0。3:3)
set(ax（1),'ytick'，0:0.1：1）
set（ax（2)，’ytick'，0：0。4:4)
set(h1，’linewidth'，1。5)
set（h2,’linewidth',1.5）
set（get(ax（1），’xlabel’），’string'，'无量纲轴向距离（x)’）
set(get（ax(1）,’ylabel’),’string’,’无量纲密度(M）’）
set（get（ax(2），’ylabel')，'string','马赫数（Ma）')
title（’无量纲密度和马赫数定常分布数值解’）
legend('无量纲密度数值解'，’马赫数数值解’）