冰彎曲破壞的彈塑性近場動力學(xué)模型

張 媛，王 超，郭春雨，葉禮裕，劉 正

(哈爾濱工程大學(xué) 船舶工程學(xué)院，哈爾濱 150001)

近些年來，在氣候變暖、科學(xué)技術(shù)發(fā)展、經(jīng)濟(jì)效益驅(qū)動以及一些政治和軍事領(lǐng)域的顯著價(jià)值需求下，極地航線開發(fā)、極地資源的勘探以及極地環(huán)境利用成為一種必然且迫切的趨勢[1]。這些變化極大地促進(jìn)了更加頻繁的極地地區(qū)海洋結(jié)構(gòu)物的作業(yè)活動，除此之外，多數(shù)內(nèi)陸結(jié)冰水域仍舊常年進(jìn)行各種海上作業(yè)活動。然而冰的存在能夠?qū)Y(jié)構(gòu)物造成嚴(yán)重的損害，甚至造成結(jié)構(gòu)物的破壞而無法正常作業(yè)。由此研究冰力學(xué)性能和破壞模式及其與結(jié)構(gòu)物相互作用過程是支持極地以及冰區(qū)結(jié)構(gòu)物發(fā)展的重要手段之一，而建立精確的冰本構(gòu)模型是這其中的關(guān)鍵技術(shù)之一。

基于數(shù)值計(jì)算的冰力學(xué)特性、冰-結(jié)構(gòu)物作用過程等研究近幾年快速發(fā)展：最初，學(xué)者們通常采用基于試驗(yàn)獲得的經(jīng)驗(yàn)數(shù)據(jù)和數(shù)學(xué)物理分析的方法來計(jì)算冰的力學(xué)特性并且用于冰與結(jié)構(gòu)物相互作用過程的冰載荷預(yù)報(bào)以及船舶運(yùn)動性能預(yù)報(bào)[2-4]。然而上述文獻(xiàn)中對于海冰性質(zhì)研究的方法，極大地依賴于經(jīng)驗(yàn)數(shù)據(jù)，不足以發(fā)展成完整成熟的數(shù)值模型[5]。后來，學(xué)者們紛紛采用基于傳統(tǒng)連續(xù)介質(zhì)力學(xué)方法的數(shù)值手段來研究海冰的特性，例如有限元方法FEM[6-8]，改進(jìn)有限元方法[9-10]，以及其他考慮冰細(xì)微參數(shù)的方法等[11]。然而傳統(tǒng)介質(zhì)力學(xué)方法在處理大變形或者斷裂問題具有一定的局限性，不能準(zhǔn)確地模擬冰破壞過程的裂紋擴(kuò)展路徑。為了克服上述的科學(xué)問題，學(xué)者開始嘗試建立基于無網(wǎng)格方法的冰本構(gòu)模型，這些方法里發(fā)展比較成型的有：離散元方法(Discrete element method,DEM)[12-13]，光滑粒子動力學(xué)(Smoothed particle hydrodynamics,SPH)[14-16]以及近場動力學(xué)方法(Peridynamics,PD)[17-19]等。其中DEM方法由于其顆粒狀的模型形式具有模擬浮碎冰和結(jié)構(gòu)物碰撞、浮碎冰之間碰撞的優(yōu)勢。而SPH和PD方法的粒子離散形式使其在海冰破碎以及裂紋擴(kuò)展方面更具有更突出的優(yōu)勢，不管是哪種粒子數(shù)值模型，雖然能夠處理斷裂問題，但是并不能完整精確地模擬冰力學(xué)性能和破壞特性，由此這些方法仍然處于發(fā)展前期，還需要投入大量的研究工作來建立適用于冰的準(zhǔn)確數(shù)值模型。

近場動力學(xué)方法(PD)，是近二十年新興的一種非局部、無網(wǎng)格粒子、積分形式的計(jì)算方法，在斷裂問題域和非斷裂問題域上都可以連續(xù)求解[20]。該方法最早由Silling等[21-22]提出并建立推導(dǎo)了其控制方程，給出了基本的數(shù)值實(shí)現(xiàn)程序案例。PD方法主要兩種表現(xiàn)形式[21-23]，經(jīng)推導(dǎo)和比較[24]兩種表達(dá)方式大同小異且具有相類似的計(jì)算精度。PD發(fā)展至今主要有鍵型、常規(guī)狀態(tài)型和非常規(guī)狀態(tài)型3種形式。且該方法已經(jīng)成功應(yīng)用于材料的彈性、塑性、熱力學(xué)等多種物理問題計(jì)算分析領(lǐng)域[25-27]。然而PD應(yīng)用于冰計(jì)算領(lǐng)域剛剛起步，多數(shù)文獻(xiàn)采用基于鍵型的彈性模型來模擬冰的破壞特性，例如，冰與剛體柱的作用過程、冰槳作用過程等，忽略了泊松比的真實(shí)狀態(tài)和海冰具有塑性屈服過程的力學(xué)特性。由此本文主要的目的是建立基于狀態(tài)型PD方法的冰彈塑性本構(gòu)模型，消除了鍵型對泊松比的限制，加入了冰破壞過程塑性變形過程，為后續(xù)的冰與結(jié)構(gòu)物的作用過程提供支撐。

1 基于常規(guī)近場動力學(xué)的彈塑性本構(gòu)模型

1.1 基于常規(guī)近場動力學(xué)的彈塑性理論

1.1.1 近場動力學(xué)理論

近場動力學(xué)方法將連續(xù)介質(zhì)離散為均勻的物質(zhì)點(diǎn)。在參考坐標(biāo)系下，其運(yùn)動方程為

t′(u-u′,x′-x,))dH+b(x,t)

(1)

上式的離散形式為

t′(u(k)-u(j),x(k)-x(j),t))dV(j)…+b(k)

(2)

式中，每個(gè)粒子的信息用其位置坐標(biāo)表示，即x(k)，該粒子占據(jù)的空間大小為V(k)且密度為ρ(k)。在笛卡爾坐標(biāo)中，材料點(diǎn)受到外界作用的位移為u(k)，此時(shí)的位置向量為y(k)。由此，材料點(diǎn)k在未變形狀態(tài)下的位置為x(k)，其位移和外力向量分別用u(k)(x(k),t)和b(k)(x(k),t)表示，對k點(diǎn)有作用的域(近場域)表示為Hx(k)。

1.1.2 塑性理論

本文將海冰視作各向同性的均勻介質(zhì)材料。下述本構(gòu)模型的建立基于該假設(shè)成立。根據(jù)經(jīng)典連續(xù)介質(zhì)力學(xué)中的塑性定義，當(dāng)材料在加載/卸載條件下的應(yīng)力水平超過屈服應(yīng)力時(shí)，材料在完全卸載時(shí)會發(fā)生永久(或塑性)變形。與彈性變形不同，塑性變形取決于加載歷史。因此，塑性變形更加關(guān)注每一時(shí)間步力密度增量和拉伸增量的關(guān)系。用PD進(jìn)行塑性分析時(shí)，需要對等效應(yīng)力進(jìn)行描述，這是因?yàn)樵谒苄苑治鲋?，彈性區(qū)域受屈服面限制。此外，屈服面在每個(gè)加載步驟中的增長和移動都需要確定，因此需要將等效塑性拉伸描述為硬化內(nèi)變量。

為了用參量的增量形式表達(dá)PD中的塑性變形，將外界作用下兩個(gè)材料點(diǎn)之間的拉伸增量Δs(k)(j)以及力密度增量Δt(k)(j)分解成彈性項(xiàng)和塑性項(xiàng)如下:

(3)

(4)

本文的塑性理論中假設(shè)塑性變形與體積變形無關(guān)，由此塑性拉伸的過程體積膨脹為0，即Δθp=0。PD方法中靜水壓力分別在力密度函數(shù)和應(yīng)變能密度的彈性項(xiàng)中，用p=-kθ表述，其中k為冰的體積模量。本文模型中塑性變形與靜水壓力無關(guān)。同樣地，每一步的應(yīng)變能密度增量分解為彈性應(yīng)變能密度增量和塑性應(yīng)變能密度增量如下:

(5)

1.1.3 屈服函數(shù)

傳統(tǒng)的塑性理論中，塑性變形由屈服面來判定。參照傳統(tǒng)塑性理論屈服面定義為[28]

(6)

(7)

依據(jù)Von-Mises 塑性屈服準(zhǔn)則，塑性屈服發(fā)生在偏斜應(yīng)變部分的應(yīng)變能密度達(dá)到極限值的時(shí)候，對于單軸壓縮來說，這個(gè)極限值等于材料達(dá)到初始屈服應(yīng)力σ0時(shí)的應(yīng)變能密度，即

(8)

式中J2為應(yīng)力偏張量的第二不變量。為了在一般荷載情況下使用該方程，需要定義等效應(yīng)力。等效應(yīng)力可定義為[29]

(9)

1.1.4 流動法則和硬化規(guī)律

(10)

式中，C(k)為一個(gè)正比例常數(shù)，上式作為常態(tài)條件表述了塑性伸長增量的流動法則。

在式(7)中，理想塑性(無硬化)狀態(tài)下，所有加載/卸載過程中，σy是恒定的。然而，塑性變形在不同的狀態(tài)下與等效塑性拉伸具有相關(guān)性，這種相關(guān)性有多種情況，例如各向同性硬化關(guān)系、運(yùn)動硬化和混合硬化關(guān)系等[28]。本文采用線性各向同性硬化模型[29]，即式(6)可改寫為

(11)

其中

在簡單的單軸拉伸情況下，無膨脹項(xiàng)，可以確定塑性增量拉伸引起的畸變應(yīng)變能密度，并由此提取參數(shù)A0。參數(shù)A0已在文獻(xiàn)[30]中針對不同的物體維度推導(dǎo)得出，本文未列出。

1.2 邊界條件

由于PD運(yùn)動方程不包含任何空間導(dǎo)數(shù)，因此一般不需要約束條件來求解積分微分方程。與局部理論不同，邊界條件是通過非零體積的虛擬邊界層施加的?；跀?shù)值實(shí)驗(yàn)，Macek等[30]建議虛擬邊界層的范圍等于領(lǐng)域尺寸δ，以確保施加的規(guī)定約束在真實(shí)區(qū)域中得到準(zhǔn)確反映。位移邊界條件可以通過對虛擬層中的材料點(diǎn)進(jìn)行約束來施加，使邊界曲面上的條件得到明確滿足。因此，虛擬層中的位移值可以基于實(shí)域值和邊界條件指定值的線性外推來近似。類似地，也可以通過近似虛擬區(qū)域中的牽引力值來施加牽引力邊界條件，從而使真實(shí)區(qū)域和虛擬層中的牽引力變化恢復(fù)邊界表面上施加的牽引力。圖1表示了位移邊界條件的施加方式，圖中Rf表示虛擬粒子邊界，寬度為δ，R為真實(shí)粒子區(qū)域，位移和速度邊界條件為U*和V*。則邊界條件可以表達(dá)為：

圖1 位移和速度約束條件[31]Fig.1 Boundary condition of displacement and velocity[31]

uf(xf,yf,t+Δt)=2U*(x*,y*,t+Δt)-u(x,y,t)

vf(xf,yf,t+Δt)=2V*(x*,y*,t+Δt)-v(x,y,t)

(12)

1.3 失效準(zhǔn)則

基于PD方法的彈性變形中，失效的判斷標(biāo)準(zhǔn)為極限伸長量，當(dāng)變形后粒子間的伸長量達(dá)到極限伸長量s0時(shí)，則判斷為失效，此過程不可逆。極限伸長量可以通過材料的極限能量釋放率得到，粒子間的伸長量和力密度的關(guān)系呈線性。然而在塑性變形中，粒子間力密度和伸長量的關(guān)系呈典型的非線性，如圖2所示，因此彈塑性本構(gòu)模型中的失效準(zhǔn)則不能再簡單地使用最終狀態(tài)的伸長量來判斷，必須對整個(gè)變形過程伸長量對應(yīng)的力密度積分來判斷材料的失效。圖2曲線下的面積代表彈塑性變形微勢w(k)(j)，該函數(shù)代表變形的程度，可由下式得到:

圖2 塑性理論中的變形本構(gòu)關(guān)系Fig.2 Constitutive relationship in plastic deformation

(13)

如果一個(gè)新的裂紋產(chǎn)生，則穿越這個(gè)裂紋的所有粒子間的作用都必須消失。同理，消除兩個(gè)粒子間作用所需的極限能量釋放率作為判別粒子間作用是否失效的標(biāo)準(zhǔn)，則可以表達(dá)為

(14)

(15)

(16)

式中：Gc為材料的極限能量釋放率，Nc為PD的對應(yīng)參數(shù)，在一維時(shí)取值為12，二維時(shí)取值為36，三維時(shí)取值為560。

粒子作用的失效通過上述準(zhǔn)則判斷之后，失效的粒子作用之間力將被移除。此時(shí)，通過監(jiān)測粒子之間相互作用的失效與否就可以表述材料的破壞過程。因此引入局部破壞系數(shù)φ(x(k),t)來表述材料在局部位置的破壞程度，該參數(shù)的物理意義是：材料點(diǎn)k在其作用域內(nèi)失效粒子作用與總粒子作用的比值，由下式計(jì)算：

(17)

2 數(shù)值執(zhí)行過程

2.1 數(shù)值執(zhí)行策略

PD方法的控制方程主要采取積分形式，因此在數(shù)值實(shí)現(xiàn)時(shí)，計(jì)算域被分解成粒子形式。計(jì)算每個(gè)粒子在每個(gè)時(shí)間步長的物理信息，作為下一時(shí)間步的基礎(chǔ)值，每個(gè)時(shí)間步下對所有粒子的物理信息求和積分則為該步長的變形狀態(tài)。

本文參考了文獻(xiàn)[24]中給出的鍵型PD的程序基礎(chǔ)框架。本文的數(shù)值執(zhí)行過程在上述程序框架進(jìn)行補(bǔ)充完善。首先，常規(guī)型PD積分過程中增加了體積膨脹項(xiàng)，可按下式數(shù)值執(zhí)行:

(18)

時(shí)間積分內(nèi)的數(shù)值執(zhí)行過程將全部替換為塑性變形的數(shù)值過程，可參考文獻(xiàn)[28-29]中的數(shù)值實(shí)現(xiàn)方法。對于失效準(zhǔn)則的積分過程，式(13)將采用如下積分過程實(shí)現(xiàn):

t+Δtw(k)(j)=tw(k)(j)+Δw(k)(j)

(19)

(20)

局部破壞系數(shù)的數(shù)值執(zhí)行形式為

(21)

數(shù)值框架流程如圖3所示。

圖3 數(shù)值實(shí)現(xiàn)程序框圖Fig.3 Flowchart of numerical implementation

2.2 數(shù)值方法驗(yàn)證

為了驗(yàn)證文章模型的準(zhǔn)確性和程序的可靠性，首先計(jì)算了帶孔二維板的變形問題。二維板的尺寸選取長和寬分別為L=H=1.0 m，板的厚度為h=0.01 m，孔的尺寸為D=0.3 m。彈性模量為E=200 MPa，密度為ρ=4 428 kg/m3，泊松比為ν=1/3，剛度模量為K=50 GPa。離散狀態(tài)的粒子間距為dx=0.002 5 m。同時(shí)采用有限元分析軟件Abaqus計(jì)算了板的有限元變形，有限元分析中，網(wǎng)格大小與粒子間距相同，網(wǎng)格數(shù)有8 910個(gè)。板的左右兩側(cè)采用位移約束，其位移變形約束表示為：ux(x=0,y,t)=-0.001 m和uy(x=L,y,t)=0.001 m。板的工況示意圖如圖4所示。

圖4 二維板拉伸變形示意圖Fig.4 Sketch map of 2D plate under tension

正如預(yù)期的那樣，塑性變形在高應(yīng)力集中區(qū)域開始。圖5中給出了帶孔板的等效應(yīng)力和位移的變化，與有限元結(jié)果對比完全一致，該案例驗(yàn)證了本文程序的準(zhǔn)確性。

圖5 有孔板拉伸狀態(tài)下的等效應(yīng)力、x方向位移以及y方向位移Fig.5 Variations of effective stress, displacement in x direction, and displacement in y direction of plate with a hole under tension

3 冰四點(diǎn)彎曲數(shù)值計(jì)算

冰的彎曲強(qiáng)度是計(jì)算冰-船作用的關(guān)鍵參數(shù)之一，四點(diǎn)彎曲試驗(yàn)可以用來測量冰的彎曲強(qiáng)度。為了驗(yàn)證本文冰模型的準(zhǔn)確性，選用Ehlers和Kujala進(jìn)行的試驗(yàn)[32]作為對比驗(yàn)證，該試驗(yàn)?zāi)Ｐ图凹虞d工況示意如圖6所示。

圖6 冰四點(diǎn)彎曲試驗(yàn)示意圖Fig.6 Schematic diagram of four-point bending test of ice

該試驗(yàn)中，冰梁的尺寸為L=4.25 m,H=0.4 m和B=0.365 m。上方的固定支撐位于中點(diǎn)向兩邊2 m處，底部位移壓頭位于0.5 m處。上方的支撐在空間中固定，下方的支撐以Vsupport=0.003 580 m/s勻速向上移動。海冰參數(shù)的選取參考試驗(yàn)和文獻(xiàn)中對冰力學(xué)性能相關(guān)分析設(shè)定[32-34]，冰材料的極限能量釋放率采用了文獻(xiàn)[35]中的參考值，見表1。

表1 試驗(yàn)工況設(shè)定和海冰參數(shù)Tab.1 Test conditions and parameters of ice

文中采用2D進(jìn)行計(jì)算，因此彈性模量轉(zhuǎn)換成2D參數(shù)為E/(1-v2)[34]。結(jié)果如圖7、8所示。冰的彎曲載荷和時(shí)間的關(guān)系曲線如圖9所示。

從計(jì)算結(jié)果可以發(fā)現(xiàn)，冰梁在下部位移支撐和上層約束支撐的共同作用，發(fā)生微小變形，如圖7(a)、(b)所示，兩個(gè)移動支撐之間的冰梁位移最大，且這段部分各個(gè)位置的位移量基本相同。冰梁上越遠(yuǎn)離移動支撐部位的縱向位移變形越小，這樣的變形一直持續(xù)到0.31 s之前。此時(shí)海冰未發(fā)生破壞，由于固定支撐和加載支撐的約束原因，冰x方向和y方向的位移符合物理規(guī)律。還可以通過等效應(yīng)力分布圖觀測到：高應(yīng)力分布在拉伸和壓縮集中的部位，例如冰梁的上側(cè)邊緣和下側(cè)邊緣，如圖5(c)所示，這兩側(cè)的位置是冰梁擠壓和拉伸最為嚴(yán)重的位置。該應(yīng)力現(xiàn)象與文獻(xiàn)[35]中通過有限元方法計(jì)算的結(jié)果具有很高的一致性。從圖8可以觀測到：冰梁彎曲后分裂成3部分，斷裂位置對應(yīng)于底部壓頭加載位置。該現(xiàn)象和試驗(yàn)現(xiàn)象具有極高的一致性。另外為了進(jìn)一步驗(yàn)證計(jì)算結(jié)果的準(zhǔn)確性，對比數(shù)值計(jì)算和試驗(yàn)中的力-時(shí)間曲線(圖9)可知：數(shù)值計(jì)算結(jié)果和試驗(yàn)結(jié)果吻合度極高。

圖7 t=0.31 s,Vsupport=0.003 580 m/s(未發(fā)生破壞之前)冰四點(diǎn)彎曲試驗(yàn)數(shù)值計(jì)算結(jié)果Fig.7 Numerical results of four-point bending test of ice (before damage) (t=0.31 s,Vsupport=0.003 580 m/s)

圖8 t=1.00 s,Vsupport=0.003 580 m/s(破壞之后)冰四點(diǎn)彎曲試驗(yàn)數(shù)值計(jì)算結(jié)果Fig.8 Numerical results of four-point bending test of ice (after damage) (t=1.00 s,Vsupport=0.003 580 m/s)

圖9 Vsupport=0.003 580 m/s冰四點(diǎn)彎曲載荷-時(shí)間曲線圖Fig.9 Force-time curve of four-point bending test of ice (Vsupport=0.003 580 m/s)

為了驗(yàn)證本文的數(shù)值模型能夠模擬不同時(shí)間工況下冰梁的彎曲破壞，分別選取了Vsupport=0.002 750 m/s和Vsupport=0.003 147 m/s加載工況進(jìn)行對比驗(yàn)證，載荷-時(shí)間變化關(guān)系計(jì)算結(jié)果如圖10所示，數(shù)值結(jié)果和試驗(yàn)結(jié)果基本保持一致。由試驗(yàn)結(jié)果可知：加載速度為Vsupport=0.003 580 m/s時(shí)的彎曲最大載荷最大，斷裂時(shí)間最短。表2列出了3種加載速度下數(shù)值計(jì)算和試驗(yàn)對應(yīng)的極值載荷數(shù)值極其發(fā)生時(shí)間，從表中可以看出試驗(yàn)的結(jié)果和計(jì)算值的誤差在10%以內(nèi)，說明本文建立的無網(wǎng)格法冰本構(gòu)模型能夠模擬冰的彎曲破壞過程。

圖10 冰四點(diǎn)彎曲載荷-時(shí)間曲線圖Fig.10 Force-time curve of four-point bending test of ice

表2 數(shù)值計(jì)算和試驗(yàn)中3種速度工況對應(yīng)的計(jì)算結(jié)果Tab.2 Numerical simulation and test results under three velocity conditions

4 結(jié) 論

1)本文建立的模型能夠準(zhǔn)確模擬冰在拉伸時(shí)的斷裂現(xiàn)象和裂紋擴(kuò)展模式。

2)冰的四點(diǎn)彎曲過程中，裂紋出現(xiàn)在固定剛性支撐的位置。

3)冰梁彎曲后分裂成3部分，斷裂位置對應(yīng)于底部壓頭加載位置。該現(xiàn)象和試驗(yàn)現(xiàn)象具有極高的一致性。

4)對于冰的四點(diǎn)彎曲過程，加載速度為Vsupport=0.003 580 m/s時(shí)的彎曲最大載荷最大，斷裂時(shí)間最短。