蘇海東

?

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

蘇海東

（長江科學院材料與結構研究所，武漢430010）

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

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

doi：10．11988/ckyyb．20150343

1　研究背景

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

式中N1，N2分別為對應于結點1和結點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構形更新至第n構形，附著在材料體上的數(shù)學網(wǎng)格結點隨材料體移動至新的結點位置x1n和x2n，網(wǎng)格內的任一點由xn-1移至xn。若按原始的高階流形法速度公式，結點1和結點2的第n步速度覆蓋函數(shù)分別為其中，在構形更新過程中的多項式系數(shù)保持不變，即?？紤]到剛體運動中的形函數(shù)N1和

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

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

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

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

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

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

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

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

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

式中ain至fi

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

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

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

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

式中a-i

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

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

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

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

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

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

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

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

及形函數(shù)L1

n，L2n，L3

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

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

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

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

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

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

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

修正后的速度荷載為

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

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

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

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

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

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

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

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

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

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

5結語

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

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

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

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

［1］ 石根華．數(shù)值流形方法與非連續(xù)變形分析［M］．裴覺民，譯．北京：清華大學出版社，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)ィ聽帢s，呂文閣．數(shù)值流形方法研究及應用進展［J］．力學進展，2010，40（1）：1-12．

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

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

［7］蘇海東，崔建華，謝小玲．高階數(shù)值流形方法的初應力公式［J］．計算力學學報，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ù)值流形方法在結構靜力分析中的應用研究［J］．長江科學院院報，2005，22（5）：74-77．

［10］林紹忠，祁勇峰，蘇海東．數(shù)值流形法中覆蓋函數(shù)的改進形式及其應用［J］．長江科學院院報，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

國家自然科學基金項目（10772034）

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