流固耦合破壞分析的多分辨率PD-SPH 方法1)

姚學(xué)昊 陳丁 武立偉 黃丹

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

引言

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

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

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

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

1 數(shù)值模擬方法

1.1 PD 方法

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

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

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

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

1.2 SPH 方法

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

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

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

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

對(duì)于不可壓縮流體,粒子xi受到的物理黏性力可簡(jiǎn)化為

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

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

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

2.1 多分辨率空間離散

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

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

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

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

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

2.2 多分辨率時(shí)間離散

考慮到PD 與SPH 的完全顯式計(jì)算,分別定義結(jié)構(gòu)和流體的計(jì)算時(shí)間步ΔtS和ΔtF,其中,固體計(jì)算的時(shí)間步長(zhǎng)應(yīng)滿足[31]

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

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

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

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

3 算例驗(yàn)證

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

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

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

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

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

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

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

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

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

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

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

圖8 不同工況的計(jì)算時(shí)間Fig.8 Computational times for different cases

3.2 潰壩流體沖擊彈性閘門

Antoci等[35]對(duì)潰壩水流沖破彈性閘門問(wèn)題進(jìn)行了實(shí)驗(yàn)研究,該問(wèn)題的初始條件示意圖如圖9 所示.實(shí)驗(yàn)在寬L=0.1 m、內(nèi)部水深H1=0.14 m 的水箱中進(jìn)行.箱體左側(cè)剛性壁面下部放置了上端固定下端自由的彈性橡膠閘門,其寬度為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,總模擬時(shí)長(zhǎng)0.4 s.針對(duì)該問(wèn)題的多分辨率工況設(shè)置如表2 所示.

表2 潰壩流體沖破彈性閘門問(wèn)題工況布置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 中不同時(shí)刻的彈性閘門在時(shí)變水壓作用下的變形以及自由液面變化過(guò)程,并與實(shí)驗(yàn)觀察結(jié)果[35]進(jìn)行比較.由圖10可以看出,PD-SPH 所得彈性閘門與自由液面形狀與實(shí)驗(yàn)結(jié)果基本吻合.在t=0.08s 時(shí),閘門自由端在水壓力作用下向x軸負(fù)方向移動(dòng),形成水流出口.在t=0.16 s 時(shí),閘門自由端位移與水流出口顯著增大.隨后,由于水位的降低,水壓力減小,水流出口逐漸減小.值得注意的是,由于本文開展的是二維模擬,實(shí)驗(yàn)中的流體飛濺現(xiàn)象不會(huì)出現(xiàn).

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

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

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

圖12 3 種工況的計(jì)算總時(shí)長(zhǎng)Fig.12 Total calculation time in three cases

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

3.3 重力壩水力劈裂

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

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

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

圖14 Koyna 大壩裂縫擴(kuò)展過(guò)程Fig.14 Crack propagation process in Koyna dam

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

3.4 水流沖擊作用下的混凝土板崩塌問(wèn)題

通過(guò)上述驗(yàn)證,本節(jié)將提出的方法應(yīng)用于工程常見的水流沖擊作用下混凝土板崩塌這一復(fù)雜流-固耦合結(jié)構(gòu)破壞問(wèn)題模擬.初始條件如圖16 所示,水體經(jīng)由寬為w=0.33 m 的水流入口并以5 m/s的速度流入;混凝土板長(zhǎng)L=2 m,高H=0.3 m,其左右兩端頂點(diǎn)受到固定約束,中部含一條長(zhǎng)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,兩相計(jì)算時(shí)間步長(zhǎng)為Δ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為水流沖擊作用下混凝土板崩塌過(guò)程,圖18則給出了混凝土板上表面中部A點(diǎn)(如圖16 所示)的豎向位移時(shí)程曲線.在t=0.12 s 時(shí),部分水體自水流出口流出,并在重力作用下加速運(yùn)動(dòng);t=0.212 s 時(shí),水流開始沖擊混凝土板上部;t=0.215 s時(shí),初始裂縫開始擴(kuò)展,此時(shí)A點(diǎn)的豎向位移較小,僅為 - 9.18×10-5m.隨著時(shí)間推移,裂縫逐漸向上擴(kuò)展,并在t=0.222 s 時(shí),擴(kuò)展至混凝土板上表面形成貫穿裂縫,A點(diǎn)位移也于此刻產(chǎn)生巨大變化.t=0.31 s時(shí),左右兩塊混凝土板在自重及水流沖擊作用下分別繞固定點(diǎn)旋轉(zhuǎn),使得裂縫左右表面呈倒V 形,此時(shí)A點(diǎn)的豎向位移為 -0.11 m ;水流則因混凝土板的阻擋作用形成橫向射流.最終,水流完全沖開了兩塊獨(dú)立運(yùn)動(dòng)的混凝土板,并自由向下運(yùn)動(dòng);A點(diǎn)在t=0.45 s 時(shí)的豎向位移為 -0.57 m.本文所得混凝土板斷裂過(guò)程與文獻(xiàn)[37]一致,表明本文提出的多分辨率PD-SPH 混合方法適用于復(fù)雜FSI 問(wèn)題求解,能夠自然實(shí)現(xiàn)流體運(yùn)動(dòng)與結(jié)構(gòu)變形破壞全過(guò)程連續(xù)仿真.

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

圖18 A 點(diǎn)豎向位移時(shí)程曲線Fig.18 Time history of vertical displacement at point A

4 小結(jié)

本文針對(duì)流-固耦合破壞問(wèn)題提出一種多時(shí)間和空間分辨率的PD-SPH 混合方法,并將其應(yīng)用于典型問(wèn)題模擬,得到以下結(jié)論.

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

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

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