多边形邻域元胞自动机

唐 强，袁卫锋

（西南科技大学制造科学与工程学院制造过程测试技术省部共建教育部重点实验室，四川 绵阳 621010）

1 引言

自然界中，许多学科领域内的问题可用代数、微分或积分等数学方程描述。少数简单或能简化的问题有解析解。多数问题由于方程的特征或求解域复杂的几何形状而无法得到解析解，这类问题主要用数值的方法来近似求解。应用广泛且解决了大批有重大意义问题的有限单元法是一种实用的数值求解工具[1]。有限单元法的思想是将连续的求解域离散成有限多个且按一定方式联结的单元，在每个单元上分片地近似待求解问题，用最小位能原理将问题离散成线性方程组，求解得到整个问题域上的近似解。有限元是基于网格近似的数值方法，在求解大变形、网格畸变或需要网格重构时求解精度降低甚至求解失败。针对有限元法的缺陷，计算力学界兴起了对无网格法的研究。无网格法是将连续的求解域离散成有限多个节点，对每个节点给定局部影响域，在每个局部影响域上近似待求解问题，用不同的加权余量法将问题离散成线性方程组。文献[2]用移动最小二乘法（MLS）构造近似函数，将由随机点离散成的问题域用规则的网格覆盖，在规则的背景网格上采用伽辽金法将问题离散。为了避免用背景网格，文献[3]在局部域上用MLS构造近似函数，用局部弱式伽辽金法将问题离散，形成了无网格局部Petrov-Galerkin（MLPG）法。MLPG法是完全不需要网格的方法。目前，已经提出了近30余种无网格法。无网格法在断裂力学[4]、大变形[5]等网格变形大或需要网格重构的领域得到了应用。无网格法在用加权余量法时，通常用数值积分形成线性方程组，增加了计算量。文献[6]提出的无网格元胞自动机算法，基于元胞自动机思想建立一种新的无网格法并应用于二维弹性力学问题。无网格元胞自动机算法构造简单，数值模拟时不需要数值积分，较传统无网格法减少了计算量。

元胞自动机（CellularAutomaton，简称CA）是文献[7]在研究生物自繁殖现象时提出来的。元胞自动机是一组规则排布在元胞空间上且具有有限多个状态量的元胞根据自身及其邻居上一时刻的状态量在离散的时间步上按照指定的演变规则进行状态更新的系统。在某一时刻，元胞的状态量为其有限多个状态量中的一种，例如阀门在某时刻非“关”即“开”，该时刻任意元胞的状态可表示为：

式中：s—元胞状态；i—第i个元胞；t—时间步；N—邻居。

元胞自动机根据元胞的特点建立演变规则，能更真实地反映问题的本质。目前，元胞自动机已广泛应用于应急疏散[8]、热传导[9]。文献[10]将元胞自动机思想与有限元法结合形成了元胞单元法。元胞单元法将节点及单元同时考虑成元胞，并将节点和位移均考虑成元胞的状态，元胞状态的更新规则按照力与位移相互转化来确定。提出的多边形邻域元胞自动机算法，将元胞自动机思想应用于无网格法，形成一种简单的无网格法。该方法将问题域离散成随机分布的场节点，仅将位移量考虑成元胞的状态量，借用有限元思想建立元胞与邻居间的演变规则，用并行运算及遗传因素提升计算效率，具有求解大变形、网格畸变问题的应用前景。

2 算法构造

二维弹性力学问题在用多边形邻域元胞自动机算法求解时，将问题域离散成随机分布的场节点，将任意一个离散的场节点考虑成一个元胞（Cell），将待求解的问题域考虑成元胞空间，如图1所示。

图1 离散模型Fig.1 Discrete Model

元胞仅会“感知”与其邻近的元胞及直接作用在该元胞上的外力。元胞的状态量为矢量，分别为x、y方向的状态量u、v。元胞对邻近元胞的“感知”用邻居元胞状态量的加权求和近似表示为式（1）。元胞的状态根据式（1）指定的演变规则在离散时间步上演化，直到元胞空间内的元胞都达到稳定状态。

式中：fix、fiy—第 i个元胞 x、y 方向上的节点力；wi、ωi—第 i个元胞x、y方向节点力权重；φ、ψ—邻居元胞状态量权重；n—包含第i个元胞在内的邻居个数。

根据模型中元胞的分布情况指定合适的“感知”半径，确定出元胞的邻居。因元胞的分布情况不同，每个元胞的邻居个数可能存在差异。根据邻居个数的不同，建立边数不同的正多边形邻域。正四边形邻域与正六边形邻域，如图2所示。

图2 正多边形邻域Fig.2 Regular Polygon Influence Zone

借用有限元思想建立式（1）的演变规则。以正四边形为例说明演变规则的构造过程：

（1）将正四边形邻域考虑成一个局部有限元模型，借助有限元求解方程Ka=P建立中心元胞i与顶点A、B、C、D间的关系。只考虑邻居元胞对中心元胞的影响，选取有限元求解方程中以中心点i为对角元的两个线性方程：

式中：k1x=［k13，k15，k17，k19］、k2x=［k23，k25，k27，k29］、k1y=［k14，k16，k18，k110］、k2y=［k24，k26，k28，k210］—局部总刚元素组合，up=［uA，uB，uC，uD］T—正多边形顶点 x 方向的状态量；vp=［vA，vB，vC，vD］T—正多边形顶点y方向的状态量。

（2）分别在正多边形中的每个三角形内插值出对应邻居元胞的状态量，如式（3）表示。每个邻居点x方向的状态量，y方向状态量只需将式（3）中的u替换成v。

利用各三角形内插值函数和为1，并将式（3）改写成矩阵形式：

式中：us=［u1，u2，u3，u4］T—邻居元胞 x 方向的状态量；vs=［v1，v2，v3，v4］T—邻居元胞 y 方向的状态量，E=［1，1，1，1］T。

（3）用式（4）求出正多边形邻域顶点的状态量由邻居元胞状态量表示：

（4）式（5）代入到式（2）中简化得到邻居元胞对中心元胞的影响：

式（6）表示的演化规则中，并未考虑元胞自身上一迭代步对应方向上的影响。将其右侧部分考虑成非自身因素的影响，用、表示。对（6）式加入自身因素的影响，可称之为遗传因素，如式（7）所示。

式中：α—遗传系数，该系数对演变步数有较大影响。

图2中1、2、3、4元胞与正四边形邻域顶点重合时：

简化式（7）得到该特殊情况下邻居元胞对中心元胞影响的表达式，就是式（3）中直接将多边形顶点状态量改写成邻居元胞的状态量：

式（9）说明当邻居元胞与多边形邻域顶点重合时，本算法的演变规则直接由局部有限元总刚获得。当离散域中每个元胞都选择有限元网格划分后与之相联结的单元上的元胞作为邻居，多边形邻域元胞自动机算法将完全退化成有限元。

传统无网格法在施加本质边界条件时，通常用特殊的处理办法，例如罚函数法。本算法借助有限元插值思想建立演变规则，有限元是本算法的特殊情况，两种算法能自然耦合。求解问题时，本算法求解域边界上元胞的邻居与多边形邻域顶点重合，直接施加本质边界条件。边界层采用有限元而模型内部采用多边形邻域元胞自动机算法会形成的过渡区用元胞自动机算法。

该算法求解二维弹性力学问题时的步骤为：

（1）将问题域离散成随机分布且总数为NCA的元胞，给出材料信息、边界条件、误差ε、初始元胞状态u=0且v=0、遗传因素α、迭代步 t=0；

（2）指定元胞“感知”半径，搜索邻居元胞，确定多边形邻域，计算式（7）中各权重系数；

（3）按随机序列更新元胞状态；按照式（7）的演变规则先更新第i个元胞x方向的状态量ui，若指定了该点x方向的边界位移为u0，则令ui=u0；然后以相同方式更新vi；若NCA个元胞都更新完毕，则迭代步t=t+1，进行下一步；

3 算例

3.1 带孔方板

有限大的四方受压带孔方板，如图3所示。其解析解可用无限大（即 l＞＞a）带孔方板近似：

图3 受压带孔方板Fig.3 Pressurized Perforated Square Plate

建立受压方板1/4模型，如图4所示。模型中均布载荷q=-2000N/m2，杨氏模量 E=2×105Pa，泊松比 υ=0.3，模型尺寸 l=4.0m，a=0.8m，板厚1.0 m，α=-0.42，元胞总数1721个。虚线表示模型初始位置，位移放大了10倍。

图4 受压方板1/4模型Fig.4 Model of 1/4 Pressurized Perforated Square Plate

模型求解是对求解域中的每个元胞进行状态更新。求解是一个动态演化的过程。选取x=y=l处的点，作出了该点y向的状态量的演变趋势，如图5所示。随着迭代步数的增加，该点的状态量逐渐趋于稳定。经验证，在-1＜α＜1内，随着α的减小，模型演化所需的求解步数减小。

显示CA数值结果与有限元结果基本一致，而与解析解存在较大的误差。原因是解析解是l＞＞a带孔方板的解，如图6所示。

图5 位移变化趋势Fig.5 Variation Trend of Displacement

图6 带孔方板位移结果Fig.6 Displacement of Pressurized Perforated Square Plate

3.2 悬臂梁

悬臂梁，如图7所示。自由端x=l处施加抛物线型表面力为梁截面的剪应力在该处的取值

图7 悬臂梁Fig.7 Cantilever Beam

该问题存在解析解：

式中：p—表面分布力的积分。

建立悬臂梁离散模型，如图8所示。载荷p=-1000N/m2，杨氏模量 E=2×105Pa，泊松比 υ=0.3，模型长 l=4.0m，模型宽 d=1.0m，板厚1.0m，α=-0.42。问题域离散成分布相对均匀且总数为1874个元胞，由于元胞分布相对均匀，因此大多数元胞确定为了正六边形邻域。位移缩小为0.25。

图8 悬臂梁模型Fig.8 Model of Cantilever Beam

悬臂梁同样具有图5所示的变化趋势。该算法与有限元、解析解结果基本吻合，如图9所示。

图9 悬臂梁位移结果Fig.9 Displacement of Cantilever Beam

3.3 并行计算

多边形邻域元胞自动机算法求解问题时，仅需存储邻居元胞的相关信息，不需要合成总刚。这较有限元及传统的无网格法节约了存储空间。本算法虽节省了存储空间，但用迭代时间步求解，存在迭代步数较多，演化时间长的问题。本算法元胞的演变仅与其邻近元胞的状态相关，元胞的演变具有高度的并行性。按照并行运算思路，将前述元胞状态更新过程修改如下：

（1）离散问题域，设置初始条件；（2）选择邻居，计算式（7）中的权重；（3）将NCA个元胞随机分配到m个处理器中，每个处理器同时按式（7）更新分配到的元胞的状态；若所有处理器中的元胞状态都更新完毕，则迭代步t=t+1；（4）若满足收敛条件，则求解结束；否则，重复第（3）步。

用单核、双核及四核更新元胞数为7267个的悬臂梁模型，双核求解用时约为单核求解用时的50%，如图10所示。四核求解用时约为单核用时的25%。

图10 并行计算Fig.10 Parallel Computation

4 结论

（1）多边形邻域元胞自动机算法用有限元插值思想建立元胞间的演变规则，演变规则构建简单计算量小且能与有限元自然耦合，本质边界条件施加简单。

（2）本算法求解二维弹性力学问题时，将问题域离散成随机分布的元胞，用元胞状态的演变来获得稳定的数值解。算例表明，该算法在元胞分布相对均匀时能获得收敛且误差较小的数值解。

（3）逐个更新元胞的状态量，整体问题求解时间较长。根据元胞自动机高度并行的特点，采用并行运算大大提升了该算法的计算效率。考虑遗传因素影响，同样提升了计算效率。

（4）该算法求解问题时不需要网格，是一种简单的无网格法。算法在求解大变形、网格畸变等等领域具有应用前景。