Electrochemical processes and porous media: mathematical and - - PowerPoint PPT Presentation

Weierstrass Institute for Applied Analysis and Stochastics

Electrochemical processes and porous media: mathematical and numerical modeling

Jürgen Fuhrmann Alfonso Caiazzo, Klaus Gärtner, Hartmut Langmach, Alexander Linke, Hong Zhao

Mohrenstrasse 39 · 10117 Berlin · Germany · Tel. +49 30 20372 0 · www.wias-berlin.de RICAM Workshop · Linz · 2011-10-05

Electrochemistry

Corrosion Electroplating Biological processes (heart beat etc.) Batteries Conversion of chemical energy stored in compounds within the cell into

electrical energy

Different variants (low/high temperature, solid/liquid electrolyte ...) Fuel cells Invented ≈ 1840 by Schönbein, Groves Conversion of chemical energy stored in hydrogen, methanol,

carbohydrates ... into electrical energy

Continuous supply of reactants, removal of products Different variants (low/high temperature, solid/liquid electrolyte ...) Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 2 (46)

Polymer electrolyte fuel cell - membrane electrode assembly (MEA) Void space Polymer electrolyte Porous matrix Catalyst particles

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 3 (46)

Polymer electrolyte fuel cell - membrane electrode assembly (MEA) Anode channel Porous layer Reaction zone Membrane Cathode channel Porous layer Reaction zone Carbon fiber + teflon Carbon fiber + teflon + nafion + catalyst Nafion (proton conducting polymer) Carbon fiber + teflon Carbon fiber + teflon + nafion + catalyst Void space Polymer electrolyte Porous matrix Catalyst particles

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 3 (46)

Polymer electrolyte fuel cell - membrane electrode assembly (MEA) Hydrogen fuel cell (H2-PEMFC) Anode channel Porous layer Reaction zone Membrane Cathode channel Porous layer Reaction zone Void space Polymer electrolyte Porous matrix Catalyst particles H2 2H2 → 4H+ +4e−

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 3 (46)

Polymer electrolyte fuel cell - membrane electrode assembly (MEA) Hydrogen fuel cell (H2-PEMFC) Anode channel Porous layer Reaction zone Membrane Cathode channel Porous layer Reaction zone Void space Polymer electrolyte Porous matrix Catalyst particles H2 2H2 → 4H+ +4e− e− e−

– +

Load H+

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 3 (46)

Polymer electrolyte fuel cell - membrane electrode assembly (MEA) Hydrogen fuel cell (H2-PEMFC) Anode channel Porous layer Reaction zone Membrane Cathode channel Porous layer Reaction zone Void space Polymer electrolyte Porous matrix Catalyst particles H2 2H2 → 4H+ +4e− e− e−

– +

Load H+ O2,N2

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 3 (46)

Polymer electrolyte fuel cell - membrane electrode assembly (MEA) Hydrogen fuel cell (H2-PEMFC) Anode channel Porous layer Reaction zone Membrane Cathode channel Porous layer Reaction zone Void space Polymer electrolyte Porous matrix Catalyst particles H2 2H2 → 4H+ +4e− e− e−

– +

Load H+ O2,N2 4H+ +4e− +O2 → 2H2O

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 3 (46)

Polymer electrolyte fuel cell - membrane electrode assembly (MEA) Hydrogen fuel cell (H2-PEMFC) Anode channel Porous layer Reaction zone Membrane Cathode channel Porous layer Reaction zone Void space Polymer electrolyte Porous matrix Catalyst particles H2 2H2 → 4H+ +4e− e− e−

– +

Load H+ O2,N2 4H+ +4e− +O2 → 2H2O H2O,N2

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 3 (46)

Polymer electrolyte fuel cell - membrane electrode assembly (MEA) Hydrogen fuel cell (H2-PEMFC) Anode channel Porous layer Reaction zone Membrane Cathode channel Porous layer Reaction zone Void space Polymer electrolyte Porous matrix Catalyst particles H2 2H2 → 4H+ +4e− e− e−

– +

Load H+ O2,N2 4H+ +4e− +O2 → 2H2O H2O,N2 Overall reaction: 4H2 +O2 → 2H2O+energy Several MEA “sandwiches” combined into a stack.

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 3 (46)

Polymer electrolyte fuel cell - membrane electrode assembly (MEA) Direct methanol fuel cell (DMFC) Anode channel Porous layer Reaction zone Membrane Cathode channel Porous layer Reaction zone Void space Polymer electrolyte Porous matrix Catalyst particles CH3OH,H2O 2CH3OH +2H2O → 2CO2 +12H+ +12e− CO2 e− e−

– +

Load H+ O2,N2 12H+ +12e− +3O2 → 6H2O H2O,N2 Overall reaction: 2CH3OH +3O2 → 4H2O+2CO2 +energy Several MEA “sandwiches” combined into a stack.

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 3 (46)

Some physical effects

Electrolyte flow (free/porous media) Transport + diffusion of dissolved species Charge transport in electric field Electrostatic potential distribution (multistep) Reactions with electron transfer Intercalation, heat transport, swelling ... Aging, ripening ...

Simplifications: high velocity asymptotics, ideal mixing, lumped reactions General case: ⇒ numerical modeling

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 4 (46)

• J. Newman, K. Thomas-Aleya (2004):Electrochemical systems;
• A. Kulikovsky (2010): Analytical Fuel Cell Modeling

Charge transport Charge transport is considered in electrodes and electrolytes.

Electrodes: Solid (metal, carbon, semiconductor ...) Mostly electronic conductors: charge carriers are electrons (e−) Electrolytes Liquid or solid (aquatic solution, molten salt, polymer membranes) Ionic conductors: electrons are blocked, charge carriers are ions of different

type, e.g. protons (H+).

Mixed conductors show properties of both Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 5 (46)

Charge transport: Nernst Planck Poisson system

Transport of i-th dissolved species (i = 1...n)

due to diffusion, electromigration, advection in dilute solution ⇒ Nernst-Planck equation:

• Ni = −

electromigration

• ziuiFci∇φ

−

diffusion

Di∇ci +

advection

• ci

v ∂tci +∇· Ni = ji

Distribution of charged species

⇒ self-consistent electric field ⇒ Poisson equation: −∇·ε∇φ = F

n

∑

i=1

zici Variables n number of species φ electrostatic potential ci species concentration

• Ni

molar flux zi charge ui mobility Di diffusion coefficient ε electrostatic permeability F Faraday constant ji Reaction

• v

Substrate velocity

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 6 (46)

Assumptions behind Nernst-Planck Dilute solution theory:

Ignore interactions (collisions) between different dissolved species ci.

Otherwise: Stefan-Maxwell terms, “concentrated solution theory”

Velocity field not influenced by moving ions.

Otherwise: contribution to momentum balance

Fluid density not influenced by concentration changes

Otherwise: variable density flow Special case: semiconductor device equations.

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 7 (46)

Bulk electroneutrality

ε << F ⇒ electroneutrality in bulk (away from interfaces and boundaries) n

∑

i=1

zici = 0

Sum up equations multiplied by zi:

∂t

• n

∑

i=1

zici

• −∇·
• n

∑

i=1

z2

i uiFci∇φ + n

∑

i=1

ziDi∇ci + v

n

∑

i=1

zici

• =

n

∑

i=1

zi ji

Express e.g. c1:

z1D1c1 = −

n

∑

i=2

D1zici

No bulk reactions ⇒

∇·

• n

∑

i=1

Fz2

i uici∇φ + n

∑

i=2

zi(Di −D1)∇ci

• = 0

Equal diffusion coefficients or small concentration gradients ⇒

Ohm’s law: ∇·κ∇φ = 0

• κ = F

n

∑

i=1

z2

i uici : conductivity

• Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 8 (46)

Nernst Planck Ohm system

Transport of i-th dissolved species (i = 2...n)

due to diffusion, electromigration, advection in dilute solution ⇒ Nernst-Planck equation:

• Ni = −

electromigration

• ziuiFci∇φ

−

diffusion

Di∇ci +

advection

• ci

v ∂tci +∇· Ni = ji

Self-consistent electric field from Ohm’s Law:

∇·κ∇φ = 0 Variables n number of species φ electrostatic potential ci species concentration

• Ni

molar flux zi charge ui mobility Di diffusion coefficient κ conductivity F Faraday constant ji Reaction

• v

Fluid velocity

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 9 (46)

Special case: solid electrolyte/electrode c1: mobile charge carriers c2: immobile charges in solid lattice (u2 = D2 = 0)

• v = 0 ⇒ c2 = const

Electroneutrality,z1 = −z2 ⇒ c1 = c2 small ∇c2 ⇒ κ = z2

1u1Fc1 Electrodes (graphite, metal): c1 ↔ free e− Polymer electrolytes in fuel cells: c1 ↔ free H+. Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 10 (46)

Inert electrode-electrolyte interface

• +

+ + + + + + + + + + + + +

• +
• +
• +
• +
• +
• +
• +

Electrode Electrons Fixed charges at lattice Electrolyte Free ions

Electrodes: electron concentration as c1 ⇒ “electron potential” φs (Acidic) electrolytes: concentration of protons (H+) as c1 ⇒ “proton potential” φl Local electroneutrality violated in boundary layer Potential jump at interface Electrons and protons attracting each other may accumulate on both sides of

the interface creating a double layer

Simple model: double layer charge Qdl = Cdl(φs −φl) Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 11 (46)

Charge transfer reactions at Interfaces

• +

+ + + + + + + + + + + + +

• +
• +
• +
• +
• +
• +
• +

Electrode Electrons Fixed charges at lattice Electrolyte Free ions Interface reaction with electron transfer

Reaction rates in mass action law depend on potential difference n

∑

i=1

aiAi ⇋

n

∑

i=1

biBi +ne− +nH+ r = k+e

nα(φs−φl )F RT

n

∏

i=1

[Ai]ai −k−e

(1−α)(φs−φl )F RT

n

∏

i=1

[Bi]bi (φs: electron potential in electrode, φl: proton potential in electrolyte)

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 12 (46)

Catalytic interface reactions

Most reactions need initial supply of energy in order to be started. With a cataylist, this inital energy barrier may be much lower. Adsorption to electrode surface modeled as reaction with free catalyst sites (electrochemical) Reaction of adsorbates with lower activation energy Desorption from electrode surface Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 13 (46)

Hydrogen Oxidation Reaction (HOR) Anodic reaction in H2-PEMFC s – catalyst site, ·ad – adsorbed species. H2 +2s ⇋ 2Had, r1 = k+

1

cH2 cre f

H2

¯ θ2 −k−

1 θ2 H

Had ⇋ H+ +e− +s, r2 = k+

2 eαδφF/RT θH −k− 2 e(α−1)δφF/RT aH+ ¯

θ cH2: hydrogen concentration θH: Share of catalyst sites occupied by hydrogen ¯ θ = 1−θH: Share of free catalyst sites

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 14 (46)

Hydrogen Oxidation Reaction (HOR) Anodic reaction in H2-PEMFC s – catalyst site, ·ad – adsorbed species. H2 +2s ⇋ 2Had, r1 = k+

1

cH2 cre f

H2

¯ θ2 −k−

1 θ2 H

Had ⇋ H+ +e− +s, r2 = k+

2 eαδφF/RT θH −k− 2 e(α−1)δφF/RT aH+ ¯

θ cH2: hydrogen concentration θH: Share of catalyst sites occupied by hydrogen ¯ θ = 1−θH: Share of free catalyst sites Typical way to include these processes into global model (for given δφ).

• NH2 = D∇cH2 +cH2

v advection-diffusion in Ω ∂tcH2 +∇· NH2 = 0 continuity in Ω

• NH2 ·

n+r1 = 0 normal flux equals reaction on Γ ⊂ ∂Ω ∂tθH −r1 +r2 = 0 evolution of catalyst coverage on Γ ⊂ ∂Ω

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 14 (46)

Oxygen Reduction Reaction (ORR) Multistep cathodic reaction in H2-PEMFC and DMFC Oxygen reaction is split into several steps: O2 +s

k±

1

⇋ O2,ad O2,ad +H+ +e−

k±

2

⇋ HO2,ad HO2,ad +H+ +e−

k±

3

⇋ H2O2,ad H2O2,ad +2H+ +2e−

k±

4

⇋ 2H2O+s

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 15 (46)

Methanol Oxidation Reaction (MOR) Anodic reaction in DMFC CH3OH +2s1

k±

1

⇋ (CH2 −OH)ad +Had (CH2 −OH)ad +s1

k±

2

⇋ (CH −OH)ad +Had (CH −OH)ad +s1

k±

3

⇋ (C −OH)ad +Had (C −OH)ad +s1

k±

4

⇋ (C −O)ad +Had (C −O)ad +OHad

k±

5

⇋ (COOH)ad +s2 (COOH)ad +OHad

k±

6

⇋ CO2 +H2O+s1 +s2 Had

k±

7

⇋ H+ +s1 +e− H2O+s2

k±

8

⇋ OHad +H+ +e−

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 16 (46)

Open questions in catalysis modeling Essentially, the way of modeling these reactions is guided by heuristics. The proposed reaction chains follow one particuala hypothesis about the pathways. Some missing effects:

Catalyst surface defects are preferred adsorption sites Different crystal directions have different reaction speeds Multisite adsorption processes Restructuring of catalyst surface due to reaction Elementary processes at electrodes may be more complex Pt catalyst dissolution Ripening of catalyst particles

More problems

“Sticky” intermediates (e.g. CO) block catalyst sites Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 17 (46)

Porous electrodes Current proportional to interface area ⇒ increase interface area per volume ⇒ porous electrode

homogenization ⇒ two coupled continua: conducting matrix electrolyte in pore space Solid-liquid interface is distributed in the volume Surface reaction terms → volume reaction terms Fluid flow → Darcy law Nernst-Planck equation → reactive transport in porous

medium

Newman/Tiedeman: Porous electrodes in batteries, double layer capacitors Neu/Krassowska, Pennachio/Savare/Colli/Franzone: “Bidomain equation”

(extracellular+intercellular space in cardiac tissue)

Moyne/Murad/Bennethum/Cushman: clays, concrete (coupled with swelling) Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 18 (46)

Potential distribution in porous electrodes −∇·ksσ∇φs =−aF jF −∂tQdl −∇·klκ∇φl =aF jF +∂tQdl jF =i0 F eα F

RT (φs−φl)

− i0 F e(α−1) F

RT (φs−φl)

Qdl =Cdla(φs −φl) φs solid matrix potential φl electrolyte potential Cdl double layer capacity Qdl double layer charge ks,kl pore space dependent factors σ conductivity of pure solid κ conductivity of pure electrolyte a interface area per volume i0 exchange current density Model limitations/extensions:

Pore size should exceed thickness of boundary layers Low concentration gradients Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 19 (46)

Proton Exchange Membrane Nafion conductivity (Tohoku Univ.) MEA by Gore Associates e.g. Nafion c DuPont:

“Solid acid”: sulfonic acid groups fixed at

PTFE backbone

Dissociation in presence of liquid water ⇒

free H3O+ ions available for conduction

⇒ conductivity depends on water content ⇒ water management needs to guarantee its

presence Problems associated with membrane

Reactant crossover Swelling free radical attacks from reaction

intermediates

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 20 (46)

Two phase flow

Cathode supply of O2 (gaseous, low solubility) remove H2O in gaseous and liquid form Anode (DMFC) supply of reactant dissolved in water removal of CO2 in dissolved and gaseous form Liquid H2O is needed in order to maintain membrane conductivity Phases need to move in opposite directions Pores kept open for gas flow by admixture of teflon Model by standard ansatz with modified capillary pressure/saturation curve Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 21 (46)

Mixed wettability: Bernoulli function based ansatz for s(pc)

0.2 0.4 0.6 0.8 1

• 40
• 30
• 20
• 10

10 20 30 40 Effective Saturation se Capillary Pressure pc 1:2 1:1 Bernoulli

Measured values for drying and wetting branches of se(pc) for quartz sand/teflon mixture with different compositions (triangles); Least squares fit to Bernoulli function based ansatz (lines) Bernoulli function: B(x) =

x ex−1

pδ

c =α(pc − p 1

2 )

se(pc) =−B′ B(pδ

c )βw +B(−pδ c )βn

• sw(pc) =sres

w +(1−sres g −sres w )se(pc) Ustohal et. al. (1998): measured saturation

curves for quartz sand/teflon mixtures of varying composition

Divisek et. al. (2003): Bernoulli function based

parametrization fitting well to these measurements

Measurements for fuel cell PTL ? Hysteresis ? Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 22 (46)

Ustohal/Stauffer/Dracos, J. Cont. Hydrol. 1998; Divisek/Fuhrmann/Gärtner/Jung, JES 2003

Stefan-Maxwell: two phase flow with gas mixture Stefan–Maxwell mean transport pore model for a mixture of ng gases empirically combined with two phase flow ∂t

• φ(1−sw) pgi

RT

• +∇·Ngi = Rgchem

i

+Rgev

i

Ngi = −krg

ng

∑

j=1

DSM

ij (p1,..., png)∇pg j

pgi partial pressure of ith component Ngi molar flux of ith component R gas constant T temperature krg(sw)

• rel. permeability

DSM

i j

Stefan Maxwell Diffusion matrix Rgev

i

evaporation reaction Rgchem

i

chemical reactions

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 23 (46)

Stefan–Maxwell Gas Transport – Diffusion Matrix Diffusion matrix DSM

ij

implicitely defined by RT

• Ngi

DK

i

+

ng

∑

i=1

yiNg j −yjNgi Dbm

ij

• = ptot∇yi+yi
• Bi

DK

i

+

ng

∑

i=1

Bi Dm

ij

yj

• 1− B j

Bi

• ∇ptot.

ptot = ∑ pgi total gas pressure yi = pgi/ptot molar fractions DK

i

Knudsen diffusions coeff. Dm

i j

binary diffusion coeff. λi mean free path gas i, ¯ r mean pore radius Ki = λi/¯ r Knudsen number Bi effective permeability

Definition

• f

capillary pressure: pc = ptot − pw

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 24 (46)

Further effects Evaporation + Condensation

Volumetric reaction between dissolved and gas phase Henry’s law gives equilibrium Kinetic constant ? Among other effects, the model contains reaction CO2|dissolved ⇌ CO2|gas

Heat transport Drag effect: force term in Darcy’s law form ion movement

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 25 (46)

Coupling between supply channels and porous electrodes Challenge: two phase coupling Droplets in H2 PEMFC CO2 bubbles in DMFC

Kandlikar et al., Mech Engr. 2009

• N. Paust et al. IMTEK, Freiburg Univ.

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 26 (46)

Mathematical/Numerical Model of a DMFC MEA 1D/2D/3D- Numerical model of a DMFC MEA including all the processes described above in cooperation with FZ Jülich

MEA processes taken into account in

temporal and spatial resolution.

11 coupled nonlinear partial differential

equations.

12 nonlinear algebraic equations per

discretization node.

Resolution of catalytic reaction chains. Two-phase flow model taking into account

mixed wettability.

0.5 1 500 1000 1500 Voltage (V) Current Density (A/m2)

• Meas. 0.5 mol/l
• Meas. 1 mol/l
• Meas. 2 mol/l

Sim., 0.5 mol/l Sim., 1 mol/l Sim., 2 mol/l

Measured and calculated polarization curves at 60oC.

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 27 (46)

Divisek/Fuhrmann/Gärtner/Jung, J. Electrochem. Soc. (2003)

Influence of anodic saturation curve

0.2 0.4 0.6 0.8 1

• 2.5 -2 -1.5 -1 -0.5 0

0.5 1 1.5 2 2.5 Saturation (m3/m3) Capillary Pressure (bar) p1/2=0.5 bar p1/2=1.0 bar p1/2=1.5 bar 0.5 1 500 1000 1500 Voltage (V) Current Density (A/m2)

• Meas. 0.5 mol/l

p1/2=0.5 bar p1/2=1.0 bar p1/2=1.5 bar

Capillary pressure vs. saturation polarization curves Best anode performance for more wettable material

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 28 (46)

• J. Fuhrmann, K Gärtner, in: Device and Materials Modeling in PEM Fuel Cells, vol. 113 of Topics in Applied Physics, 2009

Influence of cathodic saturation curve

0.2 0.4 0.6 0.8 1

• 2.5 -2 -1.5 -1 -0.5 0

0.5 1 1.5 2 2.5 Saturation (m3/m3) Capillary Pressure (bar) p1/2=0.5 bar p1/2=1.0 bar p1/2=1.5 bar 0.5 1 500 1000 1500 Voltage (V) Current Density (A/m2) measured, 0.5 mol/l most wettable medium wettable least wettable

Capillary pressure vs. saturation polarization curves Best cathode performance for less wettable material

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 29 (46)

• J. Fuhrmann, K Gärtner, in: Device and Materials Modeling in PEM Fuel Cells, vol. 113 of Topics in Applied Physics, 2009

Flow cell modeling Model simple xperimental devices to improve understanding of partial processes.

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 30 (46)

Flow cell modeling Model simple xperimental devices to improve understanding of partial processes. A+s B+s C +s Dissolved in fluid k±

A ⇃↾

k±

B ⇃↾

k±

C ⇃↾

Adsorption/Desorption Aads

k±

AB

⇀ ↽ Bads

k±

BC

⇀ ↽ Cads Adsorbed at catalyst surface ց ց H+ +e− H+ +e− Released ions Thin-layer flow cell

cAin qin Γin cAout cBout cCout Γout Mass Spec I

+ –

catalyst surface Ω Anode Cathode

How does the fraction of the desorbed intermediate B depend on the catalyst concentration ?

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 30 (46)

Model multistep surface reaction in a flow cell Rate expressions: rA = k+

A cAθf−k− A θA

rAB = k+

ABθA−k− ABθB

rB = k+

B cBθf −k− B θB

rBC = k+

BCθB−k− BCθC

rC = k+

C cCθf −k− C θC

Algebraic conditions for adsorbed species: dθA/dt−rA −rAB = 0 dθB/dt−rB +rAB −rBC = 0 dθC/dt−rC +rBC = 0

θX – fraction of catalyst sites occupied by adsorbed species θf = 1−θA −θB −θC – fraction of free catalyst sites Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 31 (46)

Solute transport coupled with surface reaction

Stationary diffusion and transport of species X = A,B,C in velocity field

v: ∂cX ∂t −div(DX gradcX − vcX) = 0

Boundary conditions:

cX = cX,in

• n Γin

Dirichlet (DX gradcX −cX v)· n = −cX v· n

• n Γout

Outflow (DX gradcX −cX v)· n = ce f f

cat rX

• n Γcat

Electrode reaction ∂cX ∂ n = 0

• n Γ\(Γin ∪Γout ∪Γcat)

No flow

Outflow of species X:

Xout =

• Γout

(DX gradcX −cX v)· n ds

Relative amount of reaction intermediates:

Iout = Bout/(Bout +Cout)

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 32 (46)

H2O2 yield in catalytic oxygen reduction

0.2 0.4 0.6 0.8 1 catalyst coverage 5 10 15 20 25 30 relative amount of reaction intermediate/%

O2(sol) H2O2(sol) H2O O2(ads) H2O2(ads) H2O(ads) Catalyst loading effect: H2O2 yield decreases with increasing coverage of active Pt nanodisks on planar glassy carbon substrate (nanostructured Pt/GC model electrodes) Top: Experimental, bottom: simulation.

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 33 (46)

Fuhrmann/Zhao/Langmach/Seidel/Jusys/Behm, Fuel Cells 2011

H2O2 yield in catalytic oxygen reduction

5 10 15 20 25 30 flow velocity /µL s

• 1

2 4 6 8 relative amount of reaction intermediate/%

O2(sol) H2O2(sol) H2O O2(ads) H2O2(ads) H2O(ads) H2O2 yield increases with increasing electrolyte flow velocity for

two different densities of similarly

sized Pt nanodisks

fully Pt covered electrode.

Top: Experimental, bottom: simulation.

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 33 (46)

Fuhrmann/Zhao/Langmach/Seidel/Jusys/Behm, Fuel Cells 2011

Choice of numerical methods Challenges:

Nonlinear coupling Dominant convection Heterogeneous species distributions

Wishes:

Maximum principle (no unphysical oscillations) Preservation of a priori bounds (positivity) Convergence, Accuracy, Efficiency

Choices:

Voronoi finite volume method: stability but limited convergence order

e.g. Fuhrmann/Linke/Langmach, Appl. Num. Math. 2011

Stabilized FEM: high convergence order, but no preservation of a priori bounds. John/Knobloch, Comput. Methods Appl. Mech. Engrg. 2008 Discontinuous Galerkin: high flexibility wrt. grids, high convergence order, but no

preservation of a priori bounds, complicated discrete systems

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 34 (46)

Diffusion-Convection-Reaction system N partial differential equations of convection/diffusion/reaction type in a one-, two or three-dimensional domain with appropriate boundary conditions. ∂tb( x,u)+∇· j( x,u)+r( x,u) = 0 u( x,t) = (u1( x,t)...uN( x,t)): vector valued function b( x,u) = (b1( x,u)...bN( x,u)): vector of storage terms r( x,u) = (r1( x,u)...rN( x,u)): vector of reaction terms

• j(

x,u) = ( j1( x,u)... jN( x,u)): vector of flux terms

• j(

x,u): sum of convective and diffusive fluxes.

Porous electrodes Two phase flow Density driven flow Charge transport in semiconductor devices Reaction-Diffusion systems ... Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 35 (46)

Two Point Flux Voronoi Finite Volume Method

• xK
• xL
• xM

K L

Integral over space-time REV control volume

0 = 1 tn −tn−1

tn

• tn−1
• K
• ∂tb(

x,u)+∇· j( x,u)+r( x,u)

• dxdt

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 36 (46)

Macneal, Quart. Math. Appl. 1953

Two Point Flux Voronoi Finite Volume Method

• xK
• xL
• xM

K L

Integral over space-time REV control volume Newton/Leibniz, Gauss, quadrature

0 = 1 tn −tn−1

tn

• tn−1
• K
• ∂tb(

x,u)+∇· j( x,u)+r( x,u)

• dxdt

0 = |K|

• b(

xK,un

K)−b(

xK,un−1

K

) tn −tn−1 +r( xK,un

K)

• −
• ∂K
• j(

x,un)· nds

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 36 (46)

Macneal, Quart. Math. Appl. 1953

Two Point Flux Voronoi Finite Volume Method

• xK
• xL
• xM

K L

Integral over space-time REV control volume Newton/Leibniz, Gauss, quadrature Approximation of inter-volume fluxes

0 = 1 tn −tn−1

tn

• tn−1
• K
• ∂tb(

x,u)+∇· j( x,u)+r( x,u)

• dxdt

0 = |K|

• b(

xK,un

K)−b(

xK,un−1

K

) tn −tn−1 +r( xK,un

K)

• −
• ∂K
• j(

x,un)· nds 0 = |K|

• b(

xK,un

K)−b(

xK,un−1

K

) tn −tn−1 +r( xK,un

K)

• −

∑

L nb. of K

|∂K ∩∂L| | xK − xL| g( xK, xL,un

K,un L) Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 36 (46)

Macneal, Quart. Math. Appl. 1953

Exponential fitting flux for linear problems Scalar, linear convection diffusion flux in given velocity field v( x): j = D∇u+uv( x) Flux projection onto control volume normal xK xL: vKL := 1 |∂K ∩∂L|

• ∂K∩∂L v(s)·(

xK − xL)ds Upwind numerical flux: g(uK,uL) := D

• B

vKL D

• uK −B
• −vKL

D

• uL
• (B(ξ) =

ξ e−ξ −1: Bernoulli function) Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 37 (46)

Allen/Southwell, Quart. J. Mech. Appl. Math. 1955; Il’in, Mat. Zametki 1969; Scharfetter/Gummel, IEEE Trans. El. Dev. 1969

Generalization to nonlinear problems Scalar, nonlinear convection diffusion flux in given velocity field v( x): j = D(u)∇u+F(u)v( x) Define g(uK,uL) := G from the solution of a local Dirichlet problem for the projection

• f the equation onto the edge

xK xL: Let w(ξ) := u( xL +ξ( xK − xL)):      D(w)w′ +F(w)vKL = G w(0) = uL w(1) = uK

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 38 (46)

Eymard/Fuhrmann/Gärtner, Numer. Math. 2006

Maximum principle Stationary convection diffusion problem: j = D∇u+uv( x), ∇· j = 0 Fundamental property of solution: “no unphysical oscillations”

Divergence free velocity field (∇·v = 0) ⇒ continuous maximum principle:

solution in a given point x is bounded by its values in a surroundig of x.

Discretely divergence free velocity projection vKL

∑

L nb. of K

|∂L∩∂K| | xK − xL| vKL = 0 ⇒ discrete maximum principle for discrete problem: solution in point xi is bounded by values in neigboring points xj.

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 39 (46)

Divergence free discrete fluxes How to guarantee discrete divergence free condition when coupling to flow problems?

From pointwise divergence free velocity field v(x) by exact calculation of

vKL = 1 |K ∩L|

• K∩L v(s)·(

xK − xL)ds

Analytical expressions (Hagen - Poiseuille, von Karman ...) Pointwisde divergence free finite elements Finite volume scheme for flow problem leading to discrete divergence free fluxes Darcy (porous media) flow: v = K∇p, vKL = K(pK − pL) Compatible finite volume solution for Navier-Stokes: work in progress Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 40 (46)

Divergence free finite elements

div(velocity space) ⊂ pressure space Lowest order Scott Vogelius elements:

(P

d,P −(d−1)) (d:space dimension) Stable on macro triangulations Maintain two independent discretizations for transport (FV) and for flow (FE) (FV): For every triangle T, calculate contributions σKL;T = ∂K ∩∂L∩T to

∂K ∩∂L

(FE): Calculate velocity projections vKL;T from continuous FE velocity field (FV): Assemble vKL from vKL;T Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 41 (46)

Burman/Linke, App. Num. Math 2008; Fuhrmann/Linke/Langmach, Appl. Num. Math. 2011

Compatible finite volume method After Nicolaides - extension of MAC scheme: Choice of unknowns:

Pressures pK on simplex mesh nodes Velocities qKL along simplex mesh edges ≡ normal to Voronoi box faces

Discrete operators divh,gradh,roth, vector calculus (up to boundary terms): divh = gradT

h

divh gradh = rotT

h roth −gradh divh

rotT

h gradh = 0

divh roth = 0 Discrete Stokes system:

• rotT

h roth v−gradh p

= 0 divh v = 0

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 42 (46)

Eymard/Fuhrmann/Linke, Proc. FVCA6, 2011

pdelib2/fvsys API Problem description

Grid Accumulation terms b(

xK,uK)

Reaction terms r(

xK,uK)

Fluxes between control volumes g(

xK, xL,uK,uL)

Species sets varying between subdomains and boundary parts

Solution strategy

Adaptive time step selection Nonlinear solver: Newton’s method Linear solver: PARDISO (direct), AMG with point block smoothers

Implementation

OpenMP based parallelization Portable to Linux/Unix/Mac/Windows OpenGL online graphics, FLTK GUI, Lua scripting language Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 43 (46)

www.pdelib.org

Advantages of Voronoï Methods

Focus on conservation character of mathematical model Straigthforward use of two point finite difference formulae, upwinding Convergence results for many cases Simplex based assembly loop as in FEM: data assembled from simplicial

contributions ⇒ known techniques for parallelization etc. apply

Heterogeneities represented by mesh

⇒ no averaging across interior boundaries

M-property + discrete maximum principle for Laplacian and properly upwinded

convection problems in 2/3D ⇒ conservation of qualitative properties (positivity, dissipativity) ⇒ physically meaningful a priori bounds

Implementation tolerant to Delaunay violations:

Algebraic expression for “|∂K ∩∂L|” may become negative, but still can be calculated

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 44 (46)

Challenges for Voronoï Methods Simple scheme but considerable efforts in mesh generation. No higher order version ⇒ high accuracy solutions need very fine meshes. Challenges

Reliable boundary conforming Delaunay meshing for general domains Anisotropic problems −

→ mesh alignment to anisotropy direction

Boundary + interior layer resolution −

→ anisotropic meshes aligned with layers

Storage terms may be overestimated at heterogeneity boundaries

− → “double layer” meshing at boundary

Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 45 (46)

Outlook Current activities:

CO Electrooxidation on Pt surfaces (Experiments of K. Krischer, Munich) Navier-Stokes/Porous media coupling with special emphasis on porous

electrodes (Poster of A. Caiazzo)

Corrosion modeling (with C.Chainais (Lille)) Nernst-Planck-Poisson-Navier-Stokes coupling using finite volume methods Electrochemistry and porous media · RICAM Workshop · Linz · 2011-10-05 · Page 46 (46)