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

氣泡及其破碎興波對浮動沖擊平臺影響探究*

2019-10-17 07:35:06王志凱姚熊亮楊娜娜
爆炸與沖擊 2019年9期

王志凱,周 鵬,孫 波,姚熊亮,楊娜娜

(1. 哈爾濱工程大學船舶工程學院,黑龍江 哈爾濱 150001;2. 哈爾濱工程大學機電工程學院,黑龍江 哈爾濱 150001)

艦艇生命力及其作戰效能是作戰艦艇能夠發揮其攻防效能的關鍵,而艦載設備系統抗爆抗沖擊能力作為其中十分重要的技術指標,直接影響艦艇海上作戰能力,各海軍強國對此均非常重視,并投入大量資源進行研究[1]。但實船實驗周期長、耗費成本高、實施難度大、危險性高、且具有不可重復性。故目前世界各國主要是運用不同實驗裝備代替實船進行考察為主,數值模擬計算為輔進行考察[2]。小型艦載設備的實驗裝備主要是陸上沖擊機,大型艦載設備主要采用浮動沖擊平臺(floating shock platform)作為其標準實驗裝備[3](以下簡稱平臺或FSP)。

實驗時,炸藥產生的沖擊波傳播到水面和海底時會形成反射波即水面截斷效應,同時,自由面的空泡和氣泡脈動,氣泡、自由面、平臺相互耦合,過程十分復雜。氣泡在經過脈動后上浮到自由面處會發生破碎,在破碎的自由面附近產生明顯的初始擾動,產生較為特殊的興波現象;在自由面處產生十幾米甚至幾十米的水冢,水冢形成時可能會直接沖擊平臺,改變平臺狀態。當水冢上升到最高位置時,由于重力作用,會產生自由落體運動,引起自由面的擾動,在落水中心處產生向外部擴散的波浪,對周圍的平臺產生力的作用,在一定程度上改變平臺運動狀態。而氣泡脈動過程中,因為氣泡載荷的低頻性,水中物體不僅會被水下爆炸氣泡載荷造成嚴重損傷,還會出現比較大的剛體位移,而在以往的理論和數值研究中,這種現象都被忽略[4-6]。同時,Zhang 等[7]發現當長寬比較小的艦船結構承受水下爆炸載荷時,艦船剛體運動占整個運動比例較大,不可忽略。當氣泡上浮到自由面處時發生破裂,在破裂的自由面附近產生明顯的初始擾動,進一步形成較為特殊的興波現象。這一現象已經探討得比較深入[8-9],特別是Wang[10]的工作中,運用邊界元方法發現氣泡破碎后產生的興波包含很大一部分能量。而沖擊環境則是單自由度系統受到外界環境激勵時,質量的相對位移或速度或加速度最大的幅值隨頻率變化的關系。對于固定系統,當不同激勵作用時,造成系統的沖擊環境不同,而整個水下爆炸過程中,平臺承受的外界激勵種類較多(沖擊波、氣泡脈動與水冢推動、后期氣泡興波等),而在以往沖擊環境計算中幾乎都是只考慮沖擊波作用影響[11-12],而忽略氣泡脈動、水冢、涌流、后期波浪的影響。

本文中利用LS-DYNA 軟件模擬整個氣泡、自由面、平臺的耦合過程,考慮氣泡脈動及其破碎興波、水冢、涌流等作用,探討水下爆炸過程中平臺運動變化與沖擊環境變化情況。

1 近水面近壁面混合邊界氣泡動力學特性分析及模型驗證

1.1 有限元模型

實際沖擊平臺考察設備時,當考察點距離平臺較近時,氣泡不僅受到自由面的強烈作用還會受到平臺的作用,實驗數據來自文獻[13],4 g PETN 混合炸藥(可等效為5.2 g TNT 炸藥)在2 m 長的鋼制正方體水箱中進行水下爆炸實驗,爆心位于水下0.26 m,距離剛性壁面0.19 m。有限元模型中炸藥、水、空氣采用多物質ALE 算法,空氣域與水域均為邊長2 m 的立方體,采用單元類型為歐拉六面體的漸變網格,藥包附近網格尺寸最小,為0.9 mm。流域邊緣網格尺寸最大,為0.6 m。流域網格總數為160 萬;TNT 為當量5.2 g 的球形藥包,布放位置與對應實驗工況相同;歐拉域四周與底面設置剛性壁面,頂部添加大氣壓力。工況示意圖如圖1 所示。

圖 1 混合邊界條件時工況示意圖及有限元圖Fig. 1 Illustration of working and finite element view in complex boundary conditions

1.1.1 空氣

空氣的狀態方程為理想氣體狀態方程:

式中:pg是空氣壓力,ρ 是空氣密度,ρ=1.225 kg/m3,γ 是絕熱指數,γ=1.4,E 是空氣內能,E=2.11×105J/kg。

1.1.2 水

水的狀態方程采用Grüneusen 狀態方程,在LS-DYNA 中通過*MAT_NALL 材料模型實現:

式中:pw為水的壓力;μ為水的壓縮系數,μ=ρ/ρ0-1;水的初始密度ρ0=1 000 kg/m3;c 為水中聲速,c=1 500 m/s;S1、S2、S3為曲線擬合參數,S1=2.56,S2=1.986,S3=1.226 8;γ0為Grüneusen 參數,γ0=0.5;α 為μ和γ0的一階體積修正量,α=0;E 為水的體積內能,單位為J/m3。

1.1.3 TNT 炸藥

通過JWL 狀態方程對TNT 炸藥的爆炸狀態進行描述,其表現形式為:

式中: pT為爆轟產物壓力,η 為爆炸產物與初始炸藥密度ρ0的比值,ρ0=1 630 kg/m3;A、B、R1、R2、ω 為狀態方程參數,A=371.2 GPa、B=3.231 GPa、R1=4.15、R2=0.95、ω=0.3;E 為爆轟產物內能,E=8.0 GJ/kg。

1.2 計算結果與實驗結果對比

圖2 為4 g PETN 混合炸藥在位于水下0.26 m、距離剛性壁面0.19 m 爆炸時,氣泡的實驗結果[13]與計算結果,剛性壁面位于氣泡的右側。

圖 2 混合邊界條件時爆炸氣泡實驗結果與計算結果對比Fig. 2 Experimental and numerical results of bubble in complex boundary conditions

通過比較實驗結果和數值結果,發現兩者在氣泡形態以及產生相應形態所對應的時刻具有高度的一致性,由此可以認為本次使用的ALE 方法能夠較好地模擬近水面近壁面混合邊界條件下,爆炸氣泡的膨脹、收縮、射流、潰滅等復雜現象,可以進行下一步更深層次的工作。

2 計算模型及工況設置

2.1 平臺有限元模型

按照美國中型平臺尺寸建立有限元模型[14],長12.2 m、寬6.1 m、總高度4.2 m、雙層底高1 m,雙層底之間通過縱桁和實肋板加強,空載吃水1.2 m,空載重心高1.01 m。在保證主尺寸、重心、吃水量等影響搖蕩性能的參數不變前提下對模型進行部分簡化。王軍等[11]發現剛體運動能量約占整個外來能量的60%~70%以上,平臺工作時,剛體運動占絕大部分。同時由于氣泡脈動以及后期水冢落水產生波浪的時間較長,為縮短計算時間,本文中將平臺材料簡化為剛體進行探究。具體結構示意圖如圖3 所示。

圖 3 平臺有限元模型Fig. 3 Finite element model of floating shock platform

2.2 工況設置

平臺有限元模型引用美軍標的設計標準,因此實驗時同樣采用美軍標相對應的工況進行計算。為提高安全標準,將原標準中40.8 kg HBX-1 炸藥替換成60 kg TNT 球形炸藥進行考核。并參照美軍標的設計要求,選擇距離自由液面垂向距離為8.23 m,沿平臺寬度方向設置多個水平爆點,具體示意圖如圖4(a)~(b)所示,工況參數如表1 所示。本文中基于平面波假定,從結構遮擋的沖擊波能量相等角度定義無量綱沖擊因子C,具體計算如下:

式中: W 為炸藥當量,kg; R 為總爆距,m。

3 數值計算結果分析

3.1 平臺運動形態變化

圖5 為工況2 的自由面和平臺運動圖,此工況下爆距較近。可以看出,在圖5(a)到圖5(b)過程中,氣泡一直膨脹,并推動平臺發生順時針橫搖運動,最大橫搖角約3°。在t=0.340 s (圖5(b))時刻氣泡體積膨脹到最大值,同時由于邊界條件的不連續性,平臺附近流場出現高壓區,在平臺兩側出現自由面突起。同時氣泡內壓遠小于周圍流體,周圍流體受到輻射的負壓,因為平臺的存在,此負壓從氣泡一直延伸到平臺的下表面,故在氣泡收縮過程中,自由面突起與平臺交界點之間的自由面迅速向內凹陷,產生空腔,使得平臺左側失去浮力的支撐作用,從而產生逆時針的恢復力矩,平臺開始向左下方運動,橫搖角度減小。直到t=0.640 s (圖5(c))時刻氣泡體積達到最小值,氣泡第一周期結束,進入第二周期并再次開始膨脹,并出現明顯的上浮運動,平臺再次向右上方移動。在t=1.010 s (圖5(d))時氣泡體積膨脹到最大,由于氣泡與平臺之間距離過近,平臺阻止了氣泡的完全膨脹,同時氣泡也推動平臺,使平臺開始發生大角度的順時針橫搖運動,此時平臺橫搖角達到12.3°,下一時刻平臺開始隨著氣泡的收縮而做逆時針運動,自由面突起由于慣性繼續上升。在t=1.240 s (圖5(e))時,氣泡縮小到最小,開始下一周期的膨脹,平臺橫搖角繼續增加。在t=1.520 s (圖5(f))時,氣泡體積第三次膨脹到最大,并已經開始產生皇冠圍裙,此后水冢開始向周圍及上方擴散,由于滯后流以及上升水冢的作用,平臺一直發生大角度的橫搖運動,直到t=2.320 s(圖5(h))達到最大值,此時橫搖角度已達到40°,嚴重威脅平臺的穩性,此時氣泡已經破碎,后續討論將氣泡移除。由于下落的水冢直接抨擊平臺左側底部,對平臺重心產生逆時針的力矩作用,平臺發生逆時針運動,橫搖角開始減小,同時部分水冢落在自由面上,自由面產生孤波向右傳播。推動平臺一直向右平移移動。在t=3.620 s (圖5(j))時刻,平臺重心回復到平衡位置,同時由于慣性作用及自由面波浪的作用,平臺繼續進行逆時針運動。在t=4.350 s (圖5(k))時刻,橫搖角達到最大(約7°),此時波浪已經很小,此后平臺由于波浪的作用和慣性開始發生向右平移和小角度的往復橫搖運動。

圖 5 工況2 水冢形態演變及平臺運動響應Fig. 5 Water spike's shape and evolution process and platform's motion response in working condition 2

圖6 為工況6 的自由面和平臺運動圖,此工況下爆距較遠。圖6(a)~(f)時,在整個上浮氣泡脈動過程中,平臺受到水冢的作用力很小,運動響應較小。平臺未對氣泡運動特性及水冢演變特征產生影響。在圖6(g)時氣泡開始破碎產生興波,中心水柱在慣性作用下繼續上升,皇冠圍裙則繼續向四周飛濺,并推動平臺產生順時針的橫搖運動。在圖6(h)時可以看到部分水冢由于重力已經開始下落并抨擊水面,產生向兩側擴散的孤立波,平臺由于波浪作用繼續進行順時針橫搖運動。在圖6(i)時,整個皇冠圍裙已經全部落入水中,中心水柱也已開始下落,平臺橫搖角達到最大值(約5°),此時平臺由于恢復力矩開始做逆時針橫搖運動,并且平臺上方水冢落下拍擊平臺上表面。在圖6(j)時平臺逆時針橫搖達到最大值,由于始終有孤立波的作用,此時橫搖角為-0.35°,由于波浪及恢復力矩,平臺繼續做順時針橫搖運動,并且此時發現水冢已經完全融入水面,平臺上表面水花開始向上飛落。在圖6(k)時發現,水面由于水冢的抨擊作用,產生大量水花,平臺上表面水花飛濺得更加劇烈,此時,平臺橫搖角達到最大值(約5°)。由于波浪和回復力矩的作用,平臺開始不斷進行小角度有阻尼的升沉運動和橫搖運動。

圖 6 工況6 水冢形態演變及平臺運動響應Fig. 6 Water spike's shape and evolution process and platform's motion response in working condition 6

3.2 平臺三自由度剛性位移

下面給出不同水平距離下平臺的三自由度曲線,升沉自由度曲線如圖7 所示。

由圖7 可以發現,在上浮氣泡脈動階段,由于氣泡脈動載荷基本是以水中聲速(約1 500 m/s)向外傳播,平臺上升位移最大與最小值時刻嚴格對應氣泡體積最大與最小值時刻,無明顯的延遲現象。同時發現升沉運動位移隨著水平位移的增大而減小。這是由于隨著水平位移變遠,總的氣泡作用力減小,同時用于抬升平臺的分作用力也減小。后期當氣泡破碎以及自由面水柱跌落階段,可以發現仍然是升沉運動位移隨著水平位移的增大而減小,而對應最大值的時刻越早。這是因為平臺上升到最大值后,由于重力作用,勢能轉化為動能再轉化為勢能,上升階段距離近的勢能大,所以后期轉化的勢能大。距離變遠時能量總量較小,能量相互轉化完成所需時間較短,使得對應最大值的時刻較早。同時也說明在爆炸后期氣泡破碎以及自由面水柱跌落形成的興波,對平臺升沉作用影響不大。

圖 7 不同水平位移下平臺升沉運動響應Fig. 7 Platform's heave motion response in different horizontal displacements

由圖8 中可以看出,在上浮氣泡脈動階段,在C=1.04 和C=0.94 兩個工況下第二次氣泡膨脹時平臺的橫搖角特別大,達到40°,已經喪失安全性和使用性。由圖5 對應時刻可以看出,此時氣泡在自由面附近直接膨脹,產生的皇冠水冢圍裙直接作用到平臺左側底部,推動平臺一直發生大角度的橫搖運動。而距離更近的工況由于力臂太小,使橫搖力矩很小,基本只造成升沉運動。距離太遠則由于力太小,使得橫搖力矩很小,造成小角度的橫搖運動。同時發現在氣泡破碎以及自由面水柱跌落階段,距離較近時對平臺的影響僅為由于慣性引起的橫搖角相對較小的來回擺動。而距離變大時,橫搖角是先增大后減小。原因為當距離特別近時,本應該用于興波的大部分能量在氣泡脈動階段直接造成平臺抬升和橫搖,轉化為了平臺的動能,從而減弱了后期波浪的能量。由圖5 的對應時刻可以看出,距離較近時,水冢落水時有很大一部分直接砸落在平臺上表面上,造成的更多的是升沉運動,而抨擊水面產生的興波和氣泡破碎產生的興波也由于距離過近,并未發展完全;在最后階段,平臺抬升使得干舷明顯增加,此時波浪的抨擊作用在平臺左側的面積較大,波浪的能量更多地轉化為橫傾運動,而不是橫搖。而當距離較大時,橫搖角是先增大后減小,而且可能第二個橫搖角最大值略大于第一個橫搖角最大值。這是因為波浪作用具有一個從無到有的發展區間,而且水冢落水主要分成兩部份,一部份是皇冠圍裙落水較早且距離較近,另一部份是中心水柱落水,時間較晚距離相對較遠,兩者在自由面上產生的波浪有一個疊加作用。

圖 8 不同水平位移下平臺橫搖運動響應Fig. 8 Platform's roll motion response in different horizontal displacements

由圖9 可以發現,在前期的氣泡上浮階段,各工況下引起的橫蕩位移值相差很小,這說明橫蕩位移主要是由后期氣泡破碎興波及水冢抨擊水所產生的波浪對平臺進行沖擊引起的。工況1 后期橫蕩發生很小主要是因為:由于爆距較近,此時氣泡能量幾乎完全使平臺進行升沉運動,對平臺橫蕩運動的貢獻不大。而剩余工況隨著水平距離增大,橫蕩位移變化較慢。這是因為波浪的產生及其對浮動的作用不是瞬間完成的,需要一定的傳播時間。距離較近,波浪作用的時間較早,故橫蕩位移在相同的時間內變化較為明顯。繼續計算使橫搖第二次達到最大值,此刻波浪作用極為有限,并且由中心水柱產生的波浪也已經作用到平臺上或仍在傳播,且由于能量損失造成波高已經很小,僅對平臺造成橫蕩運動,對安全性和使用性不會有影響。

圖 9 不同水平位移下平臺橫蕩運動響應探究Fig. 9 Platform's traverse motion response in different horizontal displacements

3.3 平臺沖擊響應變化

對于固定系統,當不同激勵作用時,造成系統的沖擊環境不同,而整個平臺考察設備過程中,平臺承受的外界激勵種類較多,為此本節分別進行沖擊波僅峰值作用(作用時間0.05 s)、沖擊波完全作用(作用時間0.5 s)、沖擊波氣泡脈動作用(作用時間2 s)、沖擊波氣泡脈動與后期波浪共同作用(作用時間7.5 s)總計4 種情況的計算,以探討氣泡及其破碎興波對平臺沖擊影響。工況2 在4 種情況下的沖擊加速度時歷曲線如圖10 所示。并分別畫出各激勵載荷作用下,較近爆距工況2 和較遠爆距工況6 沖擊響應譜曲線對比。

從圖11 可以發現,當僅沖擊波峰值對平臺作用時,此時處于沖擊波初始響應階段,由于外界激勵比較簡單,作用時間較短(此時作用時間0.05 s,只有大于20 Hz 的頻率值才達到完整一個周期,5~20 Hz 時達到1/4 周期以上,更低頻不可能達到峰值),低頻段此時出現明顯的等位移譜、等速度譜、等加速度譜。當沖擊波完全作用時,對比發現,此時等速度譜與等加速度譜幾乎不變,等位移譜值在較低頻時增加,這是因為沖擊波完全作用時,時間較長,達到0.5 s,此時2 Hz 以內的頻率值都已達到完整一個周期,0.5~2 Hz 值也已達到1/4 周期以上,故此時外界激勵基本涵蓋峰值段,等位移譜值在較低頻時增加。

從圖12 可知,當有無氣泡脈動作用時,工況2 與工況6 兩者的高頻部分幾乎不變,即譜加速度值不變,而在低頻中頻部分則出現較多共振點,且爆距較近的工況2 數值影響程度更大,這說明氣泡脈動作用與水冢直接沖擊作用對沖擊譜值影響主要是中低頻性,影響譜速度值、譜位移值,對譜加速度值沒有影響。從圖13 可以發現,在考慮波浪作用的影響時,工況2 與工況6 兩者各自有無波浪的原始譜基本完全相同,僅在8 Hz 處有明顯突起,也就是有一個周期0.125 s 左右的外載荷激勵。圖14 為波浪作用下平臺加速度的時歷曲線。觀察圖14 發現波浪周期約為0.12 s,此激勵正是后期波浪造成的,即水冢抨擊水面和氣泡興波對平臺形成周期性的激勵,而且只是通過在此激勵對應頻率從而引起共振來改變沖擊環境。

圖 10 工況2 沖擊加速度時歷曲線Fig. 10 Shock acceleration time history curve of working condition 2

圖 11 沖擊波作用時間對沖擊響應譜影響對比圖Fig. 11 Effect of shock wave's action time on shock spectrum in comparison

圖15 為其他工況有無波浪的沖擊譜計算結果,當沖擊因子小于1 時,譜速度和譜加速度的變化較緩,幾乎是一條直線,而大于1 之后則急速上升,可發現是因為距離過近時,水冢直接推動平臺引起的。可發現考慮氣泡后期波浪對整個計算結果的結果極為有限,故在后期計算不同沖擊因子對平臺沖擊的響應結果時,為了節省計算時間和計算成本,只計算到氣泡脈動載荷截止即可。將忽略氣泡破碎興波和波浪影響的平臺在不同沖擊因子C 下譜速度和譜加速度曲線擬合成公式,如下表2 所示,方便以后的直接應用。

圖 12 氣泡脈動對沖擊響應譜影響對比圖Fig. 12 Effect of bubble pulsation on shock spectrum in comparison

圖 13 后期波浪對沖擊響應譜影響對比圖Fig. 13 Effect of later wave on shock spectrum in comparison

圖 14 波浪作用對平臺加速度響應影響圖Fig. 14 Later wave's effect on shock acceleration response time

圖 15 波浪對比時沖擊譜參數隨沖擊因子變化關系Fig. 15 Relationship between shock spectrum parameters and shock factor in wave contrast

表 2 沖擊環境參數與沖擊因子關系式Table 2 Relationship between shock response spectrum parameters and shock factor

4 結 論

本文利用LS-DYNA 軟件模擬整個氣泡、自由面、浮動沖擊平臺的耦合過程,考慮氣泡脈動及其破碎興波、水冢、涌流等作用,探討水下爆炸過程中浮動沖擊平臺運動變化與沖擊環境變化情況,得到以下結論:(1)在整個水下爆炸過程中,氣泡、自由面、浮動沖擊平臺會發生強烈的耦合作用,爆炸氣泡破碎前后遵循不同規律:氣泡破碎前,浮動沖擊平臺的剛性位移主要與氣泡脈動力、水冢直接推動有關,主要引起浮動沖擊平臺的橫搖運動和升沉運動。氣泡破碎后,浮動沖擊平臺的剛性位移由氣泡破碎興波和水冢拍擊水面產生的孤立波引起,主要造成浮動沖擊平臺的橫搖運動和橫蕩運動,影響浮動沖擊平臺的安全性和使用性。(2)在水下爆炸過程中,沖擊波是影響艦船沖擊環境的主要因素,由于氣泡的低頻性,氣泡脈動及水冢對浮動沖擊平臺直接沖擊作用小幅度增加譜速度值、譜位移值,對譜加速度值幾乎無影響。(3)水冢抨擊水面和氣泡興波形成的波浪,對浮動沖擊平臺造成的激勵載荷呈周期性,而且只是通過在此激勵對應頻率引起共振來改變沖擊環境,且由于波浪載荷很小,對沖擊環境影響較小,在計算時可以忽略不計。

主站蜘蛛池模板: 亚洲V日韩V无码一区二区| 亚洲伊人久久精品影院| 亚洲欧美另类中文字幕| www.91中文字幕| 无码专区国产精品第一页| 综合社区亚洲熟妇p| 亚洲人成人无码www| 亚洲无码视频喷水| 中文字幕第4页| 亚洲性视频网站| 亚洲色无码专线精品观看| 国内老司机精品视频在线播出| 波多野结衣国产精品| 波多野结衣一二三| 亚洲人成亚洲精品| 欧美.成人.综合在线| 亚洲国产av无码综合原创国产| 欧洲亚洲欧美国产日本高清| 日韩av无码DVD| 欧美激情成人网| 中国毛片网| 91麻豆精品国产91久久久久| 98精品全国免费观看视频| 激情网址在线观看| 国产综合日韩另类一区二区| 伊人天堂网| 亚洲91精品视频| 一级香蕉视频在线观看| 国产永久无码观看在线| 伊人丁香五月天久久综合| 青青草一区二区免费精品| 国产精品人人做人人爽人人添| 亚洲视频色图| 国产白浆一区二区三区视频在线| 精品久久久久久中文字幕女| 欧美性精品| 亚洲AⅤ无码国产精品| 欧美日韩国产一级| 三上悠亚在线精品二区| 日本中文字幕久久网站| 在线看片免费人成视久网下载| 这里只有精品在线| 99久久国产综合精品2020| 精品久久国产综合精麻豆| 精品国产三级在线观看| 操国产美女| 色亚洲成人| 麻豆精品在线| 欧美精品亚洲精品日韩专区| 日韩精品资源| 中文无码影院| 国产成人精品免费av| 青青极品在线| 欧美国产在线看| 亚洲国产中文欧美在线人成大黄瓜| 制服丝袜在线视频香蕉| 国内熟女少妇一线天| 久久国语对白| 毛片在线播放网址| 凹凸国产熟女精品视频| 国模粉嫩小泬视频在线观看| 在线一级毛片| 国产视频a| 十八禁美女裸体网站| 欧美性精品不卡在线观看| 久久精品aⅴ无码中文字幕| 99视频有精品视频免费观看| 日韩av无码精品专区| 成人午夜视频在线| 人妻出轨无码中文一区二区| 亚洲看片网| 精品久久久久久久久久久| 久久美女精品国产精品亚洲| 国产视频自拍一区| 国产无码精品在线| 久青草网站| 亚洲美女一级毛片| 成人另类稀缺在线观看| 国模在线视频一区二区三区| 久久综合婷婷| 一级毛片免费不卡在线 | 国产手机在线观看|