引言

混凝土材料由于具有抗压强度高及适应性强等优点普遍应用于工程结构中. 作为一种人造石材，其组构成分多样，因而具有较为复杂的应力应变行为特性.对于混凝土处于复杂应力状态下的情况，在一般的工程结构中都能找到对应的应力工作状态.比如，桥梁结构受到往复循环加载的作用，梁顶和梁底分别对应循环压和拉应力路径.David等^[1]开展了混凝土循环压路径的试验研究，发现在等幅循环压缩作用下，混凝土切线模量逐渐减小，同时强度衰减剧烈.由于混凝土在地震载荷下的延性效应，学者们对于峰值后应变软化段的力学行为特性更为关注.上述基本特性已被国内外大量的试验结果证 实^[2-9]. 在循环加载下，混凝土的力学性能均出现不同程度的变化，随着应力幅值的增大以及循环次数的增加，割线模量均降低，强度也减小.国内外众多学者针对其本构关系开展了较为持久与深入的研究，也取得了一系列的成果.路德春等^[10]利用广义非线性准则作为屈服面，以塑性剪应变作为硬化参量，构造出能反映应变软化特性的三维弹塑性静力模型，模型具有简单易用的特点.郑付刚等^[11]针对混凝土裂纹破坏模式，在有效应力空间中利用受拉与受剪损伤因子来反映微观损伤的影响，所建立的模型可反映多轴应力状态下的非线性力学行为特性.王怀亮等^[12]利用Perzyna黏塑性连续理论建立了一个适用于碾压混凝土层面动态剪切断裂行为的应力应变模型，使用Carol率相关性界面方程作为屈服判据来描述率相关性，采用塑性断裂力学理论来模拟剪切面上的断裂及摩擦活动行为，能较好地模拟碾压混凝土单调加载下的力学行为.

采用损伤力学与塑性理论相结合的建模方法可较为合理地考虑混凝土由于微裂缝以及不可逆变形等产生的非线性特性.但采用有效应力空间需先确定真实应力，这增加了一定复杂性.若采用名义应力，通过压损伤以及拉损伤分别定义相应的硬化参量，再通过构造较为适用于混凝土屈服特性的屈服面，则可在弹塑性框架体系内简单方便地模拟混凝土的多轴非线性特性以及由于裂纹所造成的强度降低以及刚度衰化等主要特性.由Bazant等^[13-14]提出的微平面模型基于材料点各方向的微平面上应变可由应变分量向各方向的平面投影得到，这就使得微平面上的应力应变关系能摆脱应力路径的影响，微平面上应力由体积应力、偏斜应力及切向应力组成，与面接触相关的力学行为比如摩擦、滑移等可直接反映.然而构建宏观与微平面微观应力应变的关系采用动态约束来表达，而动态约束中微平面模型中平面上应力应变分量更新相对独立，某些面上往往与宏观应力应变不一致，这就使得积分后得到宏观应力应变与试验产生较大的误差.Hsieh等^[15]在$I_{1}-J_{2}$坐标系中采用抛物线作为边界面，能反映出单轴压、拉及双轴压和三轴压等的强度特性.然而，在复杂加载情况下采用随动硬化法则表现循环加载及复杂加载特性.随动硬化法则无法从机理上揭示拉压两种不同加载方式所造成的混凝土损伤，因而也就无法准确反映出复杂的非线性特性.此外，还存在以下两个问题：(1) 在中低幅往复循环载荷作用下，例如对于单轴循环压缩情形(单轴抗压强度为$f_{\rm c}$)，当载荷为0.1$f_{\rm c}$左右时，即为小幅循环载荷，在小幅循环载荷作用下，滞回环趋于消失，混凝土的累积塑性变形达到一定值后几乎不再增大，形成类似弹性变形特性.Sparks^[16]的研究结果表明：在低应力水平下的循环加载，试验试件的割线模量在整个试验过程中基本没有降低.林燕清^[17]通过对混凝土单轴及双轴进行等幅循环加载疲劳试验，也得到了相同的规律. (2)由于边界面采用固定几何形态的抛物线，在静水压力或者广义剪应力显著变化时的静水压力效应无法得到合理反映.关于静水压力对于素混凝土的强度影响特点，Tasuji等^[18]在双向应力状态下的破坏特性进行了归纳，而Mills等^[19], Launay 等^[20], Gerstle等^[21]对于三向应力状态的破坏结果进行了总结，都表明在偏平面上破坏面横断面的轨迹在拉力或者小的压缩应力作用下为近似三角形形状曲线，而在逐渐增大的静水压力作用下，则三角形形状像圆形形状逐渐过渡，表明在较大的静水压力作用下，各向异性强度性质减弱. 鉴于此，根据上述3个主要问题，本文对Hsieh的弹塑性模型进行了改进.(1)采用边界面模型的思路，即假设外部存在一个边界面，内部应力点始终在边界面内部或者边界面上移动，而始终无法逾越出边界面，且边界面内部设置一个弹性域，弹性域边界与边界面几何相似，在弹性域内部，应力应变关系假设为线弹性关系，当应力点达到弹性域上则开始进入弹塑性工作状态，此时，应力应变关系按照边界面弹塑性关系进行插值计算，而当应力点最终达到边界面上时，则按照相关联流动法则和一致性条件进行弹塑性计算分析.如图 1所示，图中位于原点的内部封闭虚线为初始屈服面，屈服面内部为纯弹性域，而初始屈服面与后续屈服面之间为屈服面的弹塑性变形阶段. 假设初始屈服面在向外扩张时候，屈服面形态与进入后续阶段的屈服面几何相似.在初始屈服面与后续屈服面之间的塑性模量按照边界面的插值公式来修正.(2)基于广义非线性强度准则，并利用变换应力方法，将模型拓展为三维弹塑性本构模型. (3)混凝土按照其硬化成因及变形破坏机理可以相应的分为3种硬化机制：①压裂型硬化及破坏进程；②拉裂型硬化及破坏进程；③上述两种的混合型硬化过程及破坏模式.相应地根据上述加载模式以及混凝土的硬化受损机理，建议了受压及受拉作用下的损伤硬化参数.

1 混凝土模型 1.1 弹性域纯弹性阶段

在$\sigma _{1}$-$\sigma_{2}$坐标系下，在初始屈服面以内，所有的应力应变关系可假设为纯弹性变形，可用广义胡克定律来表示.

1.2 弹塑性阶段 1.2.1 偏平面上塑性模量 1.2.1.1 首次加载条件

根据Kupfer等^[22]的研究成果，在二维平面应力状态下，在首次加载时，混凝土处于初始屈服面与后续临界状态之间的弹塑性阶段，偏平面上的剪切模量可表述为

$ H_{\rm p} = H_{\rm p0} \left( \delta \right)^\omega = H_{\rm p0} \left( {1 - \dfrac{\tau _0 }{\tau _{\rm ou} }} \right)^\omega $

(1)

式中，$H_{\rm p0}$为初始塑性八面体剪切模量; $\delta $为$\sigma _1$-$\sigma _2 $坐标系下当前应力点与边界面上像点距离与像点到原点距离的比值; $\tau _0 $为八面体剪应力

$ \tau _0 = \sqrt { \dfrac{1}{3} s_{ij} s_{ij} } $

(2)

其中，$s_{ij}$为偏应力分量；$\tau _{\rm ou}$为$\tau _0$在边界面上的投影像点对应的八面体剪应力；$\omega $为像点剪应力的函数，为了能统一描述单轴拉到单轴压的单调比例加载路径的像点偏应变，可采用如下函数表示

$ \omega = n\left( {\dfrac{\tau _{\rm ou} }{f_{\rm c} }} \right)^m $

(3)

根据单调比例压缩路径，可得到像点偏应变的公式

$ \gamma _{\rm ou} = \dfrac{\tau _{\rm ou} }{H_0 \left[{1 - n\left( {\dfrac{\tau _{\rm ou} }{f_{\rm c} }} \right)^m} \right]} $

(4)

其中，$H_{0}$为初始八面体剪切模量,$ f_{\rm c}$为混凝土单轴抗压强度.

1.2.1.2 卸载条件

可采用Gerstle^[23]建议的公式，来表示卸载阶段的剪切模量

$ H_{\rm p} = H_{\rm p1} \left( {\delta ^\ast } \right)^\omega $

(5)

其中，$H_{\rm p1}$为卸载开始时的初始塑性八面体剪切模量,$\delta^{\ast}$为当前点到边界面的相反边界面上像点距离与像点和反边边界面像点距离之比值,确定H_p1可采用如下方程

$ H_{\rm p1} = \dfrac{\left( {\delta _{\max }^\ast } \right)^{1 - \omega } - \left( {\delta _{\min }^\ast } \right)^{1 - \omega }}{1 - \delta _{\min }^\ast } \cdot \\ \dfrac{\tau _{\rm ou} }{\left( {1 - \omega } \right)k_{\max } \gamma _{\rm ou}^{\rm p} \left( {1 - \eta } \right)}$

(6)

其中

$ k = \dfrac{\gamma _0^{\rm p} }{\gamma _{\rm ou}^{\rm p} } \quad $

(7)

$ k_{\max } = \max (k_i ) $

(8)

$ k_{\max }$表示第$i$次循环中达到的损伤参量最大值. $\eta=k_{\min}/k_{\max}$表示损伤参量最小与最大值之比.

1.2.1.3 再加载条件

根据Karsan等^[24]的建议，可采用如下公式来描述再加载阶段的剪切模量

$ H_{\rm p} = H_{\rm p2} \left( \delta \right)^\omega $

(9)

$H_{\rm p2}$为再加载开始时的初始塑性八面体剪切模量

$ H_{\rm p2} = \left( {1 - \delta _{\rm c}^{1 - \omega } } \right)\dfrac{\tau _{\rm ou} } {\left( {1 - \omega } \right)k_{\max } \gamma _{\rm ou}^p \left( {\xi - \eta } \right)} $

(10)

其中,$\xi=k_{\rm c}/k_{\max}$,表示公共点对应的损伤变量$k_{\rm c}$与损伤参量最大值之比.

$\eta =0.199{{k}_{\max }}+0.089,{{k}_{\max }}0.686,{{k}_{\max }}>3$

(11)

$\xi = 0.95 $

(12)

$\delta _c = 1 - x + x\delta _{\min } ,{x = 0.894} $

(13)

1.2.2 球应力作用下塑性体积模量 1.2.2.1 首次加载条件

根据Chen等^[25]的研究成果，初次加载下球应力作用的塑性体积模量可用如下公式表示

${{K}_{\text{p}}}=\frac{{{K}_{\text{p0}}}}{1+1.86{{\left( \frac{{{\sigma }_{0}}}{{{f}_{\text{c}}}} \right)}^{1.5}}}$

(14)

1.2.2.2 卸载及再加载条件

假设混凝土在卸载以及再加载条件下的体积塑性模量按照线性规律变化，则可用如下公式表示

$ K_{\rm p} = K_{\rm p0} = \dfrac{ - 0.013 68f_{\rm c}^2 + 1.78f_{\rm c} - 4.04}{\varepsilon _0 } $

(15)

1.2.3 偏应力作用下耦合体变系数

试验数据表明，纯剪切加载作用下也会引起体应变，顺序为初始时发生体积收缩，随后发生体积膨胀.根据Winnicki等^[26]的成果，初始加载时可采用如下简单函数来表述

$ \beta \left( \delta \right) = a\delta + b\sqrt \delta + c $

(16)

其中，$a$,$b$,$c$均为模型参数，可取值为$a =-2.403 3$,$b =3.163 3$,$c =-0.76$.

对于卸载或者再加载条件，则假定偏应力不引起体变，则

$ \beta \left( \delta \right) = 0 $

(17)

上述公式如用于描述首次加载阶段及卸载、再加载阶段中的应力应变关系，都是针对应变软化前的阶段，也就是用于描述应力峰值前的特性，如刚度衰减以及循环加载下的滞回环特性.

1.3 后续临界状态阶段 1.3.1 边界面方程

根据Hsieh等^[15]的研究成果，可修改为如下的边界面方程来表示

$ A\left( {\dfrac{\sqrt {J_2 } }{f_{\rm c}^\ast }} \right)^2 + B\dfrac{\sqrt {J_2 } }{f_{\rm c}^\ast } - C\dfrac{\sigma _3 }{f_{\rm c}^\ast } - D\dfrac{I_1 }{f_{\rm c}^\ast } - 1 = 0 $

(18)

其中，$I_1 $为应力张量第一不变量，$J_2 $为偏应力张量第二不变量，$\sigma _3 $为最小的主应力，$f_{\rm c}^\ast $为混凝土单轴压强度.

图 3为根据式(18) 得到的三维主应力空间中的边界面，由于混凝土自身的特性，边界面具备如下几个特点：

(1) 由于存在拉伸强度，因而边界面的圆锥顶点位于三维主应力空间的负半轴一侧，即满足：$\sigma_{1}<0$,$\sigma_{2}<0$,$\sigma_{3} <0$；

(2) 当应力罗德角$\theta =0^\circ$及$\theta =60^\circ$时，对应子午面上的边界面形状为幂函数曲线，因此能够反映随球应力水平的增加，广义剪应力增量减小现象；

(3) 随着球应力的增大，在不同偏平面上，对应的边界面由较尖锐的曲边三角形趋近于较圆钝化的曲边三角形，能一定程度的反映了偏平面上的静水压力效应，但由于边界面解析式为固定的方程式，因而无法反映不同应力路径下对静水压力效应的影响；

(4) 由于边界面在静水压力轴增大方向为开口形式，因而无法考虑沿着静水压力轴路径下的塑性变形特征.

采用Hsieh等^[15]的结果，用$A$,$B$,$D$来表示$C$，其中，$A =2.010 8$，$B =0.971 4$，$D =0.231 2$.

$ C = \dfrac{1}{3}A\dfrac{f_{\rm t}^\ast }{f_{\rm c}^\ast } - \dfrac{B}{\sqrt 3 } - D - \dfrac{f_{\rm c}^\ast }{f_{\rm t}^\ast } $

(19)

其中，$f_{\rm c}^\ast $和$f_{\rm t}^\ast $为损伤参量的函数. 可表示为

$f_{\text{c}}^{*}={{f}_{\text{c}}}{{f}_{1}}\left( {{k}_{\max }} \right)\text{ }$

(20)

$f_{t}^{*}={{f}_{t}}{{f}_{2}}\left( {{k}_{\max }} \right)\text{ }$

(21)

式中，$f_{1}$和$f_{2}$分别为后续临界状态压和拉强度降低函数.

1.3.2 压硬化损伤强度衰减函数

压损伤强度衰减函数如图 4所示.

对于压压状态导致的压强度降低函数可建议采用如下函数

${{f}_{1}}\left( {{k}_{\max }} \right)=1,{{k}_{\max }}\frac{1}{{{\beta }_{1}}\left( {{k}_{{{\max }^{2}}}}-1 \right)+{{k}_{\max }}},{{k}_{\max }}>1$

(22)

对于压压导致的拉强度降低函数可认为与式(22)完全相同. 上述公式基于Karsan 等^[24]的试验结果，对于一般混凝土的压碎，能够一定程度反映最大骨料粒径对于压缩损伤的影响. 一般随最大骨料粒径的增大，断裂能增大.可先确定出断裂能，然后经由断裂能来确定最大骨料粒径.

根据对于单轴压路径所得到的有限元,压分裂所释放的断裂能$G_{\rm f}$对主压应力与主压应变积分得

${{G}_{\text{f}}}=-{{h}_{\text{f}}}\int_{-\infty }^{{{\varepsilon }_{\text{c}}}}{{{\sigma }_{11}}}d\varepsilon _{11}^{\text{p}}=-{{h}_{\text{f}}}\int_{-\infty }^{{{\varepsilon }_{\text{c}}}}{{{f}_{\text{c}}}}{{f}_{1}}\left( {{k}_{\max }} \right)d\varepsilon _{11}^{\text{p}}$

(23)

将式(22)代入式(23)并推导，可得到最终的断裂能表达式

$ {G}_{\rm f} = \dfrac{ - \left( {\beta _1 + 2} \right)h_{\rm f} f_{\rm c} \gamma _{\rm ou}^{\rm p} }{\sqrt {2\left( {1 + 4\beta _1^2 } \right)} }\cdot \\ \left( {\ln \dfrac{\dfrac{\sqrt 2 \varepsilon _{\rm c} }{3\gamma _{\rm ou}^{\rm p} } + \dfrac{1}{2\beta _1 } - \sqrt {1 + \dfrac{1}{4\beta _1^2 }} }{\dfrac{\sqrt 2 \varepsilon _{\rm c} }{3\gamma _{\rm ou}^{\rm p} } + \dfrac{1}{2\beta _1 } + \sqrt {1 + \dfrac{1}{4\beta _1^2 }} } - 1} \right) $

(24)

其中，$h_{\rm f}$表示有限元的特征长度，$\varepsilon_{\rm c}$表示对应于强度$f_{\rm c}$的主压应变. 采用幂函数来表示断裂能是最大粒径的函数，则可表示为

$ {G}_{\rm f}= c_1 D_{\max }^{c_2 } $

(25)

其中，$c_{1}$和$c_{2}$经由试验来确定拟合.

1.3.3 拉硬化损伤强度衰减函数

由单轴拉拉状态或者拉压状态导致的拉强度降低函数可由以下公式描述. 图 5为形状参数$\beta _{2}$对拉硬化损伤衰减函数形状的影响.

${{f}_{2}}\left( {{k}_{\max }} \right)=1{{k}_{\max }}1{{\text{e}}^{-{{\beta }_{2}}\left( {{k}_{\max }}-1 \right)\gamma _{\text{ou}}^{\text{p}}}},{{k}_{\max }}>1$

(26)

由单轴拉拉路径及拉压路径导致的压强度降低函数可由下列公式表述

$ f_1 \left( {k_{\max } } \right) = - 0.4\left[{f_2 \left( {k_{\max } } \right) - 1} \right]^2 + 1 $

(27)

在后续临界状态，上述公式用于描述混凝土在峰值后应变软化阶段的特点，如强度减小，变形增大，所形成的滞回环的割线模量逐渐减小.上述峰值后的现象通过采用边界面的相关联流动法则来实现，此时边界面即退化为屈服面，并由压强度降低函数或拉硬化损伤强度衰减函数来分别控制边界面的流动过程.

2 基于广义非线性强度准则的变换应力方法 2.1 变换应力公式

由姚仰平等^[27-31]提出的广义非线性强度准则(general ronlinear strength theory,GNST)，可广泛地适用于描述砂土、黏土、岩石以及混凝土等材料的强度特性，基于该准则的变换应力方法则可简单合理的将岩土类材料由二维实现为三维应力应变关系模型. 其基本变换思路为将三维主应力空间的破坏面转换为变换应力空间中的圆锥面.将广义非线性强度准则变换为扩展Mises准则，变换分为两步：(1)将三轴压缩子午面上幂函数形式表示的破坏函数转换为过渡空间中的直线形式；(2)将偏平面上的曲面三角形状转换为变换应力空间中的规则圆形曲线. 见图 6和图 7.

变换应力公式为

$ \bar {p} = \dfrac{\bar {\sigma }_{ii} }{3} = p_{\rm r} \left( {\dfrac{p + \sigma _0 }{p_{\rm r} }} \right)^n ,i=1,2,3 $

(28)

式中，$p$，$\sigma _{0} $表示有效球应力和三向拉伸强度. $p_{\rm r}$为参考球应力，反映在此应力值下，子午面上破坏曲线的割线斜率可保证式(28)括号中的比值量纲为1，使等式左右量纲相同. 对于散粒体材料，$p_{\rm r}$通常取一个标准大气压值.

在偏平面上，剪应力可表示为

$ \bar {q}_\alpha ^ * = \alpha \sqrt {\bar {I}_{1}^{2} - {3}\bar {I}_{2} } + \dfrac{ 2(1 - \alpha ) \bar {I}_{1} }{{3}\sqrt {(\bar {I}_{1} \bar {I}_{2} - \bar {I}_{3} ) / (\bar {I}_{1} \bar {I}_{2} - {9}\bar {I}_{3} )} - {1}} $

(29)

其中

${{\bar{I}}_{1}}={{\bar{\sigma }}_{1}}+{{\bar{\sigma }}_{2}}+{{\bar{\sigma }}_{3}}{{\bar{I}}_{2}}={{\bar{\sigma }}_{1}}{{\bar{\sigma }}_{2}}+{{\bar{\sigma }}_{2}}{{\bar{\sigma }}_{3}}+{{\bar{\sigma }}_{3}}{{\bar{\sigma }}_{1}}{{\bar{I}}_{3}}={{\bar{\sigma }}_{1}}{{\bar{\sigma }}_{2}}{{\bar{\sigma }}_{3}}$

(30)

$\bar {\sigma }_i = \sigma _i + \left[{p_{\rm r} \left( {\dfrac{p + \sigma _{0} }{p_{\rm r} }} \right)^n - p} \right] , n \in [0, 1] $

(31)

$\bar {\sigma }_i $ 对应的应力空间为过渡应力空间.

一般化变换应力空间中的应力，则可得到变换应力空间中的变换应力张量$\tilde {\sigma }_{ij} $

$ \tilde {\sigma }_{ij} = \bar {p}\delta _{ij} + \dfrac{\bar {q}_\alpha ^\ast }{\bar {q}}\left( {\bar {\sigma }_{ij} - \bar {p}\delta _{ij} } \right) $

(32)

2.2 静水压力效应

试验结果表明，混凝土材料具备静水压力相关特性，即剪切强度会受到约束压力的影响，同时广义剪应力的增大也会使材料的各向异性程度增强.

图 8为在偏平面上的GNST强度线，随着静水压力的增大，强度线逐渐由较为尖锐的曲边三角形向曲边圆形过渡.图 9为相同的静水压力下随着广义剪应力增大的GNST曲线，可见，当剪应力增大时，则强度线形态由近似圆形向较尖锐的曲边三角形过渡.上述两个特性说明了GNST能够反映静水压力与广义剪应力的交叉耦合影响. 图 10所示为原边界面随静水压力增大的曲线族，由图可见，几乎呈现几何相似形变化，表明原边界面受静水压力影响很小.另外，由图 3和图 10可见，由于原边界面几何形态为由三个空间曲面拼接而成，因此，在拼接处存在着突变线，在图 10中则投影为与三条主应力轴重合的线. 按照相关联流动法则，当应力点位于此线上时则会造成屈服面法线方向不唯一，因此会产生奇异性.而GNST强度线则为连续光滑过渡曲线，不存在上述问题.

3 三维弹塑性混凝土本构模型 3.1 纯弹性域阶段本构方程

在纯弹性域范围内，则应力应变关系服从广义胡克定律

$d{{\sigma }_{ij}}=D_{ijkl}^{\text{e}}d{{\varepsilon }_{kl}}$

(33)

其中，$D_{ijkl}^{\rm e} $弹性刚度张量.

$ D_{ijkl}^{\rm e} = \dfrac{E}{2\left( {1 + \nu } \right)}\left( {\dfrac{2\nu }{1 - 2\nu }\delta _{ij} \delta _{kl} + \delta _{ij} \delta _{ik} \delta _{jl} + + \delta _{ik} \delta _{jl} } \right) $

(34)

3.2 弹塑性域阶段本构方程

在未进入后续临界状态阶段时，则在弹塑性状态阶段，假设变形都随变换应力在变换应力空间仍然产生. 式(1) $\sim$式(18)中的应力量分别被变换应力空间中的对应应力量所替换，在此不赘述.

弹塑性本构方程为

$d{{\sigma }_{ij}}=D_{ijkl}^{\text{ep}}d{{\varepsilon }_{kl}}$

(35)

其中，$D_{ijkl}^{\rm ep} $为弹塑性本构张量，可表示为

$\begin{align} & D_{ijkl}^{\text{ep}}=D_{ijkl}^{\text{e}}-\frac{{{{\tilde{K}}}_{\text{e}}}}{1+{{{\tilde{K}}}_{\text{p}}}/{{{\tilde{K}}}_{\text{e}}}}{{J}_{ij}}{{J}_{kl}}- \\ & \left( D_{ijkl}^{\text{e}}{{B}_{klst}}{{{\tilde{\sigma }}}_{st}}+\frac{\beta {{{\tilde{\tau }}}_{0}}{{{\tilde{K}}}_{\text{p}}}}{1+{{{\tilde{K}}}_{\text{p}}}/{{{\tilde{K}}}_{\text{e}}}}{{J}_{ij}} \right)\cdot \left( \frac{{{{\tilde{\sigma }}}_{st}}{{B}_{stmn}}D_{mnkl}^{\text{e}}}{3{{{\tilde{H}}}_{\text{p}}}\tilde{\tau }_{0}^{2}+{{{\tilde{\sigma }}}_{st}}{{B}_{stmn}}D_{mnkl}^{\text{e}}{{B}_{klpq}}{{{\tilde{\sigma }}}_{pq}}} \right) \\ \end{align}$

(36)

$J_{ij} = \delta _{ij} $

(37)

$B_{ijkl} = \dfrac{1}{3}\left( {6\delta _{ik} \delta _{jl} - 3\delta _{ij} \delta _{ik} \delta _{jl} - \delta _{ij} \delta _{kl} } \right) $

(38)

3.3 后续临界状态本构方程

边界面方程

$ A\left( {\dfrac{\sqrt {\tilde {J}_2 } }{\tilde {f}_{\rm c}^\ast }} \right)^2 + B\dfrac{\sqrt {\tilde {J}_2 } }{\tilde {f}_{\rm c}^\ast } - C\dfrac{\tilde {\sigma }_3 }{\tilde {f}_{\rm c}^\ast } - D\dfrac{\tilde {I}_1 }{\tilde {f}_{\rm c}^\ast } - 1 = 0 $

(39)

根据图 11所示，当对原边界面方程三维化后，所得到的边界面表达式为式(39)，其在图 11中的偏平面上的边界面为大曲边三角形，而对于变换应力空间中的真三轴等比加载任意路径，如沿着原点到$CAB$点路径，则其对应的广义偏应力$q =2 \sqrt 3 J_{2} /3$，对应着变换空间边界面与小圆的交点$A$点，由于变换应力中，将任意路径由三轴压缩路径来表达，因而，小圆上任意一点的应力状态均与$A$点相同，而根据变换应力的变换思路，变换应力空间中的小圆可与真实应力空间中的小曲边三角形存在一一映射关系，因而，变换应力空间的$A$点与真实应力空间中的$C$点相对应. 当调整变换空间中屈服面形态的参数$\alpha $变化时，则存在如下3种情况：(1)当变换应力空间中的参数$\alpha = 1$时，则变换应力空间中的边界面与真实应力空间中的边界面相同，均在偏平面上表示为一个圆，此时，变换空间中边界面与真实应力空间中边界面相重合，即$C$点退化为$A$点. (2)而当$\alpha =0$时，则变换应力空间中的边界面与真实应力空间中的边界面不同，此时，真实空间边界面退化为空间滑动面(spacial mobilized plane,SMP)强度准则屈服面. (3)当$0 < \alpha <1$时，则此时变换空间中的对应每个应力点所对应的圆，其在普通应力空间中的映射关系则对应于介于圆与SMP曲边三角形之间的屈服面形态.

由于对应的是增量型本构关系，因而对应每一个增量步，都会对应一个与之不同的曲面三角形，即图 11中的真实空间边界面，而由于此真实空间边界面是Mises圆与SMP准则之间的插值表达式，而SMP准则可以考虑由于应力路径的不同所带来的静水压力效应的影响，因而，这种三维化方法能够考虑由于应力路径的差异所导致的静水压力效应的影响.

式(39)经过整理可得

$ A\tilde {J}_2 + \left( {B\sqrt {\tilde {J}_2 } - C\tilde {\sigma }_3 - D\tilde {I}_1 } \right)\tilde {f}_{\rm c}^\ast - \tilde {f}_{\rm c}^{\ast 2} = 0 $

(40)

令$E = \left( {B\sqrt {\tilde {J}_2 } - C\tilde {\sigma }_3 - D\tilde {I}_1 } \right)$，代入式(40)，求解得

$ \tilde {f}_{\rm c}^\ast = \dfrac{1}{2}\left( {E + \sqrt {E^2 + 4A\tilde {J}_2 } } \right) $

(41)

由一致性条件，对式(41)进行全微分，可得到全微分方程，并联立下式

$d{{\sigma }_{ij}}=D_{ijkl}^{\text{e}}d\varepsilon _{kl}^{\text{e}}=D_{ijkl}^{\text{e}}\left( d{{\varepsilon }_{kl}}-d\varepsilon _{kl}^{\text{p}} \right)$

(42)

(1) 对于拉后续临界状态

$\begin{align} & {{Q}_{mn}}=(\frac{1.5B{{{\tilde{\tau }}}_{0}}}{\sqrt{{{{\tilde{J}}}_{2}}}}\frac{\partial {{{\tilde{\tau }}}_{0}}}{\partial {{{\tilde{\sigma }}}_{uv}}}\frac{\partial {{{\tilde{\sigma }}}_{uv}}}{\partial {{\sigma }_{mn}}}-C\frac{\partial {{{\tilde{\sigma }}}_{\min }}}{\partial {{{\tilde{\sigma }}}_{uv}}}\frac{\partial {{{\tilde{\sigma }}}_{uv}}}{\partial {{\sigma }_{mn}}}- \\ & D\frac{\partial {{{\tilde{I}}}_{1}}}{\partial {{{\tilde{\sigma }}}_{uv}}}\frac{\partial {{{\tilde{\sigma }}}_{uv}}}{\partial {{\sigma }_{mn}}}) \\ \end{align}$

(43)

$X_{mn} = \dfrac{6A\tilde {\tau }_0 \dfrac{\partial \tilde {\tau }_0 }{\partial \tilde {\sigma }_{uv} }\dfrac{\partial \tilde {\sigma }_{uv} }{\partial \sigma _{mn} }}{\sqrt {E^2 + 4A\tilde {J}_2 } } + \dfrac{1.6\tilde {f}_{\rm c} \left( {f_2 - 1} \right)f_2 \beta _2 \dfrac{\partial \tilde {\tau }_0 }{\partial \tilde {\sigma }_{uv} }\dfrac{\partial \tilde {\sigma }_{uv} }{\partial \sigma _{mn} }}{H_{\rm p0} \left( {1 - \omega } \right)\left( {1 - \tilde {\delta }} \right)}$

(44)

$M_{kl} = \left[{\left( {1 + \dfrac{E}{\sqrt {E^2 + 4A\tilde {J}_2 } }} \right)Q_{mn} + X_{mn} } \right]D_{mnkl}^{\rm e} $

(45)

$ Q_3 = \sqrt {\left( {\dfrac{\partial \tilde {f}}{\partial \tilde {\sigma }_{ij} } - \dfrac{1}{3}\dfrac{\partial \tilde {f}}{\partial \tilde {\sigma }_{ij} }\delta _{ij} } \right)\left( {\dfrac{\partial \tilde {f}}{\partial \tilde {\sigma }_{ij} } - \dfrac{1}{3}\dfrac{\partial \tilde {f}}{\partial \tilde {\sigma }_{ij} }\delta _{ij} } \right)} $

(46)

$d{{\sigma }_{ij}}=\left( D_{ijkl}^{\text{e}}-\frac{D_{ijuv}^{\text{e}}\frac{\partial \tilde{f}}{\partial {{{\tilde{\sigma }}}_{uv}}}{{M}_{kl}}}{{{M}_{kl}}\frac{\partial \tilde{f}}{\partial {{{\tilde{\sigma }}}_{kl}}}+\frac{1.6}{\sqrt{3}}{{{\tilde{f}}}_{\text{c}}}\left( {{f}_{2}}-1 \right){{f}_{2}}{{\beta }_{2}}{{Q}_{3}}} \right)d{{\varepsilon }_{kl}}$

(47)

(2) 对于压后续临界状态

$Y_{mn} = \dfrac{6A\tilde {\tau }_0 \dfrac{\partial \tilde {\tau }_0 }{\partial \tilde {\sigma }_{uv} }\dfrac{\partial \tilde {\sigma }_{uv} }{\partial \sigma _{mn} }}{\sqrt {E^2 + 4A\tilde {J}_2 } } $

(48)

$Z_{mn} = \dfrac{2\left( {2\beta _1 k_{\max } + 1} \right)k_{\max } \tilde {f}_{\rm c} \dfrac{\partial \tilde {\tau }_0 }{\partial \tilde {\sigma }_{uv} }\dfrac{\partial \tilde {\sigma }_{uv} }{\partial \sigma _{mn} }}{\left[{\beta _1 \left( {k_{\max }^2 - 1} \right) + k_{\max } } \right]^2\gamma _{\rm ou}^{\rm p} H_{\rm p0} \left( {1 - \omega } \right)\left( {1 - \tilde {\delta }} \right)} $

(49)

$N_{kl} = \left[{\left( {1 + \dfrac{E}{\sqrt {E^2 + 4A\tilde {J}_2 } }} \right)Q_{mn} + Y_{mn} - Z_{mn} } \right]D_{mnkl}^{\rm e} $

(50)

$Q_4 = \dfrac{2\left( {2\beta _1 k_{\max } + 1} \right)\tilde {f}_{\rm c} }{\sqrt 3 \left[ {\beta _1 \left( {k_{\max }^2 - 1} \right) + k_{\max } } \right]^2\gamma _{\rm ou}^{\rm p} } $

(51)

$d{{\sigma }_{ij}}=\left( D_{ijkl}^{\text{e}}-\frac{D_{ijuv}^{\text{e}}\frac{\partial \tilde{f}}{\partial {{{\tilde{\sigma }}}_{uv}}}{{N}_{kl}}}{{{N}_{kl}}\frac{\partial \tilde{f}}{\partial {{{\tilde{\sigma }}}_{kl}}}-{{Q}_{3}}{{Q}_{4}}} \right)d{{\varepsilon }_{kl}}$

(52)

4 模型预测

采用Kupfer的试验结果与模型预测进行对比，如图 12和图 13所示，双轴压缩工况下的双轴应力比值分别为1$:$0,1$:$1和1$:$0.52.

对于三轴压缩情况，采用Imran等^[32]对混凝土在侧向约束围压为0$\sim $17.2 MPa条件下的三轴应力应变关系结果进行计算模拟，如图 14和图 15所示.

采用Karsan 等^[24]的试验数据，对混凝土在循环载荷作用下的应力应变响应进行研究，分别如图 16和图 17所示.上述4种混凝土的基本参数见表 1 $\sim$表 3所示.

图 12和图 13分别为采用所提模型及各向异性增量本构模型对双轴压缩进行的预测结果对比.各向异性增量模型为Elwi和Murray所建议的一个三维增量模型，以下简称为EM模型.由对比可发现，所提模型对于双轴不同应力比下的压缩能较好的模拟，对于应变软化，强度衰化以及不同应力比下的变形特性有比较准确的模拟. 而EM模型对于应力比为1/0.52情况下的第二主应变在软化阶段与试验数据趋势相违背，在软化段，第二主应变仍然处于压缩阶段，但EM模型却计算出了剪胀特点.而所提模型预测结果与试验规律相一致.

图 14和图 15分别为两模型在三轴压缩条件下的对比情况.图 14中所提模型预测规律与试验结果基本一致，能反映出随不同围压条件下的第一主应变与第二主应变的变化规律，应变软化特点.随围压的约束作用，压缩的强度也相应的提高，预测结果显示其具备反映静水压力相关特性.而图 15中EM模型的对比结果表明，对于围压的约束作用，压缩强度的提高很小.

由图 16可见，随着循环次数的增加，呈现出如下规律：

(1) 初始剪切模量逐渐降低，由图 16试验数据可知，初次加载时剪切模量呈现出逐渐减小的规律，但初次加载曲线呈现出外凸性，而卸载后再加载曲线则出现微弱的外凹特点，随着循环次数的增加，再加载阶段曲线呈现出明显的外凸性，即剪切模量体现出先减小后增大到软化段再减小的特征. 从预测曲线来看，能反映出随着循环次数的增加，剪切模量逐渐减小的规律.(2)从卸载曲线分析，随着循环次数的增大，卸载时的剪切模量也逐步减小，而预测曲线也能反映出这种特征.(3)随着循环次数的增加，压应变的增大导致了压损伤的出现，随之导致了应变软化现象. 从预测曲线来看，也能明显的反映应变软化特征. (4)从图 16可见，其强度软化包络线能完整的表现出与单调加载时应变软化曲线相一致的现象.

对比单轴循环拉的试验数据以及预测曲线，从中可发现除了与上述循环压相同的规律外，由于不同的硬化以及损伤机理，因此循环拉作用下的应力应变曲线与循环压存在较大差异. 还存在以下特点：

(1) 单轴拉的强度较单轴压的强度值小很多，通常为单轴压路径的0.1倍左右.而随着循环次数的增大，单轴拉的强度值亦逐步降低，在8次单轴拉循环后，最终的残余强度值约为初始时的0.2倍左右.预测曲线较好的反映了这种现象.

(2) 由试验可知，在再加载阶段曲线中，始终遵循着剪切模量逐渐减小的特点，并未出现单轴循环压路径下的先减小后增大到软化段再减小的现象. 即其再加载曲线始终是外凸的，并未出现外凹特点.

(3) 由于单轴循环拉的损伤特点，其对应变的变化更为敏感，在轴应变为0.35%时，强度值即退化为初始强度值的0.2倍，而对于单轴循环压路径下，在轴应变为0.35%时，强度值退化为初始强度值的0.8倍左右. 因此循环拉造成的应变软化现象更为显著.

5 结论

通过对Hsieh的模型进行修正，使其克服了如下缺陷：

(1) 在边界面框架下，引入初级屈服面，在初级屈服面以内是弹性域.当处于循环加载的应力点位于弹性域内部时，则可按照纯弹性解来计算，能够反映较低应力水平下的重复循环加载特性.

(2) 基于广义非线性强度准则，利用变换应力方法对Hsieh模型实现三维化，使三维化后的Hsieh模型能够考虑不同应力路径对于当前空间边界面的形态影响，即一方面可通过子午面上的幂函数曲线来反映子午面上的静水压力效应，另一方面可通过偏平面上的广义非线性强度准则来反映偏平面上的广义剪应力与静水压力交叉耦合的影响. 采用变换应力方法，避免了原边界面存在应力奇异点的问题.

(3) 建议了压拉强度损伤函数，可用于描述对应于压拉两种不同加载模式下对混凝土造成的应变软化及强度衰减等行为.

通过对比双轴、三轴以及单轴循环压和拉等多种应力路径下的试验和预测结果，在静力载荷下，修正Hsieh模型可反映出应变软化性质，双轴作用下的静水约束特性，三轴作用下抗压强度的围压约束相关性，表现出强度的硬化特点以及应变硬化规律.在循环载荷下则可简单、合理的反映出混凝土的循环刚度衰减规律、强度的应变软化包络线，以及非线性特点.可比较方便的应用于一般应力路径下混凝土的应力应变关系模拟.