A HOMOGENIZED MODEL ACCOUNTING FOR DISPERSION, INTERFACES AND SOURCE POINTS FOR TRANSIENT WAVES IN 1D PERIODIC MEDIA

. A homogenized model is proposed for linear waves in 1D microstructured media. It combines second-order asymptotic homogenization (to account for dispersion) and interface correctors (for transmission from or towards homogeneous media). A new bound on a second-order effective coefficient is proven, ensuring well-posedness of the homogenized model whatever the microstructure. Based on an analogy with existing enriched continua, the evolution equations are reformulated as a dispersive hyperbolic system. The efficiency of the model is illustrated via time-domain numerical simulations. An extension to Dirac source terms is also proposed.


Introduction
Understanding, modelling and controlling wave propagation in heterogeneous media is of major interest in numerous engineering domains, e.g.seismic simulation and protection, sound and vibration attenuation, nondestructive testing, etc.The particular case of periodically varying media, see e.g.[28] for sonic and phononic crystals, gathered much attention.Indeed the interaction between waves and the periodic medium produces pronounced dispersive features, including band-gaps ("forbidden" range of frequencies for which waves decay exponentially) at higher frequencies.On the modelling side, robust methods exploit the periodicity to establish effective models linking the salient macroscopic features of the wave propagation to the microstructure.These effective models enable efficient numerical simulation, and they also provide a solid ground to implement topological optimization algorithms of architected materials [2,14].
In this paper, the focus is on the long wavelength regime, compared to the characteristic size of the microstructure.In this regime (prior to the first band-gap), higher-order models accounting for dispersion effects are studied since [35], notably thanks to double-scale asymptotic homogenization methods [9], see e.g.[2,14] for models of scalar waves in two-dimensional media.One-dimensional situations are particularly well documented: [1,4,5,18,19,27] among other computed the higher-order homogenized coefficients and designed well-posed models thanks to the "Boussinesq trick" [3] that enables the partial or total permutation of space and time derivatives without loss of the asymptotic order of approximation.This trick was generalized in [41] that estab-lished a whole family of models.Introducing an additional modelling tool, [36] recently proposed a reformulation of these homogenized wave equations into hyperbolic systems well-suited for time-domain simulation [23], based on a phenomenological stress-gradient model presented in [21].
On the other hand, using these models in realistic bounded domains, i.e. addressing either boundary conditions or transmission conditions towards another homogeneous or microstructured medium, is an issue that is still a major research topic.In two or three dimensions, oscillating boundary layers that appear near these boundaries or interfaces, can be accounted for using corrector functions to complement the inner expansion and to obtain convergence estimates, see e.g.[34] for eigenvalue problems in bounded domains, [26] for elliptic problems with interfaces and [10,11,20] and the references therein for time-harmonic transmission problems.However, (i) these boundary or transmission correctors are complex objects whose theoretical study is quite involved [6,22] and (ii) their direct computation is often as costly as the one of the original problem posed on the microstructured domain.
Fewer works address the "practical" implementation of these correctors, i.e. the design of enriched boundary or transmission conditions suitable for numerical implementation and simulation.A notable exception is the contribution of S. Fliss and collaborators [8,20,40].These authors consider two-dimensional time-harmonic problems in half-planes or strips and address interface neighborhoods thanks to matched asymptotic expansions or "enriched" double-scale expansions accounting explicitly for local correctors.Also using matched expansions, effective transmission conditions are built for scattering problems by layered "slabs" of materials with interfaces either orthogonal [31] or aligned [32] with the layers, and in the context of water waves by [33].For purely onedimensional propagation, [15] proposes boundary and transmission conditions for wave problems homogenized up to second order, still for time-harmonic loads and fields.This methodology is also applied to non-uniformly oscillating media in [37].
Finally, the time-domain simulation of these homogenized wave models and associated boundary and transmission conditions is almost inexistant, to the best or our knowledge, although a proposal was made recently in [8,Chap. 5] for a half-plane with a Dirichlet boundary condition.In this context, the main purpose of the present paper is to build a well-posed "total" model to simulate 1D transient wave propagation through microstructured domains, combining (i) the family of second-order models from [41] and their reformulation as hyperbolic systems from [36], and (ii) the effective transmission conditions proposed in [15], adapted to the time-domain context to ensure well-posedness.An extension of the proposed tools to Dirac source points immersed in the microstructure is also given.
In Section 2, dispersive models for one-dimensional waves obtained from second-order homogenization are recalled, and a new result on the second-order homogenized coefficient is given (Prop.1).Then the proposed hyperbolic system is presented and the conditions for its well-posedness are established (Prop.2).Section 3 gives relevant corrections for transmission conditions between homogeneous and homogenized media, and states the stability conditions for the total model incorporating the hyperbolic system and these transmission conditions (Prop.3).Section 4 presents the treatment of source points.The improvements brought by the proposed models are illustrated by a set of time-domain simulations in Section 5, and Section 6 finally summarizes the findings of the papers and proposes several extensions.Auxiliary definitions and proofs are gathered in Appendix A.

Second-order homogenized model in unbounded space
This section focuses on free waves in an unbounded periodic elastic medium, characterized by density  ℓ () = (/ℓ) and Young's modulus  ℓ () = (/ℓ), in terms of 1-periodic functions (, ) and the periodicity length ℓ.The material displacement is denoted  ℓ (, ), and  ℓ =    ℓ and  ℓ are the velocity and stress fields.Combining the momentum balance equation   ( ℓ  ℓ ) =    ℓ and the elastic constitutive law  ℓ =  ℓ    ℓ , the source-free wave motion is then described equivalently by the wave equation satisfied by  ℓ , or the first-order system satisfied by ( ℓ ,  ℓ ): The main results obtained by second-order homogenization at low frequencies and wavenumbers, and moderate variations of material parameters (, ) (so that the ratio between period and wavelength is the only small scale of the problem), are summarized below.Since the derivation of these results is now classical, they are not justified: the reader is referred to [41] and the references therein for more details on the formal expansions, to [15] for error estimates in the time-harmonic case, and to [1,27] for error estimates in the transient case.
Remark 1. Rigorous asymptotic analysis should be done in terms of a non-dimensional parameter  = ℓ/, where  is a "reference" wavelength that serves as a spatial scale.For time-harmonic problems, this wavelength can be defined, see e.g.[15, Remark 1].In the present work, the focus is on modelling rather than on asymptotic analysis, and transient waves with broad frequency -and therefore wavelength-spectra are considered, hence all the expansions and estimates below are formally given in terms of ℓ for simplicity: (1), (ℓ) and (ℓ 2 ) remainders indicate leading-, first-and second-order approximations.Discussions on the choice of reference wavelength  and the influence of ratio  are postponed to numerical illustrations in Section 5.1.

Formal expansions, cell functions and macroscopic fields
The two-scale formal expansion begins with the following ansatz on the solution: where the spatial dependency is split between the "slow" space variable  and its "fast" counterpart  = /ℓ.Injecting this expression into the wave equation ( 1) and solving iteratively the resulting cascade of equations, see [41], one finds separated variable solutions of the form: written in terms of slowly varying macroscopic fields ū and fastly oscillating correctors involving 1-periodic cell functions   , modulated by the successive space derivatives of the previous macroscopic fields.The cell functions are solutions of static cell problems posed on the scaled periodicity cell [0,1], as recalled in Appendix A.1.
The formal process also provides wave equations satisfied by the macroscopic fields ū , whose source terms depend on the previously determined fields {ū  } < .Here we prefer to work with total macroscopic fields   instead, that concatenate all the contributions of order up to .In what follows, we accordingly define: Section 2.2 below provides and discusses the wave equations satisfied by these fields.
Finally, (i) defining the total macroscopic velocities   =     , and stresses   =  0     , where  0 is the effective Young's modulus to be defined later, (ii) summing the contributions of equation ( 2), (iii) using replacements e.g.ū0 =  1 + (1) in the corrector terms and (iv) neglecting higher-order contributions, the total th order approximations of ( ℓ ,  ℓ ,  ℓ ) are given by: where specific stress cell functions   are introduced [15], see Appendix A.1 and Remark 2 below.
Remark 2. The literature often focuses on the displacement  ℓ and on the wave equation it satisfies (left of ( 1)).The stress field  ℓ =  ℓ    ℓ , particularly relevant when considering boundary or transmission conditions, is then treated as a byproduct of the analysis, whose approximation accuracy suffers from the differentiation.On the other hand, the system featuring the two fields ( ℓ ,  ℓ ) in (1) emphasizes the similar roles played by these fields, that should be reflected by equally accurate approximations.This is the main motivation for the functions   associated to the stress field  ℓ , as proposed in this 1D waves context by [15].The two-fields formalism is also used for deriving boundary and transmission correctors [31][32][33] and providing associated error estimates in higher dimensions, see [10,11] and the references therein.

Homogenized wave models up to second order
Effective equations satisfied by the macroscopic fields (  ,   ,   ), also provided by the two-scale approach, are now discussed at orders 0, 1 and 2. Homogenized model at leading order.This model is given similarly to (1) by a wave equation on  0 or a system on ( 0 ,  0 ): The constant homogenized coefficients at leading order  0 and  0 are given in 1D by the usual formula: i.e. the arithmetic mean and the harmonic mean of their oscillating counterparts (, ).The celerity  0 is defined naturally as: First-order homogenized model.As shown by [13,Section 2.3] or [15, Lemma 1], the first-order macroscopic fields ( 1 ,  1 ,  1 ) also satisfy the homogenized equations (4), i.e. there is no first-order contribution to the wave equation.
Second-order homogenized model.The model ( 4) is valid in the quasistatic limit, i.e. for wavelengths much longer than ℓ.For shorter wavelengths, dispersive effects (i.e.frequency-dependent wave velocity) cannot be neglected and must be accounted for in the homogenized models.This is done by incorporating second-order terms into the wave equation, i.e. using the equations satisfied by ( 2 ,  2 ,  2 ).In [41], it was shown that second-order expansions lead to a family of enriched wave equations for the macroscopic displacement  2 : where the homogenization process provides the values of the homogenized velocity  0 given by ( 6), and the second-order coefficient , for which a new expression is provided here: Proposition 1.The second-order homogenized coefficient  can be computed in terms of the material coefficients (, ), the homogenized coefficients ( 0 ,  0 ) given by (5) and the first cell functions ( 1 ,  1 ): As a consequence,  is non-negative.
The proof of this result is given in Appendix A.2 and relies on reciprocity identities on the cell functions ( 1 ,  Remark 3. The family of models (7), parametrized by the coefficients (  ,   ,   ), is obtained thanks to the wave equation    2 =  2 0    2 + (ℓ 2 ) at leading-order, which enables to exchange second-order space and time derivatives in higher-order terms.This possibility, sometimes called "Boussinesq trick" [3], was used by [19] to exhibit a well-posed model with only "mixed" derivatives (  =   = 0,   = − < 0), to replace ill-posed models with fourth-order space derivatives (  =  > 0,   =   = 0) that arise naturally from the formal homogenization process (when  = 1).Later, [27] also studied the "mixed" model, derived an equivalent Korteweg-de-Vries equation and proved error estimates for long-time propagation.The study [1] also proved well-posedness and established error estimates for a two-parameter family of models (  −   = ,   = 0).The family (7) introduced in [41] is a generalization of these models which provides two degrees of freedom to choose a particular model.The model can therefore be optimized e.g. by fitting the dispersion curve of the periodic medium, see [36], while still satisfying the well-posedness requirements.
In higher dimensions, terms involving fourth-order space derivatives, associated with a high-order tensorvalued coefficient, are in general impossible to remove.However, the Boussinesq trick still can be performed to obtain well-posed models [3].
Finally, in this work as in [15], only the cases where   = 0 are considered because (i) a fourth-order space derivative in the wave equation would come along with additional boundary or transmission conditions that are not easy to define (although a proposal is made by [8,Chap. 5] for a Dirichlet boundary condition) and (ii) the analysis performed in [36] showed that this term leads to unstable models when   > 0, involving an additional modelling constraint.The so-called () model, featuring fourth-order mixed and time derivatives, is therefore used hereafter: and only one degree of freedom is left to choose the couple (  ,   ).
To investigate the properties of this wave equation, an equivalent system is now presented.Using this formalism enables to reuse theoretical results specific to hyperbolic systems (stability, existence of solutions), as well as dedicated numerical methods (integration schemes, discretization of interfaces).An overview may be found in the reference books [23,29].

Stress-gradient formalism
The upcoming system comes from a stress-gradient model introduced in [21] and studied in [36].This phenomenological model features four fields: the usual displacement  and stress , but also a microdisplacement  and a microstress  to account for microstructural effects.These fields satisfy the constitutive relations and equilibrium equations: written in terms of density and Young's modulus ( sg ,  sg ) and two micro-stiffness and micro-inertia coefficients (, ).All these parameters are strictly positive real numbers, and  should not be confused with a differentiation operator.
Relationships with the second-order homogenized model.Combining the relations (9), as done in [21, Section 5.3], provides a fourth-order enriched wave equation for : which is formally the same equation than the () homogenized equation given by (8).More precisely, the equivalence is obtained by identifying the coefficients ( sg ,  sg , , ) as: Moreover, identifying the terms in the two wave equations ( 8) and (10), and using the relations ( 9) and ( 11), the fields of the stress-gradient model satisfy the following relations with the macroscopic displacement  2 and stress  2 =  0    2 : Several remarks may be done from the above relations: -The microdisplacement  is a second-order term, while the microstress  is in fact a leading-order contribution: from (9) one has  =    + (1).-The parameter  controls the choice of () model:  = 1 corresponds to the "time" model () (with   = 0,   = −), and  = 0 corresponds to the "mixed" model () (with   = −,   = 0) in the wave equation (10).However  = 0 makes no sense in the original stress-gradient system (9) that therefore cannot describe the () model.
Equivalent system.To pursue towards the derivation of a hyperbolic system equivalent to the wave equation ( 8), the total macroscopic velocity , incorporating the microscopic fluctuations, and the microscopic velocity  = (ℓ) are introduced: Then, slightly reformulating the equations ( 9), the four fields (, , , ) are found to satisfy the system: where the convenient parameter  controls the choice of () model: In [36] a similar system was introduced, featuring the classical velocity  =    instead of .The present choice slightly simplifies the notation when introducing interfaces in Section 3.
Total fields approximations.Finally, approximations of the total fields ( ℓ ,  ℓ ) are given by the expansions (3) in terms of the macroscopic fields ( 2 ,  2 ).From the relations ( 12), these macroscopic fields can be recovered a posteriori from the stress-gradient velocity and stress (, ) and the auxiliary fields (, ) as: However, a simpler way to compute the total fields, without computing the fields ( 2 ,  2 ), is to use the following approximations: They are formally justified by the second-order amplitude of the fluctuations  2 −  and  2 − .

Model properties
The system (14) above can be recast in matricial form: where  := (, , , ) T and: Hyperbolicity.The matrices  and  have eigenvalues: The eigenspaces of the matrix  are: For  ̸ = 0, all these vectors are linearly independent and therefore they form a basis of R 4 .To obtain a hyperbolic system, the eigenvalues of  need to be real.To ensure that the null solution is stable, the eigenvalues of  need to be imaginary.These two conditions imply that the coefficients (,   ,   ) should satisfy: {︂ Hyperbolicity:  > 0 Stability: As already pointed out by [36] and omitted in previous studies [15], these conditions complement the condition −  −   =  coming from the homogenization process.For latter use, the following exponential matrix, needed in the numerical integration of ( 17), is introduced: where: Energy conservation.Introducing the symmetrizer matrix one obtains: i.e.  •  is symmetric, and  •  is skew-symmetric.Multiplying the system (17) by  T •  and making use of these (skew)symmetries, one obtains: Considering compactly supported initial conditions or source terms, and integrating over a space-time domain Ω × [0,  ], where Ω is chosen large enough so that  (•,  ) = 0 on Ω, one obtains the energy conservation: The volume energy ℰ vol is the sum of a kinetic and elastic energy given as: In the second expression, the macroscopic velocity  =    is introduced and 1 −  = /  is used to emphasize that the energy is positive since  > 0 (Prop.1) and   > 0 (as required by the stability condition ( 18)).
Finally, the results of this part are summarized in the following proposition: Proposition 2. If the parameters (  ,   ) satisfying the relation −  −   =  are chosen so that   < 0 and   > 0, the system (14) is hyperbolic, the null solution is stable, and the associated positive volume energy ℰ vol defined by (21) is conserved.

First-order transmission conditions and total model
Now that an effective model is proposed for unbounded microstructured media, this section focuses on how to use this model for bounded domains made of microstructured materials.In the same way as higher-order correctors were introduced in the wave equation, higher-order interface correctors will be introduced.These correctors will be defined as in the previous study [15], that focused on time-harmonic problems and in which rigorous error estimates are proven to support this proposition, following earlier works on boundary correctors, see [10,11].In the present work, an extension is made to the transient case without rigorous justification: our purpose is to establish a well-posed and practical model rather than to prove asymptotic accuracy (that is nevertheless observed in numerical examples).
To fix ideas, the transmission of waves from a homogeneous domain Ω − = { < 0} to a microstructured domain Ω + = { > 0} is first addressed in detail.Then the extension to a "slab" Ω =]0, [ surrounded at both extremities by homogeneous media is given.
For future use, let us introduce the jump  0 and the mean ⟨ ⟩ 0 of a function  across an interface at  = 0 as: and we give the following relations:

Microstructured transmission problem
The original problem to be approximated is defined as follows: the fields ( ℓ ,  ℓ ) satisfy the system (1) with constant coefficients ( − ,  − ) in Ω − = { < 0} and with oscillating coefficients ( ℓ ,  ℓ ) in Ω + = { > 0}.A perfect interface is considered at  = 0, i.e. the velocity and stress are continuous: To complete this problem, initial conditions are finally given: where  ⋆ and  ⋆ are entirely supported in the homogeneous domain Ω − to simplify the analysis.Indeed, initial conditions supported in the microstructured domain Ω + should be addressed in the homogenization process, which is outside the scope of the present work: again the reader is referred to [27] for insights on these issues in the case of an unbounded domain.In Section 4, we will study the case of a Dirac source point immersed in the microstructured medium.
In the following parts, the focus will be on the interface conditions at  = 0 when the wave propagation in the right domain Ω + is described by an effective model.

The homogenized model at leading order
The model at leading order is first presented.The macroscopic fields are noted ( 0 ,  0 ) and the correctors in the approximation (3) are ignored, i.e. the approximations ( ℓ ,  ℓ ) ≈ ( 0 ,  0 ) are used.Using the effective model at the leading order (4) in Ω + , these fields satisfy (i) the systems: (ii) the interface conditions at  = 0 inherited from the conditions (23) (see [10] and [32, Appendix 1.A.]): 0 0 = 0, and  0 0 = 0, (26) and (iii) the initial conditions ( 24) with ( ℓ ,  ℓ ) replaced by ( 0 ,  0 ).To prepare the ensuing analysis for higher-order models, the energy associated with this problem is now studied.Similarly to what is done in Section 2, one obtains from the systems (25): where the volume energy ℰ vol 0 is now: and the interface term is: Then, a sufficient condition for the stability of the problem is this term to be the time derivative of a positive interface energy: if there exists ℰ int 0 > 0 such that  int 0 =   ℰ int 0 , then (i) the total energy is conserved from ( 27) and (ii) the the volume term ℰ vol 0 is bounded: because the initial conditions on ( 0 ,  0 ) ensure ℰ int 0 (0) = 0 and ℰ int 0 ≥ 0 by assumption.This property then ensures that the fields ( 0 ,  0 ) remain bounded in time.
In the leading-order case, the interface term  int 0 is identically canceled by the transmission conditions ( 26), and therefore from (27) the energy ℰ vol 0 is conserved.

First-order homogenized model: correcting the transmission conditions
First-order approximations from (3) are now considered: where ( 1 ,  1 ) = (0, 0) in the homogeneous domain Ω − .As seen in Section 2.2, the first-order macroscopic fields ( 1 ,  1 ) satisfy the same non-dispersive systems (25) as at leading order in each of the two half-spaces.On the other hand, replacing the fields ( ℓ ,  ℓ ) by their approximations (30) above in the transmission conditions (23) and neglecting the remainders, as proposed by [15], one obtains imperfect non-symmetrical conditions on the jump of the macroscopic fields: i.e. a first-order correction compared to the conditions at the leading order (26).
A similar energy analysis than at the leading order results in where the total volume energy ℰ vol 1 is the same as ℰ vol 0 given by ( 28) with ( 0 ,  0 ) replaced by ( 1 ,  1 ).The same analogy holds for the interface term  int 1 compared with  int 0 in (29).Given the conditions (31), one has: There is no guarantee that  int 1 is the derivative of a positive interface energy.The negative case may result in an ill-posed problem.The following paragraphs explain how to deal with this issue, through symmetrization and stabilization of the jump conditions.
Spring-mass transmission conditions.To symmetrize the relations (31), the limit values e.g.   + 1 need to be written in terms of mean values across the interfaces, e.g.⟨   1 ⟩ 0 .The process introduces additional approximations, up to residuals that should be at least of second order (i.e. in (ℓ)) to preserve the overall first-order approximation.The main tools to do so are the bulk equations that link the time and space derivatives of ( 1 ,  1 ) in Ω + : Then, using also the jump-mean relations (22), the following approximations are found at the leading order: , where the jump conditions (31) were used to neglect the first-order jumps  1 0 and  1 0 .Combining these approximations with (31) and neglecting second-order terms, one finally obtains: This form is better known as spring-mass transmission conditions, that may also model e.g. a thin layer of elastic material between the considered media [30,App. 1].A new interface term  int 1,sm then replaces  int 1 in ( 32)- (33), conveniently written as the derivative of an interface energy ℰ int 1,sm : However, there is still no guarantee that this interface energy remains positive for all time, except in the particular cases where  1 (0) ≤ 0 and  1 (0) ≤ 0. It is why a last modelling step is necessary, as addressed now.
Stabilization introducing an enlarged interphase.To obtain centered transmission conditions that are associated to a positive interface energy, a thick interface   = [−ℓ, ℓ] of thickness 2ℓ, centered on the point  = 0, is introduced, and the transmission conditions are written across this interface.This method has been first used for microstructured interface homogenization [17], but also for a 2D transmission problem similar to the present case in [20,40].Using Taylor expansions of the macroscopic fields in the homogeneous and homogenized media, e.g. 1 (0 − ) =  1 (−ℓ) + ℓ    1 (−ℓ) + (ℓ) and keeping only first-order terms in ℓ, the conditions (31) become: Then, using the same first-order approximations (35) than above, one obtains equivalent (up to (ℓ)) spring-mass conditions: where •  and ⟨•⟩  are the jump and mean values of fields across the interface   .Similarly to (37), the interface term in the energy is then written: A sufficient condition for the interface energy ℰ int 1, to be positive is that the two coefficients ( 1 ,  1 ) are positive.This is achieved by choosing the interface parameter  such that: As already noticed, when  1 (0) ≤ 0 and  1 (0) ≤ 0, there is no need for a thick interface: one can choose  =  min = 0 and the conditions (36) are retrieved.In the other cases, choosing  =  min will cancel one of the factors ( 1 ,  1 ) in the relations (39), i.e. one has either  1  = 0 or  1  = 0.
Remark 5. Instead of introducing a thick interface, another way to design stable interface conditions from (36) is to slightly modify the homogenized model by choosing non-zero means to the cell functions ( 1 ,  1 ) to ensure  1 (0) ≤ 0 and  1 (0) ≤ 0. This method was introduced in [1] to design stable wave models in unbounded media, and applied in [8] to stabilize enriched boundary conditions at the edge of a half-plane.It was implemented by the authors of the present paper, and is indeed a good alternative to treat one interface.But since it requires a global change of the model (in the whole microstructured domain), it cannot be easily extended of the case of a slab treated below in Section 3.5, or any case where more than one interface needs to be addressed.The thick interface method, that introduces a local approximation, was therefore preferred.
Remark 6.When  =  0 is constant, i.e. when  is the only varying parameter, one has  1 = 0 (see Appendix A.1) and therefore there is no first-order correction on the velocity jump in the spring-mass conditions (36), but still a correction on the stress.This was shown already in [32] for electromagnetic waves ( being the electric field and  being the permittivity).However, depending on the sign of  1 (0), stabilizing the interface condition may still be necessary, and artificial corrections on the velocity jump may appear e.g. if the "thick" interface method is chosen, resulting in conditions (39).

Total "hybrid" model incorporating second-order dispersion
Finally, the second-order model studied in Section 2 is implemented in the microstructured domain to account for dispersion.Using the four fields (, , , ) of the stress-gradient system in Ω + , while ( ℓ ,  ℓ ) ≈ (, ) in Ω − , these fields satisfy: At this stage, an important remark is that the two auxiliary fields (, ) satisfy ordinary differential equations in time, and therefore no boundary condition are needed for these fields (but initial conditions must be added).The systems (41) must therefore be complemented by transmission conditions on (, ) only, as in the leading-order and first-order cases, and by initial conditions for all fields.Since we choose initial conditions only supported in the left domain Ω − , we set: (, 0) =  ⋆ (), (, 0) =  ⋆ (), (, 0) = 0 and (, 0) = 0. (42) Then, combining the results on energies of the stress-gradient system collected in Section 2.3, and the additional term coming from the interface as specified above, one obtains: where the total energy ℰ vol is now: and the interface term  int is again given by: This term is formally identical to the ones that appear for the first-order model (33).
At this stage, it would be natural to use the second-order approximations in (3) to design transmission conditions for (, ).However, for technical reasons discussed in Remark 7 below, only first-order transmission conditions are designed to obtain a "hybrid" model, using the terminology introduced in [32].
As in the previous section, the chosen total fields approximations are: and the same analysis can be performed.Indeed, the links between the fields derivatives given by ( 34) can be replaced by: The second relation in ( 45) is obtained by taking the difference between the first and third equation of (41) in Ω + , and by using  = (ℓ).The second-order remainder (ℓ) in (45) has no influence on the approximations at the leading order (35).Consequently, the same first-order spring-mass transmission conditions (39), written across a thick interface, are applied to (, ): The various asymptotic formal expansions are still valid and indicate that the transmission conditions (46) should result in an overall first-order approximation.In this way, the overall problem including (i) the systems (41), (ii) the initial conditions (42) and (iii) the transmission conditions (46), is proven to be stable using the same arguments than in the previous section.It is summarized in the next proposition.
Proposition 3. We consider as "total model" the system (41) with the symmetrized jump conditions (46) across a thick interface, where the interface parameter  satisfies (40).Associated with this system is a constant energy over time ℰ = ℰ vol + ℰ int  , where ℰ vol ≥ 0 is given by (43), and ℰ int  is given by Remark 7. Imposing transmission conditions to the second-order approximation ( 16) of the total fields, as done in the time-harmonic case in [15], would lead to where  ± = (0 ± , ) and similarly for (, , ) and their derivatives.The same ideas than in Section 3.3 can be applied to obtain symmetrized jump conditions on a thick interface of the form: with ( 1 ,  1 ) given by ( 46) and new coefficients ( 2 ,  2 ), generalizing (46).However, these second-order terms feature (i) (non-classical) second-order time derivatives of (, ) and (ii) the traces of auxiliary fields (,   ), and our attempts to prove the stability of the resulting model were unsuccessful so far.Moreover, the "hybrid" model described in Proposition 3 already provides a good compromise between implementation easiness (classical spring-mass transmission conditions are used) and both qualitative and quantitative performances, as illustrated in Section 5 below.

Microstructured slab bounded by two homogeneous domains
To end this section, one extends the previous model to the case of a microstructured medium bounded at both extremities by homogeneous domains.It is representative of many real experiments of wave propagation.The microstructured domain is Ω  =]0, [, while Ω + = { > } now denotes the right homogeneous domain, characterized by coefficients ( + ,  + ).Additional transmission conditions must be designed for the second interface at  = .The continuity of the first-order approximations (44) is written: where   = /ℓ.Then, a similar analysis as in Section 3.3 leads to a reformulation of these conditions involving the jumps •  ′ and means ⟨•⟩  ′ across a thick interface [ −  ′ ℓ,  +  ′ ℓ] as: where the two coefficients ( ′ 1 ,  ′ 1 ) must be positive for these conditions to be associated with a positive interface energy.The thickness  ′ of this second interface must be chosen so that: a condition similar to its counterpart (40) except for the signs.

Correctors for a Dirac source term
So far, no source term was considered in the evolution equations.The initial conditions were supported in the homogeneous domain, thus avoiding to adapt them to the microstructure [27].To complete the modelling framework, the case of a single force source point immersed in the microstructured medium is now addressed, using the tools previously developed for interfaces.

Source terms and jump conditions
A point force at point  s is modeled by a Dirac source term in the initial wave equation: where  is the Dirac distribution and  the time-dependent amplitude of the source.To simplify the analysis, and having in mind the microstructures with fixed periodicity length ℓ considered in the numerical examples, we furthermore assume that the point source position is fixed at  s ∈ [0, 1] inside a unit cell, and we note  s = ( s ) the density at this position, that does not depend on ℓ.
The above equation ( 49) is to be understood in the weak sense, i.e. : ∫︁ with a suitable choice of domain Ω =] − ,  + [ and functional space (Ω).As in the previous sections, space and time dependency are omitted in (50) and below for compactness.Such a source point has been addressed by [13], and even extended to double-couple source (modelling ruptures of seismic faults) and to higher dimensions (for non-periodic homogenization) in [12].However, the provided expression, relying on energy conservation arguments, only treat first-order homogenization.The present proposal exploits the fact that in 1D, a Dirac distribution yields a jump condition on the solution.As a consequence, the formalism developed in the previous sections can be reused here to complete the two-scale formal expansion.
Remark 8.This approach does not hold in higher dimensions, for which an analysis of the correction for a source point should be incorporated into the two-scale asymptotic expansion, see [12].However, it could be applied to line sources in 2D or area sources in 3D, modelling e.g. a distribution of charges in electromagnetism, arrays of microphones in acoustics etc.
Replacing the source point by a jump condition.In 1D, a Dirac source term entails a discontinuity in the stress  ℓ =  ℓ    ℓ (accordingly with the physical interpretation of a Dirac source term as a point force).Indeed, integrating the wave equation (49) between  s −  and  s + , one obtains: and then letting  → 0 gives the announced result.Accordingly, the wave equation ( 49) may be rewritten as the system: where • s denotes the jump at source location  s .This problem can be homogenized as the previously seen transmission problem, except that the same periodic (or homogenized) medium occupies both half-spaces, and the successive corrections will account for the discontinuity of  ℓ .

Families of models
Leading-order homogenization.At the leading order, as in Section 3.2,  ℓ and  ℓ are replaced by the constant  0 and  0 in the bulk equations of system (51), and the source term is unchanged (still featuring the value  s ).
First-order homogenization.At first order, as in Section 3.3, the bulk equations satisfied by the macroscopic fields denoted ( 1 ,  1 ) are the same as at leading order, and the transmission conditions are applied to the first-order approximations of the full fields (30), leading to:  <  s and  >  s : where  s =  s /ℓ is the source position inside the cell that contains  s .Then, the bulk equations provide approximations at the leading order of the jumps of spatial derivatives: Neglecting the (ℓ) terms in (52), one obtains first-order conditions: i.e. the system in (52) features an additional inhomogeneous jump condition on the macroscopic velocity.The resulting system may finally be written in the equivalent form: This expression was also obtained in [13, Section 2.4] using an energy conservation argument.
Application to the stress-gradient system.Finally, the stress-gradient system (20) satisfied by the macroscopic fields (, , , ) is considered.As in Section 3.4, only first-order corrections are applied to the source term for simplicity.Hence, the analysis performed above for these corrections remains valid: the main difference is that the stress-velocity relation    =  0    + (ℓ) is a first-order approximation rather than an equality, but it does not change the approximations at the leading order (53).Therefore jump conditions (54) are applied to (, ).
Coming back to Dirac notations, the stress spatial derivative    is involved in two equations of the system, and hence the source term must be accordingly distributed: As in Section 3.4, this is a hybrid model combining first-order corrections for the source term and second-order corrections in the wave equation to account for dispersive effects.If transmission correctors were derived at the second order, see Remark 7, applying them to the problem (51) would certainly enable to obtain a complete second-order model.

Numerical experiments
The proposed models are now applied to simulate the wave propagation into an example bilaminate material.The material properties, initial conditions and source terms are first given (Sect.5.1), the numerical implementation is briefly described (Sect.5.2), and two test-cases are addressed: (i) a wave transmitted from a homogeneous medium towards the bilaminate (Sect.5.3) and (ii) a wave emitted from a source point embedded in a slab of bilaminate, surrounded by homogeneous media (Sect.5.4).

Microstructured medium, initial conditions and source modelling
The numerical illustrations are provided for bilaminate materials, for which analytical expressions of the homogenized coefficients and cell problems are available, see Appendix A.3.Furthermore, an analytical formula for (  ,   ) (or, here, ) was proposed in [15] to achieve a better fit of the exact dispersion curve for low frequencies, see (A.10).This specific choice of () model was also discussed in [36] and found to compete reasonably with numerically optimized models that achieve an overall fit of this curve on the whole Brillouin zone.Finally, this choice is numerically found to be compatible with the stability conditions {  < 0,   > 0} established by (18).The bilaminate investigated here has a periodicity length ℓ = 20 m, with a phase ratio  = 0.25 (Fig. 1).The physical parameters are   = 1000 kg/m 3 ,   = 10 9 Pa, and   = 1500 kg/m 3 ,   = 6 10 9 Pa.The associated homogenized parameters are  0 = 1375 kg/m 3 ,  0 = 2.66 10 9 Pa and  = 0.01959, see Appendix A.3.The optimized parameters of the () model (A.10) are   = −0.107071and   = 0.087481.
In the following examples, two types of forcing are considered: (i) a Cauchy problem (24) with an initial right-going wave or (ii) null initial conditions with a time-dependent forcing at a Dirac source point (Sect.4.2).Both configurations involve the source function (), chosen as the following combination of sinusoids with bounded support: where   = 2 −1 , the coefficients   are  1 = 1,  2 = −21/32,  3 = 63/768,  4 = −1/512, and the amplitude factor is  ≡ 1, unless otherwise stated.It entails that  is a smooth function , see the display in Figure 2(a).
Moreover, () is a wide-band signal with a central frequency  c =  c /2.As announced in Remark 1, this frequency is used to define a dimensionless parameter characteristic of the "low-frequency" regime underlying the formal asymptotic expansions used throughout the paper: We denote ĝ() the Fourier transform in time of (), where  = 2 and  is the frequency.Figure 2(b) displays the normalized modulus |ĝ|() in terms of  = ℓ  / 0 = ℓ /(2 0 ), for fours values of  c .Note that even for the smallest central frequencies, non-negligible values of  are reached by the signal.

Numerical set-up
In the following experiments, comparisons are made between full-field simulations in the microstructured configurations and simulations involving their homogenized counterparts at different orders.All the numerical simulations are done on a uniform grid, with space and time discretization parameters ∆ and ∆.The discretization is fine enough to guarantee that the numerical artifacts are negligible.A much coarser grid can be used for simulations in homogenised media than for simulations in microstructured media (this is one of the great interests of homogenisation!).Nevertheless, as this is not the object of the present work, we use here the same grid ∆ = 1 m (to be compared with cell size ℓ = 20 m) for both types of configurations.
The reference simulations in the microstructured configurations are done using a fourth-order ADER scheme [29].This explicit two-step finite-volume scheme has a stencil of width 2. In the homogenized configurations, we prefer to use the standard Lax-Wendroff scheme, of width 1. Doing so ensures that the stencil does not cross the whole enlarged interphase, and only one ghost value is required in the interphase [39].Both ADER and Lax-Wendroff schemes are stable under the usual CFL condition  = max()∆/∆ ≤ 1; in practice,  = 0.95 is chosen.
In the case of second-order homogenization, the stress-gradient formulation yields a non-homogeneous firstorder system (17).A splitting strategy is then followed, where one solves alternatively the homogeneous part by the Lax-Wendroff scheme, and the relaxation part exactly thanks to the exponential (19).The reader is referred to [7] for details.
The discretization of the transmission conditions relies on the Explicit Simplified Interface Method.Roughly speaking, the numerical scheme is modified at grid points surrounding the interfaces, based on the jump conditions (31) or (46).Once again, the reader is referred to [30] for details.Lastly, the numerical approximations of the fields are denoted as follows: - ℎ for the microstructured velocity  ℓ (1); - 0 for the zeroth-order homogenized velocity (25) with the transmission conditions (26); - 1 for the first-order homogenized velocity (31) without (ℓ), endowed with the conditions (31); - 2 for the "total" homogenized velocity (44) without (ℓ), endowed with the conditions (46).

Transmission towards a microstructured half-space
A bounded domain [−400, 1100] is considered.An interface at  0 = 350 m separates the domain into two areas: (i) on the left, a homogeneous medium; (ii) on the right, the bilaminate described above.A right-going wave is initially supported in the homogeneous medium, with a time shift  0 = 0.18 s in (57).
To highlight the importance of correctors for the transmission conditions, the physical parameters of the homogeneous medium are chosen as  − =  0 and  − =  0 , to ensure impedance matching ( −  − =  0  0 ) between the left and right media.At leading order and first order, the whole domain has thus homogeneous physical parameters, and the (first-order) transmission condition is the only signature of the microstructure.
Figure 3 illustrates the low-frequency configuration:  c = 3 Hz ( c = 0.043).The wave is shown at the initial time (a) and at  = 0.5 (b), where a reflected wave is observed in the Cauchy medium.For the model at leading order (Sect.3.2), the medium appears as homogeneous: no reflected wave is generated (c), which clearly shows the inability of this naive approach to handle the boundary effects of a homogenized medium.Moreover, the small-scale gradients of the solution in the microstructured medium are not captured.The first-order model (Sect.3.3) fixes these deficiencies: the reflected wave is captured, and the correctors handle the cell evolution of the transmitted wave (d) (note that at the scale of the figure, the two interfaces cannot be distinguished from each other).The same observations can be done for the total model (Sect.3.4): at low frequency, the latter captures equally well the wave phenomena (e).A close-up on the wave in the transmitted medium shows that the mean field  (13) captures the slow variations of the solution, whereas the correctors in (15) describe finely the gradients of the solution (f).
A higher value of the central frequency is considered in Figure 4:  c = 9 Hz ( c = 0.123).Compared with the low-frequency case, the microstructured field  ℎ exhibits a dispersive behavior.The model at leading order does not capture any phenomena: no reflected wave, no dispersion, no small-scale gradient (a).The first-order model captures two of the phenomena (reflected wave and gradients), but the dispersion in the right domain cannot be handled (b): this model is intrinsically non-dispersive.On the contrary, the total model captures all these features (c).A close-up in the right subdomain illustrates the quality of the correctors of the transmitted wave (d).The agreement is slightly worse on the main front of the reflected wave.We do not know how to explain this difference in behaviour, which is also observed in other simulations not presented here.A loss of accuracy is also observed behind the main front of the reflected and transmitted waves.These oscillations are associated with high-frequency components that are not captured by the dispersive model.
Lastly, Figure 5 illustrates the high central frequency  c = 12 Hz ( c = 0.172), where dispersion becomes the predominant feature.The limitations of models at the leading order and first order are examplified (a and b).Even for this large  c , the total model captures finely the main fronts of the transmitted dispersive wave (c and d).The amplitude errors are larger on the reflected wave.
Other simulations have been done in cases without impedance matching (not shown here).Similar conclusions are obtained: (i) the model at leading order yields a large error in the reflected wave; (ii) the first-order model captures well the reflected wave and the small-scale variations of the solution, but it is restricted to low frequencies; (iii) the total model captures all the features, and maintains accuracy even for large values of  c .
The accuracy of the effective models is then examined quantitatively.A source point and a receiver are respectively located at  s <  0 and  r >  0 .We impose that these positions correspond to nodes of the grid to avoid discretization effects.Moreover, the relative position of the receiver with respect to the microstructure is kept constant.Finally, the distances , the measurement time  and the forcing amplitude  must be kept constant in scaled variables (denoted by bars)  =  ,  =   , and  = / 2 , with  = /.To satisfy these criteria, the indices  and  are introduced and one defines Err Figure 6(a) represents the error on velocity fields {  } =0,1,2 as a function of the asymptotic parameter  c .The expected orders of convergence are asymptotically retrieved, in particular the "total" model remains of first order due to the choice of total field and associated transmission conditions.However, accounting for dispersion using the stress-gradient model brings a clear improvement: additionally to the qualitative gains observed in  Figures 4 and 5, the prefactor in front of the error on  2 is divided by approximately 5 compared to that of the error on  1 .With reference to Remark 2, the error on stress fields is also plotted in Figure 6(b), and the same observations can be made, highlighting the identical treatment given to the velocity and stress fields.We also examine the of larger values of the interface thickness on the accuracy for two central frequencies (in previous experiments, the minimal value 2ℓ = 2 min ℓ ≈ 0.68 m was chosen).To do so, we consider integer values of ℓ, from which are deduced the interface parameters  1 and  1 (39).Figure 7(a) shows relative error as a function of the interphase thickness for  1 and  2 .In both cases, the error grows strictly with ℓ.However, this growth is moderate, which allows enlarge the interphase while keeping an almost constant accuracy.One can therefore choose a given thickness from numerical requirements (e.g.mesh size or needed "ghost" nodes inside the interface) without deteriorating the solution too much.Since we do not have such requirements with the chosen numerical methods, the minimum value  =  min is kept in the following numerical experiments.
Finally, to illustrate what may happen when the sufficient conditions of stability on the transmission condition coefficients are not satisfied, Figure 7(b) represents a left-going wave impacting a flat interface between a homogeneous and a homogenized media occupying the right and left half-spaces, respectively, i.e. the configuration described in Section 3.5.The interface is not thickened, i.e.  ′ = 0 in (47), which results in  ′ 1 = −7.81× 10 −11 < 0 and  ′ 1 = 46.875.When the wave crosses the interface, an instability appears  and quickly grows: here it reaches the same amplitude than the incident wave in a few milliseconds (while the characteristic times of propagation are in seconds in these numerical experiments).

Source point in a slab
This second example illustrates both the modelling of a Dirac source term and of a bounded laminated medium.A domain [−300, 1700] is divided into three parts: a heterogeneous slab of bilaminate [350, 1090] is surrounded by two homogeneous Cauchy media with properties  ± = 1000 kg/m 3 and  ± = 2.25 10 9 Pa.Contrary to the previous example, is no impedance matching between the homogeneous and homogenized Figure 8 shows snapshots of the microstructured solution  ℎ at several times, for various central frequencies.For the sake of clarity, the internal interfaces of the bilaminate are not shown.As time increases, one observes the waves emitted by the source point and scattered by the external interfaces, leading to transmitted waves (in the  homogeneous external medial) and reflected waves (in the internal microstructured medium).At low frequency  c = 3 Hz (Fig. 8(a)), the dispersion is negligeable.The only small scale variations are due to the internal microstructure.As the central frequency increases ((b)-(d)), the dispersive effects are becoming increasingly important.A slight asymmetry between left and right parts of the figures may be observed.This is due to the grid-induced assymmetry of the Dirac source point.
Figure 9 compare the microstructured velocity  ℎ with the velocity  2 of the total model (Sect.3.4).Once again, the internal interfaces are not shown.The velocities of the models at the leading order and first order are discarded, since they are insufficient to capture accurately dispersion.Figure 9(a) and (b) illustrates the low-frequency forcing  c = 3 Hz.The fields are shown at  = 0.25 s (a), where the emitted waves are still in the slab, and  = 0.6 s (b), where waves have been transmitted in the surrounding media.As in Figure 5, dispersive effects are not visible at this frequency.Figure 9(c) and (d) illustrate the case of a higher frequency  c = 9 Hz, where the dispersive effects become more visible in the microstructured model.The total model captures the first wavefronts and their transmission to the surrounding media.Figure 9(e) and (f) display the high-frequency case  c = 12 Hz.Some inaccuracies of the total model can be observed, especially inside the slab, but the qualitative agreement with the microstructured solution remains good, even for this large value of the parameter  c = 0.172.
For completeness, the stress fields are also displayed in Figure 10 for  c = 9 Hz corresponding to the intermediate configuration of Figure 9(c) and (d), and for thickened interfaces:  =  ′ = 4 m, so that  1 = 2.28 × 10 −10 Pa −1 ,  1 = 503.12Pa (left interface),  ′ 1 = 8.57 × 10 −11 Pa −1 and  ′ 1 = 521.87Pa (right interface).Once again, and despite these non-optimal values of  and  ′ , excellent agreement is observed between the microstructured and "total" fields, especially for the first wavefronts of the fields transmitted into the left and right homogeneous domains.
Finally, to provide another representation of these approximations, a receiver put at  r = 340 m records the field transmitted to the left homogeneous domain.Figure 11 shows the time history of the absolute errors | ℎ ( r , ) −   ( r , )|.The signal is null up to  ≈ 0.3 s, which corresponds of the travel time from the source point to the receiver.A low frequency (a), the error of the model at leading order is much bigger than that of the first-order model and total model.It emphasizes the role of the first-order transmission conditions.At higher frequencies (b)-(d), the models at leading order and first order are unable to capture the dispersive effects, which yields similar errors.On the contrary, the benefit induced by the total model is clearly observed.Lastly, the maximal amplitude of the error for each model increases with the frequency.

Conclusion and perspectives
In this paper, a consistant "hybrid" model was proposed for transient waves in periodic media, combining second-order correctors of the wave equation to account for dispersion and first-order correctors for transmission conditions from or toward homogeneous domains.Based on a reformulation as a hyperbolic system, its wellposedness was proven, and its efficiency was established through numerical simulations.
Effective coefficients.The homogenization indexes at leading order are first defined as: First cell problem.As already seen,  0 = 1 and the first cell function  1 is given by: in terms of the index   given by (A.8).
Second cell problem.The cell function  1 is: ]︂ .

Figure 1 .
Figure 1.Unit cell of a bilaminate with two phases (, ) and phase ratio .

Figure 2 .
Figure 2. Time evolution of the source (58) at  c = 3 Hz (a).Spectrum of the source for various central frequencies  c , as a function of the asymptotic parameter  c (b).

Figure 6 .
Figure 6.Interface between a homogeneous half-space and a microstructured half-space.Error measurements for the velocity (a) and stress (b) fields.Dashed lines indicate expected asymptotic slopes.

Figure 7 .
Figure 7. Interfaces between a homogeneous half-space and a microstructured half-space.(a) Error measurements on velocity fields  1 and  2 for variable interface thickness and two frequencies for the configuration of the previous figures.The minimum thickness 2 min ℓ = 0.68 m is indicated by the dashed vertical line.(b) Example of instability occurring for a leftgoing wave impacting a flat interface (dashed vertical line).In this example the homogeneous and homogenized media occupy the right and left half-spaces, respectively.