近似不可壓軟材料動力分析的完全拉格朗日物質點法1)

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

(大連理工大學工程力學系,工業裝備結構分析國家重點實驗室,遼寧大連 116024)

引言

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

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

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

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

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

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

1 TLMPM 基本理論

1.1 控制方程

對于大變形問題,其初始構型上坐標為X的點在變形后移動到位置x.可以據此計算變形梯度為

其中,1為2 階單位張量,u=x-X表示位移,?為梯度算子,下標0表示相對于初始構型的量.對于動力學問題,其拉格朗日格式的動量方程為[28]

其中,ρ為密度,ü 表示u對時間的2 階導數,P為PK1 應力,b為體力.引入位移試函數wu并應用高斯散度定理,推導可得等效積分弱形式為

1.2 離散及變量更新

與在計算過程中不斷更新構型的傳統MPM 不同,TLMPM 在求解力學問題時不參考當前構型,而是將初始構型上的求解域離散成若干個物質點,離散后物體的質量、體積等屬性集中在物質點上.同時在求解域范圍內設置背景網格,于是背景網格的等效質量可以由物質點得到,即

其中nI為背景網格節點的總數.類似地,位移、速度和加速度也同樣可以采用類似的方法進行映射.

將插值映射格式帶入弱形式方程(3),經整理可得每個節點的離散控制方程為

其中h0為邊界層的厚度.

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

根據式(10)可以得到各個節點在t+Δt時刻的速度.而對于物質點上的變量,則需要根據節點上的變量進行逆映射,其具體形式為

2 位移-壓強混合格式

2.1 壓強控制方程及其離散

在對凝膠、生物材料等軟材料進行數值求解時,由于其往往具有不可壓特性,因此會產生體積自鎖,對計算結果的精度造成很大影響.因此需要對當前的TLMPM 進行改進,以提升其精度.

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

其中C=FTF,Wdev和Wvol分別為應變能函數的偏量部分和體積部分.這里可以假設軟材料符合Neo-Hookean 本構關系,其偏量應變能密度Wdev的形式為[9]

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

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

在采用TLMPM 求解時為了避免體積自鎖,將應力分解成

其中p為靜水壓力.它與體積應變能存在如下關系

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

于是,物質點上的靜水壓力為

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

2.2 交錯求解流程

為了同時考慮位移和壓強場的求解,混合格式的TLMPM 采用交錯求解的格式,在一個時間步內依次進行位移和壓強場的求解,其流程如下.

(1) 位移場求解:

①將物質點的質量和體積映射到節點上;

②計算節點的內力、外力向量;

③求解位移場方程(7),更新節點速度;

④更新物質點速度、位移、變形梯度;

(2) 壓強場求解:

①將體積變形重映射,計算物質點等效體積比;

②求解壓強場方程(21),更新節點壓強;

③更新物質點的壓強;

(3) 應用本構關系,即式(17)更新應力.

3 數值算例

為了驗證本文提出的針對近似不可壓軟材料動力學的完全拉格朗日物質點算法,本節給出了幾個典型的數值算例.在各個算例的模擬中,為了消除時間步長對結果的影響,時間步長統一設置為

其中Δtcr為臨界時間步長,其確定方法可以參考文獻[17].

3.1 庫克膜問題

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

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

圖2 庫克膜:采用不同數量物質點離散的t=3 s 時右上角豎向位移模擬結果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 時的靜水壓力分布Fig.3 Cook’s membrane:hydrostatic pressure distribution at t=3 s

3.2 二維軟梁的彎曲

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

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

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

分別提取了4 次模擬在t=0.5,1.0s 的靜水壓力分布,如圖6 所示,可以看出線性的標準TLMPM的構型與其他3 種相差較大并且無法計算出合理的靜水壓力,而混合格式 TLMPM 在對近似不可壓材料的體應力計算上可以得到理想的結果.

圖6 二維軟梁的彎曲:不同時刻的靜水壓力分布.(a)～(b) 標準格式線性插值;(c)～(d) 標準格式二次插值;(e)～(f) 線性-常數混合格式;(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 三維柱體扭轉

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

圖7 三維柱體扭轉:構型及靜水壓力分布.(a) 幾何模型;(b) t=0.1 s,標準TLMPM;(c) t=0.3 s,標準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 三維柱體扭轉:A點z 方向位移變化曲線(參考解為文獻[30]的計算結果)Fig.8 Twisting of a 3D column:history curves of z-displacement of point A(the reference results refer to Ref.[30])

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

4 結論

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