基于邓肯-张模型的堆石坝有限元分析

2016-03-23汪天飞武汉大学水资源与水电工程科学国家重点实验室武汉430072

中国农村水利水电 2016年10期

汪天飞，徐青(武汉大学水资源与水电工程科学国家重点实验室，武汉 430072)

1 研究背景

随着中国水利水电建设事业的发展，堆石坝因具有对不同地形、地质和气候条件适应性好，造价低，施工速度快，结构简单，抗震性能强等优点而被广泛应用。近年，随着我国经济的不断发展，面板堆石坝得到了快速发展。随着堆石坝坝高的增加，坝体应力和变形分析已成为坝体设计中不可缺少的一部分[1]。因此，如何合理计算坝体变形是堆石坝设计中的一项重要工作。堆石料性质十分复杂，尤其是其强度特性，自1963年剑桥模型被提出了以后，基于大量的试验成果和理论分析，国内外学者提出了不同的堆石料本构模型。

人们通过不断研究与实践，认为非线性弹性模型和弹塑性模型可以较好地反映堆石料的变形特性。其中常用的非线性弹性模型以邓肯-张模型为代表，弹塑性模型有黄文熙提出的弹塑性模型、沈珠江提出的双屈服面弹塑性模型等[2,3]。

ANSYS软件是一款大型通用有限元软件，目前行业内普遍采用ANSYS软件进行有限元计算，但它不包含被岩土工程界所广泛采用的邓肯-张本构模型。而实际研究表明，邓肯-张模型意义明确，计算参数简单且易于通过试验获得，同时该模型能够反映堆石体的主要应力变形特性，计算结果精度能够满足工程需要。近年来，不少研究者将邓肯-张模型嵌入ANSYS、FLAC_3D和ABAQUS中，用于计算堆石坝的受力变形，但未详细说明如何在每一步加载前考虑已填筑坝体的初始应力。在堆石坝应力应变计算过程中，无论是新填筑坝体还是浇筑完成的坝体，坝体单元都具有一定的初始应力。每一荷载步的初始应力状态均为前一荷载步终了时的应力。因此，本文基于ANSYS研究平台，对ANSYS进行二次开发，通过参数化设计语言APDL语言编写程序，完成了邓肯-张本构模型的添加[4]。本文在堆石坝应力变形模拟过程中，输出每一荷载步加载完成后的应力，作为下一荷载步的初始应力，使大坝的填筑和蓄水过程更接近工程实际情况。由于考虑坝体单元初始应力时，会产生相应的附加位移，本文采用约束荷载平衡法消除坝体初始应力引起的附加位移对结果造成的影响，并将二次开发程序用于堆石坝应力变形计算。计算结果表明程序能较好地应用于工程实践。

2 土石坝非线性静力分析方法

2.1 邓肯-张E-B模型[5]

1970年邓肯等人在康纳(Kondner)提出的用双曲线拟合三轴试验曲线(σ1-σ3)-ε1的假定基础上，提出了双曲线模E-μ型，后来又提出了采用切线体积模量代替切线泊松比的E-B模型。邓肯-张E-B模型的理论及公式推导过程，本文不再赘述。

切线弹性模型：

切线体积模量：

泊松比：

回弹模量：

土体初始应力状态[6,7]：

σ1=γh

σ3=K0γh

K0=0.95-sinφ

式中：h为单元形心到土体表面的深度。

该模型同时还考虑粗粒料内摩擦角φ随围压σ3的变化，内摩擦角为：

根据上面的计算公式可知，要采用邓肯-张模型，只需要确定φ、C、K、n、Rf、Kb、m、Kur及nur9个参数，这些参数都可以根据常规三轴实验得到[2]。

2.2 非线性有限元分析方法

非线性有限元分析方法[8]主要有增量法和迭代法2种，尤以增量法更为普遍。增量法思路为：将荷载划分为若干级荷载增量，每次施加一个荷载增量，求解每级荷载增量下的应力增量和应变增量，最后求和得到全荷载作用下的总位移、总应力和总应变。增量法假定材料是线弹性，即方程是线性的，劲度矩阵是常数，而在不同的荷载增量中，劲度矩阵也可以是不相同的。这种方法是以很多条直线来近似逼近曲线的，当荷载划分较小时，能得到近似于真实情况的解。增量法分为基本增量法和中点增量法两种。考虑到基本增量法采用初始应力所对应的弹性常数，使得最后求解的结果与实际曲线有相当大的偏离。本文采用中点增量法进行非线性问题的求解。

中点增量法计算步骤为：先按基本增量法计算一次，将计算得到的应力与初始应力平均，得到该级荷载下的平均应力(中点应力)，再求一次刚度矩阵[D]，根据新的刚度矩阵计算求得应力、应变和位移增量。对第级荷载，中点增量法计算步骤如下。

(1)用前级终了时的应力，即本级的初始应力{σ}i-1，确定刚度矩阵[D]i，形成劲度矩阵[K]i。

(2)解线性方程组[K]i{Δδ}i={ΔR}i,得到位移增量{Δδ}i，则相应位移总量为{δ}i={δ}i-1+{Δδ}i。

(3)由{Δδ}i求解各单元应变增量{Δε}i和应力增量{Δσ}i，相应的应变和应力总量为{ε}i={ε}i-1+{Δε}i，{σ}i={σ}i-1+{Δσ}i。

(7)由{Δδ}i求解应变增量{Δε}i和应力增量{Δσ}i，进而得到应力和应变总量。

对各级荷载重复上述步骤就可以得到总荷载作用下结构产生的总应力、总应变和总位移。

2.3 APDL语言开发邓肯-张模型

邓肯-张模型是一种建立在增量广义胡克定律基础上的非线性弹性模型，主要通过不断改变其切线模量来实现其非线性过程。本文对ANSYS进行了二次开发，利用参数化语言APDL[4]编写程序，建立邓肯-张模型，模拟工程实际施工工序，逐步填筑坝体，施加荷载，最后计算得到堆石坝在施工期、竣工期的应力应变大小。该程序基于参数化设计，使用方便，对其他堆石坝工程应力变形分析具有实际意义。

2.3.1施工过程模拟[9-11]

利用ANSYS进行土石坝施工模拟时，先整体建模，根据填筑过程对单元进行划分，定义单元的材料参数。利用APDL语言编制土石坝施工模拟的计算程序。基本思路如下。

(1)根据施工过程和坝体的剖面图，分层分区建模，建模完成后进行单元网格划分，然后根据土体单元的初始应力状态，计算出初始材料参数，给每个单元定义一个材料号，把相应的材料赋给每一个单元。

(2)创建邓肯-张E-B模型宏命令。

(3)杀死所有坝体单元，对坝基部分进行计算，考虑坝基初始地应力的影响。

(4)激活坝体第一层单元，读入单元初始应力，根据初始材料参数进行计算，计算完成后提取每个单元的第1、第3主应力，根据邓肯-张E-B模型宏命令，采用中点增量法计算出激活单元相应的弹性模量和泊松比，修改各单元的材料属性，修改单元刚度矩阵，重新进行计算，此时计算得到的位移和应力作为该荷载步加载完成后的计算结果，同时在计算过程中输出单元的应力作为下一荷载步的初始应力，坝体的材料参数作为下一个荷载步计算的初始材料参数。

(5)依次加载，方法同步骤(4)，直到坝体填筑完成。

(6)根据计算结果进行相应的后处理，编制相应的程序，得到坝体施工完成后的应力、应变结果。施工过程流程图见图1。

图1 施工过程流程图

2.3.2运行过程模拟

运行过程的模拟是随着水位的抬高，将水荷载逐步施加在坝体上的，在进行土石坝运行过程仿真模拟时，应根据实际水位变化情况逐步施加水荷载。运行过程中坝体应力应变的计算方法与施工期一样，只是施加的荷载变为水荷载增量。运行过程流程图见图2。由于库水压力是随着库水水位的升高而升高的，在下一荷载步蓄水施加水压力时，上一荷载步已蓄水部位的水压力也要增加，具体见图3。

图2 运行过程流程图

图3 水荷载逐级施加示意图

2.3.3APDL语言相关说明

(1)坝体填筑过程模拟。为了实现对坝体填筑过程的模拟，ANSYS软件提供了“生死单元”功能，由EALIVE(激活)和EKILL(杀死)来实现。

(2)单元第1、第3主应力提取。首先利用ETABLE命令创建单元表，然后用“*get”命令提取单元主应力。例：

ETABLE,EtabS1,S,1(创建单元表)

*get,ArrS1(num),elem,Num,etab,EtabS1(提取单元第1主应力)

(3)单元初始应力的读入、输出。在参数化语言APDL中，最开始用于输出、读入初始应力文件的命令为iswrite和isfile，随着软件的更新，逐渐被inistate,write和inistate,read命令所代替。例：

isfile,read,dam,ist,,1(单元初始应力读入)

iswrite,on(单元初始应力输出)

2.3.4约束荷载平衡法

运用ANSYS进行数值模拟时，单元初始应力被当成一种“荷载”施加在每一个单元上，在仅有初始应力荷载时，读入初始应力文件后，将产生一个大小与原荷载相同，方向与原荷载相反的附加位移。为了避免这一附加位移对计算结果造成影响，本文采用了约束荷载平衡法。即施加相同的边界条件和与原荷载相同的荷载，使坝体产生的应力场与初始应力相同，而初始位移场为零(位移一般很小，可以忽略不计)[12]。

3 计算实例

本文算例选取均质坝和非均质坝分布进行计算。

3.1 均质坝

算例为一均质堆石坝，坝高100 m,坝顶宽度8 m，上下游坡角为30°，堆石坝密度为2.0 t/m3，其泊松比大致变化范围是0.31～0.34，根据工程经验，算例中泊松比取0.3。邓肯-张E-B模型所需参数见表1。

表1 邓肯-张E-B模型材料参数

程序对该均质坝坝体分10级均匀加载，即每层填筑高度为10 m。计算结果见图4～7。

图4 最小主应力等值线图(单位：Pa)

图5 最大主应力等值线图(单位：Pa)

图6 水平方向位移等值线图(单位：m)

图7 垂直沉降位移等值线图(单位：m)

计算结果分析：由图4、图5可知，在自重作用下，该均质坝大小主应力呈对称分布，大主应力最大值为1.5 MPa，发生在坝底底部中央，随着坝体高度的增加逐渐减小，小主应力最大值为0.62 MPa，分布规律与大主应力相同。由图6、图7可知，在自重作用下，坝体发生垂直沉降和水平位移，该均质坝为对称结构，坝体沉降和水平位移也呈对称分布。坝体沉降最大值发生在坝体中部偏上1/2～1/3坝高处，最大沉降值为72 cm，大约占坝高的0.72%(满足堆石坝坝体总沉降量为坝高1%以内要求)，方向竖直向下；水平位移则对称发生在坝体的中下部，上下游坝坡附近，水平位移最大值为18 cm。根据分析可知，坝体应力变形结果符合施工过程逐层加载的计算规律。

3.2 非均质坝

所选模型为某工程大坝最大断面剖图，最大坝高为88 m，坝顶宽度为8 m，坝顶高程为515 m，坝底高程为427 m，材料参数见表2，剖面图如图8所示，竣工期和运行期的计算结果见图9～16。计算结果分析如下。

竣工期计算结果见图9～图12。由图9～图10可知，坝体大、小主应力的最大值都发生在底部，大主应力最大值为1.6 MPa，小主应力最大值约为0.46 MPa。由图11可知，坝体水平方向位移为上游侧倾向上游，坝体下游侧倾向下游，由于坝体下游侧次堆石料的工程特性比主堆石料的工程特性差，坝体下游侧的位移明显比上游侧大，且最大位移发生在次堆石区部位，高程约为次堆石区的1/3处，水平方向位移最大值约为19 cm。由于下游侧坝体向下游发生变形，坝顶以下约1/4坝高的坝体也向下游变形。由图12可看出，与均质坝分布规律基本相同，坝体最大沉降为88 cm，最大值出现在次堆石区部位高程约为1/2～2/3坝高处。

表2 坝体材料力学参数

图8 最大断面剖面图

图9 竣工期最小主应力等值线图(单位：Pa)

图10 竣工期最大主应力等值线图(单位：Pa)

图11 竣工期水平方向位移等值线图(单位：m)

图12 竣工期垂直沉降位移等值线图(单位：m)

图13 运行期最小主应力等值线图(单位：Pa)

图14 运行期最大主应力等值线图(单位：Pa)

图15 运行期水平方向位移等值线图(单位：m)

图16 运行期垂直沉降位移等值线图(单位：m)

运行期计算结果见图13～图16。由图13～图14可知，大、小主应力的最大值都发生在底部，大主应力最大值为1.7 MPa，小主应力最大值约为0.56 MPa。运行期由于水压力作用，坝体整体向下游移动，上游侧坝体在水压力作用下，竣工期产生的向上游的位移被抵消，并开始向下游变形。下游侧坝体的水平位移增加，水平位移最大值约为25 cm，出现在坝体下游部位，高程约为次堆石区的1/3处。由图16可知，在水压力作用下，坝体沉降与竣工期规律基本相同，沉降最大值为88 cm，基本与竣工期没有差别，但是坝体上游侧的位移明显比竣工期大。该工程中坝体竣工期、运行期的最大位移与应力值见表3。