动力学模拟论文

2024-08-06

动力学模拟论文（共10篇）

动力学模拟论文篇1

0 引言

杜仲胶(TPI)是一种天然高分子材料,国际上习称古塔胶(Gutta-Percha)或巴拉塔胶(Balata)。它与天然橡胶(NR)的化学组成相同(C5H10),但二者的分子链构型不同,导致其性状迥然不同。天然橡胶是优良的弹性体,而杜仲胶则是硬质塑料[1]。杜仲胶单纯用作硫化弹性体时其综合性能并不理想,通常与橡胶或塑料共混并用。杜仲胶的共混主要是指杜仲胶与天然橡胶、顺丁橡胶、聚乙烯等共混。实验证明,杜仲胶与橡胶共混,可以降低混炼胶的生热性,改善其焦烧特性,提高生胶强度,增强动态疲劳性能[2],混炼温度会明显降低;通过改变配方比例,可以得到性能变化范围较宽的新型材料[3]。

聚合物与聚合物之间的相容性影响共混物的形态和力学性能。提高聚合物与聚合物各组分之间的热力学相容性,可以防止共混物在加工和使用过程中发生聚结或相分离,改善产品的性能[4,5,6]。天然橡胶与杜仲胶的相容性可以用溶解度参数和Huggins-Flory相互作用参数来判断。通常采用浊度法和溶胀法获得溶解度参数,但这些方法或费时费料或存在某些局限性。与实验研究相比,分子模拟(MD)技术是一种更直接的研究方法,近年来已有越来越多的科研小组开始应用分子模拟方法研究共混物的相容性,并取得了较大的进展。但是对于杜仲胶与橡胶共混的相容性研究尚鲜见报道。因此,本研究拟从模拟溶解度参数与Huggins-Flory相互作用参数获得杜仲胶与橡胶共混物的相容性。

分子模拟属于微观模拟方法,可模拟的空间尺度小于100Å,时间尺度为纳秒级, 而由于高分子的长链结构,使相分离形成的特征区尺寸在几百埃,相应的动力学过程中时间尺度从微秒到毫秒或更长,远远超出分子模拟的范围[7],因此要进一步考察杜仲胶与天然橡胶的相态结构还需要借助耗散粒子动力学(DPD)模拟。本研究用MD模拟方法对体系的溶度参数、内聚能密度和共混物分子间的Flory-Huggins作用参数等进行模拟计算,预测共混物的相容性,揭示TPI与NR分子间相互作用的本质,然后将MD模拟的结果作为DPD模拟的输入参数,实现TPI/NR的微、介观多尺度模拟。

利用MD模拟方法可以预测聚合物共混物的性能,从而为共混物配方设计提供理论指导。朱伟[8]用MD方法对钝感炸药TATB及其与氟橡胶所构成的高聚物粘结炸药PBX的弹性系数、模量进行模拟计算,发现在TATB中添加少量氟橡胶能有效改善其力学性能。付一政[9]模拟计算了高分子粘结剂端羟基聚丁二烯(HTPB)与增塑剂癸二酸二辛酯(DOS)、硝化甘油(NG)的力学性能,结果表明,添加DOS增塑剂后使HTPB的弹性模量(E)、体积模量(K)和剪切模量(G)下降, 材料刚性减弱, 柔性增强, 力学性能得到改善。本研究拟用MD模拟方法预测TPI/NR共混物的力学性能,以期为筛选出杜仲胶替代部分天然橡胶的配方提供理论指导。

1 模型构建与模拟方法

1.1 构建分子模型

运用Materials Studio 4.0(MS4.0)软件包中的Visualizer模块分别构建杜仲胶及天然橡胶的分子模型,如图1所示。在分子模拟中,分子链聚合度选择越大,其模拟的计算结果和统计结果越精确,但模拟所需要的时间就越长,对于计算机的性能要求也就更高,为此需要选择合适的分子链聚合度。为了确定能够真实反映聚合物的最小聚合度,改变聚合物的长度,计算聚合物不同长度时的溶度参数(δ),当溶度参数不再随着聚合度的改变而变化时,说明该聚合度可以保证计算结果的准确性[10,11,12]。

采用MS4.0中的Visualizer模块分别构建聚合度为5、10、15、20、25、30、35、40、45、50的杜仲胶与天然橡胶的单链模型[13],经最小化处理之后,依据在298K、1.01325×105Pa下杜仲胶与天然橡胶的密度,分别搭建三维周期性边界条件的无定形高分子模型,对所建立的模型进行能量最小化,消除建造模型中产生的局部不合理结构,为下一步的分子动力学模拟提供较为合理的结构。

1.2 MD模拟细节

利用Smart minimization方法对所构建的无定形分子模型首先进行结构优化,然后将优化后的模型进行MD模拟。在整个模拟过程中,温度为298K,温度控制方法采用Andersen方法;位能采用球形截断——长程加和修正法,截断半径取0.85nm,样条宽度(Spline width)取0.1nm,缓冲宽度(Buffer width)取0.05nm,时间步长为1fs。先进行50ps等温等容(NVT)系综的MD模拟,再进行50ps等温等压(NPT)系综的MD模拟;计算溶解度参数、内聚能密度等相关性能。

确定聚合度之后,杜仲胶与天然橡胶的高分子链利用Visualizer模块共混,混合物的密度由单体密度计算。混合模型经能量最小化之后,在1.01325×105Pa下,其他参数不变,进行150ps的恒温恒压(NPT)系综的MD模拟,分析计算不同比例共混物的溶解度参数、内聚能密度、力学性能等。

1.3 DPD模拟细节

采用DPD方法对其二元共混的相态进行研究,其中,最关键的步骤是建立真实分子结构和介观参数之间的联系。Groot等[14]认为,对于聚合物分子的粗粒化,一个珠子所代表的应该是其分子链的持久长度。因此,采用式(1)对真实聚合物分子链进行粗粒化,得到粗粒化模型[15]。

undefined

式中:N为聚合物的聚合度,Cn为特征比,可以通过MS 软件Synthia 模块求得TPI高分子链的Cn为5.15,NR高分子链的Cn为5.45,则得到每条TPI、NR分子链上含有的DPD粗粒化珠子数目,之后再用弹簧将相邻DPD珠子链接起来,就得到了线性TPI、NR分子的DPD粗粒化模型。

在DPD模拟中,最重要的参数是DPD珠子之间的相互排斥力参数(aii和aij)。在高分子体系中,同种粒子的排斥力参数与体系的数密度有关,为了使纯DPD珠子体系的压缩系数与真实体系中水的压缩系数相匹配,并且流体的涨落能够被正确地描述,aii选取应遵循式(2)[16]:

undefined

式中:kBT为能量单位,ρ为密度。在本研究的模拟中,选取kBT=1,ρ=3。因此,相同珠子之间的最大排斥力即aii=25。TPI与NR分子的特征比都接近于5,当本研究选择映射因子为5时,得到DPD珠子与珠子之间的χij′=χij×5[17],不同珠子之间的排斥力常数即aundefined可通过DPD珠子与珠子之间的χij′得到,即:

aij=aii+3.497χij′ (3)

所有的介观模拟都是在MS 软件中的DPD模块中计算,模拟中所用的盒子尺寸是20.0nm×20.0nm×20.0nm,体系数密度为3,因此,体系包含24000个DPD粒子。研究过程中,通过改变均聚物A、B的体积分数来研究相转变和相结构。由于每个体系的组分不一样,所以体系中包含的各种分子的个数也就不一样。为了便于模拟,将长度rc、质量m以及温度kBT的单位都设置为1。弹簧常数C=4.0,耗散因子γ=4.5,噪声强度σ=3,由于所采用的积分方法为Velocity-Verlet 算法,为保证在此积分方法下能够使温度得到很好地控制,所选用的时间步长为0.05。每模拟均从无序态开始,运行至少 5×104直至相结构达平衡。

2 结果与讨论

2.1 相容性判断

为了得到高分子共混的溶解性,可借助于溶度参数这个物理量。分子的溶度参数δ是内聚能密度的平方根,有:

undefined

根据杜仲胶、天然橡胶不同聚合度的溶解度参数,可以确定能够真实反映聚合物的最小聚合度。如图2所示,杜仲胶的聚合度从5增加到15,其δ迅速减小,从15开始缓慢回升,到了30开始下降,35、40、45、50时的δ基本稳定,达到最小值,约为16.6(J/cm3)1/2,从而可以确定N(TPI)=35;天然橡胶的聚合度从10到20大幅下降,之后较为平缓,在30、45时略有偏差,δ保持在16.2(J/cm3)1/2左右,则可以确定N(NR)=20。所得结果见表1,模拟得到的δ与文献[18]的值基本接近。

若二元共混体系[19]足够平衡时,则共混体系的能量ΔEmix可用式(5)计算:

undefined

式中:ϕA与ϕB分别是组分A、B的体积分数。

Huggins-Flory相互作用参数χ1,2 采用式(6)、(7)计算:

undefined

式中:Vmon为单体单元体积,在本研究中Vmon=68.69cm3/mol,T为绝对温度,R为气体常数,计算结果见表2。

χ1,2的正负并不能绝对地判断共混物的相容性,因此,需要引入临界Huggins-Flory相互作用参数[20](χ1,2)C,两种聚合物相容的条件是χ1,2≤(χ1,2)C。本研究中,经计算得到(χ1,2)C =0.1963,则由图3可知,不同混合比的χ1,2都在(χ1,2)C的下方,由此说明各种混合比的TPI/NR共混物均可相容。

2.2 DPD模拟结果分析

径向分布函数[21]g(r)是反映材料微观结构的特征物理量,它表示在一个分子周围距离为r的地方出现另一个分子的概率密度相对于随机分布概率密度的比值。分子内径向分布函数可以提供模拟体系有序度的信息,一般出现大于0.3nm 的峰表明分子链长程有序,属于结晶体系;而出现小于0.3nm的峰则表明分子链短程无序,属于无定形结构。图4为TPI、NR纯物质的分子内和分子间的径向分布函数。从图4可见,在0.3～0.7nm范围内两种纯物质均没有出现强的峰,说明所构建的分子结构属于无定形结构。分子间径向分布函数可以揭示非键原子间相互作用方式和本质,氢键作用范围为0.26～0.31nm,范德华作用范围为0.31～0.50nm。由图4可见,TPT、NR分子间主要作用方式为范德华作用。二者的径向分布函数曲线非常接近,根据相似相容的原理,可以更进一步说明其共混体系的相容性较好。

图5为不同比例的TPI/NR共混物的等密度图。为了更加直观地观察到体系内部的相结构,在等密度图中灰色表面包围的区域为TPI相,空白区域为NR相。从图5中可以看到,当共混体系达到平衡之后,TPI相分散比较均匀,体系并未发生同相归并和相分离,说明其相容性很好。

2.3 力学性能分析

对于MD模拟所得体系的轨迹文件,基于静态力学分析原理,利用MS 程序对体系实行多次小变型的单轴拉伸与纯剪切形变操作后,在原子水平上由维利公式求得内应力张量,于是由广义虎克定律等式两边对应变分量求一阶偏导数得到弹性系数矩阵。然后分别通过最小二乘法拟合拉伸应力应变曲线,求得拉伸应力与拉伸应变之比和泊松比横向应变与纵向应变之比的负值的平均值,进而可求得其它有效各向同性力学性能,如剪切模量、体积模量和拉梅系数。各向同性材料的模量可以用拉梅系数表示[22]:

undefined(拉伸模量)

undefined(体模量)

G=μ (剪切模量)

拉伸模量、体模量和剪切模量均可用来衡量材料的刚性强弱,体模量与剪切模量的比值即K/G可用来衡量体系的韧性。

因为聚合物材料的弹性模量可以理解为材料的弹性系数,所以形状记忆材料的热收缩性可以用材料的弹性模量E来特性化,E越大,形状记忆性能越好[23]。通过上述力学分析,可以用拉伸模量E来评价形状记忆高分子材料的优劣。

从图6(a)可以看出,随着TPI质量分数的增加,模量呈上升趋势,到了36.8%达到最大,至63.6%时有一突变,模量急剧下降,之后又略有回升。从图6(b)可以看到,36.8%时其韧性最好。此结论也符合实验结果[24,25,26],李良萍等[24]研究了杜仲胶与天然橡胶共混硫化胶的静态力学性能及动态伸张疲劳性能,结果表明,在共混硫化胶体系中,杜仲胶质量分数在40%以下能较好地保持其优良的力学性能。由此可以确定1/3(TPI/NR)共混体系的弹性模量及其韧性均优于其他组分,可以为后续的三组分共混提供依据,进一步从分子模拟的角度来研究杜仲胶/天然橡胶/低密度聚乙烯共混物的形状记忆性能[27]。

3 结论

(1)通过MD模拟杜仲胶、天然橡胶在5～50聚合度下的溶解度参数,确定能够真实反映聚合物的最小聚合度N(TPI)=35、N(NR)=20,得到与实验值较为吻合的溶解度参数δTPI=16.6(J/cm3)1/2、δNR=16.2(J/cm3)1/2。

(2)利用MD模拟获得Huggins-Flory相互作用参数,其值均小于临界Huggins-Flory相互作用参数,从微观尺度上可以判断共混物体系的相容性比较好。

(3)通过比较TPI、NR纯物质的分子内和分子间径向分布函数发现,二者分子链短程无序,属于无定形结构,同时揭示了TPT、NR非键原子间主要作用方式为范德华作用。二者的径向分布函数曲线非常接近,根据相似相容的原理,可以更进一步说明其共混体系的相容性。

(4)通过DPD模拟获得不同比例的TPI/NR共混物的等密度图,更加直观地观察到体系内部的相结构,经过分析预测了共混体系的相容性。结果表明,各种比例的TPI/NR共混物的相容性均较好,与MD模拟的结论一致。 MD和DPD模拟的结合实现了TPI/NR的微-介观多尺度模拟。

(5)通过分析比较不同比例的共混体系的静态力学性能,确定了1/3(TPI/NR)共混体系的弹性模量及其韧性均优于其他组分,其结论与实验结果一致。

致谢感谢西北工业大学高性能计算中心的支持!

动力学模拟论文篇2

磁场对水内能作用的分子动力学模拟

将水系统视为正则系综,采用分子动力学模拟方法计算了磁场条件下水的内能、比热和径向分布函数.结果表明,磁场可以影响径向分布函数g(r)的`分布,使水的结构发生变化,从而导致水的内能、比热发生改变.水的内能、比热、径向分布函数伴峰的高度与磁场强度的关系均呈多极值特征,且在常温下磁感应强度为0.25 T时,磁场作用最明显.

作者：张军周开学作者单位：石油大学应用物理系,山东东营,257061 刊名：石油大学学报(自然科学版) ISTIC EI PKU英文刊名：JOURNAL OF THE UNIVERSITY OF PETROLEUM，CHINA(EDITION OF NATURAL SCIENCE) 年，卷(期)： 26(2) 分类号：O482.5 关键词：磁场水分子动力学模拟内能比热径向分布函数

动力学模拟论文篇3

(上海海事大学商船学院，上海 201306)

低温等离子体技术净化船舶柴油机尾气的化学反应动力学模拟

孙永明，夏文虎，张勤

(上海海事大学商船学院，上海 201306)

针对目前尚没有一种净化方法能同时去除船舶柴油机尾气中NOx，SO2和炭粒的现状，提出将低温等离子体技术应用于船舶柴油机尾气净化中.利用MATLAB，根据化学反应动力学原理建立关于低温等离子体技术净化NOx和SO2的微分方程组，对船舶尾气中NOx和SO2的净化效果进行模拟.在此基础上，进一步研究尾气中的SO2含量和O2含量对尾气处理效果的影响.理论研究结果说明此方法是可行的.利用低温等离子体技术净化船舶柴油机尾气具有较好的发展前景.

低温等离子体；化学反应动力学；船舶柴油机；尾气净化

0 引言

为降低航运成本，远洋船舶主机的燃料一般为劣质燃油.船舶使用的燃油的平均硫含量为2.7%.船舶主机的三大污染物，NOx，SO2及炭粒[1]，对海洋环境造成极大的污染，其中NOx更是污染空气的主要物质之一.这一问题越来越引起国际海事组织(International Maritime Organization, IMO)和国际社会的重视.国际上对船舶柴油机尾气处理的研究也非常多，但是尚未有一种方法能同时去除尾气中的三大污染物并能实际应用于船上.等离子体技术应用于船舶尾气处理是一个新的研究方向，曾有学者提出等离子体放电是去除NOx的一种有效方法[2]，但近年来研究多采用脉冲电晕放电[3-5]和介质阻挡放电[6-8]的方法，它们都通过某种方式将NOx还原为对环境友好的N2.

关于等离子体氧化柴油机尾气中的颗粒物(炭粒)的大量文献说明：反应过程中OH和O自由基对PM的氧化起到重要作用.[14]本文根据化学反应动力学原理，以MATLAB为研究工具模拟用低温等离子体技术净化船舶柴油机尾气中NOx和SO2的过程.通过对模拟结果的分析，并参考相关文献对低温等离子体氧化炭粒的研究，发现低温等离子体技术结合海水吸收的方法不仅能够实现船舶柴油机尾气的脱硫脱硝和炭粒去除，而且只消耗电能，不需要价格昂贵的催化剂，具有成本低的特点，这正是船舶所有人追求的.而且利用该方法净化船舶尾气所需的设备体积不大，不受船舶机舱空间狭小的限制.若能成功将此方法应用于船舶柴油机尾气处理，则在将来IMO对船舶尾气排放要求越来越高的大环境下，该方法具有广阔的使用前景.

1 等离子体产生强氧化性粒子的反应机理

船舶主机基本上是大型低速柴油机，其排气中O2约占14%，H2O约占5.1%，产生强氧化性粒子的反应机理如下.

电极间加高电压，产生电晕放电，电晕放电又会产生高能电子(5～20 eV)，这些高能电子会激活或电离尾气中含量较多的H2O，O2，N2分子，具体反应如下.

(1)高能电子与分子碰撞：

(2)激发态氧原子与分子碰撞：

(3)此外，在使自由基消失的众脱除反应中有一反应：

由以上3个步骤的反应可见强氧化性粒子的产生过程：首先高能电子碰撞中性分子产生活性粒子，如(1)中反应所示；由于这些活性粒子中有极不稳定的激发态原子存在，如O(1D)，其与中性分子继续发生反应产生二次自由基，如(2)中反应所示；而O自由基又与O2分子反应生成臭氧O3，如(3)中反应所示.经过这3个步骤的反应，产生大量的强氧化性活性粒子OH，O自由基和O(1D)，O3等.

2 低温等离子体净化NOx和SO2的化学反应动力学数值模拟

在高电场下，电极间由于电晕放电产生的高能电子(5～20 eV)激活H2O及O2等，产生的OH和O等活性粒子将SO2和NOx氧化，其中有3个过程：电子与分子碰撞、激发态氧原子与分子碰撞和脱除反应自由基消失.相关的化学反应及化学反应速率常数[15-16]见表1.

2.1 反应动力学模型的建立

大型船舶低速柴油机的尾气中包括14%的O2，76.2%的N2，4.5%的CO2，5.1%的H2O，0.14%的NOx以及0.06%的SO2等.

忽略低温等离子体的磁效应，忽略中性原子和离子的运动对等离子体反应过程的影响，假设放电间隙内粒子均匀分布、电场均匀分布[17].

表1 相关的化学反应及化学反应速率常数

根据表1中所列出的反应，共有e，e*，O2，O，H2O，H，OH，O(1D)， O3，NO，NO2，HNO3，NO3，N2O5，SO2，SO3， HSO3，H2SO4，N2和N等20种粒子.yi(i=1，2，…，18)分别表示后18种粒子的浓度(每立方厘米分子数).根据反应动力学通式建立的初值常微分方程组如下.

dy1/dt=-k1Cey1-k3Cey1-k6y17y1y2-k19y1y1y2+k9y7y8+k10y7y9

dy2/dt=2k1Cey1+k3Cey1+k5y6y17-k6y17y1y2-k7y2y8-k16y2y13-k19y1y1y2

dy3/dt=-k2Cey3-k4y3y6-k12y3y12-k15y3y14

dy4/dt=k2Cey3

dy5/dt=k2Cey3+2k4y3y6-k8y5y9-k13y5y13-k14y5y15-k18y5y18

dy6/dt=k3Cey1-k4y3y6-k5y6y17

dy7/dt=k6y17y1y2+k19y1y1y2-k9y7y8-k10y7y9

dy8/dt=-k7y8y2-k9y8y7+k18y5y18

dy9/dt=k7y8y2-k8y9y5+k9y7y8-k10y9y7-k11y9y11

dy10/dt=k8y5y9+k12y3y12

dy11/dt=k10y7y9-k11y11y9

dy12/dt=k11y11y9-k12y3y12

dy13/dt=-k13y5y13-k16y13y2

dy14/dt=-k15y14y3+k16y2y13

dy15/dt=k13y5y13-k14y5y15

dy16/dt=k14y5y15+k15y3y14

dy17/dt=-k17Cey17

dy18/dt=2k17Cey17-k18y18y5

其中：O2的初始浓度为3.763×1018/cm3，N2为2.047 9×1019/cm3，H2O为1.371×1018/cm3，NO为3.829 7×1016/cm3，NO2为2.016×1015/cm3，SO2为1.612 5×1016/cm3，电子浓度Ce为1013/cm3，其余粒子的初始浓度为0.[18]

2.2 模拟结果

图1 NO浓度随时间的变化

图2 NO2浓度随时间的变化

图3 SO2浓度随时间的变化

应用MATLAB软件编写微分方程的计算程序，并用较高精度的ode45法(四五阶Runge-Kutta法)求解该初值下的常微分方程组，得到NO，NO2及SO2的浓度随时间变化的曲线，分别见图1～3.由图1～3可知：因NO的初始浓度比较高，NO2的初始浓度很低，所以刚开始NO在等离子体下反应生成NO2的速度比NO2反应生成HNO3的速度快得多，使得初始时NO2的浓度不但不下降反而快速上升；随着反应的进行，因为生成NO的速度较慢，NO的浓度下降导致生成NO2的速度也急剧下降，这时NO2的消耗速度又远大于其生成速度，所以NO2的浓度在达到顶峰后又迅速降低.结合化学反应方程式还可以知道，虽然有NO生成，但NO浓度较高时，NO被反应掉的速度大于其生成速度，最后NO的浓度保持在较低水平.由图1～3还可以发现，3种污染物被去除的速度都很快，虽然柴油机排放尾气的速度也非常快，反应装置又有体积限制(船舶上机舱空间有限)，停留在反应器中的时间非常短，但最后净化尾气中污染物的效果还是很好的.

图4 SO2的初始浓度为0时NO浓度随时间的变化

图5 SO2的初始浓度为0时NO2浓度随时间的变化

分析柴油机燃油中硫含量(尾气中SO2的初始浓度)对处理尾气中NOx的影响，当SO2的初始浓度为0时，得到NO和NO2的浓度随时间变化的曲线，分别见图4和5.由图4和5可知：当使用含硫量低的高质燃油时，NO2浓度的峰值(1.48×1016/cm3)较之前的峰值(1.85×1016/cm3)明显降低；同时，大型船舶柴油机排放尾气的速度非常快，尾气在反应器中停留的时间非常短，尾气中NO2的初始浓度并不高，刚开始时会上升，如果峰值过高就会增加NO2排放超标的风险；使用含硫量低的燃油时，NO2浓度到达较低值时所用时间明显减少，因而使用含硫量低的燃油有助于提高NOx的去除率.这是因为SO2与活性粒子反应消耗了活性粒子，使得参与NOx去除反应中的活性粒子的量降低.

分析废气中的O2浓度对尾气中NOx和SO2去除效果的影响，当O2的初始浓度为3.763×1017/cm3时，得到NO，NO2及SO2的浓度随时间变化的曲线，分别见图6～8.

图6 O2浓度较低时NO浓度随时间的变化

图7 O2浓度较低时NO2浓度随时间的变化

图8 O2浓度较低时SO2浓度随时间的变化

由图6～8可以明显看出：当O2浓度较低时，NO，NO2及SO2的去除速度都明显降低，达到与O2浓度较高时相同的净化效果所花费的时间明显增加.从图中可观察到O2浓度较低时反应变得特别慢，污染物的浓度相对于O2浓度较高时变化非常缓慢，这是由于活性粒子O，O(1D)，OH和O3等的生成都直接或间接地需要O2的参与，而这些活性粒子又是净化柴油机尾气中污染物NOx和SO2的必需物质，它们的浓度大小直接关系到净化反应的快慢.O2浓度对NO2浓度的影响最大：O2浓度降低使NO2的峰值向右移动，达到峰值的时间高了一个数量级，NO2达到较低的浓度所需的时间也高了一个数量级.前面已经分析过尾气在反应器中停留的时间不能太长，若尾气中O2浓度较低，则尾气在反应器中停留时间过短，最终通过反应器的尾气中NO2的浓度可能比初始时NO2的浓度还要高.

3 结论及展望

(1)利用低温等离子体技术净化船舶柴油机尾气中的NO，NO2和SO2，其去除速度都很快，能达到所需的净化效果.(2)在反应过程中有NO生成，但当NO浓度稍高时，NO被反应掉的速度大于其生成速度，因而NO浓度最终维持在一较低水平；NO2浓度初始时会增加，但达到峰值后会迅速降低.(3)使用含硫量较低的燃油时，NO2浓度的峰值明显降低，NO2浓度达到最低值所用的时间也明显减少.(4)尾气中O2浓度降低会使NO，NO2和SO2的去除速度明显降低，净化所需时间明显增加.

国内外对低温等离子体处理船舶柴油机尾气的研究尚很少，还没有相关的试验研究.本文用化学反应动力学数值模拟的方法，对该方法净化船舶柴油机尾气中的NOx和SO2进行研究，理论研究结果说明此方法是可行的，但尚需实验来验证，有条件的机构可做相关实验研究.

[1]裴梅香. 低温等离子体技术及其在柴油机排气处理中的应用[J]. 环境污染治理技术与设备, 2004, 5(5): 56-60.

[2]宿鹏浩. 用TiO2光催化剂与直流电晕放电结合去除NOx[J].上海海事大学学报, 2011, 32(4): 80-84.

[3]YU Chunjiang, XU Fei, LUO Zhongyang,etal. Influences of water vapor and fly ash addition on NO and SO2gas conversion efficiencies enhanced by pulsed corona discharge[J]. J Electrostatics, 2009, 67(6): 829-834.

[4]商克峰, 李国锋, 吴彦, 等. 添加剂对电晕放电尾气脱硝效率和NOx转化影响[J]. 大连理工大学学报, 2007, 47(1): 21-25.

[5]WANG Wenchun, ZHAO Zhibing, LIU Feng,etal. Study of NO/NOxremoval from flue gas contained fly ash and water vapor by pulsed corona discharge[J]. J Electrostatics, 2005, 63(2): 155-164.

[6]高旭东, 孙保民, 肖海平, 等. 介质阻挡放电脱除NOx反应器的评价方法及运行流量特性分析[J]. 中国电机工程学报, 2010, 30(1): 27-32.

[7]杜伯学, 刘弘景, 王克峰, 等. 介质阻挡放电产生低温等离子体除去NOx的实验研究[J]. 高电压技术, 2009, 35(9): 2186-2192.

[8]MARIO M S, AXEL V, ESTEBAN S,etal. Design of a DBD wire-cylinder reactor for NOxemission control: experimental and modeling approach[J]. J Cleaner Production, 2008, 16(2): 198-207.

[9]龚大国. 等离子体汽车尾气治理技术[J]. 重庆环境科学, 2003, 25(2): 28-32.

[10]郭鲁钢. 海水脱硫技术现状[J]. 海洋技术, 2006, 25(3): 10-14.

[11]李英峰. 海水尾气脱硫工艺排水对区域海水的影响[J].电力环境保护, 2000, 16(3): 10-13.

[12]宋晓东. 浅谈尾气脱硫工艺排水对海洋环境的影响[J]. 山东电力高等专科学报, 2000, 3(3): 44-48.

[13]吴来贵. 深圳西部电厂4号机组海水脱硫系统监测分析[J]. 热能动力工程, 2003, 18(2): 200-202.

[14]THOMAS S, MARTIN A, RAYBONE D,etal. Non-thermal plasma aftertreatment of particulates-theoretical limits and impact on reactor design[C].SAE Technical Paper, 2000: 1926. DOI: 10.4271/2000-01-1926.

[15]TOKUNAGA O, SUZUKI N. Radiation chemical reaction in NOxand SO2removal from flue gas[J]. Radiat Phys Chem, 1984, 24: 145-165.

[16]MATZING H. Chemical kinetics of flue gas cleaning by irradiation with electron[J]. Adv Chem Phys, 1991, 80: 315-402.

[17]张文豪. 低温等离子体技术净化柴油机尾气NOx的化学反应动力学模拟研究[J]. 拖拉机与农用运输车, 2008, 35(1): 79-81.

[18]KULIKOVSKY A A. Production of chemically active species in the air by a single positive streamer in a nonuniform field[J]. IEEE Trans Plasma Sci, 1997, 25(3): 439-446.

(编辑赵勉)

Chemical reaction kinetics simulation of marine diesel engine flue gas cleaning by low temperature plasma technology

SUN Yongming, XIA Wenhu, ZHANG Qin

(Merchant Marine College, Shanghai Maritime Univ., Shanghai 201306, China)

Considering that there exist no cleaning methods that can remove NOx, SO2and carbon particles together from marine diesel engine flue gas, it is proposed that the low temperature plasma technology is put forward for cleaning marine diesel engine flue gas. The differential equations of removing NOxand SO2by the technology are established by MATLAB according to the chemical reaction kinetics. The effects of removing NOxand SO2from marine diesel engine flue gas are simulated. On this basis, the influences of SO2and O2contents in flue gas on treatment effect are further studied. The theoretical research results show that the method is feasible. It has better development prospect to use the technology to clean marine diesel engine flue gas.

low temperature plasma; chemical reaction kinetics; marine diesel engine; flue gas cleaning

10.13340/j.jsmu.2015.02.015

1672-9498(2015)02-0079-05

2014-08-15

2014-10-28

孙永明(1962—)，男，浙江绍兴人，副教授，硕士，研究方向为船舶辅机工程，(E-mail)ymsun@shmtu.edu.cn

U664.121； TK421.5

荡秋千动力学过程的数值模拟篇4

武际可[1]、刘延柱[2]将秋千系统简化为一个摆长周期性变化的单摆模型,即具有两个自由度的单摆模型,生动解释了秋千系统越荡越高的力学机理与功能转化,贺承绪[3]利用单摆模型从能量的角度分析了秋千系统越荡越高的原因,Case等[4]将秋千系统简化为具有两个自由度的三质点模型,用一个三球哑铃模拟荡秋千者,通过让哑铃前后摆动,使人体的质心升高或降低来再现秋千的振荡过程.

本文将秋千系统分为人体和秋千两个部分,根据人体在荡秋千过程中的动作特点,将人体简化为由三个刚体构成的子系统,而将秋千分别简化为不计质量的刚性直杆和不可伸长的柔索,应用多体系统动力学理论数值模拟了荡秋千的动力学过程,证明了只有在人体上下屈伸的频率接近秋千系统摆荡频率的两倍时才能使整个系统越荡越高,此外,改变人体上下屈伸的初始相位(即起始时刻)可以得到完全相反的结果,即系统不仅不会越荡越高,反而会越荡越低.

1 秋千系统的四刚体力学模型

根据人体在荡秋千过程中屈膝下蹲和挺身直立的动作特点,将人体简化为3个刚体,将秋千简化为一个刚体,如图1所示.其中B1为小腿(包括脚),B2为大腿,B3为躯干,B4为不计质量的刚性直杆,O,O1,O2,O3均为旋转铰,O为秋千的悬挂点,O1为脚与秋千的连接点,O2,O3分别为膝关节和髋关节,O4表示躯干通过手与秋千的连接,一方面O4在秋千上的位置固定不变(代表手握秋千的位置),另一方面O4相对于躯干可以沿躯干轴线移动.假设OO1的长度为5 m,其中O1O4段长为1.3m,人体的几何和质量参数如表1所示,其中各部分的质心位于自身轴线(图1中虚线)的中点,所列转动惯量均为通过各自质心且垂直图1所在平面的轴线的转动惯量,小腿和大腿的质量和转动惯量为左右两侧之和.

选取描述系统位置的坐标为q=(φθ1θ2θ3)T,其中φ为秋千摆荡的幅角,θ1为小腿与O1O4之间的夹角,θ2和θ3分别为膝关节和髋关节控制变量.由于系统具有3个自由度,故上述坐标之间存在一个约束方程

根据多体系统动力学的理论,可得到系统如下形式的动力学微分方程[5]

其中A(q)为4×4的广义质量对称矩阵,Φq(q)为约束方程(1)的4×1的雅可比矩阵,λ为拉格朗日乘子,B(q,,M1,M2)为4×1的矩阵,其中包含了重力及膝关节和髋关节的驱动力矩M1,M2对方程的贡献.

假设系统初始静止,φ0=10°,膝关节控制变量θ2的规律为

其中,t1为膝关节开始屈伸的时刻,ω为膝关节屈伸的频率.若在荡秋千的过程中,躯干的轴线始终与O1O4段重合(通过髋关节实现),则髋关节控制变量θ3的规律为

假设t1=7π/2ω,ω=2ω0,其中ω0为秋千摆荡的频率,本文取ω0=2x/4.16.数值计算的有关结果示于图2,图3中,其中图2为平面上的相轨迹,可见相轨迹由内向外扩展,说明秋千越荡越高.图3为系统的机械能随时间的变化情况,系统的机械能总体趋势是增加的,但人体在高处屈膝下蹲使系统的机械能减少,而在低处挺身直立使系统的机械能增加,故系统的机械能是波动的,这些均与有关文献中的定性分析结论一致[1,2,3].

为了分析膝关节屈伸的频率对秋千摆荡的影响,现分别假设ω=2ω0和ω=ω0,其它条件不变,得到平面上的相轨迹如图4所示,可见只有当膝关节屈伸的频率为秋千摆荡频率的两倍时,秋千才越荡越高.事实上,作者所做的计算表明,当人体上下屈伸的频率接近秋千摆荡频率的两倍时,就能使秋千越荡越高.

进一步地,若取φ0=30°,人体开始屈伸的时刻分别为t1=7π/2ω和t1=5π/2ω,其它条件不变,得到平面上的相轨迹如图5所示,可见如果人体在秋千上反其道而行之,即在高处挺身站立,在低处屈膝下蹲,相轨迹由外向内收缩,秋千反而会越荡越低.也就是说,改变人体上下屈伸的初始相位(即起始时刻)可以得到完全相反的结果.

2 秋千系统的三刚体-柔索力学模型

与上述秋千系统的四刚体力学模型不同的是,现将秋千简化为不可伸长的柔索B4(由OO4和O1O4两段构成),如图6所示.此时选取描述系统位置的坐标为q=(φψθ1θ2θ3)T,其中ψ,为柔索O1O4段与OO4段之间的角度,其它坐标的物理意义同上.由于系统具有4个自由度,上述坐标之间存在一个约束方程,根据多体系统动力学的理论,同样可以得到类似于方程(1)和(2)的约束方程和动力学微分方程,不过其中A(q),Φq(q)和B(q,,M1,M2)分别为5×5,5×1和5×1的矩阵.计算得到的平面上的相轨迹如图7所示,由于柔索不能承受弯矩,故ψ角的变化使得秋千的摆荡不够平稳,但相轨迹的总体趋势仍是由内向外扩展的,而得到的系统机械能的时间历程与图3基本相同,这些都说明秋千仍能越荡越高.进一步的计算表明,改变人体上下屈伸的频率和初始相位,可以得到与上面四刚体力学模型完全类似的结论.

参考文献

[1]武际可．从荡秋千说开去——漫话共振．力学与实践，2003，25(2)：75～78

[2]刘延柱．再谈荡秋千——兼谈自激振动．力学与实践，2007，29(3)：92～93

[3]贺承绪等．关于荡秋千的能量分析．物理与工程，2005，15(1)：25～26

[4] Case WB等．对于荡秋千的动力学分析．数理译丛，1991，1：61～68

[5]洪嘉振．计算多体系统动力学．北京；高等教育出版社，2002

动力学模拟论文篇5

超临界NaCl水溶液的分子动力学模拟?

采用分子动力学模拟的方法对超临界NaCl水溶液的.微观结构进行了研究.模拟发现在所研究超临界条件下,密度的变化比温度的变化对超临界NaCl水溶液的微观结构影响更大.温度及密度对Cl--H2O径向分布函数的影响比对Na+-H2O径向分布函数的影响要大.超临界条件下,各gNa+-Cl-在0.261nm处出现峰值,表明Na+、Cl-之间发生了离子的缔合.超临界条件下,随温度增加,缔合作用增强;随密度增加,缔合作用减弱.本文工作为建立可适用于超临界条件下的电解质热力学模型提供了依据.

作者：周健朱宇汪文川陆小华王延儒王钧作者单位：周健,朱宇,陆小华,王延儒,王钧(南京化工大学化工学院,江苏省化学工程与技术重点实验室,南京210009)

汪文川(北京化工大学化学工程学院,北京100029)

刊名：物理化学学报 ISTIC SCI PKU英文刊名：ACTA PHYSICO-CHIMICA SINICA 年，卷(期)：2002 18(3) 分类号：O645.164 关键词：氯化钠微观结构超临界流体分子动力学分子模拟电解质溶液离子缔合

山东省科技发展系统动力学模拟篇6

近年来国内外学术界对科技发展指标的模拟和预测进行了不少研究。国外学者运用动态逼近函数的处理方式和马氏决策过程模型来研究人才的需求量[1,2]。国内学者王维、董奋义、王建玲等运用灰色模型GM(1,1)分别对科技人才的需求和科技投入与科技产出的关系进行预测[3,4,5],为区域科技规划与决策的制定提供了科学的理论依据;杨月运用回归分析方法,对大量科技统计数据进行分析处理,确定各变量之间的关系,建立回归分析数学模型,对科技资源的需求进行预测[6];杨雪、常新宇、张建勇等把BP神经网络应用到科技人才需求预测中,通过确定科技需求能力的衡量指标,建立神经网络预测模型,预测了区域科技人才需求状况[7,8,9]。以上研究大多是对单个科技发展指标进行预测,没有把科技发展作为整体系统进行预测研究,也较少考虑科技发展与经济发展之间的动态互动关系。

本文把科技发展与经济发展作纳入一个系统,从科技发展内部入手,分析各因素之间的关系,建立起因果关联、信息反馈的系统模型,明确了系统中各内部要素之间相互作用的关系以及与外部系统之间的联系,以科技和经济的协调发展为基点构造了系统动力学模型。根据山东省历史统计数据和未来国民经济发展目标确定系统模型中的参数,通过调整决策变量,采用Vensim软件进行仿真计算,模拟了“十二五”期间山东省科技发展各项主要指标的变化趋势,得到科技发展速度的高、中、低3套方案,最后根据山东省具体情况,提出了科技发展的对策和建议。

1 科技发展系统动力学模型的建立

1.1 系统模型结构及主要素分析

在科技发展的模型中,本文从区域的总体性、综合性、关联性等角度分析了变量之间的因果反馈关系,从而建立起了科技发展与经济协调发展的关联结构。在结合科技发展综合评价指标体系的基础上,把整个科技发展系统主要分为3个子系统,分别为科技投入、科技产出和经济发展子系统。

在科技投入子系统中,企业R&D投入和政府科技投入是科技总投入的主要来源,政府科技投入的主要作用是促进和刺激私人企业增加科技投入,引导产业技术创新方向从而达到促进技术创新发展和调整产业结构的目的。企业是技术创新的主体,企业R&D投入是促进国家或地区经济增长、提升竞争力的重要条件之一。同时考虑到教育的投入对科技发展的促进作用,把反映教育投入水平的财政教育经费以及教育投入比例引入子系统。另外把反映全社会科技投入强度的指标全社会R&D投入占GDP的比重以及反映一个国家对外技术依赖程度的指标对外技术依存度作为辅助变量。

科技产出子系统中,把由于科技的发展,所带来的一系列收益作为指标变量。科技贡献率反映了一个国家或地区科技的进步对经济发展的促进程度。新产品销售收入反映了由于企业增加的科技投入为本企业所带来的销售效益。科技论文数量与专利申请量和授权量分别反映一个地区科技研发能力和创新能力。获国家科技成果奖代表了一个地区的核心科技竞争力。高新技术产业产值则体现了一个地区高新技术产业化的能力,从另一个侧面也反映了该地区的科技实力。

经济发展子系统中,科技发展的最终目的就是为了促进经济又好又快的发展。在经济发展模块中选取GDP、GDP增加量、GDP增长率、地方财政收入、固定资产投资额等能够反映一个地区经济实力和和经济发展潜力的指标做为变量。

1.2 系统指标变量设置

系统中变量按性质不同化分为状态变量、速率变量和辅助变量。状态变量是描述系统的积累效应的变量,速率变量是描述系统中积累效应变化快慢的变量,辅助变量是状态变量与速率变量之间信息传递和转换过程的中间变量[12]。表1中列出了3个子系统的指标变量以及它们的单位和变量类型。

1.3 根据反馈关系建立系统流图

系统流图是因果关系图的扩展,是一个多层次、多节点、多回路的网络图,它是系统结构的直观表现,反映了不同变量之间的因果关系和相互影响的关系[12]。在分析系统内各变量之间的反馈关系的基础上,利用系统动力学原理,构建的山东省科技发展系统动力学流图,如图1所示。

1.4 系统方程式设置

通过查阅1978—2008年《山东统计年鉴》、《山东科技统计年鉴》,得到山东省科技和经济发展各项指标的统计数据,采用SPSS、EViews等专业统计软件,利用趋势拟合法建立有关变量之间的定量因果关系。通过不断测试和调整,得到以下主要系统方程式。

(1)科技投入模块

财政科技投入=INTEG(科技投入增加量),

企业R&D投入=INTEG(企业R&D投入增加量),

企业R&D投入增加量=

(-16.2601+0.022579*GDP年增加量)*

(1+企业R&D投入增长率(Time),

对外技术依存度=

技术引进经费/-(技术引进经费+

企业R&D经费支出),

R&D科学家和工程师数=

0.825*科技活动人员数-48 668.3,

科技活动人员数=0.034*财政教育经费+182585。

(2)科技产出模块

高新技术产业产值=-1 773.34+

51.955 4*R&D投入-

0.00211287*R&D科学家和工程师数,

科技贡献率=0.011 011+

LN(R&D投入/10 000)*0.070 218 9,

科技论文数=1.608*科技活动人员数-307030

新产品销售收入=1 206 540+

15.342*企业R&D投入。

(3)经济发展模块

GDP=INTEG(GDP年增加量),

GDP年增加量=GDP*GDP增长率,

固定资产投资额=GDP*固定资产投资率(Time)。

2 系统仿真预测

2.1 系统模型测试

系统动力学模型有效性检验方法可分为直观检验、运行检验、历史检验、灵敏度分析4种方法[12]。我们在分析系统、指标变量选择与绘制流图的过程中,基本上通过了直观检验与运行检验。灵敏度分析,是指通过改变模型中的参数、结构,运行模型,比较模型的输出,从而确定其影响程度,本模型已基本通过灵敏度检验。历史检验,是指选择某一历史时刻为起点,进行仿真,然后用已有历史数据与仿真结果数据进行误差、关联度等检验。这里,我们验证的起始时间是2004年,检验时间是5年,到2008年止。检验变量为GDP,检验结果,如表2所示。

从表2可以看出,GDP仿真值与实际历史值非常接近,相对误差极小,可以满足研究需要。除GDP外,其他指标变量的仿真结果也较为满意,在此就不再赘述。

2.2 系统仿真预测

本文的系统动力学仿真预测中,基本目标是研究山东省在“十二五”期间科技与经济的发展趋势。综合考虑系统各要素的影响作用,预测五年规划期内的以经济增长为核心的各指标的发展变化,考虑到模型的仿真意义以及与现实情况的符合程度,根据科技投入力度对山东省科技经济发展的影响,得出了高、中、低3套指标模拟方案。

2.2.1 科技发展高速增长方案

假定“十二五”期间,山东省的R&D投入处于高速增长阶段,全社会的R&D投入主要包括政府R&D科技投入和企业自筹R&D经费,在该方案中把企业R&D投入增长率和政府科技投入占财政收入的比重作为决策变量对整个科技经济发展系统进行仿真模拟。山东省的R&D投入继续高速增长,2009年R&D投入增长率为35%,“十二五”期间,按照“自主创新、重点跨越、支撑发展、引领未来”的科技发展指导方针,要使R&D投入持续高速增长,就要求企业R&D投入和政府科技投入共同保持高的增长率。设定“十二五”期间企业R&D投入增长率达到35%以上,政府R&D投入占财政收入的2.05%,根据此方案,预测结果,如表3所示。

注:带*的值为2008年数据,各价值量指标按当年价格计算(下同)。

2.2.2 科技发展中速增长方案

假设在未来时期,受诸多因素的制约,科技发展的速度比上述方案略有下降。“十二五”期间,政府科技投入占财政收入的比重为1.8%,企业R&D投入增长率保持在30%左右。科技投入的增长速度下降必将影响科技产出的增长速度,到2015年,R&D总投入为1 584亿元,全社会R&D投入占GDP的比重达到2.29%,科技贡献率为55%,高新技术产业产值为82 278亿元。具体预测结果,如表4所示。

2.2.3 科技发展低速增长方案

假设在未来时期,科技发展由于各方面因素,发展遭遇瓶颈,增速更加缓慢。假定“十二五”期间,政府和企业同时减少科技投入的比重,假设政府科技投入占财政收入比重为1.2%,企业R&D投入平均增长率在“十二五”期间为27%。全社会R&D投入增速减慢,到2015年,全社会R&D投入占GDP的比重为2.22%,高新技术产业产值为79 959亿元,新产品销售收入为22 335亿元,科技贡献率达到54.8%。具体预测结果,如表5所示。

2.3 预测结果分析

通过高、中、低3种指标模拟方案的对比分析可以看出,科技投入的力度与科技产出以及经济发展有着直接的关系。在高速增长方案中,R&D投入增长率远远大于GDP增长率,预测“十二五”结束后,R&D总投入将达到1 626亿元,全社会R&D投入占GDP的比重将达到2.34%。科技产出模块的高新技术产业产值、新产品销售收入、科技论文数量、专利申请与授权数等都得到显著增长;在中速发展方案中,科技投入处于中等水平,科技产出方面相比高速增长方案略有下降,预测“十二五”结束后R&D总投入比第一方案减少42亿元,全社会R&D投入占GDP的比重达到2.21%,相比第一方案下降了0.13个百分点,高新技术产业产值以及新产品销售收入等也有所减少;低速发展方案中,由于受各方面因素的影响,科技发展遭遇瓶颈,发展速度明显缓慢。在该方案中,预测“十二五”结束后,R&D总投入为1 539亿元,相比第一方案减少87亿元,相比第二方案减少45亿元。全社会R&D投入占GDP的比重达到2.03%,这个比重相比第一方案减少了0.31个百分点,相比第二方案减少了0.18个百分点。另外高新技术产业产值以及新产品销售收入等受限于科技发展的速度,相比前两个方案都明显减少。

3 山东省科技发展对策与建议

3.1 深化科技管理体制改革,完善科技创新体系

改革科技管理体制,健全有利于全社会科技进步的科技管理体制。各级政府科技行政部门要转变职能,加强对经济和社会发展各部门特别是各产业部门和企业科技进步工作的指导与支持,建立科技管理的知识产权保护制度,充分保护研发成果,激励科研人员的创新动力。完善以企业为主体,各高等院校、科研院所积极参与、产学研相结合的技术创新体系。

3.2 建立科技投入稳定增长机制,切实保障科技投入力度

针对山东省财政性科技投入强度偏低的现状,政府应该拓宽资金投入的渠道,增加资金投入的数量,建立科技投入稳定增长机制,切实提高财政性科技投入占地方财政支出的比例。政府科技投入的增加对企业起到引导作用,使企业在技术创新中的主体作用更加突出[10,11]。企业要建立符合科技活动规律和特点、推动持续创新的科技投入机制,进一步提高企业在成本中列支的研究开发经费在年销售额中的比例,增加科技经费筹集额,加大科技项目开发力度,在着力解决企业结构调整、产业升级过程中的重大技术瓶颈。

3.3 加强科技人才队伍建设

进一步深化科技人才管理机制,坚持科技人才培养与引进并重,加强政策和资金扶持,建设合理的人才梯队,切实将人才强省战略落在实处。建立科技人才使用和激励机制,对科技人才实施动态管理,确保科技人才在合适的岗位上实现自身价值,深化科技成果奖励制度,建立科学的科技成果评价机制,激发科技人员的持续创新能力。

3.4 深化区域之间科技合作与交流

鼓励省内高等院校、科研单位、企业开展与省外高水平院校、科学院、工程院等的科技交流与合作。加强与国家各专业协会的合作,掌握科技发展前沿动向,重点解决全省科技发展难题。建立与江苏、广州等较高科技发展水平地区的密切合作,吸取科技发展的经验和教训,根据山东省具体情况,因地制宜,切实促进全省科技发展水平。

3.5 建立和完善全面的科技发展绩效考评体系

建立全社会科技发展绩效考评体系,从全社会的科技投入、科技产出与科技促进经济发展三方面制定科学合理、易于操作的评价指标体系,引导企业及全社会增加对科技发展的贡献。完善政府科技发展绩效考评体系,从山东省和17地市两个层次建立政府科技发展年度绩效考评体系,制定奖惩措施,奖励有突出表现的地市,鞭策排名靠后的地市,切实将科技发展绩效考评作为政府工作业绩的重要考评内容。

参考文献

[1]ZULCH G,ROTTINGER S,VOLLSTEDTT.Asimulation approachfor planning and re-assigning of personnel in manufacturing[J].International Journal of Production Economics,2004,90(2):265-277.

[2]AHN H S,RIGHTER R,SHANTHIKUMAR J G.Staffing decisionsfor heterogeneous workers with turnover[J].Mathematical Methodsof Operations Research,2005(62):499-514.

[3]王维,李仕明,李钰.2005—2010年四川省人才需求预测[J].电子科技大学学报:社科版,2005(2):44-47.

[4]董奋义.科技投入与科技产出的关联分析及趋势预测[J].技术经济,2009(7):22-24.

[5]王建玲,刘思峰,邱广华,等.苏州市科技创新人才建设现状及供给预测研究[J].科技进步与对策,2010(12):141-144.

[6]杨月,沈进.多元线性回归分析在人才需求预测中的应用[J].商场现代化,2006(11):33-34.

[7]杨雪,李建华,成宝英.基于人工神经网络的科技需求能力测度与预测问题研究[J].科技管理,2004(12):29-32.

[8]常新宇,代宏坤.四川科技人才需求预测研究[J].四川理工学院学报:社会科学版,2006(2):67-70.

[9]张建勇,赵涛.基于BP神经网络的科技人才需求预测[J].科技管理研究,2009(8):501-502.

[10]李新运,于永信.21世纪初山东省经济增长对科技投入的需求分析[J].科学与管理,1999(6):13-16.

[11]李新运,于绥贞,于永信,等.山东省科技投入需求分析和投入体系构建研究[J].中国软科学,2001(3):34-38.

动力学模拟论文篇7

关键词：金属玻璃,胶体实验模拟,粒子扩散,团簇结构

金属玻璃是近几十年材料领域的研究热点, 与传统晶态合金相比, 金属玻璃的原子排列短程有序、长程无序, 独特的结构特性使其在强度、硬度、冲击断裂性能、耐腐蚀性以及耐磨性等诸多力学性能方面具有显著优势, 在航空航天、汽车工业、化学工业、运动器材等方面具有广泛的应用潜力［1］。由于金属玻璃不具备晶态结构的长程有序性, 所以其主要的结构特征是局域的短程-中程有序结构［2］。金属玻璃局域结构方面的研究是本领域最活跃的基础科学问题之一。借助于实验手段和计算技术的不断发展, 近几年人们对金属玻璃的局域结构逐步有了清晰的认识［3—5］, 但是到目前为止对于金属玻璃短程有序结构的描述仍存在很大分歧［3，4］。

普遍认为金属玻璃的短程-中程有序结构是由金属熔体, 特别是由过冷却液的微结构演化而来, 所以研究金属熔体的结构特征对于进一步深入了解金属玻璃的结构有很重要的意义。对于单原子的金属熔体, 经典的模型研究主要有Bernal等［6］提出的硬球随机密堆模型; Finney［7］、Hoare等［8］、Doye等［9］考虑到原子之间作用势, 在此基础上提出的软球堆积模型。实验方面, Turnbull于1952 年首次得到了Hg的过冷却液滴［10］。之后随着电磁、静电技术的进步, 也逐渐得到了包括Fe、Zr、Ni、Co、Ti等在内的多种金属的过冷却液［11—15］。通过电子衍射和中子散射, 对这些单原子金属熔滴的研究表明, 虽然在液体中普遍存在的二十面体局域有序单元, 但每种元素体系的二十面体单元又有所不同［16—20］, 所以单原子金属熔体的结构仍然存在很多疑点。造成这种情况的主要原因是实验手段的局限性, 无法直接观测到原子排列, 而衍射得到结果具有不确定性和不唯一性。

近几十年来, 胶体悬浮物已被证明是研究原子/分子体系晶体和非晶态材料物理行为的有效模型体系［21—23］。胶体模型系统的突出优势是, 它们可以通过显微镜直接观测单个颗粒的运动、排列规律, 这些颗粒可以作为介观尺度的模型原子, 为原子 ( 分子) 晶体材料的物理过程提供可视化的实验依据。同时由于尺度上的增大, 也使这些物理过程的时间尺度显著增大, 这使实际材料中一些难以观测的非平衡过程能够通过动力学控制给予冻结, 这将为实际材料的物理过程的认识提供直观的证据和启示。

本文通过胶体模拟的实验方法, 调节颗粒浓度, 制备了不同状态的胶体液态试样来模拟金属熔体。针对最简单的单原子液态体系, 从二维角度出发, 结合Voronoi分形技术分析了不同浓度下液态体系的团簇结构, 并从动力学角度探讨了粒子的运动规律和扩散规律。

1 实验

本实验采用聚苯乙烯小球模拟单个原子, 粒径尺寸为2. 9 #m, 球形度好, 单分散性小于5%。首先将聚苯乙烯小球的悬浮液用超纯水多次离心洗涤, 去除杂质备用。然后将小球充分分散在按一定比例配置的超纯水和酒精溶液里。接着将悬浮液滴加到直径为10 mm的试样盒中静置。由于重力作用, 悬浮液中的颗粒会慢慢沉降到试样盒底部的载玻片上, 以载玻片作为基片形成由聚苯乙烯小球组成的单层薄膜。载玻片在使用前依次经过铬酸洗液, 丙酮和酒精浸泡处理, 以获得良好的表面亲水性。控制不同的悬浮液浓度和滴加量, 可以得到不同致密度的薄膜。随着致密度的增加, 可以观察到胶体气态、液态、晶态和玻璃态。将静置的试样放在显微镜下, 由于聚苯乙烯颗粒疏水, 小球不会被基片粘连住, 但是由于重力作用, 小球只能在基片平面内作准二维的布朗运动。按照2. 5 s/张的采样频率进行采集, 可以得到粒子的连续运动图像。

2 结果和讨论

图1 给出了不同浓度下得到的二维胶体体系图像。调节悬浮液的滴加量, 可以得到不同的胶体状态。这里, 用颗粒占基底平面区域的面积分数 ρ 来标识不同的胶体状态。当 ρ = 0. 03 时[图1 ( a) ]得到的是胶体体系的气态, 当 ρ =0. 11 时[图1 ( b) ]得到的是胶体体系的极稀液态; 当 ρ =0. 2 ~0. 4 时[图1 ( c) ~ 图1 ( e) ]对应的是不同浓度的液态; 而当ρ = 0. 56时 ( 图1 ( f) ) 得到的是稠液态。对于实际金属而言, 由固态融化为液态, 密度降低很少, 其液态结构最大的特征就是高密度。所以 ρ = 0. 56 时, 胶体体系的粒子排布和原子体系液态的结构分布有共同的特征, 可以用于模拟胶体玻璃的液态。

在无序体系结构研究中, 为了描述粒子密度的分布, 引入径向分布函数

式 ( 1) 中, 平均数密度N为粒子总数, V为总体积, ρ0是平均粒子数密度, ( …) 表示系综平均, rij表示粒子i和j之间的距离, 式中的距离均以粒子直径的倍数来计算。

对胶体液态体系[图1 ( c) ~ 图1 ( f) ]的粒子位置进行计算, 得到的径向分布函数如图2 所示。可以看到, 对于稀液态, 径向分布函数仅有一个明显峰值, 且第一峰的峰值随着浓度递增。而对于稠液态 ( ρ =0. 56) , 可以看到类似原子体系液态的2 ~3 个峰值。这表明胶体过冷液态体系中存在着短程有序性。

二维情况下, 胶体体系的团簇结构较为简单, 可以根据单个粒子与周围粒子的位形关系, 作出Voronoi分形图, 即以单个粒子为中心, 作出Delau-nay多边形, 以此来分析粒子与周围粒子相互作用的情况［24］。一般来说, 若得到的Voronoi分形图显示该粒子位于周围粒子形成的几边形区域的中心, 则认为周围粒子形成的是一个几边形团簇结构。对液态浓度 ρ =0. 2 ~0. 6 取试样分析, 得到的每一种浓度某一时间的Voronoi分形图如图3 所示。用不同的颜色标识不同的多边形, 即红色-四边形、黄色-五边形、蓝紫色-六边形、绿色-七边形、青色-八边形。值得注意的是, 在研究的胶体体系中, 绝大部分都是以这五种情况出现, 占整个体系的99. 5% 以上, 其他边数的多边形出现概率较低, 所以这里主要统计这几种局域结构的情况。由于边缘效应, 图片四周边界的粒子不计入统计范围。从图3 中可以观察到, 浓度较低时, 如图3 ( a) 和图3 ( b) 中, 单个粒子所占的多边形面积较大, 因而粒子扩散能够活动的区域范围也大, 当浓度增加时, 图3 ( c) 和图3 ( d) 中, 粒子所占的多边形面积明显减小, 粒子能够活动的区域范围也相应缩小。在 ρ =0. 21 和 ρ =0. 30 的体系中, 可以看到有大小为5 ~6 个粒子区域的蓝色六边形区域。对比其光学显微照片, 发现这些表征着有序结构的六边形正对应着粒子聚集而形成的团簇结构。而在浓度较高的 ρ =0. 43 和 ρ =0. 56 的体系中, 黄色和绿色区域形成链状, 蓝色区域分布在链的中间。图中蓝色区域的大小为10 ~20 个粒子, 且在 ρ =0. 56 的体系中尺寸最大。这说明, 随着浓度的增加体系中的结构有序尺度递增。

( (a) ρ=0.21; (b) ρ=0.30; (c) ρ=0.43; (d) ρ=0.56)

对图3 中不同浓度液态体系多边形数目统计的结果如图4 所示。可以看到, 四边形和八边形的数目占比较小。随着浓度 ρ 从0. 30 递增到0. 56, 单元液态体系的五边形占比逐渐减小, 六边形和七边形数目增多。其中, 六边形代替五边形成为最主要的结构单元。对于 ρ = 0. 21 的反常结果, 主要原因之一是由于浓度太低, 粒子数目统计不够造成。

在液态和金属玻璃中, 也普遍存在二十面体团簇结构。二十面体团簇结构具有五次对称单元［25，26］。实际金属中, 团簇结构都有所扭曲, 所以得到的五次对称单元, 并非完美的五边形。在对二维液态情况下粒子位形拓扑结构的研究中, 发现液态的二维局域结构主要有五边形和六边形两种。四边形主要起到连接这两种结构的作用。单一组分的胶体体系和金属原子体系一样, 易于形成晶体［27］。晶态结构一般为f. c. c, 其主要二维平面结构特征就是六边形。一定程度上来说, 五边形数目和六边形数目代表了非晶态和晶态两种机制竞争存在的情况。五边形相对数目的增加表明其非晶态程度的增加。

追踪不同的液态浓度下粒子的运动轨迹。根据粒子运动的轨迹可以得到粒子的均方位移 ( MSD) 和时间t的关系:

式 ( 2) 中, ri ( t) 表示第i个粒子在t时刻的位置, ri ( o) 表示第i个粒子初始时刻的位置, ( …) 表示系综平均。

不同浓度的MSD如图5 所示。从图5 ( a) 中可以看到, 在稀液态状态下 ρ =0. 21 时, 粒子随时间运动的位移明显高于其他的液态情况。而在稠液态状态下 ρ =0. 56, 粒子随时间的均方位移变化要明显低于浓度较稀的情况。这表明, 随着浓度的增加, 粒子运动受到的周围粒子位形带来的阻力越来越大。

浓度的增加表现在显微图像上是粒子周围的间隙越来越小, 粒子越过周围的粒子形成围笼需要的能量越高, 直到粒子数目高到一定程度, 整个体系形成玻璃态, 粒子的运动大多被局域在周围粒子形成的笼子里［28］。

笼子效应将粒子自由扩散分成三个时间区段, 能够从均方位移-时间曲线的lg-lg变换中看出来。在短时间和长时间区段, 图5 ( b) 中两条曲线的斜率均为1, 表明粒子的运动是自由扩散。在中间时间区段, 斜率要小于1, 由于笼子效应, 粒子的运动进入一个交叉区域。从图中可以看到浓度高的体系 ( ρ =0. 56) , 开始能观察到微弱的由于笼子效应形成的中间时间区域。而对于浓度较低的体系 ( ρ≤0. 43) , 并不能看到明显的笼子效应。这是因为高浓度使得粒子处在更小的笼子里, 所以粒子与周围粒子撞击的几率减小, 逃离出笼子的时间增加。定义中间区域出现和结束的时刻为曲线上瞬时斜率偏离1 达到10%的时刻。面积分数为 ρ =0. 56 的单元胶体液态体系其出现结构弛豫初始时刻4 s, 结束时刻6 s, 弛豫持续时间为2 s。

短程扩散主要是由于粒子之间复杂的水动力学相互作用引起的。在短程扩散区段内, 浓度高的体系扩散速率慢, 但比长程扩散速率大, 主要原因是笼子效应会减慢长程扩散的速率。单元体系的扩散速率随着浓度的增加而减慢, 在高浓度 ( ρ = 0. 56) 观测到由于受到周围原子的限制, 形成弱的笼子效应, 粒子的扩散受到阻滞, 其长程扩散系数为0. 10#m2/ s, 短程扩散系数0. 11 #m2/ s, 两者差别不大, 说明笼子效应对粒子的长程扩散影响较小。

3 结论

( 1) 本文从胶体体系实验角度出发, 模拟研究了二维单原子金属液态的微结构; 研究结果表明在液态结构可以直接观察到短程-中程局域有序性的存在。

( 2) 对胶体体系的Voronoi分形研究表明, 在二维体系中, 五边形和六边形这两种局域结构特征竞争存在。随浓度增加, 五边形数目减少, 六边形代替五边形成为主要的特征结构单元。

动力学模拟论文篇8

随着各国环境保护意识的增强和世界能源结构在逐渐变化, 天然气成为最受欢迎的能源之一。换热器中预冷和深冷液化液化天然气重要的一个环节。这是一个涉及能量、热量和质量传递的相变过程以及气液两相间界面的追踪的过程。1981年, Hirts和Nichols[1]提出了VOF (Volume of Fluid) 方法, 并使用该方法对溃坝问题进行数值模拟, 论证方法的可行性。该方法的相界面构造基本思想为运动界面追踪问题的数值模拟起到了开创性的作用;1998年, Boris Halasz[2]从热力学角度, 以能量、动量、质量平衡为基础, 通过分析蒸发式冷凝器内传热传质及流动阻力, 总结出了当时所有类型的蒸发式冷却装置通用的数学模型。1988年, Osher和Sethian[3]提出Level Set (水平集) 方法, 较精确计算相界面曲率及相关的物理量。2015年, Li S[4]等人通过引入热平衡模型, 对亚音速蒸汽注入过冷池直接接触冷凝进行数值模拟, 结果表明在管出口轴向温度随轴向速度降低而升高和压力震荡主要受蒸汽流速、蒸汽冷凝和低温冷却水的静压力影响。

本文的主要目的对天然气在换热器管程中进行深冷液化的问题进行研究, 分析在低温环境下甲烷深冷液化的主要影响因素。本文采用CFX中多相流混合模型, 模拟低温环境下天然气深冷液化的相变传热传质过程, 分析了不同因素对液化的影响及液化后的流型。

1 数值方法

1.1 几何模型

由于在换热器中管程通常为弯管结构, 本文几何模型采用“S”型弯管 (如图1所示) 。模拟计算过程中, 入口通入高温甲烷气体;气体与管道壁面换热, 使管道中甲烷遇冷液化;出口处则为甲烷的气液混合物。

1.2 控制方程

混合物模型[5] (Mixture) 是一种简化的多相流模型, 它用平均速度的概念来模拟多相流中各相具有不同速度的情况。该模型能够求解混合相的能量、动量和连续性方程, 以及各相的体积分数、温度、压力等物理量。

对上述模型, 通过建立动量守恒方程、连续性方程、体积守恒方程以及相间热传递方程。再进行数值模型的稳态求解, 不考虑控制方程在时间上的连续性, 控制方程如下:

(1) 动量守恒方程

式中, ρ和分别表示混合相的平均密度和速度, 表示体积力。uF

(2) 连续性方程

式中, SM描述用户指定的质量源, Γab指从β相到α相的单位体积质量流量

(3) 体积守恒方程

式中, rα表示α相体积比, Np是流体域内所有相的数量。

(4) 相间热传递方程

式中, hα, Tα, λα分别表示α相静焓、温度、热导率, SEα是外部热源, Qα其他相传递到α相的内能。

1.3 网格划分

利用Hypermesh划分模型网格, 流体域采用四面体网格划分, 在壁面添加10层边界层, 第一层厚度0.01mm, 增长率1.2。保持壁面边界层参数不变, 改变网格最小单元尺寸。网格最小单元尺寸:0.5mm、0.7mm、1mm、1.5mm、2mm和2.5mm, 网格数量分别为2605506、1286105、449675、203434、93534和55231。分别对以上5中网格进行模拟计算, 提取管道出口位置温度的平均值与最大值。根据分析可以得出从网格数量449675以后, 网格数量的增加对模拟结果的影响很小。本文选定最小单元尺寸0.7、网格数量1286105作为分析对象。

1.4 边界条件

本节主要针对流体域动力学参数和流体域热力学参数进行分析。由于液化后的液体受重力影响对流动影响较大, 需要考虑重力和浮力的作用;甲烷气相和液相会形成交界面, 因此选用CFX多相流自由液面模型;流体域入口控制参数设置为总压、温度及气态甲烷的体积比;出口设置为静压;所有壁面均为无滑移壁面模型。

流体域内的热力学模型为总能量模型, 相间设置热传递系数1000W/m2k;管程和壳程之间存在热量传递, 为方便计算设定壳程流体温度均为100K, 故管外温度设置为100K, 管的传热系数为20000W/m2k。其入口压强随重力方向而改变。

2 数值结果与讨论

(1) 数值结果的验证。甲烷在换热器管程中深冷液化涉及气液相变及多相流动问题。高温甲烷气体流经壁面, 遇冷液化。同时, 甲烷气体与液化后产生的液体同时在管道中流动。气体液化主要发生在管道表面, 当重力作用可忽略不计时, 将会形成环状的流行 (如图2所示) 。图3展现的为高温甲烷气体从管道入口到管道出口气液百分比, 颜色越深液体比例越高。由图可以看出, 越靠近出口, 甲烷液体越多。与Osher S等实验结果基本符合, 证明了仿真的正确性。

(2) 驱动压力的影响。保持管道壁面的传热系数不变, 设置7组递增驱动压力的计算组, 分别为100Pa、200Pa、300Pa、400Pa、500Pa、600Pa和700Pa。通过仿真模拟计算, 获得管道出口处甲烷的冷凝量。图5为不同驱动压力时, 出口处甲烷冷凝量的变化曲线。由图4可以看出, 在100~300Pa阶段, 冷凝量随着驱动压力的增大成直线增大;300~500Pa阶段, 冷凝量成稳定状态;500~700Pa, 冷凝量又回到直线增长。

(3) 管壁传热系数的影响。保持初始温度和驱动压力不变, 设置9组管道壁面传热系数递增的计算组, 分别为1、2、6、8、10、12、12.5、15、17.5和20k W/m2k。通过模拟计算, 获得管道出口处甲烷的冷凝量。图3-4为不同管道壁面传热系数时, 出口处甲烷冷凝量的变化曲线。由图5可以看出, 随着管道壁面传热系数的增大, 出口处冷凝量呈直线增加。管道壁面传热系数增大, 管壁的热阻减少, 使气体与管壁的传热效率提高, 从而使冷凝的液体增多。

3 结论

本文主要研究驱动压力和管壁传热系数分别对目标参数 (出口处甲烷冷凝量) 的影响。通过简化管道模型, 并采用流体动力学混合物模型结合自由液面模型 (VOF) 对其进行数值模拟, 得到以下结论:

(1) 数值模拟结果与文献中的流型一致性较好, 从而证明了本文所建模型的合理性。

(2) 在优化的范围内, 若要提高甲烷的冷凝量, 可通过提高驱动压力和使用传热系数较高的管道材料。

参考文献

[1]Hirt CW, Nichols BD.Volume of Fluid (Vof) Method for the Dynamics of Free Boundaries.J Comput Phys.1981;39 (1) :201-25.doi:10.1016/0021-9991 (81) 90145-5.

[2]Halasz B.A general mathematical model of evaporative cooling devices.Rev Gen Therm.1998;37 (4) :245-55.doi:10.1016/S0035-3159 (98) 80092-5.

[3]Osher S, Sethian JA.Fronts Propagating with Curvature-Dependent Speed-Algorithms Based on Hamilton-Jacobi Formulations.J Comput Phys.1988;79 (1) :12-49.doi:10.1016/0021-9991 (88) 90002-2.

[4]Li SQ, Wang P, Lu T.Numerical simulation of direct contact condensation of subsonic steam injected in a water pool using VOF method and LES turbulence model.Prog Nucl Energ.2015;78:201-15.doi:10.1016/j.pnucene.2014.10.002.

[5]江帆, 黄鹏.Fluent高级应用.2008.

动力学模拟论文篇9

电沉积技术广泛用于制备各种纳米材料, 如纳米线[1,2]、薄膜以及多孔层材料。然而, 尽管过去的50年中电沉积的机理得到了广泛的研究, 但在电沉积过程中有关形核和生长的基本原理依然存在着争议[3]。因此, 电化学金属沉积过程长期吸引着科研工作者和工程师的关注[4,5]。

电极和溶液界面处的动力学过程对于电沉积是非常重要的。对形核和生长的深入理解有利于优化产品的微观结构和性能[6,7]。关于电沉积形核和生长的理论可以大致分为2类:扩散控制和动力学控制。大多数扩散控制行为的理论模型[3,8,9]都是基于连续介质方程和与形核率、核密度以及核分布有关的先期假设。在连续介质模型中, 吸附过程的细节通常被忽略, 因为岛被认为是半球形, 增加1个原子只是简单的增加球半径。实际上, 处理原子尺度下的薄膜生长, 连续方程模型已经不适用, 应该发展随机方法来描述金属离子还原和晶格形成过程中的近表面化学、反应机制和复杂的分子运输模式。

人们对动力学控制行为的兴趣日益高涨, 它可以在比较低的过电位下形成很小的金属团簇。一些早期的工作中, 很多学者研究了动力学控制电沉积的早期形核和生长[10,11,12]。例如, Stephens等[12]用动力学Monte Carlo数值方法模拟了动力学控制条件下形核和生长, 发现当表面扩散和沉积速率的比值较大时沉积层沿着单原子台阶边缘生长。然而, 大多数上述工作只注重早期的形核和生长而非电镀层的整体性质。本研究旨在深入讨论电沉积薄膜的整体微观结构, 材料的微观结构对它们的性质乃至应用都有重要影响。

本研究扩展了一种原子尺度下的二维横截面模型并用动力学Monte Carlo方法模拟了动力学控制条件下电沉积多晶镍薄膜的生长。二维电沉积[13,14]是一种被广泛研究的非连续、非局域、非平衡典型的生长过程。本研究详细分析了电沉积条件对镀层微观结构的影响, 为电沉积过程中结构演变的微观机制提供了深入的了解。

1 模型与算法

本研究采用动力学Monte Carlo (KMC) 方法对生长在镍 (100) 基底上镍薄膜的横截面 (三角晶格) 进行二维模拟。为简化运算, 基底和沉积金属的几何细节都采用最简单的情况, 即模拟模型是理想体系。考虑与新相物质的生长有关的2个基本过程: (1) 镍离子在阴极还原成吸附原子; (2) 吸附原子在生长表面上的扩散。在一个KMC循环中, 扩散原子只能跳跃到同层或上、下层原子的最近邻位置, 构建一个包含系统所有可能事件的列表, 然后随机选择一个事件执行。选中某个事件的概率正比于它的速率。因为模拟是在较低的应用电势条件下进行的, 由体扩散导致的镍离子消耗微乎其微, 所以假设在所有的表面位置, 在整个模拟时间内, 吸附原子的沉积率都是常数[11]。扩散速率依赖于原子位置和将要跳到位置原子的占有情况。

式中:υ0=2kBT/h是尝试频率, kB是玻尔兹曼常数, h是普朗克常数, T是溶液温度, E是扩散能垒, 其表达式见式 (2) [15]:

式中:Es是单个吸附原子表面扩散的激活能, ED是先沉积的最近邻原子与该原子的相互作用之和, Eboun是原子处于界面时的附加能量。晶格体系中每个点都分配一个整数来描述它所在晶粒的晶体学取向, 相同晶粒中点的该整数值都是相同的, 不同晶粒中点的该整数值不同。采用嵌入原子方法 (EAM) [16]计算出所有可能跳跃构形的能量势垒, 就如文献[17]中所述。

动力学Monte Carlo方法的算法描述如下: (1) 选择一个 (0, 1]之间均匀分布的随机数r; (2) 由满足undefined的整数l选取事件rl, 式中ri是事件i的速率常数; (3) 执行事件rl; (4) 更新所有由于执行上述事件而改变的速率常数rj; (5) 更新模拟时间Δt。执行完步骤 (5) 以后回到步骤 (1) , 如此循环往复。

对于即将发生的任意给定事件, 模拟时间Δt简单地用系统所有事件速率常数总和的倒数来计算, 出于数学完备性考虑, 再乘上因子-ln (ρ) , 其中ρ是0～1之间的随机数。

undefined

在本研究的模拟计算中, 镀层横截面的尺寸为180nm×18nm, 以水平方向为周期性边界条件沿垂直方向生长100层原子。随着沉积率由10nm/s到100nm/s, 温度由280K到360K变化, 比较其计算结果, 以此来研究沉积率和温度对电镀薄膜晶粒尺寸、相对密度和表面粗糙度的影响。考虑到模拟计算的随机性, 本研究中所有数据都是10次以上独立模拟计算的平均值。

2 结果与讨论

图1显示的是分别沉积10、50和100个原子层 (ML) 后多晶镍薄膜的截面模拟图。不同的颜色表示不同的晶粒。从图1中可以初步观察到, 动力学控制条件下镍薄膜生长过程中微观结构的一些基本特征。首先, 在较高的温度或较低的沉积率条件下, 会形成比较大的岛, 然后长成相对比较粗的晶粒。其次, 在较低的温度或较高的沉积率条件下, 薄膜中存在较多的缺陷, 如内部孔洞、杂质。再次, 在较低的温度或较高的沉积率条件下, 薄膜表面更粗糙, 在某些地方甚至出现了悬臂梁状的凸起。此外, 晶粒的形状为典型的柱状晶, 与文献[18]中观察到的一致。以上这些特征与实验结果是相符的, 表明KMC模型能够比较好地预测薄膜微观结构的演化。

薄膜的相对密度定义见式 (4) :

undefined

式中:N1为系统中所有的晶格点数, N2为内部孔洞所占住的晶格点数。相对密度越大意味着内部孔洞越少。晶粒尺寸用晶粒所包含的晶格点数描述。图2显示了相对密度和平均晶粒尺寸的变化。

由图2可以看出, 薄膜的平均晶粒尺寸和相对密度随着温度的升高或者沉积率的降低而增大。这是因为温度升高或沉积率降低时, 生长在表面上的吸附原子的扩散能力增强了。一方面, 在较高的温度下, 更多的吸附原子被激活并扩散到一个更稳定的低能量位置, 从而有利于晶粒的长大和内部孔洞被填充。相反地, 在较低的温度下, 吸附原子粘附在它沉积的位置, 很难移动, 从而形成很多很小的核和更多的孔洞。如同气象沉积多晶膜, 随着基底温度的升高, 观察到增强的表面移动性使晶粒横向尺寸增大[19]。另一方面, 在较低的沉积率条件下, 吸附原子在下一个沉积原子到来之前能扩散一段更长的距离以找到一个具有更多最近邻原子的位置。沉积率越高, 形成薄膜的内部孔洞越多、晶粒尺寸越小。因此, 在动力学控制条件下, 当温度较高或沉积率较低时可以获得较致密的微观结构。

表面粗糙度是薄膜高度围绕平均高度波动值的均方根:

undefined

式中:hi是第i个格点的高度, N是表面格点的总数目。图3显示了表面粗糙度在不同条件下随沉积原子层数的变化。在沉积100个原子层后, 对应表面粗糙度的变化如图4所示。从图3和图4可以看出, 表面粗糙度随着时间的延长而增加, 随着沉积率的降低而减小, 随着温度的升高而减小。界面上的岛状生长模式使生长表面粗糙, 层状生长模式则使表面平滑, 如果提高表面吸附原子的扩散能力使它们能够克服能量势垒而由团簇顶端跳跃到生长层面上, 就能减小生长表面的粗糙度。根据前文的分析可知, 可以通过升高温度或减小沉积率来提高吸附原子的扩散能力, 从而减小薄膜的表面粗糙度。因此, 金属吸附原子的移动性强烈影响着膜的微观结构和形貌[19]。上述结果与电沉积铜表面粗糙度的变化类似[20]。

一般地, 在比较短的生长时间内, 表面粗糙度满足Family等[21]给出的动力学标度特性:W∝tβ, 其中:t是沉积时间, β是生长因子。生长因子决定着表面生长标度行为的类别特征。相同β值的不同系统属于同一种类别, 具有相似方式的行为。

图5为在300K、70nm/s时生长因子的拟合结果, td为沉积一层原子所花费的时间。可见, 本研究所述模型的生长因子为β=0.46±0.04。图6列举了由各种模型得到的β值:1#和2#是MBE模型的模拟值[22,23], 3#是本研究的结果, 4#是扩散控制电沉积的模拟数值, 5#和6#分别是脉冲电沉积和直流电沉积的实验结果[24]。

本课题组计算的生长因子接近于分子束外延 (MBE) 模型的生长因子, 但比扩散控制电沉积模型的生长因子小。这个结果说明电沉积薄膜的表面比MBE膜表面粗糙, 动力学控制电沉积膜的表面比扩散控制电沉积膜表面平滑。扩散控制条件下的复杂形貌是由于不同传输机制的相互影响, 例如阳离子扩散、电迁移、对流、表面扩散等。然而, 在低应用电势的动力学控制下, 表面扩散是薄膜生长过程中的控制步骤, 因此它的生长因子接近于MBE模型的生长因子。另外图6显示, 大多数计算值都比实验值小, 这是由于边缘势垒的存在。当吸附原子从团簇顶端向下通过团簇边缘扩散到生长层上时, 它必须克服一个边缘势垒 (Schwoebel势) 。实际中的边缘势垒可能比较大, 吸附原子扩散到生长层是很困难的, 甚至在有些地方只存在向上扩散[25]。因此, 实际上电沉积薄膜大多数情况下都是岛状生长模式, 具有较大的生长因子。

3 结论

本研究用动力学Monte Carlo方法模拟了动力学控制条件下电沉积多晶镍薄膜的生长情况。重点研究了溶液温度和沉积率对薄膜微观结构的影响。吸附原子的表面扩散能力严重地影响薄膜的微观结构, 而溶液温度和沉积率的改变可以使吸附原子的扩散能力改变, 可以通过控制沉积条件来获得理想的微观结构以优化薄膜的性能。模拟结果显示, 薄膜的晶粒尺寸和相对密度随着温度的升高或者沉积率的降低而增大;恰恰相反, 在相同的情况下, 表面粗糙度却减小。对模拟数据进行拟合, 得到生长因子β=0.46±0.04。因为表面扩散为控制步骤, 所以动力学控制电沉积薄膜的生长因子接近于MBE模型的生长因子。因为实际的边缘势垒比较大, 致使电沉积薄膜大多采用岛状生长模式, 所以生长因子较模拟值大。

摘要：采用动力学MonteCarlo方法模拟动力学控制条件下电沉积多晶镍薄膜的生长。研究了表面扩散对电镀薄膜微观结构的影响。模拟结果显示, 溶液温度和沉积速率是2个非常重要的参数, 它们能改变吸附原子的扩散能力, 因而对薄膜的晶粒尺寸、相对密度、表面粗糙度等微观结构产生显著影响。此外, 还得到该模式下薄膜的生长因子, 大约为0.46。

动力学模拟论文篇10

逆流循环式糙米加湿调质是指糙米自下而上由螺旋搅龙输送到加湿圆盘,由加湿细雾喷头向旋转圆盘上的糙米喷雾通风,称之为逆流通风加湿;着湿后的糙米均匀散落在均质仓内,完成水分由外向内的水分渗透。若1次加湿不能满足最优磨米水分的需求,可以打开活动套筒,对糙米进行2次加湿。如此循环,直到满足要求为止,称为循环加湿[1]。在逆流循环式糙米加湿工艺中,加湿细雾和糙米物料以相反方向运动,因此逆流循环式加湿仓具有的特点:一是逆流循环式糙米加湿仓可将储存1年的稻谷或糙米(含水量在11%～13%)均匀加湿调质到15%～17%;二是采用逆流循环式加湿可实现糙米加湿可控、均匀与高效的设计目标;三是加湿工艺设计结构简单,操作方便,磨米加工碎米率低。逆流循环式糙米加湿仓的缺点是加湿量既不能过多也不能过少,润糙时间不能过长,否则起不到降低能耗和提高大米质量的目的。

日本的糙米调质技术研究起步较早。从1955年开始,日本经济高速增长,20世纪50年代末期,日本现代大型精米加工厂急剧增加,民间的技术研究开发投资高涨,主要以应用技术为主,并在20世纪70-80年代达到全面的发展[2,3]。近年来,韩国在糙米调质方面也做了大量研究。韩国国立庆尚大学宋大斌等人研制的连续式调质机,其调质实验结果确认整米率增加2.2%。逆流糙米加湿机具有效率高和节能等优点,但存在需要控制加湿量、润糙时间和热风温度等问题,在我国应用较少。为了克服逆流糙米加湿机加湿不均的缺点,目前逆流式糙米加湿机多采用循环式,且保证加湿均匀可控。逆流循环式糙米加湿机具有加湿均匀、能耗低、润糙与加湿相结合以及能保证大米加工整精米率等优点[4,5],因此在日本和韩国被广泛采用。我国的逆流循环式糙米加湿机应用较少,东北农业大学工程学院的贾富国教授研制了逆流循环式糙米加湿调质机的试验样机。

虽然我国在引进吸收的基础上开发了一些逆流式谷物加湿机,但缺乏深入研究逆流谷物通风加湿机的机理。其原因在于逆流循环式糙米加湿机存在两点边值问题,是所有谷物调质加湿机中较为复杂的一种机型。20世纪70年代以来,美国密执安州立大学农业工程系F.W.Bakker-Arkema教授,根据传热传质的基本理论推导了谷物干燥的理论模型。此理论模型由4个偏微分方程组成,通过此方程利用数值解法可以求出在一定的时间后薄层谷物的水分、温度、热空气的湿度和热空气的温度。目前,此模型已经用于玉米、水稻和小麦等干燥过程的模拟。此模型通过可逆处理后,可用于谷物加湿过程模拟,其优点是通用性好,可以用于顺流、逆流和横流加湿机,而且也考虑了湿空气温度和谷物温度的差别。但此模型的确定计算比较复杂,需要使用大型的计算机[6,7]。为了提高我国逆流循环式糙米加湿机的设计水平,本文运用基于加湿理论的数学模拟方法对逆流循环式糙米加湿机进行模拟研究,并开发了模拟软件。

1 糙米湿化动力学模型建立

薄层湿化方程是谷物深床湿化模拟的基础,是在深床湿化模拟的基础上发展起来的。薄层湿化方程就是描述加湿量、风量、热风相对湿度和谷物初始含水率等参数对湿化速率的影响模型。在进行深床湿化模拟分析时,通常是把深床谷物看作是由许多薄层组成,逐层进行计算,才能得到最终结果。因此,薄层湿化方程便成为谷物湿化模拟的关键。不同粮食具有不同的湿化特性,因此薄层湿化方程也各不相同。

目前,薄层湿化模型可以分为经验模型、半经验模型和理论模型3种不同的类型[8,9,10,11]。经验模型是直接通过试验建立的谷物水分与湿化时间关系的模型。半经验模型假定湿化速率正比于谷物实际水分与平衡水分的差值,因而得出MR=(M-Me)/(Mo-Me)=e-kt的关系式,也是谷物湿化中最常用的一种模型形式。随后的研究发现,利用一个系数k来描述各参数与湿化速率的关系存在一定的缺点,因而出现了MR=e-ktN以及MR=ae-kt等形式的半经验模型。理论模型是基于谷粒内部的液相扩散理论而建立起来的薄层湿化方程,是求解扩散方程而得出的理论模型。它除去假定温度平衡外,还假定谷物表面与周围空气存在水分平衡。由于理论方程比较复杂,扩散系数又受许多因素的影响,不易确定,因此大多数学者通常采用半经验模型进行湿化过程的模拟分析。

典型的薄层湿化总体方程描述为

$Μ R = \frac{Μ_{(t)} - Μ_{e}}{Μ_{0} - Μ_{e}}$ (1)

由式(1)可知,水分比MR为湿化时间t的函数。式中,M(t),Me,Mo分别表示物料任一时刻含水率、平衡含水率和初始含水率(干基)。该方法试图对湿化现象做统一描述,把所有传递性质概括为简单的函数,而不涉及湿化控制机理。

1.1 糙米薄层湿化模型的建立

谷物湿化模型通常是根据试验所得数据,选用适当数学方程,用曲线拟合法求解建立。常见的薄层湿化模型主要有指数模型和Page模型两种。许多学者经试验验证指出Page模型更适用于玉米和水稻[12,13,14]。本试验研究所得的数据也与Page方程基本吻合,因此选用Page方程来建立糙米的薄层湿化模型。Page方程的数学表达式为

$Μ R = \frac{Μ_{(t)} - Μ_{e}}{Μ_{0} - Μ_{e}} = \exp (- k t^{Ν})$ (2)

式中 MR—水分比;

M(t)—某一时刻含水率(kg/kg,db);

M0—初始含水率(kg/kg,db);

Me—平衡含水率(kg/kg,db);

T—湿化时间;

k,N—待定的常数。

利用试验数据,通过最小二乘曲线拟合法,应用数学软件Matlab对试验数据进行处理,求得

K=0.017 5exp(0.007 6Ta)

N=0.258 9+0.057v

式中 Ta—湿空气的温度(℃);

v—湿空气的速度(m/s)。

由此得出糙米的薄层湿化方程为

$Μ R = \frac{Μ_{(t)} - Μ_{e}}{Μ_{0} - Μ_{e}} = \exp (- k t^{Ν})$ (3)

K=0.017 5exp(0.007 6Ta)

N=0.258 9+0.057v

Me是由式(3)求出,即

Me=A-B×ln[-(T+C)ln(RH)]

式中 RH—空气相对湿度,小数;

T—糙米温度(℃);

Me—平衡水分(干基水分),小数;

其中,A=0.293 94,B=0.004 601 5,C=35.703。

1.2 糙米薄层湿化模型的验证

用风速为0.228m/s、风温为25℃进行糙米薄层湿化试验,并将试验所得结果与薄层湿化方程计算值进行比较,如表1所示。

由表1可见,试验值与计算值最大相差为2.89%。这表明,所建立的薄层湿化方程精度可满足计算分析的需要。试验研究表明,糙米水分随着时间的变化,明显地受热风温度、加湿量以及风速的影响。根据试验结果建立的薄层湿化方程,能够较好地描述糙米薄层湿化过程。

2 逆流谷物加湿机模拟软件开发

逆流谷物加湿模拟软件在功能上要求实现:如果输入谷物流速和其它参数,程序运行后可以得出谷物终水分以及其它加湿机性能参数;如果输入谷物终水分和其它参数,程序运行后可以算出谷物流速以及其它加湿机性能参数;输出结果既能够用文本形式表示,又能够用图形表示;软件具有较强的容错能力和帮助功能。

一个实用的商品化软件必须做到内容和形式的完美结合。具体到一个软件而言,不但要求其功能强大,而且要求具有良好的用户界面、图形功能和帮助功能。从软件工程学的观点来看,友好的用户界面成为评价软件好坏的一个重要指标。有人统计,在一个商品化软件中,界面的编程量占总的编程量的60%以上[15]。由此可见,软件界面的重要性。在软件开发前,笔者经过反复比较,决定采用基于 WINDOWS 平台的可视化编程语言Visual Basic和 Visual C++作为开发工具。前者在开发界面方面具有很大的优势,但由于是解释性语言,因而运行速度较慢。采用Visual Basic 编写软件中的计算、绘图部分运行速度较慢。对于一个具有较大计算量和反复迭代的程序来说,运行速度慢的缺点显得尤为突出。 Visual C++ 属于编译性语言,运行速度快,采用两种语言混合编程,优势互补,发挥最佳效益。具体做法是:对于程序中的计算、绘图部分要求计算速度快,采用Visual C++编写;对于用户界面和其它功能要件使用Visual Basic编写;两者的结合是通过动态连接库作为媒体的。在Visual C++编译环境下,将逆流模拟程序的计算部分编译成动态连接库(DLL),存放在指定的动态连接库中。在用Visual Basic编写的程序中,通过声明要使用的动态连接库后,就可以像调用Visual Basic本身的函数那样调用DLL。编写DLL的另一个好处是实现资源共享,其它软件业可以象调用API函数那样调用DLL。如果读者想了解如何编写自己的DLL,Visual Basic和Visual C++的有关知识,请参阅文献[16]。

在逆流模拟软件中,体现逆流谷物加湿调质机谷物和空气的流动特点(即谷物自底向上流动,湿空气自顶向下流动)显得十分重要。只有在程序中真实反映逆流谷物加湿调质机的流动特点,模拟结果才能与实际相符合。下面的C语言代码是从逆流谷物加湿模拟软件中摘录出来的,其功能就是仿真逆流谷物加湿机中的谷物和空气流动。谷物向上流动,用程序表示为

for(i=1;i<=10;i++)

{j=i+1;

wa[i]=wb[j];

Temp[i]=Temp[j];

}

其中,wa[i]和wb[j]分别表示第i层谷物和第j层谷物的水分含量;Temp[i]表示第i层的谷温。

空气向下运动,用程序表示为: hum[j]=hum[i]。其含义将第j层的空气湿含量赋予第i层,即将上一层的空气湿含量赋予它下一层的空气。这就表示空气向下流动。

模拟空气和谷物的流动以及判断谷物和空气系统何时处于稳定状态,可用下面一段程序完整表达出来。为了便于说明问题,在每行程序的开头都标上序号。下面的一段程序是逆流谷物加湿机模拟的核心部分,也是逆流谷物加湿机模拟不同于其他模拟程序的地方。

在逆流谷物加湿机模拟中,加湿过程是由不稳定状态逐渐过渡到稳定状态。模拟程序中就是通过反复迭代来模拟这一过程。采用的方法是将谷床某一空间位置相邻两次迭代计算出的谷物水分差作为判断标准。如果水分差在规定的误差范围内,说明逆流谷物加湿过程达到稳定状态。程序的第5行将上一次迭代的谷物水分值保存下来,在程序的第13行将本次迭代的谷物水分值同上次迭代的水分值进行比较,求出差值。如果水分差小于规定的值,表明加湿过程达到稳定态,模拟过程结束。

3 结束语

1)在本试验条件下,糙米的薄层湿化模型可以用 $Μ R = \frac{Μ_{(t)} - Μ_{e}}{Μ_{0} - Μ_{e}} = \exp (- k t^{Ν})$ 描述。