HTPB推进剂中增塑剂扩散系数计算①

李红霞，强洪夫，李新其，王洪伟

(1.第二炮兵工程大学，西安 710025;2.中国航天科工集团公司六院四十六所，呼和浩特 010076)

HTPB推进剂中增塑剂扩散系数计算①

李红霞1，强洪夫1，李新其1，王洪伟2

(1.第二炮兵工程大学，西安 710025;2.中国航天科工集团公司六院四十六所，呼和浩特 010076)

针对丁羟推进剂/衬层界面增塑剂的迁移问题，研究增塑剂的扩散特性。采用分子动力学方法模拟增塑剂在丁羟粘合剂体系中的运动，再通过爱因斯坦关系式求得扩散系数;采用对衬层加厚的粘接试件进行加速老化实验方法，再通过费克第二定律计算得到增塑剂的扩散系数;考察不同环境温度和不同增塑剂含量条件下癸二酸二辛酯的扩散系数的变化。分析认为，增塑剂在丁羟粘合体系高分子链段运动产生的间隙中，以“跳跃”方式发生空间位置迁移。结果表明，扩散系数模拟值和实验值基本一致，数量级为10-12m2/s;温度升高，增塑剂运动活性加强，有效活动空间增大，扩散系数增大;粘合体系与增塑剂的共容限量使得增塑剂含量大于3%时，模拟得到的扩散系数依次略有下降。分子动力学方法计算增塑剂扩散系数更具优势。

丁羟推进剂;增塑剂;扩散系数;实验;模拟

0 引言

目前，对于固体推进剂组分扩散迁移研究，多以现象分析和实验为主。如徐思羽［1］分析了二茂铁衍生物在固体推进剂中的迁移问题;尹华丽［2-3］用加速老化试验研究了NEPE推进剂界面迁移组分含量的变化，还用浸泡增重法和气相色谱仪考察了HTPB衬层、绝热层对DOS和T27的吸收以及DOS、T27和GFP的迁移，他们的结果对后续研究工作具有重要参考价值。由于推进剂成分多、贮存时间长、反应复杂，实验中简化研究对象、设置极端条件，往往只能得到定性的结果，很少揭示其迁移机理。而分子模拟方法既能省时省力省钱达到实验的目的，还能提供物质的微观信息，为机理分析提供理论依据，现已成功用于推进剂领域［4-5］。笔者曾尝试过推进剂中增塑剂扩散行为的分子模拟工作［6-7］，尽管当时缺少扩散系数的试验验证，但分子模拟方法在其他领域的成功应用，展示出了该方法的强大优势。

扩散系数是物质的一种很重要的输运系数，很难用实验方法精确测量和确定，已有很多学者通过分子模拟手段来计算扩散系数或研究扩散过程，Kotelyanskii［8］等采用分子动力学模拟方法计算了水和盐在反渗透膜中的扩散过程和扩散系数，对膜的截盐效果做了微观的解释，Hofmann［9］等模拟了 H2、O2在无定形聚酰胺等几种高分子聚合物膜中的扩散过程，并计算了扩散系数。扩散系数的模拟结果得到了实验验证［10-12］，如乙醇在聚硅氧烷(PDMS)中扩散系数的分子模拟结果4.6×10-6cm2/s与实验结果4.5×10-6cm2/s接近;用分子动力学模拟得到的甲基丙烯酸甲酯寡聚物在高分子量的聚甲基丙烯酸甲酯母体中扩散系数的模拟结果，与实验得到的数值在趋势上吻合得很好。

增塑剂的扩散特性直接影响丁羟推进剂/衬层界面增塑剂的迁移速度，其迁移速度主要取决于增塑剂在推进剂的含量和扩散系数。本文主要开展HTPB推进剂中增塑剂扩散的分子模拟研究和扩散系数的实验测定工作，旨在分析扩散机理，以期为界面增塑剂迁移问题的解决提供依据。

1 增塑剂扩散系数的分子模拟计算

1.1 模拟方法

分子模拟的总体趋势是向着基于量子理论的方向发展，但流体的扩散性质研究一般不涉及到电子的转移，采用经验的位能函数就能达到精度要求，故不需采用费时的量子力学方法。相对于Monte Carlo模拟，分子动力学方法可更容易得到扩散性质，本文采用基于分子力学的分子动力学方法。

分子动力学模拟与真实实验方法具有相似的过程，首先要选择一个模拟体系，即样品的准备，然后在一定的模拟条件(如，恒温、恒压等)下未解一系列运动方程，等待体系平衡后，开始测定物性性质及进行数据分析。分子动力学模拟是一种“计算实验方法”，用来计算经典多体体系的平衡和传递性质。它通过求解体系中粒子的运动方程获得粒子的运动速度和运动轨迹，并按照统计物理和统计热力学原理得出系统相应的宏观物理特性。在分子动力学模拟中，扩散性质通过时间平均获得，即通过模拟中的均方位移MSD随时间变化的Einstein关系［13］，求得扩散系数D。

式中r(t)为t时刻分子的坐标;r(0)为初始坐标;〈│r(t)-r(0)│2〉为均方位移。

扩散系数既是材料的属性，又是温度的函数，本文将用分子动力学方法计算不同环境温度和增塑剂含量条件下的增塑剂癸二酸二辛酯(DOS)扩散系数。

1.2 模拟内容

整个模拟工作均在分子模拟软件Materials Studio 4.3上进行，全部计算在高性能计算平台IBM Blade上完成。用Visualizer模块构建DOS、IPDI、HTPB分子模型，选用COMPASS力场，用分子力学法进行几何结构优化。用Amorphous Cell模块构造由10个DOS、5个IPDI、5个HTPB分子组成的混合物(共3 566个原子)模型，其体系密度 ρ=0.91 g/cm3，胞元模型的a、b、c都是3.29 nm。对混合物结构进行几何优化，对优化过的结构进行弛豫。优化方法采用Smart Minimization方法，分别用atom-based和Ewald方法求范德华作用和静电作用。模拟过程中位能采用球形截断法，截断半径 cutoff取0.95 nm，spline width 取 0.1 nm，缓冲宽度buffer width取0.05 nm，时间步长为1 fs，截断距离之外的分子间作用按平均密度近似方法进行校正。在真空条件下，采用Andersen控温方法和Berendsen控压方法，各分子起始速度按Maxwell-Boltzmann分布取样，MD模拟积分方法为Velocity Verlet算法。

在Discover模块下，选择恒温恒压系综(NPT)，在1标准大气压下，增塑剂含量为3%时，分别在273、298、310、323、348 K下进行MD模拟。达到平衡后，在相同条件下再进行1 000 ps的MD模拟，保存全运动轨迹，分析不同温度下DOS在混合体系中的均方位移。根据增塑剂在推进剂中的含量为3%、6%、9%、15%，温度为298 K，构建由DOS、IPDI、HTPB 不同分子数组成的混合物模型。用同样方法分析不同含量的DOS在混合体系中的均方位移。图1是温度298 K、增塑剂3%含量下DOS的MSD-t曲线(略去其他MSD-t曲线)。

1.3 扩散系数的模拟结果

由均方位移随时间的变化曲线，经拟合得到各直线斜率，计算得到 273、298、310、323、348 K 下 DOS 在粘接体系中的扩散系数分别为1.67×10-12、3.33×10-12、4.17 × 10-12、5.17 × 10-12、7.17 × 10-12m2/s，DOS在粘接体系中的扩散系数随温度的升高而逐渐增大。计算得到增塑剂在推进剂中的含量为3%、6%、9%、15%时的扩散系数分别为 3.33 ×10-12、3.00 ×10-12、2.33 × 10-12、2.00 × 10-12m2/s，随增塑剂含量增大，扩散系数反而下降。

图1 增塑剂DOS的MSD-t曲线(298 K、含量3%)Fig.1 Curve of MSD-t of plasticizer DOS(298 K，content 3%)

2 增塑剂扩散系数的加速老化实验测定

2.1 实验内容

固体火箭研制定型，需要评定其贮存寿命，在没有充分的现场试验结果时，通过加速贮存试验预测固体火箭使用寿命是一个重要的途径。本文采用对粘接试件进行加速老化试验方法，测定衬层中增塑剂的含量。按照QJ 2038.1A—2004矩形件制作方法制作丁羟推进剂/衬层/绝热层/壳体模拟粘接试件，增塑剂含量为3%、衬层配方与模拟部分中丁羟粘合系统一致(HTPB-IPDI)，衬层加厚到10 mm，丁羟推进剂加厚到200 mm。将粘接试件分别在25、50、70℃温度下贮存50 d后，将衬层取下由内向外切成厚度约1 mm左右的切片，记录衬层切片到界面的距离，用内标法测定切片中增塑剂的浓度。采用同样方法测定推进剂中增塑剂含量为2%、3%、4%的粘接试件在70℃温度下贮存50 d后衬层切片中增塑剂的浓度。

2.2 计算扩散系数

根据费克第二定律计算增塑剂扩散系数。在丁羟推进剂/衬层界面附近区域，假设增塑剂扩散符合一维情况，扩散系数D在整个扩散过程中不变，且聚合物体积不发生变化;由于推进剂相对衬层而言质量大得多，推进剂作为增塑剂扩散迁移的扩散源，可认为推进剂中增塑剂浓度始终不变;初始时衬层中增塑剂含量为零。可采用以下边界条件来求解费克第二定律的表达式［14］。

边界条件:

费克第二定律的表达式:

式中t为扩散时间，s;x为扩散距离，m;C为扩散剂浓度，mol/m3;Q为推进剂中增塑剂的表面浓度，mol/m2。

求解得

又因为D、Q均为常数，所以当扩散时间t一定时，式(3)可简化为

作lnC-x2直线，由直线斜率可换算得到增塑剂表观扩散系数D［15-16］。

2.3 扩散系数的实验结果

根据式(4)对实验数据进行处理，得到lnC-x2直线，由斜率可计算出增塑剂含量为3%，25、50、70℃温度下增塑剂扩散系数分别为 9.08×10-12、9.68×10-12、10.18 ×10-12m2/s。随温度升高，扩散系数依次升高。计算70℃温度下增塑剂含量为2%、3%、4%时，扩散系数分别为 9.68 ×10-12、10.18 ×10-12、9.37×10-12m2/s，扩散系数的变化趋势不明显，增塑剂含量为3%时，扩散系数相比而言最大。表1是增塑剂含量为3%、70℃，老化50 d后，衬层不同位置测得的增塑剂浓度和扩散系数。略去其他试验数据。

3 结果与讨论

表2是通过分子动力学模拟及加速老化实验方法得到不同温度下增塑剂DOS的扩散系数，DOS扩散系数模拟值和实验值数量级为10-12m2/s。随温度升高，DOS扩散系数均逐渐增大。增塑剂与推进剂粘接体系相容性好［17］，由亲合性拉力引起扩散运动，增塑剂不参与固化系统的反应，在高分子链段运动产生的间隙中(即自由体积孔穴之间)以“跳跃”方式发生空间位置迁移［18］。这在用分子动力学方法模拟DOS分子在HTPB-IPDI体系中的微观运动过程中得到直观演示。随温度升高，DOS扩散系数均逐渐增大的原因，经分析主要有2点:一是对于粘接体系而言，温度升高，运动单元热运动能量随之升高，聚合物链段间的活动性增加，分子间距离增加，聚合物中的自由体积膨胀、溶胀度升高，使得增塑剂DOS进行扩散的有效空间增大;二是温度升高，增塑剂分子的活动性增加，加上粘接体系的粘度下降，使得DOS的扩散行为更活跃。

表3是通过分子动力学模拟及加速老化实验方法得到不同增塑剂含量下DOS的扩散系数，增塑剂含量为2%、3%、4%时，实验得到的扩散系数变化不明显，含量为3%时的扩散系数相比而言最大;当推进剂中增塑剂的含量为3%、6%、9%、15%(已超过推进剂中的真实含量2% ～4%)时，模拟得到的扩散系数依次略有下降。

表1 衬层不同位置的增塑剂浓度和扩散系数Table 1 Plasticizer concentration in the different position of liner and diffusion coefficient

表2 不同温度下的增塑剂扩散系数(含量3%)Table 2 Plasticizer diffusion coefficient on the different temperatures(content 3%)

表3 不同增塑剂含量下的扩散系数Table 3 Plasticizer diffusion coefficient at different contents

共容增塑限量主要取决于粘合剂相的网络结构，在固体火箭发动机推进剂中通常可加入2%～4%的增塑剂，推进剂中粘合剂相的含量一般不超过10%，那么增塑剂约占其粘接体系的20% ～35%(质量比)，这大概就是粘接体系的共容增塑限量［1］。而增塑剂在衬层中其含量极低，甚至为零，界面处浓度差引起的一种拉曳力会使增塑剂从推进剂迁移到衬层中。衬层中粘合体系约占40%，根据这个共容增塑限量值，可知衬层最多可吸收约8% ～14%的增塑剂，这与试验现象(增塑剂迁移的相对含量反而减小，详见文献［19］)吻合。

在交联聚合物的溶胀过程中，存在着2种相反趋势的平衡过程:增塑剂力图渗入高聚物内部，使高聚物体积膨胀，引起分子网交联点之间分子链的伸展，降低高聚物的构象熵值;构象熵值降低又必然引起分子网的弹性收缩力，力图使分子网收缩。当粘接体系中DOS含量超过共容增塑限量后，DOS的增加对溶胀度的影响很小，几乎保持不变;过量的DOS占据大量空间，不仅束缚分子网的伸展，而且促使分子网收缩，使其自由体积减少，降低DOS的扩散活动。所以，扩散系数有所下降。

4 结论

(1)增塑剂扩散系数试验值和模拟值的数量级基本一致，都为10-12m2/s。

(2)随温度升高，DOS扩散系数均逐渐增大。

(3)增塑剂含量为2%、3%、4%时，试验得到的扩散系数变化不明显;当推进剂中增塑剂的含量为3%、6%、9%、15%，模拟得到的扩散系数依次略有下降。

(4)虽然试验过程中人为因素干扰较大、模拟计算的体系较小，且对实物作了必要简化，但从试验结果和模拟结果的对比分析可知，分子模拟可得到与实验相一致的结果，增塑剂扩散系数的试验结果可验证其分子模拟结果，分子动力学方法计算增塑剂扩散系数更具优势。

［1］徐思羽，唐大森.二茂铁衍生物在固体推进剂应用中的迁移问题［J］.推进技术，1983，28(3):40-45.

［2］尹华丽，李东峰，王玉，等.组分迁移对NEPE推进剂界面粘接性能的影响［J］.固体火箭技术，2005，28(2):126-129.

［3］尹华丽，李东峰，张纲要，等.丁羟推进剂粘接体系中的组分迁移［J］.固体火箭技术，2007，30(5):416-419.

［4］姚维尚，李倩，谭惠民.NEPE推进剂粘合剂性能的分子模拟研究［J］.含能材料，2007，15(6):650-655.

［5］杨月诚，焦东明，强洪夫，等.HTPB推进剂组分溶度参数的分子模拟研究［J］.含能材料，2008，16(2):191-195.

［6］李红霞，强洪夫，等.丁羟推进剂粘接体系中增塑剂迁移的分子模拟［J］.火炸药学报，2008，31(5):74-78.

［7］李红霞，强洪夫，王广，武文明.基于MD方法的增塑剂扩散行为的模拟研究［J］.含能材料，2009，17(1):36-41.

［8］Kotelyanskii M J，Wagner N J，Paulaitis M E.Molecular dynamics simulation study of the mechanisms of water diffusion in dydrated，amorphous polyamide［J］.Computational and Theoretical Polymer Science，1999，9:301-306.

［9］Hofmann D，Fritz L，Ulbrich J，et al.Molecular simulation of small molecule diffusion and solution in dense amorphous polysiloxanes and polyimides［J］.Computational and Theoretical Polymer Science，2000，10:419-436.

［10］Guo Xiang-dan，Teng Jia-wei，Xie Zai-ku.Investigation of water adsorbed on HZSM-5:experiment and molecular simulation［J］.Computer and Applied Chemistry，2006，6:486-490.

［11］Liu Qing-lin，Huang Yu.Transport behavior of oxygen and nitrogen through organasilicon-containing polystyrenes by molecular simulation［J］.J.Phys.Chem.B，2006，100:17375-17382.

［12］于坤千，李泽生，李志儒，等.寡聚物在高分子母体中的扩散——分子动力学模拟研究［J］.高等学校化学学报，2002，23(7):1327-1330.

［13］Allen M P，Tildesley D J.Computer simulation of liquids［M］.Oxford:Clarendon Press，1987.

［14］王国梅，万发荣.材料物理［M］.武汉:武汉理工大学出版社，2004:31.

［15］Gonzalez S，Gil M A，Hernandez J O，et al.Resistance to corrosion of galvanized steel covered with an epoxy-polyamide primer coating［J］.Progress in Organic Coatings，2001，41:3167-3170.

［16］Glavchev I，Petrova K，Devedjiev I.Investigation of solvent interaction in epoxy coatings［J］.Polymer Testing，2000，19:111-114.

［17］张景春.固体推进剂化学及工艺学［M］.长沙:国防科技大学出版社，1987:36.

［18］潘福生.高分子-无机杂化膜扩散特性的分子动力学模拟［D］.天津:天津大学，2006.

［19］李红霞.丁羟推进剂/衬层粘接界面失效的分子模拟与实验研究［D］.西安:第二炮兵工程学院，2009:92-93.

Measurement of diffusion coefficient of plasticizer in HTPB propellant

LI Hong-xia1，QIANG Hong-fu1，LI Xin-qi1，WANG Hong-wei2
(1.The Second Artillery Engineering University，Xi＇an 710025，China;2.The 46th Institute of the Sixth Academy of CASIC，Huhhot 010076，China)

Aiming at the plasticizer migration in the interface between HTPB propellant and liner，the diffusibility of plasticizer was studied.The diffusion coefficients were obtained by two methods.Firstly，plasticizer diffusivity in the bond system of HTPB propellant was simulated by the molecular dynamics method and the diffusion coefficient was obtained by Einstein equation.Secondly，taking the bonding interface of propellant and thickening liner as the object，the accelerated migration test of plasticizer was performed and the diffusion coefficient was obtained by the second law(Fick).The diffusion coefficients of dioctyl sebacate were discussed under different temperatures and at different contents of plasticizer.It was analyzed that plasticizer migrated by jump on the space situation which were produced when macromolecule chain moved in the HTPB bond system.The results indicated that the magnitude order of DOS diffusion coefficient was 10-12m2/s.The diffusion coefficients gradually increased with increasing of environmental temperature because of more intense activity and larger active space of plasticizer.Yet the diffusion coefficients decreased a little with the increasing of DOS content when content was more than 3%.It is more advantageous to measure diffusion coefficient of plasticizer by molecule dynamics.

HTPB propellant;plasticizer;diffusion coefficient;experiment;simulation

V512

A

1006-2793(2012)03-0387-04

2011-01-23;

2012-03-06。

国家863项目(2010AAJ119);国家973项目(973-61338)。

李红霞(1979—)，女，博士，研究方向为纳米材料应用与含能材料模拟。E-mail:lhxlihongxia@163.com

(编辑:刘红利)