空间差异性下沥青混凝土心墙坝应力与变形的随机有限元分析

郭 晴1,，刘东海，陈 辉

(1.河北工程大学 水利水电学院，河北 邯郸 056021； 2.天津大学 水利工程仿真与安全国家重点实验室，天津 300350 )

1 研究背景

沥青混凝土心墙坝是以堆石为主体，沥青混凝土作为坝体防渗材料的一种土石坝型[1]。沥青混凝土自身良好的防渗性、变形性和黏结性，以及沥青混凝土心墙可与坝壳材料同步施工以缩短施工工期的性能，使得沥青混凝土心墙坝在世界范围内的建设规模不断扩大，已经成为当前具有极大发展潜力的新型坝型。因此，有必要通过有限元计算，模拟大坝的运行情况，找出坝体可能存在的安全隐患，为大坝正常运行提供分析依据。涂幸等[2]、代凌辉[3]、Coroller[4]基于有限元分析软件对沥青混凝土心墙坝进行了三维有限元分析，探讨了在不同工况下坝体的应力、变形特性；张芸芸等[5]、Valstad等[6]对考虑堆石料湿化效应或受地震动荷载作用下的沥青混凝土心墙坝进行数值模拟，提出了改善坝体受力特性的有效方案。

但是，目前国内外的研究在采用有限元分析时多是假定同一分区采用相同的坝料力学参数。大坝实际施工质量在空间上存在差异性，这种差异性势必会导致坝体材料物理力学参数在空间分布上的不确定性，如果按坝体材料的设计参数计算分析大坝的应力和变形将会与实际情况存在差异。如杨鸽等[7]采用局部平均细分法模拟了二维的堆石料物理力学性质随机场，表明忽略材料的不确定性可能导致大坝的地震反应被低估。

本文采用考虑空间差异的随机有限元分析方法，研究沥青混凝土心墙坝的应力、变形情况；统计坝体应力、应变的空间分布规律，进而为大坝安全分析与设计优化提供参考依据。

2 沥青心墙坝随机有限元分析方法

2.1 坝体随机有限元计算流程

随机有限元法是将蒙特卡罗(Monte-Carlo)技术与有限元分析方法相结合，通过有限次的循环计算参数统计量，进而研究参数的分布特征，并且能够保证概率结果具有较高的置信度。在模拟的过程中，每一个仿真循环是完全独立的，而且任何一组循环与其他仿真循环结果无关，因此非常适合并行计算[8]。本文采用蒙特卡罗直接抽样法模拟坝体在施工过程中任何材料参数下的响应，一个仿真循环得到坝体在某组材料参数下的应力、变形等响应。随机有限元计算流程如图1所示。

(1)

2.2 材料物理力学参数的不确定性

在结构的随机有限元计算中，涉及的不确定性因素有：几何不确定性、初始条件不确定性、材料参数不确定性、边界条件不确定性，以及载荷不确定性等[9]。对于沥青混凝土心墙坝，它的几何模型大小、边界条件和受力情况都基本不变化。天然条件下存在地基覆盖层材料的非均匀性与复杂性，因此，本文将研究考虑坝体与覆盖层材料参数的不确定性对沥青混凝土心墙坝应力、变形的影响。

由于沥青混凝土与堆石料都是典型的散粒体材料，具有材料非线性特点，因此材料的本构模型选取邓肯张E-B模型[10]。将密度ρ与邓肯张E-B模型材料参数(黏聚力c、内摩擦角φ、破坏比Rf、弹性模量数K、体积模量指数m、弹性模量指数n、体积模量数Kb)作为随机变量进行抽样，Kur的取值是根据K的概率分布某次抽样后得到的数值，将2K赋值给Kur，即Kur=2K。然后再将这8个材料参数与密度作为该次抽样的堆石料材料参数赋值给堆石料的单元，从而进行该次有限元分析。

由于沥青混凝土心墙厚度较薄，一般为0.5～1.2 m，并且作为大坝防渗结构，其施工质量控制严格，心墙料物理与力学参数差异不大。根据实际工程中的现场检测，沥青混凝土心墙的技术指标大都符合设计控制指标[11]。因此，本文只考虑坝体除心墙外的过渡层区、堆石区和覆盖层区的密度与邓肯张E-B模型材料参数的不确定性，而心墙区的材料物理力学参数仍旧使用设计值。

2.3 模型参数批量赋值方法

为了使得有限元分析充分反映大坝施工压实质量的空间差异性，本文采用每个有限元单元对应不同的邓肯张E-B模型参数的方法。有限元分析软件ABAQUS可以通过调用全命令流式的inp文件，直接进行结构的有限元分析，因此可以通过改变inp文件中的命令来实现对大坝所有网格材料的重新赋值。批量赋值的过程采用开发的Fortran程序来实现，具体流程如下：

(1)确定赋值区域单元号E。通过ABAQUS软件建立坝体不同分区的集合，从inp文件中采集堆石区、过渡层以及覆盖层的单元号E。

(2)定义单元集合FE。在Fortran程序中读入单元号E，按照软件固定的格式循环，实现每个单元定义为一个单元集合的模式，生成任意单元集合FE，将结果输出到文本文件1.txt。

(3)定义单元材料截面属性。在Fortran程序中将邓肯张E-B模型材料属性的名称ME分配给相应的单元集合FE，按照软件固定的格式进行循环赋值，实现每个单元集合对应一个截面属性，将结果输出到文本文件2.txt。

(4)定义单元材料特性。读取2.2节每个单元随机抽样的密度ρi与邓肯张E-B模型参数ci，φ0i，Rfi，Ki，mi，ni，Kbi，Δφi，Kuri，Pa(Pa为大气压，Pa=100 kPa)。对材料名称ME、密度ρ和邓肯张E-B模型参数进行三重嵌套循环，实现对有限元模型任意单元的力学模型参数赋值，将结果写入文本文件3.txt。

(5)将文本文件1.txt、2.txt、3.txt合并为一个inp文件，通过ABAQUS软件直接调用不同抽样结果下的inp文件进行有限元计算。

3 空间差异性下的坝体应力、变形达标判断准则

在进行坝体应力、变形有限元计算时，随着坝体的填筑，荷载增量逐步施加，部分坝体单元的应力结果可能会超过其材料的极限应力而产生破坏。基于上述的坝体破坏方式和有限元计算方法，在进行坝料参数不确定时的坝体应力、变形有限元计算时，可以根据坝体材料参数在选取设计值时的确定性计算结果作为每一次循环模拟计算结果是否达标的判断标准。

根据随机有限元每次的循环计算结果，得出每次计算的坝体最大沉降、上下游水平位移最大值、大小主应力最大值,分别与材料参数设计值计算得出的相应最值进行比较，若小于等于设计阶段的计算结果，则说明达标，否则，说明不达标。判断式如式(2)所示。

(2)

根据前述判断准则确定坝体应力、变形最大值中没有超过设计情况计算结果的次数m占总有限元模拟次数M的比例，由此确定在材料参数不确定性下的坝体应力、变形达标概率P，如式(3)所示。进一步可以得出超标概率P′，如式(4)所示。

P=m/M；

(3)

P′=1-P。

(4)

4 工程实例分析

某沥青混凝土心墙坝最大坝高82.5 m，上下游坝坡比均为1∶1.8，采用渐变式沥青混凝土心墙，在心墙上下游侧设置有厚度为4.5 m的过渡层。图2中标明了大坝材料分区，包括堆石料Ⅰ区、堆石料Ⅱ区、过渡层区、沥青混凝土心墙区。坝址区河床存在厚度为42～50 m的覆盖层。

图2 大坝材料分区与筑坝顺序Fig.2 Materials partition and building sequence of dam

图3 大坝整体三维有限元网格剖分Fig.3 Three-dimensional finite element meshing of the whole dam

4.1 计算模型

利用ABAQUS有限元模拟软件提供的C3D8单元对模型进行网格剖分，共剖分27 142个单元，大坝三维整体网格剖分见图3。基础边界采用截断选取，竖直方向向下截取50 m，并在其底部施加固定位移约束；水平向截断长度为50 m，并在其截断面上施加固定位移约束。通过坝基地应力平衡和9级加载步模拟坝体的填筑过程(见图2)。计算中规定沿坝轴向从右岸到左岸为x坐标正向，沿坝体高程方向规定为y坐标正向，垂直于坝轴线从上游到下游为z坐标正向。

表1 材料随机参数统计特性Table 1 Random properties of dam material parameters

4.2 随机有限元计算结果分析

4.2.1 材料随机参数的统计特性

根据陈辉等[12]先前对堆石料材料参数不确定性的研究，结合现有土石坝填筑标准和相似工程现场检测资料，拟定不同分区邓肯张E-B模型材料参数的变异系数(即标准差/均值)如表1所示，各参数的期望值根据设计参数确定。由于现有的资料只有针对某堆石坝主堆石区变异系数的研究，并且考虑到变异系数只是反映数据的离散程度，本文想通过变异系数与该工程的均值真值求出标准差进行抽样，而且虽是不同分区但同为堆石料，变异系数差异性较小，所以不同分区采用了同一个变异系数。

由不同材料分区的8个参数统计特性，假定其服从一定的分布规律，对材料物理力学参数进行随机抽样。文献[12]统计了实际压实质量下堆石坝主堆石区的邓肯张E-B模型参数概率分布，认为邓肯张E-B模型参数服从正态或对数正态分布。

4.2.2 坝体随机有限元计算结果分析

绘制出由蒙特卡罗模拟得到的坝体最大沉降的滑动平均值随蒙特卡罗模拟次数的变化曲线，如图4所示。可知当蒙特卡罗模拟次数达到300次后计算结果已经较为平稳，因此下文随机模拟次数取足够大的600次。

图5是设计工况下坝体应力、变形等值线。结合有限元计算结果可知，坝体竣工期最大沉降为1.130 m，占坝高与覆盖层厚度的0.85%，位于坝体的心墙中下部。坝体向上游水平位移最大值为0.412 m，向下游水平位移最大值为0.383 m，分布较为对称，没有呈现明显的不均匀性。坝体大主应力从坝顶向坝基呈现出逐渐增大的趋势，其最大值为2 530 kPa，位于大坝底部。

注：频数为变形或应力在600次随机有限元计算中出现的次数。图6 坝体应力、变形概率分布Fig.6 Probability distribution of dam stress and deformation

图6为坝体应力、变形概率分布。随机有限元计算得到的坝体最大沉降概率分布如图6(a)所示。可以看出在考虑了堆石料与覆盖层的空间差异性后，与设计情况相比坝体最大沉降超标概率为50%，变异系数为9.6%(0.109/1.132×100%=9.6%)。随机统计结果的最大沉降为1.43 m，占坝高与覆盖层厚度的1.0%(坝高与覆盖层厚度为82.5 m+50 m=132.5 m，1.43/132.5×100%=1%)，虽然其发生概率仅有0.17%(通过600次随机有限元计算，统计最大沉降1.43 m出现的次数与600的比值可得)，但仍然会对坝体结构安全造成一定的威胁。对坝体沉降最大值进行正态分布K-S检验，sig值为0.744。由此可知坝体沉降最大值服从均值为1.132、标准差为0.109的正态分布。(sig值为统计显著性，由SPSS数据分析软件直接求得，在进行正态分布K-S检验时，若sig值>0.05，该样本服从正态分布，否则不服从正态分布。)沉降样本最大值1.43 m落在3σ区间内，因此在设计阶段很有必要考虑坝体沉降在3σ区间内的情况，对于本工程即考虑坝体最大沉降范围在0.803～1.457 m时坝体的应力、变形情况。

坝体的顺河向最大水平位移及大主应力最大值的概率分布如图6(b)—图6(d)所示，最终统计结果如表2。可看出考虑筑坝材料与覆盖层料的不确定性确实会使大坝应力、变形产生一定程度的离散。假如在设计阶段忽略这种不确定性，而仍用确定性的分析方法对结构进行分析，则有50%的可能性会低估坝体最大沉降，46%～47%的可能性低估上下游最大水平位移，43%的可能性低估坝体大主应力最大值，51%的可能性会低估坝体小主应力最大值。

表2 坝体应力、变形最大值随机有限元分析结果统计Table 2 Statistical results of dam stress and deformation by stochastic finite element analysis

4.2.3 心墙随机有限元计算结果分析

设计工况下心墙应力、变形分布如图7所示。结合有限元计算结果可知，心墙竣工期竖直沉降最大值为1.130 m，发生于心墙中下部。大主应力分布呈现出从顶部到底部逐渐增加的趋势，最大值为1 740 kPa，位于心墙底部位置，且为压应力，说明心墙受力状态良好。图8是心墙顺河向水平位移沿相对高程分布情况，可知设计工况下向上游水平位移最大值为0.050 m，主要位于心墙下部，向下游水平位移最大值为0.025 m，主要位于心墙上部。而随机有限元模拟向上游水平位移最大值均值为0.052 m，向下游水平位移最大值均值为0.023 m，模拟均值的分布情况与设计值相同。

图7 设计工况下心墙应力、变形等值线Fig.7 Stress and deformation of core wall in design condition

图8 心墙顺河向水平位移沿高度分布Fig.8 Horizontal displacement of core wall along the river against elevation

图9 心墙特征点位置Fig.9 Typical locations on core wall

在心墙典型剖面的上中下部位各选择一个特征点(有限元单元编号分别为536，5724，1844)对其应力、变形分布规律进行分析。特征点选取位置如图9所示。由于心墙向上下游方向水平位移较小，对大坝的正常运营不会产生影响，在此不再分析心墙顺河向水平位移的随机有限元分析结果。

图10 编号5724特征点应力、变形概率分布Fig.10 Probability distribution of stress and deformation at characteristic point No.5724

心墙上3个不同部位的特征点竖直沉降在以设计工况为判断标准时，都是有50%的超标概率，规律比较明显。相比之下，心墙的大、小主应力也有不同程度的超标，但无明显规律。从3个点应力、变形结果的整体来看， 3个特征点位置的不同，心墙上部的特征点超标概率较小，而心墙下部的特征点超标概率较大。从确定性有限元分析可知，心墙的竖直沉降与大小主应力最大区域更靠近心墙中下部，因此分析这种情况可能是由于下部点距离心墙的竖直沉降和大小主应力变化最大的区域更近导致的，所以在心墙施工时应对其中下部进行严格的质量控制。图10是编号5724特征点随机有限元计算得到的部分应力、变形概率分布。心墙特征点应力、变形随机有限元结果统计如表3所示。

表3 心墙特征点应力、变形随机有限元分析 结果统计Table 3 Statistical results of stress and deformation at characteristic points of core wall by stochastic finite element analysis

这一现象也可以从536与1844特征点的竖直沉降不服从正态分布看出。由于这2个特征点位距离心墙竖直沉降最大区域较远，所以其沉降样本点取值多等于该点沉降均值，没有产生过度的离散，因而不服从正态分布。可见虽然本构模型参数服从正态分布，但其应力、变形结果可能因研究位置的不同，概率分布也会不同。

5 结 论

本文针对大坝实际施工过程造成的材料参数空间差异对沥青混凝土心墙坝应力、变形存在一定影响的问题，给出了随机有限元分析的流程步骤和实现空间差异性计算的方法，并且结合工程实例，利用蒙特卡罗概率设计方法探讨了在考虑堆石料与覆盖层的空间差异性后，坝体和心墙特征点的应力、变形统计规律。具体结论如下：

(1)基于随机有限元分析各仿真循环相对独立的特点，提出了在有限元软件中实现材料参数空间差异性的方法，并且结合坝体破坏方式给出了判断坝体随机有限元结果达标的判断准则。

(2)考虑坝体材料的空间差异性确实会使大坝应力、变形发生一定程度的离散。而如果忽略这种差异性，仍用确定性分析方法对结构进行分析，则有50%左右的可能性会低估坝体的应力、变形。其中沉降最大值有可能会对坝体正常服役产生威胁，所以很有必要分析在随机有限元计算结果的最大值下，坝体结构是否仍然安全，从而指导工程设计。由于坝体上下游水平位移变异系数较小且大小主应力全为压应力，因此二者对坝体的安全运行威胁较小。

(3)空间差异性下沥青混凝土心墙上部主应力的超标概率小于下部主应力超标概率，不同位置处沉降的超标概率都为50%。且表征大坝或心墙的性能指标不一定都服从正态分布，即使该指标的最值服从正态分布，但有可能因为指标考察的结构位置不同而使其概率分布规律改变。