0

我在 Matlab 中有一个 1D 热扩散代码，我在 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
4

2 回答 2

1

The stability condition for an explicit scheme like FTCS is governed by $r = K dt/dx^2 < 1/2$ or $dt < dx^2/(2K)$ where K is your coefficient of diffusion. This is required in order to make the sign of the 4th order derivative leading truncation error term be negative.

If you do not want to be limited by timestep I suggest using an implicit scheme (albeit at a higher of computational cost than an explicit scheme). This can be achieved simply by using backward Euler for the diffusion term instead of forward Euler. Another option is Crank-Nicholson which is also implicit.

于 2013-03-28T19:43:28.680 回答
0

@Isopycnal 振荡是完全正确的，因为最大稳定步长在显式方案中受到限制。仅供参考，这通常称为离散傅立叶数或只是傅立叶数，可以针对不同的边界条件进行查找。以下内容也可以帮助您推导隐式或 Crank-Nicholson 方案，并提到Gerald W. Recktenwald 的热方程的稳定性有限差分近似。

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

于 2016-10-17T22:42:30.573 回答