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

基于Johnson-Cook本構模型的高強度裝甲鋼動態力學性能參數標定及驗證

2022-08-27 09:49:08張杜江趙振宇賀良任建偉強鷺升周貽來
兵工學報 2022年8期
關鍵詞:模型

張杜江, 趙振宇, 賀良, 任建偉, 強鷺升, 周貽來

(1.南京航空航天大學 機械結構力學及控制國家重點實驗室, 江蘇 南京 210016;2.南京航空航天大學 多功能輕量化材料與結構工業和信息化部重點實驗室, 江蘇 南京 210016;3.內蒙古第一機械集團有限公司 特種車輛及其傳動系統智能制造國家重點實驗室, 內蒙古 包頭 014030)

0 引言

簡易爆炸裝置(IED)和地雷等爆炸物產生的爆炸沖擊載荷嚴重威脅車輛乘員的生命安全,優化車體結構設計和選擇先進裝甲防護材料是提升車輛防護性能的主要措施,而車底防護結構的構型設計一旦確定,材料選擇成為另一研究重點。在集中、均布等外部載荷作用下,Jones理論分析了理想剛塑性均質梁/板的動態響應,給出了其無量綱撓度與無量綱沖量之間的關系。從無量綱沖量的表達式中可以看出,選擇屈服強度高的材料可降低無量綱沖量,從而降低均質板無量綱撓度,提升抗爆性能。隨著我國軍用裝備的快速發展,對裝甲鋼的綜合性能提出了更高的要求,某裝甲鋼具有較高的屈服強度(大于1 400 MPa),為我國裝甲車輛的防護結構設計提供了新的材料選擇。

通常,采用大量試驗對防護結構的構型或者材料選擇進行驗證,試驗成本較大,研發周期也較長,而采用數值仿真技術先期進行分析和優化設計,可節約試驗成本與時間。對防護結構的動態響應過程進行數值分析,首先需要確定結構所使用材料的本構模型及相關材料參數。其中,由于Johnson-Cook本構模型(以下簡稱J-C本構)考慮了金屬材料的硬化效應、應變率效應和溫度效應,并得到了充分的實驗驗證,在國內外實現了廣泛應用。胡昌明等利用溫度為25~300 ℃和應變率為10~10s的拉伸應力- 應變曲線,通過擬合獲得了45號鋼的J-C本構參數;郭子濤等研究了Q235鋼在常溫至900 ℃的準靜態和動態壓縮及拉伸力學性能,修正了J-C本構模型中的溫度軟化項,并利用Taylor撞擊實驗和數值仿真驗證了其動態本構關系;Senthil等通過12.7 mm和7.62 mm穿甲燃燒彈侵徹20 mm厚的7075鋁合金板,驗證了獲得的7075鋁合金J-C本構參數;高玉龍等通過溫度為-50~300 ℃的單軸拉伸試驗、應變率為10~10s的單軸拉伸和壓縮試驗,使用試驗和仿真相結合的方法,進行了0°剪切、45°剪切、單軸拉伸、單缺口拉伸和雙缺口拉伸試驗,得到了6008鋁合金J-C本構的硬化及損傷斷裂參數?;趯嶒灉y量,6008鋁合金、7075鋁合金、工業純鐵、304不銹鋼、45號鋼、Weldox 900E、Armox 500T、鈦合金等典型金屬材料的J-C本構參數如表1所示。此外,由于大型商業有限元軟件(如Abaqus、LS-DYNA等)均內置了J-C本構,該本構已廣泛應用于防爆結構設計、汽車耐撞性檢驗、飛機防鳥撞設計等領域。

表1 典型金屬材料的J-C本構模型參數

針對我國某裝甲鋼材料抗沖擊力學性能及J-C本構模型參數的研究,目前未見公開報道。為準確模擬該裝甲鋼在爆炸載荷下的響應過程,獲得其動態本構關系參數是必要的;一般而言,由于裝甲車的防護結構設計應確保防護結構在炸藥爆炸后不破裂,研究不涉及J-C本構模型的斷裂參數。

本文在常溫和高溫環境下開展某裝甲鋼的準靜態拉伸試驗研究,同時在常溫環境下,采用分離式霍普金森壓桿測試系統開展不同應變率下的動態壓縮試驗研究?;贘-C本構模型并結合材料實驗結果數據,擬合得到該裝甲鋼的J-C本構模型參數;再次,采用輕氣炮和泡沫鋁彈丸對該裝甲鋼均質梁開展沖擊試驗研究。分別采用J-C本構模型和理想彈塑性模型進行有限元仿真計算,并將數值結果與沖擊試驗進行對比分析。

1 本構模型

Johnson和Cook等基于材料各向同性假設,考慮材料的硬化效應、應變率效應和溫度效應的影響,提出了著名的J-C本構模型。J-C本構模型將應力- 應變關系表示為

(1)

(2)

(3)

分別開展常溫、高溫條件下的準靜態拉伸試驗以及分離式霍普金森壓桿試驗研究,然后基于試驗數據擬合獲得某裝甲鋼的J-C本構模型參數,具體步驟如下:

1)常溫條件下,開展準靜態拉伸試驗,根據獲得的真實應力- 應變數據可得屈服強度。隨后,考慮應變硬化系數和應變硬化指數,忽略應變率和溫度的影響,此時(1)式可簡化為

(4)

基于(4)式,對真實應力- 應變曲線塑性段進行數據擬合,可得應變硬化系數和應變硬化指數。

(5)

基于(5)式,對不同溫度下獲得的真實應力- 應變曲線,取多個相同塑性應變值的點進行擬合,得到多個熱軟化系數,然后取平均值。

3)考慮應變率對屈服應力影響時,無量綱溫度為0,則應力和無量綱應變率的關系簡化為

(6)

基于(6)式,對于不同應變率下獲得的真實應力- 應變曲線,取多個相同塑性應變值的點進行擬合,得到多個應變率硬化系數,然后取平均值。

2 材料性能測試

2.1 常溫下準靜態性能試驗

根據國家標準GB/T 228.1—2010金屬材料拉伸試驗第一部分室溫試驗方法,設計某裝甲鋼的拉伸試驗樣件,如圖1所示。試驗樣件處于彈性拉伸階段時,拉伸速率取為0.216 mm/min,此時應變率為1×10s;材料進入塑性拉伸階段后,拉伸速率取為3 mm/min。圖2給出室溫為20 ℃時3次重復試驗后得到的真實應力- 應變曲線。由圖2可見,某裝甲鋼在室溫條件下沒有明顯的屈服點。因此,選擇試驗樣件產生0.2%塑性應變時的真實應力作為屈服強度,將3次試驗得到的屈服應力進行平均,得到該裝甲鋼的屈服應力為1 458 MPa。同時,對圖2中真實應力- 應變曲線的彈性階段進行線性擬合,得到裝甲鋼的彈性模量為220 GPa。

圖1 拉伸試驗樣件尺寸Fig.1 Dimensions of tensile specimens

圖2 室溫(20 ℃)條件下真實應力- 應變曲線Fig.2 True stress versus true strain curve at room temperature (20 ℃)

2.2 高溫下準靜態性能試驗

在溫度為180 ℃、260 ℃、370 ℃和550 ℃條件下,對圖1的試樣進行3次準靜態拉伸試驗。圖3給出了室溫和高溫條件下的真實應力- 應變曲線,圖4所示為屈服強度隨溫度變化曲線圖。由圖3和圖4可見:從室溫20 ℃變為180 ℃時,屈服強度下降相對緩慢,溫度為180 ℃時,流動應力沒有出現通常認為的下降,而是流變應力大于室溫時的流變應力,Wang等在2015年首次將這一現象命名為“第三型應變時效”,這是由運動位錯與擴散的溶質原子的相互作用引起的。位錯在障礙前等待時,溶質原子向位錯擴散,在位錯周圍形成溶質原子氣團,對運動位錯“釘扎”,阻礙了位錯的運動,在宏觀上表現為金屬流動應力增大。溫度大于180 ℃時,該裝甲鋼的屈服強度明顯降低,且呈線性下降,如圖4所示。

圖3 不同溫度下真實應力- 應變曲線Fig.3 True tress versus true strain curve at different temperatures

圖4 屈服強度隨溫度變化曲線Fig.4 Yield stress versus temperature

2.3 動態性能測試

根據國家軍用標準GJB 8799—2015金屬材料動態壓縮試驗方法,采用尺寸為6 mm×5 mm的圓柱形樣件,開展不同應變率下的分離式霍普金森壓桿沖擊試驗,如圖5所示。

圖5 分離式霍普金森壓桿試驗樣件Fig.5 Specimens for Split-Hipkinson pressure bar tests

圖6給出了6種不同應變率下的動態壓縮應力- 應變曲線,可見該裝甲鋼的動態屈服強度大于其準靜態屈服強度,且隨著無量綱應變率的增加而增大。

圖6 不同應變率下真實應力- 應變曲線圖Fig.6 True stress versus strain curves at different strain rates

3 參數擬合

采用(4)式對準靜態真實應力- 應變曲線的塑性階段進行擬合,如圖7所示,擬合結果與試驗數據吻合較好,由此得應變硬化系數為1 408 MPa,應變硬化指數為0.405。

圖7 真實應力隨塑性應變變化曲線Fig.7 True stress versus plastic strain

根據(5)式,對塑性應變為0、0.005、0.010、0.015、0.020、0.025和0.030時的應力和無量綱溫度進行擬合,得到熱軟化參數分別為0.558、0.632、0.672、0.714、0.705、0.697、0.689,平均后為0.667。圖8給出了塑性應變為0和0.030時,應力隨無量綱溫度變化曲線。

圖8 應力隨無量綱溫度變化趨勢Fig.8 Stress versus dimensionless temperature

根據(6)式,對塑性應變為0、0.005、0.010、0.015、0.020、0.025和0.030時的應力和無量綱應變率進行擬合,得到熱軟化參數分別為0.014 0、0.007 2、0.006 1、0.005 6、0.005 4、0.005 4、0.005 5,平均后為0.007。圖9給出了塑性應變為0和0.030時,應力隨對數應變率變化曲線。

圖9 應力隨對數應變率的變化趨勢Fig.9 Stress versus log strain rate

綜上所述,某裝甲鋼的J-C本構模型參數列入表2。

表2 某裝甲鋼J-C本構模型參數

4 材料參數驗證

4.1 試驗裝置與樣件

一級輕氣炮沖擊試驗裝置主要由支撐底座、壓氣室、炮筒、測試艙、回收艙及高速攝影機組成,如圖10所示。由圖10可見;在壓氣室中充入具有一定壓力的氮氣,高壓氣體釋放后,泡沫鋁彈丸受到高壓氣體作用,在炮管中加速并獲得一定的初速度后沖擊靶板(某裝甲鋼均質梁);采用線切割工藝從閉孔泡沫鋁上切出直徑為58 mm,長度為100 mm的彈丸。圖11所示為測試艙中工裝夾具示意圖,當泡沫鋁彈丸運動到測試艙時,通過高速攝影機捕捉彈丸的運動過程和均質梁的響應過程。如圖12所示,均質梁的長度2=370 mm,寬度2=60 mm,名義厚度=5 mm (實測為5.13 mm)。試驗時,采用尼龍66扎帶將均質板樣件固定在夾具上。考慮到尼龍66的抗拉強度約為71.5 MPa,而扎帶的橫截面積較小,為1.3 mm×4.7 mm,經計算得到扎帶斷裂的臨界載荷約為0.43 kN。由于臨界載荷較低,在仿真時忽略尼龍扎帶的影響。試驗結束后,通過高速攝影機拍攝到的照片測量樣件的變形過程,并測量樣件的殘余撓度。

圖10 輕氣炮沖擊試驗裝置Fig.10 Schematic diagram of the impact test setup

圖11 測試艙內的試樣及工裝夾具Fig.11 Schematic diagram of beam specimen fixed in the test cabin

圖12 某裝甲鋼均質梁樣件示意圖Fig.12 Geometry and dimensions of the monolithic beam made of the studied armor steel

為了在后續的仿真分析中模擬泡沫鋁彈丸,開展了閉孔泡沫鋁(密度370 kg/m)準靜態單軸壓縮試驗,樣件直徑20 mm、高度50 mm,試驗時的名義應變率為1×10s;圖13給出了實測的工程應力- 應變曲線。圖13中,為泡沫鋁的名義致密應變,為平臺應力。

圖13 單軸壓縮下閉孔泡沫鋁工程應力和能量吸收率隨工程應變變化曲線Fig.13 Uniaxial compressive engineering stress and energy absorption efficiency versus engineering strain curve of closed-cell aluminum foam

為確定,定義能量效率參數如下:

(7)

(8)

對圖13的工程應力- 應變曲線進行多項式擬合,根據(7)式和(8)式確定能量效率曲線,可得=059,=57 MPa。

4.2 有限元模型

基于輕氣炮沖擊試驗,采用有限元軟件Abaqus對泡沫鋁彈丸沖擊均質梁的過程進行數值模擬。圖14給出的有限元模型包括泡沫鋁彈丸、某裝甲鋼均質梁和固支的工裝夾具。其中,均質梁長度為370 mm,寬度為60 mm,厚度為5.13 mm(試驗樣件的實際厚度)。為避免彈丸和均質梁之間在沖擊初始即出現網格穿透,建模時在彈丸和均質梁之間設置厚度為0.1 mm的間隙。為避免在后續計算過程中出現網格穿透,在有限元模型中設置通用接觸;均質梁和完全固定的工裝夾具之間的摩擦系數取為0.3。采用Crushable Foam本構模擬泡沫鋁:根據單軸壓縮下閉孔泡沫鋁工程應力- 應變變化曲線中的彈性階段進行線性擬合,擬合得到的斜率為泡沫鋁彈性模量,大小為210 MPa;由于泡沫鋁的塑性泊松比為0,體積應變=()=,其中為準靜態壓縮工程應變,為泡沫鋁的初始高度,為泡沫鋁的初始橫截面積,因此準靜態壓縮工程應變與體積應變是相等的,本構中的應力- 體積應變關系由圖13中數據給出。Deshpande等通過霍普金森壓桿試驗得到了泡沫鋁在不同應變率下的應力應變曲線,結果表明泡沫鋁對應變率的變化不敏感,高華等研究了多次沖擊下泡沫鋁的動態壓縮力學性能。在Radford等、Wang等、Yu等研究泡沫彈丸沖擊三明治梁的過程中,用準靜態應力- 應變曲線和Crushable Foam本構模型描述泡沫鋁的動態沖擊行為。因此,在仿真分析時忽略了泡沫鋁的應變率效應,并使用Crushable Foam本構模擬泡沫鋁。

圖14 泡沫彈丸沖擊某裝甲鋼均質梁有限元模型Fig.14 Finite element model of the monolithic armor steel beam subjected to foam projectile impact

為得到有限元模擬的最優網格尺寸,首先進行網格無關性分析。采用C3D8R實體網格劃分均質梁和泡沫鋁彈丸,均質梁在厚度方向分為3層網格,每層厚1.71 mm,在面內方向,網格尺寸大小依次取2 mm、3 mm、4 mm、5 mm、6 mm、7 mm、8 mm;對泡沫鋁彈丸的網格做相應劃分。圖15給出了泡沫鋁彈丸沖量=8.0 kPa·s時均質梁無量綱中點峰值撓度和無量綱殘余撓度隨網格尺寸的變化趨勢,其中,無量綱中點峰值撓度等于中點峰值撓度除以梁長度的一半,無量綱殘余撓度等于殘余撓度除以梁長度的一半;其他沖量下的變化趨勢類似,故不重復給出。當網格尺寸為2 mm時,無量綱中點峰值撓度和無量綱殘余撓度趨于穩定,并且與試驗值相接近。因此,選擇尺寸為2 mm的網格進行后續分析。此外,為比較不同本構模型對均質梁動態力學響應的影響,分別采用J-C本構模型和理想彈塑性本構模型對某裝甲鋼進行了仿真模擬。

圖15 均質梁無量綱中點峰值撓度和無量綱殘余撓度隨網格尺寸變化趨勢(I=8.0 kPa·s)Fig.15 Midpoint peak displacement and residual deflection of the monolithic beam versus mesh size (I=8.0 kPa·s)

4.3 結果與討論

表3給出了3次重復實驗條件下泡沫鋁彈丸沖擊均質梁的試驗與仿真結果對比。由表3可見:中點峰值撓度和殘余撓度的試驗和仿真結果吻合較好,相對誤差全部小于10%;隨著彈丸沖量的增加,峰值撓度和殘余撓度都有所增加,但試驗后的泡沫鋁名義應變基本不變。表3結果表明,獲得的某裝甲鋼J-C本構參數是可靠的。

表3 泡沫鋁彈丸沖擊均質梁的試驗與仿真結果對比

圖16給出了沖量為8.0 kPa·s的彈丸沖擊作用下均質梁的變形過程。由圖16可見:彈丸高速撞擊均質梁的前2 ms,泡沫鋁劇烈壓縮,均質梁獲得速度并與彈丸一起向右運動;=3.9 ms時,彈丸與均質梁中點運動到最大撓度處,二者此時的中點速度均為0 m/s,均質梁中存儲的彈性能最大;>3.9 ms期間,均質梁中存儲的彈性能釋放后轉化為均質梁和彈丸的動能,二者獲得向左的速度,開始反向運動(即回彈)。圖17給出了試驗結束后泡沫鋁彈丸的最終形貌。表3對比了泡沫鋁彈丸名義應變的試驗和仿真結果,相對誤差在10%左右。根據圖16的高速攝影照片,可對均質梁的動態響應過程做出定性分析。下文結合有限元計算結果進一步分析。

圖16 泡沫鋁彈丸沖量I=8.0 kPa·s時均質梁動態變形過程Fig.16 Dynamic deformation of the monolithic beam subjected to an impact impulse of I=8.0 kPa·s

圖17 試驗后的泡沫鋁彈丸最終形貌Fig.17 Final profiles of the aluminum foam projectiles after tests

沖量為8.0 kPa·s條件下,圖18對比了無量綱中點撓度試驗和仿真結果隨時間變化的曲線,圖19給出了相應的無量綱撓度隨時間變化曲線。圖18表明無量綱中點撓度的試驗結果和基于J-C本構模型的仿真結果總體吻合較好:前1.5 ms時,試驗和仿真曲線完全重合;>1.5 ms時,有限元結果略低于試驗結果,但變化趨勢相同;=3.9 ms時,試驗和仿真結果同時到達峰值,隨后出現下降(即均質梁反向運動,見圖16)?;贘-C本構的仿真結果與試驗結果的相對誤差為1.7%;相比較而言,基于理想彈塑性本構的計算結果與試驗值得偏離較大,相對誤差為12.2%。圖19表明,沖擊載荷下,均質梁的無量綱撓度先增后降,再升高、降低,呈現周期振動。取振動時的平衡位置作為殘余撓度,則基于 J-C本構模型的仿真結果與試驗結果的相對誤差為7.6%,而基于理想彈塑性本構模型的仿真與試驗結果的相對誤差為46.4%。相對于J-C本構模型,理想彈塑性本構模型無強化階段,故采用該本構模型進行仿真的誤差較大。

圖18 均質梁無量綱中點峰值撓度隨時間變化曲線(I=8.0 kPa·s)Fig.18 Midpoint displacement of the monolithic beam versus time (I=8.0 kPa·s)

圖19 均質梁無量綱撓度隨時間變化曲線(I=8.0 kPa·s)Fig.19 Deflection of the monolithic beam versus time (I=8.0 kPa·s)

基于J-C本構模型進一步開展仿真計算,分析泡沫鋁彈丸和均質梁系統中的能量轉化。在泡沫鋁彈丸沖量=8.0 kPa·s沖擊載荷下,圖20給出了系統中的總能量、總應變能、總動能、黏性耗散能和摩擦耗散能;圖21所示為不同能量占總能量的百分比隨時間變化曲線;圖22給出了彈丸和均質梁的動能和塑性耗散能隨時間變化曲線;圖23所示為彈丸和均質梁的動能和塑性耗散能占總能量的百分比隨時間變化曲線;圖24展示了均質梁在泡沫鋁彈丸沖擊下的變形過程。由圖20~圖24可見:=0 ms時刻,系統和彈丸的總動能最大,其他能量為0 J;0<<0.4 ms,彈丸的動能急劇下降,均質梁的動能有所上升(見圖22),但系統的總動能急劇下降,主要轉化為彈丸和均質梁的應變能、彈丸的黏性耗散能(見圖20),變形模式上表現為泡沫鋁的急劇壓縮(見圖24(a)~圖24(c));=0.4 ms時刻,泡沫鋁壓縮完畢,黏性耗散能不再變化;0.7 ms<<3.9 ms,均質梁的動能開始下降,轉化為均質梁的彈性能和塑性耗散能,均質梁出現塑性變形(見圖22);=3.9 ms時刻,均質梁的動能最低,其中點撓度達到最大(見圖24(g));>3.9 ms,均質梁中存儲的彈性能轉化為彈丸和均質梁的動能,二者開始反向運動(見圖24(g)~圖24(i))。

圖20 系統中總應變能、總動能、黏性耗散能和摩擦耗散能隨時間變化的曲線(I=8.0 kPa·s)Fig.20 Total strain energy, kinetic energy, viscous dissipation energy and frictional dissipation energy in the system versus time (I=8.0 kPa·s)

圖21 不同能量占總能量的百分比隨時間變化曲線 (I=8.0 kPa·s)Fig.21 Proportions of various energy versus time (I=8.0 kPa·s)

圖22 彈丸和均質梁的動能和塑性耗散能隨時間變化曲線 (I=8.0 kPa·s)Fig.22 Kinetic energy and plastic dissipation energy of foam projectile as well as monolithic beam versus time (I=8.0 kPa·s)

圖23 彈丸和均質梁的動能和塑性耗散能占總能量的百分比隨時間變化曲線(I=8.0 kPa·s)Fig.23 Percentage of kinetic energy and plastic dissipation energy of foam projectile as well as monolithic beam versus time (I=8.0 kPa·s)

圖24 泡沫鋁彈丸沖擊下的均質梁變形過程(I=8.0 kPa·s)Fig.24 Numerically predicted deformation process of the monolithic beam subjected to foam projectile impact (I=8.0 kPa·s)

圖21表明:0 ms<<3.9 ms期間,系統的動能在系統總能量中的占比逐漸減低,應變能和黏性耗散能的占比則升高;>3.9 ms期間,應變能中的彈性能和動能之間存在相互轉化,梁和彈丸獲得向左運動的速度;=5 ms時,系統總應變能占比為69.3%,黏性耗散能占比29.1%,總動能占比1.1%,摩擦耗散能的占比僅為0.7%。圖23表明,=5 ms時,泡沫鋁彈丸的塑性耗散能占比為41.6%,加上其黏性耗散能占比28.9%,彈丸的總耗散能為70.5%,均質梁的塑性耗散能占比則為18.2%。

5 結論

為準確模擬具有高屈服強度的某裝甲鋼的動態力學行為,本文基于J-C本構模型通過常溫和高溫環境下的準靜態拉伸實驗數據,以及常溫環境下的動態壓縮實驗數據,擬合得到了J-C本構參數;隨后采用輕氣炮和泡沫鋁彈丸對該裝甲鋼均質梁開展沖擊試驗研究,同時分別采用J-C本構模型和理想彈塑性本構模型進行有限元仿真計算,并對比分析了試驗與數值仿真結果。結果表明:某裝甲鋼具有應變率強化效應,溫度軟化效應顯著;采用J-C本構模型仿真的均質梁峰值撓度與試驗結果的相對誤差為1.7%~6.1%,殘余撓度相對誤差為0.6%~7.6%。

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 国产95在线 | 日韩性网站| 日本91视频| 一区二区影院| 最新加勒比隔壁人妻| 曰AV在线无码| 中文无码精品a∨在线观看| 亚卅精品无码久久毛片乌克兰| 久久精品亚洲中文字幕乱码| 99久久婷婷国产综合精| 国产农村1级毛片| 中文字幕在线观| 久一在线视频| 国产女人综合久久精品视| 播五月综合| 欧美全免费aaaaaa特黄在线| 国内精品伊人久久久久7777人| 亚洲精品动漫| 国产精品久久久久久搜索| 国产91导航| 国产欧美视频综合二区| 国产人在线成免费视频| 欧美成人a∨视频免费观看| av无码久久精品| 69免费在线视频| 2021国产精品自拍| 成年免费在线观看| 国产情精品嫩草影院88av| 在线综合亚洲欧美网站| 国产人妖视频一区在线观看| 久久久久久午夜精品| 国产在线自在拍91精品黑人| 亚洲最大综合网| 中文字幕乱码二三区免费| 中文字幕欧美日韩| 四虎免费视频网站| 国产成人狂喷潮在线观看2345| 免费观看男人免费桶女人视频| 国产国拍精品视频免费看 | 40岁成熟女人牲交片免费| 久久情精品国产品免费| 亚洲视频影院| 久久综合色视频| 在线国产综合一区二区三区| 亚洲日韩国产精品综合在线观看| 蜜桃臀无码内射一区二区三区| 精品一区二区三区视频免费观看| 在线免费看片a| 久久国产精品娇妻素人| 国产网站免费| 亚洲国产中文综合专区在| 色婷婷成人| swag国产精品| 免费毛片全部不收费的| 美女被操黄色视频网站| 国产中文一区二区苍井空| 999精品色在线观看| 欧美第九页| 国产精品美女在线| 久久青草热| 日韩欧美网址| 欧美人与性动交a欧美精品| 成人福利在线观看| 免费福利视频网站| 狠狠色噜噜狠狠狠狠色综合久| 自慰高潮喷白浆在线观看| 99r在线精品视频在线播放| 国产成本人片免费a∨短片| 国产一区二区三区免费| 丁香五月婷婷激情基地| 日韩亚洲综合在线| 亚洲六月丁香六月婷婷蜜芽| 四虎影视无码永久免费观看| 97se亚洲综合在线| 这里只有精品在线| 国产成人免费观看在线视频| 久久特级毛片| 精品国产成人三级在线观看| 成年女人a毛片免费视频| 日本人妻一区二区三区不卡影院| 99青青青精品视频在线| 九色综合伊人久久富二代|