水力学

畸形波生成及其影响因素的数值研究

  • 李梦雨 , 1, 2 ,
  • 吕超凡 , 1, 2 ,
  • 卢金友 1, 2 ,
  • 栾华龙 1, 2 ,
  • 朱勇辉 1, 2 ,
  • 朱家熙 1, 2, 3 ,
  • 葛建忠 4 ,
  • Makiko Iguchi 5
展开
  • 1 长江科学院 河流研究所, 武汉 430010
  • 2 长江科学院 水利部长江中下游河湖治理与防洪重点实验室, 武汉 430010
  • 3 河海大学 港口海岸与近海工程学院, 南京 210098
  • 4 华东师范大学 河口海岸全国重点实验室, 上海 201100
  • 5 日本水电综合技术研究所, 日本 大阪 5300000
吕超凡(1998-),男,江西九江人,博士,主要研究方向为河口海岸动力学。E-mail:

李梦雨(1992-),女,河南邓州人,高级工程师,博士,主要研究方向为河口海岸防灾减灾。E-mail:

Copy editor: 任坤杰

收稿日期: 2024-09-30

  修回日期: 2025-03-18

  网络出版日期: 2025-12-11

基金资助

国家重点研发计划项目(2022YFE0117500)

国家自然科学基金项目(52301310)

国家自然科学基金项目(42376166)

中国博士后科学基金项目(2023M740369)

中国博士后科学基金项目(2024T170761)

中央级公益性科研院所基本科研业务费专项资金资助项目(CKSF2023321/HL)

中央级公益性科研院所基本科研业务费专项资金资助项目(CKSF2024326/HL)

中央级公益性科研院所基本科研业务费专项资金资助项目(CKSF20231009/HL)

Numerical Study on Freak Wave Generation and Its Influencing Factors

  • LI Meng-yu , 1, 2 ,
  • LÜ Chao-fan , 1, 2 ,
  • LU Jin-you 1, 2 ,
  • LUAN Hua-long 1, 2 ,
  • ZHU Yong-hui 1, 2 ,
  • ZHU Jia-xi 1, 2, 3 ,
  • GE Jian-zhong 4 ,
  • Makiko Iguchi 5
Expand
  • 1 River Research Department, Changjiang River Scientific Research Institute, Wuhan 430010, China
  • 2 Key Laboratory of River-Lake Regulation and Flood Control in the Middle and Lower Reaches of Yangtze River of Ministry of Water Resources, Changjiang River Scientific Research Institute, Wuhan 430010, China
  • 3 College of Harbour, Coastal and Offshore Engineering,Hohai University, Nanjing 210098, China
  • 4 State Key Laboratory of Estuarine and Coastal Research, East China Normal University, Shanghai 201100, China
  • 5 Hydro Technology Institute Co., Ltd., Osaka 5300000, Japan

Received date: 2024-09-30

  Revised date: 2025-03-18

  Online published: 2025-12-11

摘要

畸形波生成及演化的影响因素众多,基于自主研发的黏性流数值波浪水槽,开展畸形波生成及其影响因素的数值研究。首先利用试验结果验证数值波浪水槽生成畸形波的准确性,然后模拟分析畸形波传播过程中的波面变化,采用谐波分离方法分析波群非线性对畸形波传播过程中的波面变形、聚焦位置和时间及内部结构变化的影响规律;之后开展数值试验,分析谱型、波群的组成波个数、频谱宽度、谱峰频率及水深等关键因素对畸形波生成的影响规律,深入探究各个因素的敏感性。结果表明:采用不同频谱类型产生的畸形波表现出不同的能量分布和聚焦时空变化。同时,组成波个数影响波群的聚焦重现周期。在有限水深情况下,畸形波的聚焦幅值并非随着谱宽的变窄而增大,同时受水深、谱峰频率等因素的影响较大。该研究有助于理解畸形波的生成机制,对实验室生成畸形波具有参考意义。

本文引用格式

李梦雨 , 吕超凡 , 卢金友 , 栾华龙 , 朱勇辉 , 朱家熙 , 葛建忠 , Makiko Iguchi . 畸形波生成及其影响因素的数值研究[J]. 长江科学院院报, 2025 , 42(12) : 75 -85 . DOI: 10.11988/ckyyb.20241029

Abstract

[Objective] Freak wave is a marine disaster characterized by extremely large wave height, strong nonlinearity, and high destructiveness. The results of wave superposition method for simulating freak waves are influenced by multiple parameters, and the sensitivity and interaction mechanisms of these factors require systematic investigation. [Methods] Based on a self-developed viscous-flow numerical wave tank, we conducted a numerical simulation on the generation of freak waves and their influencing factors. First, the reliability of the numerical model was verified against physical experimental data. Subsequently, the harmonic separation method was employed to examine the influence of wave group nonlinearity on wave surface deformation, focusing characteristics, and frequency spectrum structure. Through numerical experiments, the effects of key parameters—including spectral type, number of constituent waves, spectral bandwidth, spectral peak frequency, and water depth—were investigated. [Results] 1) During the generation of a freak wave, wave-wave nonlinear interactions caused energy to transfer from the primary frequency to both high and low frequencies, resulting in significant spectral broadening. Low-frequency free pseudo-harmonics propagated faster, leading to an actual wave height slightly larger than the theoretical value. High-frequency bound harmonics formed a tail wave, which had a minor influence on the shape of the main peak. 2) The spectral type significantly influenced the wave profile characteristics: the JONSWAP and P-M spectra, with concentrated energy, tended to generate freak waves with steep crests. The CWA spectrum produced gentle wave profiles; the CWS spectrum yielded the smallest focused amplitude. 3) The number of constituent waves affected the focusing recurrence period. An insufficient number could generate secondary focused waves. It was recommended to use 29 constituent waves to balance computational accuracy and efficiency. 4) Under finite water depth conditions, the focused amplitude reached its maximum when the spectral bandwidth was 0.7 Hz, indicating that the amplitude was co-modulated by the spectral bandwidth, water depth, and spectral peak frequency. 5) An increase in the spectral peak frequency enhanced nonlinearity, resulting in wave profile steepening. However, an excessively high frequency led to wave breaking, thereby reducing the amplitude. 6) Water depth influenced the wave profile by altering the dispersion characteristics. A greater water depth resulted in faster wave speed and a higher amplitude, whereas an excessively small water depth readily induced wave breaking. [Conclusion] The main innovations of this research include: establishing a high-precision viscous-flow numerical model capable of accurately simulating the evolution of nonlinear waves including breaking effects; employing the harmonic separation method to reveal the influence mechanism of wave group nonlinearity on wave surface structure and energy distribution; and clarifying the coupling effects of various factors under finite water depth conditions through multi-parameter sensitivity experiments. The findings of this study deepen the understanding of freak wave generation mechanisms, provide an important theoretical basis and parameter selection guidance for laboratory simulation of freak waves.

0 引言

畸形波是一种具有极大波高、强非线性的灾害性波浪[1],不仅发生在深海和浅海区域,也常发生在近岸海域[2-3]。据记载,1969—1994年,世界海域内共有22艘巨轮因遭遇畸形波而失事,造成约525人伤亡[4]。近年来随着气候变化,畸形波发生的频率也有所增加,据统计,2005—2021年,世界范围内就出现了429个畸形波事件(其中224个发生在近岸海域),共造成658人伤亡[5]。因此,深入了解畸形波的内部特征、准确地再现其时空演化过程,对于预报及预防畸形波危害、保障人类生命财产安全具有重要意义。
由于畸形波发生的环境多样,出现的时间地点随机,不同学者对其发生机理的诠释也有不同。畸形波生成可归结为一种或多种能量汇聚效应的叠加,主要受外部环境(如风场、地形)和内部动力(如波列不稳定性、波-流及波-波非线性作用)的影响。根据不同的影响因素,畸形波的生成机理可概括为空间聚焦、风场能量输入、时空聚焦、波浪调制不稳定性、波-波非线性相互作用以及非线性波群相互作用等[5-7]
线性时空聚焦理论是实验室常用的再现畸形波的方式之一[8]。根据线性波浪理论将不同振幅和周期的波浪进行叠加,并调制各个组成波的初相位,使其在目标时刻和位置达到最大峰值, 从而产生目标极端大波;由于实施方便,该方法是目前实验室产生畸形波最有效的方法之一,一些学者基于该方法开展了大量的研究工作。Brandini等[9]通过组成波波能汇聚的方法在数值水池中模拟生成了畸形波,并分析了畸形波的外部和运动学特征;黄国兴[10]采用线性时空聚焦方法,基于二阶的高阶谱方法,通过人工干预组成波的随机初相位得到了包含极端大波的波面时间序列;Zhao等[11]基于Navier-Stokes(N-S)方程,采用线性时空聚焦方法在二维水槽中模拟了聚焦波的产生及变形;Ning等[12]基于高阶边界元法引入“新波”(New Wave) 模型模拟了聚焦波的生成。裴玉国[13]基于双波列叠加的组合模型实现了畸形波在实验室的定时、定点生成;赵西增[14]基于高阶谱方法,采用波列叠加以及相位调制等模型模拟了畸形波的生成;刘赞强[15]基于Longuet-Higgins模型提出了一种模拟畸形波的相位调制方法,数值再现了含有极端大波的实测波列,并分析了极端波浪的特征参数;崔成等[16]采用三波列叠加模型在实验室内实现了畸形波的定点、定时生成,分析了畸形波生成及演化的整个过程和基本特征;扈喆[17]对传统波列叠加模型做出改进,提出了基于概率的波列叠加模型并有效模拟了日本Yura渔场的实测畸形波。Yu等[18]基于Boussinesq方程建立了模拟畸形波的相位聚焦模型,并对畸形波的演化过程和非线性特征参数进行了分析;张亚荣[7]利用2个单聚焦波群叠加模拟了双波群聚焦产生畸形波的过程,研究发现双波群间频率宽度越小,畸形波越容易产生;王磊[19]选取2个能量集中的聚焦波群进行能量叠加以模拟畸形波的生成,分析了双波群聚焦过程中的演化特性及畸形波生成过程中波浪非线性作用的影响,解释了2个波群间的相互作用引起的三阶非线性是双波群相比于单波群线性叠加出现相位滞后和振幅谱变化的原因;Wang等[20] 建立了单峰和双峰波列之间的最大峰度定量关系,并推导得出了双峰波峰度与单峰海况峰度关系的经验公式;崔成等[21]采用色散和方向聚焦方法在物理水池中模拟了三维短峰畸形波的生成和演化过程;刘震等[22]基于波列叠加模型在大型拖曳水池模拟了获得的包含畸形波的目标波列,分析了畸形波波面的时空演化过程中最大波高和有效波高及其比值、波幅谱、峰度和偏度的沿程变化等。
综上,通过调制波列相位使波群的能量叠加能够在实验室水槽中指定的位置和时刻产生目标畸形波。对于基于波列叠加理论产生的畸形波,其生成往往通过海浪谱来实现。在实验室模拟畸形波生成的过程,模拟结果受多种因素影响,目前这些因素对畸形波生成结果的影响和敏感性尚不明晰,同时波群内部结构组成还需进一步揭示。另外考虑畸形波生成及传播过程涉及强非线性自由面问题,采用势流理论的数学模型往往在模拟过程中存在一定的局限性,而基于黏性流理论的数学模型,由于方程中考虑了非线性项和黏性项,使其能够更准确地模拟波浪破碎、翻卷、波-物相互作用等强非线性自由面流动。因此,本文通过自主编程建立黏性流数值波浪水槽,基于波列聚焦理论在实验室生成畸形波,深入探究畸形波生成过程中波群内部的结构组成,并揭示关键因素对畸形波生成及演化的影响规律,为实验室生成畸形波提供参考。

1 数值波浪水槽

1.1 控制方程

基于直角笛卡尔坐标系,不可压缩黏性流体流动的质量守恒方程和动量守恒方程为:
$\frac{\partial {u}_{i}}{\partial {x}_{i}}=0 ,$
$\begin{array}{l}\\ \frac{\partial {u}_{i}}{\partial t}+{u}_{j}\frac{\partial {u}_{i}}{\partial {x}_{j}}=-\frac{1}{\rho }\frac{\partial p}{\partial {x}_{i}}+\frac{1}{\rho }\frac{\partial }{\partial {x}_{j}}\left(\mu \frac{\partial {u}_{i}}{\partial {x}_{j}}\right)+\\ {g}_{i}+{{s}^{\mathrm{m}}}_{i}+{A}_{\mathrm{b}} 。\end{array}$
式中:xi(i=1, 2, 3)为笛卡尔直角坐标系中的坐标轴;ui(i=1,2,3)为速度分量;t为时间;ρ为流体密度;p为压强;gi为体积力,如重力等;μ为动力黏性系数; ${s}_{i}^{\mathrm{m}}$为造波动量源项;Ab为阻尼消波源项。
畸形波的生成涉及液体、空气和固体的多相流动,采用体积函数 ${\varphi }_{m}$(m =1、 2、 3分别表示液体、空气和固体) 捕获交界面,即
$\frac{\partial {\varphi }_{m}}{\partial t}+{u}_{i}\frac{\partial {\varphi }_{m}}{\partial {x}_{i}}=0 。$
式中φm为体积函数,取值0~1。为了克服界面的不稳定性,模型中假定液体和固体的密度和黏度相同。每个网格内固体φ3的体积函数由拉格朗日方法确定,空气φ2的体积函数φ2= 1-φ1-φ3。每个网格的物理性质λ(如密度ρ和动力黏性系数μ) 可由下式计算
$\mathrm{\lambda }=\stackrel{3}{\sum _{m=1}}{\varphi }_{m}{\lambda }_{m} 。$

1.2 畸形波生成方法

根据时空聚焦模型[23],通过调制波群中各个组成波的初始相位,在物理和数值水槽特定的位置和时刻实现组成波能量的汇聚,产生特定波高的畸形波。数学模型基于波列聚焦理论、采用基于Boussines方程推导得出的动量源项造波方法来生成畸形波。根据Choi等[24]的研究,对于二维线性单频规则波的动量源项sm可简化为
${s}^{\mathrm{m}}=-g\left(2\beta x\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(-\beta {x}^{2}\right)\frac{D}{\omega }\mathrm{s}\mathrm{i}\mathrm{n}\left(-\omega t\right) 。$
式中:x为坐标轴;t为时间;ω为组成波的频率;βD为与造波源区宽度相关的因素。
根据波列聚焦理论,在数值水槽中通过线性叠加单频波生成畸形波的动量源项表示为
$\begin{array}{l}{s}^{\mathrm{m}}=-g\left(2\beta x\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(-\beta {x}^{2}\right)·\\ \stackrel{N}{\sum _{i=1}}\frac{{D}_{i}}{{\omega }_{i}}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{ }\left[{k}_{i}\left(x-{x}_{\mathrm{f}}\right)+{\omega }_{i}\left(t-{t}_{\mathrm{f}}\right)\right] 。\end{array}$
式中:N为波群中组成波的个数;xf为聚焦位置;tf为聚焦时间;ki为波群中第i个组成波的波数。

1.3 阻尼消波方法

本文采用动量源项阻尼消波法,即动量方程中添加耗散源项,实现阻尼消波。阻尼吸收源项Ab采用Liu等[25] 的方法,即
${A}_{\mathrm{b}}={c}_{\alpha }\frac{\mathrm{e}\mathrm{x}\mathrm{p}\left[{\left(\frac{\left|x-{x}_{\mathrm{s}\mathrm{t}}\right|}{{x}_{\mathrm{a}\mathrm{b}}}\right)}^{{n}_{\mathrm{c}}}\right]-1}{\mathrm{e}\mathrm{x}\mathrm{p}\left(1\right)-1} 。$
式中:Ab为吸收系数;xstxab分别为吸收区域的起始位置和长度;dcnc为通过数值试验确定的经验阻尼系数,本研究采用nc=10和dc=min {200; 0.9/dt},dt为计算时间步长。

1.4 数值方法及求解

基于交错网格笛卡尔坐标系统,控制方程采用有限差分法对非均匀矩形网格变量进行离散求解。首先,采用具有高精度且稳定的半拉格朗日方法—紧致差值曲线(Constrained Interpolation Profile,CIP)方法求解对流项;采用中心差分法求解扩散项、造波动量源项和阻尼项;通过投影法求解压力项,并更新整个流场。自由表面的获取采用具有斜率加权的双曲线切线界面捕捉(Thinc with Slope Weighting,THINC/SW)方法。固体边界采用浸入边界法(Immersed Boundary Method,IBM)求解。数值水槽的相关内容可参考文献[26]与文献[27]。

2 畸形波的生成

2.1 水槽设置及工况

物理水槽试验在浙江大学舟山校区近海馆大断面水槽内进行[27],水槽长75 m、宽1.8 m、高2.0 m,水槽内沿其中线布置编号为G1—G6的6 个测波仪监测自由面变化。试验中恒定水深h0为0.5 m,选用连续波幅(Continuous Wave Amplitude,CWA)谱,共取29个组成波,聚焦位置均设置在G3 处。基础试验中考虑6种不同波高和谱宽的畸形波进行分析,试验工况如表1所示。有关物理试验的相关内容参考文献[24]。
表1 基础试验工况

Table 1 Basic experimental conditions

工况
编号
频谱范围
f/Hz
频谱宽
df/
Hz
中心频
fc/
Hz
聚焦振
Af/m
波陡 色散
参数
1 0.53~1.13 0.6 0.83 0.03 0.092 [0.83,2.60]
2 0.53~1.13 0.6 0.83 0.06 0.184 [0.83,2.60]
3 0.53~1.13 0.6 0.83 0.09 0.276 [0.83,2.60]
4 0.38~1.28 0.9 0.83 0.03 0.092 [0.57,3.31]
5 0.38~1.28 0.9 0.83 0.06 0.184 [0.57,3.31]
6 0.38~1.28 0.9 0.83 0.09 0.276 [0.57,3.31]

注:波陡为kcAf,其中kc为中心波的波数,即为单位长度内波的周期数(1/m);色散参数为kh,k为波数(1/m),h为水深(m)。

基于上述物理水槽参数建立数值波浪水槽,计算域见图1。考虑时间成本,使用长度28.0 m的数值波浪水槽模拟长度75.0 m的物理波浪水槽,水槽以造波区的中心为原点(x=0),共设置6个波高仪(编号为G1—G6),其位置依次为x=6.00、10.79、11.54、12.29、13.08、15.18 m。为保证物理水槽和数值水槽中波浪演化规律一致,以物理水槽第一个浪高仪的信息确定源项造波位置处的输入信号。
图1 数值波浪水槽示意

Fig.1 Schematic diagram of numerical wave tank

2.2 网格选取及收敛性验证

采用3套不同分辨率的非均匀网格进行计算收敛性验证,见表2图2(a)为工况1中不同网格在G3浪高仪处的波群在聚焦位置的自由表面历时变化,横坐标为相对聚焦时刻t - Tf(t为波面历时,Tf为理论的聚焦时间),纵坐标η/A为相对波面高程(η为波面高程,A为理论的聚焦幅值)。由图2(a)可知,细网格和中网格的模拟结果基本相同,但对于粗网格,模拟结果相比理论解和试验数据存在一定的偏差。不同网格系统下流体质量守恒相对误差E
$\mathrm{E}=\frac{\left|\sum _{i,j}\left({\varphi }_{ij}^{n}\mathrm{d}x\mathrm{d}y\right)-\sum _{i,j}\left({\varphi }_{ij}^{0}\mathrm{d}x\mathrm{d}y\right)\right|}{\sum \left({\varphi }_{ij}^{0}\mathrm{d}x\mathrm{d}y\right)} 。$
式中: ${\varphi }_{ij}$为(i,j)网格内的体积函数;dx、dy分别为xy方向的网格长度; $\sum _{i,j}\left({\varphi }_{ij}^{0}\mathrm{d}x\mathrm{d}y\right)$为初始条件下计算域的总体积; $\sum _{i,j}\left({\varphi }_{ij}^{n}\mathrm{d}x\mathrm{d}y\right)$ 为计算n时刻后计算域的总体积。计算完成后计算域的质量守恒误差见表2。由表2可知,网格越细,质量守恒误差值越小,且质量守恒误差都在10-5数量级,表明数值模型的守恒性较好。总体上,采用中等网格即可实现计算收敛。利用中网格系统对不同时间步长(ΔtLarge=T/200、ΔtMiddle=T/500和ΔtSmall=T/1 000,其中T为计算时长)进行时间步长收敛性测试,见图2(b)。结果表明采用中间时间步长(ΔtMiddle=T/500)即可实现计算收敛。
表2 三套网格相关参数

Table 2 Parameter sets of coarse, middle, and fine grids

网格
类型
x
网格数
Nx
y
网格数
Ny
网格数
Nx×Ny
x向最
小网格
dx/m
y向最
小网格
dy/m
相对误
E
粗网格 314 89 27 946 0.07 0.007 7×10-5
中网格 424 150 63 600 0.05 0.005 4×10-5
细网格 614 190 116 660 0.03 0.003 3×10-5
图2 工况1不同网格系统和时间步长在G3浪高仪处的波面历时变化

Fig.2 Time-histories of surface elevation at wave gauge G3 for different grids and time steps in Case 1

2.3 畸形波波面演变特征分析

根据以往相关研究,选取非线性最强的工况3[26]分析畸形波传播过程中波面历时变化。由图3可知,在距离聚焦位置较远的浪高仪G1处,波浪以连续波群的形式出现;在接近聚焦位置浪高仪G2处,出现了一个极端的主波峰,主波峰两侧的波谷和二级波峰并不对称;在聚焦位置浪高仪G3 处,波群的能量瞬间聚焦,主波峰变得很高且两侧的波谷和二级波峰较对称;波群大波汇聚持续的时间很短,在聚焦位置之后浪高仪G4 处,主波峰峰值变小,主波峰两侧的波谷和二级波峰变得不再对称;在远离聚焦位置后,波群能量不再汇聚,波列逐渐分散且有小尾波出现。整体上计算结果与试验结果吻合良好,验证了数值波浪水槽模拟产生畸形波的准确性。
图3 工况 3中畸形波生成过程的波面历时变化

Fig.3 Time-histories of surface elevation during freak wave generation in Case 3

2.4 波群非线性对波群自由面变形及频谱结构的影响

波群非线性对其聚焦位置、时间及波面有很大影响[26]。为探究波群非线性对畸形波生成的影响,基于工况3模拟波峰及波谷聚焦的畸形波,并引入谐波分离方法[28],探究畸形波的内部组成结构,分析不同组成成分对波群自由面变形及频谱结构的影响。
图4(a)为工况3中畸形波在水槽传播过程中的时空演变,图4(b)为谐波分离后波群线性项的时空演化,由于谐波分离,图4(c)为分离谐波后波群非线性项的时空演化。由图4(a)可知,畸形波出现之前慢速传播的短波位于长波的前面,随着相位的发展,各个组成波在聚焦位置处能量汇聚,波群中出现一个大波。波群聚焦之后,由于频散关系,波群中各个组成波以各自的速度向前传播,波群变得分散。在图4(b)中,由于谐波分离,大波位置处的能量汇聚程度要弱于图4(a)。在图4(c)中波群非线性项分为a, b, c三部分,其中b传播速度与线性项相同,其为高阶约束波;a传播速度快于线性波,其为线性造波导致的低频误差波;c传播速度慢于线性波,其为线性造波导致的高频误差波。
图4 工况3中畸形波传播过程中的时空演变

Fig.4 Spatiotemporal evolution of freak wave propagation in Case 3

图5(a)为聚焦位置和时间处的自由面结果,分别为线性理论解、计算结果、谐波分离后的线性项和非线性项的结果,由于非线性的存在,聚焦时刻滞后。由图5(a)可知,谐波分离后的线性项非常接近线性理论解,峰值和相位产生轻微偏差的原因是本文谐波分离方法只考虑二阶项,实际中谐波分离后的线性项还包含着高阶奇频项(例如三阶项、五阶项等),而它们占主频波的比例很小。
图5 工况3中畸形波的内部结构组成

Fig.5 Internal structural composition of freak wave in Case 3

图5(b)为畸形波波面内部各个组成部分的频谱分布。由图5(b)可知,畸形波在传播过程中,由于波-波非线性相互作用,能量由主频向高频和低频转移,波谱变宽。经过谐波分离后的线性项,能量主要分布在主频范围内,另外少部分能量分布在低频和高频处,次峰主频率约为0.2 Hz和2.5 Hz,这导致图5(a)中谐波分离出的线性项与线性理论解出现偏差。经过谐波分离后的非线性项,能量主要分布在低频和二倍主频范围内。
采用高通低通滤波方法分别得到高频和低频所对应的时域信息。图6(a)为低频谐波下的自由面变化,图6(b)为自由面对应的振幅谱信息。由图4(c)可知,低频自由伪谐波(次谐波)传播较快,其以波峰形式的孤立子出现,而低频约束谐波传播较慢,其在聚焦时刻以波谷形式的孤立子出现(见图6(a)),其对应的频谱信息见图6(b)。因为自由伪谐波的振幅为正,这也是聚焦位置处的波高比实际波高稍大的原因,但是本研究中自由伪谐波的波峰为2.8×10-3 m,其占理论波高的比例约为2.4%,计算中可忽略不计其误差影响。
图6 工况3中畸形波内部低频谐波分布

Fig.6 Distribution of low-frequency harmonics inside freak wave in Case 3

图7(a)为高频谐波下的自由面变化,图7(b)为自由面对应的振幅谱信息。由图7(a)可知,波面历时变化先后出现了2个波列,依据由约束波和自由伪谐波的性质[29],在聚焦时刻(t=20 s左右)出现的为高频约束波,在聚焦时刻大概10 s后(t=30 s左右)以波群形式出现的为高频自由伪谐波(超谐波)。由于高频自由伪谐波传播速度较慢,主要以波群后面的小尾波形式出现,振幅较小(图7(b)),其并不影响波群聚焦位置处的波面形态。
图7 工况3中畸形波内部高频谐波分布

Fig.7 Distribution of high-frequency harmonics inside freak wave in Case 3

3 畸形波生成的影响因素分析

畸形波的生成与谱型、组成波的个数、频谱范围、聚焦幅值、谱峰频率、波陡及水深有关,开展数值试验分析上述因素对极端波浪生成及发展的影响规律,探究各个因素的敏感性。

3.1 不同频谱类型的影响

波浪频谱形状与其生成机理有关,主要决定于当地的风速、风距、风时、地形及水深等因素。本节选取P-M(Pierson-Moskowitz)谱,JONSWAP谱,等波幅(CWA)谱,等波陡(CWS)谱等工程实际中常用的频谱,分析不同谱型对畸形波生成的影响,选取的算例水深h=0.5 m,考虑波浪波峰聚焦情况,谱宽df=0.6 Hz,频谱范围为0.53~1.13 Hz,其覆盖了大部分波能,谱峰频率fp=0.83 Hz,波群组成波的个数N=29,聚焦振幅Af=0.09 m,输入聚焦位置xf=10.2 m,输入聚焦时间Tf=20 s。不同谱型的波能分布如图8所示,横坐标为频率f,纵坐标为能量谱S
图8 不同谱型的波能分布

Fig.8 Distribution of four selected spectral types

图9展示了不同谱型下畸形波的波面时空变化。由图9可知,采用不同谱型的波群具有不同的传播速度。波群中传播较慢的为短波,传播快的为长波。初始时刻,慢速传播的短波位于长波的前面,随着相位的发展,各个组成波之间相互作用,在聚焦时间和位置处,能量汇聚,波群中出现一个大波。在聚焦位置后,由于频散关系,波群中各个组成波之间的相互作用减弱,波群不再那么紧凑。另外,从图中可以观察到长时间传播过程中短波的衰减较严重。对于CWA谱,其波能分布较均匀,波群中低频和高频组成波所占比例较大,在传播过程中波群中组成波分散较为明显(图9(a));对于CWS 谱,由于低频组成波所占比例较小,因此波群中大部分组成波传播速度较快(图9(b));对于P-M谱和JONSWAP 谱,其能量分布较为集中,波群传播过程中尤其在聚焦位置附近整体看起来较为紧凑(图9(c)图9(d))。
图9 不同谱型下畸形波的波面时空变化

Fig.9 Spatiotemporal variation of freak wave under different frequency spectra

图10为不同谱型下畸形波的波面历时变化。由图10(a)可知,对比聚焦振幅,P-M 谱下波群聚焦振幅最大,其次为JONSWAP谱和CWA谱,而CWS谱下聚焦振幅最小;对比聚焦时刻,CWS谱聚焦时刻距离输入时刻最近,其次为CWA谱、JONSWAP谱和P-M谱。由图10(b)可知,不同谱型下畸形波聚焦波面处波形差异较大。CWS谱下畸形波的自由面结果较为胖宽,波谷较平坦,最大波峰两侧的次峰较低;而CWA谱下畸形波最大波峰两侧的次峰有一些小的波动;P-M 谱和JONSWAP谱下畸形波的聚焦波面除了聚焦幅值存在差异外,其他整体较为相似,而相较CWA谱和CWS谱,其波谷更深,波峰更尖,并且最大波峰两侧的次峰更高。
图10 不同谱型下畸形波的波面历时变化

Fig.10 Time-histories of surface elevation for freak wave under different spectral types

不同谱型下畸形波波群的聚焦振幅、相对聚焦时刻和位置见表3。CWA谱型波能分布比较均匀(图8),波群内波-波之间相互作用不及P-M和JONSWAP谱强烈,波浪的非线性特征较弱;而对于CWS谱型能量分布,低频长波能量占比大,而高频强非线性短波能量占比较小,因此其聚焦时刻与理论聚焦时刻相差较小,真实聚焦振幅相较P-M和JONSWAP谱小;而JONSWAP谱能量集中分布在谱峰频率附近,波-波之间相互作用较强,其真实聚焦振幅与输入值的偏差较大;对于P-M谱能量分布,高频强非线性短波能量占比较大,因此其与理论输入值偏差最大。由上可知不同谱型代表特定的问题,模拟过程中需根据当地海况及解决的实际工程问题选用合适的波浪频谱。
表3 不同谱型计算结果

Table 3 Computational results for different spectral types

谱型 相对聚焦波幅 相对聚焦位置/m 相对聚焦时间/s
CWA 1.177 1.00 0.40
CWS 1.103 0.80 0.28
P-M 1.252 1.50 0.76
JONSWAP 1.202 1.48 0.73

3.2 组成波个数的影响

天然海浪由无数个组成波叠加而成,如果模拟选取的组成波个数N足够大,则模拟的波面特性接近于真实海况。但在实际模拟中由于操作及效率的限制,常采用有限组成波来模拟目标海浪。为了探究组成波的个数对畸形波波面特性的影响规律,本小节选取实验室常用的JONSWAP谱,静水深h=0.5 m,谱宽df=0.6 Hz,聚焦波幅Af=0.09 m,谱峰周期Tp=1.2 s,输入聚焦位置Xf=10.2 m,输入聚焦时间Tf=20 s。参考文献[30],选取29个组成波生成聚焦波,另外选取5、15、51这3种组成波个数,开展数值研究,分析组成波个数N分别为 5、15、29及51(频率间隔Δf分别为0.2、0.04、0.02和0.01 Hz)对畸形波波面特性的影响。
图11为JONSWAP谱型波列中不同组成波个数下畸形波波面时空变化。由图11可知,组成波个数影响聚焦波的重现周期,组成波个数越多,聚焦波重现周期越长,组成波个数为5、15、29和51的聚焦波重现周期分别为5、25、50及100 s(周期为1/Δf)。
图11 JONSWAP谱型不同组成波个数下畸形波波面时空变化

Fig.11 Spatiotemporal variation of freak wave surface with different numbers of constituent waves under JONSWAP spectrum

图12为JONSWAP谱型聚焦位置处不同组成波个数的自由面变化。组成波的个数越多,聚焦波的重现周期越长;组成波个数N=5的波面与N=15、29和51相比出现了相位和幅值的偏差;说明足够大的组成波个数对畸形波特性无影响,经综合考虑,本数值模型模拟的波群组成波个数N取29。特别说明,波群中的组成波个数会影响畸形波的聚焦重现周期,组成波个数过少可能导致次聚焦波的出现,因此,在实际模拟中需要考虑模拟时间、聚焦时间及位置等因素以确定最佳组成波个数,避免选取组成波个数过少出现次聚焦波影响计算结果,组成波个数过多又会影响计算效率。
图12 JONSWAP谱型不同组成波个数下畸形波的自由面历时变化

Fig.12 Time-histories of free surface elevation for freak wave with different numbers of constituent waves under JONSWAP spectrum

3.3 频谱宽度的影响

Baldock等[30]表明在深水情况下,输入的谱宽越窄,聚焦位置处的波幅值越大。为探究在有限水深情况下是否能得到同样的规律,重点分析了频谱宽度对波浪聚焦结果的影响规律。选取实验室常用的JONSWAP谱,组成波个数N=29,谱宽df=0.1~0.9 Hz,谱峰频率fp=0.83 Hz,各组成波之间的频率间隔Δf均为0.01,聚焦振幅Af=0.09 m进行数值模拟。考虑水深h=0.5 m,输入聚焦位置xf=10.2 m,输入聚焦时间Tf=20 s。
图13为不同谱宽下JONSWAP谱型,横坐标为频率f,纵坐标为能量谱S。由图13可知,随着谱宽变窄,波群组成波能量更集中,单个组成波的振幅更大,组成波之间差异较小。
图13 不同谱宽下JONSWAP谱型

Fig.13 JONSWAP spectra with different spectral bandwidths

图14为聚焦位置处不同频谱宽度下畸形波的自由面变化。由图14可知,谱宽越宽,聚焦时刻距离输入时刻越近,且聚焦位置畸形波的波面越胖宽,波谷位置较高且平坦,次峰位置较低。
图14 聚焦位置处不同谱宽下畸形波的自由面历时变化

Fig.14 Time-histories of free surface elevation for freak wave with different spectral bandwidths at focal position

图15为不同谱宽下畸形波在聚焦位置处的振幅,随着谱宽增加,聚焦位置处振幅呈先增大后减小趋势,当df=0.7 Hz时振幅达到最大值;当df=0.1 Hz时,由于波群形态趋于规则波,聚焦振幅不符合变化规律。这表明对于有限水深情况,并非谱宽越窄聚焦振幅越大,除了与频谱宽度有关,聚焦振幅还与水深、谱峰频率等因素相关。
图15 不同谱宽下畸形波在聚焦位置处的振幅

Fig.15 Amplitude of freak wave at focal position under different spectral bandwidths

3.4 谱峰频率的影响

谱峰频率fp为谱型的一个重要参数,为探究谱峰频率fp对畸形波产生过程中自由面变化的影响,选取实验室常用的JONSWAP谱,谱宽df=0.6 Hz,频谱范围为0.53~1.13 Hz,考虑静水深h=0.5 m,组成波个数N=29,聚焦振幅Af=0.09 m,聚焦位置xf=10.2 m,聚焦时间Tf=20 s。选取谱峰频率fp=0.63、0.73、0.83、0.93及1.03 Hz(谱峰波陡kpA分别为0.19、0.23、0.27、0.33及0.40)进行数值模拟。图16为不同谱峰频率下JONSWAP谱型。
图16 不同谱峰频率下JONSWAP谱型

Fig.16 JONSWAP spectra with different peak frequencies

图17为不同谱峰频率下畸形波聚焦位置处自由面历时变化。由图17(a)可知,fp=0.93 Hz下波浪的聚焦振幅最大,其次为fp=0.83、0.73 和0.63 Hz,而fp=1.03 Hz下波浪的聚焦振幅最小。对比相对聚焦时刻,fp=0.63 Hz下波浪聚焦时刻距离输入时刻最近,其次分别为fp=0.73、0.83、1.03和0.93 Hz。由图17(b) 可知,不同谱峰频率下畸形波波面聚焦的结果差异较大。fp=0.63 Hz下聚焦位置处的自由面结果较为胖宽,波谷位置较高且较为平坦,最大波峰两侧的次峰出现的位置较远并且位置较低。随着谱峰频率的增加,聚焦位置处的自由面变得高窄,最大波峰两侧的次峰出现的位置越早,波形变得更陡。而对于fp=1.03 Hz的情况,其自由面出现了次波并且峰值较小。由图16可知,fp=0.63 Hz的波群中低频长波能量分布较多,而高频强非线性短波所占能量比例较小,波浪的非线性特征没有那么强烈,所以其聚焦时刻和聚焦振幅与理论值相差较小。随着谱峰频率的增大,高频强非线性短波所占能量比例增大,波群的非线性特征增强,因此真实聚焦振幅增大。而对于fp=1.03 Hz工况,其能量大部分集中在高频组分,组成波的非线性很强,聚焦位置处的波陡超过了有限水深波的极限破碎波陡,因此部分高频组分出现了破碎,导致波群聚焦位置处的振幅偏小、聚焦位置提前。以上说明谱峰频率较小时波群能量分布较均匀,而谱峰频率较大时波群能量集中在高频部分且内部非线性较强。
图17 不同谱峰频率下畸形波聚焦位置处自由面历时变化

Fig.17 Time-histories of free surface elevation for freak wave with different peak frequencies

3.5 水深的影响

色散性是波浪的重要特性,本小节探究不同色散系数对畸形波产生过程中自由面变化的影响,波群中考虑固定波数的组成波,通过改变水深改变波群的色散范围。选取实验室常用的JONSWAP谱,谱宽df=0.6 Hz,频谱范围为0.53~1.13 Hz,组成波个数N=29,聚焦振幅Af=0.09 m,聚焦位置xf=10.2 m,聚焦时间Tf=20 s。分别选取水深h=0.3、0.4、0.5、0.6及0.7 m的情况进行模拟研究,工况参数见表4图18为不同水深下的JONSWAP 输入谱型。
表4 工况参数

Table 4 Parameters of different cases

工况
编号
频谱范围
f/Hz
谱峰频率
fp/Hz
聚焦振幅
Af/m
水深
h/m
波陡 色散参数
1 0.53~1.13 0.83 0.09 0.3 0.32 [0.62,1.66]
2 0.53~1.13 0.83 0.09 0.4 0.29 [0.73,2.12]
3 0.53~1.13 0.83 0.09 0.5 0.28 [0.83,2.60]
4 0.53~1.13 0.83 0.09 0.6 0.27 [0.93,3.10]
5 0.53~1.13 0.83 0.09 0.7 0.26 [1.03,3.61]
图18 不同水深下的JONSWAP谱型

Fig.18 JONSWAP spectra with different water depths

图19为不同水深下畸形波的自由面变化。由图19(a) 可知,水深h=0.7 m时畸形波的聚焦振幅最大,其次分别为h=0.6、0.5、0.4和0.3 m,而水深h=0.3 m时,波群明显出现了破碎,聚焦振幅最小。对比聚焦时刻,水深h=0.7 m时波浪聚焦时刻距离输入时刻最近,其次分别为h=0.6、0.5、0.4和0.3 m。由图19(b) 可知,不同水深下畸形波的波面聚焦结果差异较小,除了聚焦位置处的峰值存在差异外,最大波峰两侧的次峰形态稍有差异,水深h=0.4 m时最大波峰两侧的次峰出现的位置较远并且峰值最大。
图19 不同水深下畸形波的自由面变化

Fig.19 Time-histories of free surface elevation for freak wave at different water depths

表4可知,工况 2—工况5谱峰波陡的差异很小,也是水深h=0.4、0.5、0.6和0.7 m时波浪聚焦位置处自由面形态差异较小的原因。对于水深h=0.3 m,波群中组成波的色散参数范围很小,这说明波群中组成波的频率较高,波列的非线性很强,聚焦位置处的波陡超过了有限水深波的极限破碎波陡,因此部分高频组分出现了破碎,导致聚焦位置处的振幅较小,聚焦位置和时间无明显规律。比较工况2—工况5,波群中组成波的色散参数范围逐渐增大,表明波群中高频组成波的比例增大,波群的非线性增强。另外由波浪的色散性可知,高频波浪传播速度较慢,低频波浪传播速度较快,这也是水深h=0.7 m情况下畸形波聚焦时刻距离输入时刻最近而聚焦位置处的振幅最大的原因。由上可知水深影响波群的色散性,不同水深下波群的传播速度不同,进而影响畸形波的波面形态。

4 结论

本文自主编程建立了黏性流无反射数值波浪水槽,分析了畸形波传播过程中自由面及内部频谱结构的演变,揭示了非线性对波群传播过程自由面变形及频谱结构的影响规律。在此基础上,开展了全面的数值试验,分析了不同参数对畸形波生成及发展的影响规律,探究了各个参数对计算结果的敏感性。主要结论如下:
(1)畸形波传播过程中,由于波-波间非线性相互作用,波群能量由主频向高频和低频转移。波群中的非线性项,能量主要分布在低频和2倍主频范围内。对于高频自由伪谐波,其传播较快且在聚焦时刻相位与主频波一致,导致波群聚焦位置处的波高大于实际波高,而高频约束谐波传播较慢,主要以波群后面的小尾波形式出现,并不影响波群聚焦位置主峰处的波面形态。
(2)不同的入射谱型导致波群内部波能分布的不同,进而影响畸形波生成的形态和特性,JONSWAP谱可产生能量更为集中的畸形波。波群中的组成波个数影响畸形波的聚焦重现周期,过少的组成波个数可能导致次聚焦波的出现。
(3)波群频谱宽度直接影响畸形波的非线性,频谱宽度越窄,波列组成波的振幅越大,波群聚焦位置处的波幅值也会增大。波群在有限水深情况下,并非完全满足谱宽越窄、聚焦位置波幅越大的规律,其还与水深、谱峰频率等因素有关。
(4)谱峰频率较小时波群能量分布较均匀,而谱峰频率较大时波群能量集中在高频部分且内部非线性更强。水深影响波群的色散性,不同水深下波群的传播速度不同,进而影响畸形波的波面形态。
综上所述,畸形波生成及演化的影响因素复杂,需要综合考虑多种因素准确模拟和预测畸形波。本文研究有助于理解畸形波的形成机制和演化规律,为实验室畸形波生成提供了参考。
[1]
DRAPER L. ‘Freak’ Ocean Waves[J]. Weather, 1966, 21(1): 2-4.

DOI

[2]
CHIEN H, KAO C C, CHUANG L Z H. On the Characteristics of Observed Coastal Freak Waves[J]. Coastal Engineering Journal, 2002, 44(4): 301-319.

DOI

[3]
DYSTHE K, KROGSTAD H E, MÜLLER P. Oceanic Rogue Waves[J]. Annual Review of Fluid Mechanics, 2008, 40: 287-310.

DOI

[4]
DIDENKULOVA E G, PELINOVSKY E N. Freak Waves in 2011—2018[J]. Doklady Earth Sciences, 2020, 491(1): 187-190.

DOI

[5]
DIDENKULOVA E, DIDENKULOVA I, MEDVEDEV I. Freak Wave Events in 2005—2021:Statistics and Analysis of Favourable Wave and Wind Conditions[J]. Natural Hazards and Earth System Sciences, 2023, 23(4):1653-1663.

[6]
PELINOVSKY E, KHARIF C. Extreme Ocean Waves[M]. Switzerland: Springer, 2008.

[7]
张亚荣. 基于双波群聚焦的畸形波形成机理分析[D]. 大连: 大连理工大学, 2019.

(ZHANG Ya-rong. Analysis of the Formation Mechanism of Rogue Waves Based on Double Wave Group Focusing[D]. Dalian: Dalian University of Technology, 2019. (in Chinese))

[8]
CHAPLIN J R. On Frequency-focusing Unidirectional Waves[J]. International Journal of Offshore and Polar Engineering, 1996, 6(2):131-137.

[9]
BRANDINI C, GRILLI S. Modeling of Freak Wave Generation in a 3D-NWT[C]// Isope International Ocean and Polar Engineering Conference. Stavanger, Norway. June 17-22, 2001:ISOPE-I.

[10]
黄国兴. 畸形波的模拟方法及基本特性研究[D]. 大连: 大连理工大学, 2002.

(HUANG Guo-xing. Simulation Methods and Fundamental Characteristics of Rogue Waves[D]. Dalian: Dalian University of Technology, 2002. (in Chinese))

[11]
ZHAO X Z, HU C H, SUN Z C. Numerical Simulation of Extreme Wave Generation Using VOF Method[J]. Journal of Hydrodynamics, Ser B, 2010, 22(4): 466-477.

[12]
NING D Z, ZANG J, LIU S X, et al. Free-surface Evolution and Wave Kinematics for Nonlinear Uni-directional Focused Wave Groups[J]. Ocean Engineering, 2009, 36(15/16): 1226-1243.

DOI

[13]
裴玉国. 畸形波的生成及基本特性研究[D]. 大连: 大连理工大学, 2008.

(PEI Yu-guo. The Generation of Freak Waves and Its Behaviors[D]. Dalian: Dalian University of Technology, 2008. (in Chinese))

[14]
赵西增. 畸形波的实验研究和数值模拟[D]. 大连: 大连理工大学, 2009.

(ZHAO Xi-zeng. Experimental and Numerical Study of Freak Waves[D]. Dalian: Dalian University of Technology, 2009. (in Chinese))

[15]
刘赞强. 畸形波模拟及其与核电取水构筑物作用探究[D]. 大连: 大连理工大学, 2011.

(LIU Zan-qiang. Freak Wave Simulation and Its Interaction with Nuclear Power Station Ingarage Structure[D]. Dalian: Dalian University of Technology, 2011. (in Chinese))

[16]
崔成, 张宁川, 郭传胜, 等. 变水深对畸形波及其时频能量谱的影响[J]. 海洋学报, 2011, 33(6): 173-179.

(CUI Cheng, ZHANG Ning-chuan, GUO Chuan-sheng, et al. Impact of Water Depth Variation on Simulated Freak Waves and Their Time-frequency Energy Spectrum[J]. Acta Oceanologica Sinica, 2011, 33(6): 173-179. (in Chinese))

[17]
扈喆. 畸形波数值模拟及其对结构物的作用[D]. 上海: 上海交通大学, 2015.

(HU Zhe. Numerical Simulation of Rogue Waves and Their Action on Structures[D]. Shanghai: Shanghai Jiao Tong University, 2015. (in Chinese))

[18]
YU Xiang-jun, LI Qing-hong, WANG Hua. Phase-Focusing Model of Freak Wave Based on Boussinesq Equation[C]// Proceedings of 2019 International Conference on Modeling, Analysis, Simulation Technologies and Applications(MASTA 2019):Hangzhou, China.May 26-27, 2019.

[19]
王磊. 不同波况下畸形波发生概率的模拟研究[D]. 大连: 大连理工大学, 2021.

(WANG Lei. Investigations on Occurrence Probability of Freak Waves under Various Wave Conditions[D]. Dalian: Dalian University of Technology, 2021. (in Chinese))

[20]
WANG L, DING K, ZHOU B, et al. Quantitative Prediction of the Freak Wave Occurrence Probability in Co-propagating Mixed Waves[J]. Ocean Engineering, 2023, 271: 113810.

DOI

[21]
崔成, 潘文博. 短峰畸形波生成、演化过程的外部特征研究[J]. 海洋学报, 2023, 45(7): 79-89.

(CUI Cheng, PAN Wen-bo. Study on the External Characteristics of the Generation and Evolution of Short-crested Freak Waves[J]. Haiyang Xuebao, 2023, 45(7): 79-89. (in Chinese))

[22]
刘震, 吴永顺, 周文俊, 等. 畸形波物理水池模拟与传播特性研究[J/OL]. 华中科技大学学报(自然科学版). 2024.

(LIU Zhen, WU Yong-shun, ZHOU Wen-jun, et al. Physical Wave Basin Simulation and Propagation Characteristics of Rogue Waves[J/OL]. Journal of Huazhong University of Science and Technology(Natural Science Edition), 2024. (in Chinese))

[23]
KHARIF C, PELINOVSKY E. Physical Mechanisms of the Rogue Wave Phenomenon[J]. European Journal of Mechanics - B/Fluids, 2003, 22(6): 603-634.

DOI

[24]
CHOI J, YOON S B. Numerical Simulations Using Momentum Source Wave-maker Applied to RANS Equation Model[J]. Coastal Engineering, 2009, 56(10):1043-1060.

DOI

[25]
LIU X, LIN P, SHAO S. ISPH Wave Simulation by Using an Internal Wave Maker[J]. Coastal Engineering, 2015, 95: 160-170.

DOI

[26]
ZHAO X, HU C. Numerical and Experimental Study on a 2-D Floating Body under Extreme Wave Conditions[J]. Applied Ocean Research, 2012, 35: 1-13.

DOI

[27]
LI M, ZHAO X, YE Z, et al. Generation of Regular and Focused Waves by Using an Internal Wave Maker in a CIP-based Model[J]. Ocean Engineering, 2018, 167: 334-347.

DOI

[28]
ADCOCK T A A, TAYLOR P H. Non-linear Evolution of Uni-directional Focussed Wave-groups on a Deep Water: A Comparison of Models[J]. Applied Ocean Research, 2016, 59: 147-152.

DOI

[29]
BALDOCK T E. Non-Linear Transient Water Waves[D]. London,UK: Imperial College of Science, Technology and Medicine, 1995.

[30]
BALDOCK T E, SWAN C, TAYLOR P H. A Laboratory Study of Nonlinear Surface Waves on Water[J]. Philosophical Transactions of the Royal Society of London Series A: Mathematical, Physical and Engineering Sciences, 1996, 354(1707): 649-676.

文章导航

/