Timestep stability in a 1D heat diffusion model

Hi,

I have a 1D heat diffusion code which I was using on a timescale of 10s of years and I am now trying to use the same code to work on a scale of millions of years. Obviously if I keep my timestep the same this will take ages to calculate but if I increase my timestep I encounter numerical stability issues.

My questions are:

How should I approach this problem?
What affects the maximum stable timestep? And how do I calculate this?

Many thanks,

Alex

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;
endHi,

My questions are:

How should I approach this problem?
What affects the maximum stable timestep? And how do I calculate this?

Many thanks,

Alex

close all
clear all

%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;

%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);

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

drawnow

t=t+dt;
end Hi,

My questions are:

How should I approach this problem?
What affects the maximum stable timestep? And how do I calculate this?

Many thanks,

Alex

close all
clear all

%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;

%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);

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

drawnow

t=t+dt;
end time step; heat diffusion MATLAB Answers — New Questions

Cart

Cart