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

可壓縮液體中氣泡的聲空化特性*

2022-01-19 04:44:12鄭雅欣那仁滿都拉
物理學報 2022年1期
關(guān)鍵詞:模型

鄭雅欣 那仁滿都拉

(內(nèi)蒙古民族大學數(shù)理學院,通遼 028043)

利用新提出的Gilmore-NASG 模型,在考慮液體可壓縮效應(yīng)的邊界條件下,研究了可壓縮液體中氣泡的聲空化特性,并與利用原有KM-VdW 模型計算得到的結(jié)果進行了比較.結(jié)果表明,相比于KM-VdW 模型,由于Gilmore-NASG 模型采用新的狀態(tài)方程來描述氣體、液體以及由可壓縮性引起的液體密度變化及聲速變化,所以用Gilmore-NASG 模型得到的空化氣泡的壓縮比更大、崩潰深度更深、溫度和壓力峰值更高.隨著驅(qū)動聲壓幅值的增大,兩種模型給出的結(jié)果差別愈加明顯,而隨著驅(qū)動頻率的增大,兩種模型給出的結(jié)果差別逐漸減小.這表明,在充分考慮泡內(nèi)氣體、周圍液體在不同溫度和壓強下共體積的變化所導(dǎo)致的介質(zhì)可壓縮特性下,氣泡內(nèi)的溫度和壓強可能達到更高值.同時,Gilmore-NASG 模型還預(yù)測出了氣泡壁處液體的密度變化、壓力變化、溫度變化以及液體中的聲速變化.因此,Gilmore-NASG 模型在研究高壓狀態(tài)下氣泡的空化特性以及周圍液體對氣泡空化特性的影響方面具有優(yōu)點.

1 引言

對于可壓縮液體中的聲空化氣泡,通常使用的氣泡動力學方程有Keller-Miksis[1](KM)方程和Gilmore[2]方程,兩者都是利用連續(xù)介質(zhì)力學的基本運動方程,在氣泡壁處液體馬赫數(shù)的一級近似下簡化得到[3-7].但兩種方程簡化時采用的條件有所不同,KM 方程中認為氣泡周圍液體密度和聲速是常量,而Gilmore 方程中認為液體密度和聲速是隨液體壓力變化的量.KM 方程應(yīng)用于聲致發(fā)光及聲化學問題的研究表明,此方程在低驅(qū)動聲壓下較為準確[8,9].因為在高驅(qū)動聲壓下,氣泡劇烈崩潰時產(chǎn)生的泡壁速度可以超過液體中的恒定聲速,導(dǎo)致氣泡壁處液體馬赫數(shù)異常增大,這就違背了KM 方程所基于的小馬赫數(shù)假設(shè).Gilmore 方程則考慮了液體密度變化以及由此導(dǎo)致的聲速變化,這有效地降低了氣泡劇烈崩潰時氣泡壁處液體馬赫數(shù)瞬間變得很大的可能,所以Gilmore 方程更適合于研究高驅(qū)動聲壓下的聲空化問題[10-12].由于KM 方程中認為氣泡壁處液體密度和聲速是恒定不變的,所以通常與van der Waals (VdW)狀態(tài)方程結(jié)合構(gòu)成KM-VdW 模型[13,14].Gilmore 方程中認為液體密度和聲速是隨液體壓力變化的量,所以通常與Tait 狀態(tài)方程結(jié)合構(gòu)成Gilmore-Tait 模型[2,15,16].Tait 狀態(tài)方程可以描述諸如水這樣的液體介質(zhì)的可壓縮性,給出了介質(zhì)密度與壓力之間的變化關(guān)系,但它不能直接預(yù)測介質(zhì)的溫度,也不能反映介質(zhì)熱容的影響,且Tait 方程中的多方指數(shù)選取得異常大,不符合熱力學關(guān)系[17].這說明,需要一種更合適的狀態(tài)方程來描述液體介質(zhì)的可壓縮性以及由此引起的密度變化和聲速變化[18-20].最近,Denner[21]把Gilmore 方程與Noble-Abel-Stiffend-Gas (NASG)狀態(tài)方程相結(jié)合,提出了一種新模型,即Gilmore-NASG 模型,并與傳統(tǒng)的Gilmore-Tait 模型進行了比較.從兩種模型對預(yù)測氣泡壁處液體溫度的結(jié)果來看,Gilmore-Tait 模型給出的氣泡壁處液體溫度高于氣泡內(nèi)的溫度,這是錯誤的結(jié)果,說明Tait 狀態(tài)方程自身有所缺欠.另外注意到,目前研究氣泡空化問題時在氣泡壁處采用的壓力邊界條件,多數(shù)情況下都是只考慮液體切變黏滯影響,未考慮液體體積黏滯影響的邊界條件.由Navier-Stokes 方程以及黏性應(yīng)力張量可知[4],當液體不可壓縮時,其速度的散度為零,此時體積黏滯系數(shù)相關(guān)的黏滯應(yīng)力為零,這導(dǎo)致體積黏滯作用對氣泡壁的運動不會產(chǎn)生影響;只有液體可壓縮時,其速度的散度不為零,體積黏滯作用才會對氣泡壁的運動有影響[22].因此,研究可壓縮液體中氣泡空化時,為體現(xiàn)液體的可壓縮性,應(yīng)使用包含切變粘滯作用和體積粘滯作用的邊界條件.

為了充分體現(xiàn)液體可壓縮性,本文將利用新提出的Gilmore-NASG 模型,在考慮液體可壓縮效應(yīng)的邊界條件下,數(shù)值研究可壓縮液體中空化氣泡的聲空化特性,并把結(jié)果與KM-VdW 模型給出的結(jié)果進行比較分析.

2 氣泡動力學模型

2.1 Gilmore-NASG 模型

Denner[21]把Gilmore 方程與NASG 狀態(tài)方程相結(jié)合,提出了Gilmore-NASG 模型[10].經(jīng)典的Gilmore 方程為

式中,R=R(t)為氣泡瞬時半徑,點表示對時間的導(dǎo)數(shù),H為液體在泡壁處和無窮遠處焓的差,Cl為液體中的聲速.此模型采用一種新的狀態(tài)方程,即NASG 狀態(tài)方程[20]來統(tǒng)一描述氣泡內(nèi)氣體、蒸氣和氣泡外液體的狀態(tài).NASG 狀態(tài)方程可表示為

式中,P為壓力;T為溫度;CV為熱容;Γ為多方指數(shù);V為比容,與密度ρ的關(guān)系為V=1/ρ;B為壓力常數(shù),表示分子吸引效應(yīng);b為共體積(分子所占體積),表示分子排斥效應(yīng).相比于其他狀態(tài)方程,NASG 狀態(tài)方程能夠統(tǒng)一描述氣相、液相以及氣液共存相,且滿足熱力學關(guān)系[18].由于考慮了介質(zhì)熱容的影響,所以能夠準確預(yù)測氣、液介質(zhì)的溫度.此外,狀態(tài)方程(2)中的共體積b和壓力常數(shù)B是根據(jù)實驗飽和曲線,在不同溫度和壓力下取點計算確定的,所以更能夠準確反映分子排斥效應(yīng)和吸引效應(yīng)[20].

利用狀態(tài)方程(2),在不考慮傳質(zhì)傳熱的絕熱情況下,把液體的焓差和液體中的聲速可表示為

而常數(shù)kl和驅(qū)動聲壓Pex(t)為

同時,泡內(nèi)壓力Pg、泡內(nèi)溫度Tg及泡壁溫度Tl可表示為

2.2 KM-van der Waals 模型

KM-VdW 模型是KM 方程與VdW 狀態(tài)方程相結(jié)合得到的,其中KM 方程為

由VdW 狀態(tài)方程得到的泡內(nèi)壓力Pg及泡內(nèi)溫度Tg的表達式為

(13)式—(15)式中,m為氣體摩爾數(shù),V為氣泡體積,其他量與Gilmore-NASG 模型中所代表的含義相同.以上各式中,下角標“l(fā)”表示液體相關(guān)量,“g”表示氣體相關(guān)量,“∞”表示無窮遠處的量,“0”表示靜態(tài)或環(huán)境量.

2.3 氣泡壁處壓力邊界條件

用KM-VdW 模型計算氣泡空化問題時,由于認為氣泡壁處液體密度和液體中的聲速都是不變的常量,所以通常使用不考慮液體可壓縮效應(yīng)的邊界條件[23,24],這也是合理的.然而,在Gilmore-NASG模型中認為氣泡壁處液體密度和液體中的聲速都是隨壓力變化的,即認為氣泡壁處液體是可壓縮的,所以采用考慮液體可壓縮效應(yīng)的邊界條件更合適.為此,本文采用(16)式表示的考慮液體可壓縮效應(yīng)的邊界條件[25]:

式中σ,η,λ分別代表液體的表面張力系數(shù)、切變黏滯系數(shù)、體積黏滯系數(shù).當?shù)忍栍疫叺淖詈笠豁椀扔诹銜r,方程(16)就變成我們熟悉的不考慮液體可壓縮效應(yīng)的邊界條件.

3 數(shù)值模擬

選擇水為液體介質(zhì),泡內(nèi)氣體為氬氣,氣泡初始半徑R0=4.5 μm,驅(qū)動聲壓振幅A=1.4 atm,頻率f=30 kHz.兩種模型相關(guān)的物理參數(shù)如表1和表2所列.Gilmore-NASG 模型中參數(shù)Bl和bl的計算方法以及液體、氣體環(huán)境參數(shù)的選取可參考文獻[17,20,26].關(guān)于泡內(nèi)氣體分子共體積,兩種模型所使用的值不同.在NASG 狀態(tài)方程中,在不同的溫度和壓力下,氣體分子共體積有不同的取值.本文使用的氬氣共體積是根據(jù)文獻[27,28]給出的計算方法,在溫度為20000 K 時計算得到的值.VdW 狀態(tài)方程中的參數(shù)a和b,通常認為是不隨溫度和壓力變化的常量,本文選取文獻[29]給出的值,具體計算時取a=0,這是為了與NASG 方程中Bg=0對應(yīng).為了方便描述,下面用模型Ⅰ和模型Ⅱ分別代表Gilmore-NASG 模型和KMVdW 模型.

表1 Gilmore-NASG 模型中水和氬氣的相關(guān)物理參數(shù)Table 1.Physical parameters of water and argon in Gilmore-NASG model.

表2 KM-VdW 模型中水和氬氣的相關(guān)物理參數(shù)Table 2.Physical parameters of water and argon in KM-VdW model.

3.1 可壓縮液體中氬氣泡的空化特性

本節(jié)研究在一定強度和一定頻率的超聲波作用下,可壓縮水中氬氣泡的空化特性.用模型Ⅰ計算了空化氣泡的半徑變化、速度變化、泡內(nèi)壓力變化、溫度變化以及氣泡壁處液體的密度變化、壓力變化、溫度變化、聲速變化、馬赫數(shù)變化,并與用模型Ⅱ得到的結(jié)果進行了比較,結(jié)果如圖1 所示.從圖1(a)可以看出,氣泡的膨脹過程是一個相對緩慢的過程,在此過程中兩種模型給出的氣泡半徑變化并無明顯差別,由此可知在膨脹階段液體可壓縮性對氣泡空化行為的影響很小.從氣泡第一次崩潰時的最小半徑的對比可看出(圖1(a)插圖),模型Ⅰ給出的氣泡最小半徑更小,即有著更大的壓縮比.出現(xiàn)這種現(xiàn)象的主要原因是,兩種模型中使用的狀態(tài)方程不同,所選用的泡內(nèi)氣體分子共體積(bg)和泡壁處液體分子共體積(bl)不同.模型Ⅰ中采用NASG 狀態(tài)方程統(tǒng)一描述氣泡內(nèi)氣體和氣泡壁處液體,充分考慮氣泡內(nèi)氣體的可壓縮性,認為不同溫度和壓強下氣體分子共體積應(yīng)取不同值,在高溫高壓下氣體分子共體積要小于VdW 狀態(tài)方程中的氣體分子共體積(bg<b),共體積小意味著氣體被壓縮的空間更大,所以氣泡的崩潰深度更深(即最小半徑更小).同時,NASG 狀態(tài)方程中充分考慮了氣泡壁處液體的可壓縮性對氣泡崩潰行為的影響,即考慮了液體分子共體積的影響.液體分子共體積大,意味著液體分子之間的排斥作用強,液體不易被壓縮,所以氣泡壁處液體對于氣泡的壓力大,使氣泡崩潰得更深.自然,崩潰深度越深的氣泡,泡內(nèi)氣體壓力越大(如圖1(c)),所以氣泡回彈的速度也越大(如圖1(b)),這就導(dǎo)致模型Ⅰ給出的氣泡第一次崩潰后的回彈半徑比模型Ⅱ給出的回彈半徑更大(如圖1(a)).這與文獻[15,26]給出的結(jié)果相同.圖1(c)和圖1(d)給出的是氣泡內(nèi)壓力(Pg)和氣泡壁處液體壓力(Pl)的變化,可以看出,模型Ⅰ給出的Pg和Pl峰值均明顯高于模型Ⅱ給出的值,這種差距主要是由兩種模型給出的氣泡崩潰深度不同引起的,模型Ⅰ給出的氣泡崩潰深度更深,所以模型Ⅰ給出的Pg和Pl峰值更大.比較Pg和Pl可知,模型Ⅰ給出的Pg峰值高于Pl峰值約0.3 MPa,模型Ⅱ給出的Pg峰值高 于Pl峰值約0.8 MPa.圖1(e)和圖1(f)給出的是氣泡壁處液體密度變化和液體中聲速變化,可看出模型Ⅰ給出的液體密度和液體中的聲速在氣泡崩潰時明顯增大,這是由于氣泡崩潰時氣泡壁處液體壓力發(fā)生變化(圖1(d))所導(dǎo)致.這一點從液體密度公式(5)和液體中的聲速公式(4)也可直接看出.與模型Ⅰ給出結(jié)果相比,模型Ⅱ給出的液體密度和液體中的聲速都是常數(shù)(圖中黑虛線表示).由此可看出,模型Ⅰ更充分體現(xiàn)了液體的可壓縮性.兩種模型預(yù)測出的泡內(nèi)溫度變化如圖1(g)所示.可以看出,模型Ⅰ預(yù)測出的泡內(nèi)溫度高于模型Ⅱ預(yù)測出的溫度,可達21250 K,該溫度值比較接近于安宇[30]預(yù)測出的水中氬氣泡的泡內(nèi)溫度.模型Ⅰ還預(yù)測出了氣泡壁處液體的溫度變化(圖1(h)),但模型Ⅱ不能預(yù)測氣泡壁處液體的溫度變化.模型Ⅰ預(yù)測出的氣泡崩潰時的液體最高溫度約為430 K,實際上在絕熱過程中氣泡壁處液體溫度不應(yīng)該達到此值(已高于環(huán)境溫度300 K),所以可理解為氣泡壁處液體溫度的上限.圖1(i)為兩種模型給出的馬赫數(shù)對比,可以看出模型Ⅰ給出的馬赫數(shù)峰值明顯小于模型Ⅱ給出的馬赫數(shù)峰值,這是因為模型Ⅰ中認為液體中的聲速是隨壓力變化的,所以氣泡崩潰瞬間氣泡壁的速度很大時液體中的聲速也變得很大,這樣就有效地降低了氣泡壁處液體馬赫數(shù)瞬間變得很大的可能,保證了模型Ⅰ的有效性.不管是KM 方程還是Gilmore 方程都是在小馬赫數(shù)假設(shè)條件下簡化得到,所以只有在小馬赫數(shù)條件下兩種方程才能保證它們的有效性.

圖1 (a)氣泡半徑變化;(b)氣泡壁速度變化;(c)氣泡內(nèi)壓力變化;(d)氣泡壁處液體壓力變化;(e)氣泡壁處液體密度變化;(f)氣泡壁處液體中聲速變化;(g)氣泡內(nèi)溫度變化;(h)氣泡壁處液體溫度變化;(i)氣泡壁處液體馬赫數(shù)變化Fig.1.(a)Change of bubble radius;(b)change of bubble wall velocity;(c)change of pressure in the bubble;(d)change of liquid pressure at the bubble wall;(e)change of liquid density at the bubble wall;(f)change of sound velocity in the liquid at the bubble wall;(g)change of temperature in the bubble;(h)change of liquid temperature at the bubble wall;(i)change of liquid Mach number at the bubble wall;.

3.2 驅(qū)動聲壓對可壓縮液體中氬氣泡空化特性的影響

本節(jié)主要考察驅(qū)動聲壓變化對可壓縮液體中氬氣泡空化特性的影響.利用模型Ⅰ和模型Ⅱ,在驅(qū)動頻率不變的情況下分別計算了氣泡最大半徑、最小半徑、崩潰速度、回彈速度、泡內(nèi)最大壓力、泡內(nèi)最高溫度以及氣泡壁處液體最大馬赫數(shù),結(jié)果如圖2 所示.可以看出,隨著驅(qū)動聲壓幅值的增加,兩種模型給出的氣泡最大半徑均逐漸增大(圖2(a)),氣泡最小半徑均逐漸減小(圖2(b)),氣泡壓縮比均呈增大趨勢,這與文獻[31]所得結(jié)果相同.相比較,兩種模型給出的氣泡最大半徑變化基本相同,最小半徑變化出現(xiàn)了明顯的差異,模型Ⅰ給出的氣泡最小半徑更小,即氣泡崩潰得更深,并且隨著驅(qū)動聲壓幅值的增大,兩種模型給出的氣泡最小半徑的差異逐漸增大,這結(jié)果與文獻[15,26]所得結(jié)果相同.從圖2(c)可以看出,隨著驅(qū)動聲壓幅值的增大,氣泡崩潰速度和回彈速度均逐漸增大.相比可知,兩種模型給出的氣泡崩潰速度差別不大,但回彈速度出現(xiàn)了明顯差別,模型Ⅰ給出的回彈速度更大,這可以導(dǎo)致氣泡回彈半徑更大.對比兩種模型給出的最大馬赫數(shù)(圖2(d))可知,隨著驅(qū)動聲壓幅值的增大,兩種模型給出的最大馬赫數(shù)都在增大,但模型Ⅱ給出的最大馬赫數(shù)相比于模型Ⅰ增大得更多.這表明,模型Ⅰ更適合于研究高驅(qū)動聲壓下的聲空化問題.兩種模型給出的泡內(nèi)最大壓力和泡內(nèi)最高溫度,隨驅(qū)動聲壓幅值的增大而變化情況,如圖2(e)和圖2(f)所示.可以看出,隨著驅(qū)動聲壓幅值的增大,泡內(nèi)最大壓力和最高溫度都在增大,比較而言,模型Ⅰ預(yù)測出的泡內(nèi)最大壓力和最高溫度增大得非常明顯.綜上可知,隨著驅(qū)動聲壓幅值的增大,兩種模型給出的結(jié)果差別愈加明顯,模型Ⅰ對高壓狀態(tài)下氣泡空化行為的描述更加顯現(xiàn).

圖2 (a)氣泡最大半徑、(b)最小半徑、(c)崩潰速度和回彈速度、(d)液體最大馬赫數(shù)、(e)泡內(nèi)最大壓力以及(f)泡內(nèi)最高溫度,隨驅(qū)動聲壓幅值的變化Fig.2.(a)Maximum bubble radius,(b)minimum radius,(c)collapse speed and rebound speed,(d)maximum liquid Mach number,(e)maximum pressure in the bubble,and (f)maximum temperature in the bubble change with amplitude of driving sound pressure.

3.3 驅(qū)動頻率對可壓縮液體中氬氣泡空化特性的影響

本節(jié)在驅(qū)動聲壓幅值不變的條件下,主要考察驅(qū)動頻率變化對可壓縮液體中氬氣泡空化特性的影響.利用模型Ⅰ和模型Ⅱ,分別計算了氣泡最大半徑、最小半徑、崩潰速度、回彈速度、泡內(nèi)最大壓力、泡內(nèi)最高溫度以及氣泡壁處液體最大馬赫數(shù),結(jié)果如圖3 所示.從圖3(a)和圖3(b)可以看出,隨著驅(qū)動頻率的增大,兩種模型給出的氣泡最大半徑減小,最小半徑增大,氣泡壓縮比均呈減小趨勢,與文獻[15,23]所得結(jié)果一致.主要原因是隨著驅(qū)動頻率的增大,氣泡膨脹及壓縮周期相對縮短[26,32],氣泡沒有足夠的時間膨脹到更大和壓縮至更小.兩種模型給出的氣泡最大半徑無明顯差別,最小半徑有明顯差別,隨著驅(qū)動頻率的增大差別逐漸減小,并出現(xiàn)高頻激勵下趨于一致的走向.兩種模型給出的氣泡崩潰速度和回彈速度(圖3(c))的差別也隨驅(qū)動頻率的增大而減小,并出現(xiàn)逐漸趨于一致的規(guī)律.隨著驅(qū)動頻率的增大,兩種模型給出的最大馬赫數(shù)(圖3(d))、泡內(nèi)最大壓力(圖3(e))和泡內(nèi)最高溫度(圖3(f))均出現(xiàn)逐漸趨于一致的規(guī)律.總體上,隨著驅(qū)動頻率的增大,兩種模型給出的結(jié)果差別逐漸減小,并出現(xiàn)高頻激勵時趨于一致的規(guī)律.

圖3 (a)氣泡最大半徑、(b)最小半徑、(c)崩潰速度和回彈速度、(d)液體最大馬赫數(shù)、(e)泡內(nèi)最大壓力 以及(f)泡內(nèi)最高溫度隨驅(qū)動頻率的變化Fig.3.(a)Maximum bubble radius,(b)minimum radius,(c)collapse speed and rebound speed,(d)maximum liquid Mach number,(e)maximum pressure in the bubble,and (f)maximum temperature in the bubble change with driving frequency.

4 結(jié)論

本文研究了可壓縮水中氬氣泡的聲空化特性.為充分體現(xiàn)氣泡內(nèi)氣體和周圍液體的可壓縮性,采用了Gilmore 方程與NASG 狀態(tài)方程相結(jié)合的新模型,并在氣泡壁處采用了考慮液體可壓縮效應(yīng)的邊界條件.NASG 狀態(tài)方程能夠統(tǒng)一描述氣泡內(nèi)氣體、蒸氣以及周圍液體在不同溫度和壓強下的狀態(tài),能夠充分體現(xiàn)介質(zhì)分子排斥效應(yīng)和吸引效應(yīng),可以準確表達氣、液介質(zhì)密度變化和壓強變化,能夠考慮到介質(zhì)的熱容效應(yīng),可以準確預(yù)測氣、液介質(zhì)的溫度變化.在不考慮氣泡內(nèi)氣體與周圍液體之間的質(zhì)量交換、熱交換和化學反應(yīng)的情況下,相比于KM-VdW 模型給出的結(jié)果,用Gilmore-NASG模型計算得到的空化氣泡的壓縮比更大,崩潰深度更深,溫度和壓力峰值更高.隨著驅(qū)動聲壓幅值的增大,兩種模型給出的結(jié)果差別愈加明顯,而隨著驅(qū)動頻率的增大,兩種模型給出的結(jié)果差別逐漸減小.這表明,在充分考慮泡內(nèi)氣體、周圍液體在不同溫度和壓強下共體積的變化所導(dǎo)致的介質(zhì)可壓縮特性下,氣泡內(nèi)的溫度和壓強可能達到更高值.Gilmore-NASG 模型還能準確地預(yù)算出氣泡壁處液體的密度變化、壓力變化、溫度變化以及液體中的聲速變化,但KM-VdW 模型不能預(yù)算這些量的變化.這一優(yōu)勢對于研究周圍液體對氣泡空化特性的影響非常重要.總體而言,Gilmore-NASG 模型在高壓狀態(tài)下對氣泡空化特性的研究以及周圍液體對氣泡空化特性影響的研究方面具有優(yōu)勢,在高強度聚焦超聲、沖擊波碎石治療以及聲化學等具體問題的研究中將會有重要應(yīng)用.

猜你喜歡
模型
一半模型
一種去中心化的域名服務(wù)本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數(shù)模型及應(yīng)用
p150Glued在帕金森病模型中的表達及分布
函數(shù)模型及應(yīng)用
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 国产永久在线视频| 国产永久在线视频| 欲色天天综合网| 狼友av永久网站免费观看| 久久久噜噜噜| 国产亚洲欧美日韩在线观看一区二区 | 69精品在线观看| P尤物久久99国产综合精品| 午夜毛片福利| 99精品在线视频观看| 女人毛片a级大学毛片免费| 美女扒开下面流白浆在线试听 | 亚洲精品制服丝袜二区| 欧美日韩国产在线人成app| 久久精品无码专区免费| 欧洲高清无码在线| 伊伊人成亚洲综合人网7777| av在线人妻熟妇| 亚洲无码高清视频在线观看| 高h视频在线| 98超碰在线观看| 91偷拍一区| 欧美成人综合在线| 午夜国产精品视频黄| 播五月综合| 国产在线自揄拍揄视频网站| 久久久亚洲色| av一区二区三区在线观看 | 国产成人高清精品免费| 亚洲人成影视在线观看| 91区国产福利在线观看午夜| 99热这里只有精品久久免费| 国内熟女少妇一线天| 毛片在线看网站| 日日摸夜夜爽无码| 亚洲综合香蕉| 美女扒开下面流白浆在线试听| 国产人成乱码视频免费观看| 亚洲第一区在线| 中文字幕亚洲电影| 国产黄网站在线观看| www中文字幕在线观看| 麻豆国产精品| 日韩小视频网站hq| 熟妇人妻无乱码中文字幕真矢织江 | 极品国产一区二区三区| 欧美激情首页| 天天做天天爱夜夜爽毛片毛片| 韩国v欧美v亚洲v日本v| 精品亚洲麻豆1区2区3区| 日本国产在线| 亚洲黄网视频| 色综合国产| 18禁高潮出水呻吟娇喘蜜芽| 婷婷亚洲视频| 无码内射中文字幕岛国片| 亚洲色婷婷一区二区| 亚洲人成网7777777国产| 制服丝袜亚洲| 国产成人超碰无码| 精品伊人久久久大香线蕉欧美| 日韩美毛片| 国产成人三级| 久久国产亚洲偷自| 国产综合欧美| 亚洲第一页在线观看| 欧美日韩国产在线播放| 国产一区二区三区在线观看视频 | 2020亚洲精品无码| 日韩在线2020专区| 亚洲大尺码专区影院| 国内a级毛片| 免费观看无遮挡www的小视频| 2022国产91精品久久久久久| 亚洲女同一区二区| 久久综合干| 亚洲人成人伊人成综合网无码| 日韩经典精品无码一区二区| 在线永久免费观看的毛片| 成人免费网站久久久| 国产成人精品一区二区免费看京| 国产一区亚洲一区|