HETVAL子程序？？

混凝土水化热温度场分析其实是相当于在混凝土的温度场分析中加入了一个热源，而这个热源的放热量是随着时间变化的。由于在Abaqus中没有直接功能来模拟随着时间变化的热源，所以需要借用HETVAL子程序来实现随着时间变化的热源功能，并将其耦合到混凝土温度场的计算之中。

见图1，HETVAL子程序用来提供传热分析模型的热源，这个热源是随着时间变化的，这个与混凝土随时间变化的水化放热曲线是一致的（图2）。另外，还可以将该热源在HETVAL子程序中定义为受结果状态变量影响而变化。这个特别重要，因为混凝土水化放热的曲线随着温度的增加会放热加快，而需要模拟这一现象就要考虑温度场计算的温度结果对混凝土水化放热的影响。

图1 abaqus帮助手册关于HETVAL子程序的解释（http://wufengyun.com:888/books/sub/default.htm）

图2 混凝土水化放热速率曲线 （图片来自Fairbairn, Eduardo M. R. , and M. Azenha . "[RILEM State-of-the-Art Reports] Thermal Cracking of Massive Concrete Structures）

HETVAL的主体程序部分如图3，其中子程序调用的形参解释如下：

图3 HETVAL主体子程序

01

子程序传递的变量

• CMNAME: 用户自定义的材料名字。

• TEMP为含2个元素的数组，TEMP(1)为当前温度，TEMP(2)为温度增量。

• TIME为含2个元素的数组，TIME(1)为增量结束时的时间步长，TIME(2)为增量结束时的总时间。

• DTIME为时间步长。

• PREDEF(*)为包含所有用户自定义场变量的数组。

• DPRED(*)为预定义场变量的增量数组。

02

子程序需要定义的变量

而FLUX(1)及其重要，间接表示定义的混凝土热源，如图4所示。FLUX(1)表示产生的热流量，r，也就是单位时间单位体积的产热量，单位为J/T/L^3，即为焦耳/(秒*立方米)。

图4 FLUX(1)表示产生的热流量

03

子程序需要更新的变量

STATEV(*)：用户自定义的结果依赖的状态变量数组。仅在热传递分析中，STATEV在增量步开始时会传入子程序，并且在增量步结束时传递回来，即更新状态变量。

图5 依赖于结果的STATEV状态变量

04

水化放热公式及实现程序

下面通过一个例子来探索下HETVAL子程序模拟混凝土水化热形成的温度场。

水化放热曲线很复杂，工程中一般采用简化的经验拟合表达式，比如混凝土水化放热的公式我们这里采用复合指数式，如式1所示，因为很多文献指出用该水化放热曲线模拟的温度场与实测的温度场模拟较好。Q0为混凝土水化放热的最大值，取为364000J/kg（一般试验或文献里给的都是kJ/kg这种单位，但Abaqus里的标准单位为J，这里要统一单位），而b、c这两个参数取值可在文献中查到，不同水泥类型的参数取值不一样，对于52.5普通硅酸盐水泥，这里b取为0.36，c为0.74。注意到式1中的t的单位为天。

式1：

然后对式1求导，可得到单位立方米混凝土单位时间的水化放热率，qv (J/m^3/hr)，见式2。这里将上式中的时间单位改为了hr，因为一般混凝土温度场监测以小时为单位进行记录。这里的qv就是HETVAL中的FLUX(1)。具体HETVAL子程序如下：

式2：

在该子程序中最重要为热源放热率FLUX(1)，并且及时用statev结果状态变量来存储，这样方便在后处理中输出查看，我称之为实时调试，也就是可以实时监控自己定义或想要输出变量在有限元计算中随时间的变化规律。这里statev(1)记录的是热源放热率flux(1)（时间单位为小时），statev(2)记录的是时变的水化放热总量（注意到时间单位也换为了小时，因为模型中的混凝土热工性能参数，比如热传导系数用的时间单位也是小时，这里统一对应起来），而statev(3)、statev(4)和statev(5)分别为Q0，b，c，这样可以验证下子程序的正确性，看下输出这些确定值是否在计算过程中发生了变化。

另外，特别需要注意的是，下列程序中的Cemass这个变量一定要出现在FLUX(1)的计算表达式中，否则计算出来的温度变化极其小。这并不是说模型是有错的，而是因为Cemass表示每立方米混凝土中的水泥质量（kg/m3），因为HETVAL子程序中出现的FLUX(1)表示单位体积单位时间的水化热源放热量，强调是单位体积，也就是单位立方米混凝土，所以这里的水泥放热量应该是单位立方米混凝土中所有水泥的累加水泥放热量，否则如果不乘以Cemass，则计算的单位质量水泥（kg）的水化放热量，因此导致计算的温度结果几乎不变化。

      SUBROUTINE HETVAL(CMNAME,TEMP,TIME,DTIME,STATEV,FLUX,     1 PREDEF,DPRED)C      INCLUDE 'ABA_PARAM.INC'C      CHARACTER*80 CMNAMEC      DIMENSION TEMP(2),STATEV(20),PREDEF(*),TIME(2),FLUX(2),     1 DPRED(*)      PARAMETER(b=0.36D0,c=0.74D0,Q0=364000D0,Cemass=500)     C     qv(t)      FLUX(1)=Cemass*Q0*b*c/24*(Time(2)/24)**(C-1)*exp(-b*(Time(2)/24)**c)   c      FLUX(1)=10000000000 C     qv(t)      STATEV(1)=FLUX(1)C     Q(t)      STATEV(2)=Cemass*Q0*(1.0D0-exp(-b*(Time(2)/24)**c))C     Q0      STATEV(3)=Q0C     b      STATEV(4)=bC     c      STATEV(5)=c      RETURN      END

05

模型材料设置

这里我采用一个简单的混凝土长方体来模拟水化热温度场，其中关于水化放热的设置见图6。heat generation这里一定要点选，否则模型不会考虑HETVAL子程序。同时，user material中需要定义结果状态变量，这里我定义20个，事实上子程序中需要输出的状态变量仅为5个，但这样设置一是不影响，而是方便后续额外增加新的状态变量，利于观察调试模型。

图6 Heat Generation设置

06

水化热温度场计算结果

图7为混凝土水化热温度场在11.6小时的温度等值线图，这里的NT11表示节点温度（节点温度也就一个自由度）。可以看出，等值线图几乎一样均匀，这是由于整个混凝土块与外界保持绝热状态，每个混凝土单元都在放热，也就是每个混凝土单元之间不存在温度梯度（图8热流量密度HFL等值线图也可证明，可看出热流量密度HFL很小）。而图9表示了混凝土水化热温度场计算结果随时间的变化曲线，可以看出来混凝土一直在升温，也就没有热传导，相当于均匀升温。图10的动画用等值线图的方式表示了混凝土水化热温度场随着时间的变化，但实际上在每个时刻的温度场等值线图是一致的。

图7 NT11温度等值线图

图8 HFL热流量密度等值线图

图9 模型任一单元的温升时程图

07

水化放热速率及累计放热量计算结果

图11可发现SDV1（混凝土水化放热速率，FLUX(1)）的等值线图，发现等值线都是一致均匀的，这是因为每个混凝土单元都在采用同样的水化放热速率曲线，水化放热速率都是一致的。另外，利用XYData和XYPlots，选择图12所示的某一混凝土单元，观察SDV1随时间的变化规律。图13中也绘制了混凝土水化放速率随着时间的变化曲线。同时图14也绘制了混凝土累计水化放热量随着时间的变化曲线。因此，子程序在计算过程中的正确性得以保证。可以看出，水化放热速率由0突然增大到最大值，然后逐渐减少，在24小时后放热速率变得较小且趋于恒定，但仍旧在放热（之前已经释放掉大部分的热量）。

图11 SDV1（FLUX1(1)）等值线图

图12 利用XYData绘制SDV1（FLUX1(1)）随时间的变化曲线

图13 SDV2（混凝土累计水化放热量）等值线图

而基于选中单元的中心（与SDV1一致），绘制SDV2(混凝土累计水化放热量)随时间的变化规律如图14所示。

图14 利用XYData绘制SDV2（混凝土累计水化放热量）随时间的变化曲线

最后，有需要欢迎通过微信公众号联系我们。

微信公众号：320科技工作室。