土石坝心墙接触面单元参数的标定与数值试验方法

徐 畅，詹程远,2，何 吉

(1. 武汉大学水利水电学院，湖北 武汉 430072； 2. 长江科学院，湖北 武汉 430010)

数值模拟方法已经成为工程设计的重要手段而被广泛采用，但是该方法在使用时经常由于材料参数的不合理选取而引起模拟结果的显著偏差，甚至造成严重错误[1]，因此如何获取合理的材料参数一直是数值模拟的重点和难点。

对于沥青混凝土心墙坝而言，坝体主要由堆石料坝壳、沥青混凝土心墙、过渡层(位于心墙两侧)等部分组成[2]。这些结构通常基于有限单元法，采用非线性弹性模型(邓肯-张模型等)[3]或者弹塑性模型(K-G模型等)[4]等本构模型进行模拟，相应的材料参数可以通过成熟的试验方法获得。但是由于心墙与过渡层的力学性能存在显著差异，因此现实中两者间存在剪切、错动等相对变形，并且伴随局部侵入等复杂现象[5]。这些复杂力学行为的合理模拟对于正确判断心墙的应力状态，避免防渗体的破坏至关重要。为了模拟心墙和过渡层之间复杂的力学行为，两者间通常需要设置接触面单元[6]。常用的接触面单元分为无厚度单元和有厚度单元两种，前者以Goodman单元[7]为代表，后者以Desai单元[8]为代表。Goodman单元在受压时接触面两侧存在相互嵌入的问题。为了减少嵌入的程度，通常需要显著增大接触面的法向刚度，但是这样极易造成迭代计算时误差的明显积累[9-10]。土石坝的数值模拟无论在施工期、蓄水期，还是运行期均涉及大量迭代计算，因此更适合采用Desai单元。

本文建立了一套标定土石坝心墙接触面单元参数的数值试验方法。该方法的主要步骤包括：

1)采用邓肯-张非线性弹性模型和Desai接触面单元，建立土石坝心墙接触面直剪试验的数值模型；

3)基于数值试验，标定增厚的接触面单元参数，本文依托茅坪溪土石坝的物理试验成果，采用上述方法，标定了心墙接触面单元在不同厚度时的参数，研究了各参数随单元厚度的变化规律，并将标定所得参数成功应用于某土石坝工程的三维数值模拟中。

1 接触面物理试验

茅坪溪土石坝是三峡水利枢纽的重要组成部分。为了研究该坝沥青混凝土心墙与砂砾石过渡层之间接触面的力学性能，张治军[11]等人采用大型叠环单剪仪进行了试验。试验设备如图1所示，图1中下部为方形剪切盒，尺寸600 mm×600 mm×300 mm(长×宽×高)，盒内装有砂砾石材料；上部为堆叠放置的10层方环，每环尺寸600 mm×600 mm×30 mm，底部两层方环内装有沥青混凝土材料，其余8层方环内装有砂砾石材料。试验中保持法向应力1.5 MPa不变，逐步增大水平推力直至试件破坏。试验结果显示，破坏面沿上接触面或下接触面随机分布，接触面剪切位移与剪切应力的关系曲线如图2所示。

图1 直剪物理试验示意图

图2 接触面剪切位移与剪切应力的关系曲线

2 数值试验及参数标定

2.1 本构模型

2.1.1 砂砾石及沥青混凝土本构模型

根据专著《长江三峡水利枢纽建筑物设计及施工技术》[12]，茅坪溪土石坝在大部分研究中沥青混凝土心墙、砂砾石过渡料以及坝体各种堆石料均推荐采用邓肯-张E-μ模型，并且提供了具体参数，如表1所示。本文在数值试验中继承了上述研究成果。

表1 心墙及过渡料邓肯-张E-μ模型参数

邓肯-张E-μ模型为非线性弹性模型，相应的加载切线变形模量Et、卸载切线变形模量Eur、切线泊松比μt分别为[12]：

(1)

(2)

(3)

式中：c、φ为抗剪强度参数；Rf为破坏比；Pa为大气压力；σ1、σ3为大、小主应力，以压为正；K、n、F、G、D为无量纲参数。

2.1.2 接触面本构模型

本文为了模拟沥青混凝土心墙与砂砾石过渡层之间的复杂力学行为，在两者间设置了Desai接触面单元。该单元的应力应变关系服从[12]

{dε}=[C]{dσ}

(4)

(5)

式中：[C]为接触面的柔度矩阵；Et、μt为接触面的切线弹性模量和切线泊松比，假定与砂砾石过渡料相同(式(1)～式(3))；Gt为接触面的切线剪切模量，其表达式为

(6)

τf=σntanφs+cs

(7)

式中：τf为接触面的抗剪强度；σn为接触面的法向正应力；φs、cs为接触面的摩擦角和黏聚力；Rfs为接触面的破坏比；Gs、ns为无量纲参数。

2.2 数值试验模型

根据图1中物理试验的试件尺寸和材料分布，建立三维数值模型(见图3)，同时施加约束和荷载如下：物理试验中方形剪切盒的约束作用，以法向约束的方式施加于数值模型下部砂砾石材料的所有外表面；物理试验的法向应力和千斤顶反力，分别以面力的形式沿铅直方向和水平方向施加于数值模型的顶面。

图3 数值试验模型图

2.3 标准厚度接触面参数反演

假定一组接触面参数(包括cs、φs、Rfs、Gs、ns)，基于数值模型(见图3)，采用有限元法模拟直剪试验。计算接触面边界之间(即图3中ab点之间)的平均水平错动量，并且绘制其与剪切应力(即单位面积的平均水平推力)的关系曲线(见图2)。采用参数反演方法，反复调整接触面的各项参数，直至该曲线与物理试验的结果相近为止(见图2)，此时接触面的参数即为反演所得参数，如表2中第一行所示。

表2 不同厚度的接触面参数

2.4 增厚接触面参数标定

为了考虑不同厚度对接触面参数的影响，将接触面厚度依次增加至20、50、150 mm。基于不同接触面厚度的数值模型，分别模拟直剪试验。计算试件顶面与底面之间(即图3中cd点之间)的平均水平错动量，并且绘制其与剪切应力(即单位面积的平均水平推力)的关系曲线(见图4)。以接触面厚度为1 mm的曲线为基准，反复调整其他厚度接触面的各项参数，直至所有厚度接触面对应的曲线与基准曲线相近为止(见图4)，此时即为标定所得参数，如表2所示。

图4 试件剪切位移与剪切应力的关系曲线

由表2可知，各参数与接触面厚度的相关性中，Gs最大，cs、φs及Rfs较小，ns基本无关；cs、φs及Gs与厚度正相关，Rfs为负相关。由于增厚以后的接触面单元不仅反应了接触面的力学性能，也反应了部分过渡层材料的性能，因此随着厚度增加，过渡层的影响逐渐增大，考虑到过渡层的力学性能强于接触面，造成接触面单元的力学性能逐渐增强。表2中各参数随厚度的变化规律符合上述原理。

3 工程应用

由于无法获取茅坪溪土石坝详尽的设计资料，不能为其建立三维数值模型，本文选取某沥青混凝土心墙坝作为工程应用。该坝坝高64.5 m，坝轴线长586 m，心墙厚0.4 m，过渡层厚2 m。茅坪溪与该坝的材料分区并非完全相同，为了避免不同分区的干扰，专注于接触面的研究，因此假定坝体仅采用一种堆石料，同时假定坝基仅含有一种基岩材料，材料分区如图5所示。

图5 坝体材料分布图

图6 三维有限元网格图

表3 大坝邓肯-张E-μ模型材料参数

表4 大坝线弹性模型材料参数

基于三维有限元网格(见图6)，采用有限单元法模拟了坝体的分层施工过程，总计16个计算步。施工结束时坝体的沉降分布如图7所示，由图7可知该坝施工期最大沉降量为591.55 mm，位于坝体中部。由于最大沉降量小于坝高1%，所以符合土石坝设计规范[13]的要求。心墙上游接触面与过渡层之间的错动量如图8所示，由图8可知心墙相对过渡层整体向下错动，即心墙的沉降量大于过渡层，坝体上部错动量较大、下部较小，最大错动量在1.8 mm以内。为了比较错误的接触面参数对错动量的影响，此处选取表2中1 mm厚的参数进行计算，得到心墙与过渡层间的错动量如图8所示。由图8可知，采用1 mm厚的参数时会明显夸大层间错动量。这是因为Gs与厚度正相关，厚度越小时接触面越容易发生变形。因此增厚以后的接触面单元不能沿用增厚之前的参数，工程中需要根据接触面单元的实际厚度标定参数，否则会对结果产生显著影响。

图7 施工期最大坝高断面沉降分布图

图8 心墙上游接触面错动量随高程分布图

4 结 语

本文针对土石坝沥青混凝土心墙接触面单元的参数取值问题，建立了一套数值试验方法，研究了接触面参数随单元厚度的变化规律，并且进行了工程应用。主要结论如下。

1)为了模拟沥青混凝土心墙与砂砾石过渡层之间复杂的力学行为(剪切、错动、局部嵌入等)，推荐在两者间设置Desai接触面单元。该单元是一种常用的有厚度接触面单元，使用方法简便。

2)在土石坝的数值建模中，特别是三维建模，为了降低建模难度、减少单元总数，通常需要选择较厚的接触面单元设置于心墙与过渡层之间。接触面单元在选取参数时必须考虑厚度的影响。

3)本文基于有限元方法和直剪物理试验，建立了一套标定不同厚度接触面单元参数的数值试验方法。依托茅坪溪土石坝的物理试验成果，验证了该方法的可行性。

4)本文采用上述数值试验方法，依托茅坪溪土石坝的物理试验成果，研究了接触面单元各项参数与厚度的相关关系。研究结果显示，相关性中Gs最大，cs、φs及Rfs较小，ns基本无关；cs、φs及Gs与厚度正相关，Rfs为负相关。究其原因，接触面单元不仅反应了接触面的力学性能，也反应了部分过渡层材料的性能，因此随着接触面单元厚度的增加，过渡层的影响也逐渐增大，考虑到过渡层的力学性能强于接触面，造成接触面单元的力学性能逐渐增强。

5)本文将数值试验标定的接触面参数，应用于某沥青混凝土心墙坝的三维数值模型中，成功模拟了心墙和过渡层之间的相互错动。研究结果表明，接触面参数对心墙与过渡层之间的错动量具有显著影响，进一步论证了工程中需要根据接触面单元的实际厚度标定参数的重要意义。