文章编号：1004-0609(2007)01-0124-07
有限元与Monte Carlo方法耦合的冷轧纯铝板再结晶模拟
申孝民1，关小军1，张继祥1，刘运腾1，麻晓飞1，赵宪明2
(1. 山东大学 材料科学与工程学院，济南 250061；
2. 东北大学 轧制技术及连轧自动化国家重点实验室，沈阳 110004)
摘 要：
利用有限元软件ANSYS模拟的冷轧铝板应力场及其相应的储能场，在假设的介观非均匀储能场基础 上，考虑冷轧铝板的宏观能量场的非均匀性，实现了金属变形的有限元方法与Monte Carlo 再结晶模拟方法的耦合，有效地模拟非均匀储能场基础上的冷轧铝板再结晶过程，结合现有的理论与实验结果，比较和验证了两个极限储能部位的介观结构拓扑组织和再结晶动力学参数，得到较为理想的结果。
关键词：
中图分类号：TG 111 文献标识码：A
Coupling of FEM with Monte Carlo for simulating recrystallization in cold rolling pure aluminum sheet
SHEN Xiao-min1, GUAN Xiao-jun1, ZHANG Ji-xiang1, LIU Yun-teng1, MA Xiao-fei1, ZHAO Xian-ming2
(1. School of Material Science and Engineering, Shandong University, Ji’nan 250061, China;
2. State Key Laboratory of Rolling and Automation, Northeastern University, Shenyang 110004, China)
Abstract: The finite element method (FEM) and Monte Carlo method (MC) were coupled to simulate the microstructure evolution of inhomogeneously deformed aluminum sheet after annealing. Finite element method was used to calculate the stored energy distribution of the cold rolling aluminum sheet. On the assumption of the mesoscale non-uniform energy distribution, taking macro-scale non-uniform energy distribution into account, the stored energy distribution obtained from FEM was then used in Monte Carlo method to simulate the microstructure evolution. The modeling results are compared with the theoretical and experimental results and an acceptable agreement is achieved.
Key words: aluminum sheet; Monte Carlo; FEM; recrystallization; stored energy distribution
金属材料的塑性变形行为对其介观结构敏感，通过有限元方法可有效模拟其金属变形过程；再结晶是一个复杂的材料组织物理变化过程，利用 Monte Carlo方法(简称MC法)，通过计算机模拟，可以定量、直观、连续地重现再结晶组织变化，计算金属组织演化和预报材料性能[1]。
根据再结晶理论，再结晶形核过程与变形材料组织内部的储存能(H)有关，且形核后新生成晶粒生长的驱动力也主要是储存能。实际上，绝大多数金属的冷变形是非均匀变形，其应变(ε)和储存能(H)的分布呈非均匀和各向异性。
然而，在传统的冷轧材料再结晶过程组织变化的模拟研究中，初始条件为一个假定的均匀储能场[2-3]，显然与实际出入较大。文献[4]提出的RSRP模型考虑了晶界及晶粒大小的影响，可以近似地反映介观冷轧储能场的非均匀性，但仍与轧件宏观储能场的均匀分布和理想的弹塑性材料假设有关。目前， 国内外已经有很多学者探讨将晶体塑性有限元方法(Crystal plasticity finite element method, CPFEM)与MC法进行耦合[2-3, 5-6]，以解决冷变形材料介观储能场的非均匀分布问题，取得较好的结果；但是，由于CPFEM实现的是介观上塑性变形行为的模拟，因此两者的耦合局限于介观尺度上再结晶的理论研究，因而不能真实反映轧件宏观储能场的非均匀分布，从而导致不能全面了解轧件整体的再结晶组织演变过程。
为了全面反映轧件宏观储能场的非均匀分布及其对再结晶模拟的影响，本文作者尝试将有限元技术与再结晶组织MC模拟方法相耦合，以有限元方法所模拟得到的冷轧材料宏观储能场代替RSRP模型[4]中的平均储能场作为初始条件，仍采用RSRP模型的择优取向转换规则[7]，在接近实际变形的条件下，进行了冷轧工业纯铝板的再结晶模拟。
1 耦合模型的建立
有限元方法与MC法的耦合主要解决以下几个问题：1) 有限元模型中流变应力→变形储能的转换；2) 变形储能移植为MC模型中宏观平均储能；3) 介观初始储能场的建立。
在金属塑性变形期间，所吸收的能量大部分用于改变金属的形状和转变为热能而散失，只有大约10%的能量以晶格畸变、点缺陷、位错、层错等形式储存在金属中，使金属内能升高而处于热力学不稳定状 态[8]。在具有高层错能的铝等冷变形材料中，位错成为主要的晶格缺陷和储能携带形式。
本文以储能密度作为结点储能的量度。根据Wilkens、Kocks和Scattergood的理论，储能密度可表示为[9]
式中 c为常数，对于刃型位错取值1/4π，对于螺型位错取值[4π(1-ν)]-1，本文综合刃型位错与螺型位错的作用，取为0.1；G为剪切模量；b为柏氏矢量；ρ为位错密度；Re为位错的平均影响范围。
而位错运动时，阻力τ和位错密度的关系如下[10]：
所以，流变应力可表示为[11]
将式(3)代入式(1)，则储能密度和流变应力的关系为
利用ANSYS中的PATH命令，可以定义路径并将有限元的计算结果(如应力应变和坐标等)映射到所定义的路径上，从而可获得整个稳定轧制区的应力应变等分布规律，进而由式(4)求出经过稳定轧制后材料内部各处的宏观储能密度分布。
在MC法模拟中，通常模拟区域很小，该区域内对应的宏观储能密度可视为均匀，表示为宏观平均储能密度，具体可由该方形区域4个顶点的储能密度通过几何平均法求得[12]:
(5)
式中 Hm为所选组织模拟区域的宏观平均储能密度；Hk为所选组织模拟区域4个顶点处的宏观储能密度。
根据以上方法，将有限元模型中的流变应力转换为储存能密度，从而得到了宏观的平均储能密度分 布。然而，由于有限元软件ANSYS是在材料各向同性的假设下，设置边界条件及初始条件进行迭代运算的，而没有考虑介观的晶粒大小及其晶界对储能分布的影响。因此，本文将宏观的平均储存能移植到Monte Carlo 模型中时，采用了RSRP模型的两个假设[4]：1) 各变形晶粒内部的平均储能密度取决于晶粒大小，晶粒愈小则其内部平均储能密度愈大；2) 晶界阻碍可导致其附近的位错塞积，因而晶界及其附近结点的储能密度较晶内大。由此得到了介观的冷轧变形晶粒组织的初始储能密度分布，即：
因此，在本文耦合模型所研究的材料体系中，任意结点i的能量为：
本文采用的再结晶形核模型系基于亚晶生长机制和以下两种假设条件所建立的。
假设1：在每个结点及其邻近区域内存在一个亚晶，其初始直径Di0为[13-14]：
假设2：在储能密度Hi的驱动下，亚晶在二维空间中各向同性生长，其生长过程可描述为[15]
一旦亚晶直径超过临界直径
，则认为该结点处发生再结晶形核。任一结点处的亚晶临界形核尺寸由下式确定[16]：
2 模拟的初始和边界条件
2.1.1 轧制初始参数
本文对压下率为50%的工业纯铝板轧制过程利用有限元模拟软件ANSYS/ STRUCTURAL/LS-DYNA进行了模拟。为了减少运算量，考虑到坯料及其受力的对称性，取坯料体积的1/2和轧辊体积的1/2建立计算模型。模拟中，取静摩擦系数0.3，动摩擦系数0.15，坯料初始速度0.85 m/s，轧辊转速为10 rad/s。具体的轧辊及坯料参数列于表1。
表1 轧辊与坯料部分参数
Table 1 Parameters of rollers and blank
2.1.2 储存能分布
为了反映冷轧纯铝板内变形储能的宏观不均匀分布，如图1所示，在对应稳定轧制区的横截面上(轧件1/2厚度模型)，本文分别沿板高向和宽向选取了1、2、3和4、5、6共6条路径来研究储能分布情况。图1中箭头所示方向即为路径方向，与对应的储能分布图2中的h、b变化方向一致。
图1 所选取路径及模拟区域示意图
Fig.1 Sketch map of selected path and simulation region
从图2可以看出，在冷轧板中，表层的变形储能变化最剧烈，心部次之，而沿1/4厚度处的储能分布最为均匀；最大储能区出现在上下表面与侧面相交处的4个对称区域，最小储能区则出现在两侧面的水平中心线处(如图1中a和b)。另外，通过模拟研究可以看出，板厚加大时，应力应变的变化会更加不均匀。这里不加赘述。
2.1.3 再结晶模拟区域的位置
图2 冷轧铝板变形储能的分布规律
Fig.2 Distribution of stored energy in cold rolling Aluminum sheet: (a) Stored energy curves through normal direction; (b) Stored energy curves through transverse direction
本文取轧件中储能为极值的两个区域作为研究对象，如图1所示，A区为储能极大区，B区为储能极小区。
本文作者将等轴晶粒组织进行网格变换，利用抽点法建立一个近似冷轧组织，使晶粒沿一方向(轧向)伸长，沿另一方向(高向)压缩。这一组织基本反映了纤维晶粒状的冷轧组织特征，能够较好地模拟再结晶初始组织形貌。
模拟采用200×200点阵，模拟视区为1 mm×1 mm，取向数Q=180，柏氏矢量b=1.65×10-10 m，剪切模量G=26.3 GPa，位错的影响范围Re=1.7×10-7 m，再结晶温度T=573 K，硬化指数n=0.204，晶界能量度J=0.324 J/m2。
3 再结晶模拟结果
初次再结晶的驱动力为塑性变形时所积累的变形储能。本文以再结晶分数x=99%作为再结晶完成的标志，如图3所示。
图3 储存能对再结晶拓扑组织的影响
Fig.3 Influence of stored energy on recrystallization microstructures (x=99%): (a) Microstructure in region A (MCs= 85); (b) Microstructure in region B (MCs=99)
1) 两区域所模拟的再结晶组织晶粒特征基本相同，均呈等轴状，且三晶界交角大部分接近于120?，基本符合再结晶组织的晶粒特征。2) 储存能越大(A区)，再结晶完成时对应视区内晶粒数越多，平均晶粒尺寸越小，再结晶过程越短；相反，储存能越小(B区)，再结晶完成时对应视区内晶粒数越少，平均晶粒尺寸越大，再结晶过程越长。
3.2.1 再结晶分数
如图4所示，两模拟区域的再结晶分数与退火时间均呈近似“S”形曲线，对数分析曲线中 ln{ln[1/(1-x)]}与lnt基本呈线性关系，与JAMK理论基本一致[15]，且更接近理论值2；储存能越大，t0.5和t0.99越小，再结晶速度越快，再结晶过程越短；反之相反。
图4 储存能对再结晶分数曲线及其对数分析曲线的影响
Fig.4 Influence of stored energy on recrystallization fraction and its logarithm analysis curves: (a) Curves of x—t and ln{ln[1/(1-x)]} —lnt in region A; (b) Curves of x—t and ln{ln[1/(1-x)]} —lnt in region B
3.2.2 形核率N
形核率N为单位时间内单位未再结晶面积内的形核数目[16]，即：
如图5所示，两模拟区域的形核率N随时间变化趋势基本一致，即初期N值较大，达到峰值后呈逐渐降低的趋势，至一定时刻不再形核。这与文献[17]实测的铝拉伸变形后再结晶形核率曲线基本一致。不同之处如下：1) 形核率N不同，储存能越大，N峰值越大，再结晶完成时晶粒数量越多；2) 形核时间不同，储存能越小，形核时间越长；3) 初期形核规律不同，高储能的A区单调增加，低储能的B区先减少后增加。由于B区储存能偏小，一部分满足形核条件的结点形核后，其余结点处的亚晶还没达到临界形核条件，因此需要一个孕育期，经过短暂的孕育期后，形核率会很快达到峰值。
图5 储存能对形核率的影响
Fig.5 Influence of stored energy on nucleation rate: (a) Nucleation rate curve in region A; (b) Nucleation rate curve in region B
3.2.3 平均晶粒尺寸
如图6所示，再结晶初期由于再结晶晶核的出 现，晶粒平均尺寸迅速降低；随后，由于再结晶形核率的起伏，晶粒平均尺寸也会随之波动；形核结束后，随着新生成无畸变晶粒的生长，平均晶粒尺寸缓慢增加。通过图6明显看出，再结晶完成后，低储能B区的平均晶粒尺寸明显大于高储能A区的平均晶粒尺寸。显然，冷轧纯铝板中变形储能的高低对应于该处的变形大小，变形储能对平均晶粒尺寸的影响规律与文献[16] 所指出的冷变形量增加造成再结晶晶粒尺寸下降的规律是一致的。
图6 储存能对平均晶粒尺寸的影响
Fig.6 Influence of stored energy on mean grain size: (a) Mean grain size curve in region A; (b) Mean grain size curve in region B
综上分析不难看出，再结晶动力学模拟统计结果与其拓扑组织直观模拟结果相同，且这些结果与现 有再结晶理论和其规律及相关文献实验结果基本一致，不仅表明再结晶拓扑组织模拟结果的必然性，而且验证了本文所采用的有限元-再结晶耦合模型是正确的。
4 结论
通过有限元模拟的流变应力场-储能场转换，可展现变形材料内部任意处的宏观平均储能分布，从而建立相应的介观初始储能分布的模拟模型，成功地实现有限元方法与Monte Carlo 方法的耦合，改进了传统的Monte Carlo再结晶模拟方法，使冷轧变形过程模拟与再结晶过程模拟相结合，实现了非均匀能量场基础上的冷轧纯铝板再结晶过程的模拟，模拟结果与理论和实验研究结果相吻合。
REFERENCES
[1] Srolovitz D J, Grest G S, Anderson M P, Rollett A B. Computer simulation of recrystallization—Ⅱ. Heterogeneous nucleation and growth [J]. Acta Metall, 1988, 36(8): 2115-2128.
[2] Rollett A D, Raabe D. A hybrid model for mesoscopic simulation of recrystallization [J]. Computational Materials Science, 2001, 21(1): 69-78.
[3] Radhakrishnan B, Sarma G B, Zacharia T. Modeling the kinetics and microstructural evolution during static recrystallization— Monte Carlo simulation of recrystallization [J]. Acta Mater, 1998, 46(12): 4415-4433.
[4] 张继祥. 基于Monte Carlo方法的材料退火过程模拟模型及计算机仿真关键技术研究[D]. 济南: 山东大学, 2006.
ZHANG Ji-xiang. Study on the Simulation Models of Annealing Process and Related Computer Emulation Key Technologies by Monte Carlo method [ D]. Ji’nan: Shandong University, 2006.
[5] Raabe D, Becker R C. Coupling of a crystal plasticity finite-element model with a probabilistic cellular automaton for simulating primary static recrystallization in aluminum [J]. Modelling and Simulation in Materials Science and Engineering, 2000, 8(4): 445-462.
[6] Lan Y J, Li D Z, Li Y Y. Mesoscale simulation of ferrite transformation from deformed austenite during continuous cooling in a C-Mn steel using a cellular automaton method [J]. Computational Materials Science, 2005, 32: 147-155.
[7] 张继祥，关小军，孙 胜. 一种改进的晶粒长大Monte Carlo模拟方法[J]. 金属学报，2004, 40(5): 457-461.
ZHANG Ji-xiang, GUAN Xiao-jun, SUN Sheng. A modified Monte Carlo method in grain growth simulation[J]. Acta Metallurgica Sinica, 2004, 40(5): 457-461.
[8] Gorelik S S. 金属和合金的再结晶[M]. 仝健民, 译. 北京: 机械工业出版社, 1985.
Gorelik S S. Recrystalization in Metals and Alloys [M]. TONG Jian-min, transl. Beijing: China Machine Press, 1985.
[9] Mohamed G, Bacroix B. Role of stored energy in static recrystallization of cold rolled copper single and multicrystals [J]. Acta Mater, 2000, 48(13): 3295-3302.
[10] 杨觉先. 金属塑性变形物理基础[M]. 北京：冶金工业出版社，1988.
YANG Jue-xian. Physics Elements of Metal Plastic Deformation[M]. Beijing: Metallurgical Industry Press, 1988.
[11] Hansen N, Huang X. Microstructure and flow stress of polycrystals and single crystals [J]. Acta Mater, 1998, 46(5): 1827-1836.
[12] 彭颖红. 金属塑性成形仿真技术[M]. 上海: 上海交通大学出版社, 1999.
PENG Ying-hong. Simulation Technology of Metal Plastic Deformation[M]. Shanghai: Shanghai Jiao Tong University Press, 1999.
[13] Godfrey A, Hughes D A. Scaling of the spacing of deformation induced dislocation boundaries [J]. Acta Materialia, 2000, 48(8): 1897-1905.
[14] Ye D, Matsuoka S, Suzuki N, Maeda Y. Further investigation of Neuber’s rule and the equivalent strain energy density (ESED) method [J]. International Journal of Fatigue, 2004, 26: 447-455.
[15] Huang Y, Humphreys F J. Subgrain growth and low angle boundary mobility in aluminium crystals of orientation {110} {001} [J]. Acta Materialia, 2000, 48(8): 2017-2030.
[16] 毛卫民, 赵新兵. 金属的再结晶与晶粒长大[M]. 北京：冶金工业出版社，1994.
MAO Wei-ming, ZHAO Xin-bing. Recrystallization and Grain Growth in Metal[M]. Beijing: Metallurgical Industry Press, 1994.
基金项目：山东省自然科学基金资助项目(Y2003F03)
收稿日期：2006-06－26；修订日期：2006-10-28
通讯作者：关小军, 教授，博士; 电话: 0531-88655799; E-mail: guanxj2003@126.com
