彈性非均勻介質應力場動網格算法改進

顧驍

【摘 要】 彈性非均勻介質的彈性應力場模擬是工程力學領域的重要應用。本文以相場微彈性理論為基礎，改進了Chen的彈性非均勻體系的動網格法，提升了30%左右的計算效率。

【關鍵詞】 彈性非均勻 相場微彈性理論 動網格法

1 導論

非均勻介質的彈性應力場研究一直是工程力學領域的一個重要問題。絕大多數的工程材料系統都是彈性非均勻的，即使單相多晶中每一個晶粒都是完美晶體，整個材料也是彈性非均勻的，因為從整體坐標系看，依據不同位相放置的單個晶粒的彈性模量張量也是不同的。而多相系統——比如多層膜、功能梯度材料等常見的人造復合工程材料——不僅是彈性非均勻，更是晶體結構非均勻的。存在晶體缺陷的系統也可以看作一類重要的彈性非均勻系統。人們希望建立一個存在外加應力的三維分析模型來研究上述彈性非均勻系統。但是其分析方法一直受到數學理論的限制[1，2]。

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

2 非均勻介質相場微彈性理論的動網格算法及其改進

在計算模型中，計算坐標和物理坐標之間存在映射，當時，計算單元是均勻的。彈性介質的應力應變關系、幾何方程和應力平衡方程都基于物理坐標，將這三個方程結合得到：

（1）

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

（2）

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

（3）

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

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

（4）

那么應變張量可以通過迭代變量的加減運算得到，所以這種替代可以節約計算成本。

3 改進效果分析

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

4 結語

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

參考文獻

[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.