Kristell Dréau
email: kristell.dreau@ec-nantes.fr

Nicolas Chevaugeon
email: nicolas.chevaugeon@ec-nantes.fr

Nicolas Moës
email: nicolas.moes@ec-nantes.fr

Studied X-FEM enrichment to handle material interfaces with higher order finite element

Keywords: eXtended Finite Element Method Level set, High order approximation, Free surface, Material interface, Curved geometries

An approach to improve the geometrical representation of surfaces with the eXtended Finite Element Method is proposed. Surfaces are implicitly represented using the level set method. The finite element approximation is enriched by additional functions through the notion of partition of unity, to track material interfaces. Optimal rate of convergence is achieved with curved geometries, using linear elements and linear level set in elements. In order to accelerate the convergence, the order of approximation shape functions is increased, while keeping the same computational mesh. The level set is represented on a finer submesh than the finite element mesh. A special attention to integration procedure is necessary. A new enrichment function is introduced to represent the behavior of curved material interfaces. Numerical examples including free surfaces and material interfaces in 2-D linear elasticity are presented to study convergence rates.

Introduction

Mesh generation of complex geometries can be very timeconsuming, within a classical finite element analysis. The main difficulty arises from the necessity of the mesh to conform to physical surfaces. Discontinuities such as holes, cracks and material interfaces may not cross mesh elements. Moreover, local refinements close to discontinuities and mesh modification to track the geometrical and topological changes in crack propagation problems for example, can be difficult. Also, when geometries evolve and history dependent models are used, robust methods to transfer the solution to the new mesh are needed. This issue is particularly significant, since computed fields defined on these discontinuities are often the most important ones.

In order to improve the geometrical representation of surfaces, and particularly with complex geometries including curves, a solution is to locally refine the mesh, close to discontinuities, i.e. to decrease element size. Several methods of adaptive h-refinements were developed (for example [START_REF] George | Mesh Generation, Application To Finite Elements[END_REF][START_REF] Shepard | Adaptive mesh generation for curved domains[END_REF][START_REF] Szabó | Finite Element Analysis[END_REF]), by subdividing elements or remeshing. The position of nodes can also be adjusted in order to improve the geometrical representation while keeping their number constant; this is called r-refinement.

Another mean to reduce geometrical errors is to introduce curved elements in the mesh. Elements with curved boundaries are used to represent curved geometries, increasing the degree p of geometrical shape functions and the number of geometrical nodes on boundaries. This technique allows the use of meshes with fewer elements. Adaptive p-refinement, also called p-method, is presented in [START_REF] Szabó | Finite Element Analysis[END_REF]. Usually, these shapes functions are constructed with Legendre polynomials. With iso-parametric element, shape functions are the same for the field approximation and the geometry. Curvilinear mesh generation encounters two main issues [START_REF] Dey | Curvilinear mesh generation in 3d[END_REF][START_REF] Sherwin | Mesh generation in curvilinear domains using high-order elements[END_REF]: the respect of the geometrical representation, and the generation of a valid mesh. For example, gaps and overlaps have to be removed from meshes. Integration of functions on these elements must be computed with caution [START_REF] Királyfalvi | Quasi-regional mapping for the p-version of the finite element method[END_REF][START_REF] Szabó | Finite Element Analysis[END_REF]. If an order p increases, oscillations may nevertheless appear.

Recent studies about p-version mesh generation develop the use of Bezier polynomials to represent curves and surfaces [START_REF] Luo | Automatic p-version mesh generation for curved domains[END_REF][START_REF] Luo | p-version mesh generation issues[END_REF]. Increasing the geometric approximation order and subdividing elements become easier, due to properties of Bezier polynomials. New mapping procedures have to be defined. However, mesh generation and geometrical approximation still have to be controlled.

Isogeometric analysis [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, nurbs, exact geometry and mesh refinement[END_REF] gets rid of Lagrange polynomials drawbacks: geometry is exactly represented by NURBS and refinements simplified. However, this geometrical benefit leads to additional difficulties in numerical integration [START_REF] Hughes | Efficient quadrature for nurbs-based isogeometric analysis[END_REF].

Several numerical methods were developed, in order to get rid of mesh difficulties. An idea was, for example, to base approximations on a set of nodes without creating any mesh. These kinds of methods are called Meshfree Methods [START_REF] Fries | Classification and Overview of Meshfree Methods[END_REF]. Therefore, hp-adaptivity becomes easier because of the absence of a mesh [START_REF] Duarte | An h-p adaptive method using clouds[END_REF]. However, difficulties arise from numerical integration and treatment of essential boundary conditions.

For sake of mesh simplifications, the eXtended Finite Element Method (X-FEM) was developed 10 years ago, initially for crack propagation [START_REF] Belytschko | Elastic crack growth in finite elements with minimal remeshing[END_REF][START_REF] Daux | Arbitrary branched and intersecting cracks with the extended finite element method[END_REF][START_REF] Moës | A finite element method for crack growth without remeshing[END_REF].

Within the X-FEM, surfaces that are not represented explicitly by mesh boundaries can be implicitly represented by the iso-zero values of a level set function. This is particularly useful for moving interfaces such as crack surfaces in crack propagation analysis [START_REF] Sukumar | Extended finite element method for three-dimensional crack modelling[END_REF]. The finite element approximation is enriched by additional functions through the notion of partition of unity [START_REF] Melenk | The partition of unity finite element method: basic theory and applications[END_REF] to represent, for example, discontinuities at interfaces [START_REF] Sukumar | Modeling holes and inclusions by level sets in the extended finite-element method[END_REF] or asymptotic expansions terms near crack tip [START_REF] Sukumar | Extended finite element method for three-dimensional crack modelling[END_REF] to improve convergence rate. X-FEM is a good alternative to the classical finite element method for solving mechanical problems with complex or evolving geometries.

Enrichment functions are usually defined with the help of the level set functions to access the distance to the interface at any given point. Although interfaces do not have to be meshed, the correct integration of the element stiffness matrix, for elements enriched by a discontinuous function along an interface, needs to be done carefully. These elements are split for integration purpose along the iso-zero level set, creating a set of integration cells called partition. A regular integration is done separately on each cells of the partition.

In the case of free surfaces [START_REF] Sukumar | Modeling holes and inclusions by level sets in the extended finite-element method[END_REF], the level set is only used to construct the partition and to determine if integration cells of this partition are in the material or in the hole. On the contrary, material interfaces and cracks need additional functions for displacement field [START_REF] Daux | Arbitrary branched and intersecting cracks with the extended finite element method[END_REF][START_REF] Moës | A computational approach to handle complex microstructure geometries[END_REF] depending on level set values on element nodes.

Difficulties to represent the surface in the classical finite element method are, in the context of X-FEM, partly shifted on the integration procedure. A method for eliminating the subdivision of elements during the integration of the stiffness matrix is presented in [START_REF] Ventura | On the elimination of quadrature subcells for discontinuous functions in the extended finite element method[END_REF]: discontinuous and non-differentiable functions are replaced with equivalent polynomials (exact equivalence for triangular and tetrahedral elements). However, it applies only for linear shape functions and a straight discontinuity across elements. Nevertheless, X-FEM gives good results with linear elements and linear level set in elements: optimal rate of convergence is achieved with curved geometries [START_REF] Sukumar | Modeling holes and inclusions by level sets in the extended finite-element method[END_REF].

In order to accelerate the energy error convergence, p-versions of X-FEM were developed. In [START_REF] Iarve | Mesh independent modelling of cracks by using higher order shape functions[END_REF][START_REF] Laborde | High-order extended finite element method for cracked domains[END_REF][START_REF] Stazi | An extended finite element method with higher-order elements for curved cracks[END_REF] for example, high order X-FEM for straight cracks also reaches an optimal rate of convergence. If the crack has not a flat surface, the error is reduced with shape functions of higher order without improving the geometrical representation, but the optimal rate of convergence is not achieved. The problem is of course related to the well known fact that, in the framework of classical finite element method, when increasing the order of the approximation field, the order of representation of the geometry must be increased accordingly to get optimal rate of convergence. Isoparametric elements have an energy error norm on regular problems with curved boundary that converges as O(h p ), where h is the element size and p the polynomial basis order (see [START_REF] Zienkiewicz | The Finite Element Method[END_REF]).

To get optimal rate of convergence when curved interfaces are represented with level set, the obvious solutions would be to increase the order of the level set representation near the iso-zero values of the level set, as in [START_REF] Unger | Simulation of concrete using the extended finite element method[END_REF]. In this case, the integration problem can become quite difficult: within an element, the iso-zero surface can have complex shape and topology. In [START_REF] Legay | Strong and weak arbitrary discontinuities in spectral finite elements[END_REF], X-FEM is developed with spectral elements, to deal with arbitrary discontinuities: optimal rate of convergence in energy norm are not reached for curved discontinuities, and integration is more delicate. The influence of the geometry description within the Finite Cell Method is investigated in [START_REF] Düster | The finite cell method for threedimensional problems of solid mechanics[END_REF] for three-dimensional problems with curved boundaries: using an implicit representation and using a voxel model that is promising for biomechanics simulations.

X-FEM is extended to high order approximations in [START_REF] Cheng | Higher-order xfem for curved strong and weak discontinuities[END_REF] for arbitrary curved strong and weak discontinuities. Elements crossed by a discontinuity are divided into subcells with one curved side. Optimal convergence is reached in the first case and a close-tooptimal one in the second case using the Corrected X-FEM [START_REF] Fries | A corrected xfem approximation without problems in blending elements[END_REF].

The alternative approach proposed here consists in representing the level set on a finer sub-mesh than the finite element mesh, while keeping the level set representation piecewise linear by element (without adding any degree of freedom to the mechanical field). This is used in [START_REF] Marchandise | A stabilized finite element method using a discontinuous level set approach for the computation of bubble dynamics[END_REF] to localize an interface. This idea was also developed in hp Generalized FEM [START_REF] Pereira | hp-generalized fem and crack surface representation for non-planar 3-d cracks[END_REF]: representation of crack surfaces is controlled by a Delaunay tetrahedralization of elements cut by the discontinuity.

X-FEM is presented in the next section, with a brief recall of linear elasticity. A new strategy is then introduced to use high order X-FEM and applied to domains with curved boundaries. Numerical applications of free surfaces and material interfaces are developed within the framework of linear elasticity with polynomial approximations up to order three.

Extended finite element method

Model problem

The formulation is written within the framework of linear elastostatic. Consider a domain Ω, with boundary Γ, representing a linear elastic, homogeneous and isotropic solid. The domain is submitted to impose displacement ū on Γ u and imposed tractions t̄on Γ t . We study the case without volumic force, in small strain. Stresses, strains and displacements are respectively denoted by σ, ε and u.

The problem is to find the admissible displacement field which satisfies the following equations:

• Equilibrium equations:

∇ ⋅ σ = 0 on Ω σ ⋅ n = t on Γ t ;with n the external normal • Kinematics equations: εðuÞ = 1 2 ∇u + ∇u ð Þ T on Ω u = u on Γ u
• Constitutive law: σ = C : ε; on Ω; where C is Hooke 0 s tensor Let U be the space of admissible displacement field, and U 0 the space of kinematically admissible test functions, vanishing on Γ u . The weak form of equilibrium elastostatic and boundary conditions can be written as:

Find u ∈ U such that

∫ Ω εðuÞ : C : εðυÞdΩ = ∫ Γ t t⋅υdΩ; ∀υ ∈ U 0 2.2.

eXtended Finite Element Method

Within the finite element method, the domain Ω is divided into elements, called finite elements and constituting a mesh. The displacement field is approximated on these elements:

uðxÞ = ∑ n i ϕ i ðxÞu i
where ϕ i (x) are the finite element shape functions, and u i the nodal displacements.

Within X-FEM, the mesh does not necessarily conform to internal discontinuities (such as holes, cracks, and material interfaces), which allows the use of simple meshes (structured or not). Moreover, the finite element approximation is enriched by additional functions (through the partition of unity [START_REF] Melenk | The partition of unity finite element method: basic theory and applications[END_REF]) that model the behavior of boundaries, as follows:

u h ðxÞ = ∑ i∈N ϕ i ðxÞu i + ∑ j∈N g ϕ j ðxÞFðxÞa j ð1Þ
where F(x) is the enrichment function; a j are additional degrees of freedom, for enriched nodes; N represents the set of the nodes of the mesh, and N g the set of enriched nodes. For example, in the case of material interfaces, only nodes belonging to elements crossed by the interface are enriched. If the mesh conforms to discontinuities, there is no enrichment on nodes, and classical finite element method is recovered.

As surfaces are not represented explicitly by mesh boundaries, a function, called level set φ(x, t) [START_REF] Osher | Level Set Methods and Dynamic Implicit Surfaces[END_REF], is introduced to track them. This function usually gives the signed distance to the interface at any point of the domain (Fig. 1). Consequently, the interface location is given by the iso-zero value of the level set: φ(x, t) = 0.

High order X-FEM with complex geometries

A new approach

For the moment, X-FEM is mostly used with linear approximation shape functions ϕ i (x) and linear level set. Optimal (or nearly optimal) rates of convergence are reached in the case of holes, material interfaces and cracks [START_REF] Laborde | High-order extended finite element method for cracked domains[END_REF][START_REF] Moës | A computational approach to handle complex microstructure geometries[END_REF][START_REF] Sukumar | Modeling holes and inclusions by level sets in the extended finite-element method[END_REF][START_REF] Sukumar | Extended finite element method for three-dimensional crack modelling[END_REF] for straight cracks).

Our purpose is to accelerate the convergence of X-FEM, using high order shape functions for the finite element approximation. However, the optimal rate of convergence can't be reached without improving the geometrical representation, geometrical errors becoming more important than approximation errors.

Instead of using an exact representation of the domain (as [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, nurbs, exact geometry and mesh refinement[END_REF] for example), the finite element method utilizes a polynomial representation of the geometry. The first idea for the representation of complex geometries with X-FEM would be to use elements with curved boundaries.

These elements can be easily generated by adding geometrical nodes on edges, with hierarchical shape functions. p-refinement does not need to be uniform throughout the domain. Energy error norm on regular problems converges as O(h p ) with iso-parametric elements, where h is the element size and p the polynomial basis order. However, the mesh generation might become quite difficult with this kind of elements. Indeed, Fig. 2 illustrates a classical problem of geometrical representation: when the edge 1 is improved, overlap can appear.

In reference [START_REF] Unger | Simulation of concrete using the extended finite element method[END_REF], level set of order two are used to model interfaces, providing a more accurate representation. In [START_REF] Stazi | An extended finite element method with higher-order elements for curved cracks[END_REF], higher order elements are applied to curved cracks problems (six-node triangular elements) with quadratic shape functions. Although the energy error is reduced, the rate of convergence is not improved. Furthermore, it is shown that the use of linear shape functions for the enrichment function gives better results, even though higher order shape functions are used for the continuous approximation.

Instead of using a level set representation of higher order, leading to complex shapes of iso-zero surfaces and integration difficulties, a strategy is presented to get closer to the real geometry of the discontinuity while keeping a linear representation of the geometry. This consists in representing the level set on a finer mesh than the finite element mesh, without additional degree of freedom for the mechanical field. Thus, representation of discontinuities such as holes boundaries, material interfaces or cracks is improved in a simple way.

Geometrical representation using level sets

The level set is interpolated on elements with finite element shape functions. These geometrical shape functions are generally the same than approximation shape functions (case of iso-parametric elements). For now, within X-FEM, the level set is linear on elements (linear shape functions). Therefore, the mesh has to be fine enough to accurately locate the interface.

The aim is to improve the geometrical representation while keeping the same computational mesh, and therefore without adding any degree of freedom. The level set is computed on a finer mesh than the computational mesh, called "geometrical grid", and its elements are called "cells". Initially, the geometrical grid is a copy of the computational mesh. If the geometrical representation has to be improved, the geometrical grid has to be refined. Thus, each cell crossed by the iso-zero is split into cells of same nature. For example, a triangle in 2D is divided in 4 triangles (quadtree refinement). Level set values are then computed on the new nodes. The procedure is done recursively up to a user defined maximum depth, forming the final geometrical grid (see Fig. 5). Each element of the analysis is then linked to an octree-like partition where the level set is represented and on which the integration procedure is applied for each sub-element. The approximation of the iso-zero surface is then improved at each level of refinement of the geometrical grid. The interface mesh representing the boundary of a hole is plot on Fig. 3, with no refinement, and a level of refinement equal to 1 and 2. Refining the geometrical grid improves the representation of the interface while keeping a level set piecewise linear on elements.

Numerical integration

Numerical integration has to be handled with care in X-FEM. Indeed, discontinuous functions can't be directly integrated on elements. Therefore, cells of the geometrical grid cut by the iso-zero value of the level set are split in integration cells (Fig. 4), creating partitions tied to mesh elements, as presented in [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF]. A regular integration is thus done in these integration cells (Gauss integration), on both sides of the discontinuity. For example, an element crossed by the boundary of a hole is split in integration cells, and integration is performed only on the cell belonging to the material (i.e. not in the void).

If the geometrical grid is refined according to the previous description, integration operations are applied to cells: integration cells are created, and gathered in partitions linked to mesh elements. This process is explained in Fig. 5.

Finally, partition cells are classified depending on their distances to the interface. For example, in the case of holes, partition cells are classified in two groups: cells in the material or in a hole. Thus, integration is performed only on partition cells in the material, those in holes being ignored.

Next sections prove the validity of this high order approach. The study of benchmark elasticity problems allows one to verify the convergence of energy norm with high order X-FEM. First, the case of curved free surfaces is observed with an infinite plate with a circular hole under uniaxial tension. The case of curved material interfaces presents a more sensitive issue and is afterward studied with circular inclusion. All results are presented using GMSH [START_REF] Geuzaine | Gmsh: A Three-Dimensional Finite Element Mesh Generator with Built-in Pre-and Post-Processing facilities[END_REF] as a post-processing tool.

Free surfaces

First, the treatment of free surfaces with X-FEM is considered. In order to check the validity of this proposed method, we are interested in studying the convergence of the solution for an infinite plate with a hole under uniaxial tension. Thus, the classical finite element method and the extended finite element method can easily be compared, since an analytical solution exists.

The case of free surfaces can be easily handled with X-FEM [START_REF] Daux | Arbitrary branched and intersecting cracks with the extended finite element method[END_REF][START_REF] Sukumar | Modeling holes and inclusions by level sets in the extended finite-element method[END_REF]. Nodes for which the support is totally inside voids are wiped out of the computation.

Infinite plate with a hole under uniaxial tension

We consider an infinite plate with a circular hole under uniaxial tension σ ∞ = 1 along axis x (see Fig. 6). This plate is modeled by a square domain of side L with a circular hole of radius a in its center. Exact tractions are imposed on boundaries of the domain.

Exact solution is shown in [START_REF] Szabó | Finite Element Analysis[END_REF]. Components of exact stress tensor are: σ xy ðr; θÞ = -

σ xx ðr; θÞ = 1- a 2 r 2 3 2 cos 2θ ð Þ + cos 4θ ð Þ ð Þ + 3 2 a 4 r 4 cos 4θ ð Þ ð2Þ
σ yy ðr; θÞ = - a 2 r 2 1 2 cos 2θ ð Þ-cos 4θ ð Þ ð Þ - 3 2
a 4 r 4 cos 4θ ð Þ ð3Þ
a 2 r 2 1 2 sin 2θ ð Þ + sin 4θ ð Þ ð Þ + 3 2 a 4 r 4 sin 4θ ð Þ ð4Þ
And exact displacements read: • μ: shear modulus • κ: Kolosov constant κ = 3 -4ν (simulations are carried out in plane strain)

Results

In order to compare precision of the finite element solution with the one computed with X-FEM, convergence of the relative error in the energy norm is studied, for square regular meshes with 8, 16, 32, 64 and 128 elements (triangular elements) by side. The side length is 2. Numerical experiments are done with polynomial approximation up to order three. Energy error norm is given by:

e = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ∫ Ω εðu h -u ex Þ : C : εðu h -u ex ÞdΩ q ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ∫ Ω εðu ex Þ : C : εðu ex ÞdΩ q ð7Þ
Fig. 7 presents the relative energy error functions of the mesh element size, convergence rates being computed from errors of the two finest meshes (with 64 and 128 elements). Using high order finite element approximation without improving the geometrical representation does not give good results, in term of convergence rate of the energy error norm. The use of approximation functions of orders 2 and 3 improves the convergence rate compared with an approximation of order 1. However, it does not reach the classical finite element rate (Fig. 7), that is O(h p ). We note that the quadratic approximation reaches a convergence rate quite close to the finite element one: average rate of 1.87. On the contrary, the cubic approximation converges with a rate of 1.62. Note that these rates correspond to a particularly penalizing case. Indeed, using a circular inclusion radius of 0.4, the iso-zero value of the level set passes over some elements nodes, which is slightly penalizing.

Using the present method, improving the geometrical representation of the boundary of the hole reduces geometrical errors and consequently global error. Fig. 8 presents the convergence of the relative energy error functions of the mesh element size, with cubic approximation shape functions, for different level of precision of the geometrical representation (recursive refinement of the geometrical grid). The more refined is the geometrical grid, the smaller is the relative error norm. Choosing an optimal refinement level, depending on the order of approximation p and on the element size h, allows to accelerate the convergence of energy error functions of the mesh element size (see Fig. 9): energy error converges with a quasi-optimal rate of 2.93.

We note that the condition number could be deteriorated if elements or part of them are located within the hole. This is why elements completely in the hole are removed from the computation, and integration is performed only in partition cells in the material. Moreover, a loop is done on mesh edges. If the intersection point between an edge and the iso-zero value of the level set is very close to one of the two nodes (distance smaller than 10 -3 h, with h the edge size), the level set value of the node is modified and equal to zero. This node is "moved" on the interface. Then, precision of computations is assured by a direct solver. Finally, condition number is not deteriorated.

Thus, we show that optimal convergence rate can be achieved with an optimal geometrical representation of the boundary of the hole.

Global error is made up of interpolation error and geometrical error. Approximation error decreases as (h a ) p , h a being the relative approximation element size, and p the degree of approximation shape functions. Numerical experiments show that geometrical error decreases as (h g ) 2 , h g being the relative geometrical cell size. With a linear approximation and a linear geometry on X-FEM elements (and also for iso-parametric elements), relative approximation and geometrical elements are the same: h a = h g . This relation is no longer correct with higher order approximation. A method to find the correct Fig. 6. Infinite plate with a hole under uniaxial tension. Fig. 7. Energy error norm for a plate with a hole, with orders of approximation 1, 2 and 3. Fig. 8. Energy error norm for a plate with a hole using cubic approximation shape functions and different refinement level of the geometrical grid. geometrical cell size, that is to say the correct refinement level, is presented.

Consider a mesh of relative element size h a . In order to find the geometrical cell size, let us consider geometrical and approximation errors orders to be equal:

h 2 g = h p a
Therefore, the relative geometrical cell size is:

h g = ðh a Þ p 2 = 1 n p 2 = h a 1 n p 2 -1 ð8Þ
For example, considering quadratic shape functions (p = 2), we have h g = h a . Indeed, the relative energy error converges with the optimal rate, without the need to improve the geometry.

With cubic shape functions (p = 3), the relation becomes:

h g = 1 n 1 2 h a ð9Þ
Fig. 10 shows theoretical division of the approximation element and the chosen refinement level for each mesh, as element size can only be divided by 2 r during geometrical refinement. The smaller is the mesh element size, the bigger is the refinement level. For example, the geometrical grid is refined twice using cubic shape functions and a square mesh with 32 elements by side.

This relation just provides h g = h g (h a ), to maintain the optimal convergence rate. However, it does not depend on the curvature of the interface. A straight interface does not need refinement, as geometry is correctly represented with linear level set. On the contrary, a highly curved interface has to be well represented on elements thanks to a fine geometrical grid. A reflexion has to be lead on a correct rule to define the refinement level depending on local curvature of interface.

Thus, dealing with free surfaces, the proposed method allows one to use higher order approximation shape functions while keeping the optimal convergence rate of classical finite element methods while avoiding the generation of curved meshes.

Material interfaces

X-FEM for material interface

The X-FEM approximation requires an enrichment function, able to represent the behavior of an interface between two materials. Indeed, as the mesh does not conform to the interface, some elements overlap the two materials. An enrichment term is added to the finite element approximation on nodes of elements cut by the interface, introducing additional degrees of freedom a j , and an enrichment function F(x).

u h ðxÞ = ∑ i∈N ϕ i ðxÞu i + ∑ j∈N g ϕ j ðxÞFðxÞa j
where u i are the degrees of freedom and ϕ i (x) the approximation shape functions. Elements cut by the discontinuity have to be split into sub-elements belonging to one or the other material, for numerical integration.

Displacement is continuous across a material interface, with a discontinuous gradient. In the case of material interface, the level set φ(x) is not only used to track the geometry of the interface, but it is also used to build enrichment functions. In [START_REF] Moës | A computational approach to handle complex microstructure geometries[END_REF], an enrichment strategy is proposed to model discontinuities along material interfaces. The enrichment function [START_REF] Geuzaine | Gmsh: A Three-Dimensional Finite Element Mesh Generator with Built-in Pre-and Post-Processing facilities[END_REF] was tested for a 2D circular and a 3D spherical inclusion under tension in an infinite domain. This function is continuous, with a discontinuous derivative across the interface (Fig. 11). It is maximal on the interface, and equal to zero on Fig. 9. Energy error norm for a plate with a hole using orders of approximation 1, 2 and 3, improving the geometrical representation. Fig. 10. Free surfaces: refinement level for a cubic approximation. element nodes, edges and faces not crossed by this interface. The convergence rate of energy error is observed as close-to-optimal one.

FðxÞ = ∑ i j φ i j ϕ i ðxÞ -j ∑ i φ i ϕ i ðxÞ j ð10Þ
with φ i the level set value on element nodes.

High order X-FEM for straight interfaces

Before studying curved material interfaces, we are interested in using high order approximation shape functions simply with straight material interfaces. In this case, the transition to high order X-FEM is obvious. Indeed, as geometry is represented by a linear function, the enrichment function has to remain the same when the geometrical grid is refined. Thus, the previous enrichment term [START_REF] Moës | A computational approach to handle complex microstructure geometries[END_REF] is still convenient.

For example, consider a plate of side L = 2, and section S with two materials 1 and 2: E 1 = 1, E 2 = 10, ν 1 = ν 2 = 0. Plain-strain condition is assumed. It is embedded on the left side and a variable volumic load is imposed on the plate (Fig. 12):

∥ → f ∥ = f x L 2 .
Exact displacements in the two materials are:

u 1 ðxÞ = f 3E 1 SL 2 L 3 x- x 4 4 x u 2 ðxÞ = f 3E 2 SL 2 L 3 x- x 4 4 + 31f L 2 192S 1 E 1 - 1 E 2 x
Computations are done on squared irregular meshes, with 5, 10, 20, 40, and 80 elements (triangular elements) by side (L = 2), with linear, quadratic and cubic polynomial approximations. Element edges of theses meshes do not align with the straight interface, as shown in Fig. 13: iso-zero value of the level set is plot on a mesh with 10 elements by side. Results are presented in Fig. 14: energy norm converges with close-to-optimal rates, when higher order approximations are used.

New enrichment strategy for curved material interfaces

Enrichment described in Eq. ( 10) is used with a linear representation of the geometry. In order to use high order approximation functions, the same strategy than free surfaces is applied for curved material interfaces. That is to say that the geometrical grid is locally refined while keeping a piecewise linear level set representation by element. The level set is computed on the geometrical grid, so that its representation is more precise than if it is computed on the original mesh (without refinement). However, the previous enrichment function is not suitable anymore because of a gap between the position of the iso-zero value of the level set and the discontinuity on the element (Fig. 15). A new enrichment strategy is developed to solve this problem and to obtain optimal convergence.

Elements that are not crossed by the iso-zero value of the level set are not enriched. As the enrichment function is continuous over the domain, it has to be equal to zero on nodes, edges and faces, whose adjacent elements are not enriched. The derivative is discontinuous across the interface, whose position is given by the iso-zero value of the level set. Moreover, it has to be valid whatever the refinement level of the geometrical grid. This new enrichment is here proposed for triangle elements in 2D. It will be presented for tetrahedron elements in 3D in a forthcoming paper.

Different configurations are considered, functions of the position of the iso-zero value of the level set in relation to the element.

To build such function, we start by something similar to the enrichment in the linear case: the enrichment term depends on both level set nodal values and its values on the geometrical grid.

FðxÞ = ∑ i∈mesh j φ i j ϕ i ðxÞ-j ∑ j∈octree φ j ϕ j ðxÞ j ð11Þ
This function [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, nurbs, exact geometry and mesh refinement[END_REF] has a discontinuous derivative across the iso-zero value of the level set, and is equal to zero on mesh element nodes. It is different to zero on edges, and particularly on edges that are not cut by the iso-zero value of the level set. However, mesh elements that are not crossed by the interface are not enriched. Therefore, enrichment function on edges that are not cut has to be equal to zero. Theses edges are thereafter called "non-cut edges". Fig. 16 shows that this function is no more canceled on non-cut edges, when the geometrical grid is refined.

A correction term has to be added in order to correct the enrichment function, and is called "edge bubble function" (Fig. 17).

FðxÞ = ∑ i∈mesh jφ i j ϕ i ðxÞ |fflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflffl} ① -j ∑ j∈octree φ j ϕ j ðxÞ j |fflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflffl} ② -Π lin φ NonÀcut edge Π edge x -φ NonÀcut edge Π edge x ðϕ 1 ðxÞ + ϕ 2 ðxÞÞ |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} ③: edge bubble function ð12Þ
This enrichment Eq. ( 12) is computed for an element "e" of the mesh (computational mesh) which is crossed by the iso-zero value of the level set. The first sum is computed on nodes i of the element "e". This element is split into cells called"e r ", creating the octree (part of the geometrical grid). The second sum is computed on nodes j of each element e r . Shape functions of the term (ϕ 1 (x) + ϕ 2 (x)) are elements shape functions. The "edge bubble function" corresponds to the level set computed on edge [START_REF] Hughes | Efficient quadrature for nurbs-based isogeometric analysis[END_REF], of which the linear part is removed, and which is "stretched" to the opposite node 3.

Example. Consider a triangle element (see Fig. 18) which is crossed by the iso-zero value of the level set. We are interested in the number of non-cut edges. The iso-zero value of the level set cut two edges, so there is only one non-cut edge: edge [START_REF] Hughes | Efficient quadrature for nurbs-based isogeometric analysis[END_REF]. It means that an edge bubble function corresponding to edge ( 12) is added is in the enrichment function.

Following notations are used. Let G be a Gauss point of the element, with local coordinates (u, v) and barycentric coordinates (λ 1 , λ 2 , λ 3 ) (Fig. 19). The projection of this point on the edge [START_REF] Hughes | Efficient quadrature for nurbs-based isogeometric analysis[END_REF] gives a distance χ = λ 2 λ 1 + λ 2 from node 1. Expressions of theses barycentric coordinates are the same than triangle shape functions (ϕ 1 (x), ϕ 2 (x), ϕ 3 (x)), so barycentric coordinates are used to simplify the enrichment formulation: FðxÞ = j φ 1 j λ 1 + j φ 2 j λ 2 + j φ 3 j λ 3 -j ∑ j∈octree φ j λ j j -j φ 1 j 1-

λ 2 λ 1 + λ 2 + j φ 2 j λ 2 λ 1 + λ 2 - j φ ð12Þ λ 2 λ 1 + λ 2 j " # ×ðλ 1 + λ 2 Þ
φ (12) is the representation of the level set on the edge ( 12) , with edge shape functions. The value depends on the projection of the Gauss point on this edge χ = λ 2 λ 1 + λ 2 (Fig. 19). Following notations are used:

• φ i : level set value on node i, • λ i (x) = λ i : barycentric coordinate of a point x in the element e, • χ: projection of a Gauss point on an edge (ij)

• φ (ij) (χ): level set value on edge (ij), depending on the projection χ of the Gauss point on this edge, All level set values are computed from the geometrical grid. Finally, we obtain:

FðxÞ = jφ 3 j λ 3 + j φ ð12Þ λ 2 λ 1 + λ 2 j ðλ 1 + λ 2 Þ-j ∑ j∈octree φ j λ j j
Enrichment function directly depends on edge bubble functions, and thus on the number of edges that the iso-zero value of the level set does not cut. Different configurations are presented in Fig. 20. We could even imagine a discontinuity located in one element, provided that the octree is enough refined: the enrichment function is enriched by three edges bubble functions.

Circular inclusion with imposed displacement

The following problem is a circular inclusion Ω 1 in the center of a material Ω 2 . A linear displacement is imposed on the boundary Γ 2 (r = b). Material parameters are: E 1 = 1, ν 1 = 0.25 (equivalent to the Lamé constants λ 1 = μ 1 = 0, 4), and E 2 = 10, ν 2 = 0, 3 (equivalent to the Lamé constants λ 2 = 5, 769 and μ 2 = 3, 846). We work with polar coordinates (r, θ).

The exact displacement solution reads:

u r ðrÞ = ( 1- b 2 a 2 α + b 2 a 2 r if 0 ≤ r ≤ a r- b 2 r α + b 2 r if a ≤ r ≤ b u θ = 0 with: α = ðλ 1 + μ 1 + μ 2 Þb 2 ðλ 2 + μ 2 Þa 2 + ðλ 1 + μ 1 Þðb 2 -a 2 Þ + μ 2 b 2
Displacement is continuous, but with a discontinuous derivative across the interface.

And the strain solution is:

ε rr ðrÞ = 1- b 2 a 2 ! α + b 2 a 2 if 0 ≤ r ≤ a 1 + b 2 r 2 ! α- b 2 r 2 if a ≤ r ≤ b 8 > > > > > < > > > > > : ε θθ ðrÞ = 1- b 2 a 2 ! α + b 2 a 2 if 0 ≤ r ≤ a 1- b 2 r 2 ! α + b 2 r 2 if a ≤ r ≤ b 8 > > > > > < > > > > > :
Stresses are defined from strains:

σ

Results

The numerical model is a square domain Ω 2 of side L = 2, with a circular inclusion Ω 1 of radius a at its center. Geometrical parameters are: b = 2 and a = 0.4. Exact tractions are imposed on the boundary on the domain. Computations are performed with the same meshes than those for the plate with a hole.

Convergence of relative error of energy norm toward elements size, when initial enrichment function ( 10) is used, is presented in Fig. 21. As expected, the error decreases when quadratic and cubic approximation shape functions are used. Whereas, the precision gain is lower than the previous case of free surfaces. Optimal rates of convergence are not reached: rates of 0.93 with p = 2, and 0.90 with p = 3 (rates computed from energy error values of the two finest meshes).

Fig. 22 shows strains on the plate, computed using the new enrichment function [START_REF] Hughes | Efficient quadrature for nurbs-based isogeometric analysis[END_REF]. Strains are presented with the same coarse mesh (8 elements by side), but with different refinement levels of the geometrical grid: first with no refinement, then with a geometrical grid refined once and twice. Improving the geometrical representation improves the strain field around the inclusion.

As displacements are continuous, the enrichment function has also to be continuous on the domain. Therefore, the refinement of the geometrical grid is no more recursive, but homogeneous per element. Fig. 23 presents integration cells of a square mesh with 8 elements by side, created from a geometrical grid twice refined.

Results using the proposed enrichment method Eq. ( 12) are plotted on Figs. 24 and 25 according to elements size, for approximation orders 1, 2 and 3. Initially, X-FEM converges with high order, but slower than using classical FEM, that is to say O(h p ), where h is the element size and p the polynomial basis order. Improving the geometrical representation reduces energy error. However, contrary to the previous application (plate with a hole), optimal rates of convergence cannot be reached with an optimal depth of refinement. Indeed, the required refinement level of the geometrical grid is more important with a material interface; but for a given mesh, energy error reaches a limit when the refinement level increases. Nevertheless, the improvement of the geometrical representation and the proposed X-FEM enrichment enable the use of higher order approximation shape functions, reducing considerably energy error.

Material interface with geometric singularities

Up to now, only smooth interfaces were considered. The proposed enrichment function can also be applied to inclusion with geometric singularities as polygonal shape. In order to compare results from classical FEM and X-FEM, a plate with a squared inclusion under uniaxial tension is studied.

The numerical model is a square domain of length L = 2, with a square inclusion of side L = 0.57 and of center (0, 0.1) (Fig. 26). The finite element solution is computed on a fine mesh with 400 triangular elements by side. Its element edges respect the material interface, and particularly the four corners. For example, a finite element mesh with 16 elements by side is shown in Fig. 27. The finite element solution is the reference solution. X-FEM meshes are square regular grid with 8, 16 (Fig. 27), 32, 64, and 128 triangular elements by side. All meshes are generated with GMSH. Corners of the material interface are contained in elements of X-FEM meshes, in order to study the influence of the enrichment function.

Fig. 28 shows enrichment function computed on a regular mesh with 16 elements by side, with a geometrical grid twice refined. Using the initial enrichment function [START_REF] Moës | A computational approach to handle complex microstructure geometries[END_REF], corners of the interface are truncated if elements do not respect this interface. By contrast, refining the geometrical grid allows the enrichment function to introduce a discontinuous derivative closer to the real position of the material interface.

X-FEM solutions are compared to the finite element reference solution. Energy error norms with linear, quadratic and cubic shape functions are shown in Figs. 29-31. This enrichment only introduce a discontinuous derivative in the displacement field. Geometrical singularities as corners of the squared interface lead to singularities in the field. The material interface enrichment cannot represent these singularities. This is why refinement only can reduce energy error norm, and is not able to improve convergence rates.

To obtain optimal rates of convergence, singular functions should be added to the finite element approximation, to be able to represent the behavior of the squared material interface. For example, to represent a crack propagation, the displacement field is enriched by a discontinuous field along a crack and by asymptotic fields near crack tips [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF].

Conclusion

This article studies the possibility to use higher order approximation shape functions with X-FEM, when geometries include curved discontinuities. In particular, we are interested by curved free surfaces and material interfaces. A geometrical improvement is done refining the level set representation. Moreover, a new enrichment function is required in the case of material interfaces. Convergences studies are done with two numerical examples in two dimensions: a plate with a hole under uniaxial tension, and a circular inclusion with imposed displacement. Results show that improving the representation of geometry allows to reduce energy error. Considering the plate with a hole, rates of convergence are close to classical FEM ones, that is to say O(h p ), with h the element size, and p the order of approximation shape function. This approach gets rid of difficulties met during mesh generation, as mesh do not have to respect discontinuities with X-FEM. Operations to represent geometries are simplified, and complex geometries can be represented on coarse meshes.

Study of curved crack representation will be presented in a forthcoming paper, along with development of the proposed approach in three dimensions.

Fig. 1 .

 1 Fig. 1. Level set φ(x).

Fig. 3 .

 3 Fig. 3. Mesh of the boundary of a hole (grey): elements cut by the hole are refined internally once and twice (without adding new degrees of freedom), improving the geometrical representation of the hole.

Fig. 4 .

 4 Fig. 4. Numerical integration with X-FEM: creation of a partition.

Fig. 2 .

 2 Fig. 2. Boundary of a hole with finite elements.

  ðr; θÞ = a 8μ r a ðκ-3Þ sinðθÞ + 2 a r ð1-κÞ sinðθÞ + sinð3θÞ-2 a 3 r 3 sinð3θÞ " # ð6Þ with: • a = 0.4 • E = 10 5 : Young Modulus • ν = 0.3: Poisson's ratio

Fig. 5 .

 5 Fig. 5. Numerical integration with high order X-FEM: creation of a partition.

Fig. 11 .Fig. 12 .

 1112 Fig. 11. Initial enrichment.

Fig. 13 .

 13 Fig. 13. Iso-zero value of the level set.

Fig. 14 .

 14 Fig. 14. Energy error norm for the plate with a straight material interface.

Fig. 15 .

 15 Fig.15. Iso-zero value of the level set computed on the initial mesh and the octree (refined once).

Fig. 16 .

 16 Fig. 16. Construction of the enrichment function.

Fig. 17 .

 17 Fig. 17. Bubble function on an element: refined once and twice.

Fig. 18 .

 18 Fig. 18. Triangle element crossed by the iso-zero value of the level set.

Fig. 19 .

 19 Fig. 19. Barycentric coordinates in a triangle.

Fig. 20 .Fig. 21 .

 2021 Fig. 20. Different ways for an interface to cross a triangle element.

Fig. 22 .

 22 Fig. 22. Circular inclusionstrain with improvement of the geometry.

Fig. 23 .

 23 Fig.[START_REF] Osher | Level Set Methods and Dynamic Implicit Surfaces[END_REF]. Homogeneous refinement of the octree when an element is crossed by the material interface.

Fig. 24 .

 24 Fig.[START_REF] Pereira | hp-generalized fem and crack surface representation for non-planar 3-d cracks[END_REF]. Relative error norm with quadratic approximation shape functions with the proper enrichment function.

Fig. 25 .

 25 Fig. 25. Relative error norm with cubic approximation shape functions with the proper enrichment function.

Fig. 26 .

 26 Fig. 26. Plate with a squared inclusion under uniaxial tension.

Fig. 27 .

 27 Fig.[START_REF] Stazi | An extended finite element method with higher-order elements for curved cracks[END_REF]. FEM mesh with 16 elements by side (left), X-FEM mesh with 16 elements by side (right).

Fig. 28 .

 28 Fig. 28. Enrichment function for a squared material interface.

Fig. 29 .

 29 Fig. 29. Squared inclusion: energy error norm with linear approximation shape functions.

Fig. 30 .

 30 Fig.[START_REF] Szabó | Finite Element Analysis[END_REF]. Squared inclusion: energy error norm with quadratic approximation shape functions with the material interface enrichment function.

Fig. 31 .

 31 Fig. 31. Squared inclusion: energy error norm with cubic approximation shape functions with the material interface enrichment function.


Emilie Marchandise

Nicolas Chevaugeon

Jean-François Remacle
email: remacle@gce.ucl.ac.be

Spatial and spectral superconvergence of discontinuous Galerkin method for hyperbolic problems

Keywords: Discontinuous Galerkin method, Superconvergence, Hyperbolic systems

In this paper, we analyze the spatial and spectral superconvergence properties of one-dimensional hyperbolic conservation law by a discontinuous Galerkin (DG) method. The analyses combine classical mathematical arguments with MATLAB experiments. Some properties of the DG schemes are discovered using discrete Fourier analyses: superconvergence of the numerical wave numbers, Radau structure of the X spatial error.

Introduction

The discontinuous Galerkin method (DGM or DG method) has become a very popular numerical technique for solving hyperbolic conservation laws [START_REF] Flaherty | Aspects of discontinuous Galerkin methods for hyperbolic conservation laws[END_REF]. Some remarkable numerical properties of the DGM scheme have been demonstrated in an abundant literature. The method has been proven to be highly accurate, spectrally and spatially:

• Spectral superconvergence: Super-accuracy of numerical wave numbers was conjectured in [START_REF] Hu | Eigensolution analysis of the discontinuous Galerkin method with nonuniform grids, i, one space dimension[END_REF] and subsequently proved in [START_REF] Ainsworth | Dispersive and dissipative behavior of high order discontinuous Galerkin finite element methods[END_REF]. • Spatial superconvergence: Superconvergence at Radau points [START_REF] Cockburn | Enhanced accuracy by postprocessing for finite element methods for hyperbolic equations[END_REF][START_REF] Krivodonova | Error estimation for discontinuous Galerkin solutions of two-dimensional hyperbolic problems[END_REF][START_REF] Adjerid | A posteriori error estimation for discontinuous Galerkin solutions of hyperbolic problems[END_REF], strong superconvergence at downwind points [START_REF] Biswas | Parallel adaptive finite element method for conservation laws[END_REF] and Radau structure of the spatial error [START_REF] Krivodonova | Error estimation for discontinuous Galerkin solutions of two-dimensional hyperbolic problems[END_REF].

The aim of this work is to present those results in a comprehensive manner such that they become accessible to a wider community than the one of mathematicians. We propose here a unified way that allows to discover both spatial and spectral superconvergence properties of DGM's using simple MATLAB experiments. When necessary, we add to the discrete analysis some parts of mathematical proofs that highlight the discussion. In the following, we do not assume any specific choice for the numerical fluxes: upwind or not.

1

At the end, we try to give some reasonable answers to questions that naturally come out like:

• Is it useful to apply upwind numerical fluxes, why not using centered ones?

• On a DGM scheme that use polynomial order p on a mesh of size h, the L 2 norm of the error is of order h p+1 while the error on numerical wave numbers converges like h 2p+2 . What is the rate of convergence to consider? • What is the structure of the spatial error of the DGM scheme and are the jumps of the solution of any importance in the post processing (e.g. for error analysis)?

Discontinuous Galerkin space discretization in one dimension

Let us consider the following model problem. Being (x, t) ∈ (0, 1) × (0, T ), find u(x, t) solution of

j t u + cj x u = 0, (1)
with initial condition

u(x, 0) = u 0 (x), (2)
and periodic boundary conditions. This is a very simple model problem for an hyperbolic PDE. In (1), c is a positive number. Waves travel from left to right, in the positive x direction. We assume that time integration is exact and analyze the spatial and spectral discretization errors. We consider a partition of the space (0, 1) in N segments of size x = 1/N . In each element, we discretize in space the unknown using Legendre polynomials of orders less or equal to p. For that, we consider a reference segment ∈ (-1, 1) and the unknown u k in segment k going from

x k = (k -1)/N to x k+1 = k/N is approximated as u k (x( ), t) = p i=0 P i ( )u k i (t), with x( ) = x k 1 - 2 + x k+1 1 + 2 . (3)
Fig. 1 shows the one dimensional discretization.

The Legendre polynomials we are using here have the following properties. They are orthogonal in (-1, 1) i.e. the mass matrix

M pq = 1 -1 P p ( )P q ( ) d = 1 2p + 1 pq (4) is diagonal.
Legendre polynomials have known values at = ±1:

P p (-1) = (-1) p and P p (1) = 1. ( 5
)
This property allows us to state that

u k (x(-1), t) = p i=0 (-1) i u k i (t) and u k (x(1), t) = p i=0 u k i (t). (6)
Fig. 1. One-dimensional discretization.

2

For our further computations, the following integral is of importance

D pq = 1 -1 P p P q d .
First, we have that D pq = 0 if p q because Legendre polynomial of order p is orthogonal to any polynomial of lower order and the derivative of Legendre polynomial of order p is a polynomial of order p -1. The other case p > q is computed using integration by parts

D pq = 1 -1 P p P q d = - 1 -1 P p P q d + P p P q | 1 -1 = -D qp =0 +P p P q | 1 -1 = 1 -(-1) (p+q) . ( 7
)
The elements of D pq are then equal to 0 or 2, regardless of polynomial orders. We construct the discontinuous Galerkin formulation b(u k , v) of our problem, by multiplying (1) by a test function v(x), integrating the result over a segment k, and using the divergence theorem to obtain:

b(u k , v) = x k+1 x k j t u k (x, t)v(x) dx - x k+1 x k cu k (x, t)j x v(x) dx + [cu(x, t)v(x)] x=x k+1 x=x k = 0 ∀v. (8)
At the interface between two elements, i.e. the end points x k and x k+1 , the flux vector f = cu(x, t) is not uniquely defined and a flux formula has to be supplied to complete the discretization process. At any interface point x k , we define the following consistent approximation of the flux

ū(x k , t) = u k-1 (x k , t) + (1 -)u k (x k , t), (9)
where is the upwind parameter ( = 1 means a full upwind scheme and = 0.5 means a centered scheme). Eq. ( 8) for element k can be rewritten as:

b(u k , v) = x k+1 x k j t u k (x, t)v(x) dx - x k+1
x k cu k (x, t)j x v(x) dx

+ c ū(x k+1 , t)v(x k+1 ) -c ū(x k , t)v(x k ) = 0 ∀v. (10)
A second integration by parts of the second term of (10) leads to the standard Lesaint and Raviart [START_REF] Lesaint | Sur la résolution des systèmes hyperboliques du premier ordre par la méthode des éléments finis[END_REF] DG form:

b(u k , v) = x k+1 x k j t u k (x, t)v(x) dx + x k+1 x k cj x u k (x, t)v(x) dx + c[u k (x k , t) -u k-1 (x k , t)]v(x k ) + (1 -)c[u k+1 (x k+1 , t) -u k (x k+1 , t)]v(x k+1 ) = 0 ∀v. (11)
Next, if we consider expansion (3) for the numerical discretization of the unknown u k (x, t) and if we choose Legendre polynomials P j up to order j p as test functions v(x), the discrete variational formulation becomes:

b p i=0 u k i P i , P j = p i=0 j t u k i (t) x 1 -1 P i P j dx + cu k i (t) 1 - 1
P i P j dx + p i=0 c[P i (-1)u k i (t) -P i (1)u k-1 i (t)]P j (-1) + p i=0 (1 -)c[P i (-1)u k+1 i (t) -P i (1)u k i (t)]P j (1) = 0, j = 0, . . . , p. (12)
In ( 12) and for k = 1, the periodic boundary condition is imposed by stating that u 0 = u N . If we call

A ij = (-1) i+j , B ij = (-1) j and I ij = 1, then (12)
, i.e. the discrete discontinuous Galerkin formulation on segment k, can be written in matrix form as:

x c ju k i (t)M ij + (D ij + A ij -(1 -)I ij )u k i (t) -B ij u k-1 i (t) + (1 -)B ji u k+1 i (t) = 0. ( 13
)
If u is a column vector of size N × (p + 1) that contains all unknowns of every segment, the discontinuous Galerkin formulation can be written in a more compact form

j t u = - c x M -1 Lu,
with M = diag(M, . . . , M) a block diagonal matrix whose condition number grows like 2p + 1 and

L = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ (D T + A -(1 -)I ) (1 -)B • • • -B T -B T (D T + A -(1 -)I ) . . . . . . (1 -)B • • • -B T (D T + A -(1 -)I ) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦
, a sparse matrix.

Fourier analysis

Fourier analysis consider wave-like solutions

u(x, t) = Re(Ce 2i (kx-f t+i t) ),
where C is a complex constant, k is the wave number, f the frequency and is the damping parameter. The term "wave number" refers to the number of complete wave cycles that exist in one meter of linear space. Wave number k is dimensional: it is the inverse of a distance. The frequency f is dimensional too (inverse of a time). It is the number of complete wave cycles that are completed in one second. Substituting in (1), we have the dispersion relation and no damping:

f = ck and = 0.
The speed of the waves is c and the solution is neither amplified nor damped. If x is the mesh size and k is the dimensional wave number, we define a non-dimensional wave number as k h = k x. The non-dimensional wave number k h is interpreted as a wave number where the length measure is taken as the mesh size x. For example, a non-dimensional wave number of k h = 1 5 correspond to a wave length 5 x i.e. of five element sizes. We consider the semi-discrete form of the model problem (1):

j t u + Au = 0, with A = c x M -1 L. ( 14
)
We seek how accurately the DGM scheme, i.e. A, is able to approach the j x operator.

It is easy to compute the spectrum of A for different values of p, N and using MATLAB. We compute AV = VD where D is a diagonal matrix, D ii is the ith complex eigenvalue of A and where the ith column of V is the corresponding ith eigenvector of A.

Fig. 2 shows the distribution of the discrete eigenvalues of the DGM space operator. The mesh was made of N = 10 equally spaced segments. On the left Fig. 2 The centered scheme seems not to produce any numerical dissipation. This can be easily shown. In order to prove the L 2 stability of the DGM, we sum Eq. ( 11) over all elements and integrate it in time from t = 0 to T :

B(u k , u k ) = T 0 N k=1 b(u k , u k ) dt = 1 2 N k=1 x k+1 x k u k (x, T ) 2 -u k (x, 0) 2 dx - c 2 N k=1 T 0 u k (x k , t) 2 -u k-1 (x k , t) 2 dt + c N k=1 T 0 [ u k (x k , t) + (1 -)u k-1 (x k , t)][u k (x k , t) -u k-1 (x k , t)] dt = 1 2 u(T ) 2 L 2 (0,1) - 1 2 u(0) 2 L 2 (0,1) + c - 1 2 N k=1 T 0 (u k (x k , t) -u k (x k-1 , t)) 2 dt = 0. (15)
Clearly, the scheme is L 2 stable if and only if the sign of c is the same as the sign of -1 2 . Taking = 1 2 leads to a scheme that exactly conserves the L 2 norm of the solution (no numerical dissipation). The full upwind scheme consists in choosing = (1 + sign(c))/2. In Eq. ( 15), we observe that the dissipation mechanism involved in the DG scheme is produced by the solution jumps.

The eigenfunction u i (x) is constructed using the DGM interpolant:

u i (x) = Re ⎛ ⎝ N j =0 P j (x)V ji ⎞ ⎠ .
The eigenfunctions u i (x) are an approximation of the exact eigenfunctions U i (x) of the space operator. Those can be written in the following form:

u i (x) U i (x) = cos(2 k i x + i ),
where k i is the wave number of the mode. In order to find the corresponding wave number k i and the phase angle i , we first perform a fast Fourier transform (FFT) of each numerical mode u i (x). We use the Fourier spectrum to determine the wave number k i of the mode: it corresponds to the maximum in the energy spectrum. The corresponding phase angle i is also given by the FFT: it is the angle of the corresponding complex Fourier coefficient.

Fig. 3 shows some of the eigenvectors of the DGM discretization. The mesh is equally spaced with N = 5, p = 4 and = 1. We see that the eigenmodes corresponding to small wave numbers are well resolved, they are close to cosine functions. Higher order modes are unresolved and jumps of increasing size are appearing at inter-element boundaries. Those will be damped in time (see Eq. ( 15)). Actually, Ainsworth showed in [START_REF] Ainsworth | Dispersive and dissipative behavior of high order discontinuous Galerkin finite element methods[END_REF] that the number of resolved waves grows like 2p + 1. The resolved wave numbers are then such that: The resolved wave numbers corresponds to relatively small kh i.e. kh that are such that 2p + 1 > (2 k x) + O(k x)1/3 . 1 In our example, the resolved wave numbers are the ones such that k x < 1.43 -O(kh) 1/3 . Now what happens with the non-resolved waves? We show in the next example that according to the used spatial scheme, those non-resolved waves behave differently. We advect a step function located at x = 0.5 with a uniform velocity to the right without any limiter. The mesh is made of N = 50 equally spaced segments. The step contains thus a very large range of wave numbers. We see in Fig. 4 that the upwind scheme dissipate the non-resolved waves while the centered scheme keep them alive. Those non-dissipated resolved waves pollute the solution. Note that the centered scheme produces spurious undamped waves that travel in the wrong direction while the upwind scheme obviously only produce waves that go from left to right. The use of an upwind scheme does not avoid the presence of oscillations. The remaining overshoots that are present at high times, when all unresolved modes are damped, correspond to resolved modes that are present in the initial solution. Those modes are going at the right speed so that they do only perturb the solution at the vicinity of the discontinuity (the solution is perfectly constant away from the discontinuity). A limiter should be able to take care of those oscillations. We now look on how the DGM scheme is able to approximate eigenvalues and eigenfunctions of the space operator. Look at eigenfunctions will tell us about space accuracy of the DGM scheme. Looking at eigenvalues will tell us about spectral accuracy of the DGM scheme.

Spectral accuracy

The dispersion error usually defines how accurately the scheme is able to predict numerical wave speeds. In other words, the dispersion error of mode i is

E disp = 2 k i -Im(D ii ). ( 16
)
The dissipation error usually defines how accurately the scheme is able to predict numerical damping. The dissipation error is the error in the real part of the eigenvalues D ii :

E diss = Re(D ii ). ( 17
)
For the solution of the wave equation, the exact dissipation is null and all eigenvalues should all be on the imaginary axis.

It has been conjectured [START_REF] Hu | Eigensolution analysis of the discontinuous Galerkin method with nonuniform grids, i, one space dimension[END_REF] and subsequently proved [START_REF] Ainsworth | Dispersive and dissipative behavior of high order discontinuous Galerkin finite element methods[END_REF] that the numerical eigenvalues D ii are superconvergent: the dispersion error is accurate to (kh) 2p+3 and the dissipation error accurate to (kh) 2p+2 . This means that numerical wave numbers are converging more rapidly that the actual formal accuracy of the DGM (O(h) p+1 ). Fig. 5 illustrates the superconvergence of numerical wave numbers with N = 15, p = 4 and a full upwind scheme = 1.

Note that the same kind of plots can be done for the centered case. In this case, there is no dissipation error. Dispersion errors are also superconvergent in the centered case.

Spatial accuracy

In this section, we look at the spatial accuracy of the DGM scheme. We look on how the scheme is able to provide an accurate approximation of the real eigenfunctions U i (x). We define the pointwise spatial error of mode i The L 2 norm of the error for mode i (wave number k i ) is computed as . This is consistent with classical interpolation theory : approximating a function using polynomials of order p produces truncation error of order p + 1.

e i (x) = u i (x) -U i (x).
E 2 i = 1 0 e 2 i (x) dx.
In Fig. 7, we plot the function e i (x) for different k i 's in the case N = 5, p = 4 and for two values of . The existence of a structure for the DGM spatial error appears very clearly for = 1: error shapes are similar for all modes and for all elements. For other values of , the spatial structure disappears. In order to find out the analytical form of the structure of the DG spatial error, we build the error functional. Taking into account that the exact solution u(x, t) is smooth, we have for one element:

b(U i , P j ) = x k+1 x k (j t U i + cj x U i )v dx = 0 ∀j . (18)
Now, we assume that the solution is exact outside the element. Then, we can compute the local error functional in the upwind case. For doing so, we take the difference between Eq. (11) with = 1 and v(x) = P j and Eq. (18). We also assume that we have the exact solution at the downwind end, i.e. u k-1 i (x k ) U i (x k ) and get finally:

b(e k i , P j ) =

We further assume that the leading term of the local discretization error is a polynomial of order p + 1, i.e. e k i ( , t) = a k i,p+1 R p+1 ( ) + O( x) p+2 . The error Eq. ( 19) becomes, for the leading term:

1

-1 j R p+1 P j d + R p+1 P j | -1 = 0 ∀j . ( 20
)
It is easy to show, using [START_REF] Krivodonova | Error estimation for discontinuous Galerkin solutions of two-dimensional hyperbolic problems[END_REF], that R p+1 = P p+1 -P p verifies (20) and is equal to 0 at = 1. Moreover, it is possible to prove [START_REF] Adjerid | A posteriori error estimation for discontinuous Galerkin solutions of hyperbolic problems[END_REF] that the remaining terms of the error, up to 2p are also Radau polynomials. For example, p+3 . This means that the downwind points are superconvergent at order (k x) 2p+1 . A more general result is the following one: the flux value ū defined by Eq. ( 9) is always more accurate than left or right values. Let define, respectively, the upwind and the Riemann error for the ith eigenvector:

e i (x) -a i,p+1 R p+1 = a i,p+2 R p+2 + O(k x)
E upw i = N k=1 |u k i (x k ) -U i (x k )|, E rie i = N k=1 | ūi (x k ) -U i (x k )|.
Fig. 8 shows the behavior of these indicators with a centered and an upwind scheme. In the upwind case, downwind values (E rie ) are superconvergent at order (k x) 2p+1 and upwind values (E upw ) are converging at the same rate (k x) p+1 as the DGM spatial error. Therefore, jumps of the solution at nodes u k (x k , t)u k-1 (x k , t) are a good image of the local error that exists at upwind points:

J = N k=1 |u k i (x k ) -u k-1 i (x k )| = N k=1 |u k i (x k ) -U i (x k ) + U i (x k ) -u k-1 i (x k )| N k=1 |u k i (x k ) -U i (x k )| + |u k-1 i (x k ) -U i (x k )| = E upw + E dwn = O(k x) p+1 .
In the centered case, E rie is a measure of the error on solution averages at nodes. Clearly, averages are better than left or right values but are not converging better. In the centered scheme, jumps are not useful for computing local errors. This is absolutely consistent with the L 2 stability result (15). In the centered case, jumps of the solution are not the cause of numerical dissipation.

Conclusion

In this paper, we presented a simple approach to study discontinuous Galerkin schemes using simple MATLAB experiments. Both the error on numerical wave numbers and on numerical wave shapes were studied.

One of the aim of the paper was to justify the choice of the upwind scheme vs. the centered scheme. We showed that the upwind scheme was the best choice for two major reasons:

• upwind schemes allow to dissipate non-resolved modes while being super accurate on resolved ones; • upwind schemes lead to an interesting structure of the spatial error: jumps are an image of the local discretization error.

In some further work, we will focus on the extension of the technique to 2D problems and to systems.

  (a), we have chosen an upwind discretization of the fluxes and in the right Fig. 2(b) a centered scheme.

Fig. 2 .

 2 Fig. 2. Plot of the numerical eigenvalues for the DGM operator: (a) full upwind scheme ( = 1); and (b) centered scheme ( = 0.5).

Fig. 3 .

 3 Fig. 3. Representation of nine normed eigenvectors of A corresponding to different non-dimensional wave numbers k x.

Fig. 4 .

 4 Fig. 4. Advection of a step function with upwind scheme (top) and centered scheme (bottom) (N = 50).

Fig. 5 .Fig. 6 .

 56 Fig. 5. Error on the numerical values of the DGM operator in the full upwind case: (a) dispersion error: O(kh 2p+3 ); and (b) dissipation error: O(kh 2p+2 ).

Fig. 6

 6 Fig. 6 plots spatial errors E as a function of adimensional wave numbers k x. The convergence rates are O(k x) p+1 .This is consistent with classical interpolation theory : approximating a function using polynomials of order p produces truncation error of order p + 1.In Fig.7, we plot the function e i (x) for different k i 's in the case N = 5, p = 4 and for two values of . The existence of a structure for the DGM spatial error appears very clearly for = 1: error shapes are similar for all modes and for all elements. For other values of , the spatial structure disappears.

Fig. 7 .

 7 Fig. 7. Spatial error e(x) of the first eigenvectors corresponding to wave numbers 1/ x = 0.2 for different values of : (a) centered scheme ( = 0.5); and (b) upwind scheme ( = 1.0).

1 - 1 (

 11 xj t e k i + cj e k )P j d + ce k i (-1, t)P j (-1) = 0 ∀j .

Fig. 8 .

 8 Fig.8. Global upwind and Riemann error as a function of the non-dimensional wave number k x computed respectively with an upwind and a centered scheme.

In Ainsworth's paper, wave numbers are defined differently i.e. they refer to the number of cycles in

meters of linear space.


C Oliver-Leblond

F Ragueneau

A Delaplace

The Third International Conference on Computational Modeling of Fracture and Failure of Materials and Structures 3D Global/Local Analysis of Cracking of Reinforced Concrete

The assessment of crack features -pattern, spacing, opening -is a requirement for the design of reinforced concrete structures. Those assessments cannot rely on formula or experiment only, the first one being too simplistic and the second one too heavy. Therefore, numerical tools are developed to tackle the problem.

On one hand, non-linear finite element analyses based on continuum damage mechanics [START_REF] Mazars | Application de la mécanique de l'endommagement au comportement nonlinéaire et la rupture du béton de structure[END_REF][START_REF] Simo | Strain and stress based continuum damage models formulation[END_REF] are efficient to obtain the global behavior of large reinforced concrete structures. RicRag [START_REF] Richard | Isotropic continuum damage mechanics for concrete under cyclic loading: Stiffness recovery, inelastic strains and frictional sliding[END_REF] is a finite element model coupling elasticity, isotropic damage and internal sliding. It can describe the local mechanisms related to concrete such as the asymmetry between the tensile behavior and the compressive behavior, the inelastic strains and the unilateral effect. Moreover, it is robust and can handle large-scale computation. However, this model does not model the discontinuities and thus makes the study of the cracks complex.

On the other hand, discrete element models [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF][START_REF] D'addetta | On the application of a discrete model to the fracture process of cohesive granular materias[END_REF] are inherently capable of representing discontinuities such as cracks. DEAP [START_REF] Delaplace | Discrete 3D model as complimentary numerical testing for anisotropic damage[END_REF] is a lattice model for which the main physical mechanisms of quasi-brittle materials rupture are recovered, such as spatial correlation, crack tortuosity or scale effects. Unfortunately, the mesh density needed for such modeling is prohibitive to treat the case of industrial structures.

A global/local analysis has been proposed [START_REF] Oliver-Leblond | Non-intrusive global/local analysis for the study of fine cracking[END_REF] to take advantage of each model at their adequate scale. This non-intrusive technique allows the use of finite element models at a structural scale and a decoupled local analysis of some interesting areas, i.e. around cracks, for which a discrete element model is used.

The presentation will focus on the extension of the global/local method for 3D analysis of reinforced concrete beams. The treatment of reinforcement and interface will be addressed.

The numerical strategy is confronted to two different experiments to show its capabilities. The experimental crack pattern and crack openings are obtained on the surface with a digital image correlation analysis. The first experiment is a laboratory three points bending test performed on a concrete beam with light reinforcement. The second experiment is an industrial four points bending test on a massive concrete beam with heavy reinforcement.