【摘 要】 弹性非均匀介质的弹性应力场模拟是工程力学领域的重要应用。本文以相场微弹性理论为基础，改进了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%的条件下将计算规模缩小为原来的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.
推荐访问： 力场 网格 介质 算法 均匀