混凝土拱坝地震非线性响应中的应变率效应研究

董 雨,王国盛,王进廷,杜修力

(1.北京工业大学 岩土与城市地下研究所,北京 100124;2.清华大学 土木水利学院,北京 100084)

拱坝是一种应用广泛的坝型,其坝型结构具有优良的力学性能,易解决枢纽泄洪,且容易施工导流。与同级别重力坝相比较,能够显著节省筑坝材料,提升经济效益[1-3]。混凝土作为拱坝主要的建筑材料,是一种典型的率相关材料,其强度受应变率影响十分显著。对于混凝土应变率效应的研究,最早可以追溯到1917年Abram[4]所进行的混凝土动态单轴抗压强度的研究,随后的一百多年里,学者进行了大量关于混凝土率相关强度的研究[5-19]。然而,当前在大坝抗震设计中,通常不考虑混凝土材料的率效应,或者只是在静强度的基础上考虑强度的整体放大。例如,应用比较广泛的是Rapheal[20]利用试验给出混凝土大坝在地震作用下动强度的建议值,其中混凝土动态抗压强度较静强度提高31%,动态直接拉伸强度较静态直拉强度提高66%,动态劈拉强度提高45%。上述建议值是在确定加载条件下获得的,即应变速率近似等于5 Hz振动产生的加载速率,但目前已被不分情况地普遍推广应用于大坝的设计[21]。如2015年NB 35047—2015《水电工程水工建筑物抗震设计规范》[22]中规定,混凝土动态强度标准值可较其静态标准值提高20%,动态弹性模量标准值可较其静态标准值可提高50%,动态抗拉强度标准值取为动态抗压强度标准值的10%。

在实际地震作用下混凝土拱坝的反应中,地震强度、坝体尺度和形式对混凝土的率相关特性均存在较大影响,并且在不同时刻、坝体不同部位混凝土的率效应差别也很大,尤其存在非线性和率效应的耦合影响。混凝土的动态强度若在静强度基础上提高20%,无法考虑大坝中混凝土强度在不同时刻和不同部位的变化,原因是没有采用能够反映率效应的混凝土本构模型。针对地震循环作用下混凝土材料的损伤特性,Lee等[23]建立了混凝土材料的塑性损伤本构模型,利用该模型开展了Koyna大坝的地震响应分析。陈健云等[24]通过在损伤张量中考虑应变率效应,建立了混凝土材料率相关损伤模型,并开展了300 m级高拱坝地震响应分析,获得了拱坝混凝土应变率和拉压损伤的分布规律。肖诗云等[25]将混凝土静态(Hsieh-Ting-Chen,HTC)模型发展为一致黏塑性模型,对某高拱坝进行了地震反应分析,结果表明应变率对拱坝的应力和塑性应变均有较大影响。其他学者[26-28]相继发展了混凝土材料率相关本构模型,并开展了高拱坝的地震反应分析。上述研究表明,由于地震荷载在坝体中产生的应变率分布不同,坝体各部位混凝土动态力学性能也会不同。在采用不同的混凝土损伤本构模型时,坝体的损伤分布有所区别,同时大坝的破坏模式与选取的地震波密切相关。在前人所做的数值模拟中,通过采用率相关本构模型开展研究,对考虑应变率和不考虑率影响下坝体特征点的加速度、位移值、坝体拉应力分布和损伤进行了对比分析。从应变率效应产生的本质上来看,地震荷载在坝体中产生的应变率随时间和空间不断变化,变化的应变率会引起混凝土动态强度变化,最终导致坝体的加速度、位移、应力分布情况、损伤值等是连续变化的过程,而非简单的一个状态。

本文采用有效硬化函数和动态S型强度准则对混凝土塑性损伤模型进行改进,利用S-CDP模型开展了大岗山拱坝有限元模拟。从坝体应变率的角度进行研究,选取坝体上具有代表性的一些单元,将特征单元的应变以及应变率随地震时间的变化曲线提取出来,通过直接对不同特征单元的应变率曲线进行分析,得到坝体不同部位应变率随时间的变化过程。再从宏观层面上对由于动强度的不同引起的损伤情况进行研究,分析地震作用下混凝土拱坝的应变率分布、相应的动态增长因子(dynamic increase factors,DIF)的变化情况以及应变率效应带来的影响。

1 混凝土损伤塑性模型的改进

ABAQUS有限元软件中内嵌的混凝土损伤塑性(concrete damage plasticity,CDP)模型[29]常被作为混凝土的材料模型使用,原始CDP模型需要用户自定义动态强度随应变率变化规律,然后将单轴拉压应力-应变曲线作为非弹性应变率的函数来考虑应变率效应。动强度定义不合理很容易导致数值计算不收敛。针对这一问题,本文对原始CDP模型进行了改进。首先采用在有效应力空间中确定的多轴硬化函数代替名义应力空间中的单轴硬化函数,保证三维塑性变形的合理计算;其次选用连续光滑、能反映真实动强度的动态强度S准则[30]考虑混凝土的应变率效应,在子程序中使用S准则来计算应变率相关应力-应变关系,提高计算的收敛性,同时避免惯性效应的多次考虑,从而减小计算结果与实际结果的偏差。具体的改进方法阐述如下。

原始CDP模型为双标量塑性损伤本构模型,一个损伤变量用于拉伸损伤,一个损伤变量用于压缩损伤。通过损伤变量建立了名义应力与有效应力的关系

(1)

(2)

式中：E0为初始无损伤弹性刚度张量;ε为总应变张量;εp为塑性应变张量。

式(1)中ω为包含压缩损伤和拉伸损伤的损伤变量

ω=1-(1-ωt)(1-ωc)

(3)

式中：ωt为拉伸损伤变量;ωc为压缩损伤变量。

有效应力表示的屈服函数是得到塑性演化的关键,CDP模型的屈服函数表达式如下

cc(κ)

(4)

(5)

在屈服函数式(4)中,cc(κ)和β(κ)为两个硬化函数。其中cc(κ)为压缩硬化,β(κ)描述压缩和拉伸的协同硬化,其表达式为

(6)

式中：ct(κ)为拉伸硬化函数;κ为硬化参数。

1.1 有效硬化函数

(7)

采用多轴应力条件下的统一硬化/软化函数[33]来描述名义应力空间中的硬化函数,它能够考虑应力状态对硬化/软化行为的影响,其方程如下

(8)

式中：A和D为形状参数;xi=εpd/εpd0为硬化参数;εpd为等效塑性剪切应变;εpd0为峰值应力下εpd的值。εpd的表达方式为

(9)

式中,εpi为塑性主应变,本文规定拉为正,压为负。

将式(8)代入式(7),得

(10)

式(10)为多轴条件下的有效硬化函数,对于单轴拉伸条件,式(9)中的εp2=εp3=0,εp1>0;对于单轴压缩条件,式(9)中的εp1=εp2=0,εp3<0。因此,将上述应变条件代入式(10),多轴硬化规律可退化为单轴拉伸和单轴压缩硬化规律,进而可以确定屈服函数式(4)中的cc(κ)与β(κ)。当CDP模型采用由有效应力空间中确定的硬化规律时,S-CDP模型能够更加合理地描述混凝土材料在有效应力空间中的屈服面演化规律。

1.2 应变率相关规律

对于混凝土材料的动态单轴强度准则,最具代表性的准则为CEB-FIP规范建议的双线性单轴动态强度准则[34],该准则的方程是分段函数,对于分段表达形式的动态单轴强度准则,准则中存在拐点,拐点位置的确定缺乏理论依据。当前大多数既有的动态强度准则都是发散的,无法反映混凝土的极限动强度,并且也没有区分真实动强度与宏观动态承载力,导致计算结果过高估计混凝土材料的真实动强度[35-36]。在本文将引入非线性动态单轴S准则,该准则是基于混凝土材料的动态物理机制建立的,能够合理反映混凝土材料的真实动强度和试验表现出的极限动强度,S准则的材料参数较少且具有明确的物理意义,其强度曲线连续光滑,不存在明显的拐点。S准则的方程如下

(11)

2 大岗山拱坝数值模型及地震作用下的损伤规律

2.1 大岗山拱坝模型参数与建模过程

大岗山拱坝坝型为混凝土双曲拱坝,最大坝高为210 m,坝顶高程为1 130 m[38]。大岗山坝体-地基有限元模型如图1所示。数值模型设置了与实际工程相同的28条坝体横缝[39],在进行分析的过程中,考虑了键槽咬合作用的高拱坝非线性行为[40]。拱坝横缝在地震荷载作用下张开、闭合的非线性力学行为引入了动接触边界以及限制接触面切向滑移的弹簧单元[41]。并采用黏弹性人工边界模拟无限地基辐射阻尼效应[42]。

(a) 坝体-地基

参照同类工程经验,本文的计算模型选取的地基范围为：下游面地基2倍坝高,上游面地基、左右岸及自坝底向下均为1.5倍坝高。有限元模型中,采用S-CDP模型作为混凝土的材料模型;坝体混凝土单元类型采用八节点六面体线性缩减积分实体单元(C3D8R),模型总的单元数量为121 122个,其中坝体单元数为79 638个,地基单元为41 484个,节点数为148 551。

分别采用“Static”分析步计算初始静态荷载的作用,采用“Dynamic/Implicit”分析步计算地震荷载的作用,初始增量步时长设置为0.01 s。坝体混凝土的材料参数为：混凝土密度2 400 kg/m3,弹性模量较其静态标准值提高50%取为36.65 GPa,泊松比0.17,膨胀角45°,α=0.14,形状参数A=2.5,D=2.0。S准则中强度参数的取值：动态拉伸条件下ξ=1,Fmax=10.2,u0=-1,F0=0.5;动态压缩条件下ξ=1.2,Fmax=3.72,u0=1.95,F0=0.5。基岩假定由两部分构成：远离坝体的线弹性基岩与靠近坝体的非线性基岩。线弹性基岩与非线性基岩相同的材料参数为：基岩密度2 650kg/m3,弹性模量20 GPa,泊松比0.25。非线弹性基岩采用Drucker Prager Harding模型,摩擦角65.1°,流变应力比取1,膨胀角取65.1°。作用荷载包括坝体自质量、上游库水静水压力、温度荷载、上游库水动水压力和地震荷载。其中,动水压力采用Westergaard附加质量法考虑[43]。

大岗山拱坝所处地区地震基本烈度为Ⅷ度,设防标准按100 a设计基准期超越概率为2%(重现期5 000 a)的基岩水平加速度,其设计加速度水平向为0.557 5g,竖向加速度取水平方向的2/3[44]。为突显应变率对拱坝地震非线性响应的影响,本文考虑地震动的不确定性[45],从实测地震数据库中选取3组地震动时程曲线：第一组为等比例放大1.271倍得到峰值地面加速度(peak ground acceleration,PGA) 0.557g的地震动时程曲线。第二组等比例放大1.821倍得到PGA=0.663g地震动时程曲线;第三组为等比例放大1.025倍得到PGA=0.836g的地震动时程曲线;XYZ3个方向的地震波加速度时程曲线如图2所示。

(a) PGA=0.557g的地震动时程曲线

2.2 3组地震作用下坝体的损伤规律

利用S-CDP模型作为混凝土的材料模型,分别计算了在第一组、第二组、第三组地震动荷载作用下坝体的动力响应,其损伤分布如图3～图5所示。由图3～图5可知,坝基面的损伤比较严重,坝体下游面出现了大面积的损伤,损伤结果与其他学者[46]论文中的强震下高拱坝的损伤开裂结果具有一定的相似性,说明了本文计算结果的合理性。在第一组地震的损伤结果中,上游面靠近坝顶处也出现了小范围的损伤。两组地震作用下坝体损伤最严重的部位均出现在坝基,原因是坝基部位的应力集中所导致。

(a)下游面

(a)下游面

(a)下游面

在3组地震作用下得到的绝对值最大的主应变如图6所示,将以符号εmpa来表示绝对值最大的主应变。选取7个特征单元,编号为：L1,L2,R1,R2,F1,F2,M1。

(a) 第一组地震动下的分布情况

坝体的εmpa主要分布在坝踵部位及坝体下游面的中部,结合损伤分布的情况来看,坝体在承受地震荷载作用时,最容易出现损伤的部位为坝基以及坝体下游面的中部区域。

第一组地震动、第二组地震动和第三组地震动荷载作用下坝体7个特征单元的εmpa时程曲线如图7所示。图7表明,选择的特征单元的拉应变要大于压应变,而坝体混凝土的抗拉强度低,因此出现的损伤都是拉伸损伤,没有压缩损伤。

(a) 第一组地震动

对图7的εmpa的时程曲线进行分析,第一组地震动下坝基处两个特征单元F1和F2的εmpa最大分别达到了0.004 5和0.006 5。除坝基外坝体的其他特征单元的εmpa远小于特征单元F1和F2的值,都分布在0.000 5以下。第二组地震动荷载作用下,坝基处的两个特征单元F1和F2的εmpa最大分别达到了0.007和0.007 5。除坝基外坝体的其它特征单元的εmpa远小于F1和F2单元的值,都分布在0.002以下。在第三组地震动荷载作用下,坝基处的两个特征单元F1和F2的εmpa最大分别达到了0.002 2和0.006,除坝基外坝体的其他特征单元的εmpa远小于特征单元F1和F2的值,都分布在0.001以下。

3 混凝土拱坝的应变率效应分析

为研究在地震荷载作用下混凝土拱坝的应变率效应,采用S-CDP模型,对大岗山混凝土拱坝进行了第一组地震动、第二组地震动与第三组地震动下的非线性地震响应分析,得到其应变率变化规律和相应的DIF变化规律。

3.1 第一组地震动作用下的应变率规律

基于DIF的定义(动强度与静强度的比值),通过应变率时程曲线获得第一组地震动作用下特征单元的DIF时程曲线,如图8所示。

(a) 拉伸动态增长因子

在第一组地震动作用的情况下,拉伸DIF变化范围为1.00～1.21,压缩DIF变化范围为1.000 0～1.003 3。混凝土的应变率效应增大了坝体动强度,以拉应变产生的动强度增长为主,DIF的数值分布表明动强度增长范围大致分布在0%～21%。

3.2 第二组地震动作用下的应变率规律

对图7中εmpa的时程曲线微分,得到特征单元拉应变率和压应变率随时间的变化规律如图9所示。图9(a)表明：在整个地震过程中拉应变率幅值最小的是R1单元的曲线,该单元位于坝体的中部偏右岸位置;拉应变率幅值最大的是F2单元的曲线,该单元位于坝体坝基部位。拉应变率最大峰值出现的时间段分布在地震荷载作用的5～7 s,该时间段也是地震加速度时程曲线最大峰值出现的区间。图9(b)表明：压应变率幅值最小的是R2单元的曲线,该单元位于坝体的中间部位;压应变率幅值最大的为F2单元的曲线。

(a) 拉伸应变率

分析图9中应变率的变化规律,坝体出现拉压应变率幅值最大的位置与损伤最大的位置是一致的。然而,拉应变率与压应变率幅值最小的单元则出现在坝体的不同部位。此外,坝体不同位置处单元应变率的幅值和出现幅值的时刻都是不同的。通过采用S-CDP模型,计算得到坝体多处单元的拉压应变率曲线,直观地展示了坝体不同部位混凝土应变率在时间和空间上的连续变化过程。

基于DIF的定义(动强度与静强度的比值),计算得到所选单元的DIF随时间的变化规律,如图10所示。

(a) 拉伸动态增长因子

在第二组地震动作用的情况下,坝体的拉应变率变化范围为0～1.6×10-2s-1,最大拉应变率只在6.2 s时出现一次,其余时刻坝体混凝土的拉应变率基本上都位于0.8×10-2s-1以下。坝体的压应变率范围为0～0.25×10-2s-1,其余时刻坝体混凝土的压应变率基本上都位于0.15×10-2s-1以下。在图10(a)中,拉伸DIF幅值最小的是R2单元的曲线,拉伸DIF幅值最大的是F2单元的曲线。在图10(b)中,压缩DIF幅值最小的是R2单元的曲线,压缩DIF幅值最大的是F2单元的曲线。拉伸DIF变化范围为1.00～1.23。压缩DIF变化范围为1.000 0～1.002 5。拉伸DIF的最大值远大于压缩DIF的最大值,但是应指出拉伸DIF最大值1.23只在6.2 s时出现了一次,出现的部位为坝基处的F2单元。

综上数据,采用S-CDP模型能够很好地考虑不同应变率下混凝土动强度的变化,且不同应变率下的动强度和已有的工程经验是一致的,证明了模型的合理性。在第二组地震动荷载作用下,混凝土的应变率效应对坝体动强度增长起到了一定的作用,以拉应变产生的动强度增长为主,DIF的数值分布表明动强度增长范围大致分布在0%～23%。

3.3 第三组地震动作用下的应变率规律

采用与3.2节相同的方法,通过应变率时程曲线获得第三组地震动作用下特征单元的DIF时程曲线,如图11所示。

在第三组地震动作用的情况下,坝体的拉应变率变化范围为0～1.1×10-2s-1,坝体的压应变率范围为0～0.27×10-2s-1,在整个地震动荷载作用的时间内,坝体混凝土的拉应变率基本上都位于0.4×10-2s-1以下,坝体混凝土的压应变率基本上都位于0.2×10-2s-1以下。在图11(a)中,拉伸DIF幅值最小的是R1单元的曲线,拉伸DIF幅值最大的是F2单元的曲线。在图11(b)中,压缩DIF幅值最小的是M1单元的曲线,压缩DIF幅值最大的是F2单元的曲线。拉伸DIF变化范围为1.00～1.19,压缩DIF变化范围为1.000 0～1.002 5。拉伸DIF的最大值远大于压缩DIF的最大值,但是最大的拉伸DIF值1.19只在5.4 s时出现了一次,出现的部位为位于坝基处的F2单元。

综上数据,在第三组地震动作用下,混凝土的应变率效应增大了坝体动强度,以拉应变产生的动强度增长为主,DIF的数值分布表明动强度增长范围大致分布在0%～19%。

4 率相关强度对坝体地震响应的影响

为了分析率相关强度对坝体地震响应的影响,本章考虑在第二组和第三组地震动荷载作用,分析3种工况：工况I,采用率无关CDP模型,按混凝土材料静态强度,弹性模量较其静态标准值提高50%计算的拱坝非线性响应;工况II,按照现行的NB 35047—2015《水电工程水工建筑物抗震设计规范》,直接将坝体混凝土材料静态强度提高20%,动态弹性模量标准值较其静态标准值提高50%,采用率无关CDP模型计算拱坝非线性响应;工况III,采用S-CDP模型作为混凝土材料模型,动态弹性模量标准值较其静态标准值提高50%,计算拱坝的动力非线性响应。3种不同的工况,均同时考虑坝体自体质量、静水压力、动水压力、温度荷载和地震荷载作用,第二组地震动荷载坝体损伤情况如图12所示,第三组地震动荷载坝体损伤情况如图13所示。

(a) 工况Ⅰ

(a) 工况Ⅰ

图12和图13中分别展示了第二组和第三组地震动作用下的3种工况坝体下游面的损伤情况。从图12、图13可知,两组地震动荷载的3种工况坝体下游面都出现了损伤情况但均未形成自下游坝面至上游坝面的贯穿损伤。对比工况I、工况II、工况III的损伤响应,可以发现,在同一地震动荷载作用下,工况Ⅰ和工况Ⅲ坝体的损伤区域基本一致,都出现在坝基以及坝体的中上部位,但是工况Ⅲ中的坝体损伤程度比工况Ⅰ的损伤程度小,这说明在地震荷载的作用下,坝体混凝土采用S-CDP模型计算得到的损伤比采用率无关CDP模型计算得到的损伤会减小,率相关强度对最终的损伤结果产生了影响。

对图12(b)和图12(c)的损伤云图及损伤因子进行对比,对图13(b)和图13(c)的损伤云图及损伤因子进行对比,可以看出工况Ⅱ的坝体损伤区域与损伤因子较工况Ⅲ整体偏小,但是在实际工程中混凝土拱坝的地震响应存在应变率效应,其坝体的损伤区域及损伤因子与在静强度的基础上考虑强度的整体放大20%的情况相比时会更大,这说明采用S-CDP模型时混凝土的动强度增长无法一直达到20%的程度,也说明在进行水工建筑物抗震计算时,若将整个坝体的混凝土的动强度提高20%,可能会造成损伤结果与实际情况有所偏差。

5 结 论

通过引入有效硬化函数和S型率相关强度准则,考虑混凝土材料的三维塑性变形行为和应变率效应,对混凝土塑性损伤本构模型进行了改进。使用S-CDP模型,对大岗山混凝土拱坝进行了非线性地震响应分析。分别计算了PGA=0.557g,PGA=0.663g,PGA=0.836g的地震动荷载作用下拱坝的动力响应,分析了坝体的应变分布规律、应变率的分布规律及相应的DIF变化规律、证明了S-CDP模型的合理性。主要结论如下：

(1) 在地震荷载作用下,拱坝混凝土的压应变率效应不明显;然而,坝基和部分坝体区域产生了较大的塑性变形,拱坝混凝土的拉应变率效应明显。此外,通过分析坝体中特征单元的应变以及应变率随时间的动态变化规律,发现采用S-CDP模型建立的数值模型能够合理反映混凝土的应变率效应。

(2) 通过在本构模型中考虑混凝土率效应,计算得到的大岗山拱坝坝体的拉应变率最大峰值为1.6×10-2s-1,相对应的拉伸DIF为1.23,仅在坝基处的某瞬间出现了一次,除此外拉应变率变化较大,并不能时刻保证坝体混凝土的动强度提高20%及以上,地震荷载作用下坝体混凝土的压应变率效应不显著。在水工建筑物抗震设计时,若将混凝土的动强度提高固定的20%,可能会造成计算的损伤结果与实际的损伤有所偏差。动强度提高固定的20%计算得到的损伤结果相较于本文中采用S-CDP模型计算得到的损伤结果偏大。因此在水工建筑物抗震设计时,推荐选用合适的率相关本构模型进行计算,以尽可能的保证计算的精确度。