祁勇峰，蘇海東，崔建華

(長江科學院a.水利部水工程安全與病害防治工程技術研究中心;b.材料與結構研究所，武漢 430010)

1 研究背景

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

在結構應力分析方面，流形法相對于有限元法的優(yōu)勢主要體現(xiàn)在兩方面:一是數(shù)學網(wǎng)格與物理網(wǎng)格相互獨立，兩者邊界可以不重合，只要求前者覆蓋后者，因此可采用規(guī)則的數(shù)學網(wǎng)格對復雜的物理區(qū)域進行切分形成流形元，切分過程僅涉及簡單的幾何運算，前處理方便;二是可以通過提高覆蓋函數(shù)階數(shù)或構造特殊覆蓋函數(shù)來提高數(shù)值解的精度。

當前的流形法研究多采用有限單元網(wǎng)格來構造數(shù)學網(wǎng)格。對任一結點，與該結點相連的所有有限元形成一個數(shù)學覆蓋，這樣，任一原始的有限元網(wǎng)格就是其結點覆蓋的重疊部分，或者說，定義于有限元結點上的覆蓋在有限單元的區(qū)域內完全重疊，本文稱之為完全重疊覆蓋的流形法。這種方法能夠充分利用有限元法的基本公式，在流形法的初步應用研究中發(fā)揮了重要的作用［2－5］。但研究表明，在人們實際比較關注的一些結構關鍵部位，如孔周邊、裂紋尖端附近等，2套網(wǎng)格的不匹配可能造成計算精度的下降，流形法的計算結果往往達不到常規(guī)有限元法的精度水平。解決上述問題的一種途徑是對數(shù)學網(wǎng)格進行人工干預，如文獻［6］將數(shù)學網(wǎng)格結點布置在上述關鍵部位，雖能起到一定效果，但加大了人工工作量及處理難度，抵消了流形法在前處理上的一些優(yōu)勢。

針對上述問題，石根華博士指出，現(xiàn)在的流形法只是一個初級版本。實際上，流形覆蓋有著豐富的內涵，是非常靈活的。數(shù)值流形方法的基本思想來自于數(shù)學上的流形，即部分重疊的非空開集之并。在數(shù)學上，對于各種復雜的幾何形狀，用統(tǒng)一的坐標去描述并不方便。因此，近代數(shù)學的一個重要分支——微分幾何提出了流形思想，用各種獨立的覆蓋反映各自區(qū)域的幾何特征。覆蓋與覆蓋之間的重疊區(qū)域只是為了保持覆蓋之間的聯(lián)系，或曰協(xié)調性(為保證積分的唯一性，需要引入各覆蓋的權函數(shù)，并使權函數(shù)之和為1，這就是單位分解思想)。因此，原始的流形思想本身就是以獨立的覆蓋(即非空開集)為主，以其重疊區(qū)域為輔。

為此，石根華博士提出一種新的研究思路——部分重疊覆蓋的流形法。該方法采用以獨立覆蓋為主的分析方式，獨立覆蓋之間僅用較小的重疊區(qū)域保持連續(xù)性。初步分析，部分重疊覆蓋的流形法至少存在以下優(yōu)點:便于在局部區(qū)域采用特殊的覆蓋函數(shù)以適應物理場的局部特征，比如在裂紋尖端附近采用特殊的解析解級數(shù)，收斂性要比通常的多項式逼近快得多，這樣就可以將不同區(qū)域的多種求解形式協(xié)調地聯(lián)系起來進行聯(lián)合求解，從而實現(xiàn)局部解析法與其周邊數(shù)值解的完美結合;對獨立的覆蓋函數(shù)進行收斂性分析，不涉及其它覆蓋，比完全重疊下的覆蓋精度分析要方便;相鄰覆蓋之間的聯(lián)系僅限于其較小的重疊部分，在保持協(xié)調性的同時又降低了自由度之間的相互影響，能有效地改善方程組的性態(tài)，有利于方程組的快速迭代求解。

作為部分重疊覆蓋流形法的首次嘗試，本文基于矩形覆蓋，給出部分重疊的覆蓋形式、部分重疊覆蓋流形法的基本公式及其實現(xiàn)方式，對該方法在結構分析中的有效性進行初步驗證，為下一步的研究打下良好的基礎。

2 完全重疊的覆蓋形式

結合一個簡單的例子來說明覆蓋、網(wǎng)格、重疊覆蓋的概念。如圖1(a)所示的橢圓形的物理區(qū)域(物理網(wǎng)格)以及矩形有限元數(shù)學網(wǎng)格(編號為a—p)，經(jīng)過切割操作后形成物理覆蓋系統(tǒng)(著色區(qū)域)。以網(wǎng)格e為例，在其4個結點中，與結點1相連的所有有限元數(shù)學網(wǎng)格(a，e)構成結點1的數(shù)學覆蓋，與結點2相連的有限元網(wǎng)格(a，b，e，f)構成結點2的數(shù)學覆蓋，結點3和結點4的數(shù)學覆蓋類似，分別見圖1(b)—圖1(e)，圖中所示的矩形網(wǎng)格中的著色區(qū)域即為物理覆蓋。這4個物理覆蓋的交集為圖1(f)的數(shù)學網(wǎng)格e的著色區(qū)域，定義為流形元，圖中可見流形元可以僅占據(jù)網(wǎng)格內的部分區(qū)域，并具有任意的形狀。

圖1 完全重疊的數(shù)學覆蓋Fig.1 Completely overlapping math cover

從這個例子中可見，每個有限單元網(wǎng)格正好就是其幾個結點上的數(shù)學覆蓋的共同區(qū)域，或者說，這幾個覆蓋在有限單元的區(qū)域內完全重疊，這就是通常所見的完全重疊的覆蓋形式。

3 部分重疊的覆蓋形式

仍考慮圖1中的橢圓形物理區(qū)域，圖2為矩形網(wǎng)格的部分重疊覆蓋系統(tǒng)示意圖。其中，圖2(a)中a—p所示的矩形網(wǎng)格為數(shù)學覆蓋，覆蓋之間為較小的重疊區(qū)域，圖2(b)的著色區(qū)域即為物理網(wǎng)格與數(shù)學網(wǎng)格(覆蓋)相互切割后形成的物理覆蓋。

圖2 部分重疊的數(shù)學覆蓋Fig.2 Partially overlapping math cover

以圖2(c)中的局部 4 個數(shù)學覆蓋 a，b，e，f為例，每個覆蓋包含一個獨立覆蓋區(qū)域(圖2(c)中較大的矩形)和3個覆蓋重疊區(qū)域，其中，P1—P4為相鄰兩個覆蓋之間的重疊區(qū)域，P5為4個覆蓋共同的重疊區(qū)域。在完全重疊覆蓋的流形法中，流形單元定義為覆蓋重疊的公共部分，而對于部分重疊覆蓋的流形法，流形元還應包括獨立覆蓋中的物理區(qū)域(物理覆蓋)。

4 部分重疊覆蓋的實現(xiàn)技術

與通常所見的完全重疊的覆蓋體系不同，部分重疊的覆蓋體系涉及到眾多的獨立覆蓋及其周邊較小的覆蓋重疊區(qū)域，本來需要編制新的程序來實現(xiàn)。但實際上，在完全重疊的矩形覆蓋體系下，僅需采用簡單的結點之間強制約束的技術，就可以很方便地實現(xiàn)部分重疊的覆蓋形式。以下結合圖2(c)中所示的局部區(qū)域說明實現(xiàn)過程，見圖3(a)。

圖3 采用自由度約束實現(xiàn)部分重疊覆蓋Fig.3 Partial overlapping cover obtained by using freedom constraint

首先考察圖3(b)的數(shù)學網(wǎng)格1－2－6－5，網(wǎng)格內的流形元位移函數(shù)為

式中:Ui為定義在結點i上的多項式覆蓋函數(shù);Wi為權函數(shù)，

式中:ε0=εiε，η0=ηiη，(ε，η)為網(wǎng)格內任意點的局部坐標;(εi，ηi)為結點的局部坐標，見圖4。權函數(shù)滿足

將結點2，5，6“強制約束”到結點1，即令U2=U3=U4=U1，考慮到式(3)，則式(1)成為

圖4 矩形數(shù)學網(wǎng)格Fig.4 Rectangular math mesh

即在數(shù)學網(wǎng)格中保持權函數(shù)為1，使原始的完全重疊覆蓋1－2－6－5成為獨立覆蓋1。其他的獨立覆蓋類似，如圖3(a)所示，結點4，7，8約束到結點3，結點10，13，14 約束到結點9，結點12，15，16約束到結點11，分別形成獨立覆蓋3，9，11。

再考察上述4個覆蓋之間的重疊區(qū)域，以圖3(c)所示的覆蓋e和f之間的交叉區(qū)域P4(結點2－3－7－6)為例。結點2和6為同一覆蓋，局部覆蓋函數(shù)為U1，結點3和7也為同一覆蓋，局部覆蓋函數(shù)為U2，則網(wǎng)格內的位移函數(shù)為

式中W1+W4=(1－ξ)/2，W2+W3=(1+ξ)/2，為一維的權函數(shù)(ξ為局部坐標，－1≤ξ≤+1)，因此P4區(qū)域正好是2個一維覆蓋的重疊區(qū)域，在獨立覆蓋1的邊界1—1上滿足ξ=－1，獨立覆蓋3的邊界3—3上滿足ξ=1。P1—P3的重疊區(qū)域的情況類似。

對于4個覆蓋的重疊區(qū)域P5，位移函數(shù)與完全重疊覆蓋的情況一致。

在完全重疊覆蓋流形法程序的基礎上，僅需加入結點自由度之間的關聯(lián)程序段，就可以實現(xiàn)部分重疊的覆蓋，實現(xiàn)過程簡單方便。在覆蓋之間的重疊區(qū)域，比如圖3(1)所示的P1—P5，數(shù)學網(wǎng)格的各結點對應于不同的獨立覆蓋，需要推導相應的流形元計算公式。

5 部分重疊覆蓋的流形單元分析

在完全重疊覆蓋的流形法中，基于多項式覆蓋函數(shù)和矩形數(shù)學網(wǎng)格，采用矩陣特殊運算推導了流形元計算公式，但各結點采用的是相同的覆蓋函數(shù)［7］。部分重疊覆蓋流形法的計算公式與前者大體相同，但可靈活實現(xiàn)不同結點采用不同的覆蓋函數(shù)，以下相應的位移函數(shù)與剛度矩陣公式。

5.1 位移函數(shù)

矩形網(wǎng)格中的位移函數(shù)可表示為

式中:權函數(shù)

其中:

式中:定義在結點i上的多項式覆蓋函數(shù)項數(shù)為m(i);uij，vij為覆蓋函數(shù)的系數(shù);Fi是權函數(shù)的系數(shù)列陣;ξi，ηi(i=1～4)為4 個結點的局部坐標，依次為(－1，－1)，(1，－1)，(1，1)，(－1，1);xc，yc為矩形網(wǎng)格的形心坐標;a，b為矩形網(wǎng)格的邊長;twξ是多項式基底，I為2階單位矩陣。

5.2 剛度矩陣

流形元應變矩陣為

式中，

剛度矩陣為

其中，

流形元剛度矩陣的子塊(2×2)也可表示為

式中:D為彈性矩陣;B為幾何矩陣;m(i)，m(k)分別表示結點覆蓋i，k的多項式項數(shù)。

6 算例分析

如圖5所示的重力壩，壩高100m，壩底長為60m，壩頂寬度20m，庫水深80m。壩體彈性模量為30 GPa，泊松比為0.167。計算上游面承受庫水壓力情況下的應力和變形。采用稠密網(wǎng)格的有限元計算結果作為對比，劃分1 080個4結點等參元，結點總數(shù)為1 183個。

圖5 計算網(wǎng)格Fig.5 Computational mesh

采用疏密程度不同的流形元網(wǎng)格進行計算，其中:疏網(wǎng)格有8個獨立覆蓋，包括其重疊區(qū)域在內共有流形單元21個;密網(wǎng)格有36個獨立覆蓋，流形單元總數(shù)為112個。計算模型的底部全約束。計算結果及與有限元的對比見表1和圖6、圖7，對比分析上游表面的x向位移及x，y向的應力，其中水位以下的x向應力應等于該點處的靜水壓力。

表1 上游面特征點x向位移對比Table 1 Comparison of horizontal displacements on the upstream surface calculated by manifold method and finite element method

表1的上游面x向位移顯示，疏網(wǎng)格的3階流形法計算結果與有限元總體上接近，但在靠近壩踵的部位(y=10m)由于約束的影響計算結果稍差。僅對上游底部的獨立覆蓋升階，保持其他覆蓋的階數(shù)不變，至6，7階時，位移與有限元解一致。密網(wǎng)格3階流形法的結果與有限元解一致。

圖6 上游面特征點x向應力對比Fig.6 Comparison of horizontal stresses on the upstream surface calculated by manifold method and finite element method

圖7 上游面特征點y向應力對比Fig.7 Comparison of vertical stresses on the upstream surface calculated by manifold method and finite element method

圖6和圖7的上游面x向及y向應力對比表明，疏網(wǎng)格采用3階流形法時，在中上部(y=40m以上)與有限元的結果接近，而在底部覆蓋處由于約束的影響難以獲得較好的結果。因此，將上游底部的獨立覆蓋升至7階后，其他覆蓋的階數(shù)保持不變，x向、y向應力比3階情況有較大改善，除了在靠近壩踵的部位(y=10m)處計算結果稍差外，其它部位與有限元解基本一致。密網(wǎng)格的2階流形法結果稍差，將上游的10個獨立覆蓋升至3階，其他覆蓋保持2階不變，計算結果與有限元解很接近。

以上算例表明，部分重疊覆蓋的流形法可以得到較好的計算結果。對于某些計算結果不理想的獨立覆蓋區(qū)域，通過提高該覆蓋中的多項式階數(shù)也可以提高計算精度，且對周邊覆蓋中的計算結果影響較小。但對于應力和變形過于復雜的獨立覆蓋區(qū)域往往需要很高的覆蓋函數(shù)階次，因此有時單純提高階次的效果不佳，而相比之下通過加密覆蓋，即使采用低階覆蓋函數(shù)也可以得到很好的計算結果。

7 結語

本文將現(xiàn)有的基于完全重疊覆蓋的流形法擴展到一般意義上的流形研究——部分重疊覆蓋的流形法，通過單個矩形數(shù)學網(wǎng)格各結點之間的強制約束方式，很方便地將完全重疊的覆蓋形式及流形法公式轉化為部分重疊的覆蓋形式及公式，算例分析初步驗證了這種新型流形法的有效性。

部分重疊覆蓋的流形法是一種以獨立覆蓋分析為主的全新計算模式，本文的工作僅僅是其研究的開始。我們在進一步的研究中開展了以下工作:在獨立覆蓋中采用特殊的覆蓋函數(shù)以適應物理場的局部特征，如應用裂紋尖端解析覆蓋計算裂紋尖端的位移場及應力強度因子［8］;根據(jù)需要進行大、小覆蓋之間的疏密布置并保持協(xié)調過渡，使其具有更大的靈活性;采用部分離散的思想初步形成基于流形思想的有限條分析方法，將以往只能用于規(guī)則形狀的有限條法擴展到分析復雜形狀。這些工作將另文介紹。

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

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

［2］彭自強，葛修潤.數(shù)值流形方法在八結點有限元網(wǎng)格上的實現(xiàn)［J］.武漢大學學報(工學版)，2004，37(2):72－76.(PENG Zi-qiang，GE Xiu-run.Numerical Manifold Method on Eight-node Isoparametric Element Meshes［J］.Engineering Journal of Wuhan University，2004，37(2):72－76.(in Chinese))

［3］林紹忠，明崢嶸，祁勇峰.用數(shù)值流形法分析溫度場及溫度應力［J］.長江科學院院報，2007，24(5):72－75.(LIN Shao-zhong，MING Zheng-rong，QI Yong-feng.Thermal Field and Thermal Stress Analysis Based on Numerical Manifold Method［J］.Journal of Yangtze River Scientific Research Institute，2007，24(5):72－75.(in Chinese))

［4］蘇海東，黃玉盈.數(shù)值流形方法在流固耦合諧振分析中的應用［J］.計算力學學報，2007，24(6):823－828.(SUHai-dong，HUANGYu-ying.Application of Numerical Manifold Method in Fluid-solid Interaction Harmonic Analysis［J］. Journal of Computational Mechanics，2007，24(6):823－828.(in Chinese))

［5］JIANG Qing-hui，ZHOU Chuang-bing，LI Dian-qing.A Three-Dimensional Numerical Manifold Method Based on Tetrahedral Meshes［J］.Computers ＆ Structures，2009，87(7):880－889.

［6］蘇海東，謝小玲，陳 琴.高階數(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))

［7］林紹忠，祁勇峰，蘇海東.基于矩陣特殊運算的高階流形元單元分析［J］.長江科學院院報，2006，23(6):36－39.(LIN Shao-zhong，QI Yong-feng，SU Hai-dong.Element Analysis of High-Order Numerical Manifold Method Based on Special Matrix Operations［J］.Journal of Yangtze River Scientific Research Institute，2006，23(6):36－39.(in Chinese))

［8］SU Hai-dong，QI Yong-feng.A Novel Numerical Method for Computing Stress Intensity Factors［J］.Applied Mechanics and Materials，2012，204－208:4486－4489.