基于矩形獨(dú)立覆蓋初步實(shí)現(xiàn)結(jié)構(gòu)靜力分析的自動(dòng)計(jì)算

蘇海東,龔亞琦,頡志強(qiáng),祁勇峰

（長(zhǎng)江科學(xué)院材料與結(jié)構(gòu)研究所,武漢 430010 )

?

基于矩形獨(dú)立覆蓋初步實(shí)現(xiàn)結(jié)構(gòu)靜力分析的自動(dòng)計(jì)算

蘇海東,龔亞琦,頡志強(qiáng),祁勇峰

（長(zhǎng)江科學(xué)院材料與結(jié)構(gòu)研究所,武漢 430010 )

摘 要：采用前期研究的基于矩形獨(dú)立覆蓋的新型數(shù)值流形法,提出結(jié)構(gòu)線彈性靜力分析的自動(dòng)計(jì)算方法,包括自動(dòng)的前處理、自適應(yīng)分析等。根據(jù)獨(dú)立覆蓋的特點(diǎn)提出幾個(gè)后驗(yàn)誤差指標(biāo):獨(dú)立覆蓋之間條形連接區(qū)域的應(yīng)變連續(xù)性指標(biāo)；邊界應(yīng)力指標(biāo)和獨(dú)立覆蓋的高階誤差指標(biāo)。利用新方法的h型網(wǎng)格加密及p型升階的方便性,選擇一種路徑嘗試h-p型的混合自適應(yīng),其中,對(duì)于矩形獨(dú)立覆蓋采用簡(jiǎn)單的二分法實(shí)現(xiàn)覆蓋加密。通過(guò)幾個(gè)二維算例驗(yàn)證了新方法實(shí)現(xiàn)自動(dòng)計(jì)算的可行性,只需人工輸入結(jié)構(gòu)外形、材料參數(shù)和邊界條件,其它工作完全交由計(jì)算機(jī)完成,最終得到滿足一定精度的計(jì)算結(jié)果。

關(guān)鍵詞：數(shù)值流形方法；獨(dú)立覆蓋；自動(dòng)計(jì)算；誤差估計(jì)；h-p混合自適應(yīng)；靜力分析

2016,33（02):144-150

1 研究背景

計(jì)算機(jī)輔助工程（Computer Aided Engineering, CAE),主要是指用計(jì)算機(jī)輔助求解復(fù)雜工程和產(chǎn)品的結(jié)構(gòu)強(qiáng)度、剛度、動(dòng)力響應(yīng)、彈塑性等力學(xué)性能的分析計(jì)算［1］。提高分析效率對(duì)CAE意義重大,未來(lái)的發(fā)展趨勢(shì)是讓計(jì)算機(jī)自動(dòng)、高效地完成分析,即CAE分析技術(shù)的自動(dòng)化。

本文提出自動(dòng)計(jì)算的思想,其目標(biāo)是基本做到無(wú)需人工參與就能獲得理想結(jié)果的分析,當(dāng)然,必要的結(jié)構(gòu)外形、材料參數(shù)和邊界條件的人工輸入除外。要達(dá)到此目標(biāo),必須解決CAE前處理的自動(dòng)化、計(jì)算收斂效率及自適應(yīng)等問(wèn)題。其中,自適應(yīng)是指在分析過(guò)程中由程序根據(jù)計(jì)算精度要求自動(dòng)地提高數(shù)值計(jì)算中的近似函數(shù)的逼近程度,包括h型（加密網(wǎng)格)和p型（提高局部近似函數(shù)階次)以及兩者的混合方式［2］。目前來(lái)看,即使對(duì)于最基本的結(jié)構(gòu)線彈性靜力分析而言,以有限元法為代表的各種數(shù)值方法離自動(dòng)計(jì)算的目標(biāo)仍有不小距離。

有限元法在CAE中應(yīng)用最為廣泛,但其前處理的網(wǎng)格劃分經(jīng)常會(huì)耗費(fèi)大量的人工工作量（據(jù)統(tǒng)計(jì),前處理在計(jì)算分析總時(shí)間中所占比例往往高達(dá)80%以上［1］)。其原因在于,有限元網(wǎng)格不僅要模擬結(jié)構(gòu)外形,還要保持合適的網(wǎng)格形狀,以及網(wǎng)格之間通過(guò)結(jié)點(diǎn)的正確連接,而實(shí)際工程結(jié)構(gòu)的復(fù)雜性使得網(wǎng)格剖分的難度較大。同時(shí),由于求解域內(nèi)物理場(chǎng)分布的復(fù)雜程度不同,需要網(wǎng)格密度的變化,比如目前廣泛使用的h型自適應(yīng)需要對(duì)精度不佳的區(qū)域進(jìn)行網(wǎng)格細(xì)分,涉及到的大小網(wǎng)格過(guò)渡通常不太方便,要保持網(wǎng)格合適形狀較為不易。在計(jì)算收斂效率方面,對(duì)于存在強(qiáng)奇異性的裂紋尖端（以下簡(jiǎn)稱裂尖)等復(fù)雜應(yīng)力場(chǎng)分布區(qū)域,用常規(guī)的有限單元插值函數(shù)逼近往往效果不佳。

廣義有限元在結(jié)點(diǎn)上引入級(jí)數(shù)提高近似函數(shù)逼近能力,但其主要問(wèn)題是高階近似函數(shù)情況下的方程組線性相關(guān)［3-4］。擴(kuò)展有限元在單元結(jié)點(diǎn)上增加附加函數(shù):加入階躍函數(shù)模擬裂紋兩邊的間斷,可使裂紋橫穿網(wǎng)格,同時(shí)在裂尖引入漸近位移場(chǎng)函數(shù)更好地模擬奇異性［5-6］。上述2種方法一般推薦采用均勻布置的規(guī)則有限元網(wǎng)格對(duì)求解域進(jìn)行覆蓋,但均勻網(wǎng)格在計(jì)算量上通常不經(jīng)濟(jì),而不均勻網(wǎng)格的布置就仍存在大小網(wǎng)格的過(guò)渡問(wèn)題。

無(wú)網(wǎng)格法直接借助離散點(diǎn)可構(gòu)造出高階連續(xù)的整體近似函數(shù),在很大程度上克服了有限元法網(wǎng)格剖分的困難［7-8］。但其主要問(wèn)題是:沒(méi)有了網(wǎng)格,近似函數(shù)的構(gòu)造要比有限元法難度大,計(jì)算成本較高,且在積分上一般需要特殊處理［8］。

綜上所述,有限元法存在網(wǎng)格剖分困難、裂尖奇異區(qū)收斂慢等不足,無(wú)網(wǎng)格法存在實(shí)現(xiàn)相對(duì)復(fù)雜和計(jì)算量大等缺點(diǎn),廣義有限元存在的線性相關(guān)等問(wèn)題,成為制約其實(shí)現(xiàn)自動(dòng)計(jì)算的瓶頸。因此,筆者提出基于獨(dú)立覆蓋的新型數(shù)值流形方法,嘗試結(jié)構(gòu)靜力分析及其自動(dòng)計(jì)算的新途徑。

2 新型數(shù)值流形方法

1991年,石根華博士［9］首次將現(xiàn)代數(shù)學(xué)“流形”思想引入工程計(jì)算,發(fā)明了數(shù)值流形方法（以下簡(jiǎn)稱流形法),其中很關(guān)鍵的一點(diǎn)是:與有限元法將求解域離散成網(wǎng)格的方式不同,流形法構(gòu)造物理場(chǎng)近似解的所謂“數(shù)學(xué)網(wǎng)格”與實(shí)際求解域分離（后者只用于定義積分區(qū)域),要求數(shù)學(xué)網(wǎng)格在空間上完全覆蓋求解域。如圖1（a)所示,圖中著色的橢圓形物理區(qū)域被矩形數(shù)學(xué)網(wǎng)格完全覆蓋。

對(duì)應(yīng)于數(shù)學(xué)流形中的“覆蓋”,數(shù)學(xué)網(wǎng)格又被稱為數(shù)學(xué)覆蓋。數(shù)值計(jì)算中,在各覆蓋上定義覆蓋函數(shù)作為局部近似函數(shù)Vi,通過(guò)單位分解函數(shù)φi聯(lián)系

n起來(lái)形成整體近似函數(shù)V＝∑＝

φiVi（n為覆蓋

i1 數(shù))［9-10］。局部近似函數(shù)Vi一般為多項(xiàng)式級(jí)數(shù),當(dāng)然也可以是其它的完備級(jí)數(shù),甚至是解析解級(jí)數(shù)。

到目前為止,流形法的數(shù)學(xué)網(wǎng)格一般是有限元網(wǎng)格［11］,并且推薦采用均勻的規(guī)則網(wǎng)格布置,在有限元結(jié)點(diǎn)上定義級(jí)數(shù)作為局部近似函數(shù),單位分解函數(shù)就是有限單元的形函數(shù)。筆者稱這種采用有限元數(shù)學(xué)網(wǎng)格的流形法為“常規(guī)流形法”,該方法實(shí)際上與廣義有限元非常相似,也存在高階函數(shù)情況下的方程組線性相關(guān)等問(wèn)題。

從覆蓋形式上看,有限元網(wǎng)格是一種完全重疊的覆蓋,各覆蓋沒(méi)有自己的獨(dú)立區(qū)域。2012年以來(lái),筆者首次提出了部分重疊覆蓋的“新型流形法”［12-13］,引入了“獨(dú)立覆蓋”,即單位分解函數(shù)φi＝1、近似函數(shù)就是給定級(jí)數(shù)Vi的覆蓋區(qū)域。

首先研究了部分重疊的矩形覆蓋。如圖1（a)所示,標(biāo)識(shí)為a—p的大網(wǎng)格即矩形的獨(dú)立覆蓋,覆蓋之間只有較小的條形重疊區(qū)域。以圖1（b)、圖1 （c)中局部4個(gè)覆蓋為例,每個(gè)覆蓋包含1個(gè)獨(dú)立區(qū)域和幾個(gè)重疊區(qū)域,圖中的大矩形為獨(dú)立區(qū)域（即獨(dú)立覆蓋),P1—P4為相鄰2個(gè)覆蓋之間的重疊區(qū)域,P5為4個(gè)覆蓋共同的重疊區(qū)域。建議條形厚度取較小尺寸,其單位分解函數(shù)φi取有限元的形函數(shù)實(shí)現(xiàn)覆蓋之間的線性過(guò)渡［12,14］。獨(dú)立覆蓋還可以擴(kuò)展到任意形狀［15］。

筆者在文獻(xiàn)［14］中進(jìn)一步明確了新型流形法是基于獨(dú)立覆蓋的分析,并指出,獨(dú)立覆蓋具有任意形狀、任意連接及任意加密的優(yōu)良特性,突破了有限元網(wǎng)格的種種限制。如圖2所示,只需對(duì)獨(dú)立覆蓋進(jìn)行簡(jiǎn)單的再分塊就可以方便地實(shí)現(xiàn)覆蓋加密［13,16］,因此在h型自適應(yīng)上具有很大優(yōu)勢(shì),便于有針對(duì)性地加密覆蓋及進(jìn)行大小覆蓋的連接過(guò)渡,從而提高求解效率。新方法在獨(dú)立覆蓋上直接給出某一階的完備級(jí)數(shù)（如完全多項(xiàng)式),因此p型升階也很方便。

圖1 部分重疊的矩形數(shù)學(xué)覆蓋Fig.1 Schematic diagrams of rectangular mathematical covers with partially overlapping

圖2 矩形獨(dú)立覆蓋的加密Fig.2 Refinement of a rectangular independent cover

新型流形法的部分重疊覆蓋形式不僅解決了廣義有限元和常規(guī)流形法在高階情況下的線性相關(guān)問(wèn)題［17］,而且還能在局部區(qū)域采用特殊覆蓋函數(shù)以適合物理場(chǎng)局部特征。以裂縫分析為例,與擴(kuò)展有限元在裂尖處的單元結(jié)點(diǎn)上增加漸近位移場(chǎng)的部分特征函數(shù)相比,新型流形法采用Williams解析解級(jí)數(shù)的獨(dú)立覆蓋函數(shù)是對(duì)裂尖位移場(chǎng)的最佳逼近,因此收斂更快［13,18］。另外,如果剛度矩陣采用數(shù)值積分方式,那么在計(jì)算公式及其程序?qū)崿F(xiàn)上并不比有限元法復(fù)雜多少［12,18］,同樣的帶狀稀疏正定矩陣適用于大型計(jì)算且計(jì)算量可控。

考慮到真實(shí)場(chǎng)函數(shù)的分布未知,靠人工布置覆蓋及設(shè)定級(jí)數(shù)階次以獲得滿意的計(jì)算結(jié)果仍然是不太容易。因此采用自動(dòng)的前處理技術(shù)并結(jié)合有效的誤差估計(jì)進(jìn)行自適應(yīng)分析是非常必要的,這就是自動(dòng)計(jì)算的主要內(nèi)容。相比其它方法而言,正是由于新型流形法具有的諸多優(yōu)勢(shì),才使得在結(jié)構(gòu)靜力分析中盡早地實(shí)現(xiàn)自動(dòng)計(jì)算成為可能。

3 誤差估計(jì)初步研究

對(duì)當(dāng)前計(jì)算結(jié)果進(jìn)行精度分析的所謂“后驗(yàn)”誤差估計(jì)主要有2種方法［19］:計(jì)算局部區(qū)域平衡方程殘值的殘值法；基于數(shù)值解的后處理結(jié)果的誤差估計(jì),如O.C.Zienkiewicz和J.Z.Zhu提出的ZZ法。前者數(shù)學(xué)理論基礎(chǔ)比較完善,但由于形式復(fù)雜而沒(méi)有得到廣泛應(yīng)用；后者具有計(jì)算簡(jiǎn)單、易于理解的優(yōu)點(diǎn),受到工程界的廣泛歡迎。本文以后者為主。

設(shè)ε為精確應(yīng)變（對(duì)于常量彈性矩陣D,研究應(yīng)力或應(yīng)變效果相同),ε＾為數(shù)值計(jì)算得出的應(yīng)變,則差值為

能量范數(shù)誤差定義為［19］

精確的能量為［19］

真實(shí)應(yīng)變?chǔ)磐ǔＪ俏粗?改用一個(gè)相對(duì)準(zhǔn)確值ε*代替,則‖E‖和‖ΔE‖分別成為‖E*‖和‖ΔE*‖。定義誤差估計(jì)的有效性指數(shù)［19］為

其值趨向于1表明有漸近精確性。

相對(duì)誤差定義為［19］

關(guān)于相對(duì)準(zhǔn)確值ε*,本文研究了以下幾種誤差指標(biāo)。

3.1 條形處的應(yīng)變連續(xù)性指標(biāo)η1

有限元中的ZZ法對(duì)相鄰單元結(jié)點(diǎn)的應(yīng)力平均,或?qū)Χ鄠€(gè)單元組成的單元片進(jìn)行最小二乘的應(yīng)力平滑,作為相對(duì)準(zhǔn)確解。前者方法簡(jiǎn)單但精度不高,后者的精度隨單元片增大而提高,但計(jì)算量上升。

新型流形法在覆蓋的獨(dú)立區(qū)域已具備高階光滑的應(yīng)變,只是在覆蓋的條形重疊區(qū)域容易發(fā)生應(yīng)變的突變,是誤差最大的地方。受ZZ法的啟發(fā),本文以條形區(qū)域的應(yīng)變?yōu)橥黄瓶?當(dāng)2個(gè)覆蓋在其條形重疊區(qū)域的應(yīng)變相差不大時(shí),表明在這2個(gè)覆蓋區(qū)域的應(yīng)力基本連續(xù)平滑；當(dāng)更多的覆蓋達(dá)到此要求后,將使成片區(qū)域以至于整個(gè)求解域的應(yīng)力平滑,這樣就可以避免采用最小二乘法進(jìn)行應(yīng)力平滑的較大計(jì)算量。

圖3 兩個(gè)相鄰覆蓋及其條形區(qū)域的應(yīng)變Fig.3 Strains of two adjacent covers and strip area

以一維問(wèn)題的2個(gè)覆蓋為例,假設(shè)應(yīng)變分布如圖3所示。圖中左、右兩邊的黑色曲線分別表示覆蓋1和覆蓋2的應(yīng)變曲線ε1和ε2,灰線為

研究實(shí)例表明,隨著覆蓋函數(shù)階次的提高和覆蓋尺寸的減小,ε1和ε2兩條應(yīng)變曲線在條形連接處逐漸接近,因此將其接近程度作為應(yīng)變連續(xù)甚至平滑的判據(jù)是很直觀的想法?？蓪⒓t線作為相對(duì)準(zhǔn)確的應(yīng)變?chǔ)?,即Δε＝ε*-ε1或Δε＝ε*-ε2,兩者相同,然后按式（5)計(jì)算η1。

3.2 獨(dú)立覆蓋內(nèi)的自由表面應(yīng)力指標(biāo)η2

在具有自由表面的獨(dú)立覆蓋內(nèi),計(jì)算

3.3 獨(dú)立覆蓋內(nèi)的高階誤差指標(biāo)η3

參考文獻(xiàn)［20］,以一維為例,精確解u按近似解的級(jí)數(shù)uh展開(kāi)到m階,即

式中: Nn和an分別為級(jí)數(shù)的單個(gè)函數(shù)及其系數(shù),兩者的差值假設(shè)為級(jí)數(shù)的第m＋1階函數(shù),即

應(yīng)變誤差為

則一般情況下的二維和三維能量范數(shù)誤差

假定當(dāng)Nm＋1am＋1加入到近似解時(shí),原來(lái)參數(shù)a1,…,am保持不變。設(shè)RV為取到m項(xiàng)時(shí)將uh代入到微分方程后的殘值,再次應(yīng)用伽遼金法逼近am＋1,得到

將式（12)代入式（11)得到

然后按式（5)計(jì)算η3,其中ε*按當(dāng)前近似解uh計(jì)算。

該指標(biāo)表明了采用高階后對(duì)計(jì)算結(jié)果的改進(jìn),但其主要問(wèn)題是［20］:Nm＋1可能與積分域內(nèi)的殘差RV正交,得到誤差估計(jì)為0；假定a1,…,am保持不變,且新解與原來(lái)解之間的任何相互影響都忽略不計(jì),而這些條件在實(shí)際計(jì)算中很難完全達(dá)到。文獻(xiàn)［20］用一維算例說(shuō)明,這種方式得到的‖ΔE‖是誤差能量模的上限。

4 自動(dòng)計(jì)算初步研究

新型流形法的前處理包含2個(gè)步驟,其一是形成如圖1所示的矩形數(shù)學(xué)網(wǎng)格；其二是切割操作,即數(shù)學(xué)網(wǎng)格與物理邊界相互切割,生成網(wǎng)格內(nèi)的物理積分區(qū)域（又稱為流形單元)。需要強(qiáng)調(diào)的是,切割只涉及簡(jiǎn)單幾何運(yùn)算［21］,可完全自動(dòng)化,規(guī)則的矩形數(shù)學(xué)網(wǎng)格（含獨(dú)立覆蓋及其條形連接)的生成也很簡(jiǎn)單,因此前處理難度大大降低。目前,上述過(guò)程已編寫了計(jì)算機(jī)程序,只需人工輸入結(jié)構(gòu)的輪廓,就能自動(dòng)完成網(wǎng)格生成和切割操作,再結(jié)合人工輸入的材料參數(shù)和邊界條件,自動(dòng)生成計(jì)算數(shù)據(jù),實(shí)現(xiàn)完全自動(dòng)化的前處理。

新型流形法具有h型覆蓋加密和p型覆蓋函數(shù)升階的雙重便利性,理論上可通過(guò)h型和p型的適當(dāng)組合達(dá)到最高的計(jì)算收斂效率,這就是h-p型自適應(yīng)的最佳路徑。本文暫不考慮計(jì)算效率問(wèn)題,而是選擇一種路徑嘗試h-p型的混合自適應(yīng),其計(jì)算步驟為:

（1)在初始的矩形覆蓋網(wǎng)格上,首先進(jìn)行獨(dú)立覆蓋函數(shù)的p型升階——對(duì)于多項(xiàng)式級(jí)數(shù),從2階開(kāi)始,最高升至4階；對(duì)于裂尖Williams解析級(jí)數(shù), 從3階開(kāi)始,最高升至5階。

（2)如果精度不滿足,再進(jìn)行h型覆蓋加密,對(duì)于本文研究的矩形覆蓋,暫時(shí)采用最簡(jiǎn)單的二分法,如圖4所示,將圖4（a)的一個(gè)獨(dú)立覆蓋沿x向和y向分別均分為二,一共劃分為圖4（b)的4個(gè)小覆蓋。然后,重新進(jìn)行前處理。

（3)再回復(fù)到步驟（1)進(jìn)行p型升階,目前仍從最低階開(kāi)始計(jì)算。

循環(huán)往復(fù)直至所有區(qū)域達(dá)到設(shè)定的誤差指標(biāo)。但對(duì)于凹角、裂尖等具有奇異性的區(qū)域,網(wǎng)格越加密,應(yīng)力越大,不可能達(dá)到預(yù)設(shè)精度,因此設(shè)定網(wǎng)格加密的最多次數(shù)適時(shí)停止計(jì)算。

覆蓋加密和升階判斷中的誤差指標(biāo)是第3節(jié)中定義的η1,η2和η3,對(duì)于η1超標(biāo)時(shí)涉及的2個(gè)覆蓋都進(jìn)行升階或加密操作,如圖4（b)所示,當(dāng)右邊水平向條形的誤差超標(biāo)時(shí),加密成圖4（c)的形式。

以上自適應(yīng)過(guò)程由計(jì)算機(jī)自動(dòng)完成。再加上自動(dòng)化的前處理,就能初步實(shí)現(xiàn)自動(dòng)計(jì)算。

圖4 矩形覆蓋加密的二分法Fig.4 Bisection method for refinement of rectangular covers

5 結(jié)構(gòu)靜力分析的自動(dòng)計(jì)算算例

圖5 重力壩的覆蓋自動(dòng)加密過(guò)程Fig.5 Procedures of automatically refining covers in a gravity dam model

5.1 算例1——上游表面頂點(diǎn)受集中力作用的重力壩

如圖5所示的重力壩,壩高100 m,壩底長(zhǎng)60 m,壩頂寬度20 m,下游面折角處距壩頂20 m。壩體彈性模量為30 GPa,泊松比為0.167。上游表面頂點(diǎn)受集中力作用F＝1 000 kN。采用積分形式的罰函數(shù)法（即硬彈簧)實(shí)現(xiàn)底部的全約束。另采用稠密網(wǎng)格的有限元計(jì)算結(jié)果作為參考解。預(yù)設(shè)誤差指標(biāo)為:η1＝0.05,η2＝0.15,η3＝0.05。

不妨從1個(gè)覆蓋開(kāi)始計(jì)算。程序在誤差指標(biāo)的控制下自動(dòng)地進(jìn)行覆蓋函數(shù)升階及覆蓋加密,如圖5所示,加密的覆蓋逐漸趨向于集中力作用點(diǎn)、壩踵及下游壩面拐點(diǎn)等應(yīng)力奇異部位。最終得到的上游面應(yīng)力見(jiàn)圖6,圖6（a)的垂直應(yīng)力與細(xì)密網(wǎng)格的有限元參考解吻合很好,相對(duì)誤差均在2%以下, 圖6（b)的自由面應(yīng)力大都在0.001 MPa以下,與理論應(yīng)力值為0非常接近。

圖6 重力壩上游面應(yīng)力Fig.6 Stresses on the upstream surface of the gravity dam

5.2 算例2——方孔

中部開(kāi)小方孔的矩形板在上、下兩端受均布拉力p＝10 kN/ m2,采用如圖7所示的1/4計(jì)算模型,左邊和底部用積分形式的罰函數(shù)法施加法向約束。同樣從1個(gè)覆蓋開(kāi)始計(jì)算（當(dāng)然也可以事先自動(dòng)劃分較多的覆蓋)。從表1可見(jiàn),經(jīng)過(guò)5次自動(dòng)加密后,孔口測(cè)量點(diǎn)的應(yīng)力值與有限元細(xì)密網(wǎng)格參考解的誤差在3%以下。加密的覆蓋也是逐漸趨向于孔口的角點(diǎn)奇異部位,而在應(yīng)力變化平緩的區(qū)域則是自動(dòng)地采用大覆蓋和較低的多項(xiàng)式階次（2階居多),在計(jì)算量上體現(xiàn)出經(jīng)濟(jì)性。

圖7 方孔的1/4計(jì)算模型Fig.7 1/4 model of a rectangular plate with a small internal square hole

5.3 算例3——I型裂縫的應(yīng)力強(qiáng)度因子

圖8 兩端受均布拉力的長(zhǎng)條平板示意圖(含邊界裂紋)Fig.8 A long plate with a boundary crack subjected to a uniform tension

如圖8所示的無(wú)限長(zhǎng)條平板兩端受均布拉力P＝3 kN/ m,平板寬w＝0.4 m,含邊界裂紋長(zhǎng)a＝0.05 m。應(yīng)力強(qiáng)度因子KI的理論值1.45 kPa/ m1/ 2［22］。

如圖9所示,經(jīng)過(guò)圖9（a)至圖9（e)的幾次自動(dòng)加密后,裂尖應(yīng)力強(qiáng)度因子KI計(jì)算值分別為1.67, 1.40,1.52,1.46,1.47,最后與理論值1.45很接近。

圖9 計(jì)算應(yīng)力強(qiáng)度因子的流形元網(wǎng)格Fig.9 NMM meshes for calculated stress intensity factor

表1 孔口測(cè)量點(diǎn)切向應(yīng)力比較Table 1 Tangential stresses of the measuring points at the orifice　0.01 MPa

本例的特點(diǎn)是在裂尖所在的局部覆蓋中采用了Williams解析級(jí)數(shù)作為覆蓋函數(shù),即圖10中的裂尖解析覆蓋,因此收斂很快,計(jì)算精度高,而且應(yīng)力強(qiáng)度因子作為級(jí)數(shù)的系數(shù)直接得到。

圖10 圖9(d)中的裂尖附近局部網(wǎng)格Fig.10 Local NMM meshes around the crack tip in Fig.9(d)

此例加密覆蓋時(shí),為避免條形區(qū)域落到裂尖附近,對(duì)裂尖所在的覆蓋在y方向采用了三等份。但實(shí)際操作中不容易使裂尖位于獨(dú)立覆蓋的中心。如圖9（c)所示,當(dāng)裂尖靠近獨(dú)立覆蓋的邊緣時(shí),計(jì)算精度會(huì)下降。

另一個(gè)發(fā)現(xiàn)是,裂尖附近并非越加密越好。這是因?yàn)?越靠近裂尖,位移場(chǎng)和應(yīng)力場(chǎng)的變化梯度越大,與解析覆蓋相連的數(shù)值覆蓋就要求達(dá)到更高的階次。因此,對(duì)于裂尖等奇異區(qū),如何合理地劃分覆蓋及設(shè)定覆蓋函數(shù)階次的問(wèn)題需要進(jìn)一步研究。

6 結(jié) 語(yǔ)

本文采用基于矩形獨(dú)立覆蓋的新型流形法,初步實(shí)現(xiàn)了二維結(jié)構(gòu)線彈性靜力分析的自動(dòng)計(jì)算:只需人工輸入結(jié)構(gòu)外形、材料參數(shù)和邊界條件,其它過(guò)程完全交由計(jì)算機(jī)自動(dòng)完成,最終能得到滿足一定精度的計(jì)算結(jié)果。但限于研究時(shí)間較短,在h-p自適應(yīng)路徑及計(jì)算效率上沒(méi)有優(yōu)化,因此只是自動(dòng)計(jì)算的初步研究?？紤]到計(jì)算效率是關(guān)系到新方法是否具有實(shí)際應(yīng)用價(jià)值的重要問(wèn)題,下一步將從以下幾方面加以改進(jìn)。

（1)將矩形覆蓋擴(kuò)展到任意形狀覆蓋。本文基于矩形覆蓋進(jìn)行的自動(dòng)計(jì)算,可以直接推廣到三維長(zhǎng)方體覆蓋。但對(duì)于復(fù)雜求解域,矩形覆蓋很可能造成邊界附近各覆蓋內(nèi)的材料體大小不均,位移和應(yīng)力的復(fù)雜程度分布不均,使自適應(yīng)分析效率下降。因此,針對(duì)求解域形狀進(jìn)行覆蓋劃分是有必要的,這就要將矩形覆蓋擴(kuò)展到任意形狀覆蓋,使覆蓋貼合結(jié)構(gòu)外形并能模擬不同材料分區(qū),同時(shí)也能做到準(zhǔn)確施加本質(zhì)邊界條件。

（2)研究覆蓋的合理布置?？紤]到新型流形法的靈活性,總會(huì)有相對(duì)合適的覆蓋形狀、大小以及連接方式可使計(jì)算效率大幅提升。從本文幾個(gè)算例可知,裂尖、孔口角點(diǎn)等處具有應(yīng)力分布的奇異性,自適應(yīng)分析后期的很多過(guò)程都圍繞它們進(jìn)行。如果能夠根據(jù)物理場(chǎng)的先驗(yàn)知識(shí),研究出合理的覆蓋及函數(shù)形式在這些局部區(qū)域預(yù)先布置,就可以減少自適應(yīng)分析的步驟,大大提高計(jì)算效率。

（3)研究h-p自適應(yīng)的最佳路徑和提高計(jì)算速度的算法。新型流形法的計(jì)算量可控,理論上可以通過(guò)h-p自適應(yīng)的最佳路徑達(dá)到最小的計(jì)算量。同時(shí),自適應(yīng)分析需要反復(fù)計(jì)算,如何利用前次計(jì)算結(jié)果以加快當(dāng)前計(jì)算速度還有待研究。因此,新方法的自動(dòng)計(jì)算還有很大的效率提升空間。

（4)誤差估計(jì)理論有待進(jìn)一步完善。目前的誤差指標(biāo)與真實(shí)誤差的差別有多大（需要研究各指標(biāo)的有效性指數(shù))？是否對(duì)所有情況都有效？如果只關(guān)注局部區(qū)域的精度,如何確定局部誤差？這些問(wèn)題還有待解決。

最后需要指出,本文提出的自動(dòng)計(jì)算思想不僅是上述的自動(dòng)前處理和h-p型混合自適應(yīng),而且是貫穿整個(gè)計(jì)算過(guò)程的自動(dòng)化:采用可具有任意形狀、任意連接、任意加密的覆蓋網(wǎng)格,基于CAD精確幾何描述進(jìn)行自動(dòng)化的前處理；針對(duì)求解域形狀特征、根據(jù)物理場(chǎng)先驗(yàn)知識(shí)進(jìn)行自動(dòng)化的覆蓋布置和覆蓋函數(shù)選??；通過(guò)最佳的h-p型混合路徑進(jìn)行高效的自適應(yīng)重復(fù)計(jì)算；最后基于高精度應(yīng)力結(jié)果進(jìn)行自動(dòng)化的后處理。最終目標(biāo)是結(jié)構(gòu)分析自動(dòng)化CAE軟件。

致謝：部分重疊覆蓋的思想來(lái)自于石根華博士。筆者的研究一直受到石根華博士的指導(dǎo)，在此表示由衷的感謝。

［1］ 百度百科.CAE［EB/ OL］.（2015-06-11)［2015-07-10］. http:∥baike.baidu.com/ view/112169.htm.

［2］ ZIENKIEWICZ O C, TAYLOR R L. The Finite Element Method（The Fifth Edition)［M］. Oxford, Great Britain: Butterworth-Heinemann, 2000.

［3］ DUARTE C A, BABUSKA I, ODEN J T. Generalized Finite Element Methods for Three-dimensional Structural Mechanics Problems［J］. Computer＆Structures, 2000, 77:215-232.

［4］ 李錄賢,劉書靜,張慧華,等.廣義有限元方法研究進(jìn)展［J］.應(yīng)用力學(xué)學(xué)報(bào), 2009, 26（1): 96-108.

［5］ MOES N, DOLBOW J, BELYTSCHKO T. A Finite Element Method for Crack Growth without Remeshing［J］. International Journal for Numerical Methods in Engineering, 1999, 46: 131-150.

［6］ 李錄賢,王鐵軍.擴(kuò)展有限元法（XFEM)及其應(yīng)用［J］.力學(xué)進(jìn)展, 2005, 35（1): 5-20.

［7］ 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.

［8］ 張 雄,劉 巖,馬 上.無(wú)網(wǎng)格法的理論及應(yīng)用［J］.力學(xué)進(jìn)展, 2009, 39（1): 1-36.

［9］ SHI G H. Manifold Method of Material Analysis［C］∥U. S. Army Research Office. Transactions of the Ninth Army Conference on Applied Mathematics and Computing. Minneapolis, Minnesota, U S A, June 18-21,1991:57-76.

［10］BABUSKA I,MELENK J M.The Partition of Unity Method ［J］. International Journal for Numerical Methods in Engineering,1997, 40:727-758.

［11］張湘?zhèn)?章?tīng)?zhēng)榮,呂文閣,等.數(shù)值流形方法研究及應(yīng)用進(jìn)展［J］.力學(xué)進(jìn)展, 2010, 40（1): 1-12.

［12］祁勇峰,蘇海東,崔建華.部分重疊覆蓋的數(shù)值流形方法初步研究［J］.長(zhǎng)江科學(xué)院院報(bào), 2013, 30（1): 65-70.

［13］SU Hai-dong, QI Yong-feng, GONG Ya-qi, et al. Preliminary Research of Numerical Manifold Method Based on Partly Overlapping Rectangular Covers［C］∥DDA Commission of International Society for Rock Mechanics. Proceedings of the 11th International Conference on Analysis of Discontinuous Deformation（ICADD11).Fukuoka, Japan, August 27-29, 2013: 341-347.

［14］蘇海東,頡志強(qiáng),龔亞琦,等.基于獨(dú)立覆蓋的流形法的收斂性及覆蓋網(wǎng)格特性［J］.長(zhǎng)江科學(xué)院院報(bào), 2016,（2):131-136.

［15］蘇海東,祁勇峰,龔亞琦,等.任意形狀覆蓋的數(shù)值流形方法初步研究［J］.長(zhǎng)江科學(xué)院院報(bào),2013,30（12): 91-96.

［16］蘇海東,祁勇峰.部分重疊覆蓋流形法的覆蓋加密方法［J］.長(zhǎng)江科學(xué)院院報(bào), 2013, 30（7): 95-100.

［17］林紹忠,蘇海東.數(shù)值流形法獨(dú)立覆蓋區(qū)域的一種自動(dòng)選取方法［J］.長(zhǎng)江科學(xué)院院報(bào),2014,31（12):88-91.

［18］蘇海東,祁勇峰,龔亞琦.裂紋尖端解析解與周邊數(shù)值解聯(lián)合求解應(yīng)力強(qiáng)度因子［J］.長(zhǎng)江科學(xué)院院報(bào), 2013, 30（6): 83-89 .

［19］AKIN J E.Finite Element Analysis with Error Estimators ［M］. Oxford, Great Britain: Elseviewer Butterworth -Heinemann,2005.

［20］O C辛克維奇,K摩根.有限元與近似法［M］.陶振宗,張述良,周之德,譯.北京:人民交通出版社,1989.

［21］石根華.數(shù)值流形方法與非連續(xù)變形分析［M］.裴覺(jué)民譯.北京:清華大學(xué)出版社, 1997 .

［22］中國(guó)航空研究院.應(yīng)力強(qiáng)度因子手冊(cè)［M］.北京:科學(xué)出版社,1993.

（編輯:劉運(yùn)飛)

Preliminary Implementation of Automatic Computation for Static Analysis of Structures Using NMM Based on Independent Rectangular Covers

SU Hai-dong, GONG Ya-qi, XIE Zhi-qiang, QI Yong-feng
（Material and Engineering Structure Department,Yangtze River Scientific Research Institute, Wuhan 430010,China)

Abstract:By means of numerical manifold method（NMM) based on independent rectangular covers proposed in previous study, we present an automatic computation method for static analysis of linear-elastic structures, including automatic pre-processing, self-adaptive analyses and so on. According to the characteristics of independent covers, we give 3 indexes for posterior error such as index of strain continuity in strip area between two covers, stress index on boundary surfaces and high-order error in an independent cover. By using convenient h-version mesh refinement and p-version order increasing in the new method, we implement h-p version self-adaptivity in a selected way to realize h-version refinement of rectangular covers by using simple bisection method. Some 2D numerical examples are given to illustrate the feasibility of automatic computation, in which all the procedures are automatically accomplished by the computer, except for necessary manual input of structural outlines, material parameters, and boundary conditions. Finally, we obtain calculated data with certain precision.

Key words:numerical manifold method（NMM)；independent covers；automatic computation；error estimation；h-p hybrid adaptivity；static analysis of structures

作者簡(jiǎn)介：蘇海東（1968-),男,湖北武漢人,教授級(jí)高級(jí)工程師,博士,主要從事水工結(jié)構(gòu)數(shù)值分析工作和計(jì)算方法研究,（電話)027-82927167（電子信箱)suhd@ mail.crsri.cn。

基金項(xiàng)目：國(guó)家自然科學(xué)基金項(xiàng)目（51409012)；中央級(jí)公益性科研院所基本科研業(yè)務(wù)費(fèi)項(xiàng)目（CKSF2013031/ CL,CKSF2014054/ CL)

收稿日期：2015-10-30;修回日期：2015-11-25

doi:10.11988/ ckyyb.20150919

中圖分類號(hào)：TB115;TV311

文獻(xiàn)標(biāo)志碼：A

文章編號(hào)：1001-5485(2016)02-0144-07