混凝土坝瞬态热传导分析的快速边界面法

何建平，张见明，李光耀，郭帅平（.湖南大学汽车车身先进设计制造国家重点实验室，湖南长沙 40082；2.中国电建集团中南勘测设计研究院有限公司，湖南长沙 4004）



混凝土坝瞬态热传导分析的快速边界面法

何建平1，2，张见明1，李光耀1，郭帅平1
（1.湖南大学汽车车身先进设计制造国家重点实验室，湖南长沙 410082；
2.中国电建集团中南勘测设计研究院有限公司，湖南长沙 410014）

摘要：采用与快速算法相结合的边界面法，对大坝浇筑过程中的瞬态热传导问题进行仿真分析，预测混凝土大坝施工过程中温度的分布及随时间的变化。采用时间卷积法计算时间积分，通过Taylor级数展开将基本解中的空间变量和时间变量分离，一次性计算并存储空间变量的积分，提高了卷积积分的计算效率。算例仿真分析表明：边界面法的计算温度分布与有限元法计算结果高度吻合，将边界面法应用于真实大坝浇筑过程瞬态热传导分析是可行的。

关键词：边界面法；有限元法；时间卷积法；瞬态热分析；混凝土坝

计算机辅助工程（CAE）对于推动产品研发和指导工程建设具有重要意义。计算机技术的快速发展推动了CAE的发展。但CAE仍面临许多难题[1]，比如如何对复杂几何形体进行离散才能进行有效计算，如何处理大规模工程问题的数值计算。陆续涌现出的有限差分法（FDM）、有限体积法（FVM）、有限元法（FEM）和边界元法[2]（BEM）都未能很好地解决这些问题。目前大多数CAE软件采用的是有限元法。有限元法需对整个求解域进行离散处理，其应力计算精度较位移精度低。

相比有限元法，另一种有吸引力的数值方法是以边界积分方程为基础的边界元法[3]，该法不仅可以减轻有限元法的网格离散困难，还具有高精度和便于处理无限域问题、边界及界面问题、奇异性问题等突出优点，因此在工程数值分析领域被公认是有限元法的一个重要补充。近年来一些关键技术的突破使边界元法又发生了革命性的变化：①稠密矩阵稀疏化技术，包括快速多极算法[4]、分级矩阵和小波算法（其中快速多极算法是20世纪十大算法之一），可以将边界元满阵方程的求解复杂度降为线性，从而使得边界元法在解题规模上优于有限元法；②对偶互易法、径向积分法[5]和奇异边界法[6]的提出大大扩展了边界元法的应用范围，使得边界元法不仅能够求解基本解存在的线弹性问题，还可以求解非均质问题、非线性问题和时域问题；③奇异积分[7-8]、近奇异积分处理技术的发展[9]，使得边界元法可以基于三维理论求解超薄型结构。然而，在传统的边界元法中，三维CAD几何模型被离散成边界元分析模型后，CAD模型的原始几何信息基本被丢掉，导致设计和分析变成了两个相互独立的过程。由于CAD模型和分析模型的分离，在边界元自适应网格细分过程中，需要反复地与CAD系统进行交互，而每个交互过程都是很烦琐的。

为了克服上述缺点，张见明等[10]在边界元法和边界点法[11]的基础上提出了边界面法（BFM）。在边界面法中，对边界的数值积分和场变量的插值都在边界曲面的二维参数空间中进行，CAE分析直接在CAD模型上进行，实现了复杂结构的CAE分析自动化。由于CAE模型与CAD几何模型融为一体，不管网格离散有多么粗糙，分析模型在几何上都是精确的，并且在自适应网格细分过程中，不再需要与CAD系统反复地进行交互，使自适应分析变得简单。目前边界面法已经取得了许多研究成果，如参数曲面内的网格自动生成算法[12]、多域稳态热传导分析[13]、弹性静力分析、温度应力分析、声场分析[14]以及近奇异积分计算[9]。张见明等[15]通过快速多级子算法、自适应交叉拟合（ACA）和分级矩阵（H-matrix）技术，成功地将计算量级从O（N2）降至O（NlgN），并且基于GPU高性能并行计算，在CUDA编程环境中实现了边界面法正则积分的并行加速，使边界面法的大规模工程应用成为可能。

混凝土大坝施工过程的温度控制对于防止大坝混凝土开裂进而保证大坝安全运行具有重要意义。温度控制的目的是保证基础温差、上下层温差、内外温差在许可范围内，以保证温度应力在可控范围内，对大坝施工期及运行期瞬时温度场的预测是温度控制的关键。笔者运用边界面法对大坝逐层浇筑施工期及运行期的瞬态温度进行数值仿真，结合自适应交叉拟合算法和基本解时间积分与空间积分分离方法解决瞬态热传导问题中存在的计算时间长、需要存储空间大等问题，对大坝自然冷却过程中经典瞬态热传导过程[16]及真实大坝冷却过程进行仿真，以验证边界面法计算大坝瞬态热传导问题的有效性。

1 瞬态热传导问题的时域卷积法

1. 1 边界积分方程

大坝浇筑及其温度场演变过程是一个包含水化热生成的瞬态热传导过程，坝体几何形状、内部温度、热流密度随时间而变化。瞬态热传导问题的数学模型为

式中：Ω为目标结构所占区域，其边界Γ＝Su∪Sq∪SR；Su、Sq、SR分别为已知温度、热流密度、对流关系的边界；n为外法向量；u、q分别为温度和热流密度；-u、-q为温度和热流密度的边界已知量；.u为温度对时间的导数；Q为热源密度函数；ua为对流表面流体的平均温度；k、ρ、c分别为材料的传热系数、密度、比热容；β为对流系数。

瞬态热传导问题的基本解[17]为其满足式中：u*为位移基本解；K为导温系数；q*为位移基本解的法向导数；y、x分别为源点和场点；t、τ分别为时间变量和待分析时刻；r为源点和场点间的距离函数；δ为狄利克雷函数。运用式（2）中的基本解可得到式（1）的等效积分形式[17]：式中θ为y点处的立体角。式（4）中的边界积分方程为实现瞬态热传导数值仿真的基础。为简单起见，暂不考虑热源及初始条件，亦即假设结构内部无热源并且结构的初始温度为零。实际上，当热源密度分布、初始温度分布满足Laplace方程时，式（4）中相应的体积分项可以通过互易原理转化为边界积分。边界积分方程可表示为

1. 2 时域卷积法

运用边界面法分析瞬态热传导问题时，时域方法计算时间积分的方案有拟初始条件方法和时间卷积法。拟初始条件法的基本原理是每一个时间步的计算结果都作为下一步计算的初始条件，因此需要计算对初始条件的体积分；时间卷积法的基本原理是每一步的计算都需要用到初始时刻到当前时间步之间所有时间步的边界温度及热流密度。采用时间卷积法，结合边界面法离散技术，将方程（5）最终转化为线性方程组进行求解。

将时间区间[0，τ]均匀分成M个长度为Δt的时间步，在每个时间步内假设温度和热流密度不随时间改变，于是式（5）通过时间离散为

使用边界面法将上述边界积分方程进行空间离散，具体过程可以参见稳态热问题的离散过程[13]，以下只列出离散之后的边界积分方程：

方程（10）组织成矩阵形式为式中：N为单元形函数；下标i、j分别为域边界号和单元号。

式（11）在考虑边界条件后重新组织可以得到式中XM为tM时边界上的未知量。

1. 3 基本解的展开

在式（16）右端积分的计算中，对于每一个M，前面时间步对应的矩阵HmM、GmM都需要进行更新，并且分别与节点温度向量、节点热流密度向量相乘。当边界条件类型以及时间步长都未发生改变时，矩阵A是不变的，因此A只需要分解求逆一次，便可存储起来供后续时间步计算使用。但是H1M、G1M在每一个时间步计算时都需要计算并存储，计算第M个时间步时就得计算并存储2M个矩阵。特别当时间步长发生改变时，所有卷积矩阵都得重新计算，这种计算非常耗时，缺乏应用的可行性。

研究发现，这种情况的根源在于基本解中时间变量和空间变量相互掺和在一起，虽然空间变量一直没发生变化，但随着时间的更新，矩阵必须不断更新，每更新一次就得计算一次边界积分。针对此问题的解决方法是将基本解通过Taylor级数展开，将时间变量从基本解中分离出来，然后取展开后无穷级数中的有限项进行近似计算。因此，空间变量的积分可以一次性计算并存储，卷积矩阵可以通过这些存储的矩阵计算出来，而不需要在每一个时间步中反复计算边界积分，从而节省大量的计算时间。

式中的误差函数有如下展开形式：将式（20）代入式（19）可得

式中pn（r）为与时间无关的量。对式（21）中的无穷级数截断取有限项，在计算边界积分时，首先将有限个pn（r）对应的影响矩阵Pn计算并存储，在后续计算中可直接使用。对式（12）中卷积项通过基本解展开可表示为为了将时间变量和空间变量分离，简记式中NT为截断项数。H（n）ij仅与空间有关，可在分析之前一次性计算并存储，在每个时间步的卷积矩阵计算中可直接使用，不需重复计算，节约了计算时间。

（24）

2 算 例

2. 1 大坝的自然冷却

图1 几何模型

以如图1所示大坝模型为例分析其自然冷却过程。模型的上部为混凝土坝体，下部为基底。假定混凝土内部初始温度为20℃，基底初温为与外界气温均为0℃，模型周边施加固定温度边界条件，其值也为0℃，基底和坝体内均无热源。研究在自然冷却条件下，域内温度随时间的变化过程。在使用有限元法对本文中的算例仿真时，均将其简化为二维模型，而在采用边界面法仿真时，仍然将其作为三维问题进行计算。

图2至图4为部分时间点（1 a，5 a，8 a）对相应的温度分布云图。可以看出，通过边界面法和有限元计算得到的温度云图高度吻合。表1列出了两种方法在截取的部分时间点的温度最大值与最小值，数值基本吻合。

图2 1 a后坝体上温度分布（单位：℃）

图3 5 a后坝体上温度分布（单位：℃）

图4 8 a后坝体上温度分布（单位：℃）

表1 有限元与边界面法自然冷却温度计算结果

本算例证明边界面法在瞬态热传导分析中具备很高的分析精度，具备在实际工程结构中应用的可行性。

2. 2 多域问题的瞬态热传导

以文献[16]中的经典算例进行比较计算，计算模型为岩基上16层混凝土浇筑块，几何模型如图5（a）所示。浇筑块长25m，厚1m。材料热导率为2. 78 W/（m·K），表面传热系数为16. 67 W/（m2· K），混凝土初温和气温均为0℃，层间间歇10 d，混凝土的绝热温升θ（τ）＝25. 0τ/（4. 50＋τ），其中θ为热生成率，τ为龄期。基岩的绝热温升θ（τ）＝0。边界面法网格模型如图5（b）所示，边界条件包括三类，如图6所示。

图5 16层浇筑块模型

图6 边界条件示意图

按照如前所述的热源和边界条件进行16层浇筑块的瞬态热分析，并将计算结果和有限元法计算结果进行对比。为了得到精确的对比数据，截取如图5中所示的A-A截面上的温度值绘成温度曲线如图7所示。边界面法在分析多域瞬态热传导问题的计算结果与文献[16]中有限元计算结果吻合较好，证明边界面法用于多域瞬态热分析计算是可行的。

翠山公园在斜坡及陡坎下部出露地层为Q3-4apl砂卵砾石，灰白-青灰色，母岩以花岗岩、凝灰岩、长石石英为主，磨圆度中等，分选性差，含极少量砂质，采样分析其d10均大于0.1 mm。斜坡上部覆盖Q3-4apl黄土状粉土，分选性好，垂直节理发育，根据现场采样分析其d60在0.02～0.04 mm。上下层间岩性差异使得其渗透系数存在较大差异，同时上部粉土层垂向节理的发育，使得粉土层地表水下渗主要以捷径式渗透为主。地表水捷径式下渗是翠山南部地表坑穴形成的原因之一。

2. 3 真实大坝的瞬态热传导

某大坝逐层浇筑过程的瞬态热传导问题模型如图8所示，上部为混凝土坝体，下部为基底。分析中考虑坝体逐层浇筑过程以及浇筑后的热传导状态，时间跨度为9 a，共1260个时间步，材料的导热系数为2. 78 W/（m·K），热扩散率为1. 11×10-6m2/ s，环境温度为0℃，绝热温升

图9为在坝体浇筑过程中的4个时刻（100 h、400 h、38 d、9 a）的瞬时温度分析结果。该算例坝体温度随浇筑过程变化的动态图在“完整实体工程分析软件中心”（http：/ / www. 5acae. com/）有具体呈现。对比表明边界面法可对逐层浇筑的大坝进行实时热传导分析，可用于混凝土大坝实际工况瞬态温度仿真。

图7 16层浇筑块温度分布对比

图8 某大坝模型（单位：m）

图9 大坝浇筑不同时刻的温度分布（单位：℃）

3 结 语

应用边界面法对大坝混凝土的浇筑过程进行瞬态温度场计算，实现了BFM与CAD的无缝连接，所有材料特性、边界条件、初始状态、热源，甚至网格自动划分等完全内置于CAD环境里，实现了大坝瞬态热传导分析的自动化。在边界面法中，由于CAE模型与CAD模型的一体化，避免了分析模型的几何误差，从而可以有效提高数值计算的精度。在运用时间卷积法计算积分时，将基本解按Taylor级数展开，

分离时间变量与空间变量，并将空间变量的积分进行一次性计算并存储，大大减少了计算量。通过对大坝热传导问题中经典算例的分析，并将分析结果与FEM结果进行对比，验证了该方法结果的可信性，且自动化的分析环境显著减少了计算时间和人的劳动量，表明该方法在大型工程应用上具有可行性。

参考文献：

[ 1 ] ZHANG Jianming，TANAKA M，ENDO M. The hybrid boundary node method accelerated by fast multipole method for 3D potential problems [ J ]. International Journal for Numerical Methods in Engineering，2005，63：660-680.

[ 3 ] ZHANG Jianming，YAO Zhenhan，LI Hong. A hybrid boundary node method [ J ]. International Journal for Numerical Methods in Engineering，2002，53：751-763.

[ 4 ] LIU Yijun，CHEN Shaohai. A new form of the hypersingular boundary integral equation for 3-D acoustics and its implementation with C0boundary element [ J ]. Computer Methods in Applied Mechanics and Engineering，1999，173：375-386.

[ 5 ] GAO Xiaowei. The radial integration method for evaluation of domain integrals [ J ]. Engineering Analysis with Boundary Elements，2002，26：905-916.

[ 6 ] CHEN Wen，WANG Fuzhang. A method of fundamental solutions without fictitious boundary [ J]. Engineering Analysis with Boundary Elements，2010，34（5）：530-532. [ 7 ] GUIGGIANI M，GIGANTE N. A general algorithm for the numerical solution of hypersingular boundary integral equations[J]. Journal of Applied Mechanics，1992，59：604-614.

[ 8 ] GAO Xiaowei. An effective method for numerical evaluation of general 2D and 3D high order singular boundary integrals [ J]. Computer Methods in Applied Mechanics and Engineering，2010，199：2856-2864.

[ 9 ] XIE Guizhong，ZHANG Jianming，QIN Xianyun，et al. New variable transformations for evaluating nearly singular integrals in 2D boundary element method[J]. Engineering Analysis with Boundary Elements，2011，35：811-817.

[10] ZHANG Jianming，QIN Xianyun，HAN Xu，et al. A boundary face method for potential problems in three dimensions [ J ]. International Journal for Numerical Methods in Engineering，2009，80：320-337.

[11]宋敏，覃先云，张见明，等.快速杂交边界点法在水坝热传导分析中的应用[J].水利水电科技进展，2010，30 （1）：27-29.（SONG Min，QIN Xianyun，ZHANG Jianming，et al. Application of fast hybrid boundary node method to heat conduction in concrete dams [ J ]. Advances in Science and Technology of Water Resources，2010，30（1）：27-29.（in Chinese））

[12] QIN Xianyun，ZHANG Jianming，LI Guanyao，et al. A finite element implementation of the boundary face method for potential problems in three dimensions [ J ]. Engineering Analysis with Boundary Elements，2010，34：934-943.

[13]张见明，李湘贺，陆陈俊，等.多域边界面法在稳态热传导问题中的应用[J].湖南大学学报（自然科学版），2014，41（5）：58-64.（ZHANG Jianming，LI Xianghe，LU Chenjun，et al. Application of multi-domain boundary face method in steady state heat conduction [ J]. Journal of Hunan University（Natural Sciences），2014，41（5）：58-64.（in Chinese））

[14] WANG Xianhui，LU Chenjun，ZHOU Fengling，et al. Acoustic problems analysis of 3D solid with small holes by fast multipole boundary face method [ J]. Engineering Analysis with Boundary Elements，2013，37：1703-1711.

[15]张见明，余列祥，刘路平.基于GPU加速的边界面法正则积分的研究[J].湖南大学学报（自然科学版），2013，40（3）：41-45.（ZHANG Jianming，YU Liexiang，LIU Luping. Research on regular integration in boundary face method based on GPU-acceleration [J]. Journal of Hunan University（Natural Sciences），2013，40（3）：41-45.（in Chinese））

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

[17] BREBBIA C A，TELLES J C F，WROBEL L C. Boundary element techniques：theory and applications in engineering [M]. Berlin：Springer，1984.

Fast boundary face method for transient heat conduction analysis of concrete dams/ /

HE Jianping1，2，ZHANG Jianming1，LI Guangyao1，GUO Shuaiping1（1. State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body，Hunan University，Changsha 410082，China；2. Power China Zhongnan Engineering Corporation Limited，Changsha 410014，China）

Abstract：The boundary face method（BFM），combined with a fast algorithm，was applied to transient heat conduction analysis of concrete dams during the pouring process，so as to predict the temperature distribution and its time history in concrete dams during the construction process. Using the time convolution method to compute the time integration and the Taylor series expansion method to separate the spatial variables and time variables in fundamental solutions，the spatial variable integrals can be calculated and stored at the same time，improving the calculation efficiency of the convolution integral. A case study shows that the temperature distribution obtained by the fast BFM agrees with the results of the finite element method，demonstrating the feasibility of application of the fast BFM to transient heat conduction analysis during the actual dam construction process.

Key words：boundary face method；finite element method；time convolution method；transient heat conduction analysis；concrete dam

收稿日期：（2014 09 25 编辑：郑孝宇）

作者简介：何建平（1966—），男，高级工程师，博士，主要从事水电、机械计算机仿真研究。E-mail：jianping_he@163. com

基金项目：国家自然科学基金（11172098）；国家重点基础研究项目（973计划）（2010CB328005）

中图分类号：O241；O414

文献标志码：A

文章编号：1006 7647（2016）01 0071 06