A multi-criteria h-adaptive finite-element framework for industrial part-scale thermal analysis in additive manufacturing processes

This work presents an h-adaptive finite-element (FE) strategy to address the numerical simulation of additive manufacturing (AM) of large-scale parts. The wire-arc additive manufacturing is chosen as the demonstrative technology for its manufacturing capabilities suitable for industrial purposes. The scanning path and processing parameters of the simulation are provided via a RS-274 (GCode) file, being the same as the one delivered to the AM machine. The approach is suitable for industrial applications and can be applied to other AM processes. To identify the location in the FE mesh of the heat affected zone (HAZ), a collision detection algorithm based on the separating axis theorem is used. The mesh is continuously adapted to guarantee the necessary mesh resolution to capture the phenomena inside and outside the HAZ. To do so, a multi-criteria adaptive mesh refinement and coarsening (AMR) strategy is used. The AMR includes a geometrical criterion to guarantee the FE size within the HAZ, and a Zienkiewicz–Zhu-based a-posteriori error estimator to guarantee the solution accuracy elsewhere. Thus, the number of active FEs is controlled and mesh manipulation by the end-user is avoided. Numerical simulations comparing the h-adaptive strategy with the (reference) fixed fine meshes are performed to prove the computational cost efficiency and the solution accuracy.


Introduction
Meeting the demands of new design strategies in terms of optimality and functionality of many industrially relevant settings involves the use of suitable additive manufacturing (AM) techniques. When complex geometries are generated via shape optimization analysis, or provided by a multi-scale analysis, where specific micro-scale is required to fulfill certain functionalities (e.g., weight reduction, thermal conduction or noise isolation) [30,42,51,56], traditional manufacturing processes often fail to meet the design requirements. Therefore, the use of AM manufacturing techniques to fabricate such geometries is preferred.
Powder-based AM technologies, as selective laser melting (SLM), electron beam melting (EBM) and selective laser sintering (SLS), produce very thin layers, ranging from 25 to 100 ( μm), while wire feeding technologies, for instance Wire-Arc additive manufacturing (WAAM), and direct energy deposition (DED) technologies, such as laser engineered net shaping (LENS) and electron beam freeform fabrication (EBF 3 ), provide layer thickness in the range of 1000-5000 ( μm). Although the former produce a final product of higher quality [6], they are limited by the actual dimension of the AM chamber. Therefore, WAAM or DED technologies are suitable when dealing with large-size components (where the height can range from tens of centimeters to meters) due their capacity to print thicker layers.

3
Together with the evolution of AM techniques, appropriate CAD tools have been developed to bridge the gap between designers and the manufacturing process, playing a crucial role in the design to production pipeline. Once the computer aided design (CAD) geometry has been prepared, the resulting STereoLithography (STL) mesh is sliced to produce the scanning sequence for each layer [52,53]. Two formats are commonly used: the common layer interface (CLI) [18] and the RS-274 (GCode). CLI contains only the geometry of the scanning path, while the GCode also includes a series of machine tools commands.
Nowadays, the calibration of the AM process parameters is a trial and error process, and experiments are costly. Therefore, a major effort is being put into the development of efficient and accurate finite-element (FE) frameworks for the simulation of AM processes. In addition, from the numerical point of view, complex geometries require several hundreds of scanning tracks to conform the geometry to be read from the CLI/GCode files. As a consequence, a large number of FEs need to be used to describe the computational domain, resulting in a high computational cost [39]. The part-scale meshing step is challenging when complex geometries are involved. Standard solutions adopt an initial FE mesh to accurately describe the geometry resulting in huge FE meshes and large pre-process data files.
Because of this, computational models have been restricted to validation, using simple geometries (e.g., thinwalled structures, cuboids, beams, bridges) with a reduced mesh size [15,29,35,50]. However, to assess the quality and performance of the printed part, high-fidelity (HF) simulations are mandatory to optimize the process parameters. For instance, the accurate prediction of the thermal history is necessary to analyze the metallurgical evolution at the micro-scale [31,32,38]. In this context, the sizes of the FEs and the laser-spot must be comparable [19]; for industrialsized components, this specification implies a dense FE mesh at the heat affected zone (HAZ).
From the computational point of view, problems involving growing domains require a strict control on the evolution of the active FEs and the proper load balancing among processors is crucial. Recently, a scanwise refinement scheme based on the laser position with a fix refinement distance ahead and behind from the HAZ has been proposed by [46,47] to reduce the computational cost. This approach proposes a refined volume exclusively based on the geometrical entities.
Alternatively, a multi-scale approach is performed by [28] using dynamic adaptive mesh refinement and coarsening (AMR) technique to reduce the computational cost. The simulation is split into the part-scale and the small-scale models. In the part-scale model, pre-activation adaptivity is applied to coarsen the FEs inside the current layer group prior to their activation, but it requires an initial fine mesh to describe the geometry. A multi-level method using hp-AMR combined with the finite cell method have been applied by [48] demanding more costly FE approximations.
Other more elaborated approaches use unfitted boundary strategies to properly describe complex geometries, some of them exhibiting the well known small-cut cell problem that yields ill-conditioned systems [7,9,13,14,43,45,48]. In addition, integrating cut FEs may require extra computational cost, stability issues and is not easy to implement [4,10,20,41]. This motivated the shifted-boundary method (SBM) and shifted-interface method (SIM), both substitute the initial domain by a surrogate domain, where the Dirichlet boundary conditions are replaced by equivalent Neumann boundary conditions in a surrogate domain. The challenges in the former approach are the computation and stability of the Neumann conditions for the surrogate domain and also the presence of edges and corners might induce numerical instability and requires a mapping strategy [5,41]. The latter approach makes use of the ideas of SBM, but requires the use of mixed formulations to deal with the jumps in the primal variable and its fluxes [34].
A scalable parallel AM framework is provided in [44] using a geometry criterion based on the Separating Axis Theorem (SAT). The criterion adopted transforms the HAZ into a cuboid and apply the SAT to check if the FEs intersect or not the HAZ [24,27]. All the FEs intersecting the HAZ cuboid are refined and the accuracy elsewhere is achieved by the imposition of the 2:1 ratio, i.e., both coarse and fine FEs must fullfil the following relation: h c = 2 × h f , where h c and h f denote the representative sizes for the coarse and fine FEs, respectively.
To sum-up the main difficulties encountered in the mentioned approaches for simulating the AM process for complex geometries are: (1) the usage of fixed fine meshes is unfeasible to simulate industrial-scale components; (2) some AMR strategies are purely based on geometric parameters without including an error-estimator based on the accuracy of the temperature of the problem; (3) the necessity of using fine meshes to accurately reproduce the boundaries of the domain; (4) the use of complex and costly techniques to integrate cut FEs in embedded approaches; and (5) the imposition of a 2:1 balance scheme is not optimal to keep the number of FEs controlled.
To overcome such difficulties, the AM pipeline for simulating large scale industrial components presented in this article includes the following features: (a) The geometry is defined by using CLI/GCode formats, both widely used in industrial applications. The geometry information is described by a sequence of straight lines, namely, hatches. The laser path describing the 1 3 geometry is the common link between the two formats. The GCode offers additional features for each hatch, such as laser speed, power source and cooling pauses. This enhances the flexibility of the simulation, allowing to modify the process parameters (i.e., scanning speed, power input, feeding rate, waiting time) as delivered to the AM machine. (b) The AMR strategy allows to keep the number of FEs under control throughout the simulation while having a very fine mesh resolution at the HAZ. Simultaneously, the SAT [24,27,44] is performed to find the intersection of the HAZ with a background mesh, and refine the intersecting FEs until the desired refinement level has been reached. In addition, to guarantee the proper use of resources, the proposed strategy periodically performs a load re-balancing process to provide a similar amount of FEs per process in parallel computing. (c) The solution accuracy is guaranteed by an a-posteriori error estimator criterion, on the current state of variables along the simulation. Based on a Based on a Zienkiewicz-Zhu (ZZ) error estimator, the temperature gradient at the super-convergent points of the FE is computed and projected back to the nodes, where the nodal error must be evaluated by the L 2 projection (see. e. g [57][58][59]).
The objective of this paper is to present a robust parallel h-adaptive FE framework applied to the thermal analysis of AM processes. The outline of this paper is as follows. Section 2 introduces the external heat source and the boundary conditions for the WAAM thermal model. Section 3 describes the spatial and time discretizations. In Sect. 4, the main ingredients of the AM pipeline are described. Section 4.1 describes the GCode input-data format. In Sect. 4.2, some aspects about the AMR strategy are detailed. In Sect. 4.3, all the criteria used to generate the coarsening and refinement flags to adapt the mesh are described. Section 5 provides a set of numerical examples comparing the multi-criteria AMR, the layer-wise AMR and the uniform fixed fine meshes. Finally, in Sect. 6, some relevant conclusions are presented.

The continuum problem: strong and weak forms
Let Ω be a growing open bounded domain in ℝ 3 with the boundary Γ . The strong form of the balance of energy equation is given by where Ḣ , ̇r and q are the rate of enthalpy, rate of heat source and the heat-flux per unit volume, respectively. The heat flux, q , is defined through Fourier's law as where T is the temperature field and k is the (temperature dependent) thermal-conductivity.
The phase-change process is generally much faster than the thermal diffusion due to the concentrated laser beam at the HAZ, being its global influence negligible [18], because melting and solidification occur within the same time interval. Thus, the enthalpy rate can be simply written as where C is the (temperature dependent) heat capacity. Introducing Eqs. (2,3) in Eq. (1), yields After applying the divergence theorem to the second term and introducing the test functions , compatible with the Dirichlet boundary conditions, the weak form of Eq. (4) is given by The first integral of the right-hand side of Eq. (5) is the external work of the thermal loads and the second integral is related to the heat losses. Thus, q cond is the heat loss because of the heat conduction between the substrate and the printed component, q rad is the radiation flux through the environment, and q conv is the equivalent heat flux by convection [18,36].
The solution to the AM thermal problem consists in finding T(t) ∈ H 1 ( (t)) in (t i , t f ] , that satisfies Eq. (5) subjected to the initial conditions T(x, t i ) = T 0 (x) , where t i and t f are the initial and final time of the AM process, and the following appropriate boundary conditions.

The external heat source
In the context of AM simulations, the external heat source is generated by the power input (e.g., electric arc, laser beam or electron beam) while following the pre-defined scanning path. To properly track the heat source during the material deposition process, two methodologies are widely used. The first one represents the power input as a Gaussian/doubleellipsoid density distribution [17,26,28,49], while the (1) H = −∇ ⋅ q +̇r, second one consider an average value of the power uniformly distributed in the HAZ [33,37,40]. In what follows, we will consider the latter approach suitable for part-scale analyses and validated with several numerical and experimental tests [18,19]. Figure 1 shows the power distribution as a function of the length-scale, Δl 1 , for the double-ellipsoid distribution, and Δl 2 , for the average distribution, where v scan is the laser scanning speed, Δt 1 and Δt 2 are the time increment for the double-ellipsoid and average approaches, respectively, and d 0 the spot size.

Boundary conditions
The boundary of the domain, Γ , is subjected to the following two heat loss mechanisms: (1) heat transfer by conduction at the contact surface between the substrate and the component and (2) heat transfer by radiation and convection through the surfaces in contact with the environment. Therefore, the boundary Γ is split into two regions, namely: the substratecomponent boundary, Γ s , and the environment boundary, Γ e . The temperature T s is the temperature at the interface between substrate and component on Γ s , thus Heat conduction: The heat loss by conduction between the substrate and component is given by Newton's law [23] as where h cond is the HTC by conduction.
Environment heat loss: The heat loss through the environment is a combination of both convection and radiation heat transfer, expressed by Newton's law as where h conv and h rad are the HTCs associated to the convection and radiation, respectively, and T e is the environment temperature. The HTC for the radiation [19] can be expressed as where is the Stefan-Boltzmann constant and is the emissivity parameter of the radiating surface.

The discrete problem
AM processes involve growing-in-time domains. As a consequence, an initial discrete domain, referred to as Ω(t 0 ) , is taken as point of departure to be continuously adapted along the simulation. This initial discrete domain depends on many factors, among them, the AM technology and the chamber, where the piece is printed, specific for every printing machine. We consider this initial discrete domain Ω(t 0 ) containing the bounding box of the component (domain where the piece is printed), including the substrate domain, where the piece is supported during the printing process, see Fig. 2. This split allows to handle complex mesh patterns generated by the AMR process in the component part, while no adaptivity processes in the substrate domain are performed along the simulation.
To obtain the required mesh density at the HAZ, the initial mesh is subjected to several AMR cycles according to geometrical criteria. Essentially, the initial bounding box is adapted to obtain a proper discretization to capture the temperature dissipation in the melt-pool, see Sect. 4.3.1. Once the mesh has been generated, the FEs are progressively activated following the laser path, as explained later in Sect. 4.1.

Spatial discretization
Consider the AM process defined in the time interval (0, T], divided into N time steps with t ∈ (0, T] . For a given time step, where the increment of time is defined as ( Δt = t n+1 − t n ), the laser beam activates a set of FEs as determined by the scanning path. Thus, at time step t n+1 the domain is decomposed into its active and inactive parts. The active part is defined as a FE partition Ω act ⊆ (Ω n+1 a ∪ Ω n a ) , where Ω n+1 a and Ω n a , are the sets of activated FEs at t n+1 and t n , respectively. The inactive domain is defined as is the domain of inactive FEs and Ω n+1 i the domain of inactive FEs that share nodes with Ω act . Figure 3 shows an example of this classification of FEs at time t n+1 , for a given scanning path in an arbitrary layer.
The domain Ω act is built as the union of FEs intersected by the scanning path. To determine how the domain Ω act evolves along the simulation, the path-length increment per time-step, Δl , is user-defined as an input argument. The strategy will activate all the FEs intersected by the laser path, until the increment of laser path reaches the value Δl , see Fig. 4a. Vertical movements are considered when the laser changes its position from one layer to the next one.
Similar to the classification of FEs, the nodes are also classified according to their status (active/inactive) and their geometric position in the component part. An active node belongs to the component and it is located either at the surface or in the interior of the component. Nodes are classified in the following three sets: (1) skin contains all nodes belonging to the surface of the active portion of the component (i.e., they belong simultaneously to Ω act and Ω n+1 i ), (2) int is the set of active nodes in the interior of the component (i. e. they belong to Ω act ), (3) ina is the set of nodes belonging to inactive FEs ( Ω n+1 i ). Figure 4b shows an example of all sets.
The volume activated between the initial and final position during the time-step, Δt = t n+1 − t n , defines the HAZ as where V n+1 HAZ , V (e) and n e are the discrete volume of the HAZ, the volume of each FE and number of activated FEs inside the HAZ, respectively. The average heat density distribution per unit of volume is computed as where denotes the efficiency of the heat absorption and W the total power input.

Time integration
The implicit backward Euler method is used to advance in time [1,2,16,22]. The time-increment is set automatically according to the metal deposition speed, v dep , and the recoating speed, v rec , which may vary from hatch to hatch. The deposition speed determines the time-increment when the laser is on ( Δt on ), while the re-coating speed allows the computation of the time-increment for the laser to be repositioned for the next hatch, allowing previous printed hatches to experience cooling effects, while the laser is off ( Δt off ); for example, when a given layer is totally printed, and the laser beam is moved to the next one. In some cases, additional re-positioning tasks may be performed, and the re-positioning time must be accounted for. For both laser status, the time-increment, for a constant (user-defined) Δl , is computed as In some cases, the over-heating of the laser require pauses during the manufacturing process. When a given pause is scheduled, the time evolves within a number of time-steps defined by the user, n pause . Then, the time-increment is computed as where no deposition occurs and the laser is set as off, allowing the component to cool down.

AM via hierarchical octree-based adaptive meshes
In this section, a hierarchical octree-based adaptive FE mesh strategy for AM is described. Its objective is to ease the tasks performed by the user, avoiding cumbersome geometry manipulation, and providing a user-friendly environment to tackle complex geometries. The approach enables to eliminate both the user interaction and the need for initial fine meshes to describe complex geometries. As an example, Fig. 5a shows a section cut of a complex geometry embedded in the initial bounding box and surrounded by inactive elements; Fig. 5b shows a top view of the component and the FEs mesh around it, Fig. 5c, d shows the evolution of the mesh following the HAZ for two different times. The first step to start the AM simulation consists of computing the maximum length increment per time-step, Δl , the minimum mesh size (the maximum level of refinement) to describe the HAZ and the maximum mesh size allowed (coarsest level) away from the HAZ. The path data and process parameters are read from the GCode. The number of steps of the simulation is computed dividing the total distance of the scanning path by Δl . This given, an initial uniform refinement (Cartesian Voxelization) is performed. Next, the thermal problem is solved and the error estimators based on the temperature field and the corresponding array of coarsening and refinement flags are also computed. Then, the AMR performs the mesh refinement and data transfer. Next, the FEs activation strategy generates the new computational domain for the next time step.
In large-scale simulations, the use of more than one CPU is required. In this context, once the variables of the problem have been projected to the new mesh and all nodes and related data are redistributed among tasks to balance the load per CPU and guarantee the proper memory usage. Algorithm 1 presents the implementation aspects of the adaptive AM simulation pipeline.
Algorithm 1: Time-stepping and FEs activation in parallel AMR simulation for AM processes. 1 Read input data and generate the initial background mesh (uniformly refined octree) 2 Compute array of time-steps based on the deposition ∆l defined by the user. 3 current time step ← 1 4 Initialize welding path. 5 Activate the FEs for the first step with deposition. 6 while current time step≤ num time steps do 7 Solve thermal problem 8 Compute error estimators based on the temperature field. 9 Compute array of coarsening and refinement flags ft. Compute the array of coarsening and refinement flags based on the geometry of the welding path fg. 13 Generate the final array of coarsening and refinement flags f . 14 Adapt the background mesh Ω b h . 15 Transfer variables to the new mesh. 16 Activation of FEs to build the updated domain. 17 Redistribute mesh and migrate variables among tasks. 18 Generate postprocess.
In Sect. 4.1, some basic features of the input data containing the scanning path geometry in GCode format are described. In Sect. 4.2, some aspects of the AMR techniques applied to hierarchical octree-based meshes are presented. Finally, in Sect. 4.3, the multi-criteria used to adapt the mesh size are outlined.

Input data: the GCode format
GCode is a format broadly used in industrial machine tooling. The advantage of using the GCode format relies on the fact that it can include additional data related to the processing parameters, such as the power input, time pauses, depositing and re-coating speeds, useful to inform the machine and the simulation too. Regarding the GCode commands: G1/G0 correspond to movement commands, where G0 is a fast straight movement (with no deposition), and G1 is a straight line movement with material feeding. Command E provides the information of the deposition process, where the attached value corresponds to the material feeding; if this value is different from zero, the laser is depositing new material in this segment, otherwise the laser is moving without deposition. Command F refers to the laser speed (in (mm/min)). The coordinates defined in a given command line corresponds to the final (x, y) coordinates of the hatch, while the initial coordinate of the hatch are defined in the previous command line. The power source can be changed along the AM process using the command G108 and introducing the P (power source) followed by S (of set) and the value of the new power source in [W]. In some cases, pauses are required 1 3 during the printing process. This time interval can be set by G4 followed by S and the value of the pause in [s]. In this format, the layer thickness is defined by the difference between two consecutive z-coordinates.
As an example, Fig. 6 shows a code snippet describing a laser beam path in GCode format, and Table 1 shows the list of representative properties stored by hatch.
Although the GCode format is most used due to its flexibility, there are other formats also supported by the present framework. This is the case of the CLI format, which only provides the geometry of the scanning path; process parameters need to be input separately.

Adaptive mesh coarsening and refinement
In this section, a brief description of the adaptivity strategy used is given; interested readers may find further details in reference [11].
The hierarchy of refinement levels is defined by an octree parental structure, where zero level FEs are parents to the one level (this way, one level FEs are children of the zero level FEs), and the one level FEs are parent to two level FEs. Figure 7 shows the hierarchical tree of parental relations of the FEs.
The refinement strategy produces a non-conforming mesh with hanging nodes. The treatment of hanging nodes adopted in this work is detailed in [12], and summarized herein; interested readers may find other hanging nodes treatment techniques in [8,9,21,43,54,55]. The modified test/shape function is written as where the modified parent test/shape function, * p , is written in terms of the original (parents) shape/test function, p , and the children shape/test function, c i , according to the spatial location of the child, where x i is the ith hanging child position, and n child is the number of children FEs.
The straight-forward application of the parent/child hanging node constraint may lead to scenarios, where a child node is active and contributes to an inactive parent, causing a singular system. To avoid this situation, the algorithm checks if the parent node is active or inactive. Figure 8 shows the change on the hanging/parent nodes classification in t = t n , Fig. 8a, and in t = t n+1 , Fig. 8b.
Let the zero level mesh (bounding box) be the initial mesh for the AMR strategy. The strategy requires an array of coarsening and refinement flags, , to perform the adaptivity process according to an established criterion, this being performed in an iterative procedure until the mesh does not require further adaptivity or once the maximum number of adaptivity cycles is reached. The multi-criteria strategy to mark the FEs is described in Sect. 4.3.
When mesh refinement is required, an isotropic octree division of each parent FEs into two new FEs in each direction is performed. Thus, each FEs is divided into 8 new FEs in one cycle. Figure 9 shows an initial mesh (just one

Multi-criteria AMR
This AMR strategy is developed in such a way that different adaptivity criteria can be used in the same simulation. This design takes advantage of the Object-Oriented programming, where a set of objects provide an array of coarsening and refinement flags, , to be compared to obtain the most critical scenario. The criteria used to adapt the mesh are described in the following. In Sect. 4.3.1, a geometrical criterion based on the intersection of the FEs mesh with the bounding box of the HAZ is presented. In Sect. 4.3.2, a criterion based on a-posteriori error estimators of the temperature is presented. Each criterion performs a search of FEs to mark them as: to refine (1), to do nothing (0) or to coarsen ( −1 ). The outcome of this process is an array that contains the coarsening and refinement flags for all FEs. This array is sent to the library in charge of performing the AMR process [11], where all related data structures are also updated.
For the temperature field, a point-to-point projection is used. This procedure is straightforward for nested octreebased meshes; temperature at the hanging nodes is interpolated from the corresponding parent elements.

Geometry-based adaptivity
The geometric criterion consists in a search algorithm, where an oriented reference volume intersects the FEs of the background mesh Ω b h . The oriented volume (OV) is a fictitious prism, containing the HAZ at t n+1 , being immersed in Ω b h . The dimensions of the prism are: the hatch length, the melt-pool width, b pool , and the HAZ height including the current layer thickness and the laser penetration (melt-pool depth), h pool .
The intersection between the OV and the FEs rely on the separating axis theorem/hyperplane separating theorem (SAT/HST) [24,27,44] which states that a pair of convex polytopes, formed by E edges and F faces each, intersect if they overlap at least in one of their projections onto their E 2 + 2F planes. When both polytopes are box-shaped, i. e. they contain parallel faces and edges, the number of planes to be tested is considerably reduced. Thus, E = 3 and F = 3 , resulting in only 15 testing planes. In case of using tetrahedra, the number of testing planes is increased to 44 ( E = 6 and F = 4).
In practice, this test compares the projected distance of the pair of polytopes centers, s, with their projected halfsides, r OV and r FE , onto the normal direction of the test plane, tp . If s > r OV + r FE , the plane is a separating plane, hence the polytopes do not overlap. If the test finds a separating plane, no further checks are required, and both OV and FEs do not intersect; otherwise a new test plane should be verified. The method chooses as testing planes, 3 independent planes describing the OV, 3 independent planes describing the FEs of the Ω b h , and the 9 cross-product between them. Figure 10 presents a visual description of this method.
If the OV and the FEs overlap, then the tested FEs belongs to Ω n+1 a , and is marked as to be refined ( f i = 1 ). If a separating axis is detected, and the tested FEs does not belong to Ω n+1 i , this FEs ∈ (Ω n a ∪ Ω n+1 i ) is marked as to be coarsened ( f i = −1 ), see Fig. 3. The case where the FEs ∈Ω n+1 i , then it is marked as to do nothing ( f i = 0 ), to prevent it from being The number of adaptivity cycles depends on the accuracy desired to represent the OV. The refinement level per FEs, referred to as Lev(FEs), is computed and compared with the minimum, min , and maximum, max , respectively. Then, the corresponding coarsening and refinement flag for the ith FEs, f i , is determined as follows In consequence, the number of computations required is bounded by the difference between the maximum and minimum refinement levels [44]. Figure 11 shows a 2D example, the OV inside Ω b h and the final mesh using the SAT test after computing 5 adaptivity cycles with and without the 2:1 balance between FEs.
Alternatively, another geometry-based strategy is the layer-wise AMR, which consists of keeping refinement in the latest set of layers from the current deposition layer to a fix distance down in the building direction. Therefore, the finest mesh size is preserved for the latest layers to avoid continuous remeshing, thus minimizing the data transfer and CPU time.

A-posteriori error-estimators adaptivity
To avoid that the coarsening process near the OV affects the accuracy of the solution, an a-posteriori error-estimator is used [57]. This section presents the gradient-based approach adopted in this work. The discretization error, , for a given discrete solution can be computed as where ∇T and ∇T h correspond to the temperature gradient obtained with analytical (exact) and the discrete FE solutions, respectively. However, exact solutions are only known in very simple cases. Thus, Zienkiewicz et al. [58,59] present a methodology to obtain an accurate error-estimate when the exact solution is not available. The authors show that there exist a set of super-convergent points within the FE discretization, where the solution can be used as reference value when computing the error estimator. These values are a-posteriori computed at the FE super-convergent points and projected to the FE nodes. Thus, Eq. (17) is modified as where ̄ is the discrete error estimator for the gradient of temperature, ∇T is the gradient of the temperature computed at the super-convergent points, P(⋅) is a projection operator that projects ∇T onto the nodes of the FE mesh, and ∇T h is the gradient of temperature computed with the FE solution. Fig. 11 a OV embedded into the initial domain; final mesh after 5 adaptivity cycles, b with 2:1 balance; c without 2:1 balance As a matter of example, in the standard hexahedra tri-linear FE, the super-convergent point is the barycenter of the FE.
Finally, the L 2 -norm of the ̄ for the ith FE, i , is computed and compared with the minimum and maximum admissible errors, as follows

Numerical examples
In the examples presented in this section, the WAAM thermal problem is simulated. Nevertheless, the multi-criteria approach is technology independent and can be adapted to other AM processes.
The objective of the presented numerical examples is to demonstrate the suitability of the multi-criteria approach to industrial part-scale simulations with regard to the computational cost-efficiency and accuracy. Figure 12 shows the temperature dependent properties of the Ti6Al4V titaniumalloy used in all the numerical simulations.
The Cartesian adaptive FE meshes are obtained from an uniform Cartesian Voxelization of the initial Bounding Box, which include the part geometry. The minimum and maximum L 2 threshold are 10 −4 and 10 −3 , respectively. The fixed fine FE meshes are obtained to match the FE size of the maximum refinement level of the adaptive mesh, which corresponds to the resolution inside the HAZ.
To assess the global error of the multi-criteria approach with respect to the fine mesh solution, a relative L 2 error norm of the temperature field is computed at the end of the simulation as is the point to point projection of the reference (fixed fine mesh) solution onto the adaptive mesh and u coarse h is the solution obtained with the adaptive refined mesh.
The performance and numerical accuracy among the multi-criteria AMR and different mesh strategies are assessed and compared. Section 5.1 compares the multi-criteria AMR with a reference fix fine mesh and the AMR using only the geometric criterion. In Sect. 5.2, the multi-criteria AMR is used and the results are compared with experimental data and with the layer-wise strategy. Finally, Sect. 5.3 simulates a complex geometry modelled automatically from the GCode provided. Fig. 12 Temperature dependent properties for the Ti6Al4V. a Density; b heat capacity; c thermal conductivity 1 3 The numerical simulations are carried out using the inhouse software FEMUSS (Finite-Element Method Using Sub-grid Scales) developed at the International Center for Numerical Methods in Engineering (CIMNE). All the computations are performed in a distributed-memory environment using an Intel i7-10700 @2.9GHz with 16 CPU's (8 cores and 2 threads-per-core) and 32 GB RAM memory. The thermal-problem is solved using the bi-conjugated gradient solver together with the Additive Schwarz Method. The postprocessing is done using Paraview [3].

Printing of a cube
In this section, the 3D printing of a cube of side a = 100 (mm) is simulated to study the accuracy and computational efficiency of the proposed multi-criteria AMR. To compare the solution of the approach, the numerical results obtained with the fixed fine mesh are compared to the AMR using only the geometric criterion. The bottom corner of the cube is chosen to compare the local accuracy of the AMR strategies with respect to the fixed fine mesh.
The domain has been discretized into hexahedral FEs, a total of 7 uniform AMR cycles are performed to obtain a proper discretization to capture the phenomena in the HAZ (this discretization corresponds to 2 7 = 128 FEs per side). In this case, the FE representative size inside the HAZ is h = 0.78125 (mm). This mesh contains a total of 2,097,152 hexahedral FEs. Figure 13 shows the fixed fine mesh, and the AMR using the multi-criteria and the geometric criterion, respectively. The difference of both AMR solutions compared to the fixed fine mesh is remarkable; also, the influence of the a-posteriori criterion can be appreciated when compared to the purely geometry-based AMR (see Fig. 15).
The convective flow inside the melt-pool is considered by increasing of the thermal conductivity above the melting temperature [25]. The heat loss by radiation and convection are considered at all external surfaces and the environment temperature is set to 20 ( • C). The radiation emissivity is = 0.7 . The HTC for convection increases linearly from h conv (20 • C) = 10 (W/(m 2 ⋅ • C)) to h conv (1000 • C) = 100 (W/ (m 2 ⋅ • C)).
The cube is built in a sequence of 16 layers, where the power source is W = 5 (kW) with a heat absorption efficiency of = 0.75 . The simulation accuracy depends on the melt-pool deposition process which is set to Δl = 12.5 (mm), defining the length of the melt-pool, and its base and height b pool = h pool = 12.5 (mm). In this experiment, a constant speed for both, deposition and re-coating is considered, v rec = v dep = 3600 (mm/min). Figure 14 shows the evolution of the number of active FEs for the length-scale adopted in the simulation in total of 2087 time-steps and the dashed lines are the average number of active FEs of the AMR simulations. Figure 14a shows the evolution of the three strategies; the difference in the final number of active FEs of the fixed fine mesh compared to the AMR strategies is remarkable. Figure 14b shows the initial 200 time-steps, where the number of active FEs of the fixed fine mesh develops at a much faster rate than the AMR strategies. Figure 14c shows that the AMR strategies increase the number of active FEs initially and stabilize around a much smaller value than the fixed fine mesh, keeping the size of the problem controlled throughout the simulation. The multi-criteria approach has an average number of active FEs of 10,500, almost twice the size of the purely geometric criterion (5400 average active FEs). This reduction represents 1.00% and 0.51% of the average total number of active FEs of the fixed fine mesh at the end of the simulation for the multi-criteria AMR and the geometry-based AMR, respectively.
The real-time run for the fixed fine mesh case is 7 (h) 24 (min). The computational time reduction for the multicriteria AMR and the purely geometric AMR are 88.57% and 91.21%, respectively. Table 2 shows the time spent by each module.   Figure 15 shows the evolution of the temperature at the bottom corner of the cube to assess the local accuracy of both AMR strategies. In the case of the thermal problem, the accuracy is enhanced by the a-posteriori error estimate at a computational cost around 3% higher than the purely geometric AMR for this example. Figure 16 shows the relative difference of the temperature field of the AMR strategies with respect to the fixed fine mesh at the end of the simulation. The local maximum relative differences observed are 2.47% and 21.33% for the multi-criteria AMR and the purely geometric AMR, respectively. Applying Eq. (19), the global error computed for the multi-criteria AMR and the geometric-based AMR are 0.81% and 6.71%, respectively.

Printing of a thin-wall
In the next example, the computational cost and accuracy of the multi-criteria AMR and the layer-wise AMR are compared with the experimental data obtained from the printing of a thin-walled component using the WAAM process.
In the experimental setting, a thin wall is built upon the deposition of 18 layers of thickness t = 4.5 (mm). Figure 17 shows the geometry of the thin wall and the substrate. The metal deposition goes from the initial to the final position with a constant speed of 100 (mm/min) and after concluding the current layer, the re-coating velocity is 514 (mm/min). In this experiment, two scenarios have been considered: a cooling interval between layers of (I) 5 min:21 s and (II) 10 min:21 s (see Fig. 18).  The GCode provides a variable value for the power source for the first 3 layers and thereafter the power is kept constant for the remaining 15 layers. Table 3 shows the power source, while the heat absorption efficiency is set to a value of = 0.75.
The chamber temperature, measured in-situ, is 17 ( • C) and a radiation emissivity of = 0.7 is considered. The HTC for convection is taken as h conv = 10 (W/(m 2 ⋅ • C)). To consider the contact of the bottom surface of the substrate with the stainless steel work-bench, the conduction HTC value is set to h cond = 100 (W/(m 2 ⋅ • C)). The time-stepping is defined by a fix movement of the heat source of Δl = 10.0 (mm), the melt-pool base, b pool = 20.0 (mm), and the height (layer thickness and HAZ penetration), h pool = 9.0 (mm).
The initial mesh consists of two bounding boxes: one to include the thin-wall and one for the substrate. The dimensions of the initial bounding box for the thin-wall and for the substrate are 220 × 80 × 81 (mm 3 ) and 220 × 80 × 20 (mm 3 ), respectively. The maximum and minimum adaptivity level considered in this case is 8 and 4, respectively. In addition, an alternative AMR strategy is forced to keep the finest mesh throughout the last 3 layers, while the mesh below it is coarsened to the minimum refinement level chosen. Figure 19 shows the initial mesh and the different adaptive meshes for each strategy at an intermediate time step.
The thermo-couples are positioned 5 (mm) from the thinwall deposition starting point, at the top of the currently deposited layer. The temperature is measured 21 (s) prior to the deposition of the next layer. Figure 20 shows the temperature evolution of the two AMR strategies compared to the experimental data at the measurement points of each layer for the two idle time scenarios. The temperature evolution from the layer-wise AMR and the multi-criteria AMR show good agreement with the experimental results. The average errors associated to the 5 (min) and 10 (min) curves are 4.8% and 4.4%, respectively. The small width of the thinwall makes the positioning of the thermocouples difficult in the experimental setting. Regarding the 5 (min) test, the thermo-couple was misplaced in the initial three layers and the data is not available.
The real-time run for the layer-wise AMR case is 10 (h) 34 (min), and the multi-criteria AMR yields a computational time reduction of 87.37%. The saving in time is a consequence of the reduced number of FEs required by the multicriteria AMR strategy (see Fig. 21).
The peak value of active FEs on the initial steps of the simulation is due to the layer-wise strategy, refining the substrate with a very dense FEs mesh. The peak number of active FEs is 8,804,331 for the layer-wise AMR and 526,541 for the multi-criteria AMR. After the refinement height of the layer-wise strategy leaves the substrate (around time step 200) the effect on the substrate is reduced and the number of active FEs remains almost constant. For the layer-wise strategy, the average number of active FEs is 1,166,199, being highly influenced by the initial steps, whereas the multicriteria strategy average number is 90.81% lower, being  close to the time reduction in the simulation obtained by the multi-criteria. The requirement of the layer-wise AMR strategy of having 3 adapted layers below the current deposition height is far from optimal when the manufacturing process is capable of depositing thick layers of material, compromising the computational cost due to the large number of active FEs. Even when the refinement height is no longer affecting the substrate this strategy presents almost 5 times the number of active FEs compared to the ones obtained with the multicriteria approach. Another shortcoming of the layer-wise approach is that the mesh remains constant during pauses, whereas the multi-criteria AMR reduces the number of FEs based on the data provided by the cooling process, avoiding excessive coarsening that may affect the accuracy of the model.

Printing of an optimized door handle
In the next example, a complex geometry taking from a topological optimization analysis, is analyzed using both the fixed fine mesh and the multi-criteria AMR.
Due to the complex geometry of the door handle, the structure walls are very thin and hollow, in consequence a balance between computational saving time using mesh coarsening techniques and a proper description of the shape must be guaranteed.
The goal of this example is to compare the performance of the fixed fine mesh and the multi-criteria AMR in a scenario, where the coarsening speed-up benefits are constrained by the complexity of the geometry.
The printed component geometry is obtained from the STereoLithography file and converted to a GCode file using the open-source software CURA. Hence, no user manipulation is required and the mesh resolution depends on the maximum refinement level adopted for the geometric-based AMR criterion. Figure 22 shows the door handle geometry and the points P1 and P2, of coordinates (37.73, 101.48, 2.99) (mm) and (120.0, 117.31, 52.02) (mm) are chosen to measure the temperature evolution for both FE meshes.
The discrete domain is built departing from a bounding box of dimensions 180 × 135 × 70 (mm 3 ) being generated with the data provided by a GCode file. The minimum and maximum refinement levels adopted are 4 and The environment temperature is 20 ( • C) and a radiation emissivity of = 0.7 is considered. The HTC by convection is fix to a constant value h conv = 10 (W/(m 2 ⋅ • C)). In this case, we considered the laser movement defining the time-marching scheme is Δl = 800.0 (mm), the melt-pool base and height are b pool = 500.0 ( μ m) and h pool = 100.0 ( μ m), respectively. Figure 23 shows the evolution of the component building and the inactive background mesh required. Figure 24 shows a cross section of the door handle for both FE models; note that the difference between the finest and the coarsest FEs in the multi-criteria AMR is of 2 refinement levels. Figure 25 shows the evolution and the average number of the active FE. Ideally, the AMR strategy exhibits a fluctuation of the number of active FEs around a constant threshold, as observed in previous examples using simpler geometries. In this case, the complexity of the geometry of the hollowed thin-walled door-handle induces a geometric constraint to the coarsening technique inhibiting higher coarsening capabilities. For the fixed fine mesh, the average number of active FEs is 764,923, whereas the multicriteria strategy average is 57.89% lower (322,099 active FEs). The real-time run for the fixed fine mesh case is 14 (h) 6 (min) obtaining a time reduction of 29.91% when the multi-criteria AMR is used. Figure 26 shows the temperature evolution at points P1 and P2, showing the agreement between the two strategies. Figure 27 graphically shows the relative percentual difference of the temperature field of the multi-criteria AMR is 0.28%. The strategy produces very accurate results at a fraction of the computational cost, even when coarsening is restricted by geometric constraints due to the complexity of the component.

Conclusions
In this work, a robust parallel h-adaptive framework to simulate AM processes focusing on industrial part-scale simulation is presented. The motivation arises from the necessity of the industry to take advantage of the capabilities of AM simulation to produce complex optimized solutions to enhance the performance of industrial components. In this sense, the use of the GCode file to define the geometry of the components is crucial, because it is a flexible format, where the process parameters may be introduced and changed throughout the simulation. The use of the GCode removes both CAD cleaning and meshing operations, very challenging when dealing with complex geometries. Therefore, the proposed strategy allows for a user-friendly environment for industrial part-scale analyses.
The scalability of the industrial use of AM processes depends on the reduction of the printing time for large-scale components. The success on the numerical simulation depends on its capability to reduce the overall simulation time without loss of accuracy. For this reason, the multicriteria approach has a geometrical-based criterion to identify the current position of the laser, creating a maximum refined area at the vicinity of the HAZ by means of the separating-axis theorem. All FEs belonging to the HAZ are activated on the current time-step and the geometric resolution depends on the user-defined maximum level of refinement. The boundaries of the domain are represented by a fitted FE mesh, which does not require cut FE integration schemes. The accuracy is assured by the a-posteriori error estimate based on the temperature gradients, avoiding the imposition of a 2:1 balance scheme.
The numerical examples show compelling results in terms of computational cost reduction and accuracy compared to alternative approaches to simulate AM processes. The multi-criteria approach is easy to implement and can be applied to several others AM processes.
The AMR strategy presented herein is now being developed into a thermo-mechanical setting. To make it applicable for industrial needs, a ROM of the thermo-mechanical problem may be incorporated to speed-up the computational model.