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）根据试验资料，计算坝体的初始参数，包括初始弹性模量Ei和初始泊松比μi，然后根据坝体材料分区分别赋予坝体各部位特定的初始材料参数。

（2）利用给定的材料参数、初始条件和边界条件进行填筑坝体第一步的计算。

（3）通过上一步的计算，确定每一单元的应力状态。根据第一层单元的应力状态，计算其切线模量Et和切线体积模量Bt。

（4）修改填筑坝体第一步的刚度矩阵，进行新填筑部分和已有坝体共同作用的非线性有限元计算。

（5）重复（3）～（5）步，进行坝体施工过程模拟，直至坝体竣工并蓄水完毕。程序流程图见图3.1。