蘇海東，董 鵬，頡志強(qiáng)

(1.長(zhǎng)江科學(xué)院 材料與結(jié)構(gòu)研究所，武漢 430010； 2.水利部水工程安全與病害防治工程技術(shù)研究中心，武漢 430010)

1 研究背景

在計(jì)算力學(xué)中，當(dāng)材料發(fā)生的位移或變形相對(duì)于材料自身的幾何尺度不是很小時(shí)，就必須考慮由大位移或大變形造成的非線性，稱為幾何非線性問題。目前，相關(guān)的研究工作大都采用拉格朗日描述法或歐拉描述法[1](以下分別簡(jiǎn)稱拉格朗日法和歐拉法)來求解，前者主要用于固體力學(xué)，后者主要用于流體力學(xué)。

拉格朗日法(又稱為拉格朗日觀點(diǎn))，跟蹤材料體上的物質(zhì)點(diǎn)，計(jì)算網(wǎng)格附著于材料體而跟隨物質(zhì)點(diǎn)運(yùn)動(dòng)，其優(yōu)點(diǎn)在于：可得到物質(zhì)點(diǎn)的物理量時(shí)間歷程，有利于求解隨時(shí)間變化的材料非線性問題；能準(zhǔn)確描述材料體的邊界隨時(shí)間的運(yùn)動(dòng)；控制方程相對(duì)簡(jiǎn)單。但其缺陷是特大變形引起網(wǎng)格扭曲，使計(jì)算精度急劇下降，甚至導(dǎo)致計(jì)算失敗。常見的解決措施是重新劃分網(wǎng)格，但有限元網(wǎng)格劃分通常并不容易，而且需要通過插值運(yùn)算在新、舊網(wǎng)格之間傳遞信息，易引入額外誤差。

歐拉法(又稱為歐拉觀點(diǎn))著眼于空間點(diǎn)，關(guān)注空間點(diǎn)上的物理量隨時(shí)間的變化，材料體在固定于空間的網(wǎng)格中運(yùn)動(dòng)。其優(yōu)點(diǎn)是計(jì)算網(wǎng)格始終保持不變，不存在特大變形下的網(wǎng)格扭曲，缺陷是：不易得到物質(zhì)點(diǎn)的物理量時(shí)間歷程；難以準(zhǔn)確描述材料邊界隨時(shí)間的運(yùn)動(dòng)；控制方程存在遷移項(xiàng)(流體力學(xué)稱為“對(duì)流項(xiàng)”)，造成處理上的困難，比如系數(shù)矩陣非對(duì)稱，又比如基于歐拉觀點(diǎn)的流體Navier-Stokes方程由對(duì)流項(xiàng)引起的非線性和計(jì)算穩(wěn)定問題等。表1總結(jié)了拉格朗日法和歐拉法的特點(diǎn)，可以明顯看出，2種方法優(yōu)勢(shì)互補(bǔ)。

表1 拉格朗日法和歐拉法的比較Table 1 Comparison between the Lagrangian method and the Eularian method

1964年Noh[2]提出任意拉格朗日-歐拉(Arbitrary Lagrange-Euler，ALE)方法，被認(rèn)為吸收了拉格朗日法和歐拉法各自的優(yōu)點(diǎn)：通過規(guī)定合適的網(wǎng)格運(yùn)動(dòng)形式(在ALE控制方程中考慮網(wǎng)格運(yùn)動(dòng)的影響)，既可以準(zhǔn)確描述材料的移動(dòng)邊界，又能維持網(wǎng)格的合理形狀。但要根據(jù)研究對(duì)象來設(shè)計(jì)網(wǎng)格運(yùn)動(dòng)方式需要經(jīng)驗(yàn)和技巧，有時(shí)是非常困難的，會(huì)給使用者帶來負(fù)擔(dān)。

作為當(dāng)今計(jì)算力學(xué)的研究熱點(diǎn)，無網(wǎng)格法通過空間分布的離散點(diǎn)構(gòu)造近似函數(shù)，從而避免了網(wǎng)格劃分，其產(chǎn)生的一個(gè)重要背景就是針對(duì)有限元等網(wǎng)格類方法的網(wǎng)格扭曲問題[3]。無網(wǎng)格法中的物質(zhì)點(diǎn)法(Material Point Method，MPM)[4]采用拉格朗日型的質(zhì)點(diǎn)離散材料區(qū)域，用固定的歐拉背景網(wǎng)格求解動(dòng)量方程，集合了2種描述的優(yōu)勢(shì)，避免了網(wǎng)格扭曲和遷移項(xiàng)。但無網(wǎng)格法的基礎(chǔ)理論仍有待完善，計(jì)算量大的問題一直沒有很好地解決。

筆者于2011年提出固定網(wǎng)格流形法[5]：利用數(shù)值流形方法的網(wǎng)格與材料分離的特性，采用在空間均勻布置的正三角形網(wǎng)格和一階多項(xiàng)式覆蓋函數(shù)，并考慮材料相對(duì)于空間固定網(wǎng)格的移動(dòng)，初步實(shí)現(xiàn)了固定網(wǎng)格中的幾何非線性分析。該方法同樣集合了拉格朗日法和歐拉法的優(yōu)勢(shì)。然而，均勻網(wǎng)格不利于不同區(qū)域物理場(chǎng)的計(jì)算精度控制，而有限元網(wǎng)格通過結(jié)點(diǎn)相連的要求也給局部網(wǎng)格加密帶來不便，同時(shí)，常規(guī)數(shù)值流形方法還存在高階情況下的系統(tǒng)方程線性相關(guān)問題。

筆者后期提出獨(dú)立覆蓋流形法，具有網(wǎng)格任意連接和任意加密以及獨(dú)立覆蓋級(jí)數(shù)升階的優(yōu)點(diǎn)，且高階情況的系統(tǒng)方程線性無關(guān)。本文在固定的獨(dú)立覆蓋網(wǎng)格中研究幾何非線性問題，為下一步的計(jì)算精度控制和自適應(yīng)求解以及固定網(wǎng)格中的多重非線性分析打下基礎(chǔ)。為簡(jiǎn)便起見，本文僅研究單純的幾何非線性問題，不考慮由大應(yīng)變引起的本構(gòu)關(guān)系變化。

2 獨(dú)立覆蓋流形法及其計(jì)算公式

2.1 獨(dú)立覆蓋流形法簡(jiǎn)介

1991年，石根華先生[6-7]發(fā)明了數(shù)值流形方法(以下簡(jiǎn)稱流形法)，將現(xiàn)代數(shù)學(xué)的流形思想引入數(shù)值計(jì)算：如圖1(a)所示，一系列相互重疊的覆蓋用于分割求解域，使復(fù)雜的物理場(chǎng)在每個(gè)覆蓋簡(jiǎn)單化，以利于局部近似函數(shù)的逼近。

進(jìn)一步提出基于獨(dú)立覆蓋的數(shù)值流形方法，簡(jiǎn)稱獨(dú)立覆蓋流形法[10]，強(qiáng)調(diào)獨(dú)立覆蓋是計(jì)算分析的主要對(duì)象：如圖1(c)所示，獨(dú)立覆蓋占據(jù)一個(gè)覆蓋的主要面積，重點(diǎn)研究獨(dú)立覆蓋內(nèi)的級(jí)數(shù)，而覆蓋之間的重疊區(qū)域成為面積盡可能小(可只占覆蓋面積的1%)的窄“條形”，僅起連接各覆蓋級(jí)數(shù)的作用，并以有限單元的線性形函數(shù)作為單位分解函數(shù)。這樣，由各覆蓋分區(qū)的級(jí)數(shù)聯(lián)合形成的整體近似函數(shù)，通過加權(quán)殘值法(如伽遼金)逼近真實(shí)解，稱為基于流形思想的“分區(qū)級(jí)數(shù)解”，可作為微分方程的最終解答。

圖1 流形覆蓋示意圖Fig.1 Manifold covers

獨(dú)立覆蓋流形法的整體收斂是基于各覆蓋自身的局部收斂而建立的[10]。以通常采用的多項(xiàng)式級(jí)數(shù)為例，相當(dāng)于在各覆蓋區(qū)域內(nèi)將真實(shí)場(chǎng)函數(shù)展開為泰勒級(jí)數(shù)，或者說，各覆蓋內(nèi)的真實(shí)場(chǎng)函數(shù)(包括其導(dǎo)數(shù)，比如位移的導(dǎo)數(shù)是應(yīng)變或應(yīng)力)被多項(xiàng)式級(jí)數(shù)無限逼近。數(shù)學(xué)上的級(jí)數(shù)逼近理論保證了該方法具有嚴(yán)格的收斂性，且這種收斂性不受覆蓋內(nèi)實(shí)際占有的材料體(即真實(shí)的積分區(qū)域)形狀的影響，也不受相鄰覆蓋在條形重疊區(qū)域是否錯(cuò)位連接的影響。

暫不考慮覆蓋的條形重疊區(qū)域(將由程序自動(dòng)添加，見下文)，新方法可應(yīng)用2種覆蓋網(wǎng)格方式：①對(duì)求解域進(jìn)行任意形狀的塊體網(wǎng)格劃分[11]，不同大小及形狀的網(wǎng)格可隨意連接；②即本文所采用的規(guī)則網(wǎng)格方式，以圖2(a)為例，用矩形網(wǎng)格(最好是正方形網(wǎng)格)覆蓋求解域，通過切割操作將求解域分割成塊體(即網(wǎng)格內(nèi)的真實(shí)積分區(qū)域)。

基于獨(dú)立覆蓋流形法的收斂特性，材料不必占滿所在網(wǎng)格的全部面積，真正有意義的是材料體所占據(jù)的積分區(qū)域，通常在材料體邊界附近的網(wǎng)格內(nèi)為任意形狀。因此，網(wǎng)格不必與物理邊界相匹配(本質(zhì)邊界條件可采用罰函數(shù)法施加)，這是數(shù)值流形方法及其發(fā)展而來的獨(dú)立覆蓋流形法的重要特性之一，突破了現(xiàn)有的大多數(shù)數(shù)值方法在網(wǎng)格劃分上的限制，同時(shí)也是在固定網(wǎng)格中準(zhǔn)確描述材料移動(dòng)邊界的基礎(chǔ)。

當(dāng)然，積分方法也與常見的映射到標(biāo)準(zhǔn)母單元的方式不同：以圖2(a)右上角帶有一條曲線邊界的塊體c為例，如圖3(a)所示，將相鄰頂點(diǎn)構(gòu)成的各直線段依次與某一外點(diǎn)P0相連，形成有向單純形(二維為三角形)[7]，在單純形內(nèi)可采用數(shù)值積分[11]或解析積分[12]，然后考慮各單純形的正、負(fù)向后，累加得到整個(gè)塊體上的積分；圖3(b)的曲線和直線所圍區(qū)域，曲線邊界可直接按精確的曲線方程參與積分計(jì)算[11,13]，實(shí)現(xiàn)材料邊界的準(zhǔn)確模擬；然后如圖2(b)所示，以塊體e和f為例，在塊體之間自動(dòng)地添加條形，進(jìn)行條形區(qū)域的積分，同時(shí)扣除重復(fù)的塊體積分[11]。

大、小覆蓋之間可任意錯(cuò)位相連，因而能對(duì)覆蓋進(jìn)行任意加密，比如圖2(c)所示的對(duì)塊體f進(jìn)行細(xì)分和再細(xì)分?？傊?，覆蓋網(wǎng)格具有任意形狀(的積分區(qū)域)、任意(在條形處錯(cuò)位)連接以及任意加密這3個(gè)“任意”特性，配合計(jì)算精度控制，有望實(shí)現(xiàn)自適應(yīng)分析甚至是無需人工參與的自動(dòng)計(jì)算[14-15]。

圖2 規(guī)則網(wǎng)格覆蓋求解域Fig.2 Domain to be solved covered with regular meshes

圖3 包含曲線邊界的塊體的單純形劃分Fig.3 A block with a curve boundary divided intoseveral simplexes

2.2 流形法計(jì)算格式及其改進(jìn)

流形法采用分時(shí)步的計(jì)算格式，在t時(shí)刻計(jì)算完畢后更新材料體的邊界，按更新后的構(gòu)形計(jì)算t+Δt時(shí)刻的增量位移u，采用“慣性優(yōu)勢(shì)平衡方程式”[7]，即

(KL+Kg)u=Q+Fg-Fσ。

(1)

式中：KL為小變形情況下的線性剛度矩陣；Q為荷載(本文考慮依賴于變形狀態(tài)的所謂“跟隨荷載”)；當(dāng)考慮動(dòng)力荷載(不計(jì)阻尼)時(shí)，將慣性力分離出來，增加由此引起的剛度矩陣Kg及荷載Fg，即

(2)

(3)

Fσ稱為初應(yīng)力荷載,即

(4)

式中：ρ為密度(本文暫時(shí)忽略其變化)；v為速度；T為形函數(shù)矩陣；BL為線性應(yīng)變矩陣；σ為線彈性應(yīng)力;Ω為積分區(qū)域。每步計(jì)算完畢后的應(yīng)力和速度按如下公式更新：

e=BLu,

(5)

t+Δtσ=tσ+De，

(6)

(7)

式中：相同變量以左上標(biāo)區(qū)分其發(fā)生的構(gòu)形；e為線性應(yīng)變?cè)隽?；D為線彈性材料矩陣。

與計(jì)算幾何非線性問題常用的更新拉格朗日(Updated Lagrangian, UL)格式[1]相比，流形法計(jì)算格式的不同之處在于：

(1)缺少了幾何剛度矩陣(也稱為初應(yīng)力矩陣)。在非線性分析中，剛度矩陣沒必要很準(zhǔn)確，只需保證正確的收斂方向。UL格式更接近于切向剛度，有利于收斂的穩(wěn)定性，而流形法的剛度矩陣在初應(yīng)力較大時(shí)有可能偏離方向而出現(xiàn)計(jì)算振蕩，因此要引入慣性，使材料保持其原有的變形和運(yùn)動(dòng)方向。即使在靜力分析中，也推薦采用動(dòng)力松弛法[7]，即在式(7)中將速度tv乘上一個(gè)系數(shù)(如0.95)，由衰減的動(dòng)力方式計(jì)算靜力問題。

(2)沒有非線性平衡迭代過程。在步長(zhǎng)很小的情況下，不進(jìn)行迭代計(jì)算而直接將失衡力計(jì)入下一步荷載是允許的。但需要控制步長(zhǎng)，否則會(huì)引起計(jì)算結(jié)果的“漂移”。

(3)應(yīng)變計(jì)算完全按照小變形公式，未考慮應(yīng)變的二階項(xiàng)，在計(jì)算帶有旋轉(zhuǎn)的大變形或大位移問題時(shí)，步長(zhǎng)稍大就會(huì)產(chǎn)生所謂的“體積膨脹”。對(duì)此改進(jìn)如下：式(6)中的線性應(yīng)變?cè)隽縠用Green-Lagrange應(yīng)變?cè)隽喀糯?，?/p>

ε=e+η。

(8)

式中η為應(yīng)變二階項(xiàng)，其分量ηij用張量形式表示為[1]

(9)

式中：uk，i為k中的下標(biāo)“i”表示位移分量uk對(duì)獨(dú)立坐標(biāo)xi求偏導(dǎo)數(shù)；指標(biāo)k重復(fù)出現(xiàn)2次，表示在該指標(biāo)的取值范圍內(nèi)(如平面問題，k=1和2)遍歷求和。

(4)其他如應(yīng)力在不同構(gòu)形中的轉(zhuǎn)換、時(shí)間積分格式等差異還需進(jìn)一步研究。

獨(dú)立覆蓋流形法沿用了經(jīng)典流形法的數(shù)值計(jì)算格式，僅僅在數(shù)值解的空間離散上采用了“分區(qū)級(jí)數(shù)解”方案，表現(xiàn)為形函數(shù)T的不同形式：獨(dú)立覆蓋i中的位移u、速度v和應(yīng)力σ等均為級(jí)數(shù)表達(dá)，比如多項(xiàng)式級(jí)數(shù),即

(10)

式中：p為多項(xiàng)式的最高階次；ank為級(jí)數(shù)系數(shù)；(x0，y0)取為獨(dú)立覆蓋的中心坐標(biāo)；lx和ly分別表示x和y這2個(gè)方向的尺度。而在2個(gè)覆蓋之間的條形重疊區(qū)域內(nèi)，有

V=φ1V1+φ2V2。

(11)

其中，單位分解函數(shù)采用一維有限元的形函數(shù)[1]，即

(12)

式中ξ為在條形厚度方向定義的局部坐標(biāo)(以條形中點(diǎn)為原點(diǎn))。

3 物質(zhì)點(diǎn)相對(duì)于固定空間的移動(dòng)

表1給出了拉格朗日法和歐拉法的優(yōu)缺點(diǎn)。對(duì)于固體計(jì)算而言，拉格朗日法跟蹤物質(zhì)點(diǎn)、控制方程簡(jiǎn)單的優(yōu)點(diǎn)很有吸引力，而引入歐拉法的固定網(wǎng)格以解決其網(wǎng)格扭曲問題的關(guān)鍵在于兩點(diǎn)：①物質(zhì)點(diǎn)相對(duì)于固定網(wǎng)格(或固定空間點(diǎn))的移動(dòng)；②固定網(wǎng)格中的材料移動(dòng)邊界的準(zhǔn)確模擬。本節(jié)主要討論第①點(diǎn)。

如圖4所示，在t時(shí)刻計(jì)算完畢后，根據(jù)位移級(jí)數(shù)解更新所有的材料點(diǎn)至t+Δt時(shí)刻的新位置，包括更新材料體的原邊界(虛線)至新的構(gòu)形邊界(實(shí)線)。仍采用2.2節(jié)的拉格朗日型計(jì)算格式，考慮到積分區(qū)域可為任意形狀的特性，在確定的形狀下(即實(shí)線邊界內(nèi))可以做到對(duì)式(1)至式(4)中的各矩陣和向量積分的準(zhǔn)確計(jì)算，其中的速度v和應(yīng)力σ的確定方式如下所述。

圖4 逆向追蹤示意圖Fig.4 Schematic ofbackward tracking

拉格朗日方程關(guān)注的是物質(zhì)點(diǎn)，關(guān)鍵是如何得到向量Fg和Fσ中關(guān)于物質(zhì)點(diǎn)的速度v和應(yīng)力σ。從有限元法來看，網(wǎng)格附著在材料體上并隨材料移動(dòng)，移動(dòng)后的網(wǎng)格結(jié)點(diǎn)攜帶著物理量，通過插值方式可獲得網(wǎng)格內(nèi)任意點(diǎn)的物理量(不考慮網(wǎng)格扭曲造成的精度下降)。如果網(wǎng)格不動(dòng)，那么固定結(jié)點(diǎn)上的物理量需要修正，以反映材料相對(duì)于固定網(wǎng)格(及其所代表的固定空間)的移動(dòng)，而有限元法不易做到準(zhǔn)確修正，特別是應(yīng)力。

獨(dú)立覆蓋流形法采用的“級(jí)數(shù)解”具有公式化的表達(dá)，為物理量的準(zhǔn)確修正創(chuàng)造了條件。本文借鑒流體計(jì)算中的半拉格朗日法[16]，提出對(duì)物質(zhì)點(diǎn)進(jìn)行“逆向追蹤”的思路：如圖4所示，當(dāng)前構(gòu)形中坐標(biāo)為(t+Δtx,t+Δty)的空間點(diǎn)A被某物質(zhì)點(diǎn)占據(jù)，可通過上一時(shí)步得到的位移級(jí)數(shù)解逆向求解該物質(zhì)點(diǎn)在上一時(shí)步的位置A′(tx,ty)，即

t+Δtx=tx+ux(tx,ty) ,t+Δty=ty+uy(tx,ty) 。

(13)

式中ux、uy分別為x、y兩個(gè)方向的位移級(jí)數(shù)，高階情況下難以顯式求解，可采用直接迭代法，以當(dāng)前位置為初值，在控制好步長(zhǎng)及整體計(jì)算穩(wěn)定的情況下，一般經(jīng)過幾次迭代就可得到(tx,ty)的收斂解。以此坐標(biāo)代入上一時(shí)步的應(yīng)力和速度的級(jí)數(shù)表達(dá)式求出σ和v，作為當(dāng)前位置A的物質(zhì)點(diǎn)所攜帶的物理量初值，代入式(3)和式(4)中的Fg和Fσ參與積分，從而實(shí)現(xiàn)當(dāng)前時(shí)步的拉格朗日方程的準(zhǔn)確計(jì)算。

以歐拉觀點(diǎn)對(duì)上述過程給出另一種解釋：歐拉描述關(guān)注固定空間點(diǎn)上的物理量的變化，即關(guān)注不同時(shí)刻經(jīng)過該空間點(diǎn)的不同物質(zhì)點(diǎn)所攜帶的物理量的變化，而遷移項(xiàng)(對(duì)流項(xiàng))所反映的正是這種由于物質(zhì)從一點(diǎn)移動(dòng)到另一點(diǎn)的空間不均勻性所引起的變化。通過“逆向追蹤”手段反映這種變化，就可消除遷移項(xiàng)，實(shí)現(xiàn)材料相對(duì)于固定空間點(diǎn)移動(dòng)的準(zhǔn)確描述。

以上討論的一般性情況，同樣適合于各獨(dú)立覆蓋。不過，針對(duì)各物質(zhì)點(diǎn)的操作，需要形成級(jí)數(shù)表達(dá)式，以便于在條形區(qū)域內(nèi)“加權(quán)”，以及式(6)中的應(yīng)力累積和式(7)中的速度更新：在“逆向追蹤”后，根據(jù)各積分點(diǎn)的速度v和應(yīng)力σ的修正數(shù)值，用最小二乘法得到修正后的級(jí)數(shù)表達(dá)式。另外，考慮應(yīng)變二階項(xiàng)計(jì)算應(yīng)變及應(yīng)力時(shí)，由于式(9)的復(fù)雜性，也要通過最小二乘法獲得與位移同階次的應(yīng)力級(jí)數(shù)，然后再代入式(6)與tσ累積得到t+Δtσ的級(jí)數(shù)表達(dá)式。

4 材料體移動(dòng)邊界的準(zhǔn)確模擬

對(duì)于變形過程中可能出現(xiàn)的復(fù)雜形狀邊界，本文采用3次樣條曲線描述，按單值性分為y=f(x)和x=f(y)2種，并可設(shè)置多條曲線用于描述復(fù)雜邊界。構(gòu)形變化時(shí)，更新各樣條曲線的控制點(diǎn)。

圖5 正方形網(wǎng)格中的積分區(qū)域Fig.5 Integral region in asquare mesh

采用正方形網(wǎng)格，與每個(gè)時(shí)步更新后的邊界輪廓進(jìn)行切割操作。本文提出形成固定網(wǎng)格內(nèi)的積分區(qū)域的新方法，如圖5所示，簡(jiǎn)要說明如下：

對(duì)于某個(gè)正方形網(wǎng)格，邊界輪廓(包括直線段或樣條曲線)依次與網(wǎng)格的4條邊求交，比如，得到圖5中左、右兩條邊上的交點(diǎn)A、B，作為有效頂點(diǎn)；通過判斷，結(jié)點(diǎn)1和結(jié)點(diǎn)2在域內(nèi)，也作為有效頂點(diǎn)；各段邊界線的起點(diǎn)或終點(diǎn)也有可能落在網(wǎng)格內(nèi)，成為有效頂點(diǎn)；按逆時(shí)針走向，從邊1-2開始，在各邊(及網(wǎng)格內(nèi)部)依次連接有效頂點(diǎn)，形成積分區(qū)域1-2-B-A-1；其他與邊界不相交的正方形網(wǎng)格，如判斷其形心在域內(nèi)，則整個(gè)網(wǎng)格作為積分區(qū)域參與積分。

在材料邊界的移動(dòng)過程中，不可避免地會(huì)形成在網(wǎng)格內(nèi)只占很小區(qū)域的所謂“小塊”，如圖6(a)所示的右上角網(wǎng)格。其中，在材料進(jìn)入新網(wǎng)格時(shí)，一定會(huì)首先出現(xiàn)小塊。設(shè)置一定的小塊比例，比如，小塊的面積為所有塊體平均面積的5%以下。

小塊的存在會(huì)影響方程性態(tài)，造成數(shù)值解的不穩(wěn)定?？紤]到獨(dú)立覆蓋流形法具有在不同大小的網(wǎng)格之間錯(cuò)位連接的特性，采用的處理方式如圖6(b)所示，將小塊所在的網(wǎng)格與其周邊大網(wǎng)格合并(目前暫時(shí)取小塊的長(zhǎng)邊所連接的大網(wǎng)格)，應(yīng)力和速度采用大網(wǎng)格中的級(jí)數(shù)。

隨著材料邊界的移動(dòng)，小塊有兩種變化的可能性：第1種如圖6(c)所示，邊界向下移動(dòng)，小塊完全消失；第2種如圖6(d)所示，邊界向上移動(dòng)，小塊逐漸占據(jù)一定比例的正方形網(wǎng)格面積(比如大于平均塊體面積的5%)，則成為獨(dú)立網(wǎng)格中的積分區(qū)域參與計(jì)算，其初始的應(yīng)力和速度繼承原大網(wǎng)格中的級(jí)數(shù)。第2種情況同時(shí)說明了材料進(jìn)入新網(wǎng)格時(shí)，如何通過小塊將信息傳遞給新網(wǎng)格。

圖6 小塊的處理Fig.6 Processing of small blocks

通過以上的各種處理，在幾何非線性問題的分步計(jì)算中，只要在各覆蓋分區(qū)內(nèi)采用足夠的多項(xiàng)式級(jí)數(shù)逼近，就能保證每步計(jì)算的每個(gè)物質(zhì)點(diǎn)(包括邊界點(diǎn))的位移精度，從而保證邊界輪廓的準(zhǔn)確更新、速度的準(zhǔn)確更新和應(yīng)力的準(zhǔn)確累積，在各時(shí)步實(shí)現(xiàn)移動(dòng)邊界的準(zhǔn)確描述。

圖7 固定的正方形網(wǎng)格中的懸臂梁變形Fig.7 Deformations of cantilever beam in fixedsquare meshes

5 算 例

算例1：計(jì)算如圖7所示的水平放置的懸臂梁，左端固定。該梁長(zhǎng)10 m，截面的高度和寬度均為1 m，彈性模量E為3×105kN/m2，泊松比μ為0.2。在梁自由端的截面中點(diǎn)作用始終垂直向下的集中力P。在固定的正方形網(wǎng)格中計(jì)算其大變形，梁的上、下邊緣分別采用3次樣條曲線。每個(gè)獨(dú)立覆蓋采用3階多項(xiàng)式級(jí)數(shù)，步長(zhǎng)為6.25 kN(數(shù)值試驗(yàn)表明，計(jì)算精度隨步長(zhǎng)減小而提高)。典型時(shí)步的自由端中點(diǎn)位移及其與理論解的比較見表2，相對(duì)誤差在1%左右。在計(jì)算中發(fā)現(xiàn)，若完全按照靜力計(jì)算，則在很大變形時(shí)(P>1 000 kN)會(huì)出現(xiàn)計(jì)算振蕩現(xiàn)象。除了引入初應(yīng)力剛度矩陣解決此問題以外，本文采用動(dòng)力松弛法求解靜力問題，也能保證計(jì)算的穩(wěn)定性。

表2 懸臂梁自由端的中點(diǎn)位移Table 2 Displacements of the mid-point at the free endof cantilever beam

算例2：圖7中的正方形網(wǎng)格內(nèi)的梁作為剛性桿件，彈性模量取大值以模擬剛性。桿件的左端中點(diǎn)為固定約束，從原位置開始旋轉(zhuǎn)，初始角速度為1 rad/s，步長(zhǎng)為0.01 s。各時(shí)步的自由端中點(diǎn)坐標(biāo)見表3，采用1階和4階多項(xiàng)式級(jí)數(shù)的計(jì)算結(jié)果相同，但是否考慮二階應(yīng)變對(duì)計(jì)算結(jié)果有影響：不計(jì)二階應(yīng)變，則桿件會(huì)在旋轉(zhuǎn)過程中逐漸“拉長(zhǎng)”；考慮二階應(yīng)變，得到與理論解幾乎一致的結(jié)果。此算例驗(yàn)證了本文方法模擬剛體運(yùn)動(dòng)的準(zhǔn)確性。

表3 勻速旋轉(zhuǎn)桿件的自由端中點(diǎn)坐標(biāo)Table 3 Coordinates of the mid-point at the freeend of rod with uniform rotation

算例3：為了更好地說明本文方法相對(duì)于常規(guī)拉格朗日法的優(yōu)勢(shì)，給出如圖8(a)所示的半圓結(jié)構(gòu)在垂直體力作用下的變形分析算例。該結(jié)構(gòu)半徑為4 m，底部法向約束，取彈性模量為100 kN/m2，泊松比為0.2?？紤]對(duì)稱性，取1/4圓建立計(jì)算模型如圖8(b)所示，從圓弧的1/2處劃分為兩段，分別采用y=f(x)和x=f(y)2種樣條曲線，左側(cè)和底部為法向約束。采用了8×8個(gè)正方形網(wǎng)格覆蓋其可能的變形區(qū)域，以及4階多項(xiàng)式級(jí)數(shù)。

圖8 在垂直體力作用下的半圓結(jié)構(gòu)及其固定網(wǎng)格Fig.8 Semicircle in fixed meshes under verticalbody loads

圖9 在不同垂直體力作用下的結(jié)構(gòu)變形Fig.9 Deformations of the semicircle under variedvertical body loads

不同體力荷載f作用下的結(jié)構(gòu)變形如圖9所示，結(jié)構(gòu)頂部豎向位移(VB)和底部外邊緣的水平向位移(UA)列入表4，并與大型有限元分析軟件MARC的結(jié)果對(duì)比。MARC模型劃分了540個(gè)有限元網(wǎng)格，考慮隨變形變化的荷載選項(xiàng)。從表4可見，在較小荷載作用時(shí)，本文計(jì)算值與MARC結(jié)果接近。隨著荷載增大，MARC的有限元網(wǎng)格因變形而扭曲，最終在f=17 kN/m3時(shí)出現(xiàn)計(jì)算不收斂的情況，提示矩陣奇異。而本文方法沒有網(wǎng)格扭曲現(xiàn)象，可一直計(jì)算下去，體現(xiàn)出方法的優(yōu)越性。

表4 A點(diǎn)水平位移與B點(diǎn)垂直位移Table 4 Horizontal displacements of point A andvertical displacements of point B

6 結(jié) 語

本文提出應(yīng)用級(jí)數(shù)對(duì)物質(zhì)點(diǎn)進(jìn)行“逆向追蹤”的思路，在空間固定的網(wǎng)格中采用拉格朗日控制方程求解幾何非線性問題，集合了拉格朗日法的跟蹤物質(zhì)點(diǎn)、控制方程簡(jiǎn)單、邊界描述準(zhǔn)確以及歐拉法的網(wǎng)格無扭曲的優(yōu)勢(shì)，避開了兩種方法相應(yīng)的缺陷。相比之下，任意拉格朗日和歐拉法需要設(shè)計(jì)網(wǎng)格運(yùn)動(dòng)算法，該算法通用性不強(qiáng)且有可能很復(fù)雜。本文采用級(jí)數(shù)公式反映物質(zhì)點(diǎn)的相對(duì)運(yùn)動(dòng)，相當(dāng)于公式化的物質(zhì)點(diǎn)法，求解精度與獨(dú)立覆蓋流形法的常規(guī)計(jì)算無異，而物質(zhì)點(diǎn)法的不足是用移動(dòng)的質(zhì)點(diǎn)作為積分點(diǎn)，相比常規(guī)有限元法的精度降低。本文的研究表明，除了無網(wǎng)格法的途徑之外，網(wǎng)格類的方法也可以很好地解決特大變形下的網(wǎng)格扭曲問題。

與物質(zhì)點(diǎn)法類似，采用固定網(wǎng)格會(huì)涉及物質(zhì)點(diǎn)在相鄰網(wǎng)格之間的穿越，這要求物理量在網(wǎng)格之間具有一定的連續(xù)性，特別是應(yīng)力，否則會(huì)引起局部計(jì)算的振蕩。因此，為保證每步計(jì)算的收斂性(包括應(yīng)力的連續(xù)性)，本文算例多采用3至4階多項(xiàng)式級(jí)數(shù)。相對(duì)于筆者前期提出的固定網(wǎng)格流形法而言，本文基于獨(dú)立覆蓋的固定網(wǎng)格具有任意連接和任意加密的特性，且覆蓋級(jí)數(shù)升階方便。目前正在引入自適應(yīng)分析技術(shù)，在一定的精度控制下，通過加密網(wǎng)格(也可配合級(jí)數(shù)升階)，可以實(shí)現(xiàn)每步的應(yīng)力“分區(qū)級(jí)數(shù)解”在一定程度上的整體連續(xù)甚至光滑。因此固定網(wǎng)格并非一成不變，可在各積分點(diǎn)處由“整體連續(xù)”的級(jí)數(shù)公式計(jì)算應(yīng)力和速度初值，從而實(shí)現(xiàn)在各時(shí)步的不同網(wǎng)格之間準(zhǔn)確傳遞信息，以替代通?；诰W(wǎng)格的傳遞方式。

本文方法跟蹤物質(zhì)點(diǎn)且邊界描述清晰準(zhǔn)確，便于求解隨時(shí)間變化的材料非線性和接觸非線性問題：通過“逆向追蹤”，使得當(dāng)前位置的物質(zhì)點(diǎn)攜帶材料狀態(tài)信息，其中當(dāng)前邊界上的物質(zhì)點(diǎn)還攜帶邊界接觸狀況。因此，將來有可能在空間固定的背景網(wǎng)格中同時(shí)進(jìn)行幾何、材料和接觸的多重非線性分析。

致謝：感謝石根華先生的指導(dǎo)。感謝中央級(jí)公益性科研院所基本科研業(yè)務(wù)費(fèi)項(xiàng)目(CKSF2019394/GC)的資助。