第三節 數學模擬

數學模擬就是根據物理學定律，建立數學物理方程，構成工程的數學模型，依照規定的初始條件和邊界條件，求解工程對外界作用的響應，并由此解決工程問題。其優點表現在：①數學模型易于構建，且能進行施工和運行過程的仿真分析；②數學模擬能突出構成建筑物本質特征的因素，便于分析了解建筑物的性能；③可以變動模型的有關因素，進行敏感性分析，了解它們對建筑物影響的程度及趨勢，為改進設計提供啟示。在計算技術飛速發展的今天，“計算”已經和“理論”、“實驗”并列，被普遍認為是人類認識自然的三大支柱之一。常用的求解方法有以下幾種。

（1）解析法。直接按模型的數學物理方程推導答案的解析式，是嚴格的理論解。但只能對較少邊界條件和初始條件簡單的典型課題才有解答，對工程適應性差。解析法成果精確，容易了解問題的規律性，可用于考核其他方法的精確性，是一個基礎性的方法。

（2）差分法。對于模型的數學物理方程，用差分運算代替微分運算，是近似的解法。差分法用于求解結構物的力學問題時優點不明顯，應用較少；但在水力學中仍是主要的算法，已逐步從科學研究的手段成為工程實用的算法。

（3）有限單元法。有限單元法是用離散的有限單元體逼近模擬連續體，在力學模型上是近似的，但在數學解法上是嚴格的。有限單元法發展迅速，不僅能求解位移場、應力場，也可以求解溫度場和滲流場；不僅能解決彈性問題，也可以解決彈塑性問題、動力學問題或巖體力學問題，在水工設計中的應用日趨廣泛。

（4）邊界單元法。邊界單元法的基本思想，即用積分方程方法解微分方程的思想可以追溯到20世紀初。邊界單元法是在綜合有限單元法和經典邊界積分方程方法的基礎上發展起來的，它把有限元的離散技巧引入經典的邊界積分方程方法中，通過一個滿足場方程的奇異函數——基本解作為權函數，將區域積分化為邊界積分，并在邊界上進行離散處理。其主要特點是：①將問題的維數降低一階，從而使得數據準備工作量及求解自由度大為減少；②由于離散僅在邊界上進行，故誤差只產生在邊界上，區域內的物理量仍由解析公式求出，因此具有較高的計算精度；③計算域內物理量時，無需一次全部求得，只需要計算給定點的值，從而避免了不必要的計算，提高了效率；④對應力集中、無限域等問題，該方法尤為適用。當然，也存在一些缺點：①得出的線性方程組的系數矩陣是一個滿的、不規則的矩陣，不便于應用已在有限元中發展成熟的處理稀疏對稱陣的線性代數方程組的一系列有效解法；②當問題的規模較大時，占內存較多，效率相對較低；③應用邊界單元法必須事先知道所求問題的控制方程的基本解，但從目前來看，非線性問題的基本解不易得出；④當物體嚴重不均質時將會大大影響邊界單元法的應用范圍和效果。

（5）DDA法。非連續變形分析方法（DDA）是平行于有限單元法的方法，但所有單元是被事先存在的不連續面所包圍的塊體。DDA法的單元或塊體可以是任意凸或凹形狀，甚至可以是帶孔的多邊形；而有限單元法限定只能用標準形狀的單元。此外，在DDA法中，當塊體接觸時，庫侖定律可用于接觸面，而聯立平衡方程式是對每一荷載或時間增量來選擇和求解的。在有限單元法的情況下，未知數是所有結點的自由度之和。在DDA法的情況下，未知數是所有塊體的自由度之和。

（6）流形方法。數值流形方法是利用現代數學——“流形”的有限覆蓋技術建立起來的一種新數值方法。有限覆蓋由物理覆蓋和數學覆蓋所組成，它可以處理連續問題和非連續問題。有限元在流形方法中只有一個單一的物理覆蓋，它覆蓋了全部數學覆蓋；DDA在流形方法中，則有許多物理覆蓋，它們各自覆蓋一部分數學覆蓋。這兩種方法在數值流形方法中只是兩個特殊的例子。在數值流形方法中，只要用不同覆蓋組合，就可以解決比有限元和DDA更普遍的復雜問題。

以上介紹了數值分析的一些方法。在構造數學模型（建模）時，應當做到：①合理的抽象。能有效地模擬建筑物的物理力學特性，明確表達各部分之間的關系。②精選主要因素。只引入能反映問題本質的核心因素組成計算模型。③適度的離散。在關鍵部位網格剖分較密，能據以評判建筑物的安全狀況，在不明顯影響關鍵部位的外圍部位，將單元剖分粗疏一些，以降低數值計算的維數。

一、彈性應力應變場的有限單元法

有限單元法是一種通用方法。利用有限單元法進行計算時可以考慮壩內的薄弱區，不同強度的混凝土分區，接縫和裂縫的影響，壩和地基的接觸區，壩和地基內滲透力的作用，復雜的地質構造和復雜的幾何形狀等。20世紀60年代以后，一些數學家發現，有限單元法實質上是變分法中里茲（Ritz）法的變種，從而使有限單元法的應用從求解應力場擴大到求解電磁場、溫度場和滲流場等，成為一種綜合能力很強的數值計算方法。隨著計算機技術和軟件工程的發展，其功能又有很大進步，如網格自動剖分、計算成果的自動整理、繪圖、屏幕顯示和光筆的應用等，可以使計算工作從過去繁瑣的勞動中擺脫出來，實現計算工作的自動化。

有限單元法的發展，借助于兩個重要的工具：在理論推導中采用了矩陣方法，在實際計算中采用了計算機。有限元離散、矩陣方法、計算機實現是相輔相成的。

目前能用于水工結構應力分析的三維有限元計算程序很多，各有特色。其中通用的大型結構分析程序，如SAP、ADINA、ANSYS等；專用的拱壩應力分析程序，如ADAP、AFED等。各程序的功能和適用范圍有所不同，可根據實際需要選擇不同的程序。

彈性問題有限單元法分析大體包含以下內容。

1.結構離散

先把連續體進行離散化，變為有限個在結點鉸接的單元組合體（圖4-1）。

最簡單的二維單元可取三角形單元。三維單元類型較多（圖4-2），其中最簡單的三維單元可取四面體單元。一般地，曲棱單元更能適應壩體體形，結點數較多，計算精度較高，但計算工作量相應地也較大（表4-1）。

此外，在拱壩分析中，還會用到殼體單元，如圖4-3和表4-2所示。殼體單元適用于拱壩壩體，尤其是較薄的拱壩。

2.位移模式

有限單元法中，最早提出并應用最廣的是位移法。取結點位移為基本未知量，求出結點位移后再求應力。

根據不同的單元型式，位移函數可以假設為一次式或高次式，次數越高，逼近真解越好，但計算工作量也越大。

圖4-4所示的平面問題的三角形單元ijm有6個結點位移分量，即

同時，三角形的三個結點上存在著結點力

在單元內部任意一點的水平位移u和垂直位移v的一般表達式為

式中：{δ}^e為單元結點的位移列向量；［N］為形狀函數矩陣，是坐標x、y的函數。

對于小變形問題，彈性理論中的幾何方程為

將式（4-8）代入幾何方程式（4-9），可以得出關系式

式中：［B］為單元的應變矩陣。

3.本構關系

本構關系表示材料的應力應變關系。

根據彈性理論中的廣義虎克定律，可以得出單元的應力與應變的關系式

｛σ｝=［σ_x σ_y τ_xy］^T

{ε}=［ε_x ε_y γ_xy］^T

式中：｛σ｝為應力列陣；｛ε｝為應變列陣；［D］為彈性矩陣。

將式（4-10）代入式（4-11），得

式中：［S］為應力矩陣。

4.虛功原理

設在單元e內產生虛位移{δ^*e}，相應的虛應變為{ε^*}，根據虛功原理，有

令｛ε^*｝=［B］｛δ^*e｝

｛ε^*｝^T=｛δ^*e｝^T［B］^T

｛σ｝=［D］［B］｛δ｝^e

代入式（4-14），得到

因虛位移{δ^*e}為任意向量，故由式（4-15）可得

式中：［K］^e為單元剛度矩陣。

求得單元剛度矩陣［K］^e和荷載向量{F}^e后，利用疊加原理可得出整個結構物的剛度矩陣［K］和荷載向量{F}，且滿足關系

解得結點位移后，代入式（4-12）即可求出應力。

二、滲流場的有限元分析

二維穩定滲流場的偏微分方程式及邊界條件可以由式（3-26）及式（3-28）～式（3-30）得到

根據變分原理，其解等價于下述泛函J（H）在滲流區D內求極值，即

滿足式（4-19）的函數H便是所求滲流場在定解條件下的解。該滲流場可通過離散化，用有限個單元的集合體近似求解。同樣以三角形單元為例，假定單元內水頭函數H為線性函數，即

設三角形單元三個結點的水頭值分別為H_i、H_j、H_m。式（4-20）中的三個常數α₁、α₂、α₃可用三個結點i、j、m的坐標及相應的水頭值H_i、H_j、H_m表示。于是，水頭函數

式中：［N］為形函數矩陣。

單元內水頭函數H的導數為

式中：A為三角形單元的面積；b_i、b_j、b_m、c_i、c_j、c_m為三角形單元結點坐標的函數。

求單元泛函的微分，經演算得

式中：［K］^e為單元傳導矩陣，與單元結點坐標及滲透系數有關。

對所有單元的泛函求得微分后疊加，并使其等于零（求最小值），得到泛函對結點水頭的微分方程組

式中：n為結點數。該線性方程組具有n個方程式。

寫成矩陣形式，得穩定滲流的計算公式為

式中：［K］為總傳導矩陣，由各個單元傳導矩陣［K］^e組成；{H}為結點水頭列陣；{Q}為結點流量列陣，由已知邊界條件得出。

求解總體方程組式（4-26），即可得出各結點的水頭值，然后得出流速、流量等其他滲流要素。利用這些要素可以整理出如圖3-11和圖3-12所示的流網圖，也可以整理出流速場，如圖4-5所示。

三、溫度場的有限元分析

1.溫度場計算

大體積混凝土結構必須研究并控制其溫度變化及溫度應力。在混凝土入倉后，由于水化熱的作用，使混凝土的溫度持續上升，直至達到最高溫度后，開始回降，在回降過程中產生的溫度應力，是引起混凝土裂縫的重要因素之一。既然溫度應力是由于溫度變化而產生的，那么分析工作可以分兩步走：第一步研究壩內溫度場的變化過程和最終穩定狀態，第二步計算溫度應力場。

比較穩定溫度場的偏微分方程（3-22）和穩定滲流場的偏微分方程（3-27），以及對應的溫度邊界條件式（3-20）、式（3-21）和滲流邊界條件式（3-28）、式（3-29）可以發現，溫度場的控制方程與滲流場相似。用有限單元法求解溫度場時，滲流場的公式完全適用。

不穩定溫度場的計算可參考《水工設計手冊》（第二版）等相關著作。

2.溫度應力計算

物體的溫度變化產生體積變形，當變形受到外部或內部約束時便產生溫度應力。設三角形單元3個結點的溫度變化為ΔT_i、ΔT_j、ΔT_m，則單元的平均溫度變化為

在無約束條件下，由ΔT引起的初應變為

式中：α為線脹系數，約為1.0×10^-5，1/℃。

在約束作用下，初應變{ε₀}不能自由產生，設實際產生的應變為{ε}=［ε_x ε_y γ_xy］^T，則單元的溫度應力為

式中：［D］為彈性矩陣；［B］為應變矩陣；｛δ｝^e為結點變位列陣。

由虛功原理可得因溫度變化產生的單元等效結點荷載為

式中：t為單元厚度，m；A為單元面積，m²。

對所有單元疊加，得到

求解由溫度引起的結點位移后，再由式（4-28）計算單元的溫度應力。

四、動力分析

水工結構除了承受靜力荷載外，往往還承受動力荷載，如爆破沖擊、波浪壓力、地震作用等。尤其是地震作用，對水工建筑物的安全影響很大。

用有限單元法求解水工結構地震動力響應問題時，首先將結構劃分成若干單元，然后根據達朗貝爾原理，建立以矩陣形式表示的結點動力平衡方程

式中：［M］為結構的質量矩陣；［C］為結構的阻尼矩陣，如假設為比例阻尼情況，則［C］=α₀［M］+α₁［K］，α₀、α₁為常數；［K］為結構剛度矩陣；{u（t）}、{u.（t）}、{u¨（t）}分別為結構的位移列陣、速度列陣、加速度列陣；{u¨_g（t）}為地震時地面加速度列陣。

可用時程分析法、振型分解法等求解方程（4-31）。

時程分析法是由結構基本運動方程輸入地震加速度紀錄進行積分，求得整個時間歷程內結構地震作用效應的方法。此法對線性、非線性、單自由度、多自由度結構都能應用，適應性好，但工作量大。

振型分解法是先求解結構對應其各階振型的地震作用效應后，再組合成結構總地震作用效應的方法。此法適用于多自由度彈性結構的動力分析。各階振型效應用時程分析法求得后直接疊加，稱為振型分解時程分析法；用反應譜法求得后再組合，稱為振型分解反應譜法。

下面就振型分解反應譜法求結構的地震動力反應作一簡單介紹。

采用反應譜法需要知道結構的自振頻率和振型。實際上，結構的自振特性（自振頻率和振型）本身也是結構抗震研究的重要內容之一，計算經驗表明，結構的阻尼對自振特性的影響很小，在求結構的自振特性時，可略去阻尼的影響。

在基本方程（4-31）中，令［C］和{u¨_g（t）}為零，可得到結構無阻尼自由振動方程，即

設結構作自由振動時的特解為

式中：｛δ｝為振幅列陣。

式（4-33）表示各質體作同頻率ω和同初相角γ的簡諧振動。求出{u¨（t）}、{u（t）}代入式（4-32），兩邊同除以sin（ωt+γ），得

式（4-35）為齊次線性代數方程組，{δ}有非零解的條件是

對于二維問題，如結構有n個結點，則將式（4-36）展開，可得ω²的2n次方程，稱為特征方程。求解該方程，可得2n個頻率ω，稱為特征值；再把這2n個ω依次代入式（4-35），便可解出2n組振幅{δ}，即特征向量。但{δ}中各個元素只表示相對值，通常設{δ}中的最大一個元素為1，從而可以定出{δ}中其他元素的值。要求解全部的特征值和特征向量，可用Jacobi法，但通常我們只需求解前幾階或十幾階自振頻率和振型，因此可用迭代法或子空間迭代法。

當求出s個結構的自振頻率 ω₁、ω₂…ω_s和s個振型{δ}₁、{δ}₂…{δ}_s后，將結構的位移按振型展開，得

式中：Y₁（t）、Y₂（t）、…、Y_s（t）為待求的主坐標。由式（4-37）得出{u.（t）}和{u¨（t）}，并且和{u（t）}一起代入式（4-31）。利用振型正交條件，可將動力方程（4-31）分解成一組關于主坐標Y_i（t）的非偶合微分方程，其中第i個方程為

式中：、為振型參與系數；ζ_i為阻尼比。

式（4-38）等號右邊反映了地震作用，該式和單質點體系強迫振動方程完全一致，其解Y_i（t）就是自振頻率為ω_i、阻尼比為ζ_i的單質點體系的地震反應。由式（4-37），結構的動力位移為

Y_i（t）的大小反映了第i階振型在強迫振動中所占成分的大小，如果地震作用與阻尼都相同，而以不同的ω_i值代入式（4-38），解出的Y_i（t）則不同，所以Y_i（t）的大小主要取決于ω_i。在一定的阻尼下，以ω為橫坐標，Y_i（t）為縱坐標，可作出一系列離散的分布，這些點的全體稱為位移反應譜。如以為縱坐標，以ω（或T=）為橫坐標，同樣可以分別作出速度反應譜和加速度反應譜。但Y（t）本身都是時間的函數，對于不同的時刻有不同的數值，由此作出的反應譜對于每個時刻都不一樣，而工程上所需要的是最大的地震反應。所以，工程上以Y_imax、為縱坐標，以ω（或T）為橫坐標，分別作出各種反應譜。有了反應譜，就不需要求解式（4-38）的常微分方程，可直接從各種反應譜中找出相應于第i階振型的反應譜值Y_imax、，再按振型疊加后，就可得到結構最大的地震反應（位移、速度、加速度等）。

DL 5073—2000《水工建筑物抗震設計規范》提供了水工建筑物動力法計算的設計反應譜曲線，如圖4-6所示。圖中β（T）為動力系數，其原始定義為單質點彈性體系在水平地震作用下水平反應絕對加速度最大值與地面最大水平加速度之比。

設計加速度反應譜是重要的地震動參數，其形狀及有關參數與所在場址的巖土類別、場址離地震震中的距離有關。圖4-6中，β_max為設計反應譜最大值的代表值，按表4-3選取；T_g為設計反應譜特征周期，與地基類別有關，巖石地基取0.2s，一般非巖性地基取0.3s，軟弱地基取0.7s。設計烈度不大于8度且基本自振周期大于1.0s的結構，特征周期宜延長0.05s；β_min為設計反應譜下限值，β_min應不小于設計反應譜最大值的代表值的20%。

根據設計反應譜，相應于第i階振型的位移為

式中：c為綜合影響系數，與式（3-70）中的地震作用效應折減系數ξ類似，既反映結構彈塑性帶來的抗震潛力和結構對地面的抑制，也包括施工質量和運行條件等因素，一般取，重要結構取；a_h為水平向設計地震加速度，可查表3-12得到；β_i可根據第i階周期T_i查圖4-6得到。

各階振型應力為

由于各階振型最大值反應一般不在同一時刻出現，因而存在各階振型反應如何組合的問題。常用平方和方根法組合（即SRSS法），得總的動力應力為

當兩個振型的頻率差的絕對值與其中一個較小的頻率之比小于0.1時，地震作用效應宜采用完全二次型方根法組合（即CQC法），這時

式中：S_E為地震作用效應；S_i、S_j分別為第i階、第j階振型的地震作用效應；s為計算采用的振型數；ρ_ij為第i階和第j階振型的相關系數；ζ_i、ζ_j分別為第i階、第j階振型的阻尼比；γ_ω為圓頻率比，γ_ω=ω_j/ω_i；ω_i、ω_j分別為第i階、第j階振型的圓頻率。

地震作用效應影響不超過5%的高階振型可略去不計。采用集中質量模型時，集中質量的個數不宜少于地震作用效應計算中采用的振型數的4倍。

對于擋水建筑物，滿庫時地震，水工建筑物與庫水一起振動，此時利用附加質量矩陣［M_P］，動力方程（4-31）可改寫為

求解過程與方法不變。