陳澤蕓 袁衛(wèi)鋒

西南科技大學(xué)制造科學(xué)與工程學(xué)院，綿陽,621010

彈性力學(xué)中無網(wǎng)格和有限元耦合的元胞自動(dòng)機(jī)算法

陳澤蕓 袁衛(wèi)鋒

西南科技大學(xué)制造科學(xué)與工程學(xué)院，綿陽,621010

結(jié)合有限元和無網(wǎng)格算法的優(yōu)勢(shì)，提出了一種元胞自動(dòng)機(jī)算法用以求解二維彈性力學(xué)問題。該算法將二維模型離散成一系列節(jié)點(diǎn)，這些節(jié)點(diǎn)被分成有限元群和無網(wǎng)格群。有限元區(qū)域被定義在問題的邊界附近，其中的任一節(jié)點(diǎn)和其周圍相鄰點(diǎn)的力學(xué)關(guān)系通過有限元單元建立；無網(wǎng)格區(qū)域定義在遠(yuǎn)離原理問題邊界處，其中的節(jié)點(diǎn)之間的關(guān)系借用有限元中的位移插值概念建立。無論處于有限元區(qū)域還是無網(wǎng)格區(qū)域，任何一個(gè)節(jié)點(diǎn)都被置于元胞自動(dòng)機(jī)的框架下進(jìn)行處理，即節(jié)點(diǎn)的位移通過元胞自動(dòng)機(jī)進(jìn)行求解。與有限元方法相比，所提出的元胞自動(dòng)機(jī)算法無需采用高斯消去法等傳統(tǒng)系統(tǒng)求解器，而是通過元胞自動(dòng)機(jī)的自動(dòng)演化解決問題。依據(jù)該算法，有限元和無網(wǎng)格方法可以實(shí)現(xiàn)無縫連接。數(shù)值算例驗(yàn)證了該算法的新穎性和正確性。

耦合算法；元胞自動(dòng)機(jī)；有限元方法；彈性力學(xué)

0 引言

有限元方法的基本思想是將連續(xù)求解區(qū)域離散成有限個(gè)按一定方式相互連接在一起的單元的組合體。有限元方法利用每一個(gè)單元內(nèi)假設(shè)的近似函數(shù)來分片地表示全求解域上待求的未知場函數(shù)[1]。但在求解某些特殊問題，如計(jì)算流固耦合、動(dòng)態(tài)裂紋擴(kuò)展等問題時(shí)，有限元網(wǎng)格可能會(huì)產(chǎn)生畸變，這種情況不僅需要不斷地進(jìn)行網(wǎng)格重構(gòu)，而且會(huì)嚴(yán)重影響解的精度，這就使得有限元方法在工程應(yīng)用中有一定的局限性。為解決有限元方法在計(jì)算中由于網(wǎng)格畸變引起的缺陷，近年來無網(wǎng)格方法成為研究熱點(diǎn)。與有限元方法不同，無網(wǎng)格法是用一組點(diǎn)來離散求解區(qū)域，直接借助于離散點(diǎn)來構(gòu)造近似函數(shù)，可以徹底或部分消除網(wǎng)格，不需要網(wǎng)格的初始劃分和重構(gòu)[2]，不僅可以確保計(jì)算的精度，而且可以降低計(jì)算的難度。現(xiàn)在無網(wǎng)格法已經(jīng)成功應(yīng)用于固體力學(xué)彈塑性分析[3]以及高速碰撞[4]、斷裂損傷力學(xué)等問題的研究中[5]，但無網(wǎng)格方法也存在基本邊界條件不能直接進(jìn)行施加[6]、計(jì)算量大等問題。將無網(wǎng)格法和有限元法耦合是無網(wǎng)格法處理邊界條件的方法之一。

元胞自動(dòng)機(jī)(cellular automata，CA)是一種時(shí)間、空間、狀態(tài)都離散的動(dòng)力系統(tǒng)模型，不同于一般的動(dòng)力學(xué)模型，它不是由嚴(yán)格定義的物理方程或函數(shù)確定的，而是由一系列模型構(gòu)造的規(guī)則構(gòu)成，凡是滿足這些規(guī)則的模型都可以算作CA模型。元胞自動(dòng)機(jī)在熱擴(kuò)散[7]、并行計(jì)算[8]、緊急疏散[9]、社會(huì)交通[10]和復(fù)雜的社會(huì)經(jīng)濟(jì)現(xiàn)象[11]等領(lǐng)域都有廣泛的應(yīng)用。

本文提出了一種將無網(wǎng)格元胞自動(dòng)機(jī)方法與有限元方法相結(jié)合的耦合算法，并將其應(yīng)用到二維彈性力學(xué)問題的求解中。

1 二維彈性力學(xué)的耦合模型

1.1方法描述

用一組隨機(jī)離散點(diǎn)對(duì)二維彈性力學(xué)模型進(jìn)行離散，如圖1所示。模型被分成兩個(gè)部分，靠近邊界的區(qū)域?yàn)橛邢拊獏^(qū)域，中間部分為元胞自動(dòng)機(jī)區(qū)域。有限元區(qū)域用有限元方法進(jìn)行網(wǎng)格劃分，利用有限元三角形單元插值函數(shù)建立這一區(qū)域內(nèi)節(jié)點(diǎn)間的相互關(guān)系。元胞自動(dòng)機(jī)區(qū)域因?yàn)闆]有網(wǎng)格劃分，所以需要先確定其節(jié)點(diǎn)的影響域，從影響域內(nèi)選擇與其相關(guān)的節(jié)點(diǎn)，再建立節(jié)點(diǎn)間的相互關(guān)系，具體步驟如下。

圖1 離散模型示意圖Fig.1 Discretization of the model

(1)確定節(jié)點(diǎn)影響域的形狀和大小。本文選擇的節(jié)點(diǎn)影響域形狀為正六邊形，如圖2所示。將正六邊形影響域分成6個(gè)大小完全相等的正三角形，定義圖2中一象限內(nèi)正三角形邊與x軸構(gòu)成的最小夾角θ為正六邊形影響域的角度。影響域的大小由模型離散的疏密程度決定，只需保證影響域每個(gè)三角形區(qū)域內(nèi)至少有一個(gè)點(diǎn)即可。影響域角度由模型離散點(diǎn)的位置決定，所以角度θ可能是0～60°內(nèi)的任一角度。

圖2 節(jié)點(diǎn)影響域的大小和位置Fig.2 The hexagon influence zone

(2)選擇影響域中與當(dāng)前節(jié)點(diǎn)相關(guān)的節(jié)點(diǎn)。從節(jié)點(diǎn)的正六邊形影響域分成的6個(gè)正三角形區(qū)域中各選一個(gè)點(diǎn)作為當(dāng)前節(jié)點(diǎn)的相關(guān)節(jié)點(diǎn)，若一個(gè)三角形區(qū)域中不止一個(gè)節(jié)點(diǎn)，則從中隨機(jī)選出一個(gè)點(diǎn)作為當(dāng)前節(jié)點(diǎn)的相關(guān)節(jié)點(diǎn)。

(3)建立所選節(jié)點(diǎn)間的相互關(guān)系。利用有限元方法中的插值概念構(gòu)造元胞自動(dòng)機(jī)區(qū)域內(nèi)任一節(jié)點(diǎn)與其相關(guān)節(jié)點(diǎn)間的相互關(guān)系，具體方法如下。①對(duì)正六邊形影響域分成的每個(gè)三角形區(qū)域，用有限元方法中的三角形插值函數(shù)建立中間點(diǎn)O與影響域頂點(diǎn)A～F之間位移與力的關(guān)系：

(1)

式中,K為正六邊形影響域形成的剛度矩陣；u為中間點(diǎn)O和影響域頂點(diǎn)A～F的位移分量；p為中間點(diǎn)O和影響域頂點(diǎn)A～F所受力的分量；下標(biāo)V代表影響域頂點(diǎn)A～F。

取方程組中計(jì)算節(jié)點(diǎn)O位移的兩個(gè)方程，得

KOVuV+KOOuO=pO

(2)

②借助有限元三角形單元插值構(gòu)造離散節(jié)點(diǎn)O及節(jié)點(diǎn)1～6和點(diǎn)A～F之間位移與力的關(guān)系，如在三角形OAB中有

(3)

同理，借助有限元三角形單元插值可以構(gòu)造其他5個(gè)三角形單元的內(nèi)節(jié)點(diǎn)與三角形頂點(diǎn)的位移函數(shù)關(guān)系。將這6個(gè)三角形單元形成的方程組寫成矩陣形式為

(4)

式中，uR、uV分別為真實(shí)節(jié)點(diǎn)1～6和虛擬點(diǎn)A～F的位移分量。

uO=λpO+γuR

(5)

至此，式(5)為所構(gòu)造的元胞自動(dòng)機(jī)區(qū)域內(nèi)任一節(jié)點(diǎn)與其影響域內(nèi)相關(guān)節(jié)點(diǎn)之間位移和力的關(guān)系，即CA的局部規(guī)則。

1.2與有限元法耦合

需要指出的是，元胞自動(dòng)機(jī)區(qū)域內(nèi)節(jié)點(diǎn)的影響域可能包含有限元部分的節(jié)點(diǎn)，其鄰居也會(huì)包含有限元部分內(nèi)的節(jié)點(diǎn)。兩個(gè)區(qū)域交界面上的點(diǎn)也用CA方法選擇其鄰居并建立交界面上點(diǎn)與其鄰居間的關(guān)系，交界面上的點(diǎn)必然會(huì)將有限元部分和元胞自動(dòng)機(jī)部分中的點(diǎn)選作鄰居。用上述方法建立節(jié)點(diǎn)間相互關(guān)系的同時(shí)建立了元胞自動(dòng)機(jī)部分和有限元部分的聯(lián)系，即無網(wǎng)格部分和網(wǎng)格部分用式(5)進(jìn)行自然連接，邊界上各物理量會(huì)通過式(5)傳遞到CA部分內(nèi)部。兩個(gè)部分間不需要增加過渡單元或過渡計(jì)算，算法簡單。

將有限元部分的求解方程組改寫成節(jié)點(diǎn)間相互關(guān)系：

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

(6)

式中,n為模型總的離散節(jié)點(diǎn)數(shù);ω為由正六邊形影響域剛度矩陣中的第x行系數(shù)；uR′為與節(jié)點(diǎn)x相關(guān)的節(jié)點(diǎn)位移分量;Kx,x為剛度矩陣K對(duì)角線元素。

由式(5)和式(6)整理可得到模型中所有節(jié)點(diǎn)與其相關(guān)節(jié)點(diǎn)間的關(guān)系，寫為

ux=apx-buR′

(7)

式中,a、b為常數(shù)，由相應(yīng)的三角形插值函數(shù)和有限元部分的剛度矩陣計(jì)算得到。

由此，式(7)作為整個(gè)模型的元胞自動(dòng)機(jī)自然演化局部規(guī)則。

1.3求解步驟

式(7)給出了模型中任意節(jié)點(diǎn)與其相關(guān)節(jié)點(diǎn)的關(guān)系，以式(7)作為元胞自動(dòng)機(jī)基本演化規(guī)則，離散模型總數(shù)n作為一次循環(huán)，從任意節(jié)點(diǎn)開始以任意順序進(jìn)行元胞自動(dòng)機(jī)演化求解。重復(fù)循環(huán)直到結(jié)果收斂，收斂條件是t+1時(shí)刻位移與t時(shí)刻位移的最大值小于給定誤差限，即

(8)

具體演化步驟為：①設(shè)置每個(gè)節(jié)點(diǎn)的位移初始值；②根據(jù)t時(shí)刻節(jié)點(diǎn)位移，計(jì)算t+1時(shí)刻的節(jié)點(diǎn)位移；③重復(fù)步驟②至結(jié)果滿足收斂條件。

因?yàn)槊總€(gè)時(shí)間步的循環(huán)中計(jì)算順序是任意的且每次順序都不同，這種計(jì)算順序的任意性為并行計(jì)算的實(shí)現(xiàn)提供有利條件。

1.4應(yīng)力計(jì)算

求得位移之后，用位移法求每個(gè)節(jié)點(diǎn)的應(yīng)力。有網(wǎng)格部分，根據(jù)已有網(wǎng)格劃分用有限元方法計(jì)算其單元應(yīng)力。元胞自動(dòng)機(jī)部分根據(jù)每個(gè)節(jié)點(diǎn)選出的6個(gè)相關(guān)節(jié)點(diǎn)依次選兩個(gè)點(diǎn)與當(dāng)前點(diǎn)構(gòu)成6個(gè)三角形單元，通過各節(jié)點(diǎn)的已知位移求單元應(yīng)力。例如，任意三角形OAB的單元應(yīng)力可寫為

(9)

式中,s為對(duì)應(yīng)節(jié)點(diǎn)應(yīng)力矩陣。

任一節(jié)點(diǎn)O的節(jié)點(diǎn)應(yīng)力為與節(jié)點(diǎn)相關(guān)的單元應(yīng)力的數(shù)值平均，即

(10)

式中，m為與節(jié)點(diǎn)O相關(guān)的單元個(gè)數(shù)。

2 算例

2.1算例1

四邊受均布?jí)毫Φ姆桨?，邊長為4 m，均布力大小q=1.0×1010Pa，如圖3所示。材料參數(shù)為彈性模量E=2.0×1011Pa，泊松比μ=0.3。因?yàn)槟Ｐ完P(guān)于x、y方向都對(duì)稱，所以取模型的1/4進(jìn)行計(jì)算。節(jié)點(diǎn)影響域的大小和角度取決于離散點(diǎn)的疏密程度，所以模型中各節(jié)點(diǎn)的影響域大小和角度均不相同，如圖4所示，這使得算法具有良好的適應(yīng)性。

圖3 四邊受壓方板模型Fig.3 Square plate subjected to distributed forces

圖4 節(jié)點(diǎn)影響域的大小和角度Fig.4 Sizes and orientations of the nodal influence zones

圖5 位移計(jì)算結(jié)果比較Fig.5 Displacements of the nodes on the boundary

模型總離散點(diǎn)數(shù)為1681，其中有限元區(qū)域節(jié)點(diǎn)數(shù)為312，交界面上節(jié)點(diǎn)數(shù)為144，CA區(qū)域節(jié)點(diǎn)數(shù)為1225。取模型下方y(tǒng)=2的41個(gè)節(jié)點(diǎn)位移計(jì)算結(jié)果與理論解比較，結(jié)果如圖5所示。x方向位移沿負(fù)方向線性增大，y方向位移均為-0.007。由圖5可以看出，理論解和數(shù)值解吻合得很好，最大相對(duì)誤差小于0.001%。

同樣的模型，有限元區(qū)域和交界面上離散點(diǎn)數(shù)不變，調(diào)整CA區(qū)域內(nèi)節(jié)點(diǎn)個(gè)數(shù)，從1225減小到841、529，模型計(jì)算結(jié)果相差不大。該結(jié)果說明，只要保證CA區(qū)域內(nèi)節(jié)點(diǎn)能夠與有限元區(qū)域建立聯(lián)系，則在一定范圍內(nèi)，CA區(qū)域內(nèi)節(jié)點(diǎn)的個(gè)數(shù)對(duì)計(jì)算精度影響不大。

2.2算例2

取與算例1模型一樣的方板，僅中間開一個(gè)半徑為0.5 m的孔。同樣取模型的1/4作為計(jì)算模型，模型總的離散點(diǎn)數(shù)為1677。取模型圓弧邊界和下邊界共51個(gè)節(jié)點(diǎn)的y方向應(yīng)力集中系數(shù)計(jì)算結(jié)果與理論結(jié)果、ANSYS模擬結(jié)果相比較，結(jié)果如圖6所示。從圖6可以看出本文算法的計(jì)算結(jié)果與理論解之間存在一定誤差，但與ANSYS的計(jì)算結(jié)果十分接近。算例2中的理論解是四邊受均布?jí)毫Φ臒o限大板中間帶小圓孔的理論解，而算例2中的方板模型是有限大的，所以求得的應(yīng)力數(shù)值解與理論解有一定的差異。而與ANSYS模擬結(jié)果十分接近則可以說明本文提出的耦合算法的正確性和可行性。

圖6 應(yīng)力集中系數(shù)計(jì)算結(jié)果比較Fig.6 Comparison of the stress concentration factors

3 結(jié)論

(1)將有限元方法與元胞自動(dòng)機(jī)相結(jié)合，提出了一種新的耦合算法。該算法利用插值概念建立任意離散點(diǎn)與其周圍節(jié)點(diǎn)之間的力和位移的關(guān)系，根據(jù)元胞自動(dòng)機(jī)的演化求解二維彈性力學(xué)問題。計(jì)算過程中只需要建立有限元和元胞自動(dòng)機(jī)的相應(yīng)演化規(guī)則，整體基于元胞自動(dòng)機(jī)的演化過程求解，兩個(gè)部分自然結(jié)合，不需要增加任何過渡計(jì)算。算例證明該算法簡單有效。

(2)雖然該算法需要一段時(shí)間更新計(jì)算結(jié)果，可能會(huì)降低計(jì)算效率，但是計(jì)算過程中任意計(jì)算順序體現(xiàn)了該算法在并行運(yùn)算中的潛力及優(yōu)勢(shì)，為進(jìn)一步提高計(jì)算效率提供了可能。

(3)目前該耦合算法只應(yīng)用于二維彈性力學(xué)，對(duì)于三維問題只需要將節(jié)點(diǎn)影響域變成正四面體或正八面體。在未來的探索中應(yīng)該突出有限元及元胞自動(dòng)機(jī)的特點(diǎn)和優(yōu)勢(shì)，建立更充分的理論計(jì)算基礎(chǔ)，以推廣其應(yīng)用到不連續(xù)問題和大變形的計(jì)算中。

[1] 王勖成,邵敏.有限單元法基本原理和數(shù)值方法[M].北京:清華大學(xué)出版社,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 程序及其在斜高速碰撞問題的應(yīng)用[J] .計(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ìn)雄,李梅娥,張紅艷，等.再生核質(zhì)點(diǎn)法研究進(jìn)展[J].力學(xué)進(jìn)展,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]. 計(jì)算力學(xué)學(xué)報(bào),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] 鄭亮, 馬壽峰, 賈寧. 基于駕駛員行為的元胞自動(dòng)機(jī)模型研究[J]. 物理學(xué)報(bào), 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] 顏偉, 唐德善. 區(qū)域經(jīng)濟(jì)聚集效應(yīng)的元胞自動(dòng)機(jī)模擬[J]. 商場現(xiàn)代化, 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

國家自然科學(xué)基金資助項(xiàng)目(11472232)；四川省教育廳教改項(xiàng)目(14JGCX07);制造過程測(cè)試技術(shù)省部共建實(shí)驗(yàn)室基金資助項(xiàng)目(13zxzk05)

TB121

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

陳澤蕓，女，1991年生。西南科技大學(xué)制造科學(xué)與工程學(xué)院碩士研究生。發(fā)表論文5篇。E-mail:464252016@qq.com。袁衛(wèi)鋒，男，1970年生。西南科技大學(xué)制造科學(xué)與工程學(xué)院研究員。