弹性力学中无网格和有限元耦合的元胞自动机算法

陈泽芸 袁卫锋

西南科技大学制造科学与工程学院，绵阳,621010

弹性力学中无网格和有限元耦合的元胞自动机算法

陈泽芸 袁卫锋

西南科技大学制造科学与工程学院，绵阳,621010

结合有限元和无网格算法的优势，提出了一种元胞自动机算法用以求解二维弹性力学问题。该算法将二维模型离散成一系列节点，这些节点被分成有限元群和无网格群。有限元区域被定义在问题的边界附近，其中的任一节点和其周围相邻点的力学关系通过有限元单元建立；无网格区域定义在远离原理问题边界处，其中的节点之间的关系借用有限元中的位移插值概念建立。无论处于有限元区域还是无网格区域，任何一个节点都被置于元胞自动机的框架下进行处理，即节点的位移通过元胞自动机进行求解。与有限元方法相比，所提出的元胞自动机算法无需采用高斯消去法等传统系统求解器，而是通过元胞自动机的自动演化解决问题。依据该算法，有限元和无网格方法可以实现无缝连接。数值算例验证了该算法的新颖性和正确性。

耦合算法；元胞自动机；有限元方法；弹性力学

0 引言

有限元方法的基本思想是将连续求解区域离散成有限个按一定方式相互连接在一起的单元的组合体。有限元方法利用每一个单元内假设的近似函数来分片地表示全求解域上待求的未知场函数[1]。但在求解某些特殊问题，如计算流固耦合、动态裂纹扩展等问题时，有限元网格可能会产生畸变，这种情况不仅需要不断地进行网格重构，而且会严重影响解的精度，这就使得有限元方法在工程应用中有一定的局限性。为解决有限元方法在计算中由于网格畸变引起的缺陷，近年来无网格方法成为研究热点。与有限元方法不同，无网格法是用一组点来离散求解区域，直接借助于离散点来构造近似函数，可以彻底或部分消除网格，不需要网格的初始划分和重构[2]，不仅可以确保计算的精度，而且可以降低计算的难度。现在无网格法已经成功应用于固体力学弹塑性分析[3]以及高速碰撞[4]、断裂损伤力学等问题的研究中[5]，但无网格方法也存在基本边界条件不能直接进行施加[6]、计算量大等问题。将无网格法和有限元法耦合是无网格法处理边界条件的方法之一。

元胞自动机(cellular automata，CA)是一种时间、空间、状态都离散的动力系统模型，不同于一般的动力学模型，它不是由严格定义的物理方程或函数确定的，而是由一系列模型构造的规则构成，凡是满足这些规则的模型都可以算作CA模型。元胞自动机在热扩散[7]、并行计算[8]、紧急疏散[9]、社会交通[10]和复杂的社会经济现象[11]等领域都有广泛的应用。

本文提出了一种将无网格元胞自动机方法与有限元方法相结合的耦合算法，并将其应用到二维弹性力学问题的求解中。

1 二维弹性力学的耦合模型

1.1方法描述

用一组随机离散点对二维弹性力学模型进行离散，如图1所示。模型被分成两个部分，靠近边界的区域为有限元区域，中间部分为元胞自动机区域。有限元区域用有限元方法进行网格划分，利用有限元三角形单元插值函数建立这一区域内节点间的相互关系。元胞自动机区域因为没有网格划分，所以需要先确定其节点的影响域，从影响域内选择与其相关的节点，再建立节点间的相互关系，具体步骤如下。

图1 离散模型示意图Fig.1 Discretization of the model

(1)确定节点影响域的形状和大小。本文选择的节点影响域形状为正六边形，如图2所示。将正六边形影响域分成6个大小完全相等的正三角形，定义图2中一象限内正三角形边与x轴构成的最小夹角θ为正六边形影响域的角度。影响域的大小由模型离散的疏密程度决定，只需保证影响域每个三角形区域内至少有一个点即可。影响域角度由模型离散点的位置决定，所以角度θ可能是0～60°内的任一角度。

图2 节点影响域的大小和位置Fig.2 The hexagon influence zone

(2)选择影响域中与当前节点相关的节点。从节点的正六边形影响域分成的6个正三角形区域中各选一个点作为当前节点的相关节点，若一个三角形区域中不止一个节点，则从中随机选出一个点作为当前节点的相关节点。

(3)建立所选节点间的相互关系。利用有限元方法中的插值概念构造元胞自动机区域内任一节点与其相关节点间的相互关系，具体方法如下。①对正六边形影响域分成的每个三角形区域，用有限元方法中的三角形插值函数建立中间点O与影响域顶点A～F之间位移与力的关系：

(1)

式中,K为正六边形影响域形成的刚度矩阵；u为中间点O和影响域顶点A～F的位移分量；p为中间点O和影响域顶点A～F所受力的分量；下标V代表影响域顶点A～F。

取方程组中计算节点O位移的两个方程，得

KOVuV+KOOuO=pO

(2)

②借助有限元三角形单元插值构造离散节点O及节点1～6和点A～F之间位移与力的关系，如在三角形OAB中有

(3)

同理，借助有限元三角形单元插值可以构造其他5个三角形单元的内节点与三角形顶点的位移函数关系。将这6个三角形单元形成的方程组写成矩阵形式为

(4)

式中，uR、uV分别为真实节点1～6和虚拟点A～F的位移分量。

uO=λpO+γuR

(5)

至此，式(5)为所构造的元胞自动机区域内任一节点与其影响域内相关节点之间位移和力的关系，即CA的局部规则。

1.2与有限元法耦合

需要指出的是，元胞自动机区域内节点的影响域可能包含有限元部分的节点，其邻居也会包含有限元部分内的节点。两个区域交界面上的点也用CA方法选择其邻居并建立交界面上点与其邻居间的关系，交界面上的点必然会将有限元部分和元胞自动机部分中的点选作邻居。用上述方法建立节点间相互关系的同时建立了元胞自动机部分和有限元部分的联系，即无网格部分和网格部分用式(5)进行自然连接，边界上各物理量会通过式(5)传递到CA部分内部。两个部分间不需要增加过渡单元或过渡计算，算法简单。

将有限元部分的求解方程组改写成节点间相互关系：

ux=(px-ωuR′)/Kx,xx=1,2,…,n

(6)

式中,n为模型总的离散节点数;ω为由正六边形影响域刚度矩阵中的第x行系数；uR′为与节点x相关的节点位移分量;Kx,x为刚度矩阵K对角线元素。

由式(5)和式(6)整理可得到模型中所有节点与其相关节点间的关系，写为

ux=apx-buR′

(7)

式中,a、b为常数，由相应的三角形插值函数和有限元部分的刚度矩阵计算得到。

由此，式(7)作为整个模型的元胞自动机自然演化局部规则。

1.3求解步骤

式(7)给出了模型中任意节点与其相关节点的关系，以式(7)作为元胞自动机基本演化规则，离散模型总数n作为一次循环，从任意节点开始以任意顺序进行元胞自动机演化求解。重复循环直到结果收敛，收敛条件是t+1时刻位移与t时刻位移的最大值小于给定误差限，即

(8)

具体演化步骤为：①设置每个节点的位移初始值；②根据t时刻节点位移，计算t+1时刻的节点位移；③重复步骤②至结果满足收敛条件。

因为每个时间步的循环中计算顺序是任意的且每次顺序都不同，这种计算顺序的任意性为并行计算的实现提供有利条件。

1.4应力计算

求得位移之后，用位移法求每个节点的应力。有网格部分，根据已有网格划分用有限元方法计算其单元应力。元胞自动机部分根据每个节点选出的6个相关节点依次选两个点与当前点构成6个三角形单元，通过各节点的已知位移求单元应力。例如，任意三角形OAB的单元应力可写为

(9)

式中,s为对应节点应力矩阵。

任一节点O的节点应力为与节点相关的单元应力的数值平均，即

(10)

式中，m为与节点O相关的单元个数。

2 算例

2.1算例1

四边受均布压力的方板，边长为4 m，均布力大小q=1.0×1010Pa，如图3所示。材料参数为弹性模量E=2.0×1011Pa，泊松比μ=0.3。因为模型关于x、y方向都对称，所以取模型的1/4进行计算。节点影响域的大小和角度取决于离散点的疏密程度，所以模型中各节点的影响域大小和角度均不相同，如图4所示，这使得算法具有良好的适应性。

图3 四边受压方板模型Fig.3 Square plate subjected to distributed forces

图4 节点影响域的大小和角度Fig.4 Sizes and orientations of the nodal influence zones

图5 位移计算结果比较Fig.5 Displacements of the nodes on the boundary

模型总离散点数为1681，其中有限元区域节点数为312，交界面上节点数为144，CA区域节点数为1225。取模型下方y=2的41个节点位移计算结果与理论解比较，结果如图5所示。x方向位移沿负方向线性增大，y方向位移均为-0.007。由图5可以看出，理论解和数值解吻合得很好，最大相对误差小于0.001%。

同样的模型，有限元区域和交界面上离散点数不变，调整CA区域内节点个数，从1225减小到841、529，模型计算结果相差不大。该结果说明，只要保证CA区域内节点能够与有限元区域建立联系，则在一定范围内，CA区域内节点的个数对计算精度影响不大。

2.2算例2

取与算例1模型一样的方板，仅中间开一个半径为0.5 m的孔。同样取模型的1/4作为计算模型，模型总的离散点数为1677。取模型圆弧边界和下边界共51个节点的y方向应力集中系数计算结果与理论结果、ANSYS模拟结果相比较，结果如图6所示。从图6可以看出本文算法的计算结果与理论解之间存在一定误差，但与ANSYS的计算结果十分接近。算例2中的理论解是四边受均布压力的无限大板中间带小圆孔的理论解，而算例2中的方板模型是有限大的，所以求得的应力数值解与理论解有一定的差异。而与ANSYS模拟结果十分接近则可以说明本文提出的耦合算法的正确性和可行性。

图6 应力集中系数计算结果比较Fig.6 Comparison of the stress concentration factors

3 结论

(1)将有限元方法与元胞自动机相结合，提出了一种新的耦合算法。该算法利用插值概念建立任意离散点与其周围节点之间的力和位移的关系，根据元胞自动机的演化求解二维弹性力学问题。计算过程中只需要建立有限元和元胞自动机的相应演化规则，整体基于元胞自动机的演化过程求解，两个部分自然结合，不需要增加任何过渡计算。算例证明该算法简单有效。

(2)虽然该算法需要一段时间更新计算结果，可能会降低计算效率，但是计算过程中任意计算顺序体现了该算法在并行运算中的潜力及优势，为进一步提高计算效率提供了可能。

(3)目前该耦合算法只应用于二维弹性力学，对于三维问题只需要将节点影响域变成正四面体或正八面体。在未来的探索中应该突出有限元及元胞自动机的特点和优势，建立更充分的理论计算基础，以推广其应用到不连续问题和大变形的计算中。

[1] 王勖成,邵敏.有限单元法基本原理和数值方法[M].北京:清华大学出版社,2003:1-48. WANG Xucheng, SHAO Min.The Basic Principle of Finite Element Method and Numerical Method[M].Beijing：Tsinghua University Press,2003：1-48.

[2] BELYTSCHKO T, KRONGAUZ Y, ORGAN D, et al. Meshless Methods: an Overview and Recent Developments[J].Computer Methods in Applied Mechanics and Engineering,1996,139:3-47.

[3] SONG K Z, ZHANG X. Meshless Method Based on Collocation for Elasto-plastic Analysis[C]// Proceedings of Internal Conference on Computational Engineering & Science. Los Angeles, 2000:1-17.

[4] 贝新源, 岳宗五.三维SPH 程序及其在斜高速碰撞问题的应用[J] .计算物理,1997,14(2):155-166. BEI Xinyuan ,YUE Zongwu. A Study on 3-Dimensional SPH[J]. Computational Physics,1997,14(2):155 -166.

[5] 周进雄,李梅娥,张红艳，等.再生核质点法研究进展[J].力学进展,2002,32(4):535-544. ZHOU Jinxiong, LI Mei’e, Zhang Hongyan, et al. Advance in Reproducing Kernel Particle Method[J].Advance in Mechanics,2002,32(4):535-544.

[6] BELYTSCHKO T, KRONGAUZ Y. Meshless Methods: An Overview and Recent Developments[J].Comput Methods Appl. Mech. Engrg.,1996,139:3-47.

[7] YUAN W F, TAN K H. Heat Transfer in Composite Beams Using Combined Cellular Automata and Fibre Model[J]. Computers, Materials & Continua,2009,13(1):49-62.

[8] 杨吉新. 并行元胞单元法[J]. 计算力学学报,2006,23(1):97-100. YANG Jixin. Parallel Cellular Element Method[J]. Computational Mechanics,2006,23(1):97-100.

[9] YUAN W F, TAN K H. A Model for Simulation of Crowd Behaviour in the Evacuation from a Smoke-filled Compartment[J]. Physica A: Statistical Mechanics and Its Applications,2011,390(23/24):4210-4218.

[10] 郑亮, 马寿峰, 贾宁. 基于驾驶员行为的元胞自动机模型研究[J]. 物理学报, 2010,59(7):4490-4498. ZHENG Liang,MA Shoufeng, JIA Ning.The Cellular Automaton Model of Traffic Flow Based on the Driving Behavior[J]. ACTA Physica Sinica,2010,59(7):4490-4498.

[11] 颜伟, 唐德善. 区域经济聚集效应的元胞自动机模拟[J]. 商场现代化, 2006(4):254. YAN Wei, TANG Deshan. Regional Economic Aggregation Effect of Cellular Automata Simulation[J]. Market Modernization,2006(4):254.

(编辑王艳丽)

ACouplingCAAlgorithmofMeshlessandFEMinElasticity

CHEN Zeyun YUAN Weifeng

School of Manufacturing Science and Engineering，Southwest University of Science and Technology，Mianyang，Sichuan，621010

A CA was developed to combine the advantages of FEM and meshless method which was used to analyze 2-dimensional elastic problems. By the CA, a 2-dimensional domain was discretized into a grid of nodes. These nodes were distributed randomly in the domain and they were divided into the FEM group and the meshless group. In the FEM zone which was mostly defined near the boundaries of the domain, FEM elements were used to establish the relationship between one node and other adjacent nodes. In the meshless zone which was normally away from the boundaries, the adjacent nodes were related by employing the concept of displacement interpolation used in FEM. However, each nodes, wherever they were in the FEM zone or meshless zone, were dealt with under the CA frame, that is, the unknown displacements of each nodes were evaluated by CA algorithm. Unlike FEM, the proposed CA solver worked out the results through automatic evolution, instead of Gaussian elimination or other conventional methods. Based on the current CA algorithm, FEM and meshless may be merged seamlessly. The novelty and the correctness of the current approach were proven by numerical examples.

coupling algorithm; cellular automaton(CA); finite element method(FEM); elasticity

2016-11-14

国家自然科学基金资助项目(11472232)；四川省教育厅教改项目(14JGCX07);制造过程测试技术省部共建实验室基金资助项目(13zxzk05)

TB121

10.3969/j.issn.1004-132X.2017.17.018

陈泽芸，女，1991年生。西南科技大学制造科学与工程学院硕士研究生。发表论文5篇。E-mail:464252016@qq.com。袁卫锋，男，1970年生。西南科技大学制造科学与工程学院研究员。