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

        ?

        基于敏度分析的框架抗震優(yōu)化

        2012-02-13 09:01:44劉齊茂燕柳斌
        振動與沖擊 2012年17期
        關(guān)鍵詞:一階二階層間

        劉齊茂,燕柳斌

        (1.廣西工學院 土木建筑系,廣西 柳州 545006;2.廣西大學 土木建筑學院,南寧 530004)

        結(jié)構(gòu)的抗震優(yōu)化設(shè)計屬于結(jié)構(gòu)動力響應優(yōu)化設(shè)計的范疇,具有計算量巨大、優(yōu)化設(shè)計問題建模非常復雜等特點,因此該領(lǐng)域的研究進展非常緩慢。根據(jù)各國規(guī)范的規(guī)定、震害經(jīng)驗和實驗研究結(jié)果及工程實例分析,在地震作用下采用多高層建筑結(jié)構(gòu)層間相對位移作為衡量結(jié)構(gòu)變形能力從而判斷結(jié)構(gòu)是否滿足建筑功能要求的指標是合理的。因此Zou等[1]建立以結(jié)構(gòu)質(zhì)量最小化為目標,同時滿足層間相對位移和設(shè)計變量約束的優(yōu)化數(shù)學模型,分別以振型分解反應譜法和時程分析法為基礎(chǔ)推導出抗震優(yōu)化設(shè)計準則,根據(jù)優(yōu)化設(shè)計準則提出優(yōu)化設(shè)計方法,是一種有效和實用的抗震優(yōu)化設(shè)計方法,但該方法也存在不足,如只要求地震響應過程中層間相對位移在某些時點(加速度峰值及鄰近的一些時點)滿足約束條件,但實際上,當在優(yōu)化過程中出現(xiàn)比較柔的結(jié)構(gòu)時,層間相對位移響應的最大值并不一定出現(xiàn)在加速度峰值及鄰近的一些時點,故該方法得到的優(yōu)化設(shè)計在地震響應過程中不一定能滿足約束條件,可能得到不可行的解。Lagarosa等[2]以時程分析法為基礎(chǔ),并考慮材料的非線性,采用進化算法求解了以結(jié)構(gòu)質(zhì)量最小化為目標,同時滿足層間相對位移和設(shè)計變量約束的優(yōu)化數(shù)學模型,但該方法由于采用隨機搜索方法,存在計算工作量過大的缺點。黃冀卓等[3]以鋼框架質(zhì)量最輕和結(jié)構(gòu)總動應變能最小為目標,基于相關(guān)的設(shè)計規(guī)范,給出了抗震鋼框架多目標優(yōu)化問題的一種合理提法,采用Pareto遺傳算法求解抗震鋼框架多目標優(yōu)化問題,得到多目標優(yōu)化問題的妥協(xié)解。結(jié)構(gòu)抗震優(yōu)化設(shè)計屬于結(jié)構(gòu)動力優(yōu)化設(shè)計領(lǐng)域,結(jié)構(gòu)動力優(yōu)化設(shè)計方法從數(shù)學的角度看可分為零階優(yōu)化方法(非梯度類算法)[4]、一階優(yōu)化方法[5]和二階優(yōu)化方法。零階優(yōu)化方法不需要計算結(jié)構(gòu)動力響應對設(shè)計變量的一階導數(shù),如粒子群方法[6]、遺傳算法[7-8]和模擬退火算法[9]。一階優(yōu)化方法需要計算結(jié)構(gòu)動力響應對設(shè)計變量的一階導數(shù)如共軛梯度法[10],一階優(yōu)化方法關(guān)鍵是計算結(jié)構(gòu)動力響應對設(shè)計變量的一階導數(shù)(也稱敏度分析),國內(nèi)外的學者已經(jīng)發(fā)展了多種敏度分析方法[11-12]。二階優(yōu)化方法不僅需要計算結(jié)構(gòu)動力響應對設(shè)計變量的一階導數(shù),而且需要計算結(jié)構(gòu)動力響應對設(shè)計變量的二階導數(shù)(也稱Hessian矩陣分析),如牛頓類算法[13],結(jié)構(gòu)動力響應對設(shè)計變量的二階導數(shù)的計算難度非常大,因此在國內(nèi)外很少見到報道。通常,二階優(yōu)化方法的計算效率高于一階優(yōu)化方法,一階優(yōu)化方法的計算效率高于零階設(shè)計方法。

        本文的主要工作是Newmark-β法的基礎(chǔ)上推導了層間相對位移對設(shè)計變量一階導數(shù)和二階導數(shù)的計算公式,建立結(jié)構(gòu)抗震優(yōu)化設(shè)計的數(shù)學模型,構(gòu)造求解結(jié)構(gòu)抗震優(yōu)化設(shè)計問題二階優(yōu)化算法。最后演示了一個三層兩跨的平面框架進行抗震優(yōu)化設(shè)計,優(yōu)化結(jié)果表明本文的抗震優(yōu)化設(shè)計方法能獲得結(jié)構(gòu)的最優(yōu)設(shè)計。

        1 優(yōu)化數(shù)學模型

        框架具有i=1,2,3,…,N個單元組,單元的橫截面為矩形,如圖1所示,橫截面的寬度(be)和高度(he)為設(shè)計變量。

        為了討論的方便,設(shè)計變量向量定義為:

        框架抗震優(yōu)化設(shè)計的目標函數(shù)為結(jié)構(gòu)的質(zhì)量,由設(shè)計變量表示為:

        其中:ρ為材料的密度;li為第i個單元組中單元的長度;ni為第i個單元組所含單元的數(shù)目。

        圖1 平面梁單元Fig.1 Plane beam element

        對高層建筑結(jié)構(gòu)而言,結(jié)構(gòu)的側(cè)移通常成為設(shè)計的控制指標。高層建筑結(jié)構(gòu)抗震設(shè)計必須限制結(jié)構(gòu)的層間相對位移,因為過大的層間相對位移會使主體結(jié)構(gòu)出現(xiàn)裂縫甚至破損,使結(jié)構(gòu)產(chǎn)生附加內(nèi)力,嚴重使會加速結(jié)構(gòu)的倒塌,使填充墻及建筑裝修出現(xiàn)裂縫或破損,使電梯軌道變形過大,使人不舒服,影響使用。《建筑抗震設(shè)計規(guī)范》[14]規(guī)定,以彈性方法計算的層間相對位移應滿足以下條件:

        式中,ΔuJ(t)為在t時刻第J層的彈性層間相對位移;uJ(t)和uJ-1(t)分別為在t時刻相鄰的第J和J-1樓層的水平位移;HJ為第J層的高度;θJ為抗震設(shè)計規(guī)范規(guī)定的第J層層間轉(zhuǎn)角;TT為地震的持續(xù)時間。

        在抗震優(yōu)化設(shè)計中,除了考慮層間相對位移的約束外,構(gòu)件的強度和延性約束也是非常重要的,但由于問題本身很復雜,在本文中只考慮層間相對位移的約束(假設(shè)結(jié)構(gòu)的側(cè)移是結(jié)構(gòu)設(shè)計的控制指標),但文中的優(yōu)化設(shè)計方法將可推廣到考慮構(gòu)件的強度和延性約束、同時考慮受到靜力和動力約束的情況,抗震優(yōu)化設(shè)計問題通過式(1)~式(4)描述為:

        本文將采用高效的牛頓類優(yōu)化算法求解優(yōu)化數(shù)學模型式(5),優(yōu)化數(shù)學模型可通過積分型罰函數(shù)轉(zhuǎn)變?yōu)橐幌盗袩o約束優(yōu)化設(shè)計問題,然后計算罰函數(shù)對設(shè)計變量的梯度和海森矩陣,從而可構(gòu)造高效的牛頓類優(yōu)化算法。由式(5)可見,要計算罰函數(shù)對設(shè)計變量的梯度和海森矩陣,必須計算層間相對位移和質(zhì)量對設(shè)計變量的一階和二階導數(shù)。質(zhì)量對設(shè)計變量的一階和二階導數(shù)的計算是很容易解決的問題,因此層間相對位移對設(shè)計變量的一階和二階導數(shù)的計算是解決問題的關(guān)鍵。

        2 一、二階導數(shù)計算

        2.1 層間位移及其對設(shè)計變量的一、二階導數(shù)

        2.1.1 層間相對位移

        初始條件:

        其中:M、C、K分別為結(jié)構(gòu)的質(zhì)量矩陣、阻尼矩陣和剛度矩陣;x(t)、(t)、(t)分別為相對于地面的結(jié)點位移、速度和加速度;Eu為一向量,在節(jié)點線位移u方向、v方向和節(jié)點角位移θ方向上分別為1、0、0。

        結(jié)構(gòu)的阻尼矩陣采用瑞利阻尼如下:

        其中,

        式中,ω1、ω2分別為結(jié)構(gòu)的第一、二階固有頻率;ζ1、ζ2為阻尼比,本文取 ζ1=ζ2=0.05。

        將采用Newmark-β法求解(6)、式(7),框架在地震作用下層間相對位移響應的計算步驟如下:

        (1)初始計算

        其中,積分參數(shù)δ≥0.5,β=0.25(0.5+δ)2;⑤ 計算有效的剛度矩陣K*:

        (2)計算t+Δt時刻的響應

        ① 計算t+Δt時刻有效的荷載向量F*:

        ② 計算t+Δt時刻的節(jié)點位移x(t+Δt):

        ③ 計算t+Δt時刻的節(jié)點加速度(t+Δt):

        ⑤ 用式(4)計算t+Δt時刻的層間相對位移ΔuJ(t+Δt);

        (3)重復(2)步驟①~⑤即可獲得各個時刻的層間相對位移。

        2.1.2層間位移對設(shè)計變量的一階和二階導數(shù)

        式(13)兩邊分別對設(shè)計變量di求導一次,得:

        其中,

        式(17)和式(18)分別通過式(11)和式(12)對設(shè)計變量di求導一次獲得。

        通過式(16)計算t+Δt時刻位移對設(shè)計變量di的一階導數(shù)后,式(14)兩邊分別對設(shè)計變量di求導一次,得:

        通過式(19)計算t+Δt時刻加速度對設(shè)計變量di的一階導數(shù)后,式(15)兩邊分別對設(shè)計變量di求導一次,得:

        式(4)兩邊分別對設(shè)計變量di求導一次,得到層間相對位移對設(shè)計變量的一階導數(shù):

        式(16)兩邊分別對設(shè)計變量dj進一步求導,得:

        其中:

        式(23)和式(24)分別由式(17)和式(18)對設(shè)計變量di求導一次獲得。

        由式(22)求得?2x(t+ Δt)/(?di?dj)后,式(19)兩邊分別對設(shè)計變量dj求導,得:

        式(21)兩邊分別對設(shè)計變量dj求導,得:

        (1)初始計算

        其中,積分參數(shù)δ≥0.5,β=0.25(0.5+δ)2;

        ⑥ 由式(11)計算K*。

        (2)計算t+Δt時刻的響應

        ① 由式(13)計算t+Δt時刻的節(jié)點位移x(t+Δt);

        ② 由式(14)計算t+Δt時刻 的節(jié)點加速度(t+Δt);

        ③ 由式(15)計算(t+Δt)時刻的節(jié)點速度(t+Δt);

        ④由式(4)計算t+Δt時刻的層間相對位移ΔuJ(t+Δt);

        ⑤ 由式(16)計算t+Δt時刻的節(jié)點位移對設(shè)計變量的一階導數(shù)?x(t+Δt)/?di;

        ⑦ 由式(20)計算t+Δt時刻的節(jié)點速度對設(shè)計變量的一階導數(shù)(t+Δt)/?di;

        ⑧ 由式(21)計算t+Δt時刻的層間相對位移對設(shè)計變量的一階導數(shù)?ΔuJ(t+Δt)/?di;

        ⑨ 由式(22)計算t+Δt時刻的節(jié)點位移對設(shè)計變量的二階導數(shù)?2x(t+ Δt)/(?di?dj);

        ⑩ 由式(25)計算t+Δt時刻的節(jié)點加速度對設(shè)計變量的二階導數(shù)?2(t+Δt)/(?di?dj);

        (11) 由式(26)計算t+Δt時刻的節(jié)點速度對設(shè)計變量的二階導數(shù)?2(t+Δt)/(?di?dj);

        (12) 由式(27)計算t+Δt時刻的層間相對位移對設(shè)計變量的二階導數(shù)?2ΔuJ(t+Δt)/(?di?dj)。

        (3)重復執(zhí)行①~(12)步即可求得地震過程中各個時刻的響應。

        2.2 結(jié)構(gòu)質(zhì)量對設(shè)計變量的一階和二階導數(shù)

        式(2)兩邊分別對設(shè)計變量求導一次得結(jié)構(gòu)質(zhì)量對設(shè)計變量的一階導數(shù)如下:

        式(28)兩邊分別對設(shè)計變量進一步求導得結(jié)構(gòu)質(zhì)量對設(shè)計變量的二階導數(shù)如下:

        3 優(yōu)化數(shù)學模型的求解

        將數(shù)學模型(5)的約束條件規(guī)則化,得到等效的數(shù)學模型如下:

        3.1 優(yōu)化設(shè)計問題的轉(zhuǎn)變

        不等式優(yōu)化數(shù)學模型(30)可通過罰函數(shù)轉(zhuǎn)變?yōu)橐幌盗袩o約束的積分形式的優(yōu)化數(shù)學模型,為了求導的方便本文采用內(nèi)點罰函數(shù),如下:

        其中:罰系數(shù)rk是一遞減數(shù)列,當rk→0,罰函數(shù)P(d,rk)的最小值逼近約束優(yōu)化問題的最小值。因此,求解不等式約束優(yōu)化問題(30)轉(zhuǎn)變?yōu)榍蠼庖幌盗袩o約束優(yōu)化問題如下:

        罰系數(shù)的初始值r1由下式計算:

        其中:d0為初始設(shè)計點,p0=1~50。在本文中p0=10。

        罰系數(shù)rk根據(jù)下式遞減:

        其中:c=10~50,本文取c=10;k為搜索過程中使用罰函數(shù)的個數(shù)。

        3.2 罰函數(shù)的梯度和海森矩陣

        罰函數(shù)的梯度和海森矩陣即罰函數(shù)對設(shè)計變量的一階和二階導數(shù),式(31)兩邊分別對設(shè)計變量di求導一次得罰函數(shù)對設(shè)計變量一階導數(shù)如下:

        式(35)兩邊對設(shè)計變量dj進一步求導得罰函數(shù)對設(shè)計變量的二階導數(shù)如下:

        通過式(35)和式(36)計算罰函數(shù)對設(shè)計變量的一階和二階導數(shù)后即可獲得罰函數(shù)的梯度和海森矩陣。

        3.3 優(yōu)化設(shè)計問題的求解

        利用梯度和海森矩陣信息,構(gòu)造求解不等式約束優(yōu)化數(shù)學模型(30)的優(yōu)化設(shè)計算法,其步驟如下:

        第1步 選擇可行的初始設(shè)計點d0,由式(33)計算罰系數(shù)的初始值,定義收斂容差ε1=10,令k=1;

        第2步 從設(shè)計點dk-1出發(fā),采用 Marquardt方法[15]求解優(yōu)化數(shù)學模型(32)獲得最優(yōu)設(shè)計dk;用Marquardt方法求解優(yōu)化數(shù)學模型(32)的步驟:

        ② 令i=0,λ(0)=105;

        ⑤i≥MI?是,執(zhí)行第(11)步;否,執(zhí)行第⑥步;

        ⑥計算

        ⑩ 令λ(i)=2λ(i),返回第⑥步;

        第3 步P(dk-1,rk)-P(dk,rk)≤ε1?是,dk為最優(yōu)設(shè)計,輸出結(jié)果,停止計算;否,執(zhí)行第④步;

        第4步 由式(34)計算罰系數(shù)rk+1,k=k+1,回到第②步。

        4 數(shù)值算例

        一個三層兩跨的平面框架如圖2所示,假設(shè)每個樓層的質(zhì)量為30 t,框架材料的彈性模量為20 GPa,阻尼比為5%,梁和柱的截面尺寸空間見圖2的上方,層間轉(zhuǎn)角限值為1/450,則第一層層間相對位移的限值為0.012 m,第二、三層的層間相對位移限值為0.008 m。

        支座1,5,9受到El Centro水平地震加速度時程的激勵,El Centro水平加速度峰值為0.349 g,記錄時長為 8 s,步長 Δt=0.02 s。

        從初始可行設(shè)計出發(fā),通過4.3節(jié)的優(yōu)化設(shè)計算法獲得一系列優(yōu)化設(shè)計如表1所示。

        由表1和表2可見,隨著罰系數(shù)逐步趨向于0,目標函數(shù)即結(jié)構(gòu)的質(zhì)量越來越小,優(yōu)化設(shè)計逐步趨向于最優(yōu)設(shè)計(r4=0.327和r5=0.032 7獲得的優(yōu)化設(shè)計),收斂時前后兩次最優(yōu)設(shè)計的質(zhì)量不變,說明收斂精度非常高。初始可行設(shè)計的質(zhì)量為156 735 kg,最優(yōu)設(shè)計的質(zhì)量為57 043 kg,結(jié)構(gòu)的質(zhì)量減輕了99 692 kg,最優(yōu)設(shè)計的質(zhì)量僅為初始可行設(shè)計的36.39%,優(yōu)化效果明顯。在優(yōu)化過程中僅僅用了5個罰系數(shù),說明優(yōu)化算法的效率非常高。由表2還可看出,即使僅使用1個罰系數(shù)(r1=327),得到的優(yōu)化設(shè)計也能得到極大的改善(質(zhì)量由156 735 kg下降為82 876 kg,優(yōu)化設(shè)計的質(zhì)量為初始可行設(shè)計的52.88%)。

        圖2 三層兩跨的平面框架Fig.2 A three-storey,two-bay planar frame

        表1 初始可行設(shè)計和優(yōu)化設(shè)計(單位:mm)Tab.1 Initial feasible design and optimal design(Units:mm)

        表2 初始可行設(shè)計和優(yōu)化設(shè)計的質(zhì)量(單位:kg)Tab.2 Weight of the initial feasible design and optimal design(Units:kg)

        初始可行設(shè)計和一系列優(yōu)化設(shè)計的質(zhì)量如表2所示。

        5 結(jié)論

        (1)提出一種計算層間相對位移對設(shè)計變量的一階導數(shù)和二階導數(shù)的方法,該方法只需要進行一次動力響應分析就可完成,具有非常高的計算效率。

        (2)提出框架的抗震優(yōu)化設(shè)計數(shù)學模型適用于結(jié)構(gòu)的側(cè)移是結(jié)構(gòu)設(shè)計的控制指標的情況(如高層建筑)。當還要考慮其他設(shè)計指標時,需要在優(yōu)化數(shù)學模型中增加響應的約束條件。

        (3)采用積分型罰函數(shù)將顯含時間參數(shù)的不等式約束優(yōu)化問題轉(zhuǎn)變?yōu)橐幌盗胁缓瑫r間參數(shù)的無約束優(yōu)化問題,并且提出積分型罰函數(shù)梯度和海森矩陣的計算方法。但由于采用內(nèi)點罰函數(shù),因此初始設(shè)計必須為可行的設(shè)計。

        (4)提出的優(yōu)化算法充分利用梯度和海森矩陣的信息確定搜索的方向和步長,不需要進行一維搜索,下降速度快,具有良好的收斂性,能獲得最優(yōu)設(shè)計。算例表明該優(yōu)化算法具有非常高的效率,并且在優(yōu)化過程中能獲得一系列可行的優(yōu)化設(shè)計供工程設(shè)計人員選用。

        (5)本文的抗震優(yōu)化設(shè)計方法雖然只輸入了一條地震波,但實際上可以同時輸入多條地震波,只需要在優(yōu)化數(shù)學模型中增加多個時程響應約束即可,計算方法基本不變。

        [1] Zou X K,Chan C M.An optimal resizing technique for seismic drift design of concrete buildings subjected to response spectrum and time history loadings[J].Computers and Structures,2005,83:1689-1704.

        [21]Lagarosa N D,F(xiàn)ragiadakisa M,Papadrakakisa M.Structural optimization:a tool for evaluating seismic design procedures[J].Engineering Structures,2006,28:1623-1633.

        [3]黃冀卓,王 湛.基于遺傳算法的抗震鋼框架多目標優(yōu)化設(shè)計[J].力學學報,2007,39(3):389-397.

        HUANG Ji-zhuo,WANG Zhan.Multiobjective optimization design of seismic steel frames using genetic algorithm[J].Chinese Journal of Theoretical and Applied Mechanics,2007,39(3):389-397.

        [4]Kocer F Y,Arora J S.Optimal design of h-frame transmission poles for earthquake loading[J].Journal of Structure Engineering,1999,125:1299-1308.

        [5] Arora J S,Cardoso J E B.A design sensitivity analysis principle and its implementation into ADINA[J].Computer and Structure,1989,32:691-705.

        [6] Pereza R E,Behdinan K.Particle swarm approach for structural design optimization[J]. Computers and Structures,2007,85:1579-1588.

        [7] Togan V,Daloglu A T.Optimization of 3d trusses with adaptive approach in genetic algorithms[J].Engineering Structures,2006,28:1019-1027.

        [8] Venanzi I,Materazzi A L.Multi-objective optimization of wind-excited structures[J].Engineering Structures,2007,29(6):983-990.

        [9] Lamberti L.An efficient simulated annealing algorithm for design optimization of truss structures[J].Computers and Structures,2008,86:1936-1953.

        [10]Fletcher R,Reeves C M.Function minimization by conjugate gradients[J].Computer Journal,1964,6:163-168.

        [11]Van Keulen F,Haftka R T,Kim N H.Review of options for structural design sensitivity analysis.part 1:linear systems[J]. ComputerMethodsin Applied Mechanicsand Engineering,2005,194:3213-3243.

        [12] Bogomolni M,Kirsch U,Sheinman I.Efficient design sensitivities of structures subjected to dynamic loading[J].International Journal of Solids and Structures,2006,43:5485-5500.

        [13] Kirsch U,Bogomolni M,Van Keulen F.Efficient finitedifference design-sensitivities[J].AIAA Journal,2005,43:399-405.

        [14] Broydon C G.Quasi-newton methods and their application to function minimization[J].Mathematics of Computation,1967,21:368-381.

        [15] GB50011-2001,建筑抗震設(shè)計規(guī)范[S].北京:中國建筑工業(yè)出版社,2001.

        [16] Reklaitis G V,Ravindran A,Ragsdell K M.Engineering optimization methods and applications[M].New York:John Wiley& Sons,1983.

        猜你喜歡
        一階二階層間
        求解Ericksen-Leslie方程的一階精度、線性穩(wěn)定的數(shù)值格式
        李超代數(shù)到Kac模的一階上同調(diào)
        基于超聲檢測的構(gòu)件層間粘接缺陷識別方法
        一類二階迭代泛函微分方程的周期解
        一類二階中立隨機偏微分方程的吸引集和擬不變集
        二階線性微分方程的解法
        基于層間接觸的鋼橋面鋪裝力學分析
        上海公路(2018年3期)2018-03-21 05:55:40
        一類二階中立隨機偏微分方程的吸引集和擬不變集
        結(jié)合多分辨率修正曲率配準的層間插值
        大段合采油井層間干擾主控因素研究
        无码熟妇人妻av影音先锋| 水蜜桃视频在线观看免费18| 亚洲日本精品一区二区三区| 中文字幕日韩人妻少妇毛片| 免费久久人人爽人人爽av| 日韩国产精品一区二区Hd| 国产精品一区成人亚洲| 成h视频在线观看免费| 国产精品久久久久高潮| 中文字幕亚洲无线码| 亚洲欧洲国无码| 亚洲精品国产综合久久 | 亚洲国产成人精品无码区在线秒播 | 99福利影院| 亚洲国产日韩一区二区三区四区| 男女猛烈拍拍拍无挡视频| 亚洲av无码乱码国产精品fc2| 被暴雨淋湿爆乳少妇正在播放| 精品福利一区二区三区| 人妻少妇久久久久久97人妻| 亚洲尺码电影av久久| 国产精品无套粉嫩白浆在线 | 免费国产黄网站在线观看| 91最新免费观看在线| 亚洲一区二区三区18| 凹凸国产熟女精品视频app| 久久这里只精品国产免费10| 成人国产精品高清在线观看| 日本免费久久高清视频| 男女啪动最猛动态图| 日本午夜福利| 亚洲精品国产二区在线观看| 久久不见久久见免费视频6 | 成人免费播放视频777777| 日本又黄又爽gif动态图| 成人无码激情视频在线观看| 国产精品成人av大片| 久久久久亚洲av无码专区首jn| 欧美国产日本精品一区二区三区 | 中文字幕有码在线亚洲| 国产精品亚洲а∨天堂2021|