第 7 章：改进的纤维取向建模
7.1 Folgar-Tucker 模型的局限性介绍
许多用于预测短纤维增强热塑性塑料中纤维取向的模拟采用的是Folgar-Tucker模型[122]，该模型在第5章中有所描述。Folgar-Tucker模型在Jeffery方程中加入了各向同性旋转扩散项，以考虑纤维-纤维相互作用。扩散系数被假设为广义剪切率的函数，这一假设被认为是有道理的。然而，模型将纤维-纤维相互作用的强度仅仅归结为一个单一的标量常数，即相互作用系数CI，作为拟合参数使用。由于没有基本依据来确定CI的值，Tucker及其同事通过调整CI使其与实验测得的纤维取向张量分量a11相匹配来确定其值。然而，其他分量并不保证同时与实验数据匹配。由于纤维取向受多种因素影响，扩散系数有时可能被不恰当地用于修正其他因素引起的误差。此外，如前所述，必须使用闭合近似将四阶张量与二阶张量相关联。CI的值也被发现会因不同的闭合近似方案而显著变化（Cintra和Tucker[62]）。Fan等人[103]和Phan-Thien等人 [290] 开发了一种直接仿真技术，无需采用旋转扩散假设或闭合近似。实际上，这种纤维级仿真技术为 Folgar-Tucker 模型的问题提供了见解，并提供了一种寻找模型本构修正的方法。上述作者意识到，仅使用单一标量相互作用系数，Folgar-Tucker 模型无法预测直接仿真结果。因此，相互作用系数应为各向异性的张量形式。他们将传统的 Folgar-Tucker 等各向同性旋转扩散模型推广为各向异性旋转扩散（ARD）模型（Phan-Thien 等人 [289]）。此外，还发现传统的 Folgar-Tucker 模型预测的芯区宽度远小于注射成型端封条的实验数据，表明模型中取向变化的动力学过程过快，与实际情况不符。Folgar-Tucker 模型中的 Jeffery 项被认为是导致这一局限性的原因。Tucker 及其同事 [296, 372, 393]、Sepehr 等人 [327] 和 Férec 等人 [116] 提出了通过修改 Jeffery 项来减缓取向变化的一些模型。还有一些其他情况表明，Folgar-Tucker 模型仍需进一步改进。
例如，实验表明，稀悬浮液（约5%质量分数）的取向行为与浓悬浮液（约30%质量分数）有显著不同，而Folgar-Tucker模型的数值预测并未反映出这些变化（As-Sultany [12]）。该理论模型在更浓的系统中（约30%质量分数）表现更好，而在较稀的系统中（<10%质量分数）表现较差。模型的不足可能是因为忽略了模型中的粘弹性溶剂效应，因为Folgar-Tucker方程是针对牛顿流体中刚性颗粒的模型。在简单的剪切流中，当液体是粘性的时，细长颗粒会经历一种称为Jeffery轨道的周期性运动。最近关于非牛顿液体中纤维悬浮液的流变行为的研究表明，细长颗粒的旋转会由于粘弹性效应偏离Jeffery轨道（Fan等 [102]）。这一话题将在第14章中讨论。
7.2 不均匀旋转扩散（ARD）模型
7.2.1 进化方程
如前所述，对于非稀悬浮液流动，Folgar和Tucker [122] 将纤维-纤维相互作用建模为随机力。换句话说，纤维之间的空间相互作用通过Jeffery理论框架中的布朗过程类比来建模。这导致纤维在其配置p空间中的旋转扩散。
在该模型中，扩散系数
Phan-Thien 等人 [289] 给出的具有各向异性扩散的取向张量演化方程形式如下：
其中，
对于
假设C是应变率张量及其二阶不变量的函数，如下所示：
其中，系数
Phelps 和 Tucker 提出了修正后的各向异性旋转扩散模型，其形式为：
Phelps 和 Tucker [296] 进一步假设 C 也可以是纤维取向张量的函数，形式如下：
其中
7.2.2 直接模拟
从 Folgar-Tucker 类型模型计算得到的纤维取向张量暗示了一种取向概率。
个体纤维在悬浮液中的运动细节并未计算。Fan等人[103]开发了一种纤维级模拟方法，称为直接模拟，在这种方法中，每根纤维在悬浮液中的状态被明确地进行数值计算。流体动力相互作用被建模为润滑力引起的短程相互作用和通过细长体近似（Batchelor [27]）引起的远程相互作用的叠加。当纤维接触时，润滑力变得无穷大。在模拟中，可以通过假设周期条件来模拟纤维悬浮液，初始时纤维在参考单元中随机分布，然后在三维空间中复制单元。在他们的论文中，Fan等人[103]在不同的剪切速率下模拟了一个单位单元中的300根纤维。结果表明，在零剪切速率下随机取向的纤维在受到应变时会被重新排序（图7.1）。（图片来自Fan等人[103]，经Elsevier许可）
图 7.1：不同应变下纤维构型的直接模拟结果，转载自 Fan 等人 [103]，经 Elsevier 许可
直接模拟由于计算量过大，无法用于建模实际的注塑成型过程；然而，它提供了对纤维-纤维相互作用的见解，因为在模拟中相互作用的效果自然发生，无需进行闭合近似。模拟数据可以适当平均，以提供宏观性质，如取向张量和模拟悬浮液的减缩粘度，因此可以用来预测相互作用系数张量的各个分量。
7.2.3 交互系数的计算
为了预测交互系数张量的各分量，我们将式(7.2)在稳态下写作如下形式：
通过使用式(7.7)，一旦从直接模拟的取样数据中获得取向张量的所有分量，就可以确定交互系数张量的六个独立分量。如果模拟使用标准的Folgar-Tucker模型，交互系数的标量值
使用这种方法进行的计算显示出与实验结果的良好一致。例如，当悬浮液的体积分数为0.184，纤维的长宽比为12.4时，模拟得到的
如第5章所述，Bay [29] 和 Tucker 和 Advani [370] 解释了 C I 减小的“笼效应”。如果这是真的，直接模拟可能会忽略这一效应，因为在模拟中并未考虑纤维直径。此外，模拟的剪切流动发生在无界空间中，而 Bay 的实验中肯定存在壁效应。
7.3 减少应变闭合（RSC）模型
实验证据表明，浓悬浮液中取向变化的动力学远慢于标准 Folgar-Tucker 模型的预测。Huynh [170] 将此归因于纤维以集群形式移动并经历的局部应变小于总体应变。他引入了一个“减少因子”到 Folgar-Tucker 方程中的速度梯度中，从而有效减缓了取向张量的变化率。
鉴于非均匀运动，Sepehr 等人 [327] 修改了标准 Folgar-Tucker 模型，引入了滑移变形。这实际上等同于 Huynh 作出的修改。在两种情况下，修改后的 Folgar-Tucker 模型不再遵循物质客观性原理——该原理断言，材料对给定经历或运动历史的响应必须独立于任何参考框架的变化 [356]。Tucker 及其同事（Tucker 等人 [372] 和 Wang 等人 [393]）提出了一种新的修改模型，该模型是客观的。
该方法基于谱分解定理，允许将Folgar-Tucker方程分解为两个独立的微分方程组，分别对应于a_ij的特征值和特征向量。然后通过将速度梯度减少一个因子κ（≤1）来修改特征值的动力学表达式，而特征向量的动力学表达式保持不变。修改后，重建了材料导数D a_ij/Dt的演化方程。该新模型可写为：
其中
且
其中κ是一个表征参数，用于减少纤维取向的速率。当κ=1时，恢复了标准的Folgar-Tucker方程。向量e_i（i=1,2,3）是二阶取向张量的特征向量，λ_i（i=1,2,3）是相应的特征值。注意，方程7.11未使用求和约定，所有求和均明确表示。将方程7.9与标准的Folgar-Tucker模型进行比较，可以看到四阶张量a_ijkl被替换为A_ijkl，且C_I乘以κ。由于a_ijkl是通过闭合近似来近似的项，这种修改等效于使用新的闭合近似。这就是为什么该模型被称为“减少应变闭合”模型。
类似模型也独立地由Férec等人[116]提出。减应变闭合（RSC）模型被整合到各向异性旋转扩散（ARD）模型中，形成了一个统一的演化方程，称为ARD-RSC模型[296]。
7.4 悬浮体流变学与纤维流相互作用
悬浮体流变学不是本书第一部分的主题，因为在注射成型流动的模拟中忽略了纤维流相互作用，采用了解耦的方法。解耦方法独立计算流动场，然后利用获得的流动场进行纤维取向计算。根据Tucker[369]的观点，这种简化是否允许取决于无量纲数Npδ²，其中Np是体积分数和长宽比变化的粒子数，可以通过以下公式估算：
表7.1 A_i的渐近值，i=1到4
| 情况 | A1 | A2 | A3 | A4 |
|---|---|---|---|---|
| a_r² | 6 ln 2a_r - 11 | 3a_r² | 2(ln 2a_r - 1.5) | 2 ln 2a_r - 0.5 (杆状) |
| a_r = 1 + ε, ε ≪ 1 | 147ε | 14ε - 588ε² (1 - 7ε + 3ε²) | 9ε | 9ε |
| a_r → 0 | 3πa_r + 9π² - 2 - 3πa_r + 1 - 9π² | 3πa_r - πa_r² (盘状) |
参数δ描述了平面外纤维取向，等于流动通道狭窄度的最大值（由典型间隙高度与典型平面维度之比表征）和C_I^(1/3)，其中C_I是相互作用系数。当Npδ² ≪ 1时，允许采用解耦方法。
在这一范围之外，需要考虑纤维流动耦合。Lipscomb等人[227]还表明，在复杂流动中，流动动力学可能会受到纤维取向的强烈影响，因此分离计算可能会严重错误。为了考虑纤维-流动相互作用，必须考虑悬浮体的流变学。一般来说，对于牛顿流体中的纤维悬浮液，其本构方程包含两部分，贡献于悬浮液的额外应力：
其中 τ(s)ij 是悬浮流体的粘性贡献，τij 是颗粒贡献的应力。各种本构模型之间的差异在于颗粒贡献应力的表达式。尽管关于纤维悬浮液的流变学研究始于20世纪60年代，但更广泛的分散体流变学主题可以追溯到1905-1906年由Einstein[94]的工作，他考虑了稀悬浮液（φ < 0.03）中刚性球体在牛顿流体中的情况。Einstein给出的有效悬浮液附加应力的结果为
其中 ηs 是溶剂粘度。Einstein方程描述的流体是各向同性的。Ericksen[97]和Hand[145]推导出了一种各向异性稀悬浮液模型，称为横向各向同性流体（TIF）模型。
方程写作如下：
其中，
其中，
对于完全对齐的取向，有
对于浓悬浮液，TIF 模型的另一种修改由范-廷恩和格雷厄姆 [292] 提出。
他们的修改理念是保留TIF模型中的主导项，并通过体积分数的函数依赖关系来替代TIF模型中体积分数的线性依赖关系。在Phan-Thien-Graham模型中，颗粒贡献的应力由下式给出：
其中，
其中，
Fan等人[104]进一步修改了Phan-Thien-Graham模型，加入了由于纤维随机运动引起的动量传输贡献（熵贡献）的额外项。
7.5 布朗运动模拟
当以二阶取向张量作为主要变量求解演化方程时，需要进行闭合近似，这会引入纤维取向计算中的误差。Fan等人[104]采用了一种不同的计算技术——布朗运动模拟。在此方法中，纤维的运动被描述为一个随机过程。
纤维n的取向矢量pi(n)的随机方程可写为：
其中，n = 1, …, N，N为纤维总数，Lij为之前定义的有效速度梯度张量，F(b)(t)为随机力，其性质为：
上述方程中，D(r)为扩散系数，δ(s)表示狄拉克δ函数，I为单位张量。当选择以下表达式作为D(r)时，随机方程的平均结果将等同于福尔加-塔克方程：
其中，
其中，wi为维纳过程。解方程7.21得到pi(n)后，可以使用集合平均计算结构张量：
其中，N为集合中的纤维总数。
一种结合布朗动力学模拟方法与有限元方法来求解守恒方程的技术被称为 CONNFFESSIT（非牛顿流体流动的有限元和随机模拟技术）[280]。该技术使得可以直接从微观模型模拟宏观流动场，而不是依赖闭式本构方程。这是一种多尺度模拟。对于许多感兴趣的流动问题，单个颗粒跟踪涉及大量的自由度，因此在计算时间上较为昂贵。Hulsen 等人[169] 将 CONNFFESSIT 方法扩展到了所谓的布朗配置场方法，以减少方差。Fan 等人[104] 采用这一思想来模拟非稀释纤维悬浮流动。他们将纤维运动的随机向量 p(n) 视为空间 x 和时间 t 的随机向量场函数，并使用欧拉时间导数来替换 p(n) 的拉格朗日时间导数。为了解决方程 7.21，引入了一个新的配置场 q(n) (x, t) (n = 1, …N)，它与 p(n) (x, t) 平行：
(7.28) 将随机力表示为白噪声，并对时间从
式中，
尽管布朗动力学模拟或布朗配置场方法比直接从演化方程计算