N Sukumar

D L Chopp

N Mo Es

T Belytschko

Modeling holes and inclusions by level sets in the extended ®nite-element method

Keywords: Extended ®nite-element method, Level set method, Voids, Material interface, Elastic inclusion

A methodology to model arbitrary holes and material interfaces (inclusions) without meshing the internal boundaries is proposed. The numerical method couples the level set method (S. Osher, J.A. Sethian, J. Comput. Phys. 79 (1) (1988) 12) to the extended ®niteelement method (X-FEM) (N. Mo es, J. Dolbow, T. Belytschko, Int. J. Numer. Methods Engrg. 46 (1) (1999) 131). In the X-FEM, the ®nite-element approximation is enriched by additional functions through the notion of partition of unity. The level set method is used for representing the location of holes and material interfaces, and in addition, the level set function is used to develop the local enrichment for material interfaces. Numerical examples in two-dimensional linear elastostatics are presented to demonstrate the accuracy and potential of the new technique.

Introduction

Defects such as pores, cracks, and inclusions are important to the structural integrity and durability of components. For example, material interfaces in composites are modeled to predict the mechanical behavior and to establish macroscopic material properties. The accurate modeling of voids and inclusions is hence of interest to both the theorist and the practitioner. In ®nite-element methods, the mesh is required to conform to internal boundaries in the model. Although mesh generation is very well established in two dimensions, meshing arbitrary number and distribution of defects and inclusions is still a time-consuming and burdensome task.

In this paper, we describe the modeling of internal boundaries from a dierent viewpoint, by using the level set method within the framework of the extended ®nite-element method (X-FEM) [START_REF] Mo Es | A ®nite element method for crack growth without remeshing[END_REF]. In the X-FEM, the ®nite-element approximation is enriched by additional functions through the notion of partition of unity [START_REF] Melenk | The partition of unity ®nite element method: Basic theory and applications[END_REF]. The recognition and use of partition of unity enrichment strategy to solve boundary-value problems with internal boundaries is due to Oden and co-workers [START_REF] Duarte | Generalized ®nite element methods for three dimensional structural mechanics problems[END_REF][START_REF] Duarte | Generalized ®nite element methods for three dimensional structural mechanics problems[END_REF][START_REF] Oden | A new cloud-based hp ®nite element method[END_REF],10] ± the numerical technique was coined as the generalized ®nite-element method (GFEM). Strouboulis et al. [START_REF] Strouboulis | The design and analysis of the generalized ®nite element method[END_REF] used local enrichment functions in the GFEM for modeling re-entrant corners in two dimensions, and in [START_REF] Strouboulis | The generalized ®nite element method: An example of its implementation and illustration of its performance[END_REF], enrichment functions for holes are proposed.

In the X-FEM, the additional functions are used to model cracks, voids or inhomogeneities, and also to improve accuracy in problems where some aspects of the functional behavior of the solution ®eld are known a priori. As opposed to the GFEM in which the analytical solution of a given boundary-value problem is embedded into the ®nite-element space, in the X-FEM the proximity to the ®nite-element method is consistently sought by carrying out only the requisite enrichment to model the internal boundary (crack or inclusion) of interest. The GFEM is amenable to accurate numerical solutions with coarse meshes and provides the capability of p-type adaptivity. In the X-FEM, less dependence on known closed-form solutions is emphasized which provides for greater ¯exibility to solve a wider range of problems. For a few applications in the above spirit, see [START_REF] Dolbow | Modeling fracture in Mindlin±Reissner plates with the extended ®nite element method[END_REF][START_REF] Mo Es | A ®nite element method for crack growth without remeshing[END_REF] for crack modeling in two dimensions, and [START_REF] Sukumar | Extended ®nite element method and fast marching method for three dimensional fatigue crack propagation[END_REF][START_REF] Sukumar | Extended ®nite element method for three-dimensional crack modeling[END_REF] for threedimensional elastostatics and fatigue crack growth modeling.

The level set method is a numerical technique for tracking moving interfaces [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton±Jacobi formulations[END_REF]. The evolving interface is represented as the zero level set of a function of one higher dimension (i.e., ux; t 0). The motion of the interface is embedded in a hyperbolic equation in terms of ux; t. In this paper we pursue twodimensional elastostatic applications, and hence the interface is static. The interface for geometric entities such as holes and inclusions is represented by the level set function. In addition, using u, the local enrichment function for material interfaces is constructed.

The outline of this paper is as follows. In the next section, we give a brief introduction to the X-FEM, and in Section 3, the level set method is introduced. In Section 4, the modeling of voids and inclusions is addressed. The strong and weak forms of the Galerkin method for two-dimensional elastostatics, along with the discrete equations for the X-FEM are given in Section 5. In Section 6, two dierent approaches for the numerical integration of the weak form are described. Numerical results are presented in Section 7, and we close with some concluding remarks in Section 8.

The X-FEM

In the ®nite-element method, the presence of ¯aws or inhomogeneities such as cracks, voids, and inclusions must be taken into account in the mesh generation process ± the edges of the element must conform to these geometric entities.

The X-FEM aims to alleviate much of the burden associated with mesh generation for problems with voids and interfaces by not requiring the ®nite-element mesh to conform to internal boundaries. The essence of the X-FEM lies in sub-dividing a model into two distinct parts: mesh generation for the domain (excluding internal boundaries); and enriching the ®nite-element approximation by additional functions that model the internal boundaries.

Consider a domain X & R 2 which is partitioned into ®nite elements, and let N fn 1 ; n 2 ; . . . ; n m g be a set of m nodes in the mesh. In addition, let X g & X be a region for which some enrichment is required. For a material interface, the domain X g is the curve representing the material interface and in the case of a void, X g is the region occupied by the void. The X-FEM displacement approximation assumes the general form:

u h x I n I PN / I xu I J n J PN g / J xwxa J u I ; a J P R 2 À Á ; 1
where the nodal set N g is de®ned as

N g n J : n J P N; x J X g È Y É : 2
In the above equation, x J supp n J is the support of the nodal shape function / J x, which consists of the union of all elements with n J as one of its vertices. The choice of the enrichment function wx for voids and material interfaces will be described in Section 4.

The level set method

The level set method [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton±Jacobi formulations[END_REF] is a numerical technique for tracking moving interfaces. It is based upon the idea of representing the interface as a level set curve of a higher-dimensional function ux; t. In keeping with the applications pursued in this paper, we consider one-dimensional curves in R 2 .

A moving interface Ct & R 2 can be formulated as the level set curve of a function u : R 2 Â R 3 R, where

Ct x P R 2 : ux; t È 0 É : 3
One important example of such a function u would be the signed distance function

ux; t AE min x C PCt kx À x C k; 4
where the sign is positive (negative) if x is outside (inside) the contour de®ned by Ct (we assume that the interface Ct is such that one can de®ne an interior and exterior to it). The evolution of the interface is then embedded in the evolution equation for u, which is given by Osher and Sethian [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton±Jacobi formulations[END_REF] u t F k$uk 0; 5a ux; 0 given; 5b where F x; t is the speed of the interface at x P Ct in the direction of the outward normal to the interface.

The key advantages of this method are that it is computed on a ®xed Eulerian mesh, handles topological changes in the interface naturally, and can easily be formulated in higher dimensions. Application of the level set method to various problems in engineering and the applied sciences can be found in [START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry[END_REF].

For our purposes in this paper, the interface is static, so we only use the level set theory for representing the interface. For circular voids, we use ux; 0 min

x i c PX i c i1;2;...;nc kx È À x i c k À r i c É ; 6
where X i c is the domain of the ith void, n c is the number of circular voids, and x i c and r i c are the center and radius of the ith void, respectively.

Consider an elliptical void with semi-major and semi-minor lengths of a and b, respectively, with the axes oriented along the local x1 and x2 axes, respectively. In addition, let ê Re be the mapping between the global and the local unit base vectors. Now, the coordinates of a point in the local orthogonal coordinate system with the origin at the center of the ellipse is given by: x Rx À x c , where x c is the center of the ellipse. Letting n n 1 ; n 2 , and n 1 x1 =a and n 2 x2 =b, the equation of the ellipse in the local coordinate system is: f n knk À 1 0. Hence a level set function for elliptical voids is given by Ux; 0 min i1;2;...;ne f n i ; 7a

f n i kn i k À 1; 7b n i xi 1 a i ; xi 2 b i 2 3 ; 7c xi R i x À x i c ; 7d
where n e is the number of elliptical voids. Note that in Eq. (7a), U is a level set function, but not a signed distance function. In order to obtain the signed distance function u to an elliptical void (interface), one applies the fast marching method to Eq. ( 9) with Gx 1. Now, consider a polygonal void, whose interface C p (C p p i1 l i ) consists of p segments: l 1 x 1 ; x 2 , l 2 x 2 ; x 3 ; . . . ; l p x p ; x 1 . The level set function for a polygonal interface is given by ux; 0 kx À x min ksgn x À x min Á n min ; 8a

sgn n 1 if n P 0; À1 if n < 0; & 8b kx À x min k min xiPli i1;2;...;p kx À x i k; 8c
where x min is the orthogonal projection of x on the interface C p and n min is the outward normal to the interface at x min . If no unique normal is de®ned at x min , the sign is positive if x À x min belongs to the cone of normals at x min and negative otherwise. In Fig. 1, the level set function for a hexagonal interface is shown. For more general interfaces given only as a parameterized curve, the signed distance function is most eciently calculated using the fast marching method [START_REF] Sethian | A marching level set method for monotonically advancing fronts[END_REF], which optimally solves an equation of the form

k$uk 1 Gx ; 9
where G : R 2 3 R is given. If we take Gx 1 and apply the fast marching method to Eq. ( 9), we can obtain the signed distance function to the interface. This method can be applied to both rectilinear meshes [START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry[END_REF] and also to general triangulated meshes [START_REF] Barth | Numerical schemes for the Hamilton±Jacobi and level set equations on triangulated domains[END_REF].

Modeling voids and inclusions

Level set methods are used to track the evolution of material interfaces. The geometric description of an interface (for instance, boundary of a hole or inclusion) is represented by the zero level set curve u ux; t 0. In essence, the physical description of the interface is converted to a discrete functional representation through the function ux; t. The geometrical degrees of freedom at a ®xed set of points x I (nodes) are used to determine u, and hence the location of the interface. In this study, each ®nite-element node is associated with a geometrical degree of freedom for the level set function. Finite-element interpolation is used to compute u at any point x in the domain ux

I / I xu I ; 10
where the summation is over all nodes in the connectivity of the element that contains x, / I x are standard ®nite-element shape functions, and u I are the nodal values of the level set function. The internal boundaries are then approximately piecewise linear for bilinear quadrilateral elements and exactly linear for threenoded triangular elements.

In the X-FEM, there is no need to mesh internal boundaries. The computational geometry issues that are associated with the ®nite-element mesh and the internal geometric entities (such as cracks or inclusions) are an important consideration in the X-FEM computations. In light of this, the level set method is an appealing choice that would greatly simplify and speed-up the geometric computations in the X-FEM. In addition, the level set function can also be used to construct enrichment functions within the X-FEM framework. This provides the motivation and seemingly natural choice to use the level set method in conjunction with the X-FEM.

Voids

The modeling of voids in the X-FEM is carried out using an enrichment function V x for nodes that intersect the boundary of the void. If the node lies outside the void V x 1, and V x 0 if the node is in the interior of the void. In [START_REF] Daux | Arbitrary cracks and holes with the extended ®nite element method[END_REF], the discretization of the boundary of the hole was independent of the ®nite-element mesh. The degrees of freedom associated with nodes whose support is entirely inside the void are removed from the system of equations. In the two-dimensional implementation, instead of using V x, elements that intersect the boundary of the hole are partitioned such that the weak form is not integrated in the void interior [START_REF] Daux | Arbitrary cracks and holes with the extended ®nite element method[END_REF].

In this work, the point of departure in the modeling of voids lies in using the level set method to represent the voids, and all geometric computations that pertain to the void are carried out using the level set function u. The level set method provides a convenient and appealing means to model voids with greater ¯exibility and better computational eciency. In addition, the ®nite-element mesh and the internal geometric boundary are linked, whereby a self-consistent representation of the internal boundary is obtained for a given ®nite-element mesh. The computation of the level set function for voids is described in Section 3. For circular voids, the geometrical degrees of freedom indicated in Eq. ( 10) are evaluated using Eq. ( 6)

u I min x i c PX i c i1;2;...;nc kx I È À x i c k À r i c É : 11

Inclusions

Inclusions are inhomogeneities in a matrix with diering material properties. The modeling of inclusions requires the satisfaction of the Hadamard condition, namely

F À F À a n ; 12
where F is the deformation gradient, n is the outward normal to the material interface, and a is an arbitrary vector in the plane. In [START_REF] Krongauz | EFG approximation with discontinuous derivatives[END_REF], a technique for incorporating discontinuities in derivatives for meshless methods was proposed. The use of the level set function u in the X-FEM provides an elegant and simpli®ed means to achieve the same end. Using wu as the local enrichment, the X-FEM displacement approximation for two-dimensional domains with inclusions assumes the form:

u h x I n I PN / I xu I J n J PN I / J xw ux a J u I ; a J P R 2 À Á ; 13
where the nodal set N I is de®ned as

N I n J : n J P N; x J C I f Yg; 14
where C I is the boundary of the inclusion (material interface), and the nodal set N is de®ned in Section 2. The enrichment function wux juxj is used in the elements that intersect the material interface, and it is smoothed away from the interface. By construction, since u ;n T 0, the displacement derivatives u i;n are discontinuous across C I . The level set functions for inclusions follows that which is described in Section 3 for voids. We associate a level set function u with every distinct material-pair interface (A±B). Hence, all inclusions with interfacepair A I ±B I are represented by the same level set function, and the number of distinct level set functions is equal to the number of distinct material-pair interfaces. The geometrical degrees of freedom u I for circular inclusions are evaluated analogous to that given in Eq. (11).

Nodal enrichment

The level set function u is used to develop the local enrichment function for material interfaces. In order to determine the nodes that need to be enriched, the level set function u is used. We ®rst loop over the elements of the mesh. For a given element, if there exists two nodes n I and n J in the connectivity of the element such that u I u J < 0, we add the element t to the set T of ®nite elements that intersect the interface. The enriched nodal set N e consists of all nodes that are in the connectivity of the elements in the set T. The above algorithm ensures that if a node belongs to two edges that both lie on the interface, then the node is not enriched. Hence the ®nite-element space is obtained (no enrichment) if the interface lies precisely on element edges.

Let wu be the enrichment function for a material interface. A ®rst and apparent choice for the enrichment function is: wu juj. A second alternative is to take wu juj only in the elements that intersect the interface. Since the enrichment functions exist over the domain of support of the nodes, these enrichment functions are required to be computed one element layer on either side of the elements that intersect the material interface. If N e is the set of nodes that are enriched, and u I are the level set function values at these nodes, then w I u I at these nodes. Let N w fn p1 ; n p2 ; . . . ; n q g be the set of nodes that belong to the adjacent element (outward and inward from the interface), such that n I P N w are sorted in increasing value of the level set function u. Initially, the value of w I is set to u I for the nodes in N w , then we compute a new value for w at these nodes based on the following algorithm: • for J p 1 to J q create the set P J , where P J n I P N e N w : jw I j È < ju J j and n I ; n J share an edge in the mesh É :

minimize the following expression to obtain w J :

min n I PP J w J À u I L I 2 ; 15
where L I is the distance from node n I to node n J .

• end for

The solution to the minimization problem in Eq. ( 15) is explicitly given by

w J n I PP J a I u I ; a I 1=L 2 I n K PP J 1=L 2 K ; 16
and hence the modi®ed enrichment function is given by

wx w 1 x I / I xw I : 17

Governing equations

In this section, we present the governing equations of linear elastostatics, together with the weak form and the discrete system for the X-FEM. The case of internal boundaries is considered in the presentation.

Strong form

Consider a body which is described by an open bounded domain

X & R 2 , with boundary C. The boundary C is composed of the sets C u , C t , C i h and C i I , such that C C u C t m i1 C i h n i1 C i I .
All the internal holes C i h are assumed to be traction-free, and traction continuity holds along the material interfaces C i I . The ®eld equations of elastostatics are:

$ Á r b 0 in X; 18a r C : e; 18b e $ s u; 18c
where $ s is the symmetric gradient operator and C is the tensor of elastic moduli for a homogeneous isotropic material.

The essential and natural boundary conditions are:

u u on C u ; 19a r Á n t on C t ; 19b r Á n i h 0 on C i h i 1; 2; . . . ; m; 19c sr Á n i I t 0 on C i I i 1; 2; . . . ; n; 19d
where n is the unit outward normal to X, and u and t are prescribed displacements and tractions, respectively.

Weak form

Let u be the displacement solution for the stated elastostatic boundary-value problem, with ru the corresponding Cauchy stress tensor. Let u P V be the displacement trial solution, and v P V 0 be any set of kinematically admissible test functions (virtual displacements). The space V H 1 X is the Sobolev space of functions with square-integrable ®rst derivatives in X, and V 0 H 1 0 X is the Sobolev space of functions with square-integrable ®rst derivatives in X and vanishing values on the essential boundary C u . The weak form of the governing equation and associated boundary conditions can be written as Find u P V such that

X ru : ev dX X b Á v dX Ct t Á v dC Vv P V 0 : 20

Discrete system

Consider the Bubnov±Galerkin implementation for the X-FEM in two-dimensional linear elasticity. In the X-FEM, ®nite-dimensional subspaces V h & V and V h 0 & V 0 are used as the approximating trial and test spaces. The weak form for the discrete problem can be stated as:

Find u h P V h & V such that X h ru h : ev h X h b Á v h dX C h t t Á v h dC Vv h P V h 0 & V 0 : 21
In a Bubnov±Galerkin procedure, the trial functions u h as well as the test functions v h are represented as linear combinations of the same shape functions. The trial and test functions are:

u h x I / I xu I J / J xw ux a J ; 22 v h x I / I xv I J / J xw ux b J ; 23
where / I x are the ®nite-element shape functions, ux is the level set function, and wux is the enrichment function for a material interface. On substituting the trial and test functions from Eq. ( 22) in Eq. ( 21), and using the arbitrariness of nodal variations, the following discrete system of linear equations is obtained:

Kd f; 24 where K IJ X h B T I CB J dX; 25a f I C h t /I t dC

Numerical integration

In standard ®nite-element methods, each ®nite element is associated with a material type; however, in the X-FEM, since internal boundaries of geometric entities do not conform to the mesh, such an association does not exist a priori. Using the level set function u, we can consider two dierent computational strategies to perform the numerical integration of the weak form. In the ®rst approach, we establish the association between an element or sub-element (partitioned element) and a material type, and use Gauss quadrature rules in each element to carry out the numerical integration. Alternatively, since the material type of a point x in the domain is known by the sign of ux, a quad-tree or octree sub-division in two-and three dimensions, respectively, in conjunction with Gauss quadrature, Simpson's rule, or other quadrature rules, is a possible candidate for the numerical integration.

In order to implement the ®rst approach in two dimensions, we partition the ®nite elements that intersect the boundary of the hole or inclusion, into triangles. The level set function u is used to obtain the onedimensional segments of the interface, where each segment sub-divides a ®nite element t P T (see Section 4.2.1) into precisely two domains. An algorithm is implemented to partition a ®nite element t P T that is cut by a segment a; b, where a and b lie on the edges of t. The edges of the element t are traversed and if for a particular edge e with nodal connectivity n e I ; n e J , u e I u e J < 0, then the intersection point x p is determined by

x p x e I nx e J À x e I ; n À u e I u e J À u e I : 27
The points a and b are found using the above relations. An illustration of the partitioning is depicted in Fig. 2 for an elliptical interface. In Fig. 2, the discretization of the elliptical interface that is realized by the level set computations is shown, along with the partitioned sub-elements on either side of the interface. If the discretized segment of an interface lies along an edge of a ®nite element, then the element is not partitioned. The above approach is accurate and feasible in two dimensions, since the discretization of the interface into linear segments and the partitioning of the ®nite elements into sub-triangles is readily performed. In three dimensions, however, the interface is represented by the union of triangles, and hence the implementation of a partitioning algorithm for the ®nite elements is non-trivial. In three-dimensional computations, the second approach we allude to has merit. By an octree sub-division of the ®nite elements that intersects the interface, and integrating on a given cell if all its vertices lie on the same side of the interface, an easy-to-implement algorithm emerges.

Numerical results

Three problems in two-dimensional elastostatics are presented to illustrate the accuracy of the new technique. We ®rst solve the problem of an in®nite plate with a circular hole under uniaxial tension, and then study two problems of a bimaterial subjected to pure displacement boundary conditions. Since all problems possess an exact solution, a convergence analysis is performed to study the accuracy and rate of convergence of the proposed method. The error in the energy norm that is used in the analysis is de®ned as

ku À u h k EX X e
À e h T Ce À e h dX 1=2 : 28

In®nite plate with a circular hole

An in®nite plate with a traction-free circular hole under uniaxial tension (r 0 1) along x 1 is considered. The exact solution to this problem is given in [START_REF] Szab O | Finite Element Analysis[END_REF]. In the numerical model, we consider a square domain of edge length L with a circular hole of radius a at its center. We impose the exact tractions on the boundary of the square domain, with appropriate constraints added to remove rigid body modes. The material properties chosen are: Young's modulus E 10 5 and Poisson's ratio m 0:3.

In polar coordinates r; h, the exact stress distribution for r 0 1 is given by

r 11 r; h 1 À a 2 r 2 3 2 cos 2h cos 4h 3 2 a 4 r 4 cos 4h; 29a r 22 r; h À a 2 r 2 1 2 cos 2h À cos 4h À 3 2 a 4 r 4 cos 4h; 29b r 12 r; h À a 2 r 2 1 2 sin 2h sin 4h 3 2 a 4 r 4 sin 4h; 29c
where a is the radius of the circular hole. The displacement components are:

u 1 r; h a 8l r a j 1 cos h 2 a r 1 j cos h cos 3h À 2 a 3 r 3 cos 3h ! ; 30a u 2 r; h a 8l r a j À 3 sin h 2 a r 1 À j sin h sin 3h À 2 a 3 r 3 sin 3h ! ; 30b
where l is the shear modulus and j (Kolosov constant) is de®ned as

j 3 À 4m plane strain; 3 À m 1 m plane stress: V X 31
In the numerical computations, L 2 and plane strain conditions are assumed. A convergence study is carried out using equi-spaced rectangular L h  L h meshes: L h 10; 20; 40; 80; 160 are used in the computations. The analyses are performed for two dierent hole radii: a 0:3 and a 0:4. In Fig. 3, the plot of the error in the energy norm versus the mesh spacing is shown on a log±log plot. The rate of convergence (R) in the energy norm is 0.96, which compares well to the optimal rate of one. In order to assess the accuracy of the use of level sets in the X-FEM for modeling voids, we compare the ®nite-element solution to that obtained by the new technique for a mesh with approximately the same number of unknowns. The radius of the hole a 0:4. The ®nite-element mesh used consists of 3088 three-noded triangular elements with 1569 nodes (hole explicitly modeled by the mesh), and the corresponding mesh for the X-FEM implementation has 3406 three-noded triangular elements and 1704 nodes (Fig. 4). In Fig. 4, the nodes that are removed from the discrete system are indicated by the dark circles. The number of unknowns in the FE computations is 3135, and that in the X-FEM computations is 3103. A contour plot of the normal stress r 11 using the two approaches is presented in Fig. 5. It is seen that the agreement between the results obtained by the two methods is excellent.

Bimaterial bar

Consider a bimaterial bar (X X 1 X 2 ) of length L 2 with material moduli E 1 in X 1 and E 2 in X 2 , where X 1 À1; n and X 2 n; 1. The interface C I is located at x 1 n. For the pure Dirichlet problem with u 1 0 at x 1 À1, u 1 1 at x 1 1, and no body forces, the exact displacement solution assumes the form:

u 1 x 1 1 x 1 a; À1 6 x 1 6 n; 1 E 1 E 2 x 1 À 1a; n 6 x 1 6 1; & 32
where In order to explore the appropriate enrichment function for the material interface, we simulate a one-dimensional bimaterial problem in a two-dimensional domain. Since the exact displacement solution is piecewise linear for the bimaterial bar in one dimension, we seek an enrichment function for the material interface that preserves the equivalence between the discrete X-FEM space and the discrete FEM space in which the interface is explicitly modeled using ®nite elements. By analyzing the numerical results for the bimaterial bar problem, we can generalize the enrichment methodology for general curvilinear interfaces in two dimensions, with an aim towards accuracy and convergence rate that are on par with ®nite elements in which the interface is explicitly modeled using elements.

a E 2 E 2 1 n À E 1 n À 1 : 33
We consider a two-dimensional domain X À1; 1  À1; 1 with the material interface C I located at n  À1; 1. Let Young's modulus and Poisson's ratio in X 1 À1; n  À1; 1 be E 1 1, m 0, and that in X 2 n; 1  À1; 1 be E 2 10, m 0. Numerical results are computed on an equi-spaced 10  10 regular ®nite-element mesh with ®ve dierent locations of the interface: n 0:01; 0:05; 0:1; 0:15; and 0:19. The enrichment functions w juj and w ŵ are considered. In Fig. 6(a), the ®nite-element mesh and the interface are shown for n 0:1, and in Fig. 6(b), the enrichment function ŵ is illustrated. The numerical results for the relative error in the energy norm are presented in Table 1. On using w juj as the enrichment function, a linear approximation is not obtained in elements that are adjacent to those intersected by the interface. This adversely aects the energy error norm results. However, by using the enrichment function ŵ shown in Fig. 6(b), a linear approximation is obtained in X which leads to very accurate results. The numerical results for this model problem indicate that the enrichment function away from the material interface needs to be a constant in order to pass the bimaterial patch test. In Fig. 7, a body composed of two dierent materials is shown. The material constants are constant in X 1 and X 2 , but there is a discontinuity in the material constants across the interface C 1 (r a). The Lam e constants in X 1 are chosen as: k 1 l 1 0:4, and those in X 2 are: k 2 5:7692, l 2 3:8461. These correspond to E 1 1, m 1 0:25, and E 2 10, m 2 0:3. We impose the linear displacement ®eld: u 1 x 1 , u 2 x 2 (u r r, u h 0) on the boundary C By considering displacement and traction continuity across the interface, the exact displacement solution can be written as

u r r 1 À b 2 a 2 a b 2 a 2 h i r; 0 6 r 6 a; r À b 2 r a b 2 r ; a 6 r 6 b; V b b X 35a u h 0; 35b where a k 1 l 1 l 2 b 2 k 2 l 2 a 2 k 1 l 1 b 2 À a 2 l 2 b 2 : 35c
The radial (e rr ) and hoop (e hh ) strains are given by:

e rr r 1 À b 2
a 2 a b 2 a 2 ; 0 6 r 6 a; where the appropriate Lam e constants are to be used in the evaluation of the normal stresses. The shear components of the stress and strain tensors are zero.

1 b 2 r 2 a À b 2 r 2 ; a 6 r 6 b; V b b X 36a
In the numerical model, we consider a square plate (L  L, L 2) with a circular inclusion (X 1 ) of radius a 0:4. On the boundary of the plate, the exact tractions using Eq. (37a) are imposed (a 0:4 and b 2:0 are chosen) with appropriate displacement constraints added to remove rigid body modes. A convergence study is conducted using equi-spaced rectangular L h  L h meshes: L h 10; 20; 40; 80. The ®nite-element computations are carried out using three-noded triangular elements.

The convergence analysis is carried out for three cases using the X-FEM: (a) no enrichment; (b) enrichment function w juj; and (c) enrichment function w w 1 , where w 1 is the modi®ed level set function based on the smoothing procedure described in Section 4.2.1. The motivation for the smoothing procedure presented in Section 4.2.1 stems from the numerical experiments conducted in Section 7.2. The enrichment functions for the 20 Â 20 mesh is shown in Fig. 8. In all the above cases, the weak form is integrated appropriately by partitioning the elements that are intersected by the interface. The degrees of freedom of the meshes used in the FEM and in the X-FEM computations are presented in Table 2. The numerical results for the convergence in energy norm are shown in Fig. 9. The accuracy and rate of convergence of the X-FEM with the smoothed form of the enrichment are signi®cantly better than the other two cases. The rate (R 0:91) using the enrichment function w 1 is also on par with the optimal rate of unity which is obtained using ®nite elements with the interface truly modeled by the mesh.

Even though the enrichment function w 1 used in this paper is appropriate from an accuracy and convergence viewpoint, a more general enrichment methodology which would be valid for a wider class of interface problems is attractive. To this end, one can consider the enrichment nodal data w I as also unknowns in the statement of the principle of virtual work (variational form). The nodal data w I are known on nodes that are in the connectivity of the elements that intersect the interface: the modulus of the level set function (juj) is an appropriate choice for this purpose. On elements that are one layer on either side (outer and inner) of these elements, the w I are to be determined by the numerical procedure. As an initial guess, the w I can be set using the smoothing procedure given in Eqs. ( 15)± [START_REF] Sukumar | Extended ®nite element method for three-dimensional crack modeling[END_REF]. Thereafter, an iterative scheme between the solution of the minimization problem with respect to fu I ; a I g and that with respect to w I would provide the minimization of the global problem. This would result in an optimum choice for w I , and an appropriate form for the enrichment function w which is best suited for a given problem.

Conclusions

A methodology to model arbitrary holes and material interfaces (inclusions) without meshing the internal boundaries was proposed. The numerical method couples the level set method [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton±Jacobi formulations[END_REF] to the X-FEM [START_REF] Mo Es | A ®nite element method for crack growth without remeshing[END_REF]. In the X-FEM, the ®nite-element space is enriched by adding additional functions to the approximation using the notion of partition of unity. The level set function u is used to represent holes and inclusions, and in addition, u is used to develop the local enrichment for material interfaces to model inclusions.

Numerical examples in two-dimensional linear elastostatics were presented to test the accuracy of the new technique. For the plate with a hole problem, a convergence study was performed for two dierent radii of the hole. Optimal rate of convergence in energy norm Oh was obtained using the new technique. A bimaterial boundary-value problem with pure displacement boundary conditions was also considered. Issues pertaining to appropriate choice of enrichment functions for material interfaces were addressed, and a suitable enrichment function was derived that yielded accurate results and near optimal rate of convergence for the two-dimensional bimaterial problem. This study demonstrated that by incorporating level set functions in the X-FEM, we obtain a powerful and accurate numerical technique for the modeling of voids and inclusions without meshing the internal boundaries. The new technique is simple and versatile, and the results indicate its potential and promise in the modeling of interfacial phenomena in solid and ¯uid mechanics.

Fig. 1 .

 1 Fig. 1. Level set function for a hexagonal interface: (a) mesh; (b) level set function.

X

  h /I b dX; 25b where /I / I for a ®nite-element displacement degree of freedom, and /I / I w for an enriched degree of freedom. In the above equations, C is the constitutive matrix for an isotropic linear elastic material, and the matrix B I is given by B I

Fig. 2 .

 2 Fig. 2. Partitioning and material typing of sub-elements.

Fig. 4 .

 4 Fig. 4. Mesh used in the X-FEM computations.

Fig. 3 .

 3 Fig. 3. Rate of convergence in energy norm for the plate with a hole problem.

Fig. 5 .

 5 Fig. 5. Contour plot of the normal stress r 11 : (a) FEM; (b) X-FEM.

Fig. 6 .

 6 Fig. 6. Bimaterial bar problem: (a) mesh; (b) enrichment function ŵ (n 0:1).

  2 r b. Navier's equation in polar coordinates reduces to

Fig. 7 .

 7 Fig. 7. Bimaterial boundary-value problem.

e hh r 1 À b 2 a 2 a b 2 a 2 ; 0 6 r 6 a; 1 À b 2 r 2 a b 2 r 2 ; a 6 r 6

 12126 Fig. 8. X-FEM enrichment for material interfaces: (a) w juj; (b) w w 1 .

Table 1

 1 Relative error in the energy norm for the bimaterial bar problem

	Interface location n	Enrichment
		juj	ŵ
	0.01	8:3 Â 10 À2	3:0 Â 10 À8
	0.05	1:6 Â 10 À1	2:8 Â 10 À8
	0.10	1:8 Â 10 À1	2:1 Â 10 À8
	0.15	1:8 Â 10 À1	3:8 Â 10 À8
	0.19	1:6 Â 10 À1	3:6 Â 10 À8

Table 2

 2 Degrees of freedom for the meshes used for the bimaterial problem

	Mesh h	FEM	X-FEM (no enrichment)	X-FEM (enrichment)
	0.2	277	239	287
	0.1	915	879	991
	0.05	3559	3359	3599
	0.025	13,693	13,119	13,615

Fig. 9. Rate of convergence in energy norm for the bimaterial problem.

Acknowledgements

The authors are grateful for the research support of the National Science Foundation through contract CMS±9732319, and to the Oce of Naval Research and the Army Research Oce.


David Gerault

Marine Minier

Christine Solnon

Using Constraint Programming to solve a Cryptanalytic Problem *

We describe Constraint Programming (CP) models to solve a cryptanalytic problem: the chosen key differential attack against the standard block cipher AES. We show that CP solvers are able to solve these problems quicker than dedicated cryptanalysis tools, and we prove that a solution claimed to be optimal in two recent cryptanalysis papers is not optimal by providing a better solution.

Introduction

Since 2001, AES (Advanced Encryption Standard) is the encryption standard for block ciphers [FIPS 197, 2001]. It guarantees communication confidentiality by using a secret key K to encode an original plaintext X into a ciphertext AES K (X), in such a way that the ciphertext can further be decoded into the original one using the same key, i.e., X = AES -1 K (AES K (X)). Cryptanalysis aims at testing whether confidentiality is actually guaranteed. In particular, differential cryptanalysis [START_REF] Biham | Differential cryptoanalysis of feal and n-hash[END_REF] evaluates whether it is possible to find the key within a reasonable number of trials by considering plaintext pairs (X, X ) and studying the propagation of the initial difference X ⊕ X between X and X while going through the ciphering process (where ⊕ is the xor operator). Today, differential cryptanalysis is public knowledge, and block ciphers such as AES have proven bounds against differential attacks. Hence, [START_REF] Biham | New types of cryptoanalytic attacks using related keys (extended abstract)[END_REF] proposed a new type of attack called related-key attack that allows an attacker to inject differences not only between the plaintexts X and X but also between the keys K and K (even if the secret key K stays unknown from the attacker).

To mount related-key attacks, the cryptanalist must find optimal related-key differentials, i.e., a plaintext difference δX and a key difference δK that maximize the probability that, for randomly chosen plaintext X and key K, the difference X ⊕ δX in the input plaintexts becomes the difference AES K (X) ⊕ AES K⊕δK (X ) in the output ciphertexts. Finding the optimal related-key differentials for AES is a highly combinatorial problem that hardly scales. Two main approaches have been proposed to solve this problem: a graph traversal approach [START_REF] Fouque | Structural evaluation of aes and chosenkey distinguisher of 9-round aes-128[END_REF], and a Branch & Bound approach [Biryukov and Nikolic, 2010]. The approach of [START_REF] Fouque | Structural evaluation of aes and chosenkey distinguisher of 9-round aes-128[END_REF] requires about 60 GB of memory when the key has 128 bits, and it has not been extended to larger keys. The approach of [Biryukov and Nikolic, 2010] only takes several megabytes of memory, but it requires several days of computation when the key has 128 bits, and several weeks when the key has 192 bits.

During the process of designing new ciphers, this search generally needs to be performed several times, so it is desirable that it can be done rather quickly. Another point that should not be neglected is the time needed to design and implement these approaches: To ensure that the computation is completed within a "reasonable" amount of time, it is necessary to reduce the branching by introducing clever reasoning. Of course, this hard task is also likely to introduce bugs, and checking the correctness or the optimality of the computed solutions may not be so easy. Finally, reproducibility may also be an issue. Other researchers may want to adapt these algorithms to other problems, with some common features but also some differences, and this may again be very difficult and time-consuming.

In [START_REF] Gerault | [END_REF], we proposed to use Constraint Programming (CP) to solve this problem. When using CP to solve a problem, one simply has to model the problem by means of constraints. Then, this model is solved by generic solvers usually based on a Branch & Propagate approach: The search space is explored by building a search tree, and constraints are propagated at each tree node in order to prune branches. CP opens new perspectives for this kind of cryptanalysis problems. First, it is very competitive with dedicated approaches: When the key has 128 bits, Chuffed [Chu and Stuckey, 2014] is able to find optimal solutions in less than one hour. Second, the CP model is easier to check or re-use than a full program that not only describes the problem to solve, but also how to solve it. Actually, CP allowed us to prove that a solution claimed to be optimal in [Biryukov and Nikolic, 2010;[START_REF] Fouque | Structural evaluation of aes and chosenkey distinguisher of 9-round aes-128[END_REF] is not optimal by providing a better solution.

2 Problem Statement AES block cipher. AES ciphers blocks of length n = 128 bits, and each block is a 4 × 4 matrix of bytes. The length of keys is l ∈ {128, 192, 256}. In this paper, we only consider Operations applied at each round i ∈ [0, r -1]: AES is an iterative process, and we note X i the ciphertext at the beginning of round i ∈ [0, r]. Each round is composed of the following operations, as displayed in Fig. 1:

Key K = K 0 (4×4 bytes) KS SK i [j][3] = S(K i [j][3]) KS Subkey K i+1 X i Plaintext X (4×4 bytes) ARK i = 0 ARK i = r i < r S SX i = S(X i ) SR SR(SX i ) M C Y i = M C(SR(SX i )) Xr SB Ciphertext SXr = AES K (X)
• SubBytes S. S is a non-linear permutation which is applied on each byte of X i separately, i.e., ∀j, k ∈ [0, 3], X i [j][k] is replaced by S(X i [j][k]
), according to a lookup table. We note SX i = S(X i ).

• ShiftRows SR. SR is a linear mapping that rotates on the left by 1 byte position (resp. 2 and 3 byte positions) the second row (resp. third and fourth rows) of SX i .

• MixColumns M C. M C is a linear mapping that multiplies each column of SR(SX i ) by a 4 × 4 fixed matrix chosen for its good properties of diffusion [START_REF] Daemen | [END_REF]. In particular, it has the Maximum Distance Separable (MDS) property: For each column, the total number of bytes which are different from 0, before and after applying M C, is either equal to 0 or strictly greater than 4. We note

Y i = M C(SR(SX i )).
• KeySchedule KS. The subkey at round 0 is the initial key, i.e.,

K 0 = K. For each round i ∈ [0, r -1], the subkey K i+1 is generated from K i by applying KS. It first replaces each byte K i [j][3] of the last column by S(K i [j][3]
) (where S is the SubBytes operator), and we note

SK i [j][3] = S(K i [j][3]).
Then, each column of K i+1 is obtained by performing a xor operation between bytes coming from K i , SK i , or K i+1 .

• AddRoundKey ARK. ARK performs a xor operation between bytes of Y i and subkey K i+1 to obtain X i+1 .

Let us note Bytes the set of all bytes (for all rounds), i.e., Mounting attacks requires finding a related-key differential characteristic, i.e. a plaintext difference δX = X ⊕ X and a key difference δK = K ⊕ K , such that δX becomes δSX r after r rounds with a probability as hight as possible. This can be achieved by tracking the propagation of the initial differences through the cipher. Once such an optimal differential characteristic is found, the cryptanalist asks for the encryption of messages satisfying the input difference. When this difference propagates as expected, he can infer informations leading to a recovery of the secret key K with less workload than exhaustive search. The difficulty of recovering K decreases as the probability of the used differential characteristic increases. The AES operators SR, M C, ARK, and KS are linear, i.e., they propagate differences in a deterministic way (with probability 1). However, the S operator is not linear: Given a byte difference B ⊕ B = δB, the probability that δB becomes S(B) ⊕ S(B ) = δSB is P r(δB → δSB). It is equal to 1 if δB = 0 (i.e., B = B ). However, if δB = 0, then P r(δB → δSB) ∈ { 2 256 , 4 256 }. The probability that δX becomes δSX r is equal to the product of all P r(δB → δSB) such that δB (resp. δSB) is a byte difference before (resp. after) passing through the S operator when ciphering X with K and X with K , i.e., the product of all P r(δX i

Bytes = {X[j][k], X i [j][k], SX i [j][k], Y i [j][k], K i [j][k], SK i [j][3] | i ∈ [0, r], j, k ∈ [0, 3]}
[j][k] → δSX i [j][k]) and all P r(δK i [j][3] → δSK i [j][3]) with i ∈ [0, r] and j, k ∈ [0, 3].
The goal is to find δX and δK that maximize this probability.

Two step solving process. To find δX and δK, we search for the values of all differential bytes in diffBytes. Both [Biryukov and Nikolic, 2010] and [START_REF] Fouque | Structural evaluation of aes and chosenkey distinguisher of 9-round aes-128[END_REF] propose to solve this problem in two steps. In Step 1, a Boolean variable ∆B is associated with every differential byte δB ∈ diffBytes such that ∆B = 0 ⇔ δB = 0 and , 255]. The goal of Step 1 is to find a Boolean solution that assigns values to Boolean variables such that the AES transformation rules are satisfied. During this first step, the SubBytes operation S is not considered. Indeed, it does not introduce nor remove differences. Therefore, we have

∆B = 1 ⇔ δB ∈ [1
∆X i [j][k] = ∆SX i [j][k] and ∆K i [j][3] = ∆SK i [j][3].
As we search for a solution with maximal probability, the goal of Step 1 is to search for a Boolean solution which minimizes the number of variables ∆X i [j][k] and ∆K i [j][3] which are set to 1.

In Step 2, the Boolean solution is transformed into a byte solution: For each differential byte δB ∈ diffBytes, if the corresponding Boolean variable ∆B is assigned to 0, then δB is also assigned to 0; otherwise, we search for a byte value in [1,255] to be assigned to δB such that the AES transformation rules are satisfied and the probability is maximized. Note that some Boolean solutions may not be transformable into byte solutions. These Boolean solutions are said to be byte-inconsistent.

Basic CP Model for Step 1

A first CP model for Step 1 may be derived from the AES transformation rules in a rather straightforward way. A CP model is defined by a set of variables, such that each variable x has a domain D(x), and a set of constraints, i.e., relations that restrict the values that may be simultaneously assigned to the variables.

Variables. For each differential byte δB ∈ diffBytes, we define a Boolean variable ∆B whose domain is D(∆B) = {0, 1}: it is assigned to 0 if δB = 0, and to 1 otherwise. XOR constraint. We first define a XOR constraint for ARK and KS. Let us consider three differential bytes δA, δB and δC such that δA ⊕ δB = δC. If δA = δB = 0, then δC = 0. If (δA = 0 and δB = 0) or (δA = 0 and δB = 0) then δC = 0. However, if δA = 0 and δB = 0, then we cannot know if δC is equal to 0 or not: This depends on whether δA = δB or not. When abstracting differential bytes δA, δB and δC with Boolean variables ∆A, ∆B and ∆C (which only model the fact that there is a difference or not), we obtain the following definition of the XOR constraint: XOR(∆A, ∆B, ∆C) ⇔ ∆A + ∆B + ∆C = 1.

AddRoundKey and KeySchedule constraints. Both ARK and KS are modeled with XOR constraints between ∆X i , ∆Y i , and ∆K i variables (for ARK), and between ∆K i and ∆SK i variables (for KS). ShiftRows and MixColumns. SR simply shift variables. The MDS property of M C is ensured by posting a constraint on the sum of some variables of ∆X i and ∆Y i variables, which must belong to the set {0, 5, 6, 7, 8}.

Objective function. We introduce an integer variable obj Step1 that must be minimized, and we post an equality constraint between obj Step1 and the sum of Boolean variables on which a non linear S operation is performed (all variables ∆X i [j][k] and ∆K i [j][3] with i ∈ [0, r] and j, k ∈ [0, 3]).

Let v be the optimal value of obj Step1 . It may happen that none of the Boolean solutions with obj Step1 = v is byte-consistent, or that the maximal probability p of Boolean solutions with obj Step1 = v is such that it is possible to have a better probability with a larger value for obj Step1 (i.e., p < ( 4 256 ) v+1 ). In this case, we need to search for new Boolean solutions, such that obj Step1 is minimal while being strictly greater than v. This is done by adding the constraint obj Step1 > v before solving again Step 1.

Limitations of the basic CP model. This basic CP model is complete, i.e., for any solution at the byte level (on δ variables), there exists a solution at the Boolean level (on ∆ variables). However, preliminary experiments have shown us that there is a huge number of Boolean solutions which are byte inconsistent. For example, when the number of rounds is r = 4, the optimal cost is obj Step1 = 11, and there are more than 90 millions of Boolean solutions with obj Step1 = 11. However, none of these solutions is byte-consistent. In this case, most of the Step 1 solving time is spent at generating useless Boolean solutions which are discarded in Step 2.

Additional constraints for step 1

When reasoning at the Boolean level, many solutions are not byte-consistent because constraints at the byte level have been ignored. To propagate some properties at the byte level, we introduce new variables and constraints that model equality relations between differential bytes. New equality variables. For each couple of differential bytes (δA, δB) ∈ diffBytes 2 , we introduce a Boolean equality variable EQ δA,δB which is equal to 1 if δA = δB, and to 0 otherwise. These variables are constrained to define an equivalence relation by adding a symmetry constraint (EQ δA,δB = EQ δ B ,δ A ) and a transitivity constraint (if EQ δA,δB = EQ δB,δC = 1 then EQ δ A ,δ C = 1). Also, EQ variables are related to ∆ variables by adding the constraints:

(EQ δA,δB = 1) ⇒ (∆A = ∆B) EQ δA,δB + ∆A + ∆B = 0

Revisiting the XOR constraint. When defining the constraint XOR(∆A, ∆B, ∆C), if ∆A = ∆B = 1, then we cannot know if ∆C is equal to 0 or 1. However, whenever ∆C = 0 (resp. ∆C = 1), we know for sure that the corresponding byte δC is equal to 0 (resp. different from 0), meaning that the two bytes δA and δB are equal (resp. different), i.e., that EQ δA,δB = 1 (resp. EQ δA,δB = 0). The same reasoning may be done for ∆A and ∆B because

(δA ⊕ δB = δC) ⇔ (δB ⊕ δC = δA) ⇔ (δA ⊕ δC = δB).
Therefore, we redefine the XOR constraint as follows: XOR(∆A, ∆B, ∆C) ⇔ ((∆A + ∆B + ∆C = 1) ∧ (EQ δA,δB = 1 -∆C) ∧ (EQ δA,δC = 1 -∆B) ∧ (EQ δB,δC = 1 -∆A))

Propagation of MDS at the byte level. The MDS property ensures that, for each column, the total number of bytes which are different from 0, before and after applying M C, is either equal to 0 or strictly greater than 4. This property also holds for any xor difference between two different columns of X and Y matrices. To propagate this property, for each pair of columns in X and Y matrices, we add a constraint on the sum of equality variables between bytes of these columns.

Constraints derived from KS. The KeySchedule mainly performs xor and S operations. As a consequence, each byte δK i [j][k] may be expressed as a xor between bytes of the original key difference δK 0 , and bytes of δSK i-1 (which are differences of key bytes that have passed through S during the previous round). Hence, for each byte δK i [j][k], we precompute the set V (i, j, k) such that V (i, j, k) only contains bytes of δK 0 and δSK i-1 and

δK i [j][k] = δA∈V (i,j,k) δA.
For each set V (i, j, k), we introduce a set variable V 1 (i, j, k) which is constrained to contain the subset of V (i, j, k) corresponding to the Boolean variables equal to 1. We use these set variables to infer that two differential key bytes that have the same V 1 set are equal. Also, if V 1 (i, j, k) is empty (resp. contains one or two elements), we infer that ∆K i [j][k] is equal to 0 (resp. a variable, or a xor between 2 variables).

CP model for Step 2

Given a Boolean solution for Step 1, Step 2 aims at searching for the byte-consistent solution with the highest probability (or proving that there is no byte-consistent solution). Hence, for each Boolean variable ∆A of Step 1, we define an integer variable δA whose domain depends on the value of ∆A: If ∆A = 0, then D(δA) = {0} (i.e., δA is also assigned to 0); otherwise, D(δA) = [1,255]. As we look for a byte-consistent solution with maximal probability, we also add an integer variable P A for each byte A that passes through the S operator, i.e., X i [j][k] and K i [j][3]: This variable corresponds to the base 2 logarithm of the probability Pr(δA → δS A ) of obtaining the S output difference δS A when the S input difference is δA. The domain of these variables depends on the value of ∆A in the Step 1 solution: If ∆A = 0, then Pr(0 → 0) = 1 and therefore D(P A ) = {0}; otherwise, Pr(δA → δS A ) ∈ { 2 256 , 4 256 } and D(P A ) = {-7, -6}.

The constraints basically follow the AES operations to relate variables, as described for Step 1, but consider the definition of the operations at the byte level, instead of the Boolean level. A main difference is that the SubBytes operation, which has no effect at the Boolean level, must be modeled at the byte level. This is done thanks to a ternary table constraint which extensively lists all triples (X, Y, P ) such that there exist two bytes B 1 and B 2 whose difference before and after passing through S is equal to X and Y , respectively, and such that P is the base 2 logarithm of the probability of this transformation. The goal is to find a byte-consistent solution with maximal differential probability. As we consider logarithms, this amounts to searching for a solution that maximizes the sum of all P A variables.

Results and Conclusion

The CP model for Step 1 was implemented with the MiniZinc modeling language [Nethercote et al., 2007], and we compared three CP solvers: Gecode [Gecode Team, 2006], Choco 4 [START_REF] Prudhomme | An introduction to choco 3.0: an open source java constraint programming library[END_REF], and Chuffed [Chu and Stuckey, 2014]. The best results were obtained with Chuffed. The basic CP model (described in Section 3) does not scale well because it generates a huge number of Boolean solutions which are not byte-consistent. When adding the additional constraints described in Section 4, most of these byte-inconsistent solutions are filtered out, and Chuffed is able to solve all instances in less than one hour. The CP model for Step 2 was implemented with Choco 3. It solves all instances (i.e., finds an optimal byte-consistent solution for each Boolean solution of Step 1 when it exists) in a few seconds. As a conclusion, CP solvers are much faster than the Branch & Bound approach of [Biryukov and Nikolic, 2010], which needs several days to solve these instances. It is also faster and much less memory consuming than the approach of [START_REF] Fouque | Structural evaluation of aes and chosenkey distinguisher of 9-round aes-128[END_REF], that needs 60GB and 30 minutes on a 12-core computer to pre-compute the graph.

New results for differential cryptanalysis. For r = 4 rounds, we have found a byte-consistent solution with obj Step1 = 12 and a probability equal to 2 -79 . This solution is better than the solution claimed to be optimal in [Biryukov and Nikolic, 2010] and [START_REF] Fouque | Structural evaluation of aes and chosenkey distinguisher of 9-round aes-128[END_REF]: In these papers, the authors say that the best byte-consistent solution has obj Step1 = 13, and a probability equal to 2 -81 .

We have shown how to extend our CP models to AES with longer keys (such that l ∈ {192, 256}) in [Gérault et al., 2017]. These models allowed us to find optimal solutions for all possible instances of AES with l ∈ {128, 192, 256} in less than 35 hours for Step 1, and in less than 6 minutes for Step 2. This is a clear improvement with respect to existing work as the approach of [START_REF] Fouque | Structural evaluation of aes and chosenkey distinguisher of 9-round aes-128[END_REF] cannot be extended to l > 128 due to its memory complexity, and the approach of [Biryukov and Nikolic, 2010] needs several weeks to solve instances with l = 192. Furthermore, we have shown that the solution proposed in [Biryukov and Nikolic, 2010] for l = 192 and r = 11 is inconsistent. We have also found better solutions when l = 256, and we have computed the actual optimal solution for AES with l = 256. Its probability is 2 -146 (instead of 2 -154 for the solution of [Biryukov et al., 2009]). Using this solution, we improved the relatedkey distinguisher and the basic related-key differential attack on the full AES-256 by a factor 2 6 and the q-multicollisions by a factor 2 (see [Gérault et al., 2017] for more details).

These cryptanalysis problems open new and exciting challenges for the CP community. In particular, these problems are not easy to model. More precisely, naive CP models such as the one described in Section 3 may not scale well. The introduction of equality constraints at the byte level drastically improves the solving process, but these constraints are not straightforward to find and implement. Hence, a challenge is to define new CP frameworks, dedicated to this kind of cryptanalysis problems, in order to ease the development of efficient CP models for these problems.

  Optimal related-key differentials. Let us note δB = B ⊕ B the difference between two bytes B and B , and diffBytes = {δB | B ∈ Bytes} the set of all differential bytes.

  Figure1: AES ciphering process. Each 4 × 4 array represents a group of 16 bytes. Before the first round, ARK is applied on X and K to obtain X 0 . Then, for each round i ∈ [0, r -1], S is applied on X i to obtain SX i , SR and M C are applied on SX i to obtain Y i , KS is applied on K i to obtain K i+1 (and during KS, S is applied onK i [j][3] to obtain SK i [j][3], ∀j ∈ [0, 3]),and ARK is applied on K i+1 and Y i to obtain X i+1 . The ciphertext SX r is obtained by applying SB on X r .

	keys of length l = 128, and keys are 4 × 4 matrices of bytes.
	Given a 4 × 4 matrix of bytes M , we note M [j][k] the byte at
	row j ∈ [0, 3] and column k ∈ [0, 3].

* This research was conducted with the support of the FEDER program of 2014-2020 and the region council of Auvergne