Homogenized Balance Equations for Nonlinear Poroelastic Composites

Abstract

Within this work, we upscale the equations that describe the pore-scale behaviour of nonlinear porous elastic composites, using the asymptotic homogenization technique in order to derive the macroscale effective governing equations. A porous hyperelastic composite can be thought of as being comprised of a matrix interacting with a number of subphases and percolated by a fluid flowing in the pores (which is chosen to be Newtonian and incompressible here). A general nonlinear macroscale model is derived and is then specified for a particular choice of strain energy function, namely the de Saint-Venant function. This leads to a macroscale system of PDEs, which is of poroelastic type with additional terms and transformations to account for the nonlinear behaviour of the material. Our new porohyperelastic-type model describes the effective behaviour of nonlinear porous composites by prescribing the stress balance equations, the conservation of mass and Darcy’s law. The coefficients of these macroscale equations encode the detailed microstructure of the material and are to be found by solving pore-scale differential problems. The model reduces to the following limit cases of (a) linear poroelastic composites when the deformation gradient approaches the identity, (b) nonlinear composites when there are no pores and (c) nonlinear poroelasticity when only the matrix–fluid interaction is considered. This model is applicable when the interactions between various hyperelastic solid phases occur at the pore-scale, as in biological tissues such as artery walls, the myocardium, lungs and liver.

1. Introduction

The theory of poroelasticity can be used to describe the effective mechanical behaviour of porous elastic materials with fluid flowing in the pores. This modelling approach can be applied to materials undergoing both linear and nonlinear deformations. The linear theory was developed by Biot in [1,2,3,4]. Due to the desire to apply poroelasticity to biological tissues, the theory was adapted to nonlinearities in [5,6,7,8]. The poroelastic modelling framework is applicable to a wide variety of physical scenarios, in particular to biological tissues, where the deformations are in general nonlinear. For example, a poroelastic approach has been taken to model organs such as lungs (see, e.g., [9]) and to consider the perfused myocardium [10,11]. Poroelasticity has also been applied to studying the artery walls (see, e.g., [12,13]). Nonlinear poroelasticity has also been applied to modelling tumour growth [14] and in imaging to locate tumours in an incompressible medium [15]. It is also of interest in porous thermoelasticity (see [16]). For a general overview of the micromechanics of porous media, we refer the reader to [17], where an overview of upscaling techniques, the linear theory of porous media and also the extension of the analysis to the nonlinear homogenization of a large range of scenarios including strength homogenization, nonsaturated microporomechanics, microporoplasticity and microporofracture and microporodamage theory were discussed.

We consider hyperelastic deformable porous media, and in general, these types of materials have a variety of structural features characterised over multiple scales. We have a scale (called the pore-scale) where the interactions between the fluid and the solid take place, and this scale is much smaller than the scale at which we describe the whole material (which is called the macroscale).

The coupled fluid–structure balance equations that describe the material on the pore-scale can be used in an upscaling process to obtain the macroscale governing equations of a poroelastic material. The upscaling process can be carried out by a variety of homogenization techniques. These homogenization techniques include effective medium theory, mixture theory, volume averaging and asymptotic homogenization. These techniques were discussed and described in [18,19].

We develop our analysis by means of the asymptotic homogenization technique [20,21,22,23]. This technique requires a sharp length-scale separation, between the pore-scale and the macroscale, to enable the decoupling of spatial variations that exist in the system. This technique expresses the relevant fields from the pore-scale fluid–structure interaction problem as power series of the ratio between the two different typical length scales. The balance equations obtained via this method describe the effective behaviour of the material in the homogenized macroscale domain. A key feature of this technique is that the coefficients of the macroscale model encode information about the microstructure of the material, and they are to be computed by solving differential problems at the pore-scale. In the literature, the asymptotic homogenization technique has been applied to various physical systems. For example, the asymptotic homogenization technique has been used to derive Biot’s poroelasticity (see [24,25]) and extended to study vascularised poroelastic materials [26], poroelastic composites [27] and modelling of the bones [28].

In the asymptotic homogenization literature, the majority of the applications focus on linearised balance equations. This is due to the fact that in the linearised case, it is possible to fully decouple the pore-scale and the macroscale (under some simplifying assumptions). This decoupling then leads to a large reduction in the computational complexity of the system. In the literature, the homogenization of systems involving nonlinear mechanics is generally carried out by other homogenization techniques, such as average field techniques (see, e.g., [29,30] and the references therein). These other techniques do not provide a precise description of the model coefficients in the way that the asymptotic homogenization technique can and, instead, provide bounds for the model coefficients. There are currently only a few examples of using asymptotic homogenization in situations where the material undergoes large deformations. In [31], the effective poroelastic model of Biot was extended to a nonlinear Biot model that includes pore-scale deformation. In [32], a system of effective equations that describe the flow, elastic deformation and transport in an active medium was derived. The authors considered the spatial homogenization of a coupled transport and fluid–structure interaction model. In [33], the asymptotic homogenization technique was applied to the equations that describe the dynamics of a heterogeneous material with an evolving microstructure to obtain a set of effective equations. The heterogeneous body is assumed to be composed of two hyperelastic materials, and the evolution of the microstructure is through plastic-like distortions.

We investigate materials that are subject to large deformations that have the underlying microstructure comprised of both a hyperelastic-fluid-filled porous matrix and a number of embedded hyperelastic subphases (fibres or inclusions). We then assume that both the matrix and the subphases interact with each other and the fluid at the pore-scale. This type of structure can be described as a poroelastic composite material that undergoes large deformations. For example, this is applicable to artery walls, which can be considered as a composite nonlinear elastic material consisting of a matrix with two families of symmetrically arranged embedded collagen and elastin fibres as well as fibroblast cells that interact with the fluid that is flowing in the pores [12,34,35,36,37]. If we wish to consider the artery walls in our framework, then the matrix can be identified with our matrix in the model, and the fibres and fibroblast cells could all be considered as the elastic subphases that are embedded in the matrix. The fluid in this setting would be the water that flows in the pores of the matrix. This modelling approach is applicable to the myocardium in the heart, which is a nonlinear elastic porous structure that consists of a matrix with muscle cells, fibroblasts, collagen fibres and embedded blood vessels [10,11,38]. Again, here, the various muscle cells, fibroblasts and fibres would all be the elastic subphases embedded in the matrix. By considering these systems as poroelastic composites, we are able to account for the mechanical contribution of each of the various phases individually.

Here, we generalise [27] to account for nonlinear deformations by using the asymptotic homogenization technique to upscale the interaction between a hyperelastic porous matrix where there is an incompressible Newtonian fluid flowing in the pores and a number of embedded hyperelastic subphases. We make the assumption that both the hyperelastic matrix and subphases interact with the fluid that is flowing in the pores. We assume that the length scale where the individual hyperelastic subphases are clearly visible from the surrounding matrix is comparable to the pore size. We therefore determine that this scale will be the pore-scale of the material, that is the distance between adjacent subphases is comparable with the size of the pores. This length is assumed to be much smaller than the size of the whole domain, which is denoted the macroscale. The upscaling process is then used taking into account the continuity of the stresses and elastic displacements across the interfaces between the matrix and the subphases, as well as the continuity of the stresses and velocities across the interfaces between the fluid and solid domains. Furthermore, an appropriate coordinate transformation is carried out on some quantities to formulate the full problem in the undeformed/reference configuration, i.e., by solely using Lagrangian coordinates. The resulting system of governing equations is of poroelastic type and is a generalization of the formulations for (a) multiphase elastoplastic composites [33] in the limit of no plastic distortions and (b) the formulations for hyperelastic porous media [31] in the limit of only one elastic phase. It is also a natural extension to the formulation for linear poroelastic composites [27]. All three of these formulations are recovered as particular cases by assuming that: (a) our matrix is not porous; (b) that no elastic phase other than the matrix is present; and (c) by performing a linearisation. The coefficients of the model encode the properties of the microstructure and are computed by solving local differential problems.

The paper is organised as follows. In Section 2, we formulate the fluid–structure interaction problem that characterises the behaviour of the hyperelastic porous matrix, the hyperelastic subphases and the fluid flowing in the pores. We also perform a change of coordinates, which allows the fluid–structure interaction problem to be formulated in the reference configuration. In Section 3, we apply the asymptotic homogenization technique to the system of PDEs that were described in Section 2 and determine the new model that describes the effective macroscale mechanical behaviour of nonlinear poroelastic composites. In Section 4, we discuss the general nonlinear macroscale model before prescribing a specific strain energy function, namely the de Saint-Venant, for the material and then give a physical interpretation of the novel terms. We also consider particular cases for our new model and are able to obtain previously known models from the literature. Section 5 concludes our work by highlighting and discussing the limitations of the current model and by providing further directions in which the model could be extended for particular biological applications.

2. Formulation of the Fluid–Structure Interaction Problem

We assume that we have a continuum body that is not subject to any surface or body forces to have a reference configuration, which we denote by ${\mathcal{B}}_0\in {\mathbb{R}}^3$. The body ${\mathcal{B}}_0$ has a periodic microstructure that consists of the union of a porous hyperelastic matrix ${\Omega }_{\mathrm{II}}^0$, an interconnected fluid compartment ${\Omega }_{\mathrm{f}}^0$ and a set ${\Omega }_{\mathrm{I}}^0$ of K disjoint hyperelastic subphases ${\Omega }_{\alpha }^0$, where:

$${\Omega }_{\mathrm{I}}^0=\bigcup _{\alpha =1}^K{\Omega }_{\alpha }^0,$$

(1)

and ${\overline{\Omega }}^0={\overline{\Omega }}_{\mathrm{I}}^0\cup {\overline{\Omega }}_{\mathrm{II}}^0\cup {\overline{\Omega }}_{\mathrm{f}}^0$. We note that when any domain, interface or normal vector has the superscript 0 then it is in the reference configuration. We provide a sketch of this structure as shown in Figure 1. The body undergoes a deformation described by the deformation function $\chi$, and the deformed body is denoted by ${\mathcal{B}}_t$. Each point $\mathbf{x}\in {\mathcal{B}}_t$ is such that $\mathbf{x}=\chi (\mathbf{X},t)$ with $\mathbf{X}\in {\mathcal{B}}_0$ being the point in the reference configuration. We assume that the periodicity of the body’s microstructure is preserved during the deformation, and we now denote the deformed solid porous matrix as ${\Omega }_{\mathrm{II}}^t$, the deformed connected fluid compartment as ${\Omega }_{\mathrm{f}}^t$ and the set of K disjoint deformed subphases ${\Omega }_{\alpha }^t$ as ${\Omega }_{\mathrm{I}}^t$ in ${\mathcal{B}}_t$. When any domain, interface or normal vector has the superscript t, this denotes the current configuration. The deformation from ${\mathcal{B}}_0$ to ${\mathcal{B}}_t$ is described by the deformation gradient $\mathbf{F}=\mathrm{Grad}\chi$. Figure 1 highlights the description of this structure pictorially.

We now wish to describe the equations for the fluid and the solid compartments and the interface conditions in our structure. These equations however are not all in the same coordinate systems. We wish to work in the reference configuration so we have to perform a pull back to some of the equations that are described in the current configuration. We begin by describing all equations in their natural coordinate systems before discussing the coordinate transformations.

The balance equations for the various hyperelastic domains ${\Omega }_{\alpha }$ and ${\Omega }_{\mathrm{II}}$ can be written as, $\forall \alpha =1,\cdots ,K$,

$${\nabla }_{\mathbf{X}}·{\sigma }_{\alpha }=0\mathrm{in}{\Omega }_{\alpha }^0,$$

(2)

and:

$${\nabla }_{\mathbf{X}}·{\sigma }_{\mathrm{II}}=0\mathrm{in}{\Omega }_{\mathrm{II}}^0,$$

(3)

where we neglect any volume forces and inertia. The tensors ${\sigma }_{\alpha }$ and ${\sigma }_{\mathrm{II}}$ are the solid stresses. Each subphase ${\Omega }_{\alpha }$ has the stress tensor ${\sigma }_{\alpha }$, and the porous matrix ${\Omega }_{\mathrm{II}}$ has solid stress tensor ${\sigma }_{\mathrm{II}}$. The matrix and the subphases are anisotropic, hyperelastic solids, and therefore, the constitutive laws for ${\sigma }_{\alpha }$ and ${\sigma }_{\mathrm{II}}$ are given in terms of strain energy functions,

$$\begin{array}{cc}\hfill {\sigma }_{\alpha }=\frac{\partial {\psi }_{\alpha }}{\partial {\mathbf{F}}_{\alpha }}& \mathrm{in}{\Omega }_{\alpha }^0,\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{II}}=\frac{\partial {\psi }_{\mathrm{II}}}{\partial {\mathbf{F}}_{\mathrm{II}}}& \mathrm{in}{\Omega }_{\mathrm{II}}^0,\hfill \end{array}$$

(5)

where ${\mathbf{F}}_{\alpha }$ and ${\mathbf{F}}_{\mathrm{II}}$ are the deformation gradients in the subphases and the matrix, respectively, and ${\psi }_{\alpha }$ and ${\psi }_{\mathrm{II}}$ are the strain energy functions in the subphases and matrix, respectively. We do not define a specific strain energy function at this stage and wait until Section 4.1 to specify it.

We also require equations for the fluid phase. The balance equation is given as:

$${\nabla }_{\mathbf{x}}·{\sigma }_{\mathrm{f}}=0\mathrm{in}{\Omega }_{\mathrm{f}}^t,$$

(6)

where we denoted the fluid stress tensor by ${\sigma }_{\mathrm{f}}$. The fluid is assumed to be an incompressible Newtonian fluid and, so, has the constitutive law:

$${\sigma }_{\mathrm{f}}=-p\mathbf{I}+\mu (\left({\nabla }_{\mathbf{x}}\mathbf{v}\right)+{\left({\nabla }_{\mathbf{x}}\mathbf{v}\right)}^{\mathrm{T}})\mathrm{in}{\Omega }_{\mathrm{f}}^t,$$

(7)

where $\mathbf{v}$ is the fluid velocity, p is the fluid pressure and $\mu$ is the fluid viscosity. As we consider an incompressible fluid, the incompressibility condition reads:

$${\nabla }_{\mathbf{x}}·\mathbf{v}=0\mathrm{in}{\Omega }_{\mathrm{f}}^t.$$

(8)

We can substitute the constitutive law (7) into the balance Equation (6), and then, by using (8), we obtain the Stokes problem:

$$0=-{\nabla }_{\mathbf{x}}p+\mu {\nabla }_{\mathbf{x}}·(\left({\nabla }_{\mathbf{x}}\mathbf{v}\right)+{\left({\nabla }_{\mathbf{x}}\mathbf{v}\right)}^{\mathrm{T}})\mathrm{in}{\Omega }_{\mathrm{f}}^t.$$

(9)

In order to set up an appropriate fluid–structure interaction problem, we require interface conditions between the fluid and the various solid phases and also interface conditions between the various solid subphases and the matrix. The conditions we impose are the continuity of velocities, the continuity of stresses and the continuity of displacements. This presents an issue since the fluid and solid equations are described in different coordinate systems. The fluid equations are currently presented in Eulerian coordinates and the solid equations in Lagrangian coordinates. It is not possible to properly express continuity on the interface between the fluid and solid whilst the governing equations are in different coordinate systems. For this reason, we only describe the interface conditions between the elastic subphases and the matrix here and wait until Section 2.1 to describe the interface conditions between the fluid and solids once we have formulated all equations in Lagrangian coordinates. We define the interface between each elastic subphase and the matrix as ${\Gamma }_{\alpha \mathrm{II}}^0:=\partial {\Omega }_{\alpha }^0\cap \partial {\Omega }_{\mathrm{II}}^0$ and impose continuity of stresses and displacements, namely:

$$\begin{array}{cc}\hfill {\sigma }_{\alpha }{\mathbf{n}}_{\alpha \mathrm{II}}^0={\sigma }_{\mathrm{II}}{\mathbf{n}}_{\alpha \mathrm{II}}^0& \mathrm{on}{\Gamma }_{\alpha \mathrm{II}}^0,\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill {\mathbf{u}}_{\alpha }={\mathbf{u}}_{\mathrm{II}}& \mathrm{on}{\Gamma }_{\alpha \mathrm{II}}^0,\hfill \end{array}$$

(11)

$\forall \alpha =1,\cdots ,K$, where we define the unit vector normal to the interface ${\Gamma }_{\alpha \mathrm{II}}^0$ as ${\mathbf{n}}_{\alpha \mathrm{II}}^0$, and it is pointing into the subphase ${\Omega }_{\alpha }^0$ and where ${\mathbf{u}}_{\alpha }$ and ${\mathbf{u}}_{\mathrm{II}}$ are the elastic displacements in each subphase and the matrix, respectively.

We describe the deformation from ${\mathcal{B}}_0$ to ${\mathcal{B}}_t$ by the deformation gradient $\mathbf{F}=\mathrm{Grad}\chi$. For convenience, we use the notation:

$$\mathbf{F}=\left\{\begin{array}{c}\begin{array}{c}{\mathbf{F}}_{\alpha }\mathrm{in}{\Omega }_{\alpha }^0,\hfill \\ {\mathbf{F}}_{\mathrm{II}}\mathrm{in}{\Omega }_{\mathrm{II}}^0,\hfill \\ {\mathbf{F}}_{\mathrm{f}}\mathrm{in}{\Omega }_{\mathrm{f}}^0.\hfill \end{array}\hfill \end{array}\right.$$

(12)

We note that our deformation gradient is in general discontinuous. It is however useful here to consider the alternative definition of the deformation gradient:

$$\mathbf{F}=\mathbf{I}+\mathrm{Grad}\mathbf{u},$$

(13)

where $\mathbf{I}$ is the identity tensor and $\mathrm{Grad}\mathbf{u}$ is the gradient operator of the elastic displacement. We use the notation Grad with a capital G to denote the gradient in Lagrangian coordinates and grad for Eulerian coordinates. We can specialise relationship (13) in each of the reference elastic subphases as follows:

$$\mathbf{F}=\left\{\begin{array}{c}\begin{array}{c}{\mathbf{F}}_{\alpha }=\mathbf{I}+\mathrm{Grad}{\mathbf{u}}_{\alpha }\mathrm{in}{\Omega }_{\alpha }^0,\hfill \\ {\mathbf{F}}_{\mathrm{II}}=\mathbf{I}+\mathrm{Grad}{\mathbf{u}}_{\mathrm{II}}\mathrm{in}{\Omega }_{\mathrm{II}}^0,\hfill \end{array}\hfill \end{array}\right.$$

(14)

We should note that the sketch here highlights a number of possible arrangements for the subphases. We can have subphases fully embedded in the matrix, fully embedded in the fluid or in contact with both the matrix and the fluid. We set up our fluid–structure interaction with the assumption that all elastic phases are in contact with each other and the fluid.

2.1. Fluid–Structure Interaction in Lagrangian Coordinates

Within this section, we apply a coordinate transformation to the equations in the fluid–structure interaction (FSI) in order to obtain a full system of PDEs that describe the structure in the reference configuration. We do this in order to preserve the local periodicity of the microstructure, as was done in [32]. We define the Piola transformation, $\mathsf{G}$, and the Jacobian, J, by:

$$J=\det \mathbf{F}\mathrm{and}\mathsf{G}=J{\mathbf{F}}^{-1},$$

(15)

where we have that:

$$\mathsf{G}=\left\{\begin{array}{c}\begin{array}{c}{\mathsf{G}}_{\alpha }=J_{\alpha }{\mathbf{F}}_{\alpha }^{-1}\mathrm{in}{\Omega }_{\alpha }^0,\hfill \\ {\mathsf{G}}_{\mathrm{II}}=J_{\mathrm{II}}{\mathbf{F}}_{\mathrm{II}}^{-1}\mathrm{in}{\Omega }_{\mathrm{II}}^0,\hfill \\ {\mathsf{G}}_{\mathrm{f}}=J_{\mathrm{f}}{\mathbf{F}}_{\mathrm{f}}^{-1}\mathrm{in}{\Omega }_{\mathrm{f}}^0.\hfill \end{array}\hfill \end{array}\right.$$

(16)

Here, we wish to make a remark regarding the continuity of the Piola transformation.

Remark 1

(Continuity of ${\mathsf{G}}^{\mathrm{T}}$). Here, we wish to explain the continuity of ${\mathsf{G}}^{\mathrm{T}}$ across the various interfaces appearing in our structure. We begin with Nanson’s formula:

$$\mathbf{n}\mathrm{da}=J{\mathbf{F}}^{-\mathrm{T}}\mathbf{N}\mathrm{dA},$$

(17)

where $\mathbf{n}$ is a general unit normal and $\mathrm{da}$ is a general area element in the current configuration and $\mathbf{N}$ is a general unit normal and $\mathrm{dA}$ is a general area element in the reference configuration. Using our notation for the Piola transform, we can rewrite this as:

$$\mathbf{n}\mathrm{da}={\mathsf{G}}^{\mathrm{T}}\mathbf{N}\mathrm{dA}.$$

(18)

By using this relationship, we are able to deduce the continuities on our various interfaces ${\Gamma }_{\alpha }$, ${\Gamma }_{\mathrm{II}}$ and ${\Gamma }_{\alpha \mathrm{II}}$, that is,

$$\begin{array}{c}{\mathsf{G}}_{\alpha }^{\mathrm{T}}{\mathbf{n}}_{\alpha }^0={\mathsf{G}}_{\mathrm{f}}^{\mathrm{T}}{\mathbf{n}}_{\alpha }^{0}on{\Gamma }_{\alpha },{\mathsf{G}}_{\mathrm{II}}^{\mathrm{T}}{\mathbf{n}}_{\mathrm{II}}^0={\mathsf{G}}_{\mathrm{f}}^{\mathrm{T}}{\mathbf{n}}_{\mathrm{II}}^{0}on{\Gamma }_{\mathrm{II}},\\ and{\mathsf{G}}_{\alpha }^{\mathrm{T}}{\mathbf{n}}_{\alpha \mathrm{II}}^0={\mathsf{G}}_{\mathrm{II}}^{\mathrm{T}}{\mathbf{n}}_{\alpha \mathrm{II}}^{0}on{\Gamma }_{\alpha \mathrm{II}}.\end{array}$$

(19)

These formulas are used in later sections of this work.

We use the following formulas to make the change of the coordinate systems. For a scalar $\xi$, a vector $\mathbf{a}$ and a tensor $\mathbf{A}$, we have:

$${\nabla }_{\mathbf{x}}\xi ={\mathbf{F}}^{-\mathrm{T}}{\nabla }_{\mathbf{X}}\xi \mathrm{and}{\nabla }_{\mathbf{x}}\mathbf{a}=\left({\nabla }_{\mathbf{X}}\mathbf{a}\right){\mathbf{F}}^{-1}.$$

(20)

By again using Nanson’s formula:

$$\mathbf{n}\mathrm{da}=J{\mathbf{F}}^{-\mathrm{T}}\mathbf{N}\mathrm{dA},$$

(21)

the transformation rule for a general volume element:

$$\mathrm{d}v=J\mathrm{d}V,$$

(22)

and by applying the divergence theorem, we obtain the coordinate changes:

$${\nabla }_{\mathbf{x}}·\mathbf{a}=\frac{1}{J}{\nabla }_{\mathbf{X}}·\left(\mathsf{G}\mathbf{a}\right)\mathrm{and}{\nabla }_{\mathbf{x}}·\mathbf{A}=\frac{1}{J}{\nabla }_{\mathbf{X}}·\left(\mathbf{A}{\mathsf{G}}^{\mathrm{T}}\right).$$

(23)

We are now able to use (20) and (23) to write our FSI problem in the reference configuration. The equations governing the fluid in the reference configuration are therefore given by:

$${\nabla }_{\mathbf{X}}·\left({\sigma }_{\mathrm{f}}{\mathsf{G}}_{\mathrm{f}}^{\mathrm{T}}\right)=0\mathrm{in}{\Omega }_{\mathrm{f}}^0,$$

(24)

with the fluid stress tensor ${\sigma }_{\mathrm{f}}$ transformed as:

$${\sigma }_{\mathrm{f}}{\mathsf{G}}_{\mathrm{f}}^{\mathrm{T}}=-p{\mathsf{G}}_{\mathrm{f}}^{\mathrm{T}}+\mu (\left({\nabla }_{\mathbf{X}}\mathbf{V}\right){\mathbf{F}}_{\mathrm{f}}^{-1}{\mathsf{G}}_{\mathrm{f}}^{\mathrm{T}}+{\mathbf{F}}_{\mathrm{f}}^{-\mathrm{T}}{\left({\nabla }_{\mathbf{X}}\mathbf{V}\right)}^{\mathrm{T}}{\mathsf{G}}_{\mathrm{f}}^{\mathrm{T}})\mathrm{in}{\Omega }_{\mathrm{f}}^0,$$

(25)

to correspond to the fluid balance equation, where $\mathbf{V}$ is the fluid velocity in the reference configuration. The incompressibility condition also transforms as follows:

$${\nabla }_{\mathbf{X}}·\left({\mathsf{G}}_{\mathrm{f}}\mathbf{V}\right)=0\mathrm{in}{\Omega }_{\mathrm{f}}^0.$$

(26)

The Stokes’ problem in the reference configuration is obtained by substituting (25) into (24) and using (26), that is,

$$0=-{\mathsf{G}}_{\mathrm{f}}^{\mathrm{T}}{\nabla }_{\mathbf{X}}p+\mu {\nabla }_{\mathbf{X}}·(\left({\nabla }_{\mathbf{X}}\mathbf{V}\right){\mathsf{G}}_{\mathrm{f}}{\mathbf{F}}_{\mathrm{f}}^{-\mathrm{T}}+{\mathsf{G}}_{\mathrm{f}}^{\mathrm{T}}{\left({\nabla }_{\mathbf{X}}\mathbf{V}\right)}^{\mathrm{T}}{\mathbf{F}}_{\mathrm{f}}^{-\mathrm{T}})\mathrm{in}{\Omega }_{\mathrm{f}}^0.$$

(27)

Now that all our governing equations are described in the Lagrangian framework, we are able to describe the interface conditions between the fluid and solid phases that we require to set up an appropriate FSI problem. We note that these interface conditions take place on the interfaces between the fluid and the elastic subphases in the reference configuration. Defining the interfaces between the fluid phase and each of the $\alpha$ subphases in the reference configuration as ${\Gamma }_{\alpha }^0:=\partial {\Omega }_{\alpha }^0\cap \partial {\Omega }_{\mathrm{f}}^0$ and defining the interface between the matrix and the fluid phase, again in the reference configuration, as ${\Gamma }_{\mathrm{II}}^0:=\partial {\Omega }_{\mathrm{II}}^0\cap \partial {\Omega }_{\mathrm{f}}^0$, we can impose the continuity of the velocities and stresses on the various interfaces. We therefore have:

$$\begin{array}{cc}\hfill \frac{\partial {\mathbf{u}}_{\alpha }}{\partial t}=\mathbf{V}& \mathrm{on}{\Gamma }_{\alpha }^0,\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{f}}{\mathsf{G}}_{\mathrm{f}}^{\mathrm{T}}{\mathbf{n}}_{\alpha }^0={\sigma }_{\alpha }{\mathbf{n}}_{\alpha }^0& \mathrm{on}{\Gamma }_{\alpha }^0,\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill \frac{\partial {\mathbf{u}}_{\mathrm{II}}}{\partial t}=\mathbf{V}& \mathrm{on}{\Gamma }_{\mathrm{II}}^0,\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{f}}{\mathsf{G}}_{\mathrm{f}}^{\mathrm{T}}{\mathbf{n}}_{\mathrm{II}}^0={\sigma }_{\mathrm{II}}{\mathbf{n}}_{\mathrm{II}}^0& \mathrm{on}{\Gamma }_{\mathrm{II}}^0,\hfill \end{array}$$

(31)

$\forall \alpha =1,\cdots ,K$. We have that $\partial {\mathbf{u}}_{\alpha }/\partial t$ and $\partial {\mathbf{u}}_{\mathrm{II}}/\partial t$ are the solid velocities in each of the subphases and the matrix, respectively. We have that ${\mathbf{n}}_{\alpha }^0$ and ${\mathbf{n}}_{\mathrm{II}}^0$ are the unit outward normals to the interfaces ${\Gamma }_{\alpha }^0$ and ${\Gamma }_{\mathrm{II}}^0$, respectively.

Our complete FSI problem in the reference configuration is therefore given by (2)–(5), (10), (11), (14) and (24)–(31).

Within the next section, we carry out our analysis by first nondimensionalising the system of partial differential equations (PDEs) that we formulate in the reference configuration within this section. We introduce two well-separated length scales that allow us to apply the two-scale asymptotic homogenization technique to the nondimensionalised PDEs. This allows the derivation of the macroscale governing equations.

3. The Asymptotic Homogenization Method

Here, we summarise the fluid–structure interaction problem in the reference configuration that we introduced in the previous section. We will then be ready to perform a multiscale analysis.

$$\begin{array}{cc}\hfill {\nabla }_{\mathbf{X}}·{\sigma }_{\alpha }=0& \mathrm{in}{\Omega }_{\alpha }^0,\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill {\nabla }_{\mathbf{X}}·{\sigma }_{\mathrm{II}}=0& \mathrm{in}{\Omega }_{\mathrm{II}}^0,\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill {\nabla }_{\mathbf{X}}·\left({\sigma }_{\mathrm{f}}{\mathsf{G}}_{\mathrm{f}}^{\mathrm{T}}\right)=0& \mathrm{in}{\Omega }_{\mathrm{f}}^0,\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill {\nabla }_{\mathbf{X}}·\left({\mathsf{G}}_{\mathrm{f}}\mathbf{V}\right)=0& \mathrm{in}{\Omega }_{\mathrm{f}}^0,\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill \frac{\partial {\mathbf{u}}_{\alpha }}{\partial t}=\mathbf{V}& \mathrm{on}{\Gamma }_{\alpha }^0,\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill \frac{\partial {\mathbf{u}}_{\mathrm{II}}}{\partial t}=\mathbf{V}& \mathrm{on}{\Gamma }_{\mathrm{II}}^0,\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{f}}{\mathsf{G}}_{\mathrm{f}}^{\mathrm{T}}{\mathbf{n}}_{\alpha }^0={\sigma }_{\alpha }{\mathbf{n}}_{\alpha }^0& \mathrm{on}{\Gamma }_{\alpha }^0,\hfill \end{array}$$

(38)

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{f}}{\mathsf{G}}_{\mathrm{f}}^{\mathrm{T}}{\mathbf{n}}_{\mathrm{II}}^0={\sigma }_{\mathrm{II}}{\mathbf{n}}_{\mathrm{II}}^0& \mathrm{on}{\Gamma }_{\mathrm{II}}^0,\hfill \end{array}$$

(39)

$$\begin{array}{cc}\hfill {\sigma }_{\alpha }{\mathbf{n}}_{\alpha \mathrm{II}}^0={\sigma }_{\mathrm{II}}{\mathbf{n}}_{\alpha \mathrm{II}}^0& \mathrm{on}{\Gamma }_{\alpha \mathrm{II}}^0.\hfill \end{array}$$

(40)

$$\begin{array}{cc}\hfill {\mathbf{u}}_{\alpha }={\mathbf{u}}_{\mathrm{II}}& \mathrm{on}{\Gamma }_{\alpha \mathrm{II}}^0,\hfill \end{array}$$

(41)

We have that the constitutive relationships for the fluid and the multiple solid phases are given as:

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{f}}{\mathsf{G}}_{\mathrm{f}}^{\mathrm{T}}=-p{\mathsf{G}}_{\mathrm{f}}^{\mathrm{T}}+\mu (\left({\nabla }_{\mathbf{X}}\mathbf{V}\right){\mathbf{F}}_{\mathrm{f}}^{-1}{\mathsf{G}}_{\mathrm{f}}^{\mathrm{T}}+{\mathbf{F}}_{\mathrm{f}}^{-\mathrm{T}}{\left({\nabla }_{\mathbf{X}}\mathbf{V}\right)}^{\mathrm{T}}{\mathsf{G}}_{\mathrm{f}}^{\mathrm{T}})& \mathrm{in}{\Omega }_{\mathrm{f}}^0,\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill {\sigma }_{\alpha }=\frac{\partial {\psi }_{\alpha }}{\partial {\mathbf{F}}_{\alpha }}& \mathrm{in}{\Omega }_{\alpha }^0,\hfill \end{array}$$

(43)

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{II}}=\frac{\partial {\psi }_{\mathrm{II}}}{\partial {\mathbf{F}}_{\mathrm{II}}}& \mathrm{in}{\Omega }_{\mathrm{II}}^0.\hfill \end{array}$$

(44)

We then use the constitutive relationship (42) with the incompressibility constraint (35) in the balance Equation (34) to obtain:

$$\begin{array}{cc}\hfill 0=-{\mathsf{G}}_{\mathrm{f}}^{\mathrm{T}}{\nabla }_{\mathbf{X}}p+\mu {\nabla }_{\mathbf{X}}·(\left({\nabla }_{\mathbf{X}}\mathbf{V}\right){\mathsf{G}}_{\mathrm{f}}{\mathbf{F}}_{\mathrm{f}}^{-\mathrm{T}}+{\mathsf{G}}_{\mathrm{f}}^{\mathrm{T}}{\left({\nabla }_{\mathbf{X}}\mathbf{V}\right)}^{\mathrm{T}}{\mathbf{F}}_{\mathrm{f}}^{-\mathrm{T}})& \mathrm{in}{\Omega }_{\mathrm{f}}^0.\hfill \end{array}$$

(45)

We also have the deformation gradients for each of the solid phases:

$$\begin{array}{cc}\hfill {\mathbf{F}}_{\alpha }=\mathbf{I}+{\nabla }_{\mathbf{X}}{\mathbf{u}}_{\alpha }& \mathrm{in}{\Omega }_{\alpha }^0,\hfill \end{array}$$

(46)

$$\begin{array}{cc}\hfill {\mathbf{F}}_{\mathrm{II}}=\mathbf{I}+{\nabla }_{\mathbf{X}}{\mathbf{u}}_{\mathrm{II}}& \mathrm{in}{\Omega }_{\mathrm{II}}^0,\hfill \end{array}$$

(47)

$\forall \alpha =1,\cdots ,K$. We assume that the system possesses two distinct length scales. The whole domain has a length scale, which we denote by L, and this is referred to as the macroscale. There also exists a second length scale, which we denote by d, and this describes the pore-scale. The length d is comparable to the intersubphase distance and the size of the pores. Therefore, to capture the true difference between the two scales, it is useful to perform a nondimensional analysis of the system of PDEs (32)–(47). This nondimensionalisation is carried out in the next section.

We should note that this fluid–structure interaction is the nonlinear counterpart to the one found in [27] for linear poroelastic composites. We also highlight that this FSI problem is formulated with the intention that all elastic phases are in contact with each other and the fluid.

3.1. Nondimensionalisation

We carry out the nondimensionalisation process by relying on the standard parabolic fluid velocity in the pores, which is quadratic in the pore-scale d and proportional to a given pressure gradient C. This scaling is the classical one, which ensures that a Newtonian fluid flowing in the pores is macroscopically governed by porous media flow equations (e.g., of Darcy’s type) (see [39], where this is discussed). There are of course alternative scalings available for the fluid velocity; however, these scalings do not account for the appropriate effective behaviour of a fluid flow in porous media.

Therefore, we choose the scalings:

$$\begin{array}{c}\mathbf{x}=L{\mathbf{x}}^{\prime },{\mathbf{u}}_{\alpha }=L{\mathbf{u}}_{\alpha }^{\prime },{\mathbf{u}}_{\mathrm{II}}=L{\mathbf{u}}_{\mathrm{II}}^{\prime },\mathbf{V}=\frac{C{d}^2}{\mu }{\mathbf{V}}^{\prime },p=CL{p}^{\prime },{\sigma }_{\alpha }=CL{{\sigma }_{\alpha }}^{\prime },\\ {\sigma }_{\mathrm{II}}=CL{{\sigma }_{\mathrm{II}}}^{\prime }.\end{array}$$

(48)

We then use (48), and the gradient operator becomes:

$$\nabla =\frac{1}{L}{\nabla }^{\prime },$$

(49)

and the nondimensionalised form of the system of the PDEs (32)–(41) is given by,

$$\begin{array}{cc}\hfill {\nabla }_{\mathbf{X}}·{\sigma }_{\alpha }=0& \mathrm{in}{\Omega }_{\alpha }^0,\hfill \end{array}$$

(50)

$$\begin{array}{cc}\hfill {\nabla }_{\mathbf{X}}·{\sigma }_{\mathrm{II}}=0& \mathrm{in}{\Omega }_{\mathrm{II}}^0,\hfill \end{array}$$

(51)

$$\begin{array}{cc}\hfill {\nabla }_{\mathbf{X}}·\left({\sigma }_{\mathrm{f}}{\mathsf{G}}_{\mathrm{f}}^{\mathrm{T}}\right)=0& \mathrm{in}{\Omega }_{\mathrm{f}}^0,\hfill \end{array}$$

(52)

$$\begin{array}{cc}\hfill {\nabla }_{\mathbf{X}}·\left({\mathsf{G}}_{\mathrm{f}}\mathbf{V}\right)=0& \mathrm{in}{\Omega }_{\mathrm{f}}^0,\hfill \end{array}$$

(53)

$$\begin{array}{cc}\hfill \frac{\partial {\mathbf{u}}_{\alpha }}{\partial t}=\mathbf{V}& \mathrm{on}{\Gamma }_{\alpha }^0,\hfill \end{array}$$

(54)

$$\begin{array}{cc}\hfill \frac{\partial {\mathbf{u}}_{\mathrm{II}}}{\partial t}=\mathbf{V}& \mathrm{on}{\Gamma }_{\mathrm{II}}^0,\hfill \end{array}$$

(55)

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{f}}{\mathsf{G}}_{\mathrm{f}}^{\mathrm{T}}{\mathbf{n}}_{\alpha }^0={\sigma }_{\alpha }{\mathbf{n}}_{\alpha }^0& \mathrm{on}{\Gamma }_{\alpha }^0,\hfill \end{array}$$

(56)

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{f}}{\mathsf{G}}_{\mathrm{f}}^{\mathrm{T}}{\mathbf{n}}_{\mathrm{II}}^0={\sigma }_{\mathrm{II}}{\mathbf{n}}_{\mathrm{II}}^0& \mathrm{on}{\Gamma }_{\mathrm{II}}^0,\hfill \end{array}$$

(57)

$$\begin{array}{cc}\hfill {\sigma }_{\alpha }{\mathbf{n}}_{\alpha \mathrm{II}}^0={\sigma }_{\mathrm{II}}{\mathbf{n}}_{\alpha \mathrm{II}}^0& \mathrm{on}{\Gamma }_{\alpha \mathrm{II}}^0,\hfill \end{array}$$

(58)

$$\begin{array}{cc}\hfill {\mathbf{u}}_{\alpha }={\mathbf{u}}_{\mathrm{II}}& \mathrm{on}{\Gamma }_{\alpha \mathrm{II}}^0.\hfill \end{array}$$

(59)

We also nondimensionalise the constitutive relationships for the fluid and the solid, and these are given by:

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{f}}{\mathsf{G}}_{\mathrm{f}}^{\mathrm{T}}=-p{\mathsf{G}}_{\mathrm{f}}^{\mathrm{T}}+{ϵ}^2(\left({\nabla }_{\mathbf{X}}\mathbf{V}\right){\mathbf{F}}_{\mathrm{f}}^{-1}{\mathsf{G}}_{\mathrm{f}}^{\mathrm{T}}+{\mathbf{F}}_{\mathrm{f}}^{-\mathrm{T}}{\left({\nabla }_{\mathbf{X}}\mathbf{V}\right)}^{\mathrm{T}}{\mathsf{G}}_{\mathrm{f}}^{\mathrm{T}})& \mathrm{in}{\Omega }_{\mathrm{f}}^0,\hfill \end{array}$$

(60)

$$\begin{array}{cc}\hfill {\sigma }_{\alpha }=\frac{\partial {\psi }_{\alpha }}{\partial {\mathbf{F}}_{\alpha }}& \mathrm{in}{\Omega }_{\alpha }^0,\hfill \end{array}$$

(61)

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{II}}=\frac{\partial {\psi }_{\mathrm{II}}}{\partial {\mathbf{F}}_{\mathrm{II}}}& \mathrm{in}{\Omega }_{\mathrm{II}}^0.\hfill \end{array}$$

(62)

We also have the nondimensionalised fluid balance equation and solid deformation gradients given by:

$$\begin{array}{cc}\hfill 0=-{\mathsf{G}}_{\mathrm{f}}^{\mathrm{T}}{\nabla }_{\mathbf{X}}p+{ϵ}^2{\nabla }_{\mathbf{X}}·(\left({\nabla }_{\mathbf{X}}\mathbf{V}\right){\mathsf{G}}_{\mathrm{f}}{\mathbf{F}}_{\mathrm{f}}^{-\mathrm{T}}+{\mathsf{G}}_{\mathrm{f}}^{\mathrm{T}}{\left({\nabla }_{\mathbf{X}}\mathbf{V}\right)}^{\mathrm{T}}{\mathbf{F}}_{\mathrm{f}}^{-\mathrm{T}})& \mathrm{in}{\Omega }_{\mathrm{f}}^0,\hfill \end{array}$$

(63)

$$\begin{array}{cc}\hfill {\mathbf{F}}_{\alpha }=\mathbf{I}+{\nabla }_{\mathbf{X}}{\mathbf{u}}_{\alpha }& \mathrm{in}{\Omega }_{\alpha }^0,\hfill \end{array}$$

(64)

$$\begin{array}{cc}\hfill {\mathbf{F}}_{\mathrm{II}}=\mathbf{I}+{\nabla }_{\mathbf{X}}{\mathbf{u}}_{\mathrm{II}}& \mathrm{in}{\Omega }_{\mathrm{II}}^0,\hfill \end{array}$$

(65)

where in (50)–(65), we drop the primes for the sake of a simpler notation, and we have the parameter:

$$ϵ=\frac{d}{L}.$$

(66)

Within the next section, we introduce the two-scale asymptotic homogenization technique, which we then use to upscale our system (50)–(65) in order to obtain our macroscale model.

3.2. The Two-Scale Asymptotic Homogenization Method

Here, we now introduce the asymptotic homogenization technique, which we will use to derive the effective macroscale governing equations for our structure. As described above, we have the pore-scale, denoted by d, and the macroscale, which is the average size of the whole material, denoted by L. We assume that these scales are well separated, i.e.,

$$ϵ=\frac{d}{L}\ll 1.$$

(67)

We can describe $ϵ$ as the scale separation parameter. We require a local scale spatial variable, which will capture the pore-scale variations of each of the fields appearing in (50)–(65), that is,

$$\overline{\mathbf{Y}}=\frac{\mathbf{X}}{ϵ}.$$

(68)

We also have the macroscale variable:

$$\overline{\mathbf{X}}=\mathbf{X}.$$

(69)

The newly introduced spatial variables $\overline{\mathbf{X}}$ and $\overline{\mathbf{Y}}$ represent the macroscale and the pore-scale respectively and are to be considered formally independent. The gradient operator also transforms, by the application of the chain rule, to become an operator of both scales, which we can write as:

$${\nabla }_{\mathbf{X}}\to {\nabla }_{\overline{\mathbf{X}}}+\frac{1}{ϵ}{\nabla }_{\overline{\mathbf{Y}}}.$$

(70)

We also make the assumption that all the fields in (50)–(65), $\forall \alpha =1,...,K$, are functions of the two spatial variables $\overline{\mathbf{X}}$ and $\overline{\mathbf{Y}}$ and that every field can be written as a power series expansion in $ϵ$:

$${\phi }^{ϵ}(\overline{\mathbf{X}},\overline{\mathbf{Y}},t)=\sum _{l=0}^{\infty }{\phi }^{\left(l\right)}(\overline{\mathbf{X}},\overline{\mathbf{Y}},t){ϵ}^l,$$

(71)

where $\phi$ is used to denote a general field in (50)–(65).

Remark 2

(Pore-scale periodicity). We assume that every field ${\phi }^{\left(l\right)}$ in our analysis is $\overline{\mathbf{Y}}$-periodic, and this allows us to focus our attention on a single periodic cell in our structure. This is a technical assumption that is related solely to the microscale. The fields can still vary with respect to the macroscale. Making this assumption means that we formulate the pore-scale differential problems on this periodic cell, and it is these problems that are to be solved to determine the macroscale model coefficients. In [40], the authors assumed that all the fields are periodic in the pore-scale variable and are able to compute the coefficients of the model of standard poroelasticity obtained via asymptotic homogenization. In comparison, in [41], the macroscale model of poroelasticity was solved. As such, although the effective coefficients used are those computed by following the methodology described in [40], the illustrated solution solely depends on the macroscale and is in general not periodic. It is possible to relax this assumption that all fields are $\overline{\mathbf{Y}}$-periodic and assume instead that all the fields are locally bounded. This assumption is less strict than local periodicity and means that all the fields are finite with respect to the pore-scale variable $\overline{\mathbf{Y}}$ when $ϵ\to 0$, but not necessarily periodic. This assumption however only permits the functional form of the macroscale model to be derived and does not in general allow for the model coefficients to be computed without further assumptions (see [21,23] for further details) in multidimensional problems [23].

We have that our periodic cell could potentially contain a variety of subphases and these subphases could possess various geometries and elastic properties. This assumption is particularly useful as it allows us to solve the differential problems obtained from the asymptotic homogenization technique on a single periodic cell instead of the whole material domain, therefore reducing the computational complexity. This periodic cell where we solve the differential problems is shown in Figure 2.

Remark 3

(Macroscopic uniformity). We know that the pore-scale structure of a material can vary with respect to the macroscale position (see [21,25,39,42,43]). In general, this dependence is neglected in the literature due to wanting to simplify the analysis. Here, we assume that the pore-scale geometry does not depend on the macroscale variable $\overline{\mathbf{X}}$, i.e., the material is macroscopically uniform. This assumption allows for simple differentiation under the integral sign to take place, that is:

$${\int }_{\Omega }{\nabla }_{\overline{\mathbf{X}}}·(•)={\nabla }_{\overline{\mathbf{X}}}·{\int }_{\Omega }(•).$$

(72)

If we do not assume macroscopic uniformity, then (72) is not satisfied, and in this case, the application of the Reynolds’ transport theorem is required. This may lead to additional terms appearing in the macroscale governing equations.

Remark 4

(Pore-scale geometry). In our description so far, we assume that there are various different subphases included within our periodic cell; however, without loss of generality, we can focus our attention on the case where each periodic cell contains only one hyperelastic subphase. This structure is shown in Figure 3. Therefore, we do not require the index α, and the notation can be adjusted as follows. We have the continuum body ${\mathcal{B}}_0$, which has a periodic microstructure. Within ${\mathcal{B}}_0$, we have many periodic cells; however, due to periodicity, we can identify the domain ${\Omega }^0$ with the periodic cell, which has a hyperelastic subphase, hyperelastic matrix and fluid portions denoted by ${\Omega }_{\mathrm{I}}^0$, ${\Omega }_{\mathrm{II}}^0$ and ${\Omega }_{\mathrm{f}}^0$, respectively. We are also able to simplify the notation we use for the different interfaces. The interface between the subphase and the fluid is ${\Gamma }_{\mathrm{I}}^0:=\partial {\Omega }_{\mathrm{I}}^0\cap \partial {\Omega }_{\mathrm{f}}^0$; the interface between the matrix and the fluid is ${\Gamma }_{\mathrm{II}}^0:=\partial {\Omega }_{\mathrm{II}}^0\cap \partial {\Omega }_{\mathrm{f}}^0$; the interface between our two hyperelastic solid phases is ${\Gamma }_{\mathrm{III}}^0:=\partial {\Omega }_{\mathrm{I}}^0\cap \partial {\Omega }_{\mathrm{II}}^0$, with corresponding unit normal vectors ${\mathbf{n}}_{\mathrm{I}}^0$, ${\mathbf{n}}_{\mathrm{II}}^0$ and ${\mathbf{n}}_{\mathrm{III}}^0$. If a specific application required multiple subphases to be contained in the periodic cell, then it would be simple to extend the formulation, as has been done in the case of elastic composites [44].

3.3. The Macroscale Results

The assumptions (70) and (71) of the asymptotic homogenization technique can be applied to Equations (50)–(65). We then obtain the following multiscale system of PDEs:

$$\begin{array}{cc}\hfill {\nabla }_{\overline{\mathbf{Y}}}·{\sigma }_{\mathrm{I}}^{ϵ}+ϵ{\nabla }_{\overline{\mathbf{X}}}·{\sigma }_{\mathrm{I}}^{ϵ}=0& \mathrm{in}{\Omega }_{\mathrm{I}}^0,\hfill \end{array}$$

(73)

$$\begin{array}{cc}\hfill {\nabla }_{\overline{\mathbf{Y}}}·{\sigma }_{\mathrm{II}}^{ϵ}+ϵ{\nabla }_{\overline{\mathbf{X}}}·{\sigma }_{\mathrm{II}}^{ϵ}=0& \mathrm{in}{\Omega }_{\mathrm{II}}^0,\hfill \end{array}$$

(74)

$$\begin{array}{cc}\hfill {\nabla }_{\overline{\mathbf{Y}}}·{\left({{\sigma }_{\mathrm{f}}}^{ϵ}({\mathsf{G}}_{\mathrm{f}}^{ϵ}\right)}^{\mathrm{T}})+ϵ{\nabla }_{\overline{\mathbf{X}}}·{\left({{\sigma }_{\mathrm{f}}}^{ϵ}({\mathsf{G}}_{\mathrm{f}}^{ϵ}\right)}^{\mathrm{T}})=0& \mathrm{in}{\Omega }_{\mathrm{f}}^0,\hfill \end{array}$$

(75)

$$\begin{array}{cc}\hfill {\nabla }_{\overline{\mathbf{Y}}}·\left({\mathsf{G}}_{\mathrm{f}}^{ϵ}{\mathbf{V}}^{ϵ}\right)+ϵ{\nabla }_{\overline{\mathbf{X}}}·\left({\mathsf{G}}_{\mathrm{f}}^{ϵ}{\mathbf{V}}^{ϵ}\right)=0& \mathrm{in}{\Omega }_{\mathrm{f}}^0,\hfill \end{array}$$

(76)

$$\begin{array}{cc}\hfill \frac{\partial {\mathbf{u}}_{\mathrm{I}}^{ϵ}}{\partial t}={\mathbf{V}}^{ϵ}& \mathrm{on}{\Gamma }_{\mathrm{I}}^0,\hfill \end{array}$$

(77)

$$\begin{array}{cc}\hfill \frac{\partial {\mathbf{u}}_{\mathrm{II}}^{ϵ}}{\partial t}={\mathbf{V}}^{ϵ}& \mathrm{on}{\Gamma }_{\mathrm{II}}^0,\hfill \end{array}$$

(78)

$$\begin{array}{cc}\hfill {{\sigma }_{\mathrm{f}}}^{ϵ}{\left({\mathsf{G}}_{\mathrm{f}}^{ϵ}\right)}^{\mathrm{T}}{\mathbf{n}}_{\mathrm{I}}^0={\sigma }_{\mathrm{I}}^{ϵ}{\mathbf{n}}_{\mathrm{I}}^0& \mathrm{on}{\Gamma }_{\mathrm{I}}^0,\hfill \end{array}$$

(79)

$$\begin{array}{cc}\hfill {{\sigma }_{\mathrm{f}}}^{ϵ}{\left({\mathsf{G}}_{\mathrm{f}}^{ϵ}\right)}^{\mathrm{T}}{\mathbf{n}}_{\mathrm{II}}^0={\sigma }_{\mathrm{II}}^{ϵ}{\mathbf{n}}_{\mathrm{II}}^0& \mathrm{on}{\Gamma }_{\mathrm{II}}^0,\hfill \end{array}$$

(80)

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{I}}^{ϵ}{\mathbf{n}}_{\mathrm{III}}^0={\sigma }_{\mathrm{II}}^{ϵ}{\mathbf{n}}_{\mathrm{III}}^0& \mathrm{on}{\Gamma }_{\mathrm{III}}^0,\hfill \end{array}$$

(81)

$$\begin{array}{cc}\hfill {\mathbf{u}}_{\mathrm{I}}^{ϵ}={\mathbf{u}}_{\mathrm{II}}^{ϵ}& \mathrm{on}{\Gamma }_{\mathrm{III}}^0,\hfill \end{array}$$

(82)

as well as the multiscale constitutive equations for ${{\sigma }_{\mathrm{f}}}^{ϵ}{\left({\mathsf{G}}_{\mathrm{f}}^{ϵ}\right)}^{\mathrm{T}}$, ${\sigma }_{\mathrm{I}}^{ϵ}$ and ${\sigma }_{\mathrm{II}}^{ϵ}$, which are given by:

$$\begin{array}{cc}\hfill {{\sigma }_{\mathrm{f}}}^{ϵ}{\left({\mathsf{G}}_{\mathrm{f}}^{ϵ}\right)}^{\mathrm{T}}=& -p^{ϵ}{\left({\mathsf{G}}_{\mathrm{f}}^{ϵ}\right)}^{\mathrm{T}}+ϵ(\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{V}}^{ϵ}\right){\left({\mathbf{F}}_{\mathrm{f}}^{ϵ}\right)}^{-1}{\left({\mathsf{G}}_{\mathrm{f}}^{ϵ}\right)}^{\mathrm{T}}+{\left({\mathbf{F}}_{\mathrm{f}}^{ϵ}\right)}^{-\mathrm{T}}{\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{V}}^{ϵ}\right)}^{\mathrm{T}}{\left({\mathsf{G}}_{\mathrm{f}}^{ϵ}\right)}^{\mathrm{T}})\hfill \end{array}$$

$$\begin{array}{cc}\hfill & +{ϵ}^2(\left({\nabla }_{\overline{\mathbf{X}}}{\mathbf{V}}^{ϵ}\right){\left({\mathbf{F}}_{\mathrm{f}}^{ϵ}\right)}^{-1}{\left({\mathsf{G}}_{\mathrm{f}}^{ϵ}\right)}^{\mathrm{T}}+{\left({\mathbf{F}}_{\mathrm{f}}^{ϵ}\right)}^{-\mathrm{T}}{\left({\nabla }_{\overline{\mathbf{X}}}{\mathbf{V}}^{ϵ}\right)}^{\mathrm{T}}{\left({\mathsf{G}}_{\mathrm{f}}^{ϵ}\right)}^{\mathrm{T}})\mathrm{in}{\Omega }_{\mathrm{f}}^0,\hfill \end{array}$$

(83)

$$\begin{array}{cc}\hfill & {\sigma }_{\mathrm{I}}^{ϵ}=\frac{\partial {\psi }_{\mathrm{I}}^{ϵ}}{\partial {\mathbf{F}}_{\mathrm{I}}^{ϵ}}\mathrm{in}{\Omega }_{\mathrm{I}}^0,\hfill \end{array}$$

(84)

$$\begin{array}{cc}\hfill & {\sigma }_{\mathrm{II}}^{ϵ}=\frac{\partial {\psi }_{\mathrm{II}}^{ϵ}}{\partial {\mathbf{F}}_{\mathrm{II}}^{ϵ}}\mathrm{in}{\Omega }_{\mathrm{II}}^0,\hfill \end{array}$$

(85)

and the fluid balance equation and the solid deformation gradients are given by:

$$\begin{array}{cc}\hfill 0=-{\left({\mathsf{G}}_{\mathrm{f}}^{ϵ}\right)}^{\mathrm{T}}(ϵ{\nabla }_{\overline{\mathbf{X}}}{p}^{ϵ}& +{\nabla }_{\overline{\mathbf{Y}}}{p}^{ϵ})+ϵ{\nabla }_{\overline{\mathbf{Y}}}·(\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{V}}^{ϵ}\right){\mathsf{G}}_{\mathrm{f}}^{ϵ}{\left({\mathbf{F}}_{\mathrm{f}}^{ϵ}\right)}^{-\mathrm{T}}\hfill \end{array}$$

$$\begin{array}{cc}\hfill & +{\left({\mathsf{G}}_{\mathrm{f}}^{ϵ}\right)}^{\mathrm{T}}{\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{V}}^{ϵ}\right)}^{\mathrm{T}}{\left({\mathbf{F}}_{\mathrm{f}}^{ϵ}\right)}^{-\mathrm{T}})+O\left({ϵ}^2\right)\mathrm{in}{\Omega }_{\mathrm{f}}^0,\hfill \end{array}$$

(86)

$$\begin{array}{cc}\hfill & ϵ{\mathbf{F}}_{\mathrm{I}}^{ϵ}=ϵ\mathbf{I}+ϵ{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_{\mathrm{I}}^{ϵ}+{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_{\mathrm{I}}^{ϵ}\mathrm{in}{\Omega }_{\mathrm{I}}^0,\hfill \end{array}$$

(87)

$$\begin{array}{cc}\hfill & ϵ{\mathbf{F}}_{\mathrm{II}}^{ϵ}=ϵ\mathbf{I}+ϵ{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_{\mathrm{II}}^{ϵ}+{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_{\mathrm{II}}^{ϵ}\mathrm{in}{\Omega }_{\mathrm{II}}^0,\hfill \end{array}$$

(88)

where the power series representation (71) is implied in Equations (73)–(88) through the use of the superscript $ϵ$. We then proceed with the technique by equating the coefficients of ${ϵ}^l$ for $l=0,1,...$, and this way, we derive the effective macroscale model in terms of the relevant zeroth-order fields. Following the asymptotic expansion, if any term still retains a dependence on the pore-scale, then we can apply the integral average formula. The integral average can be defined as:

$${\langle \phi \rangle }_i=\frac{1}{|{\Omega }^0|}{\int }_{{\Omega }_i^0}\phi (\overline{\mathbf{X}},\overline{\mathbf{Y}},t)\mathrm{d}\overline{\mathbf{Y}}i=\mathrm{f},\mathrm{I},\mathrm{II},$$

(89)

and again, $\phi$ is a general field. The integral average is performed over one representative cell due to the assumption of $\overline{\mathbf{Y}}$-periodicity, as discussed in Remark 2, so it is therefore a cell average. We have that the volume of the domain is given by $|{\Omega }^0|=|{\Omega }_{\mathrm{f}}^0|+|{\Omega }_{\mathrm{I}}^0|+|{\Omega }_{\mathrm{II}}^0|$.

We can equate the coefficients of ${ϵ}^0$ in (73)–(82) to obtain:

$$\begin{array}{cc}\hfill {\nabla }_{\overline{\mathbf{Y}}}·{\sigma }_{\mathrm{I}}^{\left(0\right)}=0& \mathrm{in}{\Omega }_{\mathrm{I}}^0,\hfill \end{array}$$

(90)

$$\begin{array}{cc}\hfill {\nabla }_{\overline{\mathbf{Y}}}·{\sigma }_{\mathrm{II}}^{\left(0\right)}=0& \mathrm{in}{\Omega }_{\mathrm{II}}^0,\hfill \end{array}$$

(91)

$$\begin{array}{cc}\hfill {\nabla }_{\overline{\mathbf{Y}}}·\left({{\sigma }_{\mathrm{f}}}^{\left(0\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}\right)=0& \mathrm{in}{\Omega }_{\mathrm{f}}^0,\hfill \end{array}$$

(92)

$$\begin{array}{cc}\hfill {\nabla }_{\overline{\mathbf{Y}}}·\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}{\mathbf{V}}^{\left(0\right)}\right)=0& \mathrm{in}{\Omega }_{\mathrm{f}}^0,\hfill \end{array}$$

(93)

$$\begin{array}{cc}\hfill \frac{\partial {\mathbf{u}}_{\mathrm{I}}^{\left(0\right)}}{\partial t}={\mathbf{V}}^{\left(0\right)}& \mathrm{on}{\Gamma }_{\mathrm{I}}^0,\hfill \end{array}$$

(94)

$$\begin{array}{cc}\hfill \frac{\partial {\mathbf{u}}_{\mathrm{II}}^{\left(0\right)}}{\partial t}={\mathbf{V}}^{\left(0\right)}& \mathrm{on}{\Gamma }_{\mathrm{II}}^0,\hfill \end{array}$$

(95)

$$\begin{array}{cc}\hfill {{\sigma }_{\mathrm{f}}}^{\left(0\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}{\mathbf{n}}_{\mathrm{I}}^0={\sigma }_{\mathrm{I}}^{\left(0\right)}{\mathbf{n}}_{\mathrm{I}}^0& \mathrm{on}{\Gamma }_{\mathrm{I}}^0,\hfill \end{array}$$

(96)

$$\begin{array}{cc}\hfill {{\sigma }_{\mathrm{f}}}^{\left(0\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}{\mathbf{n}}_{\mathrm{II}}^0={\sigma }_{\mathrm{II}}^{\left(0\right)}{\mathbf{n}}_{\mathrm{II}}^0& \mathrm{on}{\Gamma }_{\mathrm{II}}^0,\hfill \end{array}$$

(97)

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{I}}^{\left(0\right)}{\mathbf{n}}_{\mathrm{III}}^0={\sigma }_{\mathrm{II}}^{\left(0\right)}{\mathbf{n}}_{\mathrm{III}}^0& \mathrm{on}{\Gamma }_{\mathrm{III}}^0,\hfill \end{array}$$

(98)

$$\begin{array}{cc}\hfill {\mathbf{u}}_{\mathrm{I}}^{\left(0\right)}={\mathbf{u}}_{\mathrm{II}}^{\left(0\right)}& \mathrm{on}{\Gamma }_{\mathrm{III}}^0.\hfill \end{array}$$

(99)

The constitutive Equations (83)–(85) have coefficients of ${ϵ}^0$ given by:

$$\begin{array}{cc}\hfill {{\sigma }_{\mathrm{f}}}^{\left(0\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}=-p^{\left(0\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}& \mathrm{in}{\Omega }_{\mathrm{f}}^0,\hfill \end{array}$$

(100)

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{I}}^{\left(0\right)}=\frac{\partial {\psi }_{\mathrm{I}}^{\left(0\right)}}{\partial {\mathbf{F}}_{\mathrm{I}}^{\left(0\right)}}& \mathrm{in}{\Omega }_{\mathrm{I}}^0,\hfill \end{array}$$

(101)

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{II}}^{\left(0\right)}=\frac{\partial {\psi }_{\mathrm{II}}^{\left(0\right)}}{\partial {\mathbf{F}}_{\mathrm{II}}^{\left(0\right)}}& \mathrm{in}{\Omega }_{\mathrm{II}}^0,\hfill \end{array}$$

(102)

and the fluid balance equation and the deformation gradients (86)–(88) have coefficients of ${ϵ}^0$:

$$\begin{array}{cc}\hfill {\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{Y}}}{p}^{\left(0\right)}=0& \mathrm{in}{\Omega }_{\mathrm{f}}^0,\hfill \end{array}$$

(103)

$$\begin{array}{cc}\hfill {\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_{\mathrm{I}}^{\left(0\right)}=0& \mathrm{in}{\Omega }_{\mathrm{I}}^0,\hfill \end{array}$$

(104)

$$\begin{array}{cc}\hfill {\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_{\mathrm{II}}^{\left(0\right)}=0& \mathrm{in}{\Omega }_{\mathrm{II}}^0,\hfill \end{array}$$

(105)

Now, equating the coefficients of ${ϵ}^1$ in Equations (73)–(82) gives:

$$\begin{array}{cc}\hfill {\nabla }_{\overline{\mathbf{Y}}}·{\sigma }_{\mathrm{I}}^{\left(1\right)}+{\nabla }_{\overline{\mathbf{X}}}·{\sigma }_{\mathrm{I}}^{\left(0\right)}=0& \mathrm{in}{\Omega }_{\mathrm{I}}^0,\hfill \end{array}$$

(106)

$$\begin{array}{cc}\hfill {\nabla }_{\overline{\mathbf{Y}}}·{\sigma }_{\mathrm{II}}^{\left(1\right)}+{\nabla }_{\overline{\mathbf{X}}}·{\sigma }_{\mathrm{II}}^{\left(0\right)}=0& \mathrm{in}{\Omega }_{\mathrm{II}}^0,\hfill \end{array}$$

(107)

$$\begin{array}{cc}\hfill {\nabla }_{\overline{\mathbf{Y}}}·\left({{\sigma }_{\mathrm{f}}}^{\left(1\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}\right)+{\nabla }_{\overline{\mathbf{Y}}}·\left({{\sigma }_{\mathrm{f}}}^{\left(0\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(1\right)}\right)}^{\mathrm{T}}\right)+{\nabla }_{\overline{\mathbf{X}}}·\left({{\sigma }_{\mathrm{f}}}^{\left(0\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}\right)=0& \mathrm{in}{\Omega }_{\mathrm{f}}^0,\hfill \end{array}$$

(108)

$$\begin{array}{cc}\hfill {\nabla }_{\overline{\mathbf{Y}}}·\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}{\mathbf{V}}^{\left(1\right)}\right)+{\nabla }_{\overline{\mathbf{Y}}}·\left({\mathsf{G}}_{\mathrm{f}}^{\left(1\right)}{\mathbf{v}}^{\left(0\right)}\right)+{\nabla }_{\overline{\mathbf{X}}}·\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}{\mathbf{V}}^{\left(0\right)}\right)=0& \mathrm{in}{\Omega }_{\mathrm{f}}^0,\hfill \end{array}$$

(109)

$$\begin{array}{cc}\hfill \frac{\partial {\mathbf{u}}_{\mathrm{I}}^{\left(1\right)}}{\partial t}={\mathbf{V}}^{\left(1\right)}& \mathrm{on}{\Gamma }_{\mathrm{I}}^0,\hfill \end{array}$$

(110)

$$\begin{array}{cc}\hfill \frac{\partial {\mathbf{u}}_{\mathrm{II}}^{\left(1\right)}}{\partial t}={\mathbf{V}}^{\left(1\right)}& \mathrm{on}{\Gamma }_{\mathrm{II}}^0,\hfill \end{array}$$

(111)

$$\begin{array}{cc}\hfill ({{\sigma }_{\mathrm{f}}}^{\left(1\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}+{{\sigma }_{\mathrm{f}}}^{\left(0\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(1\right)}\right)}^{\mathrm{T}}){\mathbf{n}}_{\mathrm{I}}^0={\sigma }_{\mathrm{I}}^{\left(1\right)}{\mathbf{n}}_{\mathrm{I}}^0& \mathrm{on}{\Gamma }_{\mathrm{I}}^0,\hfill \end{array}$$

(112)

$$\begin{array}{cc}\hfill ({{\sigma }_{\mathrm{f}}}^{\left(1\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}+{{\sigma }_{\mathrm{f}}}^{\left(0\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(1\right)}\right)}^{\mathrm{T}}){\mathbf{n}}_{\mathrm{II}}^0={\sigma }_{\mathrm{II}}^{\left(1\right)}{\mathbf{n}}_{\mathrm{II}}^0& \mathrm{on}{\Gamma }_{\mathrm{II}}^0,\hfill \end{array}$$

(113)

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{I}}^{\left(1\right)}{\mathbf{n}}_{\mathrm{III}}^0={\sigma }_{\mathrm{II}}^{\left(1\right)}{\mathbf{n}}_{\mathrm{III}}^0& \mathrm{on}{\Gamma }_{\mathrm{III}}^0,\hfill \end{array}$$

(114)

$$\begin{array}{cc}\hfill {\mathbf{u}}_{\mathrm{I}}^{\left(1\right)}={\mathbf{u}}_{\mathrm{II}}^{\left(1\right)}& \mathrm{on}{\Gamma }_{\mathrm{III}}^0,\hfill \end{array}$$

(115)

and the coefficients of ${ϵ}^1$ in the constitutive Equations (83)–(85) are:

$$\begin{array}{cc}\hfill {{\sigma }_{\mathrm{f}}}^{\left(1\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}+{{\sigma }_{\mathrm{f}}}^{\left(0\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(1\right)}\right)}^{\mathrm{T}}=-p^{\left(1\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}-p^{\left(0\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(1\right)}\right)}^{\mathrm{T}}& \hfill \end{array}$$

$$\begin{array}{cc}\hfill +(\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{V}}^{\left(0\right)}\right){\left({\mathbf{F}}_{\mathrm{f}}^{\left(0\right)}\right)}^{-1}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}+{\left({\mathbf{F}}_{\mathrm{f}}^{\left(0\right)}\right)}^{-\mathrm{T}}{\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{V}}^{\left(0\right)}\right)}^{\mathrm{T}}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}})& \mathrm{in}{\Omega }_{\mathrm{f}}^0,\hfill \end{array}$$

(116)

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{I}}^{\left(1\right)}=\frac{\partial {\psi }_{\mathrm{I}}^{\left(1\right)}}{\partial {\mathbf{F}}_{\mathrm{I}}^{\left(1\right)}}& \mathrm{in}{\Omega }_{\mathrm{I}}^0,\hfill \end{array}$$

(117)

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{II}}^{\left(1\right)}=\frac{\partial {\psi }_{\mathrm{II}}^{\left(1\right)}}{\partial {\mathbf{F}}_{\mathrm{II}}^{\left(1\right)}}& \mathrm{in}{\Omega }_{\mathrm{II}}^0,\hfill \end{array}$$

(118)

while the fluid balance equation and the deformation gradients (86)–(88) have coefficients of ${ϵ}^1$:

$$\begin{array}{cc}\hfill 0=-{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}({\nabla }_{\overline{\mathbf{X}}}{p}^{\left(0\right)}+{\nabla }_{\overline{\mathbf{Y}}}{p}^{\left(1\right)})+{\nabla }_{\overline{\mathbf{Y}}}·(\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{V}}^{\left(0\right)}\right){\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}{\left({\mathbf{F}}_{\mathrm{f}}^{\left(0\right)}\right)}^{-\mathrm{T}}& \hfill \end{array}$$

$$\begin{array}{cc}\hfill +{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}{\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{V}}^{\left(0\right)}\right)}^{\mathrm{T}}{\left({\mathbf{F}}_{\mathrm{f}}^{\left(0\right)}\right)}^{-\mathrm{T}})& \mathrm{in}{\Omega }_{\mathrm{f}}^0,\hfill \end{array}$$

(119)

$$\begin{array}{cc}\hfill {\mathbf{F}}_{\mathrm{I}}^{\left(0\right)}=\mathbf{I}+{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_{\mathrm{I}}^{\left(0\right)}+{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_{\mathrm{I}}^{\left(1\right)}& \mathrm{in}{\Omega }_{\mathrm{I}}^0,\hfill \end{array}$$

(120)

$$\begin{array}{cc}\hfill {\mathbf{F}}_{\mathrm{II}}^{\left(0\right)}=\mathbf{I}+{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_{\mathrm{II}}^{\left(0\right)}+{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_{\mathrm{II}}^{\left(1\right)}& \mathrm{in}{\Omega }_{\mathrm{II}}^0.\hfill \end{array}$$

(121)

We can show that the Piola transformation ${\left({\mathsf{G}}^{ϵ}\right)}^{\mathrm{T}}$, where we exploit the notation (16), is divergence free and derive some useful identities. We have:

$$\begin{array}{cc}\hfill {\int }_{{\mathcal{B}}_0}{\nabla }_{\mathbf{X}}·{\mathsf{G}}^{\mathrm{T}}{\mathrm{dV}}_{\mathbf{X}}& ={\int }_{\partial {\mathcal{B}}_0}{\mathsf{G}}^{\mathrm{T}}·\mathrm{NdA}={\int }_{\partial {\mathcal{B}}_0}{\mathrm{JF}}^{-\mathrm{T}}\mathrm{NdA}={\int }_{\partial {\mathcal{B}}_t}\mathbf{I}·\mathrm{nda}\hfill \\ \hfill & ={\int }_{{\mathcal{B}}_t}{\nabla }_{\mathbf{x}}·\mathbf{I}{\mathrm{dV}}_{\mathbf{x}}=0,\hfill \end{array}$$

(122)

and so, ${\left({\mathsf{G}}^{ϵ}\right)}^{\mathrm{T}}$ is divergence free. We can also consider the expansion of ${\left({\mathsf{G}}^{ϵ}\right)}^{\mathrm{T}}$, which is:

$${\left({\mathsf{G}}^{ϵ}\right)}^{\mathrm{T}}={\left({\mathsf{G}}^{\left(0\right)}\right)}^{\mathrm{T}}+ϵ{\left({\mathsf{G}}^{\left(1\right)}\right)}^{\mathrm{T}}+o\left({ϵ}^2\right),$$

(123)

and the expansion of ${\nabla }_X·{\left({\mathsf{G}}^{ϵ}\right)}^{\mathrm{T}}$ is:

$$(ϵ{\nabla }_{\overline{\mathbf{X}}}+{\nabla }_{\overline{\mathbf{Y}}})·({\left({\mathsf{G}}^{\left(0\right)}\right)}^{\mathrm{T}}+ϵ{\left({\mathsf{G}}^{\left(1\right)}\right)}^{\mathrm{T}}+o\left({ϵ}^2\right))=0.$$

(124)

Then, equating the coefficient of ${ϵ}^0$ gives:

$${\nabla }_{\overline{\mathbf{Y}}}·{\left({\mathsf{G}}^{\left(0\right)}\right)}^{\mathrm{T}}=0,$$

(125)

and equating the coefficient of ${ϵ}^1$ gives:

$${\nabla }_{\overline{\mathbf{X}}}·{\left({\mathsf{G}}^{\left(0\right)}\right)}^{\mathrm{T}}+{\nabla }_{\overline{\mathbf{Y}}}·{\left({\mathsf{G}}^{\left(1\right)}\right)}^{\mathrm{T}}=0\Rightarrow {\nabla }_{\overline{\mathbf{X}}}·{\left({\mathsf{G}}^{\left(0\right)}\right)}^{\mathrm{T}}=-{\nabla }_{\overline{\mathbf{Y}}}·{\left({\mathsf{G}}^{\left(1\right)}\right)}^{\mathrm{T}}.$$

(126)

We can use the notation (16) to write (125) and (126) as their counterparts in each of the solid domains and the fluid domain as:

$${\nabla }_{\overline{\mathbf{Y}}}·{\left({\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}\right)}^{\mathrm{T}}=0,{\nabla }_{\overline{\mathbf{Y}}}·{\left({\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}\right)}^{\mathrm{T}}=0,{\nabla }_{\overline{\mathbf{Y}}}·{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}=0,$$

(127)

and:

$$\begin{array}{c}{\nabla }_{\overline{\mathbf{X}}}·{\left({\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}\right)}^{\mathrm{T}}=-{\nabla }_{\overline{\mathbf{Y}}}·{\left({\mathsf{G}}_{\mathrm{I}}^{\left(1\right)}\right)}^{\mathrm{T}},{\nabla }_{\overline{\mathbf{X}}}·{\left({\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}\right)}^{\mathrm{T}}=-{\nabla }_{\overline{\mathbf{Y}}}·{\left({\mathsf{G}}_{\mathrm{II}}^{\left(1\right)}\right)}^{\mathrm{T}},\\ {\nabla }_{\overline{\mathbf{X}}}·{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}=-{\nabla }_{\overline{\mathbf{Y}}}·{\left({\mathsf{G}}_{\mathrm{f}}^{\left(1\right)}\right)}^{\mathrm{T}}.\end{array}$$

(128)

Here, we also wish to perform the multiscale expansion of (19). We have:

$$\begin{array}{c}{\left({\mathsf{G}}_{\mathrm{I}}^{ϵ}\right)}^{\mathrm{T}}{\mathbf{n}}_{\mathrm{I}}^0={\left({\mathsf{G}}_{\mathrm{f}}^{ϵ}\right)}^{\mathrm{T}}{\mathbf{n}}_{\mathrm{I}}^0\mathrm{on}{\Gamma }_{\mathrm{I}},{\left({\mathsf{G}}_{\mathrm{II}}^{ϵ}\right)}^{\mathrm{T}}{\mathbf{n}}_{\mathrm{II}}^0={\left({\mathsf{G}}_{\mathrm{f}}^{ϵ}\right)}^{\mathrm{T}}{\mathbf{n}}_{\mathrm{II}}^0\mathrm{on}{\Gamma }_{\mathrm{II}},\\ {\left({\mathsf{G}}_{\mathrm{I}}^{ϵ}\right)}^{\mathrm{T}}{\mathbf{n}}_{\mathrm{III}}^0={\left({\mathsf{G}}_{\mathrm{II}}^{ϵ}\right)}^{\mathrm{T}}{\mathbf{n}}_{\mathrm{III}}^0\mathrm{on}{\Gamma }_{\mathrm{III}},\end{array}$$

(129)

and equating the coefficients of ${ϵ}^0$ gives:

$$\begin{array}{c}{\left({\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}\right)}^{\mathrm{T}}{\mathbf{n}}_{\mathrm{I}}^0={\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}{\mathbf{n}}_{\mathrm{I}}^0\mathrm{on}{\Gamma }_{\mathrm{I}},{\left({\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}\right)}^{\mathrm{T}}{\mathbf{n}}_{\mathrm{II}}^0={\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}{\mathbf{n}}_{\mathrm{II}}^0\mathrm{on}{\Gamma }_{\mathrm{II}},\\ {\left({\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}\right)}^{\mathrm{T}}{\mathbf{n}}_{\mathrm{III}}^0={\left({\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}\right)}^{\mathrm{T}}{\mathbf{n}}_{\mathrm{III}}^0\mathrm{on}{\Gamma }_{\mathrm{III}}.\end{array}$$

(130)

Similarly, equating the coefficients of ${ϵ}^1$ gives:

$$\begin{array}{c}{\left({\mathsf{G}}_{\mathrm{I}}^{\left(1\right)}\right)}^{\mathrm{T}}{\mathbf{n}}_{\mathrm{I}}^0={\left({\mathsf{G}}_{\mathrm{f}}^{\left(1\right)}\right)}^{\mathrm{T}}{\mathbf{n}}_{\mathrm{I}}^0\mathrm{on}{\Gamma }_{\mathrm{I}},{\left({\mathsf{G}}_{\mathrm{II}}^{\left(1\right)}\right)}^{\mathrm{T}}{\mathbf{n}}_{\mathrm{II}}^0={\left({\mathsf{G}}_{\mathrm{f}}^{\left(1\right)}\right)}^{\mathrm{T}}{\mathbf{n}}_{\mathrm{II}}^0\mathrm{on}{\Gamma }_{\mathrm{II}},\\ {\left({\mathsf{G}}_{\mathrm{I}}^{\left(1\right)}\right)}^{\mathrm{T}}{\mathbf{n}}_{\mathrm{III}}^0={\left({\mathsf{G}}_{\mathrm{II}}^{\left(1\right)}\right)}^{\mathrm{T}}{\mathbf{n}}_{\mathrm{III}}^0\mathrm{on}{\Gamma }_{\mathrm{III}}.\end{array}$$

(131)

Equating higher powers of epsilon leads to the continuity at higher orders also. We use these expressions in later sections.

We can see using (92), (100) and (125) that ${\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{Y}}}{p}^{\left(0\right)}=0$. This implies that $p^{\left(0\right)}$ does not depend on the pore-scale $\overline{\mathbf{Y}}$, that is:

$$p^{\left(0\right)}=p^{\left(0\right)}(\overline{\mathbf{X}},t)$$

(132)

We can deduce from (104) and (105) that ${\mathbf{u}}_{\mathrm{I}}^{\left(0\right)}$ and ${\mathbf{u}}_{\mathrm{II}}^{\left(0\right)}$ do not depend on the pore-scale $\overline{\mathbf{Y}}$, that is,

$$\begin{array}{cc}\hfill & {\mathbf{u}}_{\mathrm{I}}^{\left(0\right)}={\mathbf{u}}_{\mathrm{I}}^{\left(0\right)}(\overline{\mathbf{X}},t),\hfill \end{array}$$

(133)

$$\begin{array}{cc}\hfill & {\mathbf{u}}_{\mathrm{II}}^{\left(0\right)}={\mathbf{u}}_{\mathrm{II}}^{\left(0\right)}(\overline{\mathbf{X}},t),\hfill \end{array}$$

(134)

and since we have the boundary condition (99), we can write:

$${\mathbf{u}}^{\left(0\right)}={\mathbf{u}}_{\mathrm{I}}^{\left(0\right)}={\mathbf{u}}_{\mathrm{II}}^{\left(0\right)}.$$

(135)

We use (132) and (135) throughout the remainder of this work.

3.4. The Macroscale Fluid Flow

We can investigate the leading order velocity, which we denote by ${\mathbf{v}}^{\left(0\right)}$. We begin by defining the relative fluid–solid displacement, $\mathbf{w}$, as:

$$\mathbf{w}:={\mathbf{V}}^{\left(0\right)}-\frac{\partial {\mathbf{u}}^{\left(0\right)}}{\partial t}.$$

(136)

We can rearrange (136) to obtain:

$${\mathbf{V}}^{\left(0\right)}=\mathbf{w}+{\dot{\mathbf{u}}}^{\left(0\right)}$$

(137)

We are then able to use (137) and Equations (93)–(95), (100), (108) and (116) to form a Stokes’-type boundary value problem given by:

$$\begin{array}{cc}\hfill -{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}({\nabla }_{\overline{\mathbf{X}}}{p}^{\left(0\right)}+{\nabla }_{\overline{\mathbf{Y}}}{p}^{\left(1\right)})+{\nabla }_{\overline{\mathbf{Y}}}·(\left({\nabla }_{\overline{\mathbf{Y}}}\mathbf{w}\right){\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}{\left({\mathbf{F}}_{\mathrm{f}}^{\left(0\right)}\right)}^{-}\mathrm{T}& \hfill \end{array}$$

$$\begin{array}{cc}\hfill +{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}{\left({\nabla }_{\overline{\mathbf{Y}}}\mathbf{w}\right)}^{\mathrm{T}}{\left({\mathbf{F}}_{\mathrm{f}}^{\left(0\right)}\right)}^{-}\mathrm{T})=0& \mathrm{in}{\Omega }_{\mathrm{f}}^0,\hfill \end{array}$$

(138)

$$\begin{array}{cc}\hfill {\nabla }_{\overline{\mathbf{Y}}}·\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\mathbf{w}\right)=0& \mathrm{in}{\Omega }_{\mathrm{f}}^0,\hfill \end{array}$$

(139)

$$\begin{array}{cc}\hfill \mathbf{w}=\mathbf{0}& \mathrm{on}{\Gamma }_{\mathrm{I}}^0\cup {\Gamma }_{\mathrm{II}}^0.\hfill \end{array}$$

(140)

The boundary value problem (138)–(140) admits a solution. Exploiting linearity, the solution is given by:

$$\begin{array}{cc}\hfill & \mathbf{w}=-\mathbf{W}{\nabla }_{\overline{\mathbf{X}}}{p}^{\left(0\right)},\hfill \end{array}$$

(141)

$$\begin{array}{cc}\hfill & p^{\left(1\right)}=-\mathbf{P}{\nabla }_{\overline{\mathbf{X}}}{p}^{\left(0\right)}+c\left(\overline{X}\right),\hfill \end{array}$$

(142)

where $p^{\left(1\right)}$ is defined up to an arbitrary $\overline{\mathbf{Y}}$-constant function given by $c\left(\overline{X}\right)$. The second rank tensor $\mathbf{W}$ and the vector $\mathbf{P}$ are the solution to the cell problem given by:

$$\begin{array}{cc}\hfill {\nabla }_{\overline{\mathbf{Y}}}·\left(\left({\nabla }_{\overline{\mathbf{Y}}}\mathbf{W}\right){\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}{\left({\mathbf{F}}_{\mathrm{f}}^{\left(0\right)}\right)}^{-}\mathrm{T}+{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}{\left({\nabla }_{\overline{\mathbf{Y}}}\mathbf{W}\right)}^{\mathrm{T}}{\left({\mathbf{F}}_{\mathrm{f}}^{\left(0\right)}\right)}^{-}\mathrm{T}\right)& \hfill \end{array}$$

$$\begin{array}{cc}\hfill +{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}(\mathbf{I}-{\nabla }_{\overline{\mathbf{Y}}}\mathbf{P})=0& \mathrm{in}{\Omega }_{\mathrm{f}}^0,\hfill \end{array}$$

(143)

$$\begin{array}{cc}\hfill {\nabla }_{\overline{\mathbf{Y}}}·\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\mathbf{W}\right)=0& \mathrm{in}{\Omega }_{\mathrm{f}}^0,\hfill \end{array}$$

(144)

$$\begin{array}{cc}\hfill \mathbf{W}=0& \mathrm{on}{\Gamma }_{\mathrm{I}}^0\cup {\Gamma }_{\mathrm{II}}^0,\hfill \end{array}$$

(145)

This cell problem is to be supplemented by periodic conditions on the boundary $\partial {\Omega }_{\mathrm{f}}^0\backslash ({\Gamma }_{\mathrm{I}}^0\cup {\Gamma }_{\mathrm{II}}^0)$, and for the uniqueness of the solution, a further condition on the auxiliary variable $\mathbf{P}$ is required, i.e., ${\langle \mathbf{P}\rangle }_{\mathrm{f}}=0$. Since the quantity $\mathbf{w}$ retains a dependence on the pore-scale, we take the integral average of (141) over the fluid domain, which leads to:

$${\langle \mathbf{w}\rangle }_{\mathrm{f}}=-{\langle \mathbf{W}\rangle }_{\mathrm{f}}{\nabla }_{\overline{\mathbf{X}}}{p}^{\left(0\right)}$$

(146)

Therefore, the macroscale fluid flow is described by Darcy’s law.

We also consider the incompressibility constraint (109), and we integrate over the fluid domain to obtain:

$${\int }_{{\Omega }_{\mathrm{f}}^0}{\nabla }_{\overline{\mathbf{X}}}·\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}{\mathbf{V}}^{\left(0\right)}\right)\mathrm{d}\overline{\mathrm{Y}}+{\int }_{{\Omega }_{\mathrm{f}}^0}{\nabla }_{\overline{\mathbf{Y}}}·({\mathsf{G}}_{\mathrm{f}}^{\left(1\right)}{\mathbf{V}}^{\left(0\right)}+{\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}{\mathbf{V}}^{\left(1\right)})\mathrm{d}\overline{\mathrm{Y}}=0.$$

(147)

Applying the divergence theorem to the second integral and using (94), (95), (110) and (111) give:

$$\begin{array}{cc}\hfill & {\int }_{{\Omega }_{\mathrm{f}}^0}{\nabla }_{\overline{\mathbf{X}}}·\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}{\mathbf{V}}^{\left(0\right)}\right)\mathrm{d}\overline{\mathrm{Y}}+{\int }_{{\Gamma }_{\mathrm{I}}^0}({\mathsf{G}}_{\mathrm{f}}^{\left(1\right)}{\dot{\mathbf{u}}}^{\left(0\right)}+{\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}{\dot{\mathbf{u}}}_{\mathrm{I}}^{\left(1\right)}){\mathbf{n}}_{\mathrm{I}}^0\mathrm{dS}\hfill \\ \hfill & +{\int }_{{\Gamma }_{\mathrm{II}}^0}({\mathsf{G}}_{\mathrm{f}}^{\left(1\right)}{\dot{\mathbf{u}}}^{\left(0\right)}+{\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}{\dot{\mathbf{u}}}_{\mathrm{II}}^{\left(1\right)}){\mathbf{n}}_{\mathrm{II}}^0\mathrm{dS}=0.\hfill \end{array}$$

(148)

We are able to rewrite this as:

$$\begin{array}{cc}\hfill & {\int }_{{\Omega }_{\mathrm{f}}^0}{\nabla }_{\overline{\mathbf{X}}}·\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}{\mathbf{V}}^{\left(0\right)}\right)\mathrm{d}\overline{\mathrm{Y}}+{\int }_{{\Gamma }_{\mathrm{I}}^0}({\dot{\mathbf{u}}}^{\left(0\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(1\right)}\right)}^{\mathrm{T}}+{\dot{\mathbf{u}}}_{\mathrm{I}}^{\left(1\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}){\mathbf{n}}_{\mathrm{I}}^0\mathrm{dS}+{\int }_{{\Gamma }_{\mathrm{II}}^0}({\dot{\mathbf{u}}}^{\left(0\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(1\right)}\right)}^{\mathrm{T}}\hfill \\ \hfill & +{\dot{\mathbf{u}}}_{\mathrm{II}}^{\left(1\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}){\mathbf{n}}_{\mathrm{II}}^0\mathrm{dS}=0,\hfill \end{array}$$

(149)

and accounting for the continuity of the transpose of the Piola transformations applied to the normal (130) and (131) gives:

$$\begin{array}{cc}\hfill & {\int }_{{\Omega }_{\mathrm{f}}^0}{\nabla }_{\overline{\mathbf{X}}}·\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}{\mathbf{V}}^{\left(0\right)}\right)\mathrm{d}\overline{\mathrm{Y}}+{\int }_{{\Gamma }_{\mathrm{I}}^0}({\dot{\mathbf{u}}}^{\left(0\right)}{\left({\mathsf{G}}_{\mathrm{I}}^{\left(1\right)}\right)}^{\mathrm{T}}+{\dot{\mathbf{u}}}_{\mathrm{I}}^{\left(1\right)}{\left({\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}\right)}^{\mathrm{T}}){\mathbf{n}}_{\mathrm{I}}^0\mathrm{dS}\hfill \\ \hfill & +{\int }_{{\Gamma }_{\mathrm{II}}^0}({\dot{\mathbf{u}}}^{\left(0\right)}{\left({\mathsf{G}}_{\mathrm{II}}^{\left(1\right)}\right)}^{\mathrm{T}}+{\dot{\mathbf{u}}}_{\mathrm{II}}^{\left(1\right)}{\left({\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}\right)}^{\mathrm{T}}){\mathbf{n}}_{\mathrm{II}}^0\mathrm{dS}=0.\hfill \end{array}$$

(150)

This can be rewritten as:

$$\begin{array}{cc}\hfill & {\int }_{{\Omega }_{\mathrm{f}}^0}{\nabla }_{\overline{\mathbf{X}}}·\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}{\mathbf{V}}^{\left(0\right)}\right)\mathrm{d}\overline{\mathrm{Y}}+{\int }_{{\Gamma }_{\mathrm{I}}^0}({\mathsf{G}}_{\mathrm{I}}^{\left(1\right)}{\dot{\mathbf{u}}}^{\left(0\right)}+{\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}{\dot{\mathbf{u}}}_{\mathrm{I}}^{\left(1\right)}){\mathbf{n}}_{\mathrm{I}}^0\mathrm{dS}\hfill \\ \hfill & +{\int }_{{\Gamma }_{\mathrm{II}}^0}({\mathsf{G}}_{\mathrm{II}}^{\left(1\right)}{\dot{\mathbf{u}}}^{\left(0\right)}+{\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}{\dot{\mathbf{u}}}_{\mathrm{II}}^{\left(1\right)}){\mathbf{n}}_{\mathrm{II}}^0\mathrm{dS}=0,\hfill \end{array}$$

(151)

where we use the notation that $\dot{\left(\right)}=\partial \left(\right)/\partial t$. We wish to apply the divergence theorem again; however, to do this, we must include the terms on the interface ${\Gamma }_{\mathrm{III}}$, that is:

$$\begin{array}{cc}\hfill & {\int }_{{\Omega }_{\mathrm{f}}^0}{\nabla }_{\overline{\mathbf{X}}}·\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}{\mathbf{V}}^{\left(0\right)}\right)\mathrm{d}\overline{\mathrm{Y}}+{\int }_{{\Gamma }_{\mathrm{I}}^0}({\mathsf{G}}_{\mathrm{I}}^{\left(1\right)}{\dot{\mathbf{u}}}^{\left(0\right)}+{\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}{\dot{\mathbf{u}}}_{\mathrm{I}}^{\left(1\right)}){\mathbf{n}}_{\mathrm{I}}^0\mathrm{dS}\hfill \\ \hfill +& {\int }_{{\Gamma }_{\mathrm{II}}^0}({\mathsf{G}}_{\mathrm{II}}^{\left(1\right)}{\dot{\mathbf{u}}}^{\left(0\right)}+{\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}{\dot{\mathbf{u}}}_{\mathrm{II}}^{\left(1\right)}){\mathbf{n}}_{\mathrm{II}}^0\mathrm{dS}+{\int }_{{\Gamma }_{\mathrm{III}}^0}({\dot{\mathbf{u}}}^{\left(0\right)}{\left({\mathsf{G}}_{\mathrm{I}}^{\left(1\right)}\right)}^{\mathrm{T}}+{\dot{\mathbf{u}}}_{\mathrm{I}}^{\left(1\right)}{\left({\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}\right)}^{\mathrm{T}}){\mathbf{n}}_{\mathrm{III}}^0\mathrm{dS}\hfill \\ \hfill -& {\int }_{{\Gamma }_{\mathrm{III}}^0}({\dot{\mathbf{u}}}^{\left(0\right)}{\left({\mathsf{G}}_{\mathrm{II}}^{\left(1\right)}\right)}^{\mathrm{T}}+{\dot{\mathbf{u}}}_{\mathrm{II}}^{\left(1\right)}{\left({\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}\right)}^{\mathrm{T}}){\mathbf{n}}_{\mathrm{III}}^0\mathrm{dS}=0,\hfill \end{array}$$

(152)

where we can add these terms on ${\Gamma }_{\mathrm{III}}$ because they are effectively zero because of the continuity of the transpose of the Piola transform (130) and (131). Applying the divergence theorem again gives:

$$\begin{array}{cc}\hfill & {\nabla }_{\overline{\mathbf{X}}}·{\langle {\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}{\mathbf{V}}^{\left(0\right)}\rangle }_{\mathrm{f}}-{\int }_{{\Omega }_{\mathrm{I}}^0}{\nabla }_{\overline{\mathbf{Y}}}·({\mathsf{G}}_{\mathrm{I}}^{\left(1\right)}{\dot{\mathbf{u}}}^{\left(0\right)}+{\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}{\dot{\mathbf{u}}}_{\mathrm{I}}^{\left(1\right)})\mathrm{d}\overline{\mathrm{Y}}\hfill \\ \hfill & -{\int }_{{\Omega }_{\mathrm{II}}^0}{\nabla }_{\overline{\mathbf{Y}}}·({\mathsf{G}}_{\mathrm{II}}^{\left(1\right)}{\dot{\mathbf{u}}}^{\left(0\right)}+{\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}{\dot{\mathbf{u}}}_{\mathrm{II}}^{\left(1\right)})\mathrm{d}\overline{\mathrm{Y}}=0.\hfill \end{array}$$

(153)

Therefore, we can write this as:

$${\nabla }_{\overline{\mathbf{X}}}·{\langle {\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}{\mathbf{V}}^{\left(0\right)}\rangle }_{\mathrm{f}}={\langle {\nabla }_{\overline{\mathbf{Y}}}·({\mathsf{G}}_{\mathrm{I}}^{\left(1\right)}{\dot{\mathbf{u}}}^{\left(0\right)}+{\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}{\dot{\mathbf{u}}}_{\mathrm{I}}^{\left(1\right)})\rangle }_{\mathrm{I}}+{\langle {\nabla }_{\overline{\mathbf{Y}}}·({\mathsf{G}}_{\mathrm{II}}^{\left(1\right)}{\dot{\mathbf{u}}}^{\left(0\right)}+{\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}{\dot{\mathbf{u}}}_{\mathrm{II}}^{\left(1\right)})\rangle }_{\mathrm{II}}.$$

(154)

Since we have that $\mathbf{w}={\mathbf{V}}^{\left(0\right)}-{\dot{\mathbf{u}}}^{\left(0\right)}$ from (136), then this can be rearranged, and multiplying by ${\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}$ gives ${\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}{\mathbf{V}}^{\left(0\right)}={\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\mathbf{w}+{\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}{\dot{\mathbf{u}}}^{\left(0\right)}$. We then take the integral average over the fluid domain, which gives ${\langle {\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}{\mathbf{V}}^{\left(0\right)}\rangle }_{\mathrm{f}}={\langle {\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\mathbf{w}\rangle }_{\mathrm{f}}+{\langle {\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\rangle }_{\mathrm{f}}{\dot{\mathbf{u}}}^{\left(0\right)}$, and this can be used to replace the LHS of (154). Therefore, we have:

$$\begin{array}{cc}\hfill {\nabla }_{\overline{\mathbf{X}}}·\left({\langle {\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\mathbf{w}\rangle }_{\mathrm{f}}+{\langle {\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\rangle }_{\mathrm{f}}{\dot{\mathbf{u}}}^{\left(0\right)}\right)& ={\langle {\nabla }_{\overline{\mathbf{Y}}}·({\mathsf{G}}_{\mathrm{I}}^{\left(1\right)}{\dot{\mathbf{u}}}^{\left(0\right)}+{\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}{\dot{\mathbf{u}}}_{\mathrm{I}}^{\left(1\right)})\rangle }_{\mathrm{I}}\hfill \\ \hfill & +{\langle {\nabla }_{\overline{\mathbf{Y}}}·({\mathsf{G}}_{\mathrm{II}}^{\left(1\right)}{\dot{\mathbf{u}}}^{\left(0\right)}+{\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}{\dot{\mathbf{u}}}_{\mathrm{II}}^{\left(1\right)})\rangle }_{\mathrm{II}}.\hfill \end{array}$$

(155)

Then, we wish to expand the two terms on the RHS of (155), which gives:

$$\begin{array}{cc}\hfill & {\nabla }_{\overline{\mathbf{Y}}}·({\mathsf{G}}_{\mathrm{I}}^{\left(1\right)}{\dot{\mathbf{u}}}^{\left(0\right)}+{\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}{\dot{\mathbf{u}}}_{\mathrm{I}}^{\left(1\right)})+{\nabla }_{\overline{\mathbf{Y}}}·({\mathsf{G}}_{\mathrm{II}}^{\left(1\right)}{\dot{\mathbf{u}}}^{\left(0\right)}+{\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}{\dot{\mathbf{u}}}_{\mathrm{II}}^{\left(1\right)})\hfill \\ \hfill & =\mathrm{Tr}{\nabla }_{\overline{\mathbf{Y}}}({\mathsf{G}}_{\mathrm{I}}^{\left(1\right)}{\dot{\mathbf{u}}}^{\left(0\right)}+{\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}{\dot{\mathbf{u}}}_{\mathrm{I}}^{\left(1\right)})+\mathrm{Tr}{\nabla }_{\overline{\mathbf{Y}}}({\mathsf{G}}_{\mathrm{II}}^{\left(1\right)}{\dot{\mathbf{u}}}^{\left(0\right)}+{\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}{\dot{\mathbf{u}}}_{\mathrm{II}}^{\left(1\right)})\hfill \\ \hfill & =({\nabla }_{\overline{\mathbf{Y}}}·{\left({\mathsf{G}}_{\mathrm{I}}^{\left(1\right)}\right)}^{\mathrm{T}})·{\dot{\mathbf{u}}}^{\left(0\right)}+({\nabla }_{\overline{\mathbf{Y}}}·{\left({\mathsf{G}}_{\mathrm{II}}^{\left(1\right)}\right)}^{\mathrm{T}})·{\dot{\mathbf{u}}}^{\left(0\right)}+{\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}:{\nabla }_{\overline{\mathbf{Y}}}{\dot{\mathbf{u}}}_{\mathrm{I}}^{\left(1\right)}\hfill \\ \hfill & +{\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}:{\nabla }_{\overline{\mathbf{Y}}}{\dot{\mathbf{u}}}_{\mathrm{II}}^{\left(1\right)}+({\nabla }_{\overline{\mathbf{Y}}}·{\left({\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}\right)}^{\mathrm{T}})·{\dot{\mathbf{u}}}_{\mathrm{I}}^{\left(1\right)}+({\nabla }_{\overline{\mathbf{Y}}}·{\left({\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}\right)}^{\mathrm{T}})·{\dot{\mathbf{u}}}_{\mathrm{II}}^{\left(1\right)}\hfill \\ \hfill & +{\mathsf{G}}_{\mathrm{I}}^{\left(1\right)}:{\nabla }_{\overline{\mathbf{Y}}}{\dot{\mathbf{u}}}^{\left(0\right)}+{\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}:{\nabla }_{\overline{\mathbf{Y}}}{\dot{\mathbf{u}}}^{\left(0\right)}.\hfill \end{array}$$

(156)

We can cancel the terms in (156) due to ${\mathbf{u}}^{\left(0\right)}={\mathbf{u}}^{\left(0\right)}(\overline{X},t)$ and using (127) and then rewrite (155) using (156) as:

$$\begin{array}{cc}\hfill {\nabla }_{\overline{\mathbf{X}}}·\left({\langle {\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\mathbf{w}\rangle }_{\mathrm{f}}+{\langle {\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\rangle }_{\mathrm{f}}{\dot{\mathbf{u}}}^{\left(0\right)}\right)& =-\left({\langle {\nabla }_{\overline{\mathbf{X}}}·{\left({\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}\right)}^{\mathrm{T}}\rangle }_{\mathrm{I}}+{\langle {\nabla }_{\overline{\mathbf{X}}}·{\left({\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}\right)}^{\mathrm{T}}\rangle }_{\mathrm{II}}\right)·{\dot{\mathbf{u}}}^{\left(0\right)}\hfill \\ \hfill & +{\langle {\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}:{\nabla }_{\overline{\mathbf{Y}}}{\dot{\mathbf{u}}}_{\mathrm{I}}^{\left(1\right)}\rangle }_{\mathrm{I}}+{\langle {\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}:{\nabla }_{\overline{\mathbf{Y}}}{\dot{\mathbf{u}}}_{\mathrm{II}}^{\left(1\right)}\rangle }_{\mathrm{II}},\hfill \end{array}$$

(157)

where we also use (128) to replace the first two terms in (156). We return to this expression in Section 4.

3.5. The Macroscale Poroelastic Relationships

We require the macroscale constitutive relationship. To begin, we sum up the integral averages of Equations (106)–(108), that is,

$$\begin{array}{cc}\hfill & {\int }_{{\Omega }_{\mathrm{I}}^0}{\nabla }_{\overline{\mathbf{Y}}}·{\sigma }_{\mathrm{I}}^{\left(1\right)}\mathrm{d}\overline{\mathrm{Y}}+{\int }_{{\Omega }_{\mathrm{II}}^0}{\nabla }_{\overline{\mathbf{Y}}}·{\sigma }_{\mathrm{II}}^{\left(1\right)}\mathrm{d}\overline{\mathrm{Y}}+{\int }_{{\Omega }_{\mathrm{f}}^0}{\nabla }_{\overline{\mathbf{Y}}}·\left({{\sigma }_{\mathrm{f}}}^{\left(1\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}\right)\mathrm{d}\overline{\mathrm{Y}}\hfill \\ \hfill & +{\int }_{{\Omega }_{\mathrm{f}}^0}{\nabla }_{\overline{\mathbf{Y}}}·\left({{\sigma }_{\mathrm{f}}}^{\left(0\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(1\right)}\right)}^{\mathrm{T}}\right)\mathrm{d}\overline{\mathrm{Y}}+{\int }_{{\Omega }_{\mathrm{I}}^0}{\nabla }_{\overline{\mathbf{X}}}·{\sigma }_{\mathrm{I}}^{\left(0\right)}\mathrm{d}\overline{\mathrm{Y}}+{\int }_{{\Omega }_{\mathrm{II}}^0}{\nabla }_{\overline{\mathbf{X}}}·{\sigma }_{\mathrm{II}}^{\left(0\right)}\mathrm{d}\overline{\mathrm{Y}}\hfill \\ \hfill & +{\int }_{{\Omega }_{\mathrm{f}}^0}{\nabla }_{\overline{\mathbf{X}}}·\left({{\sigma }_{\mathrm{f}}}^{\left(0\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}\right)\mathrm{d}\overline{\mathrm{Y}}=0.\hfill \end{array}$$

(158)

We apply the divergence theorem to the first four integrals and then rearrange the last three integrals due to the assumption of macroscopic uniformity (72) to obtain:

$$\begin{array}{cc}\hfill & {\int }_{\partial {\Omega }_{\mathrm{I}}^0\backslash \left({\Gamma }_{\mathrm{I}}^0\cup {\Gamma }_{\mathrm{III}}^0\right)}{\sigma }_{\mathrm{I}}^{\left(1\right)}{\mathbf{n}}_{{\Omega }_{\mathrm{I}}^0\backslash \left({\Gamma }_{\mathrm{I}}^0\cup {\Gamma }_{\mathrm{III}}^0\right)}^0\mathrm{dS}+{\int }_{{\Gamma }_{\mathrm{I}}^0}{\sigma }_{\mathrm{I}}^{\left(1\right)}{\mathbf{n}}_{\mathrm{I}}^0\mathrm{dS}-{\int }_{{\Gamma }_{\mathrm{III}}^0}{\sigma }_{\mathrm{I}}^{\left(1\right)}{\mathbf{n}}_{\mathrm{III}}^0\mathrm{dS}+\hfill \\ \hfill & {\int }_{\partial {\Omega }_{\mathrm{II}}^0\backslash ({\Gamma }_{\mathrm{II}}^0\cup {\Gamma }_{\mathrm{III}}^0)}{\sigma }_{\mathrm{II}}^{\left(1\right)}{\mathbf{n}}_{{\Omega }_{\mathrm{II}}^0\backslash ({\Gamma }_{\mathrm{II}}^0\cup {\Gamma }_{\mathrm{III}}^0)}^0\mathrm{dS}+{\int }_{{\Gamma }_{\mathrm{II}}^0}{\sigma }_{\mathrm{II}}^{\left(1\right)}{\mathbf{n}}_{\mathrm{II}}^0\mathrm{dS}+{\int }_{{\Gamma }_{\mathrm{III}}^0}{\sigma }_{\mathrm{II}}^{\left(1\right)}{\mathbf{n}}_{\mathrm{III}}^0\mathrm{dS}+\hfill \\ \hfill & {\int }_{\partial {\Omega }_{\mathrm{f}}^0\backslash ({\Gamma }_{\mathrm{I}}^0\cup {\Gamma }_{\mathrm{II}}^0)}{\sigma }_{\mathrm{f}}^{\left(1\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}{\mathbf{n}}_{{\Omega }_{\mathrm{f}}^0\backslash ({\Gamma }_{\mathrm{I}}^0\cup {\Gamma }_{\mathrm{II}}^0)}^0\mathrm{dS}-{\int }_{{\Gamma }_{\mathrm{II}}^0}{\sigma }_{\mathrm{f}}^{\left(1\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}{\mathbf{n}}_{\mathrm{II}}^0\mathrm{dS}-{\int }_{{\Gamma }_{\mathrm{I}}^0}{\sigma }_{\mathrm{f}}^{\left(1\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}{\mathbf{n}}_{\mathrm{I}}^0\mathrm{dS}+\hfill \\ \hfill & {\int }_{\partial {\Omega }_{\mathrm{f}}^0\backslash ({\Gamma }_{\mathrm{I}}^0\cup {\Gamma }_{\mathrm{II}}^0)}{\sigma }_{\mathrm{f}}^{\left(0\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(1\right)}\right)}^{\mathrm{T}}{\mathbf{n}}_{{\Omega }_{\mathrm{f}}^0\backslash ({\Gamma }_{\mathrm{I}}^0\cup {\Gamma }_{\mathrm{II}}^0)}^0\mathrm{dS}-{\int }_{{\Gamma }_{\mathrm{II}}^0}{\sigma }_{\mathrm{f}}^{\left(0\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(1\right)}\right)}^{\mathrm{T}}{\mathbf{n}}_{\mathrm{II}}^0\mathrm{dS}-{\int }_{{\Gamma }_{\mathrm{I}}^0}{\sigma }_{\mathrm{f}}^{\left(0\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(1\right)}\right)}^{\mathrm{T}}{\mathbf{n}}_{\mathrm{I}}^0\mathrm{dS}+\hfill \\ \hfill & {\nabla }_{\overline{\mathbf{X}}}·{\int }_{{\Omega }_{\mathrm{I}}^0}{\sigma }_{\mathrm{I}}^{\left(0\right)}\mathrm{d}\overline{\mathrm{Y}}+{\nabla }_{\overline{\mathbf{X}}}·{\int }_{{\Omega }_{\mathrm{II}}^0}{\sigma }_{\mathrm{II}}^{\left(0\right)}\mathrm{d}\overline{\mathrm{Y}}+{\nabla }_{\overline{\mathbf{X}}}·{\int }_{{\Omega }_{\mathrm{f}}^0}{{\sigma }_{\mathrm{f}}}^{\left(0\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}\mathrm{d}\overline{\mathrm{Y}}=0,\hfill \end{array}$$

(159)

where ${\mathbf{n}}_{\mathrm{I}}^0$, ${\mathbf{n}}_{\mathrm{II}}^0$, ${\mathbf{n}}_{\mathrm{III}}^0$, ${\mathbf{n}}_{{\Omega }_{\mathrm{I}}^0\backslash \left({\Gamma }_{\mathrm{I}}^0\cup {\Gamma }_{\mathrm{III}}^0\right)}^0$, ${\mathbf{n}}_{{\Omega }_{\mathrm{II}}^0\backslash \left({\Gamma }_{\mathrm{II}}^0\cup {\Gamma }_{\mathrm{III}}^0\right)}^0$ and ${\mathbf{n}}_{{\Omega }_{\mathrm{f}}^0\backslash \left({\Gamma }_{\mathrm{I}}^0\cup {\Gamma }_{\mathrm{II}}^0\right)}^0$ are the unit normals corresponding to the interfaces ${\Gamma }_{\mathrm{I}}^0$, ${\Gamma }_{\mathrm{II}}^0$, ${\Gamma }_{\mathrm{III}}^0$, $\partial {\Omega }_{\mathrm{I}}^0\backslash ({\Gamma }_{\mathrm{I}}^0\cup {\Gamma }_{\mathrm{III}}^0)$, $\partial {\Omega }_{\mathrm{II}}^0\backslash ({\Gamma }_{\mathrm{II}}^0\cup {\Gamma }_{\mathrm{III}}^0)$ and $\partial {\Omega }_{\mathrm{f}}^0\backslash ({\Gamma }_{\mathrm{I}}^0\cup {\Gamma }_{\mathrm{II}}^0)$. The contributions over the external boundaries of ${\Omega }_{\mathrm{I}}^0$, ${\Omega }_{\mathrm{II}}^0$ and ${\Omega }_{\mathrm{f}}^0$ cancel due to $\overline{\mathbf{Y}}$-periodicity, so (159) becomes:

$$\begin{array}{cc}\hfill & {\int }_{{\Gamma }_{\mathrm{I}}^0}{\sigma }_{\mathrm{I}}^{\left(1\right)}{\mathbf{n}}_{\mathrm{I}}^{0}dS-{\int }_{{\Gamma }_{\mathrm{I}}^0}{\sigma }_{\mathrm{f}}^{\left(1\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}{\mathbf{n}}_{\mathrm{I}}^{0}dS-{\int }_{{\Gamma }_{\mathrm{I}}^0}{\sigma }_{\mathrm{f}}^{\left(0\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(1\right)}\right)}^{\mathrm{T}}{\mathbf{n}}_{\mathrm{I}}^{0}dS+{\int }_{{\Gamma }_{\mathrm{II}}^0}{\sigma }_{\mathrm{II}}^{\left(1\right)}{\mathbf{n}}_{\mathrm{II}}^{0}dS-\hfill \\ \hfill & {\int }_{{\Gamma }_{\mathrm{II}}^0}{\sigma }_{\mathrm{f}}^{\left(1\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}{\mathbf{n}}_{\mathrm{II}}^{0}dS-{\int }_{{\Gamma }_{\mathrm{II}}^0}{\sigma }_{\mathrm{f}}^{\left(0\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(1\right)}\right)}^{\mathrm{T}}{\mathbf{n}}_{\mathrm{II}}^{0}dS-{\int }_{{\Gamma }_{\mathrm{III}}^0}{\sigma }_{\mathrm{I}}^{\left(1\right)}{\mathbf{n}}_{\mathrm{III}}^0\mathrm{dS}+{\int }_{{\Gamma }_{\mathrm{III}}^0}{\sigma }_{\mathrm{II}}^{\left(1\right)}{\mathbf{n}}_{\mathrm{III}}^0\mathrm{dS}+\hfill \\ \hfill & {\nabla }_{\overline{\mathbf{X}}}·{\int }_{{\Omega }_{\mathrm{I}}^0}{\sigma }_{\mathrm{I}}^{\left(0\right)}\mathrm{d}\overline{\mathrm{Y}}+{\nabla }_{\overline{\mathbf{X}}}·{\int }_{{\Omega }_{\mathrm{II}}^0}{\sigma }_{\mathrm{II}}^{\left(0\right)}\mathrm{d}\overline{\mathrm{Y}}+{\nabla }_{\overline{\mathbf{X}}}·{\int }_{{\Omega }_{\mathrm{f}}^0}{{\sigma }_{\mathrm{f}}}^{\left(0\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}\mathrm{d}\overline{\mathrm{Y}}=0.\hfill \end{array}$$

(160)

The first eight integrals cancel using the continuity of stresses (112)–(114), so we obtain:

$${\nabla }_{\overline{\mathbf{X}}}·{\langle {\sigma }_{\mathrm{I}}^{\left(0\right)}\rangle }_{\mathrm{I}}+{\nabla }_{\overline{\mathbf{X}}}·{\langle {\sigma }_{\mathrm{II}}^{\left(0\right)}\rangle }_{\mathrm{II}}-{\nabla }_{\overline{\mathbf{X}}}·\left(p^{\left(0\right)}{\langle {\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}\rangle }_{\mathrm{f}}\right)=0,$$

(161)

where we use (100). Then, we have that:

$${\nabla }_{\overline{\mathbf{X}}}·{\sigma }_{\mathrm{Eff}}=0,$$

(162)

where:

$${\sigma }_{\mathrm{Eff}}={\langle {\sigma }_{\mathrm{I}}^{\left(0\right)}\rangle }_{\mathrm{I}}+{\langle {\sigma }_{\mathrm{II}}^{\left(0\right)}\rangle }_{\mathrm{II}}-p^{\left(0\right)}{\langle {\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}\rangle }_{\mathrm{f}}.$$

(163)

We can describe (162) and (163) as the average stress balance and constitutive law for our nonlinear poroelastic composite material.

We can write the following problem for ${\mathbf{u}}_{\mathrm{I}}^{\left(1\right)}$ and ${\mathbf{u}}_{\mathrm{II}}^{\left(1\right)}$ using (90), (91), (96)–(98), (100) and (115):

$$\begin{array}{cc}\hfill {\nabla }_{\overline{\mathbf{Y}}}·{\sigma }_{\mathrm{I}}^{\left(0\right)}=0& \mathrm{in}{\Omega }_{\mathrm{I}}^0,\hfill \end{array}$$

(164)

$$\begin{array}{cc}\hfill {\nabla }_{\overline{\mathbf{Y}}}·{\sigma }_{\mathrm{II}}^{\left(0\right)}=0& \mathrm{in}{\Omega }_{\mathrm{II}}^0,\hfill \end{array}$$

(165)

$$\begin{array}{cc}\hfill {\mathbf{u}}_{\mathrm{I}}^{\left(1\right)}={\mathbf{u}}_{\mathrm{II}}^{\left(1\right)}& \mathrm{on}{\Gamma }_{\mathrm{III}}^0,\hfill \end{array}$$

(166)

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{I}}^{\left(0\right)}{\mathbf{n}}_{\mathrm{III}}^0={\sigma }_{\mathrm{II}}^{\left(0\right)}{\mathbf{n}}_{\mathrm{III}}^0& \mathrm{on}{\Gamma }_{\mathrm{III}}^0,\hfill \end{array}$$

(167)

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{I}}^{\left(0\right)}{\mathbf{n}}_{\mathrm{I}}^0=-\left(p^{\left(0\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}\right){\mathbf{n}}_{\mathrm{I}}^0& \mathrm{on}{\Gamma }_{\mathrm{I}}^0,\hfill \end{array}$$

(168)

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{II}}^{\left(0\right)}{\mathbf{n}}_{\mathrm{II}}^0=-\left(p^{\left(0\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}\right){\mathbf{n}}_{\mathrm{II}}^0& \mathrm{on}{\Gamma }_{\mathrm{II}}^0.\hfill \end{array}$$

(169)

Since our solid stress tensors ${\sigma }_{\mathrm{I}}$ and ${\sigma }_{\mathrm{II}}$ are described by a constitutive law, our problem here is the general case. It is possible to state a generalised ansatz to this problem and state the cell problems. However, we wait to do this until the constitutive law has been specified. Within the next section, we state the macroscale model for this general case and then specify this problem for a specific choice of the constitutive law.

4. The Macroscale Model and Particular Cases

We begin by stating the general nonlinear macroscale model for poroelastic composite materials undergoing large deformations where the specific constitutive law for the material has not yet been specified, that is,

$$\left\{\begin{array}{c}\begin{array}{c}{\langle \mathbf{w}\rangle }_{\mathrm{f}}=-{\langle \mathbf{W}\rangle }_{\mathrm{f}}{\nabla }_{\overline{\mathbf{X}}}{p}^{\left(0\right)},\hfill \\ {\nabla }_{\overline{\mathbf{X}}}·{\sigma }_{\mathrm{Eff}}=0,\hfill \\ {\sigma }_{\mathrm{Eff}}={〈\frac{\partial {\psi }_{\mathrm{I}}^{\left(0\right)}}{\partial {\mathbf{F}}_{\mathrm{I}}^{\left(0\right)}}〉}_{\mathrm{I}}+{〈\frac{\partial {\psi }_{\mathrm{II}}^{\left(0\right)}}{\partial {\mathbf{F}}_{\mathrm{II}}^{\left(0\right)}}〉}_{\mathrm{II}}-p^{\left(0\right)}{\langle {\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}\rangle }_{\mathrm{f}},\hfill \\ {\nabla }_{\overline{\mathbf{X}}}·{\langle {\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\mathbf{w}\rangle }_{\mathrm{f}}=-\left({\langle {\nabla }_{\overline{\mathbf{X}}}·{\left({\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}\right)}^{\mathrm{T}}\rangle }_{\mathrm{I}}+{\langle {\nabla }_{\overline{\mathbf{X}}}·{\left({\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}\right)}^{\mathrm{T}}\rangle }_{\mathrm{II}}+{\langle {\nabla }_{\overline{\mathbf{X}}}·{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}\rangle }_{\mathrm{f}}\right)·{\dot{\mathbf{u}}}^{\left(0\right)}\hfill \\ +{\langle {\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}:{\nabla }_{\overline{\mathbf{Y}}}{\dot{\mathbf{u}}}_{\mathrm{I}}^{\left(1\right)}\rangle }_{\mathrm{I}}+{\langle {\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}:{\nabla }_{\overline{\mathbf{Y}}}{\dot{\mathbf{u}}}_{\mathrm{II}}^{\left(1\right)}\rangle }_{\mathrm{II}}-{\langle {\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\rangle }_{\mathrm{f}}:{\nabla }_{\overline{\mathbf{X}}}{\dot{\mathbf{u}}}^{\left(0\right)},\hfill \end{array}\hfill \end{array}\right.$$

where ${\psi }_{\mathrm{I}}$ and ${\psi }_{\mathrm{II}}$ are the strain energy functions for the material in the subphase and matrix, respectively. This model comprises Darcy’s law for the relative fluid–solid velocity. The second and third equations are the stress balance and the constitutive law for the nonlinear poroelastic composite material. The constitutive law is to be specified further by choosing a specific strain energy function for the matrix and the subphase, relevant to the intended application. The final equation is the conservation of mass equation. This equation is also be influenced by the choice of strain energy function for the matrix and the subphase. Within the next subsection, we prescribe a particular constitutive law for the elastic materials, and this leads to a specified model that reduces to three previously known models under appropriate simplifying assumptions.

4.1. Constitutive Law

We choose a very simple constitutive law for our solid material, which is the de Saint-Venant strain energy function given by:

$${\psi }_v=\frac{1}{2}{\mathbf{E}}_v:{\mathbb{C}}_v:{\mathbf{E}}_v\mathrm{in}{\Omega }_{\mathrm{v}}^0\mathrm{where}v=\mathrm{I},\mathrm{II},$$

(170)

where ${\psi }_v$ is described by different parameters depending on the solid domain that it is describing and ${\mathbf{E}}_v$ is the Green-Lagrangian strain tensor for each solid domain, that is,

$${\psi }_v=\left\{\begin{array}{c}\begin{array}{c}{\psi }_{\mathrm{I}}=\frac{1}{2}{\mathbf{E}}_{\mathrm{I}}:{\mathbb{C}}_{\mathrm{I}}:{\mathbf{E}}_{\mathrm{I}}\mathrm{in}{\Omega }_{\mathrm{I}}^0,\hfill \\ {\psi }_{\mathrm{II}}=\frac{1}{2}{\mathbf{E}}_{\mathrm{II}}:{\mathbb{C}}_{\mathrm{II}}:{\mathbf{E}}_{\mathrm{II}}\mathrm{in}{\Omega }_{\mathrm{II}}^0.\hfill \end{array}\hfill \end{array}\right.$$

(171)

We adopt the notation that:

$$\mathbb{C}=\left\{\begin{array}{c}\begin{array}{c}{\mathbb{C}}_{\mathrm{I}}\mathrm{in}{\Omega }_{\mathrm{I}}^0\hfill \\ {\mathbb{C}}_{\mathrm{II}}\mathrm{in}{\Omega }_{\mathrm{II}}^0\hfill \end{array}\hfill \end{array}\right.$$

(172)

and the subscript $v=\mathrm{I},\mathrm{II}$, where $\mathrm{I}$ is the subphase and $\mathrm{II}$ is the matrix throughout this section. We have that the expansion of ${\psi }_v$ is given by:

$${\psi }_v^{ϵ}=\frac{1}{2}{\mathbf{E}}_v^{ϵ}:\mathbb{C}:{\mathbf{E}}_v^{ϵ},$$

(173)

and we note that $\mathbb{C}$ is the fourth rank elasticity tensor with major and minor symmetries, namely:

$$C_{ijkl}=C_{jikl}=C_{ijlk}=C_{klij}.$$

(174)

We can now determine ${\mathbf{E}}_v$. We have:

$$\begin{array}{cc}\hfill {\mathbf{E}}_v& =\frac{1}{2}({\left({\mathbf{F}}_v\right)}^{\mathrm{T}}{\mathbf{F}}_v-\mathbf{I})=\frac{1}{2}({({\nabla }_{\mathbf{X}}{\mathbf{u}}_v+\mathbf{I})}^{\mathrm{T}}({\nabla }_{\mathbf{X}}{\mathbf{u}}_v+\mathbf{I})-\mathbf{I})\hfill \\ \hfill & =\frac{1}{2}({\left({\nabla }_{\mathbf{X}}{\mathbf{u}}_v\right)}^{\mathrm{T}}{\nabla }_{\mathbf{X}}{\mathbf{u}}_v+{\left({\nabla }_{\mathbf{X}}{\mathbf{u}}_v\right)}^{\mathrm{T}}+{\nabla }_{\mathbf{X}}{\mathbf{u}}_v),\hfill \end{array}$$

(175)

where the subscript $v=\mathrm{I},\mathrm{II}$ denotes the subphase and the matrix. We now apply the asymptotic homogenization technique to (175). Using the transformation of the gradient operator (70):

$${\nabla }_{\mathbf{X}}\to {\nabla }_{\overline{\mathbf{X}}}+\frac{1}{ϵ}{\nabla }_{\overline{\mathbf{Y}}},$$

(176)

we can write ${\mathbf{E}}_v^{ϵ}$ as:

$$\begin{array}{cc}\hfill {\mathbf{E}}_v^{ϵ}=& \frac{1}{2}\left[{({\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_v^{ϵ}+\frac{1}{ϵ}{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{ϵ})}^{\mathrm{T}}({\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_v^{ϵ}+\frac{1}{ϵ}{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{ϵ})+{({\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_v^{ϵ}+\frac{1}{ϵ}{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{ϵ})}^{\mathrm{T}}+({\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_v^{ϵ}+\frac{1}{ϵ}{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{ϵ})\right]\hfill \\ \hfill =& \frac{1}{2}[{\left({\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_v^{ϵ}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_v^{ϵ}+\frac{1}{ϵ}{\left({\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_v^{ϵ}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{ϵ}+\frac{1}{ϵ}{\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{ϵ}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_v^{ϵ}+\frac{1}{{ϵ}^2}{\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{ϵ}\right)}^{\mathrm{T}}\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{ϵ}\right)\hfill \\ \hfill & +{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_v^{ϵ}+\frac{1}{ϵ}{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{ϵ}+{\left({\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_v^{ϵ}\right)}^{\mathrm{T}}+\frac{1}{ϵ}{\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{ϵ}\right)}^{\mathrm{T}}].\hfill \end{array}$$

(177)

Then, using (177), we are able to write the expansion of the strain energy function ${\psi }_v^{ϵ}$ as:

$$\begin{array}{cc}\hfill {\psi }_v^{ϵ}=& \frac{1}{8}\left[\right({\left({\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_v^{ϵ}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_v^{ϵ}+\frac{1}{ϵ}{\left({\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_v^{ϵ}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{ϵ}+\frac{1}{ϵ}{\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{ϵ}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_v^{ϵ}+\frac{1}{{ϵ}^2}{\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{ϵ}\right)}^{\mathrm{T}}\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{ϵ}\right)\hfill \\ \hfill & +{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_v^{ϵ}+\frac{1}{ϵ}{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{ϵ}+{\left({\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_v^{ϵ}\right)}^{\mathrm{T}}+\frac{1}{ϵ}{\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{ϵ}\right)}^{\mathrm{T}}):\mathbb{C}:({\left({\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_v^{ϵ}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_v^{ϵ}\hfill \\ \hfill & +\frac{1}{ϵ}{\left({\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_v^{ϵ}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{ϵ}+\frac{1}{ϵ}{\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{ϵ}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_v^{ϵ}+\frac{1}{{ϵ}^2}{\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{ϵ}\right)}^{\mathrm{T}}\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{ϵ}\right)\hfill \\ \hfill & +{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_v^{ϵ}+\frac{1}{ϵ}{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{ϵ}+{\left({\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_v^{ϵ}\right)}^{\mathrm{T}}+\frac{1}{ϵ}{\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{ϵ}\right)}^{\mathrm{T}}\left)\right].\hfill \end{array}$$

(178)

Multiplying (178) by ${ϵ}^2$, we obtain:

$$\begin{array}{cc}\hfill {ϵ}^2{\psi }_v^{ϵ}=& \frac{1}{8}\left[\right({ϵ}^2{\left({\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_v^{ϵ}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_v^{ϵ}+ϵ{\left({\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_v^{ϵ}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{ϵ}+ϵ{\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{ϵ}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_v^{ϵ}+{\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{ϵ}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{ϵ}\hfill \\ \hfill & +{ϵ}^2{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_v^{ϵ}+ϵ{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{ϵ}+{ϵ}^2{\left({\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_v^{ϵ}\right)}^{\mathrm{T}}+ϵ{\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{ϵ}\right)}^{\mathrm{T}}):{ϵ}^2\mathbb{C}:({ϵ}^2{\left({\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_v^{ϵ}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_v^{ϵ}\hfill \\ \hfill & +ϵ{\left({\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_v^{ϵ}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{ϵ}+ϵ{\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{ϵ}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_v^{ϵ}+{\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{ϵ}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{ϵ}+{ϵ}^2{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_v^{ϵ}\hfill \\ \hfill & +ϵ{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{ϵ}+{ϵ}^2{\left({\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}_v^{ϵ}\right)}^{\mathrm{T}}+ϵ{\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{ϵ}\right)}^{\mathrm{T}}\left)\right].\hfill \end{array}$$

(179)

We can then equate the coefficients of ${ϵ}^0$, ${ϵ}^1$ and ${ϵ}^2$ in (179). For ${ϵ}^0$, we have:

$$0=\frac{1}{8}\left[\left({\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}^{\left(0\right)}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}^{\left(0\right)}\right):\left({\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}^{\left(0\right)}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}^{\left(0\right)}\right)\right]=0.$$

(180)

Equating the coefficients of ${ϵ}^1$ gives:

$$\begin{array}{cc}\hfill 0& =\frac{1}{8}\left[\right({\left({\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}^{\left(0\right)}+{\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}^{\left(0\right)}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}+{\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{\left(1\right)}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}^{\left(0\right)}+{\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}^{\left(0\right)}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{\left(1\right)}\hfill \\ \hfill & +{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}^{\left(0\right)}+{\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}^{\left(0\right)}\right)}^{\mathrm{T}}):({\left({\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}^{\left(0\right)}+{\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}^{\left(0\right)}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}+{\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{\left(1\right)}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}^{\left(0\right)}\hfill \\ \hfill & +{\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}^{\left(0\right)}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{\left(1\right)}+{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}^{\left(0\right)}+{\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}^{\left(0\right)}\right)}^{\mathrm{T}}\left)\right]=0,\hfill \end{array}$$

(181)

and finally, equating the coefficients of ${ϵ}^2$ gives:

$$\begin{array}{cc}\hfill {\psi }_v^{\left(0\right)}=& \frac{1}{8}\left[\right({\left({\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}+{\left({\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{\left(1\right)}+{\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{\left(1\right)}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}+{\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{\left(1\right)}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{\left(1\right)}\hfill \\ \hfill & +{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}+{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{\left(1\right)}+{\left({\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}\right)}^{\mathrm{T}}+{\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{\left(1\right)}\right)}^{\mathrm{T}}):\mathbb{C}:({\left({\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}\hfill \\ \hfill & +{\left({\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{\left(1\right)}+{\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{\left(1\right)}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}+{\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{\left(1\right)}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{\left(1\right)}+{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}+{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{\left(1\right)}\hfill \\ \hfill & +{\left({\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}\right)}^{\mathrm{T}}+{\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{\left(1\right)}\right)}^{\mathrm{T}}\left)\right]\hfill \\ \hfill =& \frac{1}{2}{\mathbf{E}}_v^{\left(0\right)}:\mathbb{C}:{\mathbf{E}}_v^{\left(0\right)},\hfill \end{array}$$

(182)

where we have that:

$$\begin{array}{cc}\hfill {\mathbf{E}}_v^{\left(0\right)}& ={\left({\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}+{\left({\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{\left(1\right)}+{\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{\left(1\right)}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}\hfill \\ \hfill & +{\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{\left(1\right)}\right)}^{\mathrm{T}}{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{\left(1\right)}+{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}+{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{\left(1\right)}+{\left({\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}\right)}^{\mathrm{T}}+{\left({\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_v^{\left(1\right)}\right)}^{\mathrm{T}}.\hfill \end{array}$$

(183)

Therefore, we have that the leading order term of the constitutive law in the subphase and the matrix respectively are:

$${\psi }_{\mathrm{I}}^{\left(0\right)}=\frac{1}{2}{\mathbf{E}}_{\mathrm{I}}^{\left(0\right)}:{\mathbb{C}}_{\mathrm{I}}:{\mathbf{E}}_{\mathrm{I}}^{\left(0\right)}\mathrm{and}{\psi }_{\mathrm{II}}^{\left(0\right)}=\frac{1}{2}{\mathbf{E}}_{\mathrm{II}}^{\left(0\right)}:{\mathbb{C}}_{\mathrm{II}}:{\mathbf{E}}_{\mathrm{II}}^{\left(0\right)}.$$

(184)

Now that we have an expression for ${\psi }_v^{\left(0\right)}$, we can use this to find the leading order term of the second Piola–Kirchoff stress, that is we take the derivatives of (184) with respect to ${\mathbf{E}}_{\mathrm{I}}^{\left(0\right)}$ and ${\mathbf{E}}_{\mathrm{II}}^{\left(0\right)}$, respectively. Therefore, we have:

$${\mathbf{S}}_v^{\left(0\right)}=\frac{\partial {\psi }_v^{\left(0\right)}}{\partial {\mathbf{E}}_v^{\left(0\right)}}=\mathbb{C}:{\mathbf{E}}_v^{\left(0\right)},$$

(185)

where we have a different second Piola–Kirchoff stress for both the subphase and the matrix, respectively, that is,

$${\mathbf{S}}_v^{\left(0\right)}=\left\{\begin{array}{c}\begin{array}{c}{\mathbf{S}}_{\mathrm{I}}^{\left(0\right)}={\mathbb{C}}_{\mathrm{I}}:{\mathbf{E}}_{\mathrm{I}}^{\left(0\right)}\mathrm{in}{\Omega }_{\mathrm{I}}^0,\hfill \\ {\mathbf{S}}_{\mathrm{II}}^{\left(0\right)}={\mathbb{C}}_{\mathrm{II}}:{\mathbf{E}}_{\mathrm{II}}^{\left(0\right)}\mathrm{in}{\Omega }_{\mathrm{II}}^0.\hfill \end{array}\hfill \end{array}\right.$$

(186)

We now wish to linearise the second Piola–Kirchoff stress ${\mathbf{S}}_v^{\left(0\right)}$. To do this, we can use the expression ${\mathbf{F}}_v^{\left(0\right)}=\mathbf{I}+\mathbf{H}$ where:

$$\mathbf{H}=\left\{\begin{array}{c}\begin{array}{c}{\mathbf{H}}_{\mathrm{I}}={\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}+{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_{\mathrm{I}}^{\left(1\right)}\mathrm{in}{\Omega }_{\mathrm{I}}^0,\hfill \\ {\mathbf{H}}_{\mathrm{II}}={\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}+{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_{\mathrm{II}}^{\left(1\right)}\mathrm{in}{\Omega }_{\mathrm{II}}^0.\hfill \end{array}\hfill \end{array}\right.$$

(187)

Therefore, carrying out the linearisation, we have that:

$$\begin{array}{cc}\hfill {\mathbf{S}}_v^{\left(0\right)}=& \mathbb{C}:{\mathbf{E}}_v^{\left(0\right)}\hfill \\ \hfill =& \mathbb{C}:\left(\frac{1}{2}({\left({\mathbf{F}}_v^{\left(0\right)}\right)}^{\mathrm{T}}{\mathbf{F}}_v^{\left(0\right)}-\mathbf{I})\right)\hfill \end{array}$$

$$\begin{array}{cc}\hfill =& \mathbb{C}:\left(\frac{1}{2}({(\mathbf{I}+\mathbf{H})}^{\mathrm{T}}(\mathbf{I}+\mathbf{H})-\mathbf{I})\right)\hfill \end{array}$$

(188)

$$\begin{array}{cc}\hfill =& \mathbb{C}:\left(\frac{1}{2}(\mathbf{H}+{\mathbf{H}}^{\mathrm{T}}+{\mathbf{H}}^{\mathrm{T}}\mathbf{H})\right).\hfill \end{array}$$

(189)

Ignoring the nonlinear terms, we have:

$${\left({\mathbf{S}}_v^{\left(0\right)}\right)}_{\mathrm{lin}}=\mathbb{C}:\left(\frac{1}{2}(\mathbf{H}+{\mathbf{H}}^{\mathrm{T}})\right)=\mathbb{C}:\mathrm{Sym}\mathbf{H}.$$

(190)

We are now able to use the second Piola–Kirchoff stress ${\mathbf{S}}_v^{\left(0\right)}$ to find the first Piola–Kirchoff stress:

$$\begin{array}{cc}\hfill {\sigma }_v^{\left(0\right)}& ={\mathbf{F}}_v^{\left(0\right)}{\mathbf{S}}_v^{\left(0\right)}\hfill \\ \hfill & =(\mathbf{I}+\mathbf{H})\mathbb{C}:{\mathbf{E}}_v^{\left(0\right)},\hfill \end{array}$$

(191)

where we can define the first Piola stress ${\sigma }^{\left(0\right)}$ in both of the solid constituents, the subphase and matrix, as:

$${\sigma }_v^{\left(0\right)}=\left\{\begin{array}{c}\begin{array}{c}{\sigma }_{\mathrm{I}}^{\left(0\right)}=(\mathbf{I}+{\mathbf{H}}_{\mathrm{I}}){\mathbb{C}}_{\mathrm{I}}:{\mathbf{E}}_{\mathrm{I}}^{\left(0\right)}\mathrm{in}{\Omega }_{\mathrm{I}}^0,\hfill \\ {\sigma }_{\mathrm{II}}^{\left(0\right)}=(\mathbf{I}+{\mathbf{H}}_{\mathrm{II}}){\mathbb{C}}_{\mathrm{II}}:{\mathbf{E}}_{\mathrm{II}}^{\left(0\right)}\mathrm{in}{\Omega }_{\mathrm{II}}^0.\hfill \end{array}\hfill \end{array}\right.$$

(192)

We can also linearise the first Piola–Kirchoff stress ${\sigma }^{\left(0\right)}$. Therefore, we have:

$${\left({\sigma }_v\right)}_{\mathrm{lin}}^{\left(0\right)}=\mathbb{C}:\mathrm{Sym}\mathbf{H},$$

(193)

and since we have that $\mathbb{C}$ is major and minor symmetric, then we can write this as:

$${\left({\sigma }_v\right)}_{\mathrm{lin}}^{\left(0\right)}=\mathbb{C}:\mathbf{H}.$$

(194)

We therefore can define the solid stresses in the subphase and matrix as:

$$\begin{array}{cc}\hfill & {\left({\sigma }_{\mathrm{I}}^{\left(0\right)}\right)}_{\mathrm{lin}}={\mathbb{C}}_{\mathrm{I}}:({\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}+{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_{\mathrm{I}}^{\left(1\right)})\mathrm{in}{\Omega }_{\mathrm{I}}^0,\hfill \end{array}$$

(195)

$$\begin{array}{cc}\hfill & {\left({\sigma }_{\mathrm{II}}^{\left(0\right)}\right)}_{\mathrm{lin}}={\mathbb{C}}_{\mathrm{II}}:({\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}+{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{u}}_{\mathrm{II}}^{\left(1\right)})\mathrm{in}{\Omega }_{\mathrm{II}}^0,\hfill \end{array}$$

(196)

respectively, for our chosen constitutive law. We can then use these expressions in the problem for ${\mathbf{u}}_{\mathrm{I}}^{\left(1\right)}$ and ${\mathbf{u}}_{\mathrm{II}}^{\left(1\right)}$ given by Equations (164)–(169), that is,

$$\begin{array}{cc}\hfill {\nabla }_{\overline{\mathbf{Y}}}·\left({\mathbb{C}}_{\mathrm{I}}{\nabla }_{\overline{\mathbf{Y}}}\left({\mathbf{u}}_{\mathrm{I}}^{\left(1\right)}\right)\right)+{\nabla }_{\overline{\mathbf{Y}}}·\left({\mathbb{C}}_{\mathrm{I}}{\nabla }_{\overline{\mathbf{X}}}\left({\mathbf{u}}^{\left(0\right)}\right)\right)=0& \mathrm{in}{\Omega }_{\mathrm{I}}^0,\hfill \end{array}$$

(197)

$$\begin{array}{cc}\hfill {\nabla }_{\overline{\mathbf{Y}}}·\left({\mathbb{C}}_{\mathrm{II}}{\nabla }_{\overline{\mathbf{Y}}}\left({\mathbf{u}}_{\mathrm{II}}^{\left(1\right)}\right)\right)+{\nabla }_{\overline{\mathbf{Y}}}·\left({\mathbb{C}}_{\mathrm{II}}{\nabla }_{\overline{\mathbf{X}}}\left({\mathbf{u}}^{\left(0\right)}\right)\right)=0& \mathrm{in}{\Omega }_{\mathrm{II}}^0,\hfill \end{array}$$

(198)

$$\begin{array}{cc}\hfill {\mathbf{u}}_{\mathrm{I}}^{\left(1\right)}={\mathbf{u}}_{\mathrm{II}}^{\left(1\right)}& \mathrm{on}{\Gamma }_{\mathrm{III}}^0,\hfill \end{array}$$

(199)

$$\begin{array}{cc}\hfill {\mathbb{C}}_{\mathrm{I}}{\nabla }_{\overline{\mathbf{Y}}}\left({\mathbf{u}}_{\mathrm{I}}^{\left(1\right)}\right){\mathbf{n}}_{\mathrm{III}}^0-{\mathbb{C}}_{\mathrm{II}}{\nabla }_{\overline{\mathbf{Y}}}\left({\mathbf{u}}_{\mathrm{II}}^{\left(1\right)}\right){\mathbf{n}}_{\mathrm{III}}^0=({\mathbb{C}}_{\mathrm{II}}-{\mathbb{C}}_{\mathrm{I}}){\nabla }_{\overline{\mathbf{X}}}\left({\mathbf{u}}^{\left(0\right)}\right){\mathbf{n}}_{\mathrm{III}}^0& \mathrm{on}{\Gamma }_{\mathrm{III}}^0,\hfill \end{array}$$

(200)

$$\begin{array}{cc}\hfill ({\mathbb{C}}_{\mathrm{I}}{\nabla }_{\overline{\mathbf{Y}}}\left({\mathbf{u}}_{\mathrm{I}}^{\left(1\right)}\right)+{\mathbb{C}}_{\mathrm{I}}{\nabla }_{\overline{\mathbf{X}}}\left({\mathbf{u}}^{\left(0\right)}\right)){\mathbf{n}}_{\mathrm{I}}^0=-\left(p^{\left(0\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}\right){\mathbf{n}}_{\mathrm{I}}^0& \mathrm{on}{\Gamma }_{\mathrm{I}}^0,\hfill \end{array}$$

(201)

$$\begin{array}{cc}\hfill ({\mathbb{C}}_{\mathrm{II}}{\nabla }_{\overline{\mathbf{Y}}}\left({\mathbf{u}}_{\mathrm{II}}^{\left(1\right)}\right)+{\mathbb{C}}_{\mathrm{II}}{\nabla }_{\overline{\mathbf{X}}}\left({\mathbf{u}}^{\left(0\right)}\right)){\mathbf{n}}_{\mathrm{II}}^0=-\left(p^{\left(0\right)}{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}\right){\mathbf{n}}_{\mathrm{II}}^0& \mathrm{on}{\Gamma }_{\mathrm{II}}^0.\hfill \end{array}$$

(202)

The problem given by (197)–(202) admits a unique solution up to a $\overline{\mathbf{Y}}$-constant function. Exploiting linearity, the solution is given as:

$$\begin{array}{cc}\hfill & {\mathbf{u}}_{\mathrm{I}}^{\left(1\right)}=A_{\mathrm{I}}{\nabla }_{\overline{\mathbf{X}}}\left({\mathbf{u}}^{\left(0\right)}\right)+{\mathbf{a}}_{\mathrm{I}}{p}^{\left(0\right)}+c_1\left(\overline{\mathbf{X}}\right),\hfill \end{array}$$

(203)

$$\begin{array}{cc}\hfill & {\mathbf{u}}_{\mathrm{II}}^{\left(1\right)}=A_{\mathrm{II}}{\nabla }_{\overline{\mathbf{X}}}\left({\mathbf{u}}^{\left(0\right)}\right)+{\mathbf{a}}_{\mathrm{II}}{p}^{\left(0\right)}{c}_2\left(\overline{\mathbf{X}}\right),\hfill \end{array}$$

(204)

where $c_1\left(\overline{\mathbf{X}}\right)$ and $c_2\left(\overline{\mathbf{X}}\right)$ are $\overline{\mathbf{Y}}$-constant functions. The third rank tensors $A_{\mathrm{I}}$ and $A_{\mathrm{II}}$ are the solutions of the pore-scale problems given by:

$$\begin{array}{cc}\hfill {\nabla }_{\overline{\mathbf{Y}}}·\left({\mathbb{C}}_{\mathrm{I}}{\nabla }_{\overline{\mathbf{Y}}}\left(A_{\mathrm{I}}\right)\right)+{\nabla }_{\overline{\mathbf{Y}}}·{\mathbb{C}}_{\mathrm{I}}=0& \mathrm{in}{\Omega }_{\mathrm{I}}^0,\hfill \end{array}$$

(205)

$$\begin{array}{cc}\hfill {\nabla }_{\overline{\mathbf{Y}}}·\left({\mathbb{C}}_{\mathrm{II}}{\nabla }_{\overline{\mathbf{Y}}}\left(A_{\mathrm{II}}\right)\right)+{\nabla }_{\overline{\mathbf{Y}}}·{\mathbb{C}}_{\mathrm{II}}=0& \mathrm{in}{\Omega }_{\mathrm{II}}^0,\hfill \end{array}$$

(206)

$$\begin{array}{cc}\hfill A_{\mathrm{I}}=A_{\mathrm{II}}& \mathrm{on}{\Gamma }_{\mathrm{III}}^0,\hfill \end{array}$$

(207)

$$\begin{array}{cc}\hfill {\mathbb{C}}_{\mathrm{I}}{\nabla }_{\overline{\mathbf{Y}}}\left(A_{\mathrm{I}}\right){\mathbf{n}}_{\mathrm{III}}^0-{\mathbb{C}}_{\mathrm{II}}{\nabla }_{\overline{\mathbf{Y}}}\left(A_{\mathrm{II}}\right){\mathbf{n}}_{\mathrm{III}}^0=({\mathbb{C}}_{\mathrm{II}}-{\mathbb{C}}_{\mathrm{I}}){\mathbf{n}}_{\mathrm{III}}^0& \mathrm{on}{\Gamma }_{\mathrm{III}}^0,\hfill \end{array}$$

(208)

$$\begin{array}{cc}\hfill \left({\mathbb{C}}_{\mathrm{I}}{\nabla }_{\overline{\mathbf{Y}}}\left(A_{\mathrm{I}}\right)\right){\mathbf{n}}_{\mathrm{I}}^0+{\mathbb{C}}_{\mathrm{I}}{\mathbf{n}}_{\mathrm{I}}^0=0& \mathrm{on}{\Gamma }_{\mathrm{I}}^0\hfill \end{array}$$

(209)

$$\begin{array}{cc}\hfill \left({\mathbb{C}}_{\mathrm{II}}{\nabla }_{\overline{\mathbf{Y}}}\left(A_{\mathrm{II}}\right)\right){\mathbf{n}}_{\mathrm{II}}^0+{\mathbb{C}}_{\mathrm{II}}{\mathbf{n}}_{\mathrm{II}}^0=0& \mathrm{on}{\Gamma }_{\mathrm{II}}^0,\hfill \end{array}$$

(210)

and the vectors ${\mathbf{a}}_{\mathrm{I}}$ and ${\mathbf{a}}_{\mathrm{II}}$ are the solution to this pore-scale problem:

$$\begin{array}{cc}\hfill {\nabla }_{\overline{\mathbf{Y}}}·\left({\mathbb{C}}_{\mathrm{I}}{\nabla }_{\overline{\mathbf{Y}}}\left({\mathbf{a}}_{\mathrm{I}}\right)\right)=0& \mathrm{in}{\Omega }_{\mathrm{I}}^0,\hfill \end{array}$$

(211)

$$\begin{array}{cc}\hfill {\nabla }_{\overline{\mathbf{Y}}}·\left({\mathbb{C}}_{\mathrm{II}}{\nabla }_{\overline{\mathbf{Y}}}\left({\mathbf{a}}_{\mathrm{II}}\right)\right)=0& \mathrm{in}{\Omega }_{\mathrm{II}}^0,\hfill \end{array}$$

(212)

$$\begin{array}{cc}\hfill {\mathbf{a}}_{\mathrm{I}}={\mathbf{a}}_{\mathrm{II}}& \mathrm{on}{\Gamma }_{\mathrm{III}}^0,\hfill \end{array}$$

(213)

$$\begin{array}{cc}\hfill \left({\mathbb{C}}_{\mathrm{I}}{\nabla }_{\overline{\mathbf{Y}}}\left({\mathbf{a}}_{\mathrm{I}}\right)\right){\mathbf{n}}_{\mathrm{III}}^0=\left({\mathbb{C}}_{\mathrm{II}}{\nabla }_{\overline{\mathbf{Y}}}\left({\mathbf{a}}_{\mathrm{II}}\right)\right){\mathbf{n}}_{\mathrm{III}}^0& \mathrm{on}{\Gamma }_{\mathrm{III}}^0,\hfill \end{array}$$

(214)

$$\begin{array}{cc}\hfill ({\mathbb{C}}_{\mathrm{I}}{\nabla }_{\overline{\mathbf{Y}}}\left({\mathbf{a}}_{\mathrm{I}}\right)+{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}){\mathbf{n}}_{\mathrm{I}}^0=0& \mathrm{on}{\Gamma }_{\mathrm{I}}^0,\hfill \end{array}$$

(215)

$$\begin{array}{cc}\hfill ({\mathbb{C}}_{\mathrm{II}}{\nabla }_{\overline{\mathbf{Y}}}\left({\mathbf{a}}_{\mathrm{II}}\right)+{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}){\mathbf{n}}_{\mathrm{II}}^0=0& \mathrm{on}{\Gamma }_{\mathrm{II}}^0.\hfill \end{array}$$

(216)

Both (205)–(210) and (211)–(216) are to be solved on the periodic cell. To ensure the uniqueness of the solution, we also require a further condition on $A_{\mathrm{I}}$, $A_{\mathrm{II}}$, ${\mathbf{a}}_{\mathrm{I}}$ and ${\mathbf{a}}_{\mathrm{II}}$, for example:

$${\langle A_{\mathrm{I}}\rangle }_{\mathrm{I}}+{\langle A_{\mathrm{II}}\rangle }_{\mathrm{II}}=0\mathrm{and}{\langle {\mathbf{a}}_{\mathrm{I}}\rangle }_{\mathrm{I}}+{\langle {\mathbf{a}}_{\mathrm{II}}\rangle }_{\mathrm{II}}=0.$$

(217)

Here, we wish to discuss in detail the cell problems (143)–(145), (205)–(210) and (211)–(216) and how they can potentially be solved. These pore-scale periodic cell problems are to be solved to determine the model coefficients of the final macroscale model. It is through these model coefficients that the complexity of the materials microstructure is encoded in the final model.

In general, with the asymptotic homogenization technique, these cell problems would only depend on the pore-scale and therefore can be solved in a straight-forward way. For example, solving the pore-scale asymptotic homogenization cell problems for linear elastic composites was carried out in [44], and for linear poroelasticity, the problems were solved in [40]. Similarly, it would be possible to solve the cell problems arising from linear poroelastic composites by combining the techniques used in both of these previous works. In the linear case, we have the problems (143)–(145), (205)–(210) and (211)–(216) with the simplification that ${\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}$ approaches the identity. This simplification means that the two scales are fully decoupled, and we can solve the fluid and the elastic-type cell problems.

However, due to the nonlinearity of the system we consider here, the two scales are coupled, meaning that the pore-scale periodic cell problems have a dependence on the macroscale and therefore cannot be easily solved. This dependence is through the quantity ${\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}$ appearing in (143)–(145) and (215)–(216). This quantity is the Piola transform, which involves the leading order deformation gradient ${\mathbf{F}}^{\left(0\right)}$. This depends on both the pore-scale and the macroscale, as can be seen in Equations (120) and (121). This means that the two scales are not fully decoupled, and therefore, this dramatically increases the computational complexity.

It is however crucial for a realistic analysis of the scenarios of interest (such as biological tissues) to be able to solve problems of this type. Despite the complexity, there are some potential emerging techniques that may mean it would be possible to solve this model numerically in the future. A recent example of a proposed method that could be potentially used to solve the types of problems arising in this work is found in [45]. This work investigated the potential of using Artificial Neural Networks (ANNs) for quick, accurate upscaling and localisation of problems. The method involves an incremental numerical approach where there is a rearrangement of the cell properties relating to the current deformation, and this means that there is a remodelling of the macroscopic model after each incremental time step. This method is applicable to finite strain and large deformation problems, whilst there will only be infinitesimal deformation within each incremental time step. Reference [45] investigated the full effects of the coupling between the macroscale and microscale for the first time in the analysis of fluid-saturated porous media. We believe that by following an approach similar to the one set out in [45], we could obtain a solution to our model numerically.

We can use our expressions (203) and (204) for ${\mathbf{u}}_{\mathrm{I}}^{\left(1\right)}$ and ${\mathbf{u}}_{\mathrm{II}}^{\left(1\right)}$ to rewrite ${\left({\sigma }_{\mathrm{I}}^{\left(0\right)}\right)}_{\mathrm{lin}}$ and ${\left({\sigma }_{\mathrm{II}}^{\left(0\right)}\right)}_{\mathrm{lin}}$. We have:

$$\begin{array}{cc}\hfill {\left({\sigma }_{\mathrm{I}}^{\left(0\right)}\right)}_{\mathrm{lin}}=& {\mathbb{C}}_{\mathrm{I}}{\nabla }_{\overline{\mathbf{Y}}}(A_{\mathrm{I}}{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}+{\mathbf{a}}_{\mathrm{I}}{p}^{\left(0\right)})+{\mathbb{C}}_{\mathrm{I}}{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}\hfill \\ \hfill =& {\mathbb{C}}_{\mathrm{I}}{\nabla }_{\overline{\mathbf{Y}}}{A}_{\mathrm{I}}{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}+{\mathbb{C}}_{\mathrm{I}}{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{a}}_{\mathrm{I}}{p}^{\left(0\right)}+{\mathbb{C}}_{\mathrm{I}}{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}\hfill \\ \hfill =& {\mathbb{C}}_{\mathrm{I}}{\mathbb{M}}_{\mathrm{I}}{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}+{\mathbb{C}}_{\mathrm{I}}{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}+{\mathbb{C}}_{\mathrm{I}}{\mathbf{Q}}_{\mathrm{I}}{p}^{\left(0\right)}\hfill \\ \hfill =& ({\mathbb{C}}_{\mathrm{I}}{\mathbb{M}}_{\mathrm{I}}+{\mathbb{C}}_{\mathrm{I}}){\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}+{\mathbb{C}}_{\mathrm{I}}{\mathbf{Q}}_{\mathrm{I}}{p}^{\left(0\right)},\hfill \end{array}$$

(218)

and:

$$\begin{array}{cc}\hfill {\left({\sigma }_{\mathrm{II}}^{\left(0\right)}\right)}_{\mathrm{lin}}=& {\mathbb{C}}_{\mathrm{II}}{\nabla }_{\overline{\mathbf{Y}}}(A_{\mathrm{II}}{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}+{\mathbf{a}}_{\mathrm{II}}{p}^{\left(0\right)})+{\mathbb{C}}_{\mathrm{II}}{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}\hfill \\ \hfill =& {\mathbb{C}}_{\mathrm{II}}{\nabla }_{\overline{\mathbf{Y}}}{A}_{\mathrm{II}}{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}+{\mathbb{C}}_{\mathrm{II}}{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{a}}_{\mathrm{II}}{p}^{\left(0\right)}+{\mathbb{C}}_{\mathrm{II}}{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}\hfill \\ \hfill =& {\mathbb{C}}_{\mathrm{II}}{\mathbb{M}}_{\mathrm{II}}{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}+{\mathbb{C}}_{\mathrm{II}}{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}+{\mathbb{C}}_{\mathrm{II}}{\mathbf{Q}}_{\mathrm{II}}{p}^{\left(0\right)}\hfill \\ \hfill =& ({\mathbb{C}}_{\mathrm{II}}{\mathbb{M}}_{\mathrm{II}}+{\mathbb{C}}_{\mathrm{II}}){\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}+{\mathbb{C}}_{\mathrm{II}}{\mathbf{Q}}_{\mathrm{II}}{p}^{\left(0\right)},\hfill \end{array}$$

(219)

where we define the pore-scale gradients of the auxiliary variables as:

$${\mathbb{M}}_{\mathrm{I}}={\nabla }_{\overline{\mathbf{Y}}}{A}_{\mathrm{I}},{\mathbb{M}}_{\mathrm{II}}={\nabla }_{\overline{\mathbf{Y}}}{A}_{\mathrm{II}},{\mathbf{Q}}_{\mathrm{I}}={\nabla }_{\overline{\mathbf{Y}}}{\mathbf{a}}_{\mathrm{I}},{\mathbf{Q}}_{\mathrm{II}}={\nabla }_{\overline{\mathbf{Y}}}{\mathbf{a}}_{\mathrm{II}}.$$

(220)

Then, we can return to (163) and use our linearised solid stresses to find the effective stress:

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{Eff}}=& {\langle {\left({\sigma }_{\mathrm{I}}^{\left(0\right)}\right)}_{\mathrm{lin}}\rangle }_{\mathrm{I}}+{\langle {\left({\sigma }_{\mathrm{II}}^{\left(0\right)}\right)}_{\mathrm{lin}}\rangle }_{\mathrm{II}}-p^{\left(0\right)}{\langle {\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}\rangle }_{\mathrm{f}}\hfill \\ \hfill =& \left({\langle {\mathbb{C}}_{\mathrm{I}}{\mathbb{M}}_{\mathrm{I}}+{\mathbb{C}}_{\mathrm{I}}\rangle }_{\mathrm{I}}+{\langle {\mathbb{C}}_{\mathrm{II}}{\mathbb{M}}_{\mathrm{II}}+{\mathbb{C}}_{\mathrm{II}}\rangle }_{\mathrm{II}}\right){\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}+({\langle {\mathbb{C}}_{\mathrm{I}}{\mathbf{Q}}_{\mathrm{I}}\rangle }_{\mathrm{I}}+{\langle {\mathbb{C}}_{\mathrm{II}}{\mathbf{Q}}_{\mathrm{II}}\rangle }_{\mathrm{II}}\hfill \\ \hfill & -{\langle {\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}\rangle }_{\mathrm{f}})p^{\left(0\right)}.\hfill \end{array}$$

(221)

As mentioned in Section 3.4, we return to the expression (157), restated here for convenience,

$$\begin{array}{cc}\hfill {\nabla }_{\overline{\mathbf{X}}}·\left({\langle {\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\mathbf{w}\rangle }_{\mathrm{f}}+{\langle {\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\rangle }_{\mathrm{f}}{\dot{\mathbf{u}}}^{\left(0\right)}\right)& =-\left({\langle {\nabla }_{\overline{\mathbf{X}}}·{\left({\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}\right)}^{\mathrm{T}}\rangle }_{\mathrm{I}}+{\langle {\nabla }_{\overline{\mathbf{X}}}·{\left({\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}\right)}^{\mathrm{T}}\rangle }_{\mathrm{II}}\right)·{\dot{\mathbf{u}}}^{\left(0\right)}\hfill \\ \hfill & +{\langle {\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}:{\nabla }_{\overline{\mathbf{Y}}}{\dot{\mathbf{u}}}_{\mathrm{I}}^{\left(1\right)}\rangle }_{\mathrm{I}}+{\langle {\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}:{\nabla }_{\overline{\mathbf{Y}}}{\dot{\mathbf{u}}}_{\mathrm{II}}^{\left(1\right)}\rangle }_{\mathrm{II}}\hfill \end{array}$$

(222)

We obtain expressions for ${\dot{\mathbf{u}}}_{\mathrm{I}}^{\left(1\right)}$ and ${\dot{\mathbf{u}}}_{\mathrm{II}}^{\left(1\right)}$ by taking the time derivative of (203) and (204), and we then substitute these expressions into (222) to obtain:

$$\begin{array}{cc}\hfill {\nabla }_{\overline{\mathbf{X}}}·\left({\langle {\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\mathbf{w}\rangle }_{\mathrm{f}}+{\langle {\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\rangle }_{\mathrm{f}}{\dot{\mathbf{u}}}^{\left(0\right)}\right)& ={\langle {\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}:{\nabla }_{\overline{\mathbf{Y}}}{A}_{\mathrm{I}}\rangle }_{\mathrm{I}}:{\nabla }_{\overline{\mathbf{X}}}{\dot{\mathbf{u}}}^{\left(0\right)}+{\langle {\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}:{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{a}}_{\mathrm{I}}\rangle }_{\mathrm{I}}{\dot{p}}^{\left(0\right)}\hfill \\ \hfill & +{\langle {\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}:{\nabla }_{\overline{\mathbf{Y}}}{A}_{\mathrm{II}}\rangle }_{\mathrm{II}}:{\nabla }_{\overline{\mathbf{X}}}{\dot{\mathbf{u}}}^{\left(0\right)}+{\langle {\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}:{\nabla }_{\overline{\mathbf{Y}}}{\mathbf{a}}_{\mathrm{II}}\rangle }_{\mathrm{II}}{\dot{p}}^{\left(0\right)}\hfill \\ \hfill & -\left({\langle {\nabla }_{\overline{\mathbf{X}}}·{\left({\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}\right)}^{\mathrm{T}}\rangle }_{\mathrm{I}}+{\langle {\nabla }_{\overline{\mathbf{X}}}·{\left({\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}\right)}^{\mathrm{T}}\rangle }_{\mathrm{II}}\right)·{\dot{\mathbf{u}}}^{\left(0\right)}.\hfill \end{array}$$

(223)

Expanding the LHS in (223) and using (220), we obtain:

$$\begin{array}{cc}\hfill {\nabla }_{\overline{\mathbf{X}}}·{\langle {\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\mathbf{w}\rangle }_{\mathrm{f}}+{\nabla }_{\overline{\mathbf{X}}}·\left({\langle {\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\rangle }_{\mathrm{f}}{\dot{\mathbf{u}}}^{\left(0\right)}\right)=& {\langle {\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}:{\mathbb{M}}_{\mathrm{I}}\rangle }_{\mathrm{I}}:{\nabla }_{\overline{\mathbf{X}}}{\dot{\mathbf{u}}}^{\left(0\right)}+{\langle {\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}:{\mathbf{Q}}_{\mathrm{I}}\rangle }_{\mathrm{I}}{\dot{p}}^{\left(0\right)}\hfill \\ \hfill & +{\langle {\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}:{\mathbb{M}}_{\mathrm{II}}\rangle }_{\mathrm{II}}:{\nabla }_{\overline{\mathbf{X}}}{\dot{\mathbf{u}}}^{\left(0\right)}+{\langle {\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}:{\mathbf{Q}}_{\mathrm{II}}\rangle }_{\mathrm{II}}{\dot{p}}^{\left(0\right)}\hfill \\ \hfill & -\left({\langle {\nabla }_{\overline{\mathbf{X}}}·{\left({\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}\right)}^{\mathrm{T}}\rangle }_{\mathrm{I}}+{\langle {\nabla }_{\overline{\mathbf{X}}}·{\left({\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}\right)}^{\mathrm{T}}\rangle }_{\mathrm{II}}\right)·{\dot{\mathbf{u}}}^{\left(0\right)}.\hfill \end{array}$$

(224)

Expanding the second term on the LHS further and rearranging, we obtain:

$$\begin{array}{cc}\hfill {\nabla }_{\overline{\mathbf{X}}}·{\langle {\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\mathbf{w}\rangle }_{\mathrm{f}}+{\langle {\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\rangle }_{\mathrm{f}}:{\nabla }_{\overline{\mathbf{X}}}{\dot{\mathbf{u}}}^{\left(0\right)}=& {\langle {\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}:{\mathbb{M}}_{\mathrm{I}}\rangle }_{\mathrm{I}}:{\nabla }_{\overline{\mathbf{X}}}{\dot{\mathbf{u}}}^{\left(0\right)}+{\langle {\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}:{\mathbf{Q}}_{\mathrm{I}}\rangle }_{\mathrm{I}}{\dot{p}}^{\left(0\right)}\hfill \\ \hfill +& {\langle {\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}:{\mathbb{M}}_{\mathrm{II}}\rangle }_{\mathrm{II}}:{\nabla }_{\overline{\mathbf{X}}}{\dot{\mathbf{u}}}^{\left(0\right)}+{\langle {\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}:{\mathbf{Q}}_{\mathrm{II}}\rangle }_{\mathrm{II}}{\dot{p}}^{\left(0\right)}-({\langle {\nabla }_{\overline{\mathbf{X}}}·{\left({\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}\right)}^{\mathrm{T}}\rangle }_{\mathrm{I}}\hfill \\ \hfill +& {\langle {\nabla }_{\overline{\mathbf{X}}}·{\left({\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}\right)}^{\mathrm{T}}\rangle }_{\mathrm{II}}+{\langle {\nabla }_{\overline{\mathbf{X}}}·{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}\rangle }_{\mathrm{f}})·{\dot{\mathbf{u}}}^{\left(0\right)}.\hfill \end{array}$$

(225)

Rearranging and collecting ${\dot{p}}^{\left(0\right)}$ terms gives:

$$\begin{array}{cc}\hfill {\dot{p}}^{\left(0\right)}& =\frac{1}{{\langle {\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}:{\mathbf{Q}}_{\mathrm{I}}\rangle }_{\mathrm{I}}+{\langle {\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}:{\mathbf{Q}}_{\mathrm{II}}\rangle }_{\mathrm{II}}}[{\nabla }_{\overline{\mathbf{X}}}·{\langle {\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\mathbf{w}\rangle }_{\mathrm{f}}+({\langle {\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\rangle }_{\mathrm{f}}-{\langle {\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}:{\mathbb{M}}_{\mathrm{I}}\rangle }_{\mathrm{I}}\hfill \\ \hfill & -{\langle {\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}:{\mathbb{M}}_{\mathrm{II}}\rangle }_{\mathrm{II}}):{\nabla }_{\overline{\mathbf{X}}}{\dot{\mathbf{u}}}^{\left(0\right)}+({\langle {\nabla }_{\overline{\mathbf{X}}}·{\left({\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}\right)}^{\mathrm{T}}\rangle }_{\mathrm{I}}+{\langle {\nabla }_{\overline{\mathbf{X}}}·{\left({\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}\right)}^{\mathrm{T}}\rangle }_{\mathrm{II}}\hfill \\ \hfill & +{\langle {\nabla }_{\overline{\mathbf{X}}}·{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}\rangle }_{\mathrm{f}})·{\dot{\mathbf{u}}}^{\left(0\right)}].\hfill \end{array}$$

(226)

We define:

$$M:=\frac{-1}{{\langle {\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}:{\mathbf{Q}}_{\mathrm{I}}\rangle }_{\mathrm{I}}+{\langle {\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}:{\mathbf{Q}}_{\mathrm{II}}\rangle }_{\mathrm{II}}},$$

(227)

and:

$$\tilde{\alpha }:={\langle {\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\rangle }_{\mathrm{f}}-{\langle {\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}:{\mathbb{M}}_{\mathrm{I}}\rangle }_{\mathrm{I}}-{\langle {\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}:{\mathbb{M}}_{\mathrm{II}}\rangle }_{\mathrm{II}},$$

(228)

and we can then use (227) and (228) to rewrite (226) as:

$$\begin{array}{cc}\hfill {\dot{p}}^{\left(0\right)}=& -M[{\nabla }_{\overline{\mathbf{X}}}·{\langle {\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\mathbf{w}\rangle }_{\mathrm{f}}+\tilde{\alpha }:{\nabla }_{\overline{\mathbf{X}}}{\dot{\mathbf{u}}}^{\left(0\right)}+({\langle {\nabla }_{\overline{\mathbf{X}}}·{\left({\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}\right)}^{\mathrm{T}}\rangle }_{\mathrm{I}}+{\langle {\nabla }_{\overline{\mathbf{X}}}·{\left({\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}\right)}^{\mathrm{T}}\rangle }_{\mathrm{II}}\hfill \\ \hfill & +{\langle {\nabla }_{\overline{\mathbf{X}}}·{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}\rangle }_{\mathrm{f}})·{\dot{\mathbf{u}}}^{\left(0\right)}].\hfill \end{array}$$

(229)

Finally, we can divide by M to obtain:

$$\begin{array}{cc}\hfill \frac{{\dot{p}}^{\left(0\right)}}{M}=& -{\nabla }_{\overline{\mathbf{X}}}·{\langle {\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\mathbf{w}\rangle }_{\mathrm{f}}-\tilde{\alpha }:{\nabla }_{\overline{\mathbf{X}}}{\dot{\mathbf{u}}}^{\left(0\right)}-({\langle {\nabla }_{\overline{\mathbf{X}}}·{\left({\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}\right)}^{\mathrm{T}}\rangle }_{\mathrm{I}}+{\langle {\nabla }_{\overline{\mathbf{X}}}·{\left({\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}\right)}^{\mathrm{T}}\rangle }_{\mathrm{II}}\hfill \\ \hfill & +{\langle {\nabla }_{\overline{\mathbf{X}}}·{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}\rangle }_{\mathrm{f}})·{\dot{\mathbf{u}}}^{\left(0\right)}.\hfill \end{array}$$

(230)

We therefore have now derived the effective macroscale governing equations for a nonlinear poroelastic composite that has the constitutive law given by the de Saint-Venant strain energy function. We state our novel macroscale model and then consider limit cases for the model where we obtain previously known results from the literature. The equations in the macroscale model represent a poroelastic-type system of PDEs. Therefore, the macroscale model is given by:

$$\begin{array}{cc}\hfill & {\langle \mathbf{w}\rangle }_{\mathrm{f}}=-{\langle \mathbf{W}\rangle }_{\mathrm{f}}{\nabla }_{\overline{\mathbf{X}}}{p}^{\left(0\right)},\hfill \end{array}$$

(231)

$$\begin{array}{cc}\hfill & {\nabla }_{\overline{\mathbf{X}}}·{\sigma }_{\mathrm{Eff}}=0,\hfill \end{array}$$

(232)

$$\begin{array}{cc}\hfill & {\sigma }_{\mathrm{Eff}}=\left({\langle {\mathbb{C}}_{\mathrm{I}}{\mathbb{M}}_{\mathrm{I}}+{\mathbb{C}}_{\mathrm{I}}\rangle }_{\mathrm{I}}+{\langle {\mathbb{C}}_{\mathrm{II}}{\mathbb{M}}_{\mathrm{II}}+{\mathbb{C}}_{\mathrm{II}}\rangle }_{\mathrm{II}}\right):{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}+({\langle {\mathbb{C}}_{\mathrm{I}}{\mathbf{Q}}_{\mathrm{I}}\rangle }_{\mathrm{I}}+{\langle {\mathbb{C}}_{\mathrm{II}}{\mathbf{Q}}_{\mathrm{II}}\rangle }_{\mathrm{II}}\hfill \\ \hfill & -{\langle {\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}\rangle }_{\mathrm{f}})p^{\left(0\right)},\hfill \end{array}$$

(233)

$$\begin{array}{cc}\hfill & {\dot{p}}^{\left(0\right)}=\frac{1}{{\langle {\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}:{\mathbf{Q}}_{\mathrm{I}}\rangle }_{\mathrm{I}}+{\langle {\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}:{\mathbf{Q}}_{\mathrm{II}}\rangle }_{\mathrm{II}}}[{\nabla }_{\overline{\mathbf{X}}}·{\langle {\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\mathbf{w}\rangle }_{\mathrm{f}}+({\langle {\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\rangle }_{\mathrm{f}}-{\langle {\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}:{\mathbb{M}}_{\mathrm{I}}\rangle }_{\mathrm{I}}\hfill \\ \hfill & -{\langle {\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}:{\mathbb{M}}_{\mathrm{II}}\rangle }_{\mathrm{II}}):{\nabla }_{\overline{\mathbf{X}}}{\dot{\mathbf{u}}}^{\left(0\right)}+({\langle {\nabla }_{\overline{\mathbf{X}}}·{\left({\mathsf{G}}_{\mathrm{I}}^{\left(0\right)}\right)}^{\mathrm{T}}\rangle }_{\mathrm{I}}+{\langle {\nabla }_{\overline{\mathbf{X}}}·{\left({\mathsf{G}}_{\mathrm{II}}^{\left(0\right)}\right)}^{\mathrm{T}}\rangle }_{\mathrm{II}}\hfill \\ \hfill & +{\langle {\nabla }_{\overline{\mathbf{X}}}·{\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}\rangle }_{\mathrm{f}})·{\dot{\mathbf{u}}}^{\left(0\right)}],\hfill \end{array}$$

(234)

where the macroscale pressure is denoted by $p^{\left(0\right)}$, ${\mathbf{u}}^{\left(0\right)}$ is the leading order elastic displacement, the solid velocity is represented by ${\dot{\mathbf{u}}}^{\left(0\right)}$ and, finally, $\mathbf{w}$ is the leading order relative fluid–solid velocity. Equation (231) in the macroscale model is Darcy’s law for $\mathbf{w}$. Equation (232) of the macroscale model is the stress balance equation with the effective stress tensor ${\sigma }_{\mathrm{Eff}}$. The constitutive equation given by (233) is of poroelastic type; however, it contains a nonlinear coordinate transformation. The effective elasticity tensor for the material is given by:

$$\tilde{\mathbb{C}}={\langle {\mathbb{C}}_{\mathrm{I}}{\mathbb{M}}_{\mathrm{I}}+{\mathbb{C}}_{\mathrm{I}}\rangle }_{\mathrm{I}}+{\langle {\mathbb{C}}_{\mathrm{II}}{\mathbb{M}}_{\mathrm{II}}+{\mathbb{C}}_{\mathrm{II}}\rangle }_{\mathrm{II}}$$

(235)

Finally, (234) is the conservation of mass for a nonlinear poroelastic composite material. The three terms in our expression for ${\dot{p}}^{\left(0\right)}$ that are the divergence of the Piola transforms in the subphase, matrix and fluid, respectively, describe the volume changes related to the deformation and can be viewed as a correction term that maintains the conservation of mass despite the nonlinear coordinate transformations. We can say that the effective mechanical behaviour of our nonlinear poroelastic composite material can be fully described by the model coefficients, that is the effective elasticity tensor $\tilde{\mathbb{C}}$, the hydraulic conductivity tensor ${\langle \mathbf{W}\rangle }_{\mathrm{f}}$, the transformed Biot’s tensor of coefficients, $\tilde{\alpha }$ and the transformed Biot’s modulus, M. We note that although, structurally, this model is similar to that of linear poroelasticity, the key novelty resides in the model coefficients that capture the nonlinear deformations and the additional terms. These model coefficients are to be obtained by solving the novel cell problems (143)–(145), (205)–(210) and (211)–(216)

Next, we consider particular cases for our model and are able to derive previously known models that were developed using the asymptotic homogenization technique.

4.2. Comparison with Linear Poroelastic Composites

We begin with the case where the poroelastic composite solids that we are considering are linearly elastic. This setting is applicable to many situations including the interstitial matrix of biological tissues and when describing hard hierarchical materials such as bones and tendons where deformations are very small. We note that to reduce to the linear elastic case, we have that ${\mathsf{G}}^{\left(0\right)}\to \mathbf{I}$. We also introduce the notation:

$${\langle {\phi }_{\mathrm{I}}\rangle }_{\mathrm{I}}+{\langle {\phi }_{\mathrm{II}}\rangle }_{\mathrm{II}}={\langle {\phi }_{\mathrm{I}}+{\phi }_{\mathrm{II}}\rangle }_{\mathrm{s}},$$

(236)

for fields $\phi$ with components ${\phi }_{\mathrm{I}}$ and ${\phi }_{\mathrm{II}}$ defined in the solid cell portions ${\Omega }_{\mathrm{I}}^0$ or ${\Omega }_{\mathrm{II}}^0$, respectively. This means that the model (231)–(234) reduces to:

$$\left\{\begin{array}{c}\begin{array}{c}{\langle \mathbf{w}\rangle }_{\mathrm{f}}=-{\langle \mathbf{W}\rangle }_{\mathrm{f}}{\nabla }_{\overline{\mathbf{X}}}{p}^{\left(0\right)},\hfill \\ {\nabla }_{\overline{\mathbf{X}}}·{\sigma }_{\mathrm{Eff}}=0,\hfill \\ {\sigma }_{\mathrm{Eff}}={\langle {\mathbb{C}}_{\mathrm{I}}{\mathbb{M}}_{\mathrm{I}}+{\mathbb{C}}_{\mathrm{I}}+{\mathbb{C}}_{\mathrm{II}}{\mathbb{M}}_{\mathrm{II}}+{\mathbb{C}}_{\mathrm{II}}\rangle }_s{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}+({\langle {\mathbb{C}}_{\mathrm{I}}{\mathbf{Q}}_{\mathrm{I}}+{\mathbb{C}}_{\mathrm{II}}{\mathbf{Q}}_{\mathrm{II}}\rangle }_s-\varphi \mathbf{I})p^{\left(0\right)},\hfill \\ {\dot{p}}^{\left(0\right)}=\frac{1}{{\langle Tr({\mathbf{Q}}_{\mathrm{I}}+{\mathbf{Q}}_{\mathrm{II}})\rangle }_s}\left({\nabla }_{\overline{\mathbf{X}}}·{\langle \mathbf{w}\rangle }_{\mathrm{f}}+(\varphi \mathbf{I}-{\langle Tr({\mathbb{M}}_{\mathrm{I}}+{\mathbb{M}}_{\mathrm{II}})\rangle }_s):{\nabla }_{\overline{\mathbf{X}}}{\dot{\mathbf{u}}}^{\left(0\right)}\right),\hfill \end{array}\hfill \end{array}\right.$$

(237)

This is identically the model for linear poroelastic composites presented in [27]. We note that the integral average over the fluid domain of the identity tensor, ${\langle \mathbf{I}\rangle }_{\mathrm{f}}$, is $\varphi \mathbf{I}$. We also have that ${\langle {\mathsf{G}}_{\beta }^{\left(0\right)}:{\mathbf{Q}}_{\beta }\rangle }_{\beta }$, where $\beta =\mathrm{I},\mathrm{II}$ becomes ${\langle \mathrm{Tr}{\mathbf{Q}}_{\beta }\rangle }_{\beta }$ in the limit ${\mathsf{G}}^{\left(0\right)}\to \mathbf{I}$. Similarly, ${\langle {\mathsf{G}}_{\beta }^{\left(0\right)}:{\mathbb{M}}_{\beta }\rangle }_{\beta }:{\nabla }_{\overline{\mathbf{X}}}{\dot{\mathbf{u}}}^{\left(0\right)}$ becomes ${\langle \mathrm{Tr}{\mathbb{M}}_{\beta }\rangle }_{\beta }:{\nabla }_{\overline{\mathbf{X}}}{\dot{\mathbf{u}}}^{\left(0\right)}$. As a result of taking the limit ${\mathsf{G}}^{\left(0\right)}\to \mathbf{I}$, we also have that ${\mathbf{F}}^{\left(0\right)}\to \mathbf{I}$, and using this in the cell problems (143)–(145), (205)–(210) and (211)–(216), we recover the same cell problems as in [27] and therefore obtain exactly the same model coefficients under this limit.

We also note, as remarked in [27], that if we consider only one elastic phase, then the model reduces to the macroscale model for a standard poroelastic material (see the no growth limit in [39], as well as [22,25]). Furthermore, in the limit of zero fluid (no pores), then this macroscale model reduces to the model for an elastic composite [46].

4.3. Comparison with Nonlinear Poroelasticity

We now wish to recover previous work on nonlinear poroelasticity. For this, we assume that our material has only one hyperelastic phase, the matrix that we denote by ${\Omega }_{\mathrm{s}}$, with fluid flowing in the pores. We wish to compare our Equation (225) under the assumption of only one elastic phase with the generalised Biot fluid equation found in [31]. We can rewrite our equation with only one elastic phase as:

$$\begin{array}{cc}\hfill -{\nabla }_{\overline{\mathbf{X}}}·{\langle {\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\mathbf{W}\rangle }_{\mathrm{f}}{\nabla }_{\overline{\mathbf{X}}}{p}^{\left(0\right)}& ={\langle {\mathsf{G}}_{\mathrm{s}}^{\left(0\right)}:{\mathbb{M}}_{\mathrm{s}}\rangle }_{\mathrm{s}}:{\nabla }_{\overline{\mathbf{X}}}{\dot{\mathbf{u}}}^{\left(0\right)}+{\langle {\mathsf{G}}_{\mathrm{s}}^{\left(0\right)}:{\mathbf{Q}}_{\mathrm{s}}\rangle }_{\mathrm{s}}{\dot{p}}^{\left(0\right)}-{\langle {\nabla }_{\overline{\mathbf{X}}}·{\left({\mathsf{G}}_{\mathrm{s}}^{\left(0\right)}\right)}^{\mathrm{T}}\rangle }_{\mathrm{s}}·{\dot{\mathbf{u}}}^{\left(0\right)}\hfill \\ \hfill & -{\nabla }_{\overline{\mathbf{X}}}·\left({\langle {\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\rangle }_{\mathrm{f}}{\dot{\mathbf{u}}}^{\left(0\right)}\right).\hfill \end{array}$$

(238)

We note that we do not consider any body forces, so the f appearing in [31] equals zero here. A general ansatz for ${\mathbf{u}}^{\left(1\right)}$ was also used in [31], but we are still able to make some identifications between the terms in our Equation (238) and the terms in their generalised Biot fluid equation. We note that the model coefficients of [31] all involve ${\left({\mathsf{G}}^{\left(0\right)}\right)}^{\mathrm{T}}$, so we should point out that due to a different choice in the definition of the Piola transformation between this work and [31], we have that ${\left({\mathsf{G}}^{\left(0\right)}\right)}^{\mathrm{T}}$ in [31] equals ${\mathsf{G}}^{\left(0\right)}$ in this work. We also use a modified (16) where ${\mathsf{G}}^{\left(0\right)}={\mathsf{G}}_{\mathrm{s}}^{\left(0\right)}$ in ${\Omega }_{\mathrm{s}}^0$ and ${\mathsf{G}}^{\left(0\right)}={\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}$ in ${\Omega }_{\mathrm{f}}^0$, and we have the continuity of ${\mathsf{G}}_{\mathrm{s}}^{\left(0\right)}$ and ${\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}$ on the interface between the phases. This means that we can identify the coefficients:

$$\begin{array}{c}{K}^{*}=-{\langle {\mathsf{G}}^{\left(0\right)}\mathbf{W}\rangle }_{\mathrm{f}},A^{*}={\langle {\nabla }_{\overline{\mathbf{Y}}}\left({\mathsf{G}}_{\mathrm{s}}^{\left(0\right)}A\right)\rangle }_{\mathrm{s}}={\langle {\mathsf{G}}_{\mathrm{s}}^{\left(0\right)}:{\mathbb{M}}_{\mathrm{s}}\rangle }_{\mathrm{s}},\\ B^{*}={\langle {\nabla }_{\overline{\mathbf{Y}}}\left({\mathsf{G}}^{\left(0\right)}\mathbf{a}\right)\rangle }_{\mathrm{s}}={\langle {\mathsf{G}}_{\mathrm{s}}^{\left(0\right)}:{\mathbf{Q}}_{\mathrm{s}}\rangle }_{\mathrm{s}},D^{*}={\langle {\mathsf{G}}^{\left(0\right)}\rangle }_{\Omega }.\end{array}$$

(239)

The difference in sign between $K^{*}$ and $\mathbf{W}$ is due to the difference in the choice of ansatz for the fluid problem. We therefore can recover the generalised Biot fluid equation from our Equation (238).

We now wish to compare the macroscale elasticity equation in [31] with our macroscale balance equation and constitutive law (232) and (233). We first should note that in our work, we are not considering any body forces, and therefore, we can assume that the f and b appearing in [31] are both zero. Using (233) in (232) and reducing to only one elastic phase, we obtain:

$${\nabla }_{\overline{\mathbf{X}}}·\left({\langle \mathbb{C}\mathbb{M}+\mathbb{C}\rangle }_{\mathrm{s}}:{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}+\left({\langle \mathbb{C}\mathbf{Q}\rangle }_{\mathrm{s}}-{\langle {\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}\rangle }_{\mathrm{f}}\right)p^{\left(0\right)}\right)=0$$

(240)

We are able to make the identifications between the equation of [31] using a generalised ansatz and our Equation (240), where a specific ansatz is used. The term $C{e}_x\left(u_0\right)$ can be identified with our term $\mathbb{C}:{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}$. Similarly, the term $C{e}_y\left(N(p_0,{\nabla }_x{u}_0)\right)$ can lead to our terms $\mathbb{C}\mathbb{M}:{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}+\mathbb{C}\mathbf{Q}{p}^{\left(0\right)}$ when using our ansatz. Finally, we need to compare the term ${\varphi }_1(p_0,{\nabla }_x{u}_0){\nabla }_x{p}_0$ appearing in [31] with the final term in (240). The final term in (240) can be expanded as:

$$-{\nabla }_{\overline{\mathbf{X}}}·\left({\langle {\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}\rangle }_{\mathrm{f}}{p}^{\left(0\right)}\right)=-{\langle {\left({\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\right)}^{\mathrm{T}}\rangle }_{\mathrm{f}}{\nabla }_{\overline{\mathbf{X}}}{p}^{\left(0\right)}-({\nabla }_{\overline{\mathbf{X}}}·{\langle {\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}\rangle }_{\mathrm{f}})p^{\left(0\right)}$$

(241)

The term ${\varphi }_1(p_0,{\nabla }_x{u}_0){\nabla }_x{p}_0$ can be identified with the first term in (241), and when we assume that ${\mathsf{G}}_{\mathrm{f}}^{\left(0\right)}$ is a constant, then we recover exactly the macroscale elasticity equation of [31].

We also wish to make a comparison between the cell problems found in [31] and our cell problems. The fluid cell problem (143)–(145) matches exactly the first cell problem found in [31]. We however do not require the second cell problem found in [31] as we do not include any forces in our formulation, so $f=0$. We also have the two elastic cell problems (205)–(210) and (211)–(216), and these cannot be directly compared with the cell problems found in [31] as these arise after the application of a specific ansatz. We can however compare the elastic problem of [31] with (165) and (169), and it is clear that if (165) and (169) are identical, up to a choice of sign, to the balance equation of [31]’ and continuity of stresses, then the same ansatz would produce the same cell problems.

4.4. Comparison with Nonlinear Elastic Composites

We now consider the case that our structure has no pores and therefore can be described as a composite comprised of two hyperelastic materials. We can instantly reduce our model (231) by removing the equations that govern the fluid. We are able to obtain the model derived in [33], where we make the assumption that plastic distortions are absent, that is assuming ${\mathbf{F}}_p=\mathbf{I}$ in [33].

When we assume that ${\mathbf{F}}_p=\mathbf{I}$ in [33], then we have that the plastic Green–Lagrange strain tensor ${\mathbf{E}}_p^{\left(0\right)}=\frac{1}{2}({\left({\mathbf{F}}_p^{\left(0\right)}\right)}^{\mathrm{T}}{\mathbf{F}}_p^{\left(0\right)}-\mathbf{I})=\frac{1}{2}(\mathbf{I}-\mathbf{I})=0$. This means the first Piola stress obtained in [33] becomes:

$$T_{\mathrm{lin}}^{\left(0\right)}={\mathcal{C}}_{\mathrm{R}}:\mathbf{H},$$

(242)

and we can make the identifications in the notation that ${\mathcal{C}}_{\mathrm{R}}=\mathbb{C}$ and $T_{\mathrm{lin}}^{\left(0\right)}={\sigma }_{\mathrm{lin}}^{\left(0\right)}$ in our case. Therefore, the first Piola stress tensor (194), that we obtained here matches the first Piola stress obtained by [33].

We should note that within this work, we use notation to specifically identify the two constituents of the composite, whereas [33] kept the different constituents implicit. We can modify our problem (197)–(202) to remove the involvement of the fluid. We therefore end up with the problem for elastic composites, which is given by:

$$\begin{array}{cc}\hfill {\nabla }_{\overline{\mathbf{Y}}}·\left({\mathbb{C}}_{\mathrm{I}}{\nabla }_{\overline{\mathbf{Y}}}\left({\mathbf{u}}_{\mathrm{I}}^{\left(1\right)}\right)\right)+{\nabla }_{\overline{\mathbf{Y}}}·\left({\mathbb{C}}_{\mathrm{I}}{\nabla }_{\overline{\mathbf{X}}}\left({\mathbf{u}}^{\left(0\right)}\right)\right)=0& \mathrm{in}{\Omega }_{\mathrm{I}}^0,\hfill \end{array}$$

(243)

$$\begin{array}{cc}\hfill {\nabla }_{\overline{\mathbf{Y}}}·\left({\mathbb{C}}_{\mathrm{II}}{\nabla }_{\overline{\mathbf{Y}}}\left({\mathbf{u}}_{\mathrm{II}}^{\left(1\right)}\right)\right)+{\nabla }_{\overline{\mathbf{Y}}}·\left({\mathbb{C}}_{\mathrm{II}}{\nabla }_{\overline{\mathbf{X}}}\left({\mathbf{u}}^{\left(0\right)}\right)\right)=0& \mathrm{in}{\Omega }_{\mathrm{II}}^0,\hfill \end{array}$$

(244)

$$\begin{array}{cc}\hfill {\mathbf{u}}_{\mathrm{I}}^{\left(1\right)}={\mathbf{u}}_{\mathrm{II}}^{\left(1\right)}& \mathrm{on}{\Gamma }_{\mathrm{III}}^0,\hfill \end{array}$$

(245)

$$\begin{array}{cc}\hfill {\mathbb{C}}_{\mathrm{I}}{\nabla }_{\overline{\mathbf{Y}}}\left({\mathbf{u}}_{\mathrm{I}}^{\left(1\right)}\right){\mathbf{n}}_{\mathrm{III}}^0-{\mathbb{C}}_{\mathrm{II}}{\nabla }_{\overline{\mathbf{Y}}}\left({\mathbf{u}}_{\mathrm{II}}^{\left(1\right)}\right){\mathbf{n}}_{\mathrm{III}}^0=({\mathbb{C}}_{\mathrm{II}}-{\mathbb{C}}_{\mathrm{I}}){\nabla }_{\overline{\mathbf{X}}}\left({\mathbf{u}}^{\left(0\right)}\right){\mathbf{n}}_{\mathrm{III}}^0& \mathrm{on}{\Gamma }_{\mathrm{III}}^0.\hfill \end{array}$$

(246)

When using the notation of [33], we can write (243)–(246) as:

$$\begin{array}{c}\hfill {\mathrm{Div}}_{\overline{Y}}{T}^{\left(0\right)}=\mathbf{0},\end{array}$$

(247)

$$\begin{array}{c}\hfill \left[\right[{\mathbf{u}}^{\left(1\right)}\left]\right]=\mathbf{0},\end{array}$$

(248)

$$\begin{array}{c}\hfill \left[\right[T^{\left(0\right)}·\mathbf{N}\left]\right]=\mathbf{0}.\end{array}$$

(249)

This matches identically the problem given by [33]. The reduced problem (243)–(246) has the ansatz:

$${\mathbf{u}}_{\mathrm{I}}^{\left(1\right)}={\mathbf{A}}_{\mathrm{I}}{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}\mathrm{and}{\mathbf{u}}_{\mathrm{II}}^{\left(1\right)}={\mathbf{A}}_{\mathrm{II}}{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)},$$

(250)

where $A_{\mathrm{I}}$ and $A_{\mathrm{II}}$ are third rank tensors. This is the ansatz (203) and (204) where we take $p^{\left(0\right)}=0$. This leads to the cell problem for $A_{\mathrm{I}}$ and $A_{\mathrm{II}}$:

$$\left\{\begin{array}{c}\begin{array}{c}{\nabla }_{\overline{\mathbf{Y}}}·\left({\mathbb{C}}_{\mathrm{I}}{\nabla }_{\overline{\mathbf{Y}}}\left(A_{\mathrm{I}}\right)\right)+{\nabla }_{\overline{\mathbf{Y}}}·{\mathbb{C}}_{\mathrm{I}}=0,\hfill \\ {\nabla }_{\overline{\mathbf{Y}}}·\left({\mathbb{C}}_{\mathrm{II}}{\nabla }_{\overline{\mathbf{Y}}}\left(A_{\mathrm{II}}\right)\right)+{\nabla }_{\overline{\mathbf{Y}}}·{\mathbb{C}}_{\mathrm{II}}=0,\hfill \\ A_{\mathrm{I}}=A_{\mathrm{II}},\hfill \\ {\mathbb{C}}_{\mathrm{I}}{\nabla }_{\overline{\mathbf{Y}}}\left(A_{\mathrm{I}}\right){\mathbf{n}}_{\mathrm{III}}^0-{\mathbb{C}}_{\mathrm{II}}{\nabla }_{\overline{\mathbf{Y}}}\left(A_{\mathrm{II}}\right){\mathbf{n}}_{\mathrm{III}}^0=({\mathbb{C}}_{\mathrm{II}}-{\mathbb{C}}_{\mathrm{I}}){\mathbf{n}}_{\mathrm{III}}^0,\hfill \end{array}\hfill \end{array}\right.$$

(251)

which is again cell problem (205)–(210) reduced under the assumption that our material has no pores. In [33], they remarked about the case of no plastic distortions occurring and stated the cell problem they obtained under those circumstances, and this is identical to the cell problem (251) where we make the identification that $\mathit{T}{\mathrm{Grad}}_{\overline{Y}}\xi ={\nabla }_{\overline{\mathbf{Y}}}A=\mathbb{M}$ in our work and where we use the implicit notation that:

$$\mathbb{M}=\left\{\begin{array}{c}\begin{array}{cc}\hfill {\mathbb{M}}_{\mathrm{I}}& \mathrm{in}{\Omega }_{\mathrm{I}}^0,\hfill \\ \hfill {\mathbb{M}}_{\mathrm{II}}& \mathrm{in}{\Omega }_{\mathrm{II}}^0.\hfill \end{array}\hfill \end{array}\right.$$

(252)

Therefore, the cell problem from [33] that matches (251) is given by:

$$\left\{\begin{array}{c}\begin{array}{c}{\mathrm{Div}}_{\overline{Y}}(\mathcal{C}+\mathcal{C}:\mathit{T}{\mathrm{Grad}}_{\overline{Y}}\xi )=\mathbf{0},\hfill \\ \left[\right[\xi \left]\right]=\mathbf{0},\hfill \\ \left[[(\mathcal{C}+\mathcal{C}:\mathit{T}{\mathrm{Grad}}_{\overline{Y}}\xi )·\mathbf{N}]\right]=\mathbf{0}.\hfill \end{array}\hfill \end{array}\right.$$

(253)

Finally, under the assumption of no fluid-filled pores, our macroscale model (231) can be reduced to:

$$\left\{\begin{array}{c}\begin{array}{c}{\nabla }_{\overline{\mathbf{X}}}·{\sigma }_{\mathrm{Eff}}=0,\hfill \\ {\sigma }_{\mathrm{Eff}}={\langle {\mathbb{C}}_{\mathrm{I}}{\mathbb{M}}_{\mathrm{I}}+{\mathbb{C}}_{\mathrm{I}}+{\mathbb{C}}_{\mathrm{II}}{\mathbb{M}}_{\mathrm{II}}+{\mathbb{C}}_{\mathrm{II}}\rangle }_{\mathrm{s}}{\nabla }_{\overline{\mathbf{X}}}{\mathbf{u}}^{\left(0\right)}.\hfill \end{array}\hfill \end{array}\right.$$

(254)

We can make the identification that ${\sigma }_{\mathrm{Eff}}=\widehat{\mathcal{C}}:{\mathrm{Grad}}_{\overline{X}}{\mathbf{u}}^{\left(0\right)}$, then we can see that (254) matches the model obtained in [33] in the absence of plastic distortions. This is given by:

$$\left\{\begin{array}{c}\begin{array}{c}{\mathrm{Div}}_{\overline{X}}{\langle \widehat{\mathcal{C}}:{\mathrm{Grad}}_{\overline{X}}{\mathbf{u}}^{\left(0\right)}\rangle }_{\mathrm{s}}=0,\hfill \\ {\langle \widehat{\mathcal{C}}:{\mathrm{Grad}}_{\overline{X}}{\mathbf{u}}^{\left(0\right)}\rangle }_{\mathrm{s}}={\langle \mathcal{C}:\mathit{T}{\mathrm{Grad}}_{\overline{Y}}\xi +\mathcal{C}\rangle }_s{\mathrm{Grad}}_{\overline{X}}{\mathbf{u}}^{\left(0\right)}.\hfill \end{array}\hfill \end{array}\right.$$

(255)

Therefore, we can conclude that our model for nonlinear poroelastic composites can reduce to the model of [33] under the assumption of no plastic distortions.

5. Concluding Remarks

We derived a novel framework consisting of partial differential equations that describe the effective mechanical behaviour of nonlinear poroelastic composites. These structures are comprised of a porous hyperelastic matrix with embedded hyperelastic subphases, both of which interact with the fluid flowing in the pores. This type of structure is applicable to many real-world situations, including modelling of soft biological tissues. We began by considering the multiphase fluid–structure interaction (FSI) problem among all the constituents. The problem is closed by appropriate interface conditions arising from the continuity of stresses, displacements and velocities. We also performed a coordinate transformation on certain quantities in the FSI problem in order to obtain a formulation in the reference configuration. We exploited the length scale separation between the pore-scale (where the pores and elastic subphases are clearly visible) and the macroscale (average size of the material domain) to apply the asymptotic homogenization technique to the nondimensionalised system of PDEs in order to obtain a macroscale system of governing equations. We were able to recover previously known results in the literature by considering particular limit cases of our model. The new model encodes the detailed properties of the microstructure in its coefficients, that is the microstructural details are encoded in the effective hydraulic conductivity tensor, the Biot modulus and Biot’s tensor of coefficients. These are computed by solving the arising differential problems on the periodic cell.

The model obtained here is a generalisation of the formulations for multiphase elastoplastic composites [33] in the limit of no plastic distortions and the formulations of hyperelastic porous media [31]. The model is also a natural extension to the formulation for linear poroelastic composites [27]. All three of these models are recovered as particular cases of our new macroscale model. The key novelty of this work is the ability to describe a scenario where the hyperelastic matrix is inhomogeneous at the pore-scale, that is we are able to account for the interactions between various hyperelastic phases and the fluid flowing in the pores. This is generally the case in biological tissues. This means that this model is applicable to a wide range of biological scenarios including modelling lungs. The lungs have previously been approached in a biphasic (tissue and air) manner [9]. However, the lung microstructure is more complex, and there exist collagen and elastin fibres embedded in the matrix and in the fluid, so it could therefore be beneficial to use a nonlinear poroelastic composite approach to modelling. Another example is [12], where the interaction between pulsatile blood flow and the arterial wall mechanics was modelled. The blood flow was modelled as an incompressible viscous fluid, confined by Biot’s equations of poroelasticity for the artery wall. This model could be generalised by considering the wall as a nonlinear poroelastic composite of the type we modelled in this work. The formulation of standard nonlinear poroelasticity is applicable when the solid phase can be approximated as a homogeneous matrix. The linear formulation of poroelastic composites is applicable to situations where the deformations are small such as in hard hierarchical material such as bones (see, e.g., [28,47]). Our novel model provides a formulation that bridges a gap that has not previously been considered and can be grouped with those found in [17] for strength homogenization, nonsaturated microporomechanics, microporoplasticity and microporofracture and microporodamage theory as an extension to the nonlinear homogenization of porous media.

There are some limitations of the current model, and there are a number possible theoretical extensions that could potentially improve its applicability to certain biological systems. At present, the macroscale model is derived by accounting for a quasistatic regime and only considering the incompressibility of the fluid. Within this work, we used the de Saint-Venant strain energy function for the sake of simplicity. To study a wider range of scenarios, we would need to select a variety of more detailed constitutive laws and then use those to formulate the problem before finding the corresponding cell problems and the homogenized macroscale model. Strain–energy functions specific to certain applications could be used in our formulation, for example the Holzapfel–Ogden Law for the myocardium. Another possible extension, which would be of particular interest, is to incorporate growth and remodelling in our framework. Growth and remodelling are of particular importance to settings such as arteries or heart subject to disease or ageing. Finally, a third theoretical extension we could consider would be the assumption that the solid matrix and the subphases are both incompressible, in addition to the fluid incompressibility, which was already assumed in this work. This would require the incompressibility constraint to be imposed when defining the strain energy function and when determining the Piola stresses in the material. This would lead to alternative cell problems and macroscale model. Moreover, it would be possible to assume that the fluid was in fact compressible, and this would lead to the appearance of the fluid bulk modulus in the resulting Biot’s modulus of our system. This modification could be particularly useful to modelling applications in lungs where acoustic properties can be used to aid the diagnosis of respiratory diseases [9,48]. It would also be possible to consider a three-scale approach where there would exist an intermediate local scale between the pore-scale and macroscale that is still well separated. To do this, we could follow the approach taken in [49,50] for fibre-reinforced composites.

There are a number of potential next steps for this work; however, potentially the most important of these is to investigate the model by numerical simulations. Recently, the numerical simulations of the cell problems arising from the asymptotic homogenization technique when studying linear elastic composites and linear poroelasticity were carried out by [40,44]. The simulations for a linear poroelastic composite could be obtained by using the techniques in [40,44] to compute the poroelastic coefficients in the linear problem. Obtaining numerical results for the nonlinear macroscale model is significantly more complex due to the coupling between the macroscale variables and the cell problems. Developing suitable computational schemes is a current active area of research. A recent example of a proposed method that could be potentially used to solve the types of problems arising in this work is found in [45]. It is important to note that the potential results of any simulations should be validated by experimental data, which could be related to biological tissues.