混凝土高坝全过程仿真分析

关键词：混凝土坝，仿真计算，复合并层单元，均质并层单元，分区异步长算法，水管冷却等效热传导方程。

（2）时间步长小。早龄期混凝土，由于弹性模量、徐变度、绝热温升都随着龄期而急剧变化，需采用较短的时间步长，以保证必需的计算精度。例如取时间步长D ?=0.5天，一年就有730步，如工期为3年，共有2190步。如果要考虑日照及日气温变化的影响，时间步长应减小到1～2小时左右，工期3年就有1.3～2.6万步。每一步都需求解两次（温度场和应力场）特大型线性方程组。

（3）冷却水管计算的困难，冷却水管的半径只有1～2cm，在水管周围必需采取密集计算网格，单元尺寸必需是cm级的，虽可向外逐步放大单元尺寸，但单元总数还是很多的，例如，水管间距1.5m×1.5m，采用线性单元，在1.5×1.5m²范围内就需要约90个结点，如浇筑层尺寸为18m×1.5m×30m，在长度方向取间距3.0m，一个浇筑层就有约1.2万个结点，设坝高300m，一个坝块就有约240万个结点，就目前计算机硬件水平来说，实际上无法计算。

经过几年的努力，笔者提出了一系列新的解法，较好地克服了以上各项困难，并已大量应用于实际工程。

有限元网格划分的基本原则是，单元尺寸的大小必须与温度梯度和应力梯度的变化相适应。例如采用线性单元，单元内温度和应力都是线性变化的，在划分计算网格时必须充分考虑这一特性。

如图1所示，当坝体逐步升高时，将坝体从上到下划分为4个区域。在上面新浇筑的区域R₁中，沿浇筑层厚度方向的温度变化和应力变化都比较快，每个浇筑层都应划分为n层有限元，例如n=4～5层。在其下面的区域R₂中，每个浇筑层在厚度方向的温度增量和应变增量已是线性分布的，可把原来每个浇筑层中的n层单元合并为一层单元，即一个浇筑层为一层有限单元。在再下面的区域R₃中，沿铅直方向的温度增量和应变增量的变化进一步平缓，可把若干个浇筑层合并为一个复合并层单元，复合单元内各浇筑层保留各自的力学特性和热学特性，如图2。再下面的区域R₄，已到晚期，可把几个浇筑层合并为均质并层单元，此时单元内各浇筑层的力学特性和热学特性已充分接近，可取其平均值，如图3。

并层后，单元内各浇筑层仍保持原来的弹性模晨E(?)、徐变度C(t,??)和绝热升温q (t )，按下列各式计算第n时段的单元刚度矩阵[K_n]、徐变引起的结点荷载增量和应力增量^[2,3,4]：

（1）

（2）

（3）

式中[B]为几何矩阵，为等效弹性矩阵，为与徐变有关的向量，为温度应变增量，为自生体积应变增量。由于单元内各浇筑层具有不同的力学特性，所以要进行分区数值积分。

2.2 均质并层单元

在晚期，当t _i ≤t ≤t _j各层合并后，由于单元内各层力学特性和热学特性已充分接近，可采用平均弹性模量和平均徐变度，得到一均质单元，如图3^[5]。

并层以后，单元数量大大减少。例如，对于常态混凝土坝，浇筑层厚1.5m，原先每个浇筑层内分为5层有限单元，并层以后，把4个浇筑层合并为一层有限单元，一个并层单元就包含了原来的20层有限单元。又如，对碾压混凝土坝，浇筑层厚0.30m，分为4层有限单元，并层后单元高度为6.0m，一个并层单元包含了原来的80层有限单元，因此，计算得到极大的简化。

对于早龄期混凝土，由于弹性模量和徐变度随着龄期而急剧变化，需采用很短的时间步长（D t =0.2～1.0）天，如果希望考虑日照和昼夜温差的影响，则应取D t =1～2小时，以保证必要的计算精度。对于晚龄期混凝土，因弹性模量和徐变度的变化时率已经很平缓，可以采用较大时间步长（D t =10～30d）。但在施工时期，坝块上部不断浇筑新混凝土，因此在整个施工时期，都要采用较小的时间步长。我们提出分区异步长算法，对早龄期混凝土，采用小时间步长；而对老混凝土采用大时间步长，可节省大量机时。

如图4所示，坝块划分为a和b两个区域，cc是分界面。从时间t_n至t_m，区域a采用一个大时间步长：D t _nm=t_m-t_n；而在区域b采用m个小时间步长：D t _n+1，D t _n+2，…，D t _n+m，并令

（4）

以便当t =t_m时上下层在时间上可互相衔接。

弹性徐变应力计算步骤如下：

第1步：固定接触面cc，以小时间步长，利用弹性徐变理论计算区域b中的应力增量及接触面cc上的反力。到时间t_m，区域b中的应力增量和接触面cc上的反力分别为：

（5）

第2步：固定接触面cc，以一个大时间步长，利用弹性徐变理论计算区域a中的应力增量和接触面cc上的反力。

第3步：释放接触面上的反力，得到应力增量和。

综合以上三步，最后得到应力增量如下：

（6）

把计算域划分为3个子域，如图5所示，R₁是温度变化剧烈区，R₂是过渡区，R₃是温度变化平缓区，W ₁和W ₂是分界面，B₁、B₂、B₃是相应的表面。

周围介质温度也相应地划分为三种：在边界B₁上介质温度为T_c1+f (t)，其中T_c1为变化平缓的温度（如年变化）或常数，f (t)为急剧变化的温度，如日变化或寒潮；在边界B₂上介质温度为T_c2；在边界B₃上介质温度为T_c3。

以混凝土坝为例，R₁是新浇筑的混凝土，需要考虑预冷骨料产生的初始温差、急剧变化的绝热温升和寒潮及气温日变化。R₂是较老的混凝土，R₃是老混凝土，在R₂+R₃中，绝热温升的变化和初始温度的分布已很平缓，寒潮和气温日变化已不必考虑，只需考虑介质温度的年变化。

求解的问题如下：

（7）

（8）

（9）

式中：a为导温系数；

_l 为导热系数；

b ₁、b ₂、b ₃为表面放热系数；

q ₁、q ₂、q ₃为绝热温升；

n 为法线；

t 为时间；

T为温度。

由于问题是线性的，可进行分解，令

T = U +V （10）

在三个区域中，绝热温升q 、初始温度及边界条件的分解见表1，表中q ₁、q ₂、q ₃、T₁(0)、T₂(0)、T₃(0)、T_c1、T_c2、T_c3、f(t)等都是坐标（x, y, z, t）的已知函数。由于问题是线性的，上述分解是严格的。引起温度急剧变化的因素如q ₁（早期水化热）、f(t)（气温急剧变化）、T₁(0)（与周围温度不同的初始温度）都放在子域R₁中，因此，温度的急剧变化局限在R₁中，可能波及过渡区R₂，到了R₂与R₃的接触面W ₂上，其变化已趋于零。即t _< t_s（t_s可由经验估计或由试算决定）时，W ₂上U=0。因此U满足下列条件：

（11）

式（11）的解为：

在R₃中： U = 0 （12）

表1 问题的分解

因此，U场只需在R₁+R₂中求解，由于U变化剧烈，要采用小时间步长。V场则需在全域R₁+R₂+R₃求解，但因V变化平缓，可采用大时间步长。通常R₁+R₂远小于R₃，因而可节省大量计算时间。

笔者在文献[8,9]中提出了用有限元计算冷却水管的方法，由于冷却水管的半径只有1cm左右，在水管附近必须采用非常密集的网格，因而很难用于混凝土高坝的仿真计算，为了克服这个困难，笔者提出了如下的等效热传导方程^[10]：

（13）

式中 （14）

（15）

（16）

其中D为冷却圆柱体的直径，L为冷却水管长度，c_w、r _w、q_w为冷却水的比热、容重和流量，，其中s₁、s₂分别是水管在水平和铅直方向的间距。

（17）

（18）

如混凝土绝热温升用任意函数f(t)表示如下：

（19）

（20）

在这里，用函数f (t)和y (t)考虑冷却水管的作用，包括水管半径、间距、长度、流量等因素的影响，在划分有限元网格时就不必采用密集的网格了。

基于上述算法，我们编制了计算程序Simu Dam，可在微机上进行混凝土高坝的仿真计算，先后为三峡、龙滩、东风、小湾、沙牌等一系列重力坝和拱坝进行了仿真计算，取得了良好的计算效果，图6～8为一个碾压混凝土拱坝的计算结果。

图6 某碾压混凝土拱坝下游立视图 图7 拱坝剖面C第352天温度场

图8 拱坝上游面第450天（已竣工）主应力

采用我们提出的这一套新的计算方法，对于新混凝土，采用密集的计算网格和很小的时间步长，因而大大提高了计算精度。对于老混凝土，采用并层单元和较大时间步长，大大提高了计算效率，因此，我们提出的这一套计算方法，既提高了计算精度，又提高了计算效率，为混凝土高坝仿真计算开辟了一条崭新的道路。

1. 朱伯芳，宋敬廷，混凝土温度场及温度徐变应力的有限元分析，见“水利水电工程应用电子计算机资料选编”，北京：水利电力出版社，1977。
2. 朱伯芳，大体积混凝土温度应力与温度控制[M]，北京：中国电力出版社，1999年3月，p.414-415。
3. 朱伯芳，许平，混凝土坝仿真应力分析[J]，混凝土坝技术，1999,2，p.11-17
4. Bofang Zhu, PingXu, Methods for stress analysis simulating the construction process of high concrete dams, Dam Engineering, Vol. XI, Issue 4, 243-260.
5. 朱伯芳，混凝土坝仿真应力分析的并层算法[J]，水力发电学报，1994.3，p.21-29。
6. 朱伯芳，弹性徐变体有限元时间域分区异步长解法，水利学报，1995年第7期，p.24-27。
7. 朱伯芳，不稳定温度场时间域分区异步长解法，水利学报，1995年第8期，p46-52。
8. 朱伯芳等，水工混凝土结构的温度应力与温度控制[M]，北京：水利电力出版社，1976。
9. 朱伯芳，蔡建波，混凝土坝水管冷却效果的有限元分析[J]，水利学报，1985,4，p.27-36。
10. 朱伯芳，考虑水管冷却效果的混凝土等效热传导方程[J]，水利学报，1991.3，p.28-34。

Zhu Bofang and Xu Ping

(China Institute of Water Resources and Hydropower Research, Beijing, 100038)

Abstract:

The stresses in the high concrete dams are influenced remarkablely by the process of construction. Methods are proposed to compute the stresses in the high concrete dams simulating the construction process.