瀝青混合料鋪筑溫度場的有限差分仿真模型

王黎明， 李 崟

(東北林業(yè)大學土木工程學院,哈爾濱 150040)

熱態(tài)瀝青混合料的溫度嚴重影響著壓實效果,碾壓質(zhì)量差會引致路面的早期損壞。這就需要道路工作者清楚在特定條件下混合料的降溫規(guī)律,以便合理安排施工。而熱瀝青混合料冷卻過程復(fù)雜,實地觀測不能得出鋪筑熱擴散的全面規(guī)律,而建立仿真模型是解決問題的關(guān)鍵。

20世紀60年代,Corlew等[1]基于有限差分的顯示格式建立了路面鋪筑溫度場的模型,但未考慮鋪筑過程的層厚和材料熱物性參數(shù)的變化,這使得該模型和服役期路面溫度場無本質(zhì)區(qū)別；1998年,美國明尼蘇達州編制了可預(yù)估瀝青混合料鋪筑降溫規(guī)律的“pavecool”程序[2]；但其未考慮對流換熱、路表面灑水換熱和鋪層比熱容的變化。中國雖然起步較晚,但眾學者對瀝青路面溫度場的研究頗為翔實。孫潔[3]基于通用CAE軟件,對鋪筑階段降溫規(guī)律進行了系統(tǒng)研究,但其研究缺少鋪層材料熱物性以及天空陰晴狀態(tài)對降溫規(guī)律影響的考量；呂得保[4]基于實測得出了瀝青混合料鋪筑過程的溫度回歸模型,但由于觀測因素的有限性,該模型的可信性有待商榷；王黎明等[5-6]運用Ansys軟件對熱態(tài)瀝青混合料鋪筑溫度場進行了系統(tǒng)研究,對材料熱物性和邊界條件在鋪筑過程的變化進行了細致研究,可實現(xiàn)高可信度的仿真,但該研究尚基于大型通用CAE軟件，因此難以推廣。現(xiàn)欲運用數(shù)值解方法,將求解混合料在鋪筑過程沿厚度方向的一維非穩(wěn)態(tài)溫度場這一連續(xù)函數(shù)轉(zhuǎn)化為確定沿厚度方向任意位置處每隔一定時間間隔的溫度值,進而可在文獻[6]的邊界條件基礎(chǔ)上,基于MATLAB、C語言等高級數(shù)學工具以實現(xiàn)混合料一維鋪筑溫度場的仿真。

1 基于數(shù)值解法的仿真模型

1.1 有限差分法

1.2 建立控制方程及定解條件

鋪筑溫度場的控制方程為

(1)

式(1)中：a為熱擴散系數(shù),m2/s，a=λ/ρc，其中，ρ、λ、c分別為材料的密度、導(dǎo)熱系數(shù)和比熱容,單位分別為kg/m3、W/(m2·℃)、J/(kg·℃)；?t/?τ為非穩(wěn)態(tài)項；?2t/?x2為擴散項；t為材料的溫度，℃；τ為時間，s；x為路面厚度，m。

定解條件包括初始條件和邊界條件。由于壓路機緊跟攤鋪機,鋪層沿厚度方向的各節(jié)點(除去與下承層的接觸面)的初始條件為攤鋪溫度,接觸面處的初始條件為下承層表面溫度；下承層沿厚度方向的各節(jié)點的初始條件按照式(23)和式(24)計算。邊界條件將在1.4節(jié)中做具體闡述。

圖1 路面結(jié)構(gòu)示意圖Fig.1 Schematic diagram of pavement structure

1.3 建立內(nèi)節(jié)點離散方程

由于隱式差分格式均無條件穩(wěn)定,步長Δx和Δτ的取值不受限制,因此下述離散方程均采用隱式差分格式。

采用泰勒級數(shù)展開法,當非穩(wěn)態(tài)項取向后差分(從p+1時刻的角度觀察),擴散項取p+1時刻的中心差分,即

(2)

(3)

將式(2)和式(3)代入式(1)可得鋪層內(nèi)節(jié)點離散方程:

(4)

下承層內(nèi)節(jié)點離散方程同理:

(5)

1.4 建立邊界節(jié)點離散方程

建立邊界節(jié)點離散方程采用熱平衡法。

由能量守恒定律得鋪層上邊界節(jié)點離散方程:

(6)

鋪層下邊界元體N是由鋪層材料和下承層材料共同組成的寬度為(Δx+ΔX)/2的元體,要將該元體的熱物性參數(shù)重新求解并代入到鋪層下邊界節(jié)點離散方程中:

(7)

下承層上邊界節(jié)點離散方程同式(7)。

由能量守恒定律得下承層下邊界節(jié)點離散方程:

(8)

1.5 建立下承層兩不同種材料接觸節(jié)點離散方程

下承層中兩不同種材料接觸節(jié)點離散方程與內(nèi)節(jié)點離散方程同理,只是將其換成對應(yīng)的熱物性參數(shù)、空間節(jié)點及步長,以第一和第二層下承層材料接觸節(jié)點為例建立離散方程:

(9)

式(9)中：ρ23為第一和第二層下承層材料接觸節(jié)點的密度,單位同前，ρ23=(ρ2+ρ3)/2；C23為第一和第二層下承層材料接觸節(jié)點的比熱容,單位同前，c23=(ρ2c2+ρ3c3)/(ρ2+ρ3)；λ23為第一和第二層下承層材料接觸節(jié)點的導(dǎo)熱系數(shù),單位同前，λ23=2λ2λ3/(λ2+λ3)。

1.6 數(shù)學仿真模型的實現(xiàn)方法

由于采用有限差分法求解一維非穩(wěn)態(tài)問題的代數(shù)方程組計算量較大,考慮到實際操作性,運用MATLAB數(shù)學軟件,建立M函數(shù)文件并合理調(diào)用,將代數(shù)方程組編寫代碼,進而解出各空間節(jié)點在各時間節(jié)點下的溫度值。

宏觀上講,此處涉及的代數(shù)方程組的系數(shù)矩陣M為三對角矩陣,由于M的特殊性,在編寫代碼時采用LU分解法求解。設(shè)未知數(shù)矩陣(溫度矩陣)為T,常數(shù)項矩陣為Q,顯然代數(shù)方程組數(shù)學表達式為MT=Q。令系數(shù)矩陣M由矩陣L和U的乘積表示,即M=LU,可推得LUT=Q,令UT=Y,則求解代數(shù)方程組MT=Q可轉(zhuǎn)化為求解LY=Q,如此,大大節(jié)省了內(nèi)存和計算時間。

2 仿真模型的各參數(shù)取值

2.1 基本假定

瀝青混合料壓實過程中鋪層上邊界暴露在空氣中,主要涉及對流傳熱和輻射傳熱；鋪層下邊界只涉及熱傳導(dǎo)。 現(xiàn)建立如下假設(shè):鋪筑溫度場為一維非穩(wěn)態(tài)溫度場；道路用材料均為均質(zhì)的各向同性材料； 無內(nèi)熱源； 層間無熱阻； 路面兩旁路緣石以及研究深度范圍的最底部絕熱。 取鋪層表面至下承層的42 cm深度處為研究深度范圍。

2.2 熱傳導(dǎo)影響因素的參數(shù)取值

2.2.1 鋪層厚度變化規(guī)律

由文獻[5]的觀測可知,瀝青混合料的鋪層厚度在三遍碾壓后將不再明顯壓縮。層厚隨壓實次數(shù)的變化按文獻[5]的計算方法推算。

Li=(1+tiRL)L

(10)

式(10)中：Li為壓實i次后的層厚,i取0～3的整數(shù),mm；L為設(shè)計厚度,mm；ti為壓實i次后層厚變化系數(shù),t0～t3分別取1、0.35、0.15、0；RL為瀝青層松鋪率,取0.25。

2.2.2 鋪層材料密度變化規(guī)律

由文獻[5]可知,推算模型如式(11)所示：

(11)

式(11)中：ρi為不同壓實次數(shù)后的密度,kg/m3；ρcom為混合料最終壓實密度,ρcom=Kρd。其中，ρd為試驗室成型密度；K為標準壓實度,%。

2.2.3 鋪層材料導(dǎo)熱系數(shù)變化規(guī)律

文獻[2]的熱物性測試結(jié)果表明,在25～75 ℃內(nèi),壓實中的密級配瀝青混合料(dense-graded asphalt mixtures, AC)的導(dǎo)熱系數(shù)介于2.3～2.5 W/(m2·℃),而瀝青瑪蹄脂碎石混合料(stone matrix asphalt, SMA)的導(dǎo)熱系數(shù)介于1.4～1.6 W/(m2·℃)。上述所有試樣對溫度均表現(xiàn)出相似的導(dǎo)熱性能:溫度從25 ℃上升到75 ℃,導(dǎo)熱系數(shù)約減少0.2 W/(m2·℃)。將上述測定值進行回歸,得到細粒和中粒式AC的導(dǎo)熱系數(shù)隨溫度變化的規(guī)律為

λ=-0.004T+2.6

(12)

式(12)中：λ為瀝青混合料導(dǎo)熱系數(shù),W/(m2·℃)；T為瀝青混合料溫度, ℃ 。

SMA的導(dǎo)熱系數(shù)隨溫度變化的規(guī)律為

λ=-0.004T+1.7

(13)

文獻[8]采用熱線法進行了不同溫度下導(dǎo)熱系數(shù)的測定,得到了不同溫度下的混合料導(dǎo)熱系數(shù),將其粗粒式AC的測定值進行回歸,得粗粒式AC的導(dǎo)熱系數(shù)一般性方程為

λ=-0.018T+3.406

(14)

2.2.4 鋪層材料比熱容變化規(guī)律

將文獻[9]混合料的不同集料的比熱容與溫度進行線性回歸,再用文獻[5]的方法建立混合料的不同集料的比熱容與溫度的關(guān)系式:

c=KT+1 000

(15)

式(15)中：K為集料比熱容溫變系數(shù),花崗巖1.581,玄武巖2.262,石灰?guī)r3.006；T為瀝青混合料溫度, ℃；c為瀝青混合料的比熱容,J/(kg·℃)。

2.2.5 下承層材料熱物性參數(shù)取值

在壓實過程中,取下承層材料熱物性參數(shù)為定值,如表1所示。

2.3 對流傳熱影響因素的參數(shù)取值

流體流過時與瀝青鋪層間的熱量傳遞過程,稱為對流傳熱,而影響對流傳熱量大小的參數(shù)是表面?zhèn)鳠嵯禂?shù),現(xiàn)采用Solaimanian等[10]建立的經(jīng)驗公式來計算表面?zhèn)鳠嵯禂?shù)h。

表1 不同下承層材料的熱物性參數(shù)Table 1 Thermophysical parameters of different underlying materials

(16)

式(16)中：v為風速,m/s；Ts為路面表面溫度, ℃；Ta為氣溫, ℃；h為表面?zhèn)鳠嵯禂?shù),W/(m2·℃)。

2.4 輻射傳熱影響因素的參數(shù)取值

2.4.1 太陽輻射條件

采用文獻[11]中的計算式來預(yù)估路面在晴天時的太陽輻射:

(17)

式(17)中：qDN為太陽直接輻射強度,W/m2；qdH為大氣散射和逆輻射強度,W/m2；A為表觀太陽輻射強度,W/m2；B、C分別為大氣的消光系數(shù)和散射因子,無量綱；HA為太陽高度角，按式(18)計算:

sinHA=sinφsinδ+cosφcosδcosω

(18)

式(18)中，φ為緯度, 0～90°；δ為太陽赤緯,rad。

(19)

式(19)中：n為日數(shù),1月1日起算；ω為太陽時角,(°)，計算公式為

ω=15(ST-12)

(20)

式(20)中：ST為真太陽時,以24時計，ST=北京時間+時差，時差=(當?shù)亟?jīng)度-120°)/15°。

用太陽輻射透過率[5]對晴天太陽輻射進行修正,以得到實際天空狀態(tài)下的太陽輻射:

k=100-8.121 9n

(21)

式(21)中：k為太陽輻射透過率,%；n為云量,云量不到天空的0.5/10時總云量為0,占全天1/10時總云量為1,其余依此類推。

取太陽輻射的吸收率α=0.95[5]。

2.4.2 有效輻射條件

晴天有效輻射強度qreff,sunny的計算[5]如式(22)所示:

(22)

式(22)中：ε為路面發(fā)射率；σ為斯忒藩-玻耳茲曼常量,其值為5.67×10-8W/(m2·K4)；tw為路面表面溫度,K；tsky為有效天空溫度,K,計算公式[12]為

(23)

式(23)中：tsurf為地面溫度,取下承層20 cm處的溫度,K,按式(24)推算[13]。

(24)

式(24)中：td為距下承層表面某深度d處的溫度, ℃；tS0為下承層表面溫度, ℃；d為距下承層表面的深度,m；τ為以小數(shù)計的時間參數(shù)。

采用式(24)推算下承層表面至20 cm深度處的溫度場,而20 cm深度以下的溫度場取定值,其值為20 cm 處的溫度[5]。

下承層表面溫度tS0可由實時氣溫Ta計算[13]:

tS0=1.1Ta+1.5+0.17e0.126Ta

(25)

由式(22)～式(25)計算晴天時的有效輻射,再根據(jù)實際的云量和云狀對進行修正[14]:

qreff=(1-CN)qreff,sunny

(26)

式(26)中：qreff為實際有效輻射,W/m2；qreff,sunny為晴天有效輻射,W/m2；N為以小數(shù)計的總云量,0～1.0；C為系數(shù),C=(Chnh+Cmnm+Clnl)/n,其中nh、nm和nl分別為高、中和低云的云量,Ch=0.15～0.2,Cm=0.5～0.6,Cl=0.7～0.8,為不同云狀系數(shù)。

對于有效輻射發(fā)射率,取定值ε=0.93[5]。

2.5 壓實過程

設(shè)共壓實8遍,其中初壓、復(fù)壓和終壓分別為2遍、4遍和2遍。取壓路機3～4.2 km/h碾壓,壓路機緊跟攤鋪機,初壓同一點壓實間隔為150 s；復(fù)壓和終壓同一點壓實間隔為240 s,共計壓實時間為1 740 s。

2.6 表面灑水散熱的參數(shù)取值

取壓路機灑水量200～400 kg/h,碾壓速度約3 km/h,碾輪寬2 m,則單位面積單次碾壓灑水量為33～67 g/m2,變成水蒸汽需熱量82.5～167.5 kJ,設(shè)這些熱量在碾壓后60 s內(nèi)耗散,得灑水散熱熱流密度為1 375～2 792 W/m2,由2.5節(jié),在仿真試驗中施加的路表面灑水散熱熱流密度在1 375～2 792 W/m2范圍內(nèi)取定值，施加灑水散熱熱流密度的時刻為第0、第3、第6、第10、第14、第18、第22和第26時間節(jié)點。

3 仿真模型檢驗

對仿真模型按照不同施工條件下鋪筑過程的溫度實測值進行驗證,施工條件如表2所示,并將驗證結(jié)果繪制于圖2。

表2 用于仿真模型檢驗的不同施工條件Table 2 Simulation model verification different construction conditions

圖2 不同施工條件下的仿真模型檢驗結(jié)果Fig.2 Test results of simulation models underdifferent construction conditions

由圖2(a)和圖2(b)可知,模型的仿真值與實測值在走勢和數(shù)值上基本吻合；層中的仿真值與實測值有很高的一致性；層表的仿真值相對實測值有一定出入,但總體趨勢和最終值有較好的準確性。圖2(c)為包含下承層結(jié)構(gòu)的全厚度仿真模型檢驗,各位置處的仿真值基本接近于實測值。

4 結(jié)論

(1)采用靜態(tài)定值氣溫、風速、太陽輻射、表面灑水量、下承層初始溫度場以及下承層材料熱物性參數(shù),動態(tài)變化值表面?zhèn)鳠嵯禂?shù)、有效輻射、鋪層材料的密度、導(dǎo)熱系數(shù)和比熱容,在MATLAB軟件中建立M函數(shù)文件并合理調(diào)用,得出了在特定條件下一維非穩(wěn)態(tài)鋪筑溫度場的數(shù)值解。

(2)基于傳熱學理論和模型參數(shù),確定了熱態(tài)瀝青混合料鋪筑溫度場的有限差分仿真模型涉及的參數(shù)取值。

(3)仿真模型可根據(jù)施工條件預(yù)估鋪筑溫度場的降溫規(guī)律。不同條件下仿真檢驗證明,該模型具有較好的準確性和可信性,可作為仿真試驗和實際應(yīng)用的基礎(chǔ)。