任意形狀熱夾雜位移場(chǎng)的三角形單元離散算法1)

侯佳卉 李璞 黎江林 金曉清

(重慶大學(xué)機(jī)械傳動(dòng)國(guó)家重點(diǎn)實(shí)驗(yàn)室,重慶 400040)

引言

隨著科學(xué)技術(shù)水平的日益發(fā)展,無(wú)論在航空航天、武器裝備領(lǐng)域,還是在汽車(chē)零部件等普通工程應(yīng)用領(lǐng)域,對(duì)于材料的綜合性能要求越來(lái)越高[1].大多數(shù)復(fù)合材料在制造裝配過(guò)程中,其內(nèi)部普遍存在熱膨脹、塑性變形、錯(cuò)配、位錯(cuò)、相變等非彈性變形,美國(guó)學(xué)者M(jìn)ura[2]將與之引起的相應(yīng)非彈性應(yīng)變統(tǒng)稱(chēng)為“本征應(yīng)變”.如果本征應(yīng)變只存在于基體的某局部區(qū)域,力學(xué)家Eshelby[3]將該區(qū)域命名為“夾雜”.Eshelby 夾雜問(wèn)題[4]是固體力學(xué)中的一個(gè)經(jīng)典課題,我國(guó)力學(xué)工作者在夾雜及其相關(guān)問(wèn)題上的研究也取得豐碩成果[5-10].

夾雜問(wèn)題在固體力學(xué)及材料科學(xué)領(lǐng)域具有廣泛應(yīng)用.在真實(shí)的工程零部件材料中,夾雜形狀各異,不規(guī)則形狀?yuàn)A雜的解析解通常不易獲得.1961 年Eshelby[3]完成了無(wú)限大基體材料中橢球夾雜問(wèn)題的開(kāi)創(chuàng)性工作,并提出經(jīng)典的Eshelby 張量.此后,橢球/橢圓形夾雜得到國(guó)內(nèi)外學(xué)者的廣泛研究[11-15].Ju和Sun[11]基于Eshelby 張量,通過(guò)引入虛擬橢球外單位法向矢量提出了外場(chǎng)點(diǎn)應(yīng)變的一種類(lèi)張量表示公式.基于Ju 和Sun 的成果,Jin 等[12]推導(dǎo)了二維平面應(yīng)變下橢圓柱夾雜外部彈性場(chǎng)的封閉解.隨后Jin 研究組[13]給出了三維橢球夾雜的彈性場(chǎng)完備解,系統(tǒng)解決了夾雜內(nèi)外場(chǎng)的位移、位移梯度、應(yīng)力、應(yīng)變等四類(lèi)Eshelby 張量顯式解析解;通過(guò)三維向二維問(wèn)題的轉(zhuǎn)化,可獲得二維橢圓柱內(nèi)外場(chǎng)位移解析解[14].

應(yīng)該指出,帶有尖銳角點(diǎn)的夾雜往往具有與橢球/橢圓迥異的性能,其相關(guān)研究也受到研究者廣泛關(guān)注.例如陶瓷基復(fù)合材料中的SiC 纖維橫截面呈六邊形形狀[16]、超導(dǎo)復(fù)合材料中的共晶呈多邊形[17]、基體中的量子點(diǎn)通常為多面體形狀[18],故多邊形夾雜更適合模擬此類(lèi)問(wèn)題.針對(duì)多邊形夾雜的研究,Chiu[19]求解了矩形夾雜在彈性層中受到均勻本征應(yīng)變的應(yīng)力場(chǎng)問(wèn)題,并獲得其封閉形式解.Rodin[20]提出了二維平面多邊形和三維空間多面體夾雜Eshelby 張量的算法.Nozaki 和Taya[21]通過(guò)將二維多邊形面積分轉(zhuǎn)化為沿單位圓的線積分,研究了具有均勻本征應(yīng)變的任意凸多邊形夾雜的彈性場(chǎng).基于解析延拓和保角映射,Ru[22]提出了任意形狀?yuàn)A雜在平面/半平面的解析方法,然而在映射函數(shù)包含無(wú)窮項(xiàng)時(shí),該方法在實(shí)際應(yīng)用中需截?cái)酁橛邢揄?xiàng),因而只能給出近似解,也無(wú)法討論在夾雜角點(diǎn)處的奇異性.周青華等[23]提出一種基于矩形單元的數(shù)值化計(jì)算方法,借助數(shù)值等效夾雜方法,解決了二維平面任意形狀雜質(zhì)問(wèn)題.

即便是二維問(wèn)題,夾雜模型通常涉及復(fù)雜的計(jì)算公式.其中各向同性受均勻熱膨脹本征應(yīng)變(文獻(xiàn)中稱(chēng)之為熱夾雜) 的問(wèn)題因其相關(guān)計(jì)算公式相對(duì)較簡(jiǎn)潔而廣受研究者關(guān)注[24-26].熱夾雜模型在實(shí)際工程與研究中也有重要應(yīng)用,例如可以用來(lái)分析量子線[27-28]與基體之間的晶格失配及應(yīng)變松弛,還可用于模擬孔隙壓力,處理相關(guān)的地質(zhì)學(xué)問(wèn)題等[29-30].

利用熱夾雜模型,Faux 等[31]基于格林函數(shù)方法,通過(guò)圍道積分獲得了線單元和圓弧單元的應(yīng)力解.鑒于位移解是完備彈性場(chǎng)解的一個(gè)重要組成部分[14],并且位移場(chǎng)解也可以作為通過(guò)求導(dǎo)計(jì)算應(yīng)變,進(jìn)而通過(guò)胡克定律得出應(yīng)力的計(jì)算起點(diǎn),因而在彈性力學(xué)分析中具有重要的地位.本文的研究可視為對(duì)Faux 等[31]的必要補(bǔ)充.Nakasone 等[32]提出了基于三角形單元的數(shù)值等效夾雜方法(NEIM),但是他們的算法在三角形夾雜頂點(diǎn)處存在應(yīng)力奇異,給數(shù)值編程計(jì)算帶來(lái)處理上的困難.相對(duì)應(yīng)力/應(yīng)變結(jié)果,本文分析表明位移場(chǎng)在多邊形角點(diǎn)處不存在奇異性,因而具有理論分析上良好的光滑連續(xù)性與編程計(jì)算上的數(shù)值穩(wěn)定性.

本文提出二維平面下任意形狀熱夾雜位移解的三角形單元離散算法.借助圍道積分,給出了夾雜邊界線單元所引起位移貢獻(xiàn)的解析解.在均勻分布熱本征應(yīng)變情況下:本文研究可直接得出任意多邊形夾雜的位移場(chǎng)封閉解,而針對(duì)不規(guī)則形狀?yuàn)A雜只需對(duì)其邊界進(jìn)行離散,計(jì)算簡(jiǎn)便且效率提高顯著.在現(xiàn)有文獻(xiàn)中,大多數(shù)算法是采用均勻的矩形單元網(wǎng)格[23,33]離散夾雜區(qū)域,因而即便對(duì)于均勻本征應(yīng)變情況,也需進(jìn)行二維離散,造成計(jì)算浪費(fèi).本文提出的三角形單元離散法可用于非均勻本征熱應(yīng)變情況.當(dāng)受非均勻分布本征應(yīng)變時(shí),夾雜區(qū)域可數(shù)值離散為足夠精細(xì)的三角形單元,進(jìn)而利用疊加原理獲得原夾雜所產(chǎn)生位移場(chǎng)的數(shù)值解.

1 任意形狀熱夾雜位移解數(shù)值化方法

1.1 均勻本征應(yīng)變下熱夾雜位移解

本文所有變量均定義在笛卡爾坐標(biāo)系中.假設(shè)在平面應(yīng)變條件下,無(wú)限大平面基體內(nèi)嵌有任意形狀?yuàn)A雜Ω,其內(nèi)部有均勻分布的本征應(yīng)變(圖1).夾雜與基體為各向同性,且具有相同彈性常數(shù)Cijkl

式中,μ為剪切模量,ν 為泊松比,δij為Kronecker delta符號(hào).

夾雜問(wèn)題可以借助“激勵(lì)?響應(yīng)”機(jī)制[34]進(jìn)行求解.在激勵(lì)點(diǎn)X′(ξ,η)施加沿j方向單位集中力時(shí),響應(yīng)點(diǎn)X(x,y)在i方向產(chǎn)生的位移定義為格林函數(shù)Gij(X?X′).針對(duì)二維全平面問(wèn)題,為方便計(jì)算其具體表達(dá)式[2]可寫(xiě)成如下形式

圖1 分布有均勻本征應(yīng)變的夾雜示意圖Fig.1 Schematic diagram of inclusion with uniformly distributed eigenstrains

式中,l是響應(yīng)點(diǎn)X(x,y) 與激勵(lì)點(diǎn)X′(ξ,η) 的相對(duì)距離:l2=其中l(wèi)1=x?ξ,l2=y?η.利用體力法[2],二維全平面夾雜問(wèn)題的位移解可以表示為

為了便于表示及計(jì)算機(jī)編程,引入基于Voigt 記號(hào)[14]的矩陣形式:位移和本征應(yīng)變分別用[u1,u2]T和表示.本文考慮本征應(yīng)變均勻分布的熱夾雜(即0),此時(shí)式(3)為

應(yīng)用格林定理,式(4)中的面積分可轉(zhuǎn)化為沿閉合曲線邊界?Ω(如圖1 所示)的圍道積分

為進(jìn)一步求解式(6),將夾雜邊界數(shù)值離散為線單元,如圖2 所示,設(shè)一個(gè)典型線單元L的兩個(gè)端點(diǎn)分別為(x1,y1),(x2,y2),任意位置(可以位于夾雜內(nèi)部或外部)響應(yīng)點(diǎn)記為(x0,y0).經(jīng)過(guò)一系列代數(shù)運(yùn)算和簡(jiǎn)化推導(dǎo),式(6)右端在該線單元上的積分結(jié)果可寫(xiě)成如下形式的封閉解

圖2 夾雜邊界離散示意圖Fig.2 Discretization of the inclusion boundary by line elements

任意形狀?yuàn)A雜的邊界都可用多邊形進(jìn)行數(shù)值離散近似,進(jìn)而基于線單元解的迭加可獲得最終的解析解,從而解決任意形狀?yuàn)A雜產(chǎn)生的位移影響.本文封閉解以線單元的頂點(diǎn)坐標(biāo)作為輸入?yún)⒘?提供了一種高效且便于計(jì)算機(jī)上實(shí)現(xiàn)的數(shù)值計(jì)算方法.

1.2 頂點(diǎn)奇異性分析

現(xiàn)有文獻(xiàn)[20,32]指出,應(yīng)變/應(yīng)力解在多邊形頂點(diǎn)處具有奇異性,從而帶來(lái)實(shí)際應(yīng)用中數(shù)值處理上的困難,此類(lèi)由于奇異性造成的計(jì)算編程的數(shù)值穩(wěn)定性問(wèn)題在文獻(xiàn)[34]中有相關(guān)討論.作為對(duì)比,現(xiàn)考察本文線單元所貢獻(xiàn)位移解的奇異性.當(dāng)響應(yīng)點(diǎn)(x0,y0)無(wú)限靠近線單元L的端點(diǎn)(x2,y2)時(shí),線單元封閉解式(7)簡(jiǎn)化為

式(8) 結(jié)果表示在線單元兩個(gè)頂點(diǎn)處位移結(jié)果是有界且連續(xù)的,說(shuō)明位移解在線單元端點(diǎn)處不產(chǎn)生奇異性(這與應(yīng)變/應(yīng)力解迥然不同).此外,與本文算法相比,常見(jiàn)的限元方法在處理角點(diǎn)(幾何不連續(xù))處的位移通常需要極其精密的網(wǎng)格,且往往精度較差.本文算法基于線單元,數(shù)值離散實(shí)施較易,在計(jì)算分析中可帶來(lái)極大便利.

1.3 三角形夾雜單元離散算法

在非均勻溫升問(wèn)題中,夾雜內(nèi)部受非均勻分布的熱本征應(yīng)變.基于本文均勻熱本征應(yīng)變下的線單元貢獻(xiàn)解,可進(jìn)一步提出如下數(shù)值算法:將待求解的夾雜區(qū)域離散為三角形基本單元,見(jiàn)圖3.當(dāng)三角形單元足夠精細(xì)時(shí),則可近似假定每個(gè)單元內(nèi)部分布均勻本征應(yīng)變(可取其重心點(diǎn)處本征應(yīng)變?yōu)橛?jì)算值).通過(guò)1.1 節(jié)結(jié)果利用線單元解的迭加,可以得出三角形單元的位移貢獻(xiàn).進(jìn)而基于夾雜問(wèn)題的線性性質(zhì),通過(guò)離散的三角形單元解的迭加,可求得受任意分布熱本征應(yīng)變?nèi)我庑螤钇矫鎶A雜問(wèn)題位移場(chǎng)的數(shù)值解.在處理此類(lèi)問(wèn)題時(shí),本文方法只需在夾雜域內(nèi)進(jìn)行離散,而有限元方法則要同時(shí)對(duì)無(wú)限大基體進(jìn)行離散,造成計(jì)算量的增加和時(shí)間成本的浪費(fèi).

圖3 夾雜區(qū)域三角形離散Fig.3 Triangular discretization of inclusion region

在實(shí)際計(jì)算中,夾雜區(qū)域的三角形離散可結(jié)合有限元軟件如ABAQUS 完成.在ABAQUS 中提取夾雜模型離散后的三角形單元節(jié)點(diǎn)信息,并將其作為輸入?yún)?shù)引入程序,利用三角形單元離散算法計(jì)算待求目標(biāo)場(chǎng)點(diǎn)處的位移值.下文中所有數(shù)值計(jì)算均通過(guò)FORTRAN 語(yǔ)言編程實(shí)現(xiàn).

2 數(shù)值化方法的實(shí)現(xiàn)及驗(yàn)證

2.1 均勻本征應(yīng)變下熱夾雜解驗(yàn)證

如圖4 所示,利用無(wú)限平面介質(zhì)中的橢圓熱夾雜模型作為實(shí)例驗(yàn)證位移解的有效性.如1.1 節(jié)所述,均勻本征應(yīng)變下只需離散邊界,通過(guò)對(duì)比分析(見(jiàn)表1),離散線單元數(shù)m取50 時(shí),既能保證計(jì)算精度(相對(duì)誤差大致處于0.2%范圍),又有較好的運(yùn)行效率,所以本節(jié)算例模型邊界均離散為50 個(gè)線單元.設(shè)a1,a2分別為橢圓形的長(zhǎng)、短半軸,定義形狀參數(shù)為γ=a2/a1.本文擬采取Ti-6Al-4V 材料為研究對(duì)象,其楊氏模量E為110GPa,泊松比ν 為0.34,取ε0=1.0×10?3,驗(yàn)證夾雜對(duì)基體位移場(chǎng)的擾動(dòng)影響.

圖4 橢圓形夾雜模型Fig.4 Elliptical inclusion model

表1 數(shù)值解求解得到位移分量與解析解的平均相對(duì)誤差Table 1 Average relative errors of displacement components obtained by the numerical method compared with those by analytical solution

圖5 給出了Ti-6Al-4V 材料在長(zhǎng)半徑a1=1 不變的情況下,形狀參數(shù)分別為γ=0.25,0.5,1 時(shí),橢圓夾雜模型引起的沿圖4 所示目標(biāo)線的位移,最終的位移結(jié)果通過(guò)除以u(píng)0=a1ε0×103/[2π(1 ?ν)]實(shí)現(xiàn)無(wú)量綱化.將數(shù)值解與解析解[15]進(jìn)行對(duì)比分析.當(dāng)場(chǎng)點(diǎn)(x,y)位于橢圓內(nèi)部時(shí)(內(nèi)場(chǎng)解),位移為

圖5 Ti-6Al-4V 材料不同橢圓形夾雜沿目標(biāo)線的位移Fig.5 The displacement along the target line for various Ti-6Al-4V inclusions in elliptical shape

而外場(chǎng)位移解析解為

式(10)中J11,J12,J21,J22,ρ1,ρ2的具體表達(dá)式見(jiàn)附錄.由圖觀察可知,含有均勻本征應(yīng)變的熱夾雜引起夾雜附近基體位移場(chǎng)的擾動(dòng),且位移場(chǎng)在夾雜內(nèi)部呈線性分布,與三維的橢球夾雜位移解[13]情況一致.因?yàn)閼?yīng)變可由位移求導(dǎo)獲得,由之可以推斷橢圓夾雜內(nèi)部總應(yīng)變?yōu)榫鶆蚍植?此即Eshelby 夾雜的經(jīng)典結(jié)論[3].隨著形狀參數(shù)γ 的增大,產(chǎn)生的位移影響增大,在內(nèi)外場(chǎng)邊界處,位移分量連續(xù);在夾雜外部,位移分量隨x軸坐標(biāo)的增大而減小,并與解析解結(jié)果吻合.

進(jìn)一步,在形狀參數(shù)γ 為0.5 的情況下,將現(xiàn)有方法求得的結(jié)果與解析解得到的沿共焦橢圓(圖6)產(chǎn)生的位移結(jié)果進(jìn)行比較,位移結(jié)果如圖7 所示.

圖6 橢圓形狀參數(shù)為0.5 時(shí)內(nèi)外共焦橢圓模型Fig.6 Inner and outer confocal elliptical models with aspect ratio γ=0.5

圖7 Ti-6Al-4V 材料在夾雜橢圓形狀參數(shù)為0.5 時(shí)內(nèi)外共焦橢圓位移解Fig.7 Displacement solution of internal and external confocal ellipse of Ti-6Al-4V material with aspect ratio γ=0.5

由于位移u1和u2分別關(guān)于x軸對(duì)稱(chēng)/反對(duì)稱(chēng),此處只給出θ 取0～π 范圍的位移結(jié)果對(duì)比.通過(guò)觀察可知,數(shù)值解與位移解兩種方法較好的保持了一致性,從而驗(yàn)證了數(shù)值法的有效性.

為進(jìn)行定量誤差分析,定義平均相對(duì)誤差如下

2.2 非均勻本征應(yīng)變下熱夾雜解驗(yàn)證

為說(shuō)明數(shù)值法對(duì)于非均勻本征應(yīng)變下熱夾雜問(wèn)題的處理能力.考慮含于無(wú)限大平面內(nèi)的圓形熱夾雜.如圖8 所示,其半徑為a,以該圓形夾雜的中心為坐標(biāo)原點(diǎn),夾雜內(nèi)部分布有線性本征應(yīng)變?chǔ)?ξ.最終的位移結(jié)果通過(guò)u0=aε0/[2π(1 ?ν)]無(wú)量綱化,圖9 給出了本文數(shù)值解與FEM 求解得到的圓形熱夾雜模型沿圖8 所示目標(biāo)線方向位移分布圖,由圖可見(jiàn),在外場(chǎng)時(shí),夾雜產(chǎn)生的位移影響隨著x軸的增大而減小,且逐漸趨于0.兩種方法所得結(jié)果吻合良好,驗(yàn)證了本文數(shù)值方法處理非均勻分布本征應(yīng)變的有效性.

本文算法對(duì)不規(guī)則形狀熱夾雜同樣適用.如圖10(a) 所示,存在于無(wú)限大平面的半橢圓熱夾雜,其橢圓形狀參數(shù)γ=0.5,初始本征應(yīng)變?cè)O(shè)為最終的位移結(jié)果通過(guò)u0=a1ε0/[2π(1 ?ν)]無(wú)量綱化,目標(biāo)線取x軸.對(duì)于不規(guī)則形狀?yuàn)A雜問(wèn)題,通常不易得到其解析解,再加上本征應(yīng)變?yōu)榫€性分布,其解析解更加難以獲得,而本文所提出的數(shù)值計(jì)算方法可有效解決這些難點(diǎn).把夾雜區(qū)域用三角形單元進(jìn)行離散,計(jì)算網(wǎng)格見(jiàn)圖10(a).利用本文算法獲取的位移結(jié)果與FEM 比較,如圖10(b)所示.通過(guò)觀察可知位移分量在夾雜邊界處連續(xù),且隨著外場(chǎng)點(diǎn)逐漸遠(yuǎn)離中心點(diǎn),位移分量漸漸趨于0.數(shù)值計(jì)算方法與FEM 結(jié)果保持了良好的一致性,由此驗(yàn)證了數(shù)值計(jì)算方法處理非均勻本征應(yīng)變下不規(guī)則形狀熱夾雜問(wèn)題的能力.

圖8 圓形夾雜模型示意圖Fig.8 Schematic diagram of circular inclusion model

圖9 圓形熱夾雜在線性本征應(yīng)變下的位移驗(yàn)證Fig.9 Verification of the displacement for thermal inclusion with linear eigenstrain

圖10 不同方法求解半橢圓形熱夾雜模型所得結(jié)果比較Fig.10 Comparative study on thermal inclusion of semi-elliptic shape solved by different methods

3 結(jié)論

提出一種數(shù)值化計(jì)算方法,用以處理二維任意形狀?yuàn)A雜位移問(wèn)題.主要得到如下結(jié)論:(1) 推導(dǎo)得到基本線單元的位移解形式,與應(yīng)力/應(yīng)變解不同,本文所提出的位移解在多邊形夾雜頂點(diǎn)處并無(wú)奇異性;(2)對(duì)于均勻本征應(yīng)變下任意形狀熱夾雜問(wèn)題,通過(guò)邊界線離散化,可基于線單元解求解位移場(chǎng);對(duì)于非均勻本征應(yīng)變問(wèn)題,可將夾雜區(qū)域離散為三角形單元進(jìn)行求解;(3)利用橢圓、圓形與半橢圓形熱夾雜的算例驗(yàn)證本文所提數(shù)值計(jì)算方法對(duì)于任意形狀熱夾雜問(wèn)題的處理能力,結(jié)果得到了較好的驗(yàn)證.

附 錄

附錄式(10)相關(guān)參數(shù)表達(dá)式

設(shè)a1,a2分別為橢圓形的長(zhǎng)、短半軸,當(dāng)場(chǎng)點(diǎn)X(x,y)位于橢圓外部時(shí)(即>1)