基于ANSYS二次开发的早龄期混凝土徐变应力计算

郭浩洋，戚　蓝，李少明，刘　勇，赵芳兴

(1.中国电建集团 西北勘测设计研究院有限公司，西安　710065；2.天津大学 水利工程仿真与安全国家重点实验室，天津　300072；3.杭州市水利规划研究与科技推广中心，杭州　310016)



基于ANSYS二次开发的早龄期混凝土徐变应力计算

郭浩洋1，戚蓝2，李少明3，刘勇2，赵芳兴2

(1.中国电建集团 西北勘测设计研究院有限公司，西安710065；2.天津大学 水利工程仿真与安全国家重点实验室，天津300072；3.杭州市水利规划研究与科技推广中心，杭州310016)

徐变是混凝土材料所具有的时变特性，对早龄期混凝土的应力影响很大，不可忽略，然而工程界主流的通用有限元程序不支持混凝土徐变的计算，不利于工程应用。为此，笔者通过对ANSYS用户子程序usermat进行二次开发，编写了自定义的材料模型，利用ANSYS程序实现了混凝土弹性徐变方程的隐式解法。在此基础上，通过算例验证了该二次开发子程序的正确性，并将其应用于实际工程，通过对比，采用弹性徐变方程隐式解法求得的计算结果与实测值的变化规律相同，数据也比较吻合，误差范围一般在0.3MPa以内，证明该子程序具有较好的精度，可以推广应用于工程早龄期混凝土温度徐变应力的预测和计算。

温度应力场；早龄混凝土；徐变；ANSYS二次开发；隐式解法

1　研究背景

徐变指的是混凝土结构在持续荷载作用下，变形随时间不断增加的现象。混凝土的徐变变形约为其弹性变形的1到2.5倍[1]，且在加载早期发展很快[2]。因此，对于主要荷载为温度荷载的早龄期混凝土，徐变对结构应力的影响不能忽略。

实际工程大多结构复杂，边界条件多变，用传统的分析方法求解早龄期混凝土的应力，很难得到准确结果，目前主要采用有限单元法来计算。相对于编制有限元计算程序，利用通用有限元软件来实现混凝土温度徐变应力场的计算大大降低了问题的难度，增加了灵活性，便于在工程上推广应用。因此，很多学者都在此方向作出了努力，取得了一定的成果，但也有一定的局限性。如高越权等[3]、张治成[4]根据调整龄期的有效模量法，在每一时段赋予材料对应等效弹模来模拟混凝土的徐变，但该方法不能考虑应力历史，且认为徐变完全可复，计算结果不够精确；王维红等[5]考虑采用ANSYS自带的金属蠕变本构关系曲线代替混凝土徐变系数曲线的方法分析混凝土的徐变应力，但该方法假定混凝土弹模为常量，不能模拟弹模增长迅速的早龄期混凝土；何继访等[6]通过开发ANSYS蠕变子程序usercreep，编写蠕变本构关系来模拟混凝土的徐变，但受限于蠕变子程序的结构，不能模拟混凝土的卸载。采用弹性徐变老化理论模型模拟混凝土的徐变，不仅能够反映应力历史的影响，还能较好地反映早期混凝土在卸载作用下徐变的部分可恢复性质。李骁春等[7]通过对ANSYS进行二次开发，将弹性徐变老化理论引入ANSYS，实现了对混凝土徐变的计算，但是实现的方法较为繁琐，联合了多个子程序且在每一个迭代步都需要重新启动求解，计算效率较低。本文通过开发Usermat子程序，编写自定义材料模型，将弹性徐变老化理论引入ANSYS，实现方法较为简单且提高了计算效率，经实例验证，计算结果精度满足工程要求。

2　混凝土弹性徐变计算原理

由于混凝土的徐变不仅与当前的应力有关，而且与历史的应力有关，故一般计算过程中必须记录应力历史，需要较大的存储容量。文献[8]中给出了混凝土弹性徐变方程的隐式解法，该算法不必记录应力历史，提高了计算精度和计算效率。

设混凝土徐变度为：

(1)

(2)

(3)

(4)

上述式(3)和式(4)中ωs，n为迭代算子。这是一组递推公式，在用有限元方法分析混凝土结构徐变时，不必记录应力历史，只需存储ωs，n，从而可以大大节省存储量。

(5)

式中E0，a，b均为材料常数。

由于弹性模量和徐变度都随时间而变化，故用增量法进行分析。把时间τ划分成一系列时段，则在Δτn时段内产生的应变增量为

(6)

混凝土徐变计算的隐式解法假定Δτn时段内应力呈线性变化，即∂σ/∂τ=常量，由积分中值定理可得弹性应变增量为

(7)

徐变应变增量式(2)已经给出，由式(2)、式(6)和式(7)可得应力增量与应变增量的关系，即

(8)

(9)

(10)

在有限单元法中，平衡方程为

(11)

式中：左边是应力增量对体积的积分，其中[B]为几何矩阵；右边是外荷载增量。将式(8)代入式(11)，得到基本方程

(12)

3　ANSYS徐变材料模型的开发

Usermat用户子程序是用户开发材料模型最有力的工具。该子程序通过定义材料的本构行为(应力-应变关系)来实现材料模型的开发，具体分为两个方面[10-11]：一是给出一致切线算子矩阵，即ANSYS中的雅克比矩阵(JacobianMatrix，D=∂Δσ/∂Δε)；二是由给定的应变增量(Δε，由ANSYS主程序传入Usermat)计算应力增量，从而得到新的应力，称作应力更新过程。因此，要实现徐变材料模型的开发，就需要在Usermat子程序中完成弹性徐变材料应力-应变关系的定义。

本文通过编写Usermat子程序来实现混凝土弹性徐变方程的隐式解法。其程序的流程见图1。

图1　Usermat子程序流程Fig.1　Flow chart of the Usermat subroutine

编译连接修改好的Usermat子程序，生成用户定制版本的ANSYS，运行程序，在前处理部分调用自定义材料模型，即可进行混凝土徐变的计算。

4　算例验证

本节采用文献[8]中的一则算例来验证徐变子程序的正确性。如图2所示，基岩上单层浇筑的混凝土块，长度L=25m，厚度h=3m，其表面与空气接触，计算所采用的参数详见文献[8]。计算混凝土块由于水化热作用及天然冷却而产生温度场及温度徐变应力场。

尺寸单位：m图2　基岩上的单层浇筑混凝土块Fig.2　Single layer of concrete block on the basement

选取对称结构，建立三维有限元模型计算混凝土块的温度场和徐变应力场。图3(a)为中央断面A-A上不同龄期的温度分布，前期由于水化热作用，温度不断升高，在τ=10d时，内部温度达到最高值13.21 °C(文献[8]中为13.20 °C)；以后逐渐下降，到τ=100d时，最高温度降为3.32 °C(文献[8]中为3.1 °C)。中央断面A-A上随龄期变化的水平向温度徐变应力σx分布规律如图3(b)。混凝土块浇筑初期由于温度升高，混凝土膨胀，在τ=9d之前，全断面受压，之后表面开始受拉，受拉范围逐步向内部扩展，到τ=50d时已经全断面受拉，τ=100d时拉应力最大，达到1.018 9MPa(文献[8]中为1.02MPa)。

文献[8]中算例采用二维平面有限元法计算，本文建立了三维模型且在网格划分、计算步长取值等方面与原算例均略有差别，但计算结果误差很小。本文方法计算出的混凝土块的温度、应力变化规律与文献[8]相同，验证了本子程序的正确性。

5　实际工程应用

天津某大型泵站工程由于工程进度安排，需要在冬季施工。作为泵站调节池重要组成部分的扶壁式挡土墙和导墙等大体积混凝土墙体在正常工况下需要承受巨大的水压力，因此结构的裂缝控制十分重要。为了研究寒冷环境条件下新浇扶壁式挡墙内部的温度变化及温度应力变化，在某一标段扶壁式挡土墙浇筑前，在墙体内部埋设了电子测温仪和钢筋计。该标段扶壁式挡土墙底板部分为先期浇筑，挡土墙的具体形式及主要尺寸见图4。

(a)中央断面A—A温度分布

(b)中央断面A—A温度徐变应力分布图3　中央断面A—A温度和温度徐变应力分布Fig.3　Temperature distribution and creep stresses onthe central section A—A

图4　挡土墙示意图Fig.4　Schematic of the retaining wall

在距离底板表面高度为0.3,2.0,5.0m的壁板中心位置分别埋设A,B,C3组电子测温仪和钢筋计，以记录壁板混凝土的温度变化和应力变化。挡土墙混凝土浇筑完毕后及时开始监测工作。

建立有限元模型，根据已有参数及监测到的环境温度，计算扶壁式挡土墙的温度场及徐变应力场。由于本工程未做徐变实验，徐变度采用文献[8]中推荐的用于初步设计的公式，即

式中：C1=0.23/E0；C2=0.52/E0；E0=1.45E(28)；E(28)为混凝土28d的弹性模量。

为验证有限元计算温度场的正确性，取挡土墙壁板内距离底板表面高度为2.0m处B监测点记录的温度数据与有限元模型相应位置节点计算的温度结果进行对比。对比结果如图5所示，可以看出有限元模拟结果与实测结果十分接近。说明有限元温度场计算结果精度满足要求。

图5　B监测点实测温度与计算温度对比Fig.5　Comparison between measured temperature andcalculated temperature at point B

图6为有限元计算所得的沿高度方向挡土墙壁板(附带底板)中面中线温度随时间的变化曲线。由图6可见：挡土墙浇筑后，内部温度迅速上升，并于τ=1.5d时达到最高温度24.8 °C，之后由于天然冷却而开始降温；到τ=7d时壁板中线大部分温度只有约0.5 °C，接近于环境温度。故主要的温度变化可以认为发生在前7d。壁板中线大部分范围温度分布比较均匀，在底板和顶部变化较为剧烈，据此可以推测，在这些温度变化剧烈的区域会产生较大的温度应力。

图6　挡土墙壁板中面中线温度随时间变化曲线Fig.6　Temperature vs. time along the central line of themiddle plane of the retaining wall

图7是壁板内3组不同高度钢筋计测量结果换算得出的混凝土水平向应力实测值与相应位置处有限元计算结果随时间的对比图。由图7可以看出，采用弹性徐变方程隐式解法求得的计算结果与实测值的变化规律相同，数据也比较吻合，误差范围一般在0.3MPa以内。据此可以认为，有限元温度徐变应力场计算结果满足工程精度要求。

(a)A监测点

(b)B监测点

(c)C监测点图7　监测点实测水平应力值与计算值对比Fig.7　Comparison between measured horizontal stressesand calculated values of measurement points

沿高度方向，壁板内部中面中线水平向温度徐变应力的有限元计算结果分布如图8所示。可以看出：温度徐变应力的主要影响区域在0.8～4m的范围内，超出此范围的部分，混凝土温度徐变应力相对较小。而该影响区正是混凝土壁板与底板连接的区域。由于底板先期浇筑，混凝土经过天然冷却，故其内部温度较低且相对稳定。壁板混凝土在浇筑初期，由于水泥水化热作用体积膨胀，此时就会受到底板的约束而产生压应力。由图8可知：在τ=1.5d时影响区断面底部仍然受压。随着时间推移，壁板混凝土温度迅速下降，同样由于底板的约束，壁板中面开始产生拉应力；在τ=2d时，影响区断面已经由受压变为受拉；之后拉应力不断增大，在τ=6d时达到最大值1.74MPa。

图8　壁板中面中线水平向温度徐变应力分布Fig.8　Horizontal distribution of creep stresses along thecentral line of the middle plane of retaining wall

在τ=6d时，相同养护条件下混凝土试块的抗压强度为15.3MPa，而壁板中线水平向拉应力达到了1.74MPa，极有可能在壁板上引起竖向的裂缝。后经现场踏勘，挡土墙壁板内外两侧均出现微裂缝，裂缝沿高度方向向两侧斜向对称发展，且中间有一条竖向裂缝作为对称线。

之所以会出现温度裂缝，一方面是因为挡土墙墙体分两期浇筑，两期混凝土之间温差较大；另一方面，挡土墙在冬季施工，且模板为钢模板，基本上起不到保温作用，造成墙内混凝土温度变化过快，温度变化幅度较大。

6　结　论

(1)本文通过重新编写ANSYS用户子程序Usermat，将混凝土弹性徐变方程的隐式解法引入ANSYS，实现方法较为简单且有较高的计算效率，并由算例验证了子程序的正确性， 实现了对混凝土徐变的计算。

(2)相对于编制专门的有限元计算程序，本文的工作依托ANSYS软件强大的前后处理和计算能力，极大减小了程序开发的工程量，且具有较好的通用性。

(3)由实际工程应用经验可知，该子程序具有较好的精度，可以推广应用于工程早龄期混凝土温度徐变应力的预测和计算。

(4)由本文实际工程部分的计算可知，寒冷环境下分期浇筑的薄壁混凝土结构容易产生温度裂缝，工期应尽量避免安排在冬季。如果需要在寒冷季节施工，为减少温度裂缝的发生，可采取以下措施：①采用低发热量的混凝土；②对先期浇筑的混凝土采用表面吹热风或内部埋水管通热水加温等措施，减小先后浇筑的两期混凝土之间的相对温差；③做好新浇薄壁混凝土的保温工作，以减小降温的速度和幅度。

[1]朱伯芳.混凝土结构徐变应力分析的隐式解法[J].水利学报，1983(5)：40-46.

[2]熊维.不同强度早龄期混凝土徐变及徐变对长期荷载作用下预应力构件的影响[D].天津：天津大学，2011.

[3]高岳权，黄浩良，王佶，等.基于等效模量法与ANSYS计算混凝土徐变[J].武汉理工大学学报，2012，32(3)：13-15.

[4]张治成.大跨度钢管混凝土拱桥的徐变分析 [J].工程力学，2007，24(5)：151-160.

[5]王维红，王运涛，陈宇峰,等.白果渡嘉陵江大桥徐变仿真分析[J].重庆交通大学学报(自然科学版)，2008，27(1)：9-12.

[6]何继访，杨文兵，苏砺锋,等.基于ANSYS平台的混凝土徐变计算[J].湖南交通科技，2006，32(3)：108-110，159.

[7]李骁春，吴胜兴.基于ANSYS的混凝土早期徐变应力仿真分析[J].系统仿真学报，2008，20(15)：3944-3947.

[8]朱伯芳.大体积混凝土温度应力与温度控制[M].北京：中国水利水电出版社，2012

[9]朱伯芳，吴龙珅，李玥,等.加热下部混凝土以防止混凝土结构裂缝的探索[J].水利水电技术，2009，40(2)：38-41.

[10]ANSYS.GuidetoANSYSUserProgrammableFeatures[K].U.S.:ANSYSInc, 2005.

[11]ANSYS.ANSYSProgrammer’sManual[K].U.S.:ANSYSInc, 2002.

(编辑：赵卫兵)

Calculation of Creep Stress of Early-age Concrete Based on FurtherDevelopment of ANSYS

GUO Hao-yang1, QI Lan2，LI Shao-ming3，LIU Yong2，ZHAO Fang-xing2

(1.PowerChinaNorthwestEngineeringCorporationLimited,Xi’an710065,China;2.StateKeyLaboratoryofHydraulicEngineeringSimulationandSafety,TianjinUniversity,Tianjin300072,China;3.WaterConservancyResearchandTechnologyPromotionCentreofHangzhou,Hangzhou310016,China)

Asatime-dependentpropertyofconcretematerial,creepimposessignificantinfluenceonthestressof

thermalstressfield;early-ageconcrete;creep;furtherdevelopmentofANSYS;implicitsolution

2015-07-20;

2015-08-26

郭浩洋(1989-)，男，河南洛阳人，助理工程师，硕士，研究方向为混凝土早期裂缝控制研究，(电话)13652117580(电子信箱)509341595@qq.com。

戚蓝(1955-)，女，天津人，教授，博士，研究方向为水利工程安全监测等，(电话)13752530685(电子信箱)qilan0000@163.com。

10.11988/ckyyb.20150622

2016,33(09):138-142，154

TV331

A

1001-5485(2016)09-0138-05

early-ageconcrete.However,itisdifficultinengineeringtocalculatethecreepstressbyusingcommonavailablefiniteelementprograms.Inviewofthis,throughfurtherdevelopmentoftheANSYSusersubroutineusermat,auser-definedmaterialmodelwasdevelopedtointroducetheimplicitsolutiontotheconcretecreepequation.Withthelargesterrorlowerthan0.3MPa,theusersubroutineusermatwasverifiedthroughcalculationexamplesandwasappliedtoengineering.Theresultsoftheimplicitconcretecreepequationareingoodagreementwithmeasuredvalues,andtheirvariationtrendsarethesame.Theusersubroutineusermatcanbeappliedtoengineeringtopredictandcalculatethethermalcreepstressofearly-ageconcrete.