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

面向滿應力準則的非線性連續結構截面尺寸優化設計理論與分析

2021-10-11 09:50:02何浩祥周鈺婧
振動與沖擊 2021年18期
關鍵詞:優化結構設計

吳 山,何浩祥,周鈺婧

(北京工業大學 工程抗震與結構診治北京市重點實驗室,北京 100124)

結構優化設計是指在給定的約束條件下,按某種優化目標(如重量最輕、成本最低、剛度最大等)獲得最優方案的設計方法。根據設計變量的類型,工程結構優化可以分為不同的層次[1]:①尺寸優化,在給定結構類型、材料、布局和幾何尺寸的情況下,確定各構件的最優截面尺寸,從而實現結構優化目標,是難度較低、應用最廣泛的優化方式[2];②形狀優化,在限定物理參數范圍的前提下,將構件的截面形狀、節點的空間位置以及連續體的形狀等作為設計變量,實現設計區域幾何形狀的改變,從而優化結構力學性能;③拓撲優化,以材料在設計域內的空間分布作為設計變量,使結構在指定條件下滿足性能需求或優化準則。結構優化求解的方法大致包括準則法、數學規劃法及智能算法等。準則法是從工程需求出發提出結構優化應滿足的力學準則(如滿應力準則、剛度準則和能量準則等)并迭代求解的方法,是最基本的優化求解方法。

大量研究表明:各層次結構優化的重要性是沒有絕對的高低之分的,應根據工程需求、優化目標、優化準則和計算效率從各層次開展協同研究才能真正實現結構優化設計。近年來,關于連續體拓撲優化的研究比較深入和豐富,并取得了一定應用。但對于實際結構而言,區域邊界需要保持光滑的連續性條件,優化的形狀邊界也需具有足夠的可制造性,而拓撲優化過于關注材料分布卻并未充分考慮實際條件的限定,加之考慮的荷載工況較單調,故其工程可行性偏弱。因此,基于工程結構的力學特征、普遍需求和優化準則,探究連續體的截面及形狀優化的理論形式和力學解析解具有重要的理論和工程價值。

此外,在諸多優化準則中,滿應力設計準則是最受關注的,其主要思想是在指定受力模式下使結構的每個構件或部分均達到承載力極限或設計容許值,從而使構件的強度得到充分利用[3]。在傳統的結構截面尺寸和形狀優化研究中,應優先考慮滿足滿應力設計準則,以避免出現因應力集中導致的局部變形過大或薄弱層破壞。在目前的研究中,通常采用容許應力法或彈性狀態假定認為結構始終處于彈性狀態而進行優化設計,在此條件下滿應力設計準則等效為等強度準則或等應力準則。關于滿足等應力準則的構件或簡單結構的優化設計研究取得了大量成果。例如,在簡單荷載下的矩形截面懸臂梁和簡支梁的等應力優化已經較成熟,基于多種拓撲演化策略的平面框架、桁架或簡單塔架的等應力優化也取得較令人滿意的效果[4]。

需要指出的是,由于工程結構通常會承受地震、風荷載和其他隨機荷載的作用,且當激勵幅值較大時結構的材料很可能進入彈塑性狀態并產生性能退變和較大變形,因此這些作用成為滿足結構設計需求的決定性荷載模式,相比簡單荷載情況,其優化設計研究更有理論意義。此外,若結構在接近臨界破壞彈塑性應力時也較均勻,其各處的損傷也將較接近,因此面向滿應力準則的結構優化結果也可為面向損傷的結構抗振設計提供初始優化結果,具有良好的工程意義。然而,目前在力學分析和優化領域,為了簡化分析通常將材料本構假定為線性形式,不考慮剛度退化等非線性特性,即并沒有真正實現滿應力優化設計。因此,對復雜劇烈荷載下的非線性結構滿應力優化設計是十分必要的。

綜上所述,面向實際荷載及材料特性的工程結構優化設計尚主要存在以下不足:①等應力和等強度設計理念通常被用于簡單荷載下的構件優化層面,對復雜荷載下的結構整體優化設計的探究和應用不夠深入和全面;②缺乏基于滿應力理念的連續懸臂梁結構最優截面尺寸分布的力學機理研究,從而無法為復雜荷載下結構優化設計提供堅實準確的精細化解析方法和公式;③缺乏對材料非線性結構的滿應力設計研究,不能滿足實際工程結構的彈塑性損傷設計需求和實用性需求。

有鑒于此,本文將諸如連續塔式結構和高層建筑結構等工程結構簡化為非線性連續變截面懸臂桿結構,通過構建不同的彎矩表達形式來求得最優截面尺寸解析解,采用通過試驗數據擬合獲得的二次多項式表征材料非線性本構,分別對圓形、矩形、圓環形和箱形截面懸臂桿在均布荷載、倒三角荷載和與慣性力相關荷載作用下的最優截面尺寸分布公式進行推導求解,得到完全解析解及最優剛度分布解析表達,之后通過有限元分析對以上理論方法和結果進行了充分驗證。

必須指出的是,由于實際高層建筑結構構件種類繁多,結構形式復雜,僅從力學優化角度難以獲得最優抗側剛度分布解析解,需要迭代算法對其進行深度優化設計,故本文基于力學理論對簡化模型進行理論分析,得到了最優截面尺寸分布和最優抗側剛度分布解析解,該結果主要可用于對諸如塔式結構等連續懸臂結構的優化設計,同時可為工程結構優化設計提供初始優化解,提高優化效率。

1 不同荷載形式下非線性結構最優截面尺寸分布解析

1.1 結構截面形式、荷載分布模式與材料本構

工程結構在外部激勵下發生的位移和損傷通常與結構的剛度分布和質量分布等因素密切相關。經過合理設計的結構剛度通常是沿高度方向連續變化的。連續體等應力優化理論主要適用于水塔、煙囪等構筑物和塔式結構,對于高層和超高層建筑,其樓板約束作用影響較小,其結構在水平荷載下以彎曲變形為主,因此也基本適用[5]。根據高層建筑物或構筑物不同的結構形式和工程需求,等效懸臂桿截面形式可包括圓環形、圓形、箱形和矩形等,其中圓形和矩形截面的最優尺寸分布公式可以由圓環形和箱形截面的公式表示,因此本文主要對圓環形和箱形截面的情況進行解析。在最優截面尺寸分布的基礎上,提出最優剪切剛度和彎曲剛度分布,對實際工程設計和規范修訂提供參考。

工程結構在不同荷載形式下的最優截面尺寸分布和最優剛度分布是不同的。風荷載和地震作用是使建筑結構發生破壞或損傷的常見不利作用。鐘振宇等[6]考慮風重耦合效應,推導了沿高度方向的等效靜力風荷載分布,結果表明等效靜力風荷載的分布形式非常接近倒三角形;當僅考慮第一振型時,地震作用可以等效為倒三角水平荷載;當考慮全部振型時,可將地震作用簡化為均布荷載或與慣性力相關荷載[7]。為了全面考慮風荷載和地震作用,本文考慮均布荷載、倒三角荷載和與慣性力相關荷載三種水平荷載形式。

以滿應力為目標的優化設計不僅要求結構在彈性階段實現等應力,更重要的是結構在塑性狀態下的應力相等或相近。目前相關領域的研究為了計算簡便,材料本構多假定為線性形式,不考慮剛度退化等非線性力學特性,實際工程多為鋼-混凝土結構或鋼結構,采用線性材料的研究往往難以應用于實際。為了確保研究在工程中的適用性和可行性,懸臂桿采用均質非線性材料。

1.2 優化問題描述

在地震動和風荷載作用下,對于高層建筑結構,在各層應力水平不超過限值的情況下,希望其應力沿結構高度方向的分布盡量均勻。通過對結構截面尺寸優化設計,以達到上述優化目標。該優化問題可以表述為

式中:Si為結構第i層截面尺寸;σi為結構第i層的應力值;[σ]為結構應力允許值;n為結構層數;t為時間;T為外荷載持時。

1.3 圓環形變截面懸臂桿最優截面尺寸分布

1.3.1 均布荷載下懸臂桿最優截面尺寸分布

大量試驗結果表明結構中的鋼筋混凝土可以等效為一種拉壓同性的復合材料[8],其應力-應變關系可以通過非線性函數表征。設該材料本構關系為

式中,σ和ε為材料受壓產生的應力和應變。受壓應力-應變全曲線方程可采用Hognestad建議的混凝土多項式表達式,其對鋼筋混凝土等復合材料也具有較好的通用性。假定等效復合材料為各向同性,即受拉和受壓時有相同的應力-應變全曲線,有

式中,a和b為待定系數。該表達式代表的材料應力-應變全曲線包括上升段和下降段,是材料真實力學性能的全面宏觀反應。

首先對均布荷載作用下圓環形變截面梁的最優截面尺寸分布進行求解,將結構等效為懸臂桿,桿件采用均質非線性材料,受壓應力-應變全曲線由式(3)表示。計算模型如圖1所示,取懸臂桿軸線為x軸,橫截面對稱軸為y軸,中性軸為z軸。

圖1 力學模型示意圖Fig.1 Mechanical model

首先,通過任意截面下由外荷載形成的彎矩與由梁截面尺寸與彎曲應變表示的彎矩相等的原則,建立彎矩平衡方程。以此原則優化后的懸臂桿結構在水平荷載作用下沿高度方向彎曲應變相等,因結構產生的應變以彎曲應變為主,故近似得Mises應變也是相等的,后文有限元驗證提取了結構的Mises應變,對該方法的有效性進行了驗證。相應的彎矩平衡方程為[9-10]

式中:M1為由梁截面尺寸與應變表示的彎矩公式;M2(x)為由外荷載表示的彎矩公式。對M1和M2(x)分別進行推導。

懸臂桿結構在小變形時,服從平面變形假定

式中:y為縱軸坐標;ρ為曲率半徑;d為圓形橫截面外直徑;d′為內直徑;設 σ1,σ2,ε1,ε2和 σ′1,σ′2,ε′1,ε′2分別為橫截面內外壁距離中性軸最遠兩處的應力和應變,即截面中外壁和內壁的最大應力和應變,則橫截面彎矩可寫成

式中,A為截面微面積。由式(5)及材料拉壓性質對稱[11]可知,式(6) 可寫為

對ε進行換元,設

則式(7)可寫為

當α=0時,截面為圓形。對式(9)和式(10)整理可得

當懸臂桿結構承受水平荷載作用時,在任意高度x處由均布荷載q引起的彎矩M2(x)為[12]

式中,H為結構總高度。將式(11)、式(12) 代入式(4)

式(13)為圓環形截面懸臂桿結構在均布荷載作用下最優截面尺寸分布解析解,截面尺寸滿足式(13)的懸臂桿結構在均布荷載作用下,各截面內最大應力和應變相等。

1.3.2 倒三角荷載下懸臂桿最優截面尺寸分布

倒三角荷載頂部的荷載集度為qmax,底部為0,當作用于懸臂桿結構時,荷載產生的彎矩為[13]

圓環形截面懸臂桿截面內彎矩為式(11),再通過與式(14)相等建立彎矩平衡方程,可以解得最優截面直徑分布為

1.3.3 與慣性力相關荷載下懸臂桿最優截面尺寸分布

當圓環形截面變截面懸臂桿結構承受與慣性力相關荷載時,荷載集度可表示為

式中:ρ為等效面密度;λ為質量相關系數,用來控制荷載幅值。由于與慣性力相關荷載集度與截面面積有關,而截面尺寸尚且未知,為了求解荷載作用產生的彎矩,進行如下假設

結構任一高度x截面內由于與慣性力相關荷載作用產生的彎矩可以表示為

圓環形截面內的彎矩可以表示為(11),與式(18)相等,可得

易知式(19) 為y″=f(x,y) 型高階微分方程,屬不可解類型,故考慮采用假設截面尺寸函數的方法進行求解。

由文獻[14]可知,懸臂桿截面形式按照下式分布有利于實現在與慣性相關荷載作用下的等位移原則,其高度分布函數為

圖2 與慣性力相關荷載下的懸臂桿模型Fig.2 Cantilever bar model under inertial force-related load

圓環形截面懸臂桿在與慣性力相關荷載的作用下,任意高度x截面內的彎矩為

圓環形截面的彎矩可以表示為

令式(22) 與式(21) 相等,得

因此,圓環形截面懸臂桿在與慣性力相關荷載作用下最優截面直徑分布為

圓環形截面懸臂桿截面尺寸按照式(25)分布,可以實現與慣性力相關荷載作用下的等應力原則。

1.4 箱形變截面懸臂桿最優截面尺寸分布

箱形最優截面尺寸分布和圓環形截面采用相同的原理,過程不再詳細描述。箱形外截面寬度B保持不變,對外截面高度h進行優化,設h(x)為截面高度隨結構高度變化的函數并求解析解,且有

式中:B′為箱形內截面寬度;h′為箱形內截面高度。

均布荷載下箱形截面結構最優截面尺寸分布函數為

倒三角荷載下,箱形截面結構最優截面尺寸分布函數為

與慣性力相關荷載作用下,箱形截面結構最優截面尺寸分布為

2 基于有限元的驗證分析

2.1 非線性材料本構

為了驗證本文所建立的理論方法和解析結果的正確性,需確定具體的材料非線性本構模型取值并進行彈塑性有限元分析。對于式(3)中的非線性材料本構參數,可根據試驗數據進行擬合確定。本文選擇文獻[15]中大尺寸鋼筋混凝土軸心受壓試驗的數據進行參數確定。試驗中S30-400試件長1 800 mm,正方形截面尺寸400 mm×400 mm,長細比4.5,混凝土強度設計等級為C30,實測棱柱體抗壓強度22.3 MPa,縱筋采用HRB335鋼筋,實測屈服強度408 MPa,縱筋等效配筋率2.51%,加載制度為擬靜力加載,對S30-400試件的應力-應變曲線進行二次多項式擬合,從而確定式(3)中的待定系數a和b,擬合效果如圖3所示。式(3)擬合結果為

圖3 材料本構曲線擬合效果Fig.3 Fitting effect of material constitutive curve

式(31)可視為普通鋼筋混凝土結構采用的等效復合材料受壓應力-應變全曲線表達式。

2.2 靜力驗證分析

上文對圓環形、箱形截面非線性懸臂桿在均布荷載、倒三角荷載和與慣性力相關荷載作用下的最優截面尺寸公式進行了推導求解,并得到解析解。但公式在結構處于彈性階段及彈塑性階段下的正確性及精度有待進一步驗證。利用ANSYS有限元平臺對截面尺寸按照上文公式分布的變截面非線性懸臂桿結構進行靜力分析。

懸臂桿結構的基本設計參數為:結構總高度H為20 m,非線性材料初始彈性模量25.94 GPa,泊松比0.3,有限元模型沿結構高度方向均勻設置了21個節點,懸臂桿底面與地面固接,每個節點處布置一個截面,截面尺寸按照上節公式的結果設置,采用Beam188梁單元,該單元適合模擬Timoshenko梁,考慮剪切變形的影響,每個相鄰節點之間設置taper變截面梁單元。按照求得的解析解,懸臂桿頂部截面面積應為0,然而ANSYS中梁單元不允許出現面積為0的截面,因此設頂部截面面積略大于0。材料本構采用KINH多線性隨動強化模型,考慮包辛格效應,設置了11個數據點,包含了上升段、峰值點和下降段,以便對結構產生不同應變的情況進行模擬,驗證公式的適用性和準確性。

對圓環形、箱形截面懸臂桿在均布荷載、倒三角荷載和與慣性力相關荷載作用下的結構響應進行分析,分為截面最大正應變分別為0.001,0.002和0.003的情況討論,一共建立了7個懸臂桿結構模型,設計參數和工況如表1所示。

表1 懸臂桿結構的設計參數和工況Tab.1 Design parameters and working conditions of cantilever structures

由于材料本構中應力與應變一一對應,因此可以用應變分布表征應力分布,同時為了對比設計正應變g(σ1)與實際應變,故提取加載后有限元模型的Mises應變云圖,如圖4所示。由結果可見:若按相應的理論公式設置截面尺寸,變截面連續懸臂桿結構的確能夠在各截面均能實現滿應力分布,并且適用于結構的彈性和塑性階段,驗證了公式中關鍵參數的正確性。應變云圖中結構頂部出現一定程度的應力集中,原因是在模型頂部設置了一個略大于0的面積后,該局部的實際值和理論值產生了誤差。

圖4 結構應變云圖Fig.4 Strain nephogram of specimens

提取加載后結構全部單元的最大應變,如圖5所示,所有構件應變云圖趨勢類似,高度15 m以下部分的應變均小于設計正應變g(σ1),這是由于有限元中提取的是Mises應變,結構橫截面正應變是影響Mises應變的主要因素,因此Mises應變雖然分布均勻但在幅值上與正應變存在差異,此外,taper變截面梁單元上下截面之間存在線性變化,這與最優截面尺寸分布解析解之間存在一定誤差。

圖5 結構應變分布Fig.5 Strain distribution of specimens

為了驗證優化設計效果,在最大應變為0.001的情況下,將均布荷載作用下未經過優化的等截面懸臂桿結構與優化過的變截面懸臂桿結構HU1進行對比,等截面結構底面直徑8 m,HU1底面直徑9.38 m,α均取0.9。有限元模型參數除截面變化規律外與優化過的結構保持相同。另外,由于HU1出現頂部應力集中現象,設置一個頂部截面增大的構件,該構件頂部圓形截面直徑0.085 m,HU1頂部截面直徑0.05 m,其余參數相同,三者對比結果如圖6所示。

圖6 未優化結構與優化結構對比結果Fig.6 Comparison of optimized structure

從圖6中可以看出,未經過優化的等截面懸臂桿結構應變呈明顯的非均勻分布,底部應變大,隨著結構高度增加,逐漸趨近于零,而優化過的變截面懸臂桿HU1雖然從15 m高度往上應變逐漸非線性增大,頂部出現一定的應力集中,但15 m以下應變基本相等,并且頂部應力集中現象可以通過增大頂部截面面積的方法解決,等截面結構體積為100.48 m3,應變分布均方差為4.5×10-4,頂部加強結構體積為85.56 m3,應變分布標準差為3.07×10-6,優化后結構體積減少了14.8%,應變分布明顯更均勻,充分證明了該優化方法的優越性。

需要指出的是,在均布荷載和倒三角荷載作用下可以實現滿應力準則的結構外觀呈“外凸”型,而在與慣性力相關荷載作用下則呈“內凹”型,兩者截然不同,懸臂桿結構模型如圖7所示。在實際工程中,均布荷載(倒三角荷載)和與慣性力相關荷載可能共同作用于結構,可以根據不同的工程需求對兩者賦予不同權重,如圓環形結構最優截面尺寸分布按式(32)計算

式中,?W和?E分別為倒三角荷載和與慣性力相關荷載的權重系數,?W+?E=1。當?W=?E=0.5時,懸臂桿模型如圖7(c)所示。

圖7 變截面懸臂桿有限元模型Fig.7 Finite element model of variable section cantilever

2.3 動力驗證分析

前述研究是在靜力學基礎上對地震動和風荷載作用下結構的截面尺寸分布進行優化。為了驗證優化結果在動力作用下的準確性,對前文中頂部加強的懸臂桿結構進行動力時程分析。將優化后的圓環形和箱形截面懸臂結構和未優化過的等截面懸臂結構進行對比,對所有結構基底輸入單向El Centro地震波,時間間隔0.02 s,持時15 s。由于在2.14 s時地震波加速度時程達到峰值,5.50 s時位移時程達到峰值,因此提取以上兩個時刻的結構各節點應變,結果如圖8和圖9所示。分析結果可發現:在經過頂部加強后,優化后圓形和箱形截面懸臂桿在地震作用下的應變分布均勻,均實現了滿應力準則;而未優化的等截面懸臂結構的應變分布則從底部到頂部非線性減小,應變分布離散性遠大于優化后的懸臂桿結構。經進一步分析結果可確認優化結構的應變分布規律在任何時刻均相似,充分實現了滿應力準則。綜上,按與慣性力相關荷載進行優化的結構,在動力作用下也實現了滿應力準則,且不同截面的優化結果均具有良好效果。

圖8 地震動加速度最大峰值時結構應變分布Fig.7 Strain distribution under peak ground acceleration

圖9 震動位移最大峰值時結構應變分布Fig.9 Strain distribution under peak ground displacement

3 結構最優剛度分析

上文推導出了均布荷載、倒三角荷載和與慣性力相關荷載作用下圓環形、箱形截面懸臂桿結構的最優截面尺寸分布公式,并基于ANSYS有限元平臺驗證了按照該公式布置截面的結構可以在相應工況下實現滿應力準則。其剪切剛度和彎曲剛度計算公式如下

按照式(33)和式(34)計算各種工況下公式對應的剛度,此結果可供實際工程設計和規范修訂參考。

均布荷載作用下,圓環形截面懸臂桿剪切剛度和彎曲剛度為

倒三角荷載作用下,圓環形截面懸臂桿剪切剛度和彎曲剛度為

與慣性力相關荷載作用下,圓環形截面懸臂桿剪切剛度和彎曲剛度為

箱形截面懸臂桿結構的剪切剛度和彎曲剛度同理,不再贅述。

本文的理論方法和優化結果適用于諸如塔式結構和煙囪等連續體的靜動力優化設計,也較適用于高細比較大的高層建筑等工程結構,相應的優化設計不必拘泥于截面固有形式,可根據結構截面尺寸或剛度的優化解析解進行適當調整,并可進一步通過調整抗側力構件配筋和局部桿件截面尺寸等方式靈活調整結構抗側剛度分布,使結構的等效參數與理論優化結果相近,實現滿應力準則。與拓撲優化相比,本文方法不以材料取舍為手段,相應結果的區域邊界滿足連續性條件,更符合實際條件需求和限制,具有更好的工程可行性。此外,與何浩祥等以相對位移相等為優化目標的方法相比,本文結果與之相近但有所差別,且從材料應力水平真正實現了滿應力分布,因此具有更嚴密的力學依據。

4 結 論

將工程結構簡化為圓環形、箱形截面懸臂桿結構,以等應力為優化目標進行優化設計。將風荷載和地震作用等激勵簡化為均布荷載、倒三角荷載和與慣性力相關荷載,建立彎矩平衡方程求解,得到最優截面尺寸解析解,并獲得最優剪切剛度和彎曲剛度分布。通過有限元模擬驗證了解析解的正確性和適用性。得到以下結論:

(1)在彈性和彈塑性階段,在指定荷載模式下的懸臂桿截面尺寸沿高度按本文獲得的理論分布函數進行設置可以實現滿應力。

(2)均布荷載和倒三角荷載作用下,按等應力理念優化設計的懸臂桿結構呈現“外凸”形狀,而與慣性力相關荷載作用下的結構呈現“內凹”形狀,兩者截然相反,實際工程中,結構往往受到多種荷載同時作用,可以根據當地的工程需求對不同荷載的結果賦予不同權重進而獲得組合荷載下的優化結果。

(3)基于ANSYS有限元平臺,對各種工況下的懸臂桿最優截面尺寸公式分別進行了靜力和動力優化驗證。在靜荷載作用下,按等應力優化設計的懸臂桿結構在彈性和彈塑性狀態下除頂部出現應力集中現象外,其余部分均處于等應力狀態,頂部應力集中問題可以通過增大頂部截面的方法解決。頂部經過加強后,按與慣性力相關荷載作用下優化結果設計的結構處于等應力狀態,因此優化結果在動力荷載作用下同樣有良好的效果。

猜你喜歡
優化結構設計
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
瞞天過海——仿生設計萌到家
藝術啟蒙(2018年7期)2018-08-23 09:14:18
設計秀
海峽姐妹(2017年7期)2017-07-31 19:08:17
有種設計叫而專
Coco薇(2017年5期)2017-06-05 08:53:16
論《日出》的結構
主站蜘蛛池模板: 2021国产精品自产拍在线观看| 免费不卡视频| 男女性色大片免费网站| 久久精品电影| 成人国产精品2021| 日本影院一区| 国产成人艳妇AA视频在线| 亚洲无码久久久久| 视频一区视频二区中文精品| 免费Aⅴ片在线观看蜜芽Tⅴ| 99精品在线看| 亚洲人成网7777777国产| 久久精品国产在热久久2019| 九九久久精品免费观看| 伊人色天堂| av色爱 天堂网| 国产成人高精品免费视频| 亚洲国产理论片在线播放| 日本午夜视频在线观看| 亚洲一级色| 99人体免费视频| 国产精品网址在线观看你懂的| 色丁丁毛片在线观看| 久久天天躁狠狠躁夜夜躁| 亚洲中文字幕97久久精品少妇| 欧美精品不卡| 亚洲激情99| 欧美精品v| 久久一本精品久久久ー99| 久久婷婷六月| 欧美亚洲国产精品久久蜜芽| 亚洲天堂视频在线观看免费| 91口爆吞精国产对白第三集| 亚洲熟女中文字幕男人总站| 国产日本欧美亚洲精品视| 亚洲丝袜第一页| 中国毛片网| 欧美国产综合色视频| 青青草久久伊人| 日韩精品无码不卡无码| 免费不卡视频| 日韩欧美国产另类| 国产九九精品视频| 免费观看亚洲人成网站| 欧美激情视频二区| 亚洲国产天堂久久综合226114| 97人人做人人爽香蕉精品| 亚洲一区二区三区在线视频| 亚洲第七页| 婷婷六月综合网| 永久免费av网站可以直接看的| 免费欧美一级| 国产成人精品一区二区秒拍1o| 亚洲三级电影在线播放| 欧美日韩久久综合| 拍国产真实乱人偷精品| 精品一区二区三区四区五区| 日韩AV无码免费一二三区| 国产精品理论片| 日本人妻丰满熟妇区| 日韩少妇激情一区二区| 免费一级全黄少妇性色生活片| 久久久无码人妻精品无码| 国产成人1024精品下载| 日韩精品专区免费无码aⅴ| 日韩在线观看网站| 99国产精品国产| 亚洲大尺码专区影院| 视频二区国产精品职场同事| 成人va亚洲va欧美天堂| 国产农村妇女精品一二区| 国产在线观看高清不卡| 无码人妻免费| 国产鲁鲁视频在线观看| 成人毛片免费在线观看| 精品国产成人高清在线| 国产精品免费入口视频| 国产乱人乱偷精品视频a人人澡| 亚洲成人精品| 欧美日韩高清| 亚洲精品第一在线观看视频| 99视频精品在线观看|