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

基于PB-NSGA-Ⅲ算法的高速離心泵葉輪性能優化研究*

2023-12-20 12:39:10賴喜德陳小明劉升波
機電工程 2023年12期
關鍵詞:優化設計

張 翔,賴喜德*,陳小明,劉升波

(1.西華大學 能源與動力工程學院,四川 成都 610039;2.成都成發科能動力工程有限公司,四川 成都 610052)

0 引 言

隨著高速離心泵在航空航天、石油化工等行業應用范圍的逐步擴大[1],人們對高速泵的性能也提出了更高的要求;但國內高速泵的設計研發起步較晚[2],現有產品的水力性能達不到使用要求。

為提高高速泵的水力性能,國內外學者針對高速泵優化設計開展了大量的研究,并已取得了一定的成果。

袁軍婭等人[3]建立了高速泵的性能預測模型,推導出了高速泵效率的數學表達式,以效率為優化目標,對高速泵進行了優化設計,優化后高速泵的水力效率提高了9%。

相較于理論計算,數值模擬計算可以更加精確地預測性能特性。由此,數值模擬結合正交試驗在高速泵的優化設計中逐步增多。

賀青等人[4]以空化余量為目標,采用正交優化方法得出優化方案,對航空燃油離心泵進行了抗汽蝕優化設計,優化后其蒸汽質量分數降低了19.98%;但其空化余量下降不明顯。陳建華等人[5]以揚程和水力效率為目標,采用正交優化方法,對高速井泵進行了優化設計,優化后高速井泵的揚程上升了36 m;但井泵的水力效率并未得到提高。

正交試驗得出的優化方案并不是在變量范圍內對數值進行全面搜索的結果,且優化后,部分性能目標的提升也不明顯。

嚴俊峰等人[6]將葉輪幾何參數作為優化變量,僅以效率為目標,采用遺傳算法在變量的約束范圍內全面搜索最優參數組合,為降低優化難度而構建出優化變量間的映射關系,對低比轉速高速泵進行了優化設計;然而其在優化中需要反復修正映射的系數,這加大了優化設計的工作量。為克服繁瑣的系數修正過程,XU Z K等人[7]結合近似模型與多島遺傳算法,以水力效率和揚程為優化目標,對高速電磁泵進行了性能優化,經過104次迭代得出了最優參數組合,優化后電磁泵的水力效率上升了6.23%;但該方法迭代次數多,計算時間長。為減少計算時長,白永明等人[8]基于NSGA-Ⅱ算法對高速泵進行了多目標優化設計,消除了揚程駝峰;但優化結果中可供設計人員選擇的方案過于單一。

因此,結合高速泵特點,有效縮減優化變量的數量以提高效率,減少迭代次數,探索優化變量和優化目標間的聯系,使優化后的性能參數組合在較大范圍內均勻分布,是掌握高速泵多目標優化設計的關鍵。

近年來,研究人員采用傅里葉幅度靈敏度檢驗擴展法(extended Fourier amplitude sensitivity test, EFAST)確定影響高速泵性能的主要幾何參數[9],并將其作為優化變量,以有效縮減變量個數;改進的智能算法在葉輪機械多目標優化方面[10-12]應用廣泛并取得了不錯的效果,這均為高速泵優化設計提供了參考。

高速泵的性能參數對幾何參數比較敏感[13],使得高速泵優化難度加大;且針對多目標優化問題,提升算法預測能力也是多目標優化的一個難點。已有研究并未利用改進遺傳算法,并結合高速泵特點對泵進行有效地優化設計;且現有研究中應用的算法預測能力較差,算法結果可供設計人員選擇的方案較少,高速泵的綜合性能提升不明顯。

針對高速泵的特點,首先,筆者應用EFAST法有效篩選優化變量;然后,采用支持向量機高效構建近似模型,在此基礎上采用PB-NSGA-Ⅲ算法[14]進行尋優;最終,形成一套針對高速泵的高效優化方法,以有效提升優化設計的能力。

1 高速泵葉輪多目標優化設計流程

某工程擬采用IS-60-50-158型高速泵。該泵主要由吸入管、誘導輪、葉輪和蝸殼等過流部件組成。但該泵的空化性能、水力效率以及揚程不能滿足使用需求。

為綜合提高該高速泵水力性能,在其他過流部件基本不變的條件下,可以通過優化葉輪來達到優化目標。筆者以IS-60-50-158型高速泵葉輪的幾何模型為基礎,將其作為初始模型進行多目標優化設計。

其中,葉輪的優化流程如圖1所示。

圖1 高速泵葉輪的優化流程Fig.1 Optimization process of high-speed pump impeller

在圖1的優化流程中,空化余量是空化性能的重要指標之一,空化余量需按現有誘導輪揚程考慮。

筆者以水力效率η、空化余量NPSHr和揚程H為優化目標,將泵性能參數組合為{η,NPSHr,H}。

優化流程分為兩步:

1)建立高速泵葉輪幾何參數與性能參數間的近似模型。其中,結合高速泵的結構特點以及設計要求,確定幾何參數的取值范圍;為提高優化效率,減少優化變量個數,通過靈敏度分析篩選出影響高速泵性能的主要幾何參數,將其作為優化變量;利用拉丁超立方試驗在樣本空間中生成優化變量的60組幾何參數組合;在CFX軟件中完成計算流體動力學(computational fluid dynamics,CFD)數值計算,獲取60組幾何參數組合對應的性能參數組合,采用支持向量機(support vector machine,SVM)建立高速泵幾何參數組合與性能參數組合間的近似函數,完成近似模型的構建。

2)利用所建立的近似函數,以PB-NSGA-Ⅲ算法為尋優基礎,找出性能最優的參數組合,完成多目標優化設計。

2 基于數值模擬的高速泵性能預測計算方法

筆者簡述了初始的高速泵模型,通過數值計算獲取初始高速泵的性能參數組合,同時確定數值計算方案,為后續的數值計算提供參考。

2.1 初始的高速泵模型

筆者的研究對象是IS-60-50-158型高速泵。該泵為單級懸臂式結構,轉速n=8 000 r/min,設計工況Qd=50 m3/h。

初始高速泵幾何參數如表1所示。

表1 初始高速泵的幾何參數

2.2 數值模擬方法

根據IS-60-50-158型高速泵的幾何參數,筆者運用CFturbo軟件建立了高速泵幾何模型并導出流體域。

高速泵的流體域如圖2所示。

圖2 高速泵流體域Fig.2 High-speed pump fluid domain

圖2中,高速泵流體域包括葉輪進出口段、葉輪和蝸殼。

為減少數值計算誤差,筆者在UG/NX中修改了流體域的幾何模型,進口段的長度為葉輪進口直徑的4倍。筆者采用MESH軟件將高速泵的流體域劃分為非結構網格,交界面和葉片頭部網格需要做加密處理,網格質量達到0.2以上;進口采用靜壓,出口采用質量流量;采用CFX進行數值模擬并分析其內部流場。

為了解初始模型水力性能,筆者通過數值模擬獲取了初始模型的性能參數,如表2所示。

表2 初始模型的性能參數

2.3 高速泵的優化設計目標參數

結合項目要求可知,高速泵的水力性能不能滿足使用需求。為此,筆者提出的設計參數如表3所示。

表3 高速泵的設計參數

對比表2和表3可以看到:初始的高速泵模型不能滿足設計要求,需探討其優化方法。

3 幾何參數與性能參數間的近似模型

筆者以提高高速泵水力性能為目標,應用PB-NSGA-Ⅲ算法尋優時,需要對目標中的性能特性值進行預測;但很難準確預測性能特性值。

因此,筆者利用近似模型建立高速泵葉輪性能參數組合與幾何參數組合的一一對應關系,預測幾何參數對應的性能參數。近似模型的本質是根據離散幾何參數樣本和性能參數構建近似函數。

筆者建立近似模型包括優化變量的選取、拉丁超立方試驗,以及利用支持向量機建立的幾何參數與性能參數間的近似函數。

3.1 優化變量的選取

為提高優化效率,減少優化變量個數,筆者結合幾何參數的取值范圍,通過靈敏度分析,篩選出影響高速泵性能的主要幾何參數,并將其作為優化變量。

由于幾何參數的約束范圍較寬,而全局敏感性分析中的EFAST法在約束范圍較寬的情形中具有優越性,其能夠準確分析幾何參數的敏感性。因此,筆者采用EFAST法對幾何參數進行全局靈敏度分析。

首先,需要明確幾何參數的約束條件。根據設計參數中的比轉速,結合速度系數法[16]的統計,可選取葉輪進口直徑Dj、葉輪直徑D2、葉輪出口寬度b2的取值范圍。袁壽其[17]在大量試驗和理論推導的基礎上,給出了葉片進口角β1、葉片出口角β2的推薦值。

筆者適度放寬取值范圍,得到約束條件,如表4所示。

表4 幾何參數的約束條件

根據水力損失模型[18],計算水力損失,可得揚程和水力效率,可計算泵空化余量[19-21],公式如下:

H=HT-Δh

(1)

(2)

(3)

式中:H為揚程;h為水力效率;HT為理論揚程;Δh為水力損失之和(水力損失包括吸入式水力損失、葉輪進口水力損失、葉輪流道摩擦損失、進口液流變向水力損失、葉輪出水口水力損失、蝸殼流道摩擦損失、蝸殼內擴散損失等);l為葉片進口稍前的絕對速度;v0為葉片進口稍前的相對速度;w0為進口相對速度。

筆者的主要目標是提高高速泵的水力效率和空化性能,其次是讓揚程達到設計要求,因為初始高速泵的揚程未達到設計要求但相差不大。

為找到影響高速泵性能的主因,全局靈敏度分析主要針對水力效率與空化余量。筆者結合幾何參數約束條件,根據水力損失模型及泵空化余量計算模型,采用EFAST法進行了靈敏度分析,得到一階和全局敏感性指數,如圖3所示。

圖3 幾何參數的敏感性指數Fig.3 Sensitivity index of geometric parameters

圖3中,從幾何參數的敏感性指數中可以篩選出主要幾何參數,將其作為優化變量。

一階敏感性指數代表單個變量對性能參數的敏感性;全局敏感性指數表示變量間的耦合作用。

在葉輪的幾何參數中,對空化余量影響最大的是葉輪進口直徑Dj,對水力效率影響最大的是葉輪出口寬度b2,其余3個幾何參數中,對水力效率以及空化余量影響最大的是葉片出口角β2。

因此,筆者選擇進口直徑Dj、葉輪出口寬度b2以及葉片出口角β2作為優化變量,記作幾何參數組合為{Dj,b2,β2}。

3.2 優化變量均勻樣本的生成

為提高數值模擬效率,優化變量的幾何參數組合需要均勻分布,從而在優化變量的約束范圍內,使建立的近似模型能夠預測對應的性能參數值。

在拉丁超立方抽樣框架下,筆者采用連續局部枚舉方法和平移傳播算法在抽樣空間內生成均勻樣本,生成的60組幾何參數組合如圖4所示。

圖4 樣本點分布圖Fig.4 Sample point distribution map

3.3 幾何參數與性能參數間近似函數的建立

筆者應用PB-NSGA-Ⅲ算法對高速泵葉輪進行多目標優化設計,迭代中需要確定各性能參數與幾何參數組合{Dj,b2,β2}間的函數表達式。但對于高速泵的各性能參數,實際上很難使用數學表達式來描述。

在CFturbo中,筆者利用60組幾何參數組合生成幾何模型,采用數值計算獲取60組幾何模型對應的性能參數值組合,構建出幾何參數組合與性能參數值組合間的一一對應關系,即建立幾何參數組合與性能參數值組合間的近似函數,便于尋優時利用近似函數關系。

為高效構建近似函數,筆者采用林智仁基于支持向量機開發的LIBSVM工具箱[22-24],其在構建小樣本問題的近似函數中具有良好的擬合預測效果。

支持向量機的應用分為訓練集和測試集兩個部分。在60組樣本中,隨機抽取50個樣本作為訓練集,剩余10個作為測試集。支持向量機構造近似函數的過程較為簡單,在此不加贅述。

為使近似模型的精度更高,減小模型參數對模型的影響,同時抑制過擬合和欠擬,筆者應用LIBSVM工具箱中的RBF核函數,采用交叉驗證策略選取懲罰因子(c)和方差(g)的最佳參數組合,利用最佳參數組合構建出訓練模型。

因為近似函數的精度對尋優有較大的影響,所以筆者采用回歸分析對比CFD計算值與支持向量機預測值,以評估近似函精度,如圖5所示。

圖5 支持向量機測試集的回歸分析Fig.5 Regression analysis of support vector machine test set

在圖5中,測試集中的性能參數組合包含支持向量機所建立的近似函數的預測結果與CFD計算結果,回歸分析中的支持向量機預測值為支持向量機預測幾何參數得出的性能參數值;回歸分析中的CFD計算值為采用CFX軟件獲取的幾何參數對應的性能參數值,兩者形成對比,得到決定系數。

圖中的直線為45°對角線,其含義為支持向量機的預測值與CFD計算的值相等。決定系數R2越大,說明兩者吻合性越好。決定系數大于0.9,說明支持向量機建立的性能參數與幾何參數間的近似函數精度較高,預測結果有很高的精度。

4 高速泵葉輪的多目標優化

筆者將構造出的近似模型與改進的智能算法相結合。其能夠找出性能最優的幾何參數組合。但改進的智能算法在3個目標以上的優化問題中較難收斂。因此,在高速泵優化中,需要采用近似模型預測幾何參數對應的性能參數,應用PB-NSGA-Ⅲ算法找到最優性能的參數組合。

由于高速泵轉速高,其性能對幾何參數很敏感,子群的適應度范圍較寬,從圖5中支持向量機的預測值可以看到揚程范圍較寬。為此,筆者采取自適應策略,避免“早熟”現象。此外,為減少迭代次數,筆者采用PB-NSGA-Ⅲ算法,并和基于參考點的非支配遺傳算法(NSGA-Ⅲ算法)作對比。

PB-NSGA-Ⅲ與NSGA-Ⅲ算法的流程如圖6所示。

圖6 PB-NSGA-Ⅲ算法與NSGA-Ⅲ算法流程Fig.6 PB-NSGA-III algorithm and NSGA-III algorithm process

從圖6中可以看到:PB-NSGA-Ⅲ算法與NSGA-Ⅲ算法不同的是其采用了雙指標準則,在程序中第一個指標為個體與帕累托前沿的距離;第二個指標為個體到理想點的距離,其中理想點為每代子群中每一個目標最小值的組合,筆者取兩個指標中的較小值作為選擇操作的依據[14]。

5 優化結果及分析

結合高速泵的特點,筆者應用多學科優化方法以提高其水力性能。為驗證該方法的正確性,筆者利用基于數學和計算機的NSGA-Ⅲ算法和PB-NSGA-Ⅲ算法,找出性能最優的參數組合。參數組合包含幾何參數組合與性能參數組合,對比分析兩種算法結果,從更優秀的算法結果中選出優化后的參數組合,將其作為優化方案。

筆者結合數值模型并運用葉輪機械相關理論,分析優化方案的參數組合,確認葉輪幾何參數組合是否合理;并進一步對比優化前后葉輪的流場,分析優化前后流場內流線、靜壓及湍動能的變化,驗證優化效果。

5.1 優化結果分析

筆者應用NSGA-Ⅲ算法和PB-NSGA-Ⅲ算法在MATLAB中完成相應程序。兩種算法均設定參考點數量為200個,種群數量為200個,迭代次數為400。

經過程序迭代得到200組參數組合,200組性能參數組合{η,NPSHr,η}在空間的分布能夠反映多目標優化的效果,如圖7所示。

圖7 帕累托最優前沿Fig.7 Pareto optimal front

從圖7中可以看到:PB-NSGA-Ⅲ算法得到的性能參數組合集中在狹長窄帶上,分布更為均勻,優化結果中最高水力效率為87.24%,最低效率為82.62%,其中180組性能參數組合的水力效率大于83.4%;NSGA-Ⅲ算法得到的結果較為集中,優化結果中最高水力效率為85.24%,最低效率為82.04%,其中114組性能參數組合的水力效率大于83.4%。

針對高速泵的優化,PB-NSGA-Ⅲ算法的多目標優化能力更強。水力效率相對較高的解,空化余量也相對較高;但其揚程相對較低。

在PB-NSGA-Ⅲ算法的200組性能參數組合中,兼顧水力效率與空化余量,同時在揚程滿足設計要求的情況下,筆者篩選出水力效率≥85.5%、空化余量≤9.6 m、揚程≥220 m的性能參數組合,共有17組,然后選擇水力效率最高的一組作為優化方案。

優化前后幾何參數如表5所示。

表5 初始方案與優化方案幾何參數對比

由表5可知:不同數值可反映優化前后幾何模型的變化。

優化前后高速泵的水力性能如表6所示。

表6 初始方案與優化方案性能對比

在表6中,筆者結合數值模擬驗證了優化方案中的參數組合。其中,水力效率提高了2.22%,空化余量降低了3.20%,揚程提高了6.60%;在設計工況下,對比支持向量機的預測值和數值計算獲取的結果,水力效率偏差為0.96%,空化余量偏差為1.06%,揚程偏差為1.77%,支持向量機的預測值和數值計算獲取的結果偏差在2%以內。

以上結果再次證明近似模型有很高的精度。

5.2 優化前后外特性分析

優化方案中的性能參數組合是基于數學和計算機知識得出的設計工況點的水力性能,需運用葉輪機械理論進一步驗證優化方案中的幾何參數組合是否合理。數值計算可獲取優化前后高速泵在各運行工況下的水力效率以及揚程。

優化前后水力效率和揚程如圖8所示。

圖8 初始方案與優化方案外特性曲線對比Fig.8 Comparison of the external characteristic curves of the initial scheme and the optimization scheme

從圖8中可以看到:最高水力效率出現在所設計的工況附近。優化前后在0.7Qd~0.9Qd工況的水力效率有幾處相等,但優化后在1.0Qd、1.2Qd、1.4Qd工況下,水力效率分別提高了2.22%、5.73%、30.04%;在0.6Qd、0.8Qd、1.0Qd、1.2Qd、1.4Qd工況下,揚程分別提高了0.96%、1.97%、6.60%、16.09%、42.83%。

所以,優化后高速泵在大流量工況下水力效率更高,高效率區更寬;優化后高速泵在各個工況點的揚程均有提高。

5.3 高速泵的流場分析

優化前后0.1倍葉高處的流線圖如圖9所示。

圖9 優化前后葉輪中間截面的流線圖Fig.9 The streamline diagram of the middle section of the impeller before and after optimization

從圖9中可以看到:優化前后大多數區域的流線方向是穩定的,不同的是,優化前靠近葉片尾部出現了漩渦,其會導致能量損失,對高速泵的穩定性以及水力效率均有一定的影響;而優化后漩渦消失,說明優化后的葉輪內部流態更為合理,流動損失減小。

優化前后設計工況點葉片靜壓分布如圖10所示。

圖10 優化前后葉片靜壓圖Fig.10 Blade static pressure diagram before and after optimization

從圖10中可以看到:優化后,葉片背面的靜壓分布更加均勻。在葉片頭部,消除了較小的(-0.514 MPa~-0.256 MPa)的靜壓,葉片頭部低壓區獲得了一定的控制。這是因為優化后葉輪進口直徑相對較小,葉輪進口、輪轂以及葉片進口邊所形成的空間更小,葉輪進口到葉片進口區域的壓差增大,葉片頭部的靜壓提高,從而使葉片頭部的低壓區得到一定控制。

優化前后設計工況葉輪中間截面湍動能分布如圖11所示。

圖11 優化前后中間截面湍動能分布圖Fig.11 The turbulent kinetic energy distribution map of the middle section before and after optimization

從圖11中可以看到:湍動能是描述流場中湍流發展的一個重要物理量,湍動能值越小,則湍流發散的程度越小,葉片對流體的做功能力越強。從葉片頭部一直到葉片尾部的流道內,湍動能值較小;但在葉輪與蝸殼交界面附近的區域湍動能較大,該區域的湍流程度較大,湍流的脈動長度以及時間長度較大。

優化后,葉輪與蝸殼交界面附近湍流值較大的區域明顯減小,該區域的湍流耗散程度減小,說明優化后葉輪對流體做功更強。

6 結束語

筆者針對高速泵葉輪進行了多目標優化設計,優化后的高速泵在水力性能上有大幅提高,其可為高速泵的優化設計提供參考。

研究結論如下:

1)筆者基于近似模型和PB-NSGA-Ⅲ算法,對高速泵葉輪進行了多目標優化設計。優化方案中,近似模型的預測值和數值計算獲取的結果偏差在2%以內,決定系數R2大于0.9,表明近似模型的預測精度很高。優化后,設計工況點的水力效率提高了2.22%,揚程滿足設計要求,在大流量工況下水力效率提升明顯,也表明該優化方法是可靠的;

2)優化后在0.1倍葉高處,靠近葉片尾部的漩渦消失,使泵在運行過程中的穩定性得到了提高。葉片頭部的靜壓有所提高,泵的空化余量下降,空化性能在一定程度上有所提高;

3)相對于其他優化設計方法,運用PB-NSGA-Ⅲ算法結合近似模型能夠在變量取值范圍內更全面地搜索最優解,優化結果中的性能參數組合在狹長窄帶上分布均勻,且沒有出現大量偏離狹長窄帶的解。因此,采用雙指標的PB-NSGA-Ⅲ算法對于高速泵的優化設計是十分有利的。

筆者研究了高速泵葉輪的優化設計,但是該研究中不包括含誘導輪的高速泵。因此,在后續研究中,筆者將對含誘導輪的高速泵進行優化設計,以便更好地提高高速泵的水力性能。

猜你喜歡
優化設計
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
由“形”啟“數”優化運算——以2021年解析幾何高考題為例
何為設計的守護之道?
現代裝飾(2020年7期)2020-07-27 01:27:42
《豐收的喜悅展示設計》
流行色(2020年1期)2020-04-28 11:16:38
瞞天過海——仿生設計萌到家
藝術啟蒙(2018年7期)2018-08-23 09:14:18
設計秀
海峽姐妹(2017年7期)2017-07-31 19:08:17
有種設計叫而專
Coco薇(2017年5期)2017-06-05 08:53:16
主站蜘蛛池模板: 一级做a爰片久久毛片毛片| 无码'专区第一页| 国产区免费精品视频| 色婷婷啪啪| 国产国产人免费视频成18| 一本大道东京热无码av | 亚洲狼网站狼狼鲁亚洲下载| 国产黄色爱视频| 亚洲视屏在线观看| 国产乱视频网站| 动漫精品啪啪一区二区三区| 日韩欧美国产中文| 亚洲IV视频免费在线光看| 欧洲高清无码在线| 伊人网址在线| 一级香蕉视频在线观看| 亚洲精品少妇熟女| 欧美精品在线看| h视频在线播放| 欧美另类视频一区二区三区| 国产91丝袜| 午夜a视频| 久久国产热| www.国产福利| 欧美激情综合| A级毛片无码久久精品免费| 九九九九热精品视频| 亚欧乱色视频网站大全| 国产精品美人久久久久久AV| 波多野结衣第一页| 国产成人一区在线播放| 国产成人亚洲毛片| 91丨九色丨首页在线播放| 免费在线成人网| www.99在线观看| 国产浮力第一页永久地址| 国产最新无码专区在线| 亚洲第一区欧美国产综合| 成人91在线| 国产18在线播放| 欧美亚洲国产精品第一页| 无码一区二区三区视频在线播放| 第一区免费在线观看| 国产在线98福利播放视频免费| 国产成人高精品免费视频| 黄色在线网| 首页亚洲国产丝袜长腿综合| 亚洲一区二区精品无码久久久| 手机精品视频在线观看免费| 国产午夜一级毛片| 女人av社区男人的天堂| 大香伊人久久| 亚洲综合亚洲国产尤物| 91欧美亚洲国产五月天| 国产97公开成人免费视频| 野花国产精品入口| 国产精品99r8在线观看| 亚洲最大综合网| 青青草原国产精品啪啪视频| 亚洲欧美不卡视频| 国产精品成人第一区| 亚洲无码高清免费视频亚洲| 激情午夜婷婷| 亚洲91精品视频| 久久成人国产精品免费软件| 婷婷午夜影院| 国产精品私拍在线爆乳| 久久久久久尹人网香蕉 | 欧美性精品| 久久精品中文无码资源站| 91在线国内在线播放老师| 亚洲无线国产观看| 一级毛片在线播放免费观看 | 亚洲国产成人自拍| 白丝美女办公室高潮喷水视频| 91免费在线看| 精品人妻系列无码专区久久| 成人av专区精品无码国产| 成年人福利视频| 国产精品伦视频观看免费| 亚洲色无码专线精品观看| 99re在线免费视频|