引 言
近水面波与水下物体相互作用的问题在海洋工程、船舶运输、军事等领域中有着广泛的应用．通过多重浸没点源分布的数值模拟,探索点源与自由表面相互作用下复杂的水动力学特性,对海上工程作业的安全与质量有重要的意义，对水翼船与水下航行器的优化研究也具有重要价值．
在过去几十年里,学者们提出了大量的理论、实验、数值方法,研究了基于源分布的水动力学问题．1932年,Havelock[1-2]引入自由表面Green函数来描述点源和流体结构间相互作用产生的自由表面波．1960年,Wehausen和Laitone[3]整理了自由面Green函数标准表达形式,推导出奇异积分方程的解,将移动点源引起的流动问题扩展至任意流体深度．由于Green函数的表达形式过于复杂,1981年,Noblesse[4]分别对Michell、Havelock、Peters函数进行改进,推导出三种可供选择的单积分表达形式,研究每种表达形式适用点源周围近场、远场等不同空间．1989年,Forbes[5]通过Newton迭代的方法,数值模拟了由浸没在均匀流体中的点源引起的非线性三维自由表面波高．Baar和Price[6]将两个互补的Neumann级数展开式与Newman[7]的非震荡Chebyshev展开式相结合,数值求解Kelvin波源势．2001年,Bal等[8]采用四阶差分格式以保证自由面数值解稳定,研究了以恒定速度运动的二维或三维空泡水翼的水动力学性态．2010年,Yeung和Bouger[9]使用B样条函数和波动特征函数来近似近场和远场波浪运动．2013年，Chen[10-12]提出用带耗散源Green函数数值方法计算奇异积分,模拟出了单点源及水翼的自由表面波高．2011年,Yang等[13]基于双流体公式研究定常和非定常自由表面流,模拟出波浪在水槽中、漫过固定的二维甲板及冲击固定三维物体等情况下的自由表面波高．2017年,Chen[14]将边界元法与Hess-Smith面板积分公式相结合,预测出三维升力体的水动力性能．2018年,钟万勰等[15]提出保辛摄动法,计算了一般水深周期行波解,并模拟出水面尖锐的周期行波解．
综上所述,以上文献考虑各种数值方法求解浸没单点源自由面波高,研究基于源分布下水动力学性态,但很少有文章考虑多重浸没点源分布下自由面波高的水动力学性态．由于多重点源间复杂的相互作用关系,不同距离分布、不同点源个数下自由面波高等内部流场机制皆产生复杂变化．鉴于此,本文基于带耗散源Green函数的去奇异化数值方法,建立了多重点源物理模型,数值模拟了多重点源自由表面波高,在多点源复杂的相互作用关系中寻找有规律和有意义的水动力学特点．
本文第1节为模型的建立；第2节为带耗散源Green函数和控制方程的推导；第3节为多重点源带耗散Green函数的离散；第4节为数值模拟及分析；第5节总结了多重点源分布水动力学性态规律．
1 物 理 模 型
假设流体不可压、无黏和无旋．在平静水面上定义直角坐标系Ox*y*z*,如图1所示,考虑在三维无限水深的均匀流中,浸没m个源强为M的点源
是静止水面到源点的垂直距离,η*是自由面波高．
流体域内总速度势Φ满足
Φ=Ux*+φ*,
(1)
其中,φ*表示流体扰动速度势,U是沿x*方向的水平均匀来流速度．
图1 多重点源物理模型
Fig. 1 The multi-point source physical model
φ*满足Laplace方程：
∇2φ*=0．
(2)
带耗散自由表面Bernoulli方程[12]可以表示为
∇Φ|2+gη*+μ′φ*=0, z*=0,
(3)
其中,g为重力加速度,μ′为耗散系数．
自由面线性运动学边界条件为
(4)
定义以下无纲量变换:
(5)
(6)
其中,h为特征长度,ε为无量纲源强．
通过线性化方法简化式(3),将结果与式(4)结合,得到线性耗散自由表面边界条件：
(7)
其中
(8)
2 带耗散源Green函数及控制方程
考虑单点源形式下的耗散自由表面Havelock源Green函数gμ．点源位于(ξ,σ,ζ),ζ<0,Green函数由Rankine源源势、像势和奇异波积分k组成,表示为
(9)
对于有限多个点源Qj(ξj,σj,ζj), j=1,2,3,…,m,将带耗散Green函数gμ推广至多重点源的形式:
(10)
其中
R=|(x,y,z)-(ξj,σj,ζj)|, R′=|(x,y,z)-(ξj,σj,-ζj)|．
由Gμ满足线性耗散自由表面边界条件(7),得
(11)
即
(12)
根据文献[12],使用Cauchy积分得到1/R和1/R′在复数域上的表达形式,将其代入式(12)，其中,点源
的无量纲坐标为Qj(ξj,σj,-1),得到
(13)
计算可知
(14)
根据式(13)、(14),多重点源正则波积分Kμ的表达式为
(15)
其中,i是虚数单位,∇2Kμ=0．
由Forbes[5]的推导可得,三维模型下扰动速度势φ满足
(16)
由式(3)、(5)、(6)、(8),在微小耗散μ下,得到无量纲化后的线性波高表达式为
(17)
将式(16)代入式(17),得到三维自由面波高公式：
(18)
3 多重点源带耗散Green函数的离散
用耗散Green函数模拟三维自由表面波高,根据式(10)和(15),求解在自由面z=0情况下的
(19)
在正则波积分Kμ的计算域中,选取以l=1,2,…,N1+1和n=1,2,…,N2+1为序列的网格点(kl,θn),对正则波积分Kμ进行网格离散,∂Kμ/∂x近似可得
(20)
通过式(18)和(20),多点源三维自由表面的波高可近似为
(21)
其中
(22)
4 数 值 算 例
4.1 数值验证
为了验证文中所采用数值方法的准确性和有效性,设点源个数m=1,且点源位于水下(0,0,-1)处．选取参数ε=2.6,μ=0.001,Fr=0.7,得到y=0平面的自由表面波高．将本文得到的单点源波高与Forbes[5]的数值解进行比较,如图2,本文波高的数值解用实线表示,Forbes[5]的数值解用符号“○”表示．从图中可以看出,结果吻合很好,表明本文采用的数值方法是有效的．
图2 本文单源波高与Forbes[5]波高的比较结果
Fig. 2 The comparison of single point wave elevations between the present method and the Forbes[5] method
4.2 数值结果与讨论
4.2.1 双点源二维波高数值模拟
本小节数值模拟浸没的双点源在不同距离下的自由表面波高．即以单点源波高作对比,在此基础上观察当增加的点源位于不同位置时,和单点源波高相比波形的变化规律．
选取参数ε=2.6,μ=0.001,Fr=0.7．考虑双点源,即m=2的情况下,固定点源Q1的位置为(0,0,-1),研究点源Q2(ξ2,0,-1)中ξ2取不同位置时自由表面波高的变化．
(a) 增加点源位于单点源波峰附近的自由表面波高
(a) The free surface wave elevations with an added point source located near the wave crest formed by a single point source
(b) 增加点源位于单点源波谷附近的自由表面波高
(b) The free surface wave elevations with an added point source located near the wave trough formed by a single point source
图3 双点源不同位置自由面波高
Fig. 3 The free surface wave elevations with double point sources at different locations
图3(a)、(b)的数值结果表明:
在图3(a)中,两个点源的源强度叠加后产生的波形,波高大于单点源产生的波高．取ξ2分别为0.9,0.6,0.3,当点源Q2的位置越靠近Q1单点源波的波峰时,双点源叠加效果越明显,波峰和波谷的绝对值越大．其中,ξ2=0.3,即在位于Q1波峰附近时,波峰波谷的振幅增加最显著．且随着双点源之间距离的增大,波高逐渐向后推移．
在图3(b)中,点源Q2位于Q1单点源波的波谷附近,ξ2分别取1.4,1.5,1.6时,双点源波高呈现较为明显的变化．前一部分波的波峰波谷振幅仍然增加,从点源Q2位置开始,沿来流方向自由面波高出现明显下降; 随距离增加,波以较小振幅震荡,且波高基本趋于自由面水平线上．当点源Q2的位置越靠近Q1单点源波的波谷时,双点源之间的干扰越明显,点源Q2后出现的第一个波峰峰值越小．
4.2.2 多点源二维波高数值模拟
本小节数值模拟多点源在不同距离下的自由表面波高．即从原点正下方开始,沿x轴方向浸没多个点源,研究当点源位于不同位置时,波高呈类似周期性变化的现象．
选取参数ε=2.6, μ=0.001, Fr=0.7．图4表示点源个数m=5,且点源Qj(ξj,0,-1)中ξ1,ξ2,ξ3,ξ4,ξ5分别取0,5.34,11,16.7,22.4时,自由表面的波高．从图4可以看出,在布置点源的区间内,从原点后的第一个波谷开始,波高呈类似周期性变化．每两个波峰波谷为一个周期,周期约为5.5,与相邻两个点源之间的距离近似相同,且每个周期内两个波峰之间的波谷振幅较小．
图4 m=5, y=0中心平面自由表面波
Fig. 4 The center plane free surface wave for m=5, y=0
图5 m=8, y=0中心平面自由表面波
Fig. 5 The center plane free surface wave for m=8, y=0
选取参数ε=2.6,μ=0.001,Fr=0.7．图5表示点源个数m=8,且点源Qj(ξj,0,-1)中ξ1,ξ2,ξ3,ξ4,ξ5,ξ6,ξ7,ξ8分别取0,2.5,5,7.78,10.52,13.2,15.99,18.6时,自由表面的波高．从图5可以看出,在布置点源后的第一个波峰开始,波高有规则地进行震荡,产生类似于正弦函数图像的周期性变化．波峰与波谷绝对值较为接近,周期约为2.5,与相邻两个点源之间的距离近似相同．
4.2.3 多点源自由面波高数值模拟
本节考虑在y轴正下方浸没多个等距离点源,研究其自由表面波高及等高线的变化规律．其中点源浸没深度h=-1,Froude数Fr=0.6,耗散系数μ=0.001．
(a) 自由表面波高模拟 (b) 自由表面等高线模拟
(a) The free surface wave elevation simulation (b) The free surface contour line simulation
图6 多点源自由表面波
Fig. 6 Multi-point source free surface waves
注 为了解释图中的颜色，读者可以参考本文的电子网页版本．
图6表示取m=4,源强ε=5.2 时产生的三维自由表面波．其中点源位于Q1(0,3,-1),Q2(0,1,-1),Q3(0,-1,-1),Q4(0,-3,-1)．此时,四个独立点源产生的波形相互叠加,自由表面形成如图6(a)所示的规则波,波峰向两侧水平拉伸．波面与Kelvin船波不同,以横波的形式沿x轴方向传播．图6(b)中,研究自由表面等高线变化,由于多点源的作用,自由面等高线交汇分布,点源间横波呈互相重合状态,散波呈均匀状态逐渐向外扩散,且中心区域等高线密集．由此，可以清晰看出波形分布规则．
5 结 语
本文将带耗散源Green函数去奇异化数值方法进行扩展,以便适用于多重点源的数值计算．数值研究浸没的多个点源在不同排列下波高及波形的变化规律．文中自由面的波高与Forbes[5]的数值结果进行比较,验证了本文所采用的数值方法的有效性和准确性．通过数值模拟可以发现:
1) 在双点源情况下，当增加的点源位置越靠近单点源波的波峰时,双点源叠加效果越明显; 当增加的点源位置越靠近单点源波的波谷时,双点源之间相互干扰,波峰波谷迅速下降到趋于自由面水平线上,且以小振幅缓慢震荡．
2) 在多点源情况下，点源以不同的距离纵向排列,y=0中心平面波高出现类似周期性变化的现象．波高的周期与点源间距离的大小基本趋于一致．
3) 在多点源情况下，点源等距横向排列时,多个点源自由表面波高相互叠加,以横波形式规则地向后传播,向两侧水平扩散．等高线交汇分布,波面清晰．
[1] HAVELOCK T H. Wave resistance[J]. Proceedings of the Royal Society of London, 1928, 118(779): 24-33.
[2] HAVELOCK T H. The theory of wave resistance[J]. Proceedings of the Royal Society A, 1932, 138(835): 339-348.
[3] WEHAUSEN J V, LAITONE E V. Surface waves[J]. Handbook of Physics, 1960, 9(68): 446-778.
[4] NOBLESSE F. Alternative integral representations for the Green function of the theory of ship wave resistance[J]. Journal of Engineering Mathematics, 1981, 15(4): 241-265.
[5] FORBES L K. An algorithm for 3-dimensional free-surface problems in hydrodynamics[J]. Journal of Computational Physics, 1989, 82(2): 330-347.
[6] BAAR J M, PRICE W G. Evaluation of the wavelike disturbance in the Kelvin wave source potential[J]. Journal of Ship Research, 1988, 32(1): 44-53.
[7] NEWMAN J N. Evaluation of the wave-resistance Green function, part 1: the double integral[J]. Journal of Ship Research, 1987, 31(2): 79-90.
[8] BAL S, KINNAS S A, LEE H. Numerical analysis of 2-D and 3-D cavitating hydrofoils under a free surface[J]. Journal of Ship Research, 2001, 45(1): 34-49.
[9] YEUNG R W, BOUGER Y C. A hybrid integral-equation method for steady two-dimensional ship waves[J]. International Journal for Numerical Methods in Engineering, 2010, 14(3): 317-336.
[10] CHEN Z M. A vortex based panel method for potential flow simulation around a hydrofoil[J]. Journal of Fluids and Structures, 2012, 28: 378-391.
[11] CHEN Z M. A dissipative free-surface approach for an underwater body in a uniform stream[J]. Journal of Marine Science and Technology, 2013, 18(2): 182-191.
[12] CHEN Z M. Harmonic function expansion for translating Green functions and dissipative free-surface waves[J]. Wave Motion, 2013, 50(2): 282-294.
[13] YANG C, HUANG F, WANG L, et al. Numerical simulations of highly nonlinear steady and unsteady free surface flows[J]. Journal of Hydrodynamics, 2011, 23(6): 683-696.
[14] CHEN Z M, LIU Q, PRICE W G. Straightforward discretisation of Green function and free-surface potential flow around a three-dimensional lifting body[J]. Journal of Marine Science and Technology, 2017, 22: 149-161.
[15] 钟万勰, 吴锋, 孙雁, 等. 保辛水波动力学[J]. 应用数学和力学, 2018, 39(8): 855-874.(ZHONG Wanxie, WU Feng, SUN Yan, et al. Symplectic water wave dynamics[J]. Applied Mathematics and Mechanics, 2018, 39(8): 855-874.(in Chinese))