考虑覆盖层地基材料空间变异性的沥青混凝土心墙坝地震响应研究

罗博华， 宋志强， 王 飞， 刘 琛， 刘云贺

(1. 西安理工大学 水利水电学院，西安 710048；2. 中铁第一勘察设计院集团有限公司，西安 710043)

随机性存在于岩土工程的许多方面。其来源可能是由于对参数和荷载难以进行精确的测量所致[1]。对于河床覆盖层，由于类型复杂、结构松散，因此物理力学等性质呈现出较大的不均匀性。在以往对于建立在深厚覆盖层上大坝的数值计算中，材料参数通常选取设计值，忽略了不确定性的存在，从而影响大坝地震响应的计算结果[2-4]。

本文应用随机场理论，在蒙特卡洛法的框架中[13]，建立了基于Cholesky分解的高斯空间随机场离散方法，讨论了覆盖层材料参数具有空间相关性和差异性时，覆盖层上沥青混凝土心墙坝地震响应对于覆盖层材料等效线性模型参数的敏感性，以概率思想对坝顶水平峰值加速度、心墙顶水平峰值加速度和竖向永久变形的数值计算结果进行统计分析，采用Anderson-Darling检验来探讨统计结果的分布形式。

1 覆盖层材料动力参数敏感性分析

参数敏感性分析可用来为随机参数的选取提供参考依据。本文土石坝静力分析时采用邓肯-张E-B非线性弹性模型；动力响应计算采用等效线性方法；永久变形采用沈珠江等[14]建议的修正等效线性黏-弹性模型[15-16]。

在静力分析中，对于邓肯-张E-B模型参数敏感性的研究较为充分，即使学者们关心的敏感性指标不同，也能够发现，土石坝静力参数ρ，K，φ，n，Rf对结果的影响较为显著[17-19]，结合以往研究中对随机静力参数的选取[20]，本文考虑密度ρ，系数K和初始摩擦角φ为随机参数，不再开展覆盖层材料静力参数敏感性分析。

在计算峰值加速度的等效线性方法中，土体剪切模量和阻尼比均与剪应变有关，如式(1)～式(4)所示

(1)

(2)

(3)

(4)

式中：G和λ分别为剪切模量和阻尼比；σ′m为平均有效应力；k1，k2，n为试验确定的材料参数；λmax为最大阻尼比；γmax为最大动剪应变；pa为大气压强。

在永久变形的计算中，地震的持续时间平均分为几个时间段，残余体应变和残余剪应变增量形式为

(5)

(6)

式中： Δεvr， Δγr分别为残余体应变、残余剪应变增量；γd为残余剪应变；S1为剪应力水平；N，ΔN为振动次数及其增量；c1～c2为试验参数。可以看出，静力计算的结果围压σ′m为动力计算的初始条件，动力计算结束后得到的动剪应变γd为永久变形的初始条件。为避免随机参数过多引起的互相干扰，本文对c1～c5取确定值。

目前对于动力参数敏感性研究不多，尚无较为普遍认可的结论，为了确定覆盖层材料的动力随机参数，采用正交试验法对动力参数进行敏感性分析。由极差差异评价各因素影响的敏感性，综合考虑较为敏感的参数作为随机参数，并根据敏感性结果拟定随机参数的离散方式和变异系数。

1.1 正交试验设计

试验方案选取各参数具有代表性的组合，根据“正交表”设计试验。试验中将对地震响应指标：坝顶水平向峰值加速度和竖向永久变形具有影响的参数称为因素，针对因素进行比较的试验条件称为水平。正交表Ln(mk)中：L为正交表；n为试验次数；m为因素的水平数；k为因素个数。

1.2 极差分析

假设j=A,B,…为对考察指标有影响的各个因素；i=1,2,…,r为各因素的水平数；Ai为因素A的第i水平；Xij为影响因素j的第i水平值。 在Xij下进行试验，得到正态分布的随机变量Yij；在Xij进行n次试验得到n个试验结果Yijk(k=1,2,…,n)，统计参数Kij为

(7)

Rj=max{K1j,K2j,…}-min{K1j,K2j,…}

(8)

式中，Rj为第j列因素的极差。

1.3 模型概况

某沥青混凝土心墙堆石坝最大坝高102.8 m、坝顶高程 652.8 m、坝顶宽 10 m。上游坝坡在高程 597.0 m以上坡度为 1∶1.7，高程 597.0 m 以下坡度为 1∶2.5；下游坝坡度为 1∶1.7，在高程 617.8 m和582.8 m处各设2 m宽马道。坝体分区由堆石区、围堰石渣区、过渡层、沥青混凝土心墙组成。第坝体施工期填筑分24步：第1～第6步填筑石渣材料；第7～第21步填筑堆石区材料；沥青心墙和堆石同时填筑，水荷载分三次于第22～第24步添加于心墙，至正常蓄水位。

沥青心墙与过渡料，防渗墙与覆盖层之间设置Goodman接触。为研究覆盖层材料参数随机对坝体响应的影响，假设大坝坐落在100 m高的覆盖层上，覆盖层分别向两侧和深度方向延伸1倍坝高。坝体材料参数来源于实际工程[21-22]，覆盖层材料参数根据相似工程拟定[23]。参数信息如表1和表2所示，防渗墙为混凝土材料，密度为2.4 g/cm3，弹性模量为28 GPa,泊松比为0.167。坝体采用平面四节点等参单元模拟。由于覆盖层底部坐落于基岩表面，基岩变形量小，因此静力计算中，将覆盖层底部的所有节点自由度全部约束，两侧节点在水平方向约束。

表1 静力材料邓肯-张E-B参数Tab.1 Static material Duncan-Zhang E-B parameters

表2 动力材料及永久变形参数Tab.2 Dynamic materials and permanent deformation parameters

坝体的模型信息及二维有限元网格，如图1所示。模型总单元个数3 048个，最大单元格长度为15.26 m，宽度为4.76 m。

图1 模型及随机场网格信息Fig.1 Model and grid information of random field

动力计算中，地震波选取采用太平洋数据库中的Kobe波，最大峰值加速度0.331 4g，地震动加速度时程曲线如图2所示。边界条件采用无质量地基固定边界。地震响应分析中，地震波从覆盖层底部输入，竖向震动峰值加速度为水平输入的2/3，持时25.25 s,输入间隔0.02 s。

图2 地震动加速度时程Fig.2 Time history of ground motion acceleration

1.4 计算及结果分析

本文设5个试验水平，对表2中覆盖层材料动力参数k1，k2，n，v，λmax分别乘以0.70，0.85，1.15和1.30，作为各参数正交试验的不同因素水平，如表3所示，其余分区材料参数采用表2中的固定值。

表3 正交试验不同因素水平取值Tab.3 The value of different factors in orthogonal test

先采用设计参数开展静力计算，获得动力计算的初始围压，作为地震响应计算的初试条件，每个单元的围压为固定值。根据正交设计，制定L25(56)正交表。分别计算各工况下坝顶水平向峰值加速度和竖向永久变形。试验结果如表4所示，对表4中对各影响因素进行极差分析，分析结果如表5所示。

表4 正交试验方案及试验结果Tab.4 Orthogonal test plan and test results

表5 试验方案及动力响应极差分析Tab.5 Experimental scheme and range analysis of dynamic response

对试验指标坝顶水平向峰值加速度和竖向永久变形的极差分析结果进行整理，按照各因素对各试验指标的极差值大小绘制极差值柱状图，如图3所示。由于极差最大的因素为敏感性最大的因素，从表5和图3看出，对坝顶水平向峰值加速度而言，k1，λmax，k2这3个参数的敏感性较高，对计算结果影响显著，而v，n对各指标的影响不显著，参数敏感性较低；对残余变形而言，k1，n，v敏感性较高，而k2，λmax对各指标的影响不显著，敏感性较低。综合考虑本文选取k1，n，v，λmax为随机参数。

图3 覆盖层材料动力参数敏感性对比结果Fig.3 Comparison result of dynamic parameter sensitivity of cover materia

2 材料参数空间变异性的模拟技术

2.1 随机场的实现方法

随机场模型中采用变异系数、相关函数、相关距离等特征来反映土体性质在空间上的特点。本文采用基于Cholesky分解的协方差分解法离散随机场。结合独立标准正态分布、对数正态实现覆盖层参数空间变异性的抽样模拟。各参数的随机场单元参数的统计特性可由均值μ和标准差σ表征，变异系数为

(9)

自相关距离是VanMarcke引入的一个指标，表明相关距离内的两种土壤性质是否具有明显的自相关关系，在相关距离之外，可以认为土壤性质基本上是不相关的。自相关函数可以用来计算土体任意两点间的相关系数。由于高斯自相关函数是比较稳定的，应用范围广，本文用高斯自相关函数来描述各单元形心处材料参数之间的空间相关性。

(10)

式中:τx为空间中两点水平方向上的距离，τx=|xi-xj|；τy为空间中两点在竖直方向上的距离，τy=|yi-yj|；δh和δv分别为水平和竖直自相关距离，可根据土体的性质进行拟定，详见2.2节描述。

离散后的随机场是一组相关随机变量，对这组相关随机变量进行抽样模拟，就实现了对材料参数空间变异性的模拟。根据单元形心坐标和自相关函数确定随机场，离散后得到的一组相关随机变量的相关系数矩阵[24]

ρn×n=[ρ11,ρ12,…,ρ1i,…,ρ1n]

(11)

(12)

(13)

(14)

(15)

(16)

式中，σln xi，μln xi分别为第i个单元变量取对数后的标准差和均值。

2.2 覆盖层材料参数随机场的建立

根据文献[26]中对于土体相关距离的统计分析，水平相关距离为10～80 m，而竖向波动距离为1～5 m，且相关距离在水平和竖直方向通常相差在一个数量级左右。由于在数值计算中，相关距离需要小于随机场单元长度，为保证计算精度，同时避免过量消耗计算机资源，结合模型信息，本文计算水平自相关距离和竖向自相关距离分别取50 m和5 m。

随机场的离散方式采用协方差分解法，各参数的分布形式采用正态分布或对数正态分布，变异系数均偏保守的取0.1～0.2。由于数值计算舍入误差和其他误差的原因,可能出现相关系数矩阵无法进行Cholesky分解的情况，这种情况下行列式进行微调整，使能够产生与实际相关系数矩阵接近的随机变量抽样序列[27]。

一次覆盖层材料参数K随机时的坝体系统材料分布示意图，如图4所示。可以看出当水平相关距离和竖直相关距离不同时，非均匀性呈现出向相关距离大的方向延展的特性。

图4 一次参数K随机时坝体系统材料分布场Fig.4 Material distribution field of dam body system when parameter K is random

3 覆盖层参数随机沥青混凝土心墙坝有限元计算

3.1 计算流程

随机有限元动力计算流程分为三个部分：首先确定随机参数及特征，通过程序生成材料随机参数库,再将随机参数赋予模型覆盖层材料，进行静动力计算。最后对随机有限元的结果进行统计分析，与材料参数确定时的结果进行比较，得出结论。计算流程如图5所示。

图5 随机有限元地震响应分析流程图Fig.5 Flow chart of random finite element seismic response analysis

为判断蒙特卡洛模拟次数N是否足够，每次模拟结束后，输出坝顶峰值加速度和竖向永久变形的结果，并根据式(17)计算平均值，绘制模拟次数和各指标平均值的曲线图，趋势平稳时则说明模拟次数足够。

(17)

式中：X为地震响应类型； 下标ave为均值；n为模拟次数；xi为第i次模拟的响应值。

3.2 参数分布方式及变异系数的选取

本文考虑密度ρ，系数k和初始摩擦角φ为随机参数，为避免夸大随机程度的影响，变异系数偏保守地选为15%，以上参数均服从对数正态分布。

根据动力参数敏感性的分析结果，综合考虑选取k1，n，v，λmax为随机参数。根据文献[28-29]，为了弱化敏感参数的影响，在进行分布方式及变异系数确定时，当参数的变异系数取值较大时，考虑分布方式采取对数正态分布。由于地震响应结果对k1，λmax相对更敏感，拟定其服从对数正态分布，变异系数取0.2，拟定n，v服从正态分布，变异系数分别取0.15和0.10。随机参数的选取及分布方式，如表6所示。

表6 覆盖层随机参数的选取及分布方式表Tab.6 Selection and distribution of random parameters of cover layer

3.3 计算工况

为了研究静动力参数在随机情况下对动力响应的影响，本文考虑三种工况进行研究：①只考虑静力参数随机，动力参数采取设计值；②只考虑动力参数随机，静力参数采取设计值；③静动力参数同时随机。不同工况下的计算流程略有差异。工况①在随机静力计算结束后，分别将围压值赋给模型作为地震响应计算的初始条件，动力参数采用设计值；工况②和前述敏感性分析类似，地震响应计算的初始围压均相同，动力参数从随机参数库中调取进行计算；工况③的随机静动力参数和工况①、工况②采用同一套，即N组静力随机参数取值和工况①一一对应，动力参数和工况②一一对应，而不是重新抽样，分解进行随机参数计算。

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

根据式(17)，绘出在工况③情况下，坝顶水平向峰值加速度和竖向永久变形平均值随MC模拟次数的变化，如图6所示。工况①和工况②的趋势和工况③类似，模拟次数达到300次后动力响应的计算结果已经较为平稳。因此下文对500次MC模拟的结果进行分析以估算总体的统计特性。

图6 坝顶地震响应随模拟次数变化趋势Fig.6 Variation trend of dam crest seismic response with simulation times

4.1 覆盖层参数随机对峰值加速度的影响

随机有限元计算得到的坝顶峰值加速度的概率分布情况，如图7所示。图7中实线表示随机参数下统计均值结果；虚线表示确定参数下的计算结果。可以看出，三种工况坝顶最大峰值加速度均值分别为7.003 m/s2，7.118 m/s2和7.123 m/s2，均大于确定参数下的计算结果6.99 m/s2，工况③的结果与均值的差距最大。

图7 坝顶水平向峰值加速度概率分布Fig.7 Probability distribution of horizontal peak acceleration of dam crest

将随机计算的地震响应结果分别与确定参数计算得到的结果进行比较，若大于确定值下的计算结果，则定义为结果超标。由此确定在材料参数不确定性下的坝体地震响应超标概率P=m/M，表示坝体地震响应中超过特征结果的次数m占总模拟次数M的比率。在累积曲线上查得与确定值下的结果对应的累积概率值，即累积曲线与图中虚线的交点，三种工况下超标概率分别为54.8%，79.4%，78%，随机计算的结果均有超过50%的可能大于确定值结果。

从图7直方图可以看出，工况①中静力参数随机的情况对坝顶加速度放大系数的影响不如工况②、工况③显著，但其均值仍然大于确定参数下的计算结果。随机有限元计算得到的最大峰值加速度出现在工况③中，为7.661 m/s2,比确定参数的计算结果增长了9.4%，最小峰值加速度出现在工况①中，为6.546，比确定值的结果减小了6.53%。

图8为对结果的Anderson-Darling正态分布检验。样本点集中于45°直线附近，且相伴概率p均在0.05以上。检验结果表示，三种工况中坝顶峰值加速度的结果均服从正态分布。

图8 坝顶水平向峰值加速度正态分布检验Fig.8 Test of normal distribution of horizontal peak acceleration of dam crest

随机有限元计算得到的心墙顶峰值加速度的概率分布情况，如图9所示。心墙顶的水平向峰值加速度略小于坝顶，这与沥青心墙节点和坝顶峰值加速度的节点相邻有关。三种工况下的均值分别为6.992 m/s2，7.021 m/s2，7.12 m/s2，均大于确定参数下的计算结果，心墙顶峰值加速度最大值同样出现在工况③，为7.772 m/s2，最小值出现在工况①，为6.727 m/s2，分别比确定值下的计算结果增大了11.3%，减小了3.8%。在累积曲线中查得超标概率分别为54%，61.6%，78%，和坝顶水平向峰值加速度的分布规律相似，工况①的超标概率最小，工况②、工况③稍大。随机有限元计算得到的心墙顶峰值加速度的Anderson-Darling正态分布检验结果，如图10所示，在工况①的检验中，相伴概率小于0.05，不服从正态分布，继续进行对数正态分布检验，相伴概率为0.758，结果符合对数正态分布，工况②、工况③的心墙顶水平向峰值加速度服从正态分布。

图9 心墙顶水平向峰值加速度概率分布Fig.9 Probability distribution of horizontal peak acceleration at the top of core wall

图10 心墙顶水平向峰值加速度正态分布检验Fig.10 Test of normal distribution of horizontal peak acceleration on top of core wall

4.2 覆盖层参数随机对坝顶竖向永久变形的影响

随机有限元计算得到的坝顶最大竖向永久变形的概率分布情况，如图11所示。三种工况结果的均值分别为36.484 cm，36.413 cm和36.716 cm，均大于确定值的计算结果36.35 cm，从累积曲线上查得超标概率分别为72.4%，59.6%，93%。从图11中实线和虚线之间的距离可以看出，工况③的结果与均值的差距最大，工况②的结果与均值的差距最小。即同时考虑覆盖层材料静动力参数随机时，计算结果超出确定值的概率最大，与确定值相差也最大。但这些差距在1%以内，相比于峰值加速度与确定值的最大差距达到11%，覆盖层材料动力参数随机对坝顶竖向永久变形的影响较小，这与随机计算中直接影响永久变形的参数c1～c5取确定值有关。但其均值也都超过了确定值的结果，说明考虑参数随机情况时，三种工况的计算结果大概率会超出确定值，使工程更偏于不安全。对坝体最大竖向永久变形的Anderson-Darling检验，如图12所示。检验结果表示，三种工况中坝体最大竖向永久变形均服从正态分布。

图11 坝顶竖向最大永久变形概率分布Fig.11 Probability distribution of vertical maximum permanent deformation of dam crest

图12 坝顶竖向最大永久变形正态分布检验Fig.12 Inspection of normal distribution of vertical maximum permanent deformation of dam crest

将各工况下动力响应及正态检验结果，如表7所示。超标概率及分布形式汇，如表8所示。从表中可以看出，对于峰值加速度和最大竖向永久变形而言，三种工况的统计结果服从正态分布或对数正态分布；从均值与确定值之间的差异来看，三种工况的随机计算结果均值均超过了确定参数的计算结果，说明考虑覆盖层材料参数随机时，不论是静、动力参数分别或同时随机，随机计算的结果大概率会大于确定参数下的计算结果；从变异系数的大小来看，峰值加速度的变异系数均大于竖向永久变形，这说明覆盖层材料参数随机对坝顶峰值加速度的影响比对永久变形的影响更为显著；从超标概率的大小来看，超标概率最大出现在计算工况③竖向永久变形时，且超过了90%，说明覆盖层材料静动力参数同时随机会有很大的概率使竖向永久变形的计算值超过参数确定下的计算值。

表7 地震响应结果Anderson-Darling正态检验Tab.7 Anderson-Darling normal test of seismic response results

表8 超标概率及分布形式汇总Tab.8 Summary of probability and distribution form

在对坝顶峰值加速度的影响显著性上，工况①<工况②<工况③；在对坝顶最大竖向永久变形的影响显著性上，工况②<工况①<工况③。总而言之，同时考虑覆盖层材料静动力参数随机时，对坝体地震响应的影响比考虑分别随机的情况更大。即使在1 500次随机计算中，地震响应与确定值的差距最大仅为11.3%，但仍然有70%的可能性会大于确定值的结果。

4.3 覆盖层-坝体系统中轴线的水平向峰值加速度

根据上述结论，工况③对峰值加速度的影响更为显著，且沥青心墙的动力响应和大坝峰值加速度分布趋势相似，因此本节提取出工况③沥青心墙和覆盖层中轴线所有节点的水平向峰值加速度，对每个节点计算500次结果的变异系数，绘制峰值加速度曲线及变异系数曲线图，如图13所示。

由图13可以看出，对于500次随机有限元计算，水平向峰值加速度的分布规律相同，但不同高程处节点峰值加速度的变异系数不同。对于心墙中轴线，变异系数最小出现在77 m高程处，随后增大再减小，在81 m处达到最大值，之后一直减小直至坝顶，坝顶水平向峰值加速度的变异系数在高程方向处于中间值。对于覆盖层中轴线，变异系数最小出现在覆盖层底部，随后沿高程扩大，在覆盖层中部偏上的位置达到最大值，再逐渐减小直至覆盖层顶部。

图13 中轴线峰值加速度及变异系数曲线图Fig.13 Curve of peak acceleration and coefficient of variation of the central axis

对沿高程方向的节点进行Anderson-Darling正态分布检验，发现并不是所有数据都服从正态分布或者对数正态分布；覆盖层的峰值加速度变异系数总是大于坝体，这与随机场的范围选在覆盖层有关，覆盖层参数随机使结构动力响应的结果发生离散，且这种离散对参数随机的区域影响较参数确定的区域更为显著。

5 结 论

本文应用随机场理论，在蒙特卡洛法的框架中，建立了基于Cholesky分解的高斯空间随机场离散方法，讨论了覆盖层材料参数具有空间相关性和差异性时，沥青混凝土心墙坝系统的地震响应会如何变化，并对统计结果进行了分布检验。结论如下：

(1) 实现了基于蒙特卡罗法考虑覆盖层材料参数的空间变异性，提出了一种ABAQUS和Python生成材料随机参数库的方法，这种方法可以自动赋值，计算和提取结果，减少了动力计算中由于迭代造成的复杂性。计算中程序和模型相互独立：是一种“非侵入式”的随机场实现方法，在其他相似计算中可以推广和应用。

(2) 采用正交分解法，对覆盖层动力参数进行敏感性分析，极差最大的因素为敏感性最大的因素。对加速度放大系数而言，阻尼比随剪应变变化的参数k1、最大动剪切模量的参数k2和最大阻尼比λmax的敏感性较高，对计算结果影响显著；对永久变形而言，阻尼比随剪应变化的参数k1、最大动剪切模量随围压变化的参数n、和泊松比v的敏感性较高。

(3) 在三种工况累计1 500次随机计算中，地震响应的统计结果有70%的概率会超出参数确定下的计算结果，这表示忽略材料的不确定性可能导致低估大坝的地震响应。考虑覆盖层材料静、动力参数同时随机时，地震响应的超标概率大于75%，分别随机时超标概率超过50%，即同时随机对地震响应的影响更为显著。心墙顶部峰值加速度最大超出确定参数结果的11.3%，竖向永久变形最大超出确定参数结果的1%，表明覆盖层材料参数随机对加速度的影响大于对永久变形的影响。

(4) 在高程方向取足够多的样本进行正态检验时，随机有限元计算的统计结果并不一定都符合正态或对数正态分布。覆盖层材料参数的不确定性使结构动力响应的结果发生离散，对于覆盖层-坝体系统中轴线的水平向峰值加速度，覆盖层中轴线的变异系数最大为12.3%，心墙中轴线的变异系数最大为5.8%，这表明随机结果的离散在覆盖层地基区域更为显著。

(5) 本文计算的随机响应结果均在大坝的安全范围之内，这与选取的分布方式及变异系数偏保守有关。实际覆盖层部分参数的变异系数可能会大于0.2，同时空间变异性在筑坝材料中同样存在，采取不同的变异系数，综合考虑覆盖层和坝料参数的空间差异性会得到更符合实际的结果，这个问题作者正在进一步研究中。