（水工结构工程专业论文）碾压混凝土坝温度应力仿真计算的非均质层合单元法及其应用研究.pdf

摘要 本文围绕碾压混凝土重力坝的温控防裂问题，对碾压混凝土非稳定温度场和应力场 的仿真计算进行了系统地学习和研究，并对马堵山碾压混凝土重力坝溢流坝段的温控问 题进行了详细地计算和分析。主要内容如下： 1 仿真计算高碾压混凝土坝温度场和应力场时，存在计算量大的问题，为了减少 计算工作，本文运用“非均质层合单元”仿真计算模型。运用此方法在数学上更为严密， 为解决高碾压混凝土坝温度应力施工仿真计算的计算量大的问题提供了很好的参考方 案。 2 根据试验测得的混凝土试件的温度值，运用遗传算法反演分析多个混凝土热学 参数。运用反演方法获得实际工程所需的计算参数，简单易行。 3 基于“非均质层合单元法”的温度场和应力场的仿真计算模型，编制了碾压混 凝土坝的三维有限元仿真计算程序。进行了立方体混凝土试块的温度理论解和有限元解 的比较，再次证明了有限元解的有效性。 4 在前人的研究基础上，本文阐述研究了诱导缝的开裂机理，裂缝模型，分析了 诱导缝应力强度因子等方面的内容。 5 以马堵山碾压混凝土坝项目为背景，深入系统地进行了典型溢流坝段的温控防 裂研究。 关键词：碾压混凝土重力坝，非均质层合单元法，仿真计算，混凝土热学参数，反 演分析，温度场，应力场，诱导缝 a b s t r a c t f o c u s e do nt h es u b j e c t so ft e m p e r a t u r ec o n t r o la n dc r a c kp r e v e n t i o n ，t h i sp a p e ri s m a i n l ya b o u tt h es t u d yo f t h es i m u l a t i v ec a l c u l a t i o no ft h eu n s t e a d yt e m p e r a t u r ea n dt h e r m a l s t r e s sf i e l di nt h er c cg r a v i t yd a m s t h er e s e a r c ho ft e m p e r a t u r ec o n t r o lo ft h em a d u s h a n r c c g r a v i t yd a m si sa l s op r e s e n t e di nt h i sp a p e r t h em a i nc o n t e n t sa r ea sf o l l o w s ： 1 i no r d e rt or e d u c et h ew o r ko fc a l c u l a t i o n ，t h em o d e lo ft h en o n h o m o g e n e o u s m u l t i 1 a m i n a t ee l e m e n ti sa p p l i e dh e r et os i m u l a t et h et e m p e r a t u r ea n ds t r e s sf i e l do ft h eh i g h r c cd a m s t h et h e o r yo ft h i sm e t h o di ss t r i c ti nm a t h e m a t i c sa n dc a nc u td o w nt h el a r g e w o r ko f t h es i m u l a t i v ec a l c u l a t i o n 2 a c c o r d i n gt ot h ev a l u e so ft e m p e r a t u r ea t t a i n e db yt h et e s t ，w eu s et h em e t h o d o fg a t or e v e r s es e v e r a lp a r a m e t e r so fc o n c r e t e i ti sc o n v e n i e n tt og e tt h ep a r a m e t e r sn e e d e di nt h e p r o j e c tb yu s i n gt h eb a c ka n a l y s i sm e t h o d 3 b a s e do nt h es i m u l a t i o nm o d e lo ft h en o n - h o m o g e n e o u sm u l t i - l a m i n a t ee l e m e n t ，a l a r g e3 - dp r o g r a mf o rh i 曲r c c dh a sb e e nd e v e l o p e d t h er e s u l to f t h ec o m p a r i s o nb e t w e e n t h ev a l u e so ft e m p e r a t u r eg a i n e db yc l a s s i c a lm e t h o da n df e mr e s p e c t i v e l ys h o w st h a tt h e l a t t e rm e t h o di sv a l i d 4 b a s e do nt h ep r e v i o u sr e s e a r c h ，t h et h e o r i e so ff a c t u r ea n ds t r e n g t ho ft h ei n d u c e d c r a c k si nr c cd a m sa r ei n t r o d u c e da n ds t u d i e dh e r ea sw e l la st h em o d e la n dt h es t r e n g t h f a c t o r so f t h ei n d u c e dc r a c k s 5 b a s e do nt h em a d u s h a np r o j e c t ，t h es u b j e c to ft e m p e r a t u r ec o n t r o la n dp r e v e n t i o n f r o mt h ec r a c k i n go fat y p i c a lo v e r f l o ws e c t i o nb l o c ko ft h er c cg r a v i t yd a mi ss t u d i e da sa c a s es t u d y k e y w o r d s ：r c cg r a v i t yd a m ，i n h o m o g e n e o u sm u l t i l a n f i n a t ee l e m e n t ，s i m u l a t i v ec a l c u l a t i o n ， c o n c r e t et h e r m a lp a r a m e t e r , b a c ka n a l y s i s ，t e m p e r a t u r ef i e l d ，s t r e s sf i e l d ，a n d i n d u c e d j o i n t 学位论文独创性声明： 本人所呈交的学位论文是我个人在导师指导下进行的研究工作及取得 的研究成果，尽我所知，除了文中特别加以标注和致谢的地方外，论文中 不包含其他人已经发表或撰写过的研究成果。与我一同工作的同事对本研 究所做的任何贡献均已在论文中作了明确的说明并表示了谢意。如不实， 本人负全部责任。 论文作者( 签名) ：垂墨鱼形年占月日 学位论文使用授权说明： 河海大学、中国科学技术信息研究所、国家图书馆、中国学术期刊( 光 盘版) 电子杂志社有权保留本人所送交学位论文的复印件或电子文档，可 以采用影印、缩印或其他复制手段保存论文。本人电子文档的内容和纸质 论文的内容相一致，除在保密期内的保密论文外。允许论文被查阅和借阅。 论文全部或部分内容的公布( 包括刊登) 授权河海大学研究生院办理。 论文作者( 签名) ： 五九 舢6 年月日 第一章绪论 第一章绪论 1 1 碾压混凝土坝的发展及其特点 1 1 1 碾压混凝土坝发展简述“朝 碾压混凝土筑坝的设计思想原创于在允许的条件下，采用土石坝的施工方法进行干 硬性混凝土的运输、摊铺、碾压，以达到快速施工的目的，这个设想早在6 0 年代初期 已被提出来，经过四十年左右的试验和研究，碾压混凝土坝以其独特的优越性成为世界 坝工建设领域中最受喜爱的新坝型之一。 日本是摄早研究使用混凝土筑坝技术的国家之一，世界上第一座用振动碾碾压坝体 混凝土的是1 9 7 8 年9 月的日本岛地川( s h i m a j i g a w a ) 坝。该坝为高8 9 m 的碾压混凝土重 力坝，体积3 1 7 0 0 立方米，于1 9 8 0 年建成，开创了被称为”r c d 工法”的日本式碾压混凝 土筑坝技术。8 0 年代末至9 0 年代初，日本修建了一批百米高的碾压式混凝土坝，如高 1 0 0 m 的玉川坝，高1 1 5 m 的境川i 坝，高8 4 m 的朝日小川坝，高1 0 4 m 的蛇尾川坝，日本式 碾压混凝土坝的筑坝技术日渐成熟，并形成了比较规范化的设计准则和全套施工工艺。 9 0 年代以来，高1 5 5 m 的宫濑坝和高1 5 6 m 的浦山坝的顺利建设，标志着日本碾压式混凝 土筑坝技术的最高水平。 我国自1 9 8 6 年成功地建成第一座碾压混凝土坝以来，据估计，已建、在建和正在 设计的碾压混凝土坝大约有5 0 座左右，是我国坝工发展有前景的坝型之一。近期已建、 在建和即将开工建设的高l o o m 以上的碾压混凝土坝有龙滩( 2 1 6 m ) 、江垭( 1 3 1 m ) 、 百色( 1 2 6 m ) 、大朝山( 1 2 1 m ) 、棉花滩( 1 l l m ) ，其中碾压混凝土量均超过整个 大坝混凝土量的6o 以上。正在设计的红水河上的龙滩水电站碾压混凝土重力坝，混 凝土量5 3 5 万肌3 ，其中碾压混凝土量占6 5 左右，施工月高峰浇筑强度超过2 5 万珊3 ，将 创造世界先进水平。我国的碾压混凝土筑坝技术，是在8 0 年代吸收国外先进技术的基础 上发展起来的，但也创立了有自己独特的经验，概括有：坝体体型结构尽量简单；碾压混 凝土配合比采用低水泥用量：尽量实现高掺粉煤灰和掺复合外加剂；坝体以碾压混凝土 自身防渗为主，并采用富胶凝材料，提高抗渗性能；施工上采用大仓面薄层连续浇筑的工 艺等。简言之，我国的碾压混凝土筑坝技术是高掺粉煤灰，低稠度、薄层、全断面、快速 短间歇连续填筑。在坝体体型简单化方面，主要是减少孔和洞、简化坝体结构，有利碾压、 第1 页 河海大学硕士学位论文 便于施工。龙滩、江垭、百色、大朝山、棉花滩等水电站大坝都是碾压混凝土重力坝， 坝体建筑物没有引水洞管和进水口，简化了坝体结构，可缩短工期提前发电。这些成就 标志着我国的碾压混凝土坝建设已经达到国际领先水平。 其他如美国、澳大利亚、法国、西班牙、巴西、摩洛哥、南非等国也相继修建了碾 压混凝土坝。迄今为止，全世界完建和在建的坝高超过1 5 m 的碾压混凝土坝己超过2 1 0 座，它们分布在5 大洲的2 8 个国家中，其中中国4 0 座，美国2 9 座，西班牙2 1 座，巴 西1 4 座和澳大利亚9 座，约占2 1 0 座的8 0 ，这七个国家在碾压混凝土坝工建设方面 居于世界领先地位。现在，碾压混凝土筑坝已成为大坝建筑的主要技术之一。 1 1 2 碾压混凝土坝的主要特点 碾压混凝土坝的迅速发展是与其优越的技术、经济特点紧密相关的。与常规混凝土 相比，碾压混凝土坝的主要特点也是它的优点“： ( 1 ) 工期短碾压混凝土坝采取汽车运料入仓，薄层机械碾压，连续浇筑上升的快 速施工方式，使混凝土浇筑工期可缩短3 0 5 0 。 ( 2 ) 造价低由于碾压混凝土采用干贫配比，水泥用量减少了3 ( ) 6 0 ：同时由于 施工简化，节省了大量劳力、木材、钢材和混凝土冷却费用，碾压混凝土坝可使总投资 降低1 5 3 0 。 ( 3 ) 温控措施相对简化由于混凝土中水泥用量减少，同时更以低热的粉煤灰代替 部分水泥，明显降低了混凝土放热量，可使混凝土施工中复杂的温控措施得到很大的简 化。 1 2 碾压混凝土坝温度场和应力场的研究方法及进展 1 2 1 碾压混凝土坝温度场和应力场的研究方法 碾压混凝土坝温度场和应力场的研究方法一般可分为理论解析法和近似方法两大 类： 理论解析方法主要是用来求解边界条件比较简单的低维问题，常用的有分离变量法 和拉普拉斯变换法。理论解析法对了解温度场和应力场的一般规律具有定的使用价 值，但是该方法通常只能解决比较简单的问题，对于实际工程没有太大的使用价值。但 是解析方法可以验证近似方法( 尤其是数值解法) 的可靠性和精确性。 近似方法主要包括数值解法、图解法、电热模拟法和水热模拟法等。其中后面几种 第2 面 第一章绪论 方法近年来用的很少，因此近似方法主要指的是数值解法。已经发展的数值分析方法可 以分为两大类。一类以有限差分法为代表，其特点是直接求解基本方程和相应定解条件 的近似解。另一类数值分析方法是首先建立和原问题基本方程及相应定解相等效的积分 提法，然后据之建立近似解法，其主要代表是有限元法，该方法的出现是数值分析方法 研究领域重大突破性的进展。实践经验表明，一维温度场的计算，采用差分法是方便的。 二维和三维温度场的计算，则以采用有限元方法为佳。与差分法相比，有限单元法具有 的优点是易于适应不规则边界，在温度梯度大的地方可以局部加密网格，容易与计算应 力的有限单元法程序配套，将温度场、应力场和徐变变形三者纳入一个统一的程序进行 计算。采用有限元计算非稳定温度场及徐变应力场时，一般采用有限元差分法，即 在空间域采用有限元离散，在时间域采用差分离散。目前在进行坝体温度场和应力场的 仿真计算中均采用有限元法，能有效地模拟坝体的整个施工过程和材料特性参数随龄期 的变化过程以及边界条件的变化过程，实现坝体在整个施工期和运行期的仿真计算，是 目前最为常用的计算方法。 1 2 2 碾压混凝土坝温度应力研究进展畸1 在碾压混凝土问世初期，由于其水泥用量大大减少，人们曾一度认为碾压混凝土坝 不存在温度应力问题，国内外对碾压混凝土坝的研究主要集中在碾压混凝土材料、力学 性能、施工工艺和质量控制等方面，对温度应力和温度控制问题关注不多。直到t a t r o s t e p h e nb 和s c h r a d e re k 发表了“t h e r m a lc o n s i d e r a t i o n sf o rr o l l e rc o m p a c t e d c o n c r e t e ”后，大仓面连续施工的坝体温度应力和温度控制问题才得到工程界的广泛重 视，尤其是对不设横缝或者少设横缝的碾压混凝土拱坝。已经完建的几座高碾压混凝土 坝，均出现了不同程度的裂缝，这使大坝工程师开始从一个更高的起点上对碾压混凝土 坝的温度应力问题进行研究。 在国外，主要以美国和日本的工作为代表。美国工程师在世界第一座碾压混凝土坝 ( 柳溪坝) 的温度应力研究工作中，尽管采用的是比较简单的一维模型，但他们的方法 是当时最先进的，是他们第一次把结构计算与旋工过程结合起来，用逐步递推的方法， 求出了不同时期坝体的温度场和温度应力场。1 9 9 2 年，在美国加利福尼亚州圣地亚哥市 第三次碾压混凝土会议上，巴瑞特等介绍了三维温度应力计算软件a n a c a p ，其创造性在 于把b a z a n t 的s m e a r e dc r a c k 开裂模型引入大坝温度应力的分析中，并计算了一个8 0 英尺高的碾压混凝土模型坝的温度应力，以及其开裂情况。限于当时计算机硬件水平， 第3 页 河海大学硕士学位论文 这种计算只是尝试性的。 日本学者考虑了混凝土徐变应力场的计算远远复杂于温度场，所以他们首先用有限 元或差分法计算坝体温度场，找到几个特征温差后，再利用a d i n a 程序计算三维应力场， 并以此预测了宫濑坝在施工期和运行期开裂的可能性。近几年来，日本学者在混凝土温 度徐变应力的物理仿真研究上，走在了世界的前列。他们通过大量试验资料证明：和大 体积紧密连接的应力计，可以方便地测出混凝土各部位的温度应力，并且认为，只要与 温度应力有关的材料参数的精度足够，其实测的温度应力的精度也足够。 中国在大体积混凝土结构的温度应力数值分析和理论研究方面，一直被认为是处于 世界前列的。这是我国坝工界，特别是中国水利科学研究、高等学校从5 0 年代起不懈 努力的结果。进入8 0 年代以来，中国水利水电科学研究院、清华大学、天津大学、河 海大学、武汉水利电力大学，四川联合大学、大连理工大学、西安理工大学等，都展开 了碾压混凝土坝温度应力的攻关研究，分别对沙溪口溢流坝、观音阁、铜街予、龙滩等 以建、在建和待建的碾压混凝土重力坝，进行了温度应力计算。上述单位通过理论研究 和工程实践，取得了一批有价值的成果。 中国水利水电科学研究院朱伯芳院士为减少碾压混凝土坝的计算工作量，提出了以 误差控制为特点的“扩网并层算法”、“分布异步长算法”；清华大学刘光廷教授等将断 裂力学的研究成果融入功能强大的仿真程序中，应用“人工短缝”成功地解决了溪柄碾 压混凝土拱坝两岸的温度拉应力问题。天津大学赵代深教授、李广远教授结合国家攻关 项目在混凝土坝全过程多因素仿真分析方面取得出色成果。河海大学在1 9 9 0 年至1 9 9 2 年间结合小浪底工程完成了大体积混凝土结构的二维、三维有限元仿真程序系统 ( t c s a p ) ，该程序系统具有较丰富的前后期处理和图形输出技术；河海大学傅作新、陈 里红、张子明等首次在温度应力仿真计算中考虑了混凝土的软化性能，在龙滩碾压混凝 土坝的温控研究中，建立了一、二、三维有限元综合分析的数值模型，取得了较为理想 的分析成果。原武汉水利电力大学王建江在“八五”攻关项且中提出了旨在减少网格数 量的“非均质单元法”，肖明提出了考虑外部温度变化效应的三维损伤开裂非线性有限 元分析方法，陈敏林提出了混凝土开裂跟踪的估算应力方法。大连理工大学在碾压混凝 土拱坝的三维应力场仿真分析中，根据大坝分层的施： 特点，将混凝土水化热项单独考 虑，提出了仿真分析的“波函数法”。西安理工大学陈尧隆等提出的浮动网格法。 这些研究结果大多是结合具体工程进行研究，因此形成了各具特色的研究方法。从 各家仿真计算结果来看，温度场的计算结果非常接近，但应力成果有一些差别。 钨4 页 第一章绪论 1 3 问题的提出和本文的主要研究内容 1 3 1 问题的提出 在碾压混凝土坝的设计中，通常要求不出现拉应力或者出现很小的拉应力，对于自 重、水压力等外荷载，要做到这一点并不困难。但是施工过程中在气温、水温、自身水 化热、水管冷却以及新老混凝土相互作用等复杂因素下，混凝土温度变化较大，引起较 大温度变形，当温度变形受到强制约束时，便会在混凝土内产生较大的温度应力，要把 这种温度变化所引起的拉应力限制在允许范围以内是非常不容易的，尤其是对碾压混凝 土高拱坝。实践证明，多数碾压混凝土坝在施工期间即己开裂或者在运行期间出现裂缝， 因此，碾压混凝土坝的温度应力和温度控制问题，即考虑施工过程影响的温度应力仿真 分析已经成为一个重要研究课题。 在国内碾压混凝土大坝的施工过程中，通常每连续碾压多层厚为0 _ 3 o 5 m 的混凝 土本体薄层后，形成一个升程，应混凝土表层散热等需要，升程与升程之间通常有多天 的施工间歇时间。由于水泥的水化硬化反应，混凝土浇筑后，其材料参数随时间而变化， 具有强非均质性，在数值计算中需要对这种非均质性进行尽可能精细的模拟。但是由于 坝体中层面之间的本体层厚很小，若按常规有限单元法，计算规模太大，无法对大坝施 工期和运行期的非稳定温度场和应力场进行有效的数值仿真计算。 针对这个问题，朱伯芳院士首先提出并层算法，采用拟均质单元。王建江提出“非 均匀单元法”，把阶梯型分布的材料参数简单近似地表示成单元局部坐标的连续函数， 但相邻两个升程之间材料参数有较大的突变，这样处理必然会引起较大的误差。刘光廷 提出基于位移等效的等效连续模型，也能够加大有限元网格尺寸。有关“虚拟层合单元” 的文章所介绍的层合单元法中，基本概念与非均质单元法相同，但是其中的坐标变换法 单元数学处理方法可使单元劲度或传导矩阵的计算精度得到提高。后朱伯芳院士又提出 了“并层复合单元法”。本文运用“非均质单元层合单元法”，在这个计算方法中，不但 单元是非均质的，而且这个单元是由几个非均质升程( 子非均质层) 组成，也可退化成 只有个浇筑层的均质单元，在整个计算域的有限元仿真计算过程中，大大减少了单元 数和结点数，有效提高了计算效率。 1 3 2 本文主要研究内容 本文主要围绕高碾压混凝土坝温度场和应力场的三维仿真分析、混凝土热学参数的 第5 页 河海大学硕士学位论文 确定以及碾压混凝土坝的温控防裂研究，主要完成了以下几方面的研究内容： 1 本文在查阅了大量国内外资料的基础上，综述了碾压混凝土坝的发展、温度应力分 析方法及有限元分析方法的现状。 2 综述了热传导的基本理论及温度场和弹性热应力问题的有限元求解方法，阐述了温 度单元和应力单元的匹配问题。 3 碾压混凝土坝为分批层层浇筑，浇筑时间不同，各浇筑块的物理参数，热学边界条 件都要发生相应的变化，为了对施工过程进行精确仿真，计算过程中需要反映材料 参数，边界条件和网格随时间和空间变化的特性。常规的有限元分析方法是难以实 现上述要求的，因此需要进行计算机仿真技术的研究和探讨。作者详细阐述大体积 混凝土分层施工过程中的温度及应力仿真计算的方法，并编写了相应的有限元计算 程序，通过算例验证了其正确性。 4 要对高碾压混凝土坝进行精确仿真计算，工作量非常巨大，对此本文运用“非均质 层合单元法”，采用分层积分法计算合并后的网格单元，进行了相应模型的计算。 5 反分析法是利用工程中的实测值，通过数值计算确定混凝土的材料参数，进而使计 算结果尽量接近实际工程情况。本文采用遗传算法反演混凝土热学参数，并结合立 方体试验进行了反演计算。 6 本文对马堵山碾压混凝土重力坝温控防裂研究这一科研项目深入开展了有关理论和 实际运用方面的研究工作。 第6 页 第二章碾压混凝土坝三维非稳定温度场和应力场的仿真分析 第二章碾压混凝土坝三维非稳定温度场和应力 场的仿真分析 根据国内外己建和在建的混凝土坝运行情况来看，几乎所有的混凝土坝都或多或少 地出现裂缝，而且这些裂缝大部分是由温度应力引起的，因此温度裂缝的控制是混凝土 拱坝技术的关键因素。为了掌握混凝土坝温度场和温度应力场的分布规律，了解大坝温 度裂缝的开裂机理，有必要对混凝土坝进行温度场和应力场的仿真计算。对于高碾压混 凝土坝，由于其特殊的施工特点，给坝体温度场和应力场的仿真计算带来了很大的困难， 因此很有必要寻找一种既能减少计算规模，又能满足精度要求的数值计算方法。 本章首先简单地介绍了热传导基本理论，接着阐述了三维稳定和不稳定温度场解 析解的求解方法和有限单元法分析的基本理论，进一步提出了适合于高碾压混凝土坝施 工特点的非均质层合单元法仿真分析模型，最后介绍了混凝土应力有限元分析的基本原 理和仿真分析方法。 2 1 热传导基本理论订1 2 1 1 热传导方程 热传导是一种特定方式的传导，这种传导方式的能量交换发生在固体或静止流体 ( 没有由于介质的宏观位移引起的对流运动) 内，依靠物体内部的温度梯度从高温区域 向低温区域传输热量。 在一般三维问题中，瞬态温度场t ( x ，y ，z ，) 在直角坐标系中应满足的微分方程是 p 。百o t 一夏o 【t 瓦o t ) 一言( 屯面o t ) 一瓦k o 出t ) - p q = 。( 在n 内) ( 2 1 ) 边界条件是 t = t ( 在r 。边界上) ( 2 2 ) t 娑r + i 。娶n + t 娶以：q ( 在r ：边界上) ( 2 3 ) t 瓦r + 砖瓦b + t i ”z 2 q 。征12 边界上 ( n 3 j 匕界边r在 d l “ = 心 卯一出 页 t 7 + 第 刀一砂 b + k 望缸 t 河海大学硕士学位论文 式中，p 是材料密度( k g l m 3 ) ；c 是材料比热容( a i ( k g 世) ) ；f 是时间( s ) ；k x , k ，k z 是材料沿物体3 个主方向( x ，y ，：方向) 的导热系数( w i ( m k ) ) ；q = q ( x ，y ，z ，) 是物 体内部的热源密度( 豫) ；q ，h ，r 是边界外法线的方向余弦；t = r ( r ， ) 是在r ，边界 上的给定温度；g = q ( r ，f ) 是在r ：边界上的给定热流密度( w m 2 ) ；h 是对流换热系数 ( w m 2 k ) ；瓦= l ( f ，f ) ，对于f ；边界，在自然对流条件上，t 是外界环境温度；在 强迫对流条件下，是边界层的绝热壁温度。 2 1 2 稳态热传导有限元的一般格式 1 ) 稳态热传导问题的变分原理 采用加权余量的伽辽金方法可以推导与三维稳态热传导微分方程相等效的变分原 理，给出变分原理中的泛函表达式如下： 兀= l 阻势始心1 ( 簪职卜j ， 一f 2 。砑一h 瓦一互1 丁p 可阻验证，从上述泛函的驻值条件占丌( r ) = o 得到泛函的欧拉方程就是三维稳态热传导 问题的微分方程和边界条件。也就是表示泛函的变分是和该微分方程及边界条件相等效 的。因此利用占r i ( r 1 = o 可以建立有限元的求解方程。 2 ) 有限元离散 将求解域q 离散为有限个单元体，在典型单元内各点的温度丁可以近似地由单元的 结点温度值z 插值得到 t = 于= m ( x ，y ，：) 互= n t 。 ( 2 6 ) i = 1 其中 。是每个单元的结点数；m ( 毛y ，z ) 是c 0 型插值函数，它亦具有下述性质： 嘶，= 0 兽箬 弦， 第8 页 ， 第二章碾压混凝土坝三维非稳定温度场和应力场的仿真分析 并且 雌 m = 1 ，= 】 ( 2 8 ) 将( 2 6 ) 代入有限元离散后的泛函，从占兀= 0 ，司以得到稳态热传导问题的有限兀求解 方程 k t = p( 2 9 ) 式中k 是热传导矩阵，它是对称矩阵，在引入给定温度条件以后，k 是对称正定的： t = 【正五正】是结点温度列阵a 矩阵k 和p 的元素分别表示如下： g i j = 7 z f ( ko 出n , 砂叶警等+ k z o 岔n io 昆n j d f ，+ 孙b u r n , 州z 1 0 ) 鼻2 n , q a r + 莩l n , h l d r + 军上。n , p q d f 2 ( 2 1 1 ) ( 2 1 0 ) 式中的第l 项是各单元对热传导矩阵的贡献，第2 项是第- 2 类热交换边界条件对 热传导矩阵的修正。( 2 1 1 ) 式中的3 项分别为给定热流、热交换以及热源引起的温度 荷载。可以看出热传导矩阵和温度荷载列阵都是由单元相应的矩阵集合而成。可将( 2 1 0 ) 及( 2 11 ) 式改写成单元集成的形式： = 蟛+ 巧 ( 2 1 2 ) p = 鬈+ + b 。e ( 2 1 3 ) 其中 巧= l 八k 。o 缸n 。砂n s 一+ b 可o n 百o n + k ：o n io 出n j i d q ( 2 1 4 ) 域e = h n , n r ( 2 - 1 5 ) 鬈= i n , q d f ( 2 - 1 6 ) 噶= f ：n h t o d f ( 2 1 7 ) p * f = n , p q d f ) ( 2 1 8 ) 以上是三维稳态热传导问题的有限元方程，在集成系统刚度矩阵以后，需要引入至 少限制“刚体运动”的给定位移条件。对于稳态热传导问题，在集成系统刚度矩阵以后， 还需要引入给定一个点温度的条件。 河海大学硕士学位论文 2 1 3 瞬态热传导问题 1 ) 瞬态热传导有限元的一般格式 瞬态温度场与稳态温度场的主要差别是瞬态温度场的场函数温度不仅是空间域q 的 函数，而且还是时间r 的函数。 首先建立三维瞬态热传导问题的微分方程和边界条件的等效积分形式，即 1 1 w 百o t 一昙c t 罢卜参c 勺万o t 卜昙c t 詈卜加 q + f 。( 丁一于妒+ f ：( t 罢q + b 孑b + t 老心一。户r + c z f ：码 吒罢n ，+ b 号b + t 警”：一n ( 瓦一丁) p r = 。 其中w ，w 2 ，是任意函数。按照伽辽金方法选择任意函数，设r 。上z f 懈t c t = 于 则= 0 ，并且不失一般性地可设 ( 2 2 0 ) 将上式代入( 2 2 0 ) 式，并对其中第一个域q 内积分的第2 4 项进行分部积分，则可 得到 l j ， p c 詈 + 警( t 罢) + 等 勺琴 + 塑o zf k 七：竺a z ) - s t p q q 。：。， 一f ：嘞d r f ，8 r h ( r 。一t ) a r = 0 利用上式可以建立瞬态温度场有限元的一般格式。首先将空间域q 离散为有限个单元 体，在典型单元内温度丁仍可以近似地用结点温度z 插值得到，但要注意此时的结点温 度是时间的函数，即 丁= 于= 芝m ( 墨y ，z ) z ( f ) ( 2 2 2 ) 插值函数m 只是空间域的函数，将上式代入( 2 2 1 ) 式，并考虑到j z 的任意性，就可 以得到用来确定 个结点温度z 的有限元求解方程 c 童+ k t ：p 第1 0 页 ( 2 2 3 ) 、 r 第二章碾压混凝土坝三维非稳定温度场和应力场的仿真分析 这是一组以时间r 为独立变量的线性常微分方程组。式中c 是热容矩阵，k 是热传导矩 ，阵，c 和k ( 在引入给定温度条件以后) 都是对称正定矩阵：p 是温度荷载列阵，t 是 结点温度列阵，t ( = 鲁 是结点温度对时间的导数列阵。矩阵k ，c 和p 的元素由单元 的相应的矩阵元素集成，即 巧= 巧+ 弼 ( 2 2 4 ) q = g ( 2 2 5 ) 只= 呓+ 鬈+ 辱 ( 2 2 6 ) 式中，巧e 是单元对热传导矩阵的贡献，蟛是单元热交换边界对热传导矩阵的修正，q 是 单元对热容矩阵的贡献，露是单元的对流换热边界的温度荷载。这些单元的矩阵元素 由下列各式给出： k ，e 八k 。a 出n c 。砂n j k ，c 3 砂n ic 3 砂n _ ，一k 。o 如n ic 3 赴n j d q ( 2 2 7 ) 蟛= l 帆啊d 1 1 ( 2 2 8 ) c ；2k p c n t n q q 2 9 ) 焉2 i n , q d f ( 2 3 0 ) 巧3 f ；n , h z d r ( 2 3 1 ) 呓2 m p 鲥q ( 2 3 2 ) 至此，以将时间域和空间域的偏微分方程问题在空间域内离散为n 个结点温度z ( r ) 的常 微分方程的初值问题。对于给定温度值的边界r ，上喝个结点，方程组( 2 2 3 ) 式中的相应 项应引入的条件是 i = 正( f = 1 ，2 , 3 一，啊) ( 2 3 3 ) 式中，i 是r 上碍个结点的编号。 2 ) 一阶常微分方程组的求解方法 瞬态热传导问题的有限元求解方程( 2 2 3 ) 式是一阶常微分方程组。对此方程组的求 箱1 i 页 河海大学硕士学位论文 解，在各种数值方法的著作中有很多通用的和特殊的解法可供选择。本文主要讨论直接 直接积分法是指求解常微分方程组之前不对它的形式进行变换，而是直接对它进行 数值积分。通常基于两个基本概念，一是将求解的时间域0 ，t 划分为若干个时间步 长出，在一定数目的出时间区域内，假设t 和t 的函数形式来近似方程的精确解；二是 以仅在相隔出的离散时间点上满足微分方程来代替时间域内任何时刻f 都满足微分方 将时间域0 丁等分为m 个时间步长，a t = t i m ，并认为t = o 时的初始温度列阵磊是 己知的。进一步假设t = 0 ，( = a t ) ，乞( = 2 a t ) ，f 。( = n a t ) 时刻的解l 已经求得。下一步要 计算的是f 。( = t + 出) 时刻的瓦。从该求解过程建立起求解所有离散时间点场函数的一 可以在两个时间点f 。和+ ，之间的出时间区域内，假漫7 1 和于采用如下的线性插值形式， t ( 乙+ 弘r ) = m l + m + i l 。 ( 2 3 4 ) p ， 弘石 ( o 茎善1 ) n ，= 1 一n + 。= 亍= 为。l + n o + 。l 。 ( 2 3 5 ) 瓤= 一石1吃+ 。= 石1 时间区域，典型的加权余量格式可以表示为如下形式： f w c ( 也t + 取+ 。l + 。) + k ( n o t o + m + 。l + ，) 一p ，f = o ( 2 3 6 ) 当求解初值问题时，如果已知一组参数l ，则可以利用( 2 3 6 ) 式近似确定另一组参数 第二章碾压混凝土坝三维非稳定温度场和应力场的仿真分析 l + 1 0 将( 2 3 4 ) 和( 2 3 5 ) 式代入( 2 3 6 ) 式就可以得到时间区域血前后结点上两组参数的 ，7 。 关系式，即 ( k f 蟛咖c f w 筹 k + ( k f w ( 1 吲咖c f w 警卜f w p 咖吣 式中可以代入不同的权函数。在以上的讨论中假定热传导矩阵k 和热容矩阵c 不随时间 t 而变化。( 2 3 7 ) 式可以表达为任何权函数都适用的一般形式，即 ( c t + k o ) t + 。+ 一c a t + k ( 1 一目) l = f ( 2 3 8 ) 式中 曰= f 喈d 善f w 蝣f = f w p d 喜f w 蟛 ( 2 3 9 ) 一种很方便的做法是假定p 采用与未知场函数t 相同的插值表达式，这时将得到 乒= 匕+ ，目+ 匕+ ，( 1 - 0 ) ( 2 4 0 ) 当己和p 都已知时，就可由( 2 3 8 ) 式求得下一时刻的e + 。这就是两点循环公式。可以 更明显地表示成 k l “= q = 0 ，1 ，2 ，m ) ( 2 4 1 ) 其中 霞= c a t + 0 k ( 2 4 2 ) q 。= c f 一( 1 一目) k l + ( 1 一目) l + 巩+ ( 2 4 3 ) 利用上式，从t = 0 出发，可以依次递推求得结点温度列阵t ( f ) 的各个瞬时值 王，l ，瓦，。 2 2 热应力的计算 当物体各部分温度发生变化时，物体由于热变形将产生线应变口( r 一矗) ，其中口是 材料的线膨胀系数，t 是弹性体内任一点现时的温度值，瓦是初始温度值。如果物体各 部分的热变形不受任何约束，则物体发生变形将不会引起应力。但是如果物体受到约束 或者物体由于温度变化不均匀，使得物体的热变形不能自由进行时，则会在物体中产生 第1 3 页 河海大学硕士学位论文 应力。物体由于温度变化引起的应力称为“热应力”或“温度应力”。当弹性体的温度 场已经求得时，就可以进一步求出弹性体各部分的热应力。 2 2 1 弹性热应力问题的有限元方程 物体由于热膨胀只产生线应变，而剪切应变为零。这种由于热变形产生的应变可以 看作是物体的初应变岛。对于三维问题，氏的表达格式是 g o = 口( ，一兀) 1 1 10 0 o l 。 ( 2 4 4 ) 其中，口是材料的热膨胀系数( 1 o c ) ；瓦是结构的初始温度场；t 是结构的稳态或者瞬 态温度场。丁可由温度场分析得到的单元结点温度f 通过插值求得，即可按( 2 , 6 ) 式计 算得到： n 。 丁= 亍= 芝m ( x ，y ，z ) z = n t 8 ( 2 6 ) = l 在物体中存在初应变的情况下，应力应变关系可表示成 仃= d ( 一o ) ( 2 4 5 ) 将( 2 4 5 ) 式代入虚位移原理的表达式，则可得到包括温度应变在内的用以求解热应力问 题的最小位能原理，它的泛函表达式如下： 兀小) = 北s 。e 。可巾q f 。u 嘞i ( 2 - 4 6 ) 将求解域q 进行有限元离散，就得到有限元求解方程为 k a = p( 2 4 7 ) 与不包含温度应变有限元求解方程的区别在于荷载向量中包括了由温度应变引起的温 度荷载。这里荷载向量表达为 p = p ，+ p 7 + p ( 2 4 8 ) 其中p ，l 分别是体积荷载和表面荷载引起的荷载项，是温度应变引起的荷载项， 匕= l b 7 d 。d f f 2 ( 2 4 9 ) 从以上各式可见，结构热应力问题和无热荷载的应力分析问题相比，除增加一项以初应 变形式出现的温度荷载匕以外，则是完全相同的。稳态温度应力计算在温度场分析后 第1 4 页 第二章碾压混凝土坝三维非稳定温度场和应力场的仿真分析 进行。至于瞬态温度应力的计算，可以在每一时间步的瞬态温度场计算后进行，也可以 在整个瞬态温度场分析完成后，再对每一时间步或指定的若干步进行，这可根据实际需 要和计算的方便与否决定。 2 2 2 温度单元和应力单元的匹配 从数据准备和有限元分析本身来说，都希望温度场分析和应力场分析采用同一个 有限元网格。由结点位移计算得到的应变占和由结点温度计算得到的温度应变s 。出现在 同一个关系式中，考虑到是占由晶引起的，它们应该是同阶次的。但是是由位移场取 导数得到的，而毛是直接由温度场得到的，所以合理的匹配，对于c 0 型单元应该是位 移场的插值函数比温度场的插值函数高一个阶次。例如二维问题，应力分析采用8 结点 二次单元，温度场分析采用4 结点双线性单元。另一种替代方案是温度场分析的单元插 值函数和应力分析的单元插值函数阶次相同，但在由结点温度计算温度应变矗时采用低 一阶次的插值函数。例如二维8 结点二次单元的情形，热应力分析时单元刚度矩阵和温 度荷载列阵如果采用2 2 积分，则e 0 实际上已采用了低一阶的双线性插值函数进行计 算。但是如果采用3 x 3 积分，则不满足占和矗同阶次的要求，从而影响计算结果的精度。 2 3 混凝土仿真应力场的有限元计算阴 根据弹性徐变理论，用增量初匝燹；h - 法计算施工瑚、运行期内由于变温等凼素而引 起的徐变应力变化规律。应变增量包括弹性应变增量、徐变应变增量、瞬态徐变应变增 量、温度应变增量、干- n 应变增量和自生体积变形增量，即 毛 = s ： + p ： + s ：。 + s ： + s ： + p ： ( 2 5 0 ) 式中： 占：) 一弹性应变增量； s ：卜_ 徐变应变增量； 占乒 - 瞬态徐变应变增量； s j 一温度应变增量； s ： 一干缩应变增量； s ：卜一自生体积变形增量。 弹性应变增量 s ：j 由下式计算 2 南豳】 吒 ( 乞2 竿，以- l - n ) ( 2 5 1 ) 河海大学硕士学位论文 1 一一0 l一“0 10 对2 ( 1 + ) 称 0 o 0 0 2 ( 1 + ) o 0 0 o o 2 ( 1 + 弹性模量e ( ) 可采用朱伯芳提出的双指数式求出 f 2 5 2 ) e ( f ) = e o ( 1 一e - o r 6 ) ( e o 为最终弹性模量) ( 2 5 3 ) 徐娈应娈增量见式( 2 7 0 ) ，瞬态徐受应受增量见式( 2 7 1 ) 温度应变增量 ：) 由非稳定温度场计算结果求得，求出温度场后可由下式求得 j = a a t ，& ，n ，t ，0 ，0 ，o ) ( 2 5 4 ) 式中：口一热膨胀系数，如果是各向异性材料，则三个方向的d 取值也不同：一温 度差。 干缩应变增量 s ： 由下式计算 s ：) = 一鼹) ( 2 5 5 ) 其中砖 按式t ( r ) = 五笔。计算。 在任一时刻a t , 内，由弹性徐变理论的基本假定可得增量形式的物理方程： 吒) = e 】 钆卜 玑 一 e 一 。： 一 s ： 一 s ： ) ( 2 5 6 ) 瓯 = e 妣e = 而黧鬲 ( 2 t 5 7 ) 根据虚功原理，不难导出计算混凝土结构在热、力作用下单元的温度徐变应力矩阵 方稗： 【k 。】。 民) 。= ”+ 破卜 裤r + 她r + 卜 砖 8 式中： 第1 6 页 r 2 5 8 ) 第二章碾压混凝土坝三维非稳定温度场和应力场的仿真分析 岛 。为单元刚度矩阵： 暇 。= f 别7 哦伸】施 占。) 。为单元结点位移增量向量； 蛾。r 为瞬态徐变引起的单元结点荷载增量 掣。r = i 捌7 哦】 s ：。细 n 。 蟛r 为基本徐变引起的单元结点荷载增量 蟛r = f 【司7 见】 s 品妞 迁 譬r 为温度引起的单元结点荷载增量 r = 例7 【见 t 汹 n 瞬y 为自生体积变形引起的单元结点荷载增量 时y = 吲7 d 。 矗0 细 皿 y 为干缩变形引起的单元结点荷载增量 蚪y = f 【b 九b 】 s ：) 地 群r 为外荷载引起的单元结点荷载增量 蚪r = 别7 见 占细 n ， 式中： 陋 = 陋。b ：b ，b