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

固體火箭發動機滾動球窩噴管強化接頭材料性能分析①

2012-07-09 09:12:46劉文芝武建新杜俊碩程佳慶
固體火箭技術 2012年4期
關鍵詞:有限元

劉文芝,武建新,杜俊碩,程佳慶

(1.內蒙古工業大學機械學院,呼和浩特 010051;2.中國航天科工集團公司六院四十一所,呼和浩特 010051)

固體火箭發動機滾動球窩噴管強化接頭材料性能分析①

劉文芝1,武建新1,杜俊碩1,程佳慶2

(1.內蒙古工業大學機械學院,呼和浩特 010051;2.中國航天科工集團公司六院四十一所,呼和浩特 010051)

為提高某固體火箭發動機滾動球窩接頭接觸承載能力,降低系統擺動力矩,以單珠承載試驗模型為計算模型,采用摩擦接觸問題的Lagrange乘子法與彈塑性耦合的有限元理論,計算分析了具有不同強化層材料性能的陽球試件與滾動體接觸位置處的接觸應力、摩擦應力、變形及破壞機理;解決了陽球不同強化層材料性能對接頭接觸承載性能的影響問題。根據計算分析結果,提出了滿足彈塑性摩擦接觸性能條件下,小的摩擦系數、小的彈性模量和較大的屈服極限相配合的陽球強化層材料性能要求;通過與單珠承載試驗對比分析,檢驗了有限元建模及計算結果的合理性。

固體火箭發動機;滾動球窩噴管;強化接頭材料性能;彈塑性摩擦接觸有限元;單珠承載試驗

0 引言

某固體火箭發動機全軸擺動滾動球窩噴管的滾動球窩接頭由大尺寸球面的陰球、陽球及數個小尺寸球形滾動體等組成。接頭既是活動體與固定體間的連接件,又是載荷支承件,其接觸承載能力是決定接頭結構設計及系統功能發揮的關鍵問題之一[1]。

接頭的陰球和陽球零件表面須進行強化處理,強化層材料性能必須使陽球表面具有高的接觸承載能力,否則在工作載荷作用下,滾動體與陽球接觸位置處,陽球表面會產生塑性變形;而系統不允許太大的變形,當變形超過許可范圍,則造成擺動力矩增大,增加推力向量控制的難度。因此,強化接頭材料必須在滿足總體對接頭結構質量及接頭結構尺寸要求的前提下,保證接頭陰、陽球能承受預定的工作載荷而不產生過量的塑性變形。

即使在簡單載荷作用下,對于考慮強化層材料性能和摩擦效應,并與彈塑性耦合的接觸問題的簡單模型,其解析建模及求解也非常復雜困難,一些數學問題亟待進一步解決[2-4]。本文采用更符合工程實際的彈塑性摩擦接觸有限元法,計算分析陽球試件經表面強化處理后,強化層具有不同材料性能時(不同彈性模量、屈服極限和滾動摩擦因數),與滾動體間的接觸應力、摩擦應力、變形及破壞機理;與單珠承載試驗對比分析,為降低系統擺動力矩,增加接觸承載能力,同時考慮大尺寸球面(φ≥300 mm)的強化工藝,解決不同強化層材料性能對接頭接觸性能的影響問題;提出更為合理的滿足接頭彈塑性摩擦接觸承載性能的強化層材料性能要求。

1 彈塑性摩擦接觸有限元計算模型

試驗中沒有觀測到滾動體的殘余變形,因此計算模型為彈性滾動體與彈塑性陽球試件的摩擦接觸模型,摩擦采用庫侖摩擦。

1.1 彈塑性摩擦接觸有限元計算模型

考慮實際問題的特點,選用計算精度高的接觸問題Lagrange乘子法。根據廣義變分原理,在計算域中,把接觸約束條件與Lagrange乘子相乘,與系統原勢能(結構材料的彈塑性勢能)一起構成一個修正泛函,把接觸問題轉換成求修正泛函的駐值問題,即無約束泛函的極值問題。離散接觸界面,滾動體與陽球試件間建立面-面接觸模型。陽球試件表面為接觸面,滾動體表面為目標面。把接觸力轉換成等效節點接觸力向量,寫入方程左端,并整合為系數矩陣,引入接觸定解條件[5-6],得到彈塑性摩擦接觸問題的有限元控制方程。

t+Δt時刻彈塑性摩擦接觸有限元控制方程為

式中t+ΔtKep為t+Δt時刻系統彈塑性總剛矩;Δa為系統的位移增量向量;t+ΔtQl為t+Δt時刻的外載荷向量;t+ΔtQi為t+ Δt時刻的內力向量;為t時刻接觸點對間距;t+Δtλ為t+Δt時刻拉格朗日乘子(力學意義是接觸力)。

設Z、X分別表示接觸體在接觸界面的法線方向和切線方向。

粘合狀態:系數矩陣Kca=。

滑動狀態:Kca≠(有摩擦滑動),Kcλ=(無摩擦滑動),t+Δtλ=t+ΔtλZ,-=-。

如果接觸狀態變化,可能發生不連續,為避免收斂太慢,在每一載荷增量步內用F.N.R迭代求解,其遞推迭代公式為

接觸計算判定條件為

式中t+ΔtR為t+Δt時刻接觸力向量;f為摩擦因數;?a為相對位移向量。

對于彈塑性摩擦接觸問題,由于2個非線性混合在一起,在一個增量步中,每給定一種接觸狀態,就要解一個彈塑性邊值。彈塑性迭代收斂后,檢查接觸狀態,若接觸狀態與求解前一致,則該增量步的解已獲得;不一致,則必須恢復到增量步開始時的狀態,否則就把非真實接觸條件下的變形歷史引入后續計算,造成誤差甚至錯誤的結果。然后,用新接觸狀態重新求解該增量步,直到彈塑性迭代收斂后,接觸狀態與求解前假定的一致,得到該增量步彈塑性摩擦接觸問題的解,即彈塑性迭代置于內循環,接觸迭代置于外循環。

1.2 材料性能及載荷路徑

陽球試件表面強化處理后,經計算和試驗檢測,強化層厚度為 1.22 mm[7],泊松比μ=0.3。本構關系為彈性區域σ=Eεe;塑性區域σ=σs+Hεp,εe、εp分別為彈、塑性應變,強化系數H=1×104MPa,表面強化處理后,其彈性模量E=2×105~3×105MPa,屈服極限σs可達到900 MPa以上。

陽球試件表面經不同強化工藝處理后,其基體材料性能變化不明顯。陽球試件基體材料的彈性模量E=2.1 ×105MPa,泊松比μ=0.3,屈服極限σs=1 100 MPa。本構關系為彈性區域σ=Eεe;塑性區域由試驗測得,見圖1。

陽球與滾動體間的滾動摩擦因數f=0.001~0.005[1]。

滾動體彈性模量E=2.54×105MPa,泊松比μ=0.3;滾動體半徑r=5 mm。

有限元計算分析中的加載、卸載路徑與單珠承載試驗中的加載、卸載路徑一致。

圖1 陽球基體材料的應力-應變曲線Fig.1 Stress-strain curve of the convex sphere materials

2 有限元計算結果與分析

算例為在500 kg載荷作用下,大尺寸球面(φ≥300 mm)的陽球,通過表面層強化后,其強化層材料具有如下屬性:滾動摩擦因數f分別為0.003和0.005;彈性模量E分別為2.1 ×105、2.54 ×105、3 ×105MPa;屈服極限σs分別為1 000、2 000、3 000 MPa。

根據單珠承載試驗中載荷特性和接觸體的幾何特性,建立滾動體與陽球試件彈塑性摩擦接觸軸對稱有限元模型(圖2)。

圖2 有限元模型Fig.2 Finite element model

2.1 接觸應力及接觸變形

圖3~圖5分別為強化層具有不同摩擦系數、彈性模量和屈服極限時,與滾動體接觸位置處,陽球試件表面接觸應力σ及彈塑性接觸半徑a分布、彈塑性接觸壓深b分布、殘余壓深b0及殘余半徑a0分布。

由圖3~圖5可知:

(1)彈性模量、屈服極限相同

摩擦因數分別為0.003和0.005時,接觸應力極值分別為 5 202.9、5 203.7 MPa;彈塑性接觸半徑分別為 0.60、0.59 mm,殘余半徑分別為 0.57、0.56 mm;彈塑性接觸壓深分別為 0.036 828、0.036 821 mm,殘余壓深分別為 0.016 400、0.016 366 mm。

圖3 接觸應力及接觸半徑分布Fig.3 Curves of contact stress and contact radius

圖4 彈塑性接觸壓深分布Fig.4 Curves of elasto-plastic indentation deepness

圖5 殘余壓深及殘余半徑分布Fig.5 Curves of residual indentation deepness and residual indentation radius

由于摩擦的存在,使得材料產生屈服時,需要更高的法向應力。同時,由于摩擦阻止邊界區域的滑動,因此摩擦因數越大,陽球與滾動體間的表面接觸應力越大,接觸變形越小。總體而言,滾動摩擦因數小,陽球試件與滾動體間的彈塑性摩擦接觸應力及接觸變形分布差異較小。但為了減小接觸應力、增加陽球接觸強度、降低系統擺動力矩,應盡量減小摩擦因數。

(2)摩擦因數、屈服極限相同

不同彈性模量的陽球試件,表面接觸應力及變形的分布差異較為明顯。彈性模量分別為2.1×105、2.54×105、3 ×105MPa時,接觸應力極值分別為5 065.4、5 202.9、5324.3 MPa;彈塑性接觸半徑分別為0.62 、0.60、0.59 mm,殘余半徑分別為 0.58、0.57、0.56 mm;彈塑性壓深分別為 0.039 335、0.036 828、0.034 825 mm,殘余壓深分別為 0.016 100、0.016 400、0.016 486mm。

彈性模量越小,彈塑性摩擦接觸應力越小。材料抵抗彈性接觸變形的能力越差;抵抗法向塑性接觸變形的能力增強,殘余壓深減小,有利于系統運動中,接頭承受更高載荷而不產生過量的塑性變形,從而減小系統擺動力矩。

(3)彈性模量、摩擦因數相同

由于屈服極限不同,其接觸應力和接觸變形分布明顯不同。屈服極限越小,材料抵抗塑性接觸變形的能力越差,塑性變形越大,接觸應力越小,接觸應力峰值越向接觸邊緣偏移。

屈服極限分別為1 000、2 000、3 000 MPa時,接觸應力極值分別為 3 446.8、5202.9、6 885.5 MPa;彈塑性接觸半徑為 0.72、0.60、0.56 mm,殘余半徑為 0.68、0.57、0.51mm;彈塑性壓深為 0.051 700、0.036 828、0.030 649 mm,殘余壓深為 0.035 500、0.016 400、0.007 320 mm。

與σs=2 000 MPa相比,材料的屈服極限減小50%,強化層材料很快屈服,接觸應力降低33.752%;彈塑性壓深增大40.382%,殘余壓深增大116.463%。屈服極限增大50%,接觸應力增大32.340%;彈塑性壓深減小16.778%,殘余壓深減小55.366%。高的屈服極限雖然使接觸應力有所增加,但可有效減小塑性接觸變形。

為滿足系統擺動力矩要求,規定強化接頭殘余壓深為 0.03 mm[1];企業設計規范要求[8]實際應用中,其安全系數為1.2,則許用殘余壓深為0.025 mm。強化層屈服極限為1 000 MPa的陽球試件,其表面殘余壓深已超過許用殘余壓深。

2.2 摩擦應力

圖6為強化層具有不同材料性能時,與滾動體接觸位置處,陽球試件表面摩擦應力σf分布。

由圖6可知:

(1)屈服極限相同

屈服極限σs皆為2 000 MPa時,由于滾動體與陽球接觸表面,接觸區整體呈粘合狀態,在接觸區域內,摩擦應力皆為切向拉應力。

彈性模量為 2.54×105MPa,摩擦因數分別為0.003、0.005 時,最大摩擦應力分別為 9.940、14.553 MPa。摩擦因數為0.003,彈性模量分別為2.1×105、2.54×105、3×105MPa時,最大摩擦應力分別為6.274、9.940、13.057 MPa。

摩擦因數、彈性模量越小,接觸表面的摩擦應力越小。

圖6 摩擦應力分布Fig.6 Curves of friction stress

(2)彈性模量、摩擦因數相同

屈服極限對摩擦應力的影響很大。與σs=2000 MPa相比,屈服極限減小50%,強化層很快屈服,在相同載荷作用下產生更大的徑向位移;接觸區域滑動,摩擦應力呈較大的壓應力狀態。當屈服極限達到3 000 MPa時,大的法向接觸壓力產生大的剪切應力,距接觸中心0.38 mm的接觸區域內,產生滑動趨勢,從而產生相對較小的切向摩擦壓應力;摩擦和周圍材料的約束使滑動趨勢逐漸減弱并阻止了滑動,在0.38 mm之外的接觸區域內呈粘合狀態,并產生小的切向摩擦拉應力。

2.3 層下應力

圖7、圖8分別為強化層具有不同材料性能時,陽球試件與滾動體接觸表面下方,沿接觸中心對稱軸z的主剪應力τ和等效Misses應力σv分布。由圖7、圖8可知,強化層材料性能不同的陽球試件,在接觸位置近表層區域內,主剪應力和等效Misses應力迅速增長。沿中心對稱軸距接觸中心z≥2.5 mm的基體材料內,各應力平緩下降,逐漸趨于一致。

(1)彈性模量、屈服極限相同

由于滾動摩擦因數總體較小,不同摩擦因數對接觸表面下方強化層和基體材料內的應力分布影響很小。摩擦因數分別為0.003和0.005時,最大主剪應力和等效Misses應力,分別產生在z=0.275 mm和z=0.272 mm位置處;摩擦使最大主剪應力和等效Misses應力向接觸表面移動。

(2)摩擦系數、屈服極限相同

彈性模量為2.1×105MPa時,最大主剪應力和等效Misses應力出現在z=0.29 mm位置處;彈性模量分別為2.54×105、3 ×105MPa,最大主剪應力和等效Misses應力皆出現在z=0.275 mm位置處。彈性模量越大,最大主剪應力和等效Misses應力越向接觸表面移動,當彈性模量超過2.54×105MPa時,其位置不再改變。

由于彈性模量不同,在強化層底部靠近基體材料的區域內,應力分布有較為明顯的不同;彈性模量越小,在強化層和基體材料交界面處,應力突變越小,越不易引起層間剝離和裂紋[9]的產生。

(3)彈性模量、摩擦系數相同

由于屈服極限不同,強化層內主剪應力和等效Misses應力分布明顯不同;強化層與基體材料的屈服極限越接近,材料交界面處的應力突變越小,但在基體材料局部區域內應力增大。

屈服極限分別為3 000、2 000、1000 MPa,最大主剪應力和等效Misses應力分別出現在z=0.26 mm、z=0.275 mm和z=0.32 mm位置處。強化層屈服極限為1 000 MPa,在z≥0.01 mm的區域內,強化層全部屈服。

(4)強化層材料性能不同的陽球試件,其基體材料均未產生塑性變形。理論上,層間不易產生初始裂紋[9]。

考慮大尺寸球面的強化和加工工藝,強化層厚度在1~1.4 mm[7],具有上述不同材料性能的陽球試件,在不同載荷作用下,與滾動體接觸位置處,應力和變形分布趨勢基本相同。

圖7 主剪應力分布Fig.7 Curves of main shear stress

圖8 等效Misses應力分布Fig.8 Curves of equivalent misses stress

因此,從接觸性能及破壞機理考慮,對于陽球試件的強化層,小的摩擦系數、小的彈性模量、較高的屈服極限,這樣的材料性能相互配合,才能充分發揮強化層和基體材料的綜合機械性能,減小系統擺動力矩。

3 試驗與計算結果的對比分析

根據有限元計算分析結論,對陽球試件的強化工藝提出要求,以保證其材料性能要求。陽球試件經強化處理,表面精磨0.1 mm后,掃描電鏡測得強化層厚度為1.22 mm。彈性模量E=2.54 ×105MPa,滾動摩擦因數f=0.003,屈服極限σs=2 000 MPa。陽球試件基體材料及滾動體的材料性能檢測結果見1.2節。

3.1 單珠承載試驗

滾動體與陽球試件接觸,載荷作用下陽球試件表面產生彈塑性變形,卸載后有不同程度的壓痕。每級加載的壓痕不重復。在工具顯微鏡下測量并記錄壓痕寬,對每級加載條件下的壓痕寬取平均值2ˉa0(表1)。

3.2 試驗與計算結果的對比分析

不同載荷作用下,壓痕的試驗測量結果與彈塑性摩擦接觸有限元計算結果見表1。

表1 試驗結果與彈塑性摩擦接觸有限元計算結果Table 1 Test results and computation results of elasto-plastic friction contact finite element method

每級載荷作用下,相對于試驗測量平均壓痕寬,殘余壓寬的數值計算結果平均誤差為6.15%。

電鏡下觀察壓痕(放大100倍),并切片檢測,試件在250~1 000 kg載荷作用下,均未發現層間剝離和裂紋。

作為算例的某發動機燃燒室最大工作壓強為pcmax=6 MPa,則推力向量控制系統接頭內部,陽球承受燃氣壓強所產生的最大單珠正壓力Pbmax為

式中Rs為陽球半徑;φC1及φ為接觸角;CF為發動機地面試車狀態推力損失系數;Rt為噴喉半徑;n為滾動體數目。

由表1可知,該噴管接頭在燃燒室最大壓強作用下,其接觸表面殘余壓深滿足設計要求。

因此,彈塑性摩擦接觸有限元算法及其結果合理;通過彈塑性摩擦接觸有限元計算分析所提出的陽球試件強化層材料性能要求可靠。

4 結論

(1)彈塑性條件下,滾動摩擦因數主要影響接觸表面的摩擦應力和接觸狀態。彈性模量使接觸表面應力和變形分布、強化層底部局部區域的主剪應力和等效Misses應力分布有較為明顯的不同。小的摩擦因數、彈性模量使彈塑性接觸變形有所增加,但可以有效降低表面接觸應力、摩擦應力和殘余壓深。不同的屈服極限對接觸表面應力、變形、接觸狀態及層下應力分布有很大影響;屈服極限越大,塑性變形越小,可以有效降低系統擺動力矩;但太大的屈服極限,使接觸強度顯著下降,最大主剪應力和等效應力越接近接觸表面。

(2)不同的彈性模量和屈服極限,都會引起強化層和基體材料交界面的應力突變;強化層彈性模量和屈服極限越接近基體材料,應力突變越小,可減小層間產生剝離和裂紋的幾率;但當屈服極限σs=1 000 MPa時,其殘余壓深已超出許可范圍,勢必增大系統擺動力矩。

(3)綜合考慮接觸表面及層下應力分布、變形及層間應力突變,對于強化層,小的摩擦系數、小的彈性模量、較高的屈服極限,這樣的材料性能相互配合,才能充分提高系統強化接頭的接觸性能,降低系統擺動力矩。

(4)通過彈塑性摩擦接觸有限元計算分析,對陽球材料強化工藝提出要求,通過實測,得到該型號滾動球窩噴管強化接頭陽球試件強化層的材料性能:彈性模量E=2.54 ×105MPa,滾動摩擦因數f=0.003,屈服極限σs=2 000 MPa;其材料性能可以承受燃燒室最大工作壓強。

(5)通過與單珠承載試驗對比分析,檢驗了對陽球試件強化層材料性能要求及有限元建模和算法的合理性。

[1]陳汝訓.固體火箭發動機設計與研究[M].北京:宇航出版社,1992.

[2]張焱,孔祥安,金學松,等.輪軌三維彈塑性接觸應力的算法研究[J].力學與實踐,2000,22(1):23-27.

[3]Pascoe S K.Contact stress analysis using the finite element method[D].Uni.Liverpool,1990.

[4]Kral E R,Komvopoulos K,Bogy D B.Elastic-plastic finite element analysis of repeated indentation of an elastic-plastic layered medium by rigid sphere[J].ASME J of Applied Mechanics,1995,62:42-48.

[5]王勖成.有限單元法[M].北京:清華大學出版社,2003.

[6]何君毅,林祥都.工程結構非線性問題的數值解法[M].北京:國防工業出版社,1994.

[7]劉文芝,戴美魁,韋廣梅,等.某固體火箭發動機推力向量控制系統接頭接觸性能分析[J].固體火箭技術,2011,34(5):573-578.

[8]中國航天科工集團公司六院41所.固體火箭發動機設計規范[S].2008.

[9]《合金鋼》編寫組.合金鋼[M].北京:機械工業出版社,1978.

Material performance analysis for the intensified joint of rolling ball joint socket nozzle in solid propellant rocket motor

LIU Wen-zhi1,WU Jian-xin1,DU Jun-shuo1,CHENG Jia-qing2
(1.Mechanical College of Inner Mongolia University of Technology,Hohhot 010051,China;2.The 41th Institute of the Sixth Academy of CASIC,Hohhot010051,China)

To reduce system swing moment and improve contact intensity of rolling ball joint socket nozzle in a certain solid propellant rocket motor's joint,contact stress,friction stress,deformation and failure mechanism between test parts convex sphere and ball were computed,which have different thickness and different material performance intensified layer.Computation model is single ball load-bearing experiment model and algorithm is the Lagrange multiplier method which is used for friction contact question and elasto-plastic coupling's finite element method.Different material performance intensified layers influencing the joint contact performance were figured out.On the basis of computation result,the requirement for material performance intensified layers of the convex sphere was presented,which satisfied friction contact performance,small friction coefficient,small elastic modulus and major yield limit.Rationality of FEM modeling and computation result has been validated by single ball load-bearing experiment and high load contact analysis.

solid propellant rocket motor;rolling ball joint socket nozzle;material performance of intensified joint;elasto-plastic friction contact finite element method;single ball load-bearing experiment

V435+.23

A

1006-2793(2012)04-0522-06

2012-02-24;

2012-05-06。

國家自然科學基金項目(50965013);內蒙古自然科學基金項目(2009BS0702)。

劉文芝(1969—),女,教授,研究方向為機械設計及理論。E-mail:wzliu0162@yahoo.com.cn

(編輯:崔賢彬)

猜你喜歡
有限元
基于擴展有限元的疲勞裂紋擴展分析
非線性感應加熱問題的全離散有限元方法
TDDH型停車器制動過程有限元分析
新型有機玻璃在站臺門的應用及有限元分析
上海節能(2020年3期)2020-04-13 13:16:16
基于I-DEAS的履帶起重機主機有限元計算
基于有限元模型對踝模擬扭傷機制的探討
10MN快鍛液壓機有限元分析
磨削淬硬殘余應力的有限元分析
基于SolidWorks的吸嘴支撐臂有限元分析
箱形孔軋制的有限元模擬
上海金屬(2013年4期)2013-12-20 07:57:18
主站蜘蛛池模板: 国产激爽大片在线播放| 熟妇无码人妻| 狂欢视频在线观看不卡| 色偷偷一区二区三区| 亚洲成人精品久久| 国产成人高精品免费视频| 国内熟女少妇一线天| 中文字幕人妻av一区二区| 婷婷亚洲最大| 国产乱人伦AV在线A| 国产97区一区二区三区无码| 国产第八页| 中文字幕波多野不卡一区| 无码精品一区二区久久久| 99热这里只有精品5| 欧美成人影院亚洲综合图| 国产男人天堂| 久久精品无码中文字幕| 欧美黄色网站在线看| 日韩精品欧美国产在线| 92午夜福利影院一区二区三区| 六月婷婷激情综合| 国产99视频精品免费视频7| 欧美成人怡春院在线激情| 国禁国产you女视频网站| 伊人精品视频免费在线| 日韩欧美国产成人| 亚洲性色永久网址| 亚洲精品动漫| 亚洲性影院| 综合社区亚洲熟妇p| 波多野结衣第一页| 免费播放毛片| 亚洲高清国产拍精品26u| 欧美人人干| 狼友视频一区二区三区| 精品视频一区二区三区在线播| 亚洲AV一二三区无码AV蜜桃| 欧美亚洲一区二区三区导航 | 亚洲精品大秀视频| 原味小视频在线www国产| 国产免费一级精品视频 | 亚洲欧美日韩天堂| 国产微拍一区二区三区四区| 伊大人香蕉久久网欧美| 亚洲无线观看| 女人18毛片久久| 40岁成熟女人牲交片免费| 国产精品不卡片视频免费观看| 精品伊人久久久香线蕉| 欧美日韩中文字幕在线| 亚洲三级电影在线播放| 亚洲欧美国产五月天综合| 国产一级裸网站| 国产人人乐人人爱| 成年人福利视频| 天天操精品| 内射人妻无套中出无码| 亚洲日本www| 91欧洲国产日韩在线人成| 欧洲av毛片| 亚洲妓女综合网995久久| 一本综合久久| 激情综合五月网| 国产色伊人| 四虎成人免费毛片| 在线观看亚洲人成网站| 亚洲综合欧美在线一区在线播放| 日韩黄色精品| 午夜国产精品视频黄| 欧美怡红院视频一区二区三区| 欧美亚洲另类在线观看| av在线5g无码天天| 99色亚洲国产精品11p| 国产精品高清国产三级囯产AV| 日韩精品免费一线在线观看| 欧美人人干| 538精品在线观看| 最新亚洲av女人的天堂| 国产91线观看| 国产精品嫩草影院视频| 国产亚洲精品97AA片在线播放|