如何解决传热 PDE
我一直在尝试使用显式和隐式欧拉方法来求解一个简单的传热偏微分方程:
∂T/∂t = alpha * (∂^2T/∂x^2)
T = 温度,x = 轴向尺寸。我使用的初始条件 (I.C.) 是 x = 0,T = 100 °C。以及计算网格末端的边界条件(BC):对于 T(x = L) = T(x = L-1),其中 L是计算网格的长度->最后一个和倒数第二个温度值相同。
显式方法(工作正常):T 的每个值都由 T1(i) + heat_coefficient*((T1(i+1)-2*T1(i)+T1(i-1))/dx^2)*dt
计算,除了第一个和最后一个值由我知道了和卑诗省,分别。
我的问题:
隐式方法:首先,我创建了一个三对角矩阵MAT
,它定义了下一个时间线(n+1)中的值之间的关系.在隐式方法的情况下,我不能写 I.C.和 BC完全正确,因此我将它们保存在“等式的右侧”,即 pS(i) = -(T2(i)*dx^2)/(heat_coefficient*dt);
这样,我在一个时间线 (n) 中表达所有值,然后在下一个时间线 (n+1) 中继续。
如果我把 I.C.或卑诗省进入“等式的右侧”,结果对位置 (dx) 和时间 (dt) 步长的任何变化都变得非常敏感。因此,温度函数的行为是错误的:温度曲线应该在一定时间后收敛到相同的值(如果时间变为无穷大),但是,我的曲线收敛到不同的值(并且根据 随机变化dx 和 dt).
如何实现在 100 °C 的温度下启动,以便在“无限”长时间后,所有温度曲线都收敛到 100 °C?
I.C. 的值是否应该?和 BC被放入一个三对角矩阵 - 例如,如果我有维度为 [N,N] 的矩阵,那么我指定 I.C.对于点 [1,1] 和 B.C.对于点 [N,N]? (很遗憾,我尝试的时候并没有正常工作)
另外,我对三对角矩阵的实现可能不是很好,但它有效。
这是我的代码:
% Heat transfer numerical solution test
% Implicit and explicit formula
function Temperature
clear all
clc
heat_coefficient = 0.03;
dt = 0.0001; % time step (s)
dx = 0.01; % positional step (m)
time = 1; % (s)
length = 0.5 % (m)
t_count = fix(time/dt) % number of time nodes
x_count = fix(length/dx) % number of positional nodes
T1 = zeros(x_count,1); % Explicit method
T2 = zeros(x_count,1); % Implicit method
pS = zeros(x_count,1); % Implicit method (right side of the equation,constant value)
MAT = zeros(x_count,x_count); % Tri-diagonal matrix - for implicit method
curve_count = 10; % Number of curves that are plotted
time_count = 0; % Time cycle count
% Tri-diagonal matrix diagonals specification
A = 1;
B = -(dx^2/(heat_coefficient*dt) + 2); % Main diagonal
C = A;
for j=1:t_count % Cycle over time nodes
for i=1:x_count % Cycle over positional nodes
if i == 1
T1(i) = 100; % Initial condition for explicit method
T2(i) = 100; % Initial condition for implicit method
MAT(i,i) = B;
MAT(i,i+1) = C;
elseif i == x_count
T1(i) = T1(i-1); % B.C. for explicit method (last equals the next-to-last)
T2(i) = T2(i-1); % B.C. for implicit method (last equals the next-to-last)
MAT(i,i-1)= A;
MAT(i,i) = B;
else
% Euler's explicit discretization
T1(i) = T1(i) + heat_coefficient*((T1(i+1)-2*T1(i)+T1(i-1))/dx^2)*dt;
% Implicit method - creation of tri-diagonal matrix
MAT(i,i-1)=A;
MAT(i,i) =B;
MAT(i,i+1)=C;
end
% Right side of the equation
pS(i) = -(T2(i)*dx^2)/(heat_coefficient*dt);
end
% PDEs,function
T2 = EQsystem(MAT,pS);
% Condition to plot only a certain number of temperature curves
if fix(j/(t_count/curve_count)) == j/(t_count/curve_count)
if time_count == 0
y_out1 = T1;
y_out2 = T2;
time_count = time_count + 1;
else
y_out1 = [y_out1 T1];
y_out2 = [y_out2 T2];
time_count = time_count + 1;
end
end
end
x = linspace(0,length,x_count); % Nodes along the axial axis
% Graph
subplot(2,1,2);
plot(x,y_out1); % Explicit
title('Explicit')
subplot(2,1);
plot(x,y_out2); % Implicit
title('Implicit')
end
function x=EQsystem(A,b)
% Algebraic EQs solver
% A ... matrix
% b ... right side vector
% x ... solution of the algebraic EQs
if det(A)==0
error('!')
end
x = A\b; % solution
% x = inv(A)*b; % Another possibility to obtain a solution
end
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。