Simulating fast electron beam melting with a parallel thermal free surface lattice Boltzmann method

Computers and Mathematics with Applications (

)

–

Contents lists available at ScienceDirect

Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa

Simulating fast electron beam melting with a parallel thermal free surface lattice Boltzmann method Regina Ammer a,∗ , Matthias Markl b , Ulric Ljungblad c , Carolin Körner b , Ulrich Rüde a a

Chair for System Simulation, Cauerstr. 11, 91058 Erlangen, Germany

b

Chair of Metals Science and Technology, Martensstr. 5, 91058 Erlangen, Germany

c

Arcam AB, Krokslätts Fabriker 27 A, 431 37 Mölndal, Sweden

article

info

Keywords: Thermal lattice Boltzmann method Electron beam melting process Free surface Parallel computing

abstract This paper introduces a three dimensional (3D) thermal lattice Boltzmann method for the simulation of electron beam melting processes. The multi-distribution approach incorporates a state-of-the-art volume of fluid free surface method to handle the complex interaction between gas, liquid, and solid phases. The paper provides a detailed explanation of the modeling of the electron beam gun properties, such as the absorption rate and the energy dissipation. Additionally, an algorithm for the construction of a realistic powder bed is discussed. Special emphasis is placed to a parallel, optimized implementation due to the high computational costs of 3D simulations. Finally, a thorough validation of the heat equation and the solid–liquid phase transition demonstrates the capability of the approach to considerably improve the electron beam melting process. © 2013 Elsevier Ltd. All rights reserved.

1. Introduction Electron beam melting (EBM) is an additive manufacturing method used to produce metallic structures layer by layer from metal powder (Fig. 1). The latter is positioned on the starting plate included in the vacuum chamber and is melted by the electron beam gun. After a short time of solidification the powder hopper pushes the next powder layer on top and the EBM process starts again. Numerous products are generated by EBM, e.g., medical implants like hip joints or artificial spinal disks or components for aerospace and automotive industry.1 The possibility to construct very complex structures which are strong and flexible is one of the most important advantages of EBM manufacturing. The target objectives of the 3D simulation of the EBM process are the acceleration of the building process and the improvement of the production accuracy. Therefore, a better understanding of the beam–powder interaction is necessary. Based on a two dimensional thermal free surface lattice Boltzmann model [1] we develop a thermal 3D free surface lattice Boltzmann (LB) implementation where powder particles are taken into account. The model includes hydrodynamic physical effects, e.g., the flow of the melt pool as well as thermal effects, e.g., beam absorption, melting and solidification. However, 3D simulations require a high computational intensity and thus a thorough parallelization and optimization of the program is required. The implementation is integrated in the waLBerla framework (widely applicable lattice Boltzmann solver from Erlangen) which is used for solving problems in Computational Fluid Dynamics. waLBerla is extended by several thermal effects necessary for the simulation of the EBM process, e.g., absorption and solidification. The generation of metal powder

∗

Corresponding author. E-mail addresses: [email protected], [email protected] (R. Ammer).

1 ArcamAB: http://www.arcam.com/. 0898-1221/$ – see front matter © 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.camwa.2013.10.001

2

R. Ammer et al. / Computers and Mathematics with Applications (

)

–

Fig. 1. EBM process.

particles, represented as rigid bodies, is supported by the physics engine pe framework. Special focus lies on the parallel implementation of all algorithms to take into account of high computational costs of 3D simulations. Carefully chosen numerical tests are used to evaluate the EBM method and gain insight in their fundamental processes. For physical validation, numerical results of special scenarios are compared with experimental results. Hence, the scenarios cover various electron beam properties, e.g., different energy distributions, beam powers and velocities of the beam movements which can be specified individually for each scenario. For the simulations of melting domains of meaningful size these scenarios are implemented in parallel and have to be computed on supercomputers. Our results show the high potential of the thermal LB approach to understand and develop complex processes like EBM depending on thermodynamic as well as on fluid dynamic considerations. The remainder of the paper is organized as follows. Section 2 explains the thermal free surface LB method where special attention is paid to cell conversions and solid–liquid phase transition. Subsequently, Section 3 focuses on the aspects concerning the modeling of the EBM process, i.e., the definition of the electron beam, the different absorption types and the generation of the powder particles. In Section 4 information about the two software frameworks waLBerla and pe is provided, followed by validation examples like the comparison of the analytical and numerical solution of the heat equation and the Stefan problem. Section 4 is closed with a numerical example of melting a spot in a powder bed where the behavior of the simulated melt pool is discussed. Finally, numerical aspects and targets of future work are described in Section 5. 2. Numerical methods The governing macroscopic equations to simulate the incompressible transport equations for the EBM process are given by

∇ · u = 0, 1 ∂u + (u · ∇) u = − ∇ p + ν∇ 2 u + g, ∂t ρ ∂E + ∇ · (uE ) = ∇ · (k∇ E ) + Φ , ∂t

(1) (2) (3)

where t denotes the time, u the velocity of the melt, p the pressure, ρ the density, ν the kinematic viscosity, k the thermal diffusivity, g the gravity and E the energy density. Moreover, Φ identifies the beam energy absorbed by the material. For the numerical discretization of Eqs. (1) and (2) we use an isothermal 3D LB method (LBM). A thermal LBM is explained for the numerical treatment of Eq. (3) based on a multi-distribution approach where the energy density is modeled as a passive scalar. We briefly discuss advantages of this multi-distribution approach. For the treatment of the moving boundary between the melted metal and the vacuum in the EBM chamber we use a volume-of-fluid free surface approach. The idea of the free surface approach is shown before we conclude Section 2 explaining the liquid–solid phase transition and the regarding rules for cell conversions. 2.1. The LBGK model for hydrodynamics The first lattice Boltzmann models were introduced as alternative numerical methods for hydrodynamic problems by [2] and can be seen as a discretization of the Boltzmann equation [3]. An improvement of LBM compared to its predecessor, lattice gas automata, is the existence of a particle distribution function (pdf). The pdf f (x, v, t ) indicates the probability of finding a particle with the macroscopic velocity v(x, t ) at position x. The basic idea of the LBM [4,5] is to solve the Boltzmann

R. Ammer et al. / Computers and Mathematics with Applications (

)

–

3

equation in the hydrodynamic limit for pdf in the physical momentum space

∂f ∂f + v · ∇f + F · =Ω (4) ∂t ∂v where Ω denotes the collision operator, v the microscopic velocity and F is an external force, e.g., gravity. The BGKapproximation of the collision operator [6] proposed by Bhatnagar, Gross and Krook, can be written as

 ∂f ∂f 1 + v · ∇f + F · = − f − f eq ∂t ∂v τ where f fluid by

eq

(5)

denotes the Maxwell equilibrium distribution [3]. τ is the relaxation time and related to the viscosity ν of the

ν = cs2 ∆t (τ − 0.5)

(6)

with the lattice speed of sound given by = 1/3∆t /∆x where ∆t and ∆x denote spatial and time resolution. We use the D3Q19 stencil for the discretization in the microscopic momentum space, i.e., a finite set of 19 discrete velocities {ei |i = 0, . . . , 18}, [3] and lattice weights {ωi |i = 0, . . . , 18} are defined by cs2

ωi =

 1     3  

i=0

1

 (0, 0, 0) ei = (±1, 0, 0), (0, ±1, 0), (0, 0, ±1) (±1, ±1, 0), (0, ±1, ±1), (±1, 0, ±1).

i = 1, . . . , 6

 18     1

(7)

i = 7, . . . , 18 36 Relating to the D3Q19 stencil a lattice cell is defined as a cubic volume of unit length ∆x centered around a lattice node. With this stencil Eq. (5) can be written as

 ∂ fi (x, t ) ∂ fi 1 + ei · ∇ fi (x, t ) + Fi = − fi (x, t ) − fieq (x, t ) ∂t ∂ ei τ where fi (x, t ) = f (x, ei , t ). The Maxwell distribution is discretized by a Taylor series expansion [5] and given by   (ei · u) (ei · u)2 u2 eq fi (x, t ) = ωi ρ 1 + + − . 2 4 2 cs

2cs

This discrete equilibrium function fi uT u

Ma :=

cs2

eq

2cs

(8)

(9)

is only valid in a small Mach number regime, i.e.,

≪ 1.

(10)

The discrete external force term Fi [7] in Eq. (8) is defined as



Fi = ωi ρ

(ei − u) cs2

+

(ei · u)ei



cs4

·g

(11)

where g denotes gravity. The macroscopic quantities density and momentum are computed by

ρ=



fi

i

and ρ u =



ei fi .

(12)

i

In this LBM the pressure is linearly related to the density which is the ideal gas law p = cs2 ρ.

(13)

This LBM can be divided into a stream and a collide step. The stream step propagates the local pdfs of each cell along the corresponding lattice link into the neighboring cells and the resting particle with zero velocity remains in the cell center. It is given by fi′ (x + ei ∆t , t + ∆t ) = fi (x, t ) .

(14) ′

In the collide-step, the relaxation towards the local equilibrium is performed for each cell for the pdfs fi by

 1  eq fi (x, t + ∆t ) = fi′ (x, t + ∆t ) + fi (ρ, u) − fi′ (x, t + ∆t ) − Fi . τ

(15)

We call Eqs. (14) and (15) LBGK-scheme. Chen [8] and Qian [9] show that the approximated BGK collision operator proposes enough freedom for the equilibrium distribution to satisfy isotropy, Galilean invariance and a velocity-independent pressure. Some more comments to the boundary conditions to conclude the discretization of the hydrodynamic part. We use a no-slip boundary condition with half-way bounce-back for the flow field. So far, this LBGK covers the continuity and momentum equations of Eqs. (1) and (2) without energy conservation and thermal effects. In order to solve Eq. (3) for simulating melting and solidification in the EBM process we use a thermal LBGK scheme which is introduced by Körner in [1] for the 2D model.

4

R. Ammer et al. / Computers and Mathematics with Applications (

)

–

2.2. The LBGK model for temperature The two main thermal LB approaches [10] are the multi-distribution method [11,12,10,13] or the multi-speed approach [14–16]. The thermal multi-speed LB approach can be seen as an extension of the isothermal LBM. It has additional discrete velocities and higher order velocity terms in the equilibrium functions. The internal temperature is regarded as an additional moment but it has been shown in the past that multi-speed approaches suffer from stability problems [17] and are limited to simulating only one Prandtl number which is the ratio between kinematic viscosity and thermal diffusivity. Using the multi-distribution approach, the thermal flow is modeled by a second particle distribution function set, i.e., the first set of pdfs fi is used for simulating density and momentum and the second one hi for the temperature or, like in our application, for the energy density E (compare Eq. (20)) which is defined as a passive diffusing scalar [11,18,19]. We use this multi-distribution ansatz because it can deal with almost arbitrary Prandtl numbers which is important for the EBM process. Regarding Eqs. (14) and (15) the thermal LBM pdf set hi can also be written as a stream-and-collide step, respectively

 h′i (x + ei ∆t , t + ∆t ) = h′i (x, t ) hi (x, t + ∆t ) = hi (x, t + ∆t ) + ′

1

τh

(heq i (ρ, u) − hi (x, t + ∆t )) + Φi ,

(16)

where also a discrete BGK-approximation for the collision operator and the D3Q19 stencil is used with the same velocity vectors and lattice weights given in Eq. (7). Φi is the energy source related to each lattice cell which can absorb energy from the beam and is defined as

Φi (x, t ) = ωi Eb (x, t ),

(17)

where Eb stands for the energy provided by the electron beam. More information on Eb can be found in Section 3.1. τh denotes the relaxation time belonging to the second set of pdfs hi and is related to thermal diffusivity k by k = cs2 ∆t (τh − 0.5).

(18)

The equilibrium distribution



eq

hi (x, t ) = ωi E 1 +

eq hi

is defined by

(ei · u) cs2



.

(19)

Here, a linear equilibrium function is sufficient to recover Eq. (3). Regarding Eq. (12) the third macroscopic quantity, the energy density, can be computed by E=



hi .

(20)

i

We use energy density instead of temperature since the temperature as a function of time is not continuously differentiable [20]. If the metal powder is heated by the electron beam the temperature increases constantly until the range between solidus Ts and liquidus Tl temperature is reached. If Ts is reached the temperature increases slower since the powder starts to melt and this process cools the system. If Tl is attained the temperature increases faster again. The function of the energy density increases constantly over time and due to numerical stability the energy density is computed as a macroscopic quantity by the LBGK scheme. Energy density E and temperature T are related by the following integral formula T



ρ cp dT + ρ ∆H ,

E=

(21)

0

where cp is the specific heat capacity at constant pressure and ∆H is the latent enthalpy of a lattice cell undergoing phase change. The metal powder in the EBM process consists of the multi-component superalloy TiAl6V4. For this alloy, ∆H is a function depending on the temperature

 L  T −T s ∆H (T ) = ·L  T − T s  l 0

T ≥ Tl Ts < T < Tl

(22)

T ≤ Ts ,

where L denotes the latent heat of the phase change. Regarding Eq. (21) the energy density can be computed by integrating over the temperature. In summary, the coupled-thermal LBGK is written as the following stream–collide algorithm

 streaming

fi′ (x + ei ∆t , t + ∆t ) = fi (x, t ) h′i (x + ei ∆t , t + ∆t ) = hi (x, t ),

 1 eq  fi (x, t + ∆t ) = fi′ (x, t + ∆t ) + (fi (ρ, u) − fi (x, t + ∆t )) + Fi τf collide  hi (x, t + ∆t ) = h′i (x, t + ∆t ) + 1 (heq (E , u) − hi (x, t + ∆t )) + Φi . τh i

(23)

(24)

R. Ammer et al. / Computers and Mathematics with Applications (

)

–

5

This thermal LBGK algorithm can be used for simulating the EBM process regarding hydrodynamic and thermal effects but the discretization does not consider a temperature dependent surface tension. There is only a one-way coupling between the two pdf sets taking place by hydrodynamic velocity u. The hydrodynamic velocity u is used to compute the linear thermal equilibrium function in Eq. (19). The surface tension is constant and independent from changes in the temperature. Some further comments to the thermal boundary conditions which are necessary for the discretization part. There exist a open thermal and a closed thermal boundary condition. By using the open thermal condition a fixed temperature can be set at the boundary and the gradient of temperature is zero at the boundary thus no heat conditions exists and the unknown pdf values are initialized by the equilibrium values of the local energy density values. Taking the closed thermal boundary condition we assume that no heat exchange with the environment is possible and thus, no heat condition in normal direction happens. This boundary condition is implemented by a mirror-reflection method with half-way bounce back. We use open thermal boundary condition for the validation examples in Sections 4.2 and 4.3 where walls are heated or cooled. Completing the discretization part, we show the treatment of the moving boundary between the vacuum and the liquid metal in the next section. 2.3. Free surface extension In literature there exist various approaches to handle moving boundaries numerically. Interfaces can be modeled by color-fluid methods [21,22] and phase-field methods, which have already been used for applications similar to EBM process, regarding solidification [23] and multi-component alloys, especially TiAl6V4-alloy [24]. Furthermore, level-set methods are introduced by [25,26] to capture interfaces and are optimized using a local approach [27]. The idea of the color-fluid method is the redistribution of a colored density at each node to separate the different phases, e.g., gas and liquid phase. This redistribution approach is expensive due to of local maxima computations and the colored redistributed pdfs can evoke anisotropic surface tension and oscillations. An improved color-fluid model is introduced by Latva-Kokko and Rothman [28] without spurious oscillations at the interface and anisotropic effects. This improved model enables only constant surface tension independent of the temperature. The phase-field method has already become widely applicable to describe phase transformations, especially for solidification of multi-component systems [29]. The interface is represented by a phase-field function that describes the phase movement. Level-set methods define the interface by a level set function which is Lipschitz continuous and has an opposite sign for the different phases. The value of the level set function at the interface is zero. The movement of the interface over time makes an update of the level set function necessary. Level-set methods for interface capturing are expensive and although there exist already faster, local approaches like those introduced in [27], it can be shown that the level-set methods are not strictly mass conserving [30]. The interface capturing method for the EBM process simulation is a volume-of-fluid free surface (FS) flow approach [31] which is explained in the following. It assumes that the liquid phase governs all necessary flow dynamics completely and consequently, the gas phase can be neglected. This assumption enables high density differences between the two phases. Especially in the EBM process high density differences exist between the gas phase in the vacuum chamber and the metal powder. All interface capturing methods described above have problems with such high density differences. In the FS approach a closed layer called the interface layer separates gas and liquid phases. The interface layer is exactly one cell. Each lattice cell has a cell type, e.g., wall, gas, liquid or liquid interface (Fig. 2). Regarding this interface layer the FS ansatz is as a special boundary condition. In relation to the volume-of-fluid idea a dynamic fill level ϕ(x, t ) is introduced. This fill level is defined as the fraction of the cell volume currently filled with mass of liquid m(x, t )

ϕ (x, t ) =

m(x, t ) , ρ (x, t ) (∆x)3

(25)

where ρ is the local density. The interface cells at the free surface can be filled partially, i.e., 0 < ϕinterface < 1 contrary to liquid cells which have to be completely filled and gas cells which have a fill level of zero. The mass exchange for these interface cells is computed directly at the stream step in Eq. (23) and corresponds to a change of the fill level. The mass exchange between two neighboring cells x and x + ei can be defined as

  0

f (x + ei , t ) − fi (x, t ) ∆mi = 1i     (ϕ(x, t ) + ϕ(x + ei , t )) · fi (x + ei , t ) − fi (x, t ) 2

if x + ei gas if x + ei liquid

(26)

if x + ei interface

where fi denotes the opposite direction of fi . Since gas phase flow dynamics are neglected, a special treatment for the interface cells is needed and only the gas pressure of the conducted gas region can be considered at the interface. The normal vector n(x) at the free surface at position x is approximated by the gradient of the fill levels within a local neighborhood of the cells. With this approximated normal vector n(x) the pdf fi with eTi n ≤ 0 (pointing towards the liquid phase) are set to eq

eq

fi′ (x, t + ∆t ) = fi (ρG (x), u(x)) + f (ρG (x), u(x)) − fi′ (x, t + ∆t ) i

(27)

introduced by Körner et al. in [31]. ρG is the gas pressure pG (see Eq. (13)), which has to be known, since it was shown in [31] that p = cs2 ρG = pG . u(x) is the current liquid velocity at the boundary. The effect of surface tension is also included in

6

R. Ammer et al. / Computers and Mathematics with Applications (

)

–

Fig. 2. Example of a deformed drop on a bottom wall boundary. Gas phase is separated from the fluid phase by an interface layer.

the FS boundary condition. If σ and κ(x, t ) are given as a constant surface tension parameter and as a local curvature at the interface respectively the gas pressure can be modified by the Laplace pressure to p′G = pG + 2σ κ.

(28)

At this part some comments regarding heat exchange between particles and atmosphere are in order. There is no heat exchange from the particles to the atmosphere by the free surface. The thermal boundary conditions are given by a reflection condition at the interface. The unknown pdfs hi coming from the gas phase are approximated by the outgoing hi from the liquid phase. This FS approach models the sharp interface between the metal powder particle and the vacuum appropriately. To guarantee a closed interface layer after mass exchanging operations we show several cell conversion rules and solid–liquid phase transitions in the next section. 2.4. Cell conversions & phase transition Before explaining different types of cell conversions we have to extend the set of possible cell types by thermal cell types since we have to differ between the liquid–gas interface and the solid–gas interface (Fig. 3, green, dotted arrow) due to the FS approach. Thus, each lattice cell in the simulation domain has one of the following types [31]: liquid, gas, solid, solid interface or liquid interface. We start with the liquid–gas cell conversions. Liquid and gas cells can change their state but an immediate change from gas to liquid and vice versa is not allowed by the FS approach. This cell conversion has to be via the status of a liquid-interface cell (blue, solid-line arrows, Fig. 3) to guarantee a closed interface layer. The conversion from liquid-interface cells to gas or liquid cells is governed by the fill level ϕ . If ϕ < 0 − ϵ the former liquid-interface cell converts to a gas cell and if ϕ > 1 + ϵ the liquid-interface cell converts to a liquid cell. Conversion from gas/liquid cells into liquid-interface cells depend on the neighborhood: if one liquid-interface neighbor converts during the mass exchange into a liquid or gas cell, the current liquid/gas cell has to convert into a liquid-interface cell. A careful consideration of the solid–liquid phase transition [32] and the related cell conversions (Fig. 3, red-dashed arrows) is necessary for the EBM process. This phase transition is controlled by the energy density E in Eq. (20). If E exceeds a material specific value the former solid cell converts to a liquid cell. The conversion goes the other way round if E falls below the other specific threshold. The change in E is governed by the computation of the thermal LBGK of Eqs. (23) and (24) in the solid cells where the fluid velocity u is set to zero. In order to guarantee a stable behavior of thermal cell conversions the boundary between two cell states is not a fixed value. It is implemented as a continuous hysteresis region, i.e., the phase transition between solid and liquid is done between a melting range between the solidus and liquidus temperature depending on the material of the simulated powder. The conversion from liquid or solid cells to liquid- or solid-interface cells does not require any additional initialization steps. For conversions from gas into liquid/solid-interface cells, however, a reinitialization of the pdfs fi and hi is necessary, as this data is undefined for the gas phase in the free surface approach. Details can be found in [33,1]. Before we start simulating the EBM process with this thermal FS LBGK method we discuss the realization of some further EBM specific properties in the next section. 3. Modeling aspects of EBM In this part modeling aspects regarding the interaction of the electron beam and the metal powder particles are described. We characterize the EB gun in our model, then two absorption types, explain the numerical generation of metal powder particles and show that numerical and real powder distribution results are highly concordant. 3.1. Modeling of the electron beam gun The electron beam gun is defined by its power P = U · I where U is the acceleration voltage and I the current and its shape. The shape of the EB gun determines how the energy is distributed on the surface of the powder particles on top of the powder bed. A typical beam has Gaussian shape and can be modeled by a general multi-dimensional normal distribution with standard deviation σ . Thus, we use a spatially discrete function of the beam energy Eb [33] with parameter P and σ

R. Ammer et al. / Computers and Mathematics with Applications (

)

–

7

dissipated energy coefficient [1/5e-6m]

Fig. 3. Possible cell conversions.

penetration depth [5e-6m] Fig. 4. Absorption types: relation between dissipated energy and penetration depth.

using a cell of unit length as defined in

 Eb (xi , yj ) = 2πσ 2 exp − 2σ 2 P

1

   2   xi + 0.5 xb     yj + 0.5 − yb  ,

(29)

where (xi , yi )T defines the current position of the beam. Eb does not depend on the z-direction since the EB gun is fixed at the top of the EBM machine and can only move in the (xy)-plane. The deviation of (0.5, 0.5)T in x- and y-direction is related to the cell-centered lattice. The energy absorption begins on the contact area on top of the powder layer where the induction of the beam energy starts. We assume a quite narrow EB and therefore, energy absorption proceeds only in z-direction. There are two different absorption types described in Section 3.2 both depending also on the fill level ϕ of the powder particles. Another essential property of the EB gun is its mobility in x- and y-direction. We define beam movements by a starting point s = (xb , yb )T , a velocity vb and the duration ∆t this movement lasts. This triple also enables us to simulate that the beam is switched off by setting the starting point outside the simulation domain and velocity to zero. Stationary cases of the beam are modeled by a zero velocity vector and line melting by a constant velocity. 3.2. Absorption types We have already discussed how energy only applied on the contact area depends on the beam shape. For the further energy transport in z-direction through the powder particles two types of absorption are available regarding different types of EBM machines. Exponential absorption is modeled for an EB gun of 60 kV and constant absorption for an EB gun of 120 kV. Fig. 4 shows the dissipated energy related to the penetration depth for both absorption types for experimental data from [34]. The amount of dissipated energy from real measured data should match with the approximated model. Hence, we approximate the exponential dissipation energy coefficient by a monotonically decreasing exponential function. In the corresponding absorption type most of the dissipated energy is absorbed by cells which are closer to the EB gun. AssumdE dE ing that the EB gun is quite narrow and all energy is absorbed in z-direction, i.e., dxb = dyb = 0, the exponential energy absorption following the Lambert–Beer law [35] can be written as dEb dz

= −λabs Eb

(30)

where λabs is the absorption coefficient. The energy deposited in each cell under beam radiation can be calculated by

∆Eb = −ϕλabs Eb , ∆x where ϕ denotes the fill level. Using this exponential law the energy source Φi from Eq. (17) can be computed.

(31)

The absorption type for the 120 kV is approximated by a constant function. Consequently, each cell in z-direction gets the same amount of energy until energy is dissipated completely. Thus, the constant absorption type has a higher penetration depth than the exponential one (cf. Fig. 4). For both types, absorption rate depends on the fill level at the surface of the

8

R. Ammer et al. / Computers and Mathematics with Applications (

)

–

0.02 0.018 0.016 0.014 0.012 0.01 0.008 0.006 0.004 0.002 0 0

50

100

150

200

250

Fig. 5. Real and approximated powder distributions.

powder particles. Particles with a smaller fill level absorb less energy at the surface than particles filled completely. The algorithm using the constant absorption rate will not be explained in detail in this work. But both absorption algorithms are parallelized and depend on the fill level. 3.3. Generation of powder Another modeling issue is the generation of randomly distributed powder particles to simulate the interaction between the electron beam and the powder. Starting from the data of three different real powder charges (Fig. 5) an appropriate density function for modeling the radii of particles is needed. We use the following procedure to get a realistic distribution of the powder particles in our simulation. Due to the asymmetric, right-skewed structure (cf. Fig. 5, charge 1–charge 3) of the data points an inverse Gaussian density function is a suitable approximation for the random generation of the powder particle radii

 f (x) =

  λ λ(x − µ)2 exp − . 2π x3 2µ2 x

(32)

The governing parameters of an inverse Gaussian density function are mean value µ, the shape factor λ and skewness factor ν . The relations between the shape parameter λ and mean value µ to the variance and the skewness factor for the inverse Gaussian density are given by

σ2 =

µ3 , λ

and

ν=3



µ . λ

(33)

The skewness parameter ν measures the asymmetry of the density function. Thus, ν should be constant for all relating density functions of bulk powder. The mean value of the density function can be calculated by skewness and the sizes of minimum and maximum radii of the given powder charges. This inverse Gaussian density function is shown in Fig. 5 by the solid (blue) curve. With the inverse Gaussian density function f (x) the cumulative density function F (x) can be computed. After the generation of a random number r, 0 ≤ r ≤ 1 the radius dr of one powder particle can be computed by the quantile of F (x) dr = F −1 (r ).

(34)

This procedure enables us to change the powder distribution by varying the mean size, the size distribution and the layer thickness to test different bulk powder and examine their influence on the melt pool. With these modeling aspects regarding the electron beam and the powder bed generation we can start with validation experiments and numerical simulations of small spots in the next section. 4. Numerical validation and simulation results In this section the implementation in the waLBerla2 framework is explained briefly and our thermal LBM is validated against the solution of the heat equation. Furthermore, the liquid–solid phase transition is validated on the basis of the 2 waLBerla: http://www10.informatik.uni-erlangen.de/en/Research/Projects/walberla/description.shtml.

R. Ammer et al. / Computers and Mathematics with Applications (

)

–

9

Stefan problem [36]. Finally, a physical meaningful simulation of melting a spot in a powder particle bed is shown and the behavior of the melt pool examined. 4.1. waLBerla & pe The implementation of the EBM process is integrated in the waLBerla framework which is a LBM based fluid flow solver [37,38]. We add thermodynamic properties regarding the EBM process to it, e.g., melting and solidification due to the absorption of energy. This is possible because of its modular architecture. For the generation of powder particles with the radius distribution explained in Section 3.3 we use the pe 3 framework which is described in [39]. The pe represents the particles as rigid bodies which are treated as boundaries to the LBGK-scheme of waLBerla . waLBerla executes the LBM steps of Eqs. (23) and (24) regarding the flow dynamics of the melt pool and can handle the free surface approach described in Section 2.3. If more powder layers are simulated the process starts again: the surface of the last layer (consisting of unmelted and melted solid particles) is approximated and then the next powder layer is generated by the pe . We extend the already parallel implemented LBM kernels provided by waLBerla for the new thermal LBGK to ensure sophisticated communication concepts using MPI. 4.2. Validation by heat equation We demonstrate the correct numerical transfer and implementation of the heat equation in our system. Therefore a small simulation example (sketched in Fig. 6) is examined. A box is filled with fluid and the initial temperature of the fluid is 2000 K. The wall in the (yz )-plane is heated up to a temperature of 3000 K. We use the boundary conditions which are explained in Sections 2.1 and 2.2. The open thermal boundary sets the temperatures at the two walls and for the other walls we use the closed thermal boundary. The no-slip boundary condition is set for the hydrodynamic flow field. The corresponding analytical solution of the heat equation is given by: g (xj , ∆T , L, k, tj ) =

∞ 2∆T  (−1)n

π

n =1

n



sin nπ

xi  L

  ∆T 2 2 tj exp −kn π 2 + xj + T0 L

L

(35)

2

with: ∆T = 1000 K, L = 625 · 10−6 m and thermal diffusivity k = 9.93 · 10−6 ms . Fig. 7 shows analytical and numerical solutions of Eq. (35) for three different times. Table 1 gives the L1-error norms of analytical and numerical solutions. They are in the range of 10−9 and indicate that modeling and implementation of heat equation is correct and the resolution is appropriate. This numerical example has shown that the model and implementation regarding the solution of the heat equation works well. In the next section we validate the correctness of the liquid–solid phase transition by the Stefan problem. 4.3. Validation of solid–liquid phase transition by the Stefan problem The original Stefan problem (SP) [36] treats the formation of ice in the Arctic Ocean and can be seen as a special solid–liquid boundary problem where the interface of these two phases can move with time. Hence, the SP is a benchmark problem for melting as well as for solidification processes. The conservation of energy governs the corresponding phase changing process. The following parameters are involved in the description of the SP, phase change material – in our application is this the metal powder – with constant density ρ , latent heat (L), liquidus temperature Tl , specific heat capacity cp and thermal conductivity k. The dimensionless Stefan number characterizes the specific phase transition problem by the ratio of St =

cp · (T − TL )

∆H

.

(36)

The validation setup is a cubic box filled with melted metal and is cooled at one wall boundary where a linear temperature profile is assumed. This cooling wall has a fixed temperature (using the open thermal boundary condition) which is below the solidification temperature of the metal. The solidification of the liquid metals begins at the interface between cooled wall and the liquid metal. The resulting analytical depth of the solidified metal is the function h(t ) depending on time. Stefan shows in [36] that h(t ) can be given by

 h=

2λat . ∆H ρ

(37)

Fig. 8 shows the analytical solution of Eq. (37) and the numerical solution of our simulation with the values given in Table 2. The minimum relative error of both solutions for h(t ) is 0.0062. This validation experiment of the Stefan problem

3

pe : http://www10.informatik.uni-erlangen.de/en/Research/Projects/pe/.

10

R. Ammer et al. / Computers and Mathematics with Applications ( Table 1 L1-error of numerical and analytical solution. Time (ms)

∥err∥1

0.88 2.63 8.77

3.1487 · 10−9 3.0684 · 10−9 3.0791 · 10−9

Fig. 6. Sketch of validation scenario.

Fig. 7. Comparison of numerical and analytical solutions.

Fig. 8. Comparison of numerical and analytical solutions.

)

–

R. Ammer et al. / Computers and Mathematics with Applications (

)

–

11

Table 2 Data for the Stefan problem. Parameter

Value

Unit

Latent enthalpy

∆H

2.9 · 10

T − Ts Heat conductivity

a

λ

65 3.0

Specific heat capacity

cp

7.0 · 10

Density

ρ

4000

Stefan number

St

0.15689

5

m2 s2

K kg·m s3 K 2

m2 s2 K kg m3

–

Table 3 Physical and simulation data to spot melting. Physical data Domain tbeam tsolidify rbeam

Simulation data

(0.96×0.96×0.48) mm

3

1.80 ms 0.90 ms 0.405 mm

192 × 192 × 96 (lattice cells) 15000 (lattice time steps) 7500 (lattice time steps) 81 (lattice cells)

indicates that our phase-transition treatment is physically correct. After validation examples we proceed to numerical simulations regarding the EBM process in the following section. 4.4. Numerical results We show a typical EBM simulation, the melting of a spot in one powder layer. This process can be repeated layer by layer and a link is produced by additive manufacturing. This simulation has been done at the LiMa4 compute cluster, RRZE (Germany). In the scenario a square powder bed (Table 3) is melted by a (60 kV, 5 mA)-electron beam. The electron beam is switched on for 7500 lattice time steps and the system has the same time to solidify. Here the pe -supported powder particle as well as the corresponding interaction of the electron beam gun can be observed. In Fig. 9(a)–(f) the solid powder particles are colored gray and the temperature is visualized by a blue–red scaling, i.e., dark blue corresponds to a temperature of 1500 K (melting temperature of the powder particles) and dark red a temperature of 4500 K. In order to gain a better view of the melting a cross-section of the powder bed and melt pool is visualized. Based on the alignment of the scenario a melt pool with a hemispherical form is expected. Fig. 9(a) shows the zeroth time step: the electron beam is off, all powder particles are solid and have the same initial temperature. Fig. 9(b) shows the melt pool after 2000 time steps; parts of the first powder particle layer are melted. 2000 time steps later (Fig. 9(c)) the melt pool is enhanced, deeper and the temperature is higher. As expected, the temperature is distributed radially and the melt pool has the shape of a hemisphere. Fig. 9(d) visualizes the last time step with working beam; the melt pool is increased and the second powder layer is melted. 3000 time steps later (Fig. 9(e)) it can be observed that the melt pool is still increased due to the stored energy although the electron beam is already switched off. Eventually, Fig. 9(f) shows the still liquid melt pool and it can be seen that a material specific wetting angle occurs at the boundary of the melt pool. This validation and the numerical simulation in this section show the physical correctness of modeling and implementation. 5. Conclusion Summarizing, this paper describes a parallel, thermal free surface LBGK method which is appropriate for simulating the EBM process. The implementation is included in the waLBerla and pe framework. The full algorithm was implemented to support parallel runs in order to make 3D real world scenarios possible on large scale machines. Our model characterizes the electron beam gun by acceleration voltage, current and shape and enables various movements of the beam gun to create complex structures layer by layer. We model two absorption types – exponential and constant – regarding different EBM machines which vary in penetration depth and dissipated energy. We validate our thermal LBGK model by the comparison of analytical and numerical solution of the heat equation and show that the corresponding L1-error is sufficiently small. In order to show the physical correctness of the solid-phase transition in our model we consider the Stefan problem to validate the solidification process of the metal. Finally, the numerical melting of a spot in a simulated powder bed of particles demonstrate the power of our model and implementation: we are able to simulate EBM melting scenarios of meaningful size to gain a better knowledge of the process which is our target objective. Research topics in the future will be validation and parameter studies regarding the size of the melt pool, its maximum temperature and the porosity of the final products. Another objective will be the implementation of a temperature

4 LiMa: http://www.rrze.uni-erlangen.de/dienste/arbeiten-rechnen/hpc/systeme/lima-cluster.shtml.

12

R. Ammer et al. / Computers and Mathematics with Applications (

(a) t = 0.

(b) t = 2000.

(c) t = 4000.

(d) t = 7500.

(e) t = 10000.

(f) t = 15000.

)

–

Fig. 9. Melt pool with a (uB = 60 kV, I = 5 mA)-beam. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

dependent surface tension. In this work the surface tension is constant and independent of the temperature but especially for higher temperatures the influence of temperature on the surface tension is significant. Acknowledgments Our work is supported by the European Union Seventh Framework Program—Research for SME’s with full title ‘‘High Productivity Electron Beam Melting Additive Manufacturing Development for the Part Production Systems Market’’ and grant agreement number 286695. The authors of this paper would like to thank all developers of waLBerla and pe . References [1] C. Körner, E. Attar, P. Heinl, Mesoscopic simulation of selective beam melting processes, Journal of Materials Processing Technology 211 (6) (2011) 978–987. [2] G. McNamara, C. Zanetti, Use of the Boltzmann equation to simulate lattice-gas automata, Physical Review Letters 61 (1988) 2332–2335. [3] X. He, L.-S. Luo, Theory of the lattice Boltzmann method: from the Boltzmann equation to the lattice Boltzmann equation, Physical Review E 56 (1997) 6811–6817. [4] S. Chen, G.D. Doolen, Lattice Boltzmann method for fluid flows, Annual Review of Fluid Mechanics 30 (1998) 329–364. [5] X. He, L.-S. Luo, A priori derivation of the lattice Boltzmann equation, Physical Review E 55 (1997) R6333–R6336. [6] P. Bhatnagar, E. Gross, M. Krook, A model for collision process in gases 1: small amplitude processes in charged and neutral one-component systems, Physical Review 50 (1954) 511–525.

R. Ammer et al. / Computers and Mathematics with Applications (

)

–

13

[7] L.-S. Luo, Theory of the lattice Boltzmann method: lattice Boltzmann models for non-ideal gases, Physical Review E 62 (2000) 4982–4996. [8] H. Chen, S. Chen, W. Matthaeus, Recovery of the Navier–Stokes equations using a lattice-gas Boltzmann method, Physical Review A 45 (1992) R5339–R5342. [9] Y. Qian, D. d’Humières, P. Lallemand, Lattice BGK models for Navier–Stokes equation, Europhysics Letters 17 (1992) 479–484. [10] Z. Guo, C. Zheng, B. Shi, T.S. Zhao, Thermal lattice Boltzmann equation for low Mach number flows: decoupling model, Physical Review E 75 (2007) 036704–036719. [11] X. Shan, Simulation of Rayleigh–Bénard convection using a lattice Boltzmann method, Physical Review E 55 (1997) 2780–2788. [12] B. Palmer, D. Rector, Lattice Boltzmann algorithm for simulating thermal flow in compressible fluids, Journal of Computational Physics 161 (2000) 1–20. [13] Y. Zu, Y.Y. Yan, W. Shi, L. Ren, Numerical method of lattice Boltzmann simulation for flow past a rotating circular cylinder with heat transfer, International Journal of Numerical Methods for Heat & Fluid Flow 18 (2008) 766–782. [14] F.J. Alexander, Lattice Boltzmann thermohydrodynamics, Physical Review E 47 (1993) R2249–R2252. [15] C. Teixera, H. Chen, D.M. Freed, Multi-speed thermal lattice Boltzmann method stabilization via equilibrium under-relaxation, Computer Physics Communications 129 (2000) 207–226. [16] M. Watari, M. Tsutahara, Possibility of constructing a multispeed Bhatnagar–Gross–Krook thermal model of the lattice Boltzmann method, Physical Review E 70 (2004) 016703–016712. [17] G. McNamara, B. Alder, Stabilization of thermal lattice Boltzmann methods, Journal of Statistical Physics 81 (1995) 395–408. [18] F. Massaioli, R. Benzi, S. Succi, Exponential tails in two-dimensional Rayleigh–Bénard convection, Europhysics Letters 21 (1993) 305. [19] X. He, S. Chen, G. Doolen, A novel thermal model for the lattice Boltzmann method in incompressible limit, Journal of Computational Physics 146 (1998) 282–300. [20] V. Alexiades, A.D. Solomon, Mathematical Modeling of Melting and Freezing Processes, Hemisphere Publications, 1993. [21] S. Chen, G.D. Doolen, Lattice Boltzmann method for fluid flows, Annual Review of Fluid Mechanics 30 (1) (1998) 329–364. [22] A.K. Gunstensen, D.H. Rothman, S. Zaleski, G. Zanetti, Lattice Boltzmann model of immiscible fluids, Physical Review A 43 (8) (1991) 4320. [23] W. Boettinger, J. Warren, C. Beckermann, A. Karma, Phase-field simulation of solidification 1, Annual Review of Materials Research 32 (1) (2002) 163–194. [24] Y. Chunyu, M.W.G. Teng, X.D.Z.J.Y. Rui, W. Yunzhi, 3D Phase-field simulation of effect of interfacial energy anisotropy on sideplate growth in Ti–6Al–4V, Acta Metallurgica Sinica 48 (2012) 148–158. [25] W. Mulder, S. Osher, J.A. Sethian, Computing interface motion in compressible gas dynamics, Journal of Computational Physics 100 (2) (1992) 209–228. [26] S. Osher, A level-set formulation for the solution of the Dirichlet problem for Hamilton–Jacobi equations, SIAM Journal on Mathematical Analysis 24 (5) (1993) 1145–1152. [27] D. Peng, B. Merriman, S. Osher, H. Zhao, M. Kang, A PDE-based fast local level-set method, Journal of Computational Physics 155 (2) (1999) 410–438. [28] M. Latva-Kokko, D.H. Rothman, Diffusion properties of gradient-based lattice Boltzmann models of immiscible fluids, Physical Review E 71 (5) (2005) 056702–056709. [29] B. Nestler, A. Choudhury, Phase-field modeling of multi-component systems, Current Opinion in Solid State and Materials Science 15 (3) (2011) 93–105. [30] U. Rüde, N. Thürey, Free surface lattice-Boltzmann fluid simulations with and without level sets, VMV 2004 (2004) 199–208. [31] C. Körner, M. Thies, T. Hofmann, N. Thürey, U. Rüde, Lattice Boltzmann model for free surface flow for modeling foaming, Journal of Statistical Physics 121 (2005) 179–196. [32] D. Chatterjee, S. Chakraborty, A hybrid lattice Boltzmann model for solid–liquid phase transition in presence of fluid flow, Physical Letters A 351 (2006) 359–367. [33] E. Attar, Simulation of selective electron beam melting processes, Ph.D. Thesis, 2011. [34] K. Kanaya, S. Okayama, Penetration and energy-loss theory of electrons in solid targets, Journal of Physics D: Applied Physics 5 (1972) 43–58. [35] K. Fuwa, B. Valle, The physical basis of analytical atomic absorption spectrometry. The pertinence of the Beer–Lambert law, Analytical Chemistry 35 (8) (1963) 942–946. [36] J. Stefan, Über die Theorie der Eisbildung, insbesondere über die Eisbildung im Polarmeere, Annalen der Physik 278 (2) (1891) 269–286. [37] C. Feichtinger, S. Donath, H. Köstler, J. Götz, U. Rüde, Walberla: HPC software design for computational engineering simulations, Journal of Computational Science 2 (2) (2011) 105–112. [38] H. Köstler, U. Rüde, The CSE software challenge—covering the complete stack, IT-Information Technology 55 (3) (2013) 91–96. [39] K. Iglberger, U. Rüde, Large-scale rigid body simulations, Multibody System Dynamics 25 (1) (2011) 81–95.