将Matlab热导模型转换为Python
我正在尝试把我的Matlab模型转换成Python，用于瞬态热传导的计算。不过，我在Python中得到的数值结果和Matlab模型的结果不一致。我使用Spyder这个软件来写我的代码。
到目前为止，我发现Matlab和Python在我的模型中有几个主要的区别：
- 在Matlab中，数组的索引是从1开始的，而在Python中，索引是从0开始的。
- 在Matlab中，解方程A*x = B的方法是用
x = A \ B，而在Python中则是x = np.linalg.solve(A,B)。 - 在Matlab中，把列向量
col = C转换成行向量row的方法是row = C'，而在Python中是row = C.T。
为了检查Matlab和Python模型的结果，我对比了A数组。结果显示，这两个数组并不相同：
Matlab...
A =
1.1411 -0.1411 0 0
-0.0118 1.0470 -0.0353 0
0 -0.0157 1.0470 -0.0313
0 0 -0.0470 1.0593
Python...
A
[[ 1.14106122 -0.14106122 0. 0. ]
[-0. 1.04702041 -0.04702041 0. ]
[ 0. -0.0235102 1.04702041 -0.0235102 ]
[ 0. 0. -0.04702041 1.05681633]]
我觉得我的Python代码中有些地方做得不对。我猜可能和Python的数组索引方式有关，但我不太确定。
所以，如果有人能给我一些建议，告诉我如何在Python中构建我的Matlab模型，我会非常感激。
这是我想在Python中复制的Matlab示例：
% parameters
% -------------------------------------------------------------------------
rho = 700; % density of wood, kg/m^3
d = 0.035e-2; % wood particle diameter, m
cpw = 1500; % biomass specific heat capacity, J/kg*K
kw = 0.105; % biomass thermal conductivity, W/m*K
h = 375; % heat transfer coefficient, W/m^2*K
Ti = 300; % initial particle temp, K
Tinf = 773; % ambient temp, K
% numerical model where b = 1 cylinder and b = 2 sphere
% -------------------------------------------------------------------------
nt = 1000; % number of time steps
tmax = 0.8; % max time, s
dt = tmax/nt; % time step, s
t = 0:dt:tmax; % time vector, s
nr = 3; % number or radius steps
%nr = 100; % number or radius steps
r = d/2; % radius of particle, m
dr = r/nr; % radius step, delta r
m = nr+1; % nodes from center to surface
b = 2 ; % run model as a cylinder (b = 1) or as a sphere (b = 2)
if b == 1
shape = 'Cylinder';
elseif b == 2
shape = 'Sphere';
end
alpha = kw/(rho*cpw); % thermal diffusivity, alfa = kw / rho*cp, m^2/s
Fo = alpha*dt/(dr^2); % Fourier number, Fo = alfa*dt / dr^2, (-)
Bi = h*dr/kw; % Biot numbmer, Bi = h*dr / kw, (-)
% creat array [TT] to store temperature values, row = time step, column = node
TT = zeros(1,m);
i = 1:m;
TT(1,i) = Ti; % first row is initial temperature of the cylinder or sphere
% build coefficient matrix [A] and initial column vector {C}
A = zeros(m); % pre-allocate [A] array
C = zeros(m,1); % pre-allocate {C} vector
A(1,1) = 1 + 2*(1+b)*Fo;
A(1,2) = -2*(1+b)*Fo;
C(1,1) = Ti;
for i = 2:m-1
A(i,i-1) = -Fo*(1 - b/(2*i)); % Tm-1
A(i,i) = 1 + 2*Fo; % Tm
A(i,i+1) = -Fo*(1 + b/(2*i)); % Tm+1
C(i,1) = Ti;
end
A(m,m-1) = -2*Fo;
A(m,m) = 1 + 2*Fo*(1 + Bi + (b/(2*m))*Bi);
C(m) = Ti + 2*Fo*Bi*(1 + b/(2*m))*Tinf;
% display [A] array and [C] column vector in console
A
C
% solve system of equations [A]{T} = {C} for column vector {T}
for i = 2:nt+1
T = A\C;
C = T;
C(m) = T(m) + 2*Fo*Bi*(1 + b/(2*m))*Tinf;
TT(i,:) = T'; % store new temperatures in array [TT]
end
% plot
% -------------------------------------------------------------------------
figure(b)
plot(t,TT(:,1),'--k',t,TT(:,m),'-k')
hold on
plot([0 tmax],[Tinf Tinf],':k')
hold off
axis([0 tmax Ti-20 Tinf+20])
ylabel('Temperature (K)')
xlabel('Time (s)')
nr = num2str(nr); nt = num2str(nt); dt = num2str(dt); h = num2str(h); Tinf = num2str(Tinf);
legend('center','surface',['T\infty = ',Tinf,'K'],'location','southeast')
title([num2str(shape),', nr = ',nr,', nt = ',nt,', \Deltat = ',dt,', h = ',h])
这是我在Python中的尝试：
# use Python 3 print function
from __future__ import print_function
# libraries and packages
import numpy as np
import matplotlib.pyplot as py
# parameters
# -------------------------------------------------------------------------
rho = 700 # density of wood, kg/m^3
d = 0.035e-2 # wood particle diameter, m
cpw = 1500 # biomass specific heat capacity, J/kg*K
kw = 0.105 # biomass thermal conductivity, W/m*K
h = 375 # heat transfer coefficient, W/m^2*K
Ti = 300 # initial particle temp, K
Tinf = 773 # ambient temp, K
# numerical model where b = 1 cylinder and b = 2 sphere
# -------------------------------------------------------------------------
nt = 1000 # number of time steps
tmax = 0.8 # max time, s
dt = tmax/nt # time step, s
t = np.arange(0,tmax+dt,dt)
nr = 3 # number or radius steps
r = d/2 # radius of particle, m
dr = r/nr # radius step, delta r
m = nr+1 # nodes from center m=0 to surface m=steps+1
b = 2 # run model as a cylinder (b = 1) or as a sphere (b = 2)
alpha = kw/(rho*cpw) # thermal diffusivity, alfa = kw / rho*cp, m^2/s
Fo = alpha*dt/(dr**2) # Fourier number, Fo = alfa*dt / dr^2, (-)
Bi = h*dr/kw # Biot numbmer, Bi = h*dr / kw, (-)
# create array [TT] to store temperature values, row = time step, column = node
TT = np.zeros((1,m))
# first row is initial temperature of the cylinder or sphere
for i in range(0,m):
TT[0,i] = Ti
# build coefficient matrix [A] and initial column vector {C}
A = np.zeros((m,m)) # pre-allocate [A] array
C = np.zeros((m,1)) # pre-allocate {C} vector
A[0, 0] = 1 + 2*(1+b)*Fo
A[0, 1] = -2*(1+b)*Fo
C[0, 0] = Ti
for i in range(1, m-1):
A[i, i-1] = -Fo*(1 - b/(2*i)) # Tm-1
A[i, i] = 1 + 2*Fo # Tm
A[i, i+1] = -Fo*(1 + b/(2*i)) # Tm+1
C[i, 0] = Ti
A[m-1, m-2] = -2*Fo
A[m-1, m-1] = 1 + 2*Fo*(1 + Bi + (b/(2*(m-1)))*Bi)
C[m-1, 0] = Ti + 2*Fo*Bi*(1 + b/(2*(m-1)))*Tinf
# print [A] and [C] to console
print('A \n', A)
print('C \n', C)
# solve system of equations [A]{T} = {C} for column vector {T}
for i in range(1, nt+1):
T = np.linalg.solve(A,C)
C = T
C[m-1, 0] = T[m-1, 0] + 2*Fo*Bi*(1 + b/(2*(m-1)))*Tinf
TT = np.vstack((TT, T.T))
# plot results
py.figure(1)
py.plot(t,TT[:, m-1])
py.plot(t,TT[:, 0])
py.grid()
py.show()
关于生成的Python图表（见下图），实线（红色和黑色）和虚线（红色和黑色）应该重合。在nr = 99的情况下运行上述代码时，Python的实线与Matlab的红色实线不匹配，但Python和Matlab的虚线是相符的。这让我觉得Python代码的最后一个for循环也有问题。也许我在Python中解A*x = B的方法不正确？
1 个回答
5
在生成 A 的循环中，索引的值变了1，所以这两行
A[i, i-1] = -Fo*(1 - b/(2*i)) # Tm-1
A[i, i+1] = -Fo*(1 + b/(2*i)) # Tm+1
应该改成
A[i, i-1] = -Fo*(1 - b/(2*(i+1))) # Tm-1
A[i, i+1] = -Fo*(1 + b/(2*(i+1))) # Tm+1
注意公式中 i 变成了 i+1（但 A 的索引没有变）。
另一方面，m 的含义没有改变，所以在计算 A 和 C 边缘的公式中，你不应该把 m 改成 m-1。也就是说，把这些行：
A[m-1, m-1] = 1 + 2*Fo*(1 + Bi + (b/(2*(m-1)))*Bi)
C[m-1, 0] = Ti + 2*Fo*Bi*(1 + b/(2*(m-1)))*Tinf
...
C[m-1, 0] = T[m-1, 0] + 2*Fo*Bi*(1 + b/(2*(m-1)))*Tinf
改成
A[m-1, m-1] = 1 + 2*Fo*(1 + Bi + (b/(2*m))*Bi)
C[m-1, 0] = Ti + 2*Fo*Bi*(1 + b/(2*m))*Tinf
...
C[m-1, 0] = T[m-1, 0] + 2*Fo*Bi*(1 + b/(2*m))*Tinf
另外，像 @eryksun 指出的那样，C = T 应该改为 C = T.copy()。在 Matlab 中，内存是通过“写时复制”来管理的，所以在这个赋值之后对 C 的原地修改不会影响到 T。而在 numpy 中，C = T 会让 C 和 T 指向同一个底层数组对象；如果你在 C 上做修改，也会改变 T。为了让它像 Matlab 那样工作，C 必须是 T 的一个副本。