摘要
针对传统单一尺度优化范式在填料拓扑调控机制、分子有序度耦合效应及跨尺度协同作用方面的理论缺失，通过“填料拓扑设计-分子有序排列-密度协调”多级协同策略，研究填料拓扑设计与树脂基体改性对环氧树脂基复合材料导热性能的协同调控机制。基于有限元均匀化理论构建纤维、片状及椭球体模型，揭示填料体积分数、长径比及空间取向对传热网络的调控规律。结合反向非平衡分子动力学（RNEMD）方法，提出通过分子有序排列与密度诱导非键相互作用协同提升树脂本征热导率的新策略。结果表明:片状填料（体积分数25%）的有效导热系数(ETC)达4.2 W/(m·K)；有序交联结构与密度诱导非键相互作用增强（非键能占比60%）使树脂本征热导率从0.3 W/(m·K)提升至2.85 W/(m·K)。多尺度协同分析显示，掺杂有序交联体系在25%填充与密度增强耦合作用下，ETC突破6.8 W/(m·K)，密度协同效应（ρ＞1.5 g/cm3）使热流密度标准差降低64%。研究建立的跨尺度模拟框架揭示了传热网络构建与分子有序排列的定量关联，为设计高导热复合材料提供了理论新范式。
Abstract
This research addresses the theoretical limitations of conventional single-scale optimization approaches in understanding filler topology regulation mechanisms, molecular alignment coupling effects, and multi-scale cooperative interactions. Through a multi-level collaborative strategy of "filler topology design-molecular ordering-density coordination", we systematically investigate the synergistic control mechanism between filler architectural design and resin matrix modification on the thermal conductive performance of epoxy-based composite materials. Fiber, sheet, and ellipsoid models are constructed based on finite element homogenization theory, revealing the regulation rules of filler volume fraction, aspect ratio, and spatial orientation on the heat transfer network. In conjunction with the reverse non-equilibrium molecular dynamics (RNEMD) method, we propose a novel strategy to enhance the intrinsic thermal conductivity of the resin through the synergistic effects of molecular ordering and density-induced non-bonded interactions. The study shows that the flake filler with 25% volume fraction exhibits an effective thermal conductivity (ETC) of 4.2 W/(m·K). The ordered cross-linked structure and enhanced density-induced non-bond interactions (where non-bonded energy accounts for 60%) increase the intrinsic thermal conductivity of the resin from 0.3 W/(m·K) to 2.85 W/(m·K). Multi-scale synergistic analysis shows that under the coupled effects of 25% filling and density enhancement, the ETC of the doped ordered cross-linked system exceeds 6.8 W/(m·K). The density synergistic effect (ρ＞1.5 g/cm3) reduces the standard deviation of heat flux density by 64%. The cross-scale simulation framework established in the study reveals the quantitative correlation between heat transfer network construction and molecular ordering, providing a new theoretical paradigm for designing high thermal conductivity composites.
随着5G通讯、航空航天等领域对高功率密度电子器件散热需求的提升，开发兼具高导热与优异加工性能的环氧树脂（EP）基复合材料成为研究热点[1-3]。环氧树脂因其优异的机械性能、化学稳定性和加工工艺性，在电子封装、航空航天等高新技术广泛应用[4-7]。然而，其本征热导率与快速发展的热管理需求间存在显著鸿沟，严重制约了高功率密度器件的可靠性[8-10]。
传统研究范式多聚焦于单一尺度优化，在填料拓扑结构对热输运网络的调控机制方面，尚未建立具有普适性的理论模型体系[23]。树脂基体的本征热导率提升受制于分子链无序排列导致的声子散射效应及弱非键相互作用的能量传递瓶颈[19-20]。传统环氧复合材料的热导率已难以满足散热需求，亟需通过同优化实现热输运性能的突破。
当前研究面临三重挑战：1）填料拓扑结构（如长径比、取向度）对热输运网络形成的影响机制尚未建立定量预测模型；2）树脂交联网络与分子链有序度的耦合效应对本征热导率的调控规律不明；3）宏观填料分布与微观分子排列的跨尺度协同作用缺乏系统研究。
本研究提出“填料拓扑优化-分子有序排列”双向协同调控新范式，通过多尺度数值模拟揭示导热增强机制。在介观尺度，构建纤维/片层/椭球体填料的代表性体积单元（RVE）模型，解析体积分数、长径比与空间取向对热流分布的影响规律；在分子尺度，采用反非平衡分子动力学（RNEMD）量化交联度、分子链有序度与非键相互作用对传输的调控机制。引入密度诱导非键能重构，突破填料-基体界面热阻限制。
1 多尺度建模方法
1.1 填料拓扑建模
研究主要针对环氧树脂本征结构的调整以及改变掺杂物结构实现导热性能的提升。椭球粒子必须满足曲面约束方程[11]：
(1)
其中，椭球体的球心位于（0，0，h）处，（x，y，z）是直线上的点。用旋转椭球坐标系（η，ξ，χ）中径向坐标值ξ表示粒子边界面的形状特征，在长旋转椭球坐标系中，L和D为旋转椭球面的长轴和短轴长。当L=D时，为球体。
生成填充物结构，如图1根据填充物设定的方向进行旋转随机2D或随机3D随机方向，依次生成其余结构。研究结构均选用碳纤维束，制作为片状填料，椭球填料等，热导率均呈现各向异性，其轴向热导率为60 W/（m·K），径向热导率为4 W/（m·K）。
图1椭球颗粒结构示意
Fig.1Schematic diagram of ellipsoidal particle structure
1.2 树脂基体分子动力学模型
如图2所示，利用分子动力学模拟研究了不同分子排列制备的交联环氧树脂的热传输性能，并构建不同分子链的结构模型。研究体系选用双酚A二缩水甘油醚（DGEBA）作为环氧树脂基体，搭配4-甲基环己烷二胺（HTDA）作为固化剂，构建分子动力学模型。
图2分子结构优化示意
Fig.2Schematic diagram of molecular structure optimization
交联过程通过识别DGEBA活性碳和HTDA氮原子建立反应位点，初始设定0.35 nm键合距离阈值、循环执行距离计算、原子成键、结构优化（600 K-NVT和1 atm-NPT）等步骤。每轮循环扩展0.05 nm键合阈值，最高至1.2 nm，直至达到目标交联度或最大阈值时终止。通过分子动力学模拟实现交联网络构建，每次迭代包含5次结构松弛，确保体系能量稳定。
1.3 热输运控制方程
通过求解以下拉普拉斯方程获得RVE的稳态温度场，顶部和底部表面分别设置为特定温度Ttop和Tbottom，其他边界是绝热的。同样，沿x和y方向的有效导热系数（ETC）计算为沿z方向的ETC。热输运路径如图3所示。对于非均质材料，边界条件和其余边界绝热被广泛采用。各向异性材料的导热微分方程为
(2)
本研究使用的树脂本征热导率为0.3 W/（m·K），边界条件如下：
(3)
在树脂材料与填料材料边界处，必须注意的是，假定填料和树脂基体之间存在完美的结合条件，不同材料界面处应满足以下条件：
(4)
图3复合材料传热示意
Fig.3Schematic diagram of heat transfer in composites
模拟计算过程所使用的PCFF力场是在CFF91力场参数上的拓展，PCFF力场参数已经成功地用于研究各种聚合物的力学性能，如聚乙烯、环氧树脂[12]，函数形式中的作用参数由从头算获得[13-14]。PCFF力场的精确公式[15-17]为
(5)
式中：r0为平衡键长，θ0为平衡键角，nτ与φn控制旋转周期与角轴变化，最后两项表示分子内不成键的原子间的范德华作用势和静电作用势，环氧树脂基质之间的范德华相互作用势即9-6 Lennard-Jones（L-J）电势。其中，rij为原子i和j之间的距离，ε和σ分别为能量和距离常数。
非键作用采用1 nm截断半径，同时运用尾端修正方法确保长程作用计算的完备性。模拟过程中同时考虑范德瓦尔斯与静电相互作用。使用大规模原子/分子质量并行模拟器（LAMMPS）[18]进行所有分子动力学模拟，并使用速度Verlet方法对运动方程进行积分。
为了验证分子模型和PCFF势函数，采用RNEMD方法计算了纯无定形分子随机混合交联环氧树脂在300 K下的热导率。如图4所示，在RNEMD模型中热源置于中心区域，热沉对称分布于系统两端。模型包含12 300个原子，模拟尺寸为3 nm×3 nm×14 nm。如图5所示，通过RNEMD方法，计算得出纯环氧树脂的热导率为（0.25±0.001）W/（m·K），这一结果与文献[19]中环氧树脂在室温下的热导率值约为0.25 W/（m·K）相符。此外，该结果与文献[20]中实验结果在300~340 K的平均值一致，均为约0.25 W/（m·K）。
图4反向非平衡分子动力学热导率计算方法示意
Fig.4Schematic diagram of reverse non-equilibrium molecular dynamics thermal conductivity calculation method
为了验证掺杂复合材料的数值计算方法，验证了填充球形和片状增强复合材料计算模型准确性，如图6所示。
Maxwell-Eucken模型[21]:
(6)
Hatta模型[22]:
(7)
其中，沿平面方向的计算时，S=πL/（4X），L和X分别为薄片填料的有效直径和厚度。沿厚度方向计算时，S=1-πL/（4X）。球形和片状填料的计算结果与理论值进行了对比。结果表明，本文计算结果与现有实验结果一致（如图6所示）。
2 结果与讨论
2.1 填料拓扑结构影响
填充物结构形状是影响复合材料ETC变化的因素之一[23]。本研究以具有形态演化特征的填充颗粒体系为对象，重点考察颗粒几何构型对复合材料热传导性能的调控机制。
结果表明，椭球形填充体在基体中的空间取向分布对其ETC有显著影响。当颗粒长轴与热流传输方向共线排列时，热导率随填充比例提升呈超线性增长特征，在0~12%体积分数内最大增幅达370%，较随机分布体系具有更优的热传导网络构建能力（详见图7）。
对比分析不同维度分布的椭球颗粒增强体系发现，在低填充量阶段（＜5%体积分数），轴向有序排列体系的ETC增幅斜率较二维随机分布体系提升约118%，较三维随机分布体系增幅达219%。值得注意的是，三维随机体系中颗粒长径比对ETC的调控敏感度较低，其归一化有效导热系数波动范围仅为0~170%。当填充量低于临界值（约8%）时，颗粒几何形貌对ETC的贡献度不足5%；但当填充量超过15%时，几何结构参数的调控效能与填充比例呈等效耦合特性。该研究为复合材料设计提供了重要启示，即通过优化颗粒空间拓扑构型可替代传统单纯提升填充量的方法。
图7归一化有效导热系数随长宽比的变化规律
Fig.7Normalized effective thermal conductivity variation with aspect ratio
如图8所示，空间内随机分布掺杂形成的复合材料在热输运效率方面表现出明显的优势。尽管热输运效率较高，整体ETC却低于平面内随机分布掺杂的情况。这种现象主要归因于空间内随机分布导致沿热流方向的总界面面积增大，从而增加了界面热阻。
图8长宽比为5的椭球体不同分布复合材料截面热图
Fig.8Thermal map of composite material cross-sections with ellipsoidal bodies with an aspect ratio of 5 in different distributions
如图9所示，所有椭球体的分布方向一致，且其长度方向被设定为温度梯度方向。结果表明，随着填料体积分数的增加，不同结构填料复合材料的归一化有效导热系数呈现不同程度的线性增加。当填料为不连续颗粒形状时，低体积分数（小于2%）下，球体、椭球体、二十面体和薄片的归一化有效导热系数存在微小差异，对热输运的影响不明显。随着体积分数的增加，这些微小差异逐渐显现，椭球体在传热方向的归一化有效导热系数增长率较大，且随体积分数增加，增长趋势更加明显。
图9填料结构归一化有效导热系数与体积分数关系
Fig.9Relationship between normalized effective thermal conductivity and volume fraction of filler structure
如图10所示，由于多片薄片在空间随机分布所形成的分布组合与球体相似，表现出与球体变化规律高度一致。对于特定方向分布的薄片增强复合材料，其归一化有效导热系数影响几乎可以忽略，因该方向为薄片厚度方向，无法形成有效的热网络。对于颗粒增强或片状增强复合材料，其ETC在体积分数处达到峰值，且此时界面热阻成为主导因素。
另一方面，在低体积分数下，相互融合的填充物所占比例小，对ETC贡献有限；在高体积分数下，非融合和融合填充物间距小，融合效果减弱，呈现一致变化趋势。相互融合的椭球体在xoy平面随机分布时，长度方向沿热流方向的分布概率大幅增加，更易形成有效热传输网络，提高ETC。归一化有效导热系数随着体积分数增加逐步提高，当体积分数超过0.03时，增强效果显著。随着颗粒数量增加，相互融合概率增大，热网络形成可能性更高。如图11所示，纤维相互融合界面呈现为连续截面，拓宽热传输通道，促进热流传递。
图10粘连与无粘连结构复合材料归一化有效导热系数的变化
Fig.10Variation diagram of normalized effective thermal conductivity of adhesive and non-adhesive composite materials
图11交叉融合纤维复合材料截面热图
Fig.11Thermal map of cross-mixed fiber composite sections
2.2 树脂基体改性机制
如图12所示，交联度对环氧树脂基复合材料的热传导性能有显著影响。低交联度（5%~20%）时，温度梯度波动大，热传导不均匀，热流路径随机。当交联度增至50%时，温度梯度趋于稳定，正向与反向梯度差异减小，表明交联网络形成降低了分子链自由度，减少声子散射，提高热传导效率。高交联度（50%以上）样本热流密度更稳定，表现出优异热管理性能。但交联度提高至70%时，热导率增益效果有限，可能因材料刚性增加限制热传导路径优化。交联度约50%时，材料表现出最佳热传导性能，既保证热流密度稳定性，又有效降低温度梯度波动，为优化环氧树脂基复合材料热管理性能提供理论依据。
图12密度1.3 g/cm3阶梯重叠结构温度梯度与热流密度的变化
Fig.12Variation of temperature gradient and heat flux density in a stepped overlapping structure with a density of 1.3 g/cm3
环氧树脂的热流密度随交联度变化呈非线性特征。低交联度（5%）材料因分子链松散，热传导依赖长链段热振动传递，热流密度低。交联度提升至20%~50%时，三维网络结构使分子链节段振动路径的均方根偏差（RMSD）降低28%，热流密度显著增加。但交联度超过60%后，分子链刚性增强引发空间位阻效应，振动能量耗散增加，热流密度增速趋缓甚至下降。分子动力学模拟揭示，低交联度（0~20%）下，热导率动态波动幅度达30%，传导路径不连续；中等交联度（40%~60%）时，刚柔平衡使热导率稳定性提升。高交联度（70%）材料分子链运动受限，温度梯度震荡幅度减小15%，热传导路径连续性受损。模拟结果表明：交联度在50%~60%时，材料维持适度柔性并缩短热传导路径，实现最优热导率。该平衡机制为环氧树脂热管理材料性能优化提供理论依据，热流密度峰值对应的临界交联度为工程应用的重要参数。
环氧树脂基复合材料的热导率受分子链结构和交联度的显著影响。如图13所示，中等交联度（50%或60%）的材料通常具有较高的热导率，而过高或过低的交联度会降低热导率。这是因为适度交联保持了分子链的有序排列，避免了过度交联导致的分子链刚性化。过度刚性限制了声子（热传递载体）传播，降低热导率。分子链结构也对热导率有影响。有序状分子链的热导率通常高于阶梯状分子链，因为有序排列减少了分子间散射，提升了声子的平均自由程，从而增强了热传导。
图13密度为1.3 g/cm3与1.5 g/cm3的阶梯重叠结构热导率计算结果
Fig.13Calculation results of thermal conductivity of stepped overlapping structures with a density of 1.3 g/cm3 and with a density of 1.5 g/cm3
径向分布函数g（r）反映了分子间相对位置的分布。径向分布函数计算基于环氧树脂交联网络中的特征原子对相互作用，具体选取双酚A二缩水甘油醚（DGEBA）环氧基中的氧原子与固化剂4-甲基环己烷二胺胺基的氢原子作为计算对象。该原子对的选择源于交联反应过程中氧与氮-氢的共价键形成机制。如图14所示，随着交联度增加，g（r）的峰值逐渐降低，表明分子流动性受限，结构变得更刚性。在低交联度（5%和20%）下，g（r）显示较高的峰值，表明分子结构更松散，具有较高的流动性。交联度为70%时，g（r）趋于平坦，r=0.121 nm附近急剧下降，r=0.225 nm后趋零，表明高度受限的网络中单体间距离分布均匀，分子运动受限。
图14不同交联度形成的有序结构的径向分布函数计算值
Fig.14Calculated radial distribution functions of long straight structures formed by different degrees of cross-linking
低交联度（5%）下，首峰（r=0.107 nm）g（r）值高达5.13，显示显著近程有序性；交联度增至20%时，首峰强度降至4.77，但第2配位层（r=0.205 nm）g（r）值升至0.86，表明交联网络增强中程协同作用。交联度50%时，首峰强度回升至5.35，第3配位层 （r=0.305 nm）g（r）值达1.03，显示分子链形成稳定堆积结构。高交联度（70%）时，首峰展宽，中程g（r）值降至1.73，表明刚性增大限制局部结构调整。适中交联度（50%~60%）下，三维网络兼具短程有序和中程协同作用，其优化的g（r）分布显著提升热振动传递效率。
2.3 密度协同效应
图15显示，不同交联度和密度条件下，环氧树脂的温度梯度随时间变化存在显著差异。在密度为1.5 g/cm3时，随着交联度从5%增加至70%，温度梯度的震荡幅度逐渐减小，表明高交联度有效抑制了热传导过程中的局部温度波动。低交联度（如5%）样品的温度梯度变化幅度较大，显示热传导过程中存在显著的局部温度不均匀性。随着交联度增加，环氧树脂内部逐渐形成三维网络结构，使热传导路径更加均匀，减少了局部热流密度波动。密度对温度梯度的影响也显著。在密度为1.3 g/cm3的体系中，温度梯度变化较为剧烈，表明分子间热耦合较弱，热传导不均匀性较大。而在密度为1.5 g/cm3的体系中，温度梯度变化较平稳，验证了高密度能够增强分子间热耦合效果，提高热传导的整体均匀性和稳定性。
图15密度为1.5 g/cm3有序结构温度与温度梯度计算结果
Fig.15Temperature and temperature gradient calculation results of long straight structure with a density of 1.5 g/cm3
密度的变化对材料热导率具有显著影响。如图16所示，在密度为1.5 g/cm3时，整体热导率值普遍高于密度为1.3 g/cm3的情况。表明材料密度的增加能够有效提高热传导性能，是由于分子间的紧密堆积减少了热传导中的空隙，从而增强了热量的传递效率。这种现象对于设计高密度、高导热性的复合材料具有重要意义。
图16密度为1.5 g/cm3不同交联度有序结构热导率计算结果
Fig.16Calculation results of thermal conductivity of ordered structures with different crosslinking degrees with density of 1.5 g/cm3
当ρ＞1.5 g/cm3时，能量分配呈现新的特征模式。键能占比增速提高至ΔEbond/Δρ=0.30，而非键能增速保持 ΔEnon-bond/Δρ=0.40。如图17所示，特别在ρ=2.0 g/cm3时，非键能占比达61.4%，远超键能占比20.8%。高密度下分子间距缩短导致非键相互作用网络形成连续通路，而化学键的短程特性使其空间分布呈现局域化特征。在ρ＜1 g/cm3区间遵循 Enon-bond∝ρ2.3的幂律关系，而在ρ＞1.5 g/cm3区间转变为Ebond∝exp（0.7ρ）的指数增长模式。这种行为表明材料可能经历从软物质特性向硬质材料的结构转变。
图17不同密度结构键能、非键能、总能量的数值
Fig.17Numerical values of bond energy, non-bond energy, and total energy for structures with different densities
3 结论
1）填料的拓扑结构对环氧树脂基复合材料的导热性能具有显著影响。片状填料在25%填充率下，复合材料的有效导热系数可达4.2 W/（m·K）。而椭球填料采用空间随机分布模式，通过多向热桥效应，有效降低了热通量的波动幅度，降幅高达42%。合理设计填料的形貌和分布方式能够优化热传导路径。
2）在树脂基体改性方面，有序交联结构能够显著提升分子链的有序排列程度，提升幅度达215%。同时，密度增强诱导的非键相互作用占比高达60%，使得树脂基体的本征热导率从0.3 W/（m·K）提升至2.85 W/（m·K）。径向分布函数分析进一步证实，适度交联（50%~60%）能够在r=0.107 nm处维持g（r）=5.35的短程有序性，并保留中程（r=0.205 nm）协同作用的动态调控能力，实现分子层面的热传导优化。
3）片状BN填料与有序交联树脂基体的协同作用进一步提升了复合材料的导热性能。在25%填充率和密度协同作用下，ETC突破6.8 W/（m·K）。密度协同分析表明，当材料密度ρ＞1.5 g/cm3时，非键能贡献占比达61.4%，实现了传热路径的宏观与微观贯通，显著增强了复合材料的整体导热能力。
