基于COMSOL Multiphysics求解KdV方程

糜凱華

（昆明理工大學(xué) 電力工程學(xué)院，云南 昆明 650500）

很多物理現(xiàn)象都可以抽象成數(shù)學(xué)模型，建立起相應(yīng)的偏微分方程即數(shù)學(xué)物理方程.然后對這些偏微分方程組求解得到真實世界的描述.Korteweg-deVires方程式由Korteweg和Vires于1895年對水波建模得到的，方程不考慮消耗并且波是永遠(yuǎn)運動的.孤立波是一種傳播過程中形狀、幅度和速度都保持不變的脈沖行波.它有內(nèi)似粒子的性質(zhì)，也叫孤立子.從物理學(xué)的觀點看，孤立波是物質(zhì)非線性效應(yīng)的一種特殊產(chǎn)物，而從數(shù)學(xué)上看，它是某些非線性偏微分方程的一類穩(wěn)定的、能量有限的非彌散解.孤立波在相互碰撞后，仍能保持各自的形狀和速度不變.不僅在淺水波中有孤立波現(xiàn)象，在光纖通信、金屬相變和神經(jīng)傳播等許多領(lǐng)域都有孤立波現(xiàn)象，即某種現(xiàn)象或信息以幾乎恒定的脈沖形態(tài)傳播.由于其具有這種特殊的性質(zhì)，所以在等離子物理學(xué)、高能電磁學(xué)、流體力學(xué)和非線性光學(xué)等領(lǐng)域中得到廣泛應(yīng)用.由于非線性科研的發(fā)展，求解非線性方程的方法層出不窮，如：有限元法[1]、影射法[2]、非線性變換法[3]、同倫分析法[4]、雙Jacobi橢圓函數(shù)展開法[5]、齊次平衡法[6]、試探函數(shù)法[7].而COMSOL Multiphysics的優(yōu)點是求解偏微分方程或方程組，本文利用基于有限單元法的COMSOL Multiphysics有限元軟件對KdV方程進(jìn)行數(shù)值求解.

1 KdV方程及仿真模型的建立

1.1 KdV方程的建立

本文研究的KdV方程由Korteweg和他的學(xué)生deVires于1895年對水波建模得到.KdV方程經(jīng)過變換可得到各種不同的形式.為了方便，本文取KdV方程表達(dá)式為

此KdV方程是一個典型的一維三階偏微分方程.假設(shè)該孤立波在Ω=[ ]-8,+8的范圍內(nèi)周期性變化，其表達(dá)式為

初始條件的表達(dá)式為

以上三式即構(gòu)成一個初邊值問題.式（1）的解析解為

將初始條件式（3）與式（4）比較，可知兩式并不一致.因此，就會在初始時刻產(chǎn)生兩個迭加在一起的不同振幅和速度的孤立波.

1.2 仿真模型的建立

根據(jù)上述所建立的KdV方程，采用基于有限單元法的COMSOL Multiphysics軟件建立有限元模型，其控制方程（也即物理場）選擇Δu數(shù)學(xué)模塊中Δu偏微分方程接口.因通常情況下，常見的有限元軟件只能計算二階偏微分方程或偏微分方程組，在COMSOL Multiphysics中，常規(guī)的應(yīng)用也是二階方程表述的模式.所以如何求解高階偏微分方程是一個很棘手的問題.值得注意的是COMSOL Multiphysics軟件具有強耦合[8]的計算功能，可以為高階偏微分方程提供一種合理的求解方法.

對于我們所研究的三階KdV方程，可以向方程中引入一些中間變量，通過變量代換的方式進(jìn)行降階，從一個三階方程降階得到一組存在耦合關(guān)系的二階偏微分方程組，減少三階偏微分方程求解的難度.為此，在選擇控制方程時需要輸入因變量個數(shù)為2，分別定義為u1和u2，并令中間變量u2=uxx，然后將u改為u1，即將所建立的KdV方程修改為一個偏微分方程組相應(yīng)的初始條件也變?yōu)槿缦路匠探M

周期性邊界條件變?yōu)槿缦卤磉_(dá)式為

1.3 幾何模型和有限元模型

所研究的問題屬于一維情況，故在COM?SOL Multiphysics中建立的幾何模型即是KdV方程中求解域Ω=[ ]-8,+8表示的一條線段.為了降低計算機的計算成本，加快模型的求解速度，同時也保證有一定的精度，劃分求解域的最大單元尺寸為0.1，單元的劃分總數(shù)為160個單元，其幾何模型和有限元網(wǎng)格如圖1和圖2所示.

圖1 KdV方程的幾何模型

圖2 KdV方程的有限元模型

1.4 初始條件和邊界條件的設(shè)定

在COMSOL Multiphysics中直接將KdV方程的初始值式（6）輸入初始值對話框.而周期性邊界條件則需要選擇幾何模型中線段的兩個端點，該條件保證兩端邊界解的連續(xù)性.所設(shè)置的初始值和周期性邊界條件對應(yīng)于上述推導(dǎo)得到的新的KdV偏微分方程組和定解條件，分別如式（5）、式（6）和式（7）所示.其控制方程描述如下：

2 數(shù)值結(jié)果與分析

通過上述所建立的有限元模型，在設(shè)定瞬態(tài)求解器的求解時間范圍為0～2 s，時間間隔為0.025 s的條件下進(jìn)行全局計算.通過求解，得到指定時間段內(nèi)的波形分布，如圖 3所示，分別表示第 0.0 s、0.3 s、0.7 s、1.1 s和1.2 s時的波形.為了更全面地了解碰撞過程，將所有時刻的波形分布用三維圖的形式表示，如圖4所示.

此外，為了更為清晰地分析數(shù)值模擬得出的孤立波傳播特性是否與行波解所表達(dá)的物理意義相符.本文另取孤立波的形式為

周期性邊界條件與式（2）一樣.式（9）的孤波解為

圖3 孤立波在時間t=0 s、0.3 s、0.7 s、1.1 s、1.2 s的分布

式中c和x0是由初始條件決定的常數(shù)，x0為孤波的初始位置，為傳播速度.

從圖3可以清晰地看出：兩個不同速度和振幅的孤子相互碰撞，相互穿過，相互之間沒有影響.而且該孤立波不會發(fā)生消散.孤立波向右傳播過程中，大波追趕小波，從而引起碰撞.

從圖4可以看出：孤立子在一定的空間內(nèi)會發(fā)生碰撞和重現(xiàn)，兩個迭加在一起的孤立波以不同的速度傳播，速度快的波穿越速度慢的波.

從圖5和圖6的比較中可發(fā)現(xiàn)，圖5中孤立波的傳播速度小于圖6，即c越大，移動速度越快，波高越大，但寬度越窄.顯然這和行波解表達(dá)的物理意義是相符的，即波高為c 2，波寬反比于.由此不難發(fā)現(xiàn)本文對KdV方程的數(shù)值模擬方法是正確的.

圖4 孤立波波形在所有時刻的三維表示

3 結(jié)論

本文針對一類三階KdV方程通過基于有限元法的COM?SOL Multiphysics軟件進(jìn)行求解，從數(shù)值計算結(jié)果可仿真模擬出孤立波的傳播特性.有限單元法在求解偏微分方程中已有廣泛的應(yīng)用，傳統(tǒng)方法都是建立數(shù)學(xué)模型進(jìn)行精確求解，計算量大而且比較復(fù)雜.相比之下采用有限元法進(jìn)行數(shù)值求解可以在很大程度上提高科研的效率，節(jié)省大量的人力、物力.對于高階偏微分方程在COMSOL Multiphysics軟件中很容易通過變量代換進(jìn)行降階的方式來求解.通過對孤立波的數(shù)值模型仿真分析證實本文的數(shù)值模擬方法是有效的.

圖5 孤立波在時間t=0.0 s、0.4 s的分布

圖6 孤立波在時間t=0.0 s、0.4 s的分布

[1]劉興霞，孫建安.四次樣條Galerkin有限元法求解KdV方程[J].西北師范大學(xué)學(xué)報，2008，44（4）：66-70.

[2]岳萍，龔倫訓(xùn).用射影法求解KdV方程[J].大學(xué)物理，2010，29（2）：14-16.

[3]Otwinow ski M,Paul R,Laidlaw W G.Exact tavelling wave solutions o f a class of nonlinear diffusion equations by reduction to a quadrature[J].Phys Lett,1988,A128:483.

[4]楊紅娟，石玉仁.用同倫分析法求解KdV方程的孤立波解[J].西北師范大學(xué)學(xué)報，2009，45（5）：39-43.

[5]錢天虹，劉中飛.雙Jacobi橢圓函數(shù)展開法及其對KdV方程求解[J].安徽教育學(xué)院學(xué)報，2004，22（3）：5-8.

[6]Fan E G,Zhang H Q.The homogeneo us balance method for solving nonlinear soliton equations[J].Acta Phys Sin,1998,47(3):53.

[7]Kudry ashow N A.Ex act so lutions of the generalized Kuramo to Sivashinsky equation[J].Phys Lett,1990,A 147:287.

[8]William B J,Zimmerman，中仿科技公司.COMSOL Multiphysics有限元法多物理場建模與分析[M].北京：人民交通出版社.