I am trying to convert my Matlab model for transient heat conduction to Python. Unfortunately the output from my numerical solution in Python is not matching the output from the Matlab model. I am using the Spyder IDE to write my code.
The main differences between Matlab and Python for my model (that I have found so far):
- Matlab indices for arrays start at 1 whereas Python index starts at 0
- the solution for A*x = B in Matlab is
x = A \ B whereas in Python it is x = np.linalg.solve(A,B)
- changing a column vector
col = C to a row vector row in Matlab is row = C' whereas in Python it is row = C.T
As a way to check the Matlab and Python model, I compare the A array. As you can see, the two arrays do not match:
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]]
There is something that I am not doing right in the Python code. My guess is that it has to do with how Python indexes arrays. But I'm not sure.
So any suggestions on how to construct my Matlab model in Python would be greatly appreciated?
Here is the Matlab example that I am trying to replicate in Python:
% 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])
And here is my attempt in 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()
In regards to the Python plot generated (see image below), the solid (red & black) lines and dashed (red & black) lines should be on top of each other. Running the above code at nr = 99 the Python solid line does not match the solid Matlab red line, but the dashed lines of the Python and Matlab plots do agree. This tells me that something is also wrong in the last for loop of the Python code. Maybe the way I'm solving A*x = B in Python is not correct?
https://stackoverflow.com/questions/22461342
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
Russian