About RVE size objectivity of multiscale analysis of porous media

Abstract

This study proposes a multi-scale model formulation for saturated porous media, centered on the concept of the Representative Volume Element (RVE). The linkage between scales is established by enforcing the equivalence of the total virtual power per unit volume at the larger scale with its corresponding volume-averaged counterpart at the smaller length scale. By employing the Principle of Multiscale Virtual Power (PMVP) along with appropriate constraints on micro-scale displacements and pore pressures, a robust variational theory is established. The formulation can be implemented using the finite element squared (FE2) strategy through spatial discretization. The theoretical evidence presented in this work reveal a pathological inconsistency in the objectivity of the macro scale response with respect to the RVE size. The primary contribution of this work is to offer an alternative solution to the aforementioned issue, aiming to restore the fundamental concept of RVE. To achieve this, a conveniently fine-scale constitutive approach is proposed, introducing useful adjustments in the micro-scale pore pressure field expansion.

Keywords:
Multiscale model; porous media; RVE

Graphical Abstract

1 INTRODUCTION

The mechanics of porous materials is a highly relevant discipline in various fields of knowledge, such as geomechanics, biomechanics, and materials science, among others. This area has undergone significant development since the pioneering work of Biot (1941) and Biot (1955). Defining constitutive relations at the macro-scale level is crucial for accurately capturing the complex interactions between material constituents and phases. From a modeling perspective, these equations can be: (i) explicitly assumed based on phenomenological considerations, or (ii) implicitly derived using consistent homogenization techniques Dormieux et al. (2006). The latter approach has proven to be highly effective in predicting complex material behaviors. This methodology derives the homogenized constitutive response by averaging quantities from fine-scale fields. The existence of a minimum micro-structural domain, known as the Representative Volume Element (RVE), from which the macro-scale response becomes size-insensitive, is a fundamental cornerstone for these homogenization techniques Hill (1963).

Many RVE-based multiscale strategies are founded on the Hill-Mandel lemma, which ensures equivalence in terms of internal strain energy between the macro and micro scales, Hill (1963), Hill (1965a) and Hill (1965b). Alongside the conceptual improvements, there have also been significant developments in the computational homogenization of these theoretical models, primarily based on the multilevel finite element strategy (FE²) Feyel and Chaboche (2000).

Concerning saturated porous media, Larsson et al. (2010a) introduced a variational homogenization method for analyzing transient uncoupled whereas Su et al. (2011) studied the coupled consolidation phenomena. They demonstrate that the linear multiscale approximation introduces a micro-scale size effect, manifested by the emergence of a second-order term that impacts the objectivity of the macro-scale response, implying that conventional homogenization models cannot be straightforwardly applied to saturated porous media without careful reconsideration. To address this issue, those authors indicated that the second-order effect diminishes as the size of the Representative Volume Element (RVE) approaches infinitesimally small dimensions.

In contrast, Khoei and Hajiabadi (2018) and Khoei and Saeedmonir (2021) put forward computational homogenization models for saturated and multiphase porous media, where the macro-scale response is decomposed into two terms: one stationary and the other dynamic, with the size effect encapsulated in the latter. This effect diminishes as the micro-scale domain approaches infinitesimal dimensions. However, as the fine-scale dimensions increase, these authors observe disparities in the homogenized results when compared to Direct Numerical Simulation (DNS) results.

Klahr et al. (2023) developed a variational RVE-based multiscale formulation for analyzing a saturated porous medium under large strains. In this context, they discovered a dependence on the size of the micro-cell domain, which manifests in the homogenized flow velocity.

All the advancements mentioned in the preceding paragraphs demonstrate a significant interest in understanding the size effect inherent in multiscale modeling of saturated porous media. Currently, this remains a challenging topic and an open line of research. The concept of preserving transient effects, aimed at avoiding limitations on predictive capabilities, underlies all of these previous contributions. The micro-scale size dependence appears to contradict the intuitive concept of Representative Volume Element (RVE) existence, which is a fundamental pillar upon which the entire homogenization theory was established Hill (1963). The primary motivation of this work is to contribute to clarifying ideas concerning this critical theoretical issue and to introduce a new modeling alternative.

In summary, the present work proposes a consistent homogenization model based on the Principle of Multiscale Virtual Power (PMVP) Blanco et al. (2014), along with its corresponding numerical implementation, designed to address the consolidation phenomenon in saturated soils, Lewis and Schrefler (1998), Mroginski et al. (2010), Beneyto et al. (2015). For the mechanical description of the porous medium, we employ the general theory of poromechanics developed by Coussy (2003). To define the internal and external power functionals at both macro and micro-scale levels, two additional pairs of power-conjugate variables are required compared to the classical case of a non-porous medium.

2 PRINCIPLE OF MULTISCALE VIRTUAL POWER

In this section, the coupling of the physical behavior between macro and micro length scales is established by appealing to the so-called Principle of Multiscale Virtual Power (PMVP) proposed by Blanco et al. (2014). This variational statement offers a formal framework to derive all basic ingredients of the multiscale formulation, but firstly some fundamental concepts should be introduced.

2.1 Internal and External power functionals at macro-scale

For the fully saturated case and under the hypothesis of infinitesimal transformations, the internal power ${\mathcal{P}}^{\mathcal{i}\mathcal{n}\mathcal{t}}$ was written by Coussy (2003) as follows

${\mathcal{P}}^{\mathcal{i}\mathcal{n}\mathcal{t}}={\int }_{\mathrm{\Omega }}\left[\mathrm{\sigma }:\dot{\mathrm{\epsilon }}-\nabla \cdot \left(\frac{\mathrm{p}}{{\mathrm{\rho }}^{\mathrm{f}}}\mathrm{w}\right)\right]\mathrm{d}\mathrm{\Omega },$ (1)

where $\dot{\epsilon }={\nabla }^{\mathrm{s}\mathrm{y}\mathrm{m}}\dot{\mathrm{u}}$is the rate of infinitesimal strain tensor of the solid skeleton, $\dot{\mathrm{u}}$is the rate of macro-displacement field, $\mathrm{ó}$ is the Cauchy stress tensor, $\mathrm{p}$the (scalar) pore pressure field of the fluid, ${\mathrm{\rho }}^f$ intrinsic mass density of the fluid, $\mathrm{w}={\mathrm{\rho }}^{\mathrm{f}}\mathrm{}\mathrm{V}$ is the relative mass flow vector of the fluid, $\mathrm{V}$ is the flux or the relative seepage velocity vector between the fluid and solid velocity vectors. Finally, $\mathrm{Ù}$ denote the porous media macro-scale domain.

After some algebraical manipulation and taking into account the mass balance equation ${\dot{\mathrm{m}}}^{\mathrm{f}}+\nabla \cdot \mathrm{w}=0$, being ${\dot{\mathrm{m}}}^{\mathrm{f}}$ the mass content of the fluid phase, assuming infinitesimal transformations, expression (1) can be rewritten as

${\mathcal{P}}^{\mathcal{i}\mathcal{n}\mathcal{t}}={\int }_{\mathrm{\Omega }}\left[\mathrm{\sigma }:\dot{\mathrm{\epsilon }}+\dot{\mathrm{\chi }}\mathrm{p}-\mathrm{V}\cdot \mathrm{\phi }\right]\mathrm{d}\mathrm{\Omega },$ (2)

where the following nomenclature is adopted hereafter: $\dot{\mathrm{\chi }}=\frac{\dot{{\mathrm{m}}^{\mathrm{f}}}}{{\mathrm{\rho }}^{\mathrm{f}}}$ and $\mathrm{ö}=\nabla \mathrm{p}$.

It can be seen that the second term in the r.h.s. of (2) represents the work rate done, in the pore space, by the fluid phase. The third duality product in (2) takes into account the viscous dissipation effect due to the relative motion of the fluid phase concerning the solid skeleton. Therefore, two additional pairs of power-conjugate variables, $\{\dot{\mathrm{\chi }};\mathrm{p}\}$ and $\{\mathrm{V};\mathrm{\phi }\}$, are necessary to introduce with respect to the classical solid mechanical problem.

The external power, ${\mathcal{P}}^{\mathcal{e}\mathcal{x}\mathcal{t}}$, exerted by the set of external agencies $\{\mathrm{f};\mathrm{t};\mathrm{q}\}$ has also been presented by Coussy (2003)

${\mathcal{P}}^{\mathcal{e}\mathcal{x}\mathcal{t}}={\int }_{\mathrm{\Omega }}\mathrm{f}\cdot \dot{\mathrm{u}}\mathrm{d}\mathrm{\Omega }+{\int }_{\mathrm{\Gamma }}\left[\mathrm{t}\cdot \dot{\mathrm{u}}-\mathrm{p}\frac{\mathrm{q}}{{\mathrm{\rho }}^{\mathrm{f}}}\right]\mathrm{d}\mathrm{\Omega },$ (3)

where $\mathrm{f}$ is the body force of the porous media, while $\mathrm{t}$ and $\mathrm{q}$ represent the systems of total external traction and normal fluid flux, respectively.

The well-known concepts of internal and external virtual powers can be easily established from the above expressions, considering the corresponding admissible virtual actions of their primal quantities, i.e., $\{\mathrm{ä}\dot{\mathrm{u}};\mathrm{\delta }\dot{\mathrm{\epsilon }};\mathrm{\delta }\mathrm{p};\mathrm{\delta }\mathrm{\phi }\}$. Then, the total virtual power, $\mathrm{ä}{\mathcal{P}}^{\mathcal{t}\mathcal{o}\mathcal{t}}$, can be expressed as

$\mathrm{ä}{\mathcal{P}}^{\mathcal{t}\mathcal{o}\mathcal{t}}=\mathrm{\delta }{\mathcal{P}}^{\mathcal{i}\mathcal{n}\mathcal{t}}-\mathrm{\delta }{\mathcal{P}}^{\mathcal{e}\mathcal{x}\mathcal{t}}={\int }_{\mathrm{\Omega }}\left(\mathrm{\sigma }:\mathrm{\delta }\dot{\mathrm{\epsilon }}+\dot{\mathrm{\chi }}\mathrm{}\mathrm{\delta }\mathrm{p}-\mathrm{V}\cdot \mathrm{\delta }\mathrm{\phi }\mathrm{}-\mathrm{f}\cdot \mathrm{\delta }\dot{\mathrm{u}}\right)\mathrm{d}\mathrm{\Omega }-{\int }_{{\mathrm{\Gamma }}_{\text{N}}^{\mathrm{u}}}\mathrm{t}\cdot \mathrm{\delta }\dot{\mathrm{u}}\mathrm{d}\mathrm{\Gamma }+{\int }_{{\mathrm{\Gamma }}_{\text{N}}^{\mathrm{p}}}\mathrm{\delta }\mathrm{p}\frac{\mathrm{q}}{{\mathrm{\rho }}^{\mathrm{f}}}\mathrm{d}\mathrm{\Gamma }.$ (4)

2.2 Weak form of the macro-scale balance equations

The balance equations in a weak form postulate that the external virtual power must be equal to the internal virtual power, for all admissible variations of primal descriptors. Thus, (4) can be rewritten as

${\int }_{\mathrm{\Omega }}\left(\mathrm{\sigma }:\mathrm{\delta }\dot{\mathrm{\epsilon }}+\dot{\mathrm{\chi }}\mathrm{\delta }\mathrm{p}-\mathrm{V}\cdot \mathrm{\delta }\mathrm{\phi }\right)\mathrm{d}\mathrm{\Omega }={\int }_{\mathrm{\Omega }}\mathrm{f}\cdot \mathrm{\delta }\dot{\mathrm{u}}\mathrm{d}\mathrm{\Omega }+{\int }_{{\mathrm{\Gamma }}_{\text{N}}^{\mathrm{u}}}\mathrm{t}\cdot \mathrm{\delta }\dot{\mathrm{u}}\mathrm{d}\mathrm{\Gamma }-{\int }_{{\mathrm{\Gamma }}_{\text{N}}^{\mathrm{p}}}\frac{\mathrm{q}}{{\mathrm{\rho }}^{\mathrm{f}}}\mathrm{\delta }\mathrm{p}\mathrm{d}\mathrm{\Gamma },\mathrm{}\forall \mathrm{}\mathrm{\delta }\dot{\mathrm{u}}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{}\mathrm{\delta }\mathrm{p}\mathrm{}\mathrm{a}\mathrm{d}\mathrm{m}\mathrm{i}\mathrm{s}\mathrm{s}\mathrm{i}\mathrm{b}\mathrm{l}\mathrm{e}.$ (5)

The underlying admissibility requirements invoked in (5), for virtual variations $\mathrm{ä}\dot{\mathrm{u}}$ and $\mathrm{ä}\mathrm{p}$, take into account proper regularity demands such that all the integral terms can be formally evaluated as well as homogeneous prescribed values for both continuous fields on ${\mathrm{\Gamma }}_{\text{N}}^{\mathrm{u}}$ and ${\mathrm{\Gamma }}_{\text{N}}^{\mathrm{p}}$ (i.e., where Dirichlet boundary conditions should be specified for $\dot{\mathrm{u}}$ and $\mathrm{p}$), respectively. Since $\mathrm{ä}\dot{\mathrm{u}}$ and $\mathrm{ä}\mathrm{p}$ are independent from each other, the variational form of equilibrium is finally described as a system of two coupled scalar equations, Lewis and Schrefler (1998) Di Rado et al. (2009).

$\begin{array}{l}\mathrm{G}\equiv {\int }_{\mathrm{\Omega }}\mathrm{\sigma }:\mathrm{\delta }\dot{\mathrm{\epsilon }}\mathrm{d}\mathrm{\Omega }\mathrm{}-{\int }_{\mathrm{\Omega }}\mathrm{f}\cdot \mathrm{\delta }\dot{\mathrm{u}}\mathrm{d}\mathrm{\Omega }\mathrm{}-{\int }_{{\mathrm{\Gamma }}_{\text{N}}^{\mathrm{u}}}\mathrm{t}\cdot \mathrm{\delta }\dot{\mathrm{u}}\mathrm{d}\mathrm{\Gamma },\mathrm{}\forall \mathrm{}\mathrm{\delta }\dot{\mathrm{u}}\mathrm{}\mathrm{a}\mathrm{d}\mathrm{m}\mathrm{i}\mathrm{s}\mathrm{s}\mathrm{i}\mathrm{b}\mathrm{l}\mathrm{e},\\ \mathrm{H}\equiv {\int }_{\mathrm{\Omega }}\left(\dot{\mathrm{\chi }}\mathrm{\delta }\mathrm{p}-\mathrm{V}\cdot \mathrm{\delta }\mathrm{\phi }\right)\mathrm{d}\mathrm{\Omega }\mathrm{}+\mathrm{}{\int }_{{\mathrm{\Gamma }}_{\text{N}}^{\mathrm{p}}}\frac{\mathrm{q}}{{\mathrm{\rho }}^{\mathrm{f}}}\mathrm{\delta }\mathrm{p}\mathrm{d}\mathrm{\Gamma },\mathrm{}\forall \mathrm{}\mathrm{\delta }\mathrm{p}\mathrm{}\mathrm{a}\mathrm{d}\mathrm{m}\mathrm{i}\mathrm{s}\mathrm{s}\mathrm{i}\mathrm{b}\mathrm{l}\mathrm{e}.\end{array}$ (6)

2.3 Principle of Multiscale Virtual Power

Before discussing the PMVP, it is required to postulate an adequate mechanism for transferring information and define the concept of admissibility that makes compatible the primitive descriptors (displacements, pore pressures and their corresponding gradients) between the two involved scales.

Henceforth, any object related to the micro-scale will be endowed with the subscript $\mathrm{ì}$, preserving the same physical meaning. The micro-scale domain is denoted ${\mathrm{\Omega }}_{\mathrm{\mu }}$ (also called Micro-Cell or simply MC), with volume $\left|{\mathrm{\Omega }}_{\mathrm{\mu }}\right|$ boundary ${\mathrm{\Gamma }}_{\mathrm{\mu }}$, while ${\mathrm{n}}_{\mathrm{\mu }}$ is the unit (outward) vector normal to ${\mathrm{\Gamma }}_{\mathrm{\mu }}$. In addition, the vector $\mathrm{y}$ is used to describe material points in a Cartesian system at the sub-scale. Without loss of generality, the origin of coordinates is located at the geometric center of ${\mathrm{\Omega }}_{\mathrm{\mu }}$, implying that

${\int }_{{\mathrm{\Omega }}_{\mathrm{\mu }}}\mathrm{y}\mathrm{}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{\mu }}=0$.(7)

2.4 Primal descriptors at the micro-scale level

As in the macro-scale problem, the Micro-Cell is constituted by a continuum saturated porous medium. Thus, the same primal descriptors of the poromechanics theory of the macro-scale are adopted. A fundamental assumption in the present modelling context is that the micro-scale displacement, ${\mathrm{u}}_{\mathrm{\mu }}$, can be expanded in terms of macro-scale quantities (${\mathrm{u}}_{\mathrm{\mu }}$ and $\mathrm{å}$) as follows (de Souza Neto and Feijóo (2006), Sánchez et al. (2013) and de Souza Neto et al. (2015))

${\mathrm{u}}_{\mathrm{\mu }}\left(\mathrm{y},\mathrm{t}\right)=\mathrm{u}\left(\mathrm{x},\mathrm{t}\right)+\mathrm{\epsilon }\left(\mathrm{x},\mathrm{t}\right)\cdot \mathrm{y}+{\tilde{\mathrm{u}}}_{\mathrm{\mu }}\left(\mathrm{y},\mathrm{t}\right)$,(8)

being ${\tilde{\mathrm{u}}}_{\mathrm{\mu }}\left(\mathrm{y},\mathrm{t}\right)$ the fluctuation displacement field of the MC, therefore, the micro-scale infinitesimal strain tensor yields

${\mathrm{\epsilon }}_{\mathrm{\mu }}\left(\mathrm{y},\mathrm{t}\right)={\nabla }_{\mathrm{y}}^{\mathrm{s}\mathrm{y}\mathrm{m}}{\mathrm{u}}_{\mathrm{\mu }}\left(\mathrm{y},\mathrm{t}\right)=\mathrm{\epsilon }\left(\mathrm{x},\mathrm{t}\right)+{\nabla }_{\mathrm{y}}^{\mathrm{s}\mathrm{y}\mathrm{m}}{\tilde{\mathrm{u}}}_{\mathrm{\mu }}\left(\mathrm{y},\mathrm{t}\right)=\mathrm{\epsilon }\left(\mathrm{x},\mathrm{t}\right)+{\tilde{\mathrm{\epsilon }}}_{\mathrm{\mu }}\left(\mathrm{y},\mathrm{t}\right)$,(9)

where operator ${\nabla }_{\mathrm{y}}^{\mathrm{s}\mathrm{y}\mathrm{m}}\left(\mathrm{∎}\right)$ represents the symmetric gradient concerning $\mathrm{y}$-coordinates. The strain field in (9) is composed as the addition of the macro-scale strain, $\mathrm{å}\left(\mathrm{x},\mathrm{t}\right)$, which is assumed to be uniformly distributed over the entire domain ${\mathrm{\Omega }}_{\mathrm{\mu }}$, and the micro-scale strain fluctuation ${\tilde{\mathrm{\epsilon }}}_{\mathrm{\mu }}\left(\mathrm{y},\mathrm{t}\right)$.

On the other hand, the micro-pore pressure, and its spatial gradient, can also be decomposed in an additive manner analogous to the displacement (Khoei and Hajiabadi (2018), Rivarola et al. (2019) and Anonis et al. (2024a))

${\mathrm{p}}_{\mathrm{\mu }}\left(\mathrm{y},\mathrm{t}\right)=\mathrm{p}\left(\mathrm{x},\mathrm{t}\right)+\mathrm{\phi }\left(\mathrm{x},\mathrm{t}\right)\cdot \mathrm{y}+{\tilde{\mathrm{p}}}_{\mathrm{\mu }}\left(\mathrm{y},\mathrm{t}\right)$,(10)

${\mathrm{\phi }}_{\mathrm{\mu }}\left(\mathrm{y},\mathrm{t}\right)={\nabla }_{\mathrm{y}}{\mathrm{p}}_{\mathrm{\mu }}\left(\mathrm{y},\mathrm{t}\right)=\mathrm{\phi }\left(\mathrm{x},\mathrm{t}\right)+{\nabla }_{\mathrm{y}}{\tilde{\mathrm{p}}}_{\mathrm{\mu }}\left(\mathrm{y},\mathrm{t}\right)=\mathrm{\phi }\left(\mathrm{x},\mathrm{t}\right)+{\tilde{\mathrm{\phi }}}_{\mathrm{\mu }}\left(\mathrm{y},\mathrm{t}\right)$,(11)

being ${\tilde{\mathrm{p}}}_{\mathrm{\mu }}\left(\mathrm{y},\mathrm{t}\right)$ the pore pressure fluctuation field on the MC and ${\tilde{\mathrm{\phi }}}_{\mathrm{\mu }}\left(\mathrm{y},\mathrm{t}\right)$ its corresponding gradient vector which is obtained from the ${\nabla }_{\mathrm{y}}\left(\mathrm{∎}\right)$ operator, related to the micro-scale coordinate system.

2.5 Formulation of the PMVP

The Principle of Multiscale Virtual Power states that the total virtual power per unit volume, at a point $x$ of the macro-scale, must be equal to the volumetric average of the total micro-scale virtual power (per unit volume) at the corresponding MC, for all admissible virtual actions (Blanco et al. (2014)).

Thus, recalling the definition of macro-scale total virtual power per unit volume (given by the first integrand term of the r.h.s. in (4) and assuming the same mathematical structure for its micro-scale counterpart, the PMVP gives us the following variational sentence

$\mathrm{ó}:\mathrm{ä}\dot{\mathrm{\epsilon }}+\dot{\mathrm{\chi }}\mathrm{\delta }\mathrm{p}-\mathrm{V}\cdot \mathrm{\delta }\mathrm{\phi }-\mathrm{f}\cdot \mathrm{\delta }\dot{\mathrm{u}}=\frac{1}{\left|{\mathrm{\Omega }}_{\mathrm{\mu }}\right|}\mathrm{}{\int }_{{\mathrm{\Omega }}_{\mathrm{\mu }}}\left({\mathrm{\sigma }}_{\mathrm{\mu }}:\mathrm{\delta }{\dot{\mathrm{\epsilon }}}_{\mathrm{\mu }}+{\dot{\mathrm{\chi }}}_{\mathrm{\mu }}\mathrm{\delta }{\mathrm{p}}_{\mathrm{\mu }}-{\mathrm{V}}_{\mathrm{\mu }}\cdot \mathrm{\delta }{\mathrm{\phi }}_{\mathrm{\mu }}-{\mathrm{f}}_{\mathrm{\mu }}\cdot \mathrm{\delta }{\dot{\mathrm{u}}}_{\mathrm{\mu }}\right)\mathrm{d}{\mathrm{\Omega }}_{\mathrm{\mu }},\forall \mathrm{\delta }\dot{\mathrm{u}},\mathrm{}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{}\mathrm{\delta }\mathrm{p}\mathrm{}\mathrm{a}\mathrm{d}\mathrm{m}\mathrm{i}\mathrm{s}\mathrm{s}\mathrm{i}\mathrm{b}\mathrm{l}\mathrm{e}\mathrm{}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{}\forall \mathrm{\delta }{\dot{\mathrm{u}}}_{\mathrm{\mu }}.$ (12)

Expression (12) can be viewed as a particular instance of the PMVP for the case of saturated porous media, at both scales of analysis.

2.6 Homogenized variables and variational forms of balance at the micro-scale

The variational identity (12) provides the natural way to obtain the homogenization formulae for the macro-scale stress-like entities$\{\mathrm{ó};\dot{\mathrm{\chi }};\mathrm{V}\}$ and body force $\mathrm{f}$, as well as the variational equilibrium equation at micro-scale. To attain this goal, descriptors in Eqs. (8)-(11) must be replaced in expression (12) for virtual actions in the MC domain, and then, the consequences are shown below.

a) Homogenized stress tensor:

$\mathrm{ó}=\frac{1}{\left|{\mathrm{\Omega }}_{\mathrm{\mu }}\right|}{\int }_{{\mathrm{\Omega }}_{\mathrm{\mu }}}\left({\mathrm{\sigma }}_{\mathrm{\mu }}-{\mathrm{f}}_{\mathrm{\mu }}\otimes \mathrm{y}\right)\mathrm{d}{\mathrm{\Omega }}_{\mathrm{\mu }},\forall \mathrm{t}$.(13)

Obtained from (12), taking $\mathrm{ä}{\dot{\tilde{\mathrm{u}}}}_{\mathrm{\mu }}=0,\mathrm{}\mathrm{\delta }\mathrm{p}=0,\mathrm{}\mathrm{\delta }{\tilde{\mathrm{p}}}_{\mathrm{\mu }}=0,\mathrm{}\mathrm{\delta }\mathrm{\phi }=0,\mathrm{\delta }\dot{\mathrm{u}}=0$ and allowing arbitrary variations of $\mathrm{ä}\dot{\mathrm{\epsilon }}$.

b) Homogenized mass content rate of fluid (per unit fluid density):

$\dot{\mathrm{\chi }}=\frac{{\dot{\mathrm{m}}}^{\mathrm{f}}}{{\mathrm{\rho }}^{\mathrm{f}}}=\frac{1}{\left|{\mathrm{\Omega }}_{\mathrm{\mu }}\right|}{\int }_{{\mathrm{\Omega }}_{\mathrm{\mu }}}{\dot{\mathrm{\chi }}}_{\mathrm{\mu }}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{\mu }},\forall \mathrm{t}\mathrm{}.$ (14)

Obtained from (12), taking $\mathrm{ä}\dot{\mathrm{\epsilon }}=0,\mathrm{\delta }{\dot{\tilde{\mathrm{u}}}}_{\mathrm{\mu }}=0,\mathrm{}\mathrm{\delta }{\tilde{\mathrm{p}}}_{\mathrm{\mu }}=0,\mathrm{\delta }\dot{\mathrm{u}}=0,\mathrm{}\mathrm{\delta }\mathrm{\phi }=0$ and allowing arbitrary variations of $\mathrm{ä}\mathrm{p}$.

c) Homogenized flux velocity vector:

$\mathrm{V}=\frac{1}{\left|{\mathrm{\Omega }}_{\mathrm{\mu }}\right|}{\int }_{{\mathrm{\Omega }}_{\mathrm{\mu }}}\left({\mathrm{V}}_{\mathrm{\mu }}-{\dot{\mathrm{\chi }}}_{\mathrm{\mu }}\mathrm{y}\right)\mathrm{d}{\mathrm{\Omega }}_{\mathrm{\mu }},\forall \mathrm{t}\mathrm{}.$ (15)

Which is achieved from (12), adopting $\mathrm{ä}\dot{\mathrm{\epsilon }}=0,\mathrm{}\mathrm{\delta }{\dot{\tilde{\mathrm{u}}}}_{\mathrm{\mu }}=0,\mathrm{}\mathrm{\delta }\mathrm{p}=0,\mathrm{}\mathrm{\delta }{\tilde{\mathrm{p}}}_{\mathrm{\mu }}=0,\mathrm{\delta }\dot{\mathrm{u}}=0$ with arbitrary variations of $\mathrm{ä}\mathrm{ö}$.

From (15) it is possible to decompose the homogenized flux velocity vector into a stationary part (${\mathrm{V}}_{\mathrm{s}\mathrm{t}\mathrm{a}})$and a dynamic or transient part (${\mathrm{V}}_{\mathrm{t}\mathrm{r}\mathrm{a}}$) as Janicke et al. (2020), Khoei and Hajiabadi (2018), Khoei and Saeedmonir (2021), Saeedmonir and Khoei (2022) (this decomposition is useful to discuss the MC size effect problem):

$\mathrm{V}={\mathrm{V}}_{\mathrm{s}\mathrm{t}\mathrm{a}}+{\mathrm{V}}_{\mathrm{t}\mathrm{r}\mathrm{a}}=\frac{1}{\left|{\mathrm{\Omega }}_{\mathrm{\mu }}\right|}{\int }_{{\mathrm{\Omega }}_{\mathrm{\mu }}}{\mathrm{V}}_{\mathrm{\mu }}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{\mu }}-\frac{1}{\left|{\mathrm{\Omega }}_{\mathrm{\mu }}\right|}{\int }_{{\mathrm{\Omega }}_{\mathrm{\mu }}}{\dot{\mathrm{\chi }}}_{\mathrm{\mu }}\mathrm{y}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{\mu }},\forall \mathrm{t}\mathrm{}.$ (16)

d) Homogenized body force field:

$\mathrm{f}=\frac{1}{\left|{\mathrm{\Omega }}_{\mathrm{\mu }}\right|}{\int }_{{\mathrm{\Omega }}_{\mathrm{\mu }}}{\mathrm{f}}_{\mathrm{\mu }}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{\mu }}=\frac{\mathrm{g}}{\left|{\mathrm{\Omega }}_{\mathrm{\mu }}\right|}{\int }_{{\mathrm{\Omega }}_{\mathrm{\mu }}}{\mathrm{\rho }}_{\mathrm{\mu }}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{\mu }},\forall \mathrm{t}\mathrm{}.$ (17)

Obtained from (12), taking $\mathrm{ä}\dot{\mathrm{\epsilon }}=0,\mathrm{}\mathrm{\delta }{\dot{\tilde{\mathrm{u}}}}_{\mathrm{\mu }}=0,\mathrm{}\mathrm{\delta }\mathrm{p}=0,\mathrm{}\mathrm{\delta }{\tilde{\mathrm{p}}}_{\mathrm{\mu }}=0,\mathrm{\delta }\mathrm{\phi }=0$ and arbitrary variations of$\delta \dot{\mathrm{u}}=0$, being $\mathrm{g}$ the gravity acceleration field.

e) Variational forms of balance in the MC (Integral equation of momentum balance):

${\int }_{{\mathrm{\Omega }}_{\mathrm{\mu }}}\left({\mathrm{\sigma }}_{\mathrm{\mu }}:\mathrm{\delta }\dot{\tilde{\mathrm{\epsilon }}}-{\mathrm{f}}_{\mathrm{\mu }}\cdot \mathrm{}\mathrm{\delta }{\dot{\tilde{\mathrm{u}}}}_{\mathrm{\mu }}\right)\mathrm{d}{\mathrm{\Omega }}_{\mathrm{\mu }}={\int }_{{\mathrm{\Omega }}_{\mathrm{\mu }}}\left({\mathrm{\sigma }}_{\mathrm{\mu }}:{\nabla }_{\mathrm{y}}^{\mathrm{s}\mathrm{y}\mathrm{m}}\mathrm{\delta }{\dot{\tilde{\mathrm{u}}}}_{\mathrm{\mu }}-{\mathrm{f}}_{\mathrm{\mu }}\cdot \mathrm{\delta }{\dot{\tilde{\mathrm{u}}}}_{\mathrm{\mu }}\right)\mathrm{d}{\mathrm{\Omega }}_{\mathrm{\mu }}=0,\mathrm{}\mathrm{}\forall \mathrm{}\mathrm{\delta }{\dot{\tilde{\mathrm{u}}}}_{\mathrm{\mu }}\in \mathrm{}{\tilde{\mathcal{U}}}_{\mathrm{\mu }},\forall \mathrm{t}\mathrm{}.$ (18)

Deduced from (12) by allowing for admissible variations of $\mathrm{ä}{\dot{\tilde{\mathrm{u}}}}_{\mathrm{\mu }}$ with $\mathrm{ä}\dot{\mathrm{\epsilon }}=0,\mathrm{}\mathrm{\delta }\mathrm{p}=0,\mathrm{}\mathrm{\delta }{\tilde{\mathrm{p}}}_{\mathrm{\mu }}=0,\mathrm{\delta }\mathrm{\phi }=0\mathrm{}$and $\mathrm{ä}\dot{\mathrm{u}}=0$.

f) Variational forms of balance in the MC (Integral mass balance equation):

${\int }_{{\mathrm{\Omega }}_{\mathrm{\mu }}}\left({\dot{\mathrm{\chi }}}_{\mathrm{\mu }}\mathrm{}\mathrm{\delta }{\tilde{\mathrm{p}}}_{\mathrm{\mu }}-{\mathrm{V}}_{\mathrm{\mu }}\cdot \mathrm{\delta }{\tilde{\mathrm{\phi }}}_{\mathrm{\mu }}\right)\mathrm{}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{\mu }}={\int }_{{\mathrm{\Omega }}_{\mathrm{\mu }}}\left({\dot{\mathrm{\chi }}}_{\mathrm{\mu }}\mathrm{}\mathrm{\delta }{\tilde{\mathrm{p}}}_{\mathrm{\mu }}-{\mathrm{V}}_{\mathrm{\mu }}\cdot {\nabla }_{\mathrm{y}}\mathrm{}\mathrm{\delta }{\tilde{\mathrm{p}}}_{\mathrm{\mu }}\right)\mathrm{}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{\mu }}=0,\mathrm{}\mathrm{}\forall \mathrm{}\mathrm{}\mathrm{\delta }{\tilde{\mathrm{p}}}_{\mathrm{\mu }}\in \mathrm{}{\tilde{\mathcal{P}}}_{\mathrm{\mu }},\forall \mathrm{t}\mathrm{}.$ (19)

Extracterd from (12), by allowing for admissible variations of $\mathrm{ä}{\tilde{\mathrm{p}}}_{\mathrm{\mu }}$ with $\mathrm{ä}\dot{\mathrm{\epsilon }}=0,\mathrm{}\mathrm{\delta }{\dot{\tilde{\mathrm{u}}}}_{\mathrm{\mu }}=0,\mathrm{}\mathrm{\delta }\mathrm{p}=0,\mathrm{\delta }\dot{\mathrm{u}}=0$ and $\mathrm{ä}{\tilde{\mathrm{\phi }}}_{\mathrm{\mu }}$.

3 CONSTITUTIVE EQUATIONS AT THE MICRO-SCALE

The only remaining ingredient in the proposed multiscale formulation is the specification of the material behavior in the smaller length scale. Each constituent of the micro-scale domain is a two-phase saturated porous medium; therefore, constitutive relations for the seepage velocity vector ${\mathrm{V}}_{\mathrm{\mu }}$, as well as for the mechanical stress-like quantities $\{{\mathrm{\sigma }}_{\mathrm{\mu }};{\dot{\mathrm{\chi }}}_{\mathrm{\mu }}\}$, are required. The way in which constitutive input arguments are considered to evaluate such material responses plays an important role for the purposes of this work. The above-mentioned micro-scale size dependence could be one of them, which is objectionable in the realm of RVE-based homogenization models because induces a lack of objectivity in the macro-scale response. This topic justifies the detailed treatment of the constitutive arguments presented in the following sections.

3.1 Constitutive laws based on Full Order Expansions (FOE) for primal descriptors

This strategy probably represents the most natural and straightforward choice to describe the micro-scale material behaviors. Indeed, it has been adopted in some previous contributions Larsson et al. (2010a), Su et al. (2011), Janicke et al. (2022) Khoei and Hajiabadi (2018), Khoei and Saeedmonir (2021), Saeedmonir and Khoei (2022), Wu et al. (2022). In this approach, the Full Order Expanded (FOE) version of primal descriptors, see definitions (8)-(11), are part of the list of arguments in the constitutive functionals.

With respect to the seepage phenomenon at the micro-scale level, the generalized Darcy's law is assumed to characterize the mean fluid velocity in the saturated porous medium Mroginski et al. (2010), Lewis and Schrefler (1998).

${\mathrm{V}}_{\mathrm{\mu }}={\widehat{\mathcal{V}}}_{\mathrm{\mu }}^{\text{FOE}}\left({\mathrm{\phi }}_{\mathrm{\mu }}^{\text{F0E}}\right)=-{\mathrm{k}}_{\mathrm{\mu }}\cdot \left[\left(\mathrm{\phi }+{\tilde{\mathrm{\phi }}}_{\mathrm{\mu }}\right)+{\mathrm{\rho }}_{\mathrm{\mu }}^{\mathrm{f}}\mathrm{g}\right]={\widehat{\mathcal{V}}}_{\mathrm{\mu }}^{\text{FOE}}\left(\mathrm{\phi },{\tilde{\mathrm{\phi }}}_{\mathrm{\mu }}\right)$,(20)

where ${\mathrm{k}}_{\mathrm{\mu }}$ is the symmetric second order permeability tensor and the hat-symbol, $\widehat{∎}$, denotes a generic constitutive functional. In case of isotropic saturated materials, it is ${\mathrm{k}}_{\mathrm{\mu }}={\mathrm{k}}_{\mathrm{\mu }}\mathrm{I};{\mathrm{k}}_{\mathrm{\mu }}={\mathrm{\kappa }}_{\mathrm{\mu }}/\left({\mathrm{\rho }}_{\mathrm{\mu }}^{\mathrm{f}}\left|\mathrm{g}\right|\right)$ being the hydraulic permeability which is a function of the hydraulic conductivity ${\mathrm{\kappa }}_{\mathrm{\mu }}$ and the specific weight of the fluid ${\mathrm{\rho }}_{\mathrm{\mu }}^{\mathrm{f}}\left|\mathrm{g}\right|$, $\left|\mathrm{g}\right|$ is the modulus of the acceleration of gravity and $\mathrm{I}$ denotes the second order identity tensor.

The same mathematical structure that describes the mechanical behavior of saturated porous medium due to Coussy (2003) is assumed valid for each constituent of the micro-scale domain. Thus, we have for $\{{\dot{\mathrm{\sigma }}}_{\mathrm{\mu }};{\dot{\chi }}_{\mu }\}$ (rate format is used)

${\dot{\mathrm{\sigma }}}_{\mathrm{\mu }}={\widehat{\dot{\mathrm{\sigma }}}}_{\mathrm{\mu }}^{\text{FOE}}\left({\dot{\mathrm{\epsilon }}}_{\mathrm{\mu }}^{\text{FOE}},{\dot{\mathrm{p}}}_{\mathrm{\mu }}^{\text{FOE}}\right)={\mathrm{C}}_{\mathrm{\mu }}:\left(\dot{\mathrm{\epsilon }}+{\dot{\tilde{\mathrm{\epsilon }}}}_{\mathrm{\mu }}\right)-{\mathrm{b}}_{\mathrm{\mu }}\left(\dot{\mathrm{p}}+\dot{\mathrm{\phi }}\cdot \mathrm{y}+{\dot{\tilde{\mathrm{p}}}}_{\mathrm{\mu }}\right)={\widehat{\dot{\mathrm{\sigma }}}}_{\mathrm{\mu }}^{\text{FOE}}\left(\dot{\mathrm{\epsilon }},{\dot{\tilde{\mathrm{\epsilon }}}}_{\mathrm{\mu }},\dot{\mathrm{p}},\dot{\mathrm{\phi }},{\dot{\tilde{\mathrm{p}}}}_{\mathrm{\mu }}\right),$ (21)

${\dot{\mathrm{\chi }}}_{\mathrm{\mu }}={\widehat{\dot{\mathrm{\chi }}}}_{\mathrm{\mu }}^{\text{FOE}}\left({\dot{\mathrm{\epsilon }}}_{\mathrm{\mu }}^{\text{FOE}},{\dot{\mathrm{p}}}_{\mathrm{\mu }}^{\text{FOE}}\right)={\mathrm{b}}_{\mathrm{\mu }}:\left(\dot{\mathrm{\epsilon }}+{\dot{\tilde{\mathrm{\epsilon }}}}_{\mathrm{\mu }}\right)+\frac{1}{{\mathrm{M}}_{\mathrm{\mu }}}\left(\dot{\mathrm{p}}+\dot{\mathrm{\phi }}\cdot \mathrm{y}+{\dot{\tilde{\mathrm{p}}}}_{\mathrm{\mu }}\right)={\widehat{\dot{\mathrm{\chi }}}}_{\mathrm{\mu }}^{\text{FOE}}\left(\dot{\mathrm{\epsilon }},{\dot{\tilde{\mathrm{\epsilon }}}}_{\mathrm{\mu }},\dot{\mathrm{p}},\dot{\mathrm{\phi }},{\dot{\tilde{\mathrm{p}}}}_{\mathrm{\mu }}\right).$ (22)

The micro-hydromechanical response of the porous medium is mainly defined by the elastic stiffness tensor of the skeleton ${\mathrm{C}}_{\mathrm{\mu }}$, the Biot tensor ${\mathrm{b}}_{\mathrm{\mu }}\mathrm{}$and the coefficient ${\mathrm{M}}_{\mathrm{\mu }}^{-1}$. For the particular case of isotropic constituents, it is ${\mathrm{b}}_{\mathrm{\mu }}={\mathrm{b}}_{\mathrm{\mu }}\mathrm{I}$, where the Biot coefficient ${\mathrm{b}}_{\mathrm{\mu }}=1-{\mathrm{K}}_{\mathrm{\mu }}/{\mathrm{K}}_{\mathrm{\mu }}^{\mathrm{s}}\mathrm{}$is given by the relationship between the bulk modulus of the skeleton ${\mathrm{K}}_{\mathrm{\mu }}$ and the volumetric modulus of the grain material ${\mathrm{K}}_{\mathrm{\mu }}^{\mathrm{s}}$. Finally, we have ${\mathrm{M}}_{\mathrm{\mu }}^{-1}={\mathrm{K}}_{\mathrm{\mu }}^{\mathrm{s}}/\left({\mathrm{b}}_{\mathrm{\mu }}-{\mathrm{n}}_{\mathrm{\mu }}\right)+{\mathrm{K}}_{\mathrm{\mu }}^{\mathrm{f}}/{\mathrm{n}}_{\mathrm{\mu }}$, where ${\mathrm{K}}_{\mathrm{\mu }}^{\mathrm{f}}\mathrm{}$is the bulk modulus of fluid Biot (1941), Biot (1955) and Lewis and Schrefler (1998).

In saturated porous media significant mechanical changes undergone by the soil skeleton are attributed to the concept of the effective stress field, here denoted as ${\mathrm{\sigma }}_{\mathrm{\mu }}^{\mathrm{\text{'}}}$. According to the poromechanics theory, ${\mathrm{\sigma }}_{\mathrm{\mu }}^{\mathrm{\text{'}}}$ depends on the strain tensor of the soil grains (${\mathrm{\epsilon }}_{\mathrm{\mu }}$), then in agreement with the previous constitutive law (21), it can be expressed

${\dot{\mathrm{\sigma }}}_{\mathrm{\mu }}^{\mathrm{\text{'}}}={\widehat{\dot{\mathrm{\sigma }\mathrm{\text{'}}}}}_{\mathrm{\mu }}^{\text{FOE}}\left({\dot{\mathrm{\epsilon }}}_{\mathrm{\mu }}^{\text{FOE}}\right)={\mathrm{C}}_{\mathrm{\mu }}:\left(\dot{\mathrm{\epsilon }}+{\dot{\tilde{\mathrm{\epsilon }}}}_{\mathrm{\mu }}\right)={\widehat{\dot{\mathrm{\sigma }\mathrm{\text{'}}}}}_{\mathrm{\mu }}^{\text{FOE}}\left(\dot{\mathrm{\epsilon }},{\dot{\tilde{\mathrm{\epsilon }}}}_{\mathrm{\mu }}\right)$ (23)

3.2 About the micro-scale size dependence on the macro-scale response

The transient component, ${\mathrm{V}}_{\mathrm{t}\mathrm{r}\mathrm{a}}$, of the homogenized velocity vector (see second integral in the r.h.s. of (16), has already been identified as responsible for introducing a size effect in the multiscale modelling of porous saturated solids (Larsson at al. (2010a), Su et al. (2011, Khoei and Hajiabadi (2018), Janicke et al. (2020), Saeedmonir and Khoei (2022) and Anonis et al. (2024a)). Some possible solutions have been proposed in the current literature. One is to extend the applied proposal in the context of heat flow (Larsson et al. (2010b), Özdemir et al. (2008)) and simply ignore the transient effect of the sub-scale by neglecting the ${\mathrm{V}}_{\mathrm{t}\mathrm{r}\mathrm{a}}$-contribution so that the micro-scale problem can be considered as “quasi-stationary” at all times. Another one is to adopt an infinitely small micro-scale domain size, which allows the dynamic term to be negligible and thus avoids the problem of fine-scale size dependence (Larsson at al. (2010a), Su et al. (2011, Khoei and Hajiabadi (2018), Janicke et al. (2020)). The possibility of choosing specific material parameters at the MC level in order to make the dynamic fluctuating response negligible compared with the stationary one was also mentioned in Janicke et al. (2020). In view of the previous background and the importance assigned to the transient component of the multiscale formulation, this term is now examined meticulously.

If the constitutive law (22) is replaced within the expression for ${\mathrm{V}}_{\mathrm{t}\mathrm{r}\mathrm{a}}$, and reordering conveniently, it yields

$\begin{array}{l}{\mathrm{V}}_{\mathrm{t}\mathrm{r}\mathrm{a}}=-\frac{1}{\left|{\mathrm{\Omega }}_{\mathrm{\mu }}\right|}{\int }_{{\mathrm{\Omega }}_{\mathrm{\mu }}}{\dot{\mathrm{\chi }}}_{\mathrm{\mu }}\mathrm{y}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{\mu }}=-\frac{1}{\left|{\mathrm{\Omega }}_{\mathrm{\mu }}\right|}{\int }_{{\mathrm{\Omega }}_{\mathrm{\mu }}}{\widehat{\dot{\mathrm{\chi }}}}_{\mathrm{\mu }}^{\text{FOE}}\left({\dot{\mathrm{\epsilon }}}_{\mathrm{\mu }}^{\text{FOE}},{\dot{\mathrm{p}}}_{\mathrm{\mu }}^{\text{FOE}}\right)\mathrm{y}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{\mu }}\\ -\underset{{\mathrm{T}}_1}{\underbrace{\frac{1}{\left|{\mathrm{\Omega }}_{\mathrm{\mu }}\right|}{\int }_{{\mathrm{\Omega }}_{\mathrm{\mu }}}\left({\mathrm{b}}_{\mathrm{\mu }}:\dot{\mathrm{\epsilon }}\right)\mathrm{y}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{\mu }}}}-\underset{{\mathrm{T}}_2}{\underbrace{\frac{1}{\left|{\mathrm{\Omega }}_{\mathrm{\mu }}\right|}{\int }_{{\mathrm{\Omega }}_{\mathrm{\mu }}}\frac{\dot{\mathrm{p}}}{{\mathrm{M}}_{\mathrm{\mu }}}\mathrm{y}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{\mu }}}},\\ -\underset{{\mathrm{T}}_3\to \mathcal{O}\left({\mathrm{y}}^2\right)}{\underbrace{\frac{1}{\left|{\mathrm{\Omega }}_{\mathrm{\mu }}\right|}{\int }_{{\mathrm{\Omega }}_{\mathrm{\mu }}}\left[{\mathrm{b}}_{\mathrm{\mu }}:{\dot{\tilde{\mathrm{\epsilon }}}}_{\mathrm{\mu }}\left(\dots ,\dot{\mathrm{\phi }},\dots \right)\right]\mathrm{y}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{\mu }}}}-\underset{{\mathrm{T}}_4\to \mathcal{O}\left({\mathrm{y}}^2\right)}{\underbrace{\frac{1}{\left|{\mathrm{\Omega }}_{\mathrm{\mu }}\right|}{\int }_{{\mathrm{\Omega }}_{\mathrm{\mu }}}\frac{{\dot{\tilde{\mathrm{p}}}}_{\mathrm{\mu }}\left(\dots ,\dot{\mathrm{\phi }},\dots \right)}{{\mathrm{M}}_{\mathrm{\mu }}}\mathrm{y}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{\mu }}}}-\underset{{\mathrm{T}}_5\to \mathcal{O}\left({\mathrm{y}}^2\right)}{\underbrace{\frac{1}{\left|{\mathrm{\Omega }}_{\mathrm{\mu }}\right|}{\int }_{{\mathrm{\Omega }}_{\mathrm{\mu }}}\frac{\dot{\mathrm{\phi }}}{{\mathrm{M}}_{\mathrm{\mu }}}\cdot \left(\mathrm{y}\otimes \mathrm{y}\right)\mathrm{d}{\mathrm{\Omega }}_{\mathrm{\mu }}}}.\end{array}$ (24)

Assuming nearly symmetric distributions for ${\mathrm{b}}_{\mathrm{\mu }}$ and ${\mathrm{M}}_{\mathrm{\mu }}$ properties with respect to the micro-cell barycenter, the terms ${\mathrm{T}}_1$ and ${\mathrm{T}}_{2\mathrm{}}$in (24) do not introduce micro-scale size dependence problems. The effect of these contributions tends to decrease as the micro-cell size enlarges including a major (and representative) number of heterogeneities. So, these terms contribute to qualify when a generic micro-cell is (or not) an RVE.

The last ${\mathrm{T}}_5$-term of (24) has a quadratic dependence on the $\mathrm{y}$-coordinate. Clearly, such a term introduces a size effect in the homogenized response because it increases monotonically for enlarging micro-cell dimensions and seriously compromises the concept of RVE existence. The gradient rate of pore pressures ($\dot{\mathrm{\phi }}$) in the neighborhood of macro-scale external loads can take large values, especially if short-term analysis for low permeability and cohesive soils is considered.

Although it is not straightforward to see the terms ${\mathrm{T}}_3\mathrm{}$and ${\mathrm{T}}_4$ also introduce size effect issues whenever the FOE-multiscale formulation is employed. This is because the solution of the fluctuating components ${\dot{\tilde{\mathrm{\epsilon }}}}_{\mathrm{\mu }}$ and ${\dot{\tilde{\mathrm{p}}}}_{\mathrm{\mu }}$ (present in ${\mathrm{T}}_3\mathrm{}$and ${\mathrm{T}}_4$, respectively) depend on the inserted macro-scale quantity $\dot{\mathrm{\phi }}$, and such kind of implicit functionality gives rise to a second-order dependence on $\mathrm{y}$-coordinate for ${\mathrm{T}}_3\mathrm{}$and ${\mathrm{T}}_4$. The previous sentence can be proved in some academic multiscale scenarios by solving analytically the coupled system of differential equations.

3.3 Constitutive laws based on Reduced Order Expansions for primal descriptors

It has been shown that the lack of objectivity in the macro-scale response stems from the type of expansion accepted to evaluate input arguments in constitutive laws. This analysis allows us to introduce minimally invasive modifications, pointing to the core of the problem, whereas all the remainder features and ingredients of the proposed (variationally consistent) multiscale formulation are preserved, in order to not disclaim modelling capabilities.

For this aim, slight variations are now introduced in the definition of the stress-like functionals ${\widehat{\dot{\mathrm{\sigma }}}}_{\mathrm{\mu }}$ and ${\widehat{\dot{\mathrm{\chi }}}}_{\mathrm{\mu }}$. Thus, different order of expansions for the micro-scale strain rate ${\dot{\tilde{\mathrm{\epsilon }}}}_{\mathrm{\mu }}$ and pore-pressure rate ${\dot{\tilde{\mathrm{p}}}}_{\mathrm{\mu }}\mathrm{}$are proposed, for evaluating ${\widehat{\dot{\mathrm{\sigma }}}}_{\mathrm{\mu }}$ and ${\widehat{\dot{\mathrm{\chi }}}}_{\mathrm{\mu }}$ constitutive functionals. In this way, the following Selective Expansion choices are proposed: a) A Full Order Expansion (FOE) for the micro-scale strain rate descriptor, denoted as ${\dot{\mathrm{\epsilon }}}_{\mathrm{\mu }}^{\text{FOE}}$ identical to expression (9); b) A Reduced Order Expansion (ROE) for the micro-scale pore-pressures rate as a result of neglecting the first-order term ($\dot{\mathrm{\phi }}\cdot \mathrm{y}$) in (10), thus, ${\dot{\mathrm{p}}}_{\mathrm{\mu }}^{\text{ROE}}=\dot{\mathrm{p}}+{\dot{\tilde{\mathrm{p}}}}_{\mathrm{\mu }}$.

It is worth mentioning at this point that the ROE for ${\dot{\mathrm{p}}}_{\mathrm{\mu }}$ (${\dot{\mathrm{p}}}_{\mathrm{\mu }}^{\text{ROE}}$) has a very limited local effect in the multiscale formulation, since it has only meaningful within the constitutive functions characterizing ${\dot{\mathrm{\sigma }}}_{\mathrm{\mu }}$ and ${\dot{\mathrm{\chi }}}_{\mathrm{\mu }}$. Therefore, our proposal can be formally viewed as a simplifying constitutive-like hypothesis for ${\dot{\mathrm{\sigma }}}_{\mathrm{\mu }}$ and ${\dot{\mathrm{\chi }}}_{\mathrm{\mu }}$. Note that the generalized Darcy's law remains unchanged as defined in (20), thus the complete field ${\dot{\mathrm{p}}}_{\mathrm{\mu }}^{\text{FOE}}$ is required at this instance, and the constitutive equation for the micro-scale effective stress, ${\dot{\mathrm{\sigma }}}_{\mathrm{\mu }}^{\mathrm{\text{'}}}$, also remains invariable if compared with the FOE multiscale scheme.

4 SOLUTION OF THE VARIATIONAL EQUATIONS AT THE MACRO-SCALE

The global numerical paradigm consists of two nested, time-evolving, finite element schemes where the connection between them is established at each macro-scale Gauss point. In the literature, such an approach is referred to as FE² strategy, Feyel and Chaboche (2000).

The time variable, $\mathrm{t}$, is discretized through a monotonically increasing sequence of time steps $\left[{\mathrm{t}}^0,{\mathrm{t}}^1,{\mathrm{t}}^2,\dots ,{\mathrm{t}}^{\mathrm{n}},{\mathrm{t}}^{\mathrm{n}+1},\dots \right]$. The $\mathrm{è}$-generalized rule (Lewis and Schrefler (1998)) is used to account for the problem evolution within any time interval.

To describe the physical domain $\mathrm{Ù}$ a finite element mesh ${\mathrm{\Omega }}_{\mathrm{h}}$ is used. Then it is possible to build global interpolation matrices for the displacement, ${\mathrm{N}}_{\mathrm{u}}$, and pore pressure, ${\mathrm{N}}_{\mathrm{p}}$, fields as follows

$\mathrm{u}={\mathrm{N}}_{\mathrm{u}}\stackrel{-}{\mathrm{u}}\mathrm{}\mathrm{}\mathrm{}\mathrm{};\mathrm{}\mathrm{}\mathrm{}\mathrm{p}={\mathrm{N}}_{\mathrm{p}}\stackrel{-}{\mathrm{p}}\mathrm{},$ (25)

where $\stackrel{-}{\mathrm{u}}$ and $\stackrel{-}{\mathrm{p}}\mathrm{}$are the vectors that collect all nodal displacements and pore pressures values, respectively. In order to satisfy the Babuska-Brezzi convergence requirements (Lewis and Schrefler (1998)), it is necessary to adopt different orders of interpolation for each primal variable. Following the Galerkin method, the same spatial approximation is used for virtual actions. Replacing (25) into the weak forms (6), and after performing standard manipulations specific to the finite element method, the discrete (time and spatial) version of balance, at time step ${\mathrm{t}}^{\mathrm{n}+\mathrm{\theta }}$, can be written (Dirichlet degrees of freedom are omitted)

$\begin{array}{l}{\dot{\mathrm{G}}}_{\mathrm{h}}^{\mathrm{n}+\mathrm{\theta }}\equiv {\int }_{{\mathrm{\Omega }}_{\mathrm{h}}}{{\mathrm{B}}_{\mathrm{u}}}^{\mathrm{T}}\mathrm{}{\dot{\mathrm{\sigma }}}^{\mathrm{n}+\mathrm{\theta }}\mathrm{d}\mathrm{\Omega }\mathrm{}-{\int }_{{\mathrm{\Omega }}_{\mathrm{h}}}{{\mathrm{N}}_{\mathrm{u}}}^{\mathrm{T}}{\dot{\mathrm{f}}}^{\mathrm{n}+\mathrm{\theta }}\mathrm{d}\mathrm{\Omega }-{\int }_{{\mathrm{\Gamma }}_{\text{N},\mathrm{h}}^{\mathrm{u}}}{{\mathrm{N}}_{\mathrm{u}}}^{\mathrm{T}}{\dot{\mathrm{t}}}^{\mathrm{n}+\mathrm{\theta }}\mathrm{d}\mathrm{\Gamma }=0,\\ {\mathrm{H}}_{\mathrm{h}}^{\mathrm{n}+\mathrm{\theta }}\equiv {\int }_{{\mathrm{\Omega }}_{\mathrm{h}}}{{\mathrm{N}}_{\mathrm{p}}}^{\mathrm{T}}{\dot{\mathrm{\chi }}}^{\mathrm{n}+\mathrm{\theta }}\mathrm{}\mathrm{d}\mathrm{\Omega }-{\int }_{{\mathrm{\Omega }}_{\mathrm{h}}}{{\mathrm{B}}_{\mathrm{p}}}^{\mathrm{T}}{\mathrm{V}}^{\mathrm{n}+\mathrm{\theta }}\mathrm{d}\mathrm{\Omega }+{\int }_{{\mathrm{\Gamma }}_{\text{N},\mathrm{h}}^{\mathrm{p}}}{{\mathrm{N}}_{\mathrm{p}}}^{\mathrm{T}}\frac{{\mathrm{q}}^{\mathrm{n}+\mathrm{\theta }}}{{\mathrm{\rho }}^{\mathrm{f}}}\mathrm{d}\mathrm{\Gamma }=0.\end{array}$ (26)

being ${\mathrm{B}}_{\mathrm{u}}={\nabla }^{\mathrm{s}\mathrm{y}\mathrm{m}}{\mathrm{N}}_{\mathrm{u}}$ the deformation-displacement matrix and ${\mathrm{B}}_{\mathrm{p}}=\nabla {\mathrm{N}}_{\mathrm{p}}$ is the matrix that relates pore pressures to their gradients. The system of equations (26) is solved by using a standard Newton-Raphson algorithm. Thus, the current displacement and pore pressure nodal values are updated (at time step ${\mathrm{t}}^{\mathrm{n}+1}$) in terms of the iterative increments $\mathrm{Ä}\stackrel{-}{\mathrm{u}}$ and $\mathrm{Ä}\stackrel{-}{\mathrm{p}}$, respectively. For a given $\mathrm{k}$-iteration in the time interval of the Newton-Raphson scheme, these increments are evaluated as (subscript $\mathrm{h}$ is omitted hereafter)

$\left[\begin{array}{c}\mathrm{\Delta }\stackrel{-}{\mathrm{u}}\\ \mathrm{\Delta }\stackrel{-}{\mathrm{p}}\end{array}\right]=-\left({\mathrm{J}}_{\mathrm{k}}^{\mathrm{n}+\mathrm{\theta }}\right)\left[\begin{array}{c}{\dot{\mathrm{G}}}_{\mathrm{h}}^{\mathrm{n}+\mathrm{\theta }}\\ {\mathrm{H}}_{\mathrm{h}}^{\mathrm{n}+\mathrm{\theta }}\end{array}\right]\mathrm{}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\mathrm{}{\mathrm{J}}^{\mathrm{n}+\mathrm{\theta }}\mathrm{}=\mathrm{}\left[\begin{array}{cc}\frac{\partial {\dot{\mathrm{G}}}_{\mathrm{h}}^{\mathrm{n}+\mathrm{\theta }}}{\partial {\stackrel{-}{\mathrm{u}}}^{\mathrm{n}+1}}& \frac{\partial {\dot{\mathrm{G}}}_{\mathrm{h}}^{\mathrm{n}+\mathrm{\theta }}}{\partial {\stackrel{-}{\mathrm{p}}}^{\mathrm{n}+1}}\\ \frac{\partial {\mathrm{H}}_{\mathrm{h}}^{\mathrm{n}+\mathrm{\theta }}}{\partial {\stackrel{-}{\mathrm{u}}}^{\mathrm{n}+1}}& \frac{\partial {\mathrm{H}}_{\mathrm{h}}^{\mathrm{n}+\mathrm{\theta }}}{\partial {\stackrel{-}{\mathrm{p}}}^{\mathrm{n}+1}}\end{array}\right]$ (27)

Each component of the macro-scale Jacobian matrix ${\mathrm{J}}_{\mathrm{k}}^{\mathrm{n}+\mathrm{\theta }}$ can be found in Khoei and Hajiabadi (2018) and Anonis et al. (2024a).

5 SOLUTION OF THE VARIATIONAL EQUATIONS AT THE MICRO-SCALE

Macro and micro spatial domains share the same temporal scale. Thus, all definitions concerning to the time discretization scheme, described in the section before are still valid at the micro-scale level.

The same finite element technology used in the context of the macro-scale problem is used to conform a mesh, ${\mathrm{\Omega }}_{\mathrm{\mu },\mathrm{h}}$, for the micro-scale domain ${\mathrm{\Omega }}_{\mathrm{\mu }}$. Then, the primal descriptors ${\dot{\tilde{\mathrm{u}}}}_{\mathrm{\mu }}$ and ${\tilde{\mathrm{p}}}_{\mathrm{\mu }}$ (and their corresponding admissible virtual variations $\mathrm{ä}{\dot{\tilde{\mathrm{u}}}}_{\mathrm{\mu }}$ and $\mathrm{ä}{\tilde{\mathrm{p}}}_{\mathrm{\mu }}$) can be approached through global shape function matrices for micro-displacement fluctuations, ${\mathrm{N}}_{{\tilde{\mathrm{u}}}_{\mathrm{\mu }}}$, and micro-pore pressure fluctuations, ${\mathrm{N}}_{{\tilde{\mathrm{p}}}_{\mathrm{\mu }}}$, as

$\begin{array}{l}{\dot{\tilde{\mathrm{u}}}}_{\mathrm{\mu }}={\mathrm{N}}_{{\tilde{\mathrm{u}}}_{\mathrm{\mu }}}{\dot{\widetilde{\stackrel{-}{\mathrm{u}}}}}_{\mathrm{\mu }},\mathrm{\delta }{\tilde{\mathrm{u}}}_{\mathrm{\mu }}={\mathrm{N}}_{{\tilde{\mathrm{u}}}_{\mathrm{\mu }}}\mathrm{\delta }{\dot{\widetilde{\stackrel{-}{\mathrm{u}}}}}_{\mathrm{\mu }},\mathrm{}\mathrm{}\text{with}{\dot{\widetilde{\stackrel{-}{\mathrm{u}}}}}_{\mathrm{\mu }}\text{and}\mathrm{\delta }{\dot{\widetilde{\stackrel{-}{\mathrm{u}}}}}_{\mathrm{\mu }}\in {\tilde{\mathcal{U}}}_{\mathrm{\mu },\mathrm{h}}^{\text{Per}},\\ {\tilde{\mathrm{p}}}_{\mathrm{\mu }}={\mathrm{N}}_{{\tilde{\mathrm{p}}}_{\mathrm{\mu }}}{\widetilde{\stackrel{-}{\mathrm{p}}}}_{\mathrm{\mu }},\mathrm{\delta }{\tilde{\mathrm{p}}}_{\mathrm{\mu }}={\mathrm{N}}_{{\tilde{\mathrm{p}}}_{\mathrm{\mu }}}\mathrm{\delta }{\widetilde{\stackrel{-}{\mathrm{p}}}}_{\mathrm{\mu }},\mathrm{}\mathrm{}\text{with}{\widetilde{\stackrel{-}{\mathrm{p}}}}_{\mathrm{\mu }}\text{and}\mathrm{\delta }{\widetilde{\stackrel{-}{\mathrm{p}}}}_{\mathrm{\mu }}\in {\tilde{\mathcal{P}}}_{\mathrm{\mu },\mathrm{h}}^{\text{Per}}.\end{array}$ (28)

where ${\widetilde{\stackrel{-}{\mathrm{u}}}}_{\mathrm{\mu }}$ and ${\widetilde{\stackrel{-}{\mathrm{p}}}}_{\mathrm{\mu }}$ denote the vectors containing all nodal micro-scale displacement fluctuations and pore pressure fluctuations, respectively, while ${\tilde{\mathcal{U}}}_{\mathrm{\mu },\mathrm{h}}^{\text{Per}}$ and ${\tilde{\mathcal{P}}}_{\mathrm{\mu },\mathrm{h}}^{\text{Per}}$ represent the finite-dimensional counterparts of sub-spaces ${\tilde{\mathcal{U}}}_{\mathrm{\mu }}^{\text{Per}}$ and ${\tilde{\mathcal{P}}}_{\mathrm{\mu }}^{\text{Per}}$ related to the Periodic multiscale model (see, namely de Souza Neto and Feijoo (2006), Sanchez et al. (2013) and de Souza Neto et al. (2015), Peric et al. (2011)).

After some standard mathematical treatment, the substitution of (28) into the variational forms Eqs. (18)-(19) leads to the spatial and time discrete approximation (at time $t^{n+\theta }$) for the balance equations in the micro-scale

$\underset{{\dot{\mathrm{G}}}_{\mathrm{\mu },\mathrm{h}}^{\mathrm{n}+\mathrm{\theta }}}{\underbrace{\left[{\int }_{{\mathrm{\Omega }}_{\mathrm{\mu },\mathrm{h}}}\left({\mathrm{B}}_{{\tilde{\mathrm{u}}}_{\mathrm{\mu }}}^{\mathrm{T}}{\dot{\mathrm{\sigma }}}_{\mathrm{\mu }}^{\mathrm{n}+\mathrm{\theta }}-{\mathrm{N}}_{{\tilde{\mathrm{u}}}_{\mathrm{\mu }}}^{\mathrm{T}}{\dot{\mathrm{f}}}_{\mathrm{\mu }}^{\mathrm{n}+\mathrm{\theta }}\right)\mathrm{d}{\mathrm{\Omega }}_{\mathrm{\mu }}\right]}}\cdot \mathrm{\delta }{\dot{\widetilde{\stackrel{-}{\mathrm{u}}}}}_{\mathrm{\mu }}=0,\forall \mathrm{\delta }{\dot{\widetilde{\stackrel{-}{\mathrm{u}}}}}_{\mathrm{\mu }}\in {\tilde{\mathcal{U}}}_{\mathrm{\mu },\mathrm{h}}^{\text{Per}},\mathrm{}{\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\mathrm{}\dot{\widetilde{\stackrel{-}{\mathrm{u}}}}}_{\mathrm{\mu }}\in {\tilde{\mathcal{U}}}_{\mathrm{\mu },\mathrm{h}}^{\text{Per}}\mathrm{}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{}{\widetilde{\stackrel{-}{\mathrm{p}}}}_{\mathrm{\mu }}\in {\tilde{\mathcal{P}}}_{\mathrm{\mu },\mathrm{h}}^{\text{Per}},$ (29)

$\underset{{\mathrm{H}}_{\mathrm{\mu },\mathrm{h}}^{\mathrm{n}+\mathrm{\theta }}}{\underbrace{\left[{\int }_{{\mathrm{\Omega }}_{\mathrm{\mu },\mathrm{h}}}\left({\mathrm{N}}_{{\tilde{\mathrm{p}}}_{\mathrm{\mu }}}^{\mathrm{T}}{\dot{\mathrm{\chi }}}_{\mathrm{\mu }}^{\mathrm{n}+\mathrm{\theta }}-{\mathrm{B}}_{{\tilde{\mathrm{p}}}_{\mathrm{\mu }}}^{\mathrm{T}}{\mathrm{V}}_{\mathrm{\mu }}^{\mathrm{n}+\mathrm{\theta }}\right)\mathrm{d}{\mathrm{\Omega }}_{\mathrm{\mu }}\right]}}\cdot \mathrm{\delta }{\widetilde{\stackrel{-}{\mathrm{p}}}}_{\mathrm{\mu }}=0,\forall \mathrm{\delta }{\widetilde{\stackrel{-}{\mathrm{p}}}}_{\mathrm{\mu }}\in {\tilde{\mathcal{P}}}_{\mathrm{\mu },\mathrm{h}}^{\text{Per}},\mathrm{}{\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\mathrm{}\dot{\widetilde{\stackrel{-}{\mathrm{u}}}}}_{\mathrm{\mu }}\in {\tilde{\mathcal{U}}}_{\mathrm{\mu },\mathrm{h}}^{\text{Per}}\mathrm{}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{}{\widetilde{\stackrel{-}{\mathrm{p}}}}_{\mathrm{\mu }}\in {\tilde{\mathcal{P}}}_{\mathrm{\mu },\mathrm{h}}^{\text{Per}},$ (30)

where ${\mathrm{B}}_{{\tilde{\mathrm{u}}}_{\mathrm{\mu }}}^{\mathrm{T}}={\nabla }_{\mathrm{y}}^{\mathrm{s}\mathrm{y}\mathrm{m}}{\mathrm{N}}_{{\tilde{\mathrm{u}}}_{\mathrm{\mu }}}$ and ${\mathrm{B}}_{{\tilde{\mathrm{p}}}_{\mathrm{\mu }}}^{\mathrm{T}}={\nabla }_{\mathrm{y}}{\mathrm{N}}_{{\tilde{\mathrm{p}}}_{\mathrm{\mu }}}$, are the micro-scale global matrices relating the primal variables with their corresponding gradients.

The numerical solution for nodal values of displacements, ${\widetilde{\stackrel{-}{\mathrm{u}}}}_{\mathrm{\mu }}^{\mathrm{n}+1}$, and pore pressure fluctuations, ${\widetilde{\stackrel{-}{\mathrm{p}}}}_{\mathrm{\mu }}^{\mathrm{n}+1}$, at current time step ${\mathrm{t}}^{\mathrm{n}+1}$, is obtained in terms of the iterative increments $∆{\widetilde{\stackrel{-}{\mathrm{u}}}}_{\mathrm{\mu }}$ and ${∆\widetilde{\stackrel{-}{\mathrm{p}}}}_{\mathrm{\mu }}$, respectively. They are computed by using a standard Newton-Raphson procedure applied to the system of Eqs. (29)-(30). In this sense, for a generic $\mathrm{k}$-iteration it can be expressed (subscript $\mathrm{h}$ is omitted hereafter)

$\left\{{\mathrm{J}}_{\mathrm{\mu },\mathrm{k}}^{\mathrm{n}+\mathrm{\theta }}\left[\begin{array}{c}∆{\widetilde{\stackrel{-}{\mathrm{u}}}}_{\mathrm{\mu }}\\ {∆\widetilde{\stackrel{-}{\mathrm{p}}}}_{\mathrm{\mu }}\end{array}\right]+\left[\begin{array}{c}{\dot{\mathrm{G}}}_{\mathrm{\mu },\mathrm{h}}^{\mathrm{n}+\mathrm{\theta }}\\ {\mathrm{H}}_{\mathrm{\mu },\mathrm{h}}^{\mathrm{n}+\mathrm{\theta }}\end{array}\right]\right\}\cdot \left[\begin{array}{c}\mathrm{\delta }{\dot{\widetilde{\stackrel{-}{\mathrm{u}}}}}_{\mathrm{\mu }}\\ \mathrm{\delta }{\widetilde{\stackrel{-}{\mathrm{p}}}}_{\mathrm{\mu }}\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right],\forall \mathrm{\delta }{\dot{\widetilde{\stackrel{-}{\mathrm{u}}}}}_{\mathrm{\mu }}\in {\tilde{\mathcal{U}}}_{\mathrm{\mu },\mathrm{h}}^{\text{Per}},\mathrm{}\forall \mathrm{\delta }{\widetilde{\stackrel{-}{\mathrm{p}}}}_{\mathrm{\mu }}\in {\tilde{\mathcal{P}}}_{\mathrm{\mu },\mathrm{h}}^{\text{Per}}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}∆{\widetilde{\stackrel{-}{\mathrm{u}}}}_{\mathrm{\mu }}\in {\tilde{\mathcal{U}}}_{\mathrm{\mu },\mathrm{h}}^{\text{Per}}\mathrm{}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{}{∆\widetilde{\stackrel{-}{\mathrm{p}}}}_{\mathrm{\mu }}\in {\tilde{\mathcal{P}}}_{\mathrm{\mu },\mathrm{h}}^{\text{Per}},$ (31)

${\mathrm{J}}_{\mathrm{\mu },\mathrm{k}}^{\mathrm{n}+\mathrm{\theta }}$ being the Jacobian operator in the micro-scale, which has the form (subscript $\mathrm{k}$ is removed to simplify the notation)

${\mathrm{J}}_{\mathrm{\mu },\mathrm{k}}^{\mathrm{n}+\mathrm{\theta }}=\left[\begin{array}{cc}\frac{\partial {\dot{\mathrm{G}}}_{\mathrm{\mu },\mathrm{h}}^{\mathrm{n}+\mathrm{\theta }}}{\partial {\widetilde{\stackrel{-}{\mathrm{u}}}}_{\mathrm{\mu }}^{\mathrm{n}+1}}& \frac{\partial {\dot{\mathrm{G}}}_{\mathrm{\mu },\mathrm{h}}^{\mathrm{n}+\mathrm{\theta }}}{\partial {\widetilde{\stackrel{-}{\mathrm{p}}}}_{\mathrm{\mu }}^{\mathrm{n}+1}}\\ \frac{\partial {\mathrm{H}}_{\mathrm{\mu },\mathrm{h}}^{\mathrm{n}+\mathrm{\theta }}}{\partial {\widetilde{\stackrel{-}{\mathrm{u}}}}_{\mathrm{\mu }}^{\mathrm{n}+1}}& \frac{\partial {\mathrm{H}}_{\mathrm{\mu },\mathrm{h}}^{\mathrm{n}+\mathrm{\theta }}}{\partial {\widetilde{\stackrel{-}{\mathrm{p}}}}_{\mathrm{\mu }}^{\mathrm{n}+1}}\end{array}\right]=\left[\begin{array}{cc}\frac{1}{∆\mathrm{t}}{\mathrm{K}}_{\mathrm{\mu }}& \frac{1}{∆\mathrm{t}}{\mathrm{Q}}_{\mathrm{\mu }}\\ \frac{1}{∆\mathrm{t}}{\left({\mathrm{Q}}_{\mathrm{\mu }}\right)}^{\mathrm{T}}& \left(\frac{1}{∆\mathrm{t}}{\mathrm{S}}_{\mathrm{\mu }}+\mathrm{\theta }{\mathrm{\kappa }}_{\mathrm{\mu }}\right)\end{array}\right]$ (32)

The matrices ${\mathrm{K}}_{\mathrm{\mu }}$, ${\mathrm{Q}}_{\mathrm{\mu }}$, ${\mathrm{S}}_{\mathrm{\mu }}$ and ${\mathrm{\kappa }}_{\mathrm{\mu }}$, that conform the Jacobian in (34), can be found in Anonis et al. (2024a) and Anonis et al. (2024b).

6 SCOPE OF CONSTITUTIVE FORMULATIONS BASED ON ROE

We present below a brief discussion on the possible scope or direct extention of our formulation when applied to multiscale modeling based on the RVE concept of different physical phenomena such as transient heat conduction (Özdemir et al. (2008), Larsson et al. (2010a), Ramos et al. (2017) and Waseem et al. (2020)) or mass conservation of ion species (coupled to other phenomena) (Kaessmair et al. (2021) and Saeedmonir et al. (2024), among others. Although the present work deals with multiscale modeling of porous media, the resemblance between the mass balance equation presented here and the balance of other physical quantities allows us to advance in this analysis.

Thus, in the author's opinion, whenever one attempts to model by means of multiscale methodologies based on RVE physical phenomena that satisfy the conditions mentioned in the preceding section 3.2, there will be a dependence on the RVE size and loss of objectivity in the macro-scale response.

Such conditions are:

i. That both the macro- and micro-scale are described by strong-form of mass governing equations of the type

$\begin{array}{l}\dot{\mathrm{\chi }}+\nabla \cdot \mathrm{w}=0,\\ {\dot{\mathrm{\chi }}}_{\mathrm{\mu }}\left({\dot{\mathrm{p}}}_{\mathrm{\mu }}\right)+{\nabla }_{\mathrm{\mu }}\cdot {\mathrm{w}}_{\mathrm{\mu }}=0.\end{array}$ (33)

ii. there is a constitutive functional dependence of the type ${\dot{\mathrm{\chi }}}_{\mathrm{\mu }}\left({\dot{\mathrm{p}}}_{\mathrm{\mu }}^{\text{FOE}}\right)={\dot{\widehat{\mathrm{\chi }}}}_{\mathrm{\mu }}\left(\dot{\mathrm{p}}+\dot{\mathrm{\phi }}\cdot \mathrm{y}+{\dot{\tilde{\mathrm{p}}}}_{\mathrm{\mu }}\right)$ (only the functional dependence of ${\dot{\widehat{\mathrm{\chi }}}}_{\mathrm{\mu }}$ is expressed because it is the relevant for the purpose of this discussion).

The reason, is that under these circumstances expressions like those given in (15) are reached, which in the presence of a full order expansion of ${\dot{\mathrm{p}}}_{\mathrm{\mu }}^{\text{FOE}}$ will lead to the defect explained at length in the preceding sections. On the other hand, if one assumes the constitutive ROE functional ${\dot{\mathrm{\chi }}}_{\mathrm{\mu }}\left({\dot{\mathrm{p}}}_{\mathrm{\mu }}^{\text{ROE}}\right)={\dot{\mathrm{\chi }}}_{\mathrm{\mu }}\left(\dot{\mathrm{p}}+{\dot{\tilde{\mathrm{p}}}}_{\mathrm{\mu }}\right)$, one recovers the objectivity of the response.

However, equations of the type (15) can be employed in the analysis of so many other physical phenomena individually or coupled as we mentioned at the beginning, for example the cases of transient heat conduction with $\mathrm{w}$ being the heat flux field and ${\dot{\mathrm{\chi }}}_{\mathrm{\mu }}$ the specific internal energy rate dependent on the temperature field or $\mathrm{w}$ the molar flux of species and ${\dot{\mathrm{\chi }}}_{\mathrm{\mu }}$ the ion concentration rate, etc.

Therefore, we argue that the ROE constitutive concept is applicable to all these cases and its scope goes beyond the one proposed here.

7 CONSOLIDATION EXAMPLES

The macro-scale domain consists of a soil column with a height of H = 1 m and a width of W = 0.1 m. At this scale, the material is treated as a homogenized medium. The boundary conditions are illustrated in Figure 1. At the ground surface, a compressive external traction is applied, increasing from 0 kPa to -100 kPa during the first day of analysis, and thereafter remaining constant. The macro-scale domain is discretized into 13 finite elements, as shown in Figure 1(b). The time discretization employed throughout the analysis follows the sequence: t = [0.001, 0.002, …, 0.009, 0.01, 0.02, …, 0.09, 0.1, 0.2, …, 0.9, 1, 2, …, 9, 10, 20, …, 70, 80] (in days). At the micro-scale, the periodic multiscale model is adopted. The heterogeneous material properties are defined as follows. For the matrix: Young’s modulus E = 100 kPa, Poisson’s ratio ν = 0.3, isotropic hydraulic conductivity k = 8.64 × 10−4 m/day, and initial void ratio e0 = 0.6. For the inclusion: Young’s modulus E = 1500 kPa, Poisson’s ratio ν = 0.3, isotropic hydraulic conductivity k = 1.296 × 10−3 m/day, and initial void ratio e0 = 2.0. Different micro-cells (MCs) of increasing size are considered to obtain homogenized responses using both FOE and ROE formulations. Specifically, square MCs of sizes 0.1 m × 0.1 m, 0.5 m × 0.5 m, and 1 m × 1 m are analyzed. Each MC is discretized with 36, 144, and 576 finite elements, respectively. For the Direct Numerical Simulation (DNS), see Figure 1(a), the spatial discretization required a mesh of approximately 5760 finite elements (mesh not drawing).

Figure 1
Schematic representation of the problem to be solved

Figure 2 presents the time-evolution curves of the vertical component of the effective stress tensor at the reference area A_Ref. The results reveal that micro-scale size sensitivity manifests in the homogenized stress response when the FOE-based multiscale procedure is employed. The FOE model systematically underestimates the effective stress compared to DNS, with discrepancies becoming significant for the two largest micro-cells. In contrast, the SOE-based homogenization yields effective stress curves that are in good agreement with the DNS results, independently of the micro-cell size.

Figure 2
Time-evolution curves of the vertical component of the effective stress tensor at the reference area A_Ref.

The size effect associated with the FOE procedure is even more pronounced for the homogenized seepage velocity. This is illustrated in Figure 3, where the temporal evolution of seepage velocity is compared among the different models at ARef. As expected, for the two largest micro-cells the FOE scheme produces results that markedly deviate from the DNS reference solution, while the SOE formulation provides consistent predictions.

Figure 3
Time-evolution curves of the vertical component of the seepage velocity at the reference area A_Ref.

8 CONCLUSIONS

A multiscale model of saturated porous media based on RVE has been presented. In the literature, multiscale models applied to saturated porous media have revealed some inconsistency with respect to the dependence of the micro-scale size. This drawback has been attributed to a second-order term. More specifically, to the seepage or flux velocity. The consistent PMVP-based homogenization model adopted in this work has allowed us to identify that the size effect of the RVE is caused only by one component of the flux velocity term. This component comes from the use of a full expansion of the micro-scale pore pressure at the constitutive level on the fine scale. Thus, a redefinition of the constitutive equations at the micro-scale has been proposed. This is based on a reduced-order expansion of the micro-scale pore pressure field to overcome the problem of dependence on the size of the RVE. But in addition, it is possible to maintain the dynamical term that arises during the homogenization process, and finally, it also preserves the principle of scale separation.

Acknowledgments

This work was partially supported by the Consejo Nacional de Investigaciones Cientficas y Tecnicas (CONICET, Argentina, grant PIP 2023-2025-GI-11220220100332CO) as well as Secretaria de General de Ciencia y Tecnica de la Universidad Nacional del Nordeste (UNNE), Argentina (Grants PI 21D003). Also, the authors would like to acknowledge the financial support from the Agencia Nacional de Promoción Científica y Tecnológica (ANPCyT, Argentina, grants PICTO-UNNE-2019-00014 and PICT 2021-GRF-TI-00488).

• Data availability statement:
  Research data is available in the body of the article

References

• Anonis, R. A., Mroginski, J. L., Sánchez, P. J, (2024a). Multiscale formulation for saturated porous media preserving the representative volume element size objectivity. Int. J. Numer. Meth. Eng., 125(3):e7381.
• Anonis, R. A., Mroginski, J. L., Sánchez, P. J, (2024b). Multiscale formulation for materials composed by a saturated porous matrix and solid inclusions. Comput. Meth. Appl. Mech., 429:117162.
• Beneyto, P. A., Di Rado, H. A., Mroginski, J. L., Awruch, A. M., (2015). A versatile mathematical approach for environmental geomechanic modelling based on stress state decomposition. Appl. Math. Model., 39(22):6880–6896.
• Biot, M. A., (1955). Theory of elasticity and consolidation for a porous anisotropic solid. Journal of Applied Physics, 26(2): 182–185.
• Biot, M. A., (1955). General theory of three-dimensional consolidation. Journal of Applied Physics, 12(2):155–164.
• Blanco, P. J., Sánchez, P. J, de Souza Neto, E., Feijóo, R. A., (2014). Variational foundations and generalized unified theory of rve-based multiscale models. Arch. Comput. Method. E., 23.
• Coussy, O., (2003). Poromechanics. John Wiley and Sons, Ltd.
• de Souza Neto, E., Feijóo, R. A., (2006). Variational foundations of multi-scale constitutive models of solid: small and large strain kinematical formulation. LNCC-MCTI Laboratorio Nacional de Computacao Científica. Internal Report no. 16/2006. LNCC, Petropolis, Brazil.
• de Souza Neto, E., Blanco, P. J., Sánchez, P. J, Feijóo, R. A., (2015). An rve-based multiscale theory of solids with micro-scale inertia and body force effects. Mech. Mater., 80:136–144.
• Di Rado, H. A., Beneyto, P. A., Mroginski, J. L., Awruch, A. M., (2009). Influence of the saturation-suction relationship in the formulation of non-saturated soils consolidation models. Math. Comput. Model., 49(5-6):1058–1070.
• Dormieux, L., Kondo, D., Ulm, F. J., (2006). Microporomechanics. Microporomechanics.
• Feyel, F., Chaboche, J. L., (2000). Fe2 multiscale approach for modelling the elastoviscoplastic behaviour of long fibre sic/ti composite materials. Comput. Meth. Appl. Mech., 183(3):309–330.
• Hill, R., (1963). Elastic properties of reinforced solids: Some theoretical principles. Journal of the Mechanics and Physics of Solids; 11(5): 357-372. doi: https://doi.org/10.1016/0022-5096(63)90036-X.
  » https://doi.org/10.1016/0022-5096(63)90036-X
• Hill, R., (1965a). A self-consistent mechanics of composite materials. Journal of the Mechanics and Physics of Solid; 13(4):213-222. doi: https://doi.org/10.1016/0022-5096(65)90010-4.
  » https://doi.org/10.1016/0022-5096(65)90010-4
• Hill, R., (1965b). Continuum micro-mechanics of elastoplastic polycrystals. Journal of the Mechanics and Physics of Solids; 13(2): 89-101. doi: https://doi.org/10.1016/0022-5096(65)90023-2.
  » https://doi.org/10.1016/0022-5096(65)90023-2
• Janicke, R., Larsson, F., Runesson K., (2020). A poro-viscoelastic substitute model of fine-scale poroelasticity obtained from homogenization and numerical model reduction. Comput. Mech., 65.
• Kaessmair, S., Runesson K., Steinmann, P., Janicke, R., Larsson, F., (2021). Variationally consistent computational homogenization of chemomechanical problems with stabilized weakly periodic boundary conditions. Int. J. Numer. Meth. Eng., 122(22):6429–6454.
• Khoei, A. R., Saeedmonir, S., Bonabi, A. M., (2023). Computational homogenization of fully coupled hydro-mechanical analysis of micro-fractured porous media. Comput. Geotech., 154:105121.
• Khoei, A. R., Hajiabadi M.R., (2018). Fully coupled hydromechanical multiscale model with microdynamic effects. Int. J. Numer. Meth. Eng., 115(3):293–327.
• Khoei, A. R., Saeedmonir, S., (2021). Computational homogenization of fully coupled multiphase flow in deformable porous media. Comput. Meth. Appl. Mech., 376:113660.
• Klahr, B., Medeiros Thiesen, J. L., Teixeira Pinto, O., Carniel, T. A., Fancello, E. A., (2023). A variational rve-based multiscale poromechanical formulation applied to soft biological tissues under large deformations. Eur. J. Mech. A-Solid, 99:104937.
• Larsson, F., Runesson K., Su, F., (2010a). Computational homogenization of uncoupled consolidation in micro-heterogeneous porous media. Int. J. Numer. Anal. Meth. Geomech., 34(14):1431–1458.
• Larsson, F., Runesson K., Su, F., (2010b). Variationally consistent computational homogenization of transient heat flow. Int. J. Numer. Meth. Eng., 81(13):1659–1686.
• Lewis, R., Schrefler, B., (1998). The Finite Element Method in the Static and Dynamic Deformation and Consolidation in Porous Media, volume 34. John Wiley and Sons, Ltd.
• Lopez Rivarola, F., Labanda, N., Etse, G., (2019). Thermodynamically consistent multiscale homogenization for thermo-poroplastic materials. Zeitschrift fur angewandte Mathematik und Physik ZAMP, 70:1–27.
• Mroginski, J. L., Di Rado, H. A., Beneyto, P. A., Awruch, A. M., (2010). A finite element approach for multiphase fluid flow in porous media. Mathematics and Computers in Simulation (MATCOM), 81(1):76–91.
• Peric, D., de Souza Neto, E. A., Feijóo, R. A., Partovi, M., Carneiro Molina, A. J., (2011). On micro-to-macro transitions for multi-scale analysis of non-linear heterogeneous materials: unified variational basis and finite element implementation. Int. J. Numer. Meth. Eng., 87(1-5):149–170.
• Ramos, G. R., dos Santos, T., Rossi, R, (2017). An extension of the hill–mandel principle for transient heat conduction in heterogeneous media with heat generation incorporating finite rve thermal inertia effects. Int. J. Numer. Meth. Eng., 111(6):553–580.
• Saeedmonir, S., Adeli, M.H., Khoei, A.R., (2024). A multiscale approach in modeling of chemically reactive porous media. Comput. Geotech., 165:105818.
• Saeedmonir, S., Khoei, A. R., (2022). Multiscale modeling of coupled thermo-hydro-mechanical analysis of heterogeneous porous media. Comput. Meth. Appl. Mech., 391:114518.
• Su, F., Larsson, F., Runesson K., (2011). Computational homogenization of coupled consolidation problems in micro-heterogeneous porous media. Int. J. Numer. Meth. Eng., 88(11):1198–1218.
• Sánchez, P. J, Blanco, P. J., Huespe, A.E., Feijóo, R. A, (2013). Failure-oriented multi-scale variational formulation: Micro-structures with nucleation and evolution of softening bands. Comput. Meth. Appl. Mech., 257:221–247.
• Waseem, A., Heuze, T., Stainier, L., Geers, M. G. D., Kouznetsova, V. G., (2019). Model reduction in computational homogenization for transient heat conduction. Comput. Mech., 65.
• Wu, W., Yang, Y., Zheng, H., Zhang, L., Zhang, M., (2022). Numerical manifold computational homogenization for hydro-dynamic analysis of discontinuous heterogeneous porous media. Comput. Meth. Appl. Mech., 388: 114254.
• Özdemir, I., Brekelmans, W. A. M., Geers, M. G. D., (2008). Computational homogenization for heat conduction in heterogeneous solids. Int. J. Numer. Meth. Eng., 73(2):185–204.

Edited by

• Editor:
  Eduardo Alberto Fancello and Paulo de Tarso Mendonça

Data availability

Research data is available in the body of the article

Publication Dates

• Publication in this collection
  01 Dec 2025
• Date of issue
  2026

History

• Received
  29 July 2025
• Reviewed
  09 Sept 2025
• Accepted
  01 Oct 2025