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

鐘 紅，賀 帥，李紅軍

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

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

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

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

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

溫度荷載；比例邊界有限元；溫度場；應力強度因子

0 引言

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

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

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

1 基本原理

1.1 溫度場求解

本文基于比例邊界有限元法（Scaled Boundary Finite Element Method，SBFEM）開展溫度場分析和結構分析。該方法的控制方程的推導和求解見文獻[5]，此處不再贅述。基于虛功原理[6]或矩陣函數法[4]進行推導，可以得到二維情況下求解溫度場的控制方程：

其中[E0t]、[E1t]、[E2t]為單元邊界上的系數矩陣，上標t代表溫度場，以便于與結構分析中的相應矩陣區分。該系數矩陣只和單元幾何構型和材料參數有關式（2）～式（4），其計算和組裝與有限元法類似，式中k為熱傳導系數。

徑向內部熱源為：

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

式中，[Z]是哈密頓系數矩陣：

比例邊界有限元法的單元剛度矩陣的組裝過程與有限元法相同，即先求得各個單元的剛度矩陣，而后將所有剛度陣進行組裝即可得到整體剛度矩陣，將方程（8）代入式（6）中求得：

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

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

1.2 位移場求解

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

應力為：

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

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

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

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

內部等效節點力為：

式中，{q0（ξ）}為由初始溫度應力產生的內部節點力：

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

式中，T為式（11）求得的溫度場；{β}是與溫度相關的系數。

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

對式（12）進行求解，可以求得位移為：

對位移u{ξ}進行插值得到全域的位移{u（ξ，η）}，并通過位移法可以求得應力強度因子[13]。

2 數值算例

為了驗證模型的有效性，下文分析了帶裂縫平板的溫度場、位移場、應力強度因子。

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

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

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

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

圖2 網格圖Fig. 2 Polygon mesh

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

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

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

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

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

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

圖5 網格圖Fig. 5 Polygon mesh

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

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

表2 絕熱裂縫板的應力強度因子對比（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，材料參數為彈性模量E=2.184×105Pa，泊松比為0.3，熱膨脹系數為1.67×10-5/℃，采用平面應變問題分析。

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

圖8 網格圖Fig. 8 Polygon mesh

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

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

3 結束語

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

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

[2] 程井，常曉林，周偉，等. 基于無單元伽遼金法的水工結構溫度應力及溫度裂紋擴展計算[J]. 四川大學學報：工程科學版，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—），女，副教授，主要研究方向：主要從事混凝土結構靜動力分析研究。E-mail：hzhong@dlut.edu.cn

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

李紅軍（1988—），男，高級工程師，主要研究方向：主要從事大壩靜動力分析研究。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學科代碼：130.1545

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

國家自然科學基金（51579033）；國家重點研發計劃（2016YFB0201000，2016YFC0401600）。