姚學昊 陳丁 武立偉 黃丹

(河海大學力學與材料學院,南京 211100)

引言

流-固耦合(fluid-structure interaction,FSI)是一種流體與固體相互作用和影響的力學問題,其廣泛存在于諸多工程領域,如水上降落、液艙晃蕩以及海浪砰擊艦船與海洋平臺等.由于流體流動、結構變形破壞以及二相間相互作用問題的復雜性,流-固耦合相關問題,特別是流-固耦合導致結構破壞問題的數(shù)值模擬一直是研究難點.已有諸多數(shù)值方法用于求解FSI 問題,其中,基于網(wǎng)格的數(shù)值方法因計算精度和效率高而被廣泛應用,如有限單元法(finite element method,FEM)[1-2]、有限體積法(finite volume method,FVM)[3]、有限差分法(finite difference method,FDM)[4]及其混合方法[5-6]等.然而,基于網(wǎng)格的數(shù)值方法在求解FSI 問題時往往需要網(wǎng)格重構、界面追蹤等特殊數(shù)值技術,增加了算法復雜性,并且在處理三維裂紋萌生、動態(tài)多裂紋擴展等強不連續(xù)問題時依然面臨挑戰(zhàn).無網(wǎng)格方法是另一類被廣泛使用于FSI 問題的數(shù)值方法,其不受網(wǎng)格約束,能在拉格朗日框架下求解流體域或固體域,在復雜流動問題模擬以及固體結構的大變形與斷裂分析等方面展現(xiàn)出獨特的優(yōu)勢,已成為研究FSI 問題的重要手段.

近場動力學(peridynamics,PD)[7-9]與光滑粒子流體動力學(smoothed particle hydrodynamics,SPH)[10]是兩種基于非局部假定的無網(wǎng)格方法,通常將計算域離散為一系列自由移動的粒子,通過求解控制每個粒子運動的積分方程或偏微分方程獲得整個系統(tǒng)的力學行為.PD 方法無需預先確定裂紋位置,可自然地描述裂紋萌生和擴展,近年來成功應用于固體斷裂[11]、沖擊破壞[12-13]等不連續(xù)力學問題模擬.針對FSI 問題,部分學者建立了不同的PD 模型,如Zhang等[14]和Zhou等[15]分別建立了用于水力壓裂模擬的PD 模型,Ren等[16]建立了模擬海冰在海浪作用下力學行為的PD 模型.然而,當前的PD 模型通常不適用于含自由表面流體的FSI 問題,涉及復雜流動的PD 相關算法仍有待進一步研究[17-18].SPH 方法適用于流體運動模擬,且能處理結構動態(tài)響應問題.近年來,Oger等[19]應用弱可壓SPH 方法對剛性楔形體入水問題進行研究,獲得了較為理想的數(shù)值結果,Sun等[20]結合δ+-SPH 與完全拉格朗日格式的SPH 方法研究了多種流體與彈性結構相互作用問題,Khayyer等[21]建立了不可壓縮SPH 與傳統(tǒng)SPH 混合方法,將其應用于波浪沖擊可變形結構問題.盡管SPH 可以單獨處理FSI 問題,但目前很少有研究考慮強耦合作用導致的固體破壞.

在PD-SPH 混合方法研究方面,Fan等[22-23]針對土壤的爆炸破碎問題,采用非常規(guī)態(tài)型PD 模型建立了一種固體粒子與流體粒子互為虛粒子的混合方法,Sun等[24]、姚學昊等[25]針對FSI 問題分別采用基于接觸力的混合方法.上述方法均使用相同的粒子分辨率離散固體域和流體域,在求解結構空間尺度遠小于流體的問題時需要較高的計算成本.Rahimi等[26]通過一種基于拉格朗日插值的PD-SPH 混合方法,實現(xiàn)了空間多分辨率離散,但未考慮多時間步長問題.

基于上述研究現(xiàn)狀,本文提出一種用于流-固耦合結構破壞模擬的多時間和空間分辨率PD-SPH 混合方法.該方法采用不同的粒子分辨率和時間步長對固體域和流體域進行計算以降低計算成本.通過典型的自由表面流體與可變形結構相互作用問題算例驗證了方法的精度和效率,并應用于流體沖擊混凝土板破壞過程分析.

1 數(shù)值模擬方法

1.1 PD 方法

傳統(tǒng)連續(xù)介質力學理論中物體內(nèi)部的相互作用由應力張量的散度 ? ·σ 表示[9],而PD 方法中固體域被離散為一系列包含質量m和密度ρ等基本信息的無拓撲關系的粒子Xi,每個粒子與其近場范圍HX(半徑為δ的球形區(qū)域)內(nèi)的其他粒子Xj間存在相互作用f

利用體積膨脹標量可將拉伸標量狀態(tài)分解為球量部分es=θx/3和偏量部分ed=e-es.此時,借鑒經(jīng)典理論中彈性應變能密度表達式,可得到適用于簡單彈性材料的常規(guī)態(tài)型PD 變形能密度函數(shù)

式中,k和α是與體積模量K和剪切模量G相關的PD 材料參數(shù),可根據(jù)PD 理論中變形能密度與傳統(tǒng)理論的應變能密度等效獲得[27].計算變形能密度函數(shù)關于拉伸標量狀態(tài)e的弗雷謝(Fréchet)導數(shù)則可得到力標量狀態(tài)t

針對斷裂問題,假定兩個粒子Xi和Xj之間的相對位置滿足一定關系時,粒子對內(nèi)粒子間的相互作用消失.可引入間斷函數(shù)表示

1.2 SPH 方法

在SPH 方法中,同樣將流體域離散為一系列粒子xi,通過核函數(shù)插值和粒子近似技術可得到離散形式的不可壓縮流體控制方程

式中,D /Dt表示物質導數(shù);xji=xj-xi為粒子在當前構形中的相對位置矢量,vij=vi-vj為相對速度矢量;h為光滑長度;p為壓力;fg,fΠ和fυ分別為體力、人工黏性力和物理黏性力;c0為參考聲速,其值應足夠大以滿足弱可壓縮假定[28]

式中,|v|max,pmax和ρ0分別為流體的最大流速、最大壓力以及參考密度.?i表示對粒子xi坐標的空間導數(shù);Wij=W(xij,h)為核函數(shù),本文選用B-樣條核函數(shù)[10]

式中,qs=|xij|/h;空間維度dim=1,2,3 分別對應αw=1,1 5/7π,3 /2π.控制方程組(15)中的連續(xù)性方程右側第二項為人工耗散項[28],用以抑制因數(shù)值噪聲產(chǎn)生的壓力振蕩,其中 δs為耗散強度控制常數(shù),常取值為0.1; 〈 ?ρ〉L為修正密度梯度

對于不可壓縮流體,粒子xi受到的物理黏性力可簡化為

式中,μ為流體動力黏度;φ=0.1h用于防止奇點的出現(xiàn).為了提高數(shù)值穩(wěn)定性,Monaghan[29]在動量方程中加入了人工黏性項fΠ以抑制數(shù)值振蕩,其表達式為

2 多分辨率PD-SPH 耦合方案

本文采用具有不同初始間距的固體粒子和流體粒子離散計算域,并通過流體粒子-虛粒子接觸耦合方案處理流-固交界面.為了更顯著地提高計算效率,利用不同的時間步長對固體和流體控制方程進行時間積分.

2.1 多分辨率空間離散

如圖1 所示,當SPH 流體粒子a靠近流-固界面時,PD 固體粒子b和c均位于SPH 粒子a的支持域Sa內(nèi),此時粒子b和c將以虛粒子b′和c′的形式加入粒子a的鄰近搜索列表,從而保證粒子a的支持域不被流-固界面截斷,實現(xiàn)SPH 邊界缺陷修正.考慮到PD 粒子與SPH 粒子不同的初始間距會導致虛粒子與流體粒子光滑長度不同,本文中令虛粒子的光滑長度等于SPH 粒子的光滑長度h,這也使得流體粒子a的支持域Sa的大小與虛粒子b′,c′的支持域Sb′,Sc′ 完全相同.此時,SPH 粒子a的控制方程變?yōu)?/p>

圖1 多分辨率PD-SPH 耦合方案Fig.1 Multi-resolution PD-SPH coupling scheme

式中,Kr和nr為用戶自定義參數(shù)[10].

需要注意的是,在該耦合方法中,虛粒子 κ′根據(jù)PD 粒子在初始構形中的位置分為表面虛粒子(如b′)和內(nèi)部虛粒子(如c′),僅被動地被流體粒子搜索,自身不進行SPH 數(shù)值積分,其位置信息來自對應的PD 粒子,密度等場變量信息則通過插值計算獲得[25].如圖1 所示,對于表面虛粒子b′,其場變量信息來自其支持域Sb′ 內(nèi)的流體粒子;內(nèi)部虛粒子c′的場變量信息則來自其支持域Sc′ 內(nèi)的流體粒子和表面虛粒子.虛粒子速度(無滑移邊界條件)和密度表達式為

對于PD 部分,界面處PD 粒子受到SPH 粒子的作用力,相當于在耦合界面處對PD 粒子施加力邊界條件.為了在流-固界面處滿足動力學條件,保證系統(tǒng)動量守恒,令SPH 粒子a對PD 粒子 κ 的作用力Fatoκ

2.2 多分辨率時間離散

考慮到PD 與SPH 的完全顯式計算,分別定義結構和流體的計算時間步ΔtS和ΔtF,其中,固體計算的時間步長應滿足[31]

通常,固體計算的時間步長ΔtS會小于流體計算的時間步長ΔtF,若采用單一時間步長會出現(xiàn)流體計算耗時過多的問題.本文將采用多時間步長方案,并通過執(zhí)行一次流體時間積分同時進行 ?=ΔtF/ΔtS次固體時間積分的方式保證流體與固體計算的時間一致性.該方案示意圖如圖2 所示.

圖2 ?=6 時多時間步長方案Fig.2 Multi-timestep scheme(?=6)

多分辨率PD-SPH 耦合計算流程如圖3 所示.

圖3 多分辨率PD-SPH 耦合計算流程Fig.3 Flow chart for multi-resolution PD-SPH coupling calculation

3 算例驗證

為驗證多分辨率PD-SPH 混合方法的精度和計算效率,首先對經(jīng)典的FSI 問題進行數(shù)值分析.隨后對流固相互作用導致的結構斷裂破壞問題進行模擬,進一步驗證耦合方案對于FSI 問題的適用性.本文算例均在Intel i7-7700CPU(主頻 3.6GHz)個人計算機上完成,所取光滑長度和近場范圍半徑分別為h=1.3dxF和δ=3.0dxS,其中dxF和dxS分別表示SPH粒子與PD 粒子的初始間距.

3.1 液柱靜壓力下的鋁板變形

首先對液柱靜壓力作用下的鋁板變形問題進行數(shù)值模擬.如圖4 所示,長L=1.0m、厚d=0.05 m,密度 ρS=2700kg/m3,彈性模量E=67.5 GPa,泊松比ν=0.34 的鋁板突然受到高H=2 m,密度ρF=1000kg/m3的水柱作用.根據(jù)文獻[33],鋁板經(jīng)過初始振蕩后達到平衡狀態(tài),其跨中撓度解析解為 - 6.85×10-5m.模擬中,參考聲速c0=60m/s,總模擬時長為 5 s.

圖4 液柱靜壓力下鋁板變形的初始條件(單位:m)Fig.4 Initial condition for the aluminum plate under hydrostatic pressure(unit:m)

針對當前問題,考慮到過大的流-固分辨率比會對計算精度和效率會產(chǎn)生一定的影響,設置了固定結構空間分辨率條件下的3 種多分辨率工況,具體參數(shù)設置如表1 所示.圖5 給出了3 種工況下鋁板跨中撓度時程,由圖可見,經(jīng)過初始振蕩,3 種工況下鋁板跨中撓度均收斂至約為 - 6.70×10-5m,與解析解相比誤差僅為2.19%.圖6 給出了t=5 s 時工況III 的流體粒子分布與壓力場以及鋁板撓度場,可以看出,此時流體壓力場與鋁板撓度場均是光滑的.上述結果表明,本文多分辨率PD-SPH 混合方法計算精度較高.

表1 液柱靜壓力下鋁板變形問題工況布置Table 1 The setup of different cases for aluminum plate deformation with hydrostatic pressure

圖5 鋁板跨中撓度時間歷程Fig.5 Time histories for mid-span deflection of aluminum plate

圖6 t=5 s 時工況III 的流體壓力場與板的撓度場Fig.6 Fluid pressure and plate deflection field in case-III at t=5 s

圖7為3 種工況下FSI 系統(tǒng)歸一化總能量[21](分別為當前和初始時刻系統(tǒng)總能量)時間歷程.由圖可見,歸一化總能量同樣存在振動變化過程,隨著FSI 系統(tǒng)進入穩(wěn)定狀態(tài),系統(tǒng)總能量最終保持穩(wěn)定.t=5 s 時,3 種工況的系統(tǒng)總能量損失分別為0.23%,0.34%以及0.49%,表明多分辨率混合方法具有較好的能量守恒特性.綜合圖5與圖7 可以發(fā)現(xiàn),鋁板位移初始振蕩的持續(xù)時間與FSI 系統(tǒng)總能量損失有關,系統(tǒng)能量損傷越少,其位移振蕩時間越長.該發(fā)現(xiàn)與Zhang等[34]所得結論一致,進一步驗證了本文的多分辨率混合方法.

圖7 FSI 系統(tǒng)歸一化總能量時間歷程Fig.7 Time histories for normalized total energy of FSI system

圖8 給出了3 種工況下全仿真過程計算時長,可見,隨著流體粒子間距增大,流體計算所需時間步長也增大,總模擬時長顯著減少.工況III 的計算速度約為工況I 的11.3 倍.由此表明,多分辨率PD-SPH混合方法能夠在保證計算精度的同時有效提高計算效率.

圖8 不同工況的計算時間Fig.8 Computational times for different cases

3.2 潰壩流體沖擊彈性閘門

Antoci等[35]對潰壩水流沖破彈性閘門問題進行了實驗研究,該問題的初始條件示意圖如圖9 所示.實驗在寬L=0.1 m、內(nèi)部水深H1=0.14 m 的水箱中進行.箱體左側剛性壁面下部放置了上端固定下端自由的彈性橡膠閘門,其寬度為s=0.005 m、高度為H=0.079 m,材料參數(shù)為:密度 ρS=1100kg/m3,彈性模量E=7.8 MPa,泊松比 ν=0.47.數(shù)值模擬中,參考聲速c0=15m/s,水體密度ρF=1000kg/m3,總模擬時長0.4 s.針對該問題的多分辨率工況設置如表2 所示.

表2 潰壩流體沖破彈性閘門問題工況布置Table 2 The setup of different cases for dam-break flow through an elastic gate

圖9 潰壩流體沖破彈性閘門的初始條件示意圖(單位:mm)Fig.9 Schematic sketch for initial conditions in the dam-break flow through an elastic gate(unit:mm)

圖10給出了工況III 中不同時刻的彈性閘門在時變水壓作用下的變形以及自由液面變化過程,并與實驗觀察結果[35]進行比較.由圖10可以看出,PD-SPH 所得彈性閘門與自由液面形狀與實驗結果基本吻合.在t=0.08s 時,閘門自由端在水壓力作用下向x軸負方向移動,形成水流出口.在t=0.16 s 時,閘門自由端位移與水流出口顯著增大.隨后,由于水位的降低,水壓力減小,水流出口逐漸減小.值得注意的是,由于本文開展的是二維模擬,實驗中的流體飛濺現(xiàn)象不會出現(xiàn).

圖10 不同時刻彈性閘門變形和自由液面演變Fig.10 Elastic gate deformation and change of fluid free surface

圖11為彈性閘門自由端水平位移時程曲線.通過與文獻[34]的SPH 模擬結果以及文獻[35]的實驗數(shù)據(jù)對比可以看出,3 種工況下所得位移結果均吻合較好.在P D-S P H 模擬中,彈性閘門于t=0.122 s時刻達到最大變形,此時其自由端水平位移約為0.043 m.圖12 給出了3 種工況的計算總時長,比較可知,僅采用多分辨空間離散方法(工況II)即可有效提高計算效率;若采用多時間與空間分辨率,即工況III,其計算速度為單一分辨率(工況I)的8.5 倍.

圖11 彈性閘門自由端水平位移時程Fig.11 Time histories of horizontal displacements at the free end of elastic gate

圖12 3 種工況的計算總時長Fig.12 Total calculation time in three cases

綜上,本文提出的PD-SPH 混合方法在不同分辨率工況下均能有效傳遞流-固相互作用信息,可用于涉及復雜流體流動和結構大變形的FSI 問題.

3.3 重力壩水力劈裂

為了驗證多分辨率PD-SPH 混合方法對于FSI 致結構破壞問題的適用性,選取經(jīng)典的帶裂縫Koyna 重力壩問題,分析其在超載作用下的斷裂響應.初始條件如圖13 所示,壩高 103 m、壩頂和壩底寬度分別為w=14.8 m和L=70m.水平裂縫位于大壩上游面 66.5 m 高程處,長度為1.93 m.大壩材料參數(shù)為:密度 ρS=2450kg/m3,彈性模量E=25 GPa,泊松比 ν=0.2.考慮大壩承受水頭超載作用,根據(jù)文獻[36],取上游面超高水頭為 10m,即大壩上游面水深為113 m.模擬中,PD 粒子間距dxS=0.2 m,時間步長ΔtS=5×10-6s ;水體密度 ρF=1000kg/m3,SPH 粒子間距dxF=0.4 m,時間步長ΔtF=5×10-5s,參考聲速c0=1500m/s.為與文獻工況保持一致,防止水流溢出,水體高于壩體部分右側添加了固壁邊界;計算中未考慮上游水滲入裂縫產(chǎn)生的縫面水壓.

圖13 Koyna 重力壩模型(單位:m)Fig.13 Koyna gravity dam model(unit:m)

圖14 給出了Koyna 壩在水頭超載作用下的裂縫擴展路徑.可以看出,水力裂縫首先在初始裂縫頂點萌生,并沿水平方向擴展約3.2 m.隨后,裂縫開始向右下方擴展,此時裂縫與水平方向夾角為33°.本文所得裂縫擴展路徑與文獻[36]結果對比如圖15所示,可以看出,兩種方法所得裂縫路徑基本一致,表明本文提出的多分辨率PD-SPH 混合方法能夠較好地再現(xiàn)流體作用下的結構斷裂過程.

圖14 Koyna 大壩裂縫擴展過程Fig.14 Crack propagation process in Koyna dam

圖15 重力壩裂縫擴展路徑Fig.15 Crack paths in the gravity dam

3.4 水流沖擊作用下的混凝土板崩塌問題

通過上述驗證,本節(jié)將提出的方法應用于工程常見的水流沖擊作用下混凝土板崩塌這一復雜流-固耦合結構破壞問題模擬.初始條件如圖16 所示,水體經(jīng)由寬為w=0.33 m 的水流入口并以5 m/s的速度流入;混凝土板長L=2 m,高H=0.3 m,其左右兩端頂點受到固定約束,中部含一條長0.1 m、寬0.015 m 的豎向裂縫,板的具體材料參數(shù)為:密度ρS=2400kg/m3,彈性模量E=35 GPa,泊松比 ν=0.2.模擬中,令PD 粒子與SPH 粒子的初始間距分別為dxS=0.005 m和dxF=0.01 m,兩相計算時間步長為ΔtS=1×10-6s,ΔtF=5×10-6s ;水體密度ρF=1000kg/m3,參考聲速c0=60m/s.

圖16 水流沖擊混凝土板模型(單位:m)Fig.16 Model of water impacting on a concrete slab(unit:m)

圖17為水流沖擊作用下混凝土板崩塌過程,圖18則給出了混凝土板上表面中部A點(如圖16 所示)的豎向位移時程曲線.在t=0.12 s 時,部分水體自水流出口流出,并在重力作用下加速運動;t=0.212 s 時,水流開始沖擊混凝土板上部;t=0.215 s時,初始裂縫開始擴展,此時A點的豎向位移較小,僅為 - 9.18×10-5m.隨著時間推移,裂縫逐漸向上擴展,并在t=0.222 s 時,擴展至混凝土板上表面形成貫穿裂縫,A點位移也于此刻產(chǎn)生巨大變化.t=0.31 s時,左右兩塊混凝土板在自重及水流沖擊作用下分別繞固定點旋轉,使得裂縫左右表面呈倒V 形,此時A點的豎向位移為 -0.11 m ;水流則因混凝土板的阻擋作用形成橫向射流.最終,水流完全沖開了兩塊獨立運動的混凝土板,并自由向下運動;A點在t=0.45 s 時的豎向位移為 -0.57 m.本文所得混凝土板斷裂過程與文獻[37]一致,表明本文提出的多分辨率PD-SPH 混合方法適用于復雜FSI 問題求解,能夠自然實現(xiàn)流體運動與結構變形破壞全過程連續(xù)仿真.

圖17 水流沖擊作用下混凝土板崩塌過程Fig.17 Collapse process of concrete slab subjected to impact of water flow

圖18 A 點豎向位移時程曲線Fig.18 Time history of vertical displacement at point A

4 小結

本文針對流-固耦合破壞問題提出一種多時間和空間分辨率的PD-SPH 混合方法,并將其應用于典型問題模擬,得到以下結論.

(1)該混合方法利用具有不同初始間距的粒子離散固體域和流體域,采用不同的時間步長對兩相的控制方程進行時間積分.通過將PD 粒子設置為具有與流體粒子相同光滑長度的虛粒子處理流-固界面,可使得界面處滿足動力學和運動學條件,且能夠保證系統(tǒng)動量守恒.

(2)通過液柱靜壓力作用下鋁板變形和潰壩流體沖破彈性閘門兩個經(jīng)典考題模擬,結果表明:多分辨率PD-SPH 混合方法具有較高的計算精度和較好的能量守恒特性,針對不同問題合理選擇流-固分辨率比以及時間步長比可有效提高計算效率.

(3)方法應用于經(jīng)典Koyna 重力壩水力劈裂和水流沖擊作用下混凝土板的崩塌問題,得到了固體結構的運動、變形與斷裂過程以及固體結構影響下流體的運動狀態(tài).所得結果驗證了多分辨率PD-SPH混合方法在分析流-固耦合破壞問題方面的適用性.在此基礎上,后續(xù)可擴展至三維流-固耦合破壞問題研究.本文工作希望為分析工程中同時涉及移動自由表面、復雜流動、流-固耦合導致結構大位移和損傷破壞問題的高效數(shù)值仿真提供參考.