Stabilising Effects of Lumped Integration Schemes for the Simulation of Metal-Electrolyte Reactions

Throughout this paper, we consider the domain Ω of the type shown in Fig. 1. This domain consists of an electrolyte region Ω_e , characterised by its primary fields the ionic species concentrations C_π and the electrolyte potential φ, and a neighbouring metal region Ω_m described by the displacement field vector u and the lattice hydrogen concentration C_L . Separating these two sub-domains is the metal-electrolyte interface Γ_int which is included through the hydrogen surface coverage θ_ads = C_ads /N_ads , where N_ads is the number of surface adsorption sites and C_ads is the surface concentration. The metal domain contains a crack, as near crack tips is where the critical content of hydrogen triggering embrittlement is attained.

Zoom In Zoom Out Reset image size

Regarding notation and units. For consistency between the metal and electrolyte domains, the ionic concentrations in the electrolyte and the hydrogen concentration within the metal are expressed in the SI units mol m⁻³. Consequently, for all constants related to these quantities, such as the reaction constants, their units follow accordingly; e.g., the water auto-ionization constant is given by ${10}^{-8}{(\mathrm{mol}/{{\rm{m}}}^{3})}^{2}$ (as opposed to the more common terminology of 10⁻¹⁴(mol/l)²). As a reference point for the electrolyte potential, the standard hydrogen electrode is used, placing the equilibrium constant of the hydrogen evolution reaction at 0 V_SHE. Finally, we use lightface italic letters for scalars, e.g. C_L , upright bold letters for vectors, e.g. u, and bold italic letters, such as σ or B , for second and higher order tensors.

Electrolyte sub-domain

For the electrolyte sub-domain, we consider the following ionic species: H⁺ and OH⁻ to represent a water-like electrolyte, Na⁺ and Cl⁻ to approximate seawater, and Fe²⁺ and FeOH⁺ to track the products produced by corrosion. We assume that the solubility of any gas phases produced by reactions is negligible, and that these gases disappear from the domain immediately.

The transport of ionic species inside the electrolyte is characterised by the Nernst-Planck mass balance, which describes changes in concentration based on advection, diffusion due to gradients in concentration, and electromigration due to gradients in electric potential. Thus, for a volume reaction rate R_π , charge of the ionic species z_π , diffusivity D_π , Faraday constant F, reference temperature T, and gas constant R, the Nernst-Planck equation reads

$$\begin{eqnarray}&&\begin{array}{c}{\dot{C}}_{\pi }+({\bf{v}}\cdot {\boldsymbol{\nabla }}){C}_{\pi }+{\boldsymbol{\nabla }}\cdot \left(-{D}_{\pi }{\boldsymbol{\nabla }}{C}_{\pi }\right)\\ \quad +\,\frac{{z}_{\pi }F}{{RT}}{\boldsymbol{\nabla }}\cdot \left(-{D}_{\pi }{C}_{\pi }{\boldsymbol{\nabla }}\varphi \right)+{R}_{\pi }=0\end{array}\end{eqnarray} \tag{ 1 }$$

where v represents the fluid velocity of the electrolyte, which we will assume equal to zero in the remainder to solely focus on reaction-induced convergence issues and oscillations, instead of including advection-based transport and the stability issues associated with it. ^28–34

In addition to the mass balance, we require the local electric charge to be conserved. Several popular formulations for this charge-conservation include tracking the total current flow and requiring it to be divergence-free, ∇ · (∑z_π j_π ) = 0, ^35,36 or using Gauss law to describe the charge conservation within the electrolyte and the electron double-layer, ∇² φ = − F/ε ∑z_π C_π . ^37,38 However, since this double layer is extremely thin and infeasible to account for through the discretisation of the electrolyte domain, ³⁹ we have chosen to adopt the electroneutrality condition because of its simplicity; such that

$$\begin{eqnarray}&&\displaystyle \sum _{\pi }{z}_{\pi }{C}_{\pi }=0\end{eqnarray} \tag{ 2 }$$

with this equation inducing changes in electrolyte potential to enforce the redistribution of species through the electrolyte potential-based diffusion term from Eq. 1.

In addition to these conservation laws, we consider several reactions to occur within the electrolyte. Since a water-like electrolyte is considered, the water auto-ionisation reaction is relevant: ⁴⁰

$$\begin{eqnarray}&&{{\rm{H}}}_{2}{\rm{O}}\underset{{k}_{w}^{{\prime} }}{\overset{{k}_{w}}{\rightleftharpoons }}{{\rm{H}}}^{+}+{\mathrm{OH}}^{-}\end{eqnarray} \tag{ 3 }$$

and it is implemented through the reaction term R_π from Eq. 1 as:

$$\begin{eqnarray}\begin{array}{rcl}{R}_{{{\rm{H}}}^{+},w} & = & {R}_{{\mathrm{OH}}^{-}}={k}_{w}{C}_{{{\rm{H}}}_{2}{\rm{O}}}-{k}_{w}^{{\prime} }{C}_{{{\rm{H}}}^{+}}{C}_{{\mathrm{OH}}^{-}}\\ & = & {k}_{{eq}}\left({K}_{w}-{C}_{{{\rm{H}}}^{+}}{C}_{{\mathrm{OH}}^{-}}\right)\end{array}\end{eqnarray} \tag{ 4 }$$

Here, we assume equilibrium between the forward reaction rate, governed through k_w , and the backward reaction rate ${k}_{w}^{{\prime} }$. This allows us to rewrite the reaction rate in terms of the equilibrium constant K_w = 10⁻⁸ mol² m⁻⁶ and a penalty term-like constant k_eq , where this penalty constant is chosen to be high enough to enforce the equilibrium reaction to be fulfilled at all times. An alternative to this penalty approach is using the equilibrium reaction to eliminate ${C}_{{\mathrm{OH}}^{-}}$ from the governing equations. However, this removes the boundary influx terms and would therefore complicate applying (electro-) chemical reactions involving this species on the metal-electrolyte boundary. It is furthermore not an approaches that can be generalised for the case of multiple equilibrium reactions. Therefore, we have chosen to use the penalty approach to enforce the equilibrium and resolve the complications with regards to stability and oscillations, allowing for an easy and straightforward integration of bulk and surface reactions involving an arbitrary amount of reaction and involved reaction species.

For the corrosion products, we consider the following reaction:

$$\begin{eqnarray}&&{\mathrm{Fe}}^{2+}+{{\rm{H}}}_{2}{\rm{O}}\underset{{k}_{{fe}}^{{\prime} }}{\overset{{k}_{{fe}}}{\rightleftharpoons }}{\mathrm{FeOH}}^{+}+{{\rm{H}}}^{+}\end{eqnarray} \tag{ 5 }$$

which in turn can react further through:

$$\begin{eqnarray}&&{\mathrm{FeOH}}^{+}+{{\rm{H}}}_{2}{\rm{O}}\mathop{\rightharpoonup }\limits^{{k}_{{feoh}}}\mathrm{Fe}{(\mathrm{OH})}_{2}+{{\rm{H}}}^{+}\end{eqnarray} \tag{ 6 }$$

These equations are governed by their reaction rates: k_fe (forward) and ${k}_{{fe}}^{{\prime} }$ (backward) for Reaction 5, and k_feoh for Reaction 6. Using these rate constants, the reaction rates can be written as: ⁴¹

$$\begin{eqnarray}&&{R}_{{\mathrm{Fe}}^{2+}}=-{k}_{{fe}}{C}_{{\mathrm{Fe}}^{2+}}+{k}_{{fe}}^{{\prime} }{C}_{{\mathrm{FeOH}}^{+}}{C}_{{{\rm{H}}}^{+}}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{R}_{{\mathrm{FeOH}}^{+}}={k}_{{fe}}{C}_{{{Fe}}^{2+}}-{C}_{{\mathrm{FeOH}}^{+}}({k}_{{feoh}}+{k}_{{fe}}^{{\prime} }{C}_{{{\rm{H}}}^{+}})\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{R}_{{{\rm{H}}}^{+},{fe}}={k}_{{fe}}{C}_{{\mathrm{Fe}}^{2+}}-{C}_{{\mathrm{FeOH}}^{+}}({k}_{{fe}}^{{\prime} }{C}_{{{\rm{H}}}^{+}}-{k}_{{feoh}})\end{eqnarray} \tag{ 9 }$$

From Eqs. 4 and 9, a hydrogen reaction rate ${R}_{{{\rm{H}}}^{+}}$ can be defined that encompasses both water ionisation and corrosion contributions: ${R}_{{{\rm{H}}}^{+}}={R}_{{{\rm{H}}}^{+},w}+{R}_{{{\rm{H}}}^{+},{fe}}$.

Metal sub-domain

The deformation of the metal sub-domain is described by the balance of linear momentum; such that, neglecting body forces and inertia terms,

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}\cdot {\boldsymbol{\sigma }}=0\end{eqnarray} \tag{ 10 }$$

where the stress tensor σ is obtained assuming linear-elastic material behavior. As it is typically the case, no direct influence of hydrogen on the elastic behavior of the solid is defined.

For the description of the hydrogen present within the metal, we start with the chemical potential of the dissolved hydrogen:

$$\begin{eqnarray}&&\mu ={\mu }_{0}+{RT}\,\mathrm{ln}\left(\displaystyle \frac{{\theta }_{L}}{1-{\theta }_{L}}\right)-{\overline{V}}_{H}{\sigma }_{H}\end{eqnarray} \tag{ 11 }$$

using the lattice site occupancy θ_L = C_L /N_L , the hydrostatic stresses ${\sigma }_{H}=\mathrm{tr}({\boldsymbol{\sigma }})/3$, the lattice sites density N_L , and the partial molar volume of hydrogen ${\overline{V}}_{H}$. This chemical potential is used to calculate the hydrogen fluxes through:

$$\begin{eqnarray}&&{{\bf{J}}}_{L}=-\displaystyle \frac{{D}_{L}{C}_{L}}{{RT}}{\boldsymbol{\nabla }}\mu \end{eqnarray} \tag{ 12 }$$

where D_L is the lattice diffusivity. These definitions result in the following mass balance for the lattice hydrogen:

$$\begin{eqnarray}&&{\dot{C}}_{L}+{\boldsymbol{\nabla }}\cdot \left(-\displaystyle \frac{{D}_{L}}{1-{C}_{L}/{N}_{L}}{\boldsymbol{\nabla }}{C}_{L}\right)+{\boldsymbol{\nabla }}\cdot \left(\displaystyle \frac{{D}_{L}{C}_{L}{\overline{V}}_{H}}{{RT}}{\boldsymbol{\nabla }}{\sigma }_{H}\right)=0\end{eqnarray} \tag{ 13 }$$

In contrast to conventional formulations, ^42–44 we do not assume a low lattice occupancy and as a result we obtain a hydrogen transport equation that depends not only on the concentration gradient but also on the concentration itself through the addition of a 1/(1 − C_L /N_L ) factor. This allows the present formulation to be valid over a larger range of parameters, increasing the local diffusivity to prevent the lattice hydrogen atoms to exceed the number of lattice sites (for instance, as caused by the presence of high hydrostatic stress gradients). A consequence of this is that the analytical expression for the hydrogen concentration at steady state is no longer estimated using ${C}_{L}={C}_{0}\,\exp \left({\overline{V}}_{H}{\sigma }_{H}/({RT})\right)$ ⁴⁵ (where C₀ is the concentration away from stress concentrators), but from setting Eq. 11 to be constant, resulting in:

$$\begin{eqnarray}\begin{array}{rcl}{C}_{L} & = & \displaystyle \frac{1-{C}_{L}/{N}_{L}}{1-{C}_{0}/{N}_{L}}{C}_{0}\,\exp \left(\displaystyle \frac{{\overline{V}}_{H}{\sigma }_{H}}{{RT}}\right)\\ & = & \displaystyle \frac{1}{1+\tfrac{{C}_{0}/{N}_{L}}{1-{C}_{0}/{N}_{L}}\exp \left(\tfrac{{\overline{V}}_{H}{\sigma }_{H}}{{RT}}\right)}\cdot {C}_{0}\,\exp \left(\displaystyle \frac{{\overline{V}}_{H}{\sigma }_{H}}{{RT}}\right)\end{array}\end{eqnarray} \tag{ 14 }$$

While Eq. 14 is not directly used within the model and thus solely provided for context, it shows that the maximum interstitial lattice concentration will tend closer to the reference concentration C₀ when these concentrations become closer to the number of interstitial lattice sites. As such, the influence of stress concentrators is less pronounced for higher C₀ values. It should also be noted that hydrogen sequestration in microstructural traps is not accounted for in the present model. While hydrogen traps can easily be included within Eq. 13, ^17,46,47 we have chosen not to do so to retain the focus on the stability and convergence behavior of the electrolyte and the metal-electrolyte interface.

Metal-electrolyte interface

At the metal-electrolyte interface, we consider the following reactions: ^19,48–51

$$\begin{eqnarray}&&\mathrm{Volmer}\,(\mathrm{acid}):\qquad {{\rm{H}}}^{+}+{\rm{M}}+{{\rm{e}}}^{-}\underset{{k}_{{Va}}^{{\prime} }}{\overset{{k}_{{Va}}}{\rightleftharpoons }}{\mathrm{MH}}_{{ads}}\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&\mathrm{Heyrovsky}\,(\mathrm{acid}):\qquad {{\rm{H}}}^{+}+{{\rm{e}}}^{-}+{\mathrm{MH}}_{{ads}}\underset{{k}_{{Ha}}^{{\prime} }}{\overset{{k}_{{Ha}}}{\rightleftharpoons }}{\rm{M}}+{{\rm{H}}}_{2}\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&\mathrm{Volmer}\,(\mathrm{base}):{{\rm{H}}}_{2}{\rm{O}}+{\rm{M}}+{{\rm{e}}}^{-}\underset{{k}_{{Vb}}^{{\prime} }}{\overset{{k}_{{Vb}}}{\rightleftharpoons }}{\mathrm{MH}}_{{ads}}+{\mathrm{OH}}^{-}\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&\mathrm{Heyrovsky}\,(\mathrm{base}):\qquad {{\rm{H}}}_{2}{\rm{O}}+{{\rm{e}}}^{-}+{\mathrm{MH}}_{{ads}}\underset{{k}_{{Hb}}^{{\prime} }}{\overset{{k}_{{Hb}}}{\rightleftharpoons }}{\rm{M}}+{{\rm{H}}}_{2}+{\mathrm{OH}}^{-}\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&\mathrm{Tafel}:\qquad 2{\mathrm{MH}}_{{ads}}\underset{{k}_{T}^{{\prime} }}{\overset{{k}_{T}}{\rightleftharpoons }}2{\rm{M}}+{{\rm{H}}}_{2}\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&\mathrm{Absorption}:\qquad {\mathrm{MH}}_{{ads}}\underset{{k}_{A}^{{\prime} }}{\overset{{k}_{A}}{\rightleftharpoons }}{\mathrm{MH}}_{{abs}}\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&\mathrm{Corrosion}:\qquad {\mathrm{Fe}}^{2+}+2{{\rm{e}}}^{-}\underset{{k}_{c}^{{\prime} }}{\overset{{k}_{c}}{\rightleftharpoons }}\mathrm{Fe}\end{eqnarray} \tag{ 21 }$$

where the Volmer reactions are the main source of adsorbed hydrogen, either through the acidic Volmer reaction for low pH environments or through the much slower basic Volmer reaction for alkaline environments. For a low surface coverage and highly negative metal potentials, both Heyrovsky reactions are the main sink of adsorbed hydrogen, whereas for higher surface coverage the Tafel reaction becomes a more prominent hydrogen sink. On the other side, the absorption reaction characterises the quantity of hydrogen that goes from being attached to the surface to entering the bulk metal. For the corrosion reaction, we assume its rate to be insufficient to significantly dissolve the metal. While it is possible to include metal dissolution due to corrosion through smeared approaches such as phase field ^52–54 or interface tracking schemes, ^55,56 we choose here to focus on the stability of the surface reactions when they are explicitly represented within the geometry. For the same reason, no passivation or protective layer development/dissolution is included within the context of the corrosion reactions. However, these phenomena can be readily included by adding or altering the corresponding electro-chemical reactions. ^57,58 Oxygen reduction reactions have also not been included. As is common within numerical simulations and experimental setups, ^36,59,60 a constant electric potential is imposed on the metal. By applying this potential, the anodic and cathodic reaction rates decouple, and electric currents are allowed to enter and leave through the metal (the corrosion reaction can produce electrons, but these do not necessarily need to be consumed by the hydrogen and oxygen reactions). This reduces the role of the oxygen reduction reaction to providing a source of OH⁻ ions, which will locally increase the pH slightly but otherwise not influence the results. As such, the effect of not including the oxygen reduction reaction is expected to be limited.

The reaction rates for the reactions 15–21 are given by: ²⁰

	Forward	Backward
Volmer (acid):	${\nu }_{{Va}}={k}_{{Va}}{C}_{{{\rm{H}}}^{+}}(1-{\theta }_{{ads}}){e}^{-{\alpha }_{{Va}}\tfrac{\eta F}{{RT}}}$	${\nu }_{{Va}}^{{\prime} }={k}_{{Va}}^{{\prime} }{\theta }_{{ads}}{e}^{(1-{\alpha }_{{Va}})\tfrac{\eta F}{{RT}}}$	[22]
Heyrovsky (acid):	${\nu }_{{Ha}}={k}_{{Ha}}{C}_{{{\rm{H}}}^{+}}{\theta }_{{ads}}{e}^{-{\alpha }_{{Ha}}\tfrac{\eta F}{{RT}}}$	${\nu }_{{Ha}}^{{\prime} }={k}_{{Ha}}^{{\prime} }(1-{\theta }_{{ads}}){p}_{{{\rm{H}}}_{2}}{e}^{(1-{\alpha }_{{Ha}})\tfrac{\eta F}{{RT}}}$	[23]
Volmer (base):	${\nu }_{{Vb}}={k}_{{Vb}}(1-{\theta }_{{ads}}){e}^{-{\alpha }_{{Vb}}\tfrac{\eta F}{{RT}}}$	${\nu }_{{Vb}}^{{\prime} }={k}_{{Vb}}^{{\prime} }{C}_{{\mathrm{OH}}^{-}}{\theta }_{{ads}}{e}^{(1-{\alpha }_{{Vb}})\tfrac{\eta F}{{RT}}}$	[24]
Heyrovsky (base):	${\nu }_{{Hb}}={k}_{{Hb}}{\theta }_{{ads}}{e}^{-{\alpha }_{{Hb}}\tfrac{\eta F}{{RT}}}$	${\nu }_{{Hb}}^{{\prime} }={k}_{{Hb}}^{{\prime} }(1-{\theta }_{{ads}}){p}_{{{\rm{H}}}_{2}}{C}_{{\mathrm{OH}}^{-}}{e}^{(1-{\alpha }_{{Hb}})\tfrac{\eta F}{{RT}}}$	[25]
Tafel:	${\nu }_{T}={k}_{T}\left|{\theta }_{{ads}}\right|{\theta }_{{ads}}$	${\nu }_{T}^{{\prime} }={k}_{T}^{{\prime} }(1-{\theta }_{{ads}}){p}_{{{\rm{H}}}_{2}}$	[26]
Absorption:	νA = kA (NL − CL )θads	${\nu }_{A}^{{\prime} }={k}_{A}^{{\prime} }{C}_{L}(1-{\theta }_{{ads}})$	[27]
Corrosion:	${\nu }_{c}={k}_{c}{C}_{{\mathrm{Fe}}^{2+}}{e}^{-{\alpha }_{c}\tfrac{\eta F}{{RT}}}$	${\nu }_{c}^{{\prime} }={k}_{c}^{\prime} {e}^{(1-{\alpha }_{c})\tfrac{\eta F}{{RT}}}$	[28]

These reactions use constant forward and backward rate constants k and $k^{\prime}$, charge transfer coefficient α, and equilibrium potential E_eq . The reaction rates depend on the ionic concentrations within the electrolyte, C_π , the hydrogen surface coverage, θ_ads , the lattice hydrogen concentration C_L , and the metal and electrolyte potentials through the overpotential, η = E_m − φ − E_eq . As a result, they depend on and provide the coupling between the electrolyte and metal domain through the reaction fluxes as:

$$\begin{eqnarray}&&{\overline{J}}_{{{\rm{H}}}^{+}}=-({\nu }_{{Va}}-{\nu }_{{Va}}^{{\prime} })-{\nu }_{{Ha}}\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&{\overline{J}}_{{\mathrm{OH}}^{-}}={\nu }_{{Vb}}-{\nu }_{{Vb}}^{{\prime} }+{\nu }_{{Hb}}\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&{\overline{J}}_{{\mathrm{Fe}}^{2+}}={\nu }_{c}^{{\prime} }-{\nu }_{c}\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}&&{\overline{J}}_{L}={\nu }_{A}-{\nu }_{A}^{{\prime} }\end{eqnarray} \tag{ 32 }$$

where we have neglected the backwards reaction rates of the Heyrovsky and Tafel reactions based on the assumption that the hydrogen gas produced disappears from the domain, resulting in a negligible reaction rate (${\nu }_{T}^{\prime} ={\nu }_{{Ha}}^{{\prime} }={\nu }_{{Hb}}^{\prime} =0$). Finally, the surface mass balance for the adsorbed hydrogen reads:

$$\begin{eqnarray}&&\begin{array}{l}{N}_{{ads}}{\dot{\theta }}_{{ads}}-({\nu }_{{Va}}-{\nu }_{{Va}}^{{\prime} })+{\nu }_{{Ha}}+2{\nu }_{T}+({\nu }_{A}-{\nu }_{A}^{{\prime} })\\ \,\,\,-\,({\nu }_{{Vb}}-{\nu }_{{Vb}}^{{\prime} })+{\nu }_{{Hb}}=0\end{array}\end{eqnarray} \tag{ 33 }$$

Discretisation

We discretise the governing equations presented in the previous sub-sections using the finite element method, interpolating the electrolyte potential and concentrations using quadratic triangular elements as:

$$\begin{eqnarray}&&\varphi =\displaystyle \sum _{{el}}{{\bf{N}}}_{\varphi }^{{el}}{{\boldsymbol{\varphi }}}^{{el}}\qquad {C}_{\pi }=\displaystyle \sum _{{el}}{{\bf{N}}}_{c}^{{el}}{{\bf{C}}}_{\pi }\end{eqnarray} \tag{ 34 }$$

Similarly, for the metal sub-domain, the displacements and lattice hydrogen concentration are interpolated as:

$$\begin{eqnarray}&&{\bf{u}}=\displaystyle \sum _{{el}}{{\boldsymbol{N}}}_{u}^{{el}}{{\bf{u}}}^{{el}}\qquad {C}_{L}=\displaystyle \sum _{{el}}{{\bf{N}}}_{L}^{{el}}{{\bf{C}}}_{L}^{{el}}\end{eqnarray} \tag{ 35 }$$

and we represent the hydrogen surface coverage using one-dimensional quadratic line elements at the metal-electrolyte interface, such that

$$\begin{eqnarray}&&{\theta }_{{ads}}=\displaystyle \sum _{{iel}}{{\bf{N}}}_{\theta }^{{el}}{{\boldsymbol{\theta }}}^{{el}}\end{eqnarray} \tag{ 36 }$$

In addition to this spatial discretisation, we perform the temporal discretisation using a backward Euler scheme, evaluating the governing equations at t + Δt and using the time derivative:

$$\begin{eqnarray}&&\dot{\square }=\displaystyle \frac{{\square }^{t+{\rm{\Delta }}t}-{\square }^{t}}{{\rm{\Delta }}t}\end{eqnarray} \tag{ 37 }$$

Using these discretisations, the governing equations for the metal domain, Eqs. 10 and 13, are transformed into their discretised weak form, with the following force components

$$\begin{eqnarray}&&{{\bf{f}}}_{u}={\int }_{{{\rm{\Omega }}}_{m}}{{\boldsymbol{B}}}_{u}^{T}{\boldsymbol{D}}{{\boldsymbol{B}}}_{u}{{\bf{u}}}^{t+{\rm{\Delta }}t}\,{\rm{d}}{{\rm{\Omega }}}_{m}-{\int }_{{{\rm{\Gamma }}}_{m}}{{\boldsymbol{N}}}_{u}^{T}{{\bf{t}}}_{{ext}}\,{\rm{d}}{{\rm{\Gamma }}}_{m}={\bf{0}}\end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray}\begin{array}{rcl}{{\bf{f}}}_{L} & = & {\displaystyle \int }_{{{\rm{\Omega }}}_{m}}\displaystyle \frac{1}{{\rm{\Delta }}t}{{\bf{N}}}_{L}^{T}{{\bf{N}}}_{L}\left({{\bf{C}}}_{L}^{t+{\rm{\Delta }}t}-{{\bf{C}}}_{L}^{t}\right){\rm{d}}{{\rm{\Omega }}}_{m}\\ & & +{\displaystyle \int }_{{{\rm{\Omega }}}_{m}}\displaystyle \frac{{D}_{L}}{1-{{\bf{N}}}_{L}{{\bf{C}}}_{L}^{t+{\rm{\Delta }}t}/{N}_{L}}{\left({\boldsymbol{\nabla }}{{\bf{N}}}_{L}\right)}^{T}{\boldsymbol{\nabla }}{{\bf{N}}}_{L}{{\bf{C}}}_{L}^{t+{\rm{\Delta }}t}\,{\rm{d}}{{\rm{\Omega }}}_{m}\\ & & -{\displaystyle \int }_{{{\rm{\Omega }}}_{m}}\displaystyle \frac{E}{3(1-2\nu )}\displaystyle \frac{{D}_{L}{\overline{V}}_{H}}{{RT}}{\left({\boldsymbol{\nabla }}{{\bf{N}}}_{L}\right)}^{T}\left({{\bf{N}}}_{L}{{\bf{C}}}_{L}^{t+{\rm{\Delta }}t}\right){{\boldsymbol{B}}}_{u}^{* }{{\bf{u}}}^{t+{\rm{\Delta }}t}\,{\rm{d}}{{\rm{\Omega }}}_{m}\\ & & -{\displaystyle \int }_{{{\rm{\Gamma }}}_{{int}}}{{\bf{N}}}_{L}^{T}\left({\nu }_{A}^{t+{\rm{\Delta }}t}-{\nu }_{A}^{{\prime} t+{\rm{\Delta }}t}\right)\,{\rm{d}}{{\rm{\Gamma }}}_{{int}}-{\displaystyle \int }_{{{\rm{\Gamma }}}_{m}}{{\bf{N}}}_{L}^{T}{j}_{L}\,{\rm{d}}{{\rm{\Gamma }}}_{m}=0\end{array}\end{eqnarray} \tag{ 39 }$$

where the boundary conditions on the external boundary Γ_m are defined as the external tractions t_ext and external hydrogen influx j_L . The displacement to strain and displacement to strain gradient mapping matrices are defined under plane-strain conditions as:

$$\begin{eqnarray}{{\boldsymbol{B}}}_{u}=\left[\begin{array}{cccccc}\tfrac{\partial {N}_{u1}}{\partial x} & \tfrac{\partial {N}_{u2}}{\partial x} & \cdots & 0 & 0 & \cdots \\ 0 & 0 & \cdots & \tfrac{\partial {N}_{u1}}{\partial y} & \tfrac{\partial {N}_{u2}}{\partial y} & \cdots \\ 0 & 0 & \cdots & 0 & 0 & \cdots \\ \tfrac{\partial {N}_{u1}}{\partial y} & \tfrac{\partial {N}_{u2}}{\partial y} & \cdots & \tfrac{\partial {N}_{u1}}{\partial x} & \tfrac{\partial {N}_{u2}}{\partial x} & \cdots \end{array}\right]\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}{{\boldsymbol{B}}}_{u}^{* }=\left[\begin{array}{cccccc}\tfrac{{\partial }^{2}{N}_{u1}}{\partial {x}^{2}} & \tfrac{{\partial }^{2}{N}_{u2}}{\partial {x}^{2}} & \cdots & \displaystyle \frac{{\partial }^{2}{N}_{u1}}{\partial x\partial y} & \tfrac{{\partial }^{2}{N}_{u2}}{\partial x\partial y} & \cdots \\ \tfrac{{\partial }^{2}{N}_{u1}}{\partial x\partial y} & \tfrac{{\partial }^{2}{N}_{u2}}{\partial x\partial y} & \cdots & \tfrac{{\partial }^{2}{N}_{u1}}{\partial {y}^{2}} & \tfrac{{\partial }^{2}{N}_{u2}}{\partial {y}^{2}} & \cdots \end{array}\right]\end{eqnarray} \tag{ 41 }$$

such that the strains are given by ε = B _u [u_x ; u_y ] and the hydrostatic stress gradient as ${\boldsymbol{\nabla }}{\sigma }_{H}=E/(3(1-2\nu )){{\boldsymbol{B}}}_{u}^{* }[{{\bf{u}}}_{x};\,{{\bf{u}}}_{y}]$.

For the electrolyte, the discretised weak forms of the mass balances, Eq. 1, and electroneutrality condition, Eq. 2, are given by:

$$\begin{eqnarray}\begin{array}{rcl}{{\bf{f}}}_{c\pi } & = & {\int }_{{{\rm{\Omega }}}_{e}}\frac{1}{{\rm{\Delta }}t}{{\bf{N}}}_{C}^{T}{{\bf{N}}}_{C}\left({{\bf{C}}}_{\pi }^{t+{\rm{\Delta }}t}-{{\bf{C}}}_{\pi }^{t}\right){\rm{d}}{{\rm{\Omega }}}_{e}\\ & & +{\int }_{{{\rm{\Omega }}}_{e}}{D}_{\pi }{\left({\boldsymbol{\nabla }}{{\bf{N}}}_{c}\right)}^{T}{\boldsymbol{\nabla }}{{\bf{N}}}_{C}{{\bf{C}}}_{\pi }^{t+{\rm{\Delta }}t}\,{\rm{d}}{{\rm{\Omega }}}_{e}\\ & & +{\int }_{{{\rm{\Omega }}}_{e}}\frac{{D}_{\pi }{z}_{\pi }F}{{RT}}{\left({\boldsymbol{\nabla }}{{\bf{N}}}_{c}\right)}^{T}\left({{\bf{N}}}_{C}{{\bf{C}}}_{\pi }^{t+{\rm{\Delta }}t}\right){\boldsymbol{\nabla }}{{\bf{N}}}_{\varphi }{{\boldsymbol{\varphi }}}^{t+{\rm{\Delta }}t}\,{\rm{d}}{{\rm{\Omega }}}_{e}+{{\bf{R}}}_{\pi }+{{\boldsymbol{\nu }}}_{\pi }\\ & & -{\int }_{{{\rm{\Gamma }}}_{e}}{{\bf{N}}}_{c}{j}_{c}\,{\rm{d}}{{\rm{\Gamma }}}_{e}={\bf{0}}\end{array}\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}&&{{\boldsymbol{f}}}_{\varphi }={\int }_{{{\rm{\Omega }}}_{e}}\displaystyle \sum _{\pi }{z}_{\pi }{{\bf{N}}}_{\varphi }^{T}{{\bf{N}}}_{c}{{\bf{C}}}_{\pi }^{t+{\rm{\Delta }}t}\,{\rm{d}}{{\rm{\Omega }}}_{e}={\bf{0}}\end{eqnarray} \tag{ 43 }$$

where the bulk reaction rates are given by:

$$\begin{eqnarray}\begin{array}{ccc}{{\bf{R}}}_{{{\rm{H}}}^{+}} & = & -{\int }_{{{\rm{\Omega }}}_{e}}{k}_{{eq}}{{\bf{N}}}_{c}^{T}\left({K}_{w}-\left({{\bf{N}}}_{c}{{\bf{C}}}_{{{\rm{H}}}^{+}}^{t+{\rm{\Delta }}t}\right)\left({{\bf{N}}}_{c}{{\bf{C}}}_{{\mathrm{OH}}^{-}}^{t+{\rm{\Delta }}t}\right)\right){\rm{d}}{{\rm{\Omega }}}_{e}\\ & & -{\int }_{{{\rm{\Omega }}}_{e}}{{\bf{N}}}_{c}^{T}\left({k}_{{fe}}{{\bf{N}}}_{c}{{\bf{C}}}_{{\mathrm{Fe}}^{2+}}^{t+{\rm{\Delta }}t}\right.\\ & & \left.-{{\bf{N}}}_{c}{{\bf{C}}}_{{\mathrm{FeOH}}^{+}}^{t+{\rm{\Delta }}t}\left({k}_{{fe}}^{^{\prime} }{{\bf{N}}}_{c}{{\bf{C}}}_{{{\rm{H}}}^{+}}^{t+{\rm{\Delta }}t}-{k}_{{feoh}}\right)\right){\rm{d}}{{\rm{\Omega }}}_{e}\end{array}\end{eqnarray} \tag{ 44a }$$

$$\begin{eqnarray}&&{{\bf{R}}}_{{\mathrm{OH}}^{-}}=-{\int }_{{{\rm{\Omega }}}_{e}}{{\bf{N}}}_{c}^{T}{k}_{{eq}}\left({K}_{w}-\left({{\bf{N}}}_{c}{{\bf{C}}}_{{{\rm{H}}}^{+}}^{t+{\rm{\Delta }}t}\right)\left({{\bf{N}}}_{c}{{\bf{C}}}_{{\mathrm{OH}}^{-}}^{t+{\rm{\Delta }}t}\right)\right){\rm{d}}{{\rm{\Omega }}}_{e}\end{eqnarray} \tag{ 44b }$$

$$\begin{eqnarray}&&{{\bf{R}}}_{{\mathrm{Fe}}^{2+}}={\int }_{{{\rm{\Omega }}}_{e}}{{\bf{N}}}_{c}^{T}\left({k}_{{fe}}{{\bf{N}}}_{c}{{\bf{C}}}_{{\mathrm{Fe}}^{2+}}-{k}_{{fe}}^{{\prime} }{{\bf{N}}}_{c}{{\bf{C}}}_{{\mathrm{FeOH}}^{+}}^{t+{\rm{\Delta }}t}{{\bf{N}}}_{C}{{\bf{C}}}_{{{\rm{H}}}^{+}}^{t+{\rm{\Delta }}t}\right)\,{\rm{d}}{{\rm{\Omega }}}_{e}\end{eqnarray} \tag{ 44c }$$

$$\begin{eqnarray}&&\begin{array}{c}{{\bf{R}}}_{{\mathrm{FeOH}}^{+}}=-{\int }_{{{\rm{\Omega }}}_{e}}{{\bf{N}}}_{c}^{T}\left({k}_{{fe}}{{\bf{N}}}_{c}{{\bf{C}}}_{{\mathrm{Fe}}^{2+}}^{t+{\rm{\Delta }}t}\right.\\ \left.-{{\bf{N}}}_{c}{{\bf{C}}}_{{\mathrm{FeOH}}^{+}}^{t+{\rm{\Delta }}t}\,\left(k{{\prime} }_{{fe}}{{\bf{N}}}_{c}{{\bf{C}}}_{{{\rm{H}}}^{+}}^{t+{\rm{\Delta }}t}+{k}_{{feoh}}\right)\right){\rm{d}}{{\rm{\Omega }}}_{e}\end{array}\end{eqnarray} \tag{ 44d }$$

$$\begin{eqnarray}&&{{\bf{R}}}_{{\mathrm{Na}}^{+}}={{\bf{R}}}_{{\mathrm{Cl}}^{-}}={\bf{0}}\end{eqnarray} \tag{ 44e }$$

and the surface reaction terms by:

$$\begin{eqnarray}&&{{\boldsymbol{\nu }}}_{{{\rm{H}}}^{+}}={\int }_{{{\rm{\Gamma }}}_{{int}}}{{\bf{N}}}_{C}^{T}\left({\nu }_{{Va}}^{t+{\rm{\Delta }}t}-\nu {{\prime} }_{{Va}}^{t+{\rm{\Delta }}t}+{\nu }_{{Ha}}^{t+{\rm{\Delta }}t}\right){\rm{d}}{{\rm{\Gamma }}}_{{int}}\end{eqnarray} \tag{ 45a }$$

$$\begin{eqnarray}&&{{\boldsymbol{\nu }}}_{{\mathrm{OH}}^{-}}=-{\int }_{{{\rm{\Gamma }}}_{{int}}}{{\bf{N}}}_{C}^{T}\left({\nu }_{{Vb}}^{t+{\rm{\Delta }}t}-\nu {{\prime} }_{{Vb}}^{t+{\rm{\Delta }}t}+{\nu }_{{Hb}}^{t+{\rm{\Delta }}t}\right){\rm{d}}{{\rm{\Gamma }}}_{{int}}\end{eqnarray} \tag{ 45b }$$

$$\begin{eqnarray}&&{{\boldsymbol{\nu }}}_{{\mathrm{Fe}}^{2+}}={\int }_{{{\rm{\Gamma }}}_{{int}}}{{\bf{N}}}_{C}^{T}\left({\nu }_{c}^{t+{\rm{\Delta }}t}-\nu {{\prime} }_{c}^{t+{\rm{\Delta }}t}\right){\rm{d}}{{\rm{\Gamma }}}_{{int}}\end{eqnarray} \tag{ 45c }$$

$$\begin{eqnarray}&&{{\boldsymbol{\nu }}}_{{\mathrm{Na}}^{+}}={\nu }_{{\mathrm{Cl}}^{-}}={\nu }_{{\mathrm{FeOH}}^{+}}={\bf{0}}\end{eqnarray} \tag{ 45d }$$

Algorithm 1. Overview of Gauss and lumped integration schemes

1: Start of element assembly

2: Set element flux vector and lumped weight vector to 0: ${\bf{q}}={\bf{W}}={\bf{0}}$

3: for integration point ${ip}=1:{n}_{{ip}}$ do

4: Calculate Gauss weights wip , and shape functions ${{\boldsymbol{N}}}_{{ip}}$

5: Calculate reaction rates ${\nu }_{{ip}}$ based on local integration point values (${C}_{{{\rm{H}}}^{+}}={{\bf{N}}}_{{ip}}{{\bf{C}}}_{{{\rm{H}}}^{+}}$, etc.)

6: Perform non-lumped integrations: ${\bf{q}}={\bf{q}}+{w}_{{ip}}{{\bf{N}}}^{T}{\nu }_{{ip}}$

7: Calculate lumped weights: ${\bf{W}}={\bf{W}}+{w}_{{ip}}{{\bf{N}}}_{{ip}}$

8: end for

9: for node ${nd}=1:{n}_{{nodes}}$ do

10: Calculate reaction rates ${\nu }_{{nd}}$ based on nodal values (using ${{\bf{C}}}_{{{\rm{H}}}^{+}}({nd})$, etc.)

11: Perform lumped integrations: ${\bf{q}}({nd})={\bf{q}}({nd})+{\bf{W}}({nd}){\nu }_{{nd}}$

12: end for

13: Go to next element

Finally, we have the interfacial mass balance for the adsorbed hydrogen, which is given in discretised form as:

$$\begin{eqnarray}\begin{array}{rcl}{{\bf{f}}}_{\theta } & = & {\int }_{{{\rm{\Gamma }}}_{{int}}}\frac{{N}_{{ads}}}{{\rm{\Delta }}t}{{\bf{N}}}_{\theta }^{T}{{\bf{N}}}_{\theta }\left({{\boldsymbol{\theta }}}^{t+{\rm{\Delta }}t}-{{\boldsymbol{\theta }}}^{t}\right){\rm{d}}{{\rm{\Gamma }}}_{{int}}\\ & - & {\int }_{{{\rm{\Gamma }}}_{{int}}}{{\bf{N}}}_{\theta }^{T}\left(({\nu }_{{Va}}-{\nu }_{{Va}}^{^{\prime} })-{\nu }_{{Ha}}-2{\nu }_{T}-({\nu }_{A}-{\nu }_{A}^{^{\prime} })\right.\\ & & \left.+({\nu }_{{Vb}}-{\nu }_{{Vb}}^{^{\prime} })-{\nu }_{{Hb}}\right){\rm{d}}{{\rm{\Gamma }}}_{{int}}=0\end{array}\end{eqnarray} \tag{ 46 }$$

The mass balances from Eqs. 39, 42 and 46, the momentum balance from Eq. 38, and the charge balance from Eq. 43 together fully describe the behavior of the electrolyte, metal and interface. To eliminate any stability issues originating from the coupling between the domains, both domains are solved at the same time in a monolithic manner. To solve the resulting system of equations, an iterative Newton-Raphson scheme is used. Details of the solution scheme are given in Appendix A, together with the formulation of the tangential matrices employed.

In the context of the finite element method, the weak form equations are typically integrated using a Gauss integration scheme; evaluating the value of the functions at set points within the element and using these values to estimate the integral. However, as discussed below, this results in spurious oscillations throughout the electrolyte. These oscillations originate from two main sources. One is the water auto-ionisation reaction—the first term in Eq. 44a and Eq. 44b. More specifically, oscillations are a consequence of the high value that must be assigned to the penalty reaction rate term k_eq to ensure equilibrium. A second source is the absorption reaction, as a result of the high magnitudes of the forward and backward reaction rates, k_A and ${k}_{A}^{{\prime} }$. These are relevant for the solution of the lattice hydrogen concentration C_L , see Eq. 39, and the surface coverage, see Eq. 46, bringing instabilities to the coupling between these two fields.

To prevent these oscillations, we use a lumped integration scheme. ^21,22 This scheme first constructs the lumped weight vector by employing a standard Gauss integration scheme on the element under consideration as:

$$\begin{eqnarray}&&{\bf{W}}={\int }_{{{\rm{\Omega }}}_{{el}}}{{\bf{N}}}^{T}\,{\rm{d}}{{\rm{\Omega }}}_{{el}}=\displaystyle \sum _{{ip}}{w}_{{ip}}{{\bf{N}}}^{T}({{\bf{x}}}_{{ip}})\end{eqnarray} \tag{ 47 }$$

where w_ip indicates the weighting factor of the integration point, and x_ip are the coordinates of the current integration point. The nodal integration weights are then used to integrate the reactions based on the values in the nodes. For example, the reaction rate of OH⁻, Eq. 44b, is reformulated as:

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{R}}}_{{{OH}}^{-}} & = & -{\displaystyle \int }_{{{\rm{\Omega }}}_{e}}{{\bf{N}}}_{c}^{T}{k}_{{eq}}\left({K}_{w}-\left({{\bf{N}}}_{c}{{\bf{C}}}_{{{\rm{H}}}^{+}}^{t+{\rm{\Delta }}t}\right)\left({{\bf{N}}}_{c}{{\bf{C}}}_{{\mathrm{OH}}^{-}}^{t+{\rm{\Delta }}t}\right)\right){\rm{d}}{{\rm{\Omega }}}_{e}\\ & = & -\displaystyle \sum _{{el}}\displaystyle \sum _{{nd}}{\bf{W}}({nd}){k}_{{eq}}\left({K}_{W}-{{\bf{C}}}_{{{\rm{H}}}^{+}}({nd}){{\bf{C}}}_{{\mathrm{OH}}^{-}}({nd})\right)\end{array}\end{eqnarray} \tag{ 48 }$$

and similarly for the water auto-ionisation term, Eq. 44a. Here, □(nd) is used to indicate the nodal component of vector □, and the resulting magnitude is allocated to the force vector index associated with the OH⁻ degree of freedom at node nd. This force vector assembly process is summarised in Algorithm 1 and the MATLAB finite element code developed is openly shared to facilitate understanding and uptake. ¹

Applying a similar scheme to the absorption reaction terms in Eqs. 39 and 46 leads to the following sum over all the nodes in the elements at the interface:

$$\begin{eqnarray}&&\begin{array}{c}{\int }_{{{\rm{\Gamma }}}_{{int}}}{{\bf{N}}}_{x}^{T}({\nu }_{{Va}}-{\nu }_{{Va}}^{{\prime} }){\rm{d}}{{\rm{\Gamma }}}_{{int}}=\displaystyle \sum _{{iel}}\displaystyle \sum _{{nd}}{\boldsymbol{W}}({nd})\left({k}_{A}({N}_{L}-{{\bf{C}}}_{L}({nd}))\right.\\ \left.\cdot \,{{\boldsymbol{\theta }}}_{{ads}}({nd})-{k}_{A}^{^{\prime} }{{\bf{C}}}_{L}({nd})(1-{{\boldsymbol{\theta }}}_{{ads}}({nd}))\right)\end{array}\end{eqnarray} \tag{ 49 }$$

where N_x represents the interpolation functions N_θ for the terms associated with the surface coverage, and N_L for the terms corresponding to the interstitial lattice hydrogen concentration. As a result, all off-diagonal terms associated with these reactions are taken out from the stiffness matrix. In this way, the interactions between neighbouring nodes that take place through these high reaction rates are eliminated.

Let us proceed to illustrate the effect of the lumped integration scheme proposed by considering a simple example where only the hydrogen absorption reaction is modelled over a single element. In such scenario, the weak forms read:

$$\begin{eqnarray}\begin{array}{ccc}{{\boldsymbol{q}}}_{\theta } & = & \int {{\bf{N}}}^{T}\left({k}_{A}\left({N}_{L}-{{\bf{N}}}_{L}{{\bf{C}}}_{L}\right){{\bf{N}}}_{\theta }{\boldsymbol{\theta }}\right.\\ & & \left.-{k}_{A}^{^{\prime} }{{\bf{N}}}_{L}{{\bf{C}}}_{L}\left(1-{{\bf{N}}}_{\theta }{\boldsymbol{\theta }}\right)\right)\,{\rm{d}}{\rm{\Gamma }}={\bf{0}}\end{array}\end{eqnarray} \tag{ 50 }$$

$$\begin{eqnarray}\begin{array}{ccc}{{\boldsymbol{q}}}_{L} & = & -\int {{\bf{N}}}^{T}\left({k}_{A}\left({N}_{L}-{{\bf{N}}}_{L}{{\bf{C}}}_{L}\right){{\bf{N}}}_{\theta }{\boldsymbol{\theta }}\right.\\ & & \left.-{k}_{A}^{^{\prime} }{{\bf{N}}}_{L}{{\bf{C}}}_{L}\left(1-{{\bf{N}}}_{\theta }{\boldsymbol{\theta }}\right)\right)\,{\rm{d}}{\rm{\Gamma }}={\bf{0}}\end{array}\end{eqnarray} \tag{ 51 }$$

where q is used to denote the weak forms solely related to the reaction fluxes. We shall now calculate the tangential matrix terms for this system using Gauss and lumped integration. To this end, a single quadratic element is used, the solution vectors are defined as C_L = θ = 0, and we set ${k}_{A}={k}_{A}^{{\prime} }=3$ and N_L = 1, rendering:

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{K}}}_{{Gauss}} & = & \left[\begin{array}{cccccc}0.6 & 0.3 & 0.1 & -0.6 & -0.3 & -0.1\\ 0.3 & 0.4 & 0.3 & -0.3 & -0.4 & -0.3\\ 0.1 & 0.3 & 0.6 & -0.1 & -0.3 & -0.6\\ -0.6 & -0.3 & -0.1 & 0.6 & 0.3 & 0.1\\ -0.3 & -0.4 & -0.3 & 0.3 & 0.4 & 0.3\\ -0.1 & -0.3 & -0.6 & 0.1 & 0.3 & 0.6\end{array}\right]\\ {{\boldsymbol{K}}}_{{Lumped}} & = & \left[\begin{array}{cccccc}1 & 0 & 0 & -1 & 0 & 0\\ 0 & 1 & 0 & 0 & -1 & 0\\ 0 & 0 & 1 & 0 & 0 & -1\\ -1 & 0 & 0 & 1 & 0 & 0\\ 0 & -1 & 0 & 0 & 1 & 0\\ 0 & 0 & -1 & 0 & 0 & 1\end{array}\right]\end{array}\end{eqnarray} \tag{ 52 }$$

It can be seen that both integration schemes result in a similar hydrogen flux (e.g., for each row, the same magnitude is obtained when adding the terms in columns 1-3). As expected, the lumped scheme results in a tangential matrix with all terms concentrated on the diagonals of each degree-of-freedom sub-matrix. This is in sharp contrast with the Gauss scheme, which results in a dense matrix. Differences are also showcased in Fig. 2, which shows the eigenvalues and eigenmodes obtained from the matrices in Eq. 52. These eigenmodes indicate the manner in which concentrations are expected to change, and the eigenvalues provide an indication of the frequency response of these changes. Since capacity terms are not included, three "free body motion" eigenvalues are present, representing the modes in which both the lattice and adsorbed hydrogen concentrations increase equally. Looking at the non-zero eigenvalues and accompanying eigenmodes, one can see clear differences between the lumped and Gauss integration schemes. The lumped scheme obtains three equal eigenvalues, corresponding to a single set of nodes transferring hydrogen from the surface to the metal lattice. As all eigenmodes are equal, this system response will allow for an oscillation-free transfer of hydrogen from the metal surface into the interstitial lattice sites. In contrast, Gaussian integration results in three different eigenvalues. The lower two of these will produce oscillations due to the way the three nodes are coupled. Furthermore, they do not transfer hydrogen from the surface to the lattice, but instead redistribute it between neighbouring nodes. Only the inclusion of the highest eigenmode allows for the transfer of hydrogen. These low eigenvalues corresponding to oscillatory results indicate that the hydrogen absorption reaction is prone to spurious oscillations when this term becomes dominant; this is not observed for the lumped scheme.

Zoom In Zoom Out Reset image size

We proceed to demonstrate the potential and robustness of the lumped integration scheme presented above by simulating the uptake of hydrogen in the electrolyte-metal domain shown in Fig. 1. First, the differences with Gauss integration are highlighted by examining the results obtained after a single iteration. Then, simulations are conducted for technologically-relevant ranges of material and environmental parameters, showcasing the ability of the lumped integration scheme-based framework to deliver predictions over previously unexplored conditions and time scales. In all cases, the electrolyte and metal sub-domains have 10 × 10 mm dimensions, and a cracked region of 5 × 0.4 mm is considered. Unless otherwise stated, the material and diffusion parameters are given in Table I, while the reaction constants are listed in Table II. As initial condition, we use an electrolyte with pH = 5 and initial concentrations ${C}_{{{\rm{H}}}^{+}}={10}^{-2}\,\mathrm{mol}/{{\rm{m}}}^{3}$, ${C}_{{\mathrm{OH}}^{-}}={10}^{-6}\,\mathrm{mol}\,{{\rm{m}}}^{-3}$, ${C}_{{\mathrm{Na}}^{+}}=599.99\,\mathrm{mol}\,{{\rm{m}}}^{-3}$, ${C}_{{\mathrm{Cl}}^{-}}=6\cdot {10}^{2}\,\mathrm{mol}\,{{\rm{m}}}^{-3}$, and ${C}_{{\mathrm{Fe}}^{2+}}={C}_{{\mathrm{FeOH}}^{+}}=0\,\mathrm{mol}\,{{\rm{m}}}^{-3}$. Together with $\overline{\varphi }=0\,{{\rm{V}}}_{\mathrm{SHE}}$, these concentrations are also prescribed on the left edge as boundary conditions throughout the simulation. For the metal sub-domain, C_L =0 and u = 0 are assumed at t = 0 and a constant vertical displacement of U_ext = 10 μm is prescribed at the top, while constraining both the horizontal and vertical displacements at the bottom. Finally, the only source of lattice hydrogen is through the metal-electrolyte interface, with all other boundaries not allowing any hydrogen flux. For the spatial discretisation, we use a finite element mesh composed of quadratic, triangular Bernstein elements with a characteristic element length of 0.1 mm near the metal-electrolyte interface and within the crack, increasing up to 0.5 mm further away. This results in a total of 3·10⁴ degrees of freedom (DOFs).

Table I. Material and ionic transport parameters used in all cases reported within this paper.

Parameter		Value
Young's Modulus	E	200 GPa
Poisson ratio	ν	0.3
H+ diffusion coefficient	${D}_{{{\rm{H}}}^{+}}$	9.3 · 10−9 m2 s−1
OH− diffusion coefficient	${D}_{{\mathrm{OH}}^{-}}$	5.3 · 10−9 m2 s−1
Na+ diffusion coefficient	${D}_{{\mathrm{Na}}^{+}}$	1.3 · 10−9 m2 s−1
Cl− diffusion coefficient	${D}_{{\mathrm{Cl}}^{-}}$	2 · 10−9 m2 s−1
Fe2+ diffusion coefficient	${D}_{{\mathrm{Fe}}^{2+}}$	1.4 · 10−9 m2 s−1
FeOH+ diffusion coefficient	${D}_{{\mathrm{FeOH}}^{+}}$	10−9 m2 s−1
Partial molar volume	${\overline{V}}_{H}$	2 · 10−6 mol m−3
Surface adsorption sites	Nads	10−3 mol m−2
Lattice sites	NL	106 mol m−3
Lattice diffusion coefficient	DL	10−9 m2 s−1
Temperature	T	293.15 K

Table II. Reaction rate constants used throughout the paper (backward k and forward $k^{\prime}$. For the cases in which a parametric sweep is performed, all reaction constants are taken from this table except those explicitly stated to be different.

Reaction	k	$k^{\prime}$	α	Eeq
νVa	1 · 10−4 m s−1	$1\cdot {10}^{-10}\,\mathrm{mol}{({{\rm{m}}}^{2}{\rm{s}})}^{-1}$	0.5	0 VSHE
νHa	1 · 10−10 m s−1	$0\,\mathrm{mol}{({{\rm{m}}}^{2}\mathrm{Pa}\,\ {\rm{s}})}^{-1}$	0.3	0 VSHE
νT	$1\cdot {10}^{-6}\,\mathrm{mol}{({{\rm{m}}}^{2}{\rm{s}})}^{-1}$	$0\,\mathrm{mol}{({{\rm{m}}}^{2}{\rm{s}}\ \,{\mathrm{Pa}}^{1/2})}^{-1}$	—	—
νA	1 · 103 m s−1	7 · 107 m s−1	—	—
νVb	$1\cdot {10}^{-8}\,\mathrm{mol}{({{\rm{m}}}^{2}{\rm{s}})}^{-1}$	1 · 10−13 m s−1	0.5	0 VSHE
νHb	$1\cdot {10}^{-10}\,\mathrm{mol}{({{\rm{m}}}^{2}{\rm{s}})}^{-1}$	0 m(Pa s)−1	0.3	0 VSHE
νFe	$1.5\cdot {10}^{-10}\,\mathrm{mol}{({{\rm{m}}}^{2}{\rm{s}})}^{-1}$	1.5 · 10−10 m s−1	0.5	−0.4 VSHE
kfe	0.1 s	10−3 m3(mol s)−1
kfeoh	10−3 s−1
keq	106 m3(mol s)−1

Effectiveness of lumped integration

To showcase the effect of lumped integration, results are obtained after one single Newton-Raphson iteration and compared to those obtained using Gauss integration. More specifically, simulations are performed using no lumped integration, lumped integration for solely the absorption reaction, and lumped integration for both the absorption and the auto-ionisation reactions. The results obtained are shown in Fig. 3 for the surface coverage and in Fig. 4 for the OH⁻ concentration. Consider first the results obtained for the surface coverage, Fig. 3. Remarkable differences are observed between the results obtained with Gaussian and lumped integration, with the former showing severe oscillations and surface coverage values that are orders of magnitude off from the expected solution. While these are still non-converged results, and thus are not expected to be correct after a single iteration, the oscillations observed are indicative of convergence issues. However, when lumped integration is used, these oscillations disappear and the results obtained are in agreement with expectations after a single iteration: Attaining the right order of magnitude, solely positive surface coverage and concentrations, and tending toward the fully converged solution. Similar conclusions can be drawn from the OH⁻ concentration results, Fig. 4. Using Gauss integration and only using lumped integration for the absorption reaction results in oscillations an order of magnitude higher than the correct results, while using lumped integration for both reaction terms results in vastly improved results containing only minor oscillations. It should be noted that the removal of these oscillations does not guarantee an improved convergence radius and rate, but it is a good indication of such improvements. As such, these results emphasise the ability of the lumped integration scheme to enable the use of larger time increments and quantify hydrogen uptake over relevant parameter ranges without encountering convergence issues.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Predictions across environmental and material conditions

Converged results are subsequently obtained to demonstrate the ability of the lumped integration-based framework to map the parameter space and gain insight into hydrogen absorption behavior over large time scales. All simulations are performed using an initial time increment of Δt = 30 s, increasing by 5% per new time step until the full simulation duration of 50 years is reached after a total of 303 time steps. Within these 50 years, all simulated cases achieve their steady-state solution, in which the hydrogen contents become constant. All results in this Section presented as steady-state results are therefore the result of performing the time-dependent simulation for the full 50 year duration. By starting with a smaller time step, we accurately capture short-term behavior (days) while exploiting the ever-increasing time-step size to efficiently obtain results for all relevant time scales in hydrogen embrittlement.More specifically, these time increments allow the rapid changes in pH and ionic concentrations to be captured, which occur on the order of minutes at the onset of the simulations. Increasing the increment afterwards allows modelling as well the slow hydrogen absorption timescales. A representative result from the converged simulations is shown in Fig. 5, where contours of pH are shown on the electrolyte sub-domain while contours of lattice hydrogen concentration are presented on the metal sub-domain. The metal absorbs hydrogen from the neighbouring electrolyte, increasing the lattice hydrogen content in the material while decreasing the amount of hydrogen ions within the electrolyte and thus the pH. Within the metal, the largest concentration of hydrogen is located around the crack tip, having diffused to this location due to the large hydrostatic stresses present in this region. This boundary value problem was simulated in Ref. 20 using the commercial finite element package COMSOL. However, due to the described stability and oscillation issues inherent to Gaussian integration, COMSOL-based simulations were restricted to applied potentials equal to E_m = −0.7 V_SHE or higher, while here results are shown for E_m = −1 V_SHE, a regime of relevance for hydrogen embrittlement and cathodic protection. Furthermore, the time step size restrictions required for stability when using Gauss integration (with time increments between Δt = 10⁻⁵ s and Δt = 10⁻² s) imposed a constraint on the time scales addressed. In particular, 82 h were needed to simulate 10 min of hydrogen uptake in an Intel i7-10700 CPU. In contrast, the described lumped integration scheme, not suffering from these limitations, is able to deliver predictions over a time scale of 50 years within 3 h, using the same system.

Zoom In Zoom Out Reset image size

Taking advantage of the robustness and stability of the lumped integration implementation, we proceed to obtain results in regimes previously unexplored yet critical for hydrogen embrittlement predictions. To this end, the absorbed hydrogen is quantified using the average and maximum lattice hydrogen concentrations, ${\overline{C}}_{L}$ and ${C}_{L}^{\max }$ respectively. The maximum concentration is always located at the crack tip, due to stress localization, whereas the average is obtained by numerically integrating the lattice hydrogen over the complete metal subdomain, and normalising it with the metal surface area. These two quantities give a good indication of the total behavior of the system; the average hydrogen concentration provides an estimate of overall hydrogen uptake, and the maximum hydrogen concentration indicates how unevenly hydrogen is distributed due to the presence of the stresses within the metal.

Application across the reaction constant space

We begin to explore the reaction constant space by varying the values of k_A and ${k}_{A}^{\prime}$ within the experimentally reported range. All other reaction rate constants are kept constant and take the values given in Table II. The experimental literature reports values for k_A spanning the range k_A = 2.4 · 10⁻¹² m s⁻¹⁶¹ to k_A = 1.2 · 10⁵ m s⁻¹. ¹⁸ Accordingly, simulations are conducted within the range k_A = [10⁻¹⁴, 10⁵] m s⁻¹. The backward reaction constant is chosen such that ${k}_{A}^{{\prime} }/{k}_{A}=7\cdot {10}^{4}$, altering the rate at which the two reactions occur but not their equilibrium. These simulations are performed at four applied potentials, E_m = −1 V_SHE, E_m = −0.5 V_SHE, E_m = 0 V_SHE, and E_m = 0.5 V_SHE, going from a hydrogen-dominated regime to a corrosion reaction dominated one. The resulting hydrogen uptake within the metal is shown in Fig. 6 using E_m = −1 V_SHE. The results are only shown up to k_A = 10⁻⁷ m s⁻¹, as higher absorption reaction constants led to almost identical predictions. Changes in the absorption reaction constant are seen to strongly affect the rate at which the lattice hydrogen concentration achieves equilibrium with the surrounding electrolyte. For unrealistically low values of k_A (∼10⁻¹⁴), very little hydrogen enters the metal within the simulated 50 years. In contrast, values of k_A = 10⁻⁹ m s⁻¹ and higher lead to the same average lattice hydrogen concentrations, indicating that the hydrogen absorption is almost instantaneous, and limited either by the diffusion within the metal or by the other surface reactions. Results obtained using positive metal potentials are shown in Fig. 6b using E_m = 0.5 V_SHE. Hydrogen uptake is hindered by the inhibiting effect of the metal potential, leading to a lower level of absorbed hydrogen and thus a shorter time needed to achieve steady state. No sensitivity to the absorption rate constant is observed for high k_A values. However, as k_A takes lower values, small differences are observed in the magnitude of the lattice hydrogen concentration at steady state. This points to a relation between the equilibrium hydrogen surface coverage and the absorption rate constant for low k_A values, as the Heyrovsky and Tafel reactions become dominant and remove hydrogen faster than it is absorbed, resulting in a lower degree of hydrogen uptake.

Zoom In Zoom Out Reset image size

The role of the Volmer and Heyrovsky reaction constants is subsequently investigated. Two reaction rate constants (k₁, k₂) are introduced to facilitate interpretation of the results and ensure consistency across the acid and basic regimes—reaction rate constants are chosen to ensure that acidic reactions are dominant below pH = 7 while basic reactions dominate above pH = 7. The reaction rate constant k₁ alters the Volmer reaction rates, while k₂ varies the Heyrovsky rate. From these two constants, we construct the reaction constants used in our model as: k_Va = k₁, ${k}_{{Va}}^{{\prime} }={k}_{1}\cdot {10}^{-6}$, k_Vb = k₁ · 10⁻⁴, ${k}_{{Vb}}^{{\prime} }={k}_{1}\cdot {10}^{-9}$, k_Ha = k₂, and k_Hb = k₂ · 10⁻⁴. The values for k₁ are varied between 10⁻¹⁰ and 1, while k₂ is varied between 10⁻¹⁰ and 10⁻², covering the complete range of values reported in literature. ²⁰ The average lattice hydrogen concentrations obtained in the metal while varying these two parameters are shown in Fig. 7. Two applied potentials are considered, E_m = −1 V_SHE and E_m = 0 V_SHE. For both cases, very low k₁ values result in negligible hydrogen uptake, independently of the value of k₂. For increased rates for the Volmer reaction (higher k₁), more hydrogen is present within the metal. This effect is strongest for the negative metal potential, where the Volmer reaction is further accelerated due to the role of the electric overpotential. In contrast, increasing the Heyrovsky reaction rates (higher k₂) reduces the amount of hydrogen available. As a result, the metal is nearly saturated when high Volmer reaction constants are combined with low Heyrovsky reaction constants. The time required to attain a hydrogen content that is 90% of the steady state value is shown in Fig. 8 as a function of k₁ and k₂. For the E_m = −1 V_SHE simulations, the longest required times correspond to the highest hydrogen contents, requiring close to 15 years to attain these steady state values for k₁ ≈ 10⁻⁵ and low k₂ values. When the Volmer reaction rates are decreased by lowering k₁, the time required to reach steady state is strongly reduced to only several days. In contrast, increasing the Volmer rates beyond k₁ = 10⁻⁵ brings in a reduction in the time needed to achieve steady state. This corresponds with the reaction rates at which the metal becomes fully saturated, and thus all additional hydrogen produced at the surface is no longer used to increase the lattice concentration, but rather spent on reaching the equilibrium state. For the E_m = 0 V_SHE simulations, less hydrogen is present in the metal lattice at equilibrium, and as such less time is required to attain this equilibrium.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

From the results presented in this section, it can be concluded that the presented scheme is stable, oscillation-free, and well-converging across the whole range of relevant reaction constants. Furthermore, it is interesting to note that the range of values reported in the literature causes a significant spread in the results, ranging from close to no lattice hydrogen to a fully saturated metal.

Application across the metal potential space

Next, we investigate the influence of the applied potential on hydrogen uptake, keeping the reaction rate parameters constant. The enhanced stability provided by the lumped integration scheme presented here enables obtaining results over the range E_m = [−1.5, 1] V_SHE, going from strong cathodic protection conditions to the corrosion dominated regime. Previous reported results were limited to the range E_m = [−0.7, 0.5] V_SHE, due to convergence issues resulting from the above discussed stability problems. ²⁰ The evolution in time of the average content of hydrogen within the metal is shown in Fig. 9. When the applied potential is −1.5 V_SHE, the metal becomes fully saturated within two days. Strong negative potentials accelerate surface reaction, leading to rapid saturation of the surface sites and the metal bulk sites near the surface. This, in turn, enhances hydrogen diffusivity within the metal through the non-linear diffusion term, see Eq. 13. This high diffusivity facilitates the distribution of hydrogen within the specimen, more readily attaining steady state conditions. On the other hand, increasing the metal potential brings in a reduction in electric overpotential, and thus decreases the forward surface reaction rates while accelerating the backward reaction rates. As a result, for theE_m = −1.3 V_SHE simulations it takes over 20 d to attain steady state, with the lattice concentration at this steady state being slightly lower compared to the E_m = −1.5 V_SHE results. Further increasing the metal potential shows a more pronounced effect, with hydrogen-producing surface reactions becoming slower and the steady state ${\overline{C}}_{L}$ becoming smaller, which results in a shorter time needed to attain this equilibrium condition.

Zoom In Zoom Out Reset image size

The relationship between the lattice hydrogen concentration and the applied potential at steady state is shown in Fig. 10, using a logarithmic scale. Both the average (${\overline{C}}_{L}$) and maximum (${C}_{L}^{\max }$) lattice hydrogen concentrations are provided. This mapping is expected to be useful for the hydrogen assisted fracture community, as it provides a first order approximation for input of chemo-mechanical models of hydrogen embrittlement aiming at delivering predictions over large time scales. The results show that at strongly negative potentials the average and maximum concentrations are closer to each other since the role of stress raisers is reduced, as discussed in the context of Eq. 14. Increasing the metal potential lowers both the average and maximum concentrations. Going into positive metal potentials initially increases the hydrogen uptake through an increase in corrosion rate and subsequent reduction in pH. Increasing the metal potential even further to E_m = 1 V_SHE prevents any of the hydrogen reactions from occurring due to the high electric overpotential, even though the pH of the electrolyte is strongly acidic. This is also seen in Fig. 11, showing the pH and electrolyte potential within the crack. At low metal potentials, the environment becomes highly basic and the electrolyte potential decreases. The opposite happens for high potentials, increasing the metal potential while lowering the pH. Here, one should note that due to the limited size of the domain and the zero electric potential boundary condition on the left boundary, the electrolyte never attains the same potential as the metal. This explains the ever-accelerating effect of decreasing the potential on the hydrogen reactions and the strongly inhibiting effect at high potentials.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Finally, the last case study aims at shedding light into the interplay between applied mechanical stresses and hydrogen uptake, and at demonstrating the performance of the lumped integration scheme presented for relevant regimes of mechanical load. To this end, a tensile rod subjected to a prescribed remote strain is considered, see Fig. 12. The metal domain has a 1 cm radius and 5 cm length, and is contained within an electrolyte domain with an outer radius of 5 cm. The metallic sample contains a defect that starts in the outer surface and penetrates up to a depth of 5 mm. As shown in Fig. 12, the role of the defect geometry is investigated, considering both a rounded notch with radius 0.4 mm and the case of a sharp crack, with an outer opening of 0.8 mm and inner tip radius 0.05 mm. Stress concentrators result in high gradients of hydrostatic stress, which lead to an accumulation of lattice hydrogen in their vicinity. These high lattice hydrogen concentrations shift the equilibrium of the adsorption reaction, increasing the surface coverage which in turn accelerates hydrogen recombination through the Heyrovsky and Tafel reactions. The material properties of the rod and the reaction constants for the electrolyte correspond to those used previously, and are given in Tables I and II. The electrolyte boundary conditions mimic those of the previous analysis; that is, we prescribe at its exterior boundary a constant concentration of the relevant ionic species and a constant electrolyte potential, using identical magnitudes to those reported in the previous study. The displacement is constrained at the bottom of the metal, while the vertical displacement at the top is prescribed based on the imposed average strain, U_ext = _ext · 5 cm. This strain is varied between _ext = 0 and _ext = 5 · 10⁻³, resulting in tensile stresses near the top and bottom of the domain of up to σ_yy = 700 MPa for the used geometries. The electric potential of the metal is varied between −1 V_SHE and 1 V_SHE. The finite element model exploits axial symmetry, allowing the three-dimensional domain to be discretised using a two-dimensional mesh, with this mesh using a minimum element size of 0.1 mm near the metal-electrolyte interface and larger elements with a size of up to 1 mm away from this interface. Near the notch tip, the rounded notch case uses an element with a size of 10 μm, while this characteristic element size equals 1 μm for the sharp crack. These meshes result in a total number of DOFs between 1.6 · 10⁵ (blunted crack) and 2.0 · 10⁵ (sharp crack). Details about the alterations to the previously described scheme due to the axisymmetric nature of the boundary value problem are given in Appendix B.

Zoom In Zoom Out Reset image size

The steady state results for the blunted crack case are given in Fig. 13, in terms of maps for the average lattice concentration ${\overline{C}}_{L}$ (Fig. 13a) and the maximum lattice concentration ${C}_{L}^{\max }$ (Fig. 13b), as a function of the applied potential E_m and the remote strain/stress. The maximum hydrogen attained is highly sensitive to the applied load as crack tip stresses increase with the applied strain and this leads to higher hydrogen contents—see Eq. 14. However, as shown in Fig. 13a, since this effect is localised at the crack tip, the impact on the domain-wide average hydrogen concentration is negligible. The results of Fig. 13a also reveal a ${\overline{C}}_{L}$ sensitive with applied potential that qualitative resembles the findings of Fig. 10. The time required to achieve 90% of the steady state concentration is shown in Figs. 14a and 14b, for the average and maximum lattice hydrogen concentrations, respectively. Since the average concentration is almost insensitive to the imposed strain, so is the time required to attain this average. In contrast, the time required to attain a ${C}_{L}^{\max }$ level that is 90% of that at steady state is dependent on the remote load. When the metal lattice is nearly saturated, increasing the strain facilitates achieving the maximum concentration faster, for instance altering the time from 3 years to 5 d when the metal potential is −1 V_SHE.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The results for the sharp crack case are given in Fig. 15, showing the maximum lattice concentration and time required to obtain this maximum concentration. Due to the increased stresses at the tip of the sharp crack, the hydrostatic stress gradient is much higher for this case compared to the rounded case. As a result, the metal lattice is saturated with hydrogen for less negative metal potentials, and for lower applied strains. Similar to the rounded case, once the metal is saturated, increasing the strains further reduces the time required to obtain the maximum concentration. Due to the larger hydrostatic stress gradient, these maxima are achieved faster, taking less than a day to obtain the equilibrium for the largest strain and lowest potential case. These results not only demonstrate the stability of the presented scheme and its ability to simulate more complex conditions, but also highlight the importance of coupled electro-chemo-mechanical effects in the prediction of hydrogen ingress. Even though the steady state average lattice hydrogen concentration does not depend on the applied loading, the concentration at the tip of the defect depends strongly on the strain level and the defect geometry. Furthermore, the time scale over which this hydrogen enters the metal also depends on these factors, with the time required to achieve steady state conditions ranging from several days to over ten years. The results provided bring new physical insight, given the challenges associated measuring hydrogen content locally near crack tips. While this also hinders direct validation, other outputs of the model (e.g., pH) can be validated against artificial crevice electrochemical cell measurements (see Refs. 20, 62).

Zoom In Zoom Out Reset image size

We have presented a new lumped integration scheme for reducing oscillations in the modelling of metal-electrolyte reactions. The robustness and capabilities of the scheme are demonstrated by addressing the paradigmatic case of hydrogen uptake in metals, a technologically-relevant phenomenon whose prediction is compromised by instabilities resulting from differences in reaction rates. Thus, the lumped integration scheme presented is coupled to a numerical framework resolving electrolyte behavior, interface reactions and hydrogen bulk diffusion in the presence of a mechanical load. The simulation of several boundary value problems of particular interest reveals that:

• The use of conventional Gauss integration schemes results in severe oscillations, decreasing the time increment allowed to retain stable simulations and limiting the range of parameters and time scales that can be simulated.
• The lumped integration scheme presented eliminates numerical oscillations and leads to smooth solutions. Furthermore, lumped integration leads to a significant improvement in tangential matrix conditioning for large time increments and high reaction rates.
• The use of a lumped integration scheme enables the simulation of hydrogen uptake in metals over a wide range of material and environmental parameters. Moreover, it also enables the use of large time increments, which is essential to deliver predictions over technologically-relevant time scales.
• The lumped integration scheme presented is shown to be stable over all relevant values of applied potential and reaction rate constants. Among others, this enables determining the time required to achieve steady state conditions for arbitrary choices of environment-material system.
• The results obtained for the tensile rod case study show that the scheme also remains stable in the presence of varying mechanical load. As a result, steady state times for hydrogen uptake near crack tips have been obtained for the first time.

The scheme presented has shown to make feasible the simulation of hydrogen ingress in a coupled electrolyte-metal domain and could be readily adopted to tackle other relevant electrochemical systems.

The authors gratefully acknowledge financial support through grant EP/V009680/1 ("NEXTGEM") from the Engineering and Physical Sciences Research Council (EPSRC) and the computational resources and support provided by the Imperial College Research Computing Service (http://doi.org/10.14469/hpc/2232). Emilio Martínez-Pañeda additionally acknowledges financial support from UKRI's Future Leaders Fellowship programme [grant MR/V024124/1].

The MATLAB code used to produce the results presented in this paper, together with documentation detailing the use of this code, are made freely available at www.imperial.ac.uk/mechanics-materials/codes and www.empaneda.com. Documentation is also provided, along with example files that enable reproduction of the results shown within the two-dimensional cases.

The discretised equations, Eqs. 38 to 46, are solved using a monolithic Newton-Raphson scheme. Thus, the system of governing equations,

$$\begin{eqnarray}{\left[\begin{array}{ccccc}{{\boldsymbol{K}}}_{{uu}} & {\bf{0}} & {\bf{0}} & {\bf{0}} & {\bf{0}}\\ {{\boldsymbol{K}}}_{{Lu}} & {{\boldsymbol{K}}}_{{LL}} & {{\boldsymbol{K}}}_{L\theta } & {\bf{0}} & {\bf{0}}\\ {\bf{0}} & {{\boldsymbol{K}}}_{\theta L} & {{\boldsymbol{K}}}_{\theta \theta } & {{\boldsymbol{K}}}_{\theta C} & {{\boldsymbol{K}}}_{\theta \varphi }\\ {\bf{0}} & {\bf{0}} & {{\boldsymbol{K}}}_{C\theta } & {{\boldsymbol{K}}}_{{CC}} & {{\boldsymbol{K}}}_{C\varphi }\\ {\bf{0}} & {\bf{0}} & {\bf{0}} & {{\boldsymbol{K}}}_{\varphi C} & {\bf{0}}\end{array}\right]}_{i}{\left[\begin{array}{c}\delta {\bf{u}}\\ \delta {{\bf{C}}}_{L}\\ \delta {\boldsymbol{\theta }}\\ \delta {{\bf{C}}}_{\pi }\\ \delta {\boldsymbol{\varphi }}\end{array}\right]}_{i+1}=-{\left[\begin{array}{c}{{\bf{f}}}_{u}\\ {{\bf{f}}}_{L}\\ {{\bf{f}}}_{\theta }\\ {{\bf{f}}}_{C\pi }\\ {{\bf{f}}}_{\varphi }\end{array}\right]}_{i}\end{eqnarray} \tag{ A1 }$$

is iteratively solved, such that the iterative increment in the solution (e.g., δ u_i+1) is added to the state vector at the end of each iteration. Here, the forces f and tangential matrices K are updated at the start of each iteration i using the then-available solution state (e.g., u_i ). The expressions for the stiffness matrices are provided below, categorised as a function of the relevant domain (metal, electrolyte, metal-electrolyte interface).

Appendix. Metal sub-domain

For the momentum balance, the definition of the force vector f_u is given in Eq. 38, and its associated tangential matrix is:

$$\begin{eqnarray}&&{{\boldsymbol{K}}}_{{uu}}={\int }_{{{\rm{\Omega }}}_{m}}{{\boldsymbol{B}}}_{u}^{T}{\boldsymbol{D}}{{\boldsymbol{B}}}_{u}\,{\rm{d}}{{\rm{\Omega }}}_{m}\end{eqnarray} \tag{ A2 }$$

Similarly, the force vector for the lattice hydrogen mass balance is given in Eq. 42, and the relevant tangential matrices are given by

$$\begin{eqnarray}&&{{\boldsymbol{K}}}_{{Lu}}=-{\int }_{{{\rm{\Omega }}}_{m}}\displaystyle \frac{E}{3(1-2\nu )}\displaystyle \frac{{D}_{L}{\overline{V}}_{H}}{{RT}}{\left({\boldsymbol{\nabla }}{{\bf{N}}}_{L}\right)}^{T}\left({{\bf{N}}}_{L}{{\bf{C}}}_{L}^{t+{\rm{\Delta }}t}\right){{\boldsymbol{B}}}_{u}^{* }\,{\rm{d}}{{\rm{\Omega }}}_{m}\end{eqnarray} \tag{ A3 }$$

$$\begin{eqnarray}\begin{array}{ccc}{{\boldsymbol{K}}}_{{LL}} & = & -{\int }_{{{\rm{\Omega }}}_{m}}\frac{E}{3(1-2\nu )}\frac{{D}_{L}{\bar{V}}_{H}}{{RT}}{\left({\boldsymbol{\nabla }}{{\bf{N}}}_{L}\right)}^{T}\left({{\boldsymbol{B}}}_{u}^{* }{{\bf{u}}}^{t+{\rm{\Delta }}t}\right){{\bf{N}}}_{L}\,{\rm{d}}{{\rm{\Omega }}}_{m}\\ & & +{\int }_{{{\rm{\Omega }}}_{m}}\frac{1}{{\rm{\Delta }}t}{{\bf{N}}}_{L}^{T}{{\bf{N}}}_{L}\,{\rm{d}}{{\rm{\Omega }}}_{m}\\ & & +{\int }_{{{\rm{\Omega }}}_{m}}\frac{{D}_{L}}{1-{{\bf{N}}}_{L}{{\boldsymbol{C}}}_{L}^{t+{\rm{\Delta }}t}/{N}_{L}}{\left({\boldsymbol{\nabla }}{{\bf{N}}}_{L}\right)}^{T}{\boldsymbol{\nabla }}{{\bf{N}}}_{L}\\ & & +\frac{{D}_{L}{\boldsymbol{\nabla }}{{\bf{N}}}_{L}{{\bf{C}}}_{L\,i}^{t+{\rm{\Delta }}t}}{{N}_{L}{\left(1-{{\bf{N}}}_{L}{{\boldsymbol{C}}}_{L}^{t+{\rm{\Delta }}t}/{N}_{L}\right)}^{2}}{\left({\boldsymbol{\nabla }}{{\bf{N}}}_{L}\right)}^{T}{{\bf{N}}}_{L}\,{\rm{d}}{{\rm{\Omega }}}_{m}\\ & & -\displaystyle \sum _{{el}}\displaystyle \sum _{{nd}}{\boldsymbol{W}}({nd})\left(-{k}_{A}{{\boldsymbol{\theta }}}_{{ads}}^{t+{\rm{\Delta }}t}({nd})\right.\\ & & \left.-{k}_{A}^{^{\prime} }(1-{{\boldsymbol{\theta }}}_{{ads}}^{t+{\rm{\Delta }}t}({nd}))\right){\boldsymbol{I}}({nd},{C}_{L})\end{array}\end{eqnarray} \tag{ A4 }$$

$$\begin{eqnarray}&&\begin{array}{l}{{\boldsymbol{K}}}_{L\theta }=-\displaystyle \sum _{{el}}\displaystyle \sum _{{nd}}{\boldsymbol{W}}({nd})\left({k}_{A}({N}_{L}-{{\bf{C}}}_{L}^{t+{\rm{\Delta }}t}({nd}))\right.\\ \left.+{k}_{A}^{{\prime} }{{\bf{C}}}_{L}^{t+{\rm{\Delta }}t}({nd})\right){\boldsymbol{I}}({nd},\theta )\end{array}\end{eqnarray} \tag{ A5 }$$

where the matrix I (nd, C_L ) is used to assign the lumped integration terms to the location within the tangential sub-matrices associated with the node nd and degree of freedom C_L .

Appendix. Electrolyte sub-domain

The force vectors for the mass balance and electroneutrality condition within the electrolyte are given by Eqs. 42 and 43. The related tangential matrix terms are given as follows. Firstly, for the mass balance one obtains:

$$\begin{eqnarray}&&{{\boldsymbol{K}}}_{{C}_{{{\rm{H}}}^{+}}\theta }=\displaystyle \sum _{{iel}}\displaystyle \sum _{{nd}}{\boldsymbol{W}}({nd})\left(\displaystyle \frac{\partial {\nu }_{{Va}}}{\partial \theta }-\displaystyle \frac{\partial {\nu }_{{Va}}^{{\prime} }}{\partial \theta }+\displaystyle \frac{\partial {\nu }_{{Ha}}}{\partial \theta }\right){\boldsymbol{I}}({nd},\theta )\end{eqnarray} \tag{ A6 }$$

$$\begin{eqnarray}&&{{\boldsymbol{K}}}_{{C}_{{\mathrm{OH}}^{-}}\theta }=-\displaystyle \sum _{{iel}}\displaystyle \sum _{{nd}}{\boldsymbol{W}}({nd})\left(\displaystyle \frac{\partial {\nu }_{{Vb}}}{\partial \theta }-\displaystyle \frac{\partial {\nu }_{{Vb}}^{{\prime} }}{\partial \theta }+\displaystyle \frac{\partial {\nu }_{{Hb}}}{\partial \theta }\right){\boldsymbol{I}}({nd},\theta )\end{eqnarray} \tag{ A7 }$$

$$\begin{eqnarray}\begin{array}{ccc}{{\boldsymbol{K}}}_{{C}_{{{\rm{H}}}^{+}}{C}_{{{\rm{H}}}^{+}}} & = & {{\boldsymbol{H}}}_{{C}_{{{\rm{H}}}^{+}}{C}_{{{\rm{H}}}^{+}}}+\displaystyle \sum _{{iel}}\displaystyle \sum _{{nd}}{\boldsymbol{W}}({nd})\\ & & \cdot \left(\frac{\partial {\nu }_{{Va}}}{\partial {C}_{{{\rm{H}}}^{+}}}+\frac{\partial {\nu }_{{Ha}}}{\partial {C}_{{{\rm{H}}}^{+}}}\right){\boldsymbol{I}}({nd},{C}_{{{\rm{H}}}^{+}})\\ & & +\displaystyle \sum _{{el}}\displaystyle \sum _{{nd}}{\bf{W}}({nd}){k}_{{eq}}{{\bf{C}}}_{{\mathrm{OH}}^{-}}({nd}){\boldsymbol{I}}({nd},{C}_{{{\rm{H}}}^{+}})\\ & & +\displaystyle \sum _{{el}}\displaystyle \sum _{{nd}}{\bf{W}}({nd}){k}_{{fe}}^{{\prime} }{{\bf{C}}}_{{\mathrm{FeOH}}^{+}}({nd}){\boldsymbol{I}}({nd},{C}_{{{\rm{H}}}^{+}})\end{array}\end{eqnarray} \tag{ A8 }$$

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{K}}}_{{C}_{{\mathrm{OH}}^{-}}{C}_{{\mathrm{OH}}^{-}}} & = & {{\boldsymbol{H}}}_{{C}_{{\mathrm{OH}}^{-}}{C}_{{\mathrm{OH}}^{-}}}-\displaystyle \sum _{{iel}}\displaystyle \sum _{{nd}}{\boldsymbol{W}}({nd})\displaystyle \frac{\partial {\nu }_{{Vb}}^{{\prime} }}{\partial {C}_{{\mathrm{OH}}^{-}}}{\boldsymbol{I}}({nd},{C}_{{\mathrm{OH}}^{-}})\\ & & +\displaystyle \sum _{{el}}\displaystyle \sum _{{nd}}{\bf{W}}({nd}){k}_{{eq}}{{\bf{C}}}_{{{\rm{H}}}^{+}}({nd}){\boldsymbol{I}}({nd},{C}_{{\mathrm{OH}}^{-}})\end{array}\end{eqnarray} \tag{ A9 }$$

$$\begin{eqnarray}&&{{\boldsymbol{K}}}_{{C}_{{\mathrm{Na}}^{+}}{C}_{{\mathrm{Na}}^{+}}}={{\boldsymbol{K}}}_{{C}_{{\mathrm{Cl}}^{-}}{C}_{{\mathrm{Cl}}^{-}}}={{\boldsymbol{H}}}_{{C}_{\pi }{C}_{\pi }}\end{eqnarray} \tag{ A10 }$$

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{K}}}_{{C}_{{\mathrm{Fe}}^{2+}}{C}_{{\mathrm{Fe}}^{2+}}} & = & {{\boldsymbol{H}}}_{{C}_{{\mathrm{Fe}}^{2+}}{C}_{{\mathrm{Fe}}^{2+}}}\\ & + & \displaystyle \sum _{{iel}}\displaystyle \sum _{{nd}}{\boldsymbol{W}}({nd})\left({k}_{c}{e}^{-{\alpha }_{c}\left({E}_{m}-{\boldsymbol{\varphi }}({nd})-{E}_{{eq},c}\right)}\right){\boldsymbol{I}}({nd},{C}_{{\mathrm{Fe}}^{2+}})\\ & + & \displaystyle \sum _{{el}}\displaystyle \sum _{{nd}}{\bf{W}}({nd}){k}_{{fe}}{\boldsymbol{I}}({nd},{C}_{{\mathrm{Fe}}^{2+}})\end{array}\end{eqnarray} \tag{ A11 }$$

$$\begin{eqnarray}&&\begin{array}{l}{{\boldsymbol{K}}}_{{C}_{{\mathrm{FeOH}}^{+}}{C}_{{\mathrm{FeOH}}^{+}}}={{\boldsymbol{H}}}_{{C}_{{\mathrm{FeOH}}^{+}}{C}_{{\mathrm{FeOH}}^{+}}}\\ +\displaystyle \sum _{{el}}\displaystyle \sum _{{nd}}{\bf{W}}({nd})\left({k}_{{fe}}^{{\prime} }{{\bf{C}}}_{{{\rm{H}}}^{+}}({nd})+{k}_{{feoh}}\right){\boldsymbol{I}}({nd},{C}_{{\mathrm{FeOH}}^{+}})\end{array}\end{eqnarray} \tag{ A12 }$$

$$\begin{eqnarray}&&{{\boldsymbol{K}}}_{{C}_{{{\rm{H}}}^{+}}{C}_{{\mathrm{OH}}^{-}}}=\displaystyle \sum _{{el}}\displaystyle \sum _{{nd}}{\bf{W}}({nd}){k}_{{eq}}{{\bf{C}}}_{{{\rm{H}}}^{+}}({nd}){\boldsymbol{I}}({nd},{C}_{{\mathrm{OH}}^{-}})\end{eqnarray} \tag{ A13 }$$

$$\begin{eqnarray}&&{{\boldsymbol{K}}}_{{C}_{{\mathrm{OH}}^{-}}{C}_{{{\rm{H}}}^{+}}}=\displaystyle \sum _{{el}}\displaystyle \sum _{{nd}}{\bf{W}}({nd}){k}_{{eq}}{{\bf{C}}}_{{\mathrm{OH}}^{-}}({nd}){\boldsymbol{I}}({nd},{C}_{{{\rm{H}}}^{+}})\end{eqnarray} \tag{ A14 }$$

$$\begin{eqnarray}&&{{\boldsymbol{K}}}_{{C}_{{{\rm{H}}}^{+}}{C}_{{\mathrm{Fe}}^{2+}}}=-\displaystyle \sum _{{el}}\displaystyle \sum _{{nd}}{\bf{W}}({nd}){k}_{{fe}}{\boldsymbol{I}}({nd},{C}_{{\mathrm{Fe}}^{2+}})\end{eqnarray} \tag{ A15 }$$

$$\begin{eqnarray}&&{{\boldsymbol{K}}}_{{C}_{{{\rm{H}}}^{+}}{C}_{{\mathrm{FeOH}}^{+}}}=\displaystyle \sum _{{el}}\displaystyle \sum _{{nd}}{\bf{W}}({nd})\left({k}_{{fe}}^{{\prime} }{{\bf{C}}}_{{{\rm{H}}}^{+}}({nd})-{k}_{{feoh}}\right){\boldsymbol{I}}({nd},{C}_{{\mathrm{FeOH}}^{+}})\end{eqnarray} \tag{ A16 }$$

$$\begin{eqnarray}&&{{\boldsymbol{K}}}_{{C}_{{\mathrm{Fe}}^{2+}}{C}_{{{\rm{H}}}^{+}}}=-\displaystyle \sum _{{el}}\displaystyle \sum _{{nd}}{\bf{W}}({nd}){k}_{{fe}}^{{\prime} }{{\bf{C}}}_{{\mathrm{FeOH}}^{+}}({nd}){\boldsymbol{I}}({nd},{C}_{{{\rm{H}}}^{+}})\end{eqnarray} \tag{ A17 }$$

$$\begin{eqnarray}&&{{\boldsymbol{K}}}_{{C}_{{\mathrm{Fe}}^{2+}}{C}_{{\mathrm{FeOH}}^{+}}}=-\displaystyle \sum _{{el}}\displaystyle \sum _{{nd}}{\bf{W}}({nd}){k}_{{fe}}^{{\prime} }{{\bf{C}}}_{{{\rm{H}}}^{+}}({nd}){\boldsymbol{I}}({nd},{C}_{{\mathrm{FeOH}}^{+}})\end{eqnarray} \tag{ A18 }$$

$$\begin{eqnarray}&&{{\boldsymbol{K}}}_{{C}_{{\mathrm{FeOH}}^{+}}{C}_{{{\rm{H}}}^{+}}}=\displaystyle \sum _{{el}}\displaystyle \sum _{{nd}}{\bf{W}}({nd}){k}_{{fe}}^{{\prime} }{{\bf{C}}}_{{\mathrm{FeOH}}^{+}}({nd}){\boldsymbol{I}}({nd},{C}_{{{\rm{H}}}^{+}})\end{eqnarray} \tag{ A19 }$$

$$\begin{eqnarray}&&{{\boldsymbol{K}}}_{{C}_{{\mathrm{FeOH}}^{+}}{C}_{{\mathrm{Fe}}^{2+}}}=-\displaystyle \sum _{{el}}\displaystyle \sum _{{nd}}{\bf{W}}({nd}){k}_{{fe}}{\boldsymbol{I}}({nd},{C}_{{\mathrm{Fe}}^{2+}})\end{eqnarray} \tag{ A20 }$$

$$\begin{eqnarray}&&\begin{array}{l}{{\boldsymbol{K}}}_{{C}_{{{\rm{H}}}^{+}}\varphi }={{\boldsymbol{H}}}_{{C}_{{{\rm{H}}}^{+}}\varphi }+\displaystyle \sum _{{iel}}\displaystyle \sum _{{nd}}{\boldsymbol{W}}({nd})\\ \times \left(\displaystyle \frac{\partial {\nu }_{{Va}}}{\partial \varphi }-\displaystyle \frac{\partial {\nu }_{{Va}}^{{\prime} }}{\partial \varphi }+\displaystyle \frac{\partial {\nu }_{{Ha}}}{\partial \varphi }\right){\boldsymbol{I}}({nd},\varphi )\end{array}\end{eqnarray} \tag{ A21 }$$

$$\begin{eqnarray}&&\begin{array}{l}{{\boldsymbol{K}}}_{{C}_{{\mathrm{OH}}^{-}}\varphi }={{\boldsymbol{H}}}_{{C}_{{\mathrm{OH}}^{-}}\varphi }-\displaystyle \sum _{{iel}}\displaystyle \sum _{{nd}}{\boldsymbol{W}}({nd})\\ \times \,\left(\displaystyle \frac{\partial {\nu }_{{Vb}}}{\partial \varphi }-\displaystyle \frac{\partial {\nu }_{{Vb}}^{{\prime} }}{\partial \varphi }+\displaystyle \frac{\partial {\nu }_{{Hb}}}{\partial \varphi }\right){\boldsymbol{I}}({nd},\varphi )\end{array}\end{eqnarray} \tag{ A22 }$$

$$\begin{eqnarray}&&{{\boldsymbol{K}}}_{{C}_{{\mathrm{Fe}}^{2+}}\varphi }={{\boldsymbol{H}}}_{{C}_{{\mathrm{Fe}}^{2+}}\varphi }+\displaystyle \sum _{{iel}}\displaystyle \sum _{{nd}}{\boldsymbol{W}}({nd})\left(\displaystyle \frac{\partial {\nu }_{c}}{\partial \varphi }-\displaystyle \frac{\partial {\nu }_{c}^{\prime} }{\partial \varphi }\right){\boldsymbol{I}}({nd},\varphi )\end{eqnarray} \tag{ A23 }$$

$$\begin{eqnarray}&&{{\boldsymbol{K}}}_{{C}_{{\mathrm{Na}}^{+}}\varphi }={{\boldsymbol{K}}}_{{C}_{{\mathrm{Cl}}^{-}}\varphi }={{\boldsymbol{K}}}_{{C}_{{\mathrm{FeOH}}^{+}}\varphi }={{\boldsymbol{H}}}_{{C}_{\pi }\varphi }\end{eqnarray} \tag{ A24 }$$

where the matrix ${{\boldsymbol{H}}}_{{C}_{\pi }{C}_{\pi }}$ is defined as:

$$\begin{eqnarray}\begin{array}{ccc}{{\boldsymbol{H}}}_{{C}_{\pi }{C}_{\pi }} & = & {\int }_{{{\rm{\Omega }}}_{e}}\frac{1}{{\rm{\Delta }}t}{{\bf{N}}}_{C}^{T}{{\bf{N}}}_{C}\,{\rm{d}}{{\rm{\Omega }}}_{e}+{\int }_{{{\rm{\Omega }}}_{e}}{D}_{\pi }{\left({\boldsymbol{\nabla }}{{\bf{N}}}_{c}\right)}^{T}{\boldsymbol{\nabla }}{{\bf{N}}}_{C}\,{\rm{d}}{{\rm{\Omega }}}_{e}\\ & & +{\int }_{{{\rm{\Omega }}}_{e}}\frac{{D}_{\pi }{z}_{\pi }F}{{RT}}{\left({\boldsymbol{\nabla }}{{\bf{N}}}_{c}\right)}^{T}\left({\boldsymbol{\nabla }}{{\bf{N}}}_{\varphi }{{\boldsymbol{\varphi }}}^{t+{\rm{\Delta }}t}\right){{\bf{N}}}_{C}\,{\rm{d}}{{\rm{\Omega }}}_{e}\end{array}\end{eqnarray} \tag{ A25 }$$

and the matrix ${{\boldsymbol{H}}}_{{C}_{\pi }\varphi }$ reads:

$$\begin{eqnarray}&&{{\boldsymbol{H}}}_{{C}_{\pi }\varphi }={\int }_{{{\rm{\Omega }}}_{e}}\displaystyle \frac{{D}_{\pi }{z}_{\pi }F}{{RT}}{\left({\boldsymbol{\nabla }}{{\bf{N}}}_{c}\right)}^{T}\left({{\bf{N}}}_{c}{{\bf{C}}}_{\pi }^{t+{\rm{\Delta }}t}\right){\boldsymbol{\nabla }}{{\bf{N}}}_{\varphi }\,{\rm{d}}{{\rm{\Omega }}}_{e}\end{eqnarray} \tag{ A26 }$$

Secondly, for the electroneutrality condition, the tangential matrix terms are given by:

$$\begin{eqnarray}&&{{\boldsymbol{K}}}_{\varphi {C}_{\pi }}={\int }_{{{\rm{\Omega }}}_{e}}{z}_{\pi }{{\bf{N}}}_{\varphi }^{T}{{\bf{N}}}_{c}\,{\rm{d}}{{\rm{\Omega }}}_{e}\end{eqnarray} \tag{ A27 }$$

Appendix. Metal-electrolyte interface

Finally, the relevant tangential matrices for the metal-electrolyte interface are provided. The force vector related to the mass balance of the adsorbed hydrogen is given in Eq. 46. Considering the case in which a lumped integration scheme is used for all surface reactions, the tangential matrix terms are given by:

$$\begin{eqnarray}&&\begin{array}{c}{{\boldsymbol{K}}}_{\theta L}=\displaystyle \sum _{{iel}}\displaystyle \sum _{{nd}}{\boldsymbol{W}}({nd})\left(-{k}_{A}{{\boldsymbol{\theta }}}_{{ads}}^{t+{\rm{\Delta }}t}({nd})\right.\\ \left.-\,{k}_{A}^{^{\prime} }(1-{{\boldsymbol{\theta }}}_{{ads}}^{t+{\rm{\Delta }}t}({nd}))\right){\boldsymbol{I}}({nd},{C}_{L})\end{array}\end{eqnarray} \tag{ A28 }$$

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{K}}}_{\theta \theta } & = & {\displaystyle \int }_{{{\rm{\Gamma }}}_{{int}}}\displaystyle \frac{{N}_{{ads}}}{{\rm{\Delta }}t}{{\bf{N}}}_{\theta }^{T}{{\bf{N}}}_{\theta }\,{\rm{d}}{{\rm{\Gamma }}}_{{int}}\\ & - & \displaystyle \sum _{{iel}}\displaystyle \sum _{{nd}}{\boldsymbol{W}}({nd})\left(\displaystyle \frac{\partial {\nu }_{{Va}}}{\partial \theta }-\displaystyle \frac{\partial {\nu }_{{Va}}^{{\prime} }}{\partial \theta }-\displaystyle \frac{\partial {\nu }_{{Ha}}}{\partial \theta }-2\displaystyle \frac{\partial {\nu }_{T}}{\partial \theta }-\displaystyle \frac{\partial {\nu }_{A}}{\partial \theta }\right.\\ & & \left.+\displaystyle \frac{\partial {\nu }_{A}^{{\prime} }}{\partial \theta }+\displaystyle \frac{\partial {\nu }_{{Vb}}}{\partial \theta }-\displaystyle \frac{\partial {\nu }_{{Vb}}^{{\prime} }}{\partial \theta }-\displaystyle \frac{\partial {\nu }_{{Hb}}}{\partial \theta }\right){\boldsymbol{I}}({nd},\theta )\end{array}\end{eqnarray} \tag{ A29 }$$

$$\begin{eqnarray}\begin{array}{ccc}{{\boldsymbol{K}}}_{\theta {C}_{{{\rm{H}}}^{+}}} & = & -\displaystyle \sum _{{iel}}\displaystyle \sum _{{nd}}{\boldsymbol{W}}({nd})\left({k}_{{Va}}{e}^{-{\alpha }_{{Va}}\left({E}_{m}-{{\boldsymbol{\varphi }}}^{t+{\rm{\Delta }}t}({nd})-{E}_{{eq},{Va}}\right)\tfrac{F}{{RT}}}\right.\\ & & \left.-{k}_{{Ha}}{e}^{-{\alpha }_{{Ha}}\left({E}_{m}-{{\boldsymbol{\varphi }}}^{t+{\rm{\Delta }}t}({nd})-{E}_{{eq},{Ha}}\right)\tfrac{F}{{RT}}}\right){\boldsymbol{I}}({nd},{C}_{{{\rm{H}}}^{+}})\end{array}\end{eqnarray} \tag{ A30 }$$

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{K}}}_{\theta {C}_{{\mathrm{OH}}^{-}}} & = & \displaystyle \sum _{{iel}}\displaystyle \sum _{{nd}}{\boldsymbol{W}}({nd}){k}_{{Vb}}^{{\prime} }{{\boldsymbol{\theta }}}^{t+{\rm{\Delta }}t}({nd})\\ & & \cdot {e}^{(1-{\alpha }_{{Vb}})\left({E}_{m}-{{\boldsymbol{\varphi }}}^{t+{\rm{\Delta }}t}({nd})-{E}_{{eq},{Vb}}\right)\tfrac{F}{{RT}}}{\boldsymbol{I}}({nd},{C}_{{\mathrm{OH}}^{-}})\end{array}\end{eqnarray} \tag{ A31 }$$

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{K}}}_{\theta \varphi } & = & -\displaystyle \sum _{{iel}}\displaystyle \sum _{{nd}}{\boldsymbol{W}}({nd})\left(\displaystyle \frac{\partial {\nu }_{{Va}}}{\partial \varphi }-\displaystyle \frac{\partial {\nu }_{{Va}}^{{\prime} }}{\partial \varphi }-\displaystyle \frac{\partial {\nu }_{{Ha}}}{\partial \varphi }\right.\\ & & \left.+\displaystyle \frac{\partial {\nu }_{{Vb}}}{\partial \varphi }-\displaystyle \frac{\partial {\nu }_{{Vb}}^{{\prime} }}{\partial \varphi }-\displaystyle \frac{\partial {\nu }_{{Hb}}}{\partial \varphi }\right){\boldsymbol{I}}({nd},\varphi )\end{array}\end{eqnarray} \tag{ A32 }$$

Axial symmetry is exploited in the tensile rod cases to simulate the three-dimensional domain depicted in Fig. 12. Thus, a domain defined in the coordinate system (r, θ, z) is instead evaluated within a two-dimensional coordinate system x = (r, z), assuming the results are constant in the θ direction. The main notable difference resulting from this transformation is the change to the integration scheme. For instance, the absorption reaction, Eq. 49, is integrated as follows when considering axisymmetry and lumped integration:

$$\begin{eqnarray}&&\begin{array}{c}{\int }_{{{\rm{\Gamma }}}_{{int}}}{{\bf{N}}}_{\theta }^{T}({\nu }_{{Va}}-{\nu }_{{Va}}^{{\prime} }){\rm{d}}{{\rm{\Gamma }}}_{{int}}(r,\theta ,z)\\ \quad ={\int }_{{{\rm{\Gamma }}}_{{int}}}2\pi r{{\bf{N}}}_{\theta }^{T}({\nu }_{{Va}}-{\nu }_{{Va}}^{{\prime} }){\rm{d}}{{\rm{\Gamma }}}_{{int}}(r,z)\\ \quad =\displaystyle \sum _{{iel}}\displaystyle \sum _{{nd}}{\boldsymbol{W}}({nd})\left({k}_{A}({N}_{L}-{{\bf{C}}}_{L}({nd})){{\boldsymbol{\theta }}}_{{ads}}({nd})\right.\\ \quad \left.-{k}_{A}^{^{\prime} }{{\bf{C}}}_{L}({nd})(1-{{\boldsymbol{\theta }}}_{{ads}}({nd}))\right)\end{array}\end{eqnarray} \tag{ B1 }$$

with the lumped integration weight vector W including the effects of the axisymmetry as:

$$\begin{eqnarray}&&{\bf{W}}={\int }_{{{\rm{\Gamma }}}_{{int}}}2\pi r{{\bf{N}}}^{T}\,{\rm{d}}{{\rm{\Gamma }}}_{{int}}(r,z)=\displaystyle \sum _{{ip}}2\pi {w}_{{ip}}{r}_{{ip}}{{\bf{N}}}^{T}({{\bf{x}}}_{{ip}})\end{eqnarray} \tag{ B2 }$$

And we proceed similarly with the volume reactions. For example, the water auto-ionisation reaction is described by

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{R}}}_{{{OH}}^{-}} & = & -{\displaystyle \int }_{{{\rm{\Omega }}}_{e}}{{\bf{N}}}_{c}^{T}\left({K}_{w}-\left({{\bf{N}}}_{c}{{\bf{C}}}_{{{\rm{H}}}^{+}}^{t+{\rm{\Delta }}t}\right)\left({{\bf{N}}}_{c}{{\bf{C}}}_{{\mathrm{OH}}^{-}}^{t+{\rm{\Delta }}t}\right)\right){\rm{d}}{{\rm{\Omega }}}_{e}(r,\theta ,z)\\ & = & -{\displaystyle \int }_{{{\rm{\Omega }}}_{e}}2\pi r{{\bf{N}}}_{c}^{T}\left({K}_{w}-\left({{\bf{N}}}_{c}{{\bf{C}}}_{{{\rm{H}}}^{+}}^{t+{\rm{\Delta }}t}\right)\left({{\bf{N}}}_{c}{{\bf{C}}}_{{\mathrm{OH}}^{-}}^{t+{\rm{\Delta }}t}\right)\right){\rm{d}}{{\rm{\Omega }}}_{e}(r,z)\\ & = & -\displaystyle \sum _{{el}}\displaystyle \sum _{{nd}}{\bf{W}}({nd})\left({K}_{W}-{{\bf{C}}}_{{{\rm{H}}}^{+}}({nd}){{\bf{C}}}_{{\mathrm{OH}}^{-}}({nd})\right)\end{array}\end{eqnarray} \tag{ B3 }$$

using lumped weights:

$$\begin{eqnarray}\begin{array}{rcl}{\bf{W}} & = & {\displaystyle \int }_{{{\rm{\Omega }}}_{{el}}}{{\bf{N}}}^{T}\,{\rm{d}}{{\rm{\Omega }}}_{{el}}(r,\theta ,z)={\displaystyle \int }_{{{\rm{\Omega }}}_{{el}}}2\pi r{{\bf{N}}}^{T}\,{\rm{d}}{{\rm{\Omega }}}_{{el}}(r,z)\\ & = & \displaystyle \sum _{{ip}}2\pi {r}_{{ip}}{x}_{{ip}}{{\bf{N}}}^{T}({{\bf{x}}}_{{ip}})\end{array}\end{eqnarray} \tag{ B4 }$$

Another aspect to consider is the change in the strain-displacement matrix, which now reads

$$\begin{eqnarray}{{\boldsymbol{B}}}_{u}=\left[\begin{array}{cccccc}\tfrac{\partial {N}_{u1}}{\partial r} & \tfrac{\partial {N}_{u2}}{\partial r} & \cdots & 0 & 0 & \cdots \\ 0 & 0 & \cdots & \tfrac{\partial {N}_{u1}}{\partial z} & \tfrac{\partial {N}_{u2}}{\partial z} & \cdots \\ {N}_{u1}/r & {N}_{u2}/r & \cdots & 0 & 0 & \cdots \\ \tfrac{\partial {N}_{u1}}{\partial z} & \tfrac{\partial {N}_{u2}}{\partial z} & \cdots & \tfrac{\partial {N}_{u1}}{\partial r} & \tfrac{\partial {N}_{u2}}{\partial r} & \cdots \end{array}\right]\end{eqnarray} \tag{ B5 }$$

Accordingly, the interpolation matrix used to obtain the hydrostatic stress gradient form the displacement field, ${\boldsymbol{\nabla }}{\sigma }_{H}=E/(3(1-2\nu )){{\boldsymbol{B}}}_{u}^{* }{\bf{u}}$, is formulated as

$$\begin{eqnarray}{{\boldsymbol{B}}}_{u}^{* }=\left[\begin{array}{cccccc}\tfrac{{\partial }^{2}{N}_{u1}}{\partial {r}^{2}}+\tfrac{1}{r}\tfrac{\partial {N}_{u1}}{\partial r}-\tfrac{{N}_{u1}}{{r}^{2}} & \tfrac{{\partial }^{2}{N}_{u2}}{\partial {r}^{2}}+\tfrac{1}{r}\tfrac{\partial {N}_{u2}}{\partial r}-\tfrac{{N}_{u2}}{{r}^{2}} & \cdots & \tfrac{{\partial }^{2}{N}_{u1}}{\partial r\partial z} & \tfrac{{\partial }^{2}{N}_{u2}}{\partial r\partial z} & \cdots \\ \tfrac{{\partial }^{2}{N}_{u1}}{\partial r\partial z} & \tfrac{{\partial }^{2}{N}_{u2}}{\partial r\partial z} & \cdots & \tfrac{{\partial }^{2}{N}_{u1}}{\partial {z}^{2}} & \displaystyle \frac{{\partial }^{2}{N}_{u2}}{\partial {z}^{2}} & \cdots \end{array}\right]\end{eqnarray} \tag{ B6 }$$