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

龍卷風渦旋場氣流的動力學分析

2021-09-16 09:36:06
氣象與環境科學 2021年5期
關鍵詞:分析

武 圓

(北京人文大學,北京 101300)

引 言

龍卷風一般持續時間短,襲擊時間突然,破壞力強,常給人民的生命財產帶來巨大損失[1-5]。雖然有學者在17世紀就對龍卷風開始了研究[6],但直到目前關于龍卷風的問題依然是世界難題[7-8]。目前國內也有學者對龍卷風的渦旋結構與維持特性進行了理論研究與討論[9-12],但更多的學者則傾向使用數值模擬來研究探索龍卷風渦旋場[13-18]。由于人們對龍卷風渦流場缺乏足夠的認識與了解,以至于在對處于龍卷風渦旋場中的物體進行受力分析時缺乏必要的可靠理論依據。比如現行的二維Rankin平面渦旋模型就沒有考慮龍卷風渦旋對物體的軸向作用力,而一些三維龍卷風渦旋模型如Wen建立的半經驗三維理論風場模型對氣壓降考慮得不夠完善,Burgers三維龍卷風渦旋模型對龍卷風軸向運動氣流的描述與實際情況不一致,而其他的一些三維龍卷風渦旋模型也都與現實風場存在一定的差距[19]。對此,本文將根據新提出的龍卷等渦旋的動力學模型[20],分別對龍卷風渦旋場中旋轉與軸向運動氣流的運動與受力,以及物體在龍卷風渦旋場中徑向與軸向的受力,進行理論分析,以期使人們對龍卷風的渦旋現象有更深刻的認識與理解,并為學者在以后研究龍卷風等渦旋方面的問題[21-25]時提供必要可靠的動力學模型和理論計算方法。比如依據本文對龍卷風渦旋動力學的分析與理解,可進一步研究解決一般中氣旋形成龍卷風的理論判據,或一般龍卷風的生成條件,以在對龍卷風預警時可以提前更多的時間并提高預報精度,甚至在未來有可能將龍卷風災害消滅在生成階段,可以最大限度地將龍卷風災害影響降到最小。

1 龍卷風渦旋動力學模型及其渦旋場負壓差計算

1.1 龍卷風渦旋的動力學模型及其壓差變化梯度場的計算

龍卷風渦旋運動渦核區的氣流屬于以一定角速度ω進行渦旋運動的層流,其渦核區的渦旋運動氣流的壓差分布可依據龍卷風快速旋轉的氣流所具有的離心力Fa與旋轉氣壓梯度力Fp平衡的物理條件得出[20]。如圖1所示,其中,r為旋轉氣流位置到旋轉中心的距離,Δr為旋轉氣流層的厚度;Δz為旋轉氣流層動力學模型的高度,并取Δz充分小;ω為龍卷風風場旋轉氣流的角速度;ρ為旋轉氣流的密度。根據龍卷風渦旋風場的旋轉氣流層動力學模型,可得在龍卷風渦核區中的壓差p的變化梯度隨渦旋半徑的變化規律。

圖1 龍卷風渦旋風場r處旋轉氣流層的動力學模型

若取Δr→0,則渦旋氣流層在位置區間[r,r+Δr]處的質量Δm=2πr·Δr·Δz·ρ,旋轉氣流速度v=rω,且旋轉氣流層在位置區間[r,r+Δr]處承受的負壓強差為Δp(r),則承受負壓強差為Δp(r)的面Δs(r)=2πr·Δz,則根據Fa=Fp:

(1)

(2)

令Δr→0,則由式(2)可得:

(3)

式(3)表示的是在龍卷風的渦旋場內壓差變化梯度隨其半徑的變化關系。

1.2 龍卷風渦旋場中渦旋氣流承受的負壓差的計算

由式(3)可知,龍卷風壓差梯度場中r位置渦旋氣流承受的相對大氣常壓的負壓差p(r):

(4)

2 一般幾何形物體在龍卷風梯度場徑向的受力分析

根據龍卷風動力學模型及該模型的壓差變化的梯度場方程式(3),假設在該龍卷風的動力學模型中存在一物體Ω=Ω(x,y,z),其中Ω(x,y,z)為該物體的體函數,且該物體完全存在于龍卷風渦旋梯度場內,則在該物體中任取一點Ω(xi,yi,zi),以及該點的無限接近點Ω(xi+Δx,yi+Δy,zi+Δz),則以該兩點為對角點組成了一個邊長分別為Δx、Δy和Δz的矩形(如圖2)。其中rmin與rmax分別為龍卷風旋轉氣流的內外邊界。

圖2 物體在龍卷風梯度場的受力分析

顯然,該體元的體積ΔΩ=Δx·Δy·Δz。由龍卷風的動力學模型可知:該體元ΔΩ沿水平方向(坐標x軸和y軸方向)承受梯度壓強差,令該體元承受沿x軸方向上的梯度壓強差為Δp(xi),y軸方向上的梯度壓強差為Δp(yi),z軸方向上的梯度壓強差為Δp(zi),且龍卷風渦旋場鉛垂面上沒有旋轉氣流存在,所以z軸方向上的梯度壓強差Δp(zi)=0。

結合式(3),若令Δx→0,可得該矩形承受x軸方向上的梯度壓強差為Δp(xi):

(5)

對于該矩形ΔΩ,承受該梯度壓強差Δp(xi)的面定義為:ΔSyz=Δy·Δz,則該矩形承受的沿x軸方向上的梯度壓差力為ΔF(xi),該力指向渦旋中心,且令Δx→0、Δy→0、Δz→0,則:

ΔF(xi)=Δp(xi)·ΔSyz=ρω2xi·ΔΩ

(6)

當式(6)中ΔF(xi)的xi取值遍歷其所占體積Ω(x,y,z)的取值范圍,可得在該龍卷風的動力學模型中一個其體積為Ω物體的所承受的沿x軸方向上的壓差力F(x):

(7)

對比式(7)可得該物體承受的沿y軸方向上的壓差力F(y),同理該力指向渦旋中心,且

F(y)=?Ωρω2y·dΩ

(8)

于是可得該體積為Ω的物體在龍卷風中受到的徑向合力F(r):

(9)

以上分析與計算為一般物體在龍卷風的壓差梯度場的徑向受力分析過程與計算方法。式(7)(8)與式(9)表明,物體在龍卷風渦旋場中受到的徑向力等于該物體在龍卷風渦旋壓差梯度場的體積分,即氣壓梯度力對處于龍卷風渦旋場中物體的作用及表現形式。

3 龍卷風渦旋場軸向運動氣流的動力學分析

3.1 龍卷風渦旋場對接觸地面(或大氣)的作用力分析與計算

當龍卷風渦旋場兩端的負壓區暴露時,或者龍卷風渦旋場的底部接觸地面或大氣時,龍卷風渦旋場與地面(或大氣)存在接觸面。此時分析龍卷風渦旋場對接觸面的作用力。由圖2可知,龍卷風渦旋底部與地面(或大氣)的接觸面為旋轉氣流的內外邊界rmin與rmax所圍成的環面,且令該環面為D(0),0為龍卷風渦旋距離地面高度為0的接觸環面。

在龍卷風渦旋場與地面接觸的面上任取一點D(xi,yi,0)及該點的無限接近點D′(xi+Δx,yi+Δy,Δz),并組成渦旋場氣流體元,且在該點D(xi,yi,0)氣流體元與地面接觸,并結合式(4)可得,在龍卷風旋轉氣流渦旋場中,在D(xi,yi,0)位置點上渦旋場氣流體元存在一個相對于大氣壓的壓差p(xi,yi,0),則

(10)

則已知點D(xi,yi,0)與該點的無限接近點D′(xi+Δx,yi+Δy,zi+Δz)所構成的氣流體元ΔΩ=Δx·Δy·Δz的底部與地面接觸,根據壓差各向同性的原理,可以得出地面在渦旋場D(xi,yi,0)處也存在相對大氣壓的負壓差p(xi,yi,0),且該處與體元接觸的地面在壓差p(xi,yi,0)的作用下沿渦旋軸向的受力為ΔFi(0),且同時令Δx→0、Δy→0,則:

ΔFi(0)=p(xi,yi,0)·ΔxΔy

(11)

結合式(10),當式(11)中xi、yi的取值遍歷龍卷風渦旋底部與地面的接觸面D(0)的取值范圍,可得龍卷風渦旋場對接觸地面D(0)的軸向作用力F(0):

(12)

由于龍卷風渦旋底部與地面(或大氣)的接觸面為旋轉氣流的內外邊界rmin與rmax所圍成的環面,則式(12)中

并直接計算式(12),得F(0):

(13)

分析式(13)的結果可知,當渦旋氣流與地面接觸形成龍卷風時,其旋轉氣流層將與地面接觸,而此時龍卷風渦旋底部氣流層存在的相對于大氣壓的壓差p(x,y,0)將會與地面產生作用,并對地面產生一個沿渦旋場軸向的壓強吸力,該吸力將把龍卷風渦旋氣流經過地方的氣流、泥土、甚至樹木脫離地面,式(13)的計算結果即該吸力的值。

3.2 龍卷風風場內軸向運動氣流的受力分析

根據上文3.1章節的分析可知,當環境中的氣流體元通過渦旋氣流與地面(或大氣)的接觸面在壓差p(x,y,0)作用下而進入龍卷風渦旋場并沿渦旋場軸向運動,并同時在渦旋場中與作旋轉運動的渦旋場氣流形成螺旋向上復合運動的渦旋氣流,且當該軸向運動氣流體元形成螺旋向上復合運動的渦旋氣流后,則根據式(3)(4)(5)的分析可知,該進入龍卷風渦旋場中氣流體元也會因其具有旋轉運動而在該渦旋場中繼續產生壓差,即當氣流體元在龍卷風渦旋場中沿軸向運動一定距離Δz后,該軸向運動的氣流體元也會因其旋轉運動而在其螺旋上升的同時在該渦旋場D(x,y,0)處產生壓差p(x,y,0),且使后續在渦旋場底部的氣流體元也處在壓差p(x,y,0)的作用下而被持續地吸入龍卷風渦旋場中。

假設當環境中的氣流體元通過渦旋氣流與地面(或大氣)的接觸面進入龍卷風渦旋場中運動一定距離Δz后不再作軸向運動而僅作渦旋旋轉運動時,則后續在軸向壓差p(x,y,0)作用下而作上升運動的氣流體元以一定的速度動能(見式(17)的分析,負壓差對氣流體元的做功完全轉化為體元的動能,兩者完全相等)接觸于該體元的底部面元時,后續上升的氣流體元所具有的動能轉化為對該體元的底部面元的正壓強能,進而把該氣流體元在上升一定距離Δz后底部所受到旋轉氣流層產生的負壓作用平衡掉,使該體元的底部不再承受負壓差。即該上升一定距離Δz后的氣流體元的底部相當于受到的壓強為大氣常壓,而該氣流體元的上部面元仍會繼續受到旋轉氣流層產生的軸向負壓差作用,則使該氣流體元在龍卷風渦旋場中上升一定距離Δz后仍會在軸向壓差作用力下繼續沿其軸向運動。

若對在龍卷風渦旋場中軸向運動的氣流連續分析可知,氣流體元在渦旋場軸向運動過程中,始終承受渦旋旋轉氣流產生的軸向壓差作用,這與其在渦旋場軸向運動的高度無關,僅與從渦旋底部是否有連續進入渦旋場的氣流有關。這也是一般常見的龍卷風渦旋中的氣流總是螺旋連續上升的原因。也即當龍卷風渦旋底部有連續進入渦旋場的氣流并上升一定高度z時,則該連續進入渦旋場的氣流在渦旋軸向運動上仍會受到軸向壓差p(x,y,z)的作用,其中p(x,y,z)為在龍卷風渦旋場中(x,y,z)高度位置中存在的相對于大氣常壓的壓差值。

3.3 龍卷風渦旋場軸向運動氣流的速度分析與計算

根據功能守恒定律,針對體元可得ΔW=ΔEk+ΔEp,

(14)

化簡:

(15)

式(15)的計算結果為高度z上點D(x,y,z)處軸向壓差p(x,y,z)對氣流體元ΔΩ的做功,所以取其運動高度Δz→0,則式中ρ(z)gΔz項為一階無窮小項,予以忽略。則式(15)化為

(16)

帶入式(10),可得體元在壓差p(x,y,z)作用下的氣流軸向運動瞬時速度:

(17)

在龍卷風渦旋場中,雖然在高度方向上不存在因旋轉氣流存在而產生的梯度壓差,但是進入渦旋場中沿軸向運動的氣流,都會在軸向壓差p(x,y,z)的作用下,不斷以一瞬時速度v(x,y,z)沿龍卷風軸向運動。比如對于水龍卷而言,其渦旋場中存在的最大的壓差僅為一個大氣壓,一般一個大氣壓的壓差最多能使水抬升10 m左右的高度,但在水龍卷中,則能使水以一定速度被抬升上千米高度。

式(15)中的一階無窮小項ρ(z)gΔz也體現了龍卷風對進入其渦旋場的氣流(以及水流或其他物質)一直做功的方式,雖然在龍卷風渦旋場每點高度上其軸向壓差p(x,y,z)對軸向運動氣流克服重力做功幾乎為零,但這卻是龍卷風本身能量的一種耗散方式,也是龍卷風渦旋存在周期短的重要原因之一。

根據式(17)可知,龍卷風在接觸地面或水面的瞬間,會吸入塵土或水,使其軸向運動氣流的密度大于龍卷的水平旋轉氣流的密度,其龍卷渦旋氣流的旋轉速度會大于其軸向瞬時速度。但當龍卷風在接觸地面或水面一段時間后,其軸向運動氣流的密度會與龍卷的水平旋轉氣流的密度相等,龍卷渦旋氣流的旋轉速度也會與其軸向瞬時速度相等。

4 一般幾何形物體在龍卷風渦旋場軸向的受力分析與計算

假設在一存在連續軸向運動氣流的龍卷風渦旋場內存在一物體,且根據上文3.2與3.3章節的分析可知,該物體的底部存在連續軸向運動的氣流,則該物體的底部在后續軸向上升氣流動壓的作用下,該物體底部受到旋轉氣流層產生的負壓作用被平衡,受到的壓強為大氣常壓。且由于渦旋氣流壓差的各向同性的特點,則該物體的上表面會受到渦旋場中旋轉氣流產生的相對于大氣常壓的負壓差的作用,使該物體在渦旋場中承受軸向壓差,并沿軸向受力。圖3為物體在龍卷風渦旋場的軸向受力示意圖。

圖3 物體在龍卷風渦旋場的軸向受力示意圖

ΔF(zi)=p(xi,yi,zi)ΔSicosγi

(18)

根據高等數學中關于面積分的分析,當式(18)中ΔSi的取值遍歷該物體的表面方程S=S(x,y,z)的上表面時,取面元ΔSi→0時,可得該物體Ω=Ω(x,y,z)在龍卷風渦旋場所受到的軸向合力F(z):

(19)

式(19)即為該物體Ω=Ω(x,y,z)在龍卷風渦旋場中軸向受力的合力F(z)及其計算方法,其中Sxy為該物體的上表面方程S=S(x,y,z)在龍卷風渦旋場中x—o—y平面上的投影。

5 龍卷風渦旋場內邊界(風眼)成因分析

5.1 龍卷風上升氣流在其渦旋平面上的物質通量分析與計算

根據式(17)可知,在龍卷風渦旋場內,沿軸向運動氣流高度z上的瞬時速度v(x,y,z)與其到渦旋中心的距離r不等而不同。假設,若該龍卷風渦旋穩衡存在,則該龍卷風在高度方向上旋轉角速度ω相等,且令龍卷風渦旋場中的沿軸向運動的氣流密度與旋轉運動的氣流密度ρ相等,則根據式(17),可得沿軸向運動氣流高度z上的瞬時速度v(x,y,z)=rω,即渦旋場的上升速度與其旋轉速度相等,即渦旋氣流層在位置區間[r,r+Δr]處Δt時間內的沿軸向運動的物質通量QΔt(r):

QΔt(r)=2πrΔr·rω·Δtρ

(20)

取Δr→0,可得在龍卷風渦旋z高度平面D(z)上在時間Δt內沿軸向運動的物質通量QΔt:

(21)

則單位時間內通過龍卷風渦旋z高度平面上的物質通量Q=QΔt/Δt,即

(22)

該龍卷風渦旋z高度平面上的物質通量Q可以估算龍卷風渦旋場對沿軸向運動物質的運輸效率。

5.2 龍卷風渦旋場內邊界rmin成因分析

下面根據物質通量的守恒原理,對龍卷風內邊界rmin給出其存在機制原理,并給出其存在的條件。在龍卷風軸向上任取兩個與其軸向垂直的兩個面D(z1)與D(z2),其中z1與z2分別表示龍卷風在不同高度位置的環面(龍卷風渦旋氣流在不同高度上內外邊界rmin、rmax組成的環面,見圖2),且z2>z1,并根據一般龍卷風在實際存在中一般為“漏斗”狀,即D(z2)高度位置的渦旋外邊界半徑rz2max一般大于或等于D(z1)高度位置的渦旋外部邊界半徑rz1max,則令:rz2max=λrz1max,其中λ≥1。且根據式(17),假設在該龍卷風渦旋場內,流經面D(z1)與D(z2)氣流的密度分別為ρ(z1)、ρ(z2),即在同一渦旋面上,旋轉與軸向運動氣流的密度均相等。且令龍卷風渦旋場在D(z1)面的內邊界rz1min=0,龍卷風在高度范圍內其旋轉角速度均不變,為ω。則根據式(22),可得龍卷風渦旋場在D(z1)面沿軸向運動氣流的物質通量Q(z1):

(23)

同理依據式(22)可得龍卷風渦旋場在D(z2)面沿軸向運動氣流的物質通量Q(z2):

(24)

根據物質通量的守恒原理,可得龍卷風渦旋場在D(z1)面沿軸向運動氣流的物質通量Q(z1)與其渦旋場在D(z2)面沿軸向運動氣流的物質通量Q(z2)必然相等,即Q(z1)=Q(z2)。于是聯立式(23)與式(24),并移項化簡可得

(25)

式(25)即是龍卷風渦旋場內部邊界存在的決定式。根據式(25)可知,若流經面D(z1)與D(z2)的氣流密度滿足關系λ3ρ(z2)=ρ(z1),則龍卷風渦旋場在面D(z2)位置上不會存在一個邊界為rz2min的區域;同理,若λ3ρ(z2)-ρ(z1)>0,則龍卷風渦旋場在面D(z2)位置上一定會產生一個不為零的邊界區域rz2min。由以上分析可知,龍卷風內部邊界rmin是由物質通量的守恒原理的限制產生的。

6 結論與應用展望

6.1 結 論

(1)根據龍卷風場的動力學模型,以及由該動力學模型得出的龍卷風渦旋壓差變化梯度場gradp(r)=ρω2r·er,進一步研究分析得出一般物體在龍卷風梯度場中所承受的徑向力即等于該物體在該梯度場中的體積分;并依據龍卷風渦旋場的壓差變化梯度場,計算出渦旋氣流承受的相對于周圍大氣壓的負壓差的分布規律(式4);依據式(4),并結合理想流體中壓差的各向同性的特性,對龍卷風渦旋場中的氣流在軸向上的受力與運動進行了理論分析,并得出其氣流的軸向運動是由其渦旋場中水平旋轉氣流運動影響產生的,渦旋氣流的水平運動與軸向運動是不可獨立對待、分析的。這也是文中關于龍卷風渦旋動力學模型中,僅設定存在龍卷風渦旋氣流的水平旋轉運動,卻依然可分析其渦旋氣流的軸向運動的原因。通過對物體在龍卷風渦旋場軸向的受力分析與研究,得到物體的軸向受力表達式式(19)。

(2)論文依據對龍卷風渦旋場中軸向運動氣流的分析,對龍卷風軸向運動氣流的物質通量進行了分析與計算,并給出龍卷風內部邊界rmin(或風眼)的存在需要滿足式(25)要求的條件:λ3ρ(z2)-ρ(z1)>0,即該風眼產生的原因是龍卷風軸向運動氣流遵守物質通量守恒的結果。

6.2 應用展望

文中依據新的龍卷風渦旋的動力學模型,分析得出的物體在龍卷風渦旋場中的徑向與軸向受力計算表達式,以及相關的分析過程與計算方法,都具有重要的應用意義,可以預見的有以下三個方面的應用:

(1)建筑風工程理論研究方面。比如對國家核電站等重要或其他安全等級要求較高的建筑設施,需要考慮對龍卷風渦旋的受力校核,那就可以根據該文章所得出的式(7)(8)(9)與式(19)分別對處于龍卷風渦旋場的設施與建筑物沿渦旋徑向與軸向受力進行理論分析與計算。

(2)龍卷風的數值模擬與有限元分析方面。如果現實中的一些物體(比如一些橋梁和鐵塔)是異形的,難以在龍卷風渦旋場中對物體的外形Ω=Ω(x,y,z)及其在渦旋軸向投影Sxy進行參數化方程建模,也即無法直接應用文中式(7)(8)(9)與式(19)對處于渦旋場中的物體進行受力計算。這時就需要對物體進行三維仿真建模,并同時通過流體有限元分析軟件,對龍卷風渦旋場進行參數化定義與描述,并通過本文對物體在龍卷風渦旋場中的受力分析方法,將文中關于物體的線元、面元與體元按精度要求具體賦值,也即將物體在龍卷風渦旋場中進行網格參數化,進而對物體在龍卷風渦旋場中的徑向與軸向受力進行有限元分析與計算。

(3)航空飛行器設計方面。目前世界最先進的飛行器都注重對于翼面渦升力的利用以增加飛行器飛行時的升力,比如我國的殲20威龍戰機。具體為利用飛行器機翼前方的邊條等結構產生脫體渦旋并流經于機翼上表面,顯然該脫體渦旋的軸向與翼面平行。根據式(7)(8)(9)的分析可知,此時渦旋在流經機翼上表面而對機翼產生的徑向作用力即為渦旋對機翼產生的渦升力,也即本文相關的分析方法與計算結果在航空航天飛行器設計上也有應用與參考的價值。

猜你喜歡
分析
禽大腸桿菌病的分析、診斷和防治
隱蔽失效適航要求符合性驗證分析
電力系統不平衡分析
電子制作(2018年18期)2018-11-14 01:48:24
電力系統及其自動化發展趨勢分析
經濟危機下的均衡與非均衡分析
對計劃生育必要性以及其貫徹實施的分析
現代農業(2016年5期)2016-02-28 18:42:46
GB/T 7714-2015 與GB/T 7714-2005對比分析
出版與印刷(2016年3期)2016-02-02 01:20:11
中西醫結合治療抑郁癥100例分析
偽造有價證券罪立法比較分析
在線教育與MOOC的比較分析
主站蜘蛛池模板: 色婷婷色丁香| 亚洲第一中文字幕| 成人精品在线观看| 欧美精品高清| 亚洲国产清纯| 成人福利在线观看| 99ri国产在线| 免费在线播放毛片| 欧美日韩在线第一页| 2021精品国产自在现线看| 国产剧情伊人| 日本欧美视频在线观看| 国产色婷婷视频在线观看| 久久综合亚洲色一区二区三区| 久久窝窝国产精品午夜看片| 2022国产无码在线| 88国产经典欧美一区二区三区| 九九热视频精品在线| www.日韩三级| 久久婷婷色综合老司机| 中文字幕在线永久在线视频2020| 婷婷激情亚洲| 久久熟女AV| 亚洲 成人国产| 一本综合久久| 国产成人8x视频一区二区| 高清久久精品亚洲日韩Av| 亚洲天堂视频在线观看| AV片亚洲国产男人的天堂| 性色一区| 欧美色综合久久| 毛片久久网站小视频| 99久久亚洲精品影院| 日韩无码黄色网站| h视频在线播放| 国产精品免费p区| 亚洲精品爱草草视频在线| 亚洲中文字幕国产av| 精品视频一区二区观看| 欧美日韩国产成人高清视频| 18禁黄无遮挡免费动漫网站| 国产日韩精品欧美一区喷| 久久99国产乱子伦精品免| 亚洲第一页在线观看| 美女潮喷出白浆在线观看视频| 在线播放国产99re| 亚洲一级毛片在线播放| 亚洲大尺码专区影院| 国产精品第一区| 亚洲香蕉伊综合在人在线| 99久久人妻精品免费二区| 中文字幕永久视频| yy6080理论大片一级久久| 少妇精品久久久一区二区三区| 国产免费福利网站| 成年人国产视频| 色综合久久88| 欧美国产日韩另类| 国产成人精品高清不卡在线 | 日本欧美视频在线观看| 亚洲热线99精品视频| 夜夜操天天摸| 亚洲天堂网视频| 内射人妻无码色AV天堂| 国产精品手机视频一区二区| 亚洲精品久综合蜜| 国产二级毛片| 国产在线观看成人91| 日韩少妇激情一区二区| 性69交片免费看| 无码网站免费观看| 九九热精品视频在线| 午夜福利无码一区二区| 综合色88| 国产在线自乱拍播放| 综合色88| 国产无码精品在线播放| 亚洲欧美日韩高清综合678| 国产情侣一区二区三区| 日本在线视频免费| 久久永久免费人妻精品| 中文无码伦av中文字幕|