David Grégoire
email: david.gregoire@univ-pau.fr

Carine Malheiro

Christelle Miqueu
email: christelle.miqueu@univ-pau.fr

Estimation of adsorption-induced pore pressure and confinement in a nanoscopic slit pore by a density functional theory

Keywords: pore pressure, adsorption, microporous, density functional theory, methane, carbon

Estimation of adsorption-induced pore pressure and confinement in a nanoscopic slit pore by a density functional theory

Introduction

Following the IUPAC recommendation [1,2], the pore space in porous materials is divided into three groups according to the pore size diameters: macropores of widths greater than 50 nm, mesopores of widths between 2 and 50 nm and micropores 1 of widths less than 2 nm. Zeolites, activated carbons , tight rocks, coal rocks, source rocks, cement paste or construction materials are among microporous materials. In recent years, a major attention has been paid on these microporous materials because the surface-to-volume ratio (i.e., the specific pore surface) increases with decreasing characteristic pore size. Consequently, these materials can trap an important quantity of fluid molecules as an adsorbed phase. This is important for applications in petroleum and oil recovery, gas storage, separation, catalysis or drug delivery.

In very small pores, the molecules of fluid are confined (molecular packing). This effect induces that fluid-fluid and fluid-solid interactions sum at the pore scale and has significant consequences at the macroscale, such as instantaneous swelling or shrinkage deformation. In different contexts, these deformations may be critical. A lot of natural and synthesised porous media are composed of a double porosity: the microporosity where the fluid is trapped as an adsorbed phase and a meso or a macro porosity required to ensure the transport of fluids to and from the smaller pores. If adsorption in nanopores induces instantaneous deformations at a higher scale, the matrix swelling may close the transport porosity, reducing the global permeability of the porous system or annihilating the functionality of synthesised materials. Particularly, in situ adsorption-induced coal swelling has been identified [3,4] as the principal factor leading to a rapid decrease in CO 2 injectivity during coal bed methane production enhanced by CO 2 injection. Conversely, gas desorption can lead to matrix shrinkage and microcracking, which may help oil and gas recovery in the context of unconventional petroleum engineering [5]. The effects of adsorbent deformation on physical adsorption has also been identified [6] as the major challenge concerning gas porosimetry in nano-porous non-rigid materials (e.g. metal organic frameworks). In conclusion, there is now a consensus in the community that major attention has to be focused on the coupled effects appearing at the nanoscale within microporous media because they may have significant consequences at the macroscale [7].

Experimentally, different authors tried to combine gas adsorption results and volumetric swelling data (see e.g. [7] for a review). Pioneer works of Meehan [8] showed the effect of carbon dioxide sorption on the expansion of charcoal but only mechanical deformations were reported and no adsorption quantities were measured. Later on, different authors [9,5,10,11,12,13,14] highlighted the phenomenon on bituminous coals, because it is of utmost importance in the context of CO 2 geological sequestration and coal bed reservoirs exploitation. Recently, Perrier et al. [15,[START_REF] Perrier | Coupling between adsorption and deformation in microporous media[END_REF] performed simultaneous measurements of adsorbed quantities and induced strain on an activated carbon and showed that the volumetric swelling is directly linked to the excess adsorbed quantity.

As far as modeling is concerned, standard poromechanics [START_REF] Biot | [END_REF][START_REF] Coussy | Poromechanics[END_REF] links the macroscopic strain to the total stress and the pore fluid pressure (i.e. the bulk pressure in that case) and has led to relevant descriptions of various multiphysics phenomena in macroporous materials. However, this classical approach fails in describing adsorption-induced strain and predicts volumetric shrinkage instead of volumetric swelling for these microporous materials [START_REF] Perrier | [END_REF]. Such classical framework needs to be extended to take into account the fluid confinement effects present within the nanopores. Especially, molecular packing induces a deviation of the pore fluid pressure from its bulk value, the difference being most often denoted as solvation force or pressure [START_REF] Israelachvili | Intermolecular and surface forces[END_REF]. The effect of this modified pore fluid pressure on the pore walls has to be characterized at the pore scale and upscaled to rep-resent macroscopic swelling strain. Molecular simulations are the classical tools to represent the distribution of molecules due to fluid-fluid and fluid-solid interactions at the nanoscale. Different authors used Monte Carlo simulation tools to study the pore fluid pressure profiles within nanopores. Do and coworkers performed Grand Canonical Monte Carlo simulations for argon [START_REF] Do | [END_REF]22], methane and methanol [23] in slit-shaped micropores with movable solid layers to compute the pore pressures and corresponding deformations. Kowalczyk and coworkers studied adsorption-induced deformation of microporous carbons filled with carbon dioxide [24,25] and near-critical argon [26]. These investigations showed that pressures applied on the pore surfaces may be very high (few hundred of MPa), depending on the thermodynamic conditions and on the pore sizes. However, even if macroscopic adsorption isotherms may be reconstructed in a consistent way from molecular simulations through the material pore size distribution [27], molecular simulations are not a tractable nor an efficient tool to predict resulting deformations at a macroscale in a complex material. Hence, enhanced poromechanics framework have been developed [START_REF] Perrier | [END_REF]28,[START_REF] Perrier | Coupling between adsorption and deformation in microporous media[END_REF]29] to obtain semi-analytical efficient and tractable tools capable to predict macroscopic strain induced by adsorption in homogeneous microporous materials. In these formulations, the fluid confinement is macroscopically deduced from experimental measurements of adsorbed quantities and not from the thermodynamic state of the fluid at the nanoscale. For micro-to-macro heterogeneous porous media, the fluid confinement cannot be deduced from experimental measurement and new strategies have to be proposed.

If molecular simulations are too time consuming to be coupled with poromechanical models or macroscopic numerical simulations, an elegant and relevant alternative approach consists in estimating the fluid confinement with a molecular density functional theory (DFT) 2 . It has been already shown [START_REF] Malheiro | [END_REF][START_REF] Malheiro | The journal of physical chemistry[END_REF] that nonlocal DFT is as accurate as molecular simulations to estimate density profiles of real fluids in nanopores but with a very lower computational cost. Another key feature is that DFT is written in a continuum framework and then may be coupled more easily with continuum poromechanics theory than discrete simulations results. Balbuena et al. [34] proposed a theoretical study based on DFT to estimate the solvation pressure of Lennard-Jones fluids in slit micropores. At low pressure and low temperature conditions, Neimark and co-workers [35,36,37] or Ustinov and Do [38] used different DFT variations to estimate the solvation pressure (e.g. for porous media characterisation purposes). To our knowledge, only one article in the literature proposes to apply this strategy on a large range of temperature and pressure [39] by studying the adsorption of methane on natural coals at geological conditions.

This paper aims at using a non-local density functional theory combined with a powerful equation-of-state in order to obtain highly-resolved pore pressure maps as a function of bulk pressure and pore width at different temperatures. The proposed method allows to better characterise the complex pore pressure behaviour for large thermodynamic conditions on the whole range of micropore widths even for complex fluids. This is the first step to provide a unified theory from nano to macro scales that links adsorption of the fluid in the nanopores to the macroscopic deformation of the material.

Theory

Brief overview of NLDFT

The version of the Non-Local Density Functional Theory (NLDFT) used in this work has been fully described in Malheiro et al. [START_REF] Malheiro | [END_REF]. Therefore, the discussion on this theory will be limited to its most significant features here.

Let us consider a fluid of chemical potential µ at temperature T in a volume V . In the density functional theory (DFT) framework [START_REF] Evans | Fundamentals of Inhomogeneous Fluids[END_REF], the grand thermodynamic potential Ω of an inhomogeneous fluid in an external field in the grand canonical ensemble (µ, V, T ) is given by:

Ω[ρ] = A[ρ] - N i=1 V ρ i (x)[µ bi -V ext i (x)]dx . (1)
In Eq. ( 1), Ω is the grand thermodynamic potential, A is the Helmholtz free energy, ρ = (ρ 1 , ..., ρ N ) is the vector of the different molecular densities of the N species of the inhomogeneous fluid, which depend on the position vector x, V ext i is an external field (e.g. a fluid/solid interaction potential), µ bi is the bulk chemical potential of the specie i.

The minimization of this grand potential leads to the Euler-Lagrange equations that allow computing the density distributions of the fluid at equilibrium:

δΩ[ρ] δρ i = δA[ρ] δρ i -µ bi + V ext i = 0 ∀i ∈ [1, N ] . (2)
In Eq. ( 2), (Ω, A are the volumetric grand thermodynamic potential and the volumetric Helmholtz free energy respectively. In this work, the SAFT-VR equation of state [40] has been used as the thermodynamic model to compute both the properties of the bulk fluid phase and the chemical potential at equilibrium and as the bulk limit to build the inhomogeneous Helmholtz free energy. In this work, we deal with methane, a fluid of non-associating spherical monomers of diameter σ, with a mean-field (MF) approximation, i.e. the correlations due to the long-range attractive interactions are neglected. In that case, the Helmholtz free energy functional is given by:

A[ρ] = V A[ρ(x)]dx = A ideal [ρ] ideal part + A HS [ρ] hard sphere contribution + A lr 1 [ρ] + A sr 1 [ρ] + A 2 [ρ] dispersive part . (3)
In Eq. (3):

-A ideal is the ideal part and is exactly known as:

A ideal [ρ] = k B T N i=1 V ρ i (x) ln[ρ i (x)Λ 3 i ] -1 dx . (4)
In Eq. ( 4), k B is the Boltzmann constant and Λ i is the thermal de Broglie wavelength.

-The hard sphere contribution A HS is computed non locally with the modified fundamental measure theory (MFMT) [41,42] taking into account scalar and vectorial contributions:

A HS [ρ] = k B T V Φ HS,s [ρ(x)] + Φ HS,v [ρ(x)] dx (5) with      Φ HS,s = -n 0 ln(1 -n 3 ) + n 1 n 2 1-n 3 + n 3 2 ln(1-n 3 ) 36πn 2 3 Φ HS,v = n 4 n 5 1-n 3 - n 2 n 5 n 5 ln(1-n 3 ) 12πn 2 3 - n 2 n 5 n 5 12πn 3 ln(1-n 3 ) 2 (6) and            ∀α ∈ [0, 3], n α [ρ(x)] = N i=1 V ρ i (ξ)w α ( x -ξ )dξ ∀β ∈ [4, 5], n β [ρ(x)] = N i=1 V ρ i (ξ)w β (x -ξ)dξ (7)
where

   w 0 (r) = 1 πσ 2 δ( σ 2 -r) , w 1 (r) = 1 2πσ δ( σ 2 -r) , w 2 (r) = δ( σ 2 -r) , w 3 (r) = H( σ 2 -r) , w 4 (r) = r 2πσr δ( σ 2 -r) , w 5 (r) = r r δ( σ 2 -r) . (8)
In Eq. ( 8), δ and H are respectively the Dirac and the Heaviside functions.

-The dispersive part encompasses a first-order and a second-order monomer perturbation. The first-order monomer perturbation has been partitioned to get on the one hand the "short-range" contribution A sr 1 and on the other hand the "long-range" contribution A lr 1 that becomes the classical attractive term of the DFT framework in the mean-field approximation [43]. We get

A lr 1 [ρ] = 1 2 N i=1 N j=1 V ρ(x) V ρ(ξ)Φ att ij ( x -ξ )dξ dx . (9)
In Eq. ( 9), Φ att ij is the attractive part of the square-well potential. The "shortrange" contribution A sr 1 and the second-order monomer perturbation A 2 are estimated through a coarse-grained approach (see [START_REF] Malheiro | [END_REF] for full details) from a weighted density:

ρ(x) = 3 4πσ 3 V ρ(ξ)H σ -x -ξ dξ . (10)
In the following, we are considering a slit pore of width L made of two solid parallel graphitic surfaces and filled with methane (see Fig. 1). The external potential V ext applied on the fluid is the sum of the Steele's interaction potentials [44]: In Eq. (11), ρ s = 114 molecules/m 3 is the density of graphite, ∆ = 0.335 nm is the space between two layers of graphite constituting the wall [45]. (σ, ε) are respectively the collision diameter and the depth of the well potential. σ sf and ε sf are the solid-fluid interaction parameters determined by the Lorentz-Berthelot rules

   V ext (z) = V steele (z) + V steele (L -z) V steele (z) = 2πε sf ρ s σ 2 sf ∆ 0.4 σ sf z 10 - σ sf z 4 - σ 4 sf 3∆(z+0.61∆) 3 (11)
σ sf = σ ss +σ f f 2 and ε sf = √ ε ss ε f f . (σ f f , ε f f
) are the fluid monomer parameters in the SAFT-VR theory and (σ ss , ε ss ) are the solid parameters.

As shown in Malheiro et al. [START_REF] Malheiro | [END_REF] and reprinted in Fig. 2, molecular simulations have been used to validate the methane's density profiles obtained with the NLDFT framework presented above in the full micropores range.

Pore pressure estimation

The grand thermodynamic potential of a fluid confined in a slit-like pore of width L between two parallel plates of surface A 0 being [START_REF] Hansen | Theory of Simple Liquids[END_REF]:

dΩ = -SdT - i N i dµ i + 2γdA 0 -Π f A 0 dL , (12)
the pressure Π f is:

Π f = - 1 A 0 ∂Ω ∂L T,A 0 ,µ i . ( 13
)
Thus defined, Π f becomes equal to the bulk pressure P b when L → ∞. It is the sum of two contributions Π f = P b + f s , f s being the solvation pressure as introduced by Hansen and McDonald [START_REF] Hansen | Theory of Simple Liquids[END_REF] that vanishes when L → ∞. One can demonstrate (see Appendix A for a complete development) from Eqs. ( 1) and ( 13) that:

Π f = - L 0 ρ(z) ∂V ext (z, L) ∂L dz . (14)
In case of a Steele-type external potential (see Eq. ( 11)), we get also:

Π f = - L 0 ρ(z) ∂V steele ∂z (z)dz . (15)
Eq. ( 15) is in agreement with Eq. (6.5.5) of Hansen and McDonald's book [START_REF] Hansen | Theory of Simple Liquids[END_REF] or [34]. Thus, the pressure inside a slit-like pore can be computed with Eq. ( 15) as soon as the density profiles have been computed with the NLDFT.

Results and discussion

The pore pressure Π f has been computed in the whole range of micropores ( L < 2 nm), for pressures up to 20 MPa at three temperatures : 303, 373 and 443 K. As an example, the 3D map of Π f = f (L, P b ) at 373 K is plotted in Fig. 3 whereas the 3D maps corresponding to the other temperatures are plotted in appendix B (Fig. 10). Animated versions of these 3D maps are given as supplementary materials. These maps Π f = f (L, P b ) clearly highlight that the pore pressure is a complex function of both pore size and thermodynamic conditions and that its variation is far from being obvious. Indeed, it depends on an intricate interplay between fluid-solid and fluid-fluid interactions; the first ones being highly dependent on the distance between pore walls (given an adsorbent-adsorbate pair) and the latest ones on the thermodynamic conditions (pressure and temperature). As already observed in some molecular simulations [START_REF] Do | [END_REF]23] and DFT [35,39] calculations in idealised slit or cylindrical pores, the pore pressure oscillates as a function of pore width as a result of commensurate or incommensurate packing. As mentioned in the Introduction, this behaviour has already been presented in the literature with molecular simulations. However, the calculation speed of thermodynamic theories such as NLDFT compared to molecular simulations allows drawing complete maps of pore pressure for wide thermodynamic conditions and pore ranges. Pore pressures can reach extremely high values in the ultramicropores range. For the largest micropores, the difference between the pore and bulk pressures becomes negligible, the pores being filled with methane at a nearly bulk state, except for the layers on contact with the solid walls. Moreover, the pore pressure successively increases or decreases with bulk pressure depending on the considered pore width range, as better illustrated in the 2D plots of Fig. 4. As shown in Fig. 4 for L < 0.72 nm, the pores exhibit only positive pore pressure, and so positive adsorption stress causing expansion. This range of widths corresponds to pores able to accommodate only one fluid layer. On the contrary, for L > 0.72, i.e. for pores accommodating at least two overlapping layers, one can observe a transition from initial negative pore pressure to positive ones with increasing bulk pressures. This behaviour reflects the competition between fluid-solid interactions predominant at low pressures and fluid-fluid interactions that dominate at the highest pressures. Finally, one can see in Fig. 4 that, at a given bulk pressure, Π f successively oscillates as a function of pore width, thus reflecting packing effects for such parallel ideal solid surfaces. This behaviour is better illustrated in Fig. 5 where Π f is plotted as a function of pore width at different bulk pressures at T = 373 K.

As can be seen in Fig. 5, when L increases, the pore pressure exhibits an oscillatory decay similar to the one of the pair correlation function, giving successively repulsive or attractive force on the surface of the solid walls. Even if the oscillations of Π f reflect the layers formation in the case of infinite parallel smooth walls, the same behavior has been observed in slit micropores of finite length [START_REF] Kopanichuk | [END_REF] or of heterogeneous surface [39]. For the largest micropores (L > 1.5 nm), the packing effects become negligible and there is no more additional stress applied on the surfaces. A quantitative comparison cannot be done with the previous results of Kopanichuk et al. [START_REF] Kopanichuk | [END_REF] and Yang et al. [39], due to the differences between both the fluid and the solid representations. However, the solvation pressure evolutions pre- sented here have features in perfect qualitative agreement with the works of these authors. Indeed, on the one hand, Kopanichuk et al. performed GCMC simulations in finite length graphitic slit pores. They obtained the same oscillatory but less pronounced behavior of the solvation pressure because, as they mentioned, fluid adsorption in narrow finite slits is partially depressed in compared to the infinite ones. On the other hand, in the work of Yang et al., methane is represented in a simpler manner (LJ pair-wise interactions) than it is done with a NLDFT-SAFT formalism. The principal consequence is that the densities computed by Yang et al. are underestimated at high pressure and so is the solvation pressure. For instance, the maximum average density obtained in this work at (T = 303 K, P = 200 bar) is six times greater and consequently the solvation pressure ten times greater than the ones obtained by Yang et al. in similar conditions. This oscillatory decay has been previously observed for the solvation pressure in liquids, theoretically and experimentally by different authors [START_REF] Israelachvili | Intermolecular and surface forces[END_REF]48,49,50]. Israelachvili [START_REF] Israelachvili | Intermolecular and surface forces[END_REF] proposed an exponential oscillatory decay with a pseudo-period close to σ f f and a characteristic decay length equal to α × σ f f with α ∈ {1.2; 1.7} for liquids. Fig. 6 presents the evolution of the estimated solvation pressure with the normalised pore width and one can observe that the same type of exponential decay is obtained here at T = 373 K and P bulk = 200 bar for supercritical methane. A characteristic exponential decay factor of α = 1.16 is recovered here. Therefore, one can consider that for pore widths larger that 5 × α × σ f f ≈ 2 nm, the pore pressure does not oscillate anymore and that the bulk pressure is fully recovered. Moreover, Fig. 7 shows that the dependence of the characteristic exponential decay factor α with the temperature is quite linear. To the authors' knowledge, this was never reported in the literature but only three temperatures are investigated here and further results are needed to confirm this trend.

In addition, it is noteworthy that the extrema of the pore pressure do not exactly correspond to the extrema of the average density in the pore (see Fig. 8). Indeed, such a correspondence should only happen in the case of hard walls as observed by Henderson [48] and the slight difference of position between the extrema of pore pressure and average density is a result of the fluid-solid potential V ext . In the literature, the maxima of Π f are often attributed to the maxima of density occurring for the most ordered states in commensurate pores and the minima for 1.17 the least ordered states in non-commensurate pores [START_REF] Israelachvili | Intermolecular and surface forces[END_REF]22] maybe as a result of less complete maps of pore pressure vs pore widths. This emphasizes the usefulness of the calculation speed of DFTs. Fig. 8 highlights also the influence of temperature for the three different temperatures considered here (T ∈ {443, 373, 303} K). It has already been shown in the literature [51] that the density profiles are less structured and the mean density maxima are less pronounced when temperature increases. This behaviour may be attributed to the increase of kinetics contribution in the free energy. We can observe here exactly the same behaviour on the pore pressure maxima. With increasing temperature, they are less and less pronounced. Moreover one can notice that the shift between pore pressure and mean density profiles also decreases with increasing temperature. This means that this shift is mostly driven by short range interactions whose influence decreases when temperature and thus molecular agitation increase.

In order to have a better understanding of the relation between the pore density profiles and the induced attractive or repulsive pore pressure, it is useful to compute the work of the pore pressure as this latter allows to distinguish stable and unstable configurations, as shown in Fig. 9. Hence, it is demonstrated here that the pore pressure work can be used to identify unambiguously the pore widths that are commensurate (or non-commensurate) with the size of an adsorbent.

Concluding remarks

In this paper, we have shown that:

-A non-local density functional theory can be used to obtain highly resolved 3D maps of pore pressure versus pore widths and bulk pressure of methane in graphitic slit-shaped pores at different temperatures. -These maps clearly highlight that, even for a non polar fluid such as methane, the pore pressure reaches very high values in the narrowest nanopores and is a complex function of both pore width and thermodynamic conditions. -For ultramicropores accommodating only one fluid layer, the pores exhibit only positive pore pressure, and so positive adsorption stress causing expansion. For larger pores accommodating at least two overlapping layers, one can observe a transition from initial negative pore pressures to positive ones with increasing bulk pressures. This behaviour reflects the competition between fluid-solid interactions -predominant at low pressure -and fluid-fluid interactions -that dominate at higher pressure. -For a given bulk pressure and with increasing pore widths, an oscillatory decay of the pore pressure is recovered giving successively repulsive or attractive force on the surface of the solids walls. The dependence of the characteristic exponential decay factor is linear with temperature in the range analysed here and is the order of 1.2 × σ f f . -It has been also clearly highlighted that the pore pressure maxima differ from the mean density ones. -This oscillatory decay reflects both the molecular packing and the fluid confinement. The pore pressure work distinguishes stable and unstable configurations. Therefore, it can be used to unambiguously identify the pore widths that are commensurate (or not-commensurate) with the size of an adsorbent. These conclusions emphasise the need to have a reliable continuous thermodynamic model with a low computational cost to compute the pore fluid properties and their impact on the solid matrix in order to be coupled with a poromechanical model. A non-local DFT coupled with the SAFT equation-of-state can play this role for complex fluids on wide thermodynamic conditions.

Supplementary Materials

Animated versions of Fig. 10.

A Development of the pore pressure estimation

In the context of a slit pore of width L and surface A 0 , filled with pure methane, Eqs. ( 1) and ( 13) may be rewritten as:

Ω(ρ) = A(ρ) -A 0 L/2 -L/2 ρ(z) µ b -V ext (z) dz and Π f = -1 A 0 ∂Ω ∂L T,A 0 ,µ . (16)
Therefore,

Π f = - 1 A 0 ∂A ∂L T,µ,A 0 + ∂ ∂L L/2 -L/2 ρ(z, L) µ b -V ext (z, L) dz I 1 (17)
On one hand, considering µ as a constant during the process of the partial derivation by L, we get:

∂A ∂L T,µ,A 0 = 0 (18)
and ∂ρ ∂L T,µ,A 0 = 0

On a second hand, considering the substitution z = ϕ(u) = uL, we get:

Π f = I 1 = u= z L ∂ ∂L 1/2 -1/2 ρ(uL, L) µ b -V ext (uL, L) Ldu (20) = 1/2 -1/2 ∂ ∂L ρ(uL, L) µ b -V ext (uL, L) L du (21) = 1/2 -1/2 ∂ ∂L [ρ(uL, L)µ b L] du I 2 - 1/2 -1/2 ∂ ∂L ρ(uL, L)V ext (uL, L)L du I 3 . (22)
•

I 2 = 1/2 -1/2 ρ(uL, L)µ b du I 4 + 1/2 -1/2 µ b L ∂ρ(uL, L) ∂L du I 5 (23)
where

I 4 = z=uL L/2 -L/2 ρ(z, L)µ b L dz = µ b L L/2 -L/2 ρ(z, L)dz (24) and
where

I 7 = z=uL L/2 -L/2 ρ(z, L)V ext (z, L) L dz , (32)
For a symmetrical external potential such as the Steele-type potential (see Eq. ( 11)), we get also:

Π f = -

Fig. 1 Fig. 2

 12 Fig. 1 Schematic view of a slit pore filled with methane.

Fig. 3

 3 Fig.3Pore pressure as a function of bulk pressure and pore size at 373 K.

Fig. 4

 4 Fig.4Pore pressure computed with Eq. (15) as a function of bulk pressure at 373 K depicted for various pore ranges.

Fig. 5

 5 Fig. 5 Pore pressure evolution as a function of pore size for four different values of bulk pressure (P bulk ∈ {1, 10, 100, 200} bar) at 373 K.

Fig. 6

 6 Fig. 6 Solvation pressure evolution as a function of normalised pore widths at T = 373 K and P bulk = 200 bar.

Fig. 7

 7 Fig. 7 Evolution of the characteristic exponential decay factor with the temperature at P bulk = 200 bar.

  dz with L * = L -2e and e = 0.79 × σ sf

Fig. 8

 8 Fig.8Comparison between the pore pressure evolution and the average pore density evolution as a function of pore size for three different temperatures (T ∈ {443, 373, 303} K) at 200 bar.

FFig. 9

 9 Fig. 9 Pore pressure work evolution as a function of pore widths (top) and characteristic pore density profiles (bottom, A to F) at 373 K.

I 2 = 30 )• I 3 = 1 / 2 - 1 / 2 ρ 7 + 1 / 2 - 1 / 2 V 8 + 1 / 2 - 1 / 2 ρ

 230312127121281212 I 4 + I 5 = 0 . ((uL, L)V ext (uL, L)du I ext (uL, L)L ∂ρ(uL, L) ∂L du I (uL, L)L ∂V ext (uL, L) ∂L du I 9

2 - 1 / 2 VI 3 = I 7 + I 8 + I 9 =Π f = I 1 = I 2 -I 3 = - L/ 2 -L/ 2 ρ

 21237891222 ext (uL, L)L u ∂ρ(uL, L) ∂z + ∂ρ(uL, L) (z, L) ∂V ext (z, L) ∂L dz

Fig. 10

 10 Fig. 10 Pore pressure as a function of bulk pressure and pore size at 303 K (top), 373 K (center) and 443 K (bottom).

The term nanopores is sometimes preferred in the Mechanics community to label pores of width less than

nm. However, the 2015 IUPAC report[2] recommends that this term embraces the three categories of pores, but with an upper limit ≈ 100 nm.

A comprehensive description of DFT and its applications can be found in the review of Evans[START_REF] Evans | Fundamentals of Inhomogeneous Fluids[END_REF] or the monograph of Davis[START_REF] Davis | Statistical mechanics of phases, interfaces and thin films[END_REF].

Acknowledgements Financial supports from the Communauté de Communes de Lacq-Orthez and Conseil Départemental 64 are gratefully acknowledged. We also gratefully acknowledge Dr. Mourad Abouzaid for his advices and our different discussions, especially concerning the mathematical developments.


Rıza Alp Güler
email: riza.guler@inria.fr2g.trigeorgis

George Trigeorgis

Epameinondas Antonakos
email: antonak@amazon.com

Patrick Snape

Stefanos Zafeiriou
email: s.zafeiriou@imperial.ac.uk

Iasonas Kokkinos
email: i.kokkinos@cs.ucl.ac.uk

DenseReg: Fully Convolutional Dense Shape Regression In-the-Wild

In this paper we propose to learn a mapping from image pixels into a dense template grid through a fully convolutional network. We formulate this task as a regression problem and train our network by leveraging upon manually annotated facial landmarks "in-the-wild". We use such landmarks to establish a dense correspondence field between a three-dimensional object template and the input image, which then serves as the ground-truth for training our regression system. We show that we can combine ideas from semantic segmentation with regression networks, yielding a highly-accurate 'quantized regression' architecture.

Our system, called DenseReg, allows us to estimate dense image-to-template correspondences in a fully convolutional manner. As such our network can provide useful correspondence information as a stand-alone system, while when used as an initialization for Statistical Deformable Models we obtain landmark localization results that largely outperform the current state-of-the-art on the challenging 300W benchmark. We thoroughly evaluate our method on a host of facial analysis tasks, and demonstrate its use for other correspondence estimation tasks, such as the human body and the human ear. DenseReg code is made available at http://alpguler.com/DenseReg.html along with supplementary materials.

Introduction

Non-planar object deformations, e.g. due to facial pose or expression, result in challenging but also informative signal variations. Our objective in this paper is to recover this information in a feedforward manner by employing a discriminatively trained convolutional network. Motivated by the gap between discriminatively trained systems for detection and category-level deformable models, we propose a system that combines the merits of both.

In particular, discriminative learning-based approaches typically pursue invariance to shape deformations, for in- * E. Antonakos was with Imperial College London during this work. We introduce a fully convolutional neural network that regresses from the image to a "canonical", deformationfree parameterization of the face surface, effectively yielding a dense 2D-to-3D surface correspondence field. Once this correspondence field is available, one can effortlessly solve many image-level problems by backward-warping their canonical solution from the template coordinates to the image domain for the problems of landmark localization, semantic part segmentation, and face transfer.

stance by employing local 'max-pooling' operations to ellicit responses that are invariant to local translations. As such, these models can reliably detect patterns irrespective of their deformations through efficient, feedforward algorithms. At the same time, however, this discards useful shape-related information and only delivers a single categorical decision per position. Several recent works in deep learning have aimed at enriching deep networks with information about shape by explicitly modelling the effect of similarity transformations [START_REF] Papandreou | Modeling local and global deformations in deep learning: Epitomic convolution, multiple instance learning, and sliding window detection[END_REF] or non-rigid deformations [START_REF] Jaderberg | Spatial transformer networks[END_REF]18,9]; several of these have found success in classification [START_REF] Papandreou | Modeling local and global deformations in deep learning: Epitomic convolution, multiple instance learning, and sliding window detection[END_REF], fine-grained recognition [START_REF] Jaderberg | Spatial transformer networks[END_REF], and also face detection [9]. There are works [START_REF] Lades | Distortion invariant object recognition in the dynamic link architecture[END_REF]36] that model the deformation via optimization procedures, whereas we obtain it in a feedforward manner and in a single shot. In these works, shape is treated as a nuisance, while we treat it as the goal in itself. Recent works on 3D surface correspondence [START_REF] Masci | Geodesic convolutional neural networks on riemannian manifolds[END_REF]4] have shown the merit of CNN-based unary terms for correspondence. In our case we tackle the much more challenging task of establishing a 2D to 3D correspondence in the wild by leveraging upon recent advances in semantic segmentation [START_REF] Chen | Semantic image segmentation with deep convolutional nets and fully connected crfs[END_REF]. To the best of our knowledge, the task of explicitly recovering dense correspondence in the wild has not been addressed yet in the context of deep learning.

By contrast, approaches that rely on Statistical Deformabe Models (SDMs), such as Active Appearance Models (AAMs) or 3D Morphable Models (3DMMs) aim at explicitly recovering dense correspondences between a deformation-free template and the observed image, rather than trying to discard them. This allows to both represent shape-related information (e.g. for facial expression analysis) and also to obtain invariant decisions after registration (e.g. for identification). Explicitly representing shape can have substantial performance benefits, as is witnessed in the majority of facial analysis tasks requiring detailed face information e.g. landmark localisation [START_REF] Sagonas | 300 faces in-the-wild challenge: Database and results[END_REF], 3D pose estimation, as well as 3D face reconstruction "in-the-wild" [22], where SDMs consitute the current state of the art.

However SDM-based methods are limited in two respects. Firstly they require an initialization from external systems, which can become increasingly challenging for elaborate SDMs: both AAMs and 3DMMs require at least a bounding box as initialization and 3DMMs may further require position of specific facial landmarks. Furthermore, SDM fitting requires iterative, time-demanding optimization algorithms, especially when the initialisation is far from the solution. The advent of Deep Learning has made it possible to replace the iterative optimization task with iterative regression problems [START_REF] Trigeorgis | Mnemonic descent method: A recurrent process applied for end-to-end face alignment[END_REF], but this does not alleviate the need for initialization and multiple iterations.

In this work we aim at bridging these two approaches, and introduce a discriminatively trained network to obtain, in a fully-convolutional manner, dense correspondences between an input image and a deformation-free template coordinate system.

In particular, we exploit the availability of manual facial landmark annotations "in-the-wild" in order to fit a 3D template; this provides us with a dense correspondence field, from the image domain to the 2-dimensional, U -V parameterization of the face surface. We then train a fully convolutional network that densely regresses from the image pixels to this U -V coordinate space. This provides us with dense and fine-grained correspondence information, as in the case of SDMs, while at the same time being independent of any initialization procedure, as in the case of discriminatively trained 'fullyconvolutional' networks. We demonstrate that the perfor-mance of certain tasks, such as facial landmark localisation or segmantic part segmentation, is largely improved by using the proposed network.

Even though the methodology is general, this paper is mainly concerned with human faces. The architecture for the case of human face is described in Fig. 1.

Our approach can be seen in two complementary manners: first, it provides a stand-alone, feedforward alternative to the combination of initialization with iterative fitting typically used in SDMs. This allows us to have a feedforward system that solves both the detection and correspondence problems at approximate 7 -8 frames per second for a 300 × 300 input image. Secondly, our approach can also be understood as an initialization procedure for SDMs which gets them started from a much more accurate position than the bounding box, or landmark-based initializations currently employed in the face analysis literature. When taking this approach we observe substantial gains over the current state-of-the-art systems.

We can summarize our contributions as follows:

We introduce the task of dense shape regression in the setting of CNNs, and exploit the SDM-based notion of a deformation-free UV-space to construct target ground-truth signals (Sec.2).

We propose a carefully-designed fully-convolutional shape regression system that exploits ideas from semantic segmentation and dense regression networks. Our quantized regression architecture (Sec.3) is shown to substantially outperform simpler baselines that consider the task as a plain regression problem.

We use dense shape regression to jointly tackle a multitude of problems, such as landmark localization or semantic segmentation. In particular, the template coordinates allow us to 'copy' multiple annotations constructed on a single template system, and thereby tackle multiple problems in a single go.

We use the regressed shape coordinates for the initialization of SDMs; systematic evaluations on facial analysis benchmarks show that this yields substantial performance improvements on tasks ranging from landmark localization to semantic segmentation.

We demonstrate the generic nature of the method by applying it to the task of estimating dense correspondence for human bodies and ears.

From SDMs to Dense Shape Regression

Following the deformable template paradigm [START_REF] Yuille | Deformable templates for face recognition[END_REF][START_REF] Grenander | Hands: A Pattern Theoretic Study of Biological Shapes[END_REF], we consider that object instances are obtained by deforming a prototypical object, or 'template', through dense deformation fields. This makes it possible to factor object variability within a category into variations that are associated to deformations, generally linked to the object's 2D/3D shape, and variations that are associated to appearance (or, 'texture' in graphics), e.g. due to facial hair, skin color, or illumination. This factorization largely simplifies the modelling task. SDMs use it as a stepping stone for the construction of parametric models of deformation and appearance. For instance, in AAMs a combination of Procrustes Analysis, Thin-Plate Spline warping and PCA is the standard pipeline for learning a low-dimensional linear subspace that captures category-specific shape variability [START_REF] Cootes | Active appearance models[END_REF]. Even though we have a common starting point, rather than trying to construct a linear generative model of deformations, we treat the image-to-template correspondence as a vector field that our network tries to regress.

In particular, we start from a template X = [x 1 , x 2 , ...x m ] ∈ R, where each x j ∈ R 3 is a vertex location of the mesh in 3D space. This template could be any 3D facial mesh, but in practice it is most useful to use a topology that is in correspondence with a 3D statistical shape model such as [2] or [START_REF] Paysan | A 3d face model for pose and illumination invariant face recognition[END_REF]. We compute a bijective mapping ψ, from template mesh X to the 2D canonical space U ∈ R 2×m , such that

ψ(x j ) → u j ∈ U , ψ -1 (u j ) → x j . (1)
The mapping ψ is obtained via the cylindrical unwrapping described in [START_REF] Booth | Optimal UV spaces for facial morphable model construction[END_REF]. Thanks to the cylindrical unwrapping, we can interpret these coordinates as being the horizontal and vertical coordinates while moving on the face surface:

u h j ∈ [0, 1] and u v j ∈ [0, 1].
Note that this semantically meaningful parameterization has no effect on the operation of our method.

We exploit the availability of landmark annotations "in the wild", to fit the template face to the image by obtaining a coordinate transformation for each vertex x j . We use the fittings provided by [START_REF] Zhu | Face alignment across large poses: A 3d solution[END_REF] which were fit using a modified 3DMM implementation [START_REF] Romdhani | Estimating 3d shape and texture using pixel intensity, edges, specular highlights, texture constraints and a prior[END_REF]. However, for the purpose of this paper, we require a per-pixel estimate of the location in UV space on our template mesh and thus do not require an estimate of the projection or model parameters as required by other 3D landmark recovery methods [22,[START_REF] Zhu | Face alignment across large poses: A 3d solution[END_REF].

The per-pixel UV coordinates are obtained through rasterization of the fitted mesh and non-visible vertices are culled via z-buffering. As illustrated in fig:GT, once the transformation from the template face vertices to the morphed vertices is established, the u j coordinates of each visible vertex on the canonical face can be transferred to the image space. This establishes the ground truth signal for our subsequent regression task.

Fully Convolutional Dense Shape Regression

Having described how we establish our supervision signal, we now turn to the task of estimating it through a convolutional neural network (CNN). Our aim is to estimate at any image pixel that belongs to a face region the values of u = [u h , u v ]. We need to also identify non-face pixels, e.g. by predicting a 'dummy' output.

One can phrase this problem as a generic regression task and attack it with the powerful machinery of CNNs. Unfortunately, the best performance that we could obtain this way was quite underwhelming, apparently due to the task's complexity. Our approach is to quantize and estimate the quantization error separately for each quantized value. Instead of directly regressing u, the quantized regression approach lets us solve a set of easier sub-problems, yielding improved regression results.

In particular, instead of using a CNN as a 'black box' regressor, we draw inspiration from the success of recent works on semantic part segmentation [START_REF] Tsogkas | Deep learning for semantic part segmentation with high-level guidance[END_REF][START_REF] Chen | Deeplab: Semantic image segmentation with deep convolutional nets, atrous convolution, and fully connected crfs[END_REF], and landmark classification [5,6]. These works have shown that CNNs can deliver remarkably accurate predictions when trained to predict categorical variables, indicating for instance the facial part or landmark corresponding to each pixel.

Building on these successes, we propose a hybrid method that combines a classification with a regression problem. Intuitively, we first identify a coarser face region that can contain each pixel, and then obtain a refined, region-specific prediction of the pixel's U -V field. As we will describe below, this yields substantial gains in performance when compared to the baseline of a generic regression system.

We identify facial regions by using a simple geometric approach. We tesselate the template's surface with a cartesian grid, by uniformly and separately quantizing the u h and u v coordinates into K bins, where K is a design parameter. For any image that is brought into correspondence with the template domain, this induces a discrete labelling, which can be recovered by training a CNN for classification.

On Fig. 4, the tesselations of different granularities are visualized. For a sufficiently large value of K even a plain classification result could provide a reasonable estimate of the pixel's correspondence field, albeit with some staircasing effects. The challenge here is that as the granularity of these discrete labels becomes increasingly large, the amount of available training data decreases and label complexity increases. A more detailed analysis on the effect of labelspace granularity to segmentation performance is provided in supplementary materials.

We propose to combine powerful classification results with a regression problem that will yield a refined correspondence estimate. For this, we compute the residual between the desired and quantized U -V coordinates and add a separate module that tries to regress it. We train a separate regressor per facial region, and at any pixel only penalize the regressor loss for the responsible face region. We can interpret this form as a 'hard' version of a mixture of regression experts [21]. This interpretation is further elaborated upon in the supplementary material.

The horizontal and vertical components u h , u v of the correspondence field are predicted separately. This results in a substantial reduction in computational and sample complexity -For K distinct U and V bins we have K 2 regions; the classification is obtained by combining 2 K-way classifiers. Similarily, the regression mapping involves K 2 regions, but only uses 2K one-dimensional regression units. The pipeline for quantized face shape regression is provided in Fig. 3.

We now detail the training and testing of this network; for simplicity we only describe the horizontal component of the mapping. From the ground truth construction, every position x is associated with a scalar ground-truth value u h . Rather than trying to predict u h as is, we transform it into a pair of discrete q h and continuous r h values, encoding the quantization and residual respectively:

q h = u h d , r h i = u h i -q h i d , (2)
where d = 1 K is the quantization step size (we consider u h , u v coordinates to lie in [0, 1]).

Given a common CNN trunk, we use two classification branches to predict q h , q v and two regression branches to predict r h , r v as convolution layers with kernel size 1 × 1. As mentioned earlier, we employ separate regression functions per region, which means that at any position we have K estimates of the horizontal residual vector, rh i , i = 1, . . . , K.

At test time, we let the network predict the discrete bin qh associated with every input position, and then use the respective regressor output rh qh to obtain an estimate of u:

ûh = qh d + rh qh (3)
For the q h and q v , which are modeled as categorical distributions, we use softmax followed by the cross entropy loss. For estimating rh and rv , we use a normalized version of the smooth L 1 loss [START_REF] Girshick | Fast r-cnn[END_REF]. The normalization is obtained by dividing the loss by the number of pixels that contribute to the loss.

Compared to plain regression of the coordinates, the proposed method achieves much better results. In Fig. 5 we report results of an experiment that evaluates the contribution of the q-r branches separately for different granularities. The results for the quantized branch are evaluated by transforming the discrete horzintal/vertical label into the center of the region corresponding to the quantized horizontal/vertical value respectively. The results show the merit of adopting the classification branch, as the finely quantized results(K=40,60) yield better coordinate estimates with respect to the non-quantized alternative (K=1). After K=40, we observe an increase in the failure rate for the quantized branch. The experiment reveals that the proposed quantized regression outperforms both non-quantized and the best of only-quantized alternatives. K=1, q+r (0.748) K=5, q+r (0.800) K=10, q+r (0.810) K=20, q+r (0.808) K=5, q (0.319) K=10, q (0.560) K=20, q (0.713) K=40, q (0.772) K=60, q (0.779) Figure 5: Performance of q and r, branches for various tesselation granularities, K. Areas under the curve(AUC) are reported.

Experiments

Herein, we evaluate the performance of the proposed method (referred to as DenseReg) on various face-related tasks. In the following sections, we first describe the training setup (Sec. 4.1) and then present extensive quantitative results on (i) semantic segmentation (Sec. 4.2), (ii) landmark localization on static images (Sec. 4.3), (iii) deformable tracking (Sec. 4.4), (iv) dense correspondence on human bodies (Sec. 4.5), and (v) human ear landmark localization (Sec. 4.6). Due to space constraints, we refer to the supplementary material for additional qualitative results, experiments on monocular depth estimation and further analysis of experimental results.

Training Setup

Training Databases. We train our system using the 3DDFA data of [START_REF] Zhu | Face alignment across large poses: A 3d solution[END_REF]. The 3DDFA data provides projection and 3DMM model parameters for the Basel [START_REF] Paysan | A 3d face model for pose and illumination invariant face recognition[END_REF] + Face-Warehouse [START_REF] Cao | Facewarehouse: A 3d facial expression database for visual computing[END_REF] model for each image of the 300W database. We use the topology defined by this model to define our UV space and rasterize the images to obtain per-pixel ground truth UV coordinates. Our training set consists of the LFPW trainset, Helen trainset and AFW, thus 3148 images that are captured under completely unconstrained conditions and exhibit large variations in pose, expression, illumination, age, etc. Many of these images contain multiple faces, some of which are not annotated. We deal with this issue by employing the out-of-the-box DPM face detector of Mathias et al. [START_REF] Mathias | Face detection without bells and whistles[END_REF] to obtain the regions that contain a face for all of the images. The detected regions that do not overlap with the ground truth landmarks do not contribute to the loss. For training and testing, we have rescaled the images such that their largest side is 800 pixels.

CNN Training. For the dense regression network, we adopt a ResNet101 [19] architecture with dilated convolutions (atrous) [START_REF] Chen | Semantic image segmentation with deep convolutional nets and fully connected crfs[END_REF]29], such that the stride of the CNN is 8. We use bilinear interpolation to upscale both the q and r branches before the losses. The losses are applied at the input image scale and back-propagated through interpolation. We apply a weight to the smooth L1 loss layers to balance their contribution. In our experiments, we have used a weight of 40 for quantized (d = 0.1) and a weight of 70 for non-quantized regression, which are determined by a coarse cross validation. We initialize the training with a network pre-trained for the MS COCO segmentation task [START_REF] Lin | Microsoft coco: Common objects in context[END_REF]. The new layers are initialized with random weights drawn from Gaussian distributions. Large weights of the regression losses can be problematic at initialization even with moderate learning rates. To cope with this, we use initial training with a lower learning rate for a warm start for a few iterations. We then use a base learning rate of 0.001 with a polynomial decay policy for 20k iterations with a batch size of 10 images. During training, each sample is randomly scaled with one of the ratios [0.5, 0.75, 1, 1.25, 1.5] and cropped to form a fixed 321 × 321 input image.

Semantic Segmentation

As discussed in Sec. 2, any labelling function defined on the template shape can be transferred to the image domain using the regressed coordinates. One application that can be naturally represented on the template shape is semantic segmentation of facial parts. To this end, we manually defined a segmentation mask of 8 classes (right/left eye, right/left eyebrow, upper/lower lip, nose, other) on the template shape, as shown in Fig. 1. We compare against a state-of-the-art semantic part segmentation system (DeepLab-v2) [START_REF] Chen | Deeplab: Semantic image segmentation with deep convolutional nets, atrous convolution, and fully connected crfs[END_REF] which is based on the same ResNet-101 architecture as our proposed DenseReg. We train DeepLab-v2 on the same training images (i.e. LFPW trainset, Helen trainset and AFW). We generate the ground-truth segmentation labels for both training and testing images by transferring the segmentation mask using the ground-truth deformation-free coordinates explained in Sec. 2. We employ the Helen testset [START_REF] Le | Interactive facial feature localization[END_REF] for the evaluation.

Table 1 reports evaluation results using the intersectionover-union (IoU) ratio. Additionally, Fig. 6 shows some qualitative results for both methods, along with the groundtruth segmentation labels. The results indicate that the DenseReg outperforms DeepLab-v2. The reported improvement is substantial for several parts, such as eyebrows and lips. We believe that this result is significant given that DenseReg is not optimized for the specific task-at-hand, as opposed to DeepLab-v2 which was trained for semantic segmentation. This performance difference can be justified by the fact that DenseReg was exposed to a richer label structure during training, which reflects the underlying variability and structure of the problem.

Landmark Localization on Static Images

DenseReg can be readily used for the task of facial landmark localization on static images. Given the landmarks' locations on the template shape, it is straightforward to estimate the closest points in the deformation-free coordinates on the images. The local minima of the Euclidean distance between the estimated coordinates and the landmark coordinates are considered as detected landmarks. In order to find the local minima, we simply analyze the connected components separately. Even though more sophisticated methods for covering "touching shapes" can be used, we found that this simplistic approach is sufficient for the task.

Note that the closest deformation-free coordinates among all visible pixels to a landmark point is not necessarily the correct corresponding landmark. This phenomenon is called "landmark marching" [START_REF] Zhu | High-fidelity pose and expression normalization for face recognition in the wild[END_REF] and mostly affects the jaw landmarks which are dependent on changes in head pose. It should be noted that we do not use any explicit supervision for landmark detection nor focus on ad-hoc methods to cope with this issue. Errors on jaw landmarks due to invisible coordinates and improvements thanks to deformable models can be observed in Fig. 8.

Herein, we evaluate the landmark localization performance of DenseReg as well as the performance obtained by employing DenseReg as an initialization for deformable models [START_REF] Boscaini | Learning shape correspondence with anisotropic convolutional neural networks[END_REF]44,[START_REF] Antonakos | Feature-based lucas-kanade and active appearance models[END_REF][START_REF] Trigeorgis | Mnemonic descent method: A recurrent process applied for end-to-end face alignment[END_REF] trained for the specific task. In the second scenario, we provide a slightly improved initialization with very small computational cost by reconstructing the detected landmarks with a PCA shape model that is constructed from ground-truth annotations.

We present experimental results using the very challenging 300W benchmark. This is the testing database that was used in the 300W competition [START_REF] Sagonas | 300 faces in-the-wild challenge: The first facial landmark localization challenge[END_REF][START_REF] Sagonas | 300 faces in-the-wild challenge: Database and results[END_REF] -the most important facial landmark localization challenge. The error is measured using the point-to-point RMS error normalized with the interocular distance and reported in the form of Cumulative Error Distribution (CED). Figure 7 (bottom) presents some self-evaluations in which we compare the quality of initialization for deformable modelling between DenseReg and two other standard face detection techniques (HOG-SVM [START_REF] King | Max-margin object detection[END_REF], DPM [START_REF] Mathias | Face detection without bells and whistles[END_REF]). The employed deformable models are the popular generative approach of patch-based Active Appearance Models (AAM) [START_REF] Boscaini | Learning shape correspondence with anisotropic convolutional neural networks[END_REF]44,[START_REF] Antonakos | Feature-based lucas-kanade and active appearance models[END_REF], as well as the current state-of-the-art approach of Mnemonic Descent Method (MDM) [START_REF] Trigeorgis | Mnemonic descent method: A recurrent process applied for end-to-end face alignment[END_REF]. It is interesting to notice that the performance of DenseReg without any additional deformable model on top, already outperforms even HOG-SVM detection combined with MDM. Especially when DenseReg is combined with MDM, it greatly outperforms all other combinations. [START_REF] Uřičář | Multi-view facial landmark detector learned by the structured output svm[END_REF] 0.2109 32.17

Table 2: Landmark localization results on the 300W testing dataset using 68 points. Accuracy is reported as the AUC and the Failure Rate.

Figure 7 (top) compares DenseReg+MDM with the results of the latest 300W competition [START_REF] Sagonas | 300 faces in-the-wild challenge: Database and results[END_REF]. We greatly outperform all competitors by a large margin. It should be noted that the participants of the competition did not have any restrictions on the amount of training data employed and some of them are industrial companies (e.g. Fan et al. [15]), which further illustrates the effectiveness of our approach. Finally, Table 2 reports the area under the curve (AUC) of the CED curves, as well as the failure rate for a maximum RMS error of 0.1. Apart from the accuracy improvement shown by the AUC, we believe that the reported failure rate of 3.67% is remarkable and highlights the robustness of DenseReg.

Deformable Tracking

For the challenging task of deformable face tracking on lengthy videos, we employ the testing database of the 300VW challenge [START_REF] Shen | The first facial landmark tracking inthe-wild challenge: Benchmark and results[END_REF][START_REF] Chrysos | Offline deformable face tracking in arbitrary videos[END_REF] -the only existing benchmark for deformable tracking "in-the-wild". The benchmark consists of 114 videos (∼ 218k frames in total) and includes videos captured in totally arbitrary conditions (severe occlusions and extreme illuminations). The tracking is performed based on sparse landmark points, thus we follow the same strategy as in the case of landmark localization in Sec. 4.3. We compare the output of DenseReg, as well as DenseReg+MDM which was the best performing combination for landmark localization in static images (Sec. 4.3), against the participants of the 300VW challenge. Table 3 reports the AUC and Failure Rate measures. DenseReg combined with MDM demonstrates better performance than the winner of the 300VW competition. It should be highlighted that our approach is not fine-tuned for the task-athand as opposed to the rest of the methods that were trained on video sequences and most of them make some kind of temporal modelling. Finally, similar to the 300W case, the participants were allowed to use unlimited training data (apart from the provided training seuqences), as opposed to DenseReg (and MDM) that were trained only on the 3148 images mentioned in Sec. 4.1. Please refer to the supplementary material for a more detailed presentation of the tracking results.

Dense Correspondence for the Human Body

To portray that the DenseReg system can be used for articulated shapes of complex topology, we present results on the human shape. We use the recently proposed "Unite the People" (UP) dataset [25], which provides a 3D deformable human shape model [START_REF] Loper | Smpl: A skinned multi-person linear model[END_REF] in correspondence with images from several publicly available datasets. We handle the complex geometry of the human shape by manually partitioning the surface into patches. We unwrap each patch using multidimensional scaling. The partitioning replaces the quantization and the rest of the system remains the same. Since there are no dense correspondence results between a 3D human model and image pixels in literature, we demonstrate the performance of our system through visual results from our test-set partition of the UP dataset in Fig. 9.

Ear Landmark Localization

We have also performed experiments on the human ear. We employ the 602 images and sparse landmark annotations that were generated in a semi-supervised manner by Zhou et al. [START_REF] Zhou | Estimating correspondences of deformable objects "in-the-wild[END_REF]. Due to the lack of a 3D model of the human ear, we apply Thin Plate Splines to bring the images into dense correspondence and obtain the deformationfree space. We perform landmark localization following the same procedure as in Sec. 4.3. Quantitative results are detailed in the supplementary material, where we compare DenseReg, DenseReg+AAM and DenseReg+MDM with alternative DPM detector based initializations. We observe that DenseReg results are highly accurate and clearly outperforms the DPM based alternative even without a deformable model. Examples for dense human ear correspondence estimated by our system are presented in Fig. 10.

Estimated

Groundtruth Estimated Groundtruth Estimated Groundtruth

Conclusion

We propose a fully-convolutional regression approach for establishing dense correspondence fields between objects in natural images and three-dimensional object templates. We demonstrate that the correspondence information can successfully be utilised on problems that can be geometrically represented on the template shape. Throughout the paper, we focus on face shapes, where applications are abundant and benchmarks allow a fair comparison. We show that using our dense regression method out-of-thebox outperforms a state-of-the-art semantic segmentation approach for the task of face-part segmentation, while when used as an initialisation for SDMs, we obtain the state-ofthe-art results on the challenging 300W landmark localization challenge. We demonstrate the generality of our method by performing experiments on the human body and human ear shapes. We believe that our method will find ubiquitous use, since it can be readily used for face-related tasks and can be easily integrated into many other correspondence problems.

1 Figure 1 :

 11 Figure1: We introduce a fully convolutional neural network that regresses from the image to a "canonical", deformationfree parameterization of the face surface, effectively yielding a dense 2D-to-3D surface correspondence field. Once this correspondence field is available, one can effortlessly solve many image-level problems by backward-warping their canonical solution from the template coordinates to the image domain for the problems of landmark localization, semantic part segmentation, and face transfer.

Ψ 1 Figure 2 :

 12 Figure 2: Ground-truth generation: (a) Annotated landmarks. (b) Template shape morphed based on the landmarks. (c) Deformation-free coordinates (u h and u v ), obtained by unwrapping the template shape, transferred to image domain.

Figure 3 :

 3 Figure 3: Proposed Quantized Regression Approach for the horizontal correspondence signal: The continuous signal is regressed by first estimating a grossly quantized (or, discretized) function through a classification branch. For each quantized value qh we use a separate residual regression unit's prediction, rh qh , effectively multiplexing the different residual predictions. These are added to the quantized prediction, yielding a smooth and accurate correspondence field.

Figure 4 :

 4 Figure 4: Horizontal and vertical tesselations obtained using K = 2, 4 and 8 bins.

Figure 6 :

 6 Figure 6: Exemplar semantic segmentation results.

Figure 7 :

 7 Figure 7: Landmark localization results on the 300W testing dataset using 68 points. Accuracy is reported as Cumulative Error Distribution of RMS point-to-point error normalized with interocular distance. Top: Comparison with state-of-the-art. Bottom: Self-evaluation results.

Figure 9 :

 9 Figure 9: Dense Correspondence for human body.

Figure 10 :

 10 Figure 10: Exemplar pairs of deformation-free coordinates of dense landmarks on human ear.

Table 1 :

 1 Semantic segmentation accuracy on Helen testset measured using intersection-over-union (IoU) ratio.

	Class	Methods DenseReg Deeplab-v2
	Left Eyebrow	48.35	40.57
	Right Eyebrow	46.89	41.85
	Left Eye	75.06	73.65
	Right Eye	73.53	73.67
	Upper Lip	69.52	62.04
	Lower Lip	75.18	70.71
	Nose	87.71	86.76
	Other	99.44	99.37
	Average	71.96	68.58

Table 3 :

 3 Deformable tracking results against the state-ofthe-art on the 300VW testing dataset using 68 points. Accuracy is reported as AUC and the Failure Rate.

Acknowledgements Rıza Alp Güler and Iasonas Kokkinos were supported by the European Horizon 2020 programme under grant agreement no 643666 (I-Support). G. Trigeorgis was supported by EPSRC DTA award at Imperial College London. The work of P. Snape and E. Antonakos was partially funded by an EPSRC DTA and by the European Community Horizon 2020 [H2020/2014-2020] under grant agreement no. 688520 (TeSLA). The work of S. Zafeiriou was funded by the EPSRC project EP/N007743/1 (FACER2VM).


Gilles Régnier

Célia Nicodeau

Jacques Cinquin

Jacques Verdu

Francisco Chinesta

Gilles Regnier

Eads

A Multiphysic and Multiscale Approach to Model the Continuous Welding of Thermoplastic Matrix Composites

Keywords: thermoplastic composite tape, welding, chain diffusion, degradation, crystallization, simulation Roller (Tr

published or not. The documents may come

Introduction

Thermoplastic matrix composites are more and more used in the aircraft industry for structural parts. In the case of continuous fiber-reinforced materials, a common way of processing consists to weld elementary tapes one to another to constitute the desired part. A two steps procedure was needed in the past: first a manual or automated tape placement, then a consolidation in autoclave. The aim of this study is to investigate a new process, named "Drapcocot" [START_REF] Coiffier-Colas | Future for high performance thermoplastic composites and aircraft helicopters[END_REF], developed by EADS, Dassault Aviation and Eurocopter, in which placement and consolidation are performed simultaneously in order to increase the productivity.

The process principle can be briefly described as follows (Fig. 1): the first layers are supposed to be already laid down. Then the tape under consideration is placed in its position using the tow guide. Two torches blow hot nitrogen in order to melt the material at the interface between the plies to weld. A roller applies a normal force on the tape in order to improve the adhesion between the plies. The process parameters that can be controlled are: torches temperature, head velocity, roller temperature and roller force.

The chosen material was the well known APC-2 pre-impregnated UD carbon-PEEK composite supplied by Cytec.

Fig. 1. Schematic principle of "Drapcocot" process 2 Influence of processing conditions on interfacial strength

The weld quality was characterized using peel tests on two layers samples. These latter differed by their processing conditions: torches temperature and laydown velocity. The applied roller force was kept constant and high enough to make the intimate contact time [START_REF] Lee | Model of the manufacturing of thermoplastics matrix[END_REF] negligible, compared to the welding time. The peel energy was plotted in Fig. 2 against temperature, for several values of the head velocity. The curves display a maximum of which the amplitude is a decreasing function of the head velocity. It appears thus that adhesion is under the dependence of two antagonists, thermally activated processes: the role of velocity puts in evidence the importance of kinetics aspects. The fact that, in the low temperature domain (typically T torches ≤ 700 C), adhesion is an increasing function of temperature and a decreasing function of velocity e.g. an increasing function of the time spent at high temperature, suggests that it results from a thermally activated diffusional process. To explain the decrease of adhesion at high temperature, the simplest hypothesis is, no doubt, that the polymer undergoes thermal ageing and that this latter affects adversely the welding quality. Other phenomena, for instance crystallization, could also affect the mechanical performance of welded joints. Their eventual influence has been studied [START_REF] Régnier | Welding of thermoplastic matrix composites: prediction of macromolecules diffusion at the interface[END_REF] and will be reported elsewhere.

According to this choice, the model elaboration must involve three steps: i) Determination of kinetic laws for chain diffusion. ii) Identification of thermal ageing processes susceptible to influence welding and determination of their kinetic parameters. iii) Thermokinetic mapping of the process in order to determine the material thermal history at every point, every time, in order to apply the above established kinetic models to each volume element.

Macromolecules diffusion

Once the intimate contact between plies is established, the interdiffusion of macromolecular chains at the interface occurs as long as the matrix is melt, i.e before crystallization quenches long range chain motions. The interfacial strength depends of the degree of chain interpenetration. Beyond a certain state, the interfacial zone becomes indistinguishable from the bulk polymer, the interface healing can be considered perfect. The (temperature dependant) time tw to reach this state is, indeed, a key quantity of the process modeling.

The chain motion in the melt is modeled by the reptation theory introduced by De Gennes [START_REF] Gennes | Scaling concepts in polymer physics[END_REF] and Doi and Edwards [START_REF] Doi | The theory of polymer dynamics[END_REF]. An important characteristic of the chain diffusion kinetics is the reptation time tr, e.g. the time for the chain to loose its initial configuration and to move over a distance equal to its end-to-end distance. In the frame of a welding experiment, tr is the time beyond which interface healing must be perfect. In order words the optimal time for welding tw must be of the order of the reptation time tr. This latter can be determined from rheological experiments.

Unfortunately, the PEEK grade used for APC-2 was not available so that rheometric measurements were made on a pure PEEK grade (G450) which was supposed to be the highest boundary of the APC-2 matrix in terms of molecular weight. The measurements were performed on an Ares Rheometric Scientific rheometer using a parallel plates geometry. Several isothermal tests were made under nitrogen in the 310-410 C temperature range, at 2% strain amplitude, with angular frequency ranging from 0.01 to 100 rad.s -1 . It can be recalled that the main transition temperature of PEEK are: glass transition temperature Tg = 140 C, melting point Tm= 330 C. In a programmed temperature decreasing experiment, crystallization begins at about 280 C.

The modulus of the complex viscosity η* was recorded against angular frequency ω. The transition between the newtonian regime, corresponding to the horizontal asymptote, and the rheo-thinning regime, corresponding to the oblique asymptote, is characterized by an angular frequency of the order of 1.1 rad.s-1 and thus the characteristic time λ = ω -1 of the order of 0.9s (at 360 C). At every temperature, this characteristic time, called also terminal relaxation time, and the newtonian viscosity η 0 were determined using the Carreau law [START_REF] Plummer | Nonisothermal fusion bonding in semicrystalline thermoplastics[END_REF]:

( ) [ ] ( ) α α λω η η / 1 0 1 * - + = m (1)
It has been shown that λ obeys Arrhenius law (Fig. 3), that allows to extrapolate its value at higher temperatures, inaccessible to experiments:

RT E λ λ λ exp 0 = (2)
with E λ equal to 59.7 KJ.mol -1. In the following, it will be considered that λ = t w , i.e. that the optimal welding time is the terminal relaxation time of the polymer. According to the literature [START_REF] Day | Kinetic study of the thermal decomposition of poly(aryl-ether-etherketone) (PEEK) in nitrogen[END_REF][START_REF] Jonas | Thermal stability and crystallization of poly(aryl ether ether ketone)[END_REF][START_REF] Cole | Fourier transform infra-red spectroscopic study of thermal degradation in poly(ether ether ketone)-carbon composites[END_REF] crosslinking is resulting apparently from radical processes and predominates during PEEK thermal ageing. It is expected to disfavour welding because it increases the melt viscosity.

To study crosslinking kinetics, it was chosen to perform isothermal experiments in the 360-460 C temperature range. The records of the real viscosity η' changes at low frequency, e.g. in the newtonian regime were done. Very small changes were observed, after 6 hours, at 360 C, whereas the viscosity increases by almost two orders of magnitude after about 2 hours at 460 C. At low conversion (far from the gel point), the number of crosslinking events per mass unit r is given by: [START_REF] Régnier | Welding of thermoplastic matrix composites: prediction of macromolecules diffusion at the interface[END_REF] where M n and M w are the number and weight average of molar mass and I P is the polydispersity index. It will be assumed that I P ~ I P0 and that the viscosity obeys the power law established for linear chains: [START_REF] Gennes | Scaling concepts in polymer physics[END_REF] where K is a temperature dependant material's property. More sophisticated relationships, taking into account the specific effect of branched architecture of crosslinked polymers on their viscosity are available, but it seems to us that the above equations are sufficient, in a first approach, to predict the trends of viscosity variations. Then, the crosslink concentration is thus given by: [START_REF] Doi | The theory of polymer dynamics[END_REF] and r was plotted against exposure time in Fig. 4. It varies quasi linearly with time, even at 460 C, in the first hour of exposure. The slope of the linear part gives the crosslinking rate K app . The crosslinking kinetics of this latter can be approximated by a zero order process:

app r K t = (6) 0 0 0 P P n n W W I I 1 1 r = - = - M M M M 3,4 0 W ' = K M η 0 0 1/ 3,4 P 0 W I ' r = 1- M ' η η ⎡ ⎤ ⎛ ⎞ ⎢ ⎥ ⎜ ⎟ ⎝ ⎠ ⎢ ⎥ ⎣ ⎦ exp a app E K A RT - ⎛ ⎞ = ⎜ ⎟ ⎝ ⎠
where the kinetic constant obeys to Arrhenius law (Fig. 5): [START_REF] Day | Kinetic study of the thermal decomposition of poly(aryl-ether-etherketone) (PEEK) in nitrogen[END_REF] with A = 3384 mol.g -1 .s -1 and E a = 168 KJ.mol -1 .

1.E-11

1.E-10

1.E-09

1.E-08 0.0013 0.0014 0.0015 0.0016

1/T (K -1 ) Kapp (mol.g -1 .s -1
) Fig. 5. Variation of the kinetic constant of crosslinking Indeed, the knowledge of the initial characteristics I P0 and M W0 , and this initial newtonian viscosity at the temperature under consideration, allow determining the viscosity η' from r at any time and temperature.

Thermokinetic mapping of the process

A time dependant bidimensional model has been build to predict temperature maps. The assumption that temperature is constant along the ply width is justified by the fact that the thermal conductivity is ten times higher in the fibre direction than in the transverse direction. The 2D model uses a finite volume formulation with an explicit scheme. The program was written with Matlab®.

The principle of modelling is schematised in Fig. 6. The heat transfer from hot nitrogen to the material is characterized by a heat transfer coefficient htor over a length ltor. Indeed, before the roller applies the pressure, both plies are not yet welded. The boundary condition is then convection with atmosphere. The interface is considered inexistent in the welded zone so that, thermal exchanges are due only to conduction. The roller imposes a local surface heat flow characterized by its heat transfer coefficient h rol and the contact length l rol . The global heat equation must thus contain two source terms relative to crystallization and melting. Both source terms can be identified from DSC characterizations [START_REF] Régnier | Welding of thermoplastic matrix composites: prediction of macromolecules diffusion at the interface[END_REF].

The experimental temperature mapping was made using small thermocouples inserted between the plies. Once a thermocouple is covered by the upper layer, one can be confident that it gives reliable data and thus that realistic boundary conditions for the torches can be determined.

The thermal simulations can predict the thermal history seen by the material. For example, we can see on Fig. 7 that the ply being welded has a higher temperature than the welded plaque.

Weldability predictions

It is now possible to predict, using the above models, the local extent of macromolecular diffusion and thermal ageing, and thus to appreciate, with pertinent criteria, the welding quality.

Macromolecular diffusion

To quantify the extent of reptation in the direction normal to the interface, it is possible to solve the diffusion equation in the nonisothermal case [START_REF] Yang | Healing of thermoplastic polymers at an interface under nonisothermal conditions[END_REF]. An alternative approach consists to consider the quantity N defined by: ( ) ( )

crystallisation contact t t dt N tr T t = ∫ (8)
where N, that depends of the composite thermal history and on the reptation time, describes the probability for a chain to cross the interface. When N = 1, this means that the processing time is equivalent to the reptation time in an isothermal experiment e.g. that the interface healing must theoretically be achieved.

N has been calculated, taking the origin of time when the interface reaches the melting point, for various sets of processing parameters. It was thus possible to plot the envelope of points where N = 1 in the Temperature-time diagram.

Thermal ageing

As for diffusion, it is possible to define a quantity R representing the extent of thermal ageing in a given point:

( ) exp crystallisation contact t a t E R A d t R T t ⎡ ⎤ - = ⎢ ⎥ ⎡ ⎤ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ ∫ (9)
Knowing that it is the new ply which is being laid down that will the more suffer thermal ageing, we will focus for this calculation at the lower surface of the layer n°13. This calculation was also implemented in the model so that it was possible to plot iso-ageing curves in the Temperature-velocity diagram.

Weldability predictions

Let us now consider a temperature-velocity diagram (Fig. 8) in which the curve corresponding to N=1 represents the "perfect welding". In fact, from our mechanical measurements, it was found that thermal ageing begins to hinder adhesion for R =10-7 mol.g-1. The envelope of the points corresponding to this value was also plotted in the diagram. The conditions for perfect welding would be thus: N ≥ 1 and r ≤ 10 -7 g.mol-1. According to Fig. 8, these conditions cannot be realized in the experimentally accessible temperature range (T ≤ 1000 C). In other words, there is no way to obtain, with the system under study, the theoretically maximum interfacial strength. This is confirmed by the results of peel tests: in the best case (T torches = 700 C, V= 10 mm.s-1), the peel energy is twice smaller than for autoclave consolidated samples.

It is however possible to define an "admissible welding domain" using a less severe endlife criterion for thermal ageing, for instance r = 10 -6 mol.g -1 , but bonding values are not acceptable.

Fig. 8 shows that technological improvements are necessary to reach ideal conditions in the high temperature/high velocity zone. Thus, the proposed model could be an interesting tool for the design of new processing systems.

The (T, t) diagram shows also that there is no acceptable solution in the domain of low temperature, low velocities, because the curves "isodiffusion" and "isoageing" diverge in the unfavourable sense.

A possible way towards better welding could consist to have a nonsymetric heating, the upper ply being less heated than the lower one, in order to shift the "isodiffusion" curve under the "isoageing" one. Here, also, the thermokinetic analysis suggests an interesting way for machine design.

Conclusion

The "Drapcocot" process is a new in-situ consolidation process for thermoplastic composites. As any polymer processing method in the molten state, it is characterized by a "thermal stability ceiling" which can be represented by a curve in the (time-temperature) space. Since consolidation results from chain diffusion through the interface, isodiffusion (e.g. iso-welding quality) conditions can also be represented by a curve in the (t, T) space. This curve can be called the "rheological floor". Using simple kinetic laws for diffusion and for thermal ageing and combining these laws with a thermal model, we have made temperature-head velocity maps in which boundaries corresponding to the "rheological floor" and to the "thermal stability ceeling" were plotted for the case of APC-2 (carbon/PEEK composites). This analysis showed that in the case under study, there is no "perfect welding window" because the "floor" is above the "ceiling". However a domain of "admissible welding" can be defined and corresponds well to experimental results. This model could be interesting to define technological objectives for the design of new welding systems.

Fig. 2 :

 2 Fig. 2: Influence of parameters on interfacial adhesion

Fig. 3 .

 3 Fig. 3. Variation of the terminal relaxation time versus temperature

Fig. 4 .

 4 Fig. 4. Evolution of r with time for different temperatures

Fig. 6 .

 6 Fig. 6. Principle of the 2D modeling of the Drapcocot process

Fig. 7 .

 7 Fig. 7. Moving temperature field.

Fig. 8 :

 8 Fig. 8: time-temperature diagram process