有限元数值模拟中的网格重划技术

第六章 有限元数值模拟中的网格重划技术
在用有限元方法模拟形状复杂工件的大变形过程中，随着计算过程中变形量的增加，原始定义的计算网格会逐渐畸变。

若把已经畸变的网格作为求速度增量的参考状态，会导致不精确的解，甚至无法继续进行计算。

为了使计算顺利进行，最终得到满意的解，必须严格控制单元的变形程度和单元节点的疏密布置，防止出现计算特性不好的单元。

因此，在每一个加载结束后、下一个加载开始之前，必须进行网格畸变的判断，以便于在网格变形过程中及时对计算特性不好的网格进行重划。

网格重划技术是成功模拟大变形时必须解决的关键技术，其核心内容是新旧网格之间形状和信息的准确传递，网格重划技术一直是大变形有限元计算的研究的热点之一]84~81[。

在研究网格重划技术之前，先介绍一下单元质量的评定和网格自适应技术，它们是网格重划的基础。

6.1单元质量的评定及网格自适应技术
理想的网格的单元应该是等边三角形、正方形、等边四面体和立方体。

但是对于任意的复杂的几何形状结构，试图用完全的理想的单元去离散和描述是徒劳的。

所幸的是，实际情况的要求并不如此的苛刻。

实际的单元只要与这些理想的单元形态足够的接近，就可以获得能够接受的分析结果。

评定单元几何形态质量的量化标准如下]71[：
单元边长比(Aspect Ratio)：是单元最长边与最短边之比。

理想的单元边长比是1。

可接受的单元边长比的范围是：
AR<3对线性单元，如三节点三角形、四节点四边形、四节点四面体或八节点六面体单元。

AR<10对二次单元，如六节点三角形、八节点四边形、十节点四面体或二十节点六面体。

此外，非线性分析对单元边长比的要求比线性分析高。

扭曲度（Distorsions ）：是单元在单元面内的扭转和单元的面翘曲程度的指标。

对三角形单元，扭曲度用相邻夹角与0
60之间的差别定义；对四边形单元，扭曲度用单元相邻边的角度与090之间的差别描述。

当单元面的节点不共面时，就发生面外翘曲。

网格疏密的过渡：网格疏密过渡时要求单元和节点必须匹配，用连续的网格描述单元之间的连接。

方法一：是采用单元边长的渐变，方法二：是采用节点区域的加密。

如图6.1所示
图6.1网格疏密过渡的两种方法
有限元分析的精度和效率不仅与单元的几何形态有关而且和单元的密度之间存在着密切关系。

对于每一次有限元分析，我们总希望以合理的建模和计算时间，获得最理想的计算结果。

有限元分析结果的精度与离散模型的网格划分是密切相关的。

工程问题结构形状和边界条件往往十分复杂，初始建模划分的网格并不一定保证结果计算精度和计算效率都足够高。

显然，过密的网格可能会造成计算费用的大增，而过疏的网格又无法精确描述场变量的空间变化；另外，初始预定的网格划分很难适应在不同时间点上变量的空间分布变化。

根据误差识别，能够自动调整网格疏密的网格自适应技术，成为以合理费用，提高复
杂问题计算效率，改进结果精度的有效措施。

自适应网格技术是以某种误差判据为依据的。

一旦误差准则在指定的单元中被违背，这些单元会按指定的单元细化级别在指定的载荷增量步内被细化。

常用的误差准则有]71[：
平均应变能准则。

当单元应变能大于系统平均应变能的指定倍数时细化。

Zienkiewicz_Zhu 应力误差准则。

定义计算应力与磨平应力的误差为准则。

Zienkiewicz_Zhu 应变能误差准则。

定义计算应变能与磨平应变能之差为判断准则。

在一给定区域内的节点误差准则。

落入所划区域那些节点所在的单元细化。

网格自适应技术是有限元分析中一种难度较大的技术，就目前而言，网格自适应技术还不够完善，但由于其在有限元分析中的重要地位，许多研究者都在积极探索这一领域。

目前，国外的一些有限元软件，如美国的MSC 公司Marc 提供了这种技术。

6.2塑性加工有限元法中的网格重划技术
由于变形体形状的复杂性，在塑性加工中采用的单元多是等参数单元，随着塑性加工过程的进行，初始定义的网格将会发生畸变，严重时会发生重叠，和模具发生穿透干涉，使得计算无法准确地进行下去，因此，当计算网格畸变到某种程度时必须停止计算，实施网格重划。

通常确定网格是否要进行重划的准则主要有以下几种。

1.干涉准则
假设一个单元的一边穿进了模具，当干涉量大时就会导致较差的计算结果，这时需要进行网格重划。

如图6.2，设P 是单元与模具发生干涉一边的中点，Q 是对边中点，i h 是P 点在PQ 方向上模具表面的距离即为线段PO 的长度。

e h 是PQ 之间的距离。

那么该标准可表示为
i e i C h h ≥ （7.1） 其中i C 是用户确定的干涉标准常数，一般取值在0.01～0.1之间，当比值大于该指定常数时就必须进行网格重划。

这种方法也可以适用于平面三角形单元、空间四面体、六面体单元。

Die
P
O
Q
图6.2干涉准则
2.畸变角准则
如果采用的是等参数单元，为了保证母单元计算结果映射到整体坐标系下的单元保持一一对应，必须满足
式中J 为坐标变换的雅可比(Jacobian)矩阵
在四节点四边形等参数单元中，假设i θ为四边形的任意内角，则上式等价于
即四边形单元必须保证外凸性，才能使得相应的变换维持有效性。

在实际的计算中，由于内角接近于00或0
180时计算的精度都很低，应缩小内角的取值范围，通常取
对于其他种类的四边形单元，理想的内角值为090，允许的偏差常采用的值为040。

对于三角形单元或四面体单元，理想的角度为060，允许的内角偏差为040。

3.增量步准则
按指定的增量步间隔进行网格重划分。

这种准则带有一定的盲目性。

在实际的运用中，上述的几种准则可以单独使用，也可以将其组合后使用。

6.3新网格的生成
在旧网格系统上生成新网格是网格重划技术的关键环节。

新网格的生成从总体上可分为整体重划和局部重划两种。

这两种方法各有利弊：整体重划通过网格的自适应技术有利于实现整个分析过程的自动化，使得整个分析过程在不停机的情况下自动完成，但在新旧网格之间的场量数据的转换将要花费大量的计算时间。

局部重划只是针对在变形过程中网格畸变较大的局部区域，对网格质量较好的区域不予处理，这样在数据的转换的工作量上是经济的。

但是，采用局部重划一般是采用计算机交互式绘图功能，通过人机对话在屏幕上显示，进行修改网格和产生新网格。

不难看出这种方法往往需要在分析过程中停机的情况下才能实现，难以实现分析过程的全自动化。

但是，无论采用哪种网格系统的生成方法都应使得新网格系统满足：
（1） 变形体的边界不变，以真实地反映变形过程。

（2） 网格系统下的单元形状良好。

（3） 尽量使半带宽小，以减小计算工作量。

6.4新旧网格系统的场量数据传递
重新划分网格后，为了保证分析的连续性和准确性，必须将旧网格的场量数据传递到新网格上。

信息的传递是网格重划的核心。

需要传递的场量信息有：节点速度场、变形历史积累的场变量，如：等效应力、等效应变、等效应变速率等。

数据信息传递必须准确、可靠，否则会使后续的分析计算失去意义，整个分析失败。

常见的数据转换的方法有：
在已经变形的变形体上覆盖一标准网格，作为数据传递的中间网格，把旧网格上的场量插值到中间网格上，然后再把中间网格上的信息插值到新网格上。

不难看出如采用该种方法，场量数据从旧网格到新网格将需要经过多次转换，计算量大。

该方法是把旧网格节点值直接转换到新网格节点上，最后产生新网格的场变量值。

对于二维问题，最常用的是面积加权平均法，即按相邻单元面积大小进行面积加权平均，对于新网格节点，可根据它在旧网格中包围它的单元面积即场量进行加权平均，即 j n j j n j j i A A
∑∑==1011
εε＝ （7.2）
式中j A 为就网格中包围新节点的面积，0j ε为旧网格j 单元的场变量值，i ε为新网格节点
场量值。

在变形体上设置跟踪点，一般取网格初始节点作为跟踪点，随着变形过程的进行跟踪点也随着变形体变形，当网格畸变至需要重划时，跟踪点仍存在于新网格中，新网格场变量认为是新网格单元所包含的跟踪点之场变量的平均值，而跟踪点的场变量值是通过每次变形的场变量之增量的积累求得。

上述的几种方法，各有其自身的优缺点，本文在综合一些方法的基础上总结出一种从处理上较为简单且适用范围广的方法。

有限元解题的核心思想是将研究对象离散化。

利用有限单元法离散化特征，可以将畸变后的单元进一步细化，进行分片插值，具体的方法如下：
首先，将旧网格单元的积分点的有关场量外推至旧网格节点，获得单元节点的场量。

然后对旧网格的单元进行三角形的细化处理。

也就是说将每一个二维的四边形或三角形单元都被划成更小的三角形单元；每一个三维的四面体、五面体或六面体单元都被离散成更小的四面体单元。

用经细化的旧网格的三角形角点坐标来描述新网格上任一节点的空间位置。

通过插值，不难获得新网格单元节点的节点场量和单元积分点的状态变量。

以平面四节点四边形等参元为例，说明这种方法的实现步骤。

1.在旧网格系统中插值出四节点四边形等参元形心的整体坐标及其各种场变量
不失一般性，以等效应变的分布为例，对于四节点等参元有
∑=4
1
041i i ）（＝εε （7.3） 式中ε为四节点等参元形心处场量值，0i ε为四节点等参元高斯点处场量值。

2.将旧网格积分点（高斯点）上的场量值外推至单元节点
几乎在所有已知的插值方法中，面积加权平均法最为简洁、方便，并对于四节点四边形等参元而言，由于采用双线性插值函数，故可用单元形心处的场量值来表示整体单元的场量。

如图6.3所示，节点P 由相关的单元包围,相关的单元可以这样确定：在全域内对所有单元的节点号进行搜索，如某单元中有节点P 则该单元为相关单元。

一般来说内部节点有四个相关单元，而边界处相对较少，从中也可看出这种方法在变形体内部的精度较边界上高一些。

假设所有相关单元的场量已知，节点P 点场量待求，则
图6.3 面积加权平均法求旧网格节点上场量
ip ip m i ip
p A A
εε∑==11
（7.4） 式中p ε为待求节点的P 处的等效应变值，ip ε为节点P 的第i 个相关单元的已知等效应变值，m 为节点P 相关单元的总个数，ip A 为节点P 的第i 个相关节点对节点P 的面积贡献，且定义如下 ηξηξηξd d J N dA N A i A ip ⎰⎰⎰--=1111),(),(＝ （7.5）
式中）ηξ,(N 为局部坐标系下节点P 处的形函数值。

3.判断新网格节点处于旧网格中哪个三角形单元中
在旧网格上，将每个四节点等参元绕中心划分为四个三角形单元，如图6.4所示。

图6.4 四边形单元绕形心分成三角形单元
对于三角形单元，二维平面中任意点),(y x i 在三角形三个定点的形函数的和为1，如图
6.5所示即
1=++C B A N N N （7.6）
式中 C B A N N N ,,分别为点i 在三角形三个顶点的形函数值。

图6.5 判断新网格节点),(y x i 是否在三角形中
若点i 位于三角形内部，则这三个形状函数的值均介于[0,1]之间；否则形状函数的值或大于1，或小于0。

根据点于ABC ∆中形状函数的特点，可以将
[][][]101010，，，∈∈∈C B A N N N （7.7） 作为判断一点是否在三角形中的判据。

当节点确定的三角形三个顶点的形函数值均介于[0,1]之间时，则节点在三角形中或在三角形的边上。

在实施过程中为了减少程序中采用过多的判断语句，将判据等价为 5.05.0≤-K N C B A K ,,= （7.8）
在编程判断时，由于计算机浮点计算，它的0值并非精确的0值，而可能是一个绝对值非常接近于0的正负数，如100.1-±。

因而判断节点i 是否处在三角形中的判据为 500001.05.0≤-K N C B A K ,,= （7.9）
4.将旧网格节点上场量插值至新网格节点
当节点i 处于ABC ∆中时，则可以通过三角形单元内插值出节点i 的场变量值。

Cy C By B Ay A iy u N u N u N u ++= （7.10）
用上式插值出新网格的所有节点的速度场之后，整个新网格的场量分布也就知道了。

新网格的任一节点的场量都可以用这样插值得到。

上述的方法不仅适用于二维平面问题的各种单元，对于轴对称问题即相当于把子午面看成平面来处理；对于三维的问题也同样适用，只不过是将空间单元的形状离散为四面体单元，将面积加权变成体积加权，其余的过程同上相似。