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

近自由面水下爆炸氣泡與結構相互作用數(shù)值計算研究

2016-01-15 05:24:16李健,潘力,林賢坤
振動與沖擊 2015年18期

第一作者李健男,博士,教授,1980年生

通信作者林賢坤男,博士,副教授,1976年生

近自由面水下爆炸氣泡與結構相互作用數(shù)值計算研究

李健1,潘力1,林賢坤1,榮吉利2,項大林2

(1. 廣西科技大學廣西汽車零部件與整車技術重點實驗室,廣西柳州545006; 2.北京理工大學宇航學院,北京100081)

摘要:近自由面水下爆炸氣泡在自由面與結構共同作用下會呈現(xiàn)出非常復雜的運動特性,利用體積加速度模型確定氣泡運動的初始半徑及速度,針對流場特點,通過自行開發(fā)的子程序實現(xiàn)了流場初始條件與邊界條件的定義,基于MSC.DYTRAN有限元軟件的流-固耦合處理方法對近水面水下爆炸氣泡與結構相互作用全過程進行數(shù)值模擬,通過與現(xiàn)有實驗數(shù)據(jù)的對比,驗證有限元模型建立和子程序開發(fā)的正確性。以此為基礎,研究氣泡與自由面、氣泡與結構的距離以及結構初始角度等參數(shù)對氣泡運動特性、射流角度、射流速度的影響規(guī)律,計算模型、方法及結果對相關的工程研究和計算具有一定參考價值。

關鍵詞:水下爆炸;氣泡脈動;近自由面;近剛性面;射流

基金項目:國家自然科學基金資助項目(51209042,11272057)

收稿日期:2014-07-04修改稿收到日期:2014-08-22

中圖分類號:O38文獻標志碼:A

Numerical study on interaction between bubble and structure near free surface in underwater explosion

LIJian1,PANLi1,LINXian-kun1,RONGJi-li2,XIANGDa-lin2(1. Guangxi key laboratory of Automobile Components and Vehicle technology, Guangxi University of Science and Technology, Liuzhou 545006, China2. School of Aerospace Engineering, Beijing Institute of Technology, Beijing 100081, China)

Abstract:Considering the complex phenomenon of the interaction between bubble and structure near free surface in underwater explosion, a volume-acceleration model was introduced to determine the initial condition of bubbles. Subroutines defining the initial value and boundary conditions of fluid field were developed by using MSC.DYTRAN finite element analysis software. The dynamic behaviors of bubbles near free surface and structure, including the jet of bubbles and dome of free surface were simulated and analyzed. The computational results of the behaviors of bubbles were compared with the experimental results.. The interactions between bubble and structure were studied systematically and summarized. The results show that the dynamic behaviors of bubbles have close relation to the distance parameter and initial angle of structure. The study is of referential value to relevant theoretical research and engineering calculation.

Key words:underwater explosion; bubble pulse; near free surface; near rigid wall; jet

氣泡在流體中運動及其與結構相互作用在物理、化學、醫(yī)學、工業(yè)生產(chǎn)以及兵器領域中廣泛存在并占據(jù)著重要的地位。氣泡在流體中運動及其與結構相互作用屬于典型的氣-液耦合與流-固耦合問題,在Bjerknes作用力的影響下,當氣泡在結構表面附近脈動時,膨脹階段會被結構表面輕微地排斥開,而在坍塌階段被結構表面強烈地吸引,這時氣泡坍塌出現(xiàn)在遠離結構一側,沿著坍塌方向會形成一股水射流,射流會以極高的速度穿透氣泡并最終作用在結構上,導致結構破壞。氣泡脈動在軍事上的應用非常廣泛。據(jù)美國、澳大利亞、英國、瑞典等國的24名專家對韓國“天安號”艦艇沉沒原因進行調查,結果表明導致艦艇沉沒的主要原因是魚雷在船體附近爆炸后形成的氣泡脈動以及氣泡坍塌后所產(chǎn)生的高速射流。正是由于氣泡在水中運動具有的這一特性,使其無論在民用還是軍事上都有極為廣泛的應用。所以,針對近自由面水下爆炸氣泡與結構的相互作用的研究對彈藥毀傷威力及水面艦艇目標易損性的評估都具有重要的意義。

目前國內外對近水面水下爆炸氣泡與結構相互作用研究方面主要以模型實驗或數(shù)值仿真分析為主。王詩平等[1-2]分別基于相似理論與邊界積分方法對水面艦艇模型在氣泡脈動作用下的特性進行實驗與數(shù)值仿真研究。初文華等[3]采用電極起泡方法進行了近自由面氣泡與結構相互作用的實驗研究。張弩等[4]采用雙漸進(DAA)方法研究了水下爆炸作用下水面艦船船體模型不同測點位移、速度、加速度特性。牟金磊等[5]利用MSC.DYTRAN軟件對固支方板水下爆炸氣泡作用下的動態(tài)響應進行仿真分析。劉云龍等[6]利用邊界元法進行了近自由面氣泡與傾斜壁面相互作用的數(shù)值計算。李健等[7]通過數(shù)值計算方法研究了近自由面水下爆炸氣泡與自由面水冢效應。Pearson等[8]采用邊界積分方法分別對單個及兩個氣泡與自由水面相互作用進行了研究,計算結果與實驗結果具有較好的一致性。Chahine等[9-10]通過實驗觀察到近水面爆炸氣泡與自由面相互作用的運動特性。Chew等[11]通過實驗測試方法研究了雙氣泡在剛性壁面附近的運動特性。Michael等[12]通過對大量實驗數(shù)據(jù)的總結分析,提出了近自由面水下爆炸水柱高度的經(jīng)驗公式。

從目前水下爆炸氣泡實驗與仿真分析研究成果看,國內外學者大多專注于采用邊界元法進行氣泡與結構相互作用的數(shù)值計算,邊界元方法僅在氣泡表面進行離散,單元數(shù)量較少,計算效率高,精度有保證,但在處理如射流擊穿氣泡等間斷問題時還存在一定的局限性。有限體積法(FVM)雖然計算量較大,但在處理間斷或多相流方面有一定優(yōu)勢,為此,本文采用基于FVM方法的MSC.DYTRAN軟件,以體積加速度模型為基礎,以水下爆炸實驗為依據(jù),結合多物質歐拉求解方法及流-固耦合技術,通過開發(fā)定義流場初始條件與邊界條件的子程序,對近自由面水下爆炸氣泡與結構相互作用過程進行仿真分析,研究氣泡與自由面、結構距離以及結構初始角度對氣泡運動規(guī)律以及射流特性的影響。

1理論分析

1.1氣泡運動初始條件

要模擬氣泡的運動,需要先求出氣泡的初始半徑。采用Hunter[13]提出的氣泡體積加速度模型,其主要是分別建立了藥包參數(shù)和流場壓力以及氣泡體積加速度與流場壓力之間的關系,并利用壓力等效關系最終確定氣泡的初始半徑與初始膨脹速度:

(1)

(2)

圖1為根據(jù)式(1)、式(2),55 gTNT炸藥在3.5 m水深起爆,時間t在0~0.45 ms時,氣泡半徑及徑向速度的時程曲線。炸藥引爆后,壓力瞬間增加,故初始時刻氣泡徑向速度變化率非常明顯,之后的沖擊波傳播過程中氣泡膨脹速度變化率逐漸平緩。針對這一規(guī)律,從圖1可知,與0~0.05 ms氣泡徑向速度梯度變化相比,在t>0.15 ms時,速度梯度變化較小,即可認為初始氣泡形成。取t=0.25 ms,對應的氣泡膨脹速度為200.1 m/s,初始半徑為0.1 m。

圖1 氣泡半徑及徑向速度時間歷程曲線 Fig.1 Time-history curves for radius and velocity of bubble

2有限元模型

為驗證本文所建立數(shù)值計算模型在模擬剛性壁面附近三維氣泡運動的有效性,根據(jù)文獻[14]中的實驗工況,即將55 g炸藥置于水下3.5 m處,并在其水平距離0.6 m處垂直放置一剛性平板,通過高速相機捕捉各時刻氣泡運動特性(實驗結果見圖2)。

根據(jù)實驗測試工況,同時考慮到水下爆炸氣泡運動特性,本文所建立的有限元模型為1.5 m× 1.5 m× 1.5 m的歐拉區(qū)域,整個計算區(qū)域全部用六面體歐拉單元離散,單元總數(shù)274999個,炸藥置于水下3.5 m處,水平方向距爆心位置0.6 m處垂直放置一剛性平板,采用流-固耦合方法定義流體與結構的相互作用,建立近剛性壁面三維氣泡運動模型。

根據(jù)流場初值及邊界條件特點,特別是流場存在靜水壓力的特點,利用MSC.DYTRAN軟件提供的定義初始條件和邊界條件的子程序接口EXINIT及EXFLOW2,借助FORTRAN語言開發(fā)了定義流場初值和邊界條件的子程序,最終可實現(xiàn)流場靜水壓力以及具有透射邊界條件的定義。其主要思路為:

(1)通過已劃分歐拉單元位置確定各節(jié)點坐標值

(2)根據(jù)計算工況,確定每個單元涉及節(jié)點所處位置的靜水壓力

(3)結合歐拉單元編號、節(jié)點編號、節(jié)點靜水壓值、狀態(tài)方程等,通過迭代算法確定最終單元的靜水壓力值。

(4)通過EXINIT接口將每個單元靜水壓值賦予相應編號的單元。

邊界條件的定義(其接口為EXFLOW2)與初始條件類似,只是僅需要對流場邊界的單元進行操作,且每個計算時間步長均需要重新確定并更新流場的邊界條件。

3結果與討論

3.1模型驗證

氣泡與壁面相互作用過程的計算結果見圖3。通過對比圖2與圖3發(fā)現(xiàn),氣泡在各時刻的運動均具有比較好的一致性,特別是在t=95 ms時,可以通過仿真結果清晰地看到氣泡坍塌會發(fā)生于遠離結構一側,形成由右下指向左上的射流,射流最終擊穿氣泡而作用在結構上。

圖2 55 g炸藥在水下3.5 m處爆炸氣泡與結構相互作用實驗測試結果 Fig.2 Evolutions of the bubble for 55 g detonation at a depth of 3.5 m from the experiment

圖3 55 g炸藥在水下3.5 m處爆炸氣泡與結構相互作用數(shù)值計算結果 Fig.3 Evolutions of the bubble for 55 g detonation at a depth of 3.5 m from the calculation

3.2特征參數(shù)對氣泡運動的影響

為研究剛性壁面附近氣泡的運動規(guī)律,本節(jié)將以近自由面水下爆炸氣泡與垂直放置剛性壁面相互作用數(shù)值計算模型為基礎,研究氣泡與剛性面距離、自由面距離、剛性壁面角度等參數(shù)對氣泡運動特性的影響規(guī)律。由于氣泡坍塌、射流方向、射流速度等物理現(xiàn)象的產(chǎn)生與爆心初始位置、氣泡脈動周期、氣泡最大半徑等關鍵參數(shù)有密切關系,因此以氣泡的最大半徑Rm為特征長度,中點處靜水壓力pamb為特征壓力,則特征時間為t=Rm(ρ/pamb)0.5,氣泡與剛性壁面水平距離為l=L/Rm,與自由面垂直距離為h=H/Rm,速度的無量綱表達式為v/(Rm/t),氣泡射流角度定義為射流方向與豎直向上(Y軸正向)的夾角,剛性面角度定義為其平面與水平面(X軸正向)的夾角,以此研究氣泡運動特性。

3.2.1與剛性壁面距離對氣泡運動的影響

圖4為氣泡與自由面距離保持在Rm,即h≡1.0,氣泡與剛性面距離變化時,氣泡與自由面的運動過程。從圖中可以看出,當氣泡與剛性面距離l=0.5時,由于此時自由面與剛性面Bjerknes效應均較為顯著,氣泡最終形成指向左下方的射流,與豎直方向夾角約為154°;而當特征距離增加至l=1.5,此時剛性面對氣泡運動影響減弱,表現(xiàn)為射流方向逐漸向豎直方向移動,當特征距離為l=2.5時,射流基本沿著豎直向下的方向。

圖4 氣泡與剛性面距離變化氣泡、自由面運動過程(h≡1.0) Fig.4 Evolutions of bubble contraction with different standoff distance that between bubble and structure (h≡1.0)

圖5 氣泡-剛性面距離與射流角度的關系(h≡1.0) Fig.5 Relationship betweenthe angle of jet and the standoff distance that between bubble and structure (h≡1.0)

圖5為氣泡與自由面距離h≡1.0時,氣泡與剛性面距離l變化與射流角度的關系曲線。從圖5可知,當h≡1.0時,射流角度會隨著l的增加而增加。其中當l=1.0時,射流角約為156°,當l增加至3.0~5.0時,射流角度約為178°~180°,表明此時剛性面對氣泡影響較小,氣泡坍塌于其頂端,形成背離自由面的向下射流。

圖6為氣泡與自由面距離h≡2.0時,氣泡與剛性面距離l變化與射流角度的關系曲線。從圖中可以看到,射流角度會隨著l的增加而增加。當l≥2.5 時,射流角度增加明顯,當l≥3.5時,射流方向基本垂直指向下方,表明此時剛性面對氣泡運動影響較小。

圖6 氣泡-剛性面距離與射流角度的關系(h≡2.0) Fig.6 Relationship between the angle of jet and the standoff distance that between bubble and structure (h≡2.0)

3.2.2與自由面距離對氣泡的影響

圖7為氣泡與剛性面距離L保持在Rm時,即l≡1.0時,氣泡與自由面距離h對射流角度的影響曲線。從圖中可以看到,隨著h的增加,氣泡射流角度將逐漸減小。其中當h=1.0時,此時由于自由面與垂直剛性壁面的Bjerknes力影響均較為顯著,射流角度較大,與豎直方向夾角約為156°,當h在3.0~5.0范圍內時,自由面對氣泡運動影響明顯較弱,剛性面Bjerknes作用力占優(yōu),射流角度為85°~91°,基本是垂直射向剛性面,此時,自由面效應對氣泡運動的影響較小,可忽略不計。

圖8為氣泡與剛性面距離L保持在2倍Rm時,即l≡2.0時,氣泡與自由面距離h對射流角度的影響曲線。從圖中可以看到,與前述類似,射流角度會隨著h的增加而減小。對比圖7、圖8可知,當h=1.0時,自由面Bjerknes力的影響均非常顯著,此時,隨著l的增加,剛性面對氣泡運動影響力逐漸減弱,射流角度隨之增大;而當氣泡距自由面較遠時(h≥3.0),自由面對氣泡影響較小,剛性面Bjerknes力占主導因素,表現(xiàn)為隨著l的增加,射流角度減小。

圖7 氣泡-自由面距離與射流角度的關系(l≡1.0) Fig.7 Relationship betweenthe angle of jet and the standoff distance that between bubble and free surface (l≡1.0)

圖8 氣泡-自由面距離與射流角度的關系(l≡2.0) Fig.8 Relationship betweenthe angle of jet and the standoff distance that between bubble and free surface(l≡2.0)

3.2.2剛性壁面夾角對氣泡運動的影響

由于特征距離h與l均為1.0時對氣泡運動影響較大,本節(jié)研究l=1.0且h=1.0時,剛性面角度變化對氣泡射流、射流速度的影響。圖9為l=1.0,h=1.0,剛性面角度分別為30°與150°時氣泡與自由面的運動過程。從圖中可以看到,當剛性面角度為30°時,自由面水冢將沿著結構平面方向運動,氣泡首先沿著剛性面方向坍塌,最終形成指向剛性面的射流,射流與豎直方向夾角約為40°;當剛性面與水平面夾角為150°時,自由面在氣泡脈動作用下形成沿豎直向上的高而細的水柱,氣泡在其頂端發(fā)生坍塌,最終形成指向剛性面的射流,其與豎直方向夾角約為170°,基本沿著豎直向下方向。

無量綱時間: (1)0.578, (2)0.866, (3)1.190, (4)1.443, (5)1.587, (6)1.659, (7)1.840, (8) 2.020 圖9 l=1.0, h=1.0時剛性壁面角度變化時 氣泡坍塌與射流形成過程 Fig.9 Evolutions of bubble contraction with different angle of structure (l=1.0, h=1.0)

圖10、圖11分別為剛性面角度變化與氣泡射流角度、無量綱速度的變化曲線。從圖中可知,隨著剛性面角度的增加,氣泡射流角度也會隨之增加,特別是剛性面角度在90°(垂直放置)以內時,氣泡射流角度增加明顯,其主要原因是此工況下,自由面與剛性面對氣泡運動、射流影響效果不同,即自由面導致氣泡產(chǎn)向豎直向的下坍塌與射流,而結構導致氣泡產(chǎn)生沿著兩者連線方向的射流,當剛性面角度在90°內的變化時,自由面與剛性面對氣泡沿豎直方向作用效果相反,因此,當剛性面角度在90°內變化時,射流角度變化較為明顯,而當剛性面與水平夾角大于90°時,由于剛性面與自由面對氣泡影響效果基本相同,均導致氣泡產(chǎn)生指向下的射流,因此,此工況下氣泡角度變化較小,速度增加相對較快。

圖10 l=1.0, h=1.0時射流角度與剛性面角度的關系 Fig.10 Relationship between the angle of jet and the angle of structure (l=1.0, h=1.0)

圖11 l=1.0, h=1.0時射流速度與剛性面角度的關系 Fig.11 Relationship between the velocity of jet and the angle of structure (l=1.0, h=1.0)

4結論

采用體積加速度模型確定氣泡的初始條件,對氣泡與結構相互作用全過程進行了仿真分析,通過與實驗結果的對比,驗證了氣泡初始條件確定方法、子程序開發(fā)與模型建立的正確性,以此模型為基礎,通過仿真計算,研究了近自由面水下爆炸氣泡與結構相互作用過程中氣泡與自由面、氣泡與剛性面以及剛性面角度對氣泡運動特性的影響。分析總結相關規(guī)律,獲得了如下一些主要結論:

(1)在特征距離較近時,無論是自由面還是結構均會對氣泡運動產(chǎn)生影響,一般情況下,在結構附近脈動的氣泡會形成指向結構的射流,而自由面附近脈動的氣泡會形成豎直向下射流,隨著特征距離的變化,二者對氣泡運動、射流方向、射流速度均有不同程度的影響。

(2)當氣泡與剛性面距離在3倍Rm以內時,剛性面的Bjerknes效應較明顯,氣泡射流角度隨著氣泡與剛性面距離的增加而增加,當氣泡與剛性面距離超過3.5倍Rm時,剛性面對氣泡運動影響明顯減弱,基本可以忽略不計。

(3)當氣泡與自由面距離在3倍Rm以內時,自由面的Bjerknes效應較明顯,氣泡射流角度隨著氣泡與自由面距離的增加而減小,當氣泡與自由面距離超過3倍Rm時,自由面對氣泡運動影響明顯減弱,基本可忽略不計。

(4)當氣泡與自由面、剛性面距離均較近(如Rm)時,剛性面角度在90°以內時,隨著剛性面與水平面夾角的增加,氣泡射流角度增加明顯,射流速度增加較慢,當剛性面角度大于90°時,氣泡射流角度增加較慢,射流速度增加明顯。

參考文獻

[1]王詩平, 朱楓, 吳超, 等. 爆炸氣泡與艦船相互作用的相似性方法研究[J]. 中國造船, 2012, 53(2):30-38.

WANG Shi-ping, ZHU Feng, WU Chao, et al. Similarity method in scaled experiment of interaction between underwater explosion bubbles and ships[J]. Shipbuilding of China, 2012, 53(2):30-38.

[2]王詩平, 孫士麗, 張阿漫, 等. 沖擊波和氣泡作用下艦船結構動態(tài)響應的數(shù)值模擬[J]. 爆炸與沖擊, 2011, 31(4): 367-372.

WANG Shi-ping, SUN Shi-li, ZHANG A-man, et al. Numerical simulation of dynamic response of warship structures subjected to underwater explosion shockwaves and bubbles [J]. Explosoin and Shock Waves, 2011, 31(4): 367-372.

[3]初文華, 張阿漫, 王詩平. 壁面與自由液面聯(lián)合作用下氣泡動態(tài)特性實驗研究[J]振動與沖擊,2013, 32(13):112-117.

CHU Wen-hua, ZHANG A-man, WANG Shi-ping. Experimental study on bubble pulse features under combined action of wall and free surface[J]. Journal of Vibration and Shock, 2013, 32(13):112-117.

[4]張弩, 宗智, 張文鵬, 等.基于雙漸進方法的水下爆炸氣泡載荷作用下艦船的動態(tài)響應分析[J]. 振動與沖擊, 2012, 31(23): 50-56.

ZHANG Nu,ZONG Zhi,ZHANG Wen-peng, et al. Dynamic response of a ship hull stucture subjected to an underwater explosion bubble based on doubly asymptotic approximation method[J]. Journal of Vibration and Shock, 2012, 31(23): 50-56.

[5]牟金磊, 李海濤. 固支方板在水下爆炸氣泡射流作用下動態(tài)響應[J]. 噪聲與振動控制, 2012, 6: 125-129.

MU Jin-lei, LI Hai-tao.Dynamic response of stiffened square plate subjected to jets load due to underwater explosion bubbles[J]. Noise and Vibration Control, 2012, 6: 125-129.

[6]劉云龍, 汪玉, 張阿漫. 有傾角的豎直壁面附近氣泡與自由面相互作用研究[J]. 物理學報,2013, 62(21):267-276.

LIU Yun-long, WANG Yu, ZHANG A-man. Interaction between bubble and free surface near vertical wall with inclination[J]. Acta Physica Sinica, 2013, 62(21):267-276.

[7]Li J, Rong J L. Bubble and free surface dynamics in shallow underwater explosion[J]. Ocean Engineering, 2011, 38: 1861-1868.

[8]PearsonA, Cox E, Blake J R, et al. Bubble interactions near a free surface[J]. Engineering Analysis with Boundary Elements, 2004, 28(4):295-313.

[9]Chahine G L. Interaction between an oscillating bubble and a free surface[J]. ASME, Journal of Fluids Engineering, 1977, 99(4): 709-716.

[10]Gibson D C. Cavitation adjacent to place boundaries[R]. Proc 3rd Conf Hydraulic Fluid Mech. Sydney, P.R. Austr Inst Eng. 1968, 210-214.

[11]Chew L W, Klaseboer E, Ohl S W. Interaction of two oscillating bubbles near a rigid boundary[J]. Experimental Thermal and Fluid Science, 2013, 44: 108-113.

[12]Michael M S. Explosion effects and properties: Part2 explosion effects in water[R]. ADA-056694 Maryland, USA: Naval Surface Weapons Center, 1978.

[13]Hunter K S. Global-shape-function models of an underwater explosion bubble[D]. Faculty of Graduate School of the University of Colorado, Department of Mechanical Engineering, 2001.

[14]Klaseboer E, Khoo B C, Hung K C. Dynamics of an oscillating bubble near a floating structure[J]. Journal of Fluids and Structures, 2005, 21:395-412.

主站蜘蛛池模板: 99免费视频观看| 国产成人亚洲精品无码电影| 直接黄91麻豆网站| 中国国产一级毛片| 日本午夜影院| 久久这里只有精品8| 伊人久久大香线蕉影院| 国产又大又粗又猛又爽的视频| 区国产精品搜索视频| 国产区精品高清在线观看| 国产激情无码一区二区三区免费| 欧美一级黄片一区2区| 国产乱人伦AV在线A| 国产肉感大码AV无码| 成人一区在线| 久久久精品久久久久三级| 欧美中文字幕一区| 麻豆国产在线不卡一区二区| 日本www色视频| 国产精品成人AⅤ在线一二三四| 日本91视频| 欧美日本在线一区二区三区| 欧美午夜久久| 亚洲国模精品一区| 曰韩免费无码AV一区二区| 超碰aⅴ人人做人人爽欧美| 中国精品久久| 成人福利视频网| 欧美yw精品日本国产精品| a天堂视频| 91探花在线观看国产最新| 午夜激情婷婷| 国产精品久久自在自2021| 亚洲成人在线网| 欧美日韩国产成人高清视频| 国产手机在线小视频免费观看 | 99精品国产高清一区二区| 青草精品视频| 55夜色66夜色国产精品视频| 亚洲人人视频| 日日碰狠狠添天天爽| 国产在线自乱拍播放| 99re热精品视频国产免费| 国产一区在线视频观看| 999精品色在线观看| 青青草91视频| 婷婷午夜天| 国产亚洲精品自在线| 99无码中文字幕视频| 天天视频在线91频| h视频在线观看网站| 国产麻豆aⅴ精品无码| 日韩国产亚洲一区二区在线观看| 午夜a视频| 爱爱影院18禁免费| 一级毛片在线播放| 51国产偷自视频区视频手机观看| 美女裸体18禁网站| 国产拍在线| 日韩福利在线观看| 国产在线视频自拍| 亚洲码在线中文在线观看| 日韩高清成人| 精品乱码久久久久久久| 久久久精品久久久久三级| 亚洲成a人在线观看| 国产一二视频| 免费日韩在线视频| 毛片视频网| 精品成人一区二区三区电影| 精品国产Av电影无码久久久| 久久综合丝袜日本网| 不卡午夜视频| 亚洲精品欧美日韩在线| 欧美国产日韩另类| 五月激情综合网| 日韩精品成人网页视频在线 | 91久久天天躁狠狠躁夜夜| 九九九精品成人免费视频7| 国产精品视频公开费视频| 亚洲成人在线网| 精品综合久久久久久97|