聚酰亚胺(PI)具有独特的酰亚胺环结构，较好的耐热性、电气性能、力学性能及化学稳定性等优点，目前已被广泛应用于电气绝缘、微电子、水利电力、机械制造等领域[1-3]。但随着PI应用领域技术的快速发展，纯PI制品已逐渐不能满足某些功能化需求。为解决这一问题，相关领域的国内外学者集中对PI进行特定的功能化改性。而目前通过1,3-双（3-氨基丙基）-1,1,3,3-四甲基二硅氧烷(GAPD)对PI进行硅氧功能化改性的研究不多。丁孟贤[4]研究表明：相同含量下体系中引入的硅氧烷链过长，导致其酰亚胺链段过长，分相过于明显。因此，选用二硅氧烷进行硅氧功能化改性，可以得到几乎没有分相的含Si—O结构的PI。丁梓桉等[5]利用GAPD对PI分子结构进行改性以实现硅氧功能化。但宏观的实验研究不能明确解释GAPD对PI材料热力学性能的影响机制和原理。不同含量的GAPD对PI改性的分子动力学模拟的研究不充分，尤其在PI/GAPD体系的热力学性能及结构参数方面。相比传统的实验研究，分子模拟技术可以从微观层面模拟分子的结构和运动，计算体系的各项重要参数，分析其理化性质。在新材料的设计方面，不仅可以缩短新材料的研制周期，还能够降低研究成本[6]。林家齐等[7]利用多尺度建模方法构建PI/钽铌酸钾纳米颗粒复合物模型，通过分子动力学模拟，研究不同尺寸钽铌酸钾纳米颗粒对复合材料的结构、弹性模量和相互作用能的影响规律。雷环宇等[8]综述近年来分子模拟方法在PI中的应用，并基于分子模拟的方法研究PI的热力学性能及结构参数。本实验基于分子动力学原理，采用Material Studio(MS)软件搭建硅氧功能化PI分子模型，探究PI在不同含量的GAPD分子改性下其热力学性能参数及结构参数的变化，对比各个体系参数的变化规律，为研究PI在硅氧改性方面提供一定的参考。1反应原理及模型构建1.1反应原理图1为GAPD对PI改性反应机理。10.15925/j.cnki.issn1005-3360.2022.10.015.F001图1GAPD对PI改性反应机理Fig.1Reaction mechanism of GAPD modified PI物质的量n、m0，4,4-二氨基二苯醚(ODA)和GAPD作为反应物，将一定量的N,N-二甲基乙酰胺(DMA)作为反应催化剂，同时加入1.02(n+m)的均苯四甲酸酐(PMDA)参与反应，生成相应的PI酸(PAA)溶液。将PAA溶液进行脱水缩合可形成由GAPD分子改性的硅氧功能化PI薄膜。1.2模型构建为了对比研究不同含量硅氧基对PI性能的影响，考虑纳米复合材料改性后性能一般在较低含量时达到最优，选择的研究对象为工程实际常采用的均苯型PI和硅氧基接枝密度分别为5%和10%的GAPD改性PI。根据GAPD改性PI分子的反应机理，图2为PI、GAPD和GAPD改性PI的单体分子模型。图2中灰色为C原子、白色为H原子、红色为O原子、蓝色为N原子、黄色为Si原子。研究表明模拟过程中选取5~20个重复单元建立聚合物模型即可满足精度要求。综合考虑实际情况和计算机的模拟时长，最终确定每个体系为包含两条聚合度为20的聚合物单链的无定形区模型。以结构优化后的PI单体和接枝硅氧基的GAPD改性PI单体模型使用Amorphous cell(AC)模块分别建立纤维素聚合度(DP)为20、初始密度为0.1 g/cm3的PI模型，以及硅氧基接枝密度分别为5%、10%的GAPD改性PI模型，分别记为PI、PI/5%GAPD、PI/10%GAPD。图2PI、GAPD和GAPD改性PI的单体分子模型Fig.2Monomer molecular model of PI, GAPD and modification of PI by GAPD10.15925/j.cnki.issn1005-3360.2022.10.015.F2a1(a)capton型PI单体分子模型10.15925/j.cnki.issn1005-3360.2022.10.015.F2a2(b)GAPD单体分子模型10.15925/j.cnki.issn1005-3360.2022.10.015.F2a3(c)PI/GAPD改性单体分子模型2分子动力学模拟2.1弛豫处理对于已经构建好的模型，只有在其近似于实际材料的情况下才能进行下一步研究，因此需要对已在盒子中封装好的分子模型进行弛豫处理，使模型分子体系的整体密度处于稳定均衡。首先，选择Forcite/COMPASS力场对所构建的无定形区模型进行结构优化，最大迭代次数为5×104。对模型进行一系列分子动力学优化：先在等温等压(NPT)系综下以高压(0.5 GPa)进行10 ps的MD压缩模型体系；在正则(NVT)系综下，以598 K和298 K温度对模型进行热平衡处理，每个模拟时长20 ps；在常压(1.013×105 Pa)下采用NPT系综模拟100 ps释放体系并检查密度波动，若密度-时间曲线不稳定，可将模拟时间适当延长；最后在NVT系综下进行200 ps的动力学模拟对体系做进一步动力学平衡。经过以上处理，得到合理构象的无定形模型。其中，优化处理后的PI模型密度达到1.353 g/cm3，而实际未改性PI分子密度在1.26~1.36 g/cm3，说明构建的模型与实际数据吻合较好[9]，验证了所建立模型的可靠性。然后以相同方式得到PI/GAPD改性分子的无定形模型，图3为不同PI体系分子链模型。10.15925/j.cnki.issn1005-3360.2022.10.015.F003图3不同PI体系分子链模型Fig.3Molecular chain model of different PI systems图4为PI/5%GAPD体系的密度、能量随时间的变化。从图4可以看出，体系的密度在9 ps后稳定在1.26 g/cm3左右，而其各项能量也在几乎同一时刻由剧烈地波动渐变为稳定，表明此时改性体系各项指标已经达到平衡状态，满足进行下一步研究的条件。图4PI/5%GAPD体系的密度、能量随时间的变化Fig.4Change of density and energy of PI/5%GAPD system with time10.15925/j.cnki.issn1005-3360.2022.10.015.F4a110.15925/j.cnki.issn1005-3360.2022.10.015.F4a22.2退火优化为得到弛豫处理后模型密度及其体积随温度的变化规律，对体系的模型进行650~350 K的退火模拟，设置温差为25 K，降温速率控制为25 K/100 ps，体系在所有温度下依次在NVT、NPT系综下分别进行50 ps的分子动力学模拟，压强设定为1.013×105 Pa，每个过程均记录体积、密度等参数。在退火全过程中，力场选择COMPASS[10]，分别采用Nose控温[11]，Berendsen控压法[12]，范德华作用采用Atom Based方法[13]，静电作用采用Ewald法[14]。3热力学参数计算3.1玻璃化转化温度玻璃化转化温度(Tg)是聚合物一个重要的物理参量，表示高分子链段能够运动的临界温度。低于此温度聚合物呈玻璃态，分子链几乎不存在运动能力，而随着温度升高至Tg，聚合物开始由玻璃态向高弹态转变，温度继续升高，聚合则完全呈现高弹态，在高弹态时，分子间的相互作用力变小，分子链的运动能力提升，从而更容易发生整体运动。而分子链的运动能力越强，运动越剧烈，其热稳定性就越差。目前虽然可以通过各种实验测得聚合物的Tg，但缺少普遍适用的理论解释聚合物玻璃化转变的微观机理，软件模拟预测Tg是一种便捷的途径。预测Tg常用方法为密度-温度曲线法[15]，以50 K为一个节点计算相应体系平均密度，密度-温度的散点图中拐点两侧分别拟合一条直线，延长两条直线使其相交，产生交点的横坐标为体系Tg。图5为线性拟合密度-温度曲线。10.15925/j.cnki.issn1005-3360.2022.10.015.F005图5线性拟合密度-温度曲线Fig.5Linear fitting density-temperature curve刘文兴等[16]构建了均苯型PI，并通过分子动力学模拟测算出的Tg为606.5 K。从图5可以得出，本研究中纯PI的Tg模拟值为604.12 K，与文献值相近，说明本文构建的分子模型具有较高的可靠性。以相同方式测算得出PI/5%GAPD和PI/10%GAPD的Tg值并与纯PI组对比，图6为不同PI体系的Tg。从图6可以看出，PI/5%GAPD的Tg比纯PI提高了30.11 K，PI/10%GAPD的Tg比纯PI提高了22.42 K，说明不同浓度GAPD对PI的结构改性均能够在一定程度上提高PI材料的热学性能。但PI/10%GAPD体系的Tg没有继续提升，比PI/5%GAPD的Tg降低了7.69 K。10.15925/j.cnki.issn1005-3360.2022.10.015.F006图6不同PI体系的TgFig.6Tg of different PI systems3.2热膨胀系数热膨胀系数(CTE)能影响聚合物加工过程中产生的残余应力。CTE计算公式为：βCTE=13V0(∂V∂T)p （1）式（1）中：βCTE为热膨胀系数，K-1；V0为体系的初始体积，cm3；V为体系的体积，cm3；T为体系的温度，K。本研究从分子动力学模拟的过程中得到了纯PI体系和不同浓度的PI/GAPD的改性体系的体积与温度关系图，提取出其中的数据进行拟合，通过式（1）得到三种体系分别在玻璃态和橡胶态时的CTE，表1为不同体系的CTE。10.15925/j.cnki.issn1005-3360.2022.10.015.T001表1不同体系的CTETab.1CTE of different systems体系玻璃态CTE/10-5高弹态CTE/10-4纯PI2.291855.24985PI/5%GAPD1.698694.92817PI/10%GAPD2.033175.11916K-1K-1从表1可以看出，3种体系在玻璃态时的CTE全都小于其在橡胶态时的CTE，其原因是，当温度处于Tg以下时，整个体系处于玻璃态，其占有体积变化很小，而自由体积几乎不发生变化，故此时体系的CTE较小。而当温度升高到Tg或更高时，体系的自由体积也在温度的影响下而快速膨胀变化，此时体系体积的变化同时受到占有体积和自由体积的变化影响，因此CTE较大。利用GAPD对PI进行结构改性，可以减小PI体系的CTE，且GAPD含量为5%浓度的体系比10%浓度的体系效果更明显，这可能与GAPD分子本身较小的CTE有关。3.3弹性模量弹性模量为聚合物重要的力学性能参数，其计算方法通常采用静态常应变方法进行计算。即对于已经优化好的各体系分子模型施加一个应变，使材料模型产生一个较小的形变。此时分子模型将会处于一个新的状态，然后对其重新进行优化，使这个应变可以施加到不同的方向，以此就可以计算出各体系模型的刚度矩阵，计算公式为：Cij=1V×∂2U∂εi∂εj=∂σi∂εi∂εj=σ+-σ-2εj （2）式（2）中：V为体系总体积，cm3；ε为体系应变；i和j表示应变方向，通常取值均为1、2和3；σ为应力，σ+表示正应力，说明体系处于拉伸状态；σ-表示负应力，说明体系处于压缩状态。通过式（2）计算3种体系模型的刚度矩阵，发现其均可视为各向同性材料，可以通过刚度矩阵计算出弹性常数λ和μ，计算公式为：λ=13C11+C22+C33-23C44+C55+C66μ=13C44+C55+C66 （3）而体系的弹性模量E、剪切模量G、体积模量B的计算公式为：E=μ×3λ+μλ+μG=μB=λ+23μ （4）图7为不同PI体系的力学性能。从图7可以看出，PI/GAPD体系的力学性能明显要优于纯PI体系。PI/5%GAPD体系比纯PI体系弹性模量提高16.1%，剪切模量提高7.25%，体积模量提高4.34%。由于GAPD分子中所含有的硅氧键本身就拥有良好的力学性能。且硅氧键具有一定的柔顺性，过高含量的引入导致改性后的材料提升不明显，低含量的GAPD改性体系的力学性能更优异。10.15925/j.cnki.issn1005-3360.2022.10.015.F007图7不同PI体系的力学性能Fig.7Mechanical properties of different PI system4结构参数计算4.1径向分布函数径向分布函数(RDF)通常用g(r)表示，其可以反映材料微观结构的特征物理量，表示在距离某个原子某个值的地方出现另一个其他原子的概率密度相对于随机分布概率密度的比值，可以反映非键合原子间的相互作用[17]。gAB=nBV4πr2drNB （5）式（5）中：nB为A原子以r为半径的范围内B原子的数量；NB为体系中B原子的总数；V为总体积，nm3。图8为不同PI体系的全原子的RDF。从图8可以看出，对比不同含量GAPD对PI改性前后体系的RDF曲线，发现3种体系的RDF在整体上非常相似，并无明显差异。但在1.5~2 nm处，PI/5%GAPD以及PI/10%GAPD体系均出现细微的波动，但纯PI体系在此区间没有出现，且PI/10%GAPD的波动略微强于PI/5%GAPD体系，这可能是由于GAPD改性子中的Si原子引起的，且随着GAPD含量的增加，整个体系的Si原子数量增多，导致波动更加显著。10.15925/j.cnki.issn1005-3360.2022.10.015.F008图8不同PI体系的全原子的RDFFig.8Whole atom RDF of different PI systems4.2内聚能密度影响PI热力学性能的重要参数之一为其体系中分子链段的相互作用力大小，若其相互作用力越小，其体系的Tg就越低，结构就越松散，弹性模量也就越小。而分子链段的相互作用力的大小通常使用内聚能密度ρCED表示，值越大，其所包含的内聚能越大，则其相互作用力越大，其计算公式为：ρCED=EcohV=EinterV=Etotal-EintraV （6）式（6）中：Ecoh为分子链段的内聚能，J；V为体积，cm3；Einter为不同分子链段间的能量，J；Etotal为总能量，J；Eintra为分子链段内的能量，J；“ ”表示取在NVT或NPT系综下的平均值。温度为300 K时，对三种体系进行ρCED的计算，表2为计算结果。10.15925/j.cnki.issn1005-3360.2022.10.015.T002表2不同PI体系的内聚能密度Tab.2Cohesive energy density of different PI systems体系内聚能密度纯PI400PI/5%GAPD456PI/10%GAPD413J‧cm-3J‧cm-3从表2可以看出，与纯PI体系相比，使用GAPD分子改性的硅氧功能化PI体系的ρCED都有所提高，并且GAPD浓度为5%的改性体系的ρCED最大，对PI材料的ρCED提升最为明显。所得的ρCED模拟计算结果与前文中的Tg以及弹性模量的结果相似整体，所以可以得出，引入GAPD对PI改性后，体系的ρCED与其Tg和弹性模量成正相关。4.3自由体积占比自由体积理论最早由Fox和Flory提出[18]，聚合物材料中的体积VT由两部分组成：一部分是被分子占据的体积V0（即实际占有体积），另一部分是未被占据的体积Vf（即自由体积），即:VT=V0+Vf （7）自由体积通常指体系分子之间的细微间隙，在实际材料中常以“孔穴”的形式存在。由于自由体积的存在使得体系分子拥有了部分的活动空间，所以整个体系分子具有了活动的可能性。本研究采用硬球探针法计算体系模型的自由体积参数，即利用半径为1.0 nm的球形探针探测聚合物模型[19]。在没有统一各体系总体积的情况下，各体系的自由体积通常不具有可比性，而自由体积占比(FFV)是以各体系自身总体积为基准，反映了各体系自由体积的相对大小，具有可比性，所以将其引入。FFV的表达式为：φFFV=VfV0+Vf×100% （8）式（8）中：Vf为自由体积，nm3；V0为实际占有体积，nm3。对3种体系加以计算，表3为不同PI体系的FFV。从表3可以看出，加入GAPD分子改性后，所得PI体系的自由体积均降低，这是由于GAPD改性分子的引入，阻碍了分子链段的运动及堆砌，导致分子链段运动能力下降，故复合体系的自由体积占比降低。同时，PI/GAPD复合体系的自由体积占比随着GAPD改性分子的含量增大而增大，这可能是因为GAPD改性分子与体系中其他分子间的相互作用力更强，进而限制了分子链段的运动。10.15925/j.cnki.issn1005-3360.2022.10.015.T003表3不同PI体系的FFVTab.3FFV of different PI systems体系V0/nm3Vf/nm3φFFV/%纯PI15.839934.4306921.86PI/5%GAPD15.805583.7159219.04PI/10%GAPD15.352234.0343620.814.4均方位移为了探究GAPD对PI链运动能力的影响，选用均方位移(MSD)表征不同体系下PI链热运动的位移，MSD的值越大，代表物体热运动越剧烈。有N个原子的体系中的MSD的计算公式为：MSD=ri(t)-ri(0)2 （9）式（9）中：ri(0)为零时刻i原子位置的向量；ri(t)为在t时刻i原子位置向量。分别在298、343、363 K温度下对3种体系进行100 ps的NPT分子动力学模拟，利用forcite模块计算得到各体系的MSD值。图9为298 K下不同PI体系的MSD。从图9可以看出，随着时间的增加，3种体系的MSD均呈现上升趋势。GAPD改性体系比纯PI的MSD小，PI/5%GAPD体系MSD最小。说明引入的改性分子对体系的运动能力具有一定的减弱作用。10.15925/j.cnki.issn1005-3360.2022.10.015.F009图9298 K下不同PI体系的MSDFig.9MSD of different PI system at 298 K5结论（1）通过使用GAPD对PI进行结构改性，可以有效提高PI材料的Tg和弹性模量，降低其热膨胀系数，在一定程度上增强了PI的热稳定性以及力学性能。（2）硅氧功能化PI的热力学性能并不会随着GAPD含量的增大而一直提高，当参与反应的GAPD浓度为5%时，对PI材料在热力学性能方面改性提升更加有效。相比纯PI，PI/5%GAPD的Tg提高30.11 K，弹性模量提高16.1%，剪切模量提高7.25%，体积模量提高4.34%。（3）引入GAPD对PI进行结构改性改变了PI材料的微观结构。对比纯PI组，其提高了内聚能密度，略微降低了自由体积占比以及均方位移。5%浓度的GAPD对PI材料微观结构的影响更加显著。而在材料全原子的径向分布函数方面，利用GAPD对PI改性并没有产生很大影响。