Pierre Joyot
email: p.joyot@estia.frj.trunzler@estia.frfrancisco.chinesta@paris.ensam.fr

Jean Trunzler

Fransisco Chinesta

Enriched Reproducing Kernel Approximation: Reproducing Functions with Discontinuous Derivatives

In this paper we propose a new approximation technique within the context of meshless methods able to reproduce functions with discontinuous derivatives. This approach involves some concepts of the reproducing kernel particle method (RKPM), which have been extended in order to reproduce functions with discontinuous derivatives. This strategy will be referred as Enriched Reproducing Kernel Particle Approximation (E-RKPA). The accuracy of the proposed technique will be compared with standard RKP approximations (which only reproduces polynomials).

Introduction

Meshless methods are an appealing choice to develop functional approximations (with different degrees of consistency and continuity) without a mesh support. Thus, this kind of techniques seem to be specially appropriated for treating 3D problems involving large displacements, due to the fact that the approximation is constructed only from the cloud of nodes whose position evolve during the material deformation. In this manner neither remeshing nor fields projections are a priori required.

In order to adapt the approximation for introducing some known information associated to the searched solution different possibilities exist. The first one lies in the enrichment of standard approximations in the framework of the partition of unity (PU), originally proposed by Babuška and Melenk [START_REF] Babuška | The partition of unity method[END_REF], and which is at the origin of the partition of unity finite elements (PUFEM) [START_REF] Melenk | The partition of unity finite element method: Basic theory and applications[END_REF]. Generalized finite elements (G-FEM) [START_REF] Strouboulis | The generalized finite elment method[END_REF] or extended finite elements (X-FEM) [START_REF] Moës | Extended finite element method for cohesive crack growth[END_REF][START_REF] Sukumar | Modeling holes and inclusions by level sets in the extended finite-element method[END_REF][START_REF] Ventura | A vector level set method and new discontinuity approximations for crack growth by efg[END_REF][START_REF] Ji | A hybrid extended finite/level set method for modeling phase transformtions[END_REF] are two members of this family. In the extended finite elements the approximation is locally enhanced by introducing appropriate functions describing the known behavior of the problem solution.

The second possibility lies in the introduction of some information related to the problem solution within the approximation functional basis. For this purpose, different reproduction conditions are enforced in the construction of the approximation functions. This approach has been widely used in the context of the moving least squares approximations used in the diffuse meshless techniques [START_REF] Nayrolles | Generalizing the finite element method: Diffuse approximation and diffuse elements computational mechanics[END_REF] as well as in the element free Galerkin method [START_REF] Belytschko | Meshless methods: an overview and recent developments[END_REF]. Thus, all the known information related to the problem solution can be introduced in the functional approximation. Very accurate results were obtained for example in fracture mechanics by introducing the crack tip behavior into the approximation basis [START_REF] Fleming | Enriched element-free galerkin methods for crack tip fields[END_REF].

In this work we propose a numerical strategy, based on the reproducing kernel particle techniques, able to construct approximation functions with discontinuous derivatives on fixed or moving interfaces. This problem was treated in the context of the partition of unity by Kronggauz et al. [START_REF] Kronggauz | Efg approximation with discontinuous derivatives[END_REF]. In our approach the size of the discrete system of equations remains unchanged because no additional degrees of freedom are introduced related to the enrichment. However, the fact of enriching the approximation implies, as shown later, a bigger moment matrix, with worse properties, but in some cases this enrichment is only local, and in any case, the moment matrix has a low dimension.

The starting point of our development is the reproducing kernel particle approximation (RKPA). The RKP approximation was introduced by Liu et al. [START_REF] Liu | Reproducing kernel particle methods for structural dynamics[END_REF] for enforcing some degree of consistency to standard smooth particle approximations, i.e. they proved that starting from a SPH (smooth particle hydrodynamics) approximation [START_REF] Gingold | Smoothed particle hydrodynamics: Theory and application to non-spherical stars[END_REF] it is possible to enhance the kernel function for reproducing a certain degree of polynomials. We have extended or generalized this procedure in order to reproduce any function, and more concretely, functions involving discontinuous derivatives.

Moreover, the strategy here described, can be easily coupled with a level set description of the interface location [START_REF] Sethian | Evolution, implementation and application of level set and fast marching methods for advancing fronts[END_REF] (recently introduced in the framework of the extended finite elements [START_REF] Ventura | A vector level set method and new discontinuity approximations for crack growth by efg[END_REF]), which allows to capture its position when it evolves in time. Thus, both the enrichment and the interface tracking are made in a transparent way for the user.

This paper consists of 6 sections, in which we summarize the reproducing kernel approximation (section 2) which allows to treat in section 3 the reproduction of a function with discontinuous derivatives. In sections 4 and 5 we describe and illustrate the properties of the resulting approximation. Finally, in section 6 we focus on a test problem which allows to conclude about the accuracy and potentiality of the proposed strategy.

Enriched Reproducing Kernel Particle Approximation

Let Ω be the domain where the problem is defined. The points within this domain will be noted by x or s. For the sake of simplicity, from now on, we will focus in the 1D case, but all the results have a direct 2D or 3D counterpart.

Reproduction Conditions

The approximation f a (x) of f (x) is built from the integral convolution integral

f a (x) = Ω Φ(x -s, h)f (s)dΩ, (2.1)
where Φ(xs, h) is the kernel function and h a parameter defining the size of the approximation support.

The main aim of this work is to enforce the reproduction of a general function that we can write in the form of a polynomial plus another function noted by f e (x):

f a (x) = a 0 + a 1 x + . . . + a n x n + f e (x). (2.2)
In the following paragraphs we analyze the required properties of the kernel function Φ(xs, h) for reproducing a function expressed by (2.2). From Eq. (2.1), the reproduction of a constant function a 0 is given by

Ω Φ(x -s, h)a 0 dΩ = a 0 , (2.3)
which implies

Ω Φ(x -s, h)dΩ = 1, (2.4)
which constitutes the partition of unity property. Now, the required condition to reproduce a linear function

f a (x) = a 0 + a 1 x is Ω Φ(x -s, h)(a 0 + a 1 s)dΩ = a 0 + a 1 x. (2.5)
By using the partition of unity (2.4), Eq. (2.5) can be rewritten as

Ω Φ(x -s, h)dΩ = 1 Ω Φ(x -s, h)sdΩ = x , ( 2.6)
which implies the linear consistency of the approximation. Repeating this reasoning, we can write the n-order consistency as

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ Ω Φ(x -s, h)dΩ = 1 Ω Φ(x -s, h)sdΩ = x . . . Ω Φ(x -s, h)s n dΩ = x n (2.7)
and consequently, the reproduction of the function given by (2.2) implies

Ω Φ(x -s, h)(a 0 + a 1 s + . . . + a n s n + f e (s))dΩ = = a 0 + a 1 x + . . . + a n x n + f e (x), (2.8) from which it results ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ Ω Φ(x -s, h)dΩ = 1 Ω Φ(x -s, h)sdΩ = x . . . Ω Φ(x -s, h)s n dΩ = x n Ω Φ(x -s, h)f e (s)dΩ = f e (x)
.

(2.9)

Chen et al. [START_REF] Chen | A lagrangian reproducing kernel particle method for metal forming[END_REF] applied the procedure proposed by Liu et al. [START_REF] Liu | Reproducing kernel particle methods for structural dynamics[END_REF] for enforcing n-order consistency, however their procedure does not allow to enforce directly the reproduction condition associated with f e (x).

The Moment Matrix

We will note by f r (x) the approximation function verifying the conditions (2.9). Usually a cubic spline is considered as kernel function, and consequently the conditions given by Eq. (2.9) are not satisfied. Liu et al. [START_REF] Liu | Reproducing kernel particle methods for structural dynamics[END_REF] propose the introduction of a correction function C(x, x-s) for satisfying the reproduction conditions. In our case we consider the more general form C(x, s, xs) whose pertinence will be discussed later. Thus f r (x) will be expressed by

f r (x) = Ω C(x, s, x -s)Φ(x -s, h)f (s)dΩ, (2.10)
where C(x, s, xs) is assumed to have the following form

C(x, s, x -s) = H T (x, s, x -s)b(x), (2.11)
where H T (x, s, x-s) represents the vector containing the functions considered in the approximation basis, and b(x) is a vector containing unknown functions that will be determined for satisfying the reproduction conditions. Thus, Eq.

(2.9) can be rewritten as

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ Ω H T (x, s, x -s)b(x)Φ(x -s, h)dΩ = 1 Ω H T (x, s, x -s)b(x)Φ(x -s, h)sdΩ = x . . . Ω H T (x, s, x -s)b(x)Φ(x -s, h)s n dΩ = x n Ω H T (x, s, x -s)b(x)Φ(x -s, h)f e (s)dΩ = f e (x)
.

(2.12)

Remark: Usual RKP approximations consider H T (xs), but we retain the more general form H T (x, s, x-s) in order to reproduce more general functions.

In fact, the reproduction conditions must be enforced in a discrete form. For this purpose we consider N P points (also refereed as nodes) which allow to compute the discrete form of Eq. (2.12), i.e.

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ NP I=1 H T (x, x I , x -x I )b(x)Φ(x -x I , h)Δx I = 1 NP I=1 H T (x, x I , x -x I )b(x)Φ(x -x I , h)x I Δx I = x . . . NP I=1 H T (x, x I , x -x I )b(x)Φ(x -x I , h)x n I Δx I = x n NP I=1 H T (x, x I , x -x I )b(x)Φ(x -x I , h)f e (x I )Δx I = f e (x)
, (2.13) that in a matrix form results

NP I=1 R(x I )H T (x, x I , x -x I )Φ(x -x I , h)Δx I b(x) = R(x), (2.14)
where R(x) is the reproduction vector

R T (x) = [1, x, . . . , x n , f e (x)] . ( 2
(x) = M (x) -1 R(x), (2.16)
where the moment matrix M (x) is defined by

M (x) = NP I=1 R(x I )H T (x, x I , x -x I )Φ(x -x I , h)Δx I .
(2.17)

This moment matrix differs of the usual moment matrix proposed in [START_REF] Liu | Reproducing kernel particle methods for structural dynamics[END_REF], and in fact it becomes non symmetric.

Discrete Form of the Approximation Function

The discrete form f r (x) of f a (x) derives from Eqs. (2.10), (2.11) and (2.16)

f r (x) ∼ = NP I=1 H T (x, x I , x -x I )M (x) -1 R(x)Φ(x -x I , h)f (x I )Δx I = NP I=1 ψ I (x)f I , (2.18)
where ψ I is the enriched RKP approximation shape function

ψ I (x) = H T (x, x I , x -x I )M (x) -1 R(x)Φ(x -x I , h)Δx I . (2.19)
As in the classical RKPM we take Δx I = 1. Different quadrature rules exist and they have been tested in [START_REF] Aluru | A point collocation method based on reproducing kernel approximations[END_REF] without a significant incidence on the reproducing condition accuracy.

The Case of a Function with Discontinuous Derivatives

Let Ω be the domain where the problem is defined and Γ d a point, curve or surface (in 1D, 2D and 3D respectively) where the normal function derivative becomes discontinuous. We assume that this discontinuity curve splits the domain in two subdomains Ω 0 and Ω 1 (see Figure 3

Ω 0 ∪ Ω 1 ∪ Γ d = Ω, Ω 0 ∩ Ω 1 = ∅,
where Γ d represents the curve along which the normal derivative of some approximation functions become discontinuous.

The enrichment function f e (x) which will be introduced in the reproduction vector R(x) must satisfy the following conditions

f e (x) ∈ C 0 (x ∈ Ω), f e (x) ∈ C 1 (x ∈ Ω 0 ), f e (x) ∈ C 1 (x ∈ Ω 1 ), Grad f e (x) • n ∈ C -1 (x ∈ Γ d ), (3.1)
where n denotes the unit outward vector defined on the curve Γ d .

To locate the discontinuity curve Γ d we have use of a level set function Θ(x) defined as the signed distance from x to the interface Γ d . Thus

Θ(x) = ⎧ ⎨ ⎩ Θ(x) < 0 if x ∈ Ω 0 Θ(x) > 0 if x ∈ Ω 1 Θ(x) = 0 if x ∈ Γ d . ( 3.2)
Thus, for satisfying the conditions (3.1) the enrichment function f e (x) could be assumed in the form

f e (x) = H 0 (Θ(x))Θ(x), (3.3)
where H 0 (x) represents the usual Heaviside function such that now the domain Ω consists of four subdomains Ω c0 , Ω e0 , Ω c1 and Ω e1 , where classical and enriched reproducing kernel approximation will be defined with continuity along Γ 0 and Γ 1 as proved later. Thus, we start from the relations

H 0 (Θ(x)) = 1 if Θ(x) ≥ 0 H 0 (Θ(x)) = 0 if Θ(x) < 0 . (3.4) Ω c0 Ω c1 Γ d Γ 0 Γ 1 Ω e0 Ω e1
Ω 0 = Ω c0 ∪ Ω e0 ∪ Γ 0 , Ω 1 = Ω c1 ∪ Ω e1 ∪ Γ 1 .
Obviously, Γ 0 and Γ 1 define the points, curves or surfaces (1D, 2D and 3D respectively) in Ω 0 and Ω 1 respectively, through which the determinant of the moment matrix changes from zero to a non-null value.

We define R T c (x) = [1, x, . . . , x n ], R T (x) = R T c (x), H 0 (Θ(x))Θ(x) , and in a similar manner H T c (x-x I ) = [1, xx I , . . . , (xx I ) n ], and H T (x, x I , x-

x I ) = H T c (xx I ), H 0 (Θ(x I ))Θ(x I ) . The following properties of the moment matrix can be deduced: Property 1: The inverse of the moment matrix exists in the domains Ω e0 and Ω e1 , i.e. in the subdomains whose boundary involves the discontinuity curve Γ d . Then, the moment matrix results

M (x) = NP I=1 Rc(xI )H T c (x-xI ) Rc(xI )H0(Θ(xI ))Θ(xI ) H0(Θ(xI ))Θ(xI )H T c (x-xI ) H0(Θ(xI )) 2 Θ(xI ) 2 Φ(x-xI ,h)ΔxI .
Property 2: In the domain Ω c0 , H 0 (Θ(x I )) = 0, ∀I. Thus, the moment matrix becomes singular

M (x) = NP I=1 R c (x I )H T c (x -x I ) 0 0 0 Φ(x -x I , h)Δx I .
In this case we propose to replace this moment matrix for the one derived from a standard (non-enriched) approximation

M c (x) = NP I=1 R c (x)H T c (x -x I )Φ(x -x I , h)Δx I . (4.1)
Property 3: Within the domain Ω c1 , H 0 (Θ(x I )) = 1, ∀I. Thus, M (x) results

M (x) = NP I=1 R c (x I )H T c (x -x I ) R c (x I )Θ(x I ) Θ(x I )H T c (x -x I ) Θ(x I ) 2 Φ(x -x I , h)Δx I .
The moment matrix is singular in 1D (as will be illustrated in section 5.1) as well as in some 2D and 3D cases. In this case we replace again the enriched moment matrix for the non-enriched one. Effectively, if the enriched moment matrix is not singular everywhere in Ω c1 , then Γ 1 does not exist, i.e.

Ω c1 = ∅ et Ω e1 = Ω 1 .
Property 4: When we approach the curve Γ 0 the shape functions computed from both the enriched and the non-enriched approximations are connected with continuity, i.e. lim

x∈Ωc0→Γ0 ψ I (x) = lim x∈Ωe0→Γ0 ψ I (x), ∀I.
And in a similar way when we approach the curve Γ 1 (when it exists)

lim x∈Ωc1→Γ1 ψ I (x) = lim x∈Ωe1→Γ1 ψ I (x), ∀I.
This property has important consequences because due to the continuity of the shape functions within the whole domain, the numerical integration of variational formulations does not require a specific treatment. The following results can be stated:

The curves Γ 0 and Γ 1 belong to the union of the support boundary ζ I of some nodes. In the 1D case Γ 0 and Γ 1 reduce to two points (as illustrated in section 5). In higher dimensions, both boundaries are composed by a series of arcs of circle (sphere in 3D) (when the support of the kernel function is assumed circular -spherical in 3D-).

About the Resulting Shape Function and its Derivatives

In this section we focus in a one-dimensional problem in order to illustrate some of the previous properties. The signed distance is used for defining the level set function

Θ(x) = x -x d , (5.1)
where x d represents the coordinate of the point where the discontinuity in the function derivative takes place.

Moment Matrix

For the sake of clarity, from now on, we restrict the consistency requirement to a first order. Thus, the moment matrix results

M (x) = NP I=1 M I (x)Φ(x -x I , h) with M I (x) = ⎡ ⎣ 1 x-xI H0(xI -x d )(xI -x d ) xI xI (x-xI ) xI H0(xI -x d )(xI -x d ) H0(xI -x d )(xI -x d ) H0(xI -x d )(xI -x d )(x-xI ) (H0(xI -x d )(xI -x d )) 2 ⎤ ⎦ .
For illustrating the third property, we consider x ∈ Ω c1 , which leads to the following expression of the moment matrix

M (x) = NP I=1 ⎡ ⎣ 1 x-xI (xI -x d ) xI xI (x-xI ) xI (xI -x d ) (xI -x d ) (xI -x d )(x-xI ) (xI -x d ) 2 ⎤ ⎦ Φ(x-xI ,h).
We can notice that the third row is a linear combination of the two first. Thus, the determinant vanishes and the moment matrix becomes singular. This result does not depend on the degree of consistency assumed. Now we are going to illustrate the position of the boundaries Γ 0 and Γ 1 through which the approximation changes from enriched to non-enriched.

In Figure 5.3 the vertical line indicates the position where the function derivative is expected to be discontinuous, and the horizontal one represents the enrichment domain Ω e1 ∪Ω e0 . The distance between two consecutive nodes (assumed regularly spaced) is h = 1 and the support radius is taken in the form R = a × h (where the value of a is indicated in each subfigure).

It is well known that for a standard RKP approximation (without any enrichment), the linear consistency implies that each point must be contained by more than one nodal support (in other case the moment matrix becomes singular) [START_REF] Han | Error anlysis of the reproducing kernel particle method[END_REF]. In our case we need to add to this constraint the properties listed in the previous section.

From Figure 5.3(b) we can notice:

• Point A is in the domain Ω c0 , i.e. A ∈ Ω c0 because any node whose support contains the point A is not in Ω 1 (Property 2 ). • D ∈ Ω c1 because any node whose support contains the point D is not in Ω 0 (Property 3 ). similar results when N P = 9 and consequently a node is placed just at the discontinuity point.

Remarks:

• We have noticed that the shape functions are perfectly continuous through the interfaces between the enriched and non-enriched domains (Property 4 ). Moreover the shape function derivatives become discontinuous only along the interface x d = 0. • The shape functions are quite different when a node is or not located at the discontinuity point. In the first case (node located at x d ) the shape functions are very close to those obtained by Krongauz et al. [START_REF] Kronggauz | Efg approximation with discontinuous derivatives[END_REF] using a PU enrichment.

Numerical Example

In this section we consider a simple numerical test: the steady heat conduction problem (with a constant source term g = 2) defined in the interval ]-1, 1[ which involves two homogeneous materials with different thermal conductivities

d dx k(x) dT dx = g, (6.1)
where T is the temperature field and k the thermal conductivity defined as follows

k(x) = k 0 = 1 if x ∈ Ω 0 = ]-1, 0[ , k(x) = k 1 = 10 if x ∈ Ω 1 = ]0, 1[ . ( 6.2)
The boundary conditions are defined by

T (x = -1) = 0 T (x = 1) = 1 . (6.3)
The problem is solved using the weak-form of Eq. (6.1), where the integrals were computed using a background mesh. Four integration cells are placed between two consecutive nodes, with 5 gauss integration points in each cell.

Computed temperatures with and without approximation enrichment as well as the exact temperature field are depicted in figure 6.7. In both cases (enriched and non-enriched) a first order consistency has been enforced, been the support radius R = 2× h. The related derivatives are shown in Figure 6.8. From both figures we can conclude about the better accuracy of the solution computed by using an enriched approximation, mainly in the temperature derivative in the discontinuity neighbourhood.

Convergence Analysis

In this section we compare the convergence of both enriched and non-enriched approximations, in the test problem just described. Two norms are considered, the L 2 and the H 1 norms. Moreover, the analysis has been performed for three different conductivity ratios k 01 = k1 k0 = 10, 100, 1000. Figure 6.9 illustrates these results.

In the case of non-enriched approximations (standard RKPA) we can notice from figure 6.9 that the order of the method is 1 and 0.5 (using the L 2 and H 1 norms respectively). Moreover, the error increases as the conductivity ratio increases.

On the contrary, when an enriched approximation is considered the method results of order 2.4 and 1.5 (using the L 2 and H 1 norms respectively). Moreover, the error is independent of the conductivity ratio.

In conclusion, both the order of convergence and the solution error obtained by using an enriched approximation are significantly better that the ones obtained by using the standard non-enriched formulation.

Conclusion

In this work we have proposed a new strategy to enrich reproducing kernel particle approximations in order to reproduce functions involving discontinuous derivatives through some interfaces. The computed results seems to be very accurate (from the point of view of the order of convergence as well as from the associated level of error) and it is very easy to implement. In some cases the moment matrix becomes singular far from the discontinuity location (in the order of the support size), and consequently the enrichment becomes local in space. In any case, the size of the global discrete system is not affected by the enrichment process, in contrast with the vast majority of enrichment methods making use of the partition of unity.

The case of moving interfaces involving field discontinuities can be easily treated coupling the technique described in this work with a level set description of the interface movement.

The extension to higher dimension problems involving complex discontinuity geometries is a work in progress.

. 15 )

 15 Eq. (2.14) allows the computation of vector b(x), b

Figure 3 . 1 .

 31 Figure 3.1. Problem domain containing an interface with a discontinuous normal derivative.

Figure 4 . 2 .

 42 Figure 4.2. Enriched and non-enriched domains.

Property 5 :

 5 Let ζ I be the boundary of the support related to the node x I .

Property 6 :

 6 The boundary ζ I associated with the node x I ∈ Ω 1 (respectively x I ∈ Ω 0 ) closest to the discontinuity interface Γ d defines the boundary Γ 0 (respectively Γ 1 ).

Figure 5 . 4 .Figure 5 . 5 .

 5455 Figure 5.4. Shape functions for a discontinuous derivative located at x d = 0.

11 -Figure 5 . 6 .

 1156 Figure 5.6. Shape functions for a discontinuous derivative located at x d = 0 when a node is placed just at the discontinuity point.

Figure 6 . 7 .Figure 6 . 8 .

 6768 Figure 6.7. Exact and numerical temperature fields computed with and without approximation enrichment.

Figure 6 . 9 .

 69 Figure 6.9. Convergence analysis using the L2 and the H1 norms.

Properties of the Moment MatrixIf the moment matrix is computed as described in section 2 everywhere, it can become singular, as we will prove later. We can define two boundaries in the discontinuity curve neighborhood Γ 0 and Γ 1 (as illustrated in figure4.2)

• B ∈ Ω e0 and C ∈ Ω e1 (Property 1 ).

• The position of Γ 1 (respectively Γ 0 ) is defined by the support of the node within Ω 0 (respectively Ω 1 ) closest to the discontinuity interface Γ d (Properties 5-6 ).

If a < 1.5 there are points, candidates to belong to Ω e0 ∪ Ω e1 , such that they belong only to the support of two nodes. In this way the approximation is only defined if a > 1.5 (Figure 5.3(a)).

When there is a node located at the discontinuity point this problem is avoided as depicted in Figure 5.3(c). Thus, the approximation is defined only if a > 1, like in the standard (non-enriched) RKP approximation.

Related Shape Functions


Julien Yvonnet
email: julien.yvonnet@paris.ensam.fr

David Ryckelynck

Philippe Lorong

Francisco Chinesta

Treating Moving Interfaces in

published or not. The documents may come

Introduction

Phase boundaries represent material interfaces across several fields which may exhibit sharp gradients, and even discontinuities. A wide range of numerical methods have been developed for treating these problems according to the pertinent physics and assumptions about the interface [START_REF] Lewis | Finite element simulation of metal casting[END_REF]. When a sharp interface is considered, its motion is governed by the jump in the temperature gradient normal to the phase boundary and is accompanied by the latent heat effects (Stefan condition). In order to satisfy these conditions the most common approach lies in tracking explicitly the interface motion. Within the interface tracking approach two main alternatives exist: the moving mesh methods and the mixed Eulerian-Lagrangian methods. Moving finite element mesh methods conform element boundaries to the interfaces as it evolves. Although these methods are very accurate, they are limited by the severe mesh distortion. Thus, frequent remeshing is needed, with the associated field projections between successive meshes. Moreover, remeshing is, even today, a delicate task in the 3D case. To alleviate remeshing efforts a number of Eulerian-Lagrangian methods have been developed recently that track the interface while solving the equations on a fixed grid [START_REF] Udaykumar | Computation of solid-liquid phase fronts in the sharp interface limit on fixed grids[END_REF]. Many of these methods effectively smear the discontinuity over a few grid cells, and are therefore not capable of representing the true discontinuity across the interface.

A new approach for representing localized behaviours has recently emerged in the field of the finite element method, known as the partition of unity method [START_REF] Melenk | The partition of unity finite element method: basic theory and applications[END_REF]. The main idea is to extend the classical approximation considering the product of the standard shape functions and local enrichment functions. The extended finite element method (X-FEM) is a variation on this framework. Recently, the X-FEM has been coupled to the Level Set Method [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flows[END_REF] to represent interface topologies [START_REF] Sukumar | Modeling holes and inclusions by level sets in the extended finte element method[END_REF]. In this way, the discontinuity evolution can be properly represented on a fixed background mesh, just by adding an appropriate enrichment in the functional approximation in the elements concerned by the moving discontinuity [START_REF] Ji | A hybrid finite element/level set method for modeling phase transformations[END_REF]. However, when the material in which the interface is moving, is subjected to large displacements, an updated Lagrangian description could be a better choice. When the background mesh evolves, remeshing will be required also to avoid too high background mesh distortions.

To alleviate dependance to the mesh and to provide smoother shape functions, the use of the meshfree or meshless methods is investigated. The meshless methods discretize a continuum body by a finite number of particles (or nodes) and the field of interest is interpolated under these nodes without the aid of an explicit mesh. Many meshless methods have been proposed, including the smooth particle hydrodynamics (SPH) [START_REF] Lucy | A numerical approach to the testing of fusion process[END_REF], the diffuse element method (DEM) [START_REF] Nayroles | Generalizing the finite element method: diffuse approximation and diffuse elements[END_REF], the element free Galerkin (EFG) [START_REF] Belytschko | Element-free Galerkin methods[END_REF], the reproducing kernel particle method (RKPM) [START_REF] Liu | Reproducing Kernel Particle Methods[END_REF], the HP clouds DOD96 and the partition of unity method (PUM) [START_REF] Melenk | The partition of unity finite element method: basic theory and applications[END_REF].

The introduction of moving discontinuities in these meshless methods can present difficulties for the following reasons: (i) The quality of the approximation as well as the conditioning of the global system is pathologically dependant from the size of the support of the shape functions; (ii) Imposition of essential boundary conditions needs particular treatments; (iii) Integration is not accurately defined, and (iv) The physical discontinuity across the interfaces must be introduced accurately. In order to overcome these different problems, we propose the use of the constrained natural element method (C-NEM) [YRa03] [YRb03] for treating thermal models involving moving interfaces. This approach is an extension of the natural element method [START_REF] Sukumar | The natural elements method in solid mechanics[END_REF] in which both trial and test functions are constructed on the basis of the Voronoi based interpolants [START_REF] Sibson | A vector Identity for the Dirichlet tesselations[END_REF] [HSU02], verifying the Kronecker delta property and whose support is defined by the union of the Delaunay spheres passing through the visible nodes. In the C-NEM, the introduction of a visibility criterion and its related constrained Voronoi diagram preserves the appealing properties of the NEM in any geometry (convex or not) and allows easily the introduction of material discontinuities. In a former paper the C-NEM has been sucessfully applied in some problems involving non convex domains, discontinuities and cracks [START_REF] Yvonnet | A new extension of the natural element method for non convex and dscontinuous problems, the constrained natural element method (C-NEM)[END_REF].

The structure of this work is as follow: in section 2, we introduce a simple mathematical model of a thermal problem involving a moving interface. Actually, it turns out that it results the standard Stefan problem [START_REF] Vuik | Some historical notes about the Stefan problem[END_REF]. In section 3, the C-NEM is summarized, and it will be applied in section 4 to discretize the weak formulation of the non-linear Stefan problem. A numerical benchmark is presented in section 5, which allows us to conclude about the accuracy of the proposed technique by comparing the numerical results with the exact solution of the problem.

Problem Formulation

Let Ω ⊂ 2 be a bounded domain and T the temperature field. On the domain boundary Γ ≡ ∂Ω the temperature or the thermal flux are prescribed. We will denote by Γ 1 the domain boundary where the temperature is known T (x ∈ Γ 1 , t) = T (x, t) and by Γ 2 the domain boundary where the heat flux q is imposed. The thermal model is defined in the time interval [0, t max ]. The initial temperature T (x, t = 0) = T 0 , where T 0 is assumed to be higher than the material melting temperature T m . At time, t = 0, a part of the domain boundary Γ 1 is suddenly submitted to a temperature T 1 < T m . A moving solidification front Γ I is then generated, whose position evolves in time, i.e. Γ I (t), dividing the domain Ω in two regions Ω 1 (t) (containing the solid phase at time t) and Ω 2 (t) (which contains the liquid phase) as shown in Figure 2.1. For a sake of simplicity we will consider, from now on, a homogeneous and isotropic thermal model in both phases. The heat transfer model is defined in each phase, neglecting volumetric source terms, by:

c 1 ∂T (x,t) ∂t = ∇ • (k 1 ∇T ) in Ω 1 (t) c 2 ∂T (x,t) ∂t = ∇ • (k 2 ∇T ) in Ω 2 (t) (2.1)
where c 1 and c 2 are the volumetric heat capacities of both phases, being k 1 and k 2 their thermal conductivities. The associated initial and boundary conditions are:

⎧ ⎨ ⎩ T (x, t = 0) = T 0 ∀x ∈ Ω T (x, t) = T (x, t) ∀x ∈ Γ 1 , ∀t ∈ [0, t max ] -k∇T (x, t) • n = q(x, t) ∀x ∈ Γ 2 , ∀t ∈ [0, t max ] (2.2)
The evolution of the interface Γ I (t) is described by a Stefan condition:

V(x ∈ Γ I (t)) = |[q]| L n I (x) (2.3)
where V is the interface velocity, L is the volumetric latent heat of fusion, n I (x) is the normal vector to the interface at point x which is assumed to point into the liquid phase, and [[q]] the thermal flux jump across the interface Γ I (t), i.e.

|[q]| = k 1 ∇T | Γ - I (t) -k 2 ∇T | Γ + I (t) n I (2.4)
The additional constraint prescribed on the interface Γ I (t) is:

T (x, t) = T m ; ∀x ∈ Γ I (t) (2.5)
where T m is the melting temperature.

The Constrained Natural Element Method (C-NEM)

In this section, the utility of the C-NEM to describe moving interfaces and discontinuities in a fixed cloud of nodes is discussed. After a brief review of the Voronoi-based interpolants, we introduce the constrained Voronoi diagram which is used for computing the shape functions in any domain.

Natural Neighbor Interpolation

We briefly touch upon the foundation of Sibson's natural neighbor coordinates (shape functions) that are used in the natural element method. For a more in-depth discussion on the Sibson interpolant and its application for solving second-order partial differential equations, the interested reader can refer to Braun and Sambridge [START_REF] Sambridge | Geophisical parameterization and interpolation of irregular data using natural neighbors[END_REF], and Sukmar et al. [START_REF] Sukumar | The natural elements method in solid mechanics[END_REF]. The NEM interpolant is constructed on the underlying Voronoi diagram. The Delaunay tesselation is the topological dual of the Voronoi diagram. Consider a set of nodes S = {n 1 , n 2 , . . . , n N } in 2 . The Voronoi diagram is the subdivision of 2 into regions T i (Voronoi cells) defined by:

T i = {x ∈ 2 : d(x, x i ) < d(x, x j ), ∀j = i}, ∀ i (3.1)
The Sibson coordinates of x with respect to a natural neighbor n i (see Figure 3.2) is defined as the ratio of the overlap area (volume in 3D) of their Voronoi cells to the total area (volume in 3D) of the Voronoi cell related to point x:

φ i (x) = Area(af ghe) Area(abcde) (3.2)
If the point x coincides with the node n i , i.e. (x = x i ), φ i (x i ) = 1, and all other shape functions are zero, i.e. φ j (x i ) = δ ij (δ ij being the Kroenecker's delta). The properties of positivity, interpolation, and partition of unity are then verified [START_REF] Sukumar | The natural elements method in solid mechanics[END_REF]:

⎧ ⎨ ⎩ 0 ≤ φ i (x) ≤ 1 φ i (x j ) = δ ij n i=1 φ i (x) = 1 (3.3)
The natural neighbor shape functions also satisfy the local coordinate property [START_REF] Sibson | A vector Identity for the Dirichlet tesselations[END_REF], namely:

x = n i=1 φ i (x)x i (3.4)
which combined with Eq. (3.3), implies that the natural neighbor interpolant spans the space of linear polynomials (linear completeness).

Natural neighbor shape functions are C ∞ at any point except at the nodes, where they are only C 0 , and on the boundary of the Delaunay circles (spheres in 3D) where they are only C 1 , because of the discontinuity in the neighbors nodes across these boundaries.

Another important property of this interpolant is the ability to reproduce linear functions over the boundary of convex domains. The proof can be found in Sukumar et al. [START_REF] Sukumar | The natural elements method in solid mechanics[END_REF]. An illustration is depicted in Figure 3.2 (b): as the areas associated to points on the boundary become infinite, the contribution of internal points vanish in the limit when the point approaches the convex boundary, and the shape functions associated with nodes n 1 and n 2 become linear on the segment (n 1n 2 ). This is not true in the case of non convex boundaries, and the next section focuses in an approach to circumvert this difficulty.

Consider an interpolation scheme for a scalar function T (x) : Ω ⊂ 2 → 2 , in the form:

T h (x) = n i=1 φ i (x) T i (3.5)
where T i are the nodal temperatures at the n natural neighbor nodes, and φ i (x) are the shape functions associated with each neighbor node. It is noted that Eq. (3.5) defines a local interpolation scheme. Thus, the trial and test functions used in the discretization of the variational formulation describing the thermal problem treated in this paper take the form of Eq. (3.5).

The Constrained Natural Element Method

Constrained Voronoi Diagram. It was proved in [START_REF] Yvonnet | A new extension of the natural element method for non convex and dscontinuous problems, the constrained natural element method (C-NEM)[END_REF] [SMB98] and [START_REF] Cueto | On the imposition of essential boundary conditions in Natural neighbor Galerkin Methods[END_REF] that spurious influences between "non-visible" nodes and lost of linearity in the interpolation along boundaries of non convex domains appear in the framework of the NEM. In order to avoid this drawback and to recover all properties of the method for any geometry (including non convex domains containing cracks or involving field discontinuities), a visibility criterion is introduced in order to restrict influent nodes among natural neighbors. The computation of the shape functions is done on the basis of the so-called constrained (or extended) Voronoi diagram (CVD), introduced by Seidel in [START_REF] Seidel | Constrained Delaunay triangulations and Voronoi diagrams with obstacles In: "1978-1988 Ten Years IIG[END_REF]. The constrained Voronoi cells are defined formally by:

T C i = {x ∈ n : d(x, x i ) < d(x, x j ), ∀j = i, S x→ni ∩ Γ = ∅, S x→nj ∩ Γ = ∅} (3.6)
where Γ is the domain boundary, composed with segments l i ∈ L, L being a set of segment in the plane, and S a→b denotes the segment between the points a and b.

A generalization of the constrained Delaunay triangulation to 3D doesn't exist without adding new nodes, as shown in [START_REF] Schönhardt | Uber die zerlegung von dreieckspolyedern in tetraeder[END_REF]. Nevertheless, some techniques for constructing 3D constrained Delaunay tesselations are available and provided in [START_REF] Sukumar | The natural elements method in solid mechanics[END_REF], [START_REF] Shewchuck | Sweep algorithms for constructing higher-dimensional constrained Delaunay triangulations[END_REF] by addition of Steiner points.

The Constrained Natural Element Approximation. In order to solve partial differential equations defined in non convex domains, or to reproduce some functional discontinuities, we consider the following approximation of both the trial and the test functions:

T h (x) = V i=1 φ C i (x)T i (3.7)
where V is the number of natural neighors visible from point x and φ C i is the constrained natural neighbor shape function related to the i-th node at point x. The computation of the C-n-n (constrained natural neighbor) shape functions is similar to the natural neighbor shape function, when one proceed using the constrained Voronoi diagram introduced previously. It was shown in [START_REF] Yvonnet | Interpolation naturelle sur les domaines non convexes par l'utilisation du diagramme de Voronoi contraint -M Őthode des Ől Őments C-Naturels[END_REF] and [START_REF] Yvonnet | A new extension of the natural element method for non convex and dscontinuous problems, the constrained natural element method (C-NEM)[END_REF] that the use of the constrained Voronoi diagram does not affect the properties of the NEM interpolation, allowing the extension of the linearity of the shape functions on the convex domains boundaries, to any geometry, convex or not.

The ability of the C-NEM for treating problems involving cracks has been illustrated in [START_REF] Yvonnet | A new extension of the natural element method for non convex and dscontinuous problems, the constrained natural element method (C-NEM)[END_REF]. In the present paper, we focus on its application in the context of a moving interface defining two domains with different thermal properties. Thus, defining at time t two CVD (constrained Voronoi diagrams) of Ω 1 (t) and Ω 2 (t), both with respect to the interface Γ I (t), it can be proved that the interpolated temperature field is C 1 everywhere, except at the nodes and on the interface Γ I (t) where it is only C 0 . Thus, this interpolation seems to be appropriate to simulate the Stefan problem considered in this paper.

To illustrate this behavior, we consider the situation depicted in Figure 3.3, where the point x moves from Ω 1 to Ω 2 . If x is in Ω 1 , the interpolated field is constructed from Eq. (3.7) using the neighbor visible nodes from point x (Γ I is assumed opaque). If x is on Γ I , according to the previous discussion, the interpolated field is strictly linear because it only depends on the two neighbor nodes located on Γ I . Finally, when x is in Ω 2 , the interpolated field is defined using the visible neighbor and visible nodes from point x (Γ I being opaque). The continuity of the interpolated field is then guaranteed, but a discontinuity appears in the field derivatives, because of a sudden change in the neighbor nodes across the interface. We can then reproduce the temperature field continuity, as well as the expected flux discontinuity on the interface.

C-NEM Discretization

Assuming that in our example Γ 2 ≡ Γ I , the weak formulation associated with Eq. (2.1) results:

Find T ∈ H 1 (Ω) verifying T = T on Γ 1 such that: Ω c ∂T ∂t δT dΩ = - Ω k∇T • ∇δT dΩ + ΓI (t) |[q]|δT dΓ, ∀ δT ∈ H 1 0 (Ω) (4.1)
Where H 1 (Ω) and H 1 0 (Ω) are the usual Sobolev functional spaces. Substituting the trial and test functions (both approximated in the C-NEM framework) in the above equation and using the arbitrariness of the field δT , the following system of equations is obtained:

C Ṫ + KT = F (4.2)
where T is the vector containg the unknown nodal temperatures. For time discretization, we consider the solution on the time interval [0, t max ], partitioned into steps as [t n , t n+1 ] and the generalized trapezoidal time stepping algorithm characterized by the parameter α:

∂T n+1 ∂t = T n+1 -T n -(1 -α)Δt ∂T n ∂t αΔt (4.3)
where ∂T 0 ∂t is initialized by setting:

Ω c ∂T 0 ∂t δT dΩ = - Ω k∇T 0 • ∇δT dΩ ∀ δT ∈ H 1 0 (Ω) (4.4)
which leads, after applying the discretization scheme described in Eq. (3.5), to:

C ∂T 0 ∂t = -KT 0 (4.5)
The stabilized conforming nodal integration proposed by Chen et al. in [START_REF] Chen | A stabilized conforming nodal integration for Galerkin mesh-free methods[END_REF] is employed for the numerical integration of K (see our former work [START_REF] Yvonnet | A new extension of the natural element method for non convex and dscontinuous problems, the constrained natural element method (C-NEM)[END_REF] for more details). This technique, based on the assumed strain method, has been used to reduce significantly the integration errors, and allows to satisfy the patch test exactly in the context of the natural neighbor interpolation, which is not the case if a standard Gauss quadrature scheme is used [START_REF] Cueto | On the imposition of essential boundary conditions in Natural neighbor Galerkin Methods[END_REF]. In the context of the C-NEM, the representative domains related to each node used in such techniques to define the assumed gradient, are the constrained Voronoi cells depicted in figure 3.3, which are accurately defined everywhere, and especially in the interface neighborhood.

A lumped mass matrix C is computed making use of the constrained Voronoi cells areas as nodal weights.

The iteration procedure is defined as: Knowing T n and |[q]| n at time t n , the non-linear problem associated with Eq. (4.1) results in finding T n+1 and |[q]| n+1 such that Eqs. (2.5) and (4.1) are satisfied. For this purpose we proceed as follows:

1. Compute the interface velocity V n (x) using Eq. (2.3) and update the interface position at time t n+1 using the forward Euler formula:

x n+1 J = x n J + ΔtV n (x n J ) (4.6)
where x J are the nodes defining the interface. 2. Update locally the constrained Voronoi diagram and the shape functions associated with integration points in the interface neighborhood. Then, we compute Cn+1 and K n+1 . 3. Solve Eq. (4.1) using a Newton-Raphson procedure where the tangent matrix is computed numerically. 4. Repeat while t n+1 < t max .

An alternative scheme using the Latin method [START_REF] Ladevèze | Non Linear Computational Structural Mechanics[END_REF] in the extended finite element framework can be found in Merle and Dolbow [START_REF] Dolbow | Natural neighbor Galerkin methods. Solving thermal and phase change with the extended finite element method[END_REF].

Numerical Example

In this section, we illustrate the potentiality of the proposed technique simulating a two-phases Stefan problem. The problem is essentially one-dimensional, but we solve it here in two dimensions to underline the outstanding features of the method.

The Stefan problem models the one-dimensional freezing of a semi-infinite domain (x ≥ 0). The initial temperature T 0 is assumed constant in the whole domain, being higher than the melting temperature T m . At time t = 0 the temperature at the left boundary x = 0 is suddenly prescribed to a value T 1 lower than the melting point, originating a solidification front that progresses from the boundary x = 0 in the x direction. The exact flow front position x f (t) is given by:

x f (t) = 2λ β s t (5.1)
where β s = k s /c s is the thermal diffusivity of the solid phase, and the constant λ satisfies the following relationship:

e -λ 2 erf (λ) = k l √ η(T 0 -T m )e -ηλ 2 k s (T m -T 1 )erf c(λ √ η) + λL √ π c s (T m -T 1 )
(5.2) with η = β s /β l being the ratio of the thermal diffusivities and where k l represents the liquid phase conductivity. The temperature field in the solid phase 0 ≤ x ≤ x f (t) is then:

T (x, t) = T 1 + T m -T 1 erf (λ) erf x 2 √ β s t (5.3)
and in the liquid phase x ≥ x f (t):

T (x, t) = T 0 - T 0 -T m erf c(λ √ η) erf c x 2 √ β l t (5.4)
In the present investigation, we use the water-satured sand thermal properties provided in [START_REF] Lynch | Continuously deforming finite elements for the solution of parabolic problems, with and without phase change[END_REF] and listed in Table I. T 1 and T 0 were set to -10 and 4.0 0 C respectively (λ = 0.3073). We simulate the evolution of the temperature field in Ω = [0, 1] × [0, 0.5] cm. In order to use the inifinite domain solution as reference solution, the temperature at x = 1 is prescibed to its expected value according to Eq. (5.4)

Properties

Solid Liquid

Volumetric heat capacity (cal. 0 C -1 cm -3 ) 0.49 0.62 Thermal conductivity (cal.cm -1 s -1 . 0 C -1 ) 9.6.10 -3 6.9.10 -3 Melting temperature ( 0 C) 0 .0 Volumetric latent heat of fusion (cal.cm -3 ) 19.2

In a first test, we consider in the domain Ω a 20 × 10 uniform grid and Δt = 2s. Figure 5.4 compares the computed interface position and the exact one. An excellent accuracy can be noticed, as depicted in figure 5.5, where the error in the front position is represented. Figure 5.6 shows the temperature profile at different times. We can point out that the discontinuity in the temperature gradient is accurately obtained, being in excellent agreement with the exact solution. The external boundary conditions are given by: T (t) = T 1 on (x = 0), T (t) = T ex given by Eq. (5.4) on (x = 1) and -k∇T • n = 0 on (y=0) and (y=0.5). The initial condition is given by T (t = 0) = T 0 in Ω.

In a second test, we consider the domain Ω containing 200 nodes distributed at random. The purpose of the present test is to investigate the meshless feature of the technique, in which due to its meshless character no geometrical restrictions concerning the relative nodal positions are involved. Thus, neither the background nodal distribution nor the relatinve position of the nodes defining the moving interface with respect to the background nodes, induce a lack of accuracy when high distortions in the Delaunay mesh, used to compute the Voronoi diagram, takes place. This is the main difference between the proposed strategy and the standard finite element method whose accuracy depends significantly on the geometrical quality of the mesh. Moreover, this test approaches the situations encountered when the material is also moving, inducing highly irregular nodal densities and high background mesh distortions. Figure 5.7 depicts the cloud of nodes and the interface position as well as the associated constrained Voronoi cells. Remarkably, despite of the very irregular nodal distribution and density, we can notice that the interface is not distorted as it moves through the domain. From Figs. 5.5 and 5.8 we can conclude that the accuracy is not affected significantly by the resgularity on the nodal distribution. In Figure 5.9 some temperature profiles along the line y = 0.25 are depicted, from which an excellent accuracy can be noticed.

Conclusion

In this paper, the salient features of the C-NEM meshless method are used for treating thermal problems involving moving interfaces. The use of the constrained Voronoi diagram allows to introduce the desired discontinuities in the gradient of the solution without any enrichment. In the C-NEM framework, essential boundary conditions can be imposed directly, due to the strict linearity of the shape functions over the boundaries (convex or not) and the respect of the Kronecker delta property. No user-defined parameter is involved in the shape functions support size, the support being defined like the union of the Delaunay spheres passing through a node and its neighbor visible nodes. A stabilized conforming nodal integration is computed over the constrained Voronoi cells to enhance accuracy. The most outstanding quality of this technique is the ability for introducing discontinuities located on a line (surface in 3D) defined by another set of nodes that move through a fixed (or also moving) cloud of nodes which constitutes a background scattered of nodes. The meshless character of this technique allows to proceed without remeshing (in a finite element sense), even for highly irregular nodal distributions. The method is relatively simple and delivers similar level of accuracy than r-adaptative or the partition of unity finite element schemes. The technique seems promising for the simulation of arbitrary and moving discontinuities over an unstructred fixed or also moving set of nodes.

Figure 2 . 1 .

 21 Figure 2.1. Two phases problem.

ΓFigure 3 . 2 .

 32 Figure 3.2. Construction of the Sibson shape functions.

Figure 3 . 3 .

 33 Figure 3.3. Reproducing discontinuous derivatives using the constrained Voronoi diagram.

Figure 5 . 4 .

 54 Figure 5.4. Computed C-NEM front position versus the exact solution using a 20 × 10 regular grid.

Figure 5 . 5 .Figure 5 . 6 .

 5556 Figure 5.5. Error in the front position.

Figure 5 . 7 .

 57 Figure 5.7. Computed position irregular cloud of nodes: (a) Cloud of nodes and interface position; (b) Constrained Voronoi cells.

Figure 5 Figure 5 . 9 .

 559 Figure 5.8. C-NEM interface position versus the exact solution using an irregular grid.

Table 5 . 1 .

 51 Thermal properties of the water saturated sand.


Keywords: Experimental fabric forming. Finite element simulations, Unbalanced fabric, Fabric mechanical behaviour

The deep drawing of a very unbalanced fabric on a hemispherical shape is performed both experimentally and by a numerical simulation. The shape obtained depends mainly on the mechanical behaviour of the fabric and is highly non-symmetric. The simulation of the forming process is performed using finite elements based on the tensile energy of each yarn. This approach leads to very simple expression of the interior load vectors used in the explicit dynamic approach. The result of the simulation correctly describes the shape after forming and its non-symmetry. Such a result could not have been obtained by a fishnet method classically used for fabric draping.

The use of fibre fabrics as reinforcements for composite materials is increasing because it gives the possibility to obtain complex shapes in a single operation, for a lighter final product. For instance the R.T.M. process consists of a drawing operation of the fabric before a thermoset resin is injected [CAR 96][BIC 97] [LON 98][POT 99]. The forming modes of fabrics are specific to this material and related to its woven constitution. A simulation code permits to answer two questions at the conception stage: Firstly, is the requested shape reachable by forming the fibre fabric? Secondly, where are the reinforcement directions after forming? Several codes based on geometrical approaches have been developed [MAC 56] [ VAN 91][LON 94][BIC 97][WAN 99]. These methods, called fishnet algorithms, are very fast and efficient especially for the simulation of hand operated draping. Nevertheless in the case of forming with punch and die the static boundary conditions may be very important. So is the influence of the loads on the blank holder. In addition the specifications of the fabric, i.e. the material for the warp and the weft, their respective size or the way it is woven, should be included in the analysis. An objective of this paper is to show that the mechanical properties of the fabric can be very important for the result of the shaping process. The mechanical behaviour of dry fabrics such as used in the forming stage of the R.T.M. process depends on the yarn behaviour and on the weave pattern, leading to interactions between warp and weft directions and non-linear behaviour. These geometrical nonlinearities due to changes in yarn crimp under tension have been studied from biaxial tests on cross-shaped specimens [BUE 01][BOI 01b]. The fabric studied in the present paper is very unbalanced. The rigidity in warp direction is much larger than in the weft direction.

An experimental forming test is presented that shows the very large effect of the non-symmetry of the fabric behaviour on the forming result. This forming process is simulated, within an explicit dynamic approach, using simplified finite elements in which the interior load vectors are calculated from the tensile energy of each yarn. This approach gives a computed fabric shape after forming in good agreement with the experimental result. The test performed experimentally and simulated clearly shows that it is necessary in this case to use a method that accounts for the fabric mechanical behaviour (and for the boundary condition) to obtain a good simulation result.

Deep drawing tests

Punch and die geometry

The tests were performed in the composites laboratory of the School of MMMEM at the University of Nottingham [RUD 99]. Cold press mouldings were held with a hemispherical punch-die couple, mounted on a classical 25kN-Instron machine (Figure 1).

A 6 kg ring was used as a blank-holder, to avoid premature wrinkles in the reinforcement. The preforms were impregnated with a translucent polyester resin after forming to freeze their movements. A grid of 10mmx10mm elements was drawn on the fabric prior to forming to measure the fibre angle variations and the elongations in each direction. The tool geometry (including the blank-holder) is shown in Figure 2.

Very unbalanced fabric

The fabric studied here is a 2x2 nylon twill, used in automotive parts as a reinforcement in elastomer composites. Each tow contains approximately 40 fibres. The peculiarity of this reinforcement lays in the weft tows. They are textured, so that they have a double crimped path: the normal weaving crimp and their own waviness due to their texturizing (Figure 3). This strong difference between warp and weft tows results in a very unbalanced mechanical behaviour, as shown in Figure 4. Uni-axial traction tests exhibit a nonlinear behaviour. But, for a given range of elongation, it is possible to define an equivalent linear modulus in each direction (given in Table 1). Table 1. Equivalent modulus.

Forming experiment

The size of the initially square elements of the grid changes when forming. The elongation ratio between warp and weft direction can be measured easily along a line of elements (Figure 5a). At the top of the hemisphere, the weft elongation is maximum, and the warp elongation remains very low. Thus, the ratio reaches 1/1.8, which is very specific to this unbalanced fabric. In the meantime, warp and weft deformed edges of the blank (Figure 5b) exhibit strong differences. The shear angle has been measured on one quadrant (Figure 5b) by an optical method in use at the University of Nottingham. The deformed grid (Figure 6) exhibits a highest shear angle of 39°, at the base of the hemisphere. It is not symmetrical, this angle is maximum at the approximate meridian direction = 30°, from the warp direction (so-called 0° direction). On a balanced fabric, due to the symmetry, the maximum shear angle direction is located at = 45° [AUB 89][RUD 99][YE 97].

Simulation of the forming process

Simplified dynamic equation

Fabrics are made of fibres assembled in yarns which are woven. This internal structure permits some motions between the fibres. Consequently most stiffnesses are very weak in comparison to tensile stiffness. Of course, if these possible motions (or nil rigidities) are an advantage during the forming stage they are not acceptable in the service life of the part, and it is one of the main roles of the resin that is injected after the preforming of the fabric to avoid the relative motion of the fibres. However in the forming stage the fabric is "dry" and the stress state can be assumed to be in the form :

1 + 1 1 2 h h h h 2 2 [1]
is the Cauchy stress tensor and and are unit vectors in the current warp and weft directions. Defining the tension tensor by f and t are exterior prescribed volume and surface loads.

1 h 2 h 1 + 1 1 2 T T T h h h h [2] (a) (b) (c)

Explicit finite element analysis

The simplified dynamic equation [4] is used in an explicit dynamic finite element approach. (An implicit quasi static approach can be found in [BOI 97]). If the domain is divided in finite elements, the dynamic equation can be written as : The finite element briefly described below is based on the computation of the interior loads given by equation [8] with a bilinear interpolation scheme.

Interior loads of a bilinear element

The four node element presented Figure 7 is composed of woven yarns. The numbers of yarns in the warp and weft directions are N c and N t . From the natural coordinates and , the following vectors are defined :

x g g h g g . g = h . h = 11 [9]
In the present element the interactions between warp and weft yarns are not taken into account. This is justified for the fabric presented in section 2.2. because the warp yarns are very much stiffer than the weft yarns. For other fabrics the effect of the weaving and the interaction between warp and weft yarns can be important. In this case the tensile behaviour must be obtained by biaxial tensile tests [

T ( ) T ds [10]

The components of the gradient tensor are considered in h 1 , h 2 and g 1 , g 2 bases:

s h h g g [11]
The component of the gradient 11 is related to the interpolation functions :

k 11 1 k 1 11s s N .g . B g [12]
The virtual interior work can be written as N C and N T are the number of yarns in warp and weft directions respectively. The quantities B 11s , T 11 , g 1 (and similarly B 22s , T 22 , g 2 ) are evaluated at positionp for p=1 and p for p=2, (and similarly for p ) with

2 c p 2 c N 1 1 3 N 2 t p 2 t N 1 3 N 1 [17]
The interior loads components are given explicitly as the sum of four scalars. Consequently the element is numerically efficient. Within a dynamic explicit approach such as described in section 3.2. it is classical to use reduced integration. In that case the terms of equation [16] are only computed at the center of the element. Such an approach must avoid the zero energy mode by using a specific method [BEL 83b]. This element is included in the PAM FORM code [DEL 98].

. Simulation of the forming of a very unbalanced fabric

The experimental forming operation described in section 2 is simulated using the four node fabric element described in section 3. The initial mesh used for the fabric, the punch, the die and the blank-holder are shown in Figure 8.

Figure 8. Initial mesh of the fabric, the punch, the die and the blank-holder

The deformed shape obtained at the end of the punch displacement is shown in Figure 9 and compared to the experimental shape. The global shape is correctly described by the computation. The main feature of this deformed shape is that it is very different in the warp and weft directions. The weft direction is the weak direction and it is horizontal on Figure 9. The stretching in this direction is very large in the central part due to the blank holder. The displacements of points on the edge and on the horizontal axis are very small. On the contrary the stretching in the warp (vertical) direction is small and consequently the points on the edge and on the vertical axis exhibit large displacements towards the centre of the hemisphere. The computed deformed shape does not contain any wrinkles. This is due to the lack of instability in the solutions obtained from the present formulation that only considers tensile deformation energy. Taking into account the in-plane shear stiffness, especially after the limit angle would be a source of wrinkles. Some works are now in progress at LMSP in that goal. To correctly describe the wrinkle shape, it is also necessary to add some bending stiffness of the fabric in the formulation. Nevertheless many of the wrinkles in the flat part are due to removal of the blankholder and there are no wrinkles in the hemispherical part. In this region, the global shape and the in plane strain of the fabric after forming are correctly simulated by the above presented element. In the hemispherical part, the initial squares drawn on the fabric are strongly deformed and become rectangular. At the top of the hemisphere the ratio of the side length in weft and warp directions is 1.8 (Figure 10). The simulation is in good agreement with this value and gives 1.87 for this ratio. If this forming process is simulated using a classical geometrical approach (fishnet algorithm) this ratio would be equal to 1 and the computed shape would be symmetrical which is very far from the experimental results.

Figure 11 shows the in plane shear deformation i.e. the angles between warp and weft yarns after forming. The position of the high shear zone is in correct agreement with the experiment even if the wrinkles in the experimental part modify the value of the shear angle. The angles between warp and weft yarns are not symmetrical and depend mainly on the imbalance of the fabric.

Figure 13. Tensions in the warp direction (N/yarn)

Finally Figures 12 and13 give the tensions in the warp and weft directions. The computation of these tensions is an important advantage of finite element approach. For instance the maximal possible load on the blank holder can be determined by these yarn tensions. The computed tensions and tensile strains are large in the hemispherical part especially for the weft yarns which are not very far from the ultimate value. The tensile strain is 0.87 while the tensile strain at failure is 1.1.

Conclusions

An hemispherical forming experiment performed on a very unbalanced fabric has shown that the result is highly non-symmetric. This result cannot be obtained by classical geometric draping algorithms. A simplified finite element based on the tension energy of each yarn has been used with an explicit dynamic approach. The computed shape after forming is in correct agreement with the experimental one. The wrinkles are not modelled at present. With this in mind elements with in plane shear and bending stiffness are now being developed.

This example clearly shows that the mechanical behaviour of the fabric can be important in a forming process. In this case a mechanical approach is necessary to account for this aspect. The static boundary conditions i.e. the loads on the tools can also be very important for the final result of forming. Their modelling also needs a finite element approach.

Acknowledgement

This work was performed with the support from PSA, EADS and ESI companies. This support is gratefully acknowledged by the authors.

Forming of a very

  unbalanced fabric. Experiment and simulation. J.L. Daniel * -D. Soulat*-F. Dumont * -B. Zouari * -P. Boisse * -

Figure 1 .

 1 Figure 1. Punch-die couple

Figure 2 .Figure 3 .Figure 4 .

 234 Figure 2. Tool geometry

Figure 5 .Figure 6 .

 56 Figure 5. Shape of the unbalanced fabric after forming

Figure 7 .

 7 Figure 7. Bilinear fabric element

Figure 9 .Figure 10 .

 910 Figure 9. Deformed shape. Comparison between experiment (left) and computation (right).

Figure 11 .Figure 12 .

 1112 Figure 11. Angles between warp and weft yarns after forming

  M is the mass matrix, F ext and F int are exterior and interior nodal loads. Subscript n refers to nodal quantities and superscript e refers to an element. In the central

										Mu	n			e x t F	i n t F	[6]
	difference scheme the solution	u	i n	1	on a time step t i , from t i to t i+1 is calculated
	i from u by : n
			u	i 1 n	u	i n	u	i 1/ 2 n		1 t 2	i 1	t	i	M	-1 D	i e x t F	i i n t F	t	i	[7]
	M D is a diagonal matrix calculated from the mass matrix [ZIE 89]. The
	scheme [7] is stable only if the time step is small enough [BEL 83a]. From the
	simplified dynamic equation [4], a particular form is obtained for equation [6] where
	the virtual interior work (and consequently the internal loads) is in the form :
		int T ( )		n cell	e	e n	T	e int F					e n cell	p	11	p 1 1p T L 1	p	22	p 2 2p T L	2	[8]
					elt	p=1							elt			p=1
	with	int F	F int e A elt	where	elt A	is the assembling operator over the elements.

  BUE 01] [BOI 01b]. Finite elements including the interaction due to weaving have been developed in [BOI 97][BOI 01a]. Let us consider firstly the yarns in direction g 1 .For one element, the interior virtual work can be written as:

	e
	int	yarns	l	11


Apoloniusz Tyszka

A conjecture which

Keywords:

implies that there are infinitely many primes of the form n! + 1

It is conjectured that n! + 1 is prime for infinitely many positive integers n, see [START_REF] Caldwell | On the primality of n! ± 1 and 2 × 3 × 5 × • • • × p ± 1[END_REF]. In this note, we propose a conjecture which implies this. Lemma 1. For every integers x and y greater than 1, x! • y = y! if and only if x + 1 = y.

For an integer n 6, let V n denote the following system of equations:

                   ∀i ∈ {1, . . . , n -1} \ {1, 3, 5} x i ! = x i+1 x 1 ! = x 6 x 1 • x 2 = x 3 x 3 • x 4 = x 5 x 1 • x 5 = x 6
The diagram in Figure 1 illustrates the construction of the system V n .

x 2 ! x 3 x 4 ! x 5 x 1 ! x 6 +1 +1 ! x 7 . . . x n-1 ! x n x 1 • x 2 = x 3 x 3 • x 4 = x 5 x 1 • x 5 = x 6
Fig. 1 Construction of the system V n Let f (6) = 720, and let f (n + 1) = f (n)! for every integer n 6.

Lemma 2. For every integer n 6, the system V n has exactly one solution in integers greater than 3, namely 6, 4, 24, 5, 120, f (6), . . . , f (n) .

Proof. By Lemma 1,

x 2 + 2 = x 1 . Hence, (x 2 + 2) • x 2 = x 1 • x 2 = x 3 = x 2 !. Therefore, x 2 = 4.
The rest of the proof follows from the diagram in Figure 1.

Let T n = x i ! = x j : 1 i < j n ∪ x i • x j = x j+1 : 1 i < j n -1
For an integer n 6, let Λ n denote the following statement: if a system S ⊆ T n has at most finitely many solutions in integers x 1 , . . . , x n greater than 3, then each such solution (x 1 , . . . ,

x n ) satisfies x 1 , . . . , x n f (n).
Conjecture. The statements Λ 6 , . . . , Λ 9 are true.

We present a heuristic reasoning that leads to the Conjecture. For every integer n ∈ {6, 7, 8, 9}, we consider all subsystems of the system T n . We conjecture that the largest known solution is indeed the largest possible.

Theorem 1. For every statement Λ n , the bound f (n) cannot be decreased.

Proof. It follows from Lemma 2 because V n ⊆ T n .

Theorem 2. Every statement Λ n holds true with an integer bound that depends on n.

Proof. Indeed, for every integer n 6, the system T n has a finite number of subsystems.

Lemma 3. If a prime number x is grater than 3, then

(x -1)! + 1 x > 3.
Lemma 4. (Wilson's theorem, [4, p. 89]). For every integer x 2, x is prime if and only if x divides (x -1)! + 1.

Let A denote the following system of equations:

                                   x 1 ! = x 2 x 2 ! = x 3 x 3 ! = x 4 x 5 ! = x 6 x 8 ! = x 9 x 3 • x 5 = x 6 x 4 • x 8 = x 9 x 5 • x 7 = x 8
Lemma 1 and the diagram in Figure 2 explain the construction of the system A.

x 1 ! x 2 +1 x 5 ! x 6 ! x 3 ! x 4 +1 x 8 ! x 9 x 5 • x 7 = x 8 x 3 • x 5 = x 6 x 4 • x 8 = x 9
Fig. 2 Construction of the system A Lemma 5. For every integer x 1 > 3, the system A is solvable in integers x 2 , . . . , x 9 greater than 3 if and only if x 1 ! + 1 is prime. In this case, the integers x 2 , . . . , x 9 are uniquely determined by the following equalities:

x 2 = x 1 ! x 3 = (x 1 !)! x 4 = ((x 1 !)!)! x 5 = x 1 ! + 1 x 6 = (x 1 ! + 1)! x 7 = (x 1 !)! + 1 x 1 ! + 1 x 8 = (x 1 !)! + 1 x 9 = ((x 1 !)! + 1)! Proof.
By Lemmas 1 and 3, for every integer x 1 > 3, the system A is solvable in integers x 2 , . . . , x 9 greater than 3 if and only if x 1 ! + 1 divides (x 1 !)! + 1. Hence, the claim of Lemma 5 follows from Lemma 4. Theorem 3. The statement Λ 9 implies that there are infinitely many primes of the form n! + 1.

Proof. Harvey Dubner proved that 872! + 1 is prime, see [START_REF] Caldwell | The List of Largest Known Primes Home Page[END_REF], [START_REF] Caldwell | On the primality of n! ± 1 and 2 × 3 × 5 × • • • × p ± 1[END_REF], [5, p. 7], and [START_REF] Sloane | The On-Line Encyclopedia of Integer Sequences, A002981[END_REF]. Let x 1 = 872. By Lemma 5, there exists a unique tuple (x 2 , . . . , x 9 ) of integers greater than 3 such that the tuple (x 1 , x 2 , . . . , x 9 ) solves the system A. Hence,

x 9 = ((x 1 !)! + 1)! > ((720!)! + 1)! > ((720!)!)! = (( f (6)!)!)! = f (9)
The statement Λ 9 and the inequality x 9 > f (9) imply that the system A has infinitely many solutions in integers x 1 , . . . , x 9 greater than 3. This conclusion and Lemma 5 imply that x 1 ! + 1 is prime for infinitely many integers x 1 > 3.

Let C denote the following system of equations:

                   x 1 ! = x 4 x 2 ! = x 5 x 3 ! = x 6 x 2 • x 3 = x 4 x 3 • x 5 = x 6
Lemma 1 and the diagram in Figure 3 explain the construction of the system C.

x 2 ! x 5 x 3 ! x 6 +1 x 3 • x 5 = x 6 x 1 ! x 4 x 2 • x 3 = x 4
Fig. 3 Construction of the system C Lemma 6. For every integers x 1 , x 2 greater than 3, the system C is solvable in integers x 3 , x 4 , x 5 , x 6 greater than 3 if and only if x 1 ! = x 2 (x 2 + 1). In this case, the integers x 3 , x 4 , x 5 , x 6 are uniquely determined by the following equalities:

x 3 = x 2 + 1 x 4 = x 1 ! x 5 = x 2 ! x 6 = (x 2 + 1)! Theorem 4. If the equation x 1 ! = x 2 (x 2 + 1)
has at most finitely many solutions in positive integers, then the statement Λ 6 guarantees that each such solution (x 1 , x 2 ) belongs to the set {(2, 1), (3, 2)}.

Proof. Suppose that the antecedent holds. Then, the equation x 1 ! = x 2 (x 2 + 1) has at most finitely many solutions (x 1 , x 2 ) ∈ (N \ {0, 1, 2, 3}) 2 . By Lemma 6, the system C is solvable in integers x 3 , x 4 , x 5 , x 6 greater than 3. Since C ⊆ T 6 , the statement Λ 6 implies that x 1 ! = x 4 f (6) = 720 = 6!. Hence, x 1 ∈ {1, 2, 3, 4, 5, 6}. For every integer x 1 ∈ {1, 2, 3, 4, 5, 6}, x 1 ! is a product of two consecutive integers if and only if x 1 ∈ {2, 3}.

The question of solving the equation y! = x(x + 1) was posed by P. Erdös, see [START_REF] Berend | On polynomial-factorial Diophantine equations[END_REF]. F. Luca proved that the abc conjecture implies that the equation y! = x(x + 1) has only finitely many solutions in positive integers, see [START_REF] Luca | The Diophantine equation P(x) = n! and a result of M. Overholt[END_REF].