钢管混凝土（concrete-filled steel tube, CFST）结构因承载能力较高，塑性、抗震性能较好等优点，在大跨桥梁、高层及超高层建筑、地铁站等结构的受压构件中得到广泛运用，同时因其充分发挥了钢与混凝土材料的力学性能，表现出良好的经济效果^[1-3].

大量研究结果均表明，混凝土徐变会引起结构变形和内力重分布，对CFST的长期结构行为有重要影响^[4-6]. 试验研究和理论分析均表明，当核心混凝土应力超过一定水平时（0.6f_cm，f_cm为圆柱体抗压强度平均值），徐变与应力之间表现出明显的非线性特征^[7-9]. 付学宝等^[9]分析计算了厦门港务大厦、昆明大成金融商务中心等超高层建筑的CFST柱核心混凝土应力水平，发现部分已达0.69f_cm. 高应力下的非线性徐变效应不仅会使CFST结构出现过大的变形，还会进一步影响结构的动力特性与响应^[10-11]. 为准确分析高应力水平下的CFST长期结构行为，需要建立非线性徐变分析理论.

对于非线性徐变理论，现有的研究方法主要有3类：一类是直接通过高应力下的单轴受力徐变试验，建立模拟徐变特性的流变模型^[12-13]；另一类是以低应力下的单轴徐变特性为基础，建立徐变应变和应力之间的非线性函数关系^[14-17]；还有一类是采用损伤力学的方法来建立高应力下混凝土的非线性徐变本构方程^[18-19]. 前2类方法均是从单轴常应力作用下的徐变试验出发，很难将其运用到复杂空间应力状态下的徐变效应分析；第三类方法把非线性徐变产生的原因归结为高应力下的材料损伤演化，在实际工程中得到广泛运用，但已有研究成果也存在一定不足. 如李兆霞等^[18]运用连续介质损伤力学建立了弹性徐变损伤分析模型，并未涉及高应力下材料塑性行为的影响. 黄海东等^[19]运用经典的单轴应力应变关系（Hongnestad模型）探讨了混凝土非线性徐变效应的分析方法，而单轴应力应变曲线并不能完整描述混凝土材料的非线性行为.

除此之外，轴压圆CFST短柱核心混凝土应力状态具有明显的三向受压特征，因此，其徐变特性还具有明显的三维特性^[20-23].

为准确分析高应力下轴压圆CFST短柱徐变效应，本文基于混凝土弹塑性理论和连续损伤介质力学，提出一种新的高应力水平下的轴压圆CFST柱三维非线性徐变效应分析理论框架. 该框架首先包括均匀约束条件下的混凝土塑性损伤本构模型，用以模拟高应力作用下混凝土塑性和损伤演化；其次是CFST三维非线性徐变分析模型和相应的数值分析方法. 最后，本研究基于徐变预测B4模型^[24]，结合通用有限元软件ABAQUS，二次开发计算程序，通过不同应力水平下的CFST短柱徐变试验，验证所提方法的可靠性.

Lubliner等^[25-26]提出的混凝土塑性损伤本构模型因为具有概念清晰、简便易用等优点，在实际工程中得到广泛应用，该本构模型主要包括Lubliner屈服准则、Drucker-Prager塑性流动准则、损伤演化与硬化准则等，主要适用于低围压应力下的混凝土结构弹塑性行为分析. 对轴压圆CFST柱而言，环向围压已经影响到混凝土的力学行为，对低围压下的塑性损伤本构模型进行合理修正，以适用于均匀约束条件下混凝土的弹塑性行为分析，不失为一种便捷有效的方法.

Lubliner屈服准则的表达式如式（1）所示，应力空间屈服面在偏平面上的迹线如图1所示.

(1) $\begin{split} F\left( {{\bar{ {\sigma }}},{{\tilde \varepsilon }^{{\rm{p}}}}} \right) =& \displaystyle \frac{1}{{1 - \alpha }}\left( {\sqrt {3{J_2}} + \alpha {I_1} + \beta \left\langle \;\;{{{\hat {\bar \sigma} }_{\max }}} \right\rangle - } \right. \\ & \left. {\gamma \left\langle { - \;\;{{\hat {\bar \sigma }}_{\max }}} \right\rangle } \right) - c \leqslant {\rm{0}}. \end{split}$

式中： $\tilde\varepsilon^{\rm p} $为等效塑性应变，〈〉为Macaulay括号，〈x〉=(|x|+x)/2；I₁、J₂分别为第一应力张量不变量和第二偏应力张量不变量； ${{{\bar{{\sigma }}}}_{\max }}$ 为有效应力 ${{\bar \sigma }}$ 特征值的最大值；c为混凝土黏聚力；参数α、β由混凝土单轴拉压强度和双轴等压强度确定，具体过程可参见文献[26]；在文献[25]、[26]中，参数γ通过拉子午线（tensile meridian, TM）与压子午线（compression meridian, CM）来确定，这种方法对于轴压CFST柱而言，存在着无法避免的缺陷^[27]，详细阐述如下.

拉、压子午线在Haigh-Westergaard应力空间上的表达式如下：

(2) ${\rm{TM:}}\; \left( {\gamma + 3\alpha } \right)\xi + \frac{{2\gamma + 3}}{{\sqrt 2 }}{\rho _{{\rm_{tm}}}} = \sqrt 3 \left( {1 - \alpha } \right)c,$

(3) ${\rm{CM:}}\;\left( {\gamma + 3\alpha } \right)\xi + \frac{{\gamma + 3}}{{\sqrt 2 }}{\rho _{{\rm_{cm}}}} = \sqrt 3 \left( {1 - \alpha } \right)c.$

式中： ξ为静水压力，ξ=−I₁/3; ρ为Mises等效应力， $\rho = \sqrt {3{{J}}_2} $. 当ξ为定值时，

(4) ${K_{\rm{c}}} = \frac{{{\rho _{\rm_{tm}}}}}{{{\rho _{\rm_{cm}}}}}=\frac{{\gamma {\rm{ + 3}}}}{{{\rm{2}}\gamma {\rm{ + 3}}}}.$

K_c值决定了偏平面迹线的形状，如图1所示，也说明参数γ决定了屈服面的开口形状.

拉、压子午线对静水压力负半轴的斜率可表示为

(5) $ \tan {\varphi _{{\rm{tm}}}} = \frac{{\sqrt 2 \left( {\gamma + 3\alpha } \right)}}{{2\gamma + 3}}, \; \tan {\varphi _{{\rm{cm}}}} = \frac{{\sqrt 2 \left( {\gamma + 3\alpha } \right)}}{{\gamma + 3}}. $

式中：φ为混凝土内摩擦角. 斜率的数值大小表征了屈服面的开口角度，由参数α与γ共同决定，如图2所示.

由上述分析可以得出，参数γ同时决定了应力空间的屈服面的开口形状和开口角度，而参数γ仅与参数K_c有关. 因此，在选取合适的参数K_c来确定屈服面的形状后，就可能使屈服面无法具有合适的开口角度，而开口形状往往对于混凝土的弹塑性计算而言更加重要^[27-28]，这就是Lubliner屈服准则的主要缺陷所在.

对于轴压圆CFST短柱而言，参数γ应由tanφ_cm来确定^[27]. 因核心混凝土处于三向受压，混凝土抗压强度显著提高^[28]：

(6) ${f_{{\rm{cc}}}} = {f_{\rm{c}}} + k\;{\sigma _{\rm{L}}}.$

式中：f_cc为均匀约束混凝土抗压强度，σ_L为侧压力，k为强度提高系数. 当核心混凝土压溃时，

${I_1}{\rm{ = }}{f_{{\rm{cc}}}}{\rm{ - }}2{\sigma _{\rm{L}}},\;{J_2} = \frac{1}{3}{\left( f_{\rm cc}- {\sigma _{\rm{L}}} \right)^2}$

代入式（1），可得

(7) $\gamma {\rm{ = }}\left( {1 - \alpha } \right)k - \left( {1 - 2\alpha } \right).$

由Mohr-Coulomb准则^[29]可得

(8) $\frac{{1 + \sin \varphi }}{{1 - \sin \varphi }}{\sigma _{\max }} - {\sigma _{\min }} = {f_{\rm{c}}}.$

对均匀约束混凝土而言，σ_max=–σ_L，σ_min=–f_cc，代入式（8），联立式（6）、（7）可得

(9) $\gamma {\rm{ = }}\left( {1 - \alpha } \right)\frac{{1 + \sin \varphi }}{{1 - \sin \varphi }} - \left( {1 - 2\alpha } \right).$

内摩擦角φ可按文献[30]、[31]的方法计算：

(10) $\varphi {\rm{ = }}13.84^\circ \exp\; \left( { - 21.73{{\tilde \varepsilon }^{\rm{p}}}} \right) + 22.41^\circ .$

流动准则用于确定混凝土进入塑性状态后的塑性变形，是混凝土本构模型的基本组成部分之一. 本研究选用应用广泛的非关联流动准则——Drucker-Prager流动准则，其表达式如下：

(11) $G = \sqrt {{{\left( {\zeta {\sigma _{\rm{t}}}_0\tan \psi } \right)}^2} + 3{J_2}} + \frac{1}{3}{I_1}\tan \psi .$

式中：ζ为偏心率，通常取0.1；σ_t0为混凝土单轴抗拉强度；ψ为膨胀角.

准确分析混凝土的塑性变形，关键在于获取合理的膨胀角^[32]. Wang等^[33]研究了超高性能混凝土的塑性膨胀角与环向套箍应力之间的关系，但并未考虑材料损伤演化的影响，显著低估了构件的塑性变形. 本文在文献[33]的基础上，考虑材料损伤演化的影响，重新推导膨胀角的计算方法.

混凝土应变包括弹性应变（ ${ \varepsilon _{}}^{\rm{e}}$）和塑性应变（ ${\varepsilon _{}}^{\rm{p}}$）两部分：

(12) ${ \varepsilon _{}} = { \varepsilon _{}}^{\rm{e}} + { \varepsilon _{}}^{\rm{p}}.$

根据弹塑性力学基本理论，塑性应变增量表达式如下：

(13) ${\dot { \varepsilon} _{}}^{\rm{p}}{\rm{ = }}\dot \lambda \frac{{\partial G}}{{\partial {{ {\sigma}} _{}}}}.$

式中：λ为非负标量.

环向塑性应变（ε₁^p）、径向塑性应变（ε₂^p）和纵向塑性应变（ε₃^p）的表达式如下：

(14) ${\varepsilon _1}^{\rm{p}} = {\varepsilon _2}^{\rm{p}} = {\varepsilon _1} - \frac{1}{{{E_{\rm{c}}}\left( {1 - d} \right)}}\left( {{\sigma _1} - 2\upsilon {\sigma _3}} \right), $

(15) ${\varepsilon _3}^{\rm{p}} = {\varepsilon _3} - \frac{1}{{{E_{\rm{c}}}\left( {1 - d} \right)}}\left( {{\sigma _3} - 2\upsilon {\sigma _1}} \right).$

式中：d为损伤因子，E_c为等效塑性应变的函数，υ为混凝土弹性泊松比， ${{\rm{\sigma }}_{^1}} $、 $ {{\rm{\sigma }}_2}$和 ${{{\sigma }}_3} $分别为核心混凝土环向、径向和竖向应力，ε₁、ε₂和ε₃分别为核心混凝土环向、径向和竖向应变.

把式（11）代入式（13），可得3个纵向与环向塑性应变增量表达式如下：

(16) ${\dot \varepsilon _1}^{\rm{p}} = \dot \lambda \left[ {\frac{1}{3}\tan \psi - \frac{{{\sigma _2} - 2{\sigma _1} + {\sigma _3}}}{{2\sqrt {{{\left( {\zeta {\sigma _{{\rm{t}}0}}\tan \psi } \right)}^2} + 3{J_2}} }}} \right], $

(17) ${\dot \varepsilon _2}^{\rm{p}} = \dot \lambda \left[ {\frac{1}{3}\tan \psi - \frac{{{\sigma _1} - 2{\sigma _2} + {\sigma _3}}}{{2\sqrt {{{\left( {\zeta {\sigma _{{\rm{t}}0}}\tan \psi } \right)}^2} + 3{J_2}} }}} \right], $

(18) ${\dot \varepsilon _3}^{\rm{p}} = \dot \lambda \left[ {\frac{1}{3}\tan \psi - \frac{{{\sigma _1} - 2{\sigma _3} + {\sigma _2}}}{{2\sqrt {{{\left( {\zeta {\sigma _{{\rm{t}}0}}\tan \psi } \right)}^2} + 3{J_2}} }}} \right].$

因ζ通常取0.1^[26]， ${\left[ {{{\left( {\zeta {\sigma _{{\rm{t}}0}}\tan \psi } \right)}^2} + 3{J_2}} \right]^{1/2}} \approx \sqrt {3{J_2}} $，又有

(19) ${\dot \varepsilon _1}^{\rm{p}} + {\dot \varepsilon _2}^{\rm{p}} + {\dot \varepsilon _3}^{\rm{p}} = \dot \lambda \tan \psi ,$

(20) ${\dot \varepsilon _3}^{\rm{p}} - 0.5\left({\dot \varepsilon _1}^{\rm{p}} + {\dot \varepsilon _2}^{\rm{p}}\right) = - \frac{3}{2}\dot \lambda .$

联立式（19）、（20），又dε₁^p=dε₂^p，可得

(21) $\tan \;\psi {\rm{ = }} - \frac{3}{2}\frac{{{{\dot \varepsilon }_3}^{\rm{p}} + 2{{\dot \varepsilon }_1}^{\rm{p}}}}{{{{\dot \varepsilon }_3}^{\rm{p}} - {{\dot \varepsilon }_1}^{\rm{p}}}} = \frac{{\sqrt 3 }}{2}\dot I_1^{\prime {\rm{p}}}{\left( {{{\dot J}_{\rm{2}}}^{\rm{p}}} \right)^{{\rm{ - 1/2}}}},$

(22) ${\dot I_1^{'{\rm{p}}}} = {\dot \varepsilon _1}^{\rm{p}} + {\dot \varepsilon _2}^{\rm{p}} + {\dot \varepsilon _3}^{\rm{p}},$

(23) $ {\dot J_2^{'{\rm{p}}}} = \frac{1}{6}\left[ {{{\left( {{{\dot \varepsilon }_1}^{\rm{p}} - {{\dot \varepsilon }_2}^{\rm{p}}} \right)}^2} + {{\left( {{{\dot \varepsilon }_1}^{\rm{p}} - {{\dot \varepsilon }_3}^{\rm{p}}} \right)}^2}} \right. + \left. {{{\left( {{{\dot \varepsilon }_2}^{\rm{p}} - {{\dot \varepsilon }_3}^{\rm{p}}} \right)}^2}} \right]. $

当tan ψ<0时，表示核心混凝土发生塑性收缩；当tanψ>0时，表示核心混凝土发生塑性膨胀.

在无可靠试验数据的前提下，可采用文献[32]建议的塑性膨胀角计算公式：

(24) $\psi {\rm{ = }}\left\{ \begin{array}{l} 56.3\left( {1 - {\xi _{\rm{c}}}} \right),\qquad\quad\;\;\;\;{\xi _{\rm{c}}} \leqslant 0.5;\\ 6.672{{\rm exp}{\;\left(\displaystyle\frac{{7.4}}{{4.64 + {\xi _{\rm{c}}}}}\right)\;}},\;{\xi _{\rm{c}}} > 0.5. \end{array} \right.$

式中：ξ_c为套箍系数，

(25) ${\xi _{\rm{c}}}{\rm{ = }}\frac{{{A_{\rm{s}}}{f_{\rm{y}}}}}{{{A_{\rm{c}}}{f_{{\rm{ck}}}}}}.$

其中，A_s、A_c分别为钢管与核心混凝土面积，f_y、f_ck分别为钢管屈服强度和核心混凝土圆柱体抗压强度.

损伤演化和硬化准则用以描述混凝土加载过程中的刚度退化与强度软化，是混凝土本构模型的另一基本组成部分.

根据文献[26]，引入损伤因子d，应力张量σ与应变张量ε的关系可表达如下：

(26) ${{\sigma }} = \left( {1 - d} \right){{{D}}_0}^{{\rm{el}}}:\left( {{{\varepsilon }} - {{{\varepsilon }}^{\rm{p}}}} \right).$

式中： ${{{D}}_0}^{{\rm{el}}}$ 为混凝土初始弹性矩阵. 损伤因子d表达式如下：

(27) ${1 - d} {\rm{ = }}\left( {1 - {s_{\rm{t}}}{d_{\rm{c}}}} \right)\left( {1 - {s_{\rm{c}}}{d_{\rm{t}}}} \right).$

式中：d_t和d_c分别表示混凝土拉伸损伤和压缩损伤，是等效拉伸塑性应变（ ${\tilde \varepsilon _{\rm{t}}}^{\rm{p}}$）和等效压缩塑性应变（ ${\tilde \varepsilon _{\rm{c}}}^{\rm{p}}$）的函数，

(28) ${d_{\rm{t}}}{\rm{ = }}{d_{\rm{t}}}\left( {{{\tilde \varepsilon }_{\rm{t}}}^{\rm{p}}} \right),\;{d_{\rm{c}}}{\rm{ = }}{d_{\rm{c}}}\left( {{{\tilde \varepsilon }_{\rm{c}}}^{\rm{p}}} \right),$

其中，

${{{\tilde{{\textit{ε}}}}}^{\rm{p}}}{\rm{ = }}{\left[ {{{\tilde \varepsilon }_{\rm{t}}}^{\rm{p}}\; {{\tilde \varepsilon }_{\rm{c}}}^{\rm{p}}} \right]^{\rm{T}}},$ ${\dot {\tilde \varepsilon} _{\rm{t}}}^{\rm{p}}{\rm{ = }}r\left( {{\hat{ \sigma }}} \right){\dot {\tilde \varepsilon}^{\rm{p}}}_{\max },$ ${\dot {\tilde \varepsilon} _{\rm{c}}}^{\rm{p}}{\rm{ = }} - \left( {1 - r\left( {{\hat{ \sigma }}} \right)} \right){\dot {\hat \varepsilon} ^{\rm{p}}}_{\min };$

s_t和s_c分别为损伤因子恢复函数，用以实现应力反向时的刚度恢复：

(29) $\left\{ \begin{array}{l} {s_{\rm{t}}} = 1 - {\omega _{\rm{t}}}r\left( {{\hat{ \sigma }}} \right),\;\;\;\;\;\quad\quad0 \leqslant {\omega _{\rm{t}}} \leqslant 1;\\ {s_{\rm{c}}} = 1 - {\omega _{\rm{c}}}\left[ {1{\rm{ - }}r\left( {{\hat{ \sigma }}} \right)} \right],\quad0 \leqslant {\omega _{\rm{c}}} \leqslant 1. \end{array} \right.$

(30) $ r\left( {{\hat{ \sigma }}} \right) = \frac{{\displaystyle\sum\nolimits_{i = 1}^3 {\left\langle \,{{{\hat {\bar \sigma} }_i}} \right\rangle } }}{{\displaystyle\sum\nolimits_{i = 1}^3 {\left| \,{{{\hat {\bar \sigma }}_i}} \right|} }},\quad\;\;0 \leqslant r\left( {{\hat{ \sigma }}} \right) \leqslant 1. $

其中， ${\hat{ \bar \sigma} _i}$ 为等效应力张量的主应力，w_t、w_c分别为拉伸和压缩刚度恢复权重因子.

由于0>σ₁=σ₂>σ₃，由式（30）可得， $r\left( {{\hat{ \sigma }}} \right)$=0；另外，根据文献[26]，通常 ${s_{\rm{t}}} = 1$， ${s_{\rm{c}}} = 0$；代入式（27）可得

(31) $\left( {1 - d} \right){\rm{ = }}1 - {d_{\rm{c}}}\left( {{{\tilde \varepsilon }_{\rm{c}}}^{\rm{p}}} \right).$

式中：d_c的计算方法可参见文献[34]. 根据文献[26]，混凝土的屈服面演化规律（强化准则），依然通过等效塑性应变 ${{\tilde{ \varepsilon }}^{{\rm p}}}$ 驱动：

(32) $d{\tilde \varepsilon ^{\rm{p}}}{\rm{ = }}d{\tilde \varepsilon _{\rm{c}}}^{\rm{p}} = - d{\hat \varepsilon ^{{\rm{p}}}}_{\min }.$

在确定了均匀约束条件下的混凝土塑性损伤本构关系后，还需要建立相匹配的非线性徐变分析模型.

目前关于多轴不同应力水平作用下的徐变规律研究还不成熟，没有直接可供借鉴的预测模型，本文拟通过引入徐变泊松比的概念，把单轴徐变规律拓展到多轴应力状态下，文献[20]、[34]采用了类似的方法.

单轴应力下的徐变泊松比是指在轴向应力作用下，横向徐变与轴向徐变之比. 当在多轴应力状态下，不同主应力方向的泊松比μ_{cp, i}还与应力状态有关^[35-37]，可表达为

(33) ${\mu _{{\rm{cp,}}}}_i = 0.160 - 0.074\frac{{{\sigma _i}}}{{{\sigma _j} + {\sigma _k}}} + 0.028{\left( {\frac{{{\sigma _i}}}{{{\sigma _j} + {\sigma _k}}}} \right)^2}.$

式中: σ_i为第i个主应力 $ ;i,j,k = 1,2,3;i \ne j \ne k$. 考虑徐变泊松比影响的转换矩阵A可表达

(34) ${{A}} = \left[ {\begin{array}{*{20}{c}} 1&{ - {\mu _{{\rm{cp}},2}}}&{ - {\mu _{{\rm{cp}},3}}}&0&0&0 \\ { - {\mu _{{\rm{cp}},1}}}&1&{ - {\mu _{{\rm{cp}},3}}}&0&0&0 \\ { - {\mu _{{\rm{cp}},1}}}&{ - {\mu _{{\rm{cp}},2}}}&1&0&0&0 \\ 0&0&0&{2 + {\mu _{{\rm{cp}},2}} + \;{\mu _{{\rm{cp}},3}}}&0&0 \\ 0&0&0&0&{2 + {\mu _{{\rm{cp}},1}} + {\mu _{{\rm{cp}},3}}}&0 \\ 0&0&0&0&0&{2 + {\mu _{{\rm{cp}},1}} + \;{\mu _{{\rm{cp}},2}}} \end{array}} \right].$

基于连续损伤介质力学，通过引入混凝土材料损伤因子，三维非线性徐变本构关系可表示为

(35) $\begin{split} {{\varepsilon }}(t) =& J\left( {t,\tau ,{{\sigma }}} \right){{A}} \cdot {{\sigma }}\left( \tau \right) + \int_\tau ^t {J\left( {t,\tau ,{{\sigma }}} \right){{A}} \cdot \displaystyle\frac{{{\rm d}{{\sigma }}(\tau )}}{{{\rm d}\tau }}\rm d} \tau + \\ &{{{\varepsilon }}_{\rm{s}}}(t) + {{{\varepsilon }}_{{\rm{tem}}}}(t). \end{split}$

式中：τ为混凝土加载龄期；ε_s为收缩应变；ε_tem为温度应变；J(t，τ，σ)为徐变函数，其表达式如下：

(36) $J\left( {t,\tau ,{{\sigma }}} \right) = C\left( {t,\tau ,{{\sigma }}} \right) + \frac{1}{{E\left( \tau \right)\left( {1 - d\left( \tau \right)} \right)}}.$

式中：E(τ)为混凝土在加载时刻的弹性模量；C(t，τ)为徐变度函数，需转换成用Dirichlet级数表达，以避免对全部历程的应力历史积分，如式（37）所示；d(τ)为损伤因子，其值在[0，1.0].

(37) $C\left( {t,\tau ,{{\sigma }}} \right) = \frac{1}{{\left( {1 - d\left( \tau \right)} \right)}}C\left( {t,\tau } \right).$

(38) $C\left( {t,\tau } \right){\rm{ = }}\sum\limits_{i = 1}^m \;{\bigg\{ {{a_i}(\tau )\Big[ {1 - {\rm{exp}}\;\big[ { - {\lambda _{\rm{i}}}\left( {t - \tau } \right)} \big]} \Big]} \bigg\}} .$

参数a_i(τ)与λ_i可按如下步骤^[21]确定：

1）选择合适的m与λ_i，一般情况取m=4已达到足够的精度，λ₁=1，λ₂=0.1，λ₃=0.01，λ₄=0.001；

2）确定混凝土的养护龄期τ₀和加载龄期τ；

3）如需计算t时刻的徐变效应，需将t−τ分成n段， $\left( {{t_1},{t_2},{\rm{ }} \cdots ,{t_i}, \cdots ,{t_n}\;,{t_i} > {t_0},{t_n} = t} \right)$ 确定计算时间点t_i，i=1，2，· · · ，n，t_i≥τ₀；

4）在不同加载的加载龄期，先根据相应的徐变预测模型计算C(t_i , τ)，再代入式（38），可得

(39) ${{{N}}_{m \times n}} \times {{{a}}_{m \times 1}} = {{{b}}_{n \times 1}}.$

(40) ${{{N}}_{m \times n}} = \left[ { \begin{array}{*{20}{c}} {1 - \exp \left( { - {\lambda _1}({t_1} - \tau )} \right)}&{1 - \exp \left( { - {\lambda _{\rm{2}}}({t_1} - \tau )} \right)}&{...}&{1 - \exp \left( { - {\lambda _m}({t_1} - \tau )} \right)} \\ {1 - \exp \left( { - {\lambda _1}({t_{\rm{2}}} - \tau )} \right)}&{1 - \exp \left( { - {\lambda _2}({t_2} - \tau )} \right)}&{...}&{1 - \exp \left( { - {\lambda _m}({t_2} - \tau )} \right)} \\ \vdots & \vdots &\; & \vdots \\ {1 - \exp \left( { - {\lambda _1}({t_n} - \tau )} \right)}&{1 - \exp \left( { - {\lambda _2}({t_n} - \tau )} \right)}&{...}&{1 - \exp \left( { - {\lambda _m}({t_n} - \tau )} \right)} \end{array}} \right],$

(41) $ \begin{array}{l} \\[-5pt] {{ a}_{m \times 1}} = \left[ {\begin{array}{*{20}{l}} {{a_1}\left( \tau \right)}\\ {{a_2}\left( \tau \right)}\\ {\;\;\;\; \vdots }\\ {{a_m}(\tau )} \end{array}} \right],\;{{ b}_{n \times 1}} = \left[ {\begin{array}{*{20}{l}} {C({t_1},\tau )}\\ {C({t_2},\tau )}\\ {\;\;\;\;\;\; \vdots }\\ {C({t_n},\tau )} \end{array}} \right]. \end{array}$

5）由最小二乘法，可求得待定系数为

(42) ${{{a}}_{m \times 1}} = {\left( {{{{N}}^{\rm{T}}}{{N}}} \right)^{ - 1}}\!\!\!\!\!\!\!\!_{m \times m}{\left( {{{{N}}^{\rm{T}}}{{b}}} \right)_{m \times 1}}.$

混凝土的收缩徐变会使CFST柱产生内力重分配. 在众多的分析方法中，逐步递推分析方法能够较为方便的处理这一复杂过程，且能与现有通用有限元程序相结合^[21-22]，其关键在于任意时间步Δt_i徐变应变增量的计算.

由式（35）可知，t_n、t_n+1时刻以及Δt_n+1内的徐变应变ε_c可表示为如式（43）所示，Δt_n+1内的徐变应变增量如式（44）所示.

(43) $\left. \begin{aligned} &\begin{aligned} {{{\varepsilon }}_{\rm{c}}}({t_n}) = {{A}}\left[ {\Delta {{{\sigma }}_0}C\left( {{t_n},{t_0},{{{\sigma }}_{{t_0}}}} \right)} \right. + \cdots +\\ \;\left. {C\left( {{t_n},{t_{n - 1}},{{{\sigma }}_{{t_{n - 1}}}}} \right)\Delta {{{\sigma }}_{n - 1}}} \right],\;\;\;\;\end{aligned}\\ &\begin{aligned} {{{\varepsilon }}_{\rm{c}}}({t_{n + 1}}) =& {{A}}\left[ {\Delta {{{\sigma }}_0}C\left( {{t_{n + 1}},{t_0},{{{\sigma }}_{{t_0}}}} \right)} \right. + \cdots +\\ & \left. {C\left( {{t_{n + 1}},{t_n},{{{\sigma }}_{{t_n}}}} \right)\Delta {{{\sigma }}_n}} \right].\;\;\;\;\;\;\end{aligned} \end{aligned} \right\}$

(44) $\begin{split} &\Delta {{{\varepsilon }}_{\rm{c}}}({t_{n + 1}}) = \left[ {C\left( {{t_{n + 1}},{t_0},{{{\sigma }}_{{t_0}}}} \right) - C\left( {{t_n},{t_0},{{{\sigma }}_{{t_0}}}} \right)} \right]{{A}}\Delta {{{\sigma }}_0} + \cdots +\\ & \quad \left[ {C\left( {{t_{n + 1}},{t_{n - 1}},{{{\sigma }}_{{t_{n - 1}}}}} \right) - C\left( {{t_n},{t_{n - 1}},{{{\sigma }}_{{t_{n - 1}}}}} \right)} \right]\;{{A}}\Delta {{{\sigma }}_{n - 1}} + \\ &\quad\left[ {C\left( {{t_{n + 1}},{t_n},{{{\sigma }}_{{t_n}}}} \right)} \right]{{A}}\Delta {{{\sigma }}_n}. \end{split}$

式中：Δσ_n为Δt_n增量步内的应力张量增量.

(45) $\begin{array}{l} \!\!\!\!C\left( {{t_{n + 1}},{t_j},{\sigma _j}} \right) - C\left( {{t_n},{t_j},{\sigma _j}} \right) = \;{\left( {1 - d( {{t_j}} )} \right)^{ - {\rm{1}}}} \times \;\\ \!\!\!\!\displaystyle\sum\limits_{i = 1}^m {} \bigg\{ {{a_i}({t_j}){\rm{exp}}\;\;\left[ { - {\lambda _i}({t_n} - {t_j})} \right]\Big[ {1 - {\rm{exp}}\;( - {\lambda _i}\Delta {t_{n + 1}})} \Big]} \bigg\}. \end{array} $

将式（45）代入式（43），可得Δt_n时刻的徐变应变增量的显式表达式：

(46) $\Delta \varepsilon {\left( {{t_n}} \right)^c} = \sum\limits_{i = 1}^m \; \bigg\{ {{ \omega} \left( {i,n - 1} \right)\Big[ {1 - {\rm{exp}}\;\left( { - {\lambda _i}\Delta {t_n}} \right)} \Big]} \bigg\}.$

式中：ω(i，n−1)可由下式递推关系求得,

(47) $\left. \begin{aligned} &\begin{aligned} {{\omega }}\left( {i,n} \right) =& {\left( {1 - d\left( {{t_{n - 1}}} \right)} \right)^{ - {\rm{1}}}}{{\omega }}\left( {i,n - 1} \right)\exp \;\left( { - {\lambda _i}\Delta {t_{n - 1}}} \right) + \\ &{\left( {1 - d\left( {{t_n}} \right)} \right)^{ - {\rm{1}}}}{{A}}\Delta {{{\sigma }}_n}{a_i}\left( {{t_n}} \right),\end{aligned} \\ &\begin{aligned}{{\omega }}\left( {i,1} \right) = {\left( {1 - d\left( {{t_0}} \right)} \right)^{ - {\rm{1}}}}{{A}}\Delta {{{\sigma }}_0}{a_i}\left( {{t_0}} \right). \end{aligned}\end{aligned} \right\}$

根据前文提出的三维非线性徐变本构关系，推导任意时间步徐变增量的显式表达式，再将其转换为附加应变荷载，施加在结构体系上；同时结合已经建立的均匀约束条件下的混凝土塑性损伤本构模型，实现徐变过程中的塑性修正和损伤修正. 根据前述的混凝土三维收缩徐变模型，结合ABAQUS软件强大的二次开发功能，通过各子程序协同工作，即可较为方便地实现前述全部分析过程. 各子程序的功能和相互调用过程大致如下：用户子程序USDFLD定义预定义场变量，可以实现混凝土弹性模量随时间的变化；功能子程序GETVRM通过提取每一时刻材料点变量值，包括应力和应变值，并将其传递到子程序UEXPAN中；用户子程序UEXPAN通过自编算法实现每一时步收缩、徐变应变增量的计算，并将其返回到主程序中；在主程序中引入混凝土塑性损伤本构，进行塑性和损伤修正.

付学宝等^{[9, 38]}分别进行了不同应力水平作用下的轴压圆CFST短柱徐变试验，为便于描述，简记为试验I与试验II. 各试件尺寸、加载荷载、应力水平等如表1所示，D为圆柱体试件外径；T为钢管厚度；h为时间高度；F为试验荷载大小；n_c表示混凝土应力水平，为核心混凝土的初始应力和圆柱体抗压强度平均值之比.

核心混凝土和钢管材料参数如表2所示. f_ck表示Φ150×300 mm圆柱体试块抗压强度标准值；f_y表示钢管屈服强度；E_c代表混凝土弹性模量，E_s代表钢管弹性模量. 试验在试验室内进行，室内温度相对恒定. 应变测点均位于试件中部. 试验加载方案如图3所示，均通过千斤顶与力传感器配合加载，当徐变引起的千斤顶卸载超过1%时进行补载. 另外，在文献[38]中，为避免杠杆加载引起荷载偏心，在试件底面预垫了薄钢片.

收缩徐变预测模型均选用Bazant教授提出的B4模型，一是因为该模型考虑因素较为全面，具有较高的精度^[39]；二是因为该模型把徐变分为干燥徐变与基本徐变，收缩分为干燥收缩与自收缩，而干燥徐变与干燥收缩有关，当不计干燥徐变与干燥收缩时，能够契合CFST基本不与外界发生湿度交换的特殊使用环境. 同时，对于自收缩变形的计算，因为混凝土自生水化反应会消耗一部分水分，混凝土湿度取85%^[40]. 值得注意的是，数据观测是从28 d龄期开始，而此时自收缩变形已基本完成，因此模型计算结果要剔除自收缩引起的徐变变形后，才能与试验数据对比.

本文所涉及到的参数计算结果如表3所示. 参数q₂~q₄的取值与混凝土配合比、养护温度等相关，计算方法详见文献[24].

选用ABAQUS软件对试件进行有限元模拟，混凝土和钢管分别采用实体单元C3D8、C3D8I，钢管与混凝土交界面法向采用“硬接触”模拟，切向采用摩擦罚函数模拟，界面摩擦系数取0.25^[41]. 假定垫板为不变形刚体，垫板与钢和混凝土的连结采用绑定约束. 由于无该试件受压全过程荷载位移曲线实测数据，膨胀角按式（24）进行计算，试验I与试验II试件的膨胀角分别为27.5°与34.1°，其他参数计算参见文献[34]. 有限元网格模型如图4所示. 文献[9]与[38]的实测应变与分析结果对比如图5所示，图中T表示试验组.

分析计算结果可知，对核心混凝土应力水平较小时（<0.6），采用线性徐变理论（L）与非线性徐变理论（NL）的计算结果差别不大，且均与试验值较为接近. 但对于应力水平较大的试件而言，采用线性与非线性徐变理论的计算结果差异明显，在加载110 d后，对试验I与试验II而言，采用非线性徐变理论的计算结果分别较线性徐变理论大约84%和72%。

为准确分析轴压圆CFST短柱非线性徐变效应，本文提出了一种新的三维非线性徐变效应分析理论框架，并建立了相应的数值分析方法. 将不同应力水平下的徐变试验结果与分析结果进行对比，结果表明：在高水平应力作用下，徐变与混凝土塑性损伤演化之间存在耦合特性，采用线性徐变理论会显著低估构件的长期徐变变形；当采用非线性徐变理论时，徐变应变计算值与试验值更为接近，且变化趋势也与试验值保持良好的一致性，验证了所提算法的可靠性.

本文所提方法主要适用于轴心受压圆CFST短柱的徐变效应分析，对于轴心受压长柱或偏心受压长短柱而言，由于核心混凝土不再满足均匀约束条件，其非线性徐变效应分析方法还有待进一步研究.