亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        結(jié)構(gòu)拓?fù)鋬?yōu)化局部性能約束下輕量化問題的互逆規(guī)劃解法

        2021-09-07 08:34:22隋允康彭細(xì)榮葉紅玲李宗翰
        計算力學(xué)學(xué)報 2021年4期
        關(guān)鍵詞:優(yōu)化結(jié)構(gòu)方法

        隋允康, 彭細(xì)榮*2, 葉紅玲, 李宗翰

        (1.北京工業(yè)大學(xué) 材料與制造學(xué)部,北京 100022; 2.湖南城市學(xué)院 土木工程學(xué)院,益陽 413000)

        1 引 言

        錢令希[1]指出力學(xué)性能約束可劃分為局部性約束和整體性約束,這是就截面層次優(yōu)化的劃分,同樣也適用結(jié)構(gòu)拓?fù)鋬?yōu)化層次[2-5]。

        力學(xué)性能局部性與整體性的區(qū)別直接影響建模與求解的差別。局部性力學(xué)性能約束下的結(jié)構(gòu)優(yōu)化問題中,出現(xiàn)了建模困難。因約束太多,其個數(shù)是單元或子域的個數(shù)與載荷工況的乘積。不管應(yīng)力和局部穩(wěn)定還是疲勞壽命約束,相關(guān)性能在每個單元或子域上都必須滿足。而整體性力學(xué)性能約束,如位移或頻率約束,約束數(shù)目相比局部性能約束少很多,相應(yīng)的敏度計算量也因為各種處理如采用伴隨法而變得比局部性能約束少很多。

        局部性能與整體性能的本質(zhì)區(qū)別在于單元對于性能的作用不同,局部性約束要求在所有單元里相關(guān)的性能都不能違背,故稱為單元滿足型;整體性約束只要求結(jié)構(gòu)總體不能違背相關(guān)的性能,故稱為性能的單元貢獻(xiàn)型。為了克服局部性能約束優(yōu)化問題的困難,以往研究已提出過分部[1,2-5]、化整[3-5]及集成[2,3,9-13]3種解法,分別簡稱為P解法(Partition approach)、G解法(Globalization approach)和I解法(Integration approach)。本文將闡述與之相關(guān)新提出的3種交融(Blending)即B解法,分部-集成(Partition-Integration approach)即 P -I 解法、化整-集成(Globalization-Integration approach)即G -I解法以及集成-集成(Integration-Integration approach) 即I -I解法。

        2 局部性能約束優(yōu)化問題四種解法的演化

        2.1 P解法

        在Schmit[6]奠定結(jié)構(gòu)優(yōu)化成為獨立學(xué)科基礎(chǔ)之前,結(jié)構(gòu)設(shè)計的滿應(yīng)力準(zhǔn)則其實就是分部優(yōu)化解法。Shanley等[7,8]提出的同步失效準(zhǔn)則,為解決航空結(jié)構(gòu)受壓構(gòu)件設(shè)計需要確定每個截面一組尺寸集,本質(zhì)上還是分部優(yōu)化解法。文獻(xiàn)[1]除了敘述桁架滿應(yīng)力設(shè)計,還敘述了鋼框架和鋼筋混凝土框架的構(gòu)件優(yōu)化設(shè)計,也歸于分部優(yōu)化解法。分部優(yōu)化解法的前提是在超靜定結(jié)構(gòu)中引進(jìn)內(nèi)力暫時不變的處理,或者說引入了每輪設(shè)計中結(jié)構(gòu)靜定化假設(shè)。文獻(xiàn)[2-5]在連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化的一系列研究中,也沿用了結(jié)構(gòu)靜定化假設(shè)。

        2.2 G解法

        能否在建模時把各單元的相互關(guān)聯(lián)和影響考慮進(jìn)去?該想法導(dǎo)致產(chǎn)生了應(yīng)力全局化方法,即化整優(yōu)化方法[3-5]。

        2.3 I解法

        Yang等[9]最先提出了用K-S函數(shù)將眾多應(yīng)力約束集成為一個,進(jìn)行結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計。文獻(xiàn)[10,11]也采用凝聚函數(shù)包括K-S函數(shù)集成應(yīng)力約束的結(jié)構(gòu)拓?fù)鋬?yōu)化研究,并作了詳細(xì)闡述[2,3]。Paris等[12]及Luo等[13]先后引用K-S函數(shù)將結(jié)構(gòu)拓?fù)鋬?yōu)化應(yīng)力約束集成化。將眾多約束化為一個,是化整優(yōu)化方法和集成優(yōu)化方法的共性,但是前者借助于力學(xué)原理,后者則借助于數(shù)學(xué)公式。

        2.4 B解法

        B解法是新提出的交融優(yōu)化方法,細(xì)分為3種(P -I,G -I和I -I)解法,皆植根于互逆規(guī)劃的定理2[14,15]。3種交融優(yōu)化方法在每一次迭代求解中,都包括單目標(biāo)方和多目標(biāo)方兩部分,交替求解局部性約束的優(yōu)化模型,較之局部性約束的P,G和I解法,3種B解法分別提高了求解效率。

        3 疲勞性能的分部特點和化整變換

        3.1 疲勞性能約束同應(yīng)力約束的類比

        疲勞約束如同應(yīng)力約束,都屬于單元性能滿足型約束,相應(yīng)的結(jié)構(gòu)優(yōu)化模型表示為

        (1)

        若是靜定結(jié)構(gòu),其內(nèi)力不隨設(shè)計變量變化而改變,同應(yīng)力約束問題類似,可求解N個單元疲勞壽命子問題,求解后Li l(ti)不變,一次就求得結(jié)果。連續(xù)體屬于超靜定結(jié)構(gòu),不可能像靜定結(jié)構(gòu)那樣結(jié)構(gòu)分析一次就得到最優(yōu)解,但是可以借助靜定化假設(shè),迭代逼近式(1)的整體最優(yōu)解。仿滿應(yīng)力準(zhǔn)則,可以稱此為滿疲勞設(shè)計準(zhǔn)則。

        3.2 從單元疲勞壽命化整出結(jié)構(gòu)疲勞壽命

        既然應(yīng)力約束優(yōu)化可借助應(yīng)變能化整為全局約束,那么對于疲勞約束,也能做類似的工作。為此,考慮循環(huán)載荷下的應(yīng)力壽命S -N曲線,按巴士昆(Basquin)公式,有

        σ=AL- m

        (2)

        式中σ為疲勞循環(huán)下受拉實驗桿的材料極限應(yīng)力,L為對應(yīng)的疲勞壽命,A與m為材料給定的相關(guān)量,皆為正數(shù)。

        由式(2)可推得能量的表達(dá)式為

        (3)

        式中e為實驗桿的應(yīng)變能,v為桿的體積,E為楊氏模量。

        把式(3)推廣到單元復(fù)雜應(yīng)力狀態(tài)下的疲勞問題,可以得到

        (4)

        為了推導(dǎo)出疲勞約束全局化的公式,須把式(4)用于結(jié)構(gòu)疲勞問題的單元上,稍加整理,然后求和可以得到結(jié)構(gòu)疲勞壽命的表達(dá)式為

        (5)

        式中ei l與Li l分別為l工況下i單元的單元應(yīng)變能與單元疲勞壽命,vi為i單元體積。右邊表示l工況下i單元的壽命之和,故可以稱為疲勞載荷l工況下的結(jié)構(gòu)壽命。

        3.3 疲勞過濾函數(shù)及其在疲勞約束顯式化中的作用

        沿襲以往的研究,基于ICM方法,利用不同的過濾函數(shù),將結(jié)構(gòu)拓?fù)鋬?yōu)化模型近似顯式化。為此,類比對應(yīng)力約束的處理措施[3-5],給出疲勞約束對應(yīng)的有關(guān)公式。

        在應(yīng)力約束的拓?fù)鋬?yōu)化中,涉及到P解法、G解法和I解法,當(dāng)中相應(yīng)的過濾函數(shù)有不同的選擇。

        P解法的應(yīng)力約束若取應(yīng)力為約束[3,4],有

        (6)

        P解法的疲勞約束,類比應(yīng)力約束式(6)得

        (7)

        若取應(yīng)變能表示P解法的應(yīng)力約束,單元應(yīng)力約束和全局化的結(jié)構(gòu)應(yīng)力約束可表示為

        (8a)

        (8b)

        式(8a)左端為結(jié)構(gòu)應(yīng)變能,其單元應(yīng)變能為

        (9)

        式中σi l,εi l與Pi l為l載荷工況下i單元的應(yīng)力、應(yīng)變與節(jié)點力向量,Ki與Ωi分別為i單元的剛度矩陣及子域。

        (10)

        (11)

        (12)

        至此,將式(12)代入式(5)的第i項得

        (13)

        于是,得到G解法疲勞性能化整之后全局化的結(jié)構(gòu)壽命約束為

        (14)

        3.4 疲勞約束與應(yīng)力約束的統(tǒng)一規(guī)格

        為了疲勞約束優(yōu)化問題的求解與應(yīng)力約束優(yōu)化問題的求解一致,引入了倒壽命的概念,

        I=1/L

        (15)

        在倒壽命概念下,由式(7,13,14)得

        (16a)

        (16b)

        將式(12)代入式(8)得

        (17a)

        (17b)

        式(16,17)表明疲勞與應(yīng)力統(tǒng)一了約束規(guī)格。

        順便歸納一下過濾函數(shù)的取法,(1) 應(yīng)力約束時,P解法和I解法皆按式(6,17a)??;G解法按 式(17b)取。(2) 疲勞約束時,P解法和I解法按 式(16a)?。籊解法按式 (16b)取。

        4 局部性能約束單目標(biāo)模型三種解法

        本節(jié)補充和完整了已有的局部性能約束解法,屬于單目標(biāo)模型分別為分部、化整和集成解法三種,下面分別用P模型、G模型和I模型表示。

        4.1 P解法(Partition approach)

        (18)

        式中性能符號Ψ為應(yīng)力或疲勞壽命,當(dāng)為應(yīng)力時,按式(6)或式(17a)取;當(dāng)為疲勞壽命時,按式(16a)取。依靜定化假設(shè)求得

        (19)

        式(19)尚需按變量的上下限予以修正。

        4.2 G解法(Globalization approach)

        (20)

        其中性能符號Ψ分別由式(17b)或式(16b)取得。可采用DSQP或DP -EM方法[16]求解。

        4.3 I解法(Integration approach)

        (21)

        (22)

        (23)

        接著,將式(21)轉(zhuǎn)化到x空間里,通過目標(biāo)二階近似,用Lagrange乘子法求得顯式解,迭代逼近,得

        (24)

        回代到t空間里得式(21)的解。篇幅所限,詳細(xì)推導(dǎo)過程從略。

        5 局部性能約束的單目標(biāo)與多目標(biāo)模型交替迭代的三種交融解法

        5.1 三種交融解法

        本節(jié)突破和開拓出了局部性能約束的新解法,屬于單目標(biāo)與多目標(biāo)模型交替迭代的交融解法,共三種:(1) P -I 解法(Partition-Integration approach); (2) G -I 解法(Globalization-Integration approach); (3) I -I 解法(Integration-Integration approach)。

        標(biāo)識三種解法的前一個字母,代表第4節(jié)的三種單目標(biāo)模型;后一個字母,代表第5.3節(jié)即將介紹的多目標(biāo)的集成模型,唯有一種。

        5.2 交融解法的迭代尋優(yōu)算法

        (1) 迭代次數(shù)賦初值ν=1。

        (2) 實施模塊1,即進(jìn)行結(jié)構(gòu)分析。

        (3) 實施模塊2,建立4.1~4.3三種模型之一,求解得最優(yōu)重量W*。

        (6) 結(jié)束迭代。

        5.3 交融解法中的多目標(biāo)集成模型

        多目標(biāo)模型集成模型可寫為

        (25)

        式中目標(biāo)函數(shù)的含義與第4.1節(jié)P模型的約束相同,此處不再贅述。用Lagrange乘子法推導(dǎo)出規(guī)劃(25)的顯式解,

        (26)

        迭代求解中,當(dāng)t(v + 1)與t(v )充分接近后收斂,每次小循環(huán)迭代都用0≤ti≤1修改迭代解。限于篇幅,具體推導(dǎo)從略。

        6 數(shù)值算例

        6.1 MBB梁的應(yīng)力約束拓?fù)鋬?yōu)化

        基結(jié)構(gòu)是寬為300 mm、高為50 mm簡支梁,集中載荷F=200 N作用于跨中頂點,因?qū)ΨQ性取一半結(jié)構(gòu)和一半載荷進(jìn)行計算,計算簡圖和支座約束如圖1所示,采用150×50網(wǎng)格。彈性模量E= 2.1×105MPa,泊松比為0.3。為避免應(yīng)力集中,F(xiàn)分散作用在相鄰的4個結(jié)點上,右下角豎向約束也分散在2個結(jié)點上。許用應(yīng)力[σ]=100 MPa。

        圖1 MBB梁問題取一半結(jié)構(gòu)的分析及優(yōu)化模型

        本算例分別以單純集成方法(I解法)及集成-集成交融方法(I -I 解法)求解應(yīng)力約束體積極小化問題。I解法經(jīng)過61次迭代收斂,最優(yōu)體積比為 0.548,最大Mises應(yīng)力為100.373 MPa。I -I 解法經(jīng)過47次迭代收斂,最優(yōu)體積比為0.553,最大Mises應(yīng)力為99.900 MPa。兩種方法得到的最優(yōu)拓?fù)淙鐖D2所示,對應(yīng)的Mises應(yīng)力分布如圖3所示,目標(biāo)及約束的迭代歷史曲線對比如圖4所示。相比I解法,I -I解法的迭代次數(shù)減少了14次,收斂更快。

        圖2 單純集成方法與集成-集成交融方法得到的最優(yōu)拓?fù)?/p>

        圖3 單純集成方法與集成-集成交融方法的最優(yōu)時Mises應(yīng)力分布

        圖4 集成方法與集成-集成交融方法優(yōu)化過程目標(biāo)及約束對比

        6.2 兩端固支梁的應(yīng)力約束拓?fù)鋬?yōu)化

        如圖5所示,基結(jié)構(gòu)尺寸為120 mm×60 mm×1 mm的平面體,采用200×100網(wǎng)格。材料彈性模量為E=2.1×105MPa,泊松比ν=0.3。左右邊界固定。工況1,集中載荷F1=300 N作用于上邊界中點;工況2, 集中載荷F2=300 N作用于下邊界中點。劃分為200×100個矩形單元。應(yīng)力約束為100 MPa。

        圖5 兩工況算例分析及優(yōu)化的模型

        I解法經(jīng)過64次迭代收斂,最優(yōu)體積比為 0.170,最大Mises應(yīng)力為102.590 MPa。I -I 解法經(jīng)過57次迭代收斂,最優(yōu)體積比為0.142,最大Mises應(yīng)力為96.58 MPa。兩種方法得到的最優(yōu)拓?fù)淙鐖D6所示,對應(yīng)的兩個工況下的Mises應(yīng)力分布如圖7和圖8所示,目標(biāo)及約束的迭代歷史曲線對比如圖9所示。經(jīng)對比可以看出,兩種方法得到的最優(yōu)拓?fù)錁?gòu)型不僅相似,對應(yīng)的Mises應(yīng)力也是相似的。I -I 解法得到的最優(yōu)體積比略小。相比I解法,I -I 解法的迭代次數(shù)減少了7次,收斂 更快。

        圖6 集成方法與集成-集成交融方法得到的最優(yōu)拓?fù)?/p>

        圖7 集成方法與集成-集成交融方法的最優(yōu)時工況1的Mises應(yīng)力分布

        圖8 集成方法與集成-集成交融方法的最優(yōu)時工況2的Mises應(yīng)力分布

        圖9 集成方法與集成-集成交融方法優(yōu)化過程目標(biāo)及約束對比

        6.3 五次超靜定梁的疲勞約束拓?fù)鋬?yōu)化

        如圖10所示,材料性能為E=2.1×105MPa,泊松比ν=0.3,材料對應(yīng)于循環(huán)次數(shù)1.0×106的疲勞極限為σs=250 MPa,基結(jié)構(gòu)為 800 mm×200 mm×10 mm的平面體,劃分為400×100個矩形單元。基結(jié)構(gòu)的四個角點皆為雙向約束鉸支座。工況為正弦形式循環(huán)載荷F=21000 N且均值為0,相位角為0,作用于圖10基結(jié)構(gòu)的中心點,為了避免應(yīng)力集中將荷載施加在五個節(jié)點上。疲勞壽命約束為大于或等于1.0×105次。Basquin公式中常數(shù)值取m=0.10,A=995.3 MPa。

        圖10 五次超靜定梁結(jié)構(gòu)的優(yōu)化模型

        Fig.10 Model of structure for analysis and optimization of the statically indeterminate beam

        本算例分別以單純分部方法(P解法)及分部-集成交融方法(P -I 解法)求解疲勞壽命約束體積極小化問題。P解法經(jīng)過36次迭代收斂,最優(yōu)體積比為0.488,最大疲勞壽命為99141次。P -I 解法經(jīng)過26次迭代收斂,最優(yōu)體積比為0.483,最大疲勞壽命為100675次。兩種方法得到的最優(yōu)拓?fù)淙?圖11 所示,對應(yīng)的疲勞壽命以10為底的對數(shù)分布如圖12所示,目標(biāo)及約束的迭代歷史曲線對比如 圖13 所示。相比P解法,P -I 解法的迭代次數(shù)減少了10次,收斂更快。

        圖11 單純分部方法與分部-集成交融方法得到的最優(yōu)拓?fù)?/p>

        圖12 分部方法與分部-集成交融方法的最優(yōu)時疲勞壽命分布

        圖13 分部方法與分部-集成交融方法優(yōu)化過程目標(biāo)及約束對比

        7 結(jié) 論

        類比對應(yīng)力約束的處理措施,給出疲勞約束對應(yīng)的過濾函數(shù)。仿應(yīng)力約束問題的滿應(yīng)力準(zhǔn)則,介紹了滿疲勞設(shè)計準(zhǔn)則。仿應(yīng)力約束問題的化整方法的全局化應(yīng)力約束,給出了疲勞約束化整方法的全局化疲勞約束?;诘箟勖母拍睿瑢崿F(xiàn)了疲勞約束與應(yīng)力約束的規(guī)格統(tǒng)一。對結(jié)構(gòu)疲勞約束與應(yīng)力約束的局部性能約束問題,梳理了3種單目標(biāo)優(yōu)化方法,闡述了新提出的由單目標(biāo)與多目標(biāo)交替迭代的3種交融優(yōu)化方法。數(shù)值算例表明,交融解法的求解效率,比單目標(biāo)解法更加有效,這是每次結(jié)構(gòu)分析之后做了兩次尋優(yōu)的必然結(jié)果。

        致謝:碩士生張譯文同學(xué)參加了算例計算。

        參考文獻(xiàn)(References):

        [1] 錢令希.工程結(jié)構(gòu)優(yōu)化設(shè)計[M].北京:水利電力出版社,1983.(QIAN Ling-xi.OptimizationDesignforEngineeringStructures[M].Beijing:Water Conservancy and Electric Power Press,1983.(in Chinese))

        [2] 隋允康.建?!ぷ儞Q·優(yōu)化:結(jié)構(gòu)綜合方法新進(jìn)展[M].大連:大連理工大學(xué)出版社,1996.(SUI Yun-kang.ModelingTransformaionandOptimzation:NewDevelopmentofStructuralSsynthesisMethod[M].Dalian:Dalian University of Technology Press,1996.(in Chinese))

        [3] 隋允康,葉紅玲.連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化的ICM方法[M].北京:科學(xué)出版社,2013.(SUI Yun-kang,YE Hong-ling.ContinuumTopologyOptimizationICMMethod[M].Beijing:Science Press,2013.(in Chinese))

        [4] Sui Y K,Peng X R.Modeling,SolvingandApplicationforTopologyOptimizationofContinuumStructures,ICMMethodBasedonStepFunction[M].Elsevier,2018.

        [5] 隋允康,葉紅玲,彭細(xì)榮.應(yīng)力約束全局化策略下的連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化[J].力學(xué)學(xué)報,2006,38(3):364-370.(SUI Yun-kang,YE Hong-ling,PENG Xi-rong.Topologicaal optimization of continuum structure under the strategy of globalization of stress constraints [J].ChineseJournalofTheoreticalandAppliedMechanics,2006,38(3):364-370.(in Chinese))

        [6] Schmit L A.Structural design by systematic synthesis[A].Proceeding of Conference on Electronic Computation,ASCE[C].New York, 1960.

        [7] Shanley F R.Weight-StrengthAnalysisofAircraftStructures[M].New York,1960.

        [8] Spunt L.OptimumStructuralDesign[M].New Jersey,1971.

        [9] Yang R J,Chen C J.Stress -based topology optimization[J].StructuralOptimization,1996,12(2):98-105.

        [10] 隋允康,張學(xué)勝,龍連春.應(yīng)力約束處理為應(yīng)變能集成的連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化[J].計算力學(xué)學(xué)報,2007,24(5):602-608.(SUI Yun-kang,ZHANG Xue -sheng,LONG Lian-chun.ICM method of the topology optimization for continuum structures with stress constraints approached by the integration of strain energies[J].ChineseJournalofComputationalMecha-nics,2007,24(5):602-608.(in Chinese))

        [11] 隋允康,鐵 軍.結(jié)構(gòu)拓?fù)鋬?yōu)化ICM顯式化與拋物型凝聚函數(shù)對于應(yīng)力約束的集成化[J].工程力學(xué),2010,27(2):124-133.(SUI Yun-kang,TIE Jun.The ICM explicitation approach to the structural topology optimization and the integrating approach to stress constraints based on the parabolic aggregation function [J].EngineeringMechanics,2010,27(s2):124-133.(in Chinese))

        [12] París J,Navarrina F,Colominas I,et al.Improvements in the treatment of stress constraints in structural topology optimization problems [J].JournalofComputationalandAppliedMathematics,2010,234(7):2231-2238.

        [13] Luo Y J,Wang M,Kang Z.An enhanced aggregation method for topology optimization with local stress constraints[J].ComputerMethodsinAppliedMechanicsandEngineering,2013,254(2):31-41.

        [14] 隋允康,彭細(xì)榮,葉紅玲,等.互逆規(guī)劃理論及其用于建立結(jié)構(gòu)拓?fù)鋬?yōu)化的合理模型[J].力學(xué)學(xué)報,2019,51(6):1-9.(SUI Yun-kang,PENG Xi-rong,YE Hong-ling,et al.Reciprocal programming theory and its application to establish a reasonable model of structural topology optimization[J].ChineseJournalofTheoreticalandAppliedMechanics,2019,51(6):1940-1948.(in Chinese))

        [15] 鐵 軍,隋允康,彭細(xì)榮.互逆規(guī)劃的拓寬和深化及其在結(jié)構(gòu)拓?fù)鋬?yōu)化的應(yīng)用.力學(xué)學(xué)報,2020,52(6):1822-1837.(TIE Jun,SUI Yun-kang,PENG Xi-rong.Widening and deepening of reciprocal programming and its application to structural topology optimization[J].ChineseJournalofTheoreticalandAppliedMe-chanics,2020,52(6):1822-1837.(in Chinese))

        [16] Sui Y K,Peng X R.Explicit model of dual programming and solving method for a class of separable convex programming problems[J].EngineeringOptimization,2019,51(9):1604-1625.

        猜你喜歡
        優(yōu)化結(jié)構(gòu)方法
        超限高層建筑結(jié)構(gòu)設(shè)計與優(yōu)化思考
        《形而上學(xué)》△卷的結(jié)構(gòu)和位置
        民用建筑防煙排煙設(shè)計優(yōu)化探討
        關(guān)于優(yōu)化消防安全告知承諾的一些思考
        一道優(yōu)化題的幾何解法
        論結(jié)構(gòu)
        中華詩詞(2019年7期)2019-11-25 01:43:04
        可能是方法不對
        論《日出》的結(jié)構(gòu)
        用對方法才能瘦
        Coco薇(2016年2期)2016-03-22 02:42:52
        四大方法 教你不再“坐以待病”!
        Coco薇(2015年1期)2015-08-13 02:47:34
        亚洲伊人久久一次| 亚洲av中文无码乱人伦在线观看| 成人无码一区二区三区| 极品少妇一区二区三区四区| 俺来也俺去啦久久综合网| 国产黑色丝袜在线观看视频| 一区二区三区夜夜久久| 国产精品无套一区二区久久 | AV无码人妻一区二区三区牛牛| 久久精品亚洲一区二区三区画质| 男女性杂交内射女bbwxz| 一本加勒比hezyo无码人妻 | 粗大挺进尤物人妻一区二区| 久久精品国产福利亚洲av| 免费观看91色国产熟女| 国产精品综合一区二区三区| 国产成人亚洲综合一区| 亚洲视频综合在线第一页| 中文字幕精品人妻在线| 精品丰满人妻无套内射| 97se亚洲国产综合自在线图片| 蜜桃视频色版在线观看| 国产剧情av麻豆香蕉精品| 又色又爽又高潮免费视频国产| 亚洲熟妇av乱码在线观看 | 亚洲av成人综合网| 美女熟妇67194免费入口| 国产毛片视频一区二区三区在线| 丰满熟妇人妻av无码区| 色噜噜狠狠色综合成人网 | 免费的黄网站精品久久| 日韩不卡的av二三四区| 亚洲avav天堂av在线网爱情| av中文字幕不卡无码| 麻豆夫妻在线视频观看| 香蕉久久一区二区不卡无毒影院| 亚洲成av人片在线观看ww| 日韩久久久久中文字幕人妻| 大陆少妇一区二区三区| 国产精品久久久久9999赢消| 国产精品天堂avav在线|