本文以求解拟一维喷管流动为例,比较两者在科学计算中的区别。
感受:matlab矩阵实验室在求解矩阵方面具有得天独厚的优势,尤其是在矩阵之间的运算方面、求解方程过程中,能够明显感觉到编程给人带来的快感,调试过程中实时显示运算结果是当今任何一款编程软件都无法比拟的。而近年来非常火的python,同样有其优势(可能在大数据处理,人工智能领域吧),在本次实验中表现差强人意(勉强使人满意)。
本次实验采用二阶精度的MacCormack差分格式求解CFD问题。具体分析流程详见Anderrson。
matlab代码如下(搬运):
%从亚音速到超音速的转变,非守恒形式
%实现从t时刻到t+dt时刻的转换,rou,V,T
%rou rou2 分别表示rou(i),rou(i+1)
%A,gamma,dt,dx
clear,clc
NUM=31;
C=0.5;%计算时间dt
dx=3/(NUM-1);
step=1000;%时间模拟的步数
Time=0;
gamma=1.4;
x=0:dx:3;
A=1+2.2*(x-1.5).^2;
rou_ar=zeros(step+1,NUM);V_ar=zeros(step+1,NUM);T_ar=zeros(step+1,NUM);
Ma_ar=zeros(step+1,NUM);flux_m=zeros(step+1,NUM);P_ar=zeros(step+1,NUM);
% rou=zeros(1,NUM);V=zeros(1,NUM);T=zeros(1,NUM);
drou_t=zeros(1,NUM);dV_t=zeros(1,NUM);dT_t=zeros(1,NUM);
drou_t2=zeros(1,NUM);dV_t2=zeros(1,NUM);dT_t2=zeros(1,NUM);
drou_t3=zeros(1,NUM);dV_t3=zeros(1,NUM);dT_t3=zeros(1,NUM);
rou_e=zeros(1,NUM);V_e=zeros(1,NUM);T_e=zeros(1,NUM);
rou_r=zeros(1,NUM);V_r=zeros(1,NUM);T_r=zeros(1,NUM);
dt_ar=zeros(1,step);
%初始化
rou=1-0.3146*x;
T=1-0.2314*x;
V=(0.1+1.09*x).*T.^0.5;
%
rou_ar(1,:)=rou;
V_ar(1,:)=V;
T_ar(1,:)=T;
Ma_ar(1,:)=V./sqrt(T);
flux_m(1,:)=rou.*V.*A;
P_ar(1,:)=rou.*T;
for k=1:1:step %模拟的时间步数
rou=rou_ar(k,:);
V=V_ar(k,:);
T=T_ar(k,:);
%确定本次模拟的dt时间间隔
dt_list=C*dx./(T(2:NUM-1).^0.5+V(2:NUM-1));
dt=min(dt_list);
%
%预测步
for i=2:1:(NUM-1)
drou_t(i)=-rou(i)*(V(i+1)-V(i))/dx-rou(i)*V(i)*(log(A(i+1))-log(A(i)))/dx-V(i)*(rou(i+1)-rou(i))/dx;
dV_t(i)=-V(i)*(V(i+1)-V(i))/dx-1/gamma*((T(i+1)-T(i))/dx+T(i)/rou(i)*((rou(i+1)-rou(i))/dx));
dT_t(i)=-V(i)*(T(i+1)-T(i))/dx-(gamma-1)*T(i)*((V(i+1)-V(i))/dx+V(i)*(log(A(i+1))-log(A(i)))/dx);
rou_e(i)=rou(i)+drou_t(i)*dt;
V_e(i)=V(i)+dV_t(i)*dt;
T_e(i)=T(i)+dT_t(i)*dt;
end
%BC
V_e(1)=2*V_e(2)-V_e(3);
rou_e(1)=1;
T_e(1)=1;
rou_e(NUM)=2*rou_e(NUM-1)-rou_e(NUM-2);
V_e(NUM)=2*V_e(NUM-1)-V_e(NUM-2);
T_e(NUM)=2*T_e(NUM-1)-T_e(NUM-2);
%%%%%%%%%%%%
%预测步结束
%修正步
for i=2:1:(NUM-1)
drou_t2(i)=-rou_e(i)*(V_e(i)-V_e(i-1))/dx-rou_e(i)*V_e(i)*(log(A(i))-log(A(i-1)))/dx-V_e(i)*(rou_e(i)-rou_e(i-1))/dx;
d