一种新型多参数影响下早期混凝土性态变化研究方法与流程
本发明涉及建筑工程技术领域，尤其涉及一种新型多参数影响下早期混凝土性态变化研究方法。
背景技术：
混凝土是当今使用最广、用量最大的建筑材料。由于混凝土具有工程所需要的强度和耐久性，且原料易得，造价较低，特别是能耗较低，因此用途极为广泛。相比混凝土材料的抗压强度，其抗拉强度较低，是一种容易开裂的材料，混凝土内部温度和应力、应变机理引起了众多行业的关注。混凝土浇筑完成初期，混凝土内部温度急剧上升，混凝土的弹性模量随龄期变化，易产生热胀变形。当混凝土产生的热量小于释放的热量时，混凝土温度开始降低，在此阶段会产生收缩变形；当应力超过混凝土的抗拉强度时，混凝土便会产生裂缝。早期混凝土性态的发展对施工期混凝土裂缝的产生有较大影响。
混凝土的早期热、力学特征和裂缝控制仍未能得到很好的解决，主要是因为早期混凝土材料的相关性能较为复杂，测试上也存在难度。目前已有部分学者对早期混凝土性态变化进行了一系列研究，但多数仅考虑常规因素(如水化热、应变、拉伸徐变等)影响下结构性态的变化规律。现阶段少有学者同时考虑变化的热膨胀系数和混凝土徐变、收缩变形引起的干缩应力以及塑性理论对早期混凝土结构性态变化的影响，有待进一步深入研究，为工程施工及科学研究提供理论参考。
并且，目前现有技术中已有的实时模拟方法往往在控制参数变化方面考虑不足，导致计算结果不准确。
技术实现要素：
鉴于上述问题，本发明提供一种新型多参数影响下早期混凝土性态变化研究方法，其方法实现方便、可操作性强且使用效果好、实验条件限制小，可有效节约研究成本，且计算结果精度高。
本发明提供一种新型多参数影响下早期混凝土性态变化研究方法，包括：
(1)依据工程所需选定合适的试验材料；
(2)依据相关混凝土试验规范，确定混凝土配合比；
(3)设计试验方案，并获取相关的混凝土物理参数，所述混凝土物理参数包括：导热系数、密度、比热和线性膨胀系数；
(4)用温度传感器与振弦传感器分别测量混凝土试件内部的温度与应变；
(5)结合二次开发改进后遗传算法对所述试验方案得出的热力学参数进行反演；所述热力学参数包括：导温系数、导热系数、比热和密度；
(6)基于反演参数及abaqus的二次开发平台，开发了适用于早期混凝土温度场、湿度场计算的子程序，根据所述工程建立有限元模型进行早期混凝土温、湿度场数值模拟计算；
(7)依托abaqus二次开发平台，通过调用所述温度场、湿度场和开发的应力场子程序，实现了早期黏弹塑性混凝土的热湿力学耦合模型和数值求解过程；再结合所述线性膨胀系数，研究早期混凝土的性态变化规律，实现早期混凝土热湿力的时变性耦合分析。
进一步地，早期黏弹塑性混凝土的热湿力学耦合模型和数值求解过程，包括：
采用模拟方法计算结构温度场；
根据所述温度场子程序来模拟结构相对湿度场分布；
将所述温度场和所述湿度场施加于结构进行应力场的计算。
进一步地，所述湿度场扩散控制方程为：
式中，h为混凝土结构的相对湿度，用百分数表示；t为时间；dx(h)、dy(h)、dz(h)分别为x、y、z三个方向的与相对湿度相关的湿度扩散系数；
进一步地，所述湿度场的初始条件，包括：
当湿度场模拟从浇筑开始后就开始模拟时，所述湿度场的初始条件可以假定为100％，即：
h＝h(x,y,z,t)＝h(x,y,z,t0)＝1.0(3)
当混凝土湿度场模拟从浇筑中间开始时，通过试验测得计算时刻的初始湿度值，并定为初始湿度值。
进一步地，所述湿度场的边界条件，包括：
a)第一类边界条件
混凝土结构表面相对湿度h为时间的已知函数，即有
h＝h(x,y,z,t)(4)
b)第二类边界条件
混凝土表面为绝湿边界，即
式中：n为表面法线方向；
c)第三类边界条件
当混凝土与空气接触时，用对流边界条件来表示该种边界条件，即混凝土的表面湿度交换量和表面湿度hs与大气相对湿度he之差成正比，其数学表达式为：
式中，f为表面湿度扩散系数，单位为m2/h；
d)第四类边界条件
当两固体接触时，若接触良好，则在两物体接触面上的湿度和湿度量应该都是连续的，边界条件为：
h1＝h2(7)
进一步地，所述应力场计算基于abaqus二次开发平台，在求解非线性问题时采用newton-raphson算法。
进一步地，所述newton-raphson算法为增量步算法；
为了求得在t时刻的基本徐变变形，已知徐变函数j(t,τ)，则应变张量ε(t,τ)表达为：
式中，[d]为弹性矩阵；σ(τ)为龄期τ处应力张量；ε0(t)为与应力无关的应变张量；dσ(τ)为σ(τ)的导数；
在时间增量δt＝ti+1-ti内，相应的应力张量和应变张量增量分别为：
[δσ]＝[σ(ti+1)]-[σ(ti)](31)
和
[δε]＝[ε(ti+1)]-[ε(ti)](32)
利用公式(30)，则应变增量可以表示为：
进一步地，0所述应力场子程序与所述abaqus的ti主程序之间的协同计算过程如下：
第一步：在平衡时刻tn，主程序向子程序提供总时间，时间增量δt，总应变εtotal(tn)和应变增量δεtotal(tn)，以及根据平衡时刻应力σ(tn)，分别传入子程序用来计算新的柯西应力张量σ(tn+δt)；
第二步：根据上一步更新的状态变量，计算本增量步的粘性应力
第三步：abaqus主程序更新应变张量εtotal(tn+δt)＝εtotal(tn)+δεtotal(tn)，产生新的总应变εtotal(tn+1)；
第四步：abaqus根据新的应力、应变与外荷载进行平衡迭代计算，默认最大迭代次数为16次；在最大迭代次数内达到收敛之后，会进行下一个增量步的计算；反之，在运行情况下主程序会自动减小增量步大小，再次进行计算分析，默认最大增量步次数为5次，若还是不收敛计算将会终止并报错。
进一步地，所述步骤(4)包括：
①混凝土浇筑前，在浇筑模具内表面接缝处涂抹约2mm厚硅胶；
②将聚苯泡沫板贴在模具内表面，聚苯板接缝处密封；
③振弦传感器和温度传感器固定在模具内部，其中一振弦传感器记为vws-t；另一个记为vws-b；温度传感器从上往下依次记为ts-1、ts-2、ts-3；
④称量所需材料的质量，采用搅拌机搅拌，搅拌均匀后将混凝土浇筑在模具内部，立即使用振捣棒振捣使试件振捣密实，以排除气泡；
⑤开始采集试件的温度与应变数据。
进一步地，所述二次开发改进后的遗传算法包括：
通过交叉概率和变异概率随适应度的大小及分散程度进行自适应调整；
通过式(50)、(51)实现，如下所示:
式中：pc表示交叉概率；pm表示最大变异概率；fmax为群体中最大的适应度值；favg为每代群体的平均适应度值；f′为要交叉的两个个体中较大的适应度值；f为要变异的个体的适应度值；k1、k2、k3、k4为自适应控制参数，取为k1＝k2＝0.9，k3＝k4＝0.1。
本发明实施例提供的上述技术方案的有益效果至少包括：
本发明实施例提供的一种新型多参数影响下早期混凝土性态变化研究方法，能够实时获得早期混凝土温度、湿度，该方法主要通过数值仿真、室内试验与遗传算法相结合的方式来对数值仿真模型相关参数进行实时更新，在保证精确模拟模型温度场、湿度场的情况下，进行混凝土结构物较为精确的应力场控制与模拟。
本发明的其它特征和优点将在随后的说明书中阐述，并且，部分地从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点可通过在所写的说明书、权利要求书、以及附图中所特别指出的结构来实现和获得。
下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。
附图说明
附图用来提供对本发明的进一步理解，并且构成说明书的一部分，与本发明的实施例一起用于解释本发明，并不构成对本发明的限制。在附图中：
图1为本发明提供的遗传算法的流程图；
图2为本发明提供的温度场、湿度场和应力场三场顺序耦合示意图；
图3a为本发明提供的加载龄期的影响的曲线图；
图3b为本发明提供的持载时间的影响的曲线图；
图4为本发明提供的增量法示意图；
图5为本发明提供的umat子程序计算流程图；
图6为本发明提供的umatht子程序流程图；
图7为本发明提供的混凝土放热与升温曲线示意图；
图8为本发明提供的ta测点反演计算值与测量值的曲线图；
图9为本发明提供的tb测点反演计算值与测量值的曲线图；
图10为本发明提供的tc测点反演计算值与测量值的曲线图；
图11a为本发明提供的混凝土与传感器位置的示意图；
图11b为本发明提供的图11a沿a-a、的剖面图；
图12为本发明提供的参数优化反演程序的流程图；
图13为本发明提供的jeong实测模拟应变对比的曲线图。
具体实施方式
下面将参照附图更详细地描述本公开的示例性实施例。虽然附图中显示了本公开的示例性实施例，然而应当理解，可以以各种形式实现本公开而不应被这里阐述的实施例所限制。相反，提供这些实施例是为了能够更透彻地理解本公开，并且能够将本公开的范围完整的传达给本领域的技术人员。
本发明实施例提供了一种新型多参数影响下早期混凝土性态变化研究方法，其方法实现方便、可操作性强且使用效果好、实验条件限制小，可有效节约研究成本，且计算结果精度高。采用本发明提出的多参数影响下早期混凝土性态变化研究方法进行相关分析，如下步骤：
(1)依据工程所需选定合适的试验材料，如水泥、砂、石子等；
(2)依据相关规范混凝土试验规范及实际工程试验要求，确定混凝土配合比；
(3)设计有效的试验方案获取相关的混凝土物理参数，如导热系数、密度、比热、线性膨胀系数等。
其中：有效试验方案的设计需遵循以下几点：
①试验方案应严格按照相关规范，并结合实际所需测试参数来确定；
②混凝土的拌合应严格按照真实配合比进行拌合；
③混凝土物理参数测试过程中，传感器的布置必须真实、有效；
④试验过程中测得是混凝土物理参数应真实、有效。
(4)用温度传感器与振弦传感器分别测量混凝土试件内部的温度与应变；
针对采用埋入式振弦应变计进行应变测试的相关附加说明如下：
1)若非特指，本试验选用的是江苏泰斯特电子设备制造有限公司制造的动态信号测试系统tst5912。tst5912动态信号测试分析系统可对应变应力、荷重、速度、加速度、位移、扭矩等物理量进行精确测量和分析，可以配套使用所有电压、电阻、电荷输出型传感器，温度传感器采用江苏泰斯特电子设备制造有限公司制造的电阻式传感器，而应变计是江苏泰斯特电子设备制造有限公司制造的埋入式振弦应变计。
2)如无特别说明，混凝土初凝阶段之前，振弦传感器所测得应变值与实际应变值有较大差异，需要对所得测量值做出一定的修正。认为初凝前测得的应变需要从实测应变中剔除，对于不同型号不同配比的混凝土来说，真实的初凝时间较难确定，再加之滞后变形对测量结果的影响，完全减掉这部分应变可能会造成较大偏差，甚至导致测量值错误，因此本专利中用实测应变减去初凝之前的应变作为下界，将实测应变作为上界，如果数值模拟结果位于上下界之间，则认为结果符合实际。
(5)结合二次开发改进后遗传算法对试验得出的热力学参数进行反演，从而获得较为准确且符合实际情况的特征参数。并将反演后的相关参数与试验测得数值进行对比验证，确定反演结果的有效性；反演的热力学参数主要包括导温系数、导热系数、比热、密度；改进后遗传算法的反演过程如下：参照图1所示：
①随机产生初始种群，个体数目一定，每个个体表示为染色体的基因编码；
②计算个体的适应度，并判断是否符合优化准则，如果符合，输出最佳个体和其代表的最优解，并结束计算；否则转向第三步；
③依据适应度选择再生个体，适应度高的个体被选中的概率高，适应度低的个体可能被淘汰；
④按照一定的交叉概率和交叉方法，生成新的个体；
⑤按照一定的变异概率和变异方法，生成新的个体；
⑥由交叉和变异产生新一代的种群，返回到第二步。
(6)基于abaqus的二次开发平台，开发了适用于早期混凝土温度场、湿度场计算的子程序，结合工程实际工程资料建立相应的结构有限元模型进行早期混凝土温、湿度场数值模拟计算，达到实时模拟的效果。混凝土温度场、湿度场计算的子程序分别为：
温度场计算的子程序为：
①在瞬时温度场分析过程中，umatht(usersubroutinetodefineamaterial’sthermalbehavior)子程序可以用来定义温度场的本构行为以及材料内部的生热过程，可以实现温度场参数随时间或其他变量的变化过程。子程序结构如下：
subroutineumatht(u,dudt,dudg,flux,dfdt,dfdg,
1statev,temp,dtemp,dtemdx,time,dtime,predef,dpred,
2cmname,ntgrd,nstatv,props,nprops,coords,pnewdt,
3noel,npt,layer,kspt,kstep,kinc)
c
include'aba_param.inc'
c
character*80cmname
dimensiondudg(ntgrd),flux(ntgrd),dfdt(ntgrd),
1dfdg(ntgrd,ntgrd),statev(nstatv),dtemdx(ntgrd),
2time(2),predef(1),dpred(1),props(nprops),coords(3)
usercodingtodefineu,dudt,dudg,flux,dfdt,dfdg,
andpossiblyupdatestatev,pnewdt
return
end
下面对umatht中一些用到的关键变量的含义进行说明：
1)需要定义的变量
u，增量步结束时物体单位质量的内能。该变量在增量步开始时被主程序传入到子程序中，且必须在增量步结束时更新该变量。
dudt，增量步结束时，单位质量物体内能对温度的变化率。
dudg(ntgrd)，增量步结束时，单位质量物体内能对温度空间梯度的变化率。该向量的大小取决于ntgrd的大小。在传统经典温度分析中，该变量为零。
flux(negrd)，增量步结束时的热流失量。该变量同样是开始时刻传入子程序且结束时必须被更新。
dfdt(negrd)，增量步结束时，热流失量对温度的变化率。
dfdg(ntgrd,ntgrd)，增量步结束时，热流矢量对空间梯度的变化率，该向量大小取决于ntgrd的大小。
statev(nstatv)，该向量包含全部与求解过程相关的状态变量。在非耦合热分析过程中，状态变量在增量步的开始由主程序传入umatht子程序中。但是，当用户还使用usdfld子程序时，在子程序usdfld中更新的状态变量将会传到umatht中，因为usdfld子程序的调用在umatht之前。另外，在温度～应力场完全耦合分析中，在其他定义单元力学行为的子程序(creep、umat等)中更新的状态变量，也会传递到umatht中，原因也是这些子程序在umatht之前调用。状态变量矩阵的维数是通过材料参数设置窗口中的关键字“*depvar”进行定义，输入数据的大小即为状态变量矩阵的维数。
2)传入的变量
props(nprops)，材料常数矩阵。矩阵中元素的数值对应于关键字“*usermaterial”下面的数据行中的数据。材料常数的个数，等于关键字“*usermaterial”中“constants”常数设定的值。注意，此时“usermaterialtype”应选择“thermal”。
temp，增量步开始时刻的温度值；dtemp，增量步内的温度增值；time(1)，增量步开始时刻的荷载步时间；time(2)，增量步开始时刻的总时间；dtime，时间增量；ntgrd，空间方向的个数；nstatev，状态变量的个数；noel，单元编号；kinc，增量步数。
②在瞬态热分析问题中，用户子程序film(usersubroutinetodefinenonumiformfilmcoefficientandassociatedsinktemperaturesforheattransferanalysis)可以用来定义基于点、单元或面的非均匀表面热交换系数(可以模拟结构拆模或风速对热交换系数的影响)，可以用来定义变化的环境温度，该温度可以是空间位置、时间或温度的函数。针对早期混凝土温度场分析，film子程序是用来模拟计算中温度场的第三类边界条件。子程序结构如下：
subroutinefilm(h,sink,temp,kstep,kinc,time,noel,
1npt,coords,jltyp,field,nfield,sname,node,area)
c
include'aba_param.inc'
c
dimensionh(2),time(2),coords(3),field(nfield)
character*80sname
usercodingtodefineh(1),h(2),andsink
return
end
下面对film中一些用到的关键变量的含义进行说明：
1)需要定义的变量
h(1)，积分点的表面热交换系数，单位jt-1l-2θ-1。如果该参数没有被定义，h(1)将被初始化为零。
h(2)，表面热交换系数相对于积分点上温度的变化率。单位是jt-1l-2θ-2。在某增量步内，非线性方程组的求解过程中的收敛速度可通过定义该值得以改善，特别是当热交换系数依赖于温度。
sink，大气温度。如果大气温度没有定义，大气温度将被默认为零。
2)传入的变量
time(1)，当前荷载步的时间；time(2)，当前的总时间；noel，单元编号。
③混凝土的湿度场模拟可以采用与温度场模拟同样的方法，即利用umatht来模拟湿度场的控制方程，以及湿度扩散系数的变化。此时应当注意，因为湿度场中的湿度扩散系数d与温度场中的导温系数a对应，但是umatht输入的热学参数为导热系数，因此在计算湿度场中，需要将混凝土的密度和比热都设置为常数1.0，因这两个参数在湿度场分析过程中没有任何意义，所以不影响软件对与湿度场的模拟。利用film子程序同样来模拟湿度场中的第三类边界条件，与温度场模拟类似，其中需要定义的材料参数为混凝土表面湿度转移系数h(1)以及环境的相对湿度sink。通过以上分析设置后，原来温度场结果就成为了相对湿度的结果，这样就得到了不同时刻早期混凝土结构的湿度场分布，为计算干缩变形及干缩应力提供了基础。因此，早期混凝土湿度场的模拟同样用到了uamtht和film两个子程序，且其编写过程与温度场基本一致，此处不再重复说明。
(7)依托abaqus二次开发平台，通过调用温度场、湿度场和应力场(温度场及湿度场的子程序同(6))子程序，实现了早期黏弹塑性混凝土的热-湿-力学耦合模型和数值求解过程；再结合热膨胀系数，研究早期混凝土的性态变化规律，实现早期混凝土热湿力的时变性耦合分析。
实用程序sprinc(应力场)：
软件提供的实用程序是软件自带的可以实现某种功能的单独小程序，可以实现对应力、应变不变量等变量的获取。本文在实现对于早期混凝土多场耦合分析及开裂风险计算过程中，用到了sprinc实用程序来获取结构每个时刻的主应力值，因此本文仅对该实用程序进行使用说明，其他实用程序用户可以查阅相应的用户手册。
其具体的使用格式为：
callsprinc(s,ps,lstr,ndi,nshr)
变量说明：
(1)需要提供给实用程序的变量说明：
s，应力或应变张量；
lstr，一个标识符，当其为1时，表明变量s为应力张量，当其为2时，表明变量s未应变张量；
ndi，正应力个数；
nshr，剪应力个数。
(2)通过调用实用程序可以获得的变量：
ps(i)，i的值可以为1、2或3，分别标示第一、第二和第三主应力。
子程序的调用过程如下：
在有限元分析的每个增量步内，abaqus提供的用户子程序至少被调用一次，然而不同子程序被调用的次数并不一样。针对于用于定义材料、单元或者接触行为的用户子程序，为了获得初始刚度矩阵，在每个增量步的初次迭代过程中，每个积分点子程序被调用两次，而在此增量步后续的迭代过程中每个积分点调用一次子程序。
abaqus的用户自定义子程序是按照abaqus软件提供的相应接口，遵守fortran语言编写要求，用户自己编写的程序代码。子程序是一个独立的程序单元，允许单独进行保存和编译，同时可被其它程序单元引用。它的一般结构形式是：
subroutinea(x1，x2，……，xn)
include‘aba_param.inc’(abaqus/standard用户子程序)
orinclude‘vaba_param.inc’(abaqus/explicit用户子程序)
dimension(a1，a2，…..，an)
用户编写代码段
return
end
上述代码中x1，x2，……，xn是abaqus提供给用户子程序a的接口参数，其中有些参数是abaqus主程序传到用户子程序中的，例如subroutineumat中的temp、coords等；有些是需要用户自己定义的，例如f。文件aba_param.inc和vaba_param.inc随着abaqus软件的安装而包含在操作系统中，其中含有重要的参数，帮助实现abaqus主程序对用户子程序进行编译、链接。a1，a2，…..，an是用户子程序用到的与主程序交互的变量标示符号。
在一个有限元分析中，用户可同时应用多个子程序，但是需要注意的是应该把用到的多个子程序编写在用一个以.for为后缀名的文件中。调用用户子程序有两种方法：一种是在abaquscae中运行，打开editjob菜单，通过general子菜单中的usersubroutinefile对话框选择相应的用户子程序文件即可；第二种方法是在abaquscommand窗口中运行，调用语句如下：
abaqusjob＝job-nameuser＝{source-file|object-file}
针对定义单元本构行为的子程序来说，当计算模型中含有多种材料且它们的本构并不相同时，用户可以使用cmname变量实现不同材料对其相应本构子程序的调用。现以umatht子程序为例进行说明：
if(cmname(1：4).eq.'mat1')then
callumatht_mat1(argument_list)
elseif(cmname(1：4).eq.'mat2')then
callumatht_mat2(argument_list)
endif
可以看出，umatht_mat1和umatht_mat2才是分别针对材料mat1和mat2的真正含有定义材料力学行为的子程序，其结构形式与单独的一个umatht子程序完全一致，此处umatht子程序只起到了一个指导选择子程序的作用。
黏弹塑性混凝土的热-湿-力学耦合模型和数值求解过程
结合前述研究，继续以abaqus二次开发接口为平台，研究早期混凝土温度场、湿度场和应力场的三场耦合问题，旨在更加真实准确的反映早期混凝土性态。采用各个物理场之间顺序耦合的方法，实现了早期混凝土三场的时变性耦合。
多场耦合说明：
研究中采用顺序耦合方法实现了温度场、湿度场和应力场三场耦合，其耦合过程为将温度场和湿度场的结果以荷载的形式施加于结构应力场，而同时考虑到结构温度对湿度场的湿度扩散系数有重要影响，因此，本文顺序耦合的顺序是先采用精确模拟方法计算结构温度场，然后利用温度场子程序来模拟结构相对湿度场分布，最后将两者施加于结构进行应力场的计算。
三场顺序耦合的基本流程如图2所示，本文研究的三场耦合分析总是单向的，其具体的实施方案为：
①分析温度场。设置温度场计算参数，如导热系数等，选择温度求解模块并给结构模型施加相应的初始温度，设置合理的温度场求解边界条件，最后调用子程序进行温度场模拟计算。当计算模拟结束时，读取并储存计算结果。读取结果的存储方式是文本文件，且以相应单元的编号给文件命名。
②分析湿度场。如前所述，利用温度场求解程序进行湿度场的模拟计算，将密度、比热设置为1.0，其他参数输入相应湿度场的计算参数。在湿度场计算的子程序文件中，通过自定义应用程序实现了计算某节点读取以该节点编号命名的文本文件的功能，并将节点温度值存储成向量形式，为计算温度对湿度扩散系数的影响做准备。湿度场模拟结束后，与温度场结果一样，将湿度场结果读取并储存。
③分析应力场。当温度场及湿度场分析完成之后，最后利用耦合方法模拟计算应力场。因为，应力场中的材料力学参数，如弹性模量、抗拉强度等与时间和温度均相关，因此温度场和应力场的耦合还是采用传统的方法实现。湿度场和应力场的耦合采用与(2)中温度场耦合到湿度场中的方法类似，在获取计算节点增量步内湿度变化后，可以根据相应公式计算干缩应变，同时因为本文开发的umat子程序的功能之一就是更新应力，因此可以直接将干缩应变叠加到总应变中，实现对干缩应变的考虑。模拟计算结束后，即得到结构的温度干缩应力。
经过以上几步可以实现对三场顺序耦合的模拟，可以实现对混凝土结构早期应力的精确模拟计算。图2涉及的控制方程、边界条件以及粘弹性本构的理论分析如下：
1)控制方程
利用推导温度场热传导方程思想同样可以推导湿度扩散的控制方程。对于混凝土材料，普遍认为是一种多孔复合材料，湿度满足fick第二定律，即认为物体内部可蒸发水流量qx与物体内部湿度梯度成正比。同时，利用质量守恒定律，可以得到混凝土非稳定湿度场的扩散控制方程为：
式中，h为混凝土结构的相对湿度，用百分数表示；t为时间，h；dx(h)、dy(h)、dz(h)分别为x、y、z三个方向的与相对湿度相关的湿度扩散系数，且一般认为三个方向的湿度扩散系数相等；
对于普通硅酸盐水泥混凝土，一般水灰比较大，从而结构内部湿度分布主要是由于湿度扩散引起的，所以水化反应自干燥而引起的湿度损失
式中，d(h)为混凝土湿度扩散系数；其他参数同前。
与温度场理论相似，湿度场控制方程为典型的非线性抛物线偏微分方程，无约束条件下满足方程的解有无数个，因此为了确定特定情况下的湿度场分布情况需要对方程施加一定的初始条件和边界条件。
2)初始条件与边界条件
通过前面的湿度场控制方程可知，湿度场的初始条件和边界条件与温度场的初始条件和边界条件同样应该具有相似性。
①初始条件
混凝土湿度场初始条件，因混凝土开始模拟计算的时间与浇筑时间的对应关系而不同。当湿度场模拟从浇筑开始后就开始模拟时，此时湿度场的初始条件可以假定为100％，即：
h＝h(x,y,z,t)＝h(x,y,z,t0)＝1.0(3)
当混凝土湿度场模拟不是从浇筑后即开始时，应根据实际情况或通过试验测得计算时刻的初始湿度值，并定为初始湿度值。
②边界条件
因混凝土结构各个表面所处的环境不同，湿度场同温度场类似同样也具有四类边界条件。
湿度场计算的边界条件可以分为以下四种：
a)第一类边界条件
混凝土结构表面相对湿度h为时间的已知函数，即有
h＝h(x,y,z,t)(4)
b)第二类边界条件
混凝土表面为绝湿边界，即
式中：n为表面法线方向。
c)第三类边界条件
当混凝土与空气接触时，可用对流边界条件来表示该种边界条件，即混凝土的表面湿度交换量和表面湿度hs与大气相对湿度he之差成正比，其数学表达式为：
式中，f为表面湿度扩散系数，单位为m2/h。
d)第四类边界条件
当两固体接触时，若接触良好，则在两物体接触面上的湿度和湿度量应该都是连续的，边界条件为：
3)湿度场参数的确定
通过以上内容可知，湿度场在模拟分析过程中，需要确定的参数为湿度扩散系数d、表面湿度转移系数f、大气环境相对湿度变化he。
①湿度扩散系数d
湿度扩散系数是描述混凝土结构湿度扩散能力的一个最重要的参数，受环境温度、湿度等因素影响较大，一般需要通过试验确定。在实际模拟过程中，研究人员主要是用ceb-fip(1990)规范中提出的计算公式，该公式是相对湿度h的函数：
式中，d1为湿度扩散系数的最大值，此时h＝1.0；α＝d0/d1，d0为h为零时湿度扩散系数的最小值；n为曲方程的拟合系数；hc为湿度扩散系数为最大值的一半时的相对湿度。建议，α、hc、n可以分别取0.05、0.80和15。d1也可由下式确定。
式中，d1.0取3.6×10-6m2/h；fck0取10.0mpa；混凝土的抗压强度fck可以由平均抗压强度fcm来估算，可取fck＝fcm－8.0mpa。
同时，由于早期混凝土温度场变化比较剧烈，因此必须考虑温度对湿度扩散系数的影响，同时温度对参数的影响满足arrhenius方程，所以为了考虑温度对于湿度扩散系数的影响，引入基于arrhenius方程的温度影响函数。同时，因为早期混凝土温度变化幅度大，为了更好的模拟实际情况将时间利用等效龄期代替。综上，温度场和湿度场共同作用下的湿度扩散系数表示为：
其中，f(te)为考虑龄期影响函数，其表达式见式(12)；其他参数同前。
式中，ead为混凝土扩散活化能，本研究中取值为35kj/mol。
②表面湿度转移系数
表面湿度转移系数也称为表面水分交换系数，影响该参数的因素主要是水灰比、湿度梯度、温度和风速。表面湿度扩散系数可以采用更加符合实际的改进menzel表达式：
f(h-he)＝a(0.253+0.06va)(h-he)(13)
式中，a为经验系数，取决于混凝土水灰比；va为平均风速，m/s；h、he分别为当时结构表面和大气的相对湿度值。
4)应力场计算基本理论
混凝土的强度增长与水泥的水化度直接相关，而水泥水化度不仅依赖于龄期，同样依赖于温度。温度越高，水化速率越快，弹性模量和强度增长就快；反之，则弹性模量和强度增长缓慢。同时，混凝土由于自身水化温升受环境的影响，使得各点温度值在同一时刻存在一定的差异，这就导致了相同龄期而不同部位混凝土强度的差异。通过引用基于等效龄期概念的成熟度方法，同时考虑龄期和温度共同对混凝土成熟度的影响，能够全面反映混凝土强度随龄期的变化情况，从而弥补了传统方法的不足。
混凝土结构对荷载的反应是非常复杂的，表现为三种基本类型的变形：弹性、塑性和粘性或者是它们的组合，如弹塑性、粘弹性。从时间方面考虑，可以将混凝土的变形分为两大类：与时间无关的变形和与时间相关的变形。同样地，从应力方面考虑，可以将变形分为应力相关和与应力无关两类变形形式。瞬时变形(在很短时间发生)是指与时间无关的变形，它代表了荷载施加给混凝土结构的瞬时变形反应。徐变变形是指随时间变化的变形，而不是荷载施加后结构瞬时产生的变形。
混凝土结构在时间τ的总应变ε(τ)可表示成：
ε(τ)＝εe(τ)+εc(τ)+εt(τ)+εs(τ)(14)
式中：εe(τ)为应力引起的瞬时应变，一般可认为是线弹性的；εc(τ)为混凝土徐变应变，与应力、加载龄期及持载时间等有关；εt(τ)为混凝土由温度变化而引起的变形；εs(τ)为混凝土的收缩变形，分干燥收缩和自收缩两部分。
式(14)所表示的混凝土结构的总应变中，εe(τ)和εc(τ)是由应力引起的，而εt(τ)、εs(τ)及εg(τ)的产生与当时的应力状态无关。
由以上可知，混凝土结构在荷载作用下，除了产生瞬时弹性变形外，还会产生徐变变形和收缩变形等。同时，注意到影响混凝土收缩和徐变的因素是共同的，如水泥类型、水泥细度、骨料种类、水灰比、养护方式等。因此，本发明在对混凝土结构早期应力分析过程中一并考虑了收缩、温度和徐变变形对结构应力的影响。
5)早期混凝土力学参数
在混凝土浇筑后的几天时间内，其力学参数的发展、徐变、干缩以及温度变化等之间存在复杂的相互作用。现阶段对于早期混凝土主要通过试验以及模型计算的手段来研究其强度及弹性模量的增长模式。
要想更加精确的预测早期混凝土应力增量以及开裂风险情况，准确的力学参数输入是必不可少的。
在t时刻，混凝土强度与水泥类型、温度以及浇筑温度等条件相关。在平均气温为20℃环境下，ceb-fip1990规范给出了混凝土抗压强度和弹性模量随时间变化的计算公式，分别见式(15)和式(16)。
式中，ecm(t)、fcm(t)分别为t时刻混凝土的弹性模量和抗压强度，(pa)；e28、fc28分别为混凝土第28天的弹性模量和抗压强度，(pa)；t1，为一天，d；s为与水泥品种相关的材料系数，取值在0.20～0.38之间。
为了使混凝土强度增长模型更加适合于早期混凝土，kanstad等基于混凝土早期试验及等效龄期成熟度概念，对ceb-fip1990规范中提出的混凝土强度增长公式进行了修正。修正的方式是通过引入t0参数来表征混凝土强度开始增长的时刻，修正后的混凝土强度增长公式见式(17)～(19)，修正后的强度增长公式与早期混凝土试验吻合较好。从公式中可以看出，不论是弹性模量还是抗压(拉)强度，其t0参数是相同的，这样为由试验结果获取该参数提供了方便。如果水泥水化反应过程采用之前介绍的水化度来进行描述，并定义与时刻对应的水化度为α0，那么t0参数的含义就是水泥必须首先进行部分的水化反应之后混凝土才会产生可以被“测量”的力学属性。
式中，e28、ft28和fc28分别是28d混凝土弹性模量(pa)、抗拉强度值(pa)和抗压强度值(pa)，t0为混凝土强度开始发展的时刻(h)，672即28d的小时数(h)，参数s、nt、ne与混凝土的类型有关，可通过试验或参考相关文献取值。
对于混凝土材料来说，泊松比被认为是一个相对稳定的参数，一般假定泊松比为0.2适应于几乎所有品种的混凝土。但是，deschutter建议对于早期处于硬化过程中的混凝土采用下式确定其泊松比：
式中，υ为混凝土泊松比；α为当前水化度。
从上式中可以看出，当混凝土水化度为零时，混凝土泊松比为0.5，这与新浇混凝土变形不会导致体积变化的事实是一致的。本文在计算过程中，为简便计算将混凝土泊松比取为定值0.2。
6)线性徐变理论
虽然早期混凝土的徐变现象本质上是一个非线性变化过程，但是国内外文献中建议采用的徐变模型基本上都是建立在徐变过程线性变化的假设之上，即在恒定应力下的徐变应变被认为与应力水平是线性相关的。在混凝土结构正常使用范围内，早期混凝土可以被认为是线性粘弹性材料。与线性的假设相一起，叠加原理被应用于早期混凝土徐变计算过程中。
混凝土线性叠加原理的含义是：在任意时刻t，混凝土由在t'时刻应力增量(拉伸或压缩)产生的应变不受任何应力的影响的效果，无论是应用较早或较晚的时间t'，但仅仅依赖于t'。该原理已经被广泛的应用到混凝土结构设计和研究过程中，并被认为是有效且方便的假设。
因此，单轴应力下结构时间t之前因所有的应力增量产生的应变之和，可以写成以下的徐变变化规律：
式中，t为混凝土浇筑后的时间，d(h)；τ，为加载龄期，d(h)；dσ(τ)，应力增量，pa；ε(t)为t时刻的总应变；ε0(t)为与应力无关的应变；j(t,τ)为徐变函数。
7)dpl徐变模型
如前所述，对于早期混凝土其徐变变形占有很大比例，因此混凝土的应力分析必须包含一个能够真实预测其应变和应力场的徐变模型，从而做到对开裂风险的较准确估计。混凝土徐变由基本徐变(basiccreep)和干燥徐变(dryingcreep)组成。前者是混凝土与周围环境之间无水分交换下的徐变；后者是结构承受荷载后，因结构处于干燥环境而产生的徐变。因为在早龄期干燥蠕变只在混凝土表面1-2cm深处发生，因此对干燥蠕变结构应力的影响不大，可以忽略。因此，本发明在研究早期混凝土过程中只考虑了基本徐变部分。
通过国内外学者的研究，现阶段已经提出的很多从混凝土本构角度来预测其徐变变形的模型，并且根据它们的形式可以大体分为两大类徐变模型。在第一类模型中，徐变函数被看作是加载时间和荷载持续时间的函数，这一类徐变模型函数的一般表达形式为：
式中，e(τ)表示加载时刻的混凝土弹性模量；
另一类徐变函数模型是基于混凝土变形是由两部分(或更多部分)组合而成的思想建立，即由一个可恢复的(延迟弹性)的变形和一个不可恢复的流动变形组成，其一般表达式见下式：
为了模拟特定混凝土在一外界环境下的徐变及收缩特性，国内外学者已经开发研究开发了多种徐变模型，它们之间的不同在于其精确性、简易性以及适应的条件，且通常某个模型是其他模型的修正结果。20世纪80年代，bazant和panula开发了双幂徐变模型(doublepowerlaw，dpl)，它是bazant-panula的简化形式，该模型基于线性叠加理论，可以很好的模拟早期混凝土的基本徐变，其表达式见下式：
式中：t为混凝土龄期；τ为混凝土加载龄期；j(t,τ)为徐变函数；e0为渐近弹性模量；α为长期徐变模型参数；q、d、p为徐变模型参数。
式(24)被称为双幂徐变函数是因为模型中关于龄期和持载时间都是幂函数的形式。e0为混凝土材料参数，称为渐近弹性模量，它比通常的弹性模量数值要大很多。参数q，d和p是在参考温度(23℃)下的基本徐变模型材料参数。这些系数可以相对简单地通过最优化技术依据试验数据确定。同时，该表达式中的时间t和τ都以d作为单位。
从表达式(24)中可以看出，τ-d部分对基本徐变的贡献方面应该注意以下几点：因为双幂函数中时间以天作为单位，因此当加载龄期为1天时参数d对于徐变没有影响。对于其他的加载龄期，当t小于一天时，τ-d随着d的增大而增大，当t大于1天时，τ-d随着d的增大而减小，如图3a所示。从图中可以看出，当加载龄期小于1天时，该部分对徐变的影响很大，当时间大于5天时，该影响基本为线性的。
参照图3b展示了参数p的大小对于(t-τ)p项大小的影响。不同参数p在持载时间大于1天时对徐变的影响差别非常明显，当持载时间小于1天时，该影响变小。图3a、3b的两个曲线图说明了d、p参数变化对于徐变发展的影响。式(24)更加适合于时间跨度较长的徐变模拟，其较简单的一种变化形式更加适合于模拟早期混凝土徐变效应，并公认可以得到较好的模拟结果，其表达式见式(25)。最近，一些关于早期混凝土压缩和拉伸徐变行为的实验结果与dpl计算模拟都显示了良好一致性。此外，bazant和prasannan认为dpl本身含有可以模拟徐变持续时间较小下的变形。因此，本发明利用dpl模拟早期混凝土徐变过程是合理的。
式中，τ为加载龄期，d；e-1(τ)为加载时刻的弹性模量，pa；t-τ为持载时间，d；q、d、p为材料参数，可通过混凝土试验或参考相应规范获得。
8)收缩变形与收缩应力
混凝土由收缩应变在存在约束情况下而引起的收缩应力可以使结构应力增大，存在导致混凝土开裂的可能性。因此，在早期混凝土应力场分析中，为了分析更加精确和真实，有必要考虑收缩应力。
(1)收缩变形
混凝土收缩是指混凝土由于所含水分变化、化学反应等因素而引起的体积缩小。根据产生收缩的机理及条件，混凝土的收缩变形主要分为浇筑初期(终凝前)的塑性收缩、硬化过程中的干燥收缩、自生收缩以及因碳化反应而引起的碳化收缩。
a)塑性收缩
塑性收缩是指在混凝土开始搅拌后，由于此时水泥的水化反应而出现的沁水和水分丢失现象，一般发生在混凝土搅拌后3～15h。因该塑性收缩发生在混凝土初凝到终凝之间，而混凝土此时所处状态为塑性状态，故称这种收缩变形为塑性收缩。塑性收缩变形的大小约为体积的1％，其量值较大，然而该阶段混凝土处于有塑性到硬化的变化过程，弹性模量刚刚开始发展，其值较低，因此在早期混凝土应力场应力及开裂分析中，由塑性收缩而引起的变形可以忽略。但是，当混凝土养护不良时，会导致结构表面由于塑性收缩而引起的开裂现象，对结构整体安全性构成隐患，因此对于混凝土结构，尤其是厚度较薄的结构(如地下室墙体)，必须采取足够的工程施工措施，否则塑性收缩引起表面裂缝的可能性很大。
b)自生收缩
自生收缩是指混凝土在恒温和绝湿条件下由于水化反应消耗水分而引起的体积缩小变形。这种收缩一般发生在大体积混凝土的内部、钢管混凝土等。在已经硬化后的水泥浆体中，还未发生水化反应的水泥继续进行水化反应是自生收缩产生的主要原因。混凝土自缩与水泥品种、水泥用量及掺用混合材料种类有关，其值可能为正，也可能为负。对于普通硅酸盐水泥混凝土其值一般为正值，因此对于混凝土抗裂不利，但其值一般较小。
c)干燥收缩
干燥收缩是指混凝土结构浇筑后，在不饱和的空气中失去混凝土内部毛细孔水、吸附水等而发生的不可逆收缩变形，其属于物理收缩变形。混凝土从刚刚开始浇注的塑性状态到初凝再到硬化阶段过程中一直伴随着干燥收缩的发生，收缩变形一般在(200～1000)×10-6范围内，基本上是自生收缩变形的10倍，干燥收缩变形在总收缩变形中所占比例最大。
d)碳化收缩
碳化收缩是指大气中的二氧化碳，在有水分湿度适宜(约50％)的条件下与混凝土中水泥的水化物(如硅酸钙凝胶等)发生化学反应产生caco3和游离水等而引起收缩变形，这种反应称之为碳化反应。由于碳化作用而产生的体积收缩就叫碳化收缩，具有不可逆性，且其一般仅发生在混凝土表面较浅的范围内。碳化收缩变形在一般环境的分析计算中可不予考虑，只有混凝土结构在特定环境中的特久强度和表面开裂研究分析中才加以考虑。
(2)混凝土收缩应力计算
对于收缩应变与相对湿度之间的关系，不同学者对其研究的结果并不完全相同。但是，收缩变形与空隙相对湿度关系密切。因此，可以通过建立收缩变形与相对湿度的关系来求解由于湿度变化而引起的收缩应变。bazant认为干缩引起的应变与相对湿度的损失之间存在一定的关系，其表达式为：
式中，ksh为收缩系数；
由以上分析，可以得出混凝土结构收缩应变变化增量与混凝土结构相对湿度变化之间的关系表达式为：
δεs＝kshδh(28)
其中，
一般工程模拟收缩应力过程中，采用将收缩变形换算成当量温差(式(29))的形式加以考虑，但是本文在考虑收缩变形时，直接将单元的收缩变形增量“累加”到主程序传入到子程序umat的应变增量中，从而较方便地考虑收缩变形对应力的影响。
式中，δt为收缩当量温差，℃；α为混凝土的热膨胀系数。
9)粘弹性徐变子程序
因为混凝土在浇注早期应力变化相对较大，所以dpl的数值计算的实现是通过将dpl展开成泰勒级数的形式，进而在时间增量步内获得方程(25)的近似解。此处将徐变函数展开成泰勒级数，而不是采用通常的做法将函数展开成狄利克雷级数，其原因是展开成狄利克雷级数的方法通常应用于时间跨度较长的徐变效应模拟，时间跨度可以是几年或几十年，而展开成泰勒级数的方法在模拟短时间徐变效应具有优势。另外，这两种展开方法都有一个明显的特点，就是可以通过在每一步更新的某个定义的状态变量，这样可以使下一步的应力增量δσ显示地表达为应变增量δε和应力(应变)历史的函数，这一点使得有限元模拟成为了现实。
(1)公式推导
本文对早期混凝土的温度应力研究是基于abaqus二次开发平台，在求解非线性问题时采用newton-raphson算法，即采用增量步的方法来求解非线性问题，如图4所示。考虑到混凝土早期的粘弹性行为，如上所述将徐变函数展开成级数形式，并推导应力增量的显示表达式，以及每一步需要更新的状态变量方程，其具体推导过程如下：
为了求得在t时刻的基本徐变变形，已知徐变函数j(t,τ)，则应变张量ε(t,τ)可以表达为：
式中，[d]——为弹性矩阵；σ(τ)——龄期τ处应力张量；ε0(t)——与应力无关的应变张量。
在以下对公式(30)进行变形推导过程中，假定徐变满足线性叠加原理，在每个时间增量步ti至ti+1内，假定应力为线性变化，即应力的导数为应力增量与时间增量的比值。每一个增量步内混凝土弹性模量采用等效弹性模量代替。此外，利用积分中值定理简化积分方程，进而可以推导出双幂徐变函数泰勒展开式的应力增量本构方程。
在时间增量δt＝ti+1-ti内，相应的应力张量和应变张量增量分别为：
[δσ]＝[σ(ti+1)]-[σ(ti)](31)
和
[δε]＝[ε(ti+1)]-[ε(ti)](32)
利用公式(30)，则应变增量可以表示为：
为了将应力增量显示地表达为应变增量和应力(应变)历史的函数，假定应力在增量步内为线性变化，即应力的导数为应力增量与时间增量的比值，利用中值定理(ti≤t*≤ti+1)，定义等效弹性模量为：
利用式
令
式(36)就是早期混凝土的粘性应力部分。至此，我们获得了有限元计算中使用的线性粘弹性本构方程，其表达式如下：
如前边所述，为了更好的模拟早期混凝土的徐变效应，将徐变函数展开为泰勒级数形式为：
假设n＝5，则有
h0＝a0+a1(t-td)+a2(t-td)2+a3(t-td)3+a4(t-td)4+a5(t-td)5(39)
h1＝-a1-2a2(t-td)-3a3(t-td)2-4a4(t-td)3-5a5(t-td)4(40)
以此类推通过系数a系列以及函数(t-τ)p在τ＝td处的展开式可以求得参数h2、h3、h4、h5。其中，
因为，
此时，通过带入展开后的徐变函数方程，该增量步内的等效弹性模量表达式变为：
相应地，粘性应力部分变为：
式(43)中，每个增量步结束时都需要更新的核心变量(状态变量)
将式(42)至式(44)代入式(37)就可以得到适应于早期混凝土的泰勒级数形式的粘弹性本构模型。
(2)子程序编写注意事项及其调用过程
粘弹性徐变本构模型采用用户子程序umat进行编写，其编写的核心语句是定义材料的应力应变关系、状态变量的定义和更新以及应力更新。在编写过程中，一定要注意主程序提供给子程序的应变张量为柯西应变张量，所以定义雅克比矩阵的时候应予以注意。要对状态变量的更新本质过程有清楚的了解，abaqus主程序在用户没有特定定义状态变量数值的前提下，为每个状态变量默认初始值为零，并且主程序会给产生一个状态变量的备份文件，每次运行子程序都会对状态变量进行计算更新，但是如果增量步中该迭代步没有收敛，在下一次迭代中，abaqus主程序会利用之前储存的变量备份进行计算。应力更新要注意6个应变分量的排列顺序，以及循环的编写。除了以上提到的情况之外，用户应该熟悉abaqus求解的过程，子程序被调用的频率和次数，这样才能全面的了解子程序运行的环境，有助于程序的正确编写。
umat子程序和abaqus主程序之间的交互过程为：从平衡时刻tn开始，主程序会在下一个增量步的第一次迭代中调用子程序两次，第一次的目的是为了获得当前的刚度矩阵，然后主程序会根据当前增量步的荷载增量大小以及刚度矩阵提供给子程序一个应变增量δε，子程序会根据这个应变增量，通过本构关系计算并提供给主程序新的柯西应力增量σ(tn+δt)，然后主程序会根据最新的应力验证是否与当前外荷载平衡。如果计算的应力结果收敛，那么程序会进入下一个增量步进行计算，并根据收敛情况选择下一个增量步的步长(在未固定步长的情况下)，否则主程序会首先对应变增量进行调整，进行再一次的平衡迭代计算。
图5说明了umat子程序于主程序是如何协同进行计算的，计算过程大体分为四步：
第一步：在平衡时刻tn，主程序向子程序提供总时间，时间增量δt，总应变εtotal(tn)和应变增量δεtotal(tn)，当然平衡时刻应力σ(tn)已由上一步计算求得。这些变量信息将传入子程序用来计算新的柯西应力张量σ(tn+δt)。
第二步：根据上一步更新的状态变量，计算本增量步的粘性应力
第三步：abaqus主程序更新应变张量εtotal(tn+δt)＝εtotal(tn)+δεtotal(tn)，产生新的总应变εtotal(tn+1)。
第四步：abaqus根据新的应力、应变与外荷载进行平衡迭代计算，默认最大迭代次数为16次。在最大迭代次数内达到收敛之后，会进行下一个增量步的计算。反之，在运行情况下主程序会自动减小增量步大小，再次进行计算分析，默认最大增量步次数为5次，若还是不收敛计算将会终止并报错。
本专利实现过程中步骤(4)、(5)的具体实施方案如下：
1)早期混凝土试件内部的温度与应变测试方法为：
①混凝土浇筑前，在浇筑模具内表面接缝处涂抹约2mm厚硅胶，防止浇筑完成后混凝土内水分的流失；
②将聚苯泡沫板贴在模具内表面，聚苯板接缝处密封。后期浇筑完成之后，试件顶部同样用聚苯板密封，起到保温、保湿的作用，从而减小早期混凝土因干燥收缩和塑性收缩对结构应力场的影响；
③振弦传感器和温度传感器固定在模具内部，其中一振弦传感器记为vws-t；另一个记为vws-b(详见附图11a、11b)；温度传感器从上往下依次记为ts-1、ts-2、ts-3(详见附图11a、11b)；
④称量所需材料的质量，采用搅拌机搅拌，搅拌均匀后将混凝土浇筑在模具内部，立即使用振捣棒振捣使试件振捣密实，以排除气泡；
⑤通过启动测试系统开始采集试件的温度与应变数据。
2)早期混凝土温度场计算参数反演流程图，详见图12。反演过程中相关实施过程及实例验证如下：
调用子程序进行abaqus有限元计算主要涉及的方法及理论有：
对于早期混凝土温度场的精确模拟需要用到umatht和film两个用户子程序。其中，umatht子程序用来模拟混凝土的水化放热过程，以及考虑早期混凝土比热和导热系数的变化过程；film子程序用来模拟混凝土与空气接触的第三类边界条件，可以用来模拟表面拆模或温度骤降等。
1)uamtht子程序
在对早期混凝土温度场进行模拟分析时，存在以下几个abaqus本身无法模拟实现的典型特点：
(a)早期混凝土受水泥水化反应影响，会自身产生大量的热量，导致温度升高，而混凝土生热率不仅与时间有关还与温度有关，即与水化度相关；
(b)比热和导热系数同样依赖于水化度而变化，不仅仅是时间的函数；
(c)水化度的计算依赖于等效龄期，而归根结底依赖于当前的温度和时间。
因此，基于以上特点在利用abaqus进行早期混凝土温度场模拟时，需要用到软件的二次开发接口。abaqus一共提供了两种可以模拟水化反应过程的子程序，分别是umatht和hetval。hetval子程序可以实现对水化反应过程的模拟，但是其无法模拟早期混凝土比热等参数的变化，而umatht作为定义温度场本构的子程序是可以实现所需的全部功能。因此，本文在模拟早期混凝土温度场时，采用umatht子程序。该子程序的工作流程见图6。
uamtht子程序编写的核心语句就是定义增量步结束时刻混凝土单位质量的内能u和热流失量flux。该子程序可以考虑变化的比热和热传导系数，从而可以更加精确的模拟温度场。其简化形式可以模拟绝热温升过程，当采用等效方法考虑冷却水管时，子程序也可以模拟冷却水管的作用。在编写过程中，应特别注意time变量的含义，如在本子程序中time(2)表示当前增量步开始时的总时间。当有多个求解步时，注意区分当前step的时间和总时间的表达方式。
在温度场有限元模拟计算过程中，t时刻每个积分点的实际温度tpj可以通过umatht子程序中的temp关键词获得。
2)film子程序
film子程序编写的核心语句就是定义混凝土表面与外界的热交换系数h(1)，以及定义大气温度随时间的变化过程方程sink。用户可以通过判断sname的名称对不同的外表面定义不同的随时间变化的热交换系数来模拟混凝土结构拆模或因其他原因引起的放热系数的变化。通过定义sink可以模拟寒潮引起的温度骤降等气候条件。
在温度场子程序编写过程中，需要注意一下几点：用户子程序umatht只能用于完全积分的单元类型，不能应用于减缩积分的单元类型。单位要统一，能量单位可以统一用千焦(kj)，也可以统一用焦(j)；特别注意几个用到绝对温度的公式；最后需要特别注意step--timeperiod中设置的参数是代表小时还是代表天的概念。
在绝热温度场模拟过程中，只需要用到uamtht子程序即可。当需要模拟实际的温度场分布时，需要编写uamtht和film两个子程序，并需要把它们保存在一个.for中，在job-general中调用即可，或使用commang窗口调用计算。
本文在分析早期混凝土不同温度场计算模型对温度场分析的异同时，也使用了实际情况下单位时间内单位体积混凝土的放热速率公式，其大致推导过程如下：
根据实际情况，图7中用曲线分别表示了绝热情况和实际情况下的混凝土放热和升温过程。在图中任取时间增量步tj到tj+1，在时间增量步δt内，设绝热情况下放热量为δqaj，实际情况下放热为δqpj。当时间步增量较小时，有：
由公式(45)，并考虑到一般混凝土试验在温度为20℃～25℃，上式可以写为：
由此可求得，
由以上推导，可以确定实际状态下混凝土生热率与绝热状态下生热率关系为：
式中，qpj和tpj分别代表实际状态下混凝土的生热率公式和绝对温度值，qaj和taj分别代表绝热状态下生热率公式和绝对温度值。
采用最小二乘法估算热力学参数，即将有限元计算的温度值与现场实测值的误差平方作为参数反演的目标函数j(x)，即
式中：t为测点温度计算值；tm为测点温度实测值；m为测点个数；n为测点次数。
遗传算法(ga)是一种全局随机搜索算法，通过借鉴生物界自然选择和自然基因机制，可以模拟自然界生物从低级向高级的进化过程。本遗传算法拥有2个突出特点:一是收敛速度快。一般的ga拥有收敛速度慢、易陷入局部最优的缺点，改进了遗传算法形成自适应的遗传算法(aga)，通过交叉概率和变异概率随适应度的大小及分散程度进行自适应调整，遗传算法收敛速度快，计算结果更接近全局最优解。主要通过式(50)、(51)实现，如下所示:
式中：pc表示交叉概率；pm表示最大变异概率；fmax为群体中最大的适应度值；favg为每代群体的平均适应度值；f′为要交叉的两个个体中较大的适应度值；f为要变异的个体的适应度值；k1、k2、k3、k4为自适应控制参数，取为k1＝k2＝0.9，k3＝k4＝0.1。
二是编码方式采用的是实数编码。实数编码比二进制编码节省遗传操作的时间，而且在变量定义时，由于实数编码的范围广且表示精度高，可以取大定义域甚至是未知定义域，提供的自适应机制更强。实数编码的实现如式(52)所示:
式中xi,max、xi,min分别为底i个设计变量可取的最大值、最小值；m为控制精度的整数；nrand为小于m的随机数。
参数反演实例验证：
基于遗传算法的方程参数优化为了实现早期混凝土温度场的真实模拟，完全对应实测过程，根据实测资料的模型尺寸和测点布置(图11)，运用数值模拟建立一个长×宽×高分别为0.6m×0.3m×0.9m的混凝土块模型，仪器监测每一时刻各测点的温度值。反演选用的绝热温升公式为式(53)，在遗传算法相关参数的定义中，种群规模设为20，初始交叉、变异概率分别为0.8、0.3，遗传代数为50代。根据相关文献的研究，反演参数θ0、a、b的范围分别取40～50、0.01～0.02、0.85～1.5。在混凝土温度场的模拟中，由于模拟时间较短，为了简便起见，取大气温度为定值24℃，假设混凝土块1次浇筑完成。混凝土浇筑温度为25.58℃，表1为优化程序参数反演的计算结果。模型的边界条件为第3类边界条件，混凝土密度2400kg/m3，比热根据式(54)计算，硬化后混凝土导热系数为9.185kj/(m·h·℃)，考虑水化度的导热系数为9.185×(1.33－0.33α)kj/(m·h·℃)，与空气接触的表面热交换系数为45.0kj/(m2·h·℃)，其余各面的等效热交换系数经相关公式计算所得，为82.64kj/(m2·h·℃)。温度场计算时间跨度为浇筑后的72h。反演中选用的绝热温升公式及比热公式分别为：
绝热温升公式：
式中：θ0为绝热温升；a、b为常数。
比热公式:
c＝(wcαccef+wc(1-α)cc+waca+wwcw)/ρ(54)
式中：ccef＝8.4tc+339，c表示当前的混凝土比热值，kj/(kg·℃)；wc,wa,ww分别为每立方米水泥骨料、骨料和水的质量，kg；cc,ca,cw分别为水泥，骨料和水的比热值，kj/(kg·℃)；ρ为混凝土密度值，kg/m3；ccef为混凝土水泥的假定比热值kj/(kg·℃)；α为混凝土水化反应水化度；tc为当前温度，℃。
表1参数反演计算结果
从反演计算结果表可以看出，当迭代次数达到50次时，在所采用的混凝土配合比下，绝热温升参数θ0和a、b分别为44.46、0.0181、1.3。图8～10分别给出了72h内ta、tb、tc测点温度的实测值与基于反演的绝热温升参数θ0和a、b所得的模拟值。
从图8～10实测值与反演值的对比可以看出，反演计算温度与实测值非常接近，验证了程序二次开发的合理性。
通过本文开发的新型多参数影响下早期混凝土性态变化研究方法，可以精确的模拟混凝土结构内外部的温度值和湿度值，进而研究在温度、湿度综合作用下，结构物的受力和变形特性。研究发现结构物在不受重力产生裂缝的过程中，温度对结构内部的贯穿性裂缝起到了主要作用，湿度对结构表面的微小裂缝其主要作用。为了防止早期混凝土裂缝的产生，通过控制内外温度和湿度梯度、改善结构约束条件及提高混凝土抗裂性能方面进行研究。
举例来说：
(1)实测数据
本文选用的测试系统10/s数据，限于篇幅，仅给出以10h为单位的温度实测值和应变实测值，传感器测量结果分别如表2、3所示。
表2不同时刻温度传感器测量结果
表3不同时刻振弦传感器测量结果
(2)试验方法与反演
混凝土的导温系数用来衡量混凝土本身对热扩散能力的一个参数，混凝土的导温系数越大，在受热或者冷却时，混凝土内各处温度越容易达到均匀，混凝土导温系数数值上等于其导热能力与储热能力的比值，反演结果见图8～10。
混凝土结构温度精确模拟后，主要考虑温度对混凝土湿度扩散系数的影响，数值仿真方面通过温度场求解模块模拟结构湿度场，试验方面通过一种埋置式水泥基材料内部相对湿度测量装置及其测试方法，该装置不仅能测量成熟混凝土的内部相对湿度，还适用于早龄期混凝土，观测和记录结构湿度的变化过程。
应力场的模拟中，采用自主开发的适用于混凝土的黏弹性力学模型，同时需要考虑变化的热膨胀系数。其中，热膨胀系数是物质在热胀冷缩效应作用下，几何特性随温度而发生变化的规律性系数，本专利黏性阶段选用适用于早期混凝土的双幂徐变函数模型和对室内试验数据进行了较好模拟的jeong变化的热膨胀系数公式，结果对比见图13。
显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。
- 还没有人留言评论。精彩留言会获得点赞！