Powder bed fusion modelling based on Discontinuous Galerkin formulation

Online since 29 March 2021

Article

Abstract

Microscopic defects may occur during powder-bed fusion processes which originate from molten pool instabilities. To investigate this issue, a numerical model focusing on the interaction between laser and powder, melt pool formation and thermodynamics effects is on development. The approach used in the model assumes an incompressible Newtonian flow based on an enthalpy-porosity method for solid/liquid metal transition. The jump conditions at the liquid-solid interface for the energy equation are included through the modification of the enthalpy which incorporates the latent heat of fusion. A Carman-Kozeny porosity method is implemented in the momentum equation to penalize the flow in the solid phase. Molten flows are driven by natural and thermocapillary convection which are modelled using Boussinesq approximation and Marangoni free-surface stresses respectively. This numerical model is implemented within a Discontinuous Galerkin (DG) finite element formulation which ensures high order convergence on unstructured mesh [1].

In this work, a sharp interface method is integrated into the three dimensional high-order DG code Argo [2] to capture the free surface of the melt pool. This method allows to keep the high-order of convergence of the DG scheme even near the interface not conforming with the mesh. This is essential to capture the thermo-hydrodynamics phenomena during powder-bed fusion with accuracy.

Current effort is dedicated to improve the model accounting for physical phenomena such as surface tension and vaporization. The main purpose is to investigate the occurrence of discontinuities in fused powder tracks depending on material properties and process parameters.

Keywords

Powder Bed Metal Additive Manufacturing Process, Melt Pool, Sharp Interface Method, High Order Method, Discontinuous Galerkin Method

Table of contents

Text

1 Introduction

Powder Bed Metal Additive Manufacturing (PBMAM) technology is widely accepted and recognized as advanced manufacturing. The design and production of high-performance components is within reach for industry. Nonetheless, it is essential to have a better understanding of the powder-based feedstock materials, processes, structures, and properties to produce reliable PBMAM parts. In PBMAM processes, physical phenomena occurring during the layer-by-layer manufacturing, microstructure evolutions and mechanical properties of printed part are closely linked. Three different scales can be identified, (i) microscale, (ii) meso scale and (iii) macroscale. In this context, it is essential to predict melt pool dynamics involved during PBMAM process since microstructure evolution and microscopic defects depend mainly on the melting of powder layer and the solidification of the melt pool.

For process modelling, physics-based approaches are developed at meso scale. Meso-scale melt pool models can predict fluid flow and temperature distribution based on a multi-phase approach. The melt-pool model is computationally intensive in comparison with heat conduction model, the former is usually restricted to a layer deposition. Several discretization schemes to treat the multiphase flow can be found in the literature. In this work, a high order Discontinuous Galerkin scheme is used to provide high order convergence on unstructured mesh. This enables to capture with high accuracy those physical phenomena for such moving interface problem. The proposed CFD melt pool model is based on enthalpy-porosity formulation. This model is dedicated to hydrodynamics simulations under thermocapillary effects and thermal effects during solid liquid metal transition.

The numerical challenges of simulating multiphase flows are strongly related to the discontinuity in properties at the interface and jump conditions to be satisfied [3]. Multiphase problems imply a discontinuity in the solution (strong discontinuity) or in its derivative (weak discontinuity) on an interface evolving with time. To capture the discontinuity, two techniques are studied in this work the diffuse and the sharp interface methods. Diffuse methods smooth the discontinuities over a transition region across the interface. Sharp methods modify the stencil of the numerical scheme locally and can thereby keep the order of convergence of the method even in the vicinity of the interface.

The multiphase numerical treatment dedicated for each interface (solid-liquid and metal-gas interface) is described in the first section. And a numerical simulation of a static laser shot applied on Ti-6Al-4V resulted from the proposed model is illustrated in the last section.

2 Powder-bed fusion model

2.1 Solid/liquid interface

For most metallic alloys, the transition from solid to liquid state is applied over a range of temperature. A diffusive method is considered to smoothly represent the solid/liquid interface. In most diffuse numerical methods to treat liquid/solid interface, the mathematical model is modified to implicitly integrate the jump conditions in the constitutive law [4]. This approach uses a macroscale model where the microstructure in the mushy zone is not captured accurately. The solid and liquid phases are modeled as one single phase and the energy and momentum conservation laws are modified to distinguish locally the presence of liquid or solid phase. In order to integrate the jump condition for the energy equations, the classical enthalpy method is used [5]. This technique can describe the latent heat evolution by integrating the latent heat absorption during phase change in the mathematical model. For the momentum equation, a porous model is applied which consists of adding a source term to the momentum equation to penalize the velocity in the solid region.

The numerical model implemented in Argo uses a classical enthalpy-porosity constitutive law accounting for hydrodynamic and thermal effects occurring during the phase transition from solid to liquid metal. The system of equations considered is the incompressible Navier-Stokes equations augmented with phase-transformation thermodynamical effects, porosity-based mushy-zone modeling, Marangoni free-surface stresses and buoyancy-induced natural convection terms [1]. The system can be written as follows:

Image 10000201000003070000005857C7E3A423BC7EFA.png

with u, the velocity vector, ρ the density, P the pressure, τ the viscous tensor, S(u) the source term described below, T the temperature and λ the thermal conductivity. The enthalpy h includes the latent heat for phase transformation weighted by the liquid fraction α:

Image 100002010000030600000026DDA298588D20BF9A.png

where cp is the specific heat capacity, Lm the latent heat and α the liquid fraction which can be expressed as:

Image 100002010000030700000055BEC7E3F096029E87.png

In practice, the dependence of the liquid fraction on the melting temperature is smoothed to ensure differentiability of the model [1] using a hyperbolic tangent:

Image 10000201000003070000003E1E8E4AB4F5A1328E.png

The liquidus temperature Tl specifies the temperature above which a material is completely liquid while solidus temperature Ts specifies the value at which the material is considered solid. Fig. 1 illustrates the change of phase and the evolution of the volume fraction and enthalpy for a solid with rising temperature.

Image 1000000000000262000000C966B0BEC2631DEB34.png

Figure 1: Phase diagram and evolution of the volume fraction and enthalpy with the temperature

The source term S(u) in Eq. 2 includes two contributions for the natural convection and the porosity penalty term to model the solid phase. The Boussinesq approximation is used to model the natural convection is expressed as follows:

Image 100002010000030700000024F80D3FAEAAC0F6C5.png

where, β is the dilation coefficient and g represents the gravity vector. The Carman-Kozeny [6] penalty term, which ensures negligible velocity in the solid phase, can be written as:

Image 10000201000003070000005E71DD8AFFD45764FC.png

where the constant ϵ prevents a division by zero and the constant C is computed based on the interdentritic distance d shown in Fig. 2.

Image 100000000000006E0000008D8871D1F3F6E92E0F.png

Figure 2: Interdentritic distance to evaluate permeability of the mushy zone

2.2 Liquid/gas interface

While for most materials studied a numerical diffuse method can be used for the first interface type since the transition from solid to liquid is smooth, the gas-liquid is a sharp transition with large jump in properties. Therefore, sharp interface techniques are in theory more appropriate but are also more complex to be implemented.

Currently, a diffuse approach based on Continuum Surface Force method has been implemented within Argo to treat multiphase flow for additive manufacturing processes [1]. A level-set variable ϕ is used to capture the gas-liquid interface Γ:

Image 10000201000003070000001F34299D0B9B1AA1CC.png

The variable ϕ represents the signed distance to the interface; hence the zero iso-value is the position of the interface and the sign defines the phase. The scalar ϕ is updated with the local velocity field by solving this advective equation:

Image 10000201000003070000002C11C03631EDAA26DB.png

Within Argo, the system of equations (9) is highly coupled with the resolution of the levelset [7]. The level-set method is illustrated in Fig. 3.a.

Image 10000000000002D7000000BEAED628696036F392.png

Figure 3: (a) Level-set method to capture the free surface (liquid/gas interface), (b) Regularization of material properties and interface conditions

The level-set variable is then used to regularize the properties across the interface using a hyperbolic tangent plotted in Fig. 3.b and expressed as:

Image 100002010000030700000028783B4AFEF48FBE9D.png

For instance, the density across the interface can be computed as follows:

Image 10000201000003060000002B03175D09718C21C6.png

Interface conditions are enforced using a smoothed Dirac based on the level-set variable, those conditions are applied in a narrow band defined by a thickness ϵ around the zero iso-value as shown in Fig. 3.b:

Image 10000201000003070000002D342C29174C7A26FA.png

For instance, to account for capillary forces at the interface, a source term is added in the momentum equation presented in system (2):

Image 100002010000030700000024AD195C49813B5ACC.png

This term includes normal and tangential contributions of the surface tension noted σ. The normal force requires the evaluation of the surface curvature κ, and the tangential capillary force, which is referred to as Marangoni effect, includes the temperature influence on the surface tension:

Image 10000201000003070000002BD25CAF334E218C8E.png

Image 1000000000000101000000DAB99F835178162082.png

Figure 4: Illustration of the Extended Discontinuous method to capture discontinuities within cells

The major drawbacks of this diffuse method are i) the stiffness introduced when jump in the properties are large and ii) the high mesh dependencies to be able to represent this smooth transition. The smearing also introduces a degradation of the order of convergence in the vicinity of the interface, where phenomena driving the hydrodynamics occur. The idea of a sharp transition seems very promising to conserve the high-order convergence properties of Argo’s numerical method and to capture this kind of phenomena. A recent implementation to treat moving boundary problems sharply within Argo showed promising results [2] to simulate the liquid-gas interface. For this method, similarly to XFEM, the element cut by the interface are enriched to be able to represent a discontinuity (see Figure 4) within the cells. The nodes on these elements are duplicated, but the support of the solution is restricted by the sign of the level-set. This makes it possible to discretize a discontinuity (strong or weak) in the solution. The eXtended Discontinuous Galerkin method is described in several papers [7, 9] and Figure 5 illustrates the difference with respect to a diffuse method. The method has been implemented in Argo and verified on academic test cases for thermal problem mostly [2, 10]. The method has been proved to keep the high order convergence even in the vicinity of the interface if the quadrature on cut elements is chosen wisely. The spatial integration over cut cells is handled using the hierarchical moment fitting technique proposed in Muller et al. [9] and small cut cells are agglomerated to ensure the stability of the scheme [11].

Image 10000000000002E2000000BFB358E2076FB9E5EC.pngFigure 5: Comparison of diffusive and sharp interface methods for non-conforming interface mesh

The numerical study for additive manufacturing problems has shown the limits of diffuse methods when large jumps are considered and the challenges to impose accurately the interface conditions presented in [3]. The extension of the immersed interface method introduced in Argo seems promising and it is proposed to evaluate the challenges of extending the current method for hydrodynamic cases.

3 Numerical simulation

Current work is dedicated to improve and enrich the numerical model with additional physical phenomena such as vaporization. The diffuse and sharp interface capturing methods described before are both studied to highlight pros and cons of each approaches for PBMAM simulation.

Image 1000020100000307000000381381693D02713728.png

where βr is the retrodiffusion coefficient which represents the proportion of vapor returning to the melt pool, M is the molar mass of the alloy, R is the universal gas constant, p0 is the atmospheric pressure and Lv is the latent heat of vaporization. The ejection of metallic vapor implies a force towards the melt pool. For the proposed model, this force fv is proportional to the saturated vapor pressure generated at the liquid/gas interface:

Image 1000020100000307000000307CB26FDDB7EECBDE.png

For the diffusive interface method, Eq. 16 and 17 are applied on the interface by using the smoothed Dirac function described in Eq. 13.

For the sharp interface method, the heat sink is directly applied on the interface as surface heat sink. The Eq. 10 is modified by including a velocity of vaporization proportional uv to the vaporization rate v:

Image 100002010000030700000079D3BE474F8C736987.png

3.2 Static laser shot applied to Ti-6Al-4V

The 2D numerical test case consists in a titanium substrate exposed to a static gaussian laser flux qL expressed as follows:

Image 10000201000003070000002AC1F5EF0F912AB7BD.png

where R is the reflection coefficient, PL is the laser power and r0 is the laser radius at focal point. In the case investigated, the power is equal to 250W with the radius of 50μm. Natural convection and capillary effect are also considered. The surface tension is not considered for this configuration and to limit the deformation of melt pool the Marangoni coefficient is reduced by two orders of magnitude. The material properties of the substrate and the gas are described in [12].

A symmetry plan is placed at the center of the laser heat flux and aligned with the laser beam direction. The domain sizes mentioned in [12] are duplicated. Adiabatic thermal condition and stick wall condition are imposed on lateral faces and bottom zone.

3.3 Results and discussions

Numerical simulations of the static laser shot presented in this section have been produced with the sharp interface approach including the vaporization model described in Eq. 16 and 18 during a simulation time of 150μs. The interface evolution is tracked at each time step and is represented by the zero iso-value level-set. The melt pool characteristics and evolution are illustrated in the Figure 6. The melt pool is bounded by the top surface and the white curve which represents the liquidus iso-value. Recirculation are constantly observed at the top edge of the melt pool layer and the flow is progressively intensified during the simulation as shown in Fig. 7. The temperature inside the melt pool is restricted and the maximum is closed to the boiling point due to the vaporization heat sink. The vaporization front progression and the melt pool thickness at the center of the heat source are measured for each simulation time. The simulation results show that the vaporization front progression is constant (increment of +18.5μm at each time step of 50μs) and the thickness of melt pool is constant during the vaporization (13.3μm). When the vaporization occurs on the surface of the melt pool, the vaporization velocity remains uniform due to a constant surface temperature near the boiling point (energy balance at the top surface) and the atmospheric pressure is fixed in this configuration.

Image 10000201000004B50000049DB9EB332E32FDBAFD.png

Figure 6: Evolution of the melt pool including temperature and velocity distribution near the heat source, vaporization length and melt pool thickness at the center of the heat source for different times: (a) 50μs, (b) 100μs and (c) 150μs

4 Conclusion & perspectives

To simulate and predict melt pool evolution during PBMAM processing, a melt pool model based on enthalpy-porosity formulation is currently in development. The main challenge remains the numerical treatment of the liquid/gas interface, the immersed interface method is proposed to succeed to the diffuse interface method.

Currently, main physical phenomena involved in PBMAM process are investigated such as vaporization. First numerical simulations have been performed with a vaporization model and have shown the important role of the vaporization in the keyhole formation.

Future work will be dedicated to surface tension application and thermal characterization of metal powder material.

Bibliography

[1] J.-S. Cagnone, K. Hillewaert and N. Poletz. A Discontinuous Galerkin Method For Muliphysics Welding Simulations, Key Engineering Materials, Vol. 611-612, pp. 1319-1326, (2014).

[2] P. Schrooyen, L. Arbaoui, J.-S. Cagnone, N. Poletz and K. Hillewaert. A High-order Extended Discontinuous Galerkin Method to Treat Hydrodynamics Problems, ECCM-ECFD Conference, June 11-15 2018, Glasgow (UK).

[3] D. Henneaux and A. Turchi. Identification of thermo-physical models experimental test derivation. Contractor Report ESA Contract no.4000125437/18/NL/RA TN-1.2, von Karman Institute for Fluid Dynamics, March 2019.

[4] A. Faghri and Y. Zhang. Transport Phenomena in Multiphase Systems. Elsevier, 2006.

[5] P.A. Nikrityuk. Computational thermo-fluid dynamics: in materials science and engineering. John Wiley and Sons, 2011.

[6] F.A.L. Dullien. Porous media: fluid transport and pore structure. Academic Press, San Diego, 1991.

[7] F. Pochet, K. Hillewaert, P. Geuzaine, J.F. Remacle, and E. Marchandise. A 3D strongly coupled implicit discontinuous Galerkin level-set based method for modeling two-phase flows. Computers and Fluids, 87(25):144{155, 2013.

[8] F. Kummer and M. Oberlack. An extension of the discontinuous galerkin method for the singular poisson equation. SIAM Journal of Scientific Computing, 35(2):603{622, 2013.

[9] B. Muller, F. Kummer, and M. Oberlack. Highly accurate surface and volume integration on implicit domains by means of moment-fitting. International Journal for Numerical Methods in Engineering, 96(8):512-528, 2013.

[10] P. Schrooyen, J.S. Cagnone, N. Poletz, and K. Hillewaert. A high-order extended discontinuous Galerkin method for moving immersed interface problems. In eXtended Discretization MethodS, Umea, Sweden, 2016.

[11] F. Kummer. Extended discontinuous galerkin methods for multiphase flows: the spatial discretization. In Annual Research Briefs, Center for Turbulence Research, Stanford, USA, 2013.

[12] A. Queva, G. Guillemot, C. Moriconi, C. Metton, M. Bellet, Numerical study of the impact of vaporisation on melt pool dynamics in LaserPowder Bed Fusion - Application to IN718 and Ti-6Al-4V, Additive Manufacturing Volume 35, October 2020, 101249.

Illustrations

References

Electronic reference

Larbi Arbaoui, Pierre Schrooyen, Nicolas Poletz and Koen Hillewaert, « Powder bed fusion modelling based on Discontinuous Galerkin formulation  », ESAFORM 2021 [Online], Online since 29 March 2021, connection on 29 March 2024. URL : https://popups.uliege.be/esaform21/index.php?id=2441

Authors

Larbi Arbaoui

Cenaero, Belgium

University of Liege, Belgium

Corresponding author: larbi.arbaoui@cenaero.be

Pierre Schrooyen

Cenaero, Belgium

Nicolas Poletz

Cenaero, Belgium

Koen Hillewaert

Cenaero, Belgium

University of Liege, Belgium

Details

Title
Powder bed fusion modelling based on Discontinuous Galerkin formulation
Language
en
Author, co-author
Larbi Arbaoui, Pierre Schrooyen, Nicolas Poletz and Koen Hillewaert,
Publication date
14 April 2021
Journal title
ESAFORM 2021
Copyright
CC-BY
DOI
10.25518/esaform21.2441
Permanent URL
https://popups.uliege.be/esaform21/index.php?id=2441

Statistics

Views : 262 (14 ULiège)
Downloads : 0 (0 ULiège)

Cite

Arbaoui, L., Schrooyen, P., Poletz, N., & Hillewaert, K. (2021). Powder bed fusion modelling based on Discontinuous Galerkin formulation . Paper presented at ESAFORM 2021. 24th International Conference on Material Forming, Liège, Belgique. doi: 10.25518/esaform21.2441