亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        基于渦動相關法的沉積物-水界面氧通量與水動力條件響應關系

        2020-05-06 07:30:26高學平郭曉雪孫博聞張袁寧
        水利學報 2020年3期
        關鍵詞:邊界層溶解氧時滯

        高學平,郭曉雪,孫博聞,張袁寧

        (天津大學 水利工程仿真與安全國家重點實驗室,天津 300072)

        1 研究背景

        沉積物-水界面(Sediment-Water Interface,SWI)作為水生態(tài)系統(tǒng)中的關鍵界面之一,是沉積物和水體之間物質(zhì)垂向交換的主要場所,關于其物質(zhì)通量的研究一直是國際上關注的熱點問題[1]。溶解氧作為評價水體水質(zhì)的常用指標,對于水生生物的生存、水體自凈功能的維持等起著關鍵作用。SWI氧通量由于被廣泛用于評估底棲生物的初級生產(chǎn)力、有機物礦化率,因此對研究水體物質(zhì)循環(huán)、富營養(yǎng)化治理、生態(tài)系統(tǒng)功能等都具有重要意義。

        目前關于SWI 氧通量的研究集中于分析通量與沉積物和水體中物質(zhì)含量的關系。研究表明,沉積物有機物質(zhì)含量、水體底部溶解氧濃度、泥沙粒徑和葉綠素含量等因素均會對SWI 氧通量產(chǎn)生影響[2-4]。潘延鑫等[5]對農(nóng)田排水溝的SWI 氧通量觀測發(fā)現(xiàn),上下游界面氧通量的差異可能與有機質(zhì)、鹽分含量及微生物活動等有關,但由于試驗過程中水體處于靜置狀態(tài),因此水動力條件同樣可能是影響界面氧傳輸?shù)闹匾蛩兀?-9]。Koopmans[10]和鄭陽華等[11]通過原位和試驗研究均發(fā)現(xiàn)SWI 氧通量隨水平流速增大相應增大,Scalo 等[12]在構(gòu)建氧通量代數(shù)模型中也將摩擦流速作為主要輸入?yún)?shù)。從水動力條件對氧通量的影響機制來看,目前有學者提出水動力條件可能通過控制擴散邊界層(DBL)厚度來實現(xiàn)對氧通量的影響[13-15],DBL 作為控制沉積物-水界面物質(zhì)垂向交換的主要瓶頸,與摩擦流速、雷諾數(shù)、Batchelor 尺度等水動力條件關系密切[9,16-18],但結(jié)論多為定性描述。

        對SWI 氧通量的測量一般可采用水底培養(yǎng)箱法、微電極剖面法和渦動相關法等。水底培養(yǎng)箱法通過分析封閉沉積物及其上覆水中溶解氧隨時間的變化規(guī)律評估氧通量,該方法影響了觀測區(qū)域與周圍水體間水流交換,Brink 等[19]通過內(nèi)部自帶的水流動力裝置模擬實際流動,依然難以還原真實水動力條件。微電極剖面法將微電極緩慢刺入沉積物內(nèi),根據(jù)溶解氧在沉積物-水界面的垂向分布得到氧通量,該方法雖然其垂直分辨率很高,但一般僅能獲取垂向梯度的氧通量信息,難以反映地形變化、生物活動等對氧通量的影響;該方法的測量結(jié)果還存在偶然性,R?y[20]等發(fā)現(xiàn)三維的氧通量測量結(jié)果比一維條件下高約10%。針對上述測量技術的不足,Berg 等[21]首次將渦動相關法應用于SWI 氧通量測量,通過直接測量靠近沉積物表面處水體中流速與溶解氧值獲得氧通量,可反映5 ~100 m2測量區(qū)域[22](也稱測量足跡)內(nèi)的氧通量信息。這一方法由于對水動力條件影響小,因此在研究氧通量與水動力條件響應關系方面顯示出獨特優(yōu)勢。

        本文基于渦動相關法理論基礎,采用溶解氧傳感器和聲學多普勒點式流速儀構(gòu)建了非侵入式渦動相關系統(tǒng)。通過室內(nèi)試驗對不同水平流速條件下溶解氧在沉積物-水界面的垂向分布進行觀測,獲得了擴散邊界層厚度;根據(jù)垂向流速與溶解氧濃度的實時測量結(jié)果得到了SWI 氧通量,并詳細介紹了氧通量求解過程及關鍵參數(shù)處理方法。將本文及相關研究中水動力條件、擴散邊界層厚度及氧通量進行擬合,得到了SWI 氧通量與不同水動力條件間的定量響應關系,成果可以為SWI 氧通量對水動力條件的響應機理研究提供參考。

        2 渦動相關法

        2.1 理論基礎渦動相關法最早可追溯到1950年代,起初該方法多用于分析物質(zhì)、動量、熱量等在陸-氣、海-氣邊界層的交換。2003年Berg 等[21]首次將其應用于沉積物-水界面的氧通量測量,此后陸續(xù)有學者在水底邊界層這一特定環(huán)境下對其理論進行改進,并開展了相關研究[15,23]。渦動相關法通過計算垂向流速與其它物理要素的協(xié)方差從而得到物質(zhì)通量,基于質(zhì)量守恒方程,溶質(zhì)在控制體中滿足:

        式中:C 為控制體內(nèi)溶質(zhì)濃度;uj=(u,v,w)為流速,j 為正交直角坐標系的三個方向;Dc為分子擴散系數(shù),m2/s;Sc為控制體內(nèi)溶質(zhì)的源匯項;為控制體內(nèi)溶質(zhì)濃度的變化速率。

        采用雷諾分解,將具有紊動特性的物理量分解為時均值和紊動值,即,代入式(1)得到:

        假設沉積物-水界面高度為z0,測點高度為zm,對式(3)兩端同時進行積分可得:

        由沉積物-水界面的物質(zhì)垂向傳輸方式(圖1)可知,在擴散邊界層(DBL)上部水體中物質(zhì)的垂向傳輸機制主要為湍流擴散,而在擴散邊界層內(nèi)分子擴散起主導作用,因此測點距離沉積物-水界面很近,因此基于上述分析可得:

        圖1 沉積物-水界面物質(zhì)垂向傳輸方式示意圖[1]

        式(5)說明沉積物-水界面處的物質(zhì)通量可以用靠近沉積物處水體中某點的垂向流速與溶質(zhì)濃度的協(xié)方差來表示。考慮到測量穩(wěn)定性,通常需進行較長時段測量,因此沉積物-水界面物質(zhì)通量:,其中n 為時段內(nèi)數(shù)據(jù)個數(shù)。

        圖2 試驗裝置及渦動相關系統(tǒng)示意圖(單位:cm)

        2.2 渦動相關系統(tǒng)構(gòu)建本文采用的渦動相關系統(tǒng)(圖2)由流速測量模塊和溶解氧測量模塊兩部分構(gòu)成。流速測量模塊選用聲學多普勒點式流速儀Vector6MHz(簡稱ADV),可以測量固定點的三維流速、流向、水壓、水溫、水深、波高等指標。ADV 內(nèi)置電源和數(shù)據(jù)存儲裝置,提供兩個模擬通道,采樣頻率為1 ~64 Hz,采樣體位于信號發(fā)射端正下方15 cm 處,為直徑15 mm,高度5 ~20 mm的近似圓柱水體;溶解氧測量模塊選用快速響應的ARO-EC 型溶解氧傳感器,可以測量固定點的溶解氧濃度及水溫。ARO-EC 材質(zhì)為鈦,尖端直徑12 mm,長164 mm,基于熒光壽命法進行溶解氧測量,在測量過程中不會引起信號漂移,可滿足長時間觀測需要。同時ARO-EC 為溶解氧溫度雙傳感器,響應時間均小于0.5 s,可實現(xiàn)溶解氧的快速矯正,提高測量準確性。ADV 與ARO-EC 通過水密電纜連接,ADV 通過電纜控制ARO-EC 工作并為其供電,ARO-EC 可將測量數(shù)據(jù)通過電纜傳輸給ADV 并進行保存,兩者協(xié)同工作,實現(xiàn)流速和溶解氧的同步測量。

        3 材料與方法

        3.1 試驗裝置試驗裝置如圖2所示,該試驗裝置為長×寬×高=9.3 m×0.8 m×1.2 m 的長方體有機玻璃水槽,水槽采用自循環(huán)系統(tǒng),分為上下兩層,下層用于蓄水,上層為試驗區(qū)域,中間用底板隔開,水體通過水泵在上下水箱中循環(huán)流動,試驗過程中無外界水流流入流出。水槽上部鋪設滑軌,并架設儀器固定支架,該支架用于固定渦動測量系統(tǒng)并可確保其在x、y、z 三個方向移動。

        試驗沉積物取自天津大學北洋園校區(qū)青年湖,為盡可能接近渦動相關法的適用環(huán)境,參照Donis等[24]的處理方法,試驗前將底泥中的植物去除后均勻平鋪于盛泥盒中,厚度約為5 cm,并用青年湖湖水進行7 天的靜置培養(yǎng)。

        3.2 試驗設計整理渦動相關法應用的流速條件[14-15,21,24],確定在0 ~10 cm/s 的流速范圍對沉積物-水界面處的溶解氧垂向分布及氧通量進行測量。試驗通過調(diào)節(jié)出水口處跌坎高度將水槽水深控制在30 cm;調(diào)節(jié)水泵閥門并讀取流量計讀數(shù)控制水槽入流流量進而計算得到相應試驗流速,在0.65 cm/s、0.96 cm/s、3.09 cm/s、5.25 cm/s、8.35 cm/s 和9.69 cm/s 共6 組水平流速條件下分別進行如下測量:

        (1)沉積物-水界面處溶解氧垂向分布測量:以5 min 為周期進行間歇采樣,每個周期包括3 min采樣時段和2 min 休眠時段,采樣頻率為64 Hz,采樣模式由ADV 控制,在休眠時段通過調(diào)節(jié)支架高度控制溶解氧傳感器的垂向位置。

        (2)沉積物-水界面處氧通量測量:采用連續(xù)測量模式,每次測量時長為30 min,采樣頻率64 Hz,參考Chotikarn 等[25]在室內(nèi)試驗中渦動相關系統(tǒng)的布置方式,本試驗測點位于沉積物上方5cm處,ARO-EC 探頭位于ADV 采樣體下游約2 cm 處,且與ADV 豎軸夾角為45°(圖2),每組水平流速平行測量4 次。

        4 結(jié)果與分析

        4.1 溶解氧在沉積物-水界面的垂向分布試驗選用的ARO-EC 型溶解氧傳感器在垂向移動過程中可能會破壞沉積物結(jié)構(gòu),為減小儀器在沉積物中的移動距離,有效降低儀器對沉積物結(jié)構(gòu)的破壞,因此采用與J?rgensen 等[26]類似做法測量溶解氧的垂向分布,進而分析不同水動力條件下DBL 厚度的變化規(guī)律。具體操作如下:

        (1)確定DBL 下邊界。從沉積物上方2 mm 處向下移動溶解氧傳感器,通過多次平行測量獲得溶解氧在沉積物-水界面較完整的垂向分布。圖3為水平流速8.44 cm/s 下溶解氧的垂向分布情況,可以看出,不同高度處的溶解氧濃度的變化情況存在明顯差異,在DBL 以上溶解氧沿垂向基本不變,當進入DBL 后溶解氧濃度隨高度減小基本呈線性減小,根據(jù)溶解氧剖面中線性分布的“拐點”得到DBL的下邊界,并記錄相應高度。

        圖3 SWI 溶解氧的垂向分布(水平流速8.44 cm/s)

        (2)測量溶解氧垂向分布。以DBL 下邊界為起始高度,向上移動溶解氧傳感器,測量不同水平流速下溶解氧在沉積物-水附近的垂向分布。圖4為溶解氧在沉積物-水界面附近不同水平流速下的垂向分布,可以看出,不同流速條件下沉積物-水界面處溶解氧的垂向分布存在明顯差異。

        圖4 不同水平流速下沉積物-水界面處溶解氧的垂向分布及DBL 厚度

        表1 不同水平流速下溶解氧垂向梯度、DBL 厚度和氧通量試驗結(jié)果

        (3)確定DBL 厚度。利用J?rgensen 等[26]提出的方法,將溶解氧線性濃度分布擬合線進行外延,外延線與溶解氧固定濃度的交點所對應的高度為DBL 厚度,相應結(jié)果列于表1??梢钥闯?,水平流速為0.65 ~9.69 cm/s 時,相應的溶解氧垂向梯度為12.18 ~59.88 mg/(L·mm),DBL 厚度為0.52 ~0.08 mm,DBL 厚度隨流速增大而減小。分析上述現(xiàn)象原因,隨著水平流速增加,水體紊動程度增強,DBL 上邊界處摻混隨之加劇,這使上覆水與DBL 內(nèi)溶解氧交換愈加充分,從而表現(xiàn)出DBL厚度相應減小。

        圖5 確定溶解氧和垂向流速滑動平均的窗口長度

        4.2 溶解氧在沉積物-水界面的通量以15 min 為一時段進行氧通量計算,首先將64 Hz 的原始測量數(shù)據(jù)降噪(8 Hz),再進行去尖峰、坐標旋轉(zhuǎn)校正、紊動值計算、時滯校正及通量計算,數(shù)據(jù)處理方法部分參照Kuwae 等[27],下面就紊動值計算和時滯校正進行說明。

        4.2.1 紊動值計算 湍流紊動值計算采用滑動平均法,該方法中滑動平均窗口長度的選取對氧通量的計算結(jié)果影響很大,當選取的窗口長度過小,會排除大尺度的紊動,造成通量低估;當選取的窗口長度過大,會引入非穩(wěn)定成分,使得氧通量出現(xiàn)波動。為確定合適的窗口長度,將初始窗口長度設為1 s,計算該窗口長度下氧通量,之后逐漸增加窗口長度,重復進行上述步驟,獲得不同窗口長度下的氧通量,當通量達到穩(wěn)定時所對應的窗口長度即為合適的窗口長度,計算過程如圖5所示??梢钥闯?,當窗口長度為100 s 時,不同水平流速下的氧通量均可保持穩(wěn)定,因此本次試驗選取100 s 作為最終窗口長度。

        4.2.2 時滯校正 ADV 與ARO-EC 空間位置分離及響應時間不同會造成溶解氧與垂向流速信號的不同步(時滯)。根據(jù)時滯產(chǎn)生原因,時滯的理論值應滿足下式:

        式中:x 為ARO-EC 尖端與ADV 采樣體的水平間距,2 cm;U 為水平流速,cm/s;tr為ARO-EC 的響應時間,0.5 s。

        由于氧通量結(jié)果對時滯取值較為敏感[28],而理論時滯與實際情況存在一定偏差[24],因此實際應用較少采用理論值。本文參照Lorrai 等[23]的方法確定實際時滯大小,首先參照相關研究中時滯大?。?9],確定時滯的取值為10 s,接著將溶解氧紊動值相對于垂向流速紊動值的時間序列進行移動,并計算兩者的相關性,時滯即為最大相關性所對應的移動時間。利用Matlab 軟件中的xcorr 函數(shù)進行時滯校正。圖6為水平流速5.25 cm/s 下垂向流速與溶解氧紊動值相關性隨移動時間的變化情況,可以看出,當移動時間為1.25 s 時,溶解氧和垂向流速紊動值的相關性最大,因此該水平流速下的時滯為1.25 s。采用相同方法,將不同水平流速下時滯的實測值與理論值進行分析(圖7),可以看出,時滯隨水平流速的增大而減小且與理論值擬合情況良好。說明理論時滯雖與實際情況存在一定偏差,但在實際應用中可以根據(jù)流速條件得到理論時滯值后,為時滯取值范圍的確定提供參考。

        圖6 垂向流速與溶解氧紊動值不同移動時間下的相關性(水平流速5.25 cm/s))

        圖7 時滯實測值與理論值對比

        圖8 累計氧通量

        圖9 氧通量在不同水平流速下的試驗結(jié)果

        為評估氧通量的數(shù)據(jù)質(zhì)量,對15 min 內(nèi)的氧通量進行累加得到累計氧通量(圖8)??梢钥闯?,累計氧通量有良好的線性趨勢,說明試驗過程的水動力條件穩(wěn)定,通量數(shù)據(jù)質(zhì)量良好。采用相同方法,得到不同水平流速下氧通量(圖9,表1)??梢钥闯觯搅魉贋?.65 ~9.69 cm/s 時,氧通量為-2.95±0.55 ~-25.12±2.64 mmol/(m2·d),氧通量均為負值,說明測量點處的的沉積物以耗氧為主,并且氧通量隨水平流速增加逐漸增大。

        4.3 氧通量、DBL 厚度和水動力條件的關系由不同水平流速下溶解氧在沉積物-水界面附近的垂向分布可以看出,DBL 厚度隨水平流速增大而減小,選取Batchelor 尺度對DBL 厚度進行參數(shù)化[30],Batchelor 尺度表達式為:

        圖10 DBL 厚度與Batchelor 尺度的關系

        為分析水動力條件對DBL 厚度的影響,將DBL 厚度與Batchelor 尺度進行線性擬合(圖10,藍色擬合線),可以看出,DBL 厚度隨Batchelor 尺度的增大線性增大,兩者呈顯著正相關關系,擬合關系式為δDBL1=0.96L*B1+0.04 ,相關系數(shù)為0.86。由于Batchelor 尺度可以表征湍擴散作用下標量波動所能保持的最小長度尺度,因此上述擬合關系式表明水動力條件對DBL 厚度影響顯著。考慮到DBL 厚度較難直接測量,而計算Batchelor 尺度的相關參數(shù)較容易獲取,因此對DBL 厚度進行參數(shù)化描述。整理相關試驗結(jié)果[30-32]并根據(jù)式(7)得到Batchelor 尺度。將相關試驗中DBL 厚度與Batchelor 尺度關系進行線性擬合(圖10,紅色擬合線),可以看出,本文和相關試驗擬合曲線較為一致,擬合關系式為δDBL2=0.89L*B2+0.08 ,相關系數(shù)為0.80。考慮到相關研究中未對影響運動黏滯系數(shù)和分子擴散系數(shù)的水溫條件進行實時觀測,因此Batchelor 尺度計算值與實際情況可能存在偏差。但整體而言,本文與相關試驗結(jié)果擬合公式中的參數(shù)非常接近,兩擬合公式均說明Batch?elor 尺度與DBL 厚度基本呈正相關關系(相關系數(shù)約為0.9),因此在實際應用中可用Batchelor 尺度近似表示DBL 厚度。

        為進一步研究DBL 厚度對氧通量的影響,將本文氧通量與DBL 厚度進行線性擬合(圖11,藍色擬合線)。可以看出,氧通量隨DBL 厚度的減小線性增大,兩者呈現(xiàn)負相關關系,擬合關系式為,相關系數(shù)為0.90,說明氧通量與DBL 厚度變化密切相關。整理相關研究[32-37]中氧通量與DBL 厚度并進行線性擬合(圖11,紅色擬合線),擬合關系式為相關系數(shù)為0.73。對比兩擬合曲線,發(fā)現(xiàn)本文試驗擬合曲線略高于相關試驗擬合曲線,這可能是由于室內(nèi)試驗過程中排除了生物活動的影響所致。但整體而言,兩條擬合曲線的變化趨勢較為一致:當DBL 厚度小于0.5 mm 左右時,DBL 厚度變化對氧通量影響更強烈,氧通量隨DBL 厚度減小迅速增大;DBL 大于0.5 mm 左右時,氧通量基本保持穩(wěn)定。本文試驗的DBL 厚度范圍為0.08 ~0.52 mm,最大DBL 厚度對應的水平流速為0.65 cm/s,結(jié)合Glud 等[13]的試驗結(jié)果,可以說明DBL 厚度大于0.5 mm左右時,氧通量基本保持不變。

        圖11 氧通量與DBL 厚度的關系

        4.4 討論通過將本文與相關研究結(jié)果對比,可以發(fā)現(xiàn)盡管不同研究中的沉積物的有機質(zhì)含量、孔隙度等條件存在差異,但水動力條件、擴散邊界層厚度、氧通量三者間的變化規(guī)律接近,即水動力條件對DBL 厚度影響顯著,Batchelor 尺度與DBL 厚度基本一致;DBL 厚度在小于0.5mm 時對氧通量影響劇烈,大于0.5 mm 時,氧通量基本保持穩(wěn)定。分析其原因一方面可能是由于物質(zhì)通過分子擴散、彌散及湍流擴散三種作用在SWI 進行傳輸,而在沉積物滲透雷諾數(shù)較高的環(huán)境下,最易受水動力條件影響的湍流擴散會成為物質(zhì)在SWI 垂向交換的主導方式[38]。另外,在其它因素較為穩(wěn)定的情況下,水動力條件會引發(fā)擴散邊界層厚度的改變,而SWI 物質(zhì)的垂向交換必須經(jīng)由擴散邊界層,使得其厚度變化較總氮、總磷、葉綠素含量等對SWI 氧通量的影響更為直接。因此表現(xiàn)出水動力條件會在某種程度上“掩蓋”其余因素對氧通量的影響,對于這一現(xiàn)象Murniati 等[7]也有類似闡述??紤]到多年來人們對于SWI 的物質(zhì)通量交換大多采用Ber?ner 等[39]于1980年建立的反應-輸運模型(Reaction-Transport Model,RTM)進行描述。在這一模型中擴散層厚度是關鍵影響參數(shù),因此本文得出的水動力條件可以通過改變擴散邊界層厚度進而影響氧通量的定量規(guī)律,也可適用于氨氮、硝酸鹽、可溶性有機碳等其它溶質(zhì)在SWI 的交換規(guī)律研究。

        本文對DBL 厚度進行參數(shù)化時依據(jù)相關研究中湍流耗散率的參考高度進行Batchelor 尺度的計算,關于Batchelor 尺度與DBL 厚度的內(nèi)在聯(lián)系尚待深入研究。此外,盡管不同環(huán)境條件下氧通量與DBL 厚度的變化趨勢一致,但是DBL 厚度對應的氧通量存在差異,這可能與水體及沉積物性質(zhì)有關,后期可加強氧通量與孔隙度、滲透性、有機質(zhì)含量、營養(yǎng)鹽濃度等靜態(tài)因素間響應關系的研究。

        5 結(jié)論

        本文圍繞沉積物氧通量與水動力條件間響應關系問題,較為系統(tǒng)地闡述了渦動相關法的原理及實現(xiàn)方法,基于此開展了不同水動力條件下擴散邊界層厚度及氧通量變化的試驗研究,主要結(jié)論如下:

        (1)隨著水平流速的增加,擴散邊界層厚度逐漸減小,氧通量逐漸增大。水平流速為0.65 ~9.69 cm/s時,擴散邊界層厚度為0.52 ~0.08 mm,氧通量為-2.95±0.55 ~-25.12±2.64 mmol/(m2·d)。

        (2)水動力條件對擴散邊界層厚度的影響明顯,擴散邊界層厚度與Batchelor 尺度呈正相關關系??紤]到擴散邊界層厚度較難直接測量,而計算Batchelor 尺度的相關參數(shù)較容易獲取,因此實際應用中可采用Batchelor 尺度近似表示擴散邊界層厚度。

        (3)當擴散邊界層厚度小于0.5 mm 左右時,擴散邊界層厚度變化對氧通量影響更強烈,氧通量隨擴散邊界層厚度減小迅速增大;當擴散邊界層厚度大于0.5 mm 左右時,氧通量基本保持穩(wěn)定。

        猜你喜歡
        邊界層溶解氧時滯
        帶有時滯項的復Ginzburg-Landau方程的拉回吸引子
        基于HIFiRE-2超燃發(fā)動機內(nèi)流道的激波邊界層干擾分析
        淺析水中溶解氧的測定
        污水活性污泥處理過程的溶解氧增益調(diào)度控制
        城市河道洲灘對水流溶解氧分布的影響
        一類具有邊界層性質(zhì)的二次奇攝動邊值問題
        一階非線性時滯微分方程正周期解的存在性
        非特征邊界的MHD方程的邊界層
        一類時滯Duffing微分方程同宿解的存在性
        鄭州市春季邊界層風氣候變化研究
        河南科技(2014年23期)2014-02-27 14:19:08
        精品国产a毛片久久久av| 青青青伊人色综合久久| 成美女黄网站18禁免费| 亚洲av第一区综合激情久久久 | 妇女自拍偷自拍亚洲精品| 国产三级精品三级在线专区2| 日韩色久悠悠婷婷综合| 久久开心婷婷综合中文| 国产av在线观看一区二区三区 | 国产一区二三区中文字幕| 国产无套一区二区三区久久| 国产自拍高清在线观看| 男女猛烈无遮挡免费视频| 国产肥熟女视频一区二区三区| 精品人妻潮喷久久久又裸又黄| 最新国产乱视频伦在线| 一本一道AⅤ无码中文字幕| 亚洲蜜桃视频在线观看| 亚洲啪啪色婷婷一区二区| 久久精品人搡人妻人少妇| 亚洲av成人无码一区二区三区在线观看 | 久久99精品久久久久久hb无码| а的天堂网最新版在线| 亚洲图文一区二区三区四区| 漂亮人妻出轨中文字幕| 亚洲一区二区女搞男| 香港台湾经典三级a视频| 国产av无码专区亚洲av琪琪| 国产a级午夜毛片| 69堂在线无码视频2020| 青青草视频免费在线播放| 在线观看国产成人自拍视频 | 亚洲中文字幕在线第二页 | 免费jjzz在线播放国产| 亚洲人成网站www| 久久精品国产亚洲av影院毛片| 丰满少妇人妻久久久久久| 久久久精品波多野结衣| 国产在线拍偷自拍偷精品| 久久夜色精品国产九色| 久久久久人妻精品一区二区三区 |