周學君，陳 丁，唐 軼

(1.黃岡師范學院 數(shù)理學院，湖北 黃岡 438000; 2.河海大學 力學與材料學院，南京 210098；3.云南民族大學 數(shù)學與計算機科學學院，昆明 650500)

一種新型SPH固壁邊界處理的排斥力模型

周學君1，2，陳 丁2，唐 軼3

(1.黃岡師范學院 數(shù)理學院，湖北 黃岡 438000; 2.河海大學 力學與材料學院，南京 210098；3.云南民族大學 數(shù)學與計算機科學學院，昆明 650500)

邊界排斥力法是光滑粒子流體動力學(SPH)固壁邊界處理的方法之一，但由于缺乏統(tǒng)一的排斥力模型而制約其廣泛應用?？紤]將近場動力學(Peridynamics，PD)中描述顆粒間接觸作用的短程排斥力引入到固壁邊界處理模型中，提出一種新型SPH方法邊界排斥力模型。通過Couette流和潰壩2個算例驗證了排斥力模型的有效性。排斥力表達式簡單，參數(shù)易于給定，為SPH方法中固壁邊界處理提供新思路。

光滑粒子流體動力學(SPH)；排斥力模型；近場動力學(PD)；固壁邊界; Couette流

1 研究背景

光滑粒子流體動力學(Smoothed Particle Hydrodynamics，SPH)是一種產(chǎn)生最早、發(fā)展最成熟的純Lagrangian形式的無網(wǎng)格計算方法。SPH方法最初是由Lucy(1977)[1]，Gingold等[2]提出用來解決三維開放空間的天體物理學問題，之后經(jīng)過不斷的改進和完善，已經(jīng)廣泛應用到流體力學、固體力學和生物力學等領域。SPH方法利用有限數(shù)量的粒子來離散問題域，粒子之間通過核函數(shù)來建立聯(lián)系，在模擬大變形問題時，不存在網(wǎng)格類算法(如FEM)中因出現(xiàn)網(wǎng)格畸變或纏繞而導致的算法精度降低和失敗等難題[3]。

SPH方法在模擬自由表面流問題時，由于粒子的Lagrangian特性能夠自動地獲得自由表面邊界，而不需要額外地添加邊界條件，非常方便。但對于固壁邊界條件的處理仍然存在著一定的困難，這是因為在邊界上或者鄰近邊界處，計算會被邊界截斷導致不能完全覆蓋整個區(qū)域，產(chǎn)生較大誤差。目前使用較為廣泛的SPH固壁邊界處理方法有邊界排斥力法、鏡像虛粒子法和耦合邊界粒子法，但各有優(yōu)缺點。

邊界排斥力法最早于1994年由Monaghan[4]提出，通過在固壁邊界上布置一定數(shù)量的虛粒子，使之對鄰近邊界的真實粒子產(chǎn)生適當?shù)呐懦饬?，從而防止鄰近邊界的真實粒子非物理穿透邊界。邊界排斥力法雖然不受固壁邊界形狀的影響，容易實現(xiàn)，但排斥力是人為給定的，沒有統(tǒng)一的排斥力模型，并且參數(shù)也不好確定。Libersky等(1993)[5]首創(chuàng)在邊界上鏡像布置虛粒子的方法來施加邊界條件；在此基礎之上，Randles等(1996)[6]對鏡像虛粒子法進行改進，要求所有鏡像虛粒子具有邊界處等值的變量值，采用粒子近似將邊界條件施加到邊界粒子上。鏡像虛粒子法雖然守恒性比較好，但鏡像粒子的生成比較麻煩，特別是固壁邊界形狀不太規(guī)則時更是不易。Liu等(2002[7],2003[8])提出耦合邊界粒子法來處理固壁邊界，在固壁處布置多層的邊界虛粒子，并對這些邊界虛粒子賦予與實粒子相同的質(zhì)量、密度、壓力等參數(shù)。在計算過程中，耦合邊界虛粒子參與到實粒子的守恒方程求解過程中，從而實現(xiàn)固壁邊界條件。耦合邊界粒子法的優(yōu)點是可以降低SPH 方法的邊界效應，但有時不足以防止實粒子非物理穿透邊界。

本文針對邊界排斥力法的不足，將近場動力學(Peridynamics，PD)中描述顆粒間接觸作用的短程排斥力引入到固壁邊界處理模型中。

PD方法是由美國Sandia國家實驗室的Silling[9]提出的一種非局部、無網(wǎng)格物質(zhì)點類方法[10-11]，已經(jīng)在固體材料和結構的靜動力變形以及破壞分析中成功應用[12-18]。在PD方法中，單顆粒視為由許多物質(zhì)點組成，單顆粒內(nèi)部物質(zhì)點間與分屬不同顆粒的物質(zhì)點間均存在非局部作用力，該作用力形式簡單，參數(shù)易于給定。邊界排斥力法中邊界虛粒子與實粒子之間的作用力，與顆粒間短程排斥力相似。本文主要工作是將短程排斥力模型引入到SPH方法的固壁邊界處理中，并通過編程和具體算例的計算分析，證明邊界排斥力法處理SPH固壁邊界的有效性。

2 SPH方法的近似和控制方程的離散

2.1 SPH方法的近似

在SPH方法中，將連續(xù)體離散成有限個粒子，通過對某個粒子的支持域內(nèi)其他粒子的加權求和獲得該點的數(shù)值解。SPH方法近似分2步進行：第1步是積分近似，第2步是粒子近似。

函數(shù)f(x)在問題域Ω內(nèi)的積分近似表示式為

(1)

式中：x，x′為包含在問題域Ω內(nèi)點的坐標向量；dx′為x處無窮小體元；h為光滑長度；W為光滑函數(shù)(核函數(shù))，本文采用三次樣條函數(shù)[3]為光滑函數(shù)。

式(1)可寫成離散化的粒子近似式，即

(2)

式中：mj，ρj分別表示粒子j的質(zhì)量和密度，j=1,2,…,N；N為在x處粒子的支持域內(nèi)的粒子總數(shù)。

粒子i處場函數(shù)的粒子近似式可寫成

(3)

(4)

2.2 控制方程的SPH離散

Lagrangian描述下的流體動力學控制方程應用SPH近似后，離散化的SPH公式為：

(5)

(6)

(7)

式中：α表示Cartesian坐標分量x，y和z；在同一項中重復出現(xiàn)的下標(i或j)表示Einstein求和約定；ρ，v，p分別為密度、速度和壓力；Πij為人工黏度項；f為體力；D/Dt代表物質(zhì)導數(shù)。

式(5)和式(6)中采用對稱結構，可以減少由于粒子不一致問題產(chǎn)生的誤差[19]；另外，式(6)中Πij的作用是消除SPH方法在模擬流體動力學問題時產(chǎn)生的數(shù)值不穩(wěn)定性，本文采用的人工黏度為文獻[20]中的形式。式(6)中壓力p是由弱可壓縮流體的狀態(tài)方程[6]得出，即

(8)

3 固壁邊界施加模型

利用邊界排斥力法來處理固壁邊界，需要先在邊界上布置一定數(shù)量的虛粒子，然后構造合理的排斥力公式來模擬真實邊界力?？紤]將PD方法中的顆粒間接觸作用的短程排斥力引入到模型之中，作為邊界虛粒子對實粒子的作用力，達到施加邊界條件的目的。

PD理論認為物質(zhì)點間的相互作用力具有非局部特征，這種作用力強化了近距離物質(zhì)點間的相互作用，弱化了遠距離物質(zhì)點的相互作用，恰好吻合邊界排斥力的分布特征。對于無黏性顆粒材料，PD方法采用短程排斥力描述分屬不同顆粒的物質(zhì)點間相互作用，本文將這種短程排斥力引入邊界排斥力模型中，并采用文獻[14]中彈性排斥力形式,即

(9)

式中：x，x′為物質(zhì)點的位置矢量；d為物質(zhì)點x和x′之間的短程排斥力作用的臨界距離，表征當不同的物質(zhì)點間距離

(10)

式中β為常量參數(shù)，取值在1.0左右。

若將布置在邊界上的虛粒子和實粒子看作不同的物質(zhì)點，d的大小等于光滑長度h的a倍，則邊界虛粒子會對在其影響域內(nèi)的實粒子沿著兩粒子中心線方向上產(chǎn)生排斥力，該力的大小會隨著兩粒子的距離增大而減小，直到兩粒子的距離超過d后作用力消失。

如圖1所示，當實粒子B成為邊界虛粒子A影響域內(nèi)的粒子時，則會在沿著兩粒子的中心線處虛粒子A對實粒子B產(chǎn)生一個作用力，即

(11)

圖1 邊界虛粒子對實粒子的排斥力Fig.1 The repulsive force of virtual particles on the boundary to real particles

傳統(tǒng)SPH方法處理固壁的邊界排斥力模型復雜且參數(shù)值的確定受人為因素影響較大，而上述改進的邊界排斥力公式形式簡單，參數(shù)值易于給定。

4 算例分析

本節(jié)中通過2個流體動力學的經(jīng)典算例——Couette流和潰壩，來驗證本文所研究的邊界處理模型的有效性。

這里值得一提的是，對于初始粒子的數(shù)量和間距的設置，需要根據(jù)問題域的幾何尺寸合理設定。過多的粒子雖然能保證精度但需要更多的計算量，過少的粒子又會造成計算精度降低。對于本文算例中初始粒子的設定，課題組均進行了不同數(shù)量粒子的數(shù)值試驗；本文所采用SPH代碼是由FORTRAN語言編寫，在CPU為Intel Core i7、內(nèi)存大小為32 G的臺式電腦上運行。通過比較數(shù)值試驗的CPU花費時間和計算精度，找到最佳的初始粒子配置。

4.1 Couette流

在Couette流模型中，初始靜止的流體夾在2塊固定且水平放置于間距為l的無限大平板中，流體由于上平板突然以恒定速度V0水平方向運動而產(chǎn)生流動，最終達到穩(wěn)定狀態(tài)。流體的尺寸是0.5 mm×1 mm，在SPH模擬中設置流體粒子的數(shù)量為20×40，初始間距為2.5×10-5m，流體密度ρ=1.0×103kg/m3，平板間距l(xiāng)=1.0×10-3m，V0=1.25×10-5m/s，見圖2。

Couette流中流體水平速度Vx與時間t相關的理論級數(shù)表達式[21]為

(12)

式中：y為流體粒子的豎向位置；μ為運動黏度，這里取1.0×10-6m2/s，Reynolds數(shù)Re=1.25×10-2。在SPH模擬中，時間步長為1×10-5s，以流體最前端一列間隔選取的20個粒子的速度為參照，經(jīng)過50 000步計算后流體的速度達到穩(wěn)定狀態(tài)。邊界排斥力公式(式(11))中，h=2.5×10-5m，a=1.0，β=1.0。平板采用非滑移邊界條件，并設定當流體沿平板切線方向運動時，邊界虛粒子具有與流體實粒子大小相等方向相反的速度。

圖3 SPH模擬與理論解的Couette流前端速度分布對比Fig.3 Comparison between the simulation solution in SPH and theory solution of front velocity distribution for Couette flow

圖3給出SPH方法和理論級數(shù)解(式(12))得到的在t=0.01，0.1，0.2 s和最終穩(wěn)定狀態(tài)t=∞時流體前端速度分布對比，可以看出兩者相當吻合，經(jīng)計算得出SPH模擬值的最大相對誤差為0.78%，表明本文研究的邊界方法中的排斥力模型能較客觀地反映真實邊界力效果。

圖4 潰壩的初始SPH粒子分布Fig.4 The initial particles distribution of dam breakin SPH simulation

圖5 H=300 mm時試驗[22]與SPH模擬在不同時刻的流場形態(tài)對比Fig.5 Comparison between experimental data[22] and SPH simulation result of flow field at different moments when H=300 mm

圖6 H=600 mm時試驗[22]與SPH模擬在不同時刻的流場形態(tài)對比Fig.6 Comparison between experimental data[22] and SPH simulation result of flow field at different moments when H=600 mm

4.2 潰 壩

圖5和圖6分別給出了H=300mm和600mm時，在選取的時間節(jié)點的試驗和SPH模擬結果的流場形態(tài)比較。兩圖中SPH模擬結果都采用壓強云圖，從整體上看，試驗和模擬結果對于自由表面運動的描述比較吻合，粒子壓強分布也符合實際情況；即使在最容易發(fā)生實粒子非物理穿透邊界的時刻[23-24]，如圖5(f)和圖6(f)，SPH模擬結果也沒有出現(xiàn)粒子穿透現(xiàn)象，且粒子在整個過程中分布有序，說明本文所研究的邊界力法能在不對流場產(chǎn)生明顯干擾的情況下有效處理固壁邊界。

對某些流場形態(tài)如水流垂直爬升高度、水流前端翻轉的空腔的位置和大小，模擬結果與試驗結果有些細節(jié)差異。其主要原因是SPH模擬中邊界粒子對水粒子的影響較大，以及近似處理的流體粒子湍流問題與真實流體存在一定的差異；另外，SPH模擬中器壁簡化為自由滑移邊界，沒有考慮氣壓等影響，這與試驗環(huán)境難以保持完全一致，也是差異產(chǎn)生的可能原因。這些差異在文獻[23-24]的SPH模擬潰壩算例中也同樣存在。

圖7 SPH模擬與試驗[22]得到的潰壩水流前端位置比較Fig.7 Comparison between experimental data[22] and SPH simulation of propagation of the surge front position of dam-break water flow

表>1時水流前端歸一化的平均速度的試驗值與模擬值比較

注：相對誤差=(試驗值-模擬值)/試驗值×100%

綜合試驗值和SPH模擬值的流場形態(tài)、水流前端位置和速度的比較，可以說明本文所研究的固壁邊界處理方法在自由表面流問題中是有效的。

5 結 論

針對SPH固壁邊界處理的邊界排斥力法，本文提出一種新的排斥力模型。該模型的思路來源于PD方法中描述顆粒間接觸作用的短程排斥力，提出的排斥力模型簡單，參數(shù)易于設定，利于SPH固壁邊界條件的施加。從Couette流和潰壩的數(shù)值算例結果來看，本文提出的排斥力模型能夠很好地解決粒子非物理穿透邊界，能較客觀地反映真實的邊界作用，粒子在排斥力的作用下運動規(guī)律，分布有序，表明該排斥力模型能夠在對流場不產(chǎn)生明顯干擾的情況下有效地處理固壁邊界難題。

[1] LUCY L B. A Numerical Approach to the Testing of the Fission Hypothesis[J]. The Astronomical Journal，1977，82(12): 1013-1024.

[2] GINGOLD R A，MONAGHAN J J. Smoothed Particle Hydrodynamics：Theory and Application to Nonspherical Stars[J]. Monthly Notices of the Royal Astronomical Society，1977，181(3): 375-389.

[3] LIU M B，LIU G R. Smoothed Particle Hydrodynamics (SPH): An Overview and Recent Developments[J]. Archives of Computational Methods in Engineering，2010，17(1):25-76.

[4]MONAGHAN J J. Simulation Free Surface Flows with SPH[J]. Journal of Computational Physics， 1994，110(2): 399-406.

[5] LIBERSKY L D，PETSCHCK A G，CARNEY T C，etal. High strain Lagrangian Hydrodynamics: A Three-dimensional SPH Code for Dynamic Material Response[J]. Journal of Computational Physics，1993，109(1): 67-75.

[6] RANDLES P W，LIBERSKY L D. Smoothed Particle Hydrodynamics: Some Recent Improvements and Applications[J]. Computer Methods in Applied Mechanics and Engineering， 1996，139(1): 375-408.

[7] LIU M B, LIU G R. Smoothed Particle Hydrodynamics: A Meshfree Particle Method[M]. Singapore: World Scientific Publishing Co. Pte. Ltd., 2003: 138-141.

[8] LIU M B，LIU G R，LAM K Y. Investigations into Water Mitigations Using a Meshless Particle Method[J]. Shock Waves，2002，12(3):181-195.

[9] SILLING S A. Reformulation of Elasticity Theory for Discontinuities and Long-range Forces[J]. Journal of the Mechanics and Physics of Solids， 2000， 48(1): 175-209.

[10]SILLING S A， EPTON M， WECKNER O，etal. Peridynamic States and Constitutive Modeling[J]. Journal of Elasticity， 2007， 88(2): 151-184.

[11]黃 丹，章 青，喬丕忠，等. 近場動力學方法及其應用[J]. 力學進展，2010，40(4)：448-459. [12]KILIC B. Peridynamic Theory for Progressive Failure Prediction in Homogeneous and Heterogeneous Materials[D]. Tucson,USA: The University of Arizona， 2008.

[13]胡祎樂，余 音，汪 海. 基于近場動力學理論的層壓板損傷分析方法[J]. 力學學報，2013，45(4): 624-628. [14]章 青， 顧 鑫， 郁楊天. 沖擊荷載作用下顆粒材料動態(tài)力學響應的近場動力學模擬[J]. 力學學報，2016，48(1):56-63.

[15]WECKNER O，ABEYARATNE R. The Effect of Long-range Forces on the Dynamics of a Bar[J]. Journal of the Mechanics & Physics of Solids，2005，53(3): 705-728.

[16]SILLING S A，ASKARI E. A Meshfree Method Based on the Peridynamic Model of Solid Mechanics[J]. Computers & Structures，2005，83(17/18):1526-1535.[17]KILIC B. Peridynamic Theory for Progressive Failure Prediction in Homogeneous and Heterogeneous Materials[D]. Tucson,USA: The University of Arizona，2008.[18]HUANG Dan，LU Guang-da，QIAO Pi-zhong. An Improved Peridynamic Approach for Quasi-static Elastic Deformation and Brittle Fracture Analysis[J]. International Journal of Mechanical Sciences，2015，(94/95):111-122.

[19]MONAGHAN J J. An Introduction to SPH[J]. Computer Physics Communications，1988，48(1): 89-96.

[20]LATTANZIO J C，MONAGHAN J J，PONGRACIC H，etal. Controlling Penetration[J]. SIAM Journal on Scientific and Statistical Computing，1986，7(2): 591-598.

[21]MORRIS J P，PATRICK J F，ZHU Y. Modeling Low Reynolds Number Incompressible Flows Using SPH[J]. Journal of Computational Physics，1997，136(1): 214-226.

[23]韓亞偉，強洪夫，趙玖玲，等. 光滑粒子流體動力學方法固壁處理的一種新型排斥力模型[J]. 物理學報，2013，62(4): 326-336.

[24]LIU Hu, QIANG Hong-fu, CHEN Fu-zhen，etal. A New Boundary Treatment Method in Smoothed Particle Hydrodynamics[J]. Acta Physica Sinica, 2015, 64(9):094701.

(編輯：黃 玲)

A Repulsive Model for Solid Boundary Treatment inSmoothed Particle Hydrodynamics

ZHOU Xue-jun1，2， CHEN Ding2，TANG Yi3

(1.College of Mathematics and Physics，Huanggang Normal University，Huanggang 438000， China；2.College of Mechanics and Materials，Hohai University，Nanjing 210098，China；3.College of Mathematics and Computer Science，Yunnan University of Nationalities，Kunming 650500，China)

Boundary repulsive method is one of the methods for solid boundary treatment in smoothed particle hydrodynamics (SPH)， but the method is difficult to be widely applied due to the lack of unified repulsive model. The short-range repulsive force which describes the acting force between granules in Peridynamic (PD) is introduced to solid boundary treatment model to build a new boundary repulsive model in the framework of SPH. The reliability of the method is verified by two numerical simulation examples including Couette flow and dam-break. Moreover， the repulsive formulation is simple and the parameters are easy to be given. Therefore, the present method provides a new alternative for solid boundary treatment in SPH.

smoothed particle hydrodynamics (SPH); repulsive model; peridynamic (PD); solid boundary; Couette flow

2016-04-20；

2016-06-22

國家自然科學基金項目(61462096)；江蘇省普通高校研究生科研創(chuàng)新計劃項目(KYZZ16_0268)；黃岡師范學院高級別培育項目(201617603)

周學君(1981-)，男，湖北蘄春人，講師，博士研究生，主要從事計算力學與工程仿真研究，(電話)13477625972(電子信箱)zhouxj@hhu.edu.cn。

10.11988/ckyyb.20160375

2017,34(7):54-59

O242

A

1001-5485(2017)07-0054-06