# Second-Order Computational Homogenization of Nonlinear Fluid Flow through Porous Media

Elten Polukhov, Marc-André Keip

<figure><img src="image_1.png" /></figure>

## To cite this version:

Elten Polukhov, Marc-André Keip. Second-Order Computational Homogenization of Nonlinear Fluid Flow through Porous Media. 2025. hal-05040937

HAL Id: hal-05040937

https://hal.science/hal-05040937v1

Preprint submitted on 20 Apr 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Second-Order Computational Homogenization of Nonlinear Fluid Flow through Porous Media

Elten Polukhov* and Marc-André Keip†

University of Stuttgart
Institute of Applied Mechanics, Chair of Materials Theory
Pfaffenwaldring 7, 70569 Stuttgart, Germany

Preprint, April 11, 2025

## Abstract

We present a second-order computational multiscale model for heterogeneous porous media, which allows for the scale bridging of transient fluid-flow processes through porous materials with non-separated scales. The formulation connects a homogenized macroscopic scale described by a local theory of grade two with heterogeneous microstructures described by a local theory of grade one. At the macroscale this leads to $C^1$-continuity requirements on the macroscopic solution field; at the microscale we need to take account of constraints on fluctuation fields that require $H(\text{div}) \cap H(\text{grad})$-conformity of microscopic solutions. Both these challenges are addressed through the development of mixed Hu-Washizu formulations that result in a variationally consistent homogenization framework with minimization structure across scales. We validate the second-order multiscale model by means of fully resolved, direct numerical simulations and provide comparisons with results of first-order FE² simulations. By considering linear Darcy and nonlinear Darcy-Forchheimer flow through two- and three-dimensional porous microstructures, we provide further insights into the framework and associated length-scale effects.

**Keywords:** Second-order homogenization; Porous media; Variational formulation; Size effects; Darcy-Forchheimer flow

## 1 Introduction

Porous materials are ubiquitous in nature. Their multi-field coupling properties furthermore make them ideal candidates for various engineering applications, for example, in the areas of artificial implants (Byrne and Preziosi, 2003), drug delivery (Pałka and Pokrowiecki, 2018), geothermal energy (Kumari

<footnote>*Corresponding author (elten.polukhov@mechbau.uni-stuttgart.de).</footnote>

<footnote>†Corresponding author (marc-andre.keip@mechbau.uni-stuttgart.de).</footnote>
et al., 2018), battery technology (Oh et al., 2014) and sub-surface fluid flow (Mercer and Cohen, 1990). To better understand the underlying multi-physical processes and to facilitate and enhance associated applications reliable modeling approaches are essential.

The continuum modeling of porous materials traces back to the seminal works of Darcy (1856), Biot (1941), and Terzaghi (1943). Since then, a wealth of models and underlying theories has emerged (see the review of Ehlers (2023) for more information). Almost all of these models have in common that they treat the porous material as a mixture of different phases, like a solid and a fluid phase. These phases are then equipped with phase-specific properties of, for example, mechanical, thermal, and chemical nature. The classical theories of porous media address situations in which the scale of application is very much larger than the underlying microstructure, so that classical models are usually built on local, first-order descriptions. More recently, however, the view has widened to scenarios in which the microscopic and macroscopic scales are no longer clearly separated. Associated situations could arise generally in two ways; either when the typical size of microstructural features is relatively large in comparison to the macroscopic scale (e.g., in case of osteoporosis, where the pore size of bone is large) or when the actual scale of application is relatively small such that it reaches the size of microstructures (e.g., in case of miniaturization, where components are built at micrometer scale). Such scenarios call for generalized continuum theories that are able to take account of the size-dependent properties of the underlying material (Truesdell et al., 2004; Forest and Sievert, 2006). One of the most prominent generalized theories is the so-called strain-gradient theory of Mindlin (1964).

An alternative route towards size-informed models is the development of multiscale formulations that take into account the microstructure of the material explicitly. Resulting models typically reside in the realm of *generalized homogenization methods* that combine generalized macroscopic theories with homogenized properties obtained from representative volume elements (RVEs). In contrast to corresponding first-order approaches, generalized homogenization methods incorporate information on the size of underlying microstructures. One of the most famous computational generalized multiscale models has been developed by Kouznetsova, Geers, and Brekelmans (2004a), who proposed a gradient-extended macroscopic formulation for solids that is coupled to a solid Cauchy continuum at the microscopic scale. The resulting framework has been coined *second-order computational homogenization* and has since been serving as an inspiration for further developments and extensions. We refer to the contributions of, for example, Kaczmarczyk, Pearce, and Bićanić (2008) on an extension with regard to generalized boundary conditions and Lüscher, McDowell, and Bronkhorst (2010), Blanco et al. (2016), and Lopes and Pires (2022a) on alternative strategies for the up-scaling of higher-order variables and associated constraints on fluctuation fields. As could be shown by Kouznetsova, Geers, and Brekelmans (2004b), the length scale of a second-order multiscale model (which are associated with the size of an RVE) can be related to the length-scale parameters of phenomenological strain-gradient models (Mindlin, 1964; Altan and Aifantis, 1997; Polizzotto, 2003). However, the specific forms of constraints on the fluctuation fields employed by Kouznetsova, Geers, and Brekelmans (2004a) and Blanco et al. (2016) can yield mesh dependent macroscopic response (Lopes and Pires, 2022b). A more recent generalized computational-homogenization approach that combines a second-order
continuum at the macroscale and a first-order continuum at the microscale has been presented by Yvonnet, Auffray, and Monchiet (2020). The key ingredient of their approach is the definition of a microstructure-dependent body force which—in conjunction with quadratic boundary conditions at the microscopic level—should enable constant strain gradients in homogeneous RVEs. While the model achieves the intended goals for microstructures with soft matrix and stiff inclusions, it does not yield the desired effects for microstructures with reciprocal phase contrast. The realization of constant strain gradients in homogeneous RVEs is a long-standing issue in second-order homogenization methods and remains to be a part of discussions and developments to date (Yuan, Tomita, and Andou, 2008; Forest and Trinh, 2011; Barboura and Li, 2018).

In the recent years, computational homogenization has been extended to multi-field problems. In the area of deformable porous media, we refer to the first-order approaches proposed by Su, Larsson, and Runesson (2011), Pollmann et al. (2020), Ricken et al. (2022), Tu et al. (2023), Klahr et al. (2023), Liupekevicius et al. (2024), and Anonis, Mroginski, and Sánchez (2024); In the field of second-order models we would like to highlight the contribution of Thiesen, Klahr, Carniel, Blanco, et al. (2024), who propose a second-order homogenization approach for finite poro-elasticity in terms of a saddle-point principle with displacements and pressure fields as primary variables; See Thiesen, Klahr, Carniel, Holzapfel, et al. (2025) for a corresponding numerical implementation in a slightly modified setting with second-order expansion of the pressure field and first-order expansion of the deformation field. Since the mentioned models are based on the formulation of Blanco et al. (2016), it could be expected that they exhibit mesh-dependent response as discussed in Lopes and Pires (2022a). An alternative approach to the multiscale modeling of porous media with non-separated scales has recently been proposed by Maike et al. (2024) in the form of a mesh-in-element method. This approach is based on the direct embedding of microstructural meshes into a macroscopic finite element discretization.

A common property of the above mentioned works on computational homogenization of deformable porous media is that fluid flow is governed by Darcy's law. Generally, Darcy's law, which postulates a linear relationship between the pressure gradient and the fluid flux, is only applicable to fluid flow at low velocities, i.e., with a low value of the Reynolds number. When the fluid velocity increases, Darcy's law is no longer applicable and should be replaced by a more suitable model like Forchheimer's or Brinkman's law. Such considerations are particularly relevant for the homogenization of porous media, where high contrasts in phase diffusivity or the presence of large pores could give rise to steep velocity gradients; see the contribution of Monchiet (2023) and the references therein for further details. We refer to Monchiet (2023) and Monchiet (2024) for the analytical homogenization of porous media with impervious inclusions using Forchheimer's law and to the contribution of Brown et al. (2015) and Monchiet (2024) on related discussions about Brinkman's law.

In the present work, we develop a multiscale framework for the second-order computational homogenization of porous media based on a variational minimization principle. In this context, the transient microscopic problem is formulated in terms of a constrained Hu–Washizu principle, which takes account of the required continuity of the tangential and normal components of the fluid-mass flux and which is readily suitable for the consideration of Darcy, Forchheimer and Brinkman flow
<figure><img src="image_3.png" /><figcaption>**Figure 1:** Illustration of a porous medium together with its global and local fields as well as the associated boundary conditions.</figcaption></figure>

in a straightforward manner. Complementary to that, we formulate the macroscopic boundary value problem in terms of a second-order Hu-Washizu variational principle and combine both formulations via an integrated FE² model.

The outline of the work is as follows. In Section 2, we review a canonical variational minimization formulation of fluid flow through porous media that serves as the foundation for the direct numerical simulations to be conducted for validation purposes. In Section 3, we present the variational formulations of the first-order microscopic and second-order macroscopic problems and discuss in detail the localization and homogenization conditions between relevant fields of different order. The theoretical model will be discretized in time and space using implicit time integration and suitable finite elements, respectively. The predictive capabilities of the multiscale model will eventually be demonstrated by means of a sequence of two- and three-dimensional numerical simulations in Section 4.

## 2 First-order phenomenological modeling of fluid flow through porous media

To provide a theoretical basis for the present contribution, we discuss a variational formulation of fluid flow through porous media along the lines of Miehe, Mauthe, and Teichtmeister (2015). This variational formulation also serves as theoretical framework for the direct numerical simulations documented in Section 4.

The fluid-mechanical state at a material point $x$ of a fluid-saturated porous medium $B \subset \mathbb{R}^d$, where $d \in [1, 2, 3]$ is the spatial dimension, can be described by the fluid-mass density $\rho_f \in \mathbb{R} \ge 0$ and the fluid-mass flux $h \in \mathbb{R}^d$ (Coussy, 2004). Both fields, $\rho_f$ (measured in kg/m³) and $h$ (measured in kg/m² s), are functions of time $t \in \mathbb{R}$ and coupled by the balance of fluid mass (Coussy, 2004; Gurtin, Fried, and Anand, 2010)

$$ \dot{\rho}_f + \operatorname{div} h = 0 \text{ in } B. \qquad (1) $$
Furthermore, the energetic and dissipative states of the porous medium can be described by an energy-density function $\psi(\mathbf{x}, t) := \hat{\psi}(\rho_f)$ (measured in J/m³) and a dissipation-potential-density function $\phi(\mathbf{x}, t) := \hat{\phi}(\mathbf{h})$ (measured in J/m³ s), respectively, such that the energy and dissipation potentials of the porous medium $\mathcal{B}$ arise as

$$E(\rho_f) = \int_{\mathcal{B}} \psi(\mathbf{x}, t) \, dV, \qquad (2)$$

$$D(\mathbf{h}) = \int_{\mathcal{B}} \phi(\mathbf{x}, t) \, dV. \qquad (3)$$

Once complemented by a potential of external loading $P(\mathbf{h})$ according to

$$P(\mathbf{h}) := \int_{\partial \mathcal{B}_{\mu}} \bar{\mu} h \, dA, \qquad (4)$$

where $\bar{\mu}$ is the prescribed fluid potential at the boundary $\partial \mathcal{B}_{\mu}$ and $h := \mathbf{h} \cdot \mathbf{n}$ with the outward unit normal vector $\mathbf{n}$ on $\mathcal{B}$ is the out-flux of fluid mass, the energetic-dissipative state of the porous medium can be described by a time- and space-continuous potential $\Pi$ of the form

$$\tilde{\Pi}(\mathbf{h}) := \frac{d}{dt} E(\rho_f) + D(\mathbf{h}) - P(\mathbf{h}). \qquad (5)$$

In consideration of the balance of fluid mass (Eq. 1) the latter potential can be reformulated to

$$\hat{\Pi}(\mathbf{h}) = \int_{\mathcal{B}} (-\partial_{\rho_f} \hat{\psi} \operatorname{div} \mathbf{h} + \hat{\phi}(\mathbf{h})) \, dV - \int_{\partial \mathcal{B}_{\mu}} \bar{\mu} \mathbf{h} \cdot \mathbf{n} \, dA \qquad (6)$$

and the associated continuous variational principle reads (Miehe, Mauthe, and Teichtmeister, 2015)

$$\mathbf{h}^* = \arg \inf_{\mathbf{h} \in W_h} \hat{\Pi}(\mathbf{h}), \qquad (7)$$

where $W_h$ denotes the admissible space of the fluid-mass flux $\mathbf{h}$ such that

$$W_h := \{ \mathbf{h} \in H(\operatorname{div}, \mathcal{B}) \mid \mathbf{h} \cdot \mathbf{n} = \bar{h} \text{ on } \partial \mathcal{B}_h \}. \qquad (8)$$

The Euler-Lagrange equations of the above variational principle are

$$\delta \mathbf{h} : \operatorname{grad}(\partial_{\rho_f} \hat{\psi}) + \partial_{\mathbf{h}} \hat{\phi} = \mathbf{0} \text{ in } \mathcal{B} \qquad (9)$$

$$\delta \mathbf{h} : -\partial_{\rho_f} \hat{\psi} + \bar{\mu} = 0 \text{ on } \partial \mathcal{B}_{\mu} \qquad (10)$$

where we identify the fluid potential $\mu = \partial_{\rho_f} \hat{\psi}$, which is related to the pore pressure $p$ via $p := \mu / \rho_f$, where $\rho_f$ is the intrinsic density of the fluid.
## 3 Second-order multiscale modeling of fluid flow through porous media

In what follows, we present a variationally consistent second-order computational-homogenization approach to fluid flow through porous media. This approach incorporates elements of second-order multiscale methods for solids (Kouznetsova, Geers, and Brekelmans, 2004a; Blanco et al., 2016; Yvonnet, Auffray, and Monchiet, 2020; Lopes and Pires, 2022a) and minimization-based formulations of fluid flow through porous media (Miehe, Mauthe, and Teichtmeister, 2015; Polukhov and Keip, 2020; Khurshid, Polukhov, and Keip, 2024). We derive the multi-scale model in the form of continuous and incremental variational principles that couple the macroscopic and microscopic boundary-value problems via Hu-Washizu variational principles at both scales. In order to guide the reader through the subsequent derivations, we provide a brief outline of the present section and comment on the used notation.

**Outline.** We begin with the discussion of the microscopic problem of fluid flow through porous media (Section 3.1). Here, we present the second-order scale transitions of the fluid-mass flux and its spatial derivatives (Sections 3.1.1 and 3.1.2) and link the different scales via a generalized Hill-Mandel condition that enables a convenient numerical implementation of the required $H(\text{div}) \cap H(\text{grad})$-conformity of the fluctuations of the fluid-mass flux (Section 3.1.3). We then transfer the theoretical framework into an algorithmic setting by means of discretizations in time and space (Section 3.1.4) with a particular focus on the computation of algorithmically consistent macroscopic driving forces and moduli (Section 3.1.5). Once we have set the microscopic boundary-value problem theoretically and algorithmically, we turn our view to the macroscopic variational problem of fluid flow through porous media (Section 3.2). Here, we propose a mixed Hu-Washizu variational principle that relaxes the otherwise required $C^1$-continuity of the macroscopic fluid-mass flux (Section 3.2.1). Finally, we discuss the algorithmic implementation of the macroscopic formulation, which mainly follows the implementation of the microscopic problem (Section 3.2.2).

**Notation.** In the context of multiscale modeling, we describe the response of a material across different scales. In the present contribution, we present the modeling of fluid flow at macroscopic and microscopic scales. Subsequently, we describe any heterogeneity of the material at the microscale and then incorporate the associated information into the macroscopic scale considering a homogenization procedure over a representative volume element (RVE). In that connection, all quantities $\diamond$ associated with the macroscopic scale will be labeled with an overline ($\overline{\diamond}$) and all microscopic quantities will have no special labeling ($\diamond$). The same applies to mathematical operators, where we define the gradient of a macroscopic field $\overline{\diamond}$ with respect to macroscopic coordinates $\bar{x}$ as '$grad\overline{\diamond}$', the divergence of $\overline{\diamond}$ with respect to $\bar{x}$ as '$div\overline{\diamond}$', and so on. Analogously, we define the gradient and the divergence of a microscopic field $\diamond$ with respect to microscopic coordinates $x$ as '$grad\diamond$' and '$div\diamond$', respectively. Furthermore, we introduce a compact notation for the averaging operators over the volume or over the
surface of the RVE as

$$ \langle \diamond \rangle_B := \frac{1}{|\mathcal{B}|} \int_{\mathcal{B}} \diamond dV, \qquad (11) $$

$$ \langle \diamond \rangle_{\partial \mathcal{B}} := \frac{1}{|\mathcal{B}|} \int_{\partial \mathcal{B}} \diamond dA, \qquad (12) $$

where $|\mathcal{B}|$ is the volume of the RVE.

## 3.1 The microscopic variational problem of fluid flow through porous media

The microscopic problem is defined on representative volume elements (RVEs) that are attached to the Gauss points of a macroscopic boundary-value problem (Miehe, Schotte, and Schröder, 1999). Based on the assumption that the length scales of the microscopic and the macroscopic problem are not clearly separated ($\overline{\ell} \approx \ell$), we perform a second-order scale transition as will become clear from the subsequent derivations.

### 3.1.1 Second-order down-scaling of the fluid-mass flux

Assuming that the microscopic state variables fluctuate around the corresponding macroscopic state variables, we can decompose the fluid-mass flux $h$ into a microscopic contribution $\tilde{h}$ and a macroscopic contribution $\tilde{h}$ (Polukhov and Keip, 2020)

$$ h(x, t; \bar{x}) = \tilde{h}(x, t) + \tilde{h}(x, \bar{x}, t), \qquad (13) $$

where $\tilde{h}$ can be expressed in terms of the macroscopic fluid-mass flux $\tilde{h}$ at a macroscopic material point $\bar{x}$ through a Taylor expansion. In the framework of second-order homogenization, the corresponding Taylor series is expanded up to the second-order terms such that

$$ \tilde{h}(x, \bar{x}, t) = \bar{h}(\bar{x}, t) + \overline{\mathrm{grad}}\tilde{h}(\bar{x}, t) \cdot x + \frac{1}{2}\overline{\mathrm{grad}}^2\tilde{h}(\bar{x}, t) : (x \otimes x - \bar{j}), \qquad (14) $$

where $\overline{j} := \langle x \otimes x \rangle_B$ is the homogenized volume momentum of order two<fnref n="1" /> and $\overline{\mathrm{grad}}^2\diamond := \overline{\mathrm{grad}}\overline{\mathrm{grad}}\diamond$. Furthermore, the microscopic coordinates $x$ are defined such that

$$ \langle x \rangle_B = 0. \qquad (15) $$

Here we relate the microscopic and macroscopic coordinates via $x = \bar{x} - \bar{x}_{\overline{\mathcal{P}}}$, where $\bar{x}_{\overline{\mathcal{P}}}$ is the radius vector pointing to the geometric center of an RVE placed at a macroscopic material point $\overline{\mathcal{P}}$.

<footnote n="1">¹As will be discussed below, the introduction of $\overline{j}$ in Eq. 14 allows a simplified constraint of zero volume average on the fluctuation field (Blanco et al., 2016).</footnote>
Based on Eq. 14 the divergence of the microscopic fluid-mass flux follows as

$$ \mathrm{div} \, \mathbf{h} = \mathrm{div} \, \tilde{\mathbf{h}} + \overline{\mathrm{div} \, \mathbf{h}} + \overline{\mathrm{grad} \, (\overline{\mathrm{div} \, \mathbf{h}})} \cdot \mathbf{x} \qquad (16) $$

To arrive at a macroscopic formulation that is comparable to the second-grade model of Sciarra, dell'Isola, and Coussy (2007), we decompose the macroscopic gradient of the fluid-mass flux into volumetric and deviatoric parts according to

$$ \overline{\mathrm{grad} \, \mathbf{h}} = \frac{1}{d} \overline{\mathrm{div} \, \mathbf{h}} \mathbf{I} + \mathbb{P} : \overline{\mathrm{grad} \, \mathbf{h}}, \qquad (17) $$

where $[\mathbb{P}]_{ijkl} := \delta_{ik}\delta_{jl} - \frac{1}{d}\delta_{ij}\delta_{kl}$ is the fourth-order deviatoric projection tensor in terms of the Kronecker delta $\delta_{ij}$ such that $\mathbb{P} : \overline{\mathrm{grad} \, \mathbf{h}} = \mathrm{dev}(\overline{\mathrm{grad} \, \mathbf{h}})$ and $[\mathbf{I}]_{ij} := \delta_{ij}$ is the second-order identity tensor. Based on that, we can write the second gradient of the macroscopic fluid-mass flux as

$$ \overline{\mathrm{grad}}^2 \mathbf{h} = \frac{1}{d} \mathbf{I} \otimes \overline{\mathrm{grad}}(\overline{\mathrm{div}} \mathbf{h}) + \mathbb{P} : \overline{\mathrm{grad}}^2 \mathbf{h}, \qquad (18) $$

where $\mathbb{P} : \overline{\mathrm{grad}}^2 \mathbf{h} = \overline{\mathrm{grad}}[\mathrm{dev}(\overline{\mathrm{grad}} \mathbf{h})]$, and reformulate the microscopic fluid-mass flux to<fnref n="2" />

$$ \begin{aligned} \mathbf{h} = \tilde{\mathbf{h}} + \bar{\mathbf{h}} + \frac{1}{d}(\overline{\mathrm{div}} \mathbf{h})\mathbf{x} + \mathrm{dev}(\overline{\mathrm{grad}} \mathbf{h}) \cdot \mathbf{x} + \frac{1}{2d}\overline{\mathrm{grad}}(\overline{\mathrm{div}} \mathbf{h}) \cdot (\mathbf{x} \otimes \mathbf{x} - \bar{\mathbf{j}}) \\ + \frac{1}{2}\overline{\mathrm{grad}}[\mathrm{dev}(\overline{\mathrm{grad}} \mathbf{h})] : (\mathbf{x} \otimes \mathbf{x} - \bar{\mathbf{j}}). \end{aligned} \qquad (23) $$

$$ [\overline{\mathrm{grad}}^2 \mathbf{h}]_{ijk} = \frac{1}{d} \delta_{ij} [\overline{\mathrm{grad}}(\overline{\mathrm{div}} \mathbf{h})]_k + [\overline{\mathrm{grad}}[\mathrm{dev}(\overline{\mathrm{grad}} \mathbf{h})]]_{ijk}. \qquad (19) $$

Together with the symmetry condition

$$ [\overline{\mathrm{grad}}(\overline{\mathrm{div}} \mathbf{h})]_k = [\overline{\mathrm{grad}}^2 \mathbf{h}]_{ijk} \stackrel{\mathrm{sym.}}{=} [\overline{\mathrm{grad}}^2 \mathbf{h}]_{iki} \qquad (20) $$

and the identity $[\overline{\mathrm{grad}}[\mathrm{dev}(\overline{\mathrm{grad}} \mathbf{h})]]_{ijk} = 0$, it follows that

$$ [\overline{\mathrm{grad}}(\overline{\mathrm{div}} \mathbf{h})]_k = \frac{d}{d-1} [\overline{\mathrm{grad}}[\mathrm{dev}(\overline{\mathrm{grad}} \mathbf{h})]]_{iki} \qquad (21) $$

As a consequence, the terms containing $\overline{\mathrm{grad}}[\mathrm{dev}(\overline{\mathrm{grad}} \mathbf{h})]$ and $\overline{\mathrm{grad}}(\overline{\mathrm{div}} \mathbf{h})$ in Eq. 23 could be merged such that the corresponding contribution of $\overline{\mathrm{grad}}(\overline{\mathrm{div}} \mathbf{h})$ to the fluid-mass flux $\mathbf{h}$ would read

$$ \mathbf{h}_{\overline{\mathrm{grad}}(\overline{\mathrm{div}} \mathbf{h})} = \frac{1}{d+1} \overline{\mathrm{grad}}(\overline{\mathrm{div}} \mathbf{h}) \cdot (\mathbf{x} \otimes \mathbf{x} - \bar{\mathbf{j}}). \qquad (22) $$

Since in our formulation $\overline{\mathrm{grad}}[\mathrm{dev}(\overline{\mathrm{grad}} \mathbf{h})]$ and $\overline{\mathrm{grad}}(\overline{\mathrm{div}} \mathbf{h})$ appear as independent macroscopic state variables, the latter merge is not carried out so that Eq. 23 remains in its given form.

<footnote n="2">²The symmetry of $\overline{\mathrm{grad}}^2 \mathbf{h}$ makes it possible to derive a relation between $\overline{\mathrm{grad}}(\overline{\mathrm{div}} \mathbf{h})$ and $\overline{\mathrm{grad}}[\mathrm{dev}(\overline{\mathrm{grad}} \mathbf{h})]$. To facilitate this derivation, we express Eq. 18 in index notation as</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2:** Illustration of a macroscopic state variable mapped on a cubic RVE. A change from blue to red color corresponds to a change in intensity of the macroscopic contribution from low to high.</figcaption></figure>

As can be seen from the latter equation, the second-order theory enriches the microscopic fluid-mass flux $h$ with higher-order contributions from the macroscopic scale. These correspond to additional modes of fluid flow and could be compared to analogous definitions arising within gradient-extended homogenization schemes of solid-mechanical problems, where bending and torsion modes complement the classical first-order approximations (Fig. 2).

### 3.1.2 Second-order up-scaling of the fluid-mass flux

In the context of homogenization theories, macroscopic state variables are typically defined as volume averages of their microscopic counterparts. In what follows, we will present corresponding conditions and consequences thereof.

The set of variables that define the macroscopic material state can be deduced from Eq. 23 and follows as

$$ \bar{S} := \{\bar{h}, \operatorname{div}\bar{h}, \operatorname{dev}(\operatorname{grad}\bar{h}), \operatorname{grad}(\operatorname{div}\bar{h}), \operatorname{grad}[\operatorname{dev}(\operatorname{grad}\bar{h})]\}. \quad (24) $$

In this set, the macroscopic fluid-mass flux $\bar{h}$ must be equal to the volume average of the microscopic fluid-mass flux $h$ such that

$$ \bar{h} = \langle h \rangle_{\mathcal{B}}. \quad (25) $$
Consequently, we deduce that the volume average of the fluctuations of the fluid-mass flux must vanish, i.e.,

$$ \langle \tilde{h} \rangle_B = 0 \qquad (26) $$

Note that the constraint given in Eq. 26 is a consequence of the presence of the volume momentum $\bar{j}$ in Eq. 14 (Blanco et al., 2016). If we do not subtract this term, the volume average of the fluctuations $\tilde{h}$ needs to be fixed to a non-zero value such that $\langle \tilde{h} \rangle_B = -\frac{1}{2}\overline{\mathrm{grad}^2\tilde{h}} : \langle \mathbf{x} \otimes \mathbf{x} \rangle_B$. While this term does not play a role in quasi-static or steady-state problems (Lopes and Pires, 2022a), it has an impact on the transient problem under consideration because the primary fields themselves are also present in the microscopic differential equations.

In a similar manner, we relate the first-order macroscopic terms to their microscopic counterparts via

$$ \overline{\mathrm{grad} \tilde{h}} = \langle \mathrm{grad} h \rangle_B, \qquad (27) $$

$$ \overline{\mathrm{div} \tilde{h}} = \langle \mathrm{div} h \rangle_B, \qquad (28) $$

so that we obtain the constraints

$$ \langle \mathrm{grad} \tilde{h} \rangle_B = 0 \quad \Leftrightarrow \quad \langle \tilde{h} \otimes n \rangle_{\partial B} = 0 \qquad (29) $$

$$ \langle \mathrm{div} \tilde{h} \rangle_B = 0 \quad \Leftrightarrow \quad \langle \tilde{h} \cdot n \rangle_{\partial B} = 0 \qquad (30) $$

where the surface integrals with outward unit normal vector $n$ are obtained by making use of Gauss' law. The given constraints are a priori fulfilled in case of periodic fluctuations of the fluid-mass flux.

The second-order macroscopic terms can be determined from Eqs. 13 and 14 as (Kouznetsova, Geers, and Brekelmans, 2004a; Luscher, McDowell, and Bronkhorst, 2010; Blanco et al., 2016; Lopes and Pires, 2022a)

$$ \overline{\mathrm{grad}^2 \tilde{h}} = \langle \mathrm{grad} \, h \otimes \mathbf{x} \cdot \bar{\mathbf{j}}^{-1} \rangle_B = \langle \mathrm{grad} \, h \otimes \mathbf{x} \rangle_B \cdot \bar{\mathbf{j}}^{-1}, \qquad (31) $$

which results in the constraint

$$ \langle \mathrm{grad} \, \tilde{h} \otimes \mathbf{x} \rangle_B = 0 \qquad (32) $$

This volume integral can be converted into a surface integral by using the product rule and Gauss' law so that we obtain

$$ \langle \tilde{h} \otimes n \otimes \mathbf{x} \rangle_{\partial B} - \langle \tilde{h} \rangle_B \otimes \mathbf{1} = 0. \qquad (33) $$
In the latter equation, we could exploit Eq. 26 so that we finally get

$$ \langle \tilde{h} \otimes n \otimes x \rangle_{\partial \mathcal{B}} = 0. \qquad (34) $$

As was shown by Lopes and Pires (2022a), this constraint can be further simplified in case of hyperrectangular RVEs to

$$ \int_{\partial \mathcal{B}_i^+} \tilde{h} \, dA = 0 \quad \text{with } i = 1, \dots, d, \qquad (35) $$

where $d$ denotes the spatial dimension and the superscript $\mathcal{B}_i^+$ signifies independent boundaries of a periodic RVE.

### 3.1.3 Thermodynamically consistent second-order scale transition

To guarantee a thermodynamically consistent scale transition, we define a generalized Hill-Mandel condition of the form (Fig. 3)

$$ \bar{\pi}(\bar{S}) = \inf_{\tilde{h} \in W_{\tilde{h}}} \sup_{\bar{\lambda} \in W_{\bar{\lambda}}} \left\{ \frac{1}{|\mathcal{B}|} \int_{\mathcal{B}} \left( d_t \widehat{\psi}(\rho_f) + \widehat{\phi}(h) \right) dV + \bar{\lambda} \cdot \frac{1}{|\mathcal{B}|} \int_{\mathcal{B}} \tilde{h} \, dV \right. \\ \left. + \bar{\eta} := \frac{1}{|\mathcal{B}|} \int_{\mathcal{B}} \mathrm{grad} \, \tilde{h} \otimes x \, dV - \frac{1}{2\kappa_{\lambda}} \|\bar{\lambda}\|^2 - \frac{1}{2\kappa_{\eta}} \|\bar{\eta}\|^2 \right\}, \qquad (36) $$

where $d_t \diamond$ denotes the time rate of $\diamond$, $\bar{\lambda} \in \mathbb{R}^d$ and $\bar{\eta} \in \mathbb{R}^{d \times d \times d}$ with the squared norms $\|\bar{\lambda}\|^2 := \bar{\lambda} \cdot \bar{\lambda}$ and $\|\bar{\eta}\|^2 := \bar{\eta} \cdot \bar{\eta}$ are perturbed$^3$ macroscopic Lagrange multipliers$^4$ that are energetically dual to the averaged integral constraints of the fluctuation field (Eqs. 26 and 32), the violation of which could be controlled by adjustment of the penalty parameters $\kappa_{\lambda}$ and $\kappa_{\eta}$, the operator $\cdot\cdot\cdot$ denotes a triple contraction, the fluid-mass density $\rho_f$ is related to the fluid-mass flux $h$ through the balance of mass (Eq. 1), and $W_{\diamond}$ are admissible function spaces of fields $\diamond$ according to$^5$

$$ W_{\tilde{h}} := \{ \tilde{h} \in H(\mathrm{div}, \mathcal{B}) \cap [H(\mathrm{grad}, \mathcal{B})]^d \mid \tilde{h}^+ = \tilde{h}^- \text{ on } \partial \mathcal{B} = \partial \mathcal{B}^+ \cup \partial \mathcal{B}^- \}, \qquad (38) $$

$$ W_{\bar{\lambda}} := \{ \bar{\lambda} \in [L^2(\bar{\mathcal{B}})]^d \}, \qquad (39) $$

$$ W_{\bar{\eta}} := \{ \bar{\eta} \in [L^2(\bar{\mathcal{B}})]^{d \times d \times d} \}. \qquad (40) $$

$$ \bar{\eta}^{\diamond} := \langle \mathrm{grad} \, \tilde{h} \otimes x \cdot \bar{j}^{-1} \rangle_{\mathcal{B}} = \bar{\eta}^{\diamond} \cdot \bar{j}^{-1} := \langle \mathrm{grad} \, \tilde{h} \otimes x \rangle_{\mathcal{B}}, \qquad (37) $$

<footnote>$^3$The notion of a perturbed Lagrangian method is discussed in the book of Wriggers (2008, §10.2).</footnote>

<footnote>$^4$The constraint associated with the Lagrange multiplier $\eta$ is equivalent to the one proposed by Lopes and Pires (2022a) since</footnote>

<footnote>where $\bar{\eta}^{\diamond}$ represents the Lagrange multiplier employed by Lopes and Pires (2022a).</footnote>

<footnote>$^5$The function spaces $W_{\tilde{h}}$ and $W_{\bar{\eta}}$ are associated with the macroscopic boundary value problem.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3:** Multiscale view on a heterogeneous porous medium. In the context of multiscale modeling, the macroscopic response at a macroscopic material point $\bar{x} \in \bar{\mathcal{B}}$ is determined via homogenization of the microscopic response over a representative volume element (RVE) defined on the domain $\mathcal{B}$. In the realm of second-order homogenization schemes, the classical concept of scale separation is relaxed such that the characteristic length scales of the RVE ($\ell$) and the macroscopic structure ($\bar{\ell}$) do not have to be clearly separated.</figcaption></figure>

The variational formulation given in Eq. 36 can be exploited by means of the fundamental lemma of variational calculus to derive the weak form

$$ \delta\bar{\pi}(\bar{S}) = \frac{1}{|\mathcal{B}|} \int_{\mathcal{B}} \left( \operatorname{div} \delta\tilde{h}(-\partial_{\rho_f}\hat{\psi}) + \delta\tilde{h} \cdot (\partial_h\hat{\phi} + \bar{\lambda}) + \operatorname{grad} \delta\tilde{h} : (\bar{\eta} \cdot \mathbf{x}) \right) \mathrm{d}V \\ + \delta\bar{\lambda} \cdot \left( \frac{1}{|\mathcal{B}|} \int_{\mathcal{B}} \tilde{h} \mathrm{d}V - \frac{1}{\kappa_{\lambda}}\bar{\lambda} \right) + \delta\bar{\eta} : \left( \frac{1}{|\mathcal{B}|} \int_{\mathcal{B}} \operatorname{grad} \tilde{h} \otimes \mathbf{x} \mathrm{d}V - \frac{1}{\kappa_{\eta}}\bar{\eta} \right) = 0, \quad (41) $$

so that we obtain the Euler-Lagrange equations of the microscopic boundary-value problem

$$ \delta\tilde{h} : \operatorname{grad}(\partial_{\rho_f}\hat{\psi}) + \partial_h\hat{\phi} + \bar{\lambda} - \operatorname{div}(\bar{\eta} \cdot \mathbf{x}) = \mathbf{0} \quad \text{in } \mathcal{B} \quad (42) $$

$$ \delta\tilde{h} : (-\partial_{\rho_f}\hat{\psi}\mathbf{1} + \bar{\eta} \cdot \mathbf{x}) \cdot \mathbf{n} = \mathbf{0} \quad \text{on } \partial\mathcal{B} \quad (43) $$

$$ \delta\bar{\lambda} : \langle \tilde{h} \rangle_{\mathcal{B}} - \frac{1}{\kappa_{\lambda}}\bar{\lambda} = \mathbf{0} \quad \text{in } \mathcal{B} \quad (44) $$

$$ \delta\bar{\eta} : \langle \operatorname{grad} \tilde{h} \otimes \mathbf{x} \rangle_{\mathcal{B}} - \frac{1}{\kappa_{\eta}}\bar{\eta} = \mathbf{0} \quad \text{in } \mathcal{B} \quad (45) $$

Note that we solve the microscopic problem for a given macroscopic state and hence do not consider variations of the macroscopic state variables in the microscopic boundary-value problem. This corresponds to solving an FE² problem in a staggered manner. We refer to, for example, Tan, Raju, and Lee (2020), Lange, Hütter, and Kiefer (2021), Hesch, Schmidt, and Schuß (2024), and Yang et al. (2024) for alternative approaches.
### 3.1.4 Algorithmic implementation of the microscopic variational principle

According to the variational formulation given in Eq. 36 the macroscopic state is defined at an equilibrium state of the microscopic boundary-value problem. In order to obtain the latter, we now discretize variational formulation in space and time.

**Discretization in time.** To obtain a time-discrete incremental potential corresponding to Eq. 36, we follow Miehe, Mauthe, and Teichtmeister (2015) and employ the implicit Euler scheme within a time interval $[t_n, t_{n+1}]$ with a time-step size $\tau = t_{n+1} - t_n$ such that

$$ \bar{\pi}^{\tau}(\bar{S}) = \inf_{\tilde{h} \in W_{\tilde{h}}} \sup_{\bar{\lambda} \in W_{\bar{\lambda}}} \sup_{\bar{\eta} \in W_{\bar{\eta}}} \left\{ \frac{1}{|\mathcal{B}|} \int_{\mathcal{B}} \left( \hat{\psi}(\rho_f) - \hat{\psi}(\rho_f(t_n)) + \tau \hat{\phi}(\mathbf{h}) \right) dV + \bar{\lambda}^{\tau} \cdot \frac{1}{|\mathcal{B}|} \int_{\mathcal{B}} \tilde{h} \, dV \right. \\ \left. + \bar{\eta}^{\tau} \cdot \frac{1}{|\mathcal{B}|} \int_{\mathcal{B}} \mathrm{grad} \, \tilde{h} \otimes \mathbf{x} \, dV - \frac{1}{2\kappa_{\lambda}^{\tau}} \|\bar{\lambda}^{\tau}\|^2 - \frac{1}{2\kappa_{\eta}^{\tau}} \|\bar{\eta}^{\tau}\|^2 \right\}, \quad (46) $$

where all quantities without an index refer to time $t_{n+1}$. Furthermore, we have introduced $\bar{\lambda}^{\tau} := \tau \bar{\lambda}$, $\bar{\eta}^{\tau} := \tau \bar{\eta}$, $\kappa_{\lambda/\eta}^{\tau} := \tau \kappa_{\lambda/\eta}$ and $\rho_f = \rho_f(t_n) - \tau \mathrm{div} \, \mathbf{h}$.

As we concluded in Eq. 38, the field $\tilde{h}$ is embedded into both $H(\mathrm{div}, \mathcal{B})$ and $H(\mathrm{grad}, \mathcal{B})$ spaces. Therefore, we cannot discretize the RVEs using standard $H(\mathrm{div}, \mathcal{B})$ finite elements, such as Raviart–Thomas or Brezzi–Douglas–Marini finite elements, which ensure the continuity of normal components but not of tangential components of discretized fields across finite elements (Brezzi and Fortin, 1991). This consideration would lead to a non-conforming discretization and non-defined gradients. While Huang–Zhang and Mardal–Tai–Winther finite elements ensure continuity of both normal and tangential components of a field, they come with challenges with respect to the implementation of the mapping from the parametric space to the actual mesh (Aznaran, Farrell, and Kirby, 2022; Sky et al., 2023). Hence, we extend the incremental variational problem provided in Eq. 46 to a Hu–Washizu formulation<fnref n="6" /> that allows for an implementation using mixed finite-element methods (Teichtmeister, Mauthe, and Miehe, 2019), according to

$$ \bar{\pi}_{\mathrm{mix}}^{\tau}(\bar{S}) = \inf_{\tilde{h} \in W_{\tilde{h}}} \inf_{\theta \in W_{\theta}} \sup_{\mu \in W_{\mu}} \sup_{\bar{\lambda} \in W_{\bar{\lambda}}} \sup_{\bar{\eta} \in W_{\bar{\eta}}} \left\{ \frac{1}{|\mathcal{B}|} \int_{\mathcal{B}} \left( \hat{\psi}(\check{\rho}_f) - \hat{\psi}(\check{\rho}_f(t_n)) + \tau \hat{\phi}(\mathbf{h}) \right. \right. \\ \left. \left. + \tau \mu(\theta - \mathrm{div} \, \mathbf{h}) - \frac{\tau}{2\kappa_{\mu}} \mu^2 \right) dV + \bar{\lambda}^{\tau} \cdot \frac{1}{|\mathcal{B}|} \int_{\mathcal{B}} \tilde{h} \, dV \right. \\ \left. + \bar{\eta}^{\tau} \cdot \frac{1}{|\mathcal{B}|} \int_{\mathcal{B}} \mathrm{grad} \, \tilde{h} \otimes \mathbf{x} \, dV - \frac{1}{2\kappa_{\lambda}^{\tau}} \|\bar{\lambda}^{\tau}\|^2 \right. \\ \left. - \frac{1}{2\kappa_{\eta}^{\tau}} \|\bar{\eta}^{\tau}\|^2 \right\}, \quad (47) $$

where $\check{\rho}_f := \check{\rho}_f(t_n) - \tau \theta$ is an auxiliary fluid-mass density with initial condition $\check{\rho}_f(t_0) = \rho_f(t_0)$, $\mu$ is a

<footnote n="6">⁶For mixed finite-element formulations of Stokes-Brinkman-Darcy flow on the basis of $H(\mathrm{div})$-conforming solution spaces for the velocity we refer to Carvalho, Devloo, and Gomes (2024).</footnote>
perturbed Lagrange multiplier associated with the constraint $\theta \equiv \operatorname{div} \mathbf{h}$, the violation of which could be controlled by the adjustment of the parameter $\kappa_{\mu} \in \mathbb{R} > 0$. The two additional fields $\theta$ and $\mu$ are embedded into the spaces

$$ W_{\theta} := \{\theta \in L^2(\mathcal{B})\}, \qquad (48) $$

$$ W_{\mu} := \{\mu \in L^2(\mathcal{B})\}. \qquad (49) $$

Again, by making use of the fundamental lemma of the variational calculus

$$ \begin{aligned} \delta \bar{\pi}_{\text{mix}}^{\tau}(\bar{S}) = \frac{1}{|\mathcal{B}|} \int_{\mathcal{B}} \left( \operatorname{div} \delta \tilde{h}(-\tau\mu) + \delta \tilde{h} \cdot (\tau \partial_h \hat{\phi} + \bar{\lambda}^{\tau}) + \operatorname{grad} \delta \tilde{h} : (\bar{\eta}^{\tau} \cdot \mathbf{x}) + \delta \theta(\tau\mu - \tau \partial_{\check{\rho}_f} \hat{\psi}) \right. \\ \left. + \delta \mu(\tau\theta - \tau \operatorname{div} \mathbf{h} - \frac{\tau}{\kappa_{\mu}}\mu) \right) dV + \delta \bar{\lambda}^{\tau} \cdot \left( \frac{1}{|\mathcal{B}|} \int_{\mathcal{B}} \tilde{h} \, dV - \frac{1}{\kappa_{\lambda}^{\tau}} \bar{\lambda}^{\tau} \right) \\ + \delta \bar{\eta}^{\tau} := \left( \frac{1}{|\mathcal{B}|} \int_{\mathcal{B}} \operatorname{grad} \tilde{h} \otimes \mathbf{x} \, dV - \frac{1}{\kappa_{\eta}^{\tau}} \bar{\eta}^{\tau} \right) = 0, \end{aligned} \qquad (50) $$

we obtain the Euler–Lagrange equations of the extended formulation

$$ \delta \tilde{h} : \operatorname{grad}(\tau\mu) + \tau \partial_h \hat{\phi} + \bar{\lambda}^{\tau} - \operatorname{div}(\bar{\eta}^{\tau} \cdot \mathbf{x}) = \mathbf{0} \text{ in } \mathcal{B} \qquad (51) $$

$$ \delta \tilde{h} : [[-\tau\mu\mathbf{1} + \bar{\eta}^{\tau} \cdot \mathbf{x}]] \cdot \mathbf{n} = \mathbf{0} \text{ on } \partial\mathcal{B} \qquad (52) $$

$$ \delta \theta : \tau\mu - \tau \partial_{\check{\rho}_f} \hat{\psi} = 0 \text{ in } \mathcal{B} \qquad (53) $$

$$ \delta \mu : \tau(\theta - \operatorname{div} \mathbf{h}) - \frac{\tau}{\kappa_{\mu}}\mu = 0 \text{ in } \mathcal{B} \qquad (54) $$

$$ \delta \bar{\lambda}^{\tau} : \langle \tilde{h} \rangle_{\mathcal{B}} - \frac{1}{\kappa_{\lambda}^{\tau}} \bar{\lambda}^{\tau} = \mathbf{0} \text{ in } \mathcal{B} \qquad (55) $$

$$ \delta \bar{\eta}^{\tau} : \langle \operatorname{grad} \tilde{h} \otimes \mathbf{x} \rangle_{\mathcal{B}} - \frac{1}{\kappa_{\eta}^{\tau}} \bar{\eta}^{\tau} = \mathbf{0} \text{ in } \mathcal{B} \qquad (56) $$

where we identify the Lagrange multiplier $\mu$ as the microscopic fluid potential defined via the relation $\mu = \partial_{\check{\rho}_f} \hat{\psi}$ in a weak sense.

**Discretization in space.** We now implement the above time-discrete variational formulation into a mixed finite-element setting. To make the implementation and its derivation compact, we employ a classical vector-matrix notation for tensors of different order in what follows.

Let us begin with the spatial discretization of the fluctuations of the fluid-mass flux $\tilde{h}$ and its spatial
derivatives. These are defined as

$$ \tilde{h} \approx \sum_{I=1}^{n_h^h} N_h^I d^I := \mathbf{N}_h d^e, \qquad (57) $$

$$ \mathrm{div} \, \tilde{h} \approx \sum_{I=1}^{n_h^h} \mathbf{d}^I \cdot \mathrm{grad} \, N_h^I := \mathbf{B}_h^{\mathrm{div}} \mathbf{d}^e, \qquad (58) $$

$$ \mathrm{grad} \, \tilde{h} \approx \sum_{I=1}^{n_h^h} \mathbf{d}^I \otimes \mathrm{grad} \, N_h^I := \mathbf{B}_h^{\mathrm{grad}} \mathbf{d}^e, \qquad (59) $$

where $N_h^I$ denotes the Lagrangian interpolation functions associated with node $I$ and field $\tilde{h}$, $n_h^h$ is the number of nodes per finite element associated with the field $\tilde{h}$, and $\mathbf{d}^I$ and $\mathbf{d}^e$ are the finite-element arrays containing the corresponding degrees of freedom per node and per finite element $e$, respectively. In the following numerical simulations, we select quadratic interpolation functions $N_h^I$ so that the order of the discrete $\tilde{h}$ coincides with the order of the second-order expansion of $\tilde{h}$ varying quadratically in $x$ (Eq. 14). To be specific, our implementation is based on nine-noded quadrilateral elements (Q9) in two spatial dimensions and ten-noded tetrahedral elements (TET10) in three spatial dimensions. For the numerical integration, we use nine Gauss points in two dimensions and four Gauss points in three dimensions. Since the fields $\theta$ and $\mu$ are defined in $L^2(\mathcal{B})$, we discretize them employing element-wise discontinuous linear interpolation functions of the form

$$ \theta = \sum_{I=1}^{n_h^\theta} N_\theta^I \alpha_\theta^I := \mathbf{N}_\theta \alpha_\theta^e, \qquad (60) $$

$$ \mu = \sum_{I=1}^{n_h^\mu} N_\mu^I \alpha_\mu^I := \mathbf{N}_\mu \alpha_\mu^e, \qquad (61) $$

where the number of internal nodes $n_h^{\theta/\mu}$ is 3 in two spatial dimensions and 4 in three spatial dimensions. The arrays $\alpha_{\theta/\mu}^I$ and $\alpha_{\theta/\mu}^e$ contain the internal degrees of freedom associated with $\theta$ and $\mu$ per internal node $I$ and finite element $e$, respectively. The corresponding interpolation functions $N_{\theta/\mu}^I$ could either be defined in terms of reference coordinates $\xi_i \in \mathcal{A}^e$ with the reference element $\mathcal{A}^e$ such that $N_{\theta/\mu}^I = 1 + \sum_{i=1}^d \xi_i$ or in terms of the physical coordinates $x \in \mathcal{B}$ such that $N_{\theta/\mu}^I = 1 + \sum_{i=1}^d x_i$ (Brink and Stein, 1996; Auricchio et al., 2017).

We again note that $\bar{\lambda}^\tau$ and $\bar{\eta}^\tau$ are not interpolated at the microstructure level, since they are constant over the RVE. Both variables will be determined at the Gauss points of the macroscopic boundary-value problem.
In consideration of the above, we arrive at the space-time discrete version of Eq. 50

$$ \delta \bar{\pi}_{\text{mix}}^{\tau,h} = \frac{1}{|\mathcal{B}|} \mathbf{A}_{e=1}^{n_e} \begin{bmatrix} \delta \mathbf{d}^e \\ \delta \alpha^e \end{bmatrix} \cdot \left[ \begin{bmatrix} \int_{\mathcal{B}^e} \partial \mathbf{d}^e \hat{\pi}_{\text{mix}}^{\tau,h} \, dV \\ \int_{\mathcal{B}^e} \partial \alpha^e \hat{\pi}_{\text{mix}}^{\tau,h} \, dV \end{bmatrix} + \begin{bmatrix} \delta \bar{\lambda}^\tau \\ \delta \bar{\eta}^\tau \end{bmatrix} \cdot \begin{bmatrix} \frac{1}{|\mathcal{B}|} \sum_{e=1}^{n_e} \int_{\mathcal{B}^e} \partial_{\bar{\lambda}^\tau} \hat{\pi}_{\text{mix}}^{\tau,h} \, dV - \frac{1}{\kappa_\lambda^\tau} \bar{\lambda}^\tau \\ \frac{1}{|\mathcal{B}|} \sum_{e=1}^{n_e} \int_{\mathcal{B}^e} \partial_{\bar{\eta}^\tau} \hat{\pi}_{\text{mix}}^{\tau,h} \, dV - \frac{1}{\kappa_\eta^\tau} \bar{\eta}^\tau \end{bmatrix} \right] = 0, \quad (62) $$

where $\mathcal{B}^e$ denotes the domain of a finite element $e$, such that $\mathcal{B} \approx \mathcal{B}^h = \cup_{e=1}^{n_e} \mathcal{B}^e$, $\alpha^e$ is the array of the internal degrees of freedom, such that $\alpha^e := [\alpha_\theta^e, \alpha_\mu^e]^T$. $\hat{\pi}_{\text{mix}}^{\tau,h}$ is a time-discrete extended potential density defined by

$$ \hat{\pi}_{\text{mix}}^{\tau,h}(\check{\rho}_f, \tilde{h}, \text{grad } \tilde{h}, \mu, \bar{\lambda}^{\tau}, \bar{\eta}^{\tau}) := \hat{\psi}(\check{\rho}_f) - \hat{\psi}(\check{\rho}_f(t_n)) + \tau \hat{\phi}(h) + \tau \mu(\theta - \text{div } h) - \frac{\tau}{\kappa_\mu} \mu^2 \\ + \bar{\lambda}^{\tau} \cdot \tilde{h} + \bar{\eta}^{\tau} \quad \therefore (\text{grad } \tilde{h} \otimes x), \quad (63) $$

and $\partial_{\mathbf{d}^e} \hat{\pi}_{\text{mix}}^{\tau,h}$ and $\partial_{\alpha^e} \hat{\pi}_{\text{mix}}^{\tau,h}$ are corresponding derivatives such that

$$ \partial_{\mathbf{d}^e} \hat{\pi}_{\text{mix}}^{\tau,h} = \begin{bmatrix} \mathbf{N}_h^T (\tau \partial_h \hat{\phi} + \bar{\lambda}^\tau) \\ (\mathbf{B}_h^{\text{div}})^T (-\tau \mu) \\ (\mathbf{B}_h^{\text{grad}})^T (\bar{\eta}^\tau \cdot x) \end{bmatrix}, \quad (64) $$

$$ \partial_{\alpha^e} \hat{\pi}_{\text{mix}}^{\tau,h} = \begin{bmatrix} \mathbf{N}_\theta^T (\tau \mu - \tau \partial_{\check{\rho}_f} \hat{\psi}) \\ \mathbf{N}_\mu^T (\tau \theta - \tau \text{div } h - \frac{\tau}{\kappa_\mu} \mu) \end{bmatrix}. \quad (65) $$

The nonlinearity of Eq. 62 requires an iterative solution, for which we employ the Newton-Raphson method. The linearization of Eq. 62 yields

$$ \Delta \delta \bar{\pi}_{\text{mix}}^{\tau,h} = \frac{1}{|\mathcal{B}|} \sum_{e=1}^{n_e} \begin{bmatrix} \delta \mathbf{d}^e \\ \delta \alpha^e \end{bmatrix} \cdot \int_{\mathcal{B}^e} \begin{bmatrix} \partial_{\mathbf{d}^e}^2 \hat{\pi}_{\text{mix}}^{\tau,h} & \partial_{\mathbf{d}^e}^2 \hat{\pi}_{\text{mix}}^{\tau,h} & \partial_{\mathbf{d}^e}^2 \hat{\pi}_{\text{mix}}^{\tau,h} \\ \partial_{\alpha^e \mathbf{d}^e}^2 \hat{\pi}_{\text{mix}}^{\tau,h} & \partial_{\alpha^e \alpha^e}^2 \hat{\pi}_{\text{mix}}^{\tau,h} & \mathbf{0} \end{bmatrix} \begin{bmatrix} \Delta \mathbf{d}^e \\ \Delta \alpha^e \\ \Delta \bar{\lambda}^\tau \\ \Delta \bar{\eta}^\tau \end{bmatrix} \, dV. \quad (66) $$

$$ + \begin{bmatrix} \delta \bar{\lambda}^\tau \\ \delta \bar{\eta}^\tau \end{bmatrix} \cdot \begin{bmatrix} \frac{1}{|\mathcal{B}|} \mathbf{A}_{e=1}^{n_e} \int_{\mathcal{B}^e} \partial_{\bar{\lambda}^\tau}^2 \hat{\pi}_{\text{mix}}^{\tau,h} \, dV & -\frac{1}{\kappa_\lambda^\tau} \mathbf{I} & \mathbf{0} \\ \frac{1}{|\mathcal{B}|} \mathbf{A}_{e=1}^{n_e} \int_{\mathcal{B}^e} \partial_{\bar{\eta}^\tau}^2 \hat{\pi}_{\text{mix}}^{\tau,h} \, dV & \mathbf{0} & -\frac{1}{\kappa_\eta^\tau} \mathbf{I} \end{bmatrix} \cdot \begin{bmatrix} \Delta \mathbf{d} \\ \Delta \bar{\lambda}^\tau \\ \Delta \bar{\eta}^\tau \end{bmatrix} $$

where $\mathbf{A}_{e=1}^{n_e}$ is the assembly operator over the $n_e$ finite elements of an RVE, $\Delta \mathbf{d} := \mathbf{A}_{e=1}^{n_e} \Delta \mathbf{d}^e$ denotes a vector containing the incremental solutions associated with all degrees of freedom of the RVE as a
result of a single microscopic iteration, and furthermore

$$ \partial_{\mathbf{d}^e}^2 \hat{\pi}_{\text{mix}}^{\tau,h} = \mathbf{N}_h^T (\tau \partial_{hh} \hat{\phi}_{\text{mix}}^{\tau,h}) \mathbf{N}_h \quad (67) $$

$$ \partial_{\mathbf{d}^e}^2 \alpha^e \hat{\pi}_{\text{mix}}^{\tau,h} = -\tau (\mathbf{B}_h^{\text{div}})^T \mathbf{N}_\mu, \quad (68) $$

$$ \partial_{\alpha^e}^2 \hat{\pi}_{\text{mix}}^{\tau,h} = \begin{bmatrix} \tau^2 \partial_{\tilde{\rho}_f \tilde{\rho}_f}^2 \hat{\psi} \mathbf{N}_\theta^T \mathbf{N}_\theta & \tau \mathbf{N}_\theta^T \mathbf{N}_\mu \\ \tau \mathbf{N}_\mu^T \mathbf{N}_\theta & -\frac{\tau}{\kappa_\mu} \mathbf{N}_\mu^T \mathbf{N}_\mu \end{bmatrix}, \quad (69) $$

$$ \partial_{\mathbf{d}^e}^2 \hat{\pi}_{\text{mix}}^{\tau,h} = \mathbf{N}_h^T(I), \quad (70) $$

$$ \partial_{\mathbf{d}^e}^2 \hat{\eta}^{\tau} \hat{\pi}_{\text{mix}}^{\tau,h} = (\mathbf{B}_h^{\text{grad}})^T (\mathbf{I} \cdot \mathbf{x}), \quad (71) $$

where $I$ and $\mathbf{I}$ are a second- and a sixth-order identity tensor implemented in vectorial forms, and '$\diamond$' designates a contraction of the third index of a tensor $\diamond$ with the first index of a tensor $\diamond$.<fnref n="7" />

Since the solution space of $\alpha^e$ is not continuous across finite elements, we can condense it out of the global system of equations at element level by using (Miehe, 1994)

$$ \Delta \alpha^e = - \int_{\mathcal{B}^e} \partial_{\alpha^e}^2 \hat{\pi}_{\text{mix}}^{\tau,h} \, dV^{-1} \left( \int_{\mathcal{B}^e} \partial_{\alpha^e} \hat{\pi}_{\text{mix}}^{\tau,h} \, dV + \int_{\mathcal{B}^e} \partial_{\alpha^e}^2 \hat{\pi}_{\text{mix}}^{\tau,h} \, dV \cdot \Delta \mathbf{d}^e \right). \quad (77) $$

As a result, we obtain a reduced version of the global system of equations given by

$$ \delta \hat{\pi}_{\text{mix}}^{\tau,h,\star} + \Delta \delta \hat{\pi}_{\text{mix}}^{\tau,h,\star} = \begin{bmatrix} \delta \mathbf{d} \\ \delta \bar{\lambda}^{\tau} \\ \delta \bar{\eta}^{\tau} \end{bmatrix} \cdot \left\{ \begin{bmatrix} \mathbf{R}_{\mathbf{d}} \\ \mathbf{R}_{\bar{\lambda}^{\tau}} \\ \mathbf{R}_{\bar{\eta}^{\tau}} \end{bmatrix} + \begin{bmatrix} \mathbf{K}_{\mathbf{d}\mathbf{d}} & \mathbf{K}_{\mathbf{d}\bar{\lambda}^{\tau}} & \mathbf{K}_{\mathbf{d}\bar{\eta}^{\tau}} \\ \mathbf{K}_{\bar{\lambda}^{\tau}\mathbf{d}} & \mathbf{K}_{\bar{\lambda}^{\tau}\bar{\lambda}^{\tau}} & \mathbf{K}_{\bar{\eta}^{\tau}\bar{\eta}^{\tau}} \\ \mathbf{K}_{\bar{\eta}^{\tau}\mathbf{d}} & \mathbf{K}_{\bar{\eta}^{\tau}\bar{\eta}^{\tau}} & \mathbf{K}_{\bar{\eta}^{\tau}\bar{\eta}^{\tau}} \end{bmatrix} \cdot \begin{bmatrix} \Delta \mathbf{d} \\ \Delta \bar{\lambda}^{\tau} \\ \Delta \bar{\eta}^{\tau} \end{bmatrix} \right\} = 0 \quad (78) $$

with the global arrays

$$ \mathbf{R}_{\mathbf{d}} := \frac{1}{|\mathcal{B}|} \mathbf{A}_{e=1}^{n_e} \left( \int_{\mathcal{B}^e} \partial_{\mathbf{d}^e} \hat{\pi}_{\text{mix}}^{\tau,h} \, dV - \int_{\mathcal{B}^e} \partial_{\mathbf{d}^e}^2 \alpha^e \hat{\pi}_{\text{mix}}^{\tau,h} \, dV \left\{ \int_{\mathcal{B}^e} \partial_{\alpha^e}^2 \hat{\pi}_{\text{mix}}^{\tau,h} \, dV \right\}^{-1} \int_{\mathcal{B}^e} \partial_{\alpha^e} \hat{\pi}_{\text{mix}}^{\tau,h} \, dV \right), \quad (79) $$

$$ \bar{\lambda} \triangleq [\bar{\lambda}_1^\tau, \bar{\lambda}_2^\tau]^T \quad (72) $$

$$ \bar{\eta}^{\tau} \triangleq [\bar{\eta}_{111}^{\tau}, \bar{\eta}_{221}^{\tau}, \bar{\eta}_{121}^{\tau}, \bar{\eta}_{211}^{\tau}, \bar{\eta}_{112}^{\tau}, \bar{\eta}_{222}^{\tau}, \bar{\eta}_{122}^{\tau}, \bar{\eta}_{212}^{\tau}]^{\mathrm{T}}, \quad (73) $$

$$ \bar{\eta}^{\tau} \cdot \mathbf{x} \triangleq [\bar{\eta}_{111}^{\tau} x_1, \bar{\eta}_{221}^{\tau} x_1, \bar{\eta}_{121}^{\tau} x_1, \bar{\eta}_{211}^{\tau} x_1, \bar{\eta}_{112}^{\tau} x_2, \bar{\eta}_{222}^{\tau} x_2, \bar{\eta}_{122}^{\tau} x_2, \bar{\eta}_{212}^{\tau} x_2]^{\mathrm{T}}. \quad (74) $$

$$ \mathbf{I} \triangleq \operatorname{diag}\{1, 1, 0, 0\} \quad (75) $$

$$ (\mathbf{I}^{\cdot 3} \mathbf{x}) = \partial_{\bar{\eta}^{\tau}} (\bar{\eta}^{\tau} \cdot \mathbf{x}) \triangleq \operatorname{diag}\{x_1, x_1, x_1, x_1, x_2, x_2, x_2, x_2\} \quad (76) $$

Our implementation in three-dimensions is based on an analogous extension of the arrays.

<footnote n="7">⁷For clarity, we provide the vector notations of some relevant terms in two spatial dimensions below:</footnote>
$$ \mathbf{R}_{\bar{\lambda}^{\tau}} := \frac{1}{|\mathcal{B}|} \sum_{e=1}^{n_c} \int_{\mathcal{B}^e} \partial_{\bar{\lambda}^{\tau}} \hat{\pi}_{\text{mix}}^{\tau,h} dV - \frac{1}{\kappa_{\bar{\lambda}}^{\tau}} \bar{\lambda}^{\tau}, \qquad (80) $$

$$ \mathbf{R}_{\bar{\eta}^{\tau}} := \frac{1}{|\mathcal{B}|} \sum_{e=1}^{n_c} \int_{\mathcal{B}^e} \partial_{\bar{\eta}^{\tau}} \hat{\pi}_{\text{mix}}^{\tau,h} dV - \frac{1}{\kappa_{\bar{\eta}}^{\tau}} \bar{\eta}^{\tau}, \qquad (81) $$

$$ \mathbf{K}_{dd} := \frac{1}{|\mathcal{B}|} \mathbf{A}_{e=1}^{n_c} \left( \int_{\mathcal{B}^e} \partial_{\mathbf{d}^e}^2 \hat{\pi}_{\text{mix}}^{\tau,h} dV - \int_{\mathcal{B}^e} \partial_{\mathbf{d}^e}^2 \alpha^e \hat{\pi}_{\text{mix}}^{\tau,h} dV \right) \left\{ \int_{\mathcal{B}^e} \partial_{\alpha^e}^2 \alpha^e \hat{\pi}_{\text{mix}}^{\tau,h} dV \right\}^{-1} \int_{\mathcal{B}^e} \partial_{\alpha^e}^2 \mathbf{d}^e \hat{\pi}_{\text{mix}}^{\tau,h} dV, \qquad (82) $$

$$ \mathbf{K}_{d\bar{\lambda}^{\tau}} = \mathbf{K}_{\bar{\lambda}^{\tau}d}^{\mathrm{T}} := \frac{1}{|\mathcal{B}|} \mathbf{A}_{e=1}^{n_c} \int_{\mathcal{B}^e} \partial_{\mathbf{d}^e \bar{\lambda}^{\tau}}^2 \hat{\pi}_{\text{mix}}^{\tau,h} dV, \qquad (83) $$

$$ \mathbf{K}_{d\bar{\eta}^{\tau}} = \mathbf{K}_{\bar{\eta}^{\tau}d}^{\mathrm{T}} := \frac{1}{|\mathcal{B}|} \mathbf{A}_{e=1}^{n_c} \int_{\mathcal{B}^e} \partial_{\mathbf{d}^e \bar{\eta}^{\tau}}^2 \hat{\pi}_{\text{mix}}^{\tau,h} dV, \qquad (84) $$

$$ \mathbf{K}_{\bar{\lambda}^{\tau}\bar{\lambda}^{\tau}} := -\frac{1}{\kappa_{\bar{\lambda}}^{\tau}} \mathbf{I}, \qquad (85) $$

$$ \mathbf{K}_{\bar{\eta}^{\tau}\bar{\eta}^{\tau}} := -\frac{1}{\kappa_{\bar{\eta}}^{\tau}} \mathbf{I}. \qquad (86) $$

We solve the System of Equations 78 iteratively in a staggered manner. This means that within one time step we repeatedly solve for the microscopic degrees of freedom **d** for fixed Lagrange multipliers $\bar{\lambda}^{\tau}$ and $\bar{\eta}^{\tau}$ until $||\mathbf{R}_{\mathrm{d}}|| \approx 0$ and then solve for $\bar{\lambda}^{\tau}$ and $\bar{\eta}^{\tau}$ for fixed **d**. The incremental solution for **d** follows as

$$ \Delta \mathbf{d} = -[\mathbf{K}_{dd}]^{-1} (\mathbf{K}_{d\bar{\lambda}^{\tau}} \cdot \Delta \bar{\lambda}^{\tau} + \mathbf{K}_{d\bar{\eta}^{\tau}} \cdot \Delta \bar{\eta}^{\tau}) \qquad (87) $$

and can be used as an elimination equation in the System of Equations 78 such that

$$ \begin{bmatrix} \mathbf{K}^{\star}_{\bar{\lambda}^{\tau}\bar{\lambda}^{\tau}} & \mathbf{K}^{\star}_{\bar{\lambda}^{\tau}\bar{\eta}^{\tau}} \\ \mathbf{K}^{\star}_{\bar{\eta}^{\tau}\bar{\lambda}^{\tau}} & \mathbf{K}^{\star}_{\bar{\eta}^{\tau}\bar{\eta}^{\tau}} \end{bmatrix} \cdot \begin{bmatrix} \Delta \bar{\lambda}^{\tau} \\ \Delta \bar{\eta}^{\tau} \end{bmatrix} = - \begin{bmatrix} \mathbf{R}_{\bar{\lambda}^{\tau}} \\ \mathbf{R}_{\bar{\eta}^{\tau}} \end{bmatrix}, \qquad (88) $$

which in turn can be solved for the incremental Lagrange multipliers $\bar{\lambda}^{\tau}$ and $\bar{\eta}^{\tau}$. In the latter system of equations the individual arrays read

$$ \begin{align} \mathbf{K}^{\star}_{\bar{\lambda}^{\tau}\bar{\lambda}^{\tau}} &:= \mathbf{K}_{\bar{\lambda}^{\tau}\bar{\lambda}^{\tau}} - \mathbf{K}_{\bar{\lambda}^{\tau}d} [\mathbf{K}_{dd}]^{-1} \mathbf{K}_{d\bar{\lambda}^{\tau}} \nonumber \\ \mathbf{K}^{\star}_{\bar{\lambda}^{\tau}\bar{\eta}} &:= -\mathbf{K}_{\bar{\lambda}^{\tau}d} [\mathbf{K}_{dd}]^{-1} \mathbf{K}_{d\bar{\eta}^{\tau}}, \tag{89} \end{align} $$

$$ \mathbf{K}^{\star}_{\bar{\eta}\bar{\lambda}^{\tau}} := -\mathbf{K}_{\bar{\eta}d} [\mathbf{K}_{dd}]^{-1} \mathbf{K}_{d\bar{\lambda}^{\tau}}, $$

$$ \mathbf{K}^{\star}_{\bar{\eta}\bar{\eta}^{\tau}} := \mathbf{K}_{\bar{\eta}\bar{\eta}^{\tau}} - \mathbf{K}_{\bar{\eta}^{\tau}d} [\mathbf{K}_{dd}]^{-1} \mathbf{K}_{d\bar{\eta}^{\tau}}. $$
We repeat the above procedure until the residua of both, **d** and {$\overline{\lambda}^{\tau}, \overline{\eta}^{\tau}$}, are within a chosen tolerance.<fnref n="8" />

### 3.1.5 Algorithmically consistent macroscopic response

Once we have solved the microscopic boundary-value problem, we are in a position to determine the effective response of the microstructure within the algorithmic setting. We accomplish this by taking the complete variation of Eq. 47 at a macroscopic material point, which yields in the space-time discrete setting the system of equations

$$ \overline{\delta \pi}_{\text{mix}}^{\tau,h} = \frac{1}{|\mathcal{B}|} \mathbf{A}_{e=1}^{n_c} \begin{bmatrix} \delta \mathbf{d}^e \\ \delta \alpha^e \end{bmatrix} \cdot \begin{bmatrix} \int_{\mathcal{B}^e} \partial \mathbf{d}^e \hat{\pi}_{\text{mix}}^{\tau,h} dV \\ \int_{\mathcal{B}^e} \partial \alpha^e \hat{\pi}_{\text{mix}}^{\tau,h} dV \end{bmatrix} + \begin{bmatrix} \delta \overline{\lambda} \\ \delta \overline{\eta} \\ \delta \overline{S} \end{bmatrix} \cdot \begin{bmatrix} \frac{1}{|\mathcal{B}|} \sum_{e=1}^{n_c} \int_{\mathcal{B}^e} \partial_{\overline{\lambda}} \hat{\pi}_{\text{mix}}^{\tau,h} dV - \frac{1}{\kappa_{\overline{\lambda}}^{\tau}} \overline{\lambda} \\ \frac{1}{|\mathcal{B}|} \sum_{e=1}^{n_c} \int_{\mathcal{B}^e} \partial_{\overline{\eta}} \hat{\pi}_{\text{mix}}^{\tau,h} dV - \frac{1}{\kappa_{\overline{\eta}}^{\tau}} \overline{\eta} \\ \frac{1}{|\mathcal{B}|} \sum_{e=1}^{n_c} \int_{\mathcal{B}^e} \partial_{\overline{S}} \hat{\pi}_{\text{mix}}^{\tau,h} dV \end{bmatrix} = 0. \quad (90) $$

Here, the macroscopic state variables encompassed in $\overline{S}$ are implemented in vector notation according to

$$ \overline{S} = \left[ \delta \overline{h}, \overline{\operatorname{div}} \delta \overline{h}, \operatorname{dev}(\overline{\operatorname{grad}} \delta \overline{h}), \overline{\operatorname{grad}}(\overline{\operatorname{div}} \delta \overline{h}), \overline{\operatorname{grad}}[\operatorname{dev}(\overline{\operatorname{grad}} \delta \overline{h})] \right]^T \quad (91) $$

so that the array of macroscopic driving forces arises as

$$ \partial_{\overline{S}} \overline{\pi}_{\text{mix}}^{\tau,h} = \frac{1}{|\mathcal{B}|} \sum_{e=1}^{n_c} \begin{bmatrix} \tau \int_{\mathcal{B}^e} \partial_h \hat{\phi} \, dV \\ \tau \int_{\mathcal{B}^e} \left(-\mu + \frac{1}{d} \partial_h \hat{\phi} \cdot \mathbf{x}\right) \, dV \\ \tau \int_{\mathcal{B}^e} \partial_h \hat{\phi} \otimes \mathbf{x} \, dV \\ \tau \int_{\mathcal{B}^e} \left[-\mu \mathbf{x} + \frac{1}{2d} \partial_h \hat{\phi} \cdot (\mathbf{x} \otimes \mathbf{x} - \overline{\mathbf{j}})\right] \, dV \\ \tau \int_{\mathcal{B}^e} \frac{1}{2} \partial_h \hat{\phi} \otimes (\mathbf{x} \otimes \mathbf{x} - \overline{\mathbf{j}}) \, dV \end{bmatrix} \quad (92) $$

The latter two sets of equations reflect physically well-motivated energy and dissipation dualities between macroscopic state variables and their individual driving forces. We observe that (i) the microscopic energy density $\hat{\psi}$ only appears in driving forces that are dual to the (gradient of the) divergence of the fluid-mass flux **h** via the relation $\mu := \partial_{\overline{\rho}_f} \hat{\psi}$ and (ii) the microscopic dissipation-potential density $\hat{\phi}$ appears in all the driving forces. In consideration of the incremental macroscopic balance of mass $\overline{\rho}_f = \overline{\rho}_f(t_n) - \tau \overline{\operatorname{div}} \overline{h}$, where $\overline{\rho}_f$ is the macroscopic fluid-mass density, we conclude that the macroscopic

<footnote n="8">⁸We solve the linear system of equations associated with the condensation of local fields such as given in Eq. 77 using LAPACK's routines (E. Anderson et al., 1999). Depending on the values of penalty parameters the associated tangent matrices could be singular. In these cases, least-squares routines based on singular-value decomposition could be used. We solve the global degrees of freedom using the PARDISO solver in the framework of Intel Math Kernel Library (Schenk et al., 2001; E. Wang et al., 2014). Furthermore, to reduce the computation time, we have parallelized the macroscopic finite-element assembly and the associated computations using OpenMP (2011).</footnote>
energy- and dissipation-potential-density functions have the forms<fnref n="9" />

$$ \bar{\psi} := \hat{\psi}(\bar{\rho}_f, \overline{\mathrm{grad}\bar{\rho}_f}), \qquad (93) $$

$$ \bar{\phi} := \hat{\phi}(\bar{h}, \mathrm{div}\bar{h}, \mathrm{dev}(\overline{\mathrm{grad}\bar{h}}), \overline{\mathrm{grad}}(\mathrm{div}\bar{h}), \overline{\mathrm{grad}}[\mathrm{dev}(\overline{\mathrm{grad}\bar{h}})]). \qquad (94) $$

Here, the macroscopic density functions $\bar{\psi}$ and $\bar{\phi}$ are defined as the volume averages of their microscopic counterparts at an equilibrium state such that

$$ \bar{\psi} = \langle \hat{\psi}(\bar{\rho}_f) \rangle_B, \qquad (95) $$

$$ \bar{\phi} = \langle \hat{\phi}(\bar{h}) \rangle_B. \qquad (96) $$

Comparing the list of arguments of the functions given in Eqs. 93 and 94 and of the corresponding functions derived by Sciarra, dell'Isola, and Coussy (2007) reveals consistency with respect to the energetic part; the dissipation-potential-density function of Sciarra, dell'Isola, and Coussy (2007) is however only expanded up to first order.<fnref n="10" /> Similar to the formulation of Sciarra, dell'Isola, and Coussy (2007), the present second-order homogenization formulation naturally extends Darcy flow to more generalized flow types even beyond Brinkman flow. We refer to the review of Ehlers (2022) for a discussion of several kinds of fluid flow through porous media.

**Macroscopic driving forces.** In the following, we will take a closer look at the individual macroscopic driving forces in Eq. 91. We start with the driving force dual to the macroscopic fluid-mass flux $\bar{h}$, which we denote as $\bar{m}$ and find

$$ \bar{m} := \partial_{\bar{h}} \bar{\phi} = \langle \partial_{\bar{h}} \hat{\phi} \rangle_B = \langle -\mathrm{grad} \, \mu - \bar{\lambda} + \mathrm{div}(\bar{\eta} \cdot \mathbf{x}) \rangle_B = \langle -\mu \mathbf{n} + (\bar{\eta} \cdot \mathbf{x}) \cdot \mathbf{n} \rangle_{\partial B} + \bar{\lambda}, \qquad (97) $$

where we have exploited the Euler-Lagrange Equation 51 and the fact that the Lagrange multipliers $\bar{\lambda}$ and $\bar{\eta}$ are constant over the RVE. By furthermore taking into account the Euler-Lagrange Equation 52, we conclude that

$$ \bar{m} = \bar{\lambda} \qquad (98) $$

and thus identify the Lagrange multiplier $\bar{\lambda}$ as the macroscopic driving force dual to the macroscopic fluid-mass flux $\bar{h}$.

Next we identify the macroscopic fluid potential $\bar{\mu}$ as the field that is dual to the divergence of the

<footnote n="9">⁹From $\mathrm{div}\bar{h} = \langle \mathrm{div}\,\mathbf{h} \rangle_B$ it follows that $\bar{\rho}_f = \langle \rho_f \rangle_B$.</footnote>

<footnote n="10">¹⁰We note that the second-grade poro-elastic model of Sciarra, dell'Isola, and Coussy (2007) is based on the theory of mixtures and the consideration of pore-scale effects. In the present work, we consider the bridging between scales with similar characteristic lengths that are both separated from the pore scale.</footnote>
fluid-mass flux $\overline{\overline{\text{div} \, h}}$ and the macroscopic apparent fluid-mass density $\overline{\rho}_f$ such that<fnref n="11" />

$$ \overline{\mu} := \partial_{\overline{\rho}_f} \overline{\psi} - \partial_{\overline{\text{div}} h} \overline{\phi} = \langle \mu - \frac{1}{d} \partial_h \widehat{\phi} \cdot \mathbf{x} \rangle_B = \langle \frac{1}{d} \mu \mathbf{n} \cdot \mathbf{x} \rangle_{\partial B} = \frac{1}{d} \langle \mu \operatorname{tr}(\mathbf{n} \otimes \mathbf{x}) \rangle_{\partial B}. \quad (99) $$

The driving force that is dual to the deviatoric part of the gradient of the fluid-mass flux $\overline{\text{dev}(\text{grad} h)}$ reads

$$ \overline{M} := \partial_{\overline{\text{dev}(\text{grad} h)}} \overline{\phi} = \langle \partial_h \widehat{\phi} \otimes \mathbf{x} \rangle_B = -\langle \text{grad} \, \mu \otimes \mathbf{x} \rangle_B = \langle \mu \mathbf{I} \rangle_B - \langle \mu \mathbf{n} \otimes \mathbf{x} \rangle_{\partial B}. \quad (100) $$

In consideration of the given deviatoric-volumetric split, the latter result can be reformulated by making use of Eq. 99 such that

$$ \overline{M} = \langle \mu - \overline{\mu} \rangle_B \mathbf{I} + \langle \mu \operatorname{dev}(\mathbf{n} \otimes \mathbf{x}) \rangle_{\partial B} = \langle \mu \operatorname{dev}(\mathbf{n} \otimes \mathbf{x}) \rangle_{\partial B}. \quad (101) $$

where in the last step we have taken into account the deviatoric nature of the macroscopic dual field such that $\overline{\text{dev}(\text{grad} h)} : \langle \mu - \overline{\mu} \rangle_B \mathbf{I} \equiv 0$.

We note that $\overline{\mu}$ and $\overline{M}$ could be interpreted as the first-order volumetric and deviatoric contributions of a corresponding *fluid-potential* (or *pressure*) tensor.<fnref n="12" />

The remaining higher-order driving forces dual to the corresponding state variables given in Eq. 91 arise as

$$ \overline{\mu} := \partial_{\overline{\text{grad} \rho_f}} \overline{\psi} - \partial_{\overline{\text{grad} \text{div} h}} \overline{\phi} = \langle \mu \mathbf{x} - \frac{1}{2d} \partial_h \widehat{\phi} \cdot (\mathbf{x} \otimes \mathbf{x} - \overline{\mathbf{j}}) \rangle_B, \quad (102) $$

$$ \overline{M} := \partial_{\overline{\text{grad} \operatorname{dev}(\text{grad} h)}} \overline{\phi} = \langle \frac{1}{2} \partial_h \widehat{\phi} \otimes (\mathbf{x} \otimes \mathbf{x} - \overline{\mathbf{j}}) \rangle_B. \quad (103) $$

We further simplify the latter equations by considering the Euler–Lagrange Equation 51 and the fact that $\overline{M}$ is contracted with the third-order tensor $\overline{\text{grad}}[\overline{\text{dev}(\text{grad} h)}]$, which is deviatoric in its first two indices. Because of the latter, the volumetric contributions associated with the first two indices of $\overline{M}$ affect neither the residuum nor the power density at the macroscopic scale. Hence it follows that

$$ \overline{\mu} = \langle \mu \mathbf{x} - \frac{1}{2d} \operatorname{grad} \mu \cdot (\mathbf{x} \otimes \mathbf{x} - \overline{\mathbf{j}}) \rangle_B, \quad (104) $$

$$ \overline{M} = \langle \frac{1}{2} \operatorname{grad} \mu \otimes (\mathbf{x} \otimes \mathbf{x} - \overline{\mathbf{j}}) - \frac{1}{2d} \mathbf{I} \otimes \operatorname{grad} \mu \cdot (\mathbf{x} \otimes \mathbf{x} - \overline{\mathbf{j}}) \rangle_B. \quad (105) $$

These driving forces are associated with the second-order contributions of the *fluid-potential tensor*, which we denote as *fluid hyper-potentials*.

In order to find an interpretation of the Lagrange multiplier $\overline{\eta}$, we first take a look at the dual field of

<footnote n="11">¹¹The macroscopic fluid potential $\overline{\mu}$ has the same form as the macroscopic (chemical) potential obtained from a first-order homogenization approach (Polukhov and Keip, 2020).</footnote>

<footnote n="12">¹²In the classical (first-order) Darcy theory only the volumetric contribution appears, so that there the fluid potential (or pressure) arises as a scalar-valued quantity (Khurshid, Polukhov, and Keip, 2024).</footnote>
the complete second gradient of fluid-mass flux $\overline{\text{grad}^2\bar{h}}$, which can be given as

$$ \text{sym}^{2,3}(\mathbf{I} \otimes \partial_{\overline{\text{grad}}\rho_f} \bar{\psi}) + \partial_{\overline{\text{grad}}^2\bar{h}} \bar{\phi} = \text{sym}^{2,3}(\mathbf{I} \otimes \bar{\mu} - \bar{\mathbf{M}}) = \bar{\eta} \cdot \frac{1}{2} \langle \mathbf{x} \otimes \mathbf{n} \otimes (\mathbf{x} \otimes \mathbf{x} - \bar{\mathbf{j}}) \rangle_{\partial\mathcal{B}}, \quad (106) $$

where $\text{sym}\diamond_{ijk} = \frac{1}{2}(\diamond_{ijk} + \diamond_{ikj})$; In the latter reformulations, we have used integration by parts, Gauss's law and the Euler-Lagrange Equation 56. By further making use of Gauss's law and the definition of the averaged second-order volume momentum $\bar{\mathbf{j}}$, we finally obtain

$$ \text{sym}(\eta^\dagger) = \text{sym}(\bar{\eta} \cdot \bar{\mathbf{j}}) = \text{sym}(\mathbf{I} \otimes \bar{\mu} - \bar{\mathbf{M}}) \quad (107) $$

where $\eta^\dagger$ has been defined in Eq. 37. Consequently, the Lagrange multiplier $\bar{\eta}$ is associated with the driving force that is dual to the second gradient of the macroscopic fluid-mass flux.

**Macroscopic incremental algorithmic moduli.** To determine the algorithmic macroscopic moduli, we first linearize Eq. 90 in a similar manner as in the previous section. Then we assemble the global system of equations and condense the internal degrees of freedom $\alpha^e = [\alpha_\theta^e, \alpha_\mu^e]^T$ at element level. Thereafter, once an equilibrium state has been obtained, we condense out the global degrees of freedom **d** and afterwards the Lagrange multipliers $\bar{\lambda}^\tau$ and $\bar{\eta}^\tau$.

The elimination equation for the internal degrees of freedom $\alpha^e$ arises as

$$ \Delta\alpha^e = - \left\{ \int_{\mathcal{B}^e} \partial_{\alpha^e}^2 \alpha^e \widehat{\pi}_{\text{mix}}^{\tau,h} dV \right\}^{-1} \left( \int_{\mathcal{B}^e} \partial_{\alpha^e}^2 \mathbf{d}^e \widehat{\pi}_{\text{mix}}^{\tau,h} dV \cdot \Delta\mathbf{d}^e + \int_{\mathcal{B}^e} \partial_{\alpha^e \bar{S}}^2 \widehat{\pi}_{\text{mix}}^{\tau,h} dV \cdot \Delta\bar{S} \right) \quad (108) $$

and includes the sensitivity with respect to the macroscopic state variables $\bar{S}$. As a result, we arrive at the reduced system of equations

$$ \delta \bar{\pi}_{\text{mix}}^{\tau,h,\star} + \Delta \delta \bar{\pi}_{\text{mix}}^{\tau,h,\star} = \begin{bmatrix} \delta \mathbf{d} \\ \delta \bar{\lambda}^{\tau} \\ \delta \bar{\eta}^{\tau} \\ \delta \bar{S} \end{bmatrix} \cdot \left\{ \begin{bmatrix} \cdot \\ \cdot \\ \cdot \\ \mathbf{R}_{\bar{S}} \end{bmatrix} + \begin{bmatrix} \mathbf{K}_{\mathrm{dd}} & \mathbf{K}_{\mathrm{d}\bar{\lambda}^{\tau}} & \mathbf{K}_{\mathrm{d}\bar{\eta}^{\tau}} & \mathbf{K}_{\mathrm{d}\bar{S}} \\ \mathbf{K}_{\bar{\lambda}^{\tau}\mathrm{d}} & \mathbf{K}_{\bar{\lambda}^{\tau}\bar{\lambda}^{\tau}} & \cdot & \cdot \\ \mathbf{K}_{\bar{\eta}^{\tau}\mathrm{d}} & \cdot & \mathbf{K}_{\bar{\eta}^{\tau}\bar{\eta}^{\tau}} & \cdot \\ \mathbf{K}_{\bar{S}\mathrm{d}} & \cdot & \cdot & \mathbf{K}_{\bar{S}\bar{S}} \end{bmatrix} \cdot \begin{bmatrix} \Delta \mathbf{d} \\ \Delta \bar{\lambda}^{\tau} \\ \Delta \bar{\eta}^{\tau} \\ \Delta \bar{S} \end{bmatrix} \right\} = 0, \quad (109) $$

with the individual global arrays provided in Eqs. 79 to 86 and Eq. 89

$$ \mathbf{R}_{\bar{S}} := \frac{1}{|\mathcal{B}|} \sum_{e=1}^{n_e} \int_{\mathcal{B}^e} \partial_{\bar{S}} \widehat{\pi}_{\text{mix}}^{\tau,h,\star} dV, \quad (110) $$

$$ \mathbf{K}_{\mathrm{d}\bar{S}} = \mathbf{K}_{\bar{S}\mathrm{d}}^{\mathrm{T}} := \frac{1}{|\mathcal{B}|} \sum_{e=1}^{n_e} \int_{\mathcal{B}^e} \partial_{\mathrm{d}^e \bar{S}}^2 \widehat{\pi}_{\text{mix}}^{\tau,h,\star} dV \quad (111) $$

$$ \mathbf{K}_{\bar{S}\bar{S}} := \frac{1}{|\mathcal{B}|} \sum_{e=1}^{n_e} \int_{\mathcal{B}^e} \partial_{\bar{S}\bar{S}}^2 \widehat{\pi}_{\text{mix}}^{\tau,h,\star} dV, \quad (112) $$
which are associated with the following individual terms per finite element

$$ \partial_{\mathbf{d}^e \bar{S}}^2 \hat{\pi}_{\text{mix}}^{\tau, h, \star} = \tau \left[ \begin{array}{c} \mathbf{N}_h^T (\partial_{hh}^2 \hat{\phi}) \\ \mathbf{N}_h^T (\frac{1}{d} \partial_{hh}^2 \hat{\phi} \cdot \mathbf{x}) \\ \mathbf{N}_h^T (\partial_{hh}^2 \hat{\phi} \otimes \mathbf{x}) \\ \mathbf{N}_h^T (\frac{1}{2d} \partial_{hh}^2 \hat{\phi} \cdot (\mathbf{x} \otimes \mathbf{x} - \mathbf{j})) \\ \mathbf{N}_h^T (\frac{1}{2} \partial_{hh}^2 \hat{\phi} \otimes (\mathbf{x} \otimes \mathbf{x} - \mathbf{j})) \end{array} \right] \quad (113) $$

$$ \partial_{\bar{S}\bar{S}}^2 \hat{\pi}_{\text{mix}}^{\tau, h, \star} = \tau \left[ \begin{array}{cccc} \partial_{hh}^2 \hat{\phi} & \frac{1}{d} \partial_{hh}^2 \hat{\phi} \cdot \mathbf{x} & \partial_{hh}^2 \hat{\phi} \otimes \mathbf{x} & \frac{1}{2d} \partial_{hh}^2 \hat{\phi} \cdot \mathbf{j} & \frac{1}{2} \tau \partial_{hh}^2 \hat{\phi} \otimes \mathbf{j} \\ \frac{1}{d^2} \mathbf{x} \cdot \partial_{hh}^2 \hat{\phi} \cdot \mathbf{x} & \frac{1}{d} \mathbf{x} \cdot \partial_{hh}^2 \hat{\phi} \otimes \mathbf{x} & \frac{1}{2d^2} \mathbf{x} \cdot \partial_{hh}^2 \hat{\phi} \cdot \mathbf{j} & \frac{1}{2d} \mathbf{x} \cdot \partial_{hh}^2 \hat{\phi} \otimes \mathbf{j} \\ & (\mathbf{x} \otimes \partial_{hh}^2 \hat{\phi})^{\frac{1}{2}, 2} \otimes \mathbf{x} & \frac{1}{2d} (\mathbf{x} \otimes \partial_{hh}^2 \hat{\phi})^{\frac{1}{2}, 2} \cdot \mathbf{j} & \frac{1}{2} (\mathbf{x} \otimes \partial_{hh}^2 \hat{\phi})^{\frac{1}{2}, 2} \otimes \mathbf{j} \\ & & \frac{1}{4d^2} \mathbf{j} \cdot \partial_{hh}^2 \hat{\phi} \cdot \mathbf{j} & \frac{1}{4d} \mathbf{j} \cdot \partial_{hh}^2 \hat{\phi} \otimes \mathbf{j} \\ & \text{sym.} & & \frac{1}{4} (\partial_{hh}^2 \hat{\phi} \otimes \mathbf{j})^{\frac{2}{4}, 2} \otimes \mathbf{j} \end{array} \right] \quad (114) $$

Above, we introduced the second-order tensor $j := \mathbf{x} \otimes \mathbf{x} - \mathbf{j}$ and the notation $\diamond T$ for the transpose of a tensor $\diamond$ with respect to the indices $i$ and $j$ (Ehlers, 2018).

At a microscopic equilibrium point with $||\mathbf{R}_d|| \approx 0$ we could condense $\mathbf{d}$ from the system of equations utilizing

$$ \Delta \mathbf{d} = -[\mathbf{K}_{dd}]^{-1} (\mathbf{K}_{d\bar{\lambda}^\tau} \cdot \Delta \bar{\lambda}^\tau + \mathbf{K}_{d\bar{\eta}^\tau} \cdot \Delta \bar{\eta}^\tau + \mathbf{K}_{d\bar{S}} \cdot \Delta \bar{S}). \quad (115) $$

so that we obtain the reduced symmetric system

$$ \left[ \begin{array}{ccc} \mathbf{K}_{\bar{\lambda}^\tau \bar{\lambda}^\tau} & \mathbf{K}_{\bar{\lambda}^\tau \bar{\eta}^\tau} & \mathbf{K}_{\bar{\eta}^\tau \bar{S}} \\ \mathbf{K}_{\bar{\eta}^\tau \bar{\eta}^\tau} & \mathbf{K}_{\bar{\eta}^\tau \bar{S}} & \mathbf{K}_{\bar{S}\bar{S}} \\ \text{sym.} & \mathbf{K}_{\bar{S}\bar{S}} & \mathbf{K}_{\bar{S}\bar{S}} \end{array} \right] \cdot \left[ \begin{array}{c} \Delta \bar{\lambda}^\tau \\ \Delta \bar{\eta}^\tau \\ \Delta \bar{S} \end{array} \right] = - \left[ \begin{array}{c} \cdot \\ \cdot \\ \mathbf{R}_F \end{array} \right], \quad (116) $$

with the individual arrays given in Eq. 89 and

$$ \mathbf{K}_{\bar{S}\bar{S}}^{*} := \mathbf{K}_{\bar{S}\bar{S}} - \mathbf{K}_{\bar{S}d} [\mathbf{K}_{dd}]^{-1} \mathbf{K}_{d\bar{S}}, \quad (117) $$

$$ \mathbf{K}_{\bar{S}\bar{\lambda}^\tau}^{*} := -\mathbf{K}_{\bar{S}d} [\mathbf{K}_{dd}]^{-1} \mathbf{K}_{d\bar{\lambda}^\tau} \quad (118) $$

$$ \mathbf{K}_{\bar{S}\bar{\eta}^\tau}^{*} := -\mathbf{K}_{\bar{S}d} [\mathbf{K}_{dd}]^{-1} \mathbf{K}_{d\bar{\eta}^\tau}. \quad (119) $$
Finally, by condensing the Lagrange multipliers, we obtain the incremental macroscopic moduli

$$ \bar{\mathbf{C}} := \partial_{\bar{S}\bar{S}}^2 \bar{\pi} = \mathbf{K}_{\bar{S}\bar{S}}^* - \begin{bmatrix} \mathbf{K}_{\bar{S}\bar{\lambda}^\tau}^* \\ \mathbf{K}_{\bar{S}\bar{\eta}^\tau}^* \end{bmatrix}^T \begin{bmatrix} \mathbf{K}_{\bar{\lambda}^\tau\bar{\lambda}^\tau}^* & \mathbf{K}_{\bar{\lambda}^\tau\bar{\eta}^\tau}^* \\ \text{sym.} & \mathbf{K}_{\bar{\eta}^\tau\bar{\eta}^\tau}^* \end{bmatrix}^{-1} \begin{bmatrix} \mathbf{K}_{\bar{\lambda}^\tau\bar{S}}^* \\ \mathbf{K}_{\bar{\eta}^\tau\bar{S}}^* \end{bmatrix}. \quad (120) $$

## 3.2 The macroscopic variational problem of fluid flow through porous media

Having determined all necessary ingredients, we are now in the position to posit the macroscopic variational principle

$$ \bar{h}^* = \inf_{\bar{h} \in W_{\bar{h}}} \left\{ \int_{\bar{B}} (d_t \bar{\psi}(\bar{\rho}_f, \overline{\text{grad}}\bar{\rho}_f) + \bar{\phi}(\bar{S})) \, dV + \int_{\partial \bar{B}_{\bar{t}_f}} \bar{t}_f \cdot \bar{h} \, dA - \int_{\partial \bar{B}_{\bar{\gamma}_f}} \bar{\gamma}_f \cdot \partial_{\bar{n}} \bar{h} \, dA \right\}, \quad (121) $$

with $\partial_{\bar{n}}\bar{h} := \overline{\text{grad}}\bar{h} \cdot \bar{n}$. Since this principle requires $C^1$-continuity of the fluid-mass flux $\bar{h}$ (Khurshid, Polukhov, and Keip, 2024), the corresponding finite-element implementation in two and three spatial dimensions is quite involved. We thus formulate an extended Hu–Washizu variational principle in the subsequent section.

### 3.2.1 Mixed variational formulation of second-order fluid flow through porous media

Complementary to the extended variational formulation at the microscopic scale, we postulate the extended macroscopic Hu–Washizu variational principle

$$ \{\bar{h}^*, \bar{\beta}^*\} = \inf_{\bar{h} \in W_{\bar{h}}} \inf_{\bar{\beta} \in W_{\bar{\beta}}} \inf_{\bar{\omega} \in W_{\bar{\omega}}} \sup_{\bar{\sigma} \in W_{\bar{\sigma}}} \left\{ \int_{\bar{B}} \left( d_t \bar{\psi}(\bar{\rho}_f, \bar{\chi}_f) + \bar{\phi}(\bar{S}_{\text{mix}}) + \bar{\sigma} : (\bar{\beta} - \bar{\omega}) \right. \right. \\ \left. + \frac{\kappa_{\bar{\beta}}}{2} \|\bar{\omega} - \overline{\text{grad}}\bar{h}\|^2 - \frac{1}{\kappa_{\bar{\sigma}}} \|\bar{\sigma}\|^2 \right) dV \right. \\ \left. + \int_{\partial \bar{B}_{\bar{t}_f}} \bar{h} \cdot \bar{t}_f \, dA + \int_{\partial \bar{B}_{\bar{\xi}_f}} \bar{\beta} \cdot \bar{\xi}_f \, dA \right\}, \quad (122) $$

in which $\bar{\beta}$ acts as an additional global field that is equal to the gradient of the macroscopic fluid-mass flux $\overline{\text{grad}}\bar{h}$ in an indirect manner through the auxiliary field $\bar{\omega}$ and the perturbed Lagrange multiplier $\bar{\sigma}$. Furthermore, we have introduced the macroscopic fluid-mass density $\bar{\rho}_f$ and its gradient $\bar{\chi}_f \equiv \overline{\text{grad}}\bar{\rho}_f$. While we determine the former employing the macroscopic mass balance $\dot{\bar{\rho}}_f = -\overline{\text{div}}\bar{h}$, we relate the latter to $\bar{\beta}$ through the relation $\dot{\bar{\chi}}_f := -\overline{\text{div}}(\bar{\beta}^T)$.
The corresponding function spaces of the primary fields are

$$
\begin{align}
\mathcal{W}_{\bar{h}} &:= \{\bar{h} \in H(\mathrm{div}, \bar{\mathcal{B}}) \cap [H(\mathrm{grad}, \bar{\mathcal{B}})]^d \mid \bar{h} = \bar{h}_D \text{ on } \partial \bar{\mathcal{B}}_{\bar{h}}\}, \nonumber \\
\mathcal{W}_{\bar{\beta}} &:= \{\bar{\beta} \in [H(\mathrm{div}, \bar{\mathcal{B}})]^d \cap [H(\mathrm{grad}, \bar{\mathcal{B}})]^{d \times d} \mid \bar{\beta} = \bar{\beta}_D \text{ on } \partial \bar{\mathcal{B}}_{\bar{\beta}}\}, \tag{123} \\
\mathcal{W}_{\bar{\omega}} &:= \{\bar{\omega} \in [L^2(\bar{\mathcal{B}})]^{d \times d}\}, \nonumber \\
\mathcal{W}_{\bar{\sigma}} &:= \{\bar{\sigma} \in [L^2(\bar{\mathcal{B}})]^{d \times d}\}. \nonumber
\end{align}
$$

with the decomposition of the boundary of the body $\bar{\mathcal{B}}$ into Dirichlet and Neumann boundaries such that

$$
\partial \bar{\mathcal{B}} = \partial \bar{\mathcal{B}}_{\bar{h}} \cup \partial \bar{\mathcal{B}}_{\bar{t}_f} = \partial \bar{\mathcal{B}}_{\bar{\beta}} \cup \partial \bar{\mathcal{B}}_{\bar{\xi}_f} \quad \text{and} \quad \partial \bar{\mathcal{B}}_{\bar{h}} \cap \partial \bar{\mathcal{B}}_{\bar{t}_f} = \partial \bar{\mathcal{B}}_{\bar{\beta}} \cap \partial \bar{\mathcal{B}}_{\bar{\xi}_f} = \emptyset. \tag{124}
$$

Associated with the extended mixed variational principle, we identify the modified macroscopic state
array

$$
\bar{S}_{\text{mix}} := \{ \bar{h}, \overline{\operatorname{div} \bar{h}}, \overline{\operatorname{div}(\bar{\beta}^T)}, \operatorname{dev} \bar{\beta}, \overline{\operatorname{grad}(\operatorname{dev} \bar{\beta})} \}. \quad (125)
$$

We obtain the weak form of the variational principle by using the fundamental lemma of variational
calculus as

$$
\delta \bar{\Pi} = \int_{\bar{\mathcal{B}}} \begin{bmatrix} \delta \bar{h} \\ \overline{\mathrm{div}} \delta \bar{h} \\ \overline{\mathrm{div}}(\delta \bar{\beta}^T) \\ \mathrm{dev} \, \delta \bar{\beta} \\ \overline{\mathrm{grad}}(\mathrm{dev} \, \delta \bar{\beta}) \\ \overline{\mathrm{grad}} \delta \bar{h} \\ \delta \bar{\beta} \\ \delta \bar{\omega} \\ \delta \bar{\sigma} \end{bmatrix} \cdot \begin{bmatrix} \partial_{\bar{h}} \bar{\phi} \\ - \partial_{\bar{\rho}_f} \bar{\psi} \\ - \partial_{\bar{\chi}_f} \bar{\psi} \\ - \partial_{\mathrm{dev} \, \bar{\beta}} \bar{\phi} \\ - \partial_{\overline{\mathrm{grad}}(\mathrm{dev} \, \bar{\beta})} \bar{\phi} \\ - \kappa_{\bar{\beta}} (\bar{\omega} - \overline{\mathrm{grad} \, \bar{h}}) \\ \bar{\sigma} \\ - \bar{\sigma} + \kappa_{\bar{\beta}} (\bar{\omega} - \overline{\mathrm{grad} \, \bar{h}}) \\ \bar{\beta} - \bar{\omega} - \frac{1}{\kappa_{\bar{\sigma}}} \bar{\sigma} \end{bmatrix} \, dV + \int_{\partial \bar{\mathcal{B}}_{\bar{t}_f}} \delta \bar{h} \cdot \bar{t}_f \, dA + \int_{\partial \bar{\mathcal{B}}_{\bar{\xi}_f}} \delta \bar{\beta} \cdot \bar{\xi}_f \, dA = 0, \quad (126)
$$

The corresponding Euler–Lagrange equations of the macroscopic initial-boundary-value problem then
read

$$
\begin{align}
\delta \bar{h} : & \partial_{\bar{h}} \bar{\phi} + \overline{\mathrm{grad}}(\partial_{\bar{\rho}_f} \bar{\psi}) + \overline{\mathrm{div}}[\kappa_{\bar{\beta}}(\bar{\omega} - \overline{\mathrm{grad}} \bar{h})] = \mathbf{0} \quad \text{in } \bar{B} \nonumber \\
\delta \bar{h} : & \partial_{\bar{\rho}_f} \bar{\psi} \bar{n} + \kappa_{\bar{\beta}}(\bar{\omega} - \overline{\mathrm{grad}} \bar{h}) \cdot \bar{n} = \bar{t}_f \quad \text{on } \partial \bar{B}_{\bar{t}_f} \nonumber \\
\delta \bar{\beta} : & \bar{\sigma} - \partial_{\mathrm{dev}} \bar{\beta} \bar{\phi} + \overline{\mathrm{grad}}(\partial_{\bar{\chi}_f} \bar{\psi}) + \overline{\mathrm{div}}(\partial_{\overline{\mathrm{grad}}(\mathrm{dev} \bar{\beta})} \bar{\phi}) = \mathbf{0} \quad \text{in } \bar{B} \nonumber \\
\delta \bar{\beta} : & \bar{n} \otimes \partial_{\bar{\chi}_f} \bar{\psi} + \partial_{\overline{\mathrm{grad}}(\mathrm{dev} \bar{\beta})} \bar{\phi} \cdot \bar{n} = \bar{\xi}_f \quad \text{on } \partial \bar{B}_{\bar{\xi}_f} \nonumber \\
\delta \bar{\omega} : & \kappa_{\bar{\beta}}(\bar{\omega} - \overline{\mathrm{grad}} \bar{h}) - \bar{\sigma} = \mathbf{0} \quad \text{in } \bar{B} \nonumber \\
\delta \bar{\sigma} : & \bar{\beta} - \bar{\omega} - \frac{1}{\kappa_{\bar{\sigma}}} \bar{\sigma} = \mathbf{0} \quad \text{in } \bar{B} \tag{127}
\end{align}
$$

where the Lagrange multiplier $\bar{\sigma}$ is a hyper-field associated with the macroscopic first- and second-order fluid-potential tensors. Utilizing Eqs. 100, 102 and 103, we can identify $\bar{\sigma}$ in terms of the macroscopic dual forces as $\bar{\sigma} = \bar{M} - \overline{\mathrm{grad}}\bar{\mu}^{\psi} - \overline{\mathrm{div}}\bar{M}$ with $\mu^{\psi}$ denoting the energetic part of $\mu$; see also Khurshid, Polukhov, and Keip, 2024.

### 3.2.2 Algorithmic implementation of macroscopic variational principle

The algorithmic implementation of the macroscopic variational principle mainly follows the one of the microscopic variational principle. In fact, once the microscopic problem has been discretized and the homogenization procedure has been implemented (Section 3), many ingredients for the algorithmic implementation of the macroscopic problem are readily available. We thus only briefly discuss the time and space discretizations using non-conforming mixed finite elements in the following.

**Time discretization of the macroscopic problem.** To arrive at an algorithmically consistent and unconditionally stable multiscale model, we employ an implicit Euler scheme for time integration. Analogous to the microscopic case, we determine the incremental potential within a time interval $[t_n, t_{n+1}]$ as

$$
\begin{align}
\bar{\Pi}^{\tau}(\bar{h}, \bar{\beta}, \bar{\omega}, \bar{\sigma}) = & \int_{\bar{B}} \bar{\psi}(\bar{\rho}_f, \bar{\chi}_f) - \bar{\psi}(\bar{\rho}_f(t_n), \bar{\chi}_f(t_n)) + \tau \bar{\phi}(\bar{S}_{\mathrm{mix}}) + \tau \bar{\sigma} : (\bar{\beta} - \bar{\omega}) \nonumber \\
& + \frac{\tau \kappa_{\bar{\beta}}}{2} \|\bar{\omega} - \overline{\mathrm{grad}} \bar{h}\|^2 - \frac{\tau}{\kappa_{\bar{\sigma}}} \|\bar{\sigma}\|^2 \, \mathrm{d}V + \int_{\partial \bar{B}_{\bar{t}_f}} \tau \bar{h} \cdot \bar{t}_f \, \mathrm{d}A + \int_{\partial \bar{B}_{\bar{\xi}_f}} \tau \bar{\beta} \cdot \bar{\xi}_f \, \mathrm{d}A \tag{128}
\end{align}
$$

where we again denote with $\tau = t_{n+1} - t_n$ the time-step size and drop the indices for variables that are evaluated at the current time $t_{n+1}$. Furthermore, and in analogy to the microscopic problem, we use the macroscopic mass balance to determine the macroscopic fluid-mass density $\bar{\rho}_f = \bar{\rho}_f(t_n) - \tau \overline{\mathrm{div}} \bar{h}$ and its gradient $\bar{\chi}_f = \bar{\chi}_f(t_n) - \tau \overline{\mathrm{div}}(\bar{\beta}^T)$.

**Finite-element discretization of the macroscopic problem.** Like in the case of the microscopic
problem we select second-order interpolation functions for the fluid-mass flux $\bar{h}$ such that

$$ 
\begin{align} 
\bar{h} &\approx \sum_{I=1}^{n_h^h} \bar{N}_h^I \bar{d}_h^I &=: \bar{N}_h \bar{d}_h^e, \nonumber \\ 
\overline{\mathrm{div}} \bar{h} &\approx \sum_{I=1}^{n_h^h} \bar{d}_h^I \cdot \overline{\mathrm{grad}} \bar{N}_h^I &=: \bar{B}_h^{\mathrm{div}} \bar{d}_h^e, \tag{129} \\ 
\overline{\mathrm{grad}} \bar{h} &\approx \sum_{I=1}^{n_h^h} \bar{d}_h^I \otimes \overline{\mathrm{grad}} \bar{N}_h^I &=: \bar{B}_h^{\mathrm{grad}} \bar{d}_h^e. \nonumber 
\end{align} 
$$

We discretize the macroscopic domain as well as the macroscopic fluid-mass flux by using eight-noded quadrilateral finite elements ($Q_8$) in two spatial dimensions and twenty-noded hexahedral finite elements (HEX$_{20}$) in three dimensions.<fnref n="13" /> For the numerical integration we use four Gauss points in two dimensions and eight Gauss points in three dimensions.

To arrive at matching orders of interpolations for the constraints of the variational principle, $\bar{\beta}$ is approximated using linear interpolation functions. The corresponding interpolations then arise in compact form as

$$ 
\begin{align} 
\bar{\beta} &\approx \sum_{I=1}^{n_h^\beta} \bar{N}_\beta^I \bar{d}_\beta^I &=: \bar{N}_\beta \bar{d}_\beta^e, \tag{130} \\ 
\mathrm{dev} \bar{\beta} &= \bar{\beta} - \frac{1}{d} \mathrm{tr} \bar{\beta} \mathbf{I} \approx: \bar{N}_\beta^{\mathrm{dev}} \bar{d}_\beta^e, \tag{131} \\ 
\overline{\mathrm{div}}(\bar{\beta}^T) &\approx \sum_{I=1}^{n_h^\beta} \bar{d}_\beta^I \cdot \overline{\mathrm{grad}} \bar{N}_\beta^I &=: \bar{B}_\beta^{\mathrm{div}} \bar{d}_\beta^e, \tag{132} \\ 
\overline{\mathrm{grad}}(\mathrm{dev} \bar{\beta}) &= \overline{\mathrm{grad}} \bar{\beta} - \frac{1}{d} \mathbf{I} \otimes \overline{\mathrm{div}}(\bar{\beta}^T) \approx: \bar{B}_\beta^{\mathrm{grad}} \bar{d}_\beta^e. \tag{133} 
\end{align} 
$$

Analogous to the approach documented in Section 3.1.4, the internal fields $\bar{\omega}$ and $\bar{\sigma}$ that are associated with the constraints are discretized considering linear, discontinuous interpolation functions such that

$$ 
\begin{align} 
\bar{\omega} &\approx \sum_{I=1}^{n_h^\omega} \bar{N}_\omega^I \bar{\alpha}_\omega^I &=: \bar{N}_\omega \bar{\alpha}_\omega \tag{134} \\ 
\bar{\sigma} &\approx \sum_{I=1}^{n_h^\sigma} \bar{N}_\sigma^I \bar{\alpha}_\sigma^I &=: \bar{N}_\sigma \bar{\alpha}_\sigma. \tag{135} 
\end{align} 
$$

<footnote n="13">¹³While we have generally observed stable response in case of two-dimensional simulations, three-dimensional problems can become unstable depending on the choice of the microscopic and macroscopic penalty parameters. We believe that these issues stem from the high number of degrees of freedom associated with Lagrange multipliers in three dimensions, which lead to ill-conditioning of local and global matrices.</footnote>
which allows us to condense the corresponding internal degrees of freedom from the global system of equations. Our implementation is mainly along the lines of the one that is presented in Section 3.1.4.

In consideration of the above derivations, the necessary condition of the space-time-discrete variational principle results for a pure Dirichlet problem in the statement

$$ \delta \bar{\Pi}_{\text{mix}}^{\tau,h} = \sum_{e=1}^{\bar{n}_e} \int_{\bar{B}^e} \begin{bmatrix} \delta \bar{\mathbf{d}}_h^e \\ \delta \bar{\mathbf{d}}_\beta^e \\ \delta \bar{\alpha}_\omega^e \\ \delta \bar{\alpha}_\sigma^e \end{bmatrix} \cdot \begin{bmatrix} \bar{\mathbf{r}}_h^e \\ \bar{\mathbf{r}}_\beta^e \\ \bar{\mathbf{r}}_\omega^e \\ \bar{\mathbf{r}}_\sigma^e \end{bmatrix} \, dV, \qquad (136) $$

where we have introduced the element-wise residual vector $\bar{\mathbf{r}}^e := [\bar{\mathbf{r}}_h^e, \bar{\mathbf{r}}_\beta^e, \bar{\mathbf{r}}_\omega^e, \bar{\mathbf{r}}_\sigma^e]^T$ with the individual entries

$$ \bar{\mathbf{r}}_h^e := (\bar{\mathbf{N}}_h)^T (\tau \partial_{\bar{h}} \bar{\phi}) - (\bar{\mathbf{B}}_h^{\text{div}})^T (\tau \partial_{\bar{\rho}_f} \bar{\psi}) - (\bar{\mathbf{B}}_h^{\text{grad}})^T [\tau \kappa_{\bar{\beta}} (\bar{\omega} - \overline{\text{grad} \bar{h}})] \quad (137) $$

$$ \bar{\mathbf{r}}_{\beta}^{e} := (\bar{\mathbf{N}}_{\beta})^{T} (\tau \bar{\sigma}) - (\bar{\mathbf{B}}_{\beta}^{\text{div}})^{T} (\tau \partial_{\bar{\chi}_{f}} \bar{\psi}) - (\bar{\mathbf{N}}_{\beta}^{\text{dev}})^{T} (\tau \partial_{\text{dev}} \bar{\beta} \bar{\phi}) - (\bar{\mathbf{B}}_{\beta}^{\text{grad}})^{T} (\tau \partial_{\overline{\text{grad}}(\text{dev} \bar{\beta})} \bar{\phi}) \quad (138) $$

$$ \bar{\mathbf{r}}_{\omega}^{e} := (\bar{\mathbf{N}}_{\omega})^{T} \tau [\bar{\sigma} + \kappa_{\bar{\beta}} (\bar{\omega} - \overline{\text{grad} \bar{h}})] \quad (139) $$

$$ \bar{\mathbf{r}}_{\sigma}^{e} := (\bar{\mathbf{N}}_{\sigma})^{T} \tau (\bar{\beta} - \bar{\omega} - \frac{1}{\kappa_{\bar{\sigma}}} \bar{\sigma}) \quad (140) $$

and took into account that $-\tau \partial_{\bar{\rho}_f} \bar{\psi} = \partial_{\overline{\text{div}} \bar{h}} \bar{\psi}$ and $-\tau \partial_{\bar{\chi}_f} \bar{\psi} = \partial_{\overline{\text{div}} (\bar{\beta})^T} \bar{\psi}$.

Consistent linearization finally leads to the statement

$$ \delta \bar{\Pi}_{\text{mix}}^{\tau,h} + \Delta \delta \bar{\Pi}_{\text{mix}}^{\tau,h} = \bar{\mathbf{A}}_{e=1}^{\bar{n}_e} \begin{bmatrix} \delta \bar{\mathbf{d}}^e \\ \delta \bar{\alpha}^e \end{bmatrix} \cdot \left\{ \int_{\bar{B}^e} \begin{bmatrix} \partial_{\bar{\mathbf{d}}^e} \bar{\pi} \\ \partial_{\bar{\alpha}^e} \bar{\pi} \end{bmatrix} dV + \int_{\bar{B}^e} \begin{bmatrix} \partial_{\bar{\mathbf{d}}^e}^2 \bar{\pi} \\ \partial_{\bar{\alpha}^e}^2 \bar{\pi} \\ \partial_{\bar{\mathbf{d}}^e} \bar{\alpha}^e \\ \partial_{\bar{\alpha}^e} \bar{\alpha}^e \end{bmatrix} dV \cdot \begin{bmatrix} \Delta \bar{\mathbf{d}}^e \\ \Delta \bar{\alpha}^e \end{bmatrix} \right\} = 0, \quad (141) $$

where $\bar{\alpha}^e := [\bar{\alpha}_\omega^e, \bar{\alpha}_\sigma^e]^T$ and where we take into account the effective properties obtained from the homogenization over RVEs as discussed in the previous sections.

Having described both the microscopic and the macroscopic variational formulation as well as the associated spatiotemporal discretizations, we are in a position to numerically simulate two-scale initial-boundary-value problems. These will be discussed in the next section.

# 4 Multiscale simulations of Darcy and Darcy-Forchheimer flow through porous media

In what follows, we discuss a sequence of computational multiscale initial-boundary-value problems associated with the fluid flow through porous media in order to analyze the response of heterogeneous
porous media across length scales. For that, we employ both, first- and second-order FE² methods under consideration of representative volume elements, as well as direct numerical simulations (DNS). The latter will be based on the formulation presented in Section 2 and discretized by means of H(div)-conforming, second-order Brezzi–Douglas–Fortin–Marini (BDFM₂) finite elements on quadrilateral meshes. For the multiscale simulations, we discretize both the microscopic and the macroscopic fields by means of Lagrangian interpolation functions of mixed type as discussed in Section 3. Both the DNS and the multiscale FE² frameworks have been implemented in FEAP (Taylor and Govindjee, 2021, Version 8.6.1).

For the modeling of the material behavior we employ energy-density and dissipation-potential-density functions. On the energetic side, we use the energy-density function (Miehe, Mauthe, and Teichtmeister, 2015; Kaessmair and Steinmann, 2018; Polukhov and Keip, 2020)

$$ \widehat{\psi}(\tilde{\rho}_f) = \frac{A}{2}(\tilde{\rho}_f - \tilde{\rho}_f(t_0))^2, \qquad (142) $$

where A is a material parameter that characterizes the compressibility of the fluid. It is related to a Biot-type modulus M and the fluid's intrinsic mass density $\rho_f$ via $A := M/\rho_f^2$. Unless otherwise mentioned, we will consider a vanishing initial fluid-mass density $\tilde{\rho}_f(t_0) = 0$ in the following studies. On the dissipative side, we take into account the convex dissipation-potential-density function<fnref n="14" />

$$ \widehat{\phi}(\mathbf{h}) = \frac{1}{2\rho_f^2 k} \|\mathbf{h}\|^2 + \frac{\gamma}{3} \|\mathbf{h}\|^3 = \frac{1}{2K} \|\mathbf{h}\|^2 + \frac{\gamma}{3} \|\mathbf{h}\|^3, \qquad (143) $$

where $k > 0$ denotes the isotropic permeability of the material and $\gamma \ge 0$ is the coefficient of the Forchheimer law (Ehlers, 2022; Monchiet, 2024); The permeability k can be expressed in terms of the intrinsic permeability $\mathbf{R}$ and the fluid viscosity $\eta_f$ via $k = \mathbf{R}/\eta_f$, where the intrinsic permeability correlates with the morphology of the pore network (Coussy, 2004). The above function is able to describe both Darcy and Darcy–Forchheimer fluid flow. Contrary to the linear Darcy law, the Forchheimer law establishes a nonlinear relation between the fluid-mass flux and the gradient of the fluid potential (pore pressure). While the Darcy law is thus only suitable for low velocities of the fluid-mass flux, the Darcy–Forchheimer law is predictive even in case of higher flow rates (Zeng and Grigg, 2006).

## 4.1 Fluid flow through porous media in two spatial dimensions

We now present results that have been obtained from multiscale numerical simulations using the second-order FE² approach documented above. The results are associated with different combinations of the spatial dimensions of the microscopic and macroscopic simulation domains. To evaluate the performance of the proposed method, we compare its results to those obtained from a corresponding

<footnote n="14">¹⁴Although not explicitly addressed in this work, we note that the proposed variational principle could seamlessly integrate Brinkman-type fluid flow in porous media by incorporating the gradient of the fluid-mass flux into the dissipation-potential-density function.</footnote>
<table><caption>**Table 1:** Values of material parameters of the microstructures shown in Figs. 5 to 7.</caption><thead><tr><th>Phases</th><th>A<br>m<sup>5</sup>/(kg s<sup>2</sup>)<br>(47), (142)</th><th>K<br>m<sup>3</sup> s/kg<br>(47), (142)</th><th>κ<sub>μ</sub><br>m<sup>5</sup>/(kg s)<br>(47)</th><th>κ<sub>λ</sub><sup>τ</sup><br>m<sup>3</sup>/kg<br>(46)</th><th>κ<sub>η</sub><sup>τ</sup><br>m<sup>3</sup>/kg<br>(46)</th><th>κ<sub>β</sub><br>m<sup>5</sup>/(kg s)<br>(128)</th><th>κ<sub>σ</sub><br>m<sup>5</sup>/(kg s)<br>(128)</th><th>γ<br>m<sup>5</sup>/kg<sup>2</sup><br>(128)</th></tr></thead><tbody><tr><td>Matrix</td><td>10</td><td>10<sup>-1</sup></td><td>10<sup>-15</sup></td><td>10<sup>15</sup></td><td>10<sup>15</sup></td><td>10<sup>6</sup></td><td>10<sup>14</sup></td><td>0</td></tr><tr><td>Inclusion</td><td>5</td><td>10<sup>-2</sup></td><td>10<sup>-15</sup></td><td>10<sup>15</sup></td><td>10<sup>15</sup></td><td>10<sup>6</sup></td><td>10<sup>14</sup></td><td>0</td></tr></tbody></table>

<table><caption>**Table 2:** Values of material parameters of the microstructure shown in Fig. 8.</caption><thead><tr><th>Phases</th><th>A<br>m<sup>5</sup>/(kg s<sup>2</sup>)<br>(47), (142)</th><th>K<br>m<sup>3</sup> s/kg<br>(47), (142)</th><th>κ<sub>μ</sub><br>m<sup>5</sup>/(kg s)<br>(47)</th><th>κ<sub>λ</sub><sup>τ</sup><br>m<sup>3</sup>/kg<br>(46)</th><th>κ<sub>η</sub><sup>τ</sup><br>m<sup>3</sup>/kg<br>(46)</th><th>κ<sub>β</sub><br>m<sup>5</sup>/(kg s)<br>(128)</th><th>κ<sub>σ</sub><br>m<sup>5</sup>/(kg s)<br>(128)</th><th>γ<br>m<sup>5</sup>/kg<sup>2</sup><br>(128)</th></tr></thead><tbody><tr><td>Matrix</td><td>10</td><td>10<sup>-1</sup></td><td>10<sup>-15</sup></td><td>10<sup>15</sup></td><td>10<sup>15</sup></td><td>10<sup>6</sup></td><td>10<sup>14</sup></td><td>0</td></tr><tr><td>Inclusion (green)</td><td>5</td><td>10<sup>-2</sup></td><td>10<sup>-15</sup></td><td>10<sup>15</sup></td><td>10<sup>15</sup></td><td>10<sup>6</sup></td><td>10<sup>14</sup></td><td>0</td></tr><tr><td>Inclusion (blue)</td><td>5</td><td>10<sup>-4</sup></td><td>10<sup>-15</sup></td><td>10<sup>15</sup></td><td>10<sup>15</sup></td><td>10<sup>6</sup></td><td>10<sup>14</sup></td><td>0</td></tr></tbody></table>

<table><caption>**Table 3:** Values of material parameters for the microstructures shown in Fig. 10.</caption><thead><tr><th>Phases</th><th>A<br>m<sup>5</sup>/(kg s<sup>2</sup>)<br>(47), (142)</th><th>K<br>m<sup>3</sup> s/kg<br>(47), (142)</th><th>κ<sub>μ</sub><br>m<sup>5</sup>/(kg s)<br>(47)</th><th>κ<sub>λ</sub><sup>τ</sup><br>m<sup>3</sup>/kg<br>(46)</th><th>κ<sub>η</sub><sup>τ</sup><br>m<sup>3</sup>/kg<br>(46)</th><th>κ<sub>β</sub><br>m<sup>5</sup>/(kg s)<br>(128)</th><th>κ<sub>σ</sub><br>m<sup>5</sup>/(kg s)<br>(128)</th><th>γ<br>m<sup>5</sup>/kg<sup>2</sup><br>(128)</th></tr></thead><tbody><tr><td>Matrix</td><td>10</td><td>10<sup>-1</sup></td><td>10<sup>-15</sup></td><td>10<sup>15</sup></td><td>10<sup>15</sup></td><td>10<sup>6</sup></td><td>10<sup>14</sup></td><td>0</td></tr><tr><td>Inclusions</td><td>5</td><td>10<sup>-3</sup></td><td>10<sup>-15</sup></td><td>10<sup>15</sup></td><td>10<sup>15</sup></td><td>10<sup>6</sup></td><td>10<sup>14</sup></td><td>0</td></tr></tbody></table>

<table><caption>**Table 4:** Values of material parameters for the microstructures shown in Fig. 11.</caption><thead><tr><th>Phases</th><th>A<br>m<sup>5</sup>/(kg s<sup>2</sup>)<br>(47), (142)</th><th>K<br>m<sup>3</sup> s/kg<br>(47), (142)</th><th>κ<sub>μ</sub><br>m<sup>5</sup>/(kg s)<br>(47)</th><th>κ<sub>λ</sub><sup>τ</sup><br>m<sup>3</sup>/kg<br>(46)</th><th>κ<sub>η</sub><sup>τ</sup><br>m<sup>3</sup>/kg<br>(46)</th><th>κ<sub>β</sub><br>m<sup>5</sup>/(kg s)<br>(128)</th><th>κ<sub>σ</sub><br>m<sup>5</sup>/(kg s)<br>(128)</th><th>γ<br>m<sup>5</sup>/kg<sup>2</sup><br>(128)</th></tr></thead><tbody><tr><td>Matrix</td><td>10</td><td>10<sup>-1</sup></td><td>10<sup>-15</sup></td><td>10<sup>15</sup></td><td>10<sup>15</sup></td><td>10<sup>6</sup></td><td>10<sup>14</sup></td><td>100</td></tr><tr><td>Inclusions</td><td>5</td><td>10<sup>-3</sup></td><td>10<sup>-15</sup></td><td>10<sup>15</sup></td><td>10<sup>15</sup></td><td>10<sup>6</sup></td><td>10<sup>14</sup></td><td>0</td></tr></tbody></table>
first-order homogenization approach (Polukhov and Keip, 2020) and direct numerical simulations (DNS).

We make several remarks regarding the following initial-boundary-value problems. Most of the two-dimensional simulations will be performed on rectangular macroscopic domains that are discretized with quadrilateral finite elements having four Gauss points. The dimensions of each macroscopic finite element are chosen to correspond to a patch of four neighboring RVEs. This approach enables the utilization of the same microstructures and finite-element meshes for both multiscale simulations and DNS simulations (Fig. 4).<fnref n="15" /> Given this discretization, it will be possible to compare microscopic and macroscopic results that have been obtained from different methods. In this context, we will study the microscopic and macroscopic fluid-flux densities $h$ and $\bar{h}$ obtained from first- and second-order FE² simulations as well as from DNS. When associated results are presented in corresponding diagrams, we refer to them in the following way:

* "DNS": Data obtained from the direct numerical simulation of a microstructure. The associated fluid-mass flux is denoted as $h$ and determined from Eq. 7.

* "FE²-RVE": Data obtained from the FE² simulation of a multiscale problem and extracted from the finite-element discretization of the microscopic scale. The associated fluid-mass flux is denoted as $h$ and determined from Eq. 47 using the second-order expansion given in Eq. 23 or determined from the first-order homogenization approach proposed by Polukhov and Keip (2020) in combination with the constraint on average fluctuations of the fluid-mass flux given in Eq. 36.

* "FE²-macro": Data obtained from the FE² simulation of a multiscale problem and extracted from the finite-element discretization of the macroscopic scale. The associated fluid-mass flux is denoted as $\bar{h}$ and determined from Eq. 122.

### 4.1.1 Multiscale behavior of differently sized macroscopic bodies with fixed-sized RVEs

In the present study, we analyze the fluid flow through a two-phase porous medium of different macroscopic size via multiscale simulations and DNS. The porous medium is assumed to have a periodic microstructure with fixed dimensions ($l_1 \times l_2$) = (1 × 1) m² that is composed of a permeable matrix and less permeable circular inclusions with radius $r$ = 0.25 m (Fig. 4 and Table 1). As macroscopic boundary conditions we prescribe

$$ \partial \mathcal{B}_{\text{left}} : \bar{h}_1(\mathbf{x}, t) = 0, \qquad (144) $$

$$ \partial \mathcal{B}_{\text{right}} : \bar{h}_1(\mathbf{x}, t) = 0, \qquad (145) $$

<footnote n="15">¹⁵This consideration facilitates for intuitive comparisons of results of the FE² approach with those of the direct numerical simulations. Nevertheless, in certain scenarios finer macroscopic meshes might become necessary, for example, when steep local gradients of fields are to be expected. In these cases, macroscopic finite elements can become smaller than the attached RVEs and, as a result, RVEs of neighboring Gauss points can overlap. We analyze such scenarios and the corresponding effects in Appendix A and observe that regardless of the resolution of a macroscopic finite-element mesh, the macroscopic response remains consistent for given RVE size.</footnote>
<figure><img src="image_6.png" /></figure>

**Figure 4:** *A representative macroscopic body with microstructure.* In the computational studies documented in the present work we consider macroscopic bodies with microstructure and determine their response by means of (i) direct numerical simulations of the fully resolved problem at the macroscopic scale, (ii) first-order FE² simulations, and (iii) second-order FE² simulations. To allow for meaningful comparisons of the results predicted by the different approaches, we take into account periodic microstructures that we either (i) discretize directly at the macroscopic scale or (ii,iii) map into representative volume elements (RVEs) of exact same size. In this way, we can compare the finely resolved results based on identical microstructures for all considered models. Nevertheless, it should be kept in mind that certain mismatches between the individual results are to be expected because the location of a macroscopic Gauss point needs to coincide with the center of an attached RVE (Eq. 15). However, as will become clear from the following studies, the associated discrepancies are rather low for the initial-boundary-value problems considered in the present work.
<figure><img src="image_7.png" /><figcaption>**Figure 5:** Profiles of the vertical component of the fluid-mass flux for particulate microstructures with 4, 8, 16 and 24 inclusions. The dimensions of the periodic RVEs of each microstructure are given as ($l_1 \times l_2$) = (1 × 1) m² such that the height of the macroscopic domain increases with the number of REVs. The profiles of the fluid-mass flux have been extracted from the material points located on the dashed lines drawn on the inset microstructure images.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 6:** Profiles of the horizontal component of the fluid-mass flux for particulate microstructures with 4, 8, 16 and 24 inclusions. The dimensions of the periodic RVEs of each microstructure are given as ($l_1 \times l_2$) = (1 $\times$ 1) m$^2$ such that the height of the macroscopic domain increases with the number of REVs. The profiles of the fluid-mass flux have been extracted from the material points located on the dashed lines drawn on the inset microstructure images.</figcaption></figure>
$$ \partial \mathcal{B}_{\text{top}} : \bar{h}_2(\mathbf{x}, t) = -1 \times 10^{-3} \, \text{kg/(m}^2 \text{ s)}, \qquad (146) $$

$$ \partial \mathcal{B}_{\text{bottom}} : \bar{h}_2(\mathbf{x}, t) = 0, \qquad (147) $$

$$ \partial \mathcal{B} : \bar{\beta}_{12}(\mathbf{x}, t) = \bar{\beta}_{21}(\mathbf{x}, t) = 0. \qquad (148) $$

where the last one is only relevant for the second-order simulations.

In Fig. 5, we plot the simulation results of the vertical component of the fluid-mass flux along specific sections as shown in the inset microstructures. We observe that the first- and second-order homogenization as well as the DNS results are in good agreement. However, when the size of the RVEs is comparable to the macroscopic length scale, the macroscopic boundary conditions at the upper and lower boundaries are not captured well within the RVEs close to the associated boundaries. This trend could be explained by the increasing distance of the macroscopic Gauss points from the boundary. Furthermore, we observe that the fluid-mass flux jumps across the boundaries of neighboring RVEs, with jumps being more pronounced in case of the first-order formulation. Given the fact that the corresponding RVEs are driven by first- and second-order expansions of state variables at different Gauss points with finite distances, this effect could be expected. Nevertheless, with increasing height of the macroscopic structure the relative distance between the RVEs with respect to the macroscopic dimensions decreases, so that the jumps become less and less severe for both, the first- and second-order models.

The impact of the higher-order expansions becomes even more evident in Fig. 6, where we observe large discrepancies between the first-order results on the one hand and the second-order and DNS results on the other, the latter two being in excellent agreement. Note that the plotted values of the mass flux vector are the projected values in both multiscale and DNS solutions.

### 4.1.2 Multiscale behavior of fixed-sized macroscopic bodies with differently sized RVEs

In the following study we consider a similar setup as in the preceding section, the only change being with respect to the dimensions of the macroscopic and microscopic domains. While in the previous study we analyzed macroscopic domains that were assembled from different numbers of RVEs of same size, we now investigate the behavior of a macroscopic domain of fixed size ($\bar{l}_1 \times \bar{l}_2$) = (2 × 2) m² assembled from quadrilateral RVEs of different dimensions. The dimensions of the considered three different RVEs (labeled with 'A', 'B' and 'C') are given by ($l_1^A \times l_2^A$) = (1 × 1) m², ($l_1^B \times l_2^B$) = (1/2 × 1/2) m² and ($l_1^C \times l_2^C$) = (1/3 × 1/3) m². The results are plotted in Fig. 7. Again, we observe a good agreement of the DNS and the multiscale simulations for the different RVE sizes.
<figure><img src="image_9.png" /><figcaption>**Figure 7:** Profiles of the vertical component of the fluid-mass flux for particulate microstructures with 4, 16 and 36 inclusions. The dimensions of the macroscopic domain is given as $(\bar{l}_1 \times \bar{l}_1) = (2 \times 2) \text{ m}^2$ such that the size of the periodic RVEs of each microstructure decrease with the number of REVs. The profiles of the fluid-mass flux have been extracted from the material points located on the dashed lines drawn on the inset microstructure images.</figcaption></figure>

### 4.1.3 Three-phase particulate microstructure

We now increase the complexity of the problem with respect to the macroscopic boundary conditions and the considered microstructure. As macroscopic boundary conditions we prescribe

$$ \partial \mathcal{B}_{\text{left}} : \bar{h}_1(\mathbf{x}, t) = 2 \times 10^{-3} \text{ kg/(m}^2 \text{ s)}, \qquad (149) $$

$$ \partial \mathcal{B}_{\text{right}} : \bar{h}_1(\mathbf{x}, t) = 0, \qquad (150) $$

$$ \partial \mathcal{B}_{\text{top}} : \bar{h}_2(\mathbf{x}, t) = -1 \times 10^{-3} \text{ kg/(m}^2 \text{ s)}, \qquad (151) $$

$$ \partial \mathcal{B}_{\text{bottom}} : \bar{h}_2(\mathbf{x}, t) = 0, \qquad (152) $$

$$ \partial \mathcal{B} : \bar{\beta}_{12}(\mathbf{x}, t) = \bar{\beta}_{21}(\mathbf{x}, t) = 0, \qquad (153) $$

from which the last one is only relevant for the second-order FE² simulations. With regard to the microstructure, we distribute two types of RVEs across the macroscopic domain as indicated in Fig. 8. The dimensions of the RVEs are set to $(l_1 \times l_2) = (\frac{1}{3} \times \frac{1}{3}) \text{ m}^2$, resulting in a macroscopic domain size of $(\bar{l}_1 \times \bar{l}_2) = (2 \times 2) \text{ m}^2$. The material parameters are listed in Table 2.

As indicated by the microstructure depicted in Fig. 8, we investigate the components of the fluid-mass flux obtained from the multiscale simulations and the DNS along a diagonal section of the microstructure. We observe that both simulation results are in good agreement with each other despite the complexity of the morphology and the loading conditions. For further information, we provide contour plots of the fluid-mass flux at two time instances in Fig. 9.
<figure><img src="image_10.png" /><figcaption>**Figure 8:** Profiles of the horizontal and vertical components of the fluid-mass flux for particulate microstructures composed of a matrix phase and two kinds of inclusions with different permeability. The dimensions of the macroscopic domain are given as $(\bar{l}_1 \times \bar{l}_2) = (2 \times 2)$ m$^2$ and the size of the periodic RVEs are $(l_1 \times l_2) = (\frac{1}{3} \times \frac{1}{3})$ m$^2$. The profiles of the fluid-mass flux have been extracted from the material points located on the dashed lines drawn on the inset microstructure images.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9:** Spatial distributions of the fluid-mass flux in the three-phase microstructure at time instances $t = 1$ s (left) and $t = 15$ s (right). The vectors of the fluid-mass flux are scaled by their magnitude, which is provided in the underlying contour plot.</figcaption></figure>
4.1.4 Non-centrosymmetric particulate microstructures

All microstructures considered thus far have been centrosymmetric. To further evaluate the performance
of the proposed multiscale model, we will now apply it to non-centrosymmetric RVEs (Yvonnet, Auffray,
and Monchiet, 2020). While we select the same boundary and initial conditions as in Section 4.1.1,
we further reduce the permeability of the inclusions (Table 3). The center of the inclusions is shifted
by **x**<sub>rel</sub> = [-0.15, 0.15]<sup>T</sup> m with respect to the center of each RVE. The radius of the inclusions is
r = 0.25 m.

In Fig. 10, we present the results obtained from FE² simulations and DNS microstructures made of (4×4) and (8×8) RVEs with dimensions ($l_1 \times l_2$) = (1×1) m². We observe that the vertical component of the fluid-mass flux are in good agreement. However, some discrepancies are observed near the upper and lower boundaries and at the interfaces between neighboring RVEs. These discrepancies are larger than in case of centrosymmetric microstructures.

**4.1.5 Darcy–Forchheimer fluid flow through anisotropic particulate microstructures**

We now consider a material that shows nonlinear flow behavior such that it exhibits a nonlinear relation
between the fluid-mass flux and the fluid-potential gradient. Such a behavior can be described by the
Darcy–Forchheimer law, which is imprinted in the dissipation-potential-density function provided in
Eq. 143. The anisotropic microstructure of the material is composed of matrix and embedded elliptic
inclusions with a volume fraction of 20 % and an aspect ratio of 3 to 1. The boundary conditions are
again the same as in Section 4.1.1 and the material parameters are given in Table 4.

Fig. 11 shows the results obtained from second-order FE² simulations and DNS for a microstructure composed of (6×6) RVEs. There we analyze the components of the fluid-mass flux along the horizontal and vertical coordinate axes at the time instances t = 4 s and t = 15 s. As becomes clear from the nonlinear macroscopic curves, neither of the results corresponds to a steady state.

A stark difference with respect to the previous results lies in the presence of a non-zero horizontal
component of the macroscopic fluid-mass flux, which is due to the anisotropy of the microstructure.

This time, we observe larger discrepancies between the fluid-mass flux predicted by the DNS and the
second-order FE² simulations. These discrepancies do however decrease once the solution converges
to the steady state. In this context we would furthermore like to note that the locations of the RVEs of
the FE² simulations do not align with the microstructure of the DNS, which is again due to the shifted
locations of the Gauss points (Fig. 4). As discussed earlier, this mismatch naturally leads to deviations
of results and these deviations will be more pronounced in case of anisotropic microstructures that
show a significant coupling of the response along different spatial directions.

In Fig. 11 we have also plotted results of the first-order multiscale simulations. We observe that the
fluid-mass flux deviates significantly from the DNS at early time steps. However, once the solution
converges to the steady state, the discrepancies become smaller. A similar effect could be observed in
earlier studies based on RVE-level simulations (Polukhov and Keip, 2020).
<figure><img src="image_12.png" /><figcaption>**Figure 10:** Profiles of the vertical component of the fluid-mass flux for non-centrosymmetric particulate microstructures with 16 and 64 inclusions. The dimensions of the periodic RVEs of each microstructure are given as ($l_1 \times l_2$) = (1 $\times$ 1) m$^2$ such that the height of the macroscopic domain increases with the number of REVs. The profiles of the fluid-mass flux have been extracted from the material points located on the dashed lines drawn on the inset microstructure images, passing either through the centers of the inclusions (a and c) or the centers of the RVEs (b and d).</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 11:** Profiles of the horizontal and vertical components of the fluid-mass flux for anisotropic particulate microstructures with 36 inclusions exhibiting Darcy–Forchheimer fluid flow. The dimensions of the periodic RVEs of each microstructure are given as ($l_1 \times l_2$) = (1 × 1) m² and the major axis of each inclusion is oriented by 45° with respect to the horizontal axis. The profiles of the fluid-mass flux have been extracted from the material points located on the dashed lines drawn on the inset microstructure images.</figcaption></figure>
### 4.1.6 Rate-dependence of macroscopic response

In the preceding sections, we examined the performance of the second-order computational-homogenization approach for various microstructures exhibiting linear and nonlinear material responses. Furthermore, we investigated size effects for different boundary-value problems. In the present section, we will move on and investigate the impact of loading rate on the macroscopic response.

Since we are dealing with transient problems, we expect the loading rate to have a substantial impact on the effective response; This impact should also depend on the RVE size (Polukhov and Keip, 2020). In the latter work, we analyzed the effective response of RVEs of different size under fast loading rates. There, we observed pronounced non-equilibrium effective responses that eventually converged to a common solution (Polukhov and Keip, 2020, Fig. 9).

To elucidate the rate-dependent behavior in the context of the second-order multiscale model, we consider a set of FE² simulations based on a homogeneous RVE (Fig. 12). The macroscopic domain has the dimensions ($\overline{l}_1 \times \overline{l}_2$) = (2 × 4) m² and is exposed to the boundary conditions

$$ \partial \mathcal{B}_{\text{left}} : \bar{h}_1(\mathbf{x}, t) = 0, \qquad (154) $$

$$ \partial \mathcal{B}_{\text{right}} : \bar{h}_1(\mathbf{x}, t) = 0, \qquad (155) $$

$$ \partial \mathcal{B}_{\text{top}} : \bar{h}_2(\mathbf{x}, t) = -1 \times 10^{-3} \zeta \text{ kg/(m}^2\text{s)}, \qquad (156) $$

$$ \partial \mathcal{B}_{\text{bottom}} : \bar{h}_2(\mathbf{x}, t) = 0, \qquad (157) $$

$$ \partial \mathcal{B}_{\text{bottom}} : \bar{\beta}_{22}(\mathbf{x}, t) = 0, \qquad (158) $$

$$ \partial \mathcal{B} : \bar{\beta}_{12}(\mathbf{x}, t) = \bar{\beta}_{21}(\mathbf{x}, t) = 0, \qquad (159) $$

where $\zeta := \hat{\zeta}(t) \in [0, 1] \subset \mathbb{R}$ is a ramp function that characterizes the loading profile such that

$$ \hat{\zeta}(t) = \begin{cases} \frac{1}{t_{\text{ramp}}} t & \text{for } t < t_{\text{ramp}} \\ 1 & \text{for } t \ge t_{\text{ramp}} \end{cases} \qquad (160) $$

In what follows, we consider three different loading rates by setting (i) $t_{\text{ramp}} = 1$ s (Fig. 12a,b), (ii) $t_{\text{ramp}} = 10$ s (Fig. 12c,d), and (iii) $t_{\text{ramp}} = 50$ s (Fig. 12e,f). In all simulations the time-step size is set to $\tau = 1$ s.

The results of the effective fluid-mass flux plotted in Fig. 12 show that the size of an RVE has a vanishing impact when the loading rate is low. Conversely, when the loading rate is high, the size of an RVE plays a major role. Such a behavior could be expected because of the transient nature of the microscopic problem as well as the second-order nature of the homogenization procedure.<fnref n="16" /> We further

<footnote n="16">¹⁶We would like to note explicitly that homogeneous RVEs do not lead to a macroscopic Cauchy continuum. This behavior is in contradiction with the considerations of Yuan, Tomita, and Andou (2008) and Forest and Trinh (2011) who argue that homogeneous RVEs should be associated with vanishing macroscopic hyper fields. While this argumentation appears plausible, we point out that Kouznetsova, Geers, and Brekelmans (2004b) could show an equivalence between a macroscopic model obtained from second-order homogenization over homogeneous, linear-elastic RVEs and Mindlin's</footnote>
observe that the rate- and size-dependence of the effective response decreases with decreasing sizes of RVEs, eventually resulting in a Cauchy response associated with the separation of length scales (Polukhov and Keip, 2020).

## 4.2 Fluid flow through porous media in three spatial dimensions

Before we apply the proposed multiscale model to complex three-dimensional scenarios, we validate the underlying three-dimensional FE² implementation by revisiting the Darcy and Darcy–Forchheimer models for fluid flow through porous media documented in Sections 4.1.1 and 4.1.5.

### 4.2.1 Validation of the three-dimensional implementation

In order to validate the three-dimensional implementation of the proposed formulation, we compare results of three-dimensional second-order FE² simulations to corresponding results of second-order FE² and direct numerical simulations in two dimensions. As benchmark problems, we consider the boundary-value problems documented in Sections 4.1.1 and 4.1.5 (see also the inset microstructures in Figs. 13 and 14). For the three-dimensional simulations we assume the existence of periodic patterns of fibers. The fibers have the same cross section as the inclusions of the two-dimensional problems and are oriented in the direction of the $x_3$ axis; We can thus map them into a corresponding pattern of RVEs whose projection on the $x_1-x_2$ plane resembles the two-dimensional microstructure. Furthermore, we assume the macroscopic boundary conditions

$$ \partial \bar{\mathcal{B}}_{\text{left}} : \bar{h}_1(\mathbf{x}, t) = 0, \qquad (161) $$

$$ \partial \bar{\mathcal{B}}_{\text{right}} : \bar{h}_1(\mathbf{x}, t) = 0, \qquad (162) $$

$$ \partial \bar{\mathcal{B}}_{\text{top}} : \bar{h}_2(\mathbf{x}, t) = -1 \times 10^{-3} \text{ kg/(m}^2 \text{ s)}, \qquad (163) $$

$$ \partial \bar{\mathcal{B}}_{\text{bottom}} : \bar{h}_2(\mathbf{x}, t) = 0, \qquad (164) $$

$$ \partial \bar{\mathcal{B}}_{\text{front}} : \bar{h}_3(\mathbf{x}, t) = 0, \qquad (165) $$

$$ \partial \bar{\mathcal{B}}_{\text{back}} : \bar{h}_3(\mathbf{x}, t) = 0, \qquad (166) $$

$$ \partial \bar{\mathcal{B}} : \bar{\beta}_{ij}(\mathbf{x}, t) = 0 \text{ for } i \neq j. \qquad (167) $$

As we can see from Figs. 13 and 14, the results of the three-dimensional simulations are in line with the results obtained from two-dimensional FE² and direct numerical simulations for both, linear and nonlinear problems.

<footnote>strain-gradient by establishing a link between RVE size and the length-scale parameter of an underlying strain-energy function.</footnote>
<figure><img src="image_14.png" /><figcaption>**Figure 12:** Profiles of the vertical component of the fluid-mass flux for homogeneous microstructures under different macroscopic loading rates according to Eqs. 154 to 159. The macroscopic bodies with dimensions ($\bar{l}_1 \times \bar{l}_2$) = (2 × 4) m² (shown in the sub-figures a to f) are composed of homogeneous RVEs (labeled with 'A', 'B' and 'C') with dimensions ($l_1^A \times l_2^A$) = (0.5 × 0.5) m², ($l_1^B \times l_2^B$) = (1 × 1) m² and ($l_1^C \times l_2^C$) = (2 × 2) m². The profiles of the fluid-mass flux have been extracted from the material points located on the dashed lines drawn on the inset microstructure images, where the profiles in the left column correspond to solutions at time instance $t = 1$ s and those in the right column to $t = 10$ s.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure 13:** Profiles of the vertical component of the fluid-mass flux for particulate microstructures with 24 inclusions predicted by two- and three-dimensional second-order FE² and direct numerical simulations in consideration of linear Darcy flow. The dimensions of the periodic RVEs of each microstructure are given as ($l_1 \times l_2$) = (1 × 1) m² such that the dimensions of the macroscopic domain are ($\bar{l}_1 \times \bar{l}_2$) = (2 × 12) m². The profiles of the fluid-mass flux have been extracted from the material points located on the dashed lines drawn on the inset microstructure images.</figcaption></figure>

<table><caption>**Table 5:** Values of material parameters for the two-phase microstructures in Figs. 15 and 16.</caption><thead><tr><th>Phases</th><th>A<br>m<sup>5</sup>/(kg s<sup>2</sup>)<br>(47), (142)</th><th>K<br>m<sup>3</sup> s/kg<br>(47), (142)</th><th>&kappa;<sub>&mu;</sub><br>m<sup>5</sup>/(kg s)<br>(47)</th><th>&kappa;<sup>&tau;</sup><sub>&lambda;</sub><br>m<sup>3</sup>/kg<br>(46)</th><th>&kappa;<sup>&tau;</sup><sub>&eta;</sub><br>m<sup>3</sup>/kg<br>(46)</th><th>&kappa;<sub>&beta;</sub><br>m<sup>5</sup>/(kg s)<br>(128)</th><th>&kappa;<sub>&sigma;</sub><br>m<sup>5</sup>/(kg s)<br>(128)</th><th>&gamma;<br>m<sup>5</sup>/kg<sup>2</sup><br>(128)</th></tr></thead><tbody><tr><td>Matrix</td><td>10</td><td>10<sup>-1</sup></td><td>10<sup>-14</sup></td><td>10<sup>15</sup></td><td>10<sup>12</sup></td><td>10<sup>4</sup></td><td>10<sup>12</sup></td><td>20</td></tr><tr><td>Inclusions</td><td>10</td><td>10<sup>-3</sup></td><td>10<sup>-14</sup></td><td>10<sup>15</sup></td><td>10<sup>12</sup></td><td>10<sup>4</sup></td><td>10<sup>12</sup></td><td>5</td></tr></tbody></table>

### 4.2.2 Darcy–Forchheimer fluid flow through particulate and triply-periodic minimal-surface microstructures

We simulate the response of three-dimensional macroscopic bodies with particulate and triply-periodic minimal-surface (TPMS) microstructures. The macroscopic bodies have the dimensions ($\bar{l}_1 \times \bar{l}_2 \times \bar{l}_3$) = (0.2 × 0.2 × 1.2) m³ and the three-dimensional RVEs have cubic shapes with dimensions ($l_1 \times l_2 \times l_3$) = (0.1 × 0.1 × 0.1) m³. Both kinds of RVEs are assumed to be composed of porous phases. In case of the particulate microstructure, we consider a porous matrix and porous spherical inclusions, each with 50% volume fraction<fnref n="17" />; In case of the TPMS microstructure we

<footnote n="17">¹⁷The particulate microstructure contains six randomly distributed spheres, half of which have the radius $r = 0.015$ m and the remaining half has the radius $r = 0.025$ m.</footnote>
<figure><img src="image_16.png" /><figcaption>**Figure 14:** Profiles of the vertical component of the fluid-mass flux for particulate microstructures with 36 inclusions predicted by two- and three-dimensional second-order FE² and direct numerical simulations in consideration of nonlinear Darcy–Forchheimer flow. The dimensions of the periodic RVEs of each microstructure are given as ($l_1 \times l_2$) = (1 × 1) m² such that the dimensions of the macroscopic domain are ($\bar{l}_1 \times \bar{l}_2$) = (6 × 6) m². The profiles of the fluid-mass flux have been extracted from the material points located on the dashed lines drawn on the inset microstructure images.</figcaption></figure>

consider Schoen’s gyroid with phase-volume fractions of 24% and 76%, respectively (Fig. 15).<fnref n="18" /> In both cases we assume the material parameters listed in Table 5. The macroscopic bodies are discretized with ($\bar{n}_1^e \times \bar{n}_2^e \times \bar{n}_3^e$) = (1 × 1 × 6) finite elements along the individual spatial directions and exposed to the boundary conditions

$$ \partial \bar{\mathcal{B}}_{\text{left}} : \bar{h}_1(\mathbf{x}, t) = 0, \qquad (168) $$

$$ \partial \bar{\mathcal{B}}_{\text{left}} : \bar{\beta}_{11}(\mathbf{x}, t) = 0, \qquad (169) $$

$$ \partial \bar{\mathcal{B}}_{\text{right}} : \bar{h}_1(\mathbf{x}, t) = 0, \qquad (170) $$

$$ \partial \bar{\mathcal{B}}_{\text{right}} : \bar{\beta}_{11}(\mathbf{x}, t) = 0, \qquad (171) $$

$$ \partial \bar{\mathcal{B}}_{\text{top}} : \bar{h}_3(\mathbf{x}, t) = -1 \times 10^{-3} \, \text{kg/(m}^2 \text{ s)}, \qquad (172) $$

$$ \partial \bar{\mathcal{B}}_{\text{bottom}} : \bar{h}(\mathbf{x}, t) = 0, \qquad (173) $$

$$ \partial \bar{\mathcal{B}}_{\text{top}} : \bar{\beta}_{12}(\mathbf{x}, t) = \bar{\beta}_{21}(\mathbf{x}, t) = \bar{\beta}_{33}(\mathbf{x}, t) = -1 \times 10^{-5} \, \text{kg/(m}^3 \text{ s)}, \qquad (174) $$

$$ \partial \bar{\mathcal{B}}_{\text{bottom}} : \bar{\beta}_{12}(\mathbf{x}, t) = \bar{\beta}_{21}(\mathbf{x}, t) = \bar{\beta}_{33}(\mathbf{x}, t) = 1 \times 10^{-5} \, \text{kg/(m}^3 \text{ s)}, \qquad (175) $$

<footnote n="18">¹⁸The RVEs have been generated by using the open-source Python libraries provided by Chemisky and Pruliere (2024) and Metsch, Hansen-Dörr, and Turan (2024).</footnote>
<figure><img src="image_17.png" /><figcaption>**Figure 15:** Profiles of the vertical component of the fluid-mass flux for particulate and TPMS microstructures exhibiting Darcy–Forchheimer fluid flow obtained from second-order FE² simulations. The macroscopic bodies with dimensions ($\bar{l}_1 \times \bar{l}_2 \times \bar{l}_3$) = (0.2 × 0.2 × 1.2) m³ are composed of RVEs with dimensions ($l_1^A \times l_2^A \times l_3^A$) = (0.1 × 0.1 × 0.1) m³. The profiles of the fluid-mass-flux have been extracted along the $\bar{x}_3$-axis from the material points located at $(\bar{x}_1, \bar{x}_2) = (0.0, 0.0)$ m.</figcaption></figure>

$$ \partial \bar{\mathcal{B}}_{\text{front}} : \bar{h}_2(\mathbf{x}, t) = 0, \qquad (176) $$

$$ \partial \bar{\mathcal{B}}_{\text{front}} : \bar{\beta}_{22}(\mathbf{x}, t) = 0, \qquad (177) $$

$$ \partial \bar{\mathcal{B}}_{\text{back}} : \bar{h}_2(\mathbf{x}, t) = 0, \qquad (178) $$

$$ \partial \bar{\mathcal{B}}_{\text{back}} : \bar{\beta}_{22}(\mathbf{x}, t) = 0. \qquad (179) $$

Fig. 15 shows the resulting macroscopic response. We observe an intuitive response for both microstructures. Because of the lower volume fraction of the less permeable phase in the particulate microstructure, the system reaches the steady state earlier than the gyroid. The corresponding streamlines on the selected RVEs are shown in Fig. 16.
<figure><img src="image_18.png" /><figcaption>**Figure 16:** Spatial distributions of the fluid-mass flux in the particulate (top row) and gyroid (bottom row) microstructures at different time instances. The distribution of the fluid-mass flux is shown in the form of streamlines, where the shift from blue to red represents qualitatively increase in the magnitude of fluid-mass flux. We plot the RVEs located at the macroscopic Gauss point with coordinates $\bar{x} = [0.05, 0.05, 0.45]^T$ m.</figcaption></figure>

# 5 Concluding remarks

We presented a second-order multiscale model of fluid flow through porous media based on macroscopic and microscopic variational minimization formulations with the fluid-mass flux as primary unknown field. To facilitate the numerical implementation associated with the increased continuity requirements at both microscopic and macroscopic scale, we proposed Hu-Washizu variational principles that enabled us to incorporate Darcy and Darcy-Forchheimer flow in a straightforward manner. We demonstrated the predictive capabilities of the formulation by means of a series of numerical examples in two- and three-dimensional domains and comparisons with direct numerical simulations. In particular, we investigated transient fluid-flow processes in two- and three-phase porous media with centrosymmetric and non-centrosymmetric as well as homogeneous and heterogeneous microstructures. We demonstrated that the size of the representative volume element (RVE) plays a significant role in the effective response of a porous medium, even when a homogeneous RVE is considered. These features make the proposed formulation useful for multiscale simulations in which a clear separation of length scales is not present
and could be exploited for the parameter identification of higher-order phenomenological models of
fluid flow through porous media along the lines of Sciarra, dell'Isola, and Coussy (2007) and Khurshid,
Polukhov, and Keip (2024).

Acknowledgements

We gratefully acknowledge the financial support of the German Research Foundation (Deutsche
Forschungsgemeinschaft; DFG) in the framework of the Collaborative Research Center 1313 *Interface-
Driven Multi-Field Processes in Porous Media—Flow, Transport, and Deformation* at the University of
Stuttgart (Project Number 327154368, Project B01).

A Appendix

Convergence study of the macroscopic finite-element discretization with fixed
size of RVEs

As was discussed in Section 4, a material with microstructure could either be simulated via direct
numerical simulations (which resolve the entire microstructure altogether) or by means of multiscale
simulations (which assume a homogeneous macroscopic body with attached RVEs that represent the
microstructure). While first-order homogenization schemes are not sensitive with respect to the size
of RVEs, second-order schemes are. In this appendix, we want to study the impact of a macroscopic
finite-element discretization on the solution of a homogeneous, second-order multiscale boundary-value
problem in consideration of heterogeneous RVEs with fixed size. While the size of each RVE depends on
the given microstructure and is thus a physically motivated parameter, the finite-element discretization
of a domain is arbitrary. It could thus happen that the domain of integration of a macroscopic Gauss
point is smaller than an RVE. To check whether or not such a circumstance has an impact of the effective
response we consider different finite-element discretizations of two macroscopic bodies (labeled with
'A' and 'B') with fixed dimensions ($\overline{l_1^A} \times \overline{l_2^A}$) = (2 × 4) m² and ($\overline{l_1^B} \times \overline{l_2^B}$) = (2 × 2) m², respectively.
Both bodies have microstructures described by RVEs of dimensions ($l_1^A \times l_2^A$) = (1 × 1) m² and
($l_1^B \times l_2^B$) = (0.5 × 0.5) m², respectively, where the RVE 'A' is modeled using the linear Darcy law and
RVE 'B' is modeled using the nonlinear Darcy-Forchheimer law. Fig. 17 shows corresponding results
of the multiscale simulations. As becomes clear, the size of macroscopic finite elements (and therewith
the size of integration domain per macroscopic Gauss point) does not have an impact on the effective
response. In fact, macroscopic finite elements could even be smaller than an RVE.
<figure><img src="image_19.png" /><figcaption>**Figure 17:** *Profiles of the vertical component of the effective fluid-mass flux $\bar{h}_2$ in consideration of a homogeneous macroscopic problem with attached RVEs as shown in the figures.* The dimensions of the macroscopic domains are (a) $(\overline{l}_1^A \times \overline{l}_2^A) = (2 \times 4) \text{ m}^2$ and (b) $(\overline{l}_1^B \times \overline{l}_2^B) = (2 \times 2) \text{ m}^2$ periodic RVEs are (a) $(l_1^A \times l_2^A) = (1 \times 1) \text{ m}^2$ and (b) $(l_1^B \times l_2^B) = (0.5 \times 0.5) \text{ m}^2$. The profiles of the fluid-mass flux have been extracted from the material points along the vertical axis at the center of macroscopic domain.</figcaption></figure>

## References

Altan, B. and E. Aifantis (1997). "On some aspects in the special theory of gradient elasticity". In: *Journal of the Mechanical Behavior of Materials* 8.3, pp. 231–282.

Anderson, E., Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen (1999). *LAPACK users' guide*. 3rd ed. Philadelphia, PA: Society for Industrial and Applied Mathematics.

Anonis, R. A., J. L. Mroginski, and P. J. Sánchez (2024). "Multiscale formulation for saturated porous media preserving the representative volume element size objectivity". In: *International Journal for Numerical Methods in Engineering* 125.3, e7381.

Auricchio, F., L. B. da Veiga, F. Brezzi, and C. Lovadina (2017). "Mixed finite element methods". In: *Encyclopedia of computational mechanics* second edition, pp. 1–53.

Aznaran, F. R., P. E. Farrell, and R. C. Kirby (2022). "Transformations for piola-mapped elements". In: *The SMAI Journal of computational mathematics* 8, pp. 399–437.

Barboura, S. and J. Li (2018). "Establishment of strain gradient constitutive relations by using asymptotic analysis and the finite element method for complex periodic microstructures". In: *International Journal of Solids and Structures* 136, pp. 60–76.

Biot, M. A. (1941). "General theory of three-dimensional consolidation". In: *Journal of Applied Physics* 12.2, pp. 155–164.
Blanco, P. J., P. J. Sánchez, E. A. De Souza Neto, and R. A. Feijóo (2016). “The method of multiscale virtual power for the derivation of a second order mechanical model”. In: *Mechanics of Materials* 99, pp. 53–67.

Brezzi, F. and M. Fortin (1991). *Mixed and hybrid finite element methods*. Springer-Verlag.

Brink, U. and E. Stein (1996). “On some mixed finite element methods for incompressible and nearly incompressible finite elasticity”. In: *Computational Mechanics* 19.1, pp. 105–119.

Brown, D. L., Y. Efendiev, G. Li, and V. Savatorova (2015). “Homogenization of high-contrast Brinkman flows”. In: *Multiscale Modeling & Simulation* 13.2, pp. 472–490.

Byrne, H. and L. Preziosi (2003). “Modelling solid tumour growth using the theory of mixtures”. In: *Mathematical medicine and biology: a journal of the IMA* 20.4, pp. 341–366.

Carvalho, P. G., P. R. Devloo, and S. M. Gomes (2024). “A semi-hybrid-mixed method for Stokes–Brinkman–Darcy flows with H (div)-velocity fields”. In: *International Journal for Numerical Methods in Engineering* 125.1, e7363.

Chemisky, Y. and E. Pruliere (2024). *3MAH/Microgen – Python library for microstructure generation and meshing*. URL: https://github.com/3MAH/microgen.

Coussy, O. (2004). *Poromechanics*. John Wiley & Sons.

Darcy, H. (1856). *Les fontaines publiques de la ville de Dijon: Exposition et application des principes a suivre et des formules a employer dans les questions de distribution d'eau; ouvrage terminé par un appendice relatif aux fournitures d'eau de plusieurs villes au filtrage des eaux et a la fabrication des tuyaux de fonte, de plomb, de tole et de bitume*. Paris: Victor Dalmont, Libraire des Corps imperiaux des ponts et chaussées et des mines.

Ehlers, W. (2018). *Vector and tensor calculus: An introduction*. Lecture Notes at the University of Stuttgart.

Ehlers, W. (2022). “Darcy, Forchheimer, Brinkman and Richards: classical hydromechanical equations and their significance in the light of the TPM”. In: *Archive of Applied Mechanics* 92.2, pp. 619–639.
— (2023). “A historical review on porous-media research”. In: *PAMM* 23.3, e202300271.

Forest, S. and R. Sievert (2006). “Nonlinear microstrain theories”. In: *International journal of solids and structures* 43.24, pp. 7224–7245.

Forest, S. and D. K. Trinh (2011). “Generalized continua and non-homogeneous boundary conditions in homogenisation methods”. In: *ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik* 91.2, pp. 90–109.

Gurtin, M. E., E. Fried, and L. Anand (2010). *The mechanics and thermodynamics of continua*. Cambridge University Press.

Hesch, C., F. Schmidt, and S. Schuß (2024). “Variational formulation and monolithic solution of computational homogenization methods”. In: *International Journal for Numerical Methods in Engineering* 125.20, e7567.

Kaczmarczyk, Ł., C. J. Pearce, and N. Bićanić (2008). “Scale transition and enforcement of RVE boundary conditions in second-order computational homogenization”. In: *International Journal for Numerical Methods in Engineering* 74.3, pp. 506–522.
Kaessmair, S. and P. Steinmann (2018). "Computational first-order homogenization in chemo-mechanics". In: *Archive of Applied Mechanics* 88.1, pp. 271–286.

Khurshid, H., E. Polukhov, and M.-A. Keip (2024). "Mixed variational formulation and finite-element implementation of second-order poro-elasticity". In: *International Journal of Solids and Structures* 305, p. 113055.

Klahr, B., J. L. M. Thiesen, O. T. Pinto, T. A. Carniel, and E. A. Fancello (2023). "A variational RVE-based multiscale poromechanical formulation applied to soft biological tissues under large deformations". In: *European Journal of Mechanics-A/Solids* 99, p. 104937.

Kouznetsova, V. G., M. G. Geers, and W. Brekelmans (2004a). "Multi-scale second-order computational homogenization of multi-phase materials: a nested finite element solution strategy". In: *Computer Methods in Applied Mechanics and Engineering* 193.48, pp. 5525–5550.

— (2004b). "Size of a representative volume element in a second-order computational homogenization framework". In: *International Journal for Multiscale Computational Engineering* 2.4.

Kumari, W., P. Ranjith, M. Perera, X. Li, L. Li, B. Chen, B. A. Isaka, and V. De Silva (2018). "Hydraulic fracturing under high temperature and pressure conditions with micro CT applications: geothermal energy from hot dry rocks". In: *Fuel* 230, pp. 138–154.

Lange, N., G. Hütter, and B. Kiefer (2021). "An efficient monolithic solution scheme for FE2 problems". In: *Computer Methods in Applied Mechanics and Engineering* 382, p. 113886.

Liupekevicius, R., J. A. van Dommelen, M. G. Geers, and V. G. Kouznetsova (2024). "Transient computational homogenization of heterogeneous poroelastic media with local resonances". In: *International Journal for Numerical Methods in Engineering* 125.18, e7505.

Lopes, I. A. R. and F. M. A. Pires (2022a). "An assessment of multi-scale models based on second-order computational homogenisation". In: *Computers & Structures* 259, p. 106679.

— (2022b). "Formulation and numerical implementation of a variationally consistent multi-scale model based on second-order computational homogenisation at finite strains for quasi-static problems". In: *Computer Methods in Applied Mechanics and Engineering* 392, p. 114714.

Luscher, D. J., D. L. McDowell, and C. A. Bronkhorst (2010). "A second gradient theoretical framework for hierarchical multiscale modeling of materials". In: *International Journal of Plasticity* 26.8, pp. 1248–1275.

Maike, S., J. Schröder, J. Bluhm, and T. Ricken (2024). "A mesh-in-element method for the theory of porous media". In: *International Journal for Numerical Methods in Engineering* 125.21, e7565.

Mercer, J. W. and R. M. Cohen (1990). "A review of immiscible fluids in the subsurface: Properties, models, characterization and remediation". In: *Journal of contaminant hydrology* 6.2, pp. 107–163.

Metsch, P., A. C. Hansen-Dörr, and O. T. Turan (2024). *gmshModel (Version 1.1.1) [Computer software] – A mesh modeling interface to the Gmsh-Python-API.* URL: https://doi.org/10.5281/zenodo.13353550.

Miehe, C., J. Schotte, and J. Schröder (1999). "Computational micro-macro transitions and overall moduli in the analysis of polycrystals at large strains". In: *Computational Materials Science* 16, pp. 372–382.
Miehe, C. (1994). "Aspects of the formulation and finite element implementation of large strain isotropic elasticity". In: *International journal for numerical methods in engineering* 37.12, pp. 1981–2004.

Miehe, C., S. Mauthe, and S. Teichtmeister (2015). "Minimization principles for the coupled problem of Darcy–Biot-type fluid transport in porous media linked to phase field modeling of fracture". In: *Journal of the Mechanics and Physics of Solids* 82, pp. 186–217.

Mindlin, R. D. (1964). "Micro-structure in linear elasticity". In: *Archive for rational mechanics and analysis* 16, pp. 51–78.

Monchiet, V. (2024). "Homogenization of non-Darcy flow in porous media containing spheroidal impervious inclusions". In: *International Journal of Solids and Structures*, p. 112867.

Monchiet, V. (2023). "Non-linear homogenization of Forchheimer law for porous media with cylindrical and spherical impervious inclusions". In: *Acta Mechanica* 234.12, pp. 6607–6627.

Oh, K.-Y., J. B. Siegel, L. Secondo, S. U. Kim, N. A. Samad, J. Qin, D. Anderson, K. Garikipati, A. Knobloch, B. I. Epureanu, et al. (2014). "Rate dependence of swelling in lithium-ion cells". In: *Journal of Power Sources* 267, pp. 197–202.

OpenMP, A. (2011). "Openmp application program interface v5.2". In: *OpenMP Architecture Review Board*, p. 649.

Pałka, K. and R. Pokrowiecki (2018). "Porous titanium implants: a review". In: *Advanced Engineering Materials* 20.5, p. 1700648.

Polizzotto, C. (2003). "Gradient elasticity and nonstandard boundary conditions". In: *International Journal of Solids and Structures* 40.26, pp. 7399–7423.

Pollmann, N., F. Larsson, K. Runesson, and R. Jänicke (2020). "Diffuse interface modeling and variationally consistent homogenization of fluid transport in fractured porous media". In: *European Journal of Mechanics-A/Solids* 84, p. 104067.

Polukhov, E. and M.-A. Keip (2020). "Computational homogenization of transient chemo-mechanical processes based on a variational minimization principle". In: *Advanced Modeling and Simulation in Engineering Sciences* 7.1, pp. 1–26.

Ricken, T., J. Schroeder, J. Bluhm, S. Maike, and F. Bartel (2022). "Theoretical formulation and computational aspects of a two-scale homogenization scheme combining the TPM and FE2 method for poro-elastic fluid-saturated porous media". In: *International Journal of Solids and Structures* 241, p. 111412.

Schenk, O., K. Gärtner, W. Fichtner, and A. Stricker (2001). "PARDISO: a high-performance serial and parallel sparse linear solver in semiconductor device simulation". In: *Future Generation Computer Systems* 18.1, pp. 69–78.

Sciarra, G., F. dell'Isola, and O. Coussy (2007). "Second gradient poromechanics". In: *International Journal of Solids and Structures* 44.20, pp. 6607–6629.

Sky, A., M. Neunteufel, J. S. Hale, and A. Zilian (2023). "A Reissner–Mindlin plate formulation using symmetric Hu–Zhang elements via polytopal transformations". In: *Computer Methods in Applied Mechanics and Engineering* 416, p. 116291.
Su, F., F. Larsson, and K. Runesson (2011). "Computational homogenization of coupled consolidation problems in micro-heterogeneous porous media". In: *International journal for numerical methods in engineering* 88.11, pp. 1198-1218.

Tan, V. B. C., K. Raju, and H. P. Lee (2020). "Direct FE2 for concurrent multilevel modelling of heterogeneous structures". In: *Computer Methods in Applied Mechanics and Engineering* 360, p. 112694.

Taylor, R. L. and S. Govindjee (2021). *FEAP – A finite element analysis program*. URL: http://www.ce.berkeley/feap.

Teichtmeister, S., S. Mauthe, and C. Miehe (2019). "Aspects of finite element formulations for the coupled problem of poroelasticity based on a canonical minimization principle". In: *Computational Mechanics* 64, pp. 685-716.

Terzaghi, K. (1943). *Theoretical soil mechanics*. John Wiley & Sons, New York.

Thiesen, J. L. M., B. Klahr, T. A. Carniel, P. J. Blanco, and E. A. Fancello (2024). "A second-order multiscale model for finite-strain poromechanics based on the method of multiscale virtual power". In: *Journal of Elasticity* 156.3, pp. 917-954.

Thiesen, J. L. M., B. Klahr, T. A. Carniel, G. A. Holzapfel, P. J. Blanco, and E. A. Fancello (2025). "Second-order computational homogenization for bridging poromechanical scales under large deformations". In: *Computer Methods in Applied Mechanics and Engineering* 433, p. 117481.

Truesdell, C., W. Noll, C. Truesdell, and W. Noll (2004). *The non-linear field theories of mechanics*. Springer.

Tu, V., F. Larsson, K. Runesson, and R. Jänicke (2023). "Variationally consistent homogenization of electrochemical ion transport in a porous structural battery electrolyte". In: *European Journal of Mechanics-A/Solids* 98, p. 104901.

Wang, E., Q. Zhang, B. Shen, G. Zhang, X. Lu, Q. Wu, Y. Wang, E. Wang, Q. Zhang, B. Shen, et al. (2014). "Intel math kernel library". In: *High-Performance Computing on the Intel® Xeon Phi™: How to Fully Exploit MIC Architectures*, pp. 167-188.

Wriggers, P. (2008). *Nonlinear finite element methods*. Springer Science & Business Media.

Yang, Y., W. Wu, H. Zheng, S. Wang, and L. Yang (2024). "An efficient monolithic multiscale numerical manifold model for fully coupled nonlinear saturated porous media". In: *Computer Methods in Applied Mechanics and Engineering* 418, p. 116479.

Yuan, X., Y. Tomita, and T. Andou (2008). "A micromechanical approach of nonlocal modeling for media with periodic microstructures". In: *Mechanics Research Communications* 35.1, pp. 126-133.

Yvonnet, J., N. Auffray, and V. Monchiet (2020). "Computational second-order homogenization of materials with effective anisotropic strain-gradient behavior". In: *International Journal of Solids and Structures* 191, pp. 434-448.

Zeng, Z. and R. Grigg (2006). "A criterion for non-Darcy flow in porous media". In: *Transport in porous media* 63, pp. 57-69.