The process of metal additive manufacturing involves a highly localized heat source which moves over a substrate. Its purpose is to change the state of the deposited metal either from powder to liquid or from solid to liquid such that the added material bonds with the substrate. The subsequent, rapid cooling of the heat affected zone induces undesired residual stresses, a process well understood in welding, see for example Goldak2005 ; Lindgren2007 . The numerous physical phenomena involved in this process range from the simplified view stated above to models including detailed weld pool dynamics Tanaka2014 , to models resolving the micro-structure evolution in the cooling phase (see e.g. Zinovieva2018 ; Basak2016 ; Korner2016 ). The necessary compromise between the complexity of the physical model and the ever increasing yet limited resources for its numerical resolution is mainly driven by the concrete question which needs to be answered by the simulation, (see e.g. Lindgren2006 for a guideline in welding).
In this context, efficient numerical discretizations are desirable of which the finite element method is the most prominent choice. A simple measure of its complexity is the number of degrees of freedom involved in the computation. These represent the unknown coefficients of piecewise polynomials spanned on a mesh of finite elements which resolve the highly localized gradients. Three basic strategies are available to control the number of degrees of freedom: a)-refinement, i.e. an increase of the number of finite elements in the mesh, b) -refinement, i.e. an increase of the polynomial degree of each finite element or c) the combination of both: -refinement. The achievable rate of convergence increases from a) to c) for problems with locally high gradients. Unfortunately, so does the complexity of the implementation of the respective method. This is especially the case in transient problems where the necessary adaptions need to be kept local to the (traveling) heat affected zone.
Numerous approaches are reported in literature for locally -refined discretizations, see e.g. Lindgren1997 as a representative of an early work on the subject in the context of welding and Schoinochoritis2017 for a recent review of the available variants applied to metal additive manufacturing including a review of commercial packages for this purpose. Refinements in beyond quadratic shape functions are scarce, but e.g. Nubel2007 demonstrates that it is possible to construct efficient discretizations for elasto-plastic problems by using high order discretizations whose boundary follows the plastic front. The application of -adaptive methods in this context is e.g. discussed in Olesky2015 where exponential convergence rates were achieved for boundary-conforming discretization. Recently, high-order refinements have also become of interest in the context of Isogeometric Analysis, see e.g. Elguedj2014 . All these publications treat boundary conforming discretizations. Extensions to non-boundary conforming finite element discretizations are presented in Abedian2014 ; Taghipour2018 in the context of the finite cell method. Therein, it is demonstrated that high-order convergence rates may also be achieved even if the physical boundaries of the domain do not coincide with the discretization of the mesh. This is a desirable feature in the application at hand where the physical domain grows with time. Additionally, highly localized gradients need to be followed dynamically though the course of the simulation.
As a remedy for an accurate resolution of transient gradients on non-boundary conforming domains we advocate a combination of both, the multi-level -method and the finite cell method. To this end, we first introduce these numerical methods in section 2 where we begin by recalling the multi-level -method. It offers a relatively simple management of the degrees of freedom for transient discretizations of -type. Additionally, we present some novel but straight-forward extensions necessary for elasto-plastic computations on transient discretizations of -type. We then present the computation of elasto-plasticity in the framework of the finite cell method Duster2008 and proceed to a short description of the thermo-elasto-plastic problem in section 2.3 .
After the basics of the numerical methods are presented, we evaluate their performance against three benchmark examples. In section 3.1 a sphere under internal pressure in which a plastic front evolves is investigated. While the exact location of the plastic front is known in this specific case, in a more general setting it would not be unknown. We, thus, chose a grid-like discretization whose boundaries neither coincide with the evolving plastic front nor with the physical boundaries of the problem setup itself. We demonstrate that it is possible to capture the plastic front and the stress states on non-boundary conforming domains using higher order refinements and obtain higher order convergence rates. We then move to the thermo-elasto-plastic setting in section 3.2 to make the point that it is also possible to capture the plastic stress states in a coupled setting accurately. The final numerical example presented in section 3.3 is chosen to demonstrate that the presented methodology is also capable of reproducing stress states of real experiments. We conclude the article by pointing out the potential and limits of the presented approach in section 4 .
2 Numerical Methods
This section serves to introduce and further develop the numerical techniques which are used and evaluated in section 3 . To this end, the multi-level -method is first presented in section 2.1 followed by a description on the finite cell method section 2.2 . They are then combined to solve problems in elasto-plasticity as introduced in section 2.3 . All these sections form the background for solving thermo-elasto-plastic problems as introduced in section 2.4 .
2.1 The multi-level hp-method
The aforementioned multi-level hp-method is a powerful scheme for performing local mesh adaptation in an efficient manner. It was first introduced in Zander2015 and aims at a simple degree of freedom management for dynamic discretizations of high order by performing local refinements based on superposition. For this purpose, an initial discretization of coarse base elements is overlayed locally with multiple layers of finer overlay elements so as to better capture the solution behavior such as locally high gradients, see fig. 1 . This is in contrast to standard -methods which perform refinement by replacement where coarse elements are replaced by a set of smaller elements.
The methodology to overlay finite elements as in the multi-level -method requires the enforcement of linear independence of the basis functions and compatibility of the ansatz space at the boundary of the refinement zone. This is achieved in a straightforward manner by leveraging the direct association of topological components (nodes, edges, faces, volumes) with the degrees of freedom. Compatibility and linear independence are ensured through the deactivation of specific topological components. This deactivation is governed by simple rule-sets, which work for different spatial dimensions Zander2015 ; Zander2016 alike. This forms a key strength of the method as compared to classic -methods as the simplicity of the rule-sets allow for an easy treatment of dynamic meshes with arbitrary levels of hanging nodes. The multi-level -method results in -continuity over the complete refinement hierarchy in base elements while -continuity is obtained in the leaf elements — elements on the highest refinement level with no children.
2.2 The finite cell method
The ultimate in the scope of this article is to provide a framework which leaves as much geometric and topological freedom for the emerging additively manufactured artifact as possible. At the same time, the effort for mesh generation should be kept as low as possible. These are two main features of immersed methods, a class of advanced discretization techniques that significantly reduce the effort of mesh generation by utilizing a non-boundary conforming domain discretization. They have, therefore, emerged as the method of choice to perform numerical simulations on bodies with a complex shape, topology or a combination of both. The finite cell method (FCM) introduced in Parvizian2007 ; Duster2017 , is a prominent representative of immersed methods whose core idea is to combine the advantages of the fictitious domain approach with the computational efficiency of discretizations of high order. Its basic idea is depicted in fig. 2 . A body of complex shape and topology defined on a physical domain is extended by a fictitious domain . Their union yields a computational domain with a simple boundary which can be discretized using well shaped finite elements. These form the support of the ansatz functions and are termed finite cells as their boundary is not conforming to the boundary of the original, physical boundary . The physical domain must then be recovered at the level of the numerical integration of the associated bi-linear and linear forms. An indicator function
is introduced to classify points belonging to the physical or fictitious domain. Points within the physical domain are assigned a value, whereas points in the fictitious domain have a value .
2.3 Elasto-plasticity with the multi-level hp method
The finite cell method has already been successfully used in the field of plasticity. The work most relevant to the paper at hand is published in Abedian2013a in which the FCM was implemented for the flow theory with nonlinear, isotropic hardening for small displacements and small strains. Further, in Abedian2014 it is demonstrated that the FCM leads to more efficient discretizations than the standard, boundary conforming -version of the finite element method delivers. The FCM has recently also been extended to nearly incompressible finite strain plasticity with complex geometries, see Taghipour2018 . All these investigations were carried out on static discretizations in the sense that no dynamic refinement was applied locally to capture transient plastic fronts. However, in layered deposition modeling the size of the traveling heat source is comparatively small w.r.t. the rest of the computational domain. In these applications, a dynamic refinement and coarsening offers a discretizational advantage because the computational effort can thereby be balanced with the accuracy needed locally. However, this requires a dynamic management of primal and internal variables.
To facilitate a comprehensive but concise presentation of this subject, the next sections start with the classic formulation of plasticity in 2.3.1 . The standard setting is then cast into the FCM formalism in section 2.3.2 before the integration of the multi-level -method into the elasto-plastic FCM framework is treated in section 2.3.3 along with the associated management of primal and internal variables.
2.3.1 Classic plasticity
The classic weak form of equilibrium in a solid body is given by
where u is the displacement field, is the test function,
the Cauchy stress tensor,the strain tensor, b are the body forces acting on the domain, , and t is the prescribed traction on the Neumann boundary, . When a nonlinear material is utilized, the stress state is not only a function of the instantaneous strain. Instead, the stress state also depends on the history of the loads the body was subjected to. Consequently, the weak form (1) becomes nonlinear and is solved incrementally for each time step . Within each time step the internal variables, , which contain the history of the material, are assumed constant. Linearization of the weak form with respect to the unknown u around , which is the solution at iteration i, is given e.g. in de2009computational and reads
where the fourth order tensor is the tangent modulus defined as
2.3.2 The elasto-plastic finite cell method
In order to apply the finite cell method, the domain integrals in eq. 2 are multiplied with the indicator function such that
Since the solution in the fictitious domain is unphysical, computing the tangent modulus or stresses in the fictitious domain causes unnecessary computational overhead. Therefore, the deformation in the fictitious domain is neglected and stresses are assumed to be zero (). Moreover, the tangent modulus is taken as the elastic tangent (). Incorporating these assumptions into eq. 4 provides the following linearized weak form of equilibrium for elasto-plastic problems
In this paper, a rate-independent von Mises plasticity model with nonlinear isotropic hardening is considered whereby the displacements and strains are assumed to be small. The classic, additive decomposition of the strain tensor into an elastic and a plastic counterpart is then applied such that
Figure 3 summarizes this model, in which the internal variables are given by
where is the equivalent plastic strain.
2.3.3 Multi-level -adaptivity and the transfer of primary and internal variables
In the case of dynamic multi-level hp-refinements the primary variables given by the displacement field, u, and internal variables given by , need to be transferred from the old discretization prior to the refinement to the new discretization after the refinement was carried out. Since the displacement field over the domain is discretized by the basis functions, a global continuous description is available which is directly transferred to the new discretization by means of a global -projection. However, in the general elasto-plastic finite element procedure, the evolution of the internal variables is computed at local integration points via the plastic flow rule and the hardening law. As no continuous discretization is readily available, the transfer of these variables from the integration points of the old discretization to the integration points of the new one is more involved.
As a remedy, several strategies have been developed in literature. The simplest is the point-wise transfer of the internal variables Peric1996 . Therein, an area is associated to each old integration point within which constant values of the history variables are assumed. This leads to a discontinuous approximation of the history variables within the concerned finite elements. Another strategy is the element-wise transfer Ortiz1991
, where the internal variables are interpolated by local functions. This strategy can result in an approximation with discontinuities over element boundaries. A cure is offered by the large group of strategies involving nodal projections. Therein the internal variables are first transferred from integration points of the old discretization to the nodal degrees of freedom of that old discretization. This is achieved inPeric1996 and Lee1994 by extrapolating the values from integration points to the nodal dofs for each element locally and later averaging these extrapolated values. Another possibility is to use the super convergent patch recovery technique as presented in Zienkiewicz1992a ; Zienkiewicz1992b ; Boroomand1999 ; Khoei2009 , for the transfer of the values from the integration points of the old mesh to its nodes. In this approach a patch consists of elements that surround a node and the internal variables are fitted to continuous polynomials over each patch by a least squares method and are then interpolated to the nodal points. Independent of how the values are obtained at the nodes of the old discretization, the next step is to interpolate them to the new nodal points in the new mesh. Finally, the values of the internal variables at the new integration points are interpolated from the values of the new nodal points. Since this strategy acts as a smoothing operator, the localized internal variables are smeared out over a larger part of the domain. This can be a drawback in elasto-plastic analysis, because plasticity is a local effect in most cases.
Unfortunately, none of these techniques are suitable for the multi-level hp-adaptivity. This is mainly due to the fact that the high order basis functions are modal functions, which are not associated to specific nodes. For this reason a modified version of the element-wise transfer Ortiz1991 is advocated in the paper at hand in which each component of the internal variables, , is approximated such that
is a vector of integrated Legendre polynomials as used inp-version of FEM Duster2001 and is the corresponding coefficient vector. It is explicitly pointed out that these polynomials only span the leaf element, where denotes the local coordinates of this element. The coefficient vector, is determined by applying a least square fit to the discrete integration point values, , inside the leaf element. This is done by minimizing the function, ,
which is carried out by differentiating with respect to as
Equation (10) leads to the following linear system of equations
which needs to be solved for each component of the internal variables.
If this projection strategy is applied to the elements that contain the plastic front, i.e. in which only a group of integration points in the element accumulated plastic strains, then the proposed least squares approximation may lead to oscillations and spurious values. In some cases it is even possible to obtain invalid negative values for the equivalent plastic strain. Figure 4 illustrates this for the polynomial order for a simple one dimensional example in which the equivalent plastic strain should be nonzero for the two plastified integration points and zero for the three elastic integration points. This function cannot be represented by polynomials and the well known Gibbs phenomenon occurs.
In order to circumvent this problem, the values in those elements containing a elastic-plastic interface are not projected with the least squares method. Instead, for these elements we advocate the following approach depicted in Figure (a)a . Four cases are distinguished and characterized by the location of the integration points w.r.t. the concerned interpolation point. The first case is given in Figure (b)b where the interpolation point is located inside the gray bounding box defined by all integration points. In this case the internal variables, at are approximated by a trilinear interpolation from the values of the eight integration points closely surrounding it. In all the other cases the interpolation point is located outside the gray bounding box. Therefore, its projection onto the corresponding surface, edge or corner of the bounding box, , is used for interpolation. Figure (c)c depicts the second case in which the projection of is surrounded by only four points. These are then used for bilinear interpolation. The third case occurs when the projected point is surrounded by only two points as demonstrated in Figure (d)d . In this case linear interpolation is applied. The last case occurs when the original interpolation point lies in the corner quadrants of the finite cell. In this case values at the closest integration point are used. We will demonstrate in section 3 that this approach is feasible for dynamic multi-level -discretizations in an elasto-plastic setting.
2.4 Coupling to thermal problems
Metal deposition is a multi-physics problem, where the highest temperature gradients as well as the phase changes occur in close vicinity of the moving laser beam. The phase changes between liquid and solid states are simulated by the model introduced by Celentano et. al. Celentano1994 which uses the discretized weak form of equation (14). Further details regarding this method and its application to thermal analysis of the selective laser melting process in the framework of multi-level hp-adaptivity can be found in Kollmannsberger2016 .
A dynamically adaptable data container (multi-level grid) is used to keep track of the physical domain during the deposition process. The container represents a dynamic octree. It stores the current material state of a voxel at all points in time. Material states are stored in this container and not associated to finite elements. This decoupling of material and discretization facilitates the use of the finite cell method.
Due to the changes in the temperature during the process, some regions in the domain expand, while others contract. This generates residual stresses in the part when the thermal load is removed and the material allowed to cool down to its initial temperature. These stresses are computed by the quasi-static mechanical model, given in Section 2.3 with the addition of the thermal strain defined as
where is the thermal expansion coefficient and is the second order identity tensor. Adding the thermal strain component extends the additive decomposition of strains given in equation (6) to
It should be noted that, the von Mises plasticity depends solely on the deviatoric part of the strain tensor and the thermal strain is purely hydrostatic. Therefore, the basic structure of the elasto-plastic model introduced in Section 2.3 remains as is. Only the additive decomposition of the strains changes along with the fact that the yield stress and the thermal expansion coefficient are now temperature dependent functions.
It is assumed that displacements are small and do not produce heat, so that only a one-directional coupling has to be taken into account, i.e. only the displacement field is affected by the changes in the temperature field. As shown in Figure 10 , a staggered approach is taken for the solution of the thermomechanically coupled problem. For each time step, the multi-level grid is first updated according to the metal deposition such that the physical domain, , is defined for both the thermal and the mechanical problems. Then, the thermal problem is solved to obtain the temperature distribution. Finally, the resulting temperature field is used to compute the thermal strains and temperature-dependent material properties used to solve the mechanical problem before the next time step is computed.
3 Numerical Examples
All presented examples have benchmark character and are chosen thoroughly to test the main aspects of the proposed methodology. The first example discussed in section 3.1 is chosen to test if higher order convergence rates are possible in elasto-plastic computations and under what circumstances they decay. To this end, a dynamic multi-level -refinement as presented in section 2.1 is carried out towards the plastic front in several load steps whereby the plastic front itself is not resolved in a boundary conforming manner but travels through finite cells. This situation is expected for engineering applications. The example also tests the finite cell formulation presented in section 2.2 in the elasto-plastic embedded domain setting as described in section 2.3 by resolving the boundary of the physical domain only on integration level. It also serves as a test for the transfer of history variables as detailed in section 2.3.3 .
The second example given in section 3.2 is chosen to verify the thermo-elasto-plastic coupling procedure laid out in section 2.4 . Both fields are discretized on separate computational grids which are refined and de-refined individually according to the requirements of the corresponding field variable.
The final example presented in section 3.3 tests the combination of the methodology on a experimental benchmark used in welding. As such, it constitutes a first step towards a verification of the proposed methodology for layer deposition modeling.
3.1 Internally pressurized spherical shell
In this example internal pressure, , is applied in increments to a spherical shell, which is composed of an elastic perfectly plastic material with a von Mises yield criterion. Figure (a)a illustrates the problem setup, where , and are inner radius, outer radius and the radius of the plastic front, respectively.
Hill Hill1998 provides an analytical solution to this problem. Equation (16) gives the relation between the location of the plastic front, , and the applied internal pressure, , wherein is the yield stress.
Due to symmetry, only an octant of the whole spherical shell is considered in the numerical model with the appropriate symmetry boundary conditions. The shell is embedded in a Cartesian grid of size elements, where finite cells located completely outside the spherical shell are removed from the computation. For the numerical simulations the spherical shell is selected to have an inner radius of mm and outer radius of mm. The Young’s modulus, , and Poisson’s ratio, , are set to GPa and , respectively, and the yield stress, , is set to MPa.
The geometry of the spherical shell is captured at the integration level by utilizing the smart octree algorithm as explained in Kudela2016a . Figure (b)b illustrates the resulting integration points and cells for base elements.
Figure 17 shows the distribution of the radial stresses on a dynamically adaptive multi-level -discretization for three load steps, where the polynomial degree of the basis functions is chosen to be and maximum refinement level is set to . In each step the elements that are close to the plastic front are refined and the elements that are further away from the plastic front are coarsened. Moreover, in each step the primary variables related to degrees of freedom and internal variables at integration points are transferred as explained in Section 2.3.3 .
It can be seen from the equations (17) and (18) that the stresses vary only in -direction and are constant in - and -directions. Therefore, it is sufficient to compare the computed stresses along the radius of the spherical shell to the analytical stresses as depicted in Figure 21 and Figure 25 . Both radial and tangential stress results obtained from the numerical simulation match their analytical counterparts for each load step. The different material behavior in elastic and plastic regions introduces a kink in the stress field, which is more blatant for the tangential stress. The introduced multi-level -refinement towards that elastic-plastic interface allows the numerical solution to closely capture this behavior.
In order to investigate the convergence properties of this elasto-plastic problem in the framework of the multi-level -adaptive finite cell method as described in sections 2.3, 2.2 and 2.1 , p-convergence studies are performed by uniformly elevating the order p of the polynomial shape functions, while keeping the size of the base elements fixed. For these studies a base mesh is used, which is recursively refined towards the elastic-plastic interface. Moreover, an internal pressure MPa is applied in one load step for each computation. To this end, the relative error in energy is monitored, which is defined as
with the exact internal energy, , and the numerical internal energy, , computed as
By using the analytical solutions in equation (20) the exact internal energy for an octant of the spherical shell is obtained as follows
Figure 26 shows the relative error in energy with respect to the number of degrees of freedom for studies where the maximum refinement level is increased from zero to three. To rule out any domain integration errors, all volumetric integrals were evaluated exactly following Kudela2016a . The curve without any multi-level refinement depicts exponential convergence in the pre-asymptotic range until the error is dominated by the kink in the solution field caused by the elastic-plastic interface. This kink is neither resolved by the boundaries of the finite cells nor is the error controlled by a refinement in the vicinity of this irregularity in the solution. Therefore, in the asymptotic range, the error decreases only algebraically when the polynomial degree is further increased. The kink in the solution field is better approximated by the basis functions as the elements are refined towards the elastic-plastic interface with hierarchical multi-level refinement. This can be seen in curves with one and two refinement levels towards the elastic-plastic interface, where the point at which the convergence levels off to its asymptotic value occurs much later. The curve with three refinement levels converges exponentially. Here, the error at the interface does not dominate the overall error.
This demonstrates that the combination of the techniques presented in sections 2.3, 2.2 and 2.1 is able to maintain the expected convergence behavior for higher-order methods even if neither the physical boundaries of an artifact nor an evolving elasto-plastic front are explicitly captured by the boundaries of the discretization.
3.2 Thermo-elasto-plastic bar example
In this section, the thermo-mechanical model introduced in Section 2.4 is investigated in the framework of multi-level hp-adaptivity. To this end, an idealized version of the metal casting process is simulated on a semi-infinite bar. Initially, the bar, which is given in Figure 27 , is liquid with a uniform temperature, . Then, the temperature at the boundary is instantaneously changed to and kept constant. is lower than the melting temperature, . Therefore, the bar starts to solidify from this boundary onwards.
The analytical solution to the thermal part of this problem is described in Weber1912a and is known in the literature as Neumanns’s method, see e.g. Hu1996 ; Hahn2012 . The position of the liquid-solid interface is given as
where is the diffusivity of the material and is the time. The constant is computed by solving the following nonlinear equation
where is the latent heat and is the heat capacity. The analytical temperature distribution is given in equation (24) for the semi-infinite bar.
By using the temperature distribution over the body obtained from the Neumanns’s method, Weiner and Boley Weiner1963 developed an analytical solution for the thermal stresses for this problem. They assumed that the body is composed of an elastic perfectly plastic material with a yield stress, , which linearly decreases as the temperature reaches the melting point, where it becomes zero:
Moreover, the following dimensionless quantities are introduced,
where , and are Young’s modulus, Poisson’s ratio and thermal expansion coefficient, respectively. The analytical solution of the dimensionless normal stress components in y- and z-directions are then given as
where the coordinates and denote the elastic-plastic interfaces. According to this solution, the solid material behaves elastic within the range and plastic in ranges and . The positions of these interfaces are obtained by solving the following nonlinear equations
The values that are selected for the material properties of this problem are provided in Figure 27 . With these values, the constant is computed as 0.330825295611989, while the coordinates and are found to be 0.45487188 and 0.21570439, respectively.
The transient thermal problem is solved on a base mesh with 32 hexahedral elements of order that are distributed in x-direction. Homogeneous Neumann boundary conditions are applied along y- and z-directions and Dirichlet boundary conditions are applied in the yz-plane at and . In order to better represent the kink in the temperature field, the elements are dynamically refined three times towards the solid-liquid interface. The time step for the backward Euler scheme is chosen to be seconds.
The simulation for the quasi-static mechanical problem is carried out on a base mesh with 16 hexahedral elements of order which are dynamically refined twice towards both the solid-liquid interface and the elastic-plastic interfaces. During the simulation the liquid part of the bar is treated as a fictitious domain as explained in Section 2.2 by multiplying with . Extended plane strain conditions are applied along y- and z-directions such that and are constant. This is achieved by constraining all degree of freedoms corresponding to displacements in y- and z-directions to be equal via Lagrange multipliers. As for the x-direction, the bar is fixed at and left free at .
Figure 31 presents the results of the thermal and the mechanical problems for the time states seconds. For the thermal problem the evolution of the temperature field is provided over the dynamically adapted thermal mesh. As for the mechanical problem the stress component in z-direction () is given over the dynamically adapted mechanical mesh. The sequence of figures demonstrates how elements are refined and coarsened during the simulation for thermal and mechanical problems independently.
Figure 35 compares the temperature distribution along the bar to the analytical solution obtained by the Neumanns’s method. The kink in the temperature field at melting temperature, which is caused by the latent heat release during the phase change from liquid to solid, is well captured in all time states as well as the general temperature profile along the bar. The stress component in z-direction () along the bar is depicted in Figure 39 along with the analytical solution provided by Weiner and Boley. The kinks at the elastio-plastic interfaces and the solid-liquid interface are well represented due to refinements. Moreover, it can be seen that the numerical results closely match their analytical counterparts along the bar at all time steps of the computation.
3.3 Applications to metal deposition
In this section the performance of the numerical method introduced in Section 2 in simulating metal deposition processes is investigated against the benchmark problem of a single bead laid down on the top surface of a plate as published in Truman2009 . This benchmark problem is produced by the European Network on Neutron Techniques Standardization for Structural Integrity (NeT), whose mission is to develop experimental and numerical techniques and standards for the reliable characterisation of residual stresses in structural welds. Partners from industry and academia participated in this benchmark by constructing an experiment and measuring the residual stresses on the plate Ohms2009 ; Ficquet2009 ; Hofmann2009 ; Pratihar2009 ; Turski2009 or predicting the residual stresses through numerical methods Shan2009 ; Ohms2009 ; Gilles2009 ; Bate2009 ; Ficquet2009 ; Bouchard2009 . Smith et. al. compares the residual stress measurements and predictions of all these groups in Smith2009 .
Figure 40 illustrates the setup of the benchmark problem where the plate has dimensions 180 mm 120 mm 17 mm in length, width and thickness, respectively. Additionally, the bead that is laid down on the top surface of the plate has a width of 5 mm, thickness of 0.8 mm and a length of 60 mm.
According to the benchmark protocol, the plate and the bead material is AISI Type 316L stainless steel. NeT recommends using slightly different temperature dependent material properties for the parent and weld metal. However, in this work the weld metal is not distinguished from the parent metal. Table 1 provides the temperature dependent specific heat, conductivity, thermal expansion coefficient and Young’s modulus. Over the melting temperature of 1400 °C these values are assumed to be constant. In addition to that the density and Poisson’s ratio are assumed to be constant values of 7966 kg m and 0.294, respectively and latent heat is taken as 260 kJ/kg. Finally, the monotonic tensile test results, which show the true stress corresponding to the percentage of plastic strain, is given in Table 2 by NeT. The isotropic hardening behaviour of the material is based on this table.
The heat input from the welding torch, which travels with the speed of 2.27 mms, is defined by NeT as 633 Jmm. The moving heat source is based on the Goldak’s Goldaknewfiniteelement1984 ellipsoidal model, which is widely used in welding simulations. In this model the body load is given as
where , and define the position of the welding torch and , and are the semi-axes of the ellipsoid. Shan et. al.Shan2009 uses the weld bead profile measurements to decide the semi-axes of the ellipsoid, which is followed in this work. According to these measurements , and are assumed to be 1.9 mm, 3.2 mm and 2.8 mm, respectively. Moreover, in equation (29) is the power of the heat source, which is computed from the heat input and the welding speed as 1437 Watt.
The heat loss from the system is modeled with radiation and convection boundary conditions on the surface of the weld bead and the plate. During the process, the surface is updated as new weld bead is deposited. The convection coefficient and the emissivity are taken as 10 WmK and 0.75, respectively and they are assumed to be independent of the temperature. Additionally, the ambient temperature is also assumed to be constant at 20 °C.
Due to symmetry, only half of the plate and the weld bead is considered as shown in Figure 41 for the numerical model with the appropriate symmetry boundary conditions. The plate and the weld bead are embedded in a Cartesian grid of size elements with polynomial order 3 for both the thermal and mechanical problems. The dimensions of the computational domain are 180 mm, 60 mm and 18.6 mm in length, width and thickness, respectively. It should be noted that the thickness of the computational domain is greater than the total thickness of the plate and the weld bead, which is 17.8 mm. The first six rows of elements in thickness are used to discretize the plate, while the weld bead is embedded in the last row of elements.
Figure 44 shows a snapshot of the dynamic meshes of the thermal and the mechanical problems with material states, where red, blue and white (hollow) denote weld bead, plate and air, respectively. It should be noted that only the mechanical mesh is refined once along the welding line before deposition starts, which is not coarsened afterwards. Whereas, both meshes are dynamically refined twice towards the weld torch when material is deposited. The weld torch is set to travel 1.5 mm in each time step. Therefore, the newly deposited domain near the weld torch can be captured by the sub-element boundaries and the multi-level grid that keeps track of the physical domain. As the weld torch moves further away, these sub-elements are coarsened while the multi-level grid is used to remember the deposited domain as explained in Section 2.4 . Figure 48 depicts the dynamic meshes with temperature and longitudinal stress () distributions during the deposition period.
Residual stress distributions after the parts are cooled down to the room temperature are given in Figure 51 over the physical domain. Gilles et. al. Gilles2009 measures the longitudinal and transverse residual stresses with a neutron diffraction technique over the lines AB and CD. These lines can be seen in Figure 41 . Figure 54 and Figure 57 compare the numerically predicted residual stresses to the measurements from experiment by Gilles et. al. along the lines AB and CD, respectively. It can be seen from these figures that the measured and the numerically computed stress profiles are in good agreement. However, utilizing the complete isotropic hardening values in Table 2 results in higher residual stress especially for the longitudinal stress component, which is also reported by Gilles et. al. Gilles2009 . Their explanation for effect is that the isotropic hardening law which does not consider stress relaxation due to viscoplastic effects or hardening recovery is not good enough to simulate cyclic loads which occur during the welding process. Therefore, the isotropic hardening model is modified in this work such that over one percent of plastic strain the material is assumed to behave elastic perfectly plastic. The results with this modified model match the experimental values as illustrated in Figure 54.
4 Summary and conclusions
The article at hand evaluates if high-order finite elements may be used in metal deposition modeling. To this end, the finite cell method, an embedded domain method of high order, was used. It was demonstrated that this method, combined with the multi-level -method, leads to high-order convergence rates even if neither, the physical boundaries nor the plastic front is resolved by the boundaries of the computational mesh. Therefore, high-order embedded domain modeling is a valid computational option. Very accurate results were also achieved in the thermo-elasto-plastic setting. Therein, the thermal and the elasto-plastic field variables were computed on their own high-order discetizations and coupled by means of a classical staggered scheme. In all of the computations the state variables were decoupled from the finite cells in an FCM sense by managing them on another computational grid. This allows for their separate coarsening and refinement and a simulation where material can be submerged into the computation on a sub-element level i.e. the addition of material is not equal to the addition of finite cells. The methods were thoroughly verified in semi-analytical benchmarks and then applied to an experimental benchmark of metal deposition. A good agreement with the experimentally measured residual strains/stresses was found.
While the overall computational methodology leads to very accurate results open research questions remain. For example, the presented de-refinements in the mechanical part of the simulation are clearly a valid option to accurately compute the stress states in the vicinity of a plastic front. However, the stress state in a manufactured part does not only consist of one plastic front and rather smooth stress states to the side of it. Instead, very complex, local stress states may be present. A straightforward coarsening procedure of the grid in these regions will cause smoothing effects for variables which are not diffusive by nature. It is beyond the scope of this paper to study the influence of strong, local coarsening in complex stress states and their effect in more involved computations. The investigation of this challenge remains open. Since the local plastic strains are not connected to the finite cells but stored in an extra grid which acts like a computational database, one option would be not to coarsen them at all. The computational grid, carrying the degrees of freedom could be coarsened independently. Such a procedure would allow for a constant number of degrees of freedom throughout the computation while detailed information on the local plastic strains would still be available if needed even if the finite cells were coarsened.
Next, it is worthwile to investigate the wall-clock time advantages of the computational approach w.r.t. conventional modeling techniques in detail. This is not easy to some extent because a fair comparison requires an efficient implementation of both techniques which is not available at present. Clearly, the modeling flexibility gained by the presented embedded domain method of high order using transient refinements does introduce a certain implementational as well as a computational overhead as compared to a classical fixed-grid approaches. The break-even point of refining and de-refining regarding wall-clock time further depends on the scale of the computations. However, to refine everywhere to the finest level required locally is also not an option. Thus, refinements and de-refinements must be used for efficient computations. How this is possible was investigated in the article at hand, but its comparison to other refinement techniques remains an open question. The current computations hint that the type of microscopic computations in the metal depositing process presented in the last example pay off in a multi-layer process whereby the refinement and de-refinement is not carried out at each time increment.
In summary, the article presents the first verified and validated steps of a novel computational methodology for the analysis of metal deposition modeling in an embedded domain sense using locally refined, transient discretizations of -type. As such, it demonstrates that this computational methodology is a valid option for the analysis of metal deposition.
The authors gratefully acknowledge the financial support of the German Research Foundation (DFG) under Grant RA 624/27-2.
- (1) J. A. Goldak, M. Akhlaghi, Computational Welding Mechanics, Springer, New York, 2005.
- (2) L. E. Lindgren, L.E. Lindgren-Computational Welding Mechanics Thermomechanical and Microstructural Simulations, Woodhead Publishing Limited and CRC Press LLC, 2007.
- (3) M. Tanaka, An introduction to physical phenomena in arc welding processes, Welding International 18 (11) (2004) 845–851. doi:10.1533/wint.2004.3342.
- (4) O. Zinovieva, A. Zinoviev, V. Ploshikhin, Three-dimensional modeling of the microstructure evolution during metal additive manufacturing, Computational Materials Science 141 (2018) 207–220. doi:10.1016/j.commatsci.2017.09.018.
- (5) A. Basak, S. Das, Epitaxy and Microstructure Evolution in Metal Additive Manufacturing, Annual Review of Materials Research 46 (1) (2016) 125–149. doi:10.1146/annurev-matsci-070115-031728.
- (6) C. Körner, Additive manufacturing of metallic components by selective electron beam melting — a review, International Materials Reviews 61 (5) (2016) 361–377. doi:10.1080/09506608.2016.1176289.
- (7) L. E. Lindgren, Numerical modelling of welding, Computer Methods in Applied Mechanics and Engineering 195 (48–49) (2006) 6710–6736, 00127. doi:10.1016/j.cma.2005.08.018.
- (8) L. E. Lindgren, H. A. Häggblad, J. M. J. McDill, A. S. Oddy, Automatic remeshing for three-dimensional finite element simulation of welding, Computer Methods in Applied Mechanics and Engineering 147 (3–4) (1997) 401–409, 00052. doi:10.1016/S0045-7825(97)00025-X.
- (9) B. Schoinochoritis, D. Chantzis, K. Salonitis, Simulation of metallic powder bed additive manufacturing processes with the finite element method: A critical review, Proceedings of the Institution of Mechanical Engineers, Part B: Journal of Engineering Manufacture 231 (1) (2017) 96–117. doi:10.1177/0954405414567522.
- (10) V. Nübel, A. Düster, E. Rank, An rp-adaptive finite element method for the deformation theory of plasticity, Computational Mechanics 39 (5) (2007) 557–574. doi:10.1007/s00466-006-0111-4.
- (11) M. Olesky, W. Cecot, Application of the fully automatic hp-adaptive FEM to elastic-plastic problems, Computer Methods in Materials Science 15 (1) (2015) 204–2012.
- (12) T. Elguedj, T. Hughes, Isogeometric analysis of nearly incompressible large strain plasticity, Computer Methods in Applied Mechanics and Engineering 268 (2014) 388–416. doi:10.1016/j.cma.2013.09.024.
- (13) A. Abedian, J. Parvizian, A. Düster, E. Rank, Finite cell method compared to h-version finite element method for elasto-plastic problems, Applied Mathematics and Mechanics 35 (10) (2014) 1239–1248, 00002. doi:10.1007/s10483-014-1861-9.
- (14) A. Taghipour, Parvizian, J., S. Heinze, A. Düster, The finite cell method for nearly incompressible finite strain plasticity problems with complex geometries, Computers & Mathematics with Applications (2018) accepted.
- (15) A. Düster, J. Parvizian, Z. Yang, E. Rank, The finite cell method for three-dimensional problems of solid mechanics, Computer Methods in Applied Mechanics and Engineering 197 (45–48) (2008) 3768–3782. doi:10.1016/j.cma.2008.02.036.
- (16) N. Zander, T. Bog, M. Elhaddad, F. Frischmann, S. Kollmannsberger, E. Rank, The multi-level hp-method for three-dimensional problems: Dynamically changing high-order mesh refinement with arbitrary hanging nodes, Computer Methods in Applied Mechanics and Engineering 310 (2016) 252–277. doi:10.1016/j.cma.2016.07.007.
- (17) N. Zander, T. Bog, S. Kollmannsberger, D. Schillinger, E. Rank, Multi-level hp-adaptivity: High-order mesh adaptivity without the difficulties of constraining hanging nodes, Computational Mechanics 55 (3) (2015) 499–517. doi:10.1007/s00466-014-1118-x.
- (18) J. Parvizian, A. Düster, E. Rank, Finite cell method, Computational Mechanics 41 (1) (2007) 121–133. doi:10.1007/s00466-007-0173-y.
- (19) A. Düster, E. Rank, B. A. Szabó, The p-version of the finite element method and finite cell methods, in: E. Stein, R. Borst, T. J. R. Hughes (Eds.), Encyclopedia of Computational Mechanics, Vol. 2, John Wiley & Sons, Chichester, West Sussex, 2017, pp. 1–35.
- (20) A. Abedian, J. Parvizian, A. Düster, E. Rank, The finite cell method for the J2 flow theory of plasticity, Finite Elements in Analysis and Design 69 (2013) 37–47, 00000. doi:10.1016/j.finel.2013.01.006.
- (21) E. de Souza Neto, P. Perić, P. Owen, Computational Methods for Plasticity: Theory and Applications, John Wiley & Sons, 2009.
- (22) D. Perić, C. Hochard, M. Dutko, D. R. J. Owen, Transfer operators for evolving meshes in small strain elasto-placticity, Computer Methods in Applied Mechanics and Engineering 137 (3) (1996) 331–344. doi:10.1016/S0045-7825(96)01070-5.
- (23) M. Ortiz, J. J. Quigley, Adaptive mesh refinement in strain localization problems, Computer Methods in Applied Mechanics and Engineering 90 (1) (1991) 781–804. doi:10.1016/0045-7825(91)90184-8.
- (24) N.-S. Lee, K.-J. Bathe, Error indicators and adaptive remeshing in large deformation finite element analysis, Finite Elements in Analysis and Design 16 (2) (1994) 99–139. doi:10.1016/0168-874X(94)90044-2.
O. C. Zienkiewicz, J. Z. Zhu, The superconvergent patch recovery and a posteriori error estimates. Part 2: Error estimates and adaptivity, International Journal for Numerical Methods in Engineering 33 (7) (1992) 1365–1382.doi:10.1002/nme.1620330703.
- (26) O. Zienkiewicz, J. Zhu, The superconvergent patch recovery (SPR) and adaptive finite element refinement, Computer Methods in Applied Mechanics and Engineering 101 (1-3) (1992) 207–224. doi:10.1016/0045-7825(92)90023-D.
- (27) B. Boroomand, O. C. Zienkiewicz, Recovery procedures in error estimation and adaptivity. Part II: Adaptivity in nonlinear problems of elasto-plasticity behaviour, Computer Methods in Applied Mechanics and Engineering 176 (1) (1999) 127–146. doi:10.1016/S0045-7825(98)00333-8.
- (28) A. R. Khoei, S. A. Gharehbaghi, Three-dimensional data transfer operators in large plasticity deformations using modified-SPR technique, Applied Mathematical Modelling 33 (7) (2009) 3269–3285. doi:10.1016/j.apm.2008.10.033.
- (29) A. Düster, H. Bröker, E. Rank, The p-version of the finite element method for three-dimensional curved thin walled structures, International Journal for Numerical Methods in Engineering 52 (7) (2001) 673–703.
- (30) D. Celentano, E. Oñate, S. Oller, A temperature-based formulation for finite element analysis of generalized phase-change problems, International Journal for Numerical Methods in Engineering 37 (20) (1994) 3441–3465.
- (31) S. Kollmannsberger, A. I. Özcan, E. Rank, Computational modeling of the SLM process: A hierarchical treatment of the transient heat equation with phase changes in preparation.
- (32) R. Hill, The Mathematical Theory of Plasticity, Oxford Classic Texts in the Physical Sciences, Oxford University Press, Oxford, New York, 1998.
- (33) L. Kudela, N. Zander, S. Kollmannsberger, E. Rank, Smart octrees: Accurately integrating discontinuous functions in 3D, Computer Methods in Applied Mechanics and Engineering 306 (Supplement C) (2016) 406–426, 00021. doi:10.1016/j.cma.2016.04.006.
- (34) H. Weber, B. Riemann, Die Partiellen Differential-Gleichungen Der Mathematischen Physik: Nach Riemann’s Vorlesungen in 5. Aufl. Bearb, no. v. 1, F. Vieweg und Sohn, 1912.
- (35) H. Hu, S. A. Argyropoulos, Mathematical modelling of solidification and melting: A review, Modelling and Simulation in Materials Science and Engineering 4 (4) (1996) 371.
- (36) D. W. Hahn, Heat Conduction, 3rd Edition, Wiley, Hoboken, N.J, 2012.
- (37) J. H. Weiner, B. A. Boley, Elasto-plastic thermal stresses in a solidifying body, Journal of the Mechanics and Physics of Solids 11 (3) (1963) 145–154, 00124.
- (38) C. E. Truman, M. C. Smith, The NeT residual stress measurement and modelling round robin on a single weld bead-on-plate specimen, International Journal of Pressure Vessels and Piping 86 (1) (2009) 1–2. doi:10.1016/j.ijpvp.2008.11.018.
- (39) C. Ohms, R. C. Wimpory, D. E. Katsareas, A. G. Youtsos, NET TG1: Residual stress assessment by neutron diffraction and finite element modeling on a single bead weld on a steel plate, International Journal of Pressure Vessels and Piping 86 (1) (2009) 63–72. doi:10.1016/j.ijpvp.2008.11.009.
- (40) X. Ficquet, D. J. Smith, C. E. Truman, E. J. Kingston, R. J. Dennis, Measurement and prediction of residual stress in a bead-on-plate weld benchmark specimen, International Journal of Pressure Vessels and Piping 86 (1) (2009) 20–30. doi:10.1016/j.ijpvp.2008.11.008.
- (41) M. Hofmann, R. C. Wimpory, NET TG1: Residual stress analysis on a single bead weld on a steel plate using neutron diffraction at the new engineering instrument ‘STRESS-SPEC’, International Journal of Pressure Vessels and Piping 86 (1) (2009) 122–125. doi:10.1016/j.ijpvp.2008.11.007.
- (42) S. Pratihar, M. Turski, L. Edwards, P. J. Bouchard, Neutron diffraction residual stress measurements in a 316L stainless steel bead-on-plate weld specimen, International Journal of Pressure Vessels and Piping 86 (1) (2009) 13–19. doi:10.1016/j.ijpvp.2008.11.010.
- (43) M. Turski, L. Edwards, Residual stress measurement of a 316l stainless steel bead-on-plate specimen utilising the contour method, International Journal of Pressure Vessels and Piping 86 (1) (2009) 126–131. doi:10.1016/j.ijpvp.2008.11.020.
- (44) X. Shan, C. M. Davies, T. Wangsdan, N. P. O’Dowd, K. M. Nikbin, Thermo-mechanical modelling of a single-bead-on-plate weld using the finite element method, International Journal of Pressure Vessels and Piping 86 (1) (2009) 110–121. doi:10.1016/j.ijpvp.2008.11.005.
- (45) P. Gilles, W. El-Ahmar, J.-F. Jullien, Robustness analyses of numerical simulation of fusion welding NeT-TG1 application: “Single weld-bead-on-plate”, International Journal of Pressure Vessels and Piping 86 (1) (2009) 3–12. doi:10.1016/j.ijpvp.2008.11.012.
- (46) S. K. Bate, R. Charles, A. Warren, Finite element analysis of a single bead-on-plate specimen using SYSWELD, International Journal of Pressure Vessels and Piping 86 (1) (2009) 73–78. doi:10.1016/j.ijpvp.2008.11.006.
- (47) P. J. Bouchard, The NeT bead-on-plate benchmark for weld residual stress simulation, International Journal of Pressure Vessels and Piping 86 (1) (2009) 31–42. doi:10.1016/j.ijpvp.2008.11.019.
- (48) M. C. Smith, A. C. Smith, NeT bead-on-plate round robin: Comparison of transient thermal predictions and measurements, International Journal of Pressure Vessels and Piping 86 (1) (2009) 96–109. doi:10.1016/j.ijpvp.2008.11.016.
- (49) J. Goldak, A. Chakravarti, M. Bibby, A new finite element model for welding heat sources, Metallurgical Transactions B 15 (2) (1984) 299–305. doi:10.1007/BF02667333.