MPS-FEM 方法模擬彈性船體在規(guī)則波中的運動1)

黃聰祎 趙偉文 萬德成

(上海交通大學 船舶海洋與建筑工程學院,船海計算水動力學研究中心,上海 200240)

引言

船舶在海上航行時,經常會遭遇波浪的作用.在波浪的作用下,船體可能會發(fā)生較大幅度的運動響應,給船上人員和生活帶來諸多不便.另外,船體的大幅度運動也可能會造成船體表面與水面的劇烈砰擊,作用在船體表面的砰擊壓力可能會破壞船體結構,危害船體和所載人員的安全.因此,準確地預報船體的運動響應,并估計船體受到的砰擊載荷,對于提高船舶的安全系數,提升操作性能都具有重要的意義.

船體在波浪中的運動屬于一典型的流固耦合問題.當船體運動幅度較小時,船體的彈性可以忽略,船體運動可以簡單的視為剛體浮體的運動.但當波高較高,波浪砰擊載荷較大時,船體可能會在波浪的作用下發(fā)生變形,此時船體的彈性則不可忽略.

近年來,國內外諸多學者都在研究浮體在波浪中的運動方面完成了許多工作.實驗方法在處理此類問題時成本較高,實驗難度也較大.在分析結構受力時,只能通過在結構表面附著壓力測點或應變片的方式來實現.這樣不能獲得結構表面連續(xù)的壓力分布,很難準確地捕捉到船體表面的最大壓力,更難以分析流場中的壓力分布情況.隨著近些年高性能計算技術的發(fā)展以及軟硬件性能的提高,數值模擬憑其高效、環(huán)保、適應性強的特點,被越來越多的研究人員選擇.

在求解流固耦合問題時,由于流體和固體的物理差異巨大,因此,往往對流體域和結構域采用不同的方法分別計算,然后通過流固界面上信息的傳遞實現二者的耦合.目前已經有很多研究人員采用耦合的方法對流固耦合問題進行模擬.王加夏等[1]通過計算流體力學分析軟件STAR-CCM+模擬流場,使用有限單元分析軟件Abaqus 模擬結構,通過二者之間的耦合,模擬了彈性船體在波浪中的運動,結果與試驗結果吻合較好.Oberhagemann等[2]使用黏流求解器模擬流場,使用有限元求解器ANSYS 模擬結構,模擬了LNG 船在波浪中的運動,并計算了船體受到的砰擊載荷.Lakshmynarayanana等[3]使用RANS 方法,在StarCCM+平臺上模擬流體,使用FEM方法使用ABAQUS 軟件求解結構,通過二者的耦合模擬了彈性船體在規(guī)則波浪中運動響應.Kim等[4]將BEM 方法與FEM 方法耦合,模擬了大型船舶在迎浪時的變形情況.Malenica等[5]將FEM 方法與格林函數法結合,模擬了集裝箱船在波浪中運動的流固耦合問題.

以上例子證明了在求解流固耦合問題時,耦合類的方法具有較好的表現.本文使用MPS-FEM 耦合方法,流體域的計算采用MPS 方法.MPS 方法作為一種無網格方法,追蹤粒子運動,具有拉格朗日特性,其在處理自由表面變形較大的問題時,具有天然的優(yōu)勢.目前,無網格粒子法也被廣泛地應用到對流固耦合問題的模擬中.Gao等[6]使用SPH 方法在二維數值水池中,模擬了規(guī)則波與結構的相互作用,分析了規(guī)則波與結構底面相互作用附近的壓力分布和速度分布;Omidvar等[7]研究了二維浮體與波浪相互作用問題.也有一些研究人員利用粒子法模擬了結構與規(guī)則波之間的相互作用.Sueyoshi等[8-9]首先模擬了二維的波浪中的浮體運動問題,計算得到的浮體的運動曲線與試驗數據相比吻合較好.隨后又增大規(guī)則波的波幅,模擬船體剖面在規(guī)則波中的劇烈運動.饒成平等[10]使用MPS 耦合方法計算了孤立波對彈性結構的砰擊載荷.Zhang等[11]使用MPSFEM 耦合方法,模擬了孤立波對彈性平板的砰擊過程.Lind等[12]考慮砰擊過程中空氣的影響,采用兩相流模型模擬了波浪對固定平板的砰擊.Sun等[13]使用兩相SPH 模型,模擬了畸形波對固定方柱的砰擊作用.Zhang等[14]使用MPS-FEM 方法,模擬了彈性浮體和彈性系泊平臺在孤立波中的受力及運動.Zhang等[15]用MPS 方法模擬了二維浮體在規(guī)則波中的運動,研究了規(guī)則波頻率對浮體的橫搖響應的影響.Xie等[16]采用MPS-DEM耦合方法,模擬了流體與彈性薄壁結構的相互作用.Zhang等[17]采用MPS-FEM 方法,模擬了三維潰壩流域彈性閘門的相互作用.Sun等[18]采用MPS-DEM耦合方法,模擬了彈性系泊平臺在遭遇孤立波時受到的砰擊壓力及變形情況.Sun等[19-20]使用MPS 方法,模擬了二維和三維浮式平臺在劇烈波浪作用下的變形和運動.

本文使用MPS-FEM 方法模擬彈性船體在規(guī)則波中的運動.在第1 節(jié)中,分別介紹了MPS 方法,FEM 方法的控制方程和求解流程,以及二者之間的耦合方法.在第2 節(jié)中,分別驗證了MPSFEMSJTU 求解器的造波消波模塊,以及模擬彈性體在規(guī)則波中的變形情況的準確性和穩(wěn)定性.在第3 節(jié)中,首先模擬了剛性船體在規(guī)則波中的運動,研究波長對剛性船體運動響應的影響;接著模擬了彈性船體在規(guī)則波中的運動,并與剛性船體的結果對比,分析船舶的彈性對運動的影響.最后,第4 節(jié)總結了本文得到的結論.

1 數值方法

本文采用MPS-FEM 耦合方法,模擬彈性船體在規(guī)則波中的運動.在流體域使用MPS 方法,在結構域使用FEM 方法.在本節(jié)中,將分別介紹MPS 方法和FEM 方法的求解流程,以及二者的耦合方法.

1.1 控制流場的MPS 方法

本文采用一種無網格方法,移動粒子半隱式方法(moving particle semi-implicit method,MPS)來模擬流體的運動.傳統(tǒng)的網格類數值模擬方法將計算域離散為一系列連續(xù)的網格,這種方法處理具有自由面大變形特性的問題時,很容易出現網格畸變的現象,導致計算的失效.而具有拉格朗日特性的無網格方法則不會出現這樣的問題.MPS 方法通過一系列均勻分布的粒子來離散計算域.每個粒子都攜帶著自身的物理信息,包括粒子類型,位置,壓力,速度等.每個粒子之間不存在固定的拓撲關系,而是通過一個權函數—核函數,來傳遞壓力或位移等物理信息.通過追蹤流場內每個粒子的運動來實現對流動的模擬.基于這一特點,無網格方法可以非常容易地追蹤自由表面,這是網格類方法無可取代的優(yōu)勢.因此,無網格方法被廣泛地應用于處理自由表面變形較劇烈的問題.

在MPS 方法中,控制方程為模擬不可壓縮黏性流體的連續(xù)性方程和N-S(Navier-Stokes)方程.拉格朗日形式的控制方程可寫為

式中,ρ為流體密度,P為壓力,V為速度向量,g是重力,v是流體的運動黏性系數.

接下來介紹MPS 方法對計算域的離散方法.MPS 方法在離散計算域時不需要網格,而是將計算域離散為一系列無固定拓撲關系的粒子,粒子的運動由其鄰居粒子的相應信息進行加權平均得到,這一過程通過核函數來實現.本文中使用的核函數如下所示

式中,r表示兩個粒子之間的距離,re表示粒子的影響半徑.式(3)所示的核函數是一種改進了的核函數,它既保持了原有核函數的簡潔形式,又不存在奇點,避免了當兩個粒子距離較近時出現的壓力較大的現象,從而增強了壓力場的穩(wěn)定性.

在MPS 方法中,使用粒子數密度來表示粒子分布的疏密程度,其定義為粒子作用域內所有粒子的核函數之和,定義式如下

對于不可壓縮流體,流體內部的粒子數密度是常數.

對于控制方程中壓力梯度的離散,采用梯度模型來完成.在MPS 方法中,保持壓力場的穩(wěn)定性至關重要.流場壓力的振蕩不僅會帶來計算誤差,而且可能會導致計算的崩潰.因此,選擇準確穩(wěn)定的壓力梯度模型至關重要.近年來,眾多研究者按照不同的思路改進壓力梯度模型[21-23],改進的總的原則就在于要保證粒子間的作用力始終為排斥力,從而保持動量守恒,保證計算的穩(wěn)定性.本文中使用的壓力梯度模型為Tanaka等[23]提出的梯度模型,表達式為

式中,D表示問題的維度,n0表示初始粒子數密度.張雨新[24]設計了一系列數值模擬試驗,驗證了式(5)所示的壓力梯度模型可以較好地保證系統(tǒng)的動量守恒,改善壓力控制方程中的速度散度項通過散度模型壓力離散,散度模型表達式為

在MPS 方法中,粒子的壓力通過求解壓力泊松方程得到.因此,壓力泊松方程模型的選擇也關系到計算的穩(wěn)定性和準確性.許多研究人員,嘗試從不同的角度來對壓力泊松方程進行改善[25-26],其中包括對拉普拉斯算子的改善以及對源項的改善.傳統(tǒng)的MPS 方法中的泊松方程的源項完全由粒子數密度構成,但當流動較為劇烈時,粒子分布不均勻,必然會造成流場壓力計算的不準確,導致流場壓力振蕩.于是又有部分學者提出,利用不可壓流動中速度散度為0的特點,使用速度散度作為泊松方程源項.以速度散度為源項的泊松方程可以得到較為穩(wěn)定光滑的壓力場,但是由于計算速度散度時誤差的積累,可能會導致流體體積的不守恒.為了綜合以上兩種泊松方程源項的優(yōu)點,Tanaka等[23]提出了一個混合源項法,即速度散度和粒子數密度在計算中各占一定的比例,如式(7)所示

上角標k,k+1和*分別表示當前時間步、下一時間步和中間時刻,γ表示粒子數密度源項中所占的比例.Lee等[27]研究得出γ取值為[0.01,0.05]區(qū)間時,流場壓力的穩(wěn)定性和流體體積的守恒性較好.

在求解壓力泊松方程時,將自由表面粒子壓力為0作為邊界條件.因此,準確地捕捉自由表面粒子對于穩(wěn)定地求解壓力泊松方程非常重要.在MPS 方法中,使用粒子數密度來判斷粒子分布的疏密程度.因此可以利用自由表面粒子周圍的鄰居較少的特征,使用粒子數密度來判斷自由表面粒子.在本文中,當 〈ni〉 ＜0.8n0時,n0為初始流體內部粒子數密度,認為該粒子位于自由表面,當 〈ni〉 ＞0.97n0,粒子則被判斷為流體內部粒子.若僅以這種方式來判斷,則在流動較為劇烈,粒子分布不均勻時,非常容易造成自由表面粒子的誤判,導致在流體內部出現壓力為0的位置,導致流場壓力的振蕩.因此,Zhang等[28]引入了一相對位置矢量函數F,利用自由表面處粒子兩側的粒子分布不均勻的特征來判斷粒子是否位于自由表面上.相對位置矢量函數的定義如下

當〈 |F|〉i>0.9|F|0時,認為粒子位于自由表面.其中,|F|0表示初始時刻自由面粒子的 |F| 值.張雨新[24]通過數值模擬,驗證了這種自由表面粒子判斷方法可以有效地避免當流動劇烈時,將流體內部的粒子誤判為自由表面粒子,從而提高了流場壓力的穩(wěn)定性.在處理固定邊界時,為了保證邊界附近粒子支持域的完整性,將邊界設成一層第一類邊界粒子和兩層第二類邊界粒子,如圖1 所示.其中,第一類邊界粒子的壓力通過和流體粒子一起參與壓力泊松方程的求解來得到.第二類邊界粒子的壓力通過第一類邊界粒子和流體粒子的壓力插值得到.這種多層邊界粒子的方式既可以避免粒子的穿透,又可以保證邊界附近流體粒子支持域的完整.

圖1 固壁邊界粒子示意圖Fig.1 Schematic diagram of wall boundary particles

1.2 求解結構的FEM 方法

在這一小節(jié)中,介紹用于計算結構受力和變形的FEM 方法.FEM 方法的控制方程為動力平衡方程,即建立節(jié)點單元的平衡條件.控制方程可以寫為

式中,M表示質量矩陣,K表示剛度矩陣,C表示阻尼矩陣.其中,為了簡化計算,阻尼矩陣可以表示為質量矩陣M和剛度矩陣K的線性疊加,可以寫為

式中,α1和α2分別表示瑞利阻尼系數.F(t)表示節(jié)點處隨時間變化的外力矢量.本文采用Newmark-β[29]方法求解結構動力平衡方程,通過對速度和位移進行泰勒展開,可以得到在下一時間步t=t+Δt的節(jié)點位移,可以寫為

式中,β和γ是保證計算穩(wěn)定性的控制參數,需要滿足 γ ≥0.5,β ≥0.25(0.5+γ)2條件,在本文模擬中,β=0.25,γ=0.5.將式(11)和式(12)代入控制方程中,就可以求解得到下一時間步的結構節(jié)點位移.求解流程表達式如式(13)～式(16)所示

1.3 MPS 與FEM 的耦合策略

在本文中,MPS 與FEM 求解模塊之間的信息傳遞通過串行交錯算法(conventional serial staggered,CSS)的分區(qū)耦合策略來實現.結構域和流體域分別在各自的區(qū)域內獨立求解,之后在流固交界面處通過插值完成信息的傳遞.當結構單元的離散形式不同時,流固界面不同,插值的方式也不同.在本文中,采用核函數插值技術,利用粒子法的核函數具有的插值功能,實現了流固異構界面上的位移和載荷信息的傳遞.圖2 中分別展示了流固界面之間的位移和載荷的插值技術.

如圖2(a),流體粒子m的位移為

圖2 流固界面插值技術Fig.2 The schematic diagram of coupling fluid-structure interface

式中,i表示以rei為插值半徑的流體粒子m的插值域內的計算結構節(jié)點,δi表示節(jié)點i的位移.節(jié)點i的載荷為

式中,i表示計算結構節(jié)點插值域內的流體粒子,Pi為流體粒子i的流體壓力,l0為流體域初始粒子間距.

2 數值方法驗證

在這一小節(jié)中,將對本文中使用的數值方法以及采用的求解器MPSFEM-SJTU 進行驗證.驗證分為造波模塊驗證和彈性體變形驗證兩個部分.

2.1 數值造波模塊驗證

首先在三維數值水池中生成三維規(guī)則波.數值波浪水池模型如圖3 所示,水池長L=12 m,寬D=1.5 m,水深h=1.8 m.水池左側設置了造波板,通過推板造波的方式生成波浪.造波板的位移可基于線性造波理論推導得出,表達式為

圖3 數值波浪水池示意圖Fig.3 Schematic diagram of numerical wave pool

式中,H表示波高,k表示波數.在水池的右側設置了消波區(qū),消波區(qū)長l=3 m.消波是通過在方程中增加源項來實現的,源項可以寫為

式中,xs是消波區(qū)的起始位置,Ls是消波區(qū)的長度,αs是一個無因次化的人工黏性系數,用于控制消波強度;u表示流體速度矢量.生成目標規(guī)則波的波長為λ=3m,通過Airy波的色散關系ω=以及T=2π/ω,求得規(guī)則波周期為T≈1.407 8 s.波高H=0.3 m,波陡 δ=H/λ=0.1.

取三種初始粒子間距對造波模塊進行收斂性驗證,分別為Δx=0.04 m,Δx=0.03 m和Δx=0.025 m,分別模擬規(guī)則波的生成和傳播過程.在水池底部分別布置了波高測點,波高測點的橫坐標分別為x=4.5 m,x=6.0m,x=7.5 m和x=11.0m.其中,x=11.0m 的測點位于消波區(qū)內,用以驗證消波區(qū)的消波效果.四個測點處的波高時歷曲線如圖4 所示.

從圖4 中可以看出,在前面三個測波點處,波高的時歷曲線均與解析解吻合較好.當Δx=0.04 m 時,模擬得到的規(guī)則波的波高不是很穩(wěn)定,與解析解相差較大.并且隨著波浪的傳播,在距離造波板較遠的位置,波高略低于理論值,這可能是波浪在船舶運動過程中會出現能量的衰減.但隨著初始粒子間距減小,時歷曲線逐漸貼近解析解,證明了求解器的造波模塊具有較好的收斂性.并且從圖4 中也可以看出,對于粒子間距為Δx=0.03 m和Δx=0.025 m的兩個算例,已經不存在明顯的波高衰減的現象.由此可以說明,減小初始粒子間距可以較好地避免出現波高衰減的現象.對于位于消波區(qū)內的測點,從圖3(d)中可以看出,波高測點的時歷曲線非常平穩(wěn),證明了求解器消波模塊的有效性.

圖4 四個測點處的波高時歷曲線Fig.4 Time history curves of wave height at four measuring points

圖5為粒子間距為Δx=0.03 m 時,在一個周期內,波浪場的波面分布.從圖5 中可以看出,生成的規(guī)則波波面具有較好的正弦特性,波面連續(xù),無明顯的非線性特征.沿著波浪的傳播方向不存在波高的衰減現象.因此可以證明本文使用求解器可以較好地實現造波和消波功能.

圖5 一個周期內規(guī)則波的波面分布情況Fig.5 Wave surface of regular wave in a period

2.2 彈性結構變形驗證

在這一節(jié)中,模擬二維圓柱繞流后彈性懸臂梁的運動,驗證求解器彈性結構在流體中的變形的準確性.

本文的計算模型按照Turek等[30]的模型實驗尺度建立,如圖6 所示.通道內有一個固定圓柱,圓柱后連接著一個彈性板.其中,通道高H2=0.41 m,通道長L2=2.5 m,圓柱直徑為D=0.1 m.以通道左下角為坐標原點建立坐標系,圓柱圓心坐標為(0.2,0.2),圓柱后的懸臂梁厚H1=0.02 m,長L1=0.35 m,懸臂梁使用梁單元來離散,懸臂梁端點的坐標為(0.6,0.2),彈性懸臂梁的彈性模量為Es=5.6×106Pa.通道內來流平均速度為U=1 m/s,方向沿著通道方向,流體的黏度為μ=1.0×10-3m2/s.粒子間距取為Δx=0.002 5 m,D/Δx=40.

圖6 計算模型示意圖Fig.6 Schematic diagram of the calculation model

來流流經靜止的圓柱后,會在圓柱后產生均勻規(guī)律的渦脫落現象.圓柱后的彈性懸臂梁受到渦脫落的影響產生受迫振動,其振動頻率與漩渦脫落頻率相同,周期約為T=0.8 s.圖7 展示了在半個周期內,彈性懸臂梁的變形情況與其他數值方法得到的模擬結果的對比.

圖7 t=t0,t=t0+T/8,t=t0+T/4,t=t0+3T/8 時刻彈性懸臂梁變形情況Fig.7 Deformation of elastic beam at t=t0,t=t0+T/8,t=t0+T/4,t=t0+3T/8

從圖7 中可以看出,圓柱后的彈性懸臂梁,表現出了較好的彈性特征,變形符合物理規(guī)律,并與基于SPH[31]方法模擬得到的結果吻合較好.為進一步驗證對于結構變形求解的準確性,記錄了兩個周期內懸臂梁尾端的位移隨時間的變化,并與文獻[30-31]得到的計算結果進行對比,如圖8 所示.從圖8 可以看出,MPSFEM-SJTU 求解器模擬得到的懸臂梁端點位移時歷曲線,與文獻[30-31]得到的數值計算結果在t=1.0T之前稍有差異,可能的原因是本文中選取的初始粒子間距較大,為D/Δx=40,而文獻[30]選用的初始粒子間距為D/Δx=100.但時歷曲線的振幅較為接近,可以證明本文使用的求解器在模擬彈性結構在流體中的變形時,具有較好的準確性和穩(wěn)定性.

圖8 懸臂梁端點位移的時歷曲線Fig.8 The time history curves of the displacement of the cantilever beam endpoint

3 數值模擬

在本小節(jié)中,模擬三維船體在規(guī)則波中的運動.首先模擬了剛性船體在規(guī)則波中的運動,研究了規(guī)則波的波長對船體運動的影響.接著模擬了彈性船體在規(guī)則波中的運動,并與剛性船體的運動進行對比,研究了船體的彈性對船體運動的影響.

本文中選用的船模為船長3.2 米的kvlcc2 船,縮尺比為1:100.kvlcc2 船的實船尺度和模型尺度如表1 所示.數值水池的尺度與2.1 節(jié)中的造波水池一致,船艏與造波板的距離為3 m.

表1 kvlcc2 船實船和模型船參數Table 1 Characteristic parameters of the kvlcc2 ship

在對船體模型進行離散時,可以離散為全船有限元模型,也可以離散為非均勻梁模型.前者需要大量的建模工作和計算資源,而采用變截面梁模型可以大大簡化計算,提高工作效率.在本文中選擇非均勻梁模型,將船體沿船長方向離散為一系列變截面梁單元,船體模型與梁單元的等效關系如圖8 所示.圖9(a)展示了船體梁模型分組插值的方法,即將位于同一縱剖面處的船體粒子分為一組,這樣一組粒子被等效為梁單元上的一個節(jié)點.從圖9(b)可以看出,通過映射關系,將船體梁單元上的節(jié)點的位移和變形信息傳遞到船體邊界粒子上,船體邊界粒子的載荷信息也可以通過這樣的映射關系傳遞到梁單元的節(jié)點上.

圖9 船體粒子與梁模型的等效關系Fig.9 Equivalent relationship between hull particles and beam model nodes

船體模型的初始粒子間距為dp=0.025 m,根據船體模型的特征參數,將船體密度設為ρs=945 kg/m3.非均勻梁的橫截面積、質量分布、橫截面慣性矩在各個截面各不相同,在本文中的計算式為

式中,A,M,I分別表示橫截面積、截面質量和截面慣性矩.下標i,m和S分別表示船體邊界粒子、梁單元的節(jié)點和船體結構.圖10為統(tǒng)計得到的船體梁截面質量和截面慣性矩沿船長方向的分布.將圖10(a)曲線進行積分,可以得到統(tǒng)計的船體梁模型總質量為311.6 kg,與表1 中數據一致.由此可以說明本文中采用的離散方式的合理性.

圖10 非均勻船體梁模型特征參數沿船長的變化Fig.10 Variation of characteristics along the length of a nonuniform hull beam model

3.1 模擬剛性船體在三維規(guī)則波中的運動

首先對船體在規(guī)則波中的運動響應進行收斂性驗證,選取三組粒子模型,粒子間距分別為:Δx=0.03 m,Δx=0.035 m和Δx=0.04 m,三組計算模型對應的粒子總數分別為101 萬和64 萬和43 萬.圖11展示了初始粒子間距Δx=0.03 m 的模型,在一個波浪周期內船體的運動形態(tài)和自由面速度分布.

從圖11 中可以看出,MPSFEM-SJTU 求解器可以較好地模擬三維船體在規(guī)則波中的運動情況.模擬得到的船舶運動形態(tài)以及自由表面粒子速度分布都符合物理規(guī)律.在規(guī)則波的作用下,船舶會發(fā)生周期性的運動,在本文中只考慮垂蕩和縱搖的耦合運動.

圖11 三維船體在波浪中的運動Fig.11 The ship motion in regular wave

圖12 是10～ 20s 內,不同的粒子間距下船體垂蕩和縱搖的時歷曲線.從圖中可以看出,在三組粒子間距下,船體的運動響應都表現出了較好的周期性,符合規(guī)則波激勵下的響應特征.但粒子間距對船體運動響應影響較大.當Δx=0.04 m 時,船體的垂蕩和縱搖運動都較小,這是因為當粒子間距較大時,粒子在運動過程中會出現較為明顯的能量衰減.但當初始粒子間距減小時,衰減程度減小,驗證了本文的算法具有較好的收斂性.關于船體運動響應時歷曲線數據的分析如表2 所示.可以看出,粒子間距增大,船體運動響應幅值也會大幅度減小.從圖4 中的分析可知,當粒子間距較大時,規(guī)則波的波高也會低于理論值.因此,此處的船體運動響應曲線振幅的差別不僅由粒子間距造成,規(guī)則波波高的不同也會影響船體的響應幅值.從圖12(b)和表2 中可以看出,粒子間距對于模擬船體的縱搖運動影響較小.與垂蕩響應曲線相比,三種粒子間距下船體縱搖的時歷曲線相差較小.

表2 船舶運動響應參數Table 2 Ship motion response parameters

圖12 船體在規(guī)則波中運動的時歷曲線Fig.12 The time histories curves of the ship motion

接下來分析規(guī)則波的波長對船體運動響應的影響.在這一小節(jié)均采用初始粒子間距Δx=0.03 m 的模型來進行分析和計算,船型和水池尺度與上一小節(jié)相同.取波長分別為λ=2 m,λ=3 m,λ=3.2 m,λ=3.5 m,λ=4 m進行實驗.已知船長L=3.2 m,即3 組實驗的波長分別為λ=0.625L,λ=0.9375L,λ=L,λ=1.09375L和λ=1.25L,規(guī)則波的波長均在船長附近.在本文中計算得到的船體的垂蕩和縱搖時歷曲線如圖13 所示.

由于本文波浪基于線性波浪理論,規(guī)則波的波長與頻率成反比,所以隨著波長的增大,規(guī)則波的頻率降低.由于受迫振動的響應頻率與激勵頻率相同,所以在相同波長規(guī)則波的激勵下,船體的垂蕩和縱搖響應頻率相同.隨著規(guī)則波波長增大,船體運動響應的頻率降低.在不同的規(guī)則波激勵頻率下,船體運動響應的幅值也不同,統(tǒng)計了圖13 中幾個周期內的船體運動響應平均幅值,繪制了船體運動響應隨波長的變化曲線,如圖14 所示.

圖13 不同波長下船體運動響應的時歷曲線Fig.13 Ship motion response at different wavelengths

圖14 不同波長下的船體運動響應平均幅值Fig.14 Motion response average amplitude at different wavelengths

從圖14 中可以看出,當規(guī)則波的波長由2 m 增大到3 m 時,縱搖平均幅值增大,而垂蕩平均幅值減小,當波長從 λ=3 m 繼續(xù)增大時,垂蕩和縱搖曲線的變化趨勢大致相同.當波長 λ=3.2 m=L時,船體的垂蕩運動和縱搖運動的幅值都最小.隨著波長增大,船體的運動幅值逐漸增大.當波長 λ=4 m=1.25L時,縱搖和垂蕩的幅度同時達到最大值,并且此時的運動響應幅值遠大于其他波長下的平均響應幅值.由以上分析可以認為,對于本文中所研究的波長范圍內(0.625L～ 1.25L),當規(guī)則波的波長大于船長并逐漸增大時,船體的運動響應也逐漸增大.當規(guī)則波的波長小于船長時,垂蕩和縱搖曲線出現極值時對應的規(guī)則波波長可能不相同.

3.2 模擬彈性船體在三維規(guī)則波中的運動

在這一小節(jié)中,模擬彈性船體在規(guī)則波中的運動,水池尺度、船體模型以及造波模式等均與3.1節(jié)相同.圖15 展示了彈性船體一個周期內的運動情況.

從圖15 中可以看出,彈性船體在規(guī)則波的作用下發(fā)生了周期性的縱搖運動,但由于彈性船體的彈性模量較大,剛性較強,因此難以觀察到明顯的變形.圖16 展示了在圖15(c)情況下,即當波峰位于船舯附近時,船體表面的壓力分布情況.

圖15 彈性船體在規(guī)則波中的運動Fig.15 Motion of elastic hull in regular waves

從圖16 中可以看出,與剛性船體相比,彈性船體在船舯部的壓力分布顏色更深,壓力更大,彈性船體的壓力大于剛性船體.彈性船體表面最大壓力為2393.52 Pa,剛性船體底部最大壓力為2202.61 Pa,約為彈性船體最大壓力的92%.并且較大的壓力主要分布在船舯附近.在這一時刻下,造成這一現象的主要原因有以下兩個,第一是波峰位于船舯附近,因此處水深較大,受到的水壓力也較大.第二是彈性船體在這種波浪條件下會發(fā)生中拱,即使在船體彈性模量較大的情況下.船體梁的彎曲也會增大船底粒子受到的壓力.

圖16 彈性和剛性船體的表面壓力分布Fig.16 Pressure distribution of elastic and rigid hull

為研究彈性船體和剛性船體在規(guī)則波中的運動響應的區(qū)別,分別繪制剛性船體和彈性船體的垂蕩和縱搖的時歷曲線,如圖17 所示.

從圖17 中可以看出,剛性船體的垂蕩和縱搖運動幅值都大于彈性船體.8 s 之后,船體的垂蕩曲線的幅值還存在著較小的波動,但船體縱搖的響應曲線幅值基本穩(wěn)定.統(tǒng)計8 s 之后圖17 所示的9 個周期內的船體運動響應的平均幅值,數據如表3 所示.

圖17 彈性和剛性船體的運動響應時歷曲線Fig.17 Time history curves of motion response of elastic and rigid hull

從表3 中的數據可以看出,剛性船體的縱搖平均幅值約為5.676°,彈性船體的縱搖幅值約為4.253°,相比剛性船體減小了25.1%.而對于垂蕩響應,9 個周期內剛性船體的垂蕩平均幅值為0.014 8 m,彈性船體的垂蕩平均幅值約為0.013 1 m,相比剛性船體減小了11.5%.由此可見剛性船體在規(guī)則波中的運動比彈性船體的運動更劇烈,不考慮船體的彈性時,計算得到的船體運動響應會略大于實際的船體運動響應.

表3 彈性剛性船體運動響應參數對比Table 3 Comparison of motion response parameters of elastic and rigid hulls

接下來沿船長方向,分別取了5 個截面,分別位于船艏,距船艏Lpp/3 處,船舯,距離船艏2Lpp/3處,以及船尾.五個截面受到的壓力隨時間的變化曲線如圖18 所示.

圖18 剛性和彈性船體不同截面處的壓力時歷曲線Fig.18 Time history curves of pressure at different sections of rigid and elastic hull

從圖18(a)中可以看出,在船艏處,彈性船體和剛性船體的受力情況有明顯的不同.雖然二者的壓力的時歷曲線都具有明顯的周期性和脈沖特性.但對于剛性船體,壓力大多數情況下都表現為正值,其主要成分應該是浮力.而對于彈性船體,壓力曲線在零附近有規(guī)律地波動,正值和負值同時存在.這種情況下,壓力的組成成分主要來自于浮力和甲板上浪引起的砰擊壓力.從圖15 中可以看出,剛性船體的縱搖幅度較大,在發(fā)生甲板上浪時,船艏已被完全淹沒.因此,剛性船體的船艏處難以捕捉到由甲板上浪帶來的砰擊壓力,這部分砰擊壓力可能會作用到船艏之后的某一截面上.從圖18(b)～ 圖18(d)可以看出,在1/3 船長、船中和2/3 船長這三個截面處,砰擊載荷同樣呈周期性變化,并且剛性船體與彈性船體的載荷曲線較為接近.對位于船尾截面受到的壓力,剛性船體受到的壓力幅值大于彈性船體.并且彈性船體和剛性船體都存在著一部分壓力為零的時刻,這是因為此時船體艏傾,船尾被抬起,離開了水面,因此此時壓力為零.仔細觀察還可以發(fā)現,彈性船體壓力為零的時間要短于剛性船體壓力為零的時間,這同樣是因為彈性船體的運動幅度更小.

4 結論

在本文中,采用MPS-FEM 耦合方法,使用自主開發(fā)的求解器MPSFEM-SJTU,模擬了剛性和彈性船體在規(guī)則波中的運動.首先在三維數值水池中進行空場造波,驗證了求解器的造波模塊和消波模塊.接著通過模擬二維圓柱繞流后彈性懸梁的振動,驗證求解器模擬彈性流固耦合問題的準確性.驗證了求解器后,本文模擬了彈性船體和剛性船體在規(guī)則波中的運動,得到了以下結論.

(1) 船體在規(guī)則波中的縱搖運動比垂蕩運動更穩(wěn)定.當剛性或彈性船體在規(guī)則波中的運動響應發(fā)展到較為穩(wěn)定的階段時,船體縱搖運動的幅度基本保持不變,但船體垂蕩運動的幅值還有輕微的波動.

(2) 彈性船體的運動幅度比剛性船體小,從模擬得到的彈性和剛性船體的運動響應時歷曲線可以看出,無論是垂蕩運動還是縱搖運動,剛性船體的運動幅度都較大.不考慮船體彈性時計算得到的運動響應可能偏大.

(3) 由于彈性船體會發(fā)生彎曲,因此在船體底部,彈性船體船舯附近的壓力大于剛性船體.而由于剛性船體的運動幅度大于彈性船體,所以剛性船體的船艏和船尾處受到的砰擊壓力大于彈性船體.