蘇海東，祁勇峰

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

1 研究背景和思路

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

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

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

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

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

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

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

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

圖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)格相連邊不匹配，需要采用粘接方法，強制密網(wǎng)格邊上的自由度符合粗網(wǎng)格邊上的位移分布。以上過渡方法均是常規(guī)有限元網(wǎng)格的做法。對于裂紋尖端等區(qū)域，位移場分布較為復雜，需要劃分非常細密的網(wǎng)格(比常規(guī)網(wǎng)格的尺寸小1至2個數(shù)量級)，用上述方式需要進行多層的過渡，操作起來很不方便。

圖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所示的例子進行說明。在圖3(a)的獨立覆蓋5中進行局部覆蓋加密，形成圖3(b)的形式，圖中顯示了與原覆蓋5相關區(qū)域的結點編號。

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

相關結點的控制原則是:

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

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

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

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

式中:ξ0= ξiξ，η0= ηiη，(ξi，ηi)為第 i結點的局部坐標，見圖4(a)。

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

容易驗證，W3滿足:在覆蓋3處為1;在其他2個覆蓋處為0;3個權函數(shù)之和恒為1。因此在此網(wǎng)格內(nèi)部滿足流形覆蓋聯(lián)系函數(shù)(權函數(shù))的所有條件，在水平方向上，由結點覆蓋1和2過渡到獨立覆蓋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ù)式(1)，用文獻［4］所述的強制約束方式，將結點4約束到結點3(即結點4和結點3的自由度相同)，就可以很方便地實現(xiàn)式(2)。

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

圖5 進一步加密覆蓋Fig.5 Further refinement

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

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

4 算例

以下通過2個算例驗證上述方法的有效性。

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

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

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

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

表1 上游表面應力Table 1 Stresses on upstream face of the dam 0.01 MPa

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

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

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

4.2 算例2——圓孔的應力集中

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

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

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

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

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

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

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

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

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

5 結語

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

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

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

［1］石根華.數(shù)值流形方法與非連續(xù)變形分析［M］.裴覺民譯.北京:清華大學出版社，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］.長江科學院院報，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)ィ聽帢s，呂文閣，等.數(shù)值流形方法研究及應用進展［J］.力學進展，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］.長江科學院院報，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ù)值流形方法在結構靜力分析中的應用研究［J］.長江科學院院報，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］祁勇峰，蘇海東.部分重疊覆蓋的流形法研究報告［R］.武漢:長江科學院，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］.北京:清華大學出版社，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)合求解應力強度因子［J］.長江科學院院報，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))