主要内容gydF4y2Ba

解偏微分方程gydF4y2Ba

本例通过使用符号数学工具箱™来求解微分方程来模拟海啸波现象。gydF4y2Ba

这个模拟是对这种现象的简化可视化,是基于Goring和Raichlen[1]的一篇论文。gydF4y2Ba

海啸模型的数学gydF4y2Ba

孤波(Korteweg-de Vries方程的孤子解)沿定深的管道以恒定速度从右向左传播。这相当于海啸在深海中行进。在运河的左端,有一个模拟大陆架的斜坡。到达斜坡后,孤波的高度开始增加。当水变得很浅时,大部分的波被反射回运河。然而,在斜坡的末端出现了一个狭窄但很高的水峰,并以减慢的速度沿着入射波的原始方向前进。这是海啸最终到达海岸,在海岸线上造成灾难性的破坏。接近海岸的波浪速度相对较小。海浪最终开始破裂。gydF4y2Ba

利用线性无色散水理论,得到高度gydF4y2Ba ugydF4y2Ba (gydF4y2Ba xgydF4y2Ba ,gydF4y2Ba tgydF4y2Ba )gydF4y2Ba 在不同深度的一维管道中,未受干扰的水位之上的自由表面波gydF4y2Ba hgydF4y2Ba (gydF4y2Ba xgydF4y2Ba )gydF4y2Ba 是以下偏微分方程的解。(参见[2])。gydF4y2Ba

ugydF4y2Ba tgydF4y2Ba tgydF4y2Ba =gydF4y2Ba ggydF4y2Ba (gydF4y2Ba hgydF4y2Ba ugydF4y2Ba xgydF4y2Ba )gydF4y2Ba xgydF4y2Ba

在这个公式中,下标表示偏导数,和gydF4y2Ba ggydF4y2Ba =gydF4y2Ba 9gydF4y2Ba .gydF4y2Ba 8gydF4y2Ba 1gydF4y2Ba 米gydF4y2Ba /gydF4y2Ba 年代gydF4y2Ba 2gydF4y2Ba 是重力加速度。gydF4y2Ba

考虑一个波穿过线性斜坡gydF4y2Ba hgydF4y2Ba (gydF4y2Ba xgydF4y2Ba )gydF4y2Ba 从一个深度恒定的区域gydF4y2Ba hgydF4y2Ba 2gydF4y2Ba 到一个深度恒定的区域gydF4y2Ba hgydF4y2Ba 1gydF4y2Ba ≪gydF4y2Ba hgydF4y2Ba 2gydF4y2Ba .傅里叶变换gydF4y2Ba tgydF4y2Ba 将水波偏微分方程转化为傅里叶模式的常微分方程gydF4y2Ba ugydF4y2Ba (gydF4y2Ba xgydF4y2Ba ,gydF4y2Ba tgydF4y2Ba )gydF4y2Ba =gydF4y2Ba UgydF4y2Ba (gydF4y2Ba xgydF4y2Ba ,gydF4y2Ba ωgydF4y2Ba )gydF4y2Ba egydF4y2Ba 我gydF4y2Ba ωgydF4y2Ba tgydF4y2Ba .gydF4y2Ba

-gydF4y2Ba ωgydF4y2Ba 2gydF4y2Ba UgydF4y2Ba =gydF4y2Ba ggydF4y2Ba (gydF4y2Ba hgydF4y2Ba UgydF4y2Ba xgydF4y2Ba )gydF4y2Ba xgydF4y2Ba

对于深度恒定的区域gydF4y2Ba hgydF4y2Ba 时,傅里叶模是沿相反方向匀速传播的行波gydF4y2Ba cgydF4y2Ba =gydF4y2Ba ggydF4y2Ba hgydF4y2Ba .gydF4y2Ba

ugydF4y2Ba (gydF4y2Ba xgydF4y2Ba ,gydF4y2Ba tgydF4y2Ba )gydF4y2Ba =gydF4y2Ba CgydF4y2Ba 1gydF4y2Ba (gydF4y2Ba ωgydF4y2Ba )gydF4y2Ba egydF4y2Ba 我gydF4y2Ba ωgydF4y2Ba (gydF4y2Ba tgydF4y2Ba +gydF4y2Ba xgydF4y2Ba /gydF4y2Ba cgydF4y2Ba )gydF4y2Ba +gydF4y2Ba CgydF4y2Ba 2gydF4y2Ba (gydF4y2Ba ωgydF4y2Ba )gydF4y2Ba egydF4y2Ba 我gydF4y2Ba ωgydF4y2Ba (gydF4y2Ba tgydF4y2Ba -gydF4y2Ba xgydF4y2Ba /gydF4y2Ba cgydF4y2Ba )gydF4y2Ba

解决方案gydF4y2Ba ugydF4y2Ba 2gydF4y2Ba (gydF4y2Ba xgydF4y2Ba ,gydF4y2Ba tgydF4y2Ba )gydF4y2Ba =gydF4y2Ba egydF4y2Ba 我gydF4y2Ba ωgydF4y2Ba (gydF4y2Ba tgydF4y2Ba +gydF4y2Ba xgydF4y2Ba /gydF4y2Ba cgydF4y2Ba 2gydF4y2Ba )gydF4y2Ba +gydF4y2Ba RgydF4y2Ba (gydF4y2Ba ωgydF4y2Ba )gydF4y2Ba egydF4y2Ba 我gydF4y2Ba ωgydF4y2Ba (gydF4y2Ba tgydF4y2Ba -gydF4y2Ba xgydF4y2Ba /gydF4y2Ba cgydF4y2Ba 2gydF4y2Ba )gydF4y2Ba 深水区是两个波的叠加:gydF4y2Ba

  • 以恒定速度向左移动的波gydF4y2Ba cgydF4y2Ba 2gydF4y2Ba =gydF4y2Ba ggydF4y2Ba hgydF4y2Ba 2gydF4y2Ba

  • 波向右传播的波,其振幅由频率相关的反射系数给出gydF4y2Ba RgydF4y2Ba (gydF4y2Ba ωgydF4y2Ba )gydF4y2Ba

这种选择gydF4y2Ba ugydF4y2Ba 2gydF4y2Ba 在深水区满足波动方程gydF4y2Ba RgydF4y2Ba (gydF4y2Ba ωgydF4y2Ba )gydF4y2Ba .gydF4y2Ba

解决方案gydF4y2Ba ugydF4y2Ba 1gydF4y2Ba (gydF4y2Ba xgydF4y2Ba ,gydF4y2Ba tgydF4y2Ba )gydF4y2Ba =gydF4y2Ba TgydF4y2Ba (gydF4y2Ba ωgydF4y2Ba )gydF4y2Ba egydF4y2Ba 我gydF4y2Ba ωgydF4y2Ba (gydF4y2Ba tgydF4y2Ba +gydF4y2Ba xgydF4y2Ba /gydF4y2Ba cgydF4y2Ba 1gydF4y2Ba )gydF4y2Ba 浅水区是一个以恒定速度向左传播的透射波gydF4y2Ba cgydF4y2Ba 1gydF4y2Ba =gydF4y2Ba ggydF4y2Ba hgydF4y2Ba 1gydF4y2Ba .这种选择gydF4y2Ba ugydF4y2Ba 1gydF4y2Ba 对于任意透射系数,在浅水区满足波动方程gydF4y2Ba TgydF4y2Ba (gydF4y2Ba ωgydF4y2Ba )gydF4y2Ba .gydF4y2Ba

对于过渡区域(斜率),使用gydF4y2Ba ugydF4y2Ba (gydF4y2Ba xgydF4y2Ba ,gydF4y2Ba tgydF4y2Ba )gydF4y2Ba =gydF4y2Ba UgydF4y2Ba (gydF4y2Ba xgydF4y2Ba ,gydF4y2Ba wgydF4y2Ba )gydF4y2Ba egydF4y2Ba 我gydF4y2Ba ωgydF4y2Ba tgydF4y2Ba .gydF4y2Ba

符号数学工具箱中海啸模型的参数金宝搏官方网站与解gydF4y2Ba

定义海啸模型的参数如下。忽略对频率的依赖gydF4y2Ba ωgydF4y2Ba 用下列表示法表示:gydF4y2Ba RgydF4y2Ba =gydF4y2Ba RgydF4y2Ba (gydF4y2Ba ωgydF4y2Ba )gydF4y2Ba ,gydF4y2Ba TgydF4y2Ba =gydF4y2Ba TgydF4y2Ba (gydF4y2Ba ωgydF4y2Ba )gydF4y2Ba ,gydF4y2Ba UgydF4y2Ba (gydF4y2Ba xgydF4y2Ba )gydF4y2Ba =gydF4y2Ba UgydF4y2Ba (gydF4y2Ba xgydF4y2Ba ,gydF4y2Ba ωgydF4y2Ba )gydF4y2Ba .gydF4y2Ba

信谊gydF4y2BalgydF4y2BaHgydF4y2BadepthratiogydF4y2BaggydF4y2Ba积极的gydF4y2Ba信谊gydF4y2BaxgydF4y2BatgydF4y2BawgydF4y2BaTgydF4y2BaRgydF4y2BaU (x)gydF4y2BaL1 = depthratio*L;L2 = l;h1 = depthratio*H;h2 = H;h(x) = x* h /L;C1 =√(g*h1);C2 =√(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(1我* w * (t + x / c2)) + R * exp(我* w * (t - x / c2));gydF4y2Ba

在线性斜坡上的过渡区域,使用gydF4y2BadsolvegydF4y2Ba来解傅里叶变换的ODEgydF4y2Ba UgydF4y2Ba 的gydF4y2Ba ugydF4y2Ba .gydF4y2Ba

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

解决方案gydF4y2Ba UgydF4y2Ba 是一个包含贝塞尔函数的复杂表达式。它包含两个任意的“常数”,依赖于gydF4y2Ba ωgydF4y2Ba .gydF4y2Ba

Const = setdiff(symvar(U), sym([depthratio,g,H,L,x,w]))gydF4y2Ba
Const =gydF4y2Ba
                  
                   
                    
                     
                      (gydF4y2Ba
                     
                      
                       
                        
                         
                          
                           
                            CgydF4y2Ba
                          
                          
                           
                            1gydF4y2Ba
                          
                         
                        
                       
                       
                        
                         
                          
                           
                            CgydF4y2Ba
                          
                          
                           
                            2gydF4y2Ba
                          
                         
                        
                       
                      
                     
                     
                      )gydF4y2Ba
                    
                   
                  

对于任何傅里叶模式,整体解必须是一个连续可微的函数gydF4y2Ba xgydF4y2Ba .因此,函数值和导数必须在接缝点匹配gydF4y2Ba lgydF4y2Ba 1gydF4y2Ba 而且gydF4y2Ba lgydF4y2Ba 2gydF4y2Ba .这提供了四个线性方程gydF4y2Ba TgydF4y2Ba ,gydF4y2Ba RgydF4y2Ba ,两个常数在gydF4y2Ba UgydF4y2Ba .gydF4y2Ba

Du1 (x) = diff(u1(x,0),x);Du2 (x) = diff(u2(x,0),x);dU(x) = diff(U(x),x);= [U(L1) == u1(L1,0), U(L2) == u2(L2,0),gydF4y2Ba...gydF4y2BadU(L1) == du1(L1), dU(L2) == du2(L2)];unknowns = [Const(1),Const(2),R,T];gydF4y2Ba

解这些方程。gydF4y2Ba

[Cvalue1, Cvalue2, R, T] = solve(eqs,未知);gydF4y2Ba

把结果代回gydF4y2Ba RgydF4y2Ba ,gydF4y2Ba TgydF4y2Ba ,gydF4y2Ba UgydF4y2Ba .gydF4y2Ba

U (x) =潜艇(U (x){常量(1)常量(2)},{Cvalue1, Cvalue2});gydF4y2Ba

你不能直接计算的解gydF4y2Ba ωgydF4y2Ba =gydF4y2Ba 0gydF4y2Ba 因为对应表达式的分子和分母都消失了。相反,找到这些表达式的低频限制。gydF4y2Ba

简化(限制(U(x), w, 0))gydF4y2Ba
ans =gydF4y2Ba

2gydF4y2Ba depthratiogydF4y2Ba +gydF4y2Ba 1gydF4y2Ba

简化(limit(R, w, 0))gydF4y2Ba
ans =gydF4y2Ba

-gydF4y2Ba depthratiogydF4y2Ba -gydF4y2Ba 1gydF4y2Ba depthratiogydF4y2Ba +gydF4y2Ba 1gydF4y2Ba

简化(limit(T, w, 0))gydF4y2Ba
ans =gydF4y2Ba

2gydF4y2Ba depthratiogydF4y2Ba +gydF4y2Ba 1gydF4y2Ba

这些限制非常简单。它们只取决于定义斜率的深度值的比值。gydF4y2Ba

用数值代替符号参数gydF4y2Ba

对于下面的计算,使用这些数值作为符号参数。gydF4y2Ba

  • 重力加速度:gydF4y2Ba ggydF4y2Ba =gydF4y2Ba 9gydF4y2Ba .gydF4y2Ba 8gydF4y2Ba 1gydF4y2Ba 米gydF4y2Ba /gydF4y2Ba 年代gydF4y2Ba egydF4y2Ba cgydF4y2Ba 2gydF4y2Ba

  • 运河深度:gydF4y2Ba HgydF4y2Ba =gydF4y2Ba 1gydF4y2Ba 米gydF4y2Ba

  • 浅部与深部深度比:gydF4y2Ba dgydF4y2Ba egydF4y2Ba pgydF4y2Ba tgydF4y2Ba hgydF4y2Ba rgydF4y2Ba 一个gydF4y2Ba tgydF4y2Ba 我gydF4y2Ba ogydF4y2Ba =gydF4y2Ba 0gydF4y2Ba .gydF4y2Ba 0gydF4y2Ba 4gydF4y2Ba

  • 斜率区域长度:gydF4y2Ba lgydF4y2Ba =gydF4y2Ba 2gydF4y2Ba 米gydF4y2Ba

G = 9.81;L = 2;H = 1;Depthratio = 0.04;h1 = depthratio*H;h2 = H;L1 = depthratio*L;L2 = l;C1 =√(g*h1);C2 =√(g*h2);gydF4y2Ba

定义入射孤子的振幅gydF4y2Ba 一个gydF4y2Ba 以匀速向左移动gydF4y2Ba cgydF4y2Ba 2gydF4y2Ba 在深水区。gydF4y2Ba

A = 0.3;孤子= @ (x, t) A *双曲正割(sqrt (3/4 * g * / H) * (x / c2 + t)) ^ 2;gydF4y2Ba

选择gydF4y2Ba NgydF4y2Ba tgydF4y2Ba 样本点gydF4y2Ba tgydF4y2Ba .时间尺度选择为入射孤子(时间)宽度的倍数。将傅里叶变换对应的离散频率存储在gydF4y2Ba WgydF4y2Ba .gydF4y2Ba

Nt = 64;TimeScale = 40*√(4/3*H/A/g);W =[0, 1:元/ 2 - 1,-(元/ 2 - 1):1]“* 2 *π/时间表;gydF4y2Ba

选择gydF4y2Ba NgydF4y2Ba xgydF4y2Ba 样本点gydF4y2Ba xgydF4y2Ba 每个区域的方向。创建样本点gydF4y2Ba XgydF4y2Ba 1gydF4y2Ba 对于浅水区,gydF4y2Ba XgydF4y2Ba 2gydF4y2Ba 对于深水区,和gydF4y2Ba XgydF4y2Ba 1gydF4y2Ba 2gydF4y2Ba 对于斜率区域。gydF4y2Ba

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

计算入射孤子在时间网格上的傅里叶变换gydF4y2Ba NgydF4y2Ba tgydF4y2Ba 等距采样点。gydF4y2Ba

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

基于傅立叶数据,构造了深水区行波解gydF4y2Ba 年代gydF4y2Ba .gydF4y2Ba

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

将深水区反射波的傅里叶模式转换为网格上的数值gydF4y2Ba (gydF4y2Ba xgydF4y2Ba ,gydF4y2Ba ωgydF4y2Ba )gydF4y2Ba 空间。将这些值与傅里叶系数相乘gydF4y2Ba 年代gydF4y2Ba 然后使用函数gydF4y2Ba传输线gydF4y2Ba来计算反射波gydF4y2Ba (gydF4y2Ba xgydF4y2Ba ,gydF4y2Ba tgydF4y2Ba )gydF4y2Ba 空间。请注意第一行的数值数据gydF4y2Ba RgydF4y2Ba 由NaN值组成,因为对符号数据进行了适当的数值计算gydF4y2Ba RgydF4y2Ba 为gydF4y2Ba ωgydF4y2Ba =gydF4y2Ba 0gydF4y2Ba 不可能。的第一行中定义值gydF4y2Ba RgydF4y2Ba 由于低频的限制。gydF4y2Ba

双(R =潜艇(潜艇(vpa(潜艇(R)), w w), x, X2));R(1) =双((1-sqrt (depthratio)) / (1 + sqrt (depthratio)));reflectedWave = real(S .* R .* exp(-1i*W*X2/c2)));gydF4y2Ba

对浅水区的透射波使用同样的方法。gydF4y2Ba

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

同样,对于斜率区域也可以使用这种方法。gydF4y2Ba

U12 =双(潜艇(潜艇(vpa(潜艇(U (x))), w w), x, X12));U12(1,:) = double(2/(1+根号下(depthratio)));U12 = real(ifft(S .* U12));gydF4y2Ba

规划解决方案gydF4y2Ba

为了更流畅的动画,使用三角插值沿绘图数据的列生成额外的样本点。gydF4y2Ba

孤子= interpft(孤子,10*Nt);reflectedWave = interpft(reflectedWave, 10*Nt);U12 = interpft(U12, 10*Nt);transmittedWave = interpft(transmittedWave, 10*Nt);gydF4y2Ba

创建在单独的图形窗口中显示的解决方案的动画情节。gydF4y2Ba

图(gydF4y2Ba“可见”gydF4y2Ba,gydF4y2Ba“上”gydF4y2Ba);情节([x_min, L1、L2 x_max], [h1, h1, h2, h2),gydF4y2Ba“颜色”gydF4y2Ba,gydF4y2Ba“黑”gydF4y2Ba)轴([x_min, x_max, -H-0.1, 0.6]gydF4y2Ba在gydF4y2Baline1 = plot(X1,transmittedWave(1,:),gydF4y2Ba“颜色”gydF4y2Ba,gydF4y2Ba“蓝”gydF4y2Ba);line12 = plot(X12,U12(1,:),gydF4y2Ba“颜色”gydF4y2Ba,gydF4y2Ba“蓝”gydF4y2Ba);line2 = plot(X2,孤子(1,:)+ reflectedWave(1,:)),gydF4y2Ba“颜色”gydF4y2Ba,gydF4y2Ba“蓝”gydF4y2Ba);gydF4y2Ba为gydF4y2BaT = 2:大小(孤子,1)*0.35线1。YData = transmittedWave(t,:);line12。YData = U12(t,:);么。YData =孤子(t,:) +反射波(t,:);暂停(0.1)gydF4y2Ba结束gydF4y2Ba

图中包含一个轴对象。axis对象包含4个line类型的对象。gydF4y2Ba

更多关于海啸gydF4y2Ba

在现实生活中,海啸的波长为数百公里,通常以每小时500公里以上的速度传播。(请注意,海洋的平均深度约为4公里,对应的速度为gydF4y2Ba ggydF4y2Ba hgydF4y2Ba ≈gydF4y2Ba 7gydF4y2Ba 0gydF4y2Ba 0gydF4y2Ba kgydF4y2Ba 米gydF4y2Ba /gydF4y2Ba hgydF4y2Ba ogydF4y2Ba ugydF4y2Ba rgydF4y2Ba )。在深海,振幅相当小,通常在0.5 m左右或更小。然而,当海啸传播到大陆架上时,它们的高度急剧增加:据报道,振幅可达30米或更高。gydF4y2Ba

一个有趣的现象是,尽管海啸通常以垂直于其行进方向的数百公里的波锋形式接近海岸线,但它们不会对海岸造成均匀的破坏。在某些地方它们会造成灾难,而在其他地方只观察到中波现象。这是由从海床到大陆架的不同坡度造成的。事实上,非常陡峭的斜坡会导致大部分海啸被反射回深水区,而小斜坡反射的波较少,从而传递出携带大量能量的窄而高的波。gydF4y2Ba

的不同值运行模拟gydF4y2Ba lgydF4y2Ba ,对应不同的斜率。斜坡越陡,传播的波就越低,强度越小。gydF4y2Ba

注意,这个模型忽略了色散和摩擦效应。在架子上,模拟失去了它的物理意义。在这里,摩擦效应很重要,造成了波浪的破碎。gydF4y2Ba

参考文献gydF4y2Ba

[1] Derek G. Goring和F. Raichlen,海啸-长波在大陆架上的传播,水道,港口,海岸和海洋工程杂志118(1),1992,pp. 41 - 63。gydF4y2Ba

[2] H.兰姆,流体力学,多佛,1932。gydF4y2Ba