Aymeric Pionteck

Baptiste Pierrat

Sébastien Gorges

Jean-Noël Albertini

Stéphane Avril

A Fast Method of Virtual Stent Graft Deployment for Computer Assisted EVAR

In this paper we introduce a new method simulating stent graft deployment for assisting endovascular repair of abdominal aortic aneurysms. The method relies on intraoperative images coupled with mechanical models. A multi-step algorithm has been developed to increase the reliability of simulations. The first step predicts the position of the stent graft within the aorta. The second step is an axisymmetric geometric reconstruction of each individual stent. The third step minimizes the rotation of each stent around its main axis. Finally, the last step virtually deploys each stent within a deployment box extracted from the preoperative CT scan. A proof of concept is performed on a patient. The accuracy is compatible with the clinical threshold of 3 mm: the average distance between target and simulated stents is 1.73 ± 0.37 mm. Fenestrations of the stent-graft are reconstructed with a maximum error of less than 2.5 mm, which is enables a secure catheterization of secondary arteries. In summary, the method is able to assist EVAR practitioners by providing all necessary information for a fast and accurate stent graft positioning, combining intraoperative data and a mechanical model in a very low cost framework.

I. INTRODUCTION

Abdominal aortic aneurysm (AAA) is a frequent asymptomatic pathology that results in abnormal local deformation of the aorta. Each year, aneurysm ruptures are responsible for 10,000 deaths in the United States [START_REF] Dua | Epidemiology of aortic aneurysm repair in the United States from 2000 to 2010[END_REF]. Clinical monitoring of the evolution of the aneurysmal sac diameter is used to decide whether an intervention is necessary [START_REF] Dua | Epidemiology of aortic aneurysm repair in the United States from 2000 to 2010[END_REF], [START_REF] Sajid | Endovascular Aortic Aneurysm Repair (EVAR) Has Significantly Lower Perioperative Mortality in Comparison to Open Repair: A Systematic Review[END_REF]. Two options are available: conventional open surgery or endovascular surgery (EVAR). Endovascular surgery is associated with a lower mortality rate (1.5%) than open surgery (4.6%), although long-term mortality is similar [START_REF] Greenhalgh | Endovascular repair of aortic aneurysm in patients physically ineligible for open repair[END_REF], [START_REF] Macdonald | Towards safer carotid artery stenting: A scoring system for anatomic suitability[END_REF].

During EVAR, the surgeon first makes a small incision in the groin to reach the femoral artery. From this incision, tools are introduced to position the stent graft (SG) launcher within the aneurysm. Then the SG is progressively deployed. The success of the intervention depends on the precise positioning of the SG in the artery. In some cases, a fenestrated SG is required if the aneurysm extends beyond the ostia of the renal arteries. In this case, the fenestrations of the SG must be positioned precisely in front of the renal ostium whose diameter is about 5-7mm. This phase is delicate but essential to allow the catheterization of the secondary arteries and avoid occlusions and post-operative complications [START_REF] Albertini | Aorfix Stent Graft for Abdominal Aortic Aneurysms Reduces the Risk of Proximal Type 1 Endoleak in Angulated Necks: Bench-Test Study[END_REF], [START_REF] Albertini | Pathophysiology of Proximal Perigraft Endoleak Following Endovascular Repair of Abdominal Aortic Aneurysms: a Study Using a Flow Model[END_REF] [START_REF] Carroccio | Predicting iliac limb occlusions after bifurcated aortic stent grafting: Anatomic and device-related causes[END_REF], [START_REF] Cochennec | Limb Graft Occlusion Following EVAR: Clinical Pattern, Outcomes and Predictive Factors of Occurrence[END_REF]. The lack of 3D information obliges the surgeon to perform a mental reconstruction of the scene using several images with different incidence angles, which considerably increases the duration of the procedure, as well as the time of exposure to X-rays and the volume of injected contrast products. Therefore, the virtual 3D representation of the tool location and particularly the SG in the aorta is a valuable aid to the surgeon. This would reduce the surgery time and the number of X-ray images required. Moreover, it would reduce the number of postoperative complications, most often related to inaccurate SG positioning. Latest generation systems enable the acquisition of 3D images of the tool during the intervention [START_REF] Akpek | Three-dimensional imaging and cone beam volume CT in C-arm angiography with flat panel detector[END_REF], [START_REF] Fagan | 3-D guidance of complex pulmonary artery stent placement using reconstructed rotational angiography with live overlay[END_REF]. Recently, efforts have been made to use biplanar fluoroscopic acquisitions to reconstruct the 3D shape of the device [START_REF] Baert | Three-Dimensional Guide-Wire Reconstruction from Biplane Image Sequences for Integrated Display in 3-D Vasculature[END_REF]- [START_REF] Wagner | 4D interventional device reconstruction from biplane fluoroscopy[END_REF]. However, all these methods are based on expensive equipment which are not commonplace in all hospitals, usually equipped with simple mobile C-arms. Another solution should therefore be available to obtain a three-dimensional representation of the inserted SG at low cost. Modelling and numerical simulation of SG deployment then appears as essential.

First studies on numerical simulation of SG deployment were based on finite element analyses to study the mechanics of stents and to simulate their deployment in arteries [START_REF] Mortier | A Novel Simulation Strategy for Stent Insertion and Deployment in Curved Coronary Bifurcations: Comparison of Three Drug-Eluting Stents[END_REF]- [START_REF] Auricchio | Carotid artery stenting simulation: From patient-specific images to finite element analysis[END_REF], integrating different types of constitutive behavior for the different materials of SGs [START_REF] Bock | Filling the void: A coalescent numerical and experimental technique to determine aortic stent graft mechanics[END_REF], [START_REF] Kleinstreuer | Fluid-Structure Interaction Analyses of Stented Abdominal Aortic Aneurysms[END_REF], [START_REF] Auricchio | Patient-specific aortic endografting simulation: From diagnosis to prediction[END_REF].

Perrin et al. [START_REF] Perrin | Patient-specific numerical simulation of stent-graft deployment: Validation on three clinical cases[END_REF]- [START_REF] Perrin | Patient-specific simulation of endovascular repair surgery with tortuous aneurysms requiring flexible stent-grafts[END_REF] developed a preoperative planning tool to predict the postoperative position of the SG from patient-specific models. Although essential for preoperative planning, these studies have two important limitations with regard to their use as real time assistance for the practitioner: (i) inappropriately long computation time, and (ii) lack of update from intraoperative images. Some studies have focused on reducing the computation time and developed algorithms to simulate stent deployment in "real time". They often rely on simplifications such as modeling vessels as generic tubes, on which the stent armatures are then mapped [START_REF] Appanaboyina | Computational fluid dynamics of stented intracranial aneurysms using adaptive embedded unstructured grids[END_REF]- [START_REF] Flórez-Valencia | 3D graphical models for vascular-stent pose simulation To cite this version : HAL Id : hal-00682926[END_REF]. The Fast Virtual Stenting (FVS) technique was proposed by Larrabide et al [START_REF] Larrabide | Fast Virtual Stenting with Deformable Meshes: Application to Intracranial Aneurysms[END_REF]. This technique is based on constrained deformable simple models and can virtually model stent deployment in vessel and aneurysm models. The FVS technique was tested and compared with experimental results [START_REF] Larrabide | Fast virtual deployment of self-expandable stents: Method and in vitro evaluation for intracranial aneurysmal stenting[END_REF] and with finite element models [START_REF] Flore | Stent deployment in aneurysmatic cerebral vessels: Assessment and quantification of the differences between Fast Virtual Stenting and Finite Element Analysis[END_REF] [START_REF] Bernardini | Deployment of self-expandable stents in aneurysmatic cerebral vessels: Comparison of different computational approaches for interventional planning[END_REF]. Alternative methods have been proposed, based on mass-spring models [START_REF] Spranger | Which spring is the best? Comparison methods for virtual stenting[END_REF], [START_REF] Spranger | Comparison and calibration of a real-time virtual stenting algorithm using Finite Element Analysis and Genetic Algorithms[END_REF] or on active contours [START_REF] Zhong | Fast Virtual Stenting with Active Contour Models in Intracranical Aneurysm[END_REF]. Although efficient and fast, most of these models are based on simplified mechanics and can be challenged by complex vascular geometries. In addition, applications focus on preoperative planning, as none of the work mentioned above considered intraoperative images.

A small number of studies integrated information from intraoperative images. For example, Demirci et al. [START_REF] Demirci | 3D stent recovery from one x-ray projection[END_REF] proposed an algorithm to automatically match a 3D model of the SG with an intraoperative 2D image of its structure. Zhou et al. [START_REF] Zhou | Real-Time 3-D Shape Instantiation From Single Fluoroscopy Projection for Fenestrated Stent Graft Deployment[END_REF] , and Zheng et al. [START_REF] Zheng | Real-time 3D Shape Instantiation for Partially-deployed Stent Segment from a Single 2D Fluoroscopic Image in Fenestrated Endovascular Aortic Repair[END_REF] introduced a real-time framework to generate the 3D shape of a fenestrated SG from a single 2D fluoroscopic image and position of added radioopaque markers. These methods have reduced computation times and can accurately represent the deployment of SGs in simple geometries. However, more complex cases cannot be addressed without the use of a mechanical model.

To our best knowledge, no studies have ever combined these different aspects into a single method. Achieving this combination is the objective of the present work, in order to propose a method that can assist EVAR practitioners by providing all necessary information for a fast and accurate SG positioning.

The details of the method are given in this book chapter, first introducing the global algorithm, then describing each step and finally showing a proof of concept for a patient case.

II. METHODS

The global algorithm of the method is summarized in Figure 1. The input data are the 2D intraoperative images from a mobile C-arm and the 3D geometry of the aorta obtained from a preoperative CT scan. The algorithm is divided into four main steps. The first two steps can be combined into a single stage called Stage 1. This stage is essential for the following steps but may reach insufficient accuracy, hence the possible following Stage 2. During the first step of Stage1, barycenters of each stent are positioned in 3D using a FEM model of the SG in the aorta. Then, the stents are geometrically reconstructed during the second step of Stage 1. If necessary, two refining steps are achieved during Stage 2, which is an updating or refining stage. These suplemental steps require a slightly longer calculation time but reach higher accuracy. The first step of Stage 2 consists in recovering the rotation of the stent around its main axis through a minimization loop. The second step of Stage 2 consists in deploying each stent individually.

Figure 1 -Schematic of the general algorithm

A. Data acquisition

In this section, we list the input data and describe how they are processed to extract relevant information and feed the simulations. Data available at the beginning of SG reconstruction include the SG model, updated 3D aorta geometry and intraoperative imaging. The SG models are obtained using the method described in [START_REF] Perrin | Deployment of stent grafts in curved aneurysmal arteries: toward a predictive numerical tool[END_REF], [START_REF] Demanget | Computational comparison of the bending behavior of aortic stent-grafts[END_REF]. Briefly, stent geometries are obtained from manufacturer specifications and are discretized into finite elements using a dedicated Matlab® routine. All simulations are constrained and guided by intraoperative imaging. To isolate the contour of the stents, a combination of Frangi filters [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF] and masks is applied to the image [START_REF] Demirci | 3D stent recovery from one x-ray projection[END_REF]. The Frangi filter is generally used to detect vessels or tubular structures in volumetric image data. The Frangi filter is available in open-source libraries and software such as ITK or ImageJ. The filter includes a measuring scale that allows the isolation of tubular structures of different sizes. By modifying the scale of the filter and combining it with masks, it is possible to extract binary tubular structures of the stents (Figure 2). Then, the convex hull of each stent is extracted. The twodimensional coordinates of stent barycenters are simply obtained from the convex hull (Figure 2). Apparent deployment diameters are also measured. For each stent, the proximal diameter dp and distal diameter dd are recorded. They will be used for the further geometric reconstruction of stents. In the following steps, the SG is virtually positioned in the 3D geometry of the aorta. The aortic geometry, including the centerline and the volume, is previously extracted from the preoperative scan. The volume of the aorta is segmented from the DICOM file of the preoperative scan with a front collision method implemented in VMTK [START_REF]The Vascular Modeling Toolkit[END_REF]. Then, the centerline is extracted using the Voronoi diagram method, also implemented in VMTK. The geometry of the aorta obtained from the preoperative scans may be slightly different of the aortic geometry at the day of the intervention. Indeed, it can be deformed, especially when stiff guidewires are inserted. The aortic geometry must be updated before simulating SG deployment. To do so, the geometry is rigidly and then non-rigidly registered on the intraoperative images. Several methods are available for this step [START_REF] Groher | Deformable 2D-3D Registration of Vascular Structures in a One View Scenario Supplementary Material Derivative of the Difference Measure[END_REF] [42] ( [START_REF] Pionteck | Finite-Element based Image Registration for Endovascular Aortic Aneurysm Repair[END_REF] under review). In addition, we assume that we know the projection matrix of the C-arm.

B. Corotational Euler-Bernoulli beam elements

In this section we describe the corotational Euler-Bernoulli beam elements that are used in the following steps to discretize the simplified geometry of the SG in the global positioning step, and then the stents in the individual deployment step. Simulations are carried out with Project Chrono [START_REF] Project | Chrono: An Open Source Framework for the Physics-Based Simulation of Dynamic Systems[END_REF], [START_REF] Tasora | Chrono: An open source multi-physics dynamics engine[END_REF]. The details of the theory and implementation of beam elements are described in [START_REF] Tasora | Euler-Bernoulli corotational beams in Chrono::Engine[END_REF]. We review here the main concepts. Among the different methods that allow simulating large deformations by finite elements, the corotational approach is one of the most versatile as it is based on classical linear finite elements. The corotational approach allows large displacements, but requires that the strains remain small (Figure 3). A floating coordinate system F follows the deformed element, so that the overall movement in the deformed CD state can be assumed to be composed of a large rigid body movement from the reference configuration C0 to the so-called floating or phantom configuration CS, times a small local deformation from CS to CD. The underlined symbols represent variables expressed in the floating reference basis F. A global tangent stiffness Ke and a global force vector fe are derived for each element e, given its local matrix K, its local force f and the rigid body motion of F in C0 to F in CS. At each time step, the position and rotation of F are updated.

Figure 3 -Schematic representation of the corotational approach

C. Stage 1 (preliminary stage)

Stage 1 combines the first two steps of the algorithm: positioning of the stent barycenters in the aorta and axisymmetric reconstruction of the stents. This step can be run in real-time as it has a marginal computational cost. However, it is based on assumptions that may not be fully satisfied in practice. Thus, the output of this preliminary stage will serve as the starting-point for the updating stage presented in the next section. The first step is to recover the global position of the SG inside the aorta. The SG is simulated with a simplified finite element model. The algorithm for positioning barycenters in 3D is summarized in Figure 4.

Figure 4 -Overview of the algorithm for 3D positioning of barycenters

Barycenter positioning

It is very challenging to find the position of the SG directly in the global reference frame from a single image. Indeed, the intrinsic nature of the C-arm conical projection and the discretization in pixels of the detectors lead to a significant incertitude along the projection axis. A pixel from the flat panel can be assimilated to a surface, therefore its back projection geometry is not a line but a pyramidal volume. It is from this volume that the uncertainty on point positioning comes from (Figure 5). This uncertainty represents the distance along which an object can be moved along the projection axis without moving in the projection plane. This uncertainty depends on the position of the object, i.e. the source-object distance but also the distance from the projection axis. When the pixel gets closer to the projection axis, the uncertainty may tend towards infinity. Close to the edges of the image, the uncertainty becomes lower. For a standard case (source-object distance = 800 mm, source-detector distance = 1300 mm, 750 pixel x 750 pixel detector), a point located 175 pixels from the projection axis has an uncertainty of 4.6 mm.

However, our objective is to reduce this uncertainty as much as possible. Accordingly, our method uses the aortic geometry as a support for the overall SG positioning. The first step is to relate each stent barycenter with the closest point of the artery centerline. This association is achieved in two dimensions. The centerline of the aorta is projected according to the same projection parameters as the intraoperative image. Each two-dimensional barycenter is simply associated with the nearest projected centerline point. This associativity is converted into three dimensions, assuming that the nearest 2D point is almost equivalent to the nearest 3D point. Each barycenter is then related to the corresponding centerline position CL (xCL, yCL, zCL).

The next step is to obtain the back-projection lines of each 2D barycenter. For each stent, we know (from previous image processing) the convex hull and the 2D position of its barycenter BIm (UIm, VIm) in the screen frame. The position of the target image in the 3D space, representing the configuration of the X-ray source and the flat panel detector, is known. Thus, each barycenter on the image can be associated with a three-dimensional position B2D (x2D, y2D, z2D). Since we know the projection matrix, we can obtain the corresponding back projection line for each barycenter. These lines are combined with the geometry of the aorta to obtain the global position of the SG in the 3D frame. Equations of the projection lines result from the elementary Cartesian geometry. The X-ray source is at the origin of the global coordinate system. Therefore, all projection lines have point O (0,0,0) in common. The projection parameters are known. Thus, each B2D target point of the image is associated with a back-projection line passing through this point and through the origin (Figure 5). The normalized vector of line 𝐯 ⃗ (𝐯 𝒙 , 𝐯 𝒚 , 𝐯 𝒛 ) is:

𝐯 ⃗ = 𝑩 𝟐𝑫 𝑶 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ ‖𝑩 𝟐𝑫 𝑶 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ ‖
(1)

Figure 5 -Back projection lines of 2D barycenter coordinates, according to the configuration of the mobile C-arm and uncertainty of positions along the projection axis at an estimated source-object distance

The SG must now be pre-positioned inside the aorta. The coordinates of each three-dimensional barycenter B3D (x3D, y3D, z3D) have then a back-projection line and have previously been associated with the closest centerline point CL (xCL, yCL, zCL). We know that each barycenter 3D coordinate is located on its back projection line, but the coordinate p (Equation 2) of B3D along its line is initially unknown.

𝐁 𝟑𝐃 (𝐱 𝟑𝐃 , 𝐲 𝟑𝐃 , 𝐳 𝟑𝐃 ) = 𝐁 𝟐𝐃 (𝐱 𝟐𝐃 , 𝐲 𝟐𝐃 , 𝐳 𝟐𝐃 ) + 𝑝 * 𝐯 ⃗ (𝐯 𝐱 , 𝐯 𝐲 , 𝐯 𝐳 )

(2)

As a first approximation, we assign to each barycenter the coordinate p from the nearest centerline point such as z3D=zCL. From equation 2, we obtain the following system of equations:

{ 𝐱 𝟑𝐃 = 𝐱 𝟐𝐃 + 𝑝 * 𝐯 𝐱 𝐲 𝟑𝐃 = 𝐲 𝟐𝐃 + 𝑝 * 𝐯 𝐲 𝐳 𝟑𝐃 = 𝐳 𝟐𝐃 + 𝑝 * 𝐯 𝐳 (3)
Hence:

𝑝 = 𝐳 𝟑𝐃 -𝐳 𝟐𝐃 𝐯 𝐳 (4)
(x3D, y3D) are calculated by solving the system of equations [START_REF] Greenhalgh | Endovascular repair of aortic aneurysm in patients physically ineligible for open repair[END_REF]. Then, we have a first approximation of the position of barycenters, based on the centerline. In some cases, this approach is insufficient and must be completed using a finite element model of the SG. Figure 6 shows what can happen in the case of a large aneurysm sac. In this case, the centerline follows the shape of the artery. With the approach described above, stents would be positioned in the sac, which is unlikely and mechanically unrealistic. Therefore, stents that are likely to be badly positioned must be separated from the other ones. In order to define the most precise boundary conditions for the finite element simulation, stents are divided into two categories: free and locked stents. The maximum diameter dSmax of the stent, i.e. the diameter of the fully deployed stent, is compared to the local diameter of the aorta dA at the associated point of the centerline. If dSmax > dA, the stent is locked. In this case, we assume that the stent is in equilibrium in the artery and that its barycenter is therefore very close to the local center of the artery, and therefore to CL. If dSmax < dA, for example if the stent is in the aneurysm sac, the stent is free. The position of the locked stents is assigned according to the centerline (see previous section). The position of the barycenters of free stents will be calculated using a simplified finite element model of the SG. From the initial three-dimensional geometry of the stent, the 3D position of each stent barycentre is extracted. Barycenters are connected between each other by corotational Euler-Bernoulli beam elements according to the initial configuration of the SG [START_REF] Tasora | Euler-Bernoulli corotational beams in Chrono::Engine[END_REF] (Figure 7). The model is set up using the Project Chrono libraries. The SG model in its initial configuration is pre-positioned in the aorta. Then, displacements are prescribed onto the locked stents and the resulting displacements of the free stents are calculated. Free stents cannot go outside the aortic lumen. A SG is a tubular structure with a high degree of mechanical inhomogeneity due to the combination of metal stents and textile graft. The SG model is very simplified, reducing the SG model to a succession of beam elements with the same mechanical behavior. The mechanical characteristics of these beams therefore have no physical reality, and have been optimized to ensure the robustness and stability of the model. As the model is subject to successive boundary conditions (back-projection lines, aortic volume), we assume that this simplified model is sufficient for our application, while allowing a very short computation time. The 3D positions of barycenters are finally determined, hence stent orientation.

Geometric stent reconstruction

From the updated 3D position of the barycenters, a geometric reconstruction of the stent is performed. The initial geometry of each stent, i.e. the metal structure, is first discretized into a set of points (Figure 8). Each point is defined as a vector 𝐕 ⃗ ⃗ (𝐕 𝐱 , 𝐕 𝐲 , 𝐕 𝐳 ), which originates from the stent barycenter B3D and is expressed in the stent local coordinate system. Initially, the local reference coordinate system is the translated global coordinate system. Therefore, the coordinates S(x,y,z) of the n points of a stent are defined by :

𝐒 𝑖:1→𝑛 = 𝐁 𝟑𝐃 𝐢 + 𝐕 𝐢 ⃗⃗⃗ (5)
The proximal deployment diameters dp and distal dd are measured during the image processing step. Here we assume that the stent deployment is axisymmetric. Thus, the diameter measured in the plane of the image is assumed to be the same in all directions. The local deployment diameter d of the stent is therefore interpolated along its main axis z', initially coinciding with the axis z of the global reference frame (Figure 8). New reconstruction vectors 𝐕 𝐝 ⃗⃗⃗⃗ (𝐕 𝐝 𝐱 , 𝐕 𝐝 𝐲 , 𝐕 𝐝 𝐳 ) are updated according to the diameter reduction rd such as:

𝑟𝑑 = 1 - 𝑑𝑆𝑚𝑎𝑥-𝑑 𝑑𝑆𝑚𝑎𝑥 (6)
{ 𝐕 𝐝 𝐱 = 𝐕 𝐱 * 𝑟𝑑 𝐕 𝐝 𝐲 = 𝐕 𝐲 * 𝑟𝑑 𝐕 𝐝 𝐳 = 𝐕 𝐳 [START_REF] Carroccio | Predicting iliac limb occlusions after bifurcated aortic stent grafting: Anatomic and device-related causes[END_REF] Reconstruction vectors 𝐕 𝐝 ⃗⃗⃗⃗ in the global reference frame is expressed in the global frame with the rotation matrix R according to:

( 𝐕 𝐝 ⃗⃗⃗⃗ 1 ) = 𝑹 ( V 𝐝 ⃗⃗⃗⃗⃗ 1 ) (8)
Finally, the new position of the stent is calculated with equation 5 from the updated vectors 𝐕 𝐝 ⃗⃗⃗⃗ . The SG is eventually reconstructed (Figure 8). Positioning and reconstruction of stents is based on two assumptions: the center of gravity of the stents is close to the centerline of the aorta and the deployment of the stents is axisymmetric. Results of Stage 1 (III.C) show that these assumptions are a source of uncertainty during stent reconstruction, which may prevent in some cases to correctly simulate stent deployment. Additional steps of individual stent modeling and deployment corrections are therefore required.

D. Stage 2 (refining stage)

Individual stent deployment may not be accurate enough and need to be improved. This stage combines two individual refining steps: minimization of rotation and individual deployment of stents. The goal of the first step is to determine the actual angle of rotation Φ of each stent, around its main axis z'. The aim of the second step is to individually simulate the deployment of the stent in a deployment box. Both have clinical applications. When the SG is deployed, the actual SG deployment can be reconstructed in 3D. When the SG is not fully deployed, the EVAR practitioner can visualize how the SG would deploy at its current position.

Minimization of the rotational difference

During deployment, the stent may be subject to rotations Φ around its main axis z'. In the case of axisymmetric stents, this rotation has little impact on its final deployment, although the surgeon may wish to improve the accuracy of reconstruction for critical stents. However, in the case of stents with fenestrations, their positioning depends on the rotation of the stent. It is therefore essential to determine these rotations. This step is performed within a minimization loop with a differential evolution algorithm.

In the case of axisymmetric stents, the value to be minimized is the difference eS. This difference is calculated in two dimensions. The 3D model of the stent geometrically reconstructed at the end of the previous step is projected according to the projection parameters of the target image. eS is the average distance between each point of the reconstructed stent and its nearest neighbor. This loop is used to determine the proper rotation Φ ± kθ, where θ is the periodic angle separating two peaks of the Z-shape axisymmetric stent and k a real integer. Considering fenestrated SG, all stents with fenestrations or scallops have radiopaque markers to guide the positioning. The new eM deviation to be minimized is calculated by considering only the distance between the radiopaque markers. In this case, the proper rotation Φ is exact and does not depend on θ. Indeed, the positioning of fenestrations is asymmetrical.

Individual stent deployment

The objective of this deployment box is to make maximum use of intraoperative image data. Thus, the deployment of the stent will be constrained not only by the geometry of the aorta, but also by the information from its convex hull on the image. From the convex hull obtained previously, we define a back-projection polyhedron (BpP). Each side of the BpP is a triangular element. The X-ray source is at the origin of the global coordinate system. Therefore, all triangular elements have a common point O (0, 0, 0). Each edge of the convex hull is defined as the edge opposite to the apex O (Figure 9). The volume is closed by the surface of the convex hull discretized in triangular elements. The Boolean intersection of the BpP with the aorta gives a volume called the deployment box. Boundaries of the volume are meshed with rigid shell elements. The stent is positioned in this rigid box according to its previously determined configuration.

Figure 9 -Deployment box extracted from the intersection of the aorta volume and the BpP

The stents are composed of corotational Euler-Bernoulli beam elements. The mesh size is refined in high curvature areas. The stents are modelled in their 3D position and orientation determined in the previous steps. Only the diameter of the stent is changed. The stent model is initialized in its deployed configuration. Then it is pre-constrained to the diameter of the SG launcher before the simulation begins. Euler-Bernoulli beam elements have an elastic linear behaviour. The stents have a diameter of 0.125mm and are made of 316L steel which mechanical characteristics are summarized in [START_REF] Demanget | Computational comparison of the bending behavior of aortic stent-grafts[END_REF].

Each stent, or a selection of stents at the practitioner's discretion, is deployed individually. The crimped stent is positioned inside the deployment box. The position of its barycenter and its orientation have already been recovered during the global positioning step. Then the stent is deployed (elastic recoil as the stent was crimped). The deployment is calculated using the Project Chrono engine [START_REF] Tasora | Chrono: An open source multi-physics dynamics engine[END_REF], with solver Math Kernel Library (MKL) from Intel ® . After the first contact, the time step is reduced to ensure the stability of the model. The contacts are modelled using the penalty algorithm implemented in Project Chrono, the Smooth-Contact (SMC) modeling approach. SMC uses penalty (in a discrete element method (DEM) [START_REF] Tasora | A compliant viscoplastic particle contact model based on differential variational inequalities[END_REF], [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF], regularizing the frictional contact forces, with "imaginary" spring-dashpot systems at each contact) and as such objects in contact will have slight interpenetration and integration time-step will likely be small. The simulations are performed on a computer with 4 CPUs, 3.40GHz, 16 GB RAM, but without parallelization. The computation time for deploying a stent is less than 6 minutes, without optimization. In addition, the complete calculation is easily parallelizable, as each stent deployment can be simulated on a separate core. The overall calculation time could therefore be compatible with clinical use.

Figure 10 -Flowchart of the validation scheme

Several parameters were used to assess the quality of the reconstruction. The first parameter was the DB distance, which was the distance between the 3D barycentres of the target stent and the reconstructed stent. This distance was used to assess the quality of stent positioning within the artery. The second parameter was the distance DM, which was the average distance between the point clouds of the target stent and the reconstructed stent. It was defined as the average of the Euclidean distances between a node of the reconstructed stent and its nearest neighbor among the points of the target stent. distance allowed reaching the quality of stent positioning and deployment to be assessed at the same time. Finally, the last parameter is CSAS, i.e. the cross-sectional area overlap of target and reconstructed stents, which allows the quantitative comparison of stent deformations at the SGs folds. Cross-section areas along the entire stent were measured after deployment in the aneurysm. After deployment, the cross-sectional areas were modified, particularly in terms of SG folds. Let AT be the area of the target cross-section ST and AR the area of the reconstructed cross-section SR. AU is the area of the SU intersection between ST and SR.

𝑆 𝑈 = 𝑆 𝑇 ∪ 𝑆 𝑅 (9)
𝐶𝑆𝐴𝑆 = 100 * 𝐴 𝑇 -|𝐴 𝑇 -𝐴 𝑈 | 𝐴 𝑇 (%) (10)

C. Results and discussion

Figure 11 shows the results of Stage 1 and Stage 2 and Table 2 shows a comparison of the results of the two stages. Concerning Stage 1, the average DB is generally lower than the clinical validation value, the positioning of stents in the artery is generally good. However, the error is too large on the contralateral limb, close to or greater than 3 mm with a maximum of 3.71 mm. About DM, the average error is less than 3 mm, however the distance map shows that a significant number of stents have an error too close to the limit. The CSAS map confirms this result, with an average superposition of about 70% and a minimum of 47%, which is much too low. Thus, an additional step of individual deployment seems necessary. Stage 2 showed a clear improvement in the quality of stent simulation. All measurements were below the threshold value, with a maximum for the DB of 1.40 mm and for the DM of 2.28 mm. The mean values were also improved, by reducing the mean error of -26.5% for the DB and -21.7% for the DM compared to Stage 1. There was also a clear improvement in the CSAS by a mean value of 25%. Figure 12 shows how fenestration positions (front view) were predicted. The size difference between the target and reconstructed fenestrations results from the segmentation of the post-operative SG. The fenestrations are surrounded by radioopaque markers that create artifacts during imaging. These artifacts make the markers appear larger than reality. The average distances between the centers of gravity of the fenestrations are summarized in Table 3. The overlapping of fenestrations is improved after the refining steps, which is visible on the Figure 12. The method, after correction, was therefore able to position the fenestrations precisely enough to allow the catheterization of secondary arteries, with a maximum distance less than 2.46 mm for the left renal fenestration.

-Summary and comparison of Stage 1 and Stage performances

These results are interesting in comparison with previous work. Indeed, preoperative finite element simulations [START_REF] Perrin | Patient-specific numerical simulation of stent-graft deployment: Validation on three clinical cases[END_REF] have an average error of about 3.5mm, which is higher than the results obtained here, with a longer computation time, as they are intended for planning purposes. It is though difficult to compare the previous fast methods with the one presented here. Indeed, the evaluation criteria are not similar. For example, in [START_REF] Spranger | Comparison and calibration of a real-time virtual stenting algorithm using Finite Element Analysis and Genetic Algorithms[END_REF], the results of the fast method (FM) were compared with finite element simulations (FE), which may themselves differ from reality. In complex cases, the average distance between FM and FE is about 3-4 mm, for a very short calculation time (<1min). However, these methods are not guided by intraoperative imaging. In [START_REF] Zhou | Real-Time 3-D Shape Instantiation From Single Fluoroscopy Projection for Fenestrated Stent Graft Deployment[END_REF], an accuracy of about 2-3 mm is reached, which is comparable to our results, with a shorter computation time. However, their method was based on additional markers bonded onto the device, which does not seem suitable. Compared to these previous studies, our method seems to achieve a promising compromise between accuracy, computation time and compatibility with current clinical conditions. The results of Stage 1 (first and second steps) show that the SG position is globally well predicted in the artery, but that the radial deployment of stents is not sufficiently accurate, with errors exceeding the 3mm threshold. However, corrections made in Stage 2 (third and fourth step) yield a reasonable accuracy. More specifically, accurate predictions of fenestration positions after Stage 2 would avoid complex catheterization of secondary arteries, which is one of the major source of complications for practitioners.

However, the method has several limitations. First, it depends on the quality of the input data. Indeed, we assume that stents can be individualized using image processing. If this assumption is usually satisfied, overlapped stent can challenge it, for example, in the case of very pronounced angulation of the proximal neck of the aneurysm. In addition, the geometry of the aorta is assumed to be known, as it was updated using a non-rigid registration method previously. However, registration errors could occur and add to the other errors mentioned in the chapter. The non-rigid registration step also involves restriction for clinical applications. As we assumethat we know the geometry of the aorta before the SG deployment simulation, the nonrigid registration must be performed each time the mobile C-arm position is changed, which implies to perform a DSA in order to visualize the aorta. The method should therefore only be used at key points in the procedure. Moreover, segmentation errors of the aorta volume are possible, but they are difficult to identify and probably have a minor influence on the results. A priori, the presence of thrombus is not a major problem for segmentation. Indeed, the aortic lumen is segmented during preoperative planning, not the aorta itself, and the SG is deployed in the lumen. Effects of thrombus on overall aortic stiffness are considered in the previous non-rigid registration step. But local variations in stiffness, due to thrombus, calcifications or surrounding tissues are not taken into account. Other limitations are related to the method itself. In order to save calculation time, the deployment of stents was simulated individually within rigid boxes. The rigid nature of the boxes is obviously a simplified mechanical behavior of the arterial wall but it did not induce significant errors when the simulations were compared with the post-operative CT scan. Finally, by simulating the individual deployment of stents, we do not consider the effect of the textile graft connecting them together. For example, a highly crimped stent can change the deployment diameter of the neighboring stent independently of the local geometry of the aorta. However, initial results show that these limitations do not hamper the accuracy more than what is compatible with clinical expectations. The influence of other assumptions should be investigated further with additional patient data. This method would be particularly suited for complex cases of thoracoabdominal aortic aneurysms that require the use of fenestrated stent grafts. Indeed, when positioning fenestrations, the surgeon can potentially encounter difficulties. In most cases, the method is able to provide assistance to the surgeon. A few exceptions can challenge the method. The method is based on data extracted in the plane of the intraoperative image. Missing information along the projection axis are provided by simulations. If the main deployment axis of the stent graft was located along the projection axis, the method may have difficulty simulating the device. In practice, this situation seems very unlikely though.

IV. CONCLUSION

We have presented in this chapter a methodology for simulating SG deployment based on intraoperative images coupled with mechanical models. The algorithm consists of a series of successive steps with increasing precision. A first step simulates the positioning of the SG within the aorta using a simplified finite element model and the centerline of the artery. The second step is an axisymmetric geometric reconstruction of the stents. The third step minimizes the rotation of the stent around its main axis. Finally, the last step consists in deploying each stent individually within a deployment box extracted from the imaging and geometry of the aorta.

The results of combined Stage 1 and Stage 2 yield a reasonable accuracy. More specifically, accurate positioning of fenestration would facilitate catheterization of secondary arteries. But the method suffers from several limitations. First, it depends on the quality of the input data and the ability of image processing to distinguish stents. In addition, the geometry of the aorta is supposed to be known. Next, the wall of the artery is considered rigid when the stents are deployed. Lastly, textiles are not simulated. The influence of these limitations on the accuracy of simulations needs to be explored further using additional patient data.

Finally, the translation of our methodology seems promising and could be generalized to all operating rooms equipped with mobile C-arms to assist SG deployment using real-time simulations in the future.

Figure 2 -

 2 Figure 2 -Example of stent deployment in a 3D printed AAA replica: original image (a), stent detection (b) and extraction (c), convex hull (d) and 2D barycenters (e)

Figure 6 -

 6 Figure 6 -Error in barycenter positioning due to a deformed centerline: normal centerline (a), centerline deformed in an aneurysm sac (b), resulting unrealistic stent positioning (blue) in the sac area and realistic positioning (green) in non-deformed sections (c), expected actual positioning (d)

Figure 7 -

 7 Figure 7 -From the simplified geometry of the SG (a), 3D coordinates of barycenters are first extracted (b), then barycenters are connected using beam elements (c)

Figure 8 -

 8 Figure 8 -Stent geometric reconstruction, with the local frame O'(x',y',z') and the global frame O(x,y,z), initial model (a) and after reconstruction (b)

Figure 11 -

 11 Figure 11 -Distance map representing the distance between the barycenters of the target and the reconstructed SG from Stage 1 (a) and Stage 2 (b), the RMS distance between the cloud points of target and the reconstructed stents from Stage 1 (c) and Stage 2 (d), and the cross-sectional surface superposition between target and reconstructed stents from Stage 1 (e) and Stage 2 (f). Each square represents a stent.

Figure 12 -

 12 Figure 12 -Comparison of the positions of the fenestrations. In red, the target, in blue the reconstructed fenestration after Stage 1, in green after Stage 2. Right renal fenestration (a), mesenteric fenestration (b) and left renal fenestration (c)

Table 3 -Distance between the center of gravity of target fenestrations and reconstructed fenestrations, before and after Stage 2

 3 CPUs, 3.40GHz, 16 GB RAM, the calculation time for Stage 1 is under 20 seconds. The calculation time for the rotation minimization is under 1 minute. The maximum calculation time for the deployment simulation is 6 minutes. Co optimization and parallelization are currently ongoing to enable applications in clinical conditions.

		Right renal	Mesenteric	Left renal
		fenestration	fenestration	fenestration
	DB Stage 1 (mm)	3.06	3.05	6.05
	DM Stage 2 (mm)	0.41	1.64	2.46
	With a computer powered by 4

III. PROOF OF CONCEPT

The method described above was applied to a patient who underwent an EVAR procedure.

A. Clinical data

Details about the device are given in Tab. 1. An additional stent was tested, including three fenestrations: right renal artery, mesenteric artery and left renal artery fenestrations. The preoperative and postoperative CT scans were acquired at the Saint-Etienne University Hospital under clinical conditions. Use of the clinical data was approved by the institutional review board and informed consent was obtained from the patient. The voxel size of the scans was 0.9395x0.9395x2 mm 3 . First, we evaluated the results of Stage 1. As the results were not precise enough in terms of radial expansion, a second step of individual stent deployment was required. We evaluated the results of the complete method, including Correction Part. The average diameter of the renal arteries was 5 to 7 mm. If the fenestration positioning error is less than 3 mm, the surgeon can catheterize the secondary arteries such as renal arteries. Above this threshold, it is considered that intraoperative complications are likely to arise. Therefore, the clinical validation value was set at 3mm in accordance with experienced clinicians.

B. Quality assessment of stent deployment

First, to isolate and test the stent reconstruction algorithm as precisely as possible, the following assumptions were made: the projection matrix was known; the 3D geometry of the aorta was assumed to be perfectly registered. This ensured that not introducing positioning errors related to aortic registration. Thus, target images and the 3D model of the aorta were generated from the postoperative scan. The actual 3D position of the SG was therefore known, and served as a reference to be compared with the simulation for the sake of validation (Figure 10).


Kevin Marchais
email: kevin.marchais@ensam.eu

Jérémie Girardot

Charlotte Metton

Ivan Iordanoff

• C Metton

• I Iordanoff

B K Marchais

A 3D DEM simulation to study the influence of material and process parameters on spreading of metallic powder in additive manufacturing

Keywords: Laser beam melting, Additive manufacturing simulation, Discrete element method, Granular behavior, Powder spreading

A 3D DEM simulation to study the influence of material and process parameters on spreading of metallic powder in additive manufacturing

Introduction

In metallic powder-bed-based additive manufacturing processes, for each layer of matter, 2 main phases occur. First, there is the powder spreading to form a powder bed, then there is the laser scanning to melt selected zones that correspond to the sections of the fabricated part. In order to produce pieces of better quality, faster and at the lowest cost, it is important to study the influence of the different process parameters to determine the best configurations of parameters on the deposited powder layer. Studies about powder deposition in additive manufacturing processes started recently, this work takes part in this domain of study. According to the grains disposition in the powder bed, laser-matter interactions can be impacted and can affect the quality of the final fabricated part. Indeed, if the powder bed is not dense enough or heterogeneous, some issues can happen such as porosity, lack of matter, balling effect, etc. Zielinski et al. [START_REF] Zielinski | Influence of powder bed characteristics on material quality in additive manufacturing[END_REF] studied the influence of the powder grains layout on laser melting with comparison between spreaded powder bed and bodycentered cubic (BCC) structured powder bed and showed that BCC powder bed produces a homogeneous melting whereas spreaded powder bed shows irregular melt pool. Then it is important to understand which parameters influence the final powder bed quality. Production time is impacted by the time of the powder spreading phase because a large number of layers is needed to produce a piece (e.g., 5000 layers of 60 µm needed for a height of 30 cm). One can want to increase spreader speed or increase the thickness of deposited layers to speed up the spreading phase but how does it affect final powder bed quality? Some machines use a blade as coating system while others use a roller or a brush. Some spreaders are rigid (metals), others are flexible (silicon). Haeri [START_REF] Haeri | Optimisation of blade type spreaders for powder bed preparation in additive manufacturing using dem simulations[END_REF] studied the optimization of the spreader geometry to enhance powder bed quality and obtain better powder bed than with a roller or a classic blade.

Fouda et al. [START_REF] Fouda | A dem study of powder spreading in additive layer manufacturing[END_REF] studied the effect of spreader speed and gap thickness on packing fraction and particles velocities with monodisperse (every particle has the same size) powder samples. They showed that the packing fraction could be enhanced by varying spreader speed and gap thickness. Chen et al studied the effect of spreading powder with counter-rolling roller as spreader and showed that increasing spreading speed decreases powder bed packing density and surface quality [START_REF] Chen | Packing quality of powder layer during counter-rolling-type powder spreading process in additive manufacturing[END_REF].

Depending on the material, the grain size or powder distribution, the final powder bed quality is different. Indeed, the material parameters such as adhesion combined with smallest particles highly deteriorates powder bed quality due to the formation of clusters and cavities as shown by Chen et al. [START_REF] Chen | Flow behavior of powder particles in layering process of selective laser melting : numerical modeling and experimental verification based on discrete element method[END_REF][START_REF] Chen | Powder-spreading mechanisms in powder-bed-based additive manufacturing: experiments and computational modeling[END_REF] and Meier et al. [START_REF] Meier | Critical influences of particle size and adhesion on the powder layer uniformity in metal additive manufacturing[END_REF] on TiAl6V4 powders. In this case, adhesion forces become very important compared to the grains weight.

Then the study of the adhesion parameter is fundamental. One of the cause of adhesion in powder is humidity. Higher is the humidity, higher is the adhesion. In LBM process, humidity level is extremely low so there are only Van der Waals' forces that play the role of adhesive force. However, it could be interesting to understand if less strict environment and powder storage conditions could lead to good results in powder bed quality in order to reduce the production costs. Then it is important to understand how does adhesion affect powder spreading in order to know in what conditions should be stored the powder and if powder spreading should be done under controlled atmosphere or not for example. Powders material parameters are complicated to indentify for these models and the best solution seems to be the measurement of the repose angle of powder samples as Meier et al. did [START_REF] Meier | Modeling and characterization of cohesion in fine metal powders with a focus on additive manufacturing process simulations[END_REF].

Several studies have been realized about Selective Laser Sintering (SLS) process which is very close to Laser Beam Melting (LBM) as a laser is used to transform a spreaded powder bed. The difference is about the material used in the process as LBM uses metallic powders whereas SLS uses polymeric powders. For example, Parteli and Pöschel [START_REF] Parteli | Particle-base simulation of powder application in additive manufacturing[END_REF] studied the spreading of PA12 powders (polymer) and have taken into account the non-sphericity of the particles by simulating complex shape with multisphere method, and Haeri et al. [START_REF] Haeri | Discrete element simulation and experimental study of powder spreading process in additive manufacturing[END_REF] studied the spreading of rod-shaped polymeric particles with rollers as spreaders.

In metal powder-based processes, the particles are highly spherical thanks to the gas/water atomization process [START_REF] Hoeges | Additive manufacturing using water atomized steel powders[END_REF]; however, Nan and Ghadiri [START_REF] Pasha | Jamming during particle spreading in additive manufacturing[END_REF] applied the multisphere method to 316L samples which presents less spherical grains than inconel or titanium-based powders.

The grains size is of the order of magnitude of the surface roughness properties of the plate or substrate. Then the surface roughness has to be taken into account. Nan and Ghadiri [START_REF] Pasha | Jamming during particle spreading in additive manufacturing[END_REF] chose to simulate surface roughness with cylinders arranged side by side.

So until now, only a few studies focus on the powder bed spreading simulation. The cited studies use the discrete element method to simulate powder behavior. This study, also based on discrete element method, aims to understand the behavior of polydisperse metallic powders with highly spherical particles and show the influence of surface roughness and spreader speed on powder bed properties. The influence of material parameters such as friction and adhesion is also studied in order to have a better knowledge on how precisely these parameters must be determined.

First, the numerical model developed to realize this study will be presented. Then, the results of preliminary studies including surface roughness influence will be presented. Finally, the parametric study and the obtained results will be discussed.

Numerical model

Numerical method

The Discrete Element Method (DEM), initially presented by Cundall [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF], is the chosen method to realize the simulations of this study as it is particularly adapted to granular media simulations. Each powder grain is represented by a discrete element, and the fundamental principle of dynamics is applied to each of them in order to get their acceleration. Then, by integration (Verlet velocity scheme), the velocity and the displacement of each element are computed. The forces that are considered in the simulations are gravity and interactions forces between the bodies. Simulations are computed in 3D with the open-source platform GranOO developed in C++ [START_REF] André | The granoo workbench, a new tool for developing discrete element simulations and its application to tribological problems[END_REF] and based on the Discrete Element Method. This method considers small interpenetrations between the particles to detect a collision between two particles. Unlike polymeric powder [START_REF] Parteli | Particle-base simulation of powder application in additive manufacturing[END_REF], metal grain powders are very spherical but can present small satellites at their surface. In this study, the particles are supposed perfectly spherical and satellites are not taken into account.

Choice of interaction forces

When a contact occurs between 2 bodies, a repulsion force is generated. There is also normal energy dissipation that is taken into account with a damping coefficient. There can also be a force that tends to maintain contact between bodies if they are adhesive. Finally, a friction force is applied if the relative velocity between the two bodies is not collinear to the contact direction. These 4 forces are detailed below:

Repulsive force

Modeling of the repulsion between two grains with a contact stiffness (Eq. ( 1)). This is a normal force.

-→ F r = kδ -→ n (1)
where k is the contact stiffness (N/m),

δ = R 1 + R 2 -d 12
is the interpenetration between the grains, with R 1 and R 2 the radii of the particles and d 12 the distance between their centres, δ < 0 when there is a contact.

Other models are used to model the repulsion between the grains such as the Hertz contact model [START_REF] Hertz | Ueber die berührung fester elastischer körper[END_REF]. It considers the slight deformation of two spheres that are in contact and take into account Young modulus and Poisson coefficient of materials (Eq. ( 2)).

-→ F H = 4 3 E * √ R * δ 3 2 -→ n (2)
However, the high value of Young modulus makes the contact stiffness very high. The time-step is depending on the contact stiffness and the grain mass according to the condition (3) (spring-dashpot model).

t < 0.2 M min k (3)
The additive manufacturing powder is constituted by grains whose diameter is the size of a few microns. For particles of this size, the mass is very low. It results with a very low time-step. The high computation cost of the Discrete Element Method do not allow to use Hertz model in this case because of high values of Young modulus. To obtain reasonable time computation, a higher time-step is needed. Studies using Hertz model reduce Young modulus as it doesn't affect the powder behavior as shown by Chen et al. [START_REF] Chen | Effect of young's modulus on dem results regarding transverse mixing of particles within a rotating drum[END_REF]. Meier et al. [START_REF] Meier | Modeling and characterization of cohesion in fine metal powders with a focus on additive manufacturing process simulations[END_REF] also showed that modifying the stiffness value doesn't affect the results of powder pile simulation. The contact stiffness is then modified because it does not impact the global granular behavior with the tested values. Nevertheless, some conditions must be respected : the interpenetration must be limited to avoid numerical issues. Then the contact stiffness must be high enough to face adhesion force and particles weight.

Energy dissipation

When there is a contact, a part of the energy is dissipated. The chosen model to represent it is described by Eq. ( 4) :

-→ F d = 2α √ k M V rel -→ n ( 4
)
where α is the dissipation coefficient, M the reduced mass of the two grains in contact, V rel the relative velocity of the bodies in contact. This model allows to link easily the dissipation coefficient to the restitution coefficient often used in granular problems with the relation α =

1 1+ π ln( ) 2

Adhesion

When 2 bodies are in contact, there is a surface energy that creates an attractive force between them. The chosen model to represent this phenomenon is described by Eq. ( 5). This one is very close to the JKR model [START_REF] Johnson | Surface energy and the contact of elastic solids[END_REF] and DMT model [START_REF] Derjaguin | Effect of contact deformations on the adhesion of particles[END_REF] that accounts for Van der Waals forces and is also close to models that allows to represent adhesion with liquid bridges when the humidity level is high enough in the powder [START_REF] Gladkyy | Comparison of different capillary bridge models for application in the discrete element method[END_REF]. JKR, DMT and liquid bridges models have the same shape of equations where a coefficient is multiplied by R * . This coefficient is constant with JKR and DMT models. So it has been chosen to take γ as the whole constant and to determine its value.

-→ F a = γ R * -→ n (when δ < 0) (5)
with γ the surface energy (J /m 2 ) and

R * = R 1 R 2 R 1 +R 2 (m)
the equivalent radius corresponding to the two grains.

Friction

When 2 bodies are in contact and in movement, a force is opposing to the tangential relative displacement for each body. This is the friction that is simulated by Coulomb law (Eq. ( 6)). This model could also allow to represent the satellites fixed around the powder grains [START_REF] Tan | An overview of powder granulometry on feedstock and part performance in the selective laser melting process[END_REF], increasing artificially the macroscopic friction between 2 grains.

-→ F t = -μ V T rel || -→ F r || -→ t (6)
The friction coefficient μ depends on the relative velocity between the 2 bodies. If the tangential component of the relative velocity is positive, then the friction coefficient is equal to + μ, else if it is negative then it is equal to -μ (Fig. 1). When the relative velocity tends to zero, the function must be regularized.

-→ t is the tangential vector to the contact. The principle of this method is represented in Figs. 2 and3 by the contact between 2 particles with given velocities. The forces generated by the contact are represented by the normal and tangential components of the total resulting force. Here, there is a regularization of the repulsive force by accepting an interpenetration δ and of the friction force by the calculation of the μ factor. This is a notable characteristic of this numerical method. On the opposite, the Non-Smooth Contact Dynamics [START_REF] Moreau | Non-smooth mechanics and applications[END_REF] manages non regularized forces.

Integration of acceleration

In order to compute particles motion, acceleration can be computed with the sum of the forces acting on the particle (interaction forces and weight) and its mass. Then velocity and position can be computed from acceleration using Verlet velocity integration scheme (Eq. ( 7)).

⎧ ⎨ ⎩ a n+1 = F/M v n+1 = 1 2 t(a n+1 + a n ) p n+1 = v n+1 t + a n+1 t 2 (7)
with a the acceleration, v the velocity and p the position, M the mass of the particle and F the sum of the forces applied on the particle. Angular acceleration and velocity are also computed to get the particles rotation with the use of quaternions which allows quick computations of rotations. Angular acceleration is given by Eq. ( 8) [23,24].

α n+1 = ω n • θn • ω n + 1 2 θ n I -1 τ -4( θn • ω n ) ∧ I .( θn • ω n ) ( 8
)
with α the angular acceleration, ω the angular velocity, θ the rotation, I the moment of inertia of the particle and τ the sum of the torques applied on the particle. Then angular velocity and rotation can be computed with Verlet velocity integration scheme (Eq. ( 9)).

ω n+1 = 1 2 t(α n+1 + α n ) θ n+1 = ω n+1 t + α n+1 t 2 (9)

Domain of study

To realize this study, all the simulations are computed with the same base geometry shown in Fig. 4. It contains a first plate where powder is put before the spreading (blue grains) and the fabrication plate on which powder is spread. The second plate is a little lower than the first in order to correspond to the wanted layer thickness. In this study, powder spreading is studied exclusively on a rigid and rough surface. There is no spreading on a previous non-melted powder layer. The spreader is a supposed infinitely rigid blade with a speed V .

In the fabrication zone, a layer is already existing (red) in order to simulate the surface roughness of the melted previous layer. The red layer contains discrete elements that cannot move. It allows to simulate plate/substrate surface roughness. A melted surface has a roughness of the order of magnitude of the grains size, so the irregular surface created by the discrete elements allows to approach the real roughness of the substrate. The plates width is 300 µm, and periodic boundary conditions are applied in the transverse direction in order to simulate infinite width. The plates length is 4 mm. The studied domain is relatively small because the numerical method is very expensive in calculation time, so the number of discrete elements must be limited. Nevertheless, this length is long enough to observe the powder global behavior during the spreading. The material taken into account for this study is a metallic powder of IN718 which powder grains are very spherical. Simulations are computed with spherical grains with a size distribution that follows a normal law with a diameter of 30 µm and a standard deviation of 10 µm that is similar to measured IN718 powder samples. For example, Table 1 compares a measure of virgin IN718 powder sample deciles reported by Nguyen et al. [START_REF] Nguyen | Characteristics of inconel powders for powder-bed additive manufacturing[END_REF] with numerical particle size distribution.

Each layer has approximately 10000 discrete elements. The parameters values used in the simulations are reported in Table 2. The contact stiffness is 0.1 N/m because it allows to have small enough interpenetrations between the grains and an acceptable calculation time (around 10 h). The overlap between 2 particles must be very small in order to keep the simulation stable. The condition fixed here is the following δ R < 10 2 with δ the overlap and R the radius of the particle. The contact stiffness must be high enough to balance with weight or adhesion force. Restitution coefficient is set to 0.1 to stabilize the simulation. Meier et al. [START_REF] Meier | Modeling and characterization of cohesion in fine metal powders with a focus on additive manufacturing process simulations[END_REF] showed that contact stiffness and restitution coefficient can be modified without affecting powder behavior. The other parameters are variable.

Figure 5 shows the powder behavior during the spreading. The spreader collects the powder on the first plate, and a pile is forming ahead of it as shown with experimental results by Chen et al. [START_REF] Chen | Powder-spreading mechanisms in powder-bed-based additive manufacturing: experiments and computational modeling[END_REF]. When the spreader reaches the fabrication zone, the powder spills under the spreader to form the new layer.

Powder bed properties

In order to compare the final powder bed for different parameters configurations, it is important to define properties that will permit to analyze its quality. The chosen measured properties are the packing fraction, the coordination number and the average grain size. Pieces produced by LBM process can present a high level of porosity if the powder bed is not dense enough. The packing fraction must then be as high as possible to minimize the porosity and allow the best melting efficiency. The study of the average grain size can be interesting according to the position in the powder bed to check segregation phenomenon that has already been observed on experimental testings. This phenomenon has already been observed in other domains such as mixers for example [START_REF] Poiron | Dynamics of size segregation and mixing of granular materials in a 3d-blender by nmr imaging investigation[END_REF]. In spreading process, the powder flow could result in the fact that small particles would rather be at the beginning of the spreading zone whereas bigger particles would be at the end. Finally, the coordination number is the average number of contacts per grain. This parameter is strongly linked to powder bed thermal conductivity since there is better thermal conduction trough solid matter than through gas. It is then important to get a high and homogeneous coordination number. A post processing zone is set and identified by the blue zone in Fig. 4 in the powder bed. This zone is subdivided (Fig. 6) to compute properties in each cell and compare it through the position in the powder bed. The packing fraction is computed by the ratio between the volume of each grain belonging to a cell and the volume of the cell.

Surface roughness

It has been chosen to take into account plate (or substrate) surface roughness because this is supposed to be of significant impact on powder spreading and final powder bed properties. To simulate the surface roughness, Nan et al. [START_REF] Pasha | Jamming during particle spreading in additive manufacturing[END_REF] chose to use clumped cylinders. In this study, the chosen method consists in disposing discrete elements regularly as shown in Fig. 7. First, it allows to control the roughness parameters by adjusting the discrete elements diameters. Then, contrary to Nan et al, the roughness is taken into account in every directions. However, these methods create a regular roughness whereas real part surface roughness is more random.

Simulations have been computed with different discrete elements diameters to highlight the surface roughness influence on powder flown and powder bed properties. Figure 8 shows the impact of surface roughness on powder spreading as the angle of the powder pile ahead of the spreader is increasing as the surface roughness increases. Table 3 lists powder bed density for each case, and it shows that the powder bed is more dense with a higher surface roughness. So surface roughness affects the powder flow and the final powder bed quality. For the following studies, the surface roughness is simulated with particles diameters of 30 µm.

Results

The following studies are part of a parametric study where the goal is to observe the influence of a material or process parameter by varying it. For all the studies, all parameters are fixed except the parameter that is being investigated. The friction coefficient is set to µ = 0.3. The adhesion energy is set to γ = 0 mJ/m 2 in order to not make appear adhesion effect when studying the effect of other parameters. The spreader speed is set to 100 mm/s. The chosen particle size distribution follows a normal law with an average diameter of 30 µm, a standard deviation of 10 µm, a minimal diameter of 7 µm and a maximal diameter of 70 µm. Especially for material parameters, the aim of these studies is to understand the influence of the parameter and determine how precisely it should be identified with further experiments.

Friction coefficient effect

The influence of friction coefficient on powder bed spreading is necessary to understand to evaluate the precision needed in this coefficient characterization. Even if for a given material, the friction coefficient is a proper material parameter, this study permits to understand the powder behavior whose friction coefficient is more or less high. As Fig. 9 shows, the powder bed seems to have globally the same properties for the studied parameters values (μ = 0, μ = 0.01, μ = 0.1, μ = 0.3, μ = 0.9, large range of metal friction coefficient) and no tendency is emerging. According to these results, the friction coefficient has not a consequent impact on the spreading process so it will be assumed to be constant for the following.

Friction has no influence on granular behavior during powder spreading so it is not necessary to identify precisely its value. Then the effect of adhesion energy, the second material parameter, will be studied.

Adhesion energy effect

As for the friction, the adhesion depends also on the powder sample material studied but also on the environment such as humidity level. It is necessary to understand the behavior difference between a non-adhesive powder and a more or less adhesive powder. Simulations have been computed with different adhesion energies (γ = 10 mJ/m 2 , γ = 5 mJ/m 2 , γ = 1 mJ/m 2 , γ = 0.5 mJ/m 2 , γ = 0.1 mJ/m 2 , γ = 0 mJ/m 2 , the maximum adhesion value is chosen after numerical flow tests that showed that powder flow was very deteriorated with γ = 10 mJ/m 2 ). Figure 10 shows that adhesion has an important impact on the final powder bed. Density and coordination number evolve in opposite directions. For a very adhesive powder, there will be particles heaps where the density will be locally very high (high coordination number) but also cavities that deteriorate the average powder bed density as shown in Figs. 10 and11.

On the contrary, the powder bed obtained with nonadhesive powder is rather dense but with a low coordination number compared to the other testings. The simulations realized with γ = 1 mJ/m 2 and γ = 0.5 mJ/m 2 provide a density similar to the ones provided with γ = 0 mJ/m 2 with a much higher coordination number.

Unlike friction coefficient, adhesion energy has a significant effect on powder bed final properties. This parameter will have to be precisely identified in order to compute powder bed spreading simulations. The following studies treat on process parameter influence beginning with spreader speed influence on powder bed properties.

Influence of the spreader speed

In order to observe the spreader speed influence on the powder bed final properties, several simulations have been computed with different speeds ((50, 75, 100, 125, 150, 175, 200, 225 and 250 mm/s). Figure 12 shows the evolution of the density, of the average grain size and of the coordination number according to position for the following velocities: 50, 100, 150, 200 and 250 mm/s. As expected, slower the spreading is, better the powder bed density is and higher the coordination number is. Indeed, the obtained density for a speed of 50 mm/s exceed 60% whereas for the higher speed (250 mm/s), the powder bed is less dense (40-50%). The values measured at the first and last positions are not the Fig. 8 Surface roughness influence most representative data since it corresponds to the lack of matter caused by the grains inertia for the first case and by the fall of the grains beyond the plate for the second. Moreover, at high velocities, the powder bed is not homogeneous. The density curve profile is irregular compared to the slowest velocities. The coordination number depends also on the spreading speed since a 50 mm/s spreading permits to obtain the best average coordination number (around 6). Finally, the average grains size slightly increase according to the position in the powder bed for the highest speeds. It could be the segregation phenomenon as smallest particles are mainly at the beginning of the plate. This study also shows that density evolves like coordination number in this case unlike with adhesion study. Coordination number can be considered as local density.

Discussions

The parametric study of material and process parameters permitted to show that some parameters have a significant impact on powder bed properties whereas some parameters seem not to have a major influence. First, surface roughness has a significant impact on powder flow and final powder bed properties. It is then an important parameter to consider when simulating powder spreading. Concerning the friction coefficient, it seems not to have an influence on the powder bed properties. On the other hand, adhesion parameter is very impacting the powder bed properties. Low adhesion results in a homogeneous and dense powder bed. Whereas high adhesion produces an heterogeneous powder bed with clusters and cavities. Even if material parameters cannot be directly modified for a given powder sample, humidity affects directly adhesion energy so dry environment are preferable for the powder. It seems to have a threshold effect as adhesion energy is 1 mJ/m 2 reasonable but it is preferable not to be more adhesive than that. According to these results, only adhesion must be identified precisely for a given powder sample and must be controlled in real process. Spreader speed is an important process parameter and its effect on powder bed quality was expected. A lower speed produces a better powder bed. Here, it also seems to have a threshold effect as for speeds above 100 mm/s , spreading generates heterogeneous powder bed. Finally, in the presented studies, the segregation phenomenon has not been clearly identified except a particular case when spreader speed is very high.

Conclusion

In order to produce the most representative simulations of real spreading, the material parameters of powder samples will have to be identified precisely. Since the adhesion energy is the only parameter to identify, the experimental testing of powder angle of repose can be realized in comparison with equivalent simulations as shown by Meier et al. [START_REF] Meier | Modeling and characterization of cohesion in fine metal powders with a focus on additive manufacturing process simulations[END_REF]. Except on the adhesion effect study, adhesion energy has been taken to 0 mJ/m 2 in order to eliminate the adhesion effect when studying effect of other parameters. However, it must be noted that powders always present at least a minimal amount of adhesion energy. More complete parametric studies will be computed in the future to investigate the effect of other process parameters and also the eventual coupling between several parameters. Finally, simulations with longer plates could be computed to try to identify segregation effect.

Fig. 1 Fig. 2

 12 Fig. 1 Regularization of the Coulomb's law

Fig. 3

 3 Fig. 3 Angular motion

Fig. 4 Fig. 5 Fig. 6

 456 Fig. 4 Domain of study

Fig. 7

 7 Fig. 7 Rough plate (top view)

Fig. 9 Fig. 10 Fig. 11

 91011 Fig. 9 Friction coefficient

Table 1

 1 IN718 powder sample deciles

	Sample	D 10 µm	D 50 µm	D 90 µm
	Real	21	31	49
	Numerical	20	31	43

Table 2

 2 Simulation parameters

	Parameter		Value
	ρ	Density (kg/m 3 )	8000
	d	Diameter of a particle (µm)	7-70
	L	Plate length (mm)	4
	V	Recoater speed (mm/s)	50-250
	k	Contact stiffness (N/m)	0.1
	e	Restitution coefficient	0.1
	γ	Adhesion energy (mJ/m 2 )	0-10
	μ	Friction coefficient	0-0.9

Table 3

 3

	Roughness effect on global powder bed density	Roughness	Density (%)
		Smooth	49.64
		7 µm	4 6 .79
		30 µm	5 3 .97
		50 µm	6 3 .09

Acknowledgements This work is a part of Additive Factory Hub (AFH) platform. Simulations have been computed with GranOO workbench : www.granoo.org.