CN110162827A - 一种混凝土结构时变效应的实体有限元计算方法 - Google Patents

Abstract

本发明涉及一种混凝土结构时变效应的实体有限元计算方法，包括如下步骤：1）依据逐步法将混凝土的应力历程离散化，推导出第n个时段内的徐变应变增量；2）运用Matlab软件将徐变系数表达式拟合成指数级数形式；3）将徐变系数拟合表达式代入徐变应变增量计算式，建立混凝土单轴徐变应变增量递推计算式；4）考虑徐变泊松效应，在实体单元材料本构中建立混凝土时变本构方程；5）运用Fortran语言在ANSYS软件的USERMAT子程序中开发混凝土的时变本构。与现有技术相比，本发明能实现空间应力状态下混凝土结构在复杂应力历程中的时变效应数值预测，与常用设计规范的徐变与收缩模型衔接良好。

Description

一种混凝土结构时变效应的实体有限元计算方法

技术领域

本发明涉及混凝土的徐变与收缩时变效应数值预测领域，特别是涉及一种混凝土结构时变效应的实体有限元计算方法。

背景技术

混凝土的徐变与收缩特性对混凝土桥梁及其它暴露于干燥环境下的大体积混凝土结构的长期性能有显著影响。以大跨度预应力混凝土桥梁为例，混凝土的时变特性使桥梁随服役年限增长而出现跨中挠度增加，预应力损失增大，混凝土抗裂安全储备降低，影响桥梁结构运营安全及耐久性。有必要提高混凝土结构的时变效应预测精度，以减小桥梁结构运营期内的损伤与病害，降低后期维护成本，避免结构破坏或垮塌等严重后果。

工程结构分析多采用杆系有限元程序，使用设计规范推荐的混凝土时变本构特性及徐变效应线性叠加算法。但杆系有限元程序不能充分考虑混凝土结构D区的受力特点，例如箱梁剪力滞效应，钢与混凝土间的应力重分布等。有必要开发适用于空间应力状态下混凝土结构时变效应分析的实体单元，以提高构造或受力复杂的混凝土结构时变效应的预测精度。

使用通用有限元软件开展混凝土时变效应实体有限元分析，仍至少存在3处困难。首先，实体单元徐变理论尚不完善，对于多轴应力工况，徐变泊松效应的影响尚不明确。其次，与工程实践所用的通用混凝土时变本构衔接难度大。最后，对于复杂应力历程和复杂应力状态下混凝土结构时变效应的预测精度有待进一步提高。

发明内容

为克服上述现有技术的不足，本发明提供了一种混凝土结构时变效应的实体有限元计算方法，以解决多轴应力状态和复杂应力历程下，混凝土结构时变效应预测精度不高，与设计规范的徐变与收缩时变本构衔接不足等技术难题。

为实现上述发明目的，本发明通过以下技术方案来实现。

一种混凝土结构时变效应的实体有限元计算方法，包括以下步骤：

步骤a：采用逐步法将混凝土的应力历程离散化，推导出第n个时段内的单轴徐变应变增量。

步骤b：采取线性最小二乘法初步拟合，非线性规划求极值法最终拟合的策略，将设计规范推荐的徐变系数表达式拟合成指数级数形式。

步骤c：将徐变系数拟合表达式代入单轴徐变应变增量计算式，建立混凝土单轴徐变应变增量递推计算式。

步骤d：考虑多轴应力工况的徐变泊松影响，在实体单元材料本构中引入混凝土时变效应递推计算式，建立混凝土实体单元时变本构方程。

步骤e：运用Fortran语言在ANSYS软件的USERMAT子程序中开发混凝土时变本构，实现空间应力状态下混凝土结构在复杂应力历程中的时变效应有限元分析。

所述的步骤a中，当结构截面应力连续变化时，将应力历程划分为若干个相互衔接的小过程，在每一个小过程内假定应力为常量，并以该应力状态发生徐变，则第n个时段内的单轴徐变应变增量如式(2.1)所示。

式中，Δε_c(t_n,t_n-1)为混凝土在t_n-1至t_n时段内的单轴徐变应变增量；Δσ_i为t_i时刻的应力增量；E(t_i)为t_i时刻混凝土的弹性模量；为加载龄期t_i的混凝土在t_n时刻的徐变系数；为加载龄期t_i的混凝土在t_n-1时刻的徐变系数。

所述的步骤b中，以中国设计规范(JTG 3362-2018)为例，混凝土的徐变系数表达式如式(2.2)所示。

式中，为加载龄期t_i的混凝土在t时刻的徐变系数；为混凝土在t_i时刻的名义徐变系数；β_H为与环境相对湿度及构件理论厚度相关的系数。

采取线性最小二乘法初步拟合，等式约束条件下非线性规划求极值最终拟合的策略，运用Matlab软件将式(2.2)数值拟合为如式(2.3)所示的指数级数形式。

式中，C_j、q_j(j∈[1,m])为拟合常数；m为展开的级数项数，一般m≥3有足够的拟合精度；其余符号意义同前。对于欧洲fib Mode Code 2010规范、美国AASHTO LRFD桥梁规范等亦能拟合为式(2.3)形式。

所述的步骤c中徐变应变增量递推计算式的推导过程如下：

将展开后的徐变系数级数表达式(2.3)代入式(2.1)，利用指数函数特点进一步合并，如式(2.4)所示。

式中，Δt_n为第n个时段的时间增量，即Δt_n＝t_n-t_n-1；其余符号意义同前。

对式(2.4)的累加项进行重组，并令：

式中，为递推变量；其余符号意义同前。

于是有

同样，第n+1个时段的徐变应变增量也可用式(2.7)表示。

其中，

比较式(2.5)和式(2.8)，可得到如式(2.9)所示的递推关系。

混凝土收缩应变与应力无关，当构件尺寸、材料特性、环境状况和收缩开始龄期已知时，收缩应变随时间发展曲线唯一确定。与徐变相似，将收缩应变曲线离散化，则第n个时段的收缩应变增量Δε_s(t_n,t_n-1)如式(2.10)所示。

Δε_s(t_n,t_n-1)＝ε_s(t_n,t_s)-ε_s(t_n-1,t_s) (2.10)

式中，t_s为混凝土收缩开始龄期；ε_s(t_n,t_s)为混凝土从t_s至t_n时段的收缩应变；ε_s(t_n-1,t_s)为混凝土从t_s至t_n-1时段的收缩应变。

利用上述递推关系，第n个时段的徐变与收缩应变增量Δε(t_n,t_n-1)计算如式(2.11)所示。

数值分析时，采用该递推计算式不必存储应力历程，大大节约了计算资源和时间消耗。

所述的步骤d中在实体单元建立时变本构方程的过程为：

设第n个时段内实体单元的应力增量矩阵Δσ_n与应变增量矩阵Δε_n如式(2.12)所示：

式(2.12)中单元的应变增量矩阵Δε_n表示为式(2.13)所示的三部分：

Δε_n＝Δε_e,n+Δε_c,n+Δε_s,n (2.13)

式中，Δε_e,n为弹性应变增量矩阵；Δε_c,n为徐变应变增量矩阵；Δε_s,n为收缩应变增量矩阵。

实体单元弹性应变增量矩阵如式(2.14)和(2.15)所示：

其中，

式中，Δε_e,n为单元在第n个时段的弹性应变增量矩阵；E(t_n-1)为材料在t_n-1时刻的弹性模量；A为弹性泊松比矩阵；Δσ_n为单元在第n个时段的应力增量矩阵；ν为材料的弹性应变泊松比。

实体单元徐变应变增量Δε_c,n如式(2.16)所示：

其中，

式中，为与应力分量对应的向量，C为徐变泊松比矩阵，ν’为材料的徐变泊松比；其余符号意义同前。

混凝土单元收缩应变增量矩阵Δε_s,n如式(2.19)所示：

Δε_s,n＝ε_s,n-ε_s,n-1 (2.19)

其中，

ε_s,n＝[ε_s,n ε_s,n ε_s,n 0 0 0]^T (2.20)

式中，ε_s,n为单元在第n个时段的收缩应变矩阵；ε_s,n-1为单元在第n-1个时段的收缩应变矩阵；其余符号意义同前。

由式(2.13)和(2.14)，应力增量矩阵可按式(2.21)计算。

Δσ_n＝D_nΔε_e,n＝D_n(Δε_n-Δε_c,n-Δε_s,n) (2.21)

其中，

D_n＝E(t_n-1)A^-1 (2.22)

式中，D_n为单元材料本构关系矩阵；其余符号意义同前。

把单元应变增量与节点位移增量关系式和式(2.16)代入式(2.21)，得到单元应力增量矩阵，如式(2.23)所示：

其中，

式中，B为单元几何矩阵；为单元在第n个时段的弹性变形增量矩阵；α_n为单元在第n个时段的徐变应变增量矩阵。

有限元法的平衡方程组可表示为式(2.25)：

式中，为单元外荷载增量矩阵。

将式(2.23)代入式(2.25)，得到徐变与收缩效应分析的基本方程如式(2.26)所示：

其中，

式中，k_n为单元刚度矩阵；为徐变变形产生的当量荷载增量矩阵；为收缩变形产生的当量荷载增量矩阵。

通过坐标转换，按节点序号将单元刚度矩阵、荷载矩阵集成结构刚度矩阵和荷载矩阵，形成结构平衡方程并可求解节点位移，然后将节点位移坐标转换后再由式(2.23)计算单元应力增量矩阵。

所述的步骤e中，在USERMAT子程序中开发混凝土徐变与收缩本构模型的主要步骤如下：

步骤e1：输入徐变系数、收缩应变计算参数(环境相对湿度、构件理论厚度、混凝土强度等级等等)，弹性模量随龄期变化函数，加载龄期t₀及收缩龄期t_s等数据；

步骤e2：主程序从时间t_n-1＝0开始计算，令t_n-1小于等于0.01天时，只计算弹性应变，得到初始应力增量并存储在状态变量数组中，为后面求解徐变应变增量提供应力计算参数；

步骤e3：当t_n-1大于0.01天时，根据式(2.11)建立徐变递推计算式，混凝土加载龄期在计算式中考虑，计算当前时间步列阵，并将列阵存储于状态变量数组中，由式(2.16)计算当前时间步的徐变应变增量Δε_c,n。

步骤e4：根据规范收缩模型计算当前时间的收缩应变ε_s,n，由式(2.19)、(2.20)计算当前时间步收缩应变增量列阵Δε_s,n，并将列阵ε_s,n存储于状态变量数组中。

步骤e5：根据式(2.15)、(2.22)由当前时间步的弹性模量和泊松比计算混凝土的本构矩阵D_n，在进行结构徐变收缩效应分析时，混凝土按线弹性材料处理，对于线弹性材料，一致切线算子矩阵D_ep与D_n相同。

步骤e6：根据式(2.21)计算混凝土应力增量列阵Δσ_n，并存储于状态变量数组，在时间步结束时更新单元应力，完成结构的徐变、收缩效应分析。

与现有技术相比，本发明具有以下优点：

1、本发明依据逐步递推算法建立适用于电算的实体单元时变本构方程，克服通用有限元软件分析混凝土时变效应难度大的劣势。与传统杆系有限元程序相比，能精细计算混凝土实体结构在空间应力状态和变化应力历程下的时变效应。

2、本发明与土木工程领域常用设计规范推荐的徐变与收缩模型能无缝衔接，而以往徐变递推算法采用的徐变模型仅针对水工大体积混凝土结构，本发明适用性更强。

3、本发明考虑了多轴应力工况下徐变泊松效应的影响，对于预测复杂应力状态下混凝土的时变效应，精度更高。

附图说明

图1为计算时段应力时程曲线；

图2为USERMAT子程序开发流程；

图3为圆柱体试件有限元模型；

图4为徐变系数拟合曲线与原始曲线比较。

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。本实施例以本发明技术方案为前提进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。

实施例

一种混凝土结构时变效应的实体有限元计算方法，以解决多轴应力状态和复杂应力历程下，混凝土结构时变效应预测精度不高，与设计规范的徐变与收缩时变本构衔接不足等技术难题。包括以下步骤：

步骤a：采用逐步法将混凝土的应力历程离散化，推导出第n个时段内的单轴徐变应变增量。

如图1，当结构截面应力连续变化时，将应力历程划分为若干个相互衔接的小过程，在每一个小过程内假定应力为常量，并以该应力状态发生徐变，则第n个时段内的单轴徐变应变增量如式(3.1)所示。

式中，Δε_c(t_n,t_n-1)为混凝土在t_n-1至t_n时段内的单轴徐变应变增量；Δσ_i为t_i时刻的应力增量；E(t_i)为t_i时刻混凝土的弹性模量；为加载龄期t_i的混凝土在t_n时刻的徐变系数；为加载龄期t_i的混凝土在t_n-1时刻的徐变系数。

步骤b：采取线性最小二乘法初步拟合，非线性规划求极值法最终拟合的策略，将设计规范推荐的徐变系数表达式拟合成指数级数形式。

以中国设计规范(JTG 3362-2018)为例，混凝土的徐变系数表达式如式(3.2)所示。

式中，为加载龄期t_i的混凝土在t时刻的徐变系数；为混凝土在t_i时刻的名义徐变系数；β_H为与环境相对湿度及构件理论厚度相关的系数。

采取线性最小二乘法初步拟合，等式约束条件下非线性规划求极值最终拟合的策略，运用Matlab软件将式(3.2)数值拟合为如式(3.3)所示的指数级数形式。

式中，C_j、q_j(j∈[1,m])为拟合常数；m为展开的级数项数，一般m≥3有足够的拟合精度；其余符号意义同前。对于欧洲fib Mode Code 2010规范、美国AASHTO LRFD桥梁规范等亦能拟合为式(3.3)形式。

步骤c：将徐变系数拟合表达式代入单轴徐变应变增量计算式，建立混凝土单轴徐变应变增量递推计算式。

将展开后的徐变系数级数表达式(3.3)代入式(3.1)，利用指数函数特点进一步合并，如式(3.4)所示。

式中，Δt_n为第n个时段的时间增量，即Δt_n＝t_n-t_n-1；其余符号意义同前。

对式(3.4)的累加项进行重组，并令：

式中，为递推变量；其余符号意义同前。

于是有

同样，第n+1个时段的徐变应变增量也可用式(3.7)表示。

其中，

比较式(3.5)和式(3.8)，可得到如式(3.9)所示的递推关系。

混凝土收缩应变与应力无关，当构件尺寸、材料特性、环境状况和收缩开始龄期已知时，收缩应变随时间发展曲线唯一确定。与徐变相似，将收缩应变曲线离散化，则第n个时段的收缩应变增量Δε_s(t_n,t_n-1)如式(3.10)所示。

Δε_s(t_n,t_n-1)＝ε_s(t_n,t_s)-ε_s(t_n-1,t_s) (3.10)

式中，t_s为混凝土收缩开始龄期；ε_s(t_n,t_s)为混凝土从t_s至t_n时段的收缩应变；ε_s(t_n-1,t_s)为混凝土从t_s至t_n-1时段的收缩应变。

利用上述递推关系，第n个时段的徐变与收缩应变增量Δε(t_n,t_n-1)计算如式(3.11)所示。

数值分析时，采用该递推计算式不必存储应力历程，大大节约了计算资源和时间消耗。

步骤d：考虑多轴应力工况的徐变泊松影响，在实体单元材料本构中引入混凝土时变效应递推计算式，建立混凝土实体单元时变本构方程。

设第n个时段内实体单元的应力增量矩阵Δσ_n与应变增量矩阵Δε_n如式(3.12)所示：

式(3.12)中单元的应变增量矩阵Δε_n可表示为式(3.13)所示的三部分：

Δε_n＝Δε_e,n+Δε_c,n+Δε_s,n (3.13)

式中，Δε_e,n为弹性应变增量矩阵；Δε_c,n为徐变应变增量矩阵；Δε_s,n为收缩应变增量矩阵。

实体单元弹性应变增量矩阵如式(3.14)和(3.15)所示：

其中，

式中，Δε_e,n为单元在第n个时段的弹性应变增量矩阵；E(t_n-1)为材料在t_n-1时刻的弹性模量；A为弹性泊松比矩阵；Δσ_n为单元在第n个时段的应力增量矩阵；ν为材料的弹性应变泊松比。

实体单元徐变应变增量矩阵Δε_c,n如式(3.16)所示：

其中，

式中，为与应力分量对应的向量，C为徐变泊松比矩阵，ν’为材料的徐变泊松比；其余符号意义同前。

混凝土单元收缩应变增量矩阵Δε_s,n如式(3.19)所示：

Δε_s,n＝ε_s,n-ε_s,n-1 (3.19)

其中，

ε_s,n＝[ε_s,n ε_s,n ε_s,n 0 0 0]^T (3.20)

式中，ε_s,n为单元在第n个时段的收缩应变矩阵；ε_s,n-1为单元在第n-1个时段的收缩应变矩阵；其余符号意义同前。

由式(3.13)和(3.14)，应力增量矩阵可按式(3.21)计算。

Δσ_n＝D_nΔε_e,n＝D_n(Δε_n-Δε_c,n-Δε_s,n) (3.21)

其中，

D_n＝E(t_n-1)A^-1 (3.22)

式中，D_n为单元材料本构关系矩阵；其余符号意义同前。

把单元应变增量与节点位移增量关系式和式(3.16)代入式(3.21)，得到单元应力增量矩阵，如式(3.23)所示：

其中，

式中，B为单元几何矩阵；为单元在第n个时段的弹性变形增量矩阵；α_n为单元在第n个时段的徐变应变增量矩阵。

有限元法的平衡方程组可表示为式(3.25)：

式中，为单元外荷载增量矩阵。

将式(3.23)代入式(3.25)，得到徐变与收缩效应分析的基本方程如式(3.26)所示：

其中，

式中，k_n为单元刚度矩阵；为徐变变形产生的当量荷载增量矩阵；为收缩变形产生的当量荷载增量矩阵。

通过坐标转换，按节点序号将单元刚度矩阵、荷载矩阵集成结构刚度矩阵和荷载矩阵，形成结构平衡方程并可求解节点位移，然后将节点位移坐标转换后再由式(3.23)计算单元应力增量矩阵。

步骤e：运用Fortran语言在ANSYS软件的USERMAT子程序中开发混凝土时变本构，实现空间应力状态下混凝土结构在复杂应力历程中的时变效应有限元分析。

结合图2，在USERMAT子程序中开发混凝土徐变与收缩本构模型的主要步骤如下：

e1)输入徐变系数、收缩应变计算参数(环境相对湿度、构件理论厚度、混凝土强度等级等等)，弹性模量随龄期变化函数，加载龄期t₀及收缩龄期t_s等数据；

e2)主程序从时间t_n-1＝0开始计算，令t_n-1小于等于0.01天时，只计算弹性应变，得到初始应力增量并存储在状态变量数组中，为后面求解徐变应变增量提供应力计算参数；

e3)当t_n-1大于0.01天时，根据式(3.11)建立徐变递推计算式，混凝土加载龄期在计算式中考虑，计算当前时间步列阵，并将列阵存储于状态变量数组中，由式(3.16)计算当前时间步的徐变应变增量Δε_c,n。

e4)根据规范收缩模型计算当前时间的收缩应变ε_s,n，由式(3.19)、(3.20)计算当前时间步收缩应变增量列阵Δε_s,n，并将列阵ε_s,n存储于状态变量数组中。

e5)根据式(3.15)、(3.22)由当前时间步的弹性模量和泊松比计算混凝土的本构矩阵D_n，在进行结构徐变收缩效应分析时，混凝土按线弹性材料处理，对于线弹性材料，一致切线算子矩阵D_ep与D_n相同。

e6)根据式(3.21)计算混凝土应力增量列阵Δσ_n，并存储于状态变量数组，在时间步结束时更新单元应力，完成结构的徐变、收缩效应分析。

以混凝土圆柱体试件的徐变、收缩试验为实施例，按照上述方法对混凝土单轴应力作用下的时变效应进行分析。

1、应用ANSYS软件建立实体有限元模型

试件直径为117.5mm，高305mm，混凝土采用20节点实体单元SOLID186，有限元模型如图3所示。混凝土按线弹性材料处理，其强度特征值为f_ck＝46.3MPa，加载龄期混凝土弹性模量E_c约为35GPa，且考虑混凝土弹性模量的时间效应，泊松比取0.167。试件底面为固结约束，试件在混凝土浇筑后第8天开始加载，第120天完全卸载，在试件顶面采用均布应力方式加载。

2、将徐变系数拟合成多项指数函数相加的形式

采用水灰比为0.375的快速硬化硅酸盐水泥振捣压实制成，试件在温度为17℃，相对湿度为93％的室内潮湿固化。采用fib MC 2010徐变模型，徐变系数随时间发展函数为：

将其拟合成形式，取m＝3，曲线拟合参数如表1所示。

表1曲线拟合参数

徐变模型拟合曲线与原始曲线比较如图4所示，以1天为间隔，取161组数据计算得出相关系数为0.9987，表明该曲线拟合精度较高，能满足工程计算精度的要求。

3、运用Fortran语言在USERMAT子程序中开发徐变、收缩计算模块，将编写好的USERMAT.F文件和ANSYS软件安装目录下的custom文件夹中的user文件夹一起拷贝到工作目录下，运行anscust.bat，自动编译USERMAT.F，编译成功后生成ANSYS.exe文件，该文件就是添加了用户自定义子程序的可执行文件，在ANSYS Mechanical APDL Product Launcher窗口的Customization/Preferences中连接ANSYS.exe文件。在前处理部分通过APDL命令TB，USER及TB，TBDATA调用USERMAT子程序。模型求解部分，时间步长设为1d，对试件加载全过程的应变进行精细化分析。

Claims (6)

1.一种混凝土结构时变效应的实体有限元计算方法，其特征在于，包括以下步骤：

步骤a：采用逐步法将混凝土的应力历程离散化，推导出第n个时段内的单轴徐变应变增量；

步骤b：采取线性最小二乘法初步拟合，非线性规划求极值法最终拟合的策略，将设计规范推荐的徐变系数表达式拟合成指数级数形式；

步骤c：将徐变系数拟合表达式代入单轴徐变应变增量计算式，建立混凝土单轴徐变应变增量递推计算式；

步骤d：考虑多轴应力工况的徐变泊松影响，在实体单元材料本构中引入混凝土时变效应递推计算式，建立混凝土实体单元时变本构方程；

步骤e：运用Fortran语言在ANSYS软件的USERMAT子程序中开发混凝土时变本构，实现空间应力状态下混凝土结构在复杂应力历程中的时变效应有限元分析。

2.根据权利要求1所述的一种混凝土结构时变效应的实体有限元计算方法，其特征在于，所述的步骤a中，当结构截面应力连续变化时，将应力历程划分为若干个相互衔接的小过程，在每一个小过程内假定应力为常量，并以该应力状态发生徐变，则第n个时段内的单轴徐变应变增量如式(1.1)所示：

式中，Δε_c(t_n,t_n-1)为混凝土在t_n-1至t_n时段内的单轴徐变应变增量；Δσ_i为t_i时刻的应力增量；E(t_i)为t_i时刻混凝土的弹性模量；为加载龄期t_i的混凝土在t_n时刻的徐变系数；为加载龄期t_i的混凝土在t_n-1时刻的徐变系数。

3.根据权利要求1所述的一种混凝土结构时变效应的实体有限元计算方法，其特征在于，所述的步骤b中，混凝土的徐变系数表达式如式(1.2)所示：

式中，为加载龄期t_i的混凝土在t时刻的徐变系数；为混凝土在t_i时刻的名义徐变系数；β_H为与环境相对湿度及构件理论厚度相关的系数；

采取线性最小二乘法初步拟合，等式约束条件下非线性规划求极值最终拟合的策略，运用Matlab软件将式(1.2)数值拟合为如式(1.3)所示的指数级数形式：

式中，C_j、q_j(j∈[1,m])为拟合常数；m为展开的级数项数，一般m≥3有足够的拟合精度；其余符号意义同前。对于欧洲fib Mode Code 2010规范、美国AASHTO LRFD桥梁规范等亦能拟合为式(1.3)形式。

4.根据权利要求1所述的一种混凝土结构时变效应的实体有限元计算方法，其特征在于，所述的步骤c中徐变应变增量递推计算式的详细推导过程如下：

将展开后的徐变系数级数表达式(1.3)代入式(1.1)，利用指数函数特点进一步合并，如式(1.4)所示：

式中，Δt_n为第n个时段的时间增量，即Δt_n＝t_n-t_n-1；其余符号意义同前；

对式(1.4)的累加项进行重组，并令：

式中，为递推变量，其余符号意义同前；

于是有

同样，第n+1个时段的徐变应变增量也可用式(1.7)表示：

其中，

比较式(1.5)和式(1.8)，可得到如式(1.9)所示的递推关系：

混凝土收缩应变与应力无关，当构件尺寸、材料特性、环境状况和收缩开始龄期已知时，收缩应变随时间发展曲线唯一确定；与徐变相似，将收缩应变曲线离散化，则第n个时段的收缩应变增量Δε_s(t_n,t_n-1)如式(1.10)所示：

Δε_s(t_n,t_n-1)＝ε_s(t_n,t_s)-ε_s(t_n-1,t_s) (1.10)

式中，t_s为混凝土收缩开始龄期；ε_s(t_n,t_s)为混凝土从t_s至t_n时段的收缩应变；ε_s(t_n-1,t_s)为混凝土从t_s至t_n-1时段的收缩应变；

利用上述递推关系，第n个时段的徐变与收缩应变增量Δε(t_n,t_n-1)计算如式(1.11)所示：

数值分析时，采用该递推计算式不必存储应力历程，大大节约了计算资源和时间消耗。

5.根据权利要求1所述的一种混凝土结构时变效应的实体有限元计算方法，其特征在于，所述的步骤d中在实体单元建立时变本构方程的过程为：

设第n个时段内实体单元的应力增量矩阵Δσ_n与应变增量矩阵Δε_n如式(1.12)所示：

式(1.12)中单元的应变增量矩阵Δε_n可表示为如式(1.13)所示的三部分：

Δε_n＝Δε_e,n+Δε_c,n+Δε_s,n (1.13)

式中，Δε_e,n为弹性应变增量矩阵；Δε_c,n为徐变应变增量矩阵；Δε_s,n为收缩应变增量矩阵。

实体单元弹性应变增量矩阵如式(1.14)和(1.15)所示：

其中，

式中，Δε_e,n为单元在第n个时段的弹性应变增量矩阵；E(t_n-1)为材料在t_n-1时刻的弹性模量；A为弹性泊松比矩阵；Δσ_n为单元在第n个时段的应力增量矩阵；ν为材料的弹性应变泊松比。

实体单元徐变应变增量矩阵Δε_c,n如式(1.16)所示：

其中，

式中，为与应力分量对应的向量，C为徐变泊松比矩阵，ν’为材料的徐变泊松比；其余符号意义同前；

混凝土单元收缩应变增量矩阵Δε_s,n如式(1.19)所示：

Δε_s,n＝ε_s,n-ε_s,n-1 (1.19)

其中，

ε_s,n＝[ε_s,n ε_s,n ε_s,n 0 0 0]^T (1.20)

式中，ε_s,n为单元在第n个时段的收缩应变矩阵；ε_s,n-1为单元在第n-1个时段的收缩应变矩阵；其余符号意义同前；

由式(1.13)和(1.14)，应力增量矩阵可按式(1.21)计算：

Δσ_n＝D_nΔε_e,n＝D_n(Δε_n-Δε_c,n-Δε_s,n) (1.21)

其中，

D_n＝E(t_n-1)A^-1 (1.22)

式中，D_n为单元材料本构关系矩阵；其余符号意义同前；

把单元应变增量与节点位移增量关系式和式(1.16)代入式(1.21)，得到单元应力增量矩阵，如式(1.23)所示：

其中，

式中，B为单元几何矩阵；为单元在第n个时段的弹性变形增量矩阵；α_n为单元在第n个时段的徐变应变增量矩阵。

有限元法的平衡方程组可表示为式(1.25)：

式中，为单元外荷载增量矩阵；

将式(1.23)代入式(1.25)，得到徐变与收缩效应分析的基本方程如式(1.26)所示：

其中，

式中，k_n为单元刚度矩阵；为徐变变形产生的当量荷载增量矩阵；为收缩变形产生的当量荷载增量矩阵；

通过坐标转换，按节点序号将单元刚度矩阵、荷载矩阵集成结构刚度矩阵和荷载矩阵，形成结构平衡方程并求解节点位移，然后将节点位移坐标转换后再由式(1.23)计算单元应力增量矩阵。

6.根据权利要求1所述的一种混凝土结构时变效应的实体有限元计算方法，其特征在于，所述的步骤e中，在USERMAT子程序中开发混凝土徐变与收缩本构模型的主要步骤如下：

步骤e1：输入徐变系数、收缩应变计算参数(环境相对湿度、构件理论厚度、混凝土强度等级等等)，弹性模量随龄期变化函数，加载龄期t₀及收缩龄期t_s等数据；

步骤e2：主程序从时间t_n-1＝0开始计算，令t_n-1小于等于0.01天时，只计算弹性应变，得到初始应力增量并存储在状态变量数组中，为后面求解徐变应变增量提供应力计算参数；

步骤e3：当t_n-1大于0.01天时，根据式(1.11)建立徐变递推计算式，混凝土加载龄期在计算式中考虑，计算当前时间步列阵，并将列阵存储于状态变量数组中，由式(1.16)计算当前时间步的徐变应变增量Δε_c,n；

步骤e4：根据规范收缩模型计算当前时间的收缩应变ε_s,n，由式(1.19)、(1.20)计算当前时间步收缩应变增量列阵Δε_s,n，并将列阵ε_s,n存储于状态变量数组中；

步骤e5：根据式(1.15)、(1.22)由当前时间步的弹性模量和泊松比计算混凝土的本构矩阵D_n，在进行结构徐变收缩效应分析时，混凝土按线弹性材料处理，对于线弹性材料，一致切线算子矩阵D_ep与D_n相同；

步骤e6：根据式(1.21)计算混凝土应力增量列阵Δσ_n，并存储于状态变量数组，在时间步结束时更新单元应力，完成结构的徐变、收缩效应分析。