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

        ?

        矩形圓弧角養(yǎng)殖池系統(tǒng)固液兩相流數值模型的構建與進徑比影響研究

        2022-01-11 10:10:36薛博茹任效忠胡藝萱畢春偉
        關鍵詞:養(yǎng)殖池顆粒物流場

        薛博茹,任效忠,胡藝萱,畢春偉

        (1.設施漁業(yè)教育部重點實驗室(大連海洋大學),遼寧大連 116023;2.海岸和近海工程國家重點實驗室,大連理工大學,遼寧大連 116024)

        工廠化循環(huán)水養(yǎng)殖將傳統(tǒng)漁業(yè)生產經驗與現代工業(yè)化大生產模式相結合,依靠技術手段及現代管理模式已突破水產養(yǎng)殖的地域環(huán)境限制,提升生產自動化、信息化水平的同時兼顧魚類福利需求,并使養(yǎng)殖產品更符合公眾對水產品質量安全的要求[1-2]。工廠化循環(huán)水養(yǎng)殖符合提供優(yōu)質蛋白高效供給和拓展我國糧食安全的目標,與“藍色糧倉”國家漁業(yè)產業(yè)發(fā)展戰(zhàn)略相契合,是我國環(huán)境友好型水產養(yǎng)殖的重要轉型方向之一[3-4]。在工廠化循環(huán)水養(yǎng)殖模式中,大型化養(yǎng)殖池型結構備受現代水產養(yǎng)殖業(yè)的青睞,其可明顯降低單位魚產量的養(yǎng)殖投入成本與勞動力成本[5-6],有利于促進工廠化循環(huán)水養(yǎng)殖經濟效益的提升。但大型化養(yǎng)殖池對池型結構及進、出水系統(tǒng)的設計要求更高,不合理設計帶來的生產困擾更大。

        國內外學者已有研究表明[7-8]:養(yǎng)殖池內流場條件不佳會嚴重影響系統(tǒng)的水力混合性、增加死區(qū)容積,導致飼料分散,對沉淀顆粒物的匯集與沖洗等廢物管理也造成不利影響。此外,還會致使養(yǎng)殖水質下降,對養(yǎng)殖生物健康和疾病易感性、微生物生長、疾病擴散等亦有較大影響。循環(huán)系統(tǒng)進水管射流輸入池內的水體能量是循環(huán)水養(yǎng)殖的動力源與池內水體流場的直接驅動力,池內水動力學特征量是流場特性的評價指標,其中流速是最基本也是最重要最直觀的水動力參數。不同的進水結構將營造出不同的流場特性,養(yǎng)殖池系統(tǒng)的集排污性能也體現出與之相應的差別。因此,優(yōu)化養(yǎng)殖池系統(tǒng)的流場條件是改善集排污性能的重要前提,而池底的流速大小和流場分布規(guī)律直接關系到能否為沉降于池底的固體顆粒物提供有效驅動力,并最終影響?zhàn)B殖池系統(tǒng)的集排污性能。進水系統(tǒng)的結構參數設計直接影響池內流場條件,而養(yǎng)殖池池底的流場形態(tài)又與池內水動力條件密切相關,它們之間相互聯(lián)系又最終影響系統(tǒng)的集排污性能。薛博茹等[9-10]研究已發(fā)現,調整進徑比參數C/B 可顯著改善池內的流場條件,其中進徑比參數C/B 設置在0.02~0.04 區(qū)間,有利于池內均勻穩(wěn)定流場的形成,是矩形圓弧角養(yǎng)殖池的最佳進水管布設區(qū)間。

        因此,本研究從優(yōu)選進徑比參數區(qū)間改善養(yǎng)殖池系統(tǒng)內(尤其是池底)的流場條件為切入點探索進徑比參數對系統(tǒng)集排污的影響。采用RNG k-ε 湍流模型與DPM 模型相結合的方式建立固-液兩相流動數值計算模型,分析進徑比對沉降式固體顆粒物的匯聚規(guī)律和排出效率的影響,為改善循環(huán)水養(yǎng)殖系統(tǒng)的集排污性能提供理論依據。

        1 材料與方法

        計算流體動力學(CFD)基于質量守恒方程、動量守恒方程和能量守恒方程等流體動力學三個基本方程,建立數值模型進行模擬仿真計算[11]。本研究針對養(yǎng)殖池系統(tǒng)內流體動力特性,假定流體有粘性但不可壓縮,不涉及能量方程求解。連續(xù)性方程和動量方程表示如下:

        連續(xù)性方程:

        動量方程:

        式中:ρ 為流體密度,kg·m-3;U 為笛卡爾坐標系中的矢量速度場,m·s-1;ui、uj為速度分量(i,j=1,2,3);p 為壓力,N·m-2;τ 為剪應力,N·m-2;g 為重力加速度,m·s-2;F 為外力,N;Δ為笛卡爾坐標系中的梯度或算子。

        1.1 湍流模型

        流體流動模型選用RNG k-ε 湍流模型,RNG k-ε 湍流模型基于湍流動能k 方程和湍流耗散率ε 方程,通過修正湍動粘度考慮了平均流動中的旋轉及旋轉流動情況。與采用標準的k-ε 湍流模型相比,RNG k-ε 湍流模型在處理應變率高以及流線彎曲程度較大的流動方面有優(yōu)勢,因而它可更精準的模擬以旋轉和沉降為主的養(yǎng)殖池內流場的流動[12-13]。

        采用RNG k-ε 湍流模型建立流體數值模型,湍流動能k 方程和湍流耗散率ε 方程的表達式如下:

        湍流動能k 方程:

        湍流耗散率ε 方程:

        式中:μ 為流體動力粘度,N·s·m-2;μt為湍流粘度,N·s·m-2;αk=αε=1.39 分別為k 和ε 的反向有效普朗特數;Gk為平均速度引起的湍動能k 的產生項。

        式中:C1ε=1.42,C2ε=1.68,Cμ=0.084 5(C1ε、C2ε、Cμ均為根據經驗所得常數)。

        1.2 離散相(DPM)模型

        本文的固-液兩相流問題,主要研究養(yǎng)殖池系統(tǒng)內水動力條件對沉降式固體顆粒物的影響。本文建立固-液兩相流數值模型,將流體視作連續(xù)介質,固體顆粒相作為離散相來處理。養(yǎng)殖池系統(tǒng)水體是主體,固體顆粒物占比總體較小,顆粒物的體積分數遠小于10%。且固體顆粒沿自身軌跡運動,適合采用拉格朗日法進行建模計算。因此,本研究選用DPM 模型進行固體顆粒物的數值計算分析,用以獲取更精準的顆粒物運動軌跡與統(tǒng)計排出效率[14-15]。

        FLUENT 軟件使用積分拉式坐標系下的微分方程對離散相顆粒軌跡進行求解。固體顆粒受力的微分方程為(以X 方向為例)[15]如下:

        式中:u 為連續(xù)相的速度,m·s-1;up為顆粒速度,m·s-1;ρ 為流體密度,kg·m-3;ρp為顆粒密度,kg·m-3;FD(u-up)為單位質量顆粒受到的曳力,N;gx(ρp-ρ )/ρp為單位質量粒子受到的重力,N;FX表示其他各個作用力在X 方向的分力,N。

        式中:dp為顆粒直徑,m;μ 為流體動力粘度,N·s·m-2;CD為曳力系數;Re 為相對雷諾數(顆粒雷諾數)。

        2 矩形圓弧角養(yǎng)殖池數值模型

        2.1 幾何模型與網格劃分

        本文采用Geometry 軟件構建矩形圓弧角養(yǎng)殖池三維幾何模型,導入Mesh 軟件并采用多區(qū)域網格劃分(Muitizone)方式進行計算區(qū)域網格劃分。養(yǎng)殖池模型結構參數設置如圖1 所示,其中長寬比L/B=1:1,相對弧寬比R/B=0.25,徑深比B/H=5:1。模型主體結構尺度參照工廠生產實踐中的養(yǎng)殖池進行約為1:8 的比尺縮放,本文對部分參數進行優(yōu)化與簡化處理,養(yǎng)殖池高與水深設為同一高度。進水管置于養(yǎng)殖池直壁中間位置(依據進徑比參數設置,調整其中心與池側壁的水平距離),進水管與養(yǎng)殖池同高設為200 mm、直徑為20 mm,出水口置于池底中心位置(直徑Dout=0.02 m)。計算工況流量恒為0.470 m3·h-1,射流速度V=1 m·s-1,進水管結構在水深范圍內均勻開12 孔。

        圖1 矩形圓弧角養(yǎng)殖池模型圖(mm)Fig.1 The schematic diagram of the rectangular arc angle model tank(mm)

        Muitizone 多區(qū)網格劃分可依據設定自行判斷網格劃分區(qū)域,主體區(qū)域生成純六面體網格單元,而幾何體較復雜區(qū)域采用更好適應邊界的非結構化網格單元。本文采用固-液兩相流模型進行數值模擬,選擇Muitizone 對養(yǎng)殖池內流場進行網格劃分兼顧提高計算效率與節(jié)約計算資源的需求,系統(tǒng)網格劃分如圖2(A)。本研究重點分析矩形圓弧角養(yǎng)殖池內的流場特性對系統(tǒng)集排污機制的影響,通過對比相同工況條件下同一監(jiān)測時間點固體顆粒物累計排出率來評估網格尺度對排污計算精度的影響,以獲取用于數值計算的合理網格數量。在3 種網格單元數下,系統(tǒng)在各個監(jiān)測時間點的固體顆粒物累計排出率如圖3-7(B),其中,mesh-C 網格與其他兩種網格存在較大差異,固體顆粒物累計排出率的模擬預測值低于mesh-A 和mesh-B 網格的模擬預測結果,而mesh-A 和mesh-B 網格模擬趨勢與數值呈現出高度一致性。此外,在同樣計算精度下,mesh-B 網格單元數小于mesh-A 網格單元數,因而mesh-B 網格即網格數量約為560 000可用于矩形圓弧角養(yǎng)殖池內固體顆粒物匯集規(guī)律與累計排出效率的計算分析,同時保持了相對合理可接受的計算時間要求。

        圖2 網格劃分與網格無關性驗證(mesh-A,B 和C 分別為248 141,5 59 214,1 080 455)Fig.2 Mesh generation and mesh independence verification(mesh-A,B and C is 248 141,5 59 214,and 1 080 455 respectively)

        2.2 數值求解方法

        研究基于CFD 仿真軟件FLUENT 16.0 平臺和戴爾服務器開展數值模擬計算,FLUENT 中離散相模型(DPM)將流體視為連續(xù)相、固體顆粒視為離散相。本研究對液相模型求解選擇穩(wěn)態(tài)(Steady)控制方式;向池內添加固體顆粒物后即固-液兩相流模式,對連續(xù)相模型(液相模型)求解仍然選取穩(wěn)態(tài)(Steady)控制方式,針對固相模型考慮固-液雙向耦合(Interaction with Continuous Phase),并開啟非定常粒子追蹤(unsteady particle tracking)方式,以利于真實模擬固體顆粒物沉降規(guī)律和運行軌跡。連續(xù)相和離散相均采用壓力隱式求解方法,壓力速度耦合方式選用SIMPLE 算法,壓力、動量基于二階迎風離散格式(second order upwind)求解,湍流動能、湍動能耗散率基于一階迎風離散格式(first order upwind)求解。程序計算參數設置中,設定壓力的亞松弛因子采用0.3,動量的亞松弛因子采用0.7,湍動能和湍流耗散率的亞松弛因子為0.8,其他項保持程序常規(guī)設置。

        2.3 邊界條件

        本研究設置初始邊界條件包括進、出口邊界,交界面處理和壁面邊界等,如表1 所示。同時對數值模型計算假設如下:①不考慮溫度的影響;②內不考慮養(yǎng)殖生物,外不考慮其他因素影響;③假設流體有粘性但不可壓縮;④不考慮顆粒物運動的質量損失,即假設顆粒間不存在碰撞的質量損失,不存在顆粒物的分解現象;⑤考慮了薩夫曼升力和壓力梯度力。

        表1 初始邊界條件設定Tab.1 Initial boundary conditions

        3 結果

        3.1 模型驗證

        本研究通過建立與數值模型主體尺寸相同的物理試驗模型如圖4(A)所示,驗證所構建固-液兩相流模型的合理性與有效性。物理試驗模型進水結構為對弧雙管,每管均勻開9 孔合計雙管共18 孔,流量設為恒定Q=0.814 m3·h-1。試驗以常興鲆鰈魚類2#沉性飼料代替固體污染物顆粒,常興鲆鰈魚類2#沉性飼料為類圓柱形固體小顆粒,其粒徑為2.5±0.2 mm,粒徑厚度為1:1.0±0.1 mm;飼料撒入水中可快速沉降池底,且實驗過程無松散分解現象。物理模型試驗中待養(yǎng)殖池內流場循環(huán)穩(wěn)定運行后(約20 min),向池內均勻撒入20 g 飼料并開始計時,每10 s 拍攝1 張照片監(jiān)測固體顆粒物分布,拍攝持續(xù)3 min 并記錄3 min 內各時刻固體顆粒物累計排出情況。典型試驗拍攝照片如圖3(A)所示,使用Matlab 軟件進行照片后處理及圖像分析,識別出拍攝圖片中的顆粒物總像素面積并依據公式(7)計算出各時刻固體顆粒物的累計排出率,且每個工況進行3 次重復試驗并取平均值。數值模擬通過檢測池內平均流速和殘差曲線,在池內流場達到穩(wěn)定運行后向養(yǎng)殖池系統(tǒng)注入與物理模型試驗同參數同數量的固體顆粒物。數值模擬過程中可直接監(jiān)測出水口處顆粒物排出個數,通過顆粒物排出個數與注入個數之比計算出固體顆粒物的累計排出率。數值模擬中固體顆粒物的屬性參照物理試驗中常興鲆鰈魚類2#沉性飼料設置。顆粒注入方式:在養(yǎng)殖池水面(hz=0.2 m)以池心為圓心,以0.45 m 為半徑的圓形區(qū)域內均勻撒布固體顆粒(如圖3(B)),固體顆粒具體參數設置如表1(經換算與計數統(tǒng)計,20 g 常興鲆鰈魚類2#沉性飼料約有900 個顆粒,數值計算注入900 個固體顆粒)。

        圖3 固體顆粒物撒布圖Fig.3 Distribution chart of solid particles

        式中:S 為監(jiān)測20 g 顆粒物撒布后總像素面積,m2;St為監(jiān)測各時刻池內剩余顆粒物像素面積,m2。

        圖4(B)給出數值模擬結果與物理模型試驗結果對比圖,圖中顯示數值模擬結果與試驗結果的固體顆粒物累計排出率趨勢吻合良好,且各個監(jiān)測時刻的累計排出率的誤差均在4%以內。模型驗證結果充分證明本文建立的固-液兩相流數值模型充分反映了系統(tǒng)流場內固體顆粒物的運行特性與排污性能,建立的數值模型合理、計算精度高,可應用于模擬系統(tǒng)內固體顆粒物的沉降規(guī)律與排出效率的分析研究,為分析系統(tǒng)的流場特性與集排污性能的關系和改進系統(tǒng)設計提供了一個有效的解決方法。

        圖4 物理試驗模型以及數值模擬結果與試驗結果比較圖Fig.4 Physical test model and comparison of numerical simulation results with experimental results

        3.2 各進徑比參數下顆粒累計排出率隨時間變化

        薛博茹等[9]已通過調整進水管布設位置與養(yǎng)殖池側壁之間的水平距離,優(yōu)選出進徑比參數C/B 位于0.02~0.04 區(qū)間有利于矩形圓弧角養(yǎng)殖池系統(tǒng)獲得最佳水動力條件,實現了有效改善池內與池底流場形態(tài)的目標。而進徑比參數C/B<0.02 或C/B>0.05 均不利于養(yǎng)殖池底部二次流場的形成,且會破壞均勻穩(wěn)定的流場形態(tài),進而給養(yǎng)殖池系統(tǒng)的集排污性能帶來不利影響。尤其是當進徑比過大(即C/B>0.13),進水系統(tǒng)射出的高速水體,受到射流口和中心出水口的共同影響,趨于做向心繞流運動,帶動養(yǎng)殖池中間區(qū)域流體高速運轉而忽略養(yǎng)殖池邊壁尤其是4 個圓弧角區(qū)域。現依據進徑比對單通道矩形圓弧角養(yǎng)殖池系統(tǒng)內流場特性的影響,選取較有代表性的4 種進徑比參數(C/B=0.01,0.03,0.09 和0.14),分析進徑比對沉降式固體顆粒物排出效率的影響。

        圖5(A)給出4 種進徑比參數工況下,沉降式固體顆粒物前180 s 的累計排出率,可見150 s 后池內固體顆粒物累計排出率趨于穩(wěn)定,進徑比參數對養(yǎng)殖池內固體顆粒物的排出效率影響顯著,其中進徑比參數C/B=0.03 時,固體顆粒物的累計排出率最高。為進一步對比說明進徑比參數對養(yǎng)殖池內流速和排污性能的影響,以進徑比參數C/B=0.01 時的池內平均流速和固體顆粒物累計排出率(180 s)為基準,對4 種進徑比參數下的池內平均流速和累計排出率(180 s)進行歸一化處理,歸一化結果如圖5(B),可見養(yǎng)殖池內的固體顆粒物排出效率與養(yǎng)殖池系統(tǒng)的平均流速密切相關。其中,Vavg-0.01定義為進徑比參數C/B=0.01 時的池內平均流速,Vavg-a定義為其他某一進徑比參數下養(yǎng)殖池內的平均流速;Vavg-0.01定義為第180 s 時刻進徑比參數C/B=0.01 時的固體顆粒物累計排出率,Vavg-a定義為第180 s 時刻其他某一進徑比參數下固體顆粒物累計排出率。進徑比參數C/B=0.03 時的池內平均流速比進徑比參數C/B=0.01 的池內平均流速提高約15%,而第180 s 時刻固體顆粒物的累計排出率提高約52%。流速是影響?zhàn)B殖系統(tǒng)內水動力特性最基本最重要的因素,沉降于養(yǎng)殖池底部的固體顆粒物能否順利排出,與養(yǎng)殖池底部流場特性尤其是底層流場流速是否可達到顆粒物的啟動速度密切相關,故基于調整進徑比參數優(yōu)化水動力條件可顯著提高養(yǎng)殖池系統(tǒng)的集排污性能。

        此外,進徑比參數C/B=0.14 的池內平均流速雖然高于進徑比參數C/B=0.03 工況,但固體顆粒物的累計排出率卻呈現出相反趨勢;第180 s 時刻,C/B=0.14 工況的固體顆粒物的累計排出率比進徑比參數C/B=0.03 工況的固體顆粒物的累計排出率降低60%。這一現象充分說明養(yǎng)殖池系統(tǒng)的集排污性能除了受到系統(tǒng)內流場整體的平均速度影響外,還與系統(tǒng)內流場的均勻性密切相關,只有養(yǎng)殖池系統(tǒng)底層流場流速分布均勻,且均達到固體顆粒物的最低啟動速度,養(yǎng)殖池系統(tǒng)才能整體上實現較好的集排污性能。

        4 討論

        圖6 給出固體顆粒物分布隨時間變化圖,由圖可見進徑比參數C/B 由0.01 調整到0.03,與養(yǎng)殖池系統(tǒng)集排污性能密切相關的底層流場的形態(tài)明顯改善;與此相對應,沉降于池底的固體顆粒物向池心排污口匯聚現象明顯,固體顆粒物的排出效率明顯提高。當進徑比參數C/B 增加到0.09 時,由于高速水體與池側壁發(fā)生直接碰撞現象逐漸劇烈導致能量消耗較高,使得用于維持養(yǎng)殖池系統(tǒng)內水體運轉的能量降低,直觀表現為系統(tǒng)平均流速降低且易出現低流速區(qū)和小漩渦區(qū),因此固體顆粒物仍然呈現分散分布狀態(tài),甚至在低流速區(qū)和小漩渦區(qū)域徘徊或啟動動力不足而無法向池心排水口匯聚,因此固體顆粒物的排出效率明顯下降。而當進徑比參數增加到0.14 時,由射流口射出的高流速水體帶動出水口附近流體高速運轉,此時養(yǎng)殖池內的固體顆粒物一部分散落在臨近池壁的低流速區(qū)域無法向池心匯聚,另一部分雖然快速匯聚到出水口附近外圍區(qū)域,卻隨流體高速旋轉,均無法順利排出;該工況養(yǎng)殖池系統(tǒng)內平均流速雖然較高但池底整體速度分布不均勻,在中心區(qū)域出現較大區(qū)域高速旋轉渦流區(qū),反而使固體顆粒物排出效率降低,無法滿足養(yǎng)殖系統(tǒng)的集排污需求。同時,久泡松散的殘餌糞便因池中心高速運轉流體中剪切力的增大,極易把殘餌糞便等固體廢棄物打散分解,使懸浮固體顆粒物的濃度大幅增加[16],進而影響?zhàn)B殖池系統(tǒng)的穩(wěn)定運轉與惡化系統(tǒng)的養(yǎng)殖水質環(huán)境。因此,除較小流場流速和較低的流場均勻性外[17-18],池中心排污口附近的高速旋轉渦流讓粒子在池中心形成高速環(huán)形運動軌跡而無法順利進入排水口,也是影響系統(tǒng)集排污性能的另一大不利因素。

        5 結論

        本研究運用計算流體動力學仿真技術,基于RNG k-ε 湍流模型與離散相DPM 模型相結合構建固-液兩相流數值計算模型,驗證其合理性和滿足計算精度的基礎上,開展了養(yǎng)殖池系統(tǒng)內沉降式固體顆粒物的集排污性能影響研究,分析了系統(tǒng)流場特性對固體顆粒物匯聚與排出的機理。結論如下:養(yǎng)殖池系統(tǒng)的集排污性能與養(yǎng)殖池內尤其是池底的流場條件密切相關;養(yǎng)殖池系統(tǒng)的集排污性能除了與流場速度有關外,還與系統(tǒng)的流場均勻性和池中心的高速渦流區(qū)強弱密切相關。進徑比對養(yǎng)殖池系統(tǒng)內固體顆粒物的分布規(guī)律與排出效率影響顯著,基于調整進徑比可改善池內流場形態(tài),尤其是與集排污性能有直接關系的池底流場形態(tài),為沉降式固體顆粒物的排出奠定基礎。進徑比參數C/B=0.03,既有利于單通道排污矩形圓弧角養(yǎng)殖池系統(tǒng)獲得最佳流場條件,又可顯著提高養(yǎng)殖池系統(tǒng)的集排污性能。此外,基于流場特性分析,進徑比參數設置在0.02~0.04 區(qū)間,系統(tǒng)也會在獲得最佳流場條件的同時兼具較好的集排污性能。

        猜你喜歡
        養(yǎng)殖池顆粒物流場
        水體消毒防病用藥稱量分配的簡捷方法
        相對弧寬比對雙通道方形養(yǎng)殖池的流場優(yōu)化研究
        持續(xù)陰雨期間鱖魚養(yǎng)殖池常見問題及管理措施
        當代水產(2021年6期)2021-08-13 08:17:34
        1989—2019年環(huán)渤海地區(qū)養(yǎng)殖池的空間變化特征
        海洋通報(2021年1期)2021-07-23 01:55:38
        大型空冷汽輪發(fā)電機轉子三維流場計算
        大電機技術(2021年2期)2021-07-21 07:28:24
        轉杯紡排雜區(qū)流場與排雜性能
        南平市細顆粒物潛在來源分析
        基于HYCOM的斯里蘭卡南部海域溫、鹽、流場統(tǒng)計分析
        錯流旋轉填料床脫除細顆粒物研究
        化工進展(2015年3期)2015-11-11 09:18:15
        基于瞬態(tài)流場計算的滑動軸承靜平衡位置求解
        国产一区曰韩二区欧美三区| 久久国产成人精品av| 婷婷色香五月综合激激情| 特级毛片a级毛片免费播放| 国产一区二区精品久久凹凸| 久草视频在线播放免费| 久久精品亚洲熟女av蜜謦| 亚洲精品无人区| 波多野结衣视频网址| 亚洲成av人片在线天堂无| 国产一区二区三区视频地址| 亚洲字幕av一区二区三区四区| 精品久久久久久国产| 一区二区三区在线视频免费观看| 尤物国产一区二区三区在线观看| 成人午夜福利视频后入| 精品推荐国产精品店| 亚洲国产综合精品久久av| 国产女同舌吻1区2区| 无码丰满熟妇一区二区| 天天狠狠综合精品视频一二三区| 日韩精品中文字幕免费人妻| 91在线视频在线视频| 一本一道av中文字幕无码| 日本手机在线| 亚洲国产一区二区视频| 久久久久88色偷偷| 无遮高潮国产免费观看| 亚洲女同系列高清在线观看| 一区二区三区在线少妇| 久久久g0g0午夜无码精品| 无码片久久久天堂中文字幕| 日本熟妇裸体视频在线| 蜜臀av在线观看| 欧美成人看片黄a免费看| 亚洲一二三四五区中文字幕| 一区二区黄色在线观看| 欧美最大胆的西西人体44| 久久99久久久无码国产精品色戒 | 久久亚洲国产精品成人av秋霞| 亚洲成人av一区二区麻豆蜜桃|