?

基于水動力學模型的城區防洪能力評估研究

2024-03-01 03:27吳濱濱喻海軍馬建明
關鍵詞:北河金堂縣金堂

吳濱濱,喻海軍,3,馬建明,孫 庚,穆 杰

(1.中國水利水電科學研究院,北京 100038;2.水利部防洪抗旱減災工程技術研究中心(水旱災害防御中心),北京 100038;3.水利部京津冀水安全保障重點實驗室 北京 100038;4.應急管理部國家自然災害防治研究院,北京 100085;5.北京中水科工程集團有限公司,北京 100048)

1 研究背景

城區河道現狀防洪能力的準確評估是保障城市防洪安全的基礎[1],也是識別防洪薄弱環節及提出有效工程治理措施的前提。隨著城市規模不斷擴張,橋梁、橡膠壩等河道阻水建筑物的增多,侵占了河道行洪斷面,不同程度地影響了河道行洪能力,需定期對河道防洪能力開展評估[2],以了解河道現狀防洪標準、河道可安全下泄的最大流量以及不同頻率洪水下城區的淹沒范圍和淹沒水深等。特別是多河匯聚、地勢低洼的城市,洪澇災害頻發,洪災損失巨大,迫切需要開展現狀防洪能力評估與提升等方面研究。隨著計算機技術的發展,水動力模型為河道現狀防洪能力的評估提供了有力的工具[3-5]。

已有河道防洪能力研究主要集中在單一河道(干流為主)防洪控制水面線和過流能力模擬計算方面。張景[6]模擬比較了汾河干流下游不同量級洪水下河道沿程水位與堤頂高程的關系,計算了主要城市所屬河段的過流能力,識別了風險點位,并評估了規劃蓄滯洪區的防洪能力提升效果。馮金鵬[2]計算了不同頻率洪水組合下河道典型斷面的特征水位以及保證流量,分析了過流能力較弱的河段。劉一平等[7]評估了水庫下游河道堤防漫溢的臨界流量,并分析了阻水建筑物對河道行洪能力的影響。周惠成等[8]模擬分析了多種工況組合下水庫下游河道的過流能力,探究了糙率對河道過流能力的影響。王宇明[9]復核了柴河水庫下游河道現狀防洪能力,計算了現狀斷面過流能力(保證流量)。方神光等[10]基于現有資料整理、實測水文資料分析和水動力模型相結合的方法,綜合分析了西江干流主要防洪城鎮及堤岸的現狀防洪能力。胡陽等[11]通過二維水動力模型評估了現狀防洪能力,研究了不同洪潮遭遇重現期下兩岸堤防高程差變化對防洪能力的影響。目前針對多河匯聚的城市/區域防洪能力評估開展研究相對較少,僅個別學者開展了河道水面線計算、淹沒范圍和淹沒水深分析、不同工程措施效果提升評估等研究[12-14],尚無針對多河匯聚區域開展過流能力分析等相關研究。

本文以沱江上游金堂縣城區為例,針對城區 “三河匯一江”的復雜交互關系,構建一二維耦合洪水分析模型,評估河道現狀防洪能力(含跨河建筑物)及城區洪水淹沒風險,識別金堂城區防洪的主要薄弱環節,為當地洪水風險管理和水安全保障提供技術支撐。

2 研究區域概況

金堂縣位于成都平原東北部,縣城位于北河、中河、毗河三河與沱江干流匯口處(表1、圖1和圖2),素有“千里沱江第一城”之稱,卻也是沱江流域受洪水威脅最嚴重的地區。城區上游集水面積大(沱江干流三皇廟水文站以上為6590 km2),多條支流呈扇形匯集,其中北河上游的綿遠河、石亭河、湔江發源于四川省三大暴雨區之一的鹿頭山暴雨區,中河、毗河為成都平原排洪河道,汛期暴雨量級大、持續時間長,洪水峰高量大、洪水組成與遭遇復雜;下游狹窄的金堂峽為川西水網僅有的三處出口之一,也是金堂縣城洪水唯一的下泄通道,汛期行洪不暢;加之金堂縣城自身防洪基礎設施薄弱、防洪標準偏低、地勢低洼等因素,一直是洪澇災害頻發區、重發區,已成為成都市乃至四川省防洪體系的顯著短板。

圖1 金堂縣城區位置示意圖Fig.1 Location of urban area of Jintang County

圖2 金堂縣城區放大圖Fig.2 Enlarged view of urban area of Jintang County

表1 金堂境內主要河道特性表Table 1 Characteristics of main rivers in Jintang County

近年來,金堂更是頻繁遭受洪水淹城災害,尤其是2013年、2018年、2020年三場大洪水,人民生命受到威脅、財產遭受巨大的損失、社會經濟發展秩序受到嚴重破壞,洪水災害成為當地發展最大的威脅和破壞因素。其中,沱江三皇廟水文站2013年洪水量級為20~50年一遇;2018年洪水量級為50~100年一遇;2020年洪水量級更是超100年一遇,且金堂縣城經歷了2020.8.12、2020.8.17兩輪洪水過程,三皇廟水文站共發生4次洪峰(12日2:15實測洪峰流量5930 m3/s,22:20實測洪峰流量6600 m3/s;16日22:20實測洪峰流量7950 m3/s;17日14:45實測洪峰流量8100 m3/s)[15],受災人口4.1萬人,轉移安置1620余人,直接經濟損失19.7億元。表2列出了金堂縣城區幾處沿河易受淹低洼地帶的高程及歷史洪水位,可見地面高程約444~447 m,2020.8.17、2018.7.11兩場洪水城區洪水位大都在447~449 m之間,遠高于地面高程。

表2 金堂城區沿河地勢低洼街道高程及歷史洪水位Table 2 Elevation and historical flood level of low-lying streets in the urban area of Jintang County

金堂縣城現有防洪保護體系以堤防、護岸為主,受地形地勢、市政排澇和城市空間布局限制,現有堤防高度均較低,防洪標準較低,難以抵御較大洪水,現有體系已不能滿足城區防洪要求,迫切需要開展防洪研究,識別防洪薄弱環節,為提出有效工程治理措施提供支撐。

3 研究方法

3.1 數據測量和收集本研究中地形及河道斷面數據(含跨河建筑物參數)為2020年1月實測數據。采用無人機搭載LiDAR航攝獲取激光點云數據,并通過航帶拼接、坐標轉換、高程檢查和點云濾波等處理,獲得金堂縣城區范圍內的高精度地形數據(1 m DEM,1 m DSM)。城區外地形數據采用1∶10000地形圖。河道斷面數據采用測深儀、RTK等設備人工測量獲得。

水文數據主要包括三皇廟水文站(沱江干流)和三水水文站(中河)的實時報汛數據、主要閘壩調度情況等,由金堂縣人民政府防汛抗旱指揮部辦公室(以下簡稱金堂縣防汛辦)提供。此外,毗河水位站(2019年設站,僅有2020年汛期數據)的實時報汛數據由四川省都江堰管理局東風渠管理處提供。

2020.8.17、2018.7.11、2013.7.9三場典型洪水的淹沒調查范圍及洪痕調查成果(淹沒水深)由金堂縣防汛辦提供,覆蓋城區重點位置。此外,本研究課題組也開展了2020.8.17、2018.7.11兩場洪水的洪痕調查,主要包括城區及上下游河道沿岸淹沒水位及水深的測量。因此本研究將雙方的調查成果進行相互校驗和復核,獲得最終2020.8.17、2018.7.11兩場洪水的洪痕調查成果。

3.2 洪水分析模型構建采用中國水利水電科學研究院洪水分析軟件IFMS(Integrated(IWHR)Flood Modeling System)中的一二維耦合水動力學模型開展金堂縣城區洪水演進模擬分析,該模型采用基于有限體積法的Godunov格式對一維圣維南方程組和二維淺水方程進行離散求解,已經成功應用于多個區域的洪水模擬,模型原理詳見相關文獻[16-18]。

3.2.1 一維河道模型 為了更加準確地反映金堂縣城區“三河一江”的相互影響以及城區洪水淹沒狀況,一維河道建模范圍將城區及主要河道適當外延,其中北河為綿遠河與石亭江交匯處至北河閘河段,總長約13 km,實測斷面34個;毗河為成都第二繞城高速至中河毗河匯合處河段,總長約13 km,實測斷面138個;中河為旌江干道至三河匯合口河段,總長約20 km,實測斷面129個;沱江干流為北河閘至九龍灘水電站河段,總長約19 km,實測斷面98個。此外,老城區北側北河與中河之間存在一條人工開挖的連通溝渠,枯水期基本干涸,洪水期根據北河和中河來水情況存在水流交互,本研究根據枯水期實測DEM提取溝渠斷面12個,長約0.35 km。

本研究實測跨河建筑物相關尺寸和高程,共32座橋梁和6個閘壩。在模型中,通過斷面概化考慮城區橋墩的阻水作用;對于洪水期間橋面本身也淹沒的情況,在模型中通過減少斷面過水面積來概化橋梁本身的阻水作用。對于城區1、2和3號橡膠壩,北河閘和盤龍寺樞紐,在模型設置中,均認為洪水期間為敞泄,只考慮底板和閘墩的阻水作用。

考慮到2018.7.11、2013.7.9洪水距離本研究實測斷面時間較久,且在此期間開展了毗河清淤、3號橡膠壩改擴建、工農大橋施工等工程,本文根據相關報告和資料[19-20]修正局部河道斷面。

3.2.2 二維地表模型 經現場調研發現金堂城區上游北河和中河存在洪水交互(除人工溝渠外),兩河之間的郊區在歷史洪水中受淹嚴重,因此本研究將二維模型建模范圍擴展至上游郊區,覆蓋三河一江百年一遇洪水淹沒范圍,面積約110 km2。

對整個二維區域進行網格剖分,將區域內公路、鐵路、綠島、韓灘雙島等作為內部約束邊界,最終共劃分成84 514個非結構網格,其中對道路和地形起伏大的區域進行加密處理,靠近河流網格平均尺寸為20 m,而其他區域網格平均尺寸約為50 m。采用高精度DEM數據進行二維高程插值,利用DSM數據提取城區房屋建筑分布,采用面積修正率考慮網格內建筑的影響。同時根據影像數據及城區局部1∶2000地形圖修正2013和2018年地面高程。

3.2.3 一二維耦合 采用側向連接進行一二維耦合設置(見圖3)。將北河、中河、毗河、連通溝渠以及沱江與建模范圍內地面全部進行耦合。

圖3 金堂縣城區一二維耦合模型建模范圍、高程及網格圖Fig.3 The range,elevation and grid of the 1D-2D coupled hydrodynamic model for the urban area of Jintang County

3.2.4 邊界條件 模型上游邊界條件為北河、中河、毗河入境金堂的實際流量過程或設計洪水。對于2020.8.17洪水,中河入流邊界采用三水水文站實測報汛數據,毗河入流邊界采用毗河水位站實測報汛數據推求,北河入流邊界則根據三皇廟水文站、中河入流邊界、毗河入流邊界推求,并通過不斷試算和率定最終確定入流邊界(圖4);對于2018.7.11洪水,中河入流邊界采用三水水文站實測報汛數據,毗河和北河入流邊界采用洪痕調查成果反推斷面洪峰流量,按洪峰流量等比例縮放三皇廟水文站的洪水過程線,并通過不斷試算和率定最終確定入流邊界;對于2013.7.9洪水,采用洪痕調查成果反推斷面洪峰流量,按洪峰流量等比例縮放三皇廟水文站的洪水過程線,并通過不斷試算和率定最終確定邊界。此外,本研究中三皇廟水文站及毗河、北河、中河河口設計洪水采用相關報告中已審批的成果(表3)。

圖4 2020.8.17典型洪水上游邊界條件Fig.4 Upstream boundary conditions of 2020.8.17 flood

表3 本研究采用的設計洪水成果Table 3 The design flood used in this study

模型下游邊界采用九龍灘電站攔河閘水位流量關系曲線(圖5)[23]。該曲線是在堰閘泄流公式計算的基礎上,經過了電站運行中實測水位、流量校正過的水位流量關系曲線,精度較高,滿足本研究要求。

圖5 沱江九龍灘電站攔河閘H~Q關系曲線Fig.5 H~Q relationship curve of Tuojiang Jiulongtan Hydropower Station dam

3.2.5 其他參數設置 本研究構建的一二維耦合模型計算時間步長取2 s,其他時間參數(時間步數、起始時間)根據實際模擬時間段設定。

河道一維模型主槽和河灘地的糙率由實測洪水數據及洪痕調查成果率定獲得。二維地表的糙率設置根據不同的土地利用類型進行初步設置,并經洪痕調查成果率定最終確定。

4 計算結果分析與城區防洪能力評估

4.1 模型率定與驗證本研究采用2020.8.17、2018.7.11、2013.7.9三場典型洪水分別對模型進行率定和驗證。率定的河道主槽糙率值n=0.026~0.035,灘地糙率值n=0.035~0.045,并對樹木和房屋分布較密、橋梁等阻水建筑及跌水明顯的局部斷面適當增大糙率。二維模型中,沱江河灘地糙率取0.035~0.04,農田糙率取0.04,林地糙率取0.04~0.06,道路(含兩側綠化帶)糙率取0.04,建筑糙率取0.06~0.08。

本研究收集了九龍灘電站2018.7.11和2013.7.9兩場典型洪水壩前實測最大水位,2020.8.17由于洪水位超過電站水尺測量范圍,未開展實時測量,本研究組于災后對壩前洪痕進行測量。模擬結果與實測水位或洪痕基本一致,誤差在±0.1 m以內,表明本次分析構建的洪水分析模型可以準確反映九龍灘水電站的實際泄流情況(表4)。

表4 九龍灘電站典型洪水中壩前水位實測值與模擬值對比Table 4 Comparison of measured and simulated water levels at Jiulongtan Hydropower Station

三皇廟水文站在三場典型洪水中總體模擬結果較好(表5和圖6),流量和水位納什效率系數均高于0.85,洪峰誤差均在1%以內,峰現時間誤差均在1.5 h以內,最大水位模擬誤差均在0.2 m以內。其中,2020.8.17洪水模擬精度更高,流量和水位納什效率系數均高于0.95,洪峰誤差僅0.1%,峰現時間誤差僅0.59 h,最大水位誤差僅-0.002 m。

圖6 三皇廟水文站三場典型洪水實測與模擬結果Fig.6 Measured and simulated results of three typical floods at Sanhuangmiao Hydrology Station

表5 三皇廟水文站典型洪水實測值與模擬值對比Table 5 Comparison of measured and simulated flood values at Sanhuangmiao Hydrology Station

2018.7.11、2020.8.17兩場典型洪水的模擬淹沒范圍與金堂縣城實際調查淹沒范圍基本一致(圖7,洪水淹沒范圍僅調查了城區及附近區域,上游郊區未調研),淹沒城區大部分區域、上下游鄉鎮及郊區;2013.7.9洪水模擬城區淹沒范圍跟當地人描述接近,主要北河沿岸老城區受淹嚴重,毗河兩岸局部低洼地段受淹。

圖7 三場典型洪水金堂城區及上下游淹沒調研及模擬情況Fig.7 Investigation and simulation of three typical floods in Jintang urban area and nearby areas

三場典型洪水模擬的最大淹沒水深也基本與實際調研洪痕相符(圖7),大部分點位的模擬淹沒水深誤差絕對值在0.2 m以內,占69%~84%。2020.8.17洪水模擬精度最高,在116個洪痕調查點位中,97個洪痕調查點誤差絕對值均在20 cm以內,約占84%,其中城區模擬精度較高,北河中河上游郊區誤差偏大;2018.7.11洪水模擬精度也較高,在82個洪痕調查點位中,62個洪痕調查點誤差均在20 cm以內,約占76%。2013.7.9洪水調查洪痕較少,僅13個沿河洪痕調查點(橋梁為主),其中9個點位的誤差絕對值在20 cm以內,占69%。模擬誤差較大的點位可能和局部地形概化精度有關,總體上,本研究經率定驗證的模型精度較高,基本可以反映河道及地表洪水的實際演進過程。

4.2 城區河道現狀防洪能力分析采用設計洪水和經率定驗證的模型開展金堂縣城區河段現狀堤防防洪能力評估,其中地形采用2020年1月實測數據,土地利用類型(糙率設置)也為2020年實際情況。由于金堂縣城區“三河匯一江”的復雜交互關系,存在多種洪水組合遭遇情況,本研究采用評估河段與沱江同頻,其他河段相應的方法,分別計算同一洪水頻率下單河段的水面線及淹沒范圍,并將毗河、中河、北河城區段堤頂高程與計算水面線進行比較,評判各河段的現狀防洪能力(堤防超高取0.7 m[24]);再采用包絡線法獲得不同頻率洪水下的淹沒范圍。評估結果見圖8和圖9。

圖8 金堂縣城區現狀防洪能力Fig.8 Current flood control capacity for the urban area of Jintang County

圖9 不同頻率設計洪水下金堂城區淹沒范圍及水深Fig.9 Inundation range and depth of Jintang urban area under different frequency flood

金堂縣現狀城區范圍內已建堤防(或護岸)整體防洪能力薄弱,大部河段防洪能力僅為5~10年一遇防洪標準,僅三條河城區上游局部段及三河匯口左岸滿足20年一遇防洪標準,甚至存在不足5年一遇防洪標準的薄弱段??梢娊鹛贸菂^實際防洪能力較低,迫切需要采取有效工程措施以提高防洪能力。

金堂縣城區2019年末城市面積約28 km2(其中建成區面積約23 km2,水域面積約5 km2)[25]?,F狀設防情況下,當分別遭遇5年、10年、20年、50年、100年一遇洪水時,金堂縣現有城區淹沒面積分別達1.35 km2、3.49 km2、5.72 km2、10.71 km2、12.95 km2,分別約占現有建成區面積的5.87%、15.17%、24.87%、46.57%、56.3%。其中5年一遇設計洪水下,淹沒范圍較小,僅涉及局部低洼區域,包括溫州商業城、生態水城、金海岸、金沙公園銀礫路段、一水廠段、北濱路段、棕櫚湖段等部分區域;10年一遇設計洪水下,現有老城區除了局部地勢相對較高地塊外,大部分均已淹沒;20年一遇設計洪水下,老城區、毗河及北河沿岸低洼地段受淹嚴重;50年、100年一遇設計洪水下,城區約一半區域已經被洪水包圍。

4.3 城區整體過流能力分析目前針對河段過流能力的評估多采用保證流量或者堤防漫溢流量作為河段最大泄量[6-9],本文采用后者作為河段過流能力,計算過程需要多次試算,以得出最接近斷面堤頂的對應流量。通過上一節河道現狀防洪能力分析得到金堂城區主要防洪薄弱河段為不足5年一遇防洪標準的河段,同時也是地勢低洼歷史洪災中淹沒較深的河段,因此本節選取各薄弱河段典型斷面計算不同洪水組合下其過流能力,同時評估金堂城區整體過流能力,以了解城區可以抵御多大的洪水。

金堂城區的總體過流能力受北河、中河、毗河三條河流洪水遭遇組合影響,同時受各條河流的防洪薄弱段限制。因此,有必要將三條河作為整體考慮,分析不同洪水遭遇組合情況下的過流能力。根據2020.8.17、2018.7.11、2013.7.9三場典型洪水的洪水組成,主要存在兩種極端情況:中河和毗河來水多(2018.7.11)、北河來水多(2013.7.9),因此采用這兩場典型洪水組成研究城區整體過流能力(見表6)。

表6 不同洪水組成下金堂城區薄弱河段過流能力Table 6 The flow capacity of weak reach in Jintang Urban area

當按2018.7.11洪水三條河的洪峰組成時(即北河∶中河∶毗河=49.5%∶25.3%∶25.2%),中河和毗河洪水較大,兩者總占比與北河洪水相當,城區總體過流能力相對較小,主要在于中河、毗河薄弱段過流能力小,限制了城區的總體過流能力。當總入流超過4200 m3/s時,城區開始局部受淹(毗河沿岸生態水城、溫州商業城);當總入流超過5000 m3/s時,毗河沿岸金沙公園、金海岸也開始受淹;當總入流超過5400 m3/s時,老城區開始受淹。同時,該洪水組成情況下,中河、毗河來水相對較大,兩者洪水對北河造成一定頂托,導致北河過流能力下降。北河城區薄弱段過流能力僅2700~3000 m3/s;中河、毗河城區薄弱段過流能力約為1100~1200 m3/s;兩河(中河、毗河)交匯河段過流能力約為2300~2400 m3/s。

當按2013.7.9洪水三條河的洪峰組成時(即北河∶中河∶毗河=72%∶14%∶14%),北河洪水占主導,來水較大,城區總體過流能力略有增大,主要在于北河薄弱段過流能力略大于中河、毗河薄弱段過流能力,有利于城區整體過流。當總入流為4700 m3/s時,城區基本不受淹;在4800~5200 m3/s時,僅有兩河交匯處溫州商業城局部被淹;超過5300 m3/s時,老城區也開始受淹。同時,該洪水組成情況下,北河城區薄弱段過流能力略有增大,約為3400~3600 m3/s;中河與北河在上游現有連通渠道處交互,中河城區薄弱段過流能力總體變化不大,約為1200 m3/s;中河、北河洪水均對毗河造成一定程度的頂托,導致毗河城區薄弱段過流能力下降,僅650~1000 m3/s;兩河交匯河段過流能力僅為1800 m3/s。

總體而言,在不同洪水組合下金堂城區整體過流能力較低,僅4200~5000 m3/s左右城區就開始受淹。金堂城區三河相互頂托作用也十分明顯,特別是毗河與中河交匯處角度約60°,兩者水流相互頂沖,出流不暢。因此,金堂城區防洪能力的提升要在最大化挖潛城區河道自身防洪能力(河道擴卡、疏挖整治、完善堤防等)的基礎上,開展上游蓄滯洪工程或者分洪工程,以減小城區段洪水流量及三河相互頂托作用。

4.4 城區跨河建筑物防洪能力評估金堂縣2035規劃城區范圍內跨河建筑物包括橋梁19座、閘壩2座。由于橋梁建設年代及設計標準各有不同,其過流能力亦有所差別(見表7)。大部分橋梁特別是近年來新建的橋梁基本滿足20年一遇的行洪要求,部分二十世紀八十、九十年代左右建設的橋梁由于年代久遠,設計標準較低,而不能滿足行洪要求。其中平安橋不足5年一遇防洪標準,其橋底高程低于5年一遇的洪水水面線0.35 m;毗河大橋不足10年一遇防洪標準,其橋底高程低于10年一遇的洪水水面線0.26 m;康寧橋、綠洲橋、綠島鳥橋、綠島北橋、毗河人工島連接橋、北河三橋不足20年一遇防洪標準。此外,北河大橋、清風橋等拱橋最高拱頂雖然高于20年一遇洪水水面線,但橋身兩側阻水依然比較嚴重。此外,北河閘不足20年一遇防洪標準,盤龍寺閘大于20年一遇防洪標準。因此,在金堂城區防洪能力提升措施中建議結合堤防改造工程,對平安大橋、毗河大橋進行改建或新建。

表7 跨河建筑物現狀防洪能力Table 7 Current flood control capacity of cross-river buildings

5 結論

(1)本文構建的一二維耦合洪水分析模型,經三場典型洪水實測及調查數據率定和驗證,結果表明模型總體精度較高,其中三皇廟水文站洪峰誤差均在1%以內,峰現時間誤差均在1.5 h以內,建模范圍內大部分洪痕點位的模擬最大淹沒水深誤差在0.2 m以內。

(2)金堂縣城區范圍內已建堤防(或護岸)整體防洪能力薄弱,大部分河段防洪能力僅為5~10年一遇防洪標準,僅三條河城區上游局部段及三河匯口左岸滿足20年一遇防洪標準,甚至存在不足5年一遇防洪標準的局部薄弱段。

(3)不同洪水組合下金堂城區整體過流能力較低,僅4200~5000 m3/s左右城區就開始受淹;城區三河相互頂托作用十分明顯,當北河洪水占主導時,城區整體過流能力相對較大。金堂城區防洪能力的提升建議在最大化挖潛城區河道自身防洪能力的基礎上,開展上游蓄滯洪工程或分洪工程,以減小城區段洪水流量及三河相互頂托作用。

(4)城區多座跨河建筑物(橋梁和閘)防洪能力不足20年一遇防洪標準,其中平安橋、毗河大橋分別不足5年、10年一遇防洪標準,在未來金堂城區防洪能力提升中建議結合堤防改造工程,對平安大橋、毗河大橋進行改建或新建。

猜你喜歡
北河金堂縣金堂
“五老”精神
提升技能促進就業 打造“金堂焊工”勞務品牌
金堂縣人力資源形勢分析與研判
四川省金堂縣職業高級中學:多元融合發展 增強技能支撐 加速培養人才
金堂縣:“春風行動”激活創業就業“一池春水”
晉源區北河下村:甜蜜產業帶來幸福生活
北 河
北河
金堂大學生的消費情況及帶來的經濟影響
Numerical simulation of cementing displacement interface stability of extended reach wells *
91香蕉高清国产线观看免费-97夜夜澡人人爽人人喊a-99久久久无码国产精品9-国产亚洲日韩欧美综合