部分重疊覆蓋的數(shù)值流形方法初步研究

祁勇峰，蘇海東，崔建華

(長(zhǎng)江科學(xué)院a.水利部水工程安全與病害防治工程技術(shù)研究中心;b.材料與結(jié)構(gòu)研究所，武漢 430010)

1 研究背景

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

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

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

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

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

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

2 完全重疊的覆蓋形式

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

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

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

3 部分重疊的覆蓋形式

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

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

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

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

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

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

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

式中:Ui為定義在結(jié)點(diǎn)i上的多項(xiàng)式覆蓋函數(shù);Wi為權(quán)函數(shù)，

式中:ε0=εiε，η0=ηiη，(ε，η)為網(wǎng)格內(nèi)任意點(diǎn)的局部坐標(biāo);(εi，ηi)為結(jié)點(diǎn)的局部坐標(biāo)，見圖4。權(quán)函數(shù)滿足

將結(jié)點(diǎn)2，5，6“強(qiáng)制約束”到結(jié)點(diǎn)1，即令U2=U3=U4=U1，考慮到式(3)，則式(1)成為

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

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

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

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

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

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

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

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

5.1 位移函數(shù)

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

式中:權(quán)函數(shù)

其中:

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

5.2 剛度矩陣

流形元應(yīng)變矩陣為

式中，

剛度矩陣為

其中，

流形元?jiǎng)偠染仃嚨淖訅K(2×2)也可表示為

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

6 算例分析

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

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

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

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

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

圖6 上游面特征點(diǎn)x向應(yīng)力對(duì)比Fig.6 Comparison of horizontal stresses on the upstream surface calculated by manifold method and finite element method

圖7 上游面特征點(diǎn)y向應(yīng)力對(duì)比Fig.7 Comparison of vertical stresses on the upstream surface calculated by manifold method and finite element method

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

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

7 結(jié)語(yǔ)

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

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

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

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

［2］彭自強(qiáng)，葛修潤(rùn).數(shù)值流形方法在八結(jié)點(diǎn)有限元網(wǎng)格上的實(shí)現(xiàn)［J］.武漢大學(xué)學(xué)報(bào)(工學(xué)版)，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ù)值流形法分析溫度場(chǎng)及溫度應(yīng)力［J］.長(zhǎng)江科學(xué)院院報(bào)，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ù)值流形方法在流固耦合諧振分析中的應(yīng)用［J］.計(jì)算力學(xué)學(xué)報(bào)，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ù)值流形方法在結(jié)構(gòu)靜力分析中的應(yīng)用研究［J］.長(zhǎng)江科學(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))

［7］林紹忠，祁勇峰，蘇海東.基于矩陣特殊運(yùn)算的高階流形元單元分析［J］.長(zhǎng)江科學(xué)院院報(bào)，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.