Turing patterns in a 3D morpho-chemical bulk-surface reaction-diffusion system for battery modeling

1.   Introduction

The formation of spatio-temporal structures in electrodeposition is a relevant physical phenomenon, as it impacts several applications, ranging from the durability and efficiency of batteries to electroplating ^[1]. The onset of spatio-temporal structures on the cathodic surface was proven to be initiated by a Turing morphogenetic mechanism, where the physics are modeled by a suitable reaction-diffusion system (RDS), called DIB model after the authors ^[1,2], whose spatial domain is the electrodic surface. In the DIB model, the spatial domain is assumed to be fixed and does not change over time, as the growth/corrosion effects are fully modeled by the dynamics of the system. By tweaking the parameters of the DIB model, it is possible to successfully simulate spatial ^[1] or spatio-temporal patterns ^[2] of various morphological classes that are experimentally observed under appropriate physical and chemical conditions; these include spatial patterns such as spots, holes, stripes, labyrinths, and spiral waves. The effectiveness of the DIB model has justified the development of extensions and ameliorations, such as the introduction of cross-diffusion ^[3] and the generalization of the spatial domain to be a curved surface ^[4,5].

As it stands, one of the limitations of the DIB model is that it does not fully accounts for the effects of non-uniform electrolyte concentration in a neighborhood of the electrode. Experimentally, such non uniform concentration is induced by the spatial structures arising on the electrode and, in turn, affects further structure development. In this regard, the electrode-electrolyte system has a two-way coupling that, in the long run, can drastically affect the resulting morphological class. In this paper we propose the bulk-surface DIB (BS-DIB) model in three space dimensions to fill this gap. In the proposed model, the surface represents the electrode (where the electrodeposition takes place), while the 3D bulk models the electrolyte. The physical two-way coupling mentioned above causes the proposed model to take the form of a coupled bulk-surface reaction-diffusion system (BS-RDS) ^[6,7,8].

For domains of general shape, different numerical methods were developed for the spatial approximation of BS-RDSs, such as the Bulk-Surface Finite Element Method (BS-FEM) ^[7], the Cut Finite Element Method ^[9], unfitted finite element methods ^[10], and meshless kernel methods ^[11], just to mention a few. In all of these methods, the spatially discrete problem takes the form of a large ODE system, whose dimension is equal to the number of spatial degrees of freedom. Thus, the high level of spatial resolution required by RDSs and BS-RDSs, together with the curse of dimensionality (3D), makes the numerical approximation of the BS-DIB model a challenging computational task. In the present context, where the bulk domain is a cube, a bespoke tensorized technique called Matrix-Oriented Finite Element Method (MO-FEM) ^[12,13] can be exploited to take advantage of the special geometry and drastically reduce the computational execution times. However, it is worth noting that the BS-DIB model can exhibit spatial patterns only in a neighborhood of the surface, hence a uniform spatial discretization is computationally inefficient. For this reason, we exploit the geometric flexibility of the bulk-surface virtual element method (BS-VEM) ^[8] to adopt a graded cubic mesh that is highly refined close to the surface and much coarser away from the surface. Such a mesh is simultaneously graded and entirely composed of cubic-shaped elements. Such a combination entails the presence of hanging nodes and edges, that are naturally handled by the BS-VEM and are not admissible in the BS-FEM. Compared to the MO-FEM, we show that the BS-VEM on such mesh exhibits a much lower number of degrees of freedom on equal level of spatial refinement in a neighborhood of the surface, where high spatial accuracy is actually required and produces patterns of the same morphological class. It needs to be noted that Turing patterns are highly sensitive to initial conditions, which are bound to be different between MO-FEM ad BS-VEM since the spatial meshes are different, hence obtaining the same morphological class with both methods is a sensible benchmark.

We present a wide range of numerical simulations for both the 2D DIB and the 3D BS-DIB models on equal parameters, to showcase the effect of the bulk-surface coupling. From the experiments we draw the following conclusions. First, the BS-DIB model appears to have a larger Turing region in the parameter space compared to the DIB model. In fact, for several choices of the parameters outside the Turing region of the DIB model, only the BS-DIB model exhibits spatial patterns. Second, when the DIB model exhibits spatial patterns, the BS-DIB model still exhibits patterns, but of different morphological class, thereby further highlighting the impact of the bulk-surface coupling. A rigorous analysis of the Turing instability for the BS-DIB model will be addressed in future work.

The structure of the paper is as follows. In Section 2 we introduce the BS-DIB model, we give its physical interpretation and we analyse the stability of a relevant spatially uniform equilibrium in the absence of diffusion. In Section 3, we recall from ^[14] the BS-VEM and we present a bespoke graded polyhedral mesh that allows the BS-VEM to outperform the MO-FEM when solving the BS-DIB model. In Section 4, we present an extensive list of numerical experiments that empirically demonstrate the effect of the bulk-surface coupling on pattern formation. Conclusions are drawn in Section 5.

2.   The bulk-surface DIB model on the cube

In this section, we present for the first time the morpho-chemical bulk-surface DIB (BS-DIB) PDE system for battery modeling. For simplicity, we consider a cube $\Omega$ for the 3D bulk that represents the electrolyte, and its bottom face $\Gamma$ as the surface representing the electrode where the electrodeposition process takes place, see Figure 1.

Hence, let $\Omega = [0, L]^3$ be a cube of edge $L > 0$, let $\Gamma : = [0, L]^2\times \{0\}$ be the bottom face of $\Omega$, let $\Gamma_T = [0, L]^2\times \{L\}$ be the top face of $\Omega$ and let $\Gamma_L = \partial \Omega \setminus (\Gamma\cup\Gamma_T)$ be the union of the lateral faces of $\Omega$. Let $T > 0$ be the final time. For the diffusion operators, $\Delta$ will indicate the 3D Laplace operator in $\Omega$, while $\Delta_\Gamma$ the Laplace-Beltrami operator on $\Gamma$ (which in this case coincides with the 2D Laplacian since $\Gamma$ is flat). We recall that the morpho-chemical DIB-model has been originally introduced in ^[1] on a 2D rectangular spatial domain, say $Q \subset \mathbb{R}^2$, in the framework of reaction-diffusion PDEs, with the original feature of coupling one equation for the morphology $\eta(x, y, t)$ with one for the surface chemistry $\theta(x, y, t)$, and takes the form:

where $\mathbf n: \partial Q \rightarrow \mathbb{R}^2$ is the outward unit normal vector on $\partial Q$. Essentially, here we consider that this system lives on the bottom face of the cube, i.e., $Q \equiv \Gamma$ and suitable modifications of the source terms are needed to describe the new physico-chemical features, as explained in the following. Concerning the unknown variables, also here we suppose that on the (flat) surface $\Gamma$, $\eta \in \mathbb{R}$ expresses the instantaneous increment of the electrodeposit thickness during an electrochemical process and $\theta \in [0, 1]$ is the surface coverage with an adsorbed chemical species, that influences the electrodepostion and corrosion processes. In the last ten years several theoretical and numerical results have been published on the 2D-DIB model (2.1), as reported in the Introduction and in the bibliography. Moreover, a wide range of experimental demonstrations of electrochemical pattern formation, has been documented. In most of the papers, the existence of the so-called Turing patterns has been proved theoretically, demostrated by numerical simulations and validated through comparison with experiments, also in the case of 3D surfaces (see ^[4]). Hence, to account for the more realistic presence of the battery electrolyte, the bulk-surface (BS)-DIB model seeks to find four functions: $b, q:\Omega\times [0, T] \rightarrow \mathbb{R}$ in the bulk and $\eta, \theta:\Gamma\times [0, T] \rightarrow \mathbb{R}$ on the surface that fulfil the following systems of PDEs:

and are coupled through the boundary conditions (BCs) for the bulk variables $b$ and $q$ given by

Moreover, zero Neumann BCs for (2.2b) for the surface variables $\eta$ and $\theta$ are considered, inherited by the original 2D-DIB model (2.1)

The initial conditions are prescribed by:

In (2.2), $d_\Omega > 0$ and $d_\Gamma > 0$ are adimensional diffusion coefficients. It is worth noting that the constants $\psi_\eta \geq 0$ and $\psi_\theta \geq 0$ are responsible of the strength of the BS coupling, in fact when they are equal to zero the two systems in (2.2) will be independent and the 2D DIB model (2.1) must be recovered by Eqs (2.2b)–(2.4).

In the bulk Eq (2.2a), $b(\mathbf x, t)$ represents the concentration of the electroactive cations (precursors of metal that is electrodeposited during the recharge cycle of the battery), present exclusively in the bulk and $q(\mathbf x, t)$ represents the bulk concentration of an additive species that is adsorbed at the cathode as a way of controlling shape change with the coverage degree expressed by the variable $\theta$ in (2.2b). Concerning the source terms for the bulk species $b$ and $q$, the kinetics $f_1, f_2$ are defined by:

where $b_0 = b_{bulk}\in \mathbb{R}$ and $q_0 = q_{bulk} \in \mathbb{R}$, represent the "bulk concentrations" and $k_b, k_q \geq 0$ are reaction rates.

For the surface variables $\eta$ and $\theta$, in (2.2b) the kinetics $f_3$ and $f_4$ are defined by:

and correspond to the source terms of the DIB model (2.1), but this time including the bulk contributions, as described in details below and such that

see ^[1]. In (2.6)–(2.9), the model parameters $k_b, k_q, k_2, k_3, \rho, \alpha, A_1, A_2, B, C, D$ are positive constants, and $\gamma \geq 0$.

2.1.   Physical meaning of source terms

In the kinetics $f_1, f_2$:

● $b_0 = b_{bulk}\in \mathbb{R}$ and $q_0 = q_{bulk} \in \mathbb{R}$, represent the "bulk concentrations", that prevail at equilibrium, when the bulk is homogeneous. The physical meaning of the terms $(b-b_0)$ and $(q-q_0)$, with $k_b, k_q \geq 0$ is first-order homogeneous reaction kinetics describing the tendency of the reagent to re-establish the equilibrium concentration. This can be considered a very simple model of a situation in which $b$, $q$ are the concentrations of the species involved in the electrodic reaction in the electroactive form, that is generated by the decomposition of some precursor (e.g., metallic ion with a ligand that keeps the ion in solution in non-electroactive form, from which the electroactive species forms by decomposition of the complexed one). A lower-than-equilibrium local concentration of $b$, $q$ (e.g., by cathodic consumption) generates new $b, q$ (by decomposition of the complexed form); while a higher-than-equilibrium concentration (e.g., by anodic injection) generates a consumption (e.g., by reaction with the ligand, yielding the non-electroactive form).

For aim of completeness, here we recall and update the meaning of the kinetics of the DIB modeling (i.e., $f, g, f_3, f_4$). The physical meaning of the terms in $f_3$ in (2.8) can be described as follows:

● As in ^[1], the Butler-Volmer type electrokinetic term ${A_1 (1-\theta) \eta}$ is the charge-transfer rate at sites free from adsorbates, but here this is scaled by the concentration of the electroactive species at the surface $b_{|\Gamma}$. This corresponds to first-order phenomenological kinetics and it includes naturally mass-trasport effects from the bulk to the surface (i.e., the mass-transport of the electroactive species $b$ present in the bulk and reacting electrochemically at the surface).

● The cubic term $-A_2\eta^3$ in the original DIB model ^[1] represents collectively all "hindrances" to $\eta$ (growth) resulting from "the establishment of high values of $\eta$". The most straightforward interpretation of such hindrances is mass-transport limitation that–in a Gileadi-type framework–can be approximately accounted for with a negative cubic correction to the Ⅰ-Ⅴ curve. As expounded above, in the new bulk-surface context mass-transport limitations can be naturally accounted for by multiplying the term linear in $\eta$ by the surface concentration $b_{|\Gamma}$ of the reactant $b$ present in the bulk and describing the electroactive species. Nevertheless, in the bulk-surface version of DIB it is worth retaining the cubic term in $\eta$, because this can account for other hindrances to metal growth (i.e., beyond mass-transport from the bulk to the reactive surface) appearing at high metal plating rates, such as cathodic passivation.

● The term $-B(\theta-\alpha)$ quantifies the effect of adsorbates on the electrodeposition rate. The parameter $0 < \alpha \leq 1$ takes into account the fact that adsorbates can have both inhibiting and enhancing effects on the growth rate.

● Concerning the equation for the chemical coverage $\theta$, in the original DIB model (2.1) the source term can be regarded as

and features adsorption (parameter $C$) and desorption (parameter $D$) terms including both chemical (expanded to second order) and electrochemical (first order) contributions. Here, the physical meaning of the modified kinetics $f_4$ in (2.2b) is that–coherently with the Langmuir adsorption model with monomolecular adsorption reaction–the adsorption term is directly proportional to the amount of the bulk species. In the case of a heterogeneous bulk phase, the relevant value of the bulk form is that in a neighborhood of the surface (commonly, referred to as the "catholyte") that is given by $q_{|_\Gamma}$.

2.2.   Physical meaning of boundary conditions

The coupling BCs for the bulk equations at the interface between the "growing surface $\Gamma$" and the "bulk $\Omega$" can thus naturally be written as the first two equations in (2.3) indicating that the flux of bulk species to the surface is opposite to their consumption rates at the surface. More specifically, the physical meaning of the form of the first two equations in (2.3) is that, even though $b$ and $q$ coming from the bulk are consumed at the interface (i.e., in correspondence of their values $b_{|\Gamma}$ and $q_{|\Gamma}$) to yield $\eta$ and $\theta$ only in a specific term of $f_3$ and $f_4$, the negative terms of these equations have the effect of injecting $b$ and $q$ into the bulk. Cases in which $\theta > \alpha$ (i.e., the adsorbate enhances electrodeposition, e.g., by resonant tunnelling effects), can be regarded as a special case of interfacial $b$ consumption, accounted for through the BCs. Thus the net formation rate of $\eta$ and $\theta$ is proportional to the fluxes of $b$ and $q$ to the surface.

$\nabla b \cdot \mathbf n, \nabla q \cdot \mathbf n$ denote the gradients normal to the boundary (surface) $\Gamma$, while $\psi_\eta, \psi_\theta$ are constants the role of which is to adjust the dimensionality of the equations and the physical meaning of which is explained below. $\psi_\eta$ converts adatoms (i.e., the surface species generated by the consumption of $b$ at the surface) into morphological units (the quantities actually described by $\eta$). Thus one can write: $\eta = \psi_\eta b$ and $\psi_\eta$ can be regarded as a constant as far as the density of morphological units (ca. step) is proportional to the adatom density. $\psi_\theta$ expresses an isotherm, since it connects a bulk concentration ($q$) into a surface density ($\theta$), hence, in as far as the isotherm can be linearised, we can write: $\theta = \psi_\theta q$.

To synthesize, the BCs for $b$ and $q$ in (2.3) are:

ⅰ) non-linear coupling BCs on the surface $\Gamma$, implying coupling with $\eta$ and $\theta$;

ⅱ) Dirichlet BCs on the face of $\Omega$ opposite to $\Gamma$ that is located far enough from the "bottom face" where reaction takes place, so that concentration gradients induced by reactivity have died out here;

ⅲ) zero Neumann BCs (zero flux) on the residual faces of $\Omega$ (see Figure 1). The values of the bulk variables are thus set to their equilibrium values $b_0, q_0$ respectively.

For the morpho-chemical unknowns $\eta$ and $\theta$, the BCs for the 3D BS-DIB model on $\partial \Gamma$ are still zero Neumann BCs in (2.4).

2.3.   Stability in the absence of diffusion

First of all, given the physical meaning of the BS-DIB model (2.2)-(2.3), it would make sense to consider spatial domains of more general shape. However, since the BS-DIB model is being introduced here for the first time, we have confined the presentation to the case of a cubic domain, in order to focus on the exploration of Turing patterns. A necessary condition to investigate diffusion-driven or Turing instability in the new BS model is that at least a homogeneous equilibrium of the PDE system exists and in the absence of diffusion it must be stable. In this section, we present this analysis, while the complete derivation of Turing conditions guaranteeing the existence of pattern formation in the presence of diffusion is deferred to another paper.

Hence, if $F = (f_1, f_2, f_3, f_4)$ accounts for all source terms involved and $\boldsymbol \xi = (b, q, \eta, \theta)^T$, it is easy to show that

when in the BS-DIB model all parameters in (2.8)-(2.9) different from $(B, C)$ are kept fixed after ^[2] as follows:

In fact, $k_b, k_q$ do not influence the result in (2.10), $(B, C)$ are still considered as bifurcation parameters as e.g., in ^[2] and the $(\eta, \theta)$ components of the equilibrium $\xi^*$ coincide with the homogeneous equilibrium of the original 2D-DIB model. Here we analyze the stability of such equilibrium in the special case $\gamma = 0$, when the condition on the parameter $D$ in (2.11) boils down to

To study the arising of the diffusion-driven or Turing instability in the BS-DIB model, it is necessary to prove that the equilibrium (2.10) is stable in absence of diffusion. In fact, the model (2.2), deprived of diffusion and linearized around the equilibrium in (2.10), is

where $J$ is the Jacobian of the kinetics evaluated at the equilibrium (2.10). The matrix $J$ has the following structure

where, if $\gamma = 0$, the blocks of $J$ are as follows:

Thanks to the diagonal structure of $J_\Omega$ it holds that

It follows that two eigenvalues of $J$ are $\lambda_1 = -k_b < 0$ and $\lambda_2 = -k_q < 0$. We are left to determine when the eigenvalues of $J_\Gamma$ are negative. This happens if and only if $Trace J_\Gamma < 0$ and $\det J_\Gamma > 0$. Now:

and

We obtain the following result.

Theorem 1. If $\gamma = 0$, the equilibrium (2.10) is stable in the absence of diffusion if and only if

In addition, if $A_1$, $k_2$, $k_3$, $\alpha$ are as in (2.11), the condition (2.21) specializes to

Of course, as far as the conditions for Turing pattern formation are concerned, a specific study has to be carried out considering the full Jacobian $J$ of Eq (2.14). This is an important study topic in its own right, that–nevertheless–has no impact on the analysis presented in this research. Since treating this problem exhaustively would be beyond the scope of the present work, we leave it to a subsequent publication.

Our present approach is to solve numerically the BS-DIB model for the diffusion coefficients $d_\Omega = 1$ and $d_\Gamma = 20$ and for a representative selection of parameter couples $(B, C)$ generating the whole set of Turing pattern morphologies, as described in ^[15,16]. In fact, this value of $d_\Gamma$ was that used in ^[2,15] to build a Turing region for the 2D-DIB model (2.1), that we will use as reference to investigate numerically the Turing pattern formation for the new BS-DIB model (2.2). Moreover, we shall fix the model parameters in the bulk as

such that

such that the spatially homogeneous equilibrium for the system (2.2)-(2.3) is

that is stable in absence of diffusion if $B > 20$ and $C > 2.5$. To approximate Turing patterns as steady state solutions of the BS-DIB model, we choose the initial conditions as:

where $r_\eta$ and $r_\theta$ are random spatial data defined as

where $\eta^*$ and $\theta^*$ are defined in (2.24) and $\texttt{rand}$ is the Matlab built-in function for generating uniform random numbers in $(0, 1)$.

In this scenario, our aim is to tune the coupling parameters $\psi_\eta, \psi_\theta$ to study from the numerical point of view the effect of coupling with the bulk on the morphological structure of the Turing patterns in the classes studied in ^[15]. We recall that the numerical approximation of RDSs in 3D is not straightforward because the pattern requires a very fine 3D mesh that provides sufficient spatial resolution and a long time integration to reach the asymptotic steady state. We shall apply the BS-VEM method studied in ^[8], for the space discretization of the 3D domain and surface and the IMEX Euler method as time solver. Hence, in the 3D case, if the cubic domain is approximated with a Cartesian grid, at least a million of unknowns at each time iteration are required. The usual implementation will thus end up to a sequence of linear systems where the coefficient matrix for each species is sparse, but prohibitively large. An efficient alternative to deal with this issue is the Matrix-Oriented Finite Element Method (MO-FEM) ^[12,13] where, thanks to the Cartesian structure of the numerical grid, the fully discrete problem is transformed to a sequence of Sylvester matrix equations, that can be solved efficiently in the spectral space. However, since the BS-DIB model produces spatial patterns only in a neighborhood of the surface $\Gamma$, we devise a tailor-made graded polyhedral bulk-surface mesh where the BS-VEM proves to be a competitive alternative, since such a graded mesh avoids unnecessary refinement (and degrees of freedom) away from the surface $\Gamma$. One of the major advantages of this choice is that the BS-VEM, thanks to the flexibility of polyhedral meshes, can still be used on domains of general shape, where MO techniques might not apply.

3.   The bulk surface virtual element method for the BS-DIB model

Formulating a BS-VEM ^[8] for the model (2.2) requires several steps. We start by rewriting the model (2.2) in such a way that the BCs lend themselves to a BS-VEM discretization.

3.1.   Step 1: rewriting the model with homogeneous boundary conditions

In the presence of non-zero Dirichlet BCs, it is well-known ^[17] that it is first necessary to rewrite the PDE problem in such a way that the Dirichlet conditions are homogeneous. To this end, we define the following auxiliary variables and kinetics:

With the above definitions, the model becomes

which, this time, is conveniently endowed with completely homogeneous BCs:

3.2.   Step 2: weak formulation

To write a discrete formulation of the auxiliary problem (3.7)-(3.8), we define the space of trivariate spatial functions that ensure the well-posedness of the model (3.7) and fulfil the BCs (3.8):

The dual space of $H^1_B(\Omega)$ will be denoted by $H^{-1}_B(\Omega)$. Following ^[8], the weak formulation of (3.7)-(3.8) is: find $\widetilde{b}, \widetilde{q} \in L^2([0, T]; H^1_B(\Omega))$ and $\eta, \theta \in L^2([0, T]; H^1(\Gamma))$ with $\dot{ \widetilde{b}}, \dot{ \widetilde{q}} \in L^2([0, T]; H^{-1}_B(\Omega))$ and $\dot{\eta}, \dot{\theta} \in L^2([0, T]; H^{-1}(\Gamma))$ such that

for all $\varphi \in L^2([0, T]; H^1_B(\Omega))$ and $\phi \in L^2([0, T]; H^1(\Gamma))$. The weak formulation (3.10) is well-posed for sufficiently short times thanks to the following result.

Lemma 1 (Well-posedness and stability estimates for the weak formulation). There exist $c, T > 0$ depending on the parameters of the model (2.2) such that the solution $(\widetilde{b}, \widetilde{q}, \eta, \theta)$ of the weak problem (3.10) exists and is unique for $t \in [0, T]$ and fulfills the estimates

Proof. The weak problem (3.10) falls in the class of BS-RDSs considered in ^[8], with the difference that in the BS-DIB model considered here, the kinetics $f_1, f_2, f_3, f_4$ are polynomials and are therefore only locally Lipschitz, rather that globally Lipschitz as in ^[8]. This lemma is thus proven exactly as [8, Lemma 4.1], with the difference that the resulting estimates (3.11)-(3.12) hold only for sufficiently small final time $T$. □

Remark 1 (Long time existence). Lemma 1 proves existence of a weak solution only for short times. When the kinetics are nonlinear as in the BS-DIB model (2.2), one way to prove long time existence is proving that the model possesses an invariant region, see ^[18]. The existence and the size of an invariant region might depend on the parameters of the considered BS-RDS model. For the BS-DIB model (2.2), a study of its invariant regions -if any- is one of our future research directions.

3.3.   Step 3: spatially discrete formulation

We will now describe the spatial formulation obtained by the BS-VEM following ^[8]. The choice of the BS-VEM to solve the model (2.2) is motivated by the possibility of using graded meshes that make the BS-VEM particularly competitive in this case by avoiding unnecessarily refinement away from the surface $\Gamma$. The choice of a convenient mesh will be illustrated in the next Section. For now, we illustrate the BS-VEM for arbitrary meshes.

Let us decompose the bulk $\Omega$ as the union of non-overlapping polyhedra,

If $\mathcal{F}_f$ is the set of the faces of $\mathcal{E}_h$ that are contained in $\Gamma$, then we can write

For a face $F\in\mathcal{F}_h$, the boundary space $\mathbb{B}(\partial F)$ is defined by

The preliminary space of a face $F$ is defined by

The $H^1$ projector on faces $\Pi^\nabla_F: \widetilde{\mathbb{V}}(F) \rightarrow \mathbb{P}_1(F)$ is defined, for any $v \in \widetilde{\mathbb{V}}(F)$ by

Then, the enhanced VEM space on the face $F$ is defined by

For a polyhedron $E\in\mathcal{E}_h$, the boundary space $\mathbb{B}(\partial E)$ is defined by

At this point, the preliminary VEM space on $E$ is defined by

The $H^1$ projector $\Pi^\nabla_E: \widetilde{\mathbb{V}}(E) \rightarrow \mathbb{P}_1(E)$ on the polyhedron $E$ is defined, for each $u\in \widetilde{\mathbb{V}}(E)$, by

Finally, the enhanced VEM space on the polyhedron $E$ is defined by

It is well-known that the degrees of freedom in $\mathbb{V}(F)$ and $\mathbb{V}(E)$ are the pointwise values on vertexes, see ^[19]. The global VEM spaces are defined by matching the degrees of freedom across elements. To this end, let $\mathbb{S}_\Gamma$ and $\mathbb{S}_\Omega$ be the 1-skeleton of $\Gamma$ and the 2-skeleton of $\Omega$, respectively, defined by

The global VEM spaces $\mathbb{V}_\Gamma$ and $\mathbb{V}_\Omega$ are then defined as

Notice that the space $\mathbb{V}_\Omega$ reflects the homogeneous Dirichlet BCs of the continuous counterpart $H^1_B(\Omega)$. To obtain a spatially discrete counterpart of the weak formulation (3.10), we need suitable discrete bilinear forms. Following ^[8], for all $F\in\mathcal{F}_h$, $E\in\mathcal{E}_h$, $v, w \in \mathbb{V}(F)$ and $u, z\in\mathbb{V}(E)$, we define

where $h_F$ and $h_E$ are the diameters of $F$ and $E$, respectively. In (3.24)–(3.27), the simplest form of stabilization was chosen, the so-called dofi-dofi stabilization ^[20], for two reasons. First, this choice is the most implementation-friendly, see ^[21]. Second, the dofi-dofi stabilization is known to perform well when the mesh is composed of elements that are not distorted and with very regular shapes ^[20], which is exactly our case as we will see.

Let $m_h^\Gamma, a_h^\Gamma : \mathbb{V}_\Gamma \times \mathbb{V}_\Gamma \rightarrow \mathbb{R}$ and $m_h^\Omega, a_h^\Omega : \mathbb{V}_\Omega \times \mathbb{V}_\Omega \rightarrow \mathbb{R}$ be the corresponding global forms. Furthermore, let $I_\Gamma: \mathcal{C}^0(\Gamma) \rightarrow \mathbb{V}_\Gamma$ and $I_\Omega: \mathcal{C}^0(\Omega) \rightarrow \mathbb{V}_\Omega$ be the Lagrangian interpolant operators. The spatially discrete formulation is finally given by: find $B, Q: \mathbb{V}_\Omega \times [0, T]$ and $\Lambda, \Theta : \mathbb{V}_\Gamma \times [0, T] \rightarrow \mathbb{R}$ such that

for all $\Phi : \mathbb{V}_\Omega \times [0, T] \rightarrow \mathbb{R}$ and $\Psi : \mathbb{V}_\Gamma \times [0, T] \rightarrow \mathbb{R}$. If $N_\Gamma : = \dim \mathbb{V}_\Gamma$ and $N_\Omega : = \dim \mathbb{V}_\Omega$, let $\{\psi_i\}_{i = 1}^{N_\Gamma}$ and $\{\varphi_i\}_{i = 1}^{N_\Omega}$ be the Lagrangian bases of $\mathbb{V}_\Gamma$ and $\mathbb{V}_\Omega$, respectively. The spatially discrete formulation (3.28) is well-posed for sufficiently short times thanks to the following result.

Lemma 2 (Well-posedness and stability estimates for the spatially discrete formulation). There exist $c, T > 0$ depending on the parameters of the model (2.2) and on the shape regularity of the mesh such that the solution $(B, Q, \Lambda, \Theta)$ of the spatially discrete problem (3.28) exists and is unique for $t \in [0, T]$ and fulfills the estimates

Proof. This result was proven in ^[8] in the case when the kinetics $f_1, f_2, f_3, f_4$ are globally Lipschitz. Here, $f_1, f_2, f_3, f_4$ are only locally Lipschitz, therefore the resulting estimates (3.29)-(3.30) hold only for sufficiently small final time $T$. □

Remark 2 (Long time existence for the spatially discrete solution). Lemma 2 proves existence of a spatially discrete solution only for short times. To prove long time existence, a sufficient condition would not only involve an invariant region as for the continuous problem (2.2), see Remark 1, but also a spatial method that preserves the invariant regions of the continuous model, see for instance ^[22] for the case of surface-only RDSs. To the best of our knowledge, a BS-VEM that preserves invariant regions under discretization is an open problem. It must be noted, however, that in our numerical experiments in Section 4, long time numerical solutions are found.

We express the numerical solution $(B, Q, \Lambda, \Theta)$ in the Lagrange bases:

where $b_i(t)$, $q_i(t)$, $\lambda_i(t)$, $\theta_i(t)$ are unknown time-dependent coefficients, which are collected in column vectors $\mathbf b(t), \mathbf q(t) \in \mathbb{R}^{N_\Omega}$, $\boldsymbol \eta(t), \boldsymbol \theta(t) \in \mathbb{R}^{N_\Gamma}$. Following ^[8], we substitute (3.31)–(3.34) into the spatially discrete formulation (3.28), and we obtain the following ODE system in vector form:

where the stiffness matrices $A_\Omega \in\mathbb{R}^{N_\Omega \times N_\Omega}$, $A_\Gamma \in\mathbb{R}^{N_\Gamma \times N_\Gamma}$, the lumped mass matrices $M_\Omega \in\mathbb{R}^{N_\Omega \times N_\Omega}$, $M_\Gamma\in\mathbb{R}^{N_\Gamma \times N_\Gamma}$ and the reduction matrix $R\in\mathbb{R}^{N_\Omega \times N_\Gamma}$ are defined as follows:

where $I_{N_\Gamma}$ is the identity of dimension $N_\Gamma$.

3.4.   Step 4: fully discrete formulation in vector form

Following ^[8], we discretize in time the ODE system (3.35) with the IMEX Euler scheme, which is first-order accurate. Let $\tau > 0$ be the time stepsize and let $N_T = \left\lceil \frac{T}{\tau}\right\rceil$ be the number of timesteps. For all $n = 0, \dots, N_T - 1$, the fully discrete solution $(\mathbf b^{(n)}, \mathbf q^{(n)}, \boldsymbol \eta^{(n)}, \boldsymbol \theta^{(n)})$ at time $t_n : = n \tau$ is found as follows:

where

The fully discrete formulation (3.39) is composed of four linear algebraic systems that can be solved independently of each other at each time step. Of these four linear systems, two have dimension $N_\Omega$, while the other two have dimension $N_\Gamma$. If the cube $\Omega$ is discretised with a Cartesian mesh with $N_x \in\mathbb{N}$ gridpoints along each dimension, then $N_\Omega = N_x^2(N_x-1)$ due to the boundary conditions, which makes the linear systems in (3.39) computationally expensive to solve. An extremely efficient approach to address this issue is the so-called Matrix-Oriented Finite Element Method (MO-FEM) ^[12,13], which exploits the Cartesian structure of the grid to translate the linear systems in (3.39) into tensor equations of much lower size. We will show a numerical solution to our problem carried out with MO-FEM in Section 4. However, we will show that, given the particular nature of the considered PDE problem, an even more efficient solver is given by the BS-VEM on a bespoke mesh.

3.5.   Bespoke BS-VEM mesh for the BS-DIB model

Since the domain of the model problem (2.2) is a cube, then it would be natural to choose an efficient numerical method that exploits the structure of Cartesian grids, such as the Matrix-Oriented Finite Element Method (MO-FEM) ^[12,13]. This gives us the opportunity to show the competitiveness of BS-VEM in solving the BS-DIB model (2.2) and thus further motivates the choice of a cube as spatial domain. For this reason, we choose the cubic spatial domain $\Omega = [0, L]^3$ with $L = 50$. We choose the model parameters, the final time and and the timestep as follows: $B = 66$, $C = 3$, $A_2 = 1$, $\gamma = 0.2$, $\psi_\eta = \psi_\theta = 0.2$, $T = 50$ and $\tau = 2$e-3. The other parameter values are as in (2.11). To demonstrate that the numerical solution tends to a stationary Turing pattern, as for instance in ^[2] we show, as an indicator, that the increment $\|\eta(t_{n+1}) - \eta(t_n)\|_{L^2(\Omega)}$ of the $\eta$ component of the numerical solution tends to zero.

In this section, we consider a Cartesian grid with $N_x = 129$, i.e., of $129\times 129 \times 129 \approx 2.15$e+6 equally spaced nodes (corresponding to $128 \times 128 \times 128$ equally spaced intervals), and we show that MO-FEM produces the numerical solution shown in Figure 2. In keeping with the physical meaning of the parameter choice, we can observe that the bulk components $(b, q)$ exhibit spatial patterns only in the proximity of the surface $\Gamma$, and become approximately constant away from $\Gamma$. This suggests that a uniform Cartesian cubic grid would be unnecessarily fine away from $\Gamma$. For this reason, we apply the BS-VEM with a graded cubic mesh that is highly refined close to $\Gamma$ and gradually becomes coarser as the distance from $\Gamma$ increases. Such graded polyhedral mesh is depicted in Figure 3(a). This grid is composed of two layers of $128\times 128$ cubic elements close to $\Gamma$, and five layers of gradually larger "false cubes", where the size of such false cubes doubles at each layer. For this reason, in the first two layers (closest to $\Gamma$) the number of cubes along each dimension is $128$, a power of $2$. One of the aforementioned false cubes is depicted Figure 3(b) and is actually an ennahedron: a polyhedron with nine faces and thirteen vertices. Specifically, such ennahedron is a cube the bottom face of which is split into four smaller square faces as shown in Figure 3(b). The proposed graded mesh provides a discretisation for the surface $\Gamma$ as refined as the Cartesian grid used for the MO-FEM ($129\times 129$ nodes on $\Gamma$) and, at the same time has much less nodes (approximately $5.56$e+4 versus $2.15e$+6), resulting in a discrete problem of much smaller size and shorter computational times (approximately 89 minutes for MO-FEM and 39 minutes for BS-VEM) when the linear systems are solved with Matlab's built-in direct solver \ (backslash). It must be noted that, if we applied BS-VEM on the same Cartesian mesh used for MO-FEM, the latter method would be much quicker as it is specifically designed to exploit the tensor structure of Cartesian grids. Specifically, BS-VEM would require prohibitively large computational times and memory, at least for our limited machine. On the other hand, MO-FEM (nor classical FEM) cannot be applied on graded or adaptive meshes with hanging nodes. All in all, the reduced number of degrees of freedom of the BS-VEM mesh outweights the efficiency of MO-FEM on Cartesian grids, resulting in shorter computational times. On domains of more general shapes, where tensorized methods like MO-FEM are not generally applicable, BS-VEM on a graded mesh would be even more competitive.

By comparing Figure 2 and Simulation D3 in Section 4 we observe that, for the BS-DIB model (2.2), MO-FEM and BS-VEM produce patterns of the same morphological class, i.e., reversed spots. The positioning of the patterns, however, is different, and that is to be expected because RDSs are highly sensitive to initial conditions. This implies that a small perturbation resulting from applying a different numerical method on equal initial conditions (as in (2.25)) is enough to significantly affect the positioning of the asymptotic spatial patterns, but not their morphological class ^[23]. Having ascertained the competitiveness of BS-VEM in solving the BS-DIB model, all the numerical experiments in Section 4 will be carried out using BS-VEM.

4.   Numerical experiments

We shall present eight numerical experiments to compare the DIB model (2.1) with the novel BS-DIB model (2.2). The eight experiments differ from each other for appropriate choices of the model parameters. The first four experiments, called T1 through T4, suggest that the BS-DIB model (2.2) has a larger Turing space than the DIB model (2.1), meaning that on equal parameters, the BS-DIB model might exhibit Turing patterns while the DIB model does not. The latter four experiments, called D1 through D4, explore the effect of BS coupling on the morphological class of Turing patterns, meaning that on equal parameters, the DIB (2.1) and BS-DIB (2.2) models can exhibit Turing patterns of different morphological classes. All the experiments are carried out on a cubic domain of edge length $L = 50$ on the polyhedral mesh described in Section 3.5. The final time $T$ and the timestep $\tau$ also differ for each experiment according to the stiffness of the problem and the timescale of the dynamics. A recap of the numerical experiments and the respective parameters is given in Table 1. As discussed in Section 3.5, we recall that the BS-VEM on a suitable graded mesh is more computationally efficient than MO-FEM in solving the BS-DIB model (2.2). For this reason, all the experiments presented in this Section will be carried out with BS-VEM. It must be noted that the dynamics of both the DIB (2.1) and BS-DIB (2.2) models are affected by the parameters, therefore the convergence towards the asymptotic steady state can be more or less quick depending on the model parameters. For this reason, also in the next experiments, we present a plot of the increment $\|\eta(t_{n+1})-\eta(t_n)\|_{L^2(\Gamma)}$ over time. We assume that if the increment at the chosen final time $T$ is substantially lower than the peak increment corresponding to the end of the reactivity transient regime ^[24], then the final time solution is a good approximation of the asymptotic steady state ^[25].

4.1.   Experiments T1–T4: exploration of the Turing space

In this Section, we compare the DIB (2.1) and BS-DIB (2.2) models in four experiments, called T1 through T4, where the model parameters $(C, B)$ are chosen outside -but close to- the Turing space of the DIB model (2.1), i.e., the region in the parameter space $(C, B)$ where the DIB model undergoes Turing instability and exhibits patterns, see ^[15]. Therefore, the DIB model does not produce patterns in these experiments as expected, and the numerical solutions early attain the stable equilibrium $\eta^* = 0, \theta^* = \alpha$ in (2.24) of which the initial condition (2.25) is a small perturbation. On the other hand, the BS-DIB model (2.2) exhibits spatial patterns in all experiments T1–T4. This suggests that the BS-DIB model has a larger Turing space than the DIB model. As mentioned before, a theoretical analysis of the Turing instability for the BS-DIB model is outside the scope of this work and is one of our current research directions. In particular, the numerical solutions $\eta$ and $\theta$ on the surface $\Gamma$ of the BS-DIB model have the following behaviors:

T1. Slow convergence to a thin worm pattern, where the worms slowly merge into longer worms over time, see Figure 4;

T2. Slow convergence to a stripe pattern. The transient solution exhibits spots that slowly merge into stripes over time, see Figure 5;

T3. Convergence to a reversed spots (holes) pattern, see Figure 6;

T4. Quick convergence to a homogeneous steady state different from (2.24), see Figure 7.

4.2.   Experiments D1–D4: effect of bulk-surface coupling on pattern formation

In this Section, we compare the DIB model (2.2) and the BS-DIB model (2.2) in four more experiments, called D1 through D4, that are devised to study the effect of coupling in the 3D model, measured by increasing the coefficients $\psi_\eta = \psi_\theta$ in (2.3). As opposed to Experiments T1–T4 in the previous Section, the parameters here are chosen inside the Turing space of the 2D-DIB model (2.1), in particular those for the reference labyrinth and reversed spots Turing patterns studied in ^[15]. The proposed simulations explore the effect of bulk-surface coupling on the morphological class of Turing patterns in the BS-DIB model, as follows.

In particular, we will show that the labyrinth breaks into more fragmented structures such as worms and holes/reversed spots: the stronger the coupling, i.e., the higher $\psi_\eta = \psi_\theta$, the more fragmented the pattern, up to the limit case when the BS-DIB pattern is composed by holes/reversed spots in Experiment D3. Moreover, Experiments D1–D3 indicate that the range and size of the bulk patterns increases with $\psi_\eta = \psi_\theta$ as illustrated below. Experiment D4 instead will show that when the 2D DIB model (2.1) exhibits holes/reversed spots, the effect of coupling is only a reduction in size of the holes.

To sum up, the outcome of Experiments D1 through D4 is as follows:

D1. For $B = 66, C = 3$ as in ^[15], $\psi_\eta = \psi_\theta = 0.1$, the DIB model (2.1) at $T = 200$ attains a labyrinth, while the 3D BS-DIB model (2.2) converges to a reversed spots & worms pattern, see Figure 8, as it is also evident by the behaviour of the $\eta$-increment in Figure 8;

D2. Same $(B, C)$ as in Experiment D1, $\psi_\eta = \psi_\theta = 0.15$, the 3D BS-DIB model (2.2) attains at $T = 50$ the same morphological structure, but with shorter worms than in Experiment D1, see Figure 9;

D3. Same $(B, C)$ as in Experiment D1, $\psi_\eta = \psi_\theta = 0.2$, the 3D BS-DIB model (2.2) converges to a holes/reversed spots pattern, see Figure 10;

D4. For $B = 30, C = 3$ as in ^[15], $\psi_\eta = \psi_\theta = 0.1$, at $T = 200$ both the 2D-DIB (2.1) and the 3D- BS-DIB (2.2) models converge to a holes/reversed spots pattern, where the holes are larger for the original DIB model, see Figure 11.

These results and more numerical details are reported in Table 1.

Moreover, for the Experiments D1–D3, we compare a $(x, z)$ section of the bulk component $b(x, y, z, T)$ at the final time of integration for $y = 12.5$. In Figure 12, we show that by increasing the coupling parameters $\psi_\eta = \psi_\theta$ (left to right), the amplitude of the bulk-pattern increases (as indicated by the colorbars) and length along the $z$-direction, that is the morphology structures venture more significantly into the electrolyte.

5.   Conclusions

We have introduced a BS-RD model in 3D, which we have called BS-DIB model, for electrodeposition and battery modeling. Compared to the previous DIB model in 2D, the new model fully accounts for the non-uniform electrolyte concentration in a neighborhood of the electrodic surface. The two-way coupling between bulk and surface substantially influences the long-term behavior of the system and in particular the morphological class of the Turing patterns obtained as asymptotic steady state solutions. Specifically, we find that the bulk-surface coupling has two main effects. First, we observe empirically that the BS-DIB model possesses a large Turing region in the parameter space, compared to the 2D DIB model. Second, when the parameters are chosen in the Turing space of the DIB model, the BS-DIB model still exhibits spatial patterns, but of a different spatial structure, i.e., the bulk-surface coupling affects the morphological class of the attained patterns. A theoretical Turing instability analysis of the BS-DIB model is beyond the scope of this work. These aspects form part of our current investigations.

The BS-DIB model is posed on a cubic domain, so it lends itself to efficient numerical solvers specifically devised for Cartesian grids, such as the MO-FEM ^[12]. Moreover, since the BS-DIB model exhibits spatial patterns only in a neighborhood of the surface, we have adopted the BS-VEM on a graded mesh that is highly refined close to the surface and much coarser away from the surface. Such a graded mesh combines the advantages of (ⅰ) being composed by equal elements of cubic shape, which significantly speeds up matrix assembly and improves matrix structure and (ⅱ) has far less degrees of freedom than a uniform Cartesian grid with the same level of refinement close to the surface. For this reason, the BS-VEM on a graded mesh proves to be more computationally efficient than the MO-FEM (see discussion in Section 3.5) and is thus the spatial method of choice throughout this work. Moreover, as opposed to the MO-FEM, which is confined to structured geometries such as Cartesian grids, the BS-VEM can handle domains of general shape, thereby facilitating the simulation of real case studies.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The work of MF was funded by Regione Puglia (Italy) through the research programme REFIN-Research for Innovation (protocol code 901D2CAA, project number UNISAL026). The work of IS is supported by PRIN2022 PNRR “BAT-MEN” (BATtery Modeling, Experiments & Numerics) Enhancing battery lifetime, Project code: P20228C2PP 001, CUP: F53D23010020001, funded by MIUR (Italian Ministry of University and Research) and European Union–NextGenerationEU. MF and IS are supported by the research project “Metodi avanzati per la risoluzione di PDEs su griglie strutturate, e non” (INdAM-GNCS project CUP E53C22001930001, U-UFMBAZ-2023-000205 21-02-2023). MF and IS are members of the Gruppo Nazionale Calcolo Scientifico-Istituto Nazionale di Alta Matematica (GNCS-INdAM).

Conflict of interest

The authors declare no conflicts of interest.