章子健 劉振海 張洪武 鄭勇剛

(大連理工大學(xué)工程力學(xué)系,工業(yè)裝備結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗(yàn)室,遼寧大連 116024)

引言

軟材料廣泛存在于自然界和工程中,如水凝膠、形狀記憶聚合物等一些非線性彈性體,以及人體或動(dòng)物的組織和器官等一些生物材料[1].軟材料往往具有高彈性、低模量、耐磨性、抗震性等優(yōu)良特性,并且在經(jīng)歷極大的變形時(shí)沒有能量耗散[2],因此在生活和工程中得到了廣泛的應(yīng)用.對(duì)于這些軟材料力學(xué)行為的分析主要包括理論[3]、實(shí)驗(yàn)[4]和數(shù)值模擬[5]3 種,其中數(shù)值模擬相較理論方法適用范圍更廣,并且不受實(shí)驗(yàn)環(huán)境和設(shè)備的制約[6].因此,針對(duì)軟材料力學(xué)行為分析的數(shù)值方法在近年來得到持續(xù)地發(fā)展和廣泛地應(yīng)用.

有限單元法FEM 作為模擬固體力學(xué)行為最常用的一種數(shù)值方法,被廣泛地運(yùn)用于模擬軟材料的各種復(fù)雜變形[7-8].由于軟材料往往具有不可壓或近似不可壓特性,因此在處理該類材料時(shí)需要對(duì)算法進(jìn)行特殊處理.常見的處理方法包括位移-壓強(qiáng)混合格式、B-bar和F-bar 方法等[9].然而這些方法大多都用于靜力問題的計(jì)算,而對(duì)于不可壓材料的顯式動(dòng)力學(xué)模擬,目前相關(guān)的算法還比較少.此外,由于軟材料在受外界刺激或載荷作用時(shí)可能會(huì)產(chǎn)生極端的大變形[10],采用FEM 進(jìn)行模擬時(shí)可能會(huì)產(chǎn)生網(wǎng)格畸變,導(dǎo)致模擬結(jié)果的誤差較大或者難以收斂,甚至可能模擬失敗[11].因此,FEM 在模擬大變形動(dòng)力學(xué)問題時(shí)仍然具有一定的局限性.

為了更準(zhǔn)確、有效地模擬大變形動(dòng)力學(xué)問題,Sulsky等[12-13]于1994 年將流體力學(xué)中的質(zhì)點(diǎn)網(wǎng)格法加以改進(jìn),提出了一種新的粒子類方法,即物質(zhì)點(diǎn)法(MPM).MPM 將求解對(duì)象離散成物質(zhì)點(diǎn)并設(shè)置獨(dú)立于物質(zhì)點(diǎn)的背景網(wǎng)格,這種處理方式結(jié)合了拉格朗日和歐拉描述的優(yōu)勢(shì),使其在避免FEM 中可能出現(xiàn)的網(wǎng)格畸變的同時(shí)也無需處理歐拉描述中的對(duì)流項(xiàng).因此近年來被廣泛應(yīng)用于高速?zèng)_擊[14]、裂紋擴(kuò)展[15]、邊坡失穩(wěn)[16]等諸多大變形強(qiáng)非線性問題的模擬中.

然而,傳統(tǒng)的MPM 在模擬大變形問題時(shí)物質(zhì)點(diǎn)會(huì)跨越背景網(wǎng)格,其所采用的線性插值函數(shù)的梯度不連續(xù)性會(huì)產(chǎn)生較大的跨網(wǎng)格誤差[17].為了提升MPM 的精度,眾多學(xué)者對(duì)MPM 進(jìn)行了一系列改進(jìn),提出了廣義插值物質(zhì)點(diǎn)法[18]、對(duì)流粒子域插值物質(zhì)點(diǎn)法[19]和B 樣條物質(zhì)點(diǎn)法[20]等方法.de Vaucorbeil等[21]基于B 樣條插值提出了一種完全拉格朗日物質(zhì)點(diǎn)法TLMPM,可以很好地消除這一誤差并減小計(jì)算量.該方法已經(jīng)被成功應(yīng)用于模擬大變形準(zhǔn)靜態(tài)[22]和動(dòng)力沖擊問題[23].

與FEM 類似,MPM 在計(jì)算不可壓材料時(shí)同樣會(huì)產(chǎn)生體積自鎖,也需要進(jìn)行特殊處理.Kularathna等[24]與Zhang等[25]將Chorin 提出的投影方法應(yīng)用到MPM 中,用于求解不可壓流體;Coombs等[26]將F-bar 方法擴(kuò)展到MPM 中用于克服不可壓材料中出現(xiàn)的體積自鎖現(xiàn)象;Iaconeta等[27]提出了一種位移-壓強(qiáng)混合格式的隱式MPM 用于求解不可壓材料的準(zhǔn)靜態(tài)大變形行為.然而目前在求解不可壓軟材料的大變形動(dòng)力學(xué)問題時(shí),顯式MPM 仍然缺乏有效的途徑克服體積自鎖.

本文在TLMPM 理論框架下,基于近似不可壓軟材料的體積應(yīng)變能形式,通過引入壓強(qiáng)控制方程并進(jìn)行離散,提出一種位移-壓強(qiáng)混合格式的顯式TLMPM,用于穩(wěn)定求解近似不可壓軟材料的大變形動(dòng)力學(xué)問題.此外還給出了若干算例,以驗(yàn)證所提出算法的準(zhǔn)確性和有效性.

1 TLMPM 基本理論

1.1 控制方程

對(duì)于大變形問題,其初始構(gòu)型上坐標(biāo)為X的點(diǎn)在變形后移動(dòng)到位置x.可以據(jù)此計(jì)算變形梯度為

其中,1為2 階單位張量,u=x-X表示位移,?為梯度算子,下標(biāo)0表示相對(duì)于初始構(gòu)型的量.對(duì)于動(dòng)力學(xué)問題,其拉格朗日格式的動(dòng)量方程為[28]

其中,ρ為密度,ü 表示u對(duì)時(shí)間的2 階導(dǎo)數(shù),P為PK1 應(yīng)力,b為體力.引入位移試函數(shù)wu并應(yīng)用高斯散度定理,推導(dǎo)可得等效積分弱形式為

1.2 離散及變量更新

與在計(jì)算過程中不斷更新構(gòu)型的傳統(tǒng)MPM 不同,TLMPM 在求解力學(xué)問題時(shí)不參考當(dāng)前構(gòu)型,而是將初始構(gòu)型上的求解域離散成若干個(gè)物質(zhì)點(diǎn),離散后物體的質(zhì)量、體積等屬性集中在物質(zhì)點(diǎn)上.同時(shí)在求解域范圍內(nèi)設(shè)置背景網(wǎng)格,于是背景網(wǎng)格的等效質(zhì)量可以由物質(zhì)點(diǎn)得到,即

其中nI為背景網(wǎng)格節(jié)點(diǎn)的總數(shù).類似地,位移、速度和加速度也同樣可以采用類似的方法進(jìn)行映射.

將插值映射格式帶入弱形式方程(3),經(jīng)整理可得每個(gè)節(jié)點(diǎn)的離散控制方程為

其中h0為邊界層的厚度.

在時(shí)間積分方面,對(duì)于顯式MPM,可以采用向前差分和中心差分兩種格式[21,29].為了簡便計(jì)算,考慮采用向前差分法進(jìn)行積分,便可得到時(shí)間、空間離散后的TLMPM 離散求解格式,即

根據(jù)式(10)可以得到各個(gè)節(jié)點(diǎn)在t+Δt時(shí)刻的速度.而對(duì)于物質(zhì)點(diǎn)上的變量,則需要根據(jù)節(jié)點(diǎn)上的變量進(jìn)行逆映射,其具體形式為

2 位移-壓強(qiáng)混合格式

2.1 壓強(qiáng)控制方程及其離散

在對(duì)凝膠、生物材料等軟材料進(jìn)行數(shù)值求解時(shí),由于其往往具有不可壓特性,因此會(huì)產(chǎn)生體積自鎖,對(duì)計(jì)算結(jié)果的精度造成很大影響.因此需要對(duì)當(dāng)前的TLMPM 進(jìn)行改進(jìn),以提升其精度.

在大變形框架下,近似不可壓軟材料的應(yīng)變能密度可以寫成[30-31]

其中C=FTF,Wdev和Wvol分別為應(yīng)變能函數(shù)的偏量部分和體積部分.這里可以假設(shè)軟材料符合Neo-Hookean 本構(gòu)關(guān)系,其偏量應(yīng)變能密度Wdev的形式為[9]

式中,μ為剪切模量,I1和I3分別為為C的第一和第三不變量.而對(duì)于近似不可壓材料,Wvol可以取

其中,J=det(F),K為體積模量.

在采用TLMPM 求解時(shí)為了避免體積自鎖,將應(yīng)力分解成

其中p為靜水壓力.它與體積應(yīng)變能存在如下關(guān)系

基于TLMPM 框架將該方程離散求解,便可對(duì)體應(yīng)力進(jìn)行修正以避免體積自鎖,從而得到改進(jìn)后的混合TLMPM 求解格式.為了離散方程(18),首先引入壓強(qiáng)試函數(shù)wq以推導(dǎo)出它的弱形式

于是,物質(zhì)點(diǎn)上的靜水壓力為

此外,為了滿足混合格式中的LBB 條件[32]以保證算法的收斂性,需要考慮采用不同階次的B 樣條基函數(shù)分別對(duì)位移和壓強(qiáng)進(jìn)行插值.B 樣條函數(shù)可以較為簡單地實(shí)現(xiàn)高階插值,并且易于改變插值的階次,用其作為插值函數(shù)可以提升MPM 的精度,實(shí)現(xiàn)混合格式的物質(zhì)點(diǎn)法時(shí)無需增加新的節(jié)點(diǎn),程序?qū)崿F(xiàn)較為簡便[33].物質(zhì)點(diǎn)法中具體的B 樣條插值實(shí)施細(xì)節(jié)可以參考文獻(xiàn)[20]的相關(guān)工作.在本文中,如無其他說明,位移和壓強(qiáng)插值分別采用二次和線性B 樣條插值,這種插值方法可以滿足LBB 條件,保證算法的準(zhǔn)確性和收斂性[34].

2.2 交錯(cuò)求解流程

為了同時(shí)考慮位移和壓強(qiáng)場的求解,混合格式的TLMPM 采用交錯(cuò)求解的格式,在一個(gè)時(shí)間步內(nèi)依次進(jìn)行位移和壓強(qiáng)場的求解,其流程如下.

(1) 位移場求解:

①將物質(zhì)點(diǎn)的質(zhì)量和體積映射到節(jié)點(diǎn)上;

②計(jì)算節(jié)點(diǎn)的內(nèi)力、外力向量;

③求解位移場方程(7),更新節(jié)點(diǎn)速度;

④更新物質(zhì)點(diǎn)速度、位移、變形梯度;

(2) 壓強(qiáng)場求解:

①將體積變形重映射,計(jì)算物質(zhì)點(diǎn)等效體積比;

②求解壓強(qiáng)場方程(21),更新節(jié)點(diǎn)壓強(qiáng);

③更新物質(zhì)點(diǎn)的壓強(qiáng);

(3) 應(yīng)用本構(gòu)關(guān)系,即式(17)更新應(yīng)力.

3 數(shù)值算例

為了驗(yàn)證本文提出的針對(duì)近似不可壓軟材料動(dòng)力學(xué)的完全拉格朗日物質(zhì)點(diǎn)算法,本節(jié)給出了幾個(gè)典型的數(shù)值算例.在各個(gè)算例的模擬中,為了消除時(shí)間步長對(duì)結(jié)果的影響,時(shí)間步長統(tǒng)一設(shè)置為

其中Δtcr為臨界時(shí)間步長,其確定方法可以參考文獻(xiàn)[17].

3.1 庫克膜問題

首先考慮庫克膜問題[34],該問題通常被用來驗(yàn)證算法處理體積自鎖的能力.如圖1 所示,模型左端固定,右端受到大小為0.25 N/m 的均勻剪力.模型采用近似不可壓材料,其彈性模量和泊松比分別為1000Pa和0.499,密度為1 kg/m3.為了驗(yàn)證算法的網(wǎng)格收斂性,模擬時(shí)采用不同尺寸的背景網(wǎng)格并將模型離散成不同數(shù)量的物質(zhì)點(diǎn),初始時(shí)每個(gè)背景網(wǎng)格內(nèi)包含4 個(gè)物質(zhì)點(diǎn),采用標(biāo)準(zhǔn)格式和混合格式的TLMPM 分別進(jìn)行模擬,計(jì)算獲得的t=3 s 時(shí)模型右上角點(diǎn)的豎直位移如圖2 所示.可以看出標(biāo)準(zhǔn)的TLMPM 的模擬結(jié)果在物質(zhì)點(diǎn)較密時(shí)逐漸趨近于混合格式的結(jié)果.圖3 給出了兩種算法在t=3 s 時(shí)的靜水壓力分布圖,從圖中可以看出混合格式的 TLMPM相比標(biāo)準(zhǔn)格式能得到更光滑的結(jié)果.這表明了本文提出的混合格式的TLMPM 可以更有效地避免體積自鎖,得到更合理的應(yīng)力結(jié)果.

圖1 庫克膜問題示意圖(單位:mm)Fig.1 Schematic plot of Cook’s membrane(unit:m)

圖2 庫克膜:采用不同數(shù)量物質(zhì)點(diǎn)離散的t=3 s 時(shí)右上角豎向位移模擬結(jié)果Fig.2 Cook’s membrane:simulation results of the vertical displacement of the top right corner at t=3 s using different numbers of material points

圖3 庫克膜:t=3 s 時(shí)的靜水壓力分布Fig.3 Cook’s membrane:hydrostatic pressure distribution at t=3 s

3.2 二維軟梁的彎曲

考慮一個(gè)二維近似不可壓軟材料梁的大變形彎曲問題.模型的尺寸、邊界條件如圖4 所示,梁的下端固定,并整體給定初速度為10m/s.采用近似不可壓Neo-Hookean 材料,對(duì)應(yīng)的彈性模量和泊松比分別17 MPa和0.499,密度為1100kg/m3.將模型采用不同密度的背景網(wǎng)格進(jìn)行離散,保持每個(gè)背景網(wǎng)格內(nèi)包含4 個(gè)物質(zhì)點(diǎn),分別采用混合格式和標(biāo)準(zhǔn)格式的TLMPM 進(jìn)行模擬并考慮不同階次的插值函數(shù),不同網(wǎng)格密度和插值函數(shù)的模擬結(jié)果分別如圖5和圖6 所示.圖5(a)是不同網(wǎng)格尺寸的A點(diǎn)的位移時(shí)程曲線,表明了混合格式的 TLMPM 相比標(biāo)準(zhǔn)TLMPM收斂性更好.對(duì)于網(wǎng)格尺寸為0.1 m 的模型,將插值函數(shù)的階次降低,分別得到線性位移插值的混合格式的 TLMPM和標(biāo)準(zhǔn)TLMPM 模擬結(jié)果與之前的模擬進(jìn)行對(duì)比,其中線性位移插值的混合格式的TLMPM為滿足LBB 條件壓強(qiáng)采用0階插值.結(jié)果如圖5(b)所示,從圖中可以看出混合格式的TLMPM即使降低了插值階次仍然可以得到和二次位移插值較為相似的結(jié)果,而線性插值的標(biāo)準(zhǔn)TLMPM 結(jié)果與其他的相差較大.表明了混合格式的 TLMPM 在進(jìn)行近似不可壓軟材料的大變形問題的模擬時(shí)比標(biāo)準(zhǔn)TLMPM 具有更好的效果.

圖4 二維軟梁的彎曲問題示意圖Fig.4 Schematic plot of bending of 2D soft beam

圖5 二維軟梁的彎曲:A點(diǎn)x 方向位移變化曲線Fig.5 Bending of 2D soft beam:history curves of x-displacement of point A

分別提取了4 次模擬在t=0.5,1.0s 的靜水壓力分布,如圖6 所示,可以看出線性的標(biāo)準(zhǔn)TLMPM的構(gòu)型與其他3 種相差較大并且無法計(jì)算出合理的靜水壓力,而混合格式 TLMPM 在對(duì)近似不可壓材料的體應(yīng)力計(jì)算上可以得到理想的結(jié)果.

圖6 二維軟梁的彎曲:不同時(shí)刻的靜水壓力分布.(a)～(b) 標(biāo)準(zhǔn)格式線性插值;(c)～(d) 標(biāo)準(zhǔn)格式二次插值;(e)～(f) 線性-常數(shù)混合格式;(g)～(h) 二次-線性混合Fig.6 Bending of 2D soft beam:hydrostatic distribution at different times.(a)～(b) Standard linear interpolation;(c)～(d) standard quadratic interpolation;(e)～(f) linear-constant mixed formulation and(g)～(h) quadratic-linear mixed formulation

3.3 三維柱體扭轉(zhuǎn)

最后,考慮一個(gè)三維柱體的扭轉(zhuǎn)問題.一個(gè)方形截面柱底端受約束,整個(gè)柱體受到與其坐標(biāo)相關(guān)的初速度而發(fā)生扭轉(zhuǎn),其構(gòu)型如圖7(a)所示,柱體尺寸為1 m×1 m×6 m,底面固定,將坐標(biāo)原點(diǎn)取在底面中心,給柱體施加與位置坐標(biāo)相關(guān)的初速度為:vx=-ysin(πz/12),vy=xsin(πz/12).材料參數(shù)與3.2 節(jié)的二維軟梁相同.將模型離散成20×20×120個(gè)物質(zhì)點(diǎn),并設(shè)置網(wǎng)格尺寸為0.1 m,即模型內(nèi)部的每個(gè)背景網(wǎng)格內(nèi)有8 個(gè)物質(zhì)點(diǎn).采用標(biāo)準(zhǔn)和混合格式的TLMPM 進(jìn)行模擬,得到0.1 s和0.3 s 的靜水壓力分布如圖7(b)～圖7(e) 所示.從圖中可以看出,混合TLMPM 在模擬三維大變形問題時(shí),靜水壓力結(jié)果在網(wǎng)格內(nèi)和網(wǎng)格之間仍然只有很小的振蕩,而標(biāo)準(zhǔn)的TLMPM 則受體積自鎖影響,無法得到合理的靜水壓力結(jié)果.同時(shí),得到A點(diǎn)的z方向位移時(shí)程曲線如圖8 所示,同時(shí)考慮了混合格式的傳統(tǒng)MPM,并將結(jié)果與改進(jìn)單元技術(shù)的混合格式FEM 結(jié)果[30]進(jìn)行對(duì)比,可以發(fā)現(xiàn)混合格式的TLMPM 與FEM 的結(jié)果相差很小,而標(biāo)準(zhǔn)TLMPM和混合格式的MPM的結(jié)果與它們有一定差距.

圖7 三維柱體扭轉(zhuǎn):構(gòu)型及靜水壓力分布.(a) 幾何模型;(b) t=0.1 s,標(biāo)準(zhǔn)TLMPM;(c) t=0.3 s,標(biāo)準(zhǔn)TLMPM;(d) t=0.1 s,混合TLMPM;(e) t=0.3 s,混合 TLMPMFig.7 Twisting of a 3D column:configurations and distribution of hydrostatic pressure.(a) Geometry;(b) t=0.1 s,standard TLMPM;(c) t=0.3 s,standard TLMPM;(d) t=0.1 s,mixed TLMPM and(e) t=0.3 s,mixed TLMPM

圖8 三維柱體扭轉(zhuǎn):A點(diǎn)z 方向位移變化曲線(參考解為文獻(xiàn)[30]的計(jì)算結(jié)果)Fig.8 Twisting of a 3D column:history curves of z-displacement of point A(the reference results refer to Ref.[30])

通過以上幾個(gè)算例,證明本文所提出的位移-壓強(qiáng)混合格式的TLMPM 算法在模擬近似不可壓軟材料的大變形動(dòng)力學(xué)問題時(shí)可以有效防止體積自鎖,具有很好的收斂性和準(zhǔn)確性.

4 結(jié)論

本文基于MPM 算法框架,發(fā)展了一種位移-壓強(qiáng)混合格式的完全拉格朗日物質(zhì)點(diǎn)法(TLMPM)用于模擬近似不可壓材料軟材料的大變形動(dòng)力學(xué)問題.該方法通過體積應(yīng)變能引入了關(guān)于靜水壓力的控制方程;進(jìn)一步基于MPM 進(jìn)行離散并采用完全拉格朗日格式以避免大變形問題中物質(zhì)點(diǎn)跨網(wǎng)格所產(chǎn)生的誤差;對(duì)位移場和壓強(qiáng)場進(jìn)行交錯(cuò)求解,并采用不同階次的B 樣條基函數(shù)對(duì)位移和壓強(qiáng)分別進(jìn)行插值.除此之外,還針對(duì)體積變形引入重映射技術(shù)以保證靜水壓力求解的準(zhǔn)確性.通過幾個(gè)二維和三維近似不可壓問題的模擬,驗(yàn)證了本文提出的混合格式的TLMPM 可以有效避免體積自鎖,具有很好的收斂性并能夠得到光滑的應(yīng)力結(jié)果,在模擬軟材料大變形動(dòng)力行為時(shí)效果相比標(biāo)準(zhǔn)TLMPM和混合格式的MPM 都更為有效.