Nathan Lauzeral
email: nathan.lauzeral@ec-nantes.fr

Domenico Borzacchiello

Michael Kugler

Daniel George

Yves Rémond

Alexandre Hostettler

Francisco Chinesta

A model order reduction approach to create patient-specific mechanical models of human liver in computational medicine applications

Keywords:

Background and objective: This paper focuses on computer simulation aspects of Digital Twin models in the medical framework. In particular, it addresses the need of fast and accurate simulators for the mechanical response at tissue and organ scale and the capability of integrating patient-specific anatomy from medical images to pinpoint the individual variations from standard anatomical models.

Methods:

We propose an automated procedure to create mechanical models of the human liver with patient-specific geometry and real time capabilities. The method hinges on the use of Statistical Shape Analysis to extract the relevant anatomical features from a database of medical images and Model Order Reduction to compute an explicit parametric solution for the mechanical response as a function of such features. The Sparse Subspace Learning, coupled with a Finite Element solver, was chosen to create lowrank solutions using a non-intrusive sparse sampling of the feature space.

Results: In the application presented in the paper, the statistical shape model was trained on a database of 385 three dimensional liver shapes, extracted from medical images, in order to create a parametrized representation of the liver anatomy. This parametrization and an additional parameter describing the breathing motion in linear elasticity were then used as input in the reduced order model. Results show a consistent agreement with the high fidelity Finite Element models built from liver images that were excluded from the training dataset. However, we evidence in the discussion the difficulty of having compact shape parametrizations arising from the extreme variability of the shapes found in the dataset and we propose potential strategies to tackle this issue.

Conclusions: A method to represent patient-specific real-time liver deformations during breathing is proposed in linear elasticity. Since the proposed method does not require any adaptation to the direct Finite Element solver used in the training phase, the procedure can be easily extended to more complex nonlinear constitutive behaviors -such as hyperelasticity -and more general load cases. Therefore it can be integrated with little intrusiveness to generic simulation software including more sophisticated and realistic models.

Introduction

For the sake of clarity, abbreviations and main mathematical notations are summarized in Tables 1 and2 . virtual information on real images of the patient to see through opaque tissues. In IR, this augmented representation of the patient is important to visually follow the targeted organs. For example, in radiotherapy such tool allows to follow the tumor and predict its displacement during the breathing [START_REF] Hostettler | A real-time predictive simulation of abdominal viscera positions during quiet free breathing[END_REF] , permitting to minimize the radiated area and consequently the trauma for the patient. Concerning the MIS, the AR is used as a back up for two major drawbacks of the technique: the reduced visibility and in the case of robot-assisted MIS the loss of haptic feedback during the operation [START_REF] Okamura | Haptic feedback in robot-assisted minimally invasive surgery[END_REF] . Thanks to AR, it is now possible to provide the surgeon with additional information such as the position of tumors [START_REF] Haouchine | Image-guided simulation of heterogeneous tissue deformation for augmented reality during hepatic surgery[END_REF][START_REF] Haouchine | Towards an accurate tracking of liver tumors for augmented reality in robotic assisted surgery[END_REF] or blood vessels [START_REF] Ferrari | Augmented reality visualization of deformable tubular structures for surgical simulation[END_REF][START_REF] Haouchine | Image-guided simulation of heterogeneous tissue deformation for augmented reality during hepatic surgery[END_REF] . Modifications in the organ topology due to cuts can even be taken into account [START_REF] Paulus | Surgical augmented reality with topological changes[END_REF] .

Results on haptic feedback technologies in AR, although not yet robust enough to be used in the ORs, are showing promising results and are currently used as training and learning tools [START_REF] Panait | The role of haptic feedback in laparoscopic simulation training[END_REF][START_REF] Van Der Meijden | The value of haptic feedback in conventional and robot-assisted minimal invasive surgery and virtual reality training: a current review[END_REF] . However, these approaches are often more challenging when applied to soft tissues because of the difficulty to continuously adapt the AR scene to the position of the features of interest. In abdominal surgery, for example, a lot of factors can interact with the organs. They can come from the internal or external environment (surgical tools, insufflation of gas in the peritoneum) or be physiological movements (heart beats, breathing). Hence the necessity to develop robust real-time biomechanical models of the tissues to integrate these movements and reproduce them with fidelity.

Scope of the current work

In most computational medicine applications the goal is to create a digital replica of a considered biophysical system that can realistically reproduce the most essential observed features. Two key aspects in this are the possibility to customize the models for specific use cases (accounting for inter-patient and inter-population variabilities) and the real-time interactive response to new assimilated data. This emerging technology, known as the Digital Twin, merges complex biophysical modeling and advanced real-time simulation techniques with data assimilation and analysis for decision support. Some model and simulation aspects were addressed in the review by [START_REF] Cueto | Real time simulation for computational surgery: a review[END_REF] in the surgical framework and the use of machine learning is investigated in [START_REF] Lorente | A framework for modelling the biomechanical behaviour of the human liver during breathing in real time using machine learning[END_REF][START_REF] Martínez-Martínez | Estimation of the elastic parameters of human liver biomechanical models by means of medical images and evolutionary computation[END_REF][START_REF] Martínez-Martínez | A finite element-based machine learning approach for modeling the mechanical behavior of the breast tissues under compression in real-time[END_REF] for parameters identification and mechanical behavior prediction. More specifically, models for visual and haptic feedack applications were developed in [START_REF] Niroomandi | Real-time deformable models of non-linear tissues by model reduction techniques[END_REF][START_REF] Niroomandi | Accounting for large deformations in real-time simulations of soft tissues based on reduced-order models[END_REF][START_REF] Niroomandi | Real-time simulation of biological soft tissues: a pgd approach[END_REF] , where the non-linearities and load parametrization issues are tackled, and in [START_REF] Niroomandi | Real-time simulation of surgery by reduced-order modeling and x-fem techniques[END_REF][START_REF] Quesada | Computational vademecums for real-time simulation of surgical cutting in haptic environments[END_REF] , where the simulation of surgical cutting is dealt with. In this paper we build on these works and focus on the aspects related to the integration of medical images in real-time interactive models to personalize the organ anatomy. We propose a numerical framework to personalize biomechanical models interactively using new anatomical data without rebuilding the models from scratch.

The specific anatomy of the patient is generally taken into account case-by-case thanks to pre-operative data such magnetic resonance images (MRI) or computed tomography (CT) scans. From these inputs, the surface and the volume of the organs of interest can be identified, reconstructed and meshed [START_REF] Bartz | Visualization and exploration of segmented anatomic structures[END_REF][START_REF] Heimann | Comparison and evaluation of methods for liver segmentation from ct datasets[END_REF][START_REF] Kugler | A new algorithm for volume mesh refinement on merging geometries: application to liver and vascularisation[END_REF][START_REF] Zhao | Direct multiphase mesh generation from 3d images using anisotropic mesh adaptation and a redistancing equation[END_REF] . This can be a time-consuming and computationally intensive task and needs to be repeated for any new patient, posing severe limitations to the use of such models within interactive simulation environments. Then, to reach real-time performance the models either need to simplify the formulation, e.g. using lumped parameters models, corotational finite elements or beam structures [START_REF] Haouchine | Image-guided simulation of heterogeneous tissue deformation for augmented reality during hepatic surgery[END_REF] , or use specific hardware set-up for fast parallel computations [START_REF] Joldes | Real-time nonlinear finite element computations on gpu-application to neurosurgical simulation[END_REF] . In the first case, over simplification in the model formulation may lead to incorrect mechanical behavior unless the lumped parameters are accurately tuned to capture the essential features of the equivalent distributed parameters model. In the second one, it involves high performance computing resources that are rarely available in a clinical environment. In both cases, computational time will strongly depend on the mesh refinement, limiting the accuracy to the hardware capacity. In this paper we follow the route of collocation-based Model Order Reduction (MOR). This approach allows to simultaneously tackle the specific anatomical representation and the real-time constraint without being limited by the simplification of the physics or the computational cost and bypassing mesh generation and model assembly costs. The reduction of computational costs also enables the use of inexpensive and simple hardware such as tablets or even smartphones, which could be easily installed in the ORs.

Reduced order modeling for organ twins models

MOR methods have seen a growing interest this last decade. They allow to reduce the computational complexity in numerical simulations by a parametrization of the solution, enabling real-time online computations whithout simplifying the underlying physics of the model. In return, a computationally intensive offline stage must be done beforehand. They are good candidates to the creation of Digital Twins because they handle data assimilation through their parametric formulation. A distinction can be done between the projection-based and the collocation-based MOR methods. The underlying idea in the former is to seek the solution for the model's governing equations in a lower-dimensional subspace that is specific for the problem at hand. This is accomplished by imposing the orthogonality of the equations residuals with respect to the new subspace. The resulting system of equations has a considerably reduced computational complexity, while the accuracy of the solutions is preserved. Based on the way the reduced basis is learned , projection-based MOR methods are usually divided in two categories: the a posteriori methods like the Proper Orthogonal Decomposition [START_REF] Chatterjee | An introduction to the proper orthogonal decomposition[END_REF] or Reduced Basis method [START_REF] Quarteroni | Reduced basis methods for partial differential equations: an introduction[END_REF] , where the reduced order model (ROM) is built from a set of training solutions, available from previous simulations (called snapshots), and the a priori methods like the Proper Generalized Decomposition (PGD) [START_REF] Chinesta | The Proper Generalized Decomposition for Advanced Numerical Simulations: A Primer[END_REF] , in which the parametric solution is assumed to respect the canonical tensor format and is then built from the equations governing the problem. These methods have been used for different applications in relation with the medical framework, and in particular computational surgery [START_REF] Cueto | Real time simulation for computational surgery: a review[END_REF] . Projection-based ROM including geometrical parameters have been addressed in [START_REF] Manzoni | Shape optimization for viscous flows by reduced basis methods and free-form deformation[END_REF] for optimization purposes and in [START_REF] González | Computational patient avatars for surgery planning[END_REF] to synthesize new shapes from manifold interpolation. When using real medical data, the first difficulty to create anatomical models encompassing the shape comes from its parametrization. This latter can be obtained through the use of Statistical Shape Models (SSMs), and, in particular, Point Distribution Models (PDMs) [START_REF] Cootes | Statistical models of appearance for computer vision[END_REF] . From a training set of real anatomical shapes the PDMs associate to a template the principal modes of shape variation. This way, the geometry of a specified organ can be statistically described by a limited set of parameters allowing image segmentation [START_REF] Heimann | Statistical shape models for 3d medical image segmentation: a review[END_REF] or the creation of Finite Element (FE) shape parametrized models [START_REF] Bredbenner | Development and validation of a statistical shape modeling-based finite element model of the cervical spine under low-level multiple direction loading conditions[END_REF][START_REF] Lauzeral | Shape parametrization of bio-mechanical finite element models based on medical images[END_REF][START_REF] Nicolella | Development of a parametric finite element model of the proximal femur using statistical shape and density modelling[END_REF] .

A key assumption for projection-based ROM is that the problem variational form is affine with respect to the parameters. This guarantees the complexity reduction in the projection step and is therefore fundamental for the performance of the method in terms of execution time. Unfortunately, affine approximations are not always easily recovered. This bottleneck is most likely to be encountered in problems involving a parametrization of the domain geometry, as in our case. To efficiently handle the representation of the shape, we opted for a collocation-based MOR method, the Sparse Subspace Learning (SSL) [START_REF] Borzacchiello | Non-intrusive sparse subspace learning for parametrized problems[END_REF] . On the contrary to projectionbased approaches, in collocation-based methods the residual is not constrained to be orthogonal to the solution's subspace but rather enforced to be exactly zero at properly chosen points in the parametric space, called the collocation points. The parametric solution is built by interpolating the solutions of the model at these points, which makes the ROM an a posteriori model. In our case, the main advantage is that the approach is fully data-driven, as it does not require the evaluation of the residual nor its projection onto the lower dimensional space. Hence, the lack of affinity of the variational form does not affect the performance of the method. As the shape parametrization is based on a PDM which results in the loss of affinity, it motivates our choice of the SSL. This latter is based on the use of hierarchical collocation to build the parametric model, exploiting sparsity and low-rank representation for the solutions. SSL aims for the sparsity of the representation using a sparse grid sampling approach while simultaneously building a reduced order representation of the solution. New out-of-training solutions can be computed by a hierarchical interpolation of the snapshots coefficients in the reduced basis. Since SSL is only based on the output of a FE solver and does not require any modification of it, it is regarded as a non-intrusive technique. However, the use of hierarchical collocation in the parametric space, in our case hierarchical polynomials, requires a fundamental condition on the regularity of the solution in the parametric space of shape features. To preserve this regularity, we adopted a morphing technique to smoothly adapt the FE mesh of the model according to the shape parameters [START_REF] Lauzeral | Shape parametrization of bio-mechanical finite element models based on medical images[END_REF] .

Previous work

In our previous paper [START_REF] Lauzeral | Shape parametrization of bio-mechanical finite element models based on medical images[END_REF] we developed a method to build shape-parametrized FE models from PDMs through a morphing technique. The main idea is to create the SSMs in such way that the generated shapes are adapted to the iFEMWRAP method [START_REF] Shontz | A robust solution procedure for hyperelastic solids with large boundary deformation[END_REF] used to reconstruct the volumetric mesh. This latter enables the reconstruction of a 3D mesh from the knowledge of an initial state plus displacements applied to the boundary nodes, i.e. the external mesh's surface. By iteratively minimizing the energy of deformation of the volume the position of the inner nodes is computed in such way that the mesh is adapted to FE computations. The main constraint is that the imposed displacements on the surface do not generate a tangled surface mesh. To this end we introduced a registration algorithm called Thin Plate Spline Parametrized Registration (TPS-PR) (see Appendix A ), which maximizes the morphing efficiency. By parametrizing the shape of a template 3D volumetric mesh through a Statistical Shape Analysis (SSA) using the TPS-PR algorithm, the iFEMWRAP can then be applied to any generated shape using the template as initial state.

In the remainder of this paper this method is going to be enhanced and adapted to the SSL in order to create the parametric ROM (pROM) for breathing simulations, as an example application. This is motivated by previous works from Hostettler et al. [START_REF] Hostettler | Bulk modulus and volume variation measurement of the liver and the kidneys in vivo using abdominal kinetics during free breathing[END_REF][START_REF] Hostettler | A real-time predictive simulation of abdominal viscera positions during quiet free breathing[END_REF] , where a method to predict the abdominal organs and tumors positions during free breathing was developed.

Overview and paper organization

This paper presents a method based on the SSL to create patient-specific mechanical models of the human liver by taking into account their anatomy. To do so, a SSA is used in order to get the shape parametrization of the human liver, allowing to recreate the shape of any new patient. By combining the SSA and the SSL, one can solve in real-time a chosen mechanical problem on any shape taken into account by the parametrization. In this paper the specific case of the breathing is developed. During the breathing, the displacement of the liver's surface can be estimated. From this input the presented method allows to reconstruct the organ inner displacements accordingly to the material's constitutive behavior used to represent it, enabling the visualization of tumors or blood vessels positions within the liver. We first present the liver shape database used for the SSA in Section 2.1 . The mechanical model is introduced in Section 2.2 and the associated breathing simulation is presented in Section 2.3 . The method to create the pROM is described in Section 2.4 , followed by the process to adapt the model to a new patient in Section 2.5 . Then, error metrics are defined and generated results are presented and commented in Sections 2.6 and 3 . Finally, the identified limits of the method are highlighted in the discussion and some potential solutions are suggested.

Data acquisition, processing and assimilation as well as more technological aspects are not covered in this work and will be object of a parallel publication.

Materials and methods

If not specified otherwise, all vectors mentioned by a bold italic letter belong to a three dimensional space ( x, y, z ).

A point cloud X ∈ R n ×3 will correspond to X = (x 1 , • • • , x n , y 1 , • • • , y n , z 1 , • • • , z n ) T .

Data

A database of T = 385 external surface meshes of livers was provided by the IRCAD (Institut de Recherche contre les Cancers de l'Appareil Digestif, Strasbourg, France). To our knowledge this is one of the biggest and diverse database used in the literature to build SSMs on a human organ. This database contains a large scale of liver shape collected on various kind of patients (male or female, healthy or unhealthy) in supine position. The shape meshes have different number of vertices (from 20 0 0 to 150,0 0 0). Due to the fact that this database has been collected over many years, different qualities of segmentation are included. Some of the shapes are represented in Fig. 1 .

A template shape close to the average target is chosen as reference (see Fig. 2 ) and meshed with gmsh [START_REF] Geuzaine | Gmsh: a 3-d finite element mesh generator with built-in pre-and post-processing facilities[END_REF] . Its surface is described by v S = 1393 vertices and 2782 triangles. The volume is meshed with 10,163 tetrahedrons for a total of v V = 2452 vertices.

Its maximum length is 260 mm . Special attention has been paid to the mesh quality using mesh optimization functions implemented in the open software gmsh. One of the interest of the procedure is to be able to display internal features to see through the opaque tissues of the liver. The position of the portal vein is deduced from the liver displacement fields computed through the mechanical equilibrium laws ( Eq. ( 1) ).

Model assumptions

The pROM of the liver is created under the assumption of the quasi-static mechanical equilibrium with a prescribed displacement on the boundaries. The governing equations read:

∇ • σ = 0 in the volume u = u d on the boundaries ( 1
)
where u is the displacement field and σ the Cauchy stress tensor. The user-prescribed motion u d can be deduced from the position of the abdomen and that of the surrounding organs. For example, for an intergration in a complete AR tool one could use the work presented in [START_REF] Hostettler | A real-time predictive simulation of abdominal viscera positions during quiet free breathing[END_REF] where the position of inner organs is estimated from the abdomen surface motion parametrization under free breathing. The stress-strain relationship is assumed to be linear and the material properties are considered homogeneous and isotropic, although, hypereleastic behavior could also be considered. Fig. 3 shows an example of geometry in which the vascularization position is updated based on the computed displacement field.

Breathing simulation application

The liver is considered as an elastic solid defined by E = 3 kPa and ν = 0 . 48 . Its movement during the breathing is modeled by applying fixed Dirichlet boundary conditions to the surface nodes:

u (b) = b × u d on the boundaries, (2)
where b ∈ [0, 1] is the breathing parameter. The displacement u d has been computed by taking the extreme positions of a real liver during the breathing in vivo, from data provided by the IRCAD. By registering the template shape to the maximum inhalation and exhalation positions using the TPS-PR algorithm, the displacement can be directly applied to the template mesh.

In this case a single parameter has been used to describe the breathing motion. Multiple parameters could be used as well possibly leading to a more accurate representation. The parametrization could then be obtained by applying a SSA to a database of images corresponding to maximum inhalation and exhalation positions of several patients.

Creation of the parametric reduced order model

In the general case, the offline procedure to build pROM can be separated in two main steps: the parametrization of the problem and the use of a separated variables representation to enable fast computations during the online stage. An illustration of the method is given in Fig. 4 .

Problem parametrization. To take into account the variability of the model according to the fields of interest, a parametric approach is used. Each field of interest is described by one or several parameters. In this paper a distinction is done between the explicit and the implicit parameters (respectively noted p = (p 1 , . . . , p D ) and ξ = (ξ 1 , • • • , ξ n p ) ). The firsts encompass all kind of parametrization where the parameters are directly deduced from the formulation of the fields of interest. For example, in the case of the material properties, the Young modulus or the Poisson ratio. The same can goes for boundary conditions. In the breathing application presented in this article p = b is deduced from Eq. ( 2). On the contrary, implicit parameters are "hidden" within the model and must be inferred through methods such as the SSA used to represent the liver's shape. This latter is detailed in Appendix A and Appendix B using the data presented in subsection 2.1 , resulting in the shape parametrization described by Eq. ( 10) and the parameters subspace given by Eq. ( 12) . Fig. 5 a-c illustrate how the shape parameters subspace boundaries are defined.

Separated variables representation.

Once the model is parametrized, the solution -here the displacement field -is expressed using a separated variables representation. More specifically, the canonical tensor format is used, which allows to overcome the exponential complexity of multi-parametric models. This reads:

u (x , ξ, p) = d i =1 α i B i 0 (x ) n p k =1 B i k (ξ k ) D l=1 B i n p + l (p l ) .
(

) 3
This formulation expresses the multidimensional field u as the truncated modal expansion in which each mode is the product of lower dimensional functions B i j expressing the parameter dependence. d is called the canonical rank of u . For the application considered in this paper, Eq. (3) becomes:

u (x , ξ, b) = d i =1 α i B i 0 ( x ) n p k =1 B i k ( ξ k ) B i n p +1 ( b ) .
(4) Fig. 4. Summary of the method to construct the patient-specific pROM. First an offline stage is done where the shape parametrization is built and the SSL is applied. Others parameters p than the shape can be added by the user, in this paper the breathing motion is added using one parameter such as p = b. This step can be computationally intensive but is easily parallelizable. Computational details and notations can be found in Appendix A, Appendix B and Appendix C . Then, for the online stage, the solution is written in a compact canonical tensor format (see Eq. ( 3) ) allowing fast computations. This second step enables the real-time visualization of the solutions through the parametric formulation.

To reach this specific formulation the SSL is used. The method is described in Appendix C . To put it briefly, the wor k mentioned in Section 1.4 allows to get the FE model from the shape parametrization through a morphing technique. The FE problem defined by Eq. ( 2) can then be solved for any shape ˜ X (ξ) (see Eq. ( 10) ). The SSL iteratively computes the solution of this FE problem on wellchosen collocation points, corresponding to specific set of parameters ( ξ, b ). By interpolating these snapshots a prediction of the solution over the whole parameters subspace can then be obtained.

At each iteration the interpolation error is reduced and the algorithm stops when it becomes smaller then a given tolerance. The solution is finally written as in Eq. ( 4) to allow real-time computations. Fig. 5 d illustrates the position of the collocation points for the three first increments of the SSL in the shape parameters subspace.

Model personalization for patient-specific anatomy

Once the parametric model is built, it can be adapted to a new patient liver anatomy thanks to the shape parametrization. To find the shape parameters associated to a specific anatomy, the idea is to embed the SSMs into the registration technique leading to a reduced complexity formulation. Indeed, with SSM-based registration a n p × n p system needs to be solved on each iteration, which is done on only few seconds. Once the shape parameters are estimated, the simulated breathing motion is readily computed as a particularization of Eq. (4) for the new set of shape parameters. Since the parametric solution has been computed offline the amount of online work to obtain this solution is fairly inexpensive and compatible with real-time constraint. Indeed, this particular solution format allows to display at a frequency of 25Hz for visualization up to 1 kHz [START_REF] Niroomandi | Real-time simulation of biological soft tissues: a pgd approach[END_REF] .

It must be noticed that a new shape can lie outside of the boundaries of the parameters subspace. In that case, the nearest shape inside of the parameters subspace boundaries -in terms of Euclidean distance -is considered.

Validation

For each collocation point of the SSL the rigidity matrix corresponding to the parametrized shape is computed as well as the corresponding displacements of the inner nodes. Considering only the shape, the breathing problem can be seen as adding a shape parameter. The main interest is then to be able to follow the motion of the whole volume when the breathing displacement is applied to a given shape. In a nutshell, without the SSL the following steps need to be applied for each new shape:

1. find the shape parameters, 2. reconstruct the interior of the mesh with the iFEMWRAP, 3. compute the rigidity matrix corresponding to the FE model, 4. solve the FE boundary conditions problem for a given value of the breathing parameter.

With the SSL once the shape parameters are found the solution for any breathing parameter is given by Eq. ( 4) . The reconstruction with the iFEMWRAP, the rigidity matrix computation and the solving are bypassed. This means that the SSL error comes from the "simplification" of these steps. Consequently its error lies in the reconstruction of the interior, either from the iFEMWRAP or from the solving of the FE breathing problem. To assess the accuracy of the SSL "leave-one-out" tests are performed. To this end, the ROM is built using T -1 liver shapes. The last shape is then reconstructed with M modes and the error between the ROM reconstruction and the complete model -i.e. the one built following the 4 steps -is evaluated. Several sources of error are identified:

1. The initial non-rigid registration onto the target shapes. This error is hard to quantify as the correspondence between vertices is not known. That is why the non rigidly registered template is considered as the ground truth afterward in order to be able to compare the vertices position in surface and volume. To estimate the registration accuracy the Volumetric

Overlap Error VOE (S, S ) = 100 × (1 -| S ∩ S | / | S ∪ S |
) was evaluated [START_REF] Heimann | Comparison and evaluation of methods for liver segmentation from ct datasets[END_REF] , where S and S represent the compared shapes. An

Table 3

Number of collocation points, i.e. computations, per number of shape parameters n p as a function of the SSL hierarchical level i . One mode being used for the breathing the total number of parameters is n p + 1 . The level i is limited to 3 in order to stay within a reasonable number of computations to do for n p = 10 shape parameters. These values could be increased with more computational ressources. Adding a hierchical level would increase the accuracy of the SSL and adding a shape parameter would improve the shape representation.

error of 100% means the shapes are completely dissociated, an error of 0% means they perfectly overlap. 2. The projection of the shape on a partial set of principal axis obtained with the SSA. Knowing the solution with the whole set of modes it is possible to compute the error as the distance between corresponding vertices as:

E p (M) = 1 v V s ∈ || s -s ||
where represent the whole set of nodes, s the points of the projected shape and s the corresponding points of the ground truth shape. The volume of both shapes is reconstructed with 30 iFEMWRAP iterations.

The interpolation done by the SSL approximation of the model.

This error can be measured by comparing the results from the iFEMWRAP plus the FE problem on one hand and the SSL on the other hand. This error reads:

E SSL (M) = 1 v V s ∈ || s -s ||
where represent the whole set of node, s the points of the SSL rebuilt volumetric mesh and s the corresponding points of the FE computed solution done on the iFEMWRAP reconstructed FE model. A distinction is done between b = 0 , where no FE solution is computed but only the iFEMWRAP, and b = 1 where both are.

The maximum global error is the sum of all this independent errors. The tests are done on 20 randomly chosen shapes from the database. The values of E p and E SSL are then averaged.

In our previous paper [START_REF] Lauzeral | Shape parametrization of bio-mechanical finite element models based on medical images[END_REF] the same liver shape database was used to create the SSMs. The quality of the SSA was then evaluated through three metrics: the compactness, the specificity and the generalization. The non-rigid registration method used was found as efficient as the classical Thin Plate Spline registration method.

Results

To perform the SSA, the TPS-PR was run with n b = 5 and λ r = 0 . 5 (see Appendix A for parameters description). 35 increments were done for each registration. To create the pROM, the SSL was limited to the third hierarchical level to stay within a reasonable number of collocation points. Table 3 resumes the number of collocation points per level in function of the number of shape parameters. With a MATLAB (The MathWorks, Inc., USA) implementation a computation in a collocation point -i.e. running the iFEMWRAP and solving the FE problem -took between 2 and 30s on an Intel Xeon E5-2680v3 2.5Ghz processor. We used the Parallel Toolbox on 24 cores to drastically reduce the computational time since all computations are independent and can be trivially parallelized.

First, the non-rigid registration VOE is presented in an histogram in Fig. 6 . Globally the shapes were well-registered, the error mainly comes from the smoothing of the surfaces by the TPS-PR method which does not particularly deteriorate the shapes representations.

Then, the mean error done by the sole reconstruction is presented in Fig. 7 for M = 1 , • • • , 100 shape modes. An example is provided in Fig. 10 a. The median value is also represented. It is shown that its value is always lower than the mean. This is due to the fact that some of the shapes are very peculiar and consequently require a high number of modes to be represented. That was expected considering that the database contains a large variety of livers, lot of them belonging to ill patients. Because of the large variability of the database both measures of error are quite important. 50 modes are required to get a median error lower than 4mm, this value is given in [START_REF] Hostettler | Bulk modulus and volume variation measurement of the liver and the kidneys in vivo using abdominal kinetics during free breathing[END_REF] as the threshold where protocol in radiotherapy could be significantly improved by such tool. Seeking to reach such accuracy can be debated, though. Indeed, there is an initial error in our source data, resulting from image segmentation Fig. 8. N exc : number of shapes excluded from the parameters subspace by the point selection. N out : number of sets of parameters taken out of the parameters subspace in the leave-one-out tests. These values depend on the considered number of modes M. N exc is approximately constant and represents around 5% of the number of shapes in the database. On the other hand N out steadily increase with M , meaning that the chance of having a new shape out of the pre-computed solution increases with the dimensionality. and non-rigid registration. Therefore, although the additional computational effort would improve the model accuracy, it would not increase its precision.

Next, the suitability of the parameters subspace was evaluated. Fig. 8 summarizes how many shapes were removed from the training set to build the n p -hyper-ellipsoid (see method in Appendix B ) and how many reconstructed shapes among the "leave-one-out" tests were outliers, i.e. they were located out of the parameters subspace boundaries. It appears that the number of shapes excluded from the training set varies around 20 and does not seem to depend directly on the number of modes. On the contrary, the number of outliers increases with the number of modes. Such sets of parameters are exceptional and automatically introduce an error as they do not belong to the pre-computed solution and need to be projected onto the n p -hyper-ellipsoid. For a standard use they should not be considered, that is why a distinction is made afterward between them and the inliers, i.e. the shapes located inside the parameters subspace boundaries.

Eventually, the error between the FE solution and the SSL prediction is shown in Fig. 9 a and b. As mentioned before a distinction was made between inliers and outliers. Moreover, the error was computed for b = 0 and b = 1 to assess the impact of the breathing parametrization on the error. The breathing approximation has a low impact on the error, around 0.1 μm . On the contrary, taking shapes out of the parameters subspace introduces a bigger error.

With inliers only the maximum error is around 1 μm against 1mm when all shapes are taken into account. In both cases the error increases with the number of modes, showing that the convergence is more difficult for higher dimensions. A visual representation of this error is provided in Fig. 10 b for b = 1 .

Discussion

Globally, the error done in the leave-one-out tests comes from the reconstruction from a limited number of modes. This error could be decreased by using more modes to represent the shape. However, more modes means more computations for the SSL and a more difficult convergence. Here it is difficult to go further than 10 shape modes for the SSL without reaching high computational times. Concerning the SSL error, the main flaw lies in the choice of the parameters subspace. As shown in the previous section the SSL error mainly comes from reconstructed shapes taken out of the parameters subspace, and their number increases with the number of modes. Nonetheless, this error stays within acceptable range and for shapes taken inside the parameters subspace a really high fidelity reconstruction is done. Moreover, a consequent speedup is achieved thanks to MOR while taking into account the same physical equations than the initial FE problem, without making additional assumptions. As noticed in Section 3 this initial model is already computed relatively fast, and with an implementation in a low-level programming language and a more powerful hardware a computation could certainly be done in near real-time ( ∼ 1s), which would be sufficient for breathing applications and would question the use of a pROM. As mentioned before the use case presented in this article is a simple example to illustrate the method and more complex situations involving non-linear behaviors should be tested. Anyhow, the pROM approach still has the advantage of requiring simple computational means, which is important to compensate for the lack of advanced hardware in clinical environments.

The first point to tackle to improve the method would be to reduce the number of modes necessary to represent the shape. As the model is based on a large liver shape database coming from medical data, there is an important shape variability. An idea to reduce this variability would be to clusterize the data beforehand. By sorting similar shapes into clusters the variability would be reduced inside each of them [START_REF] Gooya | Mixture of probabilistic principal component analyzers for shapes from point sets[END_REF] . A specific ROM would be then created for each cluster. When a new shape would be introduced, an additional step consisting in finding in which cluster it lies would be done. The next difficulty to deal with would be the parameters subspace choice. It needs to be a compromise between completeness and efficiency of the SSL while keeping a simple and parametrizable topology. On the contrary to the subspace taken in our previous paper [START_REF] Lauzeral | Shape parametrization of bio-mechanical finite element models based on medical images[END_REF] , we did not chose a κ%-boundary model [START_REF] Lu | Statistical shape analysis of clavicular cortical bone with applications to the development of mean and boundary shape models, Comput[END_REF] . Such subspace considers a normal distribution of the data along each dimension and fix the boundaries separately depending on their standard deviation. Consequently, if one wants to cover 95% of the variability along each dimension, it does not mean that 95% of the whole variation will be covered. In fact the variation covered lies in ([0 . 95 n p , 0 . 95] × 100)% where n p is the number of modes -i.e. the dimension. Hence, for an increasing number of modes there is a possible decrease of the variation coverage. That is why the method is not used here and that we use a non-statistical determination of the subspace. The right way to cover κ% of the variation would be to have a multivariate normal distribution of the data. Here, for a number of modes superior to 2 the data do not follow such distribution according to the Henze-Zirkler test [START_REF] Henze | A class of invariant consistent tests for multivariate normality[END_REF] . If the data could be clustered with a multivariate gaussian mixture model as proposed before, such choice for the parameters subspace boundaries would be straightforward.

Conclusion

A data-based pROM has been developed with the purpose to be integrated into a complete liver Digital Twin. This ROM allows to take into account the patient-specific shape and to model the mechanical deformations caused by the breathing motion. The procedure relies on a SSA in conjunction with the SSL. We used a previously developed registration method, the TPS-PR, to increase its efficiency. The principal novelty of this approach, when compared to literature, is the use of medical data for the shape parametrization of the ROM. A first application to simulate the breathing was done on a simple elastic case, but other applications are envisaged. We showed that, because of the large size of our database, the main limitation is the ability to statistically represent the initial liver shape with a small number of modes. Some leads are suggested to solve this issue, in particular the data clusterization. The whole process works smoothly and provides an efficient way to address the initial matters of patient-specific shape representation within interactive simulation environment, contributing to open the way to the creation of MOR-based tools for patient-specific real-time simulations.

A prototype to visually follow the liver movements during the breathing of any patient is currently being developed at the IRCAD using this approach. This will be the opportunity to validate in vivo the pertinence of this method.

Appendix A. Statistical shape analysis

The SSA is done using our previous work [START_REF] Lauzeral | Shape parametrization of bio-mechanical finite element models based on medical images[END_REF] . The aim of the SSA is to create the shape parametrization from a training data set of shapes. Two main steps are identified: the registration process and the dimension reduction with the Principal Component Analysis (PCA) [START_REF] Jolliffe | Principal Component Analysis[END_REF] . Here, the registration step is unsupervised, meaning there is no a priori knowledge on the correspondence of anatomical landmarks between each image. To do this, a well-meshed template shape is first selected. Then, all shapes from the databasethe targets -are rigidly registered to this reference shape. Next, the template undergo a non-rigid registration, called Thin Plate Spline Parametrized Registration (TPS-PR), to fit each one of the targets. The TPS-PR is done by minimizing the sum of two terms, a point maching term I match between the surfaces and the bending energy I tps of a set of control points describing the template through a spline interpolation: I(a ) = I match (a ) + λI t ps (a ) , [START_REF] Bredbenner | Development and validation of a statistical shape modeling-based finite element model of the cervical spine under low-level multiple direction loading conditions[END_REF] where λ is a weighting scalar. The point matching term reads:

I match (a ) = 1 2
|| A a -d|| 2 [START_REF] Bungartz | Sparse grids[END_REF] and the bending one:

I t ps (a ) = 1 2 a T L a ( 7
)
where a is the vector associated to the control points of the spline representation such as u = A a , where u are the displacements of the templates vertices and A the interpolant operator. The matrix L is the biharmonic operator of the TPS bending energy. The splines are described through a FE approach using Hermitian elements of the second order. The discretization n b characterizes the number of control points along each spatial dimension and is user-defined. These points are placed in a regular grid around the template such that they define a bounding box, which global size is fixed to +5% the one of the template. A gradient descent algorithm is used in conjunction with a golden-section search to minimize Eq. ( 5) . As it can be hard to tune the parameter λ, λ r is introduced as:

λ = λ r • max α∈ [0 , 1] I match (a ) max α∈ [0 , 1] I t ps (a ) λ r ∈ R + ( 8
)
where α is the golden-section search step.

Once the template is non-rigidly registered onto each target, a database of deformations associated to the template bounding box is defined. The PCA is applied to this database in order to extract the n p principal modes of deformations φ i . The representation of a in the reduced basis reads as:

˜ a (ξ 1 , • • • , ξ n p ) = n p i =1 ξ i φ i (9)
where ξ i are the associated shape parameters. These latter lie within specific ranges defined in Appendix B . Finally, a shape parametrization ˜ X can be written as:

˜

X (ξ 1 , • • • , ξ n p ) = X + n p i =1 ξ i A φ i ( 10
)
where X is the mean shape.

Appendix B. Shape parameters subspace boundaries

Once the shape parametrization has been obtained, the subspace where the parameters lie needs to be defined. The definition of this subspace is critical. It must encompass most of the shape variation but must not overpredict the solution. A subspace too large can lead to distorted liver shapes or even tangled meshes. In such cases, the iFEMWRAP method used afterward to rebuild the volume is not insured to converge.

The construction of this subspace is based on the knowledge of the shape parameters associated to the target shapes. Each set of n p shape parameters defines a shape ( Eq. ( 10) ) and a point in dimension n p ( Fig. 5 b). First, a point selection is done to exclude exceptional shapes from the training set. The mean Euclidean distance between each point and its k neighbors is computed. Here k = 30 is chosen to get a representative averaged distance from a point to its neighbors while keeping clusters if there are some. All values above a threshold are considered as exceptional and the corresponding points are removed. The threshold is empirically fixed to d + 2 × SD , where d is the mean distance and SD stands for the standard deviation, in order to keep around 95% of the shapes. Visually, removed points correspond to isolated and peculiar shapes ( Fig. 1 ), consequently removing them should not have an important impact for new shapes reconstruction. Then, the boundaries of the parameters subspace are defined ( Fig. 5 c). To create the more compact subspace as possible, the minimum volume enclosing ellipsoids method [START_REF] Moshtagh | Minimum volume enclosing ellipsoid[END_REF] is used. The idea is to generate the n p -hyperellipsoid with the smallest volume in which lie all the training shapes not rejected by the point selection. Each hyper-ellipsoid is defined by its center c ∈ R n p , its semi-axes V i and their associated lengths l i where i ∈ 1 , n p . As the Euclidean distance between the points depends on the dimension, different training sets will be considered for each n p -hyper-ellipsoid.

In n p dimensions the shape parameters ξ i are now described by a radius parameter r ∈ [0, 1] and n p -1 angle parameters θ i such as θ n p -1 ∈ [0 , 2 π ] and θ i ∈ [0, π ] when i < n p -1 . In the reference frame

F = (c; V 1 , • • • , V n p ) , they are related by: ξ 1 , F = rl 1 cos (θ 1 ) ξ k, F = rl k k -1 i =1 sin (θ i ) cos (θ k ) ∀ k ∈ 2 , n p -1 (11) ξ n, F = rl n n p -1 i =1 sin (θ i ) As (V 1 , • • • , V n p
) defines an orthonormal basis, the shape parameters are written in the initial reference frame thanks to the bijection:

T : [0 , 1] × [0 , π ] n p -2 × [0 , 2 π ] → n p -hyper-ellipsoid (r, θ 1 , • • • , θ n p -1 ) → ξ = V ξ F + c ( 12
)
where

V = V 1 , • • • , V n p .
Algorithm 1: SSL algorithm.

Inputs : Shape parametrization ˜ X (ξ) (see Equation 10), n p -hyper-ellipsoid associated to ξ, Additional parameters p, Convergence tolerance tol, Maximum number of iterations i max .

Output : Reduced solution u (x , ξ, p) 23 u (x , ξ, p) = P i -1 (x , ξ, p) [START_REF] Kugler | A new algorithm for volume mesh refinement on merging geometries: application to liver and vascularisation[END_REF] Write u in the canonical tensor format (Equation 3)

1 i = 1 , cv g = 0 2 U = [] //
sinc 1 2 (θ m -θ k m ) cos ( 1 2 (θ m -θ k m )) , if n c even D k (θ m ) = sinc 1 2 n S (θ m -θ k m ) sinc 1 2 (θ m -θ k m )
, if n c odd [START_REF] Gooya | Mixture of probabilistic principal component analyzers for shapes from point sets[END_REF] allowing to write the predicted solution on the n p -hyper-ellipsoid as:

P (x , r, θ 1 , • • • , θ n p -1 , p) = n S k =1 y k (x , p) L k (r) n p -1 m =1 D k (θ m ) . ( 15
)
To insure the interpolation stability the collocation points along the radius are computed using the Gauss-Chebychev-Lobatto points and the angles are equally spaced. The interpolation of parameters p = (p 1 , . . . , p D ) depends on the use case. In this paper, p = b is interpolated using the Lagrange polynomials from Eq. ( 13) . Eventually, the solution is written in a compact format using a low rank approximation as in Eq. ( 3) . As mentioned in Section 1.3 this last step can be done on the fly during the iterations. For the sake of simplicity it is done only once at the end. Algorithm 1 summarizes the different steps and a graphical representation is given in Fig. 4 .

The non-intrusivity of the method allows to use external software to run the computations on the collocation points. Moreover, parallelization can be used as collocation points are independent from one another. Table 3 gives an idea of the number of computations required in function of the number of parameters and the hierarchical level. It must be noticed that interpolating with Lagrangian polynomials and Dirichlet Kernel does not involve the same discretization along the concerned dimensions. The first option will require less collocation points, however the parameters subspace may be overestimated. Whatsoever, for a given number of parameters the number of hierarchical level may be guessed by knowing if the expected solution is regular or not in the parameters subspace. As interpolating functions are used to predict the parametric solution, if this latter has a trend similar to a 12) . Here all angles lie in [0, 2 π] hence some information are redundant. polynomial, for example, then the Lagrange polynomials will be able to quickly recover the pROM. On the contrary, if the solution is expected to vary quickly in the parameters subspace, then the interpolation will be difficult and a high level of iterations will be required. Consequently, the use of the SSL will depend on the smoothness of the solution variations, the number of parameters, the resources needed for a unique computation and the possibility to use parallelization. Some strategies to avoid the computations of all the collocations points in the SSL can be set up but it will not be detailed in this paper, more insights can be found in [START_REF] Borzacchiello | Non-intrusive sparse subspace learning for parametrized problems[END_REF] .

The PXDMF format [START_REF] Bordeu | PXDMF : a file format for separated variables problems[END_REF] is used to visualize the parametric solution in real-time via Paraview. Fig. 11 gives an example with 3 shape parameters plus the breathing one.

ξ

  Vector containing the implicit parameters ( ξ = (ξ 1 , • • • , ξn p ) ∈ R n p ). In the application of this paper they represent the liver's shape. pVector containing the explicit parameters (p = (p 1 , • • • , p D ) ∈ R D ).b Breathing parameter ( b ∈ [0, 1]). X Liver shape parametrization. n p Number of shape parameters used to described a liver shape. n b Discretization of the regular grid containing the control points used in the TPS-PR algorithm ( n b ∈ N * \ {1}), see Appendix A . λr Tuning factor used in the TPS-PR algorithm, see Appendix A . M Number of modes used for the shape reconstruction in the different error measures, see Section 2.6 .

Fig. 1 .

 1 Fig.1. 2D mapping of the liver shape population using PCA along the first and second principal axes. Some of the registered target shapes are displayed and linked to their two dimensional representations. Shapes taken from the dense area are similar, whereas isolated ones are unique. During the construction of the shape parameters boundaries ( Appendix B and Fig.5 c) some of these shapes are excluded in order to generate a more precise and compact pROM. The visualization was done on Matlab.

Fig. 2 .

 2 Fig. 2. Template liver shape. The control points associated to the spline representation used in the non-rigid registration method (see Appendix A ) are shown, the discretization used is n b = 5 . The mesh has been generated using gmsh and the display is done on Paraview.

Fig. 3 .

 3 Fig. 3. Visualization of the liver (in red) with the representation of the portal vein (in blue) for different values of u . The black arrows represent the direction of u d .One of the interest of the procedure is to be able to display internal features to see through the opaque tissues of the liver. The position of the portal vein is deduced from the liver displacement fields computed through the mechanical equilibrium laws ( Eq. (1) ).

Fig. 5 .

 5 Fig. 5. Illustration of the pROM construction for a 2D case. The steps are applied sequentially. The creation of the statistical boundary shape model and the SSL collocation points are represented in the shape parameters subspace.

Fig. 6 .

 6 Fig. 6. Histogram of the volumetric overlap error. This error quantifies the goodness of the registration between two closed surfaces by evaluating the percentage of volume not superimposed. A value of 100 means the shapes are dissociated and a value of 0 means they are perfectly superimposed. Here the mean value is around 10% and the standard deviation around 3.5%. Knowing the quality of the shapes from the database these values can be considered as good ones. The error mostly comes from the smoothing of the surfaces.

Fig. 7 .

 7 Fig. 7.Error made by the projection on the first M modes. The error bar is one standard deviation. This error shows that a relatively high number of modes ( > 50) is required to have a good median reconstruction ( < 4 mm ). Some of the database's liver shapes represented in Fig.1testify of the shape variability of our database and consequently justify the high number of modes required to have a good representation.

Fig. 9 .

 9 Fig. 9. Representation of the error made by the SSL. The value is averaged over the respective training data sets. Two data sets are compared by distinguishing inliers and outliers. The error is also represented for two extreme values of the breathing parameters b . For inliers only the error is very low ( < 1.5 μm ). When outliers are taken into account it increases up to 1 mm for M = 10 modes, which is still suitable for AR applications. In both cases the parameter b add an insignificant error ( < 0.1 μm ).

Fig. 10 .

 10 Fig. 10. Comparison of the shape and the portal vein representations produced by the different steps (the portal vein allows to visually assess the accuracy of the volume reconstruction). The ground truth shape was chosen randomly among the database. The red and blue shapes were reconstructed with n p = 10 shape modes. The breathing parameter was fixed to b = 1 . The projection on a partial set of principal axis provokes a degradation of the representation as evidenced by the slight shift between the yellow and blue shapes. However, the SSL replicates almost perfectly the projection and the difference is not discernible at naked eye. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

5 for j = 1 to n c do 6 on the volumetric mesh previously computed 8 9 end 10 Complete 12 if i > 1 then 13 err = 0 14 for j = 1

 5689101213141 Create empty matrix 3 while cv g = 0 and i ≤ i max do 4 Compute values of the collocation points P i j ( j ∈ [[1 , n c ]] ) Reconstruct the volumetric mesh of the shape given by ˜ Store the solution -i.e. the snapshot -as u j the solutions database U = [ U, u 1 , • • • , u n c ] 11 Get the prediction P i (x , ξ, p) of the solution over the whole parameters subspace by interpolating the snapshots contained in U to n c do 15 er r = max (er r, ||P i -1 (x , ξ P i j , p P i j )u j || )

Fig. 11 .

 11 Fig. 11. PXDMF Paraview plugin to visualize separated variables solutions. A solution with three shape parameters plus the breathing one is represented. The sliders at the bottoms allow to visualize in real-time the result for different set of parameters within their predefined boundaries. The radius and angles defined in Appendix B are used instead of the shapes parameters ξ, which is equivalent as they are linked by Eq. (12) . Here all angles lie in [0, 2 π] hence some information are redundant.

Table 1

 1 List of abbreviations.

	AR	Augmented Reality
	CT	Computed Tomography
	FE	Finite Element
	iFEMWRAP	Iterative FEMWRAP [43]
	IR	Interventional Radiology
	MIS	Minimally Invasive Surgery
	MOR	Model Order Reduction
	MRI	Magnetic Resonance Image
	OR	Operating Room
	PCA	Principal Component Analysis
	PDM	Point Distribution Model
	PGD	Proper General Decomposition
	pROM	Parametric Reduced Order Model
	ROM	Reduced Order Model
	SSA	Statistical Shape Analysis
	SSL	Sparse Subspace Learning
	SSM	Statistical Shape Model
	TPS	Thin Plate Spline
	TPS-PR	Thin Plate Spline Parametrized Registration
	VOE	Volumetric Overlap Error

Table 2

 2 List of notations.

	T	Number of liver shapes in the database.
	v V	Number of nodes in the liver template mesh.
	u	Displacement field associated to the liver mesh.
	u

d User-prescribed motion on the surface of the liver to represent the breathing motion.

n S (θ m -θ k m )

Acknowledgments

This work was performed by using HPC resources of the Centrale Nantes Supercomputing Centre on the cluster Liger, granted by the High Performance Computing Institute (ICI) ( https://ici. ec-nantes.fr/ ). This research was financially supported by BPI France through 3D-Surg, a project of the "Programme Investissements d'Avenir -PIA" of the French Government. ( http://www.3d-surg.eu/ ).

Conflicts of Interest Statement

The authors certify that they have NO affiliations with or involvement in any organization or entity with any financial interest, or non-financial interest in the subject matter or materials discussed in this manuscript.

Appendix C. Sparse subspace learning

The SSL [START_REF] Borzacchiello | Non-intrusive sparse subspace learning for parametrized problems[END_REF] is based on the sparse grid approach [START_REF] Bungartz | Sparse grids[END_REF] . It consists in a sparse and compact representation in tensor format of the space of the solutions. When combined with a dimensionality reduction method -e.g. the PGD -the solution can be written in canonical tensor format as in Eq. (3) . Such formulation allows to store the data in a compact way and to enable fast computations for real-time applications. A nonlinear greedy algorithm is used to compute a reduced order representation for Eq. ( 3) that is one with the smallest possible number of terms d .

The method is a constructive and iterative process. We described it here in a framework where the shape is parametrized (see Eq. ( 10) ) and a set of explicit parameters p are used to define the FE model. First, the subspace where the parameters lie is chosen. It will be the n p -hyper-ellipsoid created in Appendix B for the shape parameters. At each iteration i a point set

is defined following Smolyak's quadrature rule. Note that the value of n c depends on the iteration i . Each point of this set is a combination of n p shape parameters plus D other ones such as n c ) . Thanks to the shape parameters ξ P i j the iFEMWRAP can be used to get the volumetric meshes for each collocation point. A FE model corresponding to the problem defined by the parameters p P i j can then be solved for each P i j . The computed solutions are interpolated over the whole parameters subspace to get the prediction P i . In practice, in the case of n p -hyperellipsoids, the shape parameters are interpolated over the radius and the angle parameters and Eq. ( 12) is used to get their values in the initial reference frame. Given n S data points in the shape pa-

p) the associated solutions, the radius r is interpolated with the Lagrange polynomials described by: L k (r) = n S j =1 , j = k rr j r kr j [START_REF] González | Computational patient avatars for surgery planning[END_REF] and the angle θ m with the Dirichlet Kernel:


Remy Fabbro
email: remy.fabbro@ensam.eu.

C H Caceres

Scaling laws for the laser welding process in keyhole mode

Keywords:

This study shows that the keyhole model derived for determining the scaling laws of keyhole depths for laser welding when high power incident laser beams are used (typically in the multi-kW incident power range), can be also applied to determine the melted depths observed during the Selective Laser Melting process, where much lower incident powers of typically few hundred watts are focused on very small focal spots. The solution of the thermal analysis of this keyhole configuration is described by only three independent dimensionless parameters that are also involved for the analysis of a more general problem of heat conduction using similar input parameters. This global approach and the keyhole model describing the process of laser welding have been also validated by analyzing the melted depths generated by the Selective Laser Melting process. The dependence of the melted depths on the operating parameters of this process has been established, as well as the formation thresholds of the keyhole.

Introduction

Laser welding, which has been used since the 1970s, has become one of the most important laser processes in the industrial world. It allows the assembly of metal parts for a very wide range of thicknesses, from very thick, greater than ten millimeters, thanks to the use of power lasers delivering powers in the multi-kW range, to much smaller thicknesses, less than a millimeter using lasers of much lower power. Many experimental, analytical or numerical studies have made it possible to understand the main physical processes occurring during this very characteristic and complex welding mode, the keyhole (KH) welding, whose melt pool transverse cross sections are characterized by rather large aspect ratios R, defined as the melt pool depth/laser spot diameter [START_REF] Katayama | Introduction: fundamental of laser welding[END_REF]. The prediction or the analysis of these KH depths according to the different operating parameters and materials used is therefore of the utmost importance for this application.

In addition, metal addition manufacturing processes, and more particularly the Selective Laser Melting (SLM) of powders, have also become more recently widely used processes [START_REF] Yap | Review of selective laser melting: materials and applications[END_REF]. They produce molten zones of much smaller dimensions, but which can also present the same characteristics of the keyhole mode observed in welding (so with high aspect ratio), whereas the operating parameters used are very different (typically laser incident powers of a few hundred W, laser spot diameters less than 0.1 mm and processing speeds of the order of m/s). Therefore, the occurrence of the KH mode for SLM conditions has been suggested, but never really proven quantitatively, and the knowledge of these melted depths is also very important because the quality of the resulting densification has been shown to be directly related to the melt pool sizes [START_REF] Tang | Prediction of lack-of-fusion porosity for powder bed fusion[END_REF].

It is the purpose of this article to show that these two very different processes can both be described quantitatively by the same model describing the evolution of the aspect ratio R according to the operating parameters used and the thermo-physical parameters of the material. In a first step, the use of Buckingham's theorem allows to define the number of independent parameters controlling this thermal problem and justifies its solution thus determined, which is recalled. This model is then applied to the analysis of recently published experimental data obtained in SLM conditions, and it is shown that it is possible to reproduce precisely these experimental data as well as the observed thresholds for KH formation.

Scaling law methodology

Definition of the problem

In order to determine the scaling laws that control the KH depth e as a function of the main operating parameters and thermophysical properties of the welded material, several simplifying hypothesis have been used. The KH is assumed to be a vertical cylinder with a diameter d that is equal to the spot diameter of the incident laser beam; this KH is moving with the welding speed V, inside a material at an initial temperature T 0 . The incident laser beam is P and A is the fraction of this power absorbed inside the KH. One knows also that the KH walls must have a temperature at least equal to the vaporization temperature T v ; T So, one will consider that the wall temperature of the cylindrical KH is constant at T v . The used material is also characterized by its heat conductivity K(W/m.K) and by one of the two last thermophysical parameters: its heat capacity per unit volume ρC p (J/m 3 .K) or its diffusivity κ (m 2 /s).

So, in the frame of these hypothesis, one can consider that this problem is totally defined and closed by these previous p = 7 parameters: the resulting KH depth e(m), function of the 3 operating parameters, which are P(W), V(m/s), d(m), the 3 thermophysical parameters: K(W/m.K), ρCp (T/m 3 .K) and the KH wall temperature relative to the initial material temperature (T v -T 0 )(K). This means that there is a unique relation involving these p = 7 parameters.

Moreover, these 7 parameters are only depending of u = 4 fundamental units, which are the mass [M], the length [L], the time [T] and the temperature [K]. Therefore, the Vaschy-Buckingham π-theorem [START_REF] Buckingham | On physically similar systems. Illustrations of the use of dimensional equations[END_REF] states that under these conditions, there must exist p -u = 3 dimensionless independent parameters (π 1 , π 2 and π 3 ), derived from these initial p parameters, which must satisfy a relation f(π 1 , π 2 , π 3 )=0.

Construction of the 3 dimensionless independent parameters

From these 7 initial parameters, 3 dimensionless independent parameters have to be defined (a given dimensionless parameter cannot be derived from another one). There are several possibilities for this construction, by choosing ratios of parameters (or combination of parameters) using for example units of length, velocities, power, power per unit length or volume. A first obvious dimensionless parameter π 1 is the ratio e/d, which is none other than the usually defined aspect ratio R of the KH. A second dimensionless parameter can be defined as π 2 = VρC p d/K, which the ratio of the welding speed V and another velocity such as K/(ρC p d). In fact one can recognize that this parameter is similar to the well-known Peclet number Pe = VρC p d/2K that will be used below instead. A third dimensionless parameter can be π 3 = P a / (dK(T v -T 0 ), which is the ratio of two powers, one is the absorbed laser power P a = A P (with A and P being the absorptivity and the incident power respectively), and the other one results from the combination: d.K.(T v -T 0 ).

So the Vashy-Buckingham π-theorem says that it exists one relation between these 3 dimensionless parameters, which could be written as:

π 1 = R = F(π 3 , π 2 ).
Now, if one adds that the experimentally observed KH depths are usually proportional to the laser power P, one could finally write:

= = = = π R e d π f Pe P d K T T f Pe 1 3 . ( ) . . ( -) . ( ) a v 0 (1)
The relation (1) has been obtained only from dimensional analysis, which is rather efficient because it already gives an interesting scaling law with the parameters of this thermal problem. However, a complete dependence with all the involved parameters is not obtained, because the function f(Pe) is not defined here. We will see in Section 3 how f(Pe) can be only obtained from a complete solution of this thermal problem.

Other possibilities for defining dimensionless parameters

In the previous analysis, the dimensionless parameter π 1 = e/d = R is adapted to the process of laser welding in KH mode, because it involves deep penetrations characterized by large KH depths e compared to KH diameters d (i.e. R > 1). But for laser processes involving heat diffusion, (as surface treatments such as laser hardening or surface melting), where the thickness of the affected material has to be determined, a characteristic normalizing dimension that would be more appropriate to this process of heat diffusion, could be some diffusion length δ obtained during the time τ: δ =(κ.τ) 1/2 , where κ = K/(ρC p ) is the heat diffusivity and τ = d/2 V is a characteristic dwell time. So another dimensionless parameter can then be defined such as: π 1 ' = e/ δ. One can easily see that π 1 and π 1 ' are related through the relation:

π 1 ' = π 1 .2Pe 1/2 or e/δ = 2R.Pe 1/2
(2)

One could also consider another dimensionless parameter π 3 '= ΔH/ Δh(T), which is the ratio of two energy densities, being expressed in J/ m 3 . This approach was initially used by [START_REF] Hann | A simple methodology for predicting laser-weld properties from material and laser parameters[END_REF] for describing the melted depths evolution during laser processing. [START_REF] King | Observation of keyhole-mode laser melting in laser powder-bed fusion additive manufacturing[END_REF] showed that the energy density ΔH represents the energy absorbed A.P.τ during the dwell time τ, which is distributed inside a volume defined by the focal spot diameter and the diffusion length δ: π(d/ 2) 2 δ.

This energy density is finally written as ΔH = (2 3/2 π) 1/2 .4.AP.τ/ (π.d 2 .δ) = 2 3/4 AP/(πκV(d/2) 3 ) 1/2 [START_REF] Rubenchick | Scaling laws for additive manufacturing[END_REF]. (The term (2 3/2 π) 1/2 was introduced so that ΔH is consistent with the initial definition of [START_REF] Hann | A simple methodology for predicting laser-weld properties from material and laser parameters[END_REF] The second energy density Δh(T) can be defined by using the usual enthalpy formulation Δh(T) = ρC p .(T -T 0 ). If the laser welding process in the KH mode is analyzed, the enthalpy at vaporization Δh(T v ) should be used. But one could also use the enthalpy at the melting Δh(T m ), if the melting temperature T m is the main involved characteristic temperature, as for example in surface treatment.

Similarly to Eq. ( 2), one has the following relation between π 3 ' and π 3 :

= = ′ π H Δh T π Pe π Δ ( ) ( 2 .
) .

v 3 7/2 1/2 3 (3)
Finally, it is also possible to define another energy density ΔH O (in J/m 3 ) that is only related to the operating parameters P, V, d such as: ΔH O = AP/(V.d 2 ). Similarly to Eq. ( 3), a dimensionless parameter π 3 '' = ΔH O /Δh(T v ) could be defined. These different dimensionless parameters verify the relation:

= = ′′ ′ π Pe π π Pe π 1 2 . ( 2 .
) .

3 3 11/2 1/2 3 (4)
The different dimensionless parameters involved in this paper are summarized in Table 1: 3. Solution of this thermal model

Analysis of high power laser experiments

With the different hypothesis for the KH thermal model defined in § 2.1, it is possible to determine the KH depth resulting from these conditions and compare it with the previous scaling law of Eq. ( 1). As this determination has already been detailed in a previous publication [START_REF] Fabbro | Analysis and possible estimation of keyhole depths evolution, using laser operating parameters and material properties[END_REF], the main results will be shortly recalled here.

One considers that the absorbed laser power P a is homogenously distributed over the KH wall surface at T v , along the KH depth e. Therefore, if one knows the absorbed power per unit depth P z = dP/dz that is conducted through the KH wall necessary for maintaining the KH surface at T v , the KH depth must simply verify the relation e = P a /P z .

For determining the absorbed power per unit depth P z = dP/dz, we assume a 2D thermal field induced inside the material, because the KH aspect ratio R ≥ 1. Considering stationary conditions, it can be shown that the solution of the 2D heat equation is then only dependent of the Table 1 Possible dimensionless parameters involved in this thermal problem that depends of the 7 parameters : e(m), P(W), V(m/s), d(m), K(W/m K), ρCp(J/m 3 K), (T v -T 0 )(K). The relations between them are given by the Eqs. ( 2)-( 4).

π 1 (=R) (aspect ratio) π 2 (=Pe) (Pe number) π 3 π' 1 π' 3 π'' 3 e/d VρC p d/2K AP/(dK(T v - T 0 ) e/δ ΔH/ Δh (T v ) (AP/Vd 2 )/ Δh (T v )
Pe number, the KH diameter d and the two temperatures T v and T 0 . As the thermal field is now determined, the resulting heat flux per unit depth P z conducted through the KH wall can be computed and is then given by:

= - Pz K T T g Pe . (
) . ( ) v 0

(5)

In Eq. ( 5), g(Pe) is a function numerically determined of the Pe number. As a result, the aspect ratio R of the KH is then given by:

= = = R e d P d P P d K T T g Pe . . .( -) . ( ) a z a v 0 (6)
By using the independent dimensionless parameters defined in Section 2.3, one can see that the relation ( 6) can be rewritten as π 1 = π 3 /g(Pe), which is similar to the relation (1); but now the function g(Pe) is determined.

Remark: The absorbed power P a is related to the incident one P by the relation P a = A.P, where A is the fraction of incident power trapped inside the KH. Due to the cylindrical shape of the KH, the laser beam is absorbed by multi-reflections inside it. The [START_REF] Gouffe | Correction d'ouverture des corps noirs artificiels compte tenu des diffusions multiples internes[END_REF] formula takes into account this process and shows that the resulting absorptivity A depends of the aspect ratio R and varies from A 0 for small aspect ratio (typically A 0 ≈ 0.45, for a metallic surface, with a 1.06 μm wavelength laser) to about 0.8-0.9 for KH's with large aspect ratio R (typically for R ≥ 7-8). For the analysis of macro-experiments (shown in Fig. 1), where R is large, we will therefore consider that A≈1; but for SLM microexperiments, as R is not so important, one should rather consider that A ≈ A 0 .

In order to easily use Eq. ( 6), we have also shown the great interest in considering that within a given Peclet range, g(Pe) can be assumed to be a linear function of Pe such as: g(Pe) = m.Pe + n. The constants m and n are only depending of the considered Pe range (for example, typically m ≈ 2.5 and n ≈ 3 for 0.3 < Pe < 10, which corresponds to the main usual processing conditions for steel alloys, both for laser KH welding or SLM process conditions; for lower Pe, it is observed that m increases and n decreases).

Therefore using this linear relation for g(Pe), Eq. ( 6) can be rewritten as:

= + ( ) R R 1 V V 0 0 (7a) = - = where R A P n d K T T and V n m K d ρC . . . . ( ) 2 v p 0 0 0 (7b)
Eqs. ( 7a) and (7b) give the variation of the aspect ratio R as a function of the operating parameters and the thermophysical properties of the material. One can first remark that Eq. 7-a shows a linear dependence of 1/R with the welding speed V; so this behavior should be easily verified for experimental data.

In addition to these previous results, Eq. (7a) can also be rewritten in a more interesting general form, which can be applied with any operating parameters P, V and d and which separates the role of the operating parameters and the thermophysical one: So for a given material, Eq. ( 8) shows that for any aspect ratio R resulting of the use of different operating parameters P, V and d, the plot of the variable P/(d.R) as a function of (Vd) is a linear function whose ordinate at origin a and slope b are only depending of the thermophysical properties of the material and the range of the Peclet number (through the m and n dependence). When the range of the welding speed V (or the corresponding Pe number) decreases, as m increases and n decreases, one can see that the ordinate at the origin a decreases and the slope b of the linear function increases.

In Eq. ( 8), as previously discussed, one can admit that A ≈ 1, because large aspect ratios are considered here.

This behavior has been verified by analyzing the published data by Suder and Williams (2014) on KH depths obtained on St35 steel at 1.06 μm laser wavelength, for different incident laser powers and focal spots. On Fig. 1, the variable y = P/(d.R) has been plotted as a function of x = (V.d) for these different operating parameters. One sees that the data can be fitted by a unique linear function, y = a 1 + b 1 .x, with a 1 ≈ 3.0 10 5 W/m and b 1 ≈ 2.6 10 10 J/m 3 . By considering that T v ≈ 3100 K, T m ≈ 1800 K, m ≈ 2.5 and n ≈ 3 and using Eq. 9, one finds a mean value of K ≈ 40 W/m.K and C p ≈ 950 J/m 3 .K, which appears to be rather close to what is expected for this kind of steel [START_REF] Mills | Recommended Values of Thermophysical Properties for Selected Commercial Alloys[END_REF].

From Fig. 1 it can be seen that the differences between the experimental data and the representative line of the model using realistic thermodynamic parameters are quite small, which can be roughly estimated at about 20% for the large aspect ratios R. This discrepancy seems to increase with the welding speed (or when the aspect ratio R decreases).

One can therefore see that the scaling of the KH depths given by Eq. (7a) (7b)) (or (8)) is verified for these conditions of laser welding at high power and deep KH's. It has also been verified for Cu laser welding [START_REF] Fabbro | Analysis and possible estimation of keyhole depths evolution, using laser operating parameters and material properties[END_REF], and also for many other data obtained for different operating conditions.

Analysis of KH mode during selective laser melting additive process

Selective Laser melting (SLM) process uses rather low laser powers (typically several hundred watts) and focused laser beams (typically around 100 μm). It is believed that the depth of the molten pool is controlled by heat conduction in the powder bed and the underlying base material. One also observes for high laser powers (or low scanning speeds) a change in the cross section shape of the molten pool: its melted depth becomes greater than the focal spot diameter, which is characteristic of the KH mode.

This transition from conduction to a KH mode occurs above a threshold that has been analyzed by [START_REF] King | Observation of keyhole-mode laser melting in laser powder-bed fusion additive manufacturing[END_REF]. By plotting the melted depth normalized by the beam radius (which is twice the aspect ratio parameter R defined here) as a function of a normalized enthalpy X = ΔH/h s , they observed a first transition threshold characterizing the melt pool formation for X m = ΔH/h s ≈ 10; Then, for higher normalized enthalpies, the KH mode appeared typically for a threshold X v = ΔH/h s ≈ 25-30. For defining their normalized enthalpy ΔH/h s , they used the non-dimensional analysis of [START_REF] Hann | A simple methodology for predicting laser-weld properties from material and laser parameters[END_REF] where h s = ρC p T m (h s is similar to a "melting enthalpy" but rigorously, it differs from the usual h (T m ) melting enthalpy, defined in Section 2.3, because the reference to the initial temperature T 0 of the material is missing; this remark is important if one considers that the initial temperature T 0 may change due to some possible preheating occurrence during the manufacturing process). [START_REF] Rubenchick | Scaling laws for additive manufacturing[END_REF] have also reported a more recent analysis of the evolution of these melted depths e as a function of operating parameters for several experiments on steel, Inconel 625 and Ti-6Al-4 V alloys, realized on different machines. They plotted a different normalized depth (e/δ), as a function of the same normalized enthalpy ΔH/ h s , and they showed that all these data could be collapsed in one curve. This result is interesting because it shows that there is a unique mechanism underlying these data obtained for very different machines, operating parameters and materials. As different parameters of particle diameters and layer thicknesses were used, this means that these parameters are not relevant for controlling the melted depths evolution. This can be explained by the fact that the thickness of the powder bed typically used in SLM experiments is about 100 microns (corresponding to 2-3 particle diameters); therefore with a typical density of the powder of about 40-50% of the solid density, this should give a thickness of liquid (or solid) layer of about 40-50 microns, which is rather small compared to the melted depths of about several hundred of microns observed during this process. For the following discussion in Section 3.3, these different normalized data have been reproduced in Fig. 2.

For the analysis of these results, Rubenchick et al. ( 2018) used a thermal model of the melt pool formation, only based on a conductive process, where the laser beam is absorbed on the planar surface of the sample. With these hypothesis, they showed that the computed melted depth (normalized by δ) is a function of two dimensionless parameters, the previous normalized enthalpy ΔH/h s , and the inverse of the Peclet number Pe. Of course, this is in agreement with our analysis described in §2, which states that the (normalized) melted depth of this thermal problem only depends of 2 dimensionless parameters. Therefore, the fact that the data shown in Fig. 2 can be collapsed in one curve is rather surprising because it shows a relation between only 2 parameters (the normalized depths and enthalpies), and the dependence with the Peclet number does not seem to appear here. The reason of this behavior will be explained in Section 3.3.

Moreover, the use of their thermal model leads to another inconsistency reported by these authors: for high density experiments, leading to observed high aspect ratios characteristic of KH geometries, their simulations show that the peak surface temperature is much greater than T v . For these conditions, it is clear that a pure conductive model is no more relevant for describing the resulting melt pool shape, because for these high temperatures, the resulting induced recoil pressure strongly deforms the melt pool surface and leads to the usual elongated KH geometry. However, our approach described in §2 explains the reason for this similarity between the scaling laws derived from a pure conductive model with a planar melt pool surface, and those derived from a cylindrical KH model (and also those derived from experiments, as we will show in the next section). These are the same three independent dimensionless parameters, or others related to them, which can be used for the different solutions describing the thermal field of these two geometries. As discussed in §2, the use of other dimensionless parameters is possible because they are derived from the first three one by relations involving the Pe number.

We will see in the next Section 3.3 how the results described by the Fig. 2, can be reproduced by using a description of a cylindrical KH model that is moving at high velocity.

3.3. Discussion of the results presented in Fig. 2 Fig. 2 shows that for SLM processing conditions, the normalized melted depths Y = e/δ for 3 different alloys, appear to be a function of the normalized enthalpies X = ΔH/h s that can be roughly collapsed in one linear curve.

Following the analysis of [START_REF] King | Observation of keyhole-mode laser melting in laser powder-bed fusion additive manufacturing[END_REF], the conduction mode appears for X > 10 and the KH mode for X ≥ 25-30. One can see that it is then possible to fit these experimental points of these two regimes by a linear scaling of the form:

Y = Y 0 + α(X -X 0 ) ( 10
)
which have been drawn on Fig. 2. For the conduction mode, Y 0 ≈ 0 and X 0 = X m ≈ 10 (this value of the threshold X m can be more easily determined by using the original results of [START_REF] King | Observation of keyhole-mode laser melting in laser powder-bed fusion additive manufacturing[END_REF]. Similarly, for the KH mode, using Fig. 2, it is found that Y 0 ≈ 2 and X 0 = X v ≈ 25-30 and also the slope of this linear fit is α v ≈ 0.25.

In order to explain these results for the KH mode, we will consider that the melted shape for SLM conditions results from a KH formation whose depth is given by Eq. (7a). But, for usual SLM operating parameters, where typically V ≈ 0.5-1 m/s, one can verify that V > > V 0 (with the heat diffusivity κ ≈ 5 10 -6 m 2 s -1 , and d ≈ 100 μm, typically V 0 ≈ 0.12 m/s and Pe ≥ 5). Therefore, for these operating conditions, Eq. (7a) can be rewritten as:

R = R 0 .V 0 /V (11)
Now, if one expresses the two sides of Eq. ( 11), by introducing the dimensionless parameters e/δ and ΔH/h s used by [START_REF] Rubenchick | Scaling laws for additive manufacturing[END_REF], one obtains after some manipulations:

= = = - R e d e δ Pe Y Pe . (4 ) 2. 1/2 1/2 (12)
and: Combining Eqs. ( 12) and ( 13), one finally finds a linear relation between Y = e/δ and X= ΔH/h s , such that Y = α'.(ΔH/h s ) = α'. X, where the slope α' is given by:

= = - ⎛ ⎝ ⎞ ⎠ = - R R V V T m T T A P V d h π Pe T m T T ΔH h 2. . ( ) . . . . 1 2 . 1 2. . . ( ) .
= - ′ α π T m T T 2 . . ( ) m v 1/2 3/4 0 (14)
If one uses m≈2.5, T v ≈3100 K, T m ≈1800 K and T 0 ≈300 K, one obtains α'≈0.27 from ( 14). This value appears to be very close to the experimental one derived from the linear fit of Fig. 2, where α≈0.25.

So, it appears that it is because the Peclet numbers of the corresponding experiments reported in Fig. 2 are large that these data can be collapsed in one curve given by Eq. ( 11). It also confirms that the melted depths obtained from SLM experiments collected in Fig. 2 for various operating conditions, can be described as resulting from a KH formation, with an adequate thermal field description that considers a KH geometry (This conclusion could already be formulated because of the large aspect ratios commonly observed experimentally of these cross-sections.). Moreover, it is an additional validation of this KH model describing rather small KH's (but nevertheless with rather large aspect ratios) which are obtained for these various SLM operating parameters. Generating a KH under these SLM conditions characterized by high travel speeds is of course possible if one considers that the local intensity of the incident beam is very high (it is typically more than one order of magnitude higher than conventional laser welding). The KH is kept open, with a diameter corresponding approximately to the diameter of the laser spot, only because the very high recoil pressure resulting from this high incident laser intensity balances the surface tension that is very high for these small laser spot diameters.

In addition, because of the high travel speeds during the SLM process and due to the high hydrodynamic velocities induced inside the melt pool, it is observed that the KH no longer remains circular, but is rather elongated in the travel direction [START_REF] Gunenthiram | Experimental analysis of spatter generation and melt-pool behavior during the powder bed laser beam melting process[END_REF]. In this case it is necessary to compare the heat flux from an elongated KH (typically with an elliptical section) with that from a cylindrical KH. This calculation was performed by [START_REF] Miyazaki | Heat transfer from an elliptical cylinder moving through an infinite plate applied to electron beam welding[END_REF] show that for an elliptical KH the normalized heat flux (always defined by g(Pe)) is higher, due to the larger lateral surface area of the elongated KH, and increases linearly with the parameter (γ-1) (γ is the ratio of the major axis of the ellipse in the travel direction (which therefore increases with the welding speed) and the minor axis of the ellipse, which is about the diameter of the laser spot). Typically for γ ≈ 2 the difference between the two configurations is 17%. As the heat flow is greater with an elliptical KH, the corresponding depth should be lower. However, another effect should reduce this increase. Indeed, the surface temperature of the elongated KH is not uniform at T v : it is observed that only the front edge of the elongated KH is irradiated by the laser and only this front edge, located under the laser beam, is therefore brought to the vaporization temperature T v . The rear part of this elongated KH, which is formed by the hydrodynamics of the metallic liquid bypassing the KH front, has a surface temperature that is closer to the melting temperature T m than T v . Under these conditions of nonuniform KH wall surface temperature, it is likely that the heat fluxes from these two different KH geometries are comparable. This may explain why these experimental data can be reproduced using the model with a circular KH.

One can also add that for SLM conditions because of the used high welding speeds, the transfer mode of energy inside the SLM sample is dominated by the convective losses compared to conductive losses.

Finally, Eq. ( 11) can be rewritten in a much more practical equation where the operating parameters are more evidenced. It is easy to see that:

= = - = ′ R R V V m ρ C T T ΔH Δh T ( ) 0.5 ( ) ( ) A P V d p v v 0 0 . . 0 0 2 (15)
Eq. 15 shows that the aspect ratio R is equal to the ratio of two energy densities: ΔH 0 = AP/Vd 2 that only contains the operating parameters, and Δh(T v )' = 0.5 m.Δh(T v ) ≈ 0.5 m.ρC p (T v -T 0 ) is similar to an enthalpy at evaporation temperature (slightly modified through the convective parameter 0.5 m).

Determination of the threshold X v

The threshold X v ≈ 25-30 for KH generation has been experimentally determined by [START_REF] King | Observation of keyhole-mode laser melting in laser powder-bed fusion additive manufacturing[END_REF]. We show here how these results can be derived by using the KH model.

A first usual approach for defining the threshold of KH formation is to determine the operating conditions that induce the evaporation temperature at the surface of the sample. When a laser beam is moving at the surface of a material, the resulting peak temperature T max at the surface is determined by the following relation [START_REF] Hügel | Laser in der Fertigung[END_REF]):

= - + A P d K T T p Pe q .
. ( ). . ( 2)

max 0 1/2 (16)
In ( 16), p = π/8 1/2 and q = 4/π for a Gaussian intensity distribution.

For SLM conditions, as the Peclet numbers Pe are rather large, the relation ( 16) can be rewritten as:

= - A P V d π h Pe T T T . . 8 . . 1 . ( ) s v m 2 1 / 2 0 (17)
Using ( 13) and ( 17), the corresponding threshold X v1 can then be estimated:

= = ⎛ ⎝ ⎞ ⎠ = - X ΔH h π A P V d Pe h π T T T 2 2 . . . . 2 . ( ) v s s v m 1 2 3/4 1/2 2 1/2 1/2 1/4 0 (18)
For T v ≈3100 K, T m ≈1800 K and T 0 = 300 K, one finds X v1 ≈2.3, which is quite one order magnitude smaller than the threshold observed experimentally. This means that reaching the evaporation temperature at the surface of the sample, is not a sufficient condition for the KH generation.

One could also consider that a KH can be defined when its aspect ratio R is greater than one. Therefore, our previous result that gives the aspect ratio R as a function of the operating parameters (Eq. ( 11)) can be used for estimating another threshold X v2 . From the condition R = R 0 .V 0 /V = 1, one obtains:

= - A P Vd m ρ C T T . . . ( ) 2 p v 2 0 (19)
The corresponding threshold X v2 is then given by:

= = ⎛ ⎝ ⎞ ⎠ = - = - X ΔH h π A P V d Pe h π m ρC T T Pe h π m T T T Pe 2 2 . . . . 2. 2 . . . ( ) 2 ( ) v s s p v s v m 2 2 3/4 1/2 2 1/2 3 4 1 2 0 1 2 7 4 1 2 0 1 2 (20)
So, also for T v ≈3100 K, T m ≈1800 K, T 0 = 300 K and m ≈2.5, on finds: X v2 ≈ 7.3Pe 1/2 . As for SLM operating conditions we have seen that the Peclet number is rather large (Pe ≈ 5-10), this would correspond to X v2 ≈ 16-23. This range for X v2 has a much better agreement than X v1 with the experimental determination.

Finally, a third threshold X v3 for KH generation can be defined by using an approach based on the analysis of the KH front geometry: for an efficient beam trapping, the KH front inclination must be such that the reflected beam on this front must be redirected downwards, towards the inside of the sample [START_REF] Fabbro | Analysis and possible estimation of keyhole depths evolution, using laser operating parameters and material properties[END_REF]. This KH front inclination at threshold is therefore 45°(this is also roughly equivalent of considering that R≈1). The condition for obtaining an inclination angle of 45°has already been estimated [START_REF] Fabbro | Analysis and possible estimation of keyhole depths evolution, using laser operating parameters and material properties[END_REF], and is given by the relation:

= A P Vd π H . . 4 2 0 ( 21
)
where H 0 is a modified melting enthalpy derived from the piston model [START_REF] Semak | The role of recoil pressure in energy balance during laser materials processing[END_REF] that takes into account also some losses by heat diffusion. Using (21), the corresponding threshold X v3 is then: For the considered range of the Peclet number, it can be shown that H 0 /h s ≈ 2. So (22) gives: X v3 ≈ 6.Pe 1/2 ≈ 0.8 X v2 (23) For defining the KH formation threshold, Eq. ( 23) shows that the two very different approaches give similar results, which are quite close to the experimental threshold determination. It also means that reaching the evaporation temperature at the surface of the sample, defined by the threshold X v1 , is not a sufficient condition to generate a KH during these SLM operating conditions. The laser intensity should be increased in order to induce a significant depression of the melt pool.

= = ⎛ ⎝ ⎞ ⎠ = X ΔH

Conclusions

In this study, using an approach based on energy conservation of the absorbed power inside a cylindrical KH, an analytical model has been developed that gives the evolution of the KH depth as a function of the operating parameters and the thermophysical data of the used material. The main results can be summarized as follows:

-As a general law, it is obtained that the inverse of the KH depth appears to be a linear function of the welding speed. This point is always verified by corresponding experiments. -The effects of all the operating parameters and thermophysical data of the used material have been clarified. The resulting KH depths obtained with multi-kW lasers can be predicted with a rather good agreement by this model, and it should be considered as an interesting tool in an industrial environment. -It has also been shown that the melted depths observed during the Selective Laser machining process for different operating conditions and materials, can be reproduced using the same KH model. This implies that the involved SLM process can be quantitatively and definitively interpreted as a KH mode welding. -This KH model defines the conditions for achieving the threshold for KH formation. It is shown that the occurrence of the evaporation temperature at the sample surface is not a sufficient condition for achieving this threshold.

Declaration of interest

None

  A.P/(d.R) = a + b.(V.d) (8) where a = n.K.(T v -T 0 ) = n.(P/d) 0 , (9a) (P/d) 0 = K.(T v -T 0 ) is a characteristic power per unit length necessary for a KH, with the used material,and b = m.(ρC p ).(T v -T 0 )/2=m.Δh(T v )/2 (9b)

Fig. 1 .

 1 Fig. 1. For St35 steel and different operating parameters (incident powers and focal spot diameters), plot of P/(d.R) as a function of the parameter V.d and corresponding best linear fit (KH depths obtained by Suder and Williams, (2014) at 1.06 μm laser wavelength).

Fig. 2 .

 2 Fig.2. Normalized melt pool depth e/δ as a function of normalized enthalpy ΔH/h s (from[START_REF] Rubenchick | Scaling laws for additive manufacturing[END_REF][START_REF] King | Observation of keyhole-mode laser melting in laser powder-bed fusion additive manufacturing[END_REF]. The data for the conduction regime have been fitted with the blue line and those for the KH regime with the red one (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).

Acknowledgments

I acknowledge fruitful discussions with Patrice Peyre, Morgan Dal and Matthieu Schneider. This work was performed under the auspices of the Centre National de la Recherche Scientifique (CNRS) at the PIMM laboratory.