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

        ?

        不確定性簡諧激勵下連續(xù)體結構可靠性拓撲優(yōu)化

        2024-04-11 02:12:04時元昆
        振動與沖擊 2024年6期
        關鍵詞:振幅不確定性導數(shù)

        王 選, 時元昆, 陳 翔, 龍 凱

        (1. 合肥工業(yè)大學 土木與水利工程學院,合肥 230009; 2. 天津大學 機械工程學院,天津 300072;3. 華北電力大學 新能源電力系統(tǒng)國家重點實驗室,北京 102206)

        在航空航天、汽車、船舶等領域中,結構系統(tǒng)或機器中的振動和噪聲主要來源于旋轉部件引起的諧波力。為了保證結構系統(tǒng)或機械裝置的使用安全性和舒適性,通常對結構的動力學性能要求很高。一方面,要求結構的低階固有頻率遠離外激勵頻率,防止發(fā)生共振現(xiàn)象,造成結構破壞[1-2];另一方面,要求結構盡可能降低振動強度,將振動對裝置和人員的損害降到最低[3-4]。因此,為了改善結構的動力學性能,開展結構動力學拓撲優(yōu)化研究具有重要意義。早期關于動力學拓撲優(yōu)化問題的研究主要以結構的固有頻率為約束或目標,或者優(yōu)化給定階數(shù)的兩個相鄰高階特征頻率之間的間距等[5-6]。而考慮簡諧激勵下的動響應分析可幫助設計人員預測結構的持續(xù)性動力特性。簡諧激勵下的結構優(yōu)化設計對于抑制其振動響應具有重要意義。

        近年來,簡諧激勵下結構動響應拓撲優(yōu)化受到越來越多的關注。Ma等[7]基于均勻化方法研究了簡諧激勵下無阻尼結構的動柔度最小化問題。龍凱等[8]基于獨立連續(xù)映射方法建立了以結構體積最小化為目標,以關注自由度振幅平方和為約束的確定性優(yōu)化模型,并研究了激勵頻率、阻尼系數(shù)和激勵振幅對優(yōu)化結果的影響。Liu等[9]針對諧響應下大規(guī)模拓撲優(yōu)化問題,討論了模態(tài)位移法、模態(tài)加速度法和全分析方法的計算效率。房占鵬等[10]基于漸進結構優(yōu)化法實現(xiàn)了指定頻帶簡諧激勵下的約束阻尼結構布局優(yōu)化設計。Niu等[11]針對簡諧激勵下振動響應最小化問題,總結和比較了一些常用的目標函數(shù),包括動柔度、局部結構位移響應、結構加速度等。Zhao等[12]將模態(tài)疊加法與模型降階法相結合,研究了諧波激勵下具有比例阻尼的大尺度結構的拓撲優(yōu)化問題。徐家琪等[13]基于自然單元法研究了諧波激勵下結構動柔度最小化問題。Long等[14]基于二次規(guī)劃算法實現(xiàn)了簡諧激勵下應力約束拓撲優(yōu)化問題的求解。Montero等[15]提出了一種基于密度加權范數(shù)的目標函數(shù)解決連續(xù)體結構在諧波激勵下的拓撲優(yōu)化問題。Zhao等[16]提出了一種自適應混合展開法并驗證了其在諧波激勵下求解拓撲優(yōu)化問題的有效性。Wang等[17]采用拓撲優(yōu)化方法實現(xiàn)了諧響應下結構的保形拓撲優(yōu)化設計。Wang等[18]提出了一種結合二階Krylov子空間和多重網(wǎng)格法的降階方法,用以求解簡諧激勵下多相材料結構的拓撲優(yōu)化設計問題。

        不難看出,上述簡諧激勵下結構動響應拓撲優(yōu)化的研究都是在確定性條件下開展的[19],忽略了簡諧激勵振幅和頻率等不確定性因素對優(yōu)化結果的影響,可能會導致結構動力學性能對不確定因素的影響過于敏感,結構的可靠性水平較低。目前考慮不確定性的結構拓撲優(yōu)化可分為兩類:魯棒性拓撲優(yōu)化[20]和可靠性拓撲優(yōu)化(reliability-based topology optimization, RBTO)[21-22]。魯棒性拓撲優(yōu)化旨在增強結構抵抗不確定性因素干擾的能力??煽啃酝負鋬?yōu)化旨在設計滿足指定可靠性水平的結構。目前已有一些學者嘗試將考慮不確定性因素的魯棒性設計理論整合到諧響應拓撲優(yōu)化中,以提高動響應結構抵抗不確定性激勵干擾的能力。Zhang等[23]研究了未知但有界的簡諧激勵下的結構魯棒拓撲優(yōu)化問題,其中優(yōu)化目標是最小化最壞情況的動柔度。王棟[24]基于變密度方法研究了簡諧激勵作用位置不確定性情況下結構動態(tài)魯棒性拓撲優(yōu)化設計,以降低結構對載荷作用點隨機擾動的影響。Valentini等[25]采用分層抽樣的蒙特卡羅模擬方法對結構動力響應的期望值和標準差進行建模,研究了考慮激勵頻率不確定性的最小化結構動態(tài)響應的魯棒性設計。

        與上述魯棒性拓撲優(yōu)化工作不同,本文旨在將考慮不確定性的可靠性分析理論整合到諧響應拓撲優(yōu)化問題中,針對不確定性簡諧激勵下連續(xù)體結構的可靠性設計問題,提出了一種考慮簡諧載荷振幅大小和頻率不確定性的可靠性拓撲優(yōu)化方法,以設計滿足指定可靠性水平的動響應結構,確保結構在不確定性簡諧激勵下的失效概率小于或者等于指定的失效概率。

        1 簡諧激勵下確定性拓撲優(yōu)化

        1.1 諧響應有限元分析

        簡諧激勵下結構的動力學有限元方程為

        (1)

        結構整體剛度矩陣和質量矩陣可分別通過組裝對應的單元剛度矩陣和單元質量矩陣得到,其表達式為

        (2)

        本文采用工程中常見的比例阻尼,故阻尼矩陣可表示為剛度矩陣和質量矩陣的線性組合

        C=cmM+ckK

        (3)

        式中,cm和ck為比例阻尼系數(shù)。

        令U=ueiωt并將其代入式(1)可得

        (K-ω2M+iωC)u=F

        (4)

        進一步簡化式(4)可得

        Kdu=F

        (5)

        式中:Kd為阻尼系統(tǒng)的等效剛度矩陣,即動剛度矩陣,Kd=K-ω2M+iωC; 這里位移振幅u為一個復向量。

        1.2 簡諧激勵下確定性拓撲優(yōu)化模型

        由于簡諧激勵下結構的振幅反映了結構的振動強度,為了抑制結構在簡諧激勵下的振動行為,以所關注自由度對應的振幅平方和為約束,以結構相對于初始狀態(tài)的體積比為目標,建立簡諧激勵下確定性拓撲優(yōu)化(deterministic topology optimization,DTO)模型為

        (6)

        2 不確定性簡諧激勵下可靠性拓撲優(yōu)化

        2.1 簡諧激勵下可靠性拓撲優(yōu)化模型

        目前關于簡諧激勵下連續(xù)體結構拓撲優(yōu)化的研究基本是在確定性條件下進行的,忽略了簡諧激勵的振幅和頻率的不確定性對優(yōu)化結果的影響。為了提高諧響應下結構的可靠性,設計滿足指定可靠性水平的結構,本文在諧響應拓撲優(yōu)化中考慮載荷振幅和頻率的不確定性,建立了不確定性簡諧激勵下連續(xù)體結構可靠性拓撲優(yōu)化模型,其數(shù)學列式為

        (7)

        式中:X為不確定性簡諧激勵的振幅大小或頻率隨機變量;F(X)為隨機簡諧激勵的振幅函數(shù); 需要強調的是在式(6)定義的確定性優(yōu)化模型中,結構的動剛度矩陣Kd和位移振幅u僅與設計變量ρ相關,而在式(7)定義的可靠性優(yōu)化模型中,結構的動剛度矩陣Kd和位移振幅u為設計變量ρ和隨機變量X的函數(shù),受其二者影響;G(ρ,X)為用來度量結構可靠性的極限狀態(tài)函數(shù),定義為G(ρ,X)=Alimit-A,G(ρ,X)≤0表示結構失效;Pf為結構的失效概率,要求其小于或等于指定的失效概率Pf,t??煽啃阅P椭衅渌麉?shù)均與式(6)定義的確定性優(yōu)化模型相同。

        2.2 基于PMA的可靠性分析

        在實際問題中,基于蒙特卡羅的方法計算失效概率的成本較大,為此,一些基于Taylor展開的數(shù)值近似方法相繼被提出來處理式(7)中的可靠性約束,包括功能度量法(performance measure approach, PMA)[26-27]和可靠性指標法(reliability index approach, RIA)[28-29]。與RIA相比,PMA具有更高的效率和穩(wěn)定性,因此本文使用PMA?;赑MA的諧響應可靠性拓撲優(yōu)化列式表示為

        (8)

        式中,Gp(ρ,X)為極限狀態(tài)函數(shù)在目標可靠性曲面上的最可能失效點(most probable point, MPP)處計算的最小概率性能值,可以通過求解以下優(yōu)化模型得到

        min:G(ρ,Y)
        s.t. ‖Y‖=βt

        (9)

        式中:Y為由隨機變量X轉化成的標準正態(tài)分布的隨機變量;βt為目標可靠性指標,其與失效概率Pf,t之間關系為βt=-Φ-1(Pf,t),Φ-1為標準正態(tài)累積分布函數(shù)的反函數(shù)。通過求解式(9)得到的最可能失效點,記為Yp,則Gp=G(ρ,Yp)。

        本文采用先進均值(advanced mean value, AMV)法查找求解式(9),尋找最可能失效點,其迭代格式為

        (10)

        式中:YAMV,k+1為內層循環(huán)中使用AMV方法在第k+1次迭代中計算的隨機變量值; ?YG(ρ,YAMV,k)為極限狀態(tài)函數(shù)關于隨機變量Y的導數(shù);n(YAMV,k)為G(ρ,Y)在YAMV,k處的歸一化方向。

        可以看出本文考慮簡諧載荷振幅和頻率不確定性的諧響應可靠性拓撲優(yōu)化是嵌套的雙循環(huán)優(yōu)化算法,內層循環(huán)旨在尋找MPP,計算最小概率性能值Gp;外層循環(huán)更新設計變量ρ,形成新的結構。

        3 靈敏度分析

        由式(10)可知,尋找MPP需要計算極限狀態(tài)函數(shù)對隨機變量的導數(shù)。另外,為了能使用基于梯度的優(yōu)化算法求解外層優(yōu)化問題,需要計算極限狀態(tài)函數(shù)對設計變量的導數(shù),確定優(yōu)化搜索方向。因此,本章實施靈敏度分析,詳細推導極限狀態(tài)函數(shù)關于設計變量和隨機變量的導數(shù)。

        3.1 極限狀態(tài)函數(shù)對設計變量的導數(shù)

        極限狀態(tài)函數(shù)關于設計變量的導數(shù)可由鏈式法則推導得到

        (11)

        式中,導數(shù)項?G/?A可通過G(ρ,X)=Alimit-A直接求導得到

        (12)

        (13)

        設隨機簡諧激勵的振幅F與設計變量無關,式(5)兩邊同時對設計變量ρe求導可得

        (14)

        式(13)中導數(shù)項?u/?ρe可表示為

        (15)

        將式(15)代入式(13)可得

        (16)

        令伴隨向量λ為式(17)定義的伴隨方程的解

        (17)

        (18)

        基于單元密度的拓撲優(yōu)化方法通常會出現(xiàn)棋盤格現(xiàn)象。棋盤格現(xiàn)象是指單元密度為0或1的單元在設計域內交錯布置的現(xiàn)象,使得優(yōu)化后的拓撲結構的可制造性差[30]。為了避免棋盤格現(xiàn)象,獲得清晰的黑白設計,本文采用以下敏度過濾方法對極限狀態(tài)函數(shù)關于設計變量的導數(shù)進行過濾

        (19)

        w(xj)=max(R-‖xe-xj‖,0)

        (20)

        式中,xe和xj分別為單元e和單元j的中心坐標。

        3.2 極限狀態(tài)函數(shù)對隨機變量的導數(shù)

        極限狀態(tài)函數(shù)關于激勵振幅隨機變量z1和頻率隨機變量z2的導數(shù)由鏈式法則推導為

        (21)

        與導數(shù)項?A/?ρe的推導過程類似,導數(shù)項?A/?z1和?A/?z2可分別表示為

        (22)

        由于動剛度Kd與激勵振幅隨機變量z1無關, 載荷F與頻率隨機變量z2無關,式(5)兩邊分別對隨機變量z1和z2求導可得

        (23)

        則式(22)中導數(shù)項?u/?z1和?u/?z2可分別表示為

        (24)

        將式(24)代入式(22)可得

        (25)

        同樣令伴隨向量λ為式(26)定義的伴隨方程的解

        (26)

        (27)

        本文采用Svanberg[31]提出的移動漸進線方法來求解式(6)和式(8)定義的優(yōu)化問題。當內層迭代過程滿足|Gk-Gk-1|≤10-6(Gk為極限狀態(tài)函數(shù)在內層循環(huán)中第k次迭代的計算值)收斂條件時,內層迭代過程結束,進入外層拓撲優(yōu)化,每一個拓撲優(yōu)化迭代步都要進行內層逆可靠性分析循環(huán),當外層拓撲優(yōu)化達到最大迭代步數(shù)時,程序終止運行,本文算法的實施流程如圖1所示。

        圖1 本文算法的流程圖Fig.1 Flowchart of the proposed algorithm

        4 數(shù)值算例與討論

        本章通過兩個算例來驗證所提出的考慮簡諧激勵不確定性的諧響應可靠性拓撲優(yōu)化方法的有效性。所有算例均使用相同的材料參數(shù),假設材料彈性模量、泊松比和質量密度分別為210 GPa、0.3和7 800 kg/m3,阻尼系數(shù)cm=ck=1×10-4。簡諧激勵的振幅大小和頻率均為服從正態(tài)分布的隨機變量。設置所有設計變量ρe的初始值為1,最大迭代步數(shù)為300。

        4.1 懸臂梁

        懸臂梁結構尺寸如圖2所示,設計域離散為240×80個邊長為5 mm、厚度為1 mm的正四邊形單元,左邊界完全固支,右邊界中部受豎直向下的隨機簡諧激勵作用,振幅大小的平均值和標準差分別為1 000 N和100 N,頻率的平均值和標準差分別為90 Hz和9 Hz。本算例將簡諧載荷均勻分配到中間相鄰的3個有限元節(jié)點上,以載荷作用點豎直方向的振幅平方和為極限狀態(tài)函數(shù)來構造可靠性約束,設置約束限值為4.7 mm,過濾半徑設為25 mm,目標可靠性指標為βt=2。

        圖2 懸臂梁結構示意圖(mm)Fig.2 Illustration of cantilever beam(mm)

        為了研究簡諧激勵的振幅大小和頻率的不確定性對諧響應拓撲優(yōu)化結果的影響,討論以下4種不同情況。工況1:確定性載荷振幅大小和頻率下的DTO結果,如圖3所示,此時載荷振幅大小和頻率分別為1 000 N和90 Hz。工況2:載荷振幅大小確定(取其平均值1 000 N),頻率不確定下的RBTO結果,如圖4所示。工況3:載荷頻率確定(取其平均值90 Hz),載荷振幅大小不確定下的RBTO結果,如圖5所示。工況4:載荷振幅大小和頻率均不確定下的RBTO結果,如圖6所示。

        圖3 確定性載荷振幅和頻率下的DTO設計Fig.3 DTO design for deterministic load amplitude and frequency

        圖5 確定性載荷頻率、不確定性載荷振幅下的RBTO設計Fig.5 RBTO design for deterministic load frequency and uncertain load amplitude

        圖6 不確定性載荷振幅和頻率下的RBTO設計Fig.6 RBTO design for uncertain load amplitude and frequency

        為了驗證BRTO設計相對于DTO設計的優(yōu)勢,采用10 000個樣本點對DTO和RBTO設計進行蒙特卡洛仿真。4種情況對應的結構體積比(目標函數(shù)),及蒙特卡羅可靠性指標,如表1所示。由表1可知:上述4種情況對應的結構體積比分別為31.90%、33.14%、39.20%和39.26%;對應的蒙特卡羅可靠性指標分別為-0.011 8、2.025 7、2.002 8和1.984 5??梢钥闯?3種RBTO設計均很好地吻合了目標可靠性指標βt=2,滿足可靠性設計要求,最大相對誤差不超過1.3%。與RBTO設計相比,DTO設計的可靠性水平(-0.011 8)相對較低。

        表1 不同情況下懸臂梁的優(yōu)化結果

        比較圖3~圖6可知,簡諧激勵下DTO設計與BRTO設計在拓撲構型上存在顯著差異。與DTO設計相比,RBTO設計出現(xiàn)了更多的傳力路徑,材料消耗更多,更能滿足可靠性設計要求。對比工況2~工況4可知,考慮載荷振幅不確定性的RBTO設計與只考慮頻率不確定的RBTO設計相比,具有更大的結構尺寸和結構體積比。此外,考慮頻率確定載荷振幅不確定的RBTO設計與考慮載荷振幅和頻率均不確定的RBTO設計的拓撲構型基本相同,但后者結構體積比更大一點。

        4.2 簡支梁

        簡支梁結構尺寸如圖7所示。設計域離散為320×40個邊長為5 mm、厚度為1 mm的正四邊形單元,底邊兩端點固支,上邊緣中心位置處受到豎直向下的隨機簡諧激勵作用。隨機簡諧激勵的振幅大小的平均值和標準差分別為1 000 N和100 N,頻率的平均值和標準差分別為100 Hz和10 Hz。以載荷作用點豎直方向的振幅平方為極限狀態(tài)函數(shù)來構造可靠性約束,設置約束限值為0.5 mm,過濾半徑設為15 mm。

        圖7 簡支梁結構示意圖(mm)Fig.7 Illustration of simply-supported beam (mm)

        為了研究不同目標可靠性指標對諧響應可靠性拓撲優(yōu)化結果的影響,考慮載荷振幅和頻率不確定性的諧響應可靠性拓撲優(yōu)化方法在βt=1、2、3時的拓撲優(yōu)化結果,如圖8所示。由圖8可知,3個不同目標可靠性指標對應的拓撲構型存在明顯差異,其對應的結構體積比分別是37.14%、42.46%和49.07%。從優(yōu)化結果來看,結構體積比隨著可靠性指標取值的增大而增大,這說明諧響應下結構可靠性水平的提高需要消耗更多的材料。不同目標可靠性指標下固支梁的體積比迭代歷史,如圖9所示。由圖9可知,3條迭代曲線都實現(xiàn)了穩(wěn)定快速的收斂,證明了所提方法的穩(wěn)定性。

        圖8 不同目標可靠性指標下簡支梁的RBTO結果Fig.8 RBTO results of simply-supported beam under different target reliability indexes

        圖9 不同目標可靠性指標下簡支梁的體積比的迭代歷史Fig.9 Iterative histories of volume fraction ratio of simply-supported beam with different target reliability indexes

        為了驗證本文諧響應可靠性拓撲優(yōu)化方法的有效性,采用10 000個樣本點對不同目標可靠性指標下的RBTO設計進行蒙特卡羅仿真,結果如表2所示。由表2可知,不同目標可靠性指標(βt=1、2、3)對應的蒙特卡羅可靠性指標分別為0.994 0、2.006 5、3.035 7,很好地吻合了目標可靠性指標,最大相對誤差不超過1.2%。通過上述結果分析可知,考慮載荷振幅和頻率不確定的諧響應可靠性拓撲優(yōu)化方法能夠有效設計出滿足指定可靠性水平的拓撲結構。

        表2 不同目標可靠性指標下簡支梁的優(yōu)化結果

        4.3 L型梁

        L型梁的結構尺寸如圖10所示。設計域離散為25 600個邊長為5 mm、厚度為1 mm的正四邊形單元,上邊緣完全固支,右上端點受到豎直向下的隨機簡諧激勵作用,振幅大小和頻率的平均值分別為1 000 N和80 Hz。以載荷作用點豎直方向的振幅平方為極限狀態(tài)函數(shù)構造可靠性約束,設置約束限值為1.7 mm,過濾半徑R設為17 mm,目標可靠性指標為βt=2。

        圖10 L型梁結構示意圖(mm)Fig.10 Illustration of L-shaped beam (mm)

        為了研究隨機簡諧激勵的振幅和頻率的變異系數(shù)(coefficient of variations,COV)對諧響應可靠性拓撲優(yōu)化結果的影響,考慮2個不同的變異系數(shù)(變異系數(shù)=標準差/平均值),即COV 取0.01和0.05。當COV為0.01和0.05時的RBTO結果,如圖11所示。其對應的結構體積比分別是32.63%和35.65%??梢钥闯鲎罱K結構體積比隨著隨機變量的變異系數(shù)的增大而增大,這說明概率模型參數(shù)的變化對優(yōu)化結果有影響。

        圖11 不同變異系數(shù)下L型梁的RBTO結果Fig.11 RBTO results of L-shaped beam under different coefficient of variations

        此外,采用10 000個樣本點對L型梁在不同變異系數(shù)下的RBTO結果進行蒙特卡羅仿真,結果如表3所示。不同變異系數(shù)(COV為0.01、0.05)對應的蒙特卡羅可靠性指標分別為2.006 5和1.970 3,與目標可靠性指標βt=2的最大相對誤差不超過1.5%,驗證了RBTO結果的準確性,進一步說明了考慮載荷振幅和頻率不確定的諧響應可靠性拓撲優(yōu)化方法的有效性,可設計滿足指定可靠性水平的動響應結構。

        表3 不同變異系數(shù)下L型梁的優(yōu)化結果

        5 結 論

        為了改善簡諧激勵作用下結構的使用安全性和可靠性,本文將考慮不確定性的可靠性分析理論整合到諧響應拓撲優(yōu)化中,提出了一種有效的考慮載荷振幅和頻率不確定性的諧響應可靠性拓撲優(yōu)化方法。3個數(shù)值算例及蒙特卡羅仿真驗證了所提方法的有效性,結果表明:

        (1) 考慮不確定簡諧激勵的RBTO設計與DTO設計相比,結構擁有更多的傳力路徑,體積分數(shù)更大,可靠性水平更高,更能滿足可靠性設計要求。

        (2) 考慮載荷振幅和頻率不確定性的諧響應可靠性拓撲優(yōu)化方法能夠有效設計出滿足指定可靠度水平的結構。本文算例中,優(yōu)化后結構可靠性指標與目標可靠性指標之間的最大相對誤差不超過1.5%。

        猜你喜歡
        振幅不確定性導數(shù)
        法律的兩種不確定性
        法律方法(2022年2期)2022-10-20 06:41:56
        解導數(shù)題的幾種構造妙招
        英鎊或繼續(xù)面臨不確定性風險
        中國外匯(2019年7期)2019-07-13 05:45:04
        關于導數(shù)解法
        十大漲跌幅、換手、振幅、資金流向
        十大漲跌幅、換手、振幅、資金流向
        十大漲跌幅、換手、振幅、資金流向
        具有不可測動態(tài)不確定性非線性系統(tǒng)的控制
        滬市十大振幅
        導數(shù)在圓錐曲線中的應用
        97se亚洲国产综合自在线| 精品国产一区二区三区香| 国产香蕉视频在线播放| 激情航班h版在线观看| 欧美三级免费网站| 美女黄网站永久免费观看网站| 国产亚洲91精品色在线| 中国人妻与老外黑人| 国产成人亚洲综合无码| 国产激情视频在线| av一区二区在线网站| 亚洲成a人片在线观看无码3d| 亚洲av有码在线天堂| 人人妻人人澡av| 日韩中文字幕素人水野一区 | 成人国产精品一区二区视频| 国产精品jizz观看| 国产一区二区在线观看视频免费| 91色老久久偷偷精品蜜臀懂色 | 开心五月激情五月天天五月五月天| 日韩av无码一区二区三区| 亚洲国产另类久久久精品黑人| 日韩最新在线不卡av| 精品人妻在线一区二区三区在线| 国产乱人偷精品人妻a片| 巨熟乳波霸若妻在线播放| 极品美女尤物嫩模啪啪| 国产熟人精品一区二区| 成人无码α片在线观看不卡| 精品人妻无码中文字幕在线| 中文字幕日韩精品中文字幕| 欧洲熟妇色xxxx欧美老妇性| 日本一区午夜艳熟免费| 亚洲精品天堂在线观看| 国产亚洲精品久久情侣| 日韩一卡2卡3卡4卡新区亚洲| 亚洲αⅴ无码乱码在线观看性色| 国产伦精品一区二区三区| 精品久久久久久久久午夜福利| 亚洲中文字幕第一页在线| 青青草免费在线手机视频|