单原子Lennard-Jones体黏弹性弛豫时间

全文HTML

--> --> -->

2.研究对象及方法

2.1.对象及方法

-->

2.1.对象及方法

为方便对计算结果进行更好的分析, 本文中的物理量均采用“约化单位”^[15], 用上标“*”表示, 如${T}^{*}=Tk/\varepsilon ,\; {\rho }^{*}=\rho {\sigma }^{3}$

, $ {t}^{*}=t\sqrt{\varepsilon /m}/\sigma $

, $ {r}^{*}=r/\sigma $

, $ {p}^{*}=p{\sigma }^{3}/\varepsilon $

, $ {\eta }^{*}=\eta {\sigma }^{3}/\sqrt{m\varepsilon } $

. 其中ε和σ为L-J势参数, 分别表示势能为0时的原子距离、势能阱深度; m为原子质量; T为真实单位下的温度, 单位为K; ρ表示密度, 单位为mol/L; p为压强, 单位为MPa; $ \eta $

为黏度, 单位为μPa·s.
本文建立的单原子L-J体系统采用周期性边界条件^[16], 计算域尺寸(模拟盒子)为16 × 16 × 16单位体积, 并按照相应的密度设置盒子内的粒子数量, 如图1所示.

图 1 L-J体计算系统
Figure1. L-J fluid simulation system.

L-J体系统内分子间仅受范德瓦耳斯力的作用, 其具体势能函数^[17]为

$ {E}_{ij}^{{\rm{L}}\text{-}{\rm{J}}}=4{\varepsilon }_{ij}\bigg[{\left(\frac{{\sigma }_{ij}}{{\gamma }_{ij}}\right)}^{12}-{\left(\frac{{\sigma }_{ij}}{{\gamma }_{ij}}\right)}^{6}\bigg], $

式中r_ij为原子i与j之间的距离, ε_ij为势能阱的深度, σ_ij为两体互相作用的势能为零时的距离. 其中

$ {\varepsilon }_{ij}={\left({\varepsilon }_{ii}{\varepsilon }_{jj}\right)}^{1/2}, $

$ {\sigma }_{ij}= \left({\sigma }_{ii}+{\sigma }_{jj}\right)/2. $

范德瓦耳斯力作用截距半径设置为5.0^[18].
系统中原子的运动受牛顿第二定律及非哈密动量方程控制, 采用正则系统^[19]条件, 时间步长设为0.003. 模拟总时间步数为10⁸步, 系统经过5 × 10⁷个时间步长后达到平衡, 本文选取平衡后系统的数据进行分析, 并统计了系统内所有粒子在最后500个时间步长内的运动轨迹. 本文应用LAMMPS^[20]软件进行建模和具体计算.
2

2.2.数据处理方法

-->

2.2.数据处理方法

本文采用Green-Kubo公式^[7,8,21]计算黏度:

$ \eta =\frac{V}{{k}_{\rm{B}}T}\int_0^\infty { \left\langle {{\tau }_{\alpha \beta }\left({t}_{0}\right){\tau }_{\alpha \beta }\left(t\right)} \right\rangle {\rm{d}}t,} $

式中V表示体积; k_B为玻尔兹曼常数; T为热力学温度; t表示时间; τ_αβ表示剪切应力张量的非对角元素, 其具体表达式为

$ {\tau }_{\alpha \beta }=-\frac{1}{V}\sum _{i=1}^{N}m{v}_{i,\alpha }{v}_{i,\beta }-\frac{1}{V}\sum_{i=1}^{N-1}\sum _{j=i+1}^{N}{r}_{ij,\alpha }{F}_{ij,\beta }, $

其中$ {r_{ij, \alpha }}$

表示i和j粒子在α方向的位移分量, $ {F_{ij, \beta }}$

表示i粒子作用在j粒子上在β方向作用力的分量, N为粒子数, m是粒子质量.
黏弹性动态模量G(t)可通过应力自相关函数公式进行计算^[22]:

$ G\left(t\right)=\frac{V}{{k}_{\rm{B}}T}\left\langle {{\tau }_{\alpha \beta }\left({t}_{0}\right){\tau }_{\alpha \beta }\left(t\right)} \right\rangle . $

运用傅里叶变换公式将(6)式转换为复数型频域剪切模量函数G(ω):

$ G\left(\omega \right)={\rm{i}}\omega {\int }_{0}^{\infty }{\rm{e}}^{{\rm{i}}\omega t}G\left(t\right){\rm{d}}t, $

其中无穷大频率剪切模量^[6]定义为

$ {G}_{\infty }\equiv G\left(\omega \to \infty \right)=G\left(t\to 0\right). $

(7)式中复数模量的实部为储能模量((9)式)^[23], 其反映材料形变时的回弹能力^[24], 即弹性:

$ {G'}\left(\omega \right)=\omega {\int }_{0}^{\infty }G\left(t\right){\rm{s}}{\rm{i}}{\rm{n}}\left(\omega t\right){\rm{d}}t. $

(7)式中的虚部实质为损耗模量((10)式)^[25], 反映材料形变时内耗程度, 即黏性:

$ {G''}\left(\omega \right)={\rm{i}}\omega {\int }_{0}^{\infty }G\left(t\right){\rm{c}}{\rm{o}}{\rm{s}}\left(\omega t\right){\rm{d}}t. $

本文应用Einstein关系式^[26]计算得到系统的扩散率D为

$ D=\mathop {{\rm{lim}}}\limits_{t \to \infty } \frac{1}{6t}\left\langle {{\left|{{r}}_{i}\left(t\right)-{{r}}_{i}\left(0\right)\right|}^{2}} \right\rangle, $

其中t表示时间; r表示原子的位移; 〈·〉表示组内所有原子的平均.
2

2.3.模型验证

-->

2.3.模型验证

为保证模拟的准确性, 本文主要从物质的黏度、系统原子的排布规律两方面对所建立的模型和计算方法进行验证.
本文模拟了平衡状态情况下温度范围在$ {T}^{*}= 0.6—1.5$

, $ {\rho }^{*}= 0.9$

, $ \varepsilon = 1$

, $ \sigma = 1$

条件下的单原子L-J体系统. 不同温度时单原子L-J体的黏度模拟值与Gallie?o等^[27]的文献值对比示于图2. 由图2可知, 随着温度的升高, 单原子L-J体的黏度不断下降; 本文所计算的不同温度下单原子L-J体的黏度与文献值一致, 相对误差在5%之内, 本文的模拟方法具有一定的可靠性和准确性.

图 2 不同温度下L-J体的黏度验证
Figure2. Verification of viscosity of L-J fluid at different temperatures.

为验证模型在微观结构上是否准确, 本文从L-J体的原子排布出发, 对系统的径向分布函数进行验证. 计算了多种物性参数情况下单原子L-J体的径向分布函数, 其中两种不同温度下的单原子L-J体的径向分布函数与Morsali等^[28]经验公式的对比结果示于图3. 由图3可以看出, 本文中的模拟结果与相同温度密度条件下的文献数据^[28]十分符合, 且反映出径向分布函数曲线随温度变化的共同趋势: 在密度为1的情况下, 温度为0.9时(图3(a)), 径向分布函数曲线拥有多个波峰波谷, 数值波动明显, 在$ r^* = 1.0 $

, 2.0, 2.8, 3.7, 4.8依次出现峰值逐渐降低的波峰, 原子排布不规则, 属于黏弹性显著的固液共存态; 当温度升高至4时(图3(b)), 径向分布函数曲线发生明显变化, 曲线平滑, 数值波动减小, 波峰波谷数量显著减少, L-J体处于液态.

图 3 L-J体的径向分布函数验证　(a) ${\rho }^{*}=1,\; {T}^{*}=0.9$

; (b) ${\rho }^{*}=1, \;{T}^{*}=4.0$

Figure3. Verification of the radial distribution function of L-J fluid: (a) ${\rho }^{*}=1, \;{T}^{*}=0.9$

; (b) ${\rho }^{*}=1, \;{T}^{*}=4.0$

.

为进一步验证本文所观测到的原子排布现象, 模拟了$ {\rho }^{*}= 1$

, $ \varepsilon = 1$

, $ \sigma = 1$

, $ {T}^{*}= 0.85—5$

条件下单原子L-J体系统, 其径向分布函数曲线计算结果示于图4. 由图4可以明显看出, 随着温度的升高, 单原子L-J体的径向分布函数曲线逐渐平滑, 数值波动明显减小, 波峰波谷数量显著减少, 单原子L-J体逐渐由液固共存态转变为液态. 上述观察到的现象也与Heyes等^[24]对L-J体相态变化的温度密度范围划分结果一致, 故进一步验证了本文所选方法能较为精确地反映单原子L-J体内原子分布及排布情况, 能有效计算本文所需的相关参数内容.

图 4 ${\rho }^{*}=1, \;\varepsilon =1,\; \sigma =1$

时不同温度下L-J体的径向分布函数
Figure4. Radial distribution function of L-J fluid at ${\rho }^{*}=1, \varepsilon =1,\; \sigma =1$

with different temperatures.

3.模拟结果与讨论

3.1.各类弛豫时间讨论

-->

3.1.各类弛豫时间讨论

为开展单原子L-J体黏弹性时间的研究, 计算了黏弹性特征弛豫时间、Maxwell弛豫时间以及原子连通弛豫时间, 并对各类黏弹性弛豫时间进行了分析.
3

3.1.1.黏弹性特征弛豫时间计算

-->

3.1.1.黏弹性特征弛豫时间计算

通过应用(9)和(10)式计算了$ {T}^{*}= 0.85—5$

, ${\rho }^{*}\!=\! 0.85—1$

, $ \varepsilon\! =\! 0.97—1$

, $ \sigma\! =\! 0.8—1.3$

条件下共22组单原子L-J体的储能模量和损耗模量随频率的变化. 图5为${T}^{*}\!=\! 0.9—1.5$

, ${\rho }^{*} \!=\! 1$

, $ \varepsilon = 1$

, $\sigma = 1$

条件下储能模量$ G^{\prime *} $

和损耗模量$ G^{\prime \prime *} $

的计算结果. 由图5(a)可以看出, 当频率较小(0 < ω^* < 10)时, 储能模量小于损耗模量, 单原子L-J体主要表现为黏性特征, 储能模量和损耗模量均随频率的上升而增大, 单原子L-J体的弹性和黏性均增强; 当频率处于中间频率(10 < ω^* < 20)时, 储能模量和损耗模量接近, 单原子L-J体主要表现出黏弹性特征, 储能模量随着频率的上升逐渐接近并超过损耗模量; 当频率较大时(ω^* > 20), 储能模量大于损耗模量, 单原子L-J体主要表现为弹性特征, 储能模量随着频率的上升而不断增大并趋于稳定, 单原子L-J体的弹性增强, 损耗模量随着频率的上升而不断减小并趋于稳定, 单原子L-J体的黏性减弱. 图5(b)所示为$ {\tau }_{\rm{MD}}^{*} $

的计算方法^[4], 物质的黏弹性特征弛豫时间对应于储能及损耗模量相等时的频率$ {\omega }_{{\rm{c}}{\rm{r}}{\rm{o}}{\rm{s}}{\rm{s}}}^{*} $

, 即:

图 5 储能模量和损耗模量曲线　(a)不同温度条件写的储能模量和损耗模量曲线; (b)粘弹性特征弛豫时间$ {\tau }_{\rm{MD}}^{*} $

Figure5. Storage modulus and loss modulus: (a) Storage modulus and loss modulus at different temperatures; (b) viscoelastic characteristic relaxation time $ {\tau }_{\rm{MD}}^{*} $

$ {\tau }_{\rm{MD}}^{*}=1/{\omega }_{{\rm{c}}{\rm{r}}{\rm{o}}{\rm{s}}{\rm{s}}}^{*}. $

应用(12)式, 在储能模量和损耗模量计算结果的基础上得到单原子L-J体的黏弹性特征弛豫时间$ {\tau }_{\rm{MD}}^{*} $

, 结果示于图6. 由图6可以看出: 当$ {\rho }^{*}, \varepsilon ,\;\sigma $

保持不变时(图6(a)), 随着$ {T}^{*} $

的上升, $ {\tau }_{\rm{MD}}^{*} $

不断减小, 但减小的速度不断下降; 当$ {T}^{*} $

, $ \varepsilon ,\; \sigma $

保持不变时(图6(b)), 随着$ {\rho }^{*} $

的增大, $ {\tau }_{\rm{MD}}^{*} $

不断增大; 当$ {T}^{*} $

, $ {\rho }^{*} $

保持不变时(图6(c)和图6(d)), 势参数$ \sigma $

或$ \varepsilon $

的增大, 均会导致$ {\tau }_{\rm{MD}}^{*} $

的增大.

图 6 不同温度、密度以及L-J势参数下$ {\tau }_{\rm{Maxwell}}^{*} $

, $ {\tau }_{\rm{LC}}^{*} $

, $ {\tau }_{\rm{MD}}^{*} $

, $ {\tau }_{{\rm{m}}{\rm{o}}{\rm{d}}{\rm{e}}{\rm{l}}}^{*} $

　(a)不同温度条件; (b)不同密度条件; (c)不同势能为0时的原子距离条件; (d)不同势能阱深度条件
Figure6. $ {\tau }_{\rm{Maxwell}}^{*} $

, $ {\tau }_{\rm{LC}}^{*} $

, $ {\tau }_{\rm{MD}}^{*} $

, $ {\tau }_{{\rm{m}}{\rm{o}}{\rm{d}}{\rm{e}}{\rm{l}}}^{*} $

under different T, $ \rho $

, ε and σ: (a) Different temperatures; (b) different densities; (c) different distances between atoms when potential is 0; (d) different potential well depths.

3

3.1.2.Maxwell弛豫时间计算

-->

3.1.2.Maxwell弛豫时间计算

通过应用(4)式和(6)式计算了$ {T}^{*}= 0.85—5$

, $ {\rho }^{*}= 0.85—1$

, $ \varepsilon = 0.97—1$

, $ \sigma = 0.8—1.3$

条件下单原子L-J体的黏性$ {\eta }^{*} $

和无穷大频率剪切模量$ {G}_{\infty }^{*} $

, 并计算各条件下的Maxwell弛豫时间. Maxwell弛豫时间的计算公式^[6]为

$ {\tau }_{\rm{Maxwell}}^{*}={\eta }^{*}/{G}_{\infty }^{*}. $

将计算所得到的各条件下单原子L-J体的Maxwell弛豫时间$ {\tau }_{\rm{Maxwell}}^{*} $

与黏弹性特征弛豫时间$ {\tau }_{\rm{MD}}^{*} $

进行比较, 结果示于图6. 由图6可以看出, 在不同的温度、密度和势参数条件下, $ {\tau }_{\rm{Maxwell}}^{*} $

与$ {\tau }_{\rm{MD}}^{*} $

的变化趋势一致; 当${\rho }^{*},\; \varepsilon ,\; \sigma$

保持不变时(图6(a)), 随着$ {T}^{*} $

的下降, $ {\tau }_{\rm{Maxwell}}^{*} $

与$ {\tau }_{\rm{MD}}^{*} $

的差异逐渐增大, 这说明应用Maxwell弛豫时间表征黏弹性特征弛豫时间存在一定的局限性; 当$ {T}^{*} $

, $ \varepsilon , \;\sigma $

保持不变时(图6(b)), 随着$ {\rho }^{*} $

的增大, $ {\tau }_{\rm{Maxwell}}^{*} $

与$ {\tau }_{\rm{MD}}^{*} $

均逐渐增大, 但$ {\tau }_{\rm{Maxwell}}^{*} $

增大的速率更快; 当$ {T}^{*} $

, ${\rho }^{*},\; \varepsilon$

保持不变时(图6(c)), 随着$ \sigma $

的增大, $ {\tau }_{\rm{Maxwell}}^{*} $

不断增大, 但与$ {\tau }_{\rm{MD}}^{*} $

的差异越来越明显; 当$ {T}^{*} $

, ${\rho }^{*},\; \sigma$

保持不变时(图6(d)), 随着$ \varepsilon $

的增大, $ {\tau }_{\rm{Maxwell}}^{*} $

与$ {\tau }_{\rm{MD}}^{*} $

的差异逐渐增大.
上述计算结果表明, 当温度较低($ {T}^{*} < 1.5 $

)时, $ {\tau }_{\rm{Maxwell}}^{*} $

与$ {\tau }_{\rm{MD}}^{*} $

差异明显, 应用Maxwell弛豫时间表征单原子L-J体的黏弹性特征弛豫时间存在一定的局限性.
3

3.1.3.原子连通弛豫时间计算

-->

3.1.3.原子连通弛豫时间计算

原子连通数N_c是指在原子近邻(径向分布函数的第一处谷值位置)范围内的原子数, 原子连通弛豫时间$ {\tau }_{\rm{LC}} $

^[10]是指系统内所有原子中发生一次扩散或者汇聚的时间. 若$ {t}_{0} $

时刻原子连通数为$ {N}_{{\rm{c}}}\left({t}_{0}\right) $

, $ {t}_{0}+{\tau }_{\rm{LC}} $

时刻总原子连通数为$ {N}_{{\rm{c}}}({t}_{0}+{\tau }_{\rm{LC}}) $

, 则$ {N}_{{\rm{c}}}({t}_{0}+{\tau }_{\rm{LC}}) $

= $ {N}_{{\rm{c}}}\left({t}_{0}\right)-1 $

. 本文基于计算所得的L-J体原子运动轨迹, 通过MATLAB程序统计不同时刻原子间的距离, 从而分析得出原子的扩散或汇聚的时间.
在本文中, 假设粒子的扩散、汇聚过程是一个泊松过程, 即每一次粒子的扩散、汇聚过程是相互独立的. 本文将系统内各原子的扩散或汇聚时间进行统计, 图7所示为$\rho^* = 1,\; T^* = 1,\; \varepsilon = 1$

, σ = 1物性条件下L-J体原子扩散或汇聚时间分布图, 其他条件下的单原子L-J体原子扩散或汇聚时间分布与该条件下的情况类似. 通过观察图7可以发现, L-J体的原子扩散或汇聚时间分布符合指数函数分布的密度函数. 因此, 本文应用指数函数分布的密度函数即(14)式去拟合L-J体的原子扩散或汇聚时间分布. 其中, (14)式中的特征时间$ {\tau }_{\rm{LC}}^{*} $

即为原子连通弛豫时间:

图 7 $\rho^* = 1,\; T^* = 1$

, ε = 1, σ = 1条件下L-J体原子扩散或汇聚时间分布
Figure7. Distribution of the dissociation or association time of L-J fluid at $\rho^* = 1,\; T^* = 1$

, ε = 1, σ = 1

$ y= \dfrac1{{\tau }_{\rm{LC}}^{*}} \exp \Big(- \dfrac{t}{{\tau }_{\rm{LC}}^{*}}\Big). $

通过应用(14)式得到了${T}^{*}= 0.85—5$

, ${\rho }^{*}\!= 0.85—1$

, $ \varepsilon = 0.97—1$

, $ \sigma = 0.8—1.3$

条件下单原子L-J体的原子连通弛豫时间$ {\tau }_{\rm{LC}}^{*} $

, 并与黏弹性特征弛豫时间$ {\tau }_{\rm{MD}}^{*} $

和Maxwell弛豫时间$ {\tau }_{\rm{Maxwell}}^{*} $

进行比较, 结果示于图6. 由图6可以看出, 在不同的物性参数和势参数条件下, $ {\tau }_{\rm{LC}}^{*} $

与$ {\tau }_{\rm{MD}}^{*} $

具有相同的变化趋势, 与$ {\tau }_{\rm{Maxwell}}^{*} $

相比, $ {\tau }_{\rm{LC}}^{*} $

对$ {\tau }_{\rm{MD}}^{*} $

的预测更加准确; 当${\rho }^{*},\; \varepsilon ,\; \sigma$

保持不变时(图6(a)), 随着$ {T}^{*} $

的下降, ${\tau }_{\rm{LC}}^{*},\; {\tau }_{\rm{Maxwell}}^{*}$

与$ {\tau }_{\rm{MD}}^{*} $

的差异均逐渐增大, 但$ {\tau }_{\rm{LC}}^{*} $

与$ {\tau }_{\rm{MD}}^{*} $

的结果更加接近; 当$ {T}^{*} $

, $\varepsilon ,\; \sigma$

保持不变时(图6(b)), 随着$ {\rho }^{*} $

的增大, ${\tau }_{\rm{LC}}^{*},\; {\tau }_{\rm{Maxwell}}^{*}$

, $ {\tau }_{\rm{MD}}^{*} $

均呈现上升趋势, 但$ {\tau }_{\rm{Maxwell}}^{*} $

增大的速率最快, 而$ {\tau }_{\rm{LC}}^{*} $

与$ {\tau }_{\rm{MD}}^{*} $

的差异在减小; 当$ {T}^{*} $

, ${\rho }^{*},\; \varepsilon$

保持不变时(图6(c)), 随着$ \sigma $

的增大, $ {\tau }_{\rm{Maxwell}}^{*} $

与$ {\tau }_{\rm{MD}}^{*} $

的差异不断增大, 而$ {\tau }_{\rm{LC}}^{*} $

与$ {\tau }_{\rm{MD}}^{*} $

的差异存在减小的趋势; 当$ {T}^{*} $

, ${\rho }^{*},\; \sigma$

保持不变时(图6(d)), 随着$ \varepsilon $

的增大, $ {\tau }_{\rm{LC}}^{*} $

与$ {\tau }_{\rm{MD}}^{*} $

的结果逐渐接近.
与Maxwell弛豫时间相比, 原子连通弛豫时间在数值大小和变化趋势上均与黏弹性特征弛豫时间更加接近, 但原子连通弛豫时间的计算需要统计所有原子在500个时间步长内的运动轨迹, 分析每一个原子扩散或汇聚时间, 工作量庞大, 计算过程需要消耗大量的时间和计算资源. 因此, 若能通过分析各个原子周围的原子排布, 建立一种间接的反映原子扩散和汇聚过程的方法, 则可避免在统计各原子运动轨迹的过程中消耗大量的时间和计算资源.
2

3.2.预测方法及弛豫时间计算

-->

3.2.预测方法及弛豫时间计算

在上述计算过程中, 发现原子连通弛豫时间对黏弹性弛豫时间的表征更加准确, 但连通性弛豫时间的计算工作量大、耗时长. 本文在分析每个原子周围原子排布的基础上, 应用Kramers逃逸速率理论^[29]来描述原子的扩散和汇聚过程, 建立了一种新的表征单原子L-J体黏弹性弛豫时间的方法, 避免了统计所有原子的运动轨迹时需要消耗大量时间和计算资源.
Kramers逃逸速率理论是常见的描述布朗运动和化学反应扩散的经典方法, 图8所示为Kramers逃逸速率理论示意图, 其中g(r)表示的是原子的径向分布函数曲线, U(r)表示的是原子的能量变化曲线. 当粒子间距离由近处远离, 粒子间能量(U(r)曲线)由第一个谷值点a处上升至能量壁垒b点时, 发生扩散反应, 此时的能量壁垒为$ E_b^{\rm ass} $

; 而当粒子间距离由远处靠近, 对应能量(U(r)曲线)由第二个谷值点c处上升至能量壁垒b点时, 发生汇聚反应, 其能量壁垒为$ E_b^{\rm ass} $

. 该过程中的能量U可通过径向分布函数即$ g\left(r\right) $

曲线获得:

图 8 Kramers逃逸速率理论示意图
Figure8. Schematic diagram of Kramers’ rate theory.

$ U\left(r\right)=-{k}_{\rm{B}}T{\rm{l}}{\rm{n}}\left(g\left(r\right)\right). $

能量壁垒$ E_b^{\rm diss} $

和$ E_b^{\rm ass} $

可分别表示为

$ {E}_{b}^{{\rm{d}}{\rm{i}}{\rm{s}}{\rm{s}}}=U\left(b\right)-U\left(a\right), $

$ {E}_{a}^{{\rm{a}}{\rm{s}}{\rm{s}}}=U\left(b\right)-U\left(c\right). $

扩散反应速率k_diss和k_ass则可通过能量壁垒E_b及能量曲线波峰(壁垒)、波谷曲率相关参数ω_a, ω_b, ω_c和迁移率γ等描述:

$ {k}_{{\rm{d}}{\rm{i}}{\rm{s}}{\rm{s}}}=\frac{1}{2{\text{π}}}\left[{\left(\frac{{\gamma }^{2}}{4}+{{\omega }_{b}}^{2}\right)}^{\frac{1}{2}}-\frac{\gamma }{2}\right]{\left(\frac{{\omega }_{a}}{{\omega }_{b}}\right)}^{\frac{{N}_{{\rm{c}}}}{2}}{\rm{e}}^{-\frac{{E}_{b}^{{\rm{d}}{\rm{i}}{\rm{s}}{\rm{s}}}}{{k}_{\rm{B}}T}}, $

$ {k}_{{\rm{a}}{\rm{s}}{\rm{s}}}=\frac{1}{2{\text{π}}}\left[{\left(\frac{{\gamma }^{2}}{4}+{{\omega }_{b}}^{\gamma }\right)}^{\frac{1}{2}}-\frac{\gamma }{2}\right]{\left(\frac{{\omega }_{c}}{{\omega }_{b}}\right)}^{\frac{{N}_{{\rm{c}}}}{2}}{\rm{e}}^{-\frac{{E}_{b}^{{\rm{a}}{\rm{s}}{\rm{s}}}}{{k}_{\rm{B}}T}}, $

式中N_c为原子连通数, 可通过本文计算得到的径向分布函数得到; ω表示为

$ {\omega }^{2}={m}^{-1}\left|{U''}\right|; $

γ表示为

$ \gamma ={k}_{\rm{B}}T/mD, $

其中D为系统扩散系数, m为单个原子的质量.
假定初始状态系统内的A原子均两两汇聚, 形成A—A分子, 即初始条件为

$ {N}_{{\rm{A}}—{\rm{A}}}\left(t=0\right)=\frac{1}{2}{N}_{{\rm{t}}{\rm{o}}{\rm{t}}{\rm{a}}{\rm{l}}}, $

$ {N}_{{\rm{A}}}\left(t=0\right)=0, $

式中N表示系统内原子的数量.
在环境条件下, L-J体系将发生A—A分子与A原子的相互转化过程, 反应方程式如下:

$ {\rm{A}}—{\rm{A}} \Leftrightarrow 2{\rm{A}}. $

根据可逆反应机理, A—A分子与A原子的物质的量浓度存在以下关系:

$ \frac{{\rm{d}}\left[{\rm{A}}—{\rm{A}}\right]}{{\rm{d}}t}={k}_{{\rm{a}}{\rm{s}}{\rm{s}}}{\left[{\rm{A}}\right]}^{2}-{k}_{{\rm{d}}{\rm{i}}{\rm{s}}{\rm{s}}}\left[{\rm{A}}—{\rm{A}}\right], $

式中[ · ]表示物质的量浓度; k_ass为A原子靠近A—A分子的扩散、汇聚过程(图8)的速率常数; k_diss为A原子远离A—A分子的扩散分散过程(图8)的速率常数. 图9所示为不同温度、密度、势参数条件下, A原子和A—A分子的数量随时间的变化.

图 9 A原子和A—A分子的数量变化曲线　(a)不同温度条件; (b)不同密度条件; (c)不同势能为0时的原子距离条件; (d)不同势能阱深度条件
Figure9. The quantity curves of A and A—A: (a) Different temperatures; (b) different densities; (c) different distances between atoms when potential is 0; (d) different potential well depths.

由图9可知, 在不同的物性条件下, A原子和A—A分子的数量变化均存在相同的规律, 随着时间的增长, A原子的数量不断增加并趋于一个定值, A—A分子的数量不断降低并趋于一个定值, A原子和A—A分子的转换过程趋于平衡状态; 当${\rho }^{*},\; \varepsilon ,\; \sigma$

保持不变时(图9(a)), 随着$ {T}^{*} $

的上升, A—A分子的耗散速度和A原子的生成速度均增大; 当$ {T}^{*} $

, $\varepsilon ,\; \sigma$

保持不变时(图9(b)), 随着$ {\rho }^{*} $

的增大, A—A分子的耗散速度和A原子的生成速度均减小; 当$ {T}^{*} $

, $ {\rho }^{*} $

保持不变时(图9(c)和图9(d)), 无论$ \sigma $

的增大还是$ \varepsilon $

的增大, 均会导致A—A分子的耗散速度和A原子的生成速度减小.
进一步分析A—A分子的耗散速度和A原子的生成速度的变化与黏弹性的关系, 应用(26)式对所有条件下单原子L-J体内A原子数量随时间的变化进行拟合得到扩散和汇聚过程的弛豫时间即为本文提出的预测方法弛豫时间$ {\tau }_{{\rm{m}}{\rm{o}}{\rm{d}}{\rm{e}}{\rm{l}}}^* $