|
摘要：阐述了环件轧制三维有限元模拟中运用质量缩放方法的意义及其理论依据，并对显式动力学有限元模拟中各种质量缩放方法进行了分类，给出了环件轧制模拟中质量缩放方法的选取原则。最后基于ABAQUS/EXPLICIT操作平台，建立了铅环件的环件轧制三维有限元模型，通过具体模拟及分析比较计算结果对选取原则进行了验证，从而总结出了适合于环件轧制模拟的有效的质量缩放方法。
关键词：环件轧制 有限元模拟 质量缩放
中图分类号：TG
1 前言
大变形动态非线性问题（nonlinear problem of dynamic large deformation）[1]的计算和分析可以说代表目前有限元分析的最高水平，也是目前还在处于迅速发展的领域,其包括的方面主要有：板料成型工艺的数值模拟（汽车覆盖件），汽车安全性碰撞的全过程计算，材料的体积成形工艺（挤压、锻造、轧制）的全过程数值模拟等，环件轧制的三维有限元模拟就属于其中的一种材料体积成形工艺的全过程数值模拟。对于此类问题中存在的几何非线性和材料非线性问题，用于计算求解的方法主要有两种：静力隐式算法和动力显式算法。静力隐式算法的特征是迭代计算，在每一时间步内需反复迭代，迭代收敛性会受许多因素的影响，在计算中须调整迭代以满足收敛，而且其计算时间随模型单元数量呈指数增长，所以其计算时间长。而动力显式算法是递推计算，无需迭代，就时间步长直接进行递推计算，计算时间随单元数量呈线性增长，计算时间短，而且可以通过质量缩放来缩短计算时间，所以对于大规模计算和高度非线性问题[2]，动力显式算法总体效率要比静力隐式算法高，优势明显。近年来，随着在板料成型及碰撞冲击领域得到日益广泛的应用，对于材料的体积成形工艺领域，人们也都开始采用动力显式有限元法对其进行过程模拟研究[3]，而且取得了较好的成效。
动力显式算法在有限元分析中通常被用来解决两类典型的问题：结构动力学响应分析和准静态分析（其中包括复杂的非线性效应和接触条件），而环件轧制的模拟则是其中的一种准静态分析。由于动力显式算法采用的中心差分法是通过对时间求导进行计算，引入了质量矩阵[4]，所以在准静态分析中用于平衡方程式中的离散质量矩阵对计算效率和精度都起着至关重要的作用。从而在这类问题的模拟中，如何正确适当的运用质量缩放方法，使其能够起到在保持计算精度和稳定结果的同时提高计算效率的作用，也成为了一个至关重要的问题。目前，针对这个问题，人们也开始对其进行了相应的分析研究，东北大学在利用显式动力学有限元法模拟平板轧制过程中，通过计算实例研究了采用提高模拟轧制速度和质量缩放技术来缩短计算时间等方法的有效性，在忽略惯性惯性力影响的前提下得出了与实验结果相符合的计算结果[5]。燕山大学在用显式动力学有限元法分析板带轧制压力分布时提出了为了加快计算速度而采用的虚拟速度和虚拟质量带来的虚拟惯性效应会影响计算的结果[6]。而在环件轧制成形模拟方面，并没有人通过具体的实例分析来验证质量缩放技术在环轧模拟中运用的有效性和合理性。本文运用Abaqus 软件建立环件轧制的三维有限元模型，并采用动力显式算法对铅环件的轧制过程模拟进行了实例计算，验证了质量缩放技术在环轧模拟中运用的有效性，并通过对计算结果的分析比较，给出了适合环件轧制模拟的有效的质量缩放方法。
2 质量缩放在环轧模拟中的意义及理论依据
在Abaqus中利用Explict[7]来进行准静态模拟分析的求解，为了在较快的时间内获得可以接受的静态解，以提高求解的效率，通常通过两种方法来达到目的：一种是人为的提高成型模拟过程的速度；另一种则是通过质量缩放来提高模拟速度。对于环件轧制，模型的材料性能随应变率的改变而变化[8]，人为的提高进给速度就会人为的改变整个变形的过程，使得与实际不相符，而且在实际的环轧生产工艺中，进给速度最大值受环件咬入条件的限制[9]，人为的提高进给速度，如果超过了最大值，即使模拟能够得到稳定解，在实际轧制中由于环件不能满足咬入条件，轧制过程也是不可能正常进行的。所以人为的提高进给速度会使得模拟与实际脱离，失去了模拟用于验证和指导实践的作用，是不可取的。所以对于环轧的模拟，只能通过质量缩放来提高模拟速度。在Abaqus中，质量缩放是通过质量缩放因子来完成的，质量缩放因子取的太小达不到提高计算速度的效果，取的太大则会使得随着质量的增大而随之产生的虚拟惯性力的增大影响到计算的结果，甚至使得结果不能收敛。所以质量缩放因子的正确选取也就成为了环轧模拟中，质量缩放方法运用的关键。
动力显式算法的中心差分法是条件收敛的，其稳定的收敛条件是： （图片） 其中△t为时间步长，△tcr为临界时间步长，Tn是有限元系统的最小固有振动周期。在利用中心差分法求解具体问题时，时间步长必须小于由该问题求解方程性质所决定的临界时间步长△tcr，否则算法将是不收敛的。
临界时间步长通常可以近似认为是应力波通过任一单元网格所消耗时间的最小值，即：（图片） 式中Lmin为模型中单元的最小长度，与网格划分的疏密相关，cd为应力波在单元中传播的速度，对于金属材料可以表示为：（图片） 其中ρ为材料密度，λ和μ为拉梅常数，可以通过杨氏模量E和泊松比ν来表示： （图片） 代入上式则可以表示为成 （图片） 在Abaqus中用Explict进行求解计算时，所消耗的时间长短与增量数n有直接联系，n可以表示为n=T/△tcr ，其中T为整个模拟过程完成所需要的时间周期（对于环件轧制模拟，T即为轧制时间）。（图片） 则（图片） ，因为λ和μ与材料的物理性能相关，是常数，由此可以看出n值大小与时间周期T成正比，与单元长度L及ρ值成反比。n值越小，Explict在计算时处理这些增量所消耗的时间也就越短，计算效率也就越高。所以为了减小n值，可以通过减小时间周期和增大单元长度和密度来完成。在环轧模拟中，减小时间周期即为减小轧制时间，而轧制时间与进给速度成反比，减小轧制时间即意味着增大进给速度，这在前面已经说过是不可取的，而单元长度的大小代表着网格的疏密，增大单元长度即为疏划网格，网格划的太疏则会降低计算精度，达不到进行精确分析的效果，一般来说也是不可取的。而对于增大密度，根据塑性变形的体积不变定理，增大密度也就是增大质量，质量增大α倍，根据上面公式临界时间步长值△tcr就会增大（图片）倍，相应的n值就会减小（图片）倍。所以在Abaqus中，通过质量放大来达到等效于增大密度从而减小n值的作用，即能保证模拟的稳定性和精确度，又能有效的提高计算效率，是最可取的一种方法。
3 质量缩放方法的分类及选取原则
前面已经提过，动力显式算法主要用于解决结构动力学响应分析和准静态分析两种类型的问题，虽然都是用同一种算法，但这两类问题仍然有很大差别，Abaqus/Explict中用于动力显式算法的质量缩放方法有很多种，但这些方法对于两种类型问题的分析并不是通用的，大多数适合于准静态分析的质量缩放方法并不适合于结构动力学响应分析。由于环件轧制模拟属于一种准静态分析。所以这里只介绍适合于准静态分析的质量缩放方法。
Abaqus/Explict所提供的用于准静态分析的质量缩放方法主要有两种类型：
1.通过用户给定的质量缩放因子进行缩放。这是最简单直接的一种质量缩放方法，用户输入质量想要放大的倍数即缩放因子值，则Abaqus/Explict会将临界时间步长值放大倍，而且在计算过程的每一时间步中，都会将此时间步的DT值（Abaqus/Explict中以DT表示时间步长，以下均用DT表示）放大（图片）倍。
2.通过用户给定的DT值进行缩放。这种方法就是用户直接输入所需求的DT值，而不需要通过给出质量缩放因子进行缩放，但前提条件是所给出的DT值不能小于Abaqus/Explict根据用户所给出的材料性能参数计算出的原始DT值，否则就起不到质量放大的作用。这里需要说明的是Abaqus/Explict在进行求解计算时会对用户建模所生成的求解文件进行检查，并根据文件中给出的材料性能参数ρ，λ和μ从而估算给出一个原始DT值，也就是未进行质量缩放前的DT值，这个值其实也就是通常所说的临界步长值，但是Abaqus/Explict为了保证计算的稳定和收敛，所给出的原始DT值一般都要比实际临界步长值要小，对于三维模型，一般为实际值的（图片）/3～ 1倍。然后Abaqus/Explict再根据用户给出的质量缩放方法对DT值进行放大，如果给出的是缩放因子a，则Abaqus/Explict将DT值增大（图片）倍，如果用户直接给出了DT值，Abaqus/Explict则直接将其替换成用户给出的DT值，两个DT值的比值也就是用户运用这种缩放方法进行放大后，计算时间减少的倍数。所以用户在运用这种方法时必须保证给定的DT值不能小于原始DT值。对于这种方法，根据 Abaqus/Explict对DT值进行不同控制的操作，还可以将其分为三种方法。
a 直接根据给出的DT值进行缩放，对DT值不进行控制。这种方法其实和通过给出缩放因子进行缩放是同一原理。因为给出缩放因子，在表现方法上也是通过给出的因子值，将DT值进行相应倍数的放大。如果用户给出的DT值是原始DT值的4倍，效果就相当于给出的缩放因子值为16，则Abaqus/Explic会将每一时间步的DT值放大4倍。
b 用户给出DT值，在计算的每一时间步中，Abaqus/Explict自动对每一步的DT值进行控制，使得每一步的DT值在进行缩放后不超过用户最初给定的DT值。因为在用Abaqus建立模型时，网格划分采用的是自适应网格划分，那么在模拟过程中，随着模型形状变形，Abaqus会自动对模型进行自适应网格重划[10][11]。由于模型各个区域的变形程度不同，所以重划后的网格与最初网格相比，其各个单元的单元长度也发生了变化，有的单元的长度会减小，有的则会增大，随着模型变形的加剧，进行多次网格重划后的模型，其各个单元的单元长度变化也会越来越大，有的单元的长度会逐渐减小，有的则会逐渐增大，Abaqus/Explict在确定每一时间步的DT值时，为了保证所有单元的计算收敛，是以这一时刻的应力波穿过长度最大的一个单元所消耗的时间为这一时间步的DT值进行计算的。这样，如果这个DT值不超过临界时间步长则就能保证这个时间步计算的收敛性。所以在具体确定某一时间步的DT值时，随着有些单元的长度的逐渐增大，应力波穿过最大一个单元所消耗的时间也就越大，那么这一时间步的DT值就会增大，可能会超过原始DT值，甚至有可能超过了临界时间步长值，那么这一时间步的计算就不收敛了。为了保证计算的精确度，Abaqus/Explict中的自适应网格划分在每一时间步都会对模型网格进行重划，也就是说整个计算过程的每一时间步的DT值都是在变化的。所以这种缩放方法通过对每一时间步的DT值进行控制，使其不超过原始DT值，以避免出现DT值超出临界时间步长值，使得计算结果不收敛的现象，对于大变形的准静态分析是非常合适的。
c 用户给出DT值，在计算的每一时间步中，Abaqus/Explict自动对每一步的DT值进行控制，使得每一步的DT值在进行缩放后都等于用户最初给定的DT值。这种方法通过缩放DT值使每一时间步的DT值都恒等于最初的DT值，这样也就避免出现上面所说的随着网格重划而出现的逐渐增大的DT值超过临界时间步长值使得计算不收敛的现象。所以这种方法对于保证准静态分析过程的收敛性也是很有效的。
以上几种缩放方法都适合于准静态分析。但是，无论对于通过质量缩放因子缩放，还是给出DT值缩放，如果缩放系数太大，则随之产生的虚拟惯性力会影响到计算的精度和收敛性，缩放系数太小则达不到提高计算效率的目的。在准静态分析中[12]，为了保证模拟过程的稳定性，一般要求变形的动能值较小，且变化平稳，势能值则较大，动能与势能的比值较小，一般不能超过0.1，这样使得能量大部分用于模型的变形。所以一般在准静态分析中，质量缩放方法的基本选取原则就是：无论采用何种方法进行缩放，必须保证在整个过程中的动能值要较小，动能与势能的比值不能超过0.1。在满足这个前提条件的情况下进行缩放，一般都是可行的。
4 质量缩放方法的实例比较
本文以Abaqus软件为操作平台，以工业纯铅为材料，建立铅环件轧制的三维有限元模型。并用以动力显式算法为基础的求解器Abaqus/Explict进行求解计算。有限元模型如下：（图片）
图1 环件轧制三维有限元模型图 工业纯铅密度ρ=11340kg/m3 ,弹性摸量E=17GPa ,泊松比ν=0.3，最初网格划分的单元长度为0.01m，根据上文所给出△tcr公式，将上述值代入求得临界时间步长理论值为△tcr=4.955×10-6s，而Abaqus/Explict在检查求解文件时所给出的原始DT值为2.511×10-6s
由此可以看出为了保证计算过程的稳定性，Abaqus/Explict所给出原始DT值要比理论值小。
首先采用给出缩放因子进行缩放的方法，来比较不同数量级的缩放因子，对模拟结果的影响，本文分别取缩放因子为4、100、1000对环件轧制过程进行了模拟，环件最终的成形状况如图2所示:（图片）
图2 不同数量级的缩放因子时的环件最终的成形状况 图2中从左至右依次分别为质量缩放因子为4、100、1000时环件最终的成形网格图。通过对比可以看出，最左边的图效果最好，为圆环形；中间的图基本为圆环状，但左中部分呈明显的椭圆形，效果比要差；最右边的图效果最差，环件基本没有变形左中部分就被压扁了。由此，单从形状上看，采用1000的缩放因子是不合理的。缩放因子分别为4和100时，环件端面宽展变形，环件变形的动能历史，及动、势能比值比较如下：（图片）
图3 环件端面的宽展变形图4 环件变形的动能历史 图5 动、势能比值 图3横轴是环件的径向尺寸，从环件的内环面开始为0，纵轴是环件端面的宽展距离。从图中可以看出，缩放因子为4时，环件沿径向的宽展变形较平缓，基本符合实际轧制中两端宽展大，中间宽展小的情况，成形后的端面凹陷不会很明显。而缩放因子为100时，环件沿径向由内向外呈明显的上升趋势，外端面的宽展要比内端面大很多，那么成形后就会形成很明显的端面凹陷[13]。再比较环件变形的动能历史，从图4可以看出，缩放因子为4时，环件在整个变形过程的动能值较小，且变化平缓，符合上文所述的准静态分析中，要求有较小且变化平稳的动能值的原则，这样环件变形过程就会较稳定。缩放因子为100时，在同一时刻环件变形的动能值几乎都要比前者大很多，而且整个变形过程中动能变化很不平稳，升降幅度很大，不符和上述原则，这样环件变形过程就会不稳定，产生振动或摇摆的现象。最后比较环件的动、势能比值，从图四中可以看出，缩放因子为4时，在变形的初始阶段，环件的动、势能比值超过了0.1，最大值为0.3，且不稳定，有波动现象，到大约0.4秒时比值下降为小于0.1，随后趋于平稳。而缩放因子为100时，动、势能比值不仅明显超过了0.1，而且比值较大，最大值达到了1.3，比同一时刻的缩放因子为4时的比值要大很多。在接近1秒时比值才小于0.1趋于稳定，不稳定时间持续较长，几乎为整个过程的1/3，而在稳定阶段，比值也比缩放因子为4时高。由此，从以上3幅图的比较可以看出：无论从工艺需求的宽展变形和端面凹陷方面，还是从Abaqus/Explict本身进行准静态分析时对动能输出及动、势能比值的要求方面来看，缩放因子取4比取100要更合理。
以上通过对比分析确定了缩放因子取4时符合本文所建铅环件模型的模拟需求。下面比较直接给出缩放因子进行缩放即上文所述方法1和通过给出DT值进行缩放，并对DT值进行不同控制的缩放方法即上文所述的方法2中的方法b和方法c（前面已经讲过方法2中的方法a和方法一其实是一样的，所以在下面的比较中以方法a代替方法1）这三种方法对模拟成形的影响。模拟过程中仍以缩放因子4为前提，转化为DT值即为放大2倍。采用3种不同方法时，环件端面宽展变形，环件变形的动能历史，及动、势能比值比较如下：（图片）
图6 环件端面宽展变形 图7 环件变形的动能历史 从图6可以看出，三种方法中环件端面的宽展变形并无明显差别，由此可以看出这三中缩放方法对宽展变形的影响不大。从图7可以看出，虽然三种方法的动能变化都较平缓，但同一时刻方法b的动能值比方法a和c明显要小，方法a和c的动能值变化则基本一致。图8为三种方法的环件变形动、势能比值比较，从图中可以看出，方法a和c刚开始阶段比值都超过了0.1，但c相对a要好，最大比值为0.17，而a最大值为0.31，c的比值波动的幅度要比a小，比值降为小于0.1时所经历的时间也要比a短。方法b的效果则最好，环件变形的整个过程中动、势能比值都要小于0.1，最大值为0.095，而在随后的稳定阶段，同一时刻方法b的动、势能比值也要比方法a和c的小。图8为计算过程中每个时间步的DT值输出。由于Abaqus/Explict给出原始DT值为2.511×10-6s，而最初都是将DT值放大2倍即为5.022×10-6s，所以从图中可以看出三种方法刚开始的DT值都为5.022×10-6s。方法a不对DT值进行任何控制，从图中可以看出有些时间步的DT值明显大于最初值，虽然在本模型的分析中，这些较大的DT值并没有超过实际的临界时间步长而导致计算结果不收敛，但在其它模型的分析中，如果Abaqus/Explict所给出的原始DT值与实际临界时间步长值很接近时，那么运用方法a进行缩放时，较大的DT值就可能会超过实际临界时间步长值而使得计算结果不收敛。方法b是控制DT值，使每步的DT值都不超过最初值，从图中也可以看出，整个过程的DT值处在一个下降的趋势中，且都不超过5.022×10-6s。方法c是控制DT值，使每步的DT值都等于最初值，从图中明显可以看出，整个DT值变化为一条直线都等于5.022×10-6s。另外，对于使用三种缩放方法进行模拟时，记录CPU计算消耗时间分别为1小时39分，1小时40分和1小时42分。可见三种缩放方法计算效率几乎一样。（图片）
图8 环件变形的动、势能比值 图9 环件变形的DT值比较 从以上四幅图的分析比较中可以看出方法b既符合准静态分析对动能输出及动、势能比值的要求，又能保证有效的计算效率，所以是三种方法中最适合本模型分析的质量缩放方法。
总结以上的分析比较可以得出：准静态分析中质量缩放方法的选取原则对环件轧制的模拟是适用的。对于环轧模拟中缩放方法的运用，首先应该从环件成形的形状、宽展变形、动能变化及动、势能输出等方面进行比较，确定一个较合适的缩放因子，然后再运用方法b，通过对每一时间步的时间步长值进行控制，保证其不超过临界时间步长值。这样，即提高了计算效率又保证了计算过程的稳定性和收敛性。
5 结论
本文采用不同的质量缩放方法及放大系数对铅环件三维有限元模型进行了实例分析，通过分析比较不同缩放系数时的模拟结果，验证了在准静态分析中，为了保证模拟过程的稳定性，质量缩放方法的选取必须遵循：通过缩放后的整个模拟过程其变形的动能值应较小，且变化平稳，势能值则较大，动能与势能的比值应较小，一般不能超过0.1，这一选取原则对环件轧制模拟是适用的；通过分析比较采用不同缩放方法时的模拟结果，总结了一套适用于环件轧制模拟的最有效的质量缩放方法。
参考文献
[1] 曾攀. 有限元分析及应用. 北京：清华大学出版社，2004
[2] Belyschko T et al. Nonlinear Finite Elements for Continua and Stuctures.
[3] 刘立忠. 隐式静力和显式动力有限元在轧制过程模拟中的应用. 塑性工程学报，2001，8（4）：81～83
[4] Guyan R J. Reduction in stiffness and mass matrices. AIAA journal, 1965, 3(2):380
[5] 刘立忠, 刘相华. 利用显式动力学有限元法模拟平板轧制过程. 塑性工程学报， 2001，8（1）：51～54
[6] 张国名, 肖宏. 用显式动力学有限元法分析板带轧制压力分布. 燕山大学学报， 2004，2（24）：146～149
[7] 庄茁等译. ABAQUS/Explicit有限元软件入门指南. 北京：清华大学出版社,1998
[8] 罗洲. 环件轧制过程的显式有限元模拟分析. 塑性工程学报，2004，1（11）：68～70
[9] 华林，黄兴高，朱春东 . 环件轧制理论和技术. 北京：机械工业出版社，2001
[10] 江雄心，万平荣，游步东 . 三维有限元模拟中的网格重划. 金属成形工艺， 2002，2（20）：21～24
[11] Lee N K, Yoon J H, Yang D Y . Finite element analysis of large deformation by automatic renoding as a weak remeshing technique . Int J Mesh Sci, 1992, 34(4):255～273
[12] Abaqus Analysis User’s Manual. ABAQUS Corp.,2003
[13] 华林. 环件轧制成形原理和技术设计方法:[博士学位论文]. 西安:西安交通大学,2000
1/23/2006
|
