Simulation and multi-objective optimization to improve the final shape and process efficiency of a laser-based material accumulation process

Common goals of modern production processes are precision and efficiency. Typically, they are conflicting and cannot be optimized at the same time. Multi-objective optimization methods are able to compute a set of good parameters, from which a decision maker can make a choice for practical situations. For complex processes, the use of physical experiments and/or extensive process simulations can be too costly or even unfeasible, so the use of surrogate models based on few simulations is a good alternative. In this work, we present an integrated framework to find optimal process parameters for a laser-based material accumulation process (thermal upsetting) using a combination of meta-heuristic optimization models and finite element simulations. In order to effectively simulate the coupled system of heat equation with solid-liquid phase transitions and melt flow with capillary free surface in three space dimensions for a wide range of process parameters, we introduce a new coupled numerical 3d finite element method. We use a multi-objective optimization method based on surrogate models. Thus, with only few direct simulations necessary, we are able to select Pareto sets of process parameters which can be used to optimize three or six different performance measures.


Introduction
Modern production technologies follow a strategy of continuous improvement and, as part of it, producers permanently strive to deliver components with high quality at the lowest possible cost. Starting by these premises, it is straightforward to see this as an optimization problem with multiple objectives, considering the input variables as the manufacturing parameters that can be modified in the production line.
In practice, mapping process parameters towards component outcomes can be a difficult procedure when based on experimental tests. Sometimes, experimental runs are expen- Figure 1 Two-level cold forming process: In the master forming step, a preform is generated by melting the end of the wire (thermal upsetting) that is calibrated in a subsequent cold forming step sive and time consuming, making this approach a non-feasible one during production [8]. With appropriate models, numerical simulations can replace expensive experiments.
The combination of optimization methods with computer simulations where simulations are used to transform input parameters (controllable processing variables, or CPVs) into the relevant performance measures (PMs) is an actual engineering need [19,24,27,29]. This requires evaluating optimization functionals on a large amount of candidate solutions, making it a demanding computational task when the simulations are numerically expensive [16]. This is the case for simulating laser-based material accumulation processes (like the thermal upsetting process shown in Fig. 1), for which we have developed a suitable new 3d finite element method (FEM), presented and analyzed in detail in [18]. The numerical challenge is to be able to derive robust results for a wide range of admissible process parameters, therefore a full 3d simulation is necessary that is able to allow quite general geometrical variations of the time-dependent domain and topological changes of the time-dependent liquid subdomain. Here, the key aspects of the method are the two-phase Stefan problem for melting and solidification, an interface-capturing approach for the solid-liquid phase boundary, and the Navier-Stokes equations including a free capillary surface in an ALE formulation for the 3d time dependent domain.
To avoid running a large number of expensive simulations, we pursue a strategy to put together an optimization methodology with surrogate models (also called metamodels), similar to our implementation for milling processes in [19]. In addition to this former application, we are able to use here non-rectangular sets of admissible process parameters, which correspond to reasonable and effective experimental parameter sets. Metamodels are mathematical models (Response Surface, Kriging, Artificial Neural Networks, etc.) that mimic the behavior of the simulation model based on a limited number of observations [3,17,27]. Using metamodels, we reduce the computational effort required to evaluate PMs at a large amount of combinations of the process parameters.
Throughout this work, we present an integrated framework to find optimal process parameters for a material accumulation process using a combination of meta-heuristic optimization models and FEM simulations. In this section, we start by describing the material accumulation process and its inherent optimization needs. In Sect. 2, we present the essential concepts and several details of our 3d finite element method and the numerical evaluation of PMs. Later, Sect. 3 presents the optimization method for solving the multicriteria problems and Sect. 4 contains the setting of two different optimization problems and the solutions obtained for them, being the main difference in the amount of objec-tives pursued during the simulation-optimization cycle. We conclude giving some final comments in Sect. 5.

Laser-based material accumulation process
We consider the process to accumulate metallic material using thermal energy input from a laser to melt material to produce the desired upsetting of a volume portion. This process has been developed within the Collaborative Research Centre 747 at the University of Bremen and experimentally investigated in [5][6][7]. The complete upsetting process consists of the material accumulation followed by a forming step in a closed die or rotary swaging to obtain the final component shape. The correct volume accumulation obtained during the accumulation part is essential to make the forming process possible and the final shape of the component relies on the good quality of the accumulated material. Figure 1 shows the thermal upsetting procedure on a thin steel wire with diameter d 0 < 1 mm. At the beginning, the wire is partially molten from the bottom by a lateral laser beam moving upwards on a path of length s in a shielding gas atmosphere.
Since at the micro scale the surface tension exceeds gravitational force (shape-balance effect [26]), the melt forms a nearly perfect sphere. Due to the one-sided heat application, the solid-liquid interface as well as the sphere itself are tilted in direction of the laser beam. Later, after switching the laser off, overheating of the melt potentially leads to a continuation of the melting process which lasts ideally until a desired height l 0 . At this stage it is desirable that the melt releases the tilted shape and reduces its eccentricity before the temperature decreases and the molten drop solidifies. The latter is desired for the subsequent cold forming step in which the solidified preform is calibrated in a closed die or by rotary swaging, e.g., see Fig. 1.
Now we want to describe and demonstrate the general method applied to a particular process situation. The task is to generate a spherical material accumulation from a thin wire made of 1.4301 steel with initial diameter d 0 = 0.2 mm, for which the desired molten length is l 0 = 3.0 mm. For this process, there are even experimental results by our collaborators from engineering [5,7].
Our goal will be to optimize the shape of the material accumulation and the efficiency of the process, both measured by different PMs described in the next section.

Optimization needs
The nature of the process we are considering includes several parameters that determine the outcomes of the process. Additionally, there is a series of possible PMs that can be defined as our objectives in the optimization procedure.
From the practitioner's standpoint, the input parameters that are commonly used to control the process are the ones related to the application of the laser on the component. Hence, we define our input CPVs to be the laser power P [W], the translational velocity of the laser v [mm/s], and the length of the laser path s [mm].
As PMs we can define several quantities of practical interest and use them either independently or in any combination of them. Generally, we distinguish between PMs regarding process efficiency and quality of the product, which typically compete with each other. Quality is measured in terms of the preforms volume and geometrical shape, which should be spherical and centered below the wire's shaft, while process efficiency is measured with respect to time, energy and the time the laser is used: Geometry related PMs: • minimize the length error | l| with l = ll 0 , defined as the absolute difference between the desired length l 0 and the obtained molten length l, • minimize the distance between the barycenters of the preform and a fictional sphere segment of same volume centered below the wire's shaft, which can be devided into its radial and axial components: radial error (eccentricity), defined as the barycenters distance in radial direction, axial error, defined as the barycenters distance in axial direction, Efficiency related PMs: • minimize process time considered as the elapsed time from switching on the laser until the last molten portion of material solidifies, • minimize the energy spent during the time the laser is on, • minimize the laser usage time corresponding to the time in which the laser is switched on and delivers energy. It is important to remark that these definitions shouldn't be seen as exclusive, as other practical views might lead to different or modified versions of them. Section 2 will present more details about the simulated process and will clarify on the geometrical quantities to be minimized in the geometry related objectives mentioned above. The other three objectives result from the process combination among material, laser and the specific setting of the CPVs. Later on in Sect. 4, results are presented from two case studies with the mentioned input parameters but with three and six optimization objectives respectively.

Process model and numerical simulation
First, we give a mathematical model of the relevant physical conditions including solidliquid phase transitions and free surface melt flow. This is the basis for a finite element discretization of the coupled system of equations. In this work, we use a slight modification of a finite element approach proposed in [18] in order to simulate the material accumulation process. Since the model and simulation approach is elaborated in detail, analyzed numerically and also applied to the accumulation process within [18], the following presentation covers only the main aspects.

Process model
Continuum mechanics are used to model mass and heat transport by coupling conservation equations for mass, momentum and energy. Assuming a sharp solid-liquid interface , denote the time dependent physical domain as sketched in Fig. 2(left). With n {s,l} (t) denoting the outer normal to Γ {s,l} (t), the geometrical evolution of each outer boundary part Γ {s,l} (t) is tracked by a suitable condition for the corresponding normal velocity v n,{s,l} , whereas the interior boundary Γ ls (t) is captured based on the melting temperature T m .
The material movement of particles is described by the material velocity u(t). In the liquid part, fluid dynamics are modeled by the incompressible Navier-Stokes equations including surface tension effects and a kinematic boundary condition for the free capillary boundary. In the solid part, we assume no material movement and a no-slip condition at the solid-liquid interface: Here, p is the pressure, T is the temperature, f u (T) accounts for buoyancy in Boussinesqapproximation, σ is the stress tensor, K is the sum of principal curvatures and Re and We are the Reynolds and Weber numbers. Energy conservation is modeled by an enthalpy formulation for the two-phase Stefan problem in the whole domain Ω(t), in which the non-linear relation between enthalpy e and temperature T is given by the maximal monotone graph where the specific heat c p,{s,l} is assumed constant in each subdomain and Ste denotes the Stefan number. In general, this might yield a mushy region of material that is neither entirely solid nor liquid. To ensure a sharp solid-liquid interface and the definition of a liquid subdomain for the flow problem, we assign this region to the solid subdomain: Here, q {s,l} = κ {s,l} /κ l is a subdomain dependent constant with respect to heat conductivity κ {s,l} and T a is the ambient temperature. Outer boundary heat fluxes include laser heating with power P, a Gaussian laser distribution I L , and absorption coefficient A, thermal radiation due to the Stefan-Boltzmann law with Stefan-Boltzmann constant k SB and surface emissivity E, and finally cooling by the surrounding shielding gas with corresponding heat transfer coefficient δ. The Gaussian laser has a focus diameter of 40 μm and Fresnel absorption is assumed. Other general material and process parameters are listed in Table 1. Due to buoyancy forces, advectional heat transport and geometrical evolution, the subproblems for fluid flow, energy conservation and geometrical evolution are fully coupled. A critical aspect within the model is the geometrical evolution of the solid-liquid-gas triple junction γ (t): The evolution of each adjacent boundary part Γ {s,l,ls} (t) is determined by its corresponding normal velocity v n,{s,l,ls} which overall yields three conditions for the evolution of γ (t). Since the evolution of Γ ls (t) is non-material in contrast to the outer boundary parts, these conditions in general cannot be fulfilled at once. This incompatibility is not tackled in the model and instead addressed in the numerical approach in the following subsection.

Numerical method
For the numerical approximation of the coupled system of nonlinear PDEs, special challenges are given by the melting and solidification resulting in time-dependent liquid and solid subdomains with changing topology, the geometric condition for the capillary free boundary, which results in additional time-dependent changes of the liquid subdomain, and the interplay between both, leading to a triple junction γ (t) between solid, liquid, and outer domain, where the solid-liquid interface Γ ls and the capillary surface Γ l meet.
Several approaches are possible for approximation of the free surface flow, like volumeof-fluid [12,28] or level set approaches [20,22], see also [2]. A mathematically very elegant approach, which nicely fits to the variational FEMs for temperature and fluid flow, is the introduction of a surface-FE method for the geometric condition (2), based on a weak formulation and discretization of the surface-Laplace-Beltrami operator [9]. A weak formulation of equation (2) is derived using Kn l = -Γ x = ∇ Γ · ∇ Γ x where x is the position vector, Γ denotes the Laplace-Beltrami operator, and ∇ Γ the tangential gradient respectively divergence on Γ l (t). Thus, the weak formulation of (1)-(3) with test function v includes a surface integral Γ l (t) 1 We Boundary terms on ∂Γ l = γ from integration by parts vanish, as the position of the triple junction is interpreted as a Dirichlet boundary condition for the capillary surface, given by the position of the solid-liquid interface. For definition of the corresponding finite element spaces on Γ l , it is possible to use the surface trace of the bulk mesh (giving a (d -1)dimensional mesh from a d-dimensional one). Using the P 2 -P 1 Taylor-Hood element for discretization of velocity and pressure, based on an unstructured mesh of tetrahedra S(t), an isoparametric P 2 discretization of the geometry is conveniently usable with the corresponding order of domain approximation. A suitable semi-implicit stable coupling between flow computation and geometric capillary condition was introduced in [1]. The time-dependent domain change of Ω l (t) is covered by an ALE approach, where the motion of the capillary surface Γ l (t) is extended into a smooth parametrization of the whole domain, leading to an artificial advection term in the equation. An implementation was the basis of the FORTRAN-code NAVIER [1]. This was already applied for various applications of free-surface flows in 2d and 3d, but not yet in combination with melting and solidification.
Thus, for simulation of our material accumulation process, a numerical method for this combination of free surface flow with solid-liquid phase transitions had to be newly established: For the discretization of the two-phase Stefan problem, the enthalpy formulation (5) is discretized by adapting a nonlinear finite element method for P 1 elements from [11]. The corresponding implicit tracking of the interface Γ ls is preferable for general situations with variable topology of solid/liquid subdomains. For an interface-capturing approach, where the interface is represented by lower-dimensional mesh facets, topology changes are not natural. In each time step, the subproblems for energy, fluid flow, and geometrical evolution are decoupled from each other, see [18]: • First, one time step of the Stefan problem is solved, using for advection the old velocity field on the old domain, which results in a new energy density and temperature. • Based on the new temperature, we generate a new partition of mesh elements into solid and liquid ones. An element S is marked as liquid if the temperature is above melting temperature, T| S > T m . • In the liquid subdomain, a new velocity, pressure, and position of the capillary surface (with a corresponding ALE deformation of the whole liquid subdomain) are computed by using the weak formulation on the old domain, and using the new temperature for buoyancy forces. For the numerical solution, the nonlinearity of the convection term and the incompressibility constraint are decoupled by an operator splitting. Since we have no material movement in the solid subdomain and the solid-liquid phase boundary is tackled by an interface-capturing approach, the ALE evolution of the geometry effectively reduces to the liquid subdomain.
Full 3d simulations are needed for our application, due to the process design with lateral heating, where, e.g., a 2d rotationally symmetric model is not sufficient for an appropriate simulation. In 3d, with a 2d triangular surface mesh for the capillary free boundary, this combination of 3d bulk FEM with 2d surface FEM is technically much more involved than in 2d, where the outer boundary is just a polygon. For a 2d (or rotationally symmetric) situation with clearly separated melting and solidification times, a direct formulation and discretization of the Stefan condition with a corresponding ALE-formulation could be used during melting [18], but in 3d the situation is definitely more complicated, as an ALE formulation would technically be much more involved, and melting/solidification times are typically not separated, so switching between different methods is not possible. The latter is especially the case when process parameters are varied in a wide range, as it will be needed for the solution of our optimization problem.
In contrast to other free surface flow applications, it is especially important for the 3d simulations needed here for our process optimization, to take care against a degeneration of the tetrahedral mesh. This would typically be caused by the ALE mesh deformation during a drastic change in geometry, like during melting of a thin long cylindrical rod into a spherical drop, as sketched in Fig. 1. This results even in the ALE approach in mesh elements becoming flat and wide over time. To avoid mesh degeneration, an additional remeshing procedure based on the mesh generator TetGen [21] is used and accompanied by the edge collapsing approach proposed in [13] in order to coarsen the outer surface mesh, which TetGen itself is not capable of. Furthermore, the outer boundary of the resulting polygonal mesh is projected onto the isoparametric piecewise quadratic boundary of the old mesh to minimize the shape and volume error introduced by the remeshing. Additionally, each remeshing step needs accompanying transfer operations of data between different meshes. This is especially important for the flow variables, as the interpolation of a weakly divergence free velocity field from one mesh to another does not preserve the weak divergence condition for a new set of test functions. Thus, a subsequent projection is needed corresponding to the new finite element spaces.
At this point, we will only address the discrete evolution of the solid-liquid interface Γ ls (t) in detail due to its major influence on the optimization objectives. For convenience, we will not add additional indices to indicate discrete quantities and use the same notation as before. In the discrete problem, (7) is adapted to the triangulation S(t), as mentioned above, by such that a tetrahedron S can change its phase state within one time step. This results in a discrete evolution of Ω l (t) and Γ ls (t) and we have T| Γ ls > T m as sketched in Fig. 2(right). The definition allows for topology changes of the subdomains and interface during nucleation of an initial melt as well as solidification. Furthermore, the aforementioned overdetermination at γ (t) is resolved since we either have v n,s = v n,l = v n,ls = 0 or a singular jump of γ (t) if a neighboring tetrahedron at the boundary changes its phase. On the bad side, if this jump introduces a kink to Γ l (t), which typically occurs during melting, an artificial capillary wave due to an imbalance between surface tension and inner forces originates at the kink. As shown in [18], the strength of these waves is related to the ratio of time step size and the size of mesh entities in the vicinity of the triple junction γ (t), meaning that a smaller size of mesh entities should be accompanied with a smaller time step size. Further, the fluids viscosity can be artificially increased in order to dampen these waves, which has negligible influence on the heat distribution and final geometrical shape. Altogether, we have an effective numerical method for simulation of the material accumulation process for a large range of process parameters.
The discretization has of course an effect on the optimization procedure. Due to the discrete evolution of Ω l (t) and Γ ls (t) we have a discontinuous change of geometrical quantities. Therefore, a small change in input parameters might lead to a jump in geometry related PMs. The jump size is directly related to the size of mesh entities in the vicinity of Γ ls (t) at the end of the melting phase (see Sect. 4.1), which should be as small as possible. As mentioned before, this would also require a smaller time step size increasing the numerical effort further. In the following, we have chosen the size of mesh entities and time steps in such a way, that we have a good trade off between jump size in geometry related PMs and computing time.

Simulation results
In this subsection, we show some details on the results of the simulated process using the specific case of Run 12 taken from the initial DOE, cf. Table 2 and Table 4 for values of P, For visualization, the bottom tip of the wire is shown at different times in Figs. 3-5. Thereby, the axial dimension of the wire will be denoted by the z-axis with z = 0 marking the bottom tip of the initial configuration. The radial dimensions are denoted by the xand y-axis, in which the laser is impinging from positive y-direction. In order to judge the symmetry of the resulting preform, the symmetry axis (0, 0, z) of the initial configuration is indicated by a dashed line in each figure.
The nucleation and beginning of the melting phase is shown in Fig. 3. Each tetrahedron is shaded in light or dark gray according to it's phase state. This is accompanied by several temperature isolines {T m + Z · 138 • C}, where the melting temperature T m is marked in orange. At t 10 = 3.68 ms an initial melt is nucleated on the laser side without a significant geometrical deformation. Due to the ongoing laser irradiation, melting continues and by heat conduction the whole wire tip melts, with an inclined solid-liquid phase boundary. Hence, due to surface tension, the formation of a spherical melt drop begins that is tilted towards the side of energy input as can be seen at t 25 = 9.20 ms. This continues and results in a nearly quasi-stationary evolution, where more material melts and a growing melt drop is moving upwards. The situation is shown at t 60 = 22.1 ms and t 100 = 36.8 ms. Latter is the last time step of laser irradiation and shown in Fig. 4. The 3d visualization is supplemented by perpendicular cross sections showing the fluid's flow that is predominantly the spherical melt moving upwards. Due to setup of initial geometry and laser the process should be symmetric to the yz-plane. As can be seen from the xz cross section this is only approximately the case due to the non-symmetric unstructured tetrahedral mesh.
Even after switching the laser off, melting continues due to overheating of the melt. We define as endpoint of the melting phase. As shown in Fig. 5(left) melting continues until time t max = t 106 = 39.0 ms in this specific case. Due to heat conduction, melting predominantly continues on the side opposite of the laser which decreases the inclination angle of the solid-liquid phase boundary and hence reduces the tilt of the molten sphere. In Fig. 5(middle, right), we depict isolines of the melting temperature T M from different times. It can be seen that solidification starts at the shaft of the wire and at first, the melt solidifies at the laser side while melting on the opposite side, which further reduces the sphere's tilt. Later on, solidification also starts originating from the free capillary boundary and finally ends at t end = t 667 = 246 ms given by At the end, the preform has a nearly spherical shape which is slightly tilted towards the side of energy input and nearly symmetric to the xy plane.

Numerical evaluation of PMs
The process should ideally form a perfectly centered sphere segment below the shaft in which the molten length l equals the desired value l 0 . The definition of the molten length l is not straight forward since the solid-liquid phase boundary is typically inclined and material can solidify on the one side while melting on the other side. For simplicity, we relate the maximal melt volume |Ω l (t max )| to a cylinder of the same volume and diameter d 0 to match the wires shape. By taking the height of this cylinder, we get a mean position of the solid-liquid phase boundary which is given by and serves as an approximation for the molten length. We want to recall at this point, that Ω l (t max ) is given by a subtriangulation according to equation (9) and hence, the discrete interface Γ ls (t max ) does typically not coincide with an isoline of the melting temperature T m . This way, the molten length is underestimated in the numerical approach depending on the mesh element size as can be seen in Fig. 5(left). To compare the shape of the preform to a centered sphere segment below the shaft, we calculate the center of mass of the preform below l given by computed via a quadrature formula, and compare it to the center of massc of the sphere segment of same volume by as shown in Fig. 6 for Run 12. Due to the tilt of the sphere c y has the largest value and the preform appears smaller on the right hand side. In comparison, c x is rather small due to the symmetrical setup and could be neglected. But it also gives an idea for the numerical approximation error of the geometry. To reflect this, we decided to combine both radial components taking for the radial shape error in case study 2. c z reflects the shape error in axial direction and is only a rough measure for the vertical position due to the calculation of l. Nevertheless it is important, since there are input parameters resulting in only a one-sided melting of the wire. In this case, we have nearly no geometrical change from the initial geometry, which other objectives might not reflect. An even simpler approach chosen in case study 1 is the combination into a single shape error measure by

Multiobjective optimization method
The process we presented in Sects. 1 and 2 involves different PMs that exhibit conflicting behavior. For example, the processing conditions that provide the best quality in terms of product shape may not correspond to the highest energy efficiency. When multiple conflicting PMs are involved, optimizing a single objective can result in solutions that perform poorly for other objectives. Thus, it is not the best approach to obtain a single solution but rather the set of solutions corresponding to the best compromises. For this, we use the following definition for the concepts of non-dominated solutions and Pareto front as in [10] and [24]: Given a problem with conflicting PMs, we can focus our attention on finding the Pareto set, and then a decision maker on a particular moment of the process can select the best solution. This allows for the decision maker to give different weights to the PMs at any time, once the set of non-dominated solutions is known. Alternatively, one could estimate the uncertainty of the Pareto front due to variations in the design space; solutions with less variation are preferred. However, quantifying the uncertainty of the Pareto front requires extensive sampling over the design space, which makes it computationally intensive, [4].
A methodology to calculate the so called optimality influence range of the Pareto approximation due to variations in the design space has been proposed in [14]. The optimality influence range is a hyper-rectangle that encloses all the objective variations with an angle. Binois et al. [4] used conditional simulations of Gaussian process (Kriging metamodels) to estimate the Pareto front and a measurement of uncertainty of the approximated Pareto Front. Despite the fact that the latest uses metamodels to estimate the Pareto front and a measurement of uncertainty it is still computationally expensive.
In this work, we used the metamodel-based multiobjective simulation optimization method introduced in [23] to optimize two case studies for the material accumulation process described in previous Sections.
The method is schematically shown in Fig. 7 and starts by performing an experimental design to generate a set of initial data points, and a simulation run is performed at each point.
Then the set of best compromises between all performance measures is found using Definition 1, and it is called Incumbent Pareto Front.
For the initial experimental design, we use a central composite design, which is one of the most commonly chosen method when the interest is in a complete coverage of the searching space and a regression model is to be fitted. The initial set of points for an n dimensional space is of size 2 n + 2n + 1. The approximation obtained with this design performs very well when the optimization functional does not have strong oscillations and even better when it is a smooth functional (which is the case for the simulation results of our material accumulation process, besides the small jumps in geometry related PMs mentioned in Sect. 2).
Given the incumbent Pareto front obtained from the simulations at the initial data points, the main iteration steps are the following: 1. Use all available simulated data to fit a metamodel for each PM.
2. Use the metamodels to estimate the value of the PMs for a large set of input combinations. 3. Identify the best compromises between all PMs. Call the corresponding Pareto Front, Predicted Pareto Front. The corresponding CPVs settings are the predicted Pareto Set (P set ). 4. Evaluate the predicted Pareto Set using the simulation code. 5. Update the incumbent Pareto Front (based only on simulated data) using the new information. 6. Evaluate stopping criteria. Using these iterative steps, the metamodels are updated and are used to approximate a new Pareto Set. The updated models are able to obtain good approximations of the output responses near the Pareto Front at each iteration. The metamodels make use of all initial points and the new points added at each iteration during step 4. It might happen that the number of solutions on the predicted Pareto Set is larger than the remaining number of total simulation runs allowed (N left sim ), or it is larger than the maximum number of simulations allowed per iteration (N max sim ). If this is the case, a subset of min{N left sim , N max sim } solutions is selected based on a Maximin distance criterion using the predicted Pareto Front.
Jin et al. [15] compared several sequential sampling techniques for metamodeling based optimization. In [15] is recommended to sample 4n (if n < 6) and 3n if n ≥ 6 points, where n is the number of CPVs, at each iteration. Here, the number of points that are sample It has been shown in [23] that this multiobjective optimization method is able to approximate a set of Pareto solutions without having to evaluate a large number of simulations. 15q has been shown to be a good upper limit for the total number of simulations (N total sim ), where q = max{m, n}.

Optimization of a laser-based material accumulation process
For optimization, we consider the example process described at the end of Sect. 1.1, namely the generation of a spherical material accumulation by melting exactly 3.0 mm of a wire with diameter 0.2 mm. We present two optimization case studies of the laserbased material accumulation process. The objective is to produce parts with an optimal shape in an efficient way. Regarding the amount of objectives considered, the first case study considers three simultaneous PMs, while the second case study includes a total of six PMs. Both cases have three CPVs, namely laser power (P), laser path length (s), and laser velocity (v): • With constant laser velocity and path length, a high laser power uses more energy and generates a higher temperature, leading to more molten material, and time needed for cooling is longer. On the other hand, a low laser power may introduce not enough energy to melt the full wire diameter or even any material at all. • The laser path length is connected to the amount of molten material. But, depending on the laser power and velocity, typically a part of the wire is molten which is longer than the path length. • With constant laser power and path length, a lower laser velocity leads to a higher temperature and subsequently a longer cooling time, while a very fast velocity can lead again to the introduction of not enough energy to melt the full wire. The parameter range of the CPVs is based on experimental results in [5,7]. In both works, for a fixed diameter d 0 an optimal laser velocity v opt is determined for several laser powers and a linear regression model is formulated to gain a functional relationship. v opt is in this context optimal with respect to energy efficiency by maximizing the ratio l/E experimentally. It should be noted that v opt barely depends on the laser path length s which is neglected in the regression model. For d 0 = 0.2 mm the linear regression model yields v opt = 1.6P -14.5 (17) for the optimal laser velocity, cf. [5]. Due to v opt 's inherent energy efficiency, the molten sphere is typically barely overheated and the same applies for the solid shaft. After the laser is switched off, solidification starts almost immediately and hence we have l ≈ s but the inclination angle of the solid-liquid phase boundary as well as the tilt of the sphere decrease only slightly. This quality problem might be overcome by decreasing v opt and therefore increase the energy to achieve further melting after the laser is switched off. So in order to attain a specific molten length l 0 the laser path length must be decreased to s < l 0 while decreasing v opt appropriately. To achieve the desired molten length of l 0 = 3.0 mm we consider a laser path length as low as s = 2.0 mm. For this case we have validated by several simulations, that without changing the laser velocity nearly twice the laser power is needed compared to s = 3.0 mm, since energy which gives (17) for s = l 0 = 3.0 mm. Now, we want to build a neighbourhood to define our searching space as a band around this area. For this, we take upper and lower limits given as the ±25% bounds, it is The parameter ranges for P and s were taken as [20 W, 100 W] and [2 mm, 3 mm], respectively. The lower and upper limits of v are defined by v lb and v ub (see Equations (19) and 20), respectively. Thus, the range of v is different at distinct values of P and s and, as a consequence, the optimization search space is not a cube. The shaded region on Fig. 8 shows the search space. By selecting this process window, constructed around the area where experiments suggested an energy efficient process, we restrict our optimization methodology to a relatively small input domain, allowing for a higher precision. We also avoid including solutions far away from the interesting non-linear search space defined by equations (19) and (20) that might potentially mislead the construction of the metamodels.

Simulation setting for optimization
Since all simulations are starting with the same geometry, the cylindrical wire with diameter d 0 = 0.2 mm and length 10.0 mm is only triangulated once. The unstructured uniform tetrahedral mesh is generated using TetGen and subsequently inserted edge midpoints at the outer boundary are projected to the cylindrical shape giving an isoparametric piecewise quadratic approximation. The initial mesh consists of 9806 elements. As mentioned in Sect. 2, the spatial discretization has a discontinuous impact onto the geometry related PMs. Assuming a molten length of 3.0 mm, around 30% of all elements are molten. This results in relative volume change of 0.034% for each newly molten element according to (12), which is only a small jump in the PM for the length error | l|. If an element at the outer boundary is molten or barely not molten during the melting phase, the barycenter c of both corresponding final geometries differs slightly. The difference is related to the edge length on the outer boundary as well as the molten volume. With a typical edge length of 60 μm we observe changes in c xy around 10 μm and even smaller changes in the axial component c z . All these jumps can be reduced by choosing a finer spatial discretization.
Remeshing is performed to avoid numerical issues due to degeneration of tetrahedra. However, this has no significant impact on the number of elements. Time step size has been chosen constant and differently in each simulation, such that the duration of laser heating lasts exactly 100 time steps. This ensures a similar deformation per time step in most simulation runs.
Further process and material parameters are the same as for Run 12 presented in Sect. 2 and have already been listed in Table 1.
The running time of the simulations depends mostly on the flow problem, such that the size of the molten volume has a major impact. Using a workstation with four CPUs (Intel i7-3770, 3.4 GHz) handling up to six simulations simultaneously, runtime ranges between 30 minutes and 10 hours for a single simulation.

Case study 1: optimization of shape and energy
The first case study has three PMs and three CPVs. The optimization objectives, as described in Sect. For the multiobjective optimization algorithm we used the following parameters: as suggested in [23], the maximum number of evaluations allowed was set to N total sim = 3× 15 = 45; the maximum number of runs per iteration was set as N max sim = 3 × 3 = 9; and the lower bound for R 2 was set at 98% (ε = 0.02).
The optimization procedure is as follows: Initialization

Run initial experimental design
First an experiment is designed and run to collect an initial set of data points. As suggested in [25] a Central Composite Design (CCD) with one central point is used. The design was generated as follows: first an inscribed CCD of three factors with ranges [-1, 1] was constructed. Then, the values of P and s were mapped to their original ranges. After that, for each pair of P and s, the range of v was calculated using equations (19) and (20). Finally, the scaled value of v is mapped according to the corresponding range. The values of the CPVs and corresponding PMs are shown on Table 2. Figure 10 shows the initial DOE points as cyan dots. Figure 11 shows the corresponding PMs values using the same cyan color.

Find incumbent Pareto Front
After all data has been collected, the incumbent Pareto Front is identified. The incumbent Pareto solutions, from the initial points, are solutions 1, 2, 4, 7, 8, 11-15. Main Iteration, k = 1

Form a surrogate model per performance measure
At this step, a surrogate model is fitted for each PM using all available simulated data. However, since we know the formula for energy (objective 3), we only fitted a surrogate model for the first two objectives. The fitted models are Multiple Linear Regression (MLR) models with one degree of freedom. The coefficients of determination R 2 of the surrogate models are R 2 1 = 0.9673 ( l) and R 2 2 = 0.9999 (axial + radial error).

Evaluate surrogate models at a uniform grid of input combinations
Both surrogate models and equation (21) are evaluated at a grid of points on the feasible region. The grid of points was constructed with 100 equally spaced levels for each CPV. However, only the combinations that lied within ± 25% of Equation (18) were considered. In total 157,781 solutions were evaluated. Figure 9 shows the evaluation of the models.
The optimization method is implemented in Octave and ran in a PC (Intel i5-8250U, 1.6 GHz). Performing Step 1 and 2 of the optimization takes few seconds, so if we compared it with the time it will take to run thousands of simulations, metamodels are a very efficient way to represent simulation data for an optimization algorithm. Even for the best possible scenario of simulations taking only 30 minutes to compute, a similar global evaluation would need years of computational work to be done.

Find approximated Pareto Set and Front
Now, the Pareto Front of the predicted solutions is found, it consist of 7257 solutions. Since the maximum number of simulations allowed per iteration N max sim = 9, 9 solutions were selected using a max-min distance criteria algorithm with 1000 iterations. This is, 1000 subsets of 9 points were randomly selected out of the 7257 points and the set for which the minimum distance between two points is the maximal was selected. The distances are calculated based on the input values (CPVs). Table 3 shows the input and output values of the 9 new runs. Figures 10 and 11 show the input and output results as dark green circles.

Evaluate Stopping Criteria
Next, the stopping criteria are evaluated. The criteria used here are: (1) stop if the maximum number of simulations allowed was reached (no, 24 < 45); (2) stop if R 2 of all models is larger than 1ε = 0.98 (no, R 2 1 = 0.9673 and R 2 2 = 0.9999); (3) stop if no new Pareto solutions were found (no, new solutions were found). Since none of the stopping criteria were met, a new (main) iteration is needed. On the second iteration, new metamodels were fitted using all available data (24 simulations). The R 2 of the new models are R 2 1 = 0.9999 and R 2 2 = 0.9949 respectively. Later, the models were used to predict a new Pareto front which had 2852 solutions, however only 9 were evaluated. The corresponding input and output values are shown on Table 3 (runs 25-33). Afterwards, the incumbent Pareto front was updated and the new Pareto solutions are 1,4,7,8,12,14,15,16,18,19,21,22,23,24,[27][28][29][30][31][32][33]. Then, the stopping criteria  This concludes the simulation-optimization procedure and delivers a set of solutions that represent the best compromises between all PMs. Now, practitioners can select the

Case study 2: optimization for shape and process efficiency
To further improve the laser-based material accumulation process, we now consider six objectives. Three PMs are related to shape quality and 3 more relevant to process efficiency. Objective 1 corresponds to minimize length error (| l|) In this case study, the same CPVs of Case 1 were used. The range of the CPVs and fix variables was kept the same.
As in Case Study 1, the optimization was conducted following the flow chart in Fig. 7. The optimization parameters used here are: maximum number of simulations N total sim = 15 × 6 = 90; maximum number of simulations per iteration N max sim = 3 × 6 = 18; and the lower bound for R 2 was set again at 98% (ε = 0.02).
The initial design of experiments is the same as in Case Study 1. Table 4 shows the values of the CPVs and the corresponding PMs evaluations. The initial Pareto solutions are solutions 1-4, 6-8, 10-15. After the initial data is collected, a MLR model was fitted to estimate objectives 1 to 4 and objectives 5 and 6 were calculated using equations (21) and (22) respectively. The R 2 of the four metamodels are 0.96729, 0.99866, 0.95088, and 0.99997, respectively. Then, the PM values of 157,781 solutions were estimated (see Case Study 1, main iteration, Step 2 as reference) and the predicted Pareto Front was identified. The predicted Pareto Front at iteration k = 1 had 46,708 solutions, however only 18 were simulated. The selection was performed as in Case Study 1. Solutions 16 to 33 on Table 5 show the corresponding input and output values. The incumbent Pareto front is then updated and the new Pareto solutions are 1-4, 6-8, 10-16, 18-33. Since, none of the stopping criteria are met a new iteration began.
A total of 3 iterations and 69 simulation runs were performed until the method stopped. Table 5 shows the results of the additional runs. Solutions 16 to 33 are for iteration 1, 34 to 51 for iteration 2, and 52 to 69 from iteration 3. The methods stopped because the R 2 of all metamodels was larger that the lower limit.

Comments on the optimization results
For Case Study 1 (Sect. 4.2), we have chosen two PMs for measuring the shape quality of the produced components and one more for the used energy needed on the process. It is important to mention that a process with correct molten length can be achieved by choosing an appropriate velocity v for a laser power P and path length s given within the CPVs feasible region. Further, a fast moving laser with high power usually leads to better energy efficiency, since the time for heat transport in distant regions and boundary fluxes is shorter.
Looking at Tables 2 and 3, runs 21, 28 and 31 give the minimal values of energy. These runs are special, since the applied energy per unit length P v < 0.58 J mm is too low in order to establish the formation of a larger molten volume. Instead, only the region around the laser's impinging spot is melting and solidifying shortly afterwards. The opposite side of  the wire is never molten in this case, resulting in nearly no deformation with respect to the wire's initial geometry. Regarding the other two PMs, this is reflected by a non absolute value l < -2.37 mm and c xyz > 112.3 μm in which the latter is merely determined by the axial component. A molten sphere is formed for P v > 0.66 J mm which is the case in all other runs. First we look at solutions with l ≤ s which correspond to runs 8,12,16,22,23,29,30 and 33. In these runs l is typically a bit smaller than s since solidification starts almost immediately after the laser is switched off. This prevents the melt from releasing it's tilted shape as shown in Sect. 2.3 for Run 12.
With exception of Run 33, which is an extreme case with ls only being -0.046 mm, we have c xyz ∈ [28.42, 81.51] for the shape error indicating now the strength of the final preform's tilt. On the other hand we have the best energy efficiency in these runs with l E ≥ 1.21 mm J (Run 33), which is worse for all other runs. Further grouping regarding c xyz is difficult, since the discretization error for c xyz is typically ±10 μm for the underlying mesh but can go up to about 40 μm in worst cases as can be seen in Run 5. In general, the preform's tilt and therefore c xyz decreases with increasing molten length over the laser path length ls, which as a trade-off results in a worse l E ratio. The worst energy efficiency with l E ≤ 1.0 J mm is achieved by Runs 1, 2, 5, 9, 11 and 13 for which we have ls > 0.44 mm and, due to numerical approximation errors, c xyz < 13.1 μm except for the aforementioned Run 5. For Case 2 in Sect. 4.3, we considered the PMs of axial and radial errors separately (instead as the sum of them as in Case 1) and included two more PMs that are connected to the efficiency of the real process. In total, we have three geometrical PMs measuring the length, axial, and radial errors, complemented by efficiency PMs measuring the process time, used energy, and laser time. It is remarkable how our optimization method succeeds in finding the Pareto front using metamodels reaching values of R 2 > 0.98 and avoiding the simulative effort for combinations of CPVs far away from the Pareto set.
The three dimensional visualizations in Fig. 13 show how for some points in the Pareto front, the values for a single PM might be very large. For example, the axial error shows some extremely large values, but these points are still part of the optimal solutions as they deliver values of almost zero length error and small values on the other four PMs.
Although it might be expected that the last three PMs are similar, the simulated cases have shown that the relations between their values are not easy to correlate. This might be confirmed observing the shape of the Pareto front regarding these three PMs in the left plot from Fig. 13

Relation of metamodel output to experimental results
We want to relate our metamodel output to the experimental results mentioned already at the beginning of this section. v opt , given by (17), has been determined from several experiments in which l ≈ s holds. Hence, we should have | l| ≈ 0 for s = 3.0 mm. In the following, we will relate our results from Case study 2 close to s = 3.0 mm to v opt by considering Pareto solutions with s ∈ [2.8 mm, 3.0 mm] and metamodels for several PMs evaluated at s = 3.0 mm. The simulation's ability to reproduce experimental results is shown in the right plot of Fig. 14. As mentioned before, | l| ≈ 0 should be observed for all pairs (P, v) of simulated In the left plot of Fig. 15 the radial error | c xy | is shown for the same Pareto solutions as before, which is generally increasing for decreasing energy per unit length P/v. The only exception for this is, if P/v is too low to establish the formation of a molten sphere. In this case we have | c xy | ≈ 0 but we have | c z | 0 for the axial error as shown in the right plot of Fig. 15.
Finally, we take a look at Process time, which is the only efficiency related PM we have that requires a metamodel. Figure 16 contains the predicted process time values using the metamodel and the values of the computed process time for the simulated points belonging to the Pareto set with s > 2.8 mm. It can be seen that Process time is highly correlated to the signed distance to the v opt line. The physical interpretation for this is, that overheating the melt above the melting temperature has the most significant impact on the overall process time since melting continues after laser switch-off and subsequent solidification is slower. Furthermore, a high laser power coupled with a corresponding high velocity generally results in a short duration of laser heating and hence a shorter Process time. However, the impact is less significant compared to the effect of overheating.

Summary and conclusion
In this work we have dealt with the simulation based optimization for a material accumulation process. Starting with a general description of our simulation strategy based on physical models for the solid-liquid phase change and the corresponding changes in the geometrical shape of a three-dimensional domain, we presented the details of a new combination of numerical methods to conduct a FEM simulation solving for temperature, evolution of solid-liquid interphase, and also the moving boundaries of the time dependent 3d domain, that is capable to deal with various geometric and topological changes which may happen due to a large range of process parameters.
Further, the optimization method has been described, defining the precise steps to construct an efficient optimization procedure with results for multiple objectives by maintaining only a relatively small amount of FEM simulations. Our method is based on the use of a sequential surrogate modeling that makes use of few simulation results to emulate a general description of the mapping between CPVs and PMs.
We presented the results of applying our methodology to two case studies by considering first three and then six different PMs. The goal was to find the CPVs that optimize all of the defined PMs simultaneously. In general, the method was able to approximate a Pareto Front keeping only a modest number of computed FEM simulations, which is critical for the cases of interest where a single simulation or experimental run needs up to 10 hours to be computed.
Specially, for the second case study, we managed to solve a six dimensional variable space computing only 69 simulations. This might be considered from the perspective of the so called curse of dimensionality in which one might think that on a six-dimensional space approximated by only two points in each dimension, there is a total of 2 6 = 64 points. But describing any function in a single dimension by using only two points leaves plenty of possibilities out of the analysis. Statistically speaking, analyzing a higher dimensional space requires a sampling of the searching space that grows exponentially with the dimension. However, we managed to construct a full search of the space by using only 69 evaluations.
In both case studies we finish the process with a fixed set of optimal soutions that might be used by the practitioners to select during the setup of a real process. In the future, we would like to consider how to account uncertainty on the Pareto front and use it as an additional factor on the decision making process.
Furthermore, we were able to connect the resulting values of our metamodels with already published experimental data where a linear relation of efficient combinations of laser power and velocity was found, cf. Figs. 14 and 15. This shows not only that the simulation offers a good fit with the experiments, but also that the metamodels can be used to describe the experimental values as well.
The optimization of industrial processes via physical models and numerical methods involves typically the solution of coupled systems of (nonlinear) equations. For our application, we have to cope with a coupled system of nonlinear PDEs in 3d, where the time dependent domain is part of the solution. For the optimization procedure, a rather wide range of process parameters has to be studied and the numerical method has to be robust in order to give reasonable results for all admissible sets of parameters. We presented a numerical method that is based on variational principles, especially for the free capillary surface of the liquid subdomain, and is able to derive performance measures for all admissible parameters. Nevertheless, time-dependent 3d simulations consume quite a lot of computing time. Here, the chosen optimization method based on meta-models is able to derive solutions for the multi-objective optimization problem with only very few actual simulations necessary. The framework presented in the previous chapters demonstrates that it is possible to optimize processes involving complex scenarios and corresponding numerical simulations while using only a reasonable amount of computing time. Additionally, the experiments show quite robust results of the optimization method, even with non-smooth data from numerical (discretized) simulations.