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

面向增材制造的進氣道拓撲優化數值研究

2025-04-14 00:00:00郭家琦孫振華劉亞東白濤濤
航空兵器 2025年1期

摘 要:""""" 進氣道是沖壓發動機的重要部件之一, 其強度與剛度對于沖壓發動機穩定性至關重要。 為降低發動機消極質量并優化發動機設計, 本文以某進氣道為研究對象, 在滿足進氣道變形和承載要求的基礎上進行進氣道薄壁結構減重設計。 對進氣道進行靜力校核, 根據進氣道在4種工況下的初始強度校核結果確定優化設計域, 結合點陣結構, 同時考慮加強筋布局優化與加強筋尺寸優化, 運用變密度拓撲優化方法, 采用主級加筋結合局部點陣加筋填充的優化方案對進氣道進行了優化設計。 利用有限元軟件進行優化前后強度、 剛度、 穩定性校核, 優化后進氣道結構質量減少了8%, 穩定性增強, 整體剛度分布更好, 實現了進氣道的減重與結構設計空間的擴展, 證明本文采用的優化方法可行性。

關鍵詞:"""" 拓撲優化; 數值研究; 進氣道; 點陣結構; 薄壁加筋;" 沖壓發動機

中圖分類號:"""" TJ760; V431

文獻標識碼:""" A

文章編號:"""" 1673-5048(2025)01-0107-08

DOI: 10.12132/ISSN.1673-5048.2024.0193

0 引" 言

拓撲優化技術作為結構優化設計的重要分支, 是結構輕量化設計、 獲得高性能創新構型的有效設計方法, 現已被廣泛應用于航空航天領域中1-2。 降低消極質量是優化發動機設計的有效方法, 因此工程上正在尋求對發動機各部件進行減重設計。 薄壁加筋結構是航空航天裝備中的常見主承力構件, 進氣道也具有典型的薄壁特征。

加筋板殼結構因其良好的力學性能受到了越來越多的關注, 其結構特點是板殼結構和加強筋的結合。 目前針對該結構已有大量優化方法, 基結構法最先被提出用來解決各種加筋板殼結構優化問題3。 Zegard等4將基結構法引入到多種平面幾何結構進行拓撲優化。 牛斌等5采用基結構法作為主要工具, 將結構動柔度設定為優化的目標函數, 成功地完成了阻尼與加強筋布局的協同優化設計。 丁曉紅等6提出了自適應生長設計法, 該方法設定最小柔順度為優化目標來優化加強筋的布局和輪廓。 Li等7根據葉脈形成的原理, 提出一種仿生優化設計方法, 以結構應變能和第一階固有頻率的權重函數為目標函數, 對加筋板殼結構進行布局優化設計。 滕曉艷等8基于植物脈序分枝結構特性, 以諧振結構聲輻射功率為目標函數, 對薄板結構中加強筋布局優化設計。 Zhou等9基于漸近均勻化法和變密度法對筋條單胞優化設計, 提高加筋結構屈曲承載性能。 Feng等10基于變密度法和B樣條曲線參數化法, 以柔順度最小化為目標, 實現加強筋的布局和截面優化設計。 Huang等11提出了一個新的工程方法, 該方法包括了離散和連續設計變量, 并通過構建分支多點函數來近似優化問題, 使優化問題顯式, 從而實現了加筋板殼結構加強筋布局和尺寸優化設計。 目前, 對加筋板殼結構優化問題主要集中在對加強筋布局的優化設計, 考慮加強筋尺寸的優化設計還較少。

本文同時考慮加強筋布局優化與加強筋尺寸優化, 運用變密度法拓撲優化方法并結合點陣結構進行了優化設計, 在拓撲優化計算結果的基礎上對局部結構進行點陣加筋填充, 滿足進氣道變形和承載需求, 有效擴展了結構設計空間, 在保證原剛度和強度的基礎上進行進氣道的減重設計, 提高材料的利用率。

1 物理模型及優化方法

1.1 物理模型

圖1為典型進氣道結構示意圖," 由進氣道主體、" 上蒙皮、 下蒙皮、 導流板等組成, 其中上蒙皮、 下蒙皮布置在進氣道的上下表面, 導流板位于進氣道主體內部的流道中。

1.2 優化理論

(1) 單約束優化算法

本文采用Liu等12在傳統SIMP法基礎上發展的改進算法, 其表達式為

Ei(ρi)=Emin+ρpi(E0-Emin)(1)

式中: p為給定的懲罰因子; Ei(ρi)為映射單元的彈性模量; ρi為單元偽密度; E0為材料的彈性模量; Emin為給定剛度最小值, 通常取0.000 1, 避免剛度矩陣的奇異解問題。

基于SIMP材料懲罰模型13, 以最小化柔順度作為優化目標, 體分比為約束, 建立拓撲優化列式為

find: X={ρ1, ρ2, …, ρN}T

min: C=uTKudΩ

s.t. Vc-Vmax≤0

0≤ρi≤1(2)

式中: ρi(i=1, 2, …, N)為單元偽密度; Vc為設計域結構體積分數; Vmax為體積約束; C表示結構柔順度。

(2) 多約束優化算法

MMA算法14對拓撲優化問題有著良好的穩定性和適應性, 且在求解多約束問題有著更優秀的表現, 因此選擇數學規劃法中的MMA算法作為本文優化求解方法, 對式(2)增加約束, 具體形式為

find: X={x1, x2, …, xn}T

min: f0(X)

s.t. fi(X)≤fmaxi," i=1, 2, …, m

xminj≤xj≤xmaxj, j=1, 2, …, n(3)

式中: X為設計變量; xmaxj為設計變量上限; xminj為設計變量下限; fi(x)為設計變量x對應的約束值; fmaxi為上限約束; m為約束數量; f0(x)為目標函數。

基于上述公式, 初始化或更新單元設計變量, 將子問題解重復迭代, 直到滿足優化目標或者達到設置迭代步數后停止。

(3) 薄壁加筋結構的多約束

對于薄壁結構加筋的拓撲優化, 需要進一步的改進。

加筋結構需要兩方面的約束:

(a) 沖壓約束, 即變量單元的密度在薄壁結構的厚度方向上保持恒定。 由于厚度方向與單元的幾何信息密切相關, 為實現沖壓約束必須構建單元和單元之間的關系, 這對結構化網格是較為直觀的。 然而, 對于非結構化網格的結構, 單元的幾何信息和薄壁結構厚度方向沒有明顯的關系, 從而導致沖壓約束的可用性差。

(b) 尺寸約束, 拓撲優化結果往往使單元密度堆積在傳力路徑上, 盡管這種設計在理論上表現良好, 但是在實際生產過程中卻面臨諸多挑戰, 或者與預期的結構性能存在顯著差異。 為了獲得筋條形式的拓撲優化結果, 需要控制單元既不能生成孔洞, 也不能大量堆積, 因此需要孔洞尺寸約束以及最大尺寸約束。 最大尺寸約束的表達式為

∫Ωr(y)ρ(x)dΩ<∫Ωr(y)dΩ" (y∈Ω)(4)

式中: 等式左邊表示單元質心為半徑測試域Ω的總體積, 右側表示約束的總體積。 對于任意單元e, 其測試域的表達式為

x∈Ωer if x-x-e≤rmax(5)

孔洞尺寸約束的表達式為

∑i∈Revi(1-ρi)η≥Vemin(6)

式中: vi為單元體積; Vemin為測試域最小空隙體積約束; 參數η控制灰度單元被認定是空隙的程度。 通過孔洞尺寸約束以及最大尺寸約束, 可以保證拓撲優化結果在沖壓約束下構成加筋形式, 對后續復雜薄壁結構的加筋約束提供面內約束的技術基礎。

另外, 為了解決拓撲優化求解過程中可能遇到的棋盤格問題和網格依賴性求解問題等, 本文采用各向異性Helmholtz過濾方法15來求解算法方程, 避免了大量網格信息的獲取步驟, 可以更好地實現沖壓約束, 這對復雜薄壁結構的拓撲優化有重要優勢, 可以大大加快收斂進程, 并且對網格劃分的質量依賴更小。

密度過濾由濾波函數和密度的卷積實現, 表達式為

ρ~(x)=(Fρ)(x)=∫BRF(x-y)ρ(y)dy(7)

式中: F為線性卷積核; BR為以該網格中心點為中心半徑的圓形區域; ρ~為過濾后的密度, 即最后用于計算彈性模量的密度。 對于濾波函數有

∫BRF(x)dx=1(8)

卷積積分在離散情況下采用表達式為

ρ~j(x)=∑i∈Ne, jω(xi)viρi∑i∈Ne, jω(xi)vi(9)

式中: 權函數ω(xi)隨距離線性或指數增大。

在實際有限元計算的應用中, Helmholtz方程可寫為

ΔTKdΔρ~+ρ~=ρ(10)

式中: Kd為正定的張量, Kd=r2r2r2, r為過濾半徑尺度。

通過修改Kd的過濾半徑, 使沿薄壁結構厚度方向設置遠大于其他方向的過濾尺寸, 就可以實現在不借助相鄰坐標信息的前提上實現沖壓約束, 利用該特性, 建立復雜薄壁結構的加筋描述方法16

1.3 優化流程

進氣道優化流程如圖2所示。 首先, 確定進氣道的優化目標, 構造優化問題的目標函數與約束函數, 建立進氣道的有限元模型, 根據材料屬性、 約束邊界條件、 載荷邊界條件、 連接關系等信息對其進行有限元靜力校核, 根據計算結果確定優化設計區域, 而后進行優化設計域的拓撲優化計算工作。根據裝配關系同時參考有限元仿真結果, 從而確定優化設計區域及優化模型, 而后進行優化設計域的拓撲優化工作。 根據得到的拓撲優化結果重構幾何模型, 并對重構后的模型進行靜力校核, 比較優化前后結構的強度、 剛度與穩定性, 迭代優化, 直至滿足設計要求。

1.4 優化方案

將側壁設計區域簡化為500 mm×150 mm的平板作為驗證算例, 四邊簡支, 一面受均布壓強作用, 厚度方向劃分10層網格, 進行初始結構靜力校核分析; 屈曲分析主要用于研究結構在特定載荷下的穩定性以及確定結構失穩的臨界載荷, 以線性失穩特征值作為表征, 特征值越大說明結構越穩定。

首先計算正置正交加筋方案的結果," 設置筋條寬度為2 mm, 設計域體分比設置為0.3, 根據平板長寬比例設置正置正交加筋方案, 線性失穩特征值為10.040, 優化結果如圖3所示。

然后設置筋條寬度為2 mm、 域體分比為0.3的蜂窩加筋方案, 線性失穩特征值為-36.367, 計算結果如圖4所示。

下面進行拓撲優化加筋方案的對比。 使用1 mm基板加2 mm加筋板方案, 基板厚度設置3層網格, 加筋板厚度設置7層網格, 使用MMA求解器與各向異性Helmholtz過濾方法, 設置體分比為0.3, 設置不同的尺寸約束, 得到的優化結果如表1所示。

綜合分析剛度、 強度、 穩定性, 拓撲優化加筋方案具有較好的優化效果, 但線性失穩特征值較小, 且存在局部過薄的問題, 而蜂窩加筋方法的線性失穩特征值絕對值高達36.367, 具有良好的穩定性優化效果。 因此, 為了彌補薄壁面穩定性的不足, 在后續的優化過程中均采用主級加筋結合局部點陣加筋填充的填充方案。

2 進氣道結構拓撲優化

2.1 約束條件

進氣道主體與連接件、 蓋板、 蒙皮采用綁定連接, 與導流板則采用耦合連接。" 根據進氣道實際工作情況, 約束邊界條件分別定義進氣道頭部、 中部、 尾部連接件的位移約束," 如圖5所示, 對表面進行合理的固定約束, 對稱約束等。

2.2 載荷條件

進氣道工作承受最大載荷主要分為兩種, 一是轉級振蕩過程中的動態壓強載荷, 根據最大載荷位置確定了振蕩載荷的3種工況, 二是長時間工作在某一典型彈道的載荷, 確定了溫度壓強載荷。 進氣道工況如表2所示。

進氣道結構受溫度載荷、 壓強載荷的作用, 根據進氣道在流場仿真結果確定進氣道在各個坐標下的壓強場數據, 作為載荷的輸入條件, 從而確定進氣道結構校核的4種工況, 如圖6所示。

2.3 原始方案的強度計算結果

進氣道在振蕩載荷1作用下進氣道初始強度校核結果如圖7所示。

進氣道初始強度校核整體應力水平在材料屈服強度極限內, 局部應力超過855 MPa。 進氣道在4種工況下的計算結果如表3所示。 結果表明, 進氣道結構在4種工況下的應力水平、 位移水平和線性失穩特征值方面均滿足靜力強度要求, 振蕩載荷3和溫度壓強載荷兩種工況對進氣道結構更為嚴苛。

2.4 進氣道結構優化

拓撲優化計算過程靈敏度求解十分耗時, 設計域大小也直接影響著計算的時間, 同時考慮工程應用需求與工藝約束," 采用分開劃分設計域優化計算方案。 根據進氣道初始結構校核結果確定上流道兩側壁、上流道背壁、 下流道兩側壁、 下流道背壁、 后側流道背壁、 導流板、 底板為主要優化設計對象, 并開展筋條最大尺寸分別為5 mm, 10 mm, 15 mm," 20 mm," 25 mm, 30 mm的優化計算。 圖中空白部分為需要額外施加的筋條特征區域。

(1) 上流道兩側壁優化結果

上流道側壁優化結果如圖8所示。

分析結果發現, 尺寸約束增大到15 mm時優化結果較為清晰, 在兩面的交界處均有材料的保留以加強結構。 主側壁與副側壁的優化結果出現了局部過薄的問題, 優化時采用主級加筋與蜂窩加筋相結合的優化方案。

(2) 上流道背壁優化結果

上流道的優化結果如圖9所示。

上流道背壁的優化結果大體呈現出正置正交的筋條特征," 根據優化結果, 對筋條間距與寬度作部分調整作出正置正交的加筋優化設計方案。

(3) 下流道兩側壁優化結果

下流道兩側壁優化結果如圖10所示。

分析結果發現," 在兩面的交界處均有材料的保留以加強結構。 主側壁與副側壁的優化結果出現了局部過薄的問題, 優化時采用主級加筋與蜂窩加筋相結合的優化方案。

(4) 下流道背壁優化結果

下流道的優化結果如圖11所示。

下流道背壁的優化結果大體呈現出正置正交的筋條特征," 根據此優化結果, 對筋條間距與寬度作部分調整作出正置正交的加筋優化設計方案。

(5) 后側流道背壁優化結果

后側流道背壁優化結果如圖12所示," 選用主級加筋與點陣填充的優化設計方案。

(6) 導流板及其側壁優化結果

導流板優化結果如圖13所示, 整體呈現出對稱形式的筋條特征, 采用主級加筋優化設計方案。

(7) 底板優化結果

考慮到下蒙皮整體的承載特性, 考慮將下蒙皮與主體合并成底板進行優化。 優化結果如圖14所示," 采用主級加筋優化設計方案。

(8) 優化結束后的模型

上述給出了各個優化設計域在不同尺寸控制參數下的優化結果, 都經過50-70步迭代后趨于收斂, 結構邊界逐漸清晰、 平滑。 結合各優化結果的柔順度大小、 結構清晰度以及整體尺寸的把握, 本研究選用各部件在最大尺寸約束為15 mm約束下的優化結果。 重構模型, 得到優化后的模型如圖15所示, 采用主級加筋結合局部點陣加筋填充的優化方案, 按照各優化區域的優化設計方案完成重構設計, 以此模型進行優化后的靜力校核, 檢驗結構是否滿足設計要求。

3 優化后進氣道結構強度校核

首先將重構模型在振蕩載荷1下進行靜力校核, 得到結果如圖16所示。 可以看到, 在上、 下流道背壁筋條處出現了顯著的應力集中, 應力遠超材料的屈服強度。

為了對薄弱部位進行補充, 將筋條高度由2 mm增大到4.5 mm, 為減輕質量, 掏空筋條中段, 呈拱橋型, 并調整相鄰區域結構再次進行強度校核, 結果如圖17所示。 可以看到, 筋條的應力集中有了顯著改善。

針對此優化后模型開展4種載荷條件下有限元分析計算, 得到結果如表4所示。 優化模型在溫度穩壓載荷下, 應力水平基本不變, 最大位移略微增大2.1%, 整體保證了剛度分布; 在振蕩載荷作用下, 應力水平有所下降, 最大位移至少減小8.4%, 穩定性增強, 整體剛度分布更好; 優化后減重0.94 kg, 整體減重8%。

4 結" 論

(1) 本文以某進氣道為研究對象, 為實現進氣道薄壁結構減重設計, 同時考慮加強筋布局優化與加強筋尺寸優化, 運用變密度法拓撲優化方法并結合點陣結構方法進行了優化設計, 確定了主級加筋結合局部點陣加筋填充的優化方案。

(2) 利用有限元軟件對其進行優化前后強度、 剛度、 穩定性校核, 整體保證了進氣道結構的強度和剛度分布, 結構優化后進氣道結構質量減少了8%, 提高材料的利用率, 實現了對進氣道的減重設計。

參考文獻:

[1] 任淼," 劉晶晶," 文琳. 2023年國外空空導彈發展動態研究[J]. 航空兵器," 2024," 31(3): 32-39.

Ren Miao," Liu Jingjing," Wen Lin. Research on Foreign Air-to-Air Missiles’ Development in 2023[J]. Aero Weaponry," 2024," 31(3): 32-39. (in Chinese)

[2] 劉博宇," 王向明," 楊光," 等. 面向金屬增材制造的拓撲優化設計研究進展[J]. 中國激光," 2023," 50(12): 1202301.

Liu Boyu," Wang Xiangming," Yang Guang," et al. Research Progress on Topology Optimization Design for Metal Additive Manufacturing[J]. Chinese Journal of Lasers," 2023," 50(12): 1202301. (in Chinese)

[3] 霍澤凱," 王博," 周演," 等. 面向增材制造的雙蒙皮夾層結構加筋拓撲優化方法[J]. 計算力學學報," 2023," 40(6): 861-871.

Huo Zekai," Wang Bo," Zhou Yan," et al. Topology Optimization of Double-Skin Stiffened Sandwich Structure for Additive Manufacturing[J]. Chinese Journal of Computational Mechanics," 2023," 40(6): 861-871. (in Chinese)

[4] Zegard T," Paulino G H. GRAND-Ground Structure Based Topology Optimization for Arbitrary 2D Domains Using MATLAB[J]. Structural and Multidisciplinary Optimization," 2014," 50(5): 861-882.

[5] 牛斌," 閆家銘," 毛玉明," 等. 面向動力學性能的薄壁加筋板結構阻尼與筋條布局協同優化設計[J]. 中國機械工程," 2021," 32(16): 1912-1920.

Niu Bin," Yan Jiaming," Mao Yuming," et al. Collaborative Design Optimization of Damping Layers and Stiffener Layout of Thin-Walled Stiffened Plate Structures for Dynamics Performances[J]. China Mechanical Engineering," 2021," 32(16): 1912-1920. (in Chinese)

[6] 丁曉紅," 李國杰," 蔡戈堅," 等. 薄板結構的加強筋自適應成長設計法[J]. 中國機械工程," 2005," 16(12): 1057-1060.

Ding Xiaohong," Li Guojie," Cai Gejian," et al. Adaptive Growth Method of Rib Distribution for Thin Plate Structure[J]. China Mechanical Engineering," 2005," 16(12): 1057-1060. (in Chinese)

[7] Li B T," Hong J," Liu Z F. Stiffness Design of Machine Tool Structures by a Biologically Inspired Topology Optimization Method[J]. International Journal of Machine Tools and Manufacture," 2014," 84: 33-44.

[8] 滕曉艷," 江旭東," 鄒廣平," 等. 簡諧激勵下板殼結構加強筋仿生布局降噪方法[J]. 振動與沖擊," 2016," 35(6): 1-6.

Teng Xiaoyan," Jiang Xudong," Zou Guangping," et al. Bionic Approach for Stiffener Layout on Plate and Shell Structure under Harmonic Excitation for Noise Reduction[J]. Journal of Vibration and Shock," 2016," 35(6): 1-6. (in Chinese)

[9] Zhou Y," Tian K," Xu S L," et al. Two-Scale Buckling Topology Optimization for Grid-Stiffened Cylindrical Shells[J]. Thin-Walled Structures," 2020," 151: 106725.

[10] Feng S Q," Zhang W H," Meng L," et al. Stiffener Layout Optimization of Shell Structures with B-Spline Parameterization Method[J]. Structural and Multidisciplinary Optimization," 2021," 63(6): 2637-2651.

[11] Huang H," An H C," Ma H B," et al. An Engineering Method for Complex Structural Optimization Involving both Size and Topology Design Variables[J]. International Journal for Numerical Methods in Engineering," 2019," 117(3): 291-315.

[12] Liu K," Tovar A. An Efficient 3D Topology Optimization Code Written in Matlab[J]. Structural and Multidisciplinary Optimization," 2014,nbsp; 50(6): 1175-1196.

[13] Andreassen E," Clausen A," Schevenels M," et al. Efficient Topology Optimization in MATLAB Using 88 Lines of Code[J]. Structural and Multidisciplinary Optimization," 2011," 43(1): 1-16.

[14] Svanberg K. The Method of Moving Asymptotes-A New Method for Structural Optimization[J]. International Journal for Numerical Methods in Engineering," 1987," 24(2): 359-373.

[15] 王博," 周子童," 周演," 等. 薄壁結構多層級并發加筋拓撲優化研究[J]. 計算力學學報," 2021," 38(4): 487-497.

Wang Bo," Zhou Zitong," Zhou Yan," et al. Concurrent Topology Optimization of Hierarchical Stiffened Thin-Walled Structures[J]. Chinese Journal of Computational Mechanics," 2021," 38(4): 487-497. (in Chinese)

[16] Wang B," Zhou Y," Tian K," et al. Novel Implementation of Extrusion Constraint in Topology Optimization by Helmholtz-Type Anisotropic Filter[J]. Structural and Multidisciplinary Optimization," 2020," 62(4): 2091-2100.

Numerical Investigation of Inlet Topology Optimization

for Additive Manufacturing

Guo Jiaqi1," Sun Zhenhua1," 2*," Liu Yadong1," Bai Taotao1

(1.China Airborne Missile Academy," Luoyang 471009," China;

2.National Key Laboratory of Air-based Information Perception and Fusion," Luoyang 471009," China)

Abstract: The inlet is one of the important parts of ramjet," and its strength and stiffness are very important to the stability of ramjet. In order to reduce the negative mass of the ramjet and optimize the ramjet design," this paper takes an intake port as the research object," and carries out the weight reduction design of the thin-wall structure of the intake port on the basis of meeting the requirements of the intake port deformation and load. Static check is carried out on the inlet," and the optimal design domain is determined according to the results of the initial strength check of the inlet under four working conditions. Combined with the lattice structure and considering the optimization of reinforcement layout and reinforcement size," the inlet is optimized by using the variable density topology optimization method," and the optimization scheme of main reinforcement and local reinforcement filling is adopted. Finite element software is used to check the strength," stiffness and stability before and after optimization. After optimization," the structural mass of the inlet is reduced by 8%," the stability is enhanced," the overall stiffness distribution is better," the weight reduction of the inlet and the expansion of the structural design space are realized," and the feasibility of the optimization method adopted in this paper is proved.

Key words:" topology optimization; numerical investigation; inlet; lattice structure; thin stiffened structur;" ramjet

主站蜘蛛池模板: 91在线无码精品秘九色APP| 亚洲第一色视频| 国产黄网站在线观看| 日本国产精品| 国产精品任我爽爆在线播放6080 | 中文字幕亚洲乱码熟女1区2区| 人妻出轨无码中文一区二区| 久久成人国产精品免费软件| 欧美不卡视频一区发布| 国产免费a级片| 亚洲AⅤ无码国产精品| 国产性生交xxxxx免费| 国产女人水多毛片18| 日本国产一区在线观看| 欧美一级夜夜爽| 亚洲婷婷在线视频| 欧美在线国产| 国产成人精品免费视频大全五级| 成年女人a毛片免费视频| 美女无遮挡拍拍拍免费视频| 尤物亚洲最大AV无码网站| 看你懂的巨臀中文字幕一区二区 | 无码日韩精品91超碰| 国产中文一区a级毛片视频| 亚洲第一综合天堂另类专| 亚洲综合亚洲国产尤物| 国产精品尹人在线观看| 波多野结衣爽到高潮漏水大喷| 就去色综合| 欧美午夜理伦三级在线观看 | 露脸一二三区国语对白| h网址在线观看| 免费人成在线观看成人片| 最新加勒比隔壁人妻| 99久久99这里只有免费的精品| 91小视频在线观看免费版高清| 国产免费网址| 超碰免费91| 国产色伊人| 国产成人高清亚洲一区久久| 亚洲乱码精品久久久久..| 一边摸一边做爽的视频17国产| 国产清纯在线一区二区WWW| 久久精品国产91久久综合麻豆自制| 永久天堂网Av| 精品国产aⅴ一区二区三区| 亚洲天堂区| 国产真实乱子伦精品视手机观看| av色爱 天堂网| 国产国语一级毛片| 天堂网亚洲综合在线| 亚洲天堂免费观看| 国产精品七七在线播放| 亚洲人成网站在线播放2019| 色有码无码视频| 欧美另类精品一区二区三区| 五月婷婷综合色| 欧美日韩国产精品va| 国产成人精品优优av| 国产chinese男男gay视频网| 色135综合网| 欧美曰批视频免费播放免费| 91 九色视频丝袜| 97se亚洲综合在线天天| 欧美激情第一欧美在线| 试看120秒男女啪啪免费| 伊人久久久久久久久久| 国产精品自在在线午夜区app| 国产凹凸视频在线观看| 白浆视频在线观看| 91精品亚洲| 国产精欧美一区二区三区| 日韩av电影一区二区三区四区| 色综合久久88| 国产综合网站| 日韩欧美中文在线| 久久久久久尹人网香蕉| 欧美精品亚洲日韩a| 免费精品一区二区h| 亚洲美女一级毛片| 亚洲不卡影院| 国产精品无码在线看|