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

豎直振動激勵下顆粒毛細上升行為研究*

2022-06-04 06:25:36于天林凡鳳仙2
物理學報 2022年10期
關鍵詞:振動

于天林 凡鳳仙2)?

1) (上海理工大學能源與動力工程學院,上海 200093)

2) (上海理工大學上海市動力工程多相流動與傳熱重點實驗室,上海 200093)

豎直振動激勵下顆粒毛細上升現象為顆粒物料的提升、輸運和采集提供了一種全新的技術路徑.然而,已有顆粒毛細上升行為研究仍存在明顯不足,特別是缺乏對重力加速度、水平振動分量、顆粒粒徑分布影響的深入探究.針對這些問題,采用離散元方法,對不同操作條件下顆粒毛細上升現象開展數值模擬研究,并對顆粒最終毛細上升高度和平均毛細上升速度進行計算分析.結果表明,在低重力條件下,顆粒毛細效應仍能發生,顆粒最終毛細上升高度和平均毛細上升速度隨重力加速度均呈現先增加后減小的趨勢;顆粒最終毛細上升高度對水平振動分量的變化不敏感,而平均毛細上升速度隨水平振動分量的增大而增加;在平均粒徑相等的情況下,粒徑服從高斯分布的顆粒比單一粒徑顆粒的最終毛細上升高度最大值對應的臨界管徑更大,并且同處于堵塞效應影響的管徑區域時平均毛細上升速度也更快.

1 引言

顆粒物質是指粒徑大于1 μm 的大量離散固體顆粒組成的宏觀復雜體系,其具有對外界微小作用的敏感性、非線性響應和自組織行為等[1,2].特別是,豎直振動激勵下顆粒物質能夠表現出復雜而奇特的動力學行為,例如:豎直振動顆粒床中顆粒的隆起和對流[3-5];豎直振動U 形分割容器或U 形管中顆粒由一側向另一側的遷移[6-8];豎直振動激勵插入靜止顆粒床中的細管時顆粒毛細上升行為[9-11]等.這些行為難以通過傳統的流體力學、固體力學以及凝聚態物理理論進行解釋,因而成為當前工程熱物理、離散相動力學、凝聚態物理等多學科交叉的前沿領域和研究熱點.

對插入靜止顆粒床中的細管施加豎直方向正弦振動激勵作用,在一定的振動強度下,顆粒在管內逆重力上升,并最終形成一個穩定顆粒柱高度,鑒于這一現象與與液體的毛細效應類似,被稱為顆粒毛細效應[11].該現象為顆粒物料的提升、輸運和采集提供了一種全新思路,具有工程應用的潛力,因而備受關注.近年來,一些學者對此開展了一系列實驗觀測、理論分析和離散元模擬研究,探討了振動條件、幾何參數、顆粒粒徑對顆粒毛細上升高度和速度的影響[9-23],分析了顆粒毛細上升的“空穴填充”機理[21],發現了顆粒毛細上升的內在物理機制—顆粒對流[11-14],提出了顆粒毛細上升高度的理論模型[23]等.

然而,由于顆粒毛細上升過程機理復雜、影響因素繁多,在顆粒毛細上升行為研究方面仍存在明顯不足.首先,已有研究集中在開展地面實驗或在地面重力加速度下開展數值模擬,考察無量綱振動強度(最大振動加速度與重力加速度的比值)對顆粒毛細上升的影響規律,低重力條件下顆粒毛細效應的文獻報道僅見Kawamoto 等[20]在重力加速度g=1.6 m/s2情況下數值模擬得到的顆粒毛細上升高度隨時間的演變,尚缺少顆粒毛細上升行為受重力加速度影響的深入探究.其次,實際的豎直振動系統中,由于擾動的存在,微小的水平振動分量難以避免,已有對豎直振動U 形管中顆粒遷移行為的研究表明,水平振動分量的存在有利于顆粒由一個分支向另一個分支的遷移[8],然而水平振動分量對顆粒毛細上升行為的影響特性尚未被研究.再者,已有顆粒毛細上升行為研究采用的顆粒均為單分散顆粒,多分散顆粒的毛細上升行為有待研究.

鑒于已有研究中存在的上述問題,本文基于離散元方法,對豎直振動激勵下顆粒毛細上升行為開展數值模擬,對比分析低重力與地面重力條件下顆粒毛細上升過程,揭示顆粒最終毛細上升高度與平均毛細上升速度隨重力加速度、水平振動分量的變化規律,考察粒徑分布對顆粒毛細上升動力學的影響特性,以期為振動體系下顆粒的動力學行為研究提供參考.

2 離散元模型與數值計算方法

2.1 離散元模型

離散元模型的基本思想是:將顆粒物質視為離散顆粒的集合體,每個離散顆粒的運動都服從牛頓第二定律,通過求解各顆粒的運動過程,獲得顆粒系統的整體行為規律.針對顆粒毛細上升行為發生的系統,考慮顆粒間接觸力、顆粒受到的重力、切向接觸力產生的力矩和滾動摩擦力產生的力矩,建立描述顆粒運動的離散元模型.顆粒系統中任意顆粒i的運動方程可寫為:

式中,mi為顆粒i的質量,vi為顆粒i的速度,t為時間,N為與顆粒i相接觸顆粒的數目,Fn,ij為與顆粒i相接觸的顆粒j對顆粒i的法向作用力,Ft,ij為顆粒j對顆粒i的切向作用力,g為重力加速度,Ii為顆粒i的轉動慣量,ωi為顆粒i的角速度,Mt,ij為顆粒j對顆粒i的切向作用力產生的力矩,Mr,ij為顆粒j對顆粒i的滾動摩擦力產生的力矩.

分別采用黏彈性接觸模型與修正的Cundall-Strack 模型計算顆粒間的法向作用力[24]與切向作用力[25],依據定量恒轉矩模型計算滾動摩擦力產生的力矩[26],接觸力和力矩的表達式詳見文獻[12,13],此處不再贅述.在本文所考慮的顆粒系統中,除了顆粒間相互作用外,還存在顆粒與容器壁及管壁的相互作用.對于顆粒i與壁面的相互作用,將壁面視為半徑無限大的顆粒處理.

2.2 數值計算方法

基于上述離散元模型,借助開源顆粒系統數值模擬軟件LIGGGHTS[27]對豎直振動激勵下顆粒毛細上升行為開展數值模擬研究.圖1 給出了數值模擬采用的顆粒系統示意圖.如圖所示,內徑為D的圓柱形容器中填充有大量顆粒,將內徑為d、壁厚為δ的圓管沿容器的中心軸線豎直插入顆粒層,直至管底面與容器底面的距離為L.其中,顆粒層利用數值模擬方法生成,具體過程為:在容器內(包括管外和管內區域)隨機加入初速度為0 的模擬顆粒,使顆粒遵循(1)式、(2)式給出的運動方程發生沉降運動,顆粒系統的總動能先增大后減小,并最終形成總動能為0 的堆積狀態.

圖1 顆粒系統示意圖Fig.1.Schematic diagram of granular system.

以管的底面中心為坐標原點,建立三維坐標系.沿豎直方向(z向)對管施加振幅為Az、頻率為f的正弦振動,則管底面的位置坐標zb隨時間t的變化可寫為

當考察水平振動分量對顆粒毛細上升動力學行為的影響時,除上述豎直振動外,還考慮了振幅為Ax、頻率為f的水平方向(x向)正弦振動.由于該水平振動,管中心軸線的位置坐標xa隨時間t的變化可寫為

在計算模型和顆粒系統確定的情況下,離散元模擬的準確性與經濟性主要受時間步長的影響.為確定時間步長,首先利用無阻尼、無黏性碰撞時間估算顆粒接觸時間,即[28]

式中,meff=mimj/(mi+mj)為碰撞顆粒對i和j的有效質量,vimp為顆粒碰撞速度;ρ為彈性參數,其可表示為

式中,Y為楊氏模量,υ為泊松比,Ri與Rj分別為顆粒i與j的半徑.

數值模擬采用的參數見表1.在顆粒直徑dp=0.6 mm 條件下,取vimp=1 m/s 作為顆粒碰撞速度的參考值,結合表1 給出的顆粒物性參數,計算得到tc≈ 6×10—5s.因此,選取時間步長Δt=tc/30 ≈ 2×10—6s[29,30].研究發現,若采用更小的時間步長,計算結果不發生改變,這表明時間步長的選擇是合理的.此外,在該時間步長下,每計算2500 步保存一次計算結果,即一個振動周期T(T=1/f)內以等時間間隔保存20 次計算結果.

表1 數值模擬參數Table 1.Parameters used in numerical simulations.

3 結果與討論

3.1 不同重力加速度下顆粒毛細上升行為

在不考慮水平振動分量(Ax=0)情況下,采用表1 給出的參數,針對直徑dp=0.6 mm 的單一粒徑顆粒開展數值模擬研究,考察重力加速度在2.0—9.8 m/s2范圍變化時顆粒的毛細上升過程.圖2 給出了重力加速度g=2.0 m/s2和g=9.8 m/s2時顆粒毛細上升動態過程的對比情況.圖中,依據施加振動時刻(t=0 時)顆粒的位置對顆粒進行著色,位置矢量r滿足|r|

圖2 不同重力加速度下顆粒毛細上升動力學過程Fig.2.Dynamic process of granular capillary rising under different gravitational accelerations.

為定量表征顆粒毛細上升現象,引入毛細上升高度Hc、最終毛細上升高度和平均毛細上升速度Va三個參數.顆粒毛細上升高度定義為管內顆粒柱高度,即管內、管外(容器內)顆粒表面高度之差.獲得Hc的方法為:將管內z> 0 的區域沿z向劃分為高度Δz=5dp的連續單元,計算各單元內顆粒的填充率φ,將滿足φ(z) > 0.1 的單元的最大高度作為管內顆粒表面高度zin;類似地,將容器內z> 0 且|r| >d/2 的區域沿z向劃分為Δz=dp的連續單元,采用與確定管內顆粒表面高度相同的方法,獲得管外顆粒表面高度zout,基于此,Hc=zin—zout.鑒于管內顆粒柱高度隨著管的振動而發生周期性波動,取該階段10T時間內Hc的平均值作為顆粒最終毛細上升高度,考慮到2.2 節給出的計算結果保存情況,則有

式中,t1為穩定階段的起始時間,tk=t1+(k-1)T/20.

鑒于顆粒毛細效應動力學過程中由上升階段過渡到穩定階段是一個顆粒柱高度增加速度緩慢衰減的過程,給準確確定上升時間帶來困難.為解決這一問題,選取毛細上升高度首次達到0.95 倍最終上升高度(0.95)的時刻作為上升階段的截止時刻,用ta表示,如圖3 所示,則平均毛細上升速度的表達式為

圖3 毛細上升高度隨時間的演變Fig.3.Evolution of capillary rising height with time.

圖3 給出了g=2.0 m/s2,9.8 m/s2時顆粒毛細上升高度隨時間的演變情況.總體上來說,顆粒毛細上升高度均呈現出先迅速增大,而后增速降低,直到增速約為0 時顆粒毛細上升高度達到穩定的變化趨勢.然而,顆粒毛細上升高度隨時間演變的細節信息有所區別.低重力條件下,顆粒毛細上升高度存在幅度不一的起伏,且這種起伏沒有明顯的周期性;而地面重力下,顆粒毛細上升高度隨管的振動發生周期性波動,除去受初始效應影響的起始段(t=0—2T),該波動先與管的振動同頻(t=2T—10T),而后頻率減半(t> 10T),同時波動幅度顯著增加.這是因為:低重力條件下,管內顆粒堆積疏松(見圖2(a)),以至某些時段管內顆粒柱上部存在填充率φ<0.1 的稀疏區域,無法納入顆粒柱高度,導致顆粒毛細上升高度隨時間的演變呈現無序性;而地面重力條件下,管內顆粒堆積致密,除管內顆粒柱上部個別顆粒外,其余顆粒被納入顆粒柱高度的計算范圍,從而管內顆粒在管壁的周期性剪切作用下發生上升和下降運動,引起毛細上升高度發生周期性波動,此外在顆粒毛細上升過程中,顆粒間、顆粒與壁面間的碰撞和摩擦作用引起管內顆粒的填充率分布和堆積結構的改變,影響了顆粒的自組織行為,使得顆粒毛細上升高度波動呈現與管同頻向頻率減半轉變的現象.

圖4 給出了顆粒最終毛細上升高度與平均毛細上升速度隨重力加速度的變化關系.圖中,誤差條表示顆粒最終毛細上升高度的標準差.由圖可見,重力加速度在2.0—9.8 m/s2范圍內變化時,顆粒最終毛細上升高度和平均毛細上升速度均呈現出先增加后減小的趨勢,并均在g=6 m/s2時達到極大值,此時,=253.2 mm、Va=88.8 mm/s.Liu 和Zhao[21]根據振動輸入的能量與顆粒柱的重力勢能之間的平衡關系,推導出最終毛細上升高度與重力加速度及顆粒填充率負相關,而與能量輸入系數(顆粒柱獲得的能量與管振動動能之比)正相關.在重力加速度較小時,重力加速度增加,使得填充率增加,管內顆粒與管壁的作用增強,極大地促進了能量輸入系數的增加,顆粒最終毛細上升高度增加.在重力加速度較大時,填充率和能量輸入系數受重力加速度的影響均很小,因而隨著重力加速度增加,顆粒最終毛細上升高度降低.綜上所述,顆粒最終毛細上升高度呈現出隨重力加速度的增加而先增大后減小的趨勢.此外,管向上運動時,管內顆粒在管壁向上的剪切力作用下克服重力上升,在管底留下一個空穴,該空穴被容器內顆粒對流輸運的顆粒所填充;當管達到最大高度轉而向下運動時,顆粒在管壁向下的剪切力和重力的共同作用下減速上升,當顆粒速度減小至0 后,開始向下做加速運動,由于管內顆粒與容器內顆粒的碰撞引起的能量耗散,管內顆粒速度迅速減小為0;在管繼續運動至最低位置的過程中,填充到空穴的顆粒進入管內,并在管向上加速運動時受壁面剪切力作用而上升.重力加速度較低時,容器內顆粒的對流速度較小,不利于顆粒由容器內向管內的輸運;重力加速度較高時,管內顆粒重力較大,不利于填充到空穴的顆粒在管內的上升運動.因此,在管的振幅和頻率不變時,存在一個使顆粒平均毛細上升速度最大的重力加速度值.此外,由于重力加速度很低時,管向上加速運動時管內顆粒極易被提升至很高的位置,而在管向下加速運動時這些顆粒難以發生重力沉降,導致管內顆粒無法形成連續的顆粒柱,因而本文對于極低重力加速度(g<1 m/s2)下顆粒的毛細上升行為未進行探討.需要注意的是,鑒于空氣的存在并非是顆粒毛細上升的決定因素,本文的數值模擬中沒有考慮空氣的影響.張富翁等[10]在接近真空(—95 kPa)下和常壓下分別實驗觀測了平均粒徑為 0.3 mm的顆粒的毛細上升現象,發現接近真空下顆粒最終毛細上升高度比常壓下有所降低,二者之間的差在7%—13%范圍內.在今后的研究中,為精確模擬顆粒毛細上升動力學行為,需要在建模時考慮顆粒與空氣間的相互作用.

圖4 最終毛細上升高度和平均毛細上升速度隨重力加速度的變化Fig.4.Dependence of final capillary rising height and average capillary rising velocity on gravitational acceleration.

3.2 不同水平振動分量下顆粒毛細上升行為

在豎直振動的基礎上引入不同振幅的水平振動,在地面重力加速度(g=9.8 m/s2)條件下,針對直徑dp=0.6 mm 的單一粒徑顆粒毛細上升行為開展數值模擬.圖5 給出了顆粒最終毛細上升高度與平均毛細上升速度隨水平振動分量的變化特性.由圖可知,顆粒最終毛細上升高度對水平振動振幅的變化并不敏感,而顆粒平均毛細上升速度隨水平振動振幅的增大而升高.前者的原因是顆粒最終毛細上升高度由管內顆粒柱的重力勢能與管的振動輸入到顆粒柱的動能平衡關系決定,而地面重力加速度下,管內顆粒堆積密集,水平振動對顆粒柱動能輸入的影響較小,使得顆粒柱高度變化較小.后者與Sánchez 等[8]發現的豎直振動U 形管中顆粒遷移運動的規律類似.Sánchez 等[8]實驗研究發現,存在水平振動時,顆粒物質由U 形管一個分支向另一個分支遷移的速度隨著水平振動分量的增大而升高,但沒有對其原因進行解釋.顆粒柱上升階段管底區域顆粒平均填充率的計算結果表明,隨著水平振幅的增加,顆粒填充率趨于下降.基于此,圖5 所示的平均毛細上升速度變化規律可理解為:水平振動振幅的增加使得管口附近顆粒堆積更為松散,有利于容器內顆粒的對流運動,使得輸運到容器中心區域,進而進入管內的顆粒質量流率增大,顆粒平均毛細上升速度增加.

圖5 最終毛細上升高度和平均毛細上升速度隨水平振動分量的變化Fig.5.Dependence of final capillary rising height and average capillary rising velocity on horizontal vibrational component.

3.3 不同粒徑分布下顆粒毛細上升行為

在僅有豎直振動(Ax=0)時,對地面重力加速度下粒徑服從高斯分布的顆粒的毛細上升行為開展數值模擬,并將數值模擬結果與平均粒徑相等的單一粒徑顆粒的毛細上升行為進行對比分析.高斯分布函數可寫為

式中,p為頻率密度,其是顆粒直徑dp的函數;和σ分別為顆粒直徑的平均值和標準差.模擬中采用=0.6 mm,σ=0.1.

圖6 給出了兩種粒徑分布情況下顆粒最終毛細上升高度和平均毛細上升速度隨管徑的變化規律.由圖可見,兩種粒徑分布情況下顆粒最終毛細上升高度均隨管徑的變化呈現先增大后急劇減小的趨勢,因而存在一個使得顆粒最終毛細上升高度達到最大值的管徑,這與張華騰等[12]采用方形橫截面容器得到的結果定性一致.由圖還可知,與單一粒徑顆粒相比,高斯粒徑分布條件下顆粒最終毛細上升高度最大值及其對應的管徑值均較大.為表述方便,將顆粒最終毛細上升高度最大值對應的管徑稱為臨界管徑,并用dcr表示,可見兩種粒徑分布下dcr分別為8 mm 和10 mm.分析可知,當d≤dcr時,管內顆粒的動力學行為受到堵塞效應的影響,隨著管徑的增加,堵塞效應減弱,有利于顆粒的上升運動,使得此管徑段內最終毛細上升高度隨管徑的增加而增大;當d>dcr時,堵塞效應的影響消失,管內顆粒達到流態化,顆粒呈現出液體的行為,使得此管徑段內最終毛細上升高度呈現出與d≤dcr時不同的規律.本文重在關注顆粒粒徑分布對毛細效應動力學行為的影響,采用的管徑范圍有限,然而前期研究已充分說明,當d>dcr時,顆粒最終毛細上升高度隨管徑的增大而趨于減小,這與液體的毛細效應一致[11,12].由于顆粒毛細效應動力學的控制機制的變化,導致在管徑超過臨界管徑時顆粒最終毛細上升高度急劇下降.基于上述分析,高斯粒徑分布下,受振動的影響,粒徑較小的顆粒填充到大顆粒間的空隙中,使得管內顆粒柱的堆積密度高,阻礙了顆粒的流化,導致堵塞效應能夠在更大的管徑范圍內影響顆粒毛細上升行為,因而高斯粒徑分布下的臨界管徑更大.圖中結果還表明,當d≤ 8 mm 時,相同管徑條件下粒徑服從高斯分布的顆粒最終毛細上升高度較低,但平均毛細上升速度較大.前者是因為服從高斯粒徑分布的顆粒在管內的運動受堵塞效應的影響更大所致,后者則是由于堵塞效應使得顆粒毛細上升高度在更短的時間內達到穩定所致.當管徑d=9 mm 與d=10 mm 時,高斯粒徑分布下顆粒的平均毛細上升速度較小,這是由于顆粒毛細上升受控于不同的機制引起的.當d=11 mm 與d=12 mm 時,兩種粒徑分布下顆粒毛細上升動力學行為均受控于顆粒流化,使得顆粒平均毛細上升速度差異很小.豎直振動激勵下粒徑服從高斯分布的顆粒比單一粒徑顆粒能夠在更大的管徑下提升到更高的高度,對于顆粒物料的提升、輸運和采集具有重要意義,后續可以開展不同標準差下高斯分布顆粒以及其他粒徑分布下顆粒毛細上升動力學行為研究,為顆粒毛細效應動力學行為的應用提供依據.

圖6 不同粒徑分布下最終毛細上升高度和平均毛細上升速度隨管徑的變化Fig.6.Dependence of final capillary rising height and average capillary rising velocity on tube diameter under different particle size distribution.

4 結論

基于離散元方法對豎直振動激勵下顆粒毛細上升行為進行數值模擬研究,獲得了不同重力加速度條件下的顆粒毛細上升過程,揭示了顆粒最終毛細上升高度和平均毛細上升速度受重力加速度、水平振動分量和顆粒粒徑分布的影響規律.通過研究,得到以下結論:

1) 在低重力條件下,顆粒毛細效應仍能發生;隨著重力加速度的增加,最終毛細上升高度和平均毛細上升速度均呈現先增大后減小的變化趨勢,在重力加速度為6 m·s—2時達到極大值.

2) 最終毛細上升高度對水平振動分量不敏感,但平均毛細上升速度隨水平振動分量的增大而增加.

3) 粒徑服從高斯分布的顆粒與平均粒徑相等的單一粒徑顆粒相比,最終毛細上升高度最大值對應的臨界管徑更大,并且同處于堵塞效應影響區時平均毛細上升速度也更快.

猜你喜歡
振動
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
某調相機振動異常診斷分析與處理
大電機技術(2022年5期)2022-11-17 08:12:48
振動與頻率
天天愛科學(2020年6期)2020-09-10 07:22:44
This “Singing Highway”plays music
具非線性中立項的廣義Emden-Fowler微分方程的振動性
中立型Emden-Fowler微分方程的振動性
基于ANSYS的高速艇艉軸架軸系振動響應分析
船海工程(2015年4期)2016-01-05 15:53:26
主回路泵致聲振動分析
UF6振動激發態分子的振動-振動馳豫
計算物理(2014年2期)2014-03-11 17:01:44
帶有強迫項的高階差分方程解的振動性
主站蜘蛛池模板: 欧美成人一区午夜福利在线| 亚洲—日韩aV在线| 国产美女91呻吟求| 国产精品成人一区二区不卡| 中国一级特黄视频| 成人福利在线免费观看| 久久综合九九亚洲一区| lhav亚洲精品| 国产人成在线视频| 欧洲熟妇精品视频| 丁香五月婷婷激情基地| 伊人色综合久久天天| 日韩免费毛片| 狠狠色丁婷婷综合久久| 中文字幕第4页| 亚洲无线一二三四区男男| 亚洲无码37.| 波多野结衣的av一区二区三区| 亚洲成av人无码综合在线观看 | 美女高潮全身流白浆福利区| 老司机久久99久久精品播放| 亚洲一区二区黄色| 日韩小视频网站hq| 国产成人一级| 亚洲综合激情另类专区| 97久久超碰极品视觉盛宴| a级高清毛片| 美美女高清毛片视频免费观看| 国产主播在线观看| 国产小视频a在线观看| 视频一区视频二区中文精品| 直接黄91麻豆网站| 精品自窥自偷在线看| 中文字幕调教一区二区视频| 日韩最新中文字幕| 亚洲天堂区| 2021无码专区人妻系列日韩| 日本少妇又色又爽又高潮| 国产av剧情无码精品色午夜| 喷潮白浆直流在线播放| 亚洲欧美一区二区三区蜜芽| 国产精品专区第一页在线观看| 无码啪啪精品天堂浪潮av| 无码一区二区波多野结衣播放搜索| 亚洲av色吊丝无码| 国产无码网站在线观看| 欧美一级高清片久久99| 麻豆国产精品| 国产成人无码播放| 尤物精品视频一区二区三区| 99热国产这里只有精品无卡顿"| 亚洲日本中文综合在线| 超碰aⅴ人人做人人爽欧美| 亚洲AⅤ无码日韩AV无码网站| 欧美成人午夜影院| 99人妻碰碰碰久久久久禁片| 好吊色国产欧美日韩免费观看| 国产免费人成视频网| 2021最新国产精品网站| 国产成人午夜福利免费无码r| 精品国产香蕉在线播出| 亚洲高清无在码在线无弹窗| 好紧好深好大乳无码中文字幕| 欧美天堂久久| 国产免费久久精品99re不卡| 国国产a国产片免费麻豆| 97se亚洲综合在线韩国专区福利| 国产福利大秀91| 69精品在线观看| A级毛片高清免费视频就| 欧美不卡在线视频| 亚洲熟女偷拍| 精品国产自在在线在线观看| 毛片网站在线看| 伊人中文网| 亚洲国产成人综合精品2020| 亚洲综合婷婷激情| 日韩欧美高清视频| 呦女亚洲一区精品| 男女性色大片免费网站| 中字无码精油按摩中出视频| 在线播放91|