邓肯非线性模型的加卸荷准则对计算结果的影响分析[1]
米占宽 李国英
（南京水利科学研究院）
摘 要：Duncan非线性模型已广泛应用于土石坝的数值计算。加荷函数的选择是Duncan非线性模型的一个难题。由于卸荷模量和加荷模量相差悬殊（最大可达到2个数量级以上），计算结果随选用的加卸荷准有较大关系，从而使计算结果带有一定的任意性。目前一般采用Duncan等人于1980年提出的加荷函数，该函数在低应力和高围压条件下通常会低估体积模量，从而低估最小主应力和土体刚度，为此1984年Duncan等人又对最小体积模量进行了修正，在土石坝数值分析中通常认为E-B模型较之E-ν模型更好一些。本文以马吉水电站混凝土面板堆石坝为例，比较分析了E-B模型和E-ν模型的不同加卸荷准则对计算结果的影响，据此提出了相应的结论和建议。
关键词：Duncan非线性模型；加卸荷准则
1 引言
由于Duncan等人的双曲线模型可以反映土变形的非线性，并且能通过弹性常数的调整在一定程度上近似地反映塑性变形；同时由于它是建立在广义虎克定律的弹性理论基础上，易于为工程界所接受；加之参数不多、物理意义明确，且只需通过常规三轴试验即可确定；再者适用于堆石料、黏性土和砂土等大部分土体，目前已成为土工数值分析中最为普及的本构模型之一。
加荷函数的选择是Duncan非线性模型的一个难题。由于卸荷模量和加荷模量相差悬殊（最大可达到2个数量级以上），容易引起两个严重问题：一是计算结果随选用的卸荷准则不同而不同，从而带有一定的任意性；二是当相邻两个单元模量相差悬殊时，刚硬的单元往往由于应力集中而发生受拉或受剪破坏。Duncan等人于1980年和1984年分别提出了一个加卸荷函数，尤其是1984年对最小体积模量又做了一定限制，在土石坝数值分析中通常认为E-B模型较之E-ν模型更好一些。本文以马吉水电站混凝土面板堆石料坝为例，比较分析了E-B模型和E-ν模型的不同加卸荷准则对坝体应力变形尤其是面板应力变形的影响，在此基础上提出了相应的结论和建议。
2 Duncan非线性模型及其加卸荷准则简介
2.1 Duncan E-ν模型
Duncan-Chang（1970）建议以双曲线拟合三轴试验的偏应力（σ1-σ3）与轴向应变εa的关系以及侧向应变εr与轴向应变εa的关系，在此基础上得出切线杨氏模量Et和切线泊松比νt的按式（1）计算：
式中：c、φ、Rf、K、n、G、F、D分别为8个计算常数；Pa为大气压力。
对粗粒料来说，c=0，同时考虑到粗粒料的摩尔库仑包线往往不是直线，采用下列公式计算内摩擦角φ：
卸荷再加荷条件下的杨氏模量则按下式计算：
2.2 Duncan E-B模型
Duncan等人认为，E-ν模型关于εr与εa间的双曲线假设与实际情况相差较多，同时使用切线泊松比νt计算也有一些不便之处。1980年Duncan等人提出了E-B模型[2]，其中Et的确定与式（1）相同，只是以切线体积模量Bt代替切线泊松比νt：
Bt与νt有如下关系：
2.3 Duncan非线性模型应力应变关系
求得Et和νt或Bt后，平面问题的增量型应力应变关系表达为：
其中劲度矩阵系数为：
或者
2.4 加卸荷准则
1980年Duncan等人采用E-B模型编制的用于土石坝计算的FEADAM程序中，采用以应力水平Sl是否大于历史上的最大值Slmax作为卸荷准则。1984年版本中又做了如下修正：
（1）用应力函数
取代Sl作为卸荷准则。
（2）Ss大于历史上的最大值Ssmax时用加荷模量Et，Ss小于历史上的0.75Ssmax时用卸荷模量Eur，介于两者之间时则用Et和Eur内插。
（3）当σ3减小时，上述公式中的σ3一律用历史上的最大值
。
近年来国内很多科研院所和高等院校均编制了有限元程序用于土石坝数值分析，一般均采用上述加卸荷准则。在编制程序时，为构建E-ν模型和E-B模型统一的矩阵格式，一般采用式（11）的应力应变关系式，即在采用E-B模型进行数值分析时，在求得体积模量Bt后再通过式（9）进行转换求得切线泊松比。此时通常不对Bt进行修正，导致各家单位的计算结果有较大差别，实际上1984年Duncan等人分析认为在非常低的应力水平和高应力情况下，通常会低估体积模量，从而使得计算得到的泊松比过小，又对最小体积模量限制如下：
即限制泊松比的最小值为：
式（14）为单调递减函数，当φ=2.3°时，ν=0.49；当φ=90°时，ν=0.0；一般堆石体的φ不会超过50°，此时ν=0.19，为保证刚度矩阵为正定，需对ν限制为：0＜ν≤0.49。由此可以看出，当采用E-B模型计算时，泊松比最小值限制为0.19，而采用E-ν模型计算时，泊松比最小值限制为0.01，由此必然导致E-B模型计算得到的沉降偏小，水平位移偏大。
3 马吉水电站混凝土面板坝三维有限元应力变形数值分析
本节分别采用邓肯E-B模型和E-ν模型对马吉水电站混凝土面板坝进行了三维有限元应力变形数值模拟，比较分析了两种模型计算的坝体和面板应力变形的差异，在此基础上分析了产生差异的原因。其中面板和垫层间接触面模型均采用Desai薄层单元双曲线接触面模型。
3.1 马吉水电站混凝土面板坝工程概况
马吉水电站水库正常蓄水位1570.0m，死水位1500.00m，总库容46.56亿m3，混凝土面板堆石坝方案拦河坝最大坝高为275.5m，工程为一等大（1）型工程。
（1）坝料分区及断面设计。面板堆石坝坝顶高程1575.50m，河床段趾板底高程1300.00m，最大坝高为277.5m，坝顶长683.00m，坝顶宽14.00m；坝顶上游侧设有L形防浪墙，防浪墙顶高程1576.70m，墙底高程1572.00m。拦河坝上游坝面坡比为1∶1.5；下游坝坡1500.00m以上采用1∶1.6，下部采用1∶1.5，下游坝坡每50m高度设有宽5m马道平台，下游综合坡比1∶1.6（见图1）。
图1 河床段坝体断面图（单位：m）
（2）填筑及蓄水过程。坝体填筑施工顺序见图2。概述如下：
1）大坝分6期填筑。
2）面板浇筑：分为4期（高程分别为1390.00m、1450.00m、1505.00m、1572.50m），其中第1期面板在大坝第2期填筑时进行；第2期面板在大坝第3期填筑时进行；第3期面板在大坝第6期填筑时进行；第4期面板在大坝填筑全部完成后进行。
3）浇筑面板时的填筑超高：第1期面板，填筑超高50m；第2期面板，填筑超高55m；第3期面板，填筑超高67.5m。
4）蓄水过程：分为三个阶段，第1次蓄水在第4期坝体填筑时进行（此时，二期面板浇筑完成），水位由1335.79m上升至1420.50m；第2次蓄水在第6期坝体填筑时进行（此时，三期面板浇筑完成），水位上升至1500.00m；大坝竣工后，水库水位升至正常蓄水位1570.00m。
图2 马吉面板坝施工顺序图（单位：m）
3.2 单元离散及施工、蓄水模拟方法
图3和图4分别为平面和空间三维网格图，其中x向表示坝轴向（自左岸指向右岸为正），y向表示顺河向（自上游指向下游为正），z向表示垂直方向。面板分缝间距，河床部位为16m，两岸为8m，在进行三维有限元网格剖分时为适应边界条件的变化，沿纵向增加了部分断面，但增加的断面不进行分缝，沿坝轴向共截取了66个断面，共剖分单元数15403，节点数20013，三维实体单元一般采用8结点六面体等参单元，为适应边界条件以及坝料分区的变化，部分采用三棱体和四面体作为退化的六面体单元处理，单元编号根据坝体施工顺序进行。接触面厚度为0.1m。
图3 标准断面网格剖分图
根据大坝填筑和蓄水过程，有限元计算时模拟的填筑、蓄水顺序为：Ⅰ期坝体填筑（大坝填筑至高程1395.00m）→Ⅱ期坝体填筑（大坝填筑至高程1460.00m）→一期面板浇筑（高程1390.00m）→Ⅲ期坝体填筑（大坝填筑至高程1505.00m）→二期面板浇筑（高程1450.00m）→初次蓄水：水位上升至1420.50m→Ⅳ期坝体填筑（大坝填筑至高程1525.00m）→Ⅴ期坝体填筑（大坝填筑至高程1525.00m）→三期面板浇筑（高程1505.00m）→二次蓄水：水位上升至1500.00m→Ⅵ期坝体填筑至坝顶（高程1575.50m）→四期面板浇筑→水位上升至1570.00m。共分为50级进行模拟，其中坝体填筑分30级模拟计算，蓄水过程分20级模拟。
图4 三维网格剖分图
3.3 计算参数
需要说明的是，马吉水电站混凝土面板坝作为“300m级高面板堆石坝安全性及关键技术研究”依托工程之一，考虑到试验得到的坝料模量较大，根据试验成果对E-B计算参数进行了调整，调整方法如下：垫层区料、过渡区料、增模特碾区的K值调整为试验参数的0.80倍，Kb值调整为试验参数的0.70倍；上游堆石区、下游堆石区的K值调整为试验参数的0.85倍，Kb值调整为试验参数的0.75倍。为保证调整后E-ν模型的泊松比变化同E-B模型一致，本文将E-ν模型中有关切线泊松比参数的调整方法如下：将G、F分别调整为试验值的0.80倍（垫层区料、过渡区料、增模特碾区）和0.75倍（上游堆石区、下游堆石区）。最后采用E-B模型、E-ν模型参数见表1。
表1 E-B模型、E-ν模型参数
混凝土面板采用C35混凝土，按线弹性模型考虑，其弹性模量、泊松比和密度分别为E=31.5GPa，ν=0.167，ρ=2.40g/cm3。
3.4 计算结果对比
（1）坝体变形。E-B模型和E-ν模型计算得到的竣工期和蓄水期河床典型断面的应力变形计算结果（见表2），计算结果对比（见图5），限于篇幅仅给出蓄水期（系指蓄水至1570.00m高程）的图形。从上述表2、图5可以看出：①由于E-B模型限制的最小泊松比较之E-ν模型大，故前者计算的顺河向位移较大，而沉降较小。由此也导致后者计算得到的计算大主应力较大而小主应力较小；②在坝体上部的低应力区，两种模型计算值差别较小，但在坝体下部的高应力区，两种模型计算值差别较大，这主要是由于在高应力区域摩擦角较小，此时E-B模型所限制的泊松比较大所致；③坝体应力变形的极值变幅在10%之内，两种模型计算偏差较小，在数值分析中是可以接受的。
表2 坝料本构模型对坝体应力变形的影响
图5 蓄水期河床断面0+360变形分布图（单位：cm）
（2）面板应力变形。E-B模型和E-ν模型计算得到的竣工期和蓄水期面板应力变形特征值（见表3），计算结果对比（见图6、图7），限于篇幅仅给出蓄水期（系指蓄水至高程1570.00m）的图形。计算结果表明：①由于E-B模型限制的最小泊松比较之E-ν模型大，前者计算得到的面板坝轴向位移和顺河向位移略大，虽然前者计算得到的堆石体最大沉降较小，但前者坝体底部的沉降反而比后者大，由于面板最大挠度发生在面板的中下部，故E-B模型计算得到的面板挠度反而略大；②由于两种模型计算得到的面板坝轴向变形和挠度相差较小，由此计算得到的面板坝轴向拉压应力和顺坡向压应力差别也相差不大，差别较大的是面板的顺坡向拉应力，这主要是由于面板顺坡向拉应力发生在面板的底部，而该部位两种模型计算得到的变形差别较大所致。
表3 坝料本构模型对面板应力变形的影响表
图6 蓄水期面板变形分布图（单位：cm）
图7 蓄水期面板应力分布图（单位：MPa）
4 结论
本文以马吉水电站混凝土面板堆石坝为例，对其进行了三维有限元静力应力变形分析。面板与垫层间接触面模型分别采用Desai双曲线接触面模型，堆石料本构模型分别采用E-B模型和E-ν模型，比较分析了两种模型计算得到的坝体和面板应力变形的差异，并分析了差异的原因。结果表明：
（1）E-ν模型和E-B模型计算结果差别主要是卸荷处理不同造成的，前者假定卸荷时不变νt不变，后者假定Bt不变，但对Bt的最小值进行一定的限制，以使νt不致过小。
（2）E-ν模型和E-B模型计算得到的坝体应力变形、面板坝轴向拉压应力、面板顺坡向压应力差别约10%，在土石坝数值分析中，这种差别是可以接受的。差别较大的是发生在面板底部的顺坡向拉应力，这主要是由于两种模型在高应力区域计算得到的坝体变形差别较大所致。实测资料表明，面板底部不会发生顺坡向拉应力，因此两种模型的计算结果均是不合理的。
（3）《碾压式土石坝设计规范》（DL/T 5395—2007）条文说明10.4.5中，认为E-ν模型具有试验资料拟合不理想等缺点，应用于面板坝应力变形分析上是不合适的，同时认为E-B模型“参数确定简单，且在参数确定方面积累了比较成熟的经验，使用简便”，故推荐采用E-B模型。但实际上，E-ν模型和E-B模型均不能反映材料的剪胀特性，由于在低围压下，坝料通常会发生明显的剪胀，此时参数的确定两者具有同样的困难，加之两种模型计算结果差别较小，因此，不存在孰优孰劣的问题。由于两种模型计算得到的顺坡向拉应力均不合理，最好还是采用能合理反映堆石料剪胀特性的弹塑性模型进行数值分析。
参考文献
[1] 沈珠江.沈珠江土力学论文集.北京：清华大学出版社，2005.
[2] Duncan，J.M.，Wong，K.S.and Ozawa，Y.FEADAM：A computer program for finite element analysis of dams.University of California，Berkeley（Rep.No.UCB/GT/80-82）.
[3] Duncan J M，et al.FEADAM-84，A Computer Program for Finite Element Analysis of Dams[R].Report No.UCB/GT/University of California，Berkeley 1984.
[4] 朱百里，沈珠江.计算土力学.上海：上海科学技术出版社，1990.
[1]基金项目：国家重点研发计划项目课题经费资助（项目批准号：2017YFC0404803）；中央级公益性科研院所基本科研业务费专项资金（Y317005）。