自適應(yīng)一致性高階無(wú)單元伽遼金法1)

邵玉龍段慶林,?,,2)高 欣李錫夔張洪武

?(大連理工大學(xué)工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，大連116024)

?(成都理工大學(xué)地質(zhì)災(zāi)害防治與地質(zhì)環(huán)境保護(hù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，成都610059)

??(武漢大學(xué)水資源與水電工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢430072)

自適應(yīng)一致性高階無(wú)單元伽遼金法1)

邵玉龍?段慶林?,?,??,2)高 欣?李錫夔?張洪武?

?(大連理工大學(xué)工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，大連116024)

?(成都理工大學(xué)地質(zhì)災(zāi)害防治與地質(zhì)環(huán)境保護(hù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，成都610059)

??(武漢大學(xué)水資源與水電工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢430072)

近來(lái)提出的一致性高階無(wú)單元伽遼金法通過(guò)導(dǎo)數(shù)修正技術(shù)大幅度減少了所需積分點(diǎn)數(shù)目，并能夠精確地通過(guò)線性和二次分片試驗(yàn)，顯著改善標(biāo)準(zhǔn)無(wú)單元伽遼金法的計(jì)算效率、精度和收斂性.本文在此基礎(chǔ)之上，充分利用無(wú)單元法易于在局部區(qū)域添加節(jié)點(diǎn)的優(yōu)勢(shì)，發(fā)展了一致性高階無(wú)單元伽遼金法的h型自適應(yīng)分析方法.根據(jù)應(yīng)變能密度梯度該方法自適應(yīng)地確定需節(jié)點(diǎn)加密的區(qū)域，基于背景積分網(wǎng)格的局部多層細(xì)化要求生成新的計(jì)算節(jié)點(diǎn)，同時(shí)考慮了節(jié)點(diǎn)分布由密到疏漸進(jìn)過(guò)渡的情形.采用相鄰兩次計(jì)算的應(yīng)變能的相對(duì)誤差作為自適應(yīng)過(guò)程的停止準(zhǔn)則,將所發(fā)展自適應(yīng)無(wú)網(wǎng)格法應(yīng)用于由幾何外形、邊界外載和體力等因素造成的應(yīng)力集中問(wèn)題的計(jì)算分析.數(shù)值結(jié)果表明，所發(fā)展方法能夠自適應(yīng)地對(duì)高應(yīng)力梯度區(qū)域進(jìn)行節(jié)點(diǎn)加密，自動(dòng)給出合理的計(jì)算節(jié)點(diǎn)分布.與已有的標(biāo)準(zhǔn)無(wú)網(wǎng)格法的自適應(yīng)分析相比，所發(fā)展方法在計(jì)算效率、精度和應(yīng)力場(chǎng)光滑性等方面均展現(xiàn)出顯著優(yōu)勢(shì).與采用節(jié)點(diǎn)均勻分布的一致性高階無(wú)單元伽遼金法相比，它大幅度地減少了計(jì)算節(jié)點(diǎn)數(shù)目，有效提高了一致性高階無(wú)單元伽遼金法在分析應(yīng)力集中等存在局部高梯度問(wèn)題時(shí)的計(jì)算效率和求解精度.

無(wú)單元伽遼金法，無(wú)網(wǎng)格法，自適應(yīng)分析，應(yīng)變能密度，應(yīng)力集中

引言

無(wú)單元伽遼金(element-free Galerkin,EFG)法作為一種具代表性的無(wú)網(wǎng)格方法，由于具有易于建立高階近似函數(shù)、易于處理大變形和實(shí)現(xiàn)自適應(yīng)計(jì)算等優(yōu)勢(shì)而引發(fā)計(jì)算力學(xué)界的廣泛關(guān)注，已被成功應(yīng)用于金屬成型[1-2]、裂紋擴(kuò)展[3-4]以及板彎曲[5]等眾多問(wèn)題的數(shù)值分析與模擬中.然而，EFG方法的計(jì)算效率不高這一缺陷仍然阻礙著該方法在實(shí)際工程問(wèn)題中的廣泛應(yīng)用.導(dǎo)致其效率低下的一個(gè)主要原因是其需要較多的數(shù)值積分點(diǎn).例如，采用背景三角形積分網(wǎng)格的二階EFG方法通常要求在每個(gè)三角形積分子域內(nèi)至少使用16個(gè)積分點(diǎn)，這嚴(yán)重降低了它的計(jì)算效率.近來(lái)很多學(xué)者對(duì)此進(jìn)行了改進(jìn)[6-9]，其中由段慶林等[8-9]提出的一致性無(wú)單元伽遼金(consistent EFG,CEFG)法通過(guò)導(dǎo)數(shù)修正技術(shù)大幅度減少了所需的積分點(diǎn)數(shù)目.例如，二階的CEFG方法在每個(gè)背景三角形積分子域內(nèi)僅需3個(gè)積分點(diǎn).更重要的是，二階CEFG方法能夠精確地通過(guò)線性和二次分片試驗(yàn)，因而同時(shí)具有很高的計(jì)算精度和效率，已在熱傳導(dǎo)[10]、動(dòng)力響應(yīng)[11]、三維彈塑性[12]以及不可壓縮固體的變形分析[13-14]等多種問(wèn)題中獲得成功應(yīng)用，展現(xiàn)出良好的發(fā)展?jié)摿?

本文目的是在CEFG方法的基礎(chǔ)之上，充分利用無(wú)網(wǎng)格法便于在局部區(qū)域添加計(jì)算節(jié)點(diǎn)的優(yōu)勢(shì)，為CEFG發(fā)展自適應(yīng)分析方法，進(jìn)一步提高它求解局部高梯度問(wèn)題的計(jì)算效率.

自適應(yīng)方法能夠根據(jù)解的性態(tài)對(duì)網(wǎng)格或節(jié)點(diǎn)進(jìn)行局部自動(dòng)加密，可大幅度減少節(jié)點(diǎn)數(shù)目和減小計(jì)算規(guī)模，通過(guò)較小的計(jì)算量便能得到高精度解，因而在現(xiàn)代工程計(jì)算中得到了廣泛應(yīng)用[15-17].無(wú)網(wǎng)格法(如EFG法)近似函數(shù)的建立不依賴于網(wǎng)格單元，因而在建立節(jié)點(diǎn)形函數(shù)的意義上無(wú)須將新生成節(jié)點(diǎn)與周圍節(jié)點(diǎn)連成網(wǎng)格單元，這為局部加密提供了很大便利.因而，不少學(xué)者在無(wú)網(wǎng)格法的自適應(yīng)分析方法方面開(kāi)展了大量的研究工作.Duarte和Oden[18]首先提出了h,p和h-p型自適應(yīng)Hp無(wú)單元云團(tuán)法.劉欣等[19]在此基礎(chǔ)上，針對(duì)平面彈性問(wèn)題發(fā)展了一種顯示后驗(yàn)誤差指標(biāo)，對(duì)平面裂紋進(jìn)行了自適應(yīng)分析.數(shù)值結(jié)果表明，該自適應(yīng)方案能夠在高應(yīng)力梯度的裂紋尖端進(jìn)行自適應(yīng)網(wǎng)格加密，并具有很好的收斂性和精度.Chung和Belytschko[20]采用縮減影響域的形函數(shù)，在EFG方法的框架下提出了一種簡(jiǎn)單且易于實(shí)現(xiàn)的誤差分析方法，并成為后續(xù)多種自適應(yīng)無(wú)單元方法[21-22]的基礎(chǔ).Hiiussler-Combe和Korn[23]以及Rabczuk和Belytschko[24]將基于Taylor展開(kāi)的誤差估計(jì)應(yīng)用于線彈性斷裂的自適應(yīng)分析，發(fā)現(xiàn)基于該誤差分析的自適應(yīng)方法收斂性好，并能夠?qū)Ω邞?yīng)力梯度區(qū)域進(jìn)行網(wǎng)格加密.Liu和Tu[25]采用不同數(shù)量的積分點(diǎn)對(duì)能量進(jìn)行數(shù)值積分，將它們的差作為評(píng)價(jià)背景網(wǎng)格是否需要細(xì)化的標(biāo)準(zhǔn)，提出了基于背景網(wǎng)格的自適應(yīng)EFG方法.干年妃等[26]對(duì)該方法進(jìn)行了改進(jìn)，發(fā)展了基于RKPM方法的自適應(yīng)無(wú)網(wǎng)格法，對(duì)撞擊、擠壓成型等具有局部高梯度的問(wèn)題進(jìn)行了成功的數(shù)值模擬.Luo和Combe[27]使用平均應(yīng)變能梯度作為網(wǎng)格加密準(zhǔn)則，提出了基于應(yīng)變能密度梯度的自適應(yīng)EFG方法.張征等[28]將基于應(yīng)變能的自適應(yīng)準(zhǔn)則應(yīng)用于接觸問(wèn)題，相對(duì)于均勻網(wǎng)格大幅度提高了分析效率.由于重新構(gòu)造形函數(shù)比較耗時(shí)，Metsis等[29]提出了一個(gè)兩級(jí)h型自適應(yīng)無(wú)單元方法，僅需要計(jì)算所添加節(jié)點(diǎn)對(duì)剛度陣的貢獻(xiàn)，不需要重新計(jì)算初始剛度陣.

應(yīng)指出的是，盡管自適應(yīng)分析的研究已有如前所述的大量工作，但它們均基于標(biāo)準(zhǔn)的無(wú)網(wǎng)格法(如EFG方法)，在每個(gè)積分子域內(nèi)仍然須使用大量積分點(diǎn)，如Chung和Belytschko[20]在每個(gè)積分單元上使用了5×5的高斯積分點(diǎn)，Sethuraman等[21]則更是使用了6×6的積分點(diǎn)配置，這仍然嚴(yán)重影響了計(jì)算效率.與此不同的是，本文發(fā)展的自適應(yīng)方法是基于已大幅度減少積分點(diǎn)且能精確通過(guò)分片試驗(yàn)的一致性無(wú)單元伽遼金法,將能展現(xiàn)出更好的計(jì)算效率和精度.

1 一致性無(wú)單元伽遼金法：積分格式與導(dǎo)數(shù)修正

考慮區(qū)域?yàn)棣?、邊界為Γ的二維彈性體，平衡方程和相應(yīng)的邊界條件為

其中，σ為Cauchy應(yīng)力，b為體力，n為邊界的單位外法線矢量，分別為邊界上的固定力和位移.彈性本構(gòu)關(guān)系可寫(xiě)為

其中ε為應(yīng)變，D為彈性模量矩陣.

采用無(wú)單元伽遼金法，位移近似為

其中，uI為節(jié)點(diǎn)位移參數(shù)向量.

為形函數(shù)矩陣，NI(x)為節(jié)點(diǎn)的無(wú)網(wǎng)格形函數(shù)，它的計(jì)算過(guò)程可參見(jiàn)文獻(xiàn)[30].采用Nitsche法[31]施加固定位移邊界條件，最終的伽遼金離散方程為其中K為剛度陣，f為等效節(jié)點(diǎn)載荷向量，β為懲罰系數(shù)，KΓu,Kp,fΓu和fp為Nitsche法引入的額外項(xiàng)，具體的推導(dǎo)過(guò)程參見(jiàn)文獻(xiàn)[8].

由于EFG方法的節(jié)點(diǎn)形函數(shù)是非多項(xiàng)式的有理函數(shù)，為保證計(jì)算的穩(wěn)定性和精度，需要較多的數(shù)值積分點(diǎn)，這嚴(yán)重降低了計(jì)算效率.針對(duì)該問(wèn)題，本文采用Duan等[8]提出的一致性無(wú)單元伽遼金法.對(duì)于本文考慮的二階近似，該方法采用如圖1所示的二階一致三點(diǎn)積分格式 (quadratically consistent 3-point integration,QC3).圖1中的黑色實(shí)心點(diǎn)表示CEFG法的計(jì)算節(jié)點(diǎn)，將它們連成背景積分子域，每個(gè)子域采用3個(gè)域內(nèi)積分點(diǎn)(紫五角星)，在每條邊上使用兩個(gè)一維高斯點(diǎn)(紅花點(diǎn))用于計(jì)算邊界積分.

圖1 QC3積分方法示意圖Fig.1 Schematic diagram of the QC3 integration method

域積分點(diǎn)上形函數(shù)的空間導(dǎo)數(shù)由如下它與形函數(shù)之間的散度定理確定

其中，ΩS是背景積分子域，ΓS為積分子域的邊界，是MLS形函數(shù)的基底函數(shù)向量.本文采用二次近似，即

應(yīng)用圖1積分格式對(duì)式(7)進(jìn)行數(shù)值積分得到

求解方程式(11)，可得到3個(gè)積分點(diǎn)上的x方向的修正導(dǎo)數(shù)NI,x(x)，y方向修正導(dǎo)數(shù)NI,y(x)同理可求得.這些導(dǎo)數(shù)將被用于計(jì)算剛度陣K.

2 自適應(yīng)方案

2.1 細(xì)化區(qū)域的確定

本文采用Luo和Hiiussler-Combe[27]提出的基于應(yīng)變能密度梯度的準(zhǔn)則來(lái)確定需節(jié)點(diǎn)加密的區(qū)域，其中應(yīng)變能密度為

定義應(yīng)變能密度的平均梯度為

顯然，GSED能夠識(shí)別出應(yīng)變能密度變化劇烈的區(qū)域，如應(yīng)力集中區(qū)域.然而，在識(shí)別出的高梯度區(qū)域上是否需要進(jìn)一步細(xì)分網(wǎng)格(加密節(jié)點(diǎn))，還需考慮該區(qū)域當(dāng)前的網(wǎng)格分布情況.為此，定義網(wǎng)格密度為

其中，N為面積為A的區(qū)域內(nèi)的節(jié)點(diǎn)數(shù).為進(jìn)一步度量當(dāng)前計(jì)算網(wǎng)格關(guān)于應(yīng)變能密度梯度的密集程度，引入網(wǎng)格強(qiáng)度rd，定義為如下應(yīng)變能密度的平均梯度與網(wǎng)格密度的比值

將網(wǎng)格強(qiáng)度在全域上的最大和最小值分別表示為

可定義介于Rmin和Rmax之間的兩網(wǎng)格強(qiáng)度RC1和RC2作為網(wǎng)格細(xì)化的標(biāo)準(zhǔn)，采用兩級(jí)的方式，即當(dāng)rd＞RC1時(shí)，僅進(jìn)行一層加密；當(dāng)rd＞RC2時(shí)，進(jìn)行兩層加密.

在本文采用的方法中，為計(jì)算網(wǎng)格強(qiáng)度rd，需要計(jì)算應(yīng)變能密度的導(dǎo)數(shù)和網(wǎng)格密度.考慮一個(gè)節(jié)點(diǎn)數(shù)為Nk，面積為A的背景積分網(wǎng)格k，節(jié)點(diǎn)的應(yīng)變能密度為

可通過(guò)積分網(wǎng)格k內(nèi)Nk個(gè)節(jié)點(diǎn)間能量密度的平均變化近似計(jì)算應(yīng)變能密度的導(dǎo)數(shù)

2.2 細(xì)化方案

由以上闡述的應(yīng)變能密度梯度準(zhǔn)則可確定需要加密的背景積分單元.如圖2所示，由于無(wú)網(wǎng)格法允許懸空節(jié)點(diǎn)的存在(如圖2(b)中的點(diǎn)D)，因此可方便地對(duì)積分網(wǎng)格直接進(jìn)行局部細(xì)化，并將其邊上中點(diǎn)作為新引入的計(jì)算節(jié)點(diǎn)，以實(shí)現(xiàn)節(jié)點(diǎn)的加密.對(duì)于需二層加密的積分網(wǎng)格，如圖3中的網(wǎng)格單元abc，可使用積分網(wǎng)格邊上的四等分點(diǎn)作為新節(jié)點(diǎn)對(duì)該網(wǎng)格進(jìn)行二次細(xì)化，同時(shí)還須對(duì)共用該網(wǎng)格節(jié)點(diǎn)的臨近網(wǎng)格進(jìn)行加密以保證網(wǎng)格分布的光滑過(guò)渡.顯然，無(wú)網(wǎng)格法形函數(shù)僅依賴于節(jié)點(diǎn)，不依賴于網(wǎng)格，因而允許出現(xiàn)懸空節(jié)點(diǎn)，這為計(jì)算節(jié)點(diǎn)的局部加密提供了極大便利.

圖2 一層網(wǎng)格加密示意圖Fig.2 Schematic diagram of one-level mesh refinemen

圖3 兩層網(wǎng)格加密示意圖Fig.3 Schematic diagram of two-level mesh refinemen

應(yīng)注意的是，積分網(wǎng)格的局部加密還影響了式(7)中網(wǎng)格邊界積分的計(jì)算.如圖4所示，按式(7)計(jì)算單元acd內(nèi)3個(gè)積分點(diǎn)上的形函數(shù)導(dǎo)數(shù)時(shí)，邊ac的積分須使用圖示的4個(gè)積分點(diǎn)，而非圖1所示的2個(gè)積分點(diǎn)，這是CEFG方法應(yīng)用于自適應(yīng)分析須特別注意之處.

圖4 自適應(yīng)QC3積分方案示意圖Fig.4 Schematic diagram of the adaptive QC3 integration scheme

其中δ為介于0和1之間的一個(gè)計(jì)算參數(shù).通過(guò)選取合適的δ，當(dāng)式(24)得到滿足時(shí)，即相鄰兩次計(jì)算的應(yīng)變能相對(duì)誤差小于δ時(shí)，自適應(yīng)過(guò)程終止，此即是本文所采用的自適應(yīng)停止準(zhǔn)則，詳見(jiàn)文獻(xiàn)[27].

3 數(shù)值算例

本節(jié)將通過(guò)數(shù)值算例考察和驗(yàn)證本文發(fā)展的自適應(yīng)CEFG方法的有效性和優(yōu)越性.所有算例的參數(shù)均無(wú)量綱化，彈性模量E=1×107，泊松比ν=0.3.為作比較，對(duì)同樣采用應(yīng)變能密度梯度自適應(yīng)準(zhǔn)則的標(biāo)準(zhǔn)EFG方法也進(jìn)行了程序?qū)崿F(xiàn)，并稱之為自適應(yīng)EFG方法，它在每個(gè)三角形積分子域內(nèi)使用16個(gè)積分點(diǎn)，積分點(diǎn)上形函數(shù)的導(dǎo)數(shù)采用標(biāo)準(zhǔn)形式(即對(duì)形函數(shù)直接求導(dǎo))，而非由式(11)求得的修正導(dǎo)數(shù)，這是這兩種方法的本質(zhì)區(qū)別所在.

3.1 方板圓孔問(wèn)題

該算例是彈性力學(xué)的經(jīng)典考題之一，計(jì)算域尺寸、邊界條件和載荷如圖5所示.初始的計(jì)算節(jié)點(diǎn)分布和積分網(wǎng)格(節(jié)點(diǎn)典型間距h=0.5)如圖6所示.

自適應(yīng)方案中的兩個(gè)網(wǎng)格強(qiáng)度參數(shù)分別取為RC1=20%Rmax和RC2=80%Rmax，自適應(yīng)停止指標(biāo)δ=1×10-4.圖7顯示了本文方法經(jīng)過(guò)自適應(yīng)節(jié)點(diǎn)加密后的計(jì)算節(jié)點(diǎn)分布和積分網(wǎng)格.可以看出，本文方法能夠準(zhǔn)確地對(duì)存在應(yīng)力集中的小孔周邊區(qū)域進(jìn)行自適應(yīng)加密.而且，本文方法能夠十分精確地得出圖5中A點(diǎn)處的應(yīng)力集中系數(shù)為3.注意到，該問(wèn)題的解析解表明A點(diǎn)附近的應(yīng)力集中程度遠(yuǎn)高于B點(diǎn)附近，相應(yīng)地，本文自適應(yīng)方法得到的計(jì)算節(jié)點(diǎn)分布在A點(diǎn)的密集程度也遠(yuǎn)高于B點(diǎn).這些結(jié)果表明了本文方法的正確性及其處理局部應(yīng)力集中問(wèn)題的有效性.

圖5 方板圓孔問(wèn)題示意圖Fig.5 Schematic diagram of the plate with a hole problem

圖6 方板圓孔問(wèn)題的初始配置Fig.6 Initial set-up for the plate with a hole problem

圖7 方板圓孔問(wèn)題的自適應(yīng)結(jié)果Fig.7 Adaptivity for the plate with a hole problem

與已有的CEFG方法相比，本文方法的優(yōu)勢(shì)在于所引入的自適應(yīng)節(jié)點(diǎn)加密方案，這一點(diǎn)可由圖8看出.顯然，要取得同樣的高精度，本文發(fā)展的自適應(yīng)CEFG方法所需的計(jì)算節(jié)點(diǎn)數(shù)比CEFG方法要少得多，因而可大幅度減少CPU時(shí)間.

圖8 方板圓孔問(wèn)題的位移誤差--節(jié)點(diǎn)數(shù)曲線Fig.8 Displacement error-number of nodes curve of the plate with a hole problem

與同樣采用自適應(yīng)技術(shù)的EFG方法相比，本文方法在計(jì)算精度和計(jì)算效率方面具有顯著優(yōu)勢(shì)，這集中反映在表1中.而且，如圖9所示，本文發(fā)展的自適應(yīng)CEFG方法得到的應(yīng)力場(chǎng)要比自適應(yīng)EFG方法的好得多.這些數(shù)值結(jié)果表明本文方法的優(yōu)越性.

圖9 方板圓孔問(wèn)題σyy應(yīng)力場(chǎng)比較Fig.9 Comparison of σyystress fiel of the plate with a hole problem

表1 方板圓孔問(wèn)題的兩種自適應(yīng)方法比較Table 1 Comparison of the two adaptive methods for the plate with a hole problem

3.2 受壓半無(wú)限平面問(wèn)題

該算例也是彈性力學(xué)的經(jīng)典考題，其解析解見(jiàn)文獻(xiàn)[32].如圖10所示，計(jì)算域取為3×3區(qū)域，上邊界的固定位移和右邊界的載荷均按解析解施加.

初始的計(jì)算節(jié)點(diǎn)配置取為節(jié)點(diǎn)典型間距h=0.25的均勻分布，自適應(yīng)停止指標(biāo)δ=5×10-5.圖11顯示了本文自適應(yīng)方法得到的計(jì)算節(jié)點(diǎn)分布和相應(yīng)的積分網(wǎng)格.顯然，它能自適應(yīng)地對(duì)下邊界受壓區(qū)域進(jìn)行局部的節(jié)點(diǎn)加密，得到合理的計(jì)算節(jié)點(diǎn)分布，因而能比未采用自適應(yīng)技術(shù)的CEFG方法節(jié)省大量的計(jì)算節(jié)點(diǎn)，如圖12所示.

圖10 受壓半無(wú)限平面問(wèn)題示意圖Fig.10 Schematic diagram of the pressure-loaded half plane problem

圖11 受壓半無(wú)限平面問(wèn)題的自適應(yīng)結(jié)果Fig.11 Adaptivity for the pressure-loaded half plane problem

圖12 受壓半無(wú)限平面問(wèn)題的位移誤差--節(jié)點(diǎn)數(shù)曲線Fig.12 Displacement error-number of nodes curve of the pressure-loaded half plane problem

表2比較了本文發(fā)展的自適應(yīng)CEFG方法與已有的自適應(yīng)EFG方法計(jì)算該考題的誤差和消耗的CPU時(shí)間，顯然，本文方法具有更好的精度和效率.而且，如圖13所示，本文方法能得到更好的應(yīng)力場(chǎng).

表2 受壓半無(wú)限平面問(wèn)題的兩種自適應(yīng)方法比較Table 2 Comparison of the two adaptive methods for the pressure-loaded half plane problem

3.3 變體力板

該算例考察本文方法對(duì)于非均勻體力造成局部高梯度問(wèn)題的有效性.如圖14所示，考慮一4×4的平板受到如下的非均勻體力

由以上的解析解可以看出，隨著到中心原點(diǎn)處距離的增加，位移和應(yīng)力均急劇地下降為0，因而在中心附近造成了局部高梯度的存在.

圖13 受壓半無(wú)限平面問(wèn)題σxx應(yīng)力場(chǎng)比較Fig.13 Comparison of the σxxstress fiel of the pressure-loaded half plane problem

圖14 變體力板示意圖及初始節(jié)點(diǎn)分布Fig.14 Schematic diagram of the plate with non-constant body force problem and its initial node distribution

如圖 14所示，初始節(jié)點(diǎn)分布均勻，典型間距h=0.2，自適應(yīng)停止指標(biāo)所有邊界均按解析解施加位移固定邊界條件.

圖15顯示了本文方法得到的節(jié)點(diǎn)分布和相應(yīng)的積分網(wǎng)格，它能夠自適應(yīng)地在中心附近進(jìn)行局部節(jié)點(diǎn)加密.兩種自適應(yīng)方法的計(jì)算結(jié)果比較如表3所示.顯然，本文方法具有更好的計(jì)算精度和效率.圖16顯示了兩種自適應(yīng)方法得到的應(yīng)力場(chǎng)，也可以看出，本文自適應(yīng)CEFG方法更加準(zhǔn)確.

圖15 變體力板問(wèn)題的自適應(yīng)結(jié)果Fig.15 Adaptivity for the plate with non-constant body force problem

表3 變體力板問(wèn)題兩種自適應(yīng)方法比較Table 3 Comparison of the two adaptive methods for the plate with non-constant body force problem

圖16 變體力板問(wèn)題σxx應(yīng)力場(chǎng)比較Fig.16 Comparison of the σxxstress fiel of the plate with non-constant body force problem

4 結(jié)論

本文在已有的一致性高階無(wú)網(wǎng)格法的基礎(chǔ)之上，充分利用無(wú)網(wǎng)格法易于實(shí)現(xiàn)自適應(yīng)計(jì)算的優(yōu)勢(shì)，采用應(yīng)變能密度梯度準(zhǔn)則驅(qū)動(dòng)自適應(yīng)過(guò)程，發(fā)展了自適應(yīng)一致性高階無(wú)單元伽遼金法.數(shù)值結(jié)果表明，該方法在分析應(yīng)力集中等局部高梯度問(wèn)題時(shí)，大量減少了所需的節(jié)點(diǎn)數(shù)目，顯著提高了一致性高階無(wú)網(wǎng)格法分析該類問(wèn)題的計(jì)算效率.與已有的標(biāo)準(zhǔn)無(wú)網(wǎng)格法的自適應(yīng)分析相比，所發(fā)展方法在計(jì)算效率、精度、應(yīng)力場(chǎng)光滑度等方面也均展現(xiàn)出了顯著優(yōu)勢(shì).

應(yīng)說(shuō)明的是，對(duì)于由幾何突變?cè)斐傻膽?yīng)力集中問(wèn)題(如算例1)，容易判別出應(yīng)力集中區(qū)域因而可以事先進(jìn)行手動(dòng)加密.然而，對(duì)于從幾何外形上無(wú)法判別應(yīng)力集中區(qū)域的問(wèn)題(如算例3)，自適應(yīng)分析則十分有用.再考慮到局部高梯度區(qū)域隨計(jì)算發(fā)生演化的問(wèn)題，如裂尖高梯度應(yīng)力場(chǎng)隨裂紋擴(kuò)展而移動(dòng)等，所發(fā)展的自適應(yīng)一致性高階無(wú)單元伽遼金法在分析此類問(wèn)題時(shí)應(yīng)具有顯著的優(yōu)勢(shì).當(dāng)然，它在這些復(fù)雜問(wèn)題中的具體應(yīng)用仍有賴于深入研究.

1 黃凱,白鴻柏,路純紅等.無(wú)網(wǎng)格法及其在金屬成形中的應(yīng)用綜述.鍛壓技術(shù),2016,41(2)：12-16(Huang Kai,Bai Hongbai,Lu Chunhong,et al.Summarization on meshfree method and its application in metal forming.Forming and Stamping Technology,2016, 41(2)：12-16(in Chinese))

2 Hah ZH,Youn SK.Eulerian analysis of bulk metal forming processes based on spline-based meshfree method.Finite Elements in Analysis and Design,2015,106：1-15

3 Muthu N,Falzon BG,Maiti SK,et al.Modifie crack closure integral technique for extraction of SIFs in meshfree methods.Finite Elements in Analysis and Design,2014,78：25-39

4 馬文濤,許艷,馬海龍.修正的內(nèi)部基擴(kuò)充無(wú)網(wǎng)格法求解多裂紋應(yīng)力強(qiáng)度因子.工程力學(xué),2015,32(10)：18-24(Ma Wentao,Xu Yan,Ma Hailong.Solving stress intensity factors of multiple cracks by using a modifie intrinsic basis enriched meshless mehthod.Engineering Mechanics,2015,32(10)：18-24(in Chinese))

5 Tanaka S,Suzuki H,Sadamoto S,et al.Analysis of cracked shear deformable plates by an ef f ective meshfree plate formulation.Engineering Fracture Mechanics,2015,144：142-157

6 Wang DD,WuJC.Anefficientnestingsub-domain gradient smoothing integration algorithm with quadratic exactness for Galerkin meshfree methods.Comput Methods Appl Mech Engrg,2016,298：485-519

7 Hillman M,Chen JS.An accelerated,convergent,and stable nodal integration in Galerkin meshfree methods for linear and nonlinear mechanics.Int J Numer Meth Engng,2016,107：603-630

8 Duan QL,Li XK,Zhang HW,et al.Second-order accurate derivatives and integration schemes for meshfree methods.International Journal for Numerical Methods in Engineering,2012,92(4)：399-424

9 Duan QL,Gao X,Wang BB,et al.Consistent element-free Galerkin method.International Journal for Numerical Methods in Engineering,2014,99：79-101

10 王冰冰,高欣,段慶林.熱傳導(dǎo)的二階一致1點(diǎn)積分無(wú)網(wǎng)格法.工程力學(xué),2014,31(9)：1-6(Wang Bingbing,Gao Xin,Duan Qinglin. Meshfree method for heat conduction using quadratic consistent 1-point integration scheme.Engineering Mechanics,2014,31(9)：1-6 (in Chinese))

11 王冰冰,陳嵩濤,段慶林.動(dòng)力分析的二階一致無(wú)網(wǎng)格法.應(yīng)用力學(xué)學(xué)報(bào),2014,31(3)：353-358(Wang Bingbing,Chen Songtao, Duan Qinglin.Quadratically consistent meshfree method for dynamics.Chinese Journal of Applied Mechanics,2014,31(3)：353-358(in Chinese))

12 Duan Qinglin,Gao Xin,Wang Bingbing,et al.A four-point integrationschemewithquadraticexactnessforthree-dimensionalelementfree Galerkin method based on variationally consistent formulation.Comput Methods Appl Mech Engrg,2014,280(10)：84-116

13 Ortiz-Bernardin A,Hale JS,Cyron CJ.Volume-averaged nodal projection method for nearly-incompressible elasticity using meshfree and bubble basis functions.Comput Methods Appl Mech Engrg, 2015,285：427-451

14 Ortiz-Bernardin A,Puso MA,Sukumar N.Improved robustness for nearly-incompressible large deformation meshfree simulations on Delaunay tessellations.Comput Methods Appl Mech Engrg,2015, 293：348-374

15 陳炎,曹樹(shù)良,梁開(kāi)洪等.結(jié)合自適應(yīng)網(wǎng)格的基于特征線方程的分離算法.力學(xué)學(xué)報(bào),2009,41(5)：666-670(Chen Yan,Cao Shuliang,Liang Kaihong,et al.A new grid adaptive strategy combined characterized based split method.Chinese Journal of Theoretical and Applied Mechanics,2009,41(5)：666-670(in Chinese))

16 徐云,陳軍,蔚喜軍.多尺度材料模擬的自適應(yīng)有限元方法研究.力學(xué)學(xué)報(bào),2009,41(5)：722-729(Xu Yun,Chen Jun,Yu Xijun.A new adaptive finit element method for multiscale dynamic simulation.Chinese Journal of Theoretical and Applied Mechanics,2009, 41(5)：722-729(in Chinese))

17 Nguyen-Xuan H,Wu CT,Liu GR.An adaptive selective ES-FEM for plastic collapse analysis.European Journal of Mechanics A/Solids,2016,58：278-290

18 Duarte CA,Oden JT.An h-p adaptive method using clouds.Comput Methods Appl Mech Engrg,1996,139：237-262

19 劉欣,朱德懋,陸明萬(wàn)等.平面裂紋問(wèn)題的 h,p,hp型自適應(yīng)無(wú)網(wǎng)格方法的研究.力學(xué)學(xué)報(bào),2000,32(3)：308-318(Liu Xin,Zhu Demao,Lu Mingwan,et al.H,p,hp adaptive meshless method for plane crack problem.Acta Mechanica Sinica,2000,32(3)：308-318 (in Chinese))

20 Chung HJ,Belytschko T.An error estimate in the EFG method.Comp Mech,1998,21：91-100

21 Sethuraman R,Hormis JV.An error estimator and adaptivity based on basis coefficient-vector fielof element-free Galerkin method.International Journal of Computational Methods,2011,8(1)：91-118

22 He YQ,Yang HT,Deeks AJ.A node based error estimator for the element-free Galerkin(EFG)methods.International Journal of Computational Methods,2014,11(4)：13500591-24

23 Hiiussler-Combe U,Korn C.An adaptive approach with the element-free Galerkin method.Comput Methods Appl Mech Engrg, 1998,162：203-222

24 Rabczuk T,Belytschko T.Adaptivity for structured meshfree particle methods in 2D and 3D.Int J Numer Meth Engng,2005,63：1559-1582

25 Liu GR,Tu ZH.An adaptive procedure based on background cells for meshless methods.Comput Methods Appl Mech Engrg,2002, 191：1923-1943

26 干年妃,李光耀,鐘志華等.金屬體積成型過(guò)程的自適應(yīng)無(wú)網(wǎng)格方法.中國(guó)機(jī)械工程,2006,17(24)：2612-2617(Gan Nianfei,Li Guangyao,Zhong Zhihua,et al.Adaptive meshfree method on bulk forming process.China Mechanical Engineering,2006,17(24)：2612-2617(in Chinese))

27 Luo YH,Hiiussler-Combe U.A gradient-based adaptation procedure and its implementation in the element-free Galerkin method.Int J Numer Meth Engng,2003,56：1335-1354

28 張征,劉更,劉天祥.接觸問(wèn)題的自適應(yīng)無(wú)網(wǎng)格伽遼金方法.中國(guó)機(jī)械工程,2007,18(23)：2868-2873(Zhang Zheng,Liu Geng,Liu Tianxiang.An adaptive element-free method Galerkin model for contact problems.China Mechanical Engineering,2007,18(23)： 2868-2873(in Chinese))

29 Metsis P,Lantzounis N,Papadrakakis M.A new hierarchical partition of unity formulation of EFG meshless methods.Comput Methods Appl Mech Engrg,2015,283：782-805

30 Belytschko T,Lu Y,Gu L.Element-free Galerkin methods.International Journal for Numerical Methods in Engineering,1994,37(2)：229-256

31 Fern′andez-M′endez S,Huerta A.Imposing essential boundary conditions in mesh-free methods.Computer Methods in Applied Mechanics and Engineering,2004,193(12)：1257-1275

32 Duan QL,Belytschko T.Gradient and dilatational stabilizations for stress-point integration in the element-free Galerkin method.International Journal for Numerical Methods in Engineering,2009, 77(6)：776-798

ADAPTIVE CONSISTENT HIGH ORDER ELEMENT-FREE GALERKIN METHOD1)

Shao Yulong?Duan Qinglin?,?,??,2)Gao Xin?Li Xikui?Zhang Hongwu?

?(State Key Laboratory of Structural Analysis for Industrial Equipment，Dalian University of Technology，Dalian116024，China)

?(State Key Laboratory of Geohazard Prevention and Geoenvironment Protection，Chengdu University of Technology，Chengdu610059，China)

??(State Key Laboratory of Water Resources and Hydropower Engineering Science，Wuhan University，Wuhan430072，China)

The recently developed consistent high order element-free Galerkin(EFG)method not only dramatically reduces the number of quadrature points in domain integration but also accurately passes the linear and quadratic patch tests,and remarkably improves the computational efficiency,accuracy and convergence of the standard EFG methods.On this basis,this work presents the h-adaptive analysis for consistent high order EFG method by taking advantage of the convenience of the EFG method in adding approximation nodes locally.The proposed method adaptively determines the region which needs nodal refinemen according to the gradient of the strain energy density.The generation of the new approximation nodes is based on the multi-level local mesh refinemen of the background integration mesh.The gradualtransition between the regions with and without nodal refinemen is also considered.The relative error of the strain energy in two successive computation is adopted as the stop-criterion of the adaptive process.The proposed adaptive meshfree method is applied to the analysis of stress concentration caused by geometry,external boundary loads and body forces.Numerical results show that the developed method is able to refin the region with high stress gradient adaptively and to generate reasonable distribution of approximation nodes automatically.In comparison with the existing adaptive schemes of the standard EFG method,the proposed method shows remarkable advantages on computational efficiency, accuracy and the smoothness of the resulting stress fields In comparison with the consistent high order EFG method using uniform nodal distribution,the proposed adaptive method dramatically reduces the number of computational nodes. As a consequence,it significantl improves the computational efficiency and accuracy of the consistent high order EFG method for the analysis of problems with local high gradients such as stress concentration.

element-free Galerkin method，meshfree/meshless methods,adaptive analysis，strain energy density,stress concentration

O343.1

A doi：10.6052/0459-1879-16-252

2016-09-07收稿，2016-11-04錄用,2016-11-07網(wǎng)絡(luò)版發(fā)表.

1)國(guó)家自然科學(xué)基金 (11232003,11372066)、中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金 (DUT15LK07)、遼寧省教育廳重點(diǎn)實(shí)驗(yàn)室基礎(chǔ)研究(LZ2014002)、水資源與水電工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室開(kāi)放基金(2015SGG03)和地質(zhì)災(zāi)害防治與地質(zhì)環(huán)境保護(hù)國(guó)家重點(diǎn)實(shí)驗(yàn)室開(kāi)放基金(SKLGP2016K007)資助項(xiàng)目.

2)段慶林，副教授，博士，主要研究方向?yàn)闊o(wú)網(wǎng)格法、材料破壞分析與模擬等.E-mail：qinglinduan@dlut.edu.cn

邵玉龍,段慶林,高欣,李錫夔,張洪武.自適應(yīng)一致性高階無(wú)單元伽遼金法.力學(xué)學(xué)報(bào),2017,49(1)：105-116

Shao Yulong,Duan Qinglin,Gao Xin,Li Xikui,Zhang Hongwu.Adaptive consistent high order element-free Galerkin method.Chinese Journal of Theoretical and Applied Mechanics,2017,49(1)：105-116