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

        ?

        基于CFD技術(shù)的瓶裝啤酒隧道式巴氏殺菌過程研究

        2016-09-28 01:27:46菊安家彥董文勇王
        食品與機(jī)械 2016年8期
        關(guān)鍵詞:啤酒瓶噴口巴氏

        程 菊安家彥董文勇王 越

        (1.大連工業(yè)大學(xué)生物工程學(xué)院,遼寧 大連 116034;2. 大連市第21中學(xué),遼寧 大連 116021)

        ?

        基于CFD技術(shù)的瓶裝啤酒隧道式巴氏殺菌過程研究

        程菊1安家彥1董文勇2王越1

        (1.大連工業(yè)大學(xué)生物工程學(xué)院,遼寧 大連116034;2. 大連市第21中學(xué),遼寧 大連116021)

        通過數(shù)值模擬瓶裝啤酒隧道式巴氏殺菌的過熱階段(巴氏殺菌加熱部分的最后一個(gè)階段),確定瓶中的冷核(slowest heating zone,SHZ),并運(yùn)用L9(34)正交試驗(yàn),分析噴淋水溫度、噴口處噴淋水湍流強(qiáng)度以及瓶子的運(yùn)行速度三因素對(duì)瓶中啤酒溫度分布均勻性的影響。結(jié)果表明:冷核起先是位于瓶底,后向上移動(dòng),但不超過灌裝高度的1/2;在加熱至7 min時(shí),噴淋水溫度、噴口處噴淋水湍流強(qiáng)度影響顯著,瓶子的運(yùn)行速度則不顯著,試驗(yàn)的可能的最優(yōu)條件:噴淋水溫度為65 ℃,噴口處噴淋水湍流強(qiáng)度為3%,瓶子運(yùn)行的速度為3 mm/s。

        啤酒;巴氏殺菌;數(shù)值模擬;冷核

        隧道式巴氏殺菌是啤酒工業(yè)化生產(chǎn)過程中,殺死啤酒中的微生物以提高其生物穩(wěn)定性的最常用的方法之一[1]。隧道式巴氏殺菌機(jī)包括加熱區(qū)(包括預(yù)熱區(qū)和過熱區(qū))、保溫區(qū)和冷卻區(qū)三大溫區(qū)[2-3]。在加熱區(qū),啤酒溫度大約由2 ℃上升至60 ℃;在保溫區(qū),維持在60 ℃;經(jīng)冷卻區(qū)處理后,產(chǎn)品的溫度由60 ℃降至21~27 ℃[4]。瓶裝啤酒在巴氏殺菌過程中,啤酒以酒瓶為界面進(jìn)行傳熱:靠近瓶壁的啤酒溫度升高,密度減小,而中心部分溫度相對(duì)較低,密度相對(duì)較大,產(chǎn)生自然對(duì)流[5-6]。這就形成了冷核。冷核是在食物產(chǎn)品熱處理過程中,吸收熱量最少(或殺菌強(qiáng)度最弱)的區(qū)域[7-8]。因此,在啤酒巴氏殺菌的過程中,既要保證冷核達(dá)到滅菌要求,又要使瓶中除冷核以外的區(qū)域不會(huì)過巴氏殺菌,以免改變啤酒的感官,造成啤酒品質(zhì)的降低。這就要求瓶中啤酒溫度分布要均勻。

        目前研究瓶中啤酒的溫度分布及冷核的方法有兩種:① 在線檢測(cè)瓶中溫度和數(shù)值模擬啤酒巴氏殺菌過程。在線監(jiān)測(cè)過程中,溫度探針會(huì)影響瓶中液體的流動(dòng)和溫度分布[9-10],同時(shí)只能檢測(cè)一個(gè)位置,而冷核的溫度和位置會(huì)隨時(shí)間和包裝容器結(jié)構(gòu)變化而改變[11-12],且在線監(jiān)測(cè)的費(fèi)用較高,操作復(fù)雜。② 借助計(jì)算流體力學(xué)(computational fluid dynamics,CFD)技術(shù),數(shù)值模擬啤酒巴氏殺菌過程,從而得到包裝容器內(nèi)液體的傳熱規(guī)律,可有效控制食品的品質(zhì)。而且CFD能縮短時(shí)間,減少花費(fèi)[13-14]。目前,CFD模擬已被用于研究啤酒的巴氏殺菌過程,如:優(yōu)化啤酒巴氏殺菌機(jī)的隧道結(jié)構(gòu)[15];啤酒瓶方位的變化對(duì)巴氏殺菌單位(pasteurization unit,PU)的影響[2];瓶內(nèi)對(duì)流流動(dòng)對(duì)巴氏殺菌效果和老化效果的影響[3];瓶裝和罐裝對(duì)包裝內(nèi)啤酒溫度分布、速度分布和PU的影響[16];啤酒巴氏殺菌結(jié)束后,啤酒的出口溫度和PU值是否達(dá)標(biāo)[17];隨著加熱的進(jìn)行瓶內(nèi)溫度分布和SHZ位置的變化以及啤酒中酵母的死滅程度(用PU表示)[8]等。

        在巴氏殺菌過程中,瓶內(nèi)啤酒的冷核和溫度分布受加熱時(shí)間、瓶子間距、瓶子移動(dòng)的速度、噴淋水的溫度和噴口處噴淋水湍流強(qiáng)度的影響[18-19]。目前,已有關(guān)于瓶中啤酒冷核隨時(shí)間變化[8]及瓶中啤酒加熱速率隨噴淋水的噴淋強(qiáng)度變化的研究[19],但未見瓶子移動(dòng)的速度、噴淋水的溫度和噴口處噴淋水湍流強(qiáng)度對(duì)瓶中啤酒溫度分布的影響的研究。因此本研究擬考察隧道式巴氏殺菌的過熱階段,運(yùn)用L9(34)正交試驗(yàn),分析噴淋水溫度、噴口處噴淋水湍流強(qiáng)度以及瓶子的運(yùn)行速度三因素對(duì)瓶中啤酒溫度分布均勻性的影響,為瓶裝啤酒巴氏殺菌提供理論基礎(chǔ)。

        1 材料與方法

        1.1CFD模擬的驗(yàn)證方法

        選取高274 mm,直徑72 mm,容積為500 mL的啤酒瓶為研究對(duì)象,其模型見圖1。試驗(yàn)所用的巴氏殺菌機(jī)型號(hào)為SBW40A(廣東輕工機(jī)械二廠有限公司)。壓好蓋的瓶裝生啤酒通過巴氏殺菌機(jī)的隧道時(shí),用ZH5793 PU計(jì)(北京中慧天誠科技有限公司)分別測(cè)定隧道中不同位置上的瓶中啤酒的溫度(軸線上灌裝高度1/3處溫度),計(jì)算過熱階段瓶中啤酒不同時(shí)間酒溫的平均值。

        圖1 啤酒瓶模型

        1.2幾何模型

        因計(jì)算機(jī)的計(jì)算能力有限,動(dòng)網(wǎng)格的計(jì)算復(fù)雜,在三維條件下,整個(gè)模型很難進(jìn)行計(jì)算,而在二維條件下既能節(jié)約計(jì)算時(shí)間,又能節(jié)省人力和物力。所以整個(gè)模型是二維的。

        巴氏殺菌機(jī)模型網(wǎng)格的建立使用前處理軟件Gambit 2.3.16,模型的網(wǎng)格數(shù)為209 319,流體計(jì)算使用Fluent 12.0。噴淋口與瓶口的距離為300 mm,啤酒瓶個(gè)數(shù)為142個(gè)(設(shè)置了邊界條件的為3個(gè)),瓶子間距為16 mm。二維網(wǎng)格模型的部分示意圖見圖2。

        1.3控制方程(governing equation)

        運(yùn)用Ansys Fluent 12(2009)軟件求解模型及相關(guān)邊界條件的連續(xù)性方程、動(dòng)量方程、能量方程[20]。

        圖2 二維網(wǎng)格模型的部分示意圖

        (a) 連續(xù)性方程(continuity equation)或質(zhì)量守恒方程(conversation of mass equation):

        (1)

        式中:

        ρ——密度,kg/m3;

        t——時(shí)間,s;

        x——軸向坐標(biāo);

        u——軸向速度,m/s;

        v——徑向速度,m/s;

        r——徑向坐標(biāo);

        Sm——源項(xiàng),指加入到連續(xù)相的質(zhì)量(例如液滴的蒸發(fā)),也指其他自定義源項(xiàng),kg。

        (b) 動(dòng)量守恒方程:

        軸向的動(dòng)量守恒方程:

        (2)

        徑向的動(dòng)量守恒方程:

        (3)

        其中:

        (4)

        式中:

        P——靜壓,Pa;

        μ——黏度,kg/(m·s);

        Fx和Fr——外部體積力,N。

        (c) 能量守恒方程:

        (5)

        其中:

        (6)

        (7)

        式中:

        τij——應(yīng)力張量;

        Sh——化學(xué)反應(yīng)熱和其他體積熱,J;

        keff——有效導(dǎo)熱系數(shù),W/(m·℃);

        Jj’——組分j’擴(kuò)散通量,mol/(m2·s)。

        1.4模擬條件

        啤酒、玻璃啤酒瓶、空氣和水的熱物理性質(zhì)列于表1[8, 16, 21-23]。黏度模型采用k-epsilon中RNG模型,噴嘴及其附近的長方形區(qū)域?yàn)閯?dòng)網(wǎng)格,運(yùn)動(dòng)速度大小為5 mm/s。瓶內(nèi)的初始溫度為48 ℃,瓶內(nèi)壓力為1.55×105Pa。壓力和速度的耦合計(jì)算采用SIMPLEC算法。

        1.5L9(34)正交試驗(yàn)

        瓶中啤酒溫度分布均勻性和冷核受噴淋水溫度、噴口處噴淋水湍流強(qiáng)度和瓶子運(yùn)行速度的影響,因此確定四因素三水平的正交試驗(yàn),其試驗(yàn)方案見表2。

        表1 空氣、啤酒、玻璃和水的熱物理性質(zhì)

        表2 L9(34)正交試驗(yàn)方案

        2 結(jié)果與討論

        2.1CFD模擬的驗(yàn)證

        為了確保CFD模擬的可靠性,先對(duì)其進(jìn)行實(shí)驗(yàn)驗(yàn)證。啤酒在通過巴氏殺菌過熱區(qū)(工藝參數(shù):噴淋水溫度為61 ℃,初始溫度48 ℃,加熱9 min)時(shí),將啤酒瓶軸線上灌裝高度1/3位置處的溫度實(shí)驗(yàn)值和模擬值進(jìn)行對(duì)比,結(jié)果見圖3。由圖3可知,模擬值與實(shí)驗(yàn)值的溫度曲線走勢(shì)相同,它們之間的最大標(biāo)準(zhǔn)差為0.21,所以可以認(rèn)為它們是吻合的??梢娔M所設(shè)置各項(xiàng)參數(shù)是合理的,可以進(jìn)行接下來的數(shù)值模擬。

        2.2冷核的確定

        為了考察瓶內(nèi)啤酒在過熱區(qū)的溫度分布情況,設(shè)過熱區(qū)噴淋水的溫度為61 ℃,噴口處噴淋水湍流強(qiáng)度為1%,瓶子運(yùn)行速度為3 mm/s,按照方法1.2和1.3進(jìn)行數(shù)值模擬,結(jié)果見圖4。

        圖3 模擬溫度與實(shí)驗(yàn)溫度曲線的比較

        圖4 瓶中啤酒溫度分布云圖

        由圖4可知,1,5,9 min時(shí)瓶中啤酒有明顯溫度分布梯度,具有冷核,且冷核位置是隨著時(shí)間在變化的。1 min時(shí)冷核位于瓶底,這與傳統(tǒng)的理論相符;但隨著時(shí)間的進(jìn)行,冷核位置在不斷的變化,到5 min時(shí),冷核的位置位于啤酒灌裝高度1/4~1/2。瓶內(nèi)啤酒在加熱過程中由于局部溫度、密度等的差異,引起啤酒的自然對(duì)流,從而影響傳熱效果,溫度分布不均勻,因此有冷核存在。靠近瓶壁的啤酒先受熱,溫度較瓶中心的啤酒溫度更高,所以冷核一直位于啤酒瓶軸線附近的區(qū)域。而且隨著加熱的進(jìn)行,瓶底也慢慢被加熱,靠近瓶底部的啤酒溫度也慢慢地上升,所以冷核從瓶底向瓶底以上部分移動(dòng)。

        2.3L9(34)正交試驗(yàn)結(jié)果

        為了考察瓶裝啤酒通過巴氏殺菌機(jī)的過熱區(qū)時(shí),巴氏殺菌參數(shù)—噴淋水溫度、噴口處噴淋水湍流強(qiáng)度以及瓶子的運(yùn)行速度,對(duì)瓶中啤酒溫度分布均一性的影響,利用L9(34)正交試驗(yàn),數(shù)值模擬啤酒巴氏殺菌的過熱階段,分析啤酒冷核面積占啤酒總面積的比例。由于隨著加熱時(shí)間的延長,冷核的溫度在不斷升高,根據(jù)9組試驗(yàn)的不同時(shí)間的冷核溫度,將 1,3,5,7,9 min時(shí)冷核溫度的取值范圍分別設(shè)為小于:50.42,54.54,57.03,58.57,59.51 ℃。分別計(jì)算9組試驗(yàn)在1,3,5,7,9 min時(shí)啤酒冷核面積占啤酒總面積的比例,結(jié)果列于表3。按照所設(shè)的冷核溫度范圍,對(duì)所得的試驗(yàn)數(shù)據(jù)進(jìn)行方差分析,發(fā)現(xiàn)1,3,5,9 min時(shí)噴淋水溫度、噴口處噴淋水湍流強(qiáng)度以及瓶子的運(yùn)行速度對(duì)瓶中啤酒溫度分布沒有影響或影響不顯著。這是因?yàn)樵诩訜岬?~5 min時(shí),瓶中啤酒加熱時(shí)間短,溫度上升不顯著;加熱9 min時(shí),瓶中啤酒溫度分布已經(jīng)均勻。7 min時(shí)試驗(yàn)數(shù)據(jù)的方差分析及極差分析結(jié)果分別見表4、5。

        經(jīng)過計(jì)算及方差分析,取檢驗(yàn)水平α=0.05,由表4可知,噴淋水溫度和噴口處噴淋水湍流強(qiáng)度對(duì)瓶中啤酒溫度分布均一性的影響顯著,啤酒瓶運(yùn)行速度則影響不顯著。由表5可知,三因素的影響主次順序?yàn)椋築>A>C,即噴口處噴淋水湍流強(qiáng)度>噴淋水溫度>啤酒瓶運(yùn)行速度,結(jié)果與方差分析的相符,而可能的最優(yōu)條件組合為A3B3C1,即噴淋水溫度為65 ℃,噴口處噴淋水強(qiáng)度為3%,瓶子運(yùn)行的速度為3 mm/s。

        表3    二維啤酒瓶模型中的啤酒冷核面積占啤酒總面積的比例

        ?*表示影響顯著。

        表5 7 min時(shí)的極差分析

        3 結(jié)論

        本研究利用CFD技術(shù)模擬了500 mL瓶裝啤酒的隧道式巴氏殺菌過程的過熱階段,研究噴淋水溫度、噴淋口噴淋水湍流強(qiáng)度和瓶子運(yùn)行速度對(duì)瓶中啤酒溫度分布的影響,結(jié)論如下:

        (1) 運(yùn)用商業(yè)軟件Ansys Fluent12模擬的巴氏殺菌過熱階段溫度曲線,得到了實(shí)驗(yàn)溫度曲線的驗(yàn)證,說明該模擬方法和模擬參數(shù)設(shè)置是合理的。

        (2) 在加熱過程中,瓶中啤酒存在冷核,表明瓶中啤酒的溫度分布不均勻;冷核的位置隨著加熱的進(jìn)行向上移動(dòng),但不超過灌裝高度的1/2。

        (3) 利用L9(34)正交試驗(yàn)研究啤酒巴氏殺菌的過熱階段影響因素結(jié)果表明,噴口處噴淋水湍流強(qiáng)度和噴淋水溫度對(duì)瓶中啤酒溫度分布均一性的影響顯著,瓶子的移動(dòng)速度則影響不明顯。

        [1] BUZRUL S. A suitable model of microbial survival curves for beer pasteurization[J]. LWT - Food Science and Technology, 2007, 40(8): 1 330-1 336.

        [2] AUGUSTO P E D, PINHEIRO T F, CRISTIANINI M. Using computational fluid-dynamics (CFD) for the evaluation of beer pasteurization: effect of orientation of cans[J]. Food Science and Technology, 2010, 30(4): 980-986.

        [3] HORN C S, FRANKE M, BLAKEMORE F B, et al. Modelling and simulation of pasteurization and staling effects during tunnel pasteurization of bottled beer[J]. Food and Bioproducts Processing, 1997, 75(1): 23-33.

        [4] ENGELMAN M S, SANI R L. Finite-element simulation of an in-package pasteurization process[J]. Numerical Heat Transfer, 1983, 6(1): 41-54.

        [5] TEIXEIRA A A, DATTA A K. Numerically predicted transient temperature and velocity profiles during natural convection heating of canned liquid foods[J]. Journal of Food Science, 1988, 53(1): 191-195.

        [6] ZECHMAN L G. Location of the slowest heating zone for natural convection heating fluids in metal containers[J]. Journal of Food Science, 1989, 54(1): 205-209.

        [7] KUMAR A, BHATTACHARYA M, BLAYLOCK J. Numerical simulation of natural convection heating of canned thick viscous liquid food products[J]. Journal of Food Engineering, 1990, 55(5): 1 403-1 411.

        [8] BHUVANESWARI E, ANANDHARAMAKRISHNAN C. Heat transfer analysis of pasteurization of bottled beer in a tunnel pasteurizer using computational fluid dynamics[J]. Innovative Food Science & Emerging Technologies, 2014, 23: 156-163.

        [9] KANNAN A, SANDAKA P C G. Heat transfer analysis of canned food sterilization in a still retort[J]. Journal of Food Engineering, 2008, 88(2): 213-228.

        [10] FARID M, ABDUL GHANI A G. A new computational technique for the estimation of sterilization time in canned food[J]. Chemical Engineering and Processing: Process Intensification, 2004, 43(4): 523-531.

        [11] GHANI A G A, FARID M M, CHEN X D. Theoretical and experimental investigation of the thermal inactivation of bacillus stearothermophilus in food pouches[J]. Journal of Food Engineering, 2002, 51(3): 221-228.

        [12] VARMA M N, KANNAN A. CFD studies on natural convective heating of canned food in conical and cylindrical containers[J]. Journal of Food Engineering, 2006, 77(4): 1024-1036.

        [13] ASLAM BHUTTA M M, HAYAT N, BASHIR M H, et al. CFD applications in various heat exchangers design: a review[J]. Applied Thermal Engineering, 2012, 32: 1-12.

        [14] NORTON T, SUN Da-wen. Computational fluid dynamics (CFD)-an effective and efficient design and analysis tool for the food industry: a review[J]. Trends in Food Science & Technology, 2006, 17(11): 600-620.

        [15] DILAY E, VARGAS J V C, AMICO S C, et al. Modeling, simulation and optimization of a beer pasteurization tunnel[J]. Journal of Food Engineering, 2006, 77(3): 500-513.

        [16] 王亮, 于艷艷, 馬曉彬, 等.罐裝啤酒與瓶裝啤酒的巴氏殺菌過程數(shù)值模擬[J]. 食品與發(fā)酵工業(yè), 2014, 40(4): 42-46.

        [17] 洪曉敏, 許薔. 啤酒巴氏殺菌過程的CFD數(shù)值分析[J]. 食品與機(jī)械, 2016, 32(2): 160-164.

        [18] LEWICKI P P, WALCZAK W, GOSS B. Heat transfer in a tunnel pasteuriser. part i. factors affecting the rate of heating of liquid in a bottle[J]. Journal of Food Engineering, 1983, 2(4): 309-322.

        [19] 樊延敏, 王長偉. 淺談啤酒殺菌機(jī)PU值控制系統(tǒng)[J]. 啤酒科技, 2003(7): 36, 38.

        [20] 張凱, 王瑞金, 王剛. Fluent技術(shù)基礎(chǔ)與應(yīng)用實(shí)例[M]. 2版. 北京: 清華大學(xué)出版, 2010: 33-36.

        [21] AUGUSTO P E D, CRISTIANINI M. Numerical simulation of packed liquid food thermal process using computational fluid dynamics(CFD)[J]. International Journal of Food Engineering, 2011, 7(4): 1-22.

        [22] BAILEY R T, ELBAN W L. Thermal performance of aluminum and glass beer bottles[J]. Heat Transfer Engineering, 2008, 29(7): 643-650.

        [23] ERDOGDU F, TUTAR M. A computational study for axial rotation effects on heat transfer in rotating cans containing liquid water, semi-fluid food system and headspace[J]. International Journal of Heat and Mass Transfer, 2012, 55(13/14): 3 774-3 788.

        Study on tunnel pasteurization of bottled beer based on CFD technology

        CHENG Ju1ANJia-yan1DONGWeng-yong2WANGYue1

        (1.SchoolofBiologicalEngineering,DalianPolytechnicUniversity,Dalian,Liaoning116034,China;2.No.12MiddleSchoolofDalian,Dalian,Liaoning116021,China)

        In the industrial processing of beer, pasteurization is necessary the to guarantee its biological stability. The superheat phase (the last part of heating phase) was numerically simulated, in which the slowest heating zone (SHZ) inside the bottle was determined, and then the influence of three parameters on the uniformity of temperature distribution inside the bottled beer during pasteurization, i.e. the temperature of spray water, the turbulent intensity of spray water at the spray nozzles and the running velocity of bottles were analyzed using L9(34) orthogonal experiments. The results indicated that SHZ lay at the bottom of the bottle at the beginning and then moved upward without exceeding 1/2 of filling height of beer. Besides, after heating for 7 min, the temperature and the turbulent intensity of spray water at the spray nozzles affected the uniformity of temperature distribution inside the bottled beer obviously, while the running velocity of bottle had no significant effect on it. Our results showed that the probable optimum condition was 65 ℃ temperature of spray water and 3% turbulent intensity of it at the spray nozzles as well as 3 mm/s running velocity of bottle.

        beer; pasteurization; numerical simulation; slowest heating zones (SHZ)

        程菊,女,大連工業(yè)大學(xué)在讀研究生。

        王越(1972—),女,大連工業(yè)大學(xué)副教授,博士。

        E-mail:wydqs2004@163.com

        2016—02—16

        10.13652/j.issn.1003-5788.2016.08.024

        猜你喜歡
        啤酒瓶噴口巴氏
        噴口形狀對(duì)噴水推進(jìn)器性能的影響
        飛參數(shù)據(jù)在某型飛機(jī)噴口工作監(jiān)控中的應(yīng)用
        釋放巴氏新小綏螨可滿足對(duì)蘋果全爪螨的防治需求
        創(chuàng)意啤酒瓶蓋
        快樂語文(2017年15期)2017-08-11 13:06:10
        小蜘蛛
        高速?zèng)_床噴口壓力及位置對(duì)工件出模影響分析
        惱人的啤酒瓶蓋
        巴氏殺菌水牛奶在不同儲(chǔ)存條件下微生物增長規(guī)律的研究
        惱人的啤酒瓶蓋
        換啤酒
        内射中出无码护士在线| 日韩美女av一区二区| 欧美国产激情18| 日韩一卡2卡3卡4卡新区亚洲| 99精品国产在热久久国产乱| 亚洲av中文无码乱人伦在线咪咕| 亚洲精品中文字幕一二三 | 日本一区二区三区在线观看免费 | 中文字幕av素人专区| 亚洲av日韩av女同同性| 在线va免费看成| 最新国产成人综合在线观看| 人妻系列中文字幕av| 国产av国片精品jk制服| 亚洲av日韩av无码av| 99久久精品国产片| 青青草中文字幕在线播放| 日韩精品久久久久久免费| 精品少妇一区二区三区视频| 日韩女优中文字幕在线| 国产亚洲成人精品久久久| 国语对白做受xxxxx在| 国产人成午夜免电影观看| 国内自拍视频在线观看| 男女男精品视频网站免费看| 免费a级毛片无码a∨男男| 产国语一级特黄aa大片| 美女性色av一区二区三区| 国产精品久久久久一区二区三区 | 男人的天堂无码动漫av| 欧美中文字幕在线| av天堂手机一区在线| 国产在线无码一区二区三区视频| 欧美日韩视频无码一区二区三| 囯产精品无码一区二区三区AV| 国产在线一区二区三区四区乱码| 又色又爽又高潮免费视频国产| 国产成人午夜福利在线观看者| 久久婷婷色香五月综合激激情| 变态调教一区二区三区女同| 国产精品久久久|