第 6 章：数值求解方法
尽管对控制方程进行了简化，但得到的方程仍然复杂，需要采用数值方法求解。主要的数值方法包括有限差分法和有限元法。在本章中，我们将讨论这些方法在求解这些方程中的应用。关于数值方法的一些历史评论可以在附录A中找到。需要注意的是，在商业产品中采用的数值方法与当时的计算机性能密切相关。因此，在审视各种数值方案时，应了解过去CPU和内存的限制。有限差分法和有限元法在附录E和F中进行了介绍。这些附录旨在提供对这些方法的基本理解。一种特定的商业数值技术必须从其准确性、创建几何模型所需的工作量、获取相关材料数据的能力以及实施所需的计算资源等方面进行评估。学术代码通常仅关注这些方面之一，如果受到限制的话。计算资源的一个衡量标准是节点上的自由度（DOF），这指的是在计算域中每个感兴趣点必须由模拟计算出的变量的数量。在接下来描述的每种方法中，我们将讨论这一衡量标准。
6.1 中平面方法
在附录G中详细推导了注塑成型过程填充、填充结束和冷却阶段的有限差分和有限元方程。这通常被称为2.5维方法，因为在中平面模型定义的x和y方向上，压力在二维中求解，而温度在三维中确定。这一概念由Hieber和Shen [154] 提出，从1980年至1997年，这种方法一直是首选方法。2.5维近似被广泛采用，因为它能够与基于表面的CAD系统接口。在这一时期，三维建模并不常见，线框和表面模型是主流，而中平面方法正好适用于使用此类CAD模型。2.5维方法的优势在于使用有限元分析来确定压力，使用有限差分来确定温度，并简化了只需要在每个节点知道压力梯度和流度s²的动量方程。在该方法提出之初，它是一个极佳的解决方案，至今仍然如此。压力在模具成型过程中变化较为平缓，而温度在模具厚度方向上则有剧烈变化。因此，对于温度场的有限差分求解是理想的，因为一旦为中平面生成了有限元网格，用户可以根据温度解的节点数量来决定厚度方向上的细节程度。
因此，每个节点的自由度包括压力 p、流动性 S2 和温度 T。适当的中平面模型通常能提供良好的分析结果，可用于填充、充模和冷却分析。完成这些分析后，可以计算残余应力或应变，进而进行结构分析，提供收缩和翘曲结果。从商业角度来看，关键因素在于可以使用同一个模型进行所有分析。对今天的读者来说，这似乎显而易见。但在缺乏标准和计算资源的情况下，实现单一模型的概念变得困难。中平面方法落后于3D CAD软件的发展。在20世纪80年代，CAD领域发生了革命。由Parametric Technology 领导，定义制造对象的方法以照片级真实的3D表示形式呈现。这在CAD模型与其在注塑成型分析（或其他任何类型分析）中的表示之间造成了鸿沟。CAD行业围绕单一模型的概念展开——工具、装配、分析和生产。对于注塑成型分析，由于在CAD发展过程中滞后，主要问题是如何从真正的3D模型中提取中平面表示。这导致了三种克服这一困难的方法：
- 从3D模型中提取中平面
- 双域分析
- 全3D分析
在这三种方法中，前两种基于2.5D近似。最后一个项目，3D分析，与前两者截然不同。我们将在接下来的内容中讨论这些方法。
6.1.1 从三维模型中提取中面
从三维几何中自动生成中面并非一项简单的任务。图6.1展示了这一概念。挑战在于从图中左侧所示的三维几何中自动推导出具有厚度定义的所有元素的中面网格。虽然图6.1中的示例相对简单，但图6.2所示的更现实的部件则更为复杂。该领域中的许多学术工作集中在获取三维几何的中心轴上，这主要是由于图像分析和六面体网格生成的需求。
图 6.1：中面网格的生成
图 6.2：复杂注塑件的三维表示
许多方法源自图像处理领域，试图找到对象的中心轴[82]。这些方案的一个缺点是计算中面所需的时间通常为数小时。更为严重的问题是，生成的中面网格往往不如原始几何体平坦。虽然这对流体分析影响不大，但对用户来说却显得不美观。更糟糕的是，这种方法未能保留部件的结构特征，因此无法进行结构分析，进而无法进行翘曲分析。肯尼迪和余[198]提出了一种不同的自动中面生成方法。其出发点是用于立体光刻的原始网格，当时这是一种流行的方法，用于生产原型部件。该网格由许多平面三角形组成，这些三角形通常具有非常差的长宽比，定义了实际三维部件的外部表面。
这些网格经过重新划分以改善纵横比，然后被分组为表面。开发了算法来确定表面应坍缩的方向。在坍缩过程中，从原始位置到中平面的厚度变化被分配给结果的中平面单元。这种方法对于平面薄部件有效，但在处理如凸耳等小型特征时存在一些困难。生成的网格经常出现需要手动修复的不连续性。手动修正不具有商业可行性，因为它涉及用户交互，成本较高。商业上实用的算法依赖计算机而非人力。
6.1.2 双域分析法用于流动分析
而不是从三维几何体中确定中平面，另一种方法是将三维几何体转换为等效的二维半几何体。这导致了双域方法的产生。这种方法的起点是一个在三维几何体上的外部网格。这通常是一个立体光固化网格，如前一节中提到的自动中平面程序所进一步细化的。图6.3展示了基本思路。考虑一个在中心注入的矩形板的横截面，如图6.3(a)和(b)所示。图(c)显示了板的横截面中的流动；图(d)展示了使用连接单元来确保与图(b)中真实流动的物理一致。如果为表面单元规定一个厚度，就可以在表面网格上进行二维半分析。
图 6.3：双域流动分析，(a) 表示向矩形板中心注射，(b) 表示板横截面中的流动，(c) 表示表面网格上的流动前沿推进，(d) 表示使用连接单元以确保与 (b) 所示真实流动在物理上一致
然而，这样的分析在物理上是不一致的，因为材料会在顶表面流动，绕过边缘，然后沿着底表面在注射点下方形成焊线，如图6.3(c)所示。解决方法是在注射点将顶面和底面网格连接起来，如图6.3(d)所示。这样，材料就会同时沿着顶面和底面流动，如预期的那样。这导致了双域有限元分析（DD/FEA）这一名称。实际上，我们在表面网格的两侧分别进行了两次分析。由于我们实际上填充了两个区域，因此需要将注射点的流动速率加倍，以获得给定几何形状下使用中面网格计算的填充时间。实际的零件总是比简单的平板更复杂。考虑一个带有两个肋的平板的横截面，如图6.4所示。当注射点连接到顶面和底面时，流动会从该点发出并撞击每个肋，如图6.4(a)所示。在肋处，流动仅沿一侧继续，如图6.4(b)所示。再次需要将对面的表面连接起来，以便流动沿肋上升，并继续经过肋，如图6.4(c)所示。虽然这个想法很简单，但其实现却相当复杂。不过，基于这一想法的产品于1997年由Moldflow发布。该概念在美国[411]和欧洲[410]获得了专利，并在其他司法管辖区等待授予专利。
图 6.4：带两个肋的部件的双域流动分析
总结来说，双域流动分析包含三个基本步骤：
- 生成表面网格
- 在上下表面之间建立关系，以便定义厚度
- 添加连接元件以保持物理上现实的流动模式
详细实现细节参见 Yu 和 Thomas [411]。在节点自由度（DOF）方面，我们有压力 p、流度 S₂ 和温度 T，与中间平面情况相同。然而，由于我们需要在模具的两侧分别求解两个网格，因此总自由度数大约是中间平面情况的两倍。因此，双域方法确实需要更多的计算时间。不过，这仅仅是计算时间，相对于人工操作员准备模型进行中间平面分析的成本而言，相对较为经济。这一点澄清了学术目标与商业目标之间的区别。在学术环境中，双域方法很可能从未被考虑过。
6.1.3 双域结构分析
双域概念被工业用户迅速采纳，以至于立即对软件提供商产生了扩展该概念到模具冷却和翘曲分析的压力。模具冷却分析很容易扩展。模具的边界元方法得以保留，并与塑料材料内部的表面网格所包围的热传导分析相结合。为了允许进行翘曲分析，首先需要具备结构分析能力。
以下描述大多源自[113]。这种方法最好通过一个例子来介绍。核心思想是仅使用定义板外边界的空间网格来模拟板的结构性能，即板的弯曲和膜特性。从几何角度看，厚度为h的平板可以视为两个厚度为h/2的板完美粘合在一起。如果我们考虑这样的组合，可以发现它可以用两个壳体来建模，每个壳体的参考表面位于两个板的几何中心。然而，这存在一个问题，因为定义中面的节点偏离了外表面。我们希望在不修改的情况下使用外表面的网格。这可以通过使用偏心壳体单元来实现。如图6.6所示，可以定义一个结构壳体单元，其参考平面可以位于任何位置。从中面到参考平面的距离称为偏心距。回到这个问题，我们可以看到，在图6.5中，上板可以使用其顶表面作为参考表面的偏心壳体单元来建模。同样，下板可以使用其底表面作为参考表面的偏心壳体单元来建模。这解决了问题，允许我们使用板外表面的现有节点。
图 6.5：一个简单平板可分解为两个部分，每个部分厚度为原厚度的一半，并彼此完美粘合
图 6.6：用于结构分析的偏心壳单元
然而，为了获得正确的结构响应，厚度为 h/2 的两块板必须以某种方式“粘合”在一起。顶部和底部板的粘合涉及施加经典板或壳理论中的 Love-Kirchhoff 假设 [83]，要求板或壳的法线在变形后保持直线且长度不变。这通过使用多点约束来实现。总之，对于结构分析，双域方法涉及以下步骤：
- 对结构的外表面进行网格划分，并在顶部和底部表面之间建立元素关系以定义局部厚度；
- 使用参考表面位于定义三维对象外边界表面的壳单元；
- 使用多点约束确保顶部和底部表面的法线在变形后保持直线。
所使用的约束取决于所选择的单元类型。在过去，3 节点三角形单元在整个单元内具有恒定应变，其性能较差。然而，由于对结构优化的兴趣，开发了许多改进的单元形式。例如，具有 18 个自由度（每个节点六个——三个位移和三个旋转）的平面三角形面壳单元。
该单元通过叠加Bergan和Felippa [34]提出的局部膜公式与Batoz和Lardeur [28]提出的弯曲公式构建，并将结合方程转换到全局坐标系中。膜公式中使用了关于局部参考表面法线的钻孔旋转自由度，该自由度在局部单元坐标系中定义为：
为了定义节点n与其匹配节点p之间的自由度关系（见图6.7），我们要求变形前中面的法线在变形后保持直线。采用单元的局部坐标系，我们分别表示节点n的三个位移自由度和三个旋转自由度为
图 6.7：为双域分析配对的结构单元
其中，h是节点n与其匹配点p之间的距离。注意，式(6.7)中的关系是通过式(6.1)获得的，因此与所选择的单元类型有关。这种约束系统在模型底部（或顶部）表面上的所有节点上施加，但不包括边缘节点。
板边缘的元素被分配相邻上下面元素厚度的六分之一。在这些约束条件下，复合结构的结构性能与原始板相同。此时，复合模型可以施加适当的边界条件和载荷进行结构分析，并用于分析3D几何的翘曲。通常，上面的网格与下面的网格不重合。因此，从下面表面节点n出发的法线通常不会与上面表面的节点重合。相反，法线更有可能在上面元素内某点p与上面元素相交（见下面元素节点n，该节点在上面元素内某点p处与上面元素相交）。在这种情况下，需要进行插值。
图 6.8：顶面和底面的单元通常不重合，即底部单元节点 n 的法线会在顶部单元内的某点 p 与其相交，此时需要插值
6.1.4 双域有限元法的翘曲分析
为了使用双域方法进行翘曲分析，我们需要加载有限元模型，施加从注塑过程填充、填充和冷却阶段分析中得出的收缩应变或膜应力。如第4.8节所述。与中面模型相比，使用双域方法会带来轻微的计算开销。考虑到创建中面模型所需的时间，这种增加的求解时间可以忽略不计。因此，从总求解时间来看，该方法极其高效。
双域结构分析的专利已在美国 [114] 及其他几个司法管辖区获得授权。有关该方法实施细节的详情请参见专利文件。
6.2 三维分析
上述内容均涉及二维半近似。特别是，我们讨论了获得中平面模型或使用双域方法的需求。另一种方法是进行三维分析。这种方法避免了二维半近似的假设，原则上应是最终的模拟方法。 特别是，三维分析符合三维CAD建模对象的趋势。与双域方法使用有限差分和有限元相结合不同，大多数三维代码基于单一数值方法。有限元代码最为常见；我们将在附录H中提供元素形式化的描述。一些有限差分和有限体积代码起源于铸造模拟。Hetû等人 [152] 提供了最早的三维注塑成型分析示例。随后是Pichelin和Coupez [298] 和Talwar等人 [353] 的研究。这种分析假设较少；特别是，没有对零件厚度的限制，但需要将计算域用四面体或六面体单元进行网格划分。四面体网格划分通常更受欢迎，因为它可以自动完成。然而，由于塑料的低热导率，注塑成型零件往往较薄。
在厚度方向上，每毫米数百度的温度梯度要求在零件厚度方向使用大量单元。因此，在薄壁零件中，用于三维分析的单元数量急剧增加，导致计算时间过长和资源过多。为解决这一问题，采用了各向异性网格 [334]。不出所料，三维分析的计算需求高于任何其他方法。在每个节点上，我们有压力
6.2.1 有限体积法
迄今为止，我们的讨论主要集中在有限差分法和有限元法上。Chang 和 Yang [55] 描述了一种用于三维注塑成型模拟的有限体积法。这种方法的一个优点是可以使用不同的单元形状进行模拟。例如，可以在零件内部使用四面体单元，而在靠近壁面的地方使用矩形棱柱单元以捕捉温度变化。虽然有限体积法提供了一定的灵活性，但网格生成变得更加复杂，如 Chang 等人在 [54] 中所述。
6.2.2 半三维方法
Nakano [265, 266] 开发了一种并非真正三维的方法，但允许直接分析三维几何结构。
以笛卡尔坐标系中的连续性方程为例：
Nakano 对薄壁近似的二维情况进行了推广，并设定了以下速度分量：
其中，
需要注意的是，上述方程的右边为零。然而，如果假设材料是可压缩的，实际上它确实是可压缩的，那么右边不为零。因此，
其中，
Nakano 方法节省的经济性体现在每个节点需要确定的自由度（DOF）上。忽略约束条件，每个节点需要确定一个压力、一个温度和一个流体性
6.3 三维中的翘曲和收缩分析
由于三维域使用了三维单元（如四面体或六面体单元）进行网格划分，为了计算收缩应变和翘曲，我们需要将应力或应变输入到结构分析软件中。这并不是一个简单的任务。
一种粗略的方法是将计算出的体积收缩率的大约三分之一作为应变施加到每个单元上，并将其用作加载条件。对于真正三维的零件，如果模具材料的各向异性不明显，这种方法可能还不算太差。然而，一个显著的问题是，对于真正三维的零件，很难确定模拟FEM模型中哪些部分是由模具刚性约束的。虽然零件的某些区域可能被约束，但在材料冷却和收缩时，其他区域可能自由地从模具壁分离。因此，所使用的边界条件并不简单。如果分析模型中包含某种厚度概念，例如在使用各向异性网格时，可以引入某种方法来分配由约束和流动诱导结构产生的线性应变。当然，这些应变必须与计算出的体积收缩率一致。对于半结晶材料，问题更加复杂。真正三维的零件很少具有均匀厚度，因此在加工过程中会有不同的冷却速率和材料变形。这将导致模具成型过程中材料性质的变化，从而导致不同的收缩率。最终结果将是复杂的各向异性收缩分布，导致零件变形。
6.4 流道系统的三维分析
在本节中，我们将讨论流道系统弯头处的热量对流，并展示它如何导致多腔模具填充不平衡。
我们还将展示，三维分析对于捕捉这些效应是必要的。如第5.10节所述，对于流道的2.5维近似做出了若干重要假设。首先，假设流道是圆形的。即使流道的横截面是非圆形的，也会根据某些规则近似为圆形，通常采用水力半径。这可能会在流道中的流动和温度计算中引入显著误差。Beaumont [30] 发现，在实践中，流道中的温度对流可能会导致多腔模具中不同部件的填充不平衡，甚至在流道系统拐角处。如图6.9所示，假设流道系统中存在剪切加热，拐角处的对流会导致不同温度的熔体进入四个腔室。由于聚合物在较高温度下会更稀，较高温度的熔体会优先流入外腔，从而导致流动不平衡，进而可能导致模具部件的重量和性能不平衡。由于流道的2.5维近似假设径向温度相同且不考虑拐角修正，因此无法预测这种效应。在这种情况下，由于剪切加热导致的温度差异可能会导致尽管进料系统自然平衡，但腔室填充仍不平衡。
图 6.9：在三维分析中，温度能围绕流道内的转向正确对流。在这些情形中，剪切生热导致的温度差异可能使型腔充填不平衡，尽管进料系统本身是平衡的