注塑成型纤维取向模拟：从Jeffery方程到方向张量的工程实践

简介：面向不连续纤维复合材料流动诱导取向模拟的Matlab函数工具包，适合从事短纤维注塑成型、纤维取向演化及微观力学性能预测的研究与工程人员。代码实现Jeffrey、Folgar-Tucker、Phelps-Tucker等多类取向模型，提供四阶方向张量闭合近似、平面/三维方向分布函数重建、纤维长度分布演化，以及基于Halpin-Tsai、Mori-Tanaka等平均场模型的刚度与热膨胀预测，功能覆盖从取向求解到复合材料层合板性能计算的全链条。资源共98个文件，包括63个m函数脚本、23个mlx实时脚本、3个tex文档及部分mat数据、PDF说明等，m文件承载核心算法，mlx文件便于交互式示例演示，压缩包仅3.47MB，轻量易用。已有245人学习下载，适合需要对比不同取向模型、快速获取可修改Matlab代码的进阶用户。

1. 流动诱导纤维取向为什么难算：从 Jeffery 方程到方向张量

注塑成型一个玻纤增强支架，最头疼的往往不是流动求解，而是“熔体固化后，纤维到底朝向哪里”。Jeffery 方程能够描述单根纤维在剪切流中的旋转，但真实零件里每立方毫米就有数百万根纤维，逐根追踪在工程时间上不可行。张量（第二阶和第四阶）是描述取向状态的标准工具，Fiber-Orientation-Tools 正是围绕这一抽象层展开的 MATLAB 函数集：从 Jeffery 方程与 Folgar-Tucker 扩散，连续到方向分布函数重建、平均场刚度计算，再升级到经典层压板理论。它适合注塑仿真的工程师、复合材料 CAE 分析者和想弄清参数连接但没有时间读理论全文的研究生——每步都可以单独跑，不用为整合全流程发愁。

2. pDot 与 AdotJeffQuad：Jeffery 方程与 Folgar-Tucker 扩散的数值边界

2.1 先看 Jeffery 方程怎么被写成向量函数

先说单纤维模型，这部分对应仓库里的 pDot.m。向量形式的 Jeffery 方程写成函数，输入是运动张量和纤维指向：

function dp = pDot(p, D, W, lambda) % p : 3x1 单位向量，代表一根纤维的空间指向 % D : 3x3 对称张量，应变率张量（流动输入） % W : 3x3 反对称张量，涡量张量 % lambda : 形状因子，(re^2-1)/(re^2+1)，短纤维 re 趋于 1 dp = W*p + lambda * (D*p - (p'*D*p)*p); end

这四条输入在注塑场景里都来自 CFD。应变率 D 和涡量 W 从速度梯度分解得到，形状因子按纤维长径比取，玻璃纤维小长径比通常在 0.95~0.99。W*p 是刚性旋转项，后面括号中的第三项才是应变驱动的校准。如果 p 恰好等于一根真实纤维，输出 dp 就是当前时刻该纤维指向的变化速率。

在调用的时候，建议对照 AdotJeffQuad.m 或 JefferyTensorEqn.mlx 的数值积分方式，二者的推进顺序一致：先把 D、W 插值到当前时刻，再沿时间积分。这里有个容易踩的坑：时间步长如果直接沿用 CFD 的原始输出间隔，在应变率峰值处 p 会跳向错误方向。我一般会把流动结果重采样到应变率变化平缓的时间轴上，再喂给积分函数。

2.2 从单纤维转向群体：Folgar-Tucker 扩散项的意义

单纤维的 Jeffery 方程有三个缺陷：没有考虑纤维间相互作用、忽略扰动扩散、也不能反映成型件中取向的分散程度。Folgar-Tucker 模型在 Jeffery 方程右边加上一项随机扩散力矩：

dA2/dt = (W·A2 − A2·W) + λ(D·A2 + A2·D − 2D:A4) + 2 CI γ̇ (I − 3A2)

工程上真正需要关注的是最后一项：CI 是相互作用系数，γ̇ 是等效剪切速率。CI 越大，纤维取向越分散；CI 太小会让取向过于锐利，在转角处形成不真实的取向集中。这是调试时优先调整的参数。

提示：CI 的经验范围通常在 1e-4 ~ 0.01，越大越接近高浓度熔体行为。若仿真结果中取向张量的各分量收敛但分布图明显偏窄，优先怀疑 CI，而不是网格。

工具包里对应的函数是 Adot2.m 和 AdotPlanar.m，把上述方程封装成可接收流动张量的形式。代码层的大致形态如下：

function dAdt = AdotJeffQuad(A2, D, W, lambda, CI, gammaDot) % A2 : 3x3 对称半正定张量，当前取向状态 % CI : Folgar-Tucker 相互作用系数 % gammaDot : 标量等效剪切速率 A4 = closeA4(A2, 'hybrid'); termRot = W*A2 - A2*W; termOmega = lambda * (D*A2 + A2*D - 2*tensorContract4(D, A4)); termDiff = 2*CI*gammaDot * (eye(3) - 3*A2); dAdt = termRot + termOmega + termDiff; end

第 4 行调用了闭合近似模块，第 5 行的反对称乘法保持张量对称性。常见做法是先用一个简单剪切流验证计算收敛，再投入真实型腔流场。如果观察到 A2 的对角线元素之和偏离 1，多半是闭合近似与张量积分方式不匹配，而不是函数本身有 bug。函数没有对 A2 做迹归一化，这是有意的：归一化要放在整个时间步迭代之后做，否则会引入额外误差。

2.3 慢动力学 SRF、RSC、RPR 扩展，为什么你要关心

摘要里提到的慢速动力学模型（SRF 慢响应、RSC 减速、RPR 偏好旋转分布）是对 Folgar-Tucker 扩散项的修正。标准 FT 模型在低剪切弱流动条件下会预测过快的取向分散，而实际材料在低剪切区受纤维网络约束更强，取向松弛更慢。SRF 的思路是对扩散项打折，RSC 则对旋转项加入一个时间常数延迟。

opts.Model = 'RSC'; opts.RSC_TimeConstant = 0.1; % 单位与流动时间尺度保持一致 opts.SRF_ShearThreshold = 0.01; % 低于此剪切率，扩散项打折 [A, psi] = solvePsi3D(D, W, tspan, A20, opts);

注意：SRF 和 RSC 不要同时开启，它们的修正对象有重叠，一起用会让参数不可辨识。我的标定习惯是：先用标准 Folgar-Tucker 跑出一个基线，再微调 CI 让取向分布合理，最后只用一个慢动力学修正残余偏差。反过来做会让 CI 和慢动力学参数相互补偿，拟合出来的组合没有物理意义。

3. closeA4 与 AdotPlanar：闭合近似的选择逻辑与平面取向的简化

3.1 为什么 A4 不可测，但演化方程必须要 A4

上一节的演化方程里包含 D:A4，这是应变率张量与四阶方向张量的双点积。问题的根源在于：A2 的演化方程需要 A4，A4 的演化方程又需要 A6，这形成了无穷的层次耦合。实际处理是引入闭合近似，用 A2 的代数函数来逼近 A4，从而截断这个层次。仓库里的 closeA4.m 是这一层的统一入口，接收 3x3 的 A2 矩阵和一个模型标志字符串，返回 3x3x3x3 的 A4。

function A4 = closeA4(A2, model) switch lower(model) case 'linear' A4 = linearClosure(A2); case 'quadratic' A4 = quadraticClosure(A2); case 'hybrid' f = 1 - 27*det(A2); A4 = f*linearClosure(A2) + (1-f)*quadraticClosure(A2); case 'orthotropic' A4 = orthotropicClosure(A2); case 'ibof' A4 = ibofClosure(A2); otherwise error('closeA4:unknownModel', '%s 不是可用闭合模型', model); end end

真正需要留意的在调用端：平面的 closeA4planar.m 与三维 closeA4.m 行为完全不同。混用会导致 A2 的迹性质被破坏，方向盘上出现负的特征值，这在后处理里表现为取向分布图上的“空洞”，不是物理现象，是闭包选错了。

3.2 线性、二次、混合闭包的适用条件

线性闭合对所有 A2 都是严格线性的，在各向同性取向附近准确，在取向集中时误差变大。二次闭合适用高度取向的排列，但会把各向同性状态放大成虚假的择优方向。混合闭合在二者之间按 det(A2) 的偏离程度插值，是工程短纤维模流场景里最稳的方案。

下面的表格是这套工具中 TestPlanarClosures.mlx 暴露的典型行为，可作为快速判断依据：

在接近真实注塑的剪切流场里，混合闭包的误差通常比线性闭合小 30% 甚至更多，额外的计算开销可以忽略。正交各向异性闭包比混合稍慢，但在拉伸占比高的流场中不会放大取向振荡，这是混合闭包的已知软肋。

注意：当 det(A2) 接近 1/27 时，混合闭合自动退化为线性闭合，这是正常行为，不用当作故障处理。

3.3 平面闭包的差异：Bingham、椭圆半径与自然闭包

对薄壁注塑件，取向演化本质上发生在平面内，直接用三维四阶张量既浪费又容易触发非物理的厚度方向取向。仓库中 closeA4planar.m 配合 fitBingham2D.m、fitERdistn2D.m 和自然闭包处理这类问题。平面形式与三维形式的本质区别在于：平面闭包假定 A4 在板厚方向不传播取向信息。

我的经验法则是：壁厚与纤维长度之比小于 2 时优先使用平面闭包；超过 2 必须用三维闭包，否则厚度剪切会显著低估法向取向分量，导致后续层压板分析中的弯曲刚度失真。如果你做的是汽车门板这类大平面薄壁件，平面闭包配合 Bingham 重建取向分布，计算效率能提升一个量级。

4. solvePsi2D 与 fitARD：方向分布函数的细节量与 ARD 参数拟合

4.1 张量不够用，为什么还要方向分布函数

A2/A4 是取向分布的低阶矩，只保留统计上的前两阶或四阶信息。当你想知道特定方向的纤维占比，或者在制件表面判断皮芯层结构，就需要完整的方向分布函数 ψ(θ,φ)。仓库里 solvePsi2D.m 求解平面方向的瞬态 Fokker-Planck 方程，solvePsi3D.m 处理三维情形。

求解的核心是离散分布函数的演化：在球面网格上离散方向空间，从初始分布出发，沿 Jeffery 特征线搬运密度权重，同时计入扩散通量。核心函数的大致结构是：

function psi = solvePsi2D(theta, phi, Dp, Wp, tspan, options) % theta, phi : 平面角度网格点 % Dp, Wp : 平面投影后的应变率和涡量 % tspan : 时间区间 [t0 tf] % options.fiberLambda : 纤维形状因子 % options.model : 'Jeffery' | 'Folgar-Tucker' | 'ARD' psi = zeros(length(theta), length(phi)); % 初始分布：各向同性或指定分布 % 时间推进：特征线搬运加上扩散项 end

函数返回离散分布而非张量，后续计算二阶量时需要配合 p2A.m 把离散分布重新汇集为张量。常见的误用是在 ODF 演进过程中直接把 A2 当输入传入下一次积分，高频密度信息在这一步被丢失，之后无论用什么闭包都无法恢复。

4.2 各向异性旋转扩散与 fitARD 的参数含义

Phelps-Tucker 模型把扩散系数从标量 CI 推广为张量 Dr，fitARD.m 支持五常数、iARD、pARD、MRD、Wang 两常数等命名模型。拟合代码的输入通常是实验测得的取向演化序列和对应的流动场：

params = fitARD(A2meas, Dmeas, Wmeas, 'Model', 'ARD5'); params.CI % 各向同性扩散部分 params.b1 % 剪切相关权重 params.c % 剪切速率辅助项

拟合的逻辑是构建最小二乘目标，让模型预测的 dA2/dt 与实验变化率之差最小，求解器用 lsqnonlin。拟合数据量超过 1e5 时，注意把无量纲化做好，用 matlab 优化工具箱的典型配置就能收敛。有个绕不开的坑：A2meas 与 Dmeas 必须在同一坐标系下。若从 Moldflow 导出的数据没有对齐坐标系，即使参数正确，残差也会表现出规律性偏移，看起来像模型没选对。

实际操作中，我会先用 Asteady.m 做稳态预测，用稳态结果判断粗参数范围，再交给 fitARD 精修。否则高维参数在非凸目标函数上直接优化，容易陷入局部极小，最后拟合出的五常数和物理直觉对不上。

4.3 瞬态分布求解的实用陷阱

瞬态求解的老问题集中在三处：球坐标极点奇异性、高频振荡、数值耗散。solvePsi3D.m 采用球面网格，极点附近网格密集，需要额外的扩散通量修正。推荐以 60x30 网格起步，先确认结果对称性，再用 120x60 复核极值，不要一开始就上全分辨率。

收敛判据建议用分布熵而不是某个方向的峰值：熵达到稳定值后再进入下一步 A2F 计算。如果看到熵一直在缓慢下降而峰值持续上升，大概率是网格方向的耗散项放大了非物理取向集中，这时增加扩散系数或加密网格都能缓解，但要注意别把真实物理和数值效应混在一起。

5. halpin 与 mori 与 lielens：平均场均匀化如何将取向转换成刚度

5.1 平均场方法的核心假设

有了取向信息之后，需要把纤维和基体混合物换算为等效各向异性刚度。仓库里的 diluteEshelby.m、mori.m、lielens.m 属于平均场类别，Halpin-Tsai 属于经验校正型。共同思路是把纤维视为嵌入基体的夹杂，通过 Eshelby 张量建立纤维与基体之间的应变放大关系。

平均场方法的适用边界在于：当纤维长径比超过 20 时，dilute 模型会明显低估相互作用，而 Mori-Tanaka 与 Lielens 的合理性更好。Halpin-Tsai 适合作为快速上界估计，但前提是纤维近似平行；取向分散的场景直接用会导致刚度高估。

5.2 从单向刚度到方向平均的组合流程

工具包的使用路径是：先计算单向刚度，再通过方向平均得到宏观刚度。起始部分代码块：

Ef = 72e3; Em = 3.0e3; % 单位 MPa，E 玻璃纤维与 PP 基体 vfiber = 0.30; CisoM = iso2C(Em, 0.35); % 各向同性基体刚度矩阵 C_fiber = eng2C(Ef, 0.22); % 纤维刚度矩阵，假定横观各向同性 C_ud_mt = mori(C_fiber, CisoM, vfiber, 'aspect', 25); C_ud_lt = lielens(C_fiber, CisoM, vfiber, 'aspect', 25);

eng2C / C2eng 负责工程常数与刚度矩阵的互换，iso2C 构造各向同性矩阵。得到单向刚度 C_ud 后，再结合 A2 和 A4 做方向平均：

A4k = closeA4(A2, 'hybrid'); Cavg = oravg(C_ud_mt, A2, A4k, 'method', 'exact'); [E1, E2, G12, nu12] = C2eng(Cavg);

方向平均的内部实现是关于 A2/A4 的线性组合，对应聚合物的三阶张量表示。值得检查两点：Cavg 是否保持预期的弹性对称性；C2eng 换算出的工程常数是否落在物理合理区间。若 E1 异常偏大，大多数原因是 A4 闭包与方向平均坐标系不一致，而不是平均算法本身有问题。

提示：先用极端取向验证。把 A2 设成完全单方向（对角线为 [1 0 0]），算出的 E1 应恢复为单向刚度值。若偏差超过 1%，去查 A2 与 A4 是否同源、坐标系是否对齐。

5.3 Halpin-Tsai 的工程定位与长径比敏感性

Halpin-Tsai 需要最少输入：纤维模量、基体模量、长径比、体积分数。即使没有完整实验数据，用它粗判注塑件在纤维分散时的刚度下限，足以把设计空间收敛一个量级：

E_ht = halpin(Ef, Em, vfiber, 25, 'short');

对于短纤维，纵向的 ξ=2L/d，横向的 ξ=2。长径比低于 10 时 Halpin-Tsai 与 Mori-Tanaka 接近；长径比增大后两者差 15% 以上，这是因为 Halpin-Tsai 的横向修正过于简单，无法反映纤维端部应力集中的高阶效应。所以我的建议很明确：粗筛用 Halpin-Tsai，出报告用 Mori-Tanaka 或 Lielens，至少算两个模型交叉验证。

6. Clayer2laminate 与 LaminateTheory.pdf：板级分析的验证技巧

用 C2eng 得到单一材料点的刚度之后，复合材料板在厚度方向存在取向梯度，需要层压板理论计算宏观弯曲与拉伸耦合。Clayer2laminate.m 实现经典层压板理论的 ABD 矩阵拼装：

Cz = cell(1, nLayer); for k = 1:nLayer A2k = A2_thickness{k}; % 第 k 层取向张量 Ck = oravg(C_ud_mt, A2k, closeA4(A2k, 'hybrid')); Cz{k} = Ck; end [ABD, N, M] = Clayer2laminate(Cz, hvec);

返回的 ABD 矩阵可直接换算面内工程常数与弯曲刚度。注意不要跳过 closeA4 直接传入 A2，层板矩阵需要每个材料点的四阶各向异性修正信息。最常见的错误是各层局部坐标旋转约定不一致，导致 A11 与 A22 的比值异常。

这一步最实用的验证技巧是自检：构造单层 [0] 板和 [90] 板，确认面内刚度 A11 的比值与单层 E1/E2 吻合。偏差超过 5% 时，检查坐标旋转矩阵的符号约定和 C2eng 的系数排列顺序，修正后重跑。

参数敏感性分析时，我习惯在厚度方向取 5 个代表点而非逐层遍历。先把这 5 个点的取向张量对 ABD 矩阵的贡献算出来，定位敏感层位，再回到完整分层重新计算。层数划分有个经验界限：层数小于取向梯度实际厚度尺度的 1/3 时，面内刚度偏差可达 8%；超过 12 层后 ABD 矩阵收敛进入平台，再加密只是增加计算时间。把 LaminateTheory.pdf 与 Clayer2laminate 的输出对照阅读，可以快速确认自己定义的各层方向与经典理论的约定是否一致，这是整套流程收口前最值得做的一次核对。

本文还有配套的精品资源，点击获取