蘇海東，祁勇峰

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

1 研究背景和思路

數(shù)值流形方法［1］(以下簡稱流形法)是留美學(xué)者石根華博士提出的一種新型數(shù)值計(jì)算方法。該方法采用了覆蓋的概念和2套相互獨(dú)立的網(wǎng)格體系。覆蓋，或稱為數(shù)學(xué)覆蓋，用以定義各個區(qū)域的局部函數(shù)。兩套相互獨(dú)立的網(wǎng)格體系，即反映數(shù)值解精度的數(shù)學(xué)網(wǎng)格和表示幾何邊界和材料分區(qū)的物理網(wǎng)格，將整個研究區(qū)域劃分成有限個相互重疊的區(qū)域。這些區(qū)域被稱為物理覆蓋，在各個覆蓋上獨(dú)立定義局部覆蓋函數(shù)，通過聯(lián)系函數(shù)(或稱為權(quán)函數(shù))加權(quán)平均得到整個求解域上的總體函數(shù)。2套網(wǎng)格的交集，稱為流形元，跟有限元一樣，是基本計(jì)算單元，但可具有復(fù)雜任意的形狀。為保證任意形狀的流形元的積分精度，流形法多采用單純形精確積分方式［1－2］。

當(dāng)前的流形法研究一般采用有限單元網(wǎng)格來構(gòu)造數(shù)學(xué)網(wǎng)格［1，3］，權(quán)函數(shù)取為有限單元的形函數(shù)。對任一結(jié)點(diǎn)，與該結(jié)點(diǎn)相連的所有有限單元形成一個數(shù)學(xué)覆蓋，這樣，任一原始的有限元網(wǎng)格就是其結(jié)點(diǎn)覆蓋的重疊部分，或者說，定義于有限元結(jié)點(diǎn)上的覆蓋在有限單元的區(qū)域內(nèi)完全重疊，作者在文獻(xiàn)［4］中將這種類型的流形法稱為完全重疊覆蓋的流形法。研究表明，在人們實(shí)際比較關(guān)注的一些結(jié)構(gòu)關(guān)鍵部位，如孔的周邊、裂紋尖端附近等，2套網(wǎng)格的不匹配經(jīng)常會造成計(jì)算精度的下降，流形法的計(jì)算結(jié)果往往達(dá)不到常規(guī)有限元法的精度水平［5］。

作者在文獻(xiàn)［4，6］中提出部分重疊覆蓋的流形法，采用以獨(dú)立覆蓋為主的分析方式，便于在各區(qū)域用合適的覆蓋函數(shù)以適應(yīng)物理場的局部特征，獨(dú)立覆蓋之間僅用較小的重疊區(qū)域以保持覆蓋的連續(xù)性，從而將現(xiàn)有的基于完全重疊覆蓋的流形法擴(kuò)展到一般意義上的流形研究。作為部分重疊覆蓋流形法的首次嘗試，該文基于矩形覆蓋，給出部分重疊的覆蓋形式、基本公式及其實(shí)現(xiàn)方式，通過單個矩形數(shù)學(xué)網(wǎng)格各結(jié)點(diǎn)之間的強(qiáng)制約束方式，很方便地將完全重疊的覆蓋形式及公式轉(zhuǎn)化為部分重疊的覆蓋形式及公式。

眾所周知，孔周邊等局部區(qū)域具有不同于一般情況的復(fù)雜變形和應(yīng)力分布。如果采用較大的數(shù)學(xué)覆蓋、用多項(xiàng)式函數(shù)去逼近這樣復(fù)雜的分布，理論上是可行的，但往往需要很高階數(shù)的多項(xiàng)式級數(shù)，而且高階情況下由于方程組的性態(tài)不佳，易引入數(shù)值誤差，因此效果很差。實(shí)踐表明，在大的覆蓋中單純依靠提高覆蓋函數(shù)階次的方法往往會帶來計(jì)算結(jié)果的振蕩跳躍。反之，如果采用較小的覆蓋而用相對簡單的低階多項(xiàng)式，則可以更好地逼近實(shí)際復(fù)雜的分布情況［4］。

這些局部較小的覆蓋與其周邊較大覆蓋的連接是需要考慮的問題。文獻(xiàn)［4］采用了統(tǒng)一大小的規(guī)則矩形覆蓋，在這種情況下，為了將局部區(qū)域的覆蓋變小以提高計(jì)算精度，需要將所有的覆蓋都變小，或者至少要將其上、下、左、右4個方向的覆蓋都變小，從而增加了很多不必要的計(jì)算量。因此要完全發(fā)揮出部分重疊覆蓋流形法的優(yōu)勢，需要解決如何從小覆蓋過渡到周圍較大覆蓋的問題。本文繼續(xù)文獻(xiàn)［4］的工作，基于矩形覆蓋提出局部覆蓋加密技術(shù)，實(shí)現(xiàn)大、小覆蓋(或稱為網(wǎng)格)的協(xié)調(diào)過渡。

2 有限元網(wǎng)格的局部網(wǎng)格加密

對于完全重疊覆蓋流形法而言，實(shí)現(xiàn)大、小覆蓋(網(wǎng)格)的過渡相對不太方便。

如圖1所示的平面計(jì)算，通常在大、小網(wǎng)格之間采用三角形或任意形狀的四邊形進(jìn)行過渡，這些方式一般需要人工操作，不能完全體現(xiàn)流形法在網(wǎng)格自動布置上的優(yōu)勢，而且過渡的三角形或四邊形往往形狀不佳，對計(jì)算精度有影響。如文獻(xiàn)［6］計(jì)算裂紋尖端的應(yīng)力強(qiáng)度因子，采用基于面積坐標(biāo)的任意四邊形數(shù)學(xué)網(wǎng)格進(jìn)行大小網(wǎng)格過渡(使用面積坐標(biāo)是因?yàn)閱渭冃畏e分公式要求被積函數(shù)為多項(xiàng)式)。結(jié)果表明，在裂紋周邊覆蓋取為高階多項(xiàng)式時(shí)，計(jì)算值普遍比理論值偏大。

圖1 大、小網(wǎng)格的協(xié)調(diào)過渡(有限元網(wǎng)格)Fig.1 Transition between small meshes and big meshes when using FEM meshes(compatible approach)

另一種是不協(xié)調(diào)的過渡方式，如圖2所示，大、小網(wǎng)格相連邊不匹配，需要采用粘接方法，強(qiáng)制密網(wǎng)格邊上的自由度符合粗網(wǎng)格邊上的位移分布。以上過渡方法均是常規(guī)有限元網(wǎng)格的做法。對于裂紋尖端等區(qū)域，位移場分布較為復(fù)雜，需要劃分非常細(xì)密的網(wǎng)格(比常規(guī)網(wǎng)格的尺寸小1至2個數(shù)量級)，用上述方式需要進(jìn)行多層的過渡，操作起來很不方便。

圖2 大、小網(wǎng)格的不協(xié)調(diào)過渡(有限元網(wǎng)格)Fig.2 Transition between small meshes and big meshes when using FEM meshes(incompatible approach)

3 部分重疊覆蓋的覆蓋加密

針對部分重疊的矩形覆蓋，本文提出大小覆蓋的協(xié)調(diào)過渡方法，以下以圖3所示的例子進(jìn)行說明。在圖3(a)的獨(dú)立覆蓋5中進(jìn)行局部覆蓋加密，形成圖3(b)的形式，圖中顯示了與原覆蓋5相關(guān)區(qū)域的結(jié)點(diǎn)編號。

圖3 局部增加完全重疊覆蓋Fig.3 Cover refinement via adding totally overlapping covers

相關(guān)結(jié)點(diǎn)的控制原則是:

(1)與獨(dú)立覆蓋相連的結(jié)點(diǎn)由獨(dú)立覆蓋控制(與獨(dú)立覆蓋的自由度相同)，比如，結(jié)點(diǎn) 1，6，31，36仍分別由獨(dú)立覆蓋 1，3，7，9 控制;結(jié)點(diǎn) 2，3，4，5 由獨(dú)立覆蓋2控制，結(jié)點(diǎn)7，13，19，25由獨(dú)立覆蓋4控制，結(jié)點(diǎn)12，18，24，30 由獨(dú)立覆蓋6 控制，結(jié)點(diǎn)32，33，34，35由獨(dú)立覆蓋8控制。

(2)其它內(nèi)部結(jié)點(diǎn)則為通常的完全重疊覆蓋，如結(jié)點(diǎn)8至結(jié)點(diǎn)11，14至17，20至23，26至29均為完全重疊覆蓋的結(jié)點(diǎn)。

這實(shí)際上是一種在獨(dú)立覆蓋中增加完全重疊覆蓋的加密方式。

以下討論這種過渡的協(xié)調(diào)問題。由于內(nèi)部完全重疊覆蓋的網(wǎng)格是協(xié)調(diào)的，因此重點(diǎn)關(guān)注與獨(dú)立覆蓋相連的局部網(wǎng)格。以圖3(b)中與獨(dú)立覆蓋6相連的網(wǎng)格23－17－18－24為例，這種矩形網(wǎng)格內(nèi)部結(jié)點(diǎn)編號見圖4(a)，其權(quán)函數(shù)(即有限單元網(wǎng)格的形函數(shù))為式(1)［7］，

式中:ξ0= ξiξ，η0= ηiη，(ξi，ηi)為第 i結(jié)點(diǎn)的局部坐標(biāo)，見圖4(a)。

如圖4(b)所示，保持結(jié)點(diǎn)覆蓋1和2上的權(quán)函數(shù)為式(1)不變，將內(nèi)部結(jié)點(diǎn)3和4所在邊合成為1個覆蓋，其權(quán)函數(shù)為

容易驗(yàn)證，W3滿足:在覆蓋3處為1;在其他2個覆蓋處為0;3個權(quán)函數(shù)之和恒為1。因此在此網(wǎng)格內(nèi)部滿足流形覆蓋聯(lián)系函數(shù)(權(quán)函數(shù))的所有條件，在水平方向上，由結(jié)點(diǎn)覆蓋1和2過渡到獨(dú)立覆蓋3是協(xié)調(diào)的。再考察與網(wǎng)格23－17－18－24上、下相連的網(wǎng)格，比如17－11－12－18，兩者共用一條邊17－18，在上、下方向過渡也是協(xié)調(diào)的。因而總體上這種過渡方式是完全協(xié)調(diào)的。

圖4 矩形網(wǎng)格Fig.4 A rectangular mesh

實(shí)際上，仍采用一般的權(quán)函數(shù)式(1)，用文獻(xiàn)［4］所述的強(qiáng)制約束方式，將結(jié)點(diǎn)4約束到結(jié)點(diǎn)3(即結(jié)點(diǎn)4和結(jié)點(diǎn)3的自由度相同)，就可以很方便地實(shí)現(xiàn)式(2)。

以上方式可以根據(jù)計(jì)算精度的要求在局部提高網(wǎng)格密度，如在圖3(b)的基礎(chǔ)上進(jìn)一步加密成圖5的形式。

圖5 進(jìn)一步加密覆蓋Fig.5 Further refinement

除了上述在局部增加完全重疊覆蓋的加密方式外，還可以增加部分重疊的覆蓋，如圖6所示，在圖3(a)的獨(dú)立覆蓋5中增加獨(dú)立覆蓋10至18，原理與上述完全重疊覆蓋的情況類似。

進(jìn)一步推廣到采用完全重疊覆蓋和部分重疊覆蓋相結(jié)合的多種形式，所有形式都可以保證協(xié)調(diào)性，因此這種覆蓋加密的方式是非常靈活的。直觀上看，這種覆蓋加密方式也是很方便的，只要在需要加密的覆蓋內(nèi)增加橫線和豎線。采用文獻(xiàn)［4］所述的強(qiáng)制約束方式，按照上文所述的控制原則，可以由計(jì)算機(jī)自動生成結(jié)點(diǎn)之間的約束關(guān)系。另外，本文雖然僅討論了二維問題，但其思路同樣適合于三維的長方體覆蓋情況。

4 算例

以下通過2個算例驗(yàn)證上述方法的有效性。

4.1 算例1——重力壩頂部受集中力作用

仍以文獻(xiàn)［4］中的重力壩為例，如圖7所示，壩高100 m，壩底長為60 m，壩頂寬度20 m，壩體彈性模量為30 GPa，泊松比為0.167?？紤]壩頂?shù)纳嫌吸c(diǎn)作用F=1 000 kN集中力，關(guān)注壩體上游面的應(yīng)力，特別考察部分重疊覆蓋流形法模擬自由表面零應(yīng)力(理論上，上游面σx和τxy為0)的能力。采用稠密網(wǎng)格的有限元計(jì)算結(jié)果作為對比，劃分1 080個4結(jié)點(diǎn)等參元，結(jié)點(diǎn)總數(shù)為1 183。

圖7 重力壩及其有限元網(wǎng)格Fig.7 A gravity dam and its FEM meshes

如圖8所示的流形元1和流形元2兩種網(wǎng)格。流形元1為未采用局部加密的網(wǎng)格，計(jì)算結(jié)果見表1(同時(shí)列出了采用稠密網(wǎng)格的有限元結(jié)果作為對比)，可見，在上游頂、底部2個獨(dú)立覆蓋中的上游面應(yīng)力精度稍差。因此，在流形元2的網(wǎng)格中，局部采用完全重疊覆蓋對頂、底部2個獨(dú)立覆蓋進(jìn)行加密，并取獨(dú)立覆蓋n=4階，局部加密部分取n=2階?？梢?，流形元2與流形元1相比，應(yīng)力精度有明顯改善:σx和τxy總體上更接近于0，大部分點(diǎn)的數(shù)值在0.1以下，甚至明顯好于稠密網(wǎng)格的有限元計(jì)算結(jié)果;σy與有限元的結(jié)果更接近，大部分點(diǎn)(y=10 m至y=70 m處)與有限元結(jié)果相差在0.5%以下，其中y=70 m處，流形元1計(jì)算得到的應(yīng)力數(shù)值為27.81，與有限元解30.76相差約 10%，而流形元 2的結(jié)果為30.75，與有限元解基本相同。

表1 上游表面應(yīng)力Table 1 Stresses on upstream face of the dam 0.01 MPa

圖8 不同流形元計(jì)算網(wǎng)格Fig.8 Two NMM meshes

相對而言，高程80，90 m處的精度比其它部位差，說明頂部網(wǎng)格還有進(jìn)一步加密的必要(單純提高多項(xiàng)式階數(shù)的效果反而不好)。

進(jìn)一步在頂部覆蓋中采用部分重疊的覆蓋加密方式，如圖9所示的流形元3的網(wǎng)格，所有的獨(dú)立覆蓋均采用n=4階，頂部的覆蓋由于有集中力的作用，取n=5階。計(jì)算結(jié)果見表1。與流形元2的結(jié)果相比，計(jì)算精度進(jìn)一步提高。在頂部覆蓋處的 σx和τxy更接近于0。

4.2 算例2——圓孔的應(yīng)力集中

中部開圓孔的矩形板在上、下兩端承受均布拉力p，通過計(jì)算圓孔左右兩側(cè)的垂直向應(yīng)力與施加的均布拉力之比來求應(yīng)力集中系數(shù)，理論值為3.00，即圓孔處的應(yīng)力為矩形板兩端拉力p的3倍。

考慮了3種網(wǎng)格如圖10所示，圓孔周邊用完全重疊覆蓋的方式加密。

圖9 流形元3的網(wǎng)格Fig.9 NMM meshes 3

圖10 圓孔周邊覆蓋加密Fig.10 Cover refinement around a circular hole

首先計(jì)算沒有進(jìn)行局部覆蓋加密的網(wǎng)格a，其中，圖中的圓孔底部到結(jié)構(gòu)底部的一條豎線是為了使結(jié)構(gòu)的矩形邊界與圓孔之間形成有向環(huán)路而設(shè)置的［1］。總共有25個獨(dú)立覆蓋，在圓孔所在的獨(dú)立覆蓋中，當(dāng)多項(xiàng)式階數(shù)n=2時(shí)，應(yīng)力集中系數(shù)為1.03;n=4時(shí)，應(yīng)力集中系數(shù)為1.09;n=6時(shí)，應(yīng)力集中系數(shù)為1.18?？梢?，應(yīng)力集中系數(shù)隨著多項(xiàng)式階數(shù)而增大，但收斂很慢，可以預(yù)見，要達(dá)到理論值3.00，需要很高的階數(shù)。

再計(jì)算網(wǎng)格b，在圓孔所在的獨(dú)立覆蓋中進(jìn)行覆蓋加密。周圍的獨(dú)立覆蓋取4階，計(jì)算得到的應(yīng)力集中系數(shù)見表2，表中n表示局部加密的完全重疊覆蓋的多項(xiàng)式階數(shù)?？梢?，計(jì)算值雖然接近于理論值，但隨階數(shù)n的升高會出現(xiàn)數(shù)值振蕩，這是在網(wǎng)格密度不夠情況下的典型表現(xiàn)。實(shí)踐表明，在實(shí)際物理場分布復(fù)雜的區(qū)域，單純提高多項(xiàng)式的階次往往會帶來計(jì)算結(jié)果的振蕩跳躍。

表2 圓孔應(yīng)力集中系數(shù)Table 2 Stress concentration factors of the circular hole

再計(jì)算網(wǎng)格c，在圓孔所在的獨(dú)立覆蓋中進(jìn)一步加密覆蓋。周圍的獨(dú)立覆蓋取2階和5階兩種情況，計(jì)算結(jié)果見表2?？梢?，計(jì)算結(jié)果隨著局部重疊覆蓋階數(shù)的提高而較為穩(wěn)定地趨向于理論值。

以上算例表明了本文提出的覆蓋加密方法的有效性。

5 結(jié)語

本文針對部分重疊覆蓋的流形法提出了大、小矩形覆蓋之間的過渡方式，可以在大覆蓋中根據(jù)需要加密成為較小的覆蓋，并與周邊的大覆蓋保持協(xié)調(diào)的聯(lián)系。此方法相對于以往采用有限元網(wǎng)格的完全重疊覆蓋流形法而言，在大小網(wǎng)格的過渡上要方便得多。當(dāng)然，如何布置加密的覆蓋(網(wǎng)格)以及如何安排覆蓋函數(shù)的階次使計(jì)算達(dá)到所需的精度，是需要進(jìn)一步研究的問題。

部分重疊覆蓋流形法的一大優(yōu)勢是便于在局部區(qū)域采用特殊的覆蓋函數(shù)以適應(yīng)物理場的局部特征。比如，對于裂紋尖端附近區(qū)域，其位移往往具有解析的函數(shù)形式，采用這種解析函數(shù)的特殊覆蓋，其收斂性肯定要比通常的多項(xiàng)式逼近快得多［8］。然而，這種特殊覆蓋函數(shù)的適用區(qū)域非常小。因此，在大的覆蓋中進(jìn)行覆蓋加密，不僅可以降低每個覆蓋中實(shí)際函數(shù)分布的復(fù)雜度，而且還可以提供適用于解析解的特殊覆蓋區(qū)域。采用本文提出的部分重疊覆蓋流形法的局部加密方法，可以很方便地實(shí)現(xiàn)周邊大覆蓋到裂紋尖端局部覆蓋的協(xié)調(diào)過渡，這部分研究成果將另文介紹。

致謝:部分重疊覆蓋的流形法思想是由石根華博士提出的，本文的研究工作得到了石根華博士的指導(dǎo)，在此表示衷心的感謝!

［1］石根華.數(shù)值流形方法與非連續(xù)變形分析［M］.裴覺民譯.北京:清華大學(xué)出版社，1997.(SHI Gen-hua.Numerical Manifold Method(NMM)and Discontinuous Deformation Analysis(DDA)［M］.Translated by PEI Jue-min.Beijing:Tsinghua University Press，1997.(in Chinese))

［2］林紹忠.單純形積分的遞推公式［J］.長江科學(xué)院院報(bào)，2005，22(3):32 － 34.(LIN Shao-zhong.Recursive Formula for Simplex Integration［J］.Journal of Yangtze River Scientific Research Institute，2005，22(3):32 －34.(in Chinese))

［3］張湘?zhèn)?，章爭榮，呂文閣，等.數(shù)值流形方法研究及應(yīng)用進(jìn)展［J］.力學(xué)進(jìn)展，2010，40(1):1－12.(ZHANG Xiang-wei，ZHANG Zheng-rong，LV Wen-ge，et al.Advances and Perspectives in Numerical Manifold Method and Its Applications［J］.Advances in Mechanics，2010，40(1):1 －12.(in Chinese))

［4］祁勇峰，蘇海東.部分重疊覆蓋的數(shù)值流形方法初步研究［J］.長江科學(xué)院院報(bào)，2013，30(1):65 －70.(QI Yong-feng，SU Hai-dong.Preliminary Study of Numerical Manifold Method with Partially Overlapping Covers［J］.Journal of Yangtze River Scientific Research Institute，2013，30(1):65 －70.(in Chinese))

［5］蘇海東，謝小玲，陳 琴.高階數(shù)值流形方法在結(jié)構(gòu)靜力分析中的應(yīng)用研究［J］.長江科學(xué)院院報(bào)，2005，22(5):74 － 77.(SU Hai-dong，XIE Xiao-ling，CHEN Qin.Application of High-order Numerical Manifold Method in Static Analysis［J］.Journal of Yangtze River Scientific Research Institute，2005，22(5):74 － 77.(in Chinese))

［6］祁勇峰，蘇海東.部分重疊覆蓋的流形法研究報(bào)告［R］.武漢:長江科學(xué)院，2012.(QI Yong-feng，SU Haidong.Research of Numerical Manifold Method with Partially Overlapping Covers［R］.Wuhan:Yangtze River Scientific Research Institute，2012.(in Chinese))

［7］王勖成，邵 敏.有限單元法的基本概念和數(shù)值方法［M］.北京:清華大學(xué)出版社，1988.(WANG Xu-cheng，SHAO Min.Basic Principle of the FEM and Numerical Method［M］.Beijing:Tsinghua University Press，1988.(in Chinese))

［8］蘇海東，祁勇峰，龔亞琦.裂紋尖端解析解與周邊數(shù)值解聯(lián)合求解應(yīng)力強(qiáng)度因子［J］.長江科學(xué)院院報(bào)，2013，30(6):83 － 89.(SU Hai-dong，QI Yong-feng，GONG Ya-qi.Compute Stress Intensity Factors via Combining Analytical Solutions around Crack Tips with Surrounding Numerical Solutions［J］.Journal of Yangtze River Scientific Research Institute，2013，30(6):83 －89.(in Chinese))