传热 PDE

编程问答 2022-05-15

问题描述

我一直在尝试使用显式和隐式欧拉方法来求解一个简单的传热偏微分方程：

∂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 (将#修改为@）

mathmatlabmatlabnumerical-methodspde