徐亞飛, 肖映雄, 吳宇航

(湘潭大學(xué) 土木工程與力學(xué)學(xué)院,湘潭 411105)

1 引言

在有限元計算中,影響其結(jié)果準(zhǔn)確程度的因素很多,如單元類型的選取和網(wǎng)格的疏密等。對于給定的實際問題,當(dāng)選取了合適單元類型后,則需根據(jù)問題特點(diǎn)和計算精度設(shè)置網(wǎng)格的疏密。如對應(yīng)力集中問題,需要采用自適應(yīng)網(wǎng)格才能獲得滿足給定精度的數(shù)值解。自適應(yīng)網(wǎng)格技術(shù),包括自適應(yīng)加密或自適應(yīng)網(wǎng)格重劃,相應(yīng)的自適應(yīng)有限元方法可減輕計算精度與計算量的矛盾[1-5]。目前,在自適應(yīng)有限元方面也取得了大量的研究成果,如文獻(xiàn)[6,7]為膜結(jié)構(gòu)極小曲面找形和自由振動問題提出了相應(yīng)的自適應(yīng)有限元方法;文獻(xiàn)[8,9]設(shè)計了基于局部誤差控制標(biāo)準(zhǔn)的p型自適應(yīng)有限元法,并應(yīng)用于混凝土骨料模型的數(shù)值模擬;文獻(xiàn)[10]提出了一種基于自適應(yīng)網(wǎng)格加密和人工神經(jīng)網(wǎng)絡(luò)的有限元降階模型;文獻(xiàn)[11]提出了變截面變曲率梁振型的有限元后處理超收斂拼片恢復(fù)方法,獲得了相應(yīng)的網(wǎng)格自適應(yīng)分析方法。各類新型自適應(yīng)方法為解決實際工程問題提供了高效計算方案。

一般而言,可在Abaqus軟件平臺上通過GUI操作進(jìn)行有限元計算與模擬,但當(dāng)需要多次和反復(fù)進(jìn)行時,或者需要重復(fù)改動某些參數(shù)和改動原來在GUI中不能改動的參數(shù)時,迫切需要利用Python對Abaqus進(jìn)行二次開發(fā)。利用Python語言可以利用Abaqus腳本接口繞過Abaqus/CAE的圖形用戶界面GUI,直接對Abaqus內(nèi)核進(jìn)行操作。目前,在Abaqus二次開發(fā)及應(yīng)用方面,已有很多研究成果[12-17]。二次開發(fā)技術(shù)還運(yùn)用到了其他相關(guān)的領(lǐng)域中,如研究裂紋的擴(kuò)展[18],在金屬切削仿真中防止網(wǎng)格的畸變[19]?；赑ython編寫的腳本可以Abaqus直接運(yùn)行,可方便修改有限元模型和相關(guān)參數(shù),提高建模效率和避免返工。

本文采用以前學(xué)者發(fā)展的網(wǎng)格重劃分算法,利用Python-Abaqus進(jìn)行有限元建模和模擬,從數(shù)值上分析了誤差控制標(biāo)準(zhǔn)對計算結(jié)果的影響,實現(xiàn)了自適應(yīng)求解全過程。所編Python腳本可繞過Abaqus/CAE的圖形用戶界面GUI,直接對Abaqus內(nèi)核進(jìn)行操作,從而實現(xiàn)從幾何建模、網(wǎng)格剖分到自適應(yīng)有限元求解的自動化處理。然后,將所編Python腳本應(yīng)用于含有凹角或孔洞的非光滑物理區(qū)域、材料常數(shù)不連續(xù)和應(yīng)力集中等幾類典型問題的有限元分析中。結(jié)果表明,基于Abaqus網(wǎng)格重劃技術(shù)的自適應(yīng)方法對求解帶奇性解的實際問題是非常有效的,通過迭代可獲取最優(yōu)網(wǎng)格,近似解可達(dá)到要求的求解精度。Python二次開發(fā)自適應(yīng)計算與模擬,可方便多次修改模型和參數(shù),實現(xiàn)全流程腳本控制,提高建模效率。

2 自適應(yīng)網(wǎng)格重劃算法及其實現(xiàn)

在實際應(yīng)用中,由于無法事先決定網(wǎng)格的疏密分布,只能憑經(jīng)驗在梯度可能較大的部位進(jìn)行網(wǎng)格加密/重劃,使得網(wǎng)格加密/重劃帶有一定的盲目性。為了克服這種困難,需要采用自適應(yīng)網(wǎng)格技術(shù),其可將誤差估計與有限元網(wǎng)格的生成和改正有機(jī)結(jié)合起來,力求用較少的單元獲得預(yù)期的精度要求。對于一個特定的網(wǎng)格,在精確解未知的情況下,需要采用后驗誤差估計的方法來定量評價離散誤差的大小,以便指示下一步的網(wǎng)格劃分。在求得給定網(wǎng)格下有限元位移解和應(yīng)力解后,采用應(yīng)力恢復(fù)等方法構(gòu)造出改正的應(yīng)力解。研究表明,改正前的有限元應(yīng)力解和改正后的應(yīng)力解之差可用來估計有限元的誤差。這種基于后驗誤差估計的自適應(yīng)方法已廣泛應(yīng)用于各種實際問題的求解中,具有很好的普適性。

本節(jié)考慮基于后驗誤差估計的自適應(yīng)網(wǎng)格重劃算法,利用Python-Abaqus進(jìn)行建模和模擬分析,實現(xiàn)自適應(yīng)有限元求解全過程。

2.1 一般的自適應(yīng)網(wǎng)格重劃算法

自適應(yīng)網(wǎng)格重劃可以改善有限元計算結(jié)果的精度。一般如遇以下幾種情形則需要采用自適應(yīng)網(wǎng)格重劃技術(shù),(1) 事先不確定有限元網(wǎng)格需要怎樣的細(xì)化來達(dá)到給定的計算精度要求; (2) 很難在感興趣的區(qū)域(如應(yīng)力集中)附近事先設(shè)計一個足夠細(xì)化的網(wǎng)格; (3) 事先不知道網(wǎng)格需要細(xì)化的位置,如塑性區(qū)域的形成。下面給出自適應(yīng)網(wǎng)格重劃算法[1]的一般描述。

設(shè)Th是求解域Ω的任意四邊形網(wǎng)格剖分,記單元數(shù)為Ne,其中,h為Th上所有剖分單元的最大尺寸。對給定的線性或非線性力學(xué)問題,假設(shè)在當(dāng)前網(wǎng)格Th下,利用有限元法求解得到近似位移場{u}h、應(yīng)變場{ε}h和應(yīng)力場{σ}h。設(shè)原問題的精確位移場為{u},相應(yīng)的應(yīng)變場和應(yīng)力場分別為{ε}和{σ}。定義誤差{eu}={u}-{u}h,{eε}={ε}-{ε}h以及{eσ}={σ}-{σ}h。為了獲取無量綱的相對誤差度量,分別引入誤差和精確解的能量范數(shù),

(1)

(2)

(3)

(i=1,2,…,Ne)(4)

定義如下形式的單元誤差指標(biāo),

(i=1,2,…,Ne)(5)

若ξi≤1對所有單元均成立,則當(dāng)前網(wǎng)格已是最優(yōu),自適應(yīng)計算結(jié)束;否則,按式(6)計算新的單元尺寸,

(對ξi>1)(6)

在實際應(yīng)用中,應(yīng)力精確解{σ}一般未知,可通過磨平技術(shù)得到改善后的近似解{σ*}代替{σ},進(jìn)而計算誤差范數(shù)。為保證{σ*}與相應(yīng)的離散應(yīng)力場{σ}h的預(yù)測值一致,需要滿足

(7)

(8)

由此可得

(9)

評注對于四邊形等參元,由于Gauss積分點(diǎn)為超收斂點(diǎn),其上的應(yīng)力近似解比其他部位要具有更高的精度。在實際應(yīng)用中,可以利用精度較高的Gauss積分點(diǎn)的應(yīng)力值來改正節(jié)點(diǎn)應(yīng)力值,再利用磨平技術(shù)(如繞節(jié)點(diǎn)加權(quán)平均),得到改善后的節(jié)點(diǎn)應(yīng)力值{σ*}。

上述介紹的是計算應(yīng)力近似解{σ*}的一種簡單方法,精度更高的方法可采用文獻(xiàn)[20]提出的應(yīng)力恢復(fù)投影法(patch recovery technique of Z-Z method)或SPR方法[2]。

2.2 自適應(yīng)有限元算法Python-Abaqus實現(xiàn)

Abaqus提供了自適應(yīng)網(wǎng)格重劃算法的GUI操作,但當(dāng)需要多次、反復(fù)進(jìn)行計算與模擬,且需要一些更高級的技巧來實現(xiàn)Abaqus中原來沒有的一些功能或者需要重復(fù)改動某些參數(shù),可以利用Python對Abaqus進(jìn)行二次開發(fā),編制相應(yīng)的Python腳本,進(jìn)而可直接在腳本中進(jìn)行添加新功能或修改參數(shù)等,非常方便地實現(xiàn)自適應(yīng)網(wǎng)格重劃、完成有限元計算或模擬。

在Abaqus中,要實現(xiàn)自適應(yīng)網(wǎng)格重劃算法,需要選取下述三個重要指標(biāo)。

(1) 誤差指標(biāo)變量。在分析中選取什么誤差指標(biāo)變量來實施自適應(yīng)網(wǎng)格重劃,需要考慮以下三個因素來進(jìn)行選擇, ① 誤差指標(biāo)的特征,② 分析模型場函數(shù)的類型, ③ 載荷的性質(zhì)。

在Abaqus中,常用誤差指標(biāo)變量有單元能量密度(element energy density;ENDENERI)、密塞斯應(yīng)力(Mises stress;MISESRI)、等效塑性應(yīng)變(equivalent plastic strain;PEEQERI)、塑性應(yīng)變(plastic strain;PEERI)、蠕變(creep strain;CEERI)、熱流(Heat flux;HFLERI)和電勢梯度等。本文主要選取單元能量密度和密塞斯應(yīng)力作為誤差指標(biāo)變量。

結(jié)合2.1節(jié)的一般自適應(yīng)網(wǎng)格重劃算法和Abaqus中給定的上述三個重要指標(biāo),可以得到自適應(yīng)有限元算法流程,如圖1所示。

圖1 自適應(yīng)有限元算法流程Fig.1 Flow chart of adaptive finite element algorithm

利用Python語言,自行編制了如圖1所示自適應(yīng)有限元算法的腳本文件。Python腳本實現(xiàn)網(wǎng)格重劃與自適應(yīng)計算的優(yōu)點(diǎn), (1) 可自行修改和重復(fù)性操作,且所需代碼較少,如可靈活修改誤差變量、重劃網(wǎng)格的方法及誤差目標(biāo)等關(guān)鍵參數(shù),非常方便; (2) 可實現(xiàn)各種循環(huán)操作以及數(shù)據(jù)存儲等的自動化處理; (3) 可繞過Abaqus/CAE的GUI,直接對Abaqus內(nèi)核進(jìn)行操作,如可直接調(diào)用Abaqus的核心模塊adaptivityProcesses()實現(xiàn)自適應(yīng)網(wǎng)格重劃及后續(xù)計算,避開了在不同對話框和標(biāo)簽頁下輸入數(shù)據(jù)等繁瑣操作。

3 數(shù)值試驗及結(jié)果分析

將所編Python腳本應(yīng)用于含有凹角或孔洞的非光滑物理區(qū)域、材料常數(shù)不連續(xù)和應(yīng)力集中等幾類典型問題的有限元分析中,以考察全流程Python腳本的適應(yīng)性和網(wǎng)格重劃算法的有效性。

算例1(含圓孔的無限大薄板問題) 如圖2所示,設(shè)長L=20 cm,寬W=10 cm,圓孔的圓心位于薄板的中心,半徑R=1 cm,彈性模量E=2.1 GPa,泊松比ν=0.3,取單位厚度,均布荷載為q=10 kN/cm2。利用對稱性,取模型的1/4作為研究對象,對應(yīng)的有限元計算模型如圖3所示。

圖2 帶孔薄板Fig.2 Rectangular plate with a small hole

圖3 有限元計算模型Fig.3 Finite element model

利用所編Python腳本進(jìn)行分析,重點(diǎn)分析不同誤差控制標(biāo)準(zhǔn)對計算結(jié)果的影響。自適應(yīng)加密區(qū)域按不限定區(qū)域即自由劃分,數(shù)值結(jié)果總結(jié)如下。

圖4 圓孔上的應(yīng)力數(shù)值解與精確解比較Fig.4 Comparisons for the numerical solutions and the exact solution of the stress on the circular hole

圖5 當(dāng)=1%時對應(yīng)的初始網(wǎng)格和最終網(wǎng)格Fig.5 Initial mesh and the final mesh when =1%

圖6 當(dāng)=1%時對應(yīng)的初始網(wǎng)格和最終網(wǎng)格Fig.6 Initial mesh and the final mesh when =1%

圖7 兩種不同誤差控制標(biāo)準(zhǔn)下計算結(jié)果比較Fig.7 Comparisons of the computational results under two different error control standards

本文繼續(xù)將所編Python腳本應(yīng)用于含有凹角或孔洞的非光滑物理區(qū)域、材料常數(shù)不連續(xù)和應(yīng)力集中等問題的有限元分析中,以進(jìn)一步考察算法的有效性和Python腳本的適應(yīng)性。

圖8 開孔非規(guī)則薄板Fig.8 Irregular thin plate with holes

圖9 當(dāng)=1%時,自適應(yīng)迭代初始網(wǎng)格和最終網(wǎng)格Fig.9 Initial mesh and the final mesh when =1%

圖10 當(dāng)=1%時,開孔非規(guī)則薄板應(yīng)力云圖Fig.10 Stress nephograms of irregular thin plate with openings when is equal to 1%

(雙連拱隧道問題[22]) 考慮某雙連拱隧道,其基本結(jié)構(gòu)如圖11所示。假設(shè)圍巖及結(jié)構(gòu)支護(hù)體系材料為均勻各向同性的線彈性材料,對應(yīng)的材料參數(shù)列入表1。為較好地模擬隧道的受力情況,需要對外圍邊界條件進(jìn)行近似處理,即左右邊界取離隧道中心距離為60 m,上下邊界取距離隧道中心距離為40 m。按平面應(yīng)變問題分析隧道圍巖和襯徹等在自重作用下的變形和應(yīng)力分布。

圖11 雙連拱隧道結(jié)構(gòu)Fig.11 Double arch tunnel structure

表1 材料參數(shù)Tab.1 Material parameters

圖12 初始網(wǎng)格和自適應(yīng)迭代得到的最終網(wǎng)格重劃Fig.12 Initial mesh and the final mesh obtained by adaptive iteration

圖13 滿足精度要求的應(yīng)力云圖Fig.13 Stress nephograms of meeting the given accuracy

從上述算例的數(shù)值試驗結(jié)果可知,本文所編全流程Python腳本對求解含有凹角或孔洞的非光滑物理區(qū)域、材料常數(shù)不連續(xù)和應(yīng)力集中等問題都是非常有效的。按照不限區(qū)域方式實施自適應(yīng)網(wǎng)格重劃,對處理實際問題也是非常方便,均能在奇性解附近自動加密網(wǎng)格,進(jìn)一步驗證了算法的有效性和Python腳本的適應(yīng)性。另外,對多介質(zhì)問題,還可以對每一種介質(zhì)單獨(dú)設(shè)置自適應(yīng)網(wǎng)格重劃方式及選取不同的誤差控制標(biāo)準(zhǔn),應(yīng)用靈活方便。在Python腳本中,可直接添加新模型以及修改相關(guān)參數(shù),從而方便地實現(xiàn)對不同問題的自適應(yīng)有限元計算或模擬。

4 結(jié) 論

自適應(yīng)有限元法是求解含有凹角或孔洞等的非光滑物理區(qū)域、材料常數(shù)不連續(xù)和應(yīng)力集中等實際工程問題的有效數(shù)值方法之一。在利用Abaqus進(jìn)行自適應(yīng)網(wǎng)格重劃與模擬時,若創(chuàng)建的模型復(fù)雜或存在大量重復(fù)操作時,二次開發(fā)則更具優(yōu)勢。本文目的是將已有的網(wǎng)格重劃分算法與Python腳本有機(jī)結(jié)合進(jìn)行二次開發(fā),編寫相應(yīng)的自適應(yīng)有限元程序,通過直接對Abaqus內(nèi)核進(jìn)行操作,可實現(xiàn)從幾何建模、網(wǎng)格剖分、定義材料屬性、加載和有限元求解到查看計算結(jié)果的自動化處理,即實現(xiàn)全流程腳本控制。本文重點(diǎn)分析了誤差控制標(biāo)準(zhǔn)的選取對計算結(jié)果的影響,并通過幾類典型問題驗證了Python腳本的適應(yīng)性及方法的有效性。

本文做了有限的數(shù)值試驗,當(dāng)功能強(qiáng)大的Abaqus與Python相結(jié)合,其功能將更加強(qiáng)大。本文算例主要考慮的是二維單介質(zhì)或多介質(zhì)線彈性問題,相應(yīng)的自適應(yīng)網(wǎng)格劃分已經(jīng)非常復(fù)雜。而在實際工程應(yīng)用中,更多的是非線性問題,將自適應(yīng)網(wǎng)格重劃方法擴(kuò)展到非線性問題意義更大,也更能發(fā)揮有限元方法的通用性優(yōu)勢。針對非線性問題的仿真和分析,一方面,需要進(jìn)一步測試Python腳本的適應(yīng)性和有效性;另一方面,還需要利用Abaqus提供的其他自適應(yīng)網(wǎng)格技術(shù),如ALE自適應(yīng)網(wǎng)格進(jìn)行Python編程,以提高諸如大變形問題的仿真分析效率。