0

私はMatlabに数十年のタイムスケールで使用していた1D熱拡散コードを持っており、現在、同じコードを使用して数百万年のスケールで作業しようとしています。明らかに、タイムステップを同じに保つと、計算に時間がかかりますが、タイムステップを増やすと、数値安定性の問題が発生します。

私の質問は次のとおりです。

この問題にどのように取り組むべきですか？最大安定タイムステップに影響を与えるものは何ですか？そして、これをどのように計算しますか？

どうもありがとう、

アレックス

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

FTCSのような明示的なスキームの安定条件は、$ r = K dt / dx ^ 2<1/2$または$dt<dx ^ 2 /（2K）$によって制御されます。ここで、Kは拡散係数です。これは、4次導関数の先行打ち切り誤差項の符号を負にするために必要です。

タイムステップによって制限されたくない場合は、暗黙的なスキームを使用することをお勧めします（ただし、明示的なスキームよりも計算コストが高くなります）。これは、前方オイラーの代わりに拡散項に後方オイラーを使用することで簡単に実現できます。もう1つのオプションは、暗黙的なCrank-Nicholsonです。

于 2013-03-28T19:43:28.680 に答える
0

@Isopycnal Oscillationは、最大の安定したステップが明示的なスキームで制限されているという点で完全に正しいです。参考までに、これは通常、離散フーリエ数または単にフーリエ数と呼ばれ、さまざまな境界条件を調べることができます。また、以下は、暗黙的またはクランク-ニコルソンスキームの導出に役立つ可能性があり、ジェラルドW.レックテンヴァルトによる熱方程式の安定性有限差分近似に言及しています。

申し訳ありませんが、コメントを追加する担当者はまだいません

于 2016-10-17T22:42:30.573 に答える