張 影,楊曉東,賈子璇,章 胤,王 晶
(1.燕山大學(xué)里仁學(xué)院,河北 秦皇島 066004; 2.燕山大學(xué)理學(xué)院,河北 秦皇島 066004;3.燕山大學(xué)經(jīng)濟(jì)管理學(xué)院,河北 秦皇島 066004; 4.燕山大學(xué)區(qū)域經(jīng)濟(jì)發(fā)展研究中心,河北 秦皇島 066004)
水污染評(píng)估是研究水污染問(wèn)題的重要途徑,利用水動(dòng)力水質(zhì)數(shù)值模型模擬計(jì)算水污染的輸移擴(kuò)散過(guò)程,可為水污染風(fēng)險(xiǎn)評(píng)估和相關(guān)部門(mén)采取應(yīng)對(duì)措施提供重要參考信息。水庫(kù)作為水資源重要存儲(chǔ)地,庫(kù)區(qū)上游突發(fā)性水污染的研究一直備受關(guān)注。針對(duì)庫(kù)區(qū)上游河床特點(diǎn),往往建立一維、二維耦合水動(dòng)力水質(zhì)模型模擬污染物的擴(kuò)散,于庫(kù)區(qū)中上游河道窄長(zhǎng)區(qū)域構(gòu)建一維河網(wǎng)模型[1],下游河道開(kāi)闊水域構(gòu)建二維模型,并在交接面處進(jìn)行合理的耦合,既可較準(zhǔn)確地反映污染物濃度的時(shí)空變化規(guī)律,同時(shí)又可控制模型的復(fù)雜度以滿足實(shí)際需求[2]。
近年來(lái),眾多學(xué)者針對(duì)一維、二維耦合模型開(kāi)展了大量研究工作。諸裕良等[3]在河口一維、二維連接處通過(guò)接口斷面法傳遞水力因子,建立了一種一維、二維全隱河網(wǎng)海灣水動(dòng)力聯(lián)網(wǎng)數(shù)學(xué)模型;Lai等[4]提出了適用于大尺度水動(dòng)力模擬的一維、二維耦合方法,并對(duì)長(zhǎng)江中游流域的河流-湖泊洪水演進(jìn)進(jìn)行了模擬;陳文龍等[5]構(gòu)造并求解Riemann問(wèn)題實(shí)現(xiàn)了一維、二維模型耦合,有效克服了傳統(tǒng)堰流公式缺點(diǎn),并提出時(shí)間步長(zhǎng)自適應(yīng)匹配方法解決了一維、二維模型時(shí)間步長(zhǎng)不一致的問(wèn)題;王秀杰等[6]建立了復(fù)雜條件下天然河道漫潰堤洪水在防洪保護(hù)區(qū)的一維、二維水動(dòng)力模型;顧杰等[2]基于MIKE FLOOD建立了入海河流及近岸海域一維、二維耦合河流-海岸水動(dòng)力和水質(zhì)模型;田福昌等[7]建立了山洪溝道潰堤洪水演進(jìn)一維、二維水動(dòng)力耦合數(shù)值模擬模型以分析評(píng)估潰堤山洪淹沒(méi)風(fēng)險(xiǎn)。目前,一維、二維耦合模型研究多為基于堰流公式和數(shù)值通量的側(cè)向型聯(lián)解耦合,構(gòu)建水動(dòng)力的一維、二維耦合模型以研究漫潰堤洪水問(wèn)題,而針對(duì)一維、二維水動(dòng)力水質(zhì)耦合模型的研究較少,多為MIKE軟件的實(shí)際應(yīng)用。此外,針對(duì)耦合交接面的連接問(wèn)題,對(duì)參數(shù)及連接條件的過(guò)度概化,往往造成模擬結(jié)果的不準(zhǔn)確。
針對(duì)上述問(wèn)題,本文在傳統(tǒng)重疊投影法[8]的基礎(chǔ)上,于重疊區(qū)域的上邊界引入流速等物理量的二次橫向分布函數(shù),基于神經(jīng)網(wǎng)絡(luò)訓(xùn)練思想,采用最速下降法構(gòu)建訓(xùn)練算法,對(duì)分布函數(shù)進(jìn)行優(yōu)化訓(xùn)練,將收斂后的分布函數(shù)作為模型的連接條件,進(jìn)而建立了基于改進(jìn)重疊投影法的空間耦合模型,有效地降低了耦合處物理量的突變,實(shí)現(xiàn)了一維、二維水動(dòng)力和水質(zhì)模型精準(zhǔn)耦合。
采用描述明渠非恒定流的水動(dòng)力運(yùn)動(dòng)方程:
(1)
式中:x為沿程距離;t為時(shí)間;z為水位;Q為流量;qt為旁側(cè)入流流量;B為寬度;A為斷面面積;K為流量模數(shù);g為重力加速度。
結(jié)合河網(wǎng)汊點(diǎn)連接條件,利用Preissmann四點(diǎn)加權(quán)隱式格式[9]離散控制方程,然后采用河網(wǎng)三級(jí)聯(lián)解法求解(詳見(jiàn)文獻(xiàn)[10])。
河道污染物的一維對(duì)流擴(kuò)散方程為
(2)
式中:c為水流輸送的污染物濃度;Ex為污染物縱向擴(kuò)散系數(shù);K1為污染物綜合降解系數(shù)。
對(duì)方程(2)采用前差分離散時(shí)間項(xiàng)、隱式迎風(fēng)格式離散對(duì)流項(xiàng)、中心差分離散擴(kuò)散項(xiàng),整理得各斷面濃度遞推關(guān)系式,并結(jié)合汊點(diǎn)質(zhì)量平衡方程建立方程組求解,回代至各河段得各斷面污染物濃度(詳見(jiàn)文獻(xiàn)[11])。
二維控制方程有水流連續(xù)性方程和水流運(yùn)動(dòng)方程組:
(3)
(4)
式中:H為水深;u、v為x、y方向的流速;f為阻力系數(shù);Ω為科氏力系數(shù);ν為紊流渦黏性系數(shù);λ為風(fēng)應(yīng)力系數(shù);wx、wy為風(fēng)速在x、y方向上的分量。
采用貼體坐標(biāo)法[12-13],通過(guò)Poison方程對(duì)控制方程進(jìn)行坐標(biāo)變換,并采用交替方向隱格式法(ADI)求解(詳見(jiàn)文獻(xiàn)[14])。
二維對(duì)流擴(kuò)散方程用以描述污染物在水體中的輸移擴(kuò)散規(guī)律:
(5)
式中Ex、Ey為x、y方向的擴(kuò)散系數(shù)。
對(duì)流擴(kuò)散方程經(jīng)坐標(biāo)變換后,基于水動(dòng)力模型計(jì)算的流速、水位等結(jié)果,采用交替方向隱格式法求解(詳見(jiàn)文獻(xiàn)[14])。
采用一維、二維耦合的水動(dòng)力和水質(zhì)模型模擬污染物的擴(kuò)散時(shí),在耦合交接面的連接區(qū)域上對(duì)參數(shù)及連接條件的過(guò)度概化,往往造成模擬結(jié)果的不準(zhǔn)確。一般情況下,沿河寬方向各物理量——如水位、流速及污染物濃度等——不是均勻分布的,為體現(xiàn)連接后二維區(qū)域各物理量的真實(shí)分布情況,在傳統(tǒng)重疊投影法的基礎(chǔ)上,于一維、二維區(qū)域引入水位、流速、污染物濃度等物理量的二次橫向分布函數(shù),進(jìn)而構(gòu)建訓(xùn)練算法,對(duì)各分布函數(shù)中的未知參數(shù)進(jìn)行訓(xùn)練,將收斂后的分布函數(shù)作為模型的連接條件,建立改進(jìn)的一維、二維水動(dòng)力和水質(zhì)空間耦合模型。
考慮在一般情況下沿河寬方向各物理量不是均勻分布,在一維、二維交接面設(shè)置物理量沿河寬方向的分布函數(shù)Ψ(y),并采用二次函數(shù)形式來(lái)擬合其在邊界的分布:
Ψ(y)=az1y2+az2y+az3
(6)
式中:y為河寬方向坐標(biāo);az1、az2、az3為待定系數(shù)。將分布函數(shù)Ψ(y)于二維邊界離散化,設(shè)其在控制體i(i=1,2,…,N,N為二維邊界處控制體個(gè)數(shù))的取值為Ψi,有:
(7)
由于二維邊界物理量隨時(shí)間變化的模擬計(jì)算過(guò)于復(fù)雜,因而本文代之以適用于一定時(shí)段的二維邊界各控制體的波動(dòng)比值ωi,有:
(8)
此時(shí),通過(guò)訓(xùn)練得到波動(dòng)比值向量ω,在已知交接面處的平均值后就可得該物理量的分布情況,能準(zhǔn)確體現(xiàn)其在交接面的連接關(guān)系,從而提高耦合準(zhǔn)確度。
據(jù)此,在二維邊界控制體上設(shè)置水位、流速、污染物濃度分布函數(shù),需要注意的是重疊區(qū)域二維邊界的流速需要縱向流速分布函數(shù)u(y)以及流速矢量方向與縱向的夾角分布函數(shù)θ(y)兩個(gè)分布函數(shù)以表示其分布情況,分布函數(shù)皆采用二次函數(shù)形式擬合,通過(guò)訓(xùn)練可得收斂的水位波動(dòng)比值向量ωz、流速波動(dòng)比值向量ωu、夾角分布向量θ,進(jìn)而得到污染物濃度波動(dòng)比值向量ωc,作為污染物濃度在交接面上的連接條件。
考慮到二維區(qū)域數(shù)值波的傳播受柯朗-弗里德里西斯-列維(Courant-Fredrich-Lewy, CFL)條件限制,進(jìn)行水域延伸構(gòu)造虛擬重疊區(qū)域[8],并考慮在一維準(zhǔn)確計(jì)算水域中進(jìn)行二維區(qū)域物理量的細(xì)化分布函數(shù)的訓(xùn)練,將二維計(jì)算水域的邊界向一維計(jì)算水域延伸Courant個(gè)網(wǎng)格點(diǎn),見(jiàn)圖1。此時(shí),定義左側(cè)虛擬邊界為Ω1,右側(cè)一維、二維交接面為Ω2,其中,Ω1為一維、二維轉(zhuǎn)化平面,Ω2為比較訓(xùn)練平面;圖1中Φ為一維區(qū)域內(nèi)的物理量,Ψ為二維區(qū)域的物理量,L表示重疊區(qū)域中虛擬的一維斷面數(shù)。
圖1 重疊-投影示意圖
訓(xùn)練時(shí),面Ω2采用滯后條件,即下一時(shí)步虛擬重疊區(qū)域的物理量通過(guò)該時(shí)間步長(zhǎng)內(nèi)的迭代計(jì)算得到含參精確值。在一維準(zhǔn)確計(jì)算水域中進(jìn)行二維區(qū)域物理量的細(xì)化分布函數(shù)的訓(xùn)練,利用該精確值與原滯后值之差構(gòu)建目標(biāo)函數(shù)以訓(xùn)練分布函數(shù)相關(guān)系數(shù)。
以訓(xùn)練穩(wěn)定后的分布函數(shù)作為模型一維、二維的連接條件進(jìn)行耦合計(jì)算,進(jìn)而構(gòu)建改進(jìn)的重疊投影區(qū)域。改進(jìn)的重疊投影關(guān)系見(jiàn)圖2。
圖2 重疊區(qū)域一、二維連接關(guān)系示意圖
2.2.1參數(shù)訓(xùn)練步驟
a. 先計(jì)算水動(dòng)力模型,面Ω2的時(shí)間滯后條件取為ΔΦn+1=0(n為模型已完成迭代的次數(shù)),記此時(shí)的邊界值為ΦΩ2。
b. 求解一維隱格式,得面Ω1的物理量值ΦΩ1,并乘以波動(dòng)比值向量ωz、ωu將面Ω1上一維物理量值轉(zhuǎn)化為二維邊界條件ΨΩ1(P1)。
c. 計(jì)算二維區(qū)域。取面Ω2上流速等物理量投影得到一維精確邊界條件Φ′Ω2。
d. 確定目標(biāo)函數(shù)p(P1):
p(P1)=(Φ′Ω2-ΦΩ2)2
(9)
目標(biāo)函數(shù)越小,一維、二維連接效果越好。
2.2.2模型耦合步驟
a. 面Ω2的時(shí)間滯后條件取為ΔΦn+1=0,記此時(shí)的邊界值為ΦΩ2。
b. 將滯后邊界值ΦΩ2作為一維水動(dòng)力模型的下邊界,采用三級(jí)聯(lián)解法得到各水力要素在一維計(jì)算區(qū)域各斷面的取值。
d. 求解二維計(jì)算區(qū)域,代入二維邊界條件ΨΩ1(P1),采用交替方向隱格式法結(jié)合追趕法得水力要素于二維計(jì)算區(qū)域的分布。
e. 取面Ω2的水力要素分布,投影得到下一時(shí)刻一維邊界條件Φ′Ω2,即ΦΩ2,n+2。
f. 將水動(dòng)力模型相關(guān)水力參數(shù)輸入水質(zhì)模型。同水動(dòng)力模型訓(xùn)練耦合步驟,并由面Ω1的污染物濃度波動(dòng)比值向量ωc實(shí)現(xiàn)污染物濃度一維、二維區(qū)域的連接,最終得到整個(gè)計(jì)算區(qū)域內(nèi)的污染物濃度分布。
g. 重復(fù)上述步驟,根據(jù)實(shí)時(shí)訓(xùn)練出的面Ω1的分布函數(shù),實(shí)現(xiàn)一維、二維交接面物理量的轉(zhuǎn)化,直至模擬結(jié)束。
以秦皇島桃林口庫(kù)區(qū)上游青龍河流域?yàn)檠芯繀^(qū)域,模擬青龍縣某污水泄漏事件,對(duì)上述改進(jìn)的一維、二維耦合河網(wǎng)水動(dòng)力水質(zhì)模型進(jìn)行檢驗(yàn)和分析。
桃林口水庫(kù)位于河北省東北部的灤河主要支流青龍河上,總庫(kù)容8.59億m3,區(qū)域內(nèi)地勢(shì)北高南低,降雨主要集中于每年的7、8月。庫(kù)區(qū)上游河流以青龍河水系為主,河長(zhǎng)265 km,控制面積3 431.5 km2,河床主要為砂卵石,中游河床寬50~70 m,下游段河床較寬,為400~1 000 m,平均縱坡0.155 6%。
通過(guò)分析區(qū)域的河段組成,確定一維河網(wǎng)研究范圍為:上邊界取自小柳條溝北莊S358縣道,下邊界取至小菜峪。徑流量較大支流包括沙河、起河、都源河以及星干河。一維模型共布設(shè)了505個(gè)斷面,典型斷面如圖3所示,斷面平均距離Δx=500 m。為滿足CFL條件和控制方程的微分性質(zhì),設(shè)定一維模型的時(shí)間步長(zhǎng)Δt=120 s。
圖3 河網(wǎng)各河段典型斷面位置概化示意圖
二維模型研究范圍包括小菜峪斷面至入庫(kù)口。為適應(yīng)復(fù)雜河道邊界,本文建立貼體坐標(biāo)系將物理計(jì)算區(qū)域進(jìn)行轉(zhuǎn)化[15],最終生成了3 451個(gè)計(jì)算網(wǎng)格?;贑FL條件選擇時(shí)間步長(zhǎng)為4 s。一維、二維計(jì)算結(jié)果的輸出間隔均為120 s。
3.2.1初始條件
以青龍河水系監(jiān)測(cè)站最近(2019年7月15日)的一次檢測(cè)數(shù)據(jù)為基礎(chǔ),通過(guò)數(shù)據(jù)處理,得到各典型斷面的流量、水位、自然狀態(tài)下污染物的濃度等數(shù)據(jù)。采用一維、二維插值方法,分別得到一維、二維區(qū)域各斷面、控制體的流量、水位等數(shù)據(jù)作為水動(dòng)力模型的初始條件;各斷面、控制體的自然狀態(tài)下污染物的濃度作為水質(zhì)模型的初始條件。
3.2.2邊界條件
一維水動(dòng)力模型的上邊界條件是各干支流河段上游實(shí)際的流量、水位變化,下邊界條件為一維、二維耦合界面處二維模型計(jì)算得到的平均流量和平均水位。
二維水動(dòng)力模型的耦合界面為流量-水位邊界,其余邊界為固壁邊界,本文根據(jù)流域特點(diǎn)和模擬時(shí)長(zhǎng),固壁邊界采用閉邊界條件,沿邊界法線(nΓ0)方向流速為零。
河網(wǎng)水質(zhì)模型的邊界條件是單源點(diǎn)污染源,位于青龍縣河南村斷面的滿源污水處理廠的生活污水均勻地排入所在的河流系統(tǒng)中。
以秦皇島市水文局等部門(mén)提供的部分?jǐn)?shù)據(jù)為基礎(chǔ),通過(guò)Google Earth遙感與圖像處理技術(shù)[16-17]設(shè)計(jì)提取算法,提取河網(wǎng)各斷面和控制體的基礎(chǔ)數(shù)據(jù),并結(jié)合河流縱向擴(kuò)散系數(shù)[18]等水力學(xué)經(jīng)驗(yàn)計(jì)算公式,計(jì)算得到模型必要的基本參數(shù)。對(duì)于某些需要率定的參數(shù),按斷面采用試錯(cuò)法調(diào)試有較大任意性,河網(wǎng)規(guī)模較大時(shí)調(diào)試工作量十分巨大,本文采用分級(jí)法,即同一等級(jí)的河道按一定的參數(shù)取值,保證模型參數(shù)的相對(duì)準(zhǔn)確性及可行性,并通過(guò)實(shí)地勘察和咨詢(xún)水力學(xué)專(zhuān)業(yè)人員等途徑,綜合考量,對(duì)一維、二維每個(gè)經(jīng)典斷面和經(jīng)典控制體的水力參數(shù)進(jìn)行打分,采用插值方法計(jì)算所有斷面和控制體的水力參數(shù)。
假定滿源污水處理廠突發(fā)污水泄漏事故,有 9 400 m3的生活污水均勻地排入研究水系中,本文主要研究濃度為1.606×10-4mol/L的總磷(TP)在未來(lái)26.7 h內(nèi)的變化情況[19],并從秦皇島市水文局監(jiān)測(cè)站監(jiān)測(cè)時(shí)間(2019-07-15)開(kāi)始計(jì)時(shí)。
圖4和圖5為河網(wǎng)水動(dòng)力模型在一維區(qū)域的計(jì)算結(jié)果。在整個(gè)模擬期間,最大流量是876.23 m3/s,在一個(gè)時(shí)間步長(zhǎng)(120 s)內(nèi)變化最大為14.17 m3/s,最高水位是216.34 m,一個(gè)時(shí)間步長(zhǎng)內(nèi)變化最大為0.23 m。
圖4 斷面流量的時(shí)空變化
圖5 斷面水位的時(shí)空變化
圖6為河網(wǎng)水質(zhì)模型在一維區(qū)域的計(jì)算結(jié)果,可見(jiàn)污染事件發(fā)生后561個(gè)時(shí)間步長(zhǎng)(即18.7 h)污染物到達(dá)了二維區(qū)域(小菜峪斷面)。值得注意的是,污染物匯入汊點(diǎn)時(shí)未污染水流匯入稀釋?zhuān)瑢?dǎo)致輸出斷面的污染物濃度出現(xiàn)突變。
圖6 污染物TP濃度的時(shí)空變化
圖7為河網(wǎng)模型在二維區(qū)域的計(jì)算結(jié)果,即在污染事件發(fā)生后581、601個(gè)時(shí)間步長(zhǎng)(即19.33 h、20.23 h)時(shí),二維區(qū)域中各控制體上污染物濃度的分布,可以看出,當(dāng)以二次函數(shù)作為交接面污染物濃度分布的擬合函數(shù)時(shí),后續(xù)的二維區(qū)域同一橫截面的控制體上的TP濃度,仍然呈中間高、兩端低近似二次函數(shù)趨勢(shì)的分布,從而驗(yàn)證了用二次函數(shù)作為擬合函數(shù)的合理性。
(a) t=19.33 h
由以上計(jì)算結(jié)果可知,通過(guò)最速下降法訓(xùn)練改進(jìn)后的重疊投影法實(shí)現(xiàn)一維、二維模型的耦合,在河段推演、汊點(diǎn)連接、一維和二維耦合等過(guò)程均未出現(xiàn)解的“奇點(diǎn)”,方法的改進(jìn)并未引起模型解的突變,表明改進(jìn)后的模型穩(wěn)定性較高,達(dá)到了預(yù)期目標(biāo)。
首次訓(xùn)練穩(wěn)定后的一維和二維交接面處縱向流速、流速矢量方向與縱向的夾角、水位和污染物濃度的分布函數(shù)中的參數(shù)如表3所示,可以看到,4個(gè)分布函數(shù)的二次項(xiàng)系數(shù)均較小,表明河流同一截面上的該4個(gè)變量數(shù)值相差不會(huì)太大;夾角θ的二次系數(shù)為負(fù),河流中央夾角小,兩端夾角大,其他3個(gè)變量正好相反;各分布函數(shù)的一次項(xiàng)系數(shù)很小,表明河流中該4個(gè)水力要素近似呈對(duì)稱(chēng)分布。以上結(jié)果近似符合青龍河在二維區(qū)域的河流特征,進(jìn)而驗(yàn)證了參數(shù)訓(xùn)練算法的正確性。
表3 首次訓(xùn)練穩(wěn)定后各物理量分布函數(shù)的系數(shù)
圖8為通過(guò)傳統(tǒng)的重疊投影法、改進(jìn)的重疊投影法計(jì)算的一維、二維交接面與其下游相鄰20 m斷面處的平均縱向流速差值和TP平均濃度差值的變化,傳統(tǒng)重疊投影法計(jì)算的平均縱向流速差值的平均值為0.091 m/s,改進(jìn)的重疊投影法計(jì)算的平均值為0.040 m/s,實(shí)測(cè)的平均值為0.046 m/s,可見(jiàn)改進(jìn)的重疊投影法使一維向二維過(guò)渡時(shí)流速等水動(dòng)力參量的突變性大大減小,保證了水動(dòng)力耦合模型的精確性。由于流速和TP濃度的非均勻分布(河中央流速大處TP濃度高),改進(jìn)后的水質(zhì)耦合模型污染物TP擴(kuò)散速度稍微加快,更加符合實(shí)際情況。
(a) 平均縱向流速差值
由于傳統(tǒng)的重疊投影法在處理一維、二維數(shù)值交換時(shí),是將交接面上各物理量的一維數(shù)值賦予該截面所有二維控制體,即交接面上所有控制體上各物理量的數(shù)值均相等,而河流的實(shí)際情況是各物理量的數(shù)值在橫截面上呈某種分布,因而導(dǎo)致了河網(wǎng)模型的解在一維、二維耦合處出現(xiàn)較大的突變。針對(duì)庫(kù)區(qū)上游河流,本文通過(guò)引入分布函數(shù)的方法可實(shí)現(xiàn)減小此種突變的目的。
本文于一維、二維交接面處引入流速、水位等物理量的二次分布函數(shù),通過(guò)最速下降法對(duì)分布函數(shù)中的系數(shù)進(jìn)行訓(xùn)練,將訓(xùn)練得到的分布函數(shù)作為物理量從一維向二維過(guò)渡的連接,提高了耦合的精確性。桃林口庫(kù)區(qū)上游青龍河流域?qū)嵗?yàn)證結(jié)果表明,模型耦合的精確性得到了很大的提高,且計(jì)算的穩(wěn)定性較高。相比較傳統(tǒng)的重疊投影法,由于增加了最速下降法訓(xùn)練參數(shù)的過(guò)程,計(jì)算的時(shí)間復(fù)雜度有了較大提升,但隨著并行計(jì)算算法和超級(jí)計(jì)算機(jī)的快速發(fā)展,改進(jìn)的重疊投影法適用性會(huì)得到更大的提高。改進(jìn)的模型適用于上中游河道狹窄、下游河面開(kāi)闊,且在一維、二維交接面河道寬度變化較緩的河流發(fā)生的突發(fā)性污染場(chǎng)景。本文選用二次函數(shù)作為各物理量在交接面分布的擬合函數(shù),尋找更精確的分布函數(shù)是下一步的研究方向。