ドキュメンテーション

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

偏微分方程式の求解

この例では、Symbolic Math Toolbox™ を使って微分方程式を解くことにより、津波現象をシミュレーションします。

このシミュレーションは、この現象を単純化して可視化したものであり、Goring と Raichlen [1] の論文に基づいています。

津波モデルの数学

孤立波 (コルトヴェーグ・ド・フリース (KdV) 方程式のソリトン解) が、一定の水深の運河に沿って右から左へ一定速度で進行します。これは、深海上を津波が進行するのに対応します。運河の左端には、大陸棚をシミュレーションした勾配があります。孤立波は、勾配に到達した後、高さを増し始めます。水深が非常に浅くなると、ほとんどの波は反射して運河に戻ります。しかし、勾配の端で水の尖端は狭く、高くなり、速度を落として入射波の元の方向に進みます。これが、最終的に陸地を襲う津波で、海岸線に沿って壊滅的な破壊を引き起こします。陸地近くの波の速度は比較的低速です。波は最終的には砕け始めます。

線形無分散波理論を使用すると、深さ h(x) が変化する 1 次元の運河における静水位を上回る自由表面波の高さ u(x,t) は、次の偏微分方程式の解になります。([2] を参照してください)。

utt=g(hux)x

この式で、添字は偏導関数を示します。g=9.81m/s2 は重力加速度です。

一定の深さ h2 の領域から一定の深さ h1h2 の領域への線形勾配 h(x) を通過する波を考えます。t についてのフーリエ変換により、水の波の偏微分方程式が次のようにフーリエ モード u(x,t)=U(x,ω)eiωt の常微分方程式に変換されます。

-ω2U=g(hUx)x

一定の深さ h の領域では、フーリエ モードは逆方向に一定速度 c=gh で伝播する進行波です。

u(x,t)=C1(ω)eiω(t+x/c)+C2(ω)eiω(t-x/c)

水深が深い領域の解 u2(x,t)=eiω(t+x/c2)+R(ω)eiω(t-x/c2) は、次の 2 つの波の重ね合わせです。

  • 一定速度 c2=gh2 で左に進行する波

  • 周波数に依存する反射係数 R(ω) によって与えられる振幅で右に進行する波

この u2 の選択は、どの R(ω) でも水深が深い領域における波動方程式を満たします。

水深が浅い領域の解 u1(x,t)=T(ω)eiω(t+x/c1) は、一定速度 c1=gh1 で左に進行する伝送波です。この u1 の選択は、どの伝送係数 T(ω) でも水深が浅い領域における波動方程式を満たします。

遷移領域 (勾配) では、u(x,t)=U(x,w)eiωt を使用します。

Symbolic Math Toolbox における津波モデルのパラメーターと解

津波モデルのパラメーターを次のように定義します。次の表記法 R=R(ω)T=T(ω)U(x)=U(x,ω) における周波数 ω への依存関係は無視します。

syms L H depthratio g positive
syms x t w T R U(x)

L1 = depthratio*L; 
L2 = L;

h1 = depthratio*H; 
h2 = H;
h(x) = x*H/L;

c1 = sqrt(g*h1);
c2 = sqrt(g*h2);

u(x,t)  = U(x)*exp(1i*w*t);
u1(x,t) = T*exp(1i*w*(t + x/c1));
u2(x,t) = exp(1i*w*(t + x/c2)) + R*exp(1i*w*(t - x/c2));

線形勾配の遷移領域では、dsolve を使用して、u のフーリエ変換 U の ODE を解きます。

wavePDE(x,t) = diff(u,t,t) - g*diff(h(x)*diff(u,x),x);
slopeODE(x) = wavePDE(x,0); 
U(x) = dsolve(slopeODE);

U は、ベッセル関数を含む複雑な式です。これは、ω に依存する 2 つの任意 "定数" を含みます。

Const = setdiff(symvar(U), sym([depthratio,g,H,L,x,w]))
Const = (C3C4)

どのフーリエ モードでも、解全体が連続的に微分可能な x の関数でなければなりません。したがって、関数値と導関数は、シーム点 L1 および L2 で一致しなければなりません。これにより、TR、および U の 2 つの定数に対する 4 つの線形方程式が得られます。

du1(x) = diff(u1(x,0),x);
du2(x) = diff(u2(x,0),x);
dU(x)  = diff(U(x),x);

eqs =  [ U(L1) == u1(L1,0), U(L2) == u2(L2,0),...
        dU(L1) == du1(L1), dU(L2) == du2(L2)];
unknowns = [Const(1),Const(2),R,T];

これらの方程式を解きます。

[Cvalue1, Cvalue2, R, T] = solve(eqs, unknowns);

結果を RT および U に代入します。

U(x) = subs(U(x), {Const(1),Const(2)}, {Cvalue1,Cvalue2});

ω=0 の場合、対応する式の分子と分母の両方が消えるため、解を直接評価することはできません。代わりに、これらの式の低周波数極限を求めます。

simplify(limit(U(x), w, 0)) 
ans = 

2depthratio+1

simplify(limit(R, w, 0)) 
ans = 

-depthratio-1depthratio+1

simplify(limit(T, w, 0))
ans = 

2depthratio+1

これらの極限は非常にシンプルです。これらは、勾配を定義する深さの値の比にのみ依存します。

シンボリック パラメーターへの数値の代入

次の計算では、これらの数値をシンボリック パラメーターで使用します。

  • 重力加速度: g=9.81m/sec2

  • 運河の深さ: H=1m

  • 浅い領域と深い領域の深さの比: depthratio=0.04

  • 勾配領域の長さ: L=2m

g = 9.81;
L = 2;
H = 1;
depthratio = 0.04; 

h1 = depthratio*H;
h2 = H;

L1 = depthratio*L;
L2 = L; 

c1 = sqrt(g*h1);
c2 = sqrt(g*h2);

水深が深い領域を左へ一定速度 c2 で進行する振幅 A の入射ソリトンを定義します。

A = 0.3;
soliton = @(x,t) A.*sech(sqrt(3/4*g*A/H)*(x/c2+t)).^2;

t について Nt 個の標本点を選択します。時間スケールは、入射ソリトンの (時間) 幅の倍数として選択されています。対応するフーリエ変換の離散化周波数を W に格納します。

Nt =  64;
TimeScale = 40*sqrt(4/3*H/A/g);
W = [0, 1:Nt/2 - 1, -(Nt/2 - 1):-1]'*2*pi/TimeScale;

各領域の x 方向の標本点 Nx 個を選択します。標本点 X1 個を水深が浅い領域に、X2 個を水深が深い領域に、X12 個を勾配領域に作成します。

Nx = 100;
x_min = L1 - 4;
x_max = L2 + 12;
X12 = linspace(L1, L2, Nx);
X1  = linspace(x_min, L1, Nx);
X2  = linspace(L2, x_max, Nx);

等間隔の標本点 Nt 個の時間グリッドにおける入射ソリトンのフーリエ変換を計算します。

S = fft(soliton(-0.8*TimeScale*c2, linspace(0,TimeScale,2*(Nt/2)-1)))';
S = repmat(S,1,Nx);

S のフーリエ データに基づいて、水深が深い領域の進行波の解を作成します。

soliton = real(ifft(S .* exp(1i*W*X2/c2)));

水深が深い領域の反射波のフーリエ モードを (x,ω) 空間のグリッド上の数値に変換します。これらの値に S のフーリエ係数を乗算し、関数 ifft を使用して (x,t) 空間の反射波を計算します。ω=0 についてシンボリック データ R の適切な数値評価が不可能なため、数値データ R の 1 行目は NaN 値で構成されることに注意してください。R の 1 行目の値を低周波数極限として定義します。

R = double(subs(subs(vpa(subs(R)), w, W), x ,X2));
R(1,:) = double((1-sqrt(depthratio)) / (1+sqrt(depthratio)));
reflectedWave = real(ifft(S .* R .* exp(-1i*W*X2/c2)));

同じ方法を水深が浅い領域の伝送波に対して使用します。

T = double(subs(subs(vpa(subs(T)),w,W),x,X1));
T(1,:) = double(2/(1+sqrt(depthratio)));
transmittedWave = real(ifft(S .* T .* exp(1i*W*X1/c1)));

また、この方法を勾配領域に対して使用します。

U12 = double(subs(subs(vpa(subs(U(x))),w,W),x,X12));
U12(1,:) = double(2/(1+sqrt(depthratio)));
U12 = real(ifft(S .* U12));

解のプロット

アニメーションをより滑らかにするには、プロット データの列に沿って三角内挿を使用して追加の標本点を生成します。

soliton = interpft(soliton, 10*Nt);
reflectedWave = interpft(reflectedWave, 10*Nt);
U12 = interpft(U12, 10*Nt);
transmittedWave = interpft(transmittedWave, 10*Nt);

解のアニメーション プロットを作成し、別の Figure ウィンドウに表示します。

figure('Visible', 'on');
plot([x_min, L1, L2, x_max], [-h1, -h1, -h2, -h2], 'Color', 'black')
axis([x_min, x_max, -H-0.1, 0.6])
hold on

line1 = plot(X1,transmittedWave(1,:), 'Color', 'blue');
line12 = plot(X12,U12(1,:), 'Color', 'blue');
line2 = plot(X2,soliton(1,:) + reflectedWave(1,:), 'Color', 'blue');

for t = 2 : size(soliton, 1)*0.35
  line1.YData = transmittedWave(t,:);
  line12.YData = U12(t,:);
  line2.YData = soliton(t,:) + reflectedWave(t,:);
  pause(0.1)
end

津波の詳細

現実の世界では、津波は何百 km もの波長をもち、多くの場合、時速 500 km 以上で進みます (海洋の平均的な深さは約 4 km で、速度 gh700km/hour に対応します)。深海上では振幅はかなり小さく、たいていの場合およそ 0.5 m 以下です。しかし、津波が大陸棚を伝播する場合には、その高さは劇的に増加します。30 m 以上の振幅が報告されたことがあります。

興味深い現象の 1 つは、津波は通常、進行方向に垂直な何百 km にも広がる波面として海岸線に押し寄せますが、海岸に沿って一様の被害を与えるわけではないことです。大きな被害を与える地点がある一方、穏やかな波動現象しか観測されない場所もあります。これは、海底から大陸棚への勾配が異なるために引き起こされます。実際、非常に急な勾配の場合、ほとんどの津波は水深が深い領域に反射されます。一方、勾配が小さい場合、波はあまり反射されず、大きなエネルギーをもつ狭くて高い波が伝送されます。

さまざまな勾配に対応するさまざまな L の値でシミュレーションを実行します。勾配が急になるほど、伝送される波は低く、エネルギーは小さくなります。

このモデルでは分散と摩擦の効果を無視していることに注意してください。大陸棚では、このシミュレーションは物理的意味を失います。ここでは、砕波を引き起こす摩擦効果が重要です。

参考文献

[1] Derek G. Goring and F. Raichlen, Tsunamis - The Propagation of Long Waves onto a Shelf, Journal of Waterway, Port, Coastal and Ocean Engineering 118(1), 1992, pp. 41 - 63.

[2] H. Lamb, Hydrodynamics, Dover, 1932.