周學(xué)君，陳 丁，唐 軼

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

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

周學(xué)君1，2，陳 丁2，唐 軼3

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

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

光滑粒子流體動(dòng)力學(xué)(SPH)；排斥力模型；近場(chǎng)動(dòng)力學(xué)(PD)；固壁邊界; Couette流

1 研究背景

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

SPH方法在模擬自由表面流問(wèn)題時(shí)，由于粒子的Lagrangian特性能夠自動(dòng)地獲得自由表面邊界，而不需要額外地添加邊界條件，非常方便。但對(duì)于固壁邊界條件的處理仍然存在著一定的困難，這是因?yàn)樵谶吔缟匣蛘哙徑吔缣?，?jì)算會(huì)被邊界截?cái)鄬?dǎo)致不能完全覆蓋整個(gè)區(qū)域，產(chǎn)生較大誤差。目前使用較為廣泛的SPH固壁邊界處理方法有邊界排斥力法、鏡像虛粒子法和耦合邊界粒子法，但各有優(yōu)缺點(diǎn)。

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

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

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

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

2.1 SPH方法的近似

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

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

(1)

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

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

(2)

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

粒子i處場(chǎng)函數(shù)的粒子近似式可寫(xiě)成

(3)

(4)

2.2 控制方程的SPH離散

Lagrangian描述下的流體動(dòng)力學(xué)控制方程應(yīng)用SPH近似后，離散化的SPH公式為：

(5)

(6)

(7)

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

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

(8)

3 固壁邊界施加模型

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

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

(9)

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

(10)

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

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

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

(11)

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

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

4 算例分析

本節(jié)中通過(guò)2個(gè)流體動(dòng)力學(xué)的經(jīng)典算例——Couette流和潰壩，來(lái)驗(yàn)證本文所研究的邊界處理模型的有效性。

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

4.1 Couette流

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

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

(12)

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

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

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

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

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

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

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

表>1時(shí)水流前端歸一化的平均速度的試驗(yàn)值與模擬值比較

注：相對(duì)誤差=(試驗(yàn)值-模擬值)/試驗(yàn)值×100%

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

5 結(jié) 論

針對(duì)SPH固壁邊界處理的邊界排斥力法，本文提出一種新的排斥力模型。該模型的思路來(lái)源于PD方法中描述顆粒間接觸作用的短程排斥力，提出的排斥力模型簡(jiǎn)單，參數(shù)易于設(shè)定，利于SPH固壁邊界條件的施加。從Couette流和潰壩的數(shù)值算例結(jié)果來(lái)看，本文提出的排斥力模型能夠很好地解決粒子非物理穿透邊界，能較客觀地反映真實(shí)的邊界作用，粒子在排斥力的作用下運(yùn)動(dòng)規(guī)律，分布有序，表明該排斥力模型能夠在對(duì)流場(chǎng)不產(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]黃 丹，章 青，喬丕忠，等. 近場(chǎng)動(dòng)力學(xué)方法及其應(yīng)用[J]. 力學(xué)進(jìn)展，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]胡祎樂(lè)，余 音，汪 海. 基于近場(chǎng)動(dòng)力學(xué)理論的層壓板損傷分析方法[J]. 力學(xué)學(xué)報(bào)，2013，45(4): 624-628. [14]章 青， 顧 鑫， 郁楊天. 沖擊荷載作用下顆粒材料動(dòng)態(tài)力學(xué)響應(yīng)的近場(chǎng)動(dòng)力學(xué)模擬[J]. 力學(xué)學(xué)報(bào)，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]韓亞偉，強(qiáng)洪夫，趙玖玲，等. 光滑粒子流體動(dòng)力學(xué)方法固壁處理的一種新型排斥力模型[J]. 物理學(xué)報(bào)，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

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

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

10.11988/ckyyb.20160375

2017,34(7):54-59

O242

A

1001-5485(2017)07-0054-06