• <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
    “燃”在此刻,“燃”在潮頭
    国产午夜精品一二区理论片| 亚洲国产日韩欧美精品在线观看| 一二三四中文在线观看免费高清| 亚洲国产精品专区欧美| 精品熟女少妇av免费看| 国产麻豆成人av免费视频| 肉色欧美久久久久久久蜜桃 | 亚洲成人中文字幕在线播放| 十八禁网站网址无遮挡 | 水蜜桃什么品种好| 国产黄频视频在线观看| 久久久久国产网址| 三级毛片av免费| 精品久久久精品久久久| 激情五月婷婷亚洲| 热99在线观看视频| 免费高清在线观看视频在线观看| 日韩国内少妇激情av| 精品一区二区三卡| 亚洲精品成人av观看孕妇| 69人妻影院| 国产av不卡久久| 91av网一区二区| 国产精品一区二区三区四区免费观看| 久久久久免费精品人妻一区二区| 色尼玛亚洲综合影院| 午夜精品在线福利| 欧美性感艳星| 日本猛色少妇xxxxx猛交久久| 精华霜和精华液先用哪个| 边亲边吃奶的免费视频| 少妇高潮的动态图| 黄色欧美视频在线观看| 欧美日韩亚洲高清精品| 成人性生交大片免费视频hd| 午夜福利视频精品| 精品一区在线观看国产| 日本与韩国留学比较| 欧美日韩一区二区视频在线观看视频在线 | 亚洲精品日本国产第一区| 久久99精品国语久久久| 国产在线男女| 亚洲av电影不卡..在线观看| 大香蕉97超碰在线| 久久综合国产亚洲精品| 免费播放大片免费观看视频在线观看| 国产高潮美女av| 美女被艹到高潮喷水动态| 国产精品久久久久久久久免| 18禁裸乳无遮挡免费网站照片| 欧美日韩视频高清一区二区三区二| 99久久中文字幕三级久久日本| 亚洲乱码一区二区免费版| 国产高清不卡午夜福利| 尾随美女入室| 天堂网av新在线| 精品久久久久久成人av| 最近的中文字幕免费完整| 天美传媒精品一区二区| av又黄又爽大尺度在线免费看| 日韩中字成人| 日韩av不卡免费在线播放| 黄色日韩在线| 97超碰精品成人国产| 又黄又爽又刺激的免费视频.| 亚洲一区高清亚洲精品| 老司机影院成人| 午夜福利视频精品| 免费看美女性在线毛片视频| 国产精品一区二区三区四区久久| 亚洲欧洲国产日韩| 亚洲精品亚洲一区二区| 国产精品.久久久| 久久久精品免费免费高清| 欧美极品一区二区三区四区| 少妇人妻精品综合一区二区| 欧美xxxx黑人xx丫x性爽| 午夜亚洲福利在线播放| 亚洲内射少妇av| www.色视频.com| 免费看a级黄色片| 一级毛片黄色毛片免费观看视频| 欧美丝袜亚洲另类| 插阴视频在线观看视频| 久久久久久久久久黄片| 国产精品一及| 只有这里有精品99| 边亲边吃奶的免费视频| 欧美日韩综合久久久久久| 久久精品国产亚洲av涩爱| 日韩电影二区| 日日啪夜夜爽| 极品教师在线视频| 一级黄片播放器| 在线免费十八禁| 男插女下体视频免费在线播放| 国产成人a∨麻豆精品| 乱人视频在线观看| 色综合站精品国产| 精品久久久久久久久亚洲| 91av网一区二区| 亚洲熟女精品中文字幕| 美女大奶头视频| 在线观看一区二区三区| 中文欧美无线码| 草草在线视频免费看| 爱豆传媒免费全集在线观看| 免费高清在线观看视频在线观看| 永久免费av网站大全| 午夜福利网站1000一区二区三区| 亚洲成色77777| 免费av不卡在线播放| 免费大片黄手机在线观看| 亚洲av二区三区四区| 九九在线视频观看精品| 99re6热这里在线精品视频| 免费av观看视频| 日韩国内少妇激情av| 国产v大片淫在线免费观看| 狠狠精品人妻久久久久久综合| 精品久久久久久成人av| 在线观看人妻少妇| 国产免费又黄又爽又色| 观看免费一级毛片| 国产精品麻豆人妻色哟哟久久 | 国产一区二区在线观看日韩| 亚洲国产精品专区欧美| 国产成人福利小说| 天美传媒精品一区二区| 好男人视频免费观看在线| 大香蕉97超碰在线| 国产美女午夜福利| 午夜激情福利司机影院| 一级毛片黄色毛片免费观看视频| 国产成人freesex在线| 日本熟妇午夜| 国产精品1区2区在线观看.| 国产老妇女一区| 亚洲欧美一区二区三区黑人 | 麻豆精品久久久久久蜜桃| 男人和女人高潮做爰伦理| 久久久久精品性色| 婷婷色av中文字幕| 99热这里只有精品一区| 国产乱人偷精品视频| 国产高清三级在线| 麻豆成人av视频| 插逼视频在线观看| 最新中文字幕久久久久| 我要看日韩黄色一级片| 一个人免费在线观看电影| 免费观看无遮挡的男女| 少妇猛男粗大的猛烈进出视频 | 一个人免费在线观看电影| 91精品伊人久久大香线蕉| 男插女下体视频免费在线播放| 中文乱码字字幕精品一区二区三区 | 国产精品一二三区在线看| 如何舔出高潮| 男女那种视频在线观看| 国产淫语在线视频| 国产成人一区二区在线| 国产成人精品福利久久| 人妻夜夜爽99麻豆av| 2022亚洲国产成人精品| 91av网一区二区| 精品人妻视频免费看| 久久99热这里只有精品18| 男女边吃奶边做爰视频| 国产黄片视频在线免费观看| 亚洲综合精品二区| 秋霞在线观看毛片| 真实男女啪啪啪动态图| 日本熟妇午夜| 日韩欧美 国产精品| 少妇丰满av| 一级片'在线观看视频| 成人漫画全彩无遮挡| 欧美zozozo另类| 在现免费观看毛片| 国产乱来视频区| 亚洲国产精品国产精品| 80岁老熟妇乱子伦牲交| .国产精品久久| 午夜福利高清视频| 中文欧美无线码| 99热网站在线观看| 日韩强制内射视频| 国产精品久久久久久av不卡| 老司机影院成人| 如何舔出高潮| 国产探花在线观看一区二区| 亚洲性久久影院| 在线天堂最新版资源| 神马国产精品三级电影在线观看| 少妇丰满av| 精品久久国产蜜桃| 美女xxoo啪啪120秒动态图| 好男人在线观看高清免费视频| 欧美成人午夜免费资源| 一级毛片aaaaaa免费看小| 纵有疾风起免费观看全集完整版 | 亚洲国产最新在线播放| 一级黄片播放器| 亚洲av在线观看美女高潮| 久久久欧美国产精品| 插逼视频在线观看| 午夜福利在线在线| 国产 亚洲一区二区三区 | 欧美精品国产亚洲| 国精品久久久久久国模美| 国产精品嫩草影院av在线观看| 中文字幕免费在线视频6| av专区在线播放| 国产精品1区2区在线观看.| 亚洲欧美日韩无卡精品| 国产亚洲91精品色在线| 69人妻影院| 亚洲国产精品成人久久小说| 国产单亲对白刺激| 内地一区二区视频在线| 伊人久久国产一区二区| 午夜免费男女啪啪视频观看| 老司机影院毛片| 久久99精品国语久久久| 日韩欧美三级三区| 校园人妻丝袜中文字幕| 好男人在线观看高清免费视频| 成人二区视频| 波多野结衣巨乳人妻| 网址你懂的国产日韩在线| 久久精品人妻少妇| 精品久久久久久电影网| 成人漫画全彩无遮挡| 午夜福利视频1000在线观看| 熟女人妻精品中文字幕| 噜噜噜噜噜久久久久久91| 亚洲国产欧美人成| 18禁在线无遮挡免费观看视频| 在线免费十八禁| 尤物成人国产欧美一区二区三区| 精品久久久久久久人妻蜜臀av| 日本免费在线观看一区| 精品久久久久久久久亚洲| 成人漫画全彩无遮挡| 国产精品久久久久久精品电影小说 | 久久久久久久久中文| .国产精品久久| 男女下面进入的视频免费午夜| 丰满乱子伦码专区| 男插女下体视频免费在线播放| 久热久热在线精品观看| 18+在线观看网站| 肉色欧美久久久久久久蜜桃 | 日本与韩国留学比较| 久久久国产一区二区| 精品久久久精品久久久| 97精品久久久久久久久久精品| 真实男女啪啪啪动态图| 亚洲最大成人中文| 午夜福利在线在线| 天堂√8在线中文| 亚洲综合精品二区| 欧美成人午夜免费资源| 岛国毛片在线播放| 日产精品乱码卡一卡2卡三| 亚洲自偷自拍三级| 毛片一级片免费看久久久久| 久久久久性生活片| 可以在线观看毛片的网站| 亚洲成人av在线免费| 亚洲成人中文字幕在线播放| av一本久久久久| 好男人在线观看高清免费视频| 亚洲精品一区蜜桃| 免费看光身美女| av在线观看视频网站免费| 久久人人爽人人爽人人片va| 一本久久精品| 亚洲不卡免费看| 一个人看视频在线观看www免费| 在线观看av片永久免费下载| 国产精品国产三级国产av玫瑰| 免费观看无遮挡的男女| 超碰97精品在线观看| 欧美变态另类bdsm刘玥| 国产午夜福利久久久久久| 热99在线观看视频| 免费观看性生交大片5| 天堂俺去俺来也www色官网 | 最近中文字幕2019免费版| 日本午夜av视频| 搡老乐熟女国产| 国产精品久久久久久精品电影小说 | 熟妇人妻久久中文字幕3abv| 久久99精品国语久久久| 男女边吃奶边做爰视频| 精品一区二区三区视频在线| 狠狠精品人妻久久久久久综合| 国产精品99久久久久久久久| 国产三级在线视频| 久久久久久久午夜电影| 久99久视频精品免费| 国产一区有黄有色的免费视频 | 日韩欧美精品v在线| 免费看美女性在线毛片视频| 国产69精品久久久久777片| 亚洲人成网站在线播| 夜夜看夜夜爽夜夜摸| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 国产免费福利视频在线观看| 久久久久久久久久黄片| 一本久久精品| 精品久久久久久久久久久久久| 精品酒店卫生间| 日韩伦理黄色片| 校园人妻丝袜中文字幕| 晚上一个人看的免费电影| 久久久久久久久久成人| 午夜福利在线在线| 丰满乱子伦码专区| 狂野欧美白嫩少妇大欣赏| 天堂av国产一区二区熟女人妻| 天天躁夜夜躁狠狠久久av| 国产一区二区亚洲精品在线观看| 久久久久网色| 日韩电影二区| 亚洲精品aⅴ在线观看| 精品午夜福利在线看| 五月天丁香电影| 亚洲三级黄色毛片| 国产极品天堂在线| 亚洲av免费在线观看| 亚洲精品国产av成人精品| 一级爰片在线观看| 免费观看在线日韩| 日韩欧美三级三区| 免费观看的影片在线观看| 亚洲图色成人| 波野结衣二区三区在线| 男女下面进入的视频免费午夜| 女人十人毛片免费观看3o分钟| 蜜桃亚洲精品一区二区三区| 人妻一区二区av| 欧美日韩精品成人综合77777| 日本av手机在线免费观看| 亚洲成人一二三区av| 高清视频免费观看一区二区 | 亚洲国产最新在线播放| 久久久久久久久久久丰满| 大陆偷拍与自拍| av又黄又爽大尺度在线免费看| 午夜福利网站1000一区二区三区| 亚洲精品久久久久久婷婷小说| 欧美性猛交╳xxx乱大交人| 成年免费大片在线观看| 日日摸夜夜添夜夜添av毛片| 久久久亚洲精品成人影院| av在线天堂中文字幕| 狠狠精品人妻久久久久久综合| 亚洲欧美日韩东京热| 国产成人精品婷婷| 秋霞在线观看毛片| 在线天堂最新版资源| 一级毛片电影观看| 青春草国产在线视频| 日韩,欧美,国产一区二区三区| 国产一区二区三区av在线| 久久99热这里只频精品6学生| 中文字幕亚洲精品专区| 伊人久久精品亚洲午夜| 九九久久精品国产亚洲av麻豆| 国产精品嫩草影院av在线观看| 免费看光身美女| 久久久成人免费电影| 少妇熟女欧美另类| 少妇熟女aⅴ在线视频| 在线观看一区二区三区| 嫩草影院入口| av在线老鸭窝| 永久免费av网站大全| av一本久久久久| 免费看日本二区| 亚洲av二区三区四区| 国产精品一及| 亚洲成人久久爱视频| 国产亚洲精品av在线| 汤姆久久久久久久影院中文字幕 | 成人av在线播放网站| 亚洲性久久影院| 国产精品久久久久久精品电影| 亚洲欧美一区二区三区国产| 欧美变态另类bdsm刘玥| 国产午夜精品久久久久久一区二区三区| 嫩草影院新地址| 国产乱人偷精品视频| 人妻制服诱惑在线中文字幕| 99久久中文字幕三级久久日本| 国产片特级美女逼逼视频| 国产伦理片在线播放av一区| 亚洲国产av新网站| 亚洲精品日本国产第一区| 精品少妇黑人巨大在线播放| 久久久久久久久大av| 最近2019中文字幕mv第一页| 黄色配什么色好看| 九色成人免费人妻av| 精品99又大又爽又粗少妇毛片| 91久久精品电影网| 日日摸夜夜添夜夜爱| 一级爰片在线观看| 国产亚洲5aaaaa淫片| 一区二区三区免费毛片| 日韩视频在线欧美| 少妇人妻一区二区三区视频| 久久午夜福利片| 久久精品熟女亚洲av麻豆精品 | 国产亚洲91精品色在线| 少妇人妻精品综合一区二区| 麻豆成人av视频| 日本wwww免费看| 国产片特级美女逼逼视频| 尾随美女入室| 亚洲18禁久久av| 国产91av在线免费观看| 成人毛片a级毛片在线播放| 乱人视频在线观看| 又爽又黄无遮挡网站| 少妇高潮的动态图| 精品一区二区三区人妻视频| 亚洲国产欧美人成| 午夜精品国产一区二区电影 | 成年av动漫网址| 成人毛片60女人毛片免费| 日本av手机在线免费观看| 国产色婷婷99| 日韩不卡一区二区三区视频在线| 亚洲人与动物交配视频| 精品久久久久久成人av| 国产女主播在线喷水免费视频网站 | 男人和女人高潮做爰伦理| 色综合色国产| 亚洲熟妇中文字幕五十中出| 最近的中文字幕免费完整| 国产男女超爽视频在线观看| 91久久精品国产一区二区三区| 亚洲不卡免费看| 午夜日本视频在线| 久久午夜福利片| 尾随美女入室| 国产亚洲91精品色在线| 亚洲人成网站在线播| 日本与韩国留学比较| 高清欧美精品videossex| 日韩欧美一区视频在线观看 | 国产一区二区三区av在线| 97热精品久久久久久| a级一级毛片免费在线观看| 男女下面进入的视频免费午夜| 少妇熟女欧美另类| 黄色日韩在线| 国产熟女欧美一区二区| 国产真实伦视频高清在线观看| av女优亚洲男人天堂| 亚洲电影在线观看av| 黑人高潮一二区| 亚洲综合色惰| 亚洲人成网站在线播| 亚洲国产精品sss在线观看| 国产视频内射| 国产极品天堂在线| 狂野欧美白嫩少妇大欣赏| 亚洲精品自拍成人| 亚洲精品日本国产第一区| 免费观看av网站的网址| 欧美日韩国产mv在线观看视频 | 免费不卡的大黄色大毛片视频在线观看 | 亚洲精品一区蜜桃| 国产一区二区三区av在线| 国产精品熟女久久久久浪| 亚洲经典国产精华液单| 我要看日韩黄色一级片| 国产日韩欧美在线精品| 搡女人真爽免费视频火全软件| 秋霞在线观看毛片| 欧美区成人在线视频| 午夜亚洲福利在线播放| 联通29元200g的流量卡| 久久久久精品久久久久真实原创| 秋霞在线观看毛片| 亚洲国产高清在线一区二区三| 日本爱情动作片www.在线观看| 免费少妇av软件| 少妇被粗大猛烈的视频| 在线免费观看不下载黄p国产| 国产精品久久久久久久电影| 成人午夜高清在线视频| 免费看不卡的av| 国产男人的电影天堂91| 国产91av在线免费观看| av免费观看日本| 午夜激情欧美在线| 亚洲欧美一区二区三区黑人 | 寂寞人妻少妇视频99o| 少妇猛男粗大的猛烈进出视频 | 国产伦精品一区二区三区四那| 国产精品一及| 中文字幕免费在线视频6| 亚洲精品aⅴ在线观看| 欧美激情在线99| 中文资源天堂在线| 乱人视频在线观看| 欧美激情国产日韩精品一区| 综合色av麻豆| 80岁老熟妇乱子伦牲交| 草草在线视频免费看| 亚洲av电影不卡..在线观看| 国产三级在线视频| 午夜激情福利司机影院| 91午夜精品亚洲一区二区三区| 日本av手机在线免费观看| 97超碰精品成人国产| 色综合站精品国产| a级一级毛片免费在线观看| 国产欧美另类精品又又久久亚洲欧美| 精品久久久久久久久久久久久| 亚洲成人中文字幕在线播放| 美女高潮的动态| 久久99精品国语久久久| 午夜激情欧美在线| 日韩人妻高清精品专区| 能在线免费观看的黄片| 国产女主播在线喷水免费视频网站 | 成年免费大片在线观看| 亚洲精品日韩在线中文字幕| 午夜爱爱视频在线播放| 你懂的网址亚洲精品在线观看| 欧美激情在线99| 色5月婷婷丁香| 亚洲精品国产成人久久av| 精品一区二区三卡| av又黄又爽大尺度在线免费看| 日韩欧美国产在线观看| 日韩精品有码人妻一区| 一级av片app| 亚洲精品,欧美精品| 国产人妻一区二区三区在| 亚洲电影在线观看av| 春色校园在线视频观看| 国产精品一区二区在线观看99 | a级毛色黄片| 久久99精品国语久久久| 国产精品伦人一区二区| 久久国产乱子免费精品| 超碰97精品在线观看| 汤姆久久久久久久影院中文字幕 | 精品久久久噜噜| 国产av不卡久久| 国产探花在线观看一区二区| 精品午夜福利在线看| 国产综合懂色| 免费av毛片视频| videos熟女内射| 亚洲av福利一区| 99久久精品一区二区三区| 亚洲国产欧美人成| av专区在线播放| 丰满乱子伦码专区| 一级毛片aaaaaa免费看小| xxx大片免费视频| 一级爰片在线观看| 国产亚洲91精品色在线| 欧美潮喷喷水| 高清av免费在线| 国产成人精品一,二区| 国产一区亚洲一区在线观看| 久久久a久久爽久久v久久| 亚洲综合色惰| 中文在线观看免费www的网站| 97超视频在线观看视频| 天堂网av新在线| 欧美精品国产亚洲| 成人一区二区视频在线观看| 国产精品三级大全| 久久久欧美国产精品| 听说在线观看完整版免费高清| 欧美激情国产日韩精品一区| 日韩av在线大香蕉| .国产精品久久| 国产在视频线精品| 免费观看在线日韩| 国产淫片久久久久久久久| 纵有疾风起免费观看全集完整版 | 日韩欧美国产在线观看| 亚州av有码| a级一级毛片免费在线观看| 色哟哟·www| 精品久久久久久久人妻蜜臀av| 国产精品爽爽va在线观看网站| 亚洲精华国产精华液的使用体验| 搡老乐熟女国产| 人妻一区二区av| 久久精品久久精品一区二区三区| 免费大片黄手机在线观看| 天天躁日日操中文字幕| 国产高清不卡午夜福利| 亚洲四区av| 精品人妻偷拍中文字幕| 秋霞在线观看毛片| 自拍偷自拍亚洲精品老妇| 婷婷色av中文字幕| 少妇丰满av| 尤物成人国产欧美一区二区三区| 美女黄网站色视频| 大又大粗又爽又黄少妇毛片口|