(精品word)拟一维喷管流动的数值解法(MATLAB)代码-豆柴文库

您所在位置：网站首页 / (精品word)拟一维喷管流动的数值解法(MATLAB)代码.doc / 文档详情

在线预览结束，喜欢就下载吧，查找使用更方便

5 金币

下载文档

/ 3

下载提示

如果您无法下载资料，请参考说明：

1、部分资料下载需要金币，请确保您的账户上有足够的金币

2、已购买过的文档，再次下载不重复扣费

3、资料包下载后请先用软件解压，在使用对应软件打开

文本预览

(精品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;
end
end

图形处理代码
closeall；
a=［151101151201701］；
plot（x,M（a（1),:).*A（：)'。*V