徐 珍，費文平，彭仕麒

(1.四川大學(xué)水力學(xué)與山區(qū)河流開發(fā)保護國家重點實驗室，成都 610065；2.四川大學(xué)水利水電學(xué)院，成都 610065)

大體積混凝土是一種廣泛應(yīng)用于混凝土壩、房屋建筑、橋梁等大型土木結(jié)構(gòu)工程的材料[1]。由于其水化生熱的材料特點及熱量不易散發(fā)的結(jié)構(gòu)特點，大體積混凝土在不均勻溫度場和外部約束的條件下，會產(chǎn)生較大的溫度應(yīng)力，并導(dǎo)致混凝土的開裂，嚴重影響結(jié)構(gòu)安全和經(jīng)濟效益。目前工程界常利用有限元軟件ANSYS的熱分析模塊及耦合單元進行溫度應(yīng)力的仿真模擬。

自1927年，朱伯芳院士編制出了我國第一個不穩(wěn)定溫度場有限元程序以來，我國學(xué)者結(jié)合具體工程開展了各類仿真研究工作。武漢大學(xué)的彭詩明采用二維直接耦合方案進行三峽工程三期RCC圍堰施工期熱應(yīng)力的仿真計算[2]；天津大學(xué)的黃君寶通過三維直接耦合得到了上海長江大橋橋墩在溫度作用下應(yīng)力發(fā)展過程[3]；河海大學(xué)的陳長華為得到合理的施工方案，利用三維間接耦合的方案進行了三叉河口閘墩的溫度場和溫度應(yīng)力的計算分析[4]。實際上由于二維模型不能模擬大壩的三維效應(yīng)以及三維直接耦合分析無法反映混凝土的塑性變形，上述這些方案存在局限性，然而目前缺乏耦合方法對熱應(yīng)力計算結(jié)果影響的研究。本文以云南省曼點水庫大體積混凝土重力壩為例，利用有限元仿真程序ANSYS提供的熱分析模塊及ADPL二次開發(fā)，討論在二維直接耦合、三維直接耦合及三維間接耦合這3種耦合方案下，大體積混凝土溫度場和應(yīng)力場計算結(jié)果的差異性，并給出了各方案的適用范圍，為以后在節(jié)約人力、時間的同時，快速、準確地得到熱應(yīng)力計算的準確結(jié)果提供參考。

1 基本理論

1.1 熱傳導(dǎo)理論

1.1.1 熱傳導(dǎo)方程

混凝土內(nèi)部溫度場的計算是指熱傳導(dǎo)方程在特定邊界條件和初始條件下的求解，整個傳熱過程為瞬態(tài)傳熱。三維瞬態(tài)溫度場的熱傳導(dǎo)方程為[5,6]：

(1)

1.1.2 定解條件

初始條件：多數(shù)情況下，初始瞬時的溫度分布可以認為是常數(shù)，即當τ=0時：

T(x,y,z,0)=T0(x,y,z)=常數(shù)

(2)

邊界條件可由以下方式給出：

(1)第1類邊界條件：

(3)

當混凝土與水接觸時，表面溫度等于已知的水溫，就屬于這種邊界條件。

(2)第2類邊界條件：

(4)

(3)第3類邊界條件：

(5)

式中：τ為時間；θ為絕熱溫升；λ為導(dǎo)熱系數(shù)；β為表面放熱系數(shù)；Ta為外界環(huán)境溫度。

第3類邊界條件表示了固體與流體接觸時的傳熱條件，也就是混凝土邊界與空氣接觸處的邊界條件。

1.2 應(yīng)力場有限元理論分析

在低應(yīng)力水平下，混凝土處于熱線彈性階段。在大體積混凝土溫度應(yīng)力場分析中，普遍采用熱彈性本構(gòu)模型。應(yīng)變計算中采用疊加原理；在應(yīng)力場中采用有限元-差分法，即在空間域用有限元離散，在時間域用差分法離散，其中應(yīng)用最廣的為朱伯芳院士提出的隱式解法[5]。

在高應(yīng)力水平下，大體積混凝土處于非線性階段，即熱彈塑性階段。用有限元進行溫度應(yīng)力分析時，需要考慮塑性變形對溫度應(yīng)力的影響[6]。

1.3 熱-應(yīng)力耦合場分析

溫度應(yīng)力問題實際上是溫度和應(yīng)力2個物理場之間的相互作用[6]。ANSYS有限元分析軟件提供了2種熱分析的方法來計算溫度應(yīng)力：其一為間接耦合法，先進行熱分析，將得到的節(jié)點溫度作為荷載施加于后續(xù)的結(jié)構(gòu)分析模型中，然后進行結(jié)構(gòu)分析；其二為直接耦合法，直接耦合單元包括所有所需的自由度，可以通過一次求解得到耦合場分析結(jié)果。

其中，熱應(yīng)力直接耦合分析由于操作方便、能實現(xiàn)熱-結(jié)構(gòu)的雙向耦合，被廣泛使用。然而受商用分析軟件的限制，直接耦合單元不能考慮混凝土的塑性變形。因此，一方面可以借助三維間接耦合分析模擬大壩塑性變形；另一方面，考慮采用平面模型代替三維模型，實現(xiàn)建模簡單、計算快速且能考慮混凝土的材料特性的要求。

基于上述設(shè)想，本文分別用這3種耦合方案模擬大體積混凝土的熱應(yīng)力場，通過對比分析，探討這3種方案的差異性。

2 計算的關(guān)鍵技術(shù)

2.1 水泥水化熱的模擬

水泥水化熱是混凝土結(jié)構(gòu)的內(nèi)熱源[8]，它隨水泥的品種和用量、混凝土材料品種及澆注溫度的不同而異。實際上溫度場中用的是混凝土的絕熱溫升[9]。本文混凝土絕熱溫升選用指數(shù)型表達式：

θ=θ0(1-e-m τ)

(6)

式中：θ為混凝土絕熱溫升；θ0為混凝土最終絕熱溫升；m為常數(shù)，隨水泥品種、比表面積及澆注溫度變化，本文取為0.4；τ為時間。

2.2 混凝土的彈性模量

由應(yīng)力計算公式可知，溫度應(yīng)力的數(shù)值與彈性模量E(τ)成正比[10]。混凝土彈性模量隨時間變化的關(guān)系可取指數(shù)式、修正指數(shù)式或復(fù)合指數(shù)式，本文采用復(fù)合指數(shù)式：

E(τ)=E0(1-e-α τb)

(7)

式中：E(τ)為實時彈性模量；E0為初始彈性模量；取常數(shù)a為0.4，b為0.34。

本文通過APDL語言按齡期定義一系列的彈模數(shù)值，在計算過程中通過循環(huán)語句給每個齡期的混凝土彈模賦值。

2.3 水管冷卻的等效熱傳導(dǎo)方程

把水管看成負熱源，在平均意義上考慮水管冷卻作用，并用數(shù)值方法求出問題的近似解，可以把只考慮水管冷卻作用而計算的混凝土平均溫度作為絕熱溫升，從而得到等效熱傳導(dǎo)方程[11]：

(8)

該方法只是一個近似的方法，但實現(xiàn)計算比較方便，而且可以使用通常的計算網(wǎng)格，效率較高。

3 實例分析

以云南省曼點水庫工程為例，依次計算大壩在二維直接耦合分析、三維直接耦合分析及三維間接耦合分析下的溫度應(yīng)力場。

3.1 工程概況

曼點水庫大體積混凝土重力壩由非溢流壩段、溢流表孔、輸水底孔等建筑物組成，壩頂高程為1 090.0 m，河床建基面高程為1 012 m，最大壩高78.0m，壩軸線長225.0 m，壩頂寬8.0 m，下游壩坡1∶0.76，上游壩坡高程1 048.0 m以上鉛直，以下為1∶0.2。

3.2 有限元模型的建立和網(wǎng)格劃分

選取非溢流壩段及地基作為整體，建立有限元模型。采用笛卡爾直角坐標系，坐標系原點取在上游鉛直面與壩0+0樁號的交點，X方向為從上游指向下游，以上游鉛直面處為零點，Y方向為鉛直向上方向，以海拔零高程為零點，Z方向沿壩軸線方向，正方向由右手螺旋法則確定。在二維有限元分析中對基礎(chǔ)上下游側(cè)和底部施加絕熱邊界約束和法向位移約束，各澆筑層面臨空面施加對流邊界；在三維有限分析中，將基礎(chǔ)四周的鉛直面和底部施加絕熱邊界約束，各澆筑層面臨空面施加對流邊界。有限元結(jié)構(gòu)模型見圖1，有限元仿真分析計算范圍及模型參數(shù)見表1。

圖1 大體積混凝土有限元計算圖Fig.1 Finite element mesh of mass concrete

表1 3種方案的計算范圍及模型參數(shù)Tab.1 The calculation and model parameters in three case

3.3 計算參數(shù)與施工方案

采用Drucker-Prager屈服準則來模擬巖體和混凝土。混凝土及基巖的熱力學(xué)性能材料參數(shù)按表2取值。由于溫度應(yīng)力主要發(fā)生在施工期，本次分析主要考慮溫度應(yīng)力與自重應(yīng)力的組合，而不考慮與上游水荷載及地震荷載的組合。選定施工工況為計算工況，本工程澆筑進度及通水方案見表3。

表2 壩體材料參數(shù)Tab.2 Dam material parameters

表3 澆筑及通水方案Tab.3 Pouring and water solution

3.4 計算條件

(1)初始條件。壩基初始溫度為18 ℃；混凝土入倉溫度：1期混凝土澆筑溫度≤18 ℃；2期混凝土澆筑溫度≤20 ℃；3期混凝土澆筑溫度≤22 ℃。

(2)保溫措施。對壩體上、下游面采用濕草袋永久保溫；澆筑層頂表面無保溫材料。

(3)冷卻水管。采用HDPE管，外徑32 mm，壁厚2 mm，管材的導(dǎo)熱系數(shù)≥1.6 kJ/(m·h·℃)。每根水管的長度為200 m，水管的水平間距1.5 m，垂直間距為混凝土澆筑層厚度，冷卻水初溫為12 ℃，流量為1.0 m3/h，1期冷卻在開始澆筑混凝土?xí)r立即進行，每根水管通水時間為14 d。

3.5 計算結(jié)果及分析

根據(jù)所提供的力學(xué)參數(shù)及施工進度安排，本文對該工程進行了施工期全過程溫度應(yīng)力場的仿真計算。

3.5.1 耦合方案對溫度場影響的分析

選取第58層第9天進行溫度場分析。3種方案下的溫度場計算結(jié)果見圖2，溫度最大值、最小值及發(fā)生的部位見表4。

圖2 溫度場計算結(jié)果Fig.2 Calculation result of temperature field

耦合方案最大值/℃發(fā)生部位最小值/℃發(fā)生部位二維直接耦合28．4883壩體中部18．0924頂層混凝土兩端表層三維直接耦合27．8071頂部澆筑層中間15．7530頂層混凝土兩端表層三維間接耦合27．6803壩體中部15．7459頂層混凝土兩端表層

由表4可知3種方案下溫度場計算結(jié)果基本一致，這是因為大體積混凝土溫度應(yīng)力應(yīng)力場并非雙向耦合，結(jié)構(gòu)場的變化對溫度場的影響較小。

3.5.2 耦合方案對應(yīng)力場計算結(jié)果的影響

本文分別選取了3種耦合方案下2個典型澆筑時段(第1層第3天及第58層第1天)第1主應(yīng)力的計算結(jié)果(圖3～圖5)，它們分別表示了低應(yīng)力及高應(yīng)力水平下溫度應(yīng)力計算結(jié)果。另外，圖6給出了二維直接耦合分析及三維間接耦合分析的等效塑性應(yīng)力云圖。

圖3 二維直接耦合第1主應(yīng)力等值線Fig.3 The first principal stress contour map by 2D direct coupling

圖4 三維直接耦合第1主應(yīng)力等值線Fig.4 The first principal stress contour map by 3D direct coupling

圖5 三維間接耦合第1主應(yīng)力等值線Fig.5 The first principal stress contour map by 3D indirect coupling

圖6 等效應(yīng)力云圖Fig.6 Plastic equivalent stress nephogram

本文從模型維度及耦合方法入手，分別討論這3種方案對熱應(yīng)力計算結(jié)果的影響。對比方案見表5。

表5 3種方案下第1主應(yīng)力結(jié)果對比分析Tab.5 Analysis of the results of stress in three case

通過對比1，發(fā)現(xiàn)在低應(yīng)力水平下2者的誤差率為86.8%，這是因為二維模型不能完全反映巖石地基對大壩底部的約束作用；而在三維模型中，因為壩基的約束作用，底部拉應(yīng)力增大；在高應(yīng)力狀態(tài)下，2者第1主應(yīng)力數(shù)值及出現(xiàn)部位相差不大，且數(shù)值均小于混凝土抗拉強度，這反映了混凝土的塑性變形特性。另外由圖6可知，三維模型塑性區(qū)分布范圍更廣，且主要出現(xiàn)在壩址及壩體頂部。這一方面反映了大壩底部地基約束及應(yīng)力集中的影響，另一方面是由于表層對流邊界條件的存在，頂部澆筑塊降溫較快，溫度梯度較大，產(chǎn)生更大的溫度應(yīng)力，而這2點均無法在二維模型中體現(xiàn)。

在三維模型的條件下討論直接耦合與間接耦合的差異。對比2的結(jié)果表明，低應(yīng)力狀態(tài)下，由于混凝土保持為完全彈性體，2種方案下的第1主應(yīng)力數(shù)值、發(fā)生部位和分布情況相似，模擬效果一致。高應(yīng)力水平下，2者誤差率為107%，且應(yīng)力分布情況相差甚遠。究其原因，在高應(yīng)力狀態(tài)下混凝土已經(jīng)進入熱彈塑性體狀態(tài)，而三維直接耦合單元SOLID5卻不能模擬混凝土的這一材料特性，因此會出現(xiàn)第1主應(yīng)力值遠超過混凝土抗拉強度的情況，這顯然與混凝土實際受力狀態(tài)不符；而三維間接耦合的結(jié)構(gòu)單元SOLID45能夠模擬混凝土的塑性特性，并且在結(jié)果中反映了壩體逐漸擴大的塑性區(qū)。

針對這3種方案的差異性，本文給出了相應(yīng)的適用范圍，以期在以后研究中，學(xué)者能針對具體工程采用相應(yīng)耦合方案，以減小差異性帶來的影響。

(1)二維直接耦合方案建模簡單、計算快速且能反映混凝土的塑性特性，是一種良好的耦合分析方法。但二維直接耦合有限元并不能模擬整個壩段的實際邊界條件及受力情況，因此可用于計算各應(yīng)力水平下模型斷面形式單一且邊界條件較簡單的工程的熱應(yīng)力場。

(2)三維直接耦合分析能夠模擬復(fù)雜的地質(zhì)條件，并且由于使用方便，能提高計算的效率、節(jié)約工作量，建議用于應(yīng)力水平較低的各種工程。

(3)在三維間接耦合中，有限元程序能夠使用所有熱分析模塊的功能和結(jié)構(gòu)分析的功能，因而能很好地模擬實際受力情況并反映混凝土的材料特性，但是其計算模型單元、節(jié)點多，計算邊界較復(fù)雜、分析過程耗時較長，建議用于體型不規(guī)則、邊界及受力條件復(fù)雜的重大高壩工程。

4 結(jié) 語

本文以經(jīng)典熱傳導(dǎo)理論為基礎(chǔ)，通過大型有限元軟件ANSYS的熱分析模塊和參數(shù)化設(shè)計語言APDL，模擬了符合實際的混凝土澆筑過程及相應(yīng)的邊界條件，討論了3種方案下溫度應(yīng)力場耦合分析的差異。主要結(jié)論如下。

(1)由于混凝土結(jié)構(gòu)中結(jié)構(gòu)場對溫度場的影響較小，3種耦合方案溫度場的計算結(jié)果基本一致。

(2)3種方案的應(yīng)力場分析結(jié)果有差異：一方面與三維模型相比，二維模型不能反映整個壩段實際受力和邊界條件，其第1主應(yīng)力數(shù)值及塑性區(qū)范圍有誤差；另一方面對于三維模型，由于耦合分析單元的限制，在高應(yīng)力水平下，直接耦合的計算結(jié)果不可信。

(3)二維模型適用于模型斷面形式單一、邊界條件簡單的工程；三維直接耦合方案建議用于較低應(yīng)力水平(混凝土保持彈性狀態(tài))下的各種工程；三維間接耦合模擬效果最佳，但由于建模困難、操作復(fù)雜，建議用于體形不規(guī)則、邊界及受力條件復(fù)雜的重大高壩工程。

(4)在本例3種模型中均未考慮混凝土徐變的影響。在以后的研究中將考慮在不同耦合方法的前提下，實現(xiàn)基于ANSYS二次開發(fā)的大體積混凝土溫度徐變應(yīng)力的對比分析。

□

[1] 郭之章，傅 華.水工建筑物的溫度控制[M].北京：水利電力出版社，1990.

[2] 彭詩明.大體積混凝土結(jié)構(gòu)有限單元法應(yīng)用研究[D].武漢：武漢大學(xué)，2005.

[3] 黃君寶.大體積混凝土溫度場應(yīng)力場仿真方法研究[D].天津：天津大學(xué)，2005.

[4] 陳長華.考慮鋼筋作用的水工結(jié)構(gòu)施工期溫度場與溫度應(yīng)力分析[D].南京：河海大學(xué)，2006.

[5] 朱伯芳，宋敬廷.混凝土溫度場及溫度徐變應(yīng)力的有限元分析[M].北京：水利電力出版社，1977.

[6] 彭立海.大體積混凝土溫控與防裂[M].鄭州：黃河水利出版社，2005.

[7] 李彬彬，王社量，蘇三慶.各種因素耦合作用下的大體積混凝土溫度應(yīng)力場仿真分析[J].安徽建筑，2009,(4)：136-138.

[8] 朱伯芳，水工混凝土結(jié)構(gòu)的溫度應(yīng)力與溫度控制[M].北京：水利電力出版社，1999.

[9] 李 東，潘育耕.混凝土水化熱瞬態(tài)溫度場數(shù)值計算過程中的水化放熱規(guī)律及水化速率問題[J].西安建筑科技大學(xué)學(xué)報，1999,(3)：27-29.

[10] 朱伯芳.混凝土的彈性模量、徐變度與應(yīng)力松弛系數(shù)[C]∥ 水工結(jié)構(gòu)與固體力學(xué)論文集：朱伯芳論文集. 北京：水利電力出版社，1988.

[11] 朱伯芳.考慮水管冷卻效果的混凝土等效熱傳導(dǎo)方程[J]. 水利學(xué)報，1991,(3).