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

        ?

        考慮熱振特性的連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)

        2016-10-24 05:18:16劉遠(yuǎn)東尹益輝
        振動(dòng)與沖擊 2016年17期
        關(guān)鍵詞:振動(dòng)優(yōu)化結(jié)構(gòu)

        劉遠(yuǎn)東,莫 軍,尹益輝,徐 兵

        (中國(guó)工程物理研究院 總體工程研究所,四川 綿陽(yáng) 621900)

        ?

        考慮熱振特性的連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)

        劉遠(yuǎn)東,莫軍,尹益輝,徐兵

        (中國(guó)工程物理研究院 總體工程研究所,四川 綿陽(yáng)621900)

        為實(shí)現(xiàn)傳熱和振動(dòng)條件下連續(xù)體結(jié)構(gòu)的拓?fù)鋬?yōu)化設(shè)計(jì),以結(jié)構(gòu)散熱弱度最小化和動(dòng)態(tài)特征值最大化加權(quán)函數(shù)為目標(biāo),建立傳熱和振動(dòng)條件下連續(xù)體結(jié)構(gòu)的多目標(biāo)拓?fù)鋬?yōu)化模型,實(shí)現(xiàn)了相應(yīng)的算法和算例。方法中采用Rational Approximation of Material Properties(RAMP)方法對(duì)密度進(jìn)行懲罰,利用優(yōu)化準(zhǔn)則法控制設(shè)計(jì)目標(biāo)與材料分布,以敏度過(guò)濾技術(shù)抑制棋盤(pán)格效應(yīng),通過(guò)歸一化目標(biāo)函數(shù)有效地避免了不同性質(zhì)目標(biāo)函數(shù)的量級(jí)差異。通過(guò)算例,獲得了熱-振權(quán)重系數(shù)對(duì)結(jié)構(gòu)拓?fù)錁?gòu)型和目標(biāo)函數(shù)(宏觀(guān)結(jié)構(gòu)的散熱弱度和基頻)的影響規(guī)律,算例結(jié)果表明了該方法的有效性。

        連續(xù)體結(jié)構(gòu);熱-振多目標(biāo);拓?fù)鋬?yōu)化;優(yōu)化準(zhǔn)則法

        拓?fù)鋬?yōu)化研究大多集中在單一物理場(chǎng)條件下的設(shè)計(jì)問(wèn)題,如力場(chǎng)下的結(jié)構(gòu)設(shè)計(jì)、溫度場(chǎng)下的散熱設(shè)計(jì)等。這些單一場(chǎng)下的結(jié)構(gòu)優(yōu)化研究對(duì)工程實(shí)際結(jié)構(gòu)設(shè)計(jì)具有一定的指導(dǎo)意義[1-6]。然而,隨著研究的深入和優(yōu)化設(shè)計(jì)應(yīng)用工程的要求,有必要開(kāi)展更符合結(jié)構(gòu)實(shí)際工作情況的多場(chǎng)條件下的優(yōu)化設(shè)計(jì)研究。實(shí)際工程中,工作條件一般為力場(chǎng)、溫度場(chǎng)、電磁場(chǎng)、流場(chǎng)等多場(chǎng)耦合。其中高溫-振動(dòng)復(fù)合問(wèn)題就是一種典型的多場(chǎng)耦合問(wèn)題,例如航空發(fā)動(dòng)機(jī)、核電站、高速飛行器等。這類(lèi)結(jié)構(gòu)通常工作在高溫-振動(dòng)環(huán)境下,因此熱應(yīng)力和振動(dòng)成為影響結(jié)構(gòu)剛度、應(yīng)力、頻率和壽命的重要因素,有必要開(kāi)展熱-振作用下結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)研究。熱-振作用下結(jié)構(gòu)的拓?fù)鋬?yōu)化設(shè)計(jì)一般為:一方面需要熱控制,具有足夠的散熱性或隔熱性,以便使結(jié)構(gòu)獲得合適的溫度;另一方面又需要振動(dòng)控制,使結(jié)構(gòu)遠(yuǎn)離工作頻率,以避免共振現(xiàn)象。

        目前,對(duì)于熱激勵(lì)結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì),由于熱載荷屬于隱含力,依賴(lài)于材料的分布,導(dǎo)致熱彈性耦合結(jié)構(gòu)的拓?fù)鋬?yōu)化求解存在很大困難,相關(guān)的研究工作很少。LI等[7]基于ESO方法,以有限元離散厚度為設(shè)計(jì)變量、材料用量為約束,開(kāi)展了均勻穩(wěn)態(tài)溫度場(chǎng)下熱、力耦合拓?fù)鋬?yōu)化問(wèn)題的研究。CHO等[8]則基于靈敏度分析方法研究了穩(wěn)態(tài)溫度場(chǎng)下熱彈性弱耦合問(wèn)題的結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)。然而,上述文獻(xiàn)中熱力耦合拓?fù)鋬?yōu)化結(jié)果存在嚴(yán)重的中間密度分布。左孔天等[9]以解耦的形式,建立了散熱弱度最小和柔順度最小的優(yōu)化模型,開(kāi)展了熱力多目標(biāo)拓?fù)鋬?yōu)化設(shè)計(jì)研究。在以上工作基礎(chǔ)上,張衛(wèi)紅等[10]討論了基于最小柔順度和彈性應(yīng)變能為目標(biāo)函數(shù)的兩種熱彈性拓?fù)鋬?yōu)化設(shè)計(jì)方法,并給出了熱、彈性載荷對(duì)優(yōu)化結(jié)果的影響。ZHANG等[11-12]關(guān)于熱振拓?fù)鋬?yōu)化的工作是建立在均勻變化溫度場(chǎng)下的,即在力場(chǎng)范圍內(nèi)通過(guò)線(xiàn)脹系數(shù)和溫度差來(lái)表征熱載荷,在某種程度上可以退化為振動(dòng)拓?fù)鋬?yōu)化,因此有必要開(kāi)展考慮傳熱目標(biāo)函數(shù)的更一般形式熱-振拓?fù)鋬?yōu)化研究。

        本文針對(duì)輕質(zhì)/隔熱/減振結(jié)構(gòu)的多目標(biāo)拓?fù)鋬?yōu)化問(wèn)題,提出以散熱弱度和動(dòng)態(tài)特征值加權(quán)為綜合目標(biāo)的優(yōu)化模型。針對(duì)不同性質(zhì)的結(jié)構(gòu)響應(yīng)量在數(shù)量級(jí)上的差異,歸一化處理各單一目標(biāo)函數(shù);采用優(yōu)化準(zhǔn)則法和敏度過(guò)濾技術(shù)控制設(shè)計(jì)目標(biāo)和材料分布。比較不同熱-振權(quán)系數(shù)下的多目標(biāo)優(yōu)化結(jié)果,揭示了不同權(quán)系數(shù)對(duì)目標(biāo)函數(shù)的影響規(guī)律。結(jié)果表明同時(shí)考慮傳熱特性和振動(dòng)特性的連續(xù)體結(jié)構(gòu)多目標(biāo)拓?fù)鋬?yōu)化更有利于提高結(jié)構(gòu)的綜合性能。由于目標(biāo)函數(shù)對(duì)密度靈敏度是建立在單元熱能、應(yīng)變能和動(dòng)能的基礎(chǔ)上,可以有效的移植到商業(yè)軟件的二次程序開(kāi)發(fā)中,有利于開(kāi)展工程計(jì)算。同時(shí)算例結(jié)果表明該方法對(duì)于凸問(wèn)題是有效性。

        1 熱-振多目標(biāo)優(yōu)化設(shè)計(jì)模型

        1.1歸一化多目標(biāo)函數(shù)

        對(duì)結(jié)構(gòu)設(shè)計(jì)區(qū)域進(jìn)行有限單元離散后,采用常見(jiàn)的RAMP方法,有限單元材料的熱傳導(dǎo)系數(shù)H、彈性模量E和質(zhì)量M的插值計(jì)算式為:

        (1)

        (2)

        (3)

        式中:H0、E0和M0分別為實(shí)體材料的熱傳導(dǎo)系數(shù)、彈性模量和質(zhì)量,xi∈[0,1]為單元材料的偽密度值,α為懲罰因子,對(duì)介于0和1之間的偽密度值進(jìn)行懲罰,見(jiàn)文獻(xiàn)[13],本文后面取α=3。對(duì)結(jié)構(gòu)設(shè)計(jì)區(qū)域進(jìn)行有限單元離散后,采用敏度過(guò)濾技術(shù)[13]抑制有限元計(jì)算中的棋盤(pán)格現(xiàn)象。

        Findxi(i=1,2,…,n)

        (k=1,2,…,ndy)

        (4)

        ConstraintⅠ:P=HT

        ConstraintⅡ:(K-λkM)φk=0

        ConstraintⅣ:0

        i=1,2,…,n

        式中,βr是權(quán)重系數(shù),φ表示結(jié)構(gòu)的散熱弱度,λk是結(jié)構(gòu)第k階圓頻率的平方(特征方程的特征值)。約束Ⅰ是穩(wěn)態(tài)熱傳導(dǎo)有限元平衡方程,P為熱載荷列陣,H為熱傳導(dǎo)矩陣,T為節(jié)點(diǎn)溫度列陣;約束Ⅱ是自由振動(dòng)控制方程,φk是結(jié)構(gòu)的第k階特征向量,與λk形成特征對(duì),M為結(jié)構(gòu)的質(zhì)量矩陣;約束Ⅲ給出了實(shí)體材料用量上限,vi為單元i的體積,V為給定的實(shí)體材料用量。約束Ⅳ是密度變量的上下限約束,xmin為趨于0的正數(shù)(取10-3),以避免計(jì)算中出現(xiàn)剛度矩陣奇異,xmax是有限單元的最大密度值。式(4)同時(shí)考慮了結(jié)構(gòu)的傳熱特性和振動(dòng)特性,該拓?fù)鋬?yōu)化設(shè)計(jì)方法有利于連續(xù)體結(jié)構(gòu)有效性能的提高。

        1.2靈敏度分析

        熱振作用下連續(xù)體拓?fù)鋬?yōu)化設(shè)計(jì)的目標(biāo)函數(shù)是一個(gè)綜合考慮傳熱和振動(dòng)特性聯(lián)合優(yōu)化的問(wèn)題。在目標(biāo)函數(shù)歸一化后,對(duì)傳熱問(wèn)題和振動(dòng)問(wèn)題的靈敏度求解分別如下:

        對(duì)于熱傳導(dǎo)問(wèn)題,結(jié)構(gòu)總熱能表示為

        (5)

        式中,Tene是熱能,ti和hi分別表示與第i個(gè)單元對(duì)應(yīng)的溫度和單元熱傳導(dǎo)矩陣。

        由式(5)求得散熱弱度φ對(duì)偽密度xi的靈敏度為

        (6)

        式中:φ為結(jié)構(gòu)的散熱弱度。

        對(duì)于動(dòng)力問(wèn)題,結(jié)構(gòu)系統(tǒng)自然圓頻率與相應(yīng)模態(tài)向量的關(guān)系為

        (7)

        將質(zhì)量矩陣正則化,可得圓頻率對(duì)偽密度xi的靈敏度為

        (8)

        式中:Tenei、Dsenei和Dkenei分別表示單元i的熱能、應(yīng)變能和動(dòng)能。由于目標(biāo)函數(shù)(散熱弱度+頻率)對(duì)密度的靈敏度均建立在單元熱能、應(yīng)變能和動(dòng)能的基礎(chǔ)上,可以方便的移植到商業(yè)軟件的二次程序開(kāi)發(fā)中,因此式(4)的優(yōu)化模型計(jì)算效率高,更適于求解工程優(yōu)化問(wèn)題。

        1.3設(shè)計(jì)變量的迭代格式

        基于變分原理的泛函形式,引入Lagrange乘子構(gòu)造的Lagrange函數(shù),將有約束的問(wèn)題轉(zhuǎn)化為無(wú)約束的極值問(wèn)題,求解目標(biāo)函數(shù)滿(mǎn)足一系列約束條件的駐值問(wèn)題,基于Karush-Kuhn-Tucker(KKT)條件,結(jié)合文獻(xiàn)[15],可得如下變量更新迭代格式:

        (9)

        式中:xi為設(shè)計(jì)變量,ξ為阻尼系數(shù),文中取0.5,m為移動(dòng)極限,可以在0與1之間變化,本文取0.3,λj為約束的拉氏算子,且有

        (10)

        2 優(yōu)化步驟和方法

        對(duì)熱-振作用下連續(xù)體結(jié)構(gòu)進(jìn)行拓?fù)鋬?yōu)化設(shè)計(jì)的流程圖見(jiàn)圖1,基本步驟如下:

        1)初始化,宏觀(guān)結(jié)構(gòu)的傳熱學(xué)、動(dòng)力學(xué)分析:劃分有限單元網(wǎng)格、計(jì)算散熱弱度和基頻;

        2)按式(6)和式(8)進(jìn)行連續(xù)體結(jié)構(gòu)熱-振多目標(biāo)優(yōu)化的靈敏度分析;

        3)按式(9)進(jìn)行設(shè)計(jì)變量迭代,采用敏度過(guò)濾技術(shù)抑制棋盤(pán)格現(xiàn)象;

        4)設(shè)計(jì)變量更新;

        5)判斷是否滿(mǎn)足約束條件,如果不滿(mǎn)足,則重復(fù)步驟2)~4),如果滿(mǎn)足,則進(jìn)行后處理,輸出相關(guān)參數(shù)。

        圖1 熱-振作用下連續(xù)體結(jié)構(gòu)多目標(biāo)優(yōu)化設(shè)計(jì)流程圖Fig.1 The flowchart for thermal-vibrational multi-objective optimal design of continuum structure

        3 數(shù)值算例

        設(shè)計(jì)區(qū)域?yàn)橐徽叫伟澹鐖D2所示,其邊長(zhǎng)a=40,厚度t=5?;w材料的彈性模量E=2×105,泊松比ν=0.3,密度ρ=1×10-6,比熱容C0=875,熱導(dǎo)率K0=10。這里各物理量均為無(wú)量綱量。材料用量為初始結(jié)構(gòu)的50%。有限元網(wǎng)格為40×40。板的中心點(diǎn)受熱流率Q0=10的作用,板的右下角頂點(diǎn)受集中質(zhì)量m=2的作用。傳熱學(xué)的邊界條件為中心加熱,四邊溫度為T(mén)=0℃;力學(xué)邊界條件為左邊固支。表1給出了5種不同熱、振權(quán)系數(shù)條件下結(jié)構(gòu)的優(yōu)化設(shè)計(jì)結(jié)果。β1和β2分別表示熱、振權(quán)重系數(shù)。

        圖2 設(shè)計(jì)區(qū)域Fig.2 Design domain

        權(quán)系數(shù)β1β2散熱弱度?基頻f拓?fù)湫问焦r11042.880.12工況20.80.251.358.06工況30.50.552.28.36工況40.20.857.569.14工況5015263.39.31

        由表1可得,熱傳導(dǎo)單目標(biāo)拓?fù)鋬?yōu)化后,板的散熱弱度φ=51.35,板的基頻f=0.12;動(dòng)態(tài)單目標(biāo)拓?fù)鋬?yōu)化后,散熱弱度φ=5 263.3,板的基頻f=9.31。熱-振權(quán)系數(shù)β1=0.5、β2=0.5時(shí),拓?fù)鋬?yōu)化后,板的散熱弱度φ=52.2,基頻f=8.36??梢钥闯?,拓?fù)鋬?yōu)化后,結(jié)構(gòu)的熱、振目標(biāo)函數(shù)分別隨著熱、振權(quán)系數(shù)的增加而變得更好,隨著熱、振權(quán)系數(shù)的減少而變得較差。同時(shí)可以看出單目標(biāo)優(yōu)化的單項(xiàng)結(jié)果優(yōu)于多目標(biāo)優(yōu)化對(duì)應(yīng)的結(jié)果,但結(jié)構(gòu)優(yōu)化是在各種可能的結(jié)構(gòu)形式中取得最佳“折衷”的設(shè)計(jì)過(guò)程,因此熱-振多目標(biāo)拓?fù)鋬?yōu)化同時(shí)考慮傳熱和振動(dòng)特性,更能體現(xiàn)結(jié)構(gòu)綜合性能的提高。

        4 結(jié) 論

        本文建立了熱-振多目標(biāo)優(yōu)化設(shè)計(jì)模型,實(shí)現(xiàn)了相應(yīng)的算法和算例。通過(guò)算例,驗(yàn)證了方法的有效性,得到了熱-振權(quán)重系數(shù)對(duì)目標(biāo)函數(shù)和結(jié)構(gòu)構(gòu)型的影響規(guī)律。根據(jù)結(jié)果,單目標(biāo)優(yōu)化的單項(xiàng)結(jié)果優(yōu)于多目標(biāo)優(yōu)化對(duì)應(yīng)的結(jié)果,但多目標(biāo)優(yōu)化更好地滿(mǎn)足連續(xù)體結(jié)構(gòu)在使用中需要的綜合性能。同時(shí)將目標(biāo)函數(shù)的靈敏度均建立在單元熱能、應(yīng)變能和動(dòng)能的基礎(chǔ)上,使計(jì)算效率更高,可以方便地移植到商業(yè)軟件的二次程序開(kāi)發(fā)中,更適合于工程優(yōu)化問(wèn)題。

        [1]周克民,李俊峰,李霞.結(jié)構(gòu)拓?fù)鋬?yōu)化研究方法綜述[J].力學(xué)進(jìn)展,2005,35(1):69-76.

        ZHOU Kemin,LI Junfeng,LI Xia.A review on topology optimization of structures [J].Advances in Mechanics,2005,35(1):69-76.

        [2]羅震,陳立平,黃玉盈,等.連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)[J].力學(xué)進(jìn)展,2004,34(4):463-476.

        LUO Zhen,CHEN Liping,HUANG Yuying.Topological optimizationdesign for continuum structures [J].Advances in Mechanics,2004,34(4):463-476.

        [3]MATHIEU-POTVIN F,GOSSELIN L.Optimal conduction pathways for a heat generating body:a comparison exercise [J].International Journal of Heat Mass Transfer,2007,50(5-16):2996-3006.

        [4]LI Q,STEVEN G P,XIE Y M,et al.Evolutionary topology optimization for temperature reduction of heat conducting fields [J].International Journal of Heat and Mass Transfer,2004,47(23):5071-5083.

        [5]GERSBORG-HANSEN A,BENDSOE M P,SIGMUND O.Topology optimization of heat conduction problems using the finite volume method [J].Structural and Multidisciplinary Optimization,2006,31(4):251-259.

        [6]ZHANG Y,LIU S.Design of conducting paths based on topology optimization [J].Heat and Mass Transfer,2008,44(10):1217-1227.

        [7]LI Q,STEVEN G P,XIE Y M.Displacement minimization of thermo-elastic structures by evolutionary thickness design [J].Computer Methods in Applied Mechanics Engineering,1999,179(3/4):361-378.

        [8]CHO S,CHOI J Y.Efficient topology optimization of thermo elasticity problems using coupled field adjoint sensitivity analysis method [J].Finite Elements in Analysis and Design,2005,41(51):1481-1495.

        [9]左孔天,趙雨?yáng)|,陳立平.傳熱結(jié)構(gòu)的多目標(biāo)拓?fù)鋬?yōu)化設(shè)計(jì)研究[J].計(jì)算力學(xué)學(xué)報(bào),2007,24(5):620-627.

        ZUO Kongtian,ZHAO Yudong,CHEN Liping.Study on multiple objective topology optimization of thermal conductive structure [J].Chines Journal of Computational Mechanics,2007,24(5):620-627.

        [10]ZHANG W H,YANG J G,XU X J,et al.Topology optimization of thermoelastic structures:mean compliance minimization or elastic strain energy minimization [J].Structural and Multidisciplinary Optimization,2014,49(3):417-429.

        [11]YANG X W,LI Y M.Topology optimization to minimize the dynamic compliance of a bi-material plate in a thermal environment [J].Structural and Multidisciplinary Optimization,2013,47(3):399-408.

        [12]YANG X W,LI Y M.Structural topology optimization on dynamic compliance at resonance frequency in thermal environments [J].Structural and Multidisciplinary Optimization,2014,49(1):81-91.

        [13]SIGMUND O.Morphology-based black and white filters for topology optimization[J].Struct Multidisc Optim,2007,33 (4/5):401-424.

        [14]劉遠(yuǎn)東.材料/結(jié)構(gòu)一體化靜動(dòng)力學(xué)多目標(biāo)優(yōu)化設(shè)計(jì)和兩尺度實(shí)驗(yàn)研究[D].綿陽(yáng):中國(guó)工程物理研究院總體工程研究所,2007.

        [15]BENDSOE M P,SIGMUND O.Topology optimization theory,methods and application[J].Berlin:Springer,2003:71-75.

        Topologic optimization design for a continuum structure considerating its thermal and vibrational performances

        LIU Yuandong,MO Jun,YIN Yihui,XU Bing

        (Institute of Structural Mechanics,CAEP,P.O.Box 919-401,Mianyang 621900,China)

        For thermal and vibrational design problems,taking the minimum dissipation of heat transport potential capacity and the maximum fundamental frequency as objectives,a thermal and vibrational collaborative optimization model for a continuum structure was developed and corresponding numerical simulations were conducted.In the model,the rational approximation of material properties (RAMP)method was adopted to ensure clear topologies,the design objectives and material distribution were controlled with the optimization criteria method and the checkerboard effect was eliminated with the sensitivity filtering technique.To improve the smoothness of objective functions and avoid the singularity of numerical computation,the weighted objective functions were normalized.Based on numerical simulations,the effects of weight coefficients of thermal and vibrational optimization on the topologic configuration and objectives (the dissipation of heat transport potential capacity and fundamental frequency)were investigated and the results indicated that the proposed method is effective.

        continuum structure; thermal-vibrational multiple-objective; topologic optimization;optimization criteria method

        國(guó)家自然基金重點(diǎn)項(xiàng)目(11432011);中國(guó)工程物理研究院重點(diǎn)學(xué)科項(xiàng)目“計(jì)算固體力學(xué)”資助

        2015-04-07修改稿收到日期:2015-07-27

        劉遠(yuǎn)東 男,副研究員,博士,1978年生

        尹益輝 男,研究員,博士,1965年生

        TB303

        A DOI:10.13465/j.cnki.jvs.2016.17.019

        猜你喜歡
        振動(dòng)優(yōu)化結(jié)構(gòu)
        振動(dòng)的思考
        超限高層建筑結(jié)構(gòu)設(shè)計(jì)與優(yōu)化思考
        《形而上學(xué)》△卷的結(jié)構(gòu)和位置
        民用建筑防煙排煙設(shè)計(jì)優(yōu)化探討
        關(guān)于優(yōu)化消防安全告知承諾的一些思考
        一道優(yōu)化題的幾何解法
        振動(dòng)與頻率
        論結(jié)構(gòu)
        中立型Emden-Fowler微分方程的振動(dòng)性
        論《日出》的結(jié)構(gòu)
        国产人妖av在线观看| 欧美自拍区| 亚洲国产成人无码电影| 自拍偷拍韩国三级视频| 在线观看老湿视频福利| 国产suv精品一区二区| 91精品国产91久久久无码色戒| 视频在线亚洲视频在线| 97在线视频人妻无码| 亚洲精品无播放器在线播放| 少妇的诱惑免费在线观看| 五月婷婷丁香视频在线观看| 国产a在亚洲线播放| 国精产品一品二品国在线| 在线偷窥制服另类| 亚洲一区二区三区日韩在线观看| 日本大骚b视频在线| 国产A√无码专区| 亚洲成人激情在线影院| 中文字幕国产精品一二三四五区| 亚洲性爱视频| 本道无码一区二区久久激情| 麻豆国产精品伦理视频| 中文字幕日韩欧美一区二区三区| 久久久久国产精品免费免费搜索 | 亚洲粉嫩av一区二区黑人| 国产亚洲精品90在线视频| 国产女人水真多18毛片18精品| 亚洲区在线播放| 精品人妻一区二区视频| 无码色av一二区在线播放| 久久久久亚洲av无码专区| 中文字幕亚洲精品码专区| 男女主共患难日久生情的古言 | 国产一区国产二区亚洲精品| 欧美精品中文字幕亚洲专区| 亚洲男人天堂av在线| 精品嫩模福利一区二区蜜臀| 久久香蕉国产线看观看精品yw| 亚洲视频1区| 亚洲精品国产成人久久av盗摄|