彈性非均勻介質(zhì)應(yīng)力場(chǎng)動(dòng)網(wǎng)格算法改進(jìn)

顧驍

【摘 要】 彈性非均勻介質(zhì)的彈性應(yīng)力場(chǎng)模擬是工程力學(xué)領(lǐng)域的重要應(yīng)用。本文以相場(chǎng)微彈性理論為基礎(chǔ)，改進(jìn)了Chen的彈性非均勻體系的動(dòng)網(wǎng)格法，提升了30%左右的計(jì)算效率。

【關(guān)鍵詞】 彈性非均勻 相場(chǎng)微彈性理論 動(dòng)網(wǎng)格法

1 導(dǎo)論

非均勻介質(zhì)的彈性應(yīng)力場(chǎng)研究一直是工程力學(xué)領(lǐng)域的一個(gè)重要問題。絕大多數(shù)的工程材料系統(tǒng)都是彈性非均勻的，即使單相多晶中每一個(gè)晶粒都是完美晶體，整個(gè)材料也是彈性非均勻的，因?yàn)閺恼w坐標(biāo)系看，依據(jù)不同位相放置的單個(gè)晶粒的彈性模量張量也是不同的。而多相系統(tǒng)——比如多層膜、功能梯度材料等常見的人造復(fù)合工程材料——不僅是彈性非均勻，更是晶體結(jié)構(gòu)非均勻的。存在晶體缺陷的系統(tǒng)也可以看作一類重要的彈性非均勻系統(tǒng)。人們希望建立一個(gè)存在外加應(yīng)力的三維分析模型來研究上述彈性非均勻系統(tǒng)。但是其分析方法一直受到數(shù)學(xué)理論的限制[1，2]。

相場(chǎng)微彈性理論是一種彈性非均勻場(chǎng)模擬計(jì)算的相場(chǎng)分析方法[3]。它可以計(jì)算合金條幅分解[4]、位錯(cuò)附近的溶質(zhì)偏析[5]、馬氏體相變[6]等多種動(dòng)力學(xué)過程的彈性應(yīng)力場(chǎng)演化。相場(chǎng)方法善于模擬復(fù)雜結(jié)構(gòu)的晶粒、缺陷等，但往往需要很大的計(jì)算量。非均勻計(jì)算網(wǎng)絡(luò)可以很大程度地降低計(jì)算成本，通過網(wǎng)格非均勻化，將計(jì)算資源集中在應(yīng)力變化較大的區(qū)域，可以在保證計(jì)算結(jié)果精度的條件下很好地節(jié)約計(jì)算資源。Chen等人給出了相場(chǎng)模型計(jì)算單元非均勻化算法的修正方程，并將這一方法稱為動(dòng)網(wǎng)格法[7]。本文在他的算法基礎(chǔ)上，建立了改進(jìn)的非均勻介質(zhì)相場(chǎng)微彈性理論的動(dòng)網(wǎng)格算法。

2 非均勻介質(zhì)相場(chǎng)微彈性理論的動(dòng)網(wǎng)格算法及其改進(jìn)

在計(jì)算模型中，計(jì)算坐標(biāo)和物理坐標(biāo)之間存在映射，當(dāng)時(shí)，計(jì)算單元是均勻的。彈性介質(zhì)的應(yīng)力應(yīng)變關(guān)系、幾何方程和應(yīng)力平衡方程都基于物理坐標(biāo)，將這三個(gè)方程結(jié)合得到：

（1）

其中，是體系平均彈性模量，是局部彈性模量，是位移矢量，是應(yīng)變張量，是無應(yīng)力應(yīng)變張量，代表由于相變、晶體缺陷等因素引起的塑性應(yīng)變。設(shè)，取弗羅貝尼烏斯范數(shù)作為比例因子，保證對(duì)于任意矩陣都有，將公式（1）改寫為一個(gè)迭代方程：

（2）

其中下標(biāo)NStep和N+1 Step的變量分別代表第N和N+1迭代步的值。對(duì)公式（2）進(jìn)行傅里葉變換，得到：

（3）

所以，修正的迭代算法通過在倒空間迭代位移矢量得到位移的解，再計(jì)算系統(tǒng)應(yīng)變場(chǎng)。這一迭代算法已經(jīng)實(shí)現(xiàn)了計(jì)算網(wǎng)格的非均勻化，為了進(jìn)一步實(shí)現(xiàn)網(wǎng)格在每一次迭代后的重新劃分，首先需要確定局部區(qū)域分配的計(jì)算格點(diǎn)數(shù)目的權(quán)重。結(jié)合應(yīng)力場(chǎng)的判定公式，可以以應(yīng)力場(chǎng)梯度的平方作為權(quán)重。其次，需要確定如何將舊網(wǎng)格的數(shù)據(jù)導(dǎo)入新網(wǎng)格，三次樣條插值等插值方法是很好的選擇。通過合理的計(jì)算網(wǎng)格邊長調(diào)整方案，非均勻化計(jì)算網(wǎng)格的算法可以在應(yīng)變結(jié)果偏差小于1%的條件下將計(jì)算規(guī)?？s小為原來的10到40倍，這一倍數(shù)與模型具體的計(jì)算內(nèi)容有關(guān)。

在此基礎(chǔ)上，我們得到了改進(jìn)的非均勻計(jì)算單云修正方程。注意到公式（3）的迭代過程中需要計(jì)算應(yīng)變張量。而應(yīng)變張量需要通過位移矢量求偏導(dǎo)得到。該求偏導(dǎo)過程增加數(shù)值求偏導(dǎo)或者額外的傅里葉正逆變換步驟。如果用位移對(duì)于物理坐標(biāo)的梯度作為迭代變量，將公式（3）改寫為：

（4）

那么應(yīng)變張量可以通過迭代變量的加減運(yùn)算得到，所以這種替代可以節(jié)約計(jì)算成本。

3 改進(jìn)效果分析

改進(jìn)算法存在計(jì)算速度上的優(yōu)勢(shì)，也存在數(shù)據(jù)存儲(chǔ)空間上的劣勢(shì)。在計(jì)算速度方面，可以進(jìn)行定性的比較，假設(shè)計(jì)算單元一共有N個(gè)計(jì)算單元，通過傅里葉正逆變換來計(jì)算位移對(duì)于物理坐標(biāo)的梯度，從而計(jì)算應(yīng)變張量，并忽略數(shù)組傅里葉變換的規(guī)模效應(yīng)。注意到傅里葉正逆變換是迭代的主要計(jì)算成本。如果N個(gè)標(biāo)量進(jìn)行一次傅里葉正變換或逆變換需要t的時(shí)間，那么利用公式（3）進(jìn)行一次迭代的過程中，通過位移計(jì)算應(yīng)變需要，公式（3）中的傅里葉變換需要，總用時(shí)為。利用公式（4）進(jìn)行迭代的時(shí)間是，可見新的迭代算法對(duì)原算法的加速比為。表1是改進(jìn)算法和原算法分別計(jì)算Cu/Ni納米多層膜中單根位錯(cuò)產(chǎn)生的應(yīng)力場(chǎng)的迭代過程的實(shí)際計(jì)算時(shí)間，可見實(shí)測(cè)的加速比在1.2到1.4之間。所以，替代算法具有加速優(yōu)勢(shì)。在計(jì)算存儲(chǔ)空間方面，替代算法中的每個(gè)計(jì)算單元需要始終儲(chǔ)存位移偏導(dǎo)這個(gè)九維數(shù)組，而原算法只在部分時(shí)間會(huì)達(dá)到這一存儲(chǔ)量，如果適當(dāng)調(diào)整循環(huán)運(yùn)算順序并犧牲部分計(jì)算效率，可以減低存儲(chǔ)量，這對(duì)于計(jì)算規(guī)模相對(duì)較大的相場(chǎng)微彈性模型是有意義的。

4 結(jié)語

本文在Chen的非均勻彈性場(chǎng)動(dòng)網(wǎng)格算法的基礎(chǔ)上，重新選擇了迭代變量，并修改了迭代算法，使計(jì)算效率提升了30%左右。這一算法改進(jìn)對(duì)于以相場(chǎng)微彈性理論為基礎(chǔ)的晶體缺陷和相變過程的數(shù)值模擬具有重要意義，特別是對(duì)需要計(jì)算規(guī)模較大的研究復(fù)雜構(gòu)型的相場(chǎng)模型，可以大大減少計(jì)算時(shí)間。結(jié)合高性能計(jì)算和計(jì)算網(wǎng)絡(luò)非均勻化本身的效果，可以獲得理想的加速比。

參考文獻(xiàn)

[1]Demir I，Hirth J P，Zbib H M.The extended stress field around a cylindrical crack using the theory of dislocation pile-ups[J].International journal of engineering science，1992，30（7）：829-845.

[2]Leo P H， Lowengrub J S， Jou H J. A diffuse interface model for microstructural evolution in elastically stressed solids[J].Acta materialia，1998，46（6）： 2113-2130.

[3]Wang Y U，Jin Y M，Khachaturyan A G.Phase field microelasticity theory and modeling of elastically and structurally inhomogeneous solid[J].Journal of Applied Physics，2002，92（3）：1351-1360.

[4]Wang Y，Banerjee D，Su C C，et al.Field kinetic model and computer simulation of precipitation of L12 ordered intermetallics from fcc solid solution[J].Acta materialia， 1998，46（9）：2983-3001.

[5]Hu S Y，Chen L Q.Solute segregation and coherent nucleation and growth near a dislocation—a phase-field model integrating defect and phase microstructures[J].Acta materialia，2001，49（3）：463-472.

[6]Jin Y M，Artemev A，Khachaturyan A G. Three-dimensional phase field model of low-symmetry martensitic transformation in polycrystal：simulation of ζ′2 martensite in AuCd alloys[J].Acta materialia，2001，49（12）：2309-2320.

[7]Feng W M，Yu P，Hu S Y，et al.A Fourier spectral moving mesh method for the Cahn Hilliard equation with elasticity[J].Commun Comput Phys，2009，5：582-599.