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

X80鋼裂紋擴展與失效過程的多尺度模擬與實驗驗證

2022-01-15 04:07:40王汝恒錢家鈺鄧淞文
西南科技大學學報 2021年4期
關鍵詞:裂紋模型

盛 鷹 賈 彬 王汝恒 錢家鈺 鄧淞文

(1. 西南科技大學土木工程與建筑學院 四川綿陽 621010;2. 西南科技大學工程材料與結構沖擊振動四川省重點實驗室 四川綿陽 621010)

為滿足長輸油氣管道向大口徑、高壓力發展的要求,針對高鋼級管線鋼的組織形態、力學性能、變形破壞機理等開展研究具有重要意義[1-3]。X80管線鋼是當前應用最廣泛的高強度鋼之一,其在復雜的受力狀況和工作環境中可能面臨疲勞裂紋擴展、腐蝕損傷、局部屈曲失穩等多種破壞或失效風險[4-5]。

不少學者開展了X80鋼在不同工況下宏觀力學性能的研究。徐杰等[6]研究了不同的焊接溫度對X80鋼的斷裂韌度的影響。劉英義等[7]研究了X80鋼在土壤環境中的腐蝕行為。安騰[8]研究了高壓氫對X80鋼疲勞裂紋損傷行為的影響。然而,大部分學者重點研究X80鋼的宏觀尺度力學性能,而沒有特別關注材料的微觀缺陷的擴展對其宏觀力學性能的影響。

為研究材料在原子尺度的微缺陷對連續介質尺度的宏觀力學性能的影響,需要建立合理的多尺度耦合模型[9]。多尺度耦合計算方法的難點在于兩個方面[10]:一是如何實現各物理量在不同尺度下的傳遞和過渡;二是如何實現邊界條件與計算結果在不同尺度下的耦合。為解決這兩個難題,不少學者在這個領域開展了大量的研究工作,建立了多種在原子尺度與連續介質尺度之間進行數據傳遞的方法,如:遞進式原子-連續關聯方法[11]、并發式原子-連續關聯方法[12]、多尺度內聚力模型[13]等等。其中,多尺度內聚力模型的基本思想是:首先采用分子動力學方法獲得材料的界面層黏結力與張開位移的關系(T-S曲線),再將該T-S曲線的關鍵信息輸入到有限元模型中用于進一步研究材料在宏觀尺度上的力學行為。

為研究帶微缺陷的X80鋼的裂紋擴展與失效演化機制,本文在原子尺度建立了多組包含不同微觀缺陷的X80鋼拉伸微觀計算模型,獲得了各條包含裂紋尖端演化微觀信息的T-S曲線,根據T-S曲線上關鍵轉折點對應的原子構型圖,從原子尺度觀察分析了X80鋼中不同微觀缺陷擴展的現象和規律,并以鈍裂紋模型為例,獲得了其微裂紋的擴展軌跡;然后,在連續介質尺度建立了包含內聚力單元的單向拉伸CT試件有限元模型, 將T-S曲線中反映出的材料核心微觀特征引入內聚力模型中作多尺度分析,研究了微觀缺陷對材料宏觀斷裂參數的影響程度。通過實驗驗證了多尺度數值模擬結果的可靠性。

1 基于內聚力模型的多尺度方法

1.1 分子動力學模型及勢函數驗證

X80鋼的金相組織分析結果顯示,X80鋼主要由針狀鐵素體和少量珠光體構成:鐵素體為體心立方(BCC)晶格結構,其強度、硬度不高,但具有良好的塑性與韌性;珠光體是鐵素體與滲碳體的共析體,其力學性能介于鐵素體與滲碳體之間,強度較高,硬度適中,塑性和韌性較好;滲碳體是X80鋼中主要的強化相,是鐵碳合金按亞穩定平衡系統凝固和冷卻轉變時析出的Fe3C,具有復雜斜方晶體結構,其硬度很高,但塑性和沖擊韌性較差,易發生脆斷。

建立α-Fe/Fe3C的晶體結構,其中,α-Fe為BCC結構,晶格常數a0=b0=c0=0.28664 nm,α=β=γ=90°;Fe3C為斜方晶體結構,晶格常數a0=0.506 nm,b0=0.674 nm,c0=0.451 nm,α=β=γ=90°,其每個單胞由12個Fe原子和4個C原子構成,如圖1所示。

圖1 Fe3C單胞結構Fig.1 Fe3C crystal structure

建立α-Fe/Fe3C分子動力學計算模型,模型尺寸為28.664 nm×14.332 nm×1.146 nm,共計39 788個原子。選擇嵌入式原子勢函數EAM/FS[14],其第i個原子的勢能如式(1)所示:

(1)

式中Ak,Rk,ak,rk為常數,可由第一性原理計算得到[15]。

勢函數通常是經驗性的,在分子動力學仿真中往往存在爭議。為驗證α-Fe/Fe3C的EAM/FS勢函數及其參數,采用LAMMPS獲得了α-Fe/Fe3C的熔點。

設置初始構型的時間步長為0.01 ps,采用三維周期性邊界條件,在2.5 K的溫度下弛豫500 000步。利用Nose-Hover方法將壓力保持在0附近,并將系統溫度由T=2.5 K逐漸升高至2 500 K。在仿真過程中,每1 000步輸出一次熱動力學結果。升溫過程中原子平均體積隨模擬溫度的變化如圖2所示。

圖2 升溫過程中原子平均體積隨模擬溫度的變化Fig.2 Variation of average atomic volume with simulated temperature during heating process

升溫過程中,平均原子體積幾乎隨溫度線性增加。然而,當溫度增加到某一特定值時,平均原子體積突然增加。這個現象說明,系統正經歷著某種相變。α-Fe/Fe3C的熔點約為 1 873 K,仿真獲得的相變時的溫度值為 1 990 K。兩者之間的誤差約為 6.2%。驗證了勢函數的合理性。

1.2 T-S曲線與內聚力模型

考慮到裂紋尖端的應力不可能達到無窮大,學者們發展了內聚力模型。材料的開裂、裂紋擴展和破壞是由界面單元決定的,但是在內聚力模型中,裂紋擴展區域被簡化為一個零厚度區域,該區域包含兩個相重疊的面,其中一個面對任意材料都視為無損傷連續介質,另一個面為連續介質單元之間的黏聚界面,用于判定材料的損傷。在加載過程中,兩個面彼此分離,但內聚力單元并不會被破壞。當材料損壞發生時,黏聚界面張開并喪失剛度,這就意味著該黏聚界面連接的兩端連續介質不再相連,即材料在該處斷裂。基于該原因,在內聚力模型中,裂紋只能沿黏聚界面的單元傳播。

在使用黏聚單元模擬裂紋擴展之前,必須首先獲得黏聚界面的黏結力T與裂紋張開位移S之間的關系(簡稱T-S關系)。黏聚界面的分離程度是由位移間斷值δ來衡量的,即等于相鄰連續單元的位移u+和u-的差值,如式(2)所示:

δ=u+-u-

(2)

黏聚界面的分離量δ取決于分別作用于黏聚界面表面的正應力和剪應力值,而裂紋張開位移S即為各處黏聚界面的分離量δ之總和。判斷黏聚界面是否分開的判據有兩個:內聚力單元的臨界分離位移δ0和臨界分離應力T0。

更具體地說,在分子動力學模型中,取裂紋面上、下一定范圍的原子區域為T-S區域,分別計算這兩個區域內原子的平均位移,將兩者之差作為裂紋的張開位移S;將T-S區域內原子應力的平均作為材料的界面黏結力T。在計算過程中,只采用了T-S區域內的原子來獲得界面黏結力與裂紋張開位移,這樣能有效地避免數據的分散性。將不同時刻得到的S,T的值繪制于同一坐標系中即可得到T-S曲線。

在內聚力模型中,它們之間的黏結力和張開位移的關系服從事先指定的T-S曲線。而T-S曲線的形狀、極值、分布等信息充分反映了與裂紋尖端演化相關的微觀信息,如位錯發射、裂尖鈍化、孿晶、相變等等。運用內聚力單元模擬裂紋擴展過程必須確定3個特征量:內聚力單元的初始剛度、失效閾值、失效演化參量。而這3個特征量均可從T-S曲線中獲得。

綜上分析可得,本文建立的基于內聚力模型分析復合微觀缺陷擴展與失效的多尺度方法的基本步驟如下:

(1)建立相應的帶裂紋的分子動力學模型,確定T-S區域;

(2)統計T-S區域內的原子平均應力與裂紋張開位移,并繪制T-S曲線;

(3)根據T-S曲線獲得內聚力單元材料參數;

(4)建立標準緊湊拉伸試件和三點彎曲試件的有限元模型,并將通過T-S曲線獲得的參數賦予內聚力單元;

(5)在內聚力模型中模擬裂紋擴展,獲得相應的材料宏觀力學性能參數。

需要特別說明的是,在分子動力學模型中,T-S區域的尺寸選擇對結果有較大影響。在裂紋的有效影響區域內,T-S曲線應該由上升段和下降段構成。其中,在T-S曲線的上升階段,該區域的變形處于彈性階段,基本沒有金屬鍵的斷裂;而在T-S曲線的下降階段,裂紋尖端會出現明顯鈍化和材料損傷。受裂紋影響不大的區域,其T-S曲線只有應力上升段而沒有下降段,該部分區域的變形始終處于彈性階段,不應納入T-S曲線的綜合統計范圍,否則計算得到的界面黏結力的值將比真實值偏小,是不合理的。研究發現,裂紋的影響區域為裂紋長度約6倍的范圍內,即合理的T-S區域應取裂紋長度約6倍的范圍。

2 計算結果及分析

2.1 含不同微缺陷力學模型的T-S曲線

油氣管道在冶金制造過程中會產生微裂紋和微空洞等缺陷,在服役過程中,也會因為輸送介質具有很強的腐蝕性而造成管道內表面的腐蝕開裂,從而形成微觀缺陷。這些缺陷在管道內的高壓作用下可能逐漸擴展形成宏觀裂紋,存在使管道破裂的風險。因此,建立如圖3所示的分子動力學微觀尺度計算模型:模型Ⅰ不含微缺陷;模型Ⅱ為中心鈍裂紋(裂紋長2.5 nm,寬0.52 nm);模型Ⅲ是為中心空洞(空洞直徑0.6 nm)。

圖3 微觀尺度計算模型Fig.3 Microscale computation models

設置[100]和[010]方向為自由邊界條件(s),在[001]方向采用周期性邊界條件(p)。分子模擬的初始溫度為300 K,加載方式為沿[100]方向施加速率為0.025 nm/ps的拉伸應變。模擬中取時間步長t=1 ps,求解牛頓運動方程的積分方法選用Velocity-Verlet算法。模擬過程分兩步:在微正則系綜(NVE)下,首先對初始構型弛豫使其達到位能最低狀態,然后對弛豫后的納米單晶體沿[100]方向均勻拉伸20 000步。加載過程中利用Nose-Hoover方法調節溫度使其保持恒定。模擬過程中每200步記錄原子的位置、應力、總能量,并輸出系統的幾何構型。取距裂紋上、下表面1.0 nm范圍內的區域為T-S數據區。

圖4、圖5、圖6分別為通過分子動力學模擬獲得的模型Ⅰ、模型Ⅱ、模型Ⅲ的T-S曲線。在裂紋張開的任一時刻,T-S曲線中都能得到一個離散點并對應一種原子構型。因此,裂紋尖端的復雜變形機制能綜合反映到T-S曲線中。為了觀察微觀變化(如位錯、裂尖鈍化、二次裂紋萌生與擴展等等)對T-S曲線的影響,在T-S曲線中標出了關鍵轉折點所對應的原子構型圖。

圖4 模型Ⅰ在拉伸載荷下的T-S關系與變形機制Fig.4 T-S relationship and deformation mechanism of model Ⅰ under tensile load

圖5 模型Ⅱ在拉伸載荷下的T-S關系與變形機制Fig.5 T-S relationship and deformation mechanism of model Ⅱ under tensile load

圖6 模型Ⅲ在拉伸載荷下的T-S關系與變形機制Fig.6 T-S relationship and deformation mechanism of model Ⅲ under tensile load

為直觀比較模型Ⅰ、模型Ⅱ、模型Ⅲ的T-S曲線的差異,將3組T-S曲線放置在一張圖中,如圖7所示。

圖7 模型Ⅰ、模型Ⅱ、模型Ⅲ的T-S曲線對比Fig.7 T-S curve comparison of model Ⅰ, Ⅱ and Ⅲ

由圖4-圖7可得如下結果:

(1)圖4顯示了無缺陷模型Ⅰ在拉伸荷載作用下原子鍵斷裂、微裂紋萌生和擴展、模型斷裂的過程,其T-S曲線走勢與宏觀尺度試件拉伸的應力-應變曲線走勢相似。T-S曲線大致分為4個階段:第一階段的T-S曲線近似呈線性關系,最高點①處的原子應力為12.75 GPa;第二階段隨著裂紋張開位移的不斷增大,原子應力增長緩慢,由①和②的原子構型對比可見,處于模型中線附近的原子排列方式由穿插排列變為并行排列,裂紋附近原子呈無序排列狀態;第三階段原子應力隨裂紋張開位移的增大而繼續增大,局部出現裂紋,當應力峰值達到19.89 GPa時,裂紋完全貫通;第四階段原子應力急劇下降,模型完全斷裂。

(2)圖5顯示了模型Ⅱ的裂尖變形過程,其裂尖變形及失效機制為鈍化和位錯發射。模型Ⅱ的T-S曲線出現了兩次峰值。當T-S曲線第一次達到峰值點時,裂紋尖端應力集中嚴重,裂紋兩端與裂紋長度方向呈45° 的方位出現位錯,原子出現無序排列狀態,隨后應力小幅下降;隨著加載的繼續,裂紋鈍化程度進一步加大,位錯區域不斷擴大,裂紋前方不遠處出現二次裂紋,在主裂紋與二次裂紋連通前,應力有一次小幅上升,當主裂紋與二次裂紋間的大量連接原子都呈無序排列狀態時,應力達到第二個峰值9.52 GPa;隨著進一步加載,應力迅速下降,直至裂紋完全連通。

(3)圖6顯示了模型 Ⅲ 的裂尖變形過程,其裂尖變形及失效機制仍為鈍化和位錯發射。模型Ⅲ的T-S曲線出現了3次峰值。當T-S曲線第一次達到峰值點10.08 GPa時,空洞形狀由圓形拉伸為橢圓形,空洞附近的原子間距增大;隨著加載的繼續,空洞兩端沿垂直于加載的方向出現位錯,位錯前方原子呈無序排列狀態,當位錯導致的呈無序排列狀態的原子大量增多時,空洞裂紋擴展為條狀并逐漸鈍化,出現應力松弛,應力先小幅下降再小幅上升,并出現第二次峰值;繼續加載,裂紋繼續擴展,擴展路線呈非線性,應力小幅下降后上升至10.56 GPa,此時模型完全斷裂。

由圖4-圖7可進一步發現:

(1)由圖5與圖4對比可知,當模型完全斷裂時,模型Ⅱ裂紋張開位移約為模型Ⅰ裂紋張開位移的50%,模型Ⅱ最大原子應力也約為模型Ⅰ最大原子應力的50%。表明在微觀尺度,鈍裂紋對模型的強度具有不可忽略的影響。

(2)由圖5與圖6對比可知,當模型完全斷裂時,模型 Ⅲ 裂紋張開位移比模型Ⅱ裂紋張開位移約大20%,模型 Ⅲ 最大原子應力也比模型Ⅱ最大原子應力略大。表明在微觀尺度,空洞對模型的強度仍具有不可忽略的影響,但空洞對模型強度的影響比鈍裂紋對模型強度的影響更小。

(3)由圖7可知,帶微缺陷的X80鋼的微觀模型的破壞情況存在如下共同點:T-S曲線中的最高峰值點表明材料中開始出現損傷(裂紋起裂),裂尖開始出現明顯的位錯發射;隨著加載的繼續,應力集中導致原子間金屬鍵被破壞,當應力集中區域附近的大量原子呈無序排列狀態時,裂尖明顯鈍化,應力出現松弛,應力小幅下降。

(4)在微觀尺度,微缺陷會降低模型強度,不同形狀的微缺陷對模型強度的影響不同,在外荷載作用下不同形狀缺陷擴展的方向也不同。因此,研究不同微缺陷對微觀模型強度的影響是有意義的。

(5)T-S曲線綜合反映了裂紋尖端復雜的變形機制(如位錯發射、裂尖鈍化等等),因此,將分子動力學得到的反映微觀特性的T-S曲線作為搭建多尺度計算的橋梁是可行的。

2.2 微裂紋擴展軌跡與實驗驗證

由上述分析可知,鈍裂紋對微觀模型強度的影響比空洞更大,對X80鋼而言,鈍裂紋更容易導致材料沿滑移面破壞。故以鈍裂紋模型Ⅱ為算例,分析其微裂紋的擴展軌跡。

裂紋尖端原子的擴散常數可以通過Einstein布朗動力學模型得到:

(3)

式中:N為系統原子總數;S為模型維數,當前系統采用三維模型,故S=3;t為時間。

通過原子位置及最大應力,即可確定裂尖原子位置及位移,進而計算裂尖原子的擴散系數D0。

微觀模型下,裂紋的擴展路徑由位移函數(4)控制:

(4)

式中:(X(t),Y(t)),(X(t+Δt),Y(t+Δt))分別表示當前時刻、下一時刻的裂尖坐標;Vm為裂紋擴展臨界速度,本文取Vm=605.74 m/s,Δt=0.01 μs;G(t)是由C++程序生成的一個標準正態隨機變量。

于是,由式(4)可迭代得到裂尖原子下一時刻的位置。圖8為微觀尺度下裂紋尖端的軌跡。

圖8 微觀尺度下裂紋尖端的軌跡Fig.8 Trace of crack tip at microscale

由圖8可知,裂尖原子在微觀尺度下是劇烈震蕩的,而在宏觀尺度作有限元計算得到的結果近似為一條直線。由此可以看出,從微觀尺度研究裂紋的擴展是有必要的。

下面采用X80管件母材三點彎曲(SEB)試樣的斷裂韌性實驗驗證裂紋擴展軌跡的非線性[16]。圖9所示為試樣CTOD試驗結果的δ-Δa試驗關系及斷口照片。

觀察圖9可知,X80的沖擊斷口呈非鏡像對稱,斷口沿裂紋延長線方向的擴展呈鋸齒狀,并非為嚴格的直線。這與分子動力學計算得到“裂紋尖端的軌跡在微觀尺度下是劇烈震蕩的”結論一致,這也間接定性驗證了本文分子動力學計算結果的合理性。

圖9 CTOD試樣的δ-Δa試驗關系及斷口照片Fig.9 δ-Δa experiment curve of Sample CTOD and the photographs of fracture surface

2.3 T-S關系對宏觀斷裂的多尺度模擬

為分析原子尺度下微結構的演化對較高尺度下材料失效行為的影響,本節建立了包含內聚力單元的CT試件模型,并模擬了CT試件在位移加載條件下的裂紋擴展過程。通過對CT試件裂紋擴展的模擬,能檢驗分子動力學獲得的黏聚單元法則(T-S曲線)能否滿足線彈性斷裂力學的理論預期結果。CT試件的幾何尺寸和網格模型如圖10所示。CT試件的寬和高分別為1.250W和1.200W,其中W為銷孔中心到不含裂紋邊界的距離,取W=200 μm。初始裂紋長度為a=0.500W,銷孔半徑為 0.250W。在預定義裂紋擴展路徑上設置寬為0.1 nm的內聚力單元。在參考點RP-1和RP-2上施加位移載荷。

圖10 CT試件的幾何尺寸Fig.10 Geometric dimension of CT specimen

接下來將分子動力學模型T-S曲線的特征值輸入ABAQUS軟件中模擬內聚力模型的多尺度行為。運用ABAQUS軟件模擬內聚力模型的多尺度行為時,采用雙線性模型擬合T-S曲線,如圖11所示。在圖11中,考慮材料為各向同性,即初始剛度Kn,Ks,Kt均相等;雙線性模型所圍成的三角形面積為能量釋放率GTC。

圖11 雙線性模型示意圖Fig.11 Schematic diagram of bilinear model

對I型斷裂模型,設剛度Kn=Ks=Kt=K(K為輸入ABAQUS中的剛度),雙線性模型所圍成的三角形面積也等于斷裂能GC,即GTC=GC。則剛度K的表達式為:

(5)

(6)

通過在兩個銷孔處的位移加載,可觀察到裂紋張開的行為。當裂紋擴展發生時,內聚力區域即開始生成。一旦達到了臨界位移,裂紋便開始擴展。同樣以鈍裂紋模型Ⅱ為算例,給出不同時刻CT試件的裂紋擴展演示圖,如圖12所示。

圖12 與模型Ⅱ對應的CT試件裂紋擴展應力圖Fig.12 Crack propagation stress diagram of CT specimen corresponding to model II

CT試樣中裂紋張開位移與約束力的關系曲線如圖13所示。

圖13 CT試樣中位移與約束力的關系曲線Fig.13 Relation curve between displacement and binding force in CT sample

采用虛擬裂紋閉合法(VCCT)可計算得到宏觀應力強度因子,其解析式為:

(7)

從加載開始至裂紋起裂時的應力強度因子如圖14所示。由圖14可知,裂紋起裂時的應力強度因子為15 624.55 MPa·μm1/2。

圖14 加載開始至裂紋起裂時的應力強度因子Fig.14 Stress intensity factor from the beginning of loading to crack initiation

3 結論

(1)在原子尺度計算獲得的T-S曲線包含了與裂紋尖端演化相關的微觀信息,綜合反映了材料的各種微觀缺陷的特性。不同材料、不同微觀缺陷擴展得到的T-S曲線是不同的,故將T-S曲線融入內聚力模型中以描述不同微觀缺陷對材料宏觀特性的影響是可行的。

(2)微觀缺陷對微觀尺度模型強度的影響是不可忽略的。對X80鋼而言,鈍裂紋、空洞等微缺陷對微觀模型裂紋尖端最大應力、微觀模型強度的影響程度不同,鈍裂紋的影響大于空洞的影響。

(3)本文提供的基于內聚力模型的多尺度方法可用于分析微觀缺陷對宏觀應力強度因子、斷裂韌性等斷裂參數的影響,還能用于分析更復雜的斷裂機制問題。

猜你喜歡
裂紋模型
一半模型
裂紋長度對焊接接頭裂紋擴展驅動力的影響
一種基于微帶天線的金屬表面裂紋的檢測
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
Epidermal growth factor receptor rs17337023 polymorphism in hypertensive gestational diabetic women: A pilot study
微裂紋區對主裂紋擴展的影響
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
預裂紋混凝土拉壓疲勞荷載下裂紋擴展速率
主站蜘蛛池模板: 久久福利片| 69国产精品视频免费| 国产精品太粉嫩高中在线观看| 久久黄色免费电影| 日本久久网站| 欧美一区二区三区不卡免费| 人妻精品久久无码区| 日本高清在线看免费观看| a级毛片一区二区免费视频| 乱人伦视频中文字幕在线| 日韩欧美国产三级| 一级高清毛片免费a级高清毛片| 九色视频线上播放| 91免费国产在线观看尤物| 无码中文字幕加勒比高清| 国产在线一区视频| 四虎国产精品永久在线网址| 在线播放真实国产乱子伦| 国产日韩丝袜一二三区| 国产精品福利导航| AV在线天堂进入| 天堂网亚洲系列亚洲系列| 无码国产伊人| 国产成人你懂的在线观看| 国产黄色免费看| 色哟哟国产精品| 色偷偷一区二区三区| 中文字幕2区| 亚洲an第二区国产精品| 亚洲欧美不卡中文字幕| 久久精品视频亚洲| 77777亚洲午夜久久多人| 粉嫩国产白浆在线观看| 天天躁狠狠躁| 玖玖免费视频在线观看| 国产电话自拍伊人| 欧美午夜在线观看| 久久国产高潮流白浆免费观看| 一级黄色网站在线免费看| 手机在线国产精品| 亚洲高清在线天堂精品| 91视频区| 中国毛片网| 国产91成人| 天天视频在线91频| 九九热在线视频| 91免费国产在线观看尤物| 免费xxxxx在线观看网站| 国产在线麻豆波多野结衣| 在线亚洲精品自拍| 国产成人精品高清不卡在线 | 国产精品永久不卡免费视频| 国产精品手机在线观看你懂的| 欧美黑人欧美精品刺激| 精品国产aⅴ一区二区三区| 欧洲熟妇精品视频| 中文字幕欧美日韩高清| 国产裸舞福利在线视频合集| 综合网久久| 国产成人精品视频一区二区电影| 91精品国产91久无码网站| 久久综合九色综合97网| 久久先锋资源| 亚洲无码高清一区| 国产成人高清精品免费软件| 免费看av在线网站网址| 国产全黄a一级毛片| 一级毛片在线播放免费| 国产精品v欧美| 天天综合天天综合| 国产麻豆va精品视频| 大陆精大陆国产国语精品1024| 国产第一页屁屁影院| av天堂最新版在线| 精品夜恋影院亚洲欧洲| 看国产毛片| 国产成人精品亚洲日本对白优播| 曰AV在线无码| 欧美精品亚洲日韩a| 一级一级一片免费| 国产真实自在自线免费精品| 一本综合久久|