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

        ?

        具有傳質(zhì)傳熱及擴(kuò)散效應(yīng)的雙氣泡的相互作用*

        2023-10-30 06:50:46烏日樂格那仁滿都拉
        物理學(xué)報(bào) 2023年19期
        關(guān)鍵詞:傳質(zhì)聲壓空化

        烏日樂格 那仁滿都拉

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

        1 引言

        在聲空化研究中關(guān)于單氣泡的研究相對比較深入,許多研究者針對不同的側(cè)面問題進(jìn)行了研究.如Hilgenfeldt 等[1]用流體力學(xué)方法研究氣泡形狀穩(wěn)定性和擴(kuò)散穩(wěn)定性,給出了氣泡聲致發(fā)光相圖.Yasui[2]在考慮水蒸氣的非平衡蒸發(fā)和冷凝、氣泡內(nèi)外熱傳導(dǎo)及泡內(nèi)化學(xué)反應(yīng)的基礎(chǔ)上,研究了單泡聲致發(fā)光機(jī)制.Toegel 和Lohse[3]在考慮水蒸氣相變、熱傳導(dǎo)及化學(xué)反應(yīng)的情況下,計(jì)算出了單泡聲致發(fā)光相圖.Shen 等[4]研究了空化氣泡壁處液體溫度的空間分布等.但實(shí)際液體中的空化泡總是以泡群的形式存在[5?7],而空化雙泡顯然是研究空化泡群現(xiàn)象的一個(gè)重要的切入點(diǎn),所以許多學(xué)者開展了對空化雙泡的研究.當(dāng)液體中的氣泡受到聲場作用時(shí),氣泡會發(fā)生體積脈動,氣泡之間會發(fā)生相互吸引或排斥作用,這與氣泡之間的次Bjerknes力有關(guān).Mettin 等[8]首先給出次Bjerknes 力的計(jì)算表達(dá)式,并數(shù)值研究了強(qiáng)聲場作用下雙氣泡間的次Bjerknes 力的變化.Doinikov[9]提出一種能夠直接計(jì)算兩個(gè)相互作用的球形氣泡在強(qiáng)聲場中平移運(yùn)動的模型,研究了兩個(gè)氣泡平移運(yùn)動的規(guī)律.Sadighi-Bonabi 等[10]研究不同濃度硫酸溶液對雙泡相互作用力的影響,并認(rèn)為隨著溶液黏度的增加次Bjerknes 力減小.Zhang 等[11]重點(diǎn)分析了雙頻聲激勵(lì)下兩個(gè)空化氣泡間相互作用的次Bjerknes力.Pandey[12]考慮雙泡間的強(qiáng)非線性耦合后對相互作用項(xiàng)進(jìn)行高階修正,并通過研究指出高階修正項(xiàng)對近距離氣泡間的次Bjerknes 力有重要影響.Luo 和Niu[13]通過激光和水下低壓放電技術(shù)分別誘導(dǎo)兩個(gè)空化泡,并對不同時(shí)期誘導(dǎo)空泡坍塌時(shí)的噴射流和沖擊波進(jìn)行實(shí)驗(yàn)測量,揭示了雙泡坍塌形成的噴射流和沖擊波的一些規(guī)律.Zhang 等[14]數(shù)值研究了聲場中兩個(gè)空化氣泡的脈動和平移,并指出氣泡脈動使氣泡發(fā)生平移,不同的脈動導(dǎo)致不同的平移.李想等[15]在Rayleigh-Plesset 方程的基礎(chǔ)上,引入可壓縮流體次Bjerknes 輻射力,建立了考慮管道軸向氣泡分布的可壓縮雙氣泡耦合動力學(xué)模型,并分析了次Bjerkness 輻射力對雙氣泡的線性和非線性動力學(xué)特征的影響.Shen 等[16]研究氣泡間的相互作用對氣泡徑向脈動的影響,并給出氣泡間的相互作用不僅可以降低或抑制氣泡的膨脹比,還可以增大氣泡的膨脹比.Qin 等[17]利用一種包括徑向和平移運(yùn)動的綜合模型,數(shù)值研究了黏彈性組織中兩個(gè)相互作用空化泡在不同驅(qū)動條件下的非線性動力學(xué)和聲輻射特性.Zhang 等[18]數(shù)值考察了強(qiáng)聲場中的大泡與小泡的平移運(yùn)動和非球面振蕩對次Bjerknes 力的影響.王德鑫和那仁滿都拉[19]研究了不同泡半徑、不同泡間距、不同惰性氣體對雙泡聲空化特性的影響.王尋等[20]用數(shù)值方法研究了方波驅(qū)動下雙氣泡的動力學(xué)行為.以上研究都未涉及到同時(shí)考慮傳質(zhì)傳熱及擴(kuò)散效應(yīng)的空化雙泡間相互作用的研究,但對于球形空化泡群的研究中分別涉及到了質(zhì)量傳輸、熱量傳遞和化學(xué)反應(yīng)等效應(yīng)[21?23].

        本文考慮空化雙泡的傳質(zhì)傳熱及擴(kuò)散效應(yīng),利用修正的Keller-Miksis 方程與Noble-Abel-Stiffend-Gas 狀態(tài)方程相結(jié)合的動力學(xué)模型,在不同半徑、不同驅(qū)動聲壓和不同驅(qū)動頻率情況下,研究聲空化氣泡的歸一化最大半徑和氣泡間的次Bjerknes 力的變化,并與Mettin 等[8]給出的結(jié)果進(jìn)行比較分析.

        2 考慮傳質(zhì)傳熱及擴(kuò)散效應(yīng)的雙泡模型

        Yasui[21]給出的修正的Keller-Miksis (KM)方程為

        式中,小黑點(diǎn)表示對時(shí)間導(dǎo)數(shù),Ri為第i個(gè)氣泡的瞬時(shí)半徑,為第i個(gè)氣泡的單位時(shí)間單位面積上的蒸發(fā)和冷凝速率,Cl,i為第i個(gè)氣泡壁處液體中的聲速,ρl,i(ρl,∞) 為第i個(gè)氣泡壁處液體(無窮遠(yuǎn)處)密度,ρg,i為第i個(gè)氣泡內(nèi)氣體密度,Pg,i為第i個(gè)氣 泡內(nèi) 氣體 壓力,Ps,i=-Asin[2πf(t+Ri/Cl,i)]為作用于第i個(gè)氣泡的驅(qū)動聲壓,σ為表面張 力,μ為黏性系數(shù),dij為兩氣泡間的間距,(1)式的最后一項(xiàng)是相互作用項(xiàng).

        氣泡內(nèi)氣體及周圍液體的熱力學(xué)狀態(tài),采用Métayer 和Saurel[24]給出的Noble-Abel-Stiffened-Gas (NASG)狀態(tài)方程來描述:

        式中,P為壓力;γ為熱容比;CV為定容熱容;T為溫度;b為分子共體積;B為壓力常數(shù);v為比容,與密度的關(guān)系為v=1/ρ.根據(jù)聲速定義和狀態(tài)方程(3),第i個(gè)氣泡壁處液體中的聲速Cl,i可表示為

        通過狀態(tài)方程(3),可得到第i個(gè)氣泡壁處液體密度為

        第i個(gè)氣泡內(nèi)氣體和水蒸氣的壓力為

        式中,n=Ar,H2O 分別表示泡內(nèi)氬氣和水蒸氣(本文假設(shè)泡內(nèi)含有氬氣和水蒸氣),Nn,i為第i個(gè)氣泡內(nèi)氣體分子數(shù),NA是阿伏伽德羅常數(shù),Mn,i為第i個(gè)氣泡內(nèi)氣體摩爾質(zhì)量,Tg,i為第i個(gè)氣泡內(nèi)氣體溫度,ΔCn,i為第i個(gè)氣泡內(nèi)氣體的定壓熱容和定容熱容之差,定壓熱容和定容熱容可用文獻(xiàn)[25]給出的表達(dá)式計(jì)算.bl和Bl分別為液體分子共體積和壓力常數(shù),bg,n和Bg,n分別為氣體分子共體積和壓力常數(shù).對于氬氣和水蒸氣,利用文獻(xiàn)[26,27]給出的公式,計(jì)算出不同溫度和壓力下的bg,Ar和值,并做擬合處理后得到如下表達(dá)式:

        第i個(gè)空化氣泡內(nèi)單位時(shí)間單位面積上的蒸發(fā)和冷凝速率,可由氣體動力學(xué)理論導(dǎo)出的Hertz-Knudsen 公式計(jì)算[28]:

        式中,常數(shù)α=0.4,為分子平均速 度,為水蒸氣密度,飽和水蒸氣密度=0.0173 kg·m-3.因此,水蒸氣的蒸發(fā)和冷凝引起的分子數(shù)的變化可表示為

        氣體從第i個(gè)氣泡內(nèi)部擴(kuò)散到周圍液體所引起的氣體分子數(shù)的變化,可由Epstein-Plesset 理論中所用的氣體擴(kuò)散方程來表示[29]:

        式中,DAr為氬氣的擴(kuò)散系數(shù),cs為氣體的平衡濃度,c∞為氣體的飽和濃度.

        第i個(gè)氣泡內(nèi)外的熱交換量由下式計(jì)算[28]:

        式中,lth為熱邊界層厚度.χ為邊界層的熱擴(kuò)散率,λ為混合氣體的熱導(dǎo)率.

        由于氣泡的內(nèi)能是溫度和體積的函數(shù)[29]

        根據(jù)恒定體積下摩爾熱容的定義,下列關(guān)系成立:

        對于 NASG 狀態(tài)方程來講,(14)式中的第二項(xiàng)變?yōu)?/p>

        把方程(15)和(16)代入(14)式并進(jìn)行積分,就可得到第i個(gè)氣泡內(nèi)氣體溫度的計(jì)算公式:

        其中,Ei為第i個(gè)氣泡內(nèi)氣體內(nèi)能,和CV,Ar(T)分別為水蒸氣和氬氣的定容摩爾熱容[25].因?yàn)榕輧?nèi)溫度的計(jì)算需要內(nèi)能的表達(dá)式,根據(jù)熱力學(xué)第一定律,內(nèi)能變化表示為

        式中,等式右端第一項(xiàng)表示第i個(gè)氣泡內(nèi)壓力做功所引起的能量變化;第二項(xiàng)表示第i個(gè)氣泡周圍液體的水分子蒸發(fā)成水蒸氣和氣泡內(nèi)的水蒸氣分子凝結(jié)到液體中所攜帶的能量變化;第三項(xiàng)表示熱傳導(dǎo)引起的能量變化;第四項(xiàng)表示第i個(gè)氣泡內(nèi)氣體分子擴(kuò)散所引起的能量變化.

        兩個(gè)球形空化氣泡之間的次Bjerknes 力為[8]

        式中,〈·〉為時(shí)間平均值,和為氣泡1 和氣泡2的體積變化率.另外,為方便研究問題,定義了一個(gè)新的物理量,即次Bjerknes 力系數(shù),表示為[8]

        式中,fB的符號表明雙泡間是吸引還是排斥.當(dāng)fB>0 時(shí),說明雙泡相互吸引;當(dāng)fB<0 時(shí),說明雙泡相互排斥.

        另外,為了與沒有傳質(zhì)傳熱及擴(kuò)散效應(yīng)的雙空化氣泡間的相互作用比較,本文利用Mettin 等[8]給出的雙泡模型來計(jì)算沒有傳質(zhì)傳熱及擴(kuò)散效應(yīng)的雙氣泡間相互作用的次Bjerknes 力.

        3 數(shù)值模擬

        本文選擇水為液體介質(zhì),氣泡內(nèi)氣體為氬氣和水蒸氣.為了方便描述,下面用模型1 代表具有傳質(zhì)傳熱及擴(kuò)散效應(yīng)的KM-NASG 雙泡模型,用模型2 代表Mettin 等[8]給出的沒有傳質(zhì)傳熱及擴(kuò)散效應(yīng)的雙泡模型.主要對比分析兩種模型給出的空化氣泡的歸一化最大半徑以及氣泡間的次Bjerknes 力在一個(gè)聲周期內(nèi)的變化情況.用兩種模型計(jì)算時(shí)所用的相關(guān)參數(shù)見表1 和表2,這些參數(shù)的選取可參考文獻(xiàn)[8,28,30].

        表1 模型1 的相關(guān)物理參數(shù)Table 1.Related physical parameters of model 1.

        表2 模型2 的相關(guān)物理參數(shù)Table 2.Related physical parameters of model 2.

        3.1 總分子數(shù)的變化

        為了解釋說明有傳質(zhì)傳熱和擴(kuò)散效應(yīng)的雙氣泡間的相互作用與沒有傳質(zhì)傳熱和擴(kuò)散效應(yīng)的雙氣泡間的相互作用的差別,首先需要計(jì)算氣泡內(nèi)氣體總分子數(shù)NT隨氣泡半徑、驅(qū)動聲壓和驅(qū)動頻率的變化.由于空化氣泡內(nèi)的氣體分子數(shù)量與蒸發(fā)和冷凝以及擴(kuò)散等兩種效應(yīng)有關(guān),所以本文把兩種效應(yīng)都考慮后的泡內(nèi)含有的氣體分子數(shù)量稱為總分子數(shù).另外,氣泡內(nèi)氣體總分子數(shù)是隨著氣泡的膨脹、收縮崩潰與回彈階段而不斷變化[31],這里計(jì)算的是最大總分子數(shù)NT,max,即氣泡膨脹到最大半徑時(shí)氣泡內(nèi)所含有的總分子數(shù).

        圖1 空化氣泡內(nèi)氣體最大總分子數(shù) 的變化 (a)當(dāng) R20=5 μm,A=1.32P0,f0=20 kHz 時(shí),NT,max 隨氣泡1 初始半徑的變化;(b)當(dāng) R10=R20=5 μm,f0=20 kHz 時(shí), 隨驅(qū)動聲壓的變化;(c)當(dāng) R10=R20=5 μm,A=1.32P0時(shí),隨驅(qū)動頻率的變化Fig.1.Change of the maximum total molecular number of gas in cavitation bubble: (a) changes with the initial radius of bubble 1 at R20=5 μm,A=1.32P0,f0=20 kHz ;(b) changes with driving sound pressure at R10=R20=5 μm,f0=20 kHz ;(c) NT,max changes with driving frequency at R10=R20=5 μm,A=1.32P0.

        3.2 初始半徑對次Bjerknes 力的影響

        圖2 是驅(qū)動聲壓幅值A(chǔ)=1.32P0,驅(qū)動頻率f0=20 kHz時(shí),在氣泡2 的初始半徑恒定,氣泡1 的初始半徑變化的情況下,用兩種模型繪制的氣泡歸一化最大半徑和次Bjerknes 力系數(shù)的變化.從圖2(a)可以看出,在R20(=5 μm)恒定時(shí),隨著R10的增大,兩種模型給出的氣泡歸一化最大半徑的變化趨勢是一致的.比較發(fā)現(xiàn),模型1 給出的氣泡歸一化最大半徑始終大于模型2 給出的氣泡歸一化最大半徑.這是因?yàn)槟P? 考慮了氣泡的傳質(zhì)傳熱和擴(kuò)散效應(yīng),所以當(dāng)R10增大時(shí),氣泡內(nèi)的最大總分子數(shù)會明顯增大(圖1(a)),也就是氣泡內(nèi)氣體含量增多,導(dǎo)致模型1 給出的氣泡歸一化最大半徑較大.從圖2(b)可以看出,隨著R10的增大,兩種模型給出的次Bjerknes 力系數(shù)逐漸增大,顯然模型1 給出的次Bjerknes 力系數(shù)大于模型2 給出的次Bjerknes 力系數(shù),且差別逐漸增大.這主要是因?yàn)殡S著R10的增大,模型1 給出的氣泡歸一化最大半徑始終大于模型2 給出的歸一化最大半徑(圖2(a)).由于Bjerknes 力特別敏感依賴于空化氣泡半徑及其變化率[8],所以模型1 給出的次Bjerknes 力系數(shù)大于模型2 給出的值.由圖2(b)插圖可以看出,在到Blake 空化閾值半徑之前fB小于0,表明氣泡間次Bjerknes 力為斥力;到Blake 空化閾值半徑之后fB迅速增大為大于0,表明氣泡間次Bjerknes 力變?yōu)橐?同時(shí)也可以看出,模型1 給出的次Bjerknes 力由斥力變?yōu)橐Φ倪M(jìn)程提前了一些.

        圖2 空化氣泡歸一化最大半徑和次Bjerknes 力系數(shù)隨氣泡1 初始半徑的變化 (a)歸一化最大半徑;(b)次Bjerknes 力系數(shù)Fig.2.Maximum normalized radius of cavitation the bubbles and the secondary Bjerknes force coefficient change with the initial radius of bubble 1: (a) Maximum normalized radius;(b) the secondary Bjerknes force coefficient.

        3.3 驅(qū)動聲壓對次Bjerknes 力的影響

        圖3 是雙泡初始半徑R10=R20=5 μm,驅(qū)動頻率f0=20 kHz 時(shí),繪制的氣泡歸一化最大半徑和次Bjerknes 力系數(shù)隨驅(qū)動聲壓幅值的變化圖.從圖3(a)可以看出,隨著驅(qū)動聲壓幅值的增大,兩種模型給出的氣泡歸一化最大半徑呈現(xiàn)逐漸增大的趨勢.比較發(fā)現(xiàn),模型1 給出的氣泡歸一化最大半徑大于模型2 給出的值.這是因?yàn)槟P? 考慮了氣泡的傳質(zhì)傳熱和擴(kuò)散效應(yīng),所以當(dāng)驅(qū)動聲壓幅值增大時(shí),氣泡內(nèi)的最大總分子數(shù)會增多(圖1(b)),也就是氣泡內(nèi)氣體含量增多,導(dǎo)致模型1 給出的氣泡歸一化最大半徑較大.從圖3(b)可以看出,隨著驅(qū)動聲壓幅值的增大,兩種模型給出的次Bjerknes力系數(shù)逐漸增大,顯然模型1 給出的次Bjerknes力系數(shù)大于模型2 給出的值,且差別逐漸增大.這是因?yàn)殡S著驅(qū)動聲壓幅值的增大,模型1 給出的氣泡歸一化最大半徑逐漸增大且始終大于模型2 給出的值(圖3(a)).由于Bjerknes 力特別敏感依賴于空化氣泡半徑及其變化率,所以模型1 給出的次Bjerknes 力系數(shù)大于模型2 給出的力系數(shù).總體上,隨著驅(qū)動聲壓幅值的增大,模型1 給出的氣泡歸一化最大半徑和次Bjerknes 力系數(shù)大于模型2 給出的相應(yīng)值.

        圖3 空化氣泡歸一化最大半徑和次Bjerknes 力系數(shù)隨驅(qū)動聲壓的變化 (a)歸一化最大半徑;(b)次Bjerknes 力系數(shù)Fig.3.Maximum normalized radius of cavitation the bubbles and the secondary Bjerknes force coefficient change with driving sound pressure: (a) Maximum normalized radius;(b) the secondary Bjerknes force coefficient.

        3.4 驅(qū)動頻率對次Bjerknes 力的影響

        圖4 給出的是雙泡初始半徑為R10=R20=5 μm,驅(qū)動聲壓幅值為A=1.32P0時(shí),用兩種模型繪制的氣泡歸一化最大半徑和次Bjerknes 力系數(shù)隨驅(qū)動頻率的變化圖.從圖4(a)可以看出,隨著驅(qū)動頻率的增大,兩種模型給出的氣泡歸一化最大半徑逐漸減小.比較可看到,模型1 給出的氣泡歸一化最大半徑大于模型2 給出的值.其原因是模型1 考慮了氣泡的傳質(zhì)傳熱和擴(kuò)散效應(yīng),所以當(dāng)驅(qū)動頻率增大時(shí),雖然氣泡內(nèi)的氣體最大總分子數(shù)逐漸減少,但始終大于模型2 給出的值(圖1(c)),所以模型1 給出的氣泡歸一化最大半徑較大.從圖4(b)可以看出,隨著驅(qū)動頻率的增大,兩種模型給出的次Bjerknes 力系數(shù)也逐漸減小.相比可見,模型1 給出的次Bjerknes 力系數(shù)大于模型2 給出的力系數(shù).其原因主要還是模型1 給出的氣泡最大半徑及半徑的變化率始終大于模型2 給出的相應(yīng)值.總體上,隨著驅(qū)動頻率的增大,模型1 給出的次Bjerknes 力系數(shù)始終大于模型2 給出的次Bjerknes力系數(shù).

        圖4 空化氣泡歸一化最大半徑和次Bjerknes 力系數(shù)隨驅(qū)動頻率的變化 (a) 歸一化最大半徑;(b) 次Bjerknes 力系數(shù)Fig.4.Maximum normalized radius of cavitation the bubbles and the secondary Bjerknes force coefficient change with driving frequency: (a) Maximum normalized radius;(b) the secondary Bjerknes force coefficient.

        4 結(jié)論

        本文利用KM-NASG 模型,研究了具有傳質(zhì)傳熱及擴(kuò)散效應(yīng)的聲空化雙泡間的相互作用,分別計(jì)算了氣泡歸一化最大半徑和氣泡間相互作用的次Bjerknes 力(系數(shù)),并與沒有傳質(zhì)傳熱及擴(kuò)散效應(yīng)的雙泡系統(tǒng)中的氣泡歸一化最大半徑和次Bjerknes 力(系數(shù))進(jìn)行了比較.結(jié)果表明,相比于沒有傳質(zhì)傳熱及擴(kuò)散效應(yīng)的雙泡系統(tǒng),具有傳質(zhì)傳熱及擴(kuò)散效應(yīng)的雙泡系統(tǒng)中的氣泡歸一化最大半徑和氣泡間相互作用的次Bjerknes 力都要大.隨著雙泡系統(tǒng)中某一氣泡的初始半徑的增大,次Bjerknes 力會增大,兩種情況下的次Bjerknes 力的差別會逐漸增大;隨著驅(qū)動聲壓的增大,次Bjerknes 力也會增大,兩種情況下的次Bjerknes力的差別也會逐漸增大;隨著驅(qū)動頻率的增大,次Bjerknes 力會減小,兩種情況下的次Bjerknes 力的差別也會逐漸減小.因此傳質(zhì)傳熱和擴(kuò)散效應(yīng)對空化氣泡間相互作用的次Bjerknes 力有很明顯的影響,是不能忽略的重要因素.

        猜你喜歡
        傳質(zhì)聲壓空化
        功率超聲作用下鋼液中空化泡尺寸的演變特性
        鋼鐵釩鈦(2023年5期)2023-11-17 08:48:34
        基于嘴唇處的聲壓數(shù)據(jù)確定人體聲道半徑
        車輛結(jié)構(gòu)噪聲傳遞特性及其峰值噪聲成因的分析
        汽車工程(2018年12期)2019-01-29 06:46:36
        三維扭曲水翼空化現(xiàn)象CFD模擬
        氨基酸鹽吸收二氧化碳過程的傳質(zhì)特性
        不同運(yùn)動形式下水物相互作用空化數(shù)值模擬
        基于GIS內(nèi)部放電聲壓特性進(jìn)行閃絡(luò)定位的研究
        電測與儀表(2016年9期)2016-04-12 00:30:02
        PTFE膜吸收CO2的工藝參數(shù)對傳質(zhì)性能的影響
        清潔轉(zhuǎn)向酸H+表面?zhèn)髻|(zhì)行為實(shí)驗(yàn)研究
        煤顆粒熱解的傳熱傳質(zhì)分析
        亚洲av资源网站手机在线 | 人人妻人人澡人人爽国产| 肉丝高跟国产精品啪啪| 国内精品大秀视频日韩精品| 少妇被粗大的猛烈进出69影院一| 大香蕉国产av一区二区三区| 国产成人自拍视频在线免费| 亚洲中文字幕无码中字| 十八禁在线观看视频播放免费 | 一区二区三区无码高清视频| 日本一区二区高清视频| 国产精品国产三级国产专区5o| 免费看泡妞视频app| 99在线精品免费视频| 国产乱淫h侵犯在线观看| 亚洲va中文字幕欧美不卡| 男人无码视频在线观看| 蜜桃久久精品成人无码av| 东北女人啪啪对白| 亚洲啪啪色婷婷一区二区| 亚洲高清在线观看免费视频| 久久香蕉成人免费大片| 无码精品日韩中文字幕| 日本久久久久亚洲中字幕| 一区二区三区四区在线观看日本| 国产女主播在线免费观看| 国产三级黄色在线观看| 熟妇人妻中文字幕无码老熟妇| 中文字幕亚洲乱码熟女在线| 国产在线无码不卡影视影院| 美艳善良的丝袜高跟美腿| 亚洲中文字幕乱码在线视频| 人妻无码人妻有码不卡| 日韩久久一级毛片| 理论片87福利理论电影| 极品少妇一区二区三区四区| 欧洲乱码伦视频免费| 99久久免费看精品国产一| 草青青在线视频免费观看| 亚洲一区二区三区99区| 在线国产视频精品视频|