工程科学与技术   2020, Vol. 52 Issue (4): 117-123
循环荷载下考虑滞回效应的混凝土损伤模型
刘智1, 赵兰浩2, 吴晓彬1, 胡国平1, 周清勇1     
1. 江西省水利科学研究院,江西 南昌 330029;
2. 河海大学 水利水电学院,江苏 南京 210098
基金项目: 国家重点研发计划项目(2016YFC0401601;2018YFC0406705);江西省水利厅科技计划项目(202022YBKT04)
摘要: 随着中国二次能源结构逐渐优化,大批高混凝土坝工程在西南地区应运而生,其抗震安全问题日益突出。针对高混凝土坝地震过程中复杂材料非线性问题,建立循环荷载下考虑滞回效应的混凝土损伤模型。模型分别选取符合混凝土实际变形特性的拉伸与压缩骨架线以考虑材料拉压异性,骨架线中含有软化段系数以适应试验结果的离散性。采用不依赖于骨架线形状的滞回效应加卸载特征点表达式,搭建能够反映混凝土循环荷载作用下软化段滞回效应的卸载路径与重新加载路径,并建议设立残余应变临界值解决残余塑性应变与卸载应变比值随应变增长的持续发散问题。模型将复杂多轴问题转化至单轴等效应变空间中求解,计算参数少,数学表达式简单,并通过对比混凝土循环拉伸荷载试验和Koyna重力坝的震害模拟验证了模型在非线性问题求解上的正确性与高效性。研究为高混凝土坝的抗震分析及安全评价提供理论支撑。
关键词: 损伤模型    循环荷载    滞回效应    抗震分析    
Damage Model of Concrete Considering Hysteretic Effect Under Cyclic Loading
LIU Zhi1, ZHAO Lanhao2, WU Xiaobin1, HU Guoping1, ZHOU Qingyong1     
1. Jiangxi Provincial Inst. of Water Sci., Nanchang 330029, China;
2. College of Water Conservancy and Hydropower Eng., Hohai Univ., Nanjing 210098, China
Abstract: With the gradual optimization of China secondary energy structure, a large number of high concrete dam projects have emerged in the southwestern region, and their seismic safety issues have become increasingly prominent. A concrete damage model considering hysteresis effect under cyclic loading was developed to solve the complicated material non-linear problem during high concrete dam earthquake. The tension and compression skeleton curves that conform to the actual deformation characteristics of concrete were used in the proposed model for adapting to the anisotropy of material tension and compression. The skeleton curves contained the softening section coefficient to adapt to the dispersion of test results. The expression of characteristic points of loading and unloading hysteresis effect that did not depend on the shape of the skeleton curve was used to establish an unloading path and a reloading path that reflected the hysteresis effect of softening section under the cyclic loading of concrete. A critical value of residual strain was suggested to solve the continuous divergence problem, that the ratio of residual plastic strain to unloading strain increased with the unloading strain. The model transmitted the complex multi-axis problem to a uniaxial equivalent strain space for problem solving, with fewer calculation parameters and simple mathematical expressions. The reliability and efficiency of the proposed model were investigated for nonlinear problems solving by comparing the concrete cyclic tensile load test with the earthquake damage simulation of Koyna gravity dam. This study can provide effective theoretical support for seismic analysis and safety evaluation of high concrete dams.
Key words: damage model    circle loading    hysteresis effect    anti-seismic analysis    

随着中国西南强震地区水利工程建设如火如荼,地震荷载作为非常规荷载,在工程设计及结构安全分析中与常规荷载一样占据了十分重要的地位。非常规荷载具有强大的破坏性及独特的不可预知性,特别是其循环往复的荷载形式,使得混凝土材料的非线性特性愈发强烈。循环荷载作用下混凝土材料处于不断卸载与重新加载的循环过程中,将形成明显的滞回效应,期间的损伤积累及能量耗散所引起的刚度退化与应力重分布势必将影响材料后续非线性性能的变化过程。

目前针对混凝土循环荷载下滞回规则的研究较少。有学者依据试验数据推导数学公式,对受压循环荷载下的滞回行为进行简化的模型构建[1];有学者研究受拉循环荷载下的滞回规则[2-4]。Konstantinidis等[5]对当前混凝土循环受压荷载作用下的本构模型进行统计;Aslani[6]、过镇海[7]等对混凝土滞回行为的特征进行了总结。传统损伤模型通常将滞回效应简化为线性表达,这种通过线性形式描述非线性现象的假定无法完整反映循环荷载作用下的损伤积累过程,滞后的损伤积累无法及时反馈劣化部位的应力转移,势必将影响后续的仿真结果。李正等[8]将《混凝土结构设计规范》[9]中单轴应力应变骨架线与Yassin提出的滞回模型相结合,建立了循环荷载下用于非线性分析的滞回本构模型。但模型中滞回规则的重新加载曲线将返回卸载点,能够表现滞回期间的刚度退化,但无法准确描述其损伤积累过程,针对实际工程的应用性有待进一步完善。

在已有混凝土四参数损伤模型的基础上,结合滞回规则中的加卸载特征点及加卸载路径,构建循环荷载下考虑滞回效应的四参数损伤模型。模型包含了循环荷载下混凝土演化过程中的拉压异性、刚度退化、强度软化、塑性不可恢复变形及滞回效应等复杂非线性特性,并通过混凝土单轴循环荷载试验与Koyna重力坝的震害模拟验证了模型在非线性问题求解上的正确性。

1 混凝土四参数损伤模型

基于Hsieh–Ting–Chen的应力空间四参数破坏准则,韦未等[10]建立了一种基于应变空间的四参数破坏准则:

$ F\left( {I_1',J_2',{{\rm{\varepsilon}} _0}} \right) = A\frac{{J_2'}}{{{{\rm{\varepsilon}} _0}}} + B\sqrt {J_2'} + C{{\rm{\varepsilon}} _1} + DI_1' = 0 $ (1)

式中: $I_1' = {{\rm{\varepsilon}} _{ii}}(i = 1,2,3)$ 为应变张量第1不变量; $J_2' = $ $ {e_{ij}}{e_{ij}}/2(i,j = 1,2,3)$ ,为应变偏量第2不变量; ${{\rm{\varepsilon}} _1} = $ $ \dfrac{2}{{\sqrt 3 }}\sqrt {J_2'} \sin ({\rm{\theta}} + \dfrac{2}{3}\text{π} ) + \dfrac{1}{3}I_1'$ ,为最大主应变, ${{\theta}} = \dfrac{1}{3}\arcsin $ $\bigg( - \dfrac{{3\sqrt 3 J_3'}}{{2\sqrt {J_2^{'3}} }}\bigg){\text{,}}\left| \theta \right| \le {60^ \circ }$ $J_3' = {e_{ij}}{e_{jk}}{e_{ki}} $ $(i,j,k = 1,2,3)$ ,为应变偏量第3不变量; ${{\rm{\varepsilon}} _0} = \dfrac{{{f_{\rm{t}}}}}{{E}}$ ,为材料峰值应变。

韦未等[11]基于四参数破坏准则的基本思路,提出了一种新的四参数等效应变计算方法。假定四参数破坏准则在应变软化段内仍然适用,且ABCD4个参数保持不变,其形式与式(1)相同,式中 ${\varepsilon _0}$ 被等效应变 ${\varepsilon ^*}$ 替代:

$ {{\rm{\varepsilon}} ^ * } = A\frac{{J_2'}}{{{{\rm{\varepsilon}} ^ * }}} + B\sqrt {J_2'} + C{{\rm{\varepsilon}} _1} + DI_1' $ (2)

式中: ${{\rm{\varepsilon}} _1} = \dfrac{2}{{\sqrt 3 }}\sqrt {J_2'} \sin \bigg({{\theta}} + \dfrac{2}{3}\text{π} \bigg) + \dfrac{1}{3}I_1'$ ,为最大主应变; $I_1' = ({{\rm{\varepsilon}} _1} + {{\rm{\varepsilon}} _2} + {{\rm{\varepsilon}} _3})/3$ 为应变张量第1不变量; $J_2' = \dfrac{1}{2} $ $\left[ {{{\left( {{{\rm{\varepsilon}} _1} - {{\rm{\varepsilon}} _m}} \right)}^2} + {{\left( {{{\rm{\varepsilon}} _2} - {{\rm{\varepsilon}} _m}} \right)}^2} + {{\left( {{{\rm{\varepsilon}} _3} - {{\rm{\varepsilon}} _m}} \right)}^2}} \right]$ 为应变偏量第2不变量; ${{\theta}} = \dfrac{1}{3}\arcsin $ $\bigg( - \dfrac{{3\sqrt 3 J_3'}}{{2\sqrt {J_2^{'3}} }}\bigg){\text{,}}\left| \theta \right| \le {60^ \circ }$ $J_3' = {{\rm{\varepsilon}} _1}{{\rm{\varepsilon}} _2}{{\rm{\varepsilon}} _3}$ 为应变偏量第3不变量, ${{\rm{\varepsilon}} _1}{\text{、}}{{\rm{\varepsilon}} _2}{\text{、}}{{\rm{\varepsilon}} _3}$ 分别为三向主应变。

ABCD4个参数与破坏准则使用参数相同,求解式(2),且考虑到 ${{{\varepsilon}} ^ * } \ge 0$ ,可得多轴应力状态下等效应变 ${{\rm{\varepsilon}} ^ * }$

$ {{\rm{\varepsilon}} ^ * } = \frac{{(B\sqrt {J_2'} + C{{\rm{\varepsilon}} _1} + DI_1') + \sqrt {{{(B\sqrt {J_2'} + C{{\rm{\varepsilon}} _1} + DI_1')}^{\rm{2}}}{\rm{ + 4}}AJ_{\rm{2}}'} }}{2} $ (3)

式(3)形式简单明了,可将真实空间中复杂的多轴问题转化为等效空间中简单的单轴问题。模型已得到充分的理论验证及广泛的工程应用[12]

根据混凝土的拉压异性,四参数混凝土损伤模型骨架线将按照受拉与受压两种状态分别选取。损伤模型选取过镇海[7]提出的应力应变全曲线,该曲线已得到国内外科研工作者的认可,并纳入中国混凝土结构设计规范[9]

混凝土单轴受拉应力应变曲线表达式如下:

$ {\rm{\sigma}} = \left( {1 - {d_{\rm{t}}}} \right){E_{\rm{c}}}{\rm{\varepsilon}} $ (4)
$ {d_{\rm{t}}} = \left\{\!\!\!\! {\begin{array}{*{20}{l}} {1 - \dfrac{{{f_{{\rm{t}},{\rm{r}}}}}}{{{E_{\rm{c}}}{{\rm{\varepsilon }}_{t,{\rm{r}}}}}}(1.2 - 0.2{x^5}),\;x \le 1;}\\ {1 - \dfrac{{{f_{{\rm{t}},{\rm{r}}}}}}{{{E_{\rm{c}}}{{\rm{\varepsilon }}_{t,{\rm{r}}}}\left[ {{a_{\rm{t}}}{{(x - 1)}^{1.7}} + x} \right]}},\;x > 1} \end{array}} \right. $ (5)
$ x = \frac{{\rm{\varepsilon}} }{{{{\rm{\varepsilon}} _{{\rm{t}},{\rm{r}}}}}} $ (6)

式中: ${\rm{\sigma}} $ ${\rm{\varepsilon}} $ 分别为混凝土的应力与应变, ${f_{{\rm{t}},{\rm{r}}}}$ 为混凝土单轴抗拉强度; ${\varepsilon _{{\rm{t}},{\rm{r}}}}$ ${f_{{\rm{t}},{\rm{r}}}}$ 对应的应变; ${d_{\rm{t}}}$ 为单轴受拉损伤变量; ${E_{\rm{c}}}$ 为混凝土弹性模量; ${a_{\rm{t}}}$ 为混凝土受拉应力应变曲线软化段参数。

混凝土单轴受压应力应变曲线表达式:

$ {\rm{\sigma}} = \left( {1 - {d_{\rm{c}}}} \right){E_{\rm{c}}}{{\varepsilon}} $ (7)
$ {\;\;\;\;\;\;\;\;\;d_{\rm{t}}} = \left\{\!\!\!\! {\begin{array}{*{20}{l}} {1 - \dfrac{{n{f_{{\rm{c}},{\rm{r}}}}}}{{{E_{\rm{c}}}{{\rm{\varepsilon }}_{{\rm{c}},{\rm{r}}}}\left( {n - 1 + {x^n}} \right)}},\;x \le 1;}\\ {1 - \dfrac{{{f_{{\rm{c}},{\rm{r}}}}}}{{{E_{\rm{c}}}{{\rm{\varepsilon }}_{{\rm{c}},{\rm{r}}}}\left[ {{a_{\rm{c}}}{{(x - 1)}^2} + x} \right]}},\;x > 1} \end{array}} \right. $ (8)
$ n = \frac{{{E_{\rm{c}}}{{\rm{\varepsilon}} _{{\rm{c}},{\rm{r}}}}}}{{{E_{\rm{c}}}{{{\varepsilon}} _{{\rm{c}},{\rm{r}}}} - {f_{{\rm{c}},{\rm{r}}}}}} $ (9)
$ x = \frac{{\rm{\varepsilon}} }{{{{\rm{\varepsilon}} _{{\rm{c}},{\rm{r}}}}}} $ (10)

式中, ${f_{{\rm{c}},{\rm{r}}}}$ 为混凝土单轴抗压强度, ${{\rm{\varepsilon}} _{{\rm{c}},{\rm{r}}}}$ ${f_{{\rm{c}},{\rm{r}}}}$ 对应的应变, ${d_{\rm{c}}}$ 为单轴受压损伤变量, ${a_{\rm{c}}}$ 为混凝土受压应力应变曲线软化段参数。

当结构处于多轴状态时,需用等效应力、应变 ${{\rm{\sigma}} ^ * }$ ${\varepsilon ^ * }$ 替代真实应力应变 ${\rm{\sigma}} $ ${\rm{\varepsilon}} $ ,通过应变张量第一不变量 $I_1' $ 判断单元拉压状态。

2 应力卸载残余应变量值

残余应变是混凝土达到软化段后客观存在、不可恢复的塑性变形。由于试验设计方案或是其他原因,试验结果所拟合的经验公式难以捕捉残余应变的临界值,因此通常将其忽略,但当荷载持续时间足够长或加卸载次数足够多时,不考虑残余应变临界值的计算公式最终将会产生偏差,这种偏差表现为超出临界值后混凝土卸载弹性模量 ${E_{\rm{c}}} \le {\rm{0}}$

表1列举了经典混凝土残余应变计算公式,本文选用Vecchio与Palermo公式[2]。设 ${k_{\rm{p}}} = \dfrac{{{{\rm{\varepsilon}} _{\rm{p}}}}}{{{{\rm{\varepsilon}}_{\rm{r}}}}}$ ${k_{{\rm{u}}{\rm{n}}}} = \dfrac{{{{\rm{\varepsilon}} _{{\rm{u}}{\rm{n}}}}}}{{{{\rm{\varepsilon }}_{\rm{r}}}}}$ ,其中, ${{\rm{\varepsilon}} _{\rm{p}}}$ 为塑性残余应变, ${{\rm{\varepsilon}} _{\rm{r}}}$ 为抗拉或抗压峰值强度对应应变, ${{\rm{\varepsilon}} _{{\rm{u}}{\rm{n}}}}$ 为卸载点处应变。

表1 经典混凝土残余应变计算公式 Tab. 1 Classical formula for calculating residual strain of concrete

${k_{\rm{p}}}$ ${k_{{\rm{un}}}}$ 关系如图1所示。当 ${k_{{\rm{u}}{\rm{n}}}} \approx 5.23$ 时, ${k_{\rm{p}}} = {k_{{\rm{u}}{\rm{n}}}}$ ,即卸载刚度 $E = 0$ ,因此取临界值 ${k_{{\rm{u}}{\rm{n}}}} = 4.5$ ;当 ${k_{{\rm{u}}{\rm{n}}}} \ge $ $ 4.5$ 时, ${{\rm{\varepsilon}} _{\rm{p}}}{\rm{ = 0}}{\rm{.85}}{{\rm{\varepsilon}} _{{\rm{u}}{\rm{n}}}}$

图1 残余应变临界值 Fig. 1 Critical value of residual strain

3 滞回规则特征点及路径

通过混凝土循环荷载试验观测发现,循环荷载下的应力应变响应取决于荷载历史。同时,滞回效应不局限于完全卸载和完全重新加载,存在卸载不完全的局部卸载循环与重新加载不完全的局部重新加载循环。Otter等[13]通过观测试验数据,建立了一套基于荷载历史的数学经验公式,以推求循环荷载下的加卸载特征点。公式对素混凝土及钢筋混凝土均有良好的适用性,且不依赖于应力应变骨架线的形状。

3.1 完全加卸载循环

滞回规则中完全加卸载循环如图2所示。滞回循环由卸载路径AC和重新加载路径CB组成,滞回期间的损伤积累程度在模型中体现为骨架线上应变由卸载点A发展至重新加载点B

图2 损伤模型完全加卸载示意图 Fig. 2 Diagram of complete unloading and reloading of damage model

根据卸载点A处卸载应变 ${{\rm{\varepsilon}} _{{\rm{un}}}}$ 计算重新加载点B处应变 ${{\rm{\varepsilon}} _{{\rm{re}}}}$

$ \frac{{{{{\varepsilon}} _{{\rm{r}}{\rm{e}}}}}}{{{{\rm{\varepsilon}} _{\rm{r}}}}} = \frac{{{{{\varepsilon}} _{{\rm{u}}{\rm{n}}}}}}{{{{\rm{\varepsilon}} _{\rm{r}}}}} + {k_{\rm{r}}} $ (11)

式中, ${k_{\rm{r}}}$ 为重新加载系数,一般取建议值0.1。

滞回规则中,卸载曲线的曲率反映了刚度变化过程,割线模量由大到小不断改变。重新加载曲线可简化为线性,期间割线模量保持恒定,损伤值不发生变化[4]。滞回规则完全卸载曲线采用Sima等[14]通过试验拟合的经验表达式,式中计入损伤变量,能够反映卸载过程中损伤积累与刚度变化过程:

$ {\;\;\;\;\;\;\;\; \rm{\sigma}} = {\xi _1}{{\rm{e}}^{{\xi _2}\left( {1 - \frac{{{\rm{\varepsilon}} - {{\rm{\varepsilon}} _{\rm{p}}}}}{{{{\rm{\varepsilon}} _{{\rm{u}}{\rm{n}}}} - {{\rm{\varepsilon}} _{\rm{p}}}}}} \right)}}{E}\left( {{\rm{\varepsilon}} - {{\rm{\varepsilon}} _{\rm{p}}}} \right) $ (12)
$ {\!\!\!\!\!\! \rm{\sigma}} = \frac{{{{\varepsilon}} - {{{\varepsilon}} _{\rm{p}}}}}{{{{{\varepsilon}} _{{\rm{r}}{\rm{e}}}} - {{\rm{\varepsilon}} _{\rm{p}}}}}{{{\sigma}} _{{\rm{r}}{\rm{e}}}} $ (13)

式中, ${\xi _1} = \dfrac{{r\left( {1 - {d_{{\rm{u}}{\rm{n}}}}} \right)}}{{r - 1}}$ $r = \dfrac{{{{\rm{\varepsilon}} _{{\rm{u}}{\rm{n}}}}}}{{{{\rm{\varepsilon}} _{\rm{p}}}}}$ ${\xi _2} = \ln \left[ {\dfrac{{R\left( {1 - {d_{{\rm{u}}{\rm{n}}}}} \right)\left( {r - 1} \right)}}{r}} \right]$ ${d_{{\rm{u}}{\rm{n}}}}$ 为骨架线卸载点处损伤值, ${{\rm{\sigma}} _{{\rm{r}}{\rm{e}}}}$ 为骨架线重新加载点应力值, $R = \dfrac{{{E_{\rm{p}}}}}{{{E}}}$ ${E_{\rm{p}}}$ 为完全卸载时的割线模量。

重新加载阶段割线模量将保持恒定,残余应变点处的损伤值与重新加载点处的损伤值相同,因此,卸载产生的损伤积累变化量为 ${d_{{\rm{r}}{\rm{e}}}} - {d_{{\rm{u}}{\rm{n}}}}$ ,卸载时损伤计算公式为:

$ d = {d_{{\rm{u}}{\rm{n}}}} + \frac{{{d_{{\rm{r}}{\rm{e}}}} - {d_{{\rm{u}}{\rm{n}}}}}}{{{{{\varepsilon}} _{\rm{p}}} - {{{\varepsilon}} _{{\rm{u}}{\rm{n}}}}}}\left( {{{\varepsilon}} - {{{\varepsilon}} _{{\rm{u}}{\rm{n}}}}} \right) $ (14)

当卸载至残余应变点时, $d = {d_{{\rm{re}}}}$

3.2 局部重新加载循环

滞回规则中,局部重新加载循环如图3所示。

图3 损伤模型局部重新加载示意图 Fig. 3 Diagram of local reloading of damage model

图3可知,滞回循环由卸载路径AD和重新加载路径DE组成。此时荷载未完全卸载至残余应变点C,相应重新加载点E与完全加卸载循环中的重新加载点B不同,其应变值将是点A与点B之间的插值:

$ {\;\;\;\;\;\;\;\;{\rm{\varepsilon}} _{{\rm{r}}{\rm{x}}}} = {{\rm{\varepsilon}} _{{\rm{u}}{\rm{n}}}} + \left( {{{\rm{\varepsilon}} _{{\rm{r}}{\rm{e}}}} - {{\rm{\varepsilon}} _{{\rm{u}}{\rm{n}}}}} \right){\left( {\frac{{{{\rm{\sigma}} _{{\rm{u}}{\rm{n}}}} - {{\rm{\sigma}} _{\rm{u}}}}}{{{{\rm{\sigma}} _{{\rm{u}}{\rm{n}}}}}}} \right)^{{n_{{\rm{p}}{\rm{u}}}}}} $ (15)

式中: ${{{\varepsilon}} _{{\rm{rx}}}}$ 为重新加载点E处应变值; ${{{\sigma}} _{\rm{u}}}$ 为局部卸载最低点D处应力值; ${n_{{\rm{pu}}}}$ 为插值参数,经试验拟合与敏感性分析,取建议值8。

卸载后的局部重新加载公式与完全卸载时相同:

$ {\!\!\!\!\!\!\!\!\!\!\!\! \rm{\sigma}} = \frac{{{\rm{\varepsilon}} - {{\rm{\varepsilon}} _{\rm{u}}}}}{{{{\rm{\varepsilon}} _{{\rm{r}}{\rm{x}}}} - {{\rm{\varepsilon}} _{\rm{u}}}}}{{\rm{\sigma}} _{{\rm{r}}{\rm{x}}}} $ (16)

式中, ${{\rm{\varepsilon}} _{\rm{u}}}$ 为局部卸载最低点D处应变值, ${{\rm{\sigma}} _{{\rm{rx}}}}$ 为重新加载点E处应力值。

3.3 局部卸载循环

滞回规则中,局部卸载循环如图4所示。

图4 损伤模型局部卸载示意图 Fig. 4 Diagram of local unloading of damage model

图4可知,滞回循环由卸载路径FG和重新加载路径GH组成。图4中局部加载最高点F对应的骨架线卸载点应变值将是点A与点B之间的插值:

$ {{\rm{\varepsilon}} _{{\rm{ux}}}} = {{\rm{\varepsilon}} _{{\rm{un}}}} + \left( {{{\rm{\varepsilon}} _{{\rm{re}}}} - {{\rm{\varepsilon}} _{{\rm{un}}}}} \right){\left( {\frac{{{{\rm{\sigma}} _{\rm{x}}} - {{\rm{\sigma}} _{\rm{u}}}}}{{{{\rm{\sigma}} _{{\rm{re}}}} - {{\rm{\sigma}} _{\rm{u}}}}}} \right)^{{n_{{\rm{pr}}}}}} $ (17)

式中: ${{{\varepsilon}} _{{\rm{ux}}}}$ F点对应骨架线卸载点的应变值; ${{{\sigma}} _{\rm{x}}}$ F点应力; ${n_{{\rm{pr}}}}$ 为插值参数,取建议值8。

局部卸载循环中卸载路径与完全加卸载循环类似:

$ {{\sigma}} = {{{\eta}} _1}{{\rm{e}}^{{{{\eta}} _2}\left( {1 - \frac{{{{\varepsilon}} - {{{\varepsilon}} _{\rm{p}}}}}{{{{{\varepsilon}} _{\rm{x}}} - {{{\varepsilon}} _{\rm{p}}}}}} \right)}}E\left( {{{\varepsilon}} - {{{\varepsilon}} _{\rm{p}}}} \right) $ (18)

式中, ${{\rm{\eta}} _1} = \dfrac{{{{\rm{\sigma}} _{\rm{x}}}}}{{E{{\rm{\varepsilon}} _{\rm{x}}}\left( {1 - {r_1}} \right)}}$ ${{\rm{\eta}} _2} = \ln {\dfrac{{{R_1}}}{{{{\rm{\eta}} _1}}}}$ ${r_1} = \dfrac{{{{\rm{\varepsilon}} _{\rm{x}}}}}{{{{\rm{\varepsilon}} _{\rm{p}}}}}$ ${R_1} = \dfrac{{{E_{\rm{p}}}}}{{{E}}}$ ${{\rm{\varepsilon}} _{\rm{p}}}$ ${{\rm{\varepsilon}} _{{\rm{ux}}}}$ 对应的残余应变值。

4 损伤模型在程序中的实现

通过Fortran语言编译相应材料子程序,并嵌入大型有限元仿真软件,本构模型数值实现流程如图5所示。

图5 损伤模型数值实现流程 Fig. 5 Flow chart of numerical realization of damage model

损伤模型构建于等效应变空间中,因此模型集中于第1象限,滞回规则的复杂性决定了数值模拟程序实现的关键在于判断当前时间步单元加卸载状态。将单元加卸载状态分为:骨架线加载状态,记 $S\!\!=1$ ;骨架线卸载状态,记 $S\!\!=2$ ;重新加载状态,记 $S\!\!=3$ ;局部卸载状态,记 $S\!\!=4$ ;拉压转换状态,记 $S\!\!=5$ 。通过下标 $n$ $n-1$ 注明变量值与时间步之间关系。

1)当 ${\rm{\varepsilon}} _n^ * \ge {\rm{\varepsilon}} _{n - 1}^ * = {\rm{\varepsilon}} _{\max }^ * $ 时,则 ${S\!_n} = 1$ ,记录 ${\rm{\varepsilon}} _{{\rm{re}}}^ * $

2)当 ${\rm{\varepsilon}} _{\rm{p}}^ * \le {\rm{\varepsilon}} _n^ * \le {\rm{\varepsilon}} _{n - 1}^ * \le {\rm{\varepsilon}} _{\max }^ * $ 时,且 ${S\!_{n - 1}} = 1$ ${S\!_{n - 1}} = 2$ ,则记录 ${S\!_n} = 2$ ${\rm{\varepsilon}} _{\rm{u}}^ * = {\rm{\varepsilon}} _n^ * $

3)当 ${\rm{\varepsilon}} _{\rm{p}}^ * \le {\rm{\varepsilon}} _{n - 1}^ * \le {\rm{\varepsilon}} _n^ * \le {\rm{\varepsilon}} _{{\rm{re}}}^ *$ 时,且 ${S\!_{n - 1}} \ne 1$ ,则记录 ${S\!_n} = 3$ ${\rm{\varepsilon}} _x^ * = {\rm{\varepsilon}} _n^ * $

4)当 ${\rm{\varepsilon }}_{\rm{p}}^ * \le {\rm{\varepsilon }}_n^ * \le {\rm{\varepsilon }}_{n - 1}^ * \le {\rm{\varepsilon }}_x^ * <{\rm{\varepsilon }}_{\max }^ * $ 时,且 ${S\!_{n - 1}} = 3$ ${S\!_{n - 1}} =$ $ 4$ ,则记录 ${S\!_n} = 4$ ${\rm{\varepsilon}} _{\rm{u}}^ * = {\rm{\varepsilon}} _n^ * $

5)当 ${\rm{\varepsilon}} _n^ * \le {\rm{\varepsilon}} _{n - 1}^ * $ ,且 ${\rm{\varepsilon}} _n^ * \le {\rm{\varepsilon}} _{\rm{p}}^ *$ ${\rm{\varepsilon}} _n^ * \ge {\rm{\varepsilon}} _{n - 1}^ *$ ,且 ${S\!_{n - 1}} = 5$ 时,记录 ${S\!_n} = 5$

5 算例验证 5.1 单轴循环拉伸荷载作用数值验证

循环拉伸荷载算例的建立依据Vellore等[15]开展的循环拉伸试验,建立1 m×1 m×1 m的3维8结点等参单元模型对本文混凝土损伤模型进行数值验证,模型如图6所示。由图6可知,对单元左侧及底侧的结点分别施加法向约束,对右侧结点依照试验数据进行应变分级加载。

图6 等参单元2维平面示意图 Fig. 6 2–dimensional plane diagram of isoparametric elements

仿真计算材料参数见表2

表2 算例材料参数 Tab. 2 Material parameters of the example

图7可知,本文模型能够较好地反映混凝土软化段滞回效应,特别是与滞回环中的行为特征点拟合较好,即每个滞回环的卸载最低点与重新加载起始点,通常这两点是试验观测的真实数据采集点。

图7 循环拉伸荷载下应力应变曲线对比 Fig. 7 Comparison of stress–strain curves under cyclic tensile load

图8对比了循环拉伸荷载下不同模型与试验数据的损伤历程。由于试验设计属于分级加卸载,因此在分级加卸载特征值处将出现数据重合,但损伤积累历程存在较大差异,未考虑滞回效应将忽略卸载阶段的损伤积累,实际工程案例中结构将随损伤演化出现应力重分布,忽略卸载时段的损伤积累将导致后续的仿真结果大相径庭。

图8 循环拉伸荷载下损伤历程对比 Fig. 8 Comparison of damage history under cyclic tensile load

5.2 Koyna重力坝震况验证

选取遭受实际震害的Koyna重力坝作为研究对象进行工程结构数值验证,坝体模型如图9所示。采用时程法在静力分析的基础上进行动力分析,通过广义Newmark法确定每一时刻坝体与地基的应力分布及变形情况。采用固定人工边界作为地基边界条件,即在无质量地基模型中地基的截断侧边界上施加法向的固定约束,底边界上施加双向固定约束,惯性力只施加在坝体上,地基上不加惯性力,只考虑其刚度。地震荷载选取Koyna地震波,其归一化的加速度时程曲线如图10所示,地震波时长12.8 s,水平向峰值加速度为0.474g,竖直向峰值加速度为0.312g。采用Westgaard附加质量法考虑地震荷载下库水作用于坝体的动水压力。

图9 Koyna重力坝示意图 Fig. 9 Diagram of Koyna gravity dam

图10 Koyna地震波加速度时程曲线 Fig. 10 Acceleration time histories of Koyna seismic wave

图1112分别为本文模型模拟和文献试验中Koyna坝体损伤分布。

图11 不同时刻本文模型仿真Koyna坝体损伤分布 Fig. 11 Damage distribution of Koyna dam at different times in the simulation results of the model in this paper

图12 不同时刻文献试验中Koyna坝体损伤分布[16] Fig. 12 Damage distribution of Koyna dam at different times in literature test[16]

图11图12试验结果[16]对比可知,图11中本文仿真的Koyna重力坝损伤演化过程和最终开裂破坏模式与试验表现基本一致;坝体2.7 s左右在坝坡折角处进入损伤,随着向上游面方向延伸,3.4 s左右上游面进入损伤并向坝坡折角处发展,最终形成贯穿破坏。损伤演化速度较图13中传统弹塑性模型仿真结果[17]更快,损伤演化过程也存在一定差异。计入混凝土滞回效应后,将考虑软化段卸载时的损伤积累,其造成的结构局部应力重分布使得整体动力响应结果更符合实际工程震害情况。

图13 不同时刻Koyna坝体损伤分布[17] Fig. 13 Damage distribution diagram of Koyna dam at different times[17]

6 结 语

在混凝土四参数损伤模型的基础上搭建滞回规则,建立循环荷载下考虑滞回效应的混凝土损伤模型。模型基于应变等效原理,将复杂的多轴受力变形问题转换至简单的等效应变空间中求解,通过以等效应变为自变量的单轴损伤演化方程求解损伤变量,随后得到多轴状态下的真实应力。模型不仅能够考虑循环往复荷载作用下混凝土材料拉压异性、刚度退化、强度软化及塑性不可恢复变形等复杂的非线性特征,同时能够完整描述混凝土滞回效应,包括其导致的材料刚度退化与损伤积累。通过单轴循环拉伸荷载算例和Koyna重力坝震况验证,模型能够较好地反映混凝土循环往复荷载作用下的受力状态及变形规律,计算过程不依赖于结构形式、加载规律及材料参数,能够为进一步研究混凝土结构抗震性能提供支撑。

参考文献
[1]
Hooshang D,Farhad A.A comparative study on the cyclic constitutive models of concrete[C].4th International Conference on Construction Materials–Performance,Innovations and Structural Implications.Nagoya,2009:163–170.
[2]
Vecchio F J,Palermo D. Compression field modeling of reinforced concrete subjected to reversed loading:Verification[J]. Aci Structural Journal, 2003, 100(5): 155-164.
[3]
Mansour M,Hsu T T C. Behavior of reinforced concrete elements under cyclic shear Ⅱ:Theoretical model[J]. ASCE Structural Journal, 2005, 131(1): 54-65. DOI:10.1061/(ASCE)0733-9445(2005)131:1(54)
[4]
Sakai J,Kawashima K. Unloading and reloading stress–strain model for confined concrete[J]. Journal of Structural Engineering, 2006, 132(1): 112-122. DOI:10.1061/(ASCE)0733-9445(2006)132:1(112)
[5]
Konstantinidis D,Kappos A J,Izzuddin B A.Analytical model for unconfined and confined high strength concrete under cyclic loading[C].13th World Conference on Earthquake Engineering.Vancouver,2004.
[6]
Aslani F,Jowkarmeimandi R. Stress-strain model for concrete under cyclic loading[J]. Magazine of Concrete Research, 2012, 64(8): 673-685. DOI:10.1680/macr.11.00120
[7]
过镇海.钢筋混凝土原理[M].3版.北京:清华大学出版社,2013.
[8]
Li Zheng,Zhu Bingyin. Development of uniaxial cyclic constitutive model for concrete based on Chinese concrete code[J]. Building Structure, 2013, 43(Supp1): 726-730. [李正,朱炳寅. 基于规范的混凝土单轴循环本构模型的开发[J]. 建筑结构, 2013, 43(增刊1): 726-730. DOI:10.19701/j.jzjg.2013.s1.166]
[9]
中国建筑科学研究院.混凝土结构设计规范:GB50010—2010.[S].北京:中国建筑工业出版社,2015.
[10]
Wei Wei,Li Tongchun,Yao Weiming. Four-parameter failure criterion for concrete in strain space[J]. Advances in Science and Technology of Water Resources, 2004, 24(5): 27-29. [韦未,李同春,姚纬明. 建立在应变空间上的混凝土四参数破坏准则[J]. 水利水电科技进展, 2004, 24(5): 27-29. DOI:10.3880/j.issn.1006-7647.2004.05.007]
[11]
Wei Wei,Li Tongchun. A new four-parameter equivalent strain for isotropic damage model[J]. Engineering Mechanics, 2005, 22(6): 91-96. [韦未,李同春. 一种新的用于各向同性损伤模型的四参数等效应变[J]. 工程力学, 2005, 22(6): 91-96. DOI:10.3969/j.issn.1000-4750.2005.06.016]
[12]
Zhou Qiujing,Zhang Guoxin,Li Tongchun. Analysis on working performance of concrete dams with a dynamic multi-axis equivalent strain damage[J]. Water Power, 2014, 40(12): 26-30. [周秋景,张国新,李同春. 基于多轴等效应变动力损伤模型的混凝土坝工作性态分析[J]. 水力发电, 2014, 40(12): 26-30. DOI:10.3969/j.issn.0559-9342.2014.12.009]
[13]
Otter D E,Naaman A E. Model for response of concrete to random compressive loads[J]. Journal of Structural Engineering, 1989, 115(11): 2794-2809. DOI:10.1061/(ASCE)0733-9445(1989)115:11(2794)
[14]
Sima F,Roca P,Molins C. Cyclic constitutive model for concrete[J]. Engineering Structures, 2008, 30(3): 695-706. DOI:10.1016/j.engstruct.2007.05.005
[15]
Vellore,Surendra P S. Softening response of plain concrete in direct tension[J]. Aci Materials Journal, 1985, 82(3): 310-323. DOI:10.14359/10338
[16]
Wang Chao,Zhang Sherong,Li Man,et al. Classification of earthquake damage to gravity dams based on a damage index model[J]. Earthquake Engineering and Engineering Dynamics, 2014, 1(6): 218-226. [王超,张社荣,黎曼,等. 基于损伤指数模型的重力坝地震破坏等级划分[J]. 地震工程与工程振动, 2014, 1(6): 218-226. DOI:10.13197/j.eeev.2014.06.218.wangc.027]
[17]
Guo Shengshan.Study on seismic damage process mechanism and quantitative evaluation criteria of concrete dam-foundation system based on parallel computation[D].Beijing:China Institute of Water Resources and Hydropower Research,2013.
郭胜山.基于并行计算的混凝土坝–地基体系地震损伤破坏过程机理和定量评价准则研究[D].北京:中国水利水电科学研究院,2013.