COMSOL 6.4 - Molten Carbonate Transport

This tutorial shows how to model transport of the individual ions in a salt melt comprising two binary salts, where the transport equations are defined using concentrated solution theory. The example model defines the electrolyte ion transport of a molten carbonate fuel cell (or electrolyzer), with a 1D-model geometry consisting of a separator, one negative (hydrogen) and one positive (oxygen) porous electrode.

The model is solved using a time-dependent solver, simulating changes in the ion distribution during a 1 h potentiostatic hold.

The model parameters are taken from a paper by Bodén and others (Ref. 1).

Model Definition

Figure 1 shows the model geometry. The electrodes are porous, where the solid matrix consists of Ni, alloyed with Cr or Al, on the negative side, and lithiated NiO on the positive side, respectively. The separator consists of a ceramic LiAlO2 matrix.

Figure 1: Model geometry.

The electrolyte consists of a Li+/K+/CO32- salt melt, and the aim of the model is to investigate the ion distribution in the cell during load. The absence of a neutral solvent, and the high concentration of the CO32- ions participating in the electrode reactions, makes the problem unsuitable to solve using the Nernst–Planck equations (which are used in the Tertiary Current Distribution, Nernst–Planck interfaces), which assume a diluted electrolyte, implying the presence of a solvent present at a high concentration in relation to all other species. The model is therefore defined using the interface, solving for an electrolyte phase potential variable and two neutral salt (Li2CO3 and K2CO3) fraction variables. The electrolyte transport equations in this interface make use of concentrated electrolyte theory, including the assumption of local electroneutrality (Ref. 2).

Charge transfer occurs at the electrode-electrolyte interface within the porous matrix of the negative and positive electrodes. The local current density, iloc, of the negative electrode hydrogen oxidation charge transfer reaction

is defined as

(1)

where β is a symmetry factor and the exchange current density i0,neg is defined as

(2)

where xi denote the molar fractions of the reacting gas species in the gas phase. These molar fractions are set as constant in the model. (The i0,ref exchange current density refers to the exchange current density for an imaginary reference state of unity molar fractions.)

Correspondingly, on the positive side, the local current density of the oxygen reduction reaction

is defined as

(3)

with

(4)

In the kinetic expressions above, the overpotential η is defined as

(5)

Here, the electrode phase potential ϕs is set to 0 V in the negative electrode and 0.75 V in the positive electrode. The cell voltage equals the difference between the positive and negative electrode phase potentials

(6)

By defining the node stoichiometry to correspond to the oxidation of CO32-, the electrolyte phase potential ϕl is inherently defined to include CO32- ion concentration overpotentials.

The equilibrium potentials of the two electrode reactions are defined as

(7)

and

(8)

where the reference equilibrium potentials refer to the reference state of 1 atm for the partial pressures and the cell temperature 650°C.

Finally, the species sources, R (mol/m3s), for the CO32- ions are calculated using Faraday’s law of electrolysis according to

(9)

where Av (m2/m3) is the specific surface are of the electrodes.

A domain node is used in the model to define the charge transfer reactions according to the above equations.

As CO32- is transferred between the electrolyte and electrode phases, current starts to flow in the cell. The transport in the electrolyte is characterized by the use of composition-dependent binary Maxwell-Stefan diffusivities, these are defined on the node, which is a child node to the domain node. One diffusivity for each ion-pair (Li+-K+, Li+-CO32-, and K+-CO32-) is required to define the transport.

The Maxwell–Stefan transport model computes diffusive–migrative fluxes of the electrolyte species with respect to the mass-averaged velocity, u (m/s), of the electrolyte. Instead of computing the velocity using a Fluid Flow interface such as Darcy’s Law, the velocity is approximated to be proportional to the electrolyte current vector, Il (A/m2), times the carbonate ion molar mass according to

(10)

The electrolyte density, ρ (kg/m3), will in reality depend on composition, but is treated as a constant in this tutorial.

For this cell, where the carbonate ion is the only ion (CO32-) participating in the charge transfer reaction, the above expression for the electrolyte mass-averaged velocity is valid for a stationary current distribution. However, it often serves as a fair approximation also for time-dependent problems.

All the three domains are porous and are defined using a common node. By the use of a domain-specific porosity variable, the porosity is defined differently in each domain.

The model is solved in a comprising a step followed by a study step. The first step solves for the electrolyte phase potential dependent variable only, using a stationary solver, and the solution is then passed on as initial values the following time-dependent step.

Results and Discussion

Figure 2 shows the individual ion concentrations at 1 h. As a result of electroneutrality, the concentration of CO32- equals half of the sum of Li+ and K+.

Interestingly, the concentration gradients of the two cations (Li+ and K+) are opposite to each other, and in addition it can be seen that the smaller ion Li+, which has a higher binary diffusivity value vs CO32-, is subject to higher concentration gradient magnitudes than the larger K+. This segregation effect is related to it being energetically more favorable to have Li+ fulfill electroneutrality when gradients in the CO32- concentration are induced by the electrode reactions. Binary diffusivity coefficients can be seen an as inverse proportionality constants for the friction forces between the electrolyte species — the higher the diffusivity, the lower the friction.

An alternative way of plotting composition changes in the cell is to plot the relative proportions of the two neutral salts (Li2CO3 and K2CO3) in relation to each other. Figure 3 plots the quantity xLi+/(xLi++xK+) as a measure of the fraction of Li2CO3.

Figure 2: Individual ion concentrations.

Figure 3: Li2CO3 salt fraction.

Figure 4 shows the electrolyte conductivity profile over time. As a result of the concentration changes, the conductivity varies. This is a result of the composition-dependent diffusivities deployed in the model.

Figure 4: Electrolyte conductivity.

The electrolyte conductivity variations are fairly small. Figure 5 shows the electrolyte phase potential for all simulated times, ranging from 0 to 1 h in steps of 0.1 h. The variations in the potential profile are small and all line graphs more or less coincide.

Figure 5: Electrolyte phase potential.

Figure 6 shows the volumetric current densities in the cell stemming from the charge transfer reactions, also here ranging from 0 to 1 h in steps of 0.1 h. These source terms give rise to the corresponding electrolyte current density profiles seen in Figure 7. As for Figure 5, variations over time are small and all line graphs more or less coincide.

It is known that the local composition of the electrolyte impacts the electrolyte wetting of the porous matrix of the electrodes. An extended model could incorporate this effect by defining the local electrolyte volume fraction as a function of the composition. Such a model would likely result in larger spatial and temporal variations in the current distribution.