考虑坝体材料参数随机性的拱坝易损性分析

徐 强,姚 文 斌,陈 健 云,张 天 然,曾 祥 兴

(1.大连理工大学 海岸与近海工程国家重点实验室,辽宁 大连 116024; 2.大连理工大学 工程抗震研究所,辽宁 大连 116024)

0 引 言

中国拥有丰富的水能资源,但是其空间分布十分不均匀,多分布于西南地区,而西南地区多为河谷狭窄的山区,因此,拱坝成为该地区大坝修建的主要坝型之一,库容往往能达到数十亿甚至上百亿m3[1]。然而,中国西部地区位于喜马拉雅火山地震带上,地震较为频繁,强地震发生时,坝体材料在较大的应力作用下易失效。一旦拱坝失事,将给下游造成巨大的生命财产损失。针对大坝破坏失事方面的模拟,随着计算机水平的不断提高和数值分析方法的不断完善,有限元分析方法逐渐成为研究坝体在不同工况下损伤破坏情况的重要手段。

由于运输、施工和筑坝材料属性等因素的影响,坝体不同部位的混凝土材料参数存在一定的空间随机性,然而,传统的有限元分析中,假定大范围的筑坝材料参数为一个确定的值,无法反映出坝体实际的材料分布情况,造成分析结果与实测值有一定差距。如果能将这一空间随机性在数值分析过程中充分模拟,分析结果会更能反映坝体的实际损伤分布,从而为拱坝的结构设计和加固提供依据。

目前,对混凝土材料随机性进行模拟的理论主要有随机变量理论和随机场理论。随机变量理论假定空间内任意点之间的材料参数相互独立,没有相关性,该理论的优点在于简单易于实现,缺点是忽略了结构不同部位之间的相关性。而随机场理论则充分考虑了材料参数的随机性和空间位置上的相关性,更能反映实际的结构参数特性。Ardebili等[2]以Koyna重力坝为研究对象,将混凝土弹性模量、质量密度和抗拉强度均假定为随机场,并采用协方差矩阵分解和中点离散技术生成有限元模型展开研究,研究结果发现混凝土的非均质性影响着结构的抗震性能评价;Carl和 Padgett等[3]在对大型混凝土重力坝的可靠度分析中,同时考虑了地震动和坝体材料的随机性;陈健云等[4]在对Hardfill坝的有限元分析中同时考虑材料密度、拉伸弹性模量和抗拉强度等参数的空间随机性,并通过蒙特卡洛模拟方法进行分析,验证了材料参数的随机性对结构损伤的影响。

在研究坝体在不同地震动强度下的非线性响应时,需要对加速度时程反复进行调幅计算,计算量大且计算结果可能不收敛。国外学者Estekanchi提出了耐震时程法(ETA)[5],且详细阐述了耐震时程加速度曲线的合成方法,并验证了该方法在保证计算精度的同时可有效减少了计算量。白久林等[6]合成了基于中国抗震反应谱的ETA时程曲线,与增量动力分析方法(IDA)对比后,验证了ETA方法的可行性。

2003年,由中国学者李杰、陈建兵等[7]提出的概率密度演化方法(Probability Density-Evolution Method)近些年来逐渐发展完善,并得到很好的应用。该方法通过推导出解耦的广义概率密度演化方程,求解后得到响应全过程的概率信息[8]。概率密度演化方法对于非线性强、自由度数量高的复杂结构也有较好的适用性。陈健云等[9]结合广义概率密度演化理论,提出了一种混凝土坝全周期抗震性能分析方法,并与蒙特卡洛模拟方法的计算结果进行对比,验证了概率密度演化方法的可行性。

本文根据随机场理论建立了坝体材料参数随机的白鹤滩拱坝模型,根据标准设计反应谱生成了50组ETA人工加速度时程曲线,并作用于模型底部,分析了相关长度和参数平均值对坝体结构损伤分布的影响。选取坝顶位移、损伤体积比和横缝开度为损伤评价指标,对50组分析结果进行整理并建立其概率密度演化模型,根据概率密度曲线得出拱坝的可靠度,以评价白鹤滩拱坝在地震动作用下的安全性。

1 理论方法

1.1 混凝土材料空间变异性模拟

采用随机场对拱坝进行模拟时,首先要对随机场进行离散。合理的离散有助于提高模拟的精度。本文以断裂能和抗拉强度为随机变量,基于协方差矩阵分解法和中心点法生成随机场[2]。

假设随机场单元和有限元模型网格划分单元一致,共生成n组随机场。采用指数型自相关函数:

(1)

式中:ζ为两单元之间的距离,lcorr为材料的相关长度。坝体不同部位材料的相关性取决于相关长度lcorr,相关长度的含义是若坝体上任意两点之间的距离大于该长度,则这两点的材料参数相互独立,反之,如果结构上两点之间距离小于相关长度,则两点的材料参数值具有相关性,即相关长度越大,坝体材料越趋向于均值。

由离散化的自相关函数构造相关矩阵C:

C=(ρij)m×m

(2)

式中:m为模型单元数。

构造矩阵Q,并求其正态分布函数值:

(3)

V=F-1[Φ(Q)]

(4)

需要说明的是,V(i,j)的具体含义为假设材料参数服从标准正态分布,其值为负无穷到第i组随机场第j个单元的材料参数的分布函数值。

假设坝体材料参数服从韦伯分布,则i组随机场各单元参数的分布函数值取V(i,j),求出其逆累积分布函数,得出参数值。

1.2 耐震时程法(ETAM)

耐震时程法作为一种动力推覆过程,其优点在于保证预测结构响应精度的同时,可以很大程度减少计算量[10]。其思路是生成随时间逐渐递增的ETA加速度时程曲线,且不同周期对应的反应谱呈既定的倍数拟合关系。

ETA加速度时程的生成思路是先生成平稳加速度时程,经反复迭代优化后,使之各时段加速度对应的反应谱与目标反应谱有更好的拟合。在地震动持续时间内,不同时间段对应的目标反应谱与目标反应谱为倍数关系:

(5)

式中:TTarget为目标时间;SaC(T)为标准设计反应谱;t为拟合反应谱的时间;SaT(T,t)是0～t时刻的目标加速度反应谱。

根据加速度反应谱与位移谱的数学关系可知,目标位移反应谱为

(6)

式中:SuT(T,t)为时刻t的目标位移反应谱。

想要让目标加速度时程在每一点都能同时较好地拟合加速度反应谱和位移反应谱几乎是不可能的,但可以采用式(7)对加速度时程进行无约束优化,从而减小同目标反应谱之间的误差:

α[Su(T,t)-SuT(T,t)]2}dtdT

(7)

式中:ag是初始生成的地震动时程曲线,Tmax是地震动时间最大值;α是位移谱的权重系数,本文只考虑加速度反应谱的影响,因此α值取0;Sa(T,t)和Su(T,t)分别表示周期T下0～t时刻的加速度反应谱和位移反应谱。

图1为ETA时程和IDA时程的对比图。其中,ETA时程对0～5 s、0～10 s、0～15 s和0～20 s四个时段对应的反应谱进行优化拟合,并计算不同时段时程对应的反应谱如图2所示。得到的ETA时程曲线具有两方面的特点:① 随着时间的推移,其峰值加速度不断增大;② 不同时段对应的加速度反应谱与标准谱成倍数关系。

注:左为ETA法,右为IDA时程法。图1 ETA法与IDA法时程对比Fig.1 Comparison of ETA Method and IDA Method

图2 ETA时程加速度反应谱Fig.2 ETA time history acceleration response spectrum

1.3 概率密度演化法(PDEM)

通常,坝体在地震激励下的动力方程可以表示为

(8)

Z(t)=H(θ,t)

(9)

求导后,广义速度向量可表示为

(10)

由于没有其他随机因素的加入和消失,因此根据概率守恒原理经推导可得[12]:

(11)

当仅考虑系统中某一个分量时,式(11)的概率密度演化方程简化为

(12)

式中:pzθ(z,θ,t)为Z(t)与θ的联合概率密度函数。

根据初始条件和边界条件对偏微分方程进行求解,得到感兴趣变量的概率密度函数:

(13)

对于式(13)这种高维积分问题,很难得出解析解,因此,本文分别采用MOL算法和具有TVD性质的有限差分格式进行求解[13]。

1.4 直接概率积分法(DPIM)

若以θ和Y分别表示输入随机向量和输出随机向量,根据概率守恒原理可知:

(14)

其中,Pθ(θ)和PY(Y)分别表示θ和Y的概率密度函数,Ωθ和ΩY分别表示输入随机向量和输出随机向量的样本空间。通过推导,动力系统中任意时刻t的输出随机响应的概率密度函数为

(15)

方程(15)称为概率密度积分方程,由于该函数的求解过程中涉及到狄拉克函数以及高维度积分,导致该方程地求解极为困难。利用直接概率积分法可以高效地求解概率密度积分方程,得到感兴趣参数的概率密度函数。直接概率积分法包含两方面的重要技术:① 基于GF-偏差法选取代表点集并对概率空间进行剖分[14];② 将原方程中不连续的高斯函数替换为连续高斯函数,进行光滑化处理。运用下式计算输出随机响应概率密度函数[15]:

(16)

进一步推导后可得:

(17)

式中:θq表示第q个代表点;N表示所选代表点数量;Pq为该代表点的赋得概率;σ为光滑化参数。

2 随机参数的选取

考虑到拱坝坝体材料和地震动都具有随机性,因此以相关长度表征坝体材料的随机性,结构上距离超过相关长度的两点认为其随机材料属性互相独立,本文认为相关长度服从坝高10%～60%的均匀分布[1];以标准地震动反应谱的特征周期表征地震动的随机性,认为其服从0.20～0.35 s的均匀分布。

采用matlab编程在二维空间中选取代表点,并基于GF-偏差优化方法对选取点集进行优化。点集的GF值表征点集与理论分布的偏差程度。通过优化,点集的GF值从0.32下到0.28,优化效果良好。采用维诺图计算各代表点的赋得概率,点集分布如图3所示。

图3 代表点选取结果Fig.3 Results of representative point selection

3 拱坝有限元计算

3.1 有限元计算模型

以白鹤滩高拱坝为研究对象,建立其有限元模型如图4所示。其最大坝高为289 m,坝顶和坝底最大厚度分别为13 m和72 m,白鹤滩坝体-地基有限元模型由拱坝和具有辐射阻尼的岩石地基组成;坝体共设置了30条横缝。模型由25 776个六边形实体单元组成,其中包括4 656个坝体单元和21 120个地基单元。

图4 白鹤滩拱坝有限元模型Fig.4 Finite element model of Baihetan arch dam

模型考虑断裂能和抗拉强度服从威布尔分布,其他参数设置为定值。考虑自重影响,坝体混凝土密度为2 400 kg/m3,弹性模量为36 GPa,泊松比0.167;基岩密度为2 800 kg/m3,弹性模量26 GPa,泊松比0.24;静水压力方面,上游水位825 m,下游水位604 m;动水压力按照Westergaard附加质量进行施加,来模拟动水压力。

3.1.1混凝土本构模型

本文采用塑性损伤本构模型,其应力-应变关系如图5所示,在不同阶段,其关系式为

图5 塑性损伤本构的应力应变关系曲线Fig.5 Stress-strain curve of plastic damage constitutive

(18)

3.1.2拱坝坝体材料随机性的实现

根据不同参数对坝体抗震性能的影响程度,选取拱坝混凝土材料的抗拉强度和断裂能为随机变量,其他参数为定值。假定抗拉强度服从1.5～3.5 MPa的均匀分布,断裂能服从100～300 N/m的均匀分布,根据选取的50组代表点对应的不同相关长度,采用协方差矩阵分解法和中心点法生成50组随机场,按相关长度从小到大排列。图6(a)展示了由上述理论生成的白鹤滩拱坝随机场模型,从上至下相关长度在变化范围内逐渐增大。与图6(b)[2]所示的国外学者基于同一理论生成的Koyna重力坝随机场模型进行对比可以发现,其随机性分布较为一致。随着相关长度的增大,坝体材料分布逐渐趋向于均匀化,在较大范围内材料参数趋于一致,说明应用随机场理论能较好地对坝体材料的随机性进行模拟。

图6 不同相关长度对应的随机场Fig.6 Random fields corresponding to different correlation lengths

3.1.3横缝模型

白鹤滩拱坝共设置30条横缝,横缝分布如图7所示。在本研究中,采用接触边界模型模拟混凝土拱坝的横缝,横缝的接触关系[16]如图8所示,法向接触压力的计算公式为

图7 坝体横缝分布Fig.7 Distribution of dam transverse joints

(19)

式中:初始接触距离c取100 mm,初始压力p0取0.3 MPa。

每个横缝接触面包含主面和从面,主从面的接触关系如图8所示。

图8 主从面接触模型Fig.8 Dam transverse joint model

3.1.4黏弹性人工边界

黏弹性人工边界可以较为合理地模拟地基与结构的相互作用关系,具有多方面优点,不仅可以吸收边界辐射能量,同时可以反映地基的弹性恢复能力[17],其如图9所示。

图9 黏弹性人工边界模型Fig.9 Viscoelastic artificial boundary model

人工边界上,节点上的参数由式(20)计算给出:

(20)

式中:Kn为弹簧-阻尼元件法向刚度系数;Cn为弹簧-阻尼元件法向阻尼系数;Ks为弹簧-阻尼元件切向刚度系数;Cs弹簧-阻尼元件切向阻尼系数;Cp,Cs分别为P波和S波波速;r为结构几何中心到人工边界点的等效距离;ρ为介质密度;G为剪切模量;λ为第一拉梅常数;a,b为参数。

3.2 ETA地震动时程

根据选取的特征周期生成标准设计反应谱作为目标谱,通过与0～5 s、0～10 s、0～15 s和0～20 s等四个不同时间段内目标反应谱进行拟合,生成3条ETA人工地震动时程,作为3个方向的输入地震动。图10给出了Tg=0.322 s时,3个方向的ETA加速度时程,随着时间的增加,其加速度幅值不断增大。图2显示加速度时程对应的加速度反应谱在标准反应谱周围振荡,拟合效果良好。

图10 ETA时程Fig.10 ETA time history

值得注意的是,图10所示的ETA时程在不同时间点对应不用的峰值加速度,进行一次时程分析便可以得到不同峰值加速度下的动力响应,相较IDA分析方法,无需进行大量调幅计算,显著降低了计算量,适用于对大型结构的分析。

3.3 动力响应分析

3.3.1响应分析

通过ABAQUS软件计算后,发现坝体下游面损伤主要集中的中上部,部分工况下会形成贯穿性损伤,上游面损伤主要集中在拱端,损伤符合一般规律。考虑坝体材料的随机性后,损伤的分布更加离散,且坝体的损伤情况主要会受到坝体材料参数均值和相关长度的影响。图11(a)展示了3个具有相近均值,但相关长度不同的3个拱坝有限元模型的材料参数分布情况,图11(b)展示了峰值加速度为0.6g时,以上3种工况的坝体损伤情况。通过对比,在坝体材料均值接近的情况下,相关长度越小,坝体的损伤越严重,随着相关长度的增大,坝体的抗震能力逐渐提高。同样,图12对比了相关长度增大,但均值逐渐减小的3组工况,经对比,虽然相关长度不断增大,但随着材料参数的下降,坝体的破坏逐渐加重,说明坝体材料参数的均值也是决定坝体抗震能力的关键因素。

图11 参数均值相近,相关长度不同工况的计算结果Fig.11 Results under similar mean parameters and different correlation lengths

图12 相关长度增大,参数均值减小工况的计算结果Fig.12 Results under decreasing mean parameters and increasing correlation lengths

3.3.2概率密度分析

采用横缝开度、损伤体积比和相对位移作为损伤评价指标。计算结果显示,随着地震动峰值加速度(PGA)值增大,各指标呈现明显的增大趋势,不同工况的分析结果离散性较高,考虑为坝体材料的随机性导致。采用PDEM法和DPIM法分别计算了各指标概率演化过程。图13展示了随PGA值的增大,不同指标的变化趋势。各指标均随时间的推移逐渐增大,且分布范围逐渐变大。图14展示了横缝开度的概率分布情况。从等高线图中可以发现,各时刻的概率分布近似服从正态分布,均值约等于该时刻的平均横缝开度,且PGA越大,方差越大,概率分布曲线由“高瘦型”向“矮胖型”转变,说明PGA越大,离散型越高。

图13 各评价指标随PGA变化曲线Fig.13 Curve of evaluation index with the change of PGA

图14 横缝开度概率密度演化过程Fig.14 Evolution process of probability density of transverse joint openness

在PGA较小时,各指标概率分布范围较为集中,随着PGA的增大,分布范围趋于离散。两种概率密度方法计算结果高度一致,但DPIM法的平均计算速度是PDEM法的20倍左右,优势明显。

为对比直接概率积分法和概率密度演化方法的计算结果是否吻合,取10 s、15 s和20 s为特定时刻,对各特定时刻的概率密度分布情况以及累计概率分布情况进行对比。结果显示各指标使用两种方法计算的结果差距较小,DPIM法计算结果在局部略有振荡,但幅度较小。图15展示了以坝体横缝开度为指标的计算对比结果,各时刻均较为接近。

图15 概率密度分布及累计概率密度在特定时刻计算结果对比Fig.15 Comparison of probability density and accumulated probability density at specific time

3.3.3失效概率分析

失效概率是表征结构破坏程度的有效指标,当结构参数达到阈值后,认为结构发生破坏,反之,认为结构是安全的。首先采用概率密度方法计算结构的失效概率[18]。根据3.3.2节中所求概率密度曲线,通过式(21)积分后求得失效概率:

(21)

式中:XL为参数的阈值,阈值的选取可以参照文献[19],计算结果如图16所示。通过对比,直接概率积分法与概率密度演化法的计算结果较为吻合,再次验证了该方法的正确性。各项指标计算的失效概率快速增大的区域均集中在0.2g～0.4g之间,且互相之间较为吻合,具有一定等效关系。

图16 评价指标失效概率曲线Fig.16 Failure probability curve of evaluation index

通过对比两种方法计算的各损伤评价指标失效概率曲线发现,以损伤体积比为评价指标时,两种方法的计算结果在初期有所差异;随着PGA的增大,基于相对位移指标计算的失效概率曲线增速较其他两种指标计算的失效概率曲线更加平缓,说明不同损伤评价指标具有一定离散性;使用横缝开度作为评价指标计算的结果较为合理。

4 结 论

本文引入随机场理论,模拟拱坝混凝土材料参数的随机性。为避免偶然性,计算多组样本,并基于维诺图赋以相应概率。基于耐震时程法(ETA)理论,根据不同特征周期对应的标准谱生成ETA人工地震动时程。通过对拱坝有限元模型的分析研究了抗拉强度和断裂能随机分布对白鹤滩拱坝损伤分布的影响。将各工况结果乘以相应赋得概率并求和后,作为最终结果,并运用概率密度方法建立了各指标随时间变化的时变概率信息,得到其概率密度演化曲线,曲线的走向与参数均值保持一致。设定损伤阈值,基于概率密度曲线对超过阈值的部分进行积分得到失效概率,建立了峰值加速度与失效概率的关系。

研究结果表明:考虑混凝土抗拉强度和断裂能的随机性后,白鹤滩拱坝动力响应具有明显的离散性,不同工况的分析结果随着PGA的增大其差异逐渐增大;分析可知,当坝体材料随机参数均值相近时,随着相关长度的增大,坝体的损伤情况呈减小趋势,即坝体材料参数越趋向于均值,坝体抗震性能越好;而当坝体材料随机参数均值不同时,即使相关长度较大,参数均值较小的破坏依然严重。经对比,运用横缝开度和损伤体积比作为损伤评价指标计算结构的失效概率时,其计算结果更为接近。