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 <2分之1$或$ DT < DX管辖^ 2 /（2K）$其中K是你的扩散系数。这是为了使4阶微分前导截断误差项的符号为负值所必需的。

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

Answer 2

@Isopycnal振荡是完全正确的，因为在明确的方案中最大稳定步骤是有限的。仅供参考，通常称为离散傅立叶数或仅为傅里叶数，可查找不同的边界条件。 也可以帮助你推导隐式或曲柄尼科尔森方案，并提到稳定性Finite-Difference Approximations to the Heat Equation by Gerald W. Recktenwald。

对不起，我没有代表尚未添加注释