Raffaella Rizzoni
email: raffaella.rizzoni@unife.it

Serge Dumont

Frédéric Lebon

Elio Sacco

Higher order adhesive effects in composite beams

. Accounting for higher order terms of the asymptotic expansions in the adhesive, the proposed approach generalizes simpler models based on the classical spring-type interface law or on the case of perfect contact between the adherent layers.

The proposed methodology is used to evaluate stresses in two adhesive bonded joint configurations subjected to bending moment and transverse shear loading. Numerical simulations are produced and the results show good agreements with those obtained through finite element analysis.

Introduction

Structural bonding has become an essential technological solution in mechanical engineering. Whether in aeronautics or civil engineering [START_REF] Balakrishnan | Potential repair techniques for automotive composites: a review[END_REF][START_REF] Birman | Review of current trends in research and applications of sandwich structures[END_REF][START_REF] Shang | Review on techniques to improve the strength of adhesive joints with composite adherends[END_REF], this solution is becoming the preferred choice because of the desire to lighten structures in order to save energy and simplify handling. This increased use requires a better understanding of the mechanical behavior of adhesives.

Two aspects will be distinguished, the cohesive aspect of bonding related to the behavior of the adhesive per se and the adhesive aspect related to the behavior of the interface between the adhesive and the substrate. In this paper, we will focus on the cohesive behavior of adhesives. The objective of this paper is to study and model the behavior of bonding in a mechanical system composed of beams glued in their axial direction. Different failure criteria are available to predict the joint strength, based on the stress analysis obtainable by a finite element analysis or a closed-form model. Finite element analysis are preferable for complex geometry and material properties. However, in many practical cases closed-form models for adhesive bonded composite joints are very desirable, because they provided fast and accurate analysis. The main problem in the numerical resolution of gluing problems comes from the thinness, and sometimes the low stiffness of the glue. A complete finite element calculation requires extremely refined mesh sizes. One of the ways to reduce the computational cost of the analysis was already proposed in particular by the authors in previous works and consists of using simplified models to simulate the behavior of the glue (Lebon andRizzoni, 2010, 2011;[START_REF] Rizzoni | Higher order model for soft and hard elastic interfaces[END_REF]. There are many simplified models in the literature. The perfect interface model is obviously the simplest. Spring-type models are also widely used [START_REF] Benveniste | Imperfect soft and stiff interfaces in two-dimensional elasticity[END_REF][START_REF] Geymonat | Mathematical analysis of a bonded joint with a soft thin adhesive[END_REF][START_REF] Goland | Stresses in cemented joints[END_REF][START_REF] Hashin | Thin interphase/imperfect interface in elasticity with application to coated fiber composites[END_REF][START_REF] Klarbring | Derivation of the adhesively bonded joints by the asymptotic expansion method[END_REF]. There is a wide literature exploiting beam theories in the analysis of adhesive joints [START_REF] Ascione | Curve adhesive joints[END_REF][START_REF] Bennati | An enhanced beam-theory model of the asymmetric double cantilever beam (ADCB) test for composite laminates[END_REF][START_REF] Cheng | Analysis of adhesive bonded joints with nonidentical adherends[END_REF][START_REF] Goland | Stresses in cemented joints[END_REF][START_REF] Jiang | Two-dimensional analytical stress distribution model for unbalanced FRP composite single-lap joints[END_REF][START_REF] Kanninen | A dynamic analysis of unstable crack propagation and arrest in the DCB test specimen[END_REF][START_REF] Kondo | Analysis of double cantilever beam specimen[END_REF][START_REF] Luo | Linear and higher order displacement theories for adhesively bonded lap joints[END_REF][START_REF] Olsson | A simplified improved beam analysis of the DCB specimen[END_REF][START_REF] Qiao | Mechanics and fracture of crack tip deformable bi-material interface[END_REF][START_REF] Su | Analytical model for adhesively bonded composite panel-flange joints based on the Timoshenko beam theory[END_REF][START_REF] Whitney | Stress analysis of the double cantilever beam specimen[END_REF][START_REF] Williams | End corrections for orthotropic DCB specimens[END_REF][START_REF] Wu | Stress-function variational method for stress analysis of bonded joints under mechanical and thermal loads[END_REF][START_REF] Wu | Stress-function variational method for interfacial stress analysis of adhesively bonded joints[END_REF][START_REF] Zou | An analytical solution for the analysis of symmetric composite adhesively bonded joints[END_REF], just to cite a few. Other studies focus on obtaining composite beam theories via a rigorous asymptotic analysis, see e.g. [START_REF] Serpilli | Asymptotic analysis of a multimaterial with a thin piezoelectric interphase[END_REF]Lenci, 2008, 2012;[START_REF] Serpilli | Asymptotic analysis of a multimaterial with a thin piezoelectric interphase[END_REF].

1

In this paper, a model of imperfect interface developed in the threedimensional setting by the authors in a previous paper is applied [START_REF] Rizzoni | Higher order model for soft and hard elastic interfaces[END_REF]. Practically, the glue (a volume) is replaced by an elastic interface (a surface) that keeps in memory the main mechanical characteristics of the glue. A classical theory for obtaining this material interface is the method of matched asymptotic expansions. The standard result is that, if the material of the adhesive is soft, at the order zero of the asymptotic expansion the adhesive turns out to modeled by an elastic surface, across which the stress components are continuous and are related via a linear relationship to the jumps of the displacement vector field. In the literature, these transmission conditions are said to provide a spring-type interface model. If the material of the adhesive is hard, then the order zero of the asymptotic expansion returns the classical model of perfect interface, where both the traction and displacement vector fields are continuous across the interface.

For the interface model considered in this paper, the asymptotic expansion is performed up to the first order (the second term in the asymptotic expansion). This introduces a jump on both the displacement vector and the stress vector fields across the interface, corresponding to a model of imperfect interface (Lebon andRizzoni, 2010, 2011;[START_REF] Rizzoni | Higher order model for soft and hard elastic interfaces[END_REF]. The related transmission conditions turn out to be more complex than those of the perfect interface model or of the spring-type interface model, because they incorporate higher order terms in ε.

However, they are more general, including the two classical contact models as special cases. In addition, they avoid the need of a priori specification of soft or hard material for the adhesive. The presence of higher order terms in ε introduces, via the appearance of first and second derivatives of the displacements, a non-local character of the higher order interface model. This is expected to provide a better approximation of the behavior of a thin adhesive interphase when its thickness ε is still small but it tends to become more and more comparable with the thicknesses of the adherent beams, h ± . A numerical evidence of this fact can be found in the example of a symmetric double cantilever beam given in Section 3.2. Fig. 11 shows that, as the ratio ε/ h ± increases, the tangential stress predicted by the higher order interface model also increases. On the contrary, for a symmetric double cantilever beam, the spring-type interface model always estimates the shear stress to vanish throughout the glue line, an indication that, for a more refined stress analysis, the higher order interface model considered in this paper provides better results.

Another original result of the paper is the introduction of higher order terms in the Neumann boundary conditions. These extra terms, emerging at the first order in [START_REF] Rizzoni | Higher order model for soft and hard elastic interfaces[END_REF], represent forces distributed at the adhesive-adherent interface boundary and arise from the presence of a boundary layer, cf. [START_REF] Klarbring | Derivation of the adhesively bonded joints by the asymptotic expansion method[END_REF]. In the present paper, these forces are make explicit and accounted for. For the example of the double cantilever beam given in Section 3.2, we numerically show that these higher order terms have to be taken into account when applying the higher order interface model. Indeed, when these additional terms are neglected, the axial forces in the adherent beams fail to satisfy Saint-Venant's principle, cf. Fig. 12.

The paper is organized as follows. Section 2 is devoted to the formulation of the problem. In particular, the main equations of the problem are presented, based on the Timoshenko beam model for the adherents and on the higher order interface model for the adhesive. For simplicity, the higher order interface model is specialized to the case of an isotropic adhesive, but other material symmetries could be investigated. The equations governing the equilibrium of the composite and incorporating the higher order interface model can be easily numerically implemented in commercial programming platforms like MATLAB® or Mathematica. Nevertheless, in Section 2.4 we provide a strategy to obtain an analytical solution, allowing to obtain closed-form solutions for the interfacial (peel and shear) stresses. Closed-form solutions are very important and worthwhile being developed, because they can be useful as benchmark solutions for numerical analysis.

In the third Section, two applications are analyzed. The first application is the shear of two adherent beams made of different linear elastic materials and joined by a thin isotropic adhesive layer. Using the higher order interface model, we develop closed-form solutions for the stresses and the displacement fields and we compare them with the exact threedimensional closed-form solution, available for this configuration. As a second application, a double cantilever beam, both in the balanced and unbalanced configurations, is considered. The interfacial stresses and displacements distributions are calculated numerically and validated with those obtained by a finite element analysis. For both applications, the role of the extra terms arising from the boundary layer is discussed. Finally, in Section 4, some conclusions and perspectives are given.

Formulation of the problem and strategy solution

Consider two elastic beams of different materials, unequal thickness, h + and h -, same lateral width b and same length L (Fig. 1). The beams are joined by means of a elastic adhesive layer with thickness ε much smaller than the other dimensions. A three dimensional coordinate system is introduced with origin at the middle height of the adhesive and axis 2 aligned with the direction of the axes of the two beams. The composite structure may be subjected to distributed loads along its length, and constrained at the two extremities or, in alternative, subjected to concentrated loads at the extremities. With reference to Fig. 1 for notation, the kinematics of the adherent beams is described by the

vectors u ± = {u ± 1 , u ± 2 , u ± 3 } T with u ± 1 = 0 u ± 2 = u ± - ( x 3 ∓ h ± 2 ) φ ± u ± 3 = v ± (1)
where u ± , v ± are the axial and transversal displacements of the beams centerlines, respectively, and φ ± are the cross-sectional rotations; u ± , v ± and φ ± are function of the axial coordinate x 2 ∈ [0, L]. In the above equations and throughout the derivations, the thickness of the adhesive layers is assumed to be much smaller than the thickness of the adherents, i.e. ε/h ± << 1, thus, the distance of the centerlines of the beam with respect to the reference system has been approximated to h ± /2. With this notation, the shear and peel deformation of the adhesive layer, [u] and [v], are defined as the jumps of the displacement along the x 2 and x 3 directions at the adherents-adhesive interfaces respectively,

[u] = u + + h + 2 φ + -u -+ h - 2 φ -, [v] = v + -v -. ( 2
)
The strain tensor takes the form where i 2 , i 3 are the unit vectors of the x 2 and x 3 axes, respectively. The non vanishing strain components are:

e ± = e ± (i 2 ⊗ i 2 ) + 1 2 γ ± (i 2 ⊗ i 3 + i 3 ⊗ i 2 ) (3)
e ± = ũ± ,2 -x 3 φ ± ,2 γ ± = v ± ,2 -φ ± (4)
where the subscript ", 2" is taken to denote the spatial derivative with respect to x 2 , and

ũ± = u ± ± h ± 2 φ ± (5)
is the displacement in the x 2 direction of the points of the two beams at the interface.

Governing equations for the adherents

The adherent beams are assumed to be subjected to external horizontal and vertical distributed loads q ± 2 and q ± 3 , respectively, acting on the top and on the bottom of the composite beams. The adherents are subjected to external actions concentrated at the end at x 2 = 0, denoted F ± 2,0 , F ± 3,0 and C ± 0 , and to external actions concentrated at the end at x 2 = L, denoted F ± 2,L , F ± 3,L and C ± L , see Fig. 1. In [START_REF] Klarbring | Derivation of the adhesively bonded joints by the asymptotic expansion method[END_REF][START_REF] Rizzoni | Higher order model for soft and hard elastic interfaces[END_REF] it is shown that at the order zero of the asymptotic expansion the adhesive is modeled by a spring-type interface, i.e. an elastic layer across which the stress components are continuous. When higher order terms of the asymptotic expansions are taken into account, the stress components are no longer continuous, cf. (Lebon andRizzoni, 2010, 2011;[START_REF] Rizzoni | Higher order model for soft and hard elastic interfaces[END_REF].

The corresponding situation is represented in Fig. 2, where τ + , σ + (τ -, σ -) are taken to denote the shear and the normal stresses acting on the bottom (top) of the adherent + ( -), respectively. Based on the free body diagrams shown on the right side of Fig. 2, where infinitesimal elements at x 2 with length dx 2 each are represented, the force and moment balances result in the following equilibrium equation:

N ± ,2 + q ± 2 ∓ τ ± b = 0, T ± ,2 + q ± 3 ∓ σ ± b = 0, M ± ,2 + T ± ∓ h ± 2 ( q ± 2 ± τ ± b ) = 0, (6)
where N ± , T ± and M ± are the axial force, the shear force and the bending moment in each beam, respectively, b is the specimen width. The axial force, the shear force and the bending moment are related to the kinematic parameters by the classical constitutive equations:

N ± = A ± u ± ,2 , T ± = C ± ( v ± ,2 -φ ±
) ,

M ± = D ± φ ± ,2 . (7)
Here, A ± , C ± and D ± are the extensional, shear and bending stiffnesses of the adherents, respectively (no bending-extension coupling is present here, as each beam is homogeneous and the kinematics is described from the centerline of the two beams) [START_REF] Jones | Mechanics of Composites Materials[END_REF].

Interface laws for thin elastic adhesives

In (Lebon andRizzoni, 2010, 2011;[START_REF] Rizzoni | Higher order model for soft and hard elastic interfaces[END_REF] the equilibrium problem of two bodies joined by a "soft" or a "hard" linear elastic thin adhesive interphase is studied; a soft adhesive is defined by the condition that its elasticity coefficients rescale like the adhesive thickness, ε, while for a hard adhesive the elasticity coefficients are independent of ε. In the limit of a vanishing thickness, ε→0, a spring-type interface law is classically obtained for a soft material, while perfect contact, i.e. the continuity of the displacement and traction vectors, is found for a hard material (cf. also [START_REF] Benveniste | Imperfect soft and stiff interfaces in two-dimensional elasticity[END_REF][START_REF] Geymonat | Mathematical analysis of a bonded joint with a soft thin adhesive[END_REF][START_REF] Goland | Stresses in cemented joints[END_REF][START_REF] Hashin | Thin interphase/imperfect interface in elasticity with application to coated fiber composites[END_REF][START_REF] Klarbring | Derivation of the adhesively bonded joints by the asymptotic expansion method[END_REF]). When the thickness of the adhesive, ε, is not so small, higher order terms in the asymptotic expansions with respect to ε have to be considered in the derivation of the interface law. In (Lebon andRizzoni, 2010, 2011;[START_REF] Rizzoni | Higher order model for soft and hard elastic interfaces[END_REF], different higher order interface models for soft and hard materials are calculated using two methods: the matched asymptotic expansion technique based on the strong formulation of classical continuum mechanics equations (compatibility, constitutive and equilibrium equations) and an energy method based on the minimization of the potential energy. In the end, it can be shown that it is possible to obtain a condensed form of transmission conditions summarizing both the orders zero and one of the two cases of soft and hard interface materials in only one couple of equations. The result is an implicit formulation of transmission conditions for the adhesive interface, prescribing both the jumps of the traction and displacement vector fields at the interface, cf. Eqns. (128,129) in [START_REF] Rizzoni | Higher order model for soft and hard elastic interfaces[END_REF]. In the present context, where all fields are independent of x 1 , the implicit formulation can be simplified as follows:

[σi 3 ] = -ε ( K 32 ( K 33 ) -1 〈σ ,2 i 3 〉 + ( K 22 -K 32 ( K 33 ) -1 K 23 ) 〈u ,22 〉 ) , ( 8
) [u] = ε ( K 33 ) -1 ( 〈σi 3 〉 -K 23 < u ,2 > ) , (9)
where

[f (x) ] := (f (x, 0 + ) -f (x, 0 -) ), < f (x) >:= 1 2 (f (x, 0 + ) + f (x, 0 -) ) (10)
are taken to denote the jump and the average of the quantity f(x) across the interface, respectively. Notably, the above interfacial relations ( 8), ( 9) do not present the additional terms arising from the matching conditions (Dumout et al., 2018, Eqns. 12-23]. In fact, it has been shown in [START_REF] Dumont | Soft and hard interface models for bonded elements[END_REF] that the matching terms should not be considered when dealing with thin beams, as in the present model.

Denoted with b ijkl , i, j, k, l = 1, 2, 3, the elasticity coefficients of the adhesive material, the matrices K jl in ( 8), ( 9) are defined as:

( K jl ) ki : = b ijkl . ( 11
)
In view of the symmetry properties of the elasticity tensor b, the matrices K lj have the property that K jl = (K lj ) T , with j, l = 1, 2, 3.

In [START_REF] Lebon | Asymptotic behavior of a hard thin linear interphase: an energy approach[END_REF][START_REF] Rizzoni | Higher order model for soft and hard elastic interfaces[END_REF], it has been shown that, for an isotropic and homogeneous adhesive material with Lamé coefficients λ, μ, the matrices K jl are given by

K jj = (2μ + λ) i j ⊗ i j + μ (i l ⊗ i l + i k ⊗ i k ), j ∕ = l ∕ = k (12) K jl = μ ( i j ⊗ i l ) + λ ( i l ⊗ i j ) , j ∕ = l, (13)
which yield in particular

K 22 = ⎛ ⎝ μ 0 0 0 2μ + λ 0 0 0 μ ⎞ ⎠ , K 33 = ⎛ ⎝ μ 0 0 0 μ 0 0 0 2μ + λ ⎞ ⎠ , K 23 = ⎛ ⎝ 0 0 0 0 0 μ 0 λ 0 ⎞ ⎠ , K 32 = ⎛ ⎝ 0 0 0 0 0 λ 0 μ 0 ⎞ ⎠ . ( 14
)
In the present context of adherent beams, where only the normal and tangential stresses are taken into account, cf. Fig. 2, jumps and averages of the traction and displacement fields in (8,9) reduce to

[σi 3 ] = ⎧ ⎨ ⎩ 0 [τ] [σ] ⎫ ⎬ ⎭ , < σi 3 >= ⎧ ⎨ ⎩ 0 < τ > 〈σ〉 ⎫ ⎬ ⎭ , [u] = ⎧ ⎨ ⎩ 0 [u] [v] ⎫ ⎬ ⎭ , < u >= ⎧ ⎨ ⎩ 0 < u > 〈v〉 ⎫ ⎬ ⎭ ,
where [τ], [σ], 〈τ〉, < σ > and 〈v〉 are defined according to (10), [u] and [v] are given by (2) and < u > is defined as the average of the axial displacement across the adhesive interface as follows:

〈u〉 : = 1 2 (u + + u -) + h + 4 φ + - h - 4 φ -. ( 16
)
Substituting ( 14) and ( 15) into ( 8),( 9), one obtains the following interfacial relations:

[τ] = -ε λ (λ + 2μ) < σ ,2 > -4εμ (λ + μ) (λ + 2μ) < u ,22 > (17) [σ] = -ε < τ ,2 > (18) [u] = ε μ < τ > -ε < v ,2 > (19) [v] = ε (λ + 2μ) < σ > -ε λ (λ + 2μ) 〈u ,2 〉, (20)
where we recall that the subscript ", 2" is taken to denote the spatial derivative with respect to x 2 . To illustrate the interfacial relations ( 17)-( 20), it is convenient to solve for τ ± and σ ± , in order to rewrite them into the following form:

τ ± = μ ε [u] + μ〈v ,2 〉 ∓ λ 2 [ v ,2 ] ∓ ε 2 (λ + 2μ) < u ,22 > , ( 21
)
σ ± = (λ + 2μ) ε [v] + λ < u ,2 > ∓ μ 2 [ u ,2 ] ∓ ε μ 2 〈v ,22 〉. ( 22
)
These relations imply that, due to zero and first order terms in ε, the interfacial stresses are not continuous, i.e. τ + ∕ = τ -and σ + ∕ = σ -. Note also that the terms in ε -1 are the standard ones appearing in the spring- type interfacial model. Indeed, if the adhesive material parameters of the adhesive, λ and μ, are taken to rescale like ε, as classically assumed for a soft adhesive, λ = λε and μ = με, then the interfacial relations (21), ( 22) reduce to the spring-type interface equations, up to neglecting first and second order terms in ε :

τ ± = μ[u], (23) σ ± = ( λ + 2μ)[v]. (24)
On the other hand, a hard adhesive is characterized by λ and μ independent of ε. Substituting into ( 21) and ( 22) and taking the limit ε→0, the relations of a perfect interface are recovered:

[u] = 0, (25)
[v] = 0.

The above discussion indicated that the interfacial relations ( 21), ( 22) include the classical contact models of spring-type interface and perfect interface as special cases. In fact, they are more general, in the sense that, avoiding the need of a priori specification of soft or hard material, they describe the behavior of an elastic surface incorporating higher order terms in ε. Indeed, in ( 21) and ( 22), the zero and first order terms in ε introduce, via the presence of first and second derivatives of the displacements, a non-local character of the interfacial relations. This is expected to introduce a better approximation of the behavior of a thin interphase when its thickness ε is still small but it tends to become more and more comparable with the thicknesses of the adherent beams, h ± . A numerical evidence of this fact will be given in the example of a symmetric double cantilever beam given in Section 3.2. In particular, Fig. 11 shows that, as the ratio ε/h ± increases, the tangential stress predicted by ( 21) and ( 22) increases. On the contrary, for a symmetric double cantilever beam, the spring-type interface always estimates the shear stress to vanish throughout the glue line, an indication that, for a more refined stress analysis, the higher order interface model provided by ( 21) and ( 22) gives better results.

Emerging forces at the adhesive-adherent interface boundary

The zero order solution, classically stated as spring-type or perfect contact transmission conditions, fails to describe the occurrence of a boundary layer, which is instead correctly predicted by an asymptotic analysis taking into account higher order terms, cf. [START_REF] Klarbring | Derivation of the adhesively bonded joints by the asymptotic expansion method[END_REF][START_REF] Rizzoni | Higher order model for soft and hard elastic interfaces[END_REF]. The analysis developed in [START_REF] Rizzoni | Higher order model for soft and hard elastic interfaces[END_REF] reveals the emergence of non-equilibrated stress resultants, scaling like ε, at the lateral boundary of the adhesive (cf. (Qiao and Wang, 2004, Eqn. (119)]):

F * : = ± ε ( K 32 ( K 33 ) -1 〈σi 3 〉 + ( K 22 -K 32 ( K 33 ) -1 K 23 ) < u ,2 > ) . ( 27
)
In ( 27), the plus sign applies at x 2 = L, and the minus sign at x 2 = 0.

Note that, when ε goes to zero, the forces ( 27) vanish, indicating that these terms are lost at zero order. Substituting ( 14) and ( 15) into ( 27), one obtains

F * = ± ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 0 εb λ (λ + 2μ) < σ > + 4εbμ (λ + μ) (λ + 2μ) < u ,2 > εb〈τ〉 ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ . ( 28
)
Using the relations ( 21) and ( 22), the forces F * can be rewritten as

F * = -F +,* -F -,* at x 2 = 0, (29)
F * = + F +,* + F -,* at x 2 = L, ( 30
)
with

F ±,* = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 0 ±εbλv ± + εb 4 (λ + 2μ) ( 2u ± ,2 ± h ± φ ± ,2
) bμ

( ± u ± + h ± 2 φ ± + ε 2 v ± ,2 ) ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ . ( 31
)
As shown in [START_REF] Klarbring | Derivation of the adhesively bonded joints by the asymptotic expansion method[END_REF], in order to take into account the presence of the forces (31) on the lateral boundary of the adhesive, complementary terms should be added to the classical asymptotic expansions of the stress and displacement fields. For the model proposed in this paper, we assume that the extra forces (31) arising from the boundary layer are applied at the two extremities of the adhesive layer. This leads to modify the Neumann boundary conditions at the boundary of the adherent beams in order to equilibrate the forces (31). In particular, the Neumann boundary conditions are modified as follows:

N ± (0) = -F ± 2,0 -F ±,* 2 (0), N ± (L) = +F ± 2,L -F ±,* 2 (L), T ± (0) = -F ± 3,0 -F ±,* 3 (0), T ± (L) = +F ± 3,L -F ±,* 3 (L), M ± (0) = -C ± 0 -C ±,* (0), M ± (L) = +C ± L -C ±,* (L), (32)
where the additional terms

F ±,* 2 (0) = ± εbλv ± (0) + εb 4 (λ + 2μ) ( 2u ± ,2 (0) ± h ± φ ± ,2 (0) ) , ( 33
)
F ±,* 3 (0) = bμ ( ± u ± (0) + h ± 2 φ ± (0) + ε 2 v ± ,2 (0) ) , ( 34
) C ±,* (0) = ± F ±,* 2 (0) ( h ± 2 + ε 2 ) (35) F ±,* 2 (L) = ± εbλv ± (L) + εb 4 (λ + 2μ) ( 2u ± ,2 (L) ± h ± φ ± ,2 (L) ) , ( 36
)
F ±,* 3 (L) = bμ ( ± u ± (L) + h ± 2 φ ± (L) + ε 2 v ± ,2 (L) ) , ( 37
) C ±,* (L) = ± F ±,* 2 (L) ( h ± 2 + ε 2 ) , ( 38
)
are directly given by (31) for the components along the x 2 and x 3 directions and by their moments about the beams axes for the torques C ±,* (0) and C ±,* (l).

Dirichlet boundary conditions are unaffected by the presence of the boundary layer. In other words, the displacement boundary conditions on axial and cross-sectional displacements, and rotations are the classical ones for the beam theory, and we do not expect them to include any terms arising from (27). The extra forces (27) arising from the boundary layer have to be taken into account when considering the higher order interfacial relations ( 21) and ( 22). This will be made clear in the example of the double cantilever beam given in Section 3.2. As shown by Fig. 12, if the additional terms ( 33)-( 38) are neglected, then the axial forces N ± fail to satisfy Saint-Venant's principle. As a final remark, we note that in a three-dimensional setting conditions (29) and ( 30) have been found to be useful to write a weak formulation taking into account the imperfect interface equations ( 8) and ( 9), cf. [START_REF] Dumont | An asymptotic approach to the adhesion of thin stiff films[END_REF][START_REF] Dumont | Soft and hard interface models for bonded elements[END_REF][START_REF] Serpilli | An asymptotic derivation of a general imperfect interface law for linear multiphysics composites[END_REF].

Solution strategy

The equilibrium equation ( 6) together with the constitutive equation ( 7), the transmission conditions ( 21) and ( 22), the Neumann boundary conditions ( 32) and Dirichlet boundary conditions can be easily numerically implemented in commercial programming platforms like MATLAB® or Mathematica. Nevertheless, in this Section we provide a strategy to obtain a closed-form solutions for the interfacial (peel and shear) stresses. Closed-form solutions are very important and worthwhile being developed, because they are useful as benchmark solutions for numerical analysis. The first step of the strategy is to reduce the number of unknowns and get a smaller set of equations. In particular, following [START_REF] Bennati | An enhanced beam-theory model of the asymmetric double cantilever beam (ADCB) test for composite laminates[END_REF], we aim at reducing to a set of equations whose unknowns are the four interfacial stresses σ ± and τ ± . To do so, we substitute equation ( 7) into ( 6) to obtain u ± ,22 , φ ± ,222 and v ± ,2222 in terms of the interfacial stresses and their derivatives. The resulting equations are then substituted into ( 17)-( 20). After some transformations, omitted here for the sake of brevity, the following system of four coupled 4th-order ODEs is obtained:

η + σ + ,22 + η -σ - ,22 + ζ + σ + + ζ -σ -+ m = 0, (39) θ
( σ - ,2222 + σ + ,2222
)

+ κ + σ + ,22 + κ -σ - ,22 + ρ + σ + + ρ -σ -+ n = 0, (40) τ
- ,2 = α -σ -+ β -σ + -γ ( σ - ,22 + σ + ,22 ) -f , (41) τ + ,2 = β + σ -+ α + σ + + γ ( σ - ,22 + σ + ,22 ) + f . ( 42
)
In the above equations, the coefficients f, m and n are related to the distributed loads, as reported in the Appendix.

The other coefficients, α ± , β ± , γ, ζ ± , η ± , θ, κ ± and ρ ± , take into ac- count the geometry and stiffness of the beams (h ± , A ± , B ± , C ± ), and the geometry and elasticity of the adhesive layer (ε, λ, μ), see the Appendix.

We look for a solution of equations ( 39) and (40) in the form of the sum of the solution to the associated homogeneous system and a particular solution σ ± * . In case of spatially homogeneous material constants, the solution to the associated homogeneous system is written as

σ ± (x 2 ) = ∑ 6 i=1 F ± i exp(λ i x 2 ), (43)
where λ 1 , λ 2 , …λ 6 are constants. Substituting (43) into the homogeneous system associated to (39), ( 40), an eigenvalue problem is obtained and, consequently, the eigenvalues λ i , i = 1, …, 6, must satisfy the characteristic equation

θ(η --η + )λ 6 +(θ(ζ --ζ + )+κ + η --κ -η + )λ 4 +(κ + ζ --κ -ζ + +ρ + η --ρ -η + )λ 2 +(ρ + ζ --ρ -ζ + )=0, (44)
and the constants F ± 1 , F ± 2 , …F ± 6 are found to satisfying the six following (eigenvectors) condition:

( η + λ 2 i + ζ + ) F + i + ( η -λ 2 i + ζ -) F - i = 0, i = 1, …6. (45)
By substituting equation ( 43) into ( 41),( 42) and integrating, we obtain the solution for the tangential stress,

τ ± (x 2 ) = ∑ 6 i=1 G ± i exp(λ i x 2 ) + G ± 7 + τ ± * (x 2 ) (46)
where τ ± * is a particular solution incorporating σ ± * and the primitive of f, G ± i , i = 1, …7, are integration constants. In particular, the constants G ± i , i = 1, …6, satisfy the conditions:

λ i G + i - ( β + F - i + α -F + i ) -γλ 2 i ( F - i + F + i ) = 0, i = 1, …6, (47)
λ i G - i - ( α -F - i + β -F + i ) + γλ 2 i ( F - i + F + i ) = 0, i = 1, …6. (48)
Because λ i , i = 1, …, 6, can be calculated using the characteristic equation (2.4) and the integration constants

F - i , G ± i , i = 1, …, 6, are known if the constants F + i , i = 1, …, 6
, the stresses σ ± and τ ± are completely determined up to the calculation of the integration constants F + i , i = 1, …, 6. To calculate the latter ones, it is necessary to further integrate the equilibrium equation ( 6), and to use the constitutive equation ( 7), the interfacial equations ( 17)-( 20) and the boundary conditions (32) in a not straightforward way. The procedure is described in the following.

First, the internal forces are calculated. They can be obtained by substituting the expressions for the interfacial stresses ( 43) and ( 46) into the equilibrium equations of the two beams (6) and then by integrating them with respect to x 2 . This integration gives the analytical expressions for the internal forces, where six new integration constants, G ± 8 , G ± 9 , G ± 10 appear. Next, substituting the expressions for the internal forces into the constitutive equation ( 7) and integrating with respect x 2 , the analytical expressions for the displacements can be deduced. This process generates six more integration constants, G ± 11 , G ± 12 , G ± 13 . To summarize, up to now there are 38 integration constants to be determined (F ± i , i = 1, …6, and G ± j , j = 1, …13), but the 18 conditions (45), ( 47), ( 48) together with the twelve the boundary conditions (the natural boundary conditions (32) or analogous essential boundary conditions prescribing the displacements or mixed boundary conditions) appear insufficient because they consist of only 30 equations. In fact, by introducing the expressions for the interfacial stresses and displacements, obtained by integration, into the interfacial relations ( 17)-( 20), additional relations between the constants can be found, only 8 of which are linearly independent. These 8 relations, together with the 30 conditions cited above, allow to finally calculate the 38 integration constants and thus to completely solve the problem.

In the next Section, we illustrate the solution strategy by means of the simple example of the shear of a composite block. For this example, the analytical solution calculated with the proposed composite beam model is compared with the exact analytical solution obtained in the framework of linear elasticity.

Model validation and numerical examples

In this Section, two examples are presented to validate the model proposed in this paper.

The first example is the shear of a composite block, for which the closed form solution of the full three-dimensional problem is available and thus directly comparable with the approximated solution obtained by using the composite beam model. In [START_REF] Lebon | Asymptotic analysis of a thin interface: the case involving similar rigidity[END_REF][START_REF] Rizzoni | Higher order model for soft and hard elastic interfaces[END_REF] the closed form solution of the full three-dimensional equilibrium problem of a composite block under shear had already been compared with the approximated solution obtained by solving the three-dimensional equilibrium problem of two adherents in contact through soft and hard interface models for the adhesive. In Subsection 3.1, an original solution for a composite beam under shear is presented. The solution is given by simple shears of the two beams superimposed to a relative horizontal translation mimicking the shear of the adhesive layer. In the adherents, this solution correctly reproduces the corresponding analytical solution of the three-dimensional equilibrium problem of the composite block.

As a second example, the model is applied to the numerical analysis of a double cantilever beam, both in balanced and unbalanced configurations. As verification, interface stresses and displacements are obtained by three different methods: the present composite beam model based on the higher order interface equations ( 21) and ( 22); a composite beam model based on the classical spring-type interface equations, ( 23) and ( 24); a finite element analysis (FEA).

Shear of a composite block

The geometric and load configuration of a composite beam under shear is illustrated in Fig. 3. The upper and lower beams are loaded by uniform load of intensity q distributed along the horizontal direction. To prevent rigid body motions, the lower beam is constrained by two supports. Concentrated vertical forces of intensity qh + are applied at the two ends of the upper beams. In the boundary conditions (32), the additional terms arising from the boundary layer are initially neglected. Their effect on the equilibrium solution of the composite will be discussed at the end of this Section.

The load configuration depicted in Fig. 3 models the load distribution of a composite block under shear, whose exact analytical solution is available in the framework of linear elasticity and it is described in [START_REF] Lebon | Asymptotic analysis of a thin interface: the case involving similar rigidity[END_REF][START_REF] Rizzoni | Higher order model for soft and hard elastic interfaces[END_REF]. For the configuration of Fig. 3, we have q ± 2 = ±q and q ± 3 = 0. Thus, f, m, n = 0 and system (39)-( 42) becomes homogenous. Given the simple nature of the shear problem, we rely on the physical intuition in prescribing the natural choices τ ± = q b , σ ± = 0, (49) where forces and torques arising from the boundary layer have been neglected.

N ± = 0, (50)
Most of the boundary conditions are identically satisfied, the others give

G - 11 = 0, (64)
G - 12 = qh - C -, ( 65
)
G - 13 = 0. ( 66
)
The remaining three integration constants, G + 11 , G + 12 and G + 13 , are determined through three additional conditions, which can be obtained by substituting the expressions ( 56)-( 58) together with (49) into the interfacial relations ( 17)-( 20). We obtain the two equations

G + 11 + h + 2 G + 12 - q(h -) 2 2C -= qε μb -ε ( qh + C + + G + 12 ) , ( 67
)
( qh + C + + G + 12 ) x 2 + G + 13 = 0, (68)
which have to be satisfied for each x 2 in (0, L), giving the conditions Table 1 Geometrical parameters of the symmetric (S-FIT) and asymmetric (A-FIT) double cantilever beam configurations studied in the numerical simulations.

G + 11 = qε μb + q(h + ) 2 2C + + q(h -) 2 2C -, ( 69
G + 12 = - qh + C + , ( 70
)
G + 13 = 0. ( 71
)
Substituting the latter expressions into ( 56)-( 58) together with ( 64)-( 66) gives

u + = qε μb + q(h + ) 2 2C + + q(h -) 2 2C -, ( 72
)
u -= 0, ( 73
)
φ ± = - qh ± C ± , ( 74
) v ± = 0. ( 75
)
This solution corresponds to a vanishing deformation of the lower beam, and to a shear deformations of amount q(h -) 2 2C -+ q(h -) 2 2C -+ qε μb in the upper beam, the term qε μb being the contribution of the adhesive layer. The expressions ( 72)-( 75), substituted back into (1), yield the same solution obtained for the shear of a composite block in the more general context of three-dimensional linear elasticity, up to rigid motions and up the choice C ± = μ ± bh ± , with μ ± the shear moduli of the materials of the adherents (Lebon and Rizzoni, 2010, Eqn. (48)]. Finally, substitution of ( 72)-( 75) into the expressions ( 33)-( 38) of the forces arising from the boundary layer gives

F * 2 (0) = F * 2 (L) = 0, ( 76
)
F * 3 (0) = F * 3 (L) = qε. ( 77
)
The forces (77) correspond to the force resultants arising from the tangential stress acting on the lateral surface of the adhesive. Because the approximated solution ( 72)-( 75), calculated without taking into account the boundary layer, correctly reproduces the exact threedimensional solution in the adherents for the shear of the three-layer composite block, we are led to conclude that the boundary layer has only a local effect for this loading case.

Double cantilever beam

The double cantilever beam specimen is widely used to test the mixed mode fracture toughness of adhesive joints. Recently, [START_REF] Bui | The Flexible Initiation Test (FIT): a new experimental test to characterize fracture initiation in mode I at the free edge of bonded assemblies[END_REF] have proposed an original test, the Flexible Initiation Test (FIT), designed to characterize fracture initiation in mode I and decrease scattering in fracture initiation load measurements. Fig. 4 shows the DCB specimen in the symmetric FIT configuration and loading. In the configuration, deformable support beams are bonded to the composite body composed of adherents and adhesive. The presence of the (long) support beams has been experimentally proven to reduce the scattering of the fracture initiation load.

In this Section, six different configurations are analyzed: three symmetric configurations (S-FIT) and three asymmetric configurations (A-FIT), with geometrical and material parameters listed in Tables 1 and2. The geometrical parameters, shown in Fig. 4 and indicated in Table 1, have been chosen according to the numerical study and the experimental validation performed in [START_REF] Bui | The Flexible Initiation Test (FIT): a new experimental test to characterize fracture initiation in mode I at the free edge of bonded assemblies[END_REF]. The support beams and the adherents are chosen to be composed of the same material, aluminum.

Three different materials have been considered for the adhesive: • Scotch-Weld™ 3M 2216 B/A, a ductile epoxy adhesive, denoted as "adhesive 1", characterized by a low elastic modulus, 0.342 GPa, about two orders of magnitude lower than the elastic modulus of the adherents (70 GPa for aluminum);

• Araldite ® AV138M-1/Hardener HV998, denoted as "adhesive 2", a brittle epoxy adhesive with an elastic modulus of 4.13 GPa, about one order of magnitude lower than the elastic modulus of the adherents;

• Sikadur ®-30, denoted as "adhesive 3", an adhesive for bonding reinforcement with an elastic modulus of 12.8 GPa, comparable to the elastic modulus of the adherents.

In the simulations presented below, the thickness of the adhesive is fixed and equal to 0.1 mm.

The efficiency of the composite beam model based on higher order interface equations proposed in the present paper is tested against FEA, as the adhesive becomes stiffer (i.e. going from adhesive 1 to adhesive 3). For comparison, a composite beam model based on the classical spring-type interface equations has also been considered.

The equations of the composite beam model, ( 6), (7) coupled with the higher order transmission conditions (21), ( 22) have been implemented in Matlab™ (Release2018a), together with the following boundary conditions appropriate for the double cantilever beam model:

u ± (0) = 0, N ± (L) = F ±,* 2 (L), v ± (0) = 0, T ± (L) = ±P + F ±,* 3 (L), φ ± (0) = 0, M ± (L) = ±Pa + C ±,* (L), (78)
where F ±,* 2 (L), F ±,* 3 (L), and C ±,* (L) are the loading terms arising from the boundary layer (cf. Subsection 2.3). In addition to the beam model with the higher order transmission condition, equations ( 6) and ( 7), complemented by the boundary conditions (78), have been implemented in Matlab™ with the classical spring-type interface equations 23 and 24. In both cases, higher order or spring-type interface equations, the load applied to the cantilever arms has been taken equal to P = 10 MPa.

Note also that, in the plane stress case, the Lamé constants λ and μ utilized in the interfacial equations 23 and 24 or (21), ( 22) are related to the Young modulus and the Poisson ratio as follows:

λ = Eν (1 -ν) 2 , μ = E 2(1 + ν) . ( 79
)
The stress and displacement data calculated with the solutions of the composite beam model have been compared with the results obtained by a finite element analysis, performed using the commercial finite element software COMSOL™ Multiphysics 3.5. The Plane Stress -Structural Mechanics Module and triangular mesh elements have been utilized, with a minimum number of 8 elements through the adhesive thickness. present the interfacial stresses and displacements in adhesive 1 for the symmetric (S-FIT) and asymmetric (A-FIT) configurations, respectively. Figs. 5 and6 show that all the three methods agree very well for the whole range of overlap. The curves for the peel stress and the axial and transversal displacements calculated with classical spring-type interface model and the imperfect interface model display very small differences, except for the shear stress distribution for the symmetric configurations. Indeed, the shear stress turns out to vanish in the model with the springtype interface equations, but it is not vanishing and correctly predicted by the beam model with higher order interface equations, showing the same trend as FEA (cf. top right plots in Figs. 5 and6). In Figs. 789, 10 similar trends can be observed for the symmetric and asymmetric configures with adhesives 2 and 3, respectively. By comparison of the stresses distributions depicted in Figs. 5-10, it can be concluded that both peel and shear stresses increase slightly with increasing elastic modulus of the adhesive layer, for the parameters utilized in this study. An inverse trend can be observed for the displacements distributions, as expected. In other words, more compliant adhesive layers undergo larger deformations to accommodate the different deformations of the two adherents, therefore decreasing the interfacial stresses.

Fig. 11 illustrates the dependence of the normal and tangential interfacial stresses on the adhesive thickness ε for the case of adhesive 2 in the S-FIT configuration. The interfacial stress components are calculated using the spring-type model and the proposed higher order contact model for three different values of the adhesive thickness ε : 0.1 mm, 0.05 mm and 0.01 mm. The results shown in Fig. 11 indicate that for the normal stress there is almost no difference between the spring-type and the higher order interface model, at least for the set of data considered in the example. However, significant differences can be appreciated for the tangential stress. Contrary to the spring-type model which always predicts a vanishing shear stress along the glue line, the higher order contact model takes into account a non-vanishing distribution of the shear stress, increasing as the adhesive thickness increases. Note that, because the adherents' thicknesses h ± are kept fixed (4.6 mm), the numerical results actually illustrate the dependence of the interfacial stresses on the ratio ε/h ± . Then, one concludes that as the ratio ε/ h ± increases the spring-type interface model is no longer appropriate for a complete stress analysis, and more refined contact model like the higher order interface model presented in the manuscript should be implemented. Notably, all distributions satisfy Saint-Venant's principle, i.e. N ± , T ± and M ± become very small at sufficiently large distances from the loaded end (the right end for the double cantilever beam, cf. Fig. 4), except for the distribution of the axial forces N ± calculated using the higher order interface model implemented without taking into account the terms (33)-( 38). The fact that these N ± distributions do not satisfy Saint-Venant's principle is an indication that, when considering higher order interface model, it is mandatory to include the contribution of the boundary layer. One may wonder why the shear stress and bending moment distributions calculated with the higher order interface model without taking into account the boundary layer satisfy Saint-Venant's principle. In fact, this is due to the symmetry of the configuration, leading to the vanishing of the extra terms (33), ( 35), ( 36) and (38). In the case of an asymmetric A-FIT configuration, we have numerically verified (not reported here) that the shear forces T ± and the bending moments M ± calculated with the higher order interface model without taking into account the boundary layer terms do not satisfy Saint-Venant's principle.

Like most analytic models available in the literature, the beam model with imperfect interface is unable to predict the vanishing of the shear stress and the singularity of the peel stress at the free-edge. However, the good fitting to the interfacial stresses distributions predicted by FEA, especially for the shear stress, confirms the validity of the present model. Of course, some differences could be noted between the results obtained by the proposed model and the ones recovered by FEA. Indeed, it could be remarked that these differences could be expected a priori as they concern two different structural models, i.e. the two beam model and a 2D continuum solid model.

Conclusions

In this paper, we have propose an original composite beam model based on the Timoshenko beam model for the adherents and on the imperfect interface model proposed in [START_REF] Rizzoni | Higher order model for soft and hard elastic interfaces[END_REF] for the adhesive. We have shown that the method can successfully be used to calculate the interfacial and stresses (peel and shear stresses) distributions for the stress analysis of a double cantilever beam specimen.

The main advantages of the proposed composite beam model based on imperfect interface equations are as follows: (1) the imperfect interface equations are rigorously derived from a 3D continuum model, thus the interface stiffness constants have a clear physical meaning directly related to the material (type of material symmetry, elastic constants) and geometrical (thickness) parameters of the adhesive; (2) the model predicts correctly the shear stress distributions along the adherents/adhesive interfaces; (3) explicit closed forms are available for the peel and the shear stresses.

By modifying the boundary conditions of the set of governing equations, other joint configurations could be studied, including adhesively bonded single-lap joints, composite joints, etc. The higher order interface model developed in [START_REF] Rizzoni | Higher order model for soft and hard elastic interfaces[END_REF] and applied in the present paper is general enough to allow material symmetries for the adhesive different from the isotropic behavior, here considered for simplicity. Multiphysics material behavior could also be incorporated, based on some recent results [START_REF] Serpilli | An asymptotic derivation of a general imperfect interface law for linear multiphysics composites[END_REF].

The present work has been validated by FEA. The results obtained indicate that the present study provides an efficient method for accurate stress analyses in adhesive bonded joints and therefore can be used effectively for practical applications in engineering. The accurate evaluation of the interface stress would be particularly important when interface damage and delamination are considered. In these cases significant different structural responses could be obtained considering or neglecting higher order terms in the interface model. This aspect will be object of future researches.

Fig. 1 .

 1 Fig. 1. Geometry of the assemblage with applied loads.

Fig. 2 .

 2 Fig. 2. Free body diagrams of the two adherent beams (left), and of two infinitesimal elements of the beams (right).

Fig. 3 .

 3 Fig. 3. Composite beam under shear load.

)Fig. 4 .

 4 Fig. 4. a) double cantilever beam specimen in the symmetric FIT configuration (S-FIT), and b) in the asymmetric FIT configuration (A-FIT).

Fig. 5 .

 5 Fig. 5. Top figure: normal (left) and tangential (right) stresses at the upper and lower interfaces for a composite with adhesive 1 in the symmetric S-FIT configuration. Bottom figure: axial (left) and transversal (right) displacements at the upper and lower interfaces for a composite with adhesive 1 in the symmetric S-FIT configuration.

Fig. 6 .

 6 Fig. 6. Top figure: normal (left) and tangential (right) stresses at the upper and lower interfaces for a composite with adhesive 1 in the asymmetric A-FIT configuration. Bottom figure: axial (left) and transversal (right) displacements at the upper and lower interfaces for a composite with adhesive 1 in the asymmetric A-FIT configuration.

Fig. 7 .

 7 Fig. 7. Top figure: normal (left) and tangential (right) stresses at the upper and lower interfaces for a composite with adhesive 2 in the symmetric S-FIT configuration. Bottom figure: axial (left) and transversal (right) displacements at the upper and lower interfaces for a composite with adhesive 2 in the symmetric S-FIT configuration.

Fig. 8 .

 8 Fig. 8. Top figure: normal (left) and tangential (right) stresses at the upper and lower interfaces for a composite with adhesive 2 in the asymmetric A-FIT configuration. Bottom figure: axial (left) and transversal (right) displacements at the upper and lower interfaces for a composite with adhesive 2 in the asymmetric A-FIT configuration.

Fig. 9 .

 9 Fig. 9. Top figure: normal (left) and tangential (right) stresses at the upper and lower interfaces for a composite with adhesive 3 in the symmetric S-FIT configuration. Bottom figure: axial (left) and transversal (right) displacements at the upper and lower interfaces for a composite with adhesive 3 in the symmetric S-FIT configuration.

Fig. 10 .

 10 Fig. 10. Top figure: normal (left) and tangential (right) stresses at the upper and lower interfaces for a composite with adhesive 3 in the asymmetric A-FIT configuration. Bottom figure: axial (left) and transversal (right) displacements at the upper and lower interfaces for a composite with adhesive 3 in the asymmetric A-FIT configuration.

Fig. 11 .

 11 Fig. 11. Effect of the adhesive thickness on the normal (left) and tangential (right) stresses at the upper and lower interfaces for a composite with adhesive 2 in the symmetric S-FIT configuration.

Fig. 12 .

 12 Fig. 12. Distributions of axial forces (left), shear forces (middle) and bending moments (left) for a composite with adhesive 2 in the symmetric S-FIT configuration. The different curves refer to three different models: the spring-type interface model; the higher order interface model implemented without taking into account the boundary layer; the higher order interface model in the presence of the boundary layer.

  Fig. 12 illustrates the effect of the extra forces (33)-(38) introduced in Section 2.3 and arising from the boundary layer. The Figure shows the distributions of axial forces, shear forces and bending moments for a composite with adhesive 2 in the symmetric S-FIT configuration. The different curves have been calculated by using to three different interface models: the spring-type interface model; the higher order interface model implemented without taking into account the extra forces (33)-(38); the higher order interface model with the extra forces (33)-(38).

  with δ = 2AD(λ + 2μ) + bεμ(λ + μ) ( 4D + Ah 2 ) Lateral width of composite u ± : Axial displacements of adherent beams v ± : Transversal displacements of adherent beams φ ± : Cross-sectional rotations of adherent beams [f]: Jump of function f across the adhesive interface 〈f〉: Average of function f over the adhesive interface e ± : Strain tensors e ± , γ ± : Strain componentsũ± : Axial displacements of adherent beams at the adhesive interface q ± 2 : External horizontal distributed loads q ± 3 : External vertical distributed loads F ± 2,0 ,F ± 3,0 ,C ± 0 ,: External actions (forces and torque) acting on the adherent beams at the endx2 = 0 F ± 2,L ,F ± 3,L ,C ± L ,: External actions (forces and torque) acting on the adherent beams at the end x2 = L τ ± : Tangential interfacial stresses σ ± : Normal interfacial stresses N ± : Axial forces T ± : Shear forces M ± : Bending moments A ± : Extensional stiffnesses of adherent beams C ± : Shear stiffnesses of adherent beams D ± : Bending stiffnesses of adherent beams σ: Cauchy stress i3: Unit vector in the direction of x3 axis K jl : Reduced elasticity matrices of the adhesive material λ, μ: Lamé parameters of the adhesive material λ, μ: Rescaled Lamé parameters of a soft adhesive material F * : Force resultant acting on the lateral boundary of the adhesive and arising from boundary layer F ±,* : Boundary layer force resultants applied at the ends of the adherent beams F ±,* 2 , F ±,* 3 , C ±,* : Boundary layer actions (forces and torque) applied at the ends of the adherent beams f, m, n: Coefficients of Eqns. (39)-(42) related to external forcesα ± , β ± , γ, ζ ± , η ± , θ, κ ± , ρ ± : Coefficients of Eqns. (39)-(42) related to geometry and material propertiesF ± i , i = 1, 2, …6: Integration constants G ± i , i = 1,2, …13: Integration constants q: Distributed shear load a: Length of the support beams E: Elastic modulus of the adherents ν: Poisson's coefficient of the adherents

Table 2

 2 Adhesive and adherents elastic properties assumed in the numerical study.

	Quantity	Symbol	S-FIT	A-FIT	Unit
	Length of the support beams Width of the beams Length of adhesive Height of upper beam Height of lower beam Thickness of the adhesive	a b L h + h -ε	200.0 12.0 50.0 4.6 4.6 0.1	200.0 12.0 50.0 4.6 5.6 0.1	mm mm mm mm mm mm

Acknowledgements

RR was supported by the Italian Ministry of Education, University and Research (MIUR) under the Italian National Grant for Fundamental Research (FFABR 2017) and by the University of Ferrara, Italy, under the Research Programs FAR 2018 and FAR 2019. The research has been developed also thanks to the financial support of project PRIN2015 "Multi-scale mechanical models for the design and optimization of microstructured smart materials and metamaterials".

T ± = qh ± (51)

which identically verify system (39)-( 42) and the equilibrium equation ( 6). In terms of the solution strategy described in the previous Section, one can identify the following integration constants

Substituting the relations (49)-( 52) into the constitutive equation ( 7) and integrating, the following expressions for the displacement fields can be obtained:

where G ± 11 , G ± 12 and G ± 13 are six more integration constants. To calculate the integration constants, we take into account the boundary conditions. The presence of two supports in the lower adherent imposes the vanishing of the following displacement components:

The remaining (natural) boundary conditions are

Declaration of competing 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.

Appendix

The coefficients of system (39)-( 42) have the following form:

(82)

))

with

For symmetrical configurations, one has

and the coefficients ( 80)-( 87) simplify as follows:

] ,


Serge Dumont

Olivier Goubet
email: olivier.goubet@u-picardie.fr

Imen Manoubi

Decay of solutions to a water wave model with a nonlocal viscous term

Keywords: Nonlocal water wave equation, Decay rate of solutions Mathematics Subject Classification 35B40, 35Q35

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

In this article we focus on a single equation for a one-way wave equation that was investigated in [START_REF] Manoubi | Theoretical and numerical analysis of the decay rate of solutions to a water wave model with a nonlocal viscous dispersive term with Riemann-Liouville half derivative[END_REF], where fractional term is described by the Riemann-Liouville half derivative. This equation reads

u t + u x + 1 √ π ∂ ∂t t 0 u(s) √ t -s ds + uu x = u xx . (1.1)
Here, for the sake of convenience, the physical constants are set equal to 1. The nonlocal term provides both diffusion and dispersion to the equation [START_REF] Chen | Decay of solutions to a water wave model with a nonlocal viscous dispersive term[END_REF][START_REF] Manoubi | Theoretical and numerical analysis of the decay rate of solutions to a water wave model with a nonlocal viscous dispersive term with Riemann-Liouville half derivative[END_REF]. On the other hand, diffusion equations with non local operators are of great interest for mathematicians and also for other scientific communities (see [START_REF] Dutykh | Viscous-potential free-surface flows and long wave modelling[END_REF][START_REF] Haddar | Well-posedness of non-linear conservative systems when coupled with diffusive systems[END_REF][START_REF] Helie | Diffusive reprentations for the analysis and simulation of flared acoustic pipes with visco-thermal losse[END_REF][START_REF] Hristov | Fractional derivative with non-singular kernels, from the Caputo-Fabrizio definition and beyond: appraising analysis with emphasis on diffusion models, Chapter 10[END_REF][START_REF] Matignon | Asymptotic stability of linear conservative systems when coupled with diffusive systems[END_REF][START_REF] Montseny | Diffusion monodimensionnelle et intégration d'ordre 1/2[END_REF][START_REF] Montseny | Représentation Diffusive[END_REF]a n dt h e references therein).

The initial value problem for (1.1) was previously addressed in [START_REF] Goubet | Theoretical analysis of a water wave model with a nonlocal viscous dispersive term using the diffusive approach Advances in Nonlinear Analysis[END_REF][START_REF] Manoubi | Theoretical and numerical analysis of the decay rate of solutions to a water wave model with a nonlocal viscous dispersive term with Riemann-Liouville half derivative[END_REF]; we know that for initial data u 0 in L 2 (R), there exists a unique solution u ∈ C(R + ; L 2

x (R)) ∩ C 1/2 (R + ; H -2 x (R)). Besides, for more regular initial data u 0 ∈ S(R) in the Schwartz class, the solution belongs to C 1/2 (R + ; S(R)). The solution map t → u(t) is C 1 but at t = 0 + . In this article we are interested in the decay rate of solutions towards 0. We prove that any solution converges towards 0 with the expected decay rate; this is an improvement over [START_REF] Manoubi | Theoretical and numerical analysis of the decay rate of solutions to a water wave model with a nonlocal viscous dispersive term with Riemann-Liouville half derivative[END_REF] where a smallness assumption on the initial data was required. On the other hand, the numerical investigation of this decay rate has also received interest also in the last decade. Many authors have provided different numerical methods to approximate the nonlocal fractional operator [START_REF] Chen | Decay of solutions to a water wave model with a nonlocal viscous dispersive term[END_REF][START_REF] Galucio | The G α -scheme for approximation of fractional derivatives: application to the dynamics of dissipative systems[END_REF][START_REF] Galucio | An adaptation of the Gear scheme for fractional derivatives[END_REF][START_REF] Koh | Application of fractional derivatives to seismic analysis of base-isolated models[END_REF][START_REF] Manoubi | Theoretical and numerical analysis of the decay rate of solutions to a water wave model with a nonlocal viscous dispersive term with Riemann-Liouville half derivative[END_REF]. Some have proposed the so-called diffusive realization of fractional operators. This alternative provides a local representation of the fractional operator as a system of linear differential equation. Several applications of this approach can be found in [START_REF] Audounet | A threshold phenomenon in the propagation of a pointsource initiated flame[END_REF][START_REF] Helie | Diffusive reprentations for the analysis and simulation of flared acoustic pipes with visco-thermal losse[END_REF][START_REF] Matignon | Asymptotic stability of linear conservative systems when coupled with diffusive systems[END_REF][START_REF] Montseny | Diffusive representation for pseudodifferentially damped nonlinear systems[END_REF]. We use here this diffusive approach to prove theoretical results. Due to uniqueness, we will use in the sequel that a solution of the original equation with u 0 in L 2 (R) is solution to the diffusive representation and vice-versa. This article is organized as follows. We introduce a suitable diffusive representation for the equation, we state the main theoretical result and we proceed to the proof. We complete this article by numerical computations that illustrate our theoretical results.

Statement of the result

The aim of this article is to answer a question left open in [START_REF] Manoubi | Theoretical and numerical analysis of the decay rate of solutions to a water wave model with a nonlocal viscous dispersive term with Riemann-Liouville half derivative[END_REF]. Concerning the decay rate of solutions, the third author proved that if u 0 L 1 (R) is small enough, then the solution u satisfies the decay estimate max(1,

t 3/4 ) u(t, •) L 2 x (R) + max(t 1/2 , t) u(t, •) L ∞ x (R) ≤ C.
However, this result is proved assuming a smallness condition on the initial data. We relax here this assumption first for non negative initial data u 0 ≥ 0, and then conclude by a comparison result. Our main result states as follows

Theorem 1.1 Assume u 0 in L 1 (R) ∩ L 2 (R).
Then there exists C(u 0 ) such that for any t ≥ 0, max(1,

t 3/4 ) u(t, •) L 2 x (R) + max(t 1/2 , t) u(t, •) L ∞ x (R) ≤ C(u 0 ). (1.2)
We provide here a complete proof that relies on an ad hoc representation of the solution.

A diffusive representation of the solution

Introducing a new variable σ ≥ 0, for a given function u(t, x) that depends on t and on x (the variable x will be omitted in the remaining of this section) we introduce

ϕ(t,σ) = t 0 u(s) exp(-σ 2 (t -s))ds. (1.3)
Then the half-integral of a function u reads

Iu(t) = 1 √ π t 0 u(s) √ t -s ds = Ŵ(ϕ)(t) = 2 π +∞ 0 ϕ(t,σ)dσ. (1.4)
Hence equation (1.1) reads also

u t + u x + Ŵ(ϕ t ) + uu x = u xx , ϕ t + σ 2 ϕ = u, ( 1.5)
supplemented with initial data (u 0 , 0). We emphasize (see [START_REF] Manoubi | Theoretical and numerical analysis of the decay rate of solutions to a water wave model with a nonlocal viscous dispersive term with Riemann-Liouville half derivative[END_REF]) that there is no uniqueness for the diffusive representation. The representation above is convenient for our purpose.

Outline of the proof and mathematical framework

The proof of Theorem 1.1 relies on two arguments. The first one is a comparison principle between solutions. The second one is the proof of the main theorem for non negative solutions.

Throughout the article, the variable x is omitted for the functions u and ϕ.W eset

|u| 2 L 2 = R u 2 (x)dx,
for functions of one variable and

||ϕ|| 2 L 2 = R×[0,+∞) ϕ 2 (σ, x)dxdσ,
for functions of two variables. Here the third variable t is omitted.

2 Proof of the main theorem

A comparison result

Our first milestone is to prove a comparison result between solutions.

Proposition 2.1 Consider u 0 ≤ v 0 two ordered initial data that satisfy the assumptions of Theorem 1.1,i.e.u 0 ,v 0 in L 1 (R) ∩ L 2 (R). We have that u(t, x) ≤ v(t, x) for any t > 0 and a.e. x ∈ R.

Proof

We prove the comparison result first for initial data in the Schwartz's class S(R) and then the final result for general initial data by a limiting argument. This limiting argument is standard and then omitted for the sake of conciseness.

We now consider u 0 (x)<v 0 (x) two (strictly) ordered initial data in the Schwartz class. Consider u,vthe corresponding solutions and ϕ u ,ϕ v the associated diffusive representation of solutions. Consider w = vu and ψ = ϕ vϕ u . Introduce

T = inf t > 0; inf x w(t, x)<0 . (2.1)
We prove that T =+∞arguing by contradiction. Assume that T < +∞. There exists x 0 in R such that inf x w(T , x) = w(T , x 0 ) = 0andthenw x (T , x 0 ) = 0, w xx (T , x 0 ) ≥ 0. Since the solution is smooth away from the boundary t = 0wealsohavew t (T , x 0 ) ≤ 0.

On the other hand we infer from the first equation in (1.5)that

Ŵ(ψ t )(T , x 0 ) =-w t (T , x 0 ) + w xx (T , x 0 ) + w x (T , x 0 ) 1 + u(T , x 0 ) + v(T , x 0 ) 2 + w(T , x 0 ) u x (T , x 0 ) + v x (T , x 0 ) 2 ≥ 0. (2.2) Besides, since w(t, x) ≥ 0f o rt ≤ T ,t h e nψ(T , x,σ) ≥ 0. Therefore the second equation in (1.5) yields ψ t (T , x 0 ,σ) =-σ 2 ψ(T , x 0 ,σ) ≤ 0f o ra n yσ ≥ 0. Gathering this inequality with (2.2)w eh a v et h a tŴ(ψ t )(T , x 0 ) = 0. Therefore Ŵ(σ 2 ψ)(T , x 0 ) = 0, and then ψ(T , x 0 ,σ) = 0f o ra n yσ . This yields w(t, x 0 ) = 0f o ra n yt ≤ T then the contradiction. ⊓ ⊔ Choosing 0 ≤ u 0 in Proposition 2.1 yields Corollary 2.2 Assume u 0 in L 1 (R) ∩ L 2 (R)
and u 0 ≥ 0 then u(t, x) ≥ 0 for almost every x ∈ R and t > 0.

Decay of the

L 1 norm Proposition 2.3 Assume u 0 in L 1 (R) ∩ L 2 (R).W ehave |u(t)| L 1 ≤ min 1, 1 √ t |u 0 | L 1 . (2.3)
Proof In fact we just have to prove the result for

u 0 ≥ 0S i n c e u 0 -|u 0 | 2 ≤ u 0 ≤ u 0 +|u 0 | 2
, then the corresponding solutions starting from these three initial data remain ordered, say v ≤ u ≤ w. We prove the decay estimate for -v and w that are non negative, and conclude.

Assume then u 0 ≥ 0. We now know that the solution remains non negative for all time. Integrating in x (and in t)theEq. (1.1) we have, setting

m(t) = R u(t, x)dx that is the L 1 norm m(t) + Im(t) = m(0). (2.4)
Appealing Proposition 1 in [START_REF] Manoubi | Theoretical and numerical analysis of the decay rate of solutions to a water wave model with a nonlocal viscous dispersive term with Riemann-Liouville half derivative[END_REF] we know that the solution of (2.4) reads

m(t) = +∞ t e t-s √ s ds m(0).
This completes the proof of the Proposition. ⊓ ⊔

Revisiting the Schonbek's splitting method

Consider a solution (u,ϕ)to (1.5). Introduce the energy functional

E(t) =|u(t)| 2 L 2 + 2 π ||σϕ|| 2 L 2 . (2.5)
Multiplying the first equation in (1.5)b yu and the second one by 2 π ϕ t and summing the resulting equations lead to the energy dissipation mechanism

d dt E(t) + 2|u x (t)| 2 L 2 + 4 π ||ϕ t || 2 L 2 = 0. ( 2
tE(t) ≤ c |u 0 | 2 L 1 + t 0 |u(s)| 4 3 L 2 ds 3 2
.

(2.7)

Proof We begin with the identity

d dt (tE(t)) = E(t) -2t|u x | 2 L 2 - 4t π ||ϕ t || 2 L 2 . (2.8)
Dropping the last term, we infer from this

d dt (tE(t)) ≤ (|u(t)| 2 L 2 -2t|u x | 2 L 2 ) + 2 π ||σϕ|| 2 L 2 .
(2.9)

We now bound by above the last term in the right hand side of (2.9). Due to Minkowski's inequality, equality (1.3), and Fubini-Tonelli theorem we have

||σϕ|| L 2 ≤ t 0 |u(s)| L 2 ∞ 0 σ 2 exp(-2σ 2 (t -s))dσ 1 2 ds (2.10) using that ∞ 0 x 2 exp(-ax 2 )dx = √ π 4a 3 2
,wededucefrom(2.10)

||σϕ|| L 2 ≤ π 1 4 2 -7 4 t 0 |u(s)| L 2 (t -s) 3 4
ds.

(2.11)

We now tackle the first term in the right hand side of (2.9). Proceeding as in the classical Schonbek's argument we have, appealing Plancherel's theorem and the inverse Fourier transform

|u(t)| 2 L 2 -2t|u x | 2 L 2 = 1 2π R (1 -2t|ξ | 2 )|û(t,ξ)| 2 dξ ≤ 1 2π 2t|ξ | 2 ≤1 |û(t,ξ)| 2 dξ ≤ 4π|u(t)| 2 L 1 1 √ 2t 0 dξ.
(2.12)

Combining this with the L 1 estimate (2.3) leads to

|u(t)| 2 L 2 -2t|u x | 2 L 2 ≤ 2 √ 2π |u 0 | 2 L 1 √ t min 1, 1 t . (2.13)
It is worth to point out that the right hand side of (2.13) is integrable for t ∈[ 0, +∞).

Gathering these inequalities and integrating in time (2.9) yields [START_REF] Willem | Analyse Harmonique Réelle[END_REF]) that reads: there exists a constant c HLS (see [START_REF] Lieb | Analysis[END_REF] for the precise value of this constant) such that

tE(t) ≤ 8 √ 2π|u 0 | 2 L 1 + π -1 2 2 -5 2 t 0 τ 0 |u(s)| L 2 (τ -s)
|| R |g(s)||s -t| -3 4 ds|| L 2 t ≤ c 2 HLS ||g|| 2 L 4 3
.

(2.15)

Applying this to the function that is g(s

) =| u(s)| L 2 if s ∈ (0, t),a n dg(s) = 0 elswhere leads to t 0 τ 0 |u(s)| L 2 (τ -s) 3 4 ds 2 dτ ≤ c 2 HLS t 0 |u(s)| 4 3 L 2 ds 3 2
.

(2.16)

This completes the proof of the Lemma.

⊓ ⊔

We now prove the next milestone

Corollary 2.5 Assume that u 0 in L 1 (R) ∩ L 2 (R). There exists a numerical constant c such that √ tE(t) ≤ c|u 0 | 2 L 1 ∩L 2 .
Proof Here and in the sequel

|u 0 | L 1 ∩L 2 =|u 0 | L 1 +|u 0 | L 2 .
We begin with an a priori estimate on the L 2 norm of the solution. Appealing Gagliardo-Nirenberg inequality (see [START_REF] Brezis | Analyse Fonctionnelle[END_REF])

|u| L 2 ≤|u| 2 3 L 1 |u x | 1 3 L 2 .
(2.17)

Then using the Hölder's inequality

t 0 |u(s)| 4 3 L 2 ds ≤ t 0 |u(s)| 8 7
L 1 ds

7 9 t 0 |u x | 2 L 2 ds 2 9
.

(2.18)

Therefore, since the last term in the right hand side is bounded by

C(u 0 ) = ( |u 0 | 2 L 2
2 )

2 9 thanks to (2.6), using (2.3) )ds) Using the results in [START_REF] Manoubi | Theoretical and numerical analysis of the decay rate of solutions to a water wave model with a nonlocal viscous dispersive term with Riemann-Liouville half derivative[END_REF] we know that there exists two functions K (t)(x) and N (t)(x) such that, where * stands for the convolution in the x variable

t 0 |u| 4 3 L 2 ≤ 2 -2 9 |u 0 | 4 9 L 2 |u 0 | 8
u(t) = K (t) * u 0 + t 0 N (t -s) * u 2 (s)ds.
(2.21)

We also have that there exists c a numerical constant such that (see [START_REF] Manoubi | Theoretical and numerical analysis of the decay rate of solutions to a water wave model with a nonlocal viscous dispersive term with Riemann-Liouville half derivative[END_REF])

max(1, √ t) √ t|K (t)| L ∞ +|K (t)| L 1 + √ t|N (t)| L 1 + t|N (t)| L ∞ ≤ c. (2.22) Therefore |u(t)| L ∞ ≤|K (t)| L ∞ |u 0 | L 1 + t 0 |N (t)| L 3 |u(s)| 2 L 3 ds ≤ c|u 0 | L 1 min( 1 √ t , 1 t ) + c t 0 (t -s) -5 6 |u(s)| 4 3 L 2 |u(s)| 2 3
L ∞ ds.

(2.23)

Introducing the auxiliary function A(t) = √ t|u(t)| L ∞ we then infer from (2.23) and from the L 2 estimate in Corollary 2.5

A(t) ≤ c|u 0 | L 1 + c|u 0 | 4 3 L 1 ∩L 2 t 0 √ tds (t -s) -5 6 s -2 3 (sup s≤t A(s)) 2 3 .
(2.24)

We easily infer from (2.24)thatA(t) ≤ c|u 0 | 2 L 1 ∩L 2 that is the good decay estimate for small t, i.e. |u 0 | L 1 ∩L 2 t -1 2 . We now proceed to the decay rate for large t for the L 2 norm of the solution. Interpolating the L 1 and the L ∞ decay rate we have that

|u(t)| 2 L 2 ≤|u(t)| L 1 |u(t)| L ∞ ≤ c|u 0 | 2 L 1 ∩L 2 t -1 ; (2.25)
we improve this new estimate. We infer from (2.21)-(2.22)

|u(t)| L 2 ≤|K (t)| L 2 |u 0 | L 1 + t 0 |N (t)| L 2 |u(s)| 2 L 2 ds ≤ c|u 0 | L 1 t -3 4 + c t 0 (t -s) -3 4 |u(s)| 2 L 2 ds.
(2.26)

The first term in the right hand side of this inequality has the good decay rate. We focus on the second one. Due to Corollary 2.5 and using that |u(t)| L 2 remains bounded

t 0 (t -s) -3 4 |u(s)| 2 L 2 ds ≤c|u 0 | 2 L 1 ∩L 2 t 0 (t -s) -3 4 min
B(t) ≤ c |u 0 | L 1 + t 1 0 |N (t -s)| L ∞ |u(s)| 2 L 2 ds +t t 1 |N (t -s)| L 2 |u(s)| L 2 |u(s)| L ∞ ds .
(2.30)

The second term in the right hand side of (2.30) is bounded by above by c t t-1 |u 0 | 2 L 2 and then bounded for large t. The third term in the right hand side of (2.30) is bounded by above by c|u 0 | 2

L 1 ∩L 2 t 1 tds (t-s) 3 4 s 5 4
ds,usingtheL 2 estimate in Theorem 1.1 and the previous bound in Ct -1 2 for the L ∞ norm of u. This completes the proof of the theorem 3 Numerical results

The numerical scheme

Here we consider the generalized model

u t + u x + √ ν √ π ∂ ∂t t 0 u(s) √ t -s ds + γ uu x = αu xx , (3.1)
where ν, γ and α are non-negative parameters. Using the diffusive realization of the half-order derivative (1.5), the model (3.1) is rewritten as a PDE-ODE coupled system as follows

⎧ ⎨ ⎩ ∂ t u(t, x) + ∂ x (u + γ 2 u 2 ) =-2 √ ν π ∞ 0 (u(t, x) -σ 2 ϕ(t, x,σ))dσ +αu xx (t, x), t > 0, x ∈ R ∂ t ϕ(t, x,σ) =-σ 2 ϕ(t, x,σ)+ u(t, x), t > 0, x ∈ R,σ ≥ 0, (3.2) supplemented with the initial conditions ∀x ∈ R, ∀σ ≥ 0,ϕ ( 0, x,σ) = 0, ∀x ∈ R, u(0, x) = u 0 (x).
Following [START_REF] Dumont | Numerical analysis of a water wave model with a nonlocal viscous dispersive term using the diffusive approach[END_REF], we use the Gauss-Jacobi quadrature method with N m points to approximate the generalized integral in (3.2). We note by w i the weights and by σ i the nodes of the Gauss-Jacobi quadrature method. We get

D 1/2 u(t, x) ≃ N m i=1 w i 2 π u(t, x) -σ 2 i ϕ(t, x,σ i ) = N m i=1 w i 2 π u(t, x) -σ 2 i ϕ i (t, x) .
We note by U = (u,ϕ 1 ,...,ϕ N m ) T the vector of (N m + 1) unknowns, by

F (U ) = (u + γ 2 u 2 , 0,...,0) T ,
and finally

S = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ -2 √ ν π N m i=1 w i 2 √ ν π w 1 σ 2 1 ••• ••• 2 √ ν π w N m σ 2 N m 1 -σ 2 1 0 ••• 0 10 -σ 2 2 . . . 0 . . . . . . . . . 10 ••• ••• -σ 2 N m ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ .
Thus, the problem (3.2) is written in the following form

∂ t U + ∂ x F (U ) = S(U ) + G 1 ∂ 2 x U , (3.3)
where G 1 is a diagonal matrix of order N m + 1. Following [START_REF] Dumont | Numerical analysis of a water wave model with a nonlocal viscous dispersive term using the diffusive approach[END_REF], we use a splitting method to provide numerical results. Let t > 0, for all n ≥ 0, we recall that t n = n t and

U n (x) ≈ U (n t, x).
From equation (3.3), we consider the propagation equation

∂ t U + ∂ x F (U ) = G 1 ∂ 2 x U , (3.4)
and the diffusive equation

∂ t U = S(U ). (3.5)
We note by H a (respectively, H b ) the discrete operator of the solution of (3.4) (respectively, the solution of (3.5)). Then a Strang splitting method of order 2 (see for example [START_REF] Holden | Operator splitting for the KDV equation[END_REF][START_REF] Leveque | Numerical methods for conservation laws[END_REF]) between t n and t n+1 is used to solve respectively (3.4)and(3.5) as follows (1) ,

U (1) = H b t 2 U n , U (2) = H a ( t) U
U n+1 = H b t 2 U (2) . (3.6)
Here, the constructed operators H a and H b are stable and of order 2. Then, the scheme (3.6) provides an approximation of order 2 in time to the problem (3.3). Now, let t ∈[ 0, T ],

x ∈[ 0, L]. We consider periodic boundary conditions i.e u| x=0 = u| x=L .Let t > 0bea fixed time step and h > 0 be a fixed space step. We set for all n = 0,...,N , t n = n t and for all j = 1,...,K , x j = jh. We use the following notations: For the propagation Eq. (3.4), we use a Crank-Nicolson scheme for the linear part and Adams-Bashforth scheme for the nonlinear part. For the space discretization, we use standard Fourier methods. Hence, the scheme is given by

∀ i = 1,...,N , u n (x) ≈ u(n t, x), ψ n i (x) ≈ ψ(n t, x,σ i ).
ûn+1 -û n t = ûn+1 +û n 2 (-αξ 2 -iξ) - iγξ 2 (3 (u n ) 2 -(u n-1 ) 2 ) 4 . (3.7)
Then, for the diffusion equation (3.5), the exact solution is given by

H b t 2 U = e S t 2 U . (3.8)
We note that the scheme (3.7) as well as the Strang splitting method are of order two. It follows that (3.6)isoforder2.

Numerical results

In the sequel, we present the numerical results. To this end, we fix the parameters values to ν = 1andα = γ = 0.5. This case is theoretically investigated here and in [START_REF] Manoubi | Theoretical and numerical analysis of the decay rate of solutions to a water wave model with a nonlocal viscous dispersive term with Riemann-Liouville half derivative[END_REF]. Let us notice that in this article, contrarily to [START_REF] Dumont | Numerical analysis of a water wave model with a nonlocal viscous dispersive term using the diffusive approach[END_REF], the dispersive term βu xxx is not taken into account as in the present theoretical study and we have shown in [START_REF] Dumont | Numerical analysis of a water wave model with a nonlocal viscous dispersive term using the diffusive approach[END_REF] that the magnitude of the constant β has a strong influence on the decay rate of the solution for large times.

Unfortunately, an exact solution of (3.1) is not known, but a numerical study of the convergence of the scheme is provided in [START_REF] Manoubi | Theoretical and numerical analysis of the decay rate of solutions to a water wave model with a nonlocal viscous dispersive term with Riemann-Liouville half derivative[END_REF]. Then, in order to numerically investigate the decay rate of the solutions for large times, we choose an initial data u 0 that can provide an exact KdV soliton for well chosen parameters (see [START_REF] Chen | Exact travelling-wave solutions to bidirectional wave equations[END_REF]). Moreover, we choose an interval, both in time and in space, such that the solution remains small at the boundary in order to have a numerical solution (on a finite domain) that correctly approximates the exact solution that exists on the whole real line R. Actually, we have that the initial data u 0 is in the space L 1 (R) ∩ L 2 (R).

We denote by x 0 the middle of the space interval. Then the initial data reads u 0 (x) = 0.32 * sech 2 (0.4 * (xx 0 )).

(3.9) We run the computations for the numerical scheme (3.7)-(3.8) on the space interval [0, L] with L = 800. Moreover, we take the space step h = 0.1a n dN m = 20. In order to numerically confirm the order of the convergence in time of the scheme (3.7)-(3.8), we take T = 100, we denote by u n Ref the reference solution when t = 0.05 and by u n the numerical solution for different time steps t. We denote by E n ( t) = u n Refu n 2 the L 2 -error due to the time discretization of the solution. We illustrate in Fig. 1 the error E n in terms of t in the log-log scale. We observe that, as expected, the error E n decreases when t decreases. Moreover, the measured values are close to a straight line with slope 2.24. This test confirms that the numerical scheme (3.7)-(3.8) is of order 2 in time. In the sequel, we provide our simulations with the time step t = 0.1.

In order to check the convergence of the scheme (3. Refu n 2 the L 2 -error in terms of the quadrature points. We display the results in Fig. 2. As expected, we see that the error decreases when the number of grid points N m increases. In addition, we observe that the measured values are close to a straight line with slope -11. Moreover, we observe that the error E n is less that 10 -6 from N m = 15. Other numerical simulations are done and we made similar observations (see also [START_REF] Lombard | Numerical modeling of nonlinear acoustic waves in a tube connected with Helmholtz[END_REF][START_REF] Manoubi | Theoretical and numerical analysis of the decay rate of solutions to a water wave model with a nonlocal viscous dispersive term with Riemann-Liouville half derivative[END_REF]). Hence, we will use 15 grid points in the remaining of this article.

In the sequel, we use the numerical scheme (3.7)-(3.8) with the previous parameters values on a large time interval [0, T ]=[0, 1000]. Numerical results are displayed in Figs. 3 and4. It can be observed that the numerical solution remains non negative. This matches the results of Corollary 2.2. Moreover, the decay rates of the numerical solution provided in Fig. 4 in norm L 2 and in norm L ∞ matches the theoretical results of Theorem 1.1 very precisely. Indeed, one can observe that the relative error between the theoretical and the numerical L 2 decay rate (diffusion) is less than 0.53% whereas it is less than 0.8% the the decay rate in the L ∞ -norm (dispersion).

Fig. 1

 1 Fig. 1 Errorofthetime discretization

Fig. 2 11 Fig. 3

 2113 Fig. 2 Error in terms of grid points N m

  7)-(3.8) in terms of the Gauss-Jacobi quadrature points N m ,wedenotebyu n Ref the reference solution when N m = 20 and by u n the numerical solution for different values of N m .W edenotebyE n (N m ) = u n

Fig. 4 a

 4 Fig. 4 a L 2 -decay rate of the solution, b L ∞ -decay rate of the solution when ν = 1, α = γ = 0.5a n d h = t = 0.5

  .6) This implies that t →| u(t)| L 2 remains bounded by |u 0 | L 2 .W eno wstateandpro v eafirst result

	Lemma 2.4 There exists a numerical constant c such that for any t ≥ 0

4 Completing the proof of the theorem

  We already know that t →|u(t)| 2 L 2 is bounded by |u 0 | 2 L 2 . We know from Corollary 2.5 that for large t this function decays as t -1 4 . Before improving this decay rate, we prove a first decay estimate for the L ∞ norm of the solution.

	7 9 ≤ ( 7 3 )	7 9 t	1 3 , plugging this into Lemma 2.4 leads to

tE(t) ≤ c|u 0 | 2 L 1 ∩L 2 (1 + √ t); (2.

20) Since we know by (2.6)thatE(t) ≤| u 0 | 2 L 2 , gathering this with (2.20) completes the proof of the Corollary. ⊓ ⊔ 2.

  We now get rid of the extra term log t. From the previous extimates we have improved(2.25)a s|u(t)| L 2 ≤ C(u 0 )t -5 8 for instance. Using this in (2.27) and performing the same computations as in (2.27)-(2.29) leads to the good decay estimate (without the extra log t term). It remains to prove the large t estimate for the L ∞ norm of the solution. We set B(t) = t|u(t)| L ∞ . We infer from (2.21)-(2.22)

	Hence for t ≥ 2 there exists a numerical constant c such that
			0	1	(t -s) -3 4 ds +	1	t	(t -s) -3 4	ds s	≤ c	1 + log t 4 t 3	.	(2.29)
																		1,	1 s	ds.	(2.27)
	This last integral splits as (for t ≥ 1)
	0	1	(t -s) -3 4 ds +	1	t	(t -s) -3 4	ds s	≤	1 (t -1)	3 4	+ t -3 4	1 t	1	ds s(1 -s)	4 3	.	(2.28)

Acknowledgements These result were presented at the PACOM 2017 conference. The authors would like to thank the anonymous referees for their remarks and suggestions.