鐘 紅，賀 帥，李紅軍

（1.大連理工大學(xué)建設(shè)工程學(xué)部，遼寧省大連市 116024；2.中國水利水電科學(xué)研究院，北京海淀區(qū) 100044）

基于比例邊界有限元的溫度斷裂研究

鐘 紅1，賀 帥1，李紅軍2

（1.大連理工大學(xué)建設(shè)工程學(xué)部，遼寧省大連市 116024；2.中國水利水電科學(xué)研究院，北京海淀區(qū) 100044）

溫度荷載是導(dǎo)致結(jié)構(gòu)開裂、影響結(jié)構(gòu)穩(wěn)定性的一個(gè)重要因素。本文基于比例邊界有限元方法提出了溫度斷裂分析的模型。利用比例邊界有限元法求解結(jié)構(gòu)溫度場，結(jié)構(gòu)分析中的溫度荷載表示為徑向坐標(biāo)的函數(shù)并等效為若干項(xiàng)冪級(jí)數(shù)之和，基于線性疊加原理求解結(jié)構(gòu)響應(yīng)。該模型的特點(diǎn)包括：只需離散求解域的邊界從而降低了求解問題的維數(shù)；在求解裂縫尖端的溫度場、位移場和應(yīng)力場時(shí)均可獲得半解析解，應(yīng)力強(qiáng)度因子求解精度高。通過多個(gè)板承受溫度荷載的溫度場分析與結(jié)構(gòu)分析可以看出利用較少的網(wǎng)格即可獲得較高的精度，模擬所得的應(yīng)力強(qiáng)度因子隨著裂縫長度的增長而增長，并且本模型計(jì)算的應(yīng)力強(qiáng)度因子與現(xiàn)有文獻(xiàn)結(jié)果一致，證明了本模型的有效性。

溫度荷載；比例邊界有限元；溫度場；應(yīng)力強(qiáng)度因子

0 引言

隨著水電戰(zhàn)略的實(shí)施，高大混凝土壩不斷地建設(shè)，混凝土壩溫度斷裂問題也越來越突出，特別是大體積混凝土壩，由混凝土水化熱、外界氣溫、澆筑溫度、太陽輻射等引起的溫度荷載經(jīng)常會(huì)引起較大的拉應(yīng)力，其數(shù)值甚至可能超過混凝土的抗拉強(qiáng)度，從而導(dǎo)致混凝土產(chǎn)生裂縫，破壞大壩的整體性，降低耐久性，危及壩體的安全。國內(nèi)外學(xué)者開展了大量的研究，在溫度場與溫度裂縫擴(kuò)展的仿真模擬方面取了一些進(jìn)展，如Wilson編制了求解大體積混凝土結(jié)構(gòu)的二維溫度場的有限元程序DOT-DI-CE，并用該程序計(jì)算了德沃歇克（Dworshak）壩的溫度場；隨后美國陸軍工程師對有限元程序DOT-DICE進(jìn)行修改，并利用此程序計(jì)算大體積混凝土結(jié)構(gòu)的溫度場，將研究的成果應(yīng)用到Willow Creak大壩的溫度場分析中，得到了不同時(shí)期的溫度場數(shù)據(jù)，與實(shí)際情況十分符合[1]；程井與常曉林等人[2]采用無單元伽遼金法求解水工結(jié)構(gòu)應(yīng)力強(qiáng)度因子及仿真模擬溫度裂縫擴(kuò)展過程，在計(jì)算過程中不需要重新剖分網(wǎng)格，可以自由調(diào)整節(jié)點(diǎn)分布，能夠動(dòng)態(tài)地追蹤裂縫擴(kuò)展且單元邊界不影響裂縫擴(kuò)展路徑；Duflot M采用擴(kuò)展有限元方法求解了裂紋板的應(yīng)力強(qiáng)度因子及溫度場，得到了精確的結(jié)果[3]。

對溫度斷裂的研究包括溫度場求解和斷裂場模擬兩方面。通用的有限元法對于斷裂過程的模擬存在比較大的困難。比例邊界有限元法是一種新型的半解析的數(shù)值方法，它繼承了邊界元（BEM）和有限元（FEM）方法的優(yōu)點(diǎn)，只需對邊界進(jìn)行離散，不需要對整個(gè)區(qū)域進(jìn)行離散，與邊界元相對比，比例邊界有限元不需要基本解，而與有限元相比，其需要較少的節(jié)點(diǎn)，可以提高計(jì)算效率[4]。在處理斷裂問題時(shí)具有突出優(yōu)勢，無需基本解和奇異積分，在開裂方向上位移和應(yīng)力具有解析解，因此從奇異應(yīng)力模態(tài)直接提取應(yīng)力強(qiáng)度因子精度很高。宋崇民提出了在溫度荷載作用下廣義應(yīng)力強(qiáng)度因子的定義和計(jì)算方法，能得到整個(gè)區(qū)域準(zhǔn)確的溫度場、位移場及應(yīng)力場。

本文將比例邊界有限元方法引入混凝土結(jié)構(gòu)的溫度斷裂分析，開展了結(jié)構(gòu)的熱分析與結(jié)構(gòu)分析，并以多個(gè)算例驗(yàn)證了模型的有效性。

1 基本原理

1.1 溫度場求解

本文基于比例邊界有限元法（Scaled Boundary Finite Element Method，SBFEM）開展溫度場分析和結(jié)構(gòu)分析。該方法的控制方程的推導(dǎo)和求解見文獻(xiàn)[5]，此處不再贅述?；谔摴υ韀6]或矩陣函數(shù)法[4]進(jìn)行推導(dǎo)，可以得到二維情況下求解溫度場的控制方程：

其中[E0t]、[E1t]、[E2t]為單元邊界上的系數(shù)矩陣，上標(biāo)t代表溫度場，以便于與結(jié)構(gòu)分析中的相應(yīng)矩陣區(qū)分。該系數(shù)矩陣只和單元幾何構(gòu)型和材料參數(shù)有關(guān)式（2）～式（4），其計(jì)算和組裝與有限元法類似，式中k為熱傳導(dǎo)系數(shù)。

徑向內(nèi)部熱源為：

式（1）、式（5）可聯(lián)立表示成一階常微分方程：

式中，[Z]是哈密頓系數(shù)矩陣：

比例邊界有限元法的單元?jiǎng)偠染仃嚨慕M裝過程與有限元法相同，即先求得各個(gè)單元的剛度矩陣，而后將所有剛度陣進(jìn)行組裝即可得到整體剛度矩陣，將方程（8）代入式（6）中求得：

式中，ct是積分常數(shù)，可以由邊界溫度T（ξ=1）計(jì)算得到：

用形函數(shù)N（η）對T（ξ）插值，求出全域的溫度場T（ξ，η）：

1.2 位移場求解

比例邊界有限元法中的位移場在徑向可以獲得解析解，從相似中心O出發(fā)沿著徑向線，ξ處的位移為{u（ξ）}，徑向η是一個(gè)常量，因此只需要環(huán)向的插值形函數(shù)和{u（ξ）}即可表達(dá)域內(nèi)任意點(diǎn)的位移{u（ξ，η）}為：

應(yīng)力為：

式中，[D]是材料的彈性矩陣；B1（η）和B2（η）是應(yīng)變位移矩陣[7]。

用位移表達(dá)的比例邊界有限元方法的控制方程[8-9]，如下所示：

其中[E0]、[E1]、[E2]為單元邊界上的系數(shù)矩陣，只和單元的幾何構(gòu)型和材料參數(shù)有關(guān)，計(jì)算和組裝與傳統(tǒng)的有限元方法類似。{F（ξ）}為是外部荷載。若外部荷載只包含溫度荷載，則{F（ξ）}可表示為如下形式：

式中，d為計(jì)算溫度場時(shí)的特征值的相反數(shù)；| J（η）|為雅各比矩陣的行列式。

內(nèi)部等效節(jié)點(diǎn)力為：

式中，{q0（ξ）}為由初始溫度應(yīng)力產(chǎn)生的內(nèi)部節(jié)點(diǎn)力：

溫度引起的初始應(yīng)力{σ0}可以由下式推得[11]：

式中，T為式（11）求得的溫度場；{β}是與溫度相關(guān)的系數(shù)。

式中，v為泊松比；α為熱膨脹系數(shù)。

對式（12）進(jìn)行求解，可以求得位移為：

對位移u{ξ}進(jìn)行插值得到全域的位移{u（ξ，η）}，并通過位移法可以求得應(yīng)力強(qiáng)度因子[13]。

2 數(shù)值算例

為了驗(yàn)證模型的有效性，下文分析了帶裂縫平板的溫度場、位移場、應(yīng)力強(qiáng)度因子。

2.1 各向同性等溫裂縫板承受溫度荷載

考慮如圖1（a）所示的含中心裂縫的等溫裂縫平板，尺寸為2W×2W，中心裂縫長度2a，材料的彈性模量E=2.184×105Pa，泊松比ν=0.3，熱膨脹系數(shù)α=1.67×10-5/℃。圖中給出了板的溫度邊界條件，T代表溫度值，hq代表絕熱。為與文獻(xiàn)結(jié)果比較，采用平面應(yīng)變假定進(jìn)行分析?？紤]到結(jié)構(gòu)和荷載的對稱性，取1/2結(jié)構(gòu)進(jìn)行分析。整個(gè)模型離散為45個(gè)比例邊界有限單元（圖2）。

圖1 帶裂縫板（單位：℃）Fig. 1 Plate fracture

圖3給出了板的溫度云圖，其溫度值從裂縫處開始向結(jié)構(gòu)四周逐漸增大，在板頂部、底部與右側(cè)達(dá)到100℃，這一結(jié)果與參考文獻(xiàn)[12]采用擴(kuò)展無網(wǎng)格伽遼金法得到的結(jié)果相符。圖4為板的位移云圖，水平位移值從左側(cè)向右側(cè)兩角逐漸增大，在兩角處達(dá)到最大值為0.24cm；豎直方向位移值分別在板的左上角和左下角處達(dá)到最大和最小，其值分別為0.45cm與0.05cm。

圖2 網(wǎng)格圖Fig. 2 Polygon mesh

圖3 溫度場（單位：℃）Fig. 3 Temperature field

圖4 位移場（單位：cm）Fig. 4 Displacement field

2.2 各向同性絕熱裂縫板承受溫度荷載

考慮如圖1（b）所示的含中心裂縫的絕熱裂縫平板，板的尺寸是2W×2W，中心裂縫長度2a。材料的屬性同3.1節(jié)。圖中給出了板的溫度邊界條件，T代表溫度值，hq=0代表絕熱情況。由于結(jié)構(gòu)是對稱的，所以取半結(jié)構(gòu)進(jìn)行分析。采用平面應(yīng)變假定。整個(gè)模型離散為45個(gè)比例邊界有限單元（圖5）。

表1 等溫裂縫板的應(yīng)力強(qiáng)度因子對比（KI*）Tab. 1 Thermal insulation plate crack stress intensity factor of the contrast（KI*）

圖6給出了結(jié)構(gòu)的溫度云圖，其溫度值從板的低端向頂端逐漸增大，在板低端與頂端的溫度值分別為-100℃與100℃，并與擴(kuò)展無網(wǎng)格伽遼金法[12]（XEFG）所示溫度云圖相同；圖7給出了結(jié)構(gòu)的位移云圖，裂縫上部分板的水平位移從左下部向右上部逐漸增加，在右上角處達(dá)到最大值為0.20cm；下部分板的水平位移從左上部向右下部逐漸減小，在右下角處達(dá)到最小值為-0.20cm。板的左邊角處的豎向位移值最大，其值為0.1cm。

圖5 網(wǎng)格圖Fig. 5 Polygon mesh

圖6 溫度場（單位：℃）Fig. 6 Temperature field

圖7 位移場（單位：cm）Fig. 7 Displacement field

表2 絕熱裂縫板的應(yīng)力強(qiáng)度因子對比（KII*）Tab. 2 Thermal insulation plate crack stress intensity factor of the contrast（KII*）

2.3 各向同性等溫斜裂縫板承受溫度荷載

考慮如圖8所示的含中心裂縫的等溫斜裂縫平板，板的尺寸是2W×2L，θ1=0℃，θ2=100℃，W=1.0m，L=2.0m，a=0.3m，β=30o，材料參數(shù)為彈性模量E=2.184×105Pa，泊松比為0.3，熱膨脹系數(shù)為1.67×10-5/℃，采用平面應(yīng)變問題分析。

圖9給出了結(jié)構(gòu)的溫度云圖，其溫度值從裂縫處開始向結(jié)構(gòu)四周逐漸增大，裂縫處與四周的溫度值分別為0℃與100℃，并與擴(kuò)展無網(wǎng)格伽遼金法[12]（XEFG）中的溫度云圖相符。圖10給出了結(jié)構(gòu)的位移云圖，水平位移從板的左邊向后邊逐漸增大，在右邊端角達(dá)到最大值為0.4cm。豎直向位移從板的低邊向頂邊逐漸增大，在板底邊與定邊的豎直向位移值分別為0.05cm與0.75cm。本文方法計(jì)算的應(yīng)力強(qiáng)度因子KI*為0.296，KII*為0.047，擴(kuò)展無網(wǎng)格伽遼金法[14]（XEFG）所得的KI*為0.294，KII*為0.045，其相對誤差分別為0.6%、4.2%。

圖8 網(wǎng)格圖Fig. 8 Polygon mesh

圖9 溫度場（單位：℃）Fig. 9 Temperature field

圖10 位移場（單位：cm）Fig. 10 Displacement field

3 結(jié)束語

本文基于比例邊界有限元方法，提出了結(jié)構(gòu)在溫度荷載作用下的線彈性斷裂分析模型，得到了裂縫的擴(kuò)展過程和結(jié)構(gòu)的整體響應(yīng)。首先利用比例邊界有限元方法求解結(jié)構(gòu)的溫度場，在結(jié)構(gòu)分析中將溫度場轉(zhuǎn)化為與徑向坐標(biāo)有關(guān)的溫度荷載{F（ξ）}，代入結(jié)構(gòu)分析的比例邊界有限元方法求解格式，求出位移場、應(yīng)力場和應(yīng)力強(qiáng)度因子。通過三個(gè)承受溫度荷載的板的算例驗(yàn)證了本模型的計(jì)算精度，且可看出利用較少的網(wǎng)格即可獲得很高的精度。綜上，本文模型是求解結(jié)構(gòu)溫度斷裂分析的一種高精度高效率的數(shù)值方法。

[1] Tatro S，Schrader E. Thermal analysis for RCC—a practical approach[C]//Roller compacted concrete III. ASCE，1992：389-406.

[2] 程井，常曉林，周偉，等. 基于無單元伽遼金法的水工結(jié)構(gòu)溫度應(yīng)力及溫度裂紋擴(kuò)展計(jì)算[J]. 四川大學(xué)學(xué)報(bào)：工程科學(xué)版，2009（4）.CHENG Jing，CHANG Xiaolin，ZHOU Wer，et al. Simulation of thermal stress and crack propagation in Hydraulic structures based on EFGM[J]. Journal of Sichuan University: Applied Mechanics & Materials，2009 （4）.

[3] Duflot M. The extended finite element method in thermoelastic fracture mechanics[J]. International Journal for Numerical Methods in Engineering，2008，74（5）：827-847.

[4] Song C，Wolf J P. The scaled boundary finite-element method—alias consistent infinitesimal finite-element cell method—for elastodynamics[J]. Computer Methods in applied mechanics and engineering，1997，147（3）：329- 355.

[5] Song C. Analysis of singular stress fields at multi-material corners under thermal loading[J]. International journal for numerical methods in engineering，2006，65（5）：620-652.

[6] Deeks A J，Wolf J P. A virtual work derivation of the scaled boundary finite-element method for elastostatics[J]. Computational Mechanics，2002，28（6）：489-504.

[7] Song C. A matrix function solution for the scaled boundary finiteelement equation in statics[J]. Computer Methods in Applied Mechanics and Engineering，2004，193（23）：2325-2356.

[8] Song C，Wolf J P. Body loads in scaled boundary finiteelement method[J]. Computer methods in applied mechanics and engineering，1999，180（1）：117-135.

[9] Lin Z，Liao S. The scaled boundary FEM for nonlinear problems[J]. Communications in Nonlinear Science and Numerical Simulation，2011，16（1）：63-75.

[10] Li C，Ooi E T，Song C，et al. SBFEM for fracture analysis of piezoelectric composites under thermal load[J]. International Journal of Solids and Structures，2015，52：114-129.

[11] Song C，Tin-Loi F，Gao W. A definition and evaluation procedure of generalized stress intensity factors at cracks and multi-material wedges[J]. Engineering Fracture Mechanics，2010，77（12）：2316-2336.

[12] Bouhala L，Makradi A，Belouettar S. Thermal and thermomechanical influence on crack propagation using an extended mesh free method[J]. Engineering Fracture Mechanics，2012，88：35-48.

[13] Tsang D K L，Oyadiji S O，Leung A Y T. Two-dimensional fractal-like finite element method for thermoelastic crack analysis[J]. International Journal of Solids and Structures，2007，44（24）：7862-7876.

2017-03-23

2017-05-19

鐘 紅（1981—），女，副教授，主要研究方向：主要從事混凝土結(jié)構(gòu)靜動(dòng)力分析研究。E-mail：hzhong@dlut.edu.cn

賀 帥（1988—），男，碩士，主要研究方向：主要從事混凝土結(jié)構(gòu)溫度斷裂分析研究。E-mail：1255536943@qq.com

李紅軍（1988—），男，高級(jí)工程師，主要研究方向：主要從事大壩靜動(dòng)力分析研究。E-mail：lijunli1995@163.com

Thermal fracture analysis based on the Polygon Scaled boundary finite element method

ZHONG Hong1，HE Shuai1，LI Hongjun2
（1.Faculty of Infrastructure Engineering，Dalian 116024，China；2.China Institute of Water Resources and Hydropower Research，Beijing 100044，China）

Thermal load can lead to cracking of structure and thus seriously influence the stability of the structure. A numerical model for the thermal fracture analysis is presented based on the Scaled Boundary Finite Element Method（SBFEM）. The temperature distribution is solved via the SBFEM，then the corresponding thermal load is divided into several items in the form of power functions，and the structural response is further obtained based on the principle of linear superposition. In the presented model，only the boundary of the structure is discretized，the temperature，displacement and stress fields can be obtained semi-analytically，thus the stress intensity factors are solved with high precision.Through multiple plate under temperature load，temperature field analysis and structure analysis can be seen that using fewer grid can obtain higher precision.Simulation results of stress intensity factor as the growth of the crack length increases，and the stress intensity factor consistent with the existing literature，which demonstrates the effectiveness of the model.

thermal load； the Scaled boundary finite element method； temperature fields； stress intensity factor

TV313

A學(xué)科代碼：130.1545

10.3969/j.issn.2096-093X.2017.03.013

國家自然科學(xué)基金（51579033）；國家重點(diǎn)研發(fā)計(jì)劃（2016YFB0201000，2016YFC0401600）。