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

    基于Boussinesq型方程的錢塘江涌潮數(shù)值模擬

    2022-01-25 07:03:04黃婷張懷石耀霖
    地球物理學報 2022年1期
    關鍵詞:鹽官潮頭錢塘江

    黃婷,張懷,2*,石耀霖

    1 中國科學院大學地球與行星科學學院,計算地球動力學重點實驗室,北京 100049 2 南方海洋科學與工程廣東省實驗室 (珠海),廣東珠海 519080

    0 引言

    涌潮是一種發(fā)生在漲潮期間的潮汐現(xiàn)象,來潮沿著河流或狹窄海灣逆流而上,形成波長和頻率不斷變換的潮波.涌潮的形成基于來潮、河流條件以及河道地形之間的微妙平衡,并非所有的河口都會發(fā)生涌潮.據(jù)估計,全世界發(fā)生涌潮的河口有四百多個,分布在除南極洲以外的所有大陸,河口附近的環(huán)境、生態(tài)以及文化都將受到涌潮所帶來的影響(Chanson,2009).杭州灣是著名的涌潮地,而自澉浦以上的錢塘江出口段(圖1)是世界聞名的強涌潮河口(譚維炎等,1995),因錢塘江河床泥沙易沖易淤,為治理錢塘江河口環(huán)境,20 世紀60 年代以來,治江縮窄工程不斷推進(尤愛菊等,2010),河流邊界內(nèi)移,使得錢塘江河口處形成形態(tài)復雜多樣的涌潮.最具代表性的潮景是大缺口的“交叉潮”、鹽官的“一線潮”以及老鹽倉的“回頭潮”,每年都會吸引無數(shù)中外游客慕名而來,是我國獨特的自然旅游資源.對錢塘江涌潮的研究,不僅具有學術意義,也具有重要的實際應用價值.

    圖1 錢塘江河口地形示意圖黑色輪廓線是20世紀初錢塘江河流邊界,藍色線是江底地形等高線.Fig.1 Schematic diagram of the Qiantang EstuaryThe black outline indicates the boundary of the Qiantang River in the early 20th century.The blue line indicates the contour line of the topography of the river bottom.

    幾個世紀以來,國內(nèi)外的學者一直致力研究涌潮現(xiàn)象.Saint-Venant(1871)對一維渠道內(nèi)的水流進行觀測,提出了著名的圣維南方程,并將其應用于法國塞納河的涌潮,通過計算來潮的Froude數(shù),判斷涌潮的發(fā)生.一維模型基于理想假設,適用于流速方向易于評估,并且在整個涌潮持續(xù)時間內(nèi)固定的河道(Clamond and Dutykh,2013).趙雪華(1985)應用一維圣維南方程,有機地結合特征線法和特征差分格式,計算了錢塘江自杭州閘口至下游澉浦之間河段涌潮時的水位和流量情況.于普兵和潘存鴻(2010)應用迎風格式有限體積法求解一維圣維南方程,建立了錢塘江涌潮一維數(shù)值模型,并模擬計算從上游富春江電站至下游澉浦河段的潮位情況.

    錢塘江涌潮存在形狀良好的波狀涌潮和潮頭破碎的強涌潮等多種形態(tài),且波狀涌潮的頻散作用較大,需要在控制方程中添加頻散項(潘存鴻等,2008),以往的模擬大多針對的是涌潮的傳播過程并不考慮頻散效應.綜上所述,本文在前人研究的基礎上建立了新的二維高階Boussinesq型方程,采用有限體積法,考慮移動的干濕邊界,對錢塘江涌潮進行數(shù)值模擬研究,探討了在考慮頻散作用下錢塘江涌潮的形成與傳播特征.

    1 模型與數(shù)值方法

    1.1 控制方程

    1.1.1 Boussinesq型方程

    Madsen和S?rensen(1992)從經(jīng)典Boussinesq 方程出發(fā),提出了緩坡地形上的一種Boussinesq型方程,并通過理論分析和應用驗證該方程比經(jīng)典Boussinesq方程適用的水深范圍更廣.經(jīng)典Boussinesq方程適用于水深與波長比值h0/λ0<1/5的水深區(qū)域,而Madsen和S?rensen的方程在調(diào)節(jié)適當?shù)念l散系數(shù)后可將應用水深擴展到h0/λ0<1/2的區(qū)域.Beji和Nadaoka(1996)同樣從經(jīng)典Boussinesq方程出發(fā),通過對非線性項和頻散項的數(shù)學關系代換,修正了Madsen和S?rensen的方程,提出了一種考慮地形變化的二維高階Boussinesq型方程.我們在Beji和Nadaoka工作的基礎上,導出在考慮河道彎曲變化和強烈地形起伏情形下不可壓縮流體的二維高階Boussinesq型方程,表示為

    (1)

    (2)

    (3)

    其中n為曼寧系數(shù),F(xiàn)=F(fx,fy)為科氏力,有

    fx=fv,fy=-fu,

    (4)

    其中f=2Ωsinφ為Coriolis參數(shù),Ω為地球自轉角速度,φ為緯度值.A=A(Ax,Ay)表示水平動量擴散

    (5)

    其中Am為水平渦流黏度系數(shù),可以簡單地采用渦流黏性去模擬湍流混合和波破碎時的能量耗散(Chen et al.,2000;Kennedy et al.,2000).B=(Bx,By)稱為Boussinesq項或是頻散項,是關于時間和空間的混合高階導數(shù),是在靜水壓力平衡基礎上的一個附加項

    (6)

    其中h為靜水深,參數(shù)β調(diào)節(jié)著方程的頻散程度.

    1.1.2 方程的線性頻散性

    長波假設意味著,波的特征振幅η0與特征水深h0相比較小,即系數(shù)ε=η0/h0?1,波的傳播體現(xiàn)為波的相位傳播,且波的相速度與stokes一階理論(即Airy線性波理論)表述的一致,于是可以在長波假設下對控制方程進行線性頻散性分析.

    線性頻散關系表示為:

    (7)

    式中ω為波的頻率,k=2π/λ為波數(shù),h為靜水深.進而得到波速c與水深的線性頻散關系

    (8)

    將式(8)在kh=0處進行泰勒展開有

    (9)

    截取到O(k4h4)項,即方程的精度達到O(σ4)階.

    如前所述的控制方程(1)的頻散關系表示為(Beji and Nadaoka,1996)

    (10)

    于是系統(tǒng)的相速度可以表示為

    (11)

    圖2 非靜水壓力梯度的影響虛線表示初始波的形態(tài),實線表示當前波的形態(tài),空心箭頭表示額外水平壓力梯度方向.Fig.2 The influence of non-hydrostatic pressure gradientThe dashed line represents the shape of the initial wave,the solid line represents the shape of the current wave,and the hollow arrow represents the direction of the different horizontal pressure gradients.

    1.2 數(shù)值方法

    1.2.1 網(wǎng)格劃分

    在二維空間,采用節(jié)點中心格式設置空間交錯的均勻網(wǎng)格,變量η、u和v分別在三套網(wǎng)格中進行求解,如圖3所示.

    圖3 交錯網(wǎng)格節(jié)點與控制體積示意圖三個物理量η,u和v所在計算節(jié)點分別表示為圓形、三角形和正方形,所在控制體積分別填充正交線、左斜杠和右斜杠表示.Fig.3 Schematic diagram of node and control volume of the computational gridThe nodes of three independent variables η,u,and v are expressed as circular,triangular,and square respectively.An orthogonal line,left slash,and right slash represents the corresponding control volume.

    1.2.2 連續(xù)性方程

    對連續(xù)性方程(1),在時間上采用向前歐拉格式,在空間上使用有限體積法.于ηi,j節(jié)點處,在以節(jié)點為中心的控制體積ΔVi,j內(nèi)對連續(xù)性方程進行積分,引用高斯定理得

    +Δx[(Dv)n-(Dv)s]=0,

    (12)

    式中S為控制體積ΔVi,j的外邊界,n為邊界S外法線的單位向量,下標e、w、n、s分別指示邊界S的東、西、北、南四條邊界.Δx和Δy分別為在x和y方向上的空間步長,控制體積ΔVi,j=Δx×Δy.邊界處的水深使用相鄰節(jié)點的水深平均值,得到

    (13)

    式中Δt為時間步長,k為當前時間步.

    1.2.3 動量守恒方程

    形如上述動量守恒方程的保守部分可以表示為關于變量φ的形式

    (14)

    式中ν為擴散系數(shù).于φi,j節(jié)點處,在以節(jié)點為中心的控制體積ΔVi,j內(nèi)對方程進行積分得到

    (15)

    再在時間上采用向前歐拉格式得到:

    (16)

    式中m為S的各個邊界組分(e、w、n、s),Δm為邊界長度,φm為邊界處的變量值,整理后得

    (17)

    l表示邊界處的切線方向單位向量.

    對流項中C表示Courant數(shù),展開有

    上標“+”表示邊界處速度方向與邊界法線方向相同,上標“-”表示邊界處速度方向與邊界法線方向相反,即可以概括為

    (19)

    對于擴散項展開則有

    (20)

    對于控制體積邊界處的值可以由其上下游處的節(jié)點值插值得到,如東邊界處的值φe可以表示為如下形式:

    (21)

    式中函數(shù)Ψ稱為通量限制器,與迎風側的梯度與逆風側的梯度比值r有關,r表示為

    (22)

    式中φu表示當前邊界處的上游節(jié)點值,φu u表示φu的上游節(jié)點值,φd表示當前邊界處的下游節(jié)點值,φd d表示φd的下游節(jié)點值,如圖4所示.

    圖4 插值示意圖黑色實心圓和空心圓分別表示節(jié)點位置和相應節(jié)點處物理量的值,黑色實心三角和空心三角分別表示邊界位置和相對應的物理量的值.Fig.4 Schematic diagram of interpolationThe black solid circle and the hollow circle represent the position of the node and the value of the corresponding physical variables.The black solid triangle and the hollow triangle represent the position of the boundary and the value of the corresponding physical variables,respectively.

    根據(jù)不同的使用需要,通量限制器函數(shù)Ψ具有不同的形式,最常被使用的格式見表1(Versteeg and Malalasekera,2007).

    表1 常見的通量限制器函數(shù)表達形式(Versteeg and Malalasekera,2007)Table 1 Most popular limiter functions found in the literature (Versteeg and Malalasekera,2007)

    Sweby(1984)根據(jù)r-Ψ關系給出了格式符合TVD方案的充要條件,還提出了對達到二階以上高階精度格式的要求.本文采用的QUICK格式可以保證TVD特性,且至少具有二階精度.在均勻正交網(wǎng)格下,對于擴散項,使用TVD-QUICK格式有

    (23)

    于是得到方程保守部分的最終形式

    (24)

    對于控制方程中的Boussinesq項在空間上采用有限差分法,我們將在算法測試中對一維方程的Boussinesq項進行具體離散.

    1.2.4 穩(wěn)定性分析

    (25)

    1.3 涌潮的生長

    初始速度場分布與初始波形為

    (26)

    式中u0和a均為常數(shù).

    為簡便操作,取初始靜水深h=1,得到控制方程

    (27)

    η,u,x,t均為無量綱變換之后的變量,其中頻散項有

    (28)

    對方程(27)使用數(shù)值格式中介紹的方法進行數(shù)值離散,將Boussinesq項離散之后得到

    (29)

    我們分別計算頻散系數(shù)β=0和β=1/5的情況,并參照不考慮頻散項的靜水壓力模型,測試中選擇u0=0.1,dx=0.6,dt=0.25,a=6(Peregrine,1966),計算結果見圖5.

    圖5 波狀涌潮的生長過程.從下往上的曲線簇分別表示時間從0開始,以5增長直至20的計算結果灰色實線、紅色實線分別表示使用Peregrine 1966年淺水模型和Boussinesq模型的數(shù)值計算結果,黑色虛線、黑色實線分別表示使用文中所述方程及算法在β=0和β=1/5時所得的計算結果.Fig.5 Growth of an undular bore.The curves from bottom to top represent the calculation results of time from 0 to 20 with 5 as the interval The gray solid lines and the red solid lines represent the computation results by Peregrine 1966 using the shallow water model and the Boussinesq model respectively,while the black dashed lines and the black solid lines represent the computation results using the equation and algorithm described above when β=0 and β=1/5 respectively.

    從圖5中可以看到,在模擬時間內(nèi),僅考慮靜水壓力的數(shù)值模型計算結果(灰色實線),波前緣逐漸變陡,但波后方仍舊保持平滑,波前緣漲幅保持在0.1;使用Boussinesq型方程的計算結果在初始一段時間內(nèi)和靜水壓力理論的結果重合,而后波前緣逐漸隆起形成波峰,在波后方出現(xiàn)波谷,再現(xiàn)了波的頻散(二次自由面起伏)特征,形成了波狀涌潮.當頻散系數(shù)β=0時,使用如前所述的算法計算結果(黑色虛線)與Peregrine 1966年的測試結果(紅色實線)相當吻合,波前緣漲幅在0.15左右.而當頻散系數(shù)β=1/5時考慮了更高階的垂直壓力分布改變對水平壓力梯度變化的影響,因此計算得出的漲幅也會更大(黑色實線),波前緣漲幅在0.17左右.由測試結果可見,在涌潮模擬中對靜水壓力的校正是非常重要的,方程需要考慮頻散項的影響.

    1.4 模擬區(qū)域與初邊值設置

    海面的潮汐變化可以分解為許多簡單的、規(guī)則的簡諧振動,每個振動稱為一個分潮,大多數(shù)地區(qū)的潮汐主要成分是由月球引起的主太陰半日潮.錢塘江河口潮汐現(xiàn)象即為一種半日潮.2000年,浙江省錢塘江管理局與浙江省河口海岸研究所對錢塘江河口進行了一次大規(guī)??疾煨杂砍庇^測,錢塘江河口在地形上有兩個顯著的特點(林炳堯等,2002),一是杭州灣河口呈喇叭形;二是自乍浦往上至倉前、七堡之間,存在宏大的沙坎,其頂部比乍浦河床高出10 m左右.錢塘江潮波的變化、涌潮的形成和發(fā)展與這兩種地形特點密切相關.澉浦是我國潮差最大的海域之一,但是當潮波通過杭州灣時,變形的程度并不明顯,澉浦下游的潮位變化過程仍然接近簡諧曲線.一般情況下涌潮自大缺口至鹽官之間形成、發(fā)展,最終在七堡附近湮滅,大潮時期涌潮最遠可到達杭州以上的聞家堰至富陽河段(Pan and Lu,2011).

    根據(jù)實際涌潮的特點,我們的模型區(qū)域上游邊界設置在閘口處,下游邊界設置在乍浦下游,計算區(qū)域范圍取值為120.15°E—121.25°E/30.0°N—30.7°N(圖1),地形數(shù)據(jù)水平分辨率精度約90 m.河流東邊界給定水面高度變化,使用與實際記錄相近的正弦波形函數(shù)來模擬潮汐引起的強迫場(林炳堯等,2002),根據(jù)杭州灣的潮差分布(Fan et al.,2015)以及實測資料中的錢塘江涌潮的特性(林炳堯,2008),選擇振幅amp=2 m,周期period=12 h.河流西邊界在模擬的時間內(nèi)假定為水深恒定速度自由,即水位邊界條件為給定的恒定水深,速度邊界設置為自由傳輸邊界.模型考慮水波的上岸淹沒情況,因此網(wǎng)格的干濕情況是時刻變化的,形成移動的干濕邊界,網(wǎng)格干濕判斷條件如下,取最小水深Dmin=0.01 m,即有

    (30)

    經(jīng)各式模型的測試與實測校對表明,錢塘江河口底部的曼寧系數(shù)較小,約為0.01(潘存鴻等,2008),或者更小,取曼寧系數(shù)n=0.01,頻散系數(shù)β=1/5,初始時間步長Δt=0.5 s.

    2 模擬結果與分析

    我們模擬了潮波進入錢塘江河口,上行直至上游閘口處,形成涌潮的整個發(fā)展過程.結果主要展示涌潮經(jīng)過不同河段時的波高場以及速度場的變化,旨在從涌潮高度、涌潮傳播速度、潮景等涌潮傳播表征進行分析,探究在不同區(qū)域河段涌潮的特征.將陸地區(qū)域的名稱,作為其所對應相鄰的錢塘江河段的暫代名稱,用于對涌潮經(jīng)過此河段區(qū)域的定性分析.對錢塘江涌潮的觀測資料顯示,涌潮的形成基本在尖山一帶,自大缺口至鹽官河段發(fā)展達到最大,再往上行,涌潮高度逐漸減小,直至消失,其上界一般在七堡附近,最遠也可到達杭州以上(譚維炎等,1995;林炳堯等,2002;潘存鴻等,2008).

    因關注的重點在涌潮本身,我們將濕單元的數(shù)據(jù)進行提取并繪制成波高場圖和速度場圖.在我們的模型計算結果中,涌潮進入形成、發(fā)展和消亡三個階段的位置和實際記錄結果吻合.涌潮形成于尖山河段,該河段下游的水深在40 m左右,相對較深,此時的涌潮高度較小,潮頭并不明顯,潮波因江底沙洲和沿岸地形的影響分成兩股波,在尖山和大缺口河段間形成“交叉潮”(圖6a).隨后涌潮從大缺口進入了較為平直狹窄的鹽官甬道,大量水體涌入,表現(xiàn)出“一線潮”(圖6b),且強度隨著深入不斷增大(圖6c).在老鹽倉拐角處,因轉折角度較大,涌潮被反射形成“回頭潮”(圖6d).鹽官上游之后的河段水深較淺,基本處于10 m以下,江底摩擦產(chǎn)生的能量耗散作用明顯,再加上彎曲地形與涌潮之間的碰撞,涌潮強度開始減弱.涌潮后續(xù)經(jīng)過幾個拐角上溯至七堡河段(圖6e),之后進入閘口河段走向消亡階段(圖6f).2.2小節(jié)將會對三處典型潮景的涌潮特征進行具體分析.

    2.1 涌潮特征的沿程變化

    涌潮高度、涌潮速度以及Froude數(shù)是表征涌潮特性的三個重要指標.涌潮高度ΔH一般指涌潮引起的平均水面上漲高度,ΔH=hd-hu,hu、hd分別表示涌潮潮頭前方(即波前緣的上游方向)水深和潮頭后方(即波前緣的下游方向)水深,如圖7所示,相對漲幅表示為ΔH/hu.

    使用涌潮潮頭傳播的速度描述涌潮速度,涌潮速度c的流量守恒方程和動量守恒方程可以表示為(以指向上游方向為流速的正方向)

    c(hd-hu)=udhd-uuhu,

    (31)

    (32)

    其中uu、ud分別表示潮前流速和潮后流速(圖7).

    圖7 涌潮傳播示意圖uu、ud分別表示潮前流速和潮后流速,hu、hd表示潮前水深和潮后水深,c表示涌潮流速.Fig.7 Definition sketch of a bore propagationuu and ud are respectively the flow velocity before and after the tidal bore,hu and hd indicate the flow depths before and after the tidal bore respectively,c is the bore celerity.

    由式(31)、(32)可以分別得到

    (33)

    (34)

    聯(lián)合整理式(33)、(34)之后得到由潮前、潮后流速表示的涌潮傳播速度c為

    (35)

    (36)

    由式(36)可以看到,F(xiàn)roude數(shù)與潮前水深和潮后水深的比值相關,F(xiàn)r越大,涌潮引起的前后水面漲幅落差也就越大.而c同時可以由Fr表示為

    (37)

    可見涌潮速度c與Froude數(shù)Fr、流速以及水深有關.我們選取靠近河道中心位置處的涌潮剖面,計算涌潮上述三個特征以及漲幅的沿程變化情況,具體計算結果見表 2.

    表2 涌潮特征沿程變化Table 2 Variation of bore characteristics along the river

    在不考慮老鹽倉處的“回頭潮”的情況下,我們將表中的數(shù)據(jù)按照錢塘江涌潮上溯的傳播方向繪制了涌潮特征沿程變化,見圖8.涌潮的Fr隨ΔH自涌潮的形成階段、發(fā)展階段直至消亡階段,整體呈現(xiàn)出先增大再減小的規(guī)律.自鹽官以上河段的水深相差不大,在10 m左右,因此c也整體呈現(xiàn)先增大再減小的規(guī)律.涌潮于鹽官河段處上溯,強度一直增大,于外五工段區(qū)域河道處規(guī)模達到最大,ΔH約為3.5 m,F(xiàn)r約為1.54,c約為9 m·s-1,此時段的涌潮潮頭陡立,自由表面可能失穩(wěn),波開始破碎.

    圖8 涌潮特征沿程變化曲線方塊、三角、圓分別表示涌潮速度c、涌潮高度ΔH和涌潮的Froud數(shù)Fr.Fig.8 Curves of variation of bore characteristics along the riverThe celerity c,bore height ΔH,and Froude number Fr are expressed as square,triangular,and circular.

    涌潮經(jīng)過老鹽倉處的急轉彎,在反射和推擠作用下形成“回頭潮”,“回頭潮”的強度并不大,ΔH在0.7 m左右,F(xiàn)r接近1.繼續(xù)上溯的涌潮強度開始減弱但規(guī)模依舊較大,三工段至倉前的整個河段涌潮變化不大,ΔH約為2.5 m,F(xiàn)r接近1.4,c約為8 m·s-1.在經(jīng)過七堡至閘口河段,涌潮規(guī)模和強度都在減弱,F(xiàn)r減小到近1.2,ΔH也僅有1 m左右,涌潮處在消亡階段.閘口河段處的水深較七堡處要大,于是出現(xiàn)在閘口河段處的涌潮速度要稍微大于七堡河段處的涌潮速度,約為6 m·s-1.

    從實測資料分析中發(fā)現(xiàn),涌潮高度一般于鹽官河段達到極值,最大可達3 m以上,涌潮傳播速度在3~10 m·s-1之間,此后涌潮開始衰減,涌潮高度也逐漸減小,沿程的涌潮速度也呈現(xiàn)從小到大,再從大到小的變化規(guī)律,最大流速出現(xiàn)在鹽官河段(林炳堯,2008;潘存鴻等,2008),可見我們的模擬結果和實際觀測現(xiàn)象所得的特點是吻合的.

    2.2 涌潮的發(fā)展與潮景的形成

    潮景的形成是涌潮不同發(fā)展階段的標志,我們復演并選取了最具代表性的三個河段位置處的潮景,分別是大缺口處的“交叉潮”,鹽官處的“一線潮”以及老鹽倉處的“回頭潮”.

    2.2.1 交叉潮

    尖山河段是錢塘江河口的一個轉折處,河口近似呈“之”字形,從此處潮波由較為深廣的水域經(jīng)大缺口進入淺窄且平直的鹽官河段.大缺口附近河口呈平躺的“U”字形,因泥沙淤積,錢塘江河口下游發(fā)育一巨大江中沙洲(圖1),在大缺口下方尖山河段分汊河勢以及江中沙洲的影響下,此處常常發(fā)育“交叉潮”(潘存鴻等,2008).

    潮波在通過沙洲時被分成了兩股波,一股沿著北岸向西北傳播,另一股沿著南岸向西南方向傳播(圖9a),潮波帶來整個尖山河段水面的上漲,并未形成明顯的涌潮現(xiàn)象.兩股暗潮相互作用,使得西北向潮波得到加強,于尖山河段下游形成較為清晰整齊的波,涌潮形成(圖9b),涌潮的方向指向河流的上游.事實上,涌潮的形成并不是瞬時的,而是一個連續(xù)的過程,因此無法具體準確地判斷涌潮形成的精確位置與時間,當肉眼可分辨涌潮形態(tài)時,我們即認為涌潮形成.在模擬時間的80 min左右,從圖9c的涌潮剖面中可以看到,尖山河段潮頭漲幅可以達到1.5 m,涌潮高度僅約0.3 m,計算出Froude數(shù)在1左右,涌潮形成初期波形平緩,波列整齊.

    圖9 尖山河段涌潮圖(a—b)展示了不同時刻潮波引起的尖山河段的水面變化.紅色虛線表示涌潮的波前緣,紅色實線表示截取的涌潮剖面.圖(c)展示了圖(b)中的涌潮西北-東南方向的剖面形態(tài).Fig.9 Tidal bore in JianshanFigures (a—b)show the variation of water surface at Jianshan at different times respectively.The red dashed line indicates the wavefront of the tidal bore wave,and the red solid line indicates the selected tidal bore profile.Figure (c)shows the bore profile in figure (b)from NW to SE.

    在波狀涌潮上溯的過程中于尖山至大缺口河段形成“交叉潮”(圖10).尖山至大缺口河段是一個近似由西向北的轉彎,在潮波與南岸邊發(fā)生反射推擠作用的影響下,南岸形成了一股近似由南向北傳播的涌潮(圖10a),另一股涌潮沿著東南向西北行進,兩股涌潮的潮頭波高皆約為1.5 m,兩股波前形成一個鈍角的“人”字形.此時涌潮高度較小,波形完好,顯示出良好的波的特征.兩股涌潮繼續(xù)上行,進入大缺口河段入口,兩組波列相互干涉,波前形成了一個近似交叉的“十”字,即所謂的“交叉潮”(圖10b),波前緣的夾角依然是鈍角.在“交叉潮”后方清晰可見由干涉產(chǎn)生交錯排列的明暗紋,波峰處漲幅在1.5 m左右,波谷處漲幅在1~1.5 m之間,波后的水面上漲約為2 m.“交叉潮”向大缺口行進,波前緣夾角逐漸減小(圖10c),向北傳播的涌潮接近北岸,涌潮的形態(tài)也重新呈現(xiàn)“人”字形.

    圖10 “交叉潮”波高場圖(a—d)展示了不同時刻尖山河段的水面變化,紅色虛線表示涌潮的波前緣,涌潮從尖山河段上溯至大缺口河段的過程中形成“交叉潮”.Fig.10 Wave height field of cross-shaped tidal boreFigures (a—d)show the variation of water surface at Jianshan at different times respectively.The red dashed line indicates the wavefront of the tidal bore wave.The cross-shape tidal bore is formed when the bore goes upstream from Jianshan to Daquekou.

    “交叉潮”形成時,涌潮波的傳播速度雖然接近20 m·s-1,但流場的速度較小,約為1 m·s-1,波的傳播速度比流體速度大得多.速度方向的分布和兩股潮的行進方向一致,波峰的位置處流場速度較大,而波谷位置處的流場速度較小,流場也呈現(xiàn)出與波場類似的明暗條紋分布,如圖11所示.

    圖11 “交叉潮”速度場箭頭方向表示速度方向,箭頭長度表示速度大小,紅色虛線表示涌潮的波前緣位置,速度場因速度方向和大小的變化呈現(xiàn)出與波紋類似的明暗條紋.Fig.11 Velocity field of crossed tidal boreThe direction of the arrow indicates the direction of velocity while the magnitude indicates the strength,and the red dashed line indicates thewavefront of the tidal bore wave.Due to the variation of strength and direction,the velocity field presents light and dark fringes.

    2.2.2 一線潮

    鹽官以上的上游河段水深較淺,在10 m以下,因此涌潮在進入淺窄且平直的鹽官河段后,水面高度迅速上漲,潮頭走向幾乎垂直于岸邊,且波列以近乎平行的姿態(tài)上溯,于鹽官河段附近形成較為理想的“一線潮”潮景(圖12).“一線潮”是涌潮的發(fā)展階段,涌潮強度開始增大.鹽官河段處形成強度較大的“一線潮”,涌潮高度約為3 m.之后“一線潮”保持著相似的形態(tài),強度繼續(xù)變大,通過整個鹽官河段到達外五工段河道區(qū)域.

    圖12 “一線潮”波高場圖(a—c)展示了不同時刻潮波引起的鹽官河段的水面變化,涌潮通過鹽官河段時,潮頭處波高等值線變密,潮頭急速抬升,形成“一線潮”,紅線表示截取的涌潮剖面.Fig.12 Wave height field of thread-shape tidal boreFigures (a—c)show the variation of water surface at Yanguan at different times respectively.When the tidal bore reaches the Yanguan outlet,the wave height contour at the bore head becomes denser and the bore head rises steeply,then the thread-shape tidal bore forms.The red line indicates the selected tidal bore profile.

    從波形剖面可以清楚地看到涌潮強度不斷增大的過程.涌潮在進入鹽官河道入口時,依然可見波前緣后方一系列形態(tài)良好的波(圖13a).較為平緩的波前緣成長為陡立的潮頭,且潮頭前的水面會出現(xiàn)略微下降,在鹽官處形成明顯的間斷面(圖13b),計算得到涌潮的Froude數(shù)約為1.45.直至在鹽官河段的出口處,外五工段的附近,涌潮的潮頭一直在上漲,達到將近4 m(圖13c),F(xiàn)roude數(shù)達到1.54.

    圖13 “一線潮”涌潮剖面展示了圖12中選取的涌潮剖面,涌潮自東向西通過鹽官河段,潮頭逐漸陡立,涌潮強度增強,二次面起伏現(xiàn)象明顯,具有強頻散作用.Fig.13 Profile of thread-shape tidal boreThis figure shows the tidal bore profile selected in Fig.12.The tidal bore passes through the Yanguan outlet from east to west,the tidal head becomes steep,the intensity of the tidal bore increases,secondary free-surface undulations are reproduced,and the bore has a strong dispersion effect

    涌潮經(jīng)過鹽官河段的速度場分布近乎與岸邊平行,涌潮的波前緣存在明顯的速度分界.在鹽官河段入口附近,兩股流場交匯處的速度相近,速度大小約為1 m·s-1,波前緣處流體速度接近為零(圖14a).隨著涌潮深入,潮波進入甬道,潮后速度逐漸增大,于鹽官處形成“一線潮”時潮后的下游流場速度約為3 m·s-1(圖14b).“一線潮”上溯速度繼續(xù)增大,鹽官河段的出口處連接老鹽倉拐角,因此在外五工段附近,南岸的速度略微大于北岸(圖14c),靠近南岸的流場速度在3 m·s-1左右.經(jīng)計算的鹽官河段涌潮速度在6~9 m·s-1,與流體速度僅相差幾倍,可見涌潮已不再是簡單的淺水模式,具有很強的頻散效應.

    圖14 “一線潮”速度場圖(a—c)展示了不同時刻鹽官河段的速度分布,箭頭方向表示速度方向,箭頭長度表示速度大小,紅色虛線表示涌潮的波前緣位置.Fig.14 Velocity field of thread-shape boreFigures (a—c)show the velocity distribution at Yanguan at different times respectively.The direction of the arrow indicates the direction of velocity while the magnitude indicates the strength. The red dashed line indicates the wavefront of the tidal bore wave.

    2.2.3 回頭潮

    老鹽倉拐角處夾角接近130°,“一線潮”于拐角處發(fā)生反射,在老鹽倉形成了“回頭潮”的潮景.彎道處造成的碰撞和反射使得涌潮波與岸邊作用加強,且較淺的水深使得江底摩擦作用顯著,涌潮能量被耗散,于是在“回頭潮”出現(xiàn)之后,涌潮強度開始減弱.涌潮發(fā)生反射時,整個流場由下游來潮、反射波以及上游河流三個部分的流場組成(圖15b),外拐角處水面上漲最大可達4 m.沿著河道上行的涌潮疊加部分反射波繼續(xù)上溯,進入三工段區(qū)域河道.三工段河道與鹽官河段相似,狹長且河面寬度較為均勻,涌潮進入河道后,以“一線潮”的形態(tài)前進,其規(guī)模較鹽官河段的“一線潮”要小,涌潮強度開始減弱,涌潮高度減小到約為2 m(圖15c).反向往下游傳播的 “回頭潮”,潮頭漲幅在3 m左右,涌潮高度約為0.7 m(圖15d).涌潮經(jīng)過老鹽倉所帶來的水面上漲在3 m左右,老鹽倉內(nèi)拐角內(nèi)側岸邊處的一些洼地被水填滿.“回頭潮”發(fā)生后,上溯涌潮的強度開始減弱,在經(jīng)過七堡河段的拐角之后,涌潮進入消亡階段.

    圖15 “回頭潮”波高場圖(a—d)展示了不同時刻潮波引起的老鹽倉河段的水面變化,紅色虛線表示涌潮的波前緣,涌潮在老鹽倉處發(fā)生反射,形成自西向東傳播的“回頭潮”.Fig.15 Wave height field of returned tidal boreFigures (a—d)show the variation of water surface at Laoyancang at different times respectively.The red dashed line indicates the wavefront of the tidal bore wave.The tidal bore is reflected at Laoyancang,returned tidal bore forms and spreads from west to east.

    老鹽倉拐角處,上下兩股流場交匯,速度場的分布情況和波高場分布相似(圖16),由來潮、上游河流以及反射波三個速度場組成,涌潮的波前緣存在明顯的速度分界.反射波的速度方向仍然是總體指向上游的,拐角處的速度最小,僅在1 m·s-1左右(圖16b).指向上游的老鹽倉反射波和上溯的涌潮相疊加,以約為3 m·s-1的速度進入三工段區(qū)域河道,與上游的速度流場的交界處呈“一”字形態(tài)(圖16c).部分反射波形成“回頭潮”反向向東傳播,減慢了上溯流場的速度,“回頭潮”的強度也不斷減弱,涌潮同樣具有很強的頻散效應.

    圖16 “回頭潮”速度場圖(a—d)展示了不同時刻老鹽倉河段的水流速度分布,箭頭方向表示速度方向,箭頭長度表示速度大小,紅色虛線表示涌潮的波前緣位置.Fig.16 Velocity field of the returned tidal boreFigures (a—d)show the velocity distribution at Laoyancang at different times respectively.The direction of the arrow indicates the direction of velocity while the magnitude indicates the strength.The red dashed line indicates the wavefront of the tidal bore wave.

    3 結論與展望

    通過對錢塘江河口涌潮的數(shù)值模擬研究,我們得到以下結論:

    (1)錢塘江河口涌潮的形成和發(fā)展以及潮景的產(chǎn)生和錢塘江的地形特征密切相關,模擬結果與實際觀測中得到的結果可以很好吻合.潮波自錢塘江河口經(jīng)過江底沙洲的作用以及與江岸間的相互作用,在尖山至大缺口河段發(fā)育形成“交叉潮”;在鹽官河段,地形較為平坦,兩岸長直且近似平行,涌潮由較為深廣的水域進入淺窄的甬道時,波列逐漸均勻,形成“一線潮”;在老鹽倉拐角處,涌潮發(fā)生反射,形成“回頭潮”;“一線潮”潮景不僅發(fā)生在鹽官河段,三工段附近也出現(xiàn)了波動形態(tài)良好的波狀涌潮,同樣呈“一”字形.

    (2)本次模擬復演了最具代表性的“交叉潮”、“一線潮”以及“回頭潮”三大潮景,鹽官河段“一線潮”表現(xiàn)的頻散效應最強.潮景的形成標志著涌潮進入不同發(fā)展階段,“交叉潮”出現(xiàn)在涌潮初期,F(xiàn)roude數(shù)接近1,涌潮高度在0.5 m左右;涌潮由大缺口進入鹽官河段之后形成“一線潮”,進入涌潮發(fā)展階段,并達到最大規(guī)模,“一線潮”的Froude保持在1.45以上,涌潮高度最大可達3 m,且涌潮的頻散作用較強;隨后涌潮到達老鹽倉形成“回頭潮”,同時涌潮強度被減弱,繼續(xù)上溯的涌潮進入三工段河道區(qū)域形成規(guī)模較小的“一線潮”,F(xiàn)roude數(shù)在1.3~1.4之間,涌潮高度約為2 m;七堡之后,上溯的涌潮進入消亡階段,F(xiàn)roude數(shù)減小至1.3以下,涌潮高度約為1 m.

    (3)錢塘江涌潮過程中,F(xiàn)roude數(shù)和涌潮高度呈正相關,涌潮形成至進入消亡階段,F(xiàn)roude數(shù)隨涌潮高度呈現(xiàn)逐漸增大再逐漸減小的規(guī)律.尖山河段處涌潮Froude數(shù)接近1,鹽官河段涌潮的Froude數(shù)最大可達1.54,老鹽倉之后,涌潮的Froude數(shù)逐漸減小.

    (4)Boussinesq型方程可以很好地再現(xiàn)錢塘江河口涌潮的形成、發(fā)展和消亡階段,對于鹽官河段處這類具有較強頻散效應的涌潮,展現(xiàn)了涌潮中的二次自由面起伏現(xiàn)象.

    我們在前人工作的基礎上,進一步發(fā)展了錢塘江涌潮數(shù)值模型,使用考慮非線性、頻散以及耗散作用的Boussinesq型方程對錢塘江河口涌潮進行模擬,模擬結果與實際觀測中體現(xiàn)的特征相吻合.

    在后續(xù)的工作中,考慮使用Boussinesq型方程結合更高精度的地形數(shù)據(jù)以及實時的涌潮觀測資料,可以運用在鞏固堤防、維護河道、保護涌潮等數(shù)值模擬方面的工作,推進日后對錢塘江涌潮的保護和深入研究探索.同時在強烈風暴作用下,風浪過程對水位的影響也十分重要(張舒羽和潘存鴻,2013;Luettich and Westerink,2015),進一步工作中可以在本文使用的數(shù)值模式基礎上考慮風浪與潮流的相互作用,增強涌潮防災減災的實用性價值.

    致謝感謝中國科學院大學計算地球動力學重點實驗室提供的計算平臺.

    猜你喜歡
    鹽官潮頭錢塘江
    為什么錢塘江的浪潮格外壯觀
    張立勇:勇立潮頭 奮楫者先
    華人時刊(2022年13期)2022-10-27 08:55:30
    我在錢塘江邊長大
    青年文學家(2022年1期)2022-03-11 12:31:09
    國家在場:從清代滇南鹽官營看國家邊疆治理
    錢塘江觀潮
    小讀者(2021年2期)2021-03-29 05:03:18
    番禺“鹽官廚”釋讀
    廣州文博(2020年0期)2020-06-09 05:14:10
    子文同學的一篇發(fā)掘日記與廣東漢代“鹽官”
    廣州文博(2020年0期)2020-06-09 05:14:04
    也談“番禺鹽官”
    廣州文博(2020年0期)2020-06-09 05:14:04
    勇立潮頭唱大風
    人大建設(2018年6期)2018-08-16 07:23:26
    “燃”在此刻,“燃”在潮頭
    人妻人人澡人人爽人人| 制服人妻中文乱码| 久久狼人影院| 久久久久久久国产电影| 一级a爱视频在线免费观看| www日本在线高清视频| 91九色精品人成在线观看| 真人做人爱边吃奶动态| 91老司机精品| 亚洲 欧美一区二区三区| 久久影院123| 国产日韩一区二区三区精品不卡| 搡老乐熟女国产| 免费一级毛片在线播放高清视频 | 老司机午夜福利在线观看视频 | 一边摸一边抽搐一进一出视频| 日韩欧美国产一区二区入口| 国产免费现黄频在线看| 免费在线观看黄色视频的| 青春草亚洲视频在线观看| 下体分泌物呈黄色| av在线老鸭窝| 黄频高清免费视频| 热99re8久久精品国产| 我的亚洲天堂| 久热爱精品视频在线9| 男人添女人高潮全过程视频| 欧美xxⅹ黑人| 久久久久视频综合| 在线十欧美十亚洲十日本专区| 国产精品影院久久| 亚洲性夜色夜夜综合| 午夜老司机福利片| 午夜日韩欧美国产| 国产一区有黄有色的免费视频| 日韩欧美免费精品| 久久久国产一区二区| 热99国产精品久久久久久7| 在线观看人妻少妇| 俄罗斯特黄特色一大片| av片东京热男人的天堂| 免费在线观看视频国产中文字幕亚洲 | 蜜桃在线观看..| 中国美女看黄片| 一级毛片电影观看| 在线观看www视频免费| 一区福利在线观看| 久久性视频一级片| 真人做人爱边吃奶动态| 嫩草影视91久久| 性高湖久久久久久久久免费观看| 亚洲av男天堂| 老熟妇乱子伦视频在线观看 | 精品久久久久久电影网| 高清欧美精品videossex| 精品亚洲成国产av| 国产成人精品无人区| 人妻久久中文字幕网| 又黄又粗又硬又大视频| 一区二区日韩欧美中文字幕| 欧美精品一区二区免费开放| 久久国产亚洲av麻豆专区| 亚洲国产欧美日韩在线播放| 国产亚洲欧美在线一区二区| 久久这里只有精品19| 国产成人免费无遮挡视频| 蜜桃国产av成人99| 女性生殖器流出的白浆| 日韩制服骚丝袜av| 亚洲va日本ⅴa欧美va伊人久久 | 欧美日韩亚洲国产一区二区在线观看 | 国产成人精品久久二区二区91| 中文精品一卡2卡3卡4更新| av不卡在线播放| 丝袜喷水一区| 天天操日日干夜夜撸| 一区二区av电影网| 亚洲欧美一区二区三区久久| 国产精品麻豆人妻色哟哟久久| 中文字幕人妻丝袜制服| av一本久久久久| 国产视频一区二区在线看| 亚洲国产精品999| 亚洲五月婷婷丁香| 在线观看免费高清a一片| www.自偷自拍.com| 老司机影院毛片| 日韩制服骚丝袜av| 国产精品一区二区精品视频观看| 中文字幕人妻丝袜一区二区| 两人在一起打扑克的视频| 女人爽到高潮嗷嗷叫在线视频| 国产成人系列免费观看| 欧美亚洲 丝袜 人妻 在线| 国产一区二区 视频在线| 国产成人欧美| 性色av乱码一区二区三区2| 91麻豆av在线| 电影成人av| 黄色视频在线播放观看不卡| 欧美日韩亚洲高清精品| 免费av中文字幕在线| 天天影视国产精品| 亚洲国产成人一精品久久久| 亚洲午夜精品一区,二区,三区| 人人妻人人添人人爽欧美一区卜| 色婷婷av一区二区三区视频| 久久九九热精品免费| 亚洲av日韩精品久久久久久密| 久久国产精品男人的天堂亚洲| 黄网站色视频无遮挡免费观看| 午夜日韩欧美国产| 啦啦啦在线免费观看视频4| 首页视频小说图片口味搜索| 亚洲人成电影免费在线| 在线观看免费午夜福利视频| 国产精品 国内视频| 天天添夜夜摸| 在线精品无人区一区二区三| 一区二区日韩欧美中文字幕| 999精品在线视频| 欧美精品亚洲一区二区| 女人爽到高潮嗷嗷叫在线视频| 国产日韩一区二区三区精品不卡| 欧美中文综合在线视频| 爱豆传媒免费全集在线观看| 色婷婷av一区二区三区视频| 一区在线观看完整版| bbb黄色大片| 欧美亚洲 丝袜 人妻 在线| 亚洲国产欧美一区二区综合| 男女免费视频国产| 精品国产乱码久久久久久男人| 日韩,欧美,国产一区二区三区| 热99久久久久精品小说推荐| 美女午夜性视频免费| 天天躁夜夜躁狠狠躁躁| 国产欧美日韩一区二区三区在线| 日韩熟女老妇一区二区性免费视频| 女人被躁到高潮嗷嗷叫费观| 亚洲欧美日韩另类电影网站| 色老头精品视频在线观看| 亚洲va日本ⅴa欧美va伊人久久 | 在线观看www视频免费| 老熟妇乱子伦视频在线观看 | 成人国产一区最新在线观看| 好男人电影高清在线观看| 婷婷成人精品国产| 一边摸一边做爽爽视频免费| 狂野欧美激情性bbbbbb| 欧美日韩精品网址| 色老头精品视频在线观看| 无限看片的www在线观看| 91精品伊人久久大香线蕉| 真人做人爱边吃奶动态| 青青草视频在线视频观看| 精品第一国产精品| 国产精品一区二区免费欧美 | 国产av一区二区精品久久| 黄片小视频在线播放| 国产一卡二卡三卡精品| 一二三四社区在线视频社区8| 国产精品九九99| 国产三级黄色录像| 91成年电影在线观看| 成人免费观看视频高清| 大香蕉久久网| 免费高清在线观看日韩| 最近中文字幕2019免费版| 久久国产精品大桥未久av| 91精品三级在线观看| 国产亚洲欧美在线一区二区| 精品人妻熟女毛片av久久网站| 国产91精品成人一区二区三区 | 97在线人人人人妻| avwww免费| 少妇精品久久久久久久| 日本欧美视频一区| 日韩欧美一区二区三区在线观看 | av有码第一页| 国产一区二区激情短视频 | 欧美人与性动交α欧美精品济南到| 亚洲一卡2卡3卡4卡5卡精品中文| 精品一区二区三区四区五区乱码| 久久久国产一区二区| 真人做人爱边吃奶动态| 999久久久精品免费观看国产| 丝袜脚勾引网站| 国产精品一区二区精品视频观看| 女人爽到高潮嗷嗷叫在线视频| 亚洲精品第二区| 一区二区三区激情视频| 91成人精品电影| 久久毛片免费看一区二区三区| 国产一区二区三区综合在线观看| 亚洲国产精品一区二区三区在线| 国产欧美日韩一区二区精品| 熟女少妇亚洲综合色aaa.| 日日夜夜操网爽| 精品少妇一区二区三区视频日本电影| 国产日韩欧美视频二区| 啦啦啦免费观看视频1| 欧美激情 高清一区二区三区| 亚洲欧美一区二区三区黑人| 1024视频免费在线观看| 国产97色在线日韩免费| 国产1区2区3区精品| 少妇 在线观看| 久久国产精品男人的天堂亚洲| 亚洲av日韩精品久久久久久密| 日韩大码丰满熟妇| 中文精品一卡2卡3卡4更新| 午夜91福利影院| 超碰97精品在线观看| 国产av一区二区精品久久| 国产欧美亚洲国产| 亚洲美女黄色视频免费看| 欧美日韩亚洲国产一区二区在线观看 | 99国产精品一区二区蜜桃av | 日韩欧美国产一区二区入口| 高清视频免费观看一区二区| 啪啪无遮挡十八禁网站| 成人国产av品久久久| 欧美黄色淫秽网站| 国产免费视频播放在线视频| 少妇粗大呻吟视频| 天天影视国产精品| 国产精品一区二区在线观看99| 国产欧美日韩一区二区精品| 亚洲精品国产一区二区精华液| 在线观看免费高清a一片| 欧美成狂野欧美在线观看| 在线观看免费午夜福利视频| 国产精品一区二区在线不卡| 亚洲国产精品成人久久小说| 久久久国产欧美日韩av| 国产一卡二卡三卡精品| 中文字幕另类日韩欧美亚洲嫩草| 成人黄色视频免费在线看| 国产高清国产精品国产三级| 色老头精品视频在线观看| 天天躁狠狠躁夜夜躁狠狠躁| 精品人妻熟女毛片av久久网站| 久久天躁狠狠躁夜夜2o2o| 人妻久久中文字幕网| 人妻人人澡人人爽人人| 欧美亚洲日本最大视频资源| 亚洲第一青青草原| 一级,二级,三级黄色视频| 大片电影免费在线观看免费| 天天躁日日躁夜夜躁夜夜| 精品亚洲成国产av| 亚洲精品国产精品久久久不卡| 欧美日韩视频精品一区| 久久久精品94久久精品| 老司机在亚洲福利影院| 精品亚洲成a人片在线观看| 老司机影院毛片| 亚洲一卡2卡3卡4卡5卡精品中文| 午夜成年电影在线免费观看| 不卡av一区二区三区| 巨乳人妻的诱惑在线观看| 国产精品自产拍在线观看55亚洲 | 美女福利国产在线| 亚洲精品第二区| 一区福利在线观看| 亚洲精品中文字幕一二三四区 | 最新在线观看一区二区三区| 精品国内亚洲2022精品成人 | 美女高潮到喷水免费观看| 香蕉国产在线看| 久久ye,这里只有精品| 97精品久久久久久久久久精品| 久久精品熟女亚洲av麻豆精品| 精品久久久久久久毛片微露脸 | 亚洲精品av麻豆狂野| 激情视频va一区二区三区| 18禁国产床啪视频网站| 亚洲黑人精品在线| 亚洲综合色网址| 色婷婷久久久亚洲欧美| 国产亚洲精品第一综合不卡| 欧美日韩亚洲综合一区二区三区_| 免费不卡黄色视频| 亚洲av成人一区二区三| 女人被躁到高潮嗷嗷叫费观| 黄色视频不卡| 午夜精品久久久久久毛片777| 天堂8中文在线网| 亚洲熟女毛片儿| 黄片小视频在线播放| 国产av又大| 一级毛片精品| 成人av一区二区三区在线看 | 欧美日韩亚洲国产一区二区在线观看 | 大片免费播放器 马上看| 国产高清视频在线播放一区 | 成人18禁高潮啪啪吃奶动态图| 秋霞在线观看毛片| 国产主播在线观看一区二区| 啪啪无遮挡十八禁网站| 午夜成年电影在线免费观看| 一二三四在线观看免费中文在| 欧美老熟妇乱子伦牲交| 天天操日日干夜夜撸| √禁漫天堂资源中文www| 97精品久久久久久久久久精品| 久久亚洲国产成人精品v| 国产精品成人在线| 精品国内亚洲2022精品成人 | av免费在线观看网站| 成年av动漫网址| 老汉色av国产亚洲站长工具| 亚洲av国产av综合av卡| 国产高清国产精品国产三级| 又紧又爽又黄一区二区| 美女中出高潮动态图| 亚洲七黄色美女视频| 女人精品久久久久毛片| 精品卡一卡二卡四卡免费| 精品亚洲乱码少妇综合久久| 国产欧美日韩精品亚洲av| 天天操日日干夜夜撸| 午夜激情久久久久久久| av不卡在线播放| 国产日韩欧美在线精品| 满18在线观看网站| 日韩一卡2卡3卡4卡2021年| 成人18禁高潮啪啪吃奶动态图| 九色亚洲精品在线播放| 欧美午夜高清在线| 欧美人与性动交α欧美精品济南到| 欧美久久黑人一区二区| 超碰成人久久| 日韩电影二区| 视频在线观看一区二区三区| 男女高潮啪啪啪动态图| 国产精品久久久久成人av| 久久久精品免费免费高清| 久久久久国内视频| 欧美大码av| 成年av动漫网址| 久久久水蜜桃国产精品网| 久久九九热精品免费| 91大片在线观看| 纯流量卡能插随身wifi吗| 狠狠婷婷综合久久久久久88av| 久久中文看片网| 日本一区二区免费在线视频| 国产免费视频播放在线视频| 午夜影院在线不卡| 免费观看a级毛片全部| 我要看黄色一级片免费的| 国产精品国产av在线观看| 久久九九热精品免费| 色播在线永久视频| 国产亚洲一区二区精品| 宅男免费午夜| 69精品国产乱码久久久| 人妻人人澡人人爽人人| 欧美变态另类bdsm刘玥| 欧美日韩视频精品一区| 一区二区三区乱码不卡18| 久热这里只有精品99| 黄色片一级片一级黄色片| 50天的宝宝边吃奶边哭怎么回事| 亚洲av男天堂| 国产成人a∨麻豆精品| 亚洲精品久久成人aⅴ小说| 操出白浆在线播放| 欧美激情久久久久久爽电影 | 亚洲成av片中文字幕在线观看| 欧美日韩一级在线毛片| 精品少妇内射三级| 蜜桃在线观看..| 日本vs欧美在线观看视频| 国产亚洲精品一区二区www | 亚洲激情五月婷婷啪啪| 国产av国产精品国产| 涩涩av久久男人的天堂| 国产福利在线免费观看视频| 两个人看的免费小视频| 亚洲三区欧美一区| 啦啦啦 在线观看视频| 欧美黄色淫秽网站| 亚洲熟女毛片儿| 国产av精品麻豆| 久久午夜综合久久蜜桃| 日本五十路高清| 汤姆久久久久久久影院中文字幕| 99久久综合免费| 美女午夜性视频免费| 大香蕉久久成人网| 亚洲va日本ⅴa欧美va伊人久久 | 久久久精品94久久精品| 又黄又粗又硬又大视频| 国产精品香港三级国产av潘金莲| 两性夫妻黄色片| 精品少妇一区二区三区视频日本电影| 精品福利永久在线观看| 热99国产精品久久久久久7| 久久中文字幕一级| 咕卡用的链子| 国产精品一区二区在线观看99| 精品国内亚洲2022精品成人 | 久久久精品国产亚洲av高清涩受| 黄色 视频免费看| 亚洲五月婷婷丁香| 在线十欧美十亚洲十日本专区| 国产一区二区三区在线臀色熟女 | 亚洲精品第二区| 国产成人精品在线电影| 男女高潮啪啪啪动态图| av天堂在线播放| 免费观看a级毛片全部| 97精品久久久久久久久久精品| 免费观看av网站的网址| 亚洲色图综合在线观看| 一区二区日韩欧美中文字幕| 亚洲免费av在线视频| 老司机靠b影院| 亚洲av国产av综合av卡| 久久精品亚洲熟妇少妇任你| 女性生殖器流出的白浆| 亚洲国产中文字幕在线视频| 久久久国产欧美日韩av| 亚洲人成电影免费在线| 亚洲熟女毛片儿| 69av精品久久久久久 | 日韩 亚洲 欧美在线| 日韩欧美免费精品| 欧美激情久久久久久爽电影 | 亚洲欧美日韩高清在线视频 | e午夜精品久久久久久久| 精品国产乱码久久久久久男人| 一本综合久久免费| a级毛片在线看网站| 天堂俺去俺来也www色官网| 欧美日韩视频精品一区| 男女午夜视频在线观看| 亚洲欧洲日产国产| 一二三四社区在线视频社区8| 亚洲中文日韩欧美视频| 成人午夜高清在线视频| 动漫黄色视频在线观看| bbb黄色大片| 亚洲欧美日韩无卡精品| 国产人伦9x9x在线观看| 19禁男女啪啪无遮挡网站| 国产亚洲精品综合一区在线观看 | 国产午夜福利久久久久久| 欧美人与性动交α欧美精品济南到| 久久久久亚洲av毛片大全| 精品午夜福利视频在线观看一区| 亚洲,欧美精品.| 正在播放国产对白刺激| 女生性感内裤真人,穿戴方法视频| 欧美日韩国产亚洲二区| 免费看十八禁软件| 亚洲专区中文字幕在线| 国产v大片淫在线免费观看| 国产欧美日韩一区二区三| 国产99久久九九免费精品| 国产精品久久久人人做人人爽| 国产精品一区二区免费欧美| 嫩草影院精品99| 亚洲av五月六月丁香网| 久久九九热精品免费| 神马国产精品三级电影在线观看 | 国产免费av片在线观看野外av| 精品久久久久久久久久免费视频| 极品教师在线免费播放| 一级黄色大片毛片| 黄色片一级片一级黄色片| 日韩精品中文字幕看吧| 一个人观看的视频www高清免费观看 | 欧美zozozo另类| 成人特级黄色片久久久久久久| 一进一出抽搐动态| 亚洲国产精品999在线| 久久伊人香网站| 床上黄色一级片| 搞女人的毛片| 午夜福利视频1000在线观看| 波多野结衣高清无吗| 校园春色视频在线观看| 国产av不卡久久| 欧美日韩国产亚洲二区| 十八禁人妻一区二区| www.自偷自拍.com| 午夜成年电影在线免费观看| 国产精品九九99| 老司机深夜福利视频在线观看| 俺也久久电影网| 欧美又色又爽又黄视频| 精品久久久久久成人av| 19禁男女啪啪无遮挡网站| 99国产精品一区二区三区| 男人的好看免费观看在线视频 | 欧美另类亚洲清纯唯美| АⅤ资源中文在线天堂| 欧美日韩国产亚洲二区| 男人的好看免费观看在线视频 | 久久久久精品国产欧美久久久| 国产一区二区激情短视频| 欧美激情久久久久久爽电影| 欧美人与性动交α欧美精品济南到| 欧美av亚洲av综合av国产av| 亚洲人成网站高清观看| 久久中文字幕一级| 日本黄色视频三级网站网址| e午夜精品久久久久久久| 久久天堂一区二区三区四区| 中文字幕高清在线视频| 日本 欧美在线| 91字幕亚洲| www.精华液| 国产三级黄色录像| 国产黄片美女视频| 别揉我奶头~嗯~啊~动态视频| 亚洲成人精品中文字幕电影| 亚洲在线自拍视频| 精品免费久久久久久久清纯| 亚洲精品在线美女| 亚洲第一欧美日韩一区二区三区| 中文字幕人妻丝袜一区二区| 亚洲自偷自拍图片 自拍| 久久九九热精品免费| 亚洲欧美日韩无卡精品| 久久久精品大字幕| 亚洲欧美一区二区三区黑人| 悠悠久久av| 欧美一级毛片孕妇| 国产av不卡久久| 亚洲一区中文字幕在线| 精品熟女少妇八av免费久了| 国产不卡一卡二| 丁香六月欧美| 色精品久久人妻99蜜桃| 午夜免费成人在线视频| 九色国产91popny在线| 在线观看美女被高潮喷水网站 | 日韩高清综合在线| 免费看a级黄色片| 精品第一国产精品| 制服丝袜大香蕉在线| 亚洲中文字幕一区二区三区有码在线看 | 亚洲欧美精品综合久久99| 天堂动漫精品| 色老头精品视频在线观看| 操出白浆在线播放| 久久99热这里只有精品18| 看黄色毛片网站| 久久久久国内视频| 久久人人精品亚洲av| 久久久精品大字幕| 大型av网站在线播放| 999久久久国产精品视频| svipshipincom国产片| 欧美一区二区国产精品久久精品 | 性色av乱码一区二区三区2| 久久久久久免费高清国产稀缺| 国产v大片淫在线免费观看| 黄色视频不卡| 变态另类成人亚洲欧美熟女| 在线看三级毛片| 免费电影在线观看免费观看| 老司机深夜福利视频在线观看| 99在线视频只有这里精品首页| 一区福利在线观看| 国产亚洲精品综合一区在线观看 | 2021天堂中文幕一二区在线观| 黑人欧美特级aaaaaa片| av视频在线观看入口| 天天躁夜夜躁狠狠躁躁| 日韩欧美国产在线观看| 日韩欧美免费精品| 12—13女人毛片做爰片一| √禁漫天堂资源中文www| 最近视频中文字幕2019在线8| 在线观看www视频免费| 成年女人毛片免费观看观看9| 变态另类丝袜制服| 变态另类成人亚洲欧美熟女| 欧美3d第一页| 亚洲精品粉嫩美女一区| 亚洲精品一卡2卡三卡4卡5卡| e午夜精品久久久久久久| 一二三四社区在线视频社区8| 欧美日韩精品网址| 桃色一区二区三区在线观看| 长腿黑丝高跟| 国产精品免费视频内射| 欧美黄色淫秽网站| 床上黄色一级片| 日本成人三级电影网站| 成年免费大片在线观看| 亚洲精品美女久久av网站| 国产精品免费视频内射| 欧美中文日本在线观看视频| 国产三级黄色录像| 亚洲中文字幕日韩| 亚洲av电影不卡..在线观看| 精品久久久久久久毛片微露脸| 久久久久久免费高清国产稀缺| 美女扒开内裤让男人捅视频| 国产精品精品国产色婷婷| 人妻夜夜爽99麻豆av| 国产av一区在线观看免费| 男女之事视频高清在线观看| 91字幕亚洲| 精品久久久久久久久久免费视频| 色播亚洲综合网| 久久人人精品亚洲av| 久久 成人 亚洲|