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

槍彈動態擠進阻力理論與實驗研究

2022-10-11 01:31:46許輝黃陳磊王希闊劉坤李忠新吳志林
兵工學報 2022年9期
關鍵詞:有限元變形實驗

許輝, 黃陳磊, 王希闊, 劉坤, 李忠新, 吳志林

(1.南京理工大學 機械工程學院, 江蘇 南京 210094; 2.63850部隊, 吉林 白城 137000)

0 引言

槍械射擊時,槍管處于多物理場相互作用的耦合狀態,高溫高壓火藥燃氣高速流經其內壁,溫度會在極短時間內發生急劇變化,產生較大的非定常、非均勻熱沖擊應力。槍彈擠進作為彈- 槍相互作用的起始階段,具有強沖擊、短歷時、非線性及大變形等特點,對槍管壽命和射擊精度影響尤為顯著。由于現象復雜,影響因素較多,傳統內彈道理論中瞬時擠進假設難以適用。隨著槍械系統的迅速發展,亟待開展槍彈動態擠進理論研究。

近年來,彈- 槍相互作用受到廣泛關注,國內外學者相繼開展了槍彈擠進研究。South等進行了不同速度狀態下小口徑槍彈擠進實驗,獲得了準靜態擠進阻力、周向應變、刻痕深度及質量損失情況。Ritter等借助高速攝像等手段,對小口徑槍彈擠進截短身管初始運動進行了實驗研究,認為難以通過物理模型描述運動過程,可采用高階多項式擬合彈頭的位移- 時間曲線。劉東堯等借助動態壓力測量和高速攝像系統進行了中口徑彈藥的動態擠進實驗,采用對位移數據2階微分的方法,獲得了擠進阻力- 位移曲線。Li等對槍彈高速動態擠進過程進行了有限元仿真研究,獲得了勻速條件下彈頭動態擠進阻力曲線,分析了覆銅鋼和H90黃銅兩種被甲彈頭在不同速度下的動態擠進阻力變化規律。樊麗霞等采用動態顯式算法和網格自適應技術建立了彈頭擠進槍管的有限元模型,研究了鉛芯彈頭的擠進過程,獲得了彈頭膛線刻痕形成和材料流動情況,分析了擠進前后彈頭殼和鉛芯的變形特征。陸野等建立了不同坡膛角下考慮槍管、彈頭結構特性和本構非線性等因素的三維有限元模型,研究了槍管坡膛角度對擠進過程中坡膛受力的影響,獲得了不同坡膛角下擠進阻力隨位移的變化規律。沈超等對某大口徑槍彈擠進過程進行了數值模擬,采用Fortran子程序的方式實現了內彈道過程與顯式有限元方法的迭代數值求解,獲得了槍彈擠進過程中表面形貌和運動姿態隨內膛損傷發展的變化規律。周祥祥等設計了加熱擠進系統,進行了熱槍管準靜態擠進實驗,獲得了不同溫度下擠進阻力曲線,并對連發射擊時不同槍管溫度下的擠進過程進行了數值模擬,結果表明軸向摩擦阻力為擠進阻力的主要成分,導轉力隨槍管溫度升高占比隨之下降。金志明等基于阻尼器原理建立了動態擠進阻力關于彈頭速度、位移相關的函數,給出了動態擠進阻力數學模型。劉國慶等進行了不同坡膛工況下狙擊步槍彈的準靜態擠進實驗,獲得了坡膛錐角對擠進力的影響,采用非線性有限元方法模擬擠進過程,建立了彈頭擠進過程有限元計算模型,分析了準靜態擠進力的組成、坡膛錐角與準靜態擠進力間的關系。安俊斌等運用顯式動力學方法,對某大口徑槍彈擠進坡膛過程進行了數值仿真,獲得了擠進阻力變化規律,分析了彈頭被甲刻痕形成和應力狀態。程斌等構建了彈頭擠進身管的熱力耦合有限元分析模型,研究了彈頭自由行程和初始偏角對擠進過程的影響。蔡翹楚等采用非線性有限元方法建立了某型步槍彈擠進分析模型,分析了彈頭結構參數對擠進過程的影響。周彥煌等根據實驗與理論的分析,提出了準靜態擠進模型與動態擠進模型,分析了結構參數變化對擠進阻力的影響。上述成果為彈- 槍相互作用研究提供了數據支撐,驗證了實驗和數值仿真手段的可行性,具有一定的借鑒意義,但現有工作主要集中于槍彈準靜態擠進,對于動態擠進有待開展更深層次的研究。

本文擬通過理論分析、實驗研究和數值模擬相結合的方法,開展槍彈動態擠進阻力模型研究。采用顯式動力學有限元方法,進行槍彈擠進過程數值模擬,基于完全非彈性碰撞假設,構建動態擠進阻力理論模型,研究成果可充實槍彈設計理論,為輕武器系統優化提供理論參考。

1 動態擠進有限元模型的建立

1.1 槍彈擠進有限元模型

槍彈發射時,膛內火藥燃氣壓力不斷上升,當達到啟動壓力時,彈頭開始運動,逐漸擠進膛線,彈體受到沖擊后產生彈塑性變形,圓柱部形成刻痕,并在槍管的歸正和導轉作用下飛離膛口。根據彈頭擠進特點,有限元建模時引入如下假設:

1)忽略槍管后坐,不考慮溫度應力場;

2)彈頭各組成部分材料各向同性;

3)考慮彈頭彈塑性變形,其屈服強度服從Mises屈服準則;

4)不考慮對流換熱、輻射放熱和摩擦生熱。

彈頭擠進槍管有限元模型如圖1所示。由圖1可知,槍管內膛分為彈膛、坡膛和線膛,彈頭由銅被甲、鉛套和鋼芯組成,初始狀態時彈頭與坡膛之間存在一定間隙,經過一段自由行程后,弧形部與坡膛接觸,開始擠進。

圖1 彈頭擠進槍管有限元模型Fig.1 Finite element model of project and barrel

1.2 網格劃分

依據槍管內膛尺寸與彈頭結構參數建立三維模型,借助有限元前處理軟件HyperMesh劃分網格,槍管和彈頭均選用八節點六面體減縮積分單元(C3D8R),避免單元剪切閉鎖,且單元形狀對計算結果精度影響較小。

網格劃分過程中,槍管坡膛部分對擠進過程影響較大,對其進行網格加密處理,網格密度為整體軸向密度的4倍。彈頭圓柱部作為擠進膛線變形的主要部位,將其沿膛線螺旋方向劃分網格和加密。選取不同尺寸單元對模型進行網格劃分,驗證網格收斂性,劃分網格數量和計算結果如表1所示。由表1可知,隨著模型網格密度增加,數值計算結果趨于穩定,計算時間逐漸增大。方案2可滿足網格收斂性要求,且計算時間較短。因此,選用方案2進行網格劃分。模型劃分網格總數為88.3萬個,其中槍管劃分560 340個單元,被甲劃分223 590個單元,鉛套劃分69 855個單元,鋼芯劃分29 120個單元格,有限元網格模型如圖2所示。

表1 不同網格密度的數值計算結果Table 1 Numerical results of meshes with different densities

圖2 有限元網格模型Fig.2 Finite element meshes

采用Abaqus軟件的動態顯式算法,通過單點積分和基于黏性的沙漏控制,保證大變形計算速度和準確性。

1.3 材料模型

槍管、彈頭被甲、鉛套及鋼芯的材料分別為中碳低合金鋼30SiMn2MoVA、H90黃銅、純鉛、碳素結構鋼Q235,擠進過程受到應變率、損傷等多種因素的影響,彈頭產生彈塑性大變形,選用Johnson-Cook塑性模型進行描述。Johnson-Cook塑性模型中,Mises屈服應力為塑性應變、應變率和溫度相關的函數,可表示為

(1)

Johnson-Cook模型通過損傷參數描述損傷度,可表示為

=∑Δ

(2)

式中:為損傷參數(其值在0~1范圍內,初始時=0,材料發生失效時=1);Δ為等效塑性應變增量;為材料失效應變。

材料失效應變可表示為

(3)

擠進過程中,彈頭材料斷裂失效可通過Abaqus軟件中的單元刪除進行模擬,計算所用材料模型參數如表2所示。

表2 材料模型參數[9,20,22]Table 2 Material parameters of numerical simulation[9,20,22]

1.4 施加載荷和邊界條件

發射時,火藥燃氣作用于彈底,使彈頭向前運動,并隨著彈頭擠進過程不斷變化,通過實彈射擊測試獲得彈底壓力變化曲線。將沿槍管軸線方向的彈底壓力作為彈頭擠進數值模擬的主動載荷,槍管尾端定義為完全約束狀態,約束其全部自由度。

被甲外表面與槍管內膛之間、被甲內表面與鉛套之間、鉛套與鋼芯之間定義為面- 面接觸,擠進過程采用庫倫摩擦模型,摩擦系數設置為002。

2 實驗驗證與仿真結果分析

2.1 實驗方法

為驗證有限元模型的正確性,進行槍彈擠進實驗。實驗采用58 mm短槍管彈道槍(短槍管彈道槍的槍管長度為69 mm,常規彈道槍的槍管長度為520 mm,其余結構均相同)作為發射裝置,共射擊 5發,選取10A式58 mm步槍彈作為實驗用彈,并在彈頭弧形部裝有細長桿狀彈帽(為降低彈帽對彈頭運動狀態的影響,設計了聚合物材質的輕質彈帽,質量為016 g,占彈頭質量的35),以便記錄彈頭擠進位移。發射裝置和槍彈如圖3所示。借助高速攝像機(日本NAC公司產Memrecam ACS-1 M60E高速攝像機,分辨率1 280×384,幀頻200 000幀s)拍攝槍彈出膛姿態,壓力傳感器(瑞士Kistler公司產6215,量程600 MPa,線性度誤差045)安裝于距膛底16 mm位置以獲得彈殼體部膛壓,光電靶距離膛口2 m測量出膛速度,光幕和光源距膛口側03 m處(光源置于光幕后),高速攝像機置于距膛口側05 m處,確保高速攝像機鏡頭軸線與彈道在同一平面內,且與彈道方向垂直,擠進實驗原理和場景如圖4所示。

圖3 發射裝置和實驗用彈藥Fig.3 Launcher and bullet for the test

圖4 實驗原理和場景圖Fig.4 Schematic diagram and test setup

2.2 實驗結果與仿真計算對比分析

由于測量難度較大,無法通過實驗測得擠進阻力、變形力和滑動摩擦阻力,僅獲得膛壓、速度和彈頭刻痕。所測彈殼體部膛壓曲線如圖5所示,與長槍管彈道槍所測膛壓進行對比,彈頭飛離短槍管膛口壓力快速泄露,其壓力曲線在0445 ms處開始下降。以發射藥開始燃燒時刻為計時起點、彈頭出膛時刻為計時終點,彈頭出膛速度實測值為3488 m/s,運動過程如表3所示。

表3 實驗中不同時刻彈頭的位移照片Table 3 Photos of projectile displacement at different moments during the experiment

圖5 彈膛內的壓力時間曲線Fig.5 Pressure-time curve in the chamber

圖6給出了彈頭位移隨時間變化曲線。由圖6可見:仿真所得彈頭完全嵌入線膛時位移為234 mm,實測值為236 mm,位移誤差為02 mm,相對誤差為084。速度隨時間變化曲線如圖7所示。由圖7可知,仿真所得完全嵌入線膛時彈頭速度為2192 m/s,出膛速度為3321 m/s,出膛速度實測值為3488 m/s,出膛速度誤差為167 m/s,相對誤差為48。通過對比可知,仿真結果與實測值誤差較小,一致性較好,驗證了本文所建立的槍彈動態擠進有限元模型的正確性。

圖6 擠進過程位移- 時間仿真曲線與實驗曲線對比Fig.6 Comparison of displacement-time curves obtained from simulation and experiment

圖7 擠進過程速度- 時間仿真曲線與實驗曲線對比Fig.7 Comparison of velocity-time curves obtained from simulation and experiment

槍彈擠進過程中應力和應變云圖如圖8所示。由圖8可知,擠進時彈頭殼受到擠壓產生塑性變形,隨著擠進深度增加,彈頭被膛線擠壓導致的材料變形與失效形成了刻痕,在0356 ms時圓柱部完全擠進線膛。

圖8 擠進過程中彈頭的應力應變云圖Fig.8 Stress and strain nephograms of the projectile during engraving

彈頭完全擠進線膛后,彈頭殼刻痕數值模擬與實驗對比如圖9所示。由圖9可知,實驗回收的彈頭刻痕清晰規整,與仿真計算所得刻痕形貌一致,具體刻痕尺寸對比如表4所示(表4中實驗均值為5發彈頭刻痕測量結果的平均值)。對比可知刻痕長度、寬度實驗值與仿真值相對誤差分別為344、179,進一步證明了有限元計算模型的正確性。

表4 彈頭刻痕尺寸實驗值與仿真值對比Table 4 Comparison of notch size between experiment and simulation

圖9 彈頭刻痕對比Fig.9 Comparison of notches on the projectile

彈頭動態擠進阻力主要包括克服彈頭材料變形產生的變形力和軸向滑動摩擦力,受限于當前測試條件,無法通過實驗直接測得動態擠進阻力,可借助數值模擬獲得其變化規律。圖10給出了彈頭動態擠進阻力變化規律。由圖10可見:隨著彈頭膛內運動位移的增加,變形力迅速增大,在位移為123 mm時彈頭完全嵌入坡膛,動態擠進阻力出現最大值1 567 N,變形力達到峰值1 026 N,隨后急劇減小直至趨于0 N,滑動摩擦力緩慢上升至541 N,并趨于穩定值;位移為123~191 mm階段時,彈頭圓柱部開始從坡膛向線膛擠進,變形力出現一個緩降的平臺期,其大小與滑動摩擦力接近;在位移大于232 mm時,彈頭完全嵌入線膛,擠進阻力趨于穩定,變形力下降至0 N,滑動摩擦力起主要作用;由于數值模擬所得變形力中包含導轉側力,其最終殘余值主要體現為導轉側力。圖11為彈頭體積變化和變化率與位移關系曲線。圖11結合圖10發現,體積變化率和變形力曲線變化趨勢、拐點位置較為一致,由此可得彈頭動態擠進過程中變形力與擠進膛線的體積變化率存在強關聯性。

圖10 有限元仿真得到的動態擠進阻力曲線Fig.10 Dynamic engraving resistance curve obtained from finite element simulation

圖11 擠進過程中彈頭的體積變化Fig.11 Volume change of the projectile during engraving process

3 擠進阻力理論模型

3.1 基本假設

槍彈動態擠進具有非線性、瞬時性、大變形等特征,為準確描述彈- 槍相互作用,在數值模擬和實驗研究基礎上,構建槍彈動態擠進阻力模型。建模前,引入如下假設:

1)彈頭擠進變形過程視為完全非彈性碰撞;

2)假設使彈頭變形的接觸面為坡膛陽線處,其他接觸面的接觸應力為定值;

3)用庫倫摩擦定律描述動態擠進過程的高速滑動摩擦力,假設摩擦系數為定值;

4)忽略彈頭轉動及其導轉側力。

3.2 模型構建

彈頭動態擠進過程中,高溫高壓火藥燃氣產生的推力推動彈頭向前運動,彈頭在任意時刻所具有的速度為、加速度為,在Δ時間內的位移為Δ。彈頭因擠壓變形產生的可視體積變化為Δ,受擠壓部分的體積向其臨近空間移動,迫使臨近體積參與變形,參與變形部分的等效體積為Δ,彈頭抵抗變形所受到的變形力為,與身管內壁相互作用產生的滑動摩擦力為。擠進過程受力情況如圖12所示。圖12中,為彈頭被甲與槍管內壁的接觸面積,為變形分界面面積。擠進過程中,變形力和滑動摩擦力在彈頭徑向和周向上均存在分量,如圖13所示。圖13中,為在彈頭運動方向的分量,為在垂直于彈頭運動方向的分量,為在彈頭運動方向的分量,為垂直于彈頭運動方向的分量。

圖12 擠進過程受力情況Fig.12 Force in the engraving process

圖13 擠進過程簡化力學模型Fig.13 Simplified mechanical model of the engraving process

321 變形力

在彈頭運動方向上,將動態擠進過程簡化為質量塊與固定壁完全非彈性碰撞,如圖14所示。在時刻,參與碰撞的等效質量體與固定壁間發生完全非彈性碰撞,碰撞后等效質量體速度衰減至0 m/s,將動能轉化為內能。碰撞過程中,彈頭質量主體與等效質量體之間的變形分界面上作用力大小為彈頭殼材料的屈服強度與面積的乘積。

圖14 動態擠進變形簡化力學模型Fig.14 Simplified mechanical model of the dynamic engraving process

碰撞過程滿足動量守恒,即

(4)

式中:Δ為碰撞時間;為彈頭殼材料密度;為動態系數。

將(4)式寫成微分形式,整理可得

(5)

式中:()為彈頭被身管擠壓造成的可視形變體積,是與彈頭、身管結構尺寸相關的函數。

擠進過程中,彈頭和槍管的等效軸向接觸面積可表示為

(6)

式中:為靜態系數。

將(6)式代入(5)式,可得

(7)

動態接觸應力可表示為

(8)

式中:等式右端第1項為動態阻力項,第2項為靜態強度項。

準靜態擠進情況下,(7)式可表示為

(9)

該形式與文獻[15,19]中給出的靜態擠進阻力形式類似,適用于準靜態擠進過程。可通過準靜態擠進實驗確定其靜態系數,再結合動態擠進實驗和數值模擬確定動態系數。

322 滑動摩擦力

在彈頭的運動方向上受到的滑動摩擦力,與接觸面的正壓力相關。由于身管的坡膛角較小,根據投影關系,可將滑動摩擦過程簡化為滑塊與固定壁的平面摩擦模型,如圖15所示。

圖15 高速滑動摩擦簡化力學模型Fig.15 Simplified mechanical model of high-speed sliding friction

在時刻,彈頭位移為,彈頭與身管的接觸面積為,坡膛陽線處接觸面積為,該處法向動態接觸應力與相關,其他接觸面積的平均接觸應力為定值,彈頭擠進后的平均接觸應力與被甲材料的屈服強度有關,可表示為

(10)

式中:、分別為法向接觸應力系數和平均接觸應力系數。

若高速滑動摩擦系數為(與速度相關),則滑動摩擦力可表示為

=+·(-)

(11)

考慮到坡膛角很小,cos≈1,則

()=()·cos≈()

(12)

式中:()為擠進過程中彈頭和身管的接觸面積,與彈頭和身管的相對位置有關。

根據幾何關系可知,與之間存在對應關系,引入面積修正系數,則

(13)

聯立(10)式、(11)式、(12)式及(13)式,可得

(14)

綜上所述,彈頭動態擠進過程中擠進阻力可表示為

(15)

3.3 模型計算與分析

331 理論計算結果

理論模型中動態摩擦系數取=002,通過多次擬合計算確定模型系數,如表5所示。()、()的取值通過彈頭和身管幾何尺寸計算獲得,如圖16所示。理論模型計算所使用的-數據來源于動態擠進實驗。

表5 理論模型參數Table 5 Coefficients of the theoretical model

圖16 彈頭擠進過程Vs(x)、As(x)與x之間的關系Fig.16 Relationship between Vs(x), As(x), and x of projectile during the engraving process

理論計算所得動態擠進阻力曲線如圖17所示。由圖17可知,開始擠進時(=46 mm處),彈頭速度為689 m/s,隨著彈頭速度和體積變化率的增加,擠進阻力呈指數式上升。彈頭完全嵌入坡膛時(=123 mm處),速度為1586 m/s,動態擠進阻力、變形力分別達到最大值1 435 N、873 N,摩擦力也接近最大值。同時,彈頭圓柱部尾端開始嵌入坡膛起始位置,圓柱部前端距線膛起始位置約1 mm(圓柱部的長度小于坡膛長度)。由圖16可知,該時刻彈頭體積變化率最大,約為0763 mm,彈頭圓柱部全部進入坡膛后,體積變化率隨之出現陡降。

圖17 理論模型計算得到的彈頭擠進阻力Fig.17 Projectile engraving resistance obtained from the theoretical model

彈頭圓柱部開始進入線膛后(=133 mm處),體積變化率近似呈線性減小,由于彈頭速度繼續增大,變形力將緩慢下降,摩擦力基本趨于穩定,擠進阻力總體呈緩慢下降趨勢。彈頭完全嵌入線膛時(=236 mm處),動態擠進阻力趨于穩定,約為579 N。由于未考慮彈頭的導轉側力,變形力在擠進結束后趨于0 N,擠進阻力主要表現為摩擦力作用。

由此可見,彈頭擠進速度和擠進過程中的體積變化率是影響槍彈動態擠進阻力的關鍵因素。

332 理論計算與數值仿真結果分析

圖18給出了理論模型計算結果與仿真結果對比。由圖18可見,模型可得彈頭在位移為123 mm時,擠進阻力出現最大值1 435 N,仿真值為1 567 N,相對誤差為85。擠進阻力變化規律理論模型計算與數值仿真結果具有較好的一致性,由于理論模型忽略彈頭導轉側力,最終穩定值略有差異。

圖18 動態擠進阻力理論模型計算結果與仿真對比曲線Fig.18 Comparison of dynamic engraving resistance of the projectile obtained from simulation and theoretical model

圖19為變形力的理論計算結果與仿真值對比曲線。由圖19可知,模型可得彈頭在位移為123 mm時,變形力出現最大值873 N,仿真值為1 026 N,相對誤差為149。由于擠進過程中體積變化率的影響,在位移為123 mm處曲線出現陡降。

圖19 變形力理論模型計算結果與仿真對比曲線Fig.19 Comparison of deformation forces obtained from simulation and theoretical model

圖20為有限元仿真與理論模型計算得到的滑動摩擦力對比曲線。由圖20可知,理論計算所得在位移為173 mm時,滑動摩擦力出現最大值612 N,仿真值為615 N,相對誤差為05,動態擠進過程中的摩擦阻力主要受彈頭與身管之間接觸面積的影響,擠進結束后摩擦力趨于定值。

圖20 滑動摩擦力理論計算結果與仿真對比曲線Fig.20 Comparison of sliding friction obtained from simulation and theoretical model

有限元計算所得變形力具有一定的震蕩特性,滑動摩擦力較為平滑。因此兩種模型給出的變形力最大值計算結果存在較大誤差(149),滑動摩擦力最大值計算結果誤差較小(05)。擠進阻力最大值的誤差體現為變形力和滑動摩擦力計算結果的疊加作用,最終相對誤差(85)介于二者之間。綜上所述,通過對比可知,基于完全非彈性碰撞假設的槍彈動態擠進阻力理論模型的計算結果與仿真值誤差較小,一致性較好,可更為準確地描述彈- 槍相互作用過程。

4 結論

本文為揭示彈- 槍相互作用機理,通過理論分析、實驗研究和數值模擬相結合的方法,開展槍彈動態擠進阻力模型研究。采用顯式動力學有限元方法,借助Abaqus軟件模擬槍彈擠進過程。利用多參數同步測量技術進行了擠進實驗,驗證了數值模擬的有效性。基于完全彈性碰撞假設,考慮彈頭材料靜態強度、動態變形力和高速運動下摩擦力的作用,結合彈頭運動、彈- 槍結構參數構建了動態擠進阻力模型,分析了擠進阻力關鍵影響因素,并與數值模擬進行對比。得出以下主要結論:

1)通過分析動態擠進阻力曲線可知,動態擠進阻力主要由彈頭變形力和軸向滑動摩擦力組成,變形力與擠進膛線的體積變化率存在強關聯性。

2)隨著擠進位移增大,變形力先增大后減小,在圓柱部完全嵌入坡膛時達到最大值,摩擦力增大至峰值后趨于穩定。

3)彈頭擠進速度和擠進過程中的體積變化率為影響動態擠進阻力的關鍵因素。

4)動態擠進理論模型計算結果與數值模擬結果一致性較好,能準確描述彈- 槍相互作用過程。

猜你喜歡
有限元變形實驗
記一次有趣的實驗
談詩的變形
中華詩詞(2020年1期)2020-09-21 09:24:52
做個怪怪長實驗
“我”的變形計
例談拼圖與整式變形
會變形的餅
NO與NO2相互轉化實驗的改進
實踐十號上的19項實驗
太空探索(2016年5期)2016-07-12 15:17:55
磨削淬硬殘余應力的有限元分析
基于SolidWorks的吸嘴支撐臂有限元分析
主站蜘蛛池模板: 人妻91无码色偷偷色噜噜噜| 波多野结衣久久高清免费| 99在线观看国产| 亚洲性影院| 亚洲精品动漫| julia中文字幕久久亚洲| 国产高清无码第一十页在线观看| 在线精品亚洲一区二区古装| 丰满少妇αⅴ无码区| 国产成年女人特黄特色大片免费| 日韩AV无码免费一二三区| 国内熟女少妇一线天| av午夜福利一片免费看| 黄色网址手机国内免费在线观看| 操国产美女| 少妇精品在线| 在线欧美a| 亚洲一级无毛片无码在线免费视频 | 国产成人盗摄精品| 国产91小视频在线观看| 欧美a级完整在线观看| a级毛片免费看| 亚洲丝袜第一页| 一本大道香蕉中文日本不卡高清二区| 91www在线观看| 精品国产成人高清在线| 久久久波多野结衣av一区二区| 91蜜芽尤物福利在线观看| 国产成人AV综合久久| 国产a网站| 欧美日韩一区二区三| 成人免费午间影院在线观看| 亚洲视频色图| 国产第一页第二页| 蜜桃视频一区| 91在线国内在线播放老师| 中文字幕天无码久久精品视频免费 | 米奇精品一区二区三区| 热99精品视频| 91精品福利自产拍在线观看| 久久国产精品娇妻素人| 88国产经典欧美一区二区三区| 乱系列中文字幕在线视频| 中文字幕1区2区| 亚洲天堂精品在线| 玖玖精品视频在线观看| 真实国产精品vr专区| 免费全部高H视频无码无遮掩| 亚洲一区黄色| 白浆免费视频国产精品视频| 国产97公开成人免费视频| 亚洲中文精品久久久久久不卡| 日日碰狠狠添天天爽| 无码粉嫩虎白一线天在线观看| 毛片视频网址| 呦女亚洲一区精品| 老司机久久精品视频| 亚洲色精品国产一区二区三区| 2021国产精品自产拍在线观看| 国产丝袜无码一区二区视频| 亚洲精品自拍区在线观看| 五月婷婷精品| 欧美亚洲中文精品三区| 天堂在线视频精品| 中文精品久久久久国产网址| 国产永久无码观看在线| 91成人在线免费视频| 国产精品部在线观看| AV片亚洲国产男人的天堂| 97久久超碰极品视觉盛宴| 国产一区二区影院| 精品视频福利| 久久综合亚洲鲁鲁九月天| 无码中文字幕精品推荐| 成人年鲁鲁在线观看视频| 91色国产在线| 日本高清视频在线www色| 亚洲侵犯无码网址在线观看| 亚洲欧美不卡| 国产不卡国语在线| 亚洲精品色AV无码看| 香蕉精品在线|