二維裂紋傳播問(wèn)題的數(shù)值模擬

1、二維裂紋傳播問(wèn)題的數(shù)值模擬劉劉劉強(qiáng)強(qiáng)強(qiáng)黃黃黃忠忠忠億億億（清華大學(xué)數(shù)學(xué)科學(xué)系北京100084）摘摘摘要要要:研究了二維裂紋傳播問(wèn)題。經(jīng)典的連續(xù)模型對(duì)裂紋傳播描述不夠精確，而原子模型可以彌補(bǔ)這個(gè)缺點(diǎn)，但是完全使用原子模型計(jì)算量太大。為此，在裂紋附近使用原子模型，而在遠(yuǎn)離裂紋的地方使用連續(xù)模型，在交界面上，以吸收高頻波的條件作為出發(fā)點(diǎn)，建立匹配條件求解。這種匹配方法把計(jì)算復(fù)雜度從O(N2)降低到O(N)，而得到的數(shù)值解與完全用原子模型得到的

2、解吻合得很好。關(guān)關(guān)關(guān)鍵鍵鍵詞詞詞：：：原子連續(xù)模型，裂紋，裂紋傳播，匹配條件中中中圖圖圖分分分類(lèi)類(lèi)類(lèi)號(hào)號(hào)號(hào):O241.82NumericalSimulationfCrackPropagationinTwoDimensionalLIUQiangHUANGZhongyi(DepartmentofMathematicalSciencesTsinghuaUniversityBeijing100084China)Abstract:Inthispa

3、perwestudiedthecrackpropagationintwodimensionalcase.Continuummodelsbecomeinaccuratewhendescribingthecrackpropagationinmicroscale.Weusetheatomisticmodelswhicharemeaccuratenearthecracktipwhilecontinuummodelsareusedfarawayt

4、oavoidthehugecomputation.Attheinterfacewepresentanappropriatematchingcondition(suchthatthehighfrequencywavesareabsbed)togettheapproximatesolution.OurmethodcanreducethecomputationalcomplexityfromO(N2)toO(N).Thenumericalre

5、sultsfromourmethodmimictheresultsfromatomisticmodelquitewell.Keywds:atomisticcontinuummodelcrackpropagationmatchingcondition晶體中的裂紋傳播問(wèn)題在材料科學(xué)和許多工程領(lǐng)域中有著重要的應(yīng)用。建立適當(dāng)?shù)哪Ｐ?，研究?基基基金金金項(xiàng)項(xiàng)項(xiàng)目目目：：：國(guó)家自然科學(xué)基金資助項(xiàng)目(10301017)國(guó)家教育部留學(xué)回國(guó)人員科研啟動(dòng)金

6、資助項(xiàng)目；國(guó)家“九七三”重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃資助項(xiàng)目(2005CB321701)作作作者者者簡(jiǎn)簡(jiǎn)簡(jiǎn)介介介：：：劉強(qiáng)，男(漢)，四川，碩士研究生。通通通訊訊訊聯(lián)聯(lián)聯(lián)系系系人人人：：：黃忠億，副教授，Email:ghua.紋的出現(xiàn)條件、傳播的穩(wěn)定性和速度以及傳播過(guò)程中裂紋附近溫度的變化及影響等，是這類(lèi)問(wèn)題的幾個(gè)重點(diǎn)。當(dāng)晶體受到外界的影響后，晶體各個(gè)部分會(huì)產(chǎn)生不同的應(yīng)力，導(dǎo)致晶體發(fā)生形變，甚至出現(xiàn)裂紋。從實(shí)際經(jīng)驗(yàn)來(lái)看，裂紋的產(chǎn)生具有隨機(jī)性，形

7、狀很尖銳而且不規(guī)則，這些都給求解這類(lèi)問(wèn)題帶來(lái)極大的困難。從經(jīng)典力學(xué)原理出發(fā)，利用線(xiàn)彈性方程(組)研究裂紋傳播的規(guī)律，各向同性介質(zhì)的連續(xù)模型為ρ?2u?t2=?u(λ)?(?u)ρf式中：λ為L(zhǎng)ame彈性常數(shù)，為剪切模量。線(xiàn)彈性方程是線(xiàn)性化的近似方程，有著參數(shù)簡(jiǎn)單、求解方便的優(yōu)點(diǎn)。然而裂紋的產(chǎn)生涉及到化學(xué)鍵的斷裂、位錯(cuò)的產(chǎn)生與傳播等非線(xiàn)性因素，裂紋尖端附近受到晶體微觀(guān)結(jié)構(gòu)變化帶來(lái)的影響因此若僅用線(xiàn)彈性方程來(lái)研究裂紋傳播問(wèn)題，計(jì)算結(jié)果很不準(zhǔn)

8、確。分子動(dòng)力學(xué)可以克服上面的缺陷它從牛頓定律出發(fā)，對(duì)晶體中的每個(gè)原子建立運(yùn)動(dòng)方程，即：mri=??j?Φ(|ri?rj|)式中Φ為勢(shì)函數(shù)。分子動(dòng)力學(xué)方程最大的優(yōu)點(diǎn)在于它可以計(jì)算每個(gè)原子的運(yùn)動(dòng)以及產(chǎn)生的影響，但是眾所周知，哪怕是微米尺寸的一塊晶體包含的原子數(shù)目也達(dá)到109以上，再考慮到時(shí)間步長(zhǎng)的限制，要完全求解這樣規(guī)模的方程組幾乎是不可能的。近些年來(lái)，多尺度耦合算法逐漸興起[1?5]。為了彌補(bǔ)線(xiàn)彈性方程在裂紋附近的不準(zhǔn)確性，在這個(gè)區(qū)域用微

9、觀(guān)模型描述；而為了避免分子動(dòng)力學(xué)方程帶來(lái)的巨大的計(jì)算量，在遠(yuǎn)離裂紋、受原子振動(dòng)產(chǎn)生的非線(xiàn)性影響較小的的區(qū)域使用宏觀(guān)模型計(jì)算，這就是這類(lèi)方法的基本思想。耦合方法結(jié)合了兩種模型的優(yōu)點(diǎn)，可行且具有相當(dāng)準(zhǔn)確性。它的難點(diǎn)在于，使用微觀(guān)模型和宏觀(guān)模型交替計(jì)算時(shí)，如何建立耦合邊界處人工邊界條件，使它既保證結(jié)果的可靠性，又能夠應(yīng)用于實(shí)際計(jì)算。這種方法在計(jì)算晶1把這個(gè)模型應(yīng)用到式(2)，得到線(xiàn)性化方程組d2udτ2=18?0a2????????10000

10、0000???(u10?2u00u?10)???14?√340?√34340000???(u?11?2u00u1?1)???14√340√34340000???(u01?2u00u0?1)?????(9)取a→0的極限，得到彈性波方程組d2udτ2=?0????????8140002740000????2u?x2???2740008140000????2u?y2???0272027200000????2u?x?y?????.(10)2吸

11、吸吸收收收邊邊邊界界界條條條件件件現(xiàn)在尋找合適的匹配條件。取一個(gè)包含裂紋尖端的矩形D，把整個(gè)晶體分成MD區(qū)域，交界面和連續(xù)區(qū)域三個(gè)部分(見(jiàn)圖3)。MD區(qū)域內(nèi)部的原子振動(dòng)產(chǎn)生各種頻率的波，通過(guò)交界面?zhèn)鞒?。在連續(xù)模型的宏觀(guān)尺度下，可以通過(guò)網(wǎng)格剖分看到低頻波的傳播；而對(duì)于高頻波，在遠(yuǎn)離裂紋的區(qū)域是可以忽略的，因此適當(dāng)?shù)钠ヅ錀l件應(yīng)該使得高頻波在交界面被完全吸收掉。ccccccccccccccccccccccccccccccccccccccccc

12、ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccD裂紋尖端MD區(qū)區(qū)區(qū)域域域交交交界界界面面面連連連續(xù)續(xù)續(xù)區(qū)區(qū)區(qū)域域域D?DD0圖圖圖3耦合模型的區(qū)域分解對(duì)于原型交界面，文[6]建立了準(zhǔn)確的人工邊界條件，以球諧函數(shù)為基底，把交界面的波用Fourier級(jí)數(shù)表示，可以在得到t時(shí)刻在半徑為R的圓形交界面上的邊界條件：?u?r(Rst)=∞?l=0l?

13、m=?lKil(ulm(Rt))Yml(s)(11)式中i=012s=rx表示空間坐標(biāo)x的單位化向量.對(duì)于方形的交界面，利用極小化反射波能量的方法[1]，可以建立較為精確的邊界條件。為了減少計(jì)算的復(fù)雜度，把交界面上的波分解為低頻和高頻部分：uij=uij?uij。利用連續(xù)模型計(jì)算低頻波，有了差分格式得到差分網(wǎng)格上的數(shù)值結(jié)果Umn，交界面上的每個(gè)原子的低頻波uij由鄰近差分網(wǎng)格的解插值得到。假設(shè)高頻波?u沿著垂直于交界面的方向傳播(見(jiàn)圖3

14、)，那么在D0上，?u沿著x軸的正向傳播從而由式(6)得?2?u?τ2?3λ2?2?u?x2??u?τ=0.(12)因?yàn)楹苄?，省略最后一?xiàng)，選取沿x軸正向傳播的波，在D0上，得到近似的邊界條件：??u?τ=??λ??u?x(13)式中?λ=?3λ2。類(lèi)似，可以推導(dǎo)在D?和D的近似邊界條件，有?????????????D:??u?τ=??λ??u?yD?:??u?τ=?λ??u?yD0:??u?τ=??λ??u?x.(14)由于型III裂

15、紋產(chǎn)生的波主要沿著垂直于交接面的方向傳播，因此上面的吸收邊界條件已經(jīng)是一個(gè)較好的近似。而對(duì)于其他的傳播情形，還可以參考文[12]中的做法，得到更高階的離散邊界條件。3算算算法法法及及及算算算例例例下面給出具體算法。將計(jì)算區(qū)域分解成MD區(qū)域和連續(xù)模型區(qū)域，分別將初始條件離散化，令t0=0。對(duì)k=1K，重復(fù)1)–4)。1)在連續(xù)模型區(qū)域用差分格式求解連續(xù)方程，計(jì)算連續(xù)區(qū)域在tk時(shí)刻的位移Ukmn。2)在MD區(qū)域求解MD方程，得到MD區(qū)域內(nèi)部