第三节 数学模拟

数学模拟就是根据物理学定律，建立数学物理方程，构成工程的数学模型，依照规定的初始条件和边界条件，求解工程对外界作用的响应，并由此解决工程问题。其优点表现在：①数学模型易于构建，且能进行施工和运行过程的仿真分析；②数学模拟能突出构成建筑物本质特征的因素，便于分析了解建筑物的性能；③可以变动模型的有关因素，进行敏感性分析，了解它们对建筑物影响的程度及趋势，为改进设计提供启示。在计算技术飞速发展的今天，“计算”已经和“理论”、“实验”并列，被普遍认为是人类认识自然的三大支柱之一。常用的求解方法有以下几种。

（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）可改写为

求解过程与方法不变。