将Matlab数值热传导模型转换为Python
wig*_*ing 3 python matlab numpy numerical-methods python-3.x
我试图将我的Matlab模型转换为瞬态热传导到Python.不幸的是,我在Python中的数值解决方案的输出与Matlab模型的输出不匹配.我正在使用Spyder IDE编写我的代码.
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
蟒蛇...
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 = 99Python实线上运行上面的代码与实体Matlab红线不匹配,但Python和Matlab图的虚线确实一致.这告诉我for在Python代码的最后一个循环中也出现了问题.也许我在Python中解决A*x = B的方式不正确？
生成的循环中的索引值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并没有改变,所以你不应该改变m,以m-1在计算边缘的公式A和
C.也就是说,改变这些行:
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.