2.2 筑坝材料动力本构模型及计算方法
2.2.1 坝料动力本构模型及其特点
由于岩土的实际动力本构关系极为复杂，它在不同的荷载条件、土性条件及排水条件下会表现出极不相同的动本构特性，要建立一个能够适用于各种不同条件的动力本构模型普遍公式是比较困难的。目前，具体建立的动力本构模型已多达数十种，大致可以分为两大类，即非线性弹性模型和弹塑性模型。
非线性弹性模型可以分为物理类模型和经验类模型。物理类模型是用一系列具有不同屈服强度的滑块和不同初始刚度的弹簧来描述土的动力本构关系，它实质上是弹塑性模型中多屈服面模型的基础，可以较好地表达土的滞回特性，但计算得到的阻尼比实测值偏小。因物理类模型参数较多，在实际应用中很少采用。在试验数据基础上建立和发展起来的经验类模型可以分为等效线性模型和真非线性模型。等效线性黏弹性模型在目前土石坝动力反应分析中被广泛采用，此种分析方法比较简便，计算分析的稳定性好，能够合理地确定土体在地震过程中的加速度、剪应力和剪应变幅值。真非线性模型是根据不同加载条件、卸载-再加载条件直接给出动应力-应变的表达式。中国水利水电科学研究院以Masing准则为基础，研究了一种基于剪应力比控制的循环三轴试验的非线性粘弹塑性模型，即中国水科院真非线性模型。
弹塑性动力本构模型在理论上相对更为合理，能较好地反映土体的实际状态，并能够计算静、动力全过程的应力变形以及直接计算坝体的永久变形。Zienkiewicz和Mroz提出了广义塑性力学的基本思想，随后Pastor和Zienkiewicz对其基本框架进行了扩展并基于该理论建立了适用于黏土和砂土的Pastor-Zienkiewicz本构模型（简称广义塑性模型）。我国学者考虑应力相关性和颗粒破碎状态相关性对砂土广义塑性模型进行了改进，发展了筑坝堆石料的改进广义塑性模型并进行了验证，可应用于面板堆石坝弹塑性静、动力分析。清华大学针对循环荷载下堆石料的力学特性进行了较为系统的研究，建立了弹塑性循环本构模型。
1.等效线性黏弹性模型
在等效线性模型中，将土视为黏弹性体，不寻求滞回曲线的具体数学表达式，而是以等效剪切模量G和等效阻尼比λ作为动力特性指标进行计算。实际工程中一般根据试验曲线确定动剪切模量和阻尼比，然后根据试验结果给出等效剪切模量G和等效阻尼比λ与动剪应变的关系，即可根据动剪应变确定其对应的等效剪切模量G及等效阻尼比λ，进而进行坝体地震动力有限元分析。
式中 Geq和λeq——动剪切模量和阻尼比的等效值;
Gmax和λmax——动剪切模量和阻尼比的等效值的最大值;
γm——滞回圈动应变幅值;
γr——参考剪应变。
等效线性黏弹性模型由于是试验结果的归纳，形式上比较简单，因此在实际计算中得到广泛应用。但它客观上存在着多方面的不足，如不能考虑应力路径的影响，不能考虑土的各向异性，在应变较大时误差大，同时也不能直接考虑土石坝的地震永久变形。
2.改进的堆石料广义塑性模型
广义塑性模型具有许多优点，包括：不需要定义塑性势面函数直接确定塑性流动方向，不需要定义加载面函数直接确定加载方向，不需要依据相容性条件直接确定塑性模量，可以考虑剪胀和剪缩以及循环累计残余变形。此外，广义塑性模型框架清晰，便于在有限元程序中实现，用一套参数即可完成土工建筑物的静、动力分析过程。即广义塑性模型不仅适用于土工构筑物的施工填筑过程，也适用于地震动力响应分析，且可以直接计算地震永久变形。
广义塑性模型提出时主要针对砂土的液化问题。砂土液化分析时围压的变化范围较小，而由于高土石坝坝体内部平均主应力的变化范围较大，广义塑性模型在考虑压力相关性时其参数受平均主应力的影响较大，因此，该模型在高土石坝静、动力分析方面的应用存在一定的局限性。我国学者在弹性模量、加载模量和卸载模量方面考虑了筑坝材料的应力相关性，对广义塑性模型进行了改进。
改进的广义塑性模型共有17个参数，G0、ms为初始剪切模量参数；K0、mv为初始体积模量参数；αg、αf、Mg、Mf为塑性流动方向和塑性加载方向相关参数；H0、ml、β0、β1为塑性模量中加载及再加载的相关参数；Hu0、mu、γd、γDM、γu为卸载模量的相关参数。根据紫坪铺大坝筑坝堆石料静、动力试验结果，确定了模型参数，对固结排水剪试验和循环荷载试验关系曲线的模拟见图2.2-1和图2.2-2。可以看出，改进的广义塑性模型能够较好地反映堆石料的剪胀性、循环累计塑性应变、循环致密及滞回特性。结合不同岩性的筑坝材料在等σ1、等σ3、等p以及等应力比Kc应力路径下的大型三轴试验和侧限压缩试验成果，通过验证得出改进的广义塑性模型对复杂应力路径也具有很好的适应性。
图2.2-1 紫坪铺大坝筑坝堆石料静力试验应力-应变关系
图2.2-2 紫坪铺大坝筑坝堆石料循环荷载试验应力-应变关系
3.堆石料循环弹塑性模型
清华大学针对堆石材料循环本构关系进行了较为系统的研究，认为剪切和压缩各引起一个剪应变和一个体应变，每个剪应变和体应变又可分为可逆和不可逆的两种情形，将应变分解为8个分量，建议了应变分解方法，提出了一个16参数的弹塑性循环本构模型和对其进行简化后的9参数简化模型，主要针对高堆石坝材料动应力应变的基本规律及特点，做了以下几点改进：
（1）对剪切作用引起的不可逆体应变，采用了新的描述方法，主要是为了更合理地描述循环剪切荷载作用下堆石料体积累积收缩的特有规律，参数可根据常规动三轴试验得到的残余变形与荷载循环作用次数的关系整理获得。
（2）反映了反复压缩引起的残余体应变，这是以往的动本构模型所不能反映的。
（3）可以合理地描述堆石料的体积屈服特性。
（4）对模型提出了简化的方法和数值化方法，根据需要，可应用完整模型进行静动力联合计算，也可用简化模型在一定初始条件下仅进行动力计算。
采用系列化多种应力路径动三轴试验（包括常规循环、等p循环、等向压缩、径向同步）的成果，对本构模型的预测能力进行了验证。试验验证表明，该模型对高土石坝中常见应力路径下的循环应力应变响应具有相当好的适应性。
以古水面板堆石坝为计算实例，实现了基于该本构模型的数值积分算法，通过编制计算程序、进行验证性计算，表明了该模型及算法的有效性。古水面板堆石坝在设计地震作用下，坝体地震动力响应及残余变形分布见图2.2-3和图2.2-4。
图2.2-3 坝体的地震加速度响应（单位：m/s2）
4.真非线性动力本构模型
目前土体地震反应的真非线性动力模型大都以 Masing准则为基础进行补充和改进，其不足之处主要有：①模型中的滞回圈和骨干曲线与振动次数无关，没有反映土体应变历史的影响；②在周期荷载作用下，模型给出封闭的滞回圈，而土体的实际变形规律是不封闭的，而且模型滞回圈包围的面积比实测的面积大；③在不规则循环荷载作用下，当土体承受的剪应力比超过或等于历史上最大剪应力比时，使用骨干曲线表达土体的动力变形特性与实际情况相差很大。为此中国水科院研究了一种基于剪应力比控制的循环三轴试验的真非线性黏弹塑性模型。
图2.2-4 地震后残余变形分布（单位：m）
该三维真非线性动力本构模型的特点为：①与等效线性黏弹性模型相比，能够较好地模拟残余应变，用于动力分析可以直接计算残余变形；在动力分析中可以随时计算切线模量并进行非线性计算，这样得到的动力响应过程能够更好地接近实际情况。②与基于Masing准则的非线性模型相比，增加了初始加荷曲线，对剪应力比超过屈服剪应力比时的剪应力应变关系的描述较为合理；滞回圈是开放的，能够计算残余剪应变；考虑了振动次数和初始剪应力比等对变形规律的影响。
鉴于该非线性黏弹塑性模型的特点，为了更有效地进行真非线性动力反应分析，可以采用增量法和全量法交替进行的算法以控制增量法的误差积累。根据非线性黏弹塑性模型及有限元原理，推导出结构的增量和全量方程分别为
式中 {u}、
——节点位移、速度和加速度；
{ue}——弹性位移；
Δ——增量；
［M］——质量矩阵；
［C］t和［C］s——切线和割线阻尼矩阵；
［K］t和［K］s——切线和割线刚度矩阵；
{Fa}——地震力；
{Fe}——应力超过强度时加以修正的等价节点力（超越力）。
具体求解按增量步进行。对每一增量步，先求解式（2.2-3），然后如果为奇数增量步，则在假定 {
}不变的条件下，由式 （2.2-4）计算弹性位移{ue}；如果为偶数增量步，则在假定 {ue}不变的条件下计算加速度 {
}，并用此加速度校正式 （2.2-3）中的 {Δ
}，以减少用增量法解方程产生的误差积累。其他如孔隙水压力的消散和扩散计算等采用前述的方法。
针对一坝高超过250m的典型高面板堆石坝，分别采用真非线性与等效线性分析方法进行动力反应分析，对比真非线性分析与等效线性分析的差异。
（1）坝体加速度反应。图2.2-5和图2.2-6分别为采用真非线性模型与等效线性模型计算的大坝典型剖面的加速度反应等值线。
图2.2-5 真非线性模型算得的大坝典型剖面的反应加速度等值线（单位：m/s2）
从图中可见，采用真非线性模型与等效线性模型算得的大坝典型剖面的反应加速度在量值上较为接近，比较而言，真非线性模型算得的反应加速度要大一些，尤其是坝顶和坝坡的放大效应更为明显；而且真非线性模型算得的坝体下部的反应加速度比等效线性的反应更为“充分”，约大5%~10%。
（2）面板动应力反应。图2.2-7和图2.2-8分别为采用真非线性模型与等效线性模型计算的面板顺坡向动应力等值线图。
图2.2-6 等效线性模型算得的大坝典型剖面的反应加速度等值线（单位：m/s2）
图2.2-7 采用真非线性模型算得的面板顺坡向动应力的等值线（单位：MPa）
图2.2-8 采用等效线性模型算得的面板顺坡向动应力的等值线（单位：MPa）
从图中可见，采用真非线性模型与等效线性模型算得的面板顺坡向动应力的等值线图在分布上比较接近，但在量值上差别较大。比较而言，真非线性模型算得的面板动应力要大一些，约大10%~15%。采用真非线性模型可以在动力分析过程中直接计算地震残余变形，可以在一定程度上体现地震过程中地震永久变形对坝体地震反应，包括对面板应力的影响，从理论上更为合理。
2.2.2 地震永久变形分析方法
永久变形分析方法主要分为两类：一类是确定性永久变形分析，包括滑体变形分析和整体变形分析等；另一类是非确定性永久变形反应分析，包括滑体位移随机反应分析和整体位移随机反应分析等。目前实际工程中主要采用整体变形分析法进行土石坝永久变形分析，其按永久变形产生的机理不同可分为简化分析法、软化模量法、等价节点力法和等价惯性力法等四种方法，其中土石料动应力和残余应变关系的研究是整体变形分析法的关键。
2.2.2.1 Newmark滑体变形分析法
土石坝地震永久变形分析，首先是由Newmark基于屈服加速度的概念提出的，按刚塑性理论进行坝坡永久变形计算。认为当某一滑动体的加速度超过材料的屈服加速度时，沿破坏面就会发生滑动，向下滑的位移被认为是不可恢复的永久位移。这种方法计算出的永久位移大小与地震过程中出现短暂失稳的时间（次数）和失稳的程度有关，实质上是对坝坡动力稳定性的一种评价，没有体现土体的应力应变关系，尤其是残余体积变形的影响。该类方法只适用于孔隙水压力不会明显升高的土体，对于孔隙水压力会明显升高或强度会明显降低的土体，计算结果偏小。
2.2.2.2 整体变形分析法
图2.2-9 简化分析法流程示意
这类方法是将地震前后坝体及坝基均假定为连续体，先通过室内试验得到残余应变模型，再按照连续介质的有关理论来进行计算。
（1）简化分析法。简化分析法最早是Seed在1973年提出的，该法根据坝体动力分析结果，结合筑坝材料的动应力和残余应变关系曲线，确定坝体的平均残余剪切应变势，乘以坝高，近似估计坝顶的水平残留位移。简化分析法的基本流程见图2.2-9。
（2）软化模量法。软化模量法是基于“软化模型”概念提出的，该模型认为土体单元的残余应变势是由地震荷载作用引起的土料发生软化产生的，坝体永久变形等于坝体在动荷载作用前的静变形与坝体震动软化后的静变形之差，其力学模型见图2.2-10。
图2.2-10 软化模量法力学模型
软化模量法根据静力试验数据确定土体的初始模量Ei，并采用一次性加载的方法计算出土体地震前的应变εi。地震前土体初始应力σi与初始应变εi之间的关系式为
由于坝体中大部分为饱和土体，故地震时可假设土的体积模量保持不变，则土体剪切模量就会降低。地震后的总变形可以采用一次加载方式根据降低后的模量Ef求出。两次静力计算的位移差即为土体的永久变形，具体流程见图2.2-11。
（3）等价惯性力法。等价惯性力法由Taniguchi、Whitman和Marr提出。该方法利用地震动力反应分析得到的坝体中各节点的等效水平加速度分布，推算出坝体各节点上的地震等效水平惯性力，并将此惯性力作为静荷载施加在坝体节点上，方向分别指向坝体上游和下游，然后按循环三轴试验中动应力与残余应变在一定等效循环周数下的无量纲关系曲线进行迭代计算，将得到的两个方向的永久变形进行线性叠加，即为坝体最终永久变形，其流程见图2.2-12。
图2.2-11 软化模量法流程图
（4）等价节点力法。等价节点力法最早由Serff在1976年提出，该方法认为地震引起的永久变形等于某种等价节点力作用下所产生的附加变形。等价节点力法根据土石料的动力试验建立应力状态、动应力幅值和循环振次与土体残余应变的关系式。又根据坝体的静动力有限元分析确定坝体各单元的围压、固结比、振次和动应力情况。通过上述动力试验和有限元分析，可以确定坝体各单元在地震过程中的残余应变势。将此残余应变势等效为一种静节点力施加于有限元网格节点上，作为荷载按静力法计算坝体变形，即地震引起的永久变形。等价节点力法的流程见图2.2-13。
图2.2-12 等价惯性力法流程图
图2.2-13 等价节点力法流程图
在整体变形分析法中，残余应变势的确定是关键。目前已有典型的残余应变势模型有沈珠江五参数模型及其修正模型、中国水利水电科学研究院模型等。