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

復合圓柱殼沖擊壓縮數值模擬及穩(wěn)定性研究*

2021-12-03 09:07:12謝富佩周中玉谷卓偉
爆炸與沖擊 2021年11期
關鍵詞:結構模型

謝富佩,徐 緋,曾 卓,周中玉,谷卓偉

(1. 西北工業(yè)大學航空學院計算力學與工程應用研究所,陜西 西安 710072;2. 中國航空發(fā)動機集團有限公司湖南動力機械研究所,湖南 株洲 412002;3. 中國工程物理研究院流體物理研究所,四川 綿陽 621999)

復合圓柱殼體結構具有高比剛度、高比吸能以及輕質等優(yōu)點,在航空、航海、汽車工業(yè)及軍事等領域有廣泛應用。在強磁場物理、材料合成等領域,使用了一種利用炸藥爆轟驅動復合圓柱殼內部預置磁場以得到高能磁場的裝置,通常將此裝置稱為柱面內爆磁通量壓縮發(fā)生器(explosive magnetic flux compression generator),簡稱MC-1 裝置[1]。與常見的纖維編織復合圓柱殼不同,由于考慮到復合圓柱殼通電后產生并壓縮磁場的要求,此類復合結構中銅線纖維是連續(xù)平行排列的;在炸藥爆轟產生的沖擊載荷作用下,MC-1 裝置中復合圓柱殼容易發(fā)生界面失穩(wěn)和擾動增長而導致結構失穩(wěn)[2];另外,結構也可能在沖擊壓縮過程中發(fā)生較大的塑性變形而失效,而復合圓柱殼的結構響應及穩(wěn)定性是影響MC-1 裝置性能的重要因素。目前MC-1 裝置在高能量密度物理實驗技術研究領域應用較為廣泛[3],因此對復合圓柱殼體結構在沖擊壓縮下的結構響應模擬及穩(wěn)定性研究有十分重要的意義。

目前國內外對于MC-1 裝置中的復合圓柱殼在沖擊壓縮下的數值模擬研究較少。俄羅斯實驗物理科學研究所于2015 年公布了幾種MC-1 裝置的設計參數,但未公開其參數設計細節(jié);陸禹等[4]使用了引入解析形式狀態(tài)方程組的一維內爆磁流體力學編碼MC11D,對MC-1 裝置工作過程進行了一維條件下的數值模擬;張春波等[5]利用AUTODYN 軟件二次開發(fā)接口建立了內爆壓縮多層密繞螺線管復合圓柱殼三維模型,并實現(xiàn)了周期性邊界條件,使用光滑粒子流體動力學方法對其動態(tài)響應及不穩(wěn)定性開展了數值模擬;趙士操等[6]使用SPH 方法構建了纖維織物的細觀尺度計算模型,使用AUTODYN 軟件對超高速碰撞下纖維增強復合材料模型進行了分析計算;劉軍等[7]對爆轟驅動內壁刻有正弦擾動的金屬鋼殼與內部硅橡膠界面產生不穩(wěn)定性問題進行了數值模擬。由以上研究可以看出,對炸藥爆轟驅動作用下的復合圓柱殼動力學響應及穩(wěn)定性數值模擬問題,由于涉及炸藥爆炸、高速沖擊等過程,傳統(tǒng)網格法面臨大變形、高度非均勻性、運動物質交界面及邊界處理等諸多問題,而粒子法又有計算耗時較長的問題,因此對復合圓柱殼細節(jié)模型的數值模擬常存在較大困難,導致現(xiàn)有的研究缺乏對復合圓柱殼細節(jié)模型的分析手段,且實驗對結構失穩(wěn)現(xiàn)象未提出合理的失穩(wěn)判據,復合圓柱殼穩(wěn)定性的影響因素難以評估。

基于以上缺陷,本文以MC-1 裝置中的復合圓柱殼作為研究對象,使用光滑粒子動力學-有限元方法(smooth particle hydrodynamic-finite element method, SPH-FEM)耦合算法,構建復合圓柱殼二維細節(jié)爆轟模型并對其在炸藥爆轟沖擊下的結構響應進行分析計算,提出一種判斷圓柱殼結構在沖擊壓縮下的穩(wěn)定性評價方式,分析影響MC-1 裝置性能的因素,以期為MC-1 裝置設計及相關實驗研究提供參考。

1 數值計算模型

1.1 模型介紹

本文計算模型的結構參數參照俄羅斯實驗物理科學研究所核中心于2015 年公布的實驗參數,裝置結構如圖1(a)所示,復合圓柱殼外側包裹著一層環(huán)形炸藥,炸藥外側均勻布置48 個起爆點,在同一時刻起爆。裝置中的復合圓柱殼在制備時首先使用一定直徑的漆包銅導線按照規(guī)定的角度螺旋繞制成多匝多層的螺線管內筒體,將繞制的螺線管表面使用環(huán)氧混合劑固定并形成一層絕緣層,然后將從螺線管伸出的銅線覆蓋在絕緣層表面,覆蓋的銅線方向與復合圓柱殼軸線方向平行,最后在其表面使用環(huán)氧混合劑凝固。因此MC-1 裝置中復合圓柱殼由內向外分別為銅線層、內環(huán)氧層、銅線折返層及外環(huán)氧層,如圖1(b)所示。圖1(c)為公布的實驗裝置的結構參數示意圖,考慮到理想情況下銅線層及銅線折返層的厚度為銅線直徑的整數倍,因此計算時將銅線層的厚度調整為2 mm,進行模擬計算時裝置的參數如圖1(d)所示。

圖1 復合圓柱殼示意圖Fig. 1 Diagram of composite cylindrical shell

1.2 數值模型的建立

對在爆轟驅動作用下復合圓柱殼結構響應的模擬計算問題,由于復合圓柱殼三維模型較為復雜,且計算量較大,而復合圓柱殼結構沿軸線方向的截面形狀及尺寸無明顯變化,因此本文選用對復合圓柱殼具有代表性的二維截面模型進行計算。考慮到只使用有限元算法時裝置內圓柱殼結構可能會出現(xiàn)的網格畸變問題以及只使用SPH 算法可能導致的計算規(guī)模較大的問題,本文使用SPH-FEM 耦合算法[8]對爆轟沖擊下的復合圓柱殼進行數值模擬,即對裝置炸藥部分使用FEM 方法劃分網格,對本問題的主要研究對象復合圓柱殼部分使用SPH 方法建立模型,由于兩種方法都具有明顯的拉格朗日特性,因此只需要在軟件中設置接觸域即可實現(xiàn)數據傳遞,完成兩種算法的耦合。在建立模型時,由于復合圓柱殼結構較為復雜,直接在AUTODYN 軟件建模較為困難,本文使用導入外部文件[9]的方式建立模型。首先在CATIA 軟件中建立復合圓柱殼細節(jié)模型并導入到有限元軟件HyperMesh 中劃分三角形網格,在劃分網格時,網格尺寸即為后續(xù)計算的粒子間距,本文將粒子間距取為0.1,導出K 文件后編寫C 程序將模型節(jié)點格式轉化為AUTODYN 軟件可識別的SPH 粒子格式,將其存儲為Formatted 文件后導入AUTODYN軟件,進一步添加有限元炸藥,完成模型邊界和起爆點等的設置后即可完成建模。圖2 為最終在AUTODYN 中建立的15°銅線螺旋角下的復合圓柱殼數值模型。

圖2 15°螺旋角下的爆轟驅動復合圓柱殼模型Fig. 2 Detonation driven composite cylindrical shell model with a spiral angle of 15°

在數值模擬中,復合圓柱殼外側用于提供爆轟壓力的炸藥部件使用JWL 狀態(tài)方程[10]描述:

表1 Comp B 炸藥主要材料參數Table 1 Main material parameters of Comp B explosive

復合圓柱殼中銅線類型為高導無氧銅,其材料參數[11]如表2 所示,塑性部分使用Johnson-Cook 材料模型,材料的屈服應力可定義為:

表2 銅線主要材料參數Table 2 Main material parameters of copper wire

圓柱殼的環(huán)氧部分的材料模型則使用Shock 狀態(tài)方程描述:

式中:U為激波速度;up為粒子速度;c0和s為常數,參數取值參照AUTODYN 材料庫,c0取為2.73 km/s,s取為1.493。

數值計算模型中圓柱殼與炸藥之間的接觸定義為Lagrange/Lagrange 接觸,炸藥外側設置48 個均布起爆點,均在0 μs 時刻同時起爆。

2 結構響應與穩(wěn)定性分析

2.1 數值模型結果驗證與分析

為了驗證建立的數值模型計算結構響應的準確性,在AUTODYN 中建立炸藥爆轟驅動金屬環(huán)氧復合套筒模型進行數值模擬,套筒的內徑為94 mm,外徑為100 mm,沿徑向由內向外分別為厚度為1.5 mm鋼套筒和1.5 mm 的環(huán)氧套筒,外部環(huán)形炸藥的厚度為55 mm。將計算的套筒內側半徑變化曲線與使用已經過驗證的Triangels-MHD 磁流體力學程序[12]計算的結果進行對比,如圖3 所示。

圖3 套筒內側半徑變化曲線對比Fig. 3 Comparison of radius change curves within the sleeve

根據半徑變化對比曲線可以很明顯地看出,除了曲線后半段外,使用兩種方法計算得到的套筒內壁位移響應曲線大部分階段都比較接近,而曲線后半段出現(xiàn)一定差異是因為TriAngels-MHD程序主要是基于磁流體力學方程建立的,在求解過程中考慮到了磁力的影響,在壓縮后期磁力會由于半徑減小而突增,使用AUTODYN 軟件求解時則暫未考慮到后期磁力影響。但是由于磁力影響只占套筒壓縮過程中的極小一部分,因此可以認為本文使用的二維數值計算模型能較好地模擬圓柱殼在炸藥沖擊壓縮下的結構響應。

而由于金屬環(huán)氧套筒的計算模型中未引入初始擾動,整個爆轟過程中各個起爆點的起爆時間相同,起爆點分布也較為密集,使得套筒外側的爆轟載荷沿環(huán)向分布較為均勻,因此在壓縮過程中套筒無較大的環(huán)向擾動,基本保持為一個圓環(huán),如圖4(a)所示。同樣地,復合圓柱殼的變形也有此特點,以15°銅線螺旋角的復合圓柱殼模型為例,在爆轟作用下其不同時刻的變形如圖4(b)所示,從圖中可以看出,在壓縮過程中,復合圓柱殼的環(huán)氧層及銅線折返層有較大的環(huán)向變形,而決定裝置整體性能的銅線層則在壓縮前中期基本保持為環(huán)向相對光滑的狀態(tài),直至壓縮后期才有輕微的環(huán)向擾動。

圖4 圓柱殼在不同時刻的變形Fig. 4 Deformation patterns of cylindrical shell at different times

2.2 穩(wěn)定性判據

在MC-1 裝置的實際應用中,復合圓柱殼受到爆轟載荷后向內收縮是一種高壓高應變率條件下的大變形行為,在此階段中極易發(fā)生界面不穩(wěn)定增長,而復合圓柱殼失穩(wěn)是影響MC-1 裝置性能的重要因素,因此在數值模擬中需引入一種失穩(wěn)判據,以對復合圓柱殼在沖擊壓縮過程中的穩(wěn)定性進行評價。

首先以金屬環(huán)氧套筒為研究對象,圖5(a)為計算的金屬環(huán)氧套筒在沖擊壓縮過程中內壁粒子的速度變化曲線,曲線可大致分為3 段:(1) 炸藥爆轟波尚未到達套筒表面,結構尚無響應,速度為0;(2) 爆轟波到達套筒表面,套筒在爆轟波作用下開始向內收縮,內壁粒子速度急劇增加后逐漸趨于平穩(wěn);(3) 結構出現(xiàn)擾動增長,內壁粒子速度進一步增加并出現(xiàn)不穩(wěn)定增長。根據套筒內壁粒子速度變化歷史曲線可大致地判斷出結構在第3 階段出現(xiàn)失穩(wěn)。為了獲得更加精確的失穩(wěn)時間和失穩(wěn)時對應的壓縮率等數據,對第3 階段的速度曲線進行擬合并對擬合曲線進行求導,圖5(a)左上圖為對速度曲線第3 階段擬合曲線求一階導的導數曲線,導數曲線出現(xiàn)明顯轉折的時間即為套筒結構的失穩(wěn)時間。經過曲線處理后,可判斷出金屬環(huán)氧套筒結構的失穩(wěn)時間大致在15.20 μs 左右,此時套筒的內半徑為18.02 mm,整個結構壓縮率大致為61.66%。

圖5 金屬環(huán)氧套筒失穩(wěn)時間判斷曲線Fig. 5 Judgment curves of instability time of metal epoxy sleeve

為了證明提出的失穩(wěn)判據的合理性,沿金屬環(huán)氧套筒徑向由外向內等間距取點,并將其標記為Node 1~Node 5,提取其沖擊壓縮作用下的速度變化曲線,如圖5(b)所示,可以看到在炸藥爆轟載荷的作用下,金屬環(huán)氧套筒徑向各層粒子在壓縮前期變化規(guī)律基本一致,而壓縮至15 μs 左右時各層粒子的速度曲線開始出現(xiàn)較為明顯的差異,內層粒子的速度增加,而外層粒子速度減小,這種速度差異勢必會導致結構徑向的失穩(wěn),而此失穩(wěn)時間也與前面提出的使用內壁粒子速度曲線變化判斷的失穩(wěn)時間十分接近,由此可以證明提出的失穩(wěn)判據較為合理。

將上述失穩(wěn)判據推廣到復合圓柱殼結構,可以發(fā)現(xiàn)相同的規(guī)律。以15°銅線螺旋角的復合圓柱殼為例,圖6(a)為復合圓柱殼結構內壁粒子的速度變化曲線,也有明顯的三階段特征,圖6(a)左上圖為將復合圓柱殼內壁粒子速度變化曲線的第3 階段曲線擬合求導后的曲線,可以看到曲線在24.76 μs 時開始有突增的趨勢,而此時刻與圖6(b)所示的復合圓柱殼沿徑向各層粒子速度開始出現(xiàn)明顯差異的時間相近,圖6(b)中Node 1~Node 11 指沿復合套筒徑向由外向內等距標記的11 個粒子。綜合以上分析,可以證明本文提出的判斷失穩(wěn)時間的判據對復合圓柱殼結構的穩(wěn)定性分析有一定的參考價值。在后面的研究中,將繼續(xù)使用根據結構內壁粒子速度變化歷史判斷失穩(wěn)時刻的判據對不同條件下的復合圓柱殼結構的穩(wěn)定性進行分析討論。

圖6 復合圓柱套筒失穩(wěn)時間判斷曲線Fig. 6 Judgment curves of instability time of composite cylindrical shell

3 影響復合圓柱殼穩(wěn)定性因素分析

3.1 復合圓柱殼制備缺陷對穩(wěn)定性的影響

由復合圓柱殼的制備過程可知,銅線折返層的銅線是復合圓柱殼內層的銅線層延伸出來的,因此在理想情況下復合圓柱殼的銅線層和銅線折返層銅線數量應相同。在繞制復合圓柱殼結構的內層銅線層時,銅線層橢圓截面的長軸長b計算公式為

式中:D為銅線直徑,θ 為銅線的螺旋角度,同樣以15°螺旋角的復合圓柱殼模型為例,內層銅線層銅線截面的長軸長約為0.965 9 mm,則內層每層銅線的數量ni可表示為

式中:Ri為每層銅線的銅線中心距圓柱殼中心的距離,l為銅線層數,按下式計算

式中:h為模型中銅線層的厚度,對本文中計算的銅線直徑為0.25 mm 的模型,結合式(5)和式(6)估算出銅線層銅線的數量約為3668 根。

而理想條件下,將模型的銅線折返層全部排滿大致需要的銅線數量為5623 根,顯然在實際制備復合圓柱殼時并不能將銅線折返層完全由銅線排滿,在銅線折返層極有可能存在一定的缺陷。因此建立含有缺陷的復合圓柱殼模型來研究復合圓柱殼在制備過程中的缺陷對其穩(wěn)定性的影響,具體做法為在復合圓柱殼折返層建模時隨機刪去若干銅線,本文中含缺陷模型隨機在銅線折返層外側刪去100 根銅線,缺陷部分細節(jié)如圖7(a)所示,模型共有4 處這樣的缺陷。根據計算結果可以發(fā)現(xiàn),由于銅線折返層缺陷部分強度較小,在炸藥爆轟壓縮復合圓柱殼的過程中會導致強度較弱的部分銅線折返層及內環(huán)氧層侵入銅線層,并且隨著爆轟波的運動,侵入部分的深度增加,影響銅線層的變形形態(tài),進而會影響整個結構的穩(wěn)定性,不同時刻含缺陷模型的復合圓柱殼局部形態(tài)如圖7(b)所示,圖7(c)為含缺陷模型復合圓柱殼在24.58 μs 的形態(tài),可以看到其結構內壁環(huán)向不再保持為圓形,圖7(d)為失穩(wěn)時間判斷曲線,含缺陷結構的失穩(wěn)時間為22.63 μs,裝置的穩(wěn)定性明顯降低。以上對含缺陷結構的復合圓柱殼的計算表明:復合圓柱殼結構折返層的制備缺陷會影響結構穩(wěn)定性。

圖7 含缺陷復合圓柱殼變形圖及失穩(wěn)時刻判斷曲線Fig. 7 Deformation diagrams of defective composite cylindrical shell and the judgement curve of instability time

3.2 銅線螺旋角對穩(wěn)定性的影響

復合圓柱殼是非密實、非均勻的結構,銅線的螺旋角度可能是影響其穩(wěn)定性的重要因素[13]??紤]到在實際應用中復合圓柱殼中銅線螺旋角通常在20°以內,故建立螺旋角度分別為0°和10°的復合圓柱殼模型,不同銅線螺旋角的圓柱殼細節(jié)如圖8 所示。除螺旋角度外其余參數及條件均與前面的15°螺旋角模型相同,建模時不同角度模型的粒子數也相近。比較不同銅線螺旋角度下復合圓柱殼的結構響應,使用本文第2 節(jié)提到的速度判斷失穩(wěn)判據對模型的失穩(wěn)時間進行判斷,結合計算的復合圓柱殼位移變化結果,不同螺旋角模型的失穩(wěn)時間及失穩(wěn)時對應的壓縮率比較如表3 所示。

表3 不同螺旋角復合圓柱殼失穩(wěn)時間及壓縮率對比Table 3 Comparison of instability time and compression ratio of composite cylindrical shell with different spiral angles

圖8 不同螺旋角復合圓柱殼細節(jié)模型Fig. 8 Detailed composite cylindrical shell models at different spiral angle

實際上,在建立幾何模型時,可以推算出復合圓柱殼銅線部分的截面積S與銅線的螺旋角度無關的結論:

因此,當銅線直徑相同時,不同螺旋角度銅線在復合圓柱殼結構中的截面積相同,建立的計算模型中粒子數目相近,因此可得到螺旋角度對復合圓柱殼影響不大的結論。

3.3 銅線直徑對穩(wěn)定性的影響

由式(7) 可知復合圓柱殼銅線部分的截面積與銅線直徑有關,另外當銅線直徑發(fā)生變化時,相應地,復合圓柱殼各層的厚度也有可能發(fā)生改變,因此本節(jié)討論銅線直徑對復合圓柱殼的影響,主要討論0.25 和0.50 mm 兩種直徑,銅線螺旋角度均為15°。當銅線直徑為0.50 mm時,可以估算出將2 mm 的銅線層排滿需917 根銅線,這些銅線在銅線折返層只需排列一層,圖9為0.50 mm 銅線直徑下復合圓柱殼模型的結構參數及細節(jié)示意圖。使用前文提到的失穩(wěn)判據,可以判斷對0.50 mm 銅線直徑的復合圓柱殼模型的失穩(wěn)時間大致為23.62 μs,如圖10 所示,此時對應的壓縮率為66.41%,比0.25 mm 直徑的模型失穩(wěn)時對應的壓縮率小。另外當選用0.50 mm直徑的銅線時,由于在復合圓柱殼只有一層銅線折返層,制備時極易出現(xiàn)折返層缺陷而導致裝置穩(wěn)定性變差,圖11 為0.50 mm 銅線直徑下復合圓柱殼銅線折返層僅有一根銅線缺失時在23.20 μs 時的變形形態(tài),可看到僅缺失一根銅線就會導致復合圓柱殼內壁在壓縮過程中有較大的環(huán)向擾動。

圖9 0.50 mm 銅線直徑下復合圓柱殼示意圖Fig. 9 Schematic diagrams of composite cylindrical shell under 0.50 mm copper wire diameter

圖10 0.50 mm 直徑下失穩(wěn)時間判斷曲線Fig. 10 Judgment curve of 0.50 mm model instability time

圖11 t=23.20 μs 時含缺陷復合圓柱殼 變形Fig. 11 Deformation of defective cylindrical shell at t=23.20μs

4 結 論

本文利用AUTODYN 軟件建立了復合圓柱殼二維細節(jié)模型,基于SPH-FEM 耦合算法,對在炸藥爆轟沖擊壓縮下的復合圓柱殼的結構響應過程進行了數值模擬,并根據圓柱殼內壁粒子速度歷史提出了一種判斷失穩(wěn)時間的方法。通過對金屬環(huán)氧套筒爆轟模型的計算,驗證了數值模型計算的準確性,另外由復合圓柱殼及金屬環(huán)氧套筒沿徑向不同位置的速度出現(xiàn)明顯差異的時間與使用提出的失穩(wěn)判據判斷的失穩(wěn)時間一致,證明了提出的失穩(wěn)判據的合理性。

通過對不同數值模型的計算對比,可以發(fā)現(xiàn):在復合圓柱殼制備過程中存在的折返層缺陷對其結構響應及穩(wěn)定性影響較大,因此在制備過程中需要盡可能提高工藝,保證銅線折返層的完整性及均勻性;通過對0°、10°和15°等3 種螺旋角度復合圓柱殼模型的計算,可以發(fā)現(xiàn)復合圓柱殼的螺旋角度對其穩(wěn)定性影響不大;而由于銅線直徑會直接影響復合圓柱殼各層材料的厚度,因此對穩(wěn)定性影響較大,是進行裝置設計及實驗需重點考慮的參數。

猜你喜歡
結構模型
一半模型
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
新型平衡塊結構的應用
模具制造(2019年3期)2019-06-06 02:10:54
論《日出》的結構
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
創(chuàng)新治理結構促進中小企業(yè)持續(xù)成長
主站蜘蛛池模板: 国产噜噜噜| 日韩午夜福利在线观看| 手机精品视频在线观看免费| 青青青视频91在线 | 国产成人亚洲综合A∨在线播放| 免费看a级毛片| 青青草原偷拍视频| 黄网站欧美内射| …亚洲 欧洲 另类 春色| 国产三级视频网站| 久久 午夜福利 张柏芝| 色婷婷色丁香| 天堂岛国av无码免费无禁网站| 久热这里只有精品6| 中文字幕乱码二三区免费| 黄色网站不卡无码| 国产网站免费观看| 好吊妞欧美视频免费| 国产精品亚洲天堂| 456亚洲人成高清在线| 91精品人妻互换| 欧美、日韩、国产综合一区| 91成人在线免费视频| 中文字幕天无码久久精品视频免费| 欧美一级99在线观看国产| 国产欧美日韩免费| 一区二区三区精品视频在线观看| 欧美日韩高清| 国产免费人成视频网| 自慰网址在线观看| 国产91无码福利在线| 91欧美在线| h视频在线播放| 波多野结衣中文字幕一区二区| 久久五月天综合| 欧美黑人欧美精品刺激| 国产第一页亚洲| 91精品啪在线观看国产60岁 | 日韩午夜片| 亚洲av日韩av制服丝袜| 色婷婷在线影院| 国产精品成人观看视频国产 | 天天综合亚洲| 欧美不卡视频在线| 精品精品国产高清A毛片| 蜜芽国产尤物av尤物在线看| 日韩 欧美 国产 精品 综合| 亚洲丝袜第一页| 国产91在线|中文| 在线看片中文字幕| 成人在线综合| 亚洲AⅤ波多系列中文字幕| 91精品免费高清在线| 人妻精品久久久无码区色视| 欧美在线三级| 在线va视频| 国产第一页第二页| 无码中文字幕乱码免费2| a色毛片免费视频| 亚洲欧美极品| 天堂网亚洲系列亚洲系列| 精品亚洲欧美中文字幕在线看 | 色网站在线免费观看| 亚洲欧美另类专区| 欧洲熟妇精品视频| 无码内射中文字幕岛国片| 亚洲色婷婷一区二区| 夜夜操天天摸| 91综合色区亚洲熟妇p| 欧美国产日韩在线观看| 午夜毛片免费观看视频 | 一级毛片在线免费看| 91精品国产无线乱码在线| 国内老司机精品视频在线播出| 亚洲欧美成人影院| 精品国产91爱| h网址在线观看| 久夜色精品国产噜噜| 久久黄色影院| 欧美一级高清片久久99| 无码精品一区二区久久久| 国产99热|