在固定的獨立覆蓋網(wǎng)格中求解幾何非線性問題

蘇海東，董 鵬，頡志強

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

1 研究背景

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

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

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

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

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

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

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

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

2 獨立覆蓋流形法及其計算公式

2.1 獨立覆蓋流形法簡介

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

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

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

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

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

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

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

大、小覆蓋之間可任意錯位相連，因而能對覆蓋進行任意加密，比如圖2(c)所示的對塊體f進行細分和再細分。總之，覆蓋網(wǎng)格具有任意形狀(的積分區(qū)域)、任意(在條形處錯位)連接以及任意加密這3個“任意”特性，配合計算精度控制，有望實現(xiàn)自適應(yīng)分析甚至是無需人工參與的自動計算[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 流形法計算格式及其改進

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

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

(1)

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

(2)

(3)

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

(4)

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

e=BLu,

(5)

t+Δtσ=tσ+De，

(6)

(7)

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

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

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

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

(3)應(yīng)變計算完全按照小變形公式，未考慮應(yīng)變的二階項，在計算帶有旋轉(zhuǎn)的大變形或大位移問題時，步長稍大就會產(chǎn)生所謂的“體積膨脹”。對此改進如下：式(6)中的線性應(yīng)變增量e用Green-Lagrange應(yīng)變增量ε代替，即

ε=e+η。

(8)

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

(9)

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

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

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

(10)

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

V=φ1V1+φ2V2。

(11)

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

(12)

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

3 物質(zhì)點相對于固定空間的移動

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

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

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

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

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

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

(13)

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

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

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

4 材料體移動邊界的準確模擬

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

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

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

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

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

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

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

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

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

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

5 算 例

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

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

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

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

算例3：為了更好地說明本文方法相對于常規(guī)拉格朗日法的優(yōu)勢，給出如圖8(a)所示的半圓結(jié)構(gòu)在垂直體力作用下的變形分析算例。該結(jié)構(gòu)半徑為4 m，底部法向約束，取彈性模量為100 kN/m2，泊松比為0.2?？紤]對稱性，取1/4圓建立計算模型如圖8(b)所示，從圓弧的1/2處劃分為兩段，分別采用y=f(x)和x=f(y)2種樣條曲線，左側(cè)和底部為法向約束。采用了8×8個正方形網(wǎng)格覆蓋其可能的變形區(qū)域，以及4階多項式級數(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é)果對比。MARC模型劃分了540個有限元網(wǎng)格，考慮隨變形變化的荷載選項。從表4可見，在較小荷載作用時，本文計算值與MARC結(jié)果接近。隨著荷載增大，MARC的有限元網(wǎng)格因變形而扭曲，最終在f=17 kN/m3時出現(xiàn)計算不收斂的情況，提示矩陣奇異。而本文方法沒有網(wǎng)格扭曲現(xiàn)象，可一直計算下去，體現(xiàn)出方法的優(yōu)越性。

表4 A點水平位移與B點垂直位移Table 4 Horizontal displacements of point A andvertical displacements of point B

6 結(jié) 語

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

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

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

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