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

基于Chaboche硬化模型的304SS全壽命循環力學行為仿真分析

2022-07-07 13:14:54劉士杰劉繼超梁國柱
火箭推進 2022年3期
關鍵詞:模型

劉士杰,王 召,劉繼超,梁國柱

(1.北京航天動力研究所 低溫液體推進技術實驗室,北京 100076; 2.北京航天動力研究所,北京 100076; 3.北京機床研究所有限公司,北京 100102; 4.北京航空航天大學 宇航學院,北京 102206)

0 引言

304SS在航空航天、石油化工等行業被大量應用,由304SS制造的結構在非對稱循環載荷作用下會表現出多種形式的力學特性,比如蠕變和棘輪效應。Chaboche本構模型可以較好地描述金屬材料的循環力學行為,因此,它被廣泛地用于不銹鋼結構疲勞壽命評估程序。但對于工程人員而言,Chaboche混合型硬化模型參數的選取比較困難,而且所選取的模型參數是否可以代表全壽命循環的力學行為的研究鮮有報道,這在一定程度上限制了人們對模型的理解和應用。因此,本文開展了基于Chaboche硬化模型的低循環應變載荷下304SS全壽命循環力學響應仿真的研究。

1956年,Prager首先提出了一個線性硬化模型,它可以模擬包辛格效應,但這個模型卻因為在拉伸和壓縮的塑性階段具有相同的硬化模量而無法模擬棘輪效應。鑒于此,1998年,Frederick和Armstrong通過在線性Prager硬化模型中添加非線性動態恢復項提出了Armstrong-Frederick(AF)硬化模型,該模型實現了對應變路徑瞬態記憶效應的模擬,可以仿真材料的棘輪應變行為,但AF硬化模型僅能模擬穩態棘輪應變。為了提高棘輪應變的預測能力,在AF模型的基礎上提出了大量的利用非線性微分方程描述隨動硬化變量演化的硬化模型。其中Chaboche非線性隨動硬化模型(CHK-,是模型包含的AF模型的個數)最為經典,該模型由多個AF模型疊加而來,每個AF模型在整個應力應變響應階段起不同的作用。CHK-是率無關的且能很好地模擬包辛格效應與棘輪應變行為,該模型的一大優勢是,經過修改后它適用于材料在各種復雜載荷場合下復雜力學行為的模擬。因此,該模型得到了大量的工程應用,如Aguis等利用多目標遺傳算法優化獲得Chaboche模型參數,提高了P-3C航空部件疲勞壽命預測的準確性。同時,研究者對優化參數進行了敏感性分析,給出了基于數據的結構疲勞壽命分析指南。Badnava等通過考慮不同加載條件,參考松弛和多級應變控的疲勞試驗,獲得了兩個適應性函數,然后利用遺傳算法優化獲得了Chaboche混合硬化模型參數。Chaboche等建議Chaboche隨動硬化模型至少應包含3個AF硬化律,這可以提高不同應變范圍下模型的預測精度。2017年,Liu等采用Levenberg-Marquardt非線性優化算法對304SS的CHK-3 Chaboche隨動硬化模型參數進行了識別,得到的模型參數較好地模擬了304SS的穩定遲滯環。ASME 蒸汽鍋爐和壓力管道設計要求對核能使用高溫零件的非彈性力學行為進行分析,但缺乏本構模型作為指導,針對該問題,Phan等研究了316H不銹鋼非彈性本構建模過程,以實驗數據平均量為依據提出的本構模型對蠕變和率相關的高溫塑性行為進行了表示,這可為高溫構件的本構方程選取和建模提供參考。為了研究304SS的棘輪應變演化過程,2019年,Liu等利用Chaboche彈塑性本構模型對304SS應變控和應力控的力學行為進行了優化分析,得到的模型參數不僅可以仿真其應力控制的大部分棘輪行為,而且還可以模擬304SS的應變控制變形過程,但對前四分之一個循環的預測精度較差。除此以外,針對具體材料專有力學行為預測的改進Chaboche硬化模型方法得到了研究。2021年,Li等建立了10%Cr馬氏體鋼循環載荷下的初次蠕變再生的統一彈黏塑性本構模型,理論分析與試驗表明該模型可以可靠地描述初次蠕變再生現象和不同載荷的蠕變敏感性。Zhou等提出了9Cr鋼基于物理的高溫疲勞裂紋萌生預估方法,該方法融合了滑移帶形成的Tanaka-Mura模型和統一循環黏塑性本構模型,包含了沉淀硬化、晶界強化等關鍵材料硬化機理,成功地預測了9Cr鋼的循環軟化行為。同時,該項研究工作可為焊接結構的壽命預測提供參考。Meyer等研究了塑性各向異性演化機理,對現有扭曲硬化模型的分析揭示了應力驅動和應變驅動模型的差異,并提出了熱力學自洽和保證屈服面凸性的應力驅動模型,與之前模型相比,該模型在不需要額外參數的情況下更好地吻合試驗數據。

綜上可知,Chaboche硬化模型得到了廣泛研究和應用,可以同時實現對304SS應變控和應力控力學行為特性的模擬,但對全壽命循環,尤其前四分之一個循環的預測精度較差。為此,本文對基于Chaboche硬化模型的低應變循環載荷下(即±0.9%應變范圍內)304SS全壽命疲勞曲線仿真的可行性進行了研究。本文可以為其他金屬材料Chaboche彈塑性本構模型的研究提供一定的參考。

1 研究思路

低循環載荷下304SS全壽命循環力學響應的研究流程如圖1所示。

圖1 材料循環力學特性研究流程Fig.1 Study flowchart of cyclic material mechanical properties

由圖1可知,本文所指的低循環載荷下304SS力學特性研究主要由3部分組成。

1)問題分析與原因定位。根據試驗結果,對304SS材料的Masing效應、屈服應力演化和平臺屈服等基本力學行為進行分析,找出引起前四分之一個循環屈服平臺仿真誤差較大的原因。

2)解決方案與算法實現。研究Chaboche隨動和混合硬化模型模擬穩定遲滯環的算法實現流程,根據屈服平臺仿真誤差較大的原因,研究適用于前四分之一個循環力學行為仿真的模型,并針對本文的研究目的提出基于兩套硬化參數的分段仿真算法。

3)仿真驗證。對304SS±0.8%應變控制循環載荷下的全壽命力學響應曲線進行仿真,驗證算法對屈服平臺效應仿真的適應性。

2 仿真誤差分析與原因定位

2.1 Chaboche混合硬化模型

文獻[10]利用Chaboche混合硬化模型來表達304SS變形過程中的應力應變行為,其中,描述材料受力后狀態的Von Mises屈服準則為

(1)

式中:為偏應力張量;為背應力張量;為初始屈服應力;為描述等向硬化的阻應力,表示為

(2)

式中:和為材料常數,由應變控制的疲勞試驗獲取;為累積塑性應變。當阻應力的初始值為0時,對式(2)積分可以得到

=(1-e-)

(3)

由式(3)可知,值的符號決定了阻應力的正負,即決定了材料的循環硬化、循環軟化特性。

Chaboche隨動硬化模型是由多個AF硬化模型疊加而來,即

(4)

式中:為第個背應力分量所對應的材料常數;d和d分別表示塑性應變增量和累積塑性應變增量;表示使用的AF硬化模型的數目。

(5)

由式(5)可知,當→∞,應力逼近但不超過飽和值

(6)

(7)

式(5)~式(7)以及|-|-=0 (=3/2)是研究Chaboche隨動硬化模型的基礎。

由式(4)~式(7)可知,Chaboche硬化模型的應用是一種逼近方法,因此,它的精度與選取的硬化模型數量有關。Chaboche首先使用3個AF模型疊加來研究材料的單調拉伸曲線,每個AF模型所擔任的角色為

1) 第一個AF模型的和較大,它主要模擬初始階段(小應變區)的應力應變響應;

2) 第二個AF模型的和也較大,此時第一個AF模型已穩定,它主要模擬中間階段(中等應變區)的應力應變響應;

3) 第三個AF模型的和不算大,此時第一和第二個AF模型已穩定,它主要模擬最后階段(大應變區)的應力應變響應,以往的研究稱第三個AF模型為棘輪效應控制模型。

關于Chaboche硬化模型的應力控與應變控仿真算法參見文獻[10]。

2.2 304SS循環力學行為分析

基于文獻[8]和文獻[10]的研究,本文開展了多種應變控制狀態下304SS的循環力學行為試驗,以此分析304SS的循環力學行為。304SS材料中各化學元素質量百分比如表1所示。

表1 304SS化學成分質量百分比Tab.1 Chemical compositions of 304SS 單位:%

為了使304SS材料組織均勻以獲得預期的機械性能,對其采取的熱處理工藝是:在1 160 ℃下將鋼錠鍛造成棒材;然后進行退火處理,加熱到720 ℃,隨后非常緩慢地冷卻到室溫。

2.2.1 低應變載荷下的Masing效應

圖2是304SS在±0.9%對稱循環應變范圍以內的應變控試驗曲線。由圖2可知,控制載荷在低對稱應變范圍內304SS表現出Masing效應。與之類似,文獻[16-17]對304LN的Masing效應進行了分析,結論如下:室溫下,在±0.85%應變控制范圍內,材料表現出明顯的Masing效應,這與本文研究結果相近。文獻[18]對這一現象的機理進行了解釋:室溫下,馬氏體相變和位錯是造成這種材料在不同對稱應變范圍控制下表現出Masing和Non-Masing效應的根本原因。

圖2 應變范圍±0.9%以內應變控制疲勞試驗曲線Fig.2 Fatigue testing curves of symmetric controlling strain within±0.9%

對于大應變控制的情形,即對稱應變范圍大于±0.9%,304SS不僅表現出Non-Masing效應,而且它循環階段的屈服應力在不斷變化,這一現象可以通過疊放拉伸段重合部分來進行觀察,如圖3所示。

圖3 應變范圍±1.0%、±1.1%和±1.2%拉伸段疊放的εp-σ曲線Fig.3 Overlapping of εp-σ tensile curves under strain range of±1.0%,±1.1% and±1.2%

由圖3可知,在大應變控制下,304SS的屈服應力變化比較明顯,在不改變屈服應力的條件下,很難利用一組Chaboche硬化模型模擬這種演化行為。為此,本文只對在±0.9%應變范圍內表現Masing效應的304SS循環力學行為進行研究。

2.2.2 屈服平臺效應

屈服平臺效應是材料在拉伸或者壓縮階段表現出來的應變增加但應力近似不變的力學現象,尤其在前四分之一個循環中表現得尤其突出,這會顯著降低硬化模型參數對材料疲勞全壽命循環曲線的預測精度。

圖4是±0.8%應變控制下304SS的循環應力應變試驗結果。

圖4 304SS的循環應力應變試驗Fig.4 Cyclic stress-strain test for 304SS

由圖4可知,304SS在前四分之一個循環中表現出了明顯的屈服平臺效應。除此以外,304SS表現出等向硬化特性,經歷大約14個循環應力-應變曲線后基本處于穩定遲滯狀態。同時,由圖4(b)可知:①循環載荷下,304SS無明顯應變強化現象;②材料力學屈服強度對應的應力要比圖示轉折點的應力稍大。

2.3 仿真誤差分析

2.3.1 初始屈服應力分析

工程中,一般稱下屈服應力為材料的屈服強度,它對應著標準試樣測試標距內發生了0.2%的塑性應變,而材料力學本構建模中的初始屈服應力可以理解為發生某一數值的塑性應變時的應力,二者的區別要特別注意,這直接決定了分析結果的正確性。在文獻[19]中,初始屈服應力對應的塑性應變取為0.01%,如圖5所示。

圖5 初始屈服應力與屈服強度示意圖Fig.5 Diagram of the first yielding stress and yield strength

綜上可知,初始屈服應力和屈服強度的定義明確,但初始屈服應力的確定方法比較“模糊”。實際上,由圖3可知,屈服應力是一個變化的量。文獻[19]利用對數硬化模型建立了屈服應力的演化方程,較為準確地模擬了材料的后繼屈服應力。本文低循環應變載荷下304SS屈服應力的演化不太明顯,因此,此處假定初始屈服應力和后繼屈服應力相等。對于這種情況,一般有兩種方法確定初始屈服應力:①基于圖3或圖4(b)中的塑性應力—應變曲線轉折的方法,即根據轉折點確定出2倍的初始屈服應力,由此可得初始屈服應力滿足2≈450 MPa,即≈225 MPa;②取塑性應變為定值的方法,如文獻[20]取0.002 5%的塑性應變對應的應力為初始屈服應力。這兩種方法確定的應力差別較小。

同時,由圖4(b)可知,304SS的彈性模量=183.5 GPa,屈服強度≈400 MPa。顯然,初始屈服應力遠低于屈服強度。而在材料循環力學行為分析中,為了使用一套Chaboche硬化模型參數模擬絕大多數(除前四分之一個循環)循環曲線,需要基于Chaboche硬化律的屈服準則使用初始屈服應力來構造屈服面的演化方程,這樣的方法很難捕捉到前四分之一個循環準確的曲線“轉彎”狀態。

2.3.2 仿真誤差分析

文獻[10]利用發展的偽貢獻數法和試錯法,得到了可以同時近似模擬應力控制和應變控制下304SS應力應變響應的Chaboche混合硬化參數,見表2。

表2 Chaboche混合硬化模型參數Tab.2 Parameters of Chaboche combined hardening model

圖6(a)是文獻[10]利用表2得到的棘輪應變模擬結果,圖6(b)是應變控制下循環應力應變行為仿真結果。

圖6 文獻[10]仿真結果Fig.6 Simulation results from Ref.[10]

由圖6可知,文獻[10]給出的Chaboche混合硬化模型參數實現了對304SS棘輪效應和應變控制下應力應變行為的模擬,但對前四分之一個循環的模擬能力較差,這主要是由初始屈服應力和屈服強度的差異造成的。必須使用初始屈服應力構造屈服函數,才可以模擬前四分之一個循環后的后繼屈服循環曲線,而為了模擬前四分之一個循環就需要使用屈服強度來構造屈服函數。除此以外,前四分之一個循環的變化趨勢與后繼循環曲線明顯不同,無法用一套Chaboche隨動硬化模型來對二者進行模擬,因此,有必要討論能對前四分之一個循環進行仿真的可行模型。

3 全循環力學行為仿真分析

3.1 前四分之一個循環的模擬

3.1.1 基于Ramberg-Osgood模型的模擬

前四分之一個循環表示的是材料的單軸拉伸過程,可以利用Ramberg-Osgood方程[式(8)]實現對該段的模擬。

(8)

式中、ε、均為Ramberg-Osgood方程常數。利用Levenberg-Marquardt (L-M)非線性最小二乘優化算法得到如表3所示的前四分之一個循環的Ramberg-Osgood方程常數。同時,為了驗證這些參數對大于±0.9%應變控試驗前四分之一個循環的模擬能力,采用試錯法得到了±3.0%試驗中前四分之一個循環的Ramberg-Osgood方程常數,見表3。

表3 前四分之一個循環的Ramberg-Osgood常數Tab.3 Ramberg-Osgood constants of the first quarter tensile cycle

由表2和表3可知,模擬前四分之一個循環和其他循環的主要區別在于曲線“轉彎”點的應力值發生了改變。更接近材料的屈服強度,這符合2.3.2節所述的前四分之一個循環拉伸段仿真誤差的原因。

為了考察表3模型參數對前四分之一個循環模擬的適用性,圖7(a)給出了利用表3中的參數對±0.8%試驗前四分之一個循環的仿真結果。圖7(b)給出了修正參數對±3.0%試驗前四分之一個循環的仿真結果。由圖7可知,表3的數據較好地實現了對±0.8%控制下前四分之一個循環的仿真,但從±0.8%前四分之一個循環得到Ramberg-Osgood方程常數并不適用于±3.0%的循環響應曲線。二者之間的差別主要在于形狀控制參數和比例參數不同,它們共同決定了Ramberg-Osgood屈服點附近的走勢。

圖7 基于Ramberg-Osgood方程的前四分之一段仿真Fig.7 Simulation for the first quarter cycle based on Ramberg-Osgood equation

需要說明的是,±0.8%應變控的疲勞試驗條件為室溫、三角波、試驗頻率0.25 Hz;±3.0%應變控的疲勞試驗條件為室溫、三角波、試驗頻率0.05 Hz。本文的分析并未考慮率效應對結果的影響,這是以后有待分析的一個問題。

3.1.2 基于Chaboche隨動硬化模型的模擬

前四分之一個循環拉伸段具有指數函數的變化趨勢,因此,可以使用Chaboche硬化模型對該段進行模擬。

表4是基于前四分之一個循環優化得到的Chaboche隨動硬化模型參數。

表4 基于前四分之一個循環的Chaboche隨動硬化模型參數Tab.4 Parameters of Chaboche kinematic hardening model based on the first quarter cycle

利用表4中的數據對前四分之一個循環進行了仿真驗證,并與穩定遲滯環進行了圖示的對比分析,如圖8所示。

仿真使用屈服強度=380 MPa(與前文的400 MPa略有差別,此處是經試錯法調整得到)構造屈服函數,等向硬化模型參數同文獻[10]。

由圖8可知:

圖8 前四分之一個循環以及與穩定遲滯環的比較Fig.8 Simulation of the first quarter cycle and stabilized hysteresis loop

1)不能通過修改初始屈服應力為屈服強度來提高文獻[10]中的模型參數對±0.8%試驗中前四分之一個循環的模擬的準確度。

2)利用表4中的Chaboche隨動硬化模型參數,并選取較為接近前四分之一個循環宏觀“轉折點”(屈服強度附近)的應力構造屈服函數,可以較好地模擬前四分之一個循環的應力應變響應。

基于以上分析發現,可以較好地分別利用Chaboche硬化模型實現對前四分之一個循環和穩定遲滯環的仿真。但可否利用Chaboche硬化模型實現對全壽命循環曲線的仿真是值得討論的問題。

作為對比,本文給出文獻[10]中304SS混合硬化模型參數,如表5所示。

表5 文獻[10]中的Chaboche隨動硬化模型參數Tab.5 Chaboche kinematic hardening parameters from Ref.[10]

3.2 全循環周期模擬

3.2.1 算法

由第3.1節的分析可知,本文給出了兩種適合模擬前四分之一個循環曲線的方法,它們的優缺點概括如下:

1)Ramberg-Osgood模型適于模擬單調拉伸、循環應力應變曲線,以及疲勞試驗的前四分之一個循環曲線,但對于整個疲勞循環曲線的仿真難度較大,這是由循環硬(軟)化和包辛格效應等因素共同決定的;

2)Chaboche混合硬化模型不僅可以模擬前四分之一個循環的后繼屈服力學響應曲線,而且還可以較準確地模擬前四分之一個循環的曲線,只是無法使用一套模型參數來同時模擬這兩種循環力學狀態。

除此以外,Chaboche硬化模型更便于程序實現,而且具有成熟的算法(如基于徑向回退算法的應力應變仿真算法)和程序(如ANSYS/ABAQUS內置本構模型)。圖9是利用表5的參數對304SS穩定遲滯環進行的有限元仿真結果。

圖9 試棒仿真分析Fig.9 Simulation analysis of test specimen

基于以上分析,本文提出的低循環載荷下304SS全壽命周期循環力學響應曲線的模擬方案為:①在前四分之一個循環采用表4的數據仿真前四分之一個循環,如圖8所示;②在后繼循環中,尤其對穩定遲滯環,使用表5的參數進行仿真。兩個階段的計算都是Chaboche硬化模型的基本應用,關鍵是通過如下算法(單軸各向同性狀態)實現模型參數的轉換。

3)是否第一次卸載?否,則flag=1;是,則flag=flag+1。

4)若flag=1,則參數取表3;若flag=2,則參數取表4;若flag=其他,則報錯。

5)計算如下預測應力和預測屈服函數

6)回退法公式如下

+1=+(-

+1=

模型參數的轉換是由載荷形式控制的[即步驟3)],而且僅在前四分之一個循環和后繼第一個反向段起作用。

3.2.2 結果分析

基于3.2.1節算法,利用MATLAB軟件編程實現了對304SS±0.8%應變控制下全壽命循環力學行為的仿真,其中,模型參數取自表4,仿真結果見圖10。

圖10 304SS±0.8%應變控制下全壽命循環力學行為模擬Fig.10 Full-life cycle simulation of 304SS mechanical behavior under±0.8% strain controlled

盡管圖10所示的前四分之一個循環和后繼曲線(除前四分之一個循環后的壓縮段)較好地得到了仿真,但這個仿真是不可取的,原因是表4和表5參數分別控制了仿真過程中前四分之一個循環和后繼曲線中背應力的演化形式,前四分之一個循環的隨動硬化背應力只有圖10所示的40 MPa左右,而圖6(b)所示的背應力高達205 MPa,如果按照表4的參數進行仿真,那么圖10所示的前四分之一循環后的卸載段初始點設置背應力為205 MPa,其他參數保持不變,則得到的仿真結果如圖11所示。

圖11 調整前四分之一段尾背應力的全壽命循環力學行為模擬Fig.11 Simulation of life cycle mechanical behavior by adjusting the tail back-stress of the first quarter cycle

顯然,這樣的仿真是不合理的。通過以上分析可知,無法結合表6利用兩套Chaboche硬化模型參數來模擬304SS的全壽命循環曲線。

4 結論

通過本文的研究主要得出如下結論。

1)低循環載荷下,304SS表現出明顯的Masing效應和循環硬化特性,±0.8%應變控制下經歷14個應力循環后,應力應變曲線基本處于穩定遲滯狀態。

2)±0.8%應變控試驗中前四分之一個循環對應的Ramberg-Osgood模型常數是:=34.713,=0.002 24,=430 MPa;Chaboche硬化模型參數是:=744 639 MPa,=155 193,=71 633 MPa,=3 014,=20 608 MPa,=1 051,對應的初始屈服應力=380 MPa。

3)Chaboche硬化模型可以分別對前四分之一個循環和后繼循環曲線實現準確仿真,但利用兩套Chaboche硬化模型參數實現對304SS低循環載荷下全壽命循環曲線的仿真方案是不可行的。

本文的研究可以為其他金屬材料Chaboche型本構模型參數識別和仿真研究提供參考。

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 好紧太爽了视频免费无码| 国产精品护士| 色天天综合| 国产精品久久久久久久久久久久| 欧亚日韩Av| 欧美中文字幕第一页线路一| 亚洲中文字幕av无码区| 91娇喘视频| 有专无码视频| AV熟女乱| 亚洲无线观看| 91丨九色丨首页在线播放| 尤物精品国产福利网站| 日本成人不卡视频| 精品自窥自偷在线看| 人妻91无码色偷偷色噜噜噜| 亚洲人成网站在线观看播放不卡| 九九九久久国产精品| 日韩东京热无码人妻| 亚洲a免费| 国产二级毛片| 波多野结衣AV无码久久一区| 久久毛片网| 91年精品国产福利线观看久久 | 亚洲第一中文字幕| 国产精品午夜福利麻豆| 欧美啪啪精品| 免费国产不卡午夜福在线观看| 美女高潮全身流白浆福利区| 免费无码AV片在线观看国产| 免费又爽又刺激高潮网址| 99ri精品视频在线观看播放| 国产靠逼视频| 亚洲无码A视频在线| 不卡网亚洲无码| 国产无吗一区二区三区在线欢| 高清欧美性猛交XXXX黑人猛交 | 免费高清a毛片| 丰满人妻久久中文字幕| 亚洲精品无码AV电影在线播放| 色偷偷综合网| 免费高清a毛片| 国产又大又粗又猛又爽的视频| 国产国模一区二区三区四区| 亚洲国产高清精品线久久| 国产成人三级| 国产本道久久一区二区三区| 91久久国产热精品免费| 无码啪啪精品天堂浪潮av| 中文字幕亚洲专区第19页| 国产成人1024精品| 亚洲国产日韩一区| 制服丝袜无码每日更新| 国产香蕉一区二区在线网站| 国产无码在线调教| 99精品视频播放| 亚洲香蕉伊综合在人在线| 香蕉综合在线视频91| 99视频在线免费看| 欧美日韩成人| 国产9191精品免费观看| 激情网址在线观看| 国产欧美网站| 国产91九色在线播放| 午夜三级在线| 丰满的少妇人妻无码区| 亚洲天堂啪啪| 男人天堂伊人网| 综合亚洲网| 亚洲第一av网站| 亚洲精品图区| 国产日本一区二区三区| 午夜精品福利影院| 青草91视频免费观看| 少妇露出福利视频| 最新国产麻豆aⅴ精品无| 国产一级毛片网站| 99久久国产综合精品2020| 亚洲男人在线天堂| 欧美性久久久久| 四虎AV麻豆| 永久在线播放|