蘇海東12

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

1 研究背景

工程設計離不開工程計算。在計算機出現(xiàn)以前，工程計算采用的是對結構進行簡化的材料力學或結構力學方法，并結合工程經(jīng)驗，形成設計規(guī)范。隨著計算機硬件和軟件技術的發(fā)展，以有限元法為主的數(shù)值計算已成為工程設計不可或缺的手段。水利工程是最早采用有限元法的領域之一，早在20世紀60年代，我國水利工作者就開始應用當時非常先進的有限元分析技術[1]。有限元法能夠應對各種復雜結構形狀、實際荷載作用、真實材料性質，從而得到比傳統(tǒng)的材料力學和結構力學方法更準確和更符合實際的計算成果。從重力壩、拱壩、廠房等大型水工結構，到閘墩、孔洞、襯砌等局部結構，以及鋼閘門、鋼梁板等薄壁結構，力學分析無所不在，有限元的主要分析技術都會涉及到。時至今日，從常規(guī)設計荷載作用下的線彈性靜力、動力分析，已發(fā)展到能細致模擬施工、運行過程和考慮開裂破壞、接觸問題等多種非線性的結構仿真分析[2-6]。

高壩大庫的建設使工程人員面臨復雜結構、復雜環(huán)境和條件[7]，因此，作為設計重要手段的水工結構分析技術也需要提升[2]。具體到有限元法，從分析的3個過程來看：前處理中，邊界附近的網(wǎng)格要準確模擬復雜的結構外形，這與網(wǎng)格形狀盡可能規(guī)則的要求相矛盾，再加上網(wǎng)格必須通過結點相連，這給不同區(qū)域的網(wǎng)格連接帶來不便，因此，有限元的網(wǎng)格剖分常常是一項非常繁重的工作；在計算過程中，有限元法對計算精度特別是應力精度的控制尚不如人意，比如，大壩的壩踵等局部應力復雜區(qū)域難以準確地模擬；后處理面臨大量的結果數(shù)據(jù)，整理和分析的工作量巨大。上述問題不僅給水工結構的有限元分析帶來很大難度，而且經(jīng)常導致計算結果的不穩(wěn)定，即使對于常規(guī)的線彈性靜力分析而言，不同的分析人員也很難獲得一致的計算結果(特別是應力)。結構仿真分析的難度更大，比如，施工澆筑塊之間的網(wǎng)格匹配，復雜環(huán)境條件下的荷載變化帶來的計算精度控制，裂縫擴展造成邊界變化而引起的網(wǎng)格重構等問題，都不易處理。因此，雖然有限元法的模擬精度高，但由于計算結果的不穩(wěn)定及方法運用上的困難，至今仍未以定量的方式進入設計規(guī)范[2]。

長江科學院在前期研究中提出一種新的數(shù)值計算方法——獨立覆蓋流形法[8]：采用多項式等完備級數(shù)的“分區(qū)級數(shù)解”逼近真實物理場；網(wǎng)格具有任意形狀、任意連接和任意加密的特性，可以使前處理自動化；位移及應力隨級數(shù)階次的增加而自然收斂，精度容易控制；分區(qū)的級數(shù)解答為后處理的數(shù)據(jù)自動提取和整理提供了方便。在此基礎上發(fā)展了一系列新的分析技術，展現(xiàn)出一些獨特的性能和優(yōu)勢。本文從計算方法的角度，討論獨立覆蓋流形法在水工結構分析中的應用前景。

2 獨立覆蓋流形法簡介

2.1 基礎理論簡介

1991年，石根華先生發(fā)明了數(shù)值流形方法[9-10]，將現(xiàn)代數(shù)學的流形思想引入數(shù)值計算：如圖1(a)所示，一系列相互重疊的覆蓋用于分割求解域，使復雜的物理場在每個覆蓋簡單化，以利于局部近似函數(shù)的逼近(數(shù)值計算實際上就是用設定的近似函數(shù)逼近真實物理場)。

從2011年開始，在石根華先生的建議下，筆者在數(shù)值流形方法中首次引入“獨立覆蓋”[11-12]，即覆蓋的獨立區(qū)域，如圖1(b)所示,每個覆蓋中的“藍色”區(qū)域就是獨立覆蓋，其近似函數(shù)通常采用多項式等完備級數(shù)(該級數(shù)隨階次升高能無限逼近任意的連續(xù)函數(shù))。在覆蓋之間的重疊區(qū)域內(nèi)(圖中的黃色區(qū)域)，采用有限元的線性形函數(shù)將各覆蓋中的級數(shù)聯(lián)系起來，以保證近似函數(shù)的連續(xù)性。

圖1 流形覆蓋示意圖Fig.1 Sketch of manifold covers

筆者進一步提出基于獨立覆蓋的數(shù)值流形方法，簡稱獨立覆蓋流形法，強調獨立覆蓋是計算分析的主要對象。獨立覆蓋占據(jù)一個覆蓋的主要面積，重點研究獨立覆蓋內(nèi)的級數(shù)，而覆蓋之間的重疊區(qū)域成為面積盡可能小的窄“條形”，只要不對線性方程組的性態(tài)造成嚴重影響，經(jīng)初步試算，可以小至僅占覆蓋面積的1%。由此得到整個求解域內(nèi)的各覆蓋分區(qū)的級數(shù)解，稱為“分區(qū)級數(shù)解”。

獨立覆蓋流形法的整體收斂是基于各覆蓋自身的局部收斂而建立的[8]。以通常采用的多項式級數(shù)為例，這相當于在各覆蓋區(qū)域內(nèi)將真實場函數(shù)展開為泰勒級數(shù)，或者說，各覆蓋內(nèi)的真實場函數(shù)(包括其導數(shù)，比如，位移的導數(shù)是應變或應力)被多項式級數(shù)無限逼近。數(shù)學上的級數(shù)逼近理論保證了新方法具有嚴格的收斂性，且這種收斂性不受覆蓋形狀的影響，也不受相鄰覆蓋在條形重疊區(qū)域是否錯位連接的影響。

2.2 網(wǎng)格劃分方式

新方法可應用任意網(wǎng)格劃分和規(guī)則網(wǎng)格覆蓋2種網(wǎng)格方式。

2.2.1 任意網(wǎng)格劃分

如圖2(a)所示，對求解域進行塊體劃分[13]。將塊體看作網(wǎng)格：各塊體具有任意形狀，即網(wǎng)格可有任意的邊數(shù)(含曲線邊)，各邊之間可為任意夾角；相鄰塊體的邊不必在線段的端點處對接，不同大小及形狀的塊體可隨意過渡，實現(xiàn)網(wǎng)格的任意連接；各塊體還可以根據(jù)計算精度的需要任意細分為更小的塊體，實現(xiàn)網(wǎng)格的任意加密。在圖2(b)中，自動生成塊體之間的條形連接。

圖2 求解域劃分為任意的塊體網(wǎng)格Fig.2 A domain divided into arbitrary blocks

2.2.2 規(guī)則網(wǎng)格覆蓋

另一種方式是用規(guī)則的網(wǎng)格去覆蓋求解域[11-12]。如圖3(a)所示的由矩形網(wǎng)格組成的獨立覆蓋系統(tǒng)，大的矩形網(wǎng)格為獨立覆蓋，細長的小網(wǎng)格為條形。將網(wǎng)格與物理邊界(圖中以橢圓表示)相互切割，形成網(wǎng)格內(nèi)的任意形狀的“流形元”。如圖3(b)所示，矩形覆蓋(圖中的f和g)還可以任意細分成更小的覆蓋，大、小覆蓋在條形處同樣任意搭接。

圖3 由規(guī)則的矩形網(wǎng)格組成的獨立覆蓋系統(tǒng)Fig.3 An independent cover system with regular rectangular meshes

圖2和圖3中，求解域的曲線邊界無需像通常有限元法那樣離散成折線，而是直接基于精確的幾何邊界(按曲線方程)進行計算[13-14]，這不僅實現(xiàn)了結構外形的準確模擬，而且為計算機輔助設計(CAD)和計算機輔助工程分析(CAE)的融合提供了基礎條件。

3 大體積水工結構分析的自動計算

在“分區(qū)級數(shù)解”的“分區(qū)”意義上，獨立覆蓋流形法的覆蓋網(wǎng)格突破了有限元法對網(wǎng)格的各種限制?！叭我庑螤睢边m應于復雜結構外形，加上“任意連接”的特性，適合于不同區(qū)域的網(wǎng)格連接，在此基礎上的“任意加密”使計算精度易于隨網(wǎng)格細化而提升。憑借著覆蓋網(wǎng)格的優(yōu)良特性，嘗試了完全自動的前處理，以及基于網(wǎng)格自動細分的h型自適應計算。

按照任意網(wǎng)格劃分方式，采用逐層分塊方法，將求解域逐步地由大塊網(wǎng)格細分為小塊網(wǎng)格。圖4以一個重力壩模型為例演示了這一過程：圖4(a)的單獨塊體，劃分為圖4(b)的2個塊體；每個塊體再作為單獨塊體，按同樣方式繼續(xù)細分；這樣在計算機自動操作下，一步步細分為圖4(g)的很多小塊網(wǎng)格。在圖4(g)的上游面作用靜水壓力，僅采用2階多項式，計算得到的表面壓應力與荷載作用值非常接近，表明這種任意形狀的網(wǎng)格在足夠的密度下可具有很高的計算精度。該方法可以模擬材料分區(qū)，不同的材料分區(qū)首先要劃分成不同的塊體，然后再進一步細化。

圖4 重力壩模型的任意網(wǎng)格細分Fig.4 Arbitrary mesh subdivision in a gravity dam model

每個塊體是否細分，要根據(jù)當前計算結果的誤差分析。獨立覆蓋流形法采用多項式等完備級數(shù)，不僅在理論上保證了隨級數(shù)階次穩(wěn)定逼近的p型收斂性，而且提供了誤差分析的方便性。初步提出了3種后驗誤差指標[15]：相鄰覆蓋的應力在條形重疊區(qū)域的連續(xù)性指標η1；自由表面上的應力計算值與實際值差異的邊界應力指標η2；獨立覆蓋內(nèi)部的應力誤差指標η3。研究表明，上述誤差指標能有效控制計算精度，甚至有望做到逐點的應力精度控制[16]。

圖5 在CAD中輸入結構 外形、邊界條件和計算 參數(shù)Fig.5 Input of structural outlines,boundary conditions and calculation parameters in CAD

圖6 CAE自動分析中的覆蓋網(wǎng)格加密過程及變形Fig.6 Procedures of cover refinements and deformation diagrams in CAE automatic computation

基于網(wǎng)格細分及誤差分析技術，嘗試了h-p型混合自適應，初步實現(xiàn)大體積水工結構二維分析的自動計算，及CAD與CAE的融合[17]。圖5至圖7按照規(guī)則網(wǎng)格覆蓋方式，以廠房混凝土結構某個剖面的靜力分析為例，演示了從CAD精確幾何建模和輸出，到CAE的自動前處理、自動分析和自動后處理的過程。所有的人工操作(包括結構幾何外形、材料參數(shù)、邊界條件的輸入)僅限于CAD中，而CAE分析過程無需人工參與，就可獲得滿足設定精度的計算結果。這種基于CAD精確建模和CAE自動分析的融合方式，代表著CAE未來發(fā)展方向。

圖7 自動輸出的應力渲染圖Fig.7 Stress rendering diagrams exportedautomatically

計算結果的級數(shù)表達，為后處理的數(shù)據(jù)自動提取與整理提供了條件。目前已能夠通過級數(shù)表達式找到最大應力點，下一步的研究將搜索應力超標范圍及可能開裂的區(qū)域，并根據(jù)設計需求進行斷面上的應力積分以獲取內(nèi)力，直至按照規(guī)范自動計算配筋量，使后處理的自動化程度更高。

關于水工結構三維分析的自動計算還在研究中，三維計算程序已基本完成，主要工作是前處理的三維塊體劃分或切割。

圖8 重力壩模型中的幾個奇異區(qū)Fig.8 Several singularregions in a gravity dammodel

4 裂縫和壩踵等奇異區(qū)的分析

自動計算技術可以保證結構大部分區(qū)域的計算精度。但對于應力奇異區(qū)，根據(jù)彈性力學原理，應力奇異點趨向于無窮大。在奇異區(qū)的網(wǎng)格越密，應力結果越大，因此不能直接用

應力作為強度判斷標準。圖8以重力壩模型為例，給出水工結構分析可能遇到的幾種奇異區(qū)：壩踵、壩趾以及邊界轉折處等<180°的凹角；如果出現(xiàn)裂縫，則奇異區(qū)位于裂縫尖端(簡稱裂尖，見圖8中的上游面附近)。如何準確模擬奇異區(qū)及整理計算結果，目前還沒有標準方法。

首先討論裂尖。根據(jù)線彈性斷裂力學(適用于混凝土等脆性開裂的材料)，裂縫的應力強度因子是反映應力奇異性的客觀指標，裂縫繼續(xù)擴展的條件是應力強度因子大于材料的斷裂韌度。

圖9 裂縫尖端附近的覆蓋網(wǎng)格Fig.9 Cover meshesaround a crack tip

獨立覆蓋流形法在“分區(qū)級數(shù)解”的“級數(shù)解”意義上也具有非常靈活的選擇。圖9給出裂尖附近的覆蓋布置，其中的獨立覆蓋5包含裂尖，采用Williams解析級數(shù)解，稱為解析覆蓋，周邊其他覆蓋仍采用多項式級數(shù)，稱為數(shù)值覆蓋，其中，獨立覆蓋4被裂縫分割成上、下2個獨立覆蓋41和42，以模擬裂縫兩邊的獨立運動。解析覆蓋和數(shù)值覆蓋通過條形的簡單連接方式，實現(xiàn)解析解和數(shù)值解的聯(lián)合求解[12,18]。

Williams解析級數(shù)是對裂尖位移場的最佳逼近，從而達成最快收斂，因此可以采用比有限元法至少大一個數(shù)量級的網(wǎng)格，同時，由級數(shù)系數(shù)可直接求得應力強度因子，不會引入額外誤差。該方法已應用于三維裂縫分析，包括Ⅰ型(張開型)、Ⅱ型(滑開型)和Ⅲ型(撕開型)裂縫的應力強度因子計算[19]。下一步的研究將模擬裂縫擴展過程，邊界形狀隨裂尖位置而改變，因此局部采用任意加密技術進行網(wǎng)格重構以保證計算精度。

再討論對重力壩而言非常重要的壩踵。以往同樣由于奇異性問題，有限元計算的壩踵應力無法直接利用，頂多將壩踵附近受拉區(qū)的范圍作為一種指標。工程上還采用“等效應力法”[20]：由建基面上游至下游的豎向應力積分而得的內(nèi)力(包括軸力和彎矩)，依據(jù)材料力學法計算壩踵應力。筆者認為上述方法仍帶有一定的經(jīng)驗性，因為基于受拉范圍或整個建基面合力的方式，很難真實反映壩踵點的局部奇異性。

事實上，依據(jù)彈性力學，凹角附近具有解析級數(shù)的表達式，且與凹角的角度相關，裂縫只是0°凹角情況。如圖10(a)所示，壩踵處的混凝土壩體和基巖，反映的是2種不同材料分區(qū)所夾的通常為90°的凹角情況，其解析級數(shù)的某些系數(shù)也是反映奇異性的客觀指標，稱為應力奇異性指數(shù)[21]。與裂尖類似，下一步將采用解析覆蓋與數(shù)值覆蓋聯(lián)合求解的方式，基于奇異性指數(shù)開展壩踵強度的研究，應該比現(xiàn)有方法更有說服力。

如圖10(b)所示，進一步推廣到任意凹角，計算圖8中的壩趾以及邊界轉角等，從而將所有奇異區(qū)的計算和成果整理納入標準模式，再針對奇異區(qū)開展誤差分析研究，從而使結構任何部位的計算精度都能得到有效控制。

圖10 具有解析級數(shù)的凹角區(qū)域Fig.10 Regions with concave angles where analyticalseries exist

5 拱壩分析與體形優(yōu)化

拱壩的結構分析與體形優(yōu)化，是拱壩設計最基本、最重要的工作。目前設計人員一般采用拱梁分載法[22-23]進行拱壩分析，優(yōu)點是計算速度快。但該方法將連續(xù)的拱壩用拱和梁離散，存在一定的近似，而且各家程序由于假設略有不同，計算成果也有一定差異。采用三維實體有限元建模，雖然模擬精度高，但過程相對復雜，計算量大，不利于需要反復計算的體形優(yōu)化工作，而且建基面上、下游的應力奇異性不好處理[24]。采用殼體有限單元建模是不錯的選擇，因其模擬精度比拱梁分載法高，又沒有實體有限元的計算量大及奇異性問題，但體形變化需要更改殼體單元的空間坐標甚至重構網(wǎng)格，而且只能做到對拱壩真實幾何形狀的近似。

獨立覆蓋流形法采用實體計算模式，應用實體計算程序直接求解梁板殼，僅通過挑選多項式級數(shù)項的簡單方式，就可以將梁板殼基本變形假設固化在近似函數(shù)中[25]。更重要的是，在曲梁和曲殼上根據(jù)中面方程建立局部坐標系，主要通過計算局部坐標關于整體坐標的導數(shù)，就能實現(xiàn)曲梁和曲殼的精確幾何，避免了建模中的幾何誤差[26-27]。

將新的曲殼計算方法應用于拱壩分析，中面方程定義為z=f(x,y)的多項式，厚度變化也用多項式表示。計算某實際拱壩，如圖11所示，僅采用40個獨立覆蓋，位移解與劃分4 500個殼單元的有限元解吻合得很好[28]，表明用很少的計算量就可以獲得高精度的計算結果。

圖11 拱壩順流向位移及與有限元結果的對比Fig.11 Downstream displacements of an arc damcompared with FEM results

圖12 應用獨立覆蓋分析和全局優(yōu)化進行拱壩體形優(yōu)化Fig.12 Optimization of the arc dam shape based onindependent cover analysis and global optimization

在拱壩分析的基礎上，提出了新的優(yōu)化分析思路[29]：只需一次性地在河谷平面內(nèi)建立二維計算模型，中面位置及厚度由多項式方程確定，無需在體形變化時更改空間坐標；需要優(yōu)化的參數(shù)就是這2個方程中為數(shù)不多的系數(shù)；基于完全公式化的幾何描述(結構外形及約束方程)和應力表達，并應用粒子群等算法實現(xiàn)全局優(yōu)化。圖12結合上述實際拱壩初步嘗試了體形優(yōu)化工作，優(yōu)化后，體積減少12%。下一步將考慮基巖或模擬基巖分層彈模變化，完善拱壩分析與體形優(yōu)化方法。

除了拱壩以外， 梁板殼的獨立覆蓋分析方法還能應用于水利工程中大量的薄壁結構， 特別是基于精確幾何的曲殼分析方法， 非常適合于弧形鋼閘門、 水電站壓力鋼管、 蝸殼鋼板等曲面鋼結構的分析。

6 水工結構仿真分析

水工結構設計一般只考慮常規(guī)荷載作用，如自重、水壓等主要荷載，而結構的實際受力情況要更為復雜，其中，由混凝土水化熱作用和環(huán)境溫度變化引起的溫度荷載是結構產(chǎn)生裂縫的重要原因，再加上水工結構如大壩是分層、分塊澆筑的，考慮結構外形變化及各時步荷載作用下的累積應力，與通?；谝淮涡越２⑹┘映Ｒ?guī)荷載的情況可能有較大差別。因此，需要模擬實際施工與運行過程進行仿真分析。

首先針對混凝土通水冷卻的溫度場展開研究[30]。在獨立覆蓋內(nèi)可布置水管，實現(xiàn)如圖13所示的水管周邊溫度分布的精細化計算。目前在水管所在的獨立覆蓋內(nèi)采用多項式級數(shù)，下一步將應用水管周邊的解析級數(shù)以達成更快收斂。

圖13 混凝土通水冷卻溫度場分析Fig.13 Temperature field analysis of water coolingfor concrete

然后針對一個澆筑塊開展瞬態(tài)溫度場的研究?？紤]到水化熱的空間分布是均勻的，而澆筑塊受環(huán)境溫度(氣溫、水溫)的影響，在距表面一定范圍內(nèi)呈現(xiàn)出相對復雜的變化，需要在此處劃分網(wǎng)格，這樣采用高階多項式級數(shù)計算出的澆筑塊溫度場具有很高的精度。

在此基礎上，考慮大壩的分層分塊澆筑，自動進行大壩的網(wǎng)格劃分，其中，上、下層的澆筑塊可單獨分塊，通過條形錯位搭接，如圖14(a)所示。由此進行了模擬重力壩澆筑過程的三維溫度場仿真分析的自動計算[31]：人工輸入結構外形、材料參數(shù)、大壩澆筑資料等必要信息后，整個溫度場仿真分析將自動完成，如圖14(b)所示。澆筑完畢后某個澆筑層沿順流向的溫度分布見圖14(c)，計算結果與采用細密網(wǎng)格的有限元解吻合良好，其中可見上、下游附近一定范圍內(nèi)受表面氣溫的影響。下一步將研究溫度徐變應力仿真分析的自動計算[32]。

圖14 重力壩溫度場仿真分析的自動計算Fig.14 Automatic simulation analysis for temperaturefield of a gravity dam

研究發(fā)現(xiàn)，溫度場受表面影響的范圍會隨時間而改變，因此到一定階段需要調整局部網(wǎng)格。同樣，在涉及材料非線性、接觸非線性、動力時程分析、開裂分析等分時步計算中，由于荷載、材料性質或結構外形變化造成原網(wǎng)格不能滿足當前計算精度的要求，也需要適時調整局部網(wǎng)格。獨立覆蓋流形法在前一時步的分區(qū)級數(shù)解答將作為當前時步的初始狀態(tài)，同時兼具網(wǎng)格劃分或調整的靈活性，這為仿真分析技術提供了新思路。

7 結 語

針對水工結構有限元分析中存在的網(wǎng)格剖分困難、計算精度不易控制、成果整理工作量大等問題，采用獨立覆蓋流形法發(fā)展了一系列新的分析技術：應用“分區(qū)級數(shù)解”的思想，基于嚴格的級數(shù)收斂理論和靈活的級數(shù)類型選擇，有望做到結構任意部位(包括應力奇異區(qū))的計算精度控制；再結合任意形狀、任意連接和任意加密的網(wǎng)格劃分以及成果自動整理，有望實現(xiàn)無需人工參與的自動計算，并為水工結構工程設計與工程分析的融合打下基礎；基于獨立覆蓋的梁板殼分析，也為精確幾何的拱壩分析與體形優(yōu)化、水工薄壁結構分析開辟了新的路徑；嘗試了模擬重力壩澆筑過程的溫度場仿真分析的自動計算，將來會擴展到水工結構全生命周期內(nèi)的自動仿真分析，并在分時步計算中提出隨時步更新局部網(wǎng)格的新思路。

當然，獨立覆蓋流形法提出的時間不長，本文介紹了部分已完成的工作，討論更多的還是一些新想法及其在水工結構分析中的應用前景。該方法的收斂理論還需完善以實現(xiàn)所有計算的精度控制，相應的算法還要改進以提升計算效率，應用領域也需擴展到各種非線性分析。希望借助新方法，提供高精度和穩(wěn)定可靠的工程計算成果，并最大程度地減少人工工作量，實現(xiàn)水工結構分析的標準化、規(guī)范化和自動化。

致謝：獨立覆蓋流形法的研究源于石根華先生提出的“將有限元網(wǎng)格拉開成獨立覆蓋及覆蓋重疊區(qū)域”的思想，文中引用了長江科學院數(shù)值流形方法研究團隊林紹忠教授級高級工程師、頡志強高級工程師、祁勇峰高級工程師、龔亞琦高級工程師、陳積瞻碩士、袁笑晨碩士、周朝碩士、李義碩士的工作成果，在此一并表示感謝。