Homogenization based two-scale modelling of ionic transport in fluid saturated deformable porous media

TL;DR

This paper develops a homogenized two-scale model for ionic transport in deformable porous media, capturing electrochemical interactions and fluid-structure coupling, with numerical implementation demonstrating microstructure influence on effective properties.

Contribution

It introduces a homogenization approach for coupled ionic transport and deformation in porous media, including a numerical implementation in finite element software.

Findings

Homogenized coefficients depend on microstructure porosity.

The model separates electrolyte flow from solid deformation.

Microscopic responses can be reconstructed from macroscopic solutions.

Abstract

The paper deals with the homogenization of deformable porous media saturated by two-component electrolytes. The model relevant to the microscopic scale describes steady states of the medium while reflecting essential physical phenomena, namely electrochemical interactions in a dilute Newtonian solvent under assumptions of a small external electrostatic field and slow flow. The homogenization is applied to a linearized micromodel, whereby the thermodynamic equilibrium represents the reference state. Due to the dimensional analysis, scaling of the viscosity and electric permitivity is introduced, so that the limit model retains the characteristic length associated with the pore size and the electric double layer thickness. The homogenized model consists of two weakly coupled parts: the flow of the electrolyte can be solved in terms of a global pressure and streaming potentials of the two ions, independently of then the solid phase deformations which is computed afterwards for the fluid stress acting on pore walls. The two-scale model has been implemented in the SfePy finite element software. The numerical results show dependence of the homogenized coefficients on the microstructure porosity. By virtue of the corrector result of the homogenization, microscopic responses in a local representative cell can be reconstructed from the macroscopic solutions.

Tables

Table 1. Table 1 : Description of used parameters, source [ 4 ] .

Symbol	Quantity	Value	Unit
$e$	Electron charge	$1.6\times {10}^{-19}$	C
$k_B$	Boltzmann constant	$1.38\times {10}^{-23}$	J/K
$T$	Absolute temperature	$298$	K
$ℰ$	Dielectric constant	$6.93\times {10}^{-10}$	c/(mV)
${\eta }_f$	Dynamic viscosity of fluid	$1\times {10}^{-3}$	kg/(ms)
$D_1$	Diffusivity of 1st ionic species	$13.33\times {10}^{-10}$	m2/s
$D_2$	Diffusivity of 2nd ionic species	$20.32\times {10}^{-10}$	m2/s
$l$	Characteristic pore size	$1.0\times {10}^{-7}$	m
$c_c$	Characteristic concentration	$(6.02\times {10}^{24},6.02\times {10}^{26})$	particles/m3
${\Sigma }_c$	Surface charge density	$-0.129$	C/m2
$\Lambda$	Young modulus	$7.3\times {10}^9$	Pa

Equations

$${\mu }_{\alpha }={{\mu }_{\alpha }}^0+k_{B}T\ln c_{\alpha }+e{z}_{\alpha }\Psi ,$$

$$\frac{\partial c_{\alpha }}{\partial t}+\nabla \cdot \left({\mathbf{j}}_{\alpha }+\text{w}{c}_{\alpha }\right)=0\text{ in }{\Omega }_f,\alpha =1,2,$$

$${\mathbf{j}}_{\alpha }=-\frac{c_{\alpha }{D}_{\alpha }}{k_{B}T}\nabla {\mu }_{\alpha }=-\frac{c_{\alpha }{D}_{\alpha }}{k_{B}T}(\frac{k_{B}T}{c_{\alpha }}\nabla c_{\alpha }+e{z}_{\alpha }\nabla \Psi ),\alpha =1,2,$$

$${\mathbf{j}}_{\alpha }\cdot \text{n}=0\text{ on }\Gamma ,\alpha =1,2.$$

$$\mathcal{E}\nabla \cdot \text{E}=e\sum _{\beta =1}^2{z}_{\beta }{c}_{\beta }\mathrm{in}{\Omega }_f,\text{E}=-\nabla \Psi ,$$

$$-\nabla \cdot {\mathbf{\sigma }}_f=\text{f}\mathrm{in}{\Omega }_f,$$

$${\mathbf{\sigma }}_f=-p\text{I}+2{\eta }_f\mathbf{e}(\text{v})+{\mathbf{\tau }}_M,$$

$$\nabla \cdot \text{v}=0\mathrm{in}{\Omega }_f.$$

$$\nabla p-{\eta }_f{\nabla }^2\text{v}=\text{f}-e\sum _{\beta =1}^2{z}_{\beta }{c}_{\beta }\nabla \Psi \mathrm{in}{\Omega }_f.$$

$$-\nabla \cdot {\mathbf{\sigma }}_s=\text{f}\mathrm{in}{\Omega }_s,{\mathbf{\sigma }}_s=\text{A}\mathbf{e}(\text{u}),$$

$$\mathcal{E}\text{E}\cdot \text{n}=\Sigma \mathrm{on}\Gamma .$$

$$\text{w}=\text{v}.$$

$$\text{w}=0\mathrm{on}\Gamma ,$$

$${\mathbf{\sigma }}_f\cdot \text{n}={\mathbf{\sigma }}_s\cdot \text{n}\mathrm{on}\Gamma .$$

$$\begin{array}{r}Y=Y_s\cup Y_f\cup {\Gamma }_Y,Y_s=Y\setminus \overline{Y_f},{\Gamma }_Y=\overline{Y_s}\cap \overline{Y_f}.\end{array}$$

$${\nabla }^{\prime }=({\partial }_{x^{\prime }})=L_c({\partial }_x)=L_c\nabla .$$

$$p^{\epsilon }=\frac{p}{p_c},{\text{w}}^{\epsilon }=\frac{\text{w}}{v_c},{\Psi }^{\epsilon }=\frac{\Psi }{{\Psi }_c},c_{\alpha }^{\epsilon }=\frac{c_{\alpha }}{c_c},{\text{u}}^{\epsilon }=\frac{\text{u}}{u_c},$$

$${\Psi }^{\mathrm{ext},\ast }=\frac{{\Psi }^{\text{ext}}(x)}{{\Psi }_c},{\text{E}}^{\ast }=\nabla {\Psi }^{\mathrm{ext},\ast },{\Sigma }^{\ast }=\frac{\Sigma }{{\Sigma }_c},{\mathbf{f}}^{\ast }=\frac{\text{f}{L}_c}{p_c}.$$

$${\lambda }_D=\sqrt{\mathcal{E}\frac{{\Psi }_c}{e{c}_c}{\left(\sum _{\beta =1}^2{z}_{\beta }^2\right)}^{-1}}=\sqrt{\mathcal{E}{k}_{B}T{\left(e^2{c}_c\sum _{\beta =1}^2{z}_{\beta }^2\right)}^{-1}}=,$$

$${\epsilon }^2{\nabla }^2{\text{w}}^{\epsilon }-\nabla p^{\epsilon }$$

$$\nabla \cdot {\text{w}}^{\epsilon }$$

$${\epsilon }^2{\nabla }^2{\Psi }^{\epsilon }$$

$$\nabla \cdot \left({\text{Pe}}_{\alpha }{\text{w}}_{\alpha }^{\epsilon }{c}_{\alpha }^{\epsilon }+{\mathbf{j}}_{\alpha }^{\epsilon }\right)$$

$$-\nabla \cdot ({\text{A}}^{\ast }\mathbf{e}({\text{u}}^{\epsilon }))$$

$${\text{w}}^{\epsilon }$$

$${\mathbf{j}}_{\alpha }^{\epsilon }\cdot \text{n}$$

$$\epsilon \nabla {\Psi }^{\epsilon }\cdot \text{n}$$

$${\text{A}}^{\ast }\mathbf{e}({\text{u}}^{\epsilon })\cdot \text{n}$$

$${\mathbf{j}}_{\alpha }^{\epsilon }=$$

$${\mathbf{\sigma }}_f^{\epsilon }=$$

$$\begin{array}{rl}{c}_{\alpha }^{\epsilon }(x)=c_{\alpha }^{\mathrm{eq},\epsilon }(x)+\delta c_{\alpha }^{\epsilon }(x),& {\Psi }^{\epsilon }(x)={\Psi }^{\mathrm{eq},\epsilon }(x)+\delta {\Psi }^{\epsilon }(x),\\ {\text{w}}^{\epsilon }(x)={\text{w}}^{\mathrm{eq},\epsilon }(x)+\delta {\text{w}}^{\epsilon }(x),& p^{\epsilon }(x)=p^{\mathrm{eq},\epsilon }(x)+\delta p^{\epsilon }(x),\\ {\text{u}}^{\epsilon }(x)={\text{u}}^{\mathrm{eq},\epsilon }(x)+\delta {\text{u}}^{\epsilon }(x),& \end{array}$$

$$\begin{array}{rl}{c}_{\alpha }^{\mathrm{eq},\epsilon }(x)& =c_{\alpha }^b\exp (-z_{\beta }{\Psi }^{\mathrm{eq},\epsilon }(x)),\\ p^{\mathrm{eq},\epsilon }(x)& =\sum _{\beta =1}^2{c}_{\beta }^{\mathrm{eq},\epsilon }(x),\end{array}$$

$$\begin{array}{rl}{\epsilon }^2{\nabla }^2{\Psi }^{\mathrm{eq},\epsilon }=& \gamma \sum _{\beta =1}^2{z}_{\beta }{c}_{\beta }^b\exp (-z_{\beta }{\Psi }^{\mathrm{eq},\epsilon })\mathrm{in}{\Omega }_f^{\epsilon },\\ \epsilon \nabla {\Psi }^{\mathrm{eq},\epsilon }\cdot \text{n}=& -N_{\sigma }{\Sigma }^{\ast }\mathrm{on}{\Gamma }^{\epsilon },\end{array}$$

$$\sum _{\beta =1}^2{z}_{\beta }{c}_{\beta }^b=0.$$

$${\Psi }^{\mathrm{eq},\epsilon }(x)={\Psi }^{\mathrm{eq}}(y),c_{\alpha }^{\mathrm{eq},\epsilon }(x)=c_{\alpha }^{\mathrm{eq}}(y),$$

Full text

1 Introduction

Modelling the transport of an electrolyte solution through a porous medium (TEPM) is a multiscale nonlinear problem with obvious multiphysics features. The problem is of interest for a wide range of science fields, including geosciences, environmental engineering, physiology and tissue biomechanics, material science and namely chemical engineering. Moreover, there are challenging industrial applications related to energy storage (batteries), extraction of renewable energy from salinity differences, or corrosion of reinforced concrete structures. The technology-related areas of engineering and scientific research require a quantitative analysis of the TEPM using computational modelling which allows to capture influences of microstructure related phenomena on the macroscopic properties and behavior at the macroscopic level. In this context, the homogenization of periodic, or locally periodic structures is one of the most relevant modelling approaches which lead to efficient computational algorithms. On one hand, the upscaling procedure enables to compute macroscopic tensorial coefficients respecting a given microstructure, on the other hand the downscaling procedure enables to interpret the macroscopic response at the microscopic level. For both these procedures, the so-called characteristic responses which are obtained as solutions of the representative volume elements (RVE) are needed.

The aim of this paper is to develop a two-scale computational model for the quasi-static transport of a two-component electrolyte solution trough a deformable porous medium, such that the upscaling and downscaling procedures allow for studying global and local effects in a response to the microstructure and material properties.

During the last decade, there appeared a significant body of literature devoted to the modelling of the TEPM. Here we comment only on those publications which, as we believe, are the most relevant and tightly related to the present work

As stated above, a remarkable part of the related research on the matter concerns the geosciences. One of the most recognized work in this field is the paper [19] which relies on the homogenization procedure to derive a macroscopic model of expansive clays composed of a charged solid phase saturated by an electrolyte solution. This microscopic model includes equations describing electro-hydrodynamics coupled with the equation governing the flow of the electrolyte solution, ion electrodiffusion and electric potential distribution. Then the asymptotic homogenization is used to derive a two-scale model of electrokinetic phenomena, such as the electro-osmotic flow driven by the streaming potential gradient, the electrophoretic motion of mobile charges and the swelling induced by the osmosis. This model was later revisited in paper [20] with a more focus on the rigorous homogenization procedure and its analysis. It should be noted that a similar problem, *i.e. *the transport of an N-component electrolyte solution through porous rigid body subjected to a static electric field, was also studied in [17], although no assumptions about the electric double layer were considered.

In biomechanical modelling, authors of [13] and [15], use a similar approach to study a bone fluid flows at two porosity levels in the cortical bone tissue. It is worth to note a possible application of that model in the studies of the mechanosensing, cf. [14], and bone remodeling. These issues were treated in [21], where Biot’s poroelastic theory applied to the three-dimensional anisotropic media in order to account for deformation induced fluid flows in osteonal matrix under the harmonic loading. Homogenization of the ionic exchange between the charged porous medium and the electrolyte solution was elaborated in [14], being motivated by mechanosensing. A two-scale one-dimensional model for horizontal electro-osmotic flows in a number of thin horizontal slits was proposed in [5]. Therein, the pressure gradient and a horizontal electrical field were recognized as the flow driving forces. Although this work is focused on one specific case and disregards any deformation of the solid part, it provides a useful insight to the homogenization of the electroosmotic law for different types of multi-component electrolytes.

Most of the works devoted to TEPM assume that an electrolyte saturates a rigid porous medium, see [2, 3], the deformation phenomenon or evolving porous structures were considered in a number of papers, see e.g. [23, 4]. Moreover, some other works, e.g. [6, 27, 29, 26] treating the fluid-structure interaction without any electro-osmotic, or electromechanical coupling established useful platforms for extensions of those particular models to account for the phenomena featuring the transport of electrolytes. Such an extension was reported in [4], which is motivated by the study of nuclear waste disposal. In this work, the well known system of equations governing the ionic transport and extend it by elasticity of the solid part was introduced. The coupling between fluid motion and deformation of the solid matrix was also explored earlier in the work [18]. Therein the authors show that by a suitable choice of time scale, the deformation of porous medium becomes only weakly coupled to the electrokinetic system, which is advantageous for the model implementation and numerical simulations.

For completeness, let us note, that a non-stationary case model of ionic transport consisting of Stokes, Nerst-Planc and Poisson systems of equations was reported in papers [24, 30] and [10], where the upscaling procedure was treated using the two-scale convergence. However, since the time response at multiple scales introduces further difficulties in the modelling, in our study, we account for steady state problems.

Although the ionic transport in porous structures is well-known problem, there are still some challenging issues deserving more attention. Most of the papers cited above concern the theoretical issues of the mathematical modelling, without numerical simulations. On contrary, our interest lies in the implementation of a physically correct homogenized model, such that the upscaling and downscaling procedures are available for 3D microstructures. For this purpose we consider the model treated in [4] and provide essential ingredients of the modelling, starting from the model definition at the microscopic level, pursuing the linearization which allows for using the two-scale homogenization.The derived the macroscopic model involves the effective medium parameters, such as poroelasticity coefficients, permeability, diffusivity and other coupling coefficients which satisfy the Onsager reciprocity relationships.

The paper is organized as follows. In Section 2 we provide a brief introduction in the physical phenomena and their mathematical descriptions which constitute the mathematical model of the two-component electrolyte solution transported in the deforming elastic skeleton. As the next step, in Section 3, we introduce the dimensionless form of the mathematical model which is subject to the linearization procedure. The homogenization is reported in Section 4; therein the principal results are presented, namely the local problems for the so-called corrector functions, formulae for computing the homogenized equations, and, finally, the macroscopic model. In Section 5, the effective coefficients relevant to both mathematical models are quantified for varying porosity the microstructure geometry. Furthermore, this section introduces the numerical solution of the macroscopic model and its recovery at the microscopic level. Section 6 summarizes the results of this paper. A is devoted to the unfolding operator. B clarifies the introduction of the scale parameter into the model.

Basic notations

Through the paper we shall adhere to the following notation. The position $x$ in the medium is specified through the coordinates $(x_{1},x_{2},x_{3})$ with respect to a Cartesian reference frame. We shall also use the microscopic (dilated) Cartesian reference system of coordinates $(y_{1},y_{2},y_{3})$. The gradients are employed, $\nabla_{x}=(\partial/\partial x_{i})$ and $\nabla_{y}=(\partial/\partial y_{i})$ alternatively. As usually, the vectors will be denoted by bold letters, for instance, ${\textbf{{u}}}(x)$ denotes the solid matrix displacement vector field depending on the spatial variable $x$. Moreover, the components of this vector will be denoted by $u_{i}$ for $i=1,...,3$, thus ${\textbf{{u}}}=(u_{i})$. The Einstein summation convention is used which stipulates implicitly that repeated indices are summed over. By ${\mathbb{R}}$ the real number set is denoted. The differential volume and surface elements are denoted by $\mathrm{\,dV}$ and $\mathrm{\,dS}$, respectively. Function spaces are introduced subsequently in the text.

2 Mathematical model

The porous medium occupies an open bounded domain $\Omega\subset{\mathbb{R}}^{d}$, where $d=3$. Without loss of generality, we may assume that the domain represent a specimen shaped as a block $\Omega=]0,L[^{d}$, which will enable us to impose periodic boundary conditions on $\partial\Omega$. According to the phases, $\Omega$ splits into the fluid $\Omega_{f}$ and solid $\Omega_{s}$ parts, whereby both $\Omega_{s}$ and $\Omega_{f}$ are connected domains. By $\Gamma=\partial\Omega_{s}\cap\partial\Omega_{f}$ we refer to the solid-fluid interface, and n designates the unit normal vector on $\Gamma$, being outward to $\Omega_{f}$. The subscripts $\sqcup_{s}$ and $\sqcup_{f}$ will be used through the rest of the text also to denote the constants and variables belonging to the respective phases.

2.1 Electrical double layer and the electrostatic potential of a phase

Due to the surface effects of charged skeleton $\Omega_{s}$, in domain $\Omega_{f}$ occupied by the electrolyte, the so-called electrical double layer (EDL) can be distinguished in the proximity of charged pore surfaces, were the ionic charge distribution is perturbed from the bulk. The EDL splits into two sub-layers, the Stern layer and the diffuse layer, The thickness of the EDL is related to the Debye parameter $\lambda_{D}$ which will be specified below. Within the EDL, the attraction is strong enough to influence particle movement. For a deeper physical insight we refer to [11].

While in the solid phase the electrostatic potential $\Psi_{s}$ is constant, the total electric potential $\Psi_{f}$ of the fluid phase is associated with the distribution of ions in the electrolyte. Considering the effects of the EDL, the electrical potential $\Psi_{f}$ in the fluid varies strongly with the distance from to pore surface, since the ions show the tendency to arrange themselves to minimize their free energy. These effects, as illustrated in Fig. 2, result in the Poisson-Boltzmann distribution of the electrical potential $\Psi_{\textrm{EDL}}$. In the bulk, *i.e. *away from the solid-fluid interface, the electrical potential $\Psi_{f}=\Psi_{\textrm{bulk}}$ attains a constant value. We can also introduce the electrostatic part $\Psi$ of the total potential in the fluid as the difference $\Psi=\Psi_{\textrm{EDL}}-\Psi_{\textrm{bulk}}$.

The fluid flow perturbs the electrostatic distribution of ions caused by the EDL and drags some ions in direction of the flow. This produces variations of the ionic concentrations along the flow direction, so that the so-called streaming potential $\Phi_{\alpha}$ can be defined for every ${\alpha}$-th ionic species.

If an external electrical field $\mathbf{E}^{\rm{ext}}$ is present, another type of electrical potential is distinguished. By a construction, the so-called exterior (affine) potential $\Psi^{\rm{ext}}=-x\cdot\mathbf{E}^{\rm{ext}}$ can be introduced, such that $\mathbf{E}^{\rm{ext}}=-\nabla\Psi^{\rm{ext}}$. The imposed field $\mathbf{E}^{\rm{ext}}$ can usually be considered small, as compared with the fields in the EDL, so that a linearization of the non-linear problem can be employed. We may summarize that the total potential $\Psi_{f}$ is given by the sum of electrostatic potential $\Psi$, streaming potentials $\Phi_{\alpha}$, and exterior potential $\Psi^{\rm{ext}}$, thus $\Psi_{f}=\Psi+\Phi_{\alpha}+\Psi^{\rm{ext}}$.

2.2 Processes in the fluid phase

The fluid is an electrolyte solution of the solvent and two ionic species with different valencies $z_{\alpha}$, ${\alpha=1,2},z_{1}=-1,z_{2}=+1$. The ionic transport in the fluid is driven by three phenomena: 1) convection of the solvent which is determined by a convective velocity w, 2) diffusion of the ${\alpha}$-th ionic species in the solvent characterized by the diffusivity $D_{\alpha}$ of the ${\alpha}$-th ionic species, and 3) motion of ions due to the electrical field. In what follows, we introduce the system of equations describing these processes in the pore space $\Omega_{f}$ filled by the electrolyte solution. All the electrochemical and mechanical constants involved in the model are defined in Tab. 1.

Each species (labeled by $\alpha$) dissolved in the electrolyte is associated with the electrochemical potential ${\mu_{\alpha}}$, which depends on the concentration $c_{\alpha}$, the electrostatic potential $\Psi$, and temperature $T$, such that , see [11],

[TABLE]

where ${\mu_{\alpha}}^{0}$ is the standard electrochemical potential expressed at infinite dilution, $k_{B}$ is Boltzmann constant, and $e$ is the elementary charge. The transport is restricted to the Eulerian mass conservation law,

[TABLE]

where w stands for convective velocity and the effects of diffusion and migration caused by an external electrical field are expressed by the migration-diffusion flux $\mathbf{j}_{\alpha}$ given by

[TABLE]

whereby, the no-flux condition holds on the solid-fluid interface $\Gamma$,

[TABLE]

The electrokinetics of the fluid phase is characterized by distribution of the electrostatic potential $\Psi$ which satisfies the Gauss-Poisson problem,

[TABLE]

where $\mathcal{E}$ is dielectric coefficient of the solvent (assumed to be constant), and $\mathbf{E}$ is the electric field.

Obviously, the fluid flow is governed by the Navier-Stokes equations involving the fluid velocity v, the hydrostatic pressure $p$. Since we consider slow flows only, the convective term can be neglected as well as the inertia term, thus, the fluid stress tensor $\boldsymbol{\boldsymbol{\sigma}}_{f}$ satisfies the equilibrium equation,

[TABLE]

where f is the external body force. The stress in the fluid phase is extended by the Maxwell 2nd order stress tensor $\boldsymbol{\tau}_{M}=\mathcal{E}\left({\textbf{{E}}}\otimes{\textbf{{E}}}-\frac{1}{2}|{\textbf{{E}}}|^{2}{\textbf{{I}}}\right)$ so that

[TABLE]

where $\mathbf{e}({{\textbf{{v}}}})=\frac{1}{2}\left(\nabla{\textbf{{v}}}+(\nabla{\textbf{{v}}})^{\textrm{T}}\right)$, $p$ is the fluid pressure and $\eta_{f}$ is the dynamic viscosity of the electrolyte. The flow is assumed to be incompressible, hence

[TABLE]

Upon substituting (7) into (6), the Poisson equation (5) allows to rewrite the term $\nabla\cdot{\boldsymbol{\tau}_{M}}$ in terms of $\nabla\Psi$, such that the modified Stokes problem is obtained,

[TABLE]

Let us note that if the porous medium is considered to be rigid, the influence of solid matrix on fluid phase can be expressed only by boundary conditions on $\Gamma$. However, in case of the deformable porous body, the model needs to be expanded by equations describing processes in solid phase.

2.3 Fluid solid interaction

Model of the solid phase

We consider solid skeleton of the porous medium is constituted by an elastic conducting material. Therefore, as stated above, we assume the electric potential $\Psi_{s}$ is constant in the whole skeleton $\Omega_{s}$, whereas a constant surface charge $\Sigma$ is distributed on its surface $\partial\Omega_{s}$. Under standard assumptions of linear elasticity, small displacement and deformation, the displacement field u satisfies

[TABLE]

where $\mathbf{e}({{\textbf{{u}}}})=\frac{1}{2}\left(\nabla{\textbf{{u}}}+(\nabla{\textbf{{u}}})^{\textrm{T}}\right)$ is the strain tensor, and the elasticity tensor ${\textbf{{A}}}=(A_{ijkl})$ is symmetric and positive definite, $A_{ijkl}=A_{jikl}=A_{klij}$.

Interface conditions

To complete the model of the fluid saturated porous medium, the interface conditions must be prescribed on $\Gamma$. It should be pointed out that the transmission condition on $\Gamma$ concerning the electric field or electric potential may be defined using two methods. While considering no-slip hydrodynamic condition on the solid surface, i.e. ${\textbf{{v}}}=\dot{\textbf{{u}}}$, the electric potential on $\Gamma$ can be introduced using the $\zeta$-potential. This parameter is widely used in definition of the EDL, but it is more related to the electrochemical properties of the system, [11]. However, to treat more general situations, a surface charge density $\Sigma$ proportional to the normal derivative of $\Psi$ is often used instead of prescribing the $\zeta$-potential. Since, in our setting, potential $\Psi_{s}$ is a constant, the surface charge $\Sigma$ given at $\Gamma$ constitutes the boundary condition,

[TABLE]

Concerning the mechanical interaction on $\Gamma$, in general, the convective velocity w takes into an account the solid deformation extending to the fluid part and ${\textbf{{w}}}={\textbf{{v}}}-\dot{\textbf{{u}}}$. In this paper, we restrict to stationary problems so that the fluid velocity is equivalent to the convective velocity,

[TABLE]

Moreover, by virtue of the linear kinematics (small displacements and deformations) only “one-way” interaction can be considered: the solid is loaded by traction forces due to the fluid stress, however the walls of the pores $\Omega_{f}$ are considered as rigid for the flow model. Therefore, the following interface conditions ensure the no-slip of the flow and the continuity of the normal stresses,

[TABLE]

Model restrictions

In this paper, we restrict the modelling to stationary problems, so that the time derivative $\partial_{t}c_{\alpha}$ in (2) vanishes. As the result, the fluid-structure mechanical interaction simplifies and also the electrochemical interactions in the fluid can be treated using the assumption of the electroneutrality in the bulk electrolyte. Finally, to define boundary conditions on the exterior boundary $\partial\Omega$. we shall assume that $(\Psi+\Psi^{\rm{ext}}),c_{i},{\textbf{{u}}},{\textbf{{w}}}$ and $p$ are $L-$periodic, where $L$ is the side of the cube $\Omega=]0,L[^{d}$.

3 Linearization and decomposition into subproblems

In order to apply the homogenization method, the complex nonlinear problem involving electrochemical and fluid-structure interactions is linearized. We introduce the equilibrium state which enables to establish a linearized problem for perturbation fields. The electrostatic problem is related to the heterogeneous periodic structure at the pore level.

3.1 Periodic structure of the porous medium

The porous medium is characterized by parameter $\ell^{\textrm{mic}}$ characterizing the size of the micropores. By virtue of the upscaling, the characteristic length $\ell^{\textrm{mic}}$ is related to a given macroscopic characteristic length $L$, such that the scale parameter introduced by $\varepsilon=\ell^{\textrm{mic}}/L$.

The porous medium is generated as a periodic lattice by repeating the representative volume element (RVE) occupying domain $Y^{\varepsilon}=\varepsilon Y$, see A. According to the decomposition introduced in Section 2, the zoomed cell $Y=\Pi_{i=1}^{3}]0,\bar{y}_{i}[\subset{\mathbb{R}}^{3}$ splits into the solid part occupying domain $Y_{s}$ and the complementary fluid part $Y_{f}$, thus

[TABLE]

Note that, roughly speaking, the relation between the micro- and macroscopic coordinates is $y=x/\varepsilon$, see the precise decomposition ansatz in (85). For a given scale $\varepsilon>0$, $\ell_{i}=\varepsilon\bar{y}_{i}$ is the characteristic size associated with the $i$-th coordinate direction, whereby also $\varepsilon\approx\ell_{i}/L$, hence $\ell_{i}\approx\ell^{\textrm{mic}}$ (for all $i=1,2,3$) specifies the microscopic characteristic length $\ell^{\textrm{mic}}$. Below we introduce two-scale functions depending on $x\in\Omega$ and $y\in Y$ using the unfolding operator

The ratio between the characteristic dimension of macroscopic domain $L_{c}$ (this may be associated with the specimen size $L$) and characteristic dimension of microscopic periodic structure $l$ is a scale parameter $\varepsilon=\frac{l}{L_{c}},0<\varepsilon\ll 1$ and represents the smallest zoom, by which the microstructure becomes visible from the macroscopic point of view.

3.2 Dimensionless problem

In this section we shall introduce a non-dimensional form of the equations from Section 2.2 and 2.3, following the approach from [4]. The macroscopic coordinate is rescaled by $x^{\prime}=x/L_{c}$ using the characteristic length $L_{c}$. Accordingly, the dimensionless operator $\nabla^{\prime}$ is defined by

[TABLE]

In the following text, however, we shall drop the prime ′, to simplify the notation, so that by $x$ we refer to the rescaled macroscopic coordinates. The dimensionless variables are expressed in terms of the characteristic quantities denoted by subscript $\sqcup_{c}$. As they are related to scale parameter $\varepsilon$, we denote them all by superscript $\sqcup^{\varepsilon}$.

[TABLE]

where $p_{c}$ is the characteristic pressure, $v_{c}$ the characteristic velocity, $c_{c}$ characteristic concentration, $\Psi_{c}=k_{B}T/e$ the characteristic potential, $u_{c}$ the characteristic displacement. Further, the dimensionless forcing terms labeled by $\sqcup^{*}$ are defined,

[TABLE]

The thickness of the electric double layer (EDL) is represented by the Debye length $\lambda_{D}$,

[TABLE]

where the last expression is due to the special type of the electrolyte with $z_{1}=-1,z_{2}=+1$. Apart of this parameter, other two dimensionless parameters are employed in the the dimensionless form of the system (2)-(8); these are the Peclect number $\textrm{Pe}_{\alpha}=\frac{l^{2}k_{B}Tc_{c}}{\eta_{f}D_{\alpha}}$, and the ratio between electrical and thermal energy $N_{\sigma}=\frac{el\Sigma_{c}}{\mathcal{E}k_{B}T}$.

Nondimesionalized problem

Using the dimensionless form of the system (2)-(8) we can state the following problem. Given ${\textbf{{f}}}^{*},{\textbf{{E}}}^{*},\Sigma^{*}$, find $c_{\alpha}^{\varepsilon},\Psi^{\varepsilon},p^{\varepsilon},{\textbf{{w}}}^{\varepsilon}$ and ${\textbf{{u}}}^{\varepsilon}$ which satisfy the following set of equations,

[TABLE]

with interface conditions

[TABLE]

where $\gamma=\frac{l^{2}}{\lambda_{D}^{2}\sum_{{\beta=1}}^{2}z_{\beta}^{2}}$ and

[TABLE]

Concentrations $c_{\beta}^{\varepsilon}$, sum of potentials $\left(\Psi^{\varepsilon}+\Psi^{\rm{ext},*}\right)$, convective velocity ${\textbf{{w}}}^{\varepsilon}$ and pressure $p^{\varepsilon}$ are $L$-periodic. We also recall that the given surface charge density $\Sigma^{*}$ is constant, as the consequence of the assumed constant potential $\Psi_{s}$ in the solid conductor. Introduction of the parameter $\varepsilon$ into (20),(22) and (27) is the natural consequence of the adimensional choices and dimensional analysis of the system, see B.

3.3 Linearization

To apply the homogenization efficiently, the system (20)-(28) must be linearized. For this, the assumption of sufficiently small applied fields $\Psi^{\rm{ext},*}$ and ${\textbf{{f}}}^{*}$ is needed. As the consequence, the state variables are only slightly perturbed from equilibrium, which justifies the linearization.

Following the linearization procedure in [3], any unknown $a^{\varepsilon}$ can be decomposed into its equilibrium part $a^{\rm{eq},\varepsilon}$ and its perturbation $\delta a^{\varepsilon}$, thus, we consider

[TABLE]

The equilibrium quantities labeled superscript $\sqcup^{\textrm{eq}}$, are solutions of the system (20)-(28) for ${\textbf{{f}}}{{}^{*}}=0$, $\Psi^{\rm{ext},*}=0$, ${\textbf{{w}}}^{\rm{eq},\varepsilon}=0$ and, by the consequence, zero diffusive fluxes $\mathbf{j}^{\varepsilon}_{\alpha}=0$; the last statement follows from (2), (4), and (29). Since the convective velocity vanishes at the equilibrium, we can state ${\textbf{{w}}}^{\varepsilon}(x)=\delta{\textbf{{w}}}^{\varepsilon}(x)$.

3.3.1 Equilibrium state quantities

Obviously, the equilibrium solution defines the reference state of the electrolyte such that the linearized state problem governs the perturbations. The existence of equilibrium solution $({c_{\alpha}^{\rm{eq},\varepsilon}(x)},{\Psi^{\rm{eq},\varepsilon}},{p^{\rm{eq},\varepsilon}},{\textbf{{u}}}^{\rm{eq},\varepsilon})$ was shown in [2], whereby the three fields ${c_{\alpha}^{\rm{eq},\varepsilon}(x)},{\Psi^{\rm{eq},\varepsilon}}$ and ${p^{\rm{eq},\varepsilon}}$ satisfy the following two relationships,

[TABLE]

where by ${c_{\alpha}^{b}}$ is the characteristic concentration in the bulk which represents the concentration of the ${\alpha}$-th ionic species is an infinite pore. As the consequence of (32), ${\Psi^{\rm{eq},\varepsilon}}$ determines ${c_{\alpha}^{\rm{eq},\varepsilon}(x)}$, then ${p^{\rm{eq},\varepsilon}}$.

To compute ${\Psi^{\rm{eq},\varepsilon}}$, the asymptotic analysis of the dimensionless Poisson-Boltzmann equation given by (22) and (27) has been treated in [2]. Therefore, for the sake of completeness, here we only provide the resulting expressions. Substituting equilibrium concentration (32)1 into the Poisson-Boltzmann equation (22) and (27), one gets

[TABLE]

where the $L$-periodicity of ${\Psi^{\rm{eq},\varepsilon}}$ is prescribed on the external boundary $\partial_{\rm{ext}}\Omega_{f}^{\varepsilon}$. The solvability of (33) for the zero Neumann condition, $\Sigma^{*}=0$, requires that the r.h.s. integrated in $\Omega_{f}^{\varepsilon}$ must vanish. This is satisfied provided the so-called electroneutrality condition in bulk holds,

[TABLE]

We adhere this condition, hence the existence of a unique solution ${\Psi^{\rm{eq},\varepsilon}}\in H^{1}_{\#}(\Omega^{\varepsilon}_{f})$ is guaranteed. From the physical point of view it ensures that ${\Psi^{\rm{eq},\varepsilon}}$ vanishes for the zero surface charge.

Although we assume $\Sigma^{*}$ to be a constant defined on interface $\Gamma^{\varepsilon}$, even for a periodic distribution of charges $\mathcal{T}_{\varepsilon}\left({\Sigma^{*}}\right)=\tilde{\Sigma}^{*}(y)$, $y\in\Gamma_{Y}$, the problem (33) yields $\varepsilon Y$-periodic solutions ${\Psi^{\rm{eq},\varepsilon}}$ in $\Omega_{f}^{\varepsilon}$, recalling the “macroscopic” $L$-periodicity on $\partial_{\rm{ext}}\Omega_{f}^{\varepsilon}$. This property allows us to consider only the local problem in the zoomed RVE represented by cell $Y_{f}$. Then,

[TABLE]

where concentrations ${c_{\alpha}^{\rm{eq}}}(y),{\alpha}=1,2$ obey the form of the Boltzmann distribution

[TABLE]

Potential ${\Psi^{\rm{eq}}}(y)\in H^{1}_{\#}(Y_{f})$ is a solution of the Poisson-Boltzmann equation (33) imposed in $Y_{f}$, in particular

[TABLE]

To conclude, by virtue of (32), the unfolded equilibrium concentrations $\mathcal{T}_{\varepsilon}\left({{c_{\alpha}^{\rm{eq},\varepsilon}}}\right)$ and the unfolded pressure field $\mathcal{T}_{\varepsilon}\left({{p^{\rm{eq},\varepsilon}}}\right)$ are $Y$-periodic functions. Moreover, by the consequence, the unfolded displacements are also $Y$-periodic functions, whereby the macroscopic strains vanish. Therefore, in this paper, we neglect any influence of the equilibrium displacements field on the reference configuration associated with the linearization procedure considered in what follows. Note that, the equilibrium pore geometry might be perturbed due to the local strains in $Y_{s}$.

3.3.2 Perturbed state quantities

As the further step in the linearization, the total electrostatic potential $\Psi_{f}^{\varepsilon}$ is decomposed according to phenomena which participate in the total electric field. This can be considered as a superposition $\Psi_{f}^{\varepsilon}=\Psi^{\varepsilon}+{\Phi_{\alpha}^{\varepsilon}}+\Psi^{\rm{ext}}$ of local particular electric fields associated with potentials $\Psi^{\varepsilon}$, ${\Phi_{\alpha}^{\varepsilon}}$, and $\Psi^{\rm{ext}}$

potential $\Psi^{\rm{ext}}$ which yields the external electrical field ${\textbf{{E}}}^{*}=-\nabla\Psi^{\rm{ext}}$ is imposed and independent of $\varepsilon$;

• 2.

potential $\Psi^{\varepsilon}$ reflects only the effects of the EDL on the ion distribution; in the equilibrium state, $\Psi^{\varepsilon}$ is given by problem (33);

• 3.

ionic potentials $\Phi_{\alpha}$ (often referred to as the streaming potential) represents the electric field produced by motion of ${\alpha}-$th ionic species. In the equilibrium, ${\Phi_{\alpha}^{\varepsilon}}$ vanishes since both the convection ${\textbf{{w}}}^{\varepsilon}$ and the ionic flux $\mathbf{j}^{\varepsilon}_{\alpha}$ vanish. Thus, the ionic potential is identified by its perturbation only, ${\Phi_{\alpha}^{\varepsilon}}(x)={\delta\Phi_{\alpha}^{\varepsilon}}$.

To summarize the decomposition, the potential in equilibrium is given only by the electrokinetic potential, so that $\Psi^{\rm{tot,eq}}={\Psi^{\rm{eq}}}$. The decomposition of the potential in equilibrium and under the flow is also illustrated in Fig.2.

Ionic concentrations can be expressed in the context of Boltzmann distribution, so that

[TABLE]

The introduction of ionic potentials will help to eliminate the boundary condition (27) from the system, because the mobility of particles is not influenced by the choice of boundary conditions, see [22]. The boundary condition (33)2 is considered only to define $\Psi^{\rm{eq},\varepsilon}$, the potential distribution in equilibrium. It will be shown later that potential ${\Phi_{\alpha}^{\varepsilon}}$ proves useful in the decoupling the electrokinetic system.

The linearization of (38) by the first-order Taylor expansion yields

[TABLE]

Further, following the work [20], it is convenient to introduce the so-called global pressure $P^{\varepsilon}$,

[TABLE]

which consists of the hydrodynamic pressure perturbation $\delta p^{\varepsilon}$ and the osmotic pressure, see [15].

Finally, the decomposition of unknowns fields (31) is substituted into the dimensionless problem (20)-(28) and (39) is employed to express $\delta c_{\alpha}^{\varepsilon}(x)$. Note, that products of the small quantities such as $\delta c_{\alpha}^{\varepsilon}\delta\Psi^{\varepsilon}$ or $(\delta c_{\alpha}^{\varepsilon})^{2}$ are neglected. Due to the linearization and the use of the global pressure (40), the nondimensionalized problem splits into three subproblems which can be solved subsequently.

Linearized electrokinetic system

Given the body forces ${\textbf{{f}}}^{*}$ and potential $\Psi^{\rm{ext},*}$, find $L$-periodic functions $(\delta{\textbf{{w}}}^{\varepsilon},P^{\varepsilon},{\delta\Phi_{\alpha}^{\varepsilon}},\delta\Psi^{\varepsilon})$ and $\delta{\textbf{{u}}}^{\varepsilon}$ which solve the following three subproblems:

Electrokinetic problem: $(\delta{\textbf{{w}}}^{\varepsilon},P^{\varepsilon},{\delta\Phi_{\alpha}^{\varepsilon}})$ satisfy

[TABLE] 2. 2.

Electrostatic EDL problem: $\delta\Psi^{\varepsilon}$ satisfies

[TABLE] 3. 3.

Deformation problem: $\delta{\textbf{{u}}}^{\varepsilon}$ satisfies

[TABLE]

with the linearized fluid stress given by

[TABLE]

It is worth to note that, due to the linearization and decoupling in three subproblems, (44)-(45), (47) and (49) present stanard boundary conditions on interface $\Gamma^{\varepsilon}$, rather than transmission conditions, as in problem (20)-(28).

4 Homogenization

The unfolding homogenization method [9] has been used for the asymptotic analysis $\varepsilon\rightarrow 0$ of weak formulation (1) arising from the linear system (41)-(49). Since the resulting system of the limit two-scale equations corresponds to the one obtained in paper [3], in the next sections, we report only briefly the on the upscaling procedure. The main purpose is to explain the structure of the homogenized model, namely to define local problems for the so-called characteristic responses, to give formulae for computing the homogenized coefficients, and to formulate the macroscopic problem.

4.1 Convergence results

The convergence analysis is derived for the weak formulation of the linearized problem (41)-(49). The following functional spaces will be employed.

[TABLE]

where $H^{1}(\Omega_{p}^{\varepsilon})$ is the Sobolev space $W^{1,2}(\Omega_{p}^{\varepsilon})$, subscript $p=s,f$.

Weak formulation of the linearized electrokinetic system

Given ${\textbf{{f}}}^{*}\in{\textbf{{L}}}^{2}(\Omega^{\varepsilon}_{f})$ and ${\textbf{{E}}}^{*}\in{\mathbb{R}}^{3}$, whereby $\Psi^{\rm{ext},*}=x\cdot{\textbf{{E}}}^{*}$.

Find $(\delta{\textbf{{w}}}^{\varepsilon},P^{\varepsilon},{\delta\Phi_{\alpha}^{\varepsilon}})\in{\mathbf{H}}^{1}_{\#0}(\Omega^{\varepsilon}_{f})\times[H^{1}_{\#}(\Omega^{\varepsilon}_{f})]^{3}$, such that

[TABLE] 2. 2.

Find $\delta\Psi^{\varepsilon}\in H^{1}_{\#}(\Omega^{\varepsilon}_{f})$, such that

[TABLE] 3. 3.

Find ${\textbf{{u}}}^{\varepsilon}\in H^{1}_{\#}(\Omega_{s})^{d}$, such that

[TABLE]

where $\boldsymbol{\sigma}^{\varepsilon}_{f}$ is given by (50).

The two-scale limit problem can be obtained due to the weak convergences in the unfolded domain $\Omega\times Y$; the two-scale convergence method was used in [1] and [4]. The unfolded equations of the weak formulation are obtained using the unfolding operator $\mathcal{T}_{\varepsilon}$ defined in A, see [9]. Due to the a priori estimates on the solutions of (1)-(55), according to [4] the following convergence result for $\varepsilon\rightarrow 0$ can be proved: There exist limit fields $(\delta{\textbf{{w}}}^{0},P^{0})\in L^{2}(\Omega;H^{1}_{\#}(Y_{f})^{d})\times L^{2}(\Omega)$, $\left\{{\delta\Phi}^{0}_{\beta},{\delta\Phi}^{1}_{\beta}\right\}_{{\beta=1},\dots,N}\in(H^{1}(\Omega)\times L^{2}(\Omega;H^{1}_{\#}(Y_{f})))^{N}$, ${\delta\Psi}^{0}\in L^{2}(\Omega;H^{1}_{\#}(Y_{f}))$, $\delta{\textbf{{u}}}^{0}\in H^{1}_{\#}(\Omega_{s}^{\varepsilon})^{d}$ and $\delta{\textbf{{u}}}^{1}\in L^{2}(\Omega;H^{1}_{\#}(\Omega_{s}^{\varepsilon})^{d})$ such that following convergences hold

[TABLE]

for ${\beta=1,2}$ and where

As the consequence of the convergences (56), truncated asymptotic expansions of the unfolded unknown fields can be introduced which satisfy the same convergence result. These constitute the recovery sequences in subdomains of $\Omega_{f}^{\varepsilon}$ and $\Omega_{s}^{\varepsilon}$ (the respective characteristic functions are $\chi_{f}^{\varepsilon}$ and $\chi_{s}^{\varepsilon}$) which read

[TABLE]

Limit functions, $(\delta{\textbf{{w}}}^{0},P^{0},\{{\delta\Phi}^{0}_{\beta},{\delta\Phi}^{1}_{\beta}\},\delta\Psi^{0}(x,y),\delta{\textbf{{u}}}^{0},\delta{\textbf{{u}}}^{1})$ are solutions of the corresponding two-scale limit problems. These are not presented in this paper; we only introduce the multiplicative splits of the two-scale functions which allow us to establish local autonomous problems for characteristic responses.

4.2 Scale separation formulas

The local problems, relevant to the microscopic scale, can be derived from the limit problem, usually by letting all the components of the test functions, which are not relevant to the microscopic scale. Thanks to the linearization of the problem, it is possible to determine scale separation formulae as a linear combination of macroscopic fluxes and corrector base functions, otherwise called the characteristic responses. Two different macroscopic fluxes in the limit two-scale problem are recognized, namely $(\nabla_{x}{\delta\Phi_{\beta}}^{0}+{\textbf{{E}}}^{*})$ and $(\nabla_{x}P^{0}+{\textbf{{f}}}^{*})$. Therefore, we introduce scale decomposition formulae of the limits $\delta{\textbf{{w}}}^{0},P^{1},{\delta\Phi_{\beta}}^{1},{\beta=1,2}$, which read

[TABLE]

where two families of corrector base functions $(\boldsymbol{\omega}^{0,k},\pi^{0,k},\theta_{\beta}^{0,k})$ and $(\boldsymbol{\omega}^{{\alpha},k},\pi^{{\alpha},k},\theta_{\beta}^{{\alpha},k}),{\beta=1,2},$ were introduced, indexed by $k\in\{1,\dots,d\}$; note $d=3$ is the spatial dimension of the problem.

The scale separation formulae for the potential and the displacements attain the following forms:

[TABLE]

where ${\textbf{{w}}}^{ij},{\textbf{{w}}}^{P}$ and ${\textbf{{w}}}^{{\beta}}$ are corrector base functions.

All corrector functions involved in (58) and (59) are obtained as the solution of three local cell problems given in the next section.

4.3 Cell problems

Corrector basis functions introduced in (58)-(59) satisfy the local autonomous problems (cell problems) defined in subdomains of the representative periodic cell $Y$, therefore. By virtue of the linearization which decomposes the problem (41)-(49) in three subproblems, the cell problems are decoupled, thus, the following three groups are distinguished:

Group 1: two cell problems related to the electrokinetic system with responses $(\boldsymbol{\omega}^{0,k},\pi^{0,k},\theta_{\beta}^{0,k})$ and $(\boldsymbol{\omega}^{\alpha,k},\pi^{\alpha,k},\theta_{\beta}^{\alpha,k}),\alpha,\beta=1,2$, $k=1,\dots,d$.

• 2.

Group 2: one cell problem related to the electrostatic potential in the EDL with solution $\varpi^{\alpha},{\alpha=1,2}$.

• 3.

Group 3: three cell problems related to the poroelasticity with solutions $,{\textbf{{w}}}^{ij},{\textbf{{w}}}^{P}$ and ${\textbf{{w}}}^{i}$

While denoting $\mathbf{e}^{k}$ canonical basis of ${\mathbb{R}}^{d}$, $k=1,...,d$ and $\delta_{ij}$ the Kronecker symbol, all the six cell problems are introduced below.

4.3.1 Cell problems: Group 1

The first autonomous cell problem of this group is related to the macroscopic pressure gradient: Find $(\boldsymbol{\omega}^{0,k},\pi^{0,k},\theta_{\beta}^{0,k})\in{\mathbf{H}}^{1}_{\#0}(Y_{f})$, $k=1,2,3,\ {\beta=1,2}$:

[TABLE]

for all test functions $\tilde{\boldsymbol{\vartheta}}\in{\mathbf{H}}^{1}_{\#0}(Y_{f})$, $\tilde{q}\in L^{2}(Y_{f})$, $\tilde{\varphi}\in H^{1}_{\#}(Y_{f})^{N}$.

The second autonomous cell problem, corresponding to the macroscopic diffusive flux, for each species ${\alpha=1,2}$ reads: Find $(\boldsymbol{\omega}^{{\alpha},k},\pi^{{\alpha},k},\theta_{\beta}^{{\alpha},k})\in{\mathbf{H}}^{1}_{\#0}(Y_{f})$, $k=1,2,3,\ {\beta=1,2}$:

[TABLE]

for all test functions $\tilde{\boldsymbol{\vartheta}}\in{\mathbf{H}}^{1}_{\#0}(Y_{f}),\tilde{q}\in L^{2}(Y_{f}),\tilde{\varphi}\in H^{1}_{\#}(Y_{f})^{N}$.

4.3.2 Cell problems: Group 2

The cell problem associated with the macroscopic ionic potential, for each species ${\alpha=1,2}$ reads: Find corrector base functions $\varpi^{\alpha}\in H^{1}_{\#}(Y_{f}),{\alpha=1,2}$, such that

[TABLE]

for all test functions $\tilde{\varphi}\in H^{1}_{\#}(Y_{f})$.

4.3.3 Cell problems: Group 3

The following three cell problems are relevant to the homogenization of displacement perturbation. One can realize, that the first two cell problems are identical to the ones occurring in derivation of Biot’s poroelasticity equation. The first cell problem reads: find ${\textbf{{w}}}^{ij}\in{\bf{H}}_{\#}^{1}(Y_{s})^{d},\sim\kern-11.00008pt\int\limits_{Y_{s}}{\textbf{{w}}}^{ij}\mathrm{\,dV}=0$ such that

[TABLE]

for any test function $\tilde{{\textbf{{v}}}}\in{\bf{H}}_{\#}^{1}(Y_{s})^{d}$. Symbol $\boldsymbol{\Pi}^{ij}$ denotes the so-called transformation vectors $\boldsymbol{\Pi}^{ij}=(\Pi^{ij}_{k}),i,j,k=1,\dots,d$, which enable to establish local displacements defined in $Y$ generated by affine transformation of the macroscopic strains ${\textbf{{e}}}_{x}(\delta{\textbf{{u}}}^{0})$ defined in $\Omega$; it holds that ${\textbf{{e}}}_{y}(\boldsymbol{\Pi}^{ij}e_{ij}^{x}(\delta{\textbf{{u}}}^{0})={\textbf{{e}}}_{x}(\delta{\textbf{{u}}}^{0})$, where

[TABLE]

The second cell problem reads: find ${\textbf{{w}}}\in{\bf{H}}_{\#}^{1}(Y_{s}),\sim\kern-11.00008pt\int\limits_{Y_{s}}{\textbf{{w}}}^{P}\mathrm{\,dV}=0$ such that

[TABLE]

for any test function $\tilde{{\textbf{{v}}}}\in{\bf{H}}_{\#}^{1}(Y_{s})^{d}$.

Finally, the third cell problem connecting displacement perturbation and ionic potentials is needed. It reads: Find ${\textbf{{w}}}^{i}\in{\bf{H}}_{\#}^{1}(Y_{s})^{d},\sim\kern-11.00008pt\int\limits_{Y_{s}}{\textbf{{w}}}^{\alpha}\mathrm{\,dV}=0$ such that

[TABLE]

for any test function $\boldsymbol{\varphi}\in{\bf{H}}_{\#}^{1}(Y_{s})^{d}$.

4.4 Macroscopic model

By virtue of the homogenization method, the limit two-scale equations arising from (1)-(55) involve cell integrals of the two-scale functions which can be expressed in terms of the corrector basis functions. Below we list expressions of the homogenized coefficients which constitute the effective material parameters of the upscaled porous medium, [2].

The first group of the corrector functions define the following homogenized coefficients

[TABLE]

whereby $\boldsymbol{\mathcal{K}}=(\mathcal{K}_{lk})$ is the permeability tensor, $\boldsymbol{\mathcal{D}}^{{\alpha}{\beta}}=(\mathcal{D}^{{\alpha}{\beta}}_{lk})$ are diffusivity tensors; in particular $DD^{{\alpha}{\beta}}$ describes diffusion of species $\alpha$ due to the streaming potential gradient of the species $\beta$. Tensors $\boldsymbol{\mathcal{J}}^{\alpha}=(\mathcal{J}^{{\alpha}}_{lk})$ is related to the flow driven by electric fields and $\boldsymbol{\mathcal{L}}^{\alpha}=(\mathcal{L}^{{\alpha}}_{lk})$ , also known as the coupling tensor, expresses the diffusivity of species $\alpha$ due to the global pressure gradient.

The second and the third group of the corrector functions constitute poroelastic coefficients modified by the presence of the streaming potentials and the external electric field,

[TABLE]

Above the tensor $\boldsymbol{\mathcal{A}}=(\mathcal{A}_{ijkl})$ is the fourth-order positive definite effective elasticity tensor of drained skeleton, $\boldsymbol{\mathcal{B}}^{H}=(\mathcal{B}_{ij})$ is the Biot’s coupling tensor related to the pressure, while $\boldsymbol{\mathcal{C}}^{\alpha}=(\mathcal{C}^{\alpha}_{kl})$ is the tensor related to ionic potentials. These effective coefficients are sometimes referred to as the Biot poroelasticity coefficients. For convenience we may introduce coefficient

[TABLE]

Macroscopic model

We present the macroscopic model in its dimensional form, *i.e. *using quantities with physical dimensions. The macroscopic variables obtained as the limits of the oscillating solutions of the nondimensionalized problem (1)-(55) can be presented by their dimensional macroscopic counterparts – these variables will be denoted by superscript $\sqcup^{\mathrm{eff}}$. Dimensional form of the macroscopic problem reads: Find $({{P^{\rm{eff}}},{\Phi_{\alpha}^{\rm{eff}}},{{\textbf{{u}}}^{\rm{eff}}}})\in(L^{2}(\Omega)\times H^{1}(\Omega)^{d})$, such that

[TABLE]

for all test functions $q\in L^{2}(\Omega),\varphi\in L^{2}(\Omega)$ and ${\textbf{{v}}}\in H^{1}(\Omega)^{d}$ and where dimensionalized permeability tensor is denoted by $\boldsymbol{\mathcal{K}}^{D}=l^{2}\boldsymbol{\mathcal{K}}$ and $\boldsymbol{\mathcal{A}}^{G}=\Lambda\boldsymbol{\mathcal{A}}$ is Gassmann elasticity tensor. The equation (77) is so-called extended Biot equation and it is weakly coupled to the electrokinetic system (4.4) and (76) through coefficients $\boldsymbol{\mathcal{C}}^{\alpha}$.

As reported in the next section, the weak formulation (4.4)-(77) is discretized using the finite elements to obtain the numerical solutions. For the sake of completeness, we also introduce the macroscopic homogenized model in its differential form:

[TABLE]

for ${\alpha=1,2}$ and completed by periodic boundary conditions. From the (79) and (80) we can distinguish the fluid seepage and the ionic diffusion fluxes,

[TABLE]

5 Numerical simulation

The aim of this section is to explore and illustrate properties of the homogenized two-scale model described in preceding sections. For this purpose, we present numerical simulation of the ionic transport through a porous medium occupying a simple-shaped macroscopic domain, with a simple periodic microstructure. We perform a parametric study, which illustrates influence of a varying microstructure on the homogenized material properties.

The two-scale homogenized model was implemented in SfePy, a software for solving problems with coupled partial differential equations (PDEs) in weak forms by means of the finite element method (FEM) for 2D and 3D problems, [8]. SfePy is based on the Python programming language and its packages NumPy and SciPy, [12].

5.1 Algorithm of numerical implementation

In the presented examples do not consider volume forces and also disregard effects of an external electric field, thus we put ${\textbf{{f}}}=\mathbf{0}$ and ${\textbf{{E}}}=\mathbf{0}$. The used quantities and parameters are in Tab. 1. All the computations are performed for the given pore size $l$.

The numerical simulation of the problem can be divided into several steps:

Solve the potential distribution in equilibrium ${\Psi^{\rm{eq}}}$ on cell $Y$ as a solution of (37). 2. 2.

Compute concentrations ${c_{\beta}^{\rm{eq}}},{\beta=1,2}$ from (36). 3. 3.

Compute corrector functions $(\boldsymbol{\omega}^{0,k},\pi^{0,k},\theta_{\beta}^{0,k})$ and $(\boldsymbol{\omega}^{{\alpha},k},\pi^{{\alpha},k},\theta_{\beta}^{{\alpha},k})$, $k=1,\dots,d$, ${\alpha=1,2}$, ${\beta=1,2}$, related to electrokinetic system as a solution of local problems (60) and (61). 4. 4.

Compute effective coefficients relevant to the decoupled electrokinetic system from (67)-(70). 5. 5.

Compute corrector functions $\varpi^{{\alpha}},{\alpha=1,2}$ related to the potential perturbation as a solution of local problems (62). 6. 6.

Compute corrector functions $({\textbf{{w}}}^{ij},{\textbf{{w}}}^{P},{\textbf{{w}}}^{\alpha})$, ${\alpha=1,2}$, $i,j=1,\dots,d$, related to the displacement perturbation as a solution of local problems (63)-(4.3.3). 7. 7.

Compute effective coefficients relevant to the Biot poroelasticity from (71)-(73). 8. 8.

Compute solution to the macroscopic homogenized system of equations (4.4)-(77).

Since the major part of the equations were presented in their dimensionless form, all the homogenization results will be presented in the dimensionless form as well, unless stated otherwise.

5.2 Geometrical representation of microstructure

We aim to study the dependency of the effective coefficients on a change of the microstructure. For this purpose we choose only a simple geometry representation in the form of three interconnected canals aligned with $y_{1}$-, $y_{2}$- and $y_{3}$-directions trough a continuous matrix, see Fig. 3(a). The cross-section of the canals is a square with size $a$. Changing the parameter $a$ leads to a change in porosity $\hat{\phi}_{f}=\frac{|Y_{f}|}{|Y|}$.

The mesh representing the cell $Y$ was generated by a mesh generation script, which forms a part of the SfePy software. For meshing purposes, the linear hexahedron elements were used.

5.3 Semi-discretized macroscopic problem

We introduce the semi-discretized form of the macroscopic problem, which can be used in the finite element (FE) model. By $\mathbf{P}^{\rm{eff}},\boldsymbol{\Phi}^{\rm{eff}}_{\alpha}{\alpha=1,2}$ and ${\textbf{{u}}}^{\rm{eff}}$ we refer to the column vectors incorporating all degrees of freedom of FE mesh nodes associated with partitioning of the macroscopic domain $\Omega$. We need the following approximations of the terms involved in the macroscopic problem (4.4) - (77):

[TABLE]

Using the notations just introduced, we can write the linear macroscopic problem in the matrix form

[TABLE]

From the macroscopic problem (78)-(80) is immediately evident, that the electrokinetic system can be solved separately from the problem of poroelasticity, thus the solution $(\mathbf{P}^{\rm{eff}},\boldsymbol{\Phi}^{\rm{eff}}_{\alpha}),{\alpha=1,2}$ is obtained. Than, the macroscopic displacement can be found as

[TABLE]

5.3.1 Potential distribution in equilibrium

The first step in obtaining effective coefficients is to compute the distribution of potential ${\Psi^{\rm{eq}}}$ in equilibrium on the microscopic scale. All used electrokinetic quantities can be found in Tab. 1. The distribution of the dimensionless equilibrium potential ${\Psi^{\rm{eq}}}$ is shown in Fig. 3(b). The potential ${\Psi^{\rm{eq}}}$ has its maximum on the solid-fluid interface, where the surface charge $\Sigma^{*}$ is prescribed. The potential ${\Psi^{\rm{eq}}}$ gradually decreases with increasing distance from the interface. This meets our general expectation about Poisson-Boltzmann potential distribution near the solid-fluid interface. The resulting potential distribution ${\Psi^{\rm{eq}}}$ is needed for computation of concentration ${c_{\beta}^{\rm{eq}}}$ and for subsequential calculations.

5.3.2 Influence of varying microstructure on effective tensors

We study the dependency of effective coefficients on the change in porosity $\hat{\phi}_{f}$, caused by variation in $y_{1}$-direction canal size parameter $a$. These dependencies of the diagonal components of electrokinetic and poroelasticity coefficients (dimensionless) can be seen in Fig. 4 and Fig. 5. In accordance with the choice of the symmetric microstructure, the components of effective coefficients related to $y_{2}$- and $y_{3}$-direction are expected to be equal.

Fig. 4 shows the dependency of the effective coefficient related to poroelasticity on the porosity of the microstructure. The upper left graph shows components of dimensionless poroelasticity tensor $\boldsymbol{\mathcal{A}}$ and its dependency on the microstructure porosity. As expected, all the components of the poroelasticity decrease with the increasing porosity. In the upper right graph, the Biot coefficient $\boldsymbol{\mathcal{B}}$ increases with the porosity. Finally, the lower half of the figure shows ionic potential tensors $\boldsymbol{\mathcal{C}}^{1}$ and $\boldsymbol{\mathcal{C}}^{2}$. The components of $\boldsymbol{\mathcal{C}}^{1}$ related to all three direction are equal. This property applies to the components of $\boldsymbol{\mathcal{C}}^{2}$ as well. However, the tensor $\boldsymbol{\mathcal{C}}^{1}$ decrease and the tensor $\boldsymbol{\mathcal{C}}^{2}$ increase nonlinearly with the increasing porosity.

Similar nonlinear behavior is obtained for the other electrokinetic tensors as seen in Fig. 5. We observe the decrease in the components of the migration-diffusion tensors $\boldsymbol{\mathcal{J}}^{1}$ and $\boldsymbol{\mathcal{J}}^{2}$ related to anions and cations, respectively. The components of coupling tensors $\boldsymbol{\mathcal{L}}^{1}$ and $\boldsymbol{\mathcal{L}}^{2}$ decrease with the increasing porosity.

In the lower part of Fig. 5 we observe increase in permeability $\boldsymbol{\mathcal{K}}$, as expected. The last part of this figure depicts the decreasing components of diffusivity tensors $\boldsymbol{\mathcal{D}}^{11}$ and $\boldsymbol{\mathcal{D}}^{22}$. The diagonal components of $\boldsymbol{\mathcal{D}}^{12}$ and $\boldsymbol{\mathcal{D}}^{21}$ are identic with those of $\boldsymbol{\mathcal{D}}^{22}$.

5.3.3 Solution of macroscopic problem

For the purpose of numerical simulations we propose a simple experiment, see Fig. 6, where a small cuboid specimen occupied by a porous medium is placed between two ionic reservoirs and separated by semipermeable membranes. These membranes enable ionic exchange, but prevent fluid flow.

To describe the boundary conditions, we refer to the faces of the porous specimen by the intuitive notation, such that $\Gamma_{E}$ stands for the “east side” with the normal vector aligned with $x_{1}$-axis, whereas the “north side” $\Gamma_{N}$ has its normal aligned with $x_{2}$-axis. Then $\Gamma_{T}$ and $\Gamma_{B}$ refer to the top and bottom sides, respectively, see Fig. 6. On the “top” and “bottom” boundaries $\Gamma_{UT}$ and $\Gamma_{UB}$, respectively, the porous specimen is clamped.

This experiment is focused on the observation of the displacement and pressure distribution under the ionic potentials change. To this aim, we propose four macroscopic problems with varying boundary conditions related to the ionic potentials. The homogenized macroscopic problem is given by the system of equations (4.4)-(77) completed by its respective boundary conditions. In what follows we define the two sets of boundary conditions, thus obtaining two boundary problems. The two boundary value problems (BVP) are defined in the following part.

In order to prevent numerical errors, the problem is computed in its dimensionless form. Then, using the dimensionless choices from Section 3.2, we recover dimensional form of respective macroscopic quantities.

5.3.4 Boundary value problems

The boundary value problem is defined by (4.4)-(77) and by the following boundary conditions:

on $\Gamma_{UT}$ and $\Gamma_{UB}$: ${{\textbf{{u}}}^{\rm{eff}}}=\mathbf{0}$

• 2.

on $\Gamma_{W}$ and $\Gamma_{E}$: ${P^{\rm{eff}}}=0$

• 3.

on $\Gamma_{T}$: ${\Phi_{1}^{\rm{eff}}}=\bar{\Phi}_{1},\quad\boldsymbol{\mathcal{J}}^{2}\nabla{\Phi_{2}^{\rm{eff}}}\cdot{\textbf{{n}}}=\bar{g}$

• 4.

on $\Gamma_{B}$: ${\Phi_{1}^{\rm{eff}}}=\bar{\Phi}_{1},\quad\boldsymbol{\mathcal{J}}^{2}\nabla{\Phi_{2}^{\rm{eff}}}\cdot{\textbf{{n}}}=b\bar{g}$

The values of boundary conditions are $\bar{\Phi}_{1}=0.1$ and $\bar{g}=0.001$. By the choice of parameter $b$ we distinguish two BVPs.

The first boundary problem (BVP I) is defined by the choice $b=1$, so that the boundary conditions are symmetric on boundaries $\Gamma_{T}$ and $\Gamma_{B}$. Therefore, the symmetric distribution of macroscopic quantities is obtained correspondingly, as seen in Fig. 7. The deformed shape is visualized by the wire-frame, whereby the displacement field is enlarged by factor $2\cdot 10^{5}$. The swelling of the macroscopic body occurs mainly in the region, where ${\Phi_{2}^{\rm{eff}}}$ attains the lowest values.

By taking $b=5$ we get the second boundary value problem (BVP II). In this case, we increased the influx of ${\Phi_{2}^{\rm{eff}}}$ on the boundary $\Gamma_{B}$. This leads to the contraction of the porous specimen near this surface, as seen in Fig. 8. Naturally, the non-symmetrical boundary conditions result in a non-symmetric distribution of the macroscopic quantities. This effect is most visible on the swelling of the macroscopic body, which tends to react to the distribution of ${\Phi_{2}^{\rm{eff}}}$, being slightly more swelled where ${\Phi_{2}^{\rm{eff}}}$ is slightly lower. The deformed shape is visualized by the wire-frame enlarged by factor $2\cdot 10^{4}$ only, in this case.

5.4 Reconstruction of macroscopic solution on microstructure level

One of the most remarkable advantages of the homogenization method is the possibility to reconstruct the solution at the microscopic scale. After computing the global (dimensionless) responses $\left\{{\textbf{{u}}}^{0},P^{0},{\Phi^{0}}_{\beta}\right\}$, it is possible to reconstruct the associated microscopic quantities. This process is also called downscaling in contrast to the upscaling process leading to the macroscopic model. Let us now briefly introduce the reconstruction relations.

We consider a given $\varepsilon>0$ corresponding to a real size of the microstructure. This enables us to apply the decomposed forms of the fluctuating two-scale functions in (58)-(59) ${\textbf{{u}}}^{0},P^{0},{\Phi^{0}_{\alpha}},{\alpha=1,2}$ defined in domain $\Omega$ with the local characteristic responses (corrector basis functions) defined in $Y$. The local microscopic fields are given by the so-called folding mapping ([28]), such that

[TABLE]

The folding operator $\mathcal{F}^{\varepsilon}$ combines corrector basis functions defined in $Y$ with interpolated macroscopic responses transformed to the zoomed cell $Y$ by the operator $\mathcal{Q}^{\varepsilon}$. The operator $\mathcal{Q}^{\varepsilon}$ is average operator over the recovery cell.

Using this approach, the microscopic fields ${\textbf{{u}}}^{{\textrm{mic}},\varepsilon},{\textbf{{w}}}^{{\textrm{mic}},\varepsilon},P^{{\textrm{mic}},\varepsilon},\Phi_{\beta}^{{\textrm{mic}},\varepsilon},{\beta=1,2}$ can be reconstructed as follows:

[TABLE]

For the purpose of numerical modeling, we recover the dimensionless macroscopic quantities on the microscopic cell. The macroscopic quantities were taken from the solution of BVP I, the reconstruction was made once again in SfePy. As the recovery cell was taken the cell in the center of the porous specimen.

As it can be seen in Fig. 9 in 2D a 3D views (depending on the specimen orientation w.r.t. the observer), the swelling occurs inward to the canal on the microscopic level. The reconstructed velocity field ${\textbf{{w}}}^{{\textrm{mic}},\varepsilon}$ shows the electrolyte passing mainly through $y_{1}$- and $y_{3}$-directions of the canal. The reconstructed pressure $P^{{\textrm{mic}},\varepsilon}$ is shown in the top half of Fig. 10, where 3D view is illustrated using slices through all three axes of the canal. The pressure is lowest in the canal center, while highest near the solid-fluid interface. This could be traced down to the response to the solid part swelling, but is somewhat unsymmetrical due to the connections to the other macroscopic quantities. Finally, the visualization of the reconstructed potential ${\delta\Psi}^{{\textrm{mic}},\varepsilon}$ can be found in the lower half of Fig. 10. As was mentioned earlier, the ionic potentials are locally constant. And truly, reconstructed $\Phi_{\beta}^{{\textrm{mic}},\varepsilon}$ are constant on the whole fluid part of the microscopic cell $Y_{f}$ with values $\Phi_{1}^{{\textrm{mic}},\varepsilon}=5.003\times 10^{-13}$ and $\Phi_{2}^{{\textrm{mic}},\varepsilon}=1.543\times 10^{-12}$.

6 Conclusion

We applied the two-scale homogenization of the deformable porous medium saturated by two-component electrolyte and implemented the resulting two-scale model in our open source finite element code SfePy. The upscaled model for this type of media was derived in [4], however, without any computational analysis, or illustrative examples. Therefore, up to our knowledge, the computational study presented in our paper provides first quantitative analysis of the considered medium. Its behavior is illustrated in terms of the example which mimics an experiment. In the computational study, dependence of the homogenized coefficients on the fluid volume fraction was presented, whereby a symmetric geometry of the reference periodic cell generating the porous medium was employed; influences of anisotropy and other geometry-related features will be studied in our further research.

The model describes the steady state of the electrolyte flow in the solid skeleton made of an elastic electric conductor. The homogenization is applied to the linearized model is obtained for the equilibrium reference state which is defined under the assumptions of zero fluxes (fluid and solid velocities and electroneutrality in bulk of the electrolyte). Moreover, the surface electric charge is a given constant on the solid-fluid interface.

It is important to note that the homogenized model, thereby values of the homogenized coefficients are specific to a given microstructure size, as the result of the scale-dependent effects associated with the viscous flow and the electric double layer. Thus, a given characteristic microstructure size $\ell$ determines the Peclet number and another coefficient influencing the equilibrium electric potential.

By virtue of the “downscaling” procedure of the homogenization, the macroscopic fields can be used to reconstruct responses at the microscopic level. In particular, the corrector functions combined with local values of the macroscopic fields provide two-scale functions representing the $\varepsilon$-order fluctuating parts of the global pressure, the streaming potentials of the species, and of the displacements. In analogy, the reconstruction provides also the fluid velocity and the electric potential associated with the double layer potential which both are purely two-scale functions relevant to the microscopic cell. The perturbations of ionic concentrations from the electroneutrality state are recovered using the perturbed electric potentials.

There are further extensions of the present two-scale model and its computational implementation. To treat non-steady flows, the fluid-structure interaction will be more involved, thus, leading to a strong coupling between the flow, ionic concentrations and deformation. The scale decoupling procedure will be more complicated and will lead to fading memory effects of the macroscopic responses, as the homogenized coefficients will serve for time convolution kernels, cf. [7, 25].

Acknowledgment

This research is supported by part by project GACR 16-03823S of the Scientific Foundation of the Czech Republic and by the project LO 1506 of the Czech Ministry of Education, Youth and Sports. The work was also supported from European Regional Development Fund-Project „Application of Modern Technologies in Medicine and Industry” (No. CZ.02.1.01/0.0/0.0/17_048/0007280).

Appendix A Unfolding homogenization

The unfolding homogenization method is based on the properties of unfolding operator $\mathcal{T}_{\varepsilon}$ which is similar to the dilatation operator. By virtue of the coordinate decomposition into ”coarse” and ”fine” parts, any function $\psi=\psi(x)$ can be unfolded into a function of $x$ and $y$. The convergence results in the unfolded domains $\Omega\times Y$ can be found in [9]. By virtue of its definition, for specific subsequences of $\varepsilon$, domain $\Omega=]0,L[^{d}$ contains the “entire” periods $\varepsilon Y$, thus

[TABLE]

For all $z\in{\mathbb{R}}^{3}$, let $[z]$ be the unique integer such that $z-[z]\in Y$. We may write $z=[z]+\{z\}$ for all $z\in{\mathbb{R}}^{3}$, so that for all $\varepsilon>0$, we get the unique decomposition

[TABLE]

Based on this decomposition, the periodic unfolding operator $\mathcal{T}_{\varepsilon}\,{}:L^{2}(\Omega;{\mathbb{R}})\rightarrow L^{2}(\Omega\times Y;{\mathbb{R}})$ is defined as follows: for any function $v\in L^{1}(\Omega;{\mathbb{R}})$, extended to $L^{1}({\mathbb{R}}^{3};{\mathbb{R}})$ by zero outside $\Omega$, i.e. $v=0$ in ${\mathbb{R}}^{3}\setminus\Omega$,

[TABLE]

Unfolding operator $\mathcal{T}_{\varepsilon}$ has the following three important properties: For all functions $\psi$ and $\chi$:

[TABLE]

By $\mathcal{M}_{Y}(\cdot)$ we denote the average operator over $Y$, if $\mathcal{T}_{\varepsilon}\left({w^{\varepsilon}}\right)\rightharpoonup\hat{w}$ weakly in $L^{p}(\Omega\times Y)$, then ${w^{\varepsilon}}\rightharpoonup\mathcal{M}_{Y}(\hat{w})$ weakly in $L^{p}(\Omega)$. For any $D\subset Y$, ${\small\sim}\kern-9.49997pt\int_{D}=\frac{1}{|Y|}\int_{D}$; the analogical notation is employed for any $A\subset Z$, thus ${\small\sim}\kern-9.49997pt\int_{A}=\frac{1}{|Z|}\int_{A}$. Further, for any $D\subset Y$, ${\bf{H}}_{\#}^{1}(D)$ is the Sobolev space ${\textbf{{W}}}^{1,2}(Y)={\bf{H}}^{1}(Y)$ of vector-valued Y-periodic functions (indicated by the subscript $\#$).

Importantly, the unfolding operator also transforms the integration in domain $\Omega$ to $\Omega\times Y$, so that standard means of the weak convergence in Lebesgue spaces $L^{q}(\Omega\times Y)$ can be employed. For more details see [9].

Appendix B Introduction of parameter $\varepsilon$ into the dimensionless system

Here we clarify, how the scale parameter $\varepsilon$ (see Section 3.2) is introduced into the system of equations (2)-(8) along with the derivation of the dimensionless form (20)-(28).

Modified Stokes problem

The characteristic pressure $p_{c}$ is expressed using the ideal gas law,

[TABLE]

Then, by inserting the dimensionless quantities (17) and the dimensionless operator (16) into (9), we get

[TABLE]

hence (18) introduces the dimensionless force $\mathbf{f}^{*}=\frac{L_{c}}{p_{c}}{\textbf{{f}}}$. Upon substituting expressions $\Psi_{c}=k_{B}T/e$ and (89) into (90), we get $\frac{ec_{c}\Psi_{c}}{p_{c}}=1$. Further, according to [16], the ration between the velocity and pressure magnitudes $\lambda_{c}:=\frac{v_{c}\eta_{f}}{L_{c}p_{c}}$ is obtained by the dimensional analysis of the Darcy law which can also represent the viscous flow in pores. This yields

[TABLE]

where $k$ denotes the intrinsic permeability (units [m2]) depending only on the size of the micropores, $k\sim l^{2}$, hence holds and

[TABLE]

Consequently from (89)-(91), the dimensionless form (9) reads,

[TABLE]

Electrostatics

Upon substituting (17) and (16) in the Gauss-Poisson equation (5), it yields

[TABLE]

Using the Debye length definition (19) and parameter $\gamma=l^{2}(\lambda_{D}\sum_{{\beta=1}}^{2}z_{\beta}^{2})^{-1}$, we may express the characteristic concentration $c_{c}$ as

[TABLE]

By substituting $c_{c}$ and $\Psi_{c}$ into (93), we get

[TABLE]

so that $c_{c}L_{c}=\varepsilon^{2}$, hence (93) reads

[TABLE]

Similarly, by inserting (16)-(18) into (11), we get

[TABLE]

After a few easy adjustments we get its dimensionless form as follows

[TABLE]

where $N_{\sigma}=\frac{el\Sigma_{c}}{\mathcal{E}k_{B}T}$ is the ratio between electrical and thermal energy and it is usually of order $\mathcal{O}(1)$ in $\varepsilon$, [19].