Farshid Dabaghi

Adrien Petrov

Jérôme Pousin

Yves Renard

A robust finite element redistribution approach for elastodynamic contact problems

Keywords: Numerical solution, mass redistribution method, variational inequality, unilateral contact, energy conservation. 2008 MSC: 35L05, 35L85, 49A29, 65N10, 65N30, 74M15

This paper deals with a one-dimensional elastodynamic contact problem and aims to highlight some new numerical results. A new proof of existence and uniqueness results is proposed. More precisely, the problem is reformulated as a differential inclusion problem, the existence result follows from some a priori estimates obtained for the regularized problem while the uniqueness result comes from a monotonicity argument. An approximation of this evolutionary problem combining the finite element method as well as the mass redistribution method which consists on a redistribution of the body mass such that there is no inertia at the contact node, is introduced. Then two benchmark problems, one being new with convenient regularity properties, together with their analytical solutions are presented and some possible discretizations using different time-integration schemes are described. Finally, numerical experiments are reported and analyzed.

Introduction

This paper aims to give some new numerical results for the dynamical evolution with impact of a linearly elastic body and a rigid obstacle by using Email addresses: farshid.dabaghi@insa-lyon.fr (Farshid Dabaghi), apetrov@math.univ-lyon1.fr (Adrien Petrov), jerome.pousin@insa-lyon.fr (Jérôme Pousin), Yves.Renard@insa-lyon.fr (Yves Renard) the mass redistribution method introduced in [START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF]. The situations involving contact abound in industrial or biomedical problems. For this reason a considerable engineering and mathematical literature has been devoted to contact problems. Usually, with the hypotheses of small deformations, the contact phenomena is modelled by using the so-called Signorini's boundary conditions in displacement, which are based on a linearization of the physically meaningful non penetrability of masses. The elastodynamic contact problems are typically stated as hyperbolic type variational inequalities in Sobolev spaces, and a few existence of solutions results has been established. The reader is referred to [START_REF] Schatzman | A hyperbolic problem of second order with unilateral constraints: The vibrating string with a concave obstacle[END_REF][START_REF] Lebeau | A wave problem in a half-space with a unilateral constraint at the boundary[END_REF][START_REF] Kim | A boundary thin obstacle problem for a wave equation[END_REF] and also the comprehensive monograph [START_REF] Eck | Unilateral contact problems : variational methods and existence theorems[END_REF]. The uniqueness result was successfully investigated for a wave equation in a half-space with a unilateral contact at the boundary in [START_REF] Lebeau | A wave problem in a half-space with a unilateral constraint at the boundary[END_REF]. In particular, the uniqueness still remains an open problem in other frameworks. The fundamental mathematical difficulties are related to the intrinsic nonsmoothness of the problem coming from the Signorini boundary conditions. Finally, notice that the existence of a strong solution for viscoelastodynamic problems with unilateral constraints is proved in [START_REF] Petrov | Mathematical results on existence for viscoelastodynamic problems with unilateral constraints[END_REF].

On the other hand, another challenging task consists to elaborate efficient numerical methods for solving dynamical contact problems. Newmark's method is one of the most popular numerical solvers which is usually combined with finite elements methods. Unfortunately, for almost all values of the parameters of this method, unphysical energy blow-ups occur during the time integration as well as numerical instabilities at the dynamical contact boundaries. To overcome these difficulties, some numerical methods based on the Newmark scheme for solving impact problems and which are energy dissipative have been introduced in [START_REF] Carpenter | Lagrange constraints for transient finite element surface contact[END_REF]. Although stable, those methods lead to an important energy loss at impact which do not vanish when the time step decreases. Laursen et al. in [START_REF] Laursen | Improved implicit integrators for transient impact problems-geometric admissibility within the conserving framework[END_REF][START_REF] Laursen | Design of energy conserving algorithms for frictionless dynamic contact problems[END_REF][START_REF] Chawla | Energy consistent algorithms for frictional contact problems[END_REF][START_REF] Laursen | Computational contact and impact mechanics. Fundamentals of modeling interfacial phenomena in nonlinear finite element analysis[END_REF] have designed time integration schemes of Newmark type which are energy conserving. However, these schemes are unable to prevent some spurious oscillations of the displacement and of the contact stress on the contact boundary. These unphysical oscillations are avoided by the numerical methods developed in [START_REF] Deuflhard | A contact-stabilized Newmark method for dynamical contact problems[END_REF] (see also the reference therein), but these methods are still energy dissipative. Another approach consisting to remove the mass at the contact nodes is investigated in [START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF]. This mass redistribution method prevents the oscillations at the contact boundaries mentioned above, and it can be proven that the solution is unique and energy conserving, see [START_REF] Dabaghi | Convergence of mass redistribution method for the one-dimensional wave equation with a unilateral constraint at the boundary[END_REF]. Furthermore, some different variants has been proposed in [START_REF] Hager | A stable energy-conserving approach for frictional contact problems based on quadrature formulas[END_REF] with an estimate of the space approximation for linear elasticity. It has also been extended for thin structures in [START_REF] Renard | The singular dynamic method for constrained second order hyperbolic equations: application to dynamic contact problems[END_REF][START_REF] Pozzolini | Some energy conservative schemes for vibroimpacts of a beam on rigid obstacles[END_REF][START_REF] Tkachuk | Hybrid-mixed discretization of elasto-dynamic contact problems using consistent singular mass matrices[END_REF] considering different discretizations for the displacement and the velocity.

The paper is organised as follows. In Section 2, the mathematical models are presented and a new proof of existence and uniqueness results is proposed. A space semi-discretization based on a redistribution of mass is introduced in Section 3 and some time integration method are advertised. In Sections 4 and 5, two elastodynamic contact benchmark problems with their analytical solutions are exhibited. Note that an analytical solution associated to the first benchmark problem is expressed by using Taylor expansion (see [START_REF] Wriggers | Computational contact mechanics[END_REF][START_REF] Doyen | Time-integration schemes for the finite element dynamic Signorini problem[END_REF]) and it is compared to the approximated solution obtained thanks to Crank-Nicolson method. However, the numerical experiments highlight that the comparison between the analytical and approximated solutions do not converge. Therefore a new analytical solution displayed for the second benchmark is exposed in Section 5 and some numerical experiments problem are reported like convergence rates for the displacement, Lagrange multiplier and energy evolution showing the efficiency of mass redistribution method.

Description of the models

We consider a model for an elastic bar of length L vibrating longitudinally. One end of this bar is free to move, as long as it does not hit a material obstacle, while the other end can be clamped or free to move. The obstacle constrains the displacement of the extremity to be non negative. The material of the bar is supposed to be homogeneous and the theory of small deformations is considered. Let us denote by x the spatial coordinate along the bar, with the origin at the material obstacle, let u(x, t) be the displacement at time t ∈ [0, T ], T > 0 of the material point of spatial coordinate x ∈ [0, L]. Let f (x, t) denotes a density of external forces, depending on time and space. In the case where the bar is clamped at one end (see Figure 5), the mathematical problem is formulated as follows:

u tt (x, t) -u xx (x, t) = f (x, t), (x, t) ∈ (0, L) × (0, T ), (1)
with Cauchy initial data

u(x, 0) = u 0 (x) and u t (x, 0) = v 0 (x), x ∈ (0, L), (2)
and Signorini and Dirichlet boundary conditions at x = 0 and x = L, respectively,

0 ≤ u(0, t) ⊥ u x (0, t) ≤ 0 and u(L, t) = 0, t ∈ [0, T ]. (3)
While in the case where the bar is free to move at one end (see Figure 1), the motion is governed by the momentum equilibrium equation ( 1) together with Cauchy initial data (2) and Signorini and Neumann boundary conditions at x = 0 and x = L, respectively,

0 ≤ u(0, t) ⊥ u x (0, t) ≤ 0 and u x (L, t) = 0, t ∈ [0, T ]. (4)
Here u t def = ∂u ∂t and u x def = ∂u ∂x . The orthogonality has the natural meaning; namely if we have enough regularity, it means that at the boundary x = 0 the product u(0, •)u x (0, •) vanishes almost everywhere in time. If it is not the case, the above inequality is integrated on an appropriate set of test functions, leading to a weak formulation for the unilateral condition. Observe that from mathematical viewpoint, Signorini's condition means that when the bar touches the obstacle at x = 0, its reaction can be only upwards, so that u x (0, •) ≤ 0 on the set {t : u(0, •) = 0}. While in the case where the bar does not touch the obstacle, its end is free to move. More precisely, we have u x (0, •) = 0 on the set {t : u(0, •) > 0}. We suppose that the initial displacement u 0 belongs to the H 1 (0, L) and satisfies the compatibility conditions, i.e. u 0 (L) = 0 and u 0 (0) ≥ 0, the initial velocity v 0 belongs to L 2 (0, L) and the density of external forces f belongs to L 2 (0, T ; L 2 (0, L)).

We describe now the weak formulations of both the two problems. For that purpose, it is convenient to introduce the following notations:

V Dir def = {u ∈ H 1 (0, L) : u(L) = 0}, V Neu def = H 1 (0, L) and H def = L 2 (0, L).
We denote by K the following convex set:

K def = {u ∈ H 2 : u(•, t) ∈ K for almost every t}, where H 2 def = {u ∈ L 2 (0, T ; V) : u t ∈ L 2 (0, T ; H)}, K def = {u ∈ V : u(0) ≥ 0} and V = V Dir or V = V Neu depending on the conditions taken in x = L.
Therefore the weak formulations associated to (1)-(3) and to (1)-(2) together with (4) are obtained by multiplying (1) by vu and by integrating formally this result over

Q T def = (0, L) × (0, T ) to get                  find u ∈ K such that - L 0 v 0 (v(•, 0)-u 0 ) dx- Q T u t (v t -u t ) dx dt + Q T u x (v x -u x ) dx dt≥ Q T f (v-u) dx dt for all v ∈ K for which there exists ζ > 0 with v = u for t ≥ T -ζ. (5)
Notice that the sole difference of the weak formulation associated to (1)-(3) and to (1)-( 2), (4) comes from the definition of the convex set K.

It should also be noted that existence and uniqueness results are obtained for a similar situation of a vibrating string with concave obstacle in one dimensional space in [START_REF] Schatzman | A hyperbolic problem of second order with unilateral constraints: The vibrating string with a concave obstacle[END_REF] and also for a wave equations with unilateral constraint at the boundary in a half space of R N in [START_REF] Lebeau | A wave problem in a half-space with a unilateral constraint at the boundary[END_REF]. An existence result for a wave equation in a C 2 -regular bounded domain constrained by an obstacle at the boundary in R N for N ≥ 2 is proven in [START_REF] Kim | A boundary thin obstacle problem for a wave equation[END_REF].

Theorem 2.1 (Existence and uniqueness results). Assume that

u 0 ∈ V, v 0 ∈ H and f ∈ L 2 (0, T ; H). Then Problem (1)-(3) admits a unique solution u ∈ L ∞ (0, T ; V) ∩ W 1,∞ (0, T ; H).
Proof. The first step consists to introduce an auxiliary problem. Namely, let ū be the solution to Problem (1) with initial data (2) and boundary conditions (3) where Signorini's boundary condition is replaced by Dirichlet's one, i.e. ū(0, t) = 0, t ∈ [0, T ]. Note that the auxiliary problem possesses a unique solution ū belonging to C 0 ([0, T ]; [22,[START_REF] Brezis | Analyse fonctionnelle[END_REF]). Then we denote by v def = uū a solution of v tt (x, t)v xx (x, t) = 0, (x, t) ∈ (0, L) × (0, T ), [START_REF] Deuflhard | A contact-stabilized Newmark method for dynamical contact problems[END_REF] with Cauchy initial data v(x, 0) = 0 and v t (x, 0) = 0, x ∈ (0, L),

H 2 (0, L) ∩ H 1 0 (0, L)) ∩ C 1 ([0, T ]; H 1 0 (0, L)) ∩ C 2 ([0, T ]; H) (see
and Signorini and Dirichlet boundary conditions at x = 0 and x = L, respectively,

0 ≤ v(0, t) ⊥ v x (0, t) + ūx (0, t) ≤ 0 and v(L, t) = 0, t ∈ [0, T ]. ( 8
)
The second step consists to rewrite ( 6)-( 8) as a differential inclusion problem by using the characteristic method. To this aim, it is convenient to introduce the characteristic coordinates ξ def = x + t and η def = xt. Therefore the chain rule gives v xx = v ξξ + 2v ξη + v ηη and v tt = v ξξ -2v ξη + v ηη . According to (1), v ξη vanishes and it follows that v(ξ, η) = φ(ξ) + ψ(η), where φ and ψ are two differentiable functions such that v(x, t) = φ(x+t) + ψ(x-t). In particular, taking t = 0 and using the Cauchy initial data, we get φ(x) + ψ(x) = 0 and φ (x)ψ (x) = 0, which by integration implies that there exist two constants C φ and C ψ such that φ(ξ) = C φ and ψ(η) = C ψ for all ξ and η belonging to [0, L]. It comes that the boundary conditions (3) can be rewritten as follows: 0 ≤ φ(t) + ψ(-t) ⊥ φ (t) + ψ (-t) + g(t) ≤ 0 and φ(L+t) + ψ(L-t) = 0 (9) for all t belonging to [0, T ]. Thanks to the above identity, we may extend φ(t) for all t ∈ [L, 2L], i.e. we have

φ(L+t) = -ψ(L-t)
for all t belonging to [0, L]. If we choose t = L+t, we get φ(t ) = -ψ(2L-t ). We already have the solution for ψ(t) with 0 ≤ t ≤ L and if L ≤ t ≤ 2L, we can obtain φ(t ) by observing that 0 ≤ 2Lt ≤ L and by using ψ(2L-t), it comes that φ(t) = C ψ for all t belonging to [L, 2L].

Let the indicator function, namely χ (-∞,0] (x) is equal to 0 if x ∈ (-∞, 0] and +∞ otherwise. Obviously, χ (-∞,0] is a lower semi-continuous and convex function, for further details the reader is referred to [START_REF] Brézis | Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert[END_REF]. Its subdifferential is a multivalued function J N : R → P( R) defined by

J N [x] def =      {0} if x < 0, [0, +∞) if x = 0, ∅ if x > 0,
where P( R) is the set of all subsets of R. Then, the inequalities in (9) can be rewritten as follows

φ (t) + ψ (-t) + g(t) ∈ -J N [-φ(t)-ψ(-t)]. (10)
Note that at this stage, ψ(η), η ∈ [-L, L], is the unique unknown of [START_REF] Grosu | Stability of semidiscrete formulations for elastodynamics at small time steps[END_REF].

Let us define now w(t

) def = -φ(t) -ψ(-t). (11)
We insert [START_REF] Hager | A stable energy-conserving approach for frictional contact problems based on quadrature formulas[END_REF] into [START_REF] Grosu | Stability of semidiscrete formulations for elastodynamics at small time steps[END_REF] to get

w (t) ∈ -J N [w(t)] -g(t).
Finally, we find the following Cauchy problem

w (t) ∈ -J N [w(t)] -g(t) a.e. t ∈ (0, L), (12a) w(0) = 0. ( 12b
)
Then we look for a solution of the following problem:

w (t) + h(t) = -g(t) a.e. t ∈ (0, L), (13a) h(t) ∈ J N [w(t)] a.e. t ∈ (0, L), (13b) w(0) = 0. ( 13c
)
We observe that ( 13) is equivalent to [START_REF] Eck | Unilateral contact problems : variational methods and existence theorems[END_REF], the verification is left to the reader. We regularize the problem (13), we construct a solution of the regularized problem, and we show the existence of a solution by passing to the limit with respect to regularity parameter.

We regularize the problem (13) by using Yosida's regularization. More precisely, we introduce first a sequence {g ε } ε>0 such that g ε belongs to C 0 ([0, L]) for all ε > 0 and

g ε → g in H as ε 0. ( 14
)
We may approximate ( 13) by the following problem:

w ε (t) + J N,ε [w ε (t)] = -g ε (t) a.e. t ∈ (0, L), (15a)
w ε (0) = 0. ( 15b
)
where

J N,ε def = 1 ε (1-J ε )
denotes the Yosida approximation with the resolvant defined by J ε def = (1+εJ N ) -1 . Note that J N,ε is a monotone and Lipschitz continuous mapping defined on all H. Clearly the contraction principle step by step in time allows us to prove that (13) possesses a unique solution w ε belonging to C 1 ([0, L]).

We establish now some a priori estimates, which later will enable to infer the existence of a weak solution. To this aim, we multiply (15a) by w ε (t), and we integrate this expression over (0, L) to get

L 0 |w ε (t)| 2 dt + L 0 J N,ε [w ε (t)]w ε (t) dt = - L 0 g ε (t)w ε (t) dt. ( 16
)
Let us define now

I ε (-∞,0] (w) def = min z∈H 1 2ε |w-z| 2 + I (-∞,0] (z)
for all ε > 0 and w ∈ H. Since J N is the subdifferential of a proper, convex and lower semicontinuous function

I (-∞,0] , it follows that I ε (-∞,0] is convex,
Fréchet differentiable in H and its subdifferential ∂I ε (-∞,0] coincides with J N,ε . Furthermore we have

∀z ∈ H : I ε (-∞,0] (z) I (-∞,0] (z) as ε 0.
Then we may deduce from ( 16) that

1 2 w ε 2 H + I ε (-∞,0] (w ε (L)) ≤ 1 2 g ε 2 H . ( 17
)
We observe that by definition I ε (-∞,0] (w ε (L)) ≥ 0 which follows from ( 17) that there exists C 1 > 0, independent of ε > 0, such that

w ε H ≤ C 1 . (18)
Furthermore we have

|w ε (τ )| 2 ≤ 2L w ε 2 H (19)
for all τ ∈ (0, L). Then we deduce from ( 17) and ( 19) that there exists

C 2 > 0, independent of ε > 0, such that |w ε (τ )| 2 + w ε 2 H ≤ C 2 (20)
for all τ ∈ (0, L). Since J ε is a contraction on all H (see [START_REF] Brézis | Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert[END_REF]), it follows from (20) that there exists C 3 > 0, independent of ε > 0, such that

w ε V + J ε w ε V ≤ C 3 . (21)
We observe that (15a) implies that

J N,ε [w ε ] 2 H ≤ 2 w ε 2 H + 2 g ε 2 H . (22)
The above a priori estimates allow us to infer that there exists w and h belonging to V ∩ L ∞ (0, L), respectively, and passing to subsequences, if necessary, we find

J ε w ε w in L ∞ (0, L) weak * , (23a) J N,ε [w ε ] h in H weak, (23b)
as ε tends to 0. Since V → H with compact embedding, it follows that

J ε w ε → w in H, (24a) w ε → w in H and w ε w in V weak, (24b)
as ε tends to 0. It remains to prove that w and h satisfy [START_REF] Hilber | Improved numerical dissipation for time integration algorithms in structural dynamics[END_REF]. We integrate (15a) over (0, τ ), we get

τ 0 w ε (t)+J N,ε [w ε (t)]+g ε (t) dt = 0 ( 25
)
for all τ ∈ (0, L). Thanks to ( 23) and ( 24), we may pass to the limit in all the terms of ( 25), namely we find

τ 0 w (t)+h(t)+g(t) dt = 0 for all τ ∈ (0, L), which is equivalent to (13a). Let us note that J N,ε [w ε ] ∈ J N [J ε w ε ],
for some technical details, the reader is referred to [START_REF] Brézis | Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert[END_REF]. Then (23b) and (24a) enable us to deduce that lim sup

ε 0 L 0 J N,ε [w ε (t)]J ε w ε (t) dt ≤ L 0 h(t)w(t) dt, (26)
which implies that (13b) holds. Then we may conclude that there exists a solution v to Problem ( 6)-( 8).

The uniqueness result comes from the monotonicity of J N , namely J N is the subdifferential of a convex, lower semi-continuous and proper function (see [START_REF] Brézis | Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert[END_REF] for further details).

The existence and uniqueness results to Problem (1)-(2) together with (4) is established by using analogous approach developed above for Problem (1)-( 3). The main difference comes from the choice of the auxiliary problem. Let ū be the unique solution to Problem (1) with initial data (2) and boundary conditions (4) where the Signorini boundary conditions are replaced by Dirichlet boundary conditions, i.e. ū(0, t) = 0, t ∈ [0, T ]. Hence the rest of the proof is quite a routine, the verification is left to the reader.

Finally, we consider the following energy which is associated to ( 5):

E(t) = 1 2 L 0 |u t (x, t)| 2 +|u x (x, t)| 2 dx.
This energy is constant with respect to time t when the density forces f vanishes.

Note that formulations (5) and ( 13) are equivalent (in the sense that they have the same weak solutions), since the proof given in [START_REF] Dabaghi | Convergence of mass redistribution method for the one-dimensional wave equation with a unilateral constraint at the boundary[END_REF] can be straightforwardly adapted.

Discretization

Classical finite element discretization

We introduce now a semi-discrete problem in space associated to [START_REF] Dabaghi | Convergence of mass redistribution method for the one-dimensional wave equation with a unilateral constraint at the boundary[END_REF] . To this aim, we choose a set of parameters h def = L m (mesh size) having in mind the limit m → +∞, where m is an integer, and let V Neu

h def = {v h ∈ C 0 ([0, L]) : v h | [a i ,a i+1 ] ∈ P 1 , i = 0, . . . , m -1} and V Dir h def = {v h ∈ V Neu h : v h (L) = 0}.
Here, a i def = ih, i = 0, 1, . . . , m and P 1 is the space of polynomials of degree less than or equal to 1. A classical basis of V h is given by the sequence of shape functions ϕ i ∈ V h for i = 0, 1, . . . , J, defined by

ϕ i (x) def = 1 -|x-a i | h if x ∈ [a max(i-1,0) , a min(i+1,m) ], 0 otherwise,
where J = m -1 in the case where V h = V Dir h and J = m in the case where

V h = V Neu h . Observe that ϕ i (a j ) = δ ij , j = 0, 1, . . . , m (δ is Kronecker's sym- bol).
We approximate the solution u belonging to V to the weak formulation (5) by

u h (x, t) = J j=0 u j (t)ϕ j (x).
Consequently, we have u i = u h (a i ), i = 0, 1, . . . , J. We shall consider two strategies to approximate the dynamic contact problem. The first one is a classical finite element semi-discretization which uses a multiplier. It reads as follows:

             find u h : [0, T ] → V h and λ : [0, T ] → R such that for all v h ∈ V h L 0 u h,tt v h dx + L 0 u h,x v h,x dx = -λv h (0) + L 0 f v h dx a.e. t ∈ [0, T ], 0 ≤ u h (0, •) ⊥ λ ≤ 0 a.e. t ∈ [0, T ], u h (•, 0) = u 0 h and u h,t (•, 0) = v 0 h , (27)
where u 0 h and v 0 h belong to V h and λ is the Lagrange multiplier being here the contact force. Furthermore, the approximated problem [START_REF] Pozzolini | Some energy conservative schemes for vibroimpacts of a beam on rigid obstacles[END_REF] can be rewritten in an algebraic formulation form as follows:

         find U h : [0, T ] → R n and λ : [0, T ] → R such that M U h,tt + SU h = -λe 0 + F a.e. t ∈ [0, T ], 0 ≤ u 0 ⊥ λ ≤ 0 a.e. t ∈ [0, T ], U h (0) = U 0 h and U h,t (0) = V 0 h , (28)
where M and S denote the mass and stiffness matrices, respectively, and

U def = (u 0 , . . . , u n-1 ) T , e 0 def = (1, 0, . . . , 0) T and F def = ( L 0 f ϕ 0 dx, . . . , L 0 f ϕ n-1 dx).

Finite element discretization with the mass redistribution method

An alternative to the standard discretization is to consider a trivial mass redistribution method which consists to replace the mass matrix M in ( 28) by a modified mass matrix M mod defined by

M mod ij def = L h ϕ i ϕ j dx for all i, j ∈ [0, n -1].
Let us observe that there are more sophisticated mass redistribution methods, see [START_REF] Hager | Analysis of a space-time discretization for dynamic elasticity problems based on mass-free surface elements[END_REF]. The modified mass matrix has the form

M mod def = 0 0 0 M ,
where M admits an unmodified part Mij = M i+1,j+1 , for all i, j = 1, . . . , n-2.

As a consequence and in a certain sense, the node at the contact boundary evolves in a quasi-static way. It has been established in [START_REF] Dabaghi | Convergence of mass redistribution method for the one-dimensional wave equation with a unilateral constraint at the boundary[END_REF] in a slightly different context that the approximated solution using the mass redistribution method converges to the unique solution to [START_REF] Dabaghi | Convergence of mass redistribution method for the one-dimensional wave equation with a unilateral constraint at the boundary[END_REF].

Finally, note that the discrete energy associated to problem ( 28) is given by

E h (t) = 1 2 (U T h,t M U h,t +U T h SU h -U T h,t F )(t). ( 29
)
We introduce now the time discretization. In order to fix the notations, let the time interval [0, T ] be divided by n + 1 discrete time-points such that 0 = t 0 < t 1 < . . . < t n = T . Furthermore the discrete quantities U n h , U n h,t , U n h,tt and λ n are assumed to be given by algorithmic approximations of the displacement U h (t n ), the velocity U h,t (t n ) and the acceleration U h,tt (t n ) and the Lagrange multiplier λ(t n ), respectively. Some time-stepping schemes allowing to obtain an approximated solution to Problem (1)-( 4) are introduced below and their efficiency is discussed and analyzed in the next sections.

The most wide-spread time-stepping scheme for solving contact problems is the family of classical Newmark methods proposed by Newmark in [START_REF] Newmark | A method of computational for structural dynamics[END_REF]. The algorithms form a subset of the Hilber-Hughes-Taylor (HHT) family of temporal integrators, the reader is referred to [START_REF] Hilber | Improved numerical dissipation for time integration algorithms in structural dynamics[END_REF] for further details. The underlying concept of the discretizations are Taylor expansions of displacements and velocities neglecting terms of higher order. For the contact problem introduced above, the discrete evolution is described by the finite difference equations:

U n+1 h = U n h + ∆tU n h,t + 1 2 -β ∆t 2 U n h,tt + β∆t 2 U n+1 h,tt , (30a)
U n+1 h,t = U n h,t + (1-γ)∆tU n h,tt + γ∆tU n+1 h,tt , (30b)
M U n+1 h,tt + SU n+1 h = -λ n+1 e 0 + F n+1 , (30c)
0 ≤ u n+1 0 ⊥ λ n+1 ≤ 0, ( 30d
)
where ∆t is a given time step and (β, γ) are the algorithmic parameters. Furthermore U 0 h , U 0 h,t and λ 0 are given and U 0 h,tt is evaluated by using (30c). Notice that the Newmark family contains many well-known and widely-used algorithms which correspond to different choices of the algorithmic parameters (β, γ) (see [START_REF] Hughes | The finite element method[END_REF][START_REF] Laursen | Computational contact and impact mechanics. Fundamentals of modeling interfacial phenomena in nonlinear finite element analysis[END_REF]). This method is unconditionally stable for linear elastodynamic problem for γ ≥ 1 2 and β ≥ 1 4 ( 1 2 +γ) 2 (see [START_REF] Hughes | The finite element method[END_REF][START_REF] Krenk | Energy conservation in Newmark based time integration algorithms[END_REF]). One of the most used method to study dynamic problems is the Crank-Nicolson method, also called trapezoidal or average acceleration method, which is obtained by setting the parameters (β, γ) = ( 14 , 1 2 ). This method is second-order consistent and unconditionally stable in the unconstrained case. Furthermore, the total energy of the discrete evolution is preserved, for the purely elastic case, see [START_REF] Laursen | Computational contact and impact mechanics. Fundamentals of modeling interfacial phenomena in nonlinear finite element analysis[END_REF]. That is the reason why the Crank-Nicolson method is commonly used in the community of computational mechanics. However, the situation should be examined in the case of contact constraints. Indeed the order of accuracy is degraded, for further details, the reader is referred to [START_REF] Hughes | The finite element method[END_REF][START_REF] Krenk | Energy conservation in Newmark based time integration algorithms[END_REF][START_REF] Grosu | Stability of semidiscrete formulations for elastodynamics at small time steps[END_REF].

Let us define now the energy evolution by ∆E

n h def = E n+1 h -E n h
, where E n h is assumed to be given by an algorithmic approximation of the energy E h (t n ). In particular, the energy evolution associated to [START_REF] Schatzman | A hyperbolic problem of second order with unilateral constraints: The vibrating string with a concave obstacle[END_REF]

for (β, γ) = ( 1 4 , 1 2 ) is given by ∆E n h = -(V n h ) T Λ n γ e 0 , (31)
where

V n h def = U n+1 h -U n h and Λ n γ def = (1-γ)λ n + γλ n+1 .
Let us observe that in the presence of permanent contact, the energy is strictly conserved, while the release of an existing contact decreases the energy and the detection of a new contact increases the energy. However, regain and loss of energy do not balance and there is no guarantee that the energy will stay bounded during the time integration.

Another approach to study elastodynamic contact problems consists to use the θ-method. More precisely, the discrete evolution is described by the finite difference equations:

U n+1 h = U n h + ∆t((1-θ)U n h,t +θU n+1 h,t ), (32a)
U n+1 h,t = U n h,t + ∆t((1-θ)U n h,tt +θU n+1 h,tt ), (32b)
M U n+1 h,tt + SU n+1 h = -λ n+1 e 0 + F n+1 , ( 32c
) 0 ≤ u n+1 0 ⊥ λ n+1 ≤ 0, (32d)
with θ ∈ [0, 1]. Furthermore U 0 h , U 0 h,t and λ 0 are given and U 0 h,tt is evaluated by using (32c). The energy evolution associated to [START_REF] Wriggers | Computational contact mechanics[END_REF] is given by

∆E n h = 1 2 -θ (V n h,t ) T M V n h,t + 1 2 -θ (V n h ) T SV n h -(V n h ) T Λ n θ e 0 . (33)
Notice that [START_REF] Wriggers | Computational contact mechanics[END_REF] with θ = 0 and θ = 1 are called the forward and backward Euler method, respectively, while with θ = 1 2 , we obtain Crank-Nicolson's method.

Finally, we focus on the so-called Paoli-Schatzman method that consists to fix the contact constraint at an intermediate time step. Indeed the method proposed below is a slight modification of Paoli-Schatzman method (see [START_REF] Paoli | Time discretization of vibro-impact[END_REF][START_REF] Paoli | A numerical scheme for impact problems. I ,II. The multidimensional case[END_REF]) which takes into account the kernel of the modified mass matrix. A simple application of Paoli-Schatzman method based on Newmark scheme to Problem [START_REF] Renard | The singular dynamic method for constrained second order hyperbolic equations: application to dynamic contact problems[END_REF]

with γ = 1 2 leads to                      find U n+1 h : [0, T ] → R n and λ n : [0, T ] → R such that: M (U n+1 h -2U n h +U n-1 h ) ∆t 2 + S(βU n+1 h +(1-2β)U n h +βU n-1 h ) = -λ n e 0 for all n ≥ 2, 0 ≤ u n,e 0 = u n+1 0 + eu n-1 0 1 + e ⊥ λ n ≤ 0, U 0 and U 1 given. ( 34
)
Here e belongs to [0, 1] and is aimed to be interpreted as a restitution coefficient. Note that U 0 h and U 1 h are given data and U 1 h can be evaluated by a one step scheme. We may observe that taking M = M mod in (34), we are not able to resolve the problem on the kernel of M mod . That is the reason why, SU n-1 h as well as SU n h are projected on the orthogonal of the kernel of M . The energy evolution associated to (34) is given for β = 1 4 by

∆E n h = 1+e 2 λ n u n-1 0 . ( 35
)
Notice that the numerical simulations were performed by employing the finite element library Getfem++ (see [START_REF] Renard | An Open Source generic C++ library for finite element methods[END_REF]). For the reader convenience, the energy evolutions ( 31), ( 33) and ( 35) are justified in the Appendix.

The wave equation with Signorini and Neumann boundary conditions

Analytical solution

We consider in this section the analytical solution proposed in [START_REF] Wriggers | Computational contact mechanics[END_REF]. It is the solution to the equation of motion (1) together with Cauchy initial data (2) and boundary conditions (4) in the case where the density of external forces f (x, t) = -(g+u x (0, t)) with g > 0 being the gravity acceleration. Let us describe the solution of this problem (see Figure 1). Before the impact, the bar is undeformed and its initial velocity v 0 vanishes. The bar reaches the rigid obstacle at time t 1 = 2u 0 g with the velocity 2u 0 g. After the impact, the bar stays in contact as long as the shock wave travels from bottom to the top of the bar and vice versa and then it takes off. Since the velocity of the shock wave c is assumed to be equal to 1, the bar stays in contact with the rigid obstacle a time L. Finally, notice that impacts occur at t 4k+1 = 3L + 16kL, t 4k+2 = t 4k+1 + 2L, t 4k+3 = t 4k+1 + 8L and t 4k+4 = t 4k+1 + 10L. We may deduce that the analytical solution is given by

u(x, t) =                  u 0 -1 2 gt 2 if t ≤ t 1 , h 1 (x, t-t 4k+1 ) if t 4k+1 < t ≤ t 4k+2 , h 2 (x, t-t 4k+2 ) if t 4k+2 < t ≤ t 4k+3 , h 1 (x, t 4k+4 -t) if t 4k+3 < t ≤ t 4k+4 , u 0 -1 2 g t-t 4k+4 -2u 0 g 2 if t 4k+4 < t ≤ t 4(k+1)+1 . (36)
x = 0 t = 0 L Figure 1: An elastic bar with both ends free to move where we have used the following notations

h 1 (x, t) def = - 2u 0 g min(x, L-|t-L|) + +∞ n=1 2g Lν 3 n (cos(ν n t)-1) sin(ν n x), h 2 (x, t) def = u 0 - 1 2 g t- 2u 0 g 2 - 2gL 2 3 + +∞ n=1 4g λ 2 n cos(λ n t) cos(λ n x), (37) with ν n def = (n-1 2 ) π L and λ n def = n π L .

Numerical simulations

The parameters used in the numerical simulations are L = 1, g = 0.11, u 0 (x) = 0.5, v 0 (x) = 0 and u x (L, •) = 0. We evaluate the approximated solution (U n h , λ n ) to Problem (1)-( 3) by using Crank-Nicolson method that consists to choose (β, γ) = ( 14 , 1 2 ) in [START_REF] Schatzman | A hyperbolic problem of second order with unilateral constraints: The vibrating string with a concave obstacle[END_REF] for the standard and modified mass matrices. Then we compare the solutions obtained by employing Crank-Nicolson method to the analytical one (u, λ) in different points of the bar (see Figure 2). In the case where the constraint is active, the numerical experiments highlighted some spurious oscillations occurring for the approximated solution (U n h , λ n ) obtained by solving [START_REF] Schatzman | A hyperbolic problem of second order with unilateral constraints: The vibrating string with a concave obstacle[END_REF] with the standard mass matrix M (see Figure 2 (left)) while these oscillations disappeared when the standard mass matrix is replaced by the modified mass matrix, i.e. M = M mod (see Figure 2 (right)). However the approximated solutions obtained by both methods converge to the analytical solution when both the space ∆x and time ∆t steps tend to 0 in L p (0, T ; H) norms with p = 2, +∞, see Figure 3 (a constant ratio ∆x ∆t = 1 has been used). Surprisingly, the error curves in the norms u h -u L p (0,T ;V) , with p = 2, +∞ shown in the same figure are both diverging for the standard and modified mass matrices. The explanation of this divergence is illustrated in Figure 4 where the space derivative of the displacement is shown for both the analytical and the approximated solution. Indeed, Figure 4 (left) illustrates the fact that the series of the derivative corresponding to (37) do not converge. This is why it seems not to be possible to use this analytical solution to estimate the error in L p (0, T ; V) norm. For that reason, another analytical solution is exhibited in the next section.

(a) u h -u L ∞ (0,T ;H) (b) u h -u L 2 (0,T ;H) (c) u h -u L ∞ (0,T ;V) (d) u h -u L 2 (0,T ;V)

The wave equation with Signorini and Dirichlet boundary conditions

Analytical solution

In this section, we propose a new analytical piecewise affine and periodic solution to the problem composed by the momentum equilibrium equation ( 1) together with Cauchy initial data (2) and boundary conditions (3) in the case where the density of external forces f (x, t) vanishes. Suppose that the initial displacement and velocity are given by u 0 (x) = 1 2 -x 2 and v 0 (x) = 0, respectively. The length of bar is L = 1, namely, the bar is compressed at t = 0. Let us describe the solution of this problem (see Figure 5). Before the impact, the bar, which is clamped at one end, elongates until it reaches the rigid obstacle at t 1 = 1. After the impact, the bar stays in contact for the time t ∈ (t 1 , t 2 ) and it takes off at time t 2 = 2. We use characteristics method to determine the solution of the above problem. To this aim, we split each domain defined by (0, L)×(t i , t i+1 ), i = 0, 1, 2, corresponding to the phases before the impact, during the impact and after the impact, respectively, into four regions according to Figure 6. We choose below L = 1 and t i = i, i = 0, . . . , 3. Before the impact, according to characteristics lines ξ def = x + t and η def = xt, we divide the domain (0, 1)×(0, 1) into four regions as it is drawn on Figure 6. Then we observe that ξ ∈ [0, 1] and η ∈ [0, 1] in the region I, ξ ∈ [0, 1] and η ∈ [-1, 0] in the region II, ξ ∈ [1, 2] and η ∈ [0, 1] in the region III and ξ ∈ [1, 2] and η ∈ [-1, 0] in the region IV. Therefore we may conclude that the solution is given by

x = 0 L t = 0 t 2 = 2 t 3 = 3 t 1 = 1
v = 0 u = 1-x 2 t 0 = 0 v = -1/2 u = 0 t 1 = 1 v = 1/2 u = 0 t 2 = 2 v = 0 u = 1-x 2 t 3 = 3
u(x, t) = 1 2 (1-x) in the regions I, III, (38a) u(x, t) = 1 2 (1-t) in the regions II, IV. ( 38b
)
Once again by using the characteristics lines, we divide the domain (1, 2) × (0, 1) into four regions. Then the solution (38b) calculated in the region IV allows us to evaluate the initial conditions at time t 1 = 1, namely we found u(•, 1) = 0 and u t (•, 1) = -1 2 . Therefore we obtain u(x, t) = -1 2 (t 1 -t) in the region I, (39a) u(x, t) = -x 2 in the region II, (39b) u(x, t) = 1 2 (x-1) in the region III, (39c) u(x, t) = 1 2 (t-t 1 -1) in the region IV. (39d)

After the impact, we divide the domain (2, 3) × (0, 1) into four regions. It follows from (39d) that u(•, 2) = 0 and u t (•, 2) = 1 2 which implies that u(x, t) = 1 2 (t-t 2 ) in the regions I, II, u(x, t) = 1 2 (1-x) in the regions III, IV. Furthermore we observe that u(•, 3) = u 0 and u t (•, 3) = v 0 . Then we conclude that the solution is periodic of period 3. Note that the Lagrange multiplier λ is equal to the contact force u x (0, •).

Numerical simulations

We evaluate below (U n h , λ n ) by using different time-stepping methods, namely the Newmark method (30), the backward Euler method (32) with θ = 1 and the Paoli-Schatzman method (34). More precisely, for each timestepping method, the approximated solutions (U n h , λ n ) are obtained by solving [START_REF] Schatzman | A hyperbolic problem of second order with unilateral constraints: The vibrating string with a concave obstacle[END_REF] for different values of the couple (β, γ), (34) for β = 1 4 and different values for the parameter e and (32) in the particular case where θ = 1, respectively. The computation are also performed for both the standard mass matrix M and the modified one M = M mod and they are compared to the analytical solution (u, λ) exhibited in Section 5.1 (see Figures 7,[START_REF] Hager | A stable energy-conserving approach for frictional contact problems based on quadrature formulas[END_REF][START_REF] Hilber | Improved numerical dissipation for time integration algorithms in structural dynamics[END_REF][START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF][START_REF] Krenk | Energy conservation in Newmark based time integration algorithms[END_REF]19 and 21). Furthermore, the energy evolution associated to different methods mentioned above are presented in Figures 8,[START_REF] Hager | Analysis of a space-time discretization for dynamic elasticity problems based on mass-free surface elements[END_REF][START_REF] Hughes | The finite element method[END_REF][START_REF] Kim | A boundary thin obstacle problem for a wave equation[END_REF][START_REF] Laursen | Computational contact and impact mechanics. Fundamentals of modeling interfacial phenomena in nonlinear finite element analysis[END_REF]20 and 22. Note that the numerical experiments were still done for different space and time steps with a constant ratio ∆x ∆t = 10. Some error curves for Crank-Nichoslon's method are presented for standard matrix as well as for modified mass matrix. Notice that the slope of error curves for nonlinear problems reflect both the regularity of solution and the ability of numerical schemes to reproduce some technical inequalities. In the Table 1, we used the following notation:

u h -u L p (H) = u h -u L p (0,T ;H) with p = 2, +∞.
Finally, the convergence rates for different time-stepping methods, analyzed in details above, are summarized in Tables 1 and2 for the reader convenience. These convergence rates were obtained for both standard (STD) and modified (MOD) mass matrices.

One can see in Table 1, that the rate of convergence of order one methods (Newmark I and backward Euler) in L p (0, T ; H) norms of the displacement (with p = 2, +∞) are limited to approximatively 1 2 for both the standard and modified mass matrices. Concerning the other methods, which are order two methods at least for the linear part of the problem, the rate of convergence is still approximatively 1 2 for the standard mass matrix but is improved to approximatively 1 for the modified mass matrix. Note that the Paoli-Schatzman scheme is not a fully order two method, except for e = 1, due to the discretization of the contact condition.

Scheme Method u h -u L ∞ (H) u h -u L 2 (H) u h -u L ∞ (V) u h -u L 2 (V) Crank-
Scheme Method λ h -λ L 2 (0,T ) Oscillations E h -E L ∞ (0,T ) E h -E L 2 (0,T ) Crank-Nicolson STD -0.
The advantage is less pronounced in L p (0, T ; V) norms of the displacement, but it is even more important in L 2 (0, T ) norm of the contact stress (see Table 2) since convergence does not occur for all the methods with the standard mass matrix.

An important remark is that, despite the low regularity of the exact solution, order two methods perform better than order one methods and should be preferred to approximate elastodynamic problem with impact.

Concerning the efficiency of the mass redistribution method, we can conclude that, even though it represents an additional approximation compared to the standard approximation, the convergence rates on the displacement and on the contact force are improved and the method reduces the potential spurious oscillations. The mass redistribution method appears then to be a robust method to approximate elastodynamic problem with impact.

The Crank-Nicolson method

We focus here on Newmark's method [START_REF] Schatzman | A hyperbolic problem of second order with unilateral constraints: The vibrating string with a concave obstacle[END_REF] in the particular case where (β, γ) = ( 14 , 1 2 ). The analytical solution (u, λ) exhibited in Section 5.1 and the approximated solutions (U n h , λ n ) obtained for different time and space steps are represented on Figure 7. We observe some spurious oscillations on Figure 7 (left) for the solution (U n h , λ n ) to Problem [START_REF] Schatzman | A hyperbolic problem of second order with unilateral constraints: The vibrating string with a concave obstacle[END_REF] after the contact takes place while these oscillations disappear when the standard mass matrix is replaced by the modified one, see Figure 7 (right).

It can easily be seen on Figure 8 that the scheme with the standard mass matrix is unstable with a rapidly growing energy while the one with the modified mass matrix is almost conservative as the space ∆x and time ∆t steps tend to 0. The error curves for (U n h , λ n ) highlighted that the norms u h -u L p (0,T ;H) , with p = 2, +∞, and λ n -λ H converge to 0 when M = M mod as n tends to +∞, which is unfortunately not the case for the standard mass matrix (see Figure 10). Finally, Figure 9 shows that the norm E h -E L p (0,T ) , p = 2, +∞, converges when M = M mod while this norm diverges when the standard mass matrix is considered.

a) E h -E L ∞ (0,T ) (b) E h -E L 2 (0,T ) (c) u h -u L ∞ (0,T ;H) (d) u h -u L 2 (0,T ;H) (e) u h -u L ∞ (0,T ;V) (f) u h -u L 2 (0,T ;V)

The Newmark method I

We deal now with Newmark's method [START_REF] Schatzman | A hyperbolic problem of second order with unilateral constraints: The vibrating string with a concave obstacle[END_REF] in the particular case where (β, γ) = ( 12 , 1). We observe some oscillations on Figure 11 (left) for the Lagrange multiplier λ n associated to Problem [START_REF] Schatzman | A hyperbolic problem of second order with unilateral constraints: The vibrating string with a concave obstacle[END_REF] after the contact takes place while these oscillations are much less preeminent when the standard mass matrix is replaced by the modified one, see Figure 11 (right). Furthermore, the present scheme is dissipative and stable.

The Newmark method II

We consider Newmark's method [START_REF] Schatzman | A hyperbolic problem of second order with unilateral constraints: The vibrating string with a concave obstacle[END_REF] in the case where (β, γ) = ( 1 2 , 1 2 ). Some oscillations on Figure 13 (left) for the Lagrange multiplier λ n associated to Problem (30) can be observed after the contact takes place while these oscillations are much less important when the modified mass matrix is considered in [START_REF] Schatzman | A hyperbolic problem of second order with unilateral constraints: The vibrating string with a concave obstacle[END_REF], see Figure 13 (right). Note that here the scales for the Lagrange multiplier associated to standard and modified mass matrices are different, see Figure 13.

The backward Euler method

We deal here with backward Euler's method (i.e. ( 32) with θ = 1). Since the present scheme is dissipative and stable (∆E n h < 0, see (33)), the approximated solutions (U n h , λ n ) obtained by using the standard mass as well as modified mass matrices do not oscillate as it can be observed in the previous numerical simulations.

The Paoli-Schatzman method I

We are interested here in Paoli-Schatzman's method with (β, e) = ( 1 4 , 0). Note that the stability result straightforwardly follows from [START_REF] Dumont | Vibrations of a beam between obstacles. Convergence of a fully discretized approximation[END_REF].

The Paoli-Schatzman method II

We focus here on Paoli-Schatzman's method with (β, e) = ( 1 4 , 1 2 ). Observe that the stability result is still an open problem for e = 0. Some spurious oscillations can be observed for the approximated Lagrange multiplier as well as for the approximated energy.

The Paoli-Schatzman method III

We consider now Paoli-Schatzman's method with (β, e) = ( 1 4 , 1). Note that this is a second order scheme. Acknowledgments. The support of the AVCR-CNRS "Mathematical and numerical analysis of contact problems for materials with memory" grant is gratefully acknowledged.

Appendix

The aim of the present section consists to establish the energy evolutions associated to schemes [START_REF] Schatzman | A hyperbolic problem of second order with unilateral constraints: The vibrating string with a concave obstacle[END_REF], [START_REF] Wriggers | Computational contact mechanics[END_REF] and (34). Let us observe that the energy evolution, when the external force f = 0, is given by

∆E n h = 1 2 (U n+1 h,t -U n h,t ) T M (U n+1 h,t +U n h,t ) + 1 2 (U n+1 h -U n h ) T S(U n+1 h +U n h ). ( 40
)
We prove now that the energy evolution associated to the Crank-Nicolson method is given by [START_REF] Tkachuk | Hybrid-mixed discretization of elasto-dynamic contact problems using consistent singular mass matrices[END_REF]. To this aim, we observe that (30a) and (30b) lead to Clearly we have ∆t((1-θ)(-SU n h -λ n e 0 ) + θ(-SU n+1 h -λ n+1 e 0 )) = M (U n+1 h,t -U n h,t ), which allows us to conclude that (33) holds.

M (U n+1 h -U n h -∆tU n h,t ) = -∆t 2 1 2 -β SU n h +βSU n+1 h -∆t 2 1 2 -β λ n e 0 +βλ n+1 e 0 , (41) and
Finally we establish the energy evolution associated to (34) in the particular case where β = 1 4 . To this aim, let us introduce the energy evolution as follows: ).

∆E n h = 1 2 U T h,t M U h,t + 1 2 U T h SU h n+ 1 2 n-
Then it comes that

U n+ 1 2 h +U n-1 2 h = 1 2 (U n+1 h +2U n h +U n-1 h
) and U

n+ 1 2 h -U n-1 2 h = 1 2 (U n+1 h -U n-1 h ),
and also h,t , so we have

U n+1 h = U n h + ∆tU
U n+ 1 2 h,t -U n-1 2 h,t = U n+1 h -2U n h +U n-1 h ∆t and U n+ 1 2 h,t + U n-1 2 h,t = U n+1 h -U n-1 h ∆t .
Obviously, we get

∆E n h = 1 2 (U n+1 h -U n-1 h ) T M U n+1 h -2U n h +U n-1 h ∆t 2 + S U n+1 h +2U n h +U n-1 h 4 .
According to (34) together with β = 1 4 , we deduce that

∆E n h = 1 2 (U n+1 h -U n-1 h ) T (-λ n+1 e 0 ) = 1

Figure 2 :

 2 Figure 2: First analytical solution: comparison of the analytical solution (u, λ) and approximated solutions (u n h , λ n ) by using standard mass matrix (left) and modified mass matrix (right) in the contact with Crank-Nicolson method.

Figure 3 :

 3 Figure 3: First analytical solution: comparison of the error curves obtained by using standard and modified mass matrices with Crank-Nicolson method. A constant ratio ∆x ∆t = 1 has been used.

Figure 4 :

 4 Figure 4: First analytical solution: representation of the derivative of displacement with respect to x of analytical (left) and approximated (right) solutions. The analytical solution is computed with the 10 3 first terms of the series.

Figure 5 :

 5 Figure 5: An elastic bar with one end fixed and the other end free to move

Figure 6 :

 6 Figure 6: The regions allowing to determine the value of u.

Figure 7 :

 7 Figure 7: Comparison of the analytical solution (u, λ) and the approximated solutions (U n h , λ n ) by using the standard (left) and modified (right) mass matrices in the contact node with Crank-Nicolson method.

Figure 8 :

 8 Figure 8: Comparison of the energy associated to the analytical solution and the energy associated to the approximated ones for the standard (left) and modified (right) mass matrices with Crank-Nicolson method.

(

Figure 9 :Figure 10 :

 910 Figure 9: Comparison of the error curves of the energy obtained by using the standard and modified mass matrices with Crank-Nicolson method.

Figure 11 :

 11 Figure 11: Comparison of the analytical solution (u, λ) and the approximated solutions (U n h , λ n ) by using the standard (left) and modified (right) mass matrices in the contact with Newmark method I.

Figure 12 :

 12 Figure 12: Comparison of the energy associated to the analytical solution and the energy associated to the approximated ones for standard (left) and modified (right) mass matrices with Newmark method I.

Figure 13 :

 13 Figure 13: Comparison of the analytical solution (u, λ) and the approximated solutions (U n h , λ n ) by using the standard (left) and modified mass matrices (right) in the contact with Newmark method II.

Figure 14 :

 14 Figure 14: Comparison of the energy associated to the analytical solution and the energy associated to the approximated ones for the standard (left) and modified (right) mass matrices with Newmark method II.

Figure 15 :

 15 Figure 15: Comparison of the analytical solution (u, λ) and the approximated solutions (U n h , λ n ) by using the standard (left) and modified (right) mass matrices in the contact with backward Euler method.

Figure 16 :

 16 Figure 16: Comparison of the energy associated to the analytical solution and the energy associated to the approximated ones for the standard (left) and modified (right) mass matrices with backward Euler method.

Figure 17 :

 17 Figure 17: Comparison of the analytical solution (u, λ) and the approximated solutions (U n h , λ n ) by using the standard (left) and modified (right) mass matrices in the contact with Paoli-Schatzman method I.

Figure 18 :

 18 Figure 18: Comparison of the energy associated to the analytical solution and the energy associated to the approximated ones for the standard (left) and modified (right) mass matrices with Paoli-Schatzman method I.

Figure 19 :

 19 Figure 19: Comparison of the analytical solution (u, λ) and the approximated solutions (U n h , λ n ) by using the standard (left) and modified (right) mass matrices in the contact with Paoli-Schatzman method II.

Figure 20 :

 20 Figure 20: Comparison of the energy associated to the analytical solution and the energy associated to the approximated ones for the standard (left) and modified (right) mass matrices with Paoli-Schatzman method II.

Figure 21 :

 21 Figure 21: Comparison of the analytical solution (u, λ) and the approximated solutions (U n h , λ n ) by using standard (left) and modified (right) mass matrices in the contact with Paoli-Schatzman method III.

Figure 22 :

 22 Figure 22: Comparison of the energy associated to the analytical solution and the energy associated to the approximated ones for the standard (left) and modified (right) mass matrices with Paoli-Schatzman method III.

Table 1 :

 1 Convergence rates for the displacement

	Nicolson	STD MOD	0.48045 0.88075	0.47113 0.97113	-0.30273 0.38624	-0.31299 0.36192
	Newmark I	STD	0.51444	0.51768	0.26104	0.24000
	(β, γ) = ( 1 2 , 1)	MOD	0.49236	0.48034	0.28219	0.25943
	Newmark II	STD	0.57503	0.61856	0.34411	0.30100
	(β, γ) = ( 1 2 , 1 2 )	MOD	0.95055	1.00090	0.39734	0.33789
	Backward Euler	STD MOD	0.51804 0.50425	0.51771 0.49024	0.26784 0.27753	0.24847 0.25465
	Paoli-Schatzman I	STD	0.66064	0.66015	0.31204	0.28406
	e = 0	MOD	0.98009	1.00690	0.39570	0.41096
	Paoli-Schatzman II	STD	0.66064	0.66014	0.33914	0.30116
	e = 1 2	MOD	0.97659	1.00170	0.41159	0.42823
	Paoli-Schatzman III	STD	0.58876	0.61110	0.39286	0.35512
	e = 1	MOD	0.97658	1.00200	0.46380	0.46993

Table 2 :

 2 Convergence rates for the Lagrange multiplier and energy

  n+1 e 0 -λ n e 0 ) . T (1γ)λ n e 0 + γλ n+1 e 0 ) T (λ n+1 e 0λ n e 0 ). By setting the parameters (β, γ) = (14 , 1 2 ) (31) holds. Let us establish the energy evolution associated to the θ-method. We infer from (32a) and (32c) thatM (U n+1 h,t -U n h,t ) = (1-θ)∆t(-SU n h -λ n e 0 ) + θ∆t(-SU n+1 h -λ n+1 e 0 ).(44)

	Since M is a symmetric matrix, it follows from (42) that
	∆E =(	1 2	-γ)(U n+1 h	-U n h ) T S(U n+1 h	-U n h )
	-(U n+1 h h ) +∆t(β --U n γ 2 )(U n+1 h,t -U n h,t ) T S(U n+1 h +∆t(β -γ 2 )(U n+1 h,t -U n h,t Then using (44) in (40), we find	h ) -U n
	∆E n h = -+ 1 2 (U n+1 ∆t 2 h -U n (U n+1 h,t +U n h,t ) T ((1-θ)SU n h +θSU n+1 h h ) T S(U n+1 h +U n h ).	+((1-θ)λ n +θλ n+1 )e 0 )
	M (U n+1 h,t -U n h,t ) = -∆t((1-γ)SU n h +γSU n+1 Then we multiply (42) by ∆t 2 and we subtract this expression to (41), we get (42) Note that from (32a) ∆t 2 (U n+1 h,t +U n h,t ) = U n+1 h -U n h + ∆t 1 2 -θ (U n+1 h,t ), h,t -U n it follows that
	M (U n+1 h,t +U n h,t ) = ∆E n h = 1 2 -θ (U n+1 h -U n 2 ∆t h ) T S(U n+1 M (U n+1 h -U n h -U n h ) -(U n+1 h h ) + 2∆t β-γ 2 S(U n+1 h -U n h ) + 2∆t β-γ -U n h ) T (((1-θ)λ n +θλ n+1 )e 0 ) 2 (λ (43) + ∆t 1 2 -θ (U n+1 h,t -U n h,t )
	By using (43) in (40), it follows that
	∆E n h =
	1 ∆t + ∆t β-(U n+1 h,t -U n h,t ) T M (U n+1 h γ 2 (U n+1 h,t -U n h,t ) T λ n+1 e 0 -λ n e 0 ) + -U n h ) + ∆t β-γ 2	(U n+1 h,t -U n h,t ) T S(U n+1 h 1 2 (U n+1 h -U n h ) T S(U n+1 -U n h ) h +U n h ).

h

) -∆t((1-γ)λ n e 0 +γλ n+1 e 0 ). T ((1-θ)(-SU n h -λ n e 0 ) + θ(-SU n+1 h -λ n+1 e 0 )).


Yanan Cao

Cungen Cao

Liangjun Zang

Shi Wang

Dongsheng Wang

Extracting Comparative Commonsense from the Web

Commonsense acquisition is one of the most important and challenging topics in Artificial Intelligence. Comparative commonsense, such as "In general, a man is stronger than a woman", denotes that one entity has a property or quality greater or less in extent than that of another. This paper presents an automatic method for acquiring comparative commonsense from the World Wide Web. We firstly extract potential comparative statements from related texts based on multiple lexico-syntactic patterns. Then, we assess the candidates using Web-scale statistical features. To evaluate this approach, we use three measures: coverage of the web corpora, precision and recall which achieved 79.2%, 76.4% and 83.3%, respectively in our experiments. And the experimental results show that this approach profits significantly when the semantic similarity relationships are involved in the commonsense assessment.

Introduction

Commonsense knowledge plays an important role in various areas such as natural language understanding (NLU), information retrieval (IR), question answering (QA), etc. For decades, there has been a thirst in artificial intelligence research community for large-scale commonsense knowledge bases. Since the hand-coded knowledge base suffers from semantic gaps and noises, and the human effort should still be involved in the maintenance [START_REF] Singh | The Public Acquisition of Commonsense Knowledge[END_REF], much recent work focuses on automatic strategies to acquire commonsense knowledge.

As information people use every day, commonsense knowledge encompass many aspects of life. It relates to an individual object, phenomena and activity, or relations between different concepts and a sequence of events. This paper concentrates on a commonsensical concept-level relation, which is defined as follows.

Comparative Relation We say that concept c 1 and c 2 are comparative if and only if the natural language description "In general, c 1 is relatively * than c 2 " is reasonable and acceptable, where "*" is an adjective expressing the comparative degree. This statement denotes that the concept c 1 has a property or quality greater or less in extent than that of c 2 . We use the triple (c 1 , c 2 , p) to express it formally. For a specific adjective in this statement, we call it a comparative property of the given pair of concepts. For example, "In general, a man is relatively stronger than a woman" describes a comparative relation between the concepts man and woman, of which stronger is the comparative property. The formalization of this statement is (man, woman, stronger ).

In this paper, we aim to acquire comparative commonsense about given pairs of concepts from the Web, which is a vast source of information. Firstly, we submit instantiated lexico-syntactic patterns as query terms to the search engine, such as "men are * than women", and extract potential comparative statements from the retrieved texts. Then, we assess the candidates based on multiple Webscale statistical features. The semantic similarity relationships, including synonym and antonym, are involved in the computing of these statistics. This improves the effectiveness of our approach. In the experiments, we selected twenty pairs of concepts as the input data and the result set of comparative statements had a precision of 76.4% and a recall of 83.3%. Besides, we show the high coverage of the comparative information in a test web corpus, which achieved 79.2%.

The remainder of this paper is organized as follows. In Section 2, we review some related work on commonsense knowledge acquisition from texts. Section 3 and Section 4 describe the two-stage method for extracting and assessing comparative commonsense in detail. In Section 5, we evaluate our approaches and discuss the experimental results. Finally, we conclude in Section 6.

Related Work

To acquire commonsense knowledge, a number of efforts tap on textual data sources. The pioneering work mentioned in [START_REF] Schubert | Can We Derive General World Knowledge from Texts[END_REF] recognized text as a source of both explicit and implicit commonsense knowledge. Capitalizing on interpretive rules, they obtained a broad range of general propositions from noun phrases and clauses. And they tried to derive stronger generalizations based on the statistical distribution of the obtained claims. However, to design abundant rules is indeed time-cost and the extraction process was not implemented automatically.

[3] extracted commonsensical inference rules from coordinated verb phrases in Japanese text. This work was based on the assumption that if two events shared a common participant, then the two events would probably be a probablyfollow relation. They used an unsupervised approach to select the inference rules based on lexical co-occurrence information. The results were evaluated by multisubjects, and just 45% of the selected rules were accepted by 80% of the human subjects. The low precision showed that co-occurrence is important but not adequate to assess the specific causal relation between two events, because there are other event-level relations such as IsPreferableTo in [START_REF] Popescu | Information Extraction from Unstructured Web Text[END_REF].

[4] used the Web data to clean up and augment existent commonsense knowledge bases. For specific binary event-level relations such as causality, they used complete or partial event descriptions and lexico-syntactic patterns to mine potential instances of the relations, and assessed these candidates based on Web-scale statistics. For higher-level predicates, the relational definitions were involved in the assessment, which extended the work in [START_REF] Torisawa | An Unsupervised Learning Method for Commonsensical Inference Rules on Events[END_REF]. The advantage of this work is that some implicit instances of relations can be assessed at the level of the Web rather than found in one particular sentence.

In our recent work [START_REF] Cao | Acquiring Commonsense Knowledge about Properties of Concepts from Text[END_REF], we acquired commonsense knowledge about properties of concepts by analyzing how adjectives are used with nouns in everyday language. We firstly mined a large scale corpus for potential concept-property pairs using lexico-syntactic patterns and then filtered noises based on heuristic rules and statistical approaches. For each concept, we automatically selected the commonsensical properties and evaluated their applicability. The precision of the result-set achieved 67.7%. In the following, we continue our research on the relations between concepts and properties, from the view of comparing.

Extracting Candidate Comparative Statements

This section describes the first phase of acquiring comparative commonsense: we extract candidate commonsense from the Web using linguistically-motivated patterns [START_REF] Hearst | Automatic acquisition of hyponyms from large text corpora[END_REF]. In the text, comparative relations are expressed in various forms. Table1 lists several patterns followed by corresponding instances, in which NP denotes Noun Phrase.

Table 1. Examples of patterns designed to extract comparative commonsense

Pattern Instance

N P1 be * than N P2 men are stronger than women N P1 be more * than N P2 men are more powerful than women N P2 be not * than N P1

women are not stronger than men N P2 be less * than N P1

women are less powerful than men to compare N P1 and N P2, N P1 be * to compare men and women, men are stronger compared with N P2, N P1 be * compared with women, men are stronger

When we perform the extraction process, these patterns are automatically instantiated with given pairs of concepts to generate query terms. For example, to acquire comparative statements about man and woman, we instantiate the first pattern as "men are * than women" and the last one as "compared with women, men are *". Issuing the query terms, we take advantage of the search engine to retrieve relevant texts from the Web and obtain the hit count of each query string. And subsequently, we extract a comparative property from each snippet matched with the "*" wildcard, which constitutes a candidate triple with its corresponding concept-pair.

There are two explicit constrains on the extraction. The first is that the matching snippet should be an adjective phrase, and we just extract the head word no matter it has a modifier or a complement. For example, from the text "men are much stronger than women", we extracted the comparative property "stronger", while "much" is discarded. We implement this constraint by analyzing the combination of POS in adjective phrases. Second, we extracted the comparative statements in a context-independent manner. That means we don't extract the specific context even if a statement is not acceptable anymore without it. For example, in the sentence "From Shanghai to Lhasa, a train is faster than an airplane", the extracted statement (train, airplane, faster ) is dependent on the context "from Shanghai to Lhasa", or else it violates the general knowledge. To avoid obvious errors induced by this strategy, we restrict the matching snippet to be a single sentence or to have a common modifier such as "generally speaking", "as we all know", etc.

Assessing Candidate Comparative Statements

During the extraction process, we acquire all potential comparative statements. However, the Web contains massive unreliable or biased information [START_REF] Kosala | Web mining research: a survey[END_REF], and some knowledge we obtained is inconsistent with the commonsensible facts. So, we verify each candidate statement based on multiple Web-scale statistics.

Statistical Features for commonsense assessment

Feature 1 Occurrence Frequency. Because a frequently mentioned piece of information is typically more likely to be correct than an item mentioned once or twice, high frequency is an important feature for truth. Given a candidate statement (c 1 , c 2 , p), we could easily obtain its occurrence frequency F req(c 1 , c 2 , p) in the web corpora using hit counts of instantiated queries, which are returned by the search engine:

F req(c 1 , c 2 , p) = pt∈P T hits(pt(c 1 , c 2 , p)), (1)
where PT is the set of predefined patterns (referred in Table 1) and hits(pt(c 1 , c 2 , p)) represents the hit count for a query which is an instance of the pattern pt.

Feature 2 Confidence. Although frequency is an important feature, it's not a guarantee of truth. Given a pair of concepts c 1 and c 2 , we use confidence to indicate the probability of a property p to be their comparative property. More specifically:

Conf (c 1 , c 2 , p) = F req(c 1 , c 2 , p)/ F req(c 1 , c 2 , p) (2)
In this formula, the denominator is the number of potential comparative relations which c 1 and c 2 participant in, and "*" is any potential comparative property of c 1 and c 2 .

Feature 3 the Number of the Matched Patterns. The last feature is the number of different expression forms of the same comparative statement. It's motivated by the assumption that if a potential statement is extracted by multi-patterns, it seems more reliable.

Semantic Similarity Relationships in the Candidate Set

We note that, the candidate statements are not independent, and there are semantic similarity relationships among them. For example, (woman, man, weaker ) is consistent with the statement (man, woman, stronger ), while (man, woman, weaker ) is opposed to it. We make use of these semantic similarity relationships in the computing of the statistics. To see this, we first begin by introducing several fundamental notations. Given a potential comparative statement (c 1 , c 2 , p), we use SY N p = {x|x is a synonym of p} ∪ {p} to denote the set of all synonyms of the property p, and another set AN T p = {y|y is a synonym of p}, which contains all antonyms of p. Based on these two sets, we divide other statements into the following categories according to the relationship between it and (c 1 , c 2 , p).

• A supporting example, which has the same meaning of the statement (c 1 , c 2 , p). For instance, (man, woman, more powerful ) and (woman, man, weaker ) are both supporting examples of (man, woman, stronger ). It belongs to the set

SupSet(c 1 , c 2 , p) = {(c 1 , c 2 , p )|p ∈ SY N p } ∪ {(c 2 , c 1 , p |p ∈ AN T p }
• A counterexample, which is opposed to the statement (c 1 , c 2 , p) in semantic. Using above instance (man, woman, stronger ), (man, woman, weaker ) and (woman, man, stronger ) are both its counterexamples. The set of the counterexamples is

CntSet(c 1 , c 2 , p) = {(c 1 , c 2 , p )|p ∈ AN T p } ∪ {(c 2 , c 1 , p |p ∈ SY N p }
• An irrelative example is the statement that is neither a supporting example nor a counterexample. For example, (man, woman, braver ) and (man, woman, stronger ) are irrelative to each other.

Then, we use Support(c 1 , c 2 , p) to denote the frequency of supporting examples of the given triple (c 1 , c 2 , p). More specifically:

Support (c 1 , c 2 , p) = (c1,c2,p )∈SupSet(c1,c2,p) F req (c 1 , c 2 , p ) (3)
And we use Counter(c 1 , c 2 , p) to denote the frequency of counterexamples of the given triple (c 1 , c 2 , p). It is computed as follows:

Counter (c 1 , c 2 , p) = (c1,c2,p )∈CntSet(c1,c2,p) F req (c 1 , c 2 , p ) (4)
Intuitively, if supporting examples of the statement (c 1 , c 2 , p) are much more than its counterexamples, this statement is more likely to be true. So, we use ExtF req(c 1 , c 2 , p) instead of F req(c 1 , c 2 , p) to assess potential comparative statements:

ExtF req (c 1 , c 2 , p) = Support(c 1 , c 2 , p) -Counter(c 1 , c 2 , p) (5)
And confidence is accordingly computed as follows:

ExtConf (c 1 , c 2 , p) = ExtF req(c 1 , c 2 , p)/ F req(c 1 , c 2 , * ) (6)
Although the features frequency and confidence are not novel, we combine Web-scale statistics and semantic relationships, and this improves the performance of our approach.

Experiments

In this section, we report some experimental results on comparative commonsense acquisition from Chinese web corpora. We computed the estimated coverage to prove the Web a proper resource for comparative commonsense acquisition. And we show that using similarity relationships in the assessment improves the precision and recall of our approach.

Experimental Settings

We selected twenty pairs of our familiar nouns as the test concept-set, which were averagely divided into four classes: human, plant, animal (except for human) and artificiality. These concepts were all basic level categories [START_REF] Zhang | Cognitive Linguistics and Chinese Noun Phrases[END_REF] (such as "plane"), which are neither too abstract (such as "vehicle") nor too concrete (such as "jet plane"). This is because the concepts on basic level are more likely to be proper for commonsense acquisition [9]. Intuitively, not arbitrary pairs of concepts are comparable. They always have properties in common, but one's property is greater or less in extent than that of the other. To satisfy this assumption, we picked out a concept-pair from the same class, such as "man" and "woman", "watermelon" and "apple", rather than "man" and "apple". Next, we issue each concept-pair to the Google Search Engine and downloaded relevant texts the search engine returned from the Web. We call it the Test Corpus in the following experiments.

Results: Coverage of the Web Resource

Because an average person is assumed to possess commonsense knowledge, it is not communicated most of the time. That means commonsense knowledge is usually used implicit in the texts. However, it seems infeasible to extract knowledge implicit in texts just using the pattern-matching method. To investigate what amounts of comparative knowledge we might extract from explicit expressions on the Web, we propose the coverage of comparative information on the web:

Coverage = the number of comparative statements explicit on the W eb the number of comparative statements

In fact, we don't know the accurate quantity of commonsense possessed by people or that of information distributed on the whole web. So we estimated these values based on several human subjects and the Test Corpus. More specifically: EstCoverage = the number of comparative statements in the T est Corpus the number of comparative statements f rom human subjects

We asked five human subjects to contribute their comparative commonsense on the test concept-set. The subjects' knowledge was gathered in an iterative way. Firstly, they were limited to derive this knowledge just from their brains. After the subjects submitted their outcome for the first time, we presented them the Test Corpus. In the second phase, the subjects manually extracted comparative commonsense from the corpus, which may provide them new clues. Meanwhile, the subjects were allowed to modify their primary submission. Table 2 shows the quantity of the twice submissions of the subjects and knowledge automatically mined from the web corpus, respectively. We note that, the subjects supplemented 46 records according to the web resource, which were not referred in the first submission. And the average coverage achieved 79.2%, which demonstrated the Web a good resource containing adequate explicit comparative commonsense.

Results: Effectiveness of Our Approach

In this experiment, we used the feature Freq and Conf to assess potential statements as the baseline method, and used ExtFreq and ExtConf as the extended method. Precision and recall are used to evaluate the performance of these approaches.

Before our extraction, we pre-processed the related texts in the test corpus including HTML label analyzing, word segmentation and POS tagging [10]. We extracted 429 potential comparative statements using the heuristic rules (referred in Section 3.1). Then, we expressed these acquired triples in the natural language "In general, c 1 is relatively * than c 2 ", and asked the five subjects to judge whether each description was acceptable. If a potential statement was accepted by 60% of the human subjects, we consider it a piece of comparative commonsense. Finally, we assessed each potential statement using the statistical features in both baseline method and extended method, and the results were compared with that from the human judgments. shows the performance of both the baseline method and extended method. Given the minimum threshold of confidence and the number of matched patterns which are 0.057 and 3 respectively, we ranked the comparative statements according to their frequencies. This graph plots the average precision and the top n acceptable statements ranked for a pair of concepts. It's obvious that the extended method outperformed the baseline scores. The average precision and recall achieved 76.4% and 83.3%, respectively.

Comparative commonsense describes the concept-level relations from the view of comparing, which denotes that one entity has a property or quality greater or less in extent than that of another. In this paper, we propose an automatic method for acquiring comparative commonsense from the World Wide Web. We use multiple lexico-syntactic patterns to extract potential comparative statement from related texts retrieved by the search engine. And then, we assess the candidates by combining Web-scale statistics and their semantic similarity relationships including synonym and antonym.

In the experiments, we evaluated the coverage of the explicit comparative information in the web corpus, which demonstrated the Web a good resource for our extraction. The experimental results also showed that the use of semantic similarity relationships in assessment significantly improves precision and recall of our approach when the statements have close dependency.

This work is based on given comparative pairs of concepts. In the future work, we will research on the characteristics of the comparative concepts, and automatically select them as the input data.

Fig. 1 .

 1 Fig. 1. Comparison of two methods

Table 2 .

 2 Quantity of comparative commonsense obtained from 5 human subjects and that in the test web corpus

	Noun Class	Human Subjects primary modified	Test Corpus Coverage
	Human	99	128	102	79.7%
	Plant	27	33	30	90.9%
	Animal	57	57	31	54.4%
	artificiality	50	61	58	95.1%
	total	233	279	221	79.2%

Acknowledgements. This work is supported by the National Natural Science Foundation of China under Grant No. 60773059.

This article was processed using the L A T E X macro package with LLNCS style


Guang Jiang

Cao Cungen

Sui Yuefei

Han Lu

Shi Wang

A General Approach to Extracting Full Names and Abbreviations for Chinese Entities from the Web

Identifying Full names/abbreviations for entities is a challenging problem in many applications, e.g. question answering and information retrieval. In this paper, we propose a general extraction method of extracting full names/abbreviations from Chinese Web corpora. For a given entity, we construct forward and backward query items and commit them to a search engine (e.g. Google), and utilize search results to extract full names and abbreviations for the entity. To verify the results, filtering and marking methods are used to sort all the results. Experiments show that our method achieves precision of 84.7% for abbreviations, and 77.0% for full names.

Introduction

Named Entity Recognition (NER) is a basic task in text mining; it is significant for information extraction, machine translation etc. in nature language processing (NLP). In 1998, MUC (Message Understanding Conference) defined seven categories of named entity task belong to three subtasks: entity names (organizations, persons, locations), times (dates, times), and quantities (monetary values, percentages). Named Entity Recognition is a challenging topic in NLP because named entities are huge and increase as time goes on.

Among the seven categories above, organizations, persons and locations are three most import types, therefore identification of these three categories become hot points of research [START_REF] Maynard | Identifying Terms by Their Family and Friends[END_REF][START_REF] Luo | Two-Character Chinese Word Extraction Based on Hybrid of Internal and Contextual Measures[END_REF][START_REF] Wang | Learning Concepts from Text Based on the Inner-Constructive Model[END_REF]. In view of full names and abbreviations, there are two kinds of entities: full name entity and abbreviation entity. Identification of full name and abbreviation for person names and locations are easier than which of organizations, because they embody more regular laws: person name begins with a family name, e.g. (Manager Zhang), while location names are finite and their abbreviations usually are fixed usually. E.g., (Jing) is an abbreviation of (Beijing). However, organizations are more difficult to identify because they are often long, loose and flexible or sometimes contain locations.

Identification of named entities is prepared for extracting relation of different entities, which is significant in several fields such as question answering (QA), multi-word expression (MWE) and information retrieval (IR). Automatic content extraction (ACE)3 evaluation plan organized by the National institute of standards and technology (NIST) had defined seven types of entity relation: ART (artifact), GEN-AFF (Gen-affiliation), METONYMY, ORG-AFF (Orgaffiliation), PART-WHOLE, PRE-SOC (person-social) and PHYS (physical). They also can be divided into much more subtypes.

Related Research

Many researchers had been focused on the identification of named entity relations [START_REF] Li | Automated Entity Relation Tuple Extraction Using Web Mining[END_REF][START_REF] Liang | Information Extraction System Based on Event Frame[END_REF]. A general route is to convert the problem into classification, which means relations of entities are categorized as several types and a classifier is constructed. Different machine learning algorithms are applied to solve this question [START_REF] Che | Automatic Entity Relation Extraction[END_REF][START_REF] Miller | A Novel Use of Statistical Parsing to Extract Information from Text[END_REF]Kambhalta et al., 2004). As to Chinese entity relation, [START_REF] Che | Automatic Entity Relation Extraction[END_REF] attempted to extract relations of entities automatically using Winnow and SVM, with F-score 73.08% and 73.27% respectively. [START_REF] Liu | Implementation of a Kernel-Based Chinese Relation Extraction System[END_REF] used HowNet to acquire semantic knowledge, and combine both semantic sequence kernel function and KNN classifier to extract relation, the accuracy achieve about 88%. [START_REF] Dong | Chinese Automatic Entity Relation Extraction[END_REF] divided entity relation into categories: embedding relations and non-embedding relations, and construct different syntactic features respectively. Their experiment proved the method can improve performance of Chinese entity relation extraction task.

In this paper, we investigate to extract Chinese full names and abbreviations for given entities. Entities discussed in this paper are not limited to named entities, more precisely, are lexical nominal concepts, which are discriminative and exist independently; they can be either physical (e.g. Peking University) or abstract (e.g. scientific technique). In the following, we denote an entity with a Chinese full name as full name entity, and that with an abbreviated Chinese name as abbreviation entity. Different from other foreign languages, Chinese words and Chinese characters are different, we will refer to word as Chinese word and Chinese character as character itself.

In our paper, we don't study document-specific abbreviation, which means the abbreviation is just used in a specific document or context. E.g. In sentence . . . , the abbreviation (Wuhan subway) appears just as a concise anaphora in the document, but not a common fixed abbreviation.

3 Extraction of full name entity and abbreviation entity from the Web

Comparison & analysis of full name entity and abbreviation entity

There are many common architectural features for a full name entity and its corresponding abbreviation entities. We sum up several points as follows:

1. Full name entity is consisted of several words, and its abbreviation entity is composed by one Chinese character from each word. E.g. (Women's Federation) and (Fulian). 2. Full name entity is consisted of several morphemes, delete some of them and constitute its abbreviation entity. E.g. (Tsinghua University) and (Tsinghua). 3. Full name entity is consisted of coordinating morphemes; first Chinese character of each morpheme and common morphemes constitute its abbreviation entity. E.g. (agriculture and industry) and (Gongnongye). 4. Full name entity is consisted of coordinating morphemes; numeral and the common morphemes constitute its abbreviation entity. E.g. (modernization of agriculture, industry, national defense and science and technology) and (the four modernizations). 5. There are no special laws for full name entities and abbreviation entities of countries, provinces, and transliterations etc. E.g. (Shanxi province) and (Jin); (Kentucky Fried Chicken) and KFC

In the following, we introduce our method of extracting full name entity as example, as extracting abbreviation entity is analogous.

Obtaining corpus

A big problem of extracting relation of entities is OOV (out of vocabulary), which is hard to identify. However, we found that full name entity and according abbreviation entity co-exist in a sentence in many cases: AB(B is short for A) and BA(A is full name of B). Inspired by this phenomenon, for a given entity, we can construct query items to search from the Web and obtain relevant corpus.

There are two kinds of query items: forward query item Ent, and backward query item *Ent, in which Ent represents an entity. The query item contains double quotation marks to assure Ent and "/" could appear in the query result consecutively. The queries are committed to Google search engine 4 , we get only the summary of search results, not all the web pages which the result links to. In our experiment we collect top 200 search result summaries for each query item and delete the HTML tags.

F = f 1 f 2 . . . f m , f j (1 ≤ j ≤ m
) is a word, and an abbreviation entity A = a 1 a 2 . . . a n , a i (1 ≤ i ≤ n) is a Chinese character. We define a similarity mark, which is not a strict measure, but used to decide the boundary of candidate from a sentence.

SimMark(A, F )= n i=1 m j=1 j * issub(a i , f j ) (1)
,in which issub(a i , f j ) = 1 if a i appears in f j 0 else For example, we want to get the full names of entity (USTC), we extract a sentence (with segment and Pos-tagging) :

/n /n /j /w /d /v /a /a /n /w /v /w /u /n /w (Central South University Patriotic Health Campaign Committee (full name of "Ai Wei Hui") issue notice about . . . )

We found that SimM ark(, )≤SimM ark(,) ≤SimM ark(,). Therefore, is more possible become a candidate than . The descriptions of algorithm are as follows.

Input Return.

Algorithm 1: Full name/ Abbreviation Extract Algorithm (FAEA)

Filtering & sorting candidates

After using algorithm FAEA, we may get several candidates; however, the Web is so open that some sentences may be irrelevant or even wrong. Therefore, it is necessary to validate and sort all the candidates. We summarize some heuristic rules according to the commonalities between the abbreviation and full name of an entity, and they are used to filter out the candidates as follows.

Denote the set of candidates as Candidates = {Cand 1 , Cand 2 , . . . , Cand n }, Cand ∈ Candidates, five filtering rules are defined as follows: (if one of the rules is satisfied, then Cand is thought as an error, and thus be filtered out.) 1. Candis a single Chinese character; 2. len(Ent) ≥ len(Cand), where len(Ent) and len(Cand) are the numbers of characters in Cand and Ent respectively; 3. There are no common Chinese characters between Cand and Ent; 4. segnum(Cand) > segnum(Ent)+3, where segnum(Cand) and segnum(Ent) are the numbers of words of Cand and Ent with segmentation respectively; 5. Cand contains some meaningless interrogatives words;

The first four rules are straightforward, and we explain the fifth one a little. We may obtain some interrogative sentences which interrogate for the full name or abbreviation of an entity, but which do not end with any question mark, such as (What is the full name of Tsinghua). In this case, Cand = (what), which is actually an error. The fifth rule can identify such errors, and filter them out.

Attributed to the fact that an entity could have more than one full name or abbreviation entities, it's necessary to sort all the candidates using the statistical information from Sents, the candidate with rank one is most common full name or abbreviation. We define a sort comparison function to sort all the candidates. Denote the set of relevant sentences as Sents, and define

C i >C j , meaning C i precedes C j , iff 1. SubstrF req(C i ) ≥ SubstrF req(C j ) 2. LD(C i ,Ent) ≤ LD(C j , Ent), if SubstrF req(C i ) = SubstrF req(C j ) 3. SubseqF req(C i ) ≥ SubseqF req(C j ), if SubstrF req(C i ) = SubstrF req(C j )
and LD(C i ,Ent) =LD(C j , Ent).

In which, SubstrF req(Cand)is the number of sentences which Cand appear in the sentences of Sent as a substring; LD(Cand, Ent) is the Levenshtein distance of Cand and Ent; SubseqF req(Cand) is the number of sentences which Cand appear in the sentences of Sents as a subsequence.

Experiments and Discussion

We collect 300 pairs of full name entity and abbreviation entity as test data, in which about 80% are named entity, e.g. ((Peking University), (Beida)), 20). We compute abbreviation entities and full name entities respectively for each pair. The following tables illustrate our results. The precision of top k in table 3, table 4 and table 5 for full name is defined follows, the recall is either.

P recsion top k = Count top k Count F ull (2)
In which, Count top k represents the number of all correct full names extracted in top k, while Count F ull represents the number of all full names. Table 3 shows the performance when only pattern 1 is used. We can found that performance of extracting abbreviations is higher; partly because the pattern is more efficient for abbreviations, moreover, more boundary information (such as parenthesis and quotation marks) could be supplied when extracting abbreviations than full names. Performance of top 3 is higher than which of top 1; which illustrates effect of our ranking strategy; also prove that some entities have more than one abbreviation or full name. The recall of full names is only 58.4%, we found the reason is that pattern 1 cannot obtain sufficient corpus.

Table 4 shows that after pattern 2 is used, the performance is improved in both precision and recall, especially the recall of full names. The results also confirm us that the patterns accord with the most common expressions when full name and abbreviation entities co-appear. Table 5 shows the performance of our method for named entities, we can find that our method is also efficient. Many organization entities end with a suffix and we called it "suffix abbreviation", because it usually follows many different entities in sentences. E.g. (association),(university). We found that pattern 2 is very efficient for entities end with suffix abbreviation.

Table 6 shows some results with their sort compare function value; we can see that our three-level sort strategy is effective. For example, two abbreviation candidates of (Chinese Academy of Social Sciences) are (CASS) and (Siweiliangyuan); Although their SubstrFreq value is the same, we can sort them using Levenshtei Distance value correctly.

Table 7 illustrates the new pairs of full name entity and corresponding abbreviations found when we implement our algorithm for an entity (see the third column). In most cases, they are extracted when encountering with suffix abbreviation. Some pairs are partially correct and could be amended moreover. It also illustrates that our method can extend itself and get more pairs iteratively.

Furthermore, we summarize some difficulties of extracting full names or abbreviations as follows:

1. Abbreviation entities are usually OOV and difficult to segment. Therefore, the segments of abbreviation candidate supplied us are not reliable. 2. Obtain least but most useful corpora: we should have a tradeoff between more corpora and less Web search. We believe that the recall of extracting full name entities may be improved if we introduce more query items. 3. The prefix of a correct full name entity or abbreviation entity is difficult to identify when extracting candidates. In addition, the length and constitutes of prefix words are hard to be determined when we figure out the left boundary of the candidate.

Conclusions

In this paper, we aim at extracting full names and abbreviations for a given entity from the Web. We propose a method of combining both patterns and rules to solve the problem. In addition, new pairs of full name entity and abbreviation entity can be extracted simultaneously. Our experiment shows our method is efficient. However, there is still more future research to improve this work. For example, how to construct more query items to obtain more corpora, how to use named entity identification method to validate candidates.

  Fig 1 illustrates our overall framework.

Table 1 .

 1 Different actions according to P re

	Entity "E"	"E for short"	Search	Top N results	Clean HTML
	"Full name for E"				Extract Sentences
						Extract Algorithm
	Abbreviations
			Filter Rules
	/Full names		Mark Function	Evaluate	Candidate1
	New Full & Abbr.				Candidate2
	pairs found: (NFull, NAbbr)	Sort Function		Candidate3

Fig. 1. Extract Abbreviations/Full names of Entity E from the Web

3.3 Relevant sentence analysis

For the corpora obtained above, we use two patterns for the two query items respectively to extract sentences which contain the query items Pattern 1 Ent(common or bracket or . . . )(colon or bracket or is or . . . ) We call the sentence matching pattern 1 as a relevant sentence, denoted as S, it can be represented as P re+Ent+P un++Cand+P un, in which P re represent prefix words preceding entity Ent, P un represent punctuation such as comma, and Cand represent preliminary candidate.

Pattern 2 (common or bracket or . . . )(colon or bracket or is or . . . ). . . Ent Likely, relevant sentence S matching pattern 2 can be represented as Cand+P un++P re+P un+Ent. Different strategies are designed according to P re(show in table

1

). We list some examples in table 2.

Table 2 .

 2 Some examples of relevant sentence

	Entity (Peking) (at-	Query Item	Relevant Sentence "". . . . . .	P re= P re = NULL NULL √ √
	tached
	primary
	school) (four	*	1963129, ""	√
	moderni-
	zations) (primary	*		√
	school)
	3.4 Extracting candidates

Pattern 1 and Pattern 2 are analogous, so we describe our method aiming at pattern 1 in the following. We segment 5 relevant sentence P re+Ent +P un ++ Cand+ P un as P 1 P 2 . . . P m EntP unC 1 C 2 . . . C n P un . So our question converted into how to find two boundaries i, j,

1 ≤ i ≤ m, 1 ≤ j ≤ n, the full name entity candidate is C j+1 C j+2 . . . C n , meanwhile, new possible full name & abbreviation entity pair found is ( P i+1 P i+2 . . . P m , C 1 C 2 . . . C j ).

For a full name entity with segmentation

:

  Sentence S(with segment and Pos-tagging): P1P2 . . . PmEntP un C 1 C 2 . . . C n P un Entity: Ent; Relation: ""(Full name) Output: Full name Candidate of Ent: EntF ull; New Full name-abbreviation pair (N ewF ull, N ewAbbr); 1 Initialization: i=m, j=n; 2 If P1P2 . . . Pi has at least a common character with C 1 C 2 . . . C n , then i=i-1; Else go to (4); 3 If Pi is an auxiliary or preposition word, then go to (4); Else i=i-1; 4 If Pi is a geographical entity name (discerning from Pos-tagging) or SimM ark(P i P i+1 . . . P m Ent, C 1 C 2 . . . C n )> SimM ark(P i+1 P i+2 . . . P m Ent ,C 1 C 2 . . . C n ), then i=i-1, go to (2); 5 If CjCj+1 . . . Cn has no common character with Ent , then j=j-1; 6 If SimM ark(Ent,C j C j+1 . . . C n )+ SimM ark( P i+1 P i+2 . . . P m , C 1 C 2 . . . C j-1 )> SimM ark(Ent ,C j+1 C j+2 . . . C n )+ SimM ark( P i+1 P i+2 . . . P m , C 1 C 2 . . . C j ), then j=j-1; Else go to (7); 7 EntF ull=C j+1 C j+2 . . . C n , N ewF ull=Pi+1Pi+2 . . . Pm, N ewAbbr=C 1 C 2 . . . C j ;

Table 3 .

 3 Performance with only P attern 1

			Precision	Recall(of Top 1)	F-measure
	Extract Abbr	Top 1 87.1%	Top 3 92.0%	71.7%	Top 1 78.7%	Top 3 80.6%
	Extract Full	70.5%	79.3%	58.4%	63.9%	67.3%

Table 4 .

 4 Performance with pattern 1 & pattern 2

			Precision	Recall(of Top 1)	F-measure
	Extract Abbr	Top 1 84.7%	Top 3 93.6%	90.2%	Top 1 87.4%	Top 3 91.9%
	Extract Full	77.0%	82.7%	85.3%	80.9%	84.0%

Table 5 .

 5 Performance with pattern 1 & pattern 2 for named entity

			Precision	Recall(of Top 1)	F-measure
	Extract Abbr	Top 1 89.5%	Top 3 95.3%	91.1%	Top 1 90.3%	Top 3 93.2%
	Extract Full	82.3%	87.8%	85.9%	84.1%	86.8%

Table 6 .

 6 Some sort comparison function values of abbreviations extracted

	Entity	Abbreviation Substr LD Candidates Freq	Subseq Rank Freq	Result correct?
	(	17	4	36	One	Yes
	Chinese Academy of	3	2	30	Tow	Yes
	Social Sciences)	3	6	3	Three	No
		7	2	34	One	Yes
	(Doctoral candidate)	6	3	36	Two	No

Table 7 .

 7 New full names and abbreviations identified

	Full name	Abbreviation	Trigger Entity	Result correct?
	(Peking University)		(attached primary	Yes
			school)
	(China)		(mobile)	Yes
	(Beijing's Second)		(court)	No
	(China Association of		(customer service)	No
	Trade in Service)

The nist ace evaluation website: http://www.nist.gov/speech/tests/ace/ace07/.

http://www.google.cn

We use Chinese segment tool ICTCLAS. http://www.ictclas.org/

Acknowledgements This work is supported by the National Natural Science Foundation of China under Grant No.60496326, 60573063, and 60773059; the National High Technology Research and Development Program of China under Grant No. 2007AA01Z325; the Education Commission Program of Beijing under Grant No. KM201010009004.