J.-L Adia

J Yvonnet
email: julien.yvonnet@univ-paris-est.fr

Q.-C He

N C Tran

J Sanahuja

A combined Lattice-Boltzmann-Finite Element approach to modeling unsaturated poroelastic behavior of heterogeneous media

Keywords: Porous media, capillary films, concrete, shrinkage, drying, Lattice-Boltzmann method, Unsaturated poroelasticity

In this paper, we propose a procedure combining Lattice-Boltzmann and finite element simulations to model the effects of capillary pressure in porous microstructures. Starting from an explicit geometry of the microstructure, the Lattice-Boltzmann method is used to simulate the condensation from vapor phase to liquid and predict the geometry of capillary liquid films and liquid phases for arbitrary geometry of the solid skeleton. Then, a procedure is provided to prescribe surface tension and pressure due to the capillary liquid films and liquid phases. The solid skeleton can deform elastically under the action of the liquid and its deformation can be captured using the finite element method. Finally, a procedure to extract the shrinkage strain and the effective macroscopic behavior of the material, taking into account the explicit geometry of the microstructure and of the liquid capillary phases for a given saturation is described. Several numerical examples are provided to validate the methodology and show its potential.

Introduction

Shrinkage can be defined as the change of volume in porous materials (geomaterials, cementitious materials, wood, etc.), which are subjected to a drying environment. For example, in the context of nuclear power plant installations, the last barrier against accidental release of radioactive fission products into the environment is a pre-stressed concrete containment building. Such structure cannot be replaced and the issue of aging management must be treated with particular care for long-term operations [START_REF] Sanahuja | Material properties prediction for long term operation of nuclear power plants civil engineering structures: challenges at edf[END_REF]. During the life of the concrete containment vessel, it is subjected to a drying environment. Then, the concrete loaded by the pre-stressed tendons and subjected to drying exhibits shrinkage and creep deformations. Due to these deformations, the pre-stress in the concrete containment vessel which is designed to withstand various accidental situations decreases over time. Therefore, evaluating the lifetime of such structures requires a fine understanding of the mechanisms of shrinkage and creep. It is now recognized that creep and shrinkage are related phenomena and depend both on adsorption and desorption of water [START_REF] Bažant | Creep and hygrothermal effects in concrete structures[END_REF]. Understanding shrinkage is prerequisite to understanding creep. The origin of shrinkage has been primarily attributed to the three most prominent mechanisms: surface free energy [START_REF] Powers | Mechanisms of shrinkage and reversible creep of hardened portland cement paste[END_REF][START_REF] Wittmann | Interaction of hardened cement paste and water[END_REF], capillary tension [START_REF] Powers | Mechanisms of shrinkage and reversible creep of hardened portland cement paste[END_REF], and disjoining pressure [START_REF] Bažant | Thermodynamics of interacting continua with surfaces and creep analysis of concrete structures[END_REF][START_REF] Wittmann | Interaction of hardened cement paste and water[END_REF].

The principle of defining an effective stress in saturated soil has been initiated by [6]. Then in [START_REF] Biot | General theory of three-dimensional consolidation[END_REF][START_REF] Biot | The elastic coeff cients of the theory of consolidation[END_REF], Biot proposed a framework for isotropic saturated poroelasticity. These works have been extended to unsaturated porous media in [START_REF] Biot | The use of pore water coefficients in practice[END_REF][START_REF] Biot | The principle of effective stress[END_REF]. A generalized thermodynamic approach for the saturated and unsaturated porous media can be found in [START_REF] Coussy | The equivalent pore pressure and the swelling and shrinkage of cement-based materials[END_REF][START_REF] Coussy | Poromechanics[END_REF][START_REF] Tabbal | Modelling the drying shrinkage of porous materials by considering both capillary and adsorption effects[END_REF]. In the unsaturated context, a thermo-hydro-mechanical coupling associated with surface tension, or the energy related to each fluid-fluid or fluid-solid interface is required. It induces macroscopic effects such as drying shrinkage. Only recently, attempts to take into account the effects of the microstructure on the unsaturated behavior have been proposed in [START_REF] Chateau | Micromechanics of saturated and unsaturated porous media[END_REF][START_REF] Dormieux | Microporomechanics[END_REF].

To improve understanding of shrinkage mechanisms, the present paper proposes a numerical homogenization method to identify the macroscopic constitutive law of an unsaturated poroelastic medium to extend the classical analytical approach [START_REF] Dormieux | Microporomechanics[END_REF] to arbitrary microstructures. A hybrid method combining Lattice-Boltzmann method (LBM) and Finite Element Method (FEM) is proposed here for the first time to achieve the purpose.

The Lattice Boltzmann Method (LBM) (see e.g. [START_REF] Perumal | A review on the development of lattice boltzmann computation of macro fluid flows and heat transfer[END_REF][START_REF] Huang | Multiphase Lattice Boltzmann Methods: Theory and Application[END_REF]) has been introduced as an efficient numerical tool for simulating fluid problems. Owing to its ease of implementation, its kinetic theory basis and its constitutive versatility, LBM plays an essential role as a simulation tool for understanding micro and macro fluid flows. In addition, LBM can be efficiently solved on parallel comput-ers. LBM has been used for macroscopic transport properties [START_REF] Gunstensen | A Lattice-Boltzmann studies of immiscible two-phase flow through porous media[END_REF][START_REF] Martys | Simulation of multicomponent fluids in complex three-dimensional geometries by the Lattice Boltzmann method[END_REF], relative permeability [START_REF] Zhao | Relative permeability of two immiscible fluids flowing through porous media determined by Lattice Boltzmann method[END_REF], liquid/vapor transitions [START_REF] Montessori | Lattice Boltzmann modeling of complex flows for engineering applications[END_REF], ionic diffusivity, nano-and micro-scale simulations in the presence of complex porous material, heterogeneous catalysis [START_REF] Falcucci | Heterogeneous catalysis in pulsed-flow reactors with nanoporous gold hollow spheres[END_REF][START_REF] Montessori | Extended friction elucidates the breakdown of fast water transport in graphene oxide membranes[END_REF][START_REF] Vesselin | Simulating engineering flows through complex porous media via the lattice boltzmann method[END_REF][START_REF] Matthew M Montemore | Effect of nanoscale flows on the surface structure of nanoporous catalysts[END_REF], or for permittivity calculations [START_REF] Zhang | Modeling of ionic diffusivity in nonsaturated cement-based materials using Lattice Boltzmann method[END_REF][START_REF] Guihard | Assessing the permittivity of an unsaturated sand by combining a Lattice Boltzmann method simulation, electromagnetic homogenization models and measurements[END_REF]. For unsaturated porous media, LBM was used to determine the capillary curve of porous microstructures in [START_REF] Pan | Lattice-Boltzmann simulation of two-phase flow in porous media[END_REF][START_REF] Sukop | Lattice Boltzmann method for modeling liquidvapor interface configurations in porous media[END_REF] .

The LBM is here used to describe the distribution of the capillary liquid water within porous microstructure, while the FEM is employed to compute the elastic deformation of the solid skeleton under the action of the fluid phase. Previous studies combining the Lattice Boltzmann method and the finite element method can be found e.g. in [START_REF] Macmeccan | Simulating deformable particle suspensions using a coupled Lattice-Boltzmann and finiteelement method[END_REF][START_REF] De Rosis | Hydroelastic analysis of hull slamming coupling lattice boltzmann and finite element methods[END_REF]. The present procedure aims at modeling the effects of capillary pressure in porous media at the microscale. The LBM is used to simulate the condensation from vapor to liquid and predict the geometry of capillary liquid films and liquid phases for arbitrary geometry of the microstructure. Then, a procedure is provided to prescribe surface tension and pressure due to the capillary liquid films and liquid phases. The FEM is used to determine the resulting elastic deformation of the solid skeleton. We elaborate a procedure to numerically extract the shrinkage strain and the effective macroscopic behavior of the material.

In section 2, the overall LBM-FEM methodology is presented. In section 3, the Lattice Boltzmann Method (LBM) is reviewed and the relevant algorithm implemented in this paper is described. In section 4, we present the continuum model used at the microscale, which is constructed from explicit morphologies and explicit description of liquid phases forming capillary films, obtained from the LBM simulations. The homogenized model is detailed in section 5, together with the different computational methodologies, to extract the effective quantities of the resulting poroelastic medium from the RVE model. Finally, applications are presented in section 6 to validate the present model and to apply it to 3D porous microstructures to show the potential of this approach.

Overall methodology

The objective of this work is to determine the effective behavior of an unsaturated porous medium, while explicitly taking into account the formation and pressure of capillary water on the solid skeleton, which is assumed to be a linearly elastic phase. The different steps of the methodology are described in Fig. 1. First, given an arbitrary geometry of the dry porous skeleton (Fig. 1 the fluid distribution by liquid, vapor phases and thin liquid films by imperfect interfaces modeling the surface tension at solid-vapor and liquid/vapor interfaces (Fig. 1 (c)). Given this model, the pressure distribution on the solid skeleton is extracted and a homogenized poroelastic model as well as shrinkage strains can be defined and computed with finite elements. The different steps of the procedure are described in the following.

3 Lattice-Boltzmann Method (LBM) for prediction of capillary liquid distribution from condensation of vapor

A brief review of LBM for flows with phase change

In this section we briefly review the basics of LBM and its extension to flows with phase change. In the present work, we propose to describe explicitly the porous microstructure geometry to estimate the macroscopic isothermal and unsaturated poroelastic properties of a porous material for a given saturation ratio. It is then essential to accurately determine the distribution of the capillary liquid water within the porous network, and to obtain a realistic repartition of liquid and vapor in the RVE. This can be performed by multiphase fluid flow simulations. On one hand, macroscopic descriptions (Computational Fluid Dynamics, CFD) of multiphase flows require specific treatments to track the interfaces [START_REF] Tryggvason | A front-tracking method for the computations of multiphase flow[END_REF] and imply tedious computational implementation. On the other hand, microscopic modeling based on Molecular Dynamics [START_REF] Binder | Molecular dynamics simulations[END_REF] is lim-ited by small volumes of description. In this work, we adopt a mesoscopic approach to the multiphase flow using LBM [START_REF] Shan | Kinetic theory representation of hydrodynamics: a way beyond the Navier-Stokes equation[END_REF]. LBM is a discretized form of the Boltzmann equation [START_REF] Cercignani | The Boltzmann equation and its applications[END_REF], which models the evolution of probability to find a fluid molecule at a given position with a given velocity. LBM has advantages as compared to the macro and micro approaches mentioned above. First, due to its kinetic theory basis, it has a more physical foundation than CFD in the case of multiphase flows. Second, it does not require tracking the interfaces between liquid and vapor phases. Liquid and vapor phases emerge from the physical nature of the interactions between fluid molecules at the microscopic level, which are governed by a mean-field interaction potential [START_REF] Rowlinson | Molecular theory of capillarity[END_REF]. Finally, LBM allows efficient parallel computations.

There are several popular multiphase LBM models, including: the color-gradient (RK) model [START_REF] Gunstensen | Lattice Boltzmann model of immiscible fluids[END_REF], based on the immiscible lattice gas of [START_REF] Rothman | Immiscible cellular-automata: Simple models of complex hydrodynamics[END_REF], the Shan-Chen (SC) model [START_REF] Shan | Lattice Boltzmann model for simulating flows with multiple phases and components[END_REF][START_REF] Shan | Simulation of nonideal gases and liquid-gas phase transitions by the Lattice Boltzmann equation[END_REF] which is based on incorporating an attractive or repulsive force leading to phase separation, the free-energy (FE) model proposed by Swift-Osborn-Yeomans [START_REF] Swift | Lattice Boltzmann simulation of nonideal fluids[END_REF] and the He-Chen-Zhang (HCZ) model [START_REF] He | A Lattice Boltzmann scheme for incompressible multiphase flow and its application in simulation of Rayleigh-Taylor instability[END_REF]. In this paper, we have implemented a modified version of the SC model. The SC model is a good compromise between ease of implementation, computing resource and model accuracy and has been widely used during the past decades to model multiphase flow in porous media [START_REF] Martys | Simulation of multicomponent fluids in complex three-dimensional geometries by the Lattice Boltzmann method[END_REF][START_REF] Pan | Lattice-Boltzmann simulation of two-phase flow in porous media[END_REF][START_REF] Sukop | Lattice Boltzmann method for modeling liquidvapor interface configurations in porous media[END_REF].

To describe the two-phase SC LBM, the BBGKY hierarchy [START_REF] Bogoliubov | Kinetics equations[END_REF] formalism is adopted here and the mechanisms at the origin of the phase transitions and wettability phenomena are introduced following [START_REF] Benzi | Mesoscopic modeling of a two-phase in the presence of boundaries: the contact angle[END_REF]. The phase transitions and wettability in this model are induced by the interaction forces between the fluid and/or the solid molecules. The interested reader can find more details about theoretical foundations of LBM e.g. in [START_REF] Succi | The lattice Boltzmann equation: for complex states of flowing matter[END_REF].

Numerical implementation

The Lattice Boltzmann equations are derived by discretizing the Boltzmann-BGK equation in space and time through a projection of probability distribution function on a Hermite polynomial basis [START_REF] Grad | Note on n-dimensional Hermite polynomials[END_REF][START_REF] Grad | On the kinetic theory of rarefied gases[END_REF][START_REF] Shan | Discretization of the velocity space in the solution of the boltzmann equation[END_REF][START_REF] Shan | Kinetic theory representation of hydrodynamics: a way beyond the Navier-Stokes equation[END_REF] according to a Taylor series technique, i.e:

f k (x + ξ ξ ξ k ∆t, t + ∆t) streaming = f k (x, t) - ∆t τ c f k (x, t) -f eq k ρ(x, t), v (x, t) collision , (1)
where the quadrature formula in the velocity space discretization is used here in 3D, containing 19 quadrature nodes with an algebraic prediction of 5 order (D 3 Q 19 discretization scheme, see Fig. 2).

ξ ξ ξ k = ∆x ∆t        0 1 -1 0 0 0 0 1 -1 1 -1 1 -1 1 -1 0 0 0 0 0 0 0 1 -1 0 0 1 1 -1 -1 0 0 0 0 1 -1 1 -1 0 0 0 0 0 1 -1 0 0 0 0 1 1 -1 -1 1 1 -1 -1        . (2)
The discrete equilibrium distribution function is given by

f eq k (ρ, v ) = ω k ρ 1 + ξ ξ ξ k • v c 2 s + (ξ ξ ξ k • v ) 2 2c 4 s - v 2 2c 2 s , (3)
where ω k is the weight factor associated with the k th -direction given as [START_REF] Gunstensen | A Lattice-Boltzmann studies of immiscible two-phase flow through porous media[END_REF]. The local fluid density and momentum of each lattice node are related to the discrete distribution function by

ω k = 1/3(k = 0), ω k = 1/18(k = 1, • • • , 6) and ω k = 1/36(k = 7, • • • ,
ρ = 18 k=0 f k , v = 1 ρ 18 k=0 f k ξ ξ ξ k , and v = v + τ c f tot ρ . (4)
It can be shown that these quantities evolve according to the Navier-Stokes equations of fluid-dynamics in the incompressible and low Mach number limits, using the Chapmann-Enskog procedure [START_REF] Chapman | The mathematical theory of nonuniform gases[END_REF]:

∂ρ ∂t + ∇ x • (ρv) = 0, (5)
∂(ρv) ∂t + ∇ x • (ρv ⊗ v) = -∇ x • P 0 + f tot + ∇ x • ρη(∇ x v + ∇ T x v) , (6)
where the relaxation time τ c is related to the kinematic viscosity by η = c 2 s (τ -0.5)∆t, τ = ∆t τc . Above, P 0 is the ideal part of the second order fluid pressure tensor given by the perfect gas equation of state, P 0 = c 2 s ρ1. We have also introduced, f tot = f coh + f ads + f g , the total body force term including the self-consistent one, which takes into account the non-ideal effects. More precisely, f coh is responsible of phases separation, f ads for wettability phenomena and f g = ρg for the other possible external forces such as gravitational, electrical forces etc. The discrete forms of the fluid-fluid f coh and fluid-solid f ads interaction forces are given by [START_REF] Shan | Lattice Boltzmann model for simulating flows with multiple phases and components[END_REF][START_REF] Shan | Simulation of nonideal gases and liquid-gas phase transitions by the Lattice Boltzmann equation[END_REF][START_REF] Martys | Simulation of multicomponent fluids in complex three-dimensional geometries by the Lattice Boltzmann method[END_REF]:

f coh (x, t) = -Gψ(x, t) 18 k=1 ω k ψ(x + ξ ξ ξ k ∆t, t)χ f (x + ξ ξ ξ k ∆t)ξ ξ ξ k , (7)
f ads (x, t) = -Gψ(x, t) 18 k=1 ω k ψ(ρ s )χ s (x + ξ ξ ξ k ∆t)ξ ξ ξ k , (8)
where ρ s represents a fictitious density of the solid-fluid interface and χ α is a characteristic function associated to the domain Ω α , α = s, f , such that χ α (x) = 1 if x ∈ Ω α and χ α (x) = 0 otherwise. Above, ψ(x, t) is the lattice version of the mean-field potential. We define the pseudo-potential ψ using the method proposed in [START_REF] Yuan | Equations of state in a Lattice Boltzmann model[END_REF] by incorporating the non-ideal gas equation of state (EOS) [START_REF] Carnahan | Equation of state for nonattracting rigid spheres[END_REF]:

ψ(ρ) = 2 Gc 2 s ∆t p(ρ) -ρc 2 s , (9)
p(ρ) = ρRT 1 + bρ/4 + (bρ/4) 2 -(bρ/4) 3 (1 -bρ/4) 3 , (10)
where a and b are the constants of the Carnahan-Starling equation of state.

The previous EOS incorporation procedure to take into account the physical non-ideal gas EOS reduces some drawbacks of the original SC LBM, such as spurious currents, low liquid-vapor density ratio and introduces explicitly the temperature in our SC LBM. In this approach the phase transition takes place as soon as the temperature T becomes smaller than the critical temperature T crit , where the fluid starts to exhibit two coexisting liquid and vapor phases. At the critical fluid state where there is no distinction between the liquid and vapor phases (ρ l = ρ v = ρ crit ), we can show that a = 0.4963R 2 T 2 crit /p crit and b = 0.18727RT crit /p crit , with p crit the critical pressure.

The non-slip boundary condition at the solid-fluid Γ sf interface is taken into account by the half-way bounce-back scheme (Fig. 3-a),

f k (x f , t + ∆t) = f k (x f , t) - ∆t τ c f k (x f , t) -f eq k (x f , t) , (11)
where ξ ξ ξ k = -ξ ξ ξ k . The periodicity condition on ∂Ω f is ensured (Fig. 3-b) by

f in k (x -, t + ∆t) = f out k (x + , t) - ∆t τ c f out k (x + , t) -f eq,out k (x + , t) (12)
with x -and x + some nodes on the opposite parts of ∂Ω f . The overall algorithm, involving the previously mentioned SC LBM model, is summarized below. The subscript ∞ indicates converged quantities at the thermodynamic equilibrium state .

Algorithm 1: Shan-Chen Lattice Boltzmann algorithm.

Input. ρ 0 (x), v 0 (x), τ c , ∆t, ∆x, g, T /T crit , G

INITIALIZATION. ERR = ERR 0 WHILE. ERR ≤ T OL, given ρ i (x) and v i (x)
(1) Compute body and interaction forces f tot i (x) and shift velocity v i (x).

(2) Collision step:

f i k (x) = f i k (x) -∆t τc f i k (x) -f i,eq k (ρ i (x), v i (x)) . (3) Streaming step: f i+1 k (x + ξ ξ ξ k ∆t) = f i k (x). (4) Prescribe boundary conditions. (5) Update fluid variables ρ i+1 (x) and v i+1 (x), p i+1 (x), σ σ σ i+1 (x). (6) Compute ERR (7) GO to (1). Output. ρ ∞ (x), v ∞ (x), p ∞ (x), σ σ σ ∞ (x).
In the above, the error is defined as:

ERR = i [v i (x, t) -v i (x, t -t f )] 2 i v 2 i (x, t) (13)
where in practice we have used t f = 50.

4 Unsaturated porous medium RVE model with imperfect interfaces

Equations of the model with imperfect interfaces

Let us consider a porous microstructure in unsaturated conditions in a domain denoted by Ω ⊂ R 3 , associated with a Representative Volume Element, as depicted schematically in Fig. 4 (a). For one saturation level, the vapor condensates into liquid phases which distribute within the porous space due to capillary effects (we assume that gravity effects can be safely neglected due to the small scales considered). Along surfaces of the solid skeleton, the liquid phase can form thin capillary liquid films or liquid phases. A simplified model is associated with this configuration, as described in Fig. 4 (b). To avoid meshing the thin films, they are replaced by imperfect interfaces taking into account surface tension due to capillary effects. We then define the solid domain Ω s associated with the solid phase, and denote by Γ s , Γ sv and Γ v the interfaces between liquid and solid phases, solid and vapor phases, and liquid and vapor phases, respectively. The full interface between the solid phase and the pore is denoted by Γ s = Γ s ∪ Γ sv .

For simplicity we adopt ∇ = ∇ x in the sequel. We assume that the solid matrix is linear elastic and consider small perturbations in the solid phase. The equilibrium equations in the solid phase are classically given by

∇ • σ σ σ s = 0 in Ω s , (14)
where ∇ • (.) denotes divergence operator and σ s denotes the Cauchy stress.

The subscript s refers to the solid phase. The constitutive relation is given by the Hooke law,

σ σ σ s = C s : ε ε ε s in Ω s , (15)
where C s denotes the solid phase elastic tensor and ε s the linearized strain tensor. The equilibrium of the solid/vapor interphase is given by:

∇ T • τ τ τ = -t = (σ σ σ s -σ σ σ v ) • n v on Γ sv , (16)
where τ τ τ is the membrane stress tensor, . denotes the jump of a quantity across Γ sv , t denotes the traction vector and n v denotes the unit normal vector to Γ sv pointing into the vapor domain Ω v . Above, ∇ T • (.) is the surface divergence operator, which for any differentiable second order tensor A is given by:

∇ T • (A) = ∇ • (A) : 1 T (17)
with 1 T , the second order identity tensor of the plane T tangent to the interface defined by

1 T = 1 -n v ⊗ n v . ( 18
)
According to the coherent imperfect interface model, the displacement vector is continuous, u = 0 across the interface whereas the stress vector suffers a jump across the same interface t = 0 according to [START_REF] Perumal | A review on the development of lattice boltzmann computation of macro fluid flows and heat transfer[END_REF]. In particular the surface divergence operator applied to the tangent unity tensor ∇ T • 1 T = -(1 T : ∇n)n = cn gives rise to the signed local interface curvature c. We assume that γ sv is uniform on Γ sv , then:

(σ σ σ s -σ σ σ v ) • n v = -γ sv cn v on Γ sv . (19)
Above, we have introduced the membrane stress

τ = γ sv 1 T (20)
due to the presence of surface tensions. It is worth noting that the surface tensions γ sv and γ s are related by the Young's equation which defines the equilibrium of the contact line between solid-liquid-vapor domains,

γ sv -γ s = γ v cos(θ), ( 21
)
where θ is the contact angle between three interfaces. The above equations are completed by the action of pressure in the liquid and vapor phases on solid skeleton interface Γ s , which correspond finally to boundary conditions as:

σ s • n = -p n on Γ s , (22)
σ s • n v = -(p v + γ sv c)n v on Γ sv , ( 23
)
where p is the pressure exerted by the fluid on the solid, and p v is the pressure exerted by the vapor on the solid. Here, we have neglected the linear force effect on the triple line due to the surface tension γ v effects on the liquidvapor interface Γ v . We also consider a liquid phase completely wetting on the solid phase (θ = 0 o , γ sv = γ lv and γ s = 0) and the description was made by omitting the reference completely saturated state at atmospheric pressure to simplify the writing of the equations. For interested readers, a more detailed description can be found in [START_REF] Dormieux | Microporomechanics[END_REF].

Definition of pressure and surface tension from LBM simulations

In this section we define the values of p , p v , γ v and γ sv from the liquid density provided by the LBM solution.

From the SC LBM model described above, we can deduce in absence of body force, i.e. f g = 0, and at thermodynamical equilibrium state the local fluid stress tensor σ σ σ f (x) (subscript f refers to fluid), by using:

∇ • σ σ σ f (x) = -∇ • P 0 (x) + f tot (x) (24)
where

f tot (x) = -Gψ(x) 18 k=1 ω k ψ(x + ξ ξ ξ k ∆t)ξ ξ ξ k because f coh and f ads have the same form, with ψ(x + ξ ξ ξ k ∆t) = ψ(ρ α (x + ξ ξ ξ k ∆t)) if x + ξ ξ ξ k ∆t ∈ Ω α , α = s, f . We recall that P 0 (x) = ρ(x)c 2 s 1 .
By expanding in Taylor series the pseudopotential force f tot up to order 2 and integrating [START_REF] Vesselin | Simulating engineering flows through complex porous media via the lattice boltzmann method[END_REF], we get:

σ f ij (x) = -c 2 s ρ(x) + ∆t 2 c 2 s Gψ 2 (x bulk pressure )+ ∆t 3 2 c 4 s Gψ(x)∆ψ(x) + ∆t 3 4 c 4 s G|∇ψ(x)| 2 δ ij + ∆t 3 2 c 4 s G ∂ψ(x) ∂x i ∂ψ(x) ∂x j interfaces effect . (25)
At the thermodynamic equilibrium the density field is homogeneous in each bulk phase and evolves smoothly along a thick interface (δ = 0) between the phases. Therefore, the above fluid stress tensor which depends only on the local density field, contains two main contributions namely, the bulk pressure part in each bulk liquid, vapor phases

p(ρ) = c 2 s ρ + ∆t 2 c 2 s Gψ 2 (ρ) (26)
and the interfaces effects due to the presence of surface tensions γ αβ , αβ = s , sv, v between solid, liquid and vapor phases. In fact, without the bulk pressure part, all the other terms in the local fluid stress definition given in [START_REF] Matthew M Montemore | Effect of nanoscale flows on the surface structure of nanoporous catalysts[END_REF], describe the development of an interface profile with its own surface tension where the isotropy of the fluid stress tensor is violated. Indeed, along the interface thickness the fluid stress tensor becomes anisotropic with a mismatch between the transverse components and normal component (see Fig. 5).

According to [START_REF] Rowlinson | Molecular theory of capillarity[END_REF][START_REF] Benzi | Mesoscopic modeling of a two-phase in the presence of boundaries: the contact angle[END_REF], the solid, liquid and vapor surface tensions γ αβ , αβ = s , sv, v are defined as the integrals along the normal coordinate to the interface of the mismatch between normal σ f nn and transverse σ f tt components of the above fluid stress tensor. Therefore, the surface tension at the interface Γ αv with α = s, is,

γ αv = - Γ αv (σ f nn -σ f tt )dδ. ( 27
)
From the previous surface tensions definitions, we can recast the fluid pressure and surface tension into:

p = c 2 s ρ + ∆t 2 c 2 s Gψ 2 (ρ ), (28)
p v = c 2 s ρ v + ∆t 2 c 2 s Gψ 2 (ρ v ), (29)
γ αv = - 1 2 c 2 s G Γ αv | ∂ψ ∂δ | 2 dδ, α = s or , (30)
where ∂ ∂δ denotes derivative relative to the coordinate normal to the interface Γ v and dδ refers to differential along the normal coordinate of the interface.

Numerical definition of surfaces for prescribing pressure and surface tension

First, the geometry of the microstructure, including the solid phase and the voids, are meshed with tetrahedral elements (see Fig. 6 (a)). Then, a regular grid of voxels is defined on the domain Ω to solve the water phase equilibrium by LBM. The voxels of solid or void phases are set according to their belonging to corresponding elements. Then, the LBM simulations are conducted in voxels belonging to the void domain. At equilibrium, gas and liquid voxels are set according to their final density. Then, a map of distance function is constructed: for each fluid voxel, the minimal distance to any vapor voxel is computed. This distance function defined in each voxel is then interpolated at the middle point x g of each interface finite element between the solid and the porous phase (see Fig. 6 (c)). For each interface element, the minimal distance to the vapor phase d min (x g ) is then given. If d min (x g ) is lower than a threshold δ lim , then the interface element is set to belong to the imperfect interface associated to the thin film between the solid and the gas phase Γ sv ; otherwise the interface element belongs to Γ s . The algorithm is summarized in Algorithm 2 box below.

Algorithm 2: Construction of the liquid-solid interface Γ s and of the solid-vapor interface Γ sv . Input. Film thickness threshold δ lim , density field ρ(x) in Ω f , Γ s mesh (1) Detection of voxels in Ω and those in Ω v using the threshold ρ m = ρ l +ρ v 2 .

(2) For each fluid voxel, compute minimal distance d min (x) to any vapor voxel.

(3) Interpolate the obtained distance field at the center x g of interface finite elements between the solid and the porous domains.

(5) Set each Γ s element with d min (x g ) ≤ δ lim as belonging to Γ sv and the others to Γ s ;

Output. Interface elements belonging to Γ sv and Γ s .

The film thickness threshold δ lim is chosen to be equal to half of the liquidvapor interface thickness given by the SC LBM. In the SC model the interfaces are diffuse. From the model force balances it is delicate to adjust the interface thickness explicitly and independently of other parameters. The interface thickness depends not only on the temperature and τ = ∆t τc , but also on the state equation used to construct the interaction potential. For example, parameters a, R and b in the Carnahan and Starling's equation of state [START_REF] Carnahan | Equation of state for nonattracting rigid spheres[END_REF] may affect the interface thickness and numerical stability. In the present paper, the chosen parameters give us approximatively 5∆x (∆x denotes the SC LBM lattice spacing in the voxel grid) for the interface thickness for the most T and τ . Therefore, we have chosen δ lim = 2.5∆x.

In the coherent imperfect interface model, the mechanical impact of the presence of the surface tension γ sv = γ v at the solid-vapor interface is summed up by an overpressure whose intensity depends on the local signed mean curvature. Then, we propose to model the mechanical loading on the solid skeleton induced by the surface tension by applying this overpressure. The main difficulty to quantify this overpressure lies in the evaluation of the local signed mean curvature at each Gauss point of the solid-vapor interface mesh defined previously. For that purpose, we propose to compute the local signed mean curvature following the algorithm given as follows.

Algorithm 3: Solid-vapor interface Γ sv elements curvature algorithm.

Input. Interface solid-vapor Γ sv elements normal n e vectors (1) compute the elements area S e using the nodes coordinates.

(2) Compute mean Gauss coordinate x g of the element. Output. Interface Γ sv element signed mean curvatures c e .

The presence of two different grids, at the SC LBM grid and FEM mesh interface introduces a parameter F EM -LBM = F EM ∆x which is defined as the ratio of average solid-fluid interface FEM mesh edge length F EM to the SC LBM grid spacing ∆x. In this work, in order to manage between precision of the interface Γ s segmentation and local signed mean curvature calculation and computing time, we have chosen F EM -LBM = 1.

Construction of the homogenized model

Localization problem

We recall here the main equations at the microscale and define the so-called localization problem to be solved by FEM, given the microstructure mesh and the interfaces Γ sv and Γ s . The localization problem is given as follows.

Given a macroscopic strain tensor ε, a vapor pressure p v , a liquid pressure p and a surface tension γ sv , find the local strain field ε(x) in Ω such that:

∇ • σ s = 0 in Ω, (31) ε s (x) = ε (32) with σ s = C s : ε s in Ω s , ( 33
) and σ s • n v = -(p v + γ sv c)n v on Γ sv , (34)
σ s • n l = -p n l on Γ s . ( 35
) p = p v + γ sv c (γ sv = γ v ) on Γ v . ( 36
)
Condition ( 32) can be classically satisfied by the boundary condition (see e.g. [START_REF] Yvonnet | Computational Homogenization of Heterogeneous Materials with Finite Elements[END_REF])

u(x) = ε • x + ũ(x) on ∂Ω, (37)
where ũ(x) is a periodic displacement fluctuation function over Ω. Problem ( 31)-( 36) can be solved by classical finite elements.

Macroscopic model: effective quantities

Using the superposition principle, we can express the strain field solution ε(x) in Ω s as:

ε(x) = A(x) : ε + ε v (x)p v + ε (x)p + ε sv (x)γ sv ( 38
)
where A is the classical localization tensor relating local strain solution to macroscopic strains such that

A ijkl = ε (kl) ij , where ε (kl) ij
is the solution of ( 31)-(36) for p l = p v = γ sv = 0 with ε kl = 1 2 (e k ⊗ e l ), and where e i , i = 1, 2, 3 are basis vectors. From the averaging theorem (see e.g. [START_REF] Yvonnet | Computational Homogenization of Heterogeneous Materials with Finite Elements[END_REF]), we have:

σ = 1 V Ω σ(x)dΩ - 1 V Γ σ(x)n ⊗ xdΓ = 1 V Ω s C s : A(x) : ε + ε v (x)p v + ε (x)p + ε sv (x)γ sv - 1 V Ω v p v 1dΩ - 1 V Ω p 1dΩ - 1 V Γ sv ∪Γ lv cn v ⊗ xγ sv dΓ,
which can be re-written as:

σ = C : ε -B v p v -B p -B sv γ sv , (39)
with

C = (1 -φ)C s : 1 V s Ω s A(x)dΩ, (40)
B v = (φ -1)C s : 1 V s Ω s ε v (x)dΩ + φ(1 -S r )1 (41)
B = (φ -1)C s : 1 V s Ω s ε (x)dΩ + φS r 1 (42)
B sv = (φ -1)C s : 1 V s Ω s ε sv (x)dΩ - 1 V Γ sv ∪Γ lv cn v ⊗ xdΓ, (43)
where φ = V pore /V is the porosity with V pore being the volume of the pore and S r = V l /V pore is the saturation ratio. The Biot's moduli can be derived as [START_REF] Dormieux | Microporomechanics[END_REF]:

1 N v = (B v -φ1) : (C s ) -1 : 1, (44) 1
N = B -φ1 : (C s ) -1 : 1, (45)
1

N sv = (B sv -φ1) : (C s ) -1 : 1. (46)
In the particular macroscopic quantities B sv due to surface tension γ sv = γ v for θ = 0 o , averaging operation in Eq. ( 43) includes the contribution of the liquid-vapor interface. So, a similar algorithm like the solid-fluid interface segmentation Algorithm 2 have been proposed and implemented to capture the liquid-vapor Γ v interface effects. The singularity of the algorithm proposed here, is that we have first to reconstruct the mesh of the interface given by the fluid density threshold ρ m (see Algorithm 4). This reconstruction is done by smoothing the fluid density field obtained by the SC LBM and then applying the Marching-Cube algorithm [START_REF] Lorensen | Marching cubes: A high resolution 3d surface construction algorithm[END_REF].

Algorithm 4: liquid-vapor interface Γ lv mesh construction algorithm Input. Film thickness threshold δ im , density field ρ(x) in Ω f

(1) Smoothing the density ρ(x) field using a Gaussian filter;

(2) Detection of voxels in Ω and those Ω v using the threshold ρ m = ρ +ρ v 2 ;

(3) Extract the interface Γ v ∪Γ sv mesh given by the ρ m threshold using a marching-cube algorithm [START_REF] Lorensen | Marching cubes: A high resolution 3d surface construction algorithm[END_REF]. Γ sv is the part of the interface which is near of the solid interface; (4) Compute mean Gauss coordinate x g for each interface element;

(5) Compute minimal distance d min (x) to solid domain for each interface voxel in Ω and Ω v ; (6) Project the minimal distance to interface Γ v ∪ Γ sv mesh d min (x g ); (7) Each Γ v ∪ Γ sv element with d min (x) ≤ δ min is in Γ sv and the other in Γ v ;

Output. Interface Γ v mesh.

Prediction of shrinkage

The shrinkage strain ε sh can be obtained by using [START_REF] Shan | Lattice Boltzmann model for simulating flows with multiple phases and components[END_REF], neglecting the effect of vapor pressure p v and expressing ε for zero macroscopic stress. It results in:

ε sh = C -1 : B p + B sv γ sv . (47)
In the case where B = b 1 and B sv = b sv 1 due to the symmetries of the microstructure, then it is interesting to express the hydrostatic shrinkage strain which reduces in that case to:

ε sh = T r C -1 : 1 3 b p + b sv γ sv . ( 48
)

Numerical examples

As a classical usage in SC LBM simulations, the units in this paper are in lattice units, i.e., the velocity is given in lu/ts, the density in mu/lu 3 , the pressure in mu/(lu ts 2 ), temperature in tu and the surface tension in mu/ts 2 , etc. Here lu is the length unit, mu the mass unit and ts the time unit. Obtaining the unsaturated poroelastic properties in classical physical units can be performed by: x LBM lu → x (length), ρ LBM mu/lu 3 → ρ (mass) and γ lv LBM mu/ts 2 → γ lv (surface tension). The conversion mapping is completely defined by fixing the lattice length unit equal to the lattice physical spacing lu = ∆x (m), prescribing the physical liquid density and liquid-vapor surface tension, for example at a temperature of 20 o C we have ρ p = 1000 kg/m 3 and γ v p = 72.8 • 10 -3 N/m. Then, the SC LBM simulations are carried out with ∆t = 1 ts, ∆x = 1 lu. To maintain the stability of the simulations we have chosen τ c = ∆t in the sequel.

Validation of the LBM simulation model

Liquid and vapor phase separation

In this first example, we illustrate the capability of the SC LBM to model phase transition, e.g. from vapor to liquid. We recall that the critical state ( dp dρ = 0 and d 2 p dρ 2 = 0) of the fluid is given by T crit = 0.0943 tu, ρ crit = 0.1136 mu/lu 3 and p crit = 0.0044 mu/(lu ts 2 ). Below this critical temperature there is a density range for which dp dρ is negative, which yields the initiation of the presence of a transition in the liquid and vapor phases. The simulations are performed on a 128 × 128 × 4 lu 3 voxel grid with periodic boundary conditions. The initial density distribution is random but with an initial volume fraction S r = 4.32% of liquid density ρ l > ρ crit . Fig. 7 shows the time evolution of liquid and vapor phases separation from our SC LBM model implementation. We can appreciate the phase transition from vapor (white phase) into liquid (red phase). Thermodynamic equilibrium of the system has been reached for about t = 30000 ts.

Droplet test

In this test, we validate our LBM implementation and the evaluation of pressure by verifying the Laplace equation during the condensation of a liquid droplet. For that purpose, several circular droplets with different radii have been simulated. The computations have been carried out on a 128 × 128 × 4 lu 3 voxels grid at various reduced temperatures with an initial density field given by:

ρ(x, y, z) = ρ l + ρ v 2 - ρ l -ρ v 2 [tanh( 2( (x -x 0 ) 2 + (y -y 0 ) 2 -R 0 ) w )], (49)
with (x 0 , y 0 , z 0 ), R 0 and w being the droplet center coordinates, the droplet initial radius and the initial liquid vapor interface thickness, respectively (see Fig. 8). Periodic boundary condition have been prescribed. Fig. 9-a shows that the Laplace law can be reproduced by our LBM implementation. In Fig. 9-b the surface tension has been computed from Laplace law linear fit and by the integral formula given in [START_REF] Guihard | Assessing the permittivity of an unsaturated sand by combining a Lattice Boltzmann method simulation, electromagnetic homogenization models and measurements[END_REF]. We can note a general good agreement between both surface tension definitions even though there are discrepancies at the low reduced temperatures. As a last test to validate the SC LBM implementation, classical simulations for obtaining different contact angles (Fig. 10) have been carried out on a 128 × 64 × 4 lu 3 voxels grid at reduced temperature T /T c = 0.585. In these calculations, we have defined an initially semicircular static droplet with 15 lu of radius on a flat homogeneous solid surface. The simulations are performed with the top and bottom boundaries as solid wall and the left and right boundaries as periodic boundaries. Different contact angles have been obtained by changing the fictitious solid-fluid interface density ρ s as shown by figure 11-a. In fact, setting ρ s = ρ leads to θ = 0 o and setting ρ s = ρ v results in θ = 180 o . Changing ρ s in the range of (ρ v , ρ ) the contact angle can span the range (0 o , 180 o ). To retrieve the Young's equation, we proposed to compare the contact angles obtained by [START_REF] Martys | Simulation of multicomponent fluids in complex three-dimensional geometries by the Lattice Boltzmann method[END_REF] with the one provided by the geometric relationship:

r = (4h 2 + b 2 ) 8h , (50)
with

θ =      arctan[ b 2(r-h) ] θ ≤ 90 o π + arctan[ b 2(r-h) ] θ > 90 o , (51)
where b and h are the measured base and height of the droplet, after the simulation has converged. Fig. 11-a shows a good agreement between our SC LBM results and the analytical solution. In Fig. 11-b, we depict the surface tensions γ s , γ sv and γ v by changing ρ s . Only γ s and γ sv depend on ρ s and we satisfy γ s = 0 for ρ s = ρ which corresponds to the perfectly wetting liquid situation and γ sv = 0 for ρ s = ρ v which is the opposite situation. According to the Young's equation, we retrieve the fact that in perfectly wetting liquid assumption the liquid-vapor surface tension is equal to solid-vapor surface tension γ v = γ sv and the solid-liquid surface tension is then γ s = 0. In this section we propose a validation of the procedure described in Algorithm 3 box, which allows us to prescribe the force induced by the surface tension [START_REF] Shan | Kinetic theory representation of hydrodynamics: a way beyond the Navier-Stokes equation[END_REF]. We consider the case of a 3D spherical interface with its center located in the solid phase. For an elastic sphere blocked at its center, and prescribing overpressure (34) (for p v = 0) on its whole surface, the analytical expression for the radial displacement is given by [START_REF] Dormieux | Microporomechanics[END_REF]:

Validation of surface tension model and interface detection

u ref (r) = 2γ sv 3λ + 2µ r R , (52)
where λ and µ denote the Lamé's coefficients, and r and R the distance to the sphere origin and the radius of the sphere, respectively. The FEM computations have been carried out under Code Aster [54] for γ sv = 1. A mesh of tetrahedral elements has been constructed for several mesh densities, as shown in Fig. 12, where h denotes the element average size. A symmetry boundary condition on the internal faces of the mesh is prescribed and the center of the sphere is fixed.

Fig. 13 shows the radial displacement as a function of the distance to the origin of the sphere. We can note that our FEM calculations for different mesh densities are in good agreement with the reference solution as shown in Fig. 13 -a). In Fig. 13-b, the difference between the reference solution and the FEM calculations with respect to the distance to the center of curvature is depicted. Distance to origin r [-] Convergence of the solution with respect to the mesh density is shown in Fig. 14.

|u r -u ref r | [-] h = 0.25 h = 0.5 h = 1. h = 2. h = 4. (a) (b)

6.3

Validation of the complete model: capillary pressure in prolate pore

Capillary pressure

In this example, we validate the whole finite element procedure coupled with LBM by comparing the solution of capillary pressure in a prolate pore under non-saturated condition provided by our model and an approximate analytical solution [START_REF] Dormieux | Microporomechanics[END_REF]. Usually, the so-called capillary pressure curve p c = p c (S r ) is a macroscopic material property. Here, the liquid phase is assumed to be perfectly wet. The contact angle is θ = 0 and γ sv = γ v . The prolate pore has an aspect ratio α = 0.5 with radius R x = R y = 32 lu and R z = 64 lu. The LBM calculations have been carried out on a 66 × 66 × 132 lu 3 voxels grid at reduced temperature T /T c = 0.585 with ρ s = ρ (in order to achieve θ = 0 o , see Fig. 11), obtained in the case of flat interface for this temperature. Fig. 16 shows the fluid density field obtained for Sr = 0.054 and Sr = 0.40.

For this problem, the analytical solution under the spherical liquid meniscus assumption (Fig. 15-b) has been provided in [START_REF] Dormieux | Microporomechanics[END_REF] as:

p c = 2γ/R x 1 + (α 2 -1)u 2 (53)
with α = Rx Rz and u = z Rz , where z is defined in Fig. 15 (b). The analytical liquid saturation ratio is given by:

S r = 1 2 (1 -u ) 2 (2 + u ) -α( 1 + (α 2 -1)u 2 ) 3 + α 2 u ( 3 2 + u 2 (α 2 - 3 2 )).(54)
A comparison between both solutions is shown in Fig. 17. We observe the same order of magnitude but a shift in terms of saturation. The shift in the saturation may be due to the fact that the thickness of the fluid film is taken into account in the calculation of SC LBM liquid saturation ratio.

Homogenization calculations

In this section, homogenized quantities for a prolate microstructure are computed using the procedure described in section 5.2 and compared with available estimates. The porous medium is assumed to contain periodic prolate pores as depicted in Fig. 15. The solid phase of the microstructure has a linear and isotropic elastic behavior with a Young's modulus E s = 50 GPa and Poisson's ratio ν s = 0.25.

Liquid and vapor phases thermodynamics equilibrium calculations have already been presented in section 3.1. Firstly, for this microstructure, we propose to analyze the effective stiffness tensor C ef f , the saturated Biot tensor B and the liquid-vapor Biot's tensor in the case of empty pores in presence of surface tension B sv as a function of the porosity φ. This analysis is performed by comparing the results of the numerical approach with the classical analytical solutions1 [START_REF] Dormieux | Microporomechanics[END_REF]. Fig. 18 shows the normalized components of the effective stiffness tensor C ef f ijkl /C s ijkl as function of porosity φ. The tensor C ef f is orthotropic due to the morphology of the microstructure, but the available analytical solutions we use are only provided for transversally isotropic media. For porosity φ ≤ 10%, there is a good agreement between the different solutions. However, for porosity φ > 10%, the dilute scheme diverges, while Mori-Tanaka scheme and FEM remain close. The divergence between the dilute scheme and Mori-Tanaka solutions as the porosity increases reflects the fact that the diluted scheme does not take into account the interactions between the pores, unlike the Mori-Tanaka scheme (see e.g. [START_REF] Milton | Theory of Composites[END_REF]). In the FEM simulations, interactions between pores for arbitrary volume fractions is taken into account through the use of periodic boundary conditions. Fig. 19 shows the saturated Biot tensor components B ij as function of porosity φ. The computed saturated Biot's tensor exhibits transverse isotropic behavior and the FEM calculations remain in good agreement with the Mori-Tanaka solution.

Fig. 20 shows the solid-vapor Biot tensor components B sv ij as function of the porosity φ, in the case of empty pore with the presence of surface tension at the interface solid-pore. We note a good agreement between the different solutions for low porosity, while FEM and Mori-Tanaka solutions remain in good agreement for higher porosities.

Surface tension effects are more important when the surface-to-volume ratio of the microstructure is large (Fig. 20). These examples validate the implementations of the periodic boundary conditions on the solid phase, of the local signed curvature and membrane stress averaging algorithms.

Secondly, analysis is focused on the unsaturated materials properties namely, liquid Biot's tensor B and solid-vapor Biot's tensor B sv as function of the liquid saturation ratio S r , for a given porosity. Here porosity is fixed at φ = 0.1. The analytical methods cannot take into account a realistic morphology of the liquid and vapor phases in the calculation of the localization tensor under the action of theses phases. Then, the analytical approaches make the common assumption of isomorphology of these phases in the pore domain. This approximation can be expressed by:

B S r B. (55)
Fig. 21 shows the comparison between the components of the liquid phaseassociated Biot's tensor, calculated by the FEM-LBM methodology, and the Dilute and Mori-Tanaka schemes. According to the morphology of the porous network, we can note that the numerical tensor B exhibits the expected transverse isotropy. Nevertheless, the assumption made by the analytical estimators shows its limits since it is a roughly linear relationship between the liquid Biot's tensor and the degree of saturation unlike the LBM-FEM method. In addition, we find that analytical estimators underestimate the contribution of the liquid phase to the macroscopic behavior.

Moreover, in regard to the contribution of the surface tension effects to the macroscopic behavior, the isomorphology assumption (56) and the completely wetting liquid phase hypothesis, can be summarized by:

B sv φD : A f (56)
with

D = 1 V v Γ sv ∪Γ lv 1 T ds (57)
where D takes into account the effects on both, solid-vapor and liquid-vapor interfaces on the macroscopic behavior. The terms A f and V v denote the localization tensor average on the fluid domain and the pore volume, respectively.

It is interesting to compare the contributions of the volume effects p c T rB with the surface effects γ sv B sv to the macroscopic behavior of such a microstructure.

Given the capillary curve (Fig. 17) and all poroelastic properties in unsaturated condition of the periodic microsturcture containing prolate pores we can deduce the volume changes induced by drying for this microstructure and discuss the mechanisms involved (see Fig. 22). The variation in volume induced during drying in the absence of macroscopic stress has been given in Eq. 5.3. In this equation, shrinkage starts from an initial saturated state while here, the LBM-FEM or analytical calculations begin at the moment when the largest bubble is formed in the prolate pore.

In Fig. 22, shrinkage can be observed for this microstructure due to the negative sign of its absolute volume variation compared to the saturated reference configuration. In kinetic terms of this shrinkage (Fig. 22-a), we have an initial shrinkage at the beginning of the drying from the saturated reference state followed by a swelling up to 14% of liquid saturation ratio. Between 14% and 6% of liquid saturation ratio we again observe a shrinkage and then the material swells again between 6% and 0% of liquid saturation ratio. This shrinkage kinetics can be explained by analyzing the contribution of capillary and surface tension on the strain evolution (Fig. 22-b). From saturated reference configuration up to 14% of liquid saturation ratio, capillary pressure is the main loading compared to surface tension. However, between 0% and 14% of liquid saturation ratio, surface tension becomes the main loading. Now we can observe that the microstructure volume change induced by each loading, namely "volume effects" and "surface effects" on Fig. 22-a respectively, for capillary pressure and surface tension. For this microstructure morphology, capillary and surface tension loading respectively induce only shrinkage and swelling kinetics during drying. Finally, the initial shrinkage is due to capillary effects, the swelling kinetic up to 14% of liquid saturation ratio is due to the fact that the slope of surface tension swelling is greater than the one of the capillary pressure shrinkage. The opposite situation occurs between 14% and 6% of liquid saturation ratio. This results in apparent shrinkage kinetic and we retrieve the initial situation leading to the swelling kinetic between 6 % and 0% of liquid saturation ratio. For the contribution of the capillary bulk effects, we can mention that as the product S r p c controls the volumetric part of the shrinkage strain, we have zero shrinkage at 100% and 0%, so we necessarily go through (at least) one maximum between 100% and 40% of liquid saturation ratio.

Model microstructure: cubic centered lattice of spheres

We consider a model porous structure made of interpenetrable solid spheres arranged on a cubic centered periodic lattice, as depicted in Fig. 23. The side length of the RVE is 100 nm and the radius of the sphere is 23.57 nm. The microstructure is subjected to a drying/wetting cycle.

The drying/wetting cycle is performed through the capillary pressure paths shown in Fig. 24 and is obtained by LBM calculations on a lattice grid with 100 × 100 × 100 voxels. The porous microstructure has a porosity φ = 16.3%.

These isotherms are obtained by several calculations of LBM for a given liquid saturation ratio. For each saturation ratio, it is assumed that the maximum and minimum values of capillary pressure are on the desorption and adsorption paths, respectively. The simulations are performed at T /T crit = 0.585 reduced temperature, with completely wetting liquid phase assumption, i.e. the fictitious density of the solid phase is equal to the density of the liquid phase ρ s = ρ , with γ v = γ sv = 10.737 × 10 -3 mu/ts 2 . The elastic properties of the solid phase are, E s = 50 GP a for the Young modulus and ν s = 0.25 for the Poisson ratio.

For this nanometric porous system, Fig. 24 shows that we can reach at lower liquid saturation ratios some capillary pressures around 150 M P a. Fig. 24 shows also the drying/wetting paths defined from results of the present LBM calculations. We can observe a hysteresis of the curve of the capillary pressure

Conclusion

In this work, we have proposed a procedure combining Lattice-Boltzmann and finite element methods to model the effects of capillary pressure in porous micro structures. Combining FEM and LBM is introduced here for unsaturated porous media for the first time to our best knowledge. Starting from an explicit geometry of the microstructure, the Lattice-Boltzmann method (LBM) has been used to model the condensation from vapor phase to liquid and predict the geometry of capillary liquid films and liquid phases for arbitrary geometry of the solid skeleton. Then, a procedure has been elaborated to prescribe surface tension and pressure due to the capillary liquid films and liquid phases. The solid skeleton can deform linearly and elastically under the action of the liquid and its deformation can be captured using finite elements. Finally, a procedure to extract the shrinkage strain and the effective macroscopic behavior of the material, taking into account the explicit geometry of the microstructure and of the liquid capillary phases for a given saturation, has been provided.

The LBM solver has been validated on classical phase change and droplets wetting tests. Then, the solutions obtained by the present approach have been compared with available analytical solutions for simple case. These comparisons show a good agreement between them. Finally, the method has been applied to a more complex microstructures. Future works may include more realistic descriptions of the microstructure and more complex behaviors of the solid skeleton, as well as non-zero wetting angles. From a completely practical point of view, the proposed methodology in this paper could be used to discuss the phenomenological shrinkage models commonly used in civil engineering. Also, this work can be seen as a step towards a more global methodology allowing to discuss the potential interactions between creep of a solid poroviscoelastic matrix and its shrinkage induced by the capillary water loading. This last point corresponds to the configuration of the C-S-H gel in cementitious materials.

Aknowledgements

The authors gratefully acknowledge the financial support from EDF.

x y z L x L y L z R x R y R z Ω p Ω s Γ sp ∂Ω s α = Lx Lz = Ly Lz = Rx Rz = Ry Rz (a) r z O Γ lv Ω v Γ sv Ω l M M ′ C m R 0 Z 0 z ⋆ r ⋆ (

Figure 1 .

 1 Figure 1. Overall methodology: (a): geometry of the solid skeleton; (b) description of capillary water distribution after condensation simulated by Lattice Boltzmann Method; (c) simplified model containing solid, liquid phases and capillary films modeled as imperfect interfaces with surface tension.

Figure 2 .

 2 Figure 2. D 3 Q 19 molecule velocity space discretization scheme.

Figure 3 .

 3 Figure 3. Interior and exterior boundary conditions. (a) Half-way bounce-back scheme on Γ sf and (b) periodic boundary condition on ∂Ω f .

Figure 4 .

 4 Figure 4. (a) distribution of phases after condensation simulation; (b) simplified model with different surfaces and imperfect interfaces.

Figure 5 .

 5 Figure 5. Fluid stress tensor anisotropy (or mismatch) through the thick interface (δ = 0). Here δ is the interface thickness , dδ refers to the normal abscissa along the interface thickness, σ f nn and σ f tt are respectively the normal and transverse components of σ σ σ f .

Figure 6 .

 6 Figure 6. (a) Finite Element mesh of the microstructure; (b) voxel grid: grey-solid phase, blue-fluid phase; white-gas phase; (c) Computation of distances for fluid voxels to center of interface elements; (d) Definition of Γ s interface elements and Γ sv interface elements.

( 3 ) 1 e∈E

 31 Compute the elements tangent identity tensor 1 T = 1n e ⊗ n e ; (4) For each solid-vapor interface node, detect the set of all the elements sharing this node E = {e 1 , e 2 , • • • , e n }. (5) Compute smoothed elements normals at the node i using this average n i = Se e∈E n e S e . (6) Compute derivatives of shape function at each element mean Gauss point B(x g ). (7) Compute element normal gradient ∇n e = B(x g ) • (n 1 • • • n j ) T , for element with j nodes. (8) Compute elements signed mean curvature c e = -(1 T : ∇n e ).

Figure 7 .

 7 Figure 7. Simulation of condensation from vapor to liquid using SC LBM. White and red phases indicate vapor and liquid phases, respectively.

  (a) R 0 = 20 lu (b) R 0 = 35 lu (c) R 0 = 50 lu

Figure 8 . 2 • 10 - 2 Figure 9 .

 821029 Figure 8. Different final density fields in mu/lu 3 obtained at T /T c = 0.585 with w = 5lu and initial radius R 0 .

6. 1 . 3 Figure 10 . 1 • 10 - 2

 13101102 Figure 10. Different contact angles obtained through adjusting the fictitious solid density ρ s at T Tc = 0.585.

Figure 11 .

 11 Figure 11. Various contact angles, liquid-vapor, solid-vapor and solid-liquid surface tensions with fictitious solid density ρ s at T Tc = 0.585.

1 1 1 (Figure 12 .

 112 Figure 12. Mesh of a 1/8 sphere subjected to a surface tension.

Figure 13 .Figure 14 .

 1314 Figure 13. Elastic sphere subjected to surface tension: comparison between numerical and analytical solutions: (a) radial displacement; (b) difference between analytical and FEM solutions.

Figure 16 . 2 • 10 - 3 Figure 17 . 3333 FEM C ef f 1111 /C s 1111 FEM C ef f 2222 /C s 2222 FEM 1212 FEM C ef f 2323 /C s 2323 FEM C ef f 1313 /C s 1313 FEM 1122 FEM C ef f 2233 /C s 2233 FEM C ef f 1133 /C s 1133 FEM C ef f 1122 /C s 1122 Figure 18 . 23 Dilute B dil 13 Dilute B dil 12 Mori-Tanaka B mt 23 Mori-Tanaka B mt 13 Mori-Tanaka B mt 12 FEM B 23 FEM B 13 FEM B 12 Figure 19 . 2 33 FEM B sv 11 FEM B sv 22 FEMDilute B sv,dil 23 Dilute B sv,dil 13 Dilute B dil 12 Mori-Tanaka B sv,mt 23 Mori-Tanaka B sv,mt 13 Mori-Tanaka B sv,mt 12 FEM B sv 23 FEM B sv 13 FEM B sv 12 Figure 20 . 33 LBM-FEM B l 11 LBM-FEM B l 22 LBM 3 Dilute B l,dil 23 Dilute B l,dil 13 Dilute B l,dil 12 Mori-Tanaka B l,mt 23 Mori-Tanaka B l,mt 13 Mori-Tanaka B l,mt 12 LBM-FEM B l 23 LBM-FEM B l 13 LBM-FEM B l 12 Figure 21 .Figure 22 .Figure 25 . 4 • 10 - 2

 162103173333111122221212232313131122223311331122182313122313122313121923311222313122313122313122033112232313122313122313122122254102 Figure 16. Prolate pore in unsaturated condition: cut plane of the density field at two saturation states S r .

Figure 26 .

 26 Figure 26. Unsaturated Biot coefficient as function of liquid saturation ratio for a periodic cubic centred microstructure.

Figure 27 .

 27 Figure 27. Volume change induced by drying of a periodic cubic centred porous microstructure.

Dilute and Mori-Tanaka schemes

as function of the liquid saturation ratio. Fig. 25 shows the liquid phase distribution at the same liquid saturation ratio (S r 48.9 %) for both, drying (Fig. 25-a) and wetting (Fig. 25-b) paths. We can observe a different distribution according to the saturation evolution and get insight into the hysteresis of the capillary curve by some ink bottle effects as proposed by [START_REF] Coussy | Poromechanics[END_REF]. Indeed, for a given liquid saturation ratio we can see that the liquid phase is in the large pores with anchorage (solid-liquid-vapor system triple line) in the small pores on the drying path and the opposite situation occurs on the wetting path. Thus, here the capillary pressure is higher on the drying path as compared to the wetting one using Laplace law, because the curvature radius located in the small pores are smaller.

The poroelastic properties in unsaturated condition of this microstructure are obtained by the LBM-FEM methodology described in section 5.2. An illustration of the mesh used to compute them is provided in Fig. 23 (b). Due to the symmetry of this microstructure, the Biot tensors are isotropic. The stiffness tensor C ef f has cubic symmetry and is entirely determined by

From now on, we denote by b, b and b sv the Biot's coefficients of an effective isotropic medium associated with the pressure in the pore, the pressure of the liquid on Γ s and the surface tension on Γ sv . We obtain the saturated Biot's coefficient as b = 0.461. Fig. 26 shows the contributions of the bulk liquid phase and the interfaces to macroscopic behavior, respectively. The local distribution of menisci of both, drying and wetting paths have been considered in the calculations of these liquid and liquid-vapor Biot's coefficients. Moreover, we adopt for this cubic centered microstructure the assumption of isomorphology in the fluid phases presented previously. Indeed, we observe a quasi-linear relationship between the liquid Biot b l and liquid-vapor Biot b sv coefficients respectively, with S r and (1 -S r ).

Finally, Fig. 27 shows the strain or volume change induced by drying or wetting as well as a comparison between bulk capillary pressure and surface tension effects. Nevertheless, in terms of kinetics we have at the beginning of the drying from the saturated reference state shrinkage kinetic and then up to 0 % of liquid saturation ratio there is a swelling kinetic for this microstructure (see Fig. 27-a). These observations can be explained by analyzing the bulk capillary and surface tension effects. Fig. 27-b shows that the surface tension effects dominate as compared to bulk capillary pressure effects at lower liquid saturation ratio. The opposite situation is observed at higher liquid saturation ratio. This can be explained by the fact that even if the capillary pressure increases at low saturation, the application surface thereof is considerably reduced, which leads to low macroscopic effects at low saturations.