用于在模具填充过程的模拟中描述粒子统计取向分布的方法和装置

摘要

一种用于在过程的模拟中描述非球形粒子的统计取向分布的方法和装置，在该过程中用包含大量非球形粒子的悬浮液填充模具型腔。可以将所述方法和装置应用于制造纤维增强模制聚合物部件的注塑成型过程的分析或制造纤维增强金属产品的金属铸造过程的分析。这些分析的结果可以用于确定部件的张力和翘曲的方面，并且用于优化在制造过程中使用的过程条件。

说明书

技术领域

本发明涉及对进入模具型腔的包含粒子的悬浮液的流进行三维建模的领域，并且更具体地，涉及用于在过程的模拟中描述非球形粒子的统计取向分布的方法和装置，在该过程中用包含大量非球形粒子的悬浮液填充模具型腔。

背景技术

注塑成型过程或者金属铸造过程的真实3-D模拟涉及多达几十万个方程的系统。在过去已经取得了进展来改进模拟方法的效率以处理这些复杂的计算。使用优化的软件以及现代工作站的处理能力，可以在工作场所执行这种模拟，即以足够快的速度获取结果以适应纯科学研究领域之外的领域，并且可以由注塑成型物件的研究和开发部门、铸造厂和制造厂中的工程师应用。

当将粒子(诸如纤维)被加入到聚合物合成物，并且需要描述纤维的取向分布时，模拟和相关方程组明显变得更加复杂。由于计算时间太长或者模拟的精度不够，由于模拟的复杂性不允许在现今的工作站得到可接受的结果，因此迄今还没有在工作场所成功地引入这种过程的真实3-D模拟。

在纤维增强部件中，通常对于开发工程师来说，具有纤维取向分布的描述以便能够预测该部件的张力和翘曲方面是至关重要的。通常纤维用于改进塑料部件的机械性能。但是因此(热)机械性能(如热膨胀、强度和刚度)取决于纤维的取向。

近年来在许多行业中注塑成型塑料部件的使用已经稳定地增长。与以前相比，电子设备、消费品、医疗设备和汽车零件的制造商正用塑料制造越来越多他们的产品和在他们的产品中使用的部件。

因为注塑成型纤维增强部件提供提高的强度/重量比、耐用性、部件集成度和较低的成本，因此它们正在取代结构金属部件。

同时，竞争压力正驱使塑料注塑成型工业中的制造商去寻找优化设计的新方法以更好地使设计匹配制造过程。当在设计开发过程中较晚发现了需要修改部件或者模具构造时，用于实施必要的改变的延迟和相关成本明显高于在设计开发阶段的较早阶段。

想保证他们的部件是可生产的并且将最优地执行的公司想要使用计算机辅助工程技术来对注塑模具中的复杂流和所得到的纤维取向进行模拟和建模，以便在设计阶段早期更好地理解制造过程并且将这种认知结合到部件设计中。

当设计用于纤维增强部件的注塑模具以及要在其中制造的纤维增强部件时，应该考虑若干因素。诸如整体部件几何形状、最小和最大壁厚度、模具中的通过其注入液态聚合物和纤维悬浮液的门的数量和位置、模具中的型腔中的气体通过其排出的排气口的数量和位置、聚合物成分和性质、纤维性质和数量、收缩、容差和纤维取向分布的参数是其中的一些。由于紧密相互关联的关系，部件和模具设计不能可靠地纯粹基于终端部件的形式和功能，而是应该也考虑制造过程的影响。

计算机辅助工程模拟可以有利地用来为设计和制造工程师提供关于在注塑成型过程期间在模具型腔内可能发生的情况的视觉和数值反馈，允许他们更好地理解和预测预期部件设计的行为，以便可以基本上消除制造的传统、高成本的反复试验的方法。计算机辅助工程模拟的使用便于在设计阶段期间优化部件设计、模具设计以及制造处理参数，在该设计阶段中能够以最低的成本和对进度的最小影响来容易地实施必要的改变。

CAE模拟技术在用于纤维增强部件的工程过程中的应用包括：(i)“注塑成型”制造工艺的模拟，其包括流体流和热传递的计算；和(ii)压力和强度(还可能有耐用性)计算，这些都是在宏观水平上针对这些部件执行的，以确定它们在外部负载下的功能机械性质。两种模拟类型都需要适当的材料模型，该材料模型描述包含浸入的纤维的聚合物材料在液态以及固态下的特性。

宏观水平上的长度尺度由部件的几何形状的线性尺寸(整体大小、壁厚度等)确定，该线性尺寸通常在几mm到cm的范围内变化。计算单元的尺寸必须以足够的精度分解宏观长度尺度，因此它们通常比最小宏观尺寸小到一个量级。由于浸没在短纤维增强部件中的纤维的典型尺寸低于宏观计算单元的典型尺寸一个或者两个量级，因此通过统计方法来描述与宏观材料行为的建模有关的纤维性质。对于短纤维增强材料，相关的宏观性质是：(a)体积浓度，其关于整个零件通常(近似)恒定，和(b)在每个计算单元内的纤维取向(FO)的局部分布，其通常跨部件几何形状明显变化。(在详细描述的1.1和1.2节讨论此问题的进一步细节)。

通过相应分布函数的低阶(即二阶和四阶)矩提供局部FO的统计分布的简化的、出于实际目的的适当描述。由于它们的数学结构，将这些矩表示为(分别为二阶和四阶的)取向张量。在用于纤维增强部件的CAE模拟的构架内，需要四阶张量来预测纤维增强材料在宏观水平上的流变和机械性质，因为这些是四阶张量性质。二阶FO张量是具有单位迹的实值对称3×3矩阵，因此其9个分量中只有5个是独立的。通过(全)对称将四阶FO张量的独立分量的数量从3⁴＝81减小到15。

以其矩的形式描述FO分布的数学模型通过根据闭合关系以二阶张量的形式近似计算四阶取向张量而被显著地简化。闭合关系以函数的形式提供这种计算方案的数学描述，并且如果闭合关系仅仅在特定假设下近似有效，则将相关计算过程称为“闭合近似”。只使用二阶FO张量和某些闭合近似的方法导致“Folgar-Tucker”型的模型以模拟在模具填充过程期间二阶FO张量在时间和空间内的演化(详见详细描述的2.2到2.5节)。

文献“Glass fibre orientation within injection moulded automotivepedal Simulation and experiment studies，B.R.Whiteside等人，Plastics，Rubber and Composites，2000，卷29，No.1”公开了用于对增强热塑料物件内的纤维取向分布进行建模的方法，其用于非对称热塑料流，并且分析包含纤维取向预测算法。该软件使用由线性三角单元组成的二维有限元网格来近似三维模型。使用广义Hele-Shaw近似和Folgar-Tucker方程的变型来计算流场以计算纤维取向。在穿过每个单元的“厚度”的19个夹层(laminate)上使用有限差分技术来进行纤维取向、温度和粘度计算以产生三维解。然而，重要的是应该注意到由于该模型不能模拟在平面之外方向上的速度分量(Hele-Shaw近似的限制)，因此不能将该系统描述为真实的三维。当适用于真实三维模拟时，在该文献中描述的方法将导致不稳定以及过多处理能力消耗的模拟，该模拟不能在例如开发工程师的工作场所中使用。

发明内容

在此背景下，本发明的目的是提供用于在过程模拟中确定非球形粒子的取向分布的方法，在该过程中，用包含大量非球形粒子的悬浮液填充模具型腔，与现有技术方法相比，该方法更加稳定并且计算强度较小。

根据权利要求1，通过提供用计算机实现的方法来实现该目的，该方法用于通过使用模拟注塑成型过程的模拟模型来计算非球形粒子在宏观水平上的取向统计，在该过程中，形成模拟区域的几何形状的至少部分的模具型腔被由包含大量非球形粒子的溶剂形成的悬浮液填充或部分地填充，其中提供了模拟区域的几何形状的数字表示或者计算机模型，并且其中通过对模拟区域的至少部分进行细分或者离散化来形成具有多个计算单元的网格，该方法包括步骤：a)指定边界条件；b)设置初始条件；c)针对模拟区域的至少部分单元求解质量、动量和能量的平衡方程以获得宏观水平上的流体流、热流和质量传递；和d)至少部分基于所解的平衡方程的结果来求解非球形粒子取向动力学方程，以由此确定宏观水平上非球形粒子取向的变化作为空间和时间的函数。这里，对于步骤d)，可能仅仅针对包含悬浮液的计算单元求解粒子取向方程。

优选地，步骤c)还包括(cc)：至少部分基于所解的平衡方程的结果确定流体或者悬浮液的更新的自由表面或者流前沿(flow front)，其中该自由表面将型腔中由悬浮液填充的单元与空单元分开。优选地，步骤(cc)还包括根据所更新的流前沿更新边界条件。同样优选地，本发明的方法还包括步骤e)：通过确定模具型腔是否被悬浮液的模拟注射填满来确定是否完成了模拟注塑成型过程；和f)：重复步骤c)、cc)和d)，直到完成模拟注塑成型过程为止。

根据权利要求6，通过提供用于在过程模拟中描述非球形粒子的统计分布取向的方法，也实现了本发明的目的，在该过程模拟中，用由包含大量非球形粒子的溶剂形成的悬浮液填充模具型腔，该方法包括：(1)提供定义型腔的几何形状的三维计算机模型；(2)指定边界条件；(3)基于该模型对解域进行离散化以形成具有多个单元的网格；(4)针对解域的至少部分，求解能量和流方程；(5)计算各自单元中的流和温度条件作为时间的函数；(6)计算非球形粒子取向的变化；和(7)将在各自单元中的非球形粒子取向的统计分布描述为时间的函数。

本发明的方法以显著降低的计算强度利用悬浮液流来预测遍及部件的纤维取向分布和热机械性质。开发工程师可以使用模拟的结果，并能够由此优化与纤维取向相关的产品并且因此改进该物件的强度和形状稳定性。

根据详细描述，根据本发明的方法和装置的进一步的目标、特征、优势和性质将变得明显，该方法和装置用于确定在过程模拟中的非球形粒子的取向分布，其中在该过程中，用包含大量非球形粒子的悬浮液填充模具型腔。

附图说明

在以下本说明书的详细部分中，将参考附图所示的示例性实施例更详细地说明本发明，其中：

图1是穿过包括模具的注塑成型机的图形表示的横截面视图；

图2是概括根据本发明的实施例的注塑成型过程的基本处理步骤的上层流程图；

图3是进一步详细图解图2的流程图的步骤5的流程图；

图4是纤维增强塑料物件的显微图；

图5图解了在本发明的实施例中使用的纤维的模型；

图6是图解特征多项式的二次和三次以及线性和常数项的曲线图；

图7是图解矩阵的三个特征值的图解；

图8示出了对应于“极大对称”情况的特殊示例；

图9示出了四面体形体积的扭曲形式；

图10示出了根据本发明的实施例的FO矩阵(即2阶FO张量)的相空间集合M_FT的总体结构的图；

图11是图解根据本发明的实施例的算子分裂过程是以简化形式的流程图；

图12是图解时间步长方法的细节的流程图；和

图13是图解根据本发明的实施例的使用迹尺度改变的相空间投影过程的流程图。

具体实施方式

图1图示性地示出了注塑成型机1。该注塑成型机设置有螺杆2，其被馈送设置在漏斗3中的聚合物粒料。聚合物粒料通过螺杆2和加热单元4的作用而转变成粘滞团，该粘滞团在高压下被推动到模具6中的模具型腔5中。并且成型机和注塑成型制造周期在本领域内是众所周知的，在此不详细说明。使用成型机1，可以制造非纤维增强和纤维增强塑料部件。

可以根据图2所图解的过程在计算机上进行注塑成型过程的数值模拟。

一般所认为的模拟的主要步骤如下：

-步骤1，提供模拟区域的几何形状的数字表示；

-步骤2，网格划分(enmeshment)，其将计算区域细分成很多小单元，这些小单元是用于将微分方程离散化(利用不同求解算法)的基础，并且以此方式找出要模拟的物理现象的解；

-步骤3，将对于不同的材料域的必要物理数据附加到模拟模型中(数据库或者资料库)；

-步骤4，指定针对模拟项目的边界条件；

-步骤5，模拟注塑成型过程(以下将更详细地声明该步骤)；和

-步骤6，在计算机(诸如工作站)的显示器上将结果显示为图形或者数值表示。

在图3的流程图中图解了当模拟纤维增强物件的注塑成型过程时的步骤5的细节。在过程的这个部分，使用数值算法来解热流、流体流以及压力和张力的微分方程：

-步骤1，设置热物理材料性质和流前沿的初始条件；

-步骤2，使用质量、能量和动量方程的守恒来求解整个域的热方程和所有流体单元的流方程；

-步骤3，在此步骤中移动流前沿并且根据新的流前沿采用边界条件；

-在步骤4中，根据在之前的步骤中获得的流速计算纤维取向和输运。

以下将更详细地说明该步骤。将初始条件应用于新填充的单元。仅考虑包含流体材料的单元；

-在步骤5中，计算例如化学反应的附加量，并且验证单元是否凝固；

-在步骤6中，验证模具中的注塑成型过程是否结束；如果没有结束，则模拟继续下一时间步长并且该过程返回到步骤2；

-在步骤7中，计算在从模具脱模之后的物件的性质。

使用从注塑成型模拟获得的信息计算在模具外的热物理材料温度，即温度、收缩、翘曲等。

图4是纤维增强塑料物件的显微图，其中可以看到在注塑成型过程结束之后的纤维取向。最终产品中的纤维取向很大程度上取决于注塑成型过程期间的热塑性团的流型。不是针对每个单个纤维精确地确定纤维取向，而是通过分布函数进行描述。

在进入关于纤维取向计算的细节之前，提供短纤维热塑性熔体的基本流体机械方面的概述。

1.短纤维热塑性塑料熔体的基本流体机械方面

聚合物团包括分散于其中的大量的短纤维，使得所制造的部件将由短纤维增强热塑性塑料制成。在熔化状态中，即当温度足够高使得热塑性塑料基体是液体时，塑料熔体和分散的纤维的混合物组成复杂流体，该复杂流体在流体动力学和物理学的术语中通常被称为粒子悬浮液。一般而言，这种悬浮液由两种不同的相组成：(i)溶剂，在我们的情况中对应于没有分散纤维的熔化的塑料材料，以及(ii)粒子相，其由所有浸没在该溶剂中的纤维组成。

1.1纤维几何形状和材料性质

在下文中将拥有旋转对称轴的球形粒子(如图5所示)称为纤维，并且除了另外明确指出以外，同义地使用术语粒子和纤维。通过纤维的长度l和直径d以及从这两个量得到的纵横比r_a＝l/d来表征纤维的几何形状。分散在短纤维增强热塑性塑料材料中的纤维(例如，碳或者玻璃类型)的这些量的典型值例如是l≈0.5mm、d≈5μm以及r_a≈100。通常这些值以相应的大纵横比值在l＝0.1到1mm、d＝1到10μm的范围内变化。

对于浸没在纤维增强塑料中的大部分类型的纤维来说，纤维材料的质量密度可与悬浮液塑料熔体的密度相比。塑料熔体本身在模具填充阶段期间的典型温度下具有相当高的粘性。通过结合这两方面，当描述在根据优选实施例的方法中的塑料熔体内的纤维运动时忽略浮力效应和惯性效应。

关于纤维浓度，假定贯穿悬浮液的纤维的空间均匀分布，并且因此认为浓度不变。可以改变的效果可以是：(i)已经在浇口处存在的不均匀分布的纯粹对流输运，或者(ii)在模具填充阶段期间的剪切诱导粒子迁移效应的出现。

如果(i)或者(ii)组成显著效果，则用于描述浸没在热塑性塑料熔体中的短纤维的悬浮液流的模型有必要地需要是“两相流”模型，该模型将粒子相和溶剂相的流彼此耦合并且允许不均匀粒子浓度。然而，情况并非如此。可以将现象(ii)理解为一种扩散过程，其中扩散常数与相对粒子大小的平方成比例，相对粒子大小定义为悬浮粒子的实际大小(即球形粒子情况中的纤维长度或者半径)与流型腔的典型尺寸(例如，壁的厚度)的比。由于即便对于通过注塑成型由纤维增强塑料制造的薄壁部件而言，相对粒子大小通常为0.1到0.01的量级，因此可以完全合理地认为在这种情况中(ii)完全可以忽略。关于(i)，可以想象由于在螺杆中的熔体的准备过程，在浇口处可能存在不均匀浓度轮廓。然而，针对典型短纤维增强塑料部件的实际实验性观察的纤维浓度在部件内的几乎各个地方仅显示非常小地偏离常数值，这可能是均匀纤维分布的假设的最佳论据[20]。

1.2纤维对流变性质的影响

热塑性塑料材料显示非牛顿流行为。在没有纤维的纯塑料熔体的情况中，可以通过标量粘度函数来对材料的流变性质进行建模，标量粘度函数取决于温度以及流(广义牛顿流体)的状态变量。尽管已知它们仅覆盖非牛顿流性质(如例如剪切稀释(shear thinning))的有限范围，但是已经证明它们对于注塑成型模拟的目的是有用的并且令人满意地精确。此类模型也用于根据本发明的优选实施例的方法中。

向热塑性塑料基体材料加入纤维产生的机械性质在固态下是强烈各向异性的，并且该机械性质很大程度上取决于嵌入的纤维所指向的方向的局部分布。原则上，如果材料在液(即熔化)态下，则也存在各向异性材料行为。为了解决这种各向异性，不得不用粘度张量取代上述标量粘度。已经存在若干研究对纤维悬浮液的材料行为的两种建模方式进行比较。对于如在模具填充中遇到的大部分流情形，已经发现通过两种类型的模型(即标量和张量模型)所预测的充填模式显示出几乎没有差别(参见例如[21]、[22])。因此当热塑性塑料熔体包含纤维时，忽略各向异性粘度效应并且在这里也使用简单的广义牛顿模型。这等价于忽略分散的纤维的取向对熔体的流变性质的影响。粘度取决于纤维浓度，但是由于假设纤维浓度不变(见上)，因此这个方面仅作为先验已知的材料参数输入，该参数有助于热塑性塑料熔体的有效材料性质。

由于这种途径，根据优选实施例的方法使用流和纤维取向的计算的部分解耦：尽管流速局部地影响分散纤维的取向，但是纤维取向对流的影响可忽略。因此独立于纤维取向的计算进行流的计算。以此方式，熔体的局部流速作为一组外部系数进入用于纤维取向计算的模型。

2.Folgar-Tucker模型

2.1 Jeffry方程

尽管悬浮在聚合物熔体中的纤维在它们的纵横比r_a＝l/d较大的意义上说是细长粒子，但是它们足够短以使得在溶剂的局部流场中作用在其上的机械力不能够引起任何实质的变形。因此这里将单个纤维建模为细长、旋转对称的刚体，其取向由沿着旋转对称轴指向的单位向量p给出。两个向量±p都表示粒子的相同取向状态。

图5图解了旋转椭圆形的刚性粒子，该粒子的运动受粒子附近的速度向量场U影响。除了表征粒子的几何形状的量(长度l和直径d)之外，示出了取向向量p。由相应箭头的方向和长度指示流速U的方向和大小。尽管速度向量的方向都一样，但是其长度不同，指示粒子附近的流速不是恒定的。然而，由于从底部到顶部向量的长度有恒定的增加，因此流速变化的局部量(即速度梯度)看起来相同。在粒子表面的每个点处，认为粒子精确地以局部流速(“无滑动(no-slip)”边界条件)移动。如果粒子周围的流速是恒定的，则这将导致简单的平移运动，即纯粹对流输运。否则，在存在速度梯度的情况下，如虚线箭头所指示的，粒子还进行旋转运动。总而言之，刚性粒子进行的最普通类型的运动是以粒子周围的“平均速度”进行的平移，结合围绕其质心的旋转，该旋转由描述粒子附近的流速与其平均值的局部偏差的速度梯度驱动。

已经由Jeffery将以上所给出的定性描述实现为数学模型[2]，其假定粘性力是主导的，惯性力因此可忽略，并且在单个纤维运动期间由其扫过的区域上的局部流的变化较小。如上所说明的，所有这些假设都在分散在热塑性塑料熔体中的短纤维的情况中得到了满足。所假设的纤维周围的流速变化较小意味着局部速度梯度张量足以精确地描述该变化。该张量是3×3矩阵，根据流速向量的分量U_j关于空间中的点r＝(x₁，x₂，x₃)^T的笛卡尔坐标x_i的偏导数(即${(▿\otimes U)}_{\mathrm{ij}}=\partial U_j/\partial x_i$)计算该矩阵的元素。因此如

果r₀是粒子质心的空间坐标，则通过一阶Taylor展开式$U(r,t)\approx U(r_0,t)+{(▿\otimes U(r_0,t))}^T·(r-r_0)$来具有可忽略的误差地较好地近似在粒子附近的点r处和时间t处的流速U的值，并且可以把速度梯度看作局部恒定量，即可以假设$▿\otimes U(r,t)\approx ▿\otimes U(r_0,t)$适用于r₀附近的所有点r。

在这种情况下，使用以下Jeffery方程计算由作为空间和时间坐标的函数的其取向向量p(r，t)给出的纤维的瞬时取向状态：

$\left(1\right),\frac{D}{\mathrm{Dt}}p={\hat{M}}^T·p-(p^T·\hat{M}·p)p,$

这里以紧凑欧拉形式写出该方程。方程(1)的左侧(l.h.s.)的对流(或者材料)导数$\frac{D}{\mathrm{Dt}}p=\{\frac{\partial }{\partial t}+U·▿\}p$描述了纯粹纤维取向FO、流的局部速度场U(r，t)中的输运，而右侧(r.h.s.)对由有效局部速度梯度张量所驱动的粒子旋转运动进行建模。在无限细长纤维(即纵横比r_a→∞)的理想特殊情况下，张量等同于速度梯度张量在纤维具有有限纵横比(0＜r_a＜∞)的一般情况下，由以下项定义张量

$\left(1a\right),\hat{M}=\frac{\lambda +1}{2}(▿\otimes U)+\frac{\text{λ-1}}{2}{(▿\otimes U)}^T$

该项包含纤维几何形状参数λ＝(r_a²-1)/(r_a²+1)该参数对在旋转对称椭圆粒子的情况下粒子几何形状如何影响纤维的旋转运动进行编码。

利用唯一定义的分解将速度梯度张量分解成由剪切速率和涡度张量给出的其对称和不对称部分，以该可选方式将Jeffery方程写为

$\left(1b\right),\hat{G}=1/2[(▿\otimes U)+{(▿\otimes U)}^T]$(“剪切速率”)，

$\hat{W}=1/2[(▿\otimes U)-{(▿\otimes U)}^T]$(“涡度”)。

剪切速率和涡度张量经由以下恒等式关联到有效速度梯度张量：

$\left(1c\right),\hat{M}=\hat{W}+\lambda \hat{G}\iff \hat{G}=1/2\lambda (\hat{M}+{\hat{M}}^T),\hat{W}=1/2(\hat{M}-{\hat{M}}^T).$

将这些恒等式代入(1)并且考虑$p^T·\hat{W}·p=0$产生以下的可选形式Jeffery方程：

$\left(1d\right),\frac{D}{\mathrm{Dt}}p=-\hat{W}·p+\lambda [\hat{G}·p-(p^T·\hat{G}·p)p].$

仍然可以通过使用恒等式$-\hat{W}·p={\hat{W}}^T·p=1/2\mathrm{rotU}\times p$来重写(1d)的右侧的第一项来获得Jeffery方程的另一形式，该恒等式将涡度张量关联到流速向量场的旋度rotU(见例如[3])。

在粒子的几何形状是旋转对称但不是椭圆的情况下，方程(1)的形式以及有效速度梯度的定义保持不变，但是必须改变几何形状参数λ的公式，使得用有效纵横比$r_a^{*}\approx r_a$取代如上定义的“几何形状”纵横比r_a＝l/d，需要根据实验确定该有效纵横比的合适值。在这种意义上，可以将几何形状参数λ看作材料参数。在非对称非球形粒子的情况下，需要用包含三个几何形状参数的三个Jeffery型方程的耦合系统取代方程(1)[3]。

2.2纤维取向的宏观分布

在宏观水平上执行针对短纤维增强塑料的注塑成型模拟的情况下，因为单个纤维的尺寸小，因此包含在单个计算单元(即计算区域的体积单元)中的纤维的数量非常大。这一点由图4所示的显微图图解，该显微图表示从由短纤维增强塑料制成的部件得到的典型样本。因此通过宏观模型(即FO向量p的分布函数Ψ(p，r，t)(见[5]))描述包含在计算单元中的材料的局部FO状态。我们注意到尽管假设纤维浓度是均匀的，但是分布函数ψ仍然经由局部流速场U(r，t)及其梯度参数地依赖于(t，t)坐标。作为附加点，需要将说明纤维的相互作用、影响它们的取向的项包括在宏观水平上的模型中。

经由以下相应的Fokker-Planck方程实现将作为单独的、无相互作用的纤维的微观模型的Jeffery方程转换成产生许多相互作用纤维的FO统计的宏观模型：

$\left(2\right),\frac{D}{\mathrm{Dt}}\Psi =-{▿}_p·[\Psi \frac{D}{\mathrm{Dt}}p-D_r{▿}_p\Psi ]$

其中以FO分布函数Ψ作为因变量。这里作为Jeffery方程(1)的右侧的简写，并且和是在三维欧几里得(Euclidian)空间的单位球体S²上定义的梯度和散度算子的符号。根据Folgar和Tucker的方法[4]，与局部速度梯度的有效剪切速率${\gamma }_{e_{\mathrm{ff}}}=\sqrt{2G_{\mathrm{ij}}{G}_{\mathrm{ij}}}$成比例的扩散系数D_r＝C_Iγ_eff的引入产生在浓缩悬浮液中纤维纤维相互作用的简单模型。这里G_ij是如(1b)中定义的剪切速率张量的分量，并且平方根号下的项(使用爱因斯坦求和约定)是其自收缩$G_{\mathrm{ij}}{G}_{\mathrm{ij}}={\Sigma }_{i,j}{G_{\mathrm{ij}}}^2=\hat{G}:\hat{G}$的简写。无量纲、非负相互作用系数C_I是悬浮液的材料参数。通常对于短纤维增强塑料，该材料系数具有在10^-3到10^-2范围内的小(正)值。我们注意到在不可压缩的(以及几乎不可压缩的)流中，“随机”扩散项(与表示“确定性”Jeffery动力学的项相比)的相对弱点可能是稳定性问题的来源。

2.3纤维取向张量和Folgar-Tucker方程

通过Fokker-Planck方程(2)计算局部FO分布需要在针对流模拟区域的每个计算单元的单位球体S²上定义的PDE的数值解，这是对于“工业尺度”3D问题来说负担不起的昂贵的任务。因此，Advani和Tucker[6]提出使用纤维取向张量，其被定义为分布函数的矩，并且因此用FO张量的一族矩方程取代Fokker-Planck方程。由于FO分布关于变量p的反转对称Ψ(-p，...)＝Ψ(p，...)(这反映±p的方向对应于相同的取向状态的事实)，因此所有奇数阶的矩一致地消为零，以使得Ψ的矩展开式仅包含偶数阶的元素。因此这个展开式的第一非平凡矩是第二个，其由下式给出

(以索引表示：符号表示在球体S²的表面上进行积分，其中dS(p)是积分度量。第二矩被称为二阶FO张量(或者FO矩阵)并且根据定义是实对称3×3矩阵。由于根据标准化FO分布(也根据定义)，因此当满足$p^2={\Sigma }_{i=1}^3{p}_i^2=1$时，FO矩阵具有明显的性质，即其迹$\mathrm{Tr}\left({\hat{a}}^{\left(2\right)}\right)={\Sigma }_k{a}_{\mathrm{kk}}^{\left(2\right)}$等于1。展开式体系中的下一个非平凡矩是四阶FO张量其定义为：

(以索引表示：张量是完全对称的并且额外地拥有各种标准化性质：由于p²＝1，因此相同索引的任意对的和总是产生的相应元素(例如，${\Sigma }_k{a}_{\mathrm{ijkk}}^{\left(4\right)}=a_{\mathrm{ij}}^{\left(2\right)}$)，并且两个相同索引对的和总是等于1(例如，${\Sigma }_{j,k}{a}_{\mathrm{jjkk}}^{\left(4\right)}=1$)，因此包含关于的完整信息。

可以通过交换微分和积分运算来获得针对每阶的矩的上述矩方程族，即

用Fokker-Planck方程(2)的右侧的项取代并且解析地估算相应的积分。如果将这个过程应用于FO矩阵，则获得所谓的Folgar-Tucker方程(FTE)作为分级组织的一系列矩方程中的第一非平凡方程：

$\left(5\right),\frac{D}{\mathrm{Dt}}{\hat{a}}^{\left(2\right)}={\hat{a}}^{\left(2\right)}·\hat{M}+{\hat{M}}^T·{\hat{a}}^{\left(2\right)}-(\hat{M}+{\hat{M}}^T):{\hat{a}}^{\left(4\right)}+2D_r[\hat{1}-3{\hat{a}}^{\left(2\right)}]$

如在Jeffery方程中，(5)的左侧的对流导数描述由欧拉参考系中纤维的平移运动造成的局部FO状态(由矩展开式的该阶的FO矩阵表示)的纯粹输运，而(5)的右侧的前两项对由局部有效速度梯度所驱动的它们的旋转运动进行建模。从Fokker-Planck方程(2)中的扩散项的存在产生(5)的右侧的第三项。在FTE的水平上，其产生阻尼效果，该效果将FO矩阵拖向各向同性状态${\hat{a}}^{\left(2\right)}=1/3\hat{1}.$

在其欧拉形式(5)中，FTE是对流反应(convection-reaction)类型的一阶PDE的耦合系统。从拉格朗日观点来看，(5)是针对的分量的ODE的耦合系统。

2.4FO矩阵的基本性质

作为实对称3×3矩阵，FO矩阵拥有3个实特征值μ_k，μ_k的相应特征向量为E_k，E_k形成的标准正交基，即：${\hat{a}}^{\left(2\right)}·E_k={\mu }_k{E}_k,$且E_k·E_l＝δ_kl，其中δ_kl是克罗内克(Kronecker)符号，对于k＝l，该符号等于1，否则等于零。显而易见的是，通过构造，FO矩阵的特征值μ_k是非负的并且(当满足p²＝1时)它们的和(其与的迹相同，见上)总是等于1(即：${\Sigma }_{k=1}^3{\mu }_k={\Sigma }_{k=1}^3{a}_{\mathrm{kk}}^{\left(3\right)}=1$)。

对于任意单位向量p₀，特殊FO矩阵${\hat{a}}_0=p_0\otimes p_0$表示所谓的单轴取向状态，其中100％的纤维在±p₀给出的方向上取向。因此，p₀的符号不参与定义这个特殊FO矩阵的双积乘积。向量p₀是的特征值为1的特征向量，并且其相应的特征向量位于与p₀正交的平面中的剩余的两个特征值为零，否则为任意。

可以其谱表示的方式${\hat{a}}^{\left(2\right)}={\Sigma }_{k=1}^3{\mu }_k{E}_k\otimes E_k$来写FO矩阵，该谱表示具有由FO矩阵的特征向量形成的单轴取向状态的加权和的形式，权重由特征值给出。这允许将特征值μ_k解释为沿着相应的特征向量E_k的方向取向的纤维的局部部分。在这个意义上，FO矩阵的谱数据{μ_k，E_k}_{k＝1，2，3}表示在悬浮液的小体积内的纤维的局部宏观取向状态。

这作为FO矩阵和FTE的相空间的以下(数学形式)定义的动机[12]：定义：当且仅当实对称3×3矩阵是半正定并且其迹等于1时，该矩阵是FO矩阵。FTE的相空间M_FT是所有FO矩阵的集合。

可以示出的是，集合M_FT等于使用适当地标准化的(但是否则为任意的)分布函数从类型(3)的矩积分得到的所有实对称3×3矩阵的集合。最近由作者之一利用相空间M_FT的拓扑和几何性质给出了其数学特征化(见[12])。由于FTE的数值积分结果的可解释性需要因变量在FO计算的所有步骤(即，也包括中间步骤)过程中具有如上所述的特殊谱性质，因此确切地知道矩阵是否属于集合M_FT非常具有实际重要性。

2.5闭合问题

所谓的闭合问题源自这样的事实，即在矩展开式的每一阶，矩的DE包含作为变量的下一个更高阶的矩。尽管可以通过简单代数恒等式(即对一对两个相等索引求和，见上)来根据表示但是反过来是不可能的，因此需要将视为未知。在FTE中，由(5)的右侧的的存在而显示出了闭合问题，这妨碍了系统的可解性，除非通过将通过闭合+近似表示为的函数来闭合该系统。将闭合近似应用于FTE意味着由FO矩阵的某个合适的(一般为非线性)张量函数来取代(5)的右侧的精确(但未知的)四阶FO张量

公知的示例是由以下公式定义的混合闭合[7]：

$\left(6\right),{\hat{A}}_{\mathrm{hyb}}^{\left(4\right)}\left[{\hat{a}}^{\left(2\right)}\right]:=f_s\left({\hat{a}}^{\left(2\right)}\right)·{\hat{A}}_{\mathrm{qu}}^{\left(4\right)}\left[{\hat{a}}^{\left(2\right)}\right]+(1-f_s\left({\hat{a}}^{\left(2\right)}\right))·{\hat{A}}_{\mathrm{lin}}^{\left(4\right)}\left[{\hat{a}}^{\left(2\right)}\right]$

尽管有一些众所周知的缺点，但是由于该公式的代数简单性以及数值鲁棒性因而其是可接受的选择[7]。将闭合定义为在两个较简单类型的闭合之间的(凸面)插值：二次闭合定义为：

$\left(6a\right),{\hat{A}}_{\mathrm{qu}}^{\left(4\right)}\left[{\hat{a}}^{\left(2\right)}\right]:={\hat{a}}^{\left(2\right)}\otimes {\hat{a}}^{\left(2\right)}$(即以索引表示为${\left({\hat{A}}_{\mathrm{qu}}^{\left(4\right)}\right[{\hat{a}}^{\left(2\right)}\left]\right)}_{\mathrm{ijkl}}:=a_{\mathrm{ij}}^{\left(2\right)}{a}_{\mathrm{kl}}^{\left(2\right)}$)

其在单轴取向分布的特殊情况中产生精确的结果，并且由下式给出线性闭合

${\left({\hat{A}}_{\mathrm{lin}}^{\left(4\right)}\right[{\hat{a}}^{\left(2\right)}\left]\right)}_{\mathrm{ijkl}}:=-\frac{1}{35}({\delta }_{\mathrm{ij}}{\delta }_{\mathrm{kl}}+{\delta }_{\mathrm{ik}}{\delta }_{\mathrm{jl}}+{\delta }_{\mathrm{il}}{\delta }_{\mathrm{jk}})$

(6b)

$+\frac{1}{7}({\delta }_{\mathrm{ij}}{a}_{\mathrm{kl}}^{\left(2\right)}+{\delta }_{\mathrm{kl}}{a}_{\mathrm{ij}}^{\left(2\right)}+{\delta }_{\mathrm{ik}}{a}_{\mathrm{jl}}^{\left(2\right)}+{\delta }_{\mathrm{jl}}{a}_{\mathrm{ik}}^{\left(2\right)}+{\delta }_{\mathrm{il}}{a}_{\mathrm{jk}}^{\left(2\right)}{+{\delta }_{\mathrm{jk}}a}_{\mathrm{il}}^{\left(2\right)})$

其在各向同性的情况下是精确的。由标量取向因子提供在这些极端情况之间的插值权重

$\left(6c\right),f_s\left({\hat{a}}^{\left(2\right)}\right):=1-27\det \left({\hat{a}}^{\left(2\right)}\right).$

由于行列式det是FO矩阵的不变量，因此对于标量取向因子同样成立。

我们通过首字母缩写FTE-hyb来表示结合混合闭合近似(6)的FTE(5)的特殊变量。该变量非常具有实际意义，本发明的优选实施例也是基于Folgar-Tucker模型的FTE-hyb变量。

2.6FTE作为微分代数系统

针对任意实对称矩阵而形式上良好地定义(5)的右侧。对于(5)的右侧的性质是否可以自动将解轨迹限制到域M_FT这个问题，可以至少针对其精确形式(即无闭合近似)下的FTE来肯定地回答：如果将精确四阶FO张量代入(5)的右侧中，并且然后作为(5)的解计算则结果与FO分布的精确二阶矩相同，其通过使用具有相同的参数和D_r的Fokker-Planck方程(2)的解Ψ估算动量积分(3)来获得。通过这个论据，可以推断，使用精确四阶FO张量可以将(5)的解写为矩积分(3)，并且因此必然地将其轨迹限制到域M_FT。

然而，如果应用了闭合近似，则这个论据不再有效，如果必须在对完全分布函数没有任何先前知识的情况下数值求解FTE，则应用闭合近似总是必要的。以此方式，将FTE的解限制到相空间M_FT的问题总是出现在所有的具有实际意义的问题中。

将FTE的解轨迹限制到域M_FT的必要性给因变量加入附加限制，可以依据其不变量的代数不等式来用公式表示该限制(见[12])。这些关于因变量的限制将FTE变成微分代数系统(DAS)并且需要在用于数值积分的程序中注意该限制。

3.FTE的相空间的数学特征化

可以根据FO矩阵的特殊性质直接推导出相空间M_FT的基本拓扑性质，其中FO矩阵的特殊性质概括如下

定理1：相空间M_FT是被限制到由迹条件$\mathrm{Tr}\left(\hat{a}\right)=1$定义的五维超平面的实对称3×3矩阵的向量空间的有界凸子集。

M_FT的凸度考虑到向M_FT的投影映射的定义，其可以与任何适当的ODE积分器组合以构成用于FTE的合适的积分方法(见[9]的第IV.4章)。可以通过对实对称3×3矩阵的以下特征多项式的分析获得FO矩阵的不变量代数特征化

$P_a\left(\mu \right)=\det (\mu \hat{1}-\hat{a})={\mu }^2(\mu -S_a)+K_a\mu -D_a$

系数$S_a=\mathrm{Tr}\left(\hat{a}\right),$$K_a=\frac{1}{2}[\mathrm{Tr}{\left(\hat{a}\right)}^2-\mathrm{Tr}(\hat{a}·\hat{a})]$和$D_a=\det \left(\hat{a}\right)$是矩阵的不变量。(在文献中，有时由I₁＝S_a、I₂＝K_a和I₃＝D_a表示这些不变量。)可以根据这些不变量用公式表示FO矩阵的代数特征化，其根据

定理2：当且仅当实对称3×3矩阵的迹S_a等于1并且其不变量K_a和D_a是非负时，该矩阵是FO矩阵。

图6通过对特征多项式P_a(μ)的二次和三次项以及线性和常数项的单独分析而图解了此表述：正迹S_a明显地提供正特征值μ_k的存在，而两个不变量K_a和D_a的非负值阻止了负特征值的存在。此外，由于矩阵总是具有三个实特征值，因此条件S_a＞0，K_a≥0和D_a≥0对于的所有特征值为非负不仅是必要的而且是充分的。结合迹条件S_a＝1，这完成了以上定理的证明。

可以通过分别查看FO矩阵的对角和非对角元素来获得相空间M_FT的几何形状的概念，该相空间根据定理1是5D对象。首先我们注意到由给出的对角元素总是非负的并且满足迹条件∑_ka_kk＝1。因此，不依赖于坐标系的选择，它们被限制到三角集合(见图7)

这提供了必要条件，该必要条件将FO矩阵关于它们的对角元素特征化并且产生相空间集合M_FT的“对角部分”的不变量表示。

可以通过为实对称3×3矩阵的对角和非对角元素的三元组引入符号(x，y，z)和(u，v，w)，来获得域M_FT的“非对角”部分的正规描述，取任意(但固定)的对角三元组(x，y，z)∈Δ_a并且形式上定义属于固定对角三元组的所有“容许的”非对角三元组的集合

如关于以上定理2所说明的，在代数上可以将集合N_(x，y，z)特征化为同时满足下面一对不等式的所有非对角三元组(u，v，w)的集合：

$\left(9\right),K_a\ge 0\iff u^2+v^2+w^2\le \mathrm{xy}+\mathrm{xy}+\mathrm{zx},$

$\left(10\right),D_a\ge 0\iff u^{2}x+v^{2}y+w^{2}z-2\mathrm{uvw}\le \mathrm{xyz}.$

在图8中示出了对应于“极大对称”情况x＝y＝z＝1/3的特殊示例，其它情况(见图9)对应于这个四面体形的体积的扭曲形式。

为了补充不变量K_a和D_a为非负值的限制(9)和(10)，我们注意到对于任何具有正迹S_a＞0的实对称3×3矩阵，也根据以上由$K_a\le 1/3S_a^2$和$D_a\le 1/27S_a^3$来限制这些不变量，以便对于FO矩阵(S_a＝1)来说，总是由0≤K_a≤1/3和0≤D_a≤1/27限制它们的值。

通过引入集合Δ_a和N_(x，y，z)而获得的(形式)纤维化$M_{\mathrm{FT}}={\cup }_{(x,y,z)\in {\Delta }_a}(x,y,z)\times N_{(x,y,z)}$有助于得到M_FT的总体结构的图画(见图10)：可以通过局部地将可容许的非对角三元组的集合N_(x，y，z)附连到其“基点”(x，y，z)∈Δ_a的过程使单个“纤维”可视化，并且该基点贯穿整个三角集合Δ_a的随后的变化覆盖整个相空间。

在本节概括的结果清楚地显示了相空间M_FT是复杂的数学对象。根据定理2中叙述的严密数学结果，将被自然地定义为对称3×3矩阵的实向量空间中的演化方程的FTE的任何解轨迹限制到相空间域必然地意味着结合单位迹条件的不变量不等式(9)和(10)的满足。这个事实不可避免地将FTE转变成DAS。与单位迹条件不同，在存在任何目前公知的闭合近似的情况下，一般不保持不变量不等式(9)和(10)的有效性，而在不存在针对四阶FO张量的任何闭合近似的情况下，单位迹条件的满足被“构建到”精确FTE中，并且在对于较大类的闭合近似的相当普遍的先决条件下，单位迹条件仍然有效(见第4节)。(同样地，在“精确”FTE的情况下，可以宁可通过如2.6节给出的间接推论而不是通过FTE本身的代数结构推导(9)和(10)的有效性。)在学术工作中经常示出的简单例子中通常忽略这些数学事实。

然而，任何数值求解FTE的过程有必要包含闭合近似作为其基本元素。因此，任何适合于处理更复杂的情形(如通常在工业应用中出现的情形)的模拟过程必须通过适当的控制过程来对此进行解释，该控制过程将解轨迹限制到其理论上可容许的域。因此对于严格的工业应用来说，将FTE作为微分代数系统的适当处理是强制性的。

在7.7节中，提出了提供适当类型的不变量控制的过程。已经在优选实施例中实施了该过程。

4.迹守恒和稳定性的问题

将闭合近似应用于FTE意味着以FO矩阵的某(一般为非线性)张量函数取代(5)的右侧的精确(但是未知)四阶FO张量(即：数学表示为${\hat{a}}^{\left(4\right)}\leftarrow {\hat{A}}^{\left(4\right)}\left[{\hat{a}}^{\left(2\right)}\right],$进一步的细节见以上关于闭合问题的章节)。取决于闭合的选择，四阶张量具有精确张量的较小或者较大量的性质。

一个被所有合理的闭合近似满足的关于的对称性质的特殊要求是所谓的正交各向异性对称性，其由这组等式给出

$\left(11\right),A_{\mathrm{ijkl}}^{\left(4\right)}=A_{\mathrm{jikl}}^{\left(4\right)},A_{\mathrm{ijkl}}^{\left(4\right)}=A_{\mathrm{ijlk}}^{\left(4\right)},A_{\mathrm{ijkl}}^{\left(4\right)}=A_{\mathrm{klij}}^{\left(4\right)},$

即，要求4阶张量关于第一和第二索引对(ij)和(kl)以及两对的互换$\left(\mathrm{ij}\right)\leftrightarrow \left(\mathrm{kl}\right)$对称。如果假设是由具有正交各向异性对称性质(11)的闭合近似${\hat{a}}^{\left(4\right)}\leftarrow {\hat{A}}^{\left(4\right)}\left[{\hat{a}}^{\left(2\right)}\right]$所增广的(5)的解，则我们可以形式地推出其迹的以下微分方程(缩写：DE)(见[17])：

$\left(12\right),\frac{D}{\mathrm{Dt}}\mathrm{Tr}\left({\hat{a}}^{\left(2\right)}\right)=\underset{i,j}{\Sigma }{M}_{\mathrm{ij}}(a_{\mathrm{ij}}^{\left(2\right)}-\underset{k}{\Sigma }{A}_{\mathrm{ijkk}}^{\left(4\right)})+6D_r[1-\mathrm{Tr}\left({\hat{a}}^{\left(2\right)}\right)]$

在(12)的推导中还不使用迹条件$\mathrm{Tr}\left({\hat{a}}^{\left(2\right)}\right)=1,$因为应该通过对该方程的分析研究FTE的解的迹的稳定性和守恒。如果将闭合近似应用于FTE，则后者还满足标准化条件

$\left(13\right),a_{\mathrm{ij}}^{\left(2\right)}=\underset{k}{\Sigma }{A}_{\mathrm{ijkk}}^{\left(4\right)},$

并且迹的DE(12)简化为

$\left(14\right),\frac{D}{\mathrm{Dt}}\mathrm{Tr}\left({\hat{a}}^{\left(2\right)}\right)=6D_r·[1-\mathrm{Tr}\left({\hat{a}}^{\left(2\right)}\right)].$

由于扩散参数D_r根据定义是非负的，因此如果闭合近似${\hat{a}}^{\left(4\right)}\leftarrow {\hat{A}}^{\left(4\right)}\left[{\hat{a}}^{\left(2\right)}\right]$满足条件(11)和(13)，则可以推出(在D_r＞0的情况下)由迹条件定义的超平面是稳定积分流形或者(在D_r＝0的情况下)迹是FTE的第一积分。(应该注意的是，精确四阶FO张量根据定义是完全对称的并且满足${\Sigma }_k{a}_{\mathrm{ijkk}}^{\left(4\right)}=a_{\mathrm{ij}}^{\left(2\right)},$并且精确FO矩阵根据定义是对称的并且满足迹条件。)

以上的考虑显示了取决于近似张量函数和精确张量共同具有的性质的量，将闭合近似应用于FTE如何影响迹条件(其为施加到FO矩阵的不变量的最简单的条件)的有效性。

在混合闭合(6)的特殊情况下必须认真考虑迹稳定性，其中混合闭合仅拥有对称性质(11)但是不能满足标准化条件(13)。在这种闭合近似的情况下，迹的DE采用以下形式(另见[17])

其具有前因子

该前因子偏离对于满足(13)的所有闭合有效的较简单的表达式6D_r，这反映了这样的事实，即混合闭合未能满足(13)。

在下文中示出了在某些条件下，在前因子中出现的附加项可能使该项变为负。这意味着尽管超平面$\mathrm{Tr}\left({\hat{a}}^{\left(2\right)}\right)=1$仍然是具有混合闭合的FTE的积分流形，但是其可能变得局部不稳定。

可以通过以下论证的四阶段过程来推出前因子变为负的可能性：

(i)在不可压缩流场(即$\mathrm{Tr}\left(\hat{G}\right)=\mathrm{div}U=0$)的特殊情况下，前因子采用简化形式

(ii)在细长纤维的情况下，可以总是假设λ＝1，并且如果由描述的局部取向状态是近似单轴(即其特征值μ_k之一接近1)，则由于的行列式接近0，因此标量取向因子$f_s=1-27\det \left({\hat{a}}^{\left(2\right)}\right)$的值也近似为1。以此方式，获得简化表达式其在特定的假设下接近的精确值。

(iii)通过代入D_r＝C_Iγ_eff并且将谱表示${\hat{a}}^{\left(2\right)}={\Sigma }_k{\mu }_k{E}_k\otimes E_k$代入收缩$\hat{G}:{\hat{a}}^{\left(2\right)}={\Sigma }_{\mathrm{ij}}{G}_{\mathrm{ij}}{a}_{\mathrm{ij}}^{\left(2\right)},$获得前因子的近似表达式的经修改的形式：

(iv)在最后一步中，对于有效剪切速率${\gamma }_{\mathrm{eff}}=\sqrt{2G_{\mathrm{ij}}{G}_{\mathrm{ij}}}$标准化剪切速率张量的分量。标准化的剪切速率张量$\widehat{G^{\prime }}={\gamma }_{\mathrm{eff}}^{-1}\hat{G}$具有与相同的特征向量、零迹并且因此(与一样)具有负特征值。由于根据构造的分量是1级的，因此张量必然具有至少一个具有1级的绝对值的负特征值。

在到目前为止针对前因子的近似表达式的推导中所假设的各种环境下，这个事实是的期望的估计的关键。使用标准化的剪切速率张量可以用以下形式写出该近似表达式

考虑到在实际的模拟中，相互作用系数C_I是小正数(通常为0＜C_I＜10^-2，另见1.2节)，可以看出，如果以具有主特征值μ_j≈1的FO矩阵估算，将变为负，该特征值具有相应的特征向量E_k，该特征向量沿着剪切速度张量的最大负特征值的特征方向指向，因为在这种情况下，可以估计$\underset{k}{\Sigma }{\mu }_k(E_k^T·\widehat{G^{\prime }}·E_k)\approx -{\mu }_j\approx -1$并且因此得到

可以以严密的方式用公式表示以上论证。以此方式，可以给出数学证明，即针对任何不可压缩流和相互作用系数的值0≤C_I＜λ，具有混合闭合的FTE的解的迹在相空间M_FT的某些区域中变得局部不稳定，因为前因子的表达式(15b)在这些区域中不可避免地变为负。

因此，很清楚有必要在FTE的数值积分过程中注意迹稳定性问题。

5.FTE的数值积分：一般方面

用于FTE(包括任何类型的闭合近似)的数值积分的适合的方法的选择取决于方程类型的数学分类以及与FTE的特定代数形式相关的方面。如在以下的段落中所详细说明的，闭合FTE的一般结构显示出其构成FO矩阵的作为因变量的元素a_ij⁽²⁾的对流反应(convection-reaction)型的双曲偏微分方程(PDE)的耦合系统。

从方程(5)中给出的FTE的形式开始，可以通过明确地代入(5)的左侧的实质导数(material derivative)$\frac{D}{\mathrm{Dt}}...=\{\frac{\partial }{\partial t}+U·▿\}...$和(5)的右侧的象征闭合近似的数学符号${\hat{a}}^{\left(4\right)}\leftarrow {\hat{A}}^{\left(4\right)}\left[{\hat{a}}^{\left(2\right)}\right]$来重写这个方程，其产生闭合FTE的等价公式：

$\left(16\right),\{\frac{\partial }{\partial t}+U·▿\}{\hat{a}}^{\left(2\right)}={\hat{a}}^{\left(2\right)}·\hat{M}+{\hat{M}}^T·{\hat{a}}^{\left(2\right)}-(\hat{M}+{\hat{M}}^T):{\hat{A}}^{\left(4\right)}\left[{\hat{a}}^{\left(2\right)}\right]+2D_r[\hat{1}-3{\hat{a}}^{\left(2\right)}]$

FTE的这种形式已经揭示了其数学结构的大部分：左侧由输运或者对流类型的简单一阶偏微分算子组成，其中局部流速U(r，t)控制作为方程(16)的因变量的FO矩阵的元素a_ij⁽²⁾的解耦、纯粹对流输运。此外，(16)的右侧的各种项的代数结构显示，右侧(正如FO矩阵本身，并且如由数学一致性所要求的)组成二阶对称张量函数，在以下方程(16)的简要(但是等价)形式中将该函数表示为

$\left(17a\right),\{\frac{\partial }{\partial t}+U·▿\}{\hat{a}}^{\left(2\right)}={\hat{F}}^{\left(2\right)}({\hat{a}}^{\left(2\right)},{\hat{A}}^{\left(4\right)}[{\hat{a}}^{\left(2\right)}],\hat{M},D_r)$

尽管FO矩阵的右侧的显式依赖性是线性的，但是由闭合近似的函数形式${\hat{a}}^{\left(4\right)}\leftarrow {\hat{A}}^{\left(4\right)}\left[{\hat{a}}^{\left(2\right)}\right]$所确定的其隐式依赖性的性质一般而言是非线性的。因此除非选择了线性闭合近似，否则需要把函数作为因变量的非线性函数考虑。此外，可以认识到的是，右侧函数不可避免地导致针对FO元素a_ij⁽²⁾的单独方程的相互耦合，除非有效速度梯度是对角矩阵并且闭合函数具有非常特殊的代数结构。

在通过方程(1a)解出了对实际速度梯度张量和纤维几何形状参数λ的依赖性并且以有效剪切速率的形式用定义表达式D_r＝C_Iγ_eff重新替代旋转扩散参数之后，可以用以下形式重写方程(17a)

$\left(17b\right),\{\frac{\partial }{\partial t}+U·▿\}{\hat{a}}^{\left(2\right)}={\hat{F}}^{\left(2\right)}({\hat{a}}^{\left(2\right)},{\hat{A}}^{\left(4\right)}[{\hat{a}}^{\left(2\right)}],▿\otimes U,{\gamma }_{\mathrm{eff}},C_I,\lambda ).$

闭合FTE的这种形式显式地揭示了右侧关于恒定模型参数λ和C_I的依赖性，但是仍然在对FO矩阵的显式和隐式依赖性之间进行区分，以及在对速度梯度张量的显式依赖性和由有效剪切速率公式${\gamma }_{\mathrm{eff}}=\sqrt{2G_{\mathrm{ij}}{G}_{\mathrm{ij}}}$引起的对的隐式依赖性之间进行区分(详见关于Fokker-Planck方程的章节)。在文献中，有时使用索引表示以分量形式写方程(17b)，即

$\{\frac{\partial }{\partial t}+U·▿\}{a}_{\mathrm{ij}}^{\left(2\right)}=F_{\mathrm{ij}}^{\left(2\right)}(a_{\mathrm{ij}}^{\left(2\right)},A_{\mathrm{ijkl}}^{\left(4\right)}[a_{\mathrm{ij}}^{\left(2\right)}],{(▿\otimes U)}_{\mathrm{ij}},{\gamma }_{\mathrm{eff}},C_I,\lambda ).$

经常使用精简符号找到该方程的更加简化的形式$\frac{D}{\mathrm{Dt}}{a}_{\mathrm{ij}}^{\left(2\right)}=F_{\mathrm{ij}}^{\left(2\right)}(a_{\mathrm{ij}}^{\left(2\right)},{(▿\otimes U)}_{\mathrm{ij}}),$该精简符号隐藏了所有的参数依赖性并且不在分量函数F_ij⁽²⁾(...)对FO矩阵或速度梯度的显式和隐式依赖性之间进行区分。从索引表示回到在方程(17a/b)中使用的完全张量符号，获得了这些方程的简化形式

$\left(17c\right),\{\frac{\partial }{\partial t}+U·▿\}{\hat{a}}^{\left(2\right)}={\hat{F}}^{\left(2\right)}({\hat{a}}^{\left(2\right)},▿\otimes U)$

6.FTE的数值积分：“算子分裂”

存在各种用于对流反应型的耦合PDE系统$\{\frac{\partial }{\partial t}+U·▿\}y=F\left(y\right)$的数值积分方法。一种已经被证明尤其在非线性右侧函数F(...)的情况下鲁棒并且灵活的方法从等价地将PDE重整为$\frac{\partial }{\partial t}y=-U·▿y+F\left(y\right)$开始，并且通过将由线性微分算子和函数F(...)确定的“算子”所给出的右侧的两项分开处理而进行。在数学文献中，这种方法被称为[23-37]“分步法(method of fractional steps)”、“分裂法(splitting method)”或者简单地“算子分裂(operator splitting)”(关于可替选方法，参见[28-31])。通过“算子分裂”对具有形式$\frac{\partial }{\partial t}y=-U·▿y+F\left(y\right)$的方程进行的数值积分利用两个单独方程$\frac{\partial }{\partial }y=-U·▿y$和$\frac{\partial }{\partial t}y=F\left(y\right)$的数值(或者在某些情况下甚至是解析的)解，通过将右侧的两个算子之一设为零而从完整方程获得所述两个单独方程。第一个方程等价于齐次对流方程$\{\frac{\partial }{\partial t}+U·▿\}y=0,$第二个方程组成常微分方程(ODE)的系统，其一般而言是耦合并且非线性的。在优选实施例中使用这种类型的方法。使用方程(17a/b)的符号，将分裂方法应用于闭合FTE导致以下偏微分方程：

$\left(18a\right),\{\frac{\partial }{\partial t}+U·▿\}{\hat{a}}^{\left(2\right)}=\hat{0}$(→在右侧具有零矩阵的“对流步骤”)

和

$\left(18b\right)=\frac{\partial }{\partial t}{\hat{a}}^{\left(2\right)}={\hat{F}}^{\left(2\right)}({\hat{a}}^{\left(2\right)},{\hat{A}}^{\left(4\right)}[{\hat{a}}^{\left(2\right)}],▿\otimes U,{\gamma }_{\mathrm{eff}},C_I,\lambda )$(→“旋转步骤”)

将其视为在“算子分裂”方法的构架内的独立的子问题。由方程(18a)建模的物理效果是由于根据流速的流体质量的纯粹对流输运在模具型腔的被填充域内的FO统计(如由FO矩阵的元素所编码的)的空间再分布，而方程(18b)完全忽略对流输运的影响并且单独考虑由于由局部速度梯度驱动的纤维的旋转运动以及纤维间的相互作用造成的FO分布的局部变化速率。在下文中，描述若干变型，其示出可以如何组合两个子问题的(数值)解来产生完整方程(16)(或者其以简要表示等价形式(17a/b))的近似解。

6.1“简单算子分裂”

使用最简单的方式(通常表示为“简单算子分裂”)来获得实际上根据两个偏微分方程求解的方程的近似解，以连续顺序求解第二个方程，将从第一个方程得到的中间解作为第二个方程的输入(即初始值)。为了详细说明，我们通过使用具有标识$y\leftrightarrow {\hat{a}}^{\left(2\right)}$和$F\left(y\right)\leftrightarrow {\hat{F}}^{\left(2\right)}({\hat{a}}^{\left(2\right)},...)$的模型方程$\frac{D}{\mathrm{Dt}}y=F\left(y\right)$来简化符号。

流求解器(即对悬浮液的流进行建模的软件)通过在位于包含在覆盖模具型腔、浇注系统和浇口的总体计算区域的填充部分Ω⁽ⁿ⁾内的空间中的离散点r_m附近的所有计算单元中、在离散时间瞬时t_n处求解针对质量、动量和能量的离散输运方程来产生流体的状态变量，其中包括流速U(r，t)和其梯度张量流求解器在时间t₀＝0处开始，并且使用步长Δt_n+1：＝t_n+1-t_n从时间t_n进行到t_n+1。该流计算意味着在时间t_n+1处的域Ω⁽ⁿ⁺¹⁾的计算取决于其之前的状态Ω⁽ⁿ⁾和新的速度场U⁽ⁿ⁺¹⁾。可以通过使用流体体积(VOF)方法来完成该计算。计算新的流前沿和新的域Ω⁽ⁿ⁺¹⁾产生在时间t_n+1处的Ω⁽ⁿ⁺¹⁾的单元中的所有状态变量的新值。在相应的步骤中，从在根据之前的计算步骤已知的“旧”域Ω⁽ⁿ⁾中定义的值$y_m^{\left(n\right)}=y(r_m,t_n)$开始，应该通过PDE$\{\frac{\partial }{\partial t}+U·▿\}y=F\left(y\right)$的数值解来计算位于新填充的域Ω⁽ⁿ⁺¹⁾的点r_m附近的所有计算单元中的时间t_n+1处的新值$y_m^{(n+1)}=y(r_m,t_{n+1}).$

这通过以下三步过程完成，该过程组成了“简单算子分裂”的一种可能的变型：

1.扩展步骤：从域Ω⁽ⁿ⁾的单元中的$y_m^{\left(n\right)}=y(r_m,t_n)$开始，计算在新域Ω⁽ⁿ⁺¹⁾的所有单元中的初始值y_m⁽ⁿ⁾。

2.对流步骤：从扩展步骤所提供的初始值y_m⁽ⁿ⁾开始，使用域Ω⁽ⁿ⁺¹⁾中的由流求解器所计算的流速U(r_m，t_n+1)，通过具有步长Δt_n+1的对流方程$\{\frac{\partial }{\partial t}+U·▿\}y=0$的数值解计算中间值

3.旋转步骤：从之前的对流步骤所提供的初始值开始，使用域Ω⁽ⁿ⁺¹⁾中的由流求解器提供的速度梯度通过具有步长Δt_n+1的ODE系统$\frac{\partial }{\partial }y=F\left(y\right)$的数值解计算最终值$y_m^{(n+1)}=y(r_m,t_{n+1}).$

图11图解了这种算子分裂变型。

上述的“简单分裂”变型产生完整方程$\{\frac{\partial }{\partial t}+U·▿\}y=F\left(y\right)$的近似解，其步长具有的一阶时间离散化误差(即O(Δt_n+1))，并且其用于根据优选实施例的方法中。将在随后的节中说明关于该过程的三个步骤的细节。

“简单算子分裂”的可选变型为：

(i)以使用时间t_n处的“旧”速度梯度的关于Ω⁽ⁿ⁾的旋转步骤开始，然后

(ii)进行“中间”步骤以将(i)的结果从Ω⁽ⁿ⁾扩展到Ω⁽ⁿ⁺¹⁾，和

(iii)以使用时间t_n+1处的流速U的对流步骤结束。

在这个变型中，以相反的顺序执行对流和旋转步骤。尽管理论时间离散化误差具有与上述第一种变型的情况相同的级，但是由于使用了不同瞬时时间的流速和速度梯度，因此这种可选变型趋向于一致性较差。可以通过另一变型来避免该问题

(i)以如上所述的已知值y_m⁽ⁿ⁾从Ω⁽ⁿ⁾到Ω⁽ⁿ⁺¹⁾的扩展步骤开始，然后

(ii)执行旋转步骤，随后

(iii)进行对流步骤

该变型对于步骤(ii)和(iii)使用Ω⁽ⁿ⁺¹⁾中、时间t_n+1处由流求解器提供的U和的“新”值。

6.2“对称算子分裂”

理论上来说，可以通过“对称算子分裂”方法来获得关于时间离散化误差的更高阶的精度。这种方法的基本想法是将具有完整步长的部分方程之一的一个中间步骤夹在具有半步长的其它方程的步骤之间。该方法的一种可能的变型导致以下的四步骤过程：

1.扩展步骤：从域Ω⁽ⁿ⁾的单元中的$y_m^{\left(n\right)}=y(r_m,t_n)$开始，计算在新域Ω⁽ⁿ⁺¹⁾的所有单元中的初始值y_m⁽ⁿ⁾。

2.预旋转步骤：从由扩展步骤提供的Ω⁽ⁿ⁺¹⁾的已知值的扩展y_m⁽ⁿ⁾开始，使用由流求解器提供的域Ω⁽ⁿ⁺¹⁾中的速度梯度通过具有半步长Δt_n+1/2的ODE系统$\frac{\partial }{\partial t}y=F\left(y\right)$的数值解来计算预旋转值y_m^(pre)。

3.对流步骤：从由预旋转步骤提供的作为初始值的y_m^(pre)开始，使用由流求解器计算的域Ω⁽ⁿ⁺¹⁾中的流速U(r_m，t_n+1)，通过具有全步长Δt_n+1的对流方程$\{\frac{\partial }{\partial t}+U·▿\}y=0$的数值解来计算中间值

4.后旋转步骤：从由前面的对流步骤提供的初始值开始，使用由流求解器提供的域Ω⁽ⁿ⁺¹⁾中的速度梯度通过具有半步长Δt_n+1/2的ODE系统$\frac{\partial }{\partial t}y=F\left(y\right)$的数值解来计算最终值$y_m^{(n+1)}=y(r_m,t_{n+1}).$

如上所述的“对称分裂”变型产生完整方程$\{\frac{\partial }{\partial t}+U·▿\}y=F\left(y\right)$的近似解，其具有步长的二阶时间离散化误差(即O(Δt_n+1²))。

另一种“对称分裂”变型在理论上可能将旋转步骤夹在具有半步长的两个对流步骤之间。

(i)以将已知值y_m⁽ⁿ⁾扩展到新域Ω⁽ⁿ⁺¹⁾开始，

(ii)使用流速U(r_m，t_n+1)在新域中进行半步长的“预对流”步骤，然后

(iii)使用速度梯度在Ω⁽ⁿ⁺¹⁾上执行完整旋转步骤，以及

(iv)以使用流速U(r_m，t_n+1)的另一个半步长“后对流”步骤结束，步骤(iv)最终导致也具有二阶精度的$y_m^{(n+1)}=y(r_m,t_{n+1})$的近似。如在“简单分裂”的情况一样，还存在第三种变型，其交换初始扩展步骤和预旋转步骤，以后者开始并且从而使用Ω⁽ⁿ⁾上的速度梯度的“旧”值。

由于在模具填充中遇到的典型流情况经常包括导致纤维的显著旋转的剪切流，因此与纯粹的对流输运的影响相比，“旋转算子”部分总是扮演重要(并且在大部分情形中为支配的)角色。因此通过针对第一对称分裂变型详细描述的预旋转和后旋转步骤的对流步骤的对称夹入组成了“对称分裂”方法的优选变型[19]。然而，在大部分实际情况下，由流体流计算提供的时间步长的大小足够小，并且6.1节所提出的“简单分裂”变型在这些情况下以足够的精度工作。

6.3可选分裂变型

基于右侧函数F(...)的代数结构，多种可选分裂方法是可能的。如果该函数由项的和F(...)＝∑_kF_k(...)组成，则可以将右侧分裂成由部分函数F_k(...)确定的子方程$\frac{\partial }{\partial t}y=F_k\left(y\right)$的另外的集合。有可能某些部分函数包括线性算子，即$F_k\left(y\right)={\hat{L}}_k·y.$在这种情况下，可以考虑部分对流方程(即$\{\frac{\partial }{\partial t}+U·▿\}y={\Sigma }_k{\hat{L}}_k·y$)的右侧的这些线性项，其由此保持线性，而将真正的非线性函数F_j(...)中的一些(或者全部)留在单个ODE$\frac{\partial }{\partial t}y={\Sigma }_j{F}_j\left(y\right)$或者如上所说明的单独ODE的列表$\frac{\partial }{\partial t}y=F_j\left(y\right)$内进行处理。

查看方程(16)的右侧，可以看到三项的和，第一和第三项为线性的，而包括闭合近似的中间项可能具有关于FO矩阵的非线性依赖性。(基于闭合的选择，可以将该中间项分裂成另外的子项和。)根据以上的评论，(16)的右侧函数的项结构开启了构建各种其它分裂方案的可能性。与基于分裂成方程(18a/b)的变型相比，这些其它替选方案不对应于可相互分离的物理效应(如“对流输运”对局部“旋转动力学”的情况中)的明确的区别。因此，它们仅仅可以被当作对物理FO现象的适当模拟而言几乎没有用的“人工数学分裂”。这里提及这些“人工”分裂变型的可能性仅仅是为了完整性并且在这里将不作进一步讨论。

6.4“对流步骤”的数值方法

本节说明用于对对流方程(18a)进行积分的数值方法，需要将其作为闭合FTE的“算子分裂”方法中的子问题进行求解。基于针对新的时间步长计算的速度U⁽ⁿ⁺¹⁾，使用一阶逆风方案[32]来组合输运方程系统矩阵，并且针对域Ω⁽ⁿ⁺¹⁾上的FO矩阵的每个分量求解所得到的方程系统。

6.5“扩展步骤”的数值方法

使用数值方法来将“旧”域Ω⁽ⁿ⁾的计算单元内的FO矩阵的已知值扩展到“新”域Ω⁽ⁿ⁺¹⁾的单元中。由经修改的VOF方法执行“新”域Ω⁽ⁿ⁺¹⁾的计算，该VOF方法作为应用在根据针对流方程(即质量、动量和能量守恒方程)的数值解的优选实施例的方法中的算法的一部分。对于两个域共有的那些单元，在扩展步骤期间FO矩阵的值保持不变。一般而言，域Ω⁽ⁿ⁺¹⁾包含若干“新填充的单元”，其需要初始分配FO矩阵的元素的合理值。根据优选实施例，根据之前计算的质量流，通过对应于计算单元内来自/去往其邻近单元的质量输运的加权平均来完成该分配。

由于假定粒子浓度是均匀的，并且在宏观水平上通过体积平均过程得到FO矩阵，因此假定单元将从其邻近单元得到FO贡献，该贡献与输运到其中的流体质量的净值成比例。根据之前总结的FTE的相空间(即所有可能的FO矩阵的集合)的数学特性，后者是被限制到(五维)$\mathrm{Tr}\left[a_{\mathrm{ij}}^{\left(2\right)}\right]=1$超平面的实对称3×3矩阵的六维向量空间的有界凸子集。由于根据质量输运加权的平均组成到FO矩阵的凸组合中，并且因此总是导致有效的FO矩阵，初始化过程与FT模型的相空间的拓扑结构相容(理论上要求的)，这是该过程的重要性质。

6.6“喷泉流”效应的处理

术语“喷泉流”表征在一大类粘性、非牛顿流体的情况下的前沿处的自由表面附近的整体宏观流型。在“喷泉流”中，流前沿附近的上行粒子从核心区域转移到壁边界。由于流体的材料性质，这种效果实际上“自动地”发生并且不需要任何附加的建模，即：至少在基于具有适当的非牛顿本构定律的Navier-Stokes方程的3D流计算中，“喷泉流”是突生现象。然而，如果流计算基于简化模型(例如，如在[1]中的Crochet的文章的11.3节中所讨论的Hele-Shaw类型)，则情况并非如此。在使用根据优选实施例的方法的通道流(channel flow)模拟中，可以在粒子迹中清楚地辨认“喷泉流”模式，这显示由流求解器说明了这种现象。

针对新填充的单元的纤维取向的初始化过程，需要采取某种特殊的操作以确保垂直于自由表面的FO分量是零(与关于纤维取向的“喷泉流”的所观察效果的一致性所要求的)。垂直FO分量的“被要求的”变为零不是(有意地)由上述的“扩展步骤”过程所强制的，其原因如下：从数学观点来看，FTE模型是由双曲输运方程的右侧耦合的双曲输送方程的系统(即具有非线性反应项的对流反应类型的系统，见上)。因此，不适合于(即数学上讲：不可能)将在自由表面垂直方向上的FO矩阵的分量归零规定为针对自由表面单元处的FO矩阵的边界条件。该FO分量的归零应该实际上是从所计算的流的“喷泉流”行为自动出现的。

然而，由于由时间和/或空间离散化误差以及用于计算在模具填充过程中的自由表面的运动所引入的数值不精确性，有可能需要增加一些“校正”以确保针对新填充的单元的初始化过程与“喷泉流”现象所要求的自由表面处的FO行为一致。在所有新填充的单元上，检查纤维取向张量的迹与理论值1相差不超过1％。如果差别太大，则通过正交投影到其特征空间中的特征向量的可容许的三角来校正纤维取向张量。

6.7模拟开始时的FO矩阵的初始化

(闭合)FTE的双曲结构要求在模具填充模拟的每个时间步、在浇口中(和附近)的计算单元中的FO矩阵的初始化作为边界条件。因此需要做出这种初始FO状态的合适选择。在根据优选实施例的方法中，出于以下所说明的目的选择各向同性FO状态${\hat{a}}^{\left(2\right)}=1/3\hat{1}$(对应于随机FO分布)。

在高粘性剪切流中，FO状态很大程度上受高剪切速率影响，因为它被快速地驱动到其局部准平衡的“最终”状态。因此在进模口(这是实际熔体进入该部分的地方)处所观察的FO状态由其贯穿浇注系统的流历史完全确定，并且很大程度上不依赖于在浇口处假定的其初始状态。另一方面，FTE的解析结构及其运动学(即在相空间中可能的状态)的分析显示，鉴于在浇口附近的浇注系统中的流场的随后的适应，浇口处的随机取向的假设是最佳选择，因为随机取向状态产生速度梯度的所有分量的全耦合。

上述算子分裂过程是参考图11描述的简化形式：

步骤1：设置新填充单元的初始条件和边界条件，

步骤2：根据流场移动纤维取向，和

步骤3：根据速度梯度计算纤维取向。

7.“旋转步骤”ODE系统的数值方法

一种数值方法用来对ODE的耦合系统(18b)进行积分，需要将该耦合系统作为闭合FTE的“算子分裂”方法的一部分进行求解。已经针对根据优选实施例的方法的FO模块的实施特别地开发了这一方法。它包含与短纤维增强热塑性塑料材料的3D注塑成型模拟的典型的工业应用情况下相关的各个方面，例如：

·使用FO矩阵的6个独立元素的完整集的FTE的公式化(相对于使用迹条件以消除对角元素之一并将FTE的因变量的数量减少一个的“标准过程”)；

·在FTE的右侧加入“惩罚”控制项，其将迹条件定义的超平面转换为FTE的稳定积分流形；

·利用与速度梯度分量相关的FTE的右侧函数的特定缩放行为以构造积分方法，该方法根据局部剪切速率的大小(速度梯度的最大分量)选择时间积分方案；

·包含被限定在[0，1]区间内的标量取向因子的“稳定混合闭合”的实施；

·使用最小数量的运算的用于具有混合闭合的FTE的特殊r.h.s.函数的估算的高效过程的实施。

下面详细讨论这些方面。

7.1混合闭合的稳定化

由方程(6)和方程(6a/b)定义的混合闭合近似尤其在以下情况下遇到稳定性问题：速度梯度具有复杂的结构，即不是简单的剪切/拉伸类型而是反映复杂的3D流型，并且时间步长(如流求解器所确定的)的大小变得相当大。标量取向因子向负值的偏离是这些不稳定性的首要来源：一旦该因子变为负，数值解很快变为不稳定，并且指数地偏离到远离FO矩阵可容许的相空间集合的值。

因此对的值的控制提高FO计算过程的数值稳定性。

由于FO矩阵的行列式总被限制在[0，1/27]的区间内，由方程(6c)定义的取向因子的理论上可容许的值也被$0\le f_s\left({\hat{a}}^{\left(2\right)}\right)\le 1$限制。因此，以如下方式修改混合闭合的标准定义(6)：

${\hat{A}}_{\mathrm{hyb}}^{\left(4\right)}\left[{\hat{a}}^{\left(2\right)}\right]:={\tilde{f}}_s\left({\hat{a}}^{\left(2\right)}\right)·{\hat{A}}_{\mathrm{qu}}^{\left(4\right)}\left[{\hat{a}}^{\left(2\right)}\right]+(1-{\tilde{f}}_s\left({\hat{a}}^{\left(2\right)}\right))·{\hat{A}}_{\mathrm{lin}}^{\left(4\right)}\left[{\hat{a}}^{\left(2\right)}\right]$

(19)

${\tilde{f}}_s\left({\hat{a}}^{\left(2\right)}\right):=\min (1,\max (1-27\det \left({\hat{a}}^{\left(2\right)}\right),0))\in \left[\mathrm{0,1}\right].$

由方程(6a/b)给出的线性闭合项和二次闭合项的定义保持不变，而由方程(6c)定义的取向因子被上面的方程(19)中定义的受限取向因子取代。的该定义意味着如果项被评估为区间[0，1]内的值，则${\tilde{f}}_s\left({\hat{a}}^{\left(2\right)}\right)=f_s\left({\hat{a}}^{\left(2\right)}\right),$否则取0或1作为最小值或最大值。

定义(19)被称为稳定化混合闭合近似。数值试验已经显示将的值限制在区间[0，1]内避免在所考虑的测试例子中产生严重的不稳定。

已经在各种3D注塑成型模拟中成功地测试了稳定混合闭合。

7.2经由迹条件减少变量的数量

作为方程(16)所示的闭合FTE的因变量的二阶FO张量为对称3×3矩阵。由于其元素$a_{\mathrm{ij}}^{\left(2\right)}=a_{\mathrm{ji}}^{\left(2\right)}$先验地组成一组六个相互独立的变量，因此FTE是六个PDE的耦合系统。虽然FTE的代数结构形式上允许任意对称3×3矩阵，但是该矩阵需要满足表达为对其不变量的限制的若干附加条件以使其满足作为FO矩阵的条件。

这些不变量条件中最简单的迹条件$\mathrm{Tr}\left(\hat{a}\right)={\Sigma }_{k=1}^3{a}_{\mathrm{kk}}^{\left(2\right)}=1$开启了明显的可能性以消除FO矩阵的对角元素a_kk⁽²⁾之一，并且因此将变量的数量减少一个。大多数(如果不是全部)已发表的研究理论流情况(例如简单剪切流或拉伸流)下的纤维悬浮液流的研究论文通过对于所要消除的特定对角元素作出固定的选择来使用此方法。

然后通过将例如$a_{33}^{\left(2\right)}=1-a_{11}^{\left(2\right)}-a_{22}^{\left(2\right)}$代入到方程(16)的右侧来实现消除，该方程以这种方式变得不依赖于a₃₃⁽²⁾。由于对于利用方程(14)或(15a)的所有合理的闭合近似FTE形式上与迹条件相容，因此迹条件和对角元素a₃₃⁽²⁾的DPE一致地满足(即作为形式上的代数恒等式)，并且只考虑FO矩阵的剩下的五个分量的集合以及进一步的它们的PDE的耦合系统就足够了。此过程是纯形式上的。对要以这种方式消除的对角元素的选择是否是好的取决于流的局部速度场U及其梯度的空间特性。对于仅显示流速度的平面变化(例如在x₁-x₂平面内)的简单类型的流，上面所说明的对a₃₃⁽²⁾的消除是合理的选择，因为在这种情况下纤维的主导旋转动力被限制在流平面内，因为在这种情况下速度梯度的特殊形式只导致对正交x₃方向的“弱耦合”。根据这一理由似乎选择消去a₁₁⁽²⁾或a₂₂⁽²⁾很可能会导致数值问题(例如关系到解的稳定性)并且因此导致较不满意的结果。然而，在复杂模具几何形状中的3D注塑成型模拟中情况则非常不同，因为流速U及其梯度在空间和时间中以复杂的方式变化，因此不能假设其具有任何特别的简单形式，所以固定的选择(如上面的选择)很可能不是最优的选择。因此通过迹条件消去FO矩阵的一个对角元素被认为对于工业悬浮液流的模拟而言是不适合的。

在Shampine的一系列论文[34-36]中提出了反对通过单位迹条件消去一个对角元素的与数值稳定性相关的另一个理由，在该系列论文中研究了具有特殊代数约束(“守恒定律”)的ODE系统。在[34-36]中考虑了ODE的耦合系统$\partial /\partial \mathrm{tc}=F\left(c\right),$其中通过质量守恒的条件${\Sigma }_{k=1}^n{c}_k=1$将“物类浓度”的n维向量c＝(c₁，...，c_n)^T限制在n-1维的超平面。而这个条件形式上产生了消去一个浓度的可能性，例如将关系$c_n=1-{\Sigma }_{k=1}^{n-1}{c}_k$代入右侧函数F(...)，并由此同时精确地满足守恒定律。Shampine指出在数值积分过程中使用这个“手法”仅仅导致浓度c₁，...，c_n-1的数值误差的积累，其中浓度c₁，...，c_n-1通过(减少的)ODE系统的数值积分计算得出，其中经由守恒定律代数地获得剩下的浓度c_n。虽然守恒定律由于构建而总是被精确地满足，但是无法保证以这种方式计算的数值解是准确的。实际上，已知尤其对于刚性ODE系统来说，如果这种“直接消去方法”使用的不够小心的话，数值解可能任意远地偏离真实解。根据[34-36]，这些论据在包括权重因子的向量w＝(w₁，...，w_n)^T的一般线性守恒定律$w·c={\Sigma }_{k=1}^n{w}_k{c}_k=1$的情况下同样保持有效并且导致相同的结论，即便在考虑更为一般的具有g(c)＝const.形式的非线性限制的情况下也是如此。

由于通过“单位迹”限制(其为线性“守恒定律”)增强的Folgar-Tucker方程是[34-36]中出现的数学方程式的具体例子，Shampine的数学分析显示，当在FO矩阵的对角元素中引入不可控的偏差时，通过单位迹条件消去一个对角元素很可能导致数值不稳定。

由于上面所给的理由，与通过单位迹条件消去一个对角元素的“标准方法”不同，根据优选实施例的方法使用FO矩阵的全部六个元素以及FTE系统的方程用于FO计算。

7.3利用惩罚项的动态迹稳定

如果闭合FTE系统及其“旋转步骤”部分被认为是由等式(16)和等式(18b)给出，并且不使用迹条件从独立变量集合中消去FO矩阵的一个对角单元，则需要某些附加方法来保持数值计算的FO矩阵的迹尽可能的接近其期望值1。

在论文[34-36]中建议避免任何类型的消去方法(如前面的7.2节所讨论的)并且通过将通常的ODE积分方法与某些考虑(或精确或近似地)满足所需要的守恒定律或代数限制的投影过程相结合来处理完整系统。对于每个时间步骤，首先在不考虑代数限制的情况下计算数值解，并且随后通过向限制方程(例如线性限制情况下的超平面)所定义的代数流形的最接近点进行投影映射来进行校正。

作为此方法的替选的原则试图通过向模型方程的右侧加入适当的控制项以作为限制，数值计算的MO矩阵的迹能够被保持为尽可能接近其期望值1。控制项的存在在没有任何附加措施的情况下产生限制或守恒定律的近似满足。在许多情况下可以显示出，“硬”控制有效地导致一种投影，而“软”控制允许与所要求的限制间有较小的偏离。

在根据优选实施例的方法中，控制方法尤其用于将数值计算的FO矩阵的迹近似地保持在其所需要的单位值。稳定化是基于通过分别向等式(16)或等式(18b)的右侧方程加入适合的控制项来将迹条件$\mathrm{Tr}\left(\hat{a}\right)=1$所定义的5D超平面变形到闭合FTE系统的稳定积分流形的思想。从不同的角度来看，这等价于$\mathrm{Tr}\left(\hat{a}\right)=1$成为迹的微分方程(12)的稳定固定点的情况。

已经确定任何适当的控制项应满足的一般要求如下：

·如果迹条件已经被满足，则控制项必须一致地变为零。

·如果数值解不能满足迹条件，则控制项必须迫使数值解回到$\mathrm{Tr}\left(\hat{a}\right)=1$超平面。

除了这两个主要要求之外，合理的控制项应该与和速度梯度的尺度改变(rescaling)相关的FTE方程的右侧的缩放行为相容(关于这一点的更多细节参见关于流受控时间积分的7.4节)，并且它应该包含允许“软”控制或“硬”控制的调节的调整参数，该“硬”控制实际上与“无限硬”控制的限制中的投影相似。

在混合闭合近似(标准和稳定化变型)的情况下，对适合的控制项的特定选择将适合于特定的闭合近似的选择。从FTE$D/\mathrm{Dt}{\hat{a}}^{\left(2\right)}={\hat{F}}^{\left(2\right)}(...)$开始，通过以下恒等式获得迹的DE

$D/\mathrm{DtTr}\left({\hat{a}}^{\left(2\right)}\right)=\mathrm{Tr}(D/\mathrm{Dt}{\hat{a}}^{\left(2\right)})=\mathrm{Tr}\left({\hat{F}}^{\left(2\right)}(...)\right)=...$

在(标准)混合闭合的特殊情况下，可以将右侧函数的迹估算为其产生迹的DE的简单形式(15a)，其具有由(15b)给出的可变前因子在稳定化混合闭合的情况下，标量取向因子被等式(19)所定义的其受限的(稳定化的)变型所代替，因此对于迹获得了形式上完全相同的DE，其中前因子包含而非

在(稳定化)混合闭合的特殊情况下对要加入右侧函数的惩罚项的具体选择为：

参数α需要为正，但不必是随时间恒定的。3×3矩阵需要为具有单位迹的对称矩阵，其它方面是任意的。矩阵定义控制项(20)的“方向”。它可以是随时间恒定或可变的。后者包含将定义为的函数的可能。对矩阵的不同选择对应于控制项的不同变型。根据优选实施例的方法中使用的变型使用$\hat{b}=1/3\hat{1}$(即b_ij＝1/3δ_ij)，其对应于与$\mathrm{Tr}\left(\hat{a}\right)=1$超平面正交的方向上的控制项。一个可选的变型由以下选择给出：$\hat{b}={\hat{a}}^{\left(2\right)}/\mathrm{Tr}\left({\hat{a}}^{\left(2\right)}\right)$(见[18])。在加入控制项(20)之后，必须也考虑具有稳定化混合闭合的闭合FTE的一般变型

$\left(21\right),\frac{D}{\mathrm{Dt}}{\hat{a}}^{\left(2\right)}={\hat{F}}_{\mathrm{hyb}}^{\left(2\right)}(...)+{\hat{C}}_{\mathrm{hyb}}(\alpha ,{\hat{a}}^{\left(2\right)},▿\otimes U)$

相似地，控制项进入ODE(18b)的右侧，其由此假设等式(21)的形式，其中实质导数D/Dt由偏导数代替。由于控制项(20)的迹估算为因此对控制项的具体形式进行修改以根据下式的修改(15a/b)

$\left(22\right),\frac{D}{\mathrm{Dt}}\mathrm{Tr}\left({\hat{a}}^{\left(2\right)}\right)=\mathrm{Tr}\left({\hat{F}}_{\mathrm{hyb}}^{\left(2\right)}(...)\right)+\mathrm{Tr}\left({\hat{C}}_{\mathrm{hyb}}(...)\right)=\alpha ·[1-\mathrm{Tr}{(\hat{a}}^{\left(2\right)})].$

当α＞0时，值$\mathrm{Tr}\left({\hat{a}}^{\left(2\right)}\right)=1$对应于DE(22)的稳定固定点，并且因此相应的超平面变为具有混合闭合的FTE的经修改的(“控制的”)形式(21)的稳定积分流形。作为结果，(21)的任意数值解的迹被控制项(20)保持在所要求的值1附近，这是因为该项迫使所有具有$\mathrm{Tr}\left({\hat{a}}^{\left(2\right)}\right)>1$的解在方向上趋向较低的迹值，并且相应地迫使具有$\mathrm{Tr}\left({\hat{a}}^{\left(2\right)}\right)<1$的解趋向较大的迹值(也在方向)。

控制项(20)的“强度”是可以通过调整参数α而适当的调节的：参数α的较小值导致相对“软”的令迹值向1的校正，较大的α值导致较强的推(接近于瞬间投影的效果)向$\mathrm{Tr}\left({\hat{a}}^{\left(2\right)}\right)=1$超平面。实际上，当被显式ODE积分方法使用时，较大的α值可能导致问题，因为在这种情况下较大的迹修正以交替的方式将数值解推向超平面的两侧，并因此导致伪振荡。已经显示[18]在相对较宽的中间值范围内选择α值在控制项不引入不希望的缺陷的情况下产生适当的数值结果。

对于导致闭合FTE的相应的右侧函数的闭合近似的任意具体选择，控制项的一般定义

$\left(20a\right),{\hat{C}}_{\mathrm{cl}}(\alpha ,{\hat{a}}^{\left(2\right)},▿\otimes U)=(\alpha -\frac{\mathrm{Tr}\left({\hat{F}}_{\mathrm{cl}}^{\left(2\right)}(...)\right)}{1-\mathrm{Tr}\left({\hat{a}}^{\left(2\right)}\right)})·[1-\mathrm{Tr}\left({\hat{a}}^{\left(2\right)}\right)]·\hat{b},\alpha >0,\mathrm{Tr}\left(\hat{b}\right)=1$

一致地导致迹DE(22)，并且因此导致对于(稳定化)混合闭合的特殊情况的如上所述的$\mathrm{Tr}\left({\hat{a}}^{\left(2\right)}\right)=1$超平面的稳定化。如果我们允许对称、单位迹3×3矩阵为的任意函数，则方程(20a)构成控制项的最一般的形式。

7.4流受控时间积分

使用“算子分裂”方法，“旋转步骤”ODE系统(18b)的数值解需要在具有由流求解器提供的全局步长Δt_n+1＝t_n+1-t_n的时间区间[t_n，t_n+1]中、在填充区域Ω⁽ⁿ⁺¹⁾的所有计算单元(由空间向量r_m标记)内计算。在每个单元中，局部速度梯度的分量作为一组外部系数进入(18b)的右侧。如果这些较大，则可以预期FO矩阵元的值在该时间区间内变化显著，而很小的速度梯度会导致FO矩阵元与由先前计算步骤提供的初始值相比可以忽略的变化。

在注塑成型模拟的所有实际情况中，在填充阶段，速度梯度在模具型腔内的填充区域显著变化。已经观察到，剪切速率的值在接近腔壁处通常较高，而在壁间的空隙的中间区域处具有相当低的值。这导致受到壁附近流速方向上较强剪切流的纤维较快地取向，而中心区域的纤维需要明显较长的时间来改变其取向状态，这导致众所周知的“三明治结构”(由壁附近的高取向区域间的具有相对较低纤维取向程度的“核”中心区域组成)，通常在由短纤维增强热塑料制成的零件中观察到该结构。

用于ODE系统积分的数值方案通过适当(优选地适应的)选择步长来考虑右侧函数的“强度”。试图以这种方式避免由于以下原因造成的不准确：在较大值区间上的粗略分步或者右侧函数的剧烈变化；以及在该区域上划分过小的步长，其中右侧函数的(绝对)值相当小。通过适当的选择步长以及积分方案，可以将ODE积分的数值误差控制在期望的界限。

在“旋转步骤”ODE系统的特别情况下，不管由速度梯度的变化造成的右侧方程如何改变，时间区间[t_n，t_n+1]上具有全局步长Δt_n+1和全局数值误差ε_FO的数值积分必须对于计算区域的所有单元一致。为了这一目的，已经设计出了特别的数值方案，其通过以下选择以对右侧函数的最少数量的估算实现这一目标：

·“积分规则”(欧拉、中点或四阶Runge Kutta)以及

·(无量纲)局部步长，以及

·对于所选规则的此步长的局部步数量

其适合于由局部速度梯度值确定的局部右侧函数的“强度”。本方案与用于数值ODE积分的“适应步长控制”的已知方法不同，因为它利用与速度梯度分量相关的方程(16)或(18b)的右侧函数的缩放行为并且其特别地适合于FTE系统。该方案固有的对于积分规则和局部时间步的数量和大小的选择过程是基于所使用的积分规则的理论离散化误差的试探方法。下面具体说明该方案的两个方面。

方程(17b)或(18b)的右侧函数的具体形式可以从方程(16)获得

${\hat{F}}^{\left(2\right)}(...)={\hat{a}}^{\left(2\right)}·\hat{M}+{\hat{M}}^T·{\hat{a}}^{\left(2\right)}-(\hat{M}+{\hat{M}}^T):{\hat{A}}^{\left(4\right)}\left[{\hat{a}}^{\left(2\right)}\right]+{2D}_r[\hat{1}-{\widehat{3a}}^{\left(2\right)}].$

为了清楚地解析对速度梯度的所有依赖关系，需要代入有效速度梯度和剪切速率张量的定义(1a/b)，即

$\hat{M}=\frac{\lambda +1}{2}(▿\otimes U)+\frac{\lambda -1}{2}{(▿\otimes U)}^T,\hat{G}=1/2[(▿\otimes U)+{(▿\otimes U)}^T]$

以及定义旋转扩散参数和有效剪切速率的恒等式$\hat{M}+{\hat{M}}^T=2\lambda \hat{G}$和公式$D_r=C_I{\gamma }_{\mathrm{eff}},{\gamma }_{\mathrm{eff}}=\sqrt{2\hat{G}:\hat{G}}.$对于如下考虑，抑制函数经由闭合函数对FO矩阵的间接依赖性以及经由有效剪切速率对速度梯度的间接依赖性以及对于恒定模型参数C_I和λ的依赖性，由此将右侧函数的代数表达重写(如同已经对方程(17c)的右侧所做的)为以下的简化形式

$\left(23\right),{\hat{F}}^{\left(2\right)}({\hat{a}}^{\left(2\right)},▿\otimes U)={\hat{a}}^{\left(2\right)}·\hat{M}+{\hat{M}}^T·{\hat{a}}^{\left(2\right)}-2\lambda \hat{G}:{\hat{A}}^{\left(4\right)}\left[{\hat{a}}^{\left(2\right)}\right]+2C_I{\gamma }_{\mathrm{eff}}[\hat{1}-3{\hat{a}}^{\left(2\right)}],$

在的变元表中只留下对FO矩阵和速度梯度的依赖性，这样做将是有益的。使用右侧函数(23)的ODE系统(18b)的等效变形(对应于(17c))可以写成以下形式

$\left(18c\right),\frac{\partial }{\partial t}{\hat{a}}^{\left(2\right)}={\hat{F}}^{\left(2\right)}({\hat{a}}^{\left(2\right)},▿\otimes U).$

用某个因子σ＞0与的分量的相乘导致相应的相乘由于(23)右侧的特殊的代数结构，导致代数恒等式

$\left(24\right),{\hat{F}}^{\left(2\right)}({\hat{a}}^{\left(2\right)},\sigma ▿\otimes U)=\sigma {\hat{F}}^{\left(2\right)}({\hat{a}}^{\left(2\right)},▿\otimes U),$

表示函数与其第二变元相关的缩放行为。

如前所述，一般任务是ODE的耦合系统(18c)在时间区间[t_n，t_n+1]内的数值积分，其具有对于填充单元区域内的所有单元一致的全局数值误差ε_FO。这可以通过利用右侧函数的特殊性质(24)通过以下方式对速度梯度、右侧函数和ODE系统(18c)进行缩放来实现：

·以最大绝对值确定速度梯度分量的值，即${\gamma }_{\max }:={\left|\right|▿\otimes U\left|\right|}_{\infty }=\underset{i,j=\mathrm{1,2,3}}{\max }\left|{(▿\otimes U)}_{\mathrm{ij}}\right|.$(应注意由于，参数γ_max还依赖于：空间向量r_m，其标记局部计算单元；以及时间坐标t_n或t_n+1，速度梯度在该时间坐标进入右侧函数。)

·引入缩放速度梯度$\overline{▿\otimes U}:{=\gamma }_{\max }^{-1}▿\otimes U$以及缩放时间坐标τ：＝γ_maxt(注意这两个均为无量纲量，并且此缩放仅局部地应用于单元r_m以及时间区间[t_n，t_n+1]内)，并且使用经缩放的量和等式(24)将ODE(18c)转换为缩放形式：

$\left(25\right),\frac{\partial }{\partial t}{\hat{a}}^{\left(2\right)}={\gamma }_{\max }^{-1}\frac{\partial }{\partial t}{\hat{a}}^{\left(2\right)}={\hat{F}}^{\left(2\right)}({\hat{a}}^{\left(2\right)},{\gamma }_{\max }^{-1}▿\otimes U)={\hat{F}}^{\left(2\right)}({\hat{a}}^{\left(2\right)},\overline{▿\otimes U})$

原始ODE(18c)需要在“全局”时间区间[t_n，t_n+1]上以“物理”步长Δt_n+1＝t_n+1-t_n(例如以[s]为单位测量)积分，这对于区域Ω⁽ⁿ⁺¹⁾的所有计算单元是相同的，并且右侧函数的“强度”根据速度梯度的变化而变化，如针对区域中的不同单元的参数γ_max的不同值所指示的。与此不同，经转换的ODE(25)需要在具有无量纲大小的“局部”缩放的时间间隔[τ_n，τ_n+1]上积分

$\left(26\right),\Delta {\tau }_{n+1}^{\left(m\right)}:={\gamma }_{\max }(r_m,t_{n+1})·\Delta t_{n+1},$

但是由于缩放的速度梯度进入(25)的右侧，由于下列事实，右侧函数具有近似一致的“单位强度”：

i.缩放的速度梯度由于构建而为无量纲量，其分量的绝对值等于或小于1。

ii.FO矩阵元以及由闭合函数提供的四阶张量元素的值也被限制在1。

iii.因此当用估算时，右侧函数的所有分量都为O(1)阶，而当用估算时的最大分量具有与γ_max相当的绝对值。

使用任意矩阵范数‖...‖测量的“强度”，可以将这些陈述给以准确的数学形式：

$\left(27\right),\left|\right|{\hat{F}}^{\left(2\right)}({\hat{a}}^{\left(2\right)},▿\otimes U)\left|\right|=O\left({\gamma }_{\max }\right),\left|\right|{\hat{F}}^{\left(2\right)}({\hat{a}}^{\left(2\right)},\overline{▿\otimes U})\left|\right|=O\left(1\right)$

由于根据(27)的变形的ODE(25)的缩放的右侧函数对于区域Ω⁽ⁿ⁺¹⁾内的所有单元具有一致的“单位强度”，通过如下选择可以在局部变化大小Δτ_n+1^(m)的间隔上实现具有一致的全局误差ε_FO的数值积分

i.具有足够小的离散化误差的“积分规则”，以及

ii.适合的子步数量以覆盖经缩放的时间间隔。

根据优选实施例的方法中使用的积分方案使用一组属于一步方法类(见[38]的16.1节和[39]的7.2.1-7.2.3)的简单积分规则，通过引用将其包含于此。根据优选实施例的方法的FO模块中使用的具体积分规则为：

·简单正向欧拉规则，其为一阶准确度方法并且每步需要右侧函数的1次估算，

·“中点”或二阶Runge-Kutta(RK2)规则，其为二阶准确度方法并且每步需要2次估算右侧函数，以及

·四阶Runge-Kutta(RK4)规则，其为四阶准确度方法并且每步需要4次估算右侧函数。

这些积分规则的依赖于其阶次p的如下特性被用于构建本发明的一个实施例中所用的特殊方案：

a)p阶的方法每子步需要p次估算右侧函数。

b)如果以N个大小为h＝Δτ/N的等距子步在总大小为Δτ的间隔上对ODE(系统)进行数值积分，则在该间隔的最终子步处积累的总误差ε_tot可以被估计为：ε_tot～Δτ·h^p。

c)如果要求该总误差小于预先选择的阈值ε_max，则子步的大小由h＜(ε_max/Δτ)^1/p限定。相似地，子步的数量需要大于N＞Δτ·(Δτ/ε_max)^1/p。

d)在总间隔上具有小于ε_max的总误差的积分可以通过用p阶方法进行N个子步来执行，其中间隔大小限制为Δτ＜(ε_max·N^p)^1/(p+1)。这意味着如果满足$\mathrm{\Delta \tau }<{ϵ}_{\max }^{1/(p+1)}$则大小为h＝Δτ的单子步(即N＝1)即足够。

考虑到这些，可以构建“混合”积分方案，该方案通过以下选择在(26)中定义的大小为Δτ_n+1^(m)的缩放的时间间隔上以小于ε_max＝ε_FO的全局误差产生方程(25)的数值积分

(i)具有三种积分规则之一的单步

(ii)使用RK4规则的适合大小的多个子步

其中该选择取决于相比于所要求的误差阈值ε_FO(见下面内容)的Δτ_n+1^(m)的相对大小。由于此误差限制对于区域Ω⁽ⁿ⁺¹⁾的所有计算单元是相同的，因此该积分以不依赖于局部值Δτ_n+1^(m)(或相似的${\gamma }_{\max }(r_m,t_{n+1})=\Delta {\tau }_{n+1}^{\left(m\right)}/\Delta t_{n+1}$)的一致的误差ε_FO执行。

所提出的方法特定地选择积分规则，该积分规则以根据上面的a)到d)方面所显示的关系所估计的右侧函数估算的最少次数来达到所要求的误差限定。虽然由流确定的时间步长Δt_n+1与总的注塑成型时间(其为秒量级)相比通常相当小(Δt_n+1≈10^-2...10^-4s)，无量纲量Δτ_n+1^(m)可以为O(1)阶，由于γ_max的值可以因为小间隔尺寸以及聚合体熔体的高粘度而变得大得多(10...100s^-1量级)。由于FO矩阵的更新值也为O(1)阶，所要求的误差界限的合理选择在${ϵ}_{\mathrm{FO}}\cong {10}^{-2}...{10}^{-4}$范围内改变，因此可以假设典型应用中的子步长h≤0.1。可以显示对于步长h＜1/2，RK4规则的步长h的单步比使用“中点”(RK2)规则的半步长h/2的两步或四分之一步长h/4的四个欧拉步要更加准确，虽然对于所有变型来说数值工作量(四次估算右侧函数)是相同的。

对积分规则和(子)步数量的选择基于相对于针对值p∈{1，2，4}计算的一些列误差界限${ϵ}_p=\sqrt[p+1]{{ϵ}_{\mathrm{FO}}}$的Δτ_n+1^(m)的大小，该系列误差界限定义各种积分规则的精度次序。由于0＜ε_FO＜1(实际上通常满足ε_FO＜＜1)，误差界限的大小随着次序参数p的增大而单调增长(总满足0＜ε_p＜ε_p+1＜1)，所以与下面给出的算法的公式化相关的误差界限总是根据0＜ε_FO＜ε₁＜ε₂＜ε₄＜1排序。(对于例如ε_FO＝0.001的典型的选择，得到数值ε₁≈0.032，ε₂≈0.1以及ε₄≈0.25。)相似地，对于无量纲缩放时间间隔Δτ_n+1^(m)的大小有最小阈值ε_min，在其之下即便步长为$h=\Delta {\tau }_{n+1}^{\left(m\right)}$的单欧拉步仅产生FO矩阵元对于先前值的小到可以忽略的更新。这通常发生在速度梯度非常小的区域的所有单元中，例如在两个相邻壁之间相对较大的间隙尺寸的情况下的中心核区域中。例如ε_min＝10^-6的最小阈值对于典型模具填充应用来说是合理的选择。

考虑以上讨论的所有方面，在区域的每个单元内执行的、并且在图12的流程图中图解的算法可以通过如下方式进行表达：

1.首先估算${\gamma }_{\max }:={\left|\right|▿\otimes U\left|\right|}_{\infty }$并且计算缩放的速度梯度$\overline{▿\otimes U}:={\gamma }_{\max }^{-1}▿\otimes U$和缩放的间隔大小$\Delta {\tau }_{n+1}^{\left(m\right)}:={\gamma }_{\max }(r_m,t_{n+1})·\Delta t_{n+1}.$

2.在步骤12.1中检查是否$\Delta {\tau }_{n+1}^{\left(m\right)}\le {ϵ}_{\min }.$

a.如果是，则跳过ODE积分(即什么也不做)，保持该单元的FO矩阵的先前值作为新(更新的)值并退出该过程。

b.如果否，则转至步骤12.2。

3.检查是否${ϵ}_{\min }<\Delta {\tau }_{n+1}^{\left(m\right)}\le {ϵ}_1.$

a.如果是，则通过步长为$h=\Delta {\tau }_{n+1}^{\left(m\right)}$的单欧拉步骤使用估算右侧函数来更新该单元中的FO矩阵的先前值，并退出该过程。

b.如果否，则转至步骤12.3。

4.检查是否${ϵ}_1<\Delta {\tau }_{n+1}^{\left(m\right)}\le {ϵ}_2.$

a.如果是，则通过步长为$h=\Delta {\tau }_{n+1}^{\left(m\right)}$的单“中点”(KR2)步骤使用更新该单元内的FO矩阵的先前值，并退出该过程。

b.如果否，则转至步骤12.4。

5.检查是否${ϵ}_2<\Delta {\tau }_{n+1}^{\left(m\right)}\le {ϵ}_4.$

a.如果是，则通过步长为$h=\Delta {\tau }_{n+1}^{\left(m\right)}$的单“四阶Runge-Kutta”(RK4)步骤使用更新该单元内的FO矩阵的先前值，并退出该过程。

b.如果否，则转至步骤12.5。

6.在${ϵ}_4<\Delta {\tau }_{n+1}^{\left(m\right)}$的一般情况下，通过执行具有适合步长的若干RK4步骤来执行FO矩阵的更新。

a.达到所需要的ε_FO精度的最少步骤数量为满足不等式N＞Δτ·(Δτ/ε_FO)^1/p的最小整数N。使用函数int(...)，返回浮点数值的整数部分，以及函数max(...)，返回两个数值中较大的一个，可以通过以下公式计算整数N≥1：

N＝max(int(Δτ·(Δτ/ε_FO)^1/p)+1，1).

b.然后利用$h=\Delta {\tau }_{n+1}^{\left(m\right)}/N$计算需要的步长。

c.一旦计算出N和h，通过RK4规则的具有步长h的N步(在右侧使用)更新该单元中的FO矩阵的先前值并退出该过程。

在一个优选实施例中使用上述算法。它为“旋转步骤”ODE系统(18b)在区域Ω_n+1中的所有单元在时间区间[t_n，t_n+1]上定义“混合”积分规则，其具有小于ε_FO的一致全局误差和对(18b)的右侧函数的最少估算次数。通过在右侧函数的估算中使用缩放的速度梯度以及基于缩放的间隔长度Δτ_n+1^(m)的大小计算步长h来“就地”完成该积分规则对(18b)的缩放形式(25)的应用。

在此算法的典型应用(具有误差界限参数ε_FO＝0.001以及ε_min＝10^-6)中，观察到在区域Ω_n+1中的大多数单元(即大约80％)中，通过单欧拉步骤执行FO矩阵的“旋转步骤”更新，较少数量的单元仅仅由于$\Delta {\tau }_{n+1}^{\left(m\right)}\le {ϵ}_{\min }$(即γ_max的较小值)而被“跨过”，并且有合理数量的据推测位于接近模具型腔壁并因此具有高剪切速率值的单元由一个或若干(通常不多于20个)“四阶Runga-Kutta”步骤更新。

7.5具有动态迹稳定化的“旋转步骤”ODE积分

如果右侧函数包含控制项(混合闭合情况下的特殊类型(20)或者如方程(20a)给出的一般类型)以实现动态迹稳定化，则右侧函数的缩放行为仍然与本节描述的“混合”积分算法兼容，只要将方程(20)或(20a)的控制参数定义为α＝α₀·γ_max，其具有典型大小的无量纲控制参数α₀～O(1)。

通过查看一般形式(20a)的控制项来对此进行说明。使用包含任意闭合近似的右侧函数的缩放特性(24)，我们可以根据等式的顺序重写它的迹

$\mathrm{Tr}\left({\hat{F}}_{\mathrm{cl}}^{\left(2\right)}({\hat{a}}^{\left(2\right)},▿\otimes U)\right)=\mathrm{Tr}\left({\gamma }_{\max }{\hat{F}}_{\mathrm{cl}}^{\left(2\right)}({\hat{a}}^{\left(2\right)},{\gamma }_{\max }^{-1}▿\otimes U)\right)={\gamma }_{\max }\mathrm{Tr}\left({\hat{F}}_{\mathrm{cl}}^{\left(2\right)}({\hat{a}}^{\left(2\right)},\overline{▿\otimes U})\right)$

并且，使用控制因子α＝α₀·γ_max，得到控制项(20a)的等效表达式

${\hat{C}}_{\mathrm{cl}}(\alpha ,{\hat{a}}^{\left(2\right)},▿\otimes U)={\gamma }_{\max }·({\alpha }_0-\frac{\mathrm{Tr}\left({\hat{F}}_{\mathrm{cl}}^{\left(2\right)}({\hat{a}}^{\left(2\right)},\overline{▿\otimes U})\right)}{1-\mathrm{Tr}\left({\hat{a}}^{\left(2\right)}\right)})·[1-\mathrm{Tr}\left({\hat{a}}^{\left(2\right)}\right)]·\hat{b}$

右侧的表达式现在包含作为单独的因子的γ_max以及依赖于的项，该项由缩放的速度梯度代替进行估算，但是其它方面形式上与(20a)一致。这可以由以下等式表示

$\left(28\right),{\hat{C}}_{\mathrm{cl}}(\alpha ,{\hat{a}}^{\left(2\right)},▿\otimes U)={\gamma }_{\max }·{\hat{C}}_{\mathrm{cl}}({\alpha }_0,{\hat{a}}^{\left(2\right)},\overline{▿\otimes U}),$

该等式还示出，如果将控制参数选择为α＝α₀·γ_max，则在其一般形式(20a)下的控制项具有与利用γ_max对速度梯度的缩放相关的“未经控制的”右侧函数完全相同的缩放行为。

在(稳定化的)混合闭合的特殊情况下，可以通过将项(20)重写为以下形式来以同样的方式从前因子提取γ_max

其中通过与相同的公式使用代替计算“再缩放”的前因子这明显地显示(28)对于在混合闭合情况下假设的控制项的特殊形式也是有效的。

以上的考虑显示，针对“旋转步骤”ODE系统的“混合”积分算法在不被改变的情况下也能够应用于有控制项的情况，只要α₀用作控制参数并且用缩放的速度梯度估算依赖于右侧函数的项。简单地通过向方程(25)右侧加入缩放的控制项来得到对应于方程(21)的缩放形式的ODE，即需要考虑用ODE

$\left(29\right),\frac{\partial }{\partial \tau }{\hat{a}}^{\left(2\right)}={\hat{F}}^{\left(2\right)}({\hat{a}}^{\left(2\right)},\overline{▿\otimes U})+{\hat{C}}_{\mathrm{hyb}}({\alpha }_0,{\hat{a}}^{\left(2\right)},\overline{▿\otimes U})$

代替(25)作为积分算法的分离时间步骤的基础。

在根据优选实施例的方法的FO模块中实施的时间积分方案通过将本节中定义的“混合”积分算法应用到包括利用控制项(20)进行的动态迹稳定化的ODE系统(29)来执行“旋转步骤”ODE系统的数值积分。

7.6具有稳定化混合闭合的FTE的高效估算

右侧函数的估算是FO计算过程中消耗最多的部分，因此本方法以最高效的方式进行估算。

通过7.4节提出的“流受控时间积分”(FCTI)的算法方法已经部分地解决效率方面。使用该方法，能够以对FTE右侧函数(可能包含控制项)的最少估算达到一致精度。

直接影响FO计算的计算消耗的最重要的因素是对闭合近似的选择。稳定化混合闭合以低计算消耗产生合理的精度，并且被本方法采用。

具有混合闭合的FTE的代数变形

提高效率的第一步是将混合闭合的代数定义(6)代入FTE的右侧函数(23)，(在几个代数变形之后)产生特定代数形式的右侧函数：

(30)的右侧的前两项的和在形式上与表达式一致，但是在1(a)给出的有效速度梯度矩阵处包含矩阵

$\left(31\right),\hat{N}:=\hat{M}-2/7(1-f_s){\hat{G}}_M,{\hat{G}}_M:=\hat{M}+{\hat{M}}^T=2\lambda \hat{G}$

右侧剩下的三项都具有标量前因子与矩阵量的乘积的形式。该前因子由一组公式给出

以上给出的针对前因子的公式与(15b)中给出的公式一致，因为${\hat{G}}_M=2\lambda \hat{G}$并且$\mathrm{Tr}\left({\hat{G}}_M\right)=2\mathrm{\lambda Tr}\left(\hat{G}\right),$其中$\mathrm{Tr}\left(\hat{G}\right)=\mathrm{devU}.$

公式(32)在可压缩或不可压缩的流速场的情况下有效。这里强调，虽然理论上当假设流为不可压缩时应满足divU＝0，这不被明确地使用，而是将所有包含$\mathrm{Tr}\left(\hat{G}\right)=\mathrm{div}U$的项保持在公式(32)中。如果由于数值误差(其不可避免地发生在流计算过程中)导致$\mathrm{Tr}\left({\hat{G}}_M\right)\ne 0,$发现(32)将导致稳定行为，而当在前因子公式(32)中将与成比例的项先验地假设为零的情况下则会发生不稳定。

形式(30)下的具有量(31)和(32)的右侧函数的使用是通向高效估算的基础步骤，因为右侧函数(30)的代数结构被以这样的方式组织使得各种计算只需要进行一次(当计算前因子时)，这节省了很多运算。

若干代数运算发生在频繁出现的共同的子表达式(例如(31)和(32)中的项2D_r或(1-f_s)/7)中，只对它们进行一次计算并将它们的值存储在辅助变量中用于后面的使用，这显著地减少了计算工作量。通过接下来的步骤，即根据

重新定义矩阵并且使用矩阵代替进行右侧方程求解，省去了另外的一些运算。利用代数恒等式

省略了项的详尽计算，并且将右侧函数根据下式重新定义

利用方程(31a)计算的附加运算为：(i)1个乘法以得到以及(ii)3个“加法”以从矩阵的对角元减去所节省的运算数量等于FO矩阵的所有分量与前因子的相乘以及的减法运算。因此，如果使用方程(30a)和(31a)代替(30)和(31)用于估算右侧函数，则运算的总数量更少。

在稳定化混合闭合的情况下，根据(19)计算受限变量取向因子因此，用代替(32)中的f_s在不影响计算消耗的情况下得到“稳定化前因子”和

控制项(见等式(20))

的加入导致对于右侧估算的另外的运算需要。此控制项的一般结构由具有前因子的给出，并且因此具有如上面讨论的相同的形式如果为投影矩阵选择$\hat{b}=1/3\hat{1}$或$\hat{b}={\hat{a}}^{\left(2\right)}/\mathrm{Tr}\left({\hat{a}}^{\left(2\right)}\right),$则前因子被并入前因子或之一。在第一种情况下，项

被加入到前因子导致经修改的右侧函数

第二选择导致需要加入到前因子的附加项

右侧估算的高效算法

利用上述步骤，显示出一种方案，其根据以下计算序列以最小数量的代数运算估算(30b)：

1.输入：和参数：λ、C_I、α

辅助变量：λ₊，λ_-，f_s，φ₁，ζ，ξ₁，...，ξ₆

2.由利用(1a)计算的有效速度梯度矩阵初始化

${\hat{N}}_a\leftarrow {\lambda }_{+}(▿\otimes U)+{\lambda }_{-}{(▿\otimes U)}^T=\hat{M},$λ_±：＝1/2(λ±1)

3.计算对称矩阵${\hat{G}}_M\leftarrow {\hat{N}}_a+{\hat{N}}_a^T$和它的迹${\xi }_1\leftarrow \mathrm{Tr}\left({\hat{G}}_M\right).$

4.根据压缩$\zeta \leftarrow {\hat{G}}_M:{\hat{G}}_M$由${\xi }_2\leftarrow (C_I/\lambda )\sqrt{2\zeta }$计算ξ₂＝2D_r。

5.计算标量取向因子$f_s\leftarrow 1-27\det \left({\hat{a}}^{\left(2\right)}\right),$其括号形式${\tilde{f}}_s\leftarrow \min (1,\max (f_s,0))\in \left[\mathrm{0,1}\right]$以及项${\xi }_3\leftarrow (1-{\tilde{f}}_s)/7.$

6.计算压缩${\xi }_4\leftarrow {\hat{G}}_M:{\hat{a}}^{\left(2\right)}.$

7.通过以下公式使用存储在辅助变量ξ_k中的值根据(32)计算前因子：

8.计算项然后计算

9.使用辅助变量ξ₅←2ξ₃和计算矩阵应该通过以下运算序列进行计算：

${\hat{N}}_a\leftarrow {\hat{N}}_a-{\xi }_5{\hat{G}}_M\Rightarrow {\hat{N}}_a\leftarrow {\hat{N}}_a-{\xi }_6\hat{1}.$

10.通过计算$\hat{F}\leftarrow {\hat{a}}^{\left(2\right)}·{\hat{N}}_a+{\hat{N}}_a^T·{\hat{a}}^{\left(2\right)}$初始化右侧结果矩阵。然后相继完成右侧函数(30b)的估算：

11.结果：${\hat{F}}^{\left(2\right)}({\hat{a}}^{\left(2\right)},▿\otimes U)\leftarrow \hat{F}$

使用“简化符号”(CN)的向量形式FTE的重新表示。

右侧函数的高效估算过程的构建中的另一个要素基于以下事实：FO矩阵和右侧函数都是对称矩阵并且只有6个独立矩阵元。因此将它们作为一般3×3矩阵存储是不实用的，并且将估算所需要执行的运算作为对矩阵的运算来实现是低效的。根据本发明，通过使用以下识别方案将对称3×3矩阵作为实6分量向量进行处理：向量索引$\mu \leftrightarrow \left(\mathrm{ij}\right),$μ∈{1，...，6}，以及对称矩阵的索引对(ij)＝(ji)，其被称为“简化符号”(CN，参见例如[22])(在结构力学中，这也被称为“Voigt符号”)：

μ  1  2  3456(ij)＝(ji)  (11)  (22)  (33)(23)＝(32)(13)＝(31)(12)＝(21)

这个$\mu \leftrightarrow \left(\mathrm{ij}\right)$方案的特别选择当然只是6！＝720个等效方案中的一个。上面的表格显示通常的选择，这也在优选实施例的FO模块中采用。(可替选的变型通常选择非对角矩阵元对向量索引μ∈{4，5，6}的不同的分配。)使用CN方案产生对称3×3矩阵的元素到相应的6矢量(即的向量)的分量的双射映射如下面所示

这样，存储右侧函数估算结果的FO矩阵、矩阵和矩阵被分配到相应的向量。以同样的方式，如果四阶张量关于第一和第二索引对对称的话(即，如果满足$A_{\mathrm{ijkl}}^{\left(4\right)}=A_{\mathrm{jikl}}^{\left(4\right)},A_{\mathrm{ijkl}}^{\left(4\right)}=A_{\mathrm{jilk}}^{\left(4\right)}$)，则该张量的分量被分配到6×6矩阵的元素如果除此之外该张量还关于两个索引对的互换$\left(\mathrm{ij}\right)\leftrightarrow \left(\mathrm{kl}\right)$对称并且因此具有正交性对称(见方程(11))，则满足即矩阵本身对称并具有${\Sigma }_{\mu =1}^6\mu =21$个独立元素。通过进一步的对称性和张量分量间的代数关系，可以减少独立元素的数量。在全对称四阶张量的情况下，存在6个附加恒等式

其将的独立元素的数量减少到15。规范化条件$a_{\mathrm{ij}}^{\left(2\right)}={\Sigma }_{k=1}^3{A}_{\mathrm{ijkk}}^{\left(4\right)}={\Sigma }_{k=1}^3{A}_{\mathrm{ikjk}}^{\left(4\right)}$和迹条件${\Sigma }_{k=1}^3{a}_{\mathrm{kk}}^{\left(2\right)}=1$产生一组代数关系，将独立矩阵元的数量减小到更小的值(详见[22])。

引入CN方法导致如下的“向量形式”的FTE

对于任意矩阵3×3矩阵映射定义实对称3×3矩阵的六维向量空间中的线性映射，其可以被形式上写成由实6×6矩阵描述的中的矩阵矢量乘积其元素或者为零或者由的元素的简单代数项给出。

在这个意义上，CN方案产生标识其中

相似地，由以下公式以元素方式定义压缩运算

如CN方案所给出的，使用索引分配$\mu \leftrightarrow \left(\mathrm{ij}\right),v\leftrightarrow \left(\mathrm{kl}\right).$此运算还可以写成矩阵—向量乘积其中的矩阵元与的矩阵元经由如下形式相关

这揭示矩阵(与不同)不是对称的！根据CN方案(即)，矩阵通过分配张量的元素而产生，该张量本身通过闭合近似从FO矩阵计算得出，因此项代表CN中的闭合近似。

所提出的算法的计算消耗

针对右侧函数的高效估算而提出的过程的步骤2到步骤10中所执行的计算操作的最高效的实施通过将CN方案应用到此过程而产生。下表给出执行使用CN方法的过程的各个步骤所需要的运算(乘法和除法的数量，加法和减法的数量)的数量和函数调用的概况：

以下要点评述实施的各方面以及针对上述算法的各步骤在表中给出的运算计数：

·由于矩阵不是对称的，因此将其作为3×3矩阵存储。如果以的计算初始化则不需要的额外存储。

·用于的计算，对称矩阵需要为矩阵形式(→步骤9)，但是在最终步骤10中，其还作为6向量出现。

·存储和分配运算没有被计数，因为它们对总的计算消耗的贡献是可以忽略的。

·通过重新使用已经存在的变量(当不再需要其值时)来减少辅助变量的数量。这没有在上面提出的算法中执行，因为这会降低对算法进行的说明的清楚程度。

·对一对对称矩阵的压缩运算(如步骤4和步骤6中所做的)是使用CN通过7个运算和5个运算实现的。(在矩阵符号中，它由$\hat{m}:\hat{n}=\mathrm{Tr}({\hat{m}}^T·\hat{n})={\Sigma }_{i,j=1}^3{m}_{\mathrm{ij}}{n}_{\mathrm{ij}}$定义。)

·实对称3×3矩阵的行列式(→步骤5)是使用CN通过11个运算和5个运算计算的。

·如果使用辅助变量将作为对称3×3矩阵进行存储，则步骤9中执行的矩阵的组合需要6个运算和12个运算，外加2个乘法以计算变量ξ₅和ξ₆。

·矩阵运算$\hat{F}\leftarrow {\hat{a}}^{\left(2\right)}·{\hat{N}}_a+{\hat{N}}_a^T·{\hat{a}}^{\left(2\right)}$需要使用CN的27个运算和21个运算。

·任何加(或减)多个单位矩阵的类型的运算只需要β到对角元的三个加法。

表的最后一行示出获得右侧函数的单次估算需要96个运算、84个运算和3个函数调用的计算消耗以计算(双精度)浮点数值的平方根以及一对这种数值的最小值和最大值。在当前计算机硬件上，加法和乘法运算的计算消耗大致相同，对min或max函数的调用消耗大约1.5-2次运算，而平方根计算的消耗总共约为6-10次运算(即，平均8次运算，取决于所需要的精度)。

利用以上提出的算法实现的(30b)的单次估算总共需要190次运算的计算消耗。这是估算具有包括动态迹稳定化的(稳定化)混合闭合的FTE右侧的最高效的方法。

计算效率评估

根据与估算由方程(35)给出的向量形式的右侧函数的“标准”方法进行的比较产生对所提出过程的计算效率的评估。在以下情况下会需要采用这种方法，例如，如果编码被规定为具有特定模块结构，使得(i)以“一般”过程执行不依赖于闭合近似的所有运算，这种过程将6×6实矩阵作为输入，并且(ii)使用如上面所说明的CN方案和(35b)根据闭合函数的具体选择以单独的过程计算的矩阵元。可以通过如下考虑来估算这种“标准”方法的计算消耗：

·该“标准”过程将从计算矩阵6矢量和标量参数2D_r(见步骤2到步骤4)开始，加上平方根计算，这总共需要46次运算。

·如果事先解析地完成矩阵向量乘积并且只对所产生的针对六个分量的公式进行数值估算，则对应于的运算需要48次运算。

·矩阵量乘积以及将其从右侧减去的计算需要66+6＝72次运算。

·2D_r·1的加法和6D_r·a的减法(使用6D_r＝3/2·2D_r)需要另外3+12+1＝16次运算。

为平方根计算计8次运算，到这时该方法用于估算(35)右侧的计算消耗合计达190次运算。用于迹稳定化的控制项的加入增加了更多的运算(见步骤8)，但是通过将项3D_r＝3/2·2D_r从的对角线减去以代替从右侧直接减去6D_r·a可以节省大约同样多的运算。

这些“固定消耗”不依赖于计算矩阵的元素所需要的闭合近似，该计算需要在单独的过程中完成。

由于(6)和(6a/b/c)所定义的张量函数具有正交对称性，但由于二次项(6a)而不是完全对称的，因此在混合闭合的情况下矩阵有21个独立元素。将的定义转换到CN产生计算矩阵元时需要使用的如下一组表达式：

对称矩阵和都包含需要事先解析计算的全对称四阶张量表达式的分量。

当矩阵元为常量时，矩阵元为向量分量a_μ的线性表达式。矩阵元和的解析计算产生如下一对实对称6×6矩阵：

使用这些矩阵，根据如下过程确定组合矩阵的计算消耗：

·由有效算法的步骤5完成受限标量取向因子和辅助变量ζ₂的计算并且需要23次运算(为内嵌的min(0，max(...，1))调用算作3次运算)。由一次附加除法计算辅助变量ζ₁＝ζ₂/5，因此计算和ζ_1/2所需要的运算量为24。

·由初始化矩阵这对于上三角的21个独立矩阵元中的每个进行2次乘法，或总共42次运算。(利用对称性分配剩下的矩阵元。)

·的矩阵元由单次运算(即3ζ₁的计算)以及将ζ₁或3ζ₁到相应的非零矩阵元的若干分配来计算。相似地，整个更新运算只需要从的相应条目减去ζ₁或3ζ₁。由于只需要在的上三角元素执行这些减法，因此更新运算的总运算计数降到10次(即9次减法加上1次计算3ζ₁的运算)。

·矩阵包含12个不同的表达式，其中只有9个需要由单次加法或乘法计算。(剩下的通过分配获得。)由于12个表达式都不能被先验地假设为零，的计算需要12次乘法，因此对于该矩阵的独立元素，计算的总运算量为21次。随后通过分配得到剩下的矩阵元。

·的组合过程中的最终步骤由更新运算组成，其需要21次运算以将矩阵和的上三角相加并进行分配以获得剩下的下三角矩阵元。

根据(36)组合矩阵的完整过程按顺序执行以上所描述的步骤，其需要24+42+10+21+21＝118次运算。(应该注意到此过程也是以非常高效的方式建立的！)由于非对称矩阵是根据(35b)利用18个乘法从对称矩阵计算得出，因此用于上述计算矩阵的过程的总运算量为118+18＝136次。

当应用到混合闭合式时，以上所述的“标准过程”总共需要约326次运算。这达到所提出的高效过程的计算消耗的大约1.72倍(即超出大约70％)。使用高效过程以“标准”方法的计算消耗的大约60％产生对于具有稳定化混合闭合(包括动态迹稳定化)的FTE的右侧函数的单次估算的结果。

计算效率的这一提高主要是利用基于先前的解析计算来巧妙地重组右侧函数的具体代数结构而达到的，其中通过经由CN方案减少方程来最终达到优化的效率水平。

7.7迹尺度改变、不变量监测和相空间投影

任何FO矩阵需要具有非负特征值和单位迹(→2.4节)。这意味着对于FTE的解的一组代数限制，该FTE由此变为DAS(→2.6节)。除了针对迹的等式$S_a=\mathrm{Tr}\left(\hat{a}\right)=1,$针对特征值的非负性限制可以等效地以针对矩阵的第二和第三不变量$K_a=\frac{1}{2}[\mathrm{Tr}{\left(\hat{a}\right)}^2-\mathrm{Tr}(\hat{a}·\hat{a})]$和$D_a=\det \left(\hat{a}\right)$的一对不等式K_a≥0和D_a≥0(→方程(9)、(10))的形式进行公式化。

通过加入FTE右侧的适合的惩罚项经由动态迹稳定化(DTS)以“主动的”和计算上花费不多的方式将迹条件计算在内，其中该惩罚项自动地将数值解的迹保持在所需要的值S_a＝1附近。使用此技术仅获得近似地满足迹条件(即S_a≈1)的FTE的数值解。从实际的观点看，这几乎不影响解的正确性，因为总有可能通过乘以因子1/S_a来改变矩阵元的尺度。将这种迹尺度改变运算数学地描述为

并具有若干有利特性：

·如果μ_k为矩阵的特征值，改变尺度的矩阵的特征值由μ′_k＝μ_k/S_a给出，因此(i)如果S_a＞0，则它们不改变符号，并且(ii)如果S_a≈1(如DTS所保证的)，则它们的数值仅被轻微缩放。

·相应的特征向量不被尺度改变运算影响。

因此迹尺度改变的运算不导致计算的数值解的质的改变，而是在保留矩阵的所有本质特性的情况下对解进行轻微修正。从这个观点来看，FTE的具有非负特征值μ_k(或等效地：非负的第二和第三不变量K_a和D_a)且S_a≈1的任意数值解在实践中可以被看做是FO矩阵。

FTE的不变量控制的一般方面

根据[12]的结果(→第3节)，如果实对称3×3矩阵的迹为正并且第二和第三不变量K_a和D_a为非负，则保证了该矩阵特征值的非负性。

由于矩阵的迹是矩阵元素的简单线性函数，并且实质导数算子$\frac{D}{\mathrm{Dt}}=\frac{\partial }{\partial t}+U·▿$与迹运算对易，因此有可能以直接的方式得到迹的演化方程(→第4节)。在“精确”情况以及混合或一般正交闭合近似的情况下(见方程(14)和(15a/b))，由于FTE右侧的特殊代数结构，此演化方程采取闭合形式。由于FTE的这些特殊代数性质，有可能向FTE右侧本身加入惩罚项，其以特殊(预期的)方式(→DTS)影响迹的演化方程。

然而，与迹不同，第二和第三不变量K_a和D_a是矩阵元的非线性函数。因此不可能应用相同(或相似)的过程来获得针对这些不变量的闭合形式演化方程(通常不变量的演化方程的右侧是所有矩阵元的函数，不只是其不变的组合，并且不独立于编码于特征向量中的矩阵方向特性)。相似地，很难(如果不是不可能的话)在FTE右侧加入以预定方式影响非线性不变量演化的任何特定项。因此，通过FTE右侧适当的惩罚项以将不变量保持在非负区域内的第二和第三不变量的“主动”控制是不可用的。

一个可替选过程由“不受控”(或部分受控)积分过程与将数值解映射回其容许区域的投影运算的结合组成。这种“被动”过程导致有效的积分方案，如果数值积分的运算和随后的投影被应用在每个时间步，则该积分方案保持“不受控”方案的精度(具体讨论参见[9]的IV.4章关于投影方法的部分)。

应用到FTE，此结合的过程由以下方面组成

1.使用第5、6和7节(直到前面的7.6小节)所描述的方法进行的积分步骤，然后

2.不变量监测步骤，用来确定在先前的积分步骤中计算的矩阵是否仍为FO矩阵，以及

3.相空间投影步骤，如果监测步骤的结果为否定(即违反了不变量条件)，则将矩阵映射回相空间集M_FT。

由于相空间M_FT为凸集(→第3节中的定理1)，对于每个实对称3×3矩阵，M_FT中存在“最小距离”(以适合的度量测量，见下面)上的唯一的矩阵。这定义了上述过程的第三步骤所需要的投影映射。下面讨论相空间投影的更多细节。

高效不变量监测和特征值计算

根据上述过程，需要在每个时间步和计算区域的每个(填充)单元中执行不变量监测。因此不变量(第二和第三不变量或特征值，其本身也是不变量)的计算需要以最高效的方式实现以避免任何不必要的计算开销。使用符号(x，y，z)和(u，v，w)来表示第3节中引入的矩阵的对角元和非对角元，由下列公式将不变量作为矩阵元的多项式表达式给出

$\left(38\right),\left(\begin{array}{c}{S}_a=x+y+z\\ K_a=\mathrm{xy}+\mathrm{yz}+\mathrm{zx}-(u^2+v^2+w^2)\\ D_a=\mathrm{xyz}+2\mathrm{uvw}-(u^{2}x+v^{2}y+w^{2}z)\end{array}\right)\mathrm{for}\hat{a}=\begin{array}{c}\left(\begin{array}{ccc}x& w& v\\ w& y& u\\ v& u& z\end{array}\right)\end{array},$

其总共可以用11次加法和13次乘法估算。由于在FTE的数值积分过程中的每个时间步都要计算迹S_a，因此可以对于每个单元以每时间步另外的22次运算的消耗来完成通过检查不等式K_a≥0和D_a≥0进行的不变量监测。

对不变量K_a和D_a的检查是用来确定实对称3×3矩阵是否属于FO矩阵的M_FT集的最高效的方式，因为以条件μ_k≥0对特征值进行的检查需要对后者进行详尽的计算。

具体分析显示可以通过计算特征多项式P_a(μ)的根来最高效地完成矩阵的特征值的数值计算，其在计算作为P_a(μ)的系数的不变量所需要的24次运算之外消耗大约100次运算，因此其消耗高出大约5倍。这量化了使用特征值代替使用不变量用于监测FO矩阵特性的计算开销，并且突出了使用不变量(38)是高效监测的谱性质所选择的方法。矩阵的数值对角化是最昂贵的计算特征值的方法。因此不推荐这种过程，除非也关心特征向量(例如作为下面讨论的相空间投影过程的一部分)。

相空间投影

如果不变量监测指示在当前时间步中FTE的数值解移出了相空间集M_FT(即违反不等式K_a≥0和D_a≥0中的一个或两个)，则需要通过向相空间M_FT的投影校正该解。

对于任意3×3实矩阵通过定义“最接近”矩阵$\hat{a}\in M_{\mathrm{FT}}$来产生该投影映射，其中以Frobenius规范${\left|\right|\hat{m}\left|\right|}_F=\sqrt{\mathrm{Tr}({\hat{m}}^T·\hat{m})}$测量两个3×3实矩阵的距离$d_F(\hat{m},\hat{n})={\left|\right|\hat{m}-\hat{n}\left|\right|}_F,$其中Frobenius规范由3×3实矩阵的向量空间中的标量积$<\hat{m}|\hat{n}>=\mathrm{Tr}({\hat{m}}^T·\hat{n})$自身导致。这些量对对称矩阵的六维子空间的约束是直接的。可以显示出[18]对于任意3×3实矩阵存在唯一矩阵${\hat{a}}^{*}\in M_{\mathrm{FT}}$使得和M_FT的元素之间的距离准确地在$\hat{a}={\hat{a}}^{*}$的情况下最小化。

由矩阵${\hat{a}}^{*}\in M_{\mathrm{FT}}$给出该最小距离解，该矩阵由以下特性唯一地定义(其形式证明参见[18])：

·的特征向量与矩阵的特征向量相同。

·的三个特征值(μ₁^*，μ₂^*，μ₃^*)是三角集(见第3节中的等式(7)和图7)

的唯一点，该点具有中距离由矩阵的三个特征值给出的点(μ₁，μ₂，μ₃)最小的欧几里德距离。

在对称3×3矩阵的六维实矢量空间上定义的并且具有M_FT中的值的投影映射由如下公式给出简洁的数学符号

并且实际计算投影映射的值的算法由以下步骤组成：

1.给定实对称3×3矩阵作为输入，计算其不变量S_a、K_a和D_a。

2.检查条件S_a＝1、K_a≥0和D_a≥0。如果所有条件都满足，则$\hat{a}\in M_{\mathrm{FT}}$已经成立，并且投影映射仅是等式：

在这种情况下将输入矩阵分配给输出并退出。否则(即如果K_a＜0或D_a＜0)执行下列步骤：

3.计算该矩阵的特征值μ_k和相应的特征向量E_k(即$\hat{a}·E_k={\mu }_k{E}_k$)。

此步骤可以形式地表示为：

${\{{\mu }_k,E_k\}}_{k=\mathrm{1,2,3}}\leftarrow \hat{a}={\Sigma }_{k=1}^3{\mu }_k{E}_k\otimes E_k,$

$\left(39b\right),{\hat{R}}^T·\hat{a}·\hat{R}=\mathrm{diag}({\mu }_1,{\mu }_2,{\mu }_3),$

$R_{\mathrm{jk}}={\left(E_k\right)}_j\iff \hat{R}=(E_1,E_2,E_3)$

4.然后寻找对于中的μ＝(μ₁，μ₂，μ₃)具有最小欧几里德距离的唯一点${\mu }^{*}=({\mu }_1^{*},{\mu }_2^{*},{\mu }_3^{*})\in {\Delta }_a,$即：

$\left(39c\right),{\mu }^{*}=\arg \underset{p^{\prime }\in {\Delta }_a}{\min }{\left|\right|{\mu }^{\prime }-\mu \left|\right|}_2,{\left|\right|{\mu }^{\prime }-\mu \left|\right|}_2={\left[{\Sigma }_{k=1}^3{({\mu }_k^{\prime }-{\mu }_k)}^2\right]}^{1/2}$

5.最后通过以下公式组合投影映射的值

注意到如果$\hat{a}\in M_{\mathrm{FT}}$则(39c)产生μ^*＝μ，因此在这种情况下(39d)符合(39a)。在这个意义上讲，(39b/c/d)已经表示了投影映射的一般定义。仅对于$\hat{a}\notin M_{\mathrm{FT}}$的情况需要执行计算(39b/c/d)的结果的各个步骤，而在$\hat{a}\in M_{\mathrm{FT}}$的情况下(39a)使该定义形式上完整。

虽然上面给出的计算的算法的描述是数学上精确的，但是该算法的实际实施取迹S_a≈1的输入矩阵(如具有DTS的FTE时间积分方案所提供的)，这将略去本过程的步骤1中对迹的附加估算。另外，实际实施将并入各个算法步骤的下列变化和/或说明，见图13：

·在步骤1中，由(38)计算输入矩阵的不变量K_a和D_a。

·在步骤2中，只检查条件K_a≥0和D_a≥0，并且如果两个条件都满足，则使输入矩阵保持不变并将其分配到输出。

·通过数值矩阵对角化获得步骤3的结果，例如通过迭代循环雅可比旋转(见[38]，11.1章或[39]，6.5.2章)或单Givens/Housholder约简步骤以及随后的QR迭代(见[42]或[38]的11.2章)。

·实际上只在至少一个输入特征值μ_k为负时执行步骤4，所以求解点${\mu }^{*}=({\mu }_1^{*},{\mu }_2^{*},{\mu }_3^{*})\in {\Delta }_a$总位于三角集Δ_a的边缘(包括角点)。

使用“迹尺度改变的”特征值${\tilde{\mu }}_k={\mu }_k/S_a$(见方程(37)和随后关于此问题的叙述)，将根据(39c)确定μ^*所需要解决的最小化问题限制到由三角形Δ_a的三个角点(1，0，0)、(0，1，0)和(0，0，1)固定的平面S_a＝1。可以在中等效地解决此“平面”最小化问题，例如在(μ₁，μ₂)平面中，通过寻找位于投影三角形边缘上的唯一的对(μ₁^*，μ₂^*)

(即Δ_a到(μ₁，μ₂)平面的正交投影)，该唯一的对具有到对的最小的欧几里德距离。μ^*的第三坐标因此由${\mu }_3^{*}=1-({\mu }_1^{*}+{\mu }_2^{*})$给出。因为通过DTS，S_a≈1，解决“平面”最小化问题的算法产生对(39c)的精确解的很好的近似，而这比在中直接根据(39c)计算μ^*的算法简单的多并且能够更高效地实施。

·在步骤5中，计算矩阵的各个元素的公式由$a_{\mathrm{ij}}^{*}={\Sigma }_{k=1}^3{\mu }_k^{*}{R}_{\mathrm{ik}}{R}_{\mathrm{jk}}$给出，其以旋转矩阵$\hat{R}=(E_1,E_2,E_3)$的矩阵元R_ij的方式给出，该矩阵元已经在步骤3中事先算出。

8.实施

该过程的结果为描述物件的给定区域中纤维取向可能性的分布函数，该结果以图形或数字的形式显示在计算机工作站(未示出)的显示器上。这可以是执行计算的计算机或工作站的显示，或者可以在与执行模拟的计算机网络连接的计算机的显示器上。

模具或产品的设计者将使用该模拟的结果来改进由注塑成型过程得到的物件的质量。也可以由工程师利用该模拟的结果来降低用于制造所考虑的物件的成本。关于纤维信息的可靠信息的优点使得工程师能够在确定该物件的强度、刚性和稳定性特性之前得到更好的理解和信息，因为纤维(纤维通常比聚合物材料具有更高的强度特性)的取向对于物件的强度、刚性和稳定性特性具有决定性的影响。另外，纤维取向对于对其中悬浮有纤维的聚合物团进行注塑成型时可能发生的翘曲效应有影响。通过得知纤维取向，造成影响的翘曲和其它应力能够被更好地预测或者设计改变以避免。

该模拟的结果还可以被传输到其它应用软件，例如CAD软件。这样，模拟的结果可以用于设计软件和模拟软件之间的交互式过程。

本发明具有很多优点。不同的实施例或实施可以产生一个或更多如下优点。应该注意到，这不是完全的列举，可以有此处未描述的其它优点。本发明的一个优点是可以以显著减小的计算工作量确定产品的纤维增强物件中的纤维取向分布。本发明的优点是可以以提高的精度确定产品的纤维增强物件中的纤维取向分布。本发明的另一个优点是可以更快地确定产品的纤维增强物件中的纤维取向分布。

虽然为了说明的目的具体描述了本发明，应该理解这些细节只是为了说明的目的，本领域技术人员可以在其中做各种变型而不脱离本发明的范围。

权利要求中使用的术语“包含”并不排除其它单元或步骤。单称并不排除复指。

参考文献

S.G.Advani(Ed.)：Flow and Rheology in Polymer CompositesManufacturing，Elsevier，阿姆斯特丹(1994)

G.B.Jeffery：The motion of ellipsoidal particles immersed in a viscousfluid，Proc.R.Soc.A 102，161-179(1922)

M.Junk和R.Illner：A new derivation of Jeffery′s equation，J.Math.Fluid Mech.8，1-34(2006)

F.Folgar和C.L.Tucker III：Orientation behaviour of fibers inconcentrated suspensions，J.Reinf.Plast.Compos.3，98-119(1984)

C.L.Tucker和S.G.Advani：Processing of short-fiber systems，ch.6，pp.147in S.G.Advani(Ed.)：Flow and Rheology in Polymer CompositesManufacturing(见以上参考文献[1])

S.G.Advani和C.L.Tucker III：The use of tensors to describe andpredict fiber orientation in short fiber composites，J.Rheol.，751-784(1987)

S.G.Advani和C.L.Tucker III：Closure approximations forthree-dimensional structure tensors，J.Rheol.，367-386(1990)

J.S.Cintra和C.L.Tucker III：Orthotropic closure approximationsfor flow-induced fiber orientation，J.Rheol.39，1095-1122(1995)

E.Hairer、C.Lubich和G.Wanner：Geometric numerical integration，Springer，柏林(2002)

J.Linn、J.Steinbach、A.Reinhardt：Calculation of the 3D fiberorientation in the simulation of the injection molding process forshort-fiber reinforced thermoplasts，ECMI 2000 Conference，Palermo(2000)

J.Linn：Exploring the phase space of the Folgar-Tucker equation，SIAM-EMS Conference，柏林(2001)

J.Linn：On the frame-invariant description of the phase space of theFolgar-Tucker equation，p.327-332in A.Buikis，R.A.D.Fitt(Eds.)：Progress in Industrial Mathematics at ECMI 2002，Springer(2004)

S.Hess：Fokker-Planck equation approach to flow alignment in liquidcrystals，Z.Naturforsch.31A，1034ff.(1976)

M.Doi：Molecular dynamics and rheological properties ofconcentrated solutions of rodlike polymers in isotropic and liquidcristalline phases，J.Polym.Sci.，Polym.Phys.Ed.19，229-243(1981)

M.Grosso、P.L.Maffetone、F.Dupret：A closure approximation fornematic liquid crystals based on the canonical distribution subspacetheory，Rheol.Acta39，301-310(2000)

M.Simple models for complex nonequilibrium fluids，Phys.Rep.390，453-551(2004)

J.Linn：The Folgar-Tucker Model as a Differential Algebraic Systemfor Fiber Orientation Calculation，pp.215-224 in Y.Wang，K.Hutter：Trends in Applications of Mathematics to Mechanics，proceedings ofthe STAMM 2004 conference in Seeheim(德国)，Shaker(2005)

U.Strautins：Investigation of fiber orientation dynamics within theFolgar-Tucker model with hybrid closure，硕士论文，数学系，凯泽斯劳滕大学(2004)

J.Linn：Dreidimensionale Vorausberechnung der anisotropenMaterialeigenschaften in thermoplastischen Spritzgusserzeugnissen(AnIso-3D)，Projektbericht für die MAGMA GmbH，Teil IIa(2002)

J.Linn：Dreidimensionale Vorausberechnung der anisotropenMaterialeigenschaften in thermoplastischen Spritzgusserzeugnisse(AnIso-3D)，Projektbericht für die MAGMA GmbH，Teil IIb(2002)

B.E.Ver Weyst、C.L.Tucker、P.H.Foss、J.F.O′Gara：Fiberorientation in 3D injection moulded features：prediction andexperiment，Internat.Polymer Processing 14，409-420(1999)；

B.E.VerWeyst：Numerical predictions of flow-induced fiberorientation in three-dimensional geometries，博士论文，伊利诺伊大学厄本那-香槟分校(1998)

G.I.Marchuk：Splitting and Alternating Direction Methods，pp.197-462in P.G.Ciaret&J.L Lions(Eds.)：Handbook of NumericalAnalysis，卷I，北荷兰，Elsevier(1990)

K.W.Morton：Numerical solution of convection-diffusion problems，Chapman&Hall，伦敦(1996)

R.J.LeVeque，Numerical Methods for Conservation Laws，(1992)

G.Strang：″On the construction and comparison of differenceschemes″，SIAM Journ.Num.Anal.5，506-517(1968)

M.G.Crandall和A.Majda：″The method of fractional steps forconservation laws″，Math.Comp.34，285-314(1980)

H.V.Kojouharov、B.M.Chen：″Nonstandard methods for theconvective transport equation with nonlinear reactions″，Numer.Meth.Partial Diff.Eq.14，467-485(1998)；″Nonstandard methods for theconvective-dispersive transport equation with nonlinear reactions″inR.E.Mickens(ed.)：Applications of non-standard finite differenceschemes，minisymposium on nonstandard finite difference schemes：theory and applications，SIAM年会，亚特兰大GA，USA 1999，由新加坡：World Scientific(2000)发表

H.Wang、X.Shi和R.E.Ewing：″An ELLEM scheme formultidimensional advection-reaction equations and its optimal-ordererror estimate″，SIAM J.Numer.Anal.38，1846-1885(2001)

P.J.van der Houwen：″Note on the time integration of 3Dadvection-reaction equations″，J.Comput.Appl.Math.116，275-278(2000)

W.Hunsdorfer，J.G.Venver：″A note on splitting errors foradvection-reaction equations″，Appl.Numer. Math.18，191-199(1995)

S.V.Patankar：Numerical heat transfer andfluidflow，HemispherePubl.Corp.(1980)

C.A.J.Fletcher：Computational techniques for Fluid Dynamics，卷I：Fundamental and General Techniques(第二版)，Springer(1991)

L.F.Shampine：″Conservation laws and the numerical solution ofODEs″，Comp.Math.Applic.12B，1287-1296(1986)

L.F.Shampine：″Linear conservation laws for ODEs″，Comp.Math.Applic.35，45-53(1998)

L.F.Shampine：″Conservation laws and the numerical solution ofODEs，part II″，Comp.Math.Applic.38，61-72(1999)

J.Linn：″Entwicklung eines Software-Moduls zur Berechnung derFaserorientierung in der Spritzgieβsimulation mit SIGMASOFT″，technical Report for the MAGMA GmbH(2001)

W.H.Press、S.A.Teukolsky、W.T.Vetterling、B.P.Flannery：Nunlerical Recipes in Fortran 77：The Art of Scientific Computing(第二版)，剑桥大学出版社(1992)

J.Stoer、R.Bulirsch：Introduction to Numerical Analysis(第三版)，Springer(2002)

D.H.Chung、T.H.Kwon：″An invariant-based optimal fittingclosure approximation for the numerical prediction of flow-inducedfibre orientation″，J.Rheol.46，169-194(2002)

F.Dupret，V.Verleye，B.Languilier：″Numerical prediction ofmoulding of short-fibre composite parts″，Proc.1st ESAFORM Conf.，291-294(1998)

G.H.Golub，H.A.van der Vorst：″Eigenvalue Computation in the20th Century″，J.Comp.Appl.Math.123，35-65(2000)