3.1 有限元数值分析理论基础
3.1.1 土石料的本构模型
堆石坝筑坝堆石料是非线性材料，变形不仅随荷载的大小而变化，还与加荷的应力路径相关，其应力应变关系呈现明显的非线性特性。堆石体的计算模型通常采用非线性模型和弹塑性模型等。各参建单位不断修改和改进这类大型计算程序，使用了不同的土的应力应变模型，如邓肯⁃张模型、成都科大的K⁃G模型、清华大学弹塑性模型以及沈珠江院士提出的双屈服面弹塑性模型等。
3.1.1.1 邓肯⁃张模型
邓肯⁃张模型公式简单，参数物理意义明确。邓肯模型是堆石坝中最常用的模型之一。三轴试验研究结果表明，其对土体应力应变非线性特性亦能较好地反映。邓肯⁃张E⁃μ模型以切线弹性模量Et和切线泊松比μt作为计算参数。其中，切线弹性模量表达式为
切线泊松比为
其中
式中：pa为单位大气压力；G、F为试验常数；D为假设的轴向应变εα渐进值的导数；S为剪应力水平，反映材料强度发挥程度，表达式为
式中：（σ1－σ3）f为破坏时的偏应力，由摩尔⁃库仑破坏准则得：
上述各式中：K、c、Rf、n、φ、G、F、D为模型参数，由常规三轴试验确定。
由于式（3.1）计算的Et值偏大，且在求取G、F、D时带有不确定性。邓肯等人又建议用切线体积模量Bt代替μt，即
对于卸载情况，采用卸载回弹模量Eur进行计算：
式中：Kur为试验常数，可通过试验测定。
3.1.1.2 成都科大的K⁃G模型
土的总应变dε等于剪切应变dεs与体应变dεv之和，即
如果考虑球张量与偏张量的耦合作用，纯粹的剪切过程可以产生体积变形，平均法向应力之和的改变也会引起剪切变形，则
式中：p＝（σ′1＋2σ′3）/3和q＝σ1－σ3（轴对称情况，即ε2＝ε3）为有效平均应力和广义剪应力；εvp和εvq分别为p和q引起的体应变；εsp和εsq分别为p和q引起的剪切应变。
对于一般应力水平下黏性粗粒土或高应力水平下的粗粒土，可忽略剪胀作用，故可得到成科大简化K—G模型的切线体变模量、切线剪切模量和考虑中主应力影响的破坏准则如下：
式中：Ki、αK、Gi、αG、βG、a和b为试验参数。
3.1.1.3 沈珠江弹塑性模型
沈珠江提出的新弹塑性模型，其应力⁃应变关系具有剑桥模型形式，其有关参数则象邓肯E⁃μ模型一样从普通三轴剪切试验的应力⁃应变⁃体变关系曲线中得到。新理论只把屈服面看作弹性区域的边界，不再把它与硬化参数联系起来，具体建议的双屈服面方程如下：
其中
式中：r和s为参数。
根据塑性应变增量理论，可以求得双屈服面理论的体积应变增量和剪切应变增量分别为
在三轴剪切试验状态下，Δp＝Δσ1/3；Δq＝Δσ1；Δεs＝Δεa－Δεv/3；定义Et＝Δσ1/Δεa和μt＝Δεv/Δεa，可解出A1和A2，有
其中
该模型特点是假定三轴剪切试验中q—εa的关系曲线为双曲线，而εv—εa关系曲线为抛物线，方程如下
则Et可按E—μ模型计算，μt则是代替E—μ模型中μt的切线体积比，则两者计算公式如下
式中：εvd、εad和（σ1－σ3）d分别为最大正体应变及其对应的轴应变和主应力差；Cd、d和Rd分别为代替E—μ模型中的F、D和G的3个参数，分别代表σ3＝1Pa时的最大体应变、体应变随σ3变化的幂次和εvd发生时的应力比Rf（σ1－σ3）d/（σ1－σ3）f。Ei为初始弹性模量；Rs＝Rf（σ1－σ3）/（σ1－σ3）f；其他参数意义同前。
3.1.1.4 清华大学弹塑性模型
清华弹塑性模型是以黄文熙教授为首的研究组提出来的。其特点在于不是首先假设屈服面函数和塑性势函数，而是根据试验确定的各应力状态下的塑性应变增量的方向，然后按照相适应流动规则确定屈服面，再从试验结果确定其硬化参数。
1.弹性应变
弹性应变部分采用K⁃G模型计算。其中体变模量K从各向等压试验的卸载曲线确定；剪切模量G从常规三轴压缩试验确定，其一般形式为
2.屈服面的确定
关于清华弹塑性模型的详细介绍，可参考李广信编《高等土力学》，在此不一一赘述。
3.1.2 材料非线性问题的有限元解法
对于材料非线性问题，采用的有限元方法主要有增量法、迭代法以及混合法等。增量法的一个突出特点是可以考虑逐级施加荷载，这不仅能反映出施工过程中各阶段的应力和变形情况，而且能体现结构本身随施工过程的变化，更好的体现材料的非线性，因而更适于堆石坝逐级加荷的情况。下面介绍一下增量法的原理。
增量法是将全荷载分为若干级微小增量，逐级用有限元法进行计算。对于每一级增量，在计算时假定材料性质不变，作线性有限元计算，求出位移、应变和应力的增量。各级荷载之间，材料性质变化，刚度矩阵变化，反映了非线性的应力应变关系。该方法实际上是用分段直线来逼近曲线。增量法分为以下几种。
3.1.2.1 基本增量法
各级荷载下的材料性质是由刚度矩阵D来体现的，无论弹性矩阵还是弹塑性矩阵都决定于应力状态，基本增量法是根据每级的初始应力来确定刚度矩阵D的，其计算步骤为：
（1）用前级终了时的应力，也就是本级的初始应力｛σ｝l－1，确定刚度矩阵Dl。对于弹性非线性问题，就是确定切线弹性常数Etl和μtl。
（2）由Dl形成刚度矩阵Kl。
（3）解线性方程组Kl｛Δδ｝l＝｛ΔR｝l，求出位移增量｛Δδ｝l，相应的位移总量为｛δ｝l＝｛δ｝l－1＋｛Δδ｝l。
（4）由｛Δδ｝l求各单元应变增量｛Δε｝l和应力增量｛Δσ｝l，则｛ε｝l＝｛ε｝l－1＋｛Δε｝l，｛σ｝l＝｛σ｝l－1＋｛Δσ｝l。
（5）对各级荷载重复上述步骤，可以求出最后的解。
3.1.2.2 中点增量法
基本增量法由于用初始应力求D，每级荷载都有一定的误差。事实上，对于某一级荷载，应力从初始状态变化的终了状态，弹性常数也是变化的。如果用该级荷载下的平均应力所对应的D进行计算将会使得结果有所改善，这就是中点增量法。为了求平均应力，要做一次试算，按基本增量法先算一次，得出的应力与初始应力平均，就得到该级荷载的平均应力，再求D，重新计算一次。第l级荷载的计算步骤如下：
（1）～（4）同基本增量法。
（5）
（6）由平均应力
再形成
（7）解方程组
求出位移增量｛Δδ｝l，相应的位移总量为｛δ｝l＝｛δ｝l－1＋｛Δδ｝l。
（8）由｛Δδ｝l求应变增量｛Δε｝l和应力增量｛Δσ｝l，进而求应力和应变全量。
3.1.2.3 增量迭代法
对每一级荷载增量，用迭代法多次计算，使其收敛于真实解，再加下一级荷载。但增量迭代法计算的时间较长，从理论上讲，解的结果最精确。但对于土体而言，由于本构关系的复杂性，迭代未必都是收敛的，因此，工程中广泛应用的还是中点增量法。
图3.1 程序流程图
3.1.3 模拟施工步骤的实现过程
（1）根据试验资料，计算坝体的初始参数，包括初始弹性模量Ei和初始泊松比μi，然后根据坝体材料分区分别赋予坝体各部位特定的初始材料参数。
（2）利用给定的材料参数、初始条件和边界条件进行填筑坝体第一步的计算。
（3）通过上一步的计算，确定每一单元的应力状态。根据第一层单元的应力状态，计算其切线模量Et和切线体积模量Bt。
（4）修改填筑坝体第一步的刚度矩阵，进行新填筑部分和已有坝体共同作用的非线性有限元计算。
（5）重复（3）～（5）步，进行坝体施工过程模拟，直至坝体竣工并蓄水完毕。程序流程图见图3.1。