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

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

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

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

近場(chǎng)動(dòng)力學(xué)方法(PD)，是近二十年新興的一種非局部、無(wú)網(wǎng)格粒子、積分形式的計(jì)算方法，在斷裂問(wèn)題域和非斷裂問(wèn)題域上都可以連續(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é)等多種物理問(wèn)題計(jì)算分析領(lǐng)域[25-27]。然而PD應(yīng)用于冰計(jì)算領(lǐng)域剛剛起步，多數(shù)文獻(xiàn)采用基于鍵型的彈性模型來(lái)模擬冰的破壞特性，例如，冰與剛體柱的作用過(guò)程、冰槳作用過(guò)程等，忽略了泊松比的真實(shí)狀態(tài)和海冰具有塑性屈服過(guò)程的力學(xué)特性。由此本文主要的目的是建立基于狀態(tài)型PD方法的冰彈塑性本構(gòu)模型，消除了鍵型對(duì)泊松比的限制，加入了冰破壞過(guò)程塑性變形過(guò)程，為后續(xù)的冰與結(jié)構(gòu)物的作用過(guò)程提供支撐。

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

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

1.1.1 近場(chǎng)動(dòng)力學(xué)理論

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

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)表示，對(duì)k點(diǎn)有作用的域(近場(chǎng)域)表示為Hx(k)。

1.1.2 塑性理論

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

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

(3)

(4)

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

(5)

1.1.3 屈服函數(shù)

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

(6)

(7)

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

(8)

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

(9)

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

(10)

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

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

(11)

其中

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

1.2 邊界條件

由于PD運(yùn)動(dòng)方程不包含任何空間導(dǎo)數(shù)，因此一般不需要約束條件來(lái)求解積分微分方程。與局部理論不同，邊界條件是通過(guò)非零體積的虛擬邊界層施加的?；跀?shù)值實(shí)驗(yàn)，Macek等[30]建議虛擬邊界層的范圍等于領(lǐng)域尺寸δ，以確保施加的規(guī)定約束在真實(shí)區(qū)域中得到準(zhǔn)確反映。位移邊界條件可以通過(guò)對(duì)虛擬層中的材料點(diǎn)進(jìn)行約束來(lái)施加，使邊界曲面上的條件得到明確滿足。因此，虛擬層中的位移值可以基于實(shí)域值和邊界條件指定值的線性外推來(lái)近似。類似地，也可以通過(guò)近似虛擬區(qū)域中的牽引力值來(lái)施加牽引力邊界條件，從而使真實(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)為極限伸長(zhǎng)量，當(dāng)變形后粒子間的伸長(zhǎng)量達(dá)到極限伸長(zhǎng)量s0時(shí)，則判斷為失效，此過(guò)程不可逆。極限伸長(zhǎng)量可以通過(guò)材料的極限能量釋放率得到，粒子間的伸長(zhǎng)量和力密度的關(guān)系呈線性。然而在塑性變形中，粒子間力密度和伸長(zhǎng)量的關(guān)系呈典型的非線性，如圖2所示，因此彈塑性本構(gòu)模型中的失效準(zhǔn)則不能再簡(jiǎn)單地使用最終狀態(tài)的伸長(zhǎng)量來(lái)判斷，必須對(duì)整個(gè)變形過(guò)程伸長(zhǎng)量對(duì)應(yīng)的力密度積分來(lái)判斷材料的失效。圖2曲線下的面積代表彈塑性變形微勢(shì)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的對(duì)應(yīng)參數(shù)，在一維時(shí)取值為12，二維時(shí)取值為36，三維時(shí)取值為560。

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

(17)

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

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

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

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

(18)

時(shí)間積分內(nèi)的數(shù)值執(zhí)行過(guò)程將全部替換為塑性變形的數(shù)值過(guò)程，可參考文獻(xiàn)[28-29]中的數(shù)值實(shí)現(xiàn)方法。對(duì)于失效準(zhǔn)則的積分過(guò)程，式(13)將采用如下積分過(guò)程實(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ì)算了帶孔二維板的變形問(wèn)題。二維板的尺寸選取長(zhǎng)和寬分別為L(zhǎng)=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ū)域開(kāi)始。圖5中給出了帶孔板的等效應(yīng)力和位移的變化，與有限元結(jié)果對(duì)比完全一致，該案例驗(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)可以用來(lái)測(cè)量冰的彎曲強(qiáng)度。為了驗(yàn)證本文冰模型的準(zhǔn)確性，選用Ehlers和Kujala進(jìn)行的試驗(yàn)[32]作為對(duì)比驗(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(zhǎng)=4.25 m,H=0.4 m和B=0.365 m。上方的固定支撐位于中點(diǎn)向兩邊2 m處，底部位移壓頭位于0.5 m處。上方的支撐在空間中固定，下方的支撐以Vsupport=0.003 580 m/s勻速向上移動(dòng)。海冰參數(shù)的選取參考試驗(yàn)和文獻(xiàn)中對(duì)冰力學(xué)性能相關(guān)分析設(shè)定[32-34]，冰材料的極限能量釋放率采用了文獻(xiàn)[35]中的參考值，見(jiàn)表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è)移動(dòng)支撐之間的冰梁位移最大，且這段部分各個(gè)位置的位移量基本相同。冰梁上越遠(yuǎn)離移動(dòng)支撐部位的縱向位移變形越小，這樣的變形一直持續(xù)到0.31 s之前。此時(shí)海冰未發(fā)生破壞，由于固定支撐和加載支撐的約束原因，冰x方向和y方向的位移符合物理規(guī)律。還可以通過(guò)等效應(yīng)力分布圖觀測(cè)到：高應(yīng)力分布在拉伸和壓縮集中的部位，例如冰梁的上側(cè)邊緣和下側(cè)邊緣，如圖5(c)所示，這兩側(cè)的位置是冰梁擠壓和拉伸最為嚴(yán)重的位置。該應(yīng)力現(xiàn)象與文獻(xiàn)[35]中通過(guò)有限元方法計(jì)算的結(jié)果具有很高的一致性。從圖8可以觀測(cè)到：冰梁彎曲后分裂成3部分，斷裂位置對(duì)應(yīng)于底部壓頭加載位置。該現(xiàn)象和試驗(yàn)現(xiàn)象具有極高的一致性。另外為了進(jìn)一步驗(yàn)證計(jì)算結(jié)果的準(zhǔn)確性，對(duì)比數(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)行對(duì)比驗(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)對(duì)應(yīng)的極值載荷數(shù)值極其發(fā)生時(shí)間，從表中可以看出試驗(yàn)的結(jié)果和計(jì)算值的誤差在10%以內(nèi)，說(shuō)明本文建立的無(wú)網(wǎng)格法冰本構(gòu)模型能夠模擬冰的彎曲破壞過(guò)程。

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

表2 數(shù)值計(jì)算和試驗(yàn)中3種速度工況對(duì)應(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)彎曲過(guò)程中，裂紋出現(xiàn)在固定剛性支撐的位置。

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

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