移动载荷作用下浮冰的非线性动力学响应1)

孟洋涵 王展 ,**,2)

* (中国科学院力学研究所流固耦合实验室,北京 100190)

† (中国科学院大学工程科学学院,北京 100049)

** (中国科学院大学未来技术学院,北京 100049)

引言

作为经典水波理论的扩展,水弹性波考虑了自由面的弹性效应,可用于描述弹性薄板与流体之间的相互作用.对水弹性波的研究起源于对高纬度地区大型冰盖的利用,这些冰盖常被作为季节性的交通通道供车辆行驶及小型飞机起飞降落.当车辆或小型飞机在冰盖上行驶时,移动的压力源会激发一系列在冰层和流体界面处传播的波,称之为水弹性波,其涉及的回复力包括重力和冰层变形对流体施加的弹性.与重力毛细波不同,水弹性波的波长可达到数十米甚至上百米,因此在实验中更易观测.Takizawa[1]和Squire等[2]分别在Lake Saroma 和McMurdo Sound 对移动载荷作用下的冰层响应进行了测量,得到了不同运动速度下的冰层挠度及应变.

为了更好地理解运动载荷下的冰层响应问题,许多研究者从理论上对水弹性波进行了细致的研究,但早期的研究多基于线性理论.Kheisin[3]在研究点载荷及线载荷作用下的冰层位移时发现线载荷存在两个临界速度,在临界速度下冰层挠度趋于无穷.Nevel[4]扩展了Kheisin 的分析方法以处理分布在圆形区域上的均匀载荷,发现点载荷也存在对应的临界速度.进一步的分析表明这个临界速度就是自由水弹性波的最小相速度[5].接着文献[6]运用渐近傅里叶展开研究了匀速运动载荷激发的远场稳定波形,发现临界速度恰巧为水弹性波的群速度,由此给出了临界速度处冰层响应趋于无穷的物理解释.与此同时,他们还在不同的系统参数下得到了各类复杂的波形包括焦散和零波响应区.Babaei等[7]以及van der Sanden 和Short[8]利用卫星观测了车辆在冰面上行驶时所激发的波形,观测结果验证了文献[6]对远场波形的理论预测.Schulkes等[9]在文献[6]工作的基础上,进一步考虑了冰层压应力,均匀流以及流体分层对冰层响应的影响.对于临界速度处所出现的奇性,Kheisin[10]意识到可以从瞬态冰层动力学响应入手,分析结果表明,当线载荷以最小相速度cmin运动时水弹性波振幅与t1/2成正比,t为时间参数.Schulkes 和Sneyd[11]则在研究匀速线载荷作用下的时间依赖响应时得到了第二个临界速度(g为重力加速度,H为水深),当载荷以此速度运动时水弹性波振幅与t1/3成正比.然而Nugroho等[12]的研究结果却发现在三维情况下,当载荷运动速度为时仍然可以得到稳态有界波形,并不会出现水弹性波振幅无限增长的情况,因此并不能称之为临界速度.Miles 和Sneyd[13]探究了加速载荷作用下的冰层响应,研究结果显示将载荷从静止逐渐加速到两个临界速度可以避免冰层响应在临界速度处的奇异.

在对冰层动力学响应的实验研究中,Wilson[14]最早发现了观测点处冰层最大垂直位移对应的时刻与载荷经过观测点的时刻之间的差异.随后Beltaos[15],Takizawa[1]以及Squire等[2]在他们的实验中证实了这种滞后效应的存在.同时Takizawa[1]在控制方程中引入耗散项对实验中出现的滞后效应进行了解释,但其所考虑的仅是稳态冰层响应而无法得到冰层响应的时历曲线.Hosking等[16],Wang等[17]利用记忆函数在线性理论中引入黏弹性得到了临界速度V=cmin处的有界冰层响应以及滞后行为.

黏性项的引入能够合理地解释实验中所出现的滞后效应,但利用线性理论对临界速度附近的大振幅冰层动力学响应进行探究仍然存在较大的局限,因此非线性对于描述移动载荷作用下的冰层效应同样重要.Părău 和Dias[18]在考虑线性理论中的两个共振(奇点)情况时首次引入非线性效应并得到了临界速度附近的有界冰层位移.Dinvay等[19]考虑非线性、黏性以及惯性效应,利用Dirichilet-Neumman 算子建立了能够描述各类运动载荷下冰层响应的完全色散弱非线性模型并通过数值计算结果与实验结果之间的对比验证了模型的有效性.但是需要指出的是由于非线性项仅保留至二阶,该模型对应的非线性薛定谔方程是不准确的.在利用二阶模型计算临界速度附近的冰层动力学响应时所得到的数值结果可能会与实验结果存在较大差异.除此之外,他们所关注的实验均为水深较小的情况,而未对Squire等[2]在深水中的实验进行讨论.因此仍需要发展更为准确的非线性黏弹性理论来描述不同水深情况下临界速度cmin附近的冰层大幅度挠曲.

通过对相关的拟微分算子进行展开并将非线性项保留到三阶,本文将完全非线性二维问题转化为仅与自由面上变量相关的一维系统,即三阶截断模型.为验证三阶截断模型的准确性及其相对于二阶模型[19]的优势,本文从理论和数值上对波包型孤立波解进行了探究.理论上不考虑黏性和外加载荷的作用,采用多重尺度方法推导三阶非线性Schrödinger 方程,基于此方程预测不同水深下孤立波解的存在性以及三阶截断模型的准确性.数值上仍然暂不考虑黏性和外加载荷,计算完全欧拉方程、三阶截断模型以及二阶截断模型在两个典型水深下的孤立波解及分岔曲线.数值结果表明,三阶截断模型能够较好地与欧拉方程吻合,其精度远高于二阶截断模型.最后将黏性和外加载荷考虑在内,基于发展的三阶截断模型对两个工况中移动载荷作用下的浮冰动力学响应进行数值模拟并将数值结果与实验记录进行对比.

1 数学模型及数值方法

考虑密度为 ρ,水深为h的不可压缩、无黏流体,其顶部自由面覆盖密度为 ρi,厚度为d的无限大冰层,可通过弹性形变对流体施加回复力.假设冰层在受到移动载荷作用前没有受到压缩或拉伸.流体底部为刚性边界,不可穿透.建立二维笛卡尔坐标系,使x轴与未发生形变时的冰层重合,y轴垂直向上,与重力加速度g的方向相反,冰层位移为 η (x,t) .流体运动无旋,因此引入速度势函数 ϕ,其满足Laplace方程

自由面上的运动学及动力学条件为

其中,P为冰层弹性变形引起的回复力,p(x,t) 表示施加的外力载荷.采用Toland[20]提出的水弹性模型,同时考虑冰层的惯性及黏性,于是自由面上的法向应力平衡方程为

其中Pa为大气压,恒为常数,为冰层的弹性刚度,κ 为冰层曲率,s为弧长参数,黏性参数b＞0.本文将冰层视为各向同性材料,一方面是由于各向异性材料所涉及的弹性参数过于复杂,在建立数学模型描述法向应力突变时会存在较大的困难;另一方面,Takizawa[1]在实验中测量了冰层在横向和纵向上的弹性模量,两者之间差异不大.Squire等[2]也指出真实弹性模量与冰层响应并不直接相关,更为重要的是实验中的有效弹性模量,因此本文采用了简化的各向同性假设[6,19,21,22]并使用了实验测量得到的有效弹性模量.此外底部边界处满足不可穿透条件 ϕy=0,y=-h.

1.1 三阶截断模型

为了计算载荷作用下的冰层响应并将数值结果与已有的实验结果进行对比,引入Dirichilet-Neumman (DtN)算子对原始欧拉方程进行近似处理.利用该算子可避免求解自由边界上的Laplace 方程,将自由面处的边界条件转化为以正则变量表达的形式,消去边界条件对y坐标的依赖,将原始的二维问题转换为一维系统.同时可对惯性项所引入的二阶时间导数项进行处理,便于进行后续的时间依赖计算.令自由面处的速度势函数 ξ (x,t)=ϕ(x,η(x,t),t),定义DtN 算子

Craig 和Sulem[23]已经证明当 η 的范数小于一个确定的值,G (η) 是解析的并且可以展开为级数形式

其前三项可写为

其中 D=(-∂xx)1/2.深水情况下 G0=D,其他两项形式不变.由此,惯性项 ηtt可用DtN 算子表示

这里保留至关于 ξ 的三阶项.将用DtN 算子表示的各项代入自由面的运动学和动力学条件,可得

算子F 定义为

将1 +αG0记为K.对算子F 取逆算子可得

则动力学边界条件可进一步写为

这里将式(5)和式(7)称为三阶截断模型.

1.2 数值方法

本文将通过拟谱法在周期域上对三阶截断模型进行数值求解,所有的导数项均在傅里叶空间中计算,而非线性项则在物理空间中进行计算.首先对方程式(5)和式(7)进行傅里叶变换

由于 ξ 和 η 均为实数,关于p,q的两方程实际上是等价的,因此对三阶截断模型的求解可以简化为对下式的求解

在数值求解方程式(8) 后,可以进一步得到ξ和η

2 模型验证及数值结果

2.1 孤立波解及模型验证

本节从自由波包型孤立波解入手,通过比较完全欧拉方程、三阶截断模型和二阶截断模型的孤立波解,验证三阶截断模型的准确性,说明三阶截断模型相比于二阶截断模型[19]更为合理.在弱非线性理论中,波包型孤立波存在的条件是:群速度与相速度在非色散点相等且在此波数下对应波包方程为焦聚型.在不考虑惯性项的情况下,水弹性波所对应的三阶非线性薛定谔方程(NLS) 的性质会在临界水深hc=233发生变化[22].当水深小于临界深度时,NLS是焦聚型,波包型孤立波从振幅无限小的周期波分岔而来.而当水深大于临界深度时,NLS 变为焦散型,此时只存在有限振幅的波包型孤立波解.为了更好地理解自由波包型孤立波解的存在性并对三阶截断模型的准确性进行验证,首先令三阶截断模型中的黏性系数 β=0,在无外力作用下推导考虑冰层惯性的三阶NLS 方程.对弱非线性水弹性波,假设η～O(ε),ϕ～O(ε),这里 ε 是用来衡量波陡的小参数.由此,自由面上的势函数可在y=0 处展开为

因此自由面上的运动学和动力学条件在y=0 处展开为

为了得到波包的控制方程,引入慢变量X=εx,T=εt和 τ=ε2t.设解的形式为

这里 θ=kx-ωt,c.c.代表复共轭,i 为虚数单位.将解式(9)和式(10)代入运动学和动力学边界条件中,在O(ε) 阶得到

从上式可知 φ11=ϕ11/cosh(kh) .同时要使上述方程存在非零解需要满足可解性条件

上式为水弹性波的色散关系式.

在O(ε2) 阶,可以得到

其中

为水弹性波的群速度.在O(ε3) 阶时,由可解性条件经过繁复的计算最终可以得到波包的控制方程,即三阶非线性薛定谔方程

D0,D1和D2的具体表达式见附录.为了后文讨论方便,分别将色散项和非线性项的系数记作 λ,γ .令惯性项系数为0,得到的三阶NLS 与Milewski 和Wang[25]所推导的方程退化到二维情形下的形式是一致的.考虑深水情况h→∞,则此时色散项和非线性项的系数变为

需要指出的是,二阶截断模型[19]保留了原始欧拉方程的色散关系,对应三阶NLS 方程的色散项系数与本文所推导的是一致的,但其在对原始方程进行近似时仅将非线性项保留至二阶,因此对应的NLS 中非线性项的系数是不准确的.

为了验证三阶截断模型的准确性,本文重点关注了水深h=6.8 m和h=250 m 时的NLS 方程性质及孤立波解.在不考虑惯性项影响时,水深h=6.8 m对应焦聚型NLS,水深h=250 m 对应焦散型NLS.现在在方程中加入惯性项,在水深较小的情形下(h=6.8 m ＜hc)参考Takizawa[1]实验中冰层的物理参数取惯性系数 α=0.069,考虑惯性后的NLS 仍然为焦聚型但非线性项的系数有所增大,孤立波由振幅无限小的周期波分岔而来,理论上此时的三阶截断模型会具有较高的精度.进一步地,数值计算波包型孤立波解以证明这一点.对于完全欧拉方程,采用保形映射将原不规则物理域映射为规则的矩形区域再结合牛顿迭代求解(关于此算法的详细介绍可见文献[22]).将由数值计算得到的二阶截断模型和三阶截断模型的分岔曲线与完全欧拉方程的分岔曲线进行对比(图1).在较大振幅范围内,三阶截断模型所对应的分岔曲线均与完全欧拉方程的分岔曲线吻合得很好,充分说明了三阶截断模型在该水深情况下具有较高的精度,是一个好的近似模型.而二阶截断模型的分岔曲线在分岔点附近与完全欧拉方程的分岔曲线相比仍然具有较为明显的差异,这表明二阶截断模型的精度较差.同时对相同波速下二阶截断模型、三阶截断模型以及完全欧拉方程的孤立波解进行比较后发现,三阶截断模型的解与完全欧拉方程的解基本吻合,两者之间的相对振幅差异约为10-2,而二阶截断模型的解则与完全欧拉的解相差较大.

图1 (a) 水深 h=6.8 m,惯性系数 α=0.069 时,孤立波解分岔图.实线:完全欧拉方程的分岔曲线,短划线:三阶截断模型的分岔曲线,点线:二阶截断模型的分岔曲线.(b) 波速 c=1.287 8 时,欧拉方程、三阶截断模型和二阶截断模型的孤立波解.实线:完全欧拉方程,短划线:三阶截断模型,点线:二阶截断模型Fig.1 (a) Bifurcation curves of wavepacket solitary waves for h=6.8 m and α=0.069 .Solid line:full Euler equations,dashed line:third-truncation model,dotted line:quadratic-truncation model.(b) The typical profiles for c=1.287 8 .Solid line:full Euler equations,dashed line:third-truncation model,dotted line:quadratic-truncation model

对水深较大的情况(h=250 m ＞hc),在参照Squire等[2]实验中的冰层参数后取惯性系数α=0.069.此时由于惯性项的影响,NLS 方程由焦散型变为焦聚型,因此三阶截断模型在水深较大的情况下仍然会具有较高的准确性.同样地,数值计算了该水深下二阶截断模型、三阶截断模型及完全欧拉方程的孤立波解及分岔曲线(图2).从图中可以看出,当孤立波振幅较小时三阶截断模型与完全欧拉方程的分岔曲线基本吻合,但在中等振幅时两者会出现一定的差异,这与水深h=6.8 m 时的情况略有不同.产生这种现象的原因在于三阶NLS 方程非线性项的系数:当h=250 m 时,虽然由于惯性项的加入使得方程由焦散型变为焦聚型,但此时的非线性项系数γ是一个非常接近于0 的正数,焦聚的特点表现得不明显,这导致三阶截断模型的准确性相比于水深h=6.8 m时略有下降.但在此水深情况下,二阶截断模型与完全欧拉方程之间的差距较水深h=6.8 m 时更大,完全无法反映出完全欧拉方程在对应波速处的波高,其精度远小于三阶截断模型,因此在水深较大的情况下使用二阶截断模型是不合适的.图2(b)进一步展示了波速c=1.287 8 时,二阶截断模型、三阶截断模型和完全欧拉方程的孤立波解波形比对比.因此综合两种水深下三阶截断模型和二阶截断模型与欧拉方程之间的对比,三阶截断模型显然具有较高的准确性,是一个更为合理的简化模型.

图2 (a) 水深 h=250 m,惯性 α=0.069 时,孤立波解分岔图.实线:完全欧拉方程的分岔曲线,短划线:三阶截断模型的分岔曲线,点线:二阶截断模型的分岔曲线.(b) 波速 c=1.287 8 时,欧拉方程、三阶截断模型和二阶截断模型的孤立波解.实线:完全欧拉方程,短划线:三阶截断模型,点线:二阶截断模型Fig.2 The bifurcation curves of wavepacket solitary waves for h=250 m and α=0.069 .Solid line:full Euler equations,dashed line:third-truncation model,dotted line:quadratic-truncation model.(b) The typical profiles for c=1.287 8 .Solid line:full Euler equations,dashed line:third-truncation model,dotted line:quadratic-truncation model

2.2 运动载荷下的冰层非线性响应

基于三阶截断模型式(5)和式(7),接下来将重点对匀速载荷作用下的冰层响应进行计算,并将数值结果分别与实验结果进行对比.实验中各物理参数的单位及具体取值如表1 所示.首先关注Takizawa[1]的实验结果.该实验是在日本北海道的L a k e Saroma 完成的,根据实验中所记录的参数可确定惯性项系数 α=0.069 2,而外加载荷p(x-Ut) 则通过数值拟合静止载荷下的冰层挠度来确定.在实验中Takizawa 观察并记录了载荷经过监测点的时间与监测点处冰层最大垂直位移出现时间的差异,这种滞后效应一般考虑是由于黏性效应造成的.在实验中黏性的来源是多样的,包括冰层本身的黏性,覆盖在冰层上的积雪的影响以及冰层与流体界面处的边界层效应等.Hosking等[16]及Wang等[17]引入双参数记忆函数建立黏弹性模型对运动载荷下的浮冰动力学响应进行研究.基于V＜cmin时滞后时间基本保持不变的特点,Hosking等[16]通过对多组数据进行试验并与实验记录进行对比确定了最佳的记忆函数参数A0和a0.Dinvay等[19]假设阻尼与垂直速度成正比,利用耗散项 -bηt拟合黏性效应,并将黏性系数视为待定参数,在数值计算中通过反复试验优化滞后时间来确定b.除此之外,这类确定待定系数的方法也在毛细重力波问题中得到了应用.Cho等[26]等在考虑黏性耗散对稳定及瞬时波形的影响时也将耗散项系数作为一个待定参数并通过将数值与实验结果进行拟合对比来确定具体取值.本文采用较为简化的线性耗散模型 -βηt描述实验中的黏性效应,β 的确定与上述参数的确定方法类似[16,19,26]:选取多组黏性系数并数值计算同一黏性系数下不同载荷速度对应的波形和滞后时间,将数值结果与实验记录对比直至两者能够较好吻合从而确定最佳黏性系数 β .经过反复实验确定黏性系数 β=0.4,此时数值结果与Takizawa[1]的实验结果吻合得最好.

表1 浅水[1]及深水[2]实验中各物理参数的取值及单位Table 1 Values and units of the physical parameters in the shallow-water[1] and deep-water[2] experiments

图3 展示了不同运动速度的载荷作用下的冰层动力学响应,总体而言数值结果较为准确反映了实验结果的主要特征.在准静态情形下,即U=2.2 m/s和U=4.2 m/s,数值结果与实验结果吻合得非常好.当载荷移动速度增大至U=5.5 m/s,此时速度接近于最小相速度,载荷两侧开始有波形成的趋势.载荷速度为 6 .2 m/s 时,载荷前后出现明显的波动,位于载荷前端的波长较短而后侧的波长较长,三阶截断模型的数值计算结果较好地捕捉到了这两类波.这说明三阶截断模型能够准确地反映出以最小相速度为分界点,不同速度下冰层响应的差异.进一步增大载荷移动速度至U=8.9 m/s,此时载荷速度大于浅水中的重力波速度,载荷后侧的重力波消失,出现“声影区”[6].图3(c)～图3(e)中的冰层最大挠度与实验结果基本一致,但载荷两侧的波形与实验记录存在一些差异.这可能是因为在数值计算中所采用的黏性项较为简单而实验中黏性来源相对复杂,随着载荷运动速度的增大,-βηt与实际黏性项产生的效果有所差异.除此之外,本文在建模时假设冰层厚度均匀,但在实验中冰层厚度不均.这些因素都可能会对载荷作用下产生的水弹性波波形造成影响.除了对比不同载荷速度下的冰层垂直位移,实验中出现的滞后效应在数值模拟中也得到了关注.实验中Takizawa[1]记录了不同速度下载荷经过测量点的时刻及该时刻下载荷在y方向上的位置,在图3 中以红点表示.在数值计算中选取任意一点为观测点后可得到该点处冰层垂直位移的时历曲线,同时记录移动载荷通过该点的时间及对应的垂向位置并在图3中以蓝点表示,从图中可以看出数值计算得到的滞后时间在各个速度下都与实验结果较为吻合.

图3 三阶截断模型的数值结果与文献[1]实验中所记录的冰层形变对比.其中虚线代表Takizawa 的实验结果,红点代表实验中的运动载荷经过冰层形变测量仪的垂直位置.实线表示三阶截断模型式(5)和式(7)的数值结果,蓝点代表数值计算中载荷经过监测点的垂向位置Fig.3 Comparisons of the numerical results of Eq.(3) the numerical results of cubic-truncation model and the experimental records of Ref.[1].The dashed lines represent the experimental data,and the red dots indicates the position of the load as it passed the deflectometer.The solid lines shows the numerical results of Eqs.(5) and (7),and the blue dots indicates the z-position of the load as it passes the point where the time series are obtained

在Squire等[2]的实验中,冰层厚度d=1.6 m,冰层弹性刚度D=1.8×109N·m,移动载荷重量为2100 kg.与Takizawa[1]不同的是,Squire等[2]在实验中将多个应变计安装在载荷运动轨道两侧,测量的是移动载荷作用下的冰层应变.因此需要对数值计算得到的冰层垂直位移 η 进一步处理以得到冰层线应变

其中,d为冰层厚度,κ*为无量纲前的冰层曲率[27].

考虑要与Squire 的实验进行对比,在基于三阶截断模型进行数值计算时,取 G0=(-∂xx)1/2以对应于深水情形.图4 展示了三阶截断模型的数值计算结果与Squire等[2]实验结果之间的对比,载荷速度分别为 4 .5 m/s,1 7.5 m/s,1 8.4 m/s,2 0.8 m/s,前两个载荷速度为亚临界,后两个速度为超临界速度.在亚临界速度下,数值结果与实验结果吻合得比较好,不仅能够拟合冰层应变的幅值,而且能够反映出随着载荷速度的增大应变曲线宽度逐渐变窄的趋势.当载荷以较为接近最小相速度的速度 1 7.5 m/s 运动时,载荷两侧逐渐升高并形成凹槽.在图4(c)和图4(d)中,数值预测的最大应变值略大于实验结果,尤其是应变曲线中载荷前端的波形.造成差异的原因除了以上在浅水实验中所提及的两个因素外,也可能是由于本文所采用的是二维模型,但在实验中所测得的应变是三维情形下冰层的应变,从而导致两者所得到的结果存在差别.但在超临界速度下,数值计算仍然较为有效地捕捉到了载荷前后所形成的弹性波和重力波.在数值计算中同样观察到了最大应变值点与载荷通过点之间的时间差异,但由于Squire等[2]只是指出了滞后效应的存在,没有明确记录差异时间,因此无法对滞后的时间进行对比.

图4 三阶截断模型的数值结果与文献[2]实验中所测量的冰层应变对比,图中应变值量级为 1 0-6 .虚线代表Squire 的实验结果,实线表示三阶截断模型式(5)和式(7)的数值结果Fig.4 The comparison of microstrain between the numerical results of cubic-truncation model and numerical records of Ref.[2].The microstrain is of 1 0-6 .Dashed lines:experimental records,solid lines:numerical approximation of the cubic-truncation model Eqs.(5) and (7)

3 结论

本文主要考虑非线性,惯性及黏性效应的影响,研究了移动载荷作用下的浮冰动力学响应.通过对Dirichilet-Neumman 算子进行展开将原始的完全非线性问题简化为仅与自由面上的变量相关的三阶截断模型.接着重点关注了此二维水弹性问题中的自由孤立波解以验证三阶截断模型的准确性.在不考虑黏性和外力载荷下,由多重尺度方法推导了三阶NLS 方程并基于此方程对波包型孤立波解的存在性和三阶截断模型的准确性进行了探究.惯性项的引入使得三阶NLS 方程在浅水和深水中均表现为焦聚型,因此理论上三阶截断模型在任意水深下都能够较好地近似完全欧拉方程.另一方面,对二阶截断模型、三阶截断模型和完全欧拉方程的分岔曲线和孤立波解进行数值求解并对比.数值结果表明任意水深下三阶截断模型在一定振幅范围内都能够与完全欧拉方程较好地吻合,是一个好的近似模型,而二阶截断模型则在任意水深下都与完全欧拉方程存在较大的差异,精度较差.进一步地,利用所得到的三阶截断模型对匀速运动载荷作用下的冰层响应进行了计算,并将数值结果与Takizawa[1](浅水情况)和Squire等[2](深水情况)实验分别进行了对比,结果表明此三阶截断模型能够较好地拟合移动载荷作用下冰层的垂直位移及应变情况.

附录

在推导考虑惯性项的三阶非线性薛定谔方程中,将解式(9)和式(10)代入原始欧拉方程中,得到各阶方程.这里给出推导的大致思路以及最终得到的三阶非线性薛定谔方程的系数表达式.首先解Laplace 方程,可以得到

这里n表示阶数,j表示模态,Qnj表示更低阶的项.于是可以得到

这里Pnj和Rnj仍然表示低阶项.在O(ε) 阶,由这两个方程的可解性条件,可以得到色散关系式.

在O(ε2) 阶,可以得到群速度cg

不难验证cg=ωk.同时也可以得到A22和 φ22的表达式

这里

在O(ε3) 阶,同样利用可解性条件,在复杂的计算之后,给出非线性薛定谔的系数表达式