基于SPH 模型的2 種壓力振蕩抑制算法比較分析

孫雷，程聰，鄭博學，劉昌鳳

1 大連理工大學船舶工程學院，遼寧大連116024

2 高新船舶與深海開發裝備協同創新中心，上海200240

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

0 引 言

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

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

國內外學者針對該問題進行了深入研究，并提出了2 種可行的方法：密度正則化和數值耗散技術。

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

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

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

1 控制方程

1.1 基本控制方程

假設流體不可壓縮，考慮絕熱條件，流體運動應遵循質量守恒定律和動量守恒定律。在拉格朗日坐標系下，控制方程被簡化為如下形式的連續性方程和動量守恒方程：

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

1.2 狀態方程

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

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

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

1.3.1 Shepard 密度正則化

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

1.3.2 δ-SPH 方法的耗散項

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

2 SPH 數值計算模型

2.1 SPH 基本原理

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

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

為了對控制方程進行離散，將上述積分形式的表達式轉換成對支持域內粒子的求和形式，即粒子近似。場函數及其導數的核近似表達式的離散求和形式為：

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

2.2 物理黏性

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

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

2.3 邊界條件

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

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

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

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

2.4 時間積分

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

第1 步：預測步。

第2 步：校正步。

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

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

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

3 數值算例及結果討論

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

圖2 數值模型示意圖Fig.2 Schematic diagram of numerical model

3.1 靜水情況

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

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

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

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

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

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

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

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

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

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

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

3.2 液艙晃蕩

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

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

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

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

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

3.2.1 角速度較低情況

對于工況a，由于邊界運動的角速度較小，所以流體運動基本處于低速流動狀態。為觀察在長期模擬過程中2 種耗散算法對壓力振蕩的抑制效果，模擬的總物理時間設為200 s，圖6（a）～圖6（d）分別給出了該工況下對應于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

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

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

3.2.2 角速度較高情況

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

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

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

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

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

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

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

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

由式（31）不難發現，當粒子處于自由液面時，粒子的支持域在自由液面處被截斷，因此該式的結果是一個比1 小的數，這里取為0.9。

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

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

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

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

4 結 論

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

1）當液艙處于靜止狀態時，δ-SPH 密度耗散項算法比Shepard 密度正則化算法更能夠在長時間的模擬中穩定流場壓力，保證動能的嚴格守恒。

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

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

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