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

        ?

        翅片陣列對翼面的沖擊換熱影響數(shù)值模擬研究

        2023-01-18 01:14:14張?zhí)靷?/span>王克辰周文武
        上海交通大學學報 2023年1期
        關(guān)鍵詞:駐點塞爾翅片

        張?zhí)靷悾?王克辰, 張 栩, 周文武

        (上海交通大學 機械與動力工程學院, 上海 200240)

        當飛機穿過含有過冷水滴的低溫云層時,機翼表面尤其是前緣會形成積冰,破壞翼型的氣動外形,引起飛行性能下降甚至造成機毀人亡.Petty 等[1]總結(jié)得出:1975—1988年有超過800起由于飛機機翼結(jié)冰造成的安全事故.因此,高效的機翼防除冰技術(shù)對飛行安全至關(guān)重要.在眾多的防除冰方法中,熱氣膜防除冰是一種常見有效的防冰手段,從發(fā)動機引入的高溫空氣先經(jīng)調(diào)壓閥進入笛形管(熱氣供應(yīng)腔),而后通過噴射孔沖擊到機翼內(nèi)表面從而達到防除冰目的.整個過程中,氣膜防除冰效果與機翼沖擊換熱的性能緊密相關(guān),因此提升曲面沖擊換熱性能一直是相關(guān)領(lǐng)域的研究重點.

        國內(nèi)外學者針對平板沖擊換熱開展了大量的研究.前期研究表明:射流孔的幾何形狀、射流孔距靶面的距離、靶面的幾何結(jié)構(gòu)對射流的流動和換熱都有顯著的影響.Jambunathan等[2]總結(jié)了不同雷諾數(shù)、不同噴嘴形狀的平板局部和區(qū)域平均傳熱的實驗結(jié)果和經(jīng)驗關(guān)聯(lián)式.Lytle等[3]利用激光-多普勒測速和壁面壓力測量研究了平板射流沖擊的流動結(jié)構(gòu).Narayanan等[4]通過實驗研究了3種流態(tài)沖擊平板時的流場、表面壓力和換熱率,結(jié)果表明換熱的二次峰位置與近壁的流動脈動速度變化相對應(yīng).Loureiro等[5]考慮了11個不同的實驗數(shù)據(jù)集,總結(jié)了經(jīng)驗表達式,其中包括一個分段努塞爾數(shù)表達式.Siddique等[6]通過數(shù)值模擬確定了不同熱流輸入邊界條件下的局部努塞爾數(shù).除了平板流動和換熱機理的研究,部分學者也探索了在平板表面添加粗糙單元(如三角形的翅片[7-9]、帶矩形槽道[8-9]、針肋[10-12]、微型肋加糙[13]、球形突起[14])對射流沖擊換熱的影響.

        除了平板模型,半圓形曲面沖擊也是重要的研究對象.在改變圓形曲面的幾何結(jié)構(gòu)方面,Singh[15]利用ANSYS Fluent 17.2對射流沖擊復合凹坑和凸凹表面進行了數(shù)值模擬,研究了凹面上凹陷和凸起的交替位置對沖擊傳熱的強化作用.不僅如此,Singh等[16]還建立了一個基于人工神經(jīng)網(wǎng)絡(luò)的替代模型,以確定半圓形凸出物的最佳位置,最大限度地提高凹形表面的傳熱.Gau等[17]通過實驗研究了表面曲率對凹、凸表面沖擊冷卻流動和換熱過程的影響.Lyu等[18]測試了4個凹面(包括一個半圓柱凹面和3個具有不同寬深比的拋物線凹面)的射流沖擊換熱實驗結(jié)果,探究了曲率對凹表面單排人字形噴嘴射流沖擊換熱的影響.此外,劉錫晨等[19]進行了圓形和冠齒形射流口的射流實驗研究,實驗結(jié)果表明冠齒形噴嘴可以增大半圓形凹面的努塞爾數(shù).

        針對曲面沖擊在翼型防除冰的應(yīng)用研究,Hannat等[20]通過ANSYS CFX數(shù)值模擬了三維防冰腔內(nèi)部的流動與換熱,并與測試案例進行了結(jié)果對比.張靖周等[21]采用數(shù)值模擬的方法研究了射流孔直徑、射流孔間距等參數(shù)對凹腔表面溫度分布的影響.郭之強等[22]研究了光滑表面與具有長方體突起結(jié)構(gòu)對內(nèi)凹面的流動和換熱特性的影響,模擬結(jié)果表明凸起結(jié)構(gòu)可以更加均勻地發(fā)散壁面射流.歸曉燁[23]計算曲面槽道內(nèi)部強制對流,運用 Gnielinski 公式對機翼防冰系統(tǒng)槽道進行了分析,得到了槽道截面長度和寬度變化對換熱系數(shù)的影響.關(guān)于沖擊彎曲表面的數(shù)值模擬湍流模型的使用,Seyedein等[24-25]通過研究認為二方程的湍流模型只能模擬出部分的射流流動特性,改進k-ε和k-ωSST模型可以較好地擬合實際的射流沖擊曲面的換熱和流動特性.

        雖然,針對射流沖擊機理以及結(jié)防冰應(yīng)用開展了較多研究,然而之前研究大多以平板或半圓形曲面為主,這與真實的翼型曲面相差較大.因此,本文基于真實翼型結(jié)構(gòu),在幾乎不改動翼型設(shè)計的基礎(chǔ)上,數(shù)值探討如何在內(nèi)壁面添加微小翅片陣列實現(xiàn)曲面沖擊換熱性能提高.以平板和NACA0012對稱翼型為研究模型,首先探索在平板上添加導流翅片陣列的流動換熱性能,而后推廣到曲面翼型.相關(guān)研究成果能為進一步提高機翼防除冰性能提供方法借鑒和設(shè)計參考.

        1 射流沖擊流動結(jié)構(gòu)

        平板沖擊射流是一種被廣泛研究的流動現(xiàn)象,沖擊到平板后從駐點向外擴散的流動通常可以分為3個區(qū)域(見圖1):第1個區(qū)域為自由射流區(qū),在自由射流區(qū)段,來自周圍的流體被卷入射流,從而射流速度u0不斷降低,其垂直速度方向的射流面積逐漸增大;第2個區(qū)域為沖擊區(qū),在該區(qū)域,垂直平板方向速度快速下降并偏轉(zhuǎn),靜壓則快速上升,沖擊區(qū)在平板上的尺寸大約為兩個噴嘴直徑大??;第3個區(qū)域為壁面射流區(qū),射流完全轉(zhuǎn)向,速度平行壁面,隨著流體的發(fā)展,逐漸出現(xiàn)速度邊界層,壁面射流表現(xiàn)出比平行流動更高的傳熱水平,這是由于壁面射流和周圍空氣之間的剪切所產(chǎn)生的湍流,在換熱表面被輸送到邊界層引起的.

        圖1 射流流動分區(qū)Fig.1 Division of jet flow

        在沖擊區(qū)和壁面射流區(qū)的過渡段處增加翅片陣列,起到翅片前擾流、翅片段整流的作用,一方面增大駐點附近的湍流水平,提高了駐點區(qū)域的換熱;另一方面增大遠離駐點區(qū)域的換熱面積和流動速度,提高了遠離駐點區(qū)域的換熱.最終實現(xiàn)整體區(qū)域換熱水平的提升.

        2 平板射流沖擊

        2.1 平板模型

        圖2為平板計算域模型及網(wǎng)格劃分情況.平板沖擊模型和帶有的翅片陣列結(jié)構(gòu)如圖2(a)和2(b)所示,考慮機翼防冰的笛形腔內(nèi)的噴嘴長度與噴嘴直徑相當,因此噴嘴長度和噴嘴直徑d均為4 mm,噴嘴出口距靶面的高度為5d,平板厚度為0.3d.為保證出口流動對沖擊區(qū)域無影響,整體的計算域尺寸需足夠大,設(shè)置為60 mm×60 mm.靶面上的幾何結(jié)構(gòu)除了基準的無翅片以外,還包括8片、12片直翅片和12片彎翅片3種,其中直翅片的長度、高度和厚度分別為d、0.5d和 0.25d,離駐點的距離為1.25d.彎翅形狀為1/4圓弧,圓弧半徑為噴嘴直徑大小.

        圖2 平板計算域模型與網(wǎng)格Fig.2 Flat computing domain model and grid

        2.2 數(shù)值計算

        2.2.1計算域與邊界條件設(shè)置 平板沖擊的計算域同模型,采用ANSYS ICEM進行非結(jié)構(gòu)化網(wǎng)格劃分,包括流體域和固體域,流體域與固體域的交界面節(jié)點一一對應(yīng).從圖2(c)和2(d)中可以看出,在射流區(qū)域、翅片附近進行了網(wǎng)格加密,在靶面的流體一側(cè)增加了壁面棱柱層,棱柱層的第1層厚度為0.005 mm,對應(yīng)的無量綱壁面距離y+<1,網(wǎng)格質(zhì)量在0.3以上.生成網(wǎng)格后在Fluent中進行邊界條件的設(shè)置和湍流模型的選擇,為了保證入口雷諾數(shù)與真實防冰系統(tǒng)射流雷諾數(shù)(Re≈10 000)匹配,入口速度按照不同的雷諾數(shù)(Re=5 500,8 600,10 000,15 000)分別設(shè)置為u0=23.70,36.518,42.61,63.91 m/s,溫度為335 K.流體域的四周為壓力出口,平板固體區(qū)域的下底面為恒熱流密度面,熱流密度為 -2 000 W/m2,四周面為絕熱壁面,設(shè)置參考Dobbertean[9]的模擬過程.空氣的物性參數(shù)求解使用的是理想氣體狀態(tài)方程.流體域計算所使用的湍流模型為Realizablek-ε,Enchanted Wall Function模型,該模型作為改進的k-ε算法,能夠模擬射流撞擊,特別是圓形射流,并且為求解靶面近壁處的流動換熱,采用增強壁面函數(shù)型的近壁面處理.固體域的溫度分布求解為邊界條件已知的控制體導熱微分方程的求解,材料設(shè)置為鋁.采用二階迎風格式的Coupled算法,對整個計算域進行了穩(wěn)態(tài)的數(shù)值模擬.

        2.2.2CFD驗證 共設(shè)計3套不同數(shù)目的網(wǎng)格進行網(wǎng)格獨立性的驗證,網(wǎng)格數(shù)分別為120萬、225萬和322萬.網(wǎng)格獨立性驗證的結(jié)果表明,當網(wǎng)格數(shù)大于225萬時,壁面上的溫度隨著網(wǎng)格加密的變化在0.3 K以內(nèi),綜合考慮計算精度和模擬效率,最終選擇網(wǎng)格數(shù)目225萬的網(wǎng)格.

        為了確保網(wǎng)格和湍流模型的準確性,本文將數(shù)值模擬得到的駐點位置處的的局部努塞爾數(shù)大小與Garimella等[26]通過實驗歸納的實驗關(guān)聯(lián)式計算值進行了比較.局部努塞爾數(shù)的定義式如下:

        (1)

        式中:h0為局部的對流換熱系數(shù);λ為空氣的導熱系數(shù);qint為交界面處的熱流密度;Tint和Tj分別為交界面固體的溫度和射流的入口溫度.

        Garimella等[26]歸納的實驗關(guān)聯(lián)式如下:

        Nu0=0.462Re0.585Pr0.4(H/d)0.024

        (2)

        式中:Pr為空氣的普朗特數(shù);H為垂直方向上離靶面的高度.式(2)的適用范圍為Re=4 000~23 000,H/d=1~5.

        本文數(shù)值模擬得到4種雷諾數(shù)下平板駐點位置處的Nu0與實驗關(guān)聯(lián)式所得的結(jié)果對比如表1所示.表中:Er為本文模擬所得Nu0與實驗關(guān)聯(lián)式所得Nu0之間的相對誤差.由表可見,Er在10%以內(nèi),其中低雷諾數(shù)(Re=5 500)下數(shù)值模擬所得Nu0大于實驗關(guān)聯(lián)式計算的Nu0,高雷諾數(shù)下模擬所得Nu0小于實驗關(guān)聯(lián)式的Nu0.這是由于隨著雷諾數(shù)的增大, 流動再循環(huán)和混合更強,下游的熱空氣回流到駐點附近,導致Nu0的數(shù)值隨著表面熱流密度的變化有輕微的改變[27],而本文在流固交界面處的熱流密度與文獻[26]中實驗的熱流密度不同.

        表1 不同雷諾數(shù)下駐點位置的努塞爾數(shù)Tab.1 Nusselt numbers of stagnation point at different Reynolds numbers

        2.2.3結(jié)果分析 為了更加綜合地評估增加翅片之后的換熱效能,流固交界面處的局部努塞爾數(shù)是不適用的,參考Dobbertean[9]對有厚度平板的努塞爾數(shù)的定義,定義綜合的努塞爾數(shù):

        (3)

        式中:h為綜合定義的對流換熱系數(shù);qcon為固體外壁面的熱流密度,該值在平板各處為一個恒定值;To為固體外壁面的溫度.該式從研究對象來看將整個平板的換熱進行了考慮,包括外表面的熱量流失、交界面的對流換熱、交界面到外表面的熱傳導,因此是一個綜合的換熱效能.

        針對新定義的Nu,首先對比了4種結(jié)構(gòu)在Re=5 500 時的結(jié)果,如圖3所示.圖中:x、y分別為以射流孔圓心為原點,X和Y方向下的距離.由圖可見,加8片、12片直翅片相比于不加翅片均可以提升換熱能力,并且8片的換熱效果最好.這是因為在12片直翅片的結(jié)構(gòu)中,射流在向四周擴散時發(fā)生了堵塞,12片相較于8片周向面積減小了17.1%.加12片彎翅片相對于無翅片的換熱能力反而減弱,這是由于彎翅片不僅存在直翅片的堵塞問題,而且翅片是彎曲的,氣流在每個翅片的吸力面后緣均會出現(xiàn)明顯的分離,使得氣流更難向整個平板擴展出去,翅片換熱面積增大所增加的換熱小于因為射流無法擴散而導致的熱量堆積,所以換熱能力反而減弱了.

        針對圖3分析得到的結(jié)果,后續(xù)對不同雷諾數(shù)換熱性能的分析,僅使用無翅片和增強換熱效果最優(yōu)的8片直翅片的結(jié)果進行比較,如圖4所示.從圖4(a)、4(b)和4(c)、4(d)中最小和最大兩個雷諾數(shù)工況可以看出,安裝直翅片相較于無翅片的平板,Nu都是有提高的,這也意味著恒熱流外壁面的溫度有提升.

        圖3 4種結(jié)構(gòu)的努塞爾數(shù)分布云圖Fig.3 Nusselt number distribution of 4 kinds of structures

        圖4 不同雷諾數(shù)下的努塞爾數(shù)分布云圖Fig.4 Nusselt number distribution at different Reynolds numbers

        從(0,0)駐點位置處沿與X軸正方向22.5°方向取一條線上的Nu,得到一系列點線圖,即不同雷諾數(shù)下的Nu分布圖,如圖5所示.圖中:r為沿22.5°線上的點到靶面中心的距離;η為相對于無翅片時,有翅片Nu的提升比例.由圖可以直觀地看出換熱效能的提升效果,隨著雷諾數(shù)的增大,Nu隨之增大.相同雷諾數(shù)下,在徑向距離5d以內(nèi),有翅片的Nu均高于無翅片平板的Nu.在低雷諾數(shù)(Re=5 500)下,η最大,平均為5.5%左右,其余3個雷諾數(shù)下的η在1.3%左右.這是由于在低雷諾數(shù)下?lián)Q熱效果由翅片面積帶來的提升與翅片增強局部對流換熱系數(shù)帶來的提升都比較明顯,而隨著雷諾數(shù)的增大,后者帶來的提升逐漸占據(jù)主導地位,所以整體提升效果沒有低雷諾數(shù)時那么明顯.此外,由于翅片的存在導致向外擴散的周向面積減小,從而向外擴散的射流流速增快,換熱增強,徑向距離在翅片之外的區(qū)域換熱也有提升.

        圖5 綜合努塞爾數(shù)對比Fig.5 Comprehensive Nusselt number comparison

        3 翼型曲面射流沖擊

        3.1 翼型計算模型與網(wǎng)格

        熱氣膜防除冰沖擊的壁面是機翼前緣的內(nèi)表面,在平板上確認增加翅片對于換熱有提升,將該結(jié)構(gòu)應(yīng)用到翼型曲面上進行驗證.在翼型曲面模型中,整體計算域的寬度、射流孔距駐點的高度以及翼型曲面板的厚度保持不變,在曲面模型中研究了翅片長度對換熱的影響,在原先1倍噴嘴直徑長度的翅片基礎(chǔ)上,增加了2倍噴嘴直徑長度的翅片模型.整體網(wǎng)格如圖6所示,圖中:S為沿著外壁面上與葉片型線相切的方向.同樣在近壁面和翅片附近進行了加密,邊界條件和湍流模型的設(shè)置與平板的模型相同.

        圖6 曲面計算域模型與網(wǎng)格Fig.6 Curved surface computing domain model and grid

        3.2 翼型曲面沖擊換熱結(jié)果

        圖7為無翅片、8片1倍噴嘴直徑長度翅片和2倍噴嘴直徑長度翅片的曲面板外壁面在Re=5 500 工況下的綜合努塞爾數(shù)分布,這里將翼型曲面按照圖6所示的型線S方向進行展平.圖中:s為翼型外表面上的點沿型線S方向距翼型前緣的距離;z為以射流孔圓心為原點,Z方向的距離.增加翅片后紅色的高溫區(qū)域明顯增多,并且2倍的翅片的效果要優(yōu)于1倍的.因此,后續(xù)不同雷諾數(shù)下的流動換熱的分析在無翅片和2倍噴嘴直徑長度翅片的模型之間展開.

        圖7 3種模型的綜合努塞爾數(shù)Fig.7 Comprehensive Nusselt number distribution of three models

        翼型曲面板模型同樣研究了4種雷諾數(shù)下的換熱性能.圖8對比了不同雷諾數(shù)下的無翅片和有翅片的Nu分布,中間兩個雷諾數(shù)類似,因此選擇最大、最小兩個雷諾數(shù)的結(jié)果進行展示.在4種雷諾數(shù)的工況下,加翅片陣列模型的Nu在各個位置均要高于無翅片模型的Nu.在高雷諾數(shù)Re=15 000的時候,換熱效果的提升最為明顯,其現(xiàn)象與平板的結(jié)果相反.

        圖8 不同雷諾數(shù)下的努塞爾數(shù)分布云圖(曲面)Fig.8 Nusselt number distribution at different Reynolds numbers (curved surface)

        圖9顯示了從(s,z)=(0, 0)駐點位置處延展的兩條曲線上的綜合努塞爾數(shù)及其有2倍噴嘴直徑長度翅片時相對無翅片的提升比例η.從圖中可以看出,隨著雷諾數(shù)的增大,綜合努塞爾數(shù)相應(yīng)提升,并且?guī)С崞膿Q熱始終要優(yōu)于不帶翅片的.整體的提升比例在4%~10%左右,比平板的提升更為明顯.

        結(jié)合圖8和圖9觀察到的現(xiàn)象,在翼型曲面上,隨著雷諾數(shù)的增大,綜合努塞爾數(shù)變化規(guī)律與平板不同,平板在小雷諾數(shù)下提升效果更明顯,而翼型曲面在高雷諾數(shù)下效果更好;此外,整體提升效果也不同,平板提升比例明顯小于翼型曲面.造成這種現(xiàn)象的原因有:① 曲面的曲率改變了射流的流場結(jié)構(gòu),從而對對流換熱的效果產(chǎn)生了影響;② 在翼型模型中采用的是2倍噴嘴直徑長度的翅片,而在平板上為1倍噴嘴直徑長度的翅片,翅片的長度從增大換熱面積上有不同的提升效果.

        圖9 曲面綜合努塞爾數(shù)對比Fig.9 Curved surface comprehensive Nusselt number comparison

        3.3 翼型曲面沖擊流場

        翅片陣列對流場的影響對于平板和曲面的模型是類似的,均可以提高沖擊壁面附近的湍流度.湍流水平的提升可以導致傳熱率的增大,為了進一步說明換熱得到改善的內(nèi)在機理,對Re=8 600 的工況進行了瞬態(tài)的數(shù)值仿真計算.圖10給出了距離曲面板內(nèi)壁面0.5 mm處的近壁曲面上的湍動能強度k以及X、Y、Z方向上的速度u、v、w的均方根誤差(RMSE),RMSE越大,代表速度的波動越大,即在數(shù)值模擬計算中速度的脈動量越大.翅片的存在使得k增大了1倍多;并且主要改變了Y、Z兩個方向上的速度脈動量,這兩個方向也正是平行于曲面的兩個方向.說明翅片的結(jié)構(gòu)在駐點附近對沖擊射流的流場起到了擾流和增強湍流度的作用.

        圖11為翅片方向上的時均流線圖.圖中:|V| 為流場合速度的大小.由放大的流線圖可見,射流在到達翅片時與翅片發(fā)生了撞擊,速度方向發(fā)生了改變,翅片前緣同時也是四周擴散流動的駐點位置,存在高水平的換熱.不過整體流場的結(jié)構(gòu)沒有太大的變化,結(jié)合圖10的信息,可知翅片影響的是每一個時刻的流場,特別是平行于翼型曲面方向上的速度脈動量,而翅片對整體時均流場的影響比較小.

        圖10 貼近內(nèi)壁面的速度均方根誤差及湍動能(Re=8 600)Fig.10 RMSE of velocity and turbulent kinetic energy close to inner wall (Re=8 600)

        圖11 截面流線圖Fig.11 Cross-section streamline diagram

        4 結(jié)論

        為了提高防除冰目標曲面的射流沖擊換熱,研究了翅片結(jié)構(gòu)對平板和翼型曲面板上射流沖擊換熱的影響,對比不同雷諾數(shù)下的Nu,研究結(jié)果表明:

        (1) 直翅片陣列的結(jié)構(gòu)在平板可以提高射流沖擊的換熱強度,8片的提升效果最好,12片次之,8片1倍噴嘴直徑長度直翅片相較于無翅片可以提升1%~6%.

        (2) 翅片的結(jié)構(gòu)在翼型板上同樣可以提高射流沖擊的換熱強度,8片2倍噴嘴直徑長度長直翅片在曲面板上的提升為4%~10%,換熱水平的提升可以提高防除冰效果.雖然增加翅片在平板和翼型曲面上均可以提升換熱強度,不過,由于曲面曲率的影響,射流在曲面上的流動傳熱和平板是有區(qū)別的,后續(xù)的研究將針對真實的翼型進行探索.

        (3) 翅片結(jié)構(gòu)不僅可以增大駐點以及周邊區(qū)域的湍動能強度從而提升換熱性能(湍動能的提升達到了1倍左右),也可以從增加換熱面積的角度進一步提高換熱.

        (4) 低雷諾數(shù)下,翅片通過改變流動結(jié)構(gòu)增加的換熱量與增大換熱面積增加的換熱量都是比較明顯的;高雷諾數(shù)下,則以改變流動結(jié)構(gòu)的提升為主.

        猜你喜歡
        駐點塞爾翅片
        垂直翅片管自然對流傳熱特性的數(shù)值研究
        機械工程師(2022年6期)2022-06-21 08:44:24
        如果地球被我們吃掉了
        大功率COB-LED的翅片散熱器優(yōu)化研究
        超硬翅片滾刀加工在CNC磨床上的實現(xiàn)
        基于游人游賞行為的留園駐點分布規(guī)律研究
        中國園林(2018年7期)2018-08-07 07:07:48
        淺談點-連式ATP系統(tǒng)在埃塞爾比亞輕軌中的應(yīng)用
        大功率LED燈的新型底部開縫翅片散熱性能研究
        利用遠教站點,落實駐點干部帶學
        利用遠教站點,落實駐點干部帶學
        2300名干部進村“串戶”辦實事
        源流(2015年8期)2015-09-16 18:01:32
        性饥渴的农村熟妇| av在线男人的免费天堂| 久久伊人精品中文字幕有| 亚洲人精品午夜射精日韩| 人人妻人人澡人人爽欧美二区| 国产精品玖玖玖在线资源| 国内精品嫩模av私拍在线观看| 男女无遮挡高清性视频| av无码精品一区二区三区宅噜噜| 亚洲人成网站在线播放观看| 久久国产劲爆内射日本| 国产一区二区三区在线观看完整版 | 日本中国内射bbxx| 人妻少妇精品视中文字幕国语| 国产精品女同久久免费观看 | 中文字幕人妻系列一区尤物视频| 女同av一区二区三区| 欧美大片aaaaa免费观看| 国产乱沈阳女人高潮乱叫老| 日本女优一区二区在线免费观看| 日本一区二区三区视频免费观看| 国产成人精品一区二区三区视频 | 成熟妇女毛茸茸性视频| 天堂中文最新版在线中文| 国产黄色免费网站| 亚洲一本之道高清在线观看| 职场出轨的人妻中文字幕| 94久久国产乱子伦精品免费| 国产aⅴ天堂亚洲国产av| 国产av一区二区日夜精品剧情| 亚洲日韩精品a∨片无码加勒比 | 国产一区二区三区不卡在线播放| 精品无码久久久久久久久水蜜桃| 国产剧情麻豆女教师在线观看| 国产精品不卡无码AV在线播放| 国产精品亚洲精品一区二区 | 无码专区久久综合久中文字幕 | 一区二区三区日本视频| 国产精品久久久久高潮| 免费啪啪视频一区| 中文字幕亚洲精品一二三区|