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

    改進(jìn)蝙蝠算法在優(yōu)化PEMFC運(yùn)行參數(shù)研究中的應(yīng)用

    2024-12-13 00:00:00毛逸君徐洪濤
    太陽(yáng)能學(xué)報(bào) 2024年11期
    關(guān)鍵詞:數(shù)值分析

    摘 要:質(zhì)子交換膜燃料電池工況參數(shù)的優(yōu)化具有多維度和非線性的特點(diǎn),因此該文提出一種改進(jìn)蝙蝠優(yōu)化算法對(duì)其進(jìn)行優(yōu)化。首先,在種群初始化過程中采用Tent混沌映射策略,以提高初始種群遍歷性和多樣性。其次,通過自適應(yīng)權(quán)重策略提高蝙蝠個(gè)體速度動(dòng)態(tài)更新能力,避免算法陷入局部最優(yōu),并使用4種標(biāo)準(zhǔn)測(cè)試函數(shù)驗(yàn)證改進(jìn)算法的優(yōu)越性。然后,通過正交實(shí)驗(yàn)法和熵權(quán)法提出一種新的綜合性能評(píng)價(jià)目標(biāo),并將其作為智能算法優(yōu)化目標(biāo)。最后,將改進(jìn)后算法和綜合性能評(píng)價(jià)目標(biāo)應(yīng)用于質(zhì)子交換膜燃料電池運(yùn)行工況參數(shù)的優(yōu)化模擬。結(jié)果表明,對(duì)于質(zhì)子交換膜燃料電池運(yùn)行參數(shù)模擬優(yōu)化,改進(jìn)蝙蝠優(yōu)化算法比原始蝙蝠優(yōu)化算法的優(yōu)化結(jié)果提高了8.6%,說(shuō)明改進(jìn)蝙蝠算法具有更好的精度。

    關(guān)鍵詞:質(zhì)子交換膜燃料電池;數(shù)值分析;統(tǒng)計(jì)學(xué)方法;蝙蝠算法;自適應(yīng)權(quán)重

    中圖分類號(hào):TM911.4" " " " " " " " " " " " " " "文獻(xiàn)標(biāo)志碼:A

    0 引 言

    “十四五”規(guī)劃提出要進(jìn)一步加強(qiáng)新能源等戰(zhàn)略新興產(chǎn)業(yè)的發(fā)展。其中,質(zhì)子交換膜燃料電池(proton exchange membrane fuel cell, PEMFC)因其可在常溫下啟動(dòng)、效率較高以及零污染排放等優(yōu)點(diǎn),成為未來(lái)綠色能源研究重點(diǎn)[1-3]。由于PEMFC是一個(gè)多物理場(chǎng)耦合的復(fù)雜部件,其運(yùn)行參數(shù)對(duì)性能的影響也具有非線性的特點(diǎn)。不同的運(yùn)行參數(shù)會(huì)影響系統(tǒng)的功率密度與輸出性能,且過于極端的運(yùn)行參數(shù)會(huì)極大縮短PEMFC的使用壽命[4],這對(duì)PEMFC運(yùn)行條件的選擇是一個(gè)嚴(yán)峻考驗(yàn)。因此,在PEMFC的研究中,選擇合適的運(yùn)行參數(shù)對(duì)提高輸出功率和延長(zhǎng)系統(tǒng)壽命至關(guān)重要。

    為研究運(yùn)行參數(shù)對(duì)PEMFC輸出性能的影響,國(guó)內(nèi)外學(xué)者進(jìn)行了大量研究[5-7]。胡冬海等[8]探究了PEMFC與最佳工作溫度的映射關(guān)系,結(jié)果表明,在相同PEMFC輸出電流下,隨著工作溫度升高,PEMFC輸出電壓先升高后降低,存在最佳工作溫度使PEMFC輸出功率最大;Hamrang等[9]比較了平行蛇形流場(chǎng)在不同工況下的性能參數(shù),發(fā)現(xiàn)增加PEMFC的溫度、工作壓力和化學(xué)計(jì)量比都能相應(yīng)提高PEMFC的性能;楊子榮等[10]綜合分析了PEMFC工況參數(shù)對(duì)輸出性能的影響,結(jié)果表明,進(jìn)氣壓強(qiáng)、運(yùn)行溫度和化學(xué)計(jì)量比提高均能提高PEMFC的輸出功率,其中,化學(xué)計(jì)量比過高有可能引發(fā)“膜干”從而降低輸出性能。

    由于工況參數(shù)對(duì)PEMFC系統(tǒng)的影響較為復(fù)雜,單獨(dú)進(jìn)行數(shù)值模擬較難篩選出最優(yōu)工況參數(shù)組合。而智能算法能全局優(yōu)化數(shù)值模擬模型工況參數(shù),具有良好的適用性,被PEMFC模型研究廣泛采用。仇俊政等[11]對(duì)PID控制系統(tǒng)采用粒子群優(yōu)化算法(particle swarm optimization, PSO)進(jìn)行優(yōu)化,獲得了對(duì)運(yùn)行參數(shù)適應(yīng)性更好的PID參數(shù),但PSO本身存在難以跳出局部最優(yōu)與搜索速度不夠穩(wěn)定的缺陷;陳陽(yáng)等[12]利用花授粉優(yōu)化算法(flower pollination algorithm, FPA)對(duì)PSO算法進(jìn)行優(yōu)化,提出雙子群優(yōu)化算法(Bi-subgroup optimization algorithm, BSOA)來(lái)對(duì)PEMFC進(jìn)行優(yōu)化,但由于BSOA是根據(jù)總體適應(yīng)度對(duì)個(gè)體進(jìn)行排序分組,仍易陷入局部最優(yōu);段付德等[13]通過在烏鴉優(yōu)化算法(crow search algorithm, CSA)的基礎(chǔ)上加入柯西變異,來(lái)擴(kuò)大變異范圍并改善全局搜索能力,避免陷入局部最優(yōu),但該算法仍存在CSA缺乏最有個(gè)體引導(dǎo),存在一定盲目性的問題。

    為改善現(xiàn)有優(yōu)化算法在PEMFC模型仿真優(yōu)化研究中的不足,本文提出一種改進(jìn)蝙蝠優(yōu)化算法(improved bat algorithm, IBA)以及其在PEMFC運(yùn)行參數(shù)優(yōu)化研究中的應(yīng)用。該算法以蝙蝠優(yōu)化算法(bat algorithm, BA)為基礎(chǔ),加入帳篷(Tent)混沌映射策略以改善初始種群遍歷性,擴(kuò)大優(yōu)化算法初始搜索范圍。針對(duì)BA收斂速度快但易陷入局部最優(yōu)的問題,本文提出一種自適應(yīng)權(quán)重系數(shù)來(lái)控制蝙蝠搜索范圍,改善優(yōu)化算法全局搜索能力。本文使用IBA對(duì)4組測(cè)試函數(shù)進(jìn)行計(jì)算,結(jié)果表明,IBA不僅具有BA快速收斂的優(yōu)點(diǎn),并避免了其易陷入局部最優(yōu)的痛點(diǎn)。并且,IBA的收斂精度與求解速度均優(yōu)于PSO和BA。同時(shí),為更好地評(píng)價(jià)PEMFC綜合性能,本文通過結(jié)合正交實(shí)驗(yàn)法(orthogonal experimental method, OEM)和熵權(quán)法(entropy weight method, EWM)提出整合電流密度、陰極壓降和氧氣不均勻度的綜合性能評(píng)價(jià)函數(shù)[14],并將其作為智能優(yōu)化算法的適應(yīng)度函數(shù)。最后,本文用IBA對(duì)直流道質(zhì)子交換膜燃料電池工況參數(shù)進(jìn)行優(yōu)化模擬,并與BA優(yōu)化結(jié)果進(jìn)行比較。通過對(duì)算法收斂速度與求解精度進(jìn)行分析,相比于BA,IBA在PEMFC工況參數(shù)優(yōu)化問題上具有更好的優(yōu)化性能。

    1 數(shù)值模型

    1.1 物理模型

    PEMFC主要由流場(chǎng)板(flow field plate,F(xiàn)FP)、質(zhì)子交換膜(proton exchange membrane,PEM)、氣體擴(kuò)散層(gas diffusion layer,GDL)和催化層(catalyst layer,CL)組成,如圖1所示。表1列出了PEMFC模型的幾何參數(shù)與物理參數(shù)。本文主要研究目的是為了模擬智能算法優(yōu)化運(yùn)行參數(shù)對(duì)PEMFC輸出性能的影響,而直流道結(jié)構(gòu)可降低方法驗(yàn)證時(shí)間,有利于后續(xù)將優(yōu)化方法推廣到其他模型結(jié)構(gòu)。

    1.2 數(shù)學(xué)模型

    本文所采用的數(shù)學(xué)模型基于如下假設(shè)[15-16]:

    1)PEMFC在單相和穩(wěn)態(tài)條件下運(yùn)行。

    2)氣體混合物被視為理想氣體,流動(dòng)為不可壓縮流動(dòng)。

    3)多孔介質(zhì)被視為各向同性和均質(zhì)多孔材料,且孔隙率保持不變。

    4)反應(yīng)氣體不可在質(zhì)子交換膜中滲透。

    基于以上假設(shè),模型數(shù)學(xué)描述方程[17-18]為:

    質(zhì)量守恒方程:

    [▽·(ερu)=Sm] (1)

    式中:[ε]——多孔介質(zhì)孔隙率;[ρ]——?dú)怏w混合物的密度,kg/m3;[u]——?dú)怏w的速度矢量,m/s;Sm——質(zhì)量源項(xiàng)。

    動(dòng)量守恒方程:

    [▽·(ερuu)=-ε▽p+▽·(εμ▽u)+Sv] (2)

    式中:[μ]——?jiǎng)恿︷ざ?,Pa·s;[p]——?dú)怏w混合物壓強(qiáng),Pa;[Sv]——?jiǎng)恿吭错?xiàng)。

    組分守恒方程:

    [▽?(εcku)=▽·(Deffk▽ck)+Sk] (3)

    式中:[Deffk]——組分有效擴(kuò)散系數(shù);[ck]——?dú)怏w組分濃度,mol/m3;[Sk]——組分源項(xiàng)。

    能量守恒方程:

    [ρcpu·▽T=▽·(λ▽T)+ST] (4)

    式中:[cp]——比熱容,J/(kg·K);[T]——工作溫度,K;[λ]——導(dǎo)熱系數(shù),W/(m·K);[ST]——溫度源項(xiàng)。

    電荷守恒方程:

    [▽·(keffs▽?s)+S?s=0] (5)

    [▽·(keffm▽?m)+S?m=0] (6)

    式中:[keffm]、[keffs]——膜和固相的電導(dǎo)率,S/m;[?m]、[?s]——膜和固相的電勢(shì),V;[S?m]——離子電流源項(xiàng);[S?s]——電子電流源項(xiàng)。

    Butler-Volmer方程:

    [i=i0expαRFΔVRT-exp-α0FΔVRT] (7)

    式中:[αR]、[α0]——傳遞系數(shù);[i0]——交換電流密度,A/m2;[i]——電流密度,A/m2;[F]——法拉第常數(shù),96481 C/mol;[ΔV]——活化過電勢(shì),V;[R]——通用氣體常數(shù),8.314 J/(mol·K)。

    1.3 邊界條件

    該模型陰陽(yáng)極流道入口設(shè)置為速度入口邊界條件,如式(8)、式(9)所示[19]:

    [uina=ξai02FAm1ωinH2·RTinapina·1Ach] (8)

    [uinca=ξcai02FAm1ωinO2·RTincapinca·1Ach] (9)

    式中:[Ach]、[Am]——流道入口截面積以及活化區(qū)域面積,m2;[ξa]、[ξca]——陽(yáng)極和陰極化學(xué)計(jì)量比;[ωin]——不同組分氣體質(zhì)量分?jǐn)?shù)。

    模型出口設(shè)為壓力出口;GDL和CL兩邊平行于[yz]平面的側(cè)面均采用對(duì)稱邊界條件;陽(yáng)極電勢(shì)和陰極電勢(shì)分別設(shè)置為0與0.4 V。

    1.4 氧氣不均勻度

    PEMFC陰極側(cè)的氧氣分布會(huì)影響電化學(xué)反應(yīng)速度與效率,因此,氧氣不均勻度也被作為PEMFC性能評(píng)價(jià)指標(biāo)之一。為了定量分析GDL/CL界面上氧氣分布情況,將氧氣不均勻度定義為[20]:

    [S=(a-aave)2dAmaave2dAm] (10)

    [aave=adAmdAm] (11)

    式中:[aave]——GDL/CL界面上氧氣的平均摩爾濃度,mol/m2;[a]——GDL/CL界面上氧氣的局部摩爾濃度,mol/m2。

    1.5 網(wǎng)格無(wú)關(guān)性驗(yàn)證與模型驗(yàn)證

    本文采用13440、21500、43282和61372這4種網(wǎng)格數(shù)量對(duì)工作溫度為180 ℃、陰陽(yáng)極進(jìn)氣濕度80%、陰陽(yáng)極進(jìn)氣化學(xué)計(jì)量比分別為2.0和1.5的PEMFC直流道模型進(jìn)行網(wǎng)格無(wú)關(guān)性驗(yàn)證。網(wǎng)格數(shù)量對(duì)0.4 V下電流密度的影響如圖2所示。仿真結(jié)果表明,43282和61372網(wǎng)格數(shù)量之間的電流密度相對(duì)變化為0.02%。為平衡計(jì)算時(shí)間和計(jì)算精度,本研究采用43282網(wǎng)格數(shù)量進(jìn)行仿真。

    為確保本文模型的準(zhǔn)確性,在相同邊界條件和物理模型下,將本研究模擬獲得的極化曲線與Sezgin等[17]模型的極化曲線進(jìn)行比較。如圖3所示,兩者吻合良好。

    2 蝙蝠優(yōu)化算法

    蝙蝠算法是一種基于迭代優(yōu)化方法的群智能優(yōu)化算法,它是由理想化蝙蝠回聲定位捕食行為提出的。該算法由楊新社等[21]于2010年提出,具有模型簡(jiǎn)單和收斂速度快等優(yōu)點(diǎn),其遵循如下近似規(guī)則[22]:

    1)種群內(nèi)所有蝙蝠都使用回聲定位來(lái)感知距離,并有特定方法分辨障礙物與目標(biāo)區(qū)別;

    2)蝙蝠以速度[vi]在位置[xi]隨機(jī)飛行,頻率[f]固定,并在搜索獵物時(shí)根據(jù)目標(biāo)與自己之間的距離自動(dòng)調(diào)整波長(zhǎng)[Λ]與響度[A];

    3)響度從最大值[A0]變化到固定最小值[Amin]。

    假設(shè)該算法中蝙蝠種群的搜索空間為[D]維,則蝙蝠[t]時(shí)刻的速度[vi]和位置[xi]更新公式為:

    [fi=fmin+(fmax-fmin)·β] (12)

    [vti=vt-1i+(xti-x*)·fi] (13)

    [xti=xt-1i+vti] (14)

    式中:[fmax]——蝙蝠發(fā)出聲波的最大頻率;[fmin]——蝙蝠發(fā)出聲波的最小頻率;[x*]——當(dāng)前全局最優(yōu)位置;[β]——區(qū)間[0,1]內(nèi)的一個(gè)隨機(jī)向量。

    當(dāng)蝙蝠出現(xiàn)最優(yōu)位置時(shí),原最優(yōu)位置[xold]附近會(huì)根據(jù)式(15)隨機(jī)游走產(chǎn)生新的位置xnew。

    [xnew=xold+δAt] (15)

    式中:[δ]——區(qū)間[[-1],1]內(nèi)一個(gè)隨機(jī)數(shù);[At]——所有蝙蝠在這個(gè)時(shí)間步驟[t]內(nèi)的響度[Ati]的平均值。

    蝙蝠在最初捕食期間會(huì)以較低的脈沖發(fā)射率和較高的脈沖響度來(lái)搜尋獵物,以期獲得更大的搜索范圍;接近獵物后,蝙蝠會(huì)通過降低脈沖響度并提高脈沖發(fā)射率來(lái)保證搜索精度,獲得更高的捕食效率。若第[i]只蝙蝠在[t]時(shí)刻隨機(jī)游走后,蝙蝠的最優(yōu)位置發(fā)生改變,且該蝙蝠當(dāng)次搜索響度大于隨機(jī)響度,則將該蝙蝠位置設(shè)為當(dāng)前最優(yōu)位置[Xnew],并按照式(16)、式(17)更新其脈沖發(fā)射率以及脈沖響度;否則蝙蝠的最優(yōu)位置、脈沖發(fā)射率及脈沖響度不變。

    [rt+1i=r0i[1-exp(-γt)]] (16)

    [At+1i=αsAti] (17)

    式中:[rt+1i]——[t+1]時(shí)刻第[i]只蝙蝠的脈沖發(fā)射率;[At+1i]——[t+1]時(shí)刻第[i]只蝙蝠的脈沖響度;[γ]——脈沖速率增強(qiáng)系數(shù),[γ∈(0,+∞)];[αs]——聲波響度衰減系數(shù),[αs∈(0,1)]。

    當(dāng)?shù)螖?shù)等于最大迭代次數(shù)[G]時(shí),或脈沖響度[At+1i]趨向于0、脈沖發(fā)射率[rt+1i]趨向于[r0i]時(shí),蝙蝠找到獵物并達(dá)到最優(yōu)位置。

    3 改進(jìn)蝙蝠優(yōu)化算法

    PEMFC的運(yùn)行參數(shù)之間具有復(fù)雜的非線性關(guān)系,傳統(tǒng)的蝙蝠算法雖具有較快的收斂速度來(lái)搜索最優(yōu)位置,但由于蝙蝠個(gè)體位置優(yōu)化和隨機(jī)游走以及速度的更新方式都受到隨機(jī)值影響,因此在搜索過程中定位精度缺乏準(zhǔn)確性且易陷入局部最優(yōu)[23],這會(huì)嚴(yán)重影響優(yōu)化結(jié)果。針對(duì)上述問題,本文提出一種基于Tent混沌映射初始化的自適應(yīng)蝙蝠算法。通過Tent映射策略提高初始種群均勻性和遍歷性,并通過自適應(yīng)權(quán)重策略均衡算法局部搜索與全局搜索能力,填補(bǔ)蝙蝠算法容易陷入局部最優(yōu)的缺陷,提升搜索精度。

    3.1 Tent混沌映射初始化

    算法初始搜索能力與蝙蝠初始種群的遍歷性和多樣性密切相關(guān),初始種群位置多樣富且豐的種群可大幅提高蝙蝠算法的定位精度與收斂速度[24]。原始蝙蝠算法依賴于隨機(jī)生成對(duì)種群進(jìn)行初始化,會(huì)導(dǎo)致種群初始化后分布不均、多樣性較差,增加算法收斂時(shí)長(zhǎng),降低搜索精度。而多數(shù)群體智能算法為增加算法的隨機(jī)性和多樣性,多在種群初始化中采用混沌映射來(lái)生成混沌序列。因此,為提高初始種群的遍歷性與均勻性,本文采用具有均勻分布函數(shù)且相關(guān)性良好的Tent混沌映射來(lái)初始化種群,維持種群多樣性,提高算法全局搜索能力[25]。Tent映射公式為:

    [Xn+1=107Xn," "0≤Xn≤0.7103(1-Xn)," 0.7lt;Xn≤1] (18)

    式中:[Xn]——區(qū)間[0,1]內(nèi)隨機(jī)產(chǎn)生的隨機(jī)產(chǎn)生的[D]維向量的第[n]維分量;[Xn+1]——經(jīng)Tent映射后的變量值。

    根據(jù)Tent映射,初始化蝙蝠種群位置為:

    [xi=xmin+Xn+1·(xmax-xmin)] (19)

    式中:[xmax]、[xmin]——搜索范圍的上下邊界。

    3.2 自適應(yīng)權(quán)重速度更新策略

    初始化種群后,種群內(nèi)蝙蝠開始對(duì)目標(biāo)進(jìn)行定位。原始蝙蝠算法中,蝙蝠位置更新主要根據(jù)當(dāng)前最優(yōu)位置,搜索方式過于單一,易陷入局部最優(yōu)。為平衡該算法的局部尋優(yōu)和全局搜索能力,本文引入自適應(yīng)權(quán)重速度更新策略,優(yōu)化后速度更新公式如式(20)及式(21)所示[26]:

    [vti=ωtivt-1i+(xti-x*)?fi] (20)

    [ωti=ωmax-(ωmax-ωmin)Fti-FtminFtave-Ftmin," Fti≤Ftaveωmax," Ftigt;Ftave] (21)

    式中:[ωmax]、[ωmin]——慣性權(quán)重的上下邊界;[ωti]——[t]時(shí)刻第[i]只蝙蝠的慣性權(quán)重;[Ftave]——[t]時(shí)刻種群內(nèi)所有蝙蝠平均適應(yīng)度;[Ftmin]——[t]時(shí)刻種群內(nèi)所有蝙蝠最優(yōu)適應(yīng)度;[Fti]——[t]時(shí)刻第[i]只蝙蝠的適應(yīng)度值。

    該慣性權(quán)重會(huì)在優(yōu)化初期為整個(gè)種群賦予較大慣性權(quán)重,并在優(yōu)化后期賦予較小慣性權(quán)重,以此來(lái)保證初期全局搜索能力和算法收斂性能。

    3.3 改進(jìn)蝙蝠優(yōu)化算法執(zhí)行步驟

    改進(jìn)蝙蝠算法執(zhí)行步驟為:

    第1步:設(shè)置初始參數(shù)。蝙蝠種群數(shù)量[n],最大迭代次數(shù)[G],初始脈沖響度[A0],初始脈沖發(fā)射率[r0],脈沖增強(qiáng)系數(shù)[γ],脈沖頻率范圍[[fmin,fmax]],聲波衰減系數(shù)[αs],慣性權(quán)重范圍[[ωmin,ωmax]]。

    第2步:使用Tent混沌映射初始化種群,生成[n]個(gè)可行解作為初始蝙蝠種群個(gè)體位置[xi],并計(jì)算對(duì)應(yīng)適應(yīng)度[F(xi)]。

    第3步:篩選出當(dāng)前最優(yōu)位置[x*],按照式(12)~式(14)和式(20)、式(21)自適應(yīng)更新蝙蝠速度與蝙蝠位置。

    第4步:生成一個(gè)隨機(jī)數(shù)[δ1∈[0,1]],若[δ1gt;ri],則在當(dāng)前最優(yōu)蝙蝠個(gè)體位置[x*]附近通過式(15)增加擾動(dòng)更新位置;否則,通過式(14)更新蝙蝠位置。計(jì)算新位置對(duì)應(yīng)的適應(yīng)度。

    第5步:生成另一個(gè)隨機(jī)數(shù)[δ2∈[0,1]],若最優(yōu)位置得到更新且[δ2lt;Ai],則接受該位置,并根據(jù)式(16)、式(17)改變蝙蝠的脈沖響度以及脈沖發(fā)射率;否則,保持不變。

    第6步:當(dāng)[tlt;G]時(shí),重復(fù)第3步~第5步;否則,執(zhí)行第7步。

    第7步:當(dāng)達(dá)到最大迭代次數(shù)時(shí),輸出優(yōu)化結(jié)果。

    改進(jìn)后蝙蝠算法流程見圖4。

    4 算法性能校驗(yàn)

    為測(cè)試IBA的性能,選擇單峰函數(shù)Sphere和Schwefel、多峰函數(shù)Rastrigin和Griewank 4個(gè)標(biāo)準(zhǔn)測(cè)試函數(shù)對(duì)其進(jìn)行測(cè)試[27],測(cè)試函數(shù)如表2所示。算法測(cè)試結(jié)果越接近表中理論最優(yōu)值則說(shuō)明該算法性能越優(yōu),準(zhǔn)確性越高。并將測(cè)試結(jié)果與BA、PSO進(jìn)行比較,測(cè)試結(jié)果如表3所示。IBA和BA參數(shù)設(shè)置為:蝙蝠種群數(shù)量[n=20],初始脈沖響度[A0=1],脈沖頻率范圍[[fmin,fmax]]=[0,2],初始脈沖發(fā)射率[r0=0.01],脈沖增強(qiáng)系數(shù)[γ=0.9],聲波衰減系數(shù)[αs=0.9]。其中,IBA自適應(yīng)慣性權(quán)重范圍[[ωmin,ωmax]]=[0.4,1]。PSO參數(shù)設(shè)置為:粒子總數(shù)[n=20],慣性權(quán)重[ω=0.9],學(xué)習(xí)因子分別為[c1=1.5、c2=2]。為使結(jié)果更加完善,故將最大迭代次數(shù)均設(shè)置為500。同時(shí),對(duì)每個(gè)測(cè)試函數(shù)都使用各算法獨(dú)立運(yùn)行30次,以此來(lái)保證驗(yàn)證結(jié)果的客觀性。在維度[D=30]的情況下,計(jì)算30次運(yùn)行后目標(biāo)函數(shù)結(jié)果的平均值、最優(yōu)值和標(biāo)準(zhǔn)差。

    由表3可看出,對(duì)于多峰函數(shù)Rastrigin和Griewank以及單峰函數(shù)Sphere和Schwefle,IBA最優(yōu)值和平均值相比于BA和PSO算法都更接近于理論最優(yōu)值,說(shuō)明IBA算法經(jīng)過優(yōu)化后,準(zhǔn)確性高于PSO和BA算法,特別是對(duì)于有許多規(guī)則分布的局部最小值的Griewank函數(shù),IBA算法測(cè)試結(jié)果更加準(zhǔn)確。且IBA算法運(yùn)行30次后所得標(biāo)準(zhǔn)差均小于PSO和BA算法,說(shuō)明IBA算法穩(wěn)定性更為良好,算法綜合性能優(yōu)于PSO與BA算法,尤其是Sphere函數(shù),IBA標(biāo)準(zhǔn)差相對(duì)于BA算法提高11個(gè)數(shù)量級(jí),相對(duì)于PSO算法提高8個(gè)數(shù)量級(jí)。

    5 結(jié)果與討論

    5.1 綜合性能評(píng)價(jià)函數(shù)

    PEMFC輸出性能主要與電流密度與陰極壓降有關(guān),但陰極氧氣不均勻度也會(huì)影響反應(yīng)速率。因此,為了更好地評(píng)價(jià)PEMFC性能,本文通過正交實(shí)驗(yàn)表法與熵權(quán)法提出一種新的綜合性能評(píng)價(jià)函數(shù),并以此作為優(yōu)化算法適應(yīng)度函數(shù)。本模型運(yùn)行工況參數(shù)優(yōu)化范圍如表4所示。為盡可能全面覆蓋優(yōu)化范圍,綜合考慮后本研究選擇5研究因素4因素水平,并選擇L16(45)正交表設(shè)計(jì)仿真實(shí)驗(yàn)方案。同時(shí),本文選擇電流密度(I)、陰極壓降(Δp)和陰極氧氣不均勻度(S)作為性能研究目標(biāo),并對(duì)16組方案進(jìn)行仿真模擬,結(jié)果如表5所示。

    熵權(quán)法是一種定量加權(quán)方法,它根據(jù)目標(biāo)變異性的大小來(lái)分配客觀權(quán)重。其中,目標(biāo)的信息熵越小,即該目標(biāo)的可變性越大,提供的信息就越多,在整體評(píng)估中能發(fā)揮更大作用,因此,其相應(yīng)權(quán)重也越大。作為一種客觀賦權(quán)法,熵權(quán)法避免了人為因素導(dǎo)致的偏差,相對(duì)于主觀賦權(quán)法有更高的精度和客觀性。本文通過熵權(quán)法將電流密度、陰極壓降和陰極氧氣不均勻度性能目標(biāo)整合為一種綜合性能評(píng)價(jià)目標(biāo),熵權(quán)法數(shù)學(xué)模型如式(22)~式(28)所示[28-29]:

    1)根據(jù)研究目標(biāo)數(shù)值建立初始矩陣[A]。

    [A=a11a12…a1na21a22…a2n????am1am2…amn] (22)

    式中:[amn]——研究目標(biāo)的值;[m]——方案序號(hào);[n]——研究目標(biāo)數(shù)量。

    2)對(duì)各研究目標(biāo)的值去量綱化處理。

    [B=b11b12…b1nb21b22…b2n????bm1bm2…bmn] (23)

    [bij=aij-mini=1,2,…maijmaxi=1,2,…,maij-mini=1,2,…,maij," 正向目標(biāo)bij=maxi=1,2,…,maij-aijmaxi=1,2,…,maij-mini=1,2,…,maij," 負(fù)向目標(biāo)] (24)

    式中:[bmn]——標(biāo)準(zhǔn)化后研究目標(biāo)的值。

    其中,正向目標(biāo)是指隨著目標(biāo)數(shù)值的增加,評(píng)價(jià)對(duì)象的綜合評(píng)價(jià)值也隨著增加的目標(biāo);負(fù)向目標(biāo)則相反。本文中,電流密度為正向目標(biāo),陰極壓降和陰極氧氣不均勻度為負(fù)向目標(biāo)。

    3)計(jì)算各目標(biāo)在各方案下的比值[pij]。

    [pij=biji=1mbij] (25)

    4)計(jì)算各目標(biāo)信息熵[Ej]。

    [Ej=-1lnmi=1mpijlnpij] (26)

    5)確定各目標(biāo)權(quán)重[Wj]。

    [Wj=1-Ejj=1n(1-Ej)] (27)

    6)計(jì)算每個(gè)方案綜合得分[Zi]。

    [Zi=j=1nWjbij] (28)

    基于熵權(quán)法數(shù)學(xué)模型及上述模擬結(jié)果,計(jì)算得到電流密度、陰極壓降和陰極氧氣不均勻度的權(quán)重分別為0.44、0.30和0.26。因此,綜合性能評(píng)價(jià)目標(biāo)函數(shù)(comprehensive performance evaluation objective function, Ω)表達(dá)式為[Ω=0.44Inon+0.3Δpnon+0.26Snon],如圖5所示。式中,[Inon]、[Δpnon]和[Snon]分別為去量鋼化后的電流密度、陰極壓降和陰極氧氣不均勻度。

    5.2 算法收斂情況對(duì)比

    本文采用COMSOL和Matlab耦合的模擬方式,并根據(jù)綜合性能評(píng)價(jià)目標(biāo)對(duì)直流道PEMFC進(jìn)行建模優(yōu)化,如圖6所示。具體優(yōu)化流程為:利用Tent混沌算法初始化蝙蝠種群后,將種群個(gè)體位置賦值給工作溫度([T])、陽(yáng)極氣體相對(duì)濕度([φa])、陰極氣體相對(duì)濕度([φc])、陰極化學(xué)計(jì)量比([λc])和陽(yáng)極化學(xué)計(jì)量比([λa]),然后更新COMSOL模型并計(jì)算新解,將模擬結(jié)果輸出至IBA算法中計(jì)算適應(yīng)度函數(shù)的值,再繼續(xù)運(yùn)行IBA算法對(duì)應(yīng)步驟,并迭代輸出最優(yōu)解。

    為全面評(píng)價(jià)PEMFC性能,本文將綜合性能評(píng)價(jià)函數(shù)作為算法優(yōu)化適應(yīng)度函數(shù),同時(shí),將運(yùn)行參數(shù)尋優(yōu)范圍約束為表4。當(dāng)適應(yīng)度函數(shù)取到最大值時(shí),PEMFC綜合性能最優(yōu),對(duì)應(yīng)工況參數(shù)組合最好。

    為驗(yàn)證IBA對(duì)于PEMFC工況參數(shù)優(yōu)化的優(yōu)越性,將算法優(yōu)化結(jié)果與BA進(jìn)行比較。為平衡優(yōu)化計(jì)算消耗時(shí)間與精度,降低參數(shù)設(shè)置影響,故設(shè)置蝙蝠種群數(shù)量[n=20],脈沖頻率范圍[[fmin,fmax]]=[0,2],初始脈沖響度[A0=1],初始脈沖發(fā)射率[r0=0.01],最大迭代次數(shù)[G=50],脈沖增強(qiáng)系數(shù)[γ=0.9],聲波衰減系數(shù)[αs=0.9],慣性權(quán)重范圍[[ωmin,ωmax]]=[0.4,0.9],工作電壓為0.4 V。IBA與BA得到的直流道PEMFC工況參數(shù)優(yōu)化結(jié)果如表6所示,兩種算法的迭代收斂曲線如圖7所示。

    根據(jù)表6,IBA和BA兩種算法優(yōu)化結(jié)果對(duì)應(yīng)工況參數(shù)下適應(yīng)度函數(shù)值為0.990和0.912,IBA的適應(yīng)度值與BA相比提高了8.6%。由此可知,IBA算法優(yōu)化精度高于BA算法。由圖7可知,IBA算法雖收斂代數(shù)大于BA算法,但I(xiàn)BA算法得到結(jié)果更貼近最優(yōu)值,說(shuō)明IBA算法準(zhǔn)確性明顯高于BA算法。

    同時(shí),從IBA算法優(yōu)化過程中選擇15組不同數(shù)據(jù)繪制平行坐標(biāo)軸圖,如圖8所示。由圖8a可知,當(dāng)陰陽(yáng)極化學(xué)計(jì)量比與陰陽(yáng)極相對(duì)濕度近似一致時(shí),提高運(yùn)行溫度會(huì)強(qiáng)化氣體傳輸,提高反應(yīng)速率,從而降低濃差極化損失,提高輸出電流密度,使Ω的值增大。由圖8b可知,在近似運(yùn)行溫度和陰陽(yáng)極相對(duì)濕度下,增加陰極化學(xué)計(jì)量比會(huì)增大陰極氣體流速,提高陰極壓降并促使氧氣向GDL/CL界面擴(kuò)散,改善陰極氧氣不均勻度,影響Ω值。且陽(yáng)極化學(xué)計(jì)量比過高時(shí)可能導(dǎo)致“膜干”現(xiàn)象,導(dǎo)致輸出性能下降,使Ω值減小。由圖8c可知,當(dāng)運(yùn)行溫度、陰陽(yáng)極化學(xué)計(jì)量比近似一致時(shí),陰極相對(duì)濕度增加會(huì)使陰極入口處氣體水的摩爾分?jǐn)?shù)增加,陰極氣體流速明顯增大,陰極壓降增加,從而使Ω值減小。而陽(yáng)極相對(duì)濕度增加有利于提高電化學(xué)反應(yīng)中質(zhì)子傳遞效率,提高反應(yīng)速率,增大輸出電流密度,使Ω值增大。但陰陽(yáng)極相對(duì)濕度過高易產(chǎn)生“水淹”現(xiàn)象,反而降低燃料電池輸出功率,降低Ω值。

    6 結(jié) 論

    本文通過研究改進(jìn)蝙蝠優(yōu)化算法在質(zhì)子交換膜燃料電池運(yùn)行參數(shù)中的應(yīng)用,提出一種新的直流道PEMFC工況參數(shù)優(yōu)化方法,并提出一種基于電流密度、陰極壓降與陰極氧氣不均勻度的PEMFC綜合性能評(píng)價(jià)目標(biāo)。主要結(jié)論如下:

    1)改進(jìn)蝙蝠優(yōu)化算法是基于標(biāo)準(zhǔn)蝙蝠優(yōu)化算法,引入Tent混沌映射策略來(lái)提高種群初始化的遍歷性和多樣性;同時(shí)引入自適應(yīng)權(quán)重,根據(jù)實(shí)際情況動(dòng)態(tài)調(diào)整蝙蝠種群速度,避免優(yōu)化算法陷入局部最優(yōu)。算法函數(shù)測(cè)試結(jié)果表明,IBA在求解精度上有明顯提升,能夠更好地平衡局部搜索和全局搜索能力。

    2)通過熵權(quán)法得到質(zhì)子交換膜燃料電池運(yùn)行參數(shù)綜合性能評(píng)價(jià)目標(biāo)Ω,并利用改進(jìn)蝙蝠算法對(duì)其進(jìn)行研究,得到在規(guī)定參數(shù)范圍內(nèi)最優(yōu)運(yùn)行參數(shù)組合為工作溫度453 K、陽(yáng)極相對(duì)濕度80.5%、陰極相對(duì)濕度13.3%、陽(yáng)極化學(xué)計(jì)量比2.40和陰極化學(xué)計(jì)量比3.00,對(duì)應(yīng)綜合性能評(píng)價(jià)目標(biāo)的值為0.990。

    3)通過分析模擬結(jié)果,可得到運(yùn)行溫度越高,輸出電流密度越大,Ω值越大。陽(yáng)極相對(duì)濕度和陰極化學(xué)計(jì)量比的增加能提高PEMFC輸出性能,但過高的相對(duì)濕度會(huì)導(dǎo)致“水淹”現(xiàn)象,同時(shí),過高的陰極相對(duì)濕度也會(huì)增加陰極壓降,過高的陽(yáng)極化學(xué)計(jì)量比可能會(huì)產(chǎn)生“膜干”現(xiàn)象,降低Ω值。進(jìn)一步印證了PEMFC運(yùn)行工況參數(shù)多目標(biāo)優(yōu)化的必要性,以及合理工況設(shè)置對(duì)提高PEMFC輸出性能的重要性。

    [參考文獻(xiàn)]

    [1] PRAMUANJAROENKIJ A, KAKA? S. The fuel cell electric" vehicles:" the" highlight" review[J]." International journal of hydrogen energy, 2023, 48(25): 9401-9425.

    [2] ZHANG S Y, XU H T, QU Z G, et al. Bio-inspired flow channel designs for proton exchange membrane fuel cells: a review[J]." " Journal" "of" "power" "sources," 2022," 522: 231003.

    [3] LI N, WANG W T, XU R Y, et al. Design of a novel nautilus bionic flow field for proton exchange membrane fuel cell by analyzing performance[J]. International journal of heat and mass transfer, 2023, 200: 123517.

    [4] 羅熙. 基于電參數(shù)的燃料電池壽命測(cè)試方法研究[D]. 成都: 電子科技大學(xué), 2019.

    LUO X. Research on life test method of fuel cell based on electrical parameters[D]. Chengdu: University of Electronic Science and Technology of China, 2019.

    [5] 劉玉峰, 施魯浩, 楊來(lái)順, 等. 變工況對(duì)帶渦流發(fā)生器的PEMFC性能的影響[J]. 熱科學(xué)與技術(shù), 2022, 21(4): 347-355.

    LIU Y F, SHI L H, YANG L S, et al. Influence of variable operating conditions on performance of PEMFC with vortex generator[J]. Journal of thermal science and technology, 2022, 21(4): 347-355.

    [6] 趙鑫, 陳光, 張妍懿. 運(yùn)行工況對(duì)PEMFC性能與水含量的影響分析[J]. 汽車工程, 2022, 44(3): 379-384.

    ZHAO X, CHEN G, ZHANG Y Y. An analysis on the effects of operating conditions on the performance and water content of PEMFC[J]. Automotive engineering, 2022, 44(3): 379-384.

    [7] 朱京宇, 談金祝, 孫澳. 陰極相對(duì)濕度對(duì)質(zhì)子交換膜燃料電池電化學(xué)性能的影響[J]. 南京工業(yè)大學(xué)學(xué)報(bào)(自然科學(xué)版), 2021, 43(4): 456-460.

    ZHU J Y, TAN J Z, SUN A. Effects of cathode relative humidity on the electrochemical performance of proton exchange membrane fuel cell[J]. Journal of Nanjing Tech university (natural science edition), 2021, 43(4): 456-460.

    [8] HU D H, WANG Y T, LI J W, et al. Investigation of optimal operating temperature for the PEMFC and its tracking control for energy saving in vehicle applications[J]. Energy conversion and management, 2021, 249: 114842.

    [9] HAMRANG A, ABDOLLAHZADEH M, BILONDI A M, et al. Comparison of PEMFC performance with parallel serpentine and parallel serpentine-baffled flow fields under various operating and geometrical conditions; a parametric study[J]. International journal of hydrogen energy, 2023, 48(20): 7442-7459.

    [10] 楊子榮, 李巖, 冀雪峰, 等. 質(zhì)子交換膜燃料電池運(yùn)行工況參數(shù)敏感性分析[J]. 吉林大學(xué)學(xué)報(bào)(工學(xué)版), 2022, 52(9): 1971-1981.

    YANG Z R, LI Y, JI X F, et al. Sensitivity analysis of operating parameters for proton exchange membrane fuel cells[J]." Journal" "of" "Jilin" University" "(engineering and technology edition), 2022, 52(9): 1971-1981.

    [11] 仇俊政, 趙紅, 牟亮, 等. 基于粒子群PID的質(zhì)子交換膜燃料電池溫度控制[J]. 制造業(yè)自動(dòng)化, 2022, 44(8): 98-101.

    QIU J Z, ZHAO H, MOU L, et al. Temperature control of proton exchange membrane fuel cell based on particle swarm optimization PID[J]. Manufacturing automation, 2022, 44(8): 98-101.

    [12] CHEN Y, PI D C, WANG B, et al. Bi-subgroup optimization algorithm for parameter estimation of a PEMFC" model[J]." Expert" systems" with" applications, 2022, 196: 116646.

    [13] DUAN F D, SONG F, CHEN S N, et al. Model parameters identification of the PEMFCs using an improved design of crow search algorithm[J]. International journal of hydrogen energy, 2022, 47(79): 33839-33849.

    [14] ZHANG S Y, MAO Y J, LIU F, et al. Multi-objective optimization and evaluation of PEMFC performance based on orthogonal experiment and entropy weight method[J]. Energy conversion and management, 2023, 291: 117310.

    [15] DONG P C, XIE G N, NI M. The mass transfer characteristics and energy improvement with various partially blocked flow channels in a PEM fuel cell[J]. Energy, 2020, 206: 117977.

    [16] ZHANG S Y, QU Z G, XU H T, et al. A numerical study on the performance of PEMFC with wedge-shaped fins in the cathode channel[J]. International journal of hydrogen energy, 2021, 46(54): 27700-27708.

    [17] SEZGIN B, CAGLAYAN D G, DEVRIM Y, et al. Modeling and sensitivity analysis of high temperature PEM fuel cells by using Comsol Multiphysics[J]. International journal of hydrogen energy, 2016, 41(23): 10001-10009.

    [18] 梅楠. 梯形流道質(zhì)子交換膜燃料電池氣體流道及冷卻流道研究[D]. 鎮(zhèn)江: 江蘇大學(xué), 2022.

    MEI N. Study on gas channel and cooling channel of trapezoidal channel proton exchange membrane fuel cell[D]. Zhenjiang: Jiangsu University, 2022.

    [19] 王珂, 張拴羊, 徐洪濤, 等. 基于神經(jīng)元的PEMFC仿生流道性能模擬研究[J]. 太陽(yáng)能學(xué)報(bào), 2022, 43(6): 454-459.

    WANG K, ZHANG S Y, XU H T, et al. Simulation study of bionic channel-inspired of proton exchange membrane fuel cell based on neuron[J]. Acta energiae solaris sinica, 2022, 43(6): 454-459.

    [20] LIU H C, YANG W M, TAN J, et al. Numerical analysis of parallel flow fields improved by micro-distributor in proton exchange membrane fuel cells[J]. Energy conversion and management, 2018, 176: 99-109.

    [21] YANG X S, HOSSEIN GANDOMI A. Bat algorithm: a novel approach for global engineering optimization[J]. Engineering computations, 2012, 29(5): 464-483.

    [22] 黎成. 新型元啟發(fā)式蝙蝠算法[J]. 電腦知識(shí)與技術(shù), 2010, 6(23): 6569-6572.

    LI" "C." A" new" metaheuristic" "bat-inspired" "algorithm[J]. Computer knowledge and technology, 2010, 6(23): 6569-6572.

    [23] 謝良波, 李宇洋, 王勇, 等. 基于自適應(yīng)蝙蝠算法的室內(nèi)RFID定位算法[J]. 通信學(xué)報(bào), 2022, 43(8): 90-99.

    XIE L B, LI Y Y, WANG Y, et al. Indoor RFID localization algorithm based on adaptive bat algorithm[J]. Journal on communications, 2022, 43(8): 90-99.

    [24] 楊宇倫, 凌銘. 基于改進(jìn)雞群優(yōu)化算法的質(zhì)子交換膜燃料電池模型參數(shù)辨識(shí)[J]. 太陽(yáng)能學(xué)報(bào), 2023, 44(2): 269-278.

    YANG Y L, LING M. Parameter identification of proton exchange membrane fuel cells model based on improved chicken swarm optimization algorithm[J]. Acta energiae solaris sinica, 2023, 44(2): 269-278.

    [25] VALLE J, MACHICAO J, BRUNO O M. Chaotical PRNG based on composition of logistic and tent maps using deep-zoom[J]. Chaos, solitons amp; fractals, 2022, 161: 112296.

    [26] 張佳丹, 顧圣平, 鄭斯水, 等. 基于改進(jìn)蝙蝠算法的梯級(jí)水庫(kù)發(fā)電優(yōu)化調(diào)度[J]. 人民黃河, 2020, 42(6): 53-57.

    ZHANG J D, GU S P, ZHENG S S, et al. Optimal operation of cascade hydropower station reservoirs based on the improved bat algorithm[J]. Yellow river, 2020, 42(6): 53-57.

    [27] 王廷元, 何先波, 賀春林. 一種基于自適應(yīng)策略的混合鯨魚優(yōu)化算法[J]. 西華師范大學(xué)學(xué)報(bào)(自然科學(xué)版), 2021, 42(1): 92-99.

    WANG T Y, HE X B, HE C L. A hybrid whale optimization algorithm based on adaptive strategy[J]. Journal of China West Normal University (natural sciences), 2021, 42(1): 92-99.

    [28] KUMAR R, SINGH S, BILGA P S, et al. Revealing the benefits of entropy weights method for multi-objective optimization in machining operations: a critical review[J]. Journal of materials research and technology, 2021, 10: 1471-1492.

    [29] WANG Z, XING X G, YAN F. An abnormal phenomenon in entropy weight method in the dynamic evaluation of water quality index[J]. Ecological indicators, 2021, 131: 108137.

    APPLICATION OF IMPROVED BAT ALGORITHM IN OPTIMIZATION OF PEMFC OPERATING PARAMETERS

    Mao Yijun,Xu Hongtao

    (School of Energy and Power Engineering, University of Shanghai for Science and Technology, Shanghai 200093, China)

    Abstract:For the multi-dimensional and non-linear characteristics of the operating parameters optimization of the proton exchange membrane fuel cell (PEMFC) with the straight channel, this paper proposes an improved bat algorithm (IBA) to optimize it. Firstly, this paper initializes the population with tent map strategy to improve the diversity and ergodicity of the initial population. Secondly, the adaptive inertia weight is used to enhance the dynamic updating ability of individual bat speeds to prevent the algorithm from falling into the local optimum. Four sets of standard test functions verify the superiority of IBA. Then a new comprehensive performance evaluation objective is proposed by the orthogonal experiment method (OEM) and the entropy weight method (EWM), and it is used as the optimization objective of the algorithm. Finally, IBA is applied to simulate and optimize the operating parameters of the PEMFC with the straight channel. The simulation results show that the optimization result of IBA is improved by 8.6% compared with the original bat optimization algorithm, indicating that the IBA has higher solution accuracy for the simulation optimization of operating parameters of the PEMFC.

    Keywords:proton exchange membrane fuel cell; numerical analysis; statistical methods; bat algorithm; adaptive inertia weight

    猜你喜歡
    數(shù)值分析
    軟基上碗扣式滿堂支架數(shù)值分析與基礎(chǔ)驗(yàn)算
    軟基上碗扣式滿堂支架數(shù)值分析與基礎(chǔ)驗(yàn)算
    壓力溶腔對(duì)巖溶隧道施工安全影響的數(shù)值分析
    土與支護(hù)結(jié)構(gòu)相互作用及邊坡穩(wěn)定性分析
    探討補(bǔ)償回彈沖壓件模具設(shè)計(jì)的方法
    基于問題式學(xué)習(xí)的《數(shù)值分析》微課設(shè)計(jì)
    已建廠房室內(nèi)沉井施工對(duì)周邊環(huán)境影響的分析
    基于創(chuàng)新和應(yīng)用能力的數(shù)值分析課程教學(xué)研究與實(shí)踐
    民用飛機(jī)靜壓孔安裝位置研究
    科技視界(2015年28期)2015-10-14 10:37:52
    慕課背景下應(yīng)用型本科院校數(shù)值分析課程的教學(xué)改革實(shí)踐
    科技視界(2015年26期)2015-09-11 13:39:38
    亚洲人成网站在线播放欧美日韩| 夜夜躁狠狠躁天天躁| 亚洲人成伊人成综合网2020| 久久青草综合色| 午夜激情av网站| 美女大奶头视频| 欧美激情极品国产一区二区三区| 色播亚洲综合网| 久久这里只有精品19| 91成人精品电影| 黄色视频不卡| 亚洲一卡2卡3卡4卡5卡精品中文| 国产精品电影一区二区三区| 午夜激情av网站| 精品国产一区二区三区四区第35| 亚洲av熟女| bbb黄色大片| 日韩大尺度精品在线看网址 | 美女午夜性视频免费| 国产精品精品国产色婷婷| 在线观看日韩欧美| 一级毛片高清免费大全| 韩国精品一区二区三区| 日韩视频一区二区在线观看| 精品熟女少妇八av免费久了| 波多野结衣巨乳人妻| 身体一侧抽搐| 国产一区二区激情短视频| 在线观看午夜福利视频| 国产高清videossex| 欧美不卡视频在线免费观看 | 老司机靠b影院| 亚洲,欧美精品.| 电影成人av| 成熟少妇高潮喷水视频| 国产精品99久久99久久久不卡| 亚洲第一av免费看| 色播亚洲综合网| 欧美另类亚洲清纯唯美| 变态另类成人亚洲欧美熟女 | 日韩大尺度精品在线看网址 | 性欧美人与动物交配| 亚洲午夜理论影院| 无遮挡黄片免费观看| 天天添夜夜摸| 少妇的丰满在线观看| 在线观看舔阴道视频| 男人舔女人的私密视频| 久久中文字幕人妻熟女| 精品欧美一区二区三区在线| av免费在线观看网站| 国产主播在线观看一区二区| 日日干狠狠操夜夜爽| 亚洲人成伊人成综合网2020| 欧美国产精品va在线观看不卡| 欧美av亚洲av综合av国产av| 最近最新中文字幕大全免费视频| 亚洲天堂国产精品一区在线| 国产成人欧美| 国产精品久久久人人做人人爽| 国产精品 国内视频| www.精华液| 成年版毛片免费区| 亚洲,欧美精品.| 欧美成狂野欧美在线观看| 9色porny在线观看| 自线自在国产av| 亚洲精品国产精品久久久不卡| 高清毛片免费观看视频网站| e午夜精品久久久久久久| 午夜免费激情av| 日韩免费av在线播放| 此物有八面人人有两片| 国产亚洲精品第一综合不卡| 黄色 视频免费看| 国产三级在线视频| 精品久久久久久,| 国产精品日韩av在线免费观看 | 欧美中文综合在线视频| 两性午夜刺激爽爽歪歪视频在线观看 | 丝袜美足系列| 激情视频va一区二区三区| 久久久国产欧美日韩av| 亚洲av电影在线进入| 1024视频免费在线观看| 大香蕉久久成人网| 性少妇av在线| 91精品国产国语对白视频| 国产精品亚洲一级av第二区| 美女免费视频网站| 欧美日韩亚洲综合一区二区三区_| 日本 av在线| 久久久久久久午夜电影| 久久精品91无色码中文字幕| 欧美色欧美亚洲另类二区 | 亚洲中文字幕日韩| 国产不卡一卡二| 美女大奶头视频| 1024香蕉在线观看| 成在线人永久免费视频| 大香蕉久久成人网| 国产高清videossex| 中文字幕久久专区| 在线国产一区二区在线| 精品国产超薄肉色丝袜足j| 亚洲全国av大片| 国产在线观看jvid| 亚洲熟妇中文字幕五十中出| 免费观看人在逋| 视频在线观看一区二区三区| 中文字幕人妻熟女乱码| 精品日产1卡2卡| 精品第一国产精品| 咕卡用的链子| 欧美成人午夜精品| 亚洲国产精品sss在线观看| 美女高潮喷水抽搐中文字幕| 69av精品久久久久久| 日本一区二区免费在线视频| 亚洲午夜理论影院| avwww免费| 最近最新免费中文字幕在线| 日日摸夜夜添夜夜添小说| a级毛片在线看网站| 国产极品粉嫩免费观看在线| 丝袜美足系列| 婷婷精品国产亚洲av在线| 香蕉丝袜av| 亚洲国产中文字幕在线视频| 黑人巨大精品欧美一区二区mp4| 亚洲男人天堂网一区| 别揉我奶头~嗯~啊~动态视频| 高清在线国产一区| 久久久久国产精品人妻aⅴ院| 伊人久久大香线蕉亚洲五| 每晚都被弄得嗷嗷叫到高潮| 99久久综合精品五月天人人| 精品国产乱子伦一区二区三区| 一区二区三区国产精品乱码| 老熟妇乱子伦视频在线观看| 亚洲片人在线观看| 男女下面插进去视频免费观看| 最近最新免费中文字幕在线| 国产97色在线日韩免费| 久久人人爽av亚洲精品天堂| 激情在线观看视频在线高清| 久久久久久人人人人人| 少妇被粗大的猛进出69影院| 久久人妻福利社区极品人妻图片| 久久性视频一级片| 97超级碰碰碰精品色视频在线观看| xxx96com| 亚洲一区中文字幕在线| 亚洲中文字幕日韩| 午夜免费成人在线视频| 成人三级黄色视频| 99在线人妻在线中文字幕| 亚洲精华国产精华精| 精品无人区乱码1区二区| svipshipincom国产片| 精品卡一卡二卡四卡免费| 久久国产精品影院| 精品久久久久久,| 午夜亚洲福利在线播放| 首页视频小说图片口味搜索| 国产成人av激情在线播放| 中文字幕久久专区| 日日爽夜夜爽网站| 免费在线观看日本一区| 欧美久久黑人一区二区| 成人永久免费在线观看视频| 精品国产亚洲在线| 亚洲片人在线观看| 国产精品九九99| 桃红色精品国产亚洲av| 亚洲一区高清亚洲精品| 中文字幕高清在线视频| 激情视频va一区二区三区| 色综合亚洲欧美另类图片| 两个人视频免费观看高清| 亚洲av熟女| 首页视频小说图片口味搜索| 精品免费久久久久久久清纯| 国产精品野战在线观看| 色播在线永久视频| 人妻久久中文字幕网| 亚洲九九香蕉| 精品国产乱码久久久久久男人| 黄片播放在线免费| 午夜久久久在线观看| 国产97色在线日韩免费| 一a级毛片在线观看| 亚洲人成伊人成综合网2020| 亚洲精品中文字幕一二三四区| 亚洲欧美日韩无卡精品| 成人18禁高潮啪啪吃奶动态图| 国产欧美日韩一区二区精品| 女同久久另类99精品国产91| √禁漫天堂资源中文www| 国产精品一区二区三区四区久久 | 日韩一卡2卡3卡4卡2021年| 国产精品精品国产色婷婷| 久久国产乱子伦精品免费另类| 国产视频一区二区在线看| 九色亚洲精品在线播放| 欧美国产精品va在线观看不卡| 伦理电影免费视频| 黑人欧美特级aaaaaa片| 一卡2卡三卡四卡精品乱码亚洲| 法律面前人人平等表现在哪些方面| 成熟少妇高潮喷水视频| 高清在线国产一区| 伊人久久大香线蕉亚洲五| 亚洲精品中文字幕一二三四区| 久久国产精品人妻蜜桃| 久久欧美精品欧美久久欧美| 桃红色精品国产亚洲av| 热re99久久国产66热| 村上凉子中文字幕在线| 好男人在线观看高清免费视频 | 日韩精品青青久久久久久| 成在线人永久免费视频| 亚洲天堂国产精品一区在线| 国产男靠女视频免费网站| 亚洲国产精品成人综合色| 两个人看的免费小视频| 女人被躁到高潮嗷嗷叫费观| 日本免费a在线| 正在播放国产对白刺激| 91字幕亚洲| 久久人妻福利社区极品人妻图片| 夜夜躁狠狠躁天天躁| 国产野战对白在线观看| 一级,二级,三级黄色视频| 国产精品秋霞免费鲁丝片| 国产主播在线观看一区二区| 亚洲无线在线观看| av免费在线观看网站| 亚洲电影在线观看av| 欧美一级a爱片免费观看看 | 亚洲成av人片免费观看| aaaaa片日本免费| 久9热在线精品视频| bbb黄色大片| 欧美日韩福利视频一区二区| 午夜成年电影在线免费观看| 亚洲人成电影免费在线| 成年女人毛片免费观看观看9| 非洲黑人性xxxx精品又粗又长| 人人澡人人妻人| 久久久久久久久免费视频了| 99国产精品99久久久久| 欧美日本视频| 黑人操中国人逼视频| 91成人精品电影| 午夜福利视频1000在线观看 | 国产亚洲欧美98| 中文字幕av电影在线播放| 天天躁夜夜躁狠狠躁躁| 国产亚洲av嫩草精品影院| 天天一区二区日本电影三级 | 国产亚洲精品av在线| 搡老岳熟女国产| 香蕉国产在线看| 99国产精品一区二区三区| 欧美黑人欧美精品刺激| 一夜夜www| 国产精品影院久久| 亚洲熟妇中文字幕五十中出| 国产又爽黄色视频| 国产高清视频在线播放一区| 身体一侧抽搐| 一级毛片女人18水好多| 午夜免费鲁丝| 又大又爽又粗| 欧美日韩福利视频一区二区| 女人被狂操c到高潮| 日韩欧美国产在线观看| 看免费av毛片| 国产精品一区二区精品视频观看| 丝袜在线中文字幕| 亚洲自拍偷在线| 久久精品国产亚洲av高清一级| www国产在线视频色| 99久久综合精品五月天人人| 亚洲午夜理论影院| 欧美国产精品va在线观看不卡| 亚洲男人的天堂狠狠| 亚洲国产精品合色在线| 男女下面进入的视频免费午夜 | 777久久人妻少妇嫩草av网站| 久久久国产欧美日韩av| 国产一区二区三区在线臀色熟女| 精品熟女少妇八av免费久了| 午夜免费鲁丝| 久久精品91蜜桃| av在线播放免费不卡| 国产精品 国内视频| 中文字幕色久视频| 香蕉久久夜色| 国产又色又爽无遮挡免费看| 久久精品国产亚洲av高清一级| 精品欧美一区二区三区在线| √禁漫天堂资源中文www| 超碰成人久久| 午夜免费成人在线视频| 黄色视频不卡| 久久久国产成人免费| 制服人妻中文乱码| 日日夜夜操网爽| 精品久久久久久久人妻蜜臀av | 十八禁人妻一区二区| 99久久综合精品五月天人人| 国产激情久久老熟女| 久久久久国产精品人妻aⅴ院| 757午夜福利合集在线观看| 色哟哟哟哟哟哟| 午夜福利高清视频| 99精品欧美一区二区三区四区| 国产亚洲欧美在线一区二区| 亚洲成av人片免费观看| 亚洲男人的天堂狠狠| 欧美日本视频| 精品一区二区三区视频在线观看免费| 久久天堂一区二区三区四区| 欧美中文综合在线视频| 99re在线观看精品视频| 少妇的丰满在线观看| 好男人电影高清在线观看| 男人的好看免费观看在线视频 | 9191精品国产免费久久| or卡值多少钱| 脱女人内裤的视频| 国产91精品成人一区二区三区| 操美女的视频在线观看| 神马国产精品三级电影在线观看 | a级毛片在线看网站| 亚洲自偷自拍图片 自拍| 琪琪午夜伦伦电影理论片6080| 两个人看的免费小视频| 亚洲av美国av| 最新美女视频免费是黄的| 日韩精品免费视频一区二区三区| 免费久久久久久久精品成人欧美视频| 亚洲国产看品久久| 精品久久久久久成人av| avwww免费| 成人国语在线视频| 亚洲精品一区av在线观看| 别揉我奶头~嗯~啊~动态视频| 欧美成人午夜精品| a在线观看视频网站| 少妇熟女aⅴ在线视频| 久久中文字幕人妻熟女| 咕卡用的链子| 国产av精品麻豆| 麻豆一二三区av精品| 啦啦啦 在线观看视频| 午夜免费鲁丝| 国产精品日韩av在线免费观看 | 日本免费a在线| 亚洲一区二区三区不卡视频| 香蕉丝袜av| 搡老熟女国产l中国老女人| videosex国产| 99国产精品99久久久久| 黄色视频,在线免费观看| 亚洲av片天天在线观看| 叶爱在线成人免费视频播放| 好看av亚洲va欧美ⅴa在| 国产一区二区三区在线臀色熟女| 午夜福利影视在线免费观看| 高清黄色对白视频在线免费看| 又黄又粗又硬又大视频| 18禁美女被吸乳视频| 国产伦一二天堂av在线观看| 午夜亚洲福利在线播放| netflix在线观看网站| 国产私拍福利视频在线观看| 美女 人体艺术 gogo| www.www免费av| 窝窝影院91人妻| 精品人妻在线不人妻| 免费在线观看黄色视频的| 真人做人爱边吃奶动态| 欧美另类亚洲清纯唯美| 好看av亚洲va欧美ⅴa在| 国产精品久久视频播放| 一区二区三区国产精品乱码| 在线视频色国产色| 18禁裸乳无遮挡免费网站照片 | 少妇的丰满在线观看| www.自偷自拍.com| 好男人在线观看高清免费视频 | 老司机深夜福利视频在线观看| 国产不卡一卡二| 国产色视频综合| 久久精品亚洲熟妇少妇任你| 国产精品,欧美在线| 日韩大码丰满熟妇| 91在线观看av| 国产91精品成人一区二区三区| www国产在线视频色| 国产精品久久久久久亚洲av鲁大| 99riav亚洲国产免费| 亚洲va日本ⅴa欧美va伊人久久| 男女床上黄色一级片免费看| 日本撒尿小便嘘嘘汇集6| 亚洲午夜精品一区,二区,三区| 国产蜜桃级精品一区二区三区| 黄色片一级片一级黄色片| 香蕉久久夜色| 亚洲片人在线观看| 国产区一区二久久| 亚洲成人免费电影在线观看| 亚洲精品中文字幕一二三四区| 99国产精品99久久久久| 宅男免费午夜| 国产精品亚洲av一区麻豆| 18禁观看日本| 一边摸一边抽搐一进一小说| 精品久久久精品久久久| 国产精品免费一区二区三区在线| 成年女人毛片免费观看观看9| 成人国产一区最新在线观看| 丝袜人妻中文字幕| 色老头精品视频在线观看| 大陆偷拍与自拍| 国产精品电影一区二区三区| 精品人妻在线不人妻| 精品欧美一区二区三区在线| 久久久久亚洲av毛片大全| 午夜福利视频1000在线观看 | 久久久久久久精品吃奶| 午夜日韩欧美国产| 啦啦啦 在线观看视频| 精品国产乱子伦一区二区三区| 国产不卡一卡二| www.精华液| 久热这里只有精品99| 最新美女视频免费是黄的| 亚洲熟妇熟女久久| 成人永久免费在线观看视频| 国产精品久久久久久人妻精品电影| 在线观看66精品国产| 好男人在线观看高清免费视频 | 大陆偷拍与自拍| 在线十欧美十亚洲十日本专区| 老鸭窝网址在线观看| 久久精品91蜜桃| 国产欧美日韩综合在线一区二区| 国产成人系列免费观看| 制服人妻中文乱码| 91成年电影在线观看| 日韩国内少妇激情av| 1024视频免费在线观看| 老熟妇仑乱视频hdxx| 久久久久久亚洲精品国产蜜桃av| 午夜a级毛片| 久久热在线av| 琪琪午夜伦伦电影理论片6080| 国产视频一区二区在线看| cao死你这个sao货| 亚洲午夜理论影院| 久久精品国产亚洲av香蕉五月| 亚洲久久久国产精品| 精品欧美一区二区三区在线| 悠悠久久av| 国产av又大| 一区二区三区激情视频| 日本 欧美在线| 亚洲av五月六月丁香网| 黑丝袜美女国产一区| 波多野结衣巨乳人妻| 1024香蕉在线观看| 欧美日韩亚洲国产一区二区在线观看| 纯流量卡能插随身wifi吗| 在线观看免费日韩欧美大片| 精品福利观看| 国产精品av久久久久免费| 久久国产亚洲av麻豆专区| 亚洲熟妇熟女久久| 看免费av毛片| 国产一区二区三区综合在线观看| 9色porny在线观看| 国产91精品成人一区二区三区| 亚洲精品国产色婷婷电影| 在线观看一区二区三区| 久久久久国产一级毛片高清牌| 此物有八面人人有两片| 别揉我奶头~嗯~啊~动态视频| 国产一区二区三区视频了| www.www免费av| 麻豆成人av在线观看| a级毛片在线看网站| 欧美精品亚洲一区二区| 国产精品久久久久久人妻精品电影| 少妇熟女aⅴ在线视频| 18禁美女被吸乳视频| 怎么达到女性高潮| 美女午夜性视频免费| 午夜福利视频1000在线观看 | 国产欧美日韩一区二区精品| 一级a爱视频在线免费观看| 欧美在线一区亚洲| 精品福利观看| 亚洲黑人精品在线| 成人精品一区二区免费| 日韩欧美三级三区| 国产高清videossex| 成人国产综合亚洲| 亚洲aⅴ乱码一区二区在线播放 | 日韩视频一区二区在线观看| 日韩欧美国产在线观看| av欧美777| 人妻久久中文字幕网| 人人妻人人澡欧美一区二区 | 午夜两性在线视频| av在线天堂中文字幕| 看片在线看免费视频| 免费高清视频大片| 国内久久婷婷六月综合欲色啪| 免费av毛片视频| 午夜久久久久精精品| 久久精品国产清高在天天线| 校园春色视频在线观看| 99re在线观看精品视频| 国产激情久久老熟女| 久久婷婷人人爽人人干人人爱 | 1024香蕉在线观看| 91大片在线观看| 亚洲成av片中文字幕在线观看| 黄片播放在线免费| 可以免费在线观看a视频的电影网站| 午夜福利18| 精品欧美一区二区三区在线| 一a级毛片在线观看| 日本在线视频免费播放| 国产一卡二卡三卡精品| 精品一区二区三区av网在线观看| 人人妻人人爽人人添夜夜欢视频| 99精品在免费线老司机午夜| 色综合站精品国产| 中文字幕高清在线视频| 俄罗斯特黄特色一大片| 天天一区二区日本电影三级 | 亚洲最大成人中文| 欧美在线一区亚洲| av免费在线观看网站| 亚洲男人天堂网一区| 欧美日韩一级在线毛片| 欧美日韩黄片免| 两性夫妻黄色片| 亚洲成人久久性| 欧美日韩一级在线毛片| 大型黄色视频在线免费观看| 久久久精品国产亚洲av高清涩受| 国产三级在线视频| 亚洲av电影不卡..在线观看| 精品一品国产午夜福利视频| 国产av一区在线观看免费| 亚洲激情在线av| 人人澡人人妻人| 女性生殖器流出的白浆| 亚洲午夜理论影院| 亚洲七黄色美女视频| 黄色成人免费大全| 亚洲五月婷婷丁香| 亚洲国产欧美日韩在线播放| 久99久视频精品免费| 男女午夜视频在线观看| 欧美一级毛片孕妇| 在线永久观看黄色视频| 9色porny在线观看| 久久久久久国产a免费观看| 好看av亚洲va欧美ⅴa在| 69av精品久久久久久| 久久香蕉精品热| 久9热在线精品视频| 人人妻,人人澡人人爽秒播| 99久久综合精品五月天人人| av视频免费观看在线观看| 国产成人av激情在线播放| 两个人视频免费观看高清| 天天躁夜夜躁狠狠躁躁| 国产精品一区二区在线不卡| 一夜夜www| 亚洲国产欧美日韩在线播放| 多毛熟女@视频| 亚洲av电影在线进入| 妹子高潮喷水视频| 精品国产一区二区久久| 亚洲国产精品sss在线观看| 一卡2卡三卡四卡精品乱码亚洲| 少妇被粗大的猛进出69影院| 两个人看的免费小视频| 成人免费观看视频高清| 看片在线看免费视频| 大陆偷拍与自拍| 岛国在线观看网站| 午夜福利在线观看吧| 亚洲午夜精品一区,二区,三区| 一边摸一边抽搐一进一小说| 不卡一级毛片| 别揉我奶头~嗯~啊~动态视频| 国产男靠女视频免费网站| 久久久国产成人精品二区| 亚洲精品久久成人aⅴ小说| 久久热在线av| 久久久精品欧美日韩精品| 成人国产综合亚洲| 国产欧美日韩一区二区精品| av视频免费观看在线观看| 国产成+人综合+亚洲专区| 国产精品美女特级片免费视频播放器 |