matlab 与 python 在科学计算中的区别比较

CFD 专栏收录该内容
1 篇文章 0 订阅

本文以求解拟一维喷管流动为例,比较两者在科学计算中的区别。
感受: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;
dV_t2(i)=-V_e(i)*(V_e(i)-V_e(i-1))/dx-1/gamma*((T_e(i)-T_e(i-1))/dx+T_e(i)/rou_e(i)*((rou_e(i)-rou_e(i-1))/dx));
dT_t2(i)=-V_e(i)*(T_e(i)-T_e(i-1))/dx-(gamma-1)*T_e(i)*((V_e(i)-V_e(i-1))/dx+V_e(i)*(log(A(i))-log(A(i-1)))/dx);

%对导数求平均
drou_t3(i)=0.5*(drou_t(i)+drou_t2(i));
dV_t3(i)=0.5*(dV_t(i)+dV_t2(i));
dT_t3(i)=0.5*(dT_t(i)+dT_t2(i));

rou_r(i)=rou(i)+drou_t3(i)*dt;
V_r(i)=V(i)+dV_t3(i)*dt;
T_r(i)=T(i)+dT_t3(i)*dt;
end

%BC
V_r(1)=2*V_r(2)-V_r(3);
rou_r(1)=1;
T_r(1)=1;

rou_r(NUM)=2*rou_r(NUM-1)-rou_r(NUM-2);
V_r(NUM)=2*V_r(NUM-1)-V_r(NUM-2);
T_r(NUM)=2*T_r(NUM-1)-T_r(NUM-2);

rou_ar(k+1,:)=rou_r;
V_ar(k+1,:)=V_r;
T_ar(k+1,:)=T_r;
Ma_ar(k+1,:)=V_r./sqrt(T_r);
flux_m(k+1,:)=rou_r.*V_r.*A;
P_ar(k+1,:)=rou_r.*T_r;

dt_ar(k)=dt;
Time=Time+dt;
end  %模拟的时间步数

%理论解的计算
x_theory=0:0.01:3;
A_theory=1+2.2*(x_theory-1.5).^2;
num_x_th=size(x_theory,2);
Ma_theory=zeros(1,num_x_th);
%求解理论上的Ma解
for i=1:1:num_x_th
    F=@(ma)(1/ma^2)*((2/(gamma+1))*(1+(gamma-1)/2*ma^2))^((gamma+1)/(gamma-1))-A_theory(i)^2;
    if i<=(num_x_th-1)/2+1
    Ma_theory(i)=fzero(F,[1e-8,1]);
    else
    Ma_theory(i)=fzero(F,[1,10]);
    end
end

P_theory=(1+((gamma-1)/2)*(Ma_theory(end,:).^2)).^(-gamma/(gamma-1));
rou_theory=(1+((gamma-1)/2)*(Ma_theory(end,:).^2)).^(-1/(gamma-1));
T_theory=(1+((gamma-1)/2)*(Ma_theory(end,:).^2)).^(-1);

%模拟解
P_simu=P_ar(end,:);
rou_simu=rou_ar(end,:);
T_simu=T_ar(end,:);
Ma_simu=Ma_ar(end,:);


%后处理
%针对i节点,显示其在整个模拟过程中的变化。
xn=0:1:step;
i=(NUM-1)/2+1;%设置所要观察的节点
rou_i=rou_ar(:,i);
V_i=V_ar(:,i);
T_i=T_ar(:,i);
Ma_i=Ma_ar(:,i);

subplot(2,2,1)
plot(xn,rou_i);
title('密度');

subplot(2,2,2)
plot(xn,V_i);
title('V');

subplot(2,2,3)
plot(xn,T_i);
title('T');

subplot(2,2,4)
plot(xn,Ma_i);
title('Ma');

%质量流量的变化
figure(2)
flux_list=[0,50,150,200,700]+1;%注:所研究的最大700,不能超过step
plot(x,flux_m(flux_list(1),:),x,flux_m(flux_list(2),:),x,flux_m(flux_list(3),:),x,flux_m(flux_list(4),:),x,flux_m(flux_list(5),:));
legend('0dt','50dt','150dt','200dt','700dt');
title('质量流量');
%稳态的理论解和模拟解的对比
%理论解


%画图对比
figure(3)
subplot(2,2,1)
plot(x_theory,P_theory,'-',x,P_simu,'+');
legend('理论解','模拟解');
title('Pressure comparision');

subplot(2,2,2)
plot(x_theory,rou_theory,'-',x,rou_simu,'+');
legend('理论解','模拟解');
title('Density comparision');

subplot(2,2,3)
plot(x_theory,T_theory,'-',x,T_simu,'+');
legend('理论解','模拟解');
title('Temperature comparision');

subplot(2,2,4)
plot(x_theory,Ma_theory,'-',x,Ma_simu,'+');
legend('理论解','模拟解');
title('Mach comparision');

与之对应的python3 源码如下(原创):

import numpy
from scipy.optimize import fsolve
import matplotlib.pyplot as plt
for i in range(10):
    print(i)

NUM = 31
C = 0.5
dx = 3/(NUM-1)
step = 1000
Time = 0
gamma = 1.4
x = numpy.arange(0, 3.1, dx)
A = 1 + 2.2 * pow((x-1.5),2)

rou_ar = numpy.zeros([step+1, NUM])
V_ar = numpy.zeros([step+1, NUM])
T_ar = numpy.zeros([step+1, NUM])
Ma_ar = numpy.zeros([step+1,NUM])
flux_m =numpy.zeros([step+1,NUM])
P_ar = numpy.zeros([step+1,NUM])

drou_t = numpy.zeros([1,NUM])
dV_t =numpy.zeros([1,NUM])
dT_t = numpy.zeros([1,NUM])
drou_t2 = numpy.zeros([1,NUM])
dV_t2 = numpy.zeros([1,NUM])
dT_t2 = numpy.zeros([1,NUM])
drou_t3 = numpy.zeros([1,NUM])
dV_t3 = numpy.zeros([1,NUM])
dT_t3 = numpy.zeros([1,NUM])
rou_e =numpy.zeros([1,NUM])
V_e = numpy.zeros([1,NUM])
T_e = numpy.zeros([1,NUM])
rou_r = numpy.zeros([1,NUM])
V_r = numpy.zeros([1,NUM])
T_r = numpy.zeros([1,NUM])
dt_ar=numpy.zeros([1,step])
print(drou_t[0,1])

rou = 1 - 0.3146 * x
T = 1 - 0.2314 * x
V = (0.1 + 1.09 * x) * numpy.sqrt(T)
print(V[0])
for i in range(NUM):
    rou_ar[0,i] = rou[i]
    V_ar[0,i]=V[i]
    T_ar[0,i]=T[i]
    Ma_ar[0,i]=V[i]/numpy.sqrt(T[i])
    flux_m[0,i]=rou[i]*V[i]*A[i]
    P_ar[0,i]=rou[i]*T[i]

for k in range(step):
    for j in range(NUM):
        rou[j] = rou_ar[k,j]
        V[j] = V_ar[k,j]
        T[j] = T_ar[k,j]
    print(T)
    print(T[0])
    print(T[1:5])
    t = numpy.sqrt(T[1:])
    v = numpy.sqrt(V[1:])
    dt_list = C * dx/(numpy.sqrt(T[1:])+ V[1:])
    # print(dt_list)
    dt = min(dt_list)

    for i in range(1,NUM-1,1):
        drou_t[0,i] = -rou[i] * (V[i+1] - V[i]) / dx - rou[i] * V[i] * (numpy.log(A[ i + 1]) - numpy.log(A[i])) / dx - V[i] * (rou[i + 1] - rou[i]) / dx
        dV_t[0,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[0,i] = -V[i] * (T[i + 1] - T[i]) / dx - (gamma - 1) * T[i] * ((V[i + 1] - V[i]) / dx + V[i] * (numpy.log(A[i + 1]) - numpy.log(A[i])) / dx)
        rou_e[0,i] = rou[i] + drou_t[0,i] * dt
        V_e[0,i] = V[i] + dV_t[0,i] * dt
        T_e[0,i] = T[i] + dT_t[0,i] * dt

    V_e[0,0] = 2 * V_e[0,1] - V_e[0,2]
    rou_e[0,0] = 1
    T_e[0,0] = 1
    rou_e[0,NUM-1] = 2 * rou_e[0,NUM - 2] - rou_e[0,NUM - 3]
    V_e[0,NUM-1] = 2 * V_e[0,NUM - 2] - V_e[0,NUM - 3]
    T_e[0,NUM-1] = 2 * T_e[0,NUM - 2] - T_e[0,NUM - 3]

    for i in range(1, NUM-1, 1):
        drou_t2[0,i] = -rou_e[0,i] * (V_e[0,i] - V_e[0,i - 1]) / dx - rou_e[0,i] * V_e[0,i] * (numpy.log(A[i]) - numpy.log(A[i - 1])) / dx - V_e[0,i] * (rou_e[0,i] - rou_e[0,i - 1]) / dx
        dV_t2[0,i] = -V_e[0,i] * (V_e[0,i] - V_e[0,i - 1]) / dx - 1 / gamma * ((T_e[0,i] - T_e[0,i - 1]) / dx + T_e[0,i] / rou_e[0,i] * ((rou_e[0,i] - rou_e[0,i - 1]) / dx))
        dT_t2[0,i] = -V_e[0,i] * (T_e[0,i] - T_e[0,i - 1]) / dx - (gamma - 1) * T_e[0,i] * ((V_e[0,i] - V_e[0,i - 1]) / dx + V_e[0,i] * (numpy.log(A[i]) - numpy.log(A[i - 1])) / dx)

     
        drou_t3[0,i] = 0.5 * (drou_t[0,i] + drou_t2[0,i])
        dV_t3[0,i] = 0.5 * (dV_t[0,i] + dV_t2[0,i])
        dT_t3[0,i] = 0.5 * (dT_t[0,i] + dT_t2[0,i])

        rou_r[0,i] = rou[i] + drou_t3[0,i] * dt
        V_r[0,i] = V[i] + dV_t3[0,i] * dt
        T_r[0,i] = T[i] + dT_t3[0,i] * dt

    V_r[0,0] = 2 * V_r[0,1] - V_r[0,2]
    rou_r[0,0] = 1
    T_r[0,0] = 1
    rou_r[0,NUM - 1] = 2 * rou_r[0,NUM - 2] - rou_r[0,NUM - 3]
    V_r[0,NUM - 1] = 2 * V_r[0,NUM - 2] - V_r[0,NUM - 3]
    T_r[0,NUM - 1] = 2 * T_r[0,NUM - 2] - T_r[0,NUM - 3]

    for m in range(NUM):
       rou_ar[k+1, m] = rou_r[0,m]
       V_ar[k + 1,m]=V_r[0,m]
       T_ar[k + 1,m]=T_r[0,m]
       Ma_ar[k + 1,m]=V_r[0,m]/ numpy.sqrt(T_r[0,m])
       flux_m[k + 1,m]=rou_r[0,m]* V_r[0,m]* A[m]
       P_ar[k + 1,m]=rou_r[0,m]* T_r[0,m]
    dt_ar[0,k] = dt
    Time = Time + dt


x_theory = numpy.arange(0,3.01,0.01)
A_theory = 1 + 2.2 * pow((x_theory - 1.5),2)
# num_x_th = x_theory.shape()




x_theory = numpy.arange(0,3.01,0.01)
A_theory = 1 + 2.2 * pow((x_theory - 1.5),2)
num_x_th = numpy.shape(x_theory)
print(num_x_th[0])
a = num_x_th[0]
Ma_theory = numpy.zeros([1,num_x_th[0]])
gamma = 1.4
print(A_theory)
# fsolve(f,[1,2])
ma = 0.0978
print("ma",Ma_theory[0,1])
i=0
f=(1 / ma ** 2) * pow((2 / (gamma + 1)) * (1 + (gamma - 1) / 2 * ma ** 2), ((gamma + 1) / (gamma - 1)))- A_theory[i]**2
print("f=",f)
print(A_theory[i]**2)
print(pow(A_theory[i],2))
def f(ma,i):
    return (1 / ma ** 2) * pow(((2 / (gamma + 1)) * (1 + (gamma - 1) / 2 * ma ** 2)), ((gamma + 1) / (gamma - 1)))- A_theory[i]**2
#
# print(fsolve(f,[-2,2]))
#
#
bb = numpy.array([1e-8,1])
for i in range(301):
    if i <=(num_x_th[0] - 1)/2 + 1:
        Ma_theory[0,i]=fsolve(f,0.01,i)
        aa=Ma_theory[0,i]
        print(i,aa)
    else:
        Ma_theory[0,i] = fsolve(f, 1.5, i)
        print(i, aa)

P_theory=pow((1+((gamma-1)/2)*(Ma_theory[0,:]**2)),(-gamma/(gamma-1)))
rou_theory=pow((1+((gamma-1)/2)*(Ma_theory[0:]**2)),(-1/(gamma-1)))
T_theory=(1+((gamma-1)/2)/(Ma_theory[0,:]**2))
print("P", P_theory)


P_simu=P_ar[-1,:]
rou_simu=rou_ar[-1,:]
T_simu=T_ar[-1,:]
Ma_simu=Ma_ar[-1,:]
print(P_simu-P_ar[300,:])

xn = numpy.arange(0,step+1,1)
i = (NUM - 1) / 2 + 1
ii= int(i)
print(rou_ar[:,ii])

rou_i = rou_ar[:,ii]
print(rou_ar[:,ii])
V_i = V_ar[:, ii]
T_i = T_ar[:, ii]
Ma_i = Ma_ar[:,ii]

plt.subplot(2,2,1)
plt.plot(xn.T, rou_i)
plt.title("密度")

plt.subplot(2,2,2)
plt.plot(xn.T, V_i)
plt.title("V")

plt.subplot(2,2,3)
plt.plot(xn.T, T_i)
plt.title("T")

plt.subplot(2,2,4)
plt.plot(xn.T, Ma_i)
plt.title("Ma")

plt.show()

ps:python源码中省略了后处理的过程。

python的缺点在于调试、运行速度比较慢,对于我这样的小白,不知道如何利用高级python功能的话极易造成程序冗长,程序运行相对缓慢。

python最致命的缺点在于有些错误位置找不出来,因而修改不便。

python要注意:
x=range(10),那么x是一维数组,按索引查看元素过程中,一定是x[0], x[1], x[2];
y=numpy.zeros(1,2), 生成的y是1行2列的二维数组,引用一定是 y[0,0], y[0,1]
x=30, y=10, a=x/y 在python中可能得出的a 是浮点型,用于数组的索引时程序报错但不会和你说具体原因的,千万要转成int

matlab中的fsolve函数可以求指定区间上函数的零点,fsolve(myfun,[0,5])–求myfun函数在区间[0,5]范围内的零点。
而python中的 from scipy.optimise import fsolve, 同样计算区间[0,5]范围内的零点则会出现问题。因而python中的fsolve一般
fsolve(myfun,start_point),计算从一个点开始最接近改点的函数零点。如果在区间内存在多个零点,那么python中的fsolve如何使用,还请大神们指教。

  • 2
    点赞
  • 2
    评论
  • 2
    收藏
  • 一键三连
    一键三连
  • 扫一扫,分享海报

©️2021 CSDN 皮肤主题: 大白 设计师:CSDN官方博客 返回首页
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值