擴展鍵基近場動力學的幾何非線性擴展

朱其志，李惟簡，尤 濤，曹亞軍

（1.河海大學巖土力學與堤壩工程教育部重點實驗室，江蘇南京 210098；2.河海大學土木與交通學院，江蘇南京 210098）

傳統(tǒng)的基于經(jīng)典連續(xù)介質(zhì)力學的數(shù)值方法，如有限元法，能夠通過迭代求解線性方程組較好地處理線性連續(xù)問題。然而在處理非連續(xù)問題時，由于難以定義非連續(xù)界面的應力導數(shù)，所以這類方法的數(shù)學模型在裂紋開展處存在一定的缺陷。為了解決這個問題，學者們提出了多種替代性的數(shù)值方法，其中近場動力學方法因其非局部的鍵力描述[1］以及無網(wǎng)格離散的特性[2］而備受關(guān)注。該方法通過對鍵力密度積分的形式重構(gòu)控制方程，從而能夠有效處理位移不連續(xù)問題。通過近些年的發(fā)展，近場動力學方法已經(jīng)具備可靠的彈性預測能力和斷裂預測能力[3-5］，在此基礎(chǔ)上，許多學者開展了大量的研究使其能夠更加精確地捕捉斷裂響應[6-8］。

在原始的鍵基近場動力學公式[1］中，2個相互作用的粒子在法線方向上只有成對的中心力，導致模型的泊松比不能自由定義。為了克服固定泊松比的限制，一些學者在本構(gòu)關(guān)系上進行了相應的擴展[9-13］，其中應用最廣泛的是態(tài)基近場動力學模型。該模型根據(jù)經(jīng)典彈性力學思路將近場動力學變形態(tài)分解成凈體積變化部分以及剩余部分。Madenci等[14］通過應變能密度中偏量部分與等效應力的關(guān)系提出了一種基于態(tài)基近場動力學的塑性模型，但是其相關(guān)研究[15-18］受限于馮米塞斯屈服準則。同時，與常規(guī)態(tài)基近場動力學一樣，其測量體積應變的方式是基于小變形假設(shè)條件，所以在進行幾何非線性分析時需進行特殊處理[19-20］。

Zhu等[21］和Li等[22］開發(fā)的擴展鍵基近場動力學模型（XPD）將剪切應變項引入原始公式，解決了泊松比的限制問題，并且建立了宏觀彈性常數(shù)與微觀鍵模量之間的關(guān)系?；诰植繎兊臄?shù)值實現(xiàn)方法可以準確描述二維、三維任意泊松比的各向同性材料模型的線性彈性行為。在該擴展模型中，通過計算小變形理論下的“局部應變”來確定鍵剪應變效應，其中假設(shè)材料粒子的位移遠小于物體的任何相關(guān)尺寸，因此其幾何形狀和空間各點處材料的本構(gòu)關(guān)系不隨變形而改變。然而，這種近似對于薄壁或細長結(jié)構(gòu)，例如梁、板和殼這類容易受到顯著旋轉(zhuǎn)和位移影響的結(jié)構(gòu)，結(jié)果不可靠。

對于細長結(jié)構(gòu)或薄壁結(jié)構(gòu)，Diyaroglu等[23］引入了近場動力學模型來預測梁的彎曲和橫向剪切變形。O’Grady等[24］引入了非常規(guī)態(tài)基近場動力學模型來分析基于基爾霍夫洛夫理論的彎曲變形平板。Nguyen等[25-26］針對殼結(jié)構(gòu)開發(fā)了一種常規(guī)態(tài)基近場動力學模型，該模型可以用于預測加筋結(jié)構(gòu)的損壞。在這些薄的柔性體中觀察到大位移和大旋轉(zhuǎn)是很常見的，但是Nguyen等[19］提到上述近場動力學（PD）模型都基于線性分析中的小變形假設(shè)。

對于某種數(shù)值方法，將其從小變形預測延伸至大變形分析，是對其應用范圍的重要擴展。它需要考慮結(jié)構(gòu)的幾何非線性特性，這也是眾多工程應用中需要解決的關(guān)鍵問題。因此學者們發(fā)展了多種數(shù)值方法，如無網(wǎng)格方法[27-29］、格子方法[30］、大變形邊界元方法[31-32］等。其中無網(wǎng)格法不依賴網(wǎng)格劃分的質(zhì)量，在處理流體大變形問題中具有獨特的優(yōu)勢，然而其數(shù)值精度較難提高并且施加本質(zhì)邊界條件比較困難；格子方法類似于將材料離散為有限數(shù)量的顆粒，顆粒間通過直接接觸相互作用，雖然其在應對材料斷裂問題時有一定優(yōu)勢，但格子法模擬材料的力學特性嚴重依賴于顆粒的排布方式，從而影響了該方法的普適性。大變形有限元方法[33-36］由于其在連續(xù)變形預測方面的成熟度和可靠性，已被一些基準研究用作參考方法，但有限元方法的網(wǎng)格依賴性以及在處理非連續(xù)問題時的局限性也限制了其應用的范圍。

本文主要考慮幾何非線性影響，旨在將擴展鍵基近場動力學模型應用于大變形問題?？紤]到連續(xù)體的變形構(gòu)型與未變形構(gòu)型存在顯著差異，采用變形梯度張量來表示從初始構(gòu)型到當前構(gòu)型的映射。變形梯度張量可以表示為位移場相對于初始坐標的導數(shù)，因此利用Li等[22］提出的擬合方法計算局部位移函數(shù)，從而得到局部變形梯度場。通過極坐標分解，可以將變形梯度張量分解為2個獨立的運動過程即剛性旋轉(zhuǎn)和變形，從而正確計算出每個擴展鍵的變形部分。

1 擴展鍵基近場動力學

在經(jīng)典近場動力學模型中，局部平衡方程被描述為力密度函數(shù)f(X)的積分，其中屬于彼此近場域的2個物質(zhì)點沿著連接鍵的方向相互作用，如圖1。對于某一物質(zhì)點（所在位置記為X，HX代表物質(zhì)點X的近場域范圍），它在時刻t的運動平衡方程為

圖1 物質(zhì)點分布以及近場域示意Fig.1 Material distribution and definition of horizon

式中：ρ為材料密度；u?為加速度向量；b(X，t)為體力密度；VX′為物質(zhì)點X′所占體積。

取2個位置不同且相互作用的物質(zhì)點，X和X′分別為其位置坐標，則其相對位置向量可以定義為

假設(shè)兩物質(zhì)點在時刻t時分別產(chǎn)生了位移u(X，t)以及u(X′，t)，那么這兩點的相對位移向量η可以寫成

局部鍵本構(gòu)關(guān)系可以描述為

式中：C為鍵剛度張量。文獻[1］在鍵的拉壓變形以外引入了一種局部剪切機制，將鍵剛度張量表示為

式中：ξ=‖ξ‖，n=ξ/ξ表示鍵方向向量；c和κ分別表示鍵在法向以及切向方向的剛度，可以直接由材料的基本參數(shù)求得（楊氏模量E和泊松比ν）。近場域H的大小由δ來表示。

在三維模型中，有

1.1 有限變形理論下的鍵變形計算方法

使用X標記未變形構(gòu)型中某一位置，x(X，t)表示運動變形后構(gòu)型中對應位置，變形梯度張量F是變形后向量的每個分量對于未變形向量的每個分量的導數(shù)

式中：δij為克羅內(nèi)克符號。

極分解定理表明，任何二階張量都可以分解為純旋轉(zhuǎn)和對稱張量的乘積，因此固體的旋轉(zhuǎn)和變形可以分解為

式中：R和S分別為旋轉(zhuǎn)張量以及右伸長張量。

經(jīng)過剛性旋轉(zhuǎn)后的鍵相對位置向量可表示為

由幾何視圖（圖2）可知，有效相對位移為

圖2 鍵位移、旋轉(zhuǎn)、變形示意Fig.2 Bond displacement,rotation and deformation

則鍵的實際變形量可由式（13）、式（14）計算得到：

1.2 變形梯度張量的擬合

局部位移函數(shù)使用局部的物質(zhì)點位移來擬合，此位移函數(shù)可近似為一多元線性函數(shù)u。

式中：a為系數(shù)矩陣的組合；b為剛體位移。從而有

應用最小二乘法最小化每個位移分量的殘差總和，得到

式中：A為系數(shù)a組成的矩陣；Ulocal為局部物質(zhì)點位移組成的矩陣；T為

式中：Xlocal為物質(zhì)點坐標組成的矩陣 。

當?shù)玫阶冃翁荻葟埩縁后，便可以求得R和S。在大變形理論中，S與右柯西-格林張量Q的關(guān)系為

對于二維模型，當?shù)玫絊的特征值λ1和λ2后，可以定義

使用矩陣特征值的特性，λ21和λ22為矩陣S2的特征值。

由于S為正定對稱矩陣，所以λ1＞0，λ2＞0。

1.3 鍵斷裂準則

為了描述材料損傷，標量函數(shù)μ(u，ξ)被引入了力密度-位移函數(shù)關(guān)系f(u，ξ)=μ(u，ξ)C·η。

式中：?c和γc分別為鍵伸長及鍵剪切的臨界值。將能量釋放率G分解為徑向部分以及切向部分，有

式中：θ為鍵與球坐標系z軸的夾角；?為鍵在Oxy平面的投影與x軸的夾角。可得

式中：Gc和Gs分別為Ι型能量釋放率和ΙΙ型能量釋放率。物質(zhì)點x處的損傷值d(x)可按式（28）衡量：

2 數(shù)值求解系統(tǒng)

利用動態(tài)松弛法模擬結(jié)構(gòu)的準靜態(tài)力學行為。對于整個系統(tǒng)，引入自適應阻尼系數(shù)cn結(jié)合虛擬對角質(zhì)量矩陣M[38］，建立加速度向量U?、速度向量U?和力向量F之間的關(guān)系，如式（1）：

利用中心差分顯式積分，相鄰迭代步間的位移關(guān)系可以寫成

其中初始迭代步的速度向量為

隨后迭代步的速度向量通過式（32）計算：

對角質(zhì)量矩陣M中的主對角元素必須滿足

式中：Kij可采用文獻[22］方法組成的剛度矩陣。每一迭代步的自適應阻尼系數(shù)取決于

式中：L為對角局部剛度矩陣，通過式（35）計算：

需注意，如果遇到速度為零的情況，Lii便設(shè)為零。

3 算例驗證

進行相關(guān)基準測試以驗證所提出的XPD模型在大變形分析中的可靠性。3.1節(jié)分析承受集中荷載的懸臂梁，這是一個經(jīng)典的基準問題。3.2節(jié)引入了一個L形板來測試XPD大變形預測的有效性。最后將大變形XPD結(jié)果與大變形有限元結(jié)果進行對比，以驗證所建立模型的準確性。假設(shè)模型為平面應力情況，單位厚度。

3.1 承受集中力的懸臂梁

懸臂梁的幾何形狀、邊界條件和加載條件如圖3所示。其中，楊氏模量為1 000MPa，泊松比為零。對于XPD離散，近場域δ=0.175m，網(wǎng)格間距Δx=0.05m（與有限元分析的網(wǎng)格尺寸相同）。為了在左邊緣施加邊界條件，將物質(zhì)點的虛擬層添加到離散化模型中（見圖4左端面外側(cè)物質(zhì)點，寬度為δ，高度為1m），并且這些虛擬點的所有自由度都設(shè)置為零。

圖3 懸臂梁的幾何尺寸與邊界條件（單位：m）Fig.3 Cantilever beam:geometry and boundary conditions(unit:m)

圖4 物質(zhì)點分布與變形示意Fig.4 Material point distribution and its deformation

圖5、圖6展示了線性XPD、非線性XPD和非線性FEM預測的位移分布。通過對比可以看出，圖5結(jié)果與圖6結(jié)果有顯著差異，因為在小變形理論中，位移隨著外力的增加呈線性增長，即計算僅依賴模型的初始位置，而在大變形理論中，每對物質(zhì)點之間的相互作用力由模型剛體位移轉(zhuǎn)動后的變形確定。由圖6可知，大變形XPD預測與大變形有限元預測結(jié)果吻合較好，2個方向的位移分量都顯示出非常好的一致性。圖7呈現(xiàn)了整個加載過程中記錄的撓度荷載信息，表明大變形XPD方法與大變形有限元方法的加載路徑幾乎相同。

圖5 使用小變形理論預測的x2方向位移云圖Fig.5 Contour maps predicted by using small deformation theory

圖6 使用大變形理論預測的x2方向位移云圖Fig.6 Contour maps predicted by using large deformation theory

圖7 懸臂梁的撓度-荷載分析Fig.7 Cantilever beam:deflection-load analysis

3.2 承受豎向面力的L型板

考慮幾何非線性和材料損傷行為的共同影響。L形板在其右上側(cè)面承受豎直載荷，幾何尺寸和邊界條件如圖8所示。楊氏模量設(shè)置為1 000MPa、泊松比設(shè)置為ν=0.1。對于XPD離散化，δ=0.35m且Δx=0.1m。

圖8 L形板的幾何尺寸與邊界條件（單位：m）Fig.8 L-shaped plate:geometry and boundary conditions(unit:m)

記錄F=70 000kN時的位移分量的云圖，圖9和圖10的結(jié)果表明，這2種方法在求解結(jié)果上產(chǎn)生的差異很小。

圖9 L形板的使用大變形理論預測的x2方向位移云圖Fig.9 L-shaped plate:displacement distribution predicted by using large deformation theory

圖10 L形板的上表面位移Fig.10 L-shaped plate:x2-direction displacement on upper surface

圖11描繪了L形板上表面在2個不同加載階段的x2方向的位移。圖12記錄了板上表面3個不同點的位移-荷載曲線。從這些比較中可以看出，大變形XPD預測與大變形FEM的解幾乎一致。

圖11 L形板上3個指定點的位移-荷載曲線Fig.11 L-shaped plate:x2-direction displacementload curves of three specified points

圖12是大變形XPD的斷裂模擬。將臨界能量釋放率設(shè)為456 J·mm-2時，損傷在F=74 000kN處開始演化。如果檢測到任何鍵斷裂，加載將不會繼續(xù)，而是等待系統(tǒng)迭代收斂再檢測斷裂，直至沒有斷裂發(fā)生才會繼續(xù)加載。在大約F=131 000kN時，載荷不會增加，而L形板將繼續(xù)向左側(cè)尖端撕裂。圖12中呈現(xiàn)的變形模式與文獻[39］中的結(jié)果非常吻合。

圖12 L形板的大變形XPD預測的損傷模式Fig.12 L-shaped plate:damage patterns predicted by large deformation XPD

4 結(jié)語

基于擴展鍵基近場動力學模型提出了一種數(shù)值實現(xiàn)方法。通過應用變形梯度張量在新的求解框架中很好地解決了幾何非線性問題。在幾組基準測試研究中，通過新框架模擬結(jié)果與大變形有限元結(jié)果之間的對比驗證了模型的有效性和準確性。同時，該模型還能很好地模擬伴隨大變形的材料非線性行為，包括裂紋擴展等。與小變形模型相比，除了每一步需要更新粒子的位置外，并沒有增加額外的計算消耗。鑒于其預測變形的能力以及編碼簡單的優(yōu)勢，擴展鍵基近場動力學模型有望成為解決工程實踐中復雜問題的有效工具。

作者貢獻聲明：

朱其志：提供研究平臺以及理論基礎(chǔ)，進行研究指導以及寫作指導。

李惟簡：進行方法創(chuàng)新以及程序?qū)崿F(xiàn)。

尤 濤：提供編程幫助。

曹亞軍：提供理論及創(chuàng)新思路。