# Monolithic mechanical model incorporating solid contact for anisotropic porous materials

Li Zhang, Nicholas Anton Collins-Craft, Philipp Braun

► To cite this version:

Li Zhang, Nicholas Anton Collins-Craft, Philipp Braun. Monolithic mechanical model incorporating solid contact for anisotropic porous materials. 2026. hal-05535452

HAL Id: hal-05535452

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

Preprint submitted on 4 Mar 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Monolithic mechanical model incorporating solid contact for anisotropic porous materials

L. Zhang<sup>a,b</sup>, N.A. Collins-Craft<sup>a</sup>, and P. Braun<sup>a</sup>

<sup>a</sup>Laboratoire Navier, École Nationale des Ponts et Chaussées, Institut Polytechnique de Paris, Université Gustave Eiffel, CNRS, Champs-sur-Marne, France

<sup>b</sup>GEOMAR Helmholtz Centre for Ocean Research, Kiel, Germany

3rd March 2026

## Abstract

Understanding the mechanical behaviour of anisotropic porous materials as they deform is critical for applications such as geological carbon storage and subsurface energy systems. This study presents a comprehensive finite-element-based framework to simulate the evolution of porous solid media with frictionless solid contact, capturing anisotropic changes in pore geometry under complex multidirectional loading. The framework follows an innovative workflow with four steps: (1) a finite element formulation and simulator capable of resolving contact interactions with complex geometries; (2) a fast approximation method to the solution of the underlying contact problem; (3) a perturbation-based protocol for updating anisotropic stiffness matrices and Biot coefficients during progressive deformation; and (4) full integration of periodic boundary conditions into both the deformation analysis and perturbation formulations. The numerical framework allows homogenisation of systems that are not analytically tractable. Applied to single- and multi-pore systems, the framework demonstrates that progressive pore closure governs stiffness evolution, driving a transition from skeleton-controlled to contact-network-dominated behaviour. Young's modulus decreases along the loading direction while increasing orthogonally; shear loading reduces modulus in all directions and activates normal-shear coupling. Shear modulus initially rises due to oblique load transfer, whereas Poisson's ratio decreases under compression. However, following solid contact, these properties tend to approach those of the intact solid, but still remaining a distance under shear loading. Biot coefficients increase along the loading axis under uniaxial compression and rise in multiple directions under shear loading, but close to zero after solid contact under uniaxial compression process. This framework provides a robust and efficient tool for predicting nonlinear, anisotropic, and path-dependent poromechanical behaviour, linking pore-scale mechanisms to macroscopic material responses.

**Keywords:** Anisotropic porous structure; continuous deformation; solid contact; finite element method, perturbation-based stiffness updating.

## 1 Introduction

The mechanical behaviour of anisotropic porous materials plays an important role for the safe and efficient design of subsurface engineering applications, such as geological CO₂ sequestration, hydrogen storage, geohazard mitigation, and unconventional hydrocarbon extraction. To model and predict hydromechanical behaviour of rocks in the poroelasticity framework, macroscopic properties such as Young's modulus, shear modulus, and Biot coefficients are required. Microstructural evolution of porosity, solid-solid contact interactions, and structural anisotropy of the material collectively govern these macroscopic hydro-mechanical characteristics (Biot, 1973; Sahimi, 1993). Such evolution encompasses pore partitioning and fracture initiation (Lee et al., 2019), fracture closure and opening dynamics (Batzle et al., 1980) and pore clogging with resultant permeability alterations (Tang et al., 2020). Accurate prediction of these microscale properties, particularly their evolution under deformation that can induce changes in the contact network, is essential for reliable multiscale modelling of porous media.

Extensive research has been conducted on the continuous deformation of microscopic porous materials. Initially, Biot (1941) laid the theoretical foundation for poromechanics, enabling subsequent development of continuous-deformation models for porous materials. Coussy (2004) extended these principles to unsaturated media while incorporating both finite and infinitesimal strain formulations. Void geometry evolution ranging from collapse-induced softening to elongation-driven hardening, directly influences mechanical behaviour through porosity-structure coupling (Ponte Castañeda and Zaidman, 1994). The deformation of the microstructure and changes in the stress-transfer mechanisms
of mineral and rocks at the microscale has been shown to have macroscopic implications for the behaviour of reservoirs (Blenkinsop, 2002). Moreover, the integration of solid contact mechanics into models of porous media has gained momentum. Techniques ranging from coupled multiscale finite element and discrete element method (FEM×DEM) models with cohesion at the microscale (Wu et al., 2018), penalty-based and Lagrange-multiplier approaches in the particle finite element method (Carbonell et al., 2022) and discrete fracture network approaches (De Simone et al., 2023) have all been applied to modelling of porous media. Microscale modelling of soils and rocks is currently dominated by discrete element approaches, but these techniques are limited by their ability to accurately represent the micro-scale geometry, physically spurious stiffnesses and numerical instabilities in the case of classical DEM (Cundall and Strack, 1979), or perfect rigidity in the case of nonsmooth contact dynamics (Moreau, 1986; Jean and Moreau, 1987; Dubois et al., 2018) approaches. While approaches such as Level Set DEM (Kawamoto et al., 2016; Feldfogel et al., 2024) are able to account more accurately for the often-complex shape of grains, this comes at a high computational cost and does not remove the other difficulties associated with this technique. Current research combining the finite element method and contact mechanics techniques predominantly addresses larger scale questions such as soil-structure interactions or the performance of structural elements (Antunes et al., 2014), with limited application to the deformation and contact behaviour of evolving porous geomaterials.

Numerical methods for simulating solid contact problems include the boundary element method (Huang et al., 2021), DEM (Coetzee and Scheffler, 2023), meshless methods (Zhang et al., 2023), the material point method (Lei et al., 2022) and the finite element method (Pore et al., 2021). Among these, FEM offers the broadest applicability, as it is suitable for complex geometries, nonlinear materials, and continuous deformations, and has a large body of mathematical literature demonstrating its in approximating the solution of partial differential equations. Thus, it is the method adopted in this study. Common solution frameworks for solid contact model encompass Hertzian contact theory (Machado et al., 2012), penalty methods (Perić and Owen, 1992), Lagrange multiplier techniques (Papadopoulos and Solberg, 1998), augmented Lagrangian approaches (Burman et al., 2023), and nonsmooth mechanics formulations (Moreau, 1970, 1971; Moreau and Panagiotopoulos, 1988). Unlike penalty methods requiring artificial stiffness parameters, nonsmooth approaches directly compute contact forces with high precision and numerical stability by rigorously enforcing the non-penetration constraint through set-valued functions (Yastrebov, 2013) and are apt to be formulated in variational frameworks that demonstrate excellent energy-preservation properties (Capobianco and Eugster, 2016; Capobianco et al., 2017; Capobianco and Eugster, 2018; Acary et al., 2023; Capobianco et al., 2024). Further, nonsmooth mechanics formulations can often be transformed to favourable problem formulations arising from the optimisation community, of which the most favourable (and most frequent, as it arises in cases of contact between linear elastic bodies in geometrically linear frameworks) is the linear complementarity problem (LCP). However, while rigorous algorithms to obtain solutions of the LCP exist (such as Lemke's algorithm (Lemke and Howson, 1964; Lemke, 1965)), their computational complexity as the problem size increases can lead to high computational costs. In the context of multi-scale simulations, it thus of interest to have a solution that is sufficiently accurate for the purposes of obtaining effective material properties without imposing excessive computational demands.

In order to make calculations at the engineering scale, a constitutive model is required. Typically, these are phenomenological models that require strong assumptions on the part of the modeller, and often introduce a large set of material parameters that must be calibrated. By conducting a microscale analysis and applying a suitable homogenisation method, it is possible to derive material models that are closely tied to the material's physical properties and need minimal empirical calibration. Homogenisation techniques for determining effective material properties encompass analytical, semi-analytical and numerical approaches. Analytical solutions for idealised representative volume elements (RVEs) include stress fields and elastic properties of two-phase solid composites with isotropic spheres (Hill, 1965), ellipsoidal inclusions (Eshelby, 1957; Tsukrov and Kachanov, 2000), misfitting inclusions (Mori and Tanaka, 1973) and large-strain constitutive relations for elastic matrices with spherical inclusions (Hashin, 1985). For complex microstructures, numerical methods dominate, enabling homogenisation via perturbation-based stiffness matrices (Kouznetsova et al., 2001), and asymptotic analysis, mean-field approaches, transformation field analysis, and multiscale FEM (Kanouté et al., 2009; Geers et al., 2010). Analytical methods have been used to obtain drained elasticity tensors, Biot's coefficients and permeability as functions of solid volume fraction, the arrangement of spheres, the number of contacts as well as the intrinsic properties of the solid phase in simplified geometries (Bouhlel et al., 2010; Selvadurai and Suvorov, 2020). Analytical and numerical calculations at the microscale have also confirmed that poroelastic coefficients depend on pore architecture, global porosity, and solid-phase Poisson's ratio (Lydzba and Shao, 2000; Ahmed et al., 2019; Dehghani and Zilian, 2020). Microscale properties can be used to obtain macroscopic tangent moduli in the context of multiscale FEM (FEM²) calculations (Saeb et al., 2016), allowing the passage between highly-refined models at the microscale to effective models at the macroscale. The evolution of these tangent moduli due to a variety of microscale effects is critical for maintaining an accurate macro-scale simulation as the system deforms, and a variety of methods to obtain them is avaiable, from the perturbation-based stiffness matrix of multiscale nonlinear composites (Temizer and Wriggers, 2008), stress-tangent moduli derived from discrete boundary forces (Miehe, 2003), geometric softening/stiffening in linear comparison composites (Lopez-Pamies and Ponte Castañeda, 2004; Lopez-Pamies, 2006; Lopez-Pamies and Ponte Castañeda, 2007a,b), microstructure evolution in neo-Hookean matrices with cylindrical voids (Moraleda et al., 2007) and other non-linear elastic matrices (Miehe et al., 2002; van Dijk, 2016).
Conventional iterative stress-update approaches, however, require microscopic model evaluations at every macroscopic step, which is a computationally prohibitive practice, particularly for materials with strong nonlinearities or complex microstructures (Ghosh et al., 2001; Miehe and Koch, 2002). Perturbation-based stiffness updating addresses this limitation by reducing iterative cycles. It computes the RVE constitutive stiffness coefficients and mechanical parameters on demand, triggering updates only when local responses deviate from predefined thresholds (Charpin and Ehrlacher, 2014; Zhou et al., 2016). This approach not only facilitates efficient parameter transfer of stresses and fluxes, but also enables intrinsic mechanical property prediction of the microstructure itself (Miehe and Koch, 2002). While there has been some limited consideration of solid contact effects in microscale models (e.g. Rohan and Heczko (2020) developed a homogenised model accounting for contact between rigid and compliant phases for periodic microstructures with rigid inclusions), a critical gap persists in existing multiscale perturbation methods, which typically neglect configuration updates and solid contact effects (Nemat-Nasser and Hori, 1998; Charpin and Ehrlacher, 2014; Tan and Konietzky, 2014). Such omissions undermine their applicability to problems involving continuous deformations, where evolving geometry and interfacial interactions dominate mechanical responses. This study advances perturbation-based frameworks by integrating real-time configuration updating and contact mechanics, thereby extending their validity to systems in which the material fabric substantially evolves, as is common in geomechanical and metamaterial systems.

To overcome the limitations of existing studies in capturing porous solid behaviour with complex pore-scale interactions, this work develops a unified micromechanical framework that integrates finite deformation with solid-solid contact mechanics. The framework advances current methodologies through four key innovations: (1) a finite element formulation capable of accurately resolving solid contact interactions with complex geometries, implemented within a dedicated computational simulator; (2) a fast approximation method to the solution of the LCP; (3) a perturbation-based stiffness updating scheme that automatically evaluates anisotropic stiffness matrices and Biot coefficients via small-strain perturbations during configuration updating; and (4) seamless integration of stress, displacement, and periodic boundary conditions into both the deformation analysis and stiffness matrix computation. The framework is then applied to single-pore and multi-pore configurations to analyse the directional stiffness evolution and post-contact behaviour during multidirectional loading. The simulator's capability to dynamically update poromechanical properties while enforcing contact constraints enables quantitative assessment of nonlinear and anisotropy-dependent poromechanical performance, which is of interest for predicting geomaterials properties in underground storage systems and optimising anisotropic metamaterials.

## 2 Material model

In this section, the underlying system of mathematical governing equations for solid deformation and solid frictionless contact is introduced. Throughout this work, we consider a solid material, that surrounds an empty pore space, where the material properties of the solid are homogeneous and isotropic. We first consider the solid material as a body and introduce the appropriate state variables and powers, before defining the material model through the introduction of an energy potential for the bulk and for the surface. We then proceed to the spatial discretisation of the system using the finite element method, resulting in a semi-discretised system. Then, we treat the system of the solid material and the included pores as a single system to model a representative volume element (RVE) of a porous material by introducing period boundary conditions and a perturbation method to compute the macroscopic properties of the system.

<figure><img src="image_2.png" /><figcaption>Figure 1: (a) Microscopic porous media model with arbitrary pore structure (white for pore, grey for solid, incorporating (b) localised solid contact interactions (adapted from Batistić et al. (2022)) and (c) parallel solid contact configuration.</figcaption></figure>
## 2.1 State variables and powers

We consider a body $\Omega \in \mathbb{R}^3$ with boundaries denoted by $\Gamma$, where the initial position is defined by $\mathbf{X}$ and the current position $\mathbf{x}$. We consider that the geometry and properties of the system are invariant with respect to the $x_2$ variable (that is to say that we treat plane strain or plane stress problems) and so the dimensionality of the system can be reduced with $\mathbf{x} = (x_1, x_3)$ (and similarly for $\mathbf{X}$). The displacement is then given by $\mathbf{u}(\mathbf{x}) = \mathbf{x} - \mathbf{X}$ and the velocity by $\mathbf{v} = \dot{\mathbf{u}}$, where the overdot indicates the derivative with respect to time. On the boundaries of the system, we may define a local orthonormal reference frame denoted $(\mathbf{x}, \mathbf{n}, \mathbf{t})$ where $\mathbf{n} \in \mathbb{R}^2$ is the normal unit vector pointing outwards from the surface, and the vector $\mathbf{t} \in \mathbb{R}^2$ is the tangent direction along the surface, positive in the direction that locally respects the right-hand rule. We may consider the relative displacement between two different points, both lying on a surface by $u_N(\mathbf{x}^1, \mathbf{x}^2) = (\mathbf{u}(\mathbf{x}^1) - \mathbf{u}(\mathbf{x}^2)) \cdot \mathbf{n} \in \mathbb{R}$ where $\mathbf{n}$ is taken at $\mathbf{x}^1$. The relative normal velocity is then given by $v_N = \dot{u}_N$.

For a point lying on the surface of the body, the power of contact is given by

$$ \mathcal{P}_{\text{con}} = \int_{\Gamma} v_n r_N \, dx, \qquad (1) $$

where $r_N$ is the intensity of the normal reaction force related to the Cauchy stress $\sigma(\mathbf{x})$ at the surface of the body by $r_N = -\boldsymbol{\sigma} \cdot \mathbf{n} \cdot \mathbf{n}$. Then, we have the classical powers of the external, internal and inertial forces given by

$$ \mathcal{P}_{\text{ext}} = \int_{\Omega} \mathbf{v} \cdot \mathbf{f} \, dx - \int_{\Gamma_N} \mathbf{v} \cdot \boldsymbol{\tau} \, dx, \qquad (2) $$

$$ \mathcal{P}_{\text{int}} = - \int_{\Omega} \boldsymbol{\sigma} : \dot{\boldsymbol{\varepsilon}} \, dx + \int_{\Gamma} v_n r_N \, dx, \qquad (3) $$

$$ \mathcal{P}_{\text{acc}} = \int_{\Omega} \rho \mathbf{v} \cdot \dot{\mathbf{v}} \, dx, \qquad (4) $$

where $\mathbf{f}$ is the body force in $\Omega$, $\boldsymbol{\tau}$ is the imposed surface traction on $\Gamma_N$ (the region of the boundary that has the Neumann boundary condition applied), $\boldsymbol{\varepsilon} = \frac{1}{2}(\nabla \mathbf{u} + (\nabla \mathbf{u})^{\top})$ is the Cauchy strain, $\rho$ is the mass density and $\dot{\mathbf{v}}$ is the acceleration. The principle of virtual power ensures that for any virtual velocities $\bar{\mathbf{v}}, \dot{\bar{\mathbf{v}}}$ and $\bar{v}_N$ the following holds:

$$ \bar{\mathcal{P}}_{\text{acc}} = \bar{\mathcal{P}}_{\text{ext}} + \bar{\mathcal{P}}_{\text{int}}, \\ \int_{\Omega} \rho \bar{\mathbf{v}} \cdot \dot{\mathbf{v}} \, dx = \int_{\Omega} \bar{\mathbf{v}} \cdot \mathbf{f} \, dx - \int_{\Gamma_N} \bar{\mathbf{v}} \cdot \boldsymbol{\tau} \, dx - \int_{\Omega} \boldsymbol{\sigma} : \dot{\boldsymbol{\varepsilon}} \, dx + \int_{\Gamma} \bar{v}_n r_N \, dx. \qquad (5) $$

From here on we assume both that the system may be treated quasi-statically (that is, it evolves in time, but sufficiently slowly that the power of the inertial forces is zero), and that the relevant fields are sufficiently smooth away from contact boundaries. Localisation gives the equilibrium and boundary conditions of the system as

$$ \left\{ \begin{array}{ll} \nabla \cdot \boldsymbol{\sigma} + \mathbf{f} = 0 & \text{in } \Omega, \\ \boldsymbol{\tau} = \boldsymbol{\sigma} \cdot \mathbf{n} & \text{on } \Gamma_N, \\ \mathbf{u} = \mathbf{u}_D & \text{on } \Gamma_D, \\ r_N = -\boldsymbol{\sigma} \cdot \mathbf{n} \cdot \mathbf{n} & \text{on } \Gamma \setminus (\Gamma_N \cup \Gamma_D), \end{array} \right. \qquad (6) $$

where $\Gamma_D$ is the region of the boundary that has the Dirichlet boundary condition applied and $u_D$ is the prescribed displacement.

## 2.2 Energy potential

The free energy of the system gives us the laws that describes the reversible behaviour of the system. The free energy of the system $\Psi$ is given by the sum of the free energy of the bulk material $\Psi_e$ with the free energy of the surface $\Psi_s$:

$$ \Psi = \int_{\Omega} \Psi_{e}(\boldsymbol{\varepsilon}) \, dx + \int_{\Gamma} \Psi_{s}(u_{N}) \, dx. \qquad (7) $$

In this work we consider the bulk material to be purely linear elastic (and so the strain tensor is equivalent to the elastic strain tensor) and we may thus obtain the constitutive expression for the stresses with a classical specification of the bulk energy potential:

$$ \Psi_e(\boldsymbol{\varepsilon}) = \frac{1}{2} : \boldsymbol{\varepsilon} : \mathbf{D} : \boldsymbol{\varepsilon}, \qquad (8) $$

$$ \boldsymbol{\sigma}(\boldsymbol{\varepsilon}) = \frac{\partial \Psi_e(\boldsymbol{\varepsilon})}{\partial \boldsymbol{\varepsilon}} = \mathbf{D} : \boldsymbol{\varepsilon}, \qquad (9) $$
where $D$ is a fourth order elastic stiffness tensor. For convenience, we may transform this tensor for full 3D isotropic elasticity, using Voigt notation, into a two-dimensional tensor given by

$$ D_s = \frac{E_s}{(1 + \nu_s)(1 - 2\nu_s)} \begin{bmatrix} 1 - \nu_s & \nu_s & \nu_s & 0 & 0 & 0 \\ \nu_s & 1 - \nu_s & \nu_s & 0 & 0 & 0 \\ \nu_s & \nu_s & 1 - \nu_s & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1-2\nu_s}{2} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1-2\nu_s}{2} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1-2\nu_s}{2} \end{bmatrix}, \quad (10) $$

where $E_s$ and $\nu_s$ respectively denote the Young's modulus and Poisson's ratio of the solid material. We can likewise obtain the intensity of the reaction force by an appropriate specification of the surface energy:

$$ \Psi_s(u_N) = \mathcal{I}_{\mathbb{R}_+}(u_N), \quad (11) $$

$$ -r_N \in \partial_{u_N} \Psi_s(u_N), \quad (12) $$

where $\partial_{u_N}$ indicates the subdifferential with respect to $u_N$ and $\mathcal{I}_C$ is the indicator function of a convex set $C$. As the subdifferential of the indicator set of the positive real numbers is the normal cone to the positive reals, we may also write the reaction force intensity as

$$ -r_N \in \mathbb{N}_{\mathbb{R}_+}(u_N), \quad (13) $$

which can in turn be further simplified as a complementarity relation:

$$ 0 \le r_N \perp u_N \ge 0. \quad (14) $$

This relation expresses unilateral contact with impenetrability and no adhesion.

## 2.3 Spatially discretised systems

We consider now a finite-dimensional model of the system, where we discretise the model using the classical linear finite element method. We use Galerkin's method, so that the shape and weight functions are identical, and proceed in the standard fashion. While the method detailed in this paper can easily be generalised to higher order elements, here we further specialise in assuming that the discretisation is done using isoparametric linear triangle (T3) or quadrilateral (Q4) elements. These elements sufficiently balance accuracy and the size of the discrete system for our purposes.

Now, we denote our finite set of displacements $u \in \mathbb{R}^n$, the external set of forces $F \in \mathbb{R}^n$, and the structural stiffness matrix assembled in the standard way as $K \in \mathbb{R}^{m \times n}$. This matrix will be symmetric and positive semi-definite. We consider also a finite set of contact points that are indexed by $\alpha \in [1, m]$, and we assume that the local normal displacements at the contact (given by $u_N = \{\text{col}(u_N^\alpha), \alpha \in [1, m]\}$) can be related to the global displacement with a linear relation given by

$$ u_N = -B_c u + c_N, \quad (15) $$

where $B_c \in \mathbb{R}^{m \times n}$ is a matrix that selects and interpolates the applicable degrees of freedom and $c_N \in \mathbb{R}_+^m$ is an initial gap function. We can also construct a diagonal matrix $S \in \mathbb{R}^{m \times m}$ that contains the tributary area of each contact force intensity. For the sake of notational simplicity, from here on we refer to the product of this tributary area with the reaction force intensity as the contact force $f_c^\alpha = S^\alpha r_N^\alpha$, which can be converted to a column vector in the same way as $u_N$. As the components of $S$ are all positive and finite, and there is no interaction between the different terms, the complementarity relationship between $r_N$ and $u_N$ can be transferred to be between $f_c$ and $u_N$. Finally, this gives the spatially-discretised system

$$ \left\{ \begin{array}{l} Ku = F - B_c^\top f_c, \\ u_N = -B_c u + c_N, \\ 0 \le u_N \perp f_c \ge 0, \end{array} \right. \quad (16) $$

where the inequalities and complementarity are to be understood component-wise.

## 2.4 Periodic boundary conditions

In this work, the structure that we have discretised with the finite element is in fact the *microstructure* of the RVE of a porous material, in order to obtain its average macroscopic poromechanical properties. The Hill-Mandel principle of macro-homogeneity (Hill, 1963; Mandel, 1971) establishes a fundamental connection between macroscopic and microscopic deformation and stress. It states that the variation of work at a macroscopic point must equal the volume average of the microscopic work variation within the RVE. This principle also imposes constraints on the boundary conditions used in RVE simulations: either prescribed displacements, prescribed tractions, or periodic boundary conditions. Among these, periodic boundary conditions generally yield the most accurate homogenised properties
(Terada et al., 2000; Kanit et al., 2003; Larsson et al., 2011; Nguyen et al., 2012). To impose periodic boundaries, displacement differences between corresponding boundaries are constrained. For the sake of clarity, the body Ω with boundary Γ that we have previously considered is merely the solid part of the microstructure. We can also treat the entire system (inclusive of the voids) as a body with a boundary for the purposes of homogenisation. To this end, we use the symbol V to refer to the RVE in its entirety, and its boundary as ∂V. While Ω is completely contained in V, only the Dirichlet and Neumann parts of Γ are necessarily contained in ∂V, and in general the “internal” boundaries on which the contact problem is solved will not belong to ∂V.

The RVE boundary ∂V is decomposed into two parts: a positive part ∂V + and a negative part ∂V − , satisfying ∂V + ∪ ∂V − = ∂V, ∂V + ∩ ∂V − =  , with the outward normal n + = − n − at conjugate points x + ∈ ∂V + and x − ∈ ∂V − . This enforces periodicity of displacement fluctuations and antiperiodicity of tractions, as

$$
\begin{equation}
\begin{cases}
\tilde{u}(\mathbf{x}^{+}) - \tilde{u}(\mathbf{x}^{-}) = \varepsilon\mathbf{x} & \forall \mathbf{x}^{+} \in \partial V^{+} \text{ and matching } \mathbf{x}^{-} \in \partial V^{-}, \\
\sigma(\mathbf{x}^{+}) = -\sigma(\mathbf{x}^{-}) & \forall \mathbf{x}^{+} \in \partial V^{+} \text{ and matching } \mathbf{x}^{-} \in \partial V^{-}.
\end{cases}
\tag{17}
\end{equation}
$$

2.5 Perturbation method for macroscopic stiffness matrix

Here we use plane strain conditions to approximate a 3D material that is invariant and infinitely extended in the 2-direction. We want to determine the macroscopic elastic properties of the RVE, so in two dimensions and under plane strain conditions in direction 2 the following constitutive relation holds:

$$
\begin{bmatrix} \sigma_1 \\ \sigma_2 \\ \sigma_3 \\ \sigma_6 \end{bmatrix} = \begin{bmatrix} D_{11} & D_{12} & D_{13} & D_{16} \\ D_{21} & D_{22} & D_{23} & D_{26} \\ D_{31} & D_{32} & D_{33} & D_{36} \\ D_{61} & D_{62} & D_{63} & D_{66} \end{bmatrix} \begin{bmatrix} \epsilon_1 \\ 0 \\ \epsilon_3 \\ \epsilon_6 \end{bmatrix}, \quad (18)
$$

where σ$_{i}$ and ε$_{i}$ denote the stress and strain tensor components, respectively.

If the behaviour of the RVE is isotropic or transversely isotropic with respect to direction 2, the parameters D$_{61}$ , D$_{62}$ , D$_{63}$ and D$_{16}$ , D$_{26}$ and D$_{36}$ remain zero. Otherwise, these parameters acquire nonzero values, especially in large shear displacement. To determine the components of D , we apply the perturbation method (Elishakoff et al., 1997; Charpin and Ehrlacher, 2014), in which small perturbations are applied separately in horizontal, vertical, and shear displacements, and the corresponding stress variations δσ are computed for each perturbation direction. The average stress change is then obtained as

$$
\delta\boldsymbol{\sigma} = \frac{\int_V \delta\boldsymbol{\sigma}_{\text{Gauss}} \, dV}{\int_V dV}, \quad (19)
$$

where δσ represents the stress change of Gauss points in the volume V due to small displacement perturbations and σ$_{Gauss}$ represents the stress state measured at a Gauss point.

In two dimensional plane strain conditions, the elastic modulus and Poisson’s ratio in the x$_{2}$ -direction are func- tions of the mechanical properties of the solid matrix and porosity n (Nguyen et al., 2012; Charpin and Ehrlacher, 2014), i.e.

$$
E_2 = (1 - n)E_s, \tag{20}
$$

$$
\nu_{21} = \nu_{23} = \nu_s. \tag{21}
$$

Accordingly, the Young’s modulus in the x$_{1}$ - and x$_{3}$ -directions, E$_{1}$ and E$_{3}$ , shear modulus G$_{13}$ and Poisson’s ratio ν$_{13}$ can be expressed using the calculated stiffness parameters D$_{ij}$ (Charpin and Ehrlacher, 2014), i.e.

$$
E_1 = \frac{\left(\frac{D_{11}}{D_{31}}\right)^2 - \left(\frac{D_{13}}{D_{33}}\right)^2}{\frac{1}{D_{31}}\left(\frac{D_{11}}{D_{31}} + \frac{\nu_s^2}{E_2}\left[\left(\frac{D_{11}}{D_{31}}\right)^2 - 1\right]\right) + \frac{1}{D_{33}}\left(1 + \frac{\nu_s^2}{E_2}\left[1 - \left(\frac{D_{13}}{D_{33}}\right)^2\right]\right)}, \quad (22)
$$

$$
E_3 = \frac{\left(\frac{D_{33}}{D_{13}}\right)^2 - \left(\frac{D_{31}}{D_{11}}\right)^2}{\frac{1}{D_{11}}\left(1 + \frac{\nu_s^2}{E_2}\left[1 - \left(\frac{D_{31}}{D_{11}}\right)^2\right]\right) + \frac{1}{D_{31}}\left(\frac{D_{33}}{D_{13}} + \frac{\nu_s^2}{E_2}\left[\left(\frac{D_{33}}{D_{13}}\right)^2 - 1\right]\right)}, \quad (23)
$$

and

$$
\nu_{13} = \frac{D_{13}}{D_{33}} - \frac{E_1 \nu_s^2}{E_2} \left( \frac{D_{13}}{D_{33}} + 1 \right), \quad (24)
$$

$$
G_{13} = D_{66}. \tag{25}
$$
By incorporating these into the stiffness matrix formulation, we obtain

$$D_{13} = D_{31}, \qquad (26)$$

$$D_{23} = D_{32}, \qquad (27)$$

$$\frac{\nu_{13}}{E_1} = \frac{\nu_{31}}{E_3}, \qquad (28)$$

$$\frac{\nu_{23}}{E_2} = \frac{\nu_{32}}{E_3}. \qquad (29)$$

Thus, the last unknown parameter in the stiffness matrix, $D_{22}$, is determined (Tan and Konietzky, 2014):

$$D_{22} = \frac{E_2^2 E_3 (E_1 - \nu_{13}^2 E_3)}{E_1 E_2 E_3 - E_1^2 E_2 \nu_s^2 - \nu_{13}^2 E_2 E_3^2 - E_1 E_3^2 \nu_s^2 - 2 E_1 E_3^2 \nu_{13} \nu_s^2}. \qquad (30)$$

In addition to the classical elastic properties of the homogenised material, we are also interested in obtaining an estimate of the Biot coefficients in each direction of the material, which are given by (Tan and Konietzky, 2014)

$$b_1 = 1 - \frac{D_{11} + D_{12} + D_{13}}{3K_s}, \qquad (31)$$

$$b_2 = 1 - \frac{D_{12} + D_{22} + D_{23}}{3K_s}, \qquad (32)$$

$$b_3 = 1 - \frac{D_{13} + D_{23} + D_{33}}{3K_s}, \qquad (33)$$

where $K_s$ is the bulk modulus of the solid phase, i.e. $K_s = \frac{E_s}{3(1-2\nu_s)}$.

## 3 Numerical solution

Having obtained an analytical description of the spatially-discretised but time-continuous system, we now turn to the time-discretised system and its numerical resolution. We first describe the pre-processing of the system that allows the spatially discrete system to be imported and correctly interpreted by the numerical simulator. We then describe the time-stepping scheme and the contact detection and classification algorithm, leading to the formulation of the fully discretised system as a linear complementarity problem. We prove that this problem has solutions, and provide the condition under which the solution is unique. The central challenge in solid contact analysis lies in the presence of two unknowns ($u_N$ and $f_c$) where at least one variable must vanish at each contact point. To overcome these limitations, we then describe the algorithm developed to provide a fast approximation of the solution of this system. Finally, we describe the numerical resolution of the perturbation problem. The entire numerical solution process is integrated into an in-house simulator, with the iterative workflow illustrated in Figure 2. The simulation codes, written in Fortran 90, were run on Windows 10/11 and compiled with the Intel Fortran Compiler 2024.2 (IFORT, oneAPI HPC Toolkit).
<figure><img src="image_3.png" /><figcaption>Figure 2: FEM workflow for micromechanical deformation analysis with consideration of solid contact resolution and automated stiffness matrix updating.</figcaption></figure>

## 3.1 Pre-processing

In order to efficiently carry out the simulations, we first pre-process the system. The mesh is generated either using the commercial Gambit software or with the Python toolbox Nanomesh (Smeets et al., 2022) and then converted to Gambit format with an in-house developed script. The generated Gambit neutral files are directly read in by our code. However, nothing about our method requires use of this particular meshing software, so any suitable software that produces conforming meshes can easily be substituted.

Then, we analyse the system to identify which elements lie on the pore surface and thus are able to be in contact with each other. We loop through each edge in the mesh, and if that edge is present only in one element, we identify it as a possible contact surface. Then, for all of the contact surfaces, we construct an adjacency list that records which other surface nodes it is directly connected to by a boundary edge. Finally, all nodes that are directly connected through boundary edges are assigned to a contact group, and members of each group can only (potentially) be in contact with other members of the same group. This ensures that no computational resources are expended checking contact conditions on pairs of contact surfaces that belong to different pores.

Finally, the structural stiffness matrix **K** is generated by assembling the element stiffness matrices following the standard linear finite element methodology. The structural stiffness matrix is then modified in the standard way to impose the Dirichlet boundary conditions (making the matrix positive definite), and we may obtain the final linear equation system of LCP. In practice we do not actually invert the matrix and instead we obtain the LU decomposition, and each time $K^{-1}$ (or any other matrix inverse) appears below it should be understood to actually be the LU decomposition acting on the appropriate terms.
## 3.2 Initial iteration

We consider a time-discretised system, where the time steps are indexed by $k \in \mathbb{N}$ and the time at step $k$ is given by $t_k$. We denote the displacement and external force vector at time $t_k$ as $\mathbf{u}_k$ and $\mathbf{F}_k$ respectively. We apply a backwards-Euler time-discretisation to the system, as is common in quasi-statics, that is

$$ \int_{t_k}^{t_{k+1}} x(t) dt \approx \Delta t x(t_{k+1}). \qquad (34) $$

We begin the process of updating the system by updating the external force vector $\mathbf{F}_{k+1}$, using the user-defined load increments. Any force contributions of elements lying on Neumann boundaries are obtained, and then all of the element contributions are assembled into the global force vector. The force vector is then modified in the standard way to impose the Dirichlet boundary conditions at $t_{k+1}$. We then calculate the trial displacement assuming no contacts, which gives the first internal iterate of the system:

$$ \mathbf{u}_{k+1} = \mathbf{K}^{-1} \mathbf{F}_{k+1}, \qquad (35) $$

From this point, we may pass to the detection and classification of contacts.

## 3.3 Contact detection and classification

### 3.3.1 Contact detection

Once the first internal iterate has been obtained, we use this virtual geometry as the system configuration on which the contact detection operates. We implement Algorithm 1 to check for the presence of contacts.

<caption>**Algorithm 1** Contact detection algorithm with displacement scaling and boundary-surface verification</caption>

Require: contact groups, $\mathbf{x}^0$, $\mathbf{u}^1$, total scaling steps $n$
for group $\in$ contact groups **do**
    ▷ Check each possible contact group separately
    for element1 $\in$ group **do**
        for element2 $\in$ group\elements<element1 **do**
            contactelement1element2 $\leftarrow$ FALSE
            ▷ Assume no contact
            for $i = n$ to 1 **do**
                $\bar{\mathbf{u}}_i^1$ $\leftarrow$ $\frac{i}{n}\mathbf{u}^1$
                for node1 $\in$ contact surface of element1 **do**
                    for node2 $\in$ contact surface of element2 **do**
                        $\mathbf{u}_{N,\text{node1node2}}$ $\leftarrow$ $(\mathbf{x}_{\text{node1}}^0 + \bar{\mathbf{u}}_{i,\text{node1}}^1) - (\mathbf{x}_{\text{node2}}^0 + \bar{\mathbf{u}}_{i,\text{node2}}^1)$
                        if $\mathbf{u}_{N} \le 0$ **then**
                            in12 $\leftarrow$ any part of element1 boundary lies inside element2
                            in21 $\leftarrow$ any part of element2 boundary lies inside element1
                            if in12 = TRUE and in21 = TRUE **then**
                                contactelement1element2 $\leftarrow$ TRUE
                                ▷ Contact detected between element pair
                                break
                            end if
                        end if
                    end for
                    if contactelement1element2 = TRUE **then**
                        break
                    end if
                end for
                if contactelement1element2 = TRUE **then**
                    break
                end if
            end for
            end for
            end for
            end for
        end for
    end for
    return List of element pairs for which contactelement1element2 is TRUE

This algorithm loops through the elements in each contact group, then loops again through the elements with a higher element number in the same contact group (ensuring that a given pair of elements is checked only once), then checks each pair of nodes in these elements (in12 and in21) to determine if the "contactelement1element2" (contact between element1 and element2) is satisfied, when the normal gap function (15) is negative. However, during
the detection of contacts, only potential boundary elements are assessed, excluding internal elements. Excessive displacements may cause over-penetration of contact boundaries, as shown in Figure 3(a). In this example, Element 1 undergoes large deformation in the initial step and penetrates through Boundary Element 2 to contact Internal Element 3, which is not evaluated by the detection routine. In reality, Element 1 should first establish contact with Element 2. Consequently, some potential contact pairs may be missed. To mitigate against this, geometric overlap is re-evaluated by progressively scaling displacements from the updated virtual configuration back to the initial state, given by

$$ \bar{u}_i^1 = \frac{i}{n} u_{k+1}^1, \qquad (36) $$

where $i \in \{1, 2, ..., n\}$ and n indicates the total scaling steps. This operation is applied to all element pairs for which contact is not detected in the initial virtual geometry. Then, for each i considered, contacts are checked again. This operation is repeated until all element pairs in the same boundary group have been checked, and then repeated for each boundary group until the contacts are fully characterised. Contact groups are identified if overlap occurs at any scaled configuration.

<figure><img src="image_4.png" /><figcaption>Figure 3: (a) Displacement modification for over contact and (b) overlap boundary modification for cross contact in solid contact group recognition.</figcaption></figure>

Additionally, even if two boundary elements are found to overlap, direct contact is only confirmed if their surfaces actually interact, as shown in Figure 3(b). In this case, Element 1 overlaps with Element 3 by passing through the intermediate Boundary Element 2. However, the true contact should occur between Element 1 and Element 2, not Element 3. Therefore, contact is considered established when parts of the boundary surfaces of both interacting elements are simultaneously located within the contact region, indicating a partial overlap. If the apparent overlap involves only one element's boundary, the pair is treated as cross contact and non-contacting. If no contacts are detected for all values of i, there are no contacts in that time step and $u^1$ is accepted as the solution. If contacts are detected, we then proceed to §3.3.2. In such cases, the time step $\Delta t$ is adaptively reduced to enhance the resolution of incremental strain, stress, and displacement changes, allowing accurate identification of the transition from non-contact to contact states while minimising numerical errors. Therefore, different initial time steps are assigned to the pre-contact and post-contact phases. Before contact occurs, a relatively large initial time step is used to accelerate convergence. Once the first contact is detected, the time step is automatically reduced according to user-defined parameters (by default, it is reduced to one-tenth of the initial value). This adaptive reduction ensures smooth transition between loading stages and prevents over-penetration. After the first contact is established, the maximum allowable time step is also reset, providing tighter control on subsequent deformation increments.

### 3.3.2 Contact classification

We apply a surface-to-surface contact scheme where we assume that the contacting surfaces of two elements are parallel at the moment of contact, but are not necessarily aligned node-to-node. We assume a constant distribution of pressure along the contacting surfaces and apply a set of characterisation rules to appropriately allocate the contact forces to the nodes in the problem. For a pair of contacting elements identified by Algorithm 1, which possess contacting surfaces of length $\ell$ and $\ell'$ with $\ell' \le \ell$, the length of the overlapping part of the surfaces is given by $b$ and the distance from the node of the contacting surface with length $\ell$ to the edge of the overlapping region is given by $a$, as shown in the subfigures in Figure 4. The surface normals and the initial gap functions are evaluated from the geometry of the previous time step. The element associated with the contacting surface of length $\ell$ is assigned the outward-pointing normal vector $n$, while its counterpart with surface length $\ell'$ is assigned the opposite normal vector $-n$. We construct a contact interpolation and selection matrix $B_c \in \mathbb{R}^{m \times n}$ where $m$ is the number of element contact pairs identified by
Algorithm 1. We consider one contact problem per element contact pair. The contact type is determined based on the number of boundary nodes of one element that lie inside the other element, and the contact selection and interpolation matrix is constructed accordingly, as summarised in Figure 4.

<figure><img src="image_5.png" /><figcaption>Figure 4: The contact selection and interpolation matrix construction rules. $B_{c,1}$ and $B_{c,2}$ should be understood to refer to the entries that correspond the appropriate row for the given contact pair, and appropriate columns for the node under consideration.</figcaption></figure>

The $B_{c,1}$ and $B_{c,2}$ entries of the $\ell$ element and $\ell'$ element are then assembled into the appropriate row of the contact selection and interpolation matrix $B_c$. The remaining entries in that row are set to zero. The normal contact distance vector $c_{N,k}$ at time $t_k$ is then obtained by taking the average normal distance of the corresponding pairs of nodes belonging to the contacting surfaces of each contact pair.

## 3.4 Contact problem formulation

Having completed the contact detection and classification step, we have all of the necessary information to construct the contact problem at time $t_{k+1}$. The normal contact distance is given by

$$ \mathbf{u}_{N,k+1} = -B_c \mathbf{u}_{k+1} + c_{N,k}. \qquad (37) $$

The global discrete force equilibrium equation is given by

$$ K \mathbf{u}_{k+1} = F_{k+1} - B_c^\top f_{c,k+1}, \qquad (38) $$

where $f_{c,k+1}$ is the contact force between a pair of contacting elements, and the contact distance and contact force respect the discrete time complementarity condition

$$ 0 \le \mathbf{u}_{N,k+1} \perp f_{c,k+1} \ge 0, \qquad (39) $$

where the complementarity should be understood to act pair-wise (that is the $i$-th entry of $\mathbf{u}_{N,k+1}$ is in complementarity with the $i$-th entry of $f_{c,k+1}$). Straightforward re-arrangement and combination of (37), (38) and (39) gives the system

$$ \left\{ \begin{array}{l} \mathbf{u}_{N,k+1} = B_c K^{-1} B_c^\top f_{c,k+1} - B_c K^{-1} F_{k+1} + c_{N,k}, \\ 0 \le \mathbf{u}_{N,k+1} \perp f_{c,k+1} \ge 0. \end{array} \right. \qquad (40) $$

This is a linear complementarity problem in discrete time. For reasons of clarity and compactness we denote the quantity $B_c K^{-1} B_c^\top$ as $W$, and the quantity $-B_c K^{-1} F_{k+1} + c_{N,k}$ as $q$. The $W$ is often referred to in contact mechanics as the Delassus matrix (Brogliato, 2016). We refer to (40) compactly as LCP($W, q$).

### 3.4.1 Well-posedness of the linear complementarity problem

**Lemma 1.** The LCP($W, q$) possesses at least one solution.

**Proof** First, we must show that $W$ is at least positive semi-definite. We have that the structural stiffness matrix $K$ is symmetric by construction, and the modifications made to impose the Dirichlet boundary conditions change it from positive semi-definite to positive definite, without modifying the symmetry, and this property necessarily holds
for its inverse. Any matrix product of the form $ABA^T$ is at least positive semi-definite if $B$ is positive definite. As this is exactly the structure of $W$ with respect to $B_c$ and $K^{-1}$, we have that $W$ is at least positive semi-definite.

From Theorem 3.1.2 in Cottle et al. (2009), if the LCP($W$, $q$) with a positive semi-definite $W$ is feasible, it is also solvable. The feasibility conditions are that $Wf_{c,k+1} + q \ge 0$ and $u_{N,k+1} \ge 0$. We may consider a convex quadratic program given by

$$ \begin{array}{ll} \text{minimise} & \frac{1}{2} f_{c,k+1}^\top W f_{c,k+1} + f_{c,k+1}^\top q, \\ \text{subject to} & f_{c,k+1} \ge 0. \end{array} \qquad (41) $$

A solution of (41) exists and satisfies the optimality conditions $u_{N,k+1}^* = Wf_{c,k+1}^*+q \ge 0$, $f_{c,k+1}^* \ge 0$, $u_{N,k+1}^* \ge 0$, $u_{N,k+1}^* \ge 0$, that is to say that a feasible point exists, and therefore so does a solution.

**Lemma 2.** The LCP($W$, $q$) has a unique solution for $c_{N,k} > 0$.

**Proof** We denote the set of solutions of LCP($W$, $q$) by SOL($W$, $q$), which by Theorem 3.1.7 of Cottle et al. (2009) is a polyhedral set when $W$ is positive semi-definite. This solution set is explicitly given by

$$ \mathrm{SOL}(W, q) = \{ f_{c,k+1} \mid f_{c,k+1} \ge 0, Wf_{c,k+1} + q \ge 0, q^\top (f_{c,k+1} - \bar{f}_{c,k+1}) = 0, W(f_{c,k+1} - \bar{f}_{c,k+1}) = 0 \}, \quad (42) $$

where $\bar{f}_{c,k+1}$ is an arbitrary solution. The final condition $W(f_{c,k+1} - \bar{f}_{c,k+1}) = 0$ gives us $u_{N,k+1} - \bar{u}_{N,k+1} = W(f_{c,k+1} - \bar{f}_{c,k+1}) = 0$ from which we can immediately conclude that $u_{N,k+1}$ is unique. We also conclude that $f_{c,k+1} - \bar{f}_{c,k+1}$ is in the kernel of $W$ and thus also of $B_c^\top$, so we obtain that $B_c^\top (f_{c,k+1} - \bar{f}_{c,k+1}) = 0$. Then, from the second condition of (42), we have

$$ (-B_c K^{-1} F_{k+1} + c_{N,k})^\top (f_{c,k+1} - \bar{f}_{c,k+1}) = 0. \qquad (43) $$

Then we can rearrange this to

$$ (-K^{-1}F_{k+1})^\top B_c^\top (f_{c,k+1} - \bar{f}_{c,k+1}) + c_{N,k}^\top (f_{c,k+1} - \bar{f}_{c,k+1}) = 0. \qquad (44) $$

The first member of the left hand side is zero as $B_c^\top (f_{c,k+1} - \bar{f}_{c,k+1}) = 0$. We are left with $c_{N,k}^\top (f_{c,k+1} - \bar{f}_{c,k+1}) = 0$, which implies that $f_{c,k+1}$ is unique when $c_{N,k} > 0$.

### 3.4.2 Numerical resolution of the linear complementarity problem

It is worth remarking that the above proof only demonstrates well-posedness for systems with $c_{N,k} > 0$, that is to say systems that start a time step with all contacts open. If a contact closes at the end of a time step, for the next time step we may only conclude that solutions exist, and nothing about their uniqueness. This is a well-known feature of quasi-static contact mechanics (Brogliato, 2016), however it is not of particular concern to us as the non-uniqueness is in the contact forces $f_c$, rather than the contact normal displacements $u_N$, which are the primary variable that we are concerned with for the purposes of our numerical homogenisation method.

In general, there are many solution methods to resolve linear complementarity problems. Of particular note is Lemke's algorithm (Lemke and Howson, 1964; Lemke, 1965), which is guaranteed to find a solution should one exist (and we have demonstrated that solutions always exist for our system, so this algorithm will always find one). However, this algorithm (and other classic algorithms for the linear complementarity problem such as the projected Gauss-Seidel method (Jourdan et al., 1998)) are generally quite slow. The method presented in this paper is developed with an eye towards its future implementation in the context of multiscale finite element simulations (simulations where the classical constitutive law at a quadrature point of the “macro-scale” finite element model is replaced with the homogenised results of a simulation conducted at the “micro-scale”), wherein many thousands of finite element simulations need to be completed rapidly to obtain macro-scale results in a timely manner. As such, we present a fast method of approximating the solution of the linear complementarity problem that gives us tolerable accuracy at sufficient speed.

The approximation algorithm involves a two-step iteration of resolving a given system, and then modifying the system if need be, and is described in Algorithm 2.
<caption>**Algorithm 2** Approximation algorithm for the linear complementarity problem.</caption>

**Require:** $K$, $B_c$, $F_{k+1}$, $c_{N,k}$
$u_{N,k+1} \leftarrow 0$
$Y_f = K^{-1}F_{k+1}$ &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ▷ Solved with LU factorisation
$c_f \leftarrow B_c Y_f$
$Y = K^{-1}B_c^T$ &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ▷ Solved with LU factorisation
$W \leftarrow B_c Y$
$f_{c,k+1} = W^{-1}(c_N - c_f)$ &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ▷ Solved with LU factorisation
$f_{c,k+1} \leftarrow \max(0, f_{c,k+1})$ &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ▷ Enforce complementarity
**return** $f_{c,k+1}$, $u_{N,k+1}$, $u_{k+1}$

After first initialising the system, the algorithm first obtains the linear solution of (40), that is solved to obtain the values of $f_{c,k+1}$ coherent with $u_{N,k+1} = 0$. This is achieved by solving a sequence of linear systems, rather than explicitly forming $K^{-1}$, which would be computationally prohibitive and numerically unstable for large-scale problems. At each iteration, the linearised equations are solved using an LU decomposition (specifically the Pardiso parallel sparse direct solver (Schenk and Gärtner, 2004)). Then, as this can potentially result in negative $f_{c,k+1}$ (that is to say, tensile contact forces that do not respect the complementarity condition), the subsequent line projects the contact forces back onto the positive orthant. At the same time, the contact displacements are updated iteratively to satisfy the complementarity condition. Using these results, a new value of the global contact force for each node $B_c f_{c,k+1}$ is obtained, from which the updated global displacements are calculated based on (16).

## 3.5 Periodic problem

In order to enforce periodic boundary conditions, the RVE boundary must be discretised using a mapped node formulation (Nguyen et al., 2012). Each mapped node pair consists of a node on the negative boundary $\partial V^-$ and its counterpart on the positive boundary $\partial V^+$. Consequently, the displacement vector $u_{k+1}$ in (38) is partitioned into three components (Nguyen et al., 2012): internal displacements $u_{\text{int}}^{k+1}$, positive boundary displacements $u_{+}^{k+1}$, and negative boundary displacements $u_{-}^{k+1}$, i.e.

$$u_{k+1} = \begin{bmatrix} u_{\text{int}}^{k+1} \\ u_{+}^{k+1} \\ u_{-}^{k+1} \end{bmatrix}, \qquad (45)$$

where we have moved the $k+1$ to the superscript position for the sake of legibility. According to (17), the displacements of the paired positive and negative boundary satisfy

$$u_{+}^{k+1,j} = u_{-}^{k+1,j} + \varepsilon(x_{+}^{k+1,j} - x_{-}^{k+1,j}), \quad j \in \{1, 2, \dots, p\} \qquad (46)$$

where $p$ represents the number of mapping node pairs belonging to the positive and negative boundaries. Therefore, (46) can be restated as

$$u_{+}^{k+1} = u_{-}^{k+1} + g^{k+1}, \qquad (47)$$

where $g^{k+1}$ is calculated from the macroscopic strain tensor $\varepsilon$, i.e. $g^{k+1} = \varepsilon x^{k+1}$ as shown in (17).

Accordingly, (38) can be restated by dividing the stiffness matrix $K$, external force vector $F_{k+1}$ and contact force $f_{c,k+1}$ in the form (Nguyen et al., 2012)

$$\begin{bmatrix} K_{\text{int,int}} & K_{\text{int,+}} & K_{\text{int,-}} \\ K_{+,\text{int}} & K_{+,+} & K_{+,,-} \\ K_{-,\text{int}} & K_{-,+} & K_{-,,-} \end{bmatrix} \begin{bmatrix} u_{\text{int}}^{k+1} \\ u_{+}^{k+1} \\ u_{-}^{k+1} \end{bmatrix} = \begin{bmatrix} F_{\text{int}}^{k+1} + f_{\text{int}}^{c,k+1} \\ F_{+}^{k+1} + f_{+}^{c,k+1} \\ F_{-}^{k+1} + f_{-}^{c,k+1} \end{bmatrix}. \qquad (48)$$

Substituting (47) into (48), $u_{+}^{k+1}$ terms are replaced by $u_{-}^{k+1}$ and $g^{k+1}$, yielding (Nguyen et al., 2012)

$$\begin{bmatrix} K_{\text{int,int}} & K_{\text{int,+}} + K_{\text{int,-}} \\ K_{+,\text{int}} + K_{-,\text{int}} & K_{+,+} + K_{+,,-} + K_{-,+} + K_{-,,-} \end{bmatrix} \begin{bmatrix} u_{\text{int}}^{k+1} \\ u_{-}^{k+1} \end{bmatrix} = \begin{bmatrix} F_{\text{int}} + f_{\text{int}}^{c,k+1} - K_{\text{int,+}}g^{k+1} \\ F_{+} + f_{+}^{c,k+1} + F_{-} + f_{-}^{c,k+1} - (K_{+,+} + K_{-,+})g^{k+1} \end{bmatrix}. \qquad (49)$$

This formulation systematically reduces the number of unknowns, enhancing computational efficiency while preserving kinematic constraints. The resulting linear system is also solved using the Pardiso parallel sparse direct solver (LU decomposition) (Schenk and Gärtner, 2004).
3.6 Perturbation method solution

In the workflow of Figure 2, four computational phases are executed within each time step. First, the displacement is simulated to update geometric configurations and stresses. Subsequent phases apply perturbation displacements (horizontal, vertical, and shear) to the updated geometry, assembling the configuration-dependent stiffness matrix D through the perturbation method. This enables automated updates of mechanical properties and stiffness matrices during deformation. Boundary conditions for perturbations are automatically detected based on initial configurations, as shown in Algorithm 3. Based on the input file, the top, bottom, right, and left boundaries, i.e., $B_t$, $B_b$, $B_r$ and $B_l$ are respectively automatically defined. Accurate detection requires a rectangular, axis-aligned domain, with all boundary nodes positioned precisely (within a specified tolerance) on their corresponding edges. Node mapping is one-dimensional, ensuring a consistent top-to-bottom and right-to-left correspondence. Positive designations are assigned to the top and right boundaries and negative designations to the bottom and left boundaries, with mapping points automatically identified. Small horizontal displacement differences are applied to the right boundary and its corresponding mapping points in the second computational phase, while small vertical displacement differences are applied to the top boundary and its mapping points in the third computational phase. To induce shear strain, horizontal and vertical displacement perturbations are applied to the top and right boundaries in the fourth computational phase, respectively. Under the periodic boundary conditions used for stiffness matrix calculations, the total number of nodal variables and linear equations during Newton iteration is determined. Throughout the deformation process, boundary conditions, mapping points, and the equation system for different displacement perturbations remain consistent.
<caption>**Algorithm 3 Automated boundary condition recognition algorithm**</caption>

**Require:** Initial geometry configuration $\boldsymbol{x}^0$, node coordinates
$x_{\min} \leftarrow \min(\boldsymbol{x}_1^0), x_{\max} \leftarrow \max(\boldsymbol{x}_1^0)$
$z_{\min} \leftarrow \min(\boldsymbol{x}_2^0), z_{\max} \leftarrow \max(\boldsymbol{x}_2^0)$
$\eta_x \leftarrow 10^{-8} \cdot (x_{\max} - x_{\min})$
$\eta_z \leftarrow 10^{-8} \cdot (z_{\max} - z_{\min})$
**for** each node $i$ with coordinates $(x_i, z_i)$ **do**
> **if** $|z_i - z_{\max}| < \eta_z$ **then**
> > $B_t \leftarrow B_t \cup \text{node } i$
> **else if** $|z_i - z_{\min}| < \eta_z$ **then**
> > $B_b \leftarrow B_b \cup \text{node } i$
> **end if**
> **end for**
> Assign positive designation to $B_t$ and $B_r$
> Assign negative designation to $B_b$ and $B_l$
> **for** each node $i \in B_t$ **do**
> > $j \leftarrow \arg\min_{j \in B_b} |x_i - x_j|$ 
> > Record mapping relationship $i \leftrightarrow j$
> > $M_{tb} \leftarrow M_{tb} \cup \{(i, j)\}$
> **end for**
> **for** each node $i \in B_r$ **do**
> > $j \leftarrow \arg\min_{j \in B_l} |z_i - z_j|$ 
> > Record mapping relationship $i \leftrightarrow j$
> > $M_{rl} \leftarrow M_{rl} \cup \{(i, j)\}$
> **end for**
> pass $\leftarrow$ true
> **if** $|B_t| = 0$ or $|B_b| = 0$ or $|B_r| = 0$ or $|B_l| = 0$ **then**
> > pass $\leftarrow$ false
> **end if**
> **if** $|B_t| \neq |B_b|$ or $|B_r| \neq |B_l|$ **then**
> > pass $\leftarrow$ false
> **end if**
> Check that all matches in $M_{tb}$ and $M_{rl}$ are unique
> **if** any node in any boundary has no mapping nodes **then**
> > pass $\leftarrow$ false
> **end if**
> **if** not pass **then**
> > $\eta_x \leftarrow 10 \eta_x$, $\eta_z \leftarrow 10 \eta_z$, and repeat classification
> > **if** $\eta_x > 10^{-3}(x_{\max} - x_{\min})$ **or** $\eta_z > 10^{-3}(z_{\max} - z_{\min})$ **then**
> > **Terminate with failure:** boundary recognition tolerance exceeded.
> > **end if**
> **end if**
> **return** $B_t, B_b, B_r, B_l$ and mapping sets $M_{tb}, M_{rl}$

# 4 Numerical simulations

In this section, we first conduct several simulations on simplified geometries to ensure that the contact solution and perturbation methods function as expected. We then consider two case studies, firstly a single-pore case for which analytical solutions exist to compare our results against, and secondly a multi-pore case in which we demonstrate the ability of our method to obtain homogenised equivalent elastic properties for geometries in which no analytical solution is known.
## 4.1 Model validation

To validate the accuracy of the solid contact model, a solid contact test is simulated, as shown in Figure 5(a). Two separate solids subjected to vertical normal stress $\sigma_3$ are analysed, to simulate contact-induced internal vertical stress distributions. The test is conducted using matching elements (with aligned nodes across the contact surface) in Figure 5(b) and mismatching elements (where nodes do not in general align across the contact surface) in Figure 5(c) while applying different vertical stress boundary conditions. The results demonstrate that the internal vertical stress computed by the simulator in Figure 2 aligns with the applied stress boundary conditions across the different models.

<figure><img src="image_6.png" /><figcaption>Figure 5: (a) Solid contact verification model with comparative analysis of (b) vertical stress distribution using matching elements across different solid parts and (c) vertical stress distribution using mismatching elements for different solid parts.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: Square solid model ($L = 0.01$ m, $E_s = 54.2$ GPa and $\nu_s = 0.163$ m) undergoing continuous deformations with displacement perturbation methodology for updated stiffness matrix calculation.</figcaption></figure>

In order to further verify the updated stiffness matrix obtained via the perturbation method, a simple square solid model with side length of 10 mm, elastic modulus $E_s = 54.2$ GPa, and Poisson's ratio $\nu_s = 0.163$ is simulated. The verification procedure consists of two steps, as shown in Figure 6. First, prescribed displacements are applied to the model, leading to a modified configuration. Then, based on the updated configuration, small displacement perturbations are introduced to generate small additional strains and compute the updated stiffness matrix. It should be noted that a single displacement perturbation does not necessarily correspond to a single strain component, since shear displacements may induce asymmetry in the model. However, the numerical results show that the stiffness matrix remains symmetric under displacements in different directions. The computed stiffness components are $D_{11} = D_{33} = 57.9$ GPa, $D_{21} = D_{31} = D_{13} = D_{23} = 11.3$ GPa, $D_{66} = 23.3$ GPa and the other components are 0. The
calculated stiffness matrix is consistent with the isotropic elastic stiffness matrix derived directly from $E_s$ and $\nu_s$, validating the effectiveness of the perturbation method in updating mechanical properties.

## 4.2 Single-pore case

The methodology developed in this paper is applied on a single-pore case study, examining how pore shape and large multidirectional displacement variations influence the stiffness matrix and poromechanical characteristics of the porous material. The initial configuration of the model consists of a square with side length of $L = 10$ mm, containing a centrally embedded ellipse characterised by semi-minor axis $a$ and semi-major axis $b$. This configuration ensures a symmetric initial model geometry where all finite elements possess positive volumes and the mesh is free from distortion or non-physical inversion, as illustrated in Figure 7. Porosity is quantified as $n = \pi ab/L^2$, maintained constant at $n = 0.1$ across all $a$ and $b$ values. Under the plane strain assumption, coordinate axes are defined with the $x_1$-axis ($x$) and $x_3$-axis ($z$) representing horizontal and vertical directions, respectively, while the $x_2$-axis ($y$) is perpendicular to the plane of analysis. The origin coincides with the geometric centre of the model. The pore aspect ratio $R$ is defined as $R = a/b$. The solid matrix exhibits elastic modulus $E_s = 54.2$ GPa and Poisson's ratio $\nu_s = 0.163$.

<figure><img src="image_8.png" /><figcaption>Figure 7: The initial configuration of the model consisting of a square with side length of $L = 10$ mm and an ellipse with minor semi-axis $a$ and major semi-axis $b$.</figcaption></figure>

### 4.2.1 Comparison with analytical solution

At the initial stage, in the absence of continuous deformations, anisotropic macroscopic material properties arise due to the presence of elliptical pores. For the 2D plane-strain model, directional Young's moduli $E_1$ and $E_3$, shear modulus $G_{13}$ and Poisson's ratio $\nu_{13}$ are analytically determined as functions of porosity $n$ and semi-axes $a$ and $b$. In the constitutive relationship $\epsilon = C\sigma$ where $C$ is the elastic compliance matrix (equivalent to $D^{-1}$), we have (Kachanov
et al., 1994; Tsukrov and Kachanov, 2000)

$$C_{11} = \frac{1}{E_{s1}} + \frac{\pi L_{2D}}{A\sqrt{E_{s1}}} \left[ (a^2 - b^2) \cos^2 \alpha + b^2 + \frac{ab}{L_{2D}\sqrt{E_{s1}}} \right], \quad (50)$$

$$C_{13} = C_{31} = -\frac{\nu_{s13}}{E_{s1}} - \frac{\pi ab}{A\sqrt{E_{s1}E_{s3}}}, \quad (51)$$

$$C_{33} = \frac{1}{E_{s3}} + \frac{\pi L_{2D}}{A\sqrt{E_{s3}}} \left[ (b^2 - a^2) \cos^2 \alpha + a^2 + \frac{ab}{L_{2D}\sqrt{E_{s3}}} \right], \quad (52)$$

$$C_{66} = \frac{1}{G_{s13}} + \frac{\pi L_{2D}}{A\sqrt{E_{s1}E_{s3}}} \left[ (b^2 - a^2) \left( \sqrt{E_{s3}} - \sqrt{E_{s1}} \right) \cos^2 \alpha + a^2 E_{s3} + b^2 \sqrt{E_{s1}} + abL_{2D}\sqrt{E_{s1}E_{s3}} \right], \quad (53)$$

where

$$L_{2D} = \sqrt{\frac{1}{G_{s13}} - \frac{2\nu_{s13}}{E_{s1}} + \frac{2}{\sqrt{E_{s1}E_{s3}}}}, \quad (54)$$

and $E_{s1}$, $E_{s3}$, $G_{s13}$ and $\nu_{s13}$ denote the Young's modulus in the $x_1$ and $x_3$ directions, shear modulus and Poisson's ratio of the solid matrix without pores, respectively. $\alpha$ denotes the angle between the major axis of the ellipse and the horizontal direction, while $A$ represents the representative area of the whole volume. For an isotropic solid matrix with zero pore orientation angle ($\alpha = 0$), we obtain

$$E_1^{2D} = \frac{E_s}{1 + n(1 + 2R)}, \quad (55)$$

$$E_3^{2D} = \frac{E_s}{1 + n(1 + \frac{2}{R})}, \quad (56)$$

$$\nu_{13}^{2D} = \frac{\nu_{s13} E_1^{2D}}{E_s} + \frac{n E_1^{2D}}{E_s}, \quad (57)$$

$$G_{13}^{2D} = \frac{1}{\frac{1}{G_{s13}} + \frac{2n(a+b)^2}{E_s ab}}, \quad (58)$$

where the superscript 2D represents plane-strain adjustments. Specifically, these parameters are adjusted from their three-dimensional counterparts by constraining the out-of-plane strain component to zero ($\epsilon_2 = 0$), which modifies the stress-strain relationships and leads to the apparent stiffness change typical of plane-strain formulations.

<figure><img src="image_9.png" /><figcaption>Figure 8: Comparison of (a) Young's modulus in the $x_1$ and $x_3$ directions, shear modulus and (b) Poisson's ratio between the numerical solution and the analytical solution (Thorpe and Sen, 1985; Kachanov et al., 1994).</figcaption></figure>

Based on (56)–(58), the mechanical parameters can be derived under 2D plane-strain conditions. However, the present model employs three-dimensional material parameters, making direct comparison between the 2D and 3D formulations inconsistent. To establish a proper correspondence, the 2D plane-strain predictions are reconciled with the 3D elasticity theory, ensuring that the effective in-plane mechanical responses under plane-strain conditions remain consistent with those obtained from the full three-dimensional constitutive relations. Accordingly, the elastic constants are modified following (Thorpe and Sen, 1985; Kachanov et al., 1994)

$$E_1^{2D} = \frac{E_1}{1 - \nu_{23}\nu_{32}}, \quad (59)$$

$$E_3^{2D} = \frac{E_3}{1 - \nu_{21}\nu_{12}}, \quad (60)$$

$$\nu_{13}^{2D} = \frac{\nu_{13} + \nu_{12}\nu_{23}}{1 - \nu_{23}\nu_{32}}. \quad (61)$$
Combining with (20), we obtain

$$E_3 = \frac{E_2^2 E_3^{2D} - E_2^{2D} E_3^{2D} E_2 \nu_s^2}{E_2^2 - E_1^{2D} E_3^{2D} \nu_s^4}, \qquad (62)$$

$$E_1 = \frac{E_1^{2D} (E_2 - E_3 \nu_s^2)}{E_2}, \qquad (63)$$

$$\nu_{13} = \frac{\nu_{13}^{2D}}{E_1^{2D}} E_1 - \frac{\nu_s^2}{E_2} E_1. \qquad (64)$$

The analytical solution of the elastic constants is compared with the numerical results, as shown in Figure 8. It is noted that for pore aspect ratios $R > 0.2$, the analytical and numerical solutions exhibit strong agreement, confirming the validity of the numerical framework. However, discrepancies emerge for smaller $R$ values (especially $R = 0.15$). This deviation arises because the mesh quality deterioration and shape distortion arise with the small initial pore aperture.

### 4.2.2 Stiffness matrix updating after deformation

The methodologies are applied to analyse multidirectional deformations in the model with an initial pore aspect ratio $R = 0.8$, focusing on the evolving stiffness matrix and mechanical properties in updated configurations. Three distinct loading scenarios are implemented so the geometric but not constitutive nonlinearity is modeled. First, a horizontal displacement difference rate ($u_1 = -1 \times 10^{-2}$ mm/s) is imposed on the lateral boundaries, with the right boundary designated as negative displacement (that is to say we compress the RVE) and zero vertical displacements on the boundaries. The configuration undergoes compression until solid contact occurs between the opposite sides of the pore, accompanied by evolving internal contact stresses. Second, a vertical displacement rate ($u_3 = -1 \times 10^{-2}$ mm/s) is applied to the top and bottom boundaries, with the top boundary designated as negative displacement and zero horizontal displacement, inducing vertical strain. Third, simultaneous displacement rates ($u_{13} = 1 \times 10^{-2}$ mm/s vertically on the right and left boundaries and $u_{31} = 1 \times 10^{-2}$ mm/s horizontally on the top and bottom boundaries) generate shear displacement and strain.

After compressing and shearing the model over time in different directions, the pore in the solid matrix is compressed, resulting in a decrease in total porosity, as shown in Figure 9(a). The deformation behaviour exhibits pronounced nonlinearity, primarily due to geometric effects as the pore shape changes. For the initial $R = 0.8$, the porosity reduction under shear displacement initially lags behind that induced by horizontal and vertical displacements, but then accelerates rapidly, reaching effectively zero at a cumulative displacement of $u_{13} = 3.1$ mm. Porosity reduction under horizontal compression proceeds more slowly than under vertical compression, with pore closure occurring at $u_1 = 4.6$ mm versus $u_3 = 3.6$ mm for the vertical loading. After solid contact, the porosity remains near to zero (fully closed), confirming the consistency of the compaction behaviour.

During the compression process, the stress distribution evolves significantly, which can be divided into contactless and contact stages, as shown in Figure 9(b-c). Due to the initial pore aspect ratio $R = 0.8$, vertical compression reduces the pore aspect ratio, resulting in a more elongated pore shape, while horizontal compression initially increases the aspect ratio, first leading to a more uniform pore shape and later an elongated pore along the vertical direction with $R > 1.0$. Stress concentrations develop at the pore tips and the contacting surfaces, particularly in circumferential directions, intensifying with increasing strain. Tensile stress concentration dominates near pore midlines in the circumferential ($x$) direction, whereas compressive stress concentration prevails at pore tips in the circumferential ($z$) direction. Shear stress changes remain small compared to horizontal and vertical stresses. Post-contact, normal contact tractions on the internal interfaces prevent the overlap of internal boundaries; however, their magnitudes remain minimal and negligible relative to stresses within the solid matrix.
<figure><img src="image_10.png" /><figcaption>Figure 9: Evolution of (a) porosity, (b) horizontal normal stress $\sigma_1$, vertical normal stress $\sigma_3$, and shear stress $\sigma_6$, before contact, and (c) stress states after contact during progressive multidirectional compression loading.</figcaption></figure>

The evolving model configuration necessitates on-the-fly updates to the stiffness matrix and mechanical properties in accordance with (50)–(53). Figure 10 illustrates the evolution of Young’s modulus ($E_1$ and $E_3$), shear modulus ($G_{13}$), Poisson’s ratio ($\nu_{13}$), and Biot coefficients ($b_1$, $b_2$, and $b_3$) during multidirectional loading. Under uniaxial compression, Young’s modulus decreases in the loading direction while increasing linearly in the orthogonal direction. For instance, as the applied horizontal displacement $u_1$ increases to 4.64 mm, $E_1$ decreases linearly from 44.22 GPa to 40.95 GPa, while $E_3$ increases linearly from 48.8 GPa to 55.58 GPa, approximately reaching $E_s$ (54.2 GPa). When the applied vertical displacement $u_3$ increases to 3.6 mm, $E_3$ decreases linearly from 41.04 GPa to 37.58 GPa, while $E_1$ increases linearly from 44.28 GPa to 55.68 GPa (approximately $E_s$). This directional trend results from the geometric distortion of the elliptical pore. Compression reduces the effective load-bearing area along the loading direction while simultaneously elongating the pore laterally, which stiffens the material in the orthogonal direction. This trend is consistent with the analytical expression in (55)–(56), in which the effective Young’s modulus depend on the pore aspect ratio $R$. As $R$ increases, the modulus along the minor-axis direction increases, whereas the modulus along the major-axis direction decreases. Following solid contact, the Young’s moduli in different directions under horizontal or vertical compression progressively approach the solid stiffness $E_s$, particularly in the loading direction, where a pronounced jump increase is observed. In contrast, shear loading leads to a reduction in the Young’s modulus in all directions before solid contact. Under combined shear loading (simultaneous applied displacements of 3.1 mm on
the right and left boundaries and applied horizontal displacements of 3.1 mm on the top and bottom boundaries), both $E_1$ and $E_3$ decrease, from 44.22 GPa to 37.6 GPa and 41.05 GPa to 33.73 GPa, respectively. Post-contact, a further increase in displacement to 3.4 mm reverses this trend, driving $E_1$ and $E_3$ up to 47.58 GPa and 47.11 GPa. Nevertheless, both moduli remain lower than the stiffness of the intact solid, $E_s$.

<figure><img src="image_11.png" /><figcaption>Figure 10: Evolution of the effective system parameters under vertical, horizontal and shear deformations. (a) Young's modulus $E_1$, (b) Young's modulus $E_3$, (c) shear modulus $G_{13}$, (d) Poisson's ratio $\nu_{13}$, (e) Biot coefficient in different direction $b_1$, $b_2$ and $b_3$, and (f) Shear coupling moduli $D_{16}$, $D_{26}$ and $D_{36}$.</figcaption></figure>

As compression or shear displacement in different directions increases over time, the shear modulus $G_{13}$ increases, while Poisson's ratio $\nu_{13}$ decreases before solid contact. Horizontal compression exerts the most significant influence, increasing $G_{13}$ by 3.13 GPa until solid contact, compared to 2.30 GPa under vertical compression and 0.22 GPa under shear displacement, with the growth rate accelerating over time. However, after solid contact and without considering friction, $G_{13}$ exhibits an abrupt increase during different loading form, approaching the solid shear modulus $G_s$ during compression process but reaching about 20.8 GPa during shearing process which remains slightly lower than $G_s$.

Poisson's ratio $\nu_{13}$ exhibits inverse proportionality to displacement magnitude, with shear displacement causing the most pronounced reduction rate and vertical compression the least, as shown in Figure 10(d). During pre-contact deformation, $\nu_{13}$ decreases by 0.026, 0.068, and 0.045 until solid contact under horizontal, vertical, and shear dis-
placements, respectively. Following contact existing, the evolution of $\nu_{13}$ diverges depending on the loading mode. Under horizontal compression, $\nu_{13}$ exhibits an abrupt increase and rapidly approaches the value under vertical loading, which is close to solid Poisson's ratio $\nu_s$. In contrast, under shear loading, $\nu_{13}$ overshoots $\nu_s$ immediately after solid contact and subsequently decreases as the contacting pore surfaces undergo rotational rearrangement and sliding while maintaining contact.

Biot coefficients demonstrate directionally dependent behaviour, as presented in Figure 10(e). The transverse coefficient $b_2$ (in the $x_2$-direction) decreases monotonically across all loading modes. In contrast, the axial coefficients $b_1$ (in the $x_1$-direction) and $b_3$ (in the $x_3$-direction) exhibit opposite responses. Uniaxial compression increases the Biot coefficient in the loading direction while reducing it in the orthogonal direction, whereas shear loading elevates both $b_1$ and $b_3$. Following solid contact, the Biot coefficients exhibit a pronounced abrupt decrease, particularly under horizontal and vertical compression, where they progressively approach zero. These directional variations are inversely correlated with the evolution of Young's modulus, thereby underscoring the significance of elastic coupling effects on the evolution of poroelastic properties. In addition, progressive shear displacement induces a transition from a symmetric anisotropic configuration to an asymmetric state, as shown in Figure 9(b), manifesting as non-zero shear-coupling stiffness coefficients $D_{16}$, $D_{26}$, and $D_{36}$ in (18). These parameters, which characterise normal-shear strain coupling effects critical to constitutive behaviour of geomaterials, remain zero during purely horizontal or vertical compression where model symmetry is preserved. Figure 10(f) shows the evolution of these parameters, $D_{16}$, $D_{26}$, and $D_{36}$, during shear displacement. They exhibit near-linear increases from 0 to 6.85 GPa, 2.19 GPa, and 6.60 GPa, respectively, at a shear displacement of 3.1 mm. Following internal boundary contact, $D_{16}$, $D_{26}$, and $D_{36}$ exhibit pronounces jump decrease to less than 0.5 GPa. This shear-induced coupling, typically neglected in conventional stiffness formulations, underscores the necessity of constitutive models with nonlinear elasticity at the RVE scale for anisotropic media undergoing multidirectional loading.

However, as observed, slight oscillations appear in the evolution of the stiffness parameters after solid contact is established. This behaviour is mainly related to the numerical treatment of the contact constraints rather than to any physical instability. Since $\mathbf{W}$ in (41) is symmetric positive semidefinite, the associated quadratic programming problem is convex. Therefore, the exact LCP solution guarantees the uniqueness of the normal displacement $\mathbf{u}_N$. Nevertheless, the corresponding contact force $\mathbf{f}_c$ may not be unique when the solution lies on a non-smooth face of the feasible set, such as in situations involving nearly simultaneous multi-contact or transitions between active constraint sets. In these cases, small perturbations in the active set identification or numerical tolerances may lead to slight variations in the computed contact forces, which in turn induce minor oscillations in the tangent stiffness parameters. We emphasise that these oscillations are purely numerical artefacts. In our simulations, the deviations remained below the prescribed tolerance (generally less than $1 \times 10^{-4}$) within the numerical calculation workflow shown in Figure 2. Moreover, these small discrepancies did not propagate to the macroscopic level and had a negligible influence on the homogenised stress response.

## 4.3 Multi-pore case

The developed model and simulator are further applied to a two-dimensional square domain with side length $L = 10$ mm, in which multiple pores of random sizes and spatial distributions are embedded, as shown in Figure 11. The initial porosity is defined as the area fraction $n = \sum(\pi r_i^2)/L^2$, where $r_i$ denotes the radius of the $i$-th pore. The solid skeleton is also assumed to be isotropic and linearly elastic, characterised by a Young's modulus of $E_s = 54.2$ GPa and Poisson's ratio $\nu_s = 0.163$. The displacement loading scheme and periodic boundary conditions are implemented in the same manner as in the single-pore case. This configuration enables the simulation of pore-matrix interactions during progressive loading and provides a framework for examining how staged pore contact influences the evolution of the global stiffness matrix and Biot coefficients. By subjecting the model to vertical, horizontal, and shear displacements, the study captures not only the porosity-driven stiffening process but also the anisotropic effects arising from sequential pore contact and closure.
<figure><img src="image_12.png" /><figcaption>Figure 11: The initial configuration of the model consisting of multiple pores of random sizes and spatial distributions. All dimensions are in units of millimetres.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 12: Pore contact behaviour during progressive (a) vertical displacement, (b) horizontal displacement, and (c) shear displacement.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 13: Evolution of (a) effective porosity, (b) Young's modulus $E_1$, (c) Young's modulus $E_3$, (d) shear modulus $G_{13}$, (e) Poisson's ratio $\nu_{13}$, and (f) Biot coefficient in different direction $b_1$, $b_2$ and $b_3$, during progressive vertical, horizontal and shear displacements.</figcaption></figure>

Figure 12 and Figure 13 show the evolution process of different mechanical properties and pore contact behaviour under vertical, horizontal, and shear compression. The results reveal a strongly nonlinear and path-dependent poromechanical evolution under these loading conditions. As external displacement increases, pores undergo a staged transition from gradual flattening to localised pore wall contact and ultimately to full closure, as summarised in Table 1.

<table><caption>Table 1: Corresponding displacement (mm) associated with progressive pore contact and complete closure under different loading modes.</caption><thead><tr><th>Loading mode</th><th>1st contact</th><th>2nd contact</th><th>3rd contact</th><th>4th contact</th><th>5th contact</th><th>6th contact</th><th>Fully closed</th></tr></thead><tbody><tr><td>Vertical</td><td>2.9</td><td>3.67</td><td>3.7</td><td>4.19</td><td>4.21</td><td>4.32</td><td>5.07</td></tr><tr><td>Horizontal</td><td>3.2</td><td>4.0</td><td>4.1</td><td>4.4</td><td>4.5</td><td>4.8</td><td>5.07</td></tr><tr><td>Shearing</td><td>3.0</td><td>3.11</td><td>3.27</td><td>3.82</td><td>4.68</td><td>4.9</td><td>5.23</td></tr></tbody></table>

The transition between the different loading stages is governed by pore size and pore position, with smaller pores closing earlier and larger pores closing later. For normal loading, pores closer to the loading axis (along the normal
direction of the applied load boundary) consistently initiate contact earlier, whereas for shear loading, pores closer to the compressed angles initiate contact earlier while the middle pores close later. The effective porosity evolution therefore follows a nonlinear trajectory, with an initial acceleration due to pore flattening, followed by a deceleration as multiple pore contacts emerge, and finally a slow approach to zero porosity once the contact network saturates. This staged contact sequence highlights that porosity alone is consistent with earlier findings in micromechanics; rather, the connectivity and topology of the emerging contact skeleton exert a decisive influence on the macroscopic mechanical response.

The staged contact mechanism influences the evolution of effective stiffness and poromechanical coefficients, with responses differing markedly between normal (vertical and horizontal) compression and shear loading. These trends are broadly consistent with those identified in the above single-pore case, though additional effects emerge from the collective interactions among multiple pores. Specifically, before solid contact, the effective Young's modulus in the loading direction decreases continuously with displacement, whereas the modulus in the orthogonal direction increases as pores are progressively flattened. This behaviour is due to the geometric evolution of pores and the corresponding load transfer modes. For example, under vertical compression, the vertical dimension of pores shrinks, reducing the solid thickness around pores available to sustain axial loads. As these load-bearing paths become thinner, deformation is increasingly bending-dominated, leading to a rapid decline in vertical stiffness. Conversely, in the orthogonal direction, pore flattening increases the continuous solid fraction, enhancing stretching resistance and thereby raising the modulus. Once contacts occur, the increase in orthogonal stiffness becomes significantly attenuated, as the strengthening mechanism reaches a saturation point close to $E_s$. In contrast, the loading-direction stiffness transfers to increase and eventually converges to a value comparable to that in the orthogonal direction. Specifically, each successive staged contact of different pores induces a discrete jump in stiffness along the loading direction. In addition, shear displacement before solid contact leads to a continuous decrease in the Young's modulus in different directions. This behaviour mirrors the single-pore case, but the effect is amplified in multi-pore systems. Under shear, pores deform asymmetrically and undergo relative sliding, which disrupts load-bearing continuity across the matrix. Once solid contact is formed, however, stress can be directly transmitted across contacting regions, resulting in a pronounced increase in stiffness in all directions. However, because contact friction is neglected in the present model, the stiffness cannot fully recover to that of the intact solid, even under conditions of complete pore contact.

The evolution of the shear modulus ($G_{13}$) and Poisson's ratio ($\nu_{13}$) under multidirectional loading exhibits a distinct behaviour compared to the effective Young's moduli. In any compression mode, $G_{13}$ initially increases slowly as pores are progressively flattened, since oblique load-transfer paths within the solid skeleton become more effective in sustaining shear resistance. Further, after staged pore contact, $G_{13}$ increases markedly and continues until full contact is achieved. This enhancement arises because contacting pore walls tend to rotate relative to one another under shear loading; such rotational interactions intensify stress transmission across the contact interfaces, thereby strengthening the shear response. The response of Poisson's ratio ($\nu_{13}$) further highlights the anisotropic coupling. Under horizontal or vertical compression, $\nu_{13}$ decreases, especially in the later loading, reflecting constrained orthogonal contraction and diminished transverse-longitudinal coupling due to pore geometry. Once pore contacts are established, this trend reverses, and $\nu_{13}$ progressively approaches the solid Poisson's ratio $\nu_s$ as full contact is reached. The stiffening effect associated with solid contact is more pronounced under horizontal compression than under vertical compression. In contrast, under shear compression, $\nu_{13}$ increases rapidly upon contact formation but subsequently exhibits a slight decrease as rotational interactions and surface sliding between contacting pore walls develop. Overall, these observations indicate a fundamental reorganisation of load-transfer mechanisms, shifting from a skeleton-controlled regime to one dominated by the pore-contact network, consistent with the trends identified in the single-pore system.

In addition, the Biot coefficients exhibit analogous directional dependence. Under vertical or horizontal compression, pore flattening increases the effective solid fraction in the loading direction while reducing it in the orthogonal direction. Consequently, the Biot coefficient increases along the loading axis, indicating stronger coupling between pore pressure and normal stress, while it decreases transversely. After staged solid contact, the Biot coefficients in different directions during different-direction compression decrease greatly and tend to zero after full contact, denoting the decreasing effect of pore fluid pressure. In shear loading, however, the Biot coefficients in both vertical and horizontal directions increase. This is because diagonal pore compression enhances the effective load-bearing area, thereby strengthening the coupling between fluid pressure and stress in multiple orientations. Nevertheless, the Biot coefficients still exhibit abrupt decreases during staged pore contact and remain a non-zero residual value after full contact.

In summary, the findings highlight pore contact as the dominant mechanism governing the nonlinear and anisotropic poromechanical response of porous solids. The sequential nature of pore closure not only controls the rate of stiffness evolution but also determines its directional dependence. This behaviour reflects a fundamental transition in deformation mechanisms from a regime dominated by distributed solid-skeleton deformation to one governed by localised contact networks, accompanied by stress redistribution and progressive development of an anisotropic microstructural fabric.
# 5 Conclusion

In this study, a comprehensive FEM-based numerical framework is developed to analyse the behaviour of porous solids with frictionless solid contact in two dimensions. The framework incorporates a perturbation-based stiffness matrix updating scheme, enabling on-the-fly tracking of evolving anisotropic mechanical properties and poroelastic coefficients under multidirectional loading. The methodology is applied to both single-pore and multi-pore configurations to investigate the interplay between pore geometry, contact evolution, and macroscopic mechanical response. The main findings are summarised as follows:

1. The proposed framework develops a method of fast approximation to the solution of the underlying linear complementarity problem, with an accuracy that is sufficient for its intended applications in homogenization simulations. This framework also allows automated updates of stiffness parameters for evolving configurations, accurately capturing complex contact interactions and predicting anisotropic Young's modulus, shear modulus, Poisson's ratio, and Biot coefficients under multidirectional loading. This method enables homogenisation of geometries that lack analytical solutions due to the strong interactions between the pores.

2. In both single- and multi-pore systems, progressive pore closure dominates the evolution of macroscopic mechanical properties. Sequential contact induces a transition from a skeleton-controlled deformation regime to a contact-network-dominated behaviour, highlighting that initial porosity alone cannot fully describe the material state.

3. Before solid contact, Young's modulus along the loading direction decreases with displacement, while the orthogonal modulus generally increases due to pore flattening. Shear loading induces asymmetric deformation and reduces Young's modulus in all directions. After solid contact, the Young's moduli in different directions increase markedly, particularly under horizontal and vertical compression, where they progressively approach the Young's modulus of the intact solid.

4. Shear modulus initially increases during pore flattening due to enhanced oblique load transfer, and approaches to solid shear modulus after contact occurs. Poisson's ratio exhibits decreases with displacement under compression, but can increase post-contact, especially under shear loading, reflecting inclined contact surfaces facilitating stress transfer and anisotropic coupling effects.

5. Biot coefficients exhibit strong directional dependence. They increase along the loading axis under uniaxial compression and rise in multiple directions under shear loading. Following solid contact, the Biot coefficients decrease greatly, particularly under compressive loading where they approach zero. This behaviour underscores the critical role of pore geometry, contact sequence, and emerging contact networks in governing poroelastic behaviour

Despite the developments and results presented above, several limitations remain in the present model and numerical methodology. First, the current implementation is restricted to a 2D plane strain setting. In particular, although solid contact can be robustly handled in 2D, the extension to face-to-face contact in 3D configurations introduces significant algorithmic and geometrical challenges. These include more complicated contact detection, constraint enforcement, and active-set identification, which have not yet been developed in the present work. Second, frictional effects within the contact group are neglected, and only frictionless normal contact is considered. While this assumption simplifies the complementarity structure and improves computational efficiency, it may underestimate shear resistance and energy dissipation in granular or rough-contact systems. Third, the solid matrix is modelled as linear elastic, while it does not account for material nonlinearities such as plasticity, damage, or viscoelastic effects. In addition, the pore fluid pressure along the pore boundaries is assumed to be constant. Moreover, the computational domain is restricted to an axis-aligned rectangular RVE. This requirement facilitates the implementation of periodic boundary conditions and simplifies matrix assembly, but it reduces geometric flexibility.

From a numerical perspective, the objective of the solid contact treatment is to solve (16), which involves two challenges. First, the system contains two coupled unknown fields (i.e. the displacements and the contact forces), requiring a complex optimisation process. Second, determining the active contact set requires checking all possible contact groups, which becomes computationally expensive for large-scale problems due to the combinatorial growth of candidate contact configurations. A major contribution of this work is the significant speedup achieved by the proposed fast approximation LCP method. However, for cases involving complex and densely distributed contact groups, slight oscillations in the post-contact stiffness parameters may occur. These oscillations are mainly induced by the approximate treatment of the active set in highly constrained configurations. Although they remain small and within the prescribed numerical tolerance, this phenomenon can limit the utility of the developed method where a high degree of accuracy is required.

The framework presented in this paper can be expanded in a number of different directions in the future. Most obviously, an extension to three dimensional finite element simulations is eminently achievable, without any substantial changes needing to be made. A small caveat to the system's applicability to geomaterials is that it requires the
solid material to be fully connected (and thus incapable of rigid body motions after the Dirichlet boundary conditions have been applied). However, in realistic geometries of rocks and soils it is possible to have disconnected grains that are individually capable of rigid body motions, which ultimately leads to the LU decomposition not being obtainable and the method failing. The remedy in this case is straightforward, by simply changing from quasi-statics to dynamics and specifying an impact law (typically Newton's impact law), a classical nonsmooth contact dynamics method is obtained (Dubois et al., 2018), which is easily adapted to our framework. Staying within two dimensional analysis, we may also note that the Coulomb friction model can be written in LCP form, following the classic approach of Stewart and Trinkle (1996), allowing a straightforward extension to the more realistic case of frictional contact. Cohesive zone models for fracture propagation (possibly inclusive of friction) can also naturally be implemented in this framework (Monerie and Acary, 2001; Collins-Craft et al., 2021, 2022; Collins-Craft and Acary, 2025). It is only in the case of friction (or cohesion) in three dimensions, where the problem cannot be written in LCP form, but instead as a second-order cone complementarity problem (Acary et al., 2024; Acary and Collins-Craft, 2025), where the method as presented in this paper to approximate the solution of the contact problem will not work, and further effort is required to develop a similar fast approximation method. We note that the numerical homogenisation by perturbation method remains robust in all of these cases, with the small caveat that formally as dissipative behaviours frictional sliding and crack growth should not be homogenised into elastic constants, but instead should appear as macro-scale plasticity.

The framework presented in this article provides a robust and efficient tool for resolving complex contact interactions, capturing nonlinear, anisotropic, and path-dependent poromechanical responses that conventional FEM or continuum models cannot easily predict. By linking pore-scale mechanisms to macroscopic behaviour, it offers critical insights into stress redistribution, stiffness evolution, and poroelastic coupling, with direct relevance to both engineered porous materials and natural geological formations under multidirectional stress, laying the foundation for future studies on the design, optimisation, and mechanical stability assessment of porous systems.

## Acknowledgements

This research is supported by French National Research Agency (Grant ANR-22-CPJ1-0027-01 Geosolutions to meet the climate change).

## Statements and Declarations

The authors declare that they have no known competing financial interests that could have appeared to influence the work reported in this paper.

## References

Acary, V., F. Bourrier and B. Viano (1st Sept. 2023). "Variational Approach for Nonsmooth Elasto-Plastic Dynamics with Contact and Impacts". In: *Computer Methods in Applied Mechanics and Engineering* 414, p. 116156. ISSN: 0045-7825. DOI: 10.1016/j.cma.2023.116156. URL: https://www.sciencedirect.com/science/article/pii/S0045782523002803.

Acary, V. and N.A. Collins-Craft (2025). "On the Moreau-Jean Scheme with the Frémond Impact Law: Energy Conservation and Dissipation Properties for Elastodynamics with Contact, Impact and Friction". In: *Journal of Theoretical, Computational and Applied Mechanics*, pp. 1-29. DOI: 10.46298/jtcam.13480. URL: https://hal.science/hal-04230941.

Acary, V. et al. (16th Jan. 2024). "Second-Order Cone Programming for Frictional Contact Mechanics Using Interior Point Algorithm". In: *Optimization Methods and Software*, pp. 1-30. ISSN: 1055-6788, 1029-4937. DOI: 10.1080/10556788.2023.2296438. URL: https://www.tandfonline.com/doi/full/10.1080/10556788.2023.2296438.

Ahmed, S. et al. (Feb. 2019). "Drained Pore Modulus and Biot Coefficient from Pore-Scale Digital Rock Simulations". In: *International Journal of Rock Mechanics and Mining Sciences* 114, pp. 62-70. ISSN: 13651609. DOI: 10.1016/j.ijrmms.2018.12.019. URL: https://linkinghub.elsevier.com/retrieve/pii/S1365160918300959.

Antunes, F.V. et al. (Mar. 2014). "Numerical Study of Contact Forces for Crack Closure Analysis". In: *International Journal of Solids and Structures* 51.6, pp. 1330-1339. ISSN: 00207683. DOI: 10.1016/j.ijsolstr.2013.12.026. URL: https://linkinghub.elsevier.com/retrieve/pii/S0020768313005076.

Batistić, I., P. Cardiff and Ž. Tuković (Jan. 2022). "A Finite Volume Penalty Based Segment-to-Segment Method for Frictional Contact Problems". In: *Applied Mathematical Modelling* 101, pp. 673-693. ISSN: 0307-904X. DOI: 10.1016/j.apm.2021.09.009. URL: https://linkinghub.elsevier.com/retrieve/pii/S0307904X21004248.

Batzle, M.L., G. Simmons and R.W. Siegfried (1980). "Microcrack Closure in Rocks under Stress: Direct Observation". In: *Journal of Geophysical Research: Solid Earth* 85.B12, pp. 7072-7090. ISSN: 2156-2202. DOI: 10.1029/JB085iB12p07072. URL: https://onlinelibrary.wiley.com/doi/abs/10.1029/JB085iB12p07072.
Biot, M.A. (1st Feb. 1941). "General Theory of Three-Dimensional Consolidation". In: *Journal of Applied Physics* 12.2, pp. 155-164. ISSN: 0021-8979, 1089-7550. DOI: 10.1063/1.1712886. URL: https://pubs.aip.org/jap/article/12/2/155/141365/General-Theory-of-Three-Dimensional-Consolidation.

— (1973). "Nonlinear and Semilinear Rheology of Porous Solids". In: *Journal of Geophysical Research (1896-1977)* 78.23, pp. 4924-4937. ISSN: 2156-2202. DOI: 10.1029/JB078i023p04924. URL: https://onlinelibrary.wiley.com/doi/abs/10.1029/JB078i023p04924.

Blenkinsop, T. (2002). *Deformation Microstructures and Mechanisms in Minerals and Rocks*. Dordrecht: Kluwer Academic Publishers. ISBN: 978-0-412-73480-9. DOI: 10.1007/0-306-47543-X. URL: http://link.springer.com/10.1007/0-306-47543-X.

Bouhlel, M., M. Jamei and C. Geindreau (Mar. 2010). "Microstructural Effects on the Overall Poroelastic Properties of Saturated Porous Media". In: *Modelling and Simulation in Materials Science and Engineering* 18.4, p. 045009. ISSN: 0965-0393. DOI: 10.1088/0965-0393/18/4/045009. URL: https://doi.org/10.1088/0965-0393/18/4/045009.

Brogliato, B. (2016). *Nonsmooth Mechanics: Models, Dynamics and Control*. Ed. by A. Isidori et al. Communications and Control Engineering. Cham: Springer International Publishing. ISBN: 978-3-319-28662-4 978-3-319-28664-8. DOI: 10.1007/978-3-319-28664-8. URL: http://link.springer.com/10.1007/978-3-319-28664-8.

Burman, E., P. Hansbo and M.G. Larson (May 2023). "The Augmented Lagrangian Method as a Framework for Stabilised Methods in Computational Mechanics". In: *Archives of Computational Methods in Engineering* 30.4, pp. 2579-2604. ISSN: 1134-3060, 1886-1784. DOI: 10.1007/s11831-022-09878-6. URL: https://link.springer.com/10.1007/s11831-022-09878-6.

Capobianco, G. and S.R. Eugster (2016). "A Moreau-type Variational Integrator". In: *PAMM* 16.1, pp. 941-944. ISSN: 1617-7061. DOI: 10.1002/pamm.201610453. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/pamm.201610453.

— (2018). "Time Finite Element Based Moreau-type Integrators". In: *International Journal for Numerical Methods in Engineering* 114.3, pp. 215-231. ISSN: 1097-0207. DOI: 10.1002/nme.5741. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.5741.

Capobianco, G., S.R. Eugster and R.I. Leine (2017). "Moreau-Type Integrators Based on the Time Finite Element Discretization of the Virtual Action". In: *PAMM* 17.1, pp. 143-144. ISSN: 1617-7061. DOI: 10.1002/pamm.201710041. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/pamm.201710041.

Capobianco, G., J. Harsch and S. Leyendecker (1st Jan. 2024). "Lobatto-Type Variational Integrators for Mechanical Systems with Frictional Contact". In: *Computer Methods in Applied Mechanics and Engineering* 418, p. 116496. ISSN: 0045-7825. DOI: 10.1016/j.cma.2023.116496. URL: https://www.sciencedirect.com/science/article/pii/S0045782523006205.

Carbonell, J.M. et al. (June 2022). "Geotechnical Particle Finite Element Method for Modeling of Soil-Structure Interaction under Large Deformation Conditions". In: *Journal of Rock Mechanics and Geotechnical Engineering* 14.3, pp. 967-983. ISSN: 16747755. DOI: 10.1016/j.jrmge.2021.12.006. URL: https://linkinghub.elsevier.com/retrieve/pii/S1674775521001943.

Charpin, L. and A. Ehrlacher (1st Sept. 2014). "Estimating the Poroelastic Properties of Cracked Materials". In: *Acta Mechanica* 225.9, pp. 2501-2519. ISSN: 1619-6937. DOI: 10.1007/s00707-013-1082-0. URL: https://doi.org/10.1007/s00707-013-1082-0.

Coetzee, C. and O.C. Scheffler (July 2023). "Comparing Particle Shape Representations and Contact Models for DEM Simulation of Bulk Cohesive Behaviour". In: *Computers and Geotechnics* 159, p. 105449. ISSN: 0266352X. DOI: 10.1016/j.compgeo.2023.105449. URL: https://linkinghub.elsevier.com/retrieve/pii/S0266352X23002069.

Collins-Craft, N.A. and V. Acary (22nd Aug. 2025). *On the Formulation and Implementation of Mixed Mode I and Mode II Extrinsic Cohesive Zone Models with Contact and Friction*. URL: https://inria.hal.science/hal-04447397. Pre-published.

Collins-Craft, N.A., F. Bourrier and V. Acary (25th Nov. 2021). "An Extrinsic Cohesive Zone Model with Contact Developed in the Non-Smooth Mechanics Framework". In: 33rd Nordic Seminar on Computational Mechanics. Jönköping.

— (1st Oct. 2022). "On the Formulation and Implementation of Extrinsic Cohesive Zone Models with Contact". In: *Computer Methods in Applied Mechanics and Engineering* 400, p. 115545. ISSN: 0045-7825. DOI: 10.1016/j.cma.2022.115545. URL: https://www.sciencedirect.com/science/article/pii/S0045782522005369.

Cottle, R.W., J.-S. Pang and R.E. Stone (2009). *The Linear Complementarity Problem*. Ed. by R.E. O'Malley. Philadelphia: Society for Industrial and Applied Mathematics. 781 pp. ISBN: 978-0-89871-686-3. DOI: 10.1137/1.9780898719000. URL: http://epubs.siam.org/doi/book/10.1137/1.9780898719000.

Coussy, O. (2004). *Poromechanics*. 2nd ed. Chichester, England; Hoboken, NJ: Wiley. 298 pp. ISBN: 978-0-470-84920-0.

Cundall, P.A. and O.D.L. Strack (Mar. 1979). "A discrete numerical model for granular assemblies". In: *Géotechnique* 29.1, pp. 47-65. ISSN: 0016-8505. DOI: 10.1680/geot.1979.29.1.47. eprint: https://www.emerald.com/jgeot/article-pdf/29/1/47/2514611/geot_1979_29_1_47.pdf. URL: https://doi.org/10.1680/geot.1979.29.1.47.

De Simone, S. et al. (1st Dec. 2023). "Equivalent Biot and Skempton Poroelastic Coefficients for a Fractured Rock Mass from a DFN Approach". In: *Rock Mechanics and Rock Engineering* 56.12, pp. 8907-8925. ISSN: 1434-453X. DOI: 10.1007/s00603-023-03515-9. URL: https://doi.org/10.1007/s00603-023-03515-9.
Dehghani, H. and A. Zilian (1st Sept. 2020). "Poroelastic Model Parameter Identification Using Artificial Neural Networks: On the Effects of Heterogeneous Porosity and Solid Matrix Poisson Ratio". In: *Computational Mechanics* 66.3, pp. 625-649. ISSN: 1432-0924. DOI: 10.1007/s00466-020-01868-4. URL: https://doi.org/10.1007/s00466-020-01868-4.

Dubois, F., V. Acary and M. Jean (Mar. 2018). "The Contact Dynamics Method: A Nonsmooth Story". In: *Comptes Rendus Mécanique* 346.3, pp. 247-262. ISSN: 16310721. DOI: 10.1016/j.crme.2017.12.009. URL: https://linkinghub.elsevier.com/retrieve/pii/S163107211730236X.

Elishakoff, I., Y.J. Ren and M. Shinozuka (May 1997). "New Formulation of FEM for Deterministic and Stochastic Beams through Generalization of Fuchs' Approach". In: *Computer Methods in Applied Mechanics and Engineering* 144.3-4, pp. 235-243. ISSN: 00457825. DOI: 10.1016/S0045-7825(96)01173-5. URL: https://linkinghub.elsevier.com/retrieve/pii/S0045782596011735.

Eshelby, J.D. (20th Aug. 1957). "The Determination of the Elastic Field of an Ellipsoidal Inclusion, and Related Problems". In: *Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences* 241.1226, pp. 376-396. DOI: 10.1098/rspa.1957.0133. URL: https://royalsocietypublishing.org/doi/abs/10.1098/rspa.1957.0133.

Feldfogel, S. et al. (2024). "A Discretization-Convergent Level-Set-Discrete-Element-Method Using a Continuum-Based Contact Formulation". In: *International Journal for Numerical Methods in Engineering* 125.5, e7400. ISSN: 1097-0207. DOI: 10.1002/nme.7400. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.7400.

Geers, M.G.D., V.G. Kouznetsova and W.A.M. Brekelmans (Aug. 2010). "Multi-Scale Computational Homogenization: Trends and Challenges". In: *Journal of Computational and Applied Mathematics* 234.7, pp. 2175-2182. ISSN: 03770427. DOI: 10.1016/j.cam.2009.08.077. URL: https://linkinghub.elsevier.com/retrieve/pii/S0377042709005536.

Ghosh, S., K. Lee and P. Raghavan (Apr. 2001). "A Multi-Level Computational Model for Multi-Scale Damage Analysis in Composite and Porous Materials". In: *International Journal of Solids and Structures* 38.14, pp. 2335-2385. ISSN: 00207683. DOI: 10.1016/S0020-7683(00)00167-0. URL: https://linkinghub.elsevier.com/retrieve/pii/S0020768300001670.

Hashin, Z. (1985). "Large Isotropic Elastic Deformation of Composites and Porous Media". In: *International Journal of Solids and Structures* 21.7, pp. 711-720. ISSN: 00207683. DOI: 10.1016/0020-7683(85)90074-5. URL: https://linkinghub.elsevier.com/retrieve/pii/0020768385900745.

Hill, R. (Sept. 1963). "Elastic Properties of Reinforced Solids: Some Theoretical Principles". In: *Journal of the Mechanics and Physics of Solids* 11.5, pp. 357-372. ISSN: 00225096. DOI: 10.1016/0022-5096(63)90036-X. URL: https://linkinghub.elsevier.com/retrieve/pii/002250966390036X.

— (Aug. 1965). "A Self-Consistent Mechanics of Composite Materials". In: *Journal of the Mechanics and Physics of Solids* 13.4, pp. 213-222. ISSN: 00225096. DOI: 10.1016/0022-5096(65)90010-4. URL: https://linkinghub.elsevier.com/retrieve/pii/0022509665900104.

Huang, D. et al. (Dec. 2021). "Boundary Element Method for the Elastic Contact Problem with Hydrostatic Load at the Contact Interface". In: *Applied Surface Science Advances* 6, p. 100176. ISSN: 26665239. DOI: 10.1016/j.apsadv.2021.100176. URL: https://linkinghub.elsevier.com/retrieve/pii/S2666523921001227.

Jean, M. and J.J. Moreau (1987). "Dynamics in the Presence of Unilateral Contacts and Dry Friction: A Numerical Approach". In: *Unilateral Problems in Structural Analysis* — 2. Springer Vienna, pp. 151-196. DOI: 10.1007/978-3-7091-2967-8_10.

Jourdan, F., P. Alart and M. Jean (16th Mar. 1998). "A Gauss-Seidel like Algorithm to Solve Frictional Contact Problems". In: *Computer Methods in Applied Mechanics and Engineering* 155.1, pp. 31-47. ISSN: 0045-7825. DOI: 10.1016/S0045-7825(97)00137-0. URL: https://www.sciencedirect.com/science/article/pii/S0045782597001370.

Kachanov, M., I. Tsukrov and B. Shafiro (1st Jan. 1994). "Effective Moduli of Solids With Cavities of Various Shapes". In: *Applied Mechanics Reviews* 47 (1S), S151-S174. ISSN: 0003-6900, 2379-0407. DOI: 10.1115/1.3122810. URL: https://asmedigitalcollection.asme.org/appliedmechanicsreviews/article/47/1S/S151/401239/Effective-Moduli-of-Solids-With-Cavities-of.

Kanit, T. et al. (June 2003). "Determination of the Size of the Representative Volume Element for Random Composites: Statistical and Numerical Approach". In: *International Journal of Solids and Structures* 40.13-14, pp. 3647-3679. ISSN: 00207683. DOI: 10.1016/S0020-7683(03)00143-4. URL: https://linkinghub.elsevier.com/retrieve/pii/S0020768303001434.

Kanouté, P. et al. (1st Mar. 2009). "Multiscale Methods for Composites: A Review". In: *Archives of Computational Methods in Engineering* 16.1, pp. 31-75. ISSN: 1886-1784. DOI: 10.1007/s11831-008-9028-8. URL: https://doi.org/10.1007/s11831-008-9028-8.

Kawamoto, R. et al. (June 2016). "Level Set Discrete Element Method for Three-Dimensional Computations with Triaxial Case Study". In: *Journal of the Mechanics and Physics of Solids* 91, pp. 1-13. ISSN: 00225096. DOI: 10.1016/j.jmps.2016.02.021. URL: https://linkinghub.elsevier.com/retrieve/pii/S002250961530154X.

Kouznetsova, V.G., W.A.M. Brekelmans and F.P.T. Baaijens (1st Jan. 2001). "An Approach to Micro-Macro Modeling of Heterogeneous Materials". In: *Computational Mechanics* 27.1, pp. 37-48. ISSN: 1432-0924. DOI: 10.1007/s004660000212. URL: https://doi.org/10.1007/s004660000212.
Larsson, F. et al. (Jan. 2011). "Computational Homogenization Based on a Weak Format of Micro-Periodicity for RVE-problems". In: *Computer Methods in Applied Mechanics and Engineering* 200.1-4, pp. 11-26. ISSN: 00457825. DOI: 10.1016/j.cma.2010.06.023. URL: https://linkinghub.elsevier.com/retrieve/pii/S0045782510001908.

Lee, E.S. et al. (Apr. 2019). "Experimental Investigation of Macroscopic Material Nonlinear Behavior and Microscopic Void Volume Fraction Change for Porous Materials under Uniaxial Compression". In: *Composites Part B: Engineering* 163, pp. 130-138. ISSN: 13598368. DOI: 10.1016/j.compositesb.2018.11.019. URL: https://linkinghub.elsevier.com/retrieve/pii/S1359836818311922.

Lei, Z. et al. (Dec. 2022). "A Material Point-Finite Element (MPM-FEM) Model for Simulating Three-Dimensional Soil-Structure Interactions with the Hybrid Contact Method". In: *Computers and Geotechnics* 152, p. 105009. ISSN: 0266352X. DOI: 10.1016/j.compgeo.2022.105009. URL: https://linkinghub.elsevier.com/retrieve/pii/S0266352X22003469.

Lemke, C.E. (May 1965). "Bimatrix Equilibrium Points and Mathematical Programming". In: *Management Science* 11.7, pp. 681-689. ISSN: 0025-1909. DOI: 10.1287/mnsc.11.7.681. URL: http://pubsonline.informs.org/doi/abs/10.1287/mnsc.11.7.681.

Lemke, C.E. and J.T. Howson Jr. (June 1964). "Equilibrium Points of Bimatrix Games". In: *Journal of the Society for Industrial and Applied Mathematics* 12.2, pp. 413-423. ISSN: 0368-4245. DOI: 10.1137/0112033. URL: http://epubs.siam.org/doi/10.1137/0112033.

Lopez-Pamies, O. (20th Oct. 2006). "On the Effective Behavior, Microstructure Evolution and Macroscopic Stability of Elastomeric Composites." PhD thesis, Philadelphia & Saclay: University of Pennsylvania & École Polytechnique. 219 pp. URL: https://pastel.hal.science/pastel-00002978v1/document.

Lopez-Pamies, O. and P. Ponte Castañeda (1st Sept. 2004). "Second-Order Estimates for the Macroscopic Response and Loss of Ellipticity in Porous Rubbers at Large Deformations". In: *Journal of Elasticity* 76.3, pp. 247-287. ISSN: 1573-2681. DOI: 10.1007/s10659-005-1405-z. URL: https://doi.org/10.1007/s10659-005-1405-z.

— (Aug. 2007a). "Homogenization-Based Constitutive Models for Porous Elastomers and Implications for Macroscopic Instabilities: I—Analysis". In: *Journal of the Mechanics and Physics of Solids* 55.8, pp. 1677-1701. ISSN: 00225096. DOI: 10.1016/j.jmps.2007.01.007. URL: https://linkinghub.elsevier.com/retrieve/pii/S0022509607000270.

— (Aug. 2007b). "Homogenization-Based Constitutive Models for Porous Elastomers and Implications for Macroscopic Instabilities: II—Results". In: *Journal of the Mechanics and Physics of Solids* 55.8, pp. 1702-1728. ISSN: 00225096. DOI: 10.1016/j.jmps.2007.01.008. URL: https://linkinghub.elsevier.com/retrieve/pii/S0022509607000269.

Lydzba, D. and J.F. Shao (2000). "Study of Poroelasticity Material Coefficients as Response of Microstructure". In: *Mechanics of Cohesive-frictional Materials* 5.2, pp. 149-171. ISSN: 1099-1484. DOI: 10.1002/(SICI)1099-1484(200002)5:2<149::AID-CFM89>3.0.CO;2-F. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/%28SICI%291099-1484%28200002%295%3A2%3C149%3A%3AAID-CFM89%3E3.0.CO%3B2-F.

Machado, M. et al. (July 2012). "Compliant Contact Force Models in Multibody Dynamics: Evolution of the Hertz Contact Theory". In: *Mechanism and Machine Theory* 53, pp. 99-121. ISSN: 0094-114X. DOI: 10.1016/j.mechmachtheory.2012.02.010. URL: https://linkinghub.elsevier.com/retrieve/pii/S0094114X1200047X.

Mandel, J. (1971). *Plasticite classique et viscoplasticite*. Vol. 565. CISM International Centre for Mechanical Sciences. Cham: Springer International Publishing.

Miehe, C. (31st Jan. 2003). "Computational Micro-to-Macro Transitions for Discretized Micro-Structures of Heterogeneous Materials at Finite Strains Based on the Minimization of Averaged Incremental Energy q". In: *Computer Methods in Applied Mechanics and Engineering* 192.5-6, pp. 559-591. DOI: 10.1016/S0045-7825(02)00564-9.

Miehe, C. and A. Koch (1st July 2002). "Computational Micro-to-Macro Transitions of Discretized Microstructures Undergoing Small Strains". In: *Archive of Applied Mechanics* 72.4, pp. 300-317. ISSN: 1432-0681. DOI: 10.1007/s00419-002-0212-2. URL: https://doi.org/10.1007/s00419-002-0212-2.

Miehe, C., J. Schröder and M. Becker (Oct. 2002). "Computational Homogenization Analysis in Finite Elasticity: Material and Structural Instabilities on the Micro- and Macro-Scales of Periodic Composites and Their Interaction". In: *Computer Methods in Applied Mechanics and Engineering* 191.44, pp. 4971-5005. ISSN: 00457825. DOI: 10.1016/S0045-7825(02)00391-2. URL: https://linkinghub.elsevier.com/retrieve/pii/S0045782502003912.

Monerie, Y. and V. Acary (8th Jan. 2001). "Formulation Dynamique d'un Modèle de Zone Cohésive Tridimensionnel Couplant Endommagement et Frottement". In: *Revue Européenne des Éléments Finis* 10.2-4, pp. 489-503. ISSN: 1250-6559. DOI: 10.1080/12506559.2001.11869264. URL: https://www.tandfonline.com/doi/full/10.1080/12506559.2001.11869264.

Moraleda, J., J. Segurado and J. Llorca (11th Dec. 2007). "Finite Deformation of Porous Elastomers: A Computational Micromechanics Approach". In: *Philosophical Magazine* 87.35, pp. 5607-5627. ISSN: 1478-6435, 1478-6443. DOI: 10.1080/14786430701678930. URL: http://www.tandfonline.com/doi/abs/10.1080/14786430701678930.

Moreau, J.J. (1970). "Sur Les Lois de Frottement, de Plasticité et de Viscosité". In: *Comptes rendus de l'Académie des sciences. Série A - Sciences mathématiques* 271, pp. 608-611. URL: https://hal.science/hal-01868140/document.

— (12th July 1971). "Sur l'évolution d'un système élasto-visco-plastique". In: *Comptes rendus de l'Académie des sciences. Série A - Sciences mathématiques* 273, pp. 118-121.
Moreau, J.J. (1986). "Une Formulation Du Contact à Frottement Sec; Application Au Calcul Numérique". In: *Comptes rendus de l'Académie des sciences. Série 2, Mécanique, Physique, Chimie, Sciences de l'univers, Sciences de la Terre* 302.13, pp. 799-801. ISSN: 0764-4450. URL: https://gallica.bnf.fr/ark:/12148/bpt6k5675350v/f805.item.

Moreau, J.J. and P.D. Panagiotopoulos, eds. (1988). *Nonsmooth Mechanics and Applications*. Vienna: Springer Vienna. ISBN: 978-3-211-82066-7. DOI: 10.1007/978-3-7091-2624-0. URL: http://link.springer.com/10.1007/978-3-7091-2624-0.

Mori, T. and K. Tanaka (1st May 1973). "Average Stress in Matrix and Average Elastic Energy of Materials with Misfitting Inclusions". In: *Acta Metallurgica* 21.5, pp. 571-574. ISSN: 0001-6160. DOI: 10.1016/0001-6160(73)90064-3. URL: https://www.sciencedirect.com/science/article/pii/0001616073900643.

Nemat-Nasser, S. and M. Hori (1998). *Micromechanics: Overall Properties of Heterogeneous Materials*. Second. Vol. 37. North-Holland Series in Applied Mathematics and Mechanics. Amsterdam: North-Holland.

Nguyen, V.-D. et al. (Apr. 2012). "Imposing Periodic Boundary Condition on Arbitrary Meshes by Polynomial Interpolation". In: *Computational Materials Science* 55, pp. 390-406. ISSN: 09270256. DOI: 10.1016/j.commatsci.2011.10.017. URL: https://linkinghub.elsevier.com/retrieve/pii/S0927025611005866.

Papadopoulos, P. and J.M. Solberg (Aug. 1998). "A Lagrange Multiplier Method for the Finite Element Solution of Frictionless Contact Problems". In: *Mathematical and Computer Modelling* 28.4-8, pp. 373-384. ISSN: 0895-7177. DOI: 10.1016/s0895-7177(98)00128-9. URL: https://linkinghub.elsevier.com/retrieve/pii/S0895717798001289.

Perić, D. and D.R.J. Owen (1992). "Computational Model for 3-D Contact Problems with Friction Based on the Penalty Method". In: *International Journal for Numerical Methods in Engineering* 35.6, pp. 1289-1309. ISSN: 1097-0207. DOI: 10.1002/nme.1620350609. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.1620350609.

Ponte Castañeda, P. and M. Zaidman (Sept. 1994). "Constitutive Models for Porous Materials with Evolving Microstructure". In: *Journal of the Mechanics and Physics of Solids* 42.9, pp. 1459-1497. ISSN: 00225096. DOI: 10.1016/0022-5096(94)90005-1. URL: https://linkinghub.elsevier.com/retrieve/pii/0022509694900051.

Pore, T., S.G. Thorat and A.A. Nema (2021). "Review of Contact Modelling in Nonlinear Finite Element Analysis". In: *Materials Today: Proceedings* 47, pp. 2436-2440. ISSN: 2214-7853. DOI: 10.1016/j.matpr.2021.04.504. URL: https://linkinghub.elsevier.com/retrieve/pii/S2214785321034295.

Rohan, E. and J. Heczko (Apr. 2020). "Homogenization and Numerical Modelling of Poroelastic Materials with Self-Contact in the Microstructure". In: *Computers & Structures* 230, p. 106086. ISSN: 0045-7949. DOI: 10.1016/j.compstruc.2019.06.003. URL: https://linkinghub.elsevier.com/retrieve/pii/S0045794919300264.

Saeb, S., P. Steinmann and A. Javili (1st Sept. 2016). "Aspects of Computational Homogenization at Finite Deformations: A Unifying Review From Reuss' to Voigt's Bound". In: *Applied Mechanics Reviews* 68.5, p. 050801. ISSN: 0003-6900, 2379-0407. DOI: 10.1115/1.4034024. URL: https://asmedigitalcollection.asme.org/appliedmechanicsreviews/article/doi/10.1115/1.4034024/443653/Aspects-of-Computational-Homogenization-at-Finite.

Sahimi, M. (1st Oct. 1993). "Flow Phenomena in Rocks: From Continuum Models to Fractals, Percolation, Cellular Automata, and Simulated Annealing". In: *Reviews of Modern Physics* 65.4, pp. 1393-1534. ISSN: 0034-6861, 1539-0756. DOI: 10.1103/RevModPhys.65.1393. URL: https://link.aps.org/doi/10.1103/RevModPhys.65.1393.

Schenk, O. and K. Gärtner (Apr. 2004). "Solving Unsymmetric Sparse Systems of Linear Equations with PARDISO". In: *Future Generation Computer Systems* 20.3, pp. 475-487. ISSN: 0167739X. DOI: 10.1016/j.future.2003.07.011. URL: https://linkinghub.elsevier.com/retrieve/pii/S0167739X03001882.

Selvadurai, A.P.S. and A.P. Suvorov (3rd Nov. 2020). "The Influence of the Pore Shape on the Bulk Modulus and the Biot Coefficient of Fluid-Saturated Porous Rocks". In: *Scientific Reports* 10.1, p. 18959. ISSN: 2045-2322. DOI: 10.1038/s41598-020-75979-6. URL: https://www.nature.com/articles/s41598-020-75979-6.

Smeets, S., N. Renaud and L.J.C. Van Willenswaard (9th Oct. 2022). "Nanomesh: A Python Workflow Tool for Generating Meshes from Image Data". In: *Journal of Open Source Software* 7.78, p. 4654. ISSN: 2475-9066. DOI: 10.21105/joss.04654. URL: https://joss.theoj.org/papers/10.21105/joss.04654.

Stewart, D.E. and J.C. Trinkle (15th Aug. 1996). "An Implicit Time-Stepping Scheme for Rigid Body Dynamics with Inelastic Collisions and Coulomb Friction". In: *International Journal for Numerical Methods in Engineering* 39.15, pp. 2673-2691. ISSN: 0029-5981. DOI: 10.1002/(SICI)1097-0207(19960815)39:15<2673::AID-NME972>3.0.CO;2-I. URL: https://onlinelibrary.wiley.com/doi/10.1002/(SICI)1097-0207(19960815)39:15%3C2673::AID-NME972%3E3.0.CO;2-I.

Tan, X. and H. Konietzky (Jan. 2014). "Numerical Study of Variation in Biot's Coefficient with Respect to Microstructure of Rocks". In: *Tectonophysics* 610, pp. 159-171. ISSN: 00401951. DOI: 10.1016/j.tecto.2013.11.014. URL: https://linkinghub.elsevier.com/retrieve/pii/S0040195113006616.

Tang, Y. et al. (9th Mar. 2020). "Experiment Research on Physical Clogging Mechanism in the Porous Media and Its Impact on Permeability". In: *Granular Matter* 22.2, p. 37. ISSN: 1434-7636. DOI: 10.1007/s10035-020-1001-8. URL: https://doi.org/10.1007/s10035-020-1001-8.

Temizer, İ. and P. Wriggers (Dec. 2008). "On the Computation of the Macroscopic Tangent for Multiscale Volumetric Homogenization Problems". In: *Computer Methods in Applied Mechanics and Engineering* 198.3-4, pp. 495-510. ISSN: 00457825. DOI: 10.1016/j.cma.2008.08.018. URL: https://linkinghub.elsevier.com/retrieve/pii/S0045782508003149.
Terada, K. et al. (Apr. 2000). "Simulation of the Multi-Scale Convergence in Computational Homogenization Approaches". In: *International Journal of Solids and Structures* 37.16, pp. 2285-2311. ISSN: 00207683. DOI: 10.1016/S0020-7683(98)00341-2. URL: https://linkinghub.elsevier.com/retrieve/pii/S0020768398003412.

Thorpe, M.F. and P.N. Sen (May 1985). "Elastic moduli of two-dimensional composite continua with elliptical inclusions". In: *The Journal of the Acoustical Society of America* 77.5, pp. 1674-1680. ISSN: 0001-4966. DOI: 10.1121/1.391966. URL: https://doi.org/10.1121/1.391966.

Tsukrov, I. and M. Kachanov (11th Oct. 2000). "Effective Moduli of an Anisotropic Material with Elliptical Holes of Arbitrary Orientational Distribution". In: *International Journal of Solids and Structures* 37.41, pp. 5919-5941. DOI: 10.1016/S0020-7683(99)00244-9.

Van Dijk, N.P. (2016). "Formulation and Implementation of Stress-Driven and/or Strain-Driven Computational Homogenization for Finite Strain". In: *International Journal for Numerical Methods in Engineering* 107.12, pp. 1009-1028. ISSN: 1097-0207. DOI: 10.1002/nme.5198. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.5198.

Wu, H., N. Guo and J. Zhao (1st June 2018). "Multiscale Modeling and Analysis of Compaction Bands in High-Porosity Sandstones". In: *Acta Geotechnica* 13.3, pp. 575-599. ISSN: 1861-1133. DOI: 10.1007/s11440-017-0560-2. URL: https://doi.org/10.1007/s11440-017-0560-2.

Yastrebov, Vladislav A. (2013). *Numerical Methods in Contact Mechanics*. Numerical Methods in Engineering Series. London : Hoboken, NJ: ISTE ; Wiley. 391 pp. ISBN: 978-1-84821-519-1.

Zhang, Q. et al. (May 2023). "A Coupling Approach of the Isogeometric-Meshfree Method and Peridynamics for Static and Dynamic Crack Propagation". In: *Computer Methods in Applied Mechanics and Engineering* 410, p. 115904. ISSN: 00457825. DOI: 10.1016/j.cma.2023.115904. URL: https://linkinghub.elsevier.com/retrieve/pii/S0045782523000270.

Zhou, X.-Y. et al. (Mar. 2016). "Perturbation-Based Stochastic Multi-Scale Computational Homogenization Method for the Determination of the Effective Properties of Composite Materials with Random Properties". In: *Computer Methods in Applied Mechanics and Engineering* 300, pp. 84-105. ISSN: 00457825. DOI: 10.1016/j.cma.2015.10.020. URL: https://linkinghub.elsevier.com/retrieve/pii/S0045782515003497.