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

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

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

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

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

0 引 言

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

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

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

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

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

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

1 控制方程

1.1 基本控制方程

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

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

1.2 狀態(tài)方程

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

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

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

1.3.1 Shepard 密度正則化

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

1.3.2 δ-SPH 方法的耗散項(xiàng)

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

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

2.1 SPH 基本原理

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

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

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

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

2.2 物理黏性

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

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

2.3 邊界條件

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

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

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

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

2.4 時(shí)間積分

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

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

第2 步：校正步。

一個(gè)時(shí)間步以后，粒子最終的物理量為

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

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

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

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

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

3.1 靜水情況

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

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

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

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

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

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

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

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

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

圖5 無量綱動(dòng)能曲線Fig.5 Non-dimensional kinetic energy curves

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

3.2 液艙晃蕩

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

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

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

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

式中，T 為橫搖激勵(lì)周期，s。

3.2.1 角速度較低情況

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

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

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

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

3.2.2 角速度較高情況

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié) 論

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

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

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

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

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