下载文档原格式
交联度对酸酐固化环氧树脂热机械性能影响的分子动力学模拟律方成; 付可欣; 张磊; 赵云晓; 刘博闻; 谢庆【期刊名称】《《华北电力大学学报(自然科学版)》》【年(卷),期】2019(046)006【总页数】7页(P1-7)【关键词】酸酐固化环氧树脂体系; 分子动力学模拟; 交联度; 热机械参数【作者】律方成; 付可欣; 张磊; 赵云晓; 刘博闻; 谢庆【作者单位】华北电力大学新能源电力系统国家重点实验室北京102206; 华北电力大学河北省输变电设备安全防御重点实验室河北保定071003【正文语种】中文【中图分类】TM216.30 引言酸酐类固化剂固化环氧树脂所得到的环氧固化物具有机械性能好、易加工成型和成本低廉的优点，尤其是绝缘性能优异。
采用甲基四氢苯酐(MTHPA)固化双酚A型环氧树脂(DGEBA)所得到的环氧树脂固化物作为绝缘材料广泛应用于高压绝缘系统[1]。
随着超/特高压直流输电系统的快速发展，对环氧树脂基材的综合性能提出了更高的要求[2]。
通常采用填充改性制备复合材料来提升环氧树脂的热机械性能，然而大量的填充难以同时获得优异的热机械性能，甚至会带来绝缘特性的下降[3]。
因此，如果能制备出热机械学性能优异的环氧树脂基材，则可以大大减少填料的填充量，为获得综合性能优异的环氧树脂复合材料奠定基础。
环氧树脂基材的性能主要取决于交联网络分子结构[4-6]，交联网络分子结构可以通过改变交联度来实现[7]。
传统的实验方法可以通过实验测量环氧树脂基材的宏观性能参数，如玻璃化转化温度(Tg)、导热系数、弹性模量、拉伸强度等[8, 9]，却难以表征环氧树脂交联网络的分子结构、分子构象变化，难以定量分析交联度变化对环氧树脂基材性能的影响，难以研究环氧树脂微观结构与宏观性能间的关联关系。
近年来，在材料研究领域，大量研究表明分子模拟技术在研究材料组成、微观结构与性能关系上具有明显优势。
分子模拟方法主要包括：适用于简单分子或电子数量较少体系的量子力学方法和以经典力学为主的非量子力学计算方法。
分子动力学在沥青中的应用摘要：对于沥青再生技术而言，再生中新旧沥青的混溶过程是一个关键问题，其直接影响着合适的新沥青胶结料等级及其用量的确定，并对保障沥青再生效果及沥青路面的路用性能具有重要意义。
为此本文采用分子动力学方法，从分子尺度构建新旧沥青混溶模型，在微观模拟新旧沥青混溶过程的同时，计算分析再生沥青的宏观物理性能。
关键词：再生沥青；分子动力学；新旧沥青混溶；再生剂1.绪论1.1.研究背景和意义公路作为一项为公众提供出行便利性的基础设施，一直是国家建设的重点工程，是衡量一个国家经济水平和发展状况的标志之一。
目前我国公路的主要路面类型为沥青混凝土路面，约占我国公路路面总面积的80%，而随着我国交通运输业的不断发展及公路的不断使用，已经有很大一部分的沥青路面面临着养护维修乃至重建，甚至有些路面在短期内就已经受到了定程度的损坏，需要进行养护以延长使用寿命。
这就导致了每年将产生大量的废旧沥青混合料。
在沥青混合料再生技术中，一个关键问题是再生过程中新旧沥青混合问题，其混合程度直接影响再生混合料的性能指标，对整个沥青混合料的性能指标有直接的影响，旧沥青是否完全参与再生，对于新沥青的等级、掺量、混合料设计等一系列问题起基础性的作用。
研究新旧沥青混合问题对沥青混合料再生技术有着重大的意义，是再生理论的基础性内容之一，对再生体系的完善起促进作用。
目前关于新旧沥青混溶状态的研究大多集中在宏观实验研究，但由于沥青的分子组成十分繁复，新旧沥青混溶的过程也极为复杂，宏观实验难以较好地表征新旧沥青的混溶状态，一些微观试验手段例如原子力显微镜、扫描电镜等，虽然可以观察到微观始末的状态，但其内部结构的变化过程却难以观测。
因此，分子动力学是研究新旧沥青混溶状态一个直接有效的研究手段。
而分子动力学通过定义粒子支架相互作用，通过对体系内分子运动进行计算，来得到体系的热力学性质。
利用分子动力学研究新旧沥青混溶状态，在对分子微观行为进行研究的同时，也可以计算物质的宏观物理性质。
溶液中柔性树枝状高分子的分子动力学模拟张国梁;卢宇源;宋建晖;石彤非【摘要】采用分子动力学模拟方法研究了柔性树枝状高分子在无热溶剂中的静态和动态行为. 模拟结果表明: 在分子尺寸和回转半径R_g满足标度律R_g～N~(1/5)(G+1)~(2/5)P~(2/5)(其中N为树枝状分子的聚合度, G为代数, P为链节长度, g为子代代数)时;随着代数的增加, 树枝状分子和硬球的静态结构因子相似, 表明其内部结构发生了由"类星形"向"近球形"转化. 随着树枝状分子代数和链节长度的增加, 出现了"单元"(Monomer)密度几乎不变的区域. 树枝状分子的回折能力随着链节长度的增加而增强, 随着代数的增加而减弱. 树枝状分子各子代的运动能力不同, 与内层子代相比, 外层子代在短时间内扩散较慢, 但其松弛较快. 相对于Rouse指数, 树枝状分子"单元"运动的标度更接近Zimm指数.【期刊名称】《高等学校化学学报》【年(卷),期】2010(031)001【总页数】6页(P199-204)【关键词】分子动力学模拟;树枝状高分子;Zimm指数【作者】张国梁;卢宇源;宋建晖;石彤非【作者单位】中国科学院长春应用化学研究所,高分子物理与化学国家重点实验室,长春,130022;中国科学院研究生院,北京,100049;中国科学院长春应用化学研究所,高分子物理与化学国家重点实验室,长春,130022;中国科学院长春应用化学研究所,高分子物理与化学国家重点实验室,长春,130022;中国科学院研究生院,北京,100049;中国科学院长春应用化学研究所,高分子物理与化学国家重点实验室,长春,130022【正文语种】中文【中图分类】O631;O641树枝状高分子是通过多官能团单体逐步聚合得到的具有高度支化结构的大分子[1,2].由于树枝状分子在医药载体、基因工程、自组装超分子体系、纳米材料等领域具有广泛的应用前景[3],近年来引起人们的广泛关注[4～7].自 de Gennes等[8]的开辟性工作以来,人们通过实验[9～13]、理论和模拟[14～23]对树枝状分子进行深入研究.Bauer等[9]用小角 X射线散射(SAXS)发现:当代数变大时,树枝状分子聚酰胺2胺(PAMAM)的内部结构发生了类星形到球形的转变.Chai等[11]和Moreno等[12]用核磁共振(NMR)的方法研究了树枝状分子的动态行为,发现各子代有不同的移动能力,外层子代松弛的更快. Murat等[13]用随机力来替代真实溶剂分子的作用,发现外层子代部分向后折叠.Sheng等[16]用蒙特卡罗(MC)模拟研究了溶剂性质对树枝状分子尺寸的影响,证实了回转半径 Rg满足 Flory标度关系.结合高斯自洽场(GSCF)和蒙特卡罗模拟,Timoshenko等[17]研究了溶液中的树枝状分子,当代数变大时,在距中心不远处的“单元”密度存在极小值.Karatasos等[22]将溶剂分子作为真实的粒子,用分子动力学方法研究了树枝状分子的静态和动力学性质.目前,系统阐述代数和链节长度对静态和动力学性质影响的研究较少.由于大量的树枝状高分子为柔性链节,链节向后折叠呈现浓核结构,而树枝状分子为内部稀疏的浓壳结构,因此大量的研究工作集中于实现浓核2浓壳的结构转变.树枝状分子的向后折叠能力是浓核2浓壳结构转变的关键因素,研究链节长度和代数对回折能力的影响,对认识结构转变意义重大.本文研究了链节长度和代数对树枝状分子静态和动力学性质的影响,获得回折能力与链长和代数之间的关系,并揭示“单元”运动行为和各子代的扩散和松弛行为的不同成因,为树枝状高分子在相关领域的应用提供前瞻性探索.分子动力学是以牛顿第二定律为基础的微观模拟方法.所有粒子(溶质与溶剂)间的相互作用用截断2平移(Lennard2Jones)势表示:在模拟中采用约化单位,ε,σ,m(溶剂粒子的质量)均设置为 1,溶质粒子的质量设置为 2.其中ε决定了相互作用的强度,σ是长度的单位,时间的单位为τLJ=σ(m/ε)1/2.截断半径采用rc=21/6σ[24,25],树枝状高分子链上相邻粒子的连接作用为非线性弹性势 (FENE)[14]:式中,弹簧的弹性系数 k为 30,链节上相邻粒子的最大距离 R0为 115,链节上相邻粒子的距离为 r.粒子的运动采用蛙跳 (Leap2frog)[26]算法,步长为 01006.采用很高的粒子数密度 018638[24],通过Berendsen热浴[27],体系温度控制在KBT=112.溶液体系通过2×106步的松弛达到平衡,平衡后的112×107步中每200步抽样对体系的物理量进行计算.模拟体系的格子尺寸对树枝状分子的质心扩散系数 Dcm影响很大[23,24].根据式 (3),通过控制回转半径与格子尺寸的比值(Rg/L接近 0125)来减小格子尺寸的影响.图 1是链节长度 P=2,代数 G=1的树枝状高分子的拓扑结构示意图.其中圆内部分为 G=0子代,外层圆环内部分为 G=1子代.P3G6简化表示 P=3,G=6的树枝状高分子,其它分子也采用类似的简化方式.聚合度N与代数 G和链节长度 P的关系如下:2.1 树枝状高分子的大小和内部结构为了考察 G和 P对树枝状高分子大小和结构的影响,计算了回转半径 Rg和静态结构因子S(q).图 2(A)是代数 G不变,仅改变 P时 Rg与N的双对数曲线,其中 Rg的标度指数为3/5,与稀溶液中线性高分子的标度律相同[28].图 2(B)是 Rg/[P (G+1)]2/5与 N的双对数曲线,其中Rg/[P(G+1)]2/5的标度指数为 1/5,故得到柔性树枝状高分子在无热溶剂中的标度律:又因为N∝P[式(4)],故在 G不变时的标度律为 3/5.为了深入理解树枝状高分子的内部结构,计算了静态结构因子 S(q):图 3(A)和 (B)分别为 P=1和 G=5时树枝状分子 Kratky形式的静态结构因子.图3(A)中自上而下5条曲线的代数G依次为2,3,4,5,6.可见当G较小时,其静态结构因子与星形分子相似.随着G的增加,S(q)在横坐标为 510附近出现了明显的二级峰(如图中箭头所示),与硬球的静态结构因子[9]相似,表明树枝状分子的结构向近球形转化,该结果与Bauer等[9]用小角X光散射(SAXS)和 Scherren2 berg等[26]用小角中子散射 (SANS)得到的实验结果一致.图 3(B)中自下而上 4条曲线的 P依次为1, 2,3,6.可见随着 P的增加,静态结构因子值变大,二级峰消失,这说明树枝状分子需要更大的代数来实现其结构向近球形转化.2.2 树枝状高分子的“单元”数密度和回折能力由于树枝状高分子的“单元”数密度和回折能力决定了其大小和内部结构,定量计算“单元”数密度和回折能力对深入理解树枝状高分子的静态性质十分重要.为了定量分析代数G和链节长度 P对树枝状高分子空间分布和回折能力的影响,计算了树枝状高分子“单元”数密度和空间扩张量 S[14].图 4(A)为 P=3的树枝状分子“单元”数密度与距质心距离的曲线图.当 G=2时,数密度随 r增加而单调递减.而当 G=6时,数密度在 r=210处出现极小值,在 4<r<6范围内数密度值变化很小(密度平台区),随着 r的增加数密度逐渐减小到零.数密度平台区的出现是树枝状分子回折的结果,其范围随着G增加而变大,而平台区的密度值几乎不变.图4(B)是G=5时树枝状分子的“单元”数密度分布曲线.质心附近密度太大导致远离质心处各树枝状分子密度曲线区分不清,故图 4中忽略了质心附近密度.结果表明:当 2<r<5时,数密度随着 P的增加而减小.数密度的平台区域随着 P 增加而变大,而平台区的数密度值随着 P的增加而减小.图 5(A)中的插图为链节长度 P=4的树形分子,有箭头的粗线为子代键矢量,细线为树形分子链节.其中枝化点是树形分子分叉点,g和 g+1为相邻枝化点的代数,代数相邻的支化点连接即为子代键矢量.空间扩张量 S定义为相邻的子代键矢量 Ig+1投影到 Ig.由式(7)可知,S值越小,树枝状分子的回折能力越强.图 5为空间扩张量 S与枝化点代数 g的关系图.在图 5(A)中,代数相同枝化点的S值随链长 P增加而减小,说明各子代的回折能力随着 P的增加而增强,即树枝状分子的回折能力随着 P的增加而增强.图5(B)中树枝状分子的最外层、次外层和倒数第三外层子代的 S值随着 G的增加而单调递增,说明树枝状分子的回折能力随着 G的增加而减弱.当代数g≥1时,S 值随着 g增加而减小,表明子代回折能力随着代数的增加而增强.2.3 树枝状高分子的扩散和其它的动态行为采用 Karatasos[23]的方法,计算了各子代的均方位移.在短时间内,随着子代代数 g 的增加,子代的质心扩散明显变慢.这可以用 Stoke2Einstein[28,29]关系进行解释: 式(8)表明:粒子的尺寸越大,其扩散越慢.各子代的大小用 Rg估算,将各子代当作高分子,计算,其中下标 m,n只取该子代链节的下标值).以 P=2,G=6为例,当子代代数g从 0变化到 6时,各子代的 Rg值依次为1131,2176,3192,4187,5162,6121,6162.由于 Rg值随 g单调增加,故子代的扩散随着 g的增加而变慢.利用Dünweg等[24]的方法,计算了树枝状高分子的质心扩散系数Dcm.图 6为Dcm与聚合度N的双对数曲线图.其中Dcm近似满足 Zimm标度关系[28],标度指数为 -0157.P和 G对整体扩散的影响很小,如图 6中的 P=3,G=6和 P=6,G=5,虽然 G和 P不同,但聚合度N很接近,其Dcm几乎相同.考虑到树枝状高分子的特殊支化结构,我们通过自关联函数来深入探讨其松弛行为.其中键矢量是相邻子代枝化点间连线的单位矢量.采用 Karatasos[23]的方法,分析键矢量的自关联函数 C(t)[26]:式中,v(0)和 v(t)是 0时刻和 t时刻矢量 ^v的单位矢量.局域小尺度的动力学用键矢量自关联函数的松弛进行分析.图 7为 P=3,G=6子代键矢量的自关联函数曲线.可见子代键矢量的 C(t)值随着代数 g增加而减小,说明子代松弛随着 g的增加而变快.由于内层子代分布于“单元”密度很大的树枝状高分子内部,而外层子代大部分位于分子表层,表层密度较小,空间位阻效应较小,故外层子代松弛较快.Chai等[11]和Moreno等[12]用核磁共振 (NMR)研究了树枝状分子的松弛行为,发现外层子代松弛的更快.模拟结果不仅佐证了实验结果,而且给出了合理解释.通过动态结构因子[24]S(k,t)研究树枝状高分子“单元”的运动,S(k,t)形式如下:式中,|ri(t)-rj(0)|是 t时刻离子 i与 0时刻离子 j的距离,当 Rg-1ν kν a-1时,S(k,t)满足如下关系:Zimm和 Rouse模型中 z值分别为 310和 2+1/v.图 8(A)和(B)分别给出了 Zimm和 Rouse模型下 P=3,G=6的动态结构因子.相对于 Rouse模型,Zimm模型能更好地描述其在稀溶液中的动力学行为,其它树枝状分子的 S(k,t)也得到了类似结果.由于高分子“单元”在溶剂中运动时,拖曳了周围溶剂一起运动.Rouse模型忽略了“单元”与溶剂的黏性阻力,而 Z imm模型则考虑了其与溶剂的流体力学效应,故能更真实地描述树枝状分子“单元”在稀溶液中的动力学行为[28].综上所述,采用分子动力学模拟研究了无热溶剂中的柔性树枝状高分子,系统考察了代数 G和链节长度 P对其静态和动态性质的影响.结果表明,对于分子尺寸,当改变P时,回转半径 Rg满足 Rg～N3/5的标度律,而当 P和 G都改变时,Rg满足标度律Rg～N1/5(G+1)2/5P2/5.随着代数 G的增加,静态结构因子出现了明显的二级峰,与硬球的静态结构因子接近,证实了树枝状分子的结构向“近球形”转化.当代数 G 很小时,树枝状分子“单元”密度从质心到外层递减.但是随着 G和 P的增加,出现了“单元”密度几乎不变的平台区域.平台区域的密度值在 G增加时不变,在 P增加时变小.树枝状高分子的回折能力随着 P的增加而增强,随着 G的增加而减弱,其子代回折能力随着子代代数 g增加而增强.树枝状分子的质心扩散与“单元”的运动满足 Zimm标度关系.与内层子代相比,外层子代在较小时间尺度内其扩散变慢,但其松弛速度较快.[1] Tomalia D.A..Prog.Polym.Sci.[J],2005,30:294—324[2] BallauffM.,Likos C.N..Angew.Chem.Int.Ed.[J],2004,43:2998—3020[3] HelmsB.,Meijer E.J..Science[J],2006,313:929—930[4] ZHANGQi2Zhen(张其震),SUN Ji2Run(孙继润),WANGDa2Qing(王大庆),et al..Chem.J.Chinese Universities(高等学校化学学报)[J],1998,19(7):1175—1177 [5] LONG Fei(龙飞),FAN Yu(范瑜),D ING Hui2Jun(丁慧君),etal..Chem.J.Chinese Universities(高等学校化学学报)[J], 1999,20(3):495—497[6] Lee H.,Baker J.R.,Jr.,Larson R.G..J.Phys.Chem.B[J],2006,110(9):4014—4019[7] Brocorens P.,Lazzaroni R.,Br?dasJ.L..J.Phys.Chem.B[J],2007,111(31):9218—9227[8] de Gennes P.G.,Hervet H..J.Phys.Lett.[J],1983,44(9):L351—L360[9] Prosa T.J.,BauerB.J.,Amis E.J..Macromolecules[J],2001,34(14):4897—4906[10] Topp A.,BauerB.J.,Tomalia D.A.,.etal..Macromolecules[J],1999,32(21):7232—7237[11] ChaiM.,Niu Y.,et al..J.Am.Chem.Soc.[J],2001,123(20):4670—4678[12] Moreno K.X.,Simanek E.E..Macromolecules[J],2008,41(12):4108—4114[13] Scherrenberg R.,CoussensB.,et al..Macromolecules[J],1998,31(2):456—461[14] MuratM.,Grest G.S..Macromolecules[J],1996,29(4):1278—1285[15] ZacharopoulosN.,Economou I.G..Macromolecules[J],2002,35(5):1814—1821[16] Sheng Y.J.,Jiang S.Y.,Tsao H.K..Macromolecules[J],2002,35(21):7865—7868[17] MansfieldM.L.,Klushin L.I..Macromolecules[J],1993,26(16):4262—4268[18] Timoshenko E.G.,Klushin Y.A.,ConnollyR..J.Chem.Phys.[J],2002,117(19):9050—9062[19] MansfieldM.L..Macromolecules[J],2000,33(21):8043—8049[20] Ferla R.L..J.Chem.Phys.[J],1997,106(2):688—700[21] Ganazzoli F.,Ferla R.L.,TerragniG..Macromolecules[J],2000,33(17):6611—6620[22] Ganazzoli F.,Ferla R.L.,RaffainiG..Macromolecules[J],2001,34(12):4222—4228[23] Karatasos K.,AdolfD.B.,Davies G.R..J.Chem.Phys.[J],2001,115(11):5310—5318[24] DünwegB.,Kremer K..J.Chem.Phys.[J],1993,99(9):6983—6997[25] Fu C.L.,OuyangW.Z.,Sun Z.Y.,et al..J.Chem.Phys.[J],2007,127,044903[26] AllenM.P.,puter SimulationofLiquids[M],Oxford:Clarendon Press,1987[27] Berendsen H.J.C.,Postma J.P.M.,et al..J.Chem.Phys.[J],1984,81(8):3684—3690[28] RubinsteinM.,Colby R.H..Polymer Physics[M],New York:OxfordUniversity,2004[29] DoiM.,Edwards S.F..The Theory ofPolymerDynamics[M],Oxford:Clarendon Press,1996【相关文献】[1] Tomalia D.A..Prog.Polym.Sci.[J],2005,30:294—324[2] BallauffM.,Likos C.N..Angew.Chem.Int.Ed.[J],2004,43:2998—3020[3] HelmsB.,Meijer E.J..Science[J],2006,313:929—930[4] ZHANGQi2Zhen(张其震),SUN Ji2Run(孙继润),WANGDa2Qing(王大庆),etal..Chem.J.Chinese Universities(高等学校化学学报)[J],1998,19(7):1175—1177[5] LONG Fei(龙飞),FAN Yu(范瑜),D ING Hui2Jun(丁慧君),et al..Chem.J.Chinese Universities(高等学校化学学报)[J], 1999,20(3):495—497[6] Lee H.,Baker J.R.,Jr.,Larson R.G..J.Phys.Chem.B[J],2006,110(9):4014—4019[7] Brocorens P.,Lazzaroni R.,Br?das J.L..J.Phys.Chem.B[J],2007,111(31):9218—9227[8] de Gennes P.G.,Hervet H..J.Phys.Lett.[J],1983,44(9):L351—L360[9] Prosa T.J.,BauerB.J.,Amis E.J..Macromolecules[J],2001,34(14):4897—4906[10] Topp A.,BauerB.J.,Tomalia D.A.,.et al..Macromolecules[J],1999,32(21):7232—7237[11] ChaiM.,Niu Y.,et al..J.Am.Chem.Soc.[J],2001,123(20):4670—4678[12] Moreno K.X.,Simanek E.E..Macromolecules[J],2008,41(12):4108—4114[13] Scherrenberg R.,CoussensB.,et al..Macromolecules[J],1998,31(2):456—461[14] MuratM.,Grest G.S..Macromolecules[J],1996,29(4):1278—1285[15] ZacharopoulosN.,Economou I.G..Macromolecules[J],2002,35(5):1814—1821[16] Sheng Y.J.,Jiang S.Y.,Tsao H.K..Macromolecules[J],2002,35(21):7865—7868[17] MansfieldM.L.,Klushin L.I..Macromolecules[J],1993,26(16):4262—4268[18] Timoshenko E.G.,Klushin Y.A.,Connolly R..J.Chem.Phys.[J],2002,117(19):9050—9062[19] MansfieldM.L..Macromolecules[J],2000,33(21):8043—8049[20] Ferla R.L..J.Chem.Phys.[J],1997,106(2):688—700[21] Ganazzoli F.,Ferla R.L.,Terragni G..Macromolecules[J],2000,33(17):6611—6620[22] Ganazzoli F.,Ferla R.L.,Raffaini G..Macromolecules[J],2001,34(12):4222—4228[23] Karatasos K.,AdolfD.B.,Davies G.R..J.Chem.Phys.[J],2001,115(11):5310—5318[24] DünwegB.,Kremer K..J.Chem.Phys.[J],1993,99(9):6983—6997[25] Fu C.L.,OuyangW.Z.,Sun Z.Y.,et al..J.Chem.Phys.[J],2007,127,044903[26] AllenM.P.,puter Simulation ofLiquids[M],Oxford:Clarendon Press,1987[27] Berendsen H.J.C.,Postma J.P.M.,et al..J.Chem.Phys.[J],1984,81(8):3684—3690[28] RubinsteinM.,Colby R.H..Polymer Physics[M],New York:Oxford University,2004[29] DoiM.,Edwards S.F..The Theory of PolymerDynamics[M],Oxford:Clarendon Press,1996 (Ed.:D,Z)。
第３９卷第３期２０２０年３月硅㊀酸㊀盐㊀通㊀报ＢＵＬＬＥＴＩＮＯＦＴＨＥＣＨＩＮＥＳＥＣＥＲＡＭＩＣＳＯＣＩＥＴＹＶｏｌ.３９㊀Ｎｏ.３Ｍａｒｃｈꎬ２０２０ＳｉＣ模具高温模压石英玻璃物相接触角的分子动力学模拟吴㊀悠１ꎬ２ꎬ３ꎬ邹㊀斌１ꎬ２ꎬ３ꎬ王俊成１ꎬ２ꎬ３ꎬ黄传真１ꎬ２ꎬ３ꎬ朱洪涛１ꎬ２ꎬ３ꎬ姚㊀鹏１ꎬ２ꎬ３(１.山东大学机械工程学院先进射流工程技术研究中心ꎬ济南㊀２５００６１ꎻ２.山东大学高效洁净机械制造教育部重点实验室ꎬ济南㊀２５００６１ꎻ３.山东大学机械工程国家级实验教学示范中心ꎬ济南㊀２５００６１)摘要:针对高温下石英玻璃纳米液滴在ＳｉＣ模具表面接触角难以测量的问题ꎬ采用分子动力学方法ꎬ模拟研究了不同温度和粗糙表面面向模压的ＳｉＯ２/ＳｉＣ高温接触角以及ＳｉＯ２熔体的界面结构ꎮ应用压力张量法发现了ＭＳ￣Ｑ势函数模拟的ＳｉＯ２熔体表面张力较接近实际值ꎬ即ＳｉＯ２高温表面特性模拟可优先采用ＭＳ￣Ｑ势函数ꎮ针对ＳｉＣ模具纳米级表面的粗糙度ꎬ发现当粗糙度因子ｒ>１.５时润湿模式由Ｗｅｎｚｅｌ变为Ｃａｓｓｉｅ￣Ｂａｘｔｅｒꎬ此时Ｒａ的变化对接触角值无明显影响ꎬＲｍｒ值减小使得接触面积分数ｆ减小ꎬ接触角值随之增大ꎮ因此ꎬ保持ｒ大于１.５的同时适当减小Ｒｍｒ值有利于减小固液摩擦ꎬ降低石英玻璃工件和ＳｉＣ模具界面上的脱模力ꎮ随着温度升高ＳｉＯ２表面结构变得松散ꎬ导致其在ＳｉＣ表面接触角减小ꎮ在超过２３００Ｋ时接触角值的变化率增大ꎬ为减小工件￣模具界面的粘附ꎬ模压温度应选择２３００Ｋ以下ꎮ关键词:石英玻璃ꎻ碳化硅ꎻ接触角ꎻ分子动力学ꎻ模压中图分类号:Ｏ６４７.１ꎻＴＳ９４３.６５㊀㊀文献标识码:Ａ㊀㊀文章编号:１００１￣１６２５(２０２０)０３￣０９２３￣０９ＭｏｌｅｃｕｌａｒＤｙｎａｍｉｃｓＳｉｍｕｌａｔｉｏｎｏｆＨｉｇｈＴｅｍｐｅｒａｔｕｒｅＣｏｎｔａｃｔＡｎｇｌｅｆｏｒＭｏｌｄｅｄＳｉｌｉｃａＧｌａｓｓｉｎＳｉＣＤｉｅＷＵＹｏｕ１ꎬ２ꎬ３ꎬＺＯＵＢｉｎ１ꎬ２ꎬ３ꎬＷＡＮＧＪｕｎｃｈｅｎｇ１ꎬ２ꎬ３ꎬＨＵＡＮＧＣｈｕａｎｚｈｅｎ１ꎬ２ꎬ３ꎬＺＨＵＨｏｎｇｔａｏ１ꎬ２ꎬ３ꎬＹＡＯＰｅｎｇ１ꎬ２ꎬ３(１.ＣｅｎｔｅｒｆｏｒＡｄｖａｎｃｅｄＪｅｔＥｎｇｉｎｅｅｒｉｎｇＴｅｃｈｎｏｌｏｇｉｅｓꎬＳｃｈｏｏｌｏｆＭｅｃｈａｎｉｃａｌＥｎｇｉｎｅｅｒｉｎｇꎬＳｈａｎｄｏｎｇＵｎｉｖｅｒｓｉｔｙꎬＪｉｎａｎ２５００６１ꎬＣｈｉｎａꎻ２.ＫｅｙＬａｂｏｒａｔｏｒｙｏｆＨｉｇｈＥｆｆｉｃｉｅｎｃｙａｎｄＣｌｅａｎＭｅｃｈａｎｉｃａｌＭａｎｕｆａｃｔｕｒｅꎬＳｈａｎｄｏｎｇＵｎｉｖｅｒｓｉｔｙꎬＭｉｎｉｓｔｒｙｏｆＥｄｕｃａｔｉｏｎꎬＪｉｎａｎ２５００６１ꎬＣｈｉｎａꎻ３.ＮａｔｉｏｎａｌＤｅｍｏｎｓｔｒａｔｉｏｎＣｅｎｔｅｒｆｏｒＥｘｐｅｒｉｍｅｎｔａｌＭｅｃｈａｎｉｃａｌＥｎｇｉｎｅｅｒｉｎｇＥｄｕｃａｔｉｏｎꎬＳｈａｎｄｏｎｇＵｎｉｖｅｒｓｉｔｙꎬＪｉｎａｎ２５００６１ꎬＣｈｉｎａ)Ａｂｓｔｒａｃｔ:Ｄｕｅｔｏｔｈｅｄｉｆｆｉｃｕｌｔｉｅｓｏｆｍｅａｓｕｒｅｍｅｎｔｏｆｃｏｎｔａｃｔａｎｇｌｅｗｈｅｎｐｒｅｓｓｉｎｇｓｉｌｉｃａｇｌａｓｓｎａｎｏ￣ｄｒｏｐｌｅｔｉｎＳｉＣｄｉｅｓｕｒｆａｃｅａｔｈｉｇｈｔｅｍｐｅｒａｔｕｒｅꎬｔｈｅＳｉＯ２/ＳｉＣｈｉｇｈｔｅｍｐｅｒａｔｕｒｅｃｏｎｔａｃｔａｎｇｌｅａｔｄｉｆｆｅｒｅｎｔｔｅｍｐｅｒａｔｕｒｅａｎｄｒｏｕｇｈｓｕｒｆａｃｅａｓｗｅｌｌａｓｔｈｅｓｕｒｆａｃｅｓｔｒｕｃｔｕｒｅｏｆＳｉＯ２ｍｅｌｔｗｅｒｅｓｉｍｕｌａｔｅｄｂｙｍｏｌｅｃｕｌａｒｄｙｎａｍｉｃｓ.ＵｓｉｎｇｔｈｅｐｒｅｓｓｕｒｅｔｅｎｓｏｒｍｅｔｈｏｄꎬｉｔｉｓｆｏｕｎｄｔｈａｔｔｈｅｓｕｒｆａｃｅｔｅｎｓｉｏｎｏｆＳｉＯ２ｍｅｌｔｓｉｍｕｌａｔｅｄｂｙＭＳ￣Ｑｐｏｔｅｎｔｉａｌｆｕｎｃｔｉｏｎｉｓｃｌｏｓｅｔｏｔｈｅｅｘｐｅｒｉｍｅｎｔａｌｖａｌｕｅ.ＴｈｕｓｔｈｅＭＳ￣ＱｐｏｔｅｎｔｉａｌｃａｎｂｅｐｒｅｆｅｒｒｅｄｉｎｓｉｍｕｌａｔｉｏｎｏｆＳｉＯ２ｈｉｇｈ￣ｔｅｍｐｅｒａｔｕｒｅｓｕｒｆａｃｅｃｈａｒａｃｔｅｒ.ＦｏｒｔｈｅｒｏｕｇｈｎｅｓｓｏｆｎａｎｏｓｃａｌｅｓｕｒｆａｃｅｏｆＳｉＣｄｉｅ.Ｗｈｅｎｔｈｅｒｏｕｇｈｎｅｓｓｆａｃｔｏｒｒｅｘｃｅｅｄｓ１.５ꎬｔｈｅｗｅｔｔｉｎｇｍｏｄｅｃｈａｎｇｅｓｆｒｏｍＷｅｎｚｅｌｔｏＣａｓｓｉｅ￣Ｂａｘｔｅｒ.ＡｔｔｈｉｓｔｉｍｅꎬｔｈｅｃｈａｎｇｅｏｆＲａｈａｓｎｏｓｉｇｎｉｆｉｃａｎｔｉｍｐａｃｔｏｎｔｈｅｃｏｎｔａｃｔａｎｇｌｅꎬａｎｄｔｈｅｒｅｄｕｃｅｄＲｍｒｖａｌｕｅｒｅｓｕｌｔｓｉｎｔｈｅｄｅｃｒｅａｓｅｏｆｃｏｎｔａｃｔａｒｅａｆｒａｃｔｉｏｎｆａｎｄｔｈｅｉｎｃｒｅａｓｅｏｆｃｏｎｔａｃｔａｎｇｌｅ.Ｔｈｅｒｅｆｏｒｅꎬｋｅｅｐｉｎｇｒｇｒｅａｔｅｒｔｈａｎ１.５ｗｈｉｌｅａｐｐｒｏｐｒｉａｔｅｌｙｒｅｄｕｃｉｎｇＲｍｒｖａｌｕｅｉｓｃｏｎｄｕｃｉｖｅｔｏｒｅｄｕｃｉｎｇｓｏｌｉｄ￣ｌｉｑｕｉｄｆｒｉｃｔｉｏｎａｎｄｒｅｄｕｃｉｎｇｔｈｅｓｔｒｉｐｐｉｎｇｆｏｒｃｅｏｎｔｈｅｉｎｔｅｒｆａｃｅｂｅｔｗｅｅｎｆｕｓｅｄｓｉｌｉｃａｗｏｒｋｐｉｅｃｅａｎｄＳｉＣｄｉｅ.ＡｓｔｈｅｔｅｍｐｅｒａｔｕｒｅｉｎｃｒｅａｓｅｓꎬｔｈｅｓｕｒｆａｃｅｓｔｒｕｃｔｕｒｅｏｆＳｉＯ２ｂｅｃｏｍｅｌｏｏｓｅꎬｒｅｓｕｌｔｉｎｇｉｎｉｔｓｃｏｎｔａｃｔａｎｇｌｅｄｅｃｌｉｎｉｎｇｏｎｔｈｅＳｉＣｓｕｒｆａｃｅ.Ｔｈｅｃｈａｎｇｅｒａｔｅｏｆｔｈｅｃｏｎｔａｃｔａｎｇｌｅｉｎｃｒｅａｓｅｓｗｈｅｎｔｅｍｐｅｒａｔｕｒｅｅｘｃｅｅｄｓ２３００Ｋ.Ｉｎｏｒｄｅｒｔｏｒｅｄｕｃｅｔｈｅａｄｈｅｓｉｏｎｏｆｔｈｅｗｏｒｋｐｉｅｃｅｗｉｔｈｔｈｅｄｉｅꎬｔｈｅｍｏｌｄｉｎｇｔｅｍｐｅｒａｔｕｒｅｓｈｏｕｌｄｂｅｂｅｌｏｗ２３００Ｋ.Ｋｅｙｗｏｒｄｓ:ｓｉｌｉｃａｇｌａｓｓꎻｓｉｌｉｃｏｎｃａｒｂｉｄｅꎻｃｏｎｔａｃｔａｎｇｌｅꎻｍｏｌｅｃｕｌａｒｄｙｎａｍｉｃｓꎻｃｏｍｐｒｅｓｓｉｏｎｍｏｌｄｉｎｇ基金项目:山东省自然科学基金重大基础研究项目(ＺＲ２０１８ＺＢ０５２１)作者简介:吴㊀悠(１９９４￣)ꎬ男ꎬ硕士研究生ꎮ主要从事模具摩擦磨损方面的研究ꎮＥ￣ｍａｉｌ:ｊｏｌｍｕｇｉ＠１６３.ｃｏｍ通讯作者:邹㊀斌ꎬ教授ꎮＥ￣ｍａｉｌ:ｚｂ７８＠ｓｄｕ.ｅｄｕ.ｃｎ９２４㊀陶㊀瓷硅酸盐通报㊀㊀㊀㊀㊀㊀第３９卷０㊀引㊀言近年来ꎬ微结构光学元件因其在信息通信领域的广泛应用受到越来越多的关注ꎮ模压是制造光学玻璃元件常用的方法之一ꎮ模压过程中ꎬ一方面模具上精准的几何结构复印到了玻璃元件上ꎻ另一方面玻璃元件与模压模具之间的摩擦与粘附造成了模具的磨损ꎬ从而导致模具几何结构精度的丧失ꎮ这些均为模具与工件的固液界面相互接触㊁相对运动作用的结果ꎮ因此ꎬ有必要研究它们的接触特性ꎬ从而对模具结构和模压工艺的设计提供一定的指导ꎮ润湿是能够实现模压的先决条件ꎬ而接触角是表征润湿性的主要参数ꎮ对于理想光滑表面的接触角模型ꎬＹｏｕｎｇ[１]从固㊁液㊁气界面张力平衡的角度建立了经典的杨氏方程ꎮ对于粗糙表面和非均质表面ꎬＷｅｎｚｅｌ[２]和Ｃａｓｓｉｅ[３]等学者在杨氏方程的基础上分别建立了Ｗｅｎｚｅｌ和Ｃａｓｓｉｅ￣Ｂａｘｔｅｒ模型ꎮＷｅｎｚｅｌ润湿模式中ꎬ液体始终完全浸渍粗糙结构波谷中ꎬ即出现所谓 钉扎 现象ꎮ而在模压工艺中ꎬ这一现象增大了工件液相在模具表面流动的阻力ꎬ从而增大了固液摩擦ꎬ使得模具更易磨损ꎮ而Ｃａｓｓｉｅ￣Ｂａｘｔｅｒ润湿模式认为由于粗糙结构波谷中存在气相ꎬ液体无法浸渍其中ꎬ因此是模压比较期望的润湿模式ꎮ在光学玻璃材料中ꎬ石英玻璃的光学性能有其独特之处ꎮ它既可以透过远紫外光谱ꎬ且透射率在同类材料中最优ꎬ又可透过可见光和近红外光谱ꎮ同时ꎬ其机械性能也高于普通玻璃[４]ꎮ目前ꎬ石英玻璃光学元件的制造多采用超精密磨削或是激光刻蚀的方法ꎬ生产效率低ꎬ生产成本高ꎮ而模压方法可以提高加工效率ꎬ实现大规模生产ꎮ然而ꎬ由于石英玻璃软化点温度高达１５００~１６００ħꎬ常用的模具材料如ＷＣ等在此温度下无法正常工作ꎬ而ＳｉＣ热稳定性好ꎬ熔点达２８３０ħ[５]ꎬ故较为适合作为模压石英玻璃的模具材料ꎮ由于模压温度高ꎬ超出了目前接触角测量仪器的工作范围ꎬ并且微结构模具表面粗糙度达到纳米级ꎬ要研究粗糙度对润湿性的影响需要测量纳米尺度微液滴的接触角ꎬ因此试验测定比较困难ꎮ分子动力学方法使微纳液滴高温接触角的模拟及研究其成因机理成为可能ꎮ徐威等[６]通过改变ＬＪ势函数中的作用参数模拟了纳米水滴在不同能量表面上的铺展过程和润湿形态ꎬ模拟结果与经典润湿理论计算得到的结果呈现相似变化趋势ꎮ王龙[７]研究了铜和金液滴在石墨烯和碳纳米管等不同结构基底上的润湿和融合过程ꎬ发现金属液滴的融合受液滴形状和基底结构影响ꎮ由于势函数的限制ꎬ目前对于界面润湿的分子动力学模拟多集中于ＬＪ流体ꎬ如水和液氩等ꎬ而多元素流体较少见于报道ꎮ势函数按多体作用的复杂程度可分为对势和多体势ꎮ对于石英玻璃的模拟ꎬ多体势不能较好地计算Ｓｉ￣Ｏ键的断裂ꎬ因此不利于高温动力学性能的研究[８]ꎮＢｅｅｓｔ等[９]基于石英玻璃经典的ＢＭＨ势提出了ＢＫＳ势ꎬ并运用第一性原理方法结合实验数据拟合了参数ꎮＳｕｎｄａｒａｒａｍａｎ等[１０]使用ＢＫＳ势预测了石英玻璃的力学性能ꎬ发现当短程和长程截止半径分别为５.５Å和１０Å时模拟的石英玻璃结构更符合实际ꎮＤｅｍｉｒａｌｐ等[１１]首先将Ｍｏｒｓｅ势与电荷平衡法(ＱＥｑ)相结合ꎬ建立了ＭＳ￣Ｑ力场ꎬ并研究了石英玻璃在压力变化过程中的相变ꎮ丁元法等[１２]比较了ＢＫＳ与ＭＳ￣Ｑ模型下石英玻璃的高温扩散特性ꎬ认为计算高温下石英玻璃的扩散传输性能可优先选择ＭＳ￣Ｑ力场ꎬ但并未比较其表面特性ꎮ从以上文献可以看出ꎬＢＫＳ和ＭＳ￣Ｑ是石英玻璃分子动力学模拟中常用的对势ꎮ其中ꎬＢＫＳ势从低温到高温的较大温度范围均具有较好性能ꎬ而ＭＳ￣Ｑ势在高温下的传输性能要优于ＢＫＳ势ꎮ目前对于石英玻璃高温表面特性的模拟的报道较为少见ꎬ因此仍需要对两种势函数性能的优劣进行比较ꎮＨｏｓｓｅｉｎｉ等[１３]研究了不同形貌的疏水表面上的水滴行为ꎬ发现表面形貌㊁柱高度㊁空隙率和沉积角是影响表面疏水性的主要参数ꎮ因此ꎬ本模拟在ＳｉＣ表面构建了纳米方柱阵列ꎬ并分别使用粗糙度评定参数Ｒａ和Ｒｍｒ表示柱高和空隙率ꎬ研究了纳米级表面粗糙结构对面向高温模压的ＳｉＯ２/ＳｉＣ接触角的影响ꎮ本文使用分子动力学方法ꎬ对比了分别采用ＢＫＳ和ＭＳ￣Ｑ势函数计算的ＳｉＯ２熔体高温表面张力ꎮ模拟了不同模压温度和ＳｉＣ模具表面粗糙结构ＳｉＯ２的高温接触角ꎬ研究了微尺度下的高温界面润湿特征和ＳｉＯ２熔体表面层结构特点ꎬ本研究的结构框图如图１所示ꎮ１㊀模拟方法和模型１.１㊀界面张力的计算为了保证成形精度并减小工件与模具的粘附ꎬ光学玻璃模压温度一般选择在转变点温度到软化点温度第３期吴㊀悠:ＳｉＣ模具高温模压石英玻璃物相接触角的分子动力学模拟９２５㊀之间[１４]ꎬ故选择石英玻璃实际软化点附近温度１９００Ｋ作为模拟温度ꎮＳｉＯ２表面张力模拟系统如图２所示ꎮＳｉＯ２熔体表面模型的建立方法是在模拟盒子中随机加入８００个Ｓｉ原子和１６００个Ｏ原子ꎬ使其密度约为２.２ｇ/ｃｍ３ꎮ随后在５０００Ｋ下驰豫２００ｐｓ以消除初始构型的影响ꎬ并以１００Ｋ/２０ｐｓ的速度降温至１９００Ｋꎮ最后在ｘ轴方向模型两边各添加一厚度为２０Å的真空层ꎬ以避免周期性边界的影响ꎮ面向高温模压的ＳｉＯ２/ＳｉＣ界面张力模拟系统如图３(ｃ)所示ꎬ其中ＳｉＣ是由金刚石结构的β￣ＳｉＣ原胞排列而成ꎮ分别将ＳｉＯ２和ＳｉＣ单独在１９００Ｋ下驰豫２００ｐｓꎬ然后将它们添加进同一模拟盒子中ꎬ使ＳｉＯ２与ＳｉＣ{１００}表面接触并对整个体系进行驰豫ꎮ图１㊀基于分子动力学的ＳｉＣ模具高温模压石英玻璃的物相接触角模拟研究关系框图Ｆｉｇ.１㊀ＳｔｕｄｙｒｅｌａｔｉｏｎｓｈｉｐｏｆｈｉｇｈｔｅｍｐｅｒａｔｕｒｅｃｏｎｔａｃｔａｎｇｌｅｏｆｍｏｌｄｅｄｏｐｔｉｃａｌｇｌａｓｓｉｎＳｉＣｄｉｅｂａｓｅｄｏｎｍｏｌｅｃｕｌａｒｄｙｎａｍｉｃｓ图２㊀ＳｉＯ２熔体表面张力模拟计算体系Ｆｉｇ.２㊀ＳｉｍｕｌａｔｉｏｎｓｙｓｔｅｍｏｆＳｉＯ２ｍｅｌｔｓｕｒｆａｃｅｔｅｎｓｉｏｎ图３㊀ＳｉＯ２/ＳｉＣ高温界面张力模拟计算体系Ｆｉｇ.３㊀ＳｉｍｕｌａｔｉｏｎｓｙｓｔｅｍｏｆＳｉＯ２/ＳｉＣｈｉｇｈｔｅｍｐｅｒａｔｕｒｅｉｎｔｅｒｆａｃｅｔｅｎｓｉｏｎ㊀㊀表面与界面张力的计算均采用压力张量法[１５]ꎬγ＝１２ʏＬｘ２－Ｌｘ２[ＰＮ(ｘ)－ＰＴ(ｘ)]ｄｘ(１)式中ꎬγ表示表面张力ꎻＬｘ为模拟体系在ｘ方向的长度ꎻ因子１/２是由于模拟系统有两个界面ꎻＰＮ与ＰＴ分别为系统界面法向和切向压力张量分量ꎮＰＮ(ｘ)＝ρｋＢＴ－１Ｖðｉ<ｊｘ２ｉｊｒｉｊｄＵ(ｒｉｊ)ｄｒｉｊ[](２)９２６㊀陶㊀瓷硅酸盐通报㊀㊀㊀㊀㊀㊀第３９卷ＰＴ(ｘ)＝ρｋＢＴ－１Ｖðｉ<ｊｙ２ｉｊ－ｚ２ｉｊ２ｒｉｊｄＵ(ｒｉｊ)ｄｒｉｊ[](３)式中ꎬｋＢ为玻尔兹曼常数ꎻＶ为模拟盒子体积ꎻＵ表示势函数ꎻρ为液体密度ꎻＴ为温度ꎻｒ为原子间距ꎬｘ㊁ｙ㊁ｚ为其在三个坐标轴方向的分量ꎮ为了比较不同势函数模拟界面张力的可靠性ꎬ在模拟中ꎬＳｉＯ２原子间的相互作用先后用ＢＫＳ[１０]和ＭＳ￣Ｑ[１１]势函数进行表征如下ꎬ势函数参数见表１ꎮＵＢｉｊ＝ｑｉｑｊｅ２ｒｉｊ＋Ａｉｊｅ－ｒｉｊρ－Ｃｉｊｒ６ｉｊ(４)ＵＭｉｊ＝ｑｉｑｊｅ２ｒｉｊ＋Ｄ０[ｅ－２α(ｒｉｊ－ｒ０)－２ｅ－α(ｒｉｊ－ｒ０)](５)式中ꎬｑ表示电荷量ꎻＡ㊁Ｃ㊁Ｄ均为与相互作用强度有关的参数ꎻα为与平衡距离有关的参数ꎮ表１㊀ＳｉＯ２势函数参数Ｔａｂｌｅ１㊀ＰａｒａｍｅｔｅｒｓｏｆｐｏｔｅｎｔｉａｌｓｆｏｒＳｉＯ２ＰａｉｒＢＫＳＡｉｊ/(ｋｃａｌ ｍｏｌ－１)ρ/ÅＣｉｊ/(ｋｃａｌ ｍｏｌ－１ Å６)ＭＳ￣ＱＤ０/(ｋｃａｌ ｍｏｌ－１)α/Å－１ｒ０/ÅＳｉ￣Ｓｉ￣￣￣０.１７７３２.０４４３.７６０Ｏ￣Ｏ３２０２５.７９９８０.３６２３４０３５.５８７５０.５３６３１.３７３３.７９１Ｓｉ￣Ｏ４１５１７５.６４２９０.２０５２３０７９.４５５４４５.９９７０２.６５２１.６２８㊀㊀ＳｉＣ原子间的相互作用使用ｔｅｒｓｏｆｆ势函数表征ꎮ由于ＳｉＯ２/ＳｉＣ界面间为范德华作用ꎬ因此使用ＬＪ势函数表征ꎬ其参数来自于ＵＦＦ力场[１６]ꎬ如表２所示ꎮ表２㊀ＳｉＯ２/ＳｉＣ势函数参数Ｔａｂｌｅ２㊀ＰａｒａｍｅｔｅｒｓｏｆｐｏｔｅｎｔｉａｌｆｏｒＳｉＯ２/ＳｉＣｉｎｔｅｒａｃｔｉｏｎＰａｉｒε/(ｋｃａｌ ｍｏｌ－１)σ/ÅＳｉ￣Ｓｉ０.４０２３.８２６Ｏ￣Ｓｉ０.１５５３３.４７２３Ｓｉ￣Ｃ０.２０５３.６２８Ｏ￣Ｃ０.０７９３３.２７４５㊀㊀模拟过程中ꎬ模拟盒子三个方向均为周期性边界ꎮＢＫＳ势函数Ｓｉ和Ｏ原子所带电荷分别为＋２.４ｅ㊁－１.２ｅꎬ截断半径为５.５ÅꎻＭＳ￣Ｑ势函数Ｓｉ和Ｏ原子所带电荷分别为＋１.３１８ｅ㊁－０.６５９ｅꎬ截断半径为９Åꎮ静电力的计算采用ＰＰＰＭ(Ｐａｒｔｉｃｌｅ￣ＰａｒｔｉｃｌｅＰａｒｔｉｃｌｅ￣Ｍｅｓｈ)方法且精度为１０－５ꎬ截断半径为１０Åꎮ时间步长为１ｆｓꎮ模拟均在正则系综(ＮＶＴ)下进行ꎬ调温使用Ｎｏｓé￣Ｈｏｏｖｅｒ方法ꎮ１.２㊀接触角的模拟图４㊀面向高温模压的ＳｉＯ２/ＳｉＣ接触角模拟系统Ｆｉｇ.４㊀ＳｉｍｕｌａｔｉｏｎｓｙｓｔｅｍｏｆＳｉＯ２/ＳｉＣｃｏｎｔａｃｔａｎｇｌｅｆｏｒｈｉｇｈｔｅｍｐｅｒａｔｕｒｅｍｏｌｄｉｎｇ在模拟系统中ＳｉＯ２与ＳｉＣ分别为液相和固相ꎬ其中ꎬＳｉＣ基底为ＳｉＣ单胞排列而成的超晶胞ꎬ对其作出三点假设:(１)假设ＳｉＣ表面为不存在缺陷的理想表面ꎻ(２)由于ＳｉＣ不同晶面表面能相差较小ꎬ假设其{１００}晶面为与ＳｉＯ２相接触的表面ꎻ(３)假设ＳｉＯ２熔体液滴为理想的球形ꎮ接触角的模拟系统如图４所示ꎬ首先在半径２.３４８Å的球体区域随机添加１２００个Ｓｉ原子和２４００个Ｏ原子ꎬ随后在５０００Ｋ下驰豫２００ｐｓ并以１００Ｋ/２０ｐｓ的速度降温至１９００Ｋ得到模拟所用的ＳｉＯ２液滴模型ꎮ最后将ＳｉＯ２和ＳｉＣ添加进同一模拟盒子ꎬｘ㊁ｙ方向使用周期性边界ꎬｚ方向使用固定和镜像边界ꎬ其他第３期吴㊀悠:ＳｉＣ模具高温模压石英玻璃物相接触角的分子动力学模拟９２７㊀设置与上一节相同ꎮ图５㊀微结构阵列模具表面多尺度形貌示意图Ｆｉｇ.５㊀Ｍｕｌｔｉｓｃａｌｅｓｕｒｆａｃｅｔｏｐｏｇｒａｐｈｙｏｆｍｉｃｒｏ￣ｓｔｒｕｃｔｕｒｅａｒｒａｙｄｉｅ为了研究纳米级表面粗糙结构对接触角的影响ꎬ采用文献[１３]的方法构建了方柱形阵列ＳｉＣ壁面ꎮ其结构如图５所示ꎮＷｅｎｚｅｌ[２]和Ｃａｓｓｉｅ￣Ｂａｘｔｅｒ[３]模型可用式(６)㊁(７)表示:ｃｏｓθＷ＝ｒｃｏｓθＹ(６)ｃｏｓθＣ＝ｆｃｏｓθＹ＋ｆ－１(７)式中ꎬθＹ为本征接触角ꎻθＷ与θＣ为Ｗｅｎｚｅｌ和Ｃａｓｓｉｅ￣Ｂａｘｔｅｒ模型接触角ꎻｒ为粗糙度因子ꎬ表示表面的实际面积与表观面积之比ꎻｆ为接触面积分数ꎬ表示液滴与基底实际接触的面积和基底表观面积之比ꎮｒ与ｆ可用粗糙度评定参数表示为ｒ＝１＋８ＲｍｒＲａＳｍ(８)ｆ＝Ｒｍｒ２(９)式中ꎬＲａ为轮廓算数平均偏差ꎬＲｍｒ为轮廓支承长度率ꎬＳｍ为轮廓微观不平度的平均间距ꎮ模拟中若Ｒｍｒ㊁Ｓｍ取恒值ꎬ则Ｒａ取０.２５㊁０.５㊁１㊁１.５倍晶格参数ꎬ分别对应ｒ值为１.５㊁２.０㊁３.０㊁４.０ꎻ若Ｒａ取恒值ꎬ改变柱间距为２㊁１.５㊁１㊁０.５倍晶格参数使Ｒｍｒ取０.３３㊁０.４０㊁０.５０㊁０.６７ꎬ分别对应ｆ值为０.１１㊁０.１６㊁０.２５㊁０.４４ꎮ２㊀模拟结果与讨论２.１㊀界面张力图６为１９００Ｋ下ＢＫＳ与ＭＳ￣Ｑ模型中ＳｉＯ２熔体表面张力随时间的演化ꎮ由于使用压力张量法时系统需要较长时间稳定ꎬ模拟连续进行了２０ｎｓꎮ为了降低模拟初期结构不稳定的影响ꎬ在计算累计平均值时采用最后１０ｎｓ数据ꎮＢＫＳ与ＭＳ￣Ｑ势函数的计算结果分别为０.４１０Ｎ/ｍ和０.３９０Ｎ/ｍꎮ由文献[１７]可知ꎬＳｉＯ２在１９００Ｋ下的实际表面张力应为０.３０２Ｎ/ｍ左右ꎮ考虑到模拟本身的准确性以及实际实验过程中ꎬ氧分子等充当表面活性剂降低了测定的表面张力值ꎬ因此模拟值高出实验值０.０７~０.１１Ｎ/ｍ左右是可接受的[１８]ꎮ可以看出ＭＳ￣Ｑ势函数的模拟结果比ＢＫＳ势函数更为接近实验值ꎮ图６㊀ＢＫＳ势函数与ＭＳ￣Ｑ势函数下ＳｉＯ２熔体表面张力演化Ｆｉｇ.６㊀ＥｖｏｌｕｔｉｏｎｏｆＳｉＯ２ｍｅｌｔｓｕｒｆａｃｅｔｅｎｓｉｏｎｕｓｉｎｇＢＫＳｐｏｔｅｎｔｉａｌａｎｄＭＳ￣Ｑｐｏｔｅｎｔｉａｌ图７为两种不同势函数模拟的面向模压的ＳｉＯ２/ＳｉＣ高温界面张力演化图ꎮ界面张力系统达到稳定时间相对较短ꎬ因此模拟时间设为１５ｎｓꎬ取最后４.５ｎｓ时间段的数据计算平均值ꎮ从图中可见ꎬＢＫＳ和ＭＳ￣Ｑ势函数的模拟结果分别为０.８４６Ｎ/ｍ和０.６８２Ｎ/ｍꎮ同时ꎬ计算了两种表面模型的一维密度分布和径向分布函数(ＲａｄｉａｌＤｉｓｔｒｉｂｕｔｉｏｎＦｕｎｃｔｉｏｎꎬＲＤＦ)ꎮ图８(ａ)为模型中心至模拟盒子边缘的密度分布ꎬＢＫＳ与ＭＳ￣Ｑ模型的表面均存在类似汽液共存界面的低密度相ꎬ其厚度分别为４Å与６ÅꎮＢＫＳ模型内部的密度约为９２８㊀陶㊀瓷硅酸盐通报㊀㊀㊀㊀㊀㊀第３９卷２.６９ｇ/ｃｍ３ꎬ大于ＭＳ￣Ｑ的２.２５ｇ/ｃｍ３ꎮ这说明ＢＫＳ表面模型的密度大于ＭＳ￣Ｑ模型ꎬ因此原子间距要小于ＭＳ￣Ｑ模型ꎬ具有更强的原子间作用力ꎬ从而具有更大的表面张力ꎮ图８(ｂ)中ＲＤＦ的计算结果也支持了这一解释ꎮＲＤＦ图中前三个峰代表Ｏ￣Ｏ㊁Ｓｉ￣Ｏ㊁Ｓｉ￣Ｓｉ原子对的第一近邻ꎬ峰值点的横坐标即原子对的键长ꎬ其数值如图８(ｂ)所示ꎬＭＳ￣Ｑ模型中数量较多且强度较大的Ｓｉ￣Ｏ㊁Ｏ￣Ｏ键长均大于ＢＫＳ模型ꎬ从而导致其表面张力小于ＢＫＳ模型ꎮ综上所述ꎬ高温模压时ＳｉＯ２表面性质的模拟可优先选择ＭＳ￣Ｑ势函数ꎮ图７㊀面向高温模压的ＳｉＯ２/ＳｉＣ界面张力演化过程Ｆｉｇ.７㊀ＥｖｏｌｕｔｉｏｎｏｆＳｉＯ２/ＳｉＣｉｎｔｅｒｆａｃｅｔｅｎｓｉｏｎｆｏｒｈｉｇｈｔｅｍｐｅｒａｔｕｒｅｍｏｌｄｉｎｇ图８㊀使用ＢＫＳ与ＭＳ￣Ｑ势函数的ＳｉＯ２表面结构对比Ｆｉｇ.８㊀ＣｏｍｐａｒｉｓｏｎｏｆＳｉＯ２ｓｕｒｆａｃｅｓｔｒｕｃｔｕｒｅｕｓｉｎｇＢＫＳａｎｄＭＳ￣Ｑｐｏｔｅｎｔｉａｌ２.２㊀表面粗糙结构对接触角的影响液滴在理想光滑表面上的接触角称为本征接触角ꎮ由于微纳尺度下液滴表面存在一定的密度和压力涨落ꎬ为了获得具有统计学意义的接触角值ꎬ在处理数据时采用等密度拟合曲线法[１９]ꎮ最后得到面向高温模压的ＳｉＯ２/ＳｉＣ的本征接触角为１１９.２５ʎꎮ根据杨氏方程[１]ꎬ可以利用上节模拟得到的界面张力对本征接触角进行检验ꎮｃｏｓθ＝γｓｖ－γｓｌγｌｖ(１０)式中ꎬγｓｖ表示固￣气界面能ꎻγｓｌ与γｌｖ分别表示固￣液界面能和气￣液界面能ꎬ它们在数值上与固￣液界面张力和液体表面张力相等ꎮ由上节可知ꎬγｓｌ和γｌｖ的值分别为０.６８２Ｎ/ｍ和０.３９０Ｎ/ｍꎮ由于本模拟中采用的ＳｉＯ２/ＳｉＣ相互作用参数来自ＵＦＦ力场ꎬ因此选择文献[２０]中使用ＵＦＦ力场算得的ＳｉＣ{１００}晶面表面能数值０.５０２Ｎ/ｍ代入式(１０)ꎮ最终得到的本征接触角理论值为１１７.４９ʎꎬ与模拟数值差别不大ꎮ图９是ＳｉＯ２熔体在具有不同粗糙度因子的ＳｉＣ模具表面上接触角的模拟ꎬ从图中可见ꎬ接触角总体值在１３５ʎ左右波动ꎮ接触角变化不大则表明固液相互作用能较为稳定ꎬ但图９(ｂ)显示ꎬ在ｒ＝１.５时固液相互作用能较大ꎮ由图１０可以看出ꎬ当粗糙度因子ｒ＝１.５时ꎬＳｉＯ２熔体浸润了沟槽ꎬ出现 钉扎 现象ꎬ此时润湿第３期吴㊀悠:ＳｉＣ模具高温模压石英玻璃物相接触角的分子动力学模拟９２９㊀接近于Ｗｅｎｚｅｌ状态ꎮ 钉扎 现象增强了固液界面的摩擦ꎬ液滴的铺展需要克服更大的能垒ꎬ因此虽然固液相互作用能较大但并未使液滴接触角减小ꎻｒ>１.５时ꎬＳｉＯ２均处于Ｃａｓｓｉｅ￣Ｂａｘｔｅｒ润湿状态ꎬ因此接触角不再受粗糙度因子变化的影响ꎮ从Ｃａｓｓｉｅ￣Ｂａｘｔｅｒ状态到Ｗｅｎｚｅｌ状态的过渡称为润湿转变ꎮ润湿转变可以通过施加压力实现ꎬ但临界转变压力会随着纳米柱的高度减小而减小[２１]ꎻ当其减小到一定程度ꎬ就可能仅仅依靠系统压力实现润湿转变ꎮ本研究中ꎬ当ｒ>１.５时系统压力无法使液滴保持Ｗｅｎｚｅｌ状态ꎬ因此转变为Ｃａｓｓｉｅ￣Ｂａｘｔｅｒ状态ꎻ而此状态下ꎬ界面的固液相互作用比Ｗｅｎｚｅｌ状态更低ꎬ液滴在固体表面的扩散系数变小[２２]ꎬ说明此时ＳｉＯ２熔体在剪切作用下更易在ＳｉＣ模具表面滑移ꎬ从而导致粘性摩擦作用较小ꎮ图９㊀ＳｉＣ模具表面粗糙度对ＳｉＯ２熔体液滴接触角和两者相互作用能的影响Ｆｉｇ.９㊀ＥｆｆｅｃｔｏｆＳｉＣｄｉｅｒｏｕｇｈｎｅｓｓｏｎｃｏｎｔａｃｔａｎｇｌｅｏｆＳｉＯ２ｍｅｌｔｄｒｏｐｌｅｔａｎｄｉｎｔｅｒａｃｔｉｏｎｅｎｅｒｇｙ图１０㊀ＳｉＣ模具表面粗糙度对ＳｉＯ２熔体液滴接触状态的影响Ｆｉｇ.１０㊀ＥｆｆｅｃｔｏｆＳｉＣｄｉｅｒｏｕｇｈｎｅｓｓｏｎｃｏｎｔａｃｔｓｔａｔｅｏｆＳｉＯ２ｍｅｌｔｄｒｏｐｌｅｔ图１１㊀接触面积分数对面向高温模压的ＳｉＯ２/ＳｉＣ接触状态的影响Ｆｉｇ.１１㊀ＥｆｆｅｃｔｏｆｃｏｎｔａｃｔａｒｅａｆｒａｃｔｉｏｎｏｎＳｉＯ２/ＳｉＣｃｏｎｔａｃｔｓｔａｔｅｆｏｒｈｉｇｈｔｅｍｐｅｒａｔｕｒｅｍｏｌｄｉｎｇ图１２㊀接触面积分数对面向高温模压的ＳｉＯ２/ＳｉＣ接触角和相互作用能的影响Ｆｉｇ.１２㊀ＥｆｆｅｃｔｏｆｃｏｎｔａｃｔａｒｅａｆｒａｃｔｉｏｎｏｎＳｉＯ２/ＳｉＣｃｏｎｔａｃｔａｎｇｌｅａｎｄｉｎｔｅｒａｃｔｉｏｎｅｎｅｒｇｙｆｏｒｈｉｇｈｔｅｍｐｅｒａｔｕｒｅｍｏｌｄｉｎｇ９３０㊀陶㊀瓷硅酸盐通报㊀㊀㊀㊀㊀㊀第３９卷图１１为粗糙度因子ｒ>１.５时ＳｉＣ模具表面接触面积分数对ＳｉＯ２熔体液滴接触状态的影响ꎬ此时ＳｉＯ２始终呈现Ｃａｓｓｉｅ￣Ｂａｘｔｅｒ润湿状态ꎬ液滴随接触面积分数的增大逐渐在ＳｉＣ表面铺展ꎮ图１２(ａ)对比了模拟接触角和理论接触角ꎬ虽然两者存在一定误差ꎬ但其都随着接触面积分数的增大而减小ꎻ同时图１２(ｂ)显示接触角越小ꎬＳｉＯ２/ＳｉＣ高温模压界面具有越强的相互作用ꎮ表面粘着力和热应力是脱模力的两个组成部分[２３]ꎬ减小轮廓支承长度率即减小了ＳｉＯ２与ＳｉＣ模具界面的实际接触面积ꎬ从而能减小工件和模具之间的表面粘着力ꎮ而Ｃａｓｓｉｅ￣Ｂａｘｔｅｒ润湿模式无 钉扎 现象ꎬ减小了工件和模具的传热面积ꎬ在一定程度上减小了热应力ꎮ因此适当减小Ｒｍｒ值可以降低以及工件与模具之间的脱模力ꎬ从而减小模具磨损ꎬ提高寿命ꎮ同时也缩减了模具制造过程中抛光工序的工作量ꎮ２.３㊀温度对接触角的影响图１３㊀面向高温模压的ＳｉＯ２/ＳｉＣ接触角与温度的关系Ｆｉｇ.１３㊀ＲｅｌａｔｉｏｎｓｈｉｐｂｅｔｗｅｅｎＳｉＯ２/ＳｉＣｃｏｎｔａｃｔａｎｇｌｅａｎｄｔｅｍｐｅｒａｔｕｒｅｆｏｒｈｉｇｈｔｅｍｐｅｒａｔｕｒｅｍｏｌｄｉｎｇ图１３为接触面积分数为０.２５时温度对面向高温模压的ＳｉＯ２/ＳｉＣ接触角的影响ꎮ从图可见ꎬ接触角值随温度的增大而减小ꎬ说明温度升高ＳｉＯ２的表面张力减小ꎻ当温度高于２３００Ｋ时ꎬ接触角的变化较大ꎬ这是因为使用ＭＳ￣Ｑ势函数模型的ＳｉＯ２自扩散激活温度约为２３００Ｋ左右[１２]ꎮ自扩散激活温度可以视作石英玻璃的软化点温度ꎬ在此温度附近ꎬ石英玻璃的结构发生较大变化ꎬ松动和新生的化学键均大大增加[２４]ꎮ图１４(ａ)是不同温度下ＳｉＯ２熔体表面的结构特点ꎬ由图可见ꎬＳｉＯ２熔体表面存在一个等密度层ꎻ在等密度层外部ꎬ密度随温度升高而增大ꎻ在等密度层内部ꎬ密度随温度升高而减小ꎻ这意味着ＳｉＯ２熔体内部的原子在高温作用下逐渐向外扩散ꎬ这种密度梯度的减小使ＳｉＯ２表面结构更加松散ꎬ减小了其表面张力ꎮ图１４(ｂ)还表明ꎬ温度对ＳｉＯ２熔体表面原子化学键的键长并未产生较大影响ꎬ但随着温度升高ꎬ各化学键对应的峰值依次减小ꎬ表明表面层原子的无序程度增加ꎮ同时ꎬ对比图８(ｂ)可以发现ꎬ表面层中Ｓｉ￣Ｓｉ对应的峰值变得较为微弱甚至消失ꎬ这说明在表面层硅原子数量较少ꎬ而氧原子大量聚集ꎬＯ￣Ｏ键的强度远小于Ｓｉ￣Ｏ键强度ꎬ这可能是ＳｉＯ２表面张力随温度减小的另一个原因ꎮ图１４㊀不同温度下ＳｉＯ２熔体表面结构特点Ｆｉｇ.１４㊀ＳｕｒｆａｃｅｓｔｒｕｃｔｕｒｅｏｆＳｉＯ２ｍｅｌｔｉｎｄｉｆｆｅｒｅｎｔｔｅｍｐｅｒａｔｕｒｅ３㊀结㊀论采用分子动力学方法模拟了ＳｉＯ２熔体的界面结构ꎬ将ＳｉＣ模具纳米级表面理想化为纳米方柱阵列ꎬ研㊀第３期吴㊀悠:ＳｉＣ模具高温模压石英玻璃物相接触角的分子动力学模拟９３１究了粗糙度和温度对面向模压的ＳｉＯ２/ＳｉＣ高温接触角的影响ꎬ得到以下结论:(１)使用ＭＳ￣Ｑ势函数模拟的ＳｉＯ２熔体表面张力比ＢＫＳ势函数计算的结果与实验值更为接近ꎬ因此模拟ＳｉＯ２高温熔体的表面性质使用ＭＳ￣Ｑ势函数更为合理ꎮ(２)在１９００Ｋ的模压温度下ꎬ当粗糙度因子ｒ>１.５时ꎬＲａ的变化对接触角值无明显影响ꎮＲｍｒ值减小使得接触面积分数ｆ减小ꎬ接触角值随之增大ꎮ此时润湿模式从Ｗｅｎｚｅｌ转变为Ｃａｓｓｉｅ￣Ｂａｘｔｅｒꎬ减小了工件￣模具之间的摩擦ꎮ由于热应力和界面粘着力的减小ꎬ石英玻璃光学元件模压后的脱模力将会降低ꎮ同时ꎬ由于降低了Ｒｍｒ值的要求ꎬＳｉＣ模具加工过程中抛光工序的工作量也相应得到了减小ꎮ(３)面向高温模压的ＳｉＯ２/ＳｉＣ的接触角随模压温度升高而减小ꎮ当模压温度超过２３００Ｋ时ꎬ接触角变化率显著增大ꎮ因此模压温度选择在２３００Ｋ以下可以降低模压时模具因玻璃熔体的粘附造成的磨损ꎮ参考文献[１]㊀ＹｏｕｎｇＴＨ.Ａｎｅｓｓａｙｏｎｔｈｅｃｏｈｅｓｉｏｎｏｆｌｉｑｕｉｄｓ[Ｊ].Ｐｈｉｌ.Ｔｒａｎｓ.Ｒｏｙ.Ｓｏｃ.Ｌｏｎｄｏｎꎬ１８０５ꎬ９５:６５￣８７.[２]㊀ＷｅｎｚｅｌＲＮ.Ｒｅｓｉｓｔａｎｃｅｏｆｓｏｌｉｄｓｕｒｆａｃｅｓｔｏｗｅｔｔｉｎｇｂｙｗａｔｅｒ[Ｊ].Ｉｎｄ.Ｅｎｇ.Ｃｈｅｍ.ꎬ１９３６ꎬ２８:９８８￣９４.[３]㊀ＣａｓｓｉｅＡＢＤꎬＢａｘｔｅｒＳ.Ｗｅｔｔａｂｉｌｉｔｙｏｆｐｏｒｏｕｓｓｕｒｆａｃｅｓ[Ｊ].Ｔｒａｎｓ.Ｆａｒａｄａｙ.Ｓｏｃ.ꎬ１９４４ꎬ４０:５４６￣５１.[４]㊀陈㊀光.新材料概论[Ｍ].北京:科学出版社ꎬ２００３:４６￣４８.[５]㊀和丽芳.纳米碳化硅材料的制备及应用[Ｄ].太原:山西大学ꎬ２０１１.[６]㊀徐㊀威ꎬ兰㊀忠ꎬ彭本利ꎬ等.微液滴在不同能量表面上润湿状态的分子动力学模拟[Ｊ].物理学报ꎬ２０１６(２１):２１６８０１.[７]㊀王㊀龙.金属液滴在碳纳米材料表面的润湿与融合[Ｄ].济南:山东大学ꎬ２０１７.[８]㊀ＤｉｎｇＹꎬＺｈａｎｇＹꎬＺｈａｎｇＦꎬｅｔａｌ.ＭｏｌｅｃｕｌａｒｄｙｎａｍｉｃｓｓｔｕｄｙｏｆｔｈｅｓｔｒｕｃｔｕｒｅｉｎｖｉｔｒｅｏｕｓｓｉｌｉｃａｗｉｔｈＣＯＭＰＡＳＳｆｏｒｃｅｆｉｅｌｄａｔｅｌｅｖａｔｅｄｔｅｍｐｅｒａｔｕｒｅｓ[Ｃ].ＭａｔｅｒｉａｌｓＳｃｉｅｎｃｅＦｏｒｕｍ.２００７.[９]㊀ＢｅｅｓｔＢＷＨＶꎬＫｒａｍｅｒＧＪꎬＳａｎｔｅｎＲＡＶ.Ｆｏｒｃｅｆｉｅｌｄｓｆｏｒｓｉｌｉｃａｓａｎｄａｌｕｍｉｎｏｐｈｏｓｐｈａｔｅｓｂａｓｅｄｏｎａｂｉｎｉｔｉｏｃａｌｃｕｌａｔｉｏｎｓ[Ｊ].ＰｈｙｓｉｃａｌＲｅｖｉｅｗＬｅｔｔｅｒｓꎬ１９９０ꎬ６４(１６):１９５５.[１０]㊀ＳｕｎｄａｒａｒａｍａｎＳꎬＣｈｉｎｇＷＹꎬＨｕａｎｇＬ.Ｍｅｃｈａｎｉｃａｌｐｒｏｐｅｒｔｉｅｓｏｆｓｉｌｉｃａｇｌａｓｓｐｒｅｄｉｃｔｅｄｂｙａｐａｉｒｗｉｓｅｐｏｔｅｎｔｉａｌｉｎｍｏｌｅｃｕｌａｒｄｙｎａｍｉｃｓｓｉｍｕｌａｔｉｏｎｓ[Ｊ].ＪｏｕｒｎａｌｏｆＮｏｎ￣ＣｒｙｓｔａｌｌｉｎｅＳｏｌｉｄｓꎬ２０１６ꎬ４４５￣４４６:１０２￣１０９.[１１]㊀ＤｅｍｉｒａｌｐＥꎬＣａｇｉｎꎬＴａｈｉｒꎬＧｏｄｄａｒｄＷＡ.Ｍｏｒｓｅｓｔｒｅｔｃｈｐｏｔｅｎｔｉａｌｃｈａｒｇｅｅｑｕｉｌｉｂｒｉｕｍｆｏｒｃｅｆｉｅｌｄｆｏｒｃｅｒａｍｉｃｓ:ａｐｐｌｉｃａｔｉｏｎｔｏｔｈｅｑｕａｒｔｚ￣ｓｔｉｓｈｏｖｉｔｅｐｈａｓｅｔｒａｎｓｉｔｉｏｎａｎｄｔｏｓｉｌｉｃａｇｌａｓｓ[Ｊ].ＰｈｙｓｉｃａｌＲｅｖｉｅｗＬｅｔｔｅｒｓꎬ１９９９ꎬ８２(８):１７０８￣１７１１.[１２]㊀丁元法ꎬ张㊀跃ꎬ张凡伟ꎬ等.石英玻璃高温分子动力学模拟中的势函数[Ｊ].物理化学学报ꎬ２００８ꎬ２４(５):７８８￣７９２.[１３]㊀ＨｏｓｓｅｉｎｉＳꎬＳａｖａｌｏｎｉＨꎬＳｈａｈｒａｋｉＭＧ.Ｉｎｆｌｕｅｎｃｅｏｆｓｕｒｆａｃｅｍｏｒｐｈｏｌｏｇｙａｎｄｎａｎｏ￣ｓｔｒｕｃｔｕｒｅｏｎｈｙｄｒｏｐｈｏｂｉｃｉｔｙ:Ａｍｏｌｅｃｕｌａｒｄｙｎａｍｉｃｓａｐｐｒｏａｃｈ[Ｊ].ＡｐｐｌｉｅｄＳｕｒｆａｃｅＳｃｉｅｎｃｅꎬ２０１９ꎬ４８５:５３６￣５４６.[１４]㊀周书强.小口径非球面玻璃透镜的模压仿真及实验研究[Ｄ].长沙:湖南大学ꎬ２０１６.[１５]㊀ＨｅｒｄｅｓＣꎬＥｒｖｉｋＡꎬＭｅｊíａꎬｅｔａｌ.Ｐｒｅｄｉｃｔｉｏｎｏｆｔｈｅｗａｔｅｒ/ｏｉｌｉｎｔｅｒｆａｃｉａｌｔｅｎｓｉｏｎｆｒｏｍｍｏｌｅｃｕｌａｒｓｉｍｕｌａｔｉｏｎｓｕｓｉｎｇｔｈｅｃｏａｒｓｅ￣ｇｒａｉｎｅｄＳＡＦＴ￣γＭｉｅｆｏｒｃｅｆｉｅｌｄ[Ｊ].ＦｌｕｉｄＰｈａｓｅＥｑｕｉｌｉｂｒｉａꎬ２０１８ꎬ４７６:９￣１５.[１６]㊀ＲａｐｐｅＡＫꎬＣａｓｅｗｉｔＣＪꎬＣｏｌｗｅｌｌＫＳꎬｅｔａｌ.ＵＦＦꎬａｆｕｌｌｐｅｒｉｏｄｉｃｔａｂｌｅｆｏｒｃｅｆｉｅｌｄｆｏｒｍｏｌｅｃｕｌａｒｍｅｃｈａｎｉｃｓａｎｄｍｏｌｅｃｕｌａｒｄｙｎａｍｉｃｓｓｉｍｕｌａｔｉｏｎｓ[Ｊ].ＪｏｕｒｎａｌｏｆｔｈｅＡｍｅｒｉｃａｎＣｈｅｍｉｃａｌＳｏｃｉｅｔｙꎬ１９９２ꎬ１１４(２５):１００２４￣１００３５.[１７]㊀ＷｕＣꎬＣｈｅｎＧ.ＣｏｍｐｕｔａｔｉｏｎａｌｍｏｄｅｌｏｎｓｕｒｆａｃｅｔｅｎｓｉｏｎｏｆＣａＯ￣ＭｎＯ￣ＳｉＯ２ｓｌａｇｓｙｓｔｅｍ[Ｊ].ＨｏｔＷｏｒｋｉｎｇＴｅｃｈｎｏｌｏｇｙꎬ２０１４ꎬ４３(１９):５８￣６２. [１８]㊀ＢｅｌａｓｈｃｈｅｎｋｏＤＫꎬＯｓｔｒｏｖｓｋｉｉＯＩ.Ｃｏｍｐｕｔｅｒｓｉｍｕｌａｔｉｏｎｏｆｓｍａｌｌｎｏｎｃｒｙｓｔａｌｌｉｎｅｓｉｌｉｃａｃｌｕｓｔｅｒｓ[Ｊ].ＩｎｏｒｇａｎｉｃＭａｔｅｒｉａｌｓꎬ２００２ꎬ３８(９):９１７￣９２１. [１９]㊀王宝和ꎬ李㊀群.接触角的研究现状及其在凝胶干燥中的作用[Ｊ].干燥技术与设备ꎬ２０１４ꎬ１２(１):３９￣４６.[２０]㊀罗晓光.ＺｒＢ２/ＳｉＣ复合材料界面结构的分子动力学模拟[Ｄ].哈尔滨:哈尔滨工业大学ꎬ２００７.[２１]㊀ＬｉｕＢꎬＬａｎｇｅＦＦ.Ｐｒｅｓｓｕｒｅｉｎｄｕｃｅｄｔｒａｎｓｉｔｉｏｎｂｅｔｗｅｅｎｓｕｐｅｒｈｙｄｒｏｐｈｏｂｉｃｓｔａｔｅｓ:ｃｏｎｆｉｇｕｒａｔｉｏｎｄｉａｇｒａｍｓａｎｄｅｆｆｅｃｔｏｆｓｕｒｆａｃｅｆｅａｔｕｒｅｓｉｚｅ[Ｊ].ＪｏｕｒｎａｌｏｆＣｏｌｌｏｉｄ＆ＩｎｔｅｒｆａｃｅＳｃｉｅｎｃｅꎬ２００６ꎬ２９８(２):８９９￣９０９.[２２]㊀ＫａｌｉｎＭꎬＰｏｌａｊｎａｒＭ.Ｔｈｅｅｆｆｅｃｔｏｆｗｅｔｔｉｎｇａｎｄｓｕｒｆａｃｅｅｎｅｒｇｙｏｎｔｈｅｆｒｉｃｔｉｏｎａｎｄｓｌｉｐｉｎｏｉｌ￣ｌｕｂｒｉｃａｔｅｄｃｏｎｔａｃｔｓ[Ｊ].ＴｒｉｂｏｌｏｇｙＬｅｔｔｅｒｓꎬ２０１３ꎬ５２(２):１８５￣１９４.[２３]㊀叶㊀伟.强流脉冲电子束处理对热模压脱模性能的影响[Ｄ].重庆:重庆理工大学ꎬ２０１３.[２４]㊀ＷａｎｇＣꎬＫｕｚｕｕＮꎬＴａｍａｉＹ.Ｍｏｌｅｃｕｌａｒｄｙｎａｍｉｃｓｓｔｕｄｙｏｎｓｕｒｆａｃｅｓｔｒｕｃｔｕｒｅｏｆｈｏｔｗｏｒｋｉｎｇｔｅｃｈｎｏｌｏｇｙα￣ＳｉＯ２ｂｙｃｈａｒｇｅｅｑｕｉｌｉｂｒａｔｉｏｎｍｅｔｈｏｄ[Ｊ].ＪｏｕｒｎａｌｏｆＮｏｎ￣ＣｒｙｓｔａｌｌｉｎｅＳｏｌｉｄｓꎬ２００３ꎬ３１８:１３１￣１４１.。
聚丙烯酰胺稀溶液的分子模拟刘艳艳1陈攀科1罗健辉2周歌1,*江波1,*(1四川大学化学学院,绿色化学与技术教育部重点实验室,成都610064;2中国石油勘探开发研究院,北京100083)摘要：聚丙烯酰胺(PAM)是一类重要的线性水溶性聚合物,具有“百业助剂”之称,因此对其溶液性质的研究意义重大.在溶液质量浓度约为1g ·mL -1的基础上,分别构建了含有不同水分子数的溶液模型.采用分子动力学(MD)方法模拟分析了不同温度下非离子型的聚丙烯酰胺(PAM -H)和阴离子型的聚丙烯酰胺(HPAM)在纯水溶液及含不同质量分数NaCl 的水溶液中的回旋半径(R g ).结果发现,不同温度下PAM -H 和HPAM 的抗盐性能的模拟结果与实验数据基本吻合,水分子数不同的溶液模型所得模拟结果趋势没有明显变化,为了提高模拟效率,选取含有2000个水分子的溶液模型分析HPAM 链中氧负离子及氧原子的径向分布函数,从微观结构模拟说明了HPAM 水溶液粘度随NaCl 质量分数增加而减小,且HPAM 比PAM -H 具有较好的增粘效果及较差的抗盐性能的原因.关键词：分子动力学;聚丙烯酰胺;溶液性质;回旋半径;径向分布函数中图分类号：O645;O641Molecular Simulation of Dilute Polyacrylamide SolutionsLIU Yan -Yan 1CHEN Pan -Ke 1LUO Jian -Hui 2ZHOU Ge 1,*JIANG Bo 1,*(1Key Laboratory of Green Chemistry and Technology,Ministry of Education,College of Chemistry,Sichuan University,Chengdu610064,P.R.China ;2Research Institute of Petroleum Exploration and Development of PetroChina,Beijing 100083,P.R.China )Abstract ：Polyacrylamide (PAM)applied to various fields is an important class of linear water -soluble polymers.Therefore,it is of great significance to study the solution properties of PAM.We constructed solution models containingdifferent amounts of water molecules with a mass concentration of about 1g ·mL ing molecular dynamics (MD)simulations we calculated the radius of gyration (R g )for non -ionic PAM (PAM -H)and anionic PAM (HPAM)in pure water and in aqueous solutions with different mass fractions of NaCl.We discussed their behaviors at different temperatures.We found that the salt tolerance of the polyacrylamides from the simulation agreed with the experimental results at different temperatures.Furthermore,the simulation results for all the solution models containing a different amount of water molecules basically showed a similar trend.Considering computational efficiency,the solution model containing 2000water molecules was selected for our study.The radial distribution functions (RDF)for the oxygen ions and oxygen atoms of the HPAM chain were investigated in NaCl solution model containing 2000water molecules.The reduced viscosity of HPAM solutions with increasing NaCl mass fractions and a better thickening ability as well as poor salt tolerance compared to PAM -H were explained considering their microstructures as determined by RDF.Key Words ：Molecular dynamics;Polyacrylamide;Solution property;Radius of gyration;Radial distribution function[Article]物理化学学报(Wuli Huaxue Xuebao )Acta Phys.-Chim.Sin .,2010,26(11)：2907-2914聚丙烯酰胺(PAM)是一种线性的水溶性聚合物[1].由于其分子链的柔顺性和分子形状(即链构象)的易变性,以及具有高极性、易形成氢键和高反应活性的酰胺基,使得聚丙烯酰胺具有许多极有价值的November Received:July 5,2010;Revised:August 5,2010;Published on Web:September 16,2010.*Corresponding authors.JIANG Bo,Email:jiangbo@;Tel:+86-28-85418112.ZHOU Ge,Email:zhougekk@;Tel:+86-28-85418112.The project was supported by the National Natural Science Foundation of China (20904035).国家自然科学基金(20904035)资助项目鬁Editorial office of Acta Physico -Chimica Sinica2907Acta Phys.-Chim.Sin.,2010Vol.26应用性能,易通过接枝或交联得到支链或网状结构的多种改性.因此,聚丙烯酰胺类聚合物应用广泛,享有“百业助剂”之称.然而,对于最初的聚丙烯酰胺产品——丙烯酰胺均聚物(非离子型聚丙烯酰胺, PAM-H)和水解聚丙烯酰胺(阴离子型聚丙烯酰胺, HPAM),其性能不能完全满足应用需求.PAM-H易发生水解而具有阴离子型负电性,且它们分子上的—COO-对盐极其敏感,在高温、高矿化度、高剪切环境下,会导致水溶液粘度大幅度下降,不能适应高温、高矿化度抽藏驱油和钻井方面的要求.因此必须对最初的聚丙烯酰胺产品进行改性,来提高其耐温、耐盐、抗剪切性能.在聚电解质溶液中,离子基团的静电斥力引起链的扩张,使特性粘数和粘度增大;外加小分子电解质可屏蔽离子基团的静电作用,减小链尺寸,使特性粘数和粘度减小[1-3].因此,可以通过测定聚合物溶液的特性粘数或粘度来研究及预测丙烯酰胺类聚合物的性能(耐温、耐盐、抗剪切等). Yahaya等[4]引入N-苄基丙烯酰胺单体对聚丙烯酰胺进行疏水改性,同时测试其水溶液在不同高分子浓度、盐浓度、温度及剪切率下的表观粘度变化.结果发现,疏水基团在丙烯酰胺均聚物中的引入可以明显提高其增粘性,且随溶剂条件及剪切率的变化具有很好的稳定性,尤其在盐溶液中的特性粘数值要大于水溶液中的值.实际上,关于丙烯酰胺类聚合物改性及其溶液性能研究的实验文献国内外已有很多[2-7],从微观结构研究丙烯酰胺类聚合物的性质也备受关注[8-10],而这些传统的方法实验周期较长,且消耗较大的人力、物力、财力[11].计算机模拟技术[12]对高分子稀溶液的研究能够获得高分子溶液的热力学性质(高分子-溶剂体系的混合熵、混合热和混合自由能)、动力学性质(高分子溶液的粘度、高分子在溶液中的扩散和沉降)以及聚合物的分子量和分子量分布、高分子在溶液中的形态和尺寸、高分子链段间及链段与溶剂分子间的相互作用、高分子稀溶液的流变性质以及聚合溶液中的磁性等;同样也可以获得高分子浓溶液的性质.所以,我们可以用计算机模拟技术来研究聚合物分子结构在不同溶剂和温度下的变化,可以直观了解聚合物溶液宏观性质的本质,指导聚合物的改性及合成,缩短新材料的开发周期,弥补传统实验方法的不足,具有一定的环保意义.由于计算能力和研究尺度的限制,聚丙烯酰胺类聚合物稀溶液的模拟研究不多,但仍有一些关于高分子溶液性质模拟研究的文献报道[13-22],为我们进行此方面的研究提供了基础,如Tung等[18]通过分子动力学模拟技术研究了溶剂的种类对无规聚甲基丙烯酸甲酯(a-PMMA)链在溶液中的动力学性质的影响,构建了不同种类a-PMMA溶液,对均方末端距<R2>、均方回转半径<S2>以及聚合物链的非球面率<R2>/<S2>等聚合物链动力学性质进行了模拟分析.因此,本文进行了PAM-H和HPAM稀溶液的性质初探,分别构建了含不同水分子数的溶液模型,采用分子动力学(MD)方法首先研究了PAM-H和HPAM在不同水溶液模型中的回旋半径(R g)受NaCl 质量分数和温度的影响变化,从中确定稀溶液模型的体系,同时进行了径向分布函数(RDF)g(r)分析,从微观结构解释说明了PAM-H和HPAM的宏观性质,从而为以后更深入的研究、改性聚丙烯酰胺及其分子设计奠定基础.1计算方法的提出及模拟方法1.1特性粘数的理论计算特性粘数[η]指分子链在指定盐水中形成水动力学体积的特征参数.它与聚合物的结构、分子量大小和分子形态有关,在热、酸碱、盐等因素作用下,溶液中聚合物的分子尺寸会发生很大变化,从而导致聚合物溶液的特性粘数发生变化.理论上引入一个可以表征分子尺寸的参数———回旋半径(R g).特性粘数与回旋半径之间的关系如下[23-24]:[η]=63/2Φ(R g)3/M(1)式中,M为聚合物的摩尔质量,Φ=10πN A/3,N A为阿佛加德罗常数.本文通过Materials Studio(MS)4.2软件计算得到回旋半径(R g),由R g的变化趋势可看出[η]的变化趋势.1.2模拟方法模拟使用由美国Accelrys公司开发的Materials Studio4.2软件,采用amorphous cell模块的construction分别添加不同的NaCl分子数到充满一定数目的H2O分子和一个聚合物分子链的立方元胞中,construction中的目标密度(target density of the final configurations)是指模型的溶液密度,由模型含有的所有原子质量及其体积计算得到,立方元胞的体积由所含水的密度(1g·mL-1)和个数计算,模型体积元胞参数(cell parameters)随目标密度变化而自动生成,这样就得到了溶液模型的初始结构.关于溶液密度的确定,由于本文构建的是聚合物稀溶液模型,而实际稀溶液的密度近似认为接近2908No.11刘艳艳等：聚丙烯酰胺稀溶液的分子模拟水的密度(约为1g·mL-1),因此本文所设定的目标密度能够保证高分子链周围水的质量浓度约为1g·mL-1;同时由于构建的溶液模型浓度很低,其高分子链之间的作用力比高分子链与溶剂分子间的作用力小得多,因此可以忽略高分子链间的作用力,即构建模型时只选择一条高分子链是相对合理的.为了避免以后进行分子动力学模拟时,溶液模型的能量陷入低能势阱中而无法达到动力学平衡,需先用smartminimizer方法进行至少2次5000步能量优化(步长设为1fs),直到初始结构能量最小化,然后在一定温度(298、308、318、328或338K)下进行450ps分子动力学模拟,其中先进行50ps NPT(等温等压系综,系统的粒子数N、温度T和压力P恒定)动力学平衡(p=0.0001GPa),再进行400ps NVT(正则系综,系统的粒子数N、体积V和温度T恒定)分子动力学计算,然后对最后200ps分子动力学轨迹采用amorphous cell模块中的analysis来获得回旋半径R g,之后采用forcite模块中的analysis来进行径向分布函数分析.图1(a)显示了含原子数最多的溶液模型(298K时,NaCl质量分数为0.4%的HPAM-3500溶液模型,此溶液模型表示符号中,“-”前表示模型中的聚合物种类,“-”后表示模型中的水分子数,文中所有溶液模型符号均按此规则表示)在NPT系综下的能量随模拟时间的变化,在模拟时间约为50ps时溶液模型就达到平衡,因此上述选取50ps NPT平衡时间是合理的;然后取此模型50ps NPT动力学所保存的最后一个模型,继续采用NVT系综进行动力学计算所得的能量随模拟时间的变化始终保持平衡(如图1(b)所示）,这说明本文所采用的方法合理.在298K下的溶液模型分别进行了600ps、1ns和2 ns的分子动力学模拟,所得模拟结果与450ps的结果很相似.为了提高模拟工作效率,本文选取450ps 的分子动力学模拟.以上能量优化和分子动力学模拟是在PCFF力场[25-27]下运行的.Soldera[26]用PCFF 力场模拟分析了两条聚甲基丙烯酸甲酯链的立构规整度及聚丙烯酸甲酯的能量,Dai等[27]用PCFF力场模拟预测了甘油-水二元低温防护剂的玻璃化转变温度,两者都得到了较好的分子动力学模拟效果.而且PCFF力场[25]适用于聚合物、有机和无机材料、大约20种金属、碳水化合物、液体、核酸的内聚能、力学性质、压缩性质、热容、弹性常数等.因此,本文采用PCFF力场是合理的.同时,Bishop[28]和Pan[29]等已经研究了被溶剂包围的单链高分子溶液体系,模拟分析了低聚合度高分子链的动力学平衡性质,如相关函数、回旋半径、末端距、扩散系数等,结果发现可以通过低聚合度高分子链来研究实际高分子的重要性质.因此,本文首先分别构建一条聚合度为50的无规PAM-H分子链和一条聚合度为50的水解度为25%的无规HPAM分子链(如图2所示),之后将PAM-H分子链嵌入到含有1000、1500、2000、3000或3500个水分子的立方元胞中,并保证整个元胞成电中性的情况下,分别添加不同个数的Na+和Cl-.图3所示为298图1不同系综下,0.4%(w,质量分数)NaCl的HPAM-3500模型能量随模拟时间的变化Fig.1Energy of the model as a function of simulation time under different ensembles for the HPAM-3500 model with0.4%(w,mass fraction)NaCl y in HPAM-y model(or PAM-H-y model)is the number of watermolecules.图2PAM-H(a)和HPAM(b)的结构Fig.2Structures of PAM-H(a)and HPAM (b)2909Acta Phys.-Chim.Sin.,2010Vol.26K 时,含有2000个水分子的HPAM 溶液模型(HPAM -2000).所有构建元胞的离子或分子均采用smart minimizer 方法进行5000步的能量最小化.构建立方元胞时采用周期性边界条件(PBC)[30],即以立方元胞为中心,周围有26个相邻镜像立方元胞,以达到利用较少粒子进行宏观性质模拟的目的,这不仅使粒子所受到的力和宏观样品中本体粒子所受到的力一样,而且可以节省CPU 计算总时间.这是由于模拟体系越大,所需分子动力学模拟的CPU 总时间越长.由于实际溶液中相距较远的粒子之间相互作用很小,可忽略,因此,本文采用截断半径1nm,即忽略1nm 以外的粒子对其中心粒子的相互作用能.聚合度为50的水解度为25%的无规HPAM 分子链的溶液模型的构建及方法同上.本文对优化后的高分子溶液体系,用discover 模块运行分子动力学模拟,同时记录模拟轨迹.模拟温度分别为298、308、318、328和338K,恒温方法采用Andersen 方法[21,31],积分步长为1fs,即每隔1fs 各质点就会有新的坐标和速度,每隔0.5ps 保存一次分子轨迹和坐标,以便以后计算分析,且每个模型重复三次上述过程,所得分析结果为三次平均值(每组数据偏差在5%以内).2结果与讨论2.1盐质量分数和温度对PAM -H 和HPAM 回旋半径的影响分别通过对含1000、1500、2000、3000和3500个水分子的溶液模型所得回旋半径值的比较,发现随水分子数的增加,不同模型所得结果趋势没有很大的偏差,本文从中选取具有代表意义的含1000或3500个水分子的溶液模型进行讨论(如表1所示).2.1.1盐对PAM -H 和HPAM 回旋半径的影响线性聚合物的回旋半径是线性高聚物的特征参数,它直接反映线性高分子链的构象,而链的构象又影响聚合物固态和溶液的许多性质,因此用计算机模拟高分子回旋半径是一项有意义的工作.表1列出本文所构聚合物溶液模型经能量优化、动力学模拟以及分析后的各参数及不同温度下的聚合度均为50的PAM -H 和HPAM 在不同溶液模型中的回旋半径值.由表1发现,各溶液模型中,当没有盐(NaCl)存在时,在不同温度下各高分子链的回旋半径值基本是最大的.即有盐存在情况下,普通的聚丙烯酰胺类产品的回旋半径减小(尤其是HPAM),且随盐质量分数的增加聚丙烯酰胺类产品的回旋半径基本依次减小,HPAM 的回旋半径减小趋势更显著,由方程(1)可知,[η]与R g 3/M 成正比,即特性粘数也随盐质量分数的增加而减小;在没有盐存在情况下,HPAM的回旋半径明显比PAM -H 的高.据实验数据[1]可得,普通的聚丙烯酰胺类产品的主要水化基团—COO -抗盐能力较差,在水溶液中由于它的水化作用和静电排斥作用,使分子主链伸展,有着较大的回旋半径,故特性粘数较大;而在盐溶液中,由于盐离子特别是高价盐离子的存在很大程度上容易破坏水化层,放出自由水,减弱了—COO -基团沿分子主链的排斥作用,使其分子主链蜷曲,回旋半径减小,特性粘数明显下降[1-3].综上所述,模拟的HPAM 回旋半径在盐存在情况下是明显减小的,且随盐质量分数的增加回旋半径减小,即此现象的模拟结果与实验结果[1-3]基本相吻合;且在没有盐存在情况下,HPAM 的回旋半径明显比PAM -H 的回旋半径高,这是由于模拟的HPAM 中含有—COO -基团,而模拟的PAM -H 中不含有此基团,因此模拟的HPAM 的回旋半径比PAM -H 大,其溶液粘度比PAM -H 的溶液图3未含有(a)和含有(b)2%NaCl 的HPAM -2000溶液模型Fig.3HPAM -2000solution models without (a)and with (b)2%NaCl(a)(b)2910No.11刘艳艳等：聚丙烯酰胺稀溶液的分子模拟粘度高,即同一分子量范围,HPAM比PAM-H具有更好的增粘效果,但HPAM抗盐性能比PAM-H差,此模拟结果同样与实验结果相一致[1].2.1.2温度对PAM-H和HPAM回旋半径的影响据实验资料[1]可知,普通PAM类产品的溶液粘度随温度的升高而降低,这是由于温度升高使分子链降解,其中包括高分子主链的断裂和主链与侧基连接键的断裂,前者使产品的分子量降低,部分或全部失去高分子性质,后者降低产品的亲水性,然而这需要的活化能比较高,一般200-220℃以上才发生[1].然而对于没有发生降解的PAM产品,在溶液中由于温度升高容易引起高分子主链单键的内旋转而引起构象的变化,从而导致高分子链尺寸增大,回旋半径增大,分子量越大回旋半径增加得越显著,相反,分子量越小回旋半径增加得越平缓[23-24].总体而言,在没有降解的情况下(200℃以下),PAM-H和HPAM溶液的粘度随温度的增加会有所增加.而对于模拟情况,考虑到高分子链节数只有50,相对实际PAM-H产品分子量极小,回旋半径随温度的增加会有所增加,但会非常平缓,如表1中所示,我们发现,各溶液模型中的高分子链的回旋半径随温度的变化不是很稳定,但数据点基本平均分布在趋势线的两侧,因此,图4中均用趋势线来描述数据的发展趋势.从图4中发现,有些系列的回旋半径总体趋向是随温度的增高而增大,如图4(a,c,d)中所示;而有些系列的回旋半径总体趋势是随温度的升高而降低,如图4(a,b,d)中所示.而上述实验结论[23-24]是回旋半径随温度的升高应平缓增加,即模拟数据与实验结论有些偏差,这意味着:如果想得到精确的变化趋势平缓的模拟数据(如此处回旋半径随温度的变化数据),需优化模拟方法,这将在以后更加系统地研究PAM-H和HPAM溶液性能时进行考察.综上所述,对于不同溶液模型,水分子数增加没有明显改变模拟结果的趋势,但是从含有1000个水分子到含有3500个水分子的溶液模型,其所含有的原子个数从3502增加到11002,而原子个数越多,模拟所需时间越长,从考虑计算效率的角度来看,我们应考虑在模拟结果趋势没有明显改变的前提下,选取原子数相对较少的模型.对于含有1000或1500个水分子的溶液模型,在以后研究pH值对PAM-H 和HPAM溶液的影响时,按实际比例所能在模型中添加的H+、OH-的个数太少,对于含有3000或3500个水分子的溶液模型,其原子个数达到了104数量级,体系太大,所需时间太长.综合考虑,本文确定了溶液模型中水分子数为2000.2.2HPAM溶液模型的径向分布函数分析由上述可知,PAM-H和HPAM的回旋半径随盐质量分数的增加会有不同程度的减小,从而导致溶液粘度减小.本节将以2000个水分子的模型为例,从微观结构验证上述观点.表1PAM-H和HPAM各溶液模型的参数及其在不同温度下的回旋半径值(R g)Table1Parameters for solution models and radius of gyration(R g)for PAM-H and HPAM at2911Acta Phys.-Chim.Sin.,2010Vol.26图5298K HPAM -2000模型中不同NaCl 质量分数下O --H 2O(a)和O -H 2O(b)径向分布函数图Fig.5Radial distribution functions of O --H 2O (a)and O -H 2O (b)of HPAM -2000model with different massfractions of NaCl at 298K径向分布函数g (r )是反映流体微观结构特征的物理量[32-33],其表示在距离某一设定中心原子A 为r 处,某一设定原子B 的数目密度与B 的平均数目密度的比值,即:g (r )=d N ρ4πr2(2)式中,d N 表示与中心的距离为r →r +d r 间的分子数目,ρ表示原子B 的平均数目密度.已有研究者[32-33]通过径向分布函数来研究离子水化作用,即通过离子-水径向分布函数的第一峰的大小来判断离子水化的强弱,且一般采用阴离子与水中的氢做径向分布函数.本文通过径向分布函数研究了在298K 下的不同NaCl 质量分数(w =0%、1%、2%)的HPAM 溶液模型微观结构.图5(a)和5(b)分别是溶液模型中w Na C l(%)w Na C l(%)图4不同NaCl 质量分数下各溶液模型中的回旋半径-温度关系图Fig.4Radius of gyration for PAM -H or HPAM as a function of temperature with different mass fractions of NaClsolution model:(a)PAM -H -1000;(b)HPAM -1000;(c)PAM -H -3500;(d)HPAM -35002912No.11刘艳艳等：聚丙烯酰胺稀溶液的分子模拟高分子链上的氧负离子(O-)-水径向分布函数g O--H(r)和氧原子(O)-水径向分布函数g O-H(r).从图5(a)可知,不同NaCl质量分数下,HPAM分子链上的O--H2O 径向分布函数均在0.16nm处出现第一峰值,其峰值随NaCl质量分数的增加分别为1.46、1.37、1.09,这说明NaCl加入到HPAM的水溶液中,确实破坏了高分子链中O-的水化层,放出自由水,有去水化作用,且随着NaCl质量分数的增加去水化作用增强,从而减弱了—COO-基团沿分子主链的排斥作用,使其分子主链蜷曲,回旋半径依次减小(如图6所示),从图3中也可以看出HPAM尺寸变化;而对于高分子链中O的水化作用可以从图5(b)看出,O-H2O径向分布函数均在0.18nm处出现第一峰值,其峰位置相对O--H2O径向分布函数发生右移,且其峰值随NaCl质量分数的增加分别为0.571、0.569、0.561,这表示O的水化作用很弱,与O-相比可以忽略不计,且NaCl质量分数的增加对其影响不大,这也从微观角度说明了上述HPAM比PAM-H具有较好的增粘效果及较差的抗盐性能的原因.3结论在保证溶液模型中高分子链周围水分子的质量浓度约为1g·mL-1的基础上,构建含有不同水分子数的PAM-H和HPAM溶液模型,MD方法模拟计算了其在不同温度和不同NaCl质量分数下的回旋半径,还模拟分析了部分模型的径向分布函数,结果发现,不同温度下PAM-H和HPAM的回旋半径随NaCl质量分数变化的模拟结果与实验结果在趋势上基本一致,且随着水分子数的增加,各溶液模型的模拟结果趋势没有明显的改变,考虑到模拟工作效率,我们确定了含有2000个水分子的溶液模型体系.同时,对含有2000个水分子的溶液模型作了进一步研究,模拟分析了HPAM在不同盐质量分数水溶液的径向分布函数,从微观角度分析了HPAM水溶液粘度随盐质量分数增加而减小,且HPAM比PAM-H具有较好的增粘效果及较差的抗盐性能的内在原因,这也进一步证明了我们研究工作的合理性.因此,在以后的研究工作中,我们将在此基础上,进一步完善模拟方法,尤其对变化趋势平缓的性能研究,如模拟不同高分子回旋半径随温度的变化,需进一步优化模拟方法,从而对PAM-H和HPAM溶液进行更准确、更深入的研究,使模拟结果能更好地说明解释实验现象,以便更好地指导实验工作.References1Fang,D.B.;Guo,R.W.;Ha,R.H.Acrylamide polymer.Beijing: Chemical Industry Press,2006:1-86[方道斌,郭睿威,哈润华.丙烯酰胺聚合物.北京:化学工业出版社,2006:1-86]2Gong,H.J.;Xin,X.;Xu,G.Y.;Wang,Y.J.Colloids Surf.A-Physicochem.Eng.Asp.,2008,317:5223Fang,D.B.;Guo,R.W.;Zhou,S.G.;Zhang,X.C.CIESC Journal, 1966,47:645[方道斌,郭睿威,周少刚,张曦晨.化工学报,1966,47:645]4Yahaya,G.O.;Ahdab,A.A.;Ali,S.A.;Abu-Sharkh,B.F.;Hamad,E.Z.Polymer,2001,42:33635Volpert,E.;Selb,J.;Candau,F.Polymer,1998,39:10256Ye,L.;Luo,K.F.;Huang,R.H.Eur.Polym.J.,2000,36:17117McCormick,C.L.;Nonaka,T.;Johnson,C.B.Polymer,1988,29: 7318Zhang,Y.B.;Wu,C.;Fang,Q.;Zhang,Y.X.Macromolecules, 1996,29:24949Klucker,R.;Munch,J.P.;Schosseler,F.Macromolecules,1997, 30:383910Qiu,X.P.;Zhang,X.R.;Ding,M.T.Chem.J.Chin.Univ.,1995, 16(8):1321[邱星屏,张雪蓉,丁马太.高等学校化学学报,1995,16(8):1321]11Zhu,W.P.Plastics Science and Technology,2002,(5):23 [朱伟平.塑料科技,2002,(5):23]12Chen,P.F.;Zhang,L.H.China Adhesives,2007,6:52[陈鹏飞,张丽华.中国胶黏剂,2007,6:52]13Yuan,S.L.;Xu,G.Y.;Cai,Z.T.Colloid Polym.Sci.,2003,281(1):6614Zeng,Q.H.;Yu,A.B.;Lu,G.Q.Prog.Polym.Sci.,2008,33(2): 19115Dalakogloua,G.K.;Karatasosa,K.;Lyulin,S.V.;Lyulin,A.V.图6298K时,HPAM-2000模型中HPAM回旋半径-NaCl质量分数关系图Fig.6Radius of Gyration for HPAM as a function ofthe mass fractions of NaCl at298K in the HPAM-2000model2913Acta Phys.-Chim.Sin.,2010Vol.26Mater.Sci.Eng.B,2008,152:11416Le,T.C.;Todd,B.D.;Daivis,P.J.;Uhlherr,A.J.Chem.Phys., 2009,130:07490117Drew,P.M.;Adolf,D.B.Soft Matter,2005,1:14618Tung,K.L.;Lu,K.T.;Ruaan,R.C.;Juin,J.Y.Desalination, 2006,192:38019Kairn,T.;Daivis,P.J.;Matin,M.L.;Snook,I.K.Int.J.Thermophys.,2004,25(4):107520Liu,Y.M.;Li,G.Z.;Song,W.K.;Wang,J.J.Acta Phys.-Chim.Sin.,2006,22(12):1456[刘永明,李桂芝,宋万坤,王进军.物理化学学报,2006,22(12):1456]21Tao,C.G.;Feng,H.J.;Zhou,J.;Lü,L.H.;Lu,X.H.Acta Phys.-Chim.Sin.,2009,25(7):1317[陶长贵,冯海军,周健,吕玲红,陆小华.物理化学学报,2009,25(7):1317]22Shao,Q.;Lü,L.H.;Lu,X.H.;Wei,M.J.;Zhu,Y.D.;Shen,W.F.Acta Phys.-Chim.Sin.,2009,25(3):583[邵庆,吕玲红,陆小华,魏明杰,朱育丹,沈文枫.物理化学学报,2009,25(3):583] 23Rubinstein,M.;Colby,R.Polymer physics.Oxford:Oxford University Press,2002:Chapter2724He,M.J.;Chen,W.X.;Dong,X.X.Polymer physics.Shanghai: Fudan University Press,2007:15-73[何曼君,陈维孝,董西侠.高分子物理.上海:复旦大学出版社,2007:15-73]25Sun,H.;Mumby,S.J.;Maple,J.R.;Hagler,A.T.J.Am.Chem.Soc.,1994,116:297826Soldera,A.Polymer,2002,43:426927Li,D.X.;Liu,B.L.;Liu,Y.S.;Chen,C.L.Cryobiology,2008, 56:11428Bishop,M.;Kalos,M.H.;Frisch,H.L.J.Chem.Phys.,1979,70: 129929Pan,R.;Liu,X.K.;Zhang,A.M.;Gu,p.Mater.Sci.,2007, 39:88730Yang,X.Z.Molecular simulation and polymer materials.Beijing: Science Press,2002:43[杨小震.分子模拟与高分子材料.北京:科学出版社,2002:43]31Andersen,H.C.J.Chem.Phys.,1980,72:238432Zhou,J.;Zhu,Y.;Wang,W.C.;Lu,X.H.;Wang,Y.Y.;Shi,J.Acta Phys.-Chim.Sin.,2002,18(3):207[周健,朱宇,汪文川,陆小华,王延儒,时均.物理化学学报,2002,18(3):207] 33Zhou,J.;Lu,X.H.;Wang,Y,R.;Shi,J.CIESC Journal,2000,51: 143[周健,陆小华,王延儒,时均.化工学报,2000,51:143]2914。
聚合物粘结剂PVP、PAM和PVA结构与性能的分子动力学研究柴卫红;汪焰恩;魏庆华;杨明明;张映锋;魏生民【摘要】在粉末打印骨支架的工艺中,粘结剂的性质是影响骨支架质量的关键因素.采用分子动力学的模拟方法对三种常用高分子粘结剂的体系进行了构建和模拟,从微观分子层面研究了聚合物粘结剂PVP,PAM和PVA的部分性质,比如密度、内聚能及力学性能,并对三种粘结剂的性能参数进行了比较,对其内在关系进行了揭示.此外,还通过建立粘结剂与羟基磷灰石的界面相互作用模型,对三种粘结剂与羟基磷灰石的界面结合能进行了计算和比较,分析了影响高聚物粘结特性的根本原因.这一工作不仅对常用粘结剂的基本性质进行了预估,而且对骨支架粉末粘结工艺中粘结剂的选择提供了理论依据.%During the process of bone scaffold manufacturing with 3D printing technology,the properties of binder is a key factor that affects the quality of bone scaffold.In this paper,a molecular dynamics simulation method was applied to build and simulate the models of three commonly used polymer binders,and some properties of polymer binders PVP,PAM and PVP are investigated from a microscopic molecularlevel,such as density,cohesive energy and mechanical properties.The performance parameters of three binders are compared and studied,and then the internal relations are also revealed.Moreover,the interaction models of binders and HA surfaces are constructed,and the interfacial bonding energies are calculated and compared respectively,and the root reason affecting polymer bonding properties is also found.This work not only predicts the basic properties of commonly used binder,but alsoprovides the theoretical basis for the choice of binder for the bone scaffold manufacturing by 3D printing technology.【期刊名称】《计算力学学报》【年(卷),期】2018(035)003【总页数】6页(P344-349)【关键词】聚合物粘结剂;分子动力学;内聚能密度;结合能;机械性能【作者】柴卫红;汪焰恩;魏庆华;杨明明;张映锋;魏生民【作者单位】西北工业大学机电学院,西安710072;西北工业大学机电学院,西安710072;西北工业大学机电学院,西安710072;西北工业大学机电学院,西安710072;西北工业大学机电学院,西安710072;西北工业大学机电学院,西安710072【正文语种】中文【中图分类】O647.11;O3131 引言随着三维打印技术的发展，基于微滴喷射原理的三维粉末粘结打印技术已在很多领域得到应用，特别是在骨支架制备方面技术已比较成熟[1,2]。
