Effective Governing Equations for Viscoelastic Composites

Abstract

We derive the governing equations for the overall behaviour of linear viscoelastic composites comprising two families of elastic inclusions, subphases and/or fibres, and an incompressible Newtonian fluid interacting with the solid phases at the microscale. We assume that the distance between each of the subphases is very small in comparison to the length of the whole material (the macroscale). We can exploit this sharp scale separation and apply the asymptotic (periodic) homogenization method (AHM) which decouples spatial scales and leads to the derivation of the new homogenised model. It does this via upscaling the fluid–structure interaction problem that arises between the multiple elastic phases and the fluid. As we do not assume that the fluid flow is characterised by a parabolic profile, the new macroscale model, which consists of partial differential equations, is of Kelvin–Voigt viscoelastic type (rather than poroelastic). The novel model has coefficients that encode the properties of the microstructure and are to be computed by solving a single local differential fluid–structure interaction (FSI) problem where the solid and the fluid phases are all present and described by the one problem. The model reduces to the case described by Burridge and Keller (1981) when there is only one elastic phase in contact with the fluid. This model is applicable when the distance between adjacent phases is smaller than the average radius of the fluid flowing in the pores, which can be the case for various highly heterogeneous systems encountered in real-world (e.g., biological, or geological) scenarios of interest.

1. Introduction

Materials can be described and classified in a number of ways. For example, viscous materials have a time-dependent deformation when subjected to a strain and also resist shear flow. Another example is elastic materials that will strain when a load is applied but will returned to the same unstressed state when the load is removed. Viscoelastic materials are characterised by having both an elastic and a viscous response under deformation, and therefore can be described as exhibiting time-dependent strain [1,2].

There are a large variety of physical settings where we have materials that exhibit a viscoelastic response. For example, in the human spine, under normal body weight, the disks get shorter with time, which means people are shorter in the evening, and lying down (removing the body weight) allows the disks to return to normal length by morning [3,4]. Human skin can also be described as viscoelastic and can be useful in diagnostic techniques and scar modelling [5,6,7]. The theory of viscoelasticity has also been used to consider materials that have a composite-like structure. Some examples of this are found in several biological contexts such as cortical and trabecular bone [8,9]. Viscoelastic composites including fibres have been very useful in engineering and manufacturing processes due to the fact that the properties can be optimised [10].

These viscoelastic systems are usually characterised by multiple physical scales of interest. They possess a fine scale structure and this is where the various solid and fluid interactions are clearly visible (the microscale). This scale is considerably smaller than the length associated with the complete viscoelastic material (the macroscale).

Viscoelastic composites have recently been addressed by a multiscale modelling approach. In [11,12,13], a variety of different methods have been used to incorporate the different microstructural information in an computationally feasible manner. The effective response of a material that is based upon the properties of the individual constituents can be described by micromechanical models. These properties can be the viscoelastic moduli or geometrical arrangement and volume fraction of the different microstructural constituents.

When addressing a multiscale system, there are a variety of approaches that can be taken to transform a fluid–structure interaction (FSI) problem into a complete macroscale governing system. These types of procedures can be categorised as homogenization techniques. Examples of such techniques are mixture theory, effective medium theory, volume averaging and asymptotic homogenization. The choice as to which method will be suitable for the system under study should be made depending on the application of the model and the information you wish to be encoded or available from the macroscale model; however, each approach has its own benefits.

The techniques including mixture theory and effective medium theory are micromechanical approaches. These are useful when it is desirable to obtain estimates of the poroelastic coefficients for particular pore geometries. These specific geometries include spherical, ellipsoidal, or penny-shaped or in the case of diluted pores. This was considered in [14]. Volume averaging techniques can be used when the goal is to derive the functional form of the equations governing the macroscale. When using this approach, the macroscale model coefficients are in general not encoding the microstructure of the material and require further data such as experimental results to obtain them. For a further summary and comparison of homogenization techniques, the reader is directed to [15,16].

Previously, the asymptotic homogenization method has been a popular approach to study poroelastic materials such as in [17,18,19]. This theory was then extended to consider the case of an inviscid fluid [20]. The method has also been employed to study elastic composites [21,22,23,24] and electroactive materials [25,26,27,28]. The technique has also previously been exploited to address problems in viscoelasticity. For example, the technique allowed for analytical closed form expressions for the effective coefficients of fibrous viscoelastic composites to be found in [29,30,31]. In [32], the authors study the homogenized properties of linear viscoelastic composite materials in three dimensions by means of a semi-analytical approach combining the asymptotic homogenization method (AHM) with numerical computations performed by finite element simulations. The work by [33] addresses the calculation of the effective properties of non-aging linear viscoelastic composite materials and investigates the effective creep and relaxation behaviour for a variety of fibre and inclusion reinforced structures.

In this manuscript, we apply the asymptotic homogenization method to upscale the FSI problem between two families of solid obstacles and a Newtonian fluid phase. In this case, the medium could not be considered a porous material. Both families of elastic subphases are in contact with the fluid and each other. We assume that the length scale at which the different solid phases and the fluid are clearly resolved can be compared with the distance between each of the phases which we call the microscale. We make the assumption that this scale is much smaller than the macroscale. We then carry out an upscaling on a system of equations that describe the behaviour of each phase and is closed by the continuity of stresses and displacements on the solid–solid interfaces, and the continuity of stresses and velocities on the fluid–solid interfaces. The application of the asymptotic homogenization method leads to a set of viscoelastic-type equations, which are a generalization of linear elastic composites [22] and those found in [17] in which the authors discussed the emergence of a viscoelastic model through the homogenization process of a composite material comprising a single solid phase and the fluid.

The novel model derived in this manuscript has coefficients that encode the properties of the microstructure and are to be computed by solving a single local differential FSI problem where the solid and the fluid phases are all present and described by the one problem. The new model has a key difference from formulations of previous viscoelastic models. That is, our model encodes multiple different elastic phases interacting with the fluid and each other, in comparison to some previous viscoelastic formulations such as [17] where there is only one fluid and one solid phase, or indeed in the case where the interactions between a solid viscoelastic phase and a fluid are considered [34,35]. The addition of the extra interactions between multiple phases to the model is particularly beneficial to physical applications such as in the bones [8,9] where there are a variety of different microscale constituents, such as collagen and minerals, and in articular cartilage [36,37]. The coefficients capture the differences in elastic and mechanical properties between the phases, as well as the discrepancies in the elastic properties at different points in the microstructure. The model can also be compared with that of poroelastic composites [38] where we see the important difference that the scaling of the fluid viscosity has. The macroscale model coefficients capture the geometry, elastic properties and interactions of the microscale constituents and are determined by solving a local differential FSI problem where the solid and the fluid phases are all present and described by the one problem.

We organise this paper in the following manner. Section 2 sees the introduction of the fluid–structure interaction problem describing the interactions occurring between the two families of elastic subphases and the fluid phase. We continue the development of the work by performing a multiple scales analysis of the FSI problem in Section 3. This leads to the derivation of the macroscale governing equations for the homogenized effective mechanical behaviour of viscoelastic composites. In Section 4, we present the macroscale model and compare the results here with previously known models in the literature such as [17,38]. In Section 5, the work concludes with discussions surrounding the limitations of the current model and further directions for future development.

2. The Fluid–Structure Interaction Problem

Our problem begins with a set $\Omega \in {\mathbb{R}}^3$, where $\Omega$ comprises two families of N disjoint solid subphases ${\Omega }_{\mathrm{a}}$ and ${\Omega }_{\mathrm{b}}$ and an interconnected fluid domain ${\Omega }_{\mathrm{f}}$. We have that

$$\begin{array}{cc}\hfill & {\displaystyle {\Omega }_{\mathrm{a}}=\bigcup _{\delta =1}^N{\Omega }_{\delta },}\hfill \end{array}$$

(1a)

$$\begin{array}{cc}\hfill & {\displaystyle {\Omega }_{\mathrm{b}}=\bigcup _{\beta =1}^N{\Omega }_{\beta }.}\hfill \end{array}$$

(1b)

We can write that the domain is $\overline{\Omega }={\overline{\Omega }}_{\mathrm{a}}\cup {\overline{\Omega }}_{\mathrm{b}}\cup {\overline{\Omega }}_{\mathrm{f}}$. To illustrate this structure, we have provided Figure 1, which is a two-dimensional schematic of the domain $\Omega$.

Before describing the equations that govern each of the domains in our structure we first wish to clarify the notation that will be used throughout this manuscript.

Remark 1

(Notation). We use the following for a generic field. For a scalar we use ordinary lowercase letters, e.g., v, for a vector we use boldface, e.g., $\mathbf{v}$, and monospace font $\mathsf{V}$ is used for second rank tensor. We use uppercase calligraphic letters for third rank tensors, for instance $\mathcal{V}$, and finally blackboard bold font $\mathbb{V}$ is used for fourth rank tensors.

We first require a balance equation in each of the solid domains ${\Omega }_{\delta }$ and ${\Omega }_{\beta }$. We choose to neglect volume forces and inertia, so we can write for all $\delta =1,\dots ,N,\mathrm{and}\beta =1,\dots ,N$

$$\begin{array}{cc}\hfill & \nabla ·{\mathsf{T}}_{\delta }=0\mathrm{in}{\Omega }_{\delta },\hfill \end{array}$$

(2a)

$$\begin{array}{cc}\hfill & \nabla ·{\mathsf{T}}_{\beta }=0\mathrm{in}{\Omega }_{\beta }.\hfill \end{array}$$

(2b)

We use ${\mathsf{T}}_{\delta }$ and ${\mathsf{T}}_{\beta }$ as the solid stress tensors for each of the subphases ${\Omega }_{\delta }$ and ${\Omega }_{\beta }$, respectively. We make the choice that both families of subphases are general linear elastic solids. This means we can write the Cauchy stress tensors for ${\mathsf{T}}_{\delta }$ and ${\mathsf{T}}_{\beta }$ as

$$\begin{array}{cc}\hfill {\mathsf{T}}_{\delta }=& {\mathbb{C}}_{\delta }\nabla {\mathbf{u}}_{\delta },\hfill \end{array}$$

(3a)

$$\begin{array}{cc}\hfill {\mathsf{T}}_{\beta }=& {\mathbb{C}}_{\beta }\nabla {\mathbf{u}}_{\beta },\hfill \end{array}$$

(3b)

where we have that ${\mathbf{u}}_{\delta }$ and ${\mathbf{u}}_{\beta }$ denote the elastic displacement of each of the individual subphases from each family. The tensors ${\mathbb{C}}_{\delta }$ and ${\mathbb{C}}_{\beta }$ appearing in (3a)–(3b) are the fourth rank elasticity tensors in each subphase. These can be written in components as $C_{pqrs}^{\delta }$ and $C_{pqrs}^{\beta }$, for $p,q,r,s=1,2,3$. Each tensor ${\mathbb{C}}_{\delta }$ and ${\mathbb{C}}_{\beta }$ has right minor and major symmetries, these are

$$\begin{array}{cc}\hfill & C_{pqrs}^{\delta }=C_{pqsr}^{\delta };C_{pqrs}^{\beta }=C_{pqsr}^{\beta },\hfill \end{array}$$

(4a)

$$\begin{array}{cc}\hfill & C_{pqrs}^{\delta }=C_{rspq}^{\delta };C_{pqrs}^{\beta }=C_{rspq}^{\beta },\hfill \end{array}$$

(4b)

and the left minor symmetries can be obtain by combining (4a) and (4b). We can use the left minor symmetries to rewrite constitutive Equations (3a) and (3b) as

$$\begin{array}{cc}\hfill {\mathsf{T}}_{\delta }=& {\mathbb{C}}_{\delta }\xi \left({\mathbf{u}}_{\delta }\right),\hfill \end{array}$$

(5a)

$$\begin{array}{cc}\hfill {\mathsf{T}}_{\beta }=& {\mathbb{C}}_{\beta }\xi \left({\mathbf{u}}_{\beta }\right),\hfill \end{array}$$

(5b)

where we define

$$\xi (•):=\frac{\nabla (•)+{\left(\nabla (•)\right)}^T}{2}$$

(6)

as the symmetric part of the gradient operator.

We also require a fluid balance equation, we can write

$$\nabla ·{\mathsf{T}}_{\mathsf{f}}=0\mathrm{in}{\Omega }_{\mathrm{f}},$$

(7)

with the stress tensor for the fluid, ${\mathsf{T}}_{\mathsf{f}}$. Our fluid can be described as incompressible and Newtonian, so therefore the constitutive law is

$${\mathsf{T}}_{\mathsf{f}}=-p\mathsf{I}+2\mu \xi \left(\mathbf{v}\right),$$

(8)

where we represent the viscosity fluid by $\mu$ the fluid velocity by $\mathbf{v}$, and where p is the pressure. Since we have an incompressible fluid we require

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

(9)

Using the fluid constitutive law (8) in the fluid balance Equation (7) and together with the incompressibility constraint (9), this leads to the Stokes’ problem

$$\mu {\nabla }^2\mathbf{v}=\nabla p\mathrm{in}{\Omega }_{\mathrm{f}}.$$

(10)

To have a complete FSI problem, we close it with appropriate interface conditions. These conditions are placed to describe the interactions between the fluid and the solid phases. The interfaces can be defined as follows: between the fluid and the $\delta$-th inclusion/fibre/subphase we have ${\Gamma }_{\delta }:=\partial {\Omega }_{\delta }\cap \partial {\Omega }_{\mathrm{f}}$; and between the fluid phase and the $\beta$ inclusion/fibre/subphase as ${\Gamma }_{\beta }:=\partial {\Omega }_{\beta }\cap \partial {\Omega }_{\mathrm{f}}$. Across each ${\Gamma }_{\delta }$ and ${\Gamma }_{\beta }$ we then enforce continuity of velocities and stresses, i.e.,

$$\begin{array}{cc}\hfill \hspace{1em}\hspace{1em}{\dot{\mathbf{u}}}_{\delta }=\mathbf{v}& \mathrm{on}{\Gamma }_{\delta },\hfill \end{array}$$

(11a)

$$\begin{array}{cc}\hfill {\mathsf{T}}_{\mathsf{f}}{\mathbf{n}}_{\delta }={\mathsf{T}}_{\delta }{\mathbf{n}}_{\delta }& \mathrm{on}{\Gamma }_{\delta },\hfill \end{array}$$

(11b)

$$\begin{array}{cc}\hfill \hspace{1em}\hspace{1em}{\dot{\mathbf{u}}}_{\beta }=\mathbf{v}& \mathrm{on}{\Gamma }_{\beta },\hfill \end{array}$$

(11c)

$$\begin{array}{cc}\hfill {\mathsf{T}}_{\mathsf{f}}{\mathbf{n}}_{\beta }={\mathsf{T}}_{\beta }{\mathbf{n}}_{\beta }& \mathrm{on}{\Gamma }_{\beta },\hfill \end{array}$$

(11d)

for all $\delta =1,\dots ,N,\mathrm{and}\beta =1,\dots ,N$. We have used the notation ${\dot{\mathbf{u}}}_{\delta }$ and ${\dot{\mathbf{u}}}_{\beta }$ to describe the solid velocities in each inclusion/fibre/subphase ${\Omega }_{\delta }$ and ${\Omega }_{\beta }$, respectively. We must also note that ${\mathbf{n}}_{\delta }$ and ${\mathbf{n}}_{\beta }$ are the unit outward vectors that are normal to the interfaces ${\Gamma }_{\delta }$ and ${\Gamma }_{\beta }$. That is, these normals point into the fluid domain ${\Omega }_{\mathrm{f}}$.

Remark 2

(Frequency Domain). In this work, we embrace the approach of [17], which means that we consider time harmonic motion. As such, a time dependent field, say $\phi (\mathbf{x},t)$, can be decomposed in a solely spatially varying component and a harmonic time variation, i.e., we assume that $\phi (\mathbf{x},t)={\phi }_0\left(\mathbf{x}\right)exp\left(i\omega t\right)$. For instance, this then means that $\frac{\partial \phi }{\partial t}=i\omega t$. In the remainder of this work, we will identify each field with their spatially varying component and omit the subscript ${}_0$ for the sake of simplicity of notation. This approach can be carried out without loss of generality as in principle the problem could be formulated in the frequency domain and so, in other words, every sufficiently smooth time dependency could be taken into account by spanning over the frequency domain by means of the Fourier transform operator.

As such, assuming continuity of velocities at the interfaces, by means of Remark 2, we have

$$\begin{array}{cc}\hfill i\omega {\mathbf{u}}_{\delta }=\mathbf{v}& \mathrm{on}{\Gamma }_{\delta },\hfill \end{array}$$

(12a)

$$\begin{array}{cc}\hfill i\omega {\mathbf{u}}_{\beta }=\mathbf{v}& \mathrm{on}{\Gamma }_{\beta }.\hfill \end{array}$$

(12b)

We also have the boundary between each of the different solid phases, which we write as ${\Gamma }_{\delta \beta }:=\partial {\Omega }_{\delta }\cap \partial {\Omega }_{\mathrm{b}}$. We then impose continuity of stresses and displacements, which can be written as

$$\begin{array}{cc}\hfill {\mathsf{T}}_{\delta }{\mathbf{n}}_{\delta \beta }={\mathsf{T}}_{\beta }{\mathbf{n}}_{\delta \beta }& \mathrm{on}{\Gamma }_{\delta \beta }.\hfill \end{array}$$

(13a)

$$\begin{array}{cc}\hfill \hspace{1em}\hspace{1em}{\mathbf{u}}_{\delta }={\mathbf{u}}_{\beta }& \mathrm{on}{\Gamma }_{\delta \beta },\hfill \end{array}$$

(13b)

for all $\delta =1\dots N,\mathrm{and}\beta =1,\dots ,N$. The unit vector ${\mathbf{n}}_{\delta \beta }$ appearing in (13a) is normal to the interface ${\Gamma }_{\delta \beta }$ and is pointing into the inclusion ${\Omega }_{\delta }$.

We have now introduced all the equations necessary to carry out a multiple scales analysis. We do this by first non-dimensionalising the fluid–structure interaction problem we have just formed. This is done by introducing two distinct length scales. We use this sharp scale separation and apply the asymptotic homogenization method to the non-dimensionalised FSI problem. This leads to the derivation of the effective macroscale governing equations that describe the viscoelastic material.

3. Multiple Scales Analysis

We can now summarise the FSI problem we introduced in the previous section. We note that the fields in each of the equations here are all defined in the time domain and where the appearance of a time derivative is shown by the multiplication of the field by $i\omega$. We have

$$\begin{array}{cc}\hfill \nabla ·{\mathsf{T}}_{\delta }=0& \mathrm{in}{\Omega }_{\delta },\hfill \end{array}$$

(14a)

$$\begin{array}{cc}\hfill \nabla ·{\mathsf{T}}_{\beta }=0& \mathrm{in}{\Omega }_{\beta },\hfill \end{array}$$

(14b)

$$\begin{array}{cc}\hfill \nabla ·{\mathsf{T}}_{\mathsf{f}}=0& \mathrm{in}{\Omega }_{\mathrm{f}},\hfill \end{array}$$

(14c)

$$\begin{array}{cc}\hfill \nabla ·\mathbf{v}=0& \mathrm{in}{\Omega }_{\mathrm{f}},\hfill \end{array}$$

(14d)

$$\begin{array}{cc}\hfill \hspace{1em}i\omega {\mathbf{u}}_{\delta }=\mathbf{v}& \mathrm{on}{\Gamma }_{\delta },\hfill \end{array}$$

(14e)

$$\begin{array}{cc}\hfill \hspace{1em}i\omega {\mathbf{u}}_{\beta }=\mathbf{v}& \mathrm{on}{\Gamma }_{\beta },\hfill \end{array}$$

(14f)

$$\begin{array}{cc}\hfill {\mathsf{T}}_{\mathsf{f}}{\mathbf{n}}_{\delta }={\mathsf{T}}_{\delta }{\mathbf{n}}_{\delta }& \mathrm{on}{\Gamma }_{\delta },\hfill \end{array}$$

(14g)

$$\begin{array}{cc}\hfill {\mathsf{T}}_{\mathsf{f}}{\mathbf{n}}_{\beta }={\mathsf{T}}_{\beta }{\mathbf{n}}_{\beta }& \mathrm{on}{\Gamma }_{\beta },\hfill \end{array}$$

(14h)

$$\begin{array}{cc}\hfill {\mathsf{T}}_{\delta }{\mathbf{n}}_{\delta \beta }={\mathsf{T}}_{\beta }{\mathbf{n}}_{\delta \beta }& \mathrm{on}{\Gamma }_{\delta \beta },\hfill \end{array}$$

(14i)

$$\begin{array}{cc}\hfill \hspace{1em}\hspace{1em}{\mathbf{u}}_{\delta }={\mathbf{u}}_{\beta }& \mathrm{on}{\Gamma }_{\delta \beta }.\hfill \end{array}$$

(14j)

We can then use the constitutive relationships (5a), (5b), and (8), and the incompressibility constraint (14d), to rewrite the balance Equations (14a)–(14c) as

$$\begin{array}{cc}\hfill \nabla ·\left({\mathbb{C}}_{\delta }\xi \left({\mathbf{u}}_{\delta }\right)\right)=0& \mathrm{in}{\Omega }_{\delta }\hfill \end{array}$$

(15a)

$$\begin{array}{cc}\hfill \nabla ·\left({\mathbb{C}}_{\beta }\xi \left({\mathbf{u}}_{\beta }\right)\right)=0& \mathrm{in}{\Omega }_{\beta }\hfill \end{array}$$

(15b)

$$\begin{array}{cc}\hfill \hspace{1em}\hspace{1em}\hspace{1em}\mu {\nabla }^2\mathbf{v}=\nabla p& \mathrm{in}{\Omega }_{\mathrm{f}},\hfill \end{array}$$

(15c)

for all $\delta =1\dots N,\mathrm{and}\beta =1,\dots ,N$. To close the problem (14a)–(14j), we place appropriate external boundary conditions on $\partial \Omega$.

3.1. Non-Dimensionalisation of the FSI Equations

As the system we are considering is multiscale, we identify two typical length scales. We first associate L with the size of the whole material $\Omega$ (the macroscale). The second scale of interest we denote with d and this refers to the microscale, which we assume can be compared to the distance between each of the subphases and the fluid. In order to clearly see the distinction between the scales, we carry out a non-dimensionalisation of the system (14a)–(14j). To perform the non-dimensional analysis, we assume that the system is characterised by a reference pressure gradient C, as well as the length scales to obtain

$$\begin{array}{c}\mathbf{x}=L{\mathbf{x}}^{\prime },{\mathbb{C}}_{\delta }=CL{\mathbb{C}}_{\delta }^{\prime },{\mathbb{C}}_{\beta }=CL{\mathbb{C}}_{\beta }^{\prime },\\ {\mathbf{u}}_{\delta }=L{\mathbf{u}}_{\delta }^{\prime },{\mathbf{u}}_{\beta }=L{\mathbf{u}}_{\beta }^{\prime },\mathbf{v}=\frac{C{L}^2}{\mu }{\mathbf{v}}^{\prime },p=CL{p}^{\prime }.\end{array}$$

(16)

It is possible to choose different scalings for the fluid velocity and these account for the effective behaviour of fluid flow in porous media [17,19,38,39].

We are able to use (16), and noticing that we can write

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

(17)

we obtain the non-dimensionalised fluid–structure interaction problem (14a)–(14j),

$$\begin{array}{cc}\hfill \nabla ·{\mathsf{T}}_{\delta }=0& \mathrm{in}{\Omega }_{\delta }\hfill \end{array}$$

(18a)

$$\begin{array}{cc}\hfill \nabla ·{\mathsf{T}}_{\beta }=0& \mathrm{in}{\Omega }_{\beta }\hfill \end{array}$$

(18b)

$$\begin{array}{cc}\hfill \nabla ·{\mathsf{T}}_{\mathsf{f}}=0& \mathrm{in}{\Omega }_{\mathrm{f}}\hfill \end{array}$$

(18c)

$$\begin{array}{cc}\hfill \nabla ·\mathbf{v}=0& \mathrm{in}{\Omega }_{\mathrm{f}}\hfill \end{array}$$

(18d)

$$\begin{array}{cc}\hfill \hspace{1em}i\omega {\mathbf{u}}_{\delta }=\mathbf{v}& \mathrm{on}{\Gamma }_{\delta }\hfill \end{array}$$

(18e)

$$\begin{array}{cc}\hfill \hspace{1em}i\omega {\mathbf{u}}_{\beta }=\mathbf{v}& \mathrm{on}{\Gamma }_{\beta }\hfill \end{array}$$

(18f)

$$\begin{array}{cc}\hfill {\mathsf{T}}_{\mathsf{f}}{\mathbf{n}}_{\delta }={\mathsf{T}}_{\delta }{\mathbf{n}}_{\delta }& \mathrm{on}{\Gamma }_{\delta }\hfill \end{array}$$

(18g)

$$\begin{array}{cc}\hfill {\mathsf{T}}_{\mathsf{f}}{\mathbf{n}}_{\beta }={\mathsf{T}}_{\beta }{\mathbf{n}}_{\beta }& \mathrm{on}{\Gamma }_{\mathrm{b}}\hfill \end{array}$$

(18h)

$$\begin{array}{cc}\hfill {\mathsf{T}}_{\delta }{\mathbf{n}}_{\delta \beta }={\mathsf{T}}_{\beta }{\mathbf{n}}_{\delta \beta }& \mathrm{on}{\Gamma }_{\delta \beta } \hfill \end{array}$$

(18i)

$$\begin{array}{cc}\hfill \hspace{1em}\hspace{1em}{\mathbf{u}}_{\delta }={\mathbf{u}}_{\beta }& \mathrm{on}{\Gamma }_{\delta \beta }\hfill \end{array}$$

(18j)

For all $\delta =1,\dots ,N,\mathrm{and}\beta =1,\dots ,N$, the primes have been removed to simplify the notation. The non-dimensionalised constitutive laws (5a), (5b), and (8) are

$$\begin{array}{cc}\hfill {\mathsf{T}}_{\mathsf{f}}=-p\mathbf{I}+\xi \left(\mathbf{v}\right)\hspace{1em}& \hfill \end{array}$$

(19a)

$$\begin{array}{cc}\hfill {\mathsf{T}}_{\delta }={\mathbb{C}}_{\delta }\xi \left({\mathbf{u}}_{\delta }\right)& \hfill \end{array}$$

(19b)

$$\begin{array}{cc}\hfill {\mathsf{T}}_{\beta }={\mathbb{C}}_{\beta }\xi \left({\mathbf{u}}_{\beta }\right),& \hfill \end{array}$$

(19c)

and using these in the balance Equations (15a)–(15c) gives

$$\begin{array}{cc}\hfill \hspace{1em}\hspace{1em}\hspace{1em}{\nabla }^2\mathbf{v}=\nabla p& \mathrm{in}{\Omega }_{\mathrm{f}}\hfill \end{array}$$

(20a)

$$\begin{array}{cc}\hfill \nabla ·\left({\mathbb{C}}_{\delta }\xi \left({\mathbf{u}}_{\delta }\right)\right)=0& \mathrm{in}{\Omega }_{\delta }\hfill \end{array}$$

(20b)

$$\begin{array}{cc}\hfill \nabla ·\left({\mathbb{C}}_{\beta }\xi \left({\mathbf{u}}_{\beta }\right)\right)=0& \mathrm{in}{\Omega }_{\mathrm{b}}.\hfill \end{array}$$

(20c)

We are now ready to introduce the asymptotic homogenization method that will upscale the non-dimensional FSI problem (18a)–(20c) to macroscale governing equations by making the assumption that the two scales (micro and macro) are well separated.

3.2. The Asymptotic Homogenization Method

We now introduce the rules and assumptions associated with the asymptotic homogenization method which are then used to obtain the macroscale model from the FSI problem (18a)–(20c). We begin by making the assumption that the microscale length, which we denote by d, is much smaller than the average size of the viscoelastic material which has length L. That is,

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

(21)

We require a microscale spatial variable which will describe how each field varies on the microscale, that is

$$\mathbf{y}=\frac{\mathbf{x}}{ϵ}.$$

(22)

We have two spatial variables $\mathbf{x}$ and $\mathbf{y}$ and it is assumed that these are formally independent, with $\mathbf{x}$ representing the macroscale and $\mathbf{y}$ the microscale. We have that the gradient operator will also transform

$$\nabla \to {\nabla }_{\mathbf{x}}+\frac{1}{ϵ}{\nabla }_{\mathbf{y}}.$$

(23)

We require that all the fields in (18a)–(20c) are functions of both spatial variables and have a power series in $ϵ$ representation, i.e.,

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

(24)

where $\phi$ is a generic field occurring in our analysis.

Remark 3

(Microscale periodicity). To simplify the analysis in this work, we restrict our attention to a single subset of the domain which we call the periodic cell. This periodic cell may have a variety of different subphases, each of which can have a different geometry and elastic properties, this is depicted in Figure 2. For this to be possible, we make the assumption that every field ${\phi }^{\left(l\right)}$ arising in (18a)–(20c) is y-periodic. By making this assumption, we can solve the microscale differential problems arising from using the asymptotic homogenization method on just a finite subset of the material. This assumption need not be made and the analysis can be carried using a less strict assumption, i.e., the local boundedness of fields. This approach, however, only allows us to determine the functional form of the macroscale equations and the model coefficients are to be obtained by solving microscale problems on the whole microstructure of the material. This makes solving the model very computationally expensive when using the local boundedness of fields approach in comparison to microscale periodicity. Some examples of this are found in [17,40,41].

Remark 4

(Macroscopic uniformity). It is clear that the microscale geometry can differ depending on the macroscale point considered. This dependence is in general neglected in most works. We will make the assumption that the viscoelastic material is macroscopically uniform. That is, we assume the microscale geometry does not change with macroscale variable $\mathbf{x}$. By making this assumption, we have the following

$${\int }_{\Omega }{\nabla }_{\mathbf{x}}·\left(•\right)d\mathbf{y}={\nabla }_{\mathbf{x}}·{\int }_{\Omega }\left(•\right)d\mathbf{y}.$$

(25)

In the case that $\Omega =\Omega \left(\mathbf{x}\right)$, then Equation (25) does not hold. This leads to an application of the generalized Reynolds’ transport theorem which can give rise to additional contributions on the macroscale; see, e.g., [39,40,42].

Remark 5

(Local Geometry). In the set up of the problem so far, we have considered a fluid with many embedded subphases that are interacting. For the sake of clarity, we assume that each family provides only one subphase for each periodic cell, as shown in Figure 3. This assumption does not affect the generality of the properties of the model. If a particular application requires a variety of subphases in the periodic cell, then this can easily be extended, as in [21]. Therefore, the indexes δ and β are not necessary and we can amend the notation. We identify the domain Ω with the corresponding periodic cell which comprises two solid subphases (one from each family) and fluid, with each of the phases denoted by ${\Omega }_{\mathrm{a}}$, ${\Omega }_{\mathrm{b}}$, and ${\Omega }_{\mathrm{f}}$, respectively. We can also simplify the notation used for the interfaces between the different phases, so we have ${\Gamma }_{\mathrm{a}}:=\partial {\Omega }_{\mathrm{a}}\cap \partial {\Omega }_{\mathrm{f}}$ between the fluid and the inclusion $\mathrm{a}$, ${\Gamma }_{\mathrm{b}}:=\partial {\Omega }_{\mathrm{b}}\cap \partial {\Omega }_{\mathrm{f}}$ is the interface between the fluid and the inclusion $\mathrm{b}$, and ${\Gamma }_{\mathrm{ab}}:=\partial {\Omega }_{\mathrm{a}}\cap \partial {\Omega }_{\mathrm{b}}$, is the interface between the two solid phases. These interfaces have corresponding unit normal vectors ${\mathbf{n}}_{\mathrm{a}}$, ${\mathbf{n}}_{\mathrm{b}}$, and ${\mathbf{n}}_{\mathrm{ab}}$.

3.3. Deriving the Macroscale Model

We now apply the assumptions of the asymptotic homogenization method, i.e., (23) and (24), to Equations (18a)–(18j), accounting also for periodicity, to obtain

$$\begin{array}{cc}\hfill {\nabla }_{\mathbf{y}}·{\mathsf{T}}_{\mathrm{a}}^{ϵ}+ϵ{\nabla }_{\mathbf{x}}·{\mathsf{T}}_{\mathrm{a}}^{ϵ}=0& \mathrm{in}{\Omega }_{\mathrm{a}}\hfill \end{array}$$

(26a)

$$\begin{array}{cc}\hfill {\nabla }_{\mathbf{y}}·{\mathsf{T}}_{\mathrm{b}}^{ϵ}+ϵ{\nabla }_{\mathbf{x}}·{\mathsf{T}}_{\mathrm{b}}^{ϵ}=0& \mathrm{in}{\Omega }_{\mathrm{b}}\hfill \end{array}$$

(26b)

$$\begin{array}{cc}\hfill {\nabla }_{\mathbf{y}}·{\mathsf{T}}_{\mathsf{f}}^{ϵ}+ϵ{\nabla }_{\mathbf{x}}·{\mathsf{T}}_{\mathsf{f}}^{ϵ}=0& \mathrm{in}{\Omega }_{\mathrm{f}}\hfill \end{array}$$

(26c)

$$\begin{array}{cc}\hfill {\nabla }_{\mathbf{y}}·{\mathbf{v}}^{ϵ}+ϵ{\nabla }_{\mathbf{x}}·{\mathbf{v}}^{ϵ}=0& \mathrm{in}{\Omega }_{\mathrm{f}}\hfill \end{array}$$

(26d)

$$\begin{array}{cc}\hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}i\omega {\mathbf{u}}_{\mathrm{a}}^{ϵ}={\mathbf{v}}^{ϵ}& \mathrm{on}{\Gamma }_{\mathrm{a}}\hfill \end{array}$$

(26e)

$$\begin{array}{cc}\hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}i\omega {\mathbf{u}}_{\mathrm{b}}^{ϵ}={\mathbf{v}}^{ϵ}& \mathrm{on}{\Gamma }_{\mathrm{b}}\hfill \end{array}$$

(26f)

$$\begin{array}{cc}\hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{\mathsf{T}}_{\mathsf{f}}^{ϵ}{\mathbf{n}}_{\mathrm{a}}={\mathsf{T}}_{\mathrm{a}}^{ϵ}{\mathbf{n}}_{\mathrm{a}}& \mathrm{on}{\Gamma }_{\mathrm{a}}\hfill \end{array}$$

(26g)

$$\begin{array}{cc}\hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{\mathsf{T}}_{\mathsf{f}}^{ϵ}{\mathbf{n}}_{\mathrm{b}}={\mathsf{T}}_{\mathrm{b}}^{ϵ}{\mathbf{n}}_{\mathrm{b}}& \mathrm{on}{\Gamma }_{\mathrm{b}}\hfill \end{array}$$

(26h)

$$\begin{array}{cc}\hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{\mathsf{T}}_{\mathrm{a}}^{ϵ}{\mathbf{n}}_{\mathrm{ab}}={\mathsf{T}}_{\mathrm{b}}^{ϵ}{\mathbf{n}}_{\mathrm{ab}}& \mathrm{on}{\Gamma }_{\mathrm{ab}}\hfill \end{array}$$

(26i)

$$\begin{array}{cc}\hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{\mathbf{u}}_{\mathrm{a}}^{ϵ}={\mathbf{u}}_{\mathrm{b}}^{ϵ}& \mathrm{on}{\Gamma }_{\mathrm{ab}}.\hfill \end{array}$$

(26j)

The constitutive equations for the fluid and solid stress tensors under these assumptions are

$$\begin{array}{cc}\hfill \hspace{1em}ϵ{\mathsf{T}}_{\mathsf{f}}^{ϵ}=-ϵp^{ϵ}\mathbf{I}+{\xi }_{\mathbf{y}}\left({\mathbf{v}}^{ϵ}\right)+ϵ{\xi }_{\mathbf{x}}\left({\mathbf{v}}^{ϵ}\right)& \end{array}$$

(27a)

$$\begin{array}{cc}\hfill \hspace{1em}\hspace{1em}ϵ{\mathsf{T}}_{\mathrm{a}}^{ϵ}={\mathbb{C}}_{\mathrm{a}}{\xi }_{\mathbf{y}}\left({\mathbf{u}}_{\mathrm{a}}^{ϵ}\right)+ϵ{\mathbb{C}}_{\mathrm{a}}{\xi }_{\mathbf{x}}\left({\mathbf{u}}_{\mathrm{a}}^{ϵ}\right)& \end{array}$$

(27b)

$$\begin{array}{c}\hfill \hspace{1em}\hspace{1em}ϵ{\mathsf{T}}_{\mathrm{b}}^{ϵ}={\mathbb{C}}_{\mathrm{b}}{\xi }_{\mathbf{y}}\left({\mathbf{u}}_{\mathrm{b}}^{ϵ}\right)+ϵ{\mathbb{C}}_{\mathrm{b}}{\xi }_{\mathbf{x}}\left({\mathbf{u}}_{\mathrm{b}}^{ϵ}\right),\end{array}$$

(27c)

and the balance equations are

$$\begin{array}{cc}\hfill & {\nabla }_{\mathbf{y}}·\left({\mathbb{C}}_{\mathrm{a}}{\xi }_{\mathbf{y}}\left({\mathbf{u}}_{\mathrm{a}}^{ϵ}\right)\right)+ϵ{\nabla }_{\mathbf{y}}·\left({\mathbb{C}}_{\mathrm{a}}{\xi }_{\mathbf{x}}\left({\mathbf{u}}_{\mathrm{a}}^{ϵ}\right)\right)+ϵ{\nabla }_{\mathbf{x}}·\left({\mathbb{C}}_{\mathrm{a}}{\xi }_{\mathbf{y}}\left({\mathbf{u}}_{\mathrm{a}}^{ϵ}\right)\right)+{ϵ}^2{\nabla }_{\mathbf{x}}·\left({\mathbb{C}}_{\mathrm{a}}{\xi }_{\mathbf{x}}\left({\mathbf{u}}_{\mathrm{a}}^{ϵ}\right)\right)\hfill \\ \hfill & =0\mathrm{in}{\Omega }_{\mathrm{a}}\hfill \end{array}$$

(28a)

$$\begin{array}{cc}\hfill & {\nabla }_{\mathbf{y}}·\left({\mathbb{C}}_{\mathrm{b}}{\xi }_{\mathbf{y}}\left({\mathbf{u}}_{\mathrm{b}}^{ϵ}\right)\right)+ϵ{\nabla }_{\mathbf{y}}·\left({\mathbb{C}}_{\mathrm{b}}{\xi }_{\mathbf{x}}\left({\mathbf{u}}_{\mathrm{b}}^{ϵ}\right)\right)+ϵ{\nabla }_{\mathbf{x}}·\left({\mathbb{C}}_{\mathrm{b}}{\xi }_{\mathbf{y}}\left({\mathbf{u}}_{\mathrm{b}}^{ϵ}\right)\right)+{ϵ}^2{\nabla }_{\mathbf{x}}·\left({\mathbb{C}}_{\mathrm{b}}{\xi }_{\mathbf{x}}\left({\mathbf{u}}_{\mathrm{b}}^{ϵ}\right)\right)\hfill \\ \hfill & =0\mathrm{in}{\Omega }_{\mathrm{b}}\hfill \end{array}$$

(28b)

$$\begin{array}{cc}\hfill & {ϵ}^2{\nabla }_{\mathbf{x}}^2{\mathbf{v}}^{ϵ}+ϵ{\nabla }_{\mathbf{x}}·\left({\nabla }_{\mathbf{y}}{\mathbf{v}}^{ϵ}\right)+ϵ{\nabla }_{\mathbf{y}}·\left({\nabla }_{\mathbf{x}}{\mathbf{v}}^{ϵ}\right)+{\nabla }_{\mathbf{y}}^2{\mathbf{v}}^{ϵ}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hfill \\ \hfill & =ϵ{\nabla }_{\mathbf{y}}{p}^{ϵ}+{ϵ}^2{\nabla }_{\mathbf{x}}{p}^{ϵ}\mathrm{in}{\Omega }_{\mathrm{f}}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hfill \end{array}$$

(28c)

Since all the fields can be written in terms of a power series of the type (24) in (26a)–(28c), we can equate the coefficients of ${ϵ}^l$ for $l=0,1,\dots$. These obtained equations are then used to derive the macroscale model for the material. The model will be expressed in terms of relevant leading order fields. For any terms in the model that retain a dependence on the microscale variable $\mathbf{y}$, we apply the integral average. This average is defined as

$${\langle \phi \rangle }_i:=\frac{1}{|\Omega |}{\int }_{{\Omega }_i}\phi (\mathbf{x},\mathbf{y},t)d\mathbf{y}i=f,\mathrm{a},\mathrm{b}$$

(29)

where $\phi$ is a generic field and $|\Omega |$ is the volume of the domain and we note that $|\Omega | = |{\Omega }_{\mathrm{f}}| + |{\Omega }_{\mathrm{a}}| + |{\Omega }_{\mathrm{b}}|$. Due to the assumption of $\mathbf{y}$-periodicity, the integral average can be performed over one representative cell. Therefore, we can say that (29) represents a cell average.

We begin by equating coefficients of ${ϵ}^0$ in (26a)–(26j) to obtain

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

(30a)

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

(30b)

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

(30c)

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

(30d)

$$\begin{array}{cc}\hfill i\omega {\mathbf{u}}_{\mathrm{a}}^{\left(0\right)}={\mathbf{v}}^{\left(0\right)}& \mathrm{on}{\Gamma }_{\mathrm{a}},\hfill \end{array}$$

(30e)

$$\begin{array}{cc}\hfill i\omega {\mathbf{u}}_{\mathrm{b}}^{\left(0\right)}={\mathbf{v}}^{\left(0\right)}& \mathrm{on}{\Gamma }_{\mathrm{b}},\hfill \end{array}$$

(30f)

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

(30g)

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

(30h)

$$\begin{array}{cc}\hfill {\mathsf{T}}_{\mathrm{a}}^{\left(0\right)}{\mathbf{n}}_{\mathrm{ab}}={\mathsf{T}}_{\mathrm{b}}^{\left(0\right)}{\mathbf{n}}_{\mathrm{ab}}& \mathrm{on}{\Gamma }_{\mathrm{ab}},\hfill \end{array}$$

(30i)

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

(30j)

The ${ϵ}^0$ coefficients of the constitutive Equations (27a)–(27c) read

$$\begin{array}{cc}\hfill {\xi }_{\mathbf{y}}\left({\mathbf{v}}^{\left(0\right)}\right)=0& \mathrm{in}{\Omega }_{\mathrm{f}},\hfill \end{array}$$

(31a)

$$\begin{array}{cc}\hfill {\mathbb{C}}_{\mathrm{a}}{\xi }_{\mathbf{y}}\left({\mathbf{u}}_{\mathrm{a}}^{\left(0\right)}\right)=0& \mathrm{in}{\Omega }_{\mathrm{a}},\hfill \end{array}$$

(31b)

$$\begin{array}{cc}\hfill {\mathbb{C}}_{\mathrm{b}}{\xi }_{\mathbf{y}}\left({\mathbf{u}}_{\mathrm{b}}^{\left(0\right)}\right)=0& \mathrm{in}{\Omega }_{\mathrm{b}},\hfill \end{array}$$

(31c)

and the balance Equations (28a)–(28c) have coefficients of ${ϵ}^0$

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

(32a)

$$\begin{array}{cc}\hfill {\nabla }_{\mathbf{y}}·\left({\mathbb{C}}_{\mathrm{b}}{\nabla }_{\mathbf{y}}\left({\mathbf{u}}_{\mathrm{b}}^{\left(0\right)}\right)\right)=0& \mathrm{in}{\Omega }_{\mathrm{b}},\hfill \end{array}$$

(32b)

$$\begin{array}{cc}\hfill {\nabla }_{\mathbf{y}}^2{\mathbf{v}}^{\left(0\right)}=0& \mathrm{in}{\Omega }_{\mathrm{f}}.\hfill \end{array}$$

(32c)

We now similarly consider the coefficients of ${ϵ}^1$ in (26a)–(26j) which leads to

$$\begin{array}{cc}\hfill {\nabla }_{\mathbf{y}}·{\mathsf{T}}_{\mathrm{a}}^{\left(1\right)}+{\nabla }_{\mathbf{x}}·{\mathsf{T}}_{\mathrm{a}}^{\left(0\right)}=0& \mathrm{in}{\Omega }_{\mathrm{a}},\hfill \end{array}$$

(33a)

$$\begin{array}{cc}\hfill {\nabla }_{\mathbf{y}}·{\mathsf{T}}_{\mathrm{b}}^{\left(1\right)}+{\nabla }_{\mathbf{x}}·{\mathsf{T}}_{\mathrm{b}}^{\left(0\right)}=0& \mathrm{in}{\Omega }_{\mathrm{b}},\hfill \end{array}$$

(33b)

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

(33c)

$$\begin{array}{cc}\hfill {\nabla }_{\mathbf{y}}·{\mathbf{v}}^{\left(1\right)}+{\nabla }_{\mathbf{x}}·{\mathbf{v}}^{\left(0\right)}=0& \mathrm{in}{\Omega }_{\mathrm{f}},\hfill \end{array}$$

(33d)

$$\begin{array}{cc}\hfill i\omega {\mathbf{u}}_{\mathrm{a}}^{\left(1\right)}={\mathbf{v}}^{\left(1\right)}& \mathrm{on}{\Gamma }_{\mathrm{a}},\hfill \end{array}$$

(33e)

$$\begin{array}{cc}\hfill i\omega {\mathbf{u}}_{\mathrm{b}}^{\left(1\right)}={\mathbf{v}}^{\left(1\right)}& \mathrm{on}{\Gamma }_{\mathrm{b}},\hfill \end{array}$$

(33f)

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

(33g)

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

(33h)

$$\begin{array}{cc}\hfill {\mathsf{T}}_{\mathrm{a}}^{\left(1\right)}{\mathbf{n}}_{\mathrm{ab}}={\mathsf{T}}_{\mathrm{b}}^{\left(1\right)}{\mathbf{n}}_{\mathrm{ab}}& \mathrm{on}{\Gamma }_{\mathrm{ab}},\hfill \end{array}$$

(33i)

$$\begin{array}{cc}\hfill {\mathbf{u}}_{\mathrm{a}}^{\left(1\right)}={\mathbf{u}}_{\mathrm{b}}^{\left(1\right)}& \mathrm{on}{\Gamma }_{\mathrm{ab}}.\hfill \end{array}$$

(33j)

Then, the constitutive Equations (27a)–(27c) have coefficients of ${ϵ}^1$ given by

$$\begin{array}{cc}\hfill {{\mathsf{T}}_{\mathsf{f}}}^{\left(0\right)}=-p^{\left(0\right)}\mathbf{I}+{\xi }_{\mathbf{y}}\left({\mathbf{v}}^{\left(1\right)}\right)+{\xi }_{\mathbf{x}}\left({\mathbf{v}}^{\left(0\right)}\right)& \mathrm{in}{\Omega }_{\mathrm{f}},\hfill \end{array}$$

(34a)

$$\begin{array}{cc}\hfill {\mathsf{T}}_{\mathrm{a}}^{\left(0\right)}={\mathbb{C}}_{\mathrm{a}}{\xi }_{\mathbf{y}}\left({\mathbf{u}}_{\mathrm{a}}^{\left(1\right)}\right)+{\mathbb{C}}_{\mathrm{a}}{\xi }_{\mathbf{x}}\left({\mathbf{u}}_{\mathrm{a}}^{\left(0\right)}\right)& \mathrm{in}{\Omega }_{\mathrm{a}},\hfill \end{array}$$

(34b)

$$\begin{array}{cc}\hfill {\mathsf{T}}_{\mathrm{b}}^{\left(0\right)}={\mathbb{C}}_{\mathrm{b}}{\xi }_{\mathbf{y}}\left({\mathbf{u}}_{\mathrm{b}}^{\left(1\right)}\right)+{\mathbb{C}}_{\mathrm{b}}{\xi }_{\mathbf{x}}\left({\mathbf{u}}_{\mathrm{b}}^{\left(0\right)}\right)& \mathrm{in}{\Omega }_{\mathrm{b}},\hfill \end{array}$$

(34c)

and finally the balance Equations (28a)–(28c) have coefficients of ${ϵ}^1$ given by

$$\begin{array}{cc}\hfill {\nabla }_{\mathbf{y}}·\left({\mathbb{C}}_{\mathrm{a}}{\xi }_{\mathbf{y}}\left({\mathbf{u}}_{\mathrm{a}}^{\left(1\right)}\right)\right)+{\nabla }_{\mathbf{y}}·\left({\mathbb{C}}_{\mathrm{a}}{\xi }_{\mathbf{x}}\left({\mathbf{u}}_{\mathrm{a}}^{\left(0\right)}\right)\right)+{\nabla }_{\mathbf{x}}·\left({\mathbb{C}}_{\mathrm{a}}{\xi }_{\mathbf{y}}\left({\mathbf{u}}_{\mathrm{a}}^{\left(0\right)}\right)\right)=0& \mathrm{in}{\Omega }_{\mathrm{a}},\hfill \end{array}$$

(35a)

$$\begin{array}{cc}\hfill {\nabla }_{\mathbf{y}}·\left({\mathbb{C}}_{\mathrm{b}}{\xi }_{\mathbf{y}}\left({\mathbf{u}}_{\mathrm{b}}^{\left(1\right)}\right)\right)+{\nabla }_{\mathbf{y}}·\left({\mathbb{C}}_{\mathrm{b}}{\xi }_{\mathbf{x}}\left({\mathbf{u}}_{\mathrm{b}}^{\left(0\right)}\right)\right)+{\nabla }_{\mathbf{x}}·\left({\mathbb{C}}_{\mathrm{b}}{\xi }_{\mathbf{y}}\left({\mathbf{u}}_{\mathrm{b}}^{\left(0\right)}\right)\right)=0& \mathrm{in}{\Omega }_{\mathrm{b}},\hfill \end{array}$$

(35b)

$$\begin{array}{cc}\hfill {\nabla }_{\mathbf{y}}^2{\mathbf{v}}^{\left(1\right)}+{\nabla }_{\mathbf{x}}·\left({\nabla }_{\mathbf{y}}{\mathbf{v}}^{\left(0\right)}\right)+{\nabla }_{\mathbf{y}}·\left({\nabla }_{\mathbf{x}}{\mathbf{v}}^{\left(0\right)}\right)={\nabla }_{\mathbf{y}}{p}^{\left(0\right)}& \mathrm{in}{\Omega }_{\mathrm{f}}.\hfill \end{array}$$

(35c)

We have from (31b), (31c), (32a) and (32b) and the periodicity conditions that ${\mathbf{u}}_{\mathrm{a}}^{\left(0\right)}$ and ${\mathbf{u}}_{\mathrm{b}}^{\left(0\right)}$ (the leading order solid displacements) are independent of the microscale variable $\mathbf{y}$. We can also then describe these as rigid body motions and therefore we can write

$$\begin{array}{cc}\hfill {\mathbf{u}}_{\mathrm{a}}^{\left(0\right)}=& {\mathbf{u}}_{\mathrm{a}}^{\left(0\right)}\left(\mathbf{x}\right)\hfill \end{array}$$

(36a)

$$\begin{array}{cc}\hfill {\mathbf{u}}_{\mathrm{b}}^{\left(0\right)}=& {\mathbf{u}}_{\mathrm{b}}^{\left(0\right)}\left(\mathbf{x}\right).\hfill \end{array}$$

(36b)

Due to the continuity of elastic displacements ${\mathbf{u}}_{\mathrm{a}}^{\left(0\right)}={\mathbf{u}}_{\mathrm{b}}^{\left(0\right)}$ on ${\Gamma }_{\mathrm{ab}}$ given by (30j), we can write

$${\mathbf{u}}^{\left(0\right)}={\mathbf{u}}_{\mathrm{a}}^{\left(0\right)}={\mathbf{u}}_{\mathrm{b}}^{\left(0\right)},$$

(37)

and we use this in the remainder of this manuscript. We can also see from (30e) and (30f) that

$$\begin{array}{c}\hfill i\omega {\mathbf{u}}_{\mathrm{a}}^{\left(0\right)}={\mathbf{v}}^{\left(0\right)},\end{array}$$

(38a)

$$\begin{array}{c}\hfill i\omega {\mathbf{u}}_{\mathrm{b}}^{\left(0\right)}={\mathbf{v}}^{\left(0\right)},\end{array}$$

(38b)

this means that

$${\mathbf{v}}^{\left(0\right)}={\mathbf{v}}^{\left(0\right)}\left(\mathbf{x}\right),$$

(39)

and then using (37), we have that

$${\mathbf{v}}^{\left(0\right)}=i\omega {\mathbf{u}}^{\left(0\right)}.$$

(40)

3.4. Microscale Problem

We now wish to form a problem from the equations that we have just derived in the previous section. We will use conditions (36a), (36b), (37) and (39) and we can take balance Equations (35a)–(35c), we have incompressibility constraint (33d), continuities of velocities (33e), (33f), continuity of stresses (30g)–(30i), with leading order stresses from (34a)–(34c) and continuity of displacements (33j), to form the problem the problem

$$\begin{array}{cc}\hfill {\nabla }_{\mathbf{y}}^2{\mathbf{v}}^{\left(1\right)}-{\nabla }_{\mathbf{y}}{p}^{\left(0\right)}=0& \mathrm{in}{\Omega }_{\mathrm{f}}\hfill \end{array}$$

(41a)

$$\begin{array}{cc}\hfill {\nabla }_{\mathbf{y}}·{\mathbf{v}}^{\left(1\right)}+{\nabla }_{\mathbf{x}}·{\mathbf{v}}^{\left(0\right)}=0& \mathrm{in}{\Omega }_{\mathrm{f}}\hfill \end{array}$$

(41b)

$$\begin{array}{cc}\hfill {\mathbf{v}}^{\left(1\right)}=i\omega {\mathbf{u}}_{\mathrm{a}}^{\left(1\right)}& \mathrm{on}{\Gamma }_{\mathrm{a}}\hfill \end{array}$$

(41c)

$$\begin{array}{cc}\hfill {\mathbf{v}}^{\left(1\right)}=i\omega {\mathbf{u}}_{\mathrm{b}}^{\left(1\right)}& \mathrm{on}{\Gamma }_{\mathrm{b}}\hfill \end{array}$$

(41d)

$$\begin{array}{cc}\hfill {\nabla }_{\mathbf{y}}·\left({\mathbb{C}}_{\mathrm{a}}{\xi }_{\mathbf{y}}\left({\mathbf{u}}_{\mathrm{a}}^{\left(1\right)}\right)\right)+{\nabla }_{\mathbf{y}}·\left({\mathbb{C}}_{\mathrm{a}}{\xi }_{\mathbf{x}}\left({\mathbf{u}}^{\left(0\right)}\right)\right)=0& \mathrm{in}{\Omega }_{\mathrm{a}}\hfill \end{array}$$

(41e)

$$\begin{array}{cc}\hfill {\nabla }_{\mathbf{y}}·\left({\mathbb{C}}_{\mathrm{b}}{\xi }_{\mathbf{y}}\left({\mathbf{u}}_{\mathrm{b}}^{\left(1\right)}\right)\right)+{\nabla }_{\mathbf{y}}·\left({\mathbb{C}}_{\mathrm{b}}{\xi }_{\mathbf{x}}\left({\mathbf{u}}^{\left(0\right)}\right)\right)=0& \mathrm{in}{\Omega }_{\mathrm{b}}\hfill \end{array}$$

(41f)

$$\begin{array}{cc}\hfill ({\mathbb{C}}_{\mathrm{a}}{\xi }_{\mathbf{y}}\left({\mathbf{u}}_{\mathrm{a}}^{\left(1\right)}\right)+{\mathbb{C}}_{\mathrm{a}}{\xi }_{\mathbf{x}}\left({\mathbf{u}}_{\mathrm{a}}^{\left(0\right)}\right)){\mathbf{n}}_{\mathrm{a}}=(-p^{\left(0\right)}\mathbf{I}+{\xi }_{\mathbf{x}}\left({\mathbf{v}}^{\left(0\right)}\right)+{\xi }_{\mathbf{y}}\left({\mathbf{v}}^{\left(1\right)}\right)){\mathbf{n}}_{\mathrm{a}}& \mathrm{on}{\Gamma }_{\mathrm{a}}\hfill \end{array}$$

(41g)

$$\begin{array}{cc}\hfill ({\mathbb{C}}_{\mathrm{b}}{\xi }_{\mathbf{y}}\left({\mathbf{u}}_{\mathrm{b}}^{\left(1\right)}\right)+{\mathbb{C}}_{\mathrm{b}}{\xi }_{\mathbf{x}}\left({\mathbf{u}}_{\mathrm{b}}^{\left(0\right)}\right)){\mathbf{n}}_{\mathrm{b}}=(-p^{\left(0\right)}\mathbf{I}+{\xi }_{\mathbf{x}}\left({\mathbf{v}}^{\left(0\right)}\right)+{\xi }_{\mathbf{y}}\left({\mathbf{v}}^{\left(1\right)}\right)){\mathbf{n}}_{\mathrm{b}}& \mathrm{on}{\Gamma }_{\mathrm{b}}\hfill \end{array}$$

(41h)

$$\begin{array}{cc}\hfill {\mathbb{C}}_{\mathrm{a}}{\xi }_{\mathbf{y}}\left({\mathbf{u}}_{\mathrm{a}}^{\left(1\right)}\right){\mathbf{n}}_{\mathrm{ab}}-{\mathbb{C}}_{\mathrm{b}}{\xi }_{\mathbf{y}}\left({\mathbf{u}}_{\mathrm{b}}^{\left(1\right)}\right){\mathbf{n}}_{\mathrm{ab}}=({\mathbb{C}}_{\mathrm{b}}-{\mathbb{C}}_{\mathrm{a}}){\xi }_{\mathbf{x}}\left({\mathbf{u}}^{\left(0\right)}\right){\mathbf{n}}_{\mathrm{ab}}& \mathrm{on}{\Gamma }_{\mathrm{ab}}\hfill \end{array}$$

(41i)

$$\begin{array}{cc}\hfill {\mathbf{u}}_{\mathrm{a}}^{\left(1\right)}={\mathbf{u}}_{\mathrm{b}}^{\left(1\right)}& \mathrm{on}{\Gamma }_{\mathrm{ab}}\hfill \end{array}$$

(41j)

Since we have condition (40), we can replace the ${\mathbf{v}}^{\left(0\right)}$ terms in Equations (41b), (41g) and (41h) with $i\omega {\mathbf{u}}^{\left(0\right)}$. That is,

$$\begin{array}{cc}\hfill {\nabla }_{\mathbf{y}}^2{\mathbf{v}}^{\left(1\right)}-{\nabla }_{\mathbf{y}}{p}^{\left(0\right)}=0& \mathrm{in}{\Omega }_{\mathrm{f}}\hfill \end{array}$$

(42a)

$$\begin{array}{cc}\hfill {\nabla }_{\mathbf{y}}·{\mathbf{v}}^{\left(1\right)}+i\omega {\nabla }_{\mathbf{x}}·{\mathbf{u}}^{\left(0\right)}=0& \mathrm{in}{\Omega }_{\mathrm{f}}\hfill \end{array}$$

(42b)

$$\begin{array}{cc}\hfill {\mathbf{v}}^{\left(1\right)}=i\omega {\mathbf{u}}_{\mathrm{a}}^{\left(1\right)}& \mathrm{on}{\Gamma }_{\mathrm{a}}\hfill \end{array}$$

(42c)

$$\begin{array}{cc}\hfill {\mathbf{v}}^{\left(1\right)}=i\omega {\mathbf{u}}_{\mathrm{b}}^{\left(1\right)}& \mathrm{on}{\Gamma }_{\mathrm{b}}\hfill \end{array}$$

(42d)

$$\begin{array}{cc}\hfill {\nabla }_{\mathbf{y}}·\left({\mathbb{C}}_{\mathrm{a}}{\xi }_{\mathbf{y}}\left({\mathbf{u}}_{\mathrm{a}}^{\left(1\right)}\right)\right)+{\nabla }_{\mathbf{y}}·\left({\mathbb{C}}_{\mathrm{a}}{\xi }_{\mathbf{x}}\left({\mathbf{u}}^{\left(0\right)}\right)\right)=0& \mathrm{in}{\Omega }_{\mathrm{a}}\hfill \end{array}$$

(42e)

$$\begin{array}{cc}\hfill {\nabla }_{\mathbf{y}}·\left({\mathbb{C}}_{\mathrm{b}}{\xi }_{\mathbf{y}}\left({\mathbf{u}}_{\mathrm{b}}^{\left(1\right)}\right)\right)+{\nabla }_{\mathbf{y}}·\left({\mathbb{C}}_{\mathrm{b}}{\xi }_{\mathbf{x}}\left({\mathbf{u}}^{\left(0\right)}\right)\right)=0& \mathrm{in}{\Omega }_{\mathrm{b}}\hfill \end{array}$$

(42f)

$$\begin{array}{cc}\hfill ({\mathbb{C}}_{\mathrm{a}}{\xi }_{\mathbf{y}}\left({\mathbf{u}}_{\mathrm{a}}^{\left(1\right)}\right)+{\mathbb{C}}_{\mathrm{a}}{\xi }_{\mathbf{x}}\left({\mathbf{u}}_{\mathrm{a}}^{\left(0\right)}\right)){\mathbf{n}}_{\mathrm{a}}=(-p^{\left(0\right)}\mathbf{I}+i\omega {\xi }_{\mathbf{x}}\left({\mathbf{u}}^{\left(0\right)}\right)+{\xi }_{\mathbf{y}}\left({\mathbf{v}}^{\left(1\right)}\right)){\mathbf{n}}_{\mathrm{a}}& \mathrm{on}{\Gamma }_{\mathrm{a}}\hfill \end{array}$$

(42g)

$$\begin{array}{cc}\hfill ({\mathbb{C}}_{\mathrm{b}}{\xi }_{\mathbf{y}}\left({\mathbf{u}}_{\mathrm{b}}^{\left(1\right)}\right)+{\mathbb{C}}_{\mathrm{b}}{\xi }_{\mathbf{x}}\left({\mathbf{u}}_{\mathrm{b}}^{\left(0\right)}\right)){\mathbf{n}}_{\mathrm{b}}=(-p^{\left(0\right)}\mathbf{I}+i\omega {\xi }_{\mathbf{x}}\left({\mathbf{u}}^{\left(0\right)}\right)+{\xi }_{\mathbf{y}}\left({\mathbf{v}}^{\left(1\right)}\right)){\mathbf{n}}_{\mathrm{b}}& \mathrm{on}{\Gamma }_{\mathrm{b}}\hfill \end{array}$$

(42h)

$$\begin{array}{cc}\hfill {\mathbb{C}}_{\mathrm{a}}{\xi }_{\mathbf{y}}\left({\mathbf{u}}_{\mathrm{a}}^{\left(1\right)}\right){\mathbf{n}}_{\mathrm{ab}}-{\mathbb{C}}_{\mathrm{b}}{\xi }_{\mathbf{y}}\left({\mathbf{u}}_{\mathrm{b}}^{\left(1\right)}\right){\mathbf{n}}_{\mathrm{ab}}=({\mathbb{C}}_{\mathrm{b}}-{\mathbb{C}}_{\mathrm{a}}){\xi }_{\mathbf{x}}\left({\mathbf{u}}^{\left(0\right)}\right){\mathbf{n}}_{\mathrm{ab}}& \mathrm{on}{\Gamma }_{\mathrm{ab}}\hfill \end{array}$$

(42i)

$$\begin{array}{cc}\hfill {\mathbf{u}}_{\mathrm{a}}^{\left(1\right)}={\mathbf{u}}_{\mathrm{b}}^{\left(1\right)}& \mathrm{on}{\Gamma }_{\mathrm{ab}}\hfill \end{array}$$

(42j)

We now can exploit the linearity of the system (41a)–(41j) to write the following ansatz

$$\begin{array}{cc}\hfill & {\mathbf{v}}^{\left(1\right)}=i\omega \mathcal{A}{\xi }_{\mathbf{x}}{\mathbf{u}}^{\left(0\right)},\hfill \end{array}$$

(43a)

$$\begin{array}{cc}\hfill & {\mathbf{u}}_{\mathrm{a}}^{\left(1\right)}={\mathcal{B}}_{\mathrm{a}}{\xi }_{\mathbf{x}}{\mathbf{u}}^{\left(0\right)},\hfill \end{array}$$

(43b)

$$\begin{array}{cc}\hfill & {\mathbf{u}}_{\mathrm{b}}^{\left(1\right)}={\mathcal{B}}_{\mathrm{b}}{\xi }_{\mathbf{x}}{\mathbf{u}}^{\left(0\right)},\hfill \end{array}$$

(43c)

$$\begin{array}{cc}\hfill & p^{\left(0\right)}=i\omega \mathsf{P}:{\xi }_{\mathbf{x}}{\mathbf{u}}^{\left(0\right)}=\mathrm{Tr}\left(i\omega \mathsf{P}{\xi }_{\mathbf{x}}{\mathbf{u}}^{\left(0\right)}\right),\hfill \end{array}$$

(43d)

where $\mathcal{A}$, ${\mathcal{B}}_{\mathrm{a}}$ and ${\mathcal{B}}_{\mathrm{b}}$ are third rank tensors and $\mathsf{P}$ is a second rank tensor. The auxiliary fields $\mathcal{A}$, ${\mathcal{B}}_{\mathrm{a}}$, ${\mathcal{B}}_{\mathrm{b}}$ and $\mathsf{P}$ solve the following cell problem

$$\begin{array}{cc}\hfill {\nabla }_{\mathbf{y}}^2\mathcal{A}-{\nabla }_{\mathbf{y}}\mathsf{P}=0& \hfill \end{array}$$

(44a)

$$\begin{array}{cc}\hfill {\nabla }_{\mathbf{y}}·{\mathcal{A}}^{\mathrm{T}}+\mathsf{I}=0& \hfill \end{array}$$

(44b)

$$\begin{array}{cc}\hfill \mathcal{A}={\mathcal{B}}_{\mathrm{a}}& \hfill \end{array}$$

(44c)

$$\begin{array}{cc}\hfill \mathcal{A}={\mathcal{B}}_{\mathrm{b}}& \hfill \end{array}$$

(44d)

$$\begin{array}{cc}\hfill {\nabla }_{\mathbf{y}}·\left({\mathbb{C}}_{\mathrm{a}}{\xi }_{\mathbf{y}}\left({\mathcal{B}}_{\mathrm{a}}\right)\right)+{\nabla }_{\mathbf{y}}·{\mathbb{C}}_{\mathrm{a}}=0& \hfill \end{array}$$

(44e)

$$\begin{array}{cc}\hfill {\nabla }_{\mathbf{y}}·\left({\mathbb{C}}_{\mathrm{b}}{\xi }_{\mathbf{y}}\left({\mathcal{B}}_{\mathrm{b}}\right)\right)+{\nabla }_{\mathbf{y}}·{\mathbb{C}}_{\mathrm{b}}=0& \hfill \end{array}$$

(44f)

$$\begin{array}{cc}\hfill ({\mathbb{C}}_{\mathrm{a}}{\xi }_{\mathbf{y}}\left({\mathcal{B}}_{\mathrm{a}}\right)+{\mathbb{C}}_{\mathrm{a}}){\mathbf{n}}_{\mathrm{a}}=(-i\omega \mathsf{I}\otimes \mathsf{P}+i\omega \mathbb{I}+i\omega {\xi }_{\mathbf{y}}\mathcal{A}){\mathbf{n}}_{\mathrm{a}}& \hfill \end{array}$$

(44g)

$$\begin{array}{cc}\hfill ({\mathbb{C}}_{\mathrm{b}}{\xi }_{\mathbf{y}}\left({\mathcal{B}}_{\mathrm{b}}\right)+{\mathbb{C}}_{\mathrm{b}}){\mathbf{n}}_{\mathrm{a}}=(-i\omega \mathsf{I}\otimes \mathsf{P}+i\omega \mathbb{I}+i\omega {\xi }_{\mathbf{y}}\mathcal{A}){\mathbf{n}}_{\mathrm{b}}& \hfill \end{array}$$

(44h)

$$\begin{array}{cc}\hfill ({\mathbb{C}}_{\mathrm{a}}{\xi }_{\mathbf{y}}\left({\mathcal{B}}_{\mathrm{a}}\right)-{\mathbb{C}}_{\mathrm{b}}{\xi }_{\mathbf{y}}\left({\mathcal{B}}_{\mathrm{b}}\right)){\mathbf{n}}_{\mathrm{ab}}=({\mathbb{C}}_{\mathrm{b}}-{\mathbb{C}}_{\mathrm{a}}){\mathbf{n}}_{\mathrm{ab}}& \hfill \end{array}$$

(44i)

$$\begin{array}{cc}\hfill {\mathcal{B}}_{\mathrm{a}}={\mathcal{B}}_{\mathrm{b}}& \hfill \end{array}$$

(44j)

where we can define $\mathbb{I}={\delta }_{ij}{\delta }_{kl}$ as a fourth rank identity tensor.

Problem (44a)–(44j) is to be solved on the periodic cell and will also require periodic conditions on the external boundaries of $\Omega$. In order to obtain a unique solution, we place one further condition on the auxiliary variables $\mathcal{A}$, ${\mathcal{B}}_{\mathrm{a}}$, ${\mathcal{B}}_{\mathrm{b}}$ and $\mathsf{P}$, that is

$${\langle \mathcal{A}\rangle }_{\mathrm{f}}=0,{\langle {\mathcal{B}}_{\mathrm{a}}\rangle }_{\mathrm{a}}=0,{\langle {\mathcal{B}}_{\mathrm{b}}\rangle }_{\mathrm{b}}=0,{\langle \mathsf{P}\rangle }_{\mathrm{f}}=0$$

(45)

3.5. Balance Equation on the Macroscale

We must derive the macroscale balance equation. We apply the integral averages to Equations (33a)–(33c) and sum up to obtain

$$\begin{array}{cc}\hfill & {\int }_{{\Omega }_{\mathrm{a}}}{\nabla }_{\mathbf{y}}·{\mathsf{T}}_{\mathrm{a}}^{\left(1\right)}d\mathbf{y}+{\int }_{{\Omega }_{\mathrm{b}}}{\nabla }_{\mathbf{y}}·{\mathsf{T}}_{\mathrm{b}}^{\left(1\right)}d\mathbf{y}+{\int }_{{\Omega }_{\mathrm{f}}}{\nabla }_{\mathbf{y}}·{{\mathsf{T}}_{\mathsf{f}}}^{\left(1\right)}d\mathbf{y}+\hfill \\ \hfill & {\int }_{{\Omega }_{\mathrm{a}}}{\nabla }_{\mathbf{x}}·{\mathsf{T}}_{\mathrm{a}}^{\left(0\right)}d\mathbf{y}+{\int }_{{\Omega }_{\mathrm{b}}}{\nabla }_{\mathbf{x}}·{\mathsf{T}}_{\mathrm{b}}^{\left(0\right)}d\mathbf{y}+{\int }_{{\Omega }_{\mathrm{f}}}{\nabla }_{\mathbf{x}}·{{\mathsf{T}}_{\mathsf{f}}}^{\left(0\right)}d\mathbf{y}=0.\hfill \end{array}$$

(46)

The divergence theorem can be applied to the first three integrals and we can use macroscopic uniformity condition (25) to rearrange the last three integrals so that we obtain

$$\begin{array}{cc}\hfill & {\int }_{\partial {\Omega }_{\mathrm{a}}\setminus {\Gamma }_{\mathrm{a}}\cup {\Gamma }_{\mathrm{ab}}}{\mathsf{T}}_{\mathrm{a}}^{\left(1\right)}{\mathbf{n}}_{{\Omega }_{\mathrm{a}}\setminus {\Gamma }_{\mathrm{a}}\cup {\Gamma }_{\mathrm{ab}}}\mathrm{dS}+{\int }_{{\Gamma }_{\mathrm{a}}}{\mathsf{T}}_{\mathrm{a}}^{\left(1\right)}{\mathbf{n}}_{\mathrm{a}}\mathrm{dS}-{\int }_{{\Gamma }_{\mathrm{ab}}}{\mathsf{T}}_{\mathrm{a}}^{\left(1\right)}{\mathbf{n}}_{\mathrm{ab}}\mathrm{dS}+\hfill \\ \hfill & {\int }_{\partial {\Omega }_{\mathrm{b}}\setminus {\Gamma }_{\mathrm{b}}\cup {\Gamma }_{\mathrm{ab}}}{\mathsf{T}}_{\mathrm{b}}^{\left(1\right)}{\mathbf{n}}_{{\Omega }_{\mathrm{b}}\setminus {\Gamma }_{\mathrm{b}}\cup {\Gamma }_{\mathrm{ab}}}\mathrm{dS}+{\int }_{{\Gamma }_{\mathrm{b}}}{\mathsf{T}}_{\mathrm{b}}^{\left(1\right)}{\mathbf{n}}_{\mathrm{b}}\mathrm{dS}+{\int }_{{\Gamma }_{\mathrm{ab}}}{\mathsf{T}}_{\mathrm{b}}^{\left(1\right)}{\mathbf{n}}_{\mathrm{ab}}\mathrm{dS}+\hfill \\ \hfill & {\int }_{\partial {\Omega }_{\mathrm{f}}\setminus {\Gamma }_{\mathrm{a}}\cup {\Gamma }_{\mathrm{b}}}{\mathsf{T}}_{\mathsf{f}}^{\left(1\right)}{\mathbf{n}}_{{\Omega }_{\mathrm{f}}\setminus {\Gamma }_{\mathrm{a}}\cup {\Gamma }_{\mathrm{b}}}\mathrm{dS}-{\int }_{{\Gamma }_{\mathrm{b}}}{\mathsf{T}}_{\mathsf{f}}^{\left(1\right)}{\mathbf{n}}_{\mathrm{b}}\mathrm{dS}-{\int }_{{\Gamma }_{\mathrm{a}}}{\mathsf{T}}_{\mathsf{f}}^{\left(1\right)}{\mathbf{n}}_{\mathrm{a}}\mathrm{dS}+\hfill \\ \hfill & {\nabla }_{\mathbf{x}}·{\int }_{{\Omega }_{\mathrm{a}}}{\mathsf{T}}_{\mathrm{a}}^{\left(0\right)}d\mathbf{y}+{\nabla }_{\mathbf{x}}·{\int }_{{\Omega }_{\mathrm{b}}}{\mathsf{T}}_{\mathrm{b}}^{\left(0\right)}dy+{\nabla }_{\mathbf{x}}·{\int }_{{\Omega }_{\mathrm{f}}}{{\mathsf{T}}_{\mathsf{f}}}^{\left(0\right)}d\mathbf{y}=0,\hfill \end{array}$$

(47)

where the vectors ${\mathbf{n}}_{\mathrm{a}}$, ${\mathbf{n}}_{\mathrm{b}}$, ${\mathbf{n}}_{\mathrm{ab}}$, ${\mathbf{n}}_{{\Omega }_{\mathrm{I}}\setminus {\Gamma }_{\mathrm{a}}\cup {\Gamma }_{\mathrm{ab}}}$, ${\mathbf{n}}_{{\Omega }_{\mathrm{b}}\setminus {\Gamma }_{\mathrm{b}}\cup {\Gamma }_{\mathrm{ab}}}$ and ${\mathbf{n}}_{{\Omega }_{\mathrm{f}}\setminus {\Gamma }_{\mathrm{a}}\cup {\Gamma }_{\mathrm{b}}}$ are the unit normals to the interfaces ${\Gamma }_{\mathrm{a}}$, ${\Gamma }_{\mathrm{b}}$, ${\Gamma }_{\mathrm{ab}}$, $\partial {\Omega }_{\mathrm{a}}\setminus {\Gamma }_{\mathrm{a}}\cup {\Gamma }_{\mathrm{ab}}$, $\partial {\Omega }_{\mathrm{b}}\setminus {\Gamma }_{\mathrm{b}}\cup {\Gamma }_{\mathrm{ab}}$ and $\partial {\Omega }_{\mathrm{f}}\setminus {\Gamma }_{\mathrm{a}}\cup {\Gamma }_{\mathrm{b}}$. Due to $\mathbf{y}$-periodicity, we can cancel the terms that arise on the external boundaries of the phases ${\Omega }_{\mathrm{a}}$, ${\Omega }_{\mathrm{b}}$ and ${\Omega }_{\mathrm{f}}$. This gives

$$\begin{array}{cc}\hfill & {\int }_{{\Gamma }_{\mathrm{a}}}{\mathsf{T}}_{\mathrm{a}}^{\left(1\right)}{\mathbf{n}}_{\mathrm{a}}\mathrm{dS}+{\int }_{{\Gamma }_{\mathrm{b}}}{\mathsf{T}}_{\mathrm{b}}^{\left(1\right)}{\mathbf{n}}_{\mathrm{b}}\mathrm{dS}-{\int }_{{\Gamma }_{\mathrm{a}}}{\mathsf{T}}_{\mathsf{f}}^{\left(1\right)}{\mathbf{n}}_{\mathrm{a}}\mathrm{dS}-\hfill \\ \hfill & {\int }_{{\Gamma }_{\mathrm{b}}}{\mathsf{T}}_{\mathsf{f}}^{\left(1\right)}{\mathbf{n}}_{\mathrm{b}}\mathrm{dS}-{\int }_{{\Gamma }_{\mathrm{ab}}}{\mathsf{T}}_{\mathrm{a}}^{\left(1\right)}{\mathbf{n}}_{\mathrm{ab}}\mathrm{dS}+{\int }_{{\Gamma }_{\mathrm{ab}}}{\mathsf{T}}_{\mathrm{b}}^{\left(1\right)}{\mathbf{n}}_{\mathrm{ab}}\mathrm{dS}+\hfill \\ \hfill & {\nabla }_{\mathbf{x}}·{\int }_{{\Omega }_{\mathrm{I}}}{\mathsf{T}}_{\mathrm{a}}^{\left(0\right)}\mathrm{d}\mathbf{y}+{\nabla }_{\mathbf{x}}·{\int }_{{\Omega }_{\mathrm{b}}}{\mathsf{T}}_{\mathrm{b}}^{\left(0\right)}\mathrm{d}\mathbf{y}+{\nabla }_{\mathbf{x}}·{\int }_{{\Omega }_{\mathrm{f}}}{{\mathsf{T}}_{\mathsf{f}}}^{\left(0\right)}\mathrm{d}\mathbf{y}=0.\hfill \end{array}$$

(48)

The continuity of stresses interface conditions, Equations (33g)–(33i), can be used to cancel out the first six integrals in (48). This means that the remaining terms can be written as

$$\begin{array}{cc}\hfill & {\nabla }_{\mathbf{x}}·({\langle {\mathsf{T}}_{\mathrm{a}}^{\left(0\right)}\rangle }_{\mathrm{a}}+{\langle {\mathsf{T}}_{\mathrm{b}}^{\left(0\right)}\rangle }_{\mathrm{b}}+{\langle {\mathsf{T}}_{\mathrm{f}}^{\left(0\right)}\rangle }_{\mathrm{f}})=0.\hfill \end{array}$$

(49)

Our balance equation comprises the zero-th order solid stress tensors. Using Equations (34b) and (34c), we know that ${\mathbf{u}}_{\mathrm{a}}^{\left(1\right)}$ and ${\mathbf{u}}_{\mathrm{b}}^{\left(1\right)}$ appear in ${\mathsf{T}}_{\mathrm{a}}^{\left(0\right)}$ and ${\mathsf{T}}_{\mathrm{b}}^{\left(0\right)}$, respectively. Therefore, we write the leading order solid stresses as

$${\mathsf{T}}_{\mathrm{a}}^{\left(0\right)}={\mathbb{C}}_{\mathrm{a}}{\mathbb{M}}_{\mathrm{a}}{\xi }_{\mathbf{x}}\left({\mathbf{u}}^{\left(0\right)}\right)+{\mathbb{C}}_{\mathrm{a}}{\xi }_{\mathbf{x}}\left({\mathbf{u}}^{\left(0\right)}\right)$$

(50)

and

$${\mathsf{T}}_{\mathrm{b}}^{\left(0\right)}={\mathbb{C}}_{\mathrm{b}}{\mathbb{M}}_{\mathrm{b}}{\xi }_{\mathbf{x}}\left({\mathbf{u}}^{\left(0\right)}\right)+{\mathbb{C}}_{\mathrm{b}}{\xi }_{\mathbf{x}}\left({\mathbf{u}}^{\left(0\right)}\right)$$

(51)

making use of the notation

$$\begin{array}{c}\hfill {\mathbb{M}}_{\mathrm{a}}:={\xi }_{\mathbf{y}}\left({\mathcal{B}}_{\mathrm{a}}\right),{\mathbb{M}}_{\mathrm{b}}:={\xi }_{\mathbf{y}}\left({\mathcal{B}}_{\mathrm{b}}\right).\end{array}$$

(52)

We also want to write down the zero-th order fluid stress (34a) since we have an expression for ${\mathbf{v}}^{\left(1\right)}$

$${{\mathsf{T}}_{\mathsf{f}}}^{\left(0\right)}=-p^{\left(0\right)}\mathbf{I}+i\omega (\mathbb{L}+\mathbb{I}){\xi }_{\mathbf{x}}{\mathbf{u}}^{\left(0\right)},$$

(53)

or equivalently

$${{\mathsf{T}}_{\mathsf{f}}}^{\left(0\right)}=-p^{\left(0\right)}\mathbf{I}+(\mathbb{L}+\mathbb{I}){\xi }_{\mathbf{x}}{\mathbf{v}}^{\left(0\right)},$$

(54)

where we have defined

$$\mathbb{L}:={\xi }_{\mathbf{y}}\mathcal{A}.$$

(55)

Adding (50), (51) and (53) and applying the integral average over each of the domains gives

$$\begin{array}{cc}\hfill {\langle {\mathsf{T}}_{\mathrm{a}}^{\left(0\right)}\rangle }_{\mathrm{a}}+{\langle {\mathsf{T}}_{\mathrm{b}}^{\left(0\right)}\rangle }_{\mathrm{b}}+{\langle {\mathsf{T}}_{\mathrm{f}}^{\left(0\right)}\rangle }_{\mathrm{f}}=& ({\langle {\mathbb{C}}_{\mathrm{a}}{\mathbb{M}}_{\mathrm{a}}+{\mathbb{C}}_{\mathrm{a}}\rangle }_{\mathrm{a}}+{\langle {\mathbb{C}}_{\mathrm{b}}{\mathbb{M}}_{\mathrm{b}}+{\mathbb{C}}_{\mathrm{b}}\rangle }_{\mathrm{b}}){\xi }_{\mathbf{x}}\left({\mathbf{u}}^{\left(0\right)}\right)\hfill \\ \hfill & +{\langle (\mathbb{L}+\mathbb{I}-\left(\mathrm{Tr}P\right)\mathbb{I})\rangle }_{\mathrm{f}}{\xi }_{\mathbf{x}}\left({\mathbf{v}}^{\left(0\right)}\right)\hfill \end{array}$$

(56)

From (49) we know that

$${\nabla }_{\mathbf{x}}·{\mathsf{T}}_{\mathrm{Eff}}=0$$

(57)

with

$${\mathsf{T}}_{\mathrm{Eff}}={\langle {\mathsf{T}}_{\mathrm{a}}^{\left(0\right)}\rangle }_{\mathrm{a}}+{\langle {\mathsf{T}}_{\mathrm{b}}^{\left(0\right)}\rangle }_{\mathrm{b}}+{\langle {\mathsf{T}}_{\mathrm{f}}^{\left(0\right)}\rangle }_{\mathrm{f}}$$

(58)

We are able to describe (57) and (58) as the averaged force balance equation for our material.

4. The Macroscale Governing Equations and Limit Cases

We now have derived the equations necessary to write the macroscale governing equations for a linear viscoelastic composite material. That is

$$\begin{array}{cc}\hfill & {\nabla }_{\mathbf{x}}·{\mathsf{T}}_{\mathrm{Eff}}=0,\hfill \end{array}$$

(59a)

$$\begin{array}{cc}\hfill & {\mathsf{T}}_{\mathrm{Eff}}=\left({\langle {\mathbb{C}}_{\mathrm{a}}{\mathbb{M}}_{\mathrm{a}}+{\mathbb{C}}_{\mathrm{a}}\rangle }_{\mathrm{a}}+{\langle {\mathbb{C}}_{\mathrm{b}}{\mathbb{M}}_{\mathrm{b}}+{\mathbb{C}}_{\mathrm{b}}\rangle }_{\mathrm{b}}\right){\xi }_{\mathbf{x}}\left({\mathbf{u}}^{\left(0\right)}\right)+{\langle (\mathbb{L}+\mathbb{I}-\left(\mathrm{Tr}P\right)\mathbb{I})\rangle }_{\mathrm{f}}{\xi }_{\mathbf{x}}\left({\mathbf{v}}^{\left(0\right)}\right),\hfill \end{array}$$

(59b)

where ${\mathbf{u}}^{\left(0\right)}$ is the leading order solid displacement, ${\mathbf{v}}^{\left(0\right)}$ is the leading order fluid velocity and $p^{\left(0\right)}$ is the leading order pressure. Equation (59a) is the balance equation with the new constitutive law for viscoelastic composites given by (59b). We can see that the viscoelastic constitutive law takes exactly the form that is expected of a Kelvin–Voigt viscoelastic material comprising first the elastic constitutive relationship and the the second part is the viscous part of the relation as in Kelvin–Voigt materials. The addition of the multiple elastic phases being encoded in our model influences the elastic part of our constitutive law.

The new model we derive has an important difference from previous formulations of standard viscoelastic materials. Our model has the ability to incorporate multiple elastic phases all interplaying with the fluid and each other, whereas previous viscoelastic formulations of this kind, such as [17], only consider one fluid and one solid phase. These additional interactions between the multiple phases can be extremely useful in physical systems where the solid component is rarely homogeneous. The ability to model heterogeneous materials comes from the fact that discrepancies in the elastic and mechanical properties of each phase are accounted for by the multiple elasticity tensors ${\mathbb{C}}_{\mathrm{a}}$ and ${\mathbb{C}}_{\mathrm{b}}$ as well as the tensors ${\mathbb{M}}_{\mathrm{a}}$ and ${\mathbb{M}}_{\mathrm{b}}$ which account for the differences in the elastic properties at different points in the microstructure. In this work, we propose the novel cell problem (44a)–(44j), this is the problem from which the model coefficients are determined. The cell problem is an extension to the problem found for viscoelastic materials in [17] and comprises some of the elastic problem associated with poroelastic composites [38]. This cell problem comprises both the fluid and solid equations in one problem and therefore we do not have the decoupling of the different phases as seen in poroelastic-type cell problems. The model as it stands can be described as a comprehensive framework for Kelvin–Voigt viscoelastic materials that comprise various elastic phases.

We now wish to understand how our macroscale model for viscoelastic composites (59a) and (59b) compares and can reduce to previous models in the literature. We consider the viscoelastic model derived via asymptotic homogenization in [17] which considers only one elastic phase and the fluid, and we consider the model for poroelastic composites by [38] which addresses the interaction of a porous matrix and and elastic phase where fluid flows in the pores which are also comparable in size to the distance between the inclusions.

Remark 6

(Comparison with Burridge and Keller [17]). We now wish to compare the model we derived with the results of the [17] in the remark of effective viscoelasticity. There, the authors consider one elastic phase and one fluid phase and the interactions between them. We begin the comparison by rewriting the fluid stress ${{\mathsf{T}}_{\mathsf{f}}}^{\left(0\right)}$, from (53), using both the expression for ${\mathbf{v}}^{\left(1\right)}$ and $p^{\left(0\right)}$. That is,

$${{\mathsf{T}}_{\mathsf{f}}}^{\left(0\right)}=i\omega (\mathbb{L}+\mathbb{I}-\left(\mathrm{Tr}\mathsf{P}\right)\mathbb{I}){\xi }_{\mathbf{x}}{\mathbf{u}}^{\left(0\right)}.$$

(60)

We can now use this version of the fluid stress in (56) to obtain

$$\begin{array}{cc}\hfill {\langle {\mathsf{T}}_{\mathrm{a}}^{\left(0\right)}\rangle }_{\mathrm{a}}+{\langle {\mathsf{T}}_{\mathrm{b}}^{\left(0\right)}\rangle }_{\mathrm{b}}+{\langle {\mathsf{T}}_{\mathrm{f}}^{\left(0\right)}\rangle }_{\mathrm{f}}& =({\langle {\mathbb{C}}_{\mathrm{a}}{\mathbb{M}}_{\mathrm{a}}+{\mathbb{C}}_{\mathrm{a}}\rangle }_{\mathrm{a}}+{\langle {\mathbb{C}}_{\mathrm{b}}{\mathbb{M}}_{\mathrm{b}}+{\mathbb{C}}_{\mathrm{b}}\rangle }_{\mathrm{b}}\hfill \\ \hfill & +{\langle i\omega (\mathbb{L}+\mathbb{I}-\left(\mathrm{Tr}\mathsf{P}\right)\mathbb{I})\rangle }_{\mathrm{f}}){\xi }_{\mathbf{x}}\left({\mathbf{u}}^{\left(0\right)}\right).\hfill \end{array}$$

(61)

This is still the constitutive law of our model, it had just been presented in a different form. We then rewrite the macroscale model (59a) and (59b) as

$$\begin{array}{cc}\hfill & {\nabla }_{\mathbf{x}}·{\mathsf{T}}_{\mathrm{Eff}}=0,\hfill \end{array}$$

(62a)

$$\begin{array}{cc}\hfill & {\mathsf{T}}_{\mathrm{Eff}}=\left({\langle {\mathbb{C}}_{\mathrm{a}}{\mathbb{M}}_{\mathrm{a}}+{\mathbb{C}}_{\mathrm{a}}\rangle }_{\mathrm{a}}+{\langle {\mathbb{C}}_{\mathrm{b}}{\mathbb{M}}_{\mathrm{b}}+{\mathbb{C}}_{\mathrm{b}}\rangle }_{\mathrm{b}}+{\langle i\omega (\mathbb{L}+\mathbb{I}-\left(\mathrm{Tr}\mathsf{P}\right)\mathbb{I})\rangle }_{\mathrm{f}}\right){\xi }_{\mathbf{x}}\left({\mathbf{u}}^{\left(0\right)}\right).\hfill \end{array}$$

(62b)

The constitutive equation comprises the average of the stresses in both the solid domains $\mathrm{I}$ and $\mathrm{II}$ and fluid domain of the material. In [17], the constitutive law is written as the sum of the leading order stress in the solid and the leading order stress in the fluid. If we assumed we had indeed only one solid phase, our constitutive law will match that of [17]. That is,

$${\mathsf{T}}_{\mathrm{Eff}}=\left({\langle {\mathbb{C}}_{\mathrm{a}}{\mathbb{M}}_{\mathrm{a}}+{\mathbb{C}}_{\mathrm{a}}\rangle }_{\mathrm{a}}+{\langle i\omega (\mathbb{L}+\mathbb{I}-\left(\mathrm{Tr}\mathsf{P}\right)\mathbb{I})\rangle }_{\mathrm{f}}\right){\xi }_{\mathbf{x}}\left({\mathbf{u}}^{\left(0\right)}\right),$$

(63)

and we can identify ${\mathsf{T}}_{\mathrm{Eff}}={\langle {\mathsf{T}}_{\mathrm{a}}^{\left(0\right)}\rangle }_{\mathrm{a}}+{\langle {\mathsf{T}}_{\mathrm{f}}^{\left(0\right)}\rangle }_{\mathrm{f}}$ with the notation $\overline{\mathsf{\tau }}={\overline{\mathsf{\tau }}}_0+{\overline{\sigma }}_0$ which has been used in [17]. In this work, we find explicit forms of the leading order stresses by proposing a solution to the problem (41a)–(41j). In [17], the suggested equations to form a linear dissipative problem are the same as those we have in (41a)–(41j), with the exception that we have two solid elastic phases and they consider only one. In their work, they do not explicitly solve this problem, but propose that the leading order stresses in the solid and fluid are proportional to the macroscale gradient of the leading order solid displacement ${\mathbf{u}}^{\left(0\right)}$ as we have done here.

Remark 7

(Comparison with Poroelastic Composites [38]). We now can consider each of the differences between this model and that of [38] for poroelastic composites. In the work of [38], the authors consider a porous matrix with embedded elastic subphases that are in contact with the matrix and the fluid that flows in connected cylindrical pores. Due to the profile of the fluid flowing in the pores, [38] have that the fluid velocity scales by $\mathbf{v}=\frac{C{d}^2}{\mu }{\mathbf{v}}^{\prime }$, where d is the radius of the pores, which is the standard scaling for Stokes’ flow in porous media. In this current work we have that the pores are much larger and therefore scale by $\mathbf{v}=\frac{C{L}^2}{\mu }{\mathbf{v}}^{\prime }$. These two different choices relate to the observed microstructure and result in the appearance of an ${ϵ}^2$ coefficient in the fluid stress tensor and in the Stokes equation in porous media, but not in the case of viscoelastic media. This means that when applying the asymptotic homogenization method, the Stokes-type equation and the fluid stress tensor will have different orders and terms than those that we equate in this work. That is, the Stokes equation in [38] reads

$${\nabla }_{\mathbf{y}}^2{\mathbf{v}}^{\left(0\right)}-{\nabla }_{\mathbf{y}}{p}^{\left(1\right)}={\nabla }_{\mathbf{x}}{p}^{\left(0\right)},$$

(64)

we can see that the orders of ${\nabla }_{\mathbf{y}}^2\mathbf{v}$ and ${\nabla }_{\mathbf{y}}p$ are switched when compared with (41a) and [38] has the additional macroscale gradient of the pressure. The fluid stress of [38] is also different

$${\mathsf{T}}^{\left(0\right)}=-p^{\left(0\right)}\mathsf{I},$$

(65)

when compared with (34a), and we see that there are no micro or macroscale gradients of the fluid velocity $\mathbf{v}$.

Equations (42e), (42f), (42i) and (42j) appearing in the problem we derive are analogous to those that form part of the elastic problem in [38]. However, the difference in the scaling of the fluid stress leads to different continuity of stresses across the fluid–solid interfaces (42g) and (42h). That is, in [38], we have the continuity of stresses using (65), (34b) and (34c) instead of that found in (34a).

These small changes in the orders of the terms lead to the formation of only one linear dissipative problem in this work with no decoupling of the phases, rather than being able to separate the problem into a fluid problem and an elastic problem as in [38].

5. Conclusions

In this work, we have presented a system of PDEs governing the effective mechanical behaviour of viscoelastic composites. Our constitutive law takes the form of that of a Kelvin–Voigt material, where we have the effective elasticity tensor applied to the strain plus the viscosity applied to the time derivative of the strain. Our structure comprises multiple solid fibres/subphases and a fluid phase, all of which are in contact with each other. This structure is widely applicable to many real world scenarios including modelling of human bones [3,4,8,9] and skin [5,6,7], as well as in many engineering and manufacturing problems such as biomimetic materials [43] and polymers [44,45].

We begin our analysis by formulating the quasi-static fluid–structure interaction problem describing the behaviour of two families of linear elastic inclusions/fibres interacting with an incompressible Newtonian fluid. We then have made the assumption that both the fluid and the elastic fibres/inclusions/subphases are clearly visible and distinct on the microscale, and also assumed that the macroscale represents the average size of material we are modelling. These two scales are distinct, and their ratio leads to a very small scale separation parameter. We then enforce this distinction in length scales to upscale the non-dimensionalised FSI problem via asymptotic homogenization. The new model derived in this way retains the important properties of the microstructure such as geometry and stiffness in its coefficients, which are determined via solving the presented novel periodic cell problem.

The new model that we derive comprises a balance equation and constitutive law. The viscoelastic constitutive law takes exactly the form that is expected of a viscoelastic material by comprising first the elastic constitutive relationship and second the viscous part of the relation. The addition of the multiple elastic phases in our model influences the elastic part of our constitutive law which reads like that of an elastic composite.

Our new model has an important distinction from formulations of previous viscoelastic models. That is, our model encodes multiple different elastic phases interplaying with each other and the fluid in comparison with some previous viscoelastic formulations where there is only one fluid and solid phase. Accounting for these additional interactions between the multiple phases can be of particular benefit to physical applications. The differences in elastic and mechanical properties between the phases are accounted for by the multiple elasticity tensors. The latter appear in the constitutive law and are accompanied by tensors accounting for the discrepancies in the elastic properties at different points in the microstructure. In this work, we propose the novel cell problem (44a)–(44j) from which the model coefficients are calculated. The cell problem is an extension to the problem found for viscoelastic materials in [17] and comprises the elastic problem associated with poroelastic composites. This cell problem comprises both the fluid and solid equations in one problem and therefore we do not have the decoupling of the different phases as seen in poroelastic-type cell problems.

The addition of linearised inertia and compressibility of fluid to our model would be straightforward. In this case, we would have a system similar to that found in [17]. These additions would also result in corresponding changes to the macroscale model. The effective balance equation that we present in our macroscale model would have the addition of linearised leading order inertia. The compressibility of the fluid would have an influence on the cell problem with the fluid bulk modulus appearing in the Stokes-type problem component.

The current model assumes that both families of elastic phases are anisotropic linear elastic materials. We could, however, extend the model to assume that these phases exhibit hyperelastic behaviour. To do this, we would use a method similar to that found in [31,46,47,48]. By having hyperelastic phases, we dramatically increase the complexity of the numerical simulations that are to be carried out to compute the model coefficients and finally the macroscale model solution. The additional complexity is due to the fact that the two length scales (macro and micro) remain coupled. In the literature, methods to address the remaining coupling of the scales are emerging, see [49,50]. We also note that the present work can be extended for soft hyperelastic electro-active [51] and magneto-active composites [52]. This extension could provide many interesting and exciting applications in modern actuators, soft robotics and biomedicine.

There are many future directions in which this work could be developed. A first step could be to consider the cell problems using the Fourier transform method and then applying the inverse transform to obtaining a formulation in the time domain, see, e.g., ref. [53]. This would lead to cell problems that are computationally feasible to solve that have been parameterised based on real-world data which could be from, for example, biological tissues. Numerical computations performed by finite element simulations have been used in [32] to study the homogenized properties of linear viscoelastic composite materials in three dimensions by means of a semi-analytical approach combined with the asymptotic homogenization method. Finally, our formulation only accounts for Kelvin–Voigt viscoelastic materials at the macroscale. An interesting development of the theory resides in considering more complex constitutive relationships for the individual phases in order to obtain a more general framework for more general viscoelastic behaviours, such as those described in [54].