基于SPH方法的液體晃動數值模擬

方 政 李 煜 蘭慶琳 周玲玲

（河海大學 力學與材料學院，南京 210098）

液體晃動［1］是指容器內的液體由于外加擾動而產生的帶有自由表面波動的運動.在很多工程領域存在液體晃動問題，例如：升船機的承船廂內水體波動、水庫中水體由于地震引起的晃動以及貯液箱體內流體的晃動等.

早期研究液體晃動大多采用線性理論，即將水體簡化為單擺系統.當液體晃動幅度較小時，簡化模型的模擬效果較好.但隨著晃動幅度增加，會產生非線性效應，自由液面難以用線性理論精確地預測出來.目前研究大幅晃動問題，多采用數值模擬方法.Har－low最早提出MAC法［2］，并首次將其應用到模擬大幅晃動問題中.Liu等應用VOF法模擬了帶破碎自由液面的非線性晃動問題［3-4］.另外，邊界元法、有限元法也被廣泛的應用［5-6］.但是 MAC方法計算時需要的存儲空間大，計算時間長［7］；VOF法難以精確的得到自由液面的位置［8］.邊界元法、有限元法等在計算大變形問題時，網格會出現畸形，影響計算精度［7］.光滑粒子流體力學方法（Smoothed particle hydrodynamics，SPH）［9］是一種拉格朗日描述下的無網格粒子法.它不需要事先劃分網格，能夠適應大變形問題［10］.此外，SPH采用顯式的狀態方程計算壓力，顯式的時間積分法求解控制方程，能夠提高計算效率.因此比較適合模擬帶有自由液面的大幅晃動問題.

1 基本理論

光滑粒子流體力學（SPH）是一種無網格粒子化的數值模擬方法，其控制方程是拉格朗日型的Navier－Stokes方程組.該方法的基本思想是將計算域內流體離散成若干個粒子，并賦予每個粒子物理特性，例如，密度、質量、速度等；然后將描述這些物理特性的函數及其導數用核函數積分近似，得到這些函數的核近似方程；將核近似方程代入拉格朗日型的N-S方程組，得到粒子運動的常微分方程組；最后，采用顯示積分法求解得到粒子的各物理量.

1.1 SPH中函數的近似方法及控制方程

SPH方法中，利用核函數將任意函數f（x）進行積分近似，即核函數近似

式中，Ω為包含x的積分域，W 為核函數（也稱為光滑函數），h為定義核函數影響區域的光滑長度.

SPH方法模擬流體運動的過程中，光滑函數［11］影響著函數積分近似的準確性，決定了函數表達式的計算精度和效率.由于高斯型函數［12］充分光滑、穩定并且精確度較高，所以本文采用高斯型函數，表達式為

式中，R 為粒子間的相對距離；αd＝1／πh2.

拉格朗日型的N－S方程組包括連續性方程和動量方程，粒子化后得連續性方程

式中，ρi為粒子密度；N為粒子i支持域內的粒子總數；mj為粒子質量；vij為中心粒子與支持域內粒子的速度差，vij＝vi-vj；β為坐標分量.

動量方程

式中，vi為粒子速度；α，β 為左邊分量；pi，pj分別為中心粒子壓強和支持域內其他粒子壓強；μi，μj為粒子的粘性系數為剪應變率，具體如下式：

1.2 SPH的模擬精度

實際上，SPH模擬方法是用弱可壓縮性流體模型代替不可壓縮流體，因此加入了人工可壓縮率.這種人工可壓縮率引起的誤差是O（M2）（M為馬赫數，M＝v／c；v為流場中速度最大值；c為模型聲速）.模型中的壓力按照可壓縮流體的狀態方程計算.

式中，c為模型聲速；ρ0為參考密度；γ是常數，取為7［9］；ρ為實體粒子密度.

由上述可知模型聲速影響著模擬精度，模型聲速既要保證使人工可壓縮流體的特性與真實流體充分接近，同時也不宜過大，以保證時間步的增量在容許范圍之內.一般采用下式確定c值［9］

式中，c為模型聲速；vb為整體流速；δ為相對密度差；υ為動力粘滯系數；F為體積力；L為特征長度.

2 SPH的邊界處理方法

2.1 自由表面處理

SPH方法在處理自由表面時，由于粒子數量不足，計算得到的粒子密度會小于實際值，這樣在利用狀態方程計算壓強時，得到的結果會與實際值不符，影響計算的精確度［13］.因此，需要采取一定的方法，首先捕捉到自由表面的粒子，然后賦予這些粒子水的標準密度值，從而保證計算的穩定性.

通常，若粒子密度滿足Koshizuka條件［14］，即被看作是自由表面粒子.公式為

式中，ρ為粒子密度的計算值，β為小于1的參數，一般在0.8到0.98之間取值.對于弱沖擊性流動，在計算出的密度穩定性較好的情況下，利用（7）式追蹤自由表面粒子效率較高，模擬精度較好.由于液體有限振幅下的晃動屬于弱沖擊性流動，本文采用該方法處理自由表面.

2.2 固體邊界

SPH方法處理固體邊界問題時，比較常用的有邊界排斥力法和鏡像粒子法.

Monaghan提出了邊界排斥力方法［15］，即在固體邊壁上布置一系列的邊界粒子，這些虛擬粒子會給靠近它們的實體粒子一個沿軸向的排斥力，以達到防止實體粒子穿透邊界的目的.計算排斥力的公式

式中，D為常數，一般與流域速度最大值的平方量級相當；r0為截止距離；rij為邊界粒子指向實體粒子的矢徑；xij為邊界點與實體粒子矢徑的坐標分量；n1，n2為常數，通常分別取為12和4.實際模擬中，邊界排斥力法有時不能完全阻擋實體粒子穿透邊界，并且邊界排斥力會對實體粒子產生擾動，影響計算的準確性.鏡像粒子法［16］彌補這樣缺陷，其基本思想是：當內側實體粒子到邊界的距離小于κhi時，在邊界外側布置與之對應的虛粒子，這些虛粒子和實體粒子到邊界的距離相等，質量、密度、壓力一樣，速度大小相等，垂向速度相反，切向速度相同，如圖1所示.另外，鏡像虛粒子可以提高邊界處壓力計算的精確度，所以兩者往往結合使用.本文采取兩者結合的方法處理固體邊界.

圖1 粒子分布示意圖

3 箱內水體晃蕩模擬

采用矩形箱體，尺寸如圖2所示，箱體長度L＝0.8m，箱內水深為h＝0.3m.對其進行自由晃動和一階模態晃動的數值模擬.

3.1 自由晃動

將計算區域剖分成個80×30＝2400個實體粒子，粒子初始速度v0，粘性系數為v＝1.0×10-6m2／s.箱體不受外部激勵影響，橫向加速度為零，僅受重力影響，重力加速度g＝9.81m／s2.為防止實體粒子穿透邊界，采用排斥力邊界和鏡像粒子法，在邊界上布置183個虛粒子.時間步長Δt＝1.0×10-4s，采用蛙跳積分法.初始時刻，壓強按靜水壓強分布；然后根據狀態方程（5）反算出粒子密度，計算過程中始終用狀態方程計算壓強.

模擬時，先施加一個水平方向的加速度a′x，大小為2.4525m／s2，計算達到穩定時，形成傾斜液面，然后去掉a′x，使其只在重力作用下形成自由晃動.傾斜液面如圖2所示.

圖2 貯箱尺寸及傾斜液面示意圖

圖3給出了x＝0m和x＝0.8m處的液面波動情況，計算得到晃動周期T＝1.051s，頻率f＝0.9509Hz.根據勢流理論［17］所得固有頻率為0.9520Hz.結果相近，說明本文建立的數值模型較為合理，用來模擬液體晃動是可行的.

圖3 x＝0.0m和x＝0.8m處液面波動

3.2 一階晃動

對靜止的液體施加橫向的外部激勵：ax＝-bω2sin（ωt），b為外部激勵引起的震蕩幅值，本文取為0.062m，ω取為5.611.其他參數及計算方法均與上例相同.計算10s工況，共耗費機時535.6s（cpu time）.

圖4是不同時刻自由液面形狀，從圖中可以看到，一階工況下，各時刻的自由液面圍繞著液面中心點來回晃動，在中心點處的振幅大致為0m，在貯箱的一側產生波峰，在另一側產生波谷，這是最基本的非對稱波動.

圖4 不同時刻液面形狀示意圖

圖5為x＝0.8m處自由液面的波動時程，可以看到，SPH方法得到的數據與文獻［18］中的實驗數據吻合較好.從圖5中可以較明顯地看出振動波形關于自由液面不對稱，根據線性理論，只有一、二階晃動疊加時，才會出現這種情況.這是大幅晃動中的非線性現象之一.

圖5 x＝0.8m處自由液面的波動時程

圖6是右邊壁距離箱底0.2m處的壓強時程圖，圖中的靜水壓強是指某一時刻右邊壁相對于自由液面的靜水壓強分布.可以看出在2.2s之前，吻合較好，隨著時間推移瞬時壓強呈現出脈動的現象，并且與靜水壓強相比數值相差較大，但總體在一個合理的范圍內，這也是由于水體的大幅晃動而產生的非線性現象.此外，在某些波峰和波谷處出現了二次諧波，這是由于壁面對流體的阻滯作用，產生了局部回流的原因.

圖6 右邊壁距離箱底0.2m處的壓強時程圖

4 結 語

采用SPH方法模擬了貯箱內液體大幅晃動，計算結果與相關實驗數據吻合較好.分析了自由晃動和一階晃動的液面波動和壓強分布，說明本文數值模型設定的參數較為合理，顯示了SPH方法在模擬液體晃動方面的優越性.

［1］曾江紅，王照林.航天器液體晃動動力學的研究方法概述［J］.強度與環境，1997（4）：37－43.

［2］Harlow F H，Welch J E.Numerical Calculation of Time－Dependent Viscous Incompressible Flow of Fluid with Free Surface［J］.Physics of Fluids，1965（8）：2182－2189.

［3］Liu Dongming，Lin Pengzhi.A Numerical Study of Three－Dimensional Liquid Sloshing in Tanks［J］.Journal of ComputationalPhysics，2008，8（227）：3921－3939.

［4］郭紅民，向光明，謝 洋，等.毛家河水電站溢洪道三維數值模擬［J］.水電能源科學，2013，31（3）：81－85.

［5］Nakayama T，Washizu K.The Boundary Element Method Applied to the Analysis of Two－dimensional Nonlinear Sloshing Proble［J］.International Journal for Numerical Methods in Engineering，1981，17（11）：1631－1646.

［6］曹宗楊，張燎軍.拱壩-庫水動力流固耦合作用的有限元數值研究［J］.水電能源科學，2013，31（4）：58－61.

［7］陳海舟.不可壓自由表面流的SPH法數值模擬研究［D］.天津：天津大學，2008.

［8］Hirt C W，Nichols B D.Volume of Fluid（VOF）Method for the Dynamics of Free Boundaries［J］.Comput.Phys，1981，39：201－225.

［9］Liu G R，Liu M B.Smoothed Particle Hydrodynamics－a mesh free Particle method［M］.Hong Kong：World Publishing Company，2005.

［10］杜小弢，吳 衛，龔 凱，等.二維滑坡涌浪的SPH方法數值模擬［J］.水動力學研究與進展（A 輯），2006（5）：579－586.

［11］Gingold R A，Monaghan J J.Smoothed Particle Hydrodynamics：Theory and Application to Non－Spherical Stars［J］.Monthly Notices of the Royal Astronomical Society，1977，181：375－389.

［12］Monaghan J J，Lattanzio J C.A Refined Particle Method for Astrophysical Problems［J］.Astron.Astrophys，1985，149：135－143.

［13］何 濤，周 岱.二維自由表面流動的光滑粒子動力學方法數值模擬［J］.上海交通大學學報，2011（10）：1425－1429.

［14］Koshizuka S，Nobe A，Oka Y.Numerical Analysis of Breaking Waves Using the Moving Paticle Semi－Implicit Method［J］.International Journal for Numerical Methods in Fluids，1998，26：751－769.

［15］Monaghan J J.Simulating Free Surface Flows with SPH［J］.Journal of Computational Physics，1994，110（2）：399－406.

［16］Randles P W，Libersky L D.Smoothed Particle Hydrodynamics：Some Recent Improvements and Applications［J］.Computer methods in applied mechanics and engineering，1996，139：375－408.

［17］Faltinsen O M.A Nonlinear Theory of Sloshing in Rectangular Tanks［J］.Journal of Ship Research，1974，18（4）：224－241.

［18］Antonio Huerta，Wing Kam Liu.Viscous Flow with Large Free Surface Motion［J］.Computer Methods in Applied Mechanics and Engineering，1988，69（3）：277－324.