孫雷，程聰，鄭博學(xué)，劉昌鳳

1 大連理工大學(xué)船舶工程學(xué)院，遼寧大連116024

2 高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心，上海200240

3 大連海洋大學(xué)海洋與土木工程學(xué)院，遼寧大連116023

0 引 言

在船舶流體力學(xué)中，液艙晃蕩為典型的周期性水動力學(xué)問題［1］。海上航行的液化天然氣（liquefied natural gas，LNG）船在風(fēng)、浪、流的聯(lián)合作用下，其液艙內(nèi)會產(chǎn)生巨大的砰擊壓力，這將對艙壁造成直接的損壞［2］。研究如何準(zhǔn)確預(yù)報艙壁所受的砰擊載荷，對于船舶的結(jié)構(gòu)設(shè)計具有重要意義。

由于光滑粒子流體動力學(xué)（smoothed particle hydrodynamics，SPH）方法在模擬自由表面的流動、翻卷和破碎時無需特別處理［3-5］，對交界面捕捉優(yōu)勢明顯［6-7］，因此十分適合求解液艙晃蕩問題。然而，該方法仍存在壓力振蕩的問題，威脅到流場的穩(wěn)定性。雖然利用不可壓縮求解器［4］或?qū)撕瘮?shù)和核梯度進(jìn)行修正［8］可以在一定程度上緩解壓力振蕩，但仍未從根本上解決該問題。

國內(nèi)外學(xué)者針對該問題進(jìn)行了深入研究，并提出了2 種可行的方法：密度正則化和數(shù)值耗散技術(shù)。

第1 種方法是對粒子進(jìn)行密度正則化處理。Randles 和Libersky［9］提出了一種零階精度的密度正則化（Shepard 密度正則化）算法。Shao 等［8］基于該模型改進(jìn)了離散方法，在模擬二維矩形液艙晃蕩問題中獲得了穩(wěn)定的流場壓力；Chen 等［6］在該方法的基礎(chǔ)上改進(jìn)，提出了適用于大密度比兩相流的密度正則化算法，模擬了潰壩和瑞麗泰勒不穩(wěn)定性等問題，改善了壓力結(jié)果的穩(wěn)定性。Dilts［10-11］提出了一種一階精度的密度正則化（MLS密度正則化）算法。Colagrossi 和Landrini［12］利用這種方法緩解了流場的壓力噪聲，穩(wěn)定了大密度比情況下的界面流動問題；Chen 等［4］將MLS 密度正則化算法應(yīng)用到二維液艙晃蕩問題中，提高了流場壓力的計算精度。雖然MLS 密度正則化算法能夠一定程度上抑制壓力振蕩，但是在長時期的模擬過程中，并不能夠保持系統(tǒng)內(nèi)粒子體積的守恒［13］。因此，Shepard 密度正則化算法因其簡捷性和經(jīng)濟(jì)性，被廣泛使用。

第2 種方法是應(yīng)用數(shù)值耗散技術(shù)來抑制壓力振蕩。Monaghan［14］考慮速度的耗散，在動量方程中引入了人工黏性項，人工黏性能夠有效減小壓力噪聲，但是該方法是在整個流場施加速度耗散的，會引入極大的數(shù)值耗散。為了減小密度的非物理性振蕩，Molteni 和Colagrossi［13］通過在連續(xù)性方程中引入密度耗散項建立了δ-SPH 模型，對二維弱可壓縮自由表面流動問題進(jìn)行了模擬，獲得了光滑穩(wěn)定的壓力數(shù)值結(jié)果。Antuono 等［15］在Molteni 和Colagrossi 工作的基礎(chǔ)上改進(jìn)了密度耗散項，解決了流體體積不守恒和自由表面不收斂的問題，有效抑制了流場壓力的振蕩，使之更適用于靜力學(xué)問題。Antuono 等［16］對δ-SPH 模型的能量平衡進(jìn)行了深入分析，發(fā)現(xiàn)該模型的耗散機(jī)制比傳統(tǒng)SPH 模型要好，因此δ-SPH 模型也獲得了廣泛的應(yīng)用。

在計算一些船舶水動力學(xué)問題時，選取合適有效且耗散性比較小的方法來抑制流場壓力的振蕩尤為重要。而目前就此類問題進(jìn)行比較和探討的資料甚少。Calogrossi 等［13，15］在人工黏性修正的基礎(chǔ)上加入密度耗散項綜合比較，證實了密度耗散項加入后能獲得更光滑的壓力場，但是一些低沖擊性的水動力學(xué)問題沒必要施加人工黏性，因為這會帶來不必要的人為數(shù)值耗散。另外，人工黏性是一種速度耗散算法，而密度耗散項是一種密度耗散算法，為了體現(xiàn)實際的物理意義，本文將在僅考慮物理黏性、不加入人工黏性的情況下，比較2 種密度耗散算法對壓力振蕩的抑制效果。在實際工程領(lǐng)域，液艙晃蕩問題的運(yùn)動形式復(fù)雜，本文將考慮運(yùn)動速度不同情況下2 種算法的優(yōu)劣，以便為耗散項在復(fù)雜工況下的選取提供指導(dǎo)。此外，先前討論晃蕩問題的流動形態(tài)時都是直觀比較波面形狀，關(guān)于定量的波高測量算法沒有太多參考信息?；谝陨显?，本文將比較2 種密度耗散算法（Shepard 密度正則化算法和δ-SPH 密度耗散項算法）對壓力振蕩的抑制效果。

1 控制方程

1.1 基本控制方程

假設(shè)流體不可壓縮，考慮絕熱條件，流體運(yùn)動應(yīng)遵循質(zhì)量守恒定律和動量守恒定律。在拉格朗日坐標(biāo)系下，控制方程被簡化為如下形式的連續(xù)性方程和動量守恒方程：

式中：ρ 為流體密度，kg/m3；v 為速度矢量，m/s；P 為壓強(qiáng)，Pa ；F 為體積力矢量，m/s2；D*為黏性耗散項，m/s2。

1.2 狀態(tài)方程

SPH 方法在求解問題時常把不可壓縮流體假設(shè)為弱可壓縮［5，17］。這樣，求解壓力時可避免求解關(guān)于壓力的泊松方程，壓力可以直接由狀態(tài)方程顯式求解得到。在本文中，狀態(tài)方程采用如下形式［6］：

式中：cs為數(shù)值聲速，m/s ；ρ0為流體的參考密度，kg/m3；P0為背景壓力，Pa。為了控制流體的弱可壓縮性，需保證流體的密度變化范圍在1%以內(nèi)，數(shù)值聲速通常取為流場內(nèi)最大速度vmax的10倍，即滿足cs＞10vmax。

1.3 2 種SPH 耗散算法的密度求解方程

1.3.1 Shepard 密度正則化

在傳統(tǒng)SPH 方法中，臨近邊界的流體粒子支持域內(nèi)沒有足夠的臨近粒子來進(jìn)行插值，導(dǎo)致計算邊界時得到的物理量出現(xiàn)了振蕩。在實際計算過程中，會出現(xiàn)粒子分布不均的情況，這種非物理現(xiàn)象會嚴(yán)重影響物理量計算的精度。在計算流場壓力時，為了恢復(fù)支持域的緊支性，保持SPH 粒子近似的一致性，一種常用的辦法就是采用Shepard密度正則化［8］，對計算的密度進(jìn)行重整：

1.3.2 δ-SPH 方法的耗散項

對于弱可壓縮SPH 模型，為了過濾高頻壓力噪聲，抑制壓力振蕩，Molteni 和Colagrossi［13］基于數(shù)值耗散的原理，在連續(xù)性方程中加入了密度耗散項，其表達(dá)式為：

2 SPH 數(shù)值計算模型

2.1 SPH 基本原理

在SPH 方法中，計算空間被離散成一系列任意分布的粒子，粒子攜帶流場計算的基本物理信息，包括位置、速度、密度、壓力、質(zhì)量等，這些粒子之間通過光滑核函數(shù)建立聯(lián)系并相互作用。SPH方法的近似過程包括兩部分［3］：核近似和粒子近似。核近似是指流場內(nèi)某粒子點的場函數(shù)及其導(dǎo)數(shù)在該粒子支持域內(nèi)通過光滑核函數(shù)來近似積分的過程。場函數(shù)及其導(dǎo)數(shù)的核近似形式如下：

式中：r 為某 點的位置矢量，r?=r-r′ ；Ω 為支持域。

為了對控制方程進(jìn)行離散，將上述積分形式的表達(dá)式轉(zhuǎn)換成對支持域內(nèi)粒子的求和形式，即粒子近似。場函數(shù)及其導(dǎo)數(shù)的核近似表達(dá)式的離散求和形式為：

式中，q= ||ri-rj/h。

2.2 物理黏性

本文僅考慮動量方程的耗散項為物理黏性項，物理黏性的表達(dá)式采用Lo 和Shao［18］給出的形式，其離散形式為：

式中：υ 為流體運(yùn)動黏性系數(shù)，m2/s；rij=ri-rj；μi=ρiυ，μj=ρjυ，μi，μj為i，j 粒子的動力黏性系數(shù)，Pa ?s；η=0.1h，為防止出現(xiàn)奇點的系數(shù)。根據(jù)以上描述，控制方程可以表達(dá)成如下離散形式：

2.3 邊界條件

根據(jù)邊界粒子的布置及其場變量的離散形式，邊界條件可以分為斥力邊界條件、動力邊界條件和耦合動力學(xué)邊界條件等。斥力邊界條件通過改變斥力系數(shù)來人工調(diào)節(jié)斥力的作用范圍，但斥力粒子只負(fù)責(zé)提供斥力。另外，斥力粒子通常設(shè)置為一層，這會導(dǎo)致邊界粒子支持域的截斷，使邊界處出現(xiàn)數(shù)值振蕩。動力邊界條件在邊界外設(shè)置若干層虛粒子，虛粒子通過參與控制方程的運(yùn)算進(jìn)而提供物理斥力，但容易出現(xiàn)邊界粒子的非物理性穿透。耦合動力學(xué)邊界由Liu 等［19］提出，Chen 等［6］在此基礎(chǔ)上進(jìn)行了改進(jìn)，對二維情況下的大密度比多相流動問題進(jìn)行了模擬。本文采用Chen等［6］改進(jìn)的耦合動力學(xué)邊界條件，如圖1所示。

圖1 耦合動力學(xué)邊界條件Fig.1 Schematic diagram of coupled dynamic boundary condition

如圖，在邊界上布置一層斥力粒子，斥力粒子的屬性設(shè)為流體的屬性。另外，根據(jù)選取的核函數(shù)確定邊界以外所布置虛粒子的層數(shù)，斥力粒子的物理信息和虛粒子的物理信息都根據(jù)其所在支持域的周圍粒子插值得到。當(dāng)流體靠近邊界區(qū)域時，斥力粒子對流體施加的斥力為：

式中：Rij為流體粒子所受到的邊界對它的排斥力，N；mfluid為流體粒子的質(zhì)量，kg。

2.4 時間積分

本文的時間積分方案采用預(yù)測-校正法，該算法共分為2 步。

第1 步：預(yù)測步。

第2 步：校正步。

一個時間步以后，粒子最終的物理量為

由于該積分方法是一種顯式積分法，為使積分結(jié)果穩(wěn)定，計算時間步長的選取需要滿足柯朗-弗里德里希斯-列維（Courant-Friedrichs-Lewy，CFL）約束條件［3］，并滿足一定的準(zhǔn)則［17，20］，即

式中：Δt1為滿足聲速與光滑長度的時間步長；Δt2為滿足外力條件的時間步長；Δt3為滿足黏性耗散的時間步長；f為粒子加速度，N/m2。

3 數(shù)值算例及結(jié)果討論

針對液艙靜止和晃蕩情況，分別采用2 種SPH 耗散模型對其進(jìn)行數(shù)值模擬。本文的計算模型基于Li 等［21］的實驗，將三維液艙簡化為二維方箱，如圖2 所示。水箱的長度L=1 m，高度H=1 m，水深Dw=0.3 m。不考慮空氣的影響，水箱內(nèi)的介質(zhì)為水，水的密度取為ρ=1 000 kg/m3，水的運(yùn)動粘度取為υ=1.0×10-6m2/s，背景壓力取為P0=0 Pa。關(guān)于該模型的驗證參見文獻(xiàn)［22］，本文不再贅述。

圖2 數(shù)值模型示意圖Fig.2 Schematic diagram of numerical model

3.1 靜水情況

在靜水算例中，在距離水箱底部0.2 m 的左側(cè)壁面上設(shè)置一個壓力監(jiān)測點，用來實時監(jiān)測該處的靜水壓力。粒子初始間距為Δx0=0.005 m ，計算步長根據(jù)CFL 條件取為Δt=5×10-5s ，模擬的總物理時間為55 s，除耗散項不同以外，2 種SPH算法的其他數(shù)值模擬參數(shù)完全一致。

對于傳統(tǒng)SPH 方法，流場內(nèi)部粒子的密度由密度正則化最終得到，對流場內(nèi)部的粒子每隔30個時間步長執(zhí)行一次密度正則化。在δ-SPH 算法中，流場內(nèi)部粒子的密度由引入了密度耗散項的連續(xù)性方程得到，柯西數(shù)δ取為0.05。將監(jiān)測點處的壓力減去該點處的初始壓力P*，并做無量綱處理得到2 種SPH 算法的壓力耗散率，用Rdf表示，即Rdf=(P-P*)/P*。圖3 所示為傳統(tǒng)SPH 方法（即Shepard 密度正則化算法）與耗散系數(shù)為0.05 的δ-SPH 算法得到的壓力耗散程度隨時間變化的曲線。

圖3 2 種耗散算法的耗散率Fig.3 Dissipation rate of two diffusive algorithms

從圖中可以看出：在模擬的初期，因重力施加的初始效應(yīng)，壓力耗散率曲線呈現(xiàn)出振蕩的特點。隨著模擬時間的增加，2 種算法的耗散率都呈現(xiàn)出持續(xù)減小的趨勢，至t=55 s 時壓力耗散都趨近于理論值3.5%。事實上，Molteni和Colagrossi［13］已經(jīng)指出：密度耗散項算法是一種一階精度的耗散算法，精度較高，而密度正則化算法則是一種零階精度的空間均化方法，精度較低；在初始階段，壓力和速度在達(dá)到平衡的過程中會出現(xiàn)強(qiáng)烈的非線性振蕩，此時使用一階精度的密度耗散項算法和零階精度的密度正則化算法會出現(xiàn)不同的耗散效果，一階精度的耗散強(qiáng)度較高，速率較快，因此可以看出兩條耗散曲線存在差別；在達(dá)到穩(wěn)定階段后，壓力振蕩較小，壓力分布趨于靜水壓力分布，此時2 種耗散算法的耗散效果幾乎一致，因此2 條曲線趨于吻合。在其他參數(shù)相同的情況下，2種算法在經(jīng)歷了相同的耗散程度后趨于一致，證明2 種算法都達(dá)到收斂。

為了比較2 種SPH 耗散項算法在抑制壓力振蕩上的表現(xiàn)，分別選取2 種算法在5 個時間節(jié)點（t=0，12，25，40 和55 s）的流場壓力分布圖進(jìn)行了對比，如圖4 所示。各分圖中，左圖和右圖分別為密度正則化算法和δ-SPH 密度耗散項算法所得的流場壓力分布。

通過對比可以發(fā)現(xiàn)：在數(shù)值模擬的初期，如圖4（a）～圖4（b）所示，密度正則化算法與δ-SPH密度耗散項算法均能夠獲得較為光滑的壓力場。在數(shù)值模擬的中后期，如圖4（c）～圖4（e）所示，密度正則化算法獲得的流場壓力分布變得粗糙，不再保持壓力梯度分層過渡的趨勢，存在高頻壓力噪聲，壓力振蕩嚴(yán)重，且振蕩的程度隨著時間的增加而加劇了，最終導(dǎo)致數(shù)值模擬失效。與此相反，在相同的時間節(jié)點上，δ-SPH 密度耗散項算法獲得的流場壓力分布均勻、層次清晰、光滑過渡，雖然因為黏性數(shù)值耗散的原因，計算的壓力值呈現(xiàn)下降的趨勢，但是隨著時間的推移，依然能夠獲得穩(wěn)定的流場壓力。

圖4 靜水下幾個典型時刻2種算法的流場壓力分布比較Fig.4 Comparison of the pressure distribution in flow field at some typical moments in still water

以上現(xiàn)象表明：在靜水工況數(shù)值模擬中，密度正則化算法不能長期有效過濾流場的壓力噪聲，抑制壓力振蕩的效果較差；而δ-SPH 密度耗散項算法能夠長期有效地光滑流場的壓力，抑制壓力振蕩的效果較好。

從能量守恒的角度分別計算了2 種密度耗散算法中流體的動能變化情況（對動能進(jìn)行無量綱處理：將動能與初始時刻流場的總勢能作比），得到如圖5 所示的無量綱動能時歷曲線。

圖5 無量綱動能曲線Fig.5 Non-dimensional kinetic energy curves

分析圖5 可以發(fā)現(xiàn)：含有密度耗散項的δ-SPH算法與含有密度正則化的傳統(tǒng)SPH 算法相比，能夠更快地達(dá)到收斂，并且δ-SPH 算法能夠在初始效應(yīng)過后的整個數(shù)值模擬過程中基本保持動能為0，即滿足動能守恒；而傳統(tǒng)SPH 方法在數(shù)值模擬的中后期出現(xiàn)了速度波動，導(dǎo)致動能出現(xiàn)0.02%～0.05%的增長，動能不穩(wěn)定性增加。根據(jù)這個趨勢可以預(yù)測：在長時期的數(shù)值模擬中，密度正則化方法的動能守恒性會越來越差；而δ-SPH 算法的動能守恒性更好，證明其具有更好的魯棒性，更適用于長時期的數(shù)值模擬。

3.2 液艙晃蕩

為了比較2 種SPH 耗散項算法在周期性運(yùn)動條件下對流場壓力振蕩的抑制效果，對二維矩形水箱在橫搖狀態(tài)下的液體晃蕩過程進(jìn)行了數(shù)值模擬。在本算例中，密度耗散算法參數(shù)的選取均同靜水算例一致。水箱以底部中心為旋轉(zhuǎn)點做橫搖運(yùn)動，橫搖角度為正弦形式：

式中：θ 為橫搖角，rad；A 為橫搖角幅值，rad；ω為橫搖激勵頻率，rad/s。

液艙晃蕩流體固有頻率計算采用Faltinsen等［1］的公式：

為探討不同轉(zhuǎn)動角速度情況下2 種算法對壓力振蕩的抑制程度，根據(jù)數(shù)值實驗結(jié)果，分別選取并討論2 組工況：工況a——角速度較低情況；工況b——角速度較高情況。具體工況設(shè)置如下：

式中，T 為橫搖激勵周期，s。

3.2.1 角速度較低情況

對于工況a，由于邊界運(yùn)動的角速度較小，所以流體運(yùn)動基本處于低速流動狀態(tài)。為觀察在長期模擬過程中2 種耗散算法對壓力振蕩的抑制效果，模擬的總物理時間設(shè)為200 s，圖6（a）～圖6（d）分別給出了該工況下對應(yīng)于t = 0.25T，1.75T，2.25T 和2.75T 時刻2 種算法的流場壓力分布情況。各分圖中，左圖和右圖分別為密度正則化算法和δ-SPH 密度耗散項算法所得的流場壓力分布。

圖6 工況a 下幾個典型時刻2 種算法的流場壓力分布比較Fig.6 Comparison of the pressure distribution in flow field at some typical moments for case a

通過對比發(fā)現(xiàn)：2 種耗散算法在流動模擬的初期能夠獲得穩(wěn)定光滑的壓力場。如圖6（b）左圖所示，在t=1.75T 時刻，密度正則化算法流場的內(nèi)部產(chǎn)生了高頻壓力噪聲，流場壓力分布的分層狀態(tài)被破壞。隨著時間的增加，如圖6（c）～圖6（d）左圖所示，這種壓力振蕩的現(xiàn)象一直持續(xù)，其壓力分布規(guī)律與密度正則化算法計算靜水問題的規(guī)律一致。而在同樣時刻，利用δ-SPH 密度耗散項算法得到的流場壓力分布卻分層清晰，壓力分布連續(xù)。雖然由于數(shù)值耗散的原因，2 種方法獲得的壓力幅值都在衰減，但總體看來，δ-SPH 密度耗散項算法得到的流場壓力的光滑程度要優(yōu)于傳統(tǒng)SPH 密度正則化算法的結(jié)果。

2 種算法結(jié)果存在差別的原因在于：在轉(zhuǎn)動角速度較低的情況下，速度變化量不大，物理黏性還沒有完全起效，而密度正則化算法每隔一定的時間才會對流場壓力關(guān)于空間均化，一定程度上會造成壓力振蕩的積累，而該算法又不能每步施加，因為這將給數(shù)值結(jié)果帶來極大的耗散，因此，密度正則化算法得到的壓力場存在振蕩。而密度耗散項是加在連續(xù)性方程中，通過每一個時間步調(diào)節(jié)密度的振蕩來控制壓力的振蕩，因此密度耗散項算法的壓力場更為穩(wěn)定。

3.2.2 角速度較高情況

對于工況b，流動情況為中、高速狀態(tài)，橫搖周期較短。為便于比較長期流場壓力的變化情況，總物理時間取為50 s。另外，在距離液面以下0.1 m 的左側(cè)箱體壁面位置處設(shè)置一個壓力監(jiān)測點，用來監(jiān)測該點處水壓隨時間的變化；在距離水箱左側(cè)壁面0.05 m 的位置設(shè)置一個液面高度監(jiān)測點，用來實時記錄該點處水面高度隨時間的變化。

圖7 示出了2 種SPH 算法計算得到的監(jiān)測點壓力隨時間變化的曲線。由于初始時刻橫搖激勵的突然施加，壓力的數(shù)值結(jié)果呈現(xiàn)振蕩不穩(wěn)定性。為了截取穩(wěn)定有效的壓力結(jié)果，壓力值從t=2 s 開始記錄。從圖中可以發(fā)現(xiàn)，在整個數(shù)值模擬過程中，2 種耗散項算法得到的兩條曲線的相位分布一致，壓力的幅值近似一致，曲線的趨勢基本吻合，這表明2 種算法具有幾乎同等的有效獲得監(jiān)測點壓力的能力。

圖7 監(jiān)測點壓力時歷曲線Fig.7 Time histories of pressure at the monitoring point

圖8 所示為監(jiān)測點在t=20～30 s 的壓力曲線。從中可以很明顯地看出2 種耗散項算法在計算壓力中的穩(wěn)定性和振蕩程度，2 種算法計算的壓力均有一定程度的振蕩，但結(jié)果相似，表明在該組工況下2 種算法具有同等抑制壓力振蕩的效果。

圖8 監(jiān)測點在t=20～30 s 的壓力時歷曲線Fig.8 Time histories of pressure at the monitoring point（t=20～30 s）

為方便比較2 種耗散項算法在周期性運(yùn)動情況下對流場壓力振蕩的抑制效果，選取了t=2.75T，12.25T，17.5T 和24.25T 這幾個時間節(jié)點的流場壓力分布圖，如圖9（a）～圖9（d）所示。各分圖中，左圖和右圖分別為密度正則化算法和δ-SPH 密度耗散項算法所得的流場壓力分布。從流場壓力的對比中發(fā)現(xiàn)：2 種算法獲得的流場壓力分布均勻，隨著時間的增長，流場壓力的光滑程度沒有發(fā)生太大的變化，這表明2 種算法都能夠長時間起到抑制壓力振蕩、穩(wěn)定數(shù)值結(jié)果的作用。

2 種算法結(jié)果趨于一致的原因在于：在轉(zhuǎn)動角速度較高的情況下，速度的變化量相對較大，黏性開始起效，雖然2 種算法的精度不同，但獲得的結(jié)果相似。

此外，在數(shù)值模擬過程中數(shù)值耗散的一個直觀表現(xiàn)是自由液面高度的變化。為了定量比較2 種方法在自由液面高度上的耗散程度，需要事先標(biāo)記處在自由液面上的粒子，然后基于線性插值的方式計算監(jiān)測點位置的液面高度。自由液面粒子可用式（31）計算并標(biāo)記。

由式（31）不難發(fā)現(xiàn)，當(dāng)粒子處于自由液面時，粒子的支持域在自由液面處被截斷，因此該式的結(jié)果是一個比1 小的數(shù)，這里取為0.9。

在標(biāo)記的自由液面粒子中，基于線性插值的方式，計算監(jiān)測點處的液面高度，其計算方法［22］的表達(dá)式為

圖10（a）記錄了距離左側(cè)壁面0.05 m 處自由液面高度隨時間的變化情況。分析2 條曲線發(fā)現(xiàn)：因選取的角頻率遠(yuǎn)離水箱共振頻率，液面高度曲線呈現(xiàn)出周期性的特點，中后期基本上呈現(xiàn)出平穩(wěn)波動的規(guī)律，2 條曲線在整個計算過程中能夠達(dá)到90%以上的擬合程度，這在一定程度上表明了2 種算法在長時期的周期性水動力學(xué)問題中耗散程度基本持平。

圖10（b）所示為t=20～30 s 的自由液面高度隨時間變化的曲線，可見2 種算法得到的曲線相位和趨勢一致，這表明2 種算法得到的流動狀態(tài)同步，2 條曲線的幅值基本一致，表明2 種算法得到的自由液面形態(tài)和高度基本吻合，進(jìn)一步證明了2 種算法在較高速晃蕩運(yùn)動中能夠達(dá)到同樣的模擬效果。

圖10 自由液面高度時歷曲線Fig.10 Time histories of the height of free surface

4 結(jié) 論

本文不考慮人工黏性，僅考慮密度耗散，比較并分析了2 種SPH 密度耗散算法對流場壓力振蕩的抑制效果。對二維液艙的靜水和晃蕩問題進(jìn)行數(shù)值模擬，探究了不同角速度情況下2 種SPH 耗散項算法在數(shù)值結(jié)果上的差別。得到如下主要結(jié)論：

1）當(dāng)液艙處于靜止?fàn)顟B(tài)時，δ-SPH 密度耗散項算法比Shepard 密度正則化算法更能夠在長時間的模擬中穩(wěn)定流場壓力，保證動能的嚴(yán)格守恒。

2）當(dāng)液艙在較低角速度下晃蕩時，Shepard 密度正則化算法僅能在有限的時間內(nèi)過濾壓力噪聲，而δ-SPH 密度耗散項算法則能夠在長時間的模擬中光滑壓力場，后者更適用于長時間的數(shù)值計算。

3）當(dāng)液艙在較高角速度下晃蕩時，Shepard 密度正則化算法和δ-SPH 密度耗散項算法在數(shù)值模擬中均能夠持續(xù)發(fā)揮過濾壓力噪聲、抑制壓力振蕩的作用，兩者抑制壓力振蕩的程度相當(dāng)。

總之，當(dāng)液艙在角速度較低和角速度較高情況下晃蕩時，δ-SPH 密度耗散項算法均可有效抑制壓力振蕩，建議在SPH 模型中采用該密度耗散項算法。