边界条件算法详解

边界条件算法主要目的是实现对活性毁伤材料代表性体积单元周期性边界条件的施加。按模型维度差异，可将边界条件算法分为二维代表性体积单元边界条件算法及三维代表性体积单元边界条件算法。

1.二维代表性体积单元边界条件算法

典型二维代表性体积单元有限元模型如图6.18所示。对边界Γ12，Γ23，Γ14，Γ43上的节点和角节点1，2，3，4，根据周期性边界条件变形协调要求，其横坐标和纵坐标之间应该满足

式中，xΓij和yΓij分别代表四个边界上节点的横坐标和纵坐标；xi和yi分别代表四个角点的横坐标和纵坐标。

图6.18　典型代表性体积单元

1—颗粒；2—基体

除满足上述节点之间的周期性约束之外，还需约束代表性体积单元刚体位移。图6.19所示为典型二维代表性体积单元的约束及加载方式，左下角点约束x和y方向平动自由度，右下角点约束y方向平动自由度，左上角点约束x方向平动自由度。位移或力载荷施加于右下或者左上角点未被约束的自由度。以左上角点为例，若施加y方向载荷，则可实现对材料在单轴拉压情况下的模拟，而载荷的拉、压特性通过调整载荷正负实现。

图6.19　约束及加载方式

边界调价及载荷的施加通过算法编程实现，流程如图6.20所示。模型中顶点节点数为4个，除顶点各边上节点数为n个，则模型边界上共有节点4n＋4个。顶点约束方程数为2个，边界节点的约束方程数为4n个，约束方程总数为4n＋2个。以上数量关系对判断算法准确性具有重要意义。

图6.20　二维代表性体积单元边界条件施加方法

2.三维代表性体积单元边界条件算法

相比于二维代表性体积单元，三维代表性体积单元几何结构更加复杂。从构成元素角度看，立方体结构包含了面、棱、角三类几何特征。三维边界条件算法的实现需基于面节点、棱边节点及角节点来实现。图6.21所示为一立方体结构代表性体积单元，坐标原点位于D点，为不失一般性，设其长、宽、高分别为Wx、Wy、Wz，x、y、z三个方向的位移分别用u、v、w表示。

图6.21　三维代表性体积单元

（1）面节点约束方程。

相对面上的两个对应节点，只需约束方程满足上面三个方向的方程组，即可实现周期性边界条件的施加。

在垂直于x轴的相对面上，约束方程为(www.xing528.com)

在垂直于y轴的相对面上，约束方程为

在垂直于z轴的相对面上，约束方程为

（2）棱边节点约束方程。

对于代表性体积单元中处于平面交线或者交点处的棱边和角点，为其设定周期性约束时，要注意不能重复约束。例如在给面节点施加约束时，若包含了棱边上的节点，再给棱边上的节点施加约束时就会造成这部分节点的重复约束，在有限元计算中会造成过约束错误。棱边节点和角节点也会出现类似情况。因此较好的方法是在给面节点施加约束时，暂时剔除棱边和角节点，只对面中心的节点进行方程约束；在给棱边节点施加约束时，暂时先剔除其中包含的角节点，最后单独对角节点进行约束。对于图6.21所示的代表性体积单元，其12条棱边可分为3类，平行于x轴的AD、BC、FG和EH，平行于y轴的CD、BA、FE和GH，平行于z轴的HD、EA、FB和GC。

对于平行于x轴的4条棱边，约束方程为

对于平行于y轴的4条棱边，约束方程为

对于平行于z轴的4条棱边，约束方程为

式中，i依次为u、v、w，即x、y、z三个方向的位移均需满足此方程。

（3）角节点约束方程。

图6.21中所示的代表性体积单元有8个角节点，选择其中的D点作为基准点，可分别写出其他点关于D点的约束方程：

式中，i为除D点外的其余七个角节点。

此外，与二维代表性体积单元相类似，分析过程中也需要约束模型的刚体位移，载荷通过位移或者力的方式施加。三维代表性体积单元约束及加载方式如图6.22所示，约束方程施加算法流程如图6.23所示。

图6.22　约束及加载方式

图6.23　三维代表性体积单元约束方程施加算法流程

上述二维与三维边界条件算法仅适用于网格周期性排列的情况，角点、棱边、面之间有严格对应关系。周期性网格中，对应边和面上的网格数和分布一致，节点对应，便于通过程序识别和施加相应约束。但一般由于实际结构的复杂性和随机性，越接近材料真实结构和特征，越不利于实现有限元模型周期性网格的划分。因此，也有研究提出针对非周期性网格施加周期性边界条件的方法，即“一般周期性边界条件”，该方法的提出使复合材料微细观有限元模型网格划分不受周期性限制，更适合对复杂结构活性毁伤材料的模拟。