唐 強(qiáng)，袁衛(wèi)鋒

（西南科技大學(xué)制造科學(xué)與工程學(xué)院制造過(guò)程測(cè)試技術(shù)省部共建教育部重點(diǎn)實(shí)驗(yàn)室，四川 綿陽(yáng) 621010）

1 引言

自然界中，許多學(xué)科領(lǐng)域內(nèi)的問(wèn)題可用代數(shù)、微分或積分等數(shù)學(xué)方程描述。少數(shù)簡(jiǎn)單或能簡(jiǎn)化的問(wèn)題有解析解。多數(shù)問(wèn)題由于方程的特征或求解域復(fù)雜的幾何形狀而無(wú)法得到解析解，這類問(wèn)題主要用數(shù)值的方法來(lái)近似求解。應(yīng)用廣泛且解決了大批有重大意義問(wèn)題的有限單元法是一種實(shí)用的數(shù)值求解工具[1]。有限單元法的思想是將連續(xù)的求解域離散成有限多個(gè)且按一定方式聯(lián)結(jié)的單元，在每個(gè)單元上分片地近似待求解問(wèn)題，用最小位能原理將問(wèn)題離散成線性方程組，求解得到整個(gè)問(wèn)題域上的近似解。有限元是基于網(wǎng)格近似的數(shù)值方法，在求解大變形、網(wǎng)格畸變或需要網(wǎng)格重構(gòu)時(shí)求解精度降低甚至求解失敗。針對(duì)有限元法的缺陷，計(jì)算力學(xué)界興起了對(duì)無(wú)網(wǎng)格法的研究。無(wú)網(wǎng)格法是將連續(xù)的求解域離散成有限多個(gè)節(jié)點(diǎn)，對(duì)每個(gè)節(jié)點(diǎn)給定局部影響域，在每個(gè)局部影響域上近似待求解問(wèn)題，用不同的加權(quán)余量法將問(wèn)題離散成線性方程組。文獻(xiàn)[2]用移動(dòng)最小二乘法（MLS）構(gòu)造近似函數(shù)，將由隨機(jī)點(diǎn)離散成的問(wèn)題域用規(guī)則的網(wǎng)格覆蓋，在規(guī)則的背景網(wǎng)格上采用伽遼金法將問(wèn)題離散。為了避免用背景網(wǎng)格，文獻(xiàn)[3]在局部域上用MLS構(gòu)造近似函數(shù)，用局部弱式伽遼金法將問(wèn)題離散，形成了無(wú)網(wǎng)格局部Petrov-Galerkin（MLPG）法。MLPG法是完全不需要網(wǎng)格的方法。目前，已經(jīng)提出了近30余種無(wú)網(wǎng)格法。無(wú)網(wǎng)格法在斷裂力學(xué)[4]、大變形[5]等網(wǎng)格變形大或需要網(wǎng)格重構(gòu)的領(lǐng)域得到了應(yīng)用。無(wú)網(wǎng)格法在用加權(quán)余量法時(shí)，通常用數(shù)值積分形成線性方程組，增加了計(jì)算量。文獻(xiàn)[6]提出的無(wú)網(wǎng)格元胞自動(dòng)機(jī)算法，基于元胞自動(dòng)機(jī)思想建立一種新的無(wú)網(wǎng)格法并應(yīng)用于二維彈性力學(xué)問(wèn)題。無(wú)網(wǎng)格元胞自動(dòng)機(jī)算法構(gòu)造簡(jiǎn)單，數(shù)值模擬時(shí)不需要數(shù)值積分，較傳統(tǒng)無(wú)網(wǎng)格法減少了計(jì)算量。

元胞自動(dòng)機(jī)（CellularAutomaton，簡(jiǎn)稱CA）是文獻(xiàn)[7]在研究生物自繁殖現(xiàn)象時(shí)提出來(lái)的。元胞自動(dòng)機(jī)是一組規(guī)則排布在元胞空間上且具有有限多個(gè)狀態(tài)量的元胞根據(jù)自身及其鄰居上一時(shí)刻的狀態(tài)量在離散的時(shí)間步上按照指定的演變規(guī)則進(jìn)行狀態(tài)更新的系統(tǒng)。在某一時(shí)刻，元胞的狀態(tài)量為其有限多個(gè)狀態(tài)量中的一種，例如閥門(mén)在某時(shí)刻非“關(guān)”即“開(kāi)”，該時(shí)刻任意元胞的狀態(tài)可表示為：

式中：s—元胞狀態(tài)；i—第i個(gè)元胞；t—時(shí)間步；N—鄰居。

元胞自動(dòng)機(jī)根據(jù)元胞的特點(diǎn)建立演變規(guī)則，能更真實(shí)地反映問(wèn)題的本質(zhì)。目前，元胞自動(dòng)機(jī)已廣泛應(yīng)用于應(yīng)急疏散[8]、熱傳導(dǎo)[9]。文獻(xiàn)[10]將元胞自動(dòng)機(jī)思想與有限元法結(jié)合形成了元胞單元法。元胞單元法將節(jié)點(diǎn)及單元同時(shí)考慮成元胞，并將節(jié)點(diǎn)和位移均考慮成元胞的狀態(tài)，元胞狀態(tài)的更新規(guī)則按照力與位移相互轉(zhuǎn)化來(lái)確定。提出的多邊形鄰域元胞自動(dòng)機(jī)算法，將元胞自動(dòng)機(jī)思想應(yīng)用于無(wú)網(wǎng)格法，形成一種簡(jiǎn)單的無(wú)網(wǎng)格法。該方法將問(wèn)題域離散成隨機(jī)分布的場(chǎng)節(jié)點(diǎn)，僅將位移量考慮成元胞的狀態(tài)量，借用有限元思想建立元胞與鄰居間的演變規(guī)則，用并行運(yùn)算及遺傳因素提升計(jì)算效率，具有求解大變形、網(wǎng)格畸變問(wèn)題的應(yīng)用前景。

2 算法構(gòu)造

二維彈性力學(xué)問(wèn)題在用多邊形鄰域元胞自動(dòng)機(jī)算法求解時(shí)，將問(wèn)題域離散成隨機(jī)分布的場(chǎng)節(jié)點(diǎn)，將任意一個(gè)離散的場(chǎng)節(jié)點(diǎn)考慮成一個(gè)元胞（Cell），將待求解的問(wèn)題域考慮成元胞空間，如圖1所示。

圖1 離散模型Fig.1 Discrete Model

元胞僅會(huì)“感知”與其鄰近的元胞及直接作用在該元胞上的外力。元胞的狀態(tài)量為矢量，分別為x、y方向的狀態(tài)量u、v。元胞對(duì)鄰近元胞的“感知”用鄰居元胞狀態(tài)量的加權(quán)求和近似表示為式（1）。元胞的狀態(tài)根據(jù)式（1）指定的演變規(guī)則在離散時(shí)間步上演化，直到元胞空間內(nèi)的元胞都達(dá)到穩(wěn)定狀態(tài)。

式中：fix、fiy—第 i個(gè)元胞 x、y 方向上的節(jié)點(diǎn)力；wi、ωi—第 i個(gè)元胞x、y方向節(jié)點(diǎn)力權(quán)重；φ、ψ—鄰居元胞狀態(tài)量權(quán)重；n—包含第i個(gè)元胞在內(nèi)的鄰居個(gè)數(shù)。

根據(jù)模型中元胞的分布情況指定合適的“感知”半徑，確定出元胞的鄰居。因元胞的分布情況不同，每個(gè)元胞的鄰居個(gè)數(shù)可能存在差異。根據(jù)鄰居個(gè)數(shù)的不同，建立邊數(shù)不同的正多邊形鄰域。正四邊形鄰域與正六邊形鄰域，如圖2所示。

圖2 正多邊形鄰域Fig.2 Regular Polygon Influence Zone

借用有限元思想建立式（1）的演變規(guī)則。以正四邊形為例說(shuō)明演變規(guī)則的構(gòu)造過(guò)程：

（1）將正四邊形鄰域考慮成一個(gè)局部有限元模型，借助有限元求解方程Ka=P建立中心元胞i與頂點(diǎn)A、B、C、D間的關(guān)系。只考慮鄰居元胞對(duì)中心元胞的影響，選取有限元求解方程中以中心點(diǎn)i為對(duì)角元的兩個(gè)線性方程：

式中：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—正多邊形頂點(diǎn) x 方向的狀態(tài)量；vp=［vA，vB，vC，vD］T—正多邊形頂點(diǎn)y方向的狀態(tài)量。

（2）分別在正多邊形中的每個(gè)三角形內(nèi)插值出對(duì)應(yīng)鄰居元胞的狀態(tài)量，如式（3）表示。每個(gè)鄰居點(diǎn)x方向的狀態(tài)量，y方向狀態(tài)量只需將式（3）中的u替換成v。

利用各三角形內(nèi)插值函數(shù)和為1，并將式（3）改寫(xiě)成矩陣形式：

式中：us=［u1，u2，u3，u4］T—鄰居元胞 x 方向的狀態(tài)量；vs=［v1，v2，v3，v4］T—鄰居元胞 y 方向的狀態(tài)量，E=［1，1，1，1］T。

（3）用式（4）求出正多邊形鄰域頂點(diǎn)的狀態(tài)量由鄰居元胞狀態(tài)量表示：

（4）式（5）代入到式（2）中簡(jiǎn)化得到鄰居元胞對(duì)中心元胞的影響：

式（6）表示的演化規(guī)則中，并未考慮元胞自身上一迭代步對(duì)應(yīng)方向上的影響。將其右側(cè)部分考慮成非自身因素的影響，用、表示。對(duì)（6）式加入自身因素的影響，可稱之為遺傳因素，如式（7）所示。

式中：α—遺傳系數(shù)，該系數(shù)對(duì)演變步數(shù)有較大影響。

圖2中1、2、3、4元胞與正四邊形鄰域頂點(diǎn)重合時(shí)：

簡(jiǎn)化式（7）得到該特殊情況下鄰居元胞對(duì)中心元胞影響的表達(dá)式，就是式（3）中直接將多邊形頂點(diǎn)狀態(tài)量改寫(xiě)成鄰居元胞的狀態(tài)量：

式（9）說(shuō)明當(dāng)鄰居元胞與多邊形鄰域頂點(diǎn)重合時(shí)，本算法的演變規(guī)則直接由局部有限元總剛獲得。當(dāng)離散域中每個(gè)元胞都選擇有限元網(wǎng)格劃分后與之相聯(lián)結(jié)的單元上的元胞作為鄰居，多邊形鄰域元胞自動(dòng)機(jī)算法將完全退化成有限元。

傳統(tǒng)無(wú)網(wǎng)格法在施加本質(zhì)邊界條件時(shí)，通常用特殊的處理辦法，例如罰函數(shù)法。本算法借助有限元插值思想建立演變規(guī)則，有限元是本算法的特殊情況，兩種算法能自然耦合。求解問(wèn)題時(shí)，本算法求解域邊界上元胞的鄰居與多邊形鄰域頂點(diǎn)重合，直接施加本質(zhì)邊界條件。邊界層采用有限元而模型內(nèi)部采用多邊形鄰域元胞自動(dòng)機(jī)算法會(huì)形成的過(guò)渡區(qū)用元胞自動(dòng)機(jī)算法。

該算法求解二維彈性力學(xué)問(wèn)題時(shí)的步驟為：

（1）將問(wèn)題域離散成隨機(jī)分布且總數(shù)為NCA的元胞，給出材料信息、邊界條件、誤差ε、初始元胞狀態(tài)u=0且v=0、遺傳因素α、迭代步 t=0；

（2）指定元胞“感知”半徑，搜索鄰居元胞，確定多邊形鄰域，計(jì)算式（7）中各權(quán)重系數(shù)；

（3）按隨機(jī)序列更新元胞狀態(tài)；按照式（7）的演變規(guī)則先更新第i個(gè)元胞x方向的狀態(tài)量ui，若指定了該點(diǎn)x方向的邊界位移為u0，則令ui=u0；然后以相同方式更新vi；若NCA個(gè)元胞都更新完畢，則迭代步t=t+1，進(jìn)行下一步；

3 算例

3.1 帶孔方板

有限大的四方受壓帶孔方板，如圖3所示。其解析解可用無(wú)限大（即 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，元胞總數(shù)1721個(gè)。虛線表示模型初始位置，位移放大了10倍。

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

模型求解是對(duì)求解域中的每個(gè)元胞進(jìn)行狀態(tài)更新。求解是一個(gè)動(dòng)態(tài)演化的過(guò)程。選取x=y=l處的點(diǎn)，作出了該點(diǎn)y向的狀態(tài)量的演變趨勢(shì)，如圖5所示。隨著迭代步數(shù)的增加，該點(diǎn)的狀態(tài)量逐漸趨于穩(wěn)定。經(jīng)驗(yàn)證，在-1＜α＜1內(nèi)，隨著α的減小，模型演化所需的求解步數(shù)減小。

顯示CA數(shù)值結(jié)果與有限元結(jié)果基本一致，而與解析解存在較大的誤差。原因是解析解是l＞＞a帶孔方板的解，如圖6所示。

圖5 位移變化趨勢(shì)Fig.5 Variation Trend of Displacement

圖6 帶孔方板位移結(jié)果Fig.6 Displacement of Pressurized Perforated Square Plate

3.2 懸臂梁

懸臂梁，如圖7所示。自由端x=l處施加拋物線型表面力為梁截面的剪應(yīng)力在該處的取值

圖7 懸臂梁Fig.7 Cantilever Beam

該問(wèn)題存在解析解：

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

建立懸臂梁離散模型，如圖8所示。載荷p=-1000N/m2，楊氏模量 E=2×105Pa，泊松比 υ=0.3，模型長(zhǎng) l=4.0m，模型寬 d=1.0m，板厚1.0m，α=-0.42。問(wèn)題域離散成分布相對(duì)均勻且總數(shù)為1874個(gè)元胞，由于元胞分布相對(duì)均勻，因此大多數(shù)元胞確定為了正六邊形鄰域。位移縮小為0.25。

圖8 懸臂梁模型Fig.8 Model of Cantilever Beam

懸臂梁同樣具有圖5所示的變化趨勢(shì)。該算法與有限元、解析解結(jié)果基本吻合，如圖9所示。

圖9 懸臂梁位移結(jié)果Fig.9 Displacement of Cantilever Beam

3.3 并行計(jì)算

多邊形鄰域元胞自動(dòng)機(jī)算法求解問(wèn)題時(shí)，僅需存儲(chǔ)鄰居元胞的相關(guān)信息，不需要合成總剛。這較有限元及傳統(tǒng)的無(wú)網(wǎng)格法節(jié)約了存儲(chǔ)空間。本算法雖節(jié)省了存儲(chǔ)空間，但用迭代時(shí)間步求解，存在迭代步數(shù)較多，演化時(shí)間長(zhǎng)的問(wèn)題。本算法元胞的演變僅與其鄰近元胞的狀態(tài)相關(guān)，元胞的演變具有高度的并行性。按照并行運(yùn)算思路，將前述元胞狀態(tài)更新過(guò)程修改如下：

（1）離散問(wèn)題域，設(shè)置初始條件；（2）選擇鄰居，計(jì)算式（7）中的權(quán)重；（3）將NCA個(gè)元胞隨機(jī)分配到m個(gè)處理器中，每個(gè)處理器同時(shí)按式（7）更新分配到的元胞的狀態(tài)；若所有處理器中的元胞狀態(tài)都更新完畢，則迭代步t=t+1；（4）若滿足收斂條件，則求解結(jié)束；否則，重復(fù)第（3）步。

用單核、雙核及四核更新元胞數(shù)為7267個(gè)的懸臂梁模型，雙核求解用時(shí)約為單核求解用時(shí)的50%，如圖10所示。四核求解用時(shí)約為單核用時(shí)的25%。

圖10 并行計(jì)算Fig.10 Parallel Computation

4 結(jié)論

（1）多邊形鄰域元胞自動(dòng)機(jī)算法用有限元插值思想建立元胞間的演變規(guī)則，演變規(guī)則構(gòu)建簡(jiǎn)單計(jì)算量小且能與有限元自然耦合，本質(zhì)邊界條件施加簡(jiǎn)單。

（2）本算法求解二維彈性力學(xué)問(wèn)題時(shí)，將問(wèn)題域離散成隨機(jī)分布的元胞，用元胞狀態(tài)的演變來(lái)獲得穩(wěn)定的數(shù)值解。算例表明，該算法在元胞分布相對(duì)均勻時(shí)能獲得收斂且誤差較小的數(shù)值解。

（3）逐個(gè)更新元胞的狀態(tài)量，整體問(wèn)題求解時(shí)間較長(zhǎng)。根據(jù)元胞自動(dòng)機(jī)高度并行的特點(diǎn)，采用并行運(yùn)算大大提升了該算法的計(jì)算效率?？紤]遺傳因素影響，同樣提升了計(jì)算效率。

（4）該算法求解問(wèn)題時(shí)不需要網(wǎng)格，是一種簡(jiǎn)單的無(wú)網(wǎng)格法。算法在求解大變形、網(wǎng)格畸變等等領(lǐng)域具有應(yīng)用前景。