• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    基于代理模型方法的串列泵優(yōu)化設(shè)計(jì)

    2016-04-25 06:21:36趙宇王國玉黃彪王復(fù)峰
    關(guān)鍵詞:軸流泵

    趙宇, 王國玉, 黃彪, 王復(fù)峰

    (北京理工大學(xué) 流體機(jī)械工程研究所,北京 100081)

    ?

    基于代理模型方法的串列泵優(yōu)化設(shè)計(jì)

    趙宇, 王國玉, 黃彪, 王復(fù)峰

    (北京理工大學(xué) 流體機(jī)械工程研究所,北京 100081)

    摘要:針對葉輪機(jī)械傳統(tǒng)優(yōu)化設(shè)計(jì)周期長、優(yōu)化效率低等不足,提出了一種基于代理模型的優(yōu)化設(shè)計(jì)方法,并應(yīng)用于串列泵的優(yōu)化設(shè)計(jì),分析了關(guān)鍵設(shè)計(jì)參數(shù)對串列泵性能的影響。選擇首、次級葉輪葉片安放角和相位角作為優(yōu)化參數(shù),選擇效率和最小軸向截面平均壓力系數(shù)作為優(yōu)化設(shè)計(jì)目標(biāo)函數(shù),用來表示串列泵的外特性和空化特性。采用不同的方法建立代理模型,并采用敏感度分析方法和Pareto front方法進(jìn)行參數(shù)影響分析和最優(yōu)點(diǎn)的選取。采用數(shù)值計(jì)算的方法對優(yōu)化后的串列泵外特性和空化特性進(jìn)行驗(yàn)證,得到結(jié)論如下:提出的方法可以較好地應(yīng)用于串列泵的優(yōu)化設(shè)計(jì),優(yōu)化結(jié)果表明串列泵設(shè)計(jì)流量附近的效率和空化性能均得到提升。首、次級葉輪相位角對串列泵性能影響較??;首、次級葉輪葉片安放角對串列泵效率的敏感度之比和設(shè)計(jì)的載荷之比相同,首級葉輪葉片安放角是串列泵的空化特性的主要影響因素。

    關(guān)鍵詞:軸流泵;串列葉柵;系統(tǒng)優(yōu)化設(shè)計(jì);代理模型方法;空化特性

    串列葉柵是指多排葉柵直接相連的布置方式,具有靈活的氣動特性,目前在航空發(fā)動機(jī)[1-3]和水下螺旋槳推進(jìn)中[4-5]已經(jīng)得到應(yīng)用。串列泵是一種采用串列葉柵布置的泵,和普通的雙級泵相比,由于減少一級導(dǎo)葉而顯著減小泵級的軸向尺寸和質(zhì)量,從而使得串列泵具有結(jié)構(gòu)緊湊和功率密度高的特點(diǎn);又由于兩級動葉輪直接串聯(lián),前后級葉輪流動的相互作用會導(dǎo)致前后級葉輪流場變化,因此串列泵具有特殊的能量特性[6-7]。對于傳統(tǒng)泵來說,效率和空化性能不可兼得[8],而串列泵通過調(diào)整前后葉柵的載荷分配可以兼顧效率和空化性能[9-10]。盡管串列泵具有上述諸多優(yōu)點(diǎn),但與傳統(tǒng)泵相比,其優(yōu)化設(shè)計(jì)過程需要考慮更多的參數(shù),因此更為復(fù)雜。王立祥等[10-11]采用實(shí)驗(yàn)方法對串列泵進(jìn)行優(yōu)化設(shè)計(jì),分步順序處理不同設(shè)計(jì)參數(shù)對整體性能的影響,該優(yōu)化設(shè)計(jì)過程具有周期長,并且無法得到各設(shè)計(jì)參數(shù)的影響權(quán)重等缺點(diǎn)。

    系統(tǒng)優(yōu)化方法和敏感度分析方法的發(fā)展推動了準(zhǔn)確高效的現(xiàn)代設(shè)計(jì)方法的進(jìn)步,Shyy等[13]提出了一種基于代理模型優(yōu)化設(shè)計(jì)方法,并較好地應(yīng)用于解決航空動力學(xué)以及火箭燃料推進(jìn)問題。和傳統(tǒng)的優(yōu)化方法相比,基于代理模型的系統(tǒng)優(yōu)化設(shè)計(jì)方法具有以下優(yōu)點(diǎn)[14]:可以同時(shí)處理由多種不同方法獲得的數(shù)據(jù),例如實(shí)驗(yàn)結(jié)果或仿真結(jié)果;可以同時(shí)考慮多個(gè)設(shè)計(jì)變量對結(jié)果的影響,不用對每個(gè)變量進(jìn)行單獨(dú)分析;可以處理存在多個(gè)最優(yōu)點(diǎn)的優(yōu)化問題;提供多種準(zhǔn)則的優(yōu)化過程而且可以有效過濾實(shí)驗(yàn)或仿真結(jié)果中固有的噪聲?;诖砟P偷膬?yōu)化設(shè)計(jì)方法在航空動力學(xué)中已經(jīng)得到廣泛的應(yīng)用。John S等[15]應(yīng)用系統(tǒng)優(yōu)化方法對離心泵葉片的優(yōu)化設(shè)計(jì),結(jié)果表明系統(tǒng)優(yōu)化方法可以快速收斂并且優(yōu)化結(jié)果良好。 Zhang等[16]采用多目標(biāo)優(yōu)化方法設(shè)計(jì)軸流式多相流泵,結(jié)果表明泵級壓增上升10%的同時(shí)效率可以提升約3%。 J. Fan等[17]采用基于代理模型的方法對射流泵進(jìn)行優(yōu)化設(shè)計(jì),結(jié)果表明射流泵的效率增加4%的同時(shí)輸入功率可以減小20%以上。此外,J H Kim[18],Mustafa G?lcü等[19]也分別采用系統(tǒng)優(yōu)化設(shè)計(jì)方法對水力機(jī)械進(jìn)行優(yōu)化設(shè)計(jì)。

    本文采用了基于代理模型的系統(tǒng)優(yōu)化設(shè)計(jì)方法,對一串列泵進(jìn)行優(yōu)化設(shè)計(jì),采用計(jì)算流體動力學(xué)方法(computational fluid dynamics, CFD)來獲取數(shù)據(jù),建立了代理模型,對設(shè)計(jì)參數(shù)進(jìn)行優(yōu)化,并進(jìn)一步分析了首、次級葉輪進(jìn)口安放角度和相位角等關(guān)鍵設(shè)計(jì)參數(shù)對串列泵性能參數(shù)的影響。

    1系統(tǒng)優(yōu)化設(shè)計(jì)方法

    系統(tǒng)優(yōu)化設(shè)計(jì)方法流程主要包含實(shí)驗(yàn)設(shè)計(jì)、數(shù)值實(shí)驗(yàn)、代理模型的建立和驗(yàn)證等過程。

    1.1實(shí)驗(yàn)點(diǎn)設(shè)計(jì)

    在實(shí)驗(yàn)點(diǎn)設(shè)計(jì)中,通過設(shè)計(jì)方法選定實(shí)驗(yàn)點(diǎn)。然后采用數(shù)值計(jì)算的方法獲取每個(gè)實(shí)驗(yàn)點(diǎn)對應(yīng)的數(shù)據(jù),進(jìn)而建立代理模型。

    在建立代理模型之前,目標(biāo)函數(shù)和設(shè)計(jì)變量之間的關(guān)系是未知的,一般采用較為簡單的隨機(jī)分布方法。如果考慮到數(shù)值實(shí)驗(yàn)的計(jì)算成本,實(shí)驗(yàn)點(diǎn)的個(gè)數(shù)需要被控制,可以采用更復(fù)雜的實(shí)驗(yàn)點(diǎn)分布方法,例如拉丁超立方分布方法(Latin hypercube sampling, LHS),該方法是一種受約束的實(shí)驗(yàn)點(diǎn)設(shè)計(jì)方法,首先將設(shè)計(jì)空間的每個(gè)坐標(biāo)軸區(qū)間等分成互不重疊的子區(qū)間,然后在每個(gè)因素的每個(gè)子區(qū)間內(nèi)只進(jìn)行一次實(shí)驗(yàn)點(diǎn)分布[25]。在本次研究中,在基于LHS方法的同時(shí)采用面心分布方法(face-centered composite design, FCCD)[25]設(shè)計(jì)實(shí)驗(yàn)點(diǎn)。FCCD方法在設(shè)計(jì)空間的面心和頂點(diǎn)位置布置實(shí)驗(yàn)點(diǎn),用來保證設(shè)計(jì)空間邊界附近實(shí)驗(yàn)點(diǎn)的數(shù)目。

    1.2代理模型

    對于不同的問題,應(yīng)當(dāng)采用不同的方法建立代理模型并進(jìn)行比較,從而得到最好的方法。本文分別采用多項(xiàng)式法(polynomial response surface, PRS)[26],克里金法(Kriging, KRG)[27]和徑向神經(jīng)網(wǎng)絡(luò)法(radial-based neural network, RBNN)[13]建立代理模型,并比較不同代理模型預(yù)測結(jié)果的差異。

    1.2.1多項(xiàng)式響應(yīng)面法

    在多項(xiàng)式方法中,假設(shè)目標(biāo)函數(shù)是若干個(gè)設(shè)計(jì)變量多項(xiàng)式的線性組合:

    (1)

    式中:βi是系數(shù)矩陣。對于大多數(shù)工程問題,二階多項(xiàng)式擬合精度均可達(dá)到工程需求,同時(shí)考慮計(jì)算成本,本文中采用二階多項(xiàng)式方法。

    1.2.2克里金法

    克里金法是全局模型和局部偏差的組合,該方法假設(shè)局部偏差只與考察對象位置之間的距離有關(guān):

    (2)

    (3)

    式中:Ndv表示設(shè)計(jì)變量的維數(shù),σ是實(shí)驗(yàn)點(diǎn)響應(yīng)值的標(biāo)準(zhǔn)差,φk是衡量k方向上相關(guān)性的參數(shù)。

    1.2.3徑向神經(jīng)網(wǎng)絡(luò)法

    在徑向基函數(shù)神經(jīng)網(wǎng)絡(luò)法中,假設(shè)網(wǎng)絡(luò)為單輸出,它是一種三層(輸入層,隱含層,輸出層)前向結(jié)構(gòu)網(wǎng)絡(luò),隱含層有一組單元節(jié)點(diǎn),每一個(gè)節(jié)點(diǎn)均有一個(gè)稱作中心的參數(shù)矢量,每一個(gè)節(jié)點(diǎn)計(jì)算網(wǎng)絡(luò)輸入矢量與中心參數(shù)矢量間的歐氏距離,而后通過一非線性函數(shù)映射,將結(jié)果傳遞輸出層;而輸出層對隱含層各節(jié)點(diǎn)的輸出函數(shù)值作線性組合。

    目標(biāo)函數(shù)可以表示為

    (4)

    (5)

    式中:‖si-x‖是權(quán)重函數(shù)和輸入量之間的矢量距離,α是歸一化參數(shù)。

    1.3模型驗(yàn)證

    為了選取合適的模型來進(jìn)行問題分析,需要建立統(tǒng)一的誤差分析方法來評價(jià)模型的精度。常見的模型驗(yàn)證方法有分割驗(yàn)證,交叉驗(yàn)證,自助驗(yàn)證法等,本文采用的是交叉驗(yàn)證[28]。

    交叉驗(yàn)證建立在分割驗(yàn)證的基礎(chǔ)上將實(shí)驗(yàn)數(shù)據(jù)分成k個(gè)子集(k階),然后選定一組實(shí)驗(yàn)數(shù)據(jù),用其他的所有實(shí)驗(yàn)點(diǎn)來建立代理模型,利用該組實(shí)驗(yàn)數(shù)據(jù)進(jìn)行驗(yàn)證,然后選擇另一組實(shí)驗(yàn),重復(fù)建模-驗(yàn)證過程。該方法需要建立k次代理模型,然后對k次模型驗(yàn)證的結(jié)果進(jìn)行平均,從而得到交叉驗(yàn)證的最終結(jié)果。交叉驗(yàn)證方法的優(yōu)點(diǎn)是幾乎能提供無偏的誤差估計(jì),而且由于每個(gè)點(diǎn)都被使用一次來進(jìn)行驗(yàn)證,所以交叉驗(yàn)證能大大減小預(yù)測誤差的方差。本文采用一階交叉驗(yàn)證方法,其衡量參數(shù)PRESS(predictionerrorsumofsquares):

    (6)

    1.4敏感度分析

    采用代理模型獲取自變量和目標(biāo)函數(shù)的關(guān)系之后,可以進(jìn)行敏感度分析,敏感度分析可以直觀體現(xiàn)自變量和目標(biāo)函數(shù)之間的關(guān)系,在各工程領(lǐng)域獲得廣泛的應(yīng)用。本文采用的敏感度分析方法是Sobol方法[29]。

    (7)

    各項(xiàng)的平方在標(biāo)準(zhǔn)空間的積分為對應(yīng)的變化量:

    (8)

    定義整體變化量如下:

    (9)

    從而可以定義對于變量xi的局部敏感度如下:

    (10)

    全局敏感度如下:

    (11)

    2基于代理模型的串列泵優(yōu)化設(shè)計(jì)

    2.1問題建立

    本文研究對象是一個(gè)基于奇點(diǎn)分布法設(shè)計(jì)的串列泵[9],設(shè)計(jì)流量為0.461m3/s,設(shè)計(jì)揚(yáng)程為13.8m(首級葉輪設(shè)計(jì)揚(yáng)程為5.5m,次級葉輪設(shè)計(jì)揚(yáng)程為8.3m),設(shè)計(jì)效率為83.7%。表1給出了串列泵首、次級葉輪的基本參數(shù)。串列泵三維幾何模型如圖1所示。

    表1 串列泵首、次級葉輪基本參數(shù)

    圖1 串列泵三維幾何模型Fig. 1 3D geometry model of tandem pump

    與普通軸流泵相比,串列式軸流泵的性能影響參數(shù)增多,其中較為顯著的是葉片的安放角,同時(shí)考慮串列泵的兩級葉輪直接連接,首、次級葉片的相對周向位置可能會對串列泵的性能產(chǎn)生影響,所以選擇首、次級葉輪安放角,首、次級葉輪相位角作為設(shè)計(jì)變量。圖2分別給出了3個(gè)設(shè)計(jì)變量的示意圖。在圖2(a)、(b)中,三排葉柵分別表示首次級葉輪和導(dǎo)葉;在圖2(c)中導(dǎo)葉段被省略,兩排葉柵分別表示首次級葉輪。表2列出了本次優(yōu)化過程的設(shè)計(jì)變量及其變化范圍。

    圖2 設(shè)計(jì)變量示意圖Fig. 2 Schematic concept of design variables

    設(shè)計(jì)變量最小值/(°)初始值/(°)最大值/(°)首級葉輪安放角-50+5次級葉輪安放角-50+5相位角-300+30

    選擇效率作為串列泵能量性能的衡量參數(shù)。對于空化性能,考慮到首級葉輪內(nèi)部首先發(fā)生空化,選擇首級葉輪內(nèi)部最小軸向截面平均壓力作為空化性能衡量參數(shù),據(jù)此,定義最小軸向截面平均壓力系數(shù)如下

    (12)

    式中:pav,min是最小軸向截面平均壓力;p∞是環(huán)境壓力;U∞是串列泵進(jìn)口軸向速度。在本文的優(yōu)化設(shè)計(jì)中,原始模型的效率和最小軸面平均壓力系數(shù)的取值分別為83.7%、-1.43。

    2.2實(shí)驗(yàn)設(shè)計(jì)

    圖3 設(shè)計(jì)空間中實(shí)驗(yàn)點(diǎn)分布Fig. 3 Experimental points in design space

    考慮實(shí)驗(yàn)成本,在設(shè)計(jì)實(shí)驗(yàn)過程中,采用拉丁超立方分布和面心分布相結(jié)合的方法共選取45個(gè)實(shí)驗(yàn)點(diǎn),圖3給出了設(shè)計(jì)空間中實(shí)驗(yàn)點(diǎn)的分布情況。

    2.3數(shù)值實(shí)驗(yàn)

    本次優(yōu)化設(shè)計(jì)采用數(shù)值計(jì)算方法獲得實(shí)驗(yàn)點(diǎn)的數(shù)據(jù),即串列泵的效率和全沾濕條件下首級葉輪內(nèi)部最小截面平均壓力系數(shù),數(shù)值方法如下。

    2.3.1控制方程

    連續(xù)性方程和動量方程如下

    (13)

    (14)

    (15)

    (16)

    (17)

    式中:αnuc是空化核子的體積分?jǐn)?shù),取值為5×10-4;RB是空泡直徑,取值為1×10-6m;pv是飽和蒸汽壓,取值為3 169 Pa;Cdest、Cprod分別為凝結(jié)和蒸發(fā)系數(shù),取值為50和 0.01。

    2.3.2湍流模型

    對于全沾濕流動和空化流動,均采用濾波器湍流模型[30]封閉上述方程組,在濾波器湍流模型中,采用和標(biāo)準(zhǔn)k-ε湍流模型相一致的k方程和ε方程采用的標(biāo)準(zhǔn)方程形式。其湍流粘性μt定義如下

    (18)

    式中:F為濾波函數(shù),由濾波尺寸λ和湍流特征長度的比值決定:

    (19)

    在標(biāo)準(zhǔn)k-ε模型中加入濾波函數(shù)后,對尺度小于濾波器尺寸的湍流,采用標(biāo)準(zhǔn)k-ε模型計(jì)算來求解,對尺度大于濾波器尺寸的湍流結(jié)構(gòu),采用直接計(jì)算方法來求解。本文計(jì)算采用ANSYS-CFX 12.1商業(yè)計(jì)算軟件,通過二次開發(fā)將上述的FBM模型引入,濾波尺寸選擇最小網(wǎng)格尺寸的1.5倍。

    2.3.3網(wǎng)格和邊界條件設(shè)置

    采用全結(jié)構(gòu)化網(wǎng)格劃分串列泵內(nèi)部流場。為了提高網(wǎng)格質(zhì)量,將整個(gè)計(jì)算域分成4個(gè)部分,分別為入口段、首級葉輪、次級葉輪和導(dǎo)葉段。不同部分采用不同的網(wǎng)格劃分方法:入口段采用商業(yè)軟件ANSYS-ICEM 12.1商業(yè)軟件劃分流場網(wǎng)格,首級葉輪、次級葉輪和導(dǎo)葉段采用商業(yè)軟件ANSYS-TurboGrid 12.1劃分流場網(wǎng)格,導(dǎo)水錐和葉片周圍采用O型網(wǎng)格;近壁面進(jìn)行加密,保證y+=ρlΔyuτ/μl≈1(其中Δy為距離壁面最近一層網(wǎng)格厚度,uτ為壁面摩擦速度)。輪緣間隙內(nèi)包含10個(gè)節(jié)點(diǎn)。

    為了驗(yàn)證網(wǎng)格對結(jié)果的影響,對不同節(jié)點(diǎn)數(shù)的網(wǎng)格進(jìn)行計(jì)算。圖4給出了網(wǎng)格節(jié)點(diǎn)數(shù)目和串列泵能量參數(shù)之間的關(guān)系,從圖中可以看出,效率和揚(yáng)程系數(shù)在網(wǎng)格節(jié)點(diǎn)數(shù)為2.21×106后基本不變。綜合考慮計(jì)算的精度和經(jīng)濟(jì)性,最終確定的計(jì)算網(wǎng)格總節(jié)點(diǎn)數(shù)為2.21×106。圖中,揚(yáng)程系數(shù)和效率定義為

    (20)

    (21)

    式中:H為串列泵揚(yáng)程,P為串列泵軸功率。

    圖4 網(wǎng)格無關(guān)性驗(yàn)證Fig. 4 Effect of grid refinement

    計(jì)算域的進(jìn)口采用總壓進(jìn)口條件;出口設(shè)置質(zhì)量流量出口條件;固壁采用無滑移邊界條件;壁面函數(shù)為標(biāo)準(zhǔn)壁面函數(shù);動葉輪和靜葉輪交接面類型選擇Stage界面,動葉輪之間的交接面類型選擇Frozen Rotor界面。通過監(jiān)測各參量的殘差精度來確定計(jì)算的收斂性,收斂標(biāo)準(zhǔn)是10-4。對于空化流動計(jì)算,通過調(diào)節(jié)不同的進(jìn)口壓力獲得不同的裝置有效空化余量取值;出口設(shè)置質(zhì)量流量條件。計(jì)算采用相同的網(wǎng)格和邊界條件,以全沾濕工況計(jì)算結(jié)果為初始值,進(jìn)行迭代計(jì)算,通過監(jiān)測各參量的殘差精度來確定計(jì)算的收斂性,收斂標(biāo)準(zhǔn)是10-4。

    2.3.4計(jì)算方法驗(yàn)證

    利用上述計(jì)算方法針對一個(gè)比轉(zhuǎn)速為ns=700的單級軸流泵的外特性和設(shè)計(jì)流量下的空化特性進(jìn)行預(yù)測。圖5給出了計(jì)算結(jié)果和實(shí)驗(yàn)結(jié)果的對比情況。從圖中可以看出,本文所采用的計(jì)算方法可以很好地預(yù)測該軸流泵的外特性和空化特性。

    (a)外特性                   (b)空化特性圖5 單級泵外特性和空化特性曲線Fig. 5 Performances of a single stage pump

    2.4模型建立及驗(yàn)證

    通過數(shù)值計(jì)算可以獲得不同實(shí)驗(yàn)點(diǎn)即不同模型泵的效率和首級葉輪內(nèi)部最小軸截面平均壓力系數(shù)。根據(jù)計(jì)算所得的目標(biāo)函數(shù)值,分別采用二階多項(xiàng)式法,克里金法和徑向神經(jīng)網(wǎng)絡(luò)法建立代理模型。 圖6給出了3種代理模型的三維等值面分布圖,3個(gè)坐標(biāo)軸分別表示3個(gè)設(shè)計(jì)變量,不同顏色對應(yīng)于不同的目標(biāo)函數(shù)取值。

    從圖中可以看到,采用二階多項(xiàng)式法和克里金法建立的代理模型預(yù)測結(jié)果相近,而采用徑向神經(jīng)網(wǎng)絡(luò)法建立的代理模型和前兩種方法相差較大。

    采用交叉驗(yàn)證方法對3種模型進(jìn)行評價(jià),結(jié)果如表3所示。從表中可以看出三種模型的誤差值均超過10%,所以有必要進(jìn)一步修正設(shè)計(jì)空間,重新設(shè)計(jì)實(shí)驗(yàn)點(diǎn)分布。但可以進(jìn)行初步分析,從表中可以看出誤差相對較小的是二階多項(xiàng)式方法建立的代理模型,所以初步分析基于二階多項(xiàng)式代理模型結(jié)果。

    (b)克里金法(KRG)

    (c)徑向神經(jīng)網(wǎng)絡(luò)法(RBNN)圖6 不同代理模型預(yù)測目標(biāo)函數(shù)三維等值面分布Fig. 6 Predictions of 3D contours of objective functions using different surrogate models

    PRSKRGRBNN效率/%10.8812.5358.30最小平均壓力系數(shù)/%43.8343.2557.36

    2.5結(jié)果分析與討論

    圖7給出了基于二階多項(xiàng)式模型的敏感度分析結(jié)果。圖7(a)是3個(gè)自變量對于最小截面平均壓力系數(shù)的敏感度分析結(jié)果,圖7(b)是3個(gè)自變量對于效率的敏感度分析結(jié)果。圖中空心框表示局部敏感度,實(shí)心框表示全局敏感度。結(jié)果表明,首、次級葉輪安放角對效率和最小平均壓力系數(shù)影響較大,而相位角的影響很小,和前兩者相比可以忽略。在修正設(shè)計(jì)空間時(shí),可以去掉相位角,從而降低設(shè)計(jì)空間的維數(shù),使二次優(yōu)化過程得到簡化。

    圖8給出了基于二階多項(xiàng)式代理模型計(jì)算結(jié)果的最優(yōu)點(diǎn)選取情況,目標(biāo)函數(shù)空間中橫坐標(biāo)為最小截面平均壓力系數(shù),縱坐標(biāo)為效率。綜合考慮效率和最小截面平均壓力系數(shù)2種因素,同時(shí)考慮到代理模型預(yù)測精度較低,所以選擇目標(biāo)函數(shù)空間里的修正空間范圍較大,位置如圖中線框表示,即最優(yōu)點(diǎn)應(yīng)當(dāng)落在該范圍內(nèi)。

    圖7 敏感度分析Fig. 7 Global sensitivity analysis

    圖8 目標(biāo)空間代理模型預(yù)測結(jié)果Fig. 8 Predictions in objective space

    將上述修正空間映射到設(shè)計(jì)空間內(nèi),可以得到最優(yōu)點(diǎn)對應(yīng)的各個(gè)設(shè)計(jì)變量的取值范圍如下(已忽略相位角):首級葉輪安放角在-1~-8°,次級葉輪安放角在-1~-6°?;诖?,在此范圍內(nèi)重新建立設(shè)計(jì)空間,進(jìn)行二次優(yōu)化。

    2.6二次優(yōu)化

    二次優(yōu)化的目標(biāo)函數(shù)不變,由于相位角對串列泵性能影響較小,所以在二次優(yōu)化時(shí)忽略。表4給出了二次優(yōu)化的設(shè)計(jì)變量。

    表4 二次優(yōu)化設(shè)計(jì)變量

    二次優(yōu)化設(shè)計(jì)實(shí)驗(yàn)過程中,采用拉丁超立方分布方法在修正的設(shè)計(jì)空間內(nèi)選取20個(gè)實(shí)驗(yàn)點(diǎn)。采用數(shù)值計(jì)算方法獲取實(shí)驗(yàn)點(diǎn)數(shù)據(jù),然后分別采用二階多項(xiàng)式法,克里金法和徑向神經(jīng)網(wǎng)絡(luò)法建立代理模型。

    圖9給出了3種代理模型預(yù)測結(jié)果的三維示意圖,圖中2個(gè)水平面內(nèi)的坐標(biāo)軸分別表示二次優(yōu)化的2個(gè)設(shè)計(jì)變量,豎直坐標(biāo)軸表示目標(biāo)函數(shù),不同顏色表示不同的目標(biāo)函數(shù)取值。對3種模型進(jìn)行驗(yàn)證,結(jié)果如表5所示。從表中可以看出二次優(yōu)化的代理模型誤差大大減小,表明建立的代理模型具有足夠的精度。從表中還可以看出,二階多項(xiàng)式方法建立的代理模型誤差最小,基于該代理模型進(jìn)行最優(yōu)點(diǎn)的選取以及進(jìn)一步的分析。

    表5二次代理模型驗(yàn)證結(jié)果

    Table 5PRESS of surrogate models of the second optimization

    PRSKRGRBNN效率/%0.342.011.59最小平均壓力系數(shù)/%0.170.400.24

    基于二階多項(xiàng)式方法建立的代理模型對各設(shè)計(jì)參數(shù)進(jìn)行敏感度分析,圖10分別給出了2個(gè)自變量對于最小截面平均壓力系數(shù)和效率的敏感度分析結(jié)果。

    從圖10(a)中可以看出首次級葉輪葉片安放角對最小截面平均壓力系數(shù)的敏感度較高,表明首級葉輪葉片安放角和串列泵的空化性能關(guān)系密切,2個(gè)自變量的全局敏感度均大于局部敏感度,表明2個(gè)設(shè)計(jì)參數(shù)存在相互影響,表明次級葉輪葉片安放角也會對串列泵的空化性能產(chǎn)生影響,但影響力較小。從圖10(b)中可以看出,首、次級葉輪葉片安放角均會對效率產(chǎn)生影響,但次級葉輪葉片安放角占據(jù)主導(dǎo)地位,首、次級葉輪葉片安放角對效率的敏感度之比約為4∶6,與串列泵的設(shè)計(jì)參數(shù):首、次級載荷之比相同。下面基于二階多項(xiàng)式模型的預(yù)測結(jié)果尋找最優(yōu)工況。

    (a)二階多項(xiàng)式法(PRS)

    (b)克里金法(KRG)

    (c)徑向神經(jīng)網(wǎng)絡(luò)法(RBNN)圖9 二次優(yōu)化不同代理模型預(yù)測的三維響應(yīng)面Fig. 9 Predictions of 3D response surfaces using different surrogate models in the second optimization

    圖10 二次優(yōu)化敏感度分析Fig. 10 Global sensitivity analysis of the second optimization

    圖11給出了目標(biāo)函數(shù)空間內(nèi)的代理模型預(yù)測的目標(biāo)函數(shù)點(diǎn)分布情況。顏色加深的點(diǎn)對應(yīng)于效率最高或者最小截面壓力系數(shù)最大的工況,這些工況點(diǎn)形成帕氏最優(yōu)狀態(tài)(Pareto front)[32],從帕氏最優(yōu)狀態(tài)可以看到,效率和最小壓力系數(shù)是相互矛盾的性能參數(shù)(trade-off),即對于模型泵來說,如果效率較高,則最小截面壓力系數(shù)較小,即空化性能較差,如果最小壓力系數(shù)較高,即保證了較好的空化性能,則效率不能達(dá)到較高水平。

    綜合考慮2種因素,選擇最優(yōu)點(diǎn)如圖11所示。映射到設(shè)計(jì)空間中可以找到最優(yōu)點(diǎn)的設(shè)計(jì)變量取值為:首級葉輪-3.7°,次級葉輪葉片安放角為-2°。采用數(shù)值方法對優(yōu)化結(jié)果進(jìn)行初步驗(yàn)證可以得到效率和對小壓力系數(shù)分別為85.5%和-1.31。

    圖11 目標(biāo)函數(shù)空間內(nèi)的預(yù)測結(jié)果和帕氏最優(yōu)邊界Fig. 11 Predictions and Pareto front in objective space

    2.7優(yōu)化結(jié)果驗(yàn)證

    為了對優(yōu)化結(jié)果進(jìn)行進(jìn)一步驗(yàn)證,圖12給出了采用數(shù)值計(jì)算方法得到的優(yōu)化前和優(yōu)化后的串列泵外特性和空化特性曲線。從圖中可以看出,優(yōu)化設(shè)計(jì)后的串列泵設(shè)計(jì)點(diǎn)附近的效率有所提升,同時(shí)設(shè)計(jì)點(diǎn)的空化特性得到改善。表明該優(yōu)化方法可以很好地應(yīng)用于串列泵的優(yōu)化設(shè)計(jì)。

    圖13給出了優(yōu)化前后在設(shè)計(jì)流量,無空化條件下的串列泵首次級葉輪子午面平均壓力云圖。從圖中可以看出,首級葉輪經(jīng)過優(yōu)化后,進(jìn)口處的低壓區(qū)域面積有所減小,表明在設(shè)計(jì)流量下,進(jìn)口處流動損失減小,空化性能得到改善。次級葉輪經(jīng)過優(yōu)化后,沿葉高方向壓力梯度減小,表明優(yōu)化后的葉輪可以抑制次級葉輪內(nèi)部的二次流動,減小流動損失。

    (a)外特性               (b)空化特性(設(shè)計(jì)流量)圖12 優(yōu)化模型和原模型對比Fig. 12 Comparisons of prototype and optimum pump

    圖13 優(yōu)化前后串列泵子午面壓力圖對比(設(shè)計(jì)流量)Fig. 13 Average static pressure contour of meridian interface in front and rear impellers at design flow rate

    3結(jié)論

    本文基于代理模型對串列泵進(jìn)行優(yōu)化設(shè)計(jì),分析了首、次級葉輪進(jìn)口安放角度和相位角對串列泵性能的影響。采用敏感度分析方法和Pareto front方法進(jìn)行參數(shù)影響分析和最優(yōu)點(diǎn)的選取,并對優(yōu)化后的串列泵外特性和空化特性進(jìn)行驗(yàn)證,得到的主要結(jié)論如下:

    1) 基于代理模型的系統(tǒng)優(yōu)化設(shè)計(jì)方法可以較好地應(yīng)用于串列泵的優(yōu)化設(shè)計(jì);優(yōu)化結(jié)果可以提升串列泵設(shè)計(jì)流量附近的效率和空化性能。

    2) 首、次級葉輪相位角對串列泵性能影響較??;首、次級葉輪葉片安放角對串列泵效率的局部敏感度之比約為4∶6,和設(shè)計(jì)的載荷之比相同,首級葉輪葉片安放角是串列泵的空化特性的主要影響因素。

    參考文獻(xiàn):

    [1]SAHA U K, ROY B. Experimental investigations on tandem compressor cascade performance at low speeds[J]. Experimental thermal and fluid science, 1997, 14(3): 263-276.

    [2]FEINDT E G. Exact calculation of the flow of a staggered tandem cascade with moving second blade row[J]. Archive of applied mechanics, 2001, 71(9): 589-600.

    [3]QIU X, DANG T. 3D inverse method for turbomachine blading with splitter blades[C]//International Gas Turbine and Aeroengine Congress and Exhibition. ASME paper 2000-GT-0526, 2000: 8-11.

    [4]LIU Pengfei, ISLAM M, VEITCH B. Unsteady hydromechanics of a steering podded propeller unit[J]. Ocean engineering, 2009, 36(12/13): 1003-1014.

    [5]孫琴, 顧蘊(yùn)德, 鄭淑珍. 串列螺旋槳及其設(shè)計(jì)方法[M]. 北京: 人民交通出版社, 1985.

    SUN Qin, GU Yunde, ZHENG Shuzhen. Tandem propeller and design method[M]. Beijing: China Communication Press, 1985.

    [6]王國玉, 陳榮鑫, 余志毅, 等. 串列式雙級軸流泵性能的數(shù)值模擬[J]. 機(jī)械工程學(xué)報(bào), 2007, 43(12): 39-45.

    WANG Guoyu, CHEN Rongxin, YU Zhiyi, et al. Numerical simulation of the performance of tandem axial flow pump[J]. Journal of mechanical engineering, 2007, 43(12): 39-45.

    [7]趙宇, 王國玉, 白澤宇, 等. 串列式軸流泵首次級葉輪能量特性及相互作用分析[J]. 排灌機(jī)械工程學(xué)報(bào), 2013, 31(3): 210-214.

    ZHAO Yu, WANG Guoyu, BAI Zeyu, et al. Characteristics and interaction of front and rear impellers of axial flow tandem pump[J]. Journal of drainage and irrigation machinery engineering, 2013, 31(3): 210-214.

    [8]張克危. 流體機(jī)械原理[M]. 北京: 機(jī)械工業(yè)出版社, 2000.

    ZHANG Kewei. Theory of hydraulic machine[M]. Beijing: China Machine Press, 2000.

    [9]趙宇, 王國玉, 吳欽, 等. 基于計(jì)算流體力學(xué)的串列軸流泵空化性能分析[J]. 機(jī)械工程學(xué)報(bào), 2014, 50(6): 171-176, 184.

    ZHAO Yu, WANG Guoyu, WU Qin, et al. Analysis of cavitation performances of an axial flow tandem pump based on computational fluid dynamics[J]. Journal of mechanical engineering, 2014, 50(6): 171-176, 184.

    [10]王立祥, 傅健. 關(guān)于噴水推進(jìn)串列式軸流泵葉輪參數(shù)選擇與計(jì)算的探討[J]. 船舶, 2003(6): 47-51.

    WANG Lixiang, FU Jian. Parameter selection and calculation of tandem axial flow waterjet pump impeller[J]. Ship & boat, 2003(6): 47-51.

    [11]傅健. 噴水推進(jìn)串列式軸流泵水力性能的探索與研究[D]. 上海: 中國艦船研究院, 2004.

    FU Jian. A Study on hydrodynamic performance of tandem axial flow waterjet pump[D]. Shanghai: Marine Design and Research Institute of China, 2004.

    [12]YU Zhiyi, LIU Shuyan, WANG Guoyu. Numerical investigation of the performance of an axial-flow pump with tandem blades[J]. Journal of Beijing institute of technology, 2007, 16(4): 404-408.

    [13]SHYY W, PAPILA N, VAIDYANATHAN R, et al. Global design optimization for aerodynamics and rocket propulsion components[J]. Progress in aerospace sciences, 2001, 37(1): 59-118.

    [14]QUEIPO N V, HAFTKA R T, SHYY W, et al. Surrogate-based analysis and optimization[J]. Progress in aerospace sciences, 2005, 41(1): 1-28.

    [15]BOOKER A J, DENNIS J E Jr, FRANK P D, et al. Optimization using surrogate objectives on a helicopter test example[M]//BORGGAARD J, BURNS J, CLIFF E, et al. eds. Computational Methods for Optimal Design and Control. Boston: Birkhuser Boston, 1998: 49-58.

    [16]RAI M M, MADAVAN N K, HUBER F W. Improving the unsteady aerodynamic performance of transonic turbines using neural networks[C]//Proceedings of the 38th AIAA Aerospace Sciences Meeting and Exhibit. Reno, NV, 2000: 2000-0169.

    [17]MADSEN J I, SHYY W, HAFTKA R T. Response surface techniques for diffuser shape optimization[J]. AIAA journal, 2000, 38(9): 1512-1518.

    [18]PAPILA N, SHYY W, GRIFFIN L, et al. Shape optimization of supersonic turbines using global approximation methods[J]. Journal of propulsion and power, 2002, 18(3): 509-518.

    [19]VAIDYANATHAN R, TUCKER K, PAPILA N, et al. CFD-based design optimization for single element rocket injector[M]. Florida: University of Florida Department of Mechanical and Aerospace Engineering, 2003.

    [20]ANAGNOSTOPOULOS J S. A fast numerical method for flow analysis and blade design in centrifugal pump impellers[J]. Computers & fluids, 2009, 38(2): 284-289.

    [21]ZHANG Jinya, ZHU Hongwu, CHUN Yang, et al. Multi-objective shape optimization of helico-axial multiphase pump impeller based on NSGA-II and ANN[J]. Energy conversion and management, 2011, 52(1): 538-546.

    [22]FAN J, EVES J, THOMPSON H M, et al. Mincher. Computational fluid dynamic analysis and design optimization of jet pumps[J]. Computers & fluids, 2011, 46(1): 212-217.

    [23]KIM J H, KIM K Y. Hydrodynamic performance enhancement of a mixed-flow pump[C]//Proceedings of the 26th IAHR Symposium on Hydraulic Machinery and Systems. Beijing, China, 2012.

    [24]G?LCU M. Artificial neural network based modeling of performance characteristics of deep well pumps with splitter blade[J]. Energy conversion and management, 2006, 47(18/19): 3333-3343.

    [25]MCKAY M D, BECKMAN R J, CONOVER W J. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code[J]. Technometrics, 2000, 42(1): 55-61.

    [26]MYERS R H, MONTGOMERY D C, ANDERSON-COOK C M. Response surface methodology: process and product optimization using designed experiments[M]. 3rd ed. Wiley, 2009.

    [27] SACKS J, WELCH W J, MITCHELL T J, et al. Design and analysis of computer experiments[J]. Statistical science, 1989, 4(4): 409-423.

    [28]KOHAVI R. A study of cross-validation and bootstrap for accuracy estimation and model selection[C]//Proceedings of the 14th International Joint Conference on Artificial Intelligence. San Francisco, CA, USA, 1995, 2: 1137-1143.

    [29]SOBOL′ I M. Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[J]. Mathematics and computers in simulation, 2001, 55(1-3): 271-280.

    [30]JOHANSEN S T, WU Jiongyang, SHYY W. Filter-based unsteady RANS computations[J]. International journal of heat and fluid flow, 2004, 25(1): 10-21.

    [31]KUBOTA A, KATO H, YAMAGUCHI H. A new modelling of cavitating flows: A numerical study of unsteady cavitation on a hydrofoil section[J]. Journal of fluid mechanics, 1992, 240: 59-96.

    [32]HORN J, NAFPLIOTIS N, GOLDBERG D E. A niched Pareto genetic algorithm for multiobjective optimization[C]//Proceedings of the First IEEE Conference on Evolutionary Computation, 1994. IEEE World Congress on Computational Intelligence. Piscataway, NJ, 1994: 82-87.

    Optimization design of tandem pump based on surrogate method

    ZHAO Yu, WANG Guoyu, HUANG Biao, WANG Fufeng

    (Research Institute of Fluid Mechanical Engineering, Beijing Institute of Technology, Beijing 100081, China)

    Abstract:In this paper, we propose a global design optimization method based on a surrogate model for an axial-flow tandem pump, in order to overcome the disadvantages of conventional optimization methods, including their long design period and lower efficiency. We discuss the influences of the primary design parameters on pump performances, and used fixed blade angles in both the front and rear impellers, as well as the phase angle, as design variables. We selected efficiency and the minimum average pressure coefficient on the axial sectional surface as the objective functions, which represent the energy and cavitation performances, respectively. Surrogate models were constructed based on various methods. In addition, we used global sensitivity analysis and Pareto front methods to further analyze the design parameters and optimum point. The results show that the optimization results can enhance tradeoff performances well, with respect to both efficiency and cavitation performance, in the vicinity of the design flow of the tandem pump. The influence of phase angle on performance can be neglected relative to the influence of the other two design variables. The impact ratio on efficiency of a fixed blade angle in the two impellers is the same as in their designed loading distributions. A fixed blade angle in the front impeller has a great influence on cavitation performance.

    Keywords:axial-flow pump; tandem cascades; global design optimization; surrogate method; cavitation performances

    中圖分類號:TV131.32

    文獻(xiàn)標(biāo)志碼:A

    文章編號:1006-7043(2016)03-438-11

    doi:10.11990/jheu.201410085

    作者簡介:趙宇(1989-), 男, 博士研究生;通信作者:王國玉, E-mail: wangguoyu@bit.edu.cn.

    基金項(xiàng)目:國家自然科學(xué)基金資助項(xiàng)目(51239005, 51479002).

    收稿日期:2014-10-31.

    網(wǎng)絡(luò)出版地址:http://www.cnki.net/kcms/detail/23.1390.u.20160104.1648.020.html

    網(wǎng)絡(luò)出版日期:2016-01-04.

    王國玉(1961-), 男, 教授, 博士生導(dǎo)師.

    猜你喜歡
    軸流泵
    潛水軸流泵在工農(nóng)河北排澇站中的應(yīng)用
    大科技(2023年48期)2023-12-20 15:45:07
    潛水軸流泵運(yùn)行故障分析與排除研究
    空化對軸流泵葉輪流固耦合特性的影響
    潛水軸流泵電機(jī)運(yùn)行工況的特點(diǎn)及可靠性探討
    基于數(shù)值模擬的軸流泵效率分析
    固定槳式軸流泵外特性調(diào)節(jié)方法的研究
    對旋式軸流泵性能CFD計(jì)算分析
    無后置導(dǎo)葉軸流泵瞬態(tài)空化流場的數(shù)值分析
    中、小型軸流泵安裝方法分析研究
    軸流泵工況調(diào)節(jié)方法及其特性分析
    亚洲视频免费观看视频| 免费高清在线观看日韩| 欧美精品一区二区免费开放| 首页视频小说图片口味搜索| 啦啦啦中文免费视频观看日本| a级毛片黄视频| 99国产极品粉嫩在线观看| 80岁老熟妇乱子伦牲交| 亚洲精品中文字幕一二三四区 | 女人高潮潮喷娇喘18禁视频| 国产欧美日韩精品亚洲av| 久久久久久久久久久久大奶| 夜夜夜夜夜久久久久| 99精品在免费线老司机午夜| 精品国产超薄肉色丝袜足j| 国产成人系列免费观看| 亚洲人成77777在线视频| 欧美日韩视频精品一区| 淫妇啪啪啪对白视频| 777久久人妻少妇嫩草av网站| 免费看a级黄色片| 黄色视频不卡| 国产欧美日韩综合在线一区二区| 窝窝影院91人妻| 国产亚洲一区二区精品| 亚洲精品一二三| 亚洲国产精品一区二区三区在线| 日本撒尿小便嘘嘘汇集6| 高潮久久久久久久久久久不卡| 国产精品1区2区在线观看. | 国产午夜精品久久久久久| 国产黄色免费在线视频| 国产aⅴ精品一区二区三区波| 一本综合久久免费| 精品亚洲乱码少妇综合久久| 人人妻人人澡人人爽人人夜夜| 大香蕉久久成人网| 欧美国产精品一级二级三级| 一级a爱视频在线免费观看| 巨乳人妻的诱惑在线观看| 国产一区二区三区在线臀色熟女 | 久久久久久久久久久久大奶| 亚洲精品国产色婷婷电影| 久久久久网色| 人妻 亚洲 视频| 肉色欧美久久久久久久蜜桃| 精品午夜福利视频在线观看一区 | 国产欧美日韩精品亚洲av| 大香蕉久久成人网| 99精品在免费线老司机午夜| 嫩草影视91久久| 欧美激情极品国产一区二区三区| 欧美日韩一级在线毛片| 色综合婷婷激情| 国产欧美日韩一区二区三区在线| 国产成人精品在线电影| 桃花免费在线播放| 激情在线观看视频在线高清 | 久久人妻熟女aⅴ| 亚洲伊人久久精品综合| 欧美人与性动交α欧美软件| 久久99一区二区三区| av片东京热男人的天堂| 免费日韩欧美在线观看| 一区福利在线观看| 人成视频在线观看免费观看| 午夜久久久在线观看| 色老头精品视频在线观看| 国产真人三级小视频在线观看| 亚洲九九香蕉| 首页视频小说图片口味搜索| 老司机在亚洲福利影院| 国产高清激情床上av| 欧美在线一区亚洲| av天堂久久9| 女同久久另类99精品国产91| 亚洲成av片中文字幕在线观看| 超色免费av| 国产aⅴ精品一区二区三区波| 动漫黄色视频在线观看| 岛国毛片在线播放| 国产精品美女特级片免费视频播放器 | 高清毛片免费观看视频网站 | 免费人妻精品一区二区三区视频| 最新美女视频免费是黄的| 变态另类成人亚洲欧美熟女 | 日日摸夜夜添夜夜添小说| 老司机靠b影院| 一二三四社区在线视频社区8| 欧美日韩成人在线一区二区| 久久影院123| 少妇 在线观看| 久热这里只有精品99| 纵有疾风起免费观看全集完整版| 久久ye,这里只有精品| 久久国产精品影院| 一二三四社区在线视频社区8| 欧美精品一区二区免费开放| 99精品在免费线老司机午夜| 乱人伦中国视频| 露出奶头的视频| 极品人妻少妇av视频| 亚洲中文日韩欧美视频| 十八禁高潮呻吟视频| 黄色视频,在线免费观看| 性少妇av在线| 在线 av 中文字幕| 国产主播在线观看一区二区| 久久性视频一级片| 亚洲中文日韩欧美视频| 日韩 欧美 亚洲 中文字幕| 高清黄色对白视频在线免费看| 国产激情久久老熟女| 亚洲午夜精品一区,二区,三区| 欧美黑人精品巨大| 别揉我奶头~嗯~啊~动态视频| 99精品在免费线老司机午夜| 色在线成人网| 久久精品91无色码中文字幕| 精品少妇久久久久久888优播| 国产高清视频在线播放一区| 国产精品自产拍在线观看55亚洲 | 国产在视频线精品| 亚洲视频免费观看视频| 三级毛片av免费| 欧美乱码精品一区二区三区| 日韩欧美一区二区三区在线观看 | 极品人妻少妇av视频| 午夜福利影视在线免费观看| 国产精品国产av在线观看| 十八禁人妻一区二区| 人人澡人人妻人| 黄色 视频免费看| 岛国毛片在线播放| 国产亚洲欧美精品永久| 久久久久久免费高清国产稀缺| 捣出白浆h1v1| 久久这里只有精品19| 夜夜骑夜夜射夜夜干| 精品视频人人做人人爽| 久久香蕉激情| 岛国在线观看网站| 他把我摸到了高潮在线观看 | 老司机亚洲免费影院| 亚洲成人免费av在线播放| 老汉色av国产亚洲站长工具| 麻豆av在线久日| 日本av手机在线免费观看| 日韩免费av在线播放| 久久精品人人爽人人爽视色| 免费av中文字幕在线| 成人特级黄色片久久久久久久 | 最新美女视频免费是黄的| 女人被躁到高潮嗷嗷叫费观| 精品国产一区二区三区久久久樱花| 最新在线观看一区二区三区| 激情视频va一区二区三区| 国产在线视频一区二区| 18禁裸乳无遮挡动漫免费视频| 在线十欧美十亚洲十日本专区| 亚洲精品久久成人aⅴ小说| a在线观看视频网站| 啦啦啦中文免费视频观看日本| 午夜免费成人在线视频| 亚洲av欧美aⅴ国产| 国产精品一区二区在线不卡| 成人av一区二区三区在线看| 国产成人精品久久二区二区91| 两个人看的免费小视频| 成人国产一区最新在线观看| a级毛片在线看网站| 99精品在免费线老司机午夜| 80岁老熟妇乱子伦牲交| 国产97色在线日韩免费| 免费在线观看影片大全网站| 交换朋友夫妻互换小说| 国产精品电影一区二区三区 | 久久精品国产亚洲av高清一级| 国产成人啪精品午夜网站| 亚洲精品在线美女| 18禁美女被吸乳视频| 丰满人妻熟妇乱又伦精品不卡| 激情视频va一区二区三区| 成人亚洲精品一区在线观看| 伦理电影免费视频| 18禁观看日本| 国产亚洲精品久久久久5区| 久久久久久久精品吃奶| 黄色怎么调成土黄色| 他把我摸到了高潮在线观看 | 狠狠精品人妻久久久久久综合| 国产精品av久久久久免费| 国产精品 欧美亚洲| 婷婷成人精品国产| 日韩欧美免费精品| 天堂8中文在线网| 亚洲少妇的诱惑av| 欧美成狂野欧美在线观看| 午夜日韩欧美国产| 一区福利在线观看| 亚洲熟女精品中文字幕| 成年人午夜在线观看视频| 午夜福利,免费看| 国产一区二区在线观看av| 免费日韩欧美在线观看| 日韩熟女老妇一区二区性免费视频| av又黄又爽大尺度在线免费看| 日本黄色视频三级网站网址 | 国产成人免费观看mmmm| 欧美日韩黄片免| 国产av精品麻豆| 国产精品 欧美亚洲| 老熟妇仑乱视频hdxx| 亚洲欧美一区二区三区黑人| 日日夜夜操网爽| 在线观看一区二区三区激情| 欧美日韩福利视频一区二区| 国产精品一区二区精品视频观看| 日韩熟女老妇一区二区性免费视频| 午夜激情久久久久久久| 好男人电影高清在线观看| 丝袜美腿诱惑在线| 午夜免费鲁丝| 欧美变态另类bdsm刘玥| 久久精品国产99精品国产亚洲性色 | 亚洲国产成人一精品久久久| 最近最新免费中文字幕在线| 中文字幕精品免费在线观看视频| 人妻一区二区av| 国产97色在线日韩免费| 久久久久久久大尺度免费视频| 天堂中文最新版在线下载| 夫妻午夜视频| 成人国产一区最新在线观看| 国产色视频综合| 免费在线观看视频国产中文字幕亚洲| 午夜福利,免费看| 日韩大片免费观看网站| 天天躁夜夜躁狠狠躁躁| 亚洲一区二区三区欧美精品| 又黄又粗又硬又大视频| 日本黄色日本黄色录像| 免费av中文字幕在线| 日韩制服丝袜自拍偷拍| 夜夜爽天天搞| 9191精品国产免费久久| 亚洲久久久国产精品| 久久久国产一区二区| 日韩免费av在线播放| 久久久久久久久久久久大奶| 夜夜骑夜夜射夜夜干| 99国产精品免费福利视频| 亚洲精品乱久久久久久| 蜜桃国产av成人99| 母亲3免费完整高清在线观看| 国产精品电影一区二区三区 | 久久精品91无色码中文字幕| 精品久久久精品久久久| 在线天堂中文资源库| 久久久国产一区二区| 亚洲全国av大片| 亚洲精华国产精华精| 欧美日韩国产mv在线观看视频| 1024香蕉在线观看| 亚洲精品国产一区二区精华液| 成人国产一区最新在线观看| 麻豆国产av国片精品| 曰老女人黄片| 国内毛片毛片毛片毛片毛片| 中文字幕av电影在线播放| 欧美黑人精品巨大| 亚洲国产欧美在线一区| 午夜精品国产一区二区电影| 桃红色精品国产亚洲av| 下体分泌物呈黄色| 黄频高清免费视频| 咕卡用的链子| 国产亚洲欧美在线一区二区| 黄色怎么调成土黄色| 精品一品国产午夜福利视频| 最新在线观看一区二区三区| 日韩欧美国产一区二区入口| 久久香蕉激情| 999久久久精品免费观看国产| 一夜夜www| 久久 成人 亚洲| 久久中文看片网| 亚洲综合色网址| 丝瓜视频免费看黄片| 在线观看免费日韩欧美大片| 麻豆成人av在线观看| 欧美日韩亚洲综合一区二区三区_| 亚洲精品av麻豆狂野| 久久毛片免费看一区二区三区| 18禁裸乳无遮挡动漫免费视频| 成人av一区二区三区在线看| 欧美成狂野欧美在线观看| 曰老女人黄片| 天天操日日干夜夜撸| 精品少妇内射三级| 美女高潮喷水抽搐中文字幕| 国产精品欧美亚洲77777| av天堂在线播放| 亚洲熟女精品中文字幕| 久久久精品区二区三区| 十八禁人妻一区二区| 在线 av 中文字幕| 国产免费福利视频在线观看| 91九色精品人成在线观看| 亚洲熟妇熟女久久| 日本av手机在线免费观看| 国产av一区二区精品久久| 少妇的丰满在线观看| 嫩草影视91久久| 国产av又大| 欧美日韩精品网址| cao死你这个sao货| 热99国产精品久久久久久7| av有码第一页| av国产精品久久久久影院| 欧美老熟妇乱子伦牲交| 久久午夜综合久久蜜桃| 蜜桃在线观看..| 亚洲综合色网址| 国产亚洲精品第一综合不卡| 免费不卡黄色视频| 成人国产一区最新在线观看| 国产麻豆69| 另类精品久久| 两性夫妻黄色片| 97人妻天天添夜夜摸| 男女之事视频高清在线观看| 亚洲欧美一区二区三区久久| 亚洲精品粉嫩美女一区| 黄片播放在线免费| 考比视频在线观看| 成人18禁在线播放| 日韩欧美一区二区三区在线观看 | 一级毛片女人18水好多| 一本一本久久a久久精品综合妖精| 男女高潮啪啪啪动态图| 午夜福利乱码中文字幕| 国产亚洲午夜精品一区二区久久| h视频一区二区三区| 亚洲天堂av无毛| 侵犯人妻中文字幕一二三四区| 男女免费视频国产| 老司机在亚洲福利影院| 国产精品 国内视频| 欧美精品啪啪一区二区三区| 男女之事视频高清在线观看| 欧美 亚洲 国产 日韩一| 国产精品自产拍在线观看55亚洲 | 99久久人妻综合| 精品国产乱码久久久久久小说| 久久久精品区二区三区| 久久午夜亚洲精品久久| 久久婷婷成人综合色麻豆| 91国产中文字幕| 欧美成人午夜精品| 亚洲精品自拍成人| 黑人巨大精品欧美一区二区mp4| 天天操日日干夜夜撸| 超碰成人久久| 国产伦人伦偷精品视频| 老司机亚洲免费影院| 99国产综合亚洲精品| 久久热在线av| 亚洲成人国产一区在线观看| 在线播放国产精品三级| 宅男免费午夜| 9色porny在线观看| 久久天堂一区二区三区四区| 两个人看的免费小视频| 99国产精品免费福利视频| 黄片播放在线免费| 欧美日本中文国产一区发布| 国产成人免费无遮挡视频| 女人久久www免费人成看片| 脱女人内裤的视频| 日韩熟女老妇一区二区性免费视频| 一级毛片精品| 精品久久久久久久毛片微露脸| 久久 成人 亚洲| 黄色视频不卡| 色综合婷婷激情| 国产av精品麻豆| 午夜日韩欧美国产| 精品国产乱码久久久久久男人| 高清欧美精品videossex| 91av网站免费观看| 欧美日韩av久久| 男女边摸边吃奶| 亚洲五月色婷婷综合| 一边摸一边抽搐一进一小说 | 黄片小视频在线播放| 制服人妻中文乱码| 人妻一区二区av| 国产成人av激情在线播放| 午夜福利,免费看| 久久久久久久久免费视频了| 国产在线视频一区二区| 中文字幕高清在线视频| 亚洲九九香蕉| 欧美激情久久久久久爽电影 | 亚洲精品中文字幕在线视频| 亚洲一卡2卡3卡4卡5卡精品中文| 99国产极品粉嫩在线观看| 操美女的视频在线观看| 日韩成人在线观看一区二区三区| 日韩中文字幕视频在线看片| 日本欧美视频一区| 国产亚洲午夜精品一区二区久久| 18禁国产床啪视频网站| 中文字幕av电影在线播放| 性高湖久久久久久久久免费观看| 真人做人爱边吃奶动态| 国产成人精品在线电影| 国产在视频线精品| 最近最新中文字幕大全免费视频| 一二三四在线观看免费中文在| 久久天躁狠狠躁夜夜2o2o| 一级a爱视频在线免费观看| 久久久久国内视频| 不卡一级毛片| 欧美日韩视频精品一区| 大型黄色视频在线免费观看| 夜夜夜夜夜久久久久| 丝袜喷水一区| 欧美一级毛片孕妇| 少妇被粗大的猛进出69影院| 国产在线精品亚洲第一网站| 十八禁网站免费在线| 黄色毛片三级朝国网站| 亚洲九九香蕉| 亚洲一区中文字幕在线| 高清毛片免费观看视频网站 | 我的亚洲天堂| 国产欧美亚洲国产| 丝袜美足系列| 欧美日韩亚洲综合一区二区三区_| av不卡在线播放| 久久精品国产亚洲av高清一级| 岛国在线观看网站| 国产黄色免费在线视频| 欧美 日韩 精品 国产| 国产精品秋霞免费鲁丝片| 黑人巨大精品欧美一区二区蜜桃| 亚洲欧美激情在线| 国产精品偷伦视频观看了| 婷婷成人精品国产| 日韩一卡2卡3卡4卡2021年| 在线观看人妻少妇| 国产精品免费大片| 午夜福利,免费看| 国产av国产精品国产| videosex国产| 午夜福利影视在线免费观看| 亚洲黑人精品在线| 妹子高潮喷水视频| 无人区码免费观看不卡 | √禁漫天堂资源中文www| 欧美午夜高清在线| 久久99热这里只频精品6学生| 亚洲精品国产精品久久久不卡| 日日夜夜操网爽| 亚洲av电影在线进入| 一本久久精品| 欧美精品一区二区大全| 女同久久另类99精品国产91| 午夜精品久久久久久毛片777| 久久精品国产a三级三级三级| 一本—道久久a久久精品蜜桃钙片| 美女午夜性视频免费| 男女高潮啪啪啪动态图| 欧美日本中文国产一区发布| 一本久久精品| 久久免费观看电影| 欧美日韩黄片免| 久久久精品国产亚洲av高清涩受| 少妇精品久久久久久久| 欧美精品av麻豆av| 日韩熟女老妇一区二区性免费视频| 无限看片的www在线观看| 人成视频在线观看免费观看| 久久久久久免费高清国产稀缺| 无人区码免费观看不卡 | 天堂8中文在线网| 午夜福利免费观看在线| 岛国在线观看网站| 欧美精品一区二区免费开放| 精品人妻1区二区| 国产男女内射视频| 国产av一区二区精品久久| 怎么达到女性高潮| 91av网站免费观看| 女同久久另类99精品国产91| 黄色片一级片一级黄色片| 丁香欧美五月| 色尼玛亚洲综合影院| 捣出白浆h1v1| 我要看黄色一级片免费的| 国产精品 国内视频| 丝袜在线中文字幕| 国产区一区二久久| 国产精品免费一区二区三区在线 | 五月天丁香电影| 亚洲av国产av综合av卡| 好男人电影高清在线观看| 91麻豆av在线| 高清视频免费观看一区二区| 夜夜骑夜夜射夜夜干| 男女边摸边吃奶| 午夜视频精品福利| 成人国语在线视频| 高清视频免费观看一区二区| 中文字幕精品免费在线观看视频| 在线观看免费日韩欧美大片| 国产亚洲精品第一综合不卡| 久久人妻av系列| 国产成人精品在线电影| 久热爱精品视频在线9| 热99re8久久精品国产| 欧美黑人精品巨大| 精品高清国产在线一区| 男女边摸边吃奶| 韩国精品一区二区三区| 国产片内射在线| 国产一区二区三区在线臀色熟女 | 欧美性长视频在线观看| 青草久久国产| 精品亚洲成a人片在线观看| 伦理电影免费视频| 99精国产麻豆久久婷婷| 黑人巨大精品欧美一区二区mp4| 菩萨蛮人人尽说江南好唐韦庄| 在线永久观看黄色视频| 99九九在线精品视频| 又大又爽又粗| 免费日韩欧美在线观看| 午夜精品久久久久久毛片777| 精品一区二区三卡| 亚洲色图综合在线观看| av网站免费在线观看视频| 成人av一区二区三区在线看| 国产激情久久老熟女| 成人免费观看视频高清| 侵犯人妻中文字幕一二三四区| 国产精品秋霞免费鲁丝片| 热99久久久久精品小说推荐| 亚洲av片天天在线观看| 亚洲人成电影观看| 波多野结衣av一区二区av| av一本久久久久| 国产一区二区三区在线臀色熟女 | 五月天丁香电影| 99香蕉大伊视频| 午夜福利乱码中文字幕| 亚洲av成人一区二区三| 桃红色精品国产亚洲av| 丁香欧美五月| 日日爽夜夜爽网站| 最近最新免费中文字幕在线| av天堂在线播放| 9色porny在线观看| 亚洲国产精品一区二区三区在线| 中文字幕制服av| svipshipincom国产片| 国产精品一区二区在线观看99| 变态另类成人亚洲欧美熟女 | 少妇被粗大的猛进出69影院| 亚洲欧美一区二区三区久久| 国产精品成人在线| 国产成人系列免费观看| 亚洲精品av麻豆狂野| 男女边摸边吃奶| 国产亚洲精品第一综合不卡| 精品一区二区三区av网在线观看 | 丰满少妇做爰视频| 欧美精品一区二区免费开放| 精品国产一区二区三区四区第35| 成人精品一区二区免费| 亚洲人成电影观看| 免费少妇av软件| 丁香六月欧美| 又大又爽又粗| 国产在线观看jvid| 国产免费视频播放在线视频| 18在线观看网站| 色在线成人网| 日日爽夜夜爽网站| 又大又爽又粗| 国产日韩欧美在线精品| 少妇猛男粗大的猛烈进出视频| 中亚洲国语对白在线视频| 99久久人妻综合| 青青草视频在线视频观看| 天天操日日干夜夜撸| 99国产极品粉嫩在线观看| 精品一区二区三区av网在线观看 | av电影中文网址| 久久久国产成人免费| 国产免费av片在线观看野外av| 母亲3免费完整高清在线观看| 99久久人妻综合| 我的亚洲天堂| 亚洲国产欧美一区二区综合| 亚洲av成人不卡在线观看播放网| 男女高潮啪啪啪动态图| 在线天堂中文资源库| 麻豆成人av在线观看| 99久久精品国产亚洲精品| 国产伦理片在线播放av一区| e午夜精品久久久久久久| 国产精品免费一区二区三区在线 |