黃琴龍 吳鳴濤 蔡爵威
(1.同濟(jì)大學(xué)道路與交通工程教育部重點(diǎn)實(shí)驗(yàn)室 上海 201804; 2.上海市建筑科學(xué)研究院有限公司 上海 201108)
道面濕滑嚴(yán)重影響飛機(jī)的制動能力和轉(zhuǎn)彎性能,極易造成飛機(jī)“漂滑”現(xiàn)象[1],對飛機(jī)起降的安全性構(gòu)成極大威脅。由于飛機(jī)滑跑速度高,機(jī)場道面抗滑性能對水膜厚度的變化非常敏感,有效獲取道面水膜厚度的動態(tài)演變信息,有助于識別與預(yù)警濕滑狀態(tài)下飛機(jī)滑跑安全風(fēng)險(xiǎn),是機(jī)場安全運(yùn)營工作和“四型機(jī)場”建設(shè)的重要內(nèi)容。
近年來,國內(nèi)外許多學(xué)者通過理論推導(dǎo)和試驗(yàn)數(shù)據(jù)擬合等手段,對道面水膜厚度分布開展研究。張理等[2]采用謝才公式和曼寧公式進(jìn)行理論推導(dǎo),得到了不同降雨強(qiáng)度和道路橫縱坡條件下的水膜厚度公式。李光元等[3]采用極差分析法得到水膜厚度影響因素的主次順序,回歸分析了人工降雨試驗(yàn)數(shù)據(jù),得到了水泥混凝土道面的水膜厚度預(yù)測模型。尚彥宇等[4]通過零縱坡軸推導(dǎo)了路面徑流特征,建立了超高漸變段路面模型??傮w而言,上述模型適用于降雨強(qiáng)度穩(wěn)定不變的情況,但無法實(shí)時(shí)、動態(tài)求解道面水膜厚度。
本文基于運(yùn)動波方程和水力學(xué)基礎(chǔ)理論,擬結(jié)合道面紋理,橫、縱坡排水,以及實(shí)時(shí)降雨強(qiáng)度,提出道面水膜厚度求解的基本控制方程;采用Preissmann法和MacCormack法對模型進(jìn)行動態(tài)求解,并在恒定降雨強(qiáng)度和變化降雨強(qiáng)度2種條件下進(jìn)行模型驗(yàn)證。
道面水膜厚度的動態(tài)演變,可采用一維模型和二維模型進(jìn)行計(jì)算。二維模型采用有限元方法,精細(xì)化模擬道面的真實(shí)紋理狀態(tài)和水流的運(yùn)動,從而精確求解道面水膜厚度,但計(jì)算耗費(fèi)時(shí)間較長、概化困難,且需要高分辨率高程數(shù)據(jù)進(jìn)行支撐,二維模型過大的計(jì)算量會極大增加系統(tǒng)負(fù)擔(dān)。而一維模型根據(jù)道面高程數(shù)據(jù)和計(jì)算網(wǎng)格密度,提取道面內(nèi)水流流線,并計(jì)算各流線上水膜厚度;通過各流線的匯流關(guān)系,計(jì)算道面各網(wǎng)格的水膜厚度值,最終獲取道面水膜厚度的動態(tài)演變,具有計(jì)算效率高、計(jì)算時(shí)間短、所需邊界信息及高程數(shù)據(jù)較少等優(yōu)點(diǎn)[5]。因此本研究采用一維模型對道面水膜厚度的動態(tài)演變進(jìn)行求解。
為了獲取道面水膜厚度動態(tài)演變,首先需要確定道面徑流的流線。在網(wǎng)格單元劃分的基礎(chǔ)上,假定網(wǎng)格單元內(nèi)最大坡降方向與坡度均為定值,區(qū)域內(nèi)最大坡降方向即為流線方向,由此可計(jì)算出網(wǎng)格邊界任意點(diǎn)起始的流線軌跡。對于不同網(wǎng)格間流線軌跡的銜接,上游網(wǎng)格單元內(nèi)的流線終點(diǎn)即為下游網(wǎng)格單元內(nèi)的流線起點(diǎn),由此可將流線從跑道中線逐步推算到道肩。
通過對網(wǎng)格進(jìn)行編號(Xj,Yk)確定單元及其對應(yīng)位置,提取道面的流線,網(wǎng)格間流線組合示意圖見圖1。
圖1 道面網(wǎng)格劃分與網(wǎng)格間組合流線示意圖
流線與坡降方向平行,對于流線經(jīng)過的任意網(wǎng)格,起點(diǎn)不同的流線不會相交;對于任意起點(diǎn)在單元內(nèi)的流線終點(diǎn),均可以通過其在單元內(nèi)的位置和坡降方向信息計(jì)算得到,并代入下一單元內(nèi)繼續(xù)計(jì)算。
在任意網(wǎng)格單元(Xj,Yk)內(nèi),流線起點(diǎn)位置為(x0,y0),流線終點(diǎn)位置為(x1,y1)。當(dāng)流線起點(diǎn)在x方向時(shí)(即x0≥0,y0=0),對于不同起點(diǎn)和坡度的網(wǎng)格單元,流線終點(diǎn)可能在坡降方向所指向的2條邊界(即y方向邊界、x方向邊界)及其交點(diǎn)(x方向邊界和y方向邊界的交點(diǎn))處,其所對應(yīng)下一個(gè)單元的網(wǎng)格編號分別為(Xj+1,Yk)、(Xj,Yk+1)、(Xj+1,Yk+1)。圖1右側(cè)圖為流線起點(diǎn)在x方向時(shí)的示意圖,計(jì)算方法如式(1)所示。
(1)
式中:LX為網(wǎng)格單元在x方向的長度,m;LY為網(wǎng)格單元在y方向的長度,m;i為網(wǎng)格單元內(nèi)的綜合坡度,%;ix為網(wǎng)格單元內(nèi)x方向的坡度,%;iy為網(wǎng)格單元內(nèi)y方向的坡度,%。
當(dāng)流線起點(diǎn)在y方向時(shí)(即x0=0,y0>0),只需將式(1)中的x與y、X與Y互換即可。
水泥混凝土道面徑流可視為坡面流的匯流過程,可采用一維圣維南方程進(jìn)行表征,其基本方程包括連續(xù)方程和動力方程,分別如式(2)、式(3)所示。
(2)
(3)
式中:t為時(shí)間;q為降雨強(qiáng)度;x為距道面徑流某固定斷面沿流程的距離;h、u、Q分別為x處過水?dāng)嗝娴乃?即水膜厚度)、斷面平均流速、斷面平均流量;g為重力加速度;S0、Sf分別為坡底源項(xiàng)和摩擦源項(xiàng)。
式(2)是連續(xù)方程,第一項(xiàng)為沿程流量的變化率,第二項(xiàng)為蓄量的變化率,第三項(xiàng)為降雨強(qiáng)度,反映了流線中的水量平衡;式(3)是動力方程,第一項(xiàng)為局部加速度項(xiàng),第二項(xiàng)為對流加速度項(xiàng),第三項(xiàng)為壓力梯度項(xiàng),第四項(xiàng)為水流內(nèi)部及邊界的摩阻損失,反映了水流克服壓力和摩阻引起的能量損失而獲得加速度。
道面徑流屬于地表漫流,可將其流動視為運(yùn)動波,因此可以忽略動力方程中的局部加速度項(xiàng)、對流加速度項(xiàng)、壓力梯度項(xiàng),將式(3)調(diào)整為式(4)。
g(Sf-S0)=0
(4)
曼寧公式作為一種經(jīng)驗(yàn)公式,常用于明渠流水力計(jì)算。道面上薄層水流狀態(tài)是紊流,具有自由表面,且重力是其流動的主導(dǎo)因素,因此將曼寧公式帶入謝才公式,可獲得流速u與糙率n、水力坡度J,以及水力半徑R的關(guān)系,如式(5)所示。
(5)
式中:水力半徑R為過水?dāng)嗝婷娣eA與濕周C的比值,其值約等于水膜厚度h,即R=A/C≈h;水力坡度J為比降,即道面橫坡i,即J=i。代入式(5)可得式(6)。
(6)
式(6)兩邊同乘h,得到流量Q與水膜厚度h、糙率n及道面橫坡i的關(guān)系,如式(7)所示。
(7)
由于道面橫坡i足夠小,tani近似與i相等。因此將式(7)調(diào)整為式(8)。
(8)
式中:m= 5/3;α= (tani)1/2/n;MH/T 5036-2017 《民用機(jī)場排水設(shè)計(jì)規(guī)范》中給出的瀝青混凝土或水泥混凝土道面的糙率為0.013,取n=0.013。
聯(lián)立式(8)與式(2),得到道面水膜厚度求解的基本控制方程,如式(9)所示。
(9)
采用Preissmann四點(diǎn)隱式差分法和MacCormack差分法,將式(9)在時(shí)間和空間2個(gè)維度上,進(jìn)行離散化并求解。
對于函數(shù)f(x,t),按式(10)計(jì)算各因變量對其導(dǎo)數(shù)的離散形式。
(10)
式中:θ為加權(quán)系數(shù),0 <θ≤1。
將式(10)帶入式(9)可得
(11)
MacCormack差分法按照預(yù)測-修正的方式,對式(9)進(jìn)行求解。預(yù)測步、修正步分別如式(12)、式(13)所示。
(12)
(13)
在求解過程中,取未降雨(或剛開始降雨)的時(shí)刻為 0 時(shí)刻,且流線起點(diǎn)為道面高程的極大值點(diǎn),整個(gè)模型的邊界條件滿足式(14)。
(14)
季天劍等[6]選取了AC-25、AC-20和SMA-13 3種典型道面結(jié)構(gòu),開展了道面水膜厚度試驗(yàn),建立了計(jì)算水膜厚度的回歸方程。本文參考其試驗(yàn)條件,將其試驗(yàn)數(shù)據(jù)作為實(shí)測值,對比分析本文模型的水膜厚度計(jì)算值,結(jié)果見圖2。
圖2 SMA-13改進(jìn)型路面結(jié)構(gòu)實(shí)測值與計(jì)算值對比
圖2a)選取恒定降雨強(qiáng)度5.4 mm/min,排水長度0~8 m共9種工況,對水膜厚度實(shí)測值與計(jì)算值進(jìn)行比較,可以看出,在不同排水長度下,水膜厚度的計(jì)算值與實(shí)測值吻合度較高。圖2b)選取恒定降雨強(qiáng)度3.6,5.4,6.3 mm/min,排水長度0~8 m共27種工況,對水膜厚度實(shí)測值與計(jì)算值進(jìn)行比較,可以看出,在不同排水長度、不同降雨強(qiáng)度下,水膜厚度的計(jì)算值與實(shí)測值仍保持較高吻合度,兩者的平均絕對誤差為0.54 mm,平均精度為85.1%??紤]到文獻(xiàn)[6]中水膜厚度測量分辨率為0.5 mm,因此可以認(rèn)為,在恒定降雨強(qiáng)度下,本文提出的計(jì)算模型具有較高精度。
采用本文提出的模型,計(jì)算變化降雨強(qiáng)度及坡度情況下的水膜厚度演變情況。隨機(jī)生成30 min的降雨強(qiáng)度序列,計(jì)算結(jié)果見圖3。由圖3可見,水膜厚度演變情況與降雨強(qiáng)度變化呈現(xiàn)出較好的一致性,模型很好地反映了在降雨強(qiáng)度改變時(shí)水膜厚度的演變情況,以及降雨停止后積水的消散,避免了回歸模型中的結(jié)果突變現(xiàn)象。
圖3 變化降雨強(qiáng)度下水膜厚度演變情況
獲取流線上每一點(diǎn)處的水膜厚度時(shí)空演變情況見圖4。
圖4 水膜厚度的時(shí)空演變特征
由圖4可見,流線上每一點(diǎn)處的水膜厚度隨時(shí)間的變化情況都被非常直觀地呈現(xiàn)了出來,相當(dāng)于得到了該面域的水膜厚度時(shí)空演變情況。若將網(wǎng)格劃分地更細(xì),得到的結(jié)果精度會更高,但與此同時(shí)消耗的計(jì)算力也越大,考慮到模型預(yù)測的實(shí)時(shí)性,需要根據(jù)實(shí)際情況在計(jì)算速度和精度之間取一定的平衡。
本文采用Pressimann四點(diǎn)隱式差分法和MacCormack差分法對控制方程進(jìn)行求解。隨機(jī)生成30 min的降雨強(qiáng)度系列計(jì)算結(jié)果分別見圖5。
圖5 Preissmann法和MacCormack法的計(jì)算結(jié)果
由圖5可見,2種算法的計(jì)算結(jié)果基本一致,但在求解過程中,Preissmann四點(diǎn)隱式差分法的計(jì)算時(shí)間為64 s,MacCormack差分法的計(jì)算時(shí)間為36 s,后者耗時(shí)僅為前者56.25%??紤]到機(jī)場系統(tǒng)的數(shù)據(jù)體量巨大,MacCormack差分法更有助于快速處理并獲取道面水膜厚度的動態(tài)演變信息,有效提升風(fēng)險(xiǎn)預(yù)警的安全裕度。因此在后續(xù)研究中,可采用MacCormack差分法對模型進(jìn)行求解。
北京首都國際機(jī)場安全預(yù)警平臺布設(shè)了水膜厚度感知系統(tǒng),安裝了Vaisala MD30遙感式水膜厚度傳感器和脈沖式雨量計(jì),實(shí)時(shí)監(jiān)測跑道區(qū)域的水膜厚度和降雨強(qiáng)度。在安裝區(qū)域內(nèi),跑道橫坡和縱坡分別為1.2%和0.15%。
采用該系統(tǒng)2020年11月18日的降雨強(qiáng)度和水膜厚度監(jiān)測數(shù)據(jù),對本文模型進(jìn)行實(shí)測驗(yàn)證,計(jì)算結(jié)果見圖6。
圖6 模型計(jì)算值與實(shí)測值的結(jié)果
由圖6可見,模型的計(jì)算值與實(shí)測值具有很好的一致性,均方根誤差為0.17 mm。因此,在變化降雨強(qiáng)度下,本文模型具有較高精度,能夠精確、實(shí)時(shí)地描述道面水膜厚度的動態(tài)演變。
1) 本文基于運(yùn)動波方程,結(jié)合道面動態(tài)時(shí)空信息,建立了道面水膜厚度動態(tài)演變的求解模型,并通過實(shí)測數(shù)據(jù)驗(yàn)證了模型的可靠性。
2) 在恒定降雨強(qiáng)度下,本文模型能夠精確計(jì)算出道面水膜厚度,平均絕對誤差為0.54 mm,平均相對誤差為14.9%。
3) 在變化降雨強(qiáng)度下,本文模型能夠很好地反映道面水膜厚度的演變特征,并且模型結(jié)果與實(shí)測數(shù)據(jù)一致性較好,均方根誤差為0.17 mm。
4) 采用了Preissmann四點(diǎn)隱式差分法和MacCormack差分法進(jìn)行模型求解,后者在求解耗時(shí)僅為前者的56.25%,可在后續(xù)研究中優(yōu)先考慮。