含I 型裂紋八次對稱二維準晶動力學問題的數值模擬

王雅楠，李聯和，2，王桂霞，2，胡學佳

（1.內蒙古師范大學 數學科學學院，內蒙古 呼和浩特 010022；2.內蒙古自治區應用數學中心，內蒙古 呼和浩特 010022）

由于具有耐高溫及抗氧化的顯著特性，準晶熱障涂層在發動機等部件中應用廣泛［1-2］。發動機熱障涂層技術是通過對內燃機內部件噴涂熱障涂層使基體金屬起到絕熱屏蔽作用，從而促進燃料的完全燃燒，減少有害氣體的排放及冷卻水的消耗［3］。裂紋、孔洞、位錯及夾雜等缺陷是影響準晶材料強度的主要原因，因此關于準晶斷裂力學問題的研究成為數學、力學以及工程等相關領域科研工作者們關注的焦點。

眾所周知，對準晶斷裂動力學問題的求解比對靜力學問題的求解更困難。帶裂紋的準晶靜力學問題研究成果較多［4-9］，對于準晶動力學問題，即使增加很多限制條件，仍然很難進行解析求解。選擇數值方法進行求解，限制條件少，能更好地模擬實際問題?；诓▌印獢U散模型［10］，祝愛玉等［11］首次將求解晶體材料問題的有限差分法推廣應用到求解準晶彈性動力學問題上，研究了動態應力強度因子隨準晶材料常數的變化規律及裂紋動態擴展過程，取得了有意義的數值結果。王曉芳等［12］研究了帶Ⅰ型Griffith 裂紋的二十面體Al-Pd-Mn 準晶在沖擊載荷作用下的動態響應問題，數值求解了裂紋尖端處應力、位移以及動態應力強度因子，重點分析了相位子場波的傳播及擴散情況。王芳［13］用有限差分法求解十次對稱準晶中的動力學裂紋擴展問題，采用中心Griffith 裂紋與單邊Griffith 裂紋矩形試樣考查了裂紋尖端處的位移、應力、應變場變量隨時間的變化規律。Sladek［14］基于局部Petrov-Galerkin 無網格方法來解決十邊形準晶中的裂紋問題。楊連枝［15］通過將二十面體準晶的彈性邊值問題轉化為等價變分問題，開發了二十面體準晶的靜態彈性有限元算法，給出了具有不同聲子—相位子耦合參數二十面體準晶的解析解與數值解。Zhao 等［16］利用廣義不連續位移邊界積分—微分方程法分析了一維六方熱電彈準晶復合材料中一個任意形狀的三維界面裂紋問題。馬晴等［17］探究了含Ⅱ型單邊裂紋的八次對稱準晶動力學問題，探討了板端加載與裂紋面加載對動態應力強度因子的影響。陳曦［18］分析了含Ⅰ型單邊裂紋的八次對稱準晶動力學問題中不同參數對裂紋尖端處聲子場應力強度因子及相位子場位移分量的影響。上述結果表明了有限差分法對求解準晶動力學問題的可行性和有效性。

波動—電報模型［19］更具一般性，利用該模型研究準晶動力學問題才剛剛開始。基于波動—電報模型，Li［20］研究了十次對稱準晶動力學問題，分別推導出滿足四階偏微分方程的兩個基本解和滿足二階雙曲型擴散方程的一個基本解，周建敏等［21］數值求解了帶Ⅱ型Griffith 裂紋的十次對稱準晶動力學問題。本文將基于波動—電報模型，利用有限差分法對含I 型Griffith 裂紋的八次對稱二維準晶動力學問題進行數值模擬。

1 基本方程

1.1 運動方程

假設以z-軸為二維八次對稱準晶周期方向，xy-平面為準周期平面，裂紋沿著周期方向穿透整個材料的裂紋表面，即?/?z=0。根據準晶彈性理論，有變形幾何方程

其中：ui和wi分別代表聲子場和相位子場位移；εij和wij分別表示聲子場和相位子場應變，x1和x2表示在指標符號系統中的笛卡爾坐標x和y；ui表示uxi；wi，εij以及wij的表示方法與ui類似。及廣義Hooke 定律

其中λ，μ為Lame 常數，R為聲子場—相位子場耦合常數，ki(i=1，2，3)為彈性常數；σxx，σyy，σxy，σyx表示聲子場應力，Hxx，Hyy，Hxy，Hyx表示相位子場應力。把變形幾何方程和廣義Hooke 定律代入到波動—電報模型

得到波動—電報方程組

1.2 初始條件及邊界條件

初始條件為

邊界條件

其中P(t)=P0f(t)為動態加載函數，P0是常數，f(t)是Heaviside函數。

1.3 網格劃分

空間步長取Δx=Δy=h，時間步長取Δt=τ。首先將求解區域[0，L]×[0，H] 向四周延伸h/2，再對計算區域[-h/2，L+h/2]×[-h/2，H+h/2]進行網格劃分，如圖1。

圖1 網格劃分Fig.1 Mesh subdivision

1.4 離散格式

采用中心差分公式對方程組（1）進行離散，有

方程組（2）為方程組（1）的全離散格式。

在邊界處，對垂直于邊界的偏導數用向前差分或向后差分進行離散，對平行邊界的偏導數使用中心差分進行離散。對于角點處使用外推法計算。初始條件的離散類似。

2 誤差估計

定理1在計算區域內部，空間步長取Δx=h1，Δy=h2，時間步長取Δt=τ。假設聲子場和相位子場的位移存在高階導數，離散格式（2）的截斷誤差為

證明由Taylor 展式得

3 穩定性估計

證明使用Fourier 分析法［22-23］對方程組（2）進行穩定性分析，將

當特征值mi(i=1，2，3，4)＜1 時，即滿足

時，方程組（2）是穩定的。

4 數值實例

若無特殊說明取a=2.4 mm，板長L=10 mm，板寬H=25 mm，p0=1 MPa。

由于八次準晶未有實測數據，故取聲子場彈性模量c11=2.343×10-1g·mm-1·μs-2，c12=0.574×10-1g·mm-1·μs-2，相位子場彈性模量K1=1.220×10-1g·mm-1·μs-2，K2=0.240×10-1g·mm-1·μs-2，K3=0.200×10-1g·mm-1·μs-2，R=0.01μ，ρ=4.186×10-3g·mm-3。

裂紋尖端和非裂紋尖端（3.7，0.3）處聲子場位移和相位子場位移隨時間的變化如圖2。

圖2 位移隨時間變化圖Fig.2 Displacement versus time diagram

從圖2 可以看出聲子場位移呈波動趨勢，相位子場位移呈擴散趨勢。裂紋尖端處位移值較大。但由于對稱性，裂紋尖端處y方向位移值為0。

裂紋尖端處聲子場和相位子場的應力分量隨時間變化如圖3。由圖3 可以看出，波到達裂紋尖端的時間t0≈3.35 μs，橫波的速度c1=≈7.48 mm/μs，波傳播速度等于V=≈7.46 mm/μs，約等于橫波的速度，這說明聲子場中的波傳播起主要作用。

圖3 應力分量隨時間變化圖Fig.3 Stress component versus time chart

選取試樣的寬度為25 mm，30 mm，35 mm，聲子場動態無量綱應力強度因子如圖4。對于Ⅰ型對稱準晶，無量綱應力強度因子為，其中

由圖4 可以看出，試樣寬度H為25 mm，應力波到達裂紋尖端的時間t0≈3.35 μs；試樣寬度H為30 mm，應力波到達裂紋尖端的時間t0≈4.01 μs；試樣寬度H為35 mm，應力波到達裂紋尖端的時間t0≈4.68 μs，表明試樣寬度H越大，應力波到達裂紋面的時間越長。

圖4 試樣寬度H 對無量綱化應力強度因子的影響Fig.4 Effect of specimen width（H）on intensity factor of dimensionless stress

板寬H=20 mm 時，有效質量密度ρp和摩擦系數κ對裂紋尖端處相位子場位移wx的影響如圖5。

由圖5 可以看出，當摩擦系數遠大于有效質量密度時，相位子場波的傳播呈指數衰減趨勢，摩擦系數對波的傳播影響占主導。隨著摩擦系數逐漸減小，相位子場波的傳播呈指數衰減且出現波動趨勢，有效質量密度的改變將會影響相位子場的位移量。

圖5 有效質量密度ρp 和摩擦系數κ 對相位子場位移wx 的影響Fig.5 The effect of effective mass density（ρp）and friction coefficient（κ）on the phase field displacement（wx）

5 總結

本文基于波動—電報模型，求解了帶Ⅰ型裂紋的八次對稱準晶動力學問題，利用有限差分法構造了混合精度的離散格式，該格式的截斷誤差為利用Fourier 分析法給出穩定條件。數值實例表明：

（1）裂紋尖端處相較于非尖端處位移值較大；

（2）隨著板寬H的增加，應力波到達裂紋面所需的時間也在增加；

（3）摩擦系數遠大于有效質量密度時，相位子場波的傳播呈指數衰減趨勢，摩擦系數對波的傳播影響占主導，隨著摩擦系數逐漸減小，相位子場波的傳播呈指數衰減且出現波動趨勢，有效質量密度的改變將會影響相位子場的位移量。