許 澎, 許思傳, 唐軍英, 高 源
(1. 同濟(jì)大學(xué) 汽車(chē)學(xué)院, 上海 201804; 2. 安徽明天氫能科技股份有限公司, 安徽 六安 237000;3. 同濟(jì)大學(xué) 機(jī)械與能源工程學(xué)院, 上海 201804)
質(zhì)子交換膜燃料電池(PEMFC,proton exchange membrane fuel cell)因其零污染、能量轉(zhuǎn)換效率高、啟動(dòng)響應(yīng)迅速和燃料來(lái)源廣泛等特點(diǎn)被視為未來(lái)最有前景的車(chē)用動(dòng)力來(lái)源之一[1-2].作為動(dòng)力源的燃料電池欲實(shí)現(xiàn)大規(guī)模應(yīng)用須經(jīng)受高電位、啟停、反極、極限電流和冷啟動(dòng)等嚴(yán)苛工況.其中,零下冷啟動(dòng)問(wèn)題來(lái)源于低溫啟動(dòng)過(guò)程中冰的形成覆蓋電化學(xué)活性位點(diǎn)和阻塞反應(yīng)氣傳輸通道,造成電化學(xué)反應(yīng)啟動(dòng)失敗和耐久性降低[3-4].目前,有關(guān)燃料電池低溫啟動(dòng)過(guò)程的數(shù)值模擬研究均聚焦在宏觀尺度下的水熱傳遞、反應(yīng)氣傳質(zhì)和相變等行為,以求解一維和三維微分或偏微分方程組為主要形式.如Jiao等[5-6]采用有限體積法的思想對(duì)燃料電池內(nèi)的物理和化學(xué)現(xiàn)象建立一系列偏微分方程組,實(shí)現(xiàn)對(duì)輸出電壓和冰體積分?jǐn)?shù)等宏觀量的直接觀察;許等[7]基于分層集總參數(shù)方法對(duì)燃料電池單電池進(jìn)行一維冷啟動(dòng)數(shù)值建模,在單電池內(nèi)每一個(gè)組件中心處進(jìn)行建模和研究.但上述文獻(xiàn)[5-7]研究考察現(xiàn)象均為整片單電池宏觀尺度下的物理和化學(xué)行為.此外,目前有關(guān)燃料電池低溫啟動(dòng)實(shí)驗(yàn)研究中可挖掘數(shù)據(jù)較少,除輸出電壓、歐姆阻抗和溫度等結(jié)果參數(shù)外,若不借助X射線(xiàn)成像等高成本手段無(wú)法觀察電池內(nèi)部機(jī)理.再者,介觀模擬尺度下燃料電池多孔介質(zhì)內(nèi)低溫結(jié)冰研究和報(bào)道還未見(jiàn).
相較宏觀方法,由于演化方程簡(jiǎn)單、易并行計(jì)算和邊界條件處理簡(jiǎn)單等優(yōu)勢(shì),格子Boltzmann(LB)方法在多孔介質(zhì)兩相流等方面逐漸發(fā)揮了廣泛應(yīng)用[8-10].作為介觀模擬尺度的代表之一,格子Boltzmann方法觀察流體分子的速度分布函數(shù),通過(guò)考察分子的時(shí)空演化過(guò)程并根據(jù)分布函數(shù)與宏觀量之間的關(guān)系來(lái)獲得宏觀流體信息.格子Boltzmann方法與宏觀模型的相同點(diǎn)是兩者都為微觀分子的統(tǒng)計(jì)學(xué)量,不在乎分子對(duì)描述對(duì)象的影響,不同之處在于前者沒(méi)有連續(xù)性假設(shè);與微觀分子模型的相同點(diǎn)在于兩者都從微觀層次觀察流體分子的運(yùn)動(dòng)狀態(tài),不同之處在于后者表現(xiàn)單個(gè)分子的行為,前者是微觀分子的統(tǒng)計(jì)行為.
近年來(lái),格子Boltzmann在多孔介質(zhì)液固相變研究領(lǐng)域有些發(fā)展,主要有表征體元(REV)尺度和孔隙尺度.REV尺度的研究對(duì)象為一個(gè)控制體,求解控制體內(nèi)體積平均宏觀物理量,不關(guān)心微觀結(jié)構(gòu),僅依賴(lài)于統(tǒng)計(jì)參數(shù);孔隙尺度下的格子Boltzmann方法研究對(duì)象為流體分子團(tuán),多孔介質(zhì)的骨架為流場(chǎng)和溫度場(chǎng)的邊界.目前有關(guān)格子Boltzmann方法在多孔介質(zhì)相變領(lǐng)域的研究較少,且多數(shù)集中在REV尺度下的融化過(guò)程[11-14],如Gao等[12]對(duì)REV尺度下多孔介質(zhì)方腔內(nèi)的融化過(guò)程進(jìn)行建模研究,考察孔隙率、Darcy數(shù)和Rayleigh數(shù)對(duì)自然對(duì)流的影響.此外,最近一些學(xué)者[15-17]對(duì)孔隙尺度下多孔介質(zhì)內(nèi)融化相變問(wèn)題展開(kāi)了相關(guān)研究,如杲東彥等[15]研究了正方形有序排列組成的多孔介質(zhì)中Rayleigh數(shù)和Prandtl數(shù)對(duì)融化的影響.上述研究中所構(gòu)建的多孔介質(zhì)結(jié)構(gòu)缺乏真實(shí)物理背景;此外,由于相變儲(chǔ)能系統(tǒng)中液固兩相的熱物理參數(shù)差別較小,一些研究中假設(shè)液固兩相熱物理參數(shù)相同,但液態(tài)水對(duì)應(yīng)液相和固相的比熱容和導(dǎo)熱系數(shù)差異較大,因此上述研究缺乏對(duì)液固兩相不同熱物理性參數(shù)的考察.本文通過(guò)構(gòu)造燃料電池用真實(shí)氣體擴(kuò)散層三維微孔隙結(jié)構(gòu),嚴(yán)格考察凝固模型中液態(tài)水、固體冰和燃料電池用碳纖維的熱物理參數(shù),首次采用基于格子Boltzmann方法的液固相變模型研究了燃料電池氣體擴(kuò)散內(nèi)液態(tài)水的結(jié)冰過(guò)程,從介觀角度加深對(duì)燃料電池多孔介質(zhì)內(nèi)結(jié)冰機(jī)理的理解.
目前,基于格子Boltzmann方法的相變模型主要有相場(chǎng)法和焓方法.本文采用Huang學(xué)者[18-19]提出的固液相變總焓格子Boltzmann方法,該方法可以高效處理非線(xiàn)性潛熱源項(xiàng),避免焓值迭代計(jì)算和線(xiàn)性方程組求解,且可以實(shí)現(xiàn)相界面處變熱物性的處理.由于結(jié)冰為相變問(wèn)題,計(jì)算區(qū)域涉及液相和固相,需特殊處理以實(shí)現(xiàn)固相速度無(wú)滑移條件.網(wǎng)格點(diǎn)處的固相率實(shí)質(zhì)上為單元內(nèi)的平均固相率,因此,Huang學(xué)者提出了體積格子Boltzmann格式處理固相速度無(wú)滑移條件.鑒于LB方程的多松弛時(shí)間(MRT)處理可以使用多個(gè)松弛時(shí)間參數(shù),具備多個(gè)調(diào)節(jié)參數(shù),在數(shù)值穩(wěn)定性和基本物理原理等方面具有明顯優(yōu)勢(shì),本研究中對(duì)于速度場(chǎng)和溫度場(chǎng),都采用MRT形式.
(1)
格子Boltzmann方程一般分為碰撞過(guò)程和遷移過(guò)程.在計(jì)算過(guò)程中,碰撞步在矩空間執(zhí)行,遷移步在速度空間執(zhí)行.碰撞步在矩空間的表達(dá)式為
(2)
(3)
S=diag(s0,se,sε,sj,sq,sp,sp)
(4)
遷移步在速度空間下的表達(dá)式為
(5)
平衡態(tài)矩函數(shù)meq(x,t)為
(6)
式中:ρ0為參考密度;c為格子速度,c=δx/δt,δx為格子單位長(zhǎng)度;u為速度矢量.矩空間的離散作用力項(xiàng)為
(7)
固相的出現(xiàn)不影響任意格點(diǎn)處的質(zhì)量守恒關(guān)系,格點(diǎn)處的宏觀密度根據(jù)臨時(shí)速度分布函數(shù)計(jì)算,表達(dá)式為
(8)
由于固相不運(yùn)動(dòng),假設(shè)介觀尺度下固相對(duì)應(yīng)的速度分布函數(shù)一直處于平衡狀態(tài),固相使得臨時(shí)速度分布函數(shù)中固相對(duì)應(yīng)的部分趨于平衡態(tài).因此,格點(diǎn)處速度分布函數(shù)表示為
(9)
式中:fs為固相率;us為固相速率.
根據(jù)速度分布函數(shù)計(jì)算宏觀速度,即:
(10)
格子Boltzmann方程可以通過(guò)Chapman-Enskog多尺度分析推導(dǎo)出宏觀方程.速度分布函數(shù)的格子Boltzmann方程對(duì)應(yīng)的二階近似宏觀方程為
(11)
(12)
在Jiaung等[20-21]的固液相變LB模型中,潛熱源項(xiàng)以離散源項(xiàng)的形式添加在LB方程中,顯式LB方程變?yōu)殡[式,故而需要采用焓值迭代計(jì)算或線(xiàn)性方程組求解,計(jì)算量巨大.Huang等學(xué)者[18-19]構(gòu)建了一種新的總焓LB模型,總焓分布函數(shù)gi(x,t)的MRTLB方程為
gi(x+eiδt,t+δt)=
(13)
與速度分布函數(shù)的計(jì)算過(guò)程類(lèi)似,總焓分布函數(shù)的碰撞步亦在矩空間進(jìn)行,遷移步亦在速度空間進(jìn)行,即:
(14)
(15)
(16)
總焓平衡態(tài)矩函數(shù)neq(x,t)和參考比熱容Cp,ref、總焓H、溫度T和比熱容Cp有關(guān),表達(dá)式為
(17)
與之對(duì)應(yīng)的總焓平衡態(tài)分布函數(shù)為
(18)
總焓H的計(jì)算式為
(19)
溫度T和固相率fs由總焓H確定,即:
(20)
(21)
式(20)~(21)中:Hs、Hl分別為固態(tài)和液態(tài)的總焓.
與速度分布函數(shù)類(lèi)似,總焓MRT LB模型亦可以通過(guò)Chapman-Enskog分析推出宏觀方程,即:
(22)
宏觀熱擴(kuò)散系數(shù)為
(23)
邊界條件在數(shù)值模擬研究中至關(guān)重要.由于格點(diǎn)處的物理量通常代表網(wǎng)格單元內(nèi)的平均值,基于半步長(zhǎng)反彈思想,本文對(duì)速度無(wú)滑移邊界條件和溫度邊界條件分別采用半步長(zhǎng)反彈格式[22-23]和鏡像半步長(zhǎng)反彈格式[24-25]處理.如圖1所示,壁面格點(diǎn)xw離邊界網(wǎng)格點(diǎn)xb半個(gè)網(wǎng)格步長(zhǎng),為δx/2,以圖1中左邊邊界網(wǎng)格點(diǎn)xb為例,第5方向的未知速度分布函數(shù)fi(xb,t+δt)和未知總焓分布函數(shù)gi(xb,t+δt)可以表示為
(24)
(25)
(26)
圖1 邊界條件處理示意圖
對(duì)于多孔介質(zhì)微觀結(jié)構(gòu)的重構(gòu)主要有兩種方法:基于“圖像"的重構(gòu)方法和基于“隨機(jī)"的重構(gòu)方法[26].前者通常采用激光掃描顯微鏡和X射線(xiàn)斷層攝影術(shù)等方法對(duì)多孔介質(zhì)材料進(jìn)行掃描,通過(guò)對(duì)掃描得到的二維圖像進(jìn)行計(jì)算機(jī)處理后得到微觀結(jié)構(gòu),該方法基于實(shí)際材料,結(jié)構(gòu)較為真實(shí),但成本較高,受圖像處理技術(shù)的限制;后者采用計(jì)算機(jī)的隨機(jī)數(shù)字發(fā)生器生成材料典型的位置和方向特征,并提前設(shè)定規(guī)則來(lái)生成材料結(jié)構(gòu),該方法基于計(jì)算機(jī),成本低,速度快,雖只能接近真實(shí)材料的微觀孔隙結(jié)構(gòu),但對(duì)于目前的研究已足夠.以下將介紹本文研究中氣體擴(kuò)散層微孔隙結(jié)構(gòu)的重構(gòu)方法及過(guò)程.
在孔隙尺度下,多孔介質(zhì)可以看成由無(wú)數(shù)根半徑相同的碳纖維雜亂無(wú)章的平鋪層疊起來(lái)的.因此,在幾何上可以簡(jiǎn)化為:無(wú)數(shù)根圓柱體重疊相交;圓柱體可以通過(guò)它的軸線(xiàn)方向向量u=(u1,u2,u3)、軸線(xiàn)通過(guò)的一點(diǎn)P0(x0,y0,z0)和圓柱體半徑r唯一確定,如圖2所示;任意一根圓柱體可以定義為距離圓柱體軸線(xiàn)距離不大于圓柱體半徑r的所有點(diǎn)的集合P(x,y,z);圓柱體內(nèi)所有點(diǎn)的集合占據(jù)空間總體積比例為1與孔隙率ε的差值.
為了簡(jiǎn)化多孔介質(zhì)孔隙結(jié)構(gòu)的構(gòu)建模型,作以下假設(shè):① 碳纖維呈筆直圓柱體形狀;② 所有碳纖維半徑完全相同;③ 碳纖維存在重疊交叉的可能.圖3為本文多孔介質(zhì)隨機(jī)重構(gòu)方法流程圖.首先需要給出多孔介質(zhì)的目標(biāo)孔隙率ε、纖維半徑r、各項(xiàng)同性系數(shù)K和三維結(jié)構(gòu)的長(zhǎng)寬高Lx·Ly·Lz,并且初始化孔隙率E=1.然后,利用Matlab中的隨機(jī)函數(shù)rand隨機(jī)生成圓柱體軸線(xiàn)的方向向量和通過(guò)軸線(xiàn)的一個(gè)點(diǎn)坐標(biāo),過(guò)程如下:
圖2 PEMFC碳纖維圓柱體模型
圖3 隨機(jī)重構(gòu)氣體擴(kuò)散層(GDL)微孔隙結(jié)構(gòu)算法流程圖
Fig.3 Flowchart of stochastic reconstruction algorithm for gas diffusion layer(GDL) micro-porous structure
微尺度下碳纖維的重疊層交并不是均勻分布和各向同性的,多孔介質(zhì)材料的各向異性可以看成制造過(guò)程中對(duì)材料施加外部壓力導(dǎo)致的.如圖4所示,各向同性系數(shù)K可以定位為碳纖維最終厚度與初始厚度的比值,K=1代表完全各向同性,即最終厚度等于初始厚度.各向同性系數(shù)表達(dá)式為
(27)
在各向同性系數(shù)已知情況下,方向角度可從式(28)得到:
θ=cos-1(Kcosθ′)
(28)
式中:θ′為初始碳纖維的方向角度,可由Matlab的rand函數(shù)給出;θ為碳纖維壓縮后方向角度.
圖4 壓縮對(duì)碳纖維方向角的影響
根據(jù)方向角度得出方向向量(u1,u2,u3),再通過(guò)圓柱體軸線(xiàn)的正交基準(zhǔn)法向矢量和隨機(jī)數(shù)rand函數(shù),得到通過(guò)軸線(xiàn)的隨機(jī)點(diǎn)(x0,y0,z0).在圓柱體方向向量、通過(guò)軸線(xiàn)的隨機(jī)點(diǎn)以及碳纖維半徑都確定后,通過(guò)while循環(huán)函數(shù)對(duì)三維空間中所有點(diǎn)遍歷并判斷當(dāng)前點(diǎn)距離圓柱體軸線(xiàn)的距離是否小于碳纖維半徑r,小于作記號(hào)f=1,否則作記號(hào)f=0,當(dāng)前點(diǎn)距離圓柱軸線(xiàn)的距離公式見(jiàn)式(16).每生成一個(gè)圓柱體后,須判斷當(dāng)前孔隙率值,若E>ε,繼續(xù)執(zhí)行生成方向向量和隨機(jī)點(diǎn),反之,程序停止.
(29)
本文取碳纖維半徑r=3.5 μm,目標(biāo)孔隙率分別為0.5、0.6、0.7、0.8和0.9,各向同性系數(shù)K=0.8,三維尺寸Lx·Ly·Lz=200 μm×200 μm×200 μm,模型比例為1個(gè)網(wǎng)格長(zhǎng)度代表1 μm.圖5為對(duì)應(yīng)孔隙率分別為0.6的氣體擴(kuò)散層微孔隙三維結(jié)構(gòu)視覺(jué)圖,本文研究的二維氣體擴(kuò)散層為x方向1/2處的截取的平面.
為了驗(yàn)證本文基于格子Boltzmann方法的燃料電池GDL內(nèi)結(jié)冰模型的精確性和有效性,針對(duì)液、固和多孔介質(zhì)三種材料的不同熱物理參數(shù)進(jìn)行了一維半無(wú)限大空間凝固熱傳導(dǎo)、二維直角區(qū)域凝固和二維介質(zhì)方腔凝固的數(shù)值試驗(yàn),嚴(yán)格選取對(duì)應(yīng)真實(shí)燃料電池GDL層內(nèi)碳纖維、液態(tài)水和固體冰的熱物理參數(shù).注意,如沒(méi)有特殊聲明,本文涉及的數(shù)值皆為格子單位.
圖5 孔隙率0.6氣體擴(kuò)散層微孔隙三維結(jié)構(gòu)視圖
為了驗(yàn)證本文的凝固模型可以處理液固相變分界線(xiàn)兩側(cè)不同熱物理參數(shù)問(wèn)題,即液態(tài)水比熱容Cpl和冰比熱容Cps、液態(tài)水導(dǎo)熱系數(shù)λl和冰λs不相同,進(jìn)行了一維半無(wú)限大凝固熱傳導(dǎo)數(shù)值實(shí)驗(yàn).如圖6所示,一維半無(wú)限大空間初始狀態(tài)下充滿(mǎn)溫度Ti=1的液態(tài),零時(shí)刻左側(cè)壁面溫度Tb驟降至0,液固分界線(xiàn)上的溫度維持在Tm=0.5,液固分界線(xiàn)從左往右逐漸移動(dòng).相關(guān)參數(shù)設(shè)置如下:Cpl=1.0,Cps=2.0,Ss=0.004,Sl=0.008,λl=0.15,λs=0.6,其中液態(tài)水和冰Stefan數(shù)分別為Sl=Cpl(Ti-Tm)/L.Ss=Cps(Tm-Tb)/L.注意,本文所有的數(shù)值模擬均采用dx=dt=1,松弛因子設(shè)置為σ0=1.0,σ0=1/τg,σε=σe,σq=σj.圖6b為對(duì)應(yīng)不同時(shí)刻下的溫度θ隨坐標(biāo)x變化情況,LBM解和解析解[27]吻合良好.
為了驗(yàn)證模型的使用適用于二維問(wèn)題,進(jìn)行了二維直角區(qū)域凝固數(shù)值試驗(yàn).計(jì)算域和溫度邊界條件如圖7所示.初始狀態(tài)下,方腔內(nèi)充滿(mǎn)溫度Ti=0.3的液體,零時(shí)刻開(kāi)始,方腔的左壁面和下壁面溫度驟降至Tb=-1.0,方腔的上壁面和右壁面保持絕熱狀態(tài),相關(guān)參數(shù)設(shè)置如下:Cpl=Cps=1.0,Sl=Ss=Cpl(Tm-Tb)/L=4.0,λl=λs=1.0.與Lin等的解析解[27]對(duì)比發(fā)現(xiàn),Fo=0.25時(shí)刻的等溫線(xiàn)和液固界面線(xiàn)吻合良好.
a 示意圖
b 溫度分布
為了驗(yàn)證本文模型適用于研究多孔介質(zhì)方腔,進(jìn)行了二維介質(zhì)方腔凝固的數(shù)值試驗(yàn).注意,由于引入了多孔介質(zhì),式(17)和式(18)中參考比熱容Cp,ref定義為[28]:
Cp,ref=3CpsCplCpp/(CpsCpl+CplCpp+CpsCpp)
(30)
式中:Cpp多孔介質(zhì)比熱容.
a 示意圖
b 溫度云圖
c 液固分界線(xiàn)(Fo=0.25)
d 等溫線(xiàn)
a 示意圖
b 液固分界線(xiàn)
首先,介紹格子單位和物理單位之間的換算關(guān)系.在本文的凝固相變模型中,Stefan數(shù)、量綱一溫度θ和Fourier數(shù)F0的定義見(jiàn)式(31):
(31a)
(31b)
(31c)
式中:α為導(dǎo)熱系數(shù);t為時(shí)間.
表1 碳纖維、液態(tài)水和冰熱物理參數(shù)
在本文中,冰體積分?jǐn)?shù)定義為孔隙內(nèi)冰體積占據(jù)多孔介質(zhì)中孔隙體積百分比.圖9表示不同孔隙率氣體擴(kuò)散層內(nèi)結(jié)冰過(guò)程中冰體積分?jǐn)?shù)隨量綱一時(shí)間F0的變化關(guān)系.從圖中可以看出,在冰體積分?jǐn)?shù)小于0.1時(shí),冰體積分?jǐn)?shù)和量綱一時(shí)間F0呈現(xiàn)斜率較大的線(xiàn)性關(guān)系,大于0.1時(shí),呈現(xiàn)斜率較低的線(xiàn)性關(guān)系.孔隙率越大,初始狀態(tài)液態(tài)水體積分?jǐn)?shù)越大,相同條件下完全結(jié)冰所用時(shí)間越長(zhǎng).
圖9 不同孔隙率氣體擴(kuò)散層內(nèi)結(jié)冰過(guò)程中冰體積分?jǐn)?shù)與F0變化關(guān)系
圖10為對(duì)應(yīng)孔隙率0.6的氣體擴(kuò)散層方腔內(nèi)50 %高度處自左側(cè)壁向右側(cè)壁的溫度分布隨時(shí)間變化情況,初始狀態(tài)即F0<0時(shí),自左至右溫度都為θ=1,F0=0時(shí),左側(cè)壁溫度驟降至θ=0.溫度分布θ在不同時(shí)刻和不同區(qū)域表現(xiàn)出不同的變化趨勢(shì),比如當(dāng)F0=0.01時(shí),在θ<0.5時(shí)為固相,且溫度和位置呈現(xiàn)線(xiàn)性變化趨勢(shì),θ>0.5時(shí)為液相,溫度θ<1時(shí)表現(xiàn)出拋物線(xiàn)形狀;F0=0.72時(shí),θ<0.5和θ>0.5區(qū)域溫度都呈現(xiàn)線(xiàn)性變化趨勢(shì).此外,還可以發(fā)現(xiàn),當(dāng)某一時(shí)刻下的液相溫度接近θ=0.5時(shí),液固分界線(xiàn)右移速度放慢,原因是溫度梯度較小,熱傳導(dǎo)放緩.
圖10 孔隙率為0.6 氣體擴(kuò)散層結(jié)冰過(guò)程中方腔1/2高度處溫度分布隨時(shí)間變化情況
Fig.10 Variation of dimensionless temperatureθwith time in 1/2 height for gas diffusion layer with porosity 0.6
對(duì)于均勻充滿(mǎn)一種流體的方腔,在其左側(cè)壁受低溫?zé)醾鲗?dǎo)過(guò)程中,根據(jù)傳熱學(xué)原理,等溫線(xiàn)應(yīng)是垂直分布的.圖11為孔隙率0.6的GDL中在F0=0.72時(shí)刻的溫度分布和液固分界線(xiàn)分布.由于相變溫度θ=0.5,對(duì)應(yīng)等溫線(xiàn)中θ=0.5處為液固分界線(xiàn)點(diǎn).對(duì)于等溫線(xiàn)θ=0.5左側(cè)的固體區(qū)域來(lái)說(shuō),由于碳纖維的導(dǎo)熱系數(shù)大于固體冰的導(dǎo)熱系數(shù),導(dǎo)致碳纖維處的等溫線(xiàn)比其他區(qū)域往右偏移,因此,對(duì)應(yīng)等溫線(xiàn)θ=0.5是凹凸不平的.圖11b為液固分界線(xiàn)示意圖,分界線(xiàn)左側(cè)為結(jié)冰區(qū)域,右側(cè)為未結(jié)冰區(qū)域.
通過(guò)計(jì)算不同孔隙率的氣體擴(kuò)散層方腔內(nèi)液態(tài)水結(jié)冰過(guò)程可以發(fā)現(xiàn),對(duì)應(yīng)氣體擴(kuò)散層孔隙率分別為0.5、0.6、0.7、0.8和0.9時(shí),完全結(jié)冰所用量綱一時(shí)間F0分別為2.67、3.11、3.68、4.31和4.84.孔隙率越大,GDL完全結(jié)冰所用時(shí)間越久,因?yàn)榭紫堵试酱?液態(tài)水含量越多.此外,液固分界線(xiàn)基本呈現(xiàn)豎直狀態(tài)并伴隨凹凸不平現(xiàn)象逐漸從左向右移動(dòng).
自然對(duì)流是自然界中較為普遍的一種物理現(xiàn)象.由于流體內(nèi)部存在溫差,使得各部分流體的密度不同,溫度較高的流體密度小會(huì)上浮,溫度低的流體密度大會(huì)下沉,引起流體內(nèi)的自然對(duì)流現(xiàn)象,自然對(duì)流現(xiàn)象在受熱融化的相變問(wèn)題中較為明顯且研究較多.為了研究自然對(duì)流現(xiàn)象在燃料電池氣體擴(kuò)散層受冷凝固問(wèn)題中的影響,本文對(duì)孔隙率從0.5~0.9氣體擴(kuò)散層方腔的左側(cè)受冷凝固過(guò)程進(jìn)行了研究.
圖11 孔隙率為0.6的 氣體擴(kuò)散層在F0=0.72的等溫線(xiàn)和結(jié)冰情況
Fig.11 Isothermal curve and liquid-solid interface atF0=0.72 for gas diffusion layer with porosity 0.6
研究發(fā)現(xiàn),對(duì)應(yīng)氣體擴(kuò)散層的孔隙率分別為0.5、0.6、0.7、0.8和0.9,有自然對(duì)流情況下的結(jié)冰時(shí)間比無(wú)自然對(duì)流時(shí)分別減少0、0、0.001、0.001和0.007.原因是GDL孔隙率越大,未被碳纖維封閉的液相區(qū)域越大,自然對(duì)流現(xiàn)象也越明顯.圖12為對(duì)應(yīng)孔隙率為0.6和0.36時(shí)刻的自然對(duì)流現(xiàn)象,箭頭所指粗實(shí)線(xiàn)為液固分界線(xiàn),帶箭頭細(xì)實(shí)線(xiàn)為流線(xiàn).顯然,流線(xiàn)主要集中在孔隙較大且密集的區(qū)域.
圖12 孔隙率為0.6 氣體擴(kuò)散層內(nèi)結(jié)冰過(guò)程中在F0=0.36時(shí)刻的自然對(duì)流現(xiàn)象
Fig.12 Natural convection phenomenon atF0=0.36 for gas diffusion layer with porosity 0.6
基于介觀模擬尺度的格子Boltzmann方法,首次建立了孔隙尺度下燃料電池氣體擴(kuò)散層內(nèi)液態(tài)水凝固模型,對(duì)孔隙率分別為0.5、0.6、0.7、0.8和0.9的二維200 μm×200 μm氣體擴(kuò)散層內(nèi)微孔隙結(jié)構(gòu)中液態(tài)水結(jié)冰過(guò)程進(jìn)行數(shù)值模擬研究.主要成果如下:
(1) 引入格子Boltzmann方法到燃料電池用氣體擴(kuò)撒層內(nèi)液態(tài)水結(jié)冰問(wèn)題研究中,并且通過(guò)一維半無(wú)限大空間凝固熱傳導(dǎo)、二維直角區(qū)域凝固和二維介質(zhì)方腔凝固三組數(shù)值試驗(yàn)驗(yàn)證了本文模型可以精確處理碳纖維、液態(tài)水和冰三種物質(zhì)不同熱物理參數(shù)的問(wèn)題.
(2) 對(duì)應(yīng)孔隙率分別為0.5、0.6、0.7、0.8和0.9時(shí),氣體擴(kuò)散層孔中液態(tài)水完全結(jié)冰所用量綱一時(shí)間F0分別為2.67、3.11、3.68、4.31和4.84.氣體擴(kuò)散層孔隙率越大,完全結(jié)冰所用時(shí)間越多.
(3) 對(duì)應(yīng)孔隙率分別為0.5、0.6、0.7、0.8和0.9時(shí),有自然對(duì)流情況下的結(jié)冰時(shí)間F0比無(wú)自然對(duì)流時(shí)分別減少0、0、0.001、0.001和0.007.在燃料電池氣體擴(kuò)散層微孔隙內(nèi)液態(tài)水結(jié)冰過(guò)程中,自然對(duì)流現(xiàn)象對(duì)結(jié)冰時(shí)間影響較小.