基于加速度势的船舶波浪力数值模拟-李志富.pdf

《基于加速度势的船舶波浪力数值模拟-李志富.pdf》由会员分享，可在线阅读，更多相关《基于加速度势的船舶波浪力数值模拟-李志富.pdf（11页珍藏版）》请在得力文库 - 分享文档赚钱的网站上搜索。

1、57卷第3期(总第219期) 中 国 造 船2016年9月 SHIPBUILDING OF CHINA、，0157 No3(Serial No219)Sep2016文章编号：1000-4882(2016)03009811基于加速度势的船舶波浪力数值模拟李志富，任慧龙+，石玉云(哈尔滨工程大学船舶工程学院，哈尔滨150001)摘 要波浪中航行船舶所受波浪力的时域模拟通常在大地坐标系下进行，为得到每一时刻作用于船体表面的水动压力，需要在每一时间步求解速度势的时问偏导数。为避免采用差分法直接求解该偏导数所导致的数值不稳定，文中引入了流体加速度势的概念。其次，据速度势所满足的场方程和边界条件，推导了加

2、速度势所应满足的场方程和边界条件，并建立了求解加速度势的边界积分方程。再次，根据加速度势的物理意义，推导了由加速度势表达的Bernoulli方程，以求解每一时刻船体表面上的压力值。此外，为数值求解速度势、加速度势所满足的边界积分方程，采用了八节点高阶曲面元进行了边界积分方程的数值离散，以提高计算收敛速率。最后，以半球和WigleyI船型为例，计算了相关的水动力系数并与解析解和试验值进行了对比，发现符合良好。关 键 词：大地坐标系；格林函数；加速度势；波浪力；有限差分中图分类号：0352 文献标识码：A0 引 言准确预报船舶在波浪中航行时的运动响应和所承受的波浪载荷对于船舶的设计和营运至关重要。

3、船舶的大幅摇荡运动会严重影响船舶的舒适性，同时过大的波浪载荷会导致船体结构的破坏，甚至是船体的断裂。在势流理论框架下，通常在时域上求解航行船舶的相关水动力问题11。Lin2】在大地坐标系下，以时域自由面格林函数为积分核，建立了流体速度势所满足的边界积分方程，并对航行船舶的阻力、线性以及非线性波浪力进行了数值模拟和分析。此后，Datta31也做了类似的研究工作，但是对边界积分方程采用了连续的Bspline进行离散求解。最近，Datta【4通过对外飘网格进行直壁修正，在大地坐标系下对一系列的渔船进行了运动响应的时域分析，并详细研究了外飘网格修正对计算结果的影响。在以上的研究工作中，求解Bcmoul

4、li方程中速度势的时间偏导数，均采用差分法求解，其数值计算精度受所取时间步长的影响较大，特别在求解船舶大幅运动时，水线附近会出现干、湿网格交替的情形。针对以上问题，为了直接求解速度势的时间偏导数，可通过边界积分方程直接求解加速度势。Wu【5建立了以加速度势为未知量的边界积分方程，并推导了相关的边界条件。Huang【6】采用加速度势，基于弱散射假定，求解船舶在波浪中航行时所受到的非线性波浪载荷。Zhou7】在时域内求解Spar平台所收稿日期：20151028：修改稿收稿日期：20160830基金项目： 国家重点基础研究发展计划(2011CB3703)；国家自然科学基金(51109046)万方数据

5、57卷第3期(总第219期) 李志富，等：基于加速度势的船舶波浪力数值模拟受到的全非线性波浪力，其边界积分方程采用多节点高阶边界元进行离散。上述计算所采用的积分核均为Rankine源，且较为可靠的计算大部分为无航速的海洋平台。利用时域自由面格林函数求解船舶波浪力时，格林函数的计算直接影响了整个时域模拟的精度和效率。Beck【8】根据Green函数的变化特性，采用小参数级数展开、大参数渐进展开的方法对格林函数在不同区域上的值进行了计算。Newman【9】采用切比雪夫级数逼近的方法对格林函数进行了计算。Lin和Yue【2】综合了级数展开、渐进展开和二维级数逼近进行了计算。通过对无量纲形式格林函数形

6、式的分析，Clement oo证明了格林函数实际上为一四阶微分方程的解。基于此，Li【ll】采用一种改进的精细时程积分算法对格林函数进行了精确求解。为准确、高效地求解船舶在波浪中运动时所受到的波浪力，本文首先在大地坐标系下建立了流体速度势所满足的边界积分方程。其次，基于速度势的边界条件较为详细地推导了流体加速度势所应满足的定解条件，并建立了相应的边界积分方程。继而，采用改进的精细时程积分算法进行格林函数数值计算，并提出了相应的制表插值策略以提高计算效率。为加快计算的几何收敛速率，采用八节点高阶曲面元实现了边界积分方程的数值离散。最后，以半球和WigleyI船型为例，计算了相关的水动力系数，并与

7、解析解和试验值进行了对比，数值解与解析解和试验值符合良好。1 基础理论为合理描述船舶的运动和流场，如图1所示，定义两个坐标系：大地坐标系oxyz，oxy平面位于静水面上，z轴正向垂直向上，X轴的正向与船舶的航行方向同向，该坐标系用于描述船舶的运动和流场；固船坐标系(与船体固结的坐标系)OtXYz，坐标原点位于船舶重心，在初始时刻固船坐标系和大地坐标系坐标轴相互重合。0图1坐标系定义不意图11速度势定解条件基于势流理论的假定12】，将波物相互作用问题由求解流体速度矢量转化为了求解流体速度势的问题。由于所关心的问题为船舶在波浪中航行时所受到的波浪力，故将速度势作如下分解少(p，t)=旃+ (1)式

8、中，办为己知的入射波速度势，声为待求解的扰动速度势。在流体域内，速度势应满足Laplace方程V2(p，t)=0 in D (2)万方数据100 中 国 造 船 学术论文征目由向上，速度势满足线性化的运动学边界条件和动力学边界条件髦+g掣：o onl1辞Oo+p二=- n 、-t2。6包。7在物面上，速度势满足法向不可穿透条件娑：K一娑。n&在有限的时间内，物体的扰动在无限远处为零熙h警，-0 ror尺斗一对于时域内求解，速度势还应满足以下初始条件礼=o=。，乱=。when蜓。12速度势边界积分方程根据速度势定解条件，选取时域自由表面格林函数作为积分核13】：G(p，t；q，f)=万(fr)x

9、Go+日(ff)0Go(p，q)=1r-1r(3)(4)(5)(6)(7)(8)G(p，t；q,r)=2e坼厶(艘)磊s洒_i(ff)救 (9)式中，Go(p，q)为瞬时效应，G(p，q，tf)为记忆效应，万(f)为脉冲函数，H(t)为单位阶跃函数，p(x，Y，z)为场点，g(孝，7，f)为源点，g(孝，rl，一f)为源点关于自由面的镜像点，R2=(X-善)2+()，一77)2为场点和源点之间的水平距离，r2=R2+fzf12，“=R2+fz+12。对扰动速度势和时域自由面格林函数应用格林定理，并对时间变量在(0，f)上做积分，可得求解流体速度势的混合分布模型。由于在计算压力分布的时候需要求解

10、速度势的空间偏导数，故在物体内部4兀K()一嚣)_k叱毒(击)心+：dfk叱一筹呶一扣“仃(瓦OG慨 )13加速度势定解条件 -为了采用相同的边界积分方程求解流体加速度势，要求加速度势和速度势除物面条件以外，满足相同的边界条件和初始条件。本文以Kukkanen141的工作为基础，推导了流体加速度势满足的物面条件。万方数据57卷第3期(总第219期) 李志富，等：基于加速度势的船舶波浪力数值模拟 101物体表面上任意一点的法向速度可以写为K=U刀=(U+国，)刀式中，U和国分别为在固船坐标系下给出的平移速度和旋转角速度矢量，矢量(相对于旋转中心)，以为物面上一点的法向量。对式(4)两端在大地坐标

11、系下取绝对时间导数如下导：采(uV萌)以出锄m L “7 J上式左端可以写为旦f型1：以型+塑v出L锄 出 出由dn(it=国以以及全导数的定义可得苦(辨吐警巾V”V小(酬V通过简单直接的推导，上式最终可以写为(12)r为物体表面上一点的位置(13)(14)(15)旦，型1：旦f，型兰鲤 (16)dtananI(itan式中，缈=ddt即为加速度势。由式(12)可以将式(13)右端第一项展开为导(m)=亭(即+缈r)力娟+西，)一(国)一 )根据全导数的定义可以将式(12)右端第二项展开为未(V旃力)=av_出A旃栉+V旃岩； (。8)：(盟+(以+国，)v)v办)卅(国v办)以综合式( 、式

12、( 、式(,-3t16)17) 18)可得加速度势应满足的物面条件警叫+西，巾一(缈H)。以一(旦呈堕+(H+国x，)v)v办)以+(国v识)甩 (19)14格林函数的精细计算通过引入新的参变量，可以将时域自由表面格林函数的记忆项写为 o(pq)=2廊()(20)6()=r瓶siIl(屁)xe一舡厶兄(廊)卜 (21)式中，2=一(z+f)，7，f=(卜t)届7。Clement151发现可以通过求解以下四阶常微分方程的方式来求解形如式(21)的无穷积分：万方数据中 国 造 船 学术论文0(4)+f0(3)+ff24+4t16(2)+74f0(1)+946：0(22)、 ，“【ll】提出了一种改

13、进的精细时程积分算法实现了对格林函数的精确计算。通过类似的办法，可以求得格林函数对R的偏导数岛和对z偏导数G，。由于格林函数及其空间偏导数是通过求解形如式(22)所示的四阶常微分方程得到的，因此也会同时得到格林函数及其空间导数关于时间f的一阶偏导数，即在给定的数据点上同时已知函数值及函数的一阶时间偏导数。基于此，本文采用了双三次SplineHermite插值方式【16】实现了对格林函数精确插值计算。通过数值试验发现时间步长取Ar=0005，0，001)时=O001，00l，1)时=001，已经能够满足工程中的精度要求。15边界积分方程的离散求解为加快边界积分方程的几何收敛速度，本文选取能够适用

14、于复杂几何形体的八节点高阶曲面元进行边界积分方程的离散求解，大大减少了格林函数的计算次数和系数矩阵的维数。通过引入八节点插值形函数Ni=三(1+白孝)(1+，乃叩)(一l+彭孝+仉刁)j=1,3,5,7 (23)三(1+彭孝+仉77)(1一芎，72一矿孝2) ，=2，4，6，8面元内任意一点的位置坐标、分布源强可以表达为形函数插值的形式x,y，z)=(善，叩)(_，Yj，乃)，=l8盯(x)=(善，印)巳(24)(25)将以上两式代入边界积分方程式(11)便可以求解分布源强，通过求解边界积分方程式(10)可以进一步求解速度势。由于在求解分布源强的同时已经计算了格林函数偏导数的积分，因此本文采用

15、积分表达式来求解速度势的空间偏导数。16船体表面压力的计算根据Bernoulli方程，流场中任意一点的压力可以写为见(雕)=叫等+三啊捌 (26)由式(1)以及加速度势妒可以将上式进一步写为Pd(p，f)=一p(缈一UV+去IV函lz) (27)p。(p，f)叫警 )Ps(p，t)=一pgz (29)式(27)保留扰动速度势的平方项是为了合理地考虑航速效应【2J，项的影响。对应的波浪力为E 2LPa以ds，坞2 J岛(rxn)ds互2I，(t)PI以ds，M 2kp-(rxn)凼对于零航速线性问题，则不计入该(30)(31)万方数据57卷第3期(总第219期) 李志富，等：基于加速度势的船舶波

16、浪力数值模拟只2Lf)以即ds，丝=f(f)见(，以)ds (32)对于船舶运动和载荷的线性模拟，静水恢复力项采用静水恢复力系数进行表达，恢复力的具体形式可参见文献13】。此外，由加速度势的物面边界条件可知，在每一时间步的求解过程中，运动方程的右端含有待求解的船体加速度，需要进行迭代求解，根据数值试验可知，每一时间步内进行三次迭代计算即可满足计算精度和收敛性的要求。2模型有效性验证21算法有效性验证为验证本文所提出的基于加速度势的波浪力数值模拟方法的有效性和稳定性，本文首先以一半球为例，对半球在波浪中做强迫垂荡运动的辐射波浪力进行数值计算，并与解析解进行了对比。半球强迫垂荡运动幅值取为口R=0

17、05，运动周期对应的等效波数取为kR=10，半球网格划分节点总数为M一。R=200，时间步长取为At=T30，时间变量t无量纲化处理为tr。网格划分如图2所示，强迫垂荡运动波浪力如图3所示。由图3可知，利用本文所提出的基于加速度势的波浪力数值计算结果与Hulme解析解【l 7】吻合良好，初步验证了本文方法对于无航速情形的适用性。图2半球网格划分示意图图3半球强迫垂荡波浪力时历22算法收敛性检验针对本文所提出的基于加速度势的波浪力数值模拟算法，有必要对其收敛性进行探讨，以寻求一种最为经济、合理的计算参数选择方案。为此，本文以WigleyI船型在静水中做强迫垂荡运动为例，万方数据中 国 造 船 学

18、术论文进行了网格收敛性和时间收敛性检验，WigleyI船型的定义由Jo啪ee18】给出。对于网格数目收敛性检验，计算时间步长取为&=T30。不同节点数目下垂荡运动响应之间的对比如图5所示，其中垂荡运动仉通过波幅h。进行无量纲化，时间变量通过重力加速度g和船体垂线间长厶。进行无量纲化。图中。表示垂向划分节点数目，l表示纵向划分节点数目，典型网格划分示意图如图4所示。由图5可知，考虑计算的经济性，并兼顾计算精度的要求，对于常规船舶纵向划分40个节点、垂向划分8个节点即可。图4 WigleyI船型典型网格划分示意图图5计算网格数目收敛性检验(Fr=0O)对于时间步长收敛性检验，节点划分数目取为D=8

19、、L=40。不同时间步长下所得到垂荡运动响应之间的对比如图6所示。由图6可知，对于常规船舶的水动力计算，时间步长一般取为At=T30即可满足计算精度的需求，并且能够同时兼顾计算效率。图6计算时间步长收敛性检验(Fr=00)万方数据57卷第3期(总第219期) 李志富，等：基于加速度势的船舶波浪力数值模拟23波浪力数值模拟为进一步验证基于高阶面元离散的加速度势数值模型的正确性和有效性，针对无限水深中有航速WigleyI船型的波浪力和运动响应进行了数值计算。该船的几何特征、质量特性以及模型试验结果由Joumee【l 8】给出，试验中仅在垂荡和纵摇两个方向有自由运动。在WigleyI船型的波浪力及运

20、动响应模拟中，航速弗劳德数取为Fr=02，波长取为五k=150，波幅取为危=130cm，浪向角取为顶浪航行，船体纵向网格节点划分为T=40，垂向网格节点划分为No=8，计算时间步长取为At=Z30。其中厶为垂线间长，为航行船舶的遭遇周期。为了与Journee的试验值进行对比，在数值模拟过程中，船体仅在垂荡和纵摇两个模态有位移响应。垂荡和纵摇两个模态的波浪时历分别如图7和图8所示，其中波浪力仅包括水动力。垂荡和纵摇两个模态的位移时历分别如图9和图10所示，图中将数值计算结果与Joumee的模型试验结果进行了对比。由图9图10可知，利用本文所提出的基于高阶面元离散的加速度势波浪力数值求解算法所得到

21、的运动响应与试验值符合良好。但是其中垂荡运动模态平衡位置下移了一段位置，纵摇运动模态的平衡位置上移了一段距离，其原因在于大地坐标系下的边界积分方程在求解速度势的同时考虑了定长兴波流场对于非定常运动的影响。图7航行船舶垂荡模态波浪力(n一02)图8航行船舶纵摇模态波浪力(n=O2)万方数据中 国 造 船 学术论文3 结 论图9航行船舶垂荡运动时历(Fr=02)图10航行船舶纵摇运动时历(n=02)基于时域势流理论，建立了求解流场速度势的边界积分方程。通过引入流体加速度势有效地避免了Bernoulli方程中速度势时间偏导数的有限差分求解所带来的问题，并通过利用高阶曲面元进行边界积分方程的数值离散实

22、现了减少格林函数计算次数的目的。本文通过半球强迫垂荡运动辐射波浪力的计算和WigleyI型船舶水动力的相关计算得出以下主要结论：(1)半球强迫垂荡运动辐射波浪力的数值模拟结果与Hulme提供的解析解吻合良好，初步证明了本文所提出的基于高阶曲面元离散的加速度势的求解方法的正确性。(2)通过对给定时间步长，不同纵向和垂向网格划分所得到的WigleyI船型的运动响应对比，发现对于常规船舶的计算，船舶纵向划分40个节点、垂向划分8个节点即可满足计算收敛性的要求。(3)通过对给定船体网格划分方式，不同时间步长选取所得到的WigleyI船型的运动响应对比发现，对于常规船舶的计算，一般取时间步长At=T30

23、即可满足计算精度的需求，并且能够同时兼顾计算效率。(4)通过数值计算得到的航行船舶的波浪力和运动响应时间历程发现，在大地坐标系下建立的基于加速度势的边界积分方程求解式，能够很好地计及定常兴波流场对船舶非定常运动的影响。此外，需要说明的是对于常规的非外飘船型，可以利用本文基于加速度势的算法直接进行数值模拟；对于常规的弱外飘船型，可以采用修正水线附近单元的处理方法以避免应用本文算法时可能万方数据57卷第3期(总第219期) 李志富，等：基于加速度势的船舶波浪力数值模拟 107会出现的数值不稳定；对于常规的大外飘船型，为了避免求解时可能出现的数值发散问题，可以采用基于加速度势求解的RankineGr

24、een匹配边界元法进行模拟，此种匹配算法将作为本文进一步研究的工作。；，参考文献【l】vAsQUEZ G FONSECA N，SOARES C G Experimental and numerical study ofthe vertical motions ofa bulk carrier anda Ro-Ro ship in extreme waves明Journal ofOcean Engineering and Marine Energy,2015：1-172】LIN Woeimin，YUE DNumerical solutions for large-amplitude ship m

25、otions in the time domainCproceedings of theEighteenth SymposiumonNavalHydrodynamics，F，1991，NationalAcademicPress，Washington，DC3】DATTASEN DA B-spline solver for the forward-speed diffraction problem of a floating body in the time domain明Applied Ocean Research，2006，28(2)：147-160【4】DATTA K RODRIGUES J

26、 M，SOARES C G Study of the motions of fishing vessels by a time domain panel method【J】Ocean Eng；201 l，38(5-6)：782-7925】wU G X，EATOCK T RTransient motion of a floating body in steep water wavesC#Proceedings of the l lthInternational Workshop on Water Waves and Floating Bodies(IWWWFB)，F，199616mJANG Yi

27、fengNonlinear ship motions by a mallkille panel method【D】Massachusetts Institute ofTechnology,19977】ZHOU B Z，WU G XResonance of a tension leg platform exited by third-harmonic force in nonlinear regular wavesJ】Philosophical Transactions of the Royal Society of London A：Mathematical，Physical and Engi

28、neering Sciences，2015，373(2033)：20140105【8】BECK R FLIAPIS STransient motions offloating bodies at zero forward speedJ】Journal ofship research，1987，31(3)：164176【9】NEWMAN J NThe approximation offree-surface Green fimctionsJ】Wave Asymptotics，1992：107-135【10CLI!MENT A HAn ordinary differential equation

29、for the Green function oftime-domain free-surface hydrodynamics叨Joumal 0fEngineering Mathematics，1998，33(2)：201217【ll】U Z FREN H L，TONG X Wet a1A precise computation method oftransient free surface Green functionJOcean Eng，2015，105：318-32612】DAI Y，DUAN W Y Potential flow theory of ship motions in wa

30、vesJ】National Defense Industry Publishing House，Beijing，200813】WEHAUSEN J VLAITONE E V Surface wavesMEncyclopaedia ofPhysics，1960，V01IX【14】KUKKANEN TMATUSIAK JNonlinear hull girder loads ofa RoPax shipJOcean Eng，2014，75：11415】CLEMENT A HRecent developments of computational timedomain hydrodynamics b

31、ased on a differential approach ofthe green functionCProceedings of the Recent computational developments in steady and unsteady navalhydrodynaroics Colloquium,F199816】PRESS H，：TEUKOLSKY A，VETTERLING Tet a1Numerical recipes in C+111e art of computer programming【M】Cambridge University Press New York2

32、00217】HULME AThe wave forces acting on a floating hemisphere undergoing forced periodic oscillationsJ】Journal of FluidMechanics，1982，121：443-463【18】JOURN童E JExperiments and calculations on 4 Wigley hull forms in head waves【R】Delft University of TechnologyReport,1992万方数据中 国 造 船 学术论文Prediction of匀ve F

33、orceBased on Acceleration Potential MethodasedOn Acceleratmn OtentiaLI ZhifU REN Huilong+，SHI Yuyun(College ofShipbuilding Engineering，Harbin Engineering Uniwrsi哪,Harbin 150001，China)AbstractShips motion response to wave load is usually evaluated in the earthfixed coordinate syStem，and timepartial d

34、erivatives of velocity potential should be compmedTo avoid elrors coming from the differentialscheme，acceleration potential is introduced，and the bodys surface condition satisfied by the accelerationpotential is derivedHigher order eightnode panels are adopted for the discretization of the boundary

35、integralequationComputed hydrodynamic coefficients for a hemisphere and Wigley hull agree well withexperimental dataKeywords：earthfixed coordinate；Green function；acceleration potential；wave force；fmite differentialmethod作者简介李志富 男，1990年生，博士研究生。环境载荷与结构强度分析。任麓龙男，1965年生，教授，博士生导师。石玉云女，1990年生，博士研究生。环境载荷与结构强度分析。+通信作者：任慧龙万方数据