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

基于興趣子域動態代理模型的船舶結構可靠性優化

2021-08-31 00:45:20羅文俊王德禹
中國艦船研究 2021年4期
關鍵詞:優化模型

羅文俊,王德禹

1 中遠海運特種運輸股份有限公司,廣東 廣州 510623

2 上海交通大學 海洋工程國家重點實驗室,上海 200240

3 高新船舶與深海開發裝備協同創新中心,上海 200240

0 引 言

基于可靠性的優化設計因在滿足船舶結構經濟性的同時能夠保證足夠的安全性而極具優勢。目前,可靠性優化設計已經在航空航天、汽車和船舶等領域得到廣泛應用[1-3]。

可靠性優化設計方法主要分為3 種:雙循環法、單循環法和解耦法。雙循環策略雖然精度高,但計算成本較大、效率較低,而后兩種方法能有效解決此問題。Liang 等[4]根據庫恩-塔克(KKT)最優條件代替可靠性分析過程,通過計算得到近似最小功能點,將概率約束轉化為確定性約束,提出了單循環法。劉勤等[5]采用單循環法將優化迭代計算與可靠壽命的迭代求解同步進行,并通過傳動箱的結構可靠性優化結果驗證了方法的高效性。單循環法雖然效率較高,但只適用于線性和弱非線性問題,對于強非線性問題,無法保證精度甚至有可能出現迭代不收斂的情況。而解耦策略在保證精度的前提下,不僅能提高求解效率,同時還具有良好的適用性。Du 和Chen[6]提出的序貫優化與可靠性評估(sequential optimization and reliability assessment, SORA)法通過將可靠性分析結果轉化為優化邊界條件的偏移向量,使確定性優化與可靠性分析被解耦。李海燕和井元偉[7]將SORA 法應用于多學科可靠性優化,降低了協同優化計算的復雜度。解耦方法較前2 種集成策略具有很大的優勢,而SORA 法是解耦方法中應用最為廣泛的一種方法。

船舶結構響應屬于多參數、強非線性問題,采用代理模型替代有限元模型能極大地提高可靠性優化效率,對于強非線性和高維度工程問題,Kriging 模型能實現良好的擬合。劉瞻等[8]通過結合重要抽樣法和人工蜂群算法,建立了參數優化的Kriging 模型,有效提高了可靠度計算精度和效率。張干鋒和王德禹[9]對常規Kriging 模型進行了分區間泛化和參數融合,其將改進后的Kriging模型用于船舶結構型線優化,在預測多維度響應時具有更高的精度。代理模型分靜態代理模型和動態代理模型2 種。靜態代理模型需要在設計空間中采集足夠的樣本點才能構造高擬合精度的代理模型,計算成本較大,而動態代理模型采用序貫抽樣方法,能以較少的樣本點實現對有限元模型的高度擬合。Echard 等[10]引入了一種學習函數,其以該學習函數作為篩選準則進行序貫抽樣來構造動態Kriging 模型,以較少的樣本點實現了對極限狀態邊界的高度擬合。高月華和王希誠[11]基于Kriging 代理模型提出了多點加點準則,基于樣本集信息和預測函數特征添加新樣本集,在尋優迭代之時自適應地提升代理模型的精度。

根據結構可靠性理論,最大可能失效點(MPP)周圍區域對可靠度的貢獻最大,而結構優化的最優解通常位于約束邊界附近,因此可靠性優化必然存在一個可能存在最優解的區域,即興趣子域。Zhao 等[12]通過計算影響可靠度主要區域的大小,提出興趣子域的概念,從而構造了在MPP 附近區域實現高度近似的代理模型。龍周等[13]引入少數類合成的過采樣算法(synthetic minority oversampling technique,SMOTE)對MPP 附近進行過采樣,構建對極限狀態邊界高度擬合的BP 神經網絡模型,極大地提高了計算效率和精度。

針對船舶結構因強非線性、多參數、多響應而導致計算成本較大、難以獲取大量樣本點的問題,本文擬提出可靠性優化策略。首先,提出并確定基于SORA 法中圓弧搜索法(arc search method,ASM)的興趣子域范圍,基于信息熵函數H制定自適應空間減縮規則,對設計空間進行不斷的縮減,進而構造基于興趣子域的自適應空間減縮序貫抽樣策略,以采用盡可能少的樣本點構造對興趣子域進行高度擬合的Kriging動態代理模型,減少計算成本;其次,提出基于SORA 法的概率約束可行性檢查方法,減少不必要的可靠性評估過程,提高可靠性分析效率;最后,將Kriging 模型和多島遺傳算法(multi-island genetic algorithm,MIGA)嵌入SORA 法中進行可靠性優化,以保證可靠性優化結果的準確性。

1 基于興趣子域的自適應空間減縮及序貫抽樣策略

1.1 圓弧搜索法

本文中的SORA 法采用ASM 法作為可靠性分析方法,該方法由Du 等[14]提出,用于將逆可靠性分析法中的優化過程轉化為數學迭代計算,從而簡化計算,提高可靠性分析效率。ASM 法的原理如下:逆可靠性分析法是在期望可靠度指標的等值面(||u||= βa)上尋找概率約束性能函數G(u)(其中,G>0 為可行性條件)的最小值,約束條件為可靠度指標期望值 βa。根據KKT 條件,可以將逆可靠性分析的優化過程推導為如下迭代公式:

ASM 法的迭代原理如圖1 所示。圖中:x1,x2為設計變量;Gk為概率約束性能函數第k次迭代取值。

圖1 圓弧搜索法[15]Fig. 1 Arc search method[15]

當式(1)迭代收斂至uk+1,即功能函數與期望可靠度指標面相切時,若G(uk+1)<G(uk),說明樣本點朝著功能函數G減小的方向前進,uk+1即為可靠性分析的逆最大可能失效點;當出現G(uk+1)≥G(uk)的情況時,進行第2 階段尋優,尋優空間為期望可靠度指標面上uk+1和uk之間的圓弧。采用以下數學優化模型進行尋優:

通過第2 階段的尋優迭代策略,可以保證新的迭代收斂點的功能函數值比上一個最優點的更小,實現對原迭代策略的補充完善。

本文充分利用ASM 法的尋優過程原理,提出基于興趣子域的自適應空間減縮及序貫抽樣策略,實現了對興趣子域的高度擬合;而對于非興趣子域,則采用空間減縮技術,無需對這些區域進行精確擬合,因而減少了可靠性優化的循環次數,可在保證擬合精度的同時提高效率。

1.2 興趣子域范圍的確定

ASM 法是在半徑為可靠度指標β 的圓(球面或超球面)上尋找功能函數的最小值,如圖2 所示,分別為當功能函數為凹函數和凸函數時的搜索過程(G>0 為可行域)。

圖2 圓弧搜索法尋優過程[15]Fig. 2 Optimization process of arc search method[15]

1.3 自適應空間減縮及序貫抽樣策略

1.3.1 自適應空間減縮規則

圖3 興趣子域示意圖Fig. 3 Schematic diagram of interest subdomain

構建對興趣子域進行高度擬合的動態代理模型,可同時提高確定性優化和可靠性分析過程的計算效率。本節將采用蒙特卡羅抽樣生成整個設計空間的樣本點集S,采用Lv 等[16]提出的AK-LS方法中基于信息熵的主動學習函數H實現序貫抽樣,在提高Kriging 模型全局擬合精度的同時,采用自適應的空間減縮技術不斷刪減樣本點集S中位于非興趣子域的樣本點,從而縮小設計空間的范圍直至興趣子域。學習函數H如下:

由于初始Kriging 模型的全局擬合精度較低,因此非興趣子域中樣本點的刪減條件應當遠離興趣子域。隨著Kriging 模型擬合精度的不斷提高,非興趣子域中樣本點的刪減條件自適應地不斷逼近興趣子域邊界,即在進行序貫抽樣更新動態代理模型的同時不斷減縮設計空間。最后,只保留興趣子域及其附近的部分樣本點,通過學習函數H構造代理模型來實現對興趣子域的局部高精度擬合。考慮到隨著樣本點的不斷增加,整個設計空間中樣本點的Hmax值逐漸減小,因此,本文在序貫抽樣的同時采用Hmax作為篩選因子。

本文采用的自適應空間減縮規則如下:

其中,篩選條件中下確界的設定參考了Cox和John[17]提出的用于優化的下置信界(lcb)函數,下置信界函數如下所示:

1.3.2 自適應空間減縮及序貫抽樣流程

基于興趣子域的自適應空間減縮及序貫抽樣策略是通過不斷減縮全局設計空間,直至保留興趣子域,同時通過具備高效提高全局精度特點的學習函數H建立Kriging 動態代理模型來擬合興趣子域。

Lv 等[16]通過大量的算例結果分析指出,當H函數對應的迭代停止準則為Hmax≤0.5 時,代理模型具有足夠的擬合精度。因此,為保證Kriging模型在興趣子域內具有足夠的擬合精度,本文擬采用Hmax≤0.5 作為自適應空間減縮和序貫抽樣的迭代停止準則。

基于興趣子域的自適應空間減縮及序貫抽樣策略流程如圖4 所示。

圖4 自適應空間減縮及序貫抽樣策略Fig. 4 Adaptive spatial reduction and sequence sampling strategy

1) 在設計空間中采用拉丁超立方抽樣或蒙特卡羅抽樣生成樣本集S,S集中的所有樣本點均無需調用有限元計算。

2) 生成訓練集T。采用最優拉丁超立方抽樣生成少量樣本點,作為訓練集T。

3) 調用功能函數或有限元軟件計算,得到訓練集響應值G。

4) 根據訓練集T,同時對約束條件和目標函數使用Matlab 自帶的DACE 工具箱建立多個Kriging 模型。

5) 分別用Kriging 模型預測S集中樣本點對應的響應值。

6) 采用學習函數H篩選最優樣本點ui,具體準則見1.3.1 節,同時根據空間縮減規則選擇樣本點集D。

7) 進行序貫抽樣迭代停止準則判斷。若滿足所有功能函數的迭代停止條件,則輸出代理模型,否則,將篩選的新樣本點ui加入訓練集,轉到步驟3)。若某個功能函數對應的代理模型已達到精度要求,則無需繼續對該代理模型進行最優樣本點篩選,同時從樣本集S中剔除樣本集D中的樣本點,從而更新S集,轉到步驟5),循環直到所有動態代理模型均滿足序貫抽樣迭代停止準則。

2 基于興趣子域的可靠性優化方法

2.1 基于SORA 法的概率約束可行性檢查方法

SORA 法作為最理想的可靠性優化解耦方法,在工程實際中得到了廣泛應用。不過該方法的提出者Du 和Chen[6]指出,該方法值得改進之處在可靠性評估尋找逆最大可能失效點的過程中,并非所有的概率約束都起作用,因此研究識別概率約束可行性檢查方法,以此避免不必要的可靠性分析過程可以極大地減少計算量,提高計算效率。如圖5所示,圖中:G1,G2,G3為3 個功能函數,x1,x2為2 個設計變量,由于最優點H2的位置遠離第2 個概率約束的極限狀態邊界,在對該最優點進行可靠性分析時,第2 個概率約束是不起作用的,所以無需對第2 個概率約束進行可靠性分析。

圖5 概率約束條件偏移示意圖Fig. 5 Schematic diagram of probability constraint deviation

本節采用ASM 法對SORA 法中的可靠性進行評估,在確定性優化計算得出最優解樣本點之后,以最優解樣本點為中心進行蒙特卡羅抽樣,生成若干個隨機樣本點,然后分別計算各個概率約束條件的失效概率以及整體失效概率。若某個概率約束條件的失效概率高于設計值,則采用ASM 法計算逆最大可能失效點,從而計算偏移向量,在下一循環的確定性優化中按該偏移向量移動該概率約束條件;若某個概率約束條件的失效概率低于設計值,則直接設置偏移向量為0 的向量,即無需移動該概率約束條件。

本節提出的概率約束可行性檢查方法操作簡單、易于實現,能極大地減少不必要的可靠性評估過程,提高了計算效率。

2.2 基于興趣子域的可靠性優化流程

MIGA 算法具備遺傳算法的優點,即適用于全局的、離散的、非光滑的優化問題,同時,還克服了傳統遺傳算法易早熟、局部優化能力弱且計算耗時的缺點。本節將具有全局尋優性能的MIGA算法和基于興趣子域的Kriging 模型嵌入SORA解耦方法中,通過SORA 法,將可靠性分析解耦成了確定性優化和可靠性分析的優化過程。本節提出的基于興趣子域的可靠性優化設計流程如圖6所示,歸納如下:

圖6 基于興趣子域的可靠性優化流程Fig. 6 Reliability optimization process based on interest subdomain

1) 在確定性空間中,采用自適應空間減縮及序貫抽樣策略構造基于興趣子域的Kriging 模型,具體過程詳見1.3.2 節。

2) 設置所有概率約束條件的偏移向量初值為0 的向量。

3) 在確定性優化部分,根據偏移向量V調整約束條件位置,采用MIGA 算法作為優化算法,通過每次調用Kriging 模型計算各個功能函數的響應值,并判斷各約束響應是否滿足約束條件,然后輸出符合約束條件的最優樣本點ui。

4) 基于最優樣本點ui,根據2.1 節中的概率約束可行性檢查方法,篩選出需要進行可靠性分析的概率約束條件,然后,采用ASM 法調用Kriging模型尋找逆最大可能失效點,計算偏移向量,以作為下一次循環過程確定性優化中對應約束條件的偏移向量。

5) 根據2.1 節中計算的整體失效概率Pf判斷是否停止迭代,若得到滿足約束條件要求的最優可行解,可靠性優化結束;否則,轉到步驟3)進行下一次循環。

2.3 數學算例效果校驗

本文采用文獻[18]中具有非線性概率約束的數學算例,來驗證本文所提基于興趣子域的自適應空間減縮及序貫抽樣策略的有效性。該算例的數學模型如下:

式中:ux表示由2 個設計變量x1,x2組成的樣本點; Φ為標準正態分布的分布函數;N表示正態分布;u0為初始樣本點; μxi為均值。

根據3 個概率約束函數Gi(X),在平面空間中繪制如圖7 所示曲線。

圖7 概率約束函數曲線Fig. 7 Curves of probability constraint function

根據1.3.2 節的自適應空間減縮及序貫抽樣流程,如圖8(a)所示,橫、縱坐標分別表示2 個設計變量,采用拉丁超立方抽樣方法,在設計空間中生成10 000 個樣本點作為樣本集S,生成10 個樣本點作為訓練集T。圖中,小圓點代表樣本集S,*點代表訓練集T。

圖8(a)~圖8(f)所示為該數學算例的自適應空間縮減過程圖。由圖可知,隨著序貫抽樣的進行,代理模型精度不斷提高,同時,設計空間不斷縮小至興趣子域。其中,至圖8(f)所示在第15 次空間縮減后,共剩余3 265 個樣本點,分布區域主要集中在興趣子域中,訓練集38 個樣本點,主要分布在約束邊界上。

圖8 自適應空間縮減過程圖Fig. 8 Flowchart of adaptive space reduction

由表1 可知,采用本文提出的基于興趣子域的可靠性優化方法得到了該數學算例的全局最優解,該最優解與理論解的相對誤差為0.066 8%,證明本文方法具有足夠的計算精度。且本文方法相比其他方法極大地減少了調用功能函數的次數,其中和文獻[18] 中最優的HCC 方法相比減少了40.6%的調用次數,證明本文可靠性優化方法高效。

表1 數學算例可靠度計算結果Table 1 Reliability calculation results of the math example

3 船舶艙段的可靠性優化

3.1 船舶艙段模型介紹

本文將采用基于興趣子域的可靠性優化設計方法,對一艘618TEU 型多用途船的貨艙艙段進行可靠性優化。該船主要裝載礦石等散裝貨物和集裝箱,其主尺度等參數如表2 所示。

表2 618TEU 型多用途船主尺度及航速Table 2 Main dimensions and speed of the 618TEU multipurpose vessel

本文優化的區域為該船的中間貨艙。由于模型結構、載荷和邊界條件都是沿中縱剖面對稱分布,為節省計算成本,船寬方向僅選取艙段的1/2 作為計算模型。同時,考慮到減少邊界條件的影響,船長方向選取1/2+1+1/2 艙段作為計算模型。艙段模型的垂向范圍則選取為船體型深。該艙段有限元模型如圖9 所示。

圖9 艙段有限元模型Fig. 9 Finite element model of the cabin

艙段模型的坐標規定為:船長方向為X軸,正方向由船艉指向船艏;船寬方向為Y軸,正方向由右舷指向左舷;型深方向為Z軸,正方向由基線指向甲板。

該艙段首、尾端面載荷對稱分布,可在兩端面建立多點約束(multi-point constraints,MPC),將MPC 與兩端面上節點耦合,從而施加約束和載荷。艙段模型邊界約束如表3 所示,有限元模型上邊界約束情況如圖10 所示。表3 中:Ux表示約束X方向的平移;Uy表示約束Y方向的平移;Uz表示約束Z方向的平移;URy表示約束繞Y軸的轉動;URz表示約束繞Z軸的轉動。

圖10 艙段模型約束情況Fig. 10 Constraints of the cabin model

表3 邊界條件施加表Table 3 Introduction of boundary conditions

貨艙艙段的載荷布置情況按中國船級社的《散貨船結構強度直接計算分析指南》[19]確定。有限元模型中施加的載荷主要分為3 個部分:端面彎矩、外部水壓力和艙室內貨物壓力,如圖11 所示。

圖11 艙段模型的載荷分布情況Fig. 11 Load distribution of the cabin model

3.2 船舶艙段可靠性優化數學模型

圖12 艙段橫剖面圖及設計變量分布情況[13]Fig. 12 Cross section of the cabin and distribution of design variables[13]

次確定性優化迭代后得出的最優設計方案所對應的板厚值;而標準差 σti則設為制造加工的許用偏差,本文取為對應板材厚度均值的2%,即ti~N(μti,(0.02μti)2)。

各設計變量的編號以及對應的名稱和初始值如表4 所示。

表4 設計變量參數表Table 4 Parameter list of design variables

式中:Gj(ti)為 各個概率約束條件的功能函數;gj(ti)為按照規范分類的板和梁的應力響應。gj(ti)對應的計算區域及許用應力要求如表5 所示。通過對表5 中6 個結構分類區域分別建立分組,輸出每個分組中單元絕對值最大的von Mises應力值作為表5 中6 個變量的應力響應值。

表5 g i(x)對應的計算區域及許用應力要求Table 5 Analysis domain and allowable stress for gi(x)

3.3 船舶艙段可靠性優化計算結果

本文通過Isight 優化平臺,集成Matlab,Patran和Nastran 軟件,采用2.2 節所述的基于興趣子域的可靠性優化流程,對此艙段結構進行可靠性優化設計。即采用Isight 自帶的MIGA 優化算法,采用Matlab 軟件實現SORA 算法,然后采用Patran軟件進行參數化建模及后處理,并調用Nastran 軟件進行模型分析計算。具體參數設置和流程如下:

1) 采用最優拉丁超立方抽樣技術生成180 個初始樣本點,對6 個概率約束條件和目標函數分別構造初代Kriging 模型,進而采用基于興趣子域的自適應空間減縮及序貫抽樣策略構造7 個Kriging動態代理模型。

2) 采用SORA 算法結合MIGA 算法進行可靠性評估和優化,具體流程詳見2.2 節。并將不考慮概率約束的優化結果與可靠性優化結果進行對比,結果如表6 和表7 所示。

表6 可靠性優化設計與無概率約束優化結果對比Table 6 Comparison between reliability-based design optimization and probability-free constrained optimization

由表7 可知,對船舶艙段進行可靠性優化后,可靠度達99.68%,滿足概率約束條件;艙段總質量相比初始方案下降了6.32%,優化效果顯著;艙段總質量相比確定性優化增加了0.728%,符合基于可靠性的優化設計原理,即考慮到設計變量的不確定性,通過偏移向量將邊界條件朝著可靠域偏移,從而犧牲部分經濟性來達到可靠度指標。通過對比文獻[13]的可靠性優化結果可知,本文與文獻[13]相比少調用了94 次有限元計算次數,且艙段總質量和文獻[13]的結果相比減少了0.883%,證明了本文所提可靠性優化策略的高效性和適用性。

表7 可靠性優化設計與無概率約束優化結果對比Table 7 Comparison between reliability-based design optimization and probability-free constrained optimization

采用蒙特卡羅法對上述可靠性優化結果進行驗證。在設計空間均勻生成106個隨機樣本點,采用蒙特卡羅仿真(Monte Carlo simulation,MCS)得到可靠度Pr= 99.67%,由此可知,本文的可靠性優化方案滿足可靠度要求。為檢驗本文中基于興趣子域的自適應空間減縮及序貫抽樣策略所構建的動態代理模型的預測精度,將最終方案進行了有限元計算驗證,表8 所示為概率約束條件和優化目標對比結果。

艙段質量作為優化目標,可靠性優化設計和有限元分析兩種方式的計算結果分別為456.6 t 和452.536 t,相 對 誤 差 為0.898 1,因 此 綜 合 表8 可知,本文構建的替代艙段模型的Kriging 模型對各個約束功能函數和優化目標的預測相對誤差均在1%以內,驗證了本文所提可靠性優化方法的精度,以及其在船舶結構領域的適用性。

表8 可靠性優化結果與有限元模型結果對比Table 8 Comparison between result of reliability-based design optimization and FEM

4 結 論

本文針對船舶結構可靠性優化,采用序貫抽樣策略構建了高精度、高效率的動態代理模型,用以達到使用盡可能少的樣本點構建具有足夠擬合精度的代理模型的目的,從而減少調用耗時的有限元模型計算的次數,節約計算成本。得到的主要結論如下:

1) 通過研究本文所采用SORA 算法中ASM法的尋優原理,以及確定性優化的最優解可能存在的區域,確定了基于可靠性優化的興趣子域范圍,并提出了基于興趣子域的自適應空間減縮序貫抽樣策略,通過制定自適應空間減縮規則,在采用主動學習函數進行序貫抽樣的同時,逐步對設計空間進行了縮減,從而能采用盡可能少的樣本點構造對興趣子域進行局部高度擬合的動態代理模型。

2) 針對SORA 法在可靠性評估尋找逆最大可能失效點的過程中并不是所有的概率約束都起作用的問題,提出了概率約束可行性檢查方法,該方法操作簡單、易于實現,能極大地減少不必要的可靠性評估過程,節省了計算成本。

3) 將基于興趣子域構建的動態Kriging 代理模型與MIGA 算法相結合嵌入SORA 算法中,構建了基于興趣子域的可靠性優化方法,通過采用數學算例,驗證了該方法的計算效率和精度。

4) 將基于興趣子域的可靠性優化方法應用于船舶艙段結構可靠性優化設計中,結果顯示艙段總質量相比初始方案下降了6.32%,相比確定性優化增加了0.883%,可靠度達到了99.68%。該方案在保證結構安全性的同時提升了經濟性,證明了本文所提可靠性優化策略的高效性和適用性。

猜你喜歡
優化模型
一半模型
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
由“形”啟“數”優化運算——以2021年解析幾何高考題為例
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 国产情精品嫩草影院88av| 精品91在线| 精品色综合| 国产网友愉拍精品视频| 不卡视频国产| 国产婬乱a一级毛片多女| 国产精品分类视频分类一区| 特级aaaaaaaaa毛片免费视频 | 国产精品免费p区| 久久久久青草大香线综合精品| 欧美成人怡春院在线激情| 91探花在线观看国产最新| 国产va在线| 欧美午夜网| 无码日韩人妻精品久久蜜桃| 尤物午夜福利视频| 99久久国产综合精品2020| 国产精品伦视频观看免费| 在线日韩一区二区| 毛片免费网址| www.日韩三级| 亚洲欧美成人综合| 欧美日韩北条麻妃一区二区| 国产精品刺激对白在线| 欧美国产菊爆免费观看| 亚洲精品桃花岛av在线| 免费无码一区二区| 性视频久久| 日韩第九页| 少妇极品熟妇人妻专区视频| 国内精品小视频在线| 精品视频在线观看你懂的一区| 波多野结衣亚洲一区| 热这里只有精品国产热门精品| 国产在线专区| 国产亚洲精久久久久久无码AV| 国产精品成人免费视频99| 国产综合色在线视频播放线视| 91成人在线免费观看| 人人看人人鲁狠狠高清| 国产一级特黄aa级特黄裸毛片| 中文字幕在线日本| 全裸无码专区| 一本大道香蕉中文日本不卡高清二区 | 国产精鲁鲁网在线视频| 国产在线专区| 夜精品a一区二区三区| 白丝美女办公室高潮喷水视频| 日韩午夜福利在线观看| 日韩国产一区二区三区无码| 欧美视频在线观看第一页| 手机精品视频在线观看免费| 2021天堂在线亚洲精品专区| 国产爽歪歪免费视频在线观看| 美女被操黄色视频网站| 国产亚洲精品自在久久不卡| 十八禁美女裸体网站| 国产黄色爱视频| 免费国产无遮挡又黄又爽| 午夜性爽视频男人的天堂| 久久这里只有精品66| 欧美高清国产| AV不卡国产在线观看| 99国产在线视频| 亚洲综合极品香蕉久久网| 亚洲无线观看| AV不卡无码免费一区二区三区| 久久青草免费91线频观看不卡| 亚洲最猛黑人xxxx黑人猛交| 40岁成熟女人牲交片免费| 色综合久久无码网| 国产精品粉嫩| 国产丝袜无码精品| 91久久大香线蕉| av色爱 天堂网| 久久精品国产免费观看频道| 伊人激情综合网| 又粗又大又爽又紧免费视频| 国产区精品高清在线观看| 国模沟沟一区二区三区| www.youjizz.com久久| 中文字幕亚洲综久久2021|