Appendix: Details of the Numerical Implementation
Elimination of Surface Concentrations
To eliminate the surface concentrations from the Butler-Volmer equation, the mass transfer relation is rewritten to express surface concentrations in terms of bulk concentrations:
C R s = − i loc n F k m,R + C R C_{\text{R}}^{\text{s}} = -\frac{i_{\text{loc}}}{nFk_{\text{m,R}}} + C_{\text{R}} C R s = − n F k m,R i loc + C R
C O s = i loc n F k m,O + C O C_{\text{O}}^{\text{s}} = \frac{i_{\text{loc}}}{nFk_{\text{m,O}}} + C_{\text{O}} C O s = n F k m,O i loc + C O
Substituting into the concentration-dependent Butler-Volmer equation and dividing by i 0 C ref \frac{i_0}{C_{\text{ref}}} C ref i 0 :
C ref i 0 i loc = [ ( − i loc n F k m,R + C R ) exp ( α A F R T η ) − ( i loc n F k m,O + C O ) exp ( − α C F R T η ) ] \frac{C_{\text{ref}}}{i_0} i_{\text{loc}} = \left[\left(-\frac{i_{\text{loc}}}{nFk_{\text{m,R}}} + C_{\text{R}}\right) \exp\left(\frac{\alpha_{\text{A}} F}{RT}\eta\right) - \left(\frac{i_{\text{loc}}}{nFk_{\text{m,O}}} + C_{\text{O}}\right) \exp\left(\frac{-\alpha_{\text{C}} F}{RT}\eta\right)\right] i 0 C ref i loc = [ ( − n F k m,R i loc + C R ) exp ( RT α A F η ) − ( n F k m,O i loc + C O ) exp ( RT − α C F η ) ]
Collecting all terms containing i loc i_{\text{loc}} i loc on the left-hand side:
( C ref i 0 + 1 n F k m,R exp ( α A F R T η ) + 1 n F k m,O exp ( − α C F R T η ) ) i loc = C R exp ( α A F R T η ) − C O exp ( − α C F R T η ) \left(\frac{C_{\text{ref}}}{i_0} + \frac{1}{nFk_{\text{m,R}}} \exp\left(\frac{\alpha_{\text{A}} F}{RT}\eta\right) + \frac{1}{nFk_{\text{m,O}}} \exp\left(\frac{-\alpha_{\text{C}} F}{RT}\eta\right)\right) i_{\text{loc}} = C_{\text{R}} \exp\left(\frac{\alpha_{\text{A}} F}{RT}\eta\right) - C_{\text{O}} \exp\left(\frac{-\alpha_{\text{C}} F}{RT}\eta\right) ( i 0 C ref + n F k m,R 1 exp ( RT α A F η ) + n F k m,O 1 exp ( RT − α C F η ) ) i loc = C R exp ( RT α A F η ) − C O exp ( RT − α C F η )
Isolating i loc i_{\text{loc}} i loc yields the final expression:
i loc = i 0 [ C R C ref exp ( α A F R T η ) − C O C ref exp ( − α C F R T η ) ] 1 + i 0 n F k m,R C ref exp ( α A F R T η ) + i 0 n F k m,O C ref exp ( − α C F R T η ) i_{\text{loc}} = \frac{i_0 \left[\frac{C_{\text{R}}}{C_{\text{ref}}} \exp\left(\frac{\alpha_{\text{A}} F}{RT}\eta\right) - \frac{C_{\text{O}}}{C_{\text{ref}}} \exp\left(\frac{-\alpha_{\text{C}} F}{RT}\eta\right)\right]}{1 + \frac{i_0}{nFk_{\text{m,R}}C_{\text{ref}}} \exp\left(\frac{\alpha_{\text{A}} F}{RT}\eta\right) + \frac{i_0}{nFk_{\text{m,O}}C_{\text{ref}}} \exp\left(\frac{-\alpha_{\text{C}} F}{RT}\eta\right)} i loc = 1 + n F k m,R C ref i 0 exp ( RT α A F η ) + n F k m,O C ref i 0 exp ( RT − α C F η ) i 0 [ C ref C R exp ( RT α A F η ) − C ref C O exp ( RT − α C F η ) ]
Source Term Linearization
The approach is shown for the Φ s \Phi_{\text{s}} Φ s governing equation (the same procedure applies to Φ l \Phi_{\text{l}} Φ l ).
The source term is linearized around the known potential field of iteration n n n when solving for iteration n + 1 n+1 n + 1 :
∇ ⃗ ⋅ ( − σ eff ∇ ⃗ Φ s n + 1 ) + a ( i loc n + d d Φ s n i loc n ⋅ ( Φ s n + 1 − Φ s n ) ) = 0 \vec{\nabla} \cdot \left(-\sigma^{\text{eff}}\vec{\nabla}\Phi_{\text{s}}^{n+1}\right) + a\left(i_{\text{loc}}^{n} + \frac{d}{d\Phi_{\text{s}}^n}i_{\text{loc}}^{n} \cdot (\Phi_{\text{s}}^{n+1} - \Phi_{\text{s}}^{n})\right) = 0 ∇ ⋅ ( − σ eff ∇ Φ s n + 1 ) + a ( i loc n + d Φ s n d i loc n ⋅ ( Φ s n + 1 − Φ s n ) ) = 0
Computing d d Φ s n i loc n \frac{d}{d\Phi_{\text{s}}^n}i_{\text{loc}}^{n} d Φ s n d i loc n
Using the quotient rule:
d d Φ s n ( f g ) = g ⋅ f ′ − f ⋅ g ′ g 2 \frac{d}{d\Phi_{\text{s}}^n}\left(\frac{f}{g}\right) = \frac{g \cdot f' - f \cdot g'}{g^2} d Φ s n d ( g f ) = g 2 g ⋅ f ′ − f ⋅ g ′
where:
f = i 0 [ C R C ref exp ( α A F R T η ) − C O C ref exp ( − α C F R T η ) ] f = i_0 \left[\frac{C_{\text{R}}}{C_{\text{ref}}} \exp\left(\frac{\alpha_{\text{A}} F}{RT}\eta\right) - \frac{C_{\text{O}}}{C_{\text{ref}}} \exp\left(\frac{-\alpha_{\text{C}} F}{RT}\eta\right)\right] f = i 0 [ C ref C R exp ( RT α A F η ) − C ref C O exp ( RT − α C F η ) ]
f ′ = i 0 F R T [ α A C R C ref exp ( α A F R T η ) + α C C O C ref exp ( − α C F R T η ) ] f' = i_0 \frac{F}{RT} \left[\alpha_{\text{A}}\frac{C_{\text{R}}}{C_{\text{ref}}} \exp\left(\frac{\alpha_{\text{A}} F}{RT}\eta\right) + \alpha_{\text{C}} \frac{C_{\text{O}}}{C_{\text{ref}}} \exp\left(\frac{-\alpha_{\text{C}} F}{RT}\eta\right)\right] f ′ = i 0 RT F [ α A C ref C R exp ( RT α A F η ) + α C C ref C O exp ( RT − α C F η ) ]
g = 1 + i 0 n F k m,R C ref exp ( α A F R T η ) + i 0 n F k m,O C ref exp ( − α C F R T η ) g = 1 + \frac{i_0}{nFk_{\text{m,R}}C_{\text{ref}}} \exp\left(\frac{\alpha_{\text{A}} F}{RT}\eta\right) + \frac{i_0}{nFk_{\text{m,O}}C_{\text{ref}}} \exp\left(\frac{-\alpha_{\text{C}} F}{RT}\eta\right) g = 1 + n F k m,R C ref i 0 exp ( RT α A F η ) + n F k m,O C ref i 0 exp ( RT − α C F η )
g ′ = i 0 F α A n F k m,R C ref R T exp ( α A F R T η ) − i 0 F α C n F k m,O C ref R T exp ( − α C F R T η ) g' = \frac{i_0 F\alpha_{\text{A}}}{nFk_{\text{m,R}}C_{\text{ref}}RT} \exp\left(\frac{\alpha_{\text{A}} F}{RT}\eta\right) - \frac{i_0 F\alpha_{\text{C}}}{nFk_{\text{m,O}}C_{\text{ref}}RT} \exp\left(\frac{-\alpha_{\text{C}} F}{RT}\eta\right) g ′ = n F k m,R C ref RT i 0 F α A exp ( RT α A F η ) − n F k m,O C ref RT i 0 F α C exp ( RT − α C F η )
Substituting and simplifying yields the final expression:
d d Φ s n i loc n = F R T ( C R α A exp ( α A F R T η ) + C O α C exp ( − α C F R T η ) ) ( 1 + i 0 C ref n F ( 1 k m,R exp ( α A F R T η ) + 1 k m,O exp ( − α C F R T η ) ) ) 2 \frac{d}{d\Phi_{\text{s}}^n}i_{\text{loc}}^{n} = \frac{\frac{F}{RT}\left(C_{\text{R}}\alpha_{\text{A}} \exp\left(\frac{\alpha_{\text{A}}F}{RT}\eta\right) + C_{\text{O}}\alpha_{\text{C}} \exp\left(\frac{-\alpha_{\text{C}}F}{RT}\eta\right)\right)}{\left(1 + \frac{i_0}{C_{\text{ref}}nF}\left(\frac{1}{k_{\text{m,R}}} \exp\left(\frac{\alpha_{\text{A}}F}{RT}\eta\right) + \frac{1}{k_{\text{m,O}}} \exp\left(\frac{-\alpha_{\text{C}}F}{RT}\eta\right)\right)\right)^2} d Φ s n d i loc n = ( 1 + C ref n F i 0 ( k m,R 1 exp ( RT α A F η ) + k m,O 1 exp ( RT − α C F η ) ) ) 2 RT F ( C R α A exp ( RT α A F η ) + C O α C exp ( RT − α C F η ) )
+ i 0 C ref n F F R T ( α A + α C ) ( C R k m,O + C O k m,R ) exp ( ( α A − α C ) F R T Δ Φ ) ( 1 + i 0 C ref n F ( 1 k m,R exp ( α A F R T η ) + 1 k m,O exp ( − α C F R T η ) ) ) 2 + \frac{\frac{i_0}{C_{\text{ref}}nF}\frac{F}{RT}(\alpha_{\text{A}} + \alpha_{\text{C}})\left(\frac{C_{\text{R}}}{k_{\text{m,O}}} + \frac{C_{\text{O}}}{k_{\text{m,R}}}\right)\exp\left((\alpha_{\text{A}} - \alpha_{\text{C}})\frac{F}{RT}\Delta\Phi\right)}{\left(1 + \frac{i_0}{C_{\text{ref}}nF}\left(\frac{1}{k_{\text{m,R}}} \exp\left(\frac{\alpha_{\text{A}}F}{RT}\eta\right) + \frac{1}{k_{\text{m,O}}} \exp\left(\frac{-\alpha_{\text{C}}F}{RT}\eta\right)\right)\right)^2} + ( 1 + C ref n F i 0 ( k m,R 1 exp ( RT α A F η ) + k m,O 1 exp ( RT − α C F η ) ) ) 2 C ref n F i 0 RT F ( α A + α C ) ( k m,O C R + k m,R C O ) exp ( ( α A − α C ) RT F ΔΦ )