蘇海東

?

高階數(shù)值流形方法的速度公式

蘇海東

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

高階數(shù)值流形方法可以顯著提高結(jié)構(gòu)計算精度，但目前在涉及大位移的動力分析中往往得到精度很差、甚至不正確的速度結(jié)果。基于平面三角形數(shù)學(xué)網(wǎng)格和一階多項式覆蓋函數(shù)，通過一個剛體桿件旋轉(zhuǎn)算例探討其中的原因，得出必須考慮構(gòu)形坐標變化對速度的影響，并提出高階流形法的3種速度處理方法及相應(yīng)的高階速度公式。該方法對一些在結(jié)點處增加廣義自由度的類似方法（如廣義有限元）的幾何非線性問題分析也具有一定的參考價值。

數(shù)值流形方法；高階多項式覆蓋函數(shù)；大位移；速度公式；廣義自由度

doi：10．11988/ckyyb．20150343

1　研究背景

作為計算力學(xué)的經(jīng)典和普遍使用的方法，有限元法廣泛應(yīng)用于實際工程計算和科學(xué)研究中已長達數(shù)十年，其特點是以待求的場變量作為網(wǎng)格結(jié)點上的未知數(shù)，在結(jié)構(gòu)應(yīng)力分析中表現(xiàn)為位移自由度。近年來在有限元法的基礎(chǔ)上，數(shù)值流形方法［1］、擴展有限元［2］、廣義有限元［3］等新方法，通過在結(jié)點處增加廣義自由度來提高計算精度，在結(jié)構(gòu)的小變形應(yīng)力分析中取得了滿意的效果［4-6］。

在上述新方法中，數(shù)值流形方法率先進入了大位移或大變形的幾何非線性問題研究領(lǐng)域。研究表明，在結(jié)點上引入除常數(shù)位移場變量以外的多項式函數(shù)后（被稱為高階數(shù)值流形方法），這種廣義自由度對計算結(jié)果的影響逐漸顯現(xiàn)出來。作者在高階數(shù)值流形法的大變形靜力分析中曾指出［7］，有必要在應(yīng)力累積中考慮構(gòu)形坐標的變化，并相應(yīng)地提出了高階初應(yīng)力公式，通過懸臂梁大變形靜力分析的計算結(jié)果得到了驗證。

進一步討論，如果分析大位移條件下的高階數(shù)值流形法的動力問題，在速度計算中是否也需要進行特殊處理呢？對于此問題，目前的研究均未涉及到，本文通過一個簡單例子展開討論。

為使讀者對數(shù)值流形方法中的廣義自由度有所了解，以下首先簡要介紹高階數(shù)值流形方法。

2　高階數(shù)值流形方法簡介

數(shù)值流形方法［1］（以下簡稱流形法）是我國留美學(xué)者石根華博士利用現(xiàn)代數(shù)學(xué)中流形分析的有限覆蓋技術(shù)建立起來的一種新的數(shù)值分析方法，具有廣闊的發(fā)展前景。流形法引入2套覆蓋網(wǎng)格：一是用于構(gòu)造物理場近似解的數(shù)學(xué)網(wǎng)格；另一個是表示材料邊界、用于定義積分區(qū)域的物理網(wǎng)格。2套網(wǎng)格相互獨立，只要求前者在空間上完全包容后者。

流形法對各個覆蓋獨立定義覆蓋函數(shù)，通過聯(lián)絡(luò)函數(shù)（又稱權(quán)函數(shù)）將各個獨立覆蓋函數(shù)連起來形成總體的覆蓋函數(shù)。規(guī)定覆蓋Ui上的局部覆蓋函數(shù)為ui（x）（x∈Ui），它可以是多項式級數(shù)或其他級數(shù)，隨著級數(shù)項數(shù)的增加及階數(shù)的提高，計算精度得以提高。這些覆蓋函數(shù)用權(quán)函數(shù)聯(lián)系起來，其含義是加權(quán)平均，權(quán)函數(shù)滿足

因此總體覆蓋函數(shù)為

式中n為覆蓋總數(shù)。

現(xiàn)有的流形法多采用有限元網(wǎng)格來定義數(shù)學(xué)網(wǎng)格，與有限元任一結(jié)點相連的所有網(wǎng)格形成一個數(shù)學(xué)覆蓋，局部覆蓋函數(shù)定義于結(jié)點上形成結(jié)點覆蓋。這樣，任一原始的有限元網(wǎng)格就是各結(jié)點覆蓋的重疊部分，在此重疊區(qū)域（即有限元網(wǎng)格）內(nèi)，覆蓋的權(quán)函數(shù)取為有限單元的形函數(shù)。本文以三角形數(shù)學(xué)網(wǎng)格為例，其形函數(shù)為面積坐標表示的一階多項式Li（下文統(tǒng)一用上標表示結(jié)點序號，注意與指數(shù)相區(qū)別，用下標n-1，n，n+1表示構(gòu)形序號）。

考慮在各個結(jié)點覆蓋上采用一階多項式的局部覆蓋函數(shù)，即

式中的a至f為多項式系數(shù)，即待求未知數(shù)，可以看出，這些未知數(shù)已沒有通常意義下的位移自由度的物理含義，因此被稱為廣義自由度。

三角形網(wǎng)格內(nèi)的位移函數(shù)如式（3）所示，結(jié)點覆蓋函數(shù)用三角形網(wǎng)格的形函數(shù)連接成為二階多項式函數(shù)，即

式中Te為網(wǎng)格內(nèi)的插值函數(shù)（三角形網(wǎng)格的形函數(shù)與各單項式函數(shù)1，x，y之乘積）為以ai至fi表示的網(wǎng)格自由度向量。

與有限元法一樣，流形法通過最小勢能原理建立平衡方程來進行材料的變形和應(yīng)力分析，不同之處僅在于位移近似函數(shù)的表示方式。文獻［8］指出，流形法的動力分析采用了經(jīng)典的NewMark公式，只是將參數(shù)定為β=1/2和γ=1。其動力影響反映在慣性力引起的剛度矩陣（即通常所指的質(zhì)量矩陣）和荷載向量上，后者包含速度［1］。第n步計算完畢后，n+1步的初始速度為［1］

式中：Vn為第n步的速度；Δ un為第n步的增量位移；Δt為時間步長。

對于采用一階多項式覆蓋函數(shù)的高階流形法而言，其結(jié)點的速度覆蓋按式（4）計算后也具有式（2）的形式，同時，網(wǎng)格內(nèi)的速度分布具有式（3）的形式，只是系數(shù)不同。

高階流形法由于在結(jié)點處引入了帶有廣義自由度的多項式覆蓋函數(shù)，因而構(gòu)形坐標的變化會對材料體的特征量（如應(yīng)力、速度）產(chǎn)生影響。作者已在文獻［7］中對高階流形法的初應(yīng)力公式進行了修正，取得了很好的效果。以下針對速度問題，以剛性桿件的勻速旋轉(zhuǎn)為例進行討論。

3　高階流形法速度公式存在的問題

如圖1所示，采用二維流形元網(wǎng)格計算剛性桿件，桿件長10 m，矩形截面的高度和寬度均為1 m。按平面應(yīng)力計算，彈性模量取大值以模擬剛性。桿件的左端中點為固定約束，從圖示的位置開始按1 rad/s的角速度順時針旋轉(zhuǎn)，計算步長為0．01 s。

圖1　剛性桿件及其流形元網(wǎng)格Fig．1　A rigid bar and its NMM meshes

在文獻［9］編制的靜力分析程序的基礎(chǔ)上增加動力計算部分，分別取覆蓋函數(shù)為0階（即常規(guī)的位移自由度）和式（2）的一階多項式，計算圖1的桿件轉(zhuǎn)動，得到的桿件自由端中點的坐標如表1所示。

表1　原始方法的桿件自由端中點坐標Table 1　Coordinates of the middle point at the free end of bar（conventional method）

表1中可見，0階計算結(jié)果與理論值（x= 10cos（t），y=-10sin（t））很接近，說明對于剛體旋轉(zhuǎn)運動，用0階計算就已足夠準確。而在同樣條件下，一階計算結(jié)果與理論值沒有可比性，旋轉(zhuǎn)速度明顯快于理論值。這表明，當引入一階廣義自由度時，高階流形法的速度存在問題。

為便于分析，下面以一維問題為例進行討論。見圖2所示的一個一維網(wǎng)格。在第n-1步時結(jié)點1的初始位置為x1n-1，結(jié)點2的初始位置為x2n-1，網(wǎng)格內(nèi)任一點的坐標為xn-1。位移計算完畢后，按式（4）計算得到的結(jié)點1和結(jié)點2的速度覆蓋函數(shù)即第n步的初始速度）分別為：

式中a-1n，a-2n，b-1n，b-2n為構(gòu)形變化前的多項式系數(shù)，則一維網(wǎng)格內(nèi)任意點的速度為

式中N1，N2分別為對應(yīng)于結(jié)點1和結(jié)點2的一維網(wǎng)格形函數(shù)。

圖2　從第n-1步移至第n步的一維網(wǎng)格Fig．2　Schematic diagram of one-dimensional mesh moving from the（n-1）thstep to the nthstep

如圖2所示，在進行第n步計算前，根據(jù)計算得到的位移量Δu將材料體從第n-1構(gòu)形更新至第n構(gòu)形，附著在材料體上的數(shù)學(xué)網(wǎng)格結(jié)點隨材料體移動至新的結(jié)點位置x1n和x2n，網(wǎng)格內(nèi)的任一點由xn-1移至xn。若按原始的高階流形法速度公式，結(jié)點1和結(jié)點2的第n步速度覆蓋函數(shù)分別為其中，在構(gòu)形更新過程中的多項式系數(shù)保持不變，即?？紤]到剛體運動中的形函數(shù)N1和

作為第n步計算的初始速度。然而表1表明，這種速度處理是不正確的。

實際上，材料的構(gòu)形坐標發(fā)生了變化，即

從以上分析可知，由坐標移動產(chǎn)生的速度誤差正是造成表1中的一階速度計算結(jié)果不正確的原因，因此需要對高階流形法的速度進行修正。

4　高階流形法速度公式的修正

高階流形法速度公式的修正考慮以下3種方法。

4．1方法1——結(jié)點覆蓋上的速度修正

比較式（8）與式（6）可見，欲使Vn=V-n，則

考慮到式（11）中的Δu隨網(wǎng)格內(nèi)的材料點位置而變化，作為一種近似，在上述公式中可取結(jié)點處的Δu1和Δu2，即結(jié)點上的覆蓋多項式系數(shù)修正為

推廣到二維問題，設(shè)結(jié)點i（i=1，2，3）的第n步速度覆蓋函數(shù)為

式中ain至fi

n等為系數(shù)，三角形網(wǎng)格內(nèi)的速度分布為

其中Lin為在第n步構(gòu)形下的三角形網(wǎng)格形函數(shù)。

與一維情況的式（10）類似，相對于V-n的速度差值為

設(shè)三角形網(wǎng)格第i結(jié)點的位移增量為Δui和Δvi，相應(yīng)于一維公式（12），二維第n步的速度公式為

式中a-i

n至-fin為構(gòu)形變化前的多項式系數(shù)。

重新計算如圖1所示的剛性桿件旋轉(zhuǎn)問題，計算結(jié)果見表2，可見相對表1有很大改善，表明了覆蓋修正公式的作用。但與理論值相比還有一定的誤差，如t=1 s時，x坐標相差約7 cm，y坐標相差約5 cm。顯而易見，這是由速度覆蓋公式中近似取結(jié)點處的Δu和Δv所造成的，而準確的Δu和Δv在網(wǎng)格內(nèi)應(yīng)隨坐標變化。

4．2方法2——速度荷載修正

仍考慮如圖2所示的一維網(wǎng)格。網(wǎng)格內(nèi)每一點的位移變化為

表2　結(jié)點覆蓋速度修正的桿件自由端中點坐標Table 2　Coordinates of the middle point at the free end of bar（velocity modification of nodal covers）

其中，系數(shù) Δain和 Δbi

n是將它們均與Δ u1和Δ u2有關(guān)，這表明速度修正量+已不局限于結(jié)點覆蓋本身，而與一維網(wǎng)格內(nèi)的另一結(jié)點覆蓋有關(guān)。對于處在相鄰網(wǎng)格的交界點，在不同網(wǎng)格內(nèi)的速度修正量是不同的，這將造成修正后的速度不連續(xù)。

同理，對于二維三角形網(wǎng)格，將

及形函數(shù)L1

n，L2n，L3

n的表達式代入式（15）可得

同樣，上式中的Δain，Δbin，Δci

n與網(wǎng)格各結(jié)點的位移增量有關(guān)，對于處在相鄰網(wǎng)格的交界處，在不同網(wǎng)格內(nèi)的速度修正量也是不同的，因此，無法像方法1那樣對某一結(jié)點的速度覆蓋函數(shù)進行簡單修正，需要尋找其他方法。

考察流形法的基本控制方程，由最小勢能原理推導(dǎo)得出［1］，其中考慮慣性力的勢能表達式（仍以一維問題為例）為

式中：a=2u-2V，為加速度；ρ為密度；A為流形Δt2Δt元面積?？疾炱渲械牡?項∫ρu2VdA

AΔt，一般的做法是［1］，將網(wǎng)格內(nèi)的位移u=uTeTTe及速度V=TeVe代入，取極值后得到關(guān)于單元速度的荷載表達式為式中：Ve為網(wǎng)格結(jié)點速度的廣義自由度向量；M= ∫Aρ TTeTedA為單元質(zhì)量矩陣。

實際上，式（22）為單元的積分表達式，只要求計算Fv的積分，并未要求將V寫成包含結(jié)點速度的廣義自由度向量Ve的形式?？梢詢H將u=uTeTTe代入并取極值得這就給網(wǎng)格內(nèi)的速度修正（而不是結(jié)點速度修正）提供了依據(jù)。

修正后的速度荷載為

由于在計算速度荷載Fv時進行了速度修正，因而得到的增量位移Δ是比較準確的。但式（4）中的Vn仍采用了構(gòu)形變化前的系數(shù)，并沒有進行修正，因此計算出來的Vn+1仍然不準確，進而在計算下一步（第n+1步）的Fv時產(chǎn)生誤差而生成不平衡荷載。

下面考慮一種方法，在本步（第n步）式（24）的Fv中對速度增加一個修正量αΔV（α為系數(shù)），成為Vn-ΔV+αΔV，使得計算出的增量位移Δ un包含對Vn的修正，從而預(yù)先消除上述不平衡荷載。由于流形法計算要求每步的增量位移很小，這樣，相鄰計算步的速度變化量也很小，因此，基于Vn+1≈Vn的假設(shè)，式（4）可看成即2 V對應(yīng)于，可見，由修正量αΔV引起的

nFv荷載，其計算出的位移值Δ un，是對Vn的修正量的2倍。因此，系數(shù)α應(yīng)取為0．5。

再計算剛性桿件旋轉(zhuǎn)，計算結(jié)果見表3，可見相對表2又有很大改善，與理論值更接近。

表3　速度荷載修正的桿件自由端中點坐標Table 3　Coordinates of the middle point at the free end ofbar（modification of velocity loads）

4．3方法3——網(wǎng)格速度修正

對上述思路進一步擴展，認為同一結(jié)點在不同的網(wǎng)格中的速度不同，即速度不再連續(xù)。計算程序不再以結(jié)點的形式記錄速度，而以網(wǎng)格內(nèi)的速度信息來記錄，即使是連接相鄰網(wǎng)格的同一結(jié)點，它在各網(wǎng)格內(nèi)記錄的速度信息也是不同的。

對于二維三角形網(wǎng)格情況，按式（20）計算網(wǎng)格內(nèi)的結(jié)點速度修正量，對網(wǎng)格內(nèi)的各結(jié)點Vn直接進行修正（對于同一結(jié)點，在不同網(wǎng)格內(nèi)的修正值不一樣），然后代入式（23）計算速度荷載。該步計算完畢后，基于網(wǎng)格內(nèi)的速度值按式（4）計算新的速度。

剛性桿件旋轉(zhuǎn)的計算結(jié)果見表4，在保留小數(shù)點后三位有效數(shù)字時，與0階的計算結(jié)果完全相同，說明已經(jīng)完全排除了線性項誤差的影響。

表4　網(wǎng)格速度修正的桿件自由端中點坐標Table 4　Coordinates of the middle point at the free end of bar（modification of mesh velocities）

5結(jié)語

對于采用一階多項式覆蓋函數(shù)的高階流形法，本文提出了3種速度修正方法：

方法1（結(jié)點覆蓋修正）最簡單，但計算精度較差；方法2（速度荷載修正）的計算精度基本令人滿意，并能保持速度場的連續(xù)性，但要求相鄰計算步的速度差別小；方法3（網(wǎng)格速度修正）計算精度最高，但速度在相鄰網(wǎng)格間有間斷，不再保持連續(xù)性。3種方法各有優(yōu)缺點，在實際計算中可視情況選用。二階或更高階多項式覆蓋函數(shù)情況，可采用類似方法。另外，本文的研究是針對采用整體坐標來表示多項式覆蓋函數(shù)時所產(chǎn)生的高階速度問題，若采用局部化的覆蓋函數(shù)形式［10］，也許會減弱甚至避免上述問題，這還需進一步研究。

在結(jié)點處增加廣義自由度的各種方法中，數(shù)值流形方法已在大位移和大變形的計算中領(lǐng)先一步，本文方法對于其他具有廣義自由度的方法（如廣義有限元）也具有一定的參考價值。

致謝：感謝數(shù)值流形方法的發(fā)明者石根華博士的指導(dǎo)。感謝國家自然科學(xué)基金（10772034）的資助。

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

［2］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．

［3］DAUX C，MOES N，DOLBOW J，et al．Arbitrary Branched and Intersecting Cracks with the Extended Finite Element Method［J］．International Journal for Numerical Methods in Engineering，2000，48：1741-1760．

［4］ 張湘?zhèn)?，章爭榮，呂文閣．數(shù)值流形方法研究及應(yīng)用進展［J］．力學(xué)進展，2010，40（1）：1-12．

［5］李錄賢，王鐵軍．擴展有限元法（XFEM）及其應(yīng)用［J］．力學(xué)進展，2005，35（1）：5-20．

［6］李錄賢，劉書靜，張慧華．廣義有限元方法研究進展［J］．應(yīng)用力學(xué)學(xué)報，2009，26（1）：96-108．

［7］蘇海東，崔建華，謝小玲．高階數(shù)值流形方法的初應(yīng)力公式［J］．計算力學(xué)學(xué)報，2010，27（2）：270-274．

［8］ WANG C Y，CHUANG C C，SHENG J．Time Integration Theories for the DDA Method with Finite Element Meshes ［C］∥SALAMI M R，BANKS D．Proceedings of the First International Conference on Analysis of Discontinuous Deformation．Analysis（DDA）and Simulations of Discontinuous Media．Berkeley，CA：TSI Press，Albuquerque，NM，1996：263-287．

［9］蘇海東，謝小玲，陳琴．高階數(shù)值流形方法在結(jié)構(gòu)靜力分析中的應(yīng)用研究［J］．長江科學(xué)院院報，2005，22（5）：74-77．

［10］林紹忠，祁勇峰，蘇海東．數(shù)值流形法中覆蓋函數(shù)的改進形式及其應(yīng)用［J］．長江科學(xué)院院報，2006，23（6）：55-58．

（編輯：王慰）

Velocity Equations for High-order Numerical Manifold Method

SU Hai-dong
（Department of Material and Structure，Yangtze River Scientific Research Institute，Wuhan430010，China）

Computational accuracy of structure deformation can be improved greatly by the high-order Numerical Manifold Method（NMM）．However，poor accuracy or even incorrect velocity results were obtained in the dynamic analysis involved in large displacement．Based on 2-D triangular mathematical meshes and 1-order polynomial cover functions，the reason of the above cases is discussed through an example of rotation of a rigid bar in this paper．Three treatments and the corresponding equations for high-order velocities are presented for the first time，reflecting the change of configuration coordinates under large displacement．The high-order numerical manifold method is useful to other methods such as Generalized Finite Element Method（GFEM）which introduces generalized freedoms at nodes when solving geometric nonlinear problems．

Numerical Manifold Method（NMM）；high-order polynomial cover function；large displacement；velocity equation；generalized degree of freedom

O33；O34

A

1001-5485（2016）07-0121-05

2015-04-22；

2015-06-04

國家自然科學(xué)基金項目（10772034）

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