热力学模拟 篇1
钨和钨合金具有高强度、高密度, 良好的抗腐蚀性、导电、导热性和热膨胀系数小等一系列优良物理、力学性能。在军事领域中, 如:卫星、飞行器、航空发动机、核潜艇、导弹和舰艇等装备的一些关键部件大量采用了钨及钨合金, 它们的机构设计和安全性能问题与所采用的材料物理化学和机械特性有着十分密切的关系。
国内外许多研究学者对钨及钨合金进行了大量的宏观实验研究, Zhou等[3]采用拉伸冲击实验装置对钨合金试件在不同应变率下的工况进行了实验研究, 且在对冲击拉伸后断裂试件扫描分析, 建立了含有应变率相关参数的应力-应变模型。刘海燕等[4]基于旋转盘式间接杆-杆型冲击拉伸装置, 对不同颗粒度细化钨合金材料在不同应变率工况进行了实验研究, 通过分析钨合金材料在弹性阶段存在非线性, 判断了在此阶段材料内部已经损伤、钨颗粒减小材料屈服极限提高以及细化钨合金材料具有应变率致脆现象。李娜等[5]通过实验研究了83W和89W材料各向异性在不同温度环境的动、静加载条件下, 应力-应变关系及失效问题。Pink等[6,7]研究了钨合金材料在不同温度下的拉伸性能, 并发现钨合金材料在一定温度下出现动态应变时效现象和负应变率敏感性问题。王世良等[8]研究了钨合金的合成以及相应的物理问题, 谷怀鹏等[9]研究了镍基单晶合金物理特性问题。王建伟等[10]采用分子动力学方法研究bcc-Fe不同空位浓度对中子辐照损伤的影响。以上这些专家与学者分别对钨合金材料的实验研究、相关合金物理特性和金属材料的分子动力学模拟等方面均取得了一定的研究成果, 并促进钨及其合金材料的进一步发展。
然而, 利用扫描电子显微镜对钨合金材料和钨合金结构破坏观察研究分析, 钨合金材料主要由钨颗粒构成, 而在不同加载工况下, 钨颗粒的破坏是造成钨合金整体失效的主要原因。钨颗粒的强度决定了钨合金的物理力学性能, 但是由于钨的高强度、高硬度, 在宏观试验中很难测得其实际的物理性能的常规参数。针对这一问题, 采用分子动力学方法研究钨颗粒的材料力学性能成为了一种有效的手段, 钨颗粒尺寸对宏观材料尺寸相对无限小, 而对微纳米尺寸又可认为是无限大, 分子动力学方法它可以根据原子间的相互作用计算整个原子系统的性能, 分析材料的纳观变形机理, 再现宏观力学行为, 同时又提供大量的微观信息。
分子动力学模拟是一种用来计算一个经典多体体系的平衡和传递性质的方法, 是指对于原子核和电子组成的多体系统, 应用牛顿力学运动方程, 对原子核运动过程进行模拟, 记录运动过程中的轨迹, 最终利用统计力学方法计算一些相关物理量, 得出物质结构特性的一种理论计算方法。该模拟手段已成功应用于研究晶体熔化、晶体冲击、晶格畸变等方面[11], 并且随着计算方法的发展以及计算机性能的提高, 将被广泛应用于生物、化学、医药等领域, 有着广阔的应用前景。采用分子动力学方法, 对钨合金的拉伸塑性问题还未展开系统研究。
利用分子动力学模拟钨合金拉伸塑性力学行为, 对[100]晶向进行温度效应、尺度效应、应变率效应的模拟研究, 得到了钨合金微观结构变化图像和应力应变关系等相关塑性力学数据, 以及解决了在微纳米尺度材料实验中较难观测原子轨迹及提取相应数据测量的重要问题, 对于进一步了解微纳米尺度钨合金材料力学性能及其材料制备与应用提供了重要的依据。
1 钨合金的拉伸模型
1.1 嵌入原子势
嵌入式原子势 (EAM) 是Daw和Baskes[12]根据Hohen-kohn定理建立的势函数模型, 具有运算量小、运算效率高、适合金属间的作用等优点。
EAM原子势的总能量E可以表示为:
式中:φij表示距离为rij的原子i和原子j之间的对势能;Fi表示电子密度为ρi的原子i嵌入能;Mi (Pi) 为修正项;Pi为原子i处电子密度中高次项的贡献。
1.2 钨合金计算模型
图1为金属体心立方结构示意图, 晶格常数为a=b=c=0.3165nm[13]。图2是以体心立方结构为基础, 建立钨合金圆柱形超晶胞计算模型 (钨合金主体钨为体心立方结构) , 初始构型按理想晶格排列, x, y, z坐标轴分别对应[100], [010], [001]方向。圆柱半径为r=14a, 长为l=40a, 两端分别固定长度h=2a为夹具区域, 模型左端固定, 右端采用“velocity”加载方式进行拉伸, 共有49773个原子。为模拟有限长度的力学行为, x, y, z方向均取自由边界条件。
2 模拟方法
采用Ackland的EAM势函数[14]描述钨合金原子之间的相互作用, 它能准确的拟合出钨合金结合能、缺陷形成能、层错能等参数, 更能准确拟合钨合金的弹性性质。首先, 本模型采用Hoover-Nose控温系统, 在正则系宗 (Canonical Ensemble, NVT) , 即表示具有确定的粒子数 (N) 、体积 (V) 、温度 (T) 环境下, 对模型体系进行弛豫, 将模型控制在一定的温度使体系达到一个初步稳定的状态。再利用微正则系宗 (Microcanonical Ensemble, NVE) 。即表示具有确定的粒子数 (N) 、体积 (V) 、总能量 (E) 对体系进一步弛豫, 达到理想的加载环境。在此基础上, 采取与宏观实验相类似的加载方式, 对钨合金圆柱模型进行拉伸加载问题的分子动力学计算模拟。
3 钨合金体系模型弛豫
所谓模型弛豫, 就是体系达到一个稳定的状态, 以便获得准确的数据。图3为[100]晶向、截面半径r=14a、长度l=40a、温度300K下模型弛豫的动能与势能的变化情况, 由图3分析可知, 在Hoover-Nose控温系统下, 在体系进行弛豫初始时, 体系动能与势能上下波动幅度较大, 此时体系并有达到平衡状态。在经过30ps的弛豫以后, 体系动能在2000eV小幅波动变化, 而势能在-433000eV左右小幅度周期性波动, 体系基本上达到一个相对稳定的状态, 可以对其进行加载计算。
4 数值模拟与实验验证
4.1 钨合金塑性变形数值模拟研究
4.1.1 在拉伸过程中的体系结构变化
采用[100]晶向、截面半径r=14a、长度l=40a、初始温度300K。首先在NVT系宗下弛豫50ps, 再在NVE系宗下弛豫50ps。对模型右端加载区域施加0.25nm/ps沿x轴正向速度进行拉伸, 中间区域为自由运动, 每步输出一次热力学信息 (温度、压力、势能、总能量、x轴应力、x轴应变) 共运行100000步。图4为对钨合金试件拉伸加载作用下, 原子运动变化趋势 (基于OVITO软件[16]对计算数据进行图形绘制) 。
当0≤t<7.5ps时, 在加载速率初始阶段, 钨合金为完整晶体, 原子结构为规则的晶格排列, 拉伸变形为弹性变形。随着晶体拉伸持续进行, 原子之间产生了位错, 并伴随着滑移变形, 在模型表面出现轻微滑移台阶。
当7.5≤t<15ps时, 钨合金在拉力的作用下不断伸长, 发生了由弹性变形到塑性变形的转变。从图中可看出原子在 (110) 晶面、沿[111]晶向进行滑移, 随着拉伸的进行, 原子的运动越来越快, 开始出现部分金属键发生断裂, 当应变不断增加时, 位错不断增殖, 滑移持续发生, 中间段横截面发生改变并不断的减小, 并出现了颈缩现象。当拉伸进一步伸长时, 塑性变形主要集中在颈缩部位, 其他部分原子位置基本不变, 颈缩区域不断伸长, 导致滑移系的连续发生, 颈缩区域横截面不断减小, 颈缩现象更加明显, 大部分金属键已发生断裂。
当15≤t<28.5ps时, 随着拉伸的进一步进行, 颈缩区域原子明显减少, 原子相互作用力越来越小, 当金属键完全断裂, 晶体发生断裂, 表明纳米尺度下钨合金高应变率拉伸变形方式类似于宏观体心立方材料的拉伸变形。
4.1.2 钨合金高应变率拉伸过程变化曲线
计算原子i的应力分量采用如下公式:
式中:m表示原子质量;a, b表示x, y, z轴对称张量的六个分量;r表示两个原子之间的距离;F表示两个原子之间的作用力。第一项表示原子i的动能, 第二项是原子i周围Np近邻原子数经过n次循环的对能, 第三项是原子i周围Nb的键能, 第四项是原子i周围Na键角对应力的作用, 第五项是原子i周围Nd面角对应力的作用, 第六项是原子i周围Ni非健项对应力的作用, 最后一项是原子i内部约束力Nf对应力的作用。
图5为不同工况拉伸加载下铝合金应力-应变曲线。从图5 (a) 可以看出, 分子动力学模拟与宏观实验得到的应力-应变曲线相同的是初始为线性上升, 应力与应变成正比例关系, 材料为弹性变形阶段, 区别是分子动力学模拟得到的弹性极限到最高点后发生突降现象, 并产生塑性流动。钨合金模型随着加载速率的增加, 屈服强度逐渐增大。钨合金在0.010, 0.050, 0.100nm/ps的高应变率加载速率拉伸过程中发生了二次屈服, 而在0.001nm/ps时钨合金模型没有产生二次屈服。
在分子动力学模拟中, 不同拉伸速率下的钨合金模型采用了相同的弛豫时间, 因此, 在较低的拉伸速率下, 弛豫时间充分, 钨合金可以再次达到稳定的平衡状态, 而在较高的拉伸速率下, 弛豫时间不足以使钨合金体系达到稳定的平衡状态。因此, 在相对较高的拉伸速率下, 钨合金体系会产生二次屈服, 而在较低的拉伸速率下没有发生二次屈服。
而MD模拟弹性极限到达最高点后发生突降现象, 是由于在拉伸过程中原子之间持续产生位错, 进一步发生了滑移, 当达到钨合金屈服应力峰值时, 原子之间滑移不断加剧, 引起原子之间大量金属键发生断裂, 故导致“突降”问题。
从对图5 (b) 应力-应变曲线分析可知, 钨合金模型在温度300K、长度l=14a、[100]晶向以0.1nm/ps速率加载的条件下, 随着模型横截面半径的增大, 在弹性阶段钨合金的屈服强度也越大。这说明在高应变率拉伸过程中, 横截面的越大, 在单位时间内产生的位错密度越大, 位错之间的交互作用也就越强, 增加了位错运动的阻力, 导致了屈服强度的提高。
由图5 (c) 可以看到, 在温度300K、横截面半径r=14a、[100]晶向以0.1nm/ps加载的条件下, 随着钨合金模型长度的增加, 钨合金在弹性阶段屈服强度在增大。从图中可以看出, 不同长度钨合金模型的应力-应变曲线斜率相等, 说明弹性模量与模型长度尺寸没有关系。
从对图5 (d) 的应力-应变分析可知, 钨合金塑性变形对温度的变化具有很高的敏感性, 但应力-应变变化的趋势基本相同。首先, 是线性上升, 为钨合金的弹性变形阶段。其次, 随着应变的增加, 体系进入塑性变形阶段, 且体系应力呈现一定的波动性。而随着温度的上升, 其屈服强度逐步下降, 屈服后应力波动幅度逐渐降低, 与宏观拉伸变形相似。
这是因为, 随着温度升高, 金属原子振动更加剧烈, 从而金属键合力减小, 金属发生软化现象, 当温度升高到一定程度时, 金属会发生一定的黏性流动, 且温度越高, 则黏性流动的抗力越小, 其发生蠕变的范围也越大, 蠕变很大程度上导致金属屈服应力等非弹性变形对应的流动应力下降, 但提高了钨合金的塑性。
图5 (e) 表明, 不同晶向的高应变率拉伸加载钨合金材料具有不同的等效弹性刚度和屈服强度, [111]晶向钨合金模型屈服强度最大, [110]晶向次之, 而[100]晶向钨合金模型屈服强度最低。[100]晶向的钨合金呈现最大屈服应变, [111]与[110]晶向的屈服应变近似, 可近似为[100]晶向的三分之一。晶向对力学性能的影响这一特性对于纳米材料在不同领域的应用具有极其重要的研究价值。
4.2 钨合金拉伸实验验证
4.2.1 准静态拉伸验证
实验材料选取91钨合金细化颗粒材料为对象, 进行了不同工况的实验研究, wW∶mNi∶mFe=91.0∶6.3∶2.7。它是将10.0~15.0μm的细化预合金粉末经固相烧结得到的钨含量相同晶粒度不同的细观粒状复合材料。图6为将烧结后钨合金材料 (棒材) 经线切割加工后成扁平哑铃型试件。在实验前需要把试件实验段的圆弧部分沿纵向研磨光滑, 以减小试件拉伸段圆弧部分的应力集中现象。
图7为采用S-570扫描电子显微镜试验系统对91钨合金细化颗粒材料试件扫描后, 给出的钨合金微观组织像。由图7可知, 钨合金钨颗粒为半径不同, 但近似球形, 相对均匀的被黏结相包围。颗粒与颗粒之间的界面接触面较小, 黏结相和晶粒之间的界面在所有界面中占多数, 还可以看到有孔洞等微缺陷。
准静态单向拉伸实验在MTS810液压伺服材料实验机上进行, 采用位移控制速率方式加载, 在试件实验段安装引伸计测量应变, 通过计算机输出实验结果数据文件。表1给出了不同应变率下的钨合金材料准静态屈服强度。
图8为91在不同应变率下, 钨合金材料准静态拉伸实验的应力-应变曲线。由图8分析可知, 随着拉伸应变率的增加, 91钨合金材料的强度有一定的提高, 具有典型的应变率强化效应;在不同的拉伸应变率加载条件下, 钨合金材料应力-应变曲线主要呈现了弹性变形阶段, 且材料伸长率达到0.153%时, 试件发生断裂, 几乎没有塑性变形, 据此可以判断此颗粒度的细化钨合金材料为典型脆性材料。
对准静态91钨合金细化颗粒材料在准静态拉伸加载破坏的试件材料断口利用电镜扫描图像 (见图9) 。对图像观测分析, 可以看到大部分钨颗粒发生解理断裂, 颗粒的解理断面上是扇形花样, 解理台阶为扇形的肋。扇形花样的肋线的汇合点, 为解理裂纹的起源点。河流花样以扇形的方式向外扩展形成所谓的扇形花样。从钨合金材料断口处还可以观察到相邻颗粒的接触面在受外载时分开后留下的痕迹。说明钨合金材料中同样存在少量颗粒与颗粒之间的W-W界面。
从图9中钨合金材料断口的显微图像上较为光亮的部分是材料的黏结相, 黏结相断口呈现纤维状, 有韧窝, 是韧性断裂。钨合金材料在外受载时, 黏结相协调了变形并重新分配应力, 使强度较高的钨颗粒承担了主要载荷;继续增大外载到一定值, 钨颗粒之中产生微裂纹, 裂纹扩展, 钨颗粒被解理。钨颗粒的解理断裂是材料失效的主要原因。
4.2.2 动态拉伸实验
图10为旋转盘式间接杆-杆型冲击拉伸实验加载装置示意图。采用该装置对钨合金试件进行加卸载实验, 其中加载装置为直径1.4m, 转动线速率可以达到100m/s的飞轮。
91细化钨合金材料动态拉伸实验是在室温 (14℃) 条件下进行的, 分别在1×10-3, 200s-1和500s-1三种应变率工况下进行动态拉伸实验研究, 在每一工况下进行了四组实验。图11为在不同应变率条件下实验获取的钨合金材料动态拉伸的应力-应变关系曲线。
对图11分析可知, 不同应变率动态拉伸导致了钨合金材料应力-应变曲线有一定的差异, 反映了91钨合金颗粒材料具有应变率效应, 且随着加载速率的提高, 屈服应力相应增加, 但是由于材料自身的缺陷, 导致了测量的数值具有一定波动, 也降低了钨合金材料强度的下降。
为了更好地了解钨合金材料的力学性能及在高应变率下钨合金强度变化规律, 采用分子动力学方法对该问题进行了模拟研究 (通过研究分析模型形状对强度变化影响可以忽略) 。
表2为在计算模型截面半径r=14a、长度l=40a、温度300K, 加载方向采用[100]晶向下, 不同应变率加载相对应的钨合金屈服强度数据。
通过对表2的数据进行平均差值及数据拟合得到了钨合金拉伸的应变率及屈服强度的幂律拟合公式。如下:
图12为通过拟合的公式与钨合金的实验数据对比。对图12分析可知, 低应变率拉伸加载区域拟合的较好, 且相同的应变率条件下, MD模拟对应的屈服强度高于宏观实验的屈服强度, 分析认为, 钨合金是完美晶体没有任何缺陷, 而宏观实验材料内部很可能存在微小的空洞、微裂纹以及杂质等外在的因素从而降低了钨的屈服强度。MD模拟钨合金得到的屈服强度很高, 这一结论与Михаловский等[15]实验测得拉伸钨晶须的屈服强度24.75GPa近似, 这说明采取的钨合金计算模型、弛豫时间、加载方式是合理的, 拟合得到的幂律公式是准确的, 尤其对于在实验中很难对高应变率拉伸加载测试手段获得相应的塑性力学参数, 分子动力学方法显得具有更加重要实际意义。
5 结论
(1) 钨合金材料在相同晶向、截面尺寸、温度的条件下钨合金随着拉伸应变率增加, 屈服强度升高;且高应变率下的钨合金在屈服时对应的屈服强度更大, 并且容易发生二次屈服。
(2) 相同晶向、截面尺寸、应变率的条件下拉伸钨合金随着温度升高, 屈服强度不断降低;而相同晶向、温度、应变率的条件下拉伸钨合金截面尺寸增加, 屈服强度有所增加。
(3) 在相同温度、应变率、截面尺寸的条件下, 晶向效应的改变对钨合金材料屈服强度和伸长率具有较大的影响。
热力学模拟 篇2
随着国民经济的的不断进步和发展, 风机的产生在国民经济的生产发展中起到很大的促进作用, 风机将随着时代的发展, 不断更新技术研究, 从而能够更好的适应经济发展的需要, 传统的风机设计, 人们仅靠试验取得数据和经验公式, 试验发现问题, 改进设计。但由于试验研究方法受到各种条件的限制, 很多模拟参数的测量受到很多不良因素的影响, 给测量结果带来很大的困难, 很容易降低风机数值的实用性, 对风机数值测量的误差加大。而现阶段, 由于科学技术的不断发展, 利用商业CFD软件对风机的全三维流场进行模拟已越来越普遍, 也就是利用计算流体力学对风机进行数值模拟的研究, 给数值模拟工作带来了很大的便利, 通过对计算结果进行了分析, 模拟结果有助于理解风机内部的流动规律。
1 计算流体力学的概念分析
计算流体力学 (Computational Fluid Dynamics, 简称CFD) 起源于20世纪60年代, 当时的学科兴起跟计算机的技术发展有很大关系, 随着人们对其不断的发展和研究, 计算流体力学已经被广泛的应用, 各种商品化的CFD通用性软件开始应用这类力学研究, 同时更是对很多工业领域的生产发展起到很大的作用, 计算流体力学以计算机为基础, 利用数值的方法进行对流体力学各类问题的研究和模拟, 主要在离散格式、湍流模型与网格生成等方面进行相对的数值试验、计算机模拟和分析研究, 利用计算流体力学研发出得CFD技术, 不仅极大的克服了传统流体力学中不完善的问题, 而且还在应用领域得以全面的扩大, 很多核能、化工、建筑等领域都有其力学的涉略。风机在以上领域也有其所用之处, 为此, 计算流体力学对风机的设计和研究也有很大的作用。
2 风机的数值模拟分析
众所周知, 风机的国民经济发展的重要工具, 其在对生产过程中发出的大量湿、热、工业粉尘、甚至有害气体和蒸汽都有着有效的防护和净化处理的作用, 同时还能回收再利用, 有效的对资源进行合理的分配整合, 其中风机在纺织业的作用较为突出, 络筒机的离心风机提供了吸纱的作用, 不仅可以免去资源浪费, 还能减少纺纱机的能源消耗, 有效的提高纺纱质量, 具有更多的促进作用。在工业发展中, 风机从节能、降低噪声污染的角度来说, 尤其更大的促进作用, 因此在风机的设计原理上, 更多的要注重高效率, 但就目前市面上的风机产品, 可谓参差不齐, 很多规格和品种配套性极差, 为此在工业应用上也受到了很大的影响, 需要对已有的风机进行改造, 数字模拟其实是以电子计算机为工具, 把数学模型蕴藏的定量关系展示出来, 利用计算流体力学对风机的复杂流动问题的模拟计算, 通过数值离散求解流体运动方程, 揭示风机流体机理和流动规律, 从而研制出新的风机设计, 使整个产品从开发到运用都能够达到更为经济和省时的作用。
3 基于计算流体力学的风机数值模拟的应用
利用计算流体力学来研究风机的数值模拟, 这种方法对风机的设计提供更为依据原理, 对风机的不断完善起到促进作用, 其应用范围很广, 例如:通过对地铁专用轴流风机的设计来说, 这类风机主要应用在地铁车站和隧道区间内, 因其受都流量大、压头高和功率大等特点的制约, 试验成为了地铁轴流风机的设计检验的一般途径, 但是却在人力物力上有极大的消耗, 造成设计成本的浪费。为了克服这一弊端, 采用计算流体力学的原理, 对地铁轴流风机采用进行数值模拟, 主要是对地铁轴流风机在不同转速和安装角度进行模拟, 通过得出的最后结果进行指导设计方案, 并将模拟结果与厂家的试验数据作了对比, 酌情查处风机是否有需要改动之处, 从而提高风机的设计效率, 具有明显的应用价值和经济效益。
4 结论
以上对计算流体力学的风机数值模拟的分析和研究, 计算流体力学不仅是对风机的设计有很大的促进作用, 更大的提高风机的设计效率, 随着科学技术的进步, 其作用会越来越大, 充分了的利用计算机和数值数学的结合, 对流体力学的各类问题进行数值试验、计算机模拟和分析研究, 以解决实际问题。从而有助于人们对风机的构造设计进行深入了解和不断完善, 依靠合理的计算来优化风机的设计技术, 计算流体力学不仅是科学技术革新的依据, 更是极大满足了国民经济发展的需要, 计算流体力学进行对风机数值模拟的技术研究, 更是对设计高效率的风机具有重大意义。
摘要：风机是在国民经济发展的各个部门都被广泛使用的机器, 通常在冶金、石油、化工、纺织、电力、轻工等工矿业较为广泛, 在这些部门生产发展中起到很大的帮助作用, 随着计算机软硬件的发展水平的提高, 应用计算流体力学软件可以对风机进行数值模拟和分析, 为深入了解和分析风机数值带来巨大便利, 利用数值方法通过计算机求解描述流体运动的控制方程, 揭示流体运动的物理规律, 进而得出风机的工作原理, 本文同过对计算流体力学进行分析, 解析出风机的数值模拟, 不断完善发展风机技术。
关键词：计算流体力学,风机,数值模拟,发展前景
参考文献
[1]黄其柏.离心风机旋转频率噪声的理论与声辐射特性研究[D].西部大开发科教先行与可持续发展——中国科协年学术年会文集, 2009.
[2]姚宏, 王大枚, 雷丛林.浅圆仓五种机械通风方式比较试验[D].中国粮油学会第二届学术年会论文选集 (综合卷) , 2010.
热力学模拟 篇3
关键词：桥梁；剪力键；滑板式橡胶支座；地震响应；减隔震设计
最新的《城市桥梁抗震设计规范》1（后简称《规范》）中明确指出桥梁抗震体系包括两种类型：一种是通过桥墩部位形成塑性铰耗能的延性抗震设计2，另一种是通过采用柔性支座延长结构周期以避开地震能量集中区域，采用阻尼器等耗能装置来耗散地震能量的减隔震设计3.两者的主要差异在于所选择的耗能部件不同，前者是通过主体结构塑性变形耗能，容易使结构产生较大损伤且震后难以维修；而后者则是将易于更换的支座作为耗能部件，这显然更加切合实际.减隔震设计在国外使用较广，但国内多用于大型桥梁，并且减隔震装置多采用铅芯橡胶支座、高阻尼橡胶支座等4.对于高烈度区的中小跨径桥梁，采用上述减隔震支座往往会导致费用过高，因此选择一种既经济又合理可靠的减隔震措施尤为必要.
从2008年汶川地震桥梁震害5中可发现，一些因支座破坏导致主梁滑移错位的桥梁，其下部结构大多未出现明显损伤，这是由于支座滑动后能有效地减小主梁传递至下部结构的荷载，从而保护下部结构不受损伤，这同减隔震设计的思路是一致的.但支座破坏时并非一定出现滑移现象，还可能发生支座被“锁死”，此时就犹如固定支座一样，这将显著增大桥墩内力.由此可见，支座的滑移的确能起到保护下部结构的作用.但利用支座破坏后滑移来减隔震并不一定可靠，设置滑板支座的桥梁通过主梁滑动来减小主梁传递至下部结构的荷载，在保证桥梁抗震能力的同时可减小桥墩的尺寸.但由于滑板支座的低摩擦力，易导致主梁在常遇地震下发生纵向滑移.为改善该问题，本文提出一种外设剪力键的滑板橡胶支座系统.
在该支座系统中，剪力键作为一种“熔断保护装置”，在罕遇地震下能发生破坏，隔断上部结构传递至下部结构的惯性力，减小桥墩损伤，同时因剪力键破坏，支座产生滑移，充分发挥滑板橡胶支座摩擦耗能的能力，从而起到减隔震作用.而在常遇地震下，剪力键保持弹性状态，这既提高了支座整体刚度，也提供了足够的恢复力，防止主梁发生纵向滑移.Nielson和Choi等 6-7介绍了剪力键在OpenSees中的模拟方法；聂利英等8模拟了剪力键的力学性能，并对剪力键在异型板桥中的抗震性能进行了研究；夏修身等9对剪力键工作及破坏阶段分别进行了分析模拟，认为剪力键对位移峰值影响较小，但对墩底内力影响较大.孟兮10、陈浩11等人对剪力键进行了试验模拟和分析.以上这些研究从不同角度对剪力键抗震性能进行了分析，但均偏重于定性分析，缺少剪力键的定量分析.剪力键设置的不合理可能导致结构产生更大损伤.如何确定剪力键数量、尽量减小桥墩结构在地震中的损伤，同时体现基于性能的抗震设计思想，对改善结构的抗震性能尤为重要.
基于此，本文首先阐述了带剪力键的支座系统分阶段工作原理，根据基于性能的抗震设计思想并结合我国《规范》要求，明确地提出了剪力键的设计计算方法，采用数值分析方法给出了简易计算公式，以及支座系统在OpenSees中模拟方法.最后，结合一个具体算例，通过非线性时程分析，分析剪力键的减隔震效果，对比了不同数目的剪力键设置对减隔震效果的影响，验证了本文提出的剪力键设计方法的合理性.
1支座系统工作原理
本文所述支座系统主要包括两个部分：钢销剪力键和聚四氟乙烯滑板式橡胶支座，如图1所示.剪力键其一端嵌固于盖梁中，另一端伸入梁体但与梁体存在一定的间隙.上部结构自重产生的竖向力由聚四氟乙烯滑板式橡胶支座承担，而地震作用下产生的惯性力由支座和剪力键共同承担.在实际工程中，剪力键可采用盖梁处竖向钢筋外伸至梁底予以实现；间隙则通过主梁浇筑前，在钢筋周围包裹所需厚度的泡沫层予以实现.
带剪力键的支座系统各阶段的工作原理如下：
1）正常使用阶段：梁底圆孔和剪力键之间存在间隙，剪力键在本阶段不参与工作，因而不会影响主梁在温度、收缩和徐变作用下的自由变形.
2）常遇地震阶段：常温地震发生的可能性较大但震级较小，桥梁要求处于正常使用状态.剪力键在该阶段保持弹性，提高支座系统的水平承载力，保证主梁不发生纵向滑移.
3）罕遇地震阶段：罕遇地震发生可能性较小但震级一般较大，此时需要保证结构存在明确地耗能部位.该阶段剪力键将破坏从而保证滑板橡胶支座能发生较大滑移，充分耗散地震能量，同时减小上部结构传递至下部结构的水平力，避免下部结构出现过大损伤.
以上明确了支座系统需要两个参数：间隙gap和剪力键数量.gap大小的需要满足主梁自由伸缩的要求.剪力键数量则需要确保其在地震作用下分阶段工作；当剪力键数目较多时，桥梁响应接近固定支座桥梁；当剪力键数目较少时，桥梁响应接近普通滑板式橡胶支座桥梁.
2剪力键设计
我国《规范》采用两级设防水准的抗震设计思想，要求常规桥梁在E1作用下不受损坏或者不需修复即可继续使用，E2作用下不致倒塌或产生严重结构损伤并且经临时加固后可维持应急交通使用.同国外一些抗震设计规范12中的多性能水准相比较，我国《规范》的性能水准隐含在设防目标中，难以灵活选择，不太适合减隔震桥梁设计.因此，本文以《规范》为基础，结合基于性能的抗震设计思想，将桥梁分为3种设计类型，如表1所示，间接地分成多个性能水准.
设计类型
E1设防目标
E2设防目标
完全减隔震
主梁处于弹性范围
桥墩无损伤
主梁处于弹性范围
桥墩无损伤
部分减隔震
主梁处于弹性范围
桥墩无损伤
主梁处于弹性范围
桥墩处于部分损伤，出现塑性铰
非减隔震
主梁处于弹性范围
桥墩无损伤
主梁处于弹性范围，桥墩处于
严重损伤，保护层混凝土发生剥落
根据本文前述的3个阶段工作原理可知，剪力键gap参数可根据主梁在温度及收缩徐变作用下的位移值来确定.剪力键数目可通过式（1）计算：
3.1剪力键恢复力模型确定
3.1.1 数值模拟方法
采用ANSYS软件来精确模拟剪力键的力学特性.剪力键模型为圆柱体，底部完全固结，顶部按自由端来模拟；剪力键计算长度取为梁底至盖梁顶部的距离；荷载通过表面效应单元均匀施加在剪力键顶部；钢材本构关系采用双折线模型.分析剪力键的非线性受力全过程，得到剪力键的荷载位移关系曲线.
采用双折线简化得到的剪力键荷载位移曲线，如图2所示，其中A点为钢筋外层首次屈服点，ABC段及CD段面积相等.
3.1.2 简易计算公式
采用数值分析计算繁琐，不适用于设计.为了方便工程应用，根据数值分析结果，提出满足工程精度要求的剪力键恢复力模型的简易计算公式.
在剪应力作用下，剪力键截面并非严格意义的平截面，但为简化计算，在计算弯曲变形时仍采用平截面假定，这在剪应力不大的情况下能满足精度要求.对于圆形截面悬臂梁，弯曲正应力为三角形分布，剪应力呈二次抛物线分布，其计算公式为式（2）：
8计算得到极限位移，再考虑双折线本构关系的影响进行折减.对式9进行积分，得到极限位移的计算公式如10所示：
Δu=λ2∫Lφxdx=0.946λ2φ*yL2.10
式中：λ2取0.451 8；φ*y为采用极限强度计算得到的屈服曲率.
以长度100 mm，直径32 mm的剪力键为例，采用ANSYS分析得到的屈服、极限荷载分别是16.6 kN和23.5 kN；屈服、极限位移分别是0.72 mm和11.48 mm；而采用简易公式计算得到的屈服、极限荷载分别是16.8 kN和22.3 kN；屈服、极限位移分别为0.77 mm和10.9 mm.两者最大误差为6.5%，满足工程要求.
3.2支座系统模拟
如图3所示，支座系统可分为3个部分：剪力键、滑板橡胶支座和间隙.用OpenSees分析时，支座系统可采用零长度单元进行模拟，恢复力模型为特殊的多折线弹塑性模型，结合OpneSees材料库中几种普通材料以串、并联的方式进行模拟.
1剪力键为双折线模型，如图3a所示，采用Hysteretic材料进行模拟，材料参数包括初始刚度、屈服强度和屈服后刚度，均可通过上述的简易公式或通过对剪力键进行数值分析得到；
2间隙采用ElasticPPGap材料模拟，材料参数包括gap大小、初始刚度及屈服荷载.gap大小可通过计算结构的温度及收缩徐变下主梁变形值来确定；初始刚度和屈服荷载均取很大值例如刚度取剪力键刚度的100 000倍，且需保证gap材料的屈服位移大于剪力键极限位移，这是为了让gap模型仅提供一个间隙而不影响剪力键的刚度和承载力.由于剪力键的间隙在纵向是双向的，因此还需对两个方向的模型进行并联，最终得到图3b的双向gap模型；
3滑板橡胶支座为理想弹塑性模型，如图3c所示，采用Steel01材料进行模拟，材料参数包括初始刚度、屈服强度和屈服后刚度.为保证结构计算收敛，屈服后刚度不应取为零，而应设置为一个较小值.
最后将剪力键模型同gap模型串联，得到带间隙的剪力键模型；然后将滑板式橡胶支座模型与之并联，得到最终的支座系统模型.
为便于计算，首先视剪力键材料为理想弹塑性材料，其屈服强度取实际材料的极限强度，按式
4算例
4.1算例简介
算例为一座三跨连续梁桥，跨径布置为20 m+30 m+20 m，如图4所示.主梁为箱型截面、C50混
图3支座模拟示意图
Fig.3Schematic of bearing
凝土；桥墩采用2 m×1.2 m矩形截面、C30混凝土.
地震设防烈度为8度，Ⅱ类场地.桥墩处布置带剪力键的支座系统.为验证本文分析理论和简化分析，桥台处支座简化为理想滑动支座.
4.2模型建立
主梁在地震中一般处于弹性状态，因此本文中上部结构采用弹性梁单元模拟.桥墩在地震作用下可能进入塑性阶段，采用基于柔度法的弹塑性纤维梁单元模拟；其中无约束混凝土及约束混凝土均采用基于KentPark单轴混凝土模型模拟，纵向钢筋采用基于Menegotto和Pinto建议的模型.
4.3地震波选取
以《规范》提供的设计反应谱为标准，从美国太平洋地震工程研究中心的地震库中拟合6条地震波，其中3条对应E1作用的地震波，3条对应E2作用的地震波见表2.
5结果分析及验证
5.1分析工况确定
首先，计算得到式1中的参数.其中结构在E1和E2地震作用下墩顶水平力FE1和FE2分别为618 kN和1 869 kN.墩底对应无损伤的墩顶水平力Fex为764.7 kN，剪力键直径取32 mm，其极限承载力Fu为23.5 kN，滑板式橡胶支座所提供的最大水平力Ff为495 kN.计算合理的剪力键数目在6～11根.本文选取8根.另根据模型计算得到的温度变形值确定gap参数取为5 mm.
为对比设置剪力键数目对结构抗震性能的影响，分别计算设置过少剪力键，设置合理剪力键和设置过多剪力键3种工况，另外计算选取无剪力键滑板橡胶支座及固定支座两种情况作为参考，工况设置如表3所示.
5.2结果分析
1E1地震作用下主梁纵向位移
在E1地震作用下结构需处于弹性阶段，主梁纵向不应出现残余位移.由于墩柱在E1地震作用下均要求按弹性进行设计，结构内力大多均能满足要求.因此以下仅对比E1地震作用下结构的位移.
图5给出了wave2地震波下主梁的纵向位移时程图.从图中可看出，剪力键设置合理时，主梁位移始终在平衡位置两侧波动，不发生滑移.剪力键数量不足时，主梁位移将偏离平衡位置，且震后无法恢复，说明主梁发生滑移.以上表明若剪力键设置合理，能够有效地防止主梁发生纵向滑移，避免主梁出现震后残余位移.
图6给出了不同剪力键数目下主梁纵向位移峰值和残余位移变化趋势.从峰值位移可看出，随剪力键数目增加，主梁位移逐步减小.由残余位移可知，工况LS8和工况LS12的残余位移全部处于gap范围以内，同设置固定支座的情况接近，而LS4和LSH均产生较大的残余位移，这说明设置8根以上剪力键能保证其不发生破坏并始终处于弹性阶段，避免主梁发生滑移，而设置过少剪力键因其提前破坏导致主梁发生滑移.
以上分析表明，设置8根剪力键能满足E1地震作用下的设计要求.
2E2地震作用下墩底纵向水平力
对于完全减隔震的设计类型，E2地震作用下需保证桥墩等主体结构不发生损伤，且位移在允许范围内.采用带剪力键支座系统同其他减隔震桥梁一样会产生较大位移，可通过设置阻尼器、较宽的伸缩缝或合适的主梁搭接长度来控制或适应，以下对比E2地震作用下不同工况墩顶的水平力.
图7所示为一条地震波下墩底纵向水平力时程.由图可知，工况LSG墩底水平力可达1 400 kN，远高于墩底发生损伤所对应的水平力，这说明如果设置固定支座将导致桥墩发生较大程度损伤；其他工况在中后期的墩底水平力均被削弱至495 kN左右，说明在地震过程中剪力键均发生破坏，剪力键的破坏削弱了墩顶水平力和控制桥墩的内力.
由图8的墩底峰值内力可看出，随剪力键数目增加，墩底纵向水平力呈逐步增加的趋势，说明剪力键数目的增加会增大桥墩内力，桥墩内力可通过剪力键数目来控制.工况LS12峰值内力均超过765 kN，桥墩已进入塑性并产生部分损伤，说明剪力键设置过多，会导致桥墩内力过大而出现非预期的损伤.而设置8根剪力键均能保证桥墩处于无损伤状态.综合E1，E2两个阶段的结构响应分析可知，设置8根剪力键能满足两个阶段地震作用的设计要求，剪力键发挥了预期效果.
3支座系统滞回模型
绘制E2地震wave1作用下的剪力键滞回曲线如图9所示.从图中可看出，对于工况LSH，支座系统在水平力达到495 kN左右时开始发生滑移，此时恢复力为理想弹塑性，屈服力完全取决于摩擦力大小；对于工况LS4～LS12，支座系统屈服力随剪力键数目增加而增大，当剪力键发生破坏以后，支座系统发生滑移，后期同工况LSH一致.而工况LSG支座系统始终保持弹性不发生滑移.
前文中E2地震作用下的分析是在假定支撑长度足够且主梁不发生落梁破坏的基础上，因此有必要对这种假定的合理性进行讨论.从图9中可知，剪力键发生破坏后支座会产生较大位移，在不引入外界耗能装置的前提下，这是减隔震支座共同存在的一个缺点.然而，工况LS8在wave1下墩梁相对位移峰值为16 cm，小于《规范》抗震措施规定中对最小支承长度的要求，这表明只要按规范设置合理的支撑长度，本文算例桥梁在E2地震下不会发生落梁破坏，即前文分析的假定是合理的.
6结论
本文以高烈度区某等跨径桥梁为例，分析了带剪力键支座系统对桥梁抗震性能的影响，得到如下结论：
剪力键数目对支座系统的抗震性能有显著影响.设置过少剪力键会使其在常遇地震下过早发生破坏，无法充分发挥作用；设置过多剪力键会导致其在罕遇地震下难以发生破坏，无法起到控制下部结构内力的作用，无法充分发挥支座耗能能力.
本文所述支座系统造价低廉，施工方便，通过合理设置剪力键的滑板式橡胶支座能够弥补普通滑板支座的缺点，满足正常使用阶段、常遇地震阶段及罕遇地震阶段不同的需求.采用本文提供的剪力键恢复力模型及剪力键数量的计算方法能够准确合理地设置剪力键，达到减隔震效果.
剪力键的设置和桥墩损伤有关，为使整个结构达到合理的抗震性能，需将剪力键设计同墩柱设计结合起来.单独设置剪力键并不一定能提高结构整体的抗震性能.
参考文献
1]CJJ 166—2011 城市桥梁抗震设计规范S].北京：中国建筑工业出版社，2011：7-10.
CJJ 166—2011 Code for seismic design of urban bridgeS].Beijing： China Architecture & Building Press， 2011：7-10. In Chinese
2]范立础.桥梁延性抗震设计M].北京：人民交通出版社， 2001： 74-95.
FAN Lichu. Ductility seismic design of bridgesM]. Beijing： China Communication Press， 2001： 74-95. In Chinese
3]范立础，王志强.桥梁减隔震设计M].北京：人民交通出版社， 2001： 140-144.
FAN Lichu， WANG Zhiqiang. Seismic isolation design of bridgesM]. Beijing： China Communication Press， 2001： 140-144.In Chinese
4]庄军生.桥梁减震、隔震支座和装置M] .北京：中国铁道出版社，2012：45-75.
ZHUANG Junsheng. Bridge isolation and deviceM]. Beijing： China Railway Publishing House，2012：45-75.In Chinese
5]王克海，李冲，李茜，等. 考虑支座摩擦滑移的中小跨径桥梁抗震设计方法J]. 工程力学， 2014，316： 85-92.
WANG Kehai， LI Chong， LI Qian， et al. Seismic design method of small and medium spans bridge considering bearing friction slipping J]. Engineering Mechanics，2014，316：85-92.In Chinese
6]NIELSON B G. Analytical fragility curves for highway bridges in moderate seismic zonesD].Atlanta， Georgia： Georgia Institute of Technology， 2005：295-298.
7]CHOI E. Seismic analysis and retrofit of midAmerica bridgesD].Atlanta， Georgia： Georgia Institute of Technology， Civil and Environment Engineering，2002：47-66.
8]聂利英，李建中，胡世德，等. 抗震销的模拟及其对异型板桥抗震性能影响分析J]. 地震工程与工程振动，2008，282：108-113.
NIE Liying，LI Jianzhong，HU Shide，et al.The effects of antiseismic bolt on seismic estimation of heteromorphic bridge and its simulationJ]. Journal of Earthquake Engineering and Engineering Vibration，2008，282：108-113.In Chinese
9]夏修身，陈兴冲，王希慧，等. 剪力键对隔震桥梁地震反应的影响J]. 地震工程与工程振动，2012，326：104-109.
XIA Xiushen， CHEN Xingchong， WANG Xihui，et al.Effect of shear key on seisrnic responses of bridge using isolation bearingJ]. Journal of Earthquake Engineering and Engineering Vibration， 2012，326：104-109. In Chinese
10]孟兮，倪燕平. 减震榫设计及试验研究J]. 北京交通大学学报，2013，373：103-106.
MENG Xi，NI Yanping.Design and experimental study of shock absorberJ]. Journal of Beijing Jiaotong University， 2013，373：103-106. In Chinese
11]陈浩. 减震榫力学性能及其在桥梁减震中的应用D].北京：北京交通大学，2010：47-69.
CHEN Hao. Mechanical characteristic of absorber and its application on the bridgesD].Beijing：Beijing Jiaotong University， 2010： 47-69. In Chinese
热力学模拟 篇4
冲击作用下多孔材料热力学特征的模拟与分析
强冲击波的作用可在多孔材料中诱发复杂的时空耗散过程,在这期间系统整体处于远离平衡的`状态,对这一过程的稳定模拟和结果分析均具有较强的挑战性.使用近期针对超高速碰撞而发展起来的物质点方法对这一过程进行模拟,引入积分几何和数值图像处理中的形态学描述来处理和分析系统热力学特征量例如温度的Turing斑图,揭示出3个Minkowski泛函(白色区域相对面积、边界总长度、欧拉特征量)与系统中“高温区”所占份额、“热点”在空间的分布方式之间的对应,揭示出Minkowski泛函演化特征与冲击波及多孔材料相互作用过程之间的对应.研究了孔隙度和冲击波强度对物质状态参量的影响.
作 者：许爱国 张广财 蔚喜军 潘小飞 朱建士 Xu Aiguo Zhang Guangcai Yu Xijun Pan Xiaofei Zhu Jianshi 作者单位：北京应用物理与计算数学研究所计算物理实验室,北京,100088 刊 名：中国工程科学 ISTIC英文刊名：ENGINEERING SCIENCE 年，卷(期)： 11(9) 分类号：O347.5 关键词：冲击波 多孔材料 物质点法 形态学量度热力学模拟 篇5
当今计算机技术发展迅速, 我国的工程分析软件发展处于落后状态, 并且交叉学科的应用技术也很浅薄。程序设计作为一种重要的研究手段往往在博士学习阶段才会被重视, 并被视为一种高难度的工作。这实际上是不应该的。
在机械中通常用有限元法进行结构力学分析, 并能做到形变的模拟。此方法是将一个连续的弹性体离散化, 将其看作由弹簧单元组成的网格体, 并根据载荷—位移关系联立方程组求解。其步骤为网格划分、生成方程组、迭代计算, 这些数学过程的实现都是繁琐的。节点的数据存储、方程组的描述占用较大的存储空间, 按部就班地开发求解器也耗时耗力。所以, 针对平面弹性力学问题, 可以对有限元法的分析过程进行重新设计, 其中最主要的是利用差分法将问题大幅简化, 使其能够轻易地设计迭代程序来实现问题的求解和数值模拟。
一个连续体就可以利用最简单的差分法来进行离散。差分公式的推导多种多样, 且大部分是有针对性的, 并不通用。为了得到对大多数问题具有通用性、适合程序设计、数据结构描述方便的差分方程数学模型, 这里引用“节点领域”的平衡条件推导差分方程的方法。
1 位移平衡方程
为了方便数据对形体的描述, 将一个均匀材质的平面弹性体假设为无数个微小正方形单元组成。每个节点作为待定节点, 与周围的节点存在一定的位移平衡关系, 这也是一个差分方程所描述的意义。根据离散的情况, 节点位移平衡方程有4种类型, 其中有3种边界上的类型 (图1) 。
下面列出通过“节点领域”的平衡条件推导差分方程所得到的结果。原方程变换整理为矩阵形式。
(1) 非边界单元节点结构见图2。
每个节点有u、v, 水平、垂直两个位移分量:
undefined
其中, E为弹性模量, ν为泊松比, f为体力。
(2) 一般边界如图3所示, 有4个方向, x、y轴方向与图2相同。
其中, (a) (b) 所示方向的位移平衡方程为
undefined
其中, t为载荷。
由于对称性, 它们的差分方程相同, 但节点位置不同, 下面也同样。 (c) (d) 为
undefined
(3) 角边界见图4。
其中, (a) 、 (b) 所示方向的位移平衡方程为
(4) 内尖角边界见图5。
其中, (a) 、 (b) 所示方向的位移平衡方程为
以上就是4种情况下所需要的节点间关系的描述方程。
2 迭代计算
以梁弯曲问题为例, 见图6。有3种平衡方程类型, 其中②和③有4个方向。
首先, 描述节点的数据, 以C++为例, 代码如下:
由于网格单元的规律性, 使用一个二维数组直接描述形体:
对方程组的描述利用节点的相互位置, 如图7所示。
这样, 在程序设计上, 就直接省去了一般有限元法的网络划分的繁琐过程。将上述位移平衡方程整理如下:
undefined
其中, A为系数矩阵, i为节点编号, 用于迭代计算图7中节点[i, j]的位移。整理好方程用程序描述。下面以如图2所示的非边界方程类型为例的代码描述:
另外, 两个边界类型的方程4个方向的描述在这里省略。此例子的程序中所描述的方程共有9个, 用哪种方程进行迭代计算需加以判断。
如果某节点设为约束, 则对应的方程不参与迭代计算, 以保证约束的位移分量为0。
于是每个节点根据临近节点关系, 即平衡方程, 进行“扫描式”的迭代计算。所有节点按照顺序计算一遍完成一次迭代, 并反复进行。在没有对方程组完整描述的情况下实现Gauss-Seidel迭代法。
另外, 对于约束节点, 可将方程整理为的形式来计算约束反力。
图8为300单位长度的梁, 左侧边界5个节点设为固定约束, 右侧边界5节点设为方向向下、值为2的均匀载荷的迭代结果。
3 结束语
将任意弹性离散为只由正方形单元所组成, 简化了节点坐标数据的存储, 和一般有限元法中的网络划分的过程, 也为结果的可视化和后处理提供了方便。
对迭代求解所需要的方程组, 其描述上不需要耗费存储空间去完整地建立, 同时迭代程序简洁和直观的设计, 也降低了问题的抽象程度, 从而快速实现弹性问题的模拟计算。
通过“节点领域”的平衡条件推导差分方程的方法具有一定通用性, 特别是边界的差分方程上, 可以适用于各种弹性问题, 也为迭代程序的设计提供了便捷。
摘要：为了能够简捷、快速地得到一般平面弹性力学问题的结果, 将一个平面弹性体假想为由若干正方形组成的近似结构, 这个结构可以简单地用二维数组直接进行描述。以梁弯曲为例, 引用“节点领域”的差分法推导每个节点与周围节点的平衡关系, 利用每个节点的关系方程进行迭代计算, 从而编写程序实现了模拟梁弯曲变形的过程。总结了适合于程序编写、供各种基本边界条件下, 联立方程组迭代求解弹性体的位移分布的简单思路。
关键词：差分法,迭代法,程序设计,弹性力学
参考文献
[1]徐芝纶.弹性力学中的差分方法[M].北京:高等教育出版社, 1989.
[2]张传力, 冯定, 彭兴黔.弹性力学平面问题的位移函数及其差分格式[J].江汉石油学院学报, 1997 (3) .
热力学模拟 篇6
在对蛋白质折叠的分子模拟中, 一般认为不易形成折叠结构的主要原因与力场的选择或模拟时间有关, 为此本文将刚体动力学模型应用于研究力场因素对蛋白质多肽体系在形成螺旋结构过程中的影响, 由于丙氨酸饱和多肽链具有很高的螺旋性质, 并且在理解蛋白质结构方面起着非常重要的作用, 所以它被广泛的研究。本文选用的研究对象是17个丙氨酸的饱和多肽链.它的N端和C端分别由两个甲基 (CH3) 饱和。
1. 模型和方法
1.1 刚体动力学模型
本文把丙氨酸饱和多肽链划分为两种不同的多体 (刚体, 即原子集合) 系统, 这里以三个丙氨酸的饱和多肽链为例来说明这两种刚体结构, 如图1所示.图中所有的肽平面和侧链都被定义成为刚体, 刚体间由化学键链接.需要强调的是在结构I中, 丙氨酸饱和多肽链末端的两个甲基被单独定义成为刚体, 这是结构I和结构II间的唯一区别。
1.2 模拟方法
本文中, 笔者用TINKER软件包进行了相关的分子动力学模拟.为了研究丙氨酸饱和多肽链折叠分子动力学模拟对力场的依赖性, 本文采用了AMBER[11]和AMOEBAPRO[12]两种力场。
1.2.1 全原子模拟
17个丙氨酸的饱和多肽链的初始状态是由TINKER生成的直链结构.在隐性水条件下, 采用全原子模型进行分子动力学模拟, 分子动力学模拟在NVT系统下运行, 温度恒定为298K, 考虑到模拟体系的尺寸, 体系的预平衡时间为50ns, 而本文进行分子动力学模拟的时间为200ns, 在全原子模拟中, 时间步长取为1fs, 坐标文件每1ps保存一次。
刚体结构I中共有7个刚体, 分别用“body+1…7”表示, 刚体结构II中有5个刚体, 分别用“BODY+I…V”表示。
1.2.2 刚体动力学模拟
为了执行刚体动力学模拟, 在TINKER输入文件中, 笔者首先根据两种划分方式定义了不同的刚体结构, 17个丙氨酸的饱和多肽链的结构I由37个刚体组成, 而结构II则由35个刚体组成, 这两种刚体系统都能形成链状的拓扑结构, 为了与全原子模拟结果相比较, 刚体动力学模拟同样选定在NVT系综下运行, 温度恒定为298K, 模拟时间同样取200ns, 坐标文件每1ps保存一次, 唯一的不同是刚体动力学模拟所采用的时间步长为2fs。
2. 结果与讨论
研究分子动力学模拟效率的一个常用量是计算相对于参考结构的均方根偏差 (RMSD) 。本文仅计算了丙氨酸饱和多肽链中所有重原子 (即不考虑氢原子) 相对于参考结构随时间变化的均方根偏差, 17个丙氨酸的饱和多肽链的标准螺旋结构由SYBYL8.1生成, 本文对所有丙氨酸饱和多肽链的模拟轨迹都与参考结构进行了比较研究。
首先笔者在表1中给出了针对17个丙氨酸饱和多肽链的所有分子动力学模拟轨迹的最小均方根偏差的值.通过对分子动力学模拟轨迹均方根偏差的分析, 我们发现在结构I的刚体动力学模拟中, 在两种力场下都没有形成螺旋结构.而在结构II的刚体动力学模拟中, 采用AMBER力场形成了螺旋结构, 且在AMOEBAPRO力场中则没有形成的螺旋结构。
图2给出了基于全原子模型和刚体动力学模型在不同力场下得到的模拟结果。从图2 (a) 中, 我们发现在AMBER力场下, 基于结构II的刚体动力学模拟在50ns内就形成了稳定的折叠结构, 相对于17个丙氨酸的饱和多肽链的参考结构, 结构II的刚体动力学模拟得到了最小的均方根偏差值, 而结构I的刚体动力学模拟却得到了较大的均方根偏差值, 需要特别指出的是这种现象只出现在AMBER力场中, 对此笔者认为由于丙氨酸多肽链的末端是用两个甲基饱和的, 在刚体结构I中, 末端的甲基被单独定义为刚体, 所以系统末端的运动相对灵活。笔者认为这也就是在AMBER力场中结构I的刚体动力学模拟不能得到折叠状态的重要原因。
极化力场是可以更精确研究化学过程的新一代力场, 为了测试刚体动力学模型在极化力场中对丙氨酸饱和多肽链的折叠的分子模拟, 本文选用了AMOEBAPRO力场来研究。图2 (b) 给出了在AMOEBAPRO力场下, 无论用全原子模拟还是用刚体动力学模拟, 17个丙氨酸饱和多肽链在200ns内都不能形成螺旋结构。笔者认为这是由于AMOEBAPRO力场精细地计算分子体系内部原子间的相互作用, 若想在基于刚体动力学模型的模拟中观察到螺旋结构需要进行较长时间的分子模拟。
结论
本文基于刚体动力学模型选用了2种力场对17个丙氨酸饱和多肽链的折叠过程做了分子模拟, 并将结果同全原子模拟进行了比较。研究结果表明用刚体动力学模型对丙氨酸饱和多肽链折叠的分子动力学模拟依赖于对力场的选择, 在AMBER力场中, 刚体动力学模拟的结果要好于全原子模拟, 需要注意的是在不同的力场中, 不同刚体结构的定义也会对丙氨酸饱和多肽链的折叠产生影响。对刚体动力学模拟而言, AMBER力场是在两端的甲基没有被单独定义为刚体结构的情况下, 研究丙氨酸饱和多肽链折叠选用的最好力场。在AMORBAPRO力场下的刚体动力学模拟中没有观察到螺旋结构的主要原因是力场自身特点的不同及模拟时间不够, 最后需要指出的是在本文中笔者仅讨论力场及不同刚体结构的定义对研究蛋白质多肽体系折叠的影响, 但对较大蛋白体系的研究有待于以后开展大量的研究。
摘要：根据丙氨酸饱和多肽链的分子结构, 定义了两种不同的刚体结构, 然后基于刚体动力学模型采用AMBER和AMOE-BAPRO两种力场对丙氨酸饱和多肽链的折叠进行分子模拟, 并将结果同全原子模拟做了比较.研究结果表明利用刚体动力学模型研究丙氨酸饱和多肽链的折叠不仅与刚体结构的划分有关, 还依赖于对力场的选择。
关键词：刚体动力学模型,丙氨酸饱和多肽链,蛋白质折叠,分子力场,分子模拟
参考文献
[1]DOBSON C M.Protein folding and misfolding[J].Na-ture, 2010, 426:884-890
[2]DUAN Yong, KOLLMAN P A.Pathways to a proteinfolding intermediate observed in a 1-microsecond simulation inaqueous solution[J].Science, 1998, 282:740-744
[3]EBBINGHAUS S, DHAR A, MCDONALD J D, et al.Protein folding stability and dynamics imaged in a living cell[J].Nature Method, 2009, 7:319-325
[4]KARPLUS M, PETSKO G A.Molecular dynamics sim-ulations in biology[J].Nature, 1990, 347:631-639
[5]CLEMENTI C.Coarse-grained models of protein fold-ing:toy models or predictive tools?[J]Curr Opin Struct Biol, 2008, 18:10-15
[6]CHUN H M, PADILLA C E, CHIN D N, et al.MBO (N) D:A multibody method for long-time molecular dy-namics simulations[J].J Comput Chem, 1997, 21:159-184
[7]MUKHERJEEA R M, CROZIERB P S, PLIMPTONBS J, et al.Substructured molecular dynamics using multibody dy-namics algorithms[J].Inter J Nonlin Mech, 2008, 43:1040-1055
[8]PONDER J W, CASE D A.Force fields for protein sim-ulations[J].Adv Protein Chem, 2003, 66:27-85
[9]COUCH V A, CHENG N, NAMBIAR K, et al.Struc-tural Characterization ofα-Helices of Implicitly SolvatedPoly-alanine[J].J Phys Chem B, 2006, 110:3410-3419
[10]PONDER J W.TINKER molecular modeling package, Washington University Medical School
[11]WANG J, CIEPLAK P, KOLLMAN P A.How welldoes a restrained electrostatic potential (RESP) model perform incalcluating conformational energies of organic and biologicalmolecules?[J]J Comput Chem, 2000, 21:1049-1074
热力学模拟 篇7
沥青老化对多孔沥青混凝土的抗松散掉粒性能有很大的影响。沥青老化后在低温下变得更硬和更脆,抗疲劳性能也急剧下降,在反复行车荷载作用下,多孔沥青路面表面集料更易松散脱落 [4]。多孔沥青路面在荷兰的实际使用情况亦表明松散掉粒更容易发生在冬季寒冷季节,而在夏季高温季节,表面掉粒并不明显,永久变形是主要病害。
为了研究沥青老化对多孔沥青混凝土性能的影响,该文通过室内长期老化试验模拟多孔沥青混凝土在现场的老化过程。为了能更好地模拟现场长期老化过程,所采用的室内长期老化模拟试验综合了紫外线、热氧和湿气的共同作用。沥青及其砂浆经室内老化后,利用动态剪切试验,建立动态剪切模量和相位角主曲线,分析老化对其力学响应的影响。
1 原材料与试验
所用沥青为壳牌沥青公司提供的SBS改性沥青(Cariphalte XS),其主要性能指标为:25 ℃时针入度为50～85(0.1 mm),软化点大于65 ℃,5 ℃延度大于35 cm。所采用的沥青砂浆由沥青、填料和细集料按重量比0.34∶0.30∶0.36在175 ℃充分拌和而成。其中,填料含75%的矿粉和25%消石灰;细集料仅采用小于0.5 mm的部分。
为了更好地模拟多孔沥青路面中的现场老化过程,沥青及砂浆先后进行了短期和长期老化试验。短期老化试验方法为:将沥青或沥青砂浆铺成2 mm厚,放在175 ℃的鼓风烘箱内1.5 h来模拟现场拌和、摊铺和碾压时的老化过程。经室内短期老化后,将所得样品放入加速老化试验箱内进行长期老化模拟,其试验条件是:温度40 ℃,湿度70%,紫外光(波长300～400 nm)强度60 W/m2,老化时间1 000 h。研究表明上述老化方法可模拟多孔沥青路面早期3年的现场老化效果[4]。
采用DSR动态剪切流变仪(DSR MCR 101)分析沥青的动态剪切复合模量与相位角。所用试验平台为平行板频率扫描模式,即在某一温度下测量不同频率的力学响应。试验按应变控制,所施加的应变在线性粘弹性范围内。试验温度为-10～60 ℃,试验频率为0.1～50 Hz。当温度低于25 ℃时采用直径为8 mm的平行板,沥青层厚度为2 mm;当温度高于25 ℃时则采用直径为25 mm的平行板,沥青层厚度为1 mm。利用时温等效原理,可将不同温度测得的剪切复合模量与相位角平移,得到二者的主曲线[6]。
DSR动态剪切流变仪的平行板一般只适用于沥青或由沥青和填料组成的胶浆。对于含有小于细集料的沥青砂浆体系,细集料的尺寸效应将会影响试验结果。一般地,为了减少集料尺寸效应,试件的尺寸宜比最大集料粒径大一个数量级。为此,针对沥青砂浆的特点,在DSR动态剪切流变仪的基础上,专门开发了成套试验平台,如图1所示。砂浆柱的直径为6 mm,试样高20 mm。试件两端分别放置一个钢环,以便于安放试件。根据DSR动态剪切流变仪内部加载系统,设计了与试件相匹配夹具。试件的成型过程如图2所示。试模由硅胶材料按试件尺寸加工而成。试件一次可浇铸5个,成型时沥青砂浆的温度控制在175 ℃。试模先冷却至室温,再放入5 ℃温度下4 h后脱模。
通过上述对试件和夹具的专门设计,可利用DSR动态剪切流变仪进行沥青砂浆剪切力学试验。沥青砂浆的动态剪切试验条件与方法和前述的沥青试验相似,但由于沥青砂浆柱温度高时会严重软化,试验最高温度为35 ℃。试验采用小应变控制模式,以保证沥青砂浆的应变水平在其线性粘弹性范围内。
对于所设计的沥青砂浆试验,沥青砂浆柱所受到的应力和应变可按下式计算
undefined
式中:τ为最大剪切应力,MPa;T为加载扭矩,N·mm;R为沥青砂浆柱半径,3 mm;γ为最大剪切应变;θ为位移角,rad;heff为沥青砂浆柱有效高度,mm。为了精确确定heff值,在有限元模拟平台对沥青砂浆试件经有限元模拟计算,其值取12.742 mm [5]。
2 结果讨论与分析
图3为未老化、短期老化和长期老化沥青的剪切复合模量与相位角主曲线。主曲线建立的参考温度为10 ℃。如图3所示,剪切复合模量随着老化程度的增加而增加,而相位角则随着老化程度增加而减小。这种趋势在低频区表现更明显。如放大图所示,在高频区,即使放大后也发现老化对复合模量和相位角的影响很有限。未老化沥青和短期老化沥青的主曲线明显呈现出一平台区。这个平台区表明了SBS改性剂在沥青内形成了网络结构。然而经过1 000 h的紫外老化、热氧老化等综合作用后,这一平台区已逐渐消失,表明了SBS网络结构在长期老化模拟试验过程破坏严重 [7,8,9]。
图4为新拌沥青砂浆、短期老化和长期老化沥青砂浆的剪切复合模量与相位角主曲线。主曲线建立的参考温度为10 ℃。为了方便比较,未老化沥青的主曲线也加入图4中。如图所示,由于细集料和填料的加入,沥青砂浆的复合模量明显高于纯沥青,但二者对于相位角的影响极为有限,其原因可能是在小应变下,沥青砂浆的变形主要由内部沥青起主要作用。老化也导致沥青砂浆复合模量增加,而相位角减小。图4中的数据显示短期老化的效果略强于长期老化,这种结果的产生是沥青砂浆的材料变异性引起的。由于细集料和填料的加入,使得砂浆试件的不均匀性比纯沥青样高。如图5所示,沥青砂浆试验重复性差的情况确实存在。两个短期老化沥青砂浆测得的试验数据相差较大,而两个长期老化样品的试验数据则很接近。上述结果表明了沥青砂浆在制备过程中较纯沥青易出现材料变异。由于材料变异的影响可能会掩盖了老化的影响,因此为了更好地评价老化对力学性能的影响,宜进行平行试验加以验证。
3 结 论
a.室内长期老化模拟试验综合了紫外老化、热氧老化和湿气的影响,可更真实地反映多孔沥青混凝土现场长期老化过程。
b.通过试件尺寸和试验夹具的特殊设计,可在DSR剪切流变仪基础上开发出适用于沥青砂浆的成套试验方法。
c.SBS聚合物网络结构使得主曲线中出现明显的平台,长期老化后平台逐渐消失,表明SBS网络结构在老化过程中受到破坏,可严重降低改性效果。
e.沥青砂浆在制备过程中易出现材料变异,宜进行平行试验加以验证。
参考文献
[1]王牧,刘长军,刘明.多孔排水性沥青混合料性能的研究[J].公路,1997(5):35-46.
[2]毛爱明.排水沥青路面压实度超密现象分析[J].公路,2006(7):158-160.
[3]郑鑫,雷学坤,章建龙,等.国内外低噪声沥青路面研究状况概述[J].公路与汽运,2007(3):67-69.
[4]Hagos ET.The Effect of Aging on Binder Properties of Porous Asphalt Concrete Aging[D].Delft University of Technol-ogy,the Netherlands,2008.
[5]Mo LT.Damage Development of the Adhesive Zone and Mortar of Porous Asphalt Concrete[D].Delft University ofTechnology,the Netherlands,2010.
[6]赵延庆,吴剑,文健.沥青混合料动态模量及其主曲线的确定与分析[J].公路,2006(8):163-166.
[7]应荣华,郑健龙,陈骁.沥青老化效应的试验研究[J].公路交通科技,2007(6):20-24.
[8]陈华鑫,周燕,王秉纲.SBS改性沥青老化后的动态力学性能[J].2009(1):1-5.
热力学模拟 篇8
堵国成,陈坚,DU Guo-cheng,CHEN Jian(工业生物技术教育部重点实验室,江苏,无锡,214036)
期 刊：水资源保护 ISTICPKU Journal：WATER RESOURCES PROTECTION 年，卷(期)：, 23(2) 分类号：X788 关键词：Fenton反应 甲基橙 动力学模型
热力学模拟 篇9
1整体—局部均质化方法
均质化法是20世纪70年代基于小参数渐近展开的多尺度摄动理论发展起来的一种方法, 它结合有限元技术能够解决有周期性分布复杂结构的刚度问题和微观应力分布问题, 也可用于材料细观结构的拓扑优化, 但很少用于研究建筑结构中的钢筋混凝土材料。本文利用Luo等[10]导出的有效力学性质公式计算钢筋按不同体积比分布于混凝土材料时, 钢筋混凝土材料的力学性质。
设宏观变量xi与微观变量yi之比为摄动参数η:
包含整体与局部坐标的位移控制方程为:
ui (x) =u
利用链式微分规则可得到用特征函数表示的均质化弹性常数D
其中, 弹性模量Eijkl包含了特征体积单元上钢筋和混凝土的材料性质, 特征函数张量
2复合材料宏观弹性特性分析的微观力学方法
2.1 修正的蔡氏经验公式
复合材料弹性常数的计算有许多经验公式, 但都无法直接用于钢筋混凝土材料, 修正的蔡氏半经验半理论公式通过调整修正系数能够精确地计算各种复合材料的力学性质, 具体计算可参阅文献[6]。
2.2 经典Mori-Tanaka方法
为检验整体—局部均质化法的计算结果, 本文同时用平均Mori-Tanaka理论考察了钢筋含量对钢筋混凝土材料有效弹性模量的影响。
(C
设钢筋在复合材料中呈现横观各向同性性质, 因此, 式 (5) 中的C
C
C
其中, λ0, μ0, λ1, μ1分别为基体和夹杂物的拉梅常数;δij为Kronecker函数, 将式 (6) , 式 (7) 代入式 (5) 后得到二相纤维增强复合材料的弹性模量、剪切模量及泊松比。
3整体—局部均质化法在钢筋混凝土等效模量分析中的应用
利用整体—局部均质化理论计算了三维钢筋混凝土材料的整体力学性能, 主要计算结果如图1～图4所示。
由图1, 图2可知, 三种方法算得的EL相差很小, 工程法算得的ET相差比较大。由图3, 图4发现, 三种方法算得的剪切模量, 泊松比都很接近, 可以证明:用整体—局部均质化法计算钢筋混凝土材料的有效弹性模量是非常有效的。
采用两根尺寸相同的 (200 mm×450 mm×2 400 mm) 试件计算了两端简支梁受集中载荷和对称载荷作用时的应力和变形, 其中钢筋混凝土材料的两根预应力钢筋和两根非预应力钢筋分别在受拉侧和受压侧单排布置。各材料的力学参数如表1所示, 用整体—局部均质化法计算得到的等效材料的力学性质如表2所示。
如图5, 图6所示, 简支梁的挠度和应力均随荷载的增加呈线性增加, 同一受力情况下的原始材料和等效材料的结果非常接近。
4结语
1) 采用整体—局部均质化法计算钢筋混凝土材料的有效弹性模量非常有效, 可简化费时的建模过程, 为大体积、复杂的结构仿真分析提供一种新的简单高效的方法, 具有广泛的工程应用前景。
热力学模拟 篇10
1 MD模拟方法
1.1 模型的确立
碳纳米管的模型大体可以分为两类。第一类是刚性模型。第二类是柔性模型,如简谐弹簧势能和余弦势能函数来描述碳纳米管中碳碳原子间的各种作用,该模型的最大优点是可以精确的描述碳纳米管本身的性质。
水分子的势能模型种类有很多。如基于量子力学从头计算的MCYL模型[12],半经验型的ST2[13],SPC[14],TIP3P[15],和TIP4P[15]等。其中SPC (Simple Point Charge)模型形式简单,用来进行常温下的内能和径向分布函数等计算较为准备,故被广泛应用,本本文采用水-水分子间作用,其位能形式如下:
式中前半部分为静电作用,后半部分为短程Lennard-Jones(LJ)作用。ε为Lennard-Jones能量参数,σ为Lennard-Jones半径。u(ri,rj)为水分子间的相互作用势能,qi代表第i个原子或离子上所带有的电荷,rij为原子或离子i与j间的距离。SPC势能模型参数见表1。
1.2 MD模拟细节
本文模拟了矢量参数n,m分别为8,9,10的单臂碳纳米管,其直径根据下式计算:
其中a是C-C键长,为0.142 nm。具体参数见表2。
本文的模拟选取在NVT等温等压系踪中进行,采用Nose-Hoover控温法,控制温度保持在300 K,控制压强在1 atm。模拟过程中采用(8,8)(9,9)和(10,10)碳纳米管模型直径分别为1.085 nm,1.22 nm和1.356 nm,管长为73.79 nm,模型中分别放有36个,40个和50个甲醇分子。模拟总步数7.0×105,其中前2.0×105步为平衡步数,后4.0×105步用于统计各种性质,每隔100步保存一次轨迹文件,所采取的步长为0.5 fs。
2 结果与讨论
2.1 径向分布函数
径向分布函数g(r)是反映分子微观结构特征的物理量,可提供关于分子结构性质的信息,它的物理意义是指距离中心粒子r处出现另一粒子的概率密度相对于随机分布密度。通过对O-O和O-H径向分布函数的研究,便于了解不同状态下水分子的微观结构。我们通过分子动力学模拟得到不同管径下的水分子中的O-O, H-H 和O-H径向分布函数(见图1)。
从图1可以看出O-O和O-H径向分布函数分别在0.151 nm和0.299 nm处出现峰值。由此得出,随着碳纳米管管径的增大,径向分布函数的波峰位置影响较不明显,但是随着管径的增大,波峰峰值逐渐增加。其中能够代表甲醇分子微观结构特殊性的是氢键的径向分布函数,即若两水分子O-O间的距离小于0.35 nm、O-H之间的夹角小于30°,则认为两水分子间形成氢键。图1c所示水分子之间形成氢键的线性分布,随着管径的增大,波峰峰值逐渐增大,说明随着管径的增大,水分子之间的O-H聚集密度增大,氢键作用得到加强,而波谷随着管径的增大逐渐下降,且位置左移,峰宽逐渐变窄,说明体系结构由松散、无序逐渐变得紧密,有规律。
水分子放入(8,8)型碳纳米管,模拟结束时这些水分子形成了一条由氢键相连的S形纵列。(10,10)型碳纳米管内的水分子相对于(8,8)型碳纳米管内水分子,由于管内空间增大而导致受限效应减弱,分子的取向已完全成随机分布。延碳纳米管管轴方向观察水分子在管内的分布状态,水分子在靠近管壁的内层由氢键相连成环状.而在管轴附近和远离管壁的管间位置,基本上没有水分子。
2.2 水分子在碳纳米管中的自扩散系数的计算
图2显示了水分子在(8,8)(9,9)和(10,10)碳纳米管中扩散的MSD随时间的变化关系曲线。自扩散系数可通过两种等价的方法计算,即对速度自相关函数求积分Green-Kubo法和对均方位移求斜率的Einstein法[16]。现利用爱因斯坦关系式求得
式中:KB——玻尔兹曼常数
m——粒子质量
ri(t)、ri(0)、vi(0)、vi(t)——分子i在t时刻和0时刻的位移矢量和速度矢量
利用Einstein方法根据MSD和时间(图2)的导数求解水分子的扩散系数。
水分子在碳纳米管中扩散时,碳纳米管的孔径小于水分子之间的平均自由程,则水分子对碳纳米管壁的碰撞,较之水分子间的碰撞要频繁的多,这种扩散称之为Knudsen扩散[17]。Knudsen扩散的通量方程为:
式中:JA——组分A的扩散通量
CA——组分A的浓度
r——扩散距离
DK——Knudsen扩散系数
由气体分子运动论导出Knudsen扩散系数的计算公式:
式中:r——孔半径,cm
T——热力学温度,K
R——气体常数
M——扩散组分的分子量
表3显示的是利用MD方法得出的水分子的自扩散系数,现可利用Knudsen扩散系数的计算公式算出在不同管径的碳纳米管中水分子的自扩散系数(如表4)。
表3给出了不同管径下水分子的自扩散系数,从表3可以看出随着碳纳米管管径的增大,水分子的自扩散系数变大,表明管径越小水分子在碳纳米管中的扩散越不明显,扩散能力几乎丧失。同时由于受限的作用水分子在管内形成了有序的结构,并和管壁之间存在强烈的相互作用,这些都阻碍了水分子在管内的扩散。表4则是利用Knudsen扩散计算公式算出不同管径中水分子的扩散系数。Knudsen扩散区域中的扩散系数与碳纳米管的孔径大小成比例,这是由于随碳纳米管孔径的增大,水分子与管壁的碰撞频率减小,从而使得进入碳纳米管孔道的水分子数目增多,扩散系数增大。由此得水分子在不同管径的碳纳米管中的扩散属于Knudsen扩散。
摘要：利用分子动力学模拟研究了常温常压下受限于(8,8)(9,9)和(10,10)单壁碳纳米管中的水分子,对受限水分子平衡体系的状态和径向分布函数等静态性质进行了分析。结果显示,在不同管径的碳纳米管内部,水分子的微观结构都有高度的有序性。通过做水分子在碳纳米管内的均方位移与时间关系图,利用Einstein法算出不同管径中水分子的扩散系数,然后利用Knudsen扩散公式计算出的水分子的扩散系数,对种方法的计算结果进行比较。