
4.1 温度场数值计算原理
大坝面板混凝土通常承受两种不同性质的荷载:第一类包括水压、泥沙压、地震、渗压、风浪、冰凌以及结构自重与设备重量等;另一类是混凝土本身的体积变化所引发的荷载,包括温变、徐变、干湿变化、混凝土自生体积变形等所引发的荷载。对于前一类荷载,要保证结构不产生或只产生很小的拉应力并不困难。但在施工和运行期间,要把后一类荷载所产生的拉应力控制在允许范围内则是一件很不容易的事情。正是由于后一类荷载(其中主要是温变)的作用,在大坝混凝土结构中会由于产生过大的拉应力而出现裂缝。
温度变化对大体积混凝土的影响主要有:引起结构内力的变化,导致混凝土裂缝;对结构的应力状态引起应力重分布,不能按照设计时确定的应力状态发展。温度变化引起的应力甚至超过其他荷载应力,尤其是在结构温度急剧变化时,将产生很大的拉应力,而混凝土为脆性材料,抗拉强度只有抗压强度的1/10左右;拉伸变形能力很小,短期加载时的极限拉伸变形只有(0.6~1.0)×10-4,约相当于温度降低6~10℃的变形;长期加载时的极限拉伸应变只有(1.2~2.0)×10-4。大坝的面板混凝土通常是暴露在外面的,表面与空气或水接触,一年四季中气温和水温的变化在大体积混凝土结构中会引起相当大的拉应力,而导致混凝土的开裂破坏。
因此通过数值模拟计算对大坝面板混凝土进行抗裂分析,依据面板混凝土的浇筑进度计划,计算不同时段不同工况下面板混凝土的温度应力场与温度应力值以及变形分析。
4.1.1 方程的有限元格式
在有限元计算中,有限元格式的好坏直接决定了与被模拟物理过程的相似程度和同质性。在有限元解中,有3种误差来源:①由于域的近似性产生的误差;②由于解的近似方法产生的误差;③由于数值计算产生的误差(即在计算机中数值的积分和舍入误差)。其中前两种误差都在建立控制方程的有限元格式阶段产生,建立适合控制方程和研究域的有限元格式,能有效地减小有限元解与理论解之间的差值,是保证有限元解法收敛的必要条件。
根据热传导理论可知,混凝土的瞬态温度场应满足以下方程:

在第一类边界上:

在第二类边界上:

在第三类边界上:

在边界条件方程中φ可以代表T,α是放热系数,φa是环境温度。
将有限元模型的自由度确定为φ=[T1,T2,…,Tn,w1,w2,…,wn]T,是节点温度和节点湿度列阵。
对空间域进行离散,将空间域离散为有限个单元体,在典型单元内温度场φ仍可以近似地用节点温度φi插值得到。此时节点温度是时间的函数

对时间域进行离散,取一个典型地时间单元长度Δt,采用两点公式,单元内φ由节点值φn和φn+1插值得到

4.1.2 有限元格式中权函数的选择
在有限元计算中,常用的权函数有以下几种:配点法、子域法、最小二乘法、力矩法、伽辽金法。配点法简单地强迫余量在域内n个点上等于零,误差过大,不采用。子域法强迫余量在n个子域内的积分等于零,在非子域区域误差不可控,划分积分为零区域和积分不为零区域在一般有限元计算中无法体现作用。力矩法强迫余量的各次矩等于零,但得到的求解方程的系数矩阵往往是不对称的。从计算精度和对算法的熟悉程度来考虑权函数的取法从最小二乘法和伽辽金法当中选择。
最小二乘法是用域内误差的平方作为变分的泛函。当近似函数仍采用的线性插值时,最小二乘法的变分泛函是

假定φn已知,对φn+1变分使泛函Π为极小,经过变换后可以得到两点循环公式是


由上式可见计算工作量十分庞大。最小二乘法一般具有较高的精度,但也常常花费较多的时间。伽辽金法的精度较高,而且计算花费时间相对较少,所以本文采用伽辽金法选择权函数。
采用伽辽金法选择权函数,同时考虑到边界条件,可以建立如下的有限元格式

其中

式中 ω——时段单元内时间域的形函数。
当φn已知时,则可通过下式求得φn+1

4.1.3 有限元计算中解的稳定性
抛物线形方程用有限元法和有限差分法求解时,都存在一个稳定性的问题。这是步进积分带来的特点。在抛物型方程求解的过程中,不同的格式对稳定性会产生不同的反应,见图4.1。根据局部傅里叶准则:

图4.1 加权变系数对稳定性的影响

式中 α——物体的导温系数,m2/s;
Δt——时间步长,s;
Δx——空间网格步长,m。
曲线a将区域划分成稳定和不稳定区域两部分,曲线下方区域(1)是不稳定的,而上方是稳定的。曲线b将区域划分成振荡和不振荡两部分。曲线下方是振荡的,上方是不振荡的。由此可见,向后差分格式即θ=1才是无条件不振荡。
因为导温系数与质扩散系数相差3个数量级,因此如果采用同一套网格,则计算温度场时的傅里叶数与计算湿度场时的傅里叶数将相差3个数量级。因此在计算湿度场时要采用比计算温度场时尺寸小3个数量级的网格。
4.1.4 温度应力数值计算
根据面板混凝土材料的本构关系,第n个计算时段的位移{Δδ}满足以下有限单元法方程

式中 [K]——刚度矩阵;
——由徐变引起的荷载矩阵增量;
——由温度变化引起的荷载增量;
{ΔF}——外荷载增量。