黃焱 王建平 孫劍橋

(天津大學(xué)水利工程仿真與安全國(guó)家重點(diǎn)實(shí)驗(yàn)室,天津 300350)

引言

近年來(lái),環(huán)北極與近北極國(guó)家相繼出臺(tái)了新的極地政策或戰(zhàn)略部署[1-3],極地海域的商業(yè)航運(yùn)、資源開(kāi)發(fā)、科學(xué)考察、環(huán)境保護(hù)及軍事保障等活動(dòng)也隨之日趨白熱化.海冰是極地海域標(biāo)志性的環(huán)境條件,也是極地航行船舶和極地作業(yè)裝備在安全性保障上面對(duì)的首要挑戰(zhàn).為應(yīng)對(duì)這一挑戰(zhàn),學(xué)術(shù)界與工程界開(kāi)展了大量關(guān)于冰-結(jié)構(gòu)相互作用過(guò)程的數(shù)值模擬方法研究[4-6],以期能夠?yàn)楣こ萄b備在現(xiàn)實(shí)冰環(huán)境條件下的承載過(guò)程提供可靠的預(yù)報(bào).然而,天然海冰因其復(fù)雜的力學(xué)性質(zhì),致使當(dāng)前的數(shù)值模擬結(jié)果直接面向工程應(yīng)用時(shí)尚存在極大的阻礙.因此,如何建立可靠的冰力學(xué)性質(zhì)的數(shù)值模擬方法,已成為當(dāng)前學(xué)術(shù)界關(guān)注的焦點(diǎn)問(wèn)題之一.

已有大量的研究表明,天然冰材料在變形和破壞行為上都具有顯著的各向異性特征[7-10],這一特征成為冰-結(jié)構(gòu)相互作用過(guò)程中產(chǎn)生復(fù)雜載荷過(guò)程的關(guān)鍵誘因.因此,在數(shù)值模擬中準(zhǔn)確再現(xiàn)冰材料這一關(guān)鍵的力學(xué)特征,就成為保障數(shù)值預(yù)報(bào)結(jié)果合理性的必要條件.基于此,一些學(xué)者探索了將天然冰材料表現(xiàn)出的整體各向異性特征引入數(shù)值模擬方法的途徑,如Castelnau 等[11]探索了一種模擬多晶冰黏塑性行為各向異性的自洽方法;Hibler 和Schulson 等[12]將缺陷導(dǎo)致的各向異性引入到海冰數(shù)值模型中,并據(jù)此進(jìn)行了雙軸壓縮模擬,探討了模型中海冰的裂紋擴(kuò)展模式;季順迎和岳前進(jìn)[13]針對(duì)海冰的動(dòng)力過(guò)程模擬問(wèn)題,討論了彈塑性各向異性模型的適用性;文獻(xiàn)[14]基于原位冰芯樣本的實(shí)測(cè)數(shù)據(jù),建立了不同晶體織構(gòu)下的彈性剛度張量,并將其引入冰層力學(xué)模型之中,由此計(jì)算和分析了多晶冰織構(gòu)特征所導(dǎo)致的各向異性對(duì)冰內(nèi)地震波傳播方式的影響.但正如冰力學(xué)研究領(lǐng)域知名學(xué)者E.M.Schulson 所言,多晶冰材料(見(jiàn)圖1)的各向異性特征根源于單晶冰的各向異性[15].所以,對(duì)于天然冰材料各向異性特征的模擬,首先應(yīng)建立針對(duì)單晶冰各向異性的再現(xiàn)方法,由此構(gòu)建出模擬天然冰材料各向異性特征的基本框架,再進(jìn)一步由此框架擴(kuò)展出多晶冰材料的模擬方法,才可能全面、準(zhǔn)確地再現(xiàn)天然冰材料的特殊力學(xué)性質(zhì).本文將針對(duì)當(dāng)前學(xué)術(shù)界缺乏單晶冰材料各向異性模擬方法的局面,開(kāi)展探索性研究工作.

圖1 S2 冰的水平向與垂向晶體切片[15]Fig.1 Horizontal and vertical sections of S2 ice[15]

另一方面,在數(shù)值模擬方法上,基于非局部作用思想的近場(chǎng)動(dòng)力學(xué)方法,因其在統(tǒng)一連續(xù)介質(zhì)力學(xué)理論與粒子模擬方法上具有突出的潛質(zhì),已在學(xué)術(shù)界引起廣泛關(guān)注[16-19].在近場(chǎng)動(dòng)力學(xué)模型中,物質(zhì)點(diǎn)間的相互作用力函數(shù)(即“本構(gòu)力函數(shù)”)包含了材料的物性信息,它不再以傳統(tǒng)的應(yīng)力-應(yīng)變關(guān)系等形式出現(xiàn),也不假設(shè)位移場(chǎng)的連續(xù)性,求解時(shí)不存在空間求導(dǎo).此外,位移場(chǎng)的連續(xù)與否并不影響基本方程的求解,不連續(xù)性(如裂紋、界面)在求解運(yùn)動(dòng)方程時(shí)自然產(chǎn)生,因而不會(huì)出現(xiàn)傳統(tǒng)連續(xù)性模型求解時(shí)出現(xiàn)的病態(tài)特征[18].當(dāng)前,近場(chǎng)動(dòng)力學(xué)方法已被大量應(yīng)用于巖石、復(fù)合材料裂紋萌生擴(kuò)展過(guò)程的模擬[20-22],并在材料復(fù)雜力學(xué)行為特征的模擬上顯示出突出的優(yōu)勢(shì).鑒于此,這一方法也開(kāi)始被引入冰-結(jié)構(gòu)相互作用問(wèn)題的模擬中,并因其在冰材料斷裂破壞現(xiàn)象模擬上的良好表現(xiàn)而備受青睞[23-26].需要指出的是,目前近場(chǎng)動(dòng)力學(xué)方法在應(yīng)用于工程預(yù)報(bào)時(shí),都還是基于簡(jiǎn)單的各向同性材料假設(shè),進(jìn)而導(dǎo)致已有的模擬結(jié)果尚難以合理再現(xiàn)真實(shí)作用場(chǎng)景.

據(jù)此,本文在回顧單晶冰彈性各向異性力學(xué)表征的基礎(chǔ)上,采用近場(chǎng)動(dòng)力學(xué)方法建立單晶冰態(tài)基模型,通過(guò)對(duì)單軸壓縮過(guò)程的數(shù)值模擬,探討彈性各向異性特征的引入方式及相關(guān)的參數(shù)標(biāo)定方法,以期建立合理可靠的單晶冰彈性各向異性數(shù)值模擬方法.

1 單晶冰彈性各向異性的力學(xué)特征

在小位移條件下,應(yīng)力與應(yīng)變之間的關(guān)系為

其中,σi為應(yīng)力分量,i,j=1,2,···,6;εi為應(yīng)變分量;Sij為柔度矩陣分量.

對(duì)于單晶冰所具備的六方晶系結(jié)構(gòu),如果假設(shè)其基面內(nèi)符合力學(xué)各向同性,則可分別將晶體的C軸方向及其基面內(nèi)的任意一對(duì)正交方向分別定義為X3,X1和X2應(yīng)力主軸,相應(yīng)的柔度矩陣中只有5 個(gè)非零參數(shù)[9]

Fletcher[9]收集整理了大量的單晶冰內(nèi)彈性波傳播的實(shí)驗(yàn)測(cè)試數(shù)據(jù),并將這些數(shù)據(jù)回歸為描述單晶冰彈性各向異性的柔度矩陣表達(dá),矩陣中各項(xiàng)系數(shù)的值分別為:S11=1040 GPa?1;S12=?430 GPa?1;S13=?240 GPa?1;S33=850 GPa?1;S44=3140 GPa?1.

進(jìn)一步利用張量變換可獲得任一方向有效柔度的表達(dá)式為[9,15]

其中,Φ表示任一加載方向與C軸的夾角.而各個(gè)變形方向上的楊氏模量為該方向有效柔度的倒數(shù)

其中,Eref為參考楊氏模量.該值基于Fletcher 實(shí)驗(yàn)測(cè)試獲得的單晶冰柔度矩陣系數(shù)計(jì)算得到,并將作為本文數(shù)值模型后續(xù)標(biāo)定的參考值.

圖2 為由式(3)和式(4)繪制得到的單晶冰彈性各向異性特征表征圖:圖中的實(shí)線描繪的是Eref在與C軸平行的任意平面內(nèi)的空間分布情況,該曲線上任意一點(diǎn)的橫縱坐標(biāo)比值為該平面內(nèi)變形方向與C軸夾角Φ的正切值,該點(diǎn)距離坐標(biāo)原點(diǎn)的長(zhǎng)度為相應(yīng)方向上的楊氏模量值;圖中的虛線表征以C軸方向楊氏模量為特征值的各向同性邊界.

圖2 單晶冰彈性各向異性表征圖[9]Fig.2 Reference Young's modulus as a function of orientation[9]

由圖2 可知,楊氏模量的分布是一條關(guān)于坐標(biāo)原點(diǎn)中心對(duì)稱和坐標(biāo)軸雙對(duì)稱的閉合曲線,故對(duì)各向異性的描述可限定于第一象限內(nèi).

2 態(tài)基近場(chǎng)動(dòng)力學(xué)模型

2.1 本構(gòu)模型

在近場(chǎng)動(dòng)力學(xué)模型的發(fā)展歷經(jīng)了兩個(gè)階段,即鍵基模型和態(tài)基模型.鍵基模型的理論基礎(chǔ)會(huì)導(dǎo)致恒定泊松比的天然缺陷[18],為彌補(bǔ)這一理論缺陷,又進(jìn)一步發(fā)展出了態(tài)基模型[27].因此,態(tài)基模型已成為當(dāng)前近場(chǎng)動(dòng)力學(xué)模型應(yīng)用中的主角[28-29].下面對(duì)態(tài)基模型的基礎(chǔ)理論進(jìn)行簡(jiǎn)要介紹.

在態(tài)基近場(chǎng)動(dòng)力學(xué)模型中,參考構(gòu)型中x位置處材料點(diǎn)在t時(shí)刻的運(yùn)動(dòng)方程為

其中ρ為參考構(gòu)型中材料的密度;u為位移向量場(chǎng);Hx為x材料點(diǎn)的“近場(chǎng)范圍”;T(x′,x,t) 為材料點(diǎn)x′施加在x上的力密度向量;T(x,x′,t) 為材料點(diǎn)x施加在x′上的力密度向量;T(x′,x,t)?T(x,x′,t) 對(duì)應(yīng)于鍵基模型中的對(duì)點(diǎn)力函數(shù);dVx′ 為材料點(diǎn)x′的體積;b為體積力向量場(chǎng).

材料點(diǎn)x′和x之間的相互作用通過(guò)“鍵”來(lái)實(shí)現(xiàn).以ξ表示參考構(gòu)型中材料點(diǎn)x′和x的相對(duì)位置

在t時(shí)刻的當(dāng)前構(gòu)型中,材料點(diǎn)x′和x的相對(duì)位移為

故向量 ξ +η 表示兩材料點(diǎn)在當(dāng)前構(gòu)型中的相對(duì)位置.

定義材料點(diǎn)x′和x之間“鍵”的伸長(zhǎng)率為

其中,當(dāng)“鍵”處于拉伸狀態(tài)時(shí),伸長(zhǎng)率為正;當(dāng)“鍵”處于壓縮狀態(tài)時(shí),伸長(zhǎng)率為負(fù).

定義材料點(diǎn)x處的膨脹度為

其中d為與材料無(wú)關(guān)的參數(shù).

在線性態(tài)基近場(chǎng)動(dòng)力學(xué)固體模型中,力密度向量與當(dāng)前構(gòu)型中的相對(duì)位置向量 ξ +η 平行,表示為

其中,ω為影響函數(shù),可根據(jù)材料點(diǎn)x′,x之間的距離控制它們的相互作用大小;a和b為與材料有關(guān)的參數(shù);Λ 為

在近場(chǎng)動(dòng)力學(xué)中,認(rèn)為處于近場(chǎng)范圍內(nèi)的材料點(diǎn)均會(huì)與中心材料點(diǎn)發(fā)生相互作用,即非局部相互作用,這是近場(chǎng)動(dòng)力學(xué)與傳統(tǒng)連續(xù)介質(zhì)力學(xué)最重要的區(qū)別之一.

文獻(xiàn)[30]認(rèn)為,影響函數(shù)應(yīng)當(dāng)包含 δ 作為長(zhǎng)度尺度的對(duì)近場(chǎng)范圍內(nèi)其他材料點(diǎn)的影響,因此其提出以下影響函數(shù)

將影響函數(shù)代入式(10)可得態(tài)基近場(chǎng)動(dòng)力學(xué)的力密度向量表達(dá)

其中,k和G分別為體積模量和剪切模量,在三維情況下分別為

其中,E為楊氏模量;ν 為泊松比.

2.2 時(shí)間積分格式

近場(chǎng)動(dòng)力學(xué)運(yùn)動(dòng)方程是一個(gè)積分-微分方程,它的解通過(guò)空間和時(shí)間的數(shù)值積分來(lái)獲得.本文采用無(wú)網(wǎng)格離散化方法,將材料離散為一系列具有確定體積的粒子,離散后粒子xi的運(yùn)動(dòng)方程變?yōu)?/p>

其中N為粒子xi近場(chǎng)范圍內(nèi)的粒子數(shù),下標(biāo)m表示第m個(gè)時(shí)間步.

粒子xi膨脹度θi為

其中,下標(biāo)ij代表該物理量與粒子xi和xj相關(guān).

在獲得tm時(shí)刻的加速度后,采用以下顯式積分算法獲得下一時(shí)間步的速度和位移[30]

其中Δt為計(jì)算時(shí)間步長(zhǎng).

顯式積分算法雖然簡(jiǎn)單明了,但是其只是條件穩(wěn)定的,當(dāng)時(shí)間步長(zhǎng)大于最大容許時(shí)間步長(zhǎng)之后便不能獲得穩(wěn)定收斂解.Madenci 等利用馮·諾依曼穩(wěn)定性分析得到了臨界時(shí)間步長(zhǎng)為[30]

3 基于各向同性假設(shè)的單晶冰基礎(chǔ)態(tài)基模型

本文的近場(chǎng)動(dòng)力學(xué)模擬借助C++程序和Visual Studio 平臺(tái)實(shí)現(xiàn),有限元模擬采用ANSYS 軟件完成.

在探究如何將各向異性引入單晶冰近場(chǎng)動(dòng)力學(xué)模型的方法之前,首先以圖2 中展示的各向同性假設(shè)邊界為基準(zhǔn),對(duì)模型粒子數(shù)量(或粒子間距)取值進(jìn)行標(biāo)定.

如圖2 所示,各向同性假設(shè)下的楊氏模量為11.76 GPa,相應(yīng)地,體積模量k=8.7 GPa;剪切模量G=4.62 GPa.

不論是天然淡水湖冰還是海冰材料,其中包含的晶粒直徑大多集中在2～ 10 mm 范圍內(nèi)[31-33].據(jù)此,本文采用的單晶冰試樣模型尺寸為2.5 mm × 2.5 mm ×5.0 mm,密度為0.935 mg/mm3,如圖3 所示.

圖3 各向同性模型及載荷條件(單位:mm)Fig.3 Isotropic model and loading conditions (unit:mm)

令n為試件橫截面長(zhǎng)度與粒子間距的比值.本文中,n分別取為8,10,16,20,25,40,此時(shí)粒子間距 ?x分別為0.3125 mm,0.25 mm,0.15625 mm,0.125 mm,0.1 mm,0.0625 mm,據(jù)此建立了6 種模型,相應(yīng)的粒子數(shù)分別為1024,2000,8192,16000,31250,128000.近場(chǎng)范圍半徑δ取為3?x[34-36].Silling 和Askari[34]建議將邊界條件通過(guò)一個(gè)有限厚度帶狀區(qū)域內(nèi)的粒子定義.本節(jié)向模型頂部和底部同時(shí)施加大小相同、方向相反的恒力壓縮載荷,大小為F=100 N,通過(guò)厚度為δ的帶狀區(qū)域內(nèi)的粒子體積力定義b

其中,M為頂部或者底部δ厚度區(qū)域內(nèi)的粒子數(shù).

由于本節(jié)施加的載荷為靜力載荷,為獲得穩(wěn)定解,借助自適應(yīng)動(dòng)態(tài)松弛算法使模型快速達(dá)到穩(wěn)定狀態(tài)[37].同時(shí)采用商用有限元軟件ANSYS 建立相同尺度和力學(xué)性質(zhì)的有限元模型,施加相同載荷,并進(jìn)行靜力分析.提取近場(chǎng)動(dòng)力學(xué)模型中位于Z軸上的所有粒子的Z向位移,與有限元計(jì)算結(jié)果進(jìn)行對(duì)比,得到結(jié)果如圖4 所示.

從圖4 可以看出,當(dāng)本模型中n從25 增加到40后,計(jì)算結(jié)果與有限元結(jié)果的差異均在2%以內(nèi),然而計(jì)算量卻增加了近6 倍.因此,本文將 ?x=L/25作為離散粒子間距,可以滿足計(jì)算精度和效率需求.

圖4 各粒子間距模型中心軸Z 向位移與有限元結(jié)果的對(duì)比Fig.4 Displacement variation along the loading direction

4 各向異性特征的引入方式與參數(shù)標(biāo)定

在近場(chǎng)動(dòng)力學(xué)理論中,要體現(xiàn)材料的彈性各向異性,需要從力密度的本構(gòu)關(guān)系入手.在力密度向量的表達(dá)式(10)中,括號(hào)中含參數(shù)a的第一項(xiàng)表示體積膨脹度對(duì)力密度的貢獻(xiàn),含參數(shù)b的第二項(xiàng)可表示體積膨脹度以外的變形對(duì)力密度的貢獻(xiàn).本文認(rèn)為,材料的彈性各向異性行為受體積變形和其他變形的綜合影響,因此本文提出了基于影響函數(shù)ω的各向異性模型.

4.1 基于影響函數(shù)的彈性各向異性模型

Silling 等[27]在提出態(tài)基近場(chǎng)動(dòng)力學(xué)理論時(shí)定義影響函數(shù) ω 為:一個(gè)位于中心材料點(diǎn)近場(chǎng)范圍內(nèi)的嚴(yán)格為正的標(biāo)量值函數(shù).若影響函數(shù)僅與“鍵”長(zhǎng)有關(guān),即表示為,則影響函數(shù)為球形.目前,大量研究工作都是基于球形的影響函數(shù)[30,38-39],所表征的僅為材料的各向同性.Silling 等[27]指出,對(duì)于不規(guī)則形狀邊界的結(jié)構(gòu)體或者混合物表面的模擬,可以通過(guò)定義形狀不同的影響函數(shù)來(lái)實(shí)現(xiàn).本文即遵循這一思路,擬通過(guò)定義合適的影響函數(shù),以期使模型呈現(xiàn)出各向異性特征.

在一般的態(tài)基近場(chǎng)動(dòng)力學(xué)理論中,力密度如式(10)所示.本文初步假定材料的各向異性完全由影響函數(shù)ωij控制,模型中的其他部分與各向同性材料相同,并且各向同性部分的材料參數(shù)以單晶冰C軸方向的力學(xué)屬性為基準(zhǔn).

由于在各向同性理論中,力密度中的模量參數(shù)均與楊氏模量成正相關(guān),因此,假設(shè)在各向異性模型中,各模量參數(shù)和楊氏模量具有相同的關(guān)于角度的表達(dá),即力密度向量為

式(23)中的分母項(xiàng)可視為單晶冰楊氏模量隨加載方向變化的歸一化表征,由式(3)和式(4)得到;φ為模型中“鍵”與C軸的夾角.

模型約束端的總反力由所有與被約束粒子相關(guān)的力密度集體貢獻(xiàn),表示為

其中,R為總反力;M為被約束的粒子xi的數(shù)量;N為粒子xi近場(chǎng)范圍內(nèi)的粒子數(shù).

在一般單軸壓縮試驗(yàn)中,名義應(yīng)力定義為兩端反力與構(gòu)件橫截面積的比值,應(yīng)變定義為構(gòu)件標(biāo)距段的伸長(zhǎng)量與標(biāo)距的比值,即

其中,F為反力;A0為構(gòu)件變形前的橫截面面積;?l為標(biāo)距段的伸長(zhǎng)量;l0為標(biāo)距.

類似地,定義近場(chǎng)動(dòng)力學(xué)模型在單軸壓縮模擬中的名義應(yīng)力為底部反力R與變形前橫截面面積的比值;標(biāo)距為頂部虛粒子層底部與底部虛粒子層頂部之間的距離;名義應(yīng)變?yōu)闃?biāo)距段的伸長(zhǎng)量與標(biāo)距的比值.結(jié)構(gòu)如圖5 所示.模型尺寸與第3 節(jié)中的基礎(chǔ)模型相同,C軸始終位于XZ平面內(nèi).頂部的應(yīng)變率載荷通過(guò)設(shè)置三層虛粒子向下的速度為應(yīng)變率與試樣高度的乘積來(lái)施加,底部固定約束通過(guò)剛性固定三層虛粒子的位置來(lái)施加.本文中應(yīng)變率大小均為10 s?1.

圖5 恒定應(yīng)變率加載模型(單位:mm)Fig.5 Constant strain rate loading model (unit:mm)

在單晶冰中,C軸和垂直于C軸方向的力學(xué)屬性是非常重要的兩個(gè)屬性;同時(shí),單晶冰的楊氏模量在加載方向與C軸呈45°附近時(shí)最小,因此這三個(gè)方向的楊氏模量將是本文考慮的重點(diǎn).在載荷沿Z方向加載的前提下,本文通過(guò)在XZ平面內(nèi)順時(shí)針旋轉(zhuǎn)C軸的角度,實(shí)現(xiàn)C軸與加載方向夾角Φ的調(diào)節(jié).采用式(23)的影響函數(shù)進(jìn)行Φ=0°,45°,90°三個(gè)不同方向上的恒定應(yīng)變率壓縮變形模擬,得到名義應(yīng)力隨名義應(yīng)變的變化關(guān)系.這三個(gè)加載方向的名義應(yīng)力-名義應(yīng)變關(guān)系如圖6 所示.

圖6 名義應(yīng)力-名義應(yīng)變曲線Fig.6 Nominal stress-nominal strain curve

對(duì)于線彈性模型,名義應(yīng)力-名義應(yīng)變曲線近似為直線,其斜率為楊氏模量.采用最小二乘法對(duì)該斜率進(jìn)行擬合,即可得到數(shù)值模型的楊氏模量.在式(23)的影響函數(shù)下,三個(gè)加載方向的楊氏模量分別為:E0=8.461 GPa,E45=7.805 GPa,E90=7.877 GPa;下標(biāo)數(shù)字表示加載方向與C軸方向的夾角.

觀察模擬結(jié)果并將其與圖2 的參考楊氏模量進(jìn)行對(duì)比可以發(fā)現(xiàn),盡管將楊氏模量與角度的關(guān)系直接代入影響函數(shù)中,能夠使數(shù)值模型呈現(xiàn)一定程度的各向異性,但其結(jié)果與實(shí)驗(yàn)測(cè)試得到的參考楊氏模量存在較大差異,仍需對(duì)該影響函數(shù)進(jìn)行修正.

4.2 影響函數(shù)修正方案一:引入輔助參數(shù)α

如前所述,影響函數(shù)雖然初步繼承了楊氏模量隨加載方向變化的力學(xué)表征關(guān)系,但計(jì)算結(jié)果與參考楊氏模量仍存在數(shù)量上的差異,尤其是沿C軸加載方向.據(jù)此,本文初步推斷,若使模型能夠模擬出單晶冰沿各個(gè)加載方向的楊氏模量,需要在影響函數(shù)中加入輔助參數(shù)α的控制

式中α可從整體水平上調(diào)節(jié)影響函數(shù)的值.式(23)可視為式(27)在α=1 時(shí)的特殊情況.此時(shí)的影響函數(shù)存在以下關(guān)系

式(28)表示當(dāng)任意一根“鍵”ξij的影響函數(shù)中的α變?yōu)樵瓉?lái)的ζ倍后,影響函數(shù)的值變?yōu)樵瓉?lái)的1/ζ倍,在相同變形下,力密度也變?yōu)樵瓉?lái)的1/ζ倍.根據(jù)式(25),反力也變?yōu)樵瓉?lái)的1/ζ倍,則楊氏模量也應(yīng)變?yōu)樵瓉?lái)的1/ζ.

考慮到C軸方向的特殊性,首先調(diào)整α使沿C軸方向的楊氏模量與參考楊氏模量匹配.根據(jù)上一節(jié)計(jì)算得到的C軸方向的楊氏模量E0=8.461 GPa,令=0.7195.將該α控制下的影響函數(shù)代入數(shù)值模型,進(jìn)行0°,45°和90°三個(gè)方向的加載模擬,最終得到該三個(gè)方向的楊氏模量分別為:E0=11.84 GPa,E45=10.81 GPa,E90=10.94 GPa.

將上述楊氏模量計(jì)算結(jié)果分別與α=1 時(shí)的計(jì)算結(jié)果進(jìn)行對(duì)比,可以得到三個(gè)方向楊氏模量的比值分別為:0.7146,0.7220 和0.7200.這些比值與α=0.7195 非常接近,表明前述楊氏模量隨α增大而減小的推斷成立.

然而,盡管在α=0.7195 時(shí)沿C軸方向的數(shù)值模型楊氏模量與參考楊氏模量非常接近,但在另外兩個(gè)方向的差異較大.進(jìn)一步,繪制α=1 與α=0.7195下各加載方向的楊氏模量計(jì)算結(jié)果,并與參考楊氏模量進(jìn)行對(duì)比,如圖7 所示.可以發(fā)現(xiàn),數(shù)值模型楊氏模量隨加載方向變化的曲線形狀在兩種α值下幾乎沒(méi)有變化,即模型的各向異性程度沒(méi)有得到改善.因此,仍需進(jìn)一步修正影響函數(shù)中的角度控制參量.

圖7 兩種不同α 時(shí)的楊氏模量對(duì)比Fig.7 Comparison of Young's modulus for two different α

4.3 影響函數(shù)修正方案二:引入輔助參數(shù)α,β,γ

修正方案一的計(jì)算結(jié)果表明,僅通過(guò)單一的輔助控制參數(shù)并不能調(diào)整模型的各向異性程度.因此本文進(jìn)一步推斷,雖然影響函數(shù)與楊氏模量具有相似的關(guān)于角度的函數(shù)表達(dá),但表達(dá)式中的各個(gè)分項(xiàng)的系數(shù)并不一定相同.據(jù)此,進(jìn)一步假定影響函數(shù)形式為

其中,α,β,γ為獨(dú)立變化的三個(gè)輔助控制參數(shù).

將該影響函數(shù)代入力密度的公式,分別令φ=0°,45°和90°,得到沿這三個(gè)方向的力密度為

由式(30)可知,沿C軸方向“鍵”的力密度主要受α的控制,垂直于C軸方向主要受β的控制,其他方向?qū)⑹艿饺齻€(gè)參數(shù)的綜合影響.由于S44≥S11>S33,故與C軸呈45°夾角的“鍵”受γ影響最大.基于加載方向附近的“鍵”對(duì)該方向楊氏模量影響更大的思想,本文提出以下關(guān)于α,β,γ標(biāo)定方法,上標(biāo)(l)表示第l次標(biāo)定

(1)l=0,令α(0),β(0),γ(0)的初始值均為0.7195,針對(duì)0°,45°,90°三個(gè)加載方向上的楊氏模量容許誤差為Δ(本文設(shè)定為0.05%);

(2) 以α(l),β(l),γ(l)建立相應(yīng)的近場(chǎng)動(dòng)力學(xué)模型,進(jìn)行三個(gè)方向的恒應(yīng)變率單軸壓縮數(shù)值模擬,分別得到對(duì)應(yīng)的楊氏模量;

(3) 定義第l步時(shí)的總相對(duì)誤差為

計(jì)算此時(shí)的總相對(duì)誤差 ?E(l);

(4) 若 ?E(l)≤?,輸出α,β,γ的最終值為α(l),β(l),γ(l),反之,令

(5) 令l=l+1,重復(fù)第 (2)～ (4)步.

將α(0)=β(0)=γ(0)=0.7195 作為第一次循環(huán)的起始條件輸入,執(zhí)行上述標(biāo)定過(guò)程.圖8 展示了α,β,γ值以及三個(gè)方向上楊氏模量的相對(duì)誤差隨迭代次數(shù)的變化.可以看出,隨著迭代次數(shù)的增加,α,β,γ分別穩(wěn)定收斂于不同的值,同時(shí)三個(gè)方向的相對(duì)誤差隨著迭代次數(shù)的增加而下降,這表明本文提出的迭代標(biāo)定方法基本有效.此外,總誤差也隨著迭代次數(shù)的增加迅速下降,如圖9 所示.經(jīng)過(guò)22 次循環(huán)后,總誤差小于設(shè)定的容許誤差Δ,此時(shí)α(22)=0.5364,β(22)=0.7506,γ(22)=1.2425,與此相對(duì)應(yīng)的三個(gè)方向楊氏模量的計(jì)算結(jié)果分別為:E0=11.76 GPa,E45=8.793 GPa,E90=9.614 GPa.

圖8 α,β,γ 及三個(gè)方向楊氏模量的相對(duì)誤差隨著迭代次數(shù)的變化Fig.8 α,β,γ and the relative error of Young's modulus as a function of iteration number

圖9 總誤差隨著迭代次數(shù)的變化Fig.9 Total error as a function of iteration number

4.4 其他加載方向的適用性驗(yàn)證

影響函數(shù)修正方案二能夠較為準(zhǔn)確地模擬沿0°,45°和90°三個(gè)加載方向的彈性變形.進(jìn)一步,對(duì)其他加載方向模擬結(jié)果的準(zhǔn)確性開(kāi)展驗(yàn)證.

將第22 次循環(huán)時(shí)的α,β,γ值作為終值,計(jì)算在此終值下與C軸方向分別呈15°,30°,60°和75°的四個(gè)加載方向的楊氏模量,及其與參考楊氏模量的相對(duì)誤差,如表1 所示.圖10 從曲線形狀上展示了數(shù)值模型得到的楊氏模量與參考楊氏模量的對(duì)比.

圖10 楊氏模量隨加載方向的變化趨勢(shì)對(duì)比Fig.10 Graphic comparison of Young's modulus under different loading directions

從表1 可以看出,其他四個(gè)加載方向的數(shù)值模型楊氏模量與參考楊氏模量的相對(duì)誤差均位于5%以下,處于可接受的誤差范圍內(nèi).同時(shí),數(shù)值模型楊氏模量隨加載方向的變化趨勢(shì)與參考楊氏模量較為契合.這進(jìn)一步驗(yàn)證了,本文提出的基于影響函數(shù)的彈性各向異性近場(chǎng)動(dòng)力學(xué)模型,能夠?qū)崿F(xiàn)單晶冰彈性各向異性力學(xué)表征的有效模擬.

表1 四個(gè)加載方向的楊氏模量計(jì)算結(jié)果及其相對(duì)誤差Table 1 Young's modulus of the numerical model under different loading directions and the relative error

5 結(jié)論

本文基于態(tài)基近場(chǎng)動(dòng)力學(xué)理論,根據(jù)單晶冰楊氏模量沿不同加載方向的變化規(guī)律,推斷粒子間“鍵”的變形與力密度也存在類似的規(guī)律,提出了一種從影響函數(shù)出發(fā)的彈性各向異性模型.以前人實(shí)驗(yàn)測(cè)試得到的楊氏模量值為參考,通過(guò)開(kāi)展與C軸呈0°,45°和90°三個(gè)加載方向的單晶冰單軸壓縮數(shù)值模擬,提出了針對(duì)該影響函數(shù)的修正和輔助參數(shù)標(biāo)定方法,并且最終在15°,30°,60°和75°等其他四個(gè)加載方向進(jìn)行了驗(yàn)證.

計(jì)算結(jié)果表明,本文所提出的針對(duì)影響函數(shù)的修正與參數(shù)標(biāo)定方法,能夠較為便捷地找到數(shù)值模型楊氏模量與參考楊氏模量相一致的影響函數(shù)最優(yōu)解,即本文提出的基于影響函數(shù)的近場(chǎng)動(dòng)力學(xué)數(shù)值模擬方法,能夠合理、準(zhǔn)確地模擬單晶冰的彈性各向異性行為.

值得補(bǔ)充的是,極地海洋工程中常見(jiàn)的天然冰具有典型的多晶結(jié)構(gòu),其力學(xué)行為與晶體的排列方式、晶粒尺寸等有著密不可分的聯(lián)系.同時(shí),多晶冰的變形及裂紋成核與擴(kuò)展模式不僅與單晶冰的變形有關(guān),還受到晶界滑移等晶間作用機(jī)制的影響與控制.因此,在后續(xù)探索多晶冰的數(shù)值模擬方法時(shí),仍需在單晶冰模擬的基礎(chǔ)上,重點(diǎn)關(guān)注晶體織構(gòu)的建模、單晶之間力的傳遞過(guò)程、斷裂準(zhǔn)則的引入等,以實(shí)現(xiàn)多晶冰變形與破壞各向異性的有效模擬.