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

        ?

        大型氣候?qū)嶒?yàn)室氣流組織仿真分析

        2019-12-25 13:10:18馬建軍姜亞軍
        裝備環(huán)境工程 2019年11期
        關(guān)鍵詞:實(shí)驗(yàn)室模型

        馬建軍,姜亞軍

        (中國(guó)飛機(jī)強(qiáng)度研究所,西安 710065)

        大型氣候?qū)嶒?yàn)室可模擬極端低溫、高溫、濕熱、降雪、凍雨、太陽(yáng)輻照、淋雨等綜合氣候環(huán)境,滿足大型裝備的綜合氣候環(huán)境試驗(yàn)需求[1-2]。在這些氣候環(huán)境中,溫度是最基本的因素之一,試驗(yàn)溫度允差應(yīng)在±2 ℃范圍內(nèi)(對(duì)于大于5 m3的試件,溫度允差可以為±3 ℃)[3]。此外,不同氣候試驗(yàn)項(xiàng)目對(duì)風(fēng)速還有要求[4],如GJB 150A 中的高溫、低溫試驗(yàn)規(guī)定試件附近風(fēng)速不超過1.7 m/s,太陽(yáng)輻照試驗(yàn)風(fēng)速應(yīng)盡可能小,通常保持在0.25~1.5 m/s。實(shí)驗(yàn)室的溫度和速度場(chǎng)是氣流組織決定的。由于實(shí)驗(yàn)室尺寸較大,層高超過22 m,有效試驗(yàn)容積近10 000 m3,屬于高大空間[5],需要有效的設(shè)計(jì)氣流組織,以滿足試驗(yàn)溫度均勻性和風(fēng)速要求。

        氣流組織設(shè)計(jì)方法主要有四種:射流公式法、Zonal Model 模型、模型實(shí)驗(yàn)和CFD 數(shù)值模擬。前三種方法在適應(yīng)性、計(jì)算準(zhǔn)確度、投資、周期等方面存在問題。自1974 年丹麥P. V. Nielsen 首次將CFD 技術(shù)應(yīng)用于室內(nèi)氣流CFD 數(shù)值模擬以來,CFD 在氣流組織設(shè)計(jì)上越來越成熟[6],并且已應(yīng)用于小型實(shí)驗(yàn)箱的氣流組織設(shè)計(jì)研究[7-8]。文中采用CFD 方法對(duì)大型氣候?qū)嶒?yàn)室氣流組織進(jìn)行數(shù)值分析,以獲得合理的氣流組織設(shè)計(jì)方案。

        1 氣流組織形式

        一般環(huán)境實(shí)驗(yàn)室氣流組織形式主要有三種:側(cè)送側(cè)回、上送下回、下送上回。小型環(huán)境實(shí)驗(yàn)箱/實(shí)驗(yàn)室多采用側(cè)送側(cè)回和上送下回(孔板送風(fēng))的形式。側(cè)送側(cè)回送風(fēng)射流射程比較長(zhǎng),送風(fēng)氣流在到達(dá)工作區(qū)之前已充分混合,速度場(chǎng)和溫度場(chǎng)都趨于均勻和穩(wěn)定。對(duì)于大型氣候?qū)嶒?yàn)室,送風(fēng)速度很高。上送下回氣流組織常見的是全面孔板送風(fēng)和密布散流器送風(fēng),全面孔板送風(fēng)射流的擴(kuò)散和混合性較好,射流的混合過程很短,溫差和風(fēng)速度衰減快,不適用于高大空間實(shí)驗(yàn)室;散流器送風(fēng)則可根據(jù)降溫或加熱需求,改變送風(fēng)射流的方向,降溫時(shí)水平送風(fēng),加熱時(shí)豎直向下送風(fēng),以抵消浮力作用。國(guó)外大型氣候?qū)嶒?yàn)室如美國(guó)麥金利氣候?qū)嶒?yàn)室和韓國(guó)ADD 氣候?qū)嶒?yàn)室均采用上送下回的散流器送風(fēng)氣流組織形式,如圖1 所示。

        圖1 國(guó)外大型氣候環(huán)境實(shí)驗(yàn)

        美國(guó)麥金利氣候?qū)嶒?yàn)室寬76 m、深61 m,中心高為21 m,下沿高度為11 m,容積為93 000 m3,實(shí)驗(yàn)室頂部有24 個(gè)散流器送風(fēng),回風(fēng)口在實(shí)驗(yàn)室右側(cè)壁面下部。韓國(guó)ADD 實(shí)驗(yàn)室寬42 m,深32 m,高15 m,頂部有9 個(gè)散流器送風(fēng),回風(fēng)口在實(shí)驗(yàn)室后側(cè)壁面下部。

        大型氣候?qū)嶒?yàn)室尺寸與麥金利實(shí)驗(yàn)室相當(dāng),但其為平頂結(jié)構(gòu),有效高度超過麥金利實(shí)驗(yàn)室,如圖2 所示。氣流組織形式同樣采用上送下回的形式,采用30 個(gè)送風(fēng)口均勻分布,右側(cè)壁面下部均勻分布5 個(gè)回風(fēng)口。由于散流器送風(fēng)射流角度通常只有兩種形式(水平或豎直),旋流器卻可以通過調(diào)節(jié)動(dòng)葉片來連續(xù)調(diào)整送風(fēng)射流角度,且摻混效果比散流器好[9],因此采用旋流器作為送風(fēng)口。

        圖2 大型氣候?qū)嶒?yàn)室

        實(shí)驗(yàn)室空艙條件下的速度場(chǎng)和溫度場(chǎng)分布是最基本的氣流組織設(shè)計(jì)內(nèi)容,若氣流組織設(shè)計(jì)滿足極端低溫和極端高溫兩個(gè)極端溫度工況,可以認(rèn)為氣流組織也將滿足中間溫度工況。文中主要分析極端低溫和極端高溫兩種空艙狀態(tài)下實(shí)驗(yàn)室的氣流組織。

        2 模型建立

        2.1 數(shù)值方法

        室內(nèi)的空氣流動(dòng)和溫度分布要滿足連續(xù)性方程、動(dòng)量方程和能量方程,這些方程可以表示成式(1)所示的通用形式。

        式中:ρ 為空氣密度;t 為時(shí)間;對(duì)于連續(xù)性方程φ= 1,對(duì)于動(dòng)量方程 φ = Vj(j=1,2,3),對(duì)于能量方程φ=T ;v 是速度矢量; Γφ為擴(kuò)散系數(shù);Sφ為源項(xiàng)。

        采用FLUENTL 軟件進(jìn)行分析,湍流模型采用常用的k-e 兩方程模型[10-11],近壁面采用標(biāo)準(zhǔn)壁面函數(shù),考慮浮力的影響開啟重力選項(xiàng)。為簡(jiǎn)化計(jì)算,忽略對(duì)結(jié)果影響不大的因素,對(duì)計(jì)算模型進(jìn)行簡(jiǎn)化如下:忽略壁面、大門、頂棚的接縫和泄漏,認(rèn)為其是大平板結(jié)構(gòu);空氣采用不可壓縮理想氣體模型,僅考慮溫度對(duì)密度的影響;空氣的黏性采用薩特蘭公式;由于壁面溫度相差很小,忽略壁面間的輻射換熱。

        大空間氣流組織計(jì)算由于重力的作用,收斂較為困難,為達(dá)到快速收斂的目的,采用如下的計(jì)算策略:壓力差分采用Body Force Weighted,其余均采用二階迎風(fēng)差分格式;首先進(jìn)行穩(wěn)態(tài)計(jì)算,不開啟重力選項(xiàng),直到達(dá)到穩(wěn)定狀態(tài),獲得初步穩(wěn)定流場(chǎng)與溫度場(chǎng);然后進(jìn)行瞬態(tài)計(jì)算,開啟重力選項(xiàng),步長(zhǎng)為5 s,逐漸逼近,直到達(dá)到穩(wěn)定狀態(tài),完成計(jì)算。

        2.2 旋流器簡(jiǎn)化模型

        典型的旋流器結(jié)構(gòu)如圖3 所示,主要包括靜葉片、動(dòng)葉片(深色部分)和外殼。通過調(diào)節(jié)動(dòng)葉片,改變其與靜葉片之間的夾角α,即可實(shí)現(xiàn)送風(fēng)射流角度從水平到豎直的連續(xù)變化。

        圖3 旋流器結(jié)構(gòu)

        數(shù)值分析時(shí),如果對(duì)風(fēng)口進(jìn)行詳細(xì)幾何建模,由于風(fēng)口的尺寸與計(jì)算空間相比很小,而風(fēng)口需要較細(xì)的網(wǎng)格才能準(zhǔn)確模擬流動(dòng)狀態(tài),將導(dǎo)致網(wǎng)格繁多,計(jì)算周期大大延長(zhǎng)。通常采用風(fēng)口簡(jiǎn)化模型來模擬風(fēng)口的流動(dòng)[12-13],對(duì)于旋流器,采用指定速度法進(jìn)行模擬,如圖4 所示。旋流風(fēng)口的葉片出風(fēng)平面簡(jiǎn)化為一個(gè)圓形面,將其等分為12 個(gè)扇形區(qū)域,其中6 塊陰影部分為出風(fēng)面,其余為封閉區(qū)域。送風(fēng)特性由三個(gè)參數(shù)決定:外徑R2、內(nèi)徑R1以及切向角θ,送風(fēng)軸向速度vaxis按式(2)計(jì)算。

        式中:V˙為風(fēng)口體積流量,m3/s。送風(fēng)切向速度vtan按式(3)計(jì)算。

        圖4 旋流風(fēng)口簡(jiǎn)化模型

        根據(jù)坐標(biāo)變換得到直角坐標(biāo)系下的速度分量,再通過UDF 附加到出風(fēng)面當(dāng)作邊界條件。為驗(yàn)證簡(jiǎn)化模型的可行性,對(duì)詳細(xì)模型和簡(jiǎn)化模型進(jìn)行了數(shù)值分析對(duì)比。旋流器尺寸為B=1.5 m,R=1.9 m,R1=0.5 m,R2=1.5 m,T=0.2 m。如圖5 所示,將旋流器安裝在一個(gè)直徑40 m、高度4 0m 的空間內(nèi),風(fēng)口懸掛在空間頂部圓心處,出風(fēng)口平面與空間頂部3.5 m,計(jì)算網(wǎng)格劃分如圖5 所示。

        圖5 旋流風(fēng)口數(shù)值分析計(jì)算網(wǎng)格

        詳細(xì)模型和簡(jiǎn)化模型的送風(fēng)射流軌跡如圖6 和圖7 所示??梢姦?與α 存在對(duì)應(yīng)關(guān)系,通過改變?chǔ)?的值可以近似模擬實(shí)際的流動(dòng)狀態(tài),可以模擬α 的范圍為5°~25°。

        圖6 真實(shí)模型流動(dòng)軌跡

        簡(jiǎn)化模型與詳細(xì)模型的流動(dòng)分布對(duì)比見表1,其中流量為8 m3/s。表1 中Modelr和Models分別指真實(shí)模型和簡(jiǎn)化模型,vsmax為風(fēng)口出風(fēng)面最大風(fēng)速,L1.7、L1.0分別表示送風(fēng)射流末端速度為1.7、1.0 m/s時(shí)的射程。從表1 中可以看出,可以根據(jù)實(shí)際風(fēng)口葉片夾角α 的流動(dòng)狀態(tài)來確定簡(jiǎn)化模型的切向角θ。最大面風(fēng)速因?yàn)楹?jiǎn)化模型的先天缺陷而與詳細(xì)模型存在一定差異,但在射程上則差別很小,可以準(zhǔn)確復(fù)現(xiàn)送風(fēng)氣流在室內(nèi)的流動(dòng),簡(jiǎn)化模型是可行的。下面將采用這種簡(jiǎn)化模型對(duì)大空間氣候?qū)嶒?yàn)室的氣流組織進(jìn)行分析與優(yōu)化。

        圖7 簡(jiǎn)化模型流動(dòng)軌跡

        2.3 計(jì)算網(wǎng)格和邊界條件

        由于送風(fēng)口出風(fēng)流動(dòng)復(fù)雜,為準(zhǔn)確反應(yīng)流動(dòng)同時(shí)降低計(jì)算量,在風(fēng)口附近、地坪、天花板、大門、保溫壁面附近網(wǎng)格同樣進(jìn)行加密處理。第一層網(wǎng)格高度為10 mm,并在計(jì)算過程中動(dòng)態(tài)調(diào)整網(wǎng)格尺寸,以滿足壁面30<y+<300,避免不合理的傳熱[14-15]。最終網(wǎng)格數(shù)量為200 萬,網(wǎng)格劃分如圖8 所示。

        邊界條件的設(shè)置如下所述。

        1)旋流器:速度分布按2.2 節(jié)中的簡(jiǎn)化模型通過UDF 施加,送風(fēng)溫度極端低溫時(shí)為-60 ℃,極端高溫時(shí)為+79 ℃。

        2)回風(fēng)口:壓力出口條件。

        3)地板:地板為混凝土結(jié)構(gòu),蓄熱能力較強(qiáng),且與氣流組織耦合傳熱,很難預(yù)測(cè)其達(dá)到溫度穩(wěn)定的時(shí)間,因此地板設(shè)置為溫度邊界條件,分析地板溫度與空氣目標(biāo)溫度差10、20、30 ℃時(shí)的氣流組織。

        表1 簡(jiǎn)化模型與詳細(xì)模型流動(dòng)參數(shù)對(duì)比

        圖8 網(wǎng)格劃分

        4)壁面、大門、頂棚:均設(shè)置為對(duì)流換熱邊界,固體域?yàn)?.2 m 厚聚氨酯泡沫(導(dǎo)熱系數(shù)為0.0237 W/(m·K)),自由來流溫度在極端低溫工況時(shí)為35 ℃,極端高溫工況時(shí)為-10 ℃,對(duì)流換熱系數(shù)為 10 W/(m2·K)。

        2.4 評(píng)價(jià)方法

        為評(píng)估氣流組織的效果,引入溫度不均勻度σ,其按式(4)計(jì)算。

        式中:ti為監(jiān)測(cè)點(diǎn)溫度,℃;t 為監(jiān)測(cè)點(diǎn)溫度平均值,℃;n 為監(jiān)測(cè)點(diǎn)數(shù)量。

        監(jiān)測(cè)點(diǎn)即為工作區(qū)內(nèi)網(wǎng)格的中心點(diǎn),工作區(qū)指距離室內(nèi)壁面2 m,離地面1 m,離風(fēng)口2 m 的有效試驗(yàn)空間。

        3 結(jié)果分析

        3.1 極端低溫

        極端低溫時(shí),送風(fēng)量按換氣次數(shù)6 次/h 給定,送風(fēng)角度θ=65°。當(dāng)?shù)孛鏈囟葹?35 ℃(地面與目標(biāo)室溫差Δtf=20 ℃)時(shí)的室內(nèi)流動(dòng)跡線如圖9 所示??梢姷蜏毓r下,送風(fēng)射流一出送風(fēng)口便快速與周圍空氣混合,同時(shí)在浮力作用下到達(dá)地面,整個(gè)空間流動(dòng)較為均勻,沒有因?yàn)榛仫L(fēng)口都在側(cè)面而出現(xiàn)流動(dòng)的不均勻性。回風(fēng)口附近有較高的風(fēng)速,但其影響范圍較小。

        圖9 Δtf=20 ℃時(shí)極端低溫工況流線

        實(shí)驗(yàn)室中部(y=30 m)處的速度和溫度場(chǎng)分布情況如圖10 所示??梢钥闯觯靡嬗谛黠L(fēng)口良好的摻混效果,送風(fēng)射流速度衰減很快,大部分區(qū)域風(fēng)速在0.3~0.7 m/s 之間,滿足低溫試驗(yàn)要求的不大于1.7 m/s 的要求。溫度分布同樣比較均勻,整個(gè)工作區(qū)內(nèi)溫度介于-57.5~-56.0 ℃之間,沒有明顯的溫度梯度。垂直方向在地表附近溫度梯度較為明顯,但處于工作區(qū)外。

        圖10 Δtf =20 ℃時(shí)y=30 m 剖面速度和溫度分布

        工作區(qū)內(nèi),不同地面溫度或Δtf時(shí),寬度(x 方向)和高度(z 方向)剖面平均溫度ta和溫度不均勻度σ 如圖11 所示??梢钥闯觯孛鏈囟壬邔?dǎo)致ta和σ 同時(shí)增大,這是上送下回的氣流組織形式?jīng)Q定的。在浮力作用下,冷氣流下沉,熱氣流上升,Δtf的增大不僅會(huì)導(dǎo)致地面散熱量的增大,還會(huì)進(jìn)一步增強(qiáng)浮力作用,引起實(shí)驗(yàn)室整體溫度的上升。實(shí)驗(yàn)室整體溫度上升也意味著送、回風(fēng)溫差的變大,導(dǎo)致溫度均勻性變差。

        圖11 不同地面溫度下剖面平均溫度和溫度不均勻度分布Fig.11 Section plane average temperature and un-uniformity under different floor temperature

        地面溫度的影響見表2。從表2 中可知,地面溫度每升高10 ℃,室內(nèi)平均溫度ta升高1.54 ℃,回風(fēng)溫度tr升高1.33 ℃。同時(shí)由于浮力作用增強(qiáng),室內(nèi)平均速度va相應(yīng)增大0.1 m/s。Δtf增大到30 ℃時(shí),回風(fēng)口的溫度甚至低于室內(nèi)平均溫度,這說明送風(fēng)射流與周圍空氣換熱不充分,出現(xiàn)了氣流短路的現(xiàn)象。

        表2 極端低溫時(shí)不同地面溫度下室內(nèi)參數(shù)

        3.2 極端高溫

        極端高溫時(shí),考慮空氣密度變小和浮力作用,將送風(fēng)體積流量提高50%,換氣次數(shù)提高到9 次/h,這時(shí) 保 持 送 風(fēng) 角 度 θ=65°。 地 面 溫 度 為 54 ℃(Δtf=-20 ℃)時(shí)室內(nèi)的流動(dòng)跡線如圖12 所示,送風(fēng)射流受浮力作用而很難到達(dá)地面。

        圖12 Δtf =-20 ℃時(shí)極端高溫工況流線

        實(shí)驗(yàn)室中部(y=30 m)處的速度和溫度場(chǎng)分布情況如圖13 所示。雖然送風(fēng)體積流量加大,但工作區(qū)內(nèi)大部分區(qū)域的風(fēng)速僅為0.3 m/s 左右。離地面5 m以下的空間存在明顯的溫度梯度,達(dá)到了0.6 ℃/m;地面5 m 以上的空間溫度分布很均勻,整個(gè)區(qū)域溫差僅0.5 ℃。

        工作區(qū)內(nèi)不同地面溫度或Δtf時(shí),寬度(x 方向)和高度(z 方向)剖面平均溫度ta和溫度不均勻度σ如圖14 所示。地面溫度降低導(dǎo)致平均溫度ta的降低和σ 增大,但變化輻度遠(yuǎn)比低溫時(shí)小。地面溫度對(duì)地面5 m 以下區(qū)域高度方向的溫度影響非常明顯,對(duì)5 m 以上的空間影響不大。這是因?yàn)榕c低溫工況相反,高溫時(shí)送風(fēng)溫度較高,而呈現(xiàn)出向上流動(dòng)的趨勢(shì),地面溫度較低,導(dǎo)致地面附近空氣溫度降低而向下流動(dòng)。又因送風(fēng)速度衰減較快無法到達(dá)地面,導(dǎo)致地面溫度的影響僅限于地面附近區(qū)域。

        圖13 Δtf =-20 ℃時(shí)y=30 m 剖面速度和溫度分布

        圖14 不同地面溫度下剖面平均溫度和溫度不均勻度分布

        地面溫度的影響見表3,地面溫度對(duì)室內(nèi)平均溫度、回風(fēng)溫度和平均速度的影響并不顯著,但對(duì)地面附近區(qū)域溫度影響很大,導(dǎo)致不均勻度隨地面溫度降低而顯著增大。

        表3 極端高溫時(shí)不同地面溫度下室內(nèi)參數(shù)

        3.3 氣流組織優(yōu)化

        極端高溫工況面臨的主要問題是送風(fēng)射流無法有效抵達(dá)地面,調(diào)整送風(fēng)角度θ=75°,使送風(fēng)射流方向盡可能豎直向下。當(dāng)?shù)孛鏈囟葹?4 ℃(Δtf=-30 ℃)時(shí)室內(nèi)流動(dòng)跡線如圖15 所示,可見氣流流動(dòng)得到有效改善,送風(fēng)射流可以有效抵達(dá)地面。這也導(dǎo)致送風(fēng)射流速度衰減較慢,離出風(fēng)口9 m 處送風(fēng)射流速度才降至1.7 m/s 以下,對(duì)于高度較高的大型試件存在風(fēng)速超標(biāo)的風(fēng)險(xiǎn)。地面附近區(qū)域高度方向溫度梯度明顯改善,地面1~5 m 空間內(nèi),溫度在77.0~77.5 ℃之間,如圖16 所示。

        圖15 θ=75°和Δtf=-30 ℃時(shí)極端高溫工況流線

        圖16 θ=75°和Δtf=-30 ℃時(shí)y=30 m 剖面速度和溫度分布

        送風(fēng)角度的改變,使得送風(fēng)氣流到達(dá)地面,工作區(qū)平均風(fēng)速增大至0.38 m/s,增大了與地面的換熱量。室內(nèi)溫度(ta=77.48 ℃)和回風(fēng)溫度(tr=77.23 ℃)有所降低,但工作區(qū)溫度不均勻度降低至0.38 ℃,甚至優(yōu)于調(diào)整送風(fēng)角度前Δtf=-20 ℃時(shí)的0.40 ℃。地面附近高度方向溫度梯度明顯改善,5 m 高度以上的空間由于送風(fēng)射流摻混作用降低,不均勻度略增大,如圖17 示。

        圖17 θ=75°和Δtf=-30 ℃時(shí)剖面平均溫度和溫度不均勻度分布對(duì)比

        5 結(jié)論

        文中主要對(duì)大型氣候?qū)嶒?yàn)室的氣流組織進(jìn)行了仿真分析,設(shè)計(jì)了旋流風(fēng)口幾何模型,并提出了一種簡(jiǎn)化風(fēng)口CFD 仿真模型。根據(jù)國(guó)內(nèi)外氣候?qū)嶒?yàn)室現(xiàn)狀,設(shè)計(jì)了上送下回的氣流組織形式。對(duì)兩種最極端的試驗(yàn)工況極端低溫和極端高溫下的氣流組織進(jìn)行了仿真分析,研究了送風(fēng)量、送風(fēng)角度、地面溫度對(duì)室內(nèi)溫度場(chǎng)和速度場(chǎng)的影響,得出以下結(jié)論。

        1)采用旋流器作為送風(fēng)口和上送下回的氣流組織形式適用于高度較高的大型氣候?qū)嶒?yàn)室。

        2)受浮力影響,極端低溫和極端高溫工況應(yīng)采用不同的送風(fēng)角度,以使送風(fēng)氣流抵達(dá)地面,降低地面附近的溫度梯度,高溫時(shí)還應(yīng)提高送風(fēng)量。

        3)極端低溫工況下,地面溫度升高將導(dǎo)致室內(nèi)溫度和溫度不均勻度整體上升,甚至出現(xiàn)氣流短路現(xiàn)象。

        4)極端高溫工況下,地面溫度對(duì)地面附近區(qū)域的溫度場(chǎng)影響較大,對(duì)5 m 以上的空間影響不顯著。高溫工況時(shí),增大送風(fēng)角度θ,使氣流方向盡可能向下,將改善地面附近的溫度場(chǎng),但導(dǎo)致送風(fēng)速度衰減較慢,對(duì)于高度較高的試件,試驗(yàn)風(fēng)速有可能超標(biāo),此時(shí)應(yīng)減小θ 并延長(zhǎng)高溫持續(xù)時(shí)間,使地面溫度上升到合理值。

        猜你喜歡
        實(shí)驗(yàn)室模型
        一半模型
        重要模型『一線三等角』
        電競(jìng)實(shí)驗(yàn)室
        重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
        電競(jìng)實(shí)驗(yàn)室
        電競(jìng)實(shí)驗(yàn)室
        電競(jìng)實(shí)驗(yàn)室
        電競(jìng)實(shí)驗(yàn)室
        電競(jìng)實(shí)驗(yàn)室
        3D打印中的模型分割與打包
        一区二区三区在线乱码| 久久久无码人妻精品一区 | 三级日本午夜在线观看| 精品久久一区二区三区av制服| 国产交换精品一区二区三区| 久久精品av在线观看| 女人18毛片a级毛片| 97人人模人人爽人人少妇| 国产在线视频一区二区三区| 精品无码久久久久久久动漫| 久久这里只精品国产2| 宅男久久精品国产亚洲av麻豆| 国产精品国产三级农村妇女| 人妻少妇69久久中文字幕| 电影内射视频免费观看| 人妻献身系列第54部| 国产精品玖玖玖在线资源| 我和丰满老女人性销魂| 美女视频黄a视频全免费网站色 | 三个男吃我奶头一边一个视频| 亚洲精品综合欧美一区二区三区 | 91孕妇精品一区二区三区| 亚欧免费无码AⅤ在线观看 | 国产精品无码久久久久免费AV| 日韩在线手机专区av| 青青操视频手机在线免费观看| 日本二一三区免费在线| av中文字幕潮喷人妻系列| 18成人片黄网站www| 日韩精品区欧美在线一区| 麻豆国产VA免费精品高清在线| 亚洲国产丝袜美女在线| 久久亚洲春色中文字幕久久| 少妇无码av无码专线区大牛影院| 国精品午夜福利视频不卡| 欧美 日韩 国产 成人 在线观看| 久久精品国产99精品九九| 亚洲无码观看a| 91精品国产91综合久久蜜臀| 蜜臀av在线播放一区二区三区| 久久久久人妻精品一区蜜桃|