999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

隧道窯冷卻帶流-固耦合散熱數(shù)值模擬

2016-05-31 01:06:32許寧王雪峰張柱銀黃建軍鹽城工學(xué)院機(jī)械工程學(xué)院江蘇鹽城05湖南工業(yè)大學(xué)機(jī)械工程學(xué)院湖南株洲007中國建材西安墻體材料研究設(shè)計院陜西西安7006中國建筑材料集團(tuán)有限公司北京0006
中國陶瓷工業(yè) 2016年2期
關(guān)鍵詞:散熱

許寧,王雪峰,張柱銀,王 君,黃建軍,肖 慧,郅 曉(. 鹽城工學(xué)院機(jī)械工程學(xué)院,江蘇 鹽城 05;. 湖南工業(yè)大學(xué)機(jī)械工程學(xué)院,湖南 株洲 007;. 中國建材西安墻體材料研究設(shè)計院,陜西 西安 7006;. 中國建筑材料集團(tuán)有限公司,北京 0006)

?

隧道窯冷卻帶流-固耦合散熱數(shù)值模擬

許寧1, 2,王雪峰1,2,張柱銀2,王 君2,黃建軍1,肖 慧3,郅 曉4
(1. 鹽城工學(xué)院機(jī)械工程學(xué)院,江蘇 鹽城 224051;2. 湖南工業(yè)大學(xué)機(jī)械工程學(xué)院,湖南 株洲 412007;3. 中國建材西安墻體材料研究設(shè)計院,陜西 西安 710061;4. 中國建筑材料集團(tuán)有限公司,北京 100036)

摘 要:以某燒結(jié)磚隧道窯為研究對象,為獲得隧道窯冷卻帶散熱溫度場數(shù)據(jù),建立了隧道窯窯壁-磚塊-流體耦合傳熱物理模型與數(shù)學(xué)模型,采用FLUENT流體計算軟件對其進(jìn)行了隧道窯窯壁-磚塊-流體耦合傳熱仿真,獲得了窯內(nèi)整體比較精確直觀的溫度場分布。仿真的結(jié)果表明:風(fēng)速為5 m/s、孔隙率為0.5以及磚塊適當(dāng)稀碼時,隧道窯冷卻帶降溫效果較好,其與隧道窯燒成曲線冷卻段吻合也較好,為風(fēng)速調(diào)整、碼磚方式和孔隙率的選擇提供了理論依據(jù),提高了燒成質(zhì)量和生產(chǎn)效率。

關(guān)鍵詞:隧道窯冷卻帶;流-固耦合;散熱;數(shù)值模擬

E-mail:xuning196402@163.com

0 引 言

隧道窯是現(xiàn)代化的連續(xù)式燒成的熱工設(shè)備,廣泛用于磚塊和陶瓷產(chǎn)品的焙燒生產(chǎn),燒結(jié)磚隧道窯冷卻帶在燒成帶之后,主要的作用是使焙燒好的磚塊迅速冷卻。經(jīng)常在窯尾安裝風(fēng)機(jī),通過改變風(fēng)速的大小從而控制磚塊冷卻的速率與溫度。磚塊內(nèi)部傳熱機(jī)理是比較復(fù)雜的,涉及流相及固相的熱傳導(dǎo)和對流換熱[1]。磚塊上下溫差受冷卻風(fēng)影響較大,出磚時磚塊表面容易產(chǎn)生大的裂紋,影響磚塊的質(zhì)量。

隨著微型計算機(jī)技術(shù)的發(fā)展以及計算方法的不斷更新進(jìn)步,計算機(jī)應(yīng)用于燒結(jié)磚隧道窯冷卻帶的三維流-固耦合數(shù)值模擬成為現(xiàn)實。有關(guān)流固耦合的計算方法大致可以分為兩類:一類是解析解;一類是數(shù)值解。數(shù)值解法通常采用的是有限單元法或有限差分法。Anton is Antoniou等人計算了二維流-固耦合流場[2]。梅國棟等采用有限差分法分析三維流-固耦合問題[3]。浙江大學(xué)李迎博士建立了內(nèi)燃機(jī)活塞組-缸套-冷卻水-機(jī)體流固耦合系統(tǒng)。對柴油機(jī)穩(wěn)態(tài)工況下的耦合傳熱、周期性瞬態(tài)傳熱以及冷啟動、外特性的瞬態(tài)傳熱與流動進(jìn)行了數(shù)值仿真[4]。Ed Bendel主要通過測量排氣管道表面的溫度值作為邊界條件進(jìn)行數(shù)值模擬,并且利用Rad Therm和Star-CD的方法來分析對流、傳導(dǎo)的換熱特性[5]。劉上在FLUENT軟件平臺上,采用氣固耦合傳熱計算的方法模擬了二維擴(kuò)散燃燒實驗器中的溫度、組分分布[6]。大連理工大學(xué)周禛建立了單缸的缸蓋- 缸套 - 冷卻水耦合傳熱模型并且進(jìn)行了耦合傳熱仿真[7]。

本文運用FLUENT軟件對隧道窯冷卻帶流-固耦合散熱進(jìn)行數(shù)值模擬。模擬不同的風(fēng)速、孔隙率、氣體流動的數(shù)學(xué)模型以及磚塊碼放方式對隧道窯散熱的影響。為磚塊的散熱找到最佳的風(fēng)速、孔隙率、數(shù)學(xué)模型以及碼磚的方式。為降低磚塊的裂縫,提高磚塊的質(zhì)量尋找有效途徑。

1 隧道窯冷卻帶的物理模型

在GAMBIT中建立了隧道窯冷卻帶的物理模型。為了減少計算量,提高運算效率對模型進(jìn)行適當(dāng)?shù)暮喕x取六節(jié)窯車進(jìn)行流-固耦合散熱的數(shù)值模擬。隧道窯冷卻帶的截面寬度為Y=3600 mm,窯車頂面到窯頂之間的距離Z=1200 mm,冷卻帶的長度為X=8000 mm。對模型進(jìn)行假設(shè):冷卻帶要求保溫性良好;窯壁可以假設(shè)為絕熱;窯尾的風(fēng)機(jī)可以簡化為速度入口;出口設(shè)為壓力出口;不考慮磚塊之間的傳熱影響以及磚垛之間的輻射傳熱影響。

圖1 隧道窯冷卻帶的物理模型Fig.1 Physical model of the tunnel kiln cooling zone

隧道窯冷卻帶的物理模型分為磚塊區(qū)、流體區(qū)。根據(jù)數(shù)值模擬要求采用Tet/Hybrid混合網(wǎng)格劃分方法,即結(jié)構(gòu)網(wǎng)格和非結(jié)構(gòu)網(wǎng)格相結(jié)合的劃分方法。對磚塊區(qū)的網(wǎng)格劃分,格間距為4,流體區(qū)的網(wǎng)格劃分格間距為6。總計網(wǎng)格總數(shù)為670155。檢查網(wǎng)格質(zhì)量確保minimum volume不為負(fù)數(shù)。采用有限體積法對控制方程進(jìn)行離散化,動量方程和能量方程用一階迎風(fēng)格式,迭代直至收斂。合理的劃分網(wǎng)格有利于減少計算時間節(jié)約計算成本。

2 隧道窯冷卻帶的數(shù)學(xué)模型

本文計算中運用的控制方程為三維常物性不可壓縮流動的連續(xù)性方程及能量方程。計算域為磚塊與空氣、粉塵流動的區(qū)域。

2.1標(biāo)準(zhǔn)的k-ε模型

式中,Gk是由于平均速度梯度引起的湍動能產(chǎn)生;Gb是由于浮力影響引起的湍動能產(chǎn)生;YM為可壓縮湍流脈動膨脹對總的耗散率的影響;C1ε、C2ε、C3ε為經(jīng)驗常數(shù),F(xiàn)LUENT中默認(rèn)值為C1ε=1.44、C2ε=1.92、C3ε=0.09;σk=1.0、σε=1.3。

2.2流固耦合控制方程簡化

隧道窯冷卻帶內(nèi)氣體的流動看做粘性不可壓縮的湍流流動。它在流動與傳熱的過程中遵循三大守恒方程:質(zhì)量守恒方程、動量守恒方程(N-S方程)、能量守恒方程。

質(zhì)量守恒方程

最后在求解過程中采用非耦合隱式格式(即SIMPLE算法)。該算法占用的內(nèi)存小,計算快捷靈活。用有限容積法導(dǎo)出的離散方程可以保證具有守恒特性,而且離散方程系數(shù)的物理意義明確,是目前流動與傳熱問題的數(shù)值計算中應(yīng)用最廣的一種方法。

2.3模型邊界條件及參數(shù)設(shè)定

在對隧道窯冷卻帶三維流固耦合散熱數(shù)值模擬的研究中,將流體的進(jìn)口設(shè)置為速度入口,流體初始速度方向要與源面垂直;流體的出口設(shè)定為壓力出口[8-10];隧道窯的窯壁在整個計算的過程中設(shè)為絕熱壁;內(nèi)部窯壁與磚塊的區(qū)域設(shè)為流體區(qū)域[11],流體不可壓縮,流體中有氣相和固相,將固相中的顆粒假設(shè)為密度相同直徑均勻的顆粒[12],忽略粒子之間的相互作用。磚塊的密度為2500 kg/m3,導(dǎo)熱系數(shù)為80 w/m.k,熱生成率為0.57。

3 數(shù)值模擬結(jié)果與分析

3.1隧道窯冷卻帶流-固耦合溫度場分布

根據(jù)FLUENT軟件數(shù)值模擬計算經(jīng)驗,流體從隧道窯窯尾以5 m/s左右的速度進(jìn)入,進(jìn)入窯內(nèi)流體的溫度為常溫。

圖2是從窯長方向取得XY截面,從圖中可以看出,顏色越深代表溫度越高,氣體從窯尾進(jìn)入對磚塊進(jìn)行冷卻,磚塊運動的方向與氣體運動方向相反,散熱溫度隨窯車運動的方向逐漸降低。截面內(nèi)溫度分層情況比較均勻,磚塊的左右平均溫差不超過42 ℃。且在磚垛與窯壁之間或者磚垛與磚垛之間的間隙冷卻的速度較快,主要因為磚垛與窯壁或者磚垛與磚垛之間的阻力較小氣體流速較快。

圖2 Z=0.5m截面溫度場分布圖Fig.2 Temperature feld distribution on the cross-section of Z = 0.5 m

圖3 X=3.6m截面溫度場分布圖Fig.3 Temperature feld distribution on the cross-section of X = 3.6 m

圖4 X=0.8m 截面溫度場分布圖Fig.4 Temperature feld distribution on the cross-section of X = 0.8 m

圖3、圖4分別是沿著窯寬方向取得兩個YZ截面。圖5是沿著窯長方向取得XZ截面。從圖中都可以看到窯頂?shù)臏囟纫雀G車臺面的溫度低,溫度的變化梯度均勻。氣體從窯體進(jìn)入磚垛時速度會加大,這是由于進(jìn)入窯體的風(fēng)速為5 m/s,馬赫數(shù)小于0.2可以認(rèn)為氣體的流動是低速,不考慮氣體的壓縮性,氣體的速度與截面積成反比。氣體進(jìn)入磚垛時截面積變小則速度增大。氣體與固相邊界(即:磚垛)相互流動時還會產(chǎn)生湍流。在整個流-固耦合散熱的過程中,磚垛附近產(chǎn)生較強(qiáng)的湍流影響冷卻速度,在磚垛與窯頂之間主要是層流,冷卻的速度較快,則出現(xiàn)窯頂?shù)臏囟纫雀G車臺面附近區(qū)域的溫度低。

3.2計算結(jié)果分析

為了研究不同風(fēng)速對磚塊的平均溫度和上下溫差影響。選擇了10組不同風(fēng)速進(jìn)行數(shù)值模擬,計算殘差取為10-3,在計算迭代過程中,當(dāng)監(jiān)測點的溫度呈現(xiàn)穩(wěn)定不變或者近似穩(wěn)定不變時,即認(rèn)為收斂。對比結(jié)果如下:

從圖6可以看出,要保證隧道窯冷卻帶磚塊冷卻的溫度與速度,同時還要保持窯內(nèi)截面的溫差最小,就要控制進(jìn)入窯尾冷卻風(fēng)的風(fēng)速。風(fēng)速過低(風(fēng)速小于5 m/s)時,隨著氣體流速升高溫度逐漸降低,截面上下溫差也在逐漸降低。風(fēng)速達(dá)到5 m/s時,氣體與磚垛達(dá)到產(chǎn)生湍流的條件,此時截面的上下溫差達(dá)到最低。風(fēng)速繼續(xù)增加時,湍流的強(qiáng)度也隨著風(fēng)速增加逐漸降低,此時磚垛附近區(qū)域的散熱平均溫度下降逐漸變緩,截面的上下溫差也在不斷的增大。

從圖7可以看出,孔隙率不斷增大,磚塊碼放逐漸變得稀疏。Z=0.5 m時,XY截面的平均溫度逐漸降低。XY截面的上下溫差先減小后增大再減小的過程,主要是因為孔隙率小于0.5時,隨著孔隙率的逐漸增加流體與磚垛的相互流動時還沒有產(chǎn)生湍流,溫差還在逐漸減小。當(dāng)孔隙率在0.5-0.7之間時,湍流逐漸形成并且逐漸加強(qiáng),影響散熱溫差也在逐漸增大。當(dāng)孔隙率大于0.7時,磚塊已經(jīng)很稀疏,湍流現(xiàn)象逐漸消失,截面的溫差也在逐漸降低。但考慮生產(chǎn)的產(chǎn)量與成本選取孔隙率為0.5較為合適。

從圖8可以看出,二種不同氣體的流動模型冷卻曲線與隧道窯燒成曲線冷卻段曲線對比結(jié)果,兩條冷卻曲線都與隧道窯燒成曲線冷卻段曲線接近,但采用模型方式對磚塊進(jìn)行冷卻曲線更接近隧道窯冷卻曲線。主要因為氣體在磚垛內(nèi)流動的方式還是湍流,同時在做數(shù)值模擬時對物理模型和數(shù)學(xué)模型都進(jìn)行了適當(dāng)?shù)暮喕吔鐥l件以及參數(shù)設(shè)定與實際有點差別,計算的方法也不夠完善。

圖5 Y=1.7m 截面溫度場分布圖Fig.5 Temperature feld distribution on the cross-section of Y = 1.7 m

圖6 Z=0.5m時不同風(fēng)速對XY截面平均溫度和溫差的影響Fig.6 Effect of different draft speed on average temperature and temperature difference on XY cross-section when Z = 0.5 m

圖7 Z=0.5m時不同孔隙率對XY截面平均溫度和溫差的影響Fig. 7 Effect of different pore rate on average temperature and temperature difference on XY cross-section when Z = 0.5 m

圖8 不同氣體流動模型冷卻曲線與隧道窯燒成曲線冷卻段對比Fig.8 Contrast of the cooling curves of different fuid-fow models with the cooling segment of the tunnel kiln fring curve

圖9 不同碼放方式對磚塊散熱溫度變化的影響Fig.9 Effect of different bricks arrangement on brick heat transfer

從圖9可以看出,隨著風(fēng)速加快Z=0.5 m的XY截面平均溫度下降,在相同風(fēng)速、孔隙率下磚塊稀碼比磚塊密碼的平均溫度低。從圖中曲線下降的斜率可以看出:磚塊稀碼溫度隨著風(fēng)速下降的速度快于磚塊密碼。主要因為磚塊稀碼時,氣體與磚垛接觸面積加大,溫度降低速度也會加快,平均溫度也會相對較低。磚塊稀碼時氣體在磚垛內(nèi)部的湍流強(qiáng)度小于磚塊密碼,其散熱效果也要優(yōu)于磚塊密碼。

4 結(jié) 論

本文運用FLUENT對隧道窯冷卻帶六節(jié)窯車段進(jìn)行流-固耦合散熱數(shù)值模擬,結(jié)論如下:

(1)用數(shù)值模擬的方法對隧道窯冷卻帶磚塊的散熱進(jìn)行模擬,可以代替大量的反復(fù)實驗,節(jié)約大量的人力物力縮短研究時間。

(2)文中分別從風(fēng)速、孔隙率、碼磚方式的角度展開研究。研究表明:風(fēng)速5 m/s,孔隙率0.5時磚塊散熱時上下溫差較小,平均溫度較低;在相同風(fēng)速與孔隙率的條件下磚塊稀碼比磚塊密碼的散熱條件好。

(3)選擇合適的數(shù)學(xué)計算模型可以有效的模擬隧道窯冷卻帶散熱的過程,為以后的實驗研究提供理論依據(jù)。

(4)文中的散熱數(shù)值模擬均基于實際生產(chǎn)條件,模擬的結(jié)果與實際數(shù)據(jù)相吻合,表明文章中數(shù)值模擬的結(jié)果是可行的,接下來進(jìn)一步的工作將模擬的結(jié)果應(yīng)用于實際生產(chǎn)中,針對生產(chǎn)的實際情況做適當(dāng)?shù)恼{(diào)整。

參考文獻(xiàn):

[1] 百敏麗, 薛瑩瑩, 呂繼祖, 等. 隔熱材料導(dǎo)熱系數(shù)的數(shù)值模擬預(yù)測[J]. 熱科學(xué)與技術(shù), 2007, 6(2): 136-140.

[2] AKYUZLU K M, ANTONIOU A. Determination of regression rate in an ablation hybrid rocket solid fuel using a physics based comprehend div mathematical model. AIAA, 2002: 2002-3577. [3] 梅國棟, 王云海. 三維流固耦合數(shù)值模擬在銅鑼山隧道安全性評價中的應(yīng)用[J]. 中國安全生產(chǎn)科學(xué)技術(shù), 2009, 5(6): 57-61. [4] 李迎. 內(nèi)燃機(jī)流固耦合傳熱問題數(shù)值仿真與應(yīng)用研究[D]. 浙江大學(xué), 2006.

[5] BENDELL E. Investigation of a coupled CFD and thermal modeling methodology for prediction of vehicle underbody temperature. SAE, 2005-01-2044.

[6] 劉上, 孫得川. 固液發(fā)動機(jī)實驗燃燒器的氣固耦合傳熱計算[J]. 推進(jìn)技術(shù), 2008, 29(3): 258-261.

[7] 周禛, 薛冬新, 唐斌, 等. 某中速柴油機(jī)冷卻液流動及流固耦合傳熱計算分析[J]. 柴油機(jī), 2013, 35(1): 17-32.

[8] KAYA S, MANCHAN E. Modelling and optimization of the firing zone of a tunnel kiln to predict the optimal feed locations and mass fluxes of the fuel and secondary air. Applied Energy,2009, 86: 325-332.

[9] NICOLAU V P, DADAM A P. Numerical and experimental thermal analysis of a tunnel kiln used in ceramic production. Journal of the Brazilian Society of Mechanical Sciences & Engineering 2009, 31 (4): 297-304.

[10] SANTOS G M. Study of thermal behavior of a tunnel kiln used in red ceramic industry. Master Thesis, Federal University of Santa Catarina, Santa Catarina, Brazil, 2001.

[11] NACCACHE M F, GOMES M S P. Numerical simulation of flow and heat transfer through a tunnel kiln. 18th International Congress of Mechanical Engineering, November 6-11, 2005.

[12] OBA R, POSSAMAI T S. Thermal analysis of tunnel kiln used to produce roof tiles. Applied Thermal Engineering, 2013. 63(1): 59-65.

通信聯(lián)系人:許寧,男,博士,教授。

Received date:2015-12-18. Revised date: 2015-12-20.

Correspondent author:XU Ning, male, Ph.D., Professor.

Numerical Simulation of Fluid - Solid Coupling Heat Transfer in Tunnel Kiln Cooling Zone

XU Ning1, 2, Wang Xuefeng1, 2, ZHANG Zhuyin2, WANG Jun2, HUANG Janjun1, XIAO Hui3, ZHI Xiao4
(1. School of Mechanical Engineering, Yancheng Institute of Technology, Yancheng 224051, Jiangsu, China; 2. School of Mechanical Engineering, Hunan University of Technology, Zhuzhou 412007, Hunan, China; 3. Xi'an Research & Design Institute of Wall & Roof Materials, Xi'an 710061, Shaanxi, China; 4. China National Materials Group Corporation, Beijing 100036, China)

Abstract:In order to obtain the cooling temperature feld data in the tunnel kiln cooling zone, the tunnel kiln wall bricks-fuid coupling heat transfer physical model and mathematic model were established, taking a sintered brick tunnel kiln as the research object. The tunnel kiln wall bricks-fuid coupling heat transfer simulations were performed by commercial software FLUENT, and the overall accurate intuitive temperature feld distribution was obtained. The simulation results show: the cooling effect on the tunnel kiln cooling zone is better when the draft speed is 5 m/s, the pore rate is 0.5, and the bricks arrangement is sparse, and that is also agreeable with the tunnel kiln fring curve's cooling segment. This information provides a theoretical basis for the adjustment of the draft speed, pore rate and bricks arrangement to improve the quality of bricks fring and productivity.

Key words:tunnel kiln cooling zone; fuid - solid coupling; heat transfer; numerical simulation

基金項目:國家科技支撐項目水泥、墻材窯爐大氣污染治理關(guān)鍵技術(shù)研究與示范(2013BCA13B01);湖南省自然科學(xué)基金資助項目(13JJ9034);江蘇省新型環(huán)保重點實驗室資助項目(AE201109)。

收稿日期:2015-12-18。

修訂日期:2015-12-20。

DOI:10.13958/j.cnki.ztcg.2016.02.006

中圖分類號:TQ174.6+53

文獻(xiàn)標(biāo)志碼:A

文章編號:1006-2874(2016)02-0026-05

猜你喜歡
散熱
適用于車廂快速制冷優(yōu)化循環(huán)系統(tǒng)
一種小型電機(jī)定子扣片槽開口形式設(shè)計優(yōu)化
CPU故障的處理技巧
中文信息(2016年11期)2017-02-11 16:48:48
全密封式大功率硅整流柜的設(shè)計思考與應(yīng)用
商情(2016年43期)2016-12-26 00:00:00
功率芯片表面絕緣層厚度對石墨烯散熱效果的影響
IGBT模塊壓接技術(shù)研究
LED光源散熱不良的原因及對策
淺析FI分類號在服務(wù)器散熱領(lǐng)域中的應(yīng)用
列車車載LED燈結(jié)構(gòu)散熱分析
科技視界(2016年9期)2016-04-26 09:27:29
發(fā)動機(jī)活塞不確定缸的反復(fù)粘缸的故障排除
主站蜘蛛池模板: 国产精品久久久久鬼色| 欧美区在线播放| 91精品专区| 亚洲欧美一区二区三区图片| 亚洲男人天堂网址| 亚洲欧美在线综合一区二区三区| 亚洲午夜福利精品无码| 国产精品大白天新婚身材| 免费欧美一级| 久久精品丝袜| 久久精品国产免费观看频道| 五月婷婷综合色| 国产在线视频导航| 狠狠色噜噜狠狠狠狠奇米777 | 久久亚洲国产视频| 亚洲福利网址| 欧美成人精品高清在线下载| 久久久久国产精品嫩草影院| 久久人妻xunleige无码| 91在线播放国产| 国产又爽又黄无遮挡免费观看| 精品一區二區久久久久久久網站| 无码免费视频| 免费看的一级毛片| 亚洲国产清纯| 制服丝袜一区二区三区在线| 国产欧美日韩免费| 国产一级二级在线观看| 国产aⅴ无码专区亚洲av综合网| 欧美成人第一页| 中文字幕久久波多野结衣| 亚洲美女一级毛片| 日韩黄色在线| 98超碰在线观看| 精品国产www| 在线视频精品一区| 日韩在线永久免费播放| 国产内射一区亚洲| 伊人久久福利中文字幕| 激情五月婷婷综合网| 亚洲天堂首页| …亚洲 欧洲 另类 春色| 黄色网站在线观看无码| 欧洲亚洲一区| 国产亚洲美日韩AV中文字幕无码成人 | 国产91小视频在线观看| 亚洲精品国产综合99久久夜夜嗨| 欧美日本二区| 欧美一区二区福利视频| 91在线精品麻豆欧美在线| 中文国产成人精品久久| 日本欧美午夜| 久久免费视频6| 欧美在线视频不卡第一页| 亚洲精品国偷自产在线91正片| 国产高清又黄又嫩的免费视频网站| 色婷婷亚洲综合五月| 亚洲精品国产乱码不卡| 麻豆国产精品一二三在线观看| 亚洲无码精品在线播放| 成人国内精品久久久久影院| 欧洲欧美人成免费全部视频| 午夜国产精品视频黄| a级毛片免费网站| 一区二区午夜| 久久久久免费精品国产| 青草视频免费在线观看| 91精品啪在线观看国产91九色| 国产成人永久免费视频| 国产成人免费观看在线视频| 欧美人人干| 日韩在线欧美在线| 99九九成人免费视频精品| 欧美亚洲激情| 欧美日本在线| 精品综合久久久久久97| 欧美视频免费一区二区三区| 一区二区三区四区在线| 97久久精品人人| …亚洲 欧洲 另类 春色| 久久亚洲精少妇毛片午夜无码| 久久一日本道色综合久久|