|
1．引言
上世纪八十年代，一个由汽车制造工程师和模具设计师组成的团队开始了INDEED软件的研发，这个软件名来源于Innovative Deep Drawing（创新的深拉延）。顾名思义，这是一个专用于金属深拉延的有限元数值模拟分析软件。在当时，显式有限元法非常流行，特别是涉及到高度的非线性材料响应，有限形变，以及在由于复杂的接触问题导致边界条件发生突变的情况下，显式有限元法有其独特的优势。鉴于板料成形过程是一个准静态过程，工程师们还是决定在INDEED 中采用隐式算法；一个明显的优点是，可方便而精确地进行回弹计算。如今，尽管显式算法在解决准静态问题时有诸多的缺点和不足，但其代表性软件（比如LS−DYNA）在板料成形数值模拟市场上似乎已经成为主导。这种统治地位是基于一种广为流传的思想，那就是隐式算法无论是在计算的可靠性还是效率都比不上显式算法。然而，尽管在很多情况下，显式算法可以得出令人满意的结果，但是在板料成形的模拟应用中，使用隐式算法还是有其突出的优越性。特别是在计算冲压力和板料回弹的情况下，隐式算法得出的结果更加可信。由于像超高强度钢、双相钢和特殊铝合金这样的新材料不断被采用，其仿真计算变得越来越重要。
在INDEED 的研发中，工程师们把重点放在了针对更加复杂工况的分析结果质量上，例如精确的回弹计算和冲压力计算。因此，INDEED采用了隐式的求解方法，使用了一种独创的特殊壳单元，这种单元在厚度方向有一个附加的自由度，使用了领先的材料模型，将材料的变形梯度采用乘法分解为弹性和塑性两个部分，这些技术奠定了Indeed——高端成形模拟软件的基石。
过去，使用隐式算法求解问题时常常由于有限元模型非常大，必须要使用高性能的计算机。现在，由于采用了新型的基于域分解的有限元技术，隐式算法的计算速度已经大幅提高。对于那些需要计算大型成形模拟而又不愿购买昂贵的矢量机和大型机的小公司来讲，这是个好消息。最近几年中，INDEED的研发团队已经推出了FETI-INDEED，这是INDEED的并行版本，它基于FETI(Finite Element Tearing and Interconnecting)域分解技术。FETI-INDEED遵循MPI标准(Message Passing Interface)，可以运行在异构硬件平台上，也可以运行在像LINUX集群这样的同构平台上，而LINUX集群在汽车工业中正变得越来越重要。
正如INDEED这个名字的寓意一样，INDEED致力于深拉延问题的数字化模拟。现在，INDEED已经涵盖了从弯管、轧制、翻边、接合到碰撞、液压、液压深拉延等几乎所有的金属成形模拟领域（如图1所示）。 （图片）
图1：INDEED的应用领域 除了上述的高精度的方法外，INDEED还集成了特别设计的膜单元，这样使用户可以快速对成形过程进行大致的评估。
本文的后续部分，将主要描述INDEED软件的主要特点。重点在INDEED独创的壳单元和材料模型上，这两个独创的技术正是INDEED的特别之处。另外阐述了FETI的实现方法，进而展示了FETI-INDEED进行模拟的实例，以便进一步理解这种方法的可扩展性。
2. INDEED的特点
2.1 积分方案
在使用LS−DYNA或者PAM－STAMP等显式有限元方法做深拉延模拟时，为了节省计算时间，通常都要放大成形过程的速度和（或者）板料网格的密度。第一种方法是为了减少模拟的时间，第二种方法为了增加所谓的临界时间步，从而减少时间步的总数。这些方法必然会导致分析出现不符合实际工艺过程的动力学效果。例如，冲压力的值就会被夸大；会发现在实际工艺中只发生弹性变形的情况，而在模拟结果中发生了塑性变形。因为有这样的情况，在分析显式深拉延模拟结果时，必须额外小心。掌握这种显式算法，自然也要求掌握大量的力学背景知识并具有丰富的计算经验。这些要求对于仅在特定问题上需要进行数字模拟的模具设计师来讲要求过高。为了避免这种为提高效率而改变一些参数以便消除与实际结果不符的动力学响应情况，INDEED采用了隐式算法，在每个增量步都进行平衡迭代。
在对成形过程进行显式的模拟时，通常要有数十万个时间步，而使用INDEED隐式模拟时，大约需要100到200个增量步，当然这还要取决于求解问题的非线性程度。显式积分法避免了单元刚度阵的组装和与之相应的系统方程的求解。为了提高求解效率，必须尽可能把每一个时间步所要进行的运算次数压缩到最少。在显式算法中，最耗费时间的操作有三个，第一是本构方程的求解，第二是单元算法，第三是接触搜索。为了提高显式有限元算法的效率，软件开发者倾向于将这些过程尽可能的简化。基于这个原因，显式算法中通常只使用最简单的本构方程和最原始的单元算法，也就是Belytschko−Tsay壳单元。同样的原因，显式算法通常在壳单元的厚度方向采用很少的积分点。
在隐式算法中，与求解巨大的线性方程组相比较模型变化对计算时间的影响可以说是微不足道的。选择相对复杂的材料模型或者选择一种在有限元单元上需要大量计算的有限元模型，都不会对隐式算法的效率产生太大的影响。甚至增加厚度方向的积分点对隐式成形模拟求解速度的影响也不是很大。
2.2 单元
INDEED 集成了实体单元、壳单元和膜单元，可以涵盖所有的成形模拟。使用者可以选择比较快速的膜单元算法，或者比较准确的壳单元。使用膜单元可以得出一个大致的结果，因为膜单元的弯曲特性仅仅是通过几何关系而不是物理关系得出的。膜的特性和弯曲特性就可以相互独立，分别进行计算，因而可以采用较快的迭代方程求解器。对成形过程的第一次评估，或者在对模具的几何形状或者冲压的一些参数进行反复改进的时候，一般可以使用膜单元。使用膜单元，在减薄、材料流动和起皱预估等方面都具有足够的准确性。不过，如果需要更加精确的计算结果，例如需要计算准确的冲压力或者回弹，或者是在做最后一次的成形过程校验的时候，我们还是推荐使用壳单元。这种单元是为INDEED软件特别开发的。它满足以下的要求：
·较好的弯曲和膜特性；
·可以处理材料和几何非线性；
·可以处理双面接触。
基于传统壳理论的壳单元不适合于描述双面接触，因为这个理论采用平面应力假说（σ33 = 0），常常忽略厚度变化（ε33 = 0）。在真实的金属成形过程中，在板料和模具的接触区域，法向的应力是不为零的，并且由于大的薄膜应变导致的减薄效应对接触区的影响很大，从而影响到坯料上的摩擦力分布，因而它们不应该被忽略【1】。
为了克服传统的壳单元的缺点，SCHOOP【2】教授提出了特殊的壳单元理论－－导向壳理论，在1989年这种单元在INDEED【3】中的第一版中得到了实现。在其后的15年中，为了减少锁死效应和改进收敛性，最初的INDEED壳单元已经做了许多的改进。当需要得到准确的冲压力和预知板料回弹特性时，INDEED的这种壳单元相对于其他板料成形模拟软件使用的壳单元，特别是那些基于传统壳理论的壳单元，其优越性就显而易见了。
在新版本的INDEED (V8.0)中，又增加了一种新的壳单元。这种单元是为了满足高精度回弹分析要求而特别设计的。它基于改进的导向壳理论【4】，可以作为先前壳单元的一个替代。和以前的壳单元一样，它是一个3节点单元，有21个自由度，每个节点7个自由度。这些自由度的物理意义如下：
（图片）代表在有限元节点n处下表面的位置矢量的三个分量
ux1n, ux2n, ux3n代表在有限元节点n处上表面的位置矢量的三个分量
αn 附加的自由度，主要是考虑到单元在弯曲变形的时候，材料的中面和板料几何中面的偏移（图片）
图2：具有双面接触和厚度方向附加自由度的INDEED壳单元示意图 由于有了第7自由度α，直线假说可以表达为下面的形式（图片） 其中，x为壳单元内任一点的位置矢量，r为中面的位置矢量，d = (ux−lx )/t是在纤维方向的指示器。Z是纤维方向的坐标。t表示初始的单元厚度。从上面的等式中容易得出，几何中面和材料中面的距离等于8倍的第7自由度α。
对于传统的有限元来说，在单元体内的位置矢量是通过形函数插值得到的。对实际几何体也就是单元变形后的几何形状的位置矢量进行偏微分，相对于初始的位置或者单元参考构形可以得出变形梯度F，由此可以计算出有关应变的几个度量值。INDEED的新壳单元是一种非协调的壳单元。只有相对于单元中面的恒定变形梯度F0是通过上述方法得出的。单元中面外侧的应力值是通过插值得到的，这意味着它们不一定和给定的单元位移场（假设应变理论）相匹配。
对于这个新的壳单元，有以下的假设：
a） 在Kirchhoff壳理论的基础上，单元的曲率是线性分布的。
b） 横向剪应变是通过平衡方程式Q = t • M，根据曲率的可变部分计算得到的。
c） 相对于变形后的单元节点形成的三角形单元平面，曲率和横向剪应变的方向并不与之垂直，标量积Fo• dn是弯曲和剪切变形的离散化的度量。
d） 纤维相对于单元平面的转动角不是通过自身的单元曲面曲率，而是由形函数估算出来的。在单元的边线上，处理方法类似于Timoshenko梁的处理方法，因而相邻的单元在它们的共同边界上有很大的协调性，即沿边长的C1连续性，和垂直于边长的C0连续。然而，由于沿单元的曲面的可积分性不是强制保证的，在单元内是非协调的。
新的壳单元有以下的改进：
a） 具有非常好的收敛特性，也就是说，壳体结构的网格可以很粗。
b） 不会出现任何的剪切闭锁死现象。
c） 没有任何的沙漏模式，不需要做任何的稳定算法，只要和至少与一个其他的单元相连。
在8.0版INDEED中，开始采用这种新的壳单元。我们期待这种壳单元能够在回弹特性分析方面有很大的帮助。这里要提及的是，尽管该厚壳单元只有3个节点，其定义却类似于体单元，也就是说用户必须要定义6个节点，3个在上面，3个在下面。单元外表面相对其中面的运动变形由程序自动计算。
表1给出了INDEED中的所有单元类型。（图片）
表1：INDEED不同单元类型的应用总结 2.3 材料模型
在影响回弹和冲压力计算准确性的诸多因素当中，材料模型起关键作用。由于安全和轻型化需求，近几年来新的钢种，例如超高强度钢和双相钢在车身的设计中使用的越来越普遍。由于它们有很高的屈服强度，这些材料会产生的很大的回弹，也需要很大的冲压力，这些对模具的设计制造带来了很多的问题。为了解决这个问题，就需要更精确的模拟分析，特别是对新的钢种还没有实践经验可供借鉴，而且实际上能够描述这些新钢种复杂变形行为的材料模型并不是很多。
为了精确描述高强度钢和超高强度钢在塑性状态下复杂的材料行为，INDEED容许不同硬化模型的组合，比如各向同性，随动强化和一般变形强化模型。这就意味着，屈服面可以在偏应力空间内伸展、移动和旋转。和大多数的商用有限元软件比较，特别是显式的软件，在INDEED中，弹塑性材料模型是基于将总的变形梯度F乘法分解为弹性和塑性两个部分: F = Fel Fpl。这意味着变形过程被分为连续的两个部分，塑性和弹性过程。质点首先从初始状态弹性变形到中间状态，从中间状态到最终状态假设都是塑性变形。在许多的商业软件中，变形率也被附加分解为弹性和塑性变形率。而对F的乘法分解可以导入弹性势能，对变形率的附加分解会导出亚弹性或增量弹性模型。在亚弹性或增量弹性模型中，塑性行为是和路径相关的。就是说，从相同的初始位置出发的两条不同的路径，产生相同的总体变形，却有两种不同的应力状况。
在INDEED中质点材料状态是由一系列的变量描述的， ζ= (C, N, Z1, κ1, κ2)，它们决定了单元的外在（可以观测到的）和内在（内存中）的状态。外在状态是由右Cauchy−Green应变张量定义为C=FTF。塑性变形的度量N =(Fpl)−1。对称偏差因子Z1描述了在偏应力空间中屈服面的中心位置。两个标量硬化参数κ1 和 κ2作为内部状态变量。 Z1的引入描述了随动强化，而等向和变形强化在屈服函数中用材料四面体B(κ2)表示。内部状态变量是由偏微分方程控制的。这个特殊的方程最初是由KRAVIETZ【5】在1988年为INDEED开发的。从那以后已经做了多次的改进，并且模拟了许多工况来验证这种材料模型。为了给出INDEED的弹塑性材料模型所要求的参数，还特别研发了为双向载荷试验的十字试样。
同时可以使用基于Hadamard理论的失效准则，借此可以预测到颈缩现象。
在INDEED中，使用Newton−Raphson迭代法来求解非线性材料方程组的。对于小的塑性变形，这种方法等效于半径回归法。
对于许多常用的成形用钢，INDEED已经预先定义了其材料参数，以材料库的形式和INDEED软件一道提供给客户。
2.4 域分解法
FETI【6】方法是由FARHAT 和 ROUX【7】研发的，可以认为是不同架构计算机并行有限元软件的实施策略，在这种方法中，全局的内存分布在不同的工作站中的。
FETI方法，就是域分解方法。域分解法是建立在将全部的网格分为不同的部分也就是域的基础上的，如图3所示。全局的有限元的求解是通过迭代完成的。在每次迭代的过程中，局部问题通过与临域相连的子域的边界条件得到求解。通过迭代过程得到边界方程的右边项，从而每个子域的局部解完全满足全局解的条件。由于在每个子域内都满足局部平衡，只要相邻的子域沿其边界能够满足两个条件，全局也可以达到平衡。第一个条件是位移的连续性，这样全局的位移是几何相容的。第二个条件是，相互作用力的平衡。（图片）
图3：域分解和分布计算 域分解法是基于将子域之间的相互作用力作为一个接口参数。给定这个参数一个值，并且它保证在所有的接口处的相互作用力平衡，这样，在每个子域和其交界处的边界条件则规定为外力。由于相互作用力是平衡的，只有作用力沿着交界处是连续的，局部问题的解就是全局问题的解。如果作用力不是连续的，则设定的相互作用力就是不正确的。如果在位移上有一个跳动，则相互作用力必须减小以消除这个跳跃。如果子结构发生相互穿透，则是设定的相互作用力太大了，也必须减小。因此，沿着子域之间的边界的位移跳跃实际上是边界问题的残留部分，由其引起的相互作用力用迭代过程来修正，即采用共轭梯度算法。
在这个过程中，对于每次迭代，每个子域的信息都会传递给它的相邻子域。而非相邻子域之间的信息则不能在一次迭代中交流，它们之间的信息只有在一定次数的迭代后才能传递到对方，这个数字相当于在网络划分中它们之间的拓扑距离。由于有这个特点，当子域数量增加的时候，许多域分解法都会出现收敛速度急剧退化的缺点。而在FETI法内有一种机制，使得收敛的速度和子域的数量没有关系，这样就可以划分许多的子域。由于规定相互作用力是边界的边界条件，每个浮动的子域（其中没有一个节点是固定的）被定义为未知刚体运动部件。因而，每个域分解法的迭代过程分为两步，第一是局部步，在这一步中，每个子域在消除刚体位移后的相互作用力作为边界条件，求解局部问题。第二是全局步，在这一步中，通过调节所有浮动子域的刚体位移，来消除因此导致的位移跳动。在全局步中，需要隐式求解所谓的粗化网格问题，其未知量就是局部的刚体位移。这个问题就类似于一个单元是子域的有限元求解问题，单元的自由度就是刚体的运动。
正是由于有这个粗化网格过程，域分解法在子域数量上是可以扩充的。由于在每次迭代求解过程中，大部分的计算是局部问题的独立求解，因而它很适合于分布存储式并行计算。而数据的传输也主要是边界位移的交流。因此，这个算法有很好的粒度，并行运算的效率很高。上边的两大特性使得并行计算的子域划分性能很好。
为了将有限元网格分解为不同的子域，选择了METIS【8】包。每一个子域和一个运行在工作站上的进程相连接。边界节点就是相邻子域共有的节点。
借助于拉格郎日乘子，在边界节点的可互换性得到了加强。子域的局部奇异性通过两个阶段来解决。首先消除每个局部问题的刚体运动，然后同时对所有子域直接求解，以重新获得所有的部分解。再通过正交条件，将这些子域模式的贡献与拉格郎日乘子相关联。使用共轭梯度算法，求解出局部刚体模式和拉格郎日乘子的耦合系统方程，从而完成问题的求解。
由于数据是在集群的若干节点上分布的，一台机器必须通过网络进行通讯，以控制迭代过程，或者将迭代法得出的局部解组成全局解。进程之间必须进行的通讯方式采用了精心设计的MPI标准。
每一个进程中都包含一个顺序执行的INDEED任务，其中的数据是相对应子域中略加改动后的数据。图4是在一个时间步长上的流程图，图的右边指出了通讯的数量。（图片）
图4: 一个时间步中的流程图 为了展示FETI−INDEED的可扩张性，参照VDI ’91测试标准，采用了一个相对较小的深拉延例题，其中有700个3节点单元使用了12250个自由度。图5是使用普通的顺序执行的INDEED和FETI−INDEED分别计算出的冲压力的对比。FETI−INDEED中将网格分为4个子域。（图片）
图5：冲压力的比较 顺序执行的INDEED和FETI−INDEED得出了相同的结果。在不同阶段所需的CPU时间如表2所示。容易看出，即使对于有较少单元的模型，FETI−INDEED的表现也很优秀。表2：INDEED计算小型问题时的性能表
（图片）与采用传统的（共享内存），非基于域分解技术的并行方式相比较，可以有以下的结论：
·传统的并行成形模拟应用程序只能利用8－12颗的CPU。
·并行计算有一个较小数量的额外开销。使用3颗CPU运行FETI−INDEED，和单颗CPU运行普通INDEED版本的速度相当。
·对于像汽车侧围这样的有大量自由度的大型模型，FETI−INDEED可以支持至少100个处理器，因此其可扩充性是非常好的。
·相对便宜的LINUX集群可以当作超级计算机使用。
当使用共享内存INDEED版本的时候，数量巨大的单元导致对内存和带宽的要求很高，即使在降低并行效率的情况下，也需要64位选址机型才能满足这个要求。基于RISC的系统可以满足这个要求，但是这些都是比较昂贵的服务器系统。而采用FETI−INDEED可以在很多的处理器上进行模拟，每个进程所需的内存数量远远低于32位选址机型的限制。所以，FETI−INDEED结合了隐式算法的准确性和显式算法的高效性。
3. 结论
INDEED是一个高端的有限元模拟软件，几乎适用于所有的金属成形过程。由于采用了隐式积分法，基于改进壳理论的先进壳单元技术和为新钢种和铝合金成形特别研发的材料模型，INDEED可以精确地计算回弹和冲压力。
INDEED是第一个采用FETI技术的商业化有限元软件。FETI−INDEED对于小型有限元网格也展示了很好的可划分性。它已经在100多个处理器的平台上做过测试计算，证明是一种很有效的方法。INDEED使得小型公司以较低的成本，使用隐式解法对大型的模拟问题进行精确的求解。
现在已经发布了基于分布计算的INDEED版本――FETI−INDEED。它可以使用廉价的LINUX集群和异构平台，例如CAD工作站和PC作为强大的计算平台。
除了高端方案，INDEED还提供了一种快速方案，这种方案采用高效的膜单元，并采用迭代方程求解器。使用这种方案可以使用一台工作站或者PC，不需要很长的CPU时间和很多的内存占用，就可对成形过程做一个大致的模拟，从而对模具进行改进。最后的对模具的几何形状及工艺参数的验证，则使用基于精确的实体壳单元和双面接触的INDEED或FETI−INDEED来完成。
参考文献
[1] El Rifai Kassem, K. et al. Mathematical Modelling of the Sheet Metal Forming Process with INDEED, (J.−L. Chenot, R.D. Wood, ed.), Numerical Methods in Industri
[2] al FormingProcesses, NUMIFORM 92, Balkema, 1992
[3] Schoop, H. Schalenelement mit 6 kinematischen Freiheitsgraden bei großen Verschie−bungen, ZAMM 67 (1987) 4, T 237−239
[4] Schoop, H. Entwicklung eines Viereckelementes fu¨ r den Tiefziehprozeß, Forschungsberichtim Auftrage der INPRO, Berlin 1989
[5] Schoop, H. Dokumentation zu SHELL1, Forschungsbericht im Auftrage der INPRO, Berlin，2002
[6] Kravietz, A. Eine Stoffgleichungsstruktur fu¨ r Tiefziehbleche. Projektbericht im Auftrage derINPRO, Berlin, 1988
[7] Kessler L, Weiher J, Roux F −X, Diemer, J. Forming Simulation of High− and Ultrahigh−Strength Steel Using INDEED with the FETI−Method on a Workstation Cluster, Proceedings of the 5th International Conference and Workshop on Numerical Simulation of 3D Sheet Forming Processes, 2002, Korea
[8] C Farhat and F −X. Roux, Implicit parallel processing in structural mechanics (J. Tinsley Oden, ed.), Computational Mechanics Advances, Vol. 2, North−Holland, 1994
[9]http://www−users.cs.umn.edu/~karypis/metis
3/8/2006
|
