D Duhamel
email: denis.duhamel@enpc.fr

J.-M Mencik

Time response analysis of periodic structures via wave-based absorbing boundary conditions

Keywords: periodic structures, time response, absorbing boundary conditions, wave finite element method

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Many engineering structures are periodic, e.g., railway tracks, pipelines, ribbed plates, metamaterials, and so on. Among these are 1D periodic structures made up of cells/substructures of arbitrary shapes (e.g., 2D substructures). Usually, these structures are analyzed in the frequency domain via different methods, e.g., the semi-analytical finite element (SAFE) method developed in [START_REF] Nelson | Vibrations and waves in laminated orthotropic circular cylinders[END_REF][START_REF] Aalami | Waves in Prismatic Guides of Arbitrary Cross Section[END_REF][START_REF] Gavric | Computation of propagative waves in free rail using a finite element technique[END_REF][START_REF] Shorter | Wave propagation and damping in linear viscoelastic laminates[END_REF][START_REF] Bartoli | Modeling wave propagation in damped waveguides of arbitrary cross-section[END_REF][START_REF] Treyssède | Investigation of elastic modes propagating in multi-wire helical waveguides[END_REF][START_REF] Li | On the study of vibration of a supported railway rail using the semi-analytical finite element method[END_REF] for uniform waveguides. Within this framework, the displacement solution is separated into an analytical harmonic function, along the direction of a waveguide, and a finite element (FE) field on the cross-section. Concerning the prediction of the dynamic response of truly periodic structures with non-uniform sections under harmonic loading, the wave finite element (WFE) method is a well established numerical tool. This consists in computing wave modes (propagation constants, wave shapes) of a periodic structure from the FE model of a substructure and its related mass, damping and stiffness matrices which can be obtained from any FE software or MATLAB. Afterwards, these wave modes can be used to calculate the harmonic response of periodic structures in an efficient way, i.e., by computing small matrix systems for one substructure, or a few of them. The main steps of the method can be found in [START_REF] Mace | Finite element prediction of wave motion in structural waveguides[END_REF][START_REF] Duhamel | Finite element analysis of the vibrations of waveguides and periodic structures[END_REF][START_REF] Mencik | A wave-based model reduction technique for the description of the dynamic behavior of periodic structures involving arbitrary-shaped substructures and large-sized finite element models[END_REF][START_REF] Mencik | Multi-mode propagation and diffusion in structures through finite elements[END_REF].

Also, some of its recent extensions and applications are reported in [START_REF] Mencik | A wave finite element-based approach for the modeling of periodic structures with local perturbations[END_REF][START_REF] Mencik | New advances in the forced response computation of periodic structures using the wave finite element (WFE) method[END_REF][START_REF] Mencik | A model reduction strategy for computing the forced response of elastic waveguides using the wave finite element method[END_REF][START_REF] Mencik | A wave finite element approach for the analysis of periodic structures with cyclic symmetry in dynamic substructuring[END_REF][START_REF] Ichchou | Wave finite elements for low and mid-frequency description of coupled structures with damage[END_REF][START_REF] Hoang | Wave finite element method for waveguides and periodic structures subjected to arbitrary loads[END_REF][START_REF] Renno | On the forced response of waveguides using the wave and finite element method[END_REF][START_REF] Waki | Free and forced vibrations of a tyre using a wave/finite element approach[END_REF][START_REF] Singh | Stochastic wave finite element quadratic formulation for periodic media: 1D and 2D[END_REF][START_REF] Fan | Model reduction schemes for the wave and finite element method using the free modes of a unit cell[END_REF]. So far, these studies have been conducted in the frequency domain in which the WFE method has been originally formulated.

The time response analysis of periodic structures under time-dependent excitations -e.g., fast loadings like shocks or blast -is not so well reported. For infinite structures, the main difficulty is to solve a time differential equation on a given finite domain, and therefore, be able to express appropriate (time-dependent) absorbing/infinite conditions at the boundary of this domain. Timedependent absorbing boundary conditions (BCs) are usually formulated for homogeneous domains like acoustic or elastic media. For instance, absorbing BCs involving differential operators of different orders have been proposed in [START_REF] Engquist | Absorbing boundary conditions for numerical simulation of waves[END_REF][START_REF] Engquist | Absorbing boundary conditions for numerical simulation of waves[END_REF][START_REF] Clayton | Absorbing boundary conditions for acoustic and elastic wave equations[END_REF][START_REF] Reynolds | Boundary conditions for the numerical solution of wave propagation problems[END_REF] and in [START_REF] Bayliss | Radiation boundary conditions for wave-like equations[END_REF] using sequences of local non-reflecting BCs in spherical and cylindrical coordinates. This approach has been modified to describe absorbing BCs for waves at some discrete angles from the surface normal [START_REF] Higdon | Absorbing boundary conditions for difference approximations to the multidimensional wave equation[END_REF][START_REF] Higdon | Numerical absorbing boundary conditions for the wave equation[END_REF]. More efficient BCs can be obtained by adding supplementary variables to model an exterior surface, see [START_REF] Collino | High-order absorbing boundary conditions for wave propagation models. straight line boundary and corner cases[END_REF][START_REF] Hagstrom | A formulation of asymptotic and exact boundary conditions using local operators[END_REF][START_REF] Lee | Absorbing boundary condition for scalar-wave propagation problems in infinite media based on a root-finding algorithm[END_REF]. Another possibility is the perfectly matched layer (PML) method proposed in [START_REF] Berenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF][START_REF] Berenger | Three-dimensional perfectly matched layer for the absorption of electromagnetic waves[END_REF][START_REF] Collino | Optimizing the perfectly matched layer[END_REF][START_REF] Asvadurov | On optimal finite-difference approximation of PML[END_REF] which consists in surrounding a computational domain with one or several layers of elements in which the wave equation is analytically continued into complex coordinates. Contrary to the previous approaches which are easy to implement but of limited accuracy, the PML method does not require analytical expressions of the absorbing BCs, but instead, one or several additional FE layers with increasing and well selected absorption properties.

Although more complex than the homogeneous case, the time domain analysis of periodic media has been carried out in various ways, mostly for electromagnetic applications. In [START_REF] Baida | Finite Difference Time Domain Method For Grating Structures[END_REF][START_REF] Li | Time-Domain Finite Element Methods for Maxwell's Equations in Metamaterials[END_REF], PMLs have been used as absorbing BCs. In [START_REF] Deinega | Long-time behavior of PML absorbing boundaries for layered periodic structures[END_REF], it has been found, however, that the PML method can be prone to divergence issue for analyzing layered periodic electromagnetic media with evanescent waves. In [START_REF] Zheng | Implementation of Mur's absorbing boundaries with periodic structures to speed up the design process using finite-difference time-domain method[END_REF], approximate time domain absorbing BCs have been developed for the analysis of the electromagnetic wave propagation in periodic media. Also, in [START_REF] Gao | Analysis of time-domain scattering by periodic structures[END_REF], the time domain scattering of waves in infinite media by diffraction gratings has been conducted in a theoretical way using a Dirichlet-to-Neumann (DtN) map, without numerical results. It should be remarked that the DtN map approach involves a non-local operator which, as such, yields large sized fully-populated impedance matrices. A similar approach has been used in [START_REF] Marrocco | Transient radiation by periodic structures: Accuracy of the (time domain-floquet wave)-fdtd algorithm[END_REF][START_REF] Sirenko | Exact absorbing boundary conditions for periodic threedimensional structures: Derivation and implementation in discontinuous Galerkin time-domain method[END_REF]. In these works, the transient electromagnetic wave interaction with diffraction gratings has been analyzed via discontinuous Galerkin methods and exact non-local absorbing BCs for 3D periodic media. The procedure requires convolution integrals in space and time on a boundary.

Although well suited to describe absorbing BCs for homogeneous media, or simple periodic media, all of these approaches appear difficult to apply to complex systems like periodic elastic structures with arbitrary-shaped cells. Regarding the PML method, the issue mainly concerns the determination of the absorption properties of the layers to absorb waves emanating from an heterogeneous medium. Also, most of these approaches focus on the analysis of non-dispersive media having well defined wave propagation velocities. For dispersive media like beams [START_REF] Ruge | A comparison of infinite Timoshenko and Euler-Bernoulli beam models on Winkler foundation in the frequency-and time-domain[END_REF], or multi-layered systems [START_REF] Zhao | Stable high-order absorbing boundary condition based on new continued fraction for scalar wave propagation in unbounded multilayer media[END_REF], absorbing BCs can be formulated in terms of boundary operators involving complicated functions of the frequency -e.g., square roots of the frequency -which, as such, cannot be converted to simple functions of time after inverse Fourier transforms. Also, these approaches appear difficult to implement given that a substantial part of the time history of the solutions must be kept, which means computing time integrals at each time step as discussed in [START_REF] Weber | Rational transmitting boundaries for time-domain analysis of dam-reservoir interaction[END_REF]. Another strategy for describing absorbing BCs for beams involves changing the geometry at the ends by considering tapered beams [START_REF] Yang | A non-reflecting boundary for use in a finite element beam model of a railway track[END_REF], which, however, seems difficult to apply to more complicated structures.

Considering the potential of the WFE method for predicting the dynamic behavior of periodic structures in the frequency domain, in this work it is intended to extend this method to the time response analysis of infinite periodic structures. This, however, cannot be carried out straightforwardly. The issue is linked to the wave-based matrix equations involved in the WFE method which represent complicated functions of the frequency (i.e., iω) and which do not generally possess an analytical expression. This makes any time domain transform hard to perform. The goal of the present paper is to provide a possible answer to this question by targeting the problem of infinite 1D periodic structures subject to localized time-dependent excitations. These are structures which are frequently encountered in mechanical engineering, e.g., railway tracks or slender periodic beam-like structures. In order to stay as close as possible to a classical FE modeling, a finite FE mesh of a structure truncated at some finite extent from the excitation sources is considered (see Fig. 1). This allows analyzing the evanescent near-field effect in the vicinity of the excitations, but also, analyzing structures which are not necessarily strictly periodic in that region and which are modeled by finite elements. In this case, the coupling effect with the rest of the infinite structure -i.e., two semi-infinite periodic structures that would expand to the left and right directions at the left and right boundaries of the truncated structure -are taken into account via absorbing BCs. Although easy established in the frequency domain via the WFE method -i.e., by retaining only the waves which propagate outwards from the structure to the left or right semi-infinite structure -these absorbing BCs cannot be easily converted to the time domain for the reasons mentioned earlier.

In this paper, an original approach well suited for modeling periodic structures with absorbing BCs in the time domain is proposed. Within the FE framework, the periodic structures are usually modeled by means of vectors of displacements, velocities and accelerations. The key idea of the present work is to consider vectors of supplementary variables for modeling absorbing BCs by means of classical second-order time differential equations. In this sense, a global second-order time differential equation can be proposed for modeling a periodic structure with absorbing BCs, which can be easily integrated via the Newmark algorithm. The advantages of the proposed approach compared to the aforementioned ones are: (i) it is well adapted to the analysis of periodic structures like those made up of arbitrary-shaped substructures (e.g., 2D substructures); (ii) it can be easily implemented in an FE model and yields second-order time differential equations which are standard and easy to solve. The key steps of the proposed approach are summarized as follows. First, the absorbing BCs are expressed in the frequency domain via impedance matrices and the WFE method. Then, the impedance matrices are decomposed via rational approximations using rational functions with computed poles and matrices of residues. These rational functions can be reformulated in terms of mass, damping and stiffness matrices by using vectors of supplementary variables at the left and right boundaries of the structure in a similar way to the works in [START_REF] Collino | High-order absorbing boundary conditions for wave propagation models. straight line boundary and corner cases[END_REF][START_REF] Hagstrom | A formulation of asymptotic and exact boundary conditions using local operators[END_REF][START_REF] Lee | Absorbing boundary condition for scalar-wave propagation problems in infinite media based on a root-finding algorithm[END_REF]. Finally, by assembling these new matrices together with the FE mass, damping and stiffness of the structure, this yields a classical dynamic matrix system in the time domain which can be solved by the Newmark algorithm.

The rest of the paper is organized as follows. In Sec. 2, the main steps of the WFE method are recalled. The FE modeling of a periodic structure of finite length -i.e., with a finite number of substructures -subject to time-dependent excitations and absorbing BCs at its left and right ends is presented. The related impedance matrices, in the frequency domain, are formulated. Also, the strategy to express the absorbing BCs in the time domain is proposed. The Newmark integration scheme used to solve the equation of motion of the structure is finally detailed. In Sec. 3, the proposed approach is numerically validated to predict the time response of an infinite Euler-Bernoulli beam lying on an elastic foundation and subject to a point harmonic force, where analytical solutions exist. In Sec. 4, the proposed approach is used to describe the time response of two periodic structures involving 2D substructures. Point forces of harmonic type or Ricker wavelet type (shock) are analyzed. Numerical comparisons with an equivalent infinite full FE model are brought. Concluding remarks are finally brought in Sec. 5

Theory

General problem

The present paper investigates the dynamic response of infinite 1D periodic structures subject to time-dependent loadings (e.g., shocks), and the modeling of such infinite structures by means of structures of finite length and appropriate absorbing BCs for the waves escaping from excitations. For instance, a schematic of a periodic structure made up of N = 10 identical substructures subject to time-dependent forces and absorbing BCs is shown in Fig. 1. The substructures under concern can be of arbitrary shape -e.g., 2D substructures as shown in Fig. 2 -and are supposed to be linear, elastic, isotropic and damped (viscous damping). Also, the excitations are supposed to be confined to a part of the structure encompassing those substructures and bounded by two left and right boundaries S L and S R , as shown in Fig. 1, where absorbing BCs are considered. From a physical point of view, these absorbing BCs mean that the energy is escaping towards infinity and that no wave comes from infinity. The purpose of the present paper is to propose a FE model which is well suited for time response analysis of such truncated periodic structures via appropriate "time-based" absorbing conditions on S L and S R to describe semi-infinite structures at both sides.

WFE method

The wave propagation analysis in 1D periodic structures can be conducted with the WFE method [START_REF] Duhamel | Finite element analysis of the vibrations of waveguides and periodic structures[END_REF][START_REF] Mencik | Multi-mode propagation and diffusion in structures through finite elements[END_REF][START_REF] Zhong | On the direct solution of wave propagation for repetitive structures[END_REF]. The basics of the method are recalled hereafter. Let us consider infinite structures under harmonic disturbance e iωt which are built from identical substructures as shown in Fig. 2. Those substructures are supposed to share the same FE mesh, and are modeled by means of identical mass, damping and stiffness matrices M s , C s and K s . The related dynamic equilibrium equation is given by:

D s q s = F s , (1)
where q s and F s refer to the displacement vector and the force vector (respectively), and where D s is the dynamic stiffness matrix of the substructures (similar for all the substructures) expressed by

D s = -ω 2 M s + iωC s + K s .
The FE mesh of a substructure is shown in Fig. 2, and involves left (L) and right (R) boundaries which are described with the same number n of degrees of freedom (DOFs). Rearranging Eq. ( 1) yields the following transfer matrix relation between the right and left boundaries of the substructure:

u s R = S s u s L , (2)
where u s R and u s L are 2n × 1 state vectors expressed by:

u s R = q s R F s R , u s L = q s L -F s L . (3)
Also, S s is a symplectic 2n × 2n matrix (also called transfer matrix) expressed by: where D s * refers to the dynamic stiffness matrix of the substructure condensed on the left and right boundaries [START_REF] Mencik | New advances in the forced response computation of periodic structures using the wave finite element (WFE) method[END_REF].

S s = -(D s * LR ) -1 D s * LL -(D s * LR ) -1 D s * RL -D s * RR (D s * LR ) -1 D s * LL -D s * RR (D s * LR ) -1 , (4)
The eigenvalues and eigenvectors of the transfer matrix S s occur in pairs as (µ j , φ j ) and (µ j = 1/µ j , φ j ) with |µ j | < 1 (see [START_REF] Mencik | A wave-based model reduction technique for the description of the dynamic behavior of periodic structures involving arbitrary-shaped substructures and large-sized finite element models[END_REF] for further details about the computation of the eigensolutions of S s ). Let us denote by d the substructure length, i.e., along the main direction of the periodic structure. According to the Bloch's theorem [START_REF] Willcox | Balanced model reduction via the proper orthogonal decomposition[END_REF], the eigenvalues of S s have the meaning of wave parameters which are given by µ j = e -ik j d and µ j = e ik j d (k j being the wavenumbers); also, the eigenvectors φ j (resp. φ j ) have the meaning of wave shapes, for the waves traveling to the right and left directions (respectively) of the periodic structure. Those vectors of wave shapes are of size 2n × 1 and are usually partitioned as follows:

φ j = φ qj φ Fj , φ j = φ qj φ Fj , (5)
where φ qj and φ qj (resp. φ Fj and φ Fj ) are n × 1 vectors involving displacement (resp. force) components. The related n × n matrices of wave shapes -namely, Φ q , Φ q , Φ F and Φ F -are given by:

Φ q = [φ q1 • • • φ qn ] , Φ q = [φ q1 • • • φ qn ] , Φ F = [φ F1 • • • φ Fn ] , Φ F = [φ F1 • • • φ Fn ]. (6)

Absorbing boundary conditions in the frequency domain

Consider a periodic structure involving a finite number N of substructures which is enclosed between two left and right boundaries S L and S R where absorbing BCs are considered (see Fig. 1). Such BCs are used to describe the coupling effect between the periodic structure and two semi-infinite periodic structures, not shown here, that would expand to the left and right directions. The related impedance matrices for the left and right boundaries S L and S R -namely, Z L and Z R -can be defined, in the frequency domain, as follows:

F L = Z L q L , F R = Z R q R , (7)
where q L and q R (resp. F L and F R ) are the displacement (resp. force) vectors (size n × 1) for the periodic structure on S L and S R . Following the WFE procedure and expanding those vectors on the basis of wave shapes, this yields [START_REF] Mencik | New advances in the forced response computation of periodic structures using the wave finite element (WFE) method[END_REF]:

q L = Φ q Q L , q R = Φ q Q R , (8)
and

F L = -Φ F Q L , F R = Φ F Q R , (9)
where Q L (resp. Q R ) is the vector of wave amplitudes, at position S L (resp. S R ), for the waves traveling to the left (resp. right) direction of the structure. The fact that left-going (resp. right-going) waves are only considered at S L (resp. S R ) results from the absorbing BC, i.e., the fact that no wave comes from infinity. Eq. ( 8) yields Q L = (Φ q ) -1 q L and Q R = (Φ q ) -1 q R and, therefore, the following expressions for the impedance matrices:

Z L = -Φ F (Φ q ) -1 , Z R = Φ F (Φ q ) -1 . ( 10
)
Note that, withing the WFE framework, Φ q , Φ q , Φ F and Φ F , and therefore Z L and Z R , are matrixvalued functions of the frequency/pulsation ω (or iω). The thing is, the dependency of these matrices on iω is not straightforward, i.e., complicated functions which do not generally possess an analytical expression. The issue is linked to the numerical steps involved in the WFE method (derivation of a transfer matrix, resolution of an eigenproblem, among others) which makes these matrix-valued functions hard to figure out. To solve this problem, it is proposed to recast Eq. ( 7) as two alternative matrix equations involving new impedance matrices which are simply described in terms of secondorder polynomials of iω. This provides a simple way to express absorbing BCs, for the periodic structure, in the time domain. The strategy is explained as follows.

Let us decompose the impedance matrices Z L and Z R via rational approximations (following similar ideas as in [START_REF] Moon | Rational function interpolation of electromagnetic transfer functions of high-speed interconnect systems from discrete time-domain and frequency-domain data[END_REF][START_REF] Cottereau | Modèle dynamique équivalent de matrices d'impédance de fondation[END_REF][START_REF] Gustavsen | Computer code for rational approximation of frequency dependent admittance matrices[END_REF] for approximating transfer functions):

Z L = P k=1 R Lk iω -p Lk + K L , Z R = P k=1 R Rk iω -p Rk + K R , (11)
where p Lk and p Rk denote poles (k = 1, . . . , P ), and R Lk and R Rk denote matrices of residues. Some of these usually appear in conjugate pairs, i.e., (p Lk , p Lk ) and (p Rk , p Rk ), and (R Lk , R Lk ) and (R Rk , R Rk ).

In this case, it can be easily proven that:

R Lk iω -p Lk + R Lk iω -p Lk = 2 iω {R Lk } -{p Lk R Lk } -ω 2 -2iω {p Lk } + |p Lk | 2 , R Rk iω -p Rk + R Rk iω -p Rk = 2 iω {R Rk } -{p Rk R Rk } -ω 2 -2iω {p Rk } + |p Rk | 2 . ( 12
)
As a result, Eq. ( 11) leads to:

Z L = Q k=1 2 iω {R L(2k) } -{p L(2k) R L(2k) } -ω 2 -2iω {p L(2k) } + |p L(2k) | 2 + P k=2Q+1 R Lk iω -p Lk + K L , Z R = Q k=1 2 iω {R R(2k) } -{p R(2k) R R(2k) } -ω 2 -2iω {p R(2k) } + |p R(2k) | 2 + P k=2Q+1 R Rk iω -p Rk + K R , (13)
where 2Q < P . To remove the denominator terms in Eq. ( 13), let us introduce n × 1 vectors of supplementary variables X Lk and X Rk , and let us rewrite Eq. ( 7) by means of Eq. ( 13) as follows:

F L = Q k=1 2 iω {R L(2k) } -{p L(2k) R L(2k) } X Lk + P k=2Q+1 R Lk (iω)X L(k-Q) + K L q L , F R = Q k=1 2 iω {R R(2k) } -{p R(2k) R R(2k) } X Rk + P k=2Q+1 R Rk (iω)X R(k-Q) + K R q R , (14)
where:

-ω 2 -2iω {p L(2k) } + |p L(2k) | 2 X Lk = q L for k = 1, . . . , Q, -ω 2 -2iω {p R(2k) } + |p R(2k) | 2 X Rk = q R for k = 1, . . . , Q, (-ω 2 -iωp Lk )X L(k-Q) = q L for k = (2Q + 1), . . . , P, (-ω 2 -iωp Rk )X R(k-Q) = q R for k = (2Q + 1), . . . , P. (15)
Finally, let us introduce the following (P -Q)n × 1 vectors X L and X R defined by:

X L =            X L1 . . . X LQ X L(Q+1) . . . X L(P -Q)            , X R =            X R1 . . . X RQ X R(Q+1) . . . X R(P -Q)            . (16)
Then, Eqs. ( 14) and ( 15) yield:

K L(qq) K L(qX) K L(Xq) K L(XX) q L X L = F L 0 , K R(qq) K R(qX) K R(Xq) K R(XX) q R X R = F R 0 , (17)
where

K L(qq) = K L , K L(qX) = 2iω R L(2) • • • R L(2Q) -2 p L(2) R L(2) • • • p L(2Q) R L(2Q) iω R L(2Q+1) • • • R LP , K L(Xq) = -1 Q×1 ⊗ I n -1 (P -2Q)×1 ⊗ I n , K L(XX) =   blkdiag -ω 2 -2iω {p L(2k) } + |p L(2k) | 2 I n Q k=1 0 0 blkdiag (-ω 2 -iωp Lk )I n P k=2Q+1   , (18)
and

K R(qq) = K R , K R(qX) = 2iω R R(2) • • • R R(2Q) -2 p R(2) R R(2) • • • p R(2Q) R R(2Q) iω R R(2Q+1) • • • R RP , K R(Xq) = -1 Q×1 ⊗ I n -1 (P -2Q)×1 ⊗ I n , K R(XX) =   blkdiag -ω 2 -2iω {p R(2k) } + |p R(2k) | 2 I n Q k=1 0 0 blkdiag (-ω 2 -iωp Rk )I n P k=2Q+1   , (19)
where ⊗ denotes the Kronecker product. The block components of the matrices occurring in Eqs. ( 18) and ( 19) represent polynomials of iω of orders 0, 1 or 2, which as such can be simply and quickly convert to the time domain (see hereafter).

Absorbing boundary conditions in the time domain

By separating the terms of identical powers of iω in Eq. ( 17), and by invoking the classical timefrequency transforms q(ω) → q(t), iωq → q, -ω 2 q → q and X(ω) → X(t), iωX → Ẋ, -ω 2 X → Ẍ (where dot and double-dot notations mean single and double time derivatives, respectively), this yields:

M L qL ẌL +C L qL ẊL +K L q L X L = F L 0 , M R qR ẌR +C R qR ẊR +K R q R X R = F R 0 , (20)
where:

M L =     0 0 0 0 blkdiag {I n } Q k=1 0 0 0 blkdiag {I n } P k=2Q+1     , C L =     0 2 R L(2) • • • R L(2Q) R L(2Q+1) • • • R LP 0 blkdiag -2 {p L(2k) }I n Q k=1 0 0 0 blkdiag {-p Lk I n } P k=2Q+1     , K L =     K L -2 p L(2) R L(2) • • • p L(2Q) R L(2Q) 0 -1 Q×1 ⊗ I n blkdiag |p L(2k) | 2 I n Q k=1 0 -1 (P -2Q)×1 ⊗ I n 0 0     , (21)
and

M R =     0 0 0 0 blkdiag {I n } Q k=1 0 0 0 blkdiag {I n } P k=2Q+1     , C R =     0 2 R R(2) • • • R R(2Q) R R(2Q+1) • • • R RP 0 blkdiag -2 {p R(2k) }I n Q k=1 0 0 0 blkdiag {-p Rk I n } P k=2Q+1     , K R =     K R -2 p R(2) R R(2) • • • p R(2Q) R R(2Q) 0 -1 Q×1 ⊗ I n blkdiag |p R(2k) | 2 I n Q k=1 0 -1 (P -2Q)×1 ⊗ I n 0 0     . (22)
Let us denote by M, C and K the mass, damping and stiffness matrices of the periodic structure (N substructures), and let us write the related equation of motion as follows:

Mq + C q + Kq = F, (23)
where q = q(t) and F = F(t) are the displacement and force vectors, respectively. In this case, the force vector is expressed by

F = [F T I F T L F T R ]
where F L and F R refer to the force vectors on S L and S R (absorbing BCs), and F I refers to the force vector for the internal DOFs (I) of the structure. Also, the displacement vector is expressed by q = [q T I q T L q T R ] where q I is the displacement vector for the internal DOFs. Then, let us rewrite Eq. ( 23) as follows:

  M II M IL M IR M LI M LL M LR M RI M RL M RR     qI qL qR   +   C II C IL C IR C LI C LL C LR C RI C RL C RR     qI qL qR   +   K II K IL K IR K LI K LL K LR K RI K RL K RR     q I q L q R   =   F I F L F R   . (24)
By considering the absorbing BCs (Eq. ( 20)), this yields:

M tot       qI qL qR ẌL ẌR       + C tot       qI qL qR ẊL ẊR       + K tot       q I q L q R X L X R       =       F I 0 0 0 0       . ( 25
)
Eq. ( 25) represents a second-order differential matrix equation for the displacement vector q and the vector of supplementary variables X = [X T L X T R ] T . This indeed represents a classical dynamic equation, in the time domain, of a structure with absorbing BCs and subject to an input force vector F I = F I (t), with the only modification that supplementary DOFs are added at the boundaries. Therefore, this equation can be solved in a standard way via a time integration numerical scheme (e.g., Newmark scheme).

In [START_REF] Reynolds | Boundary conditions for the numerical solution of wave propagation problems[END_REF], the matrices M tot , C tot and K tot are given by:

M tot =       M II M IL M IR 0 0 M LI M LL M LR 0 0 M RI M RL M RR 0 0 0 0 0 M L(XX) 0 0 0 0 0 M R(XX))       , C tot =       C II C IL C IR 0 0 C LI C LL C LR -C L(qX) 0 C RI C RL C RR 0 -C R(qX) 0 0 0 C L(XX) 0 0 0 0 0 C R(XX))       , K tot =       K II K IL K IR 0 0 K LI K LL -K L(qq) K LR -K L(qX) 0 K RI K RL K RR -K R(qq) 0 -K R(qX) 0 K L(Xq) 0 K L(XX) 0 0 0 K R(Xq) 0 K R(XX)       , (26)
where:

M L(XX) = blkdiag {I n } Q k=1 0 0 blkdiag {I n } P k=2Q+1 , C L(qX) = 2 R L(2) • • • R L(2Q) R L(2Q+1) • • • R LN , C L(XX) =   blkdiag -2 {p L(2k) }I n Q k=1 0 0 blkdiag {-p Lk I n } P k=2Q+1   , K L(qq) = K L , K L(qX) = -2 p L(2) R L(2) • • • p L(2Q) R L(2Q) 0 , K L(Xq) = -1 Q×1 ⊗ I n -1 (P -2Q)×1 ⊗ I n , K L(XX) = blkdiag |p L(2k) | 2 I n Q k=1 0 0 0 , (27)
and

M R(XX) = blkdiag {I n } Q k=1 0 0 blkdiag {I n } P k=2Q+1 , C R(qX) = 2 R R(2) • • • R R(2Q) R R(2Q+1) • • • R RN , C R(XX) =   blkdiag -2 {p R(2k) }I n Q k=1 0 0 blkdiag {-p Rk I n } P k=2Q+1   , K R(qq) = K R , K R(qX) = -2 p R(2) R R(2) • • • p R(2Q) R R(2Q) 0 , K R(Xq) = -1 Q×1 ⊗ I n -1 (P -2Q)×1 ⊗ I n , K R(XX) = blkdiag |p R(2k) | 2 I n Q k=1 0 0 0 . (28)

Solution by implicit methods

Let us rewrite Eq. ( 25) in concise form as follows:

M tot ÿ + C tot ẏ + K tot y = F tot , (29)
where

y = y(t) =       q I (t) q L (t) q R (t) X L (t) X R (t)       , F tot = F tot (t) =       F I (t) 0 0 0 0       . ( 30
)
Eq. ( 29) can be solved by means of the classical Newmark algorithm. Hence, consider a time step of ∆t and define the following vector y n = y(n∆t) which represents the vector y expressed at a discrete time n∆t for n = 0, 1, 2, . . .. Following the Newmark scheme, one has:

y n+1 = ỹn+1 + β∆t 2 ÿn+1 , ẏn+1 = ỹn+1 + γ∆tÿ n+1 , (31)
where

ỹn+1 = y n + ∆t ẏn + ( 1 2 -β)∆t 2 ÿn , ỹn+1 = ẏn + (1 -γ)∆tÿ n , (32)
and

ÿn+1 = 1 β∆t 2 (y n+1 -ỹn+1 ). ( 33
)
Also, in Eq. ( 31), γ and β are the so-called Newmark parameters, e.g., γ = 1/2 and β = 1/4 (average constant acceleration rule). The residue of relation ( 29) is given by:

r(y) = M tot ÿ + C tot ẏ + K tot y -F tot . (34)
From the knowledge of the solution at n∆t, one has to find the solution at the subsequent time (n + 1)∆t such that:

r(y n+1 ) = M tot ÿn+1 + C tot ẏn+1 + K tot y n+1 -F tot (t n+1 ) = 1 β∆t 2 M tot (y n+1 -ỹn+1 ) + C tot ( ỹn+1 + γ β∆t (y n+1 -ỹn+1 )) +K tot y n+1 -F tot (t n+1 ) = 0. (35)
The solution at (n + 1)∆t follows by solving a linear system:

My n+1 = F n+1 , (36)
where:

M = 1 β∆t 2 M tot + γ β∆t C tot + K tot , F n+1 = 1 β∆t 2 M tot ỹn+1 -C tot ( ỹn+1 - γ β∆t ỹn+1 ) + F tot (t n+1 ). ( 37
)

Euler-Bernoulli beam on an elastic foundation

For validation purposes, let us analyze the dynamic response of an infinite Euler-Bernoulli beam lying on an elastic foundation (lineic stiffness k f ) as shown in Fig. 3 and subject to some "localized" lineic forces f (x, t) for which analytical solutions are available. In this case, the governing equation of motion of the beam is given by:

ρSv + EI ∂ 4 v ∂x 4 + ξ ∂ 4 v ∂x 4 + k f v = f (x, t), (38)
where v = v(x, t) represents the transverse displacement, ρ is the density, S is the cross-sectional area, E is the Young's modulus, I is the inertia moment, and ξ is a damping parameter. For harmonic disturbance of the form f (x)e iωt , Eq. ( 38) leads to:

(-ρSω 2 + k f )v + EI(1 + iωξ) ∂ 4 v ∂x 4 = f (x), (39)
where v = v(x). This simple example of a beam on an elastic foundation aims at illustrating the methodology for obtaining absorbing BCs in time. There exist plenty of works about the dynamic analysis of beams with or without elastic foundations. However, even in this simple case, obtaining time absorbing BCs is not obvious. In fact, only a few works have been conducted on this topic. These mostly address approximate (not exact) expressions of absorbing BCs for Euler-Bernoulli beams. For instance, approximate expressions of absorbing BCs in time for Euler-Bernoulli beams have been proposed in [START_REF] Delavaud | Rolling noise simulation: A new approach in time domain[END_REF]. In this case, the off-diagonal terms of the boundary impedance matrices are ignored, i.e., the coupling terms between the displacement and the moment and between the force and the rotation. Also, approximate expressions of absorbing BCs (in time) for beams under moving loads in convected coordinates have been proposed in [START_REF] Andersen | Finite element modelling of infinite Euler beams on Kelvin foundations exposed to moving loads in convected co-ordinates[END_REF]. In this case, absorbing BCs are obtained via the use of added stiffness and damping, tuned at a specific frequency. In [START_REF] Feng | Matching boundary conditions for the Euler-Bernoulli beam[END_REF], absorbing BCs have been derived for beams discretized by finite difference schemes. In [START_REF] Montseny | Fractional integro-differential boundary control of the Euler-Bernoulli beam[END_REF], integrodifferential operators have been used for boundary dynamic feedback control to get absorbing BCs. In [START_REF] Tang | Artificial boundary conditions for Euler-Bernoulli beam equation[END_REF], artificial absorbing BCs in time have been proposed, while in [START_REF] Arbabi | Propagation of waves in infinite beams: PML approach[END_REF], the PML method has been used. Finally, in [START_REF] Ruge | A comparison of infinite Timoshenko and Euler-Bernoulli beam models on Winkler foundation in the frequency-and time-domain[END_REF], absorbing BCs for beams lying on a Winkler foundation have been obtained via the use of rational approximations for low frequencies, and asymptotic analysis and fractional derivatives for high frequencies. In this case, the proposed rational approximations appear to be roughly similar to those investigated in the present paper, except that in the present paper, the proposed approach is intended to work for complex periodic structures (not only for beams) and is not restricted to any low frequency assumptions other than those raised by the FE modeling.

Wave shapes: analytical expressions

For this simple case, there exist analytical expressions of the matrices of wave shapes Φ q , Φ q , Φ F and Φ F , see Eq. ( 6). Here, the number of boundary DOFs is n = 2. The determination of these matrices involves solving Eq. ( 39) with f (x) = 0 (free wave propagation), which yields v = Q 1 e -ikx + Q 2 e -kx + Q 1 e ikx + Q 2 e kx where:

k = ρSω 2 -k f EI(1 + iωξ) 1/4 . ( 40
)
Hence, by expressing the transverse displacement v together with the rotation θ = ∂v/∂x, this yields:

v θ = Φ q e -ikx 0 0 e -kx Q 1 Q 2 + Φ q e ikx 0 0 e kx Q 1 Q 2 , (41)
where:

Φ q = 1 1 -ik -k , Φ q = 1 1 ik k . ( 42
)
Also, by expressing the shearing force V = -EI(∂ 3 v/∂x 3 ) and the bending moment M = EI(∂ 2 v/∂x 2 ), this yields:

V M = Φ F e -ikx 0 0 e -kx Q 1 Q 2 + Φ F e ikx 0 0 e kx Q 1 Q 2 , ( 43
)
where:

Φ F = EIk 2 -ik k -1 1 , Φ F = EIk 2 ik -k -1 1 . ( 44
)

Impedance matrices: analytical expressions

Analytical expressions of the impedance matrices follow as Z L = -Φ F (Φ q ) -1 and Z R = Φ F (Φ q ) -1 , see Eq. ( 10). In the present case, the matrix inverses (Φ q ) -1 and (Φ q ) -1 are expressed by:

(Φ q ) -1 = 1 1 -i 1 -1 k -i 1 k , (Φ q ) -1 = 1 1 -i 1 1 k -i -1 k . ( 45
)
Therefore, the impedance matrices are written as:

Z L = -Φ F (Φ q ) -1 = - EIk 2 1 -i 2ik -(1 + i) -(1 + i) 2 k , ( 46
)
and

Z R = Φ F (Φ q ) -1 = - EIk 2 1 -i 2ik 1 + i 1 + i 2 k . ( 47
)
These are analytical expressions of the impedance matrices that could also have been obtained numerically by means of the WFE method. For more complex structures (see next section), only those obtained with the WFE method will be available.

Numerical results

Consider an infinite beam of rectangular cross-section having the following parameters: height h = 0.001 m, width b = 0.01 m (i.e., S = bh and I = bh 3 /12), Young's modulus E = 2.2 × 10 11 Pa, density ρ = 7800 kg/m 3 , damping parameter ξ = 0.001 s. Also, the lineic stiffness of the elastic foundation is k f = 1 N/m 2 . The system is at rest at time t = 0 -i.e., v(x, 0) = 0 and v(x, 0) = 0and, for t ≥ 0, it is excited by a point harmonic force of frequency f 0 = 5 Hz (at x = 0):

f (x, t) = cos(2πf 0 t)δ(x) for t ≥ 0. ( 48
)
The harmonic solution of this equation is found by solving Eq. ( 39) analytically, i.e.,

v(x, t) = - i 4k 3 0 E 0 I e -ik 0 |x| - 1 4k 3 0 E 0 I e -k 0 |x| e iω 0 t , (49)
which should match the solution of the equation of motion [START_REF] Deinega | Long-time behavior of PML absorbing boundaries for layered periodic structures[END_REF] after a certain time, i.e., once the transient effects become negligible. Here, the wavenumber k 0 is given by Eq. ( 40) when ω = ω 0 ; also, E 0 = E(1 + iξω 0 ). On the other hand, the time-space response of the beam can be computed by considering the proposed approach. Within this framework, a beam of finite length L -i.e., x ∈ [-L/2, L/2] where, for instance, L = 5 m -excited at x = 0 (Eq. ( 48)) is considered as shown in Fig. 3. Here, the system beam-foundation is modeled by means of 500 identical substructures that represent identical two-node Hermitian beam elements of length d = 0.01 m. For f 0 = 5 Hz, this gives a wavelength of λ = 2π/ {k 0 } = 0.39 m and, therefore, 39 elements per wavelength. Additional numerical tests (not presented here) with d = 0.1 m and 4 elements per wavelength also give correct, although less accurate, results. As for the left and right boundaries x = -L/2 (S L ) and x = L/2 (S R ), they are described via absorbing BCs, see Secs. 2.3 and 2.4. The rational approximations of the impedance matrices Z L and Z R with P poles/residues (see Eqs. ( 11)) are computed with the MATLAB's rationalfit function. Especially, the components of the impedance matrix Z R issued from different orders P are plotted and compared to the analytical expressions given by Eq. ( 47) as shown in Fig. 4. Results show stable curves which closely match the analytical solutions for orders larger than P = 7. Besides, the relative errors between the proposed and analytical curves are given in Table 1 for orders up to P = 12 which is the value proposed by MATLAB and which will be used hereafter. In this case, the relative error is smaller than -80 dB, which appears to be enough to meet good convergence. The relative error e lm on the lm component of the impedance matrix is defined as:

e lm = ( i |Z rat lm (ω i ) -Z lm (ω i )| 2 ) 1/2 ( i |Z lm (ω i )| 2 ) 1/2 , (50)
with Z lm (ω i ) the lm component of the impedance matrix at the discrete circular frequency ω i , and Z rat lm (ω i ) its rational approximation. Also, the transverse displacement of the beam can be computed. This involves solving the differential matrix equation ( 29) with the Newmark algorithm where ∆t = 0.01 s, y 0 = 0 and ẏ0 = 0. The related transverse displacement field, at t = 4 s and t = 20 s, are shown in Fig. 5 along with the analytical harmonic solution, Eq. [START_REF] Zheng | Implementation of Mur's absorbing boundaries with periodic structures to speed up the design process using finite-difference time-domain method[END_REF]. At time t = 4 s, the proposed and analytical solutions appear to be slightly different which is explained since transient effects are taken into account in the proposed solution, while they are not described with the analytical solution. Those transient effects become negligible after a certain time, e.g., t = 20 s as shown in Fig. 5. In this case, the proposed solution closely matches the analytical one, as expected.

Finally, Fig. 6 shows the history of the displacement solution at position x = L/2. Again, it is seen that, after a certain time (transient period), the solution issued from the proposed approach stabilizes towards the harmonic solution. Fig. 7 presents identical results in term of velocity where the proposed solution shows better accuracy than the displacement solution. Similar trends can be found for the acceleration (not displayed here). Note that the proposed displacement solution is characterized by low frequency "ripples" (Fig. 6, left side) which are, however, less pronounced for the velocity solution. These correspond to transient phenomena of low frequencies generated by the initial conditions -i.e., the nullity of the displacement and the velocity at t = 0 s -which make the proposed solution far from the harmonic solution. The authors were able to verify through additional simulations that these ripples can be decreased by adding viscosity, e.g., by increasing ξ or by adding viscosity on the term k f . More generally, these ripples are decreasing over time as this could be observed on an interval of [0 , 40] s (not reproduced here). Periodic structures are usually made up of substructures of complex shapes, which as such cannot be modeled analytically. Such periodic structures may concern, for instance, a 2D beam with periodic distributions of holes and elastic supports (springs of stiffness K s ) as shown in Fig. 8. Here, square substructures of dimensions 2 × 2 m 2 with holes of radius 0.4 m are considered which are similar to those depicted in Fig. 2. The FE mesh of a substructure is generated on MATLAB with the DistMesh algorithm [START_REF] Persson | A Simple Mesh Generator in MATLAB[END_REF] The time response of the infinite periodic structure subject to a vertical point force F 0 (t) at x = 0 (top node) is analyzed. Two kinds of excitations are considered, i.e., a harmonic force and a Ricker wavelet of frequency 100 Hz. Within the framework of the proposed approach, a periodic structure involving N = 20 substructures and absorbing BCs is considered as shown in Fig. 8. In this case, the structure has a length of L = 40 m and is supposed to be long enough to gather useful information about the local/evanescent displacement field in the vicinity of the excitation. The related FE model has been detailed in Secs. 2.3 and 2.4. In this case, the rational approximations of the impedance matrices Z L and Z R (Eq. ( 11)) are expressed by means of P = 15 poles/residues. The time response of the structure is computed over a time range of [0 , 0.1] s by solving the differential matrix equation [START_REF] Collino | High-order absorbing boundary conditions for wave propagation models. straight line boundary and corner cases[END_REF] with the Newmark algorithm where ∆t = 10 -4 s, y 0 = 0 and ẏ0 = 0. For comparison purposes, an equivalent FE model of an "infinite" structure with a larger number of substructures ( 200) is considered and simulated over the time range [0 , 0.1] s which is supposed to be small enough to prevent wave reflections (free boundaries).

Harmonic force

Consider a harmonic point force of magnitude F 0 (t) = 10 4 cos(2π × 100t) acting at x = 0, and assume that the structure is at rest at t = 0. The time response is analyzed over a time range of [0 , 0.1] s which is supposed to be broad enough to include several oscillations (10 in this case) and cover the transient phase. Especially, the time variation of the transverse displacement at x = 20 m (right boundary, top node) can be computed as shown in Fig. 9. It is shown that the proposed solution perfectly matches the reference one over the whole time range. It is also numerically stable, i.e., a smooth curve that well predicts the oscillating nature of the signal. Additional simulations are made which concern the deformed shape of the structure at t = 0.1 s, where a scaling factor of 5 × 10 3 is used to magnify the displacement levels for a better understanding (see Fig. 10). Again, results show that the proposed solution is in perfect agreement with the reference one, which gives credit to the present work.

To further assess the accuracy of the proposed approach, the following relative error based on the displacement vector of the periodic structure (central part with 20 substructures) is investigated:

err = max n q(t n ) -q FE (t n ) 2 max n q FE (t n ) 2 , (51)
where q = [q T I q T L q T R ] T and q FE = [(q FE I ) T (q FE L ) T (q FE R ) T ] T are the displacement vectors of the periodic structure computed with the proposed approach and the FE method, respectively. Also, • 2 denotes the 2-norm. More precisely, the proposed relative error err is defined as the ratio between the maximum error q(t n )-q FE (t n ) 2 (among all the discrete times t n considered) and the maximum value of q FE (t n ) 2 (for the same discrete times). In the present case, this relative error yields err = 1.07% which, therefore, validates the proposed approach.

Ricker wavelet force pulse

Consider now a Ricker wavelet, which might represent some shock applied to the structure at x = 0: where f 0 = 100 Hz and t 0 = 5 × 10 -3 s. Again, the time history of the structure at x = 20 m can be computed as shown in Fig. [START_REF] Mencik | Multi-mode propagation and diffusion in structures through finite elements[END_REF]. In this case again, the proposed solution appears to be numerically stable and in perfect agreement with the reference FE solution. A good match between those solutions is also observed regarding the deformed shapes of the structure, e.g., at t = 0.01 s as shown in Fig. 12. Also, the relative error made for computing the displacement vector of the periodic structure with the proposed approach can be assessed via Eq. ( 51), which gives err = 1.05%. This completely validates the proposed approach.

F 0 (t) = 10 4 × 1 -2π 2 f 2 0 (t -t 0 ) 2 e -π 2 f 2 0 (t-t 0 ) 2 for t ≥ 0, (52)

Square substructures with circular rubber layers

To further highlight the relevance of the proposed approach, a second test case is considered which concerns a periodic structure made up of square substructures with circular rubber layers, as shown in Fig. 13. The material properties of the inner parts (cores) and outer parts of the substructuresi.e., inside and outside the rubber layers -are similar to those depicted in Sec. 4.1. As for the rubber layers, the properties are: inner and outer radii of 0.3 m and 0.7 m, Young's modulus of 1.5 × 10 8 Pa, Poisson's ratio of 0.48 and density of 950 kg/m 3 . The FE mesh of a substructure is shown in Fig. 13 and involves 2D linear plane stress triangular elements with n = 42 DOFs on the left and right boundaries. Each substructure therefore includes a soft layer surrounding a stiff core, and in this sense, it is likely to behave in a more complicated way than previously. Again, the time response of an infinite periodic structure consisting of an infinite number of substructures lying on periodic elastic supports (stiffness K s = 10 5 N/m) is analyzed, in a quite analogous way than previously (see Fig. 8). Here, the structure is subject to a vertical point force F 0 (t) (x = 0, top node) which represents a Ricker wavelet force pulse whose expression is similar to Eq. ( 52). Concerning the proposed approach, a periodic structure involving N = 20 substructures and absorbing BCs is considered (as before, see Fig. 8). In this case, the rational approximations of the impedance matrices Z L and Z R (Eq. ( 11)) are expressed by means of P = 19 poles/residues. Also, for comparison purposes, an equivalent FE model with 200 substructures is considered.

The time response of the structure (time range of [0 , 0.1] s) issued from the proposed approach is displayed in Fig. 14 along with the FE solution. In this case again, the proposed approach succeeds in predicting the dynamic behavior of the periodic structure over the whole time period. Fig. 14 shows the time variation of the transverse displacement at x = 20 m (right boundary, top node). The curve exhibits a main oscillation (say, between 0.02 s and 0.03 s) followed by several other oscillations of lower magnitudes which reflect the internal dynamic behavior of the substructures (motion of the cores embedded in the rubber layers). A closer look reveals that the deformed shape of the periodic structure, obtained via the proposed approach (e.g., at t = 0.01 s) also perfectly matches the reference solution, as shown in Fig. 15 where, for a better focus, only the deformations of the 10 central substructures are shown. Especially, the "core motion" phenomenon inside the substructures is well highlighted. In this case, the relative error made for computing the displacement vector of the periodic structure with the proposed approach (Eq. ( 51)) is err = 0.83%.

Conclusion

A FE procedure has been proposed to model infinite periodic structures subject to localized timedependent excitations. This involves considering a structure with a finite number of substructures encompassing the excitations, and absorbing BCs to model semi-infinite structures outside the excitation region. The WFE method has been used to express the absorbing BCs by means of impedance matrices, in the frequency domain. Those impedance matrices have been rewritten in terms of polynomials of the frequency iω up to order 2 which can be simply converted to the time domain. The procedure works by decomposing the impedance matrices in terms of rational functions, and by adding vectors of supplementary variables. As a result, a periodic structure can be modeled using a classical second-order time differential equation which can be integrated with the Newmark algorithm. Numerical experiments have been proposed which have clearly demonstrated the accuracy of the approach.

Several test cases have been considered which concern an infinite Euler-Bernoulli beam on an elastic foundation subject to a point harmonic force, and two periodic structures with 2D substructures subject to different kinds of point forces (harmonic force and Ricker wavelet). Follow-on works could include the analysis of infinite periodic structures with localized nonlinear effects, e.g., subject to fast loadings of high magnitudes (shocks, blast).

Figure 1 :

 1 Figure 1: Schematic of a periodic structure with N = 10 substructures subject to time-dependent excitations and semiinfinite BCs.

Figure 2 :

 2 Figure 2: Periodic structure of infinite length, and FE mesh of a substructure.

Figure 3 :

 3 Figure 3: Infinite flexural beam on an elastic foundation.

Figure 4 :

 4 Figure 4: Real and imaginary parts of the components of the impedance matrix ZR: (green crosses) Analytical solutions; (blue line) proposed solutions with P = 3; (red line) proposed solutions with P = 7; (orange line) proposed solutions with P = 10; (purple line) proposed solutions with P = 12.

Figure 5 :

 5 Figure 5: Transverse displacement field of the beam on the elastic foundation at t = 4 s (left) and t = 20 s (right): (blue crosses) Proposed approach; (red line) analytical harmonic solution.

Figure 6 :

 6 Figure 6: Time response at x = L/2 for t ∈ [0 , 20] s (left) and t ∈ [19 , 20] s (right). (blue line and crosses) Proposed approach; (red line) analytical theory, harmonic response.

Figure 7 :

 7 Figure 7: Time response of the velocity at x = L/2 for t ∈ [0 , 20] s (left) and t ∈ [19 , 20] s (right). (blue line and crosses) Proposed approach; (red line) analytical theory, harmonic response.

4 .

 4 Periodic structure with 2D substructures 4.1. Square substructures with holes

  , and involves 2D linear plane stress triangular elements (three nodes) with n = 26 DOFs on the left and right boundaries. Regarding the modeling of the periodic supports, a nodal stiffness of K s /2 (vertical direction) is added to the FE model of the substructures at the left and right boundaries (bottom node). Other substructure parameters are: thickness e = 0.005 m, Young's modulus E = 7 × 10 10 Pa, Poisson's ratio ν = 0.35, density ρ = 2700 kg/m 3 , and stiffness K s = 10 5 N/m. Rayleigh-type damping matrices C s = aM s + bK s (see Sec. 2.2) are also considered where a = 0.01 s -1 and b = 5 × 10 -5 s.

Figure 8 :

 8 Figure 8: Schematic of an infinite periodic structure (square substructures with holes) with periodic elastic supports.

Figure 9 :

 9 Figure 9: Periodic structure (square substructures with holes) subject to a harmonic force: time response at x = L/2 (top node, vertical displacement). (blue crosses) Proposed approach; (red line) reference FE method.

Figure 10 :

 10 Figure 10: Periodic structure (square substructures with holes) subject to a harmonic force: deformed shape at t = 0.1 s.

Figure 11 :

 11 Figure 11: Periodic structure (square substructures with holes) subject to a Ricker wavelet pulse: time response at x = L/2 (top node, vertical displacement). (blue crosses) Proposed approach; (red line) reference FE method.

Figure 12 :

 12 Figure 12: Periodic structure (square substructures with holes) subject to a Ricker wavelet pulse: deformed shape at t = 0.01 s.

Figure 13 :

 13 Figure 13: FE mesh of a substructure with a circular rubber layer (in blue).

Figure 14 :

 14 Figure 14: Periodic structure (square substructures with circular rubber layers) subject to a Ricker wavelet pulse: time response at x = L/2 (top node, vertical displacement). (blue crosses) Proposed approach; (red line) reference FE method.

Figure 15 :

 15 Figure 15: Periodic structure (square substructures with circular rubber layers) subject to a Ricker wavelet pulse: deformed shape at t = 0.01 s.

Table 1 :

 1 Relative errors between the proposed and analytical solutions (components of the impedance matrix ZR) for different orders P . Component (1,2) 0.106 4.51 × 10 -3 2.54 × 10 -4 6.83 × 10 -6Component (2,2) 0.082 6.90 × 10 -3 5.49 × 10 -4 3.47 × 10 -6

	Relative error


Ngoc Bich
email: bich-ngoc.ho@univ-pau.fr

David Ho

Camélia Matei Pino-Perez

Joseph Ghimbeu

Deneb Diaz

Cécile Peredo-Mancilla

David Hort

Bessieres

Ho

David Pino Perez
email: david.pino-perez@univ-pau.fr

Deneb Peredo Mancilla

Joseph Diaz

Cécile Hort

David Bessieres

Camélia Matei Ghimbeu

Determination of methane, ethane and propane on activated carbons by experimental pressure swing adsorption method

Keywords:

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Climate change is one of the most considerable environmental threats nowadays. Therefore, the reduction of air pollutant emission has attracted a lot of attention as a global challenge. Taking into account that petroleum and coal are naturally unclean, natural gas is thought to be a promising alternative thanks to its availability, adaptability, high energy value, and low carbon dioxide emission upon combustion. In 2019, the world annual natural gas consumption reached 3,929,000 million m 3 , making it become the third largest energy source [START_REF] Abdullah | Methane enrichment from high carbon dioxide content natural gas by pressure swing adsorption[END_REF][START_REF] Li | Selective gas adsorption and separation in metal-organic frameworks[END_REF]. Raw natural gas contains prominently up to 70-90% of methane (CH 4 ), other alkanes such as ethane (C 2 H 6 , 6.4%) and propane (C 3 H 8 , 5.3%), butane (C 4 H 10 , 1.4%) and other non-hydrocarbons such as carbon dioxide (CO 2 , 5%) [START_REF] Faramawy | Natural gas origin, composition, and processing: A review[END_REF]. However, the presence of compounds other than methane decreases the combustion heat capacity of the gas and can lead to the damage of the pipelines, they must be eliminated from natural gas before application [START_REF] Wang | Model and experimental data research of natural gas storage for vehicular usage[END_REF][START_REF] Koonaphapdeelert | Low pressure biomethane gas adsorption by activated carbon[END_REF].

Up to date, numerous physical and chemical treatments for natural gas, some regular methods such as membrane separation, water scrubbing, physical and chemical absorption have been applied in both laboratory and industrial scales and they all have proved their pros and cons [START_REF] Kadam | Recent advancement in biogas enrichment and its applications[END_REF][START_REF] Vrbová | Upgrading Biogas to Biomethane Using Membrane Separation[END_REF][START_REF] Sahota | Review of trends in biogas upgradation technologies and future perspectives[END_REF]. Water scrubbing is commonly used thanks to its simple operation and high CH 4 purity outcome, but it requires a huge amount of water use. On the other hand, chemical adsorption by using amines can achieve CH 4 of 98% purity but its needs of chemicals and regeneration energy are undesirable. In order to overcome these issues, pressure swing adsorption (PSA) is suggested due to its effectiveness and its capacity of undesired compounds elimination [START_REF] Pahinkar | A novel temperature swing adsorption process for natural gas purification: Part I, model development[END_REF]. PSA can be used to enrich the methane concentration to enhance the calorific value and to meet the gas quality requirement. Moreover, PSA is also promising in sustaining gas production and minimizing pollution problems due to gas flaring [START_REF] Angelidaki | Biogas upgrading and utilization: Current status and perspectives[END_REF][START_REF] Murad | Flared Gas Emission Control from an Oil Production Platform[END_REF].

Consequently, the exploration of optimum adsorbents for natural gas purification is in high demand. In our days, material screening is carried out not only in adsorption experiments but also in artificial intelligence such as molecular simulations and machine learning in order to accurately predict the structural relationship between adsorption capacity and structural properties of adsorbents [START_REF] Pardakhti | Machine Learning Using Combined Structural and Chemical Descriptors for Prediction of Methane Adsorption Performance of Metal Organic Frameworks (MOFs)[END_REF][START_REF] Simon | The materials genome in action: identifying the performance limits for methane storage[END_REF][START_REF] Lei | Selective adsorption of CO2 by Hex-star phosphorene from natural gas: Combining molecular simulation and real adsorbed solution theory[END_REF]. Most of absorbents studied are zeolites [START_REF] Weh | Nitrogen rejection from natural gas by dual reflux-pressure swing adsorption using activated carbon and ionic liquidic zeolite[END_REF][START_REF] Arefi Pour | Adsorption separation of CO2/CH4 on the synthesized NaA zeolite shaped with montmorillonite clay in natural gas purification process[END_REF][START_REF] Zhu | Modification of zeolite by metal and adsorption desulfurization of organic sulfide in natural gas[END_REF], silica gels [START_REF] Sui | Adsorption and desorption of binary mixture of acetone and ethyl acetate on silica gel[END_REF], metal-organic frameworks (MOFs) [START_REF] Li | Selective gas adsorption and separation in metal-organic frameworks[END_REF][START_REF] Bao | Kinetic separation of carbon dioxide and methane on a copper metal-organic framework[END_REF] and activated carbons [START_REF] Esteves | Adsorption of natural gas and biogas components on activated carbon[END_REF][START_REF] Oginni | Effect of one-step and two-step H3PO4 activation on activated carbon characteristics[END_REF][START_REF] Tang | Development of activated carbon pellets using a facile low-cost binder for effective malachite green dye removal[END_REF][START_REF] Queiroz | Activated carbon obtained from amazonian biomass tailings (acai seed): Modification, characterization, and use for removal of metal ions from water[END_REF][START_REF] Matsui | Removal of siloxane from digestion gas of sewage sludge[END_REF]. Amongst them, activated carbons have shown their advantageous result comparing to other absorbents because their production cost are less expensive than others while they still have greater adsorption capacity [START_REF] Petersson | 14 -Biogas cleaning[END_REF]. Activated carbons can be produced by pyrolysis and activation of natural sources such as rice husk [START_REF] Menya | Production and performance of activated carbon from rice husks for removal of natural organic matter from water: A review[END_REF], tree wood, bamboo, charcoal, or sometimes rare carbonaceous plant-derived as butnea pollen [START_REF] Ahmed | Investigation on activated carbon derived from biomass Butnea monosperma and its application as a high performance supercapacitor electrode[END_REF] or clays [START_REF] Pinto | Porous Materials Prepared from Clays for the Upgrade of Landfill Gas[END_REF]. They also have been used in water and air purification, dentistry, catalytic chemistry for years, which therefore is supposed to give a promising result in gas pressure swing adsorption technology [START_REF] Rezvani | Activated carbon surface modification by catalytic chemical vapor deposition of natural gas for enhancing adsorption of greenhouse gases[END_REF].

A lot of experiments were performed on mostly dominant compounds in natural gas, methane, and carbon dioxide combustion in the research on absorbents for purification of the raw natural gas [START_REF] Chen | Effect of swing on removing CO2 from offshore natural gas by adsorption[END_REF][START_REF] Jiang | Post-combustion CO2 capture from a natural gas combined cycle power plant using activated carbon adsorption[END_REF][START_REF] Carchini | A theoretical study of gas adsorption on calcite for CO2 enhanced natural gas recovery[END_REF]. However, The elimination of ethane and propane, two other compounds in natural gas, is also desirable to increase its heat combustion [START_REF] Speight | Liquid fuels from natural gas[END_REF]. Furthermore, the study of gas adsorption on activated carbon for other hydrocarbons thus allows us to comprehend their potential performance for natural gas purification, hence lead to effective solutions on an industrial scale.

Previously published works have only generally shown that the adsorption of gas molecules onto the surface of activated carbon is influenced by numerous factors such as surface chemistry, morphology, and affinity of the adsorbate with the adsorbent [START_REF] Abedini | Adsorption of propane and propylene on M-MOF-74 (M=Cu, Co): Equilibrium and kinetic study[END_REF][START_REF] Peredo-Mancilla | CO2 and CH4 Adsorption Behavior of Biomass-Based Activated Carbons[END_REF]. In this context, the objective of this study is to investigate methane, ethane, and propane isothermal adsorption on five commercial activated carbons and to deduce the relationship between the textural properties of adsorbents and their adsorption performance in both low and high working pressures. The main physical properties of these materials, experimental set-up, and the methodology are detailed in Section Materials and methods. In Results and Discussion, an analysis of the results is provided and the relationship between the structural properties of the activated carbons and the adsorption capacity is deduced. Besides, the experimental outcome was fitted to the modified Langmuir model and compared to various reported results in the literature.

Materials and method

Gas cylinders and activated carbons

Methane, ethane, propane, and helium were purchased from Linde AG (France), with a purity up to 99.999% (purity denoted 5.0) as described in Table 1. They were conducted to the system via a pressure regulator of the same provider. Propane was exceptionally pumped into a smaller gas tank under a counter-current heating for pre-gasification. Additional purification was not required. Five Norit commercial activated carbons named RX 1.5, ROx 0.8, CNR-115, CGRAN, and GAC 1240 were purchased from Cabot Corporation (USA). CGRAN and CNR-115 were chemically activated by H 3 PO 4 while RX 1.5, ROx 0.8, and GAC 1240 were activated by steam. The morphology of each activated carbon surface was analyzed by scanning electron microscopy (SEM) using an FEI model Quanta 400 analyzer.

Since the activated carbons were provided commercially, they are considered stable i.e do not lose adsorption capacity over numerous adsorption-desorption cycles. On the industrial scale, depending on the frequency of use, the service time of Norit activated carbons can last for months. Moreover, exhausted commercial activated carbons are able to be thermally or chemically reactivated and returned to 97% or greater compared to virgin carbon with reduced cost [START_REF] Corporation | Reactivation Services[END_REF][START_REF] Koehlert | Activated Carbon: Fundamentals and New Applications[END_REF]. On the laboratory scale, the durability of Norit activated carbons was reported to be constant in Plaza et al. and Aumeier et al., where 5-10 cycles of adsorption-desorption cycles were carried out for different configurations of adsorption [START_REF] Plaza | Post-combustion CO2 capture with a commercial activated carbon: Comparison of different regeneration strategies[END_REF][START_REF] Aumeier | Preliminary Study on the Application of Temperature Swing Adsorption in Aqueous Phase for Pesticide Removal[END_REF]. The adsorption capacities of methane, ethane, and propane onto the five activated carbon in this study are hence judged to remain steady.

Experimental system and measurement protocol 2.2.1. System set-up

Figure 1 represents A schematic plan of the homemade adsorption system. This system is composed of a dosing cell and an adsorption cell manufactured by Top Industries. Both of them are made from stainless steel and have similar volume (18.5 ± 0.2

3 ). These two cells were put in series and coupled with a PR 400B-S manometer of MKS Instruments, Inc. This manometer is operating in a wide range of pressure 0-3.3 MPa with an accuracy of 0.25% of reading. A heating wire and isolation foam wrap around the whole system to ensure the homogeneity of temperature. The system is also integrated with a temperature sensor and display devices for isotherm monitoring.

Sample preparation

An amount of activated carbon (0.9 ± 0.02 g) was introduced directly into the cleaned adsorption cell and put under vacuum using a vacuum pump at 90°C for at least 8 hours. This treatment at high temperature is required for activated carbons to ensure that all water and residual adsorbate molecules are removed from the pores thus guarantee the reliability of the adsorption isotherms.

Once the adsorption was finished, the desorption step was conducted under vacuum at 90°C for 8 hours again in order to regenerate the performance of the adsorbent.

Measurement principle 2.3.1. Calibration

Based on the fact that activated carbons do not adsorb noble gas, a helium expansion was undertaken from the dosing cell into the adsorption cell before each measurement to determine the adsorption volume of the system possessing an activated carbon. Knowing the dosage volume in advance by carbon dioxide calibration, by this consecutive measurement, the calculated adsorption volume is usually inferior to the one without sample.

=

(1)

where is the measured pressure of the dosing cell when V2 is closed and is the total pressure in the whole equipment when V2 is opened. The total volume of the system, denoted is hence calculated:

= + (2)

Adsorption isotherms

A dose of gas is released into the dosing cell (V1 opened, V2 closed). Once the equilibrium in the dosing cell is stabilized ( ), the gas is then released to the adsorption cell by opening V2.

For the first point of the isotherm ( = 1), the number of moles of the adsorbate gas introduced in the dosing cell 1 is calculated as in Equation 3:

1 = ⋅ 1 (3)
where 1 , which is also denoted 1 ( , ) is the molar density of the adsorbate at a particular initial pressure and constant working temperature. ( , ) is calculated by using Equation 4, which is considered as a polynomial function of at 303K and its parameters were obtained from NIST database [40]. The coefficients of different pure gases are showed in Table 2.

( , ) = ⋅ 4 + ⋅ 3 + ⋅ 2 + ⋅ + (4)
The adsorption isotherms are described through an accumulative process. The adsorption equilibrium is reached when the pressure is stabilized again ( ), the molar quantity of gas adsorbed onto the surface of the adsorbent 1 is calculated as in Equation 5:

1 = 1 - ⋅ 1 (5)
The procedure is repeated by closing V2 and introducing a new dose of gas into the dosing cell. and at the -th point are calculated as follows [START_REF] Gasparik | High-Pressure Methane Sorption Isotherms of Black Shales from The Netherlands[END_REF][START_REF] Peredo-Mancilla | Experimental Determination of the CH4 and CO2 Pure Gas Adsorption Isotherms on Different Activated Carbons[END_REF]:

= ⋅ -1 + ⋅ (6)

Adsorption isotherms fitting parameters

The Langmuir model is the most popular model used to represent experimental isotherms by assuming that the adsorbate behaves as an ideal gas in isothermal conditions and the adsorbent has an ideal chemically homogeneous solid surface [START_REF] Swenson | Langmuir's Theory of Adsorption: A Centennial Review[END_REF][START_REF] Awad | Adsorption of organic pollutants by natural and modified clays: A comprehensive review[END_REF]. In this context, the modified Langmuir model was chosen, where we assumed that during the adsorption procedure, the attachment gas molecules onto the surface of activated carbons in the form of a single-layer coverage. We denoted the excess adsorbed quantity of gas at the working pressure called P. The experimental results were then compared to the calculated estimation by using the modified Langmuir model with two parameters [START_REF] Li | Comprehensive polynomial simulation and prediction for Langmuir volume and Langmuir pressure of shale gas adsorption using multiple factors[END_REF]:

= × + 1 - ( , ) (7)
where , called the Langmuir pressure, is the pressure where the adsorbate occupied half of the sites the surface of the activated carbon and n L is the maximum amount absorbed, corresponding to the occupancy of all sites. In this study, the value of had been fixed at 423 -3 for methane, 546 -3 for ethane and 581 -3 for propane [START_REF] Ortiz Cancino | Adsorption of pure CO2 and a CO2/CH4 mixture on a black shale sample: Manometry and microcalorimetry measurements[END_REF][START_REF]PubChem, Ethane -structure, chemical names, physical and chemical properties, classification, patents, literature, biological activities, safety/hazards/toxicity information, supplier lists, and more[END_REF][START_REF]Occupational Safety and Health of the German Social Accident Assurance[END_REF].

Since the Langmuir fitting parameters are determined considering the confirmation of fitted Langmuir equation in low working pressure and its acceptance in high-pressure zone, a deviation standard was thus required and calculated as Equation 8:

Δ = 1 √ √ √ √ ∑ =1 ( ) - 2 (8)
where Δ stands for the standard deviation of the sample, N is number of values observed and is the corresponding mean value calculated by the modified Langmuir equation 7 [START_REF] Ettarh | Common descriptive and analytical statistics in investigative studies[END_REF]. Along with the SEM images, the N 2 adsorption/desorption isotherms at 77 K for 5 activated carbons as well as their pore size distributions were published previously in Peredo et al [START_REF] Peredo-Mancilla | Experimental Determination of the CH4 and CO2 Pure Gas Adsorption Isotherms on Different Activated Carbons[END_REF]. The N 2 adsorption/desorption capacities of the five ACs range widely from 300 to 600 3 -1 AC at ∕ 0 = 1. CNR-115 is found to have the highest N 2 adsorption capacity while GAC 1240 absorbed N 2 significantly less than the other activated carbons. Indeed, the isotherms on RX 1.5, ROX 0.8, CNR-115, and GAC 1240 give type I behavior, where the presence of a saturation limit indicates the completion of a single monolayer of N 2 at the material surface [START_REF] Lowell | Adsorption isotherms[END_REF]. Meanwhile, CGRAN carries type IV isotherm, representing a finite multi-layer formation corresponding to a complete filling of the capillaries, which leads to the deduction of a wide distribution of pore sizes for this AC [START_REF] Sultan | Optimization of adsorption isotherm types for desiccant air-conditioning applications[END_REF].

Results and Discussion

Material properties

In terms of pore size distribution, the particle size of activated carbons also differs widely [START_REF] Peredo-Mancilla | Experimental Determination of the CH4 and CO2 Pure Gas Adsorption Isotherms on Different Activated Carbons[END_REF]. While CGRAN and CNR-115 own a broad distribution of pore sizes from 0.7 nm to 1.2 nm in one peak, RX 1.5, ROX 0.8 have 2 distinct peaks at 0.7 nm and 1.8 nm. Meanwhile, GAC 1240 possesses a single peak of micropores at 0.6 mn, indicating that this activated carbon is more microporous than the others. These findings are consistent with the previous suggestion about CGRAN, that its pores rather have similar sizes.

On the other hand, the influence of porosity and specific surface of absorbent has been proven to affect gas adsorption capacity [START_REF] Belhachemi | Comparison of NO2 removal using date pits activated carbon and modified commercialized activated carbon via different preparation methods: Effect of porosity and surface chemistry[END_REF][START_REF] Cao | Influences of micro-pores and meso-pores on elastic and plastic properties of porous materials[END_REF]. In this case, the values of BET areas vary from 982 (GAC 1240) to 1714 2 -1 (CNR-115), belong to usual value range of activated carbons [START_REF] Oginni | Effect of one-step and two-step H3PO4 activation on activated carbon characteristics[END_REF][START_REF] Saka | SEM, iodine number analysis and preparation of activated carbon from acorn shell by chemical activation with ZnCl2[END_REF][START_REF] Hidayu | Characterization of Activated Carbon Prepared from Oil Palm Empty Fruit Bunch Using BET and FT-IR Techniques[END_REF][START_REF] Kharrazi | A novel post-modification of powdered activated carbon prepared from lignocellulosic waste through thermal tension treatment to enhance the porosity and heavy metals adsorption[END_REF][START_REF] Wang | Preparation of activated carbon from peanut shell with KOH activation and its application for H2S adsorption in confined space[END_REF]. Moreover, along with the highest BET areas, CNR-115 and RX 1.5 also own a great proportion of microporous volume (0.64 and 0.61 3 -1 respectively). CGRAN has the most important total pore volume (0.99 3 -1 ) but also the highest mesopore volume (0.54 3 -1 ) due to the high amount of H 3 PO 4 during activation, which leads into the average value of BET area [START_REF] González-García | Activated carbon from lignocellulosics precursors: A review of the synthesis methods, characterization techniques and applications[END_REF]. This result approved the two previously published experimental data of nitrogen adsorption/desorption isotherms and pore size distributions, where it was suggested that there was the presence of an important volume of mesopore in CGRAN [START_REF] Peredo-Mancilla | Experimental Determination of the CH4 and CO2 Pure Gas Adsorption Isotherms on Different Activated Carbons[END_REF]. Meanwhile, ROX 0.8 possesses a low mesopore volume (0.16 3 -1 ) but also a relatively low micropore volume ( 483 -1 ), therefore carries the similar BET areas value as CGRAN. In the case of GAC 1240, due to the low proportion of and the low total pore volume (0.56 3 -1 , this activated carbon is considered least porous having the smallest BET area (982 2 -1 ).

Adsorption isotherms

Tables 3 -5 and Figure 3 present the experimental results of the 5 above-described activated carbons on methane, ethane, and propane at 303K. All isotherms belong to type I which indicates the monolayer isotherms of physisorption as described in Langmuir's theory of adsorption. In other words, the adsorption uptake increased rapidly at low pressure and then saturated upon the pressure augmentation when the pores are gradually filled [START_REF] Ng | A Universal Isotherm Model to Capture Adsorption Uptake and Energy Distribution of Porous Heterogeneous Surface[END_REF]. All activated carbons showed a major difference between the adsorption capacities, as the adsorbed quantity of ethane and propane was significantly higher than methane. Indeed, in the case of RX 1.5, the adsorption capacity for ethane achieved 9.12 -1 at 2.41 MPa while there was only 7.25 mol of methane adsorbed per kilogram of this activated carbon. In addition, at a low working pressure (0.50 MPa), the propane adsorption capacity already leveled up to 8.09 -1 while the ones for methane and ethane are 3.68 and 7.37

(MPa) ( -1 ) (MPa) ( -1 ) (MPa) ( -1 ) (MPa) ( -1 ) (MPa) ( -
-1 respectively. Likewise, at 1.9 MPa, CNR-115 underwent an adsorption capacity of 8.79

-1 for ethane in relation to 4.17 -1 for methane whilst the same capacity for propane attained at 0.57 MPa. ROx 0.8 had a less important adsorption capacity when 1 kilogram adsorbed only 5.14 mole of methane or 6.44 mol of ethane at 1.2 MPa. The same results also happened to CGRAN, having 5.08 and 7.71

-1 as methane and ethane adsorption capacities and reached 7.26

-1 as propane capacities at 0.604 MPa. In the same way, at 1.0 MPa the measured CH 4 adsorption capacity of GAC 1240 was 5.23

-1 while for ethane and propane, it adsorbed 6.7 and 7.6

-1 respectively. This phenomenon can be attributed that at low and intermediate pressures, the interaction between propane and the surface of the adsorbent is much higher than methane. Being linear molecules, the more carbon atoms in the chain, the stronger the covalent bond created between alkane molecules and the surface of the activated carbon. For that reason, at low pressure, propane gets adsorbed the most, thus ethane adsorbs the second and methane adsorbed the least [START_REF] Ponraj | Separation of methane from ethane and propane by selective adsorption and diffusion in MOF Cu-BTC: A molecular simulation study[END_REF]. This explanation is also consistent with the argument of Pino-Perez et al., that porous materials adsorb preferentially long-chain molecules than short-chain molecules in a study on propane, pentane, and heptane adsorption by using Ecosorb activated carbon [START_REF] Pérez | Étude Expérimentale de l'Adsorption de gaz sur des solides micro ou meso-poreux à haute pression[END_REF]. In addition, the strong interactions of ethane and propane with the adsorbents can also be described by Van der Waals's attractive force based on the polarizability of the gas molecules. Indeed, since the polarizability of propane (62.9 × 10 25 3 ) is greater than ethane (44.3 × 10 25 3 ) and methane (25.93 × 10 25 3 ), propane is expected to be the most attracted onto the surface of adsorbents [START_REF] Li | Selective gas adsorption and separation in metal-organic frameworks[END_REF]. Despite the higher uptakes, a plateau was observed for C 2 H 6 and C 3 H 8 isotherms of most activated carbons i.e the maximal filling of micro-pores at 303K. Once entering a higher pressure zone, bigger pores start to be filled by the adsorbate. As propane and ethane have longer carbon chains, their rotation and translation get limited, leading to rapid saturation compared to methane. The saturation is clearly observed in the case of RX 1.5, ROX 0.8, and GAC 1240, whose volume of mesopores are low. They possess 0.20, 0.16 and 0.20 3 -1 as mesopores volume respectively. This finding is in agreement with the study of Mhaskar et al., where the C 3 H 8 adsorption capacity on silica gel in the range of 0 -1 bar is significantly higher than ethane and methane ones [START_REF] Mhaskar | Natural Gas Treatment Using Adsorptive Separation, Natural Gas -Extraction to End UsePublisher[END_REF].

In general, RX 1.5, ROx 0.8 and GAC 1240 adsorbed methane molecules more than CNR-115 and CGRAN. Indeed, the surface acidity of CNR-115 and CGRAN are reported to be higher than other commercial activated carbons due to the initial chemical activation [START_REF] Peredo-Mancilla | Comparative study of the CH4/CO2 adsorption selectivity of activated carbons for biogas upgrading[END_REF]. As methane is an inert compound, together with the pore blockage of the functional groups, the acidic surface tends to adsorb less than basic ones. This outcome is consistent with the findings of Contreras et al. that the presence of acidic functionalities disfavors the CH 4 uptakes and the basic surface is thought to be responsible for the improvement on CH 4 adsorption [START_REF] Contreras | On the methane adsorption capacity of activated carbons: in search of a correlation with adsorbent properties[END_REF]. In reference to the CO 2 adsorption capacity, it also seems to be coherent to Peredo et al. [START_REF] Peredo-Mancilla | Experimental Determination of the CH4 and CO2 Pure Gas Adsorption Isotherms on Different Activated Carbons[END_REF]. As the surface of RX 1.5 and ROx 0.8 are more basic, they adsorbed more CO 2 pure gas than CNR-115 and CGRAN. On the other hand, it was observed that not only the C 2 H 6 and C 3 H 8 isotherms, but also the CH 4 isotherm of GAC 1240 saturated considerably faster than other ACs due to its minimal micro-pores and mesopores volumes. Indeed, since GAC 1240 owns only 0.36 and 0.20 3 -1 micro-pores and mesopores respectively, at low working pressure, most of its pores filled quickly, leading to fewer spaces for gas molecules at high pressure. The same scenario happened to C 2 H 6 isotherm on RX 1.5 and ROx 0.8, the ethane uptake saturated at 2 MPa, justifying the inferior mesopores volume (0.20 3 -1 and 0.16 3 -1 respectively).

(MPa) ( -1 ) (MPa) ( -1 ) (MPa) ( -1 ) (MPa) ( -1 ) (MPa) ( -
In the case of C 3 H 8 isotherms, it is necessary to note that the interval of working pressure is reduced due to the limit of gas-liquid transition of propane. However, it was noticed that the results for C 3 H 8 are in coherence with CH 4 and C 2 H 6 , especially in the case of low mesoporous activated carbons like RX 1.5 and ROx 0.8. Here, RX 1.5 showed the highest propane uptakes and saturated at 0.7 MPa. The saturation was also witnessed clearly in ROx 0.8 at 0.5 MPa.

Langmuir fitting parameters

Langmuir fitting parameters of three adsorbates for five adsorbents are presented in Table 6. These parameters are considered accurate as the deviations Δ are relatively low. Furthermore, the relationship between the BET surface of each activated carbon and its maximum Langmuir adsorption capacity is also presented in Figure 4a. This linear regression proved an adequate correlation for methane (blue squares, 2 = 0.902), ethane (red diamonds, 2 = 0.968) and propane (yellow diamonds, 2 = 0.980). In all cases, the larger the BET surface is, the more important the maximum Langmuir becomes. The relation between n L and the pore volumes of each activated carbons including micropores, mesopores, and total pores are also established in Figure 4b -4d. The dependence of n L on V micro also show and analogous trend i.e the adsorption capacities increase along with the increment of V micro . However, there is no clear linear relation is noticed like the BET surface area. This outcome is consistent with the precedent findings of Peredo et al. and Ortiz, the higher BET surface area is responsible for the higher maximum adsorption capacity since it is related directly to the number of available physisorption sites [START_REF] Peredo-Mancilla | Experimental Determination of the CH4 and CO2 Pure Gas Adsorption Isotherms on Different Activated Carbons[END_REF][START_REF] Cancino | Etude expérimentale de l'adsorption du méthane dans des gaz de schistes colombiens et de la séparation méthane/dioxyde de carbone[END_REF]. [START_REF] Peredo-Mancilla | CO2 and CH4 Adsorption Behavior of Biomass-Based Activated Carbons[END_REF][START_REF] Vargas | Carbon dioxide and methane adsorption at high pressure on activated carbon materials[END_REF][START_REF] Shirazani | Starch-based activated carbon micro-spheres for adsorption of methane with superior performance in ANG technology[END_REF]. As observed, the methane uptakes of activated carbons range widely depending on the textural properties and surface chemistry of the adsorbent. For example, AC-H 3 PO 4 is an olive stones-based activated carbon activated by H 3 PO 4 that had inferior specific surface area (1178 2 -1 ) but superior methane adsorption performance. This outcome can be attributed to the original surface chemistry of the precursor and the activation process. In the same way, MP48 was also activated by H 3 PO 4 that had a lower BET surface (1368 2 -1 ) and micropore volume (0.37 3 -1 ) but the methane uptake of MP48 at 1 MPa was only 2.4 -1 compared to 3.1 -1 of CNR-115. On the other hand, the morphology of MCa2 was close to CNR-115 (S BET = 1700 2 -1 , V micro = 0.64 3 -1 ) but this activated carbon was activated by CaCl 2 , which may improve the affinity of its surface with CH 4 molecules and enhanced the adsorption capacity [START_REF] Vargas | Carbon dioxide and methane adsorption at high pressure on activated carbon materials[END_REF]. In the case of HS-KOH1:2, the KOH treatment is proved to increase the CH 4 adsorption by widening pores and adjusting the affinity with the adsorbate. However, the CH 4 isotherm of HS-KOH1:2 was quickly saturated, indicating that this adsorbent did not possess enough mesopores to adsorb methane at high pressure [START_REF] Shirazani | Starch-based activated carbon micro-spheres for adsorption of methane with superior performance in ANG technology[END_REF].

Sample Compound p L n L Δn (MPa) ( -1 )

Ethane adsorption

Liang et al. [START_REF] Liang | Asphalt-derived high surface area activated porous carbons for the effective adsorption separation of ethane and ethylene[END_REF][START_REF] Liang | Ethane selective adsorbent Ni(bdc)(ted)0.5 with high uptake and its significance in adsorption separation of ethane and ethylene[END_REF], Wang et al. [START_REF] Wang | Novel C-PDA adsorbents with high uptake and preferential adsorption of ethane over ethylene[END_REF] and Pires et al. [START_REF] Pires | Ethane Selective IRMOF-8 and Its Significance in Ethane-Ethylene Separation by Adsorption[END_REF] carried out previously ethane adsorption isotherms on different activated carbons from various precursors and activation agents. The experimental outcomes are shown on in Table 8. At 1 MPa, RX 1.5 and CNR-115 own superior ethane adsorption capacity over C-700-3 and C-800-3 despite the lower BET surfaces. It is because C-700-3 and C-800-3 are rather mesoporous material due to the excess use of KOH as activation agent, which resulted in the low ethane uptakes at low pressure but allows to adsorb ethane at high pressure without being saturated. Similarly, since C-PDA-3 was activated with an increase in KOH/C ratio, some micropores were transformed into mesopores, hence the ethane uptakes is reduced at 1 MPa but it still has room for adsorption at high working pressure. On the other hand, as observed, RX 1.5 and CNR-115 are more microporous hence adsorbed ethane more quickly but also had a faster saturation tendency.

In general, activated carbons have shown their promising adsorption capacity over MOFs. Liang et al. and Pires et al. [START_REF] Liang | Ethane selective adsorbent Ni(bdc)(ted)0.5 with high uptake and its significance in adsorption separation of ethane and ethylene[END_REF][START_REF] Pires | Ethane Selective IRMOF-8 and Its Significance in Ethane-Ethylene Separation by Adsorption[END_REF] reported the ethane adsorption capacities of two MOFs Ni(bdc)(ted) 0.5 and IRMOF-8 were 5.0 and 4.0 -1 respectively. Meanwhile, the specific surface area ad total pore volume of Ni(bdc)(ted) 0.5 were 1701 2 -1 and 0.79 3 -1 , which were relatively similar to the morphology of RX 1.5 and CNR-115 [START_REF] Liang | Ethane selective adsorbent Ni(bdc)(ted)0.5 with high uptake and its significance in adsorption separation of ethane and ethylene[END_REF]. Therefore, the inferior performance of MOFs is not only due to the lower BET area but also the weak affinity of adsorbate with the MOF surface. The adsorption capacity is hence thought to be resulted from various factors such as morphology, surface chemistry of the adsorbent and physico-chemical properties of the adsorbate.

Propane adsorption

Regarding the propane adsorption, RX 1.5 and CNR-115 also show advantageous capacities in comparison with the MOFs Cu-MOF-74 and Co-MOF-74 by Abedini et al. and the activated carbon Ecosorb by Pino et al. [START_REF] Abedini | Adsorption of propane and propylene on M-MOF-74 (M=Cu, Co): Equilibrium and kinetic study[END_REF][START_REF] Pino | Adsorption of Alkanes in the Dense Vapor Phase on a Microporous Activated Carbon[END_REF]. However, a quick saturation of C 3 H 8 adsorption isotherms was observed for all adsorbents, indicating the total micropore filling at low temperature and low pressure due to the high affinity between the long-chain molecules and the force field of the adsorbent.

As with the ethane uptakes, the propane adsorption capacity depends predominantly on the morphology of the adsorbent. Since the BET area and average pore size of the activated carbon Ecosorb were 1290 2 -1 and 0.46 nm, this adsorbent was more microporous than CNR-115 and RX 1.5 and did not own too much available sites for adsorption at high pressure. Consequently, at 0.5 MPa, the propane uptakes of Ecosorb was 5.7

-1 while RX 1.5 and CNR-115 reached 8.0-8.1 -1 . Similarily, in addition to the inferior BET area and micropore volume, Cu-MOF-74 and Co-MOF-74 also had weaker affinity with propane, leading into lower propane uptakes compared to activated carbons.

Conclusion

The adsorption isotherms of natural gas components which are methane, ethane, and propane were successfully carried out for a series of five commercial activated carbons. Most activated carbons presented the type I isotherm, indicating the highly porous adsorbent with monolayer adsorption. Regarding the three adsorbates, the increase in chain length leads to the preferential adsorption of gaseous molecules as reported in literature [START_REF] Abedini | Adsorption of propane and propylene on M-MOF-74 (M=Cu, Co): Equilibrium and kinetic study[END_REF][START_REF] Pino | Adsorption of Alkanes in the Dense Vapor Phase on a Microporous Activated Carbon[END_REF]. Therefore, all activated carbons showed a greater adsorption capacity on propane than ethane and methane at low pressure. However, once the experimental pressure enters the higher zone i.e micropores are filled, the adsorption of gas molecules on mesoporous adsorbents tends to saturate slower.

Among the studied activated carbons, for ethane and propane adsorption, RX 1.5 and CNR-115 showed superior adsorption capacities of the three gases because they possess not only high BET surface but also adequate volumes of micro-pores and mesopore. Hence, they have a good balance of different porosity to satisfy the low and high-pressure adsorption. This relationship was not noticed previously in the studies of smaller molecules because they can rotate and translate more freely compared to carbon chains.

The obtained data for pure gas adsorption isotherms were well fitted to the Langmuir two parameters model, indicating they are all microporous adsorbent and the surface coverage can be considered monolayer. In terms of the highest maximum Langmuir capacity, RX1.5, CNR-115, and ROx 0.8 showed the highest adsorption capacity of CH 4 , C 2 H 6 , and C 3 H 8 respectively compared to CGRAN and GAC 1240. The Langmuir adsorption capacity also was found to be influenced directly by the BET area of the adsorbent.

Figure 1 :

 1 Figure 1: Scheme of pure natural gas adsorption system.

Figure 2

 2 Figure 2 presents the morphology of five commercial activated carbons by SEM analysis. They possess two types of morphology: granular (CGRAN, GAC 1240, and ROx 0.8) or pellet (CNR 115 and RX 1.5). In terms of surface roughness, at the scale of 200 m, CGRAN has a hollow-channel form when all four other materials have a rough porous surface.

Figure 2 :

 2 Figure 2: SEM images of five studied activated carbons: a) CNR-115, b) RX 1.5, c) CGRAN, d) ROX 0.8 and e) GAC 1240 at 200 m scale.

Figure 3 : 8 Δ

 38 Figure 3: Adsorption isotherms of natural gas components at 303 K on five activated carbons: RX 1.5 (red), ROx 0.8 (green), CNR-115 (black), CGRAN (blue) and GAC 1240 (olive).

Figure 4 :

 4 Figure 4: Maximum Langmuir capacity at 303K as a function of the (a) BET surface areas, (b) micro-pores volumes, (c) meso-pores volume and (d) total pores volume.

3. 4 .

 4 Comparison with other activated carbons in literature 3.4.1. Methane adsorption Methane adsorption onto activated carbons is commonly published in the literature because this compound exists in both natural gas and raw biogas as a main component. The comparison of methane uptakes of several similar activated carbons (reported in Vargas et al., Tabatabaei Shirazani et al., and Peredo-Macilla et al.) are presented in Table 7

Table 1

 1 CAS registry number and purity of the gas cylinders.

	Component CAS Reg. No.	Purity
	Helium	7440-59-7	99.999%
	Methane	74-82-8	99.995%
	Ethane	74-84-0	99.95%
	Propane	74-98-6	99.95%

Table 2

 2 Natural gas compounds kinetic diameters and polynomial parameters[START_REF] Ismail | Gas Separation Membranes: Polymeric and Inorganic[END_REF][START_REF] Ettarh | Common descriptive and analytical statistics in investigative studies[END_REF][START_REF] Matteucci | Transport of Gases and Vapors in Glassy and Rubbery Polymers[END_REF].

	Parameter	Methane Ethane Propane
	Kinetic diameter (pm)	380	400	430
	a (⋅10 -10 )	2.10	2000	2158
	b (⋅10 -8 )	4.35	400	186
	c (⋅10 -4 )	649	4.00	6.66
	d	0.04	0.04	0.04
	e (⋅10 -4 )	0	3.00	0
	2	0.99	1	0.99

Table 3

 3 Experimental results of the adsorption of methane on five activated carbons at 303K

	RX 1.5	RX 0.8	CGRAN	CNR-115	GAC 1240

Table 4

 4 Experimental results of the adsorption of ethane on five activated carbons at 303K.

	RX 1.5	RX 0.8	CGRAN	CNR-115	GAC 1240
	(MPa) (		-1 ) (MPa) (		-1 ) (MPa) (		-1 ) (MPa) (		-1 ) (MPa) (	-1 )
	0.02	1.68	0.02	1.74	0.04	1.36	0.04	1.24	0.04	2.81
	0.09	3.84	0.10	3.75	0.12	2.63	0.11	2.71	0.27	5.31
	0.24	5.78	0.27	5.03	0.27	3.84	0.27	4.51	0.55	6.13
	0.56	7.38	0.48	5.69	0.58	5.25	0.58	6.26	0.86	6.54
	0.85	8.14	0.73	6.23	0.89	6.09	0.92	7.32	1.08	6.73
	1.49	8.95	1.22	6.88	1.55	7.18	1.49	8.40	1.90	6.89
	1.92	9.15	1.50	7.20	1.91	7.52	1.88	8.79
	2.41	9.13	1.90	7.59	2.38	7.71	2.42	8.94
			2.42	7.96

Table 5

 5 Experimental results of the adsorption of propane on five activated carbons at 303K.

	RX 1.5	RX 0.8	CGRAN	CNR-115	GAC 1240

Table 6

 6 Langmuir fitting parameters of adsorption isotherms of methane, ethane, propane on five commercial activated carbons at 303K.

		CH 4	1.09	10.79	0.008
	RX 1.5	C 2 H 6	0.18	10.29	0.016
		C 3 H 8	0.03	8.82	0.015
		CH 4	1.79	8.45	0.003
	CNR-115	C 2 H 6	0.38	10.98	0.015
		C 3 H 8	0.13	9.03	0.010
		CH 4	0.94	9.62	0.006
	ROx 0.8	C 2 H 6	0.20	11.62	0.024
		C 3 H 8	0.08	9.20	0.026
		CH 4	2.24	10.28	0.002
	CGRAN	C 2 H 6	0.38	9.31	0.016
		C 3 H 8	0.09	8.02	0.018
		CH 4	0.40	7.44	0.008
	GAC 1240	C 2 H 6	0.07	7.24	0.009
		C 3 H 8	0.05	7.75	0.021

Table 7

 7 Methane uptakes of RX 1.5 and CNR-115 in comparison with other adsorbents in literature.

	Sample	Reference Working conditions CH 4 adsorption capacity (	-1 )
	RX 1.5	This work	303K and 1 MPa	5.1
	CNR-115	This work	303K and 1 MPa	3.1
	AC-H 3 PO 4	[35]	298K and 1 MPa	3.0
	MP48	[68]	298 K and 1 MPa	2.4
	MCa2	[68]	298 K and 1 MPa	4.6
	HS-KOH1:2	[69]	298K and 1 MPa	5.2

Table 8

 8 Ethane uptakes of RX 1.5 and CNR-115 in comparison with other adsorbents in literature.

	Sample	Reference Working conditions C 2 H 6 adsorption capacity (	-1 )
	RX 1.5	This work 303K and 1.2 MPa	8.7
	CNR-115	This work 303K and 1.2 MPa	7.9
	C-700-3	[70]	298K and 1 MPa	7.2
	C-800-3	[70]	298K and 1 MPa	7.1
	C-PDA-3	[72]	298K and 1 MPa	6.6
	Ni(bdc)(ted) 0.5	[71]	298K and 1 MPa	5.0
	IRMOF-8	[73]	298K and 1 MPa	4.0

Table 9

 9 Propane uptakes of RX 1.5 and CNR-115 in comparison with other adsorbents in literature.

	Sample	Reference Working conditions C 3 H 8 adsorption capacity (	-1 )
	RX 1.5	This work 303K and 0.5 MPa	8.1
	CNR-115	This work 303K and 0.5 MPa	8.0
	Ecosorb	[74]	303K and 0.5 MPa	5.7
	Cu-MOF-74	[34]	303K and 0.5 MPa	5.6
	Co-MOF-74	[34]	303K and 0.5 MPa	5.4

Declaration of interest

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

Graphical Abstract

Determination of methane, ethane and propane on activated carbons by experimental pressure swing adsorption method

Bich Ngoc Ho,David Pino Perez,Deneb Peredo Mancilla,Joseph Diaz,Cécile Hort,David Bessieres,Camélia Matei Ghimbeu

Highlights

Determination of methane, ethane and propane on activated carbons by experimental pressure swing adsorption method

Bich Ngoc Ho,David Pino Perez,Deneb Peredo Mancilla,Joseph Diaz,Cécile Hort,David Bessieres,Camélia Matei Ghimbeu • The adsorption capacity of gas molecules was proved to depend on the material properties of each activated carbon.

• Specific surface areas play a major role in the scale of adsorption capacity.

• The influence of the carbon chain on the saturation of adsorption isotherms was verified.

• Adsorption isotherms of studied gases on commercial activated carbons were well fit by the modified Langmuir model.