有限元方法在二维纳米材料中的运用.doc

二维纳米级聚合物领域的有限元分析
摘要：
在研究二维纳米级聚合物的有限元方法的两种类型是化学键的因素和兰纳-琼斯势能函数的因素，而这两种因素构成了原子间和分子间的力场基础。纳米尺度的有限元方法可以用来模拟高分子的运动状况。这个数值过程包括三个步骤。首先要创造一个非晶格的随机运动同时伴随一个弛豫过程的聚合物邻域场。然后在这个邻域中形成有限元网格，化学键是由化学键建模的元素。如果两个非成键原子和单体之间的距离短于兰纳-琼斯吸引或者排斥的作用范围，那么它们之间就要加一个兰纳-琼斯势能因子。最终，外部载荷和边界条件的应用可以将高分子链变形过程逐步模拟出来。在绝活物变形过程中，要将不合理的兰纳-琼斯因子剔除，而新形成的兰纳-琼斯因子在每一步中被插入到聚合物邻域中。进一步,在整个变形过程中观察到聚合物链迁移, 纳米级孔隙,孔隙聚结和裂纹的产生过程。

1. 引言
在目前研究纳米结构的分析方法中，分子动力学是最受欢迎的方法之一。它被用于计算聚合物、纳米管和纳米管增强聚合物的模量和强度。运用分子动力学，对温度引起的高频分子热振动和静态变形的模拟可以同时进行。然而，分子热振动的频率范围在1015数量级。分子动力学的模拟可能仅仅提供发生在规模过快（微小或者毫微秒）的变形过程。因此，相应的应变速率会比实际工程实践的数值要大的多。
蒙特卡罗模拟是被用来研究纳米级聚合物变形的另一种方法，这种方法基于统计力学的角度。变形是指材料在受到外界压力时材料表现出的一系列应变增量的现象。压力是一个表达势能和温度的函数。实际上，每一个应变增量都会涉及到Metropolis能量最小周期。已有前人研究发现，当他们使用蒙特卡罗方法模拟非晶态聚合物的变形时，由于过程的收敛速度太慢,他们只能获取短暂时间尺度的变形信息。实验的应变率范围是108到109/s，与分子动力学的数值相近。
对于大多数工程应用,在较大的温度范围内材料的模保持为恒定的常数。这就意味着分子热运动的运动频率不是弹性变形的主要因素。那么问题就来了：从一个分子计算模拟模型的高频率热力学振动中有可能得到材料的弹性形变响应么？
1986年时，Theodorous和Suter提出了一个详细的关于研究分子振动对玻璃态聚丙烯弹性常数贡献的理论。他们得到了这样的结论，在高频率分子热振动中只有一个是可以忽略不计的贡献。他们建立了一个为研究非晶态聚合物的多维数据集的平衡原子场。每一个原子都被分配了一个基于均匀应变假设的微小初始位移。因为分配的原子场不是一个平衡场，所以会存在一个能量最小化的过程。模量是材料应对变化过程中由势能衍生出来的一个量。相比于使用分子动力学方法以及蒙特卡罗方法，这个方法节省了大量的计算过程和计算量；然而，它存在一定的局限性：确定复杂形状材料的原子场域的初始位移以及确定复杂的负荷任然是一个困难的障碍所在。
而有限元方法是一种高效能、常用的计算方法。有限元法在早期是以变分原理为基础发展起来的，所以它广泛地应用于以拉普拉斯方程和泊松方程所描述的各类物理场中（这类场与泛函的极值问题有着紧密的联系）。自从1969年以来，某些学者在流体力学中应用加权余数法中的迦辽金法(Galerkin)或最小二乘法等同样获得了有限元方程，因而有限元法可应用于以任何微分方程所描述的各类物理场中，而不再要求这类物理场和泛函的极值问题有所联系。有限元法的基本思想：由解给定的泊松方程化为求解泛函的极值问题。
在本文中纳米级的有限元模型是基于前人基础后有所发展的方法。首先，要建立纳米元素的刚度矩阵，　然后使用该模型来模拟变形和计算高分子领域的应力-应变曲线。在整个变形的聚合物中，高分子链的迁移和纳米级孔隙和裂纹的产生是观察到的。这种新方法可以在分子运动模型中对从静态分子位移从超高频率的热振动过程进行去耦合操作。高应变率的问题---不可避免的分子动力学模拟可以因此消除掉。它的作用效率可以跟Theodorous和Suter的方法差不多。并且，方法的用法并不局限于材料的形状等条件的制约。

2. 制定二维纳米级有限元
2.1原子力场
原子力场的核心是原子力学。原子力场来源于势能，并与变形的化学键合分子间作用力相互影响。原子力场的设置是整个模拟计算过程的基础条件。
2.2化学键
大部分材料的局部几何化学键的信息可以被定义为键长r，键角

θ以及二面角Ф，正如图1所示。化学键的变形包括键的延伸（化学键键长r的变化），键角的弯曲（键角

θ的变化）以及二面角的扭转（二面角Ф的变化）。有时候还需要加一些附加条件来表示分子结构和分子的变形过程。在二维数值模型中，只需要对键的延伸与键角的弯曲进行讨论。

势能的最简单的表达式与键延伸谐波函数相关。键延伸的谐波函数提出了键的伸长或者缩短与恢复力之间的线性关系。这个势能的表达式如下所示：

这里r0是两个原子之间的初始距离而kr是键的刚性程度。从而恢复力与键延伸之间的关系可以推出如下关系式：

表达势能的另一种改进形式是莫尔斯势能表达式，可以这么来表达：

总结上述的各个式子，键长增长与恢复力之间的关系可以用下列式子表示：

图2表示的是来源于莫尔斯势能的力与位移之间的关系。这是一条非线性的曲线。键额刚性程度kr是这条曲线的斜率，从这个可以体现出键的长度的变化情况。如果键长的增长比较小，恢复力-键场增加的曲线大约呈线性变化。在处理谐波势能与莫尔斯势能关系是，如果键的应变很小也可以得到类似的结论。随着键的应变的扩大，莫尔斯势能就会变得越加的准确。莫尔斯势能的表达式最初是由量子力学推导而来的，但是反过来它有成了支持量子力学的一个例证。
谐波势能通常能够用来描述分子力学中两个临近化学键的交互作用。它提供了回复过程中弯矩和弯曲角度之间的关系。正如下面式子给出的关系：

在上式中

θ0是平衡键角。弯矩和弯曲角度的关系可以写成如下关系式：

2.3分子间相互作用
兰纳-琼斯势能函数被用于描述两个原子或者单体之间不以化学键结合的相互作用关系。兰纳-琼斯势能函数可以被描述为如下方程：

在此处ε表示势阱的深度，σ表示在零势能距离时两个非键合的原子或者单体之间的距离。两个非键合的原子或者单体之间的兰纳-琼斯势能力（吸引力或者排斥力）可以表示为：

2.4二维纳米尺度下的高聚物场的有限元模型
有限元模型的设立考虑了两种类型的影响因子：包括化学键因子和兰纳-琼斯因子。
2.4.1 化学键元
在图3中显示了标准的碳氢化合物高聚物的分子结构。基于恢复力和力矩的本质作用，通过弹性连接构成了以共价键链段为模型的弹性棒。弹性棒能够传递恢复力；而弹性连接是恢复弯曲的基础所在。在图中的链被划分为化学键因子（图中以虚线标注）。每一个因子有一个弹性棒以及两个弹性连接组成，这两个弹性连接分布在每个弹性因子的两端。化学键因子的拉伸刚度由恢复力与键延伸之间的关系所确定，这在之前的等式中已经提出。同时，之前也推出了共价键的抗弯刚度与弹性接头的抗弯刚度是相等的。

在图4中显示了化学键因子的节点力和位移。在图中展示了在每一个节点有三个节点力：对于i分别是Fxi，Fyi和Mi，对于j分别是Fxj，Fyj和Mj；每一个节点有三个节点位移：对于i有ui，vi和

θi，对于j有uj，vj和

θj。那么节点力的增量和节点位移的增量在x方向的关系是：

上述kr表示延伸刚度。如果在这里引入谐波势能的概念，kr表示先前提到的延伸刚度；如果这里引入莫尔斯势能，那么kr则表示恢复力与键延伸曲线的斜率，正如图2中的曲线那样。

假设节点i和j的键刚度为k

θi和k

θj，在键中点处的键刚度为2k

θi和2k

θj。节点弯矩增量诱导在每一个加载过程中可以用以下表达式表达：

在这里Δ

θL可以表示为以下式子（其中L表示为因子的长度）：

根据静力平衡条件,节点力在y方向增量为：

结合之前推出的等式，可以得到节点力与节点增量之间的关系：

因此因子的刚度矩阵可以写成如下式子：

2.4.2 兰纳-琼斯键元
兰纳-琼斯键是一种非线性的弹性模型。图5展示了两个非键合原子之间兰纳-琼斯势能力对距离关系曲线。如果两者之间的距离小于1.244σ，这条曲线的斜率会是正的；若两者间的距离大于1.244σ，纳米这条曲线的斜率将会小于零。为了能够反映实际情况下的兰纳-琼斯键，这条曲线的斜率被定义为该因子的刚度。然而，负的因子刚度可能会导致整体的刚度矩阵的性质不是正定的。为了克服这个困难，我们可以采用如下的方法：如果力-距离曲线的斜率是正的，那么我们可以通过斜率来确定该因子所具有的刚度。但如果这个曲线斜率是负的，我们会使用一个较小的正扰动作为该因子的刚度而不是它原始的负值。从函数角度来考虑，这个兰纳-琼斯键因子可以定义为如下的函数：

上述F值表示两个非键合原子之间的 兰纳-琼斯键力，r是两个非键合原子之间的距离，Δ是一个小扰动项，用来避免刚度矩阵的奇异性。
一个小扰动项更换的负的因子刚度可能会产生错误数值计算结果。为了避免直至替代过程产生的错误，在每一个加载过程中都加上一个弛豫过程，正如图5所示。这边提到的弛豫过程在后续过程中会进行介绍。

从图5中我们可以看到，在斜率为负值的区域兰纳-琼斯相关作用随着距离的增加而迅速减小。当距离r远远大于零势能处的距离σ时兰纳-琼斯变得可以忽略不计。在数值计算过程中，非键合原子之间的距离在加载步骤就开始进行检测。如果两原子之间的距离比兰纳-琼斯键力的临界距离（2.5σ），那么久会有产生一个兰纳-琼斯键元插入在他们之间。在高聚物发生变形过程中，不合理的错误兰纳-琼斯键元将会被剔除，新的复合实际情况的兰纳-琼斯键元会被补充进来。
2.5 加载-弛豫过程
人们通常使用类似有限元方法来形成刚度矩阵。如果采用莫尔斯势能，化学键的元是一个非线性元素。兰钠-琼斯元是高度非线性的元。再者，当力-距离曲线的斜率为负值时，我们会使用一个较小的位移代替实际的力与兰纳-琼斯键的关系。这个替代过程违反了兰纳-琼斯键的原始的关系。这种处理过程产生的错误必须在数值计算过程进行剔除。
为了能纠正这个错误，我们才采用了加载-弛豫过程。假设节点位置与节点力在第k个加载过程为Xk-1和Fk-1,第k个加载过程的加载力为ΔFk，那么第k次的节点位移增量为ΔUk。
第一步：计算出节点位移增量ΔUk的式子为：

第二步：计算节点位置：

第三步：计算实际加载力ΔF'k。无论用谐波形式的能量，还是莫尔斯能，节点力都能够得以计算出来。减去轴节力之前节点力的加载步之后加载步骤,一个增量轴节的力量。减去之前的节点力，加载从加载后一步一步,一个派生实际的增量，其中一部分是兰纳-琼斯相互作用。
第四步：检查错误。如果，第k次的节点位移增量的ΔUk是可以接受的，反之则是不能接受的。必须采取别的方法
第五步：弛豫过程。令，并计算其结果。重复步骤2到5，直到满足条件：

3高聚物邻域的纳米元分析
3.1高聚物场的产生
在图六中显示了一条简单的高聚物链段，该链段是随机排列的。开始单体是随机放置在聚合物领域。化学键长度限制在范围之间在Lmin和Lmax之间，而第二单体必须位于图六所示的阴影部分处。键角的选择是随机的。第三单体也以同样方式生成。随着这个过程进行会产生高聚物链段。

聚合物的形成过程正如如下步骤：首先，单体随机分布在所选区域的范围里。其次，加入更多的单体，产生高聚物链段，并且链段发生生长，形成高聚物场。这种随机的高聚物场并不是实际意义的高聚物，必须经历弛豫过程。因此在第三步在高分子领域运用不平衡力的计算每个原子和应用，后用来进行弛豫过程。最终形成一个原子域。下一步是检查的内部压力，对高分子领域这参数与单体密度有关。如果密度太大，会形成聚合物抗压的一个状态。如果密度太小，聚合物会处在一个拉升的状态。
通过下面的过程可以得到聚合物邻域的局部应力。在聚合物邻域中进行一次假象的切割。所有连接切割两边的化学键和兰纳-琼斯键都要被消除。在化学键和兰纳-琼斯键消除之前，每一个的原子灰处在平衡位置，这就意味着应用每个原子的合力为零。但当化学键和兰纳-琼斯键移除后，单体附近就不再处于平衡位置。在每个部分的非平衡力量计算。这些相互作用的合力被分成左右两边的牵引，因此，内部压力是派生。
在本文中提供的例子,应力状态在聚合物领域的边界进行检查。单体在聚合物领域之外，生成符合周期性边界条件。我们沿着四方的虚构的削减，边界上的聚合物研究牵引作用。随后计算压力对体系的作用。通过调整聚合物实际密度，得到聚合物领域初始的零压力状态。
3.2聚合物纳米级元模型
聚合物链建模为化学键链。这些键连接附近的碳氢化合物单体，包括同一链和不同的链。在聚合物的每一个加载步骤，在聚合物变形过程中，不断地进行剔除和插入兰纳-琼斯元。这些参数已有其他研究者进行计算：

3.3聚合物场的刚度和强度
聚合物邻域在单向加载变形是模拟的。变形过程和破坏机理可以理解如下。在低压力，单体在附近不牵强附会的仓位,化学和兰纳-琼斯键近似线性力-位移关系。因此,整个聚合物领域将是一个近似线性系统。这可以解释最初的应力-应变的线性关系。聚合物链之间的距离增加与应变的增加,兰纳-琼斯元减少。因此聚合物刚度减少。兰纳-琼斯元持续减少,聚合物领域变得不那么均匀，从而出现孔隙。最终的结构聚合物领域失败。虽然没有对比实验是可能的，这个简化的二维高分子场的潜力，使聚合物纳米变形。
4. 结论
在本文中提出了纳米元处理二维材料的方法。在计算二维纳米元的时候主要考虑两类纳米元：化学键元与兰纳-琼斯元。通过模拟一个简化的二位高分子场来生成基于宏观应力与应变的相互作用和分子力场。当我们一步一步地使高分子进行变形，可以观察到纳米级孔隙、高分子的凝聚和裂缝形成现象。而实际的高聚物是三维结构的，因此不可能建立一个互相交缠的聚合物场，所以才是使用的简化二维场。为了定量分析刚度、强度和失效机理，聚合物的三维结构模型的建立是必须的，而这也是目前的研究方向。
参考文献：
【1】Frankland SJV, Brenner DW. Molecular dynamics simulations of polymer-nanotube composites. Materials Research Society Symposium Proceedings 1999:593.
【2】Frankland SJV, Caglar A, Brenner DW, Griebel M. Reinforcement mechanism in polymer nanotube composites: simulated non-bonded and cross-linked system. Materials Research Society Symposium Proceedings 2000:594.
【3】Ozaki T, Iwasa Y, Mitani T. Stiffness of single-walled carbon nanotubes under large strain. Physical Review Letters 2000;84:1712.
【4】Boyce MC, Socrate S, Llana PG. Constitutive model for the finite deformation stress-strain behavior of poly(ethylene terephthalate above the glass transition. Polymer 2000;41:2183.
【5】陈大勇. 高分子材料凝聚态多尺度贯通有限元分析[D]. 新疆大学, 2006.
【6】人宪. 有限元法基础[M]. 国防工业出版社, 2004.