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

        ?

        射流動力學邊界條件的研究*

        2019-11-05 02:05:38曹建明王德超張凱妹
        新能源進展 2019年5期
        關鍵詞:瑞利液膜泰勒

        曹建明,王德超,舒 力,張凱妹

        (長安大學 汽車學院,西安 710064)

        0 前 言

        線性穩(wěn)定性理論僅研究第一級波,因此僅應用第一級波的動力學邊界條件[1]。但線性穩(wěn)定性理論能夠研究三維擾動的瑞利-泰勒波,因此本文將推導出三維擾動的第一級瑞利波和瑞利-泰勒波的線性動力學邊界條件。非線性穩(wěn)定性理論能夠研究三級波,因此要應用第一級波、第二級波和第三級波的動力學邊界條件[2]。但非線性穩(wěn)定性理論僅研究一維擾動的瑞利波,因此本文將推導出一維擾動三級瑞利波的非線性動力學邊界條件。

        目前,對于線性穩(wěn)定性理論,平面液膜、圓射流和環(huán)狀液膜二維擾動的第一級瑞利波的動力學邊界條件以及平面液膜二維擾動的第一級瑞利-泰勒波的動力學邊界條件已有定論[3-18]。對于非線性穩(wěn)定性理論,平面液膜一維擾動三級瑞利波的動力學邊界條件和圓射流在特殊條件下(如:噴嘴出口液體流速為零,韋伯數(shù)等于 1)一維擾動三級瑞利波的動力學邊界條件已有研究[19-20]。上述研究鮮有對于詳細推導過程的論述。本文的推導結果與上述研究結果完全一致,創(chuàng)新點在于:從動力學邊界條件的普適性原始表達式入手,經(jīng)過細致全面的推導,將得到平面液膜、圓射流和環(huán)狀液膜三維擾動瑞利波的線性動力學邊界條件,以及平面液膜、圓射流和環(huán)狀液膜三維擾動瑞利-泰勒波的線性動力學邊界條件,得到圓射流和環(huán)狀液膜一維擾動三級瑞利波的非線性動力學邊界條件。給出各種動力學邊界條件的簡化和適用情況,并將全部的動力學邊界條件全面、系統(tǒng)地相互關聯(lián)在一起,為進行射流碎裂過程研究的學者提供各種物理模型下明確的動力學邊界條件表達式。

        動力學邊界條件的通式由楊-拉普拉斯方程確定[3]。

        1 液相與氣相的應力張量差

        1.1 平面液膜的應力張量差

        對于沿x、y、z方向的三維擾動,有量綱形式的液相與氣相的應力張量差為:

        式中:p為擾動壓力,Pa;μ為動力學粘度系數(shù),Pa?s;u、v、w為沿x、y、z方向的擾動速度,m/s;為直角坐標系的哈密頓算子,其中,分別為沿x、y、z方向的單位矢量。角標“l(fā)”表示液相,“g”表示氣相;“j= 1”表示平面液膜的上氣液交界面,“j= 2”表示下氣液交界面。

        由于平面液膜很薄,三維擾動可以簡化為一維擾動。一維擾動的線性穩(wěn)定性理論應力張量差為方程(7)。由于液相有粘性,因此線性穩(wěn)定性理論的液相壓力pl和氣相壓力pgj要通過連續(xù)性方程和納維-斯托克斯方程組(Navier-Stokes equations)求得。對于非線性穩(wěn)定性理論,由于假設液相和氣相均為無粘性、不可壓縮的一維擾動理想流體,因此,方程(7)可以簡化為:

        非線性穩(wěn)定性理論的液相壓力pl和氣相壓力pgj可以通過伯努利方程直接求得。

        1.2 圓射流和環(huán)狀液膜的應力張量差

        對于沿r、θ、z方向的三維擾動,有量綱形式的液相與氣相的應力張量差為:

        對于圓射流和環(huán)狀液膜,三維擾動也可以簡化為一維擾動。一維擾動線性穩(wěn)定性理論的應力張量差為方程(12)。圓射流和環(huán)狀液膜一維擾動非線性穩(wěn)定性理論的應力張量差與平面液膜的同式,為方程(8)。

        2 由表面張力引起的附加壓強

        2.1 附加壓強

        如圖1所示,曲率半徑R1和R2分別表示平面液膜氣液交界曲面上沿x方向和y方向的曲率半徑,或者分別表示圓射流和環(huán)狀液膜氣液交界曲面上沿z方向和r方向的曲率半徑。曲面DEFG為氣液交界面,過O點做曲面DEFG的法線N,沿法線N做兩個相互垂直的平面A1B1C1和A2B2C2,兩平面與曲面DEFG相交于兩條相互垂直的曲線A1B1和A2B2,曲線A1B1和A2B2的曲率半徑分別為R1=OC1和R2=OC2。R3表示圓射流和環(huán)狀液膜氣液交界曲面DEFG沿z軸扭轉一個角度θ后,到達曲面D'E'F'G'位置的曲率半徑,R3=O'C2(O'與B2在圖中重合),如圖2所示。

        如圖3所示,當曲率中心在液體內部時,氣體擾動壓力要大于液體擾動壓力,將氣液交界面壓向液體一側彎曲,那么方程(1)中等號右側應為負號,即由表面張力引起的附加壓強為負;而當曲率中心在液體外部時,液體擾動壓力要大于氣體擾動壓力,將氣液交界面壓向氣體一側彎曲,方程(1)中等號右側應為正號,由表面張力引起的附加壓強為正[4]。因此,如果沿某一方向的氣液交界面有波動,那么附加壓強為正負號。如果沿某一方向的氣液交界面沒有波動,而是一條直線,則附加壓強為零;如果沿某一方向的氣液交界面沒有波動,而是一個圓,那么附加壓強恒為負。

        圖1 平面液膜氣液交界曲面上的曲率半徑Fig.1 Curvature radius of liquid and gas interfaces of a plane liquid film

        圖2 圓射流和環(huán)狀液膜氣液交界曲面上的曲率半徑Fig.2 Curvature radius of liquid and gas interfaces of a liquid jet and an annular liquid film

        圖 3 表面波氣液交界面上由表面張力引起的附加壓強正負號的選取Fig.3 Plus and minus of addition pressure caused by surface tension of liquid and gas interfaces of the surface wave

        2.2 曲率半徑和曲率

        2.2.1 瑞利波

        瑞利波是沿液體噴射方向的單向波,曲率半徑和曲率指的是氣液交界面表面波沿射流噴射方向曲線的曲率半徑和曲率。對于直角坐標系的平面液膜,沿x方向曲面上任意一點的y方向氣液交界線為一條直線,即沿y方向沒有波動。因此,瑞利波的曲率半徑和曲率僅與x坐標有關,而與y坐標無關。對于圓柱坐標系的圓射流和環(huán)狀液膜,沿z方向曲面上任意一點的θ方向氣液交界線為一個圓,即沿θ方向沒有波動。但是,z方向曲線的曲率半徑和曲率要受到r方向和θ方向曲率半徑和曲率的影響。

        (1)平面液膜瑞利波的曲率半徑和曲率

        如圖1所示,平行于x軸的平面A1B1C1的曲率半徑即為沿x方向的曲率半徑R1=OC1。曲率為:

        式中:角標“R”表示瑞利波,后文角標“T”表示泰勒波,角標“R-T”表示瑞利-泰勒聯(lián)合波。Δα為弧線MM1上兩端點M和M1切線的夾角,當M→M1時,Δα →dα,見圖4a。Δs為弧線MM1的長度,當M→M1時,Δs →ds,見圖4b。

        式中:ξ為射流氣液交界面表面波的擾動振幅,m。

        由于表面波既有波峰又有波谷,是波動的。因此KR1有正負號。

        式中:j= 1, 2。其中:(- 1)j中的“j= 1”表示波峰位置,“j= 2”表示波谷位置。ξj,xx中的“j= 1”表示平面液膜的上氣液交界面,“j= 2”表示下氣液交界面。

        如圖1所示,平行于y軸的平面A2B2C2的曲率半徑即為沿y方向的曲率半徑R2=OC2。當僅有瑞利波,而無泰勒波時,沿x方向曲面上任意一點的y方向氣液交界線均為一條直線。由于直線的曲率半徑R2= ∞,曲率則平面液膜非線性的瑞利波曲率為:

        該式即為一條波動曲線的曲率方程。

        忽略非線性的ξ的一階偏導數(shù)項,對方程(20)線性化,得平面液膜線性化的瑞利波曲率。

        圖4 平面液膜沿x方向的曲率半徑Fig.4 x direction curvature radius of a plane liquid film

        (2)圓射流和環(huán)狀液膜瑞利波的曲率半徑和曲率

        如圖2所示,R1、R2和R3分別表示圓射流和環(huán)狀液膜沿z方向、r方向和θ方向的曲率半徑。沿z方向平面A1B1C1的曲率半徑為R1=OC1。由于A1B1為一條波動的曲線,因此有正負號。

        如圖2和圖 5所示。沿平行于r方向的平面A2B2C2的曲率半徑為R2=OC2。將弧線A2B2的曲率中心C2置于z軸底平面上,過O點做垂直于z軸、平行于r軸的輔助直線OB。過O點做輔助直角三角形ΔOMD、OC2⊥OM、ND⊥OB,則∠OMD=∠BOC2,三角形 ΔOMD與 ΔBOC2為相似三角形。當Δz→0時用于非線性穩(wěn)定性分析;當此可以近似看作是小擾動,可用于線性穩(wěn)定性分析。

        圖5 圓射流和環(huán)狀液膜瑞利波的曲率半徑R2Fig.5 R wave curvature radius R2 of a liquid jet and an annular liquid film

        當僅有瑞利波,而無泰勒波時,圍繞z軸是個圓面,則其附加壓強恒為負號。由于瑞利波的波動,圓面的半徑R2要隨著z方向的位移而變化。因此,對于線性穩(wěn)定性分析:

        對于非線性穩(wěn)定性分析:

        如圖2所示,當曲面DEFG繞z軸扭轉一個角度θ后,扭轉到曲面D'E'F'G'位置。由于表面張力的作用,曲面DEFG扭轉到曲面D'E'F'G'后會有變形,O點除了會發(fā)生沿z軸的旋轉之外,還會降低高度,即擾動振幅ξ變小,并沿z軸正向有一個微小位移。因此,曲率半徑R3相對于R2會有所變化。因此,與旋轉角度有關的曲率為:

        與旋轉角度有關的曲率KR3決定了圓射流的扭轉。

        瑞利波的曲率為方程(22)、(23)、(25)之和,即:

        忽略掉非線性的ξ一階導數(shù)項,對方程(26)進行線性化,得:

        方程(27)即為線性穩(wěn)定性分析瑞利波的曲率方程。

        對于環(huán)狀液膜,射流雖然也具有軸對稱特點,但由于液膜很薄,一旦發(fā)生扭轉,液膜會立即碎裂。因此,環(huán)狀液膜不存在曲率半徑R3和曲率KR3,在方程(27)中將曲率KR3項直接刪去即可。即

        瑞利波非線性穩(wěn)定性分析的曲率為方程(22)、(24)、(25)之和,即

        對于非線性穩(wěn)定性理論,由于圓射流僅研究常見的n= 0的正對稱波形和n= 1的反對稱波形,為單股狀,并未發(fā)生扭轉。環(huán)狀液膜一旦發(fā)生扭轉,液膜會立即碎裂。因此曲率KR3項應刪去,方程(29)變成

        方程(30)即為圓射流和環(huán)狀液膜非線性穩(wěn)定性分析瑞利波正反對稱波形的曲率方程。

        2.2.2 瑞利-泰勒波

        當既有瑞利波,又有泰勒波時,波動為沿射流噴射方向和垂直方向的雙向波,曲率半徑和曲率指的是氣液交界面表面波曲面的曲率半徑和曲率。對于直角坐標系的平面液膜,沿x方向曲面上任意一點的y方向氣液交界線不再為一條直線,而是一條波動的曲線,即沿x方向和y方向均有波動。因此,瑞利-泰勒波的曲率半徑和曲率與x坐標和y坐標均有關。對于圓柱坐標系的圓射流和環(huán)狀液膜,沿z方向曲面上任意一點的θ方向氣液交界線不再為一個圓,而是一條波動的曲線,即沿z方向和θ方向均有波動。

        從數(shù)學的微分幾何角度來說,過任意曲面上的某一點上具有無窮多個相互正交的平面,這些平面與曲面相交產(chǎn)生無窮多個相互正交的曲線,其中總會存在一條使得曲率K1成為極大值的曲線,與之正交的曲線的曲率K2為極小值。則K1和K2稱為該曲面的主曲率。曲面上兩個主曲率之和的平均值稱為該曲面的平均曲率(又稱為中曲率),即兩個主曲率之積稱為高斯曲率(又稱為總曲率或全曲率),即K=K1K2[5]。該理論是由哥丁根大學的法國數(shù)學家Sophie Germain于1831年在她的著作《Theory of Elasticity》中最早提出的[6]。在流體力學中,采用平均曲率作為曲面曲率,并且將平均曲率中的2舍去。因此,流體力學中的曲面曲率為兩個主曲率之和,即K=K1+K2[5-6]。

        泰勒波的振幅非常小,其曲率半徑趨近于無窮大,曲率則趨近于零。可以認為其曲率為極小值,與其垂直的瑞利波的曲率則為極大值。瑞利波曲線曲率與泰勒波曲線曲率之和就是瑞利-泰勒波曲面的曲率。因此,瑞利-泰勒波曲面的主曲率就是KR和KT。在對瑞利-泰勒波色散關系式的研究中發(fā)現(xiàn),當忽略泰勒波時,瑞利波仍然存在,而當忽略瑞利波時,泰勒波將不復存在。說明瑞利波是主波,在射流的碎裂過程中起主導作用。

        (1)平面液膜瑞利-泰勒波的曲率半徑和曲率

        沿x方向瑞利波的曲率KR見方程(20)。沿y方向泰勒波的曲率為:

        則瑞利-泰勒波曲面的曲率為:

        對方程(32)線性化,得平面液膜線性化的瑞利-泰勒波曲率。

        (2)圓射流和環(huán)狀液膜瑞利-泰勒波的曲率半徑和曲率

        圍繞z軸不再是個圓面,而是一個波動的曲面,曲面不但沿θ方向有波動,而且其波動的半徑還要隨著z方向的位移而變化。θ方向、r方向和z方向的附加壓強均有正負號。因此,圓射流瑞利波線性穩(wěn)定性分析的曲率為方程(27)變成

        環(huán)狀液膜瑞利波線性穩(wěn)定性分析的曲率為方程(28)變成

        圓射流和環(huán)狀液膜瑞利波非線性穩(wěn)定性分析的曲率為方程(29)變成

        對于圓射流的正反對稱波形和環(huán)狀液膜,非線性穩(wěn)定性分析的曲率方程(30)變成

        根據(jù)弧長與弧度的關系,ξ,s=ξ,(rθ)。對于旋轉角度,r可以看作常數(shù),有因此,沿θ方向泰勒波的曲率為

        沿r方向泰勒波的曲率,對于線性穩(wěn)定性分析

        對于非線性穩(wěn)定性分析

        沿z方向泰勒波的曲率,對于線性穩(wěn)定性分析

        對于非線性穩(wěn)定性分析

        泰勒波線性穩(wěn)定性分析的曲率為方程(38)、(39)、(41)之和,即

        對方程(43)線性化,得

        泰勒波非線性穩(wěn)定性分析的曲率為方程(38)、(40)、(42)之和,即

        圓射流瑞利-泰勒波線性穩(wěn)定性分析的曲率為方程(34)與(44)之和,即

        環(huán)狀液膜瑞利-泰勒波線性穩(wěn)定性分析的曲率為方程(35)與(44)之和,即

        瑞利-泰勒波非線性穩(wěn)定性分析的曲率為方程(36)與(45)之和。即

        目前,由于非線性穩(wěn)定性理論只研究到瑞利波,圓射流和環(huán)狀液膜采用的曲線曲率方程均為(30),而由于用于研究瑞利-泰勒波的曲面曲率方程(48)太過復雜,目前尚未有人嘗試使用。

        3 平面液膜的動力學邊界條件

        射流的動力學邊界條件為液相與氣相的應力張量差等于由表面張力引起的附加壓強。

        3.1 平面液膜的線性動力學邊界條件

        線性穩(wěn)定性理論能夠研究粘性流體。粘性液體射流必須由連續(xù)性方程和納維-斯托克斯方程組得到的擾動壓力差確定動力學邊界條件,為線性的動力學邊界條件。平面液膜有量綱形式的線性三維擾動瑞利波動力學邊界條件為

        方程(49)~(51)兩側同除以液相基流壓力、即液相噴射壓力Pl,并將氣液壓力比液流韋伯數(shù)液流歐拉數(shù)液流雷諾數(shù)代入,得瑞利波量綱一化的動力學邊界條件。

        對于一維擾動,則取方程(51)和方程(54)即可[7-11]。

        對于瑞利-泰勒波,有量綱形式的線性三維擾動瑞利波動力學邊界條件為

        將方程(55)~(57)量綱一化,得瑞利-泰勒波量綱一化的動力學邊界條件。

        對于一維擾動,則取方程(57)和(60)即可。

        對于無粘性流體,只需將方程(49)~(60)中含有μl和Rel的粘性項直接刪去即可。

        3.2 平面液膜的非線性動力學邊界條件

        目前,對于平面液膜的非線性穩(wěn)定性理論,研究的是無粘性、不可壓縮的理想流體。理想流體允許采用伯努利方程得到的壓力差來確定動力學邊界條件。對于理想流體射流的研究采用的是一維擾動的非線性穩(wěn)定性理論。與運動學邊界條件同樣,非線性穩(wěn)定性理論的擾動振幅有第一級波、第二級波和第三級波之分。目前的非線性穩(wěn)定性理論僅研究正對稱和反對稱模式的瑞利波。

        平面液膜的非線性動力學邊界條件為

        式中:Δp=pl-pgj(MPa)為液氣相壓力差。理想流體的動力學邊界條件又可稱為拉普拉斯方程。

        以速度勢函數(shù)?表示的伯努利方程氣液相一般形式為

        式中:t為時間;速度勢函數(shù)的定義為u=φ,x,w=φ,z。

        對于多級表面波,有

        式中:m= 1, 2, 3為表面波的級數(shù)為拉格朗日積分常數(shù)。對于液相,由于基流速度Vl=0,有對于氣相,由于是靜止氣體環(huán)境,則f(t)=0gj。gz為體力勢,也就是重力勢。由于細小的噴霧液滴受重力影響很小,弗勞德數(shù)(Froude number)Fr很大,因此重力勢可以忽略不計,即gz=0。方程(63)可以寫為

        將方程(64)和(65)代入方程(61),可得有量綱形式的平面液膜瑞利波動力學邊界條件,適用于對理想流體的非線性分析[19]。

        將方程(68)分別代入方程(66)和(67),可得平面液膜經(jīng)麥克勞林級數(shù)展開的有量綱和量綱一化的多級波動力學邊界條件。

        平面液膜線性化邊界和經(jīng)麥克勞林級數(shù)展開的非線性動力學邊界條件的邊界選取位于:有量綱形式的液相 -a≤z≤a;氣相a≤z≤∞(j= 1)和-∞≤z≤-a(j= 2)。量綱一形式的液相(j= 2)。其中:a為平面液膜在噴嘴出口處的半厚度。

        4 圓射流的動力學邊界條件

        4.1 圓射流的線性動力學邊界條件

        圓射流有量綱形式的線性三維擾動瑞利波動力學邊界條件為

        方程(71)~(73)兩側同除以Pl,得瑞利波量綱一化的動力學邊界條件。

        對于一維擾動,則取方程(71)和(74)即可[12-16]。

        對于瑞利-泰勒波,有量綱和量綱一化的線性三維擾動瑞利波動力學邊界條件為

        對于一維擾動,則取方程(77)和(80)即可。

        對于無粘性流體,只需將方程(71)~(82)中含有l(wèi)μ和Rel的粘性項直接刪去即可。

        4.2 圓射流的非線性動力學邊界條件

        圓射流一維擾動瑞利波正反對稱波形(階數(shù)n=0, 1)的非線性動力學邊界條件為拉普拉斯方程。

        用勢函數(shù)?表示的伯努利方程氣液相合式為

        對于多級表面波,有

        式中:m= 1, 2, 3為表面波的級數(shù)為拉格朗日積分常數(shù)。對于液相,由于基流速度Url=0,有對于氣相,由于是靜止氣體環(huán)境,則fgj(t)=0。gz可以忽略不計。方程(85)可以寫為

        將方程(86)和(87)代入方程(83),且非線性位移r=a+ξ,可得圓射流有量綱形式的非線性動力學邊界條件。

        圓射流線性化邊界和經(jīng)麥克勞林級數(shù)展開的非線性邊界選取位于:有量綱形式的液相0≤r≤a;氣相a≤r≤∞。量綱一化的液相氣相其中:a為圓射流在噴嘴出口處的半徑。

        5 環(huán)狀液膜的動力學邊界條件

        5.1 環(huán)狀液膜的線性動力學邊界條件

        環(huán)狀液膜有量綱形式的線性三維擾動瑞利波動力學邊界條件為

        方程(92)~(94)兩側同除以Pl,得瑞利波量綱一化的動力學邊界條件。

        對于一維擾動,則取方程(92)和(95)即可[17-18]。

        對于瑞利-泰勒波,有量綱和量綱一化的線性三維擾動動力學邊界條件為

        對于一維擾動,則取方程(98)和(101)即可。

        對于無粘性流體,只需將方程(92)~(103)中含有μl和Rel的粘性項直接刪去即可。

        5.2 環(huán)狀液膜的非線性動力學邊界條件

        環(huán)狀液膜一維擾動瑞利波的非線性動力學邊界條件為拉普拉斯方程,與方程(83)同式。

        用勢函數(shù)?表示的伯努利方程氣液相合式與方程(84)同式;多級表面波的與方程(85)同式。液相伯努利方程與方程(86)同式,氣相方程為

        將氣液相伯努利方程(104)和(86)代入方程(83),且非線性位移r=rj+ξj,可得環(huán)狀液膜有量綱形式的非線性動力學邊界條件。

        環(huán)狀液膜量綱一化的非線性動力學邊界條件為

        環(huán)狀液膜線性化邊界和經(jīng)麥克勞林級數(shù)展開的非線性邊界選取位于:有量綱形式的液相ri≤r≤ro;氣相 0≤r≤ri(j= i)和ro≤r≤ ∞(j= o)。量綱一形式的液相和其中:ri為環(huán)狀液膜在噴嘴出口處的內環(huán)半徑,ro為在噴嘴出口處的外環(huán)半徑,

        6 結 語

        對于射流的線性和非線性動力學邊界條件,學者們給出的各種表達式并不全面。本文從動力學邊界條件的普適性原始表達式入手,對平面液膜、圓射流和環(huán)狀液膜氣液交界面三維擾動瑞利表面波、以及瑞利-泰勒表面波的線性動力學邊界條件和一維擾動多級瑞利表面波的非線性動力學邊界條件進行了全面、系統(tǒng)的詳細推導,并給出射流的各種動力學邊界條件的適用條件。為進行射流碎裂過程研究的學者提供了各種明確的動力學邊界條件表達式。

        猜你喜歡
        瑞利液膜泰勒
        亞瑞利散斑場的二階累積量鬼成像
        考慮軸彎曲的水潤滑軸承液膜建模方法
        高空高速氣流下平板液膜流動與破裂規(guī)律
        液膜破裂對PCCS降膜的影響*
        一起綿羊泰勒焦蟲病的診斷治療經(jīng)過
        馬瑞利推出多項汽車零部件技術
        汽車零部件(2015年4期)2015-12-22 05:32:56
        瑞利波頻散成像方法的實現(xiàn)及成像效果對比研究
        豎直窄矩形通道內彈狀流中液膜特性研究
        泰勒公式的簡單應用
        河南科技(2014年14期)2014-02-27 14:12:08
        泰勒公式與泰勒級數(shù)的異同和典型應用
        欧美亚洲韩国国产综合五月天| 亚洲日韩中文字幕在线播放| 国产精品18久久久| 车上震动a级作爱视频| 无码中文字幕专区一二三| 欧美性xxxx极品高清| 亚洲av鲁丝一区二区三区黄| 麻豆精品久久久久久久99蜜桃 | 仙女白丝jk小脚夹得我好爽| 人妻中文字幕av有码在线| 国产区一区二区三区性色| 永久免费a∨片在线观看| 国产精品jizz视频| 加勒比黑人在线| 亚洲综合精品一区二区| av黄页网国产精品大全| 中文字幕人妻中文| 伊人网视频在线观看| 久久青青草原国产精品最新片| 久久午夜一区二区三区| 99久久无色码中文字幕人妻蜜柚| 色翁荡息又大又硬又粗又视频图片| 天天狠狠综合精品视频一二三区| 亚洲区福利视频免费看| 亚洲成人免费av影院| 国产午夜精品一区二区三区嫩草 | 成年人男女啪啪网站视频| 精品少妇一区二区三区免费| 激情内射日本一区二区三区| 3344永久在线观看视频| 国产亚洲精选美女久久久久 | 乱人伦中文无码视频在线观看| 国产高清精品自在线看| 日韩精品首页在线观看| 中美日韩在线一区黄色大片| 欧美性受xxxx狂喷水| 99久久综合精品五月天| 国产精品一区又黄又粗又猛又爽| 黑人大群体交免费视频| 精品久久人人妻人人做精品| 久久99老妇伦国产熟女高清|