考慮相變的近場動力學熱-力耦合模型及多孔介質凍結破壞模擬1)

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

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

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

引言

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

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

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

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

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

圖1 物質點間的相互作用Fig.1 Interaction between material points xpand in PD theory

1 考慮相變的近場動力學熱-力耦合方程

1.1 鍵型近場動力學運動方程

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

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

定義物質點對的伸長率達到某一臨界伸長率s0時,該物質點對發(fā)生斷裂,物質點之間不再有力的作用,用間斷函數來標記物質點對的斷裂狀態(tài),則有

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

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

1.2 含相變的鍵型近場動力學的熱傳導方程

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

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

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

由式(8)求溫度可得

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

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

式中,kT是近場動力學中的微熱導率,其與經典力學理論中的導熱系數K相關,表示為

式中,K是材料的熱傳導系數.

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

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

2 熱-力耦合模型的數值求解方法

對于耦合方程組的數值求解,一般采用單步法或交錯法進行.因本文熱-力耦合方程為單向耦合的形式,此外為避免單步法中采用小時間增量而造成計算量過大的問題,采用交錯算法對耦合方程組進行數值求解.其中,結構變形場和溫度場自然分開求解,即運動方程用于求解變形場,熱傳導方程用于求解溫度場.

數值求解時,運動方程和熱傳導方程采用顯式差分方案進行求解,分別為

式中,n是時間步編號,i和j是物質點編號,ΔtM和ΔtT分別是運動方程和熱傳導方程的時間增量步.

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

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

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

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

圖3 是熱-力耦合模型數值計算流程圖.在進行交錯求解時,對得到的位移場是否達到穩(wěn)定的判斷直接影響求解效率和求解精度,是程序實施的關鍵.本文通過物質點在當前構型和前一個時間步構型中位移差值的均方根值Er和一個足夠小量 Ω 的比較來判斷當前位移場是否達到穩(wěn)態(tài).當Er<Ω 時,認為位移場達到了穩(wěn)定狀態(tài),進而進行下一步溫度場的求解;否則,繼續(xù)進行位移場的求解,直至滿足穩(wěn)定性條件.

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

3 數值算例驗證與分析

3.1 方板角凍結過程模擬

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

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

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

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

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

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

3.2 方板熱致變形模擬

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

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

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

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

圖7 有限元(上)和近場動力學(下)模擬的t=200s 時結果比較Fig.7 Simulation result comparison of finite element and periynamics at t=200s

圖8 近場動力學和有限元模擬得到的A,B,C 點溫度和位移時間歷程比較Fig.8 Comparison of temperature and displacement histories of points A,B and C obtained from peridynamics and finite element simulation

3.3 多孔介質的低溫凍結破壞模擬

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

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

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

在近場動力學模擬中,板被均勻離散,物質點間距Δx=1 μm,近場范圍取 δ=3Δx.砂漿材料參數和3.2 節(jié)的材料一致,水和冰的參數為:熱傳導系數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,其中下標l和S分別表示水和冰.

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

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

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

4 結論

本文建立了考慮材料相變的近場動力學熱-力耦合模型及其數值求解方法,并對其有效性進行了驗證,取得的主要結論如下:

基于熱焓表達的近場動力學熱傳導模型,可以較好地模擬包含材料相變的熱傳導過程,能夠準確預測出相變界面位置和溫度場的分布,為求解相變傳熱問題提供分析方法.

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