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

陳澤蕓 袁衛(wèi)鋒

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

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

陳澤蕓 袁衛(wèi)鋒

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

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

耦合算法；元胞自動機；有限元方法；彈性力學

0 引言

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

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

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

1 二維彈性力學的耦合模型

1.1方法描述

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

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

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

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

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

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

(1)

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

取方程組中計算節(jié)點O位移的兩個方程，得

KOVuV+KOOuO=pO

(2)

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

(3)

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

(4)

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

uO=λpO+γuR

(5)

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

1.2與有限元法耦合

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

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

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

(6)

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

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

ux=apx-buR′

(7)

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

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

1.3求解步驟

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

(8)

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

因為每個時間步的循環(huán)中計算順序是任意的且每次順序都不同，這種計算順序的任意性為并行計算的實現(xiàn)提供有利條件。

1.4應力計算

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

(9)

式中,s為對應節(jié)點應力矩陣。

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

(10)

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

2 算例

2.1算例1

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

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

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

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

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

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

2.2算例2

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

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

3 結(jié)論

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

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

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

[1] 王勖成,邵敏.有限單元法基本原理和數(shù)值方法[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] 周進雄,李梅娥,張紅艷，等.再生核質(zhì)點法研究進展[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] 顏偉, 唐德善. 區(qū)域經(jīng)濟聚集效應的元胞自動機模擬[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

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

TB121

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

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