Matlab:一维热扩散模型的时间步长稳定性

Question

我在Matlab中有一个一维的热扩散代码，我在10年的时间尺度上使用了它，现在我正在尝试使用相同的代码在数百万年的尺度上工作。显然，如果我保持我的时间步长不变，这将需要很长时间来计算，但是如果我增加我的时间步长，我会遇到数值稳定性问题。

我的问题是：

我应该如何处理这个问题？什么会影响最大稳定时间步长？我该如何计算呢？

非常感谢，

亚历克斯

close all
clear all

dx = 4;    % discretization step in m
dt = 0.0000001; % timestep in Myrs
h=1000;        % height of box in m
nx=h/dx+1;
model_lenth=1; %length of model in Myrs
nt=ceil(model_lenth/dt)+1;     % number of tsteps to reach end of model
kappa = 1e-6; % thermal diffusivity
x=0:dx:0+h;     % finite difference mesh
T=38+0.05.*x;  % initial T=Tm everywhere ...
time=zeros(1,nt);
t=0;
Tnew = zeros(1,nx);

%Lower sill
sill_1_thickness=18;
Sill_1_top_position=590;
Sill_1_top=ceil(Sill_1_top_position/dx);
Sill_1_bottom=ceil((Sill_1_top_position+sill_1_thickness)/dx);

%Upper sill
sill_2_thickness=10;
Sill_2_top_position=260;
Sill_2_top=ceil(Sill_2_top_position/dx);
Sill_2_bottom=ceil((Sill_2_top_position+sill_2_thickness)/dx);

%Temperature of dolerite intrusions
Tm=1300;

T(Sill_1_top:Sill_1_bottom)=Tm; %Apply temperature to intrusion 1

% unit conversion to SI:
secinmyr=24*3600*365*1000000;   % dt in sec
dt=dt*secinmyr;

%Plot initial conditions
figure(1), clf
f1 = figure(1); %Make full screen
set(f1,'Units', 'Normalized', 'OuterPosition', [0 0 1 1]); 
plot (T,x,'LineWidth',2)
xlabel('T [^oC]')
ylabel('x[m]')
axis([0 1310 0 1000])
title(' Initial Conditions')
set(gca,'YDir','reverse');

%Main calculation
for it=1:nt

  %Apply temperature to upper intrusion
  if it==10;  
      T(Sill_2_top:Sill_2_bottom)=Tm; 
  end

  for i = 2:nx-1
      Tnew(i) = T(i) + kappa*dt*(T(i+1) - 2*T(i) + T(i-1))/dx/dx;
  end

  Tnew(1) = T(1);
  Tnew(nx) = T(nx);

  time(it) = t;

  T = Tnew; %Set old Temp to = new temp for next loop
  tmyears=(t/secinmyr);

  %Plot a figure which updates in the loop of temperature against depth
  figure(2), clf
  plot (T,x,'LineWidth',2)
  xlabel('T [^oC]')
  ylabel('x[m]')
  title([' Temperature against Depth after ',num2str(tmyears),' Myrs'])
  axis([0 1300 0 1000])
  set(gca,'YDir','reverse');%Reverse y axis

  %Make full screen     
  f2 = figure(2); 
  set(f2,'Units', 'Normalized', 'OuterPosition', [0 0 1 1]); 

  drawnow

  t=t+dt;
end

Answer 1

像FTCS这样的显式格式的稳定性条件由$r =K dt/dx^2 < 1/2$或$dt < dx^2/(2K)$控制，其中K是扩散系数。这是为了使四阶导数前导截断误差项的符号为负所必需的。

如果你不想受到时间步长的限制，我建议使用隐式方案(尽管比显式方案计算成本更高)。这可以通过对扩散项使用后向欧拉而不是前向欧拉来简单地实现。另一种选择是Crank-Nicholson，它也是隐式的。

Answer 2

@Isopycnal振荡是完全正确的，因为最大稳定步长在显式方案中是有限的。仅供参考，这通常被称为离散傅立叶数或仅仅是傅立叶数，并且可以查找不同的边界条件。此外，下面的内容可以帮助您推导隐式格式或克朗克-尼科尔森格式，并提到稳定性Finite-Difference Approximations to the Heat Equation by Gerald W. Recktenwald。

抱歉，我还没有添加评论的代表