• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    基于二維淺水方程的胖頭泡蓄滯洪區(qū)洪水演進數(shù)值模擬

    2024-01-01 00:00:00孫振宇周運浩張勤旭秦夢恩張明亮
    人民珠江 2024年6期

    摘要:探究蓄滯洪區(qū)洪水演進對蓄滯洪區(qū)的防洪工作有著重要意義。為在洪水來臨時,短時間內(nèi)對洪水運動進行預(yù)測。運用 Godunov 型有限體積法建立了1個二維淺水水動力數(shù)值模型,該模型的的優(yōu)勢在于采用非結(jié)構(gòu)性三角形網(wǎng)格,能夠較好地擬合不規(guī)則邊界,并采用 GPU 異構(gòu)并行技術(shù),有效解決了非結(jié)構(gòu)化網(wǎng)計算耗時長的問題,運用限制器減少二階精度數(shù)值的擴散和振蕩,引入負水深和干濕邊界處理技術(shù),使模型模擬的精確度得到了有效提高。模型應(yīng)用多個潰壩水流實驗算例對模型進行率定,在驗證其具有良好的干濕邊界計算能力和靜水和諧性的基礎(chǔ)上,采用嫩江流域1998年洪水數(shù)據(jù)對嫩江流域的胖頭泡蓄滯洪區(qū)的洪水演進進行模擬計算。結(jié)果顯示:模型對胖頭泡蓄滯洪區(qū)模擬的蓄水量為34.9億 m3,與實測值的相對誤差為1.7%,洪水淹沒面積為1267.27 km2,與實測值的相對誤差為9.24%。對研究區(qū)域不同的土地類型進行糙率分區(qū)的情況下,模型模擬的蓄水量結(jié)果更趨近于實測結(jié)果。并行化后的程序能使模型計算速度提升33%,因此模型能快速、精準計算復(fù)雜地形潰決水流的洪水演進過程。研究結(jié)果可以為蓄滯洪區(qū)的防洪規(guī)劃、災(zāi)情預(yù)警提供相關(guān)信息。

    關(guān)鍵詞:淺水方程;有限體積法;蓄滯洪區(qū);洪水演進

    中圖分類號:TV133 文獻標識碼:A文章編號:1001-9235(2024)06-0082-10

    Numerical Simulation of Flood Evolution in Pangtoupao Flood Storage Area Based on Two-Dimensional Shallow Water Equation

    SUN Zhenyu1, ZHOU Yunhao1, ZHANG Qinxu1, QIN Mengen1, ZHANG Mingliang1,2

    (1. School of Ocean Science and Environment, Dalian Ocean University, Dalian 116023, China;2. Liaoning Coastal EcologicalEnvironment and Disaster Protection Engineering Technology Innovation Center, Dalian 116023, China)

    Abstract: Exploring the flood evolution in flood storage areas is of great significance for flood control work in flood storage areas. To predict the dynamics of incoming floods within a short period of time, The Godunov-type finite volume method was used to establish a two-dimensional hydrodynamic numerical model of shallow water. The model had the advantages of adopting unstructured triangular mesh, which could better fit irregular boundaries, and it adopted GPU heterogeneous parallel technology, which effectively solved the problem of slow computation speed of an unstructured network. The model used a limiter to reduce the diffusion and oscillation of second-order precision values. The accuracy of the model simulation was improved effectively by using negative water depth and dry and wet boundary processing techniques. The model was calibrated using multiple dam-break flow experiments and validated with an excellent capability in calculating dry and wet boundaries and static water balance. Flood data from 1998 in the Nenjiang River Basinwas used to simulate and calculate flood evolution in the Pangtoupo flood storage area of the Nenjiang River Basin. The results show that the model simulates a storage capacity of 34.9 billion m3 for the Pangtoupo flood storage area, with a relative error of 1.7% compared to the measured value, and a flooded area of 1267.27 km2, with a relative error of 9.24% compared to the measured value. When different land types in the study area are classified based on roughness coefficients, the simulated storage capacity results of the model approach to the measured values. The parallelized program improves the model′s calculation speed by 33%. The model thus can rapidly and accurately calculate the flood evolution process of dam-break flows in complex terrain. The results provide relevant information for flood control planning and disaster warning in flood storage areas.

    Keywords: shallow water equation; finite volume method; flood storage area; flood evolution

    洪水一直以來都是全球最主要的自然災(zāi)害之一,目前很多專家和學者依據(jù)水動力學原理建立了一維和二維水動力模型,用于在洪水來臨時提供有效的預(yù)報和預(yù)警。在研究一維水動力模型方面,大部分學者[1-4]是對河道、河網(wǎng)和水庫等區(qū)域進行洪水演進模擬和驗證,但在研究大尺度水體區(qū)域過程中,如湖泊、蓄滯洪區(qū)、河口等,其水深、流動速度等要素在水平方向的改變會遠大于垂直方向上的改變,所以僅用一維模型很難精準的模擬洪水演進過程;若采用三維模型模擬,雖然可以更真實的反映水流運動特性,但對于垂向沿程變化不明顯的水體,例如陸面洪水等,其垂直方向上的水深、流速等各要素值接近平均分布[5],所以使用三維模型計算會較為復(fù)雜且計算速度效率將大幅降低,因此目前多采用二維模型對研究區(qū)域進行模擬計算?,F(xiàn)如今對二維淺水方程組的二維模型研究和應(yīng)用越來越廣泛,因其在洪水模擬中的計算效率以及計算精度遠高于簡化模型,可以更精準地描述洪水的運動情況[6-7],尤其是基于 Godunov 型有限體積法建立的淺水方程模型,其具有良好的捕捉動邊界的能力,在模擬洪水運動中的瞬時陸地流量等方面具有很大的潛力[8-9]。目前眾多國內(nèi)外學者[10-16]都基于有限體積法對二維淺水方程模型進行改進,用于洪水演進數(shù)值模擬當中。平面二維淺水方程的應(yīng)用可以較為恰當?shù)孛枋鰷\水水流運動,在滿足模擬精度的同時,計算效率也能得到相應(yīng)的提高。由于洪水水流其流速快以及強間斷性的特點,模型更適合于采用顯格式計算,而顯格式計算受到時間步長及網(wǎng)格精度的限制;在研究地形不規(guī)則的研究區(qū)域中,相對于結(jié)構(gòu)化網(wǎng)格而言,非結(jié)構(gòu)化網(wǎng)格具有更好的擬合邊界的能力,網(wǎng)格質(zhì)量相對更好,但對計算機的性能要求以及計算成本較高,需要耗費更長的時間。

    基于非結(jié)構(gòu)性三角形網(wǎng)格,并運用 GPU異構(gòu)并行技術(shù)解決運用非結(jié)構(gòu)化網(wǎng)格模型計算緩慢問題。使用 Roe 格式計算界面通量,運用 Godunov 有限體積法建立出1個二維淺水動力模型,并且引入了負水深和干濕邊界處理技術(shù),應(yīng)用模型模擬多個潰壩水流實驗室算例,對模型進行率定,驗證模型穩(wěn)定性以及準確性。利用衛(wèi)星遙感技術(shù)識別胖頭泡蓄滯洪區(qū)的地貌景觀特征,對不同的區(qū)域設(shè)置其相應(yīng)的糙率值,從而保證模型的精準度,模擬在相同或不同的糙率條件下,洪水在蓄滯洪區(qū)的運動狀況及洪水的衰減作用。

    1水動力模型

    1.1二維淺水控制方程

    遵循淺水Boussinesq假設(shè)與靜壓假定,對 Navier-Stokes方程沿垂直方向流速進行深度平均近似,在忽略風應(yīng)力、科氏力以及和二階擴散項的情況下,二維淺水控制方程最終表達式見式(1)、(2)。

    + + = S(1)

    式中:U 為守恒性變量的向量;F 為 x 的對流通量;G 為 y 的對流通量;S 為源項。

    ?0?

    U =( v(u)),F(xiàn) =( v(2)),G =( 2(v)),S = -- g(g)h(h) --τ(τ)by(bx)(2)

    式中:h 為水深;u 為 x 方向的流速,v 為 y 方向的流速;τbx 為 x 方向的底部摩阻項,τby 為 y 方向的底部摩阻項;g 為重力加速度;η為水位;n 為曼寧系數(shù)。其表達式見式(3)。

    τbx = gn2 uu2 + v2, τby = gn2 vu2 + v2(3)

    1.2有限體積離散

    用非結(jié)構(gòu)性三角網(wǎng)格對研究區(qū)域進行劃分,用 Godunov 型有限體積法對控制方程進行數(shù)值離散和求解,使淺水方程在非結(jié)構(gòu)三角網(wǎng)格上進行積分可得式(4)。

    dv +?·Eadv dv =?·Edif dv + Sdv(4)

    式中:E adv 為對流通量;Eaif為擴散通量。運用 Green 公式將方程進行積分整理后可得式(5)。

    ΔUi = (E ij(*)?nij)?lij +Sdv(5)

    式中:Ai 為第i個單元的面積;m 為單元邊的個數(shù); E*ij·nij為第i個單元第j 條邊的通量;lij為單元每條邊的長度;Vi 為控制單元;Ui 為控制單元平均值,其表達式見式(6)。

    Ui =Udv(6)

    1.3動量通量計算

    由于非結(jié)構(gòu)網(wǎng)格在每個控制單元邊界存在有間斷現(xiàn)象,從而構(gòu)成 Riemann 問題。本文采用 Roe 格式求解 Riemann 界面法向數(shù)值通量問題[17],數(shù)值通量表達式見式(7)。

    E* ? n = (F,G )R? n +(F,G )L? n -| J(?)|(UL -

    UR )](7)

    式中:J(?)為 Roe 平均的 Jacobian 矩陣;UL、UR 為單元邊界兩側(cè)的 Riemann 守恒性變量。

    1.4動邊界處理及負水深處理

    干濕界面是指在計算區(qū)域內(nèi)有水和無水交替變化或部分淹沒的區(qū)域,即動邊界。為了真實的模擬在天然河流,河口以及海灣處發(fā)生的洪水演進、潰壩水流演進等現(xiàn)象,從而引入干濕界面處理技術(shù)。針對復(fù)雜地形引入了干濕界面處理技術(shù)[18],以h0為最小限制水深,若水深低于 h0則為干邊界;高于h0則為濕邊界。圖1為4種干濕邊界情況。

    根據(jù)干濕邊界的4種情況判斷,可以把計算網(wǎng)格分為3種網(wǎng)格:網(wǎng)格的3個頂點全為濕點,且3個邊界全為濕邊界或半濕邊界組成,則此網(wǎng)格為濕網(wǎng)格;網(wǎng)格的3個邊界全為干邊界或者半濕邊界組成,則此網(wǎng)格為干網(wǎng)格;除濕網(wǎng)格與干網(wǎng)格外其余均為半濕網(wǎng)格。

    基于淺水方程的模擬過程當中,邊界通量的計算是由時間步長、單寬通量及邊界長度決定。如果時間步長設(shè)置過小,會使模型計算速度大幅降低。在較小水深的情況下,時間步長如果設(shè)置過大,會導(dǎo)致計算不穩(wěn)定,并且會造成網(wǎng)格邊界的通量高于上游網(wǎng)格的水量,從而導(dǎo)致上游網(wǎng)格形成負水深,最后導(dǎo)致質(zhì)量不守恒。針對負水深問題,在計算過程中記錄下負水深單元網(wǎng)格 Ai,并將這些單元網(wǎng)格的水深調(diào)整為 h0(最小限制水深),流速調(diào)置為0,之后將單元網(wǎng)格Ai 周圍的共邊網(wǎng)格中高于最小限制水深的水量補入單元網(wǎng)格 Ai 當中,直至負水深成為最小限制水深。

    1.5并行化處理

    由于二維水動力模型數(shù)據(jù)復(fù)雜,計算速度緩慢,因此本模型采用 GPU異構(gòu)并行技術(shù)實現(xiàn)高速運算。在計算過程中,CPU 作為主要領(lǐng)導(dǎo),負責產(chǎn)生和交付多線程任務(wù),而 GPU 作為計算任務(wù)的執(zhí)行者,將任務(wù)分配至 GPU 硬件的計算單元中進行計算,兩者進行必要的信息傳遞,完成整個計算流程。OpenACC是一種高級編程模型,可通過指令來加速C/C++和 Fortran語言中的并行部分。其原理是先在主機端的加速器上分配空間,然后將數(shù)據(jù)從 CPU 的內(nèi)存復(fù)制到 GPU 的內(nèi)存,并發(fā)送相應(yīng)的代碼。這些數(shù)據(jù)和代碼在加速器上排隊等待空閑資源來進行加速計算,其并行過程主要包括3個部分:首先是對變量參數(shù)、網(wǎng)格信息和流量等數(shù)據(jù)進行初始化;其次是通量計算和項源求解;最后輸出網(wǎng)格單元的水深、流速等信息,模型計算框架見圖2。在本模型中,使用 kernel、loop、private 和 routine 指令來完成并行計算,模型中每個附入代碼的區(qū)域都被注釋(! acc)來控制在 GPU 上的并行計算。在內(nèi)部循環(huán)計算過程中,使用 kernels loop 組合指令將計算區(qū)域加載到 GPU 上,成為一個能執(zhí)行的加速函數(shù),使用! acc routine seq 指令用來分配加速設(shè)備上的內(nèi)存。由于本模型使用的是顯格式求解控制方程,所以在網(wǎng)格計算過程中,每個單元網(wǎng)格的計算只需前一步時間步長的結(jié)果,因此使用OpenACC模式的并行算法可以實現(xiàn)對本模型的并行計算。

    2數(shù)值模擬驗證與應(yīng)用

    2.1在 90°彎道下的潰壩面模擬

    在潰壩模型應(yīng)用中,特別是大尺度區(qū)域進行潰壩水流數(shù)值模擬時,可能會遇到在山區(qū)河道中出現(xiàn)急彎的情況,為了研究本模型是否能處理地形存在突變的潰壩水流,模型模擬了上游矩形水庫潰壩實驗,模型示意見圖3。上游水庫幾何尺寸為244 cm×239 cm,下游河道為寬49.5 cm 的 L 型矩形斷面,河道上端長約400 cm,下端長約300 cm,下游河道高程比上游水庫高程高33 cm。河道與水庫內(nèi)無坡度,河道出口端為自由出口,試驗前河道設(shè)置為濕河床,水庫初始水位高于河道20 cm,初始底部粗糙值設(shè)置為0.012。試驗時將河道與水庫的連接閘門快速打開從而產(chǎn)生潰壩水流。圖4為在不同測點(G1—G6)水位的實測值和模擬值隨時間的變化圖像。

    本算例總計算時長為40 s,圖中 G1點為上游水庫當中的測點,由圖4a可以看到其測點的水位會隨著閘門打開呈現(xiàn)持續(xù)降低的趨勢,水位下降的幅度隨著時間的增加而逐步減少。G2、G3和 G4測點都位于河道上端,都呈現(xiàn)出了一種水位突然增大的現(xiàn)象,其原因是下游彎道對水波的反射作用所造成的,其水位均是到達極大值后經(jīng)過一段時間再突然增大到最大值,之后再緩慢降低;由圖4b—4d 可看出測點 G4的水位最大值出現(xiàn)得較早,其原因是 G4點相對于 G2和 G3測點離 L 形轉(zhuǎn)彎處的距離更近。 G5、G6均在彎道后方的河道下游,由于水流受到彎道的調(diào)整作用,其水位的變化過程相對于其他測點來說更平緩了些,模型計算結(jié)果與實測結(jié)果大致相同。從整體結(jié)果來看,本模型能夠模擬類似于這種地形存在突變的潰壩水流實驗。

    2.2過駝峰的潰壩波傳播模擬

    本算例是為了驗證模型的穩(wěn)定性,以及動態(tài)邊界處理的合理性與有效性。計算區(qū)域為75 m×30 m 的矩形封閉水槽,其中大壩位于 x=16 m 處,大壩的厚度忽略不計。水槽底部的曼寧系數(shù)為0.018,四周采用固壁邊界。水槽有3個駝峰地形,其中2個高度為1 m 的低駝峰分別位于(30 m,6 m)和(30 m,24 m),一個高度為3 m 的高駝峰位于(47.5 m,15 m)。駝峰高程的表達式見式(8)。

    試驗時,大壩上游初始水位為1.875 m。試驗?zāi)M計算了300 s 內(nèi)水流運動情況。t=2、6、12、24、30、300 s 的二維水深結(jié)果見圖5。由結(jié)果可看出,水流特征具有良好的對稱性,符合水流運動規(guī)律。當潰壩波傳播時高駝峰未曾被水淹沒,但高度為1 m 的低駝峰曾被水完全淹沒,低駝峰的漲水和退水現(xiàn)象比較明顯。t=2 s 時,水流通過小駝峰,并處于漲水狀態(tài);t=6 s 時,小駝峰已被完全淹沒,潰壩波傳播至高峰處;t=12 s 時,洪水繞高峰兩側(cè)流過,且繞流現(xiàn)象較為明顯;t=24、30 s 時,洪水已完全淹沒平底區(qū)域,水槽末端固壁邊界對水流有明顯的反射現(xiàn)象,t=300 s 時,由于潰壩波與邊界相互作用和能量的消散,研究區(qū)域內(nèi)的水流幾乎成靜止狀態(tài)。

    2.3胖頭泡蓄滯洪區(qū)洪水演進

    胖頭泡蓄滯洪區(qū)作為哈爾濱市重要的防洪體系之一,其作用是承擔松花江哈爾濱市發(fā)生巨大洪水時的防洪工作,從而減輕松花江的防洪壓力。利用1998年胖頭泡蓄滯洪區(qū)的歷史資料,使用本模型對1998年大洪水中胖頭泡蓄滯洪區(qū)的洪水演進過程進行了數(shù)值模擬,其中還包括有網(wǎng)格分布、地形插值、分區(qū)域糙率設(shè)置、邊界的條件以及初始條件的設(shè)置等。從淹沒面積、淹沒最大水深,蓄水量等方面,將模擬結(jié)果與1998年洪水數(shù)據(jù)進行對比,驗證模型的合理性。

    2.3.1網(wǎng)格敏感度分析

    相對于結(jié)構(gòu)化網(wǎng)格,非結(jié)構(gòu)化網(wǎng)格能更好的模擬不規(guī)則邊界,且網(wǎng)格的質(zhì)量更好,但非結(jié)構(gòu)網(wǎng)格在計算速度上會低于結(jié)構(gòu)化網(wǎng)格,使其在使用過程中對計算機的性能要求較高,且計算的成本也較大。因此本文在保證模型計算精準度的前提下,為了提高模型計算速度;節(jié)約計算成本,本文對網(wǎng)格的收斂性進行了分析。最終選出網(wǎng)格數(shù)量分別為13240、16877、20303個的3套網(wǎng)格,分別對其地理位置相同的4個測點的洪水水深和洪水淹沒面積進行分析比對,4個測點地理位置見圖6,淹沒水深變化圖和洪水淹沒面積見圖7、8。從結(jié)果來看,隨著網(wǎng)格數(shù)量的增加,模型計算時間也隨之增加,模擬時長為200 h 時,13240、16877、20303個的網(wǎng)格所用計算時間分別需為54、74、90 min,但其水深及淹沒面積的計算結(jié)果并未發(fā)生較大變化。因此在不同網(wǎng)格精度對模擬結(jié)果影響結(jié)果不大的前提下,選擇計算速度較快且計算成本較低的13240個網(wǎng)格完成模擬計算。

    2.3.2變化曼寧系數(shù)分析

    糙率系數(shù)是指研究區(qū)域中各種地物對水流阻力產(chǎn)生作用的1個綜合系數(shù),因此糙率系數(shù)的取值要根據(jù)研究區(qū)域的地形地貌、各地區(qū)土地實際利用情況、植被覆蓋情況及洪水泛濫的季節(jié)等進行初步分析,從而劃定不同土地類型的糙率取值。

    由于胖頭泡蓄滯洪區(qū)的土地情況較為復(fù)雜,故本文將地理信息系統(tǒng)(GIS)和遙感影像處理等技術(shù)應(yīng)用于胖頭泡蓄滯洪區(qū)洪水淹沒的研究當中。根據(jù)衛(wèi)星遙感影像,可將研究區(qū)域大致分為四類:水體、濕地、植被和人工建筑地區(qū)。經(jīng)過查證,胖頭泡蓄滯洪區(qū)植被主要為玉米等農(nóng)作物,而位于東北的玉米一般在每年8月份收獲,根據(jù)《1998年松花江大洪水》的記載,洪水發(fā)生的時間為八月至九月初。因此,當洪水發(fā)生時玉米正處于成熟階段,故在本文中將植被的糙率值設(shè)置為0.065。其他區(qū)域的糙率值參考了實際情況及文獻[19-22],設(shè)置如下:水體:0.03;濕地:0.04;人工建筑區(qū):0.09。由于裸地面積相對于大尺度的研究區(qū)面積極小,對研究區(qū)域洪水演進影響也非常小,同時考慮到小區(qū)域土地類型特征的設(shè)置無法在大尺度網(wǎng)格上得到滿足,故在此次模擬中忽略了裸地糙率值的設(shè)定。衛(wèi)星影像及土地利用情況見圖9。為了探究在胖頭泡蓄滯洪區(qū)設(shè)置全域糙率和分區(qū)糙率對模型模擬的影響作用,根據(jù)其地貌特點結(jié)合《松花江流域防汛資料匯編》對全域糙率設(shè)置為 n=0.055和 n=0.065。將模擬結(jié)果與糙率分區(qū)模擬結(jié)果結(jié)合《1998年松花江大洪水》中提供的分洪和退洪數(shù)據(jù)推導(dǎo)出的各個時刻蓄水量的結(jié)果進行對比,其對比結(jié)果見圖10。結(jié)果顯示,在糙率分區(qū)的情況下,模擬洪水在340 h 的蓄水量為34.9億 m3,與《1998年松花江大洪水》記錄的35.5億 m3較為接近,相對誤差為1.7%。

    2.3.3胖頭泡洪水演進模擬分析

    本模型以1998年松花江實測洪水資料為依據(jù),在糙率分區(qū)的情況下對胖頭泡蓄滯洪區(qū)進行數(shù)值模擬計算。由資料可知,本次洪水決口位于大安市上游、蓄滯洪區(qū)國營農(nóng)場附近的嫩江左岸堤防處附近,其流量過程及決口位置見圖11。本次模擬總時長為340 h,計算結(jié)果選取100、150、200和340 h 4個時刻的洪水淹沒范圍以及流場分布圖來演示洪水的演進過程,胖頭泡蓄滯洪區(qū)地形高程(Z)變化及流場分布見圖12。從圖12中可以看出,洪水從胖頭泡決口進入蓄滯洪區(qū)后,分洪初期洪水向著東和東北方向流入低洼地區(qū),后期洪水主要向著南和東南方向演進,洪水演進過程符合胖頭泡蓄滯洪區(qū)西北高、東南低的地勢地貌特征。在模擬過程中,串行計算結(jié)果與并行計算結(jié)果基本一致,在使用并行計算時,計算效率提升了33%。不同時刻淹沒面積計算結(jié)果與其他文獻[22]中的 MIKE、EFDC模型的計算結(jié)果對比圖13,以《1998年松花江大洪水》中洪水調(diào)查量算的淹沒面積為依據(jù),3種模型模擬淹沒面積的結(jié)果與實測數(shù)據(jù)的相對差值見表1。從結(jié)果數(shù)據(jù)來看,本模型與實測值1160 km2的相對誤差為9.24%,在誤差的允許范圍之內(nèi),故本模型計算結(jié)果合理。模型計算所得的最大水深分布面積結(jié)果,與其他文獻[22]中 MIKE、EFDC模型的計算結(jié)果對比結(jié)果見圖14。圖13、14可知,二維淺水控制方程模擬結(jié)果與王靜[22]的模型對1998年大洪水進行建模計算的結(jié)果比較一致??傮w來說,模型能夠模擬復(fù)雜地形上的洪水運動過程,可以為復(fù)雜地形的洪水演進進行預(yù)測。

    3結(jié)論

    基于有限體積法以 Roe格式求解 Riemann界面法向數(shù)值通量問題,模型采用了干濕邊界處理以及負水深處理技術(shù),進一步提高了模型地穩(wěn)定性以及靜水和諧性;在非結(jié)構(gòu)三角網(wǎng)格控制體上用有限體積法對方程進行離散,使模型能夠較好的擬合類似蓄滯洪區(qū)等復(fù)雜地形及不規(guī)則邊界,通過引入OpenACC模式的并行算法,成功實現(xiàn)了對該模型的加速計算,使模型計算效率提升了33%。將模型應(yīng)用于胖頭泡蓄滯洪區(qū)洪水演進數(shù)值模擬后發(fā)現(xiàn):洪水由潰口流入胖頭泡蓄滯洪區(qū)后,洪水先向著東和東北方向流入低洼地區(qū),有部分洪水順勢向南和東南方向順勢而下,后期洪水主體主要向東南和南演進,到達排水口后有少部分洪水向西擴散,分洪212 h后洪水由排水口回歸嫩江,這與胖頭泡蓄滯洪區(qū)西北高東南低的地勢相符合,其洪水流量與淹沒面積呈正相關(guān)。在分洪340 h 的過程當中,模型模擬洪水淹沒面積和蓄水量的結(jié)果與測量數(shù)據(jù)基本符合,且洪水演進過程中的流場分布合理。此模型可以對大區(qū)域復(fù)雜地形的洪水運動進行快速預(yù)測,其結(jié)果具有可靠性、合理性,本模型可以為蓄滯洪區(qū)的防洪工作、預(yù)報預(yù)警提供科學的依據(jù)。

    參考文獻:

    [1]季益柱,丁全林,王玲玲,等.三峽水庫一維水動力數(shù)值模擬及可視化研究[J].水利水電技術(shù),2012,43(11):21-24.

    [2]劉恒恒,齊鵬云,萬能勝,等.巢湖閘下河網(wǎng)洪水演進及閘門調(diào)度數(shù)值模擬[J].人民長江,2022,53(5):41-46.

    [3]胡嘉鏜,李適宇.珠江三角洲一維鹽度與三維斜壓耦合模型[J].水利學報,2008(11):1174-1182.

    [4]王盼,何洋,杜志水.基于水動力數(shù)值計算的城市設(shè)計洪水模擬研究[J].西安理工大學學報,2020,36(3):362-366.

    [5]張明亮.近海及河流環(huán)境水動力數(shù)值模擬方法與應(yīng)用[M].北京:科學出版社,2015.

    [6] FAN Y Y,AO T Q,YU H J,et al. A coupled 1D-2D hydrodynamic model for urban flood inundation[J]. Advances in Meteorology,2017. DOI:10.1155/2017/2819308.

    [7] GUAN M,WRIGHT N G,SLEIGH P A. A robust 2D shallow water model for solving flow over complex topography using homogenousfluxmethod [J].InternationalJournal for Numerical Methods in Fluids,2013,73(3):225-249.

    [8] OALBARRIEAT M,MEDINA R,GONZALEZ M,et al. A finite volume-finite difference hybrid model for tsunami propagation and run up[J].Computers and GeoSciences,2010,37(8):1003-1014.

    [9]張明亮.濱海鹽沼濕地退化機制及生態(tài)修復(fù)技術(shù)研究進展[J].大連:大連海洋大學學報,2022,37(4):539-549.

    [10]LI D M,ZHEN Z,ZHANG H Q,et al. Hydrodynamic model of Daya Bay based on finite element method[J]. Environmental Earth Sciences,2020,79(1). DOI:10.1007/S12665-020-09019-X.

    [11]GALLEGOS H A,SCHUBERT J E,SANDERS B F. Two- dimensional, high-resolutionmodelingofurbandam-break flooding:A case study of Baldwin Hills,California[J]. Advances in Water Resources,2009,32(8):1323-1335.

    [12]陳海鑫.二維淺水流動數(shù)值模擬研究及其應(yīng)用[D].大連:大連理工大學,2013.

    [13]LIANG Q H. Flood simulation using a well-balanced shallow flow model[J]. Journal of Hydraulic Engineering,2010,136(9):669-675.

    [14]宋利祥,周建中,王光謙,等.潰壩水流數(shù)值計算的非結(jié)構(gòu)有限體積模型[J].水科學進展,2011,22(3):373-381.

    [15]張大偉,程曉陶,黃金池,等.基于 Godunov格式的潰壩水流數(shù)學模型[J].水科學進展,2010,21(2):167-172.

    [16]侯精明,張兆安,馬利平,等.基于 GPU 加速技術(shù)的非結(jié)構(gòu)流域雨洪數(shù)值模型[J].水科學進展,2021,32(4):567-576.

    [17]冀永鵬,張洪興,王運濤,等.基于二維淺水方程的城市地面洪水演進數(shù)值模擬研究[J].水資源與水工程學報,2020,31(2):42-49,56.

    [18]HUANG Y X,ZHANG N C,PEI Y G. Well-balanced finite volume scheme for shallow water flooding and drying over arbitrary topography[J]. Engineering Applications of Computational Fluid Mechanics,2013(7):40-54.

    [19]王秀杰,胡冰,苑希民,等.洪水與風暴潮共同作用下的潰堤洪水一維、二維耦合模型及應(yīng)用[J].南水北調(diào)與水利科技,2017,15(5):43-49.

    [20]李大鳴,羅珊,范麗虹,等.西三洼洪水演進數(shù)值模擬及洪水風險分析[J].水利水電技術(shù),2018,49(8):78-86.

    [21]肖玉紅,李大鳴,白玲,等.蓄滯洪區(qū)洪水演進過程數(shù)值模擬及洪災(zāi)損失分析[J].中國農(nóng)村水利水電,2012(12):143-147,149.

    [22]王靜.基于 EFDC 的胖頭泡蓄滯洪區(qū)洪水演進模擬研究[D].大連:大連理工大學,2017.

    (責任編輯:李燕珊)

    日日啪夜夜撸| 好男人视频免费观看在线| 免费人妻精品一区二区三区视频| 国产精品一区二区三区四区免费观看| 午夜激情福利司机影院| 91久久精品国产一区二区三区| 久久精品国产亚洲av涩爱| 成年av动漫网址| 欧美成人a在线观看| 国产伦精品一区二区三区视频9| 夫妻性生交免费视频一级片| 亚洲美女搞黄在线观看| 国产免费一区二区三区四区乱码| 亚洲av.av天堂| 夜夜看夜夜爽夜夜摸| 夜夜骑夜夜射夜夜干| 免费播放大片免费观看视频在线观看| 91精品一卡2卡3卡4卡| 日本一二三区视频观看| 一级片'在线观看视频| 纵有疾风起免费观看全集完整版| 在线观看一区二区三区激情| 五月天丁香电影| 成人特级av手机在线观看| 亚洲av中文av极速乱| 日韩欧美一区视频在线观看 | 哪个播放器可以免费观看大片| 全区人妻精品视频| 久久这里有精品视频免费| 晚上一个人看的免费电影| 欧美3d第一页| 日韩强制内射视频| 国产一区二区在线观看日韩| 人体艺术视频欧美日本| av在线app专区| 欧美丝袜亚洲另类| 久久97久久精品| 日韩av免费高清视频| 日日啪夜夜撸| 国产白丝娇喘喷水9色精品| 亚洲av在线观看美女高潮| 日本爱情动作片www.在线观看| 午夜日本视频在线| 91狼人影院| 久久久久久人妻| 国产精品久久久久久av不卡| 亚洲内射少妇av| 男女下面进入的视频免费午夜| 国产在线免费精品| 精品一区二区三卡| 美女视频免费永久观看网站| 亚洲欧美清纯卡通| 黄色怎么调成土黄色| 亚洲av男天堂| av播播在线观看一区| 欧美日韩一区二区视频在线观看视频在线| 欧美日本视频| 极品少妇高潮喷水抽搐| a 毛片基地| 在线观看免费高清a一片| 少妇的逼好多水| 欧美日韩精品成人综合77777| 久久人妻熟女aⅴ| 免费大片黄手机在线观看| 亚洲av中文字字幕乱码综合| 久久精品国产a三级三级三级| 下体分泌物呈黄色| 一区二区三区四区激情视频| 日本爱情动作片www.在线观看| 2018国产大陆天天弄谢| 深夜a级毛片| 成人亚洲精品一区在线观看 | 少妇 在线观看| 欧美成人一区二区免费高清观看| 黄色怎么调成土黄色| 一级片'在线观看视频| 亚洲综合色惰| av播播在线观看一区| 男女啪啪激烈高潮av片| av天堂中文字幕网| 亚洲性久久影院| 久久久精品免费免费高清| 妹子高潮喷水视频| 日韩亚洲欧美综合| 欧美zozozo另类| 在线观看美女被高潮喷水网站| 国产精品欧美亚洲77777| 亚洲综合色惰| 熟女电影av网| 不卡视频在线观看欧美| 久久这里有精品视频免费| 一级爰片在线观看| 老熟女久久久| 妹子高潮喷水视频| 91久久精品国产一区二区三区| 国产精品无大码| 国产爱豆传媒在线观看| 国产成人freesex在线| 麻豆成人av视频| 国产精品国产av在线观看| 青春草国产在线视频| 国产成人a∨麻豆精品| 中国国产av一级| 亚洲伊人久久精品综合| 亚洲精品aⅴ在线观看| 亚洲国产日韩一区二区| 极品教师在线视频| 亚洲欧美中文字幕日韩二区| 日本一二三区视频观看| av国产久精品久网站免费入址| 久久97久久精品| 国产精品一及| 亚洲欧美精品专区久久| 你懂的网址亚洲精品在线观看| 色5月婷婷丁香| 蜜桃在线观看..| 中文资源天堂在线| 老司机影院成人| a级毛片免费高清观看在线播放| 久久久成人免费电影| 黄片wwwwww| 国产精品偷伦视频观看了| 久久久亚洲精品成人影院| 99久久精品一区二区三区| av线在线观看网站| 天堂俺去俺来也www色官网| a级毛色黄片| 久久人人爽人人片av| 欧美xxxx性猛交bbbb| 91精品伊人久久大香线蕉| 成人18禁高潮啪啪吃奶动态图 | 日韩av在线免费看完整版不卡| h日本视频在线播放| 久久精品国产鲁丝片午夜精品| 热re99久久精品国产66热6| 亚洲美女搞黄在线观看| 欧美一级a爱片免费观看看| 免费少妇av软件| 亚洲人成网站高清观看| 麻豆成人午夜福利视频| 97热精品久久久久久| 在线观看免费高清a一片| 91精品伊人久久大香线蕉| 亚洲高清免费不卡视频| 女人十人毛片免费观看3o分钟| 久久精品久久久久久噜噜老黄| 亚洲av中文av极速乱| 边亲边吃奶的免费视频| 22中文网久久字幕| freevideosex欧美| 亚洲av二区三区四区| 在线播放无遮挡| 免费看光身美女| 女性生殖器流出的白浆| 国产高清不卡午夜福利| 免费av不卡在线播放| 在线观看免费视频网站a站| tube8黄色片| 高清毛片免费看| 中文乱码字字幕精品一区二区三区| 精品久久久精品久久久| 一区二区三区免费毛片| av在线老鸭窝| 亚洲人成网站在线播| 啦啦啦视频在线资源免费观看| 激情五月婷婷亚洲| 日韩强制内射视频| 最近中文字幕高清免费大全6| 亚洲,一卡二卡三卡| 国产亚洲91精品色在线| 久久久久精品久久久久真实原创| 免费观看a级毛片全部| 亚洲成人中文字幕在线播放| 精品久久久久久电影网| 18禁裸乳无遮挡动漫免费视频| 精华霜和精华液先用哪个| av播播在线观看一区| 一区二区三区精品91| 欧美日韩视频精品一区| 又大又黄又爽视频免费| 最新中文字幕久久久久| 成人免费观看视频高清| 又黄又爽又刺激的免费视频.| 久久久久久久久久久丰满| 欧美97在线视频| 国产伦精品一区二区三区四那| 欧美一级a爱片免费观看看| 在线免费十八禁| 免费观看无遮挡的男女| 免费大片18禁| 视频中文字幕在线观看| 嫩草影院新地址| 中文精品一卡2卡3卡4更新| 在线观看免费日韩欧美大片 | 国产黄片美女视频| 制服丝袜香蕉在线| 久久99蜜桃精品久久| 十分钟在线观看高清视频www | 国产欧美日韩一区二区三区在线 | 国产v大片淫在线免费观看| 男人狂女人下面高潮的视频| 五月天丁香电影| 乱系列少妇在线播放| 国内精品宾馆在线| 少妇 在线观看| 又爽又黄a免费视频| 99热这里只有是精品50| 国产成人a区在线观看| av在线老鸭窝| 男女无遮挡免费网站观看| 极品少妇高潮喷水抽搐| 国产一区有黄有色的免费视频| 少妇精品久久久久久久| 国产精品麻豆人妻色哟哟久久| 午夜免费男女啪啪视频观看| 在线精品无人区一区二区三 | 岛国毛片在线播放| 成年女人在线观看亚洲视频| 亚洲人成网站高清观看| 国产日韩欧美亚洲二区| 亚洲精品一区蜜桃| 嘟嘟电影网在线观看| 人妻系列 视频| 国产乱人偷精品视频| 久久精品久久久久久久性| 又爽又黄a免费视频| 免费人成在线观看视频色| 伊人久久国产一区二区| 精品亚洲乱码少妇综合久久| 国产大屁股一区二区在线视频| 亚洲美女搞黄在线观看| 美女视频免费永久观看网站| 国产91av在线免费观看| 精品99又大又爽又粗少妇毛片| 亚洲av不卡在线观看| 久久久色成人| 成年免费大片在线观看| 国产 一区 欧美 日韩| 亚洲人成网站在线观看播放| 国产毛片在线视频| 国产熟女欧美一区二区| 这个男人来自地球电影免费观看 | 人人妻人人添人人爽欧美一区卜 | 舔av片在线| 综合色丁香网| 日本欧美国产在线视频| a级毛色黄片| 免费黄频网站在线观看国产| 国产色婷婷99| 亚洲人成网站在线观看播放| 亚洲色图av天堂| 色网站视频免费| 欧美精品人与动牲交sv欧美| 亚洲欧美日韩无卡精品| 亚洲国产欧美在线一区| 国产欧美另类精品又又久久亚洲欧美| 高清av免费在线| 欧美xxxx黑人xx丫x性爽| 国产黄片视频在线免费观看| 精品午夜福利在线看| 美女国产视频在线观看| 免费不卡的大黄色大毛片视频在线观看| 久久ye,这里只有精品| 高清日韩中文字幕在线| 舔av片在线| 少妇猛男粗大的猛烈进出视频| 午夜福利影视在线免费观看| 国产一级毛片在线| 人人妻人人爽人人添夜夜欢视频 | 国产无遮挡羞羞视频在线观看| 亚洲三级黄色毛片| 久久久久久久久久人人人人人人| 亚洲一区二区三区欧美精品| 久久久国产一区二区| 国语对白做爰xxxⅹ性视频网站| 精品99又大又爽又粗少妇毛片| 91精品一卡2卡3卡4卡| 男女边摸边吃奶| 少妇的逼好多水| 老司机影院成人| 18禁裸乳无遮挡动漫免费视频| 欧美日韩视频高清一区二区三区二| 熟女人妻精品中文字幕| 一级毛片aaaaaa免费看小| 你懂的网址亚洲精品在线观看| 国产精品av视频在线免费观看| 日韩成人av中文字幕在线观看| 国产精品国产三级国产专区5o| a级一级毛片免费在线观看| 人体艺术视频欧美日本| 一区在线观看完整版| 丰满乱子伦码专区| 春色校园在线视频观看| 蜜桃亚洲精品一区二区三区| 国产 一区 欧美 日韩| 久久久欧美国产精品| 老司机影院毛片| 最黄视频免费看| 青青草视频在线视频观看| 中文乱码字字幕精品一区二区三区| 亚洲不卡免费看| 精品视频人人做人人爽| 精品国产露脸久久av麻豆| 精品久久久久久久久亚洲| 97在线视频观看| 黄色怎么调成土黄色| 成年人午夜在线观看视频| 国产av精品麻豆| 国产精品无大码| 亚洲av电影在线观看一区二区三区| 高清欧美精品videossex| av线在线观看网站| 一个人看的www免费观看视频| 亚洲成人中文字幕在线播放| 久久人人爽av亚洲精品天堂 | 成人亚洲精品一区在线观看 | 在线观看美女被高潮喷水网站| 在线免费观看不下载黄p国产| 国产一区有黄有色的免费视频| 人体艺术视频欧美日本| 少妇 在线观看| 日日啪夜夜爽| 国产一区二区三区综合在线观看 | 亚洲国产高清在线一区二区三| 日韩伦理黄色片| 国产黄片视频在线免费观看| 国产一级毛片在线| 99国产精品免费福利视频| 国产极品天堂在线| 久久ye,这里只有精品| 熟女人妻精品中文字幕| 成人午夜精彩视频在线观看| 交换朋友夫妻互换小说| 色哟哟·www| 一级黄片播放器| 国产又色又爽无遮挡免| 婷婷色麻豆天堂久久| 成人黄色视频免费在线看| 亚洲av电影在线观看一区二区三区| 婷婷色综合www| 亚洲成人一二三区av| 99热这里只有是精品50| 成年美女黄网站色视频大全免费 | 亚洲av综合色区一区| 亚洲久久久国产精品| 欧美bdsm另类| 黑人猛操日本美女一级片| 婷婷色综合www| 国产色婷婷99| 99re6热这里在线精品视频| 国产精品福利在线免费观看| 久久韩国三级中文字幕| 又黄又爽又刺激的免费视频.| 身体一侧抽搐| 一级二级三级毛片免费看| 亚洲欧美日韩另类电影网站 | 日日啪夜夜爽| 国产在线男女| kizo精华| 永久网站在线| 99久久中文字幕三级久久日本| 免费久久久久久久精品成人欧美视频 | 午夜激情福利司机影院| 高清在线视频一区二区三区| 欧美最新免费一区二区三区| 久久久久视频综合| 欧美xxxx黑人xx丫x性爽| 亚洲精品456在线播放app| 亚洲在久久综合| 国产成人91sexporn| 嫩草影院新地址| 中文字幕免费在线视频6| 高清午夜精品一区二区三区| 亚洲自偷自拍三级| 欧美日韩国产mv在线观看视频 | 日韩中字成人| 久久综合国产亚洲精品| 韩国高清视频一区二区三区| 国产在线免费精品| av免费在线看不卡| 国精品久久久久久国模美| 亚洲精品国产av蜜桃| 国产免费福利视频在线观看| 一本—道久久a久久精品蜜桃钙片| 性色avwww在线观看| 成人综合一区亚洲| 91久久精品电影网| 97超视频在线观看视频| 国产精品久久久久成人av| 亚洲欧美成人精品一区二区| 亚洲激情五月婷婷啪啪| 一级片'在线观看视频| 秋霞在线观看毛片| 国产成人freesex在线| 成人毛片a级毛片在线播放| 一级毛片电影观看| 99热这里只有是精品在线观看| av在线app专区| 免费观看在线日韩| 国产精品一区二区在线不卡| 干丝袜人妻中文字幕| 大又大粗又爽又黄少妇毛片口| 日韩av在线免费看完整版不卡| 最近中文字幕高清免费大全6| 五月天丁香电影| 国产av一区二区精品久久 | 99热国产这里只有精品6| 国产成人免费无遮挡视频| 精品一区二区三卡| 男人舔奶头视频| 午夜激情久久久久久久| 黄色视频在线播放观看不卡| 人妻一区二区av| 日本免费在线观看一区| 亚州av有码| 免费看av在线观看网站| 欧美日本视频| 亚洲精品自拍成人| 好男人视频免费观看在线| 少妇的逼水好多| 我的老师免费观看完整版| 夫妻午夜视频| 国产在线一区二区三区精| 身体一侧抽搐| 五月玫瑰六月丁香| 麻豆成人午夜福利视频| 18+在线观看网站| 日韩一本色道免费dvd| av免费在线看不卡| 久久精品国产亚洲网站| 日本欧美视频一区| 男女免费视频国产| 成人影院久久| 中文字幕亚洲精品专区| 天堂中文最新版在线下载| 丰满迷人的少妇在线观看| 国产69精品久久久久777片| h日本视频在线播放| 国产成人免费观看mmmm| 国产精品久久久久久av不卡| av黄色大香蕉| av.在线天堂| 中国国产av一级| 精品99又大又爽又粗少妇毛片| 99视频精品全部免费 在线| 国产男女内射视频| 嘟嘟电影网在线观看| 欧美区成人在线视频| 日韩强制内射视频| 纯流量卡能插随身wifi吗| 1000部很黄的大片| 青青草视频在线视频观看| 九九爱精品视频在线观看| 我的老师免费观看完整版| 国产高清国产精品国产三级 | 欧美人与善性xxx| 久久久久久久久大av| 亚洲人与动物交配视频| 十八禁网站网址无遮挡 | 国产精品久久久久久久电影| 在线观看国产h片| 中文在线观看免费www的网站| 欧美日韩视频精品一区| 午夜福利在线观看免费完整高清在| a级一级毛片免费在线观看| 高清日韩中文字幕在线| 一级黄片播放器| 亚洲在久久综合| 亚洲成人中文字幕在线播放| 毛片女人毛片| 日韩亚洲欧美综合| 欧美一级a爱片免费观看看| 亚洲精品自拍成人| 毛片一级片免费看久久久久| 精品久久久噜噜| 亚洲av中文av极速乱| 亚洲国产最新在线播放| 国产精品人妻久久久久久| 菩萨蛮人人尽说江南好唐韦庄| 99热这里只有是精品在线观看| 男女免费视频国产| 亚洲aⅴ乱码一区二区在线播放| 人妻 亚洲 视频| www.色视频.com| 如何舔出高潮| 日韩av不卡免费在线播放| 国产黄片视频在线免费观看| 97超碰精品成人国产| 精品久久国产蜜桃| videos熟女内射| 国内揄拍国产精品人妻在线| 亚洲精品第二区| 亚洲国产高清在线一区二区三| 国产av国产精品国产| 国产精品国产av在线观看| 久久久国产一区二区| 高清在线视频一区二区三区| 少妇被粗大猛烈的视频| 久久精品国产亚洲网站| 国产精品一区二区在线观看99| 丰满乱子伦码专区| 精品亚洲成a人片在线观看 | 91精品伊人久久大香线蕉| 91精品一卡2卡3卡4卡| 一级片'在线观看视频| 亚洲精品乱码久久久v下载方式| 日产精品乱码卡一卡2卡三| 麻豆成人午夜福利视频| 久久国产精品大桥未久av | 久久ye,这里只有精品| 大又大粗又爽又黄少妇毛片口| 日本vs欧美在线观看视频 | 国产精品久久久久久精品古装| 22中文网久久字幕| 丰满少妇做爰视频| 伦精品一区二区三区| 久久久久视频综合| 精品99又大又爽又粗少妇毛片| 亚洲激情五月婷婷啪啪| 美女cb高潮喷水在线观看| 女性生殖器流出的白浆| 青春草亚洲视频在线观看| 免费少妇av软件| 亚洲国产精品专区欧美| 亚洲精品456在线播放app| 亚洲精品亚洲一区二区| 99久久精品国产国产毛片| 成人亚洲精品一区在线观看 | 秋霞在线观看毛片| 亚洲av日韩在线播放| 日日摸夜夜添夜夜爱| 久久久欧美国产精品| 国产亚洲91精品色在线| 内地一区二区视频在线| 亚洲精品视频女| 日韩亚洲欧美综合| 三级经典国产精品| 日韩强制内射视频| 亚洲人成网站高清观看| 亚洲美女视频黄频| 少妇熟女欧美另类| 成人毛片60女人毛片免费| 精品久久久久久久末码| 国产一区二区三区综合在线观看 | 大片电影免费在线观看免费| 一个人免费看片子| 最近最新中文字幕免费大全7| 久久久久久人妻| 久久99蜜桃精品久久| 99热这里只有是精品在线观看| 国产精品国产三级专区第一集| a 毛片基地| 亚洲国产高清在线一区二区三| 精品久久久久久久末码| 五月开心婷婷网| 日本爱情动作片www.在线观看| 日日摸夜夜添夜夜添av毛片| 国产亚洲最大av| 99热这里只有是精品50| 久久久久久久久大av| 高清毛片免费看| 九九久久精品国产亚洲av麻豆| 亚洲国产成人一精品久久久| 少妇人妻 视频| 亚洲精品国产av蜜桃| 观看av在线不卡| 免费看av在线观看网站| 国产精品成人在线| 国产精品一二三区在线看| 午夜日本视频在线| 久久久久性生活片| 亚洲欧美一区二区三区黑人 | 久久女婷五月综合色啪小说| 久久6这里有精品| 人妻制服诱惑在线中文字幕| 亚洲欧美日韩卡通动漫| 国产一区二区三区综合在线观看 | 国产视频首页在线观看| 80岁老熟妇乱子伦牲交| 男女边吃奶边做爰视频| a级一级毛片免费在线观看| 少妇的逼好多水| 亚洲丝袜综合中文字幕| 高清av免费在线| 成年女人在线观看亚洲视频| 日本猛色少妇xxxxx猛交久久| 欧美97在线视频| 国产亚洲欧美精品永久| 国产淫语在线视频| 国产精品三级大全| 日韩欧美一区视频在线观看 | 色婷婷av一区二区三区视频| 精品亚洲乱码少妇综合久久| av免费在线看不卡| 国产一区有黄有色的免费视频| 五月玫瑰六月丁香| 99热这里只有是精品50| 日本猛色少妇xxxxx猛交久久| 成年女人在线观看亚洲视频| 97在线视频观看| 国产日韩欧美在线精品| 免费av中文字幕在线| 国产淫语在线视频| 97在线视频观看| 久久亚洲国产成人精品v| 国产精品女同一区二区软件| 亚洲欧美日韩东京热| 免费观看无遮挡的男女| 国产精品一区二区在线观看99| 日韩欧美精品免费久久| 天堂中文最新版在线下载| 久久久久精品性色| 亚洲激情五月婷婷啪啪| 妹子高潮喷水视频| 少妇被粗大猛烈的视频|