李星 顧鑫 夏曉舟 陳愛玖 章青,

*(河海大學(xué)力學(xué)與材料學(xué)院,南京 210098)

?(華北水利水電大學(xué)土木與交通學(xué)院,鄭州 450045)

引言

多孔介質(zhì)材料組成的結(jié)構(gòu)服役于寒冷潮濕的環(huán)境時(shí),其內(nèi)部會(huì)發(fā)生傳熱傳質(zhì)和相變現(xiàn)象.多孔介質(zhì)中孔隙水的凍結(jié)將導(dǎo)致材料變形,誘發(fā)材料的損傷甚至破壞失效,如寒冷地區(qū)混凝土結(jié)構(gòu)的凍融破壞[1-2]和巖土體的凍脹破壞[3]均與材料內(nèi)部水的凍結(jié)有關(guān).孔隙水結(jié)冰時(shí)體積膨脹,將導(dǎo)致冰與孔隙壁發(fā)生相互作用,與此同時(shí),相變潛熱的釋放也將影響結(jié)構(gòu)的溫度場(chǎng)和變形場(chǎng),是一個(gè)復(fù)雜的熱-力耦合問題.發(fā)展能夠準(zhǔn)確刻畫包含材料相變過程的熱-力耦合模型及其數(shù)值計(jì)算方法,模擬真實(shí)的水結(jié)冰過程和其中的熱-力耦合效應(yīng),對(duì)于揭示含水材料的凍脹破壞機(jī)理具有重要的科學(xué)意義.

在進(jìn)行熱-力耦合問題模擬時(shí),準(zhǔn)確求解溫度場(chǎng)是獲得結(jié)構(gòu)變形場(chǎng)的基礎(chǔ).材料的相變問題具有強(qiáng)非線性(不同相之間物理性質(zhì)的快速變化)和不連續(xù)(不同相之間邊界、裂紋)的特征,在進(jìn)行固-液相變的熱傳導(dǎo)數(shù)值模擬時(shí),需要將固-液界面的演化作為解的一部分進(jìn)行確定,亦稱為“Stefan”問題,相應(yīng)的數(shù)值求解方法主要有界面跟蹤法和固定網(wǎng)格法[4-6].界面跟蹤法使用移動(dòng)網(wǎng)格數(shù)值方案來捕獲固-液移動(dòng)邊界,將整個(gè)求解域在相變界面處分為固相域和液相域分開求解,該方法需要在每一個(gè)時(shí)間步長下確定固-液兩相界面的位置且在相變界面上施加“Stefan”條件,求解不易實(shí)施,不適用于復(fù)雜的相變界面.固定網(wǎng)格法使用固定網(wǎng)格方案,將求解域中固相和液相編寫入一組控制方程進(jìn)行求解,此類方法通過對(duì)控制方程中傳輸量的躍變來考慮相變過程,不需要跟蹤固-液兩相界面的位置,代表性方法有熱焓法模型[7-8]、顯熱容模型[9-10]和等效熱容模型[10]等.

連續(xù)介質(zhì)力學(xué)理論建立了偏微分型的熱-力耦合方程,但在求解含孔洞、演化裂紋等強(qiáng)非線性和不連續(xù)問題時(shí)面臨瓶頸.近場(chǎng)動(dòng)力學(xué)(peridynamics,PD)采用空間積分方程描述力學(xué)響應(yīng),在數(shù)值求解時(shí)不存在空間求導(dǎo),可以有效避免傳統(tǒng)連續(xù)性模型求解時(shí)出現(xiàn)的病態(tài)特征[11-12],為熱致變形破壞等問題的數(shù)值模擬提供了新思路[13].

起初,近場(chǎng)動(dòng)力學(xué)主要用于研究機(jī)械載荷作用下材料和結(jié)構(gòu)的損傷破壞問題,隨后被用于模擬各種擴(kuò)散問題.Gerstle等[14]首次在近場(chǎng)動(dòng)力學(xué)框架內(nèi)建立了非局部擴(kuò)散模型,研究了一維物體中的熱傳導(dǎo)過程,為擴(kuò)散問題的PD 建模奠定了基礎(chǔ).此后,PD 理論逐步被應(yīng)用到熱傳導(dǎo)領(lǐng)域中,研究者先后構(gòu)建了鍵型 PD 多維非穩(wěn)態(tài)熱傳導(dǎo)模型[15-16]與態(tài)型PD熱傳導(dǎo)模型[17-19],并對(duì)多維結(jié)構(gòu)、多種邊界條件、含缺陷和裂紋的熱傳導(dǎo)問題進(jìn)行了模擬,證明了PD理論在求解熱傳導(dǎo)問題時(shí)的有效性.但上述PD 熱傳導(dǎo)模型都未涉及相變問題的描述,Madenci等[20]用近場(chǎng)動(dòng)力學(xué)微分算子將偏微分型的熱焓方程改寫為積分型方程,模擬了物質(zhì)的一維相變過程,模擬結(jié)果和解析解吻合較好.Nikolaev等[21-22]率先建立了包含材料相變的鍵型PD 熱傳導(dǎo)模型,求解了均質(zhì)物質(zhì)一維凍結(jié)和二維對(duì)流傳熱相變問題,對(duì)比解析解驗(yàn)證了模型的正確性,為鍵型PD 熱傳導(dǎo)模型中考慮相變提供了理論參考.

近年來,PD 理論也被廣泛應(yīng)用到熱致變形及熱致?lián)p傷破壞問題的研究中[23-24].國內(nèi)外學(xué)者建立了多種近場(chǎng)動(dòng)力學(xué)熱-力耦合模型,通過對(duì)大量的熱傳導(dǎo)問題的分析,證明其在分析熱致?lián)p傷及破壞問題時(shí)的可行性.目前建立的PD 熱-力耦合模型主要包括非耦合的模型和全耦合的模型兩類.前者只考慮溫度對(duì)結(jié)構(gòu)變形的影響而忽略結(jié)構(gòu)變形對(duì)溫度的影響,通過在PD 運(yùn)動(dòng)方程中的本構(gòu)力函數(shù)或力矢量狀態(tài)中引入溫度影響項(xiàng)構(gòu)建PD 熱-力耦合模型[25-28]后者在熱傳導(dǎo)方程中加入結(jié)構(gòu)變形加熱和冷卻項(xiàng)來反映結(jié)構(gòu)變形對(duì)溫度場(chǎng)的影響,實(shí)現(xiàn)了溫度場(chǎng)和變形場(chǎng)的雙向耦合[29-31].需要指出的是,目前已建立的PD 熱-力耦合模型在進(jìn)行溫度場(chǎng)的求解時(shí)還未涉及材料相變影響,材料相變時(shí)潛熱的吸收與釋放會(huì)顯著影響溫度場(chǎng)分布,進(jìn)而影響結(jié)構(gòu)的變形.據(jù)此,本文在近場(chǎng)動(dòng)力學(xué)框架內(nèi)發(fā)展了考慮相變的鍵型PD熱-力耦合模型及其數(shù)值求解方法,模擬了飽水多孔介質(zhì)在降溫條件下的凍結(jié)破壞過程,以期建立合理可靠的多孔介質(zhì)凍結(jié)破壞的數(shù)值模擬方法.

圖1 物質(zhì)點(diǎn)間的相互作用Fig.1 Interaction between material points xpand in PD theory

1 考慮相變的近場(chǎng)動(dòng)力學(xué)熱-力耦合方程

1.1 鍵型近場(chǎng)動(dòng)力學(xué)運(yùn)動(dòng)方程

材料的微彈性模量c可通過PD和經(jīng)典連續(xù)介質(zhì)力學(xué)中的應(yīng)變能密度等效獲得,具體為

式中,E為材料彈性模量,ν 是泊松比,h是平面問題中板的厚度,A是一維問題中桿件的截面面積.

定義物質(zhì)點(diǎn)對(duì)的伸長率達(dá)到某一臨界伸長率s0時(shí),該物質(zhì)點(diǎn)對(duì)發(fā)生斷裂,物質(zhì)點(diǎn)之間不再有力的作用,用間斷函數(shù)來標(biāo)記物質(zhì)點(diǎn)對(duì)的斷裂狀態(tài),則有

臨界伸長率s0與材料的斷裂能釋放率有關(guān),表示為

式中,GF為斷裂能釋放率.因此,材料的損傷可以通過材料點(diǎn)的斷裂鍵數(shù)與總鍵數(shù)之比來確定

1.2 含相變的鍵型近場(chǎng)動(dòng)力學(xué)的熱傳導(dǎo)方程

在使用固定網(wǎng)格方案進(jìn)行相變問題模擬時(shí),熱焓法被廣泛采用.熱焓法的主要思想是將溫度、比熱容和潛熱合并為熱傳導(dǎo)方程中的熱焓項(xiàng),在整個(gè)求解域中建立統(tǒng)一的能量方程.通過對(duì)熱焓求解即可確定兩相界面,因而不需要追蹤相變界面,可用來求解復(fù)雜邊界條件下的相變問題[4].具體地,引入熱焓H(T,x,t),將基于傅里葉定律的熱傳導(dǎo)方程重新寫成以下形式[20]

式中,H(T,x,t)是t時(shí)刻x處的熱焓,K是材料熱傳導(dǎo)系數(shù),T(x,t)是t時(shí)刻x處的溫度.熱焓和導(dǎo)熱系數(shù)分別為

式中,C1和C2分別是液相和固相的比熱容,K1和K2分別是液相和固相的熱傳導(dǎo)系數(shù),L是材料相變潛熱,Tmelt是材料相變溫度,G(T) 是Heaviside 階躍函數(shù).

由式(8)求溫度可得

基于此,可在近場(chǎng)動(dòng)力學(xué)框架內(nèi)建立熱焓形式的熱傳導(dǎo)模型,進(jìn)行相變問題的模擬.采用近場(chǎng)動(dòng)力學(xué)方法求解熱傳導(dǎo)問題不僅考慮了非局部性,而且允許不連續(xù)演化.如圖2 所示,在PD 熱傳導(dǎo)模型中,質(zhì)點(diǎn)間的非局部相互作用是伴隨熱能交換產(chǎn)生的.一個(gè)質(zhì)點(diǎn)與其鄰域內(nèi)的其他質(zhì)點(diǎn)發(fā)生熱量交換,熱量交換通過兩個(gè)質(zhì)點(diǎn)之間的熱鍵(T-bond)完成,熱鍵的相關(guān)定義與運(yùn)動(dòng)方程中物質(zhì)點(diǎn)之間鍵的定義相同.Nikolaev等[21]給出了鍵型PD 中用熱焓表達(dá)的熱傳導(dǎo)方程,表示為

圖2 PD 模型中熱鍵Fig.2 T-bonds in PD model

式中,kT是近場(chǎng)動(dòng)力學(xué)中的微熱導(dǎo)率,其與經(jīng)典力學(xué)理論中的導(dǎo)熱系數(shù)K相關(guān),表示為

式中,K是材料的熱傳導(dǎo)系數(shù).

如圖2 所示,在近場(chǎng)動(dòng)力學(xué)模擬中,為體現(xiàn)不同相之間的材料屬性差異,根據(jù)熱鍵兩端連接的物質(zhì)點(diǎn)種類的不同,把熱鍵分為液相熱鍵、固相熱鍵和液-固界面熱鍵,計(jì)算時(shí)熱鍵的熱力學(xué)參數(shù)根據(jù)熱鍵的種類進(jìn)行賦值,并根據(jù)相變狀態(tài)實(shí)時(shí)更新.

綜上所述,結(jié)合式(1)和式(11)即可得到考慮相變的鍵型PD 熱-力耦合模型,第一個(gè)方程是滿足熱彈性本構(gòu)關(guān)系的運(yùn)動(dòng)方程,第二個(gè)方程是熱焓形式的熱傳導(dǎo)方程.需要指出的是,在該耦合模型中只考慮溫度對(duì)變形的影響而忽略了結(jié)構(gòu)變形對(duì)溫度的影響,屬于單向耦合的形式.

2 熱-力耦合模型的數(shù)值求解方法

對(duì)于耦合方程組的數(shù)值求解,一般采用單步法或交錯(cuò)法進(jìn)行.因本文熱-力耦合方程為單向耦合的形式,此外為避免單步法中采用小時(shí)間增量而造成計(jì)算量過大的問題,采用交錯(cuò)算法對(duì)耦合方程組進(jìn)行數(shù)值求解.其中,結(jié)構(gòu)變形場(chǎng)和溫度場(chǎng)自然分開求解,即運(yùn)動(dòng)方程用于求解變形場(chǎng),熱傳導(dǎo)方程用于求解溫度場(chǎng).

數(shù)值求解時(shí),運(yùn)動(dòng)方程和熱傳導(dǎo)方程采用顯式差分方案進(jìn)行求解,分別為

式中,n是時(shí)間步編號(hào),i和j是物質(zhì)點(diǎn)編號(hào),ΔtM和ΔtT分別是運(yùn)動(dòng)方程和熱傳導(dǎo)方程的時(shí)間增量步.

顯式積分算法雖然簡(jiǎn)單明了,但其是條件穩(wěn)定的,求解時(shí)采用的時(shí)間增量步需要滿足穩(wěn)定性條件.可以根據(jù)von Neumann 穩(wěn)定性分析方法獲得穩(wěn)定性時(shí)間增量步ΔtM和ΔtT[33].采用交錯(cuò)算法時(shí)先進(jìn)行溫度場(chǎng)的求解,將當(dāng)前的溫度狀態(tài)看作是一個(gè)“假定穩(wěn)態(tài)”,然后進(jìn)行變形場(chǎng)的循環(huán)計(jì)算,直至此時(shí)的變形場(chǎng)達(dá)到穩(wěn)定狀態(tài).當(dāng)ΔtT和ΔtM差異很大,甚至前者高出后者幾個(gè)數(shù)量級(jí)時(shí),采用真實(shí)的ΔtM來進(jìn)行變形場(chǎng)的計(jì)算會(huì)導(dǎo)致較高的計(jì)算開銷.為提高求解效率使變形場(chǎng)快速達(dá)到穩(wěn)定狀態(tài),本文運(yùn)動(dòng)方程采用自適應(yīng)動(dòng)態(tài)松弛法(adaptive dynamic relaxation,ADR)進(jìn)行求解.

在ADR 方法中,系統(tǒng)所有質(zhì)點(diǎn)的運(yùn)動(dòng)方程通過引入虛擬慣性力和阻尼項(xiàng),寫成一系列常微分方程[33]

式中,D為虛擬對(duì)角密度矩陣,cM為阻尼系數(shù),分別通過Greschgorin 定理和Rayleigh 得到.F(u,u′,x,x′,t)是由PD 相互作用力和體力構(gòu)成的力密度矢量.下一時(shí)間步的位移和速度可使用中心差分的顯式積分方法得到

在動(dòng)態(tài)松弛法中,力矢量F是唯一的物理量,密度矩陣D、阻尼系數(shù)cM和時(shí)間步長Δt不需要是物理量.因此可以選擇合適的參數(shù)在滿足穩(wěn)定性要求的前提下,實(shí)現(xiàn)快速的收斂.

圖3 是熱-力耦合模型數(shù)值計(jì)算流程圖.在進(jìn)行交錯(cuò)求解時(shí),對(duì)得到的位移場(chǎng)是否達(dá)到穩(wěn)定的判斷直接影響求解效率和求解精度,是程序?qū)嵤┑年P(guān)鍵.本文通過物質(zhì)點(diǎn)在當(dāng)前構(gòu)型和前一個(gè)時(shí)間步構(gòu)型中位移差值的均方根值Er和一個(gè)足夠小量 Ω 的比較來判斷當(dāng)前位移場(chǎng)是否達(dá)到穩(wěn)態(tài).當(dāng)Er<Ω 時(shí),認(rèn)為位移場(chǎng)達(dá)到了穩(wěn)定狀態(tài),進(jìn)而進(jìn)行下一步溫度場(chǎng)的求解;否則,繼續(xù)進(jìn)行位移場(chǎng)的求解,直至滿足穩(wěn)定性條件.

圖3 熱-力耦合模型數(shù)值計(jì)算程序流程圖Fig.3 Flow chart of numerical program of thermomechanical coupled model

3 數(shù)值算例驗(yàn)證與分析

3.1 方板角凍結(jié)過程模擬

由單種相變材料組成的邊長為1.5 m 的正方形板,處于液相狀態(tài),初始溫度T0=1.3 K,材料相變溫度Tmelt=1.0K.板的左邊界和下邊界施加恒定的溫度邊界Ti=0K,上邊界和右邊界保持絕熱.在邊界溫度作用下,方板從邊界處開始降溫并發(fā)生從左下角向右上角的凍結(jié).為驗(yàn)證前述考慮相變的PD 熱傳導(dǎo)模型的準(zhǔn)確性,模擬了相變潛熱分別為L1=0.25 J/kg和L2=1.0J/kg 時(shí)方板的凍結(jié)過程.

在近場(chǎng)動(dòng)力學(xué)模型中,方板被均勻離散,物質(zhì)點(diǎn)間距Δx=5 mm,近場(chǎng)范圍取δ=3Δx.材料熱力學(xué)參數(shù)為,比熱容Cs=Cl=1.0J/(kg·K),密度ρs=ρl=1.0kg/m3,熱傳導(dǎo)系數(shù)Ks=Kl=1.0W/(m·K),其中下標(biāo)S和l分別表示固相和液相.

圖4 是模擬得到的潛熱L1=0.25 J/kg 的方板在不同時(shí)刻的溫度場(chǎng)分布和相變界面位置.從圖4(a)中可以看出,模型在邊界溫度作用下開始降溫,當(dāng)材料溫度低于相變溫度時(shí)發(fā)生凍結(jié),凍結(jié)過程由左下角向右上角發(fā)展.在圖4(b)中,藍(lán)色表示液相,綠色表示固相,紅色表示相變的界面.從圖中可以看出,方板在凍結(jié)過程相變界面自左下角開始向右上角移動(dòng),符合方板凍結(jié)的客觀規(guī)律,說明PD 方法能夠追蹤到凍結(jié)過程中相變界面位置,可以有效模擬出方板的角凍結(jié)過程.

圖4 方板角凍結(jié)PD 模擬結(jié)果Fig.4 PD simulation results of square plate corner freezing

圖5 給出了時(shí)間為0.1 s 時(shí)相變界面位置,并列入了Kova?evi?等[34]采用有限體積法FVM和基函數(shù)配點(diǎn)法RBFCM 得到的結(jié)果.可以明顯看出,在采用不同的相變潛熱L1=0.25 J/kg和L2=1.0J/kg 時(shí),PD 方法捕捉到的相變界面位置和文獻(xiàn)[34]用兩種方法模擬得到的結(jié)果均吻合較好,充分驗(yàn)證了本文給出的含相變的PD 熱傳導(dǎo)模型在求解含相變的熱傳導(dǎo)問題時(shí)的有效性.

圖5 相變界面位置比較Fig.5 Comparison of phase change interface position

3.2 方板熱致變形模擬

如圖6 所示,假設(shè)正方形板邊長為1 cm,初始溫度為T0=10°C,上邊和右邊受到降溫載荷Ti=1°C的作用,左邊和下邊絕熱并約束法向位移.將該模型作為平面應(yīng)變問題,分別利用PD 熱-力耦合模型和有限元軟件ABAQUS 模擬了該方板在降溫時(shí)的力學(xué)響應(yīng).

圖6 受到雙面溫度載荷作用的正方形板Fig.6 Square plate subjected to double side temperature load

在近場(chǎng)動(dòng)力學(xué)模擬中,方板被均勻離散,物質(zhì)點(diǎn)間距Δx=0.1mm,近場(chǎng)范圍取δ=3Δx.材料熱力學(xué)參數(shù)如下[35],彈性模量E=25.7 GPa,泊松比ν=0.25,熱傳導(dǎo)系數(shù)K=1.0J/(s·m·K),密度ρ=1800kg/m3,比熱容C=1688J/(kg·K),熱膨脹系數(shù)α=1.8×10-5K-1.

圖7 給出了近場(chǎng)動(dòng)力學(xué)和有限元模擬得到的板在t=200s 時(shí)的溫度場(chǎng)和位移場(chǎng)云圖.可以看出,在相同的時(shí)間下,兩種方法得到的溫度場(chǎng)和位移場(chǎng)基本一致.圖8 給出了兩種方法模擬得到的方板中A,B,C三個(gè)位置處的溫度和水平方向的位移的時(shí)間歷程曲線,隨著時(shí)間的增加,方板溫度由右上角向左下角逐漸降低,溫度降低引起材料收縮,A和B處的水平位移為負(fù)值,這與熱脹冷縮的實(shí)際情況相符.模擬表明,無論是在同一溫度下還是隨著溫度的變化,PD 模擬結(jié)果和有限元結(jié)果都吻合較好,驗(yàn)證了PD 熱-力耦合模型在求解熱力耦合問題時(shí)的準(zhǔn)確性.

圖7 有限元(上)和近場(chǎng)動(dòng)力學(xué)(下)模擬的t=200s 時(shí)結(jié)果比較Fig.7 Simulation result comparison of finite element and periynamics at t=200s

圖8 近場(chǎng)動(dòng)力學(xué)和有限元模擬得到的A,B,C 點(diǎn)溫度和位移時(shí)間歷程比較Fig.8 Comparison of temperature and displacement histories of points A,B and C obtained from peridynamics and finite element simulation

3.3 多孔介質(zhì)的低溫凍結(jié)破壞模擬

基于前述建立的熱力耦合模型及其數(shù)值計(jì)算方法,建立如下的計(jì)算模型來模擬多孔介質(zhì)的凍結(jié)破壞過程.混凝土是典型的準(zhǔn)脆性材料且其凍融破壞比較常見,因此選用混凝土材料進(jìn)行說明.雖然導(dǎo)致混凝土凍融破壞的因素很多,眾多學(xué)者也先后提出了多個(gè)理論模型來揭示混凝土凍融破壞的機(jī)理,但因孔隙水凍結(jié)體積增大而引起內(nèi)應(yīng)力從而誘發(fā)破壞已是學(xué)術(shù)界的共識(shí).水結(jié)冰后體積增大9%,在混凝土孔隙內(nèi)產(chǎn)生結(jié)晶壓力和靜水壓力.由此產(chǎn)生的內(nèi)應(yīng)力超過混凝土的材料強(qiáng)度時(shí),將誘發(fā)混凝土材料的損傷及微裂紋的萌生擴(kuò)展,經(jīng)過多次凍融循環(huán)后,最終導(dǎo)致混凝土結(jié)構(gòu)的破壞.

如圖9 所示,假設(shè)含有圓形孔隙的砂漿方板,板邊長300μm,開始時(shí)孔隙內(nèi)充滿液態(tài)水,初始溫度為T0=12°C,結(jié)冰溫度Tmelt=0°C.方板上邊界和右邊界施加溫度為Ti=-30°C 的降溫條件,板左邊界和下邊界絕熱并約束法向位移.為反映孔隙水凍結(jié)時(shí)體積膨脹效應(yīng),當(dāng)孔隙凍結(jié)時(shí)在孔隙邊界上施加沿著徑向的位移載荷,其大小與孔隙半徑相關(guān),按照體積增大9%確定,位移載荷的施加順序與孔隙凍結(jié)的順序保持一致.

圖9 含孔隙的砂漿方板Fig.9 Mortar square plate with pores

在近場(chǎng)動(dòng)力學(xué)模擬中,板被均勻離散,物質(zhì)點(diǎn)間距Δx=1 μm,近場(chǎng)范圍取 δ=3Δx.砂漿材料參數(shù)和3.2 節(jié)的材料一致,水和冰的參數(shù)為:熱傳導(dǎo)系數(shù)Kl=0.6J/(s·m·K),Ks=2.14J/(s·m·K),比熱容Cl=4182J/(kg·K),Cs=2060J/(kg·K),密度ρl=1000kg/m3,ρs=920kg/m3,相變潛熱L=334 kJ/kg,其中下標(biāo)l和S分別表示水和冰.

圖10為模擬得到的砂漿方板凍結(jié)破壞的結(jié)果.如圖10(a)所示,由于孔隙水結(jié)冰時(shí)釋放潛熱,孔隙和周圍砂漿基質(zhì)的溫度變化滯后于其他區(qū)域,板內(nèi)溫度場(chǎng)的分布明顯受到影響而呈非均勻分布.圖10(b)為方板中孔隙水的結(jié)冰過程,其中藍(lán)色表示砂漿基質(zhì),青色表示液態(tài)水,黃色表示冰,紅色表示冰水界面.圖10(c)為方板凍結(jié)過程中的損傷圖.當(dāng)t=0.006 s 時(shí),靠近方板邊界處的孔先開始結(jié)冰.水結(jié)冰的體積膨脹促使孔隙邊界向外擴(kuò)張,在較近的孔隙1和2 之間引起較高應(yīng)力集中誘發(fā)砂漿基質(zhì)的損傷,如圖10(c)所示.隨著板內(nèi)溫度的持續(xù)降低,更多的孔隙開始結(jié)冰.當(dāng)t=0.105 s 孔隙全部凍結(jié),板內(nèi)的損傷由孔隙周邊的點(diǎn)狀分布匯聚成裂紋并在基質(zhì)中擴(kuò)展.從最終的損傷分布可以看出,凍結(jié)過程中裂紋的萌生與擴(kuò)展主要集中分布在孔隙邊界上和孔隙之間.孔隙邊界上裂紋的環(huán)向擴(kuò)展引起孔隙邊界處基質(zhì)材料剝落,是引起孔隙半徑增大的主要原因,印證了在凍融實(shí)驗(yàn)中得出的孔隙半徑隨著凍融次數(shù)增加而增大的結(jié)論[1].此外,孔隙半徑的增加致使孔隙內(nèi)充填更多的水,會(huì)加劇凍融破壞的發(fā)生.孔隙之間的裂紋擴(kuò)展導(dǎo)致材料內(nèi)部孔隙連通性增加、滲透性提高,易在凍融循環(huán)條件下形成更多更大的宏觀裂紋,是引起結(jié)構(gòu)凍脹開裂、表面剝蝕和宏觀力學(xué)性能下降的主要原因.

圖10 砂漿方板的凍結(jié)破壞過程Fig.10 Frost failure process of mortar plate

通過以上分析可知,本文建立的考慮相變的鍵型PD 熱-力耦合模型及其數(shù)值求解方法能夠模擬混凝土材料凍結(jié)破壞的全過程,真實(shí)再現(xiàn)凍結(jié)破壞現(xiàn)象.在一定程度上可以解釋混凝土材料產(chǎn)生凍融破壞的原因,為進(jìn)一步研究和揭示混凝土材料凍融破壞機(jī)理提供參考.

4 結(jié)論

本文建立了考慮材料相變的近場(chǎng)動(dòng)力學(xué)熱-力耦合模型及其數(shù)值求解方法,并對(duì)其有效性進(jìn)行了驗(yàn)證,取得的主要結(jié)論如下:

基于熱焓表達(dá)的近場(chǎng)動(dòng)力學(xué)熱傳導(dǎo)模型,可以較好地模擬包含材料相變的熱傳導(dǎo)過程,能夠準(zhǔn)確預(yù)測(cè)出相變界面位置和溫度場(chǎng)的分布,為求解相變傳熱問題提供分析方法.

模擬了以砂漿為代表的含水多孔介質(zhì)的凍結(jié)破壞過程,再現(xiàn)凍結(jié)破壞現(xiàn)象,為今后繼續(xù)研究并揭示此類材料的凍結(jié)破壞機(jī)理提供了分析思路.