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

高強螺栓抗剪連接滑移數值模型

2021-01-08 08:53:34趙中偉樊雄濤
同濟大學學報(自然科學版) 2020年12期
關鍵詞:模型

趙中偉,樊雄濤,吳 剛

(1.東南大學土木工程學院,江蘇南京211189;2.遼寧工程技術大學土木工程學院,遼寧阜新123000)

有限元法已成為幾乎所有研究領域的主要分析方法,例如機械和土木工程。數值分析方法可以用來確定各種連接的潛在失效模式和極限承載能力。高強度螺栓廣泛應用于建筑結構中的梁柱連接。由于摩擦的存在,高強度螺栓抗剪連接的加載過程是一種高度非線性的力學行為,在數值分析中需要多次反復迭代,計算時間長而且很容易產生不收斂的結果。

高強度螺栓對整體節點的力學性能有顯著影響。為了精確模擬高強螺栓的力學性能,已有很多學者在節點有限元模型中建立了精細化的高強螺栓的數值模型[1-3]。馬舒淇[4]等人利用三線性模型將錨桿與巖石界面的滑移模型進行了簡化。Hwang[5]利用精細化的三維數值模型對螺栓的安裝過程進行了動態仿真,通過施加扭矩使螺栓逐步擰緊直至破壞。精細化的數值模型需要建立在足夠的計算能力上才能完成預定目標。以目前的計算機的計算能力,可以說不可能在整體結構中建立精細化的數值模型進行靜力分析,更別說滯回分析以及動力分析。另外,由于眾多接觸單元的存在,不收斂是目前精細化數值模型不可回避的一個問題。

近年來,已有很多學者對如何減少高強度螺栓精細化數值模型的計算成本進行了廣泛而深入的研究。Liu和Chen[6]通過建立多尺度有限元模型以減小包含螺栓群數值模型對計算機計算能力的要求。Bogdanovich和Kizhakkethara[7]利用子模型的方法對螺栓連接進行了應力分析,結果表明子模型為螺栓連接高梯度應力區的應力分析提供了一種有效的方法。Pearce[8]利用顯式有限元方法對螺栓連接的擬靜力及滯回性能進行了研究。以上所述的改進方法可以在一定程度上減少計算時間,但計算代價大且計算結果發散的問題并沒有得到根本的解決。李啟才[9]對高強度螺栓連接的力學性能進行了試驗研究。石永久和王元清等學者對螺栓的力學性能進行了一系列深入的研究[10-12],提出了高強度螺栓的抗剪簡化模型,及循環荷載作用下的滯回模型。

為了改善目前精細化螺栓連接數值模型收斂困難且計算時間長的問題,本文基于通用有限元軟件提出了簡化的高強度螺栓滑移數值模型。利用ANSYS同時建立了高強度螺栓的精細化數值模型和簡化的螺栓滑移數值模型。在此基礎上,對兩種模型進行了滯回性能分析對比,驗證了所提出的滑移數值模型計算結果的可靠性,通過參數化分析,對螺栓滑移數值模型的計算誤差進行了系統研究。

1 建立數值模型

建立了高強度螺栓連接的精細化數值模型及簡化的螺栓滑移數值模型。螺栓連接的幾何參數見圖1。試樣的長度和寬度分別為515mm和80mm,高強度螺栓的材料為20MnTiB,這種材料在中國GB/T 1231-2006中推薦用于高強度螺栓。螺栓的屈服強度為940MPa,拉伸強度為1 040MPa,螺栓預緊力為155kN,鋼板的材料為Q235B,其屈服強度、彈性模量、泊松比和密度分別為23.5MPa、210GPa、0.3和7 800 kg·m-3,螺栓孔直徑比螺栓桿直徑大1.5mm。

1.1 高強螺栓精細化數值模型

基于所提出的通用有限元程序ANSYS建立了現有的精細化模型,采用SOLID185單元對螺栓和鋼板進行網格劃分,SOLID185單元用于實體結構的三維建模,單元最大尺寸為2mm。該單元具有8個節點,每個節點具有三個自由度,即節點沿x、y和z方向上的平動位移。該單元具有塑性、超彈性、應力強化、蠕變、大撓度和大應變能力。用接觸單元CONTA174和TARGE170對螺栓和不同鋼板之間的接觸行為進行模擬。采用PRETS179對螺栓進行預緊力加載。PRETS179用于定義2D或3D網格結構中的預緊截面,該單元沿規定的加載方向具有一個平移自由度,螺栓連接的精細化模型如圖2所示。接觸和目標表面之間的滲透量取決于法向剛度(FKN),FKN允許的最大值為1,法向剛度越大,計算精度越高,但越不容易收斂。可以通過增加FKN的值來減少穿透,因此,本文法向接觸剛度FKN為1,其他參數采用默認值。選擇罰函數法和拉格朗日乘子作為接觸算法,允許不同鋼板之間的滑動。

圖1 螺栓連接幾何參數(單位:mm)Fig.1 Geometric parameters of the adopted specimen(unit:mm)

圖2 螺栓連接精細化數值模型Fig.2 Refined numerical model of a bolted connection

1.2 簡化高強螺栓抗剪連接滑移數值模型

螺栓連接的不同組件之間的接觸狀態表現出一種高度非線性的行為,單元劃分必須足夠小才能獲得準確的結果。建立現有的螺栓連接精細化數值模型比較繁瑣,計算成本較高。因此,非常有必要提出一種簡化的高強度螺栓連接數值模擬方法,在此基礎上,可以準確地預測螺栓連接的力學性能,并且可以降低計算成本。基于此,可以在大型的梁柱節點數值模型中建立簡化的螺栓滑移數值模型,以精確考慮局部螺栓對節點整體力學性能的影響。

提出一種適用于精確模擬高強度螺栓滑移的數值計算模型。鋼板可采用殼單元,如ANSYS中的SHELL181單元進行模擬;采用梁單元如BEAM188單元進行模擬螺桿;采用只受壓彈簧(COMBIN39)單元模擬螺桿與鋼板孔壁之間的擠壓行為,同時,COMBIN39單元也被用來模擬螺栓與鋼板之間的摩擦行為。COMBIN39單元是具有自定義荷載位移曲線能力的單向單元,該單元具有軸向(longitudinal)和扭轉(torsional)兩個功能選項,可用來對1維、2維和3維結構進行分析。當軸向功能被開啟時,該單元就是具有兩個節點的單軸拉壓單元,每個節點具有沿x、y和z方向上的平動位移。在螺栓孔壁周圍的節點同時建立非線性受壓單元(以下稱為接觸單元)和摩擦單元,即螺孔周圍每個節點位置建立一個受壓彈簧單元和一個非線性摩擦單元。將螺孔壁上所有節點的節點坐標系調整至如圖3a所示。激活MPC184單元的剛性梁選項,用以連接所有非線性彈簧單元的節點(包括接觸單元和摩擦單元)和螺桿。所提出的簡化螺栓滑移數值模型如圖3所示。

為了能夠精確模擬高強度螺栓的摩擦行為,必須對COMBIN39單元進行特殊的設置。摩擦力的大小直接決定于接觸面之間的摩擦系數以及法向壓力的大小。假定摩擦單元能夠承受的最大內力為滑動摩擦力(Fmax),當彈簧單元內力達到該值時,即使位移變化,彈簧內力也不再改變;如果彈簧內力未達到滑動摩擦力(Fmax),此時為靜摩擦力,則彈簧不會發生變形。

為了精確模擬高強度螺栓滑移后螺桿與螺孔的擠壓行為,將非線性彈簧單元的設置如下:

(1)KEYOPT(1)=0:將卸載路徑與加載路徑設置為相同;

(2)KEYOPT(2)=0:將彈簧受壓時的本構設置為預定的荷載-位移曲線;

(3)KEYOPT(3)=0:將單元的自由度設置為沿節點局部坐標系的x軸方向;

將接觸單元在受壓時的剛度和承載力設置為足夠大,以防止螺桿與鋼板之間的侵蝕,經過試算得出兩者分別設置為5×108kN·m-1和105kN已經足夠大,可以滿足要求。受拉時的剛度設置為0。COMBIN39單元的自由度的參考坐標系為節點坐標系,同時接觸單元應該允許螺桿與鋼板之間的相對滑移,滑移大小等于螺桿直徑與螺孔直徑的差。基于此,將螺孔周圍所有節點的坐標系進行旋轉,使節點坐標系的x軸方向沿螺孔的徑向,如圖3所示。與接觸單元所對應的荷載位移曲線如圖4所示。圖中d1和d2分別代表螺孔直徑和螺桿直徑,圖中所示位移為螺桿與螺孔之間的相對位移。

圖3 簡化的螺栓滑移數值模型Fig.3 Schematic of the newly proposed bolt-slip model

圖4 接觸單元的荷載位移曲線Fig.4 Force-deflection curve of the contact element

摩擦對于高強度螺栓的力學行為有至關重要的影響,對于摩擦型高強度螺栓來說,滑移意味著螺栓連接的失效。因此,在所提出的滑移數值模型中,必須精確考慮摩擦的影響。本文采用COMBIN39單元對螺栓與鋼板之間的摩擦行為進行模擬,為了精確模擬復雜的摩擦行為,將摩擦單元的荷載位移曲線設置為如圖5所示。通過增加摩擦單元的初始剛度,將角θ設置足夠大以接近于90°。同時將KEYOPT(1)設置為1,使卸載路徑與加載路徑平行。

圖5 摩擦單元的荷載位移曲線Fig.5 Force-displacement curve of friction element.

由于彈簧單元自由度的參考坐標系與節點坐標系相同,因此,簡化螺栓滑移模型的有限元模型如圖6所示,摩擦單元的內力方向如圖7所示。本簡化數值模型中不考慮高強度螺栓預緊力。通過式(1)計算螺栓滑移時所對應的摩擦力大小[13]。

圖6 簡化滑移螺栓數值模型Fig.6 Simplified numerical model of the bolted connection

圖7 摩擦單元內力方向示意Fig.7 Direction of friction around the bolt hole

式中:m是摩擦面數量;μ是滑動摩擦系數;P是螺栓的預緊力。

沿整體坐標系x軸方向(加載方向)的摩擦力的合力計算如下:

式中:Fmax是滑移前摩擦單元能夠承受的最大內力;n是沿螺孔壁的節點個數。

將式(2)代入式(1)可得滑移前摩擦單元所能承受的最大內力,即

2 試驗驗證

為了驗證所提出的數值模型的可靠性,將基于精細化模型和簡化模型所得到的計算結果與文獻[9]和文獻[14]中的試驗結果進行對比。對比結果如圖8所示。分析中將被連接鋼板的一端固定,在另一端施加拉力。圖中理論摩擦力是指根據式(1)所得到的被連接鋼板發生相對滑移時的拉力為151kN。從圖中可以看出,基于精細化模型和簡化模型所得到的滑移摩擦力分別為150kN和148kN。滑移后螺栓桿與孔壁接觸,承載力進一步提高。從對比結果可以看出,數值模型計算得到的螺栓連接的極限承載力與試驗結果基本吻合。基于試驗所得到的滑移階段對應的荷載高于數值模型,這主要是預緊力的誤差以及鋼板接觸摩擦系數與數值模型理想摩擦系數的不同所引起;在彈性階段,簡化模型與精細化模型和試驗存在誤差,這是由于所使用的單元類型不同所引起,因為簡化模型所使用的為殼單元,而精細化模型是實體單元;滑移后的誤差則是由于鋼材本構的誤差以及單元類型的不同共同引起。簡化模型完成一次計算需要5 min,精細化模型則需要2 h,由此可以看出簡化模型的優越性。

圖8 不同方法荷載--位移曲線對比Fig.8 Comparison of results derived by different methods

3 簡化數值模型滯回性能對比驗證

為了驗證所提出的螺栓滑移模型的可靠性,分別對兩種不同數值模型施加循環荷載進行滯回性能分析。材料本構采用雙線等向強化模型本構模型。將鋼板所用鋼材的屈服強度和切線模量設為256MPa和0.007E,E為鋼材在常溫下的彈性模量,高強度螺栓的屈服強度和切線模量為940MPa和0,螺桿的極限拉力為295.2kN。兩種模型所得滯回曲線如圖9所示。從對比結果可以看出,兩種模型所得滯回曲線吻合很好,兩種模型均可以準確捕捉高強度螺栓的滑移行為。所得最大滑移摩擦力與精細化模型基本一致,誤差基本控制在5%以內。因此可以驗證本文所提出的螺栓滑移數值模型的可靠性。此外,簡化滑移數值模型完成一次滯回分析所需的時間為0.4h,而精細化數值模型完成一次滯回分析所需時間為18h,計算時間減少為原來的2%。

圖9 滯回曲線對比Fig.9 Comparison of hysteretic curves

4 誤差分析

從對比結果可以看出,所提出的簡化螺栓滑移數值模型在高強度螺栓的滯回性能分析中具有很高的計算精度。但是,簡化的滑移模型畢竟未考慮螺栓的預緊力,因此簡化模型中螺桿的抗剪承載力要比實際的高。對于板厚較大,螺桿直徑較小的高強度螺栓連接可能會得到偏于危險的計算結果。基于此,對各種幾何尺寸的螺栓連接進行了滯回性能分析,以系統研究簡化滑移數值模型的適用范圍。為便于說明,特規定螺栓預緊力、螺桿直徑、螺孔直徑、邊板厚、中板厚和鋼板屈服強度的代表符號分別為F、d1、d2、t1、t2、和fy。

4.1 螺桿和螺孔直徑對計算誤差的影響

為了驗證所提出的螺栓滑移模型在計算具有不同幾何尺寸的高強度螺栓連接抗剪承載力的精確性,以螺桿和螺孔直徑為變化參數,通過改變d1與d2的具體數值大小,系統對比研究了不同螺桿直徑對簡化螺栓滑移模型計算精度的影響,鋼板和螺栓的屈服強度分別為256MPa和940MPa,邊板和中板厚度分別為8mm和16mm。同時,將螺栓預緊力分別設置為155kN和50kN,研究了螺栓預緊力對計算精度的影響。計算結果對比如圖10所示。圖10c和圖10f所示結果為將螺桿與螺孔直徑設置為相同的值,即螺桿與螺孔直徑之間沒有空隙可以允許滑移。在該情況下,簡化滑移模型與精細化模型結果吻合較好,在預緊力為50kN時,計算誤差稍大于預緊力為155kN的情況。同樣,圖10d的計算誤差大于圖10a的計算誤差。該對比結果表明,過小的預緊力會加大簡化滑移模型的計算誤差。圖10b和圖10e表明當螺栓空隙較大時,簡化的螺栓滑移模型同樣具有較高的計算精度。精細化有限元模型的收斂性遠不如簡化的滑移數值模型。精細化有限元模型均存在不收斂現象,因此,滯回曲線只能得到一部分。另外,從圖10所有的對比結果可以看出,精細化模型與簡化模型的前兩個滯回曲線的對比結果高度吻合,后續的滯回環的誤差主要是由于螺栓孔的殘余變形較大。而簡化滑移模型所用的殼單元不能很精確的模擬鋼板的殘余變形,導致在模擬螺桿與鋼板在后續的相互作用中存在誤差,該結論可以從圖10d明顯看出。對于精細化數值模型,當滑移位移較大時,鋼板由于螺栓桿的擠壓會產生較大的塑性變形,引起摩擦力的波動變化,并最終導致位移達到20mm時計算結果發散,而簡化模型的滑移變形達到25mm時依然可以收斂。不過該誤差較小,且所得結果偏于安全。

圖10 不同幾何尺寸下滯回曲線對比Fig.10 Comparison of hysteretic curves corresponding to different geometrical size

4.2 鋼板厚度對計算誤差的影響

由于簡化的滑移數值模型未考慮預緊力的影響,因此可能人為的提高了螺栓的極限強度。當鋼板取不同的厚度時,可能會導致不同的計算誤差。因此,本節系統研究了不同鋼板厚度對簡化滑移模型計算精度的影響。以中板和邊板的厚度為變化參數,通過改變t1與t2的具體數值大小,系統對比研究了不同板厚對簡化螺栓滑移模型計算精度的影響,螺栓預緊力F大小為155kN,螺桿和螺孔直徑分別為20mm和21.5mm,鋼板和螺栓的屈服強度分別為256MPa和940MPa。不同板厚下兩種模型計算結果對比如圖11所示。從圖中可以看出,不同板厚對應下的精細化數值模型依然得不到完整的滯回曲線,但當板厚較小時,收斂性得到改善,如圖11d所示。圖11d給出了兩種計算模型對應計算點的應力分布云圖。從圖中可以看出,兩種計算模型的應力分布特征基本一致,螺栓連接的破壞主要是鋼板的擠壓破壞,由于擠壓作用產生塑性變形,該塑性變形可反映在所得滯回曲線中,從而驗證了簡化模型在預測高強度螺栓失效模式方面的精確性。圖11b所示結果誤差最大,經歷大變形后,鋼板會產生塑性變形,由此導致滑移階段摩擦力的改變,由167kN提高到342kN,而簡化模型的摩擦力由148kN提高到214kN。兩者初始階段的計算誤差為11%,且該誤差發生在板厚較厚時(板厚與螺栓直徑比為1.6),當板厚較薄時,誤差很小。因此,對于厚鋼板,簡化的螺栓滑移模型依然可以保持較高的計算精度。

4.3 屈服強度對計算誤差的影響

研究了鋼板屈服強度對簡化模型計算精度的影響。將鋼板的屈服強度fy設置為不同的具體數值,切線模量保持不變。螺栓預緊力F大小為155kN,邊板和中板厚度分別為8mm和16mm,螺桿和螺孔直徑分別為20mm和21.5mm,螺桿屈服強度為960MPa。不同屈服強度下的對比結果如圖12所示。從計算結果可以看出,隨著鋼板屈服強度的提高,精細化有限元模型的收斂能力下降,但是簡化的螺栓滑移數值模型未受到影響。從對比結果可以看出,簡化的螺栓滑移數值模型的計算精度不會受到鋼板屈服強度的影響。從圖12c所示結果可以看出,精細化模型在滑移變形達到10mm時,發生不收斂,而簡化模型在變形達到25mm時依然可以收斂。

圖11 不同板厚所得結果的比較Fig.11 Comparison of the results derived by different thickness

從圖11和圖12中可以看出簡化模型和精細化模型計算結果存在系統性的偏大或偏小的情況,其原因可以歸結為兩方面,對于發生在螺栓連接有較大剪切變形時的誤差,此時鋼板的螺栓孔已發生塑性變形,螺栓孔增大,螺栓孔周圍的接觸應力也發生改變。而簡化模型不能考慮鋼板塑性變形所帶來的影響,因此會產生誤差;誤差的另一個原因是接觸單元以及螺栓桿預緊單元所施加預緊力的偏差。但是誤差總體較小,基本可以忽略。

圖12 不同屈服強度所得結果的比較Fig.12 Comparison of the results derived by different yield strength

5 簡化模型應用

為了驗證所提出的簡化高強度螺栓滑移數值模型的高效性與精確度,基于螺栓滑移模型建立了全螺栓隔板貫通節點的數值模型,如圖13所示。

圖13 全螺栓隔板貫通節點數值模型Fig.13 Finite element model of fully bolted diaphragm-through connections

該模型尺寸采用文獻[15-16]中SJ-1的尺寸。全部螺栓采用10.9級M24扭剪型高強度螺栓,螺栓孔直徑為26mm。為了減小計算代價,將應力較小的柱頂和柱底部分采用線單元建立,建立方法不再贅述,具體可參考文獻[17]。節點部分所有螺栓的建立可采用循環程序建立,進而可以減小模型建立所用時間。節點的等效應力云圖如圖14所示。從圖中可以看出,鋼梁的應力遠高于鋼柱的應力水平,節點的耗能能力由鋼梁與隔板的相對滑移以及鋼梁的塑性變形提供。從圖14中可以準確地觀察到螺栓的滑移現象,在整個計算過程中并未出現不收斂現象,而且計算效率得到了提高。將簡化模型所得滯回曲線與試驗進行對比,如圖15所示。

圖14 節點區域應力云(單位:MPa)Fig.14 Von Mises stress of nodal domain(Unit:MPa)

圖15 不同方法所得滯回曲線對比Fig.15 Comparison of results derived by different methods

橫軸測角是指柱頂水平位移與柱高度的比值。從結果可以看出,數值模型所得結果與試驗結果吻合較好,由于螺栓滑移所引起的剛度退化階段可以得到精確的反映。基于簡化模型所得到的滯回環可以精確反映由螺栓滑移所引起的“捏縮”效應。因此,本文所提出的簡化數值模型可以精確地用于大型鋼結構的抗剪連接中。在保證計算結果準確性的同時,實現了在整體結構中精確考慮螺栓滑移的影響。

6 結論

基于通用有限元軟件提出了一種簡化的螺栓滑移數值模型,并提出了利用非線性彈簧單元精確考慮高強度螺栓摩擦力的數值計算方法,推導并得到了摩擦彈簧單元實常數的計算公式。通過將簡化滑移模型的計算結果與精細化模型計算結果的對比驗證了所提出的數值模型的精確性。同時計算結果表明,通過所提出的數值模型,高強度螺栓的滑移可以被精確的模擬,結果的收斂性得到本質性的改善,計算時間減少為原來的10%。

誤差分析的研究結果表明,精細化模型與簡化模型的前兩個滯回曲線的對比結果高度吻合,后續的滯回環的誤差主要是由于螺栓孔的殘余變形較大。而簡化滑移模型所用的殼單元不能很精確的模擬鋼板的殘余變形,導致在模擬螺桿與鋼板在后續的相互作用中存在誤差。

計算誤差隨著鋼板厚度的增加而增大,當板厚是螺桿直徑的3倍左右時,簡化的螺栓滑移模型依然可以保持較高的計算精度。簡化的螺栓滑移數值模型的計算精度不會受到鋼板屈服強度的影響。

作者貢獻聲明:

趙中偉:負責建立模型,論文寫作。

樊雄濤:數值分析,結果總結,論文修改。

吳剛:負責總體理論研究方向的把控與評價。

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 在线视频亚洲色图| 四虎亚洲国产成人久久精品| 亚洲日韩精品无码专区97| 亚洲AV无码一区二区三区牲色| 免费看a级毛片| 女同国产精品一区二区| 国产人人射| 免费a级毛片视频| 内射人妻无码色AV天堂| 欧美一级高清视频在线播放| 日韩欧美中文字幕在线精品| 亚洲久悠悠色悠在线播放| 国产中文在线亚洲精品官网| 国产亚洲美日韩AV中文字幕无码成人 | 国产精品太粉嫩高中在线观看| 国产精品 欧美激情 在线播放| 午夜性刺激在线观看免费| 国产成人亚洲欧美激情| 99精品国产自在现线观看| 白浆视频在线观看| 国产高清在线丝袜精品一区| 欧美在线一级片| 91麻豆精品国产91久久久久| 91福利一区二区三区| 综合色婷婷| 欧美午夜在线播放| yy6080理论大片一级久久| 国产手机在线ΑⅤ片无码观看| 91青青草视频在线观看的| 日韩高清无码免费| 久久96热在精品国产高清 | 毛片a级毛片免费观看免下载| 日韩A级毛片一区二区三区| 亚洲日韩精品伊甸| 三级毛片在线播放| 午夜啪啪网| 亚洲日韩精品无码专区97| 18禁不卡免费网站| 黄色成年视频| 天天做天天爱天天爽综合区| 欧美三級片黃色三級片黃色1| 精品国产一区91在线| 国产欧美日韩资源在线观看| 国产精品女主播| 99视频精品全国免费品| 无码人妻热线精品视频| 亚洲第一成年网| 就去色综合| 国产流白浆视频| 99精品视频九九精品| 国产福利在线免费观看| 国产一级妓女av网站| 亚洲欧美在线综合图区| 亚洲成在线观看| 久热这里只有精品6| 无码精品一区二区久久久| 中文精品久久久久国产网址| 亚洲色精品国产一区二区三区| 亚洲成a人片| av天堂最新版在线| 国产精品第一区在线观看| 国产www网站| 久久久亚洲国产美女国产盗摄| 伊人蕉久影院| 香蕉99国内自产自拍视频| 国产精品欧美日本韩免费一区二区三区不卡 | 国产精品成人久久| 77777亚洲午夜久久多人| 欧美成人午夜在线全部免费| 日韩一级毛一欧美一国产| 一级毛片免费观看不卡视频| 九九久久精品国产av片囯产区| 四虎永久在线| 美女被操黄色视频网站| 91精品国产一区| 久久这里只精品国产99热8| 啪啪免费视频一区二区| 97国产精品视频自在拍| 亚洲中文字幕久久无码精品A| 久久精品中文字幕免费| 国产精品第页| 久热这里只有精品6|