姚學(xué)昊，黃 丹

(河海大學(xué)工程力學(xué)系，江蘇，南京 211100)

流固耦合(Fluid-structure interaction, FSI)問(wèn)題廣泛存在于科學(xué)和工程領(lǐng)域，常涉及結(jié)構(gòu)變形破壞以及流體破碎等復(fù)雜流動(dòng)現(xiàn)象[1]，往往難以獲得解析解。隨著計(jì)算力學(xué)研究發(fā)展，任意拉格朗日-歐拉(Arbitrary lagrangian-eulerian, ALE)方法[2-3]、浸入邊界法(Immersed boundary method, IBM)[4]等各種數(shù)值方法被用于此類問(wèn)題分析，然而在處理結(jié)構(gòu)大變形、準(zhǔn)確快速地追蹤自由表面和運(yùn)動(dòng)邊界等方面依然存在困難。

與網(wǎng)格類方法不同，光滑粒子動(dòng)力學(xué)(Smoothed particle hydrodynamics, SPH)方法[1,5]是一種拉格朗日無(wú)網(wǎng)格法，它適合于追蹤自由表面和運(yùn)動(dòng)邊界，且在處理大變形問(wèn)題時(shí)不會(huì)產(chǎn)生網(wǎng)格畸變。ANTOCI等[6]采用SPH 方法模擬流體與彈性結(jié)構(gòu)的相互作用問(wèn)題，獲得了較理想的結(jié)果。KHAYYER 等[7]和ZHANG 等[8]分別提出了一種不可壓縮SPH 方法以及多分辨率SPH 方法來(lái)提高流固耦合問(wèn)題的求解精度與效率。但SPH 方法在模擬固體時(shí)仍需要一些修正技術(shù)以克服拉伸不穩(wěn)定問(wèn)題[1,5]等固有缺陷。近年來(lái)，研究人員開(kāi)發(fā)了諸多SPH 與其他方法耦合的算法，如SPH 與有限元法(Finite element method, FEM)[9-10]、光滑點(diǎn)插值法(Smoothed point interpolation method, SPIM)[11]的耦合等。然而，上述算法主要針對(duì)含結(jié)構(gòu)變形的流固耦合問(wèn)題，在涉及固體損傷與破壞的相關(guān)問(wèn)題中具有局限性。近場(chǎng)動(dòng)力學(xué)(Peridynamics, PD)[12-14]是一種新興的非局部方法，采用空間積分代替?zhèn)鹘y(tǒng)微分方程，在分析不連續(xù)力學(xué)問(wèn)題時(shí)展現(xiàn)出獨(dú)特優(yōu)勢(shì)，可應(yīng)用于多種大變形和結(jié)構(gòu)破壞[15-17]、流固耦合結(jié)構(gòu)破壞[18-19]等模擬。采用PD-SPH 耦合方法來(lái)分析流固耦合問(wèn)題，能同時(shí)發(fā)揮PD 方法模擬固體破壞和SPH 方法追蹤流體運(yùn)動(dòng)界面的優(yōu)勢(shì)。

近來(lái)，部分學(xué)者開(kāi)展了PD 與SPH 方法結(jié)合的初步研究。ZHOU 等[20]結(jié)合兩種方法提出SPD(Smoothed peridynamics)方法并應(yīng)用于大變形與斷裂問(wèn)題。FAN 等[21-22]提出一種PD-SPH 耦合模型模擬爆炸引起的土壤破碎問(wèn)題；SUN 等[23]最近提出一種基于移動(dòng)虛粒子的PD-SPH 耦合方法，并應(yīng)用于流固耦合問(wèn)題。

本文建立一種能求解流固耦合及結(jié)構(gòu)破壞的新型PD-SPH 耦合模型。針對(duì)流-固界面處理，提出一種基于虛粒子和排斥力的耦合算法，既能防止流體粒子穿透交界面，又能對(duì)交界面處流體粒子進(jìn)行邊界缺陷修正，提高計(jì)算精度。通過(guò)模擬流體作用下的彈性結(jié)構(gòu)變形問(wèn)題，驗(yàn)證了耦合方法的可行性和有效性，并進(jìn)一步開(kāi)展了涉及結(jié)構(gòu)斷裂破壞的復(fù)雜流固耦合問(wèn)題全過(guò)程仿真分析。

1 數(shù)值模型

1.1 固體模型

1.1.1 常規(guī)態(tài)型PD 基本方程

如圖1 所示，將空間域R離散為一系列包含物性信息且具有相互作用f的PD 粒子，在任意時(shí)刻t，粒子X(jué)i的運(yùn)動(dòng)方程可寫(xiě)為：

圖1 PD 模型示意圖Fig. 1 Sketch for PD model

PD 理論分為鍵型(bond-based)PD 理論、常規(guī)態(tài)型(ordinary state-based)PD 理論以及非常規(guī)態(tài)型(non-ordinary state-based)PD 理論。由于鍵型與非常規(guī)態(tài)型PD 理論通常分別存在泊松比限制及數(shù)值不穩(wěn)定問(wèn)題，本文基于常規(guī)態(tài)型PD 理論進(jìn)行固體域求解。此時(shí)，點(diǎn)對(duì)相互作用f可表示為：

1.2 流體模型

1.2.1 控制方程

1.3 流-固交界面處理方案

PD 和SPH 方法均通過(guò)粒子離散計(jì)算域，故可通過(guò)基于粒子-粒子接觸的耦合方案處理流-固交界面。圖3 為本文提出的基于虛粒子和排斥力的PDSPH 耦合方案示意圖。在該方案中，PD 粒子將作為兩種類型的耦合虛粒子(邊界虛粒子和內(nèi)部虛粒子)參與SPH 計(jì)算。同時(shí)，反作用力將作用于PD粒子從而實(shí)現(xiàn)SPH 粒子對(duì)PD 粒子的影響。

該耦合方案具體計(jì)算流程如圖4 所示。

圖4 基于虛粒子和排斥力的PD-SPH 耦合方案流程Fig. 4 Flow chart for PD-SPH coupling scheme based on virtual particle and repulsive force

2 算例分析

2.1 靜水壓力作用下的鋁板變形

為了檢驗(yàn)耦合方案的有效性、穩(wěn)定性及計(jì)算精度，首先研究了靜水壓力作用下的鋁板變形問(wèn)題[9]。如圖5 所示，在寬1.0 m 的水箱中，箱內(nèi)水深H= 2.0 m，密度為 1000 kg/m3；水箱底部為厚度d= 0.05 m 的鋁板，其密度是 2700 kg/m3，彈性模量67.5 GPa，泊松比0.34。

圖5 靜水壓力下的鋁板示意圖Fig. 5 Sketch for aluminium plate under hydrostatic pressure

初始時(shí)刻，鋁板突然受到水柱的靜水壓力載荷而發(fā)生變形。經(jīng)過(guò)短時(shí)間的初始振蕩，系統(tǒng)最終達(dá)到具有特定靜態(tài)變形的平衡狀態(tài)。根據(jù)理論解[9]，2.0 m 高水柱產(chǎn)生的靜水壓力載荷作用下鋁板跨中撓度大小為 -6.85×10-5m。在模擬計(jì)算中：人工粘性參數(shù)απ= 0.03 ； 參考聲速c0=190 m/s；光滑長(zhǎng)度h=1.33Δx；近場(chǎng)范圍半徑 δ=3Δx；時(shí)間步長(zhǎng) Δt=1×10-6s；總模擬時(shí)間為1 s。流體粒子與固體粒子的初始間距一致，粒子空間分辨率分別取L/Δx=100、L/Δx=80和L/Δx=60，對(duì)應(yīng)的粒子間距為0.01 m、0.0125 m 以及0.0167 m，此時(shí)SPH 粒子總數(shù)分別為21060、13648、7836，PD 粒子總數(shù)分別為525、340、195。本算例在配有主頻為2.7 GHz 的Intel CPU 的計(jì)算機(jī)上完成，計(jì)算總耗時(shí)分別為8 h 31 min、6 h 4 min、4 h 32 min。

圖6 給出了模擬得到的t=1 s 時(shí)不同空間分辨率(粒子間距)下流體壓力云圖以及放大1000 倍后的鋁板變形圖?？梢?jiàn)，基于PD-SPH 的流固耦合求解方法在不同空間分辨率條件下均可得到光滑的流體壓力場(chǎng)和鋁板撓度場(chǎng)。

圖6 t = 1 s 時(shí)不同空間分辨率下流體壓力云圖和鋁板變形Fig. 6 The fluid pressure contours and the deflection of the aluminium plate at t = 1 s with different spatial resolutions

圖7 為三種不同空間分辨率下鋁板跨中撓度的時(shí)間歷程，表1 則給出了不同分辨率下跨中撓度數(shù)值解及其與解析解的相對(duì)誤差。結(jié)合圖7 與表1 可以看出，空間分辨率L/Δx=60時(shí)，由于粒子數(shù)較少，PD 和SPH 方法本身的計(jì)算精度較低；隨著空間分辨率的增加，PD-SPH 耦合方法的跨中撓度計(jì)算結(jié)果逐漸向靜態(tài)理論解析解收斂，并最終在L/Δx=100時(shí)，與解析解基本一致，相對(duì)誤差僅為3.07%。由此表明：本文提出的PD-SPH 耦合算法能有效、準(zhǔn)確模擬準(zhǔn)靜態(tài)流固耦合問(wèn)題。

表1 不同空間分辨率下鋁板跨中撓度計(jì)算誤差Table 1 Computational error of mid-span deflection of aluminium plate with different spatial resolutions

圖7 鋁板跨中撓度Fig. 7 Mid-span deflection of aluminium plate

2.2 潰壩水流沖擊彈性板

如圖8 所示，寬為0.584 m 的水箱左側(cè)存在高H= 0.292 m，寬L= 0.146 m 的水柱，水柱在重力作用下倒塌產(chǎn)生潰壩流，并沖擊水箱中高h(yuǎn)e=0.08 m，厚度為w= 0.012 m 的彈性板，隨后發(fā)生劇烈的流固耦合作用。該算例中，水的密度為1000 kg/m3；彈性板的密度為 2500 kg/m3，彈性模量為1.0 MPa，泊松比為0.0。模擬中，人工粘性參數(shù)設(shè)為 απ= 0.02 ； 參考聲速c0=35 m/s；粒子間距Δx= 0.002 m，對(duì)應(yīng)的SPH 粒子與PD 粒子總數(shù)分別為12744 和246；光滑長(zhǎng)度h=1.5Δx；近場(chǎng)范圍半徑 δ=3Δx； 時(shí)間步長(zhǎng)為 Δt=5×10-6s；總模擬時(shí)間為0.75 s。使用與算例2.1 相同的計(jì)算機(jī)配置，計(jì)算總耗時(shí)為2 h 20 min。

圖8 潰壩水流沖擊彈性板的初始條件示意圖Fig. 8 Sketch for initial conditions in dam-break flows impacting the elastic plate

為了對(duì)PD-SPH 結(jié)果進(jìn)行定量驗(yàn)證，在PDSPH 模擬中計(jì)算了彈性板頂端A點(diǎn)(如圖8 所示)的水平位移時(shí)間歷程，與其它已有文獻(xiàn)結(jié)果[29-32]的對(duì)比如圖9 所示。可以看出：在水流沖擊彈性板前期(0.4 s 前)，PD-SPH 方法的位移計(jì)算結(jié)果與其它方法所得結(jié)果吻合良好，且彈性板的位移最大值以及最大位移的出現(xiàn)時(shí)間基本一致，其中，PD-SPH 方法在0.244 s 時(shí)獲得彈性板水平位移最大值，約0.046 m；在0.4 s 以后，受自由表面融合現(xiàn)象的影響，不同方法得到的位移結(jié)果存在明顯差異，PD-SPH 模擬結(jié)果與粒子有限元法(PFEM)[29]、SPH-再生核粒子法(RKPM)[32]模擬結(jié)果更為接近。

圖9 A 點(diǎn)水平方向位移時(shí)程Fig. 9 Time history of horizontal displacement of point A

圖10 為不同時(shí)刻下，潰壩水流沖擊彈性板產(chǎn)生的流體飛濺與結(jié)構(gòu)變形。可以看出：0.16 s 時(shí)，水流沖擊彈性板并沿板的左側(cè)向上爬升，同時(shí)，彈性板在沖擊力作用下出現(xiàn)彎曲變形；0.26 s 時(shí)，水流完全覆蓋彈性板左側(cè)，使其產(chǎn)生較大變形；0.42 s 時(shí)，水流撞擊右側(cè)壁面形成射流，彈性板則出現(xiàn)較大幅度的回彈現(xiàn)象。同時(shí)，圖10 還給出了文獻(xiàn)中PFEM 方法[29]和SPH 方法[30]的模擬結(jié)果，比較看出：每個(gè)時(shí)刻3 種方法得到的流體自由表面形狀和彈性板變形都是相近的。此外，PD-SPH方法也計(jì)算得到了光滑的流體壓力場(chǎng)，并且與PFEM結(jié)果相似。由此表明：本文提出的PD-SPH 耦合方法對(duì)劇烈流固耦合問(wèn)題的模擬是可行的。

圖10 不同時(shí)刻彈性板變形和流體飛濺的PFEM[29]、SPH[30]和PD-SPH 結(jié)果的比較Fig. 10 Comparison of PFEM[29], SPH[30] and PD-SPH results of plate deformation and fluid splashing at different times

2.3 流體作用下的結(jié)構(gòu)破壞

為了驗(yàn)證本文PD-SPH 方法在流固耦合破壞問(wèn)題模擬方面的可行性，本節(jié)模擬結(jié)構(gòu)在潰壩流體作用下的運(yùn)動(dòng)、破壞過(guò)程。如圖11 所示，水箱長(zhǎng)0.8 m，高0.4 m，其內(nèi)部水柱和擋板的幾何尺寸、材料參數(shù)均與算例2.2 相同。為了實(shí)現(xiàn)擋板在水流作用下的斷裂破壞，假定擋板材料強(qiáng)度極限較低。本算例中，SPH 和PD 粒子總數(shù)分別為14267 和258，總模擬時(shí)間為0.8 s，計(jì)算總耗時(shí)為2 h 27 min。

圖11 流體作用下結(jié)構(gòu)破壞的初始條件示意圖Fig. 11 Sketch for initial conditions for structural failure under fluid action

圖12 為不同時(shí)刻下結(jié)構(gòu)變形破壞和流體飛濺的模擬結(jié)果，圖13 則給出了結(jié)構(gòu)質(zhì)心A(如圖11 所示)的坐標(biāo)隨時(shí)間的變化。在0.18 s 前，流體運(yùn)動(dòng)過(guò)程以及結(jié)構(gòu)變形與算例2.2 完全一致；0.18 s 時(shí)，結(jié)構(gòu)底部開(kāi)始斷裂(如圖12(a)所示)，質(zhì)心A的y坐標(biāo)時(shí)程曲線出現(xiàn)瞬時(shí)震蕩(如圖13(b)所示)；在水流的持續(xù)作用下，裂紋不斷擴(kuò)展，產(chǎn)生的應(yīng)力波也在水體內(nèi)部傳播，導(dǎo)致流體壓力場(chǎng)振蕩，0.19 s 時(shí)裂紋呈現(xiàn)明顯的V 字型，上部擋板僅剩1 個(gè)粒子與底部連接(如圖12(b)所示)；最終在0.2 s 時(shí)，上部擋板與底部完全斷開(kāi)(如圖12(c)所示)；隨后，擋板在水流作用下向水箱右側(cè)運(yùn)動(dòng)，于0.318 s 時(shí)撞擊水箱底部(如圖12(d)和圖13(b)所示)，并出現(xiàn)反彈現(xiàn)象(接觸力由式(18)計(jì)算)；經(jīng)過(guò)幾次接觸回彈，擋板開(kāi)始沿水箱底部滑動(dòng)，0.524 s時(shí)擋板撞擊水箱右側(cè)(如圖12(f)和圖13(a)所示)，并開(kāi)始向左滑動(dòng)。最終，擋板將在水箱底部保持靜止?fàn)顟B(tài)。本文模擬所得擋板運(yùn)動(dòng)過(guò)程與文獻(xiàn)[33]基本一致，由此表明：本文提出的PD-SPH 方法可應(yīng)用于流固耦合破壞問(wèn)題求解，且能夠較好地預(yù)測(cè)流體運(yùn)動(dòng)過(guò)程和結(jié)構(gòu)變形破壞以及破壞后部分結(jié)構(gòu)的運(yùn)動(dòng)。

圖12 不同時(shí)刻結(jié)構(gòu)變形破壞和流體飛濺的模擬結(jié)果Fig. 12 Simulation results of structural deformation and failure and fluid splashing at different times

圖13 結(jié)構(gòu)質(zhì)心A 的坐標(biāo)時(shí)間歷程Fig. 13 Time history of the coordinate of centroid A

3 結(jié)論

本文建立了一種新的PD-SPH 耦合模型求解流固耦合作用下的結(jié)構(gòu)破壞問(wèn)題，主要工作和結(jié)論如下：

(1)分別采用PD 方法與SPH 方法離散固體域和流體域，可充分發(fā)揮兩種方法的優(yōu)勢(shì)。通過(guò)基于虛粒子和排斥力的耦合算法處理流-固交界面，能夠?qū)α鞴探缑嫣幍牧黧w粒子施加速度邊界，修正核函數(shù)缺陷，提高計(jì)算精度。

(2)采用PD-SPH 耦合方法模擬了靜水壓力作用下的鋁板變形問(wèn)題以及潰壩水流沖擊彈性板問(wèn)題，所得結(jié)構(gòu)變形和流體運(yùn)動(dòng)過(guò)程與解析解或文獻(xiàn)結(jié)果吻合較好，驗(yàn)證了本文PD-SPH 耦合方法在模擬復(fù)雜流固耦合問(wèn)題方面的適用性。

(3)對(duì)流體沖擊作用下的結(jié)構(gòu)破壞過(guò)程分析表明：所提出的PD-SPH 耦合方法易于實(shí)現(xiàn)流固耦合作用下結(jié)構(gòu)從變形到破壞乃至破壞后部分結(jié)構(gòu)運(yùn)動(dòng)的全過(guò)程仿真模擬，可為流-固耦合結(jié)構(gòu)破壞問(wèn)題的研究提供新參考。