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

基于粗粒度并行遺傳算法的隔震層參數(shù)優(yōu)化

2023-12-15 09:46:54劉全明賀一哲
地震工程學(xué)報 2023年6期
關(guān)鍵詞:優(yōu)化結(jié)構(gòu)

黨 育, 劉全明, 賀一哲

(蘭州理工大學(xué) 土木工程學(xué)院, 甘肅 蘭州 730050)

0 引言

如何選擇合理的隔震層參數(shù)和布置隔震支座是隔震結(jié)構(gòu)設(shè)計中的關(guān)鍵問題。目前工程師通常采用試算法,該方法可保證設(shè)計滿足要求,但無法保證設(shè)計為最優(yōu)。針對此問題,國內(nèi)外學(xué)者提出了各種優(yōu)化方法進(jìn)行隔震結(jié)構(gòu)設(shè)計。Fallah等[1]采用非支配排序多目標(biāo)遺傳算法NSGA-II算法,將隔震結(jié)構(gòu)簡化為剪切型,優(yōu)化得到滑移隔震結(jié)構(gòu)的隔震層質(zhì)量、摩擦系數(shù)和阻尼比;Fan等[2]將鉛芯橡膠支座系統(tǒng)等效為線性,采用二次規(guī)劃法得到隔震層的最優(yōu)屈服后剛度和屈服力;Nigdeli等[3]采用和聲搜索算法,在近場和遠(yuǎn)場地震下優(yōu)化得到隔震層的剛度和阻尼比;葉昆等[4]將隔震結(jié)構(gòu)簡化為等效線性的兩自由度結(jié)構(gòu),以水平向減震系數(shù)和隔震層最大位移的線性組合為目標(biāo)函數(shù),得到隔震層的最優(yōu)屈服強(qiáng)度比;洪紹文等[5]以隔震結(jié)構(gòu)的投入成本比最小作為目標(biāo)函數(shù),采用無導(dǎo)數(shù)信賴域算法優(yōu)化得到隔震支座的鉛芯和橡膠直徑。以上研究均將隔震結(jié)構(gòu)簡化為剪切型模型,隔震層簡化為用2~3個參數(shù)表達(dá)的線性或非線性模型。但由于對隔震結(jié)構(gòu)的力學(xué)模型進(jìn)行了較多簡化,影響了隔震結(jié)構(gòu)動力分析結(jié)果的準(zhǔn)確性,且得到的結(jié)果也只是隔震層最優(yōu)參數(shù),與設(shè)計最終要求的具體隔震支座布置尚有差距。由此,Dang等[6]針對隔震結(jié)構(gòu)建立了三維有限元模型,采用經(jīng)典遺傳算法結(jié)合整數(shù)規(guī)劃法,分兩階段完成隔震層參數(shù)優(yōu)化和隔震支座布置,既保證了隔震結(jié)構(gòu)動力分析的準(zhǔn)確性,也可直接優(yōu)化得到隔震支座布置。但由于用經(jīng)典遺傳算法進(jìn)行優(yōu)化時,種群個體要依次進(jìn)行時程分析,時程分析由單一CPU串行平臺的結(jié)構(gòu)動力分析軟件(SAP2000、ETABS)完成,整個優(yōu)化過程所需時間大致為每次時程分析時間的累加,導(dǎo)致優(yōu)化的耗時較長,限制了該方法的使用。

目前,提高計算效率的方法大致分為3種:采用更高性能的CPU;采用圖像處理器GPU[7];采用多個CPU的并行計算[8]。針對文獻(xiàn)[6]的問題,經(jīng)典遺傳算法中的各種群個體需要分別進(jìn)行時程分析,但實際上,各種群個體的時程分析計算具有天然的并行性,同時,隨著微處理技術(shù)的迅速發(fā)展,多核CPU的計算機(jī)已普及,若在多核CPU上實現(xiàn)并行遺傳算法,就無需提高計算機(jī)硬件(采用高性能計算機(jī))或修改計算平臺(采用基于GPU的有限元分析軟件)來提高隔震結(jié)構(gòu)的優(yōu)化效率。基于此,本文提出了基于粗粒度并行遺傳算法的隔震層參數(shù)優(yōu)化方法,即采用Python-ETABS交互使用的編程技術(shù),用Python編寫粗粒度并行遺傳算法,多進(jìn)程調(diào)用ETABS程序并進(jìn)行遺傳操作,實現(xiàn)基于粗粒度并行遺傳算法的隔震層參數(shù)優(yōu)化分析。最后,通過一個隔震工程的實例進(jìn)行驗證,證明本文提出的方法既可保證優(yōu)化結(jié)果的準(zhǔn)確性,又可提高計算效率。

1 隔震層參數(shù)優(yōu)化的并行模型

對于隔震層參數(shù)優(yōu)化問題,種群中的每個個體都要進(jìn)行動力分析和遺傳操作。由于動力分析需要調(diào)用ETABS完成,并行調(diào)用的ETABS數(shù)量越多,占用的計算資源就越大。為平衡計算效率和計算資源,研究選用了粗粒度并行模型,即將種群劃分為多個子種群,各子種群獨立并發(fā),執(zhí)行時程分析和遺傳操作。

隔震層參數(shù)優(yōu)化的粗粒度并行模型如圖1所示。其中,多進(jìn)程并行操作采用Python的Multiprocessing庫[9]來實現(xiàn),各進(jìn)程之間無交互,充分利用CPU的多核計算資源,提高計算效率。由于整個優(yōu)化過程需要遺傳算法和結(jié)構(gòu)動力分析相結(jié)合,遺傳算法用Python編程完成,而動力時程分析由ETABS實現(xiàn),故采用ETABS應(yīng)用程序接口ETABS API[10]完成Python與ETABS的相互調(diào)用。整個過程為:Python編程生成各子種群,各子種群獨立地用 ETABS API中的函數(shù)來打開對應(yīng)的隔震結(jié)構(gòu)模型,定義質(zhì)量源,設(shè)置荷載工況,設(shè)定節(jié)點約束,定義隔震支座參數(shù),進(jìn)行動力分析和輸出結(jié)果;再將結(jié)果傳遞至Python,完成子種群的遺傳操作,生成新的子種群;再次進(jìn)行子種群并行計算,從所有種群中挑選出最優(yōu)個體并記錄下來;重復(fù)以上過程,當(dāng)最優(yōu)個體達(dá)到迭代終止條件后,優(yōu)化程序結(jié)束。

圖1 隔震層參數(shù)優(yōu)化的粗粒度并行模型Fig.1 Coarse-grained parallel model for parameter optimization of isolation layer

2 基于粗粒度并行遺傳算法的隔震層參數(shù)優(yōu)化

2.1 種群編碼

隔震層參數(shù)優(yōu)化時,隔震層的力學(xué)模型可簡化為雙線性模型,如圖2所示。圖中,Keq為隔震層的水平等效剛度,Xo為隔震層的屈服位移,Xm為隔震層的極限位移。

圖2 隔震層的力學(xué)模型Fig.2 Mechanical model of isolation layer

隔震層的設(shè)計變量如式(1)所示:

(1)

粗粒度并行遺傳算法是將種群劃分為多個子種群,因此,染色體種群表示如式(2)所示:

Λ=[X1,X2,…,Xn]=[(Keq1,Xo1,Xm1),
(Keq2,Xo2,Xm2),…,(Keqn,Xon,Xmn)]

(2)

種群編碼采用浮點數(shù)編碼,原因是隔震層優(yōu)化參數(shù)為連續(xù)實數(shù)域且范圍較大,采用浮點式編碼可降低編碼長度,提高搜索速度。

2.2 適應(yīng)度評估

適應(yīng)度函數(shù)是根據(jù)目標(biāo)函數(shù)確定的用于區(qū)分種群中個體優(yōu)良程度的標(biāo)準(zhǔn),是遺傳算法進(jìn)行自然選擇的唯一依據(jù)。

在隔震層參數(shù)優(yōu)化設(shè)計中,采用的目標(biāo)函數(shù)如式(3)所示:

(3)

式中:β為設(shè)防地震下隔震結(jié)構(gòu)的水平向減震系數(shù);u為罕遇地震作用下隔震支座的最大水平位移;θ為罕遇地震作用下上部結(jié)構(gòu)的最大層間位移角;βmax為水平向減震系數(shù)限值,與非隔震結(jié)構(gòu)相比,隔震工程上部結(jié)構(gòu)的水平地震作用減少量不宜超過1/4,再考慮到支座性能偏差調(diào)整系數(shù)ψ=0.8[12],因此取βmax=0.25×0.8=0.2;umax為隔震支座的水平位移限值,依據(jù)《建筑抗震設(shè)計規(guī)范》[12]的要求,umax=min(0.55D,3tr),其中,D為支座的有效直徑,tr為支座內(nèi)部橡膠總厚度;θmax為上部結(jié)構(gòu)的層間位移角限值,依據(jù)《建筑隔震設(shè)計標(biāo)準(zhǔn)》[13],取θmax=1/150。y越小,則該染色體對應(yīng)的隔震層參數(shù)就越接近最優(yōu)值。

在遺傳算法中,越是優(yōu)良的個體適應(yīng)度越大,這恰恰與目標(biāo)函數(shù)相反,因此,必須將目標(biāo)函數(shù)進(jìn)行一定的轉(zhuǎn)換,將其變?yōu)檫m應(yīng)度函數(shù)的形式,本文采用的變換如式(4)所示:

fitness=-y

(4)

另外,在隔震層設(shè)計時,還需要滿足以下約束條件[6]:

(5)

在遺傳算法中,以上約束條件是通過懲罰策略對個體的適應(yīng)度進(jìn)行修正來完成的,即個體不滿足任意一個約束條件時,就將該個體對應(yīng)的目標(biāo)函數(shù)值擴(kuò)大數(shù)倍,意味著該個體的適應(yīng)度函數(shù)值將減小數(shù)倍,這樣,該個體將在進(jìn)化過程中被淘汰。

2.3 種群遷移

種群遷移使各子種群之間可傳遞優(yōu)良個體,避免優(yōu)化結(jié)果陷入局部最優(yōu)。種群遷移主要解決種群的遷移規(guī)模、遷移策略和遷移拓?fù)淙齻€問題[11]。本文采用單向環(huán)型遷移拓?fù)渥鳛楦髯臃N群間傳遞優(yōu)良個體的拓?fù)?遷移策略采用第i子種群的最優(yōu)個體來代替第i+1子種群的最差個體,實現(xiàn)種群間的信息交流。遷移規(guī)模由遷移周期和遷移率來決定,考慮到隔震層參數(shù)優(yōu)化問題的優(yōu)化參數(shù)較少,為平衡優(yōu)化解的質(zhì)量和優(yōu)化效率,取遷移周期為2代,遷移率為0.1。

2.4 遺傳操作

遺傳操作的多進(jìn)程并行模型如圖3所示。各種群被分配到不同的進(jìn)程,遺傳選擇采用精英保留策略,每個進(jìn)程根據(jù)適應(yīng)度函數(shù)選擇出種群中適應(yīng)度最大的個體保留下來。為保證種群子代與父代個體數(shù)相等,從父代種群中隨機(jī)選擇出剩余個體,按照設(shè)定的交叉率和變異率確定是否進(jìn)行交叉和變異操作:如果不交叉和變異,則子代與父代相同,否則產(chǎn)生新的子代,這樣反復(fù)進(jìn)行,直至達(dá)到終止條件。為保證種群個體的多樣性,本文取交叉率為[0.7,0.9]的隨機(jī)數(shù),變異率為[0.01,0.05]的隨機(jī)數(shù)。

圖3 遺傳操作并行模型Fig.3 Parallel model of genetic operation

2.5 隔震層參數(shù)優(yōu)化設(shè)計流程

基于粗粒度并行遺傳算法的隔震層參數(shù)優(yōu)化流程如圖4所示,具體步驟為:

圖4 基于粗粒度并行遺傳算法的隔震層參數(shù)優(yōu)化計算流程Fig.4 Optimization calculation process of isolation layer parameters based on coarse-grained parallel genetic algorithm

(1) 預(yù)估隔震支座,給出設(shè)計變量Keq、Xo、Xm的取值范圍;

(2) 給定粗粒度并行遺傳算法的相關(guān)參數(shù),包括子種群數(shù)量、染色體個數(shù)、交叉率、變異率、遷移周期、遷移率和最大迭代次數(shù);

(3) 確定種群編碼,生成初始種群;

(4) 并行計算出各種群個體的目標(biāo)函數(shù),評估個體的適應(yīng)度;

(5) 判斷迭代代數(shù)N是否為遷移周期的K倍,如果是,則執(zhí)行遷移操作;

(6) 各子種群并行執(zhí)行遺傳操作,形成新的種群,計算新種群個體的目標(biāo)函數(shù)值;

(7) 判斷是否滿足迭代終止條件:如果是,輸出隔震層最優(yōu)參數(shù)及對應(yīng)的目標(biāo)函數(shù)值;如果否,令N=N+1,并轉(zhuǎn)到步驟(5)。

3 工程實例及優(yōu)化分析

3.1 隔震工程概況

以某實際隔震工程為例進(jìn)行分析,該工程為乙類建筑,地下1層,地上4層,建筑總高度為19.0 m;設(shè)防烈度為8度(0.3g),場地類別為Ⅱ類,設(shè)計地震分組為第三組。隔震層設(shè)置在地下室與上部結(jié)構(gòu)之間,隔震結(jié)構(gòu)的有限元模型如圖5所示,隔震支座的原平面布置圖如圖6所示。

圖5 隔震結(jié)構(gòu)的有限元模型Fig.5 Finite element model of the isolation structure

圖6 隔震支座原平面布置圖(單位:mm)Fig.6 Original plan layout of the isolation bearings (Unit:mm)

考慮到設(shè)計結(jié)果要與優(yōu)化結(jié)果進(jìn)行對比,因此優(yōu)化時采用標(biāo)準(zhǔn)人工波[6]進(jìn)行計算,得到上部結(jié)構(gòu)的水平向減震系數(shù)為0.32,且罕遇地震下各隔震支座的最大位移不超過248 mm,小于支座容許位移275 mm[12];同時,罕遇地震下各隔震支座均受壓。由此可見,原設(shè)計的隔震支座布置方案滿足設(shè)計要求。

3.2 優(yōu)化過程及結(jié)果

對該工程采用粗粒度并行遺傳算法進(jìn)行隔震層參數(shù)設(shè)計。

首先,根據(jù)隔震支座在重力荷載代表值下的豎向壓應(yīng)力限制要求,預(yù)估隔震支座直徑范圍為400~700 mm,各直徑的支座可以采用天然橡膠支座(LNR)和鉛芯橡膠支座(LRB)兩種類型。各隔震支座的力學(xué)性能如表1所列。

表1 隔震支座的力學(xué)性能參數(shù)Table 1 Mechanical performance parameters of isolation bearings

由表1可知,若所有支座均為LNR400,隔震層參數(shù)最小,若所有支座均為LRB700,隔震層參數(shù)最大,由此得到隔震層參數(shù)范圍為:Keq∈[19.8,58.49]kN/mm,Xo∈[0,11.32]mm,Xm∈[68,140]mm。

目標(biāo)函數(shù)如式(6)所示:

(6)

約束條件如式(7)所示:

(7)

其中,由于初步選定的最小隔震支座直徑為400 mm,根據(jù)表1可知,直徑為400 mm的隔震支座橡膠總厚度為68.6 mm,則umax= min(400×0.55,3×68.6) mm=205.8 mm。

將以上模型采用粗粒度遺傳算法進(jìn)行優(yōu)化,選定子種群數(shù)為7,每個子種群的染色體個數(shù)為10,總迭代次數(shù)為30代。計算機(jī)處理器為Intel Xeon Gold 5115,共有10核,內(nèi)存為32 GB。時程分析時仍選用標(biāo)準(zhǔn)人工波,最終得到隔震層的最優(yōu)參數(shù)為:Keq=50.49 kN/mm,Xo=9.34 mm,Xm=122.75 mm。

(1) 優(yōu)化方案與原設(shè)計方案比較

根據(jù)已得到的隔震層最優(yōu)參數(shù),按照文獻(xiàn)[6]中的第二階段優(yōu)化方法,確定出隔震支座的最優(yōu)布置(圖7)。該方案由26個LRB600和4個LRB400組成,優(yōu)化方案與原設(shè)計方案的對比如表2~3所列。

表2 原設(shè)計方案與優(yōu)化方案對比Table 2 Comparison between the original design scheme and the optimization scheme

圖7 優(yōu)化后的隔震支座布置圖(單位:mm)Fig.7 Plan layout of the isolation bearings after optimization (Unit:mm)

從表2可看出,與原設(shè)計方案相比,采用優(yōu)化方案后,上部結(jié)構(gòu)的水平向減震系數(shù)由0.32減小至0.30,減小幅度為6.3%;罕遇地震下各隔震支座的最大位移由248 mm減小至210 mm,減小幅度為14.7%。優(yōu)化方案除滿足設(shè)計要求外,還可同時減小上部結(jié)構(gòu)的層剪力和隔震支座位移,減震效果優(yōu)于原設(shè)計方案。

表3為優(yōu)化方案與原設(shè)計方案上部結(jié)構(gòu)在罕遇地震下的最大層間位移角比較。從表3可看出,無論優(yōu)化方案還是原設(shè)計方案,頂層的層間位移角明顯小于其他層,這是由于該結(jié)構(gòu)的頂層質(zhì)量較小,使得頂層剪力小于其他層,但該層與其余層的抗側(cè)剛度差別不大。因此,結(jié)構(gòu)頂層的層間位移角遠(yuǎn)小于其他層,但與原設(shè)計方案相比,優(yōu)化方案的頂層層間位移角與第3層差別更小,說明優(yōu)化方案的結(jié)構(gòu)側(cè)向變形更為均勻。此外,除3層Y向兩方案相同外,優(yōu)化方案各層的最大層間位移角均比原設(shè)計方案小,說明優(yōu)化設(shè)計方案在罕遇地震下對上部結(jié)構(gòu)的保護(hù)要優(yōu)于原設(shè)計方案。

表3 原設(shè)計方案與優(yōu)化方案在罕遇地震下上部結(jié)構(gòu)的層間位移角對比Table 3 Comparison between story drift ratios of superstructure from the original and optimization schemes under rare earthquakes

(2) 與文獻(xiàn)[6]方法的計算效率對比

用粗粒度并行遺傳算法和文獻(xiàn)[6]采用的經(jīng)典遺傳算法方法分別對以上隔震結(jié)構(gòu)進(jìn)行優(yōu)化,設(shè)定兩種方法的進(jìn)化終止條件均為:連續(xù)7代中本代與上代的目標(biāo)函數(shù)值之差的絕對值都小于0.01,且最大迭代次數(shù)不超過30次。記錄兩種方法達(dá)到優(yōu)化終止條件時的目標(biāo)函數(shù)最優(yōu)值、進(jìn)化代數(shù)和所需時間,并列于表4。由于遺傳算法是一種隨機(jī)算法,對于不同初始種群、不同變異率和交叉率,其計算過程各異,因此,為使文獻(xiàn)[6]和粗粒度并行遺傳算法的計算過程具有可比性,對同一優(yōu)化問題分別計算了5次。

表4 兩種方法的計算耗時和最優(yōu)解對比Table 4 Comparison between calculation time and optimal solution of the two methods

由表4可看出,5次優(yōu)化分析中,采用經(jīng)典遺傳算法求得的最優(yōu)目標(biāo)函數(shù)值,有3次與全局最優(yōu)解相同,有2次未求得全局最優(yōu)解,并且每次優(yōu)化需要迭代25代左右才終止;采用粗粒度并行遺傳算法,5次均達(dá)到全局最優(yōu)解,且每次優(yōu)化僅需迭代18代左右就可終止。這說明采用粗粒度并行遺傳算法比經(jīng)典最優(yōu)算法更易促進(jìn)最優(yōu)解的出現(xiàn),并提高全局收斂速度。同時,采用粗粒度并行遺傳算法,每次優(yōu)化所需時間都明顯少于經(jīng)典遺傳算法,平均耗時分別為3.46 h和21.39 h,加速比大致為6,說明粗粒度并行遺傳算法不僅能準(zhǔn)確求解隔震層最優(yōu)參數(shù),還可大大提高優(yōu)化計算效率。

在5次計算中,分別選取兩方法中目標(biāo)函數(shù)最優(yōu)值最小且計算耗時最少的結(jié)果進(jìn)行比較,即取經(jīng)典遺傳算法第1次的計算結(jié)果和粗粒度并行遺傳算法第2次的結(jié)果,繪制出進(jìn)化次數(shù)和目標(biāo)函數(shù)的關(guān)系曲線如圖8所示。

圖8 兩種算法的進(jìn)化曲線Fig.8 Evolution curves of the two algorithms

從圖8可看出,經(jīng)典遺傳算法直至第15代后才基本趨于穩(wěn)定,整個進(jìn)化曲線呈階梯狀下降,而粗粒度并行遺傳算法運(yùn)行至第7代后已基本趨于穩(wěn)定,整個進(jìn)化曲線呈斷崖式下降,說明粗粒度遺傳算法可快速收斂至最優(yōu)解,優(yōu)化效率更高。

4 結(jié)論

針對基于經(jīng)典遺傳算法的隔震層參數(shù)優(yōu)化方法效率不高的問題,利用Python的多進(jìn)程機(jī)制和Python與ETABS的交互,將粗粒度并行遺傳算法應(yīng)用于隔震層參數(shù)優(yōu)化,并得到如下結(jié)論:

(1) 采用粗粒度并行遺傳算法優(yōu)化后的方案,除滿足設(shè)計要求外,還可同時減小上部結(jié)構(gòu)的層剪力和隔震支座位移,減震效果優(yōu)于原設(shè)計方案,說明該方法可有效進(jìn)行隔震層參數(shù)優(yōu)化。

(2) 用10核CPU計算,與經(jīng)典遺傳算法相比,基于粗粒度并行遺傳算法的隔震層參數(shù)優(yōu)化方法既能準(zhǔn)確求解出全局最優(yōu)解,又可顯著提高優(yōu)化效率,加速比約為6,可基本滿足隔震工程設(shè)計的及時性需求,具有較好的工程應(yīng)用價值。

猜你喜歡
優(yōu)化結(jié)構(gòu)
超限高層建筑結(jié)構(gòu)設(shè)計與優(yōu)化思考
《形而上學(xué)》△卷的結(jié)構(gòu)和位置
民用建筑防煙排煙設(shè)計優(yōu)化探討
關(guān)于優(yōu)化消防安全告知承諾的一些思考
一道優(yōu)化題的幾何解法
由“形”啟“數(shù)”優(yōu)化運(yùn)算——以2021年解析幾何高考題為例
論結(jié)構(gòu)
中華詩詞(2019年7期)2019-11-25 01:43:04
新型平衡塊結(jié)構(gòu)的應(yīng)用
模具制造(2019年3期)2019-06-06 02:10:54
論《日出》的結(jié)構(gòu)
創(chuàng)新治理結(jié)構(gòu)促進(jìn)中小企業(yè)持續(xù)成長
主站蜘蛛池模板: 麻豆国产精品| 一级毛片免费高清视频| 久久大香伊蕉在人线观看热2| 欧美第九页| 亚洲精品自拍区在线观看| 久操中文在线| 国产麻豆aⅴ精品无码| 日韩av电影一区二区三区四区| 日韩欧美网址| 亚洲91精品视频| 无码'专区第一页| 国产超薄肉色丝袜网站| 98超碰在线观看| 不卡无码h在线观看| 亚洲床戏一区| 亚洲有码在线播放| 国产在线无码av完整版在线观看| 亚洲成年人片| 日韩无码黄色网站| 美女裸体18禁网站| 欧美精品在线免费| 欧美www在线观看| 国产免费黄| 亚洲第一黄片大全| 91伊人国产| 色综合中文| 精品欧美一区二区三区久久久| 伊人久久福利中文字幕| 色综合五月婷婷| 成人午夜久久| 国产精品嫩草影院视频| 国产国产人成免费视频77777| 午夜小视频在线| 国产天天色| 日本伊人色综合网| 人妻一区二区三区无码精品一区| 性视频一区| 欧美福利在线观看| 秘书高跟黑色丝袜国产91在线 | 欧美日本不卡| www.99在线观看| 亚洲欧洲国产成人综合不卡| 日韩精品欧美国产在线| 97久久免费视频| 露脸国产精品自产在线播| 黄色三级毛片网站| 国产自产视频一区二区三区| 亚洲AV人人澡人人双人| 国内精品久久九九国产精品| 一级爆乳无码av| 国产精品内射视频| 精品国产网站| 91一级片| 国产精品无码一二三视频| 亚洲视频在线青青| 青青青国产精品国产精品美女| 国产精品制服| 无码日韩视频| 欧美日韩成人在线观看 | 狠狠操夜夜爽| 国产美女精品人人做人人爽| 欧美精品v欧洲精品| 狠狠亚洲五月天| 亚洲三级网站| 午夜日b视频| 国产精品伦视频观看免费| 日本一本在线视频| 欧美亚洲激情| 久久亚洲精少妇毛片午夜无码| 67194在线午夜亚洲 | 日本一本在线视频| 国模粉嫩小泬视频在线观看| 国产一区二区三区在线观看免费| 亚洲第一中文字幕| 亚洲精品天堂自在久久77| 久久人妻xunleige无码| 91免费国产在线观看尤物| 亚洲欧美日本国产综合在线 | 国产真实自在自线免费精品| 国产国拍精品视频免费看 | 亚洲精品国产乱码不卡| 婷婷色婷婷|