doc. dr. Jaka Dujc • MODELIRANJE GORIVNIH CELIC
MODELIRANJE GORIVNIH CELIC
MODELLING OF PEM FUEL CELLS
doc. dr. Jaka Dujc, univ. dipl. inž. grad.	Znanstveni članek
jdujc@fgg.uni-lj.si	UDK 519.6:662.9 Univerza v Ljubljani, Fakulteta za gradbeništvo in geodezijo, Jamova 2, Ljubljana
Povzetek l Gorivne celice PEM bodo v prihodnjih desetletjih igrale vodilno vlogo pri prizadevanju za čistejše tehnologije v mobilnosti in stacionarnih aplikacijah. Modeliranje gorivnih celic PEM na več nivojih omogoča vpogled v zapletene procese in služi tudi kot orodje za optimizacijo komponent in materialov. Predstavljamo tipične pristope, uporabljene pri simulacijah in analizah mikrostruktur, MEA-sklada ter manjših gorivnih celic, podrobneje pa obravnavamo 2+1D-pristop za analizo velikih celic. Pri tem pristopu je gorivna celica razdeljena na 2-dimenzionalne in 1-dimenzionalne domene. Dve vzporedni 2D-do-meni predstavljata kanale bipolarnih plošč ter sosednji porozni GDL na anodni in katodni strani. 2D-domeni sta povezani z lD-domenami, ki predstavljajo MEA-sklad.
Ključne besede: gorivne celice PEM, mobilnost, numerični modeli, simulacije, metoda končnih elementov, mobilnost, povezani problemi
Summary l The main topics of this article are numerical models of PEM fuel cells. It is estimated that PEM fuel cells will play a leading role in the pursuit of cleaner technologies in mobility and in stationary applications. Multi-scale modeling of PEM fuel cells gives us an insight into complex physical and electro-chemical processes and also serves as a tool for optimizing components and materials. We present the typical approaches used in the simulations and analyses of microstructures, membrane-electrode-assembly and smaller-size fuel cells. In more detail we consider the 2 + 1D approach for the analysis of large fuel cells. In this approach, the fuel cell is divided into regions of 2D and 1D domains. Two parallel 2D domains represent the anode and cathode side flow fields and gas diffusion layers. Anode and cathode sides are connected by the 1D domains, representing the membrane electrode assembly.
Key words: polymer electrolyte membrane fuel cells, mobility, numerical models, simulation, finite element method, coupled problems
1*UVOD
Pred nedavnim so nas delegati 195 držav, ki sodelujejo v Medvladnem forumu o podnebnih spremembah [IPCC, 2018], ponovno opozorili na vpliv človeka, ki ga ima na okolje, in potrebo po zmanjšanju izpustov toplogred-nih plinov. V zadnjih letih so se na področju transporta že zgodile spremembe, ki kažejo premik k uporabi čistejših virov energije. Predvsem uporaba električnih baterijskih avtomobilov je postala nekaj vsakdanjega. Je pa število električnih vozil v primerjavi z uveljavljenimi vozili na fosilna goriva še vedno zelo majhno.
Poleg baterijskih vozil so širši javnosti na voljo manj znana električna vozila na gorivne celice. Gorivne celice PEM (ang. proton exchange membrane fuel cell) so eden glavnih kandidatov za veliko vlogo tako v transportu (slika 1) kot tudi v stacionarnih aplikacijah prihodnjih desetletij. Pri gorivnih celicah PEM se kemična energija vodika in kisika preko kemičnih reakcij pretvori v električno in toplotno energijo ter vodo. Glavna razlika med PEM-tehnologijo in baterijami je v tem, da je pri slednjih vsa energija shranjena v baterijah, medtem ko je pri PEM-tehnologiji energija
shranjena v rezervoarju vodika, ki se nato z uporabo gorivnih celic pretvori v električno energijo. Uporaba rezervoarja namesto baterij prinese v določenih primerih veliko prednost. Pri prometu se prednost PEM-tehnologije pokaže predvsem v dosegu vozil in času polnjenja. V tem pogledu so PEM-vozila povsem enakovredna vozilom na fosilna goriva.
Poleg cestnega prometa pa so gorivne celice PEM zanimive tudi pri drugih aplikacijah (za pregled glej npr. [FCS, 2017]). V septembru 2018 je na severu Nemčije začel voziti prvi potniški vlak na vodik, ki na redni liniji zamenjuje stare dizelske vlake [Guardian, 2018]. V porastu je uporaba sistemov z gorivnimi celicami na mestih, kjer je nujna neprekinjena oskrba z električno energijo. V vlogi pomožnih
268
MODELIRANJE GORIVNIH CELIC • doc. dr. Jaka Dujc
generatorjev se PEM-sistemi vgrajujejo v bolnišnice, podatkovne centre, telekomunikacijske stolpe itd. (npr. [DOE, 2017]).
Za gradbenike oziroma za gradbišča je zanimiva uporaba gorivnih celic kot primarnega vira energije. PEM-tehnologija je namreč idealna za odmaknjena področja oziroma za območja, kjer električno omrežje še ni vzpostavljeno (npr. [DOE, 2017]). Še posebej pa je PEM-teh-
nologija privlačna, če se izkoristi tudi proizvedena toplota, ki je sicer zavržena pri delovanju gorivnih celic. Pri takih sistemih se gorivne celice uporabijo za kogeneracijo električne in toplotne energije v stavbah [Ham, 2015]. Ne nazadnje lahko gorivne celice poganjajo tudi gradbeno mehanizacijo [RR, 2015].
Ko govorimo o prednostih in možnostih uporabe gorivnih celic PEM, je treba omeniti tudi
Slika 1« Osebna in tovorna vozila na gorivne celice: a) Toyota mirai (avtor slike Michal Setlak*), b) honda clarity (avtor slike Alexander Migl*), c) hyundai ix35 (avtor slike Spielvogel*), d) toyota project portal; *(licenca CC BY-SA 4.0 (https://creativecommons.org/licenses/by-sa/4.o), Wikimedia Commons).
nekaj pomanjkljivosti. Prva pomanjkljivost gre na račun (ne)obstoječe infrastrukutre za polnjenje z vodikom, a se tudi na tem področju v Sloveniji že kažejo določeni premiki [MZI, 2018]. Preostale prepreke, ki jih bo treba v prihodnjih letih premostiti, so povezane predvsem s povečanjem zmogljivosti sistema in znižanjem cene množične proizvodnje. Zaradi tega so razvoj in raziskave na tem področu v polnem teku in vedno več je potreb po natančnih numeričnih modelih gorivnih celic. Glavne ovire pri modeliranju gorivnih celic PEM so: razmerje med površino gorivne celice in debelino poroznih komponent (zaradi česar je 3D-analiza celotne celice praktično nemogoča) in številni povezani fizikalni in elektrokemični procesi, ki so vključeni v delujočo gorivno celico (zaradi česar so modeli zelo zapleteni in numerično nestabilni).
V tem delu v nadaljevanju predstavljamo pristope, ki jih uporabljamo pri preskakovanju prej omenjenih preprek pri modeliranju. V razdelku 2 predstavljamo osnovne koncepte delovanja ter sestavne dele gorivne celice. V razdelku 3 so prikazani pristopi za modeliranje od nivoja mikrostrokture pa vse do manjše gorivne celice. Največji poudarek v članku pa je na predstavitvi numeričnega 2+1D-modela za analizo velikih gorivnih celic. Ta model je predstavljen v razdelku 4, kjer poleg formulacije prikazujemo tudi nekaj rezultatov numeričnih simulacij. V razdelku 5 podajamo še nekaj zaključnih besed.
2*KAJ SO GORIVNE CELICE PEM IN KAKO DELUJEJO?
Gorivne celice PEM so naprave, ki spreminjajo kemično energijo goriva v električno energijo
Slika 2« Shema elektrokemičnih procesov pri delovanju gorivne celice.
s pomočjo kemične reakcije in oksidanta. Pri gorivnih celicah PEM se za gorivo uporablja vodik, za elektrolit pa se uporablja posebna polimerna membrana, ki je nepropustna za elektrone, dopušča pa prehod vodikovih ionov (protonov). Za oksidant se lahko uporablja čisti kisik ali pa kisik iz zraka.
2.1 Elektrokemične reakcije
Na sliki 2 je prikazan princip delovanja gorivne celice PEM. Na sredini gorivne celice je elektrolitska membrana, levo in desno od nje pa so porozne in električno prevodne plasti anode in katode. Električni krog je sklenjen s povezavo anode in katode preko porabnika (npr. elektromotorja).
Delovanje gorivne celice je razdeljeno na dva procesa, in sicer na oksidacijo na anod-ni strani in na redukcijo na katodni strani.
Pri oksidaciji na anodni strani vodik odda elektrone
H2^2H+ + 2e~,	(1)
vodikovi ioni (protoni) nato prehajajo preko elektrolitske membrane na katodno stran, elektroni pa pridejo na katodno stran preko porabnika. Na katodni strani nastane redukcija
(2)
~02+ 2H+ + 2e~ H20, kjer se elektroni, vodikovi ioni ter kisik združijo v vodo. Celotno kemično reakcijo delovanja gorivne celice lahko zapišemo kot
H20 + toplota.	(3)
Za zanimivost naj omenimo, da je delovanje gorivne celice PEM ravno obraten proces od pridobivanja vodika z elektrolizo vode in tudi sestava naprave za elektrolizo je enaka gorivni celici PEM.
2.2 Sestavni deli
Kot je bilo že omenjeno, sestavljajo gorivno celico PEM anoda, katoda in elektrolitska
Gradbeni vestnik • letnik 67 • december 2018
269
doc. dr. Jaka Dujc • MODELIRANJE GORIVNIH CELIC
Slika 3^ Prikaz delovanja gorivne celice s tipičnimi sestavnimi deli.
Slika MEA-sklad s prikazom glavnih procesov (avtor izvirne slike CFA213FCE) [licenca CC BY-SA 3.0 (https://creativecommons.org/licenses/by-sa/4.0), Wikimedia Commons].
3«OD MIKROSTRUKTURE DO MODELA MANJŠE GORIVNE CELICE
membrana, kjer so posamezna področja sestavljena iz več slojev. Tipična konfiguracija gorivne celice PEM sestoji iz 7 slojev (glej sliko 3); če začnemo z anodne strani:
•	anodna bipolarna plošča, po kateri se dovaja vodik in se prevaja električni tok,
•	porozen sloj GDL (ang. gas diffusion layer), preko katerega prehaja vodik, hkrati pa omogoča prevajanje električnega toka,
•	plast platinastega katalizatorja CL (ang. catalyst layer), ki omogoči vodiku, da odda elektrone,
•	elektrolitska membrana PEM, ki preprečuje transport elektronov, omogoča pa prehod vodikovih ionov,
•	plast platinastega katalizatorja CL, ki omogoči, da se elektroni, vodikovi ioni ter kisik združijo v vodo,
•	GDL, preko katerega prehaja kisik in se odstranjuje proizvedena voda, omogoča pa tudi prevajanje električnega toka,
•	katodna bipolarna plošča, po kateri se dovaja kisik, odvaja proizvedena voda in se prevaja električni tok.
V skupnosti, ki se ukvarja z raziskavami gorivnih celic, se je za sestav iz srednjih petih slojev (GDL + CL + PEM + CL + GDL), kjer potekajo glavni procesi (glej sliko 4), uveljavilo ime MEA-sklad (ang. membrane electrode assembly). Tako je pri posamezni gorivni celici MEA-sklad stisnjen med dve bipolarni plošči.
Modeliranje gorivnih celic je zahtevno, saj so kompleksne naprave, v katerih potekajo različni fizikalni/kemični procesi, vključenih je več različnih materialov, dodatno težavo pa povzroča dejstvo, da imamo opravka z dogajanjem na različnih skalah (slika 5). Na najmanjši skali, ko govorimo npr. o mikrostruk-turi GDL, imamo opravka s strukturami (vlakni in porami), ki so večinoma manjše od 10u m. Če bi si želeli podrobneje pogledati dogajanje v katalizatorski plasti ali pa v membrani, bi morali obravnavati še manjše strukture. Na drugi strani spektra, če se omejimo zgolj na delovanje gorivnih celic oz. celičnih sestavov (ang. cell stack), pa imamo opravka z dimenzijami
nekaj 10 cm. Narava gorivnih celic pa je taka, da ima dogajanje na najmanjšem nivoju zelo velik vpliv na dogajanje celotnega sistema.
Razkorak med temeljnim razumevanjem, povezanimi s procesi na ravni mikrostrukture in modelom celotnega sistema, ki bi s primerno zanesljivostjo opisal delovanje gorivnih celic, je velik. Vendar pa postane naša naloga obvladljiva, če uporabimo večplasten pristop (ang. multiscale approach) in premostimo vrzel v več segmentih. Rezultati z manjših skal se tipično uporabljajo za parametrizacijo naslednje skale, rezultati z večjih skal pa se uporabljajo za izboljšanje procesov na krajših skalah in za optimiranje delovanja celice.
V naslednjih razdelkih si bomo pogledali, kakšne analize in kakšni načini modeliranja se uporabljajo na posameznih skalah.
3.1 MIKRONIVO
Numerične simulacije, eksperimenti ter študije na mikronivoju so bistveni za razumevanje temeljnega dogajanja v gorivnih celicah. Z njimi lahko določamo povezave med karakteristikami mikrostrukture (npr. geodetska zvitost (ang. geodesic tortuosity), parameter zožitve (ang. constrictivity) ter koeficient poroznosti; [stenzel, 2016]) in efektivnimi transportnimi lastnostmi poroznih materialov (koeficient difuzije, prepustnost in konduktivnost ([Holzer, 2017a], [Math2Market, 2016]). Pri eksperimentalnem delu se veliko podatkov zajame z metodami slikanja (npr. rentgenska tomografija, FIB-tomografija, nevtronska radiografija),
270
MODELIRANJE GORIVNIH CELIC • doc. dr. Jaka Dujc
Slika 5* Obravnavanje gorivnih celic na več skalah in karakteristične dolžine za posamezni nivo.
dobljeni podatki pa se nadalje obravnavajo plini, upoštevani kot robni pogoji. Idealizirana z metodami mikrostrukturne analize [Holzer, slika MEA-modela je prikazana na sliki 7.
Slika 5* Simulacija vdora vode Monte Carlo na mikronovoju. Vlakna so obarvana rjavo, voda pa modro. Simuliran je vzorec velikosti 1126 m u 1126 m u 70 u m, [Capone, 2018].
20l7b) ali pa se na pridobljeni geometriji izvajajo numerične simulacije. Rezultat ene takih simulacij je prikazan na sliki 6, kjer se je z uporabo metode Monte Carlo simuliral vdor tekoče vode v GDL [Capone, 2018].
3.2 MEA-NIVO
Z rezultati in opažanji iz mikroskale lahko gradimo numerične modele na večjih skalah. Eden takih je stacionarni enodimenzionalni makrohomogni model MEA-sklada, katerega detajli so predstavljeni v [Vetter, 2018]. V tem modelu je MEA upoštevana s petimi zaporednimi domenami, medtem ko so vplivi bipolarnih plošč oz. kanalov, po katerih se dovajajo
Model je opisan z osmimi parcialnimi diferencialnimi enačbami (PDE), s katerimi opišemo transport naboja elektronov in protonov, kon-dukcijo toplote, transport kemično vezane vode, difuzijo vodne pare, difuzijo vodika, difuzijo kisika ter transport tekoče vode. V tabeli 1 so povzeti fizikalni procesi z oznakami njihovih neznank, tokov ter domen, kjer je posamezen proces prisoten. Fizikalne procese opišemo s kontinuitetnimi enačbami
Vja=Qa> oc = {(pe,(pp,T,:i,c„,cH2,c02,cl}, (4)
kjer je ja tok za posamezen proces, Qa predstavlja izvor/ponor, a pa je oznaka za proces. Vsi procesi, razen transporta kemično vezane vode, ki je odvisen tudi od gostote toka protonov, imajo naslednjo obliko gostote toka ja=-KaVa,	(5)
kjer je Ka splošna oznaka transportnih koeficientov. Enačbe (4) in (5) predstavjajo diferencialne enačbe drugega reda, s katerimi se v gradbeništvu srečamo npr. pri modelu paličja ali pa pri enodimenzionalnem modelu kondukcije toplote. Oblika diferencialnih enačb torej ni zelo kompleksna. Model postane kompleksen zaradi vrste povezav med fizikalnimi procesi, in sicer so procesi med seboj povezani tako na nivoju transportnih koeficientov, ki so v splošnem odvisni od več fizikalnih polj ( Ka = f(tp. , T,X, cv, cHi, c0i, c,) ), kakor tudi preko izvorov oz. ponorov (Qa = f (Ve ,vP,T, cv, ch 2, co2, ci)). Ko govorimo o kompleksnosti, je treba omeniti tudi pet različnih domen, v katerih potekajo različni fizikalni procesi in so sestavljene iz različnih materialov.
S takim MEA-modelom [Vetter, 2018] lahko zelo podrobno opišemo in študiramo povezave med fizikalnimi in elektrokemičnimi procesi, zavedati pa se je treba, da nam da tak model le informacijo o »povprečnem« odzivu gorivne celice, saj geometrija bipolarnih plošč
Slika 7* Idealizirana geometrija MEA-modela s petimi sloji (ni v merilu). Različni fizikalni procesi v posameznih slojih so označeni z barvami.
271
doc. dr. Jaka Dujc • MODELIRANJE GORIVNIH CELIC
Ime	Neznanka	Tok	Izvor	GDL-A	CL-A	PEM	CL-K	GDL-K
Naboj elektronov	%e		Q%c	x	x		x	x
Naboj protonov	%	j p	Q%		x	x	x	
Temperatura	T	j	Qt	x	x	x	x	x
Kemično vezana voda	X	X	Qx		x	x	x	
Vodna para	Cv	j	Qcv	x	x		x	x
Vodik	CH2	J c„ 2	QCH2	x	x			
Kisik	S	jc0 02	QC02				x	x
Tekoča voda H2O	Cl	j	QCl				x	x
Tabela 1* Povzetek fizikalnih procesov, njihovih prisoten.
s kanali ni eksplicitno upoštevana. Vsekakor pa je tak model dobrodošel za optimizacijo mikrostrukture oz. kot izhodišče za razširitev na dve oziroma tri dimenzije. Na tem mestu je treba omeniti še simulacijska orodja, v katerih implementiramo MEA-mo-dele. Izbira orodja je v veliki meri odvisna od namena uporabe. In sicer se po navadi v prvem koraku, kadar numerični model šele gradimo, naslonimo na komercialna orodja, kot je na primer COMSOL [Comsol, 2017], ki nam omogočijo, da hitro vključimo enačbe, ni se nam pa treba obremenjevati z numeričnimi podrobnostmi in metodami reševanja. V večini primerov pa je okvir komercialnih programov premajhen in z njim ne moremo zajeti vsega, kar naši modeli potrebujejo. Takrat se poslužujemo orodij, kot so Matlab [MathWorks, 2017] (uporabljen v [Vetter, 2018]), Mathematica [Wolfram, 2018] ali pa AceGen in AceFEM [Korelc, 2016], ki znotraj Mathematice omogočata hitro produkcijo in testiranje numeričnih modelov.
3.3 Makrohomogeni modeli na nivoju majhne gorivne celice
Z razširitvijo v prejšnem razdelku prikazanega modela na tri dimenzije lahko upoštevamo tudi vpliv geometrije bipolarnih plošč na delovanje gorivne celice. Smo pa pri takem
neznank, tokov in domen, kjer je posamezen proces
dejanju zelo omejeni z velikostjo območja, ki smo ga s takim modelom sposobni modelirati. Problematične so namreč zelo majhne dimenzije plasti v MEA, ki naredijo podrobno tridimenzionalno analizo na večji skali praktično nemogočo. Tako se pri analizi večjih dimenzij pogosto poslužujemo ali redukciji geometrije, kjer posamezne dele gorivne celice upoštevamo kot robne pogoje, ali pa se omejimo pri izbiri fizikalnih procesov, ki jih vključimo v model.
Eden takšnih tridimenzionalnih modelov je bil predstavljen v [Dujc, 2018], kjer je bil glavni poudarek na študiji vpliva različnih hidrofobnih lastnosti materiala na delovanje gorivne celice. V modelu in pri eksperimentalnih raziskavah [Forner-Cuenca, 2015] je bil namreč uporabljen GDL-material, ki je bil sestavljen iz hidrofobnih in hidrofilnih pasov.
Geometrija modela je omejena na GDL na katodni strani, preostale plasti pa so zajete z nelinearnim elektrokemičnim robnim pogojem na strani katalizatorja ter z robnimi pogoji na strani bipolarne plošče. Poleg nelinearnega robnega pogoja pa model upošteva tudi transport vode v tekočem in plinastem stanju z upoštevanjem evaporacije in kondenzacije ter (difuzni in konvektni) transport kisika. Za dodaten fizikalen vpliv je bila upoštevana tudi mehanska kompresija GDL-strukture ter vpliv strukture na transportne parametre. V [Dujc, 2018] so upoštevali, da gre pri kompresiji GDL vsa sprememba volumna na račun zapiranja praznin (por), in so tako povezali volumske deformacije GDL-a
£y=tr{e\	(6)
s prostorsko porazdeljeno vrednostjo efektivne poroznosti
""'u^'	(7)
kjer je s tenzor deformacij, e0 pa poroznost ne-deformiranega materiala. Efektivna poroznost £eff je v numeričnem modelu pomembna, saj skupaj s stopnjo zasičenosti por s tekočo vodo določa sposobnost difuzivnega transporta plinov skozi porozni material. Bolj ko je material stisnjen in bolj ko je zasičen z vodo, slabša je difuzija kisika, ki je potreben za kemično reakcijo.
Na sliki 8 so predstavljeni rezultati numerične simulacije. Na levi strani slike 8 lahko opazimo, da je vrednost normaliziranega difuziv-nega koeficienta najvišja na mestu kanalov, in sicer v območju hidrofobnega materiala, v katerem je manjša količina vode. Najnižja vrednost je pričakovano pod rebri bipolarne plošče na mestu, kjer je material hidrofilen. Pod rebri je material najbolj stisnjen, prav tako je v hidrofilnih območjih zasičenost z vodo največja. Podoben učinek je viden tudi na desni strani slike 8, kjer je prikazana gostota električnega toka. Ta je najvišja na hidrofob-nih mestih v kanalih in najnižja na hidrofilnih mestih pod rebri.
Normaliziran koeficient difuzije [1]
Gostota električnega toka [A/cm2]
a)
Z,
K b)
Slika 8* Simulacija materiala z različnimi hidrofobnimi lastnostmi: a) normaliziran koeficient difuzije za kisik ter b) porazdelitev gostote električnega toka pri napetosti 0.6 V.
272
letnik 67 • december 2018
MODELIRANJE GORIVNIH CELIC • doc. dr. Jaka Dujc
4* MODEL ZA ANALIZO VELIKIH GORIVNIH CELIC, ANALOGEN GRADBENIŠKEMU MODELU PLOŠČ IN STEBROV
Ko imamo opravka z analizo večje gorivne celice ali pa z analizo nekaj gorivnih celic, povezanih v sestav (ang. cell bundle), detajlna trodimenzionalna analiza ne pride v poštev. Največkrat se na tem mestu odločimo, da zreduciramo dimenzije simulacije toka plinov v plinskih kanalih, izvedemo manj zapleten model, ki opisuje procese v MEA, ali pa kombinacijo obeh. Primeri takšnih modelov so 3+0D-model, ki upošteva simulacijo pretoka plinov v treh dimenzijah z ničdimenzionalno (brez mreženja) povezavo med anodo in katodo, ali pa v nadaljevanju predstavljeni 2+1D-model, pri katerem je pretok plina simuliran v 2D, upoštevana pa je 1D MEA-po-vezava med anodno in katodno stranjo (glej sliko 9).
Osnovni princip 2+1D-modela, ki je bil prvič predstavljen v [Schumacher, 2012], je analogen numeričnemu modelu plošč in stebrov. Tridimenzionalni problem razdelimo na dve vzporedni dvodimenzionalni območji (»plošči«), ki ju povežemo z več enodimenzionalnimi območji (»stebri«). V primeru gorivnih celic predstavljata vzporedni »plošči« dogajanje v ravnini gorivne celice, in sicer se osredotočimo na procese v bipolarni plošči in pripadajočih kanalih ter na plast GDL. Eno dvodimenzionalno območje predstavlja dogajanje na anodni stani, medtem ko drugo dvodimenzionalno območje predstavlja katodno stran. lD-povezave med obema območjema niso nič drugega kot numerični model MEA, ki je bil na kratko že predstavljen (glej 3.2 ter [Vetter, 2018]). Primer mreže za tako 2+1D-analizo je prikazan na desni strani slike 9.
Omeniti je treba razširitev geometrije v 1D-kakor tudi v 2D-domenah, ki v prvotni študiji niso vsebovale GDL-slojev, uporabo nestis-ljivega Navier-Stokesovega toka tekočin, ki nadomešča potencialni tok, nadgradnjo izoter-mnega modela z vključitvijo temperaturnega polja po celotnem območju ter vključitev dodatnih fizikalnih procesov v lD-regijah, ki so bile predhodno opisane le s transporti elektronov, protonov in kemično vezane vode.
V	nadaljevanju predstavljeni numerični model sloni na simulacijskem orodju SESES [Seses, 2018]. Seses je ogrodje, s katerim se krmili celotna simulacija in v katerem se izvajajo simulacije 2D-domen. Odziv lD-domen pa se določi s klicem zunanjih samostoječih rutin, ki so zapisane v programskem jeziku C.
4.1 Modeliranje 2D-odziva - »plošče«
V	2D-domenah, s katerimi obravnavamo procese v kanalih bipolarnih plošč ter v porozni GDL-strukturi, modeliramo konveksne in difuzne tokove. Za regije kanalov predvidimo nestisljivi Navier-Stokesov tok
-/?vVv = //VVv-VP-
llH
(8)
kjer je p gostota plina, v=[v„vy] polje hitrosti, H viskoznost mešanice plinov, P je pritisk, hch pa predstavlja višino kanalov. Enačba (8) je bila izpeljana is splošne 3D-različice Navier-Stokesovih enačb, kjer smo upoštevali, da je z komponenta hitrosti enaka nič (vz=0) ter da imata ostali dve komponenti hitrosti kvadratičen Hagen Poiseuillov profil v smeri komponente. V tem smislu je polje hitrosti povprečna različica komponent hitrosti, če bi jih izračunali s 3D-analizo.
b)
Slika 9« Redukcija 3-dimenzionalnega problema na 2- in 1-dimenzionalne domene: a) osnovna shema ter b) primer mreže končnih elementov z 2D-ploščami in 1D-stebri.
V tem delu predstavljeni 2+1D-model je nadgradnja pristopa iz [Schumacher, 2012].
Med kanali v GDL upoštevamo Darcyjev zakon
M
-VP,
(9)
kjer je KGDL propoustnost GDL. Opozoriti je treba, da smo zaradi redukcije dimenzije izgubili del GDL, ki sicer leži povsod pod kanali in pod rebri bipolarne plošče. Kar se tiče procesov v ravnini gorivne celice, je GDL upoštevan le na mestu reber, hkrati pa je postavljen tudi v isto ravnino kot kanali.
2D-stacionarna oblika kontinuitetne enačbe o
ohranitvi mase ima obliko
V(/7v)=nm,	(10)
kjer nm predstavlja izvor/ponor mase. K enačbam (8)-(10) na robnem območju r pripadajo naslednji robni pogoji
MTJ = pvin, P(T0Ut) = P0Ut,	(11)
kjer je A skupen masni tok, vin je hitrost plinov na vhodu rin v gorivno celico, Pout pa je pritisk pri izhodu rout iz celice.
V	tem modelu gorivne celice kot posamezne komponente plina nastopajo ch2 vodik, co2 kisik ter ch20 vodna para. Omeniti je treba, da ta model ne upošteva vode v tekočem stanju, bi pa to lahko z nekaj spremembami spremenili. Za vsako od komponent plina a = {ch 2, co2> ch 2o } velja naslednja diferencialna enačba s pripadajočim robnim pogojem
v • (pav) = -V • ja + Tla, ja = -DaVca,
ca(Xi„) = cain.,	(12)
kjer jepa gostota posamezne komponente, ja difuzni tok komponente, določen s Fickovim difuznim koeficientom Da, na izvor/ponor komponente, določen z 1D MEA-modelom, caM pa je robna vrednost za koncentracijo.
V	2D-domenah upoštevamo tudi transport naboja pe znotraj GDL in bipolarnih plošč, ki je določen z
V-;t = ne, je = -a?2DV<pe, r.p.naanodi:
VeVtfD ) = 0. r-P"na kat0di : V^fD ) = Uce,l> (13)
kjer je je tok naboja, ne je izvor/ponor naboja, določen z 1D-modelom, a^2D je električna konduktivnost GDL oziroma reber, Uce„ je robna vrednost (napetost gorivne celice), določena na katodni strani, pa označuje srednje črte reber, kjer so predpisani robni pogoji za naboj.
Zadnji od fizikalnih procesov, ki so obravnavani v 2D-regijah, je konvektivni in konduktivni transport toplote, ki je opisan z naslednjo diferencialno enačbo
v • (-k?2Dt+pcpTv) = nT-n3D, (14)
kjer je K™ toplotna konduktivnost GDL oziroma bipolarnih plošč, T je temperatura, Cp je toplotna kapaciteta mešanice plinov, nT je toplotni tok, določen z 1D-MEA-modelom,
Gradbeni vestnik • letnik 67 • december 2018
nso pa dodaten člen, ki upošteva 3D-učinke (izgube oz. hlajenje v smeri pravokotno na ravnino gorivne celice), ki jih ni mogoče zajeti z 2D-analizo.
4.2 Primerjava 2D- in 3D-simulacij toka tekočin
4.4 Povezava med 2D- in lD-regijami
Upoštevana je močna povezava med 2D- in lD-domenami, kjer so vse neznanke problema (prostostne stopnje) 2D anodne strani, povezane z nasprotno 2D katodno domeno preko lD-modela, ki ga opišemo s številnimi ele-
a)
!o.
Comsol 3D brez GDL SeSes 3D brez GDL SeSes 2D brez GDL Empirično brez GDL Comsol 3D z GDL SeSes 3D z GDL SeSes 2DzGDL
b)
0.10
Pretok [1/min]
Slika 10* Primerjava 2D- in 3D-simulacij toka tekočin: a) geometrija testnega primera z dotokom na levi strani in odtokom na desni strani ter b) krivulje padec tlaka-pretok.
Obnašanje 2D-modela toka tekočin, enačbi (8) in (9), je bilo preverjeno na več testnih geometrijah. Na sliki 10 prikazujemo enega od testnih primerov, in sicer so na levi strani prikazani geometrija in parametri analize, na desni strani pa so prikazane krivulje padca tlaka v odvisnosti od pretoka. Teste smo naredili enkrat brez upoštevanja GDL-plasti, s tem smo preverili posamično delovanje enačbe (8), ter enkrat z upoštevanjem GDL-plasti, kjer smo preverili sočasno delovanje enačb (8) in (9). S tem testom obravnavamo le tok tekočine, zato so v simulacijah vsi preostali fizikalni procesi izključeni.
Rezultate, dobljene z 2D-modelom v SESES [Seses, 2018], smo primerjali s 3D-rezultati iz SESES in COMSOL [Comsol, 2017]. Dodatno smo za primer brez upoštevanja GDL za referenčni rezultat vzeli še Darcy-Weissbachov empirični zakon [Grote, 2007]. S slike 10b je razvidno dobro ujemanje med 2D-pristopom in referenčnimi rezultati. Pri testnem primeru se kvaliteta rezultatov z uporabo 2D-pristopa bistveno ni zmanjšala, opazno pa je znatno izboljšanje simulacijskih časov. Če se omejimo le na orodje SESES - s COMSOL so bili računski časi občutno daljši -se je pri analizi brez GDL računski čas zmanjšal s 97 s na 0,225 s, pri analizi z GDL pa se je računski čas zmanjšal s 337 s na 0,48 s.
4.3 Modeliranje lD-odziva - »stebri«
Uporabljeni MEA-model je zelo podoben že predstavljenemu v razdelku 3.2 (glej tudi [Vetter, 2018]). Edina razlika je v tem, da je opisan z le sedmimi namesto osmimi parcialnimi diferencialnimi enačbami (4), saj transport tekoče vode na katodni strani ni vključen.
menti ter neznankami. Vrednosti spremenljivk 2D-domen se uporabljajo kot Dirichletovi robni
pogoji, ki se uporabijo pri reševanju 1D-prob-
.2D,A rp2D,A 2D,A 2D,A
T ' Chj2 ' h2° , na c--- -, c ). Robni tok-
lema (na anodi Ve katodi vlD>K,T2D , ° , h ° , ovi (»reakcije«), ki jih z uporabljenimi Dirichle-
tovimi robnimi pogoji izračunamo z 1D-mod-
f ID, a f id,a f ID,A f
elom (na anodi fe ,Jt ,Jc„ ,a
1d, a
Slika 11* Shematična predstavitev povezav v 2+lD-modelu. Vrednosti spremenljivk v 2D-domenah predstavljajo robne pogoje za lD-simulacijo (na anodi
2D,A T 2D,A 2D,A 2D,A
Ve ,T , Ch , Ch° , na katodi
„-LUK T 2DK 2D,K 2D,K
Ve ,T , , V , Ch2° ). Robni
tokovi, določeni z lD-simulacijami (na
f1D,A, f1D,A, f1D,A, f1D,A
anodi *ve >f >Jc„ >
2
*1D,K
h 2°, na
1D,K ,
f1d,k f1d,k f1d,k fid,k
katodi fe , Jt , L° , Lh° ),
se uporabijo za določitev izvorov/2 ponorov v 2D-simulaciji.
j ID,K J ID,K J ID,K J ID,K
na katodi fe ,JfT ,,jh2°), pa se uporabijo za določitev izvorov/ponorov v 2D-modelu (n^nnn, v enačbah (10), (12), (13), (14)). Za lažje razlikovanje med posameznimi spremenljivkami uporabljamo oznaki 1D in 2D, ki označujeta, kateremu modelu spremenljivka pripada, ter oznaki A in K, ki označujeta anodno oziroma katodno stran. Shematična predstavitev povezave v 2+1D-modelu je prikazana na sliki 11. Pri pretvarjanju med robnimi tokovi in izvori/ ponori moramo paziti, da so slednji izraženi na prostorninsko enoto, zato moramo upoštevati tudi debelino 2D-domen t^. Izvor toplote v 2D-domeni anode je tako
f\d,a JT
U2td'a =
(15)
kjer fT1D,A določa robni toplotni tok na anodno stran, izračunan z 1D-modelom. Na podoben način določimo izvore/ponore tudi za preostale procese.
4.5 Optimizacija sheme reševanja
Naša prva implementacija 2+1D-modela je temeljila na shemi, predstavljeni v prejšnjem razdelku. Čeprav je ta model veliko hitrejši od popolnega 3D-modela, pa še vedno ni dovolj hiter za izvedbo simulacij tehničnih velikosti gorivnih celic. Še vedno je potrebno ogromno število neznank, s katerimi opišemo povezave med anodno in katodo (vsako vozlišče v 2D-domeni potrebuje vsaj 100 dodatnih prostostnih stopenj). Poleg velikosti problema pa dodatno težavo povzročajo nelinearnosti, ki so del 1D-modela. 1D-model potrebuje veliko iteracijskih korakov, da doseže rešitev, v nekaterih primerih pa model tudi divergira, in se poruši celotna shema reševanja 2+1D-modela.
V izognitev omenjenim težavam predlagamo optimizacijo postopka rešitve na naslednji način (glej sliko 12). Namesto istočasnega reševanja 2D- in 1D-simulacij naredimo to zaporedoma. V prvem koraku se ukvarjamo le z 1D-modelom, kjer opravljamo številne simulacije s korakanjem skozi razpon robnih
2D,a rp 2D,a 2D,a 2d,a
pogojev (na anodi Ve	,ch2 ,ch2°,
.. ,^2D,K rp2D,K 2D,K 2D,K v ,.
na katodi Ve ,T ,c°2 ,cH2° ), ki se pričakujejo v 2D-simulaciji. Končni rezultat take parametrične študije je vrsta mejnih tokov s povezanimi vrednostmi robnih pogojev. Ta niz podatkov nato uporabimo za definiranje robnih tokov, ki so potrebni pri simulaciji 2D-domen. Podatke lahko v 2D-analizo vključimo na več načinov. Ena možnost je, da izgradimo veliko bazo podatkov, iz katere med analizo beremo želene vrednosti, ali
274
Slika 12« Optimizacija sheme reševanja in zamenjava implicitnega modela MEA z eksplicitnim modelom.
Slika 13* Krivulje napetost-gostota električnega toka, dobljene z 1D MKE-simulacijo v AceFEM (Korelc, 2016] (polna črta) ter dobljene z nastavkom (16) (nepovezane točke) za
različne mešanice vhodnih plinov. Vrednost ch
/en
določa razmerje parcialnih tlakov
posamezne komponente na anodni ( ch
pa za drugo možnost uporabimo nevronske mreže. Mi smo se odločili za tretjo in najbolj pragmatično možnost, in sicer smo uporabili analitične nastavke za določitev robnih tokov. Pri določitvi oblike nastavkov nam je koristilo poznavanje problema, zato smo se lahko oprli tudi na intuicijo. Seveda pa nikakor ne trdimo, da so izbrani nastavki najboljša možnost v danem primeru. Na slikah 13 in 14 so prikazani rezultati parametrične študije, ki je bila opravljena v Mathematici z AceFEM [Korelc, 2016]. Za gostoto električnega toka i, ki je eden glavnih podatkov pri gorivnih celicah, smo izbrali naslednji nastavek
i =k0+ ky + k2V2 + k,V3 +	+ k5c20D2'K +
+k,Vc%< +ktV2c2HD-A lk9V2c^,(16)

kjer V določa razliko v napetostih med elektrodama (V = (p2eD-K -(p2eD-A), koeficienti k0 do k9 pa so bili v Mathematici določeni z metodo najmanjših kvadratov. S stališča
■ ) oziroma katodni (cO ) strani.
2+1D-simulacij nas zanima območje napetosti med 0,6 V in 0,9 V, zato smo, da bi dobili čim boljše ujemanje, to območje razdelili na 3 podobmočja (0,9V-0,8V, 0,8V-0,7V and 0,7V-0,6V) z različnimi vrednostmi koeficientov. Na sliki 13 so prikazane krivulje napetost-gostota električnega toka, dobljene z 1D MKE-simulacijo v AceFEM [Korelc, 2016] ter dobljene z nastavkom (16) za različne mešanice vhodnih plinov. S slike 13 je v območju fitanja razvidno dobro ujemanje med obema metodama določanja gostote električnega toka.
Izraz (16) je neposredno povezan z robnima
tokovoma naboja elektronov
f£A=-i,	(17)
Podobno pa lahko, ker je gostota električnega toka linearno povezana s porabo kisika in vodika v kemičnih reakcijah, za robna tokova vodika in kisika zapišemo
Jh2 af >JO2 2f>	(l«)
kjer je F Faradayeva konstanta. Za toplotna tokova na anodni in katodni strani (f1D-A in fT1D,K ) ter za tok vodne pare na anodni in katodni strani (-CO in fcD'2O) pa predpostavimo kvadratične polinome, odvisne od gostote električnega toka i. Primerjava rezultatov med kvadratičnimi polinomi ter MKE-rezultati je prikazana na sliki 14. Za dani nabor rezultatov omogočajo kvadratični nastavki zelo dobro
Slika 14* Primerjava med MKE-rezultati (Ace-Gen, AceFEM (Korelc, 2016]) in kvadratičnimi nastavki: a) krivulje toplotni tok-gostota električnega toka, b) krivulje tok vodne pare-gostota električnega toka.
4.6 Rezultati simulacij
Delovanje numeričnega modela, predstavljenega v prejšnjih razdelkih, je bilo preverjeno na gorivni celici, ki ima aktivno območje, veliko približno 30,25 m2; glej sliko 15. Na levi strani slike 15 je predstavljena geometrija ene bipolarne plošče. V bipolarni plošči so trije paralelni serpentinasti kanali z dovodom plina na spodnji levi strani in z odvodom na zgornji desni strani. Posamezna elektroda je bila v 2D simulirana z mrežo končnih elementov, prikazano na desni strani slike 15 - mreža vsebuje 48.400 (220 x 220) kvadratnih končnih elementov. Celoten problem je torej opisan z dvema mrežama po 48.400 končnimi elementi, ki sta v vsakem vozlišču (48.841 vozlišč) med seboj povezani z dodatnimi 1D-elementi. Na slikah 16 in 17 je prikazanih nekaj rezultatov simulacije. Na sliki 16 levo je prikazana porazdelitev koncentracije vodika, medtem ko je na desni strani prikazana porazdelitev kon-
Gradbeni vestnik • letnik 67 • december 2018
doc. dr. Jaka Dujc • MODELIRANJE GORIVNIH CELIC
centracije kisika. Pričakovano je koncentracija obeh plinov največja spodaj levo, kjer se plini dovajajo, ter najmanjša zgoraj desno, kjer se plini odvajajo iz gorivne celice. Koncentracija plinov je tipično večja v kanalih, kjer je upoštevan prost pretok plinov. Na mestih reber so koncentracije manjše, saj se do teh mest plini prenašajo s počasnejšim procesom difuzije. Na sliki 17 levo je prikazana porazdelitev potenciala elektronov. Opaziti je različne vrednosti na mestih reber, kjer je upoštevan robni pogoj, ter v kanalih. So pa razlike v vrednostih zelo majhne, saj sta tako GDL kot bipolarna plošča zelo prevodna za elektrone. Na sliki 17 desno je prikazana še porazdelitev gostote električnega toka. Ta je zelo podobna kot pri koncentracijah plinov, in sicer ima največjo vrednost v kanalih pri vhodu v celice (1,2 A/cm2) in najmanjšo vrednost pod rebri pri izhodu iz celice (0,6 A/cm2).
Slika 15* a) Geometrija gorivne celice s tremi vzporednimi serpentinastimi kanali. b) Mreža končnih elementov (220 x 220), uporabljena za simulacijo ene 2D-domene.
Slika 16* a) Porazdelitev koncentracije vodika na anodi. b) Porazdelitev koncentracije kisika na katodi.
a)

b)
Slika 17* a) Porazdelitev potenciala elektronov na katodni strani. b) Porazdelitev gostote električnega toka.
5-SKLEP
Gorivne celice PEM so na pragu komerciali-zacije v avtomobilski industriji in pri stacionarnih aplikacijah. Na tej stopnji je treba zagotoviti zanesljive modele na vseh skalah, ki se lahko uporabljajo kot orodje za optimizacijo materialov gorivnih celic, da bi dosegli večjo
učinkovitost, večjo energijsko gostoto in večjo zanesljivost. V članku so bili predstavljeni pristopi za analizo mikrostruktur, ki lahko v kombinaciji z makrohomogenimi modeli služijo kot orodje za optimiranje poroznih materialov.
Podrobneje je bil predstavljen 2+1D-model za analizo velikih celic. Model se je izkazal za zelo uporabno orodje, s katerim je mogoče relativno hitro in natančno simulirati velike gorivne celice. Z nekaj spremembami bi ta model lahko pretvorili v model za analizo celičnega sestava. Tak model bi služil kot orodje za študije in načrtovanje celotnega sistema gorivnih celic.
276
6'ZAHVALA
V tem članku predstavljeno delo je bilo opravljeno med avtorjevim delovanjem na Institute of Computational Physics, Zurich Univ. of Appl.
Sciences v Winterthuru, v Švici. Raziskovalno delo so s financiranjem omogočile naslednje institucije: Swiss National Science Founda-
tion (Energy Turnaround, NRP 70, project no. 153790, grant no. 407040_153790), Swiss Commission for Technology and Innovation (Swiss Competence Center for Energy Research - SCCER Mobility, contract no. KTI.2014.0115), Swiss Federal Office of Energy SFOE in Belenos Clean Power Holding.
7*LITERATURA
Capone, L., Marmet, P., Holzer, L., Dujc, J., Schumacher, J., Lamibrac, A., Büchi, F., Becker, J., An ensemble Monte Carlo simulation study of water distribution in porous gas diffusion layers for proton exchange membrane fuel cells, Journal of Electrochemical Energy Conversion and Storage, 15(3), 10 str., 2018.
COMSOL, COMSOL Multiphasics Reference Manual, version 5.3, https://www.comsol.com/, 2017.
DOE, The United States Department of Energy, Fuel Cells for Stationary Power Applications, https://www.energy.gov/sites/prod/files/2018/01/ f46/fcto_fc_stationary_power_apps.pdf, 2017.
Dujc, J., Marmet, P., Schumacher, J.O., Forner-Cuenca, A., Cochet, M., Boillat, P., Modelling the Effects of using Gas Diffusion Layers with Patterned Wettability for Advanced Water Management in Proton Exchange Membrane Fuel Cells, Journal of Electrochemical Energy Conversion and Storage, 15(2), 14 str., 2018.
FCS, The Fuel Cell Store, Introduction to Fuel Cell Applications, http://www.fuelcellstore.com/blog-section/intro-fuel-cell-applications, 2017. Forner-Cuenca, A., Biesdorf, J., Gubler, L., Kristiansen, P.M., Schmidt, T.J., Boillat, P., Engineered Water Highways in Fuel Cells: Radiation Grafting of Gas Diffusion Layers, Adv. Mater., 27(41), 6317-6322, 2015.
Grote, K.H., Feldhusen, J., DUBBEL Taschenbuch für den Maschinenbau, volume 22, Springer, 2007.
Guardian, Germany launches world's first hydrogen-powered train, https://www.theguardian.com/environment/2018/sep/17/germany-launch-es-worlds-first-hydrogen-powered-train, 2018.
Ham, S.W., Jo, S.Y., Dong, H.W., Jeong, J.W., A simplified PEM fuel cell model for building cogeneration applications, Energy and Buildings, 107, 213-225, 2015.
Holzer, L., Pecho, O., Schumacher, J., Marmet, P., Stenzel, O., Büchi, F.N., Lamibrac, A., Münch, B., Microstructure-property relationships in a gas diffusion layer (GDL) for Polymer Electrolyte Fuel Cells, Part I: effect of compression and anisotropy of dry GDL, Electrochim. Acta, 227, 419-434, 2017. Holzer, L., Pecho, O., Schumacher, J., Marmet, P., Büchi, F.N., Lamibrac, A., Münch, B., Microstructure-property relationships in a gas diffusion layer (GDL) for Polymer Electrolyte Fuel Cells, Part II: pressure-induced water injection and liquid permeability, Electrochim. Acta, 241, 414-432, 2017. IPCC, The Intergovernmental Panel on Climate Change, Global Warming of 1.5 °C, https://www.ipcc.ch/report/sr15/, 2018. Korelc, J., Wriggers, P., Automation of Finite Element, Springer, 2016. MathWorks, Inc., MATLAB 9.2 R2017a, Natick, Massachusetts, Združene države Amerike, 2017. Math2Market, GeoDict, http://www.geodict.com, 2016.
MZI, Ministrstvo za infrastrukturo, Akcijski program za alternativna goriva v prometu, http://www.mzi.gov.si/fileadmin/mzi.gov.si/pageuploads/ NOVICE_slike/Akcijski_program_za_alternativna_goriva.doc, 2018.
Randall-Reilly, LLC, JCB just bought a controlling stake in a hydrogen fuel cell company. Hmmm..., https://www.equipmentworld.com/jcb-just-bought-a-controlling-stake-in-a-hydrogen-fuel-cell-company-hmmm/, 2015
SESES, SESES Manual, https://www.zhaw.ch/en/engineering/institutes-centres/icp/multiphysik-modellierung/multiphysics-software-nm-seses/, August 2016.
Stenzel, O., Pecho, O., Holzer, L., Neumann, M., Schmidt, V., Predicting effective conductivities based on geometric microstructure characteristics, AIChE J., 62, 1834-1843, 2016.
Vetter, R., Schumacher, J., Free Open Reference Implementation of a Two- Phase PEM Fuel Cell Model, Computer Physics Communications, 234, 223-234, 2018.
Schumacher, J.O., Eller, J., Sartoris, G., Colinart, T., Seyfang, B., 2+1D modelling of a polymer electrolyte fuel cell with glassy-carbon micro-structures, Mathematical and Computer Modelling of Dynamical Systems, 18(4), 355-377, 2012. Wolfram Research, Inc., Mathematica, Champaign, Illinois, Združene države Amerike, 2018.
277