A coupled Cellular Automaton–Lattice Boltzmann model for grain structure simulation during additive manufacturing

Computational Materials Science 124 (2016) 37–48

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

A coupled Cellular Automaton–Lattice Boltzmann model for grain structure simulation during additive manufacturing Abha Rai ⇑, Matthias Markl, Carolin Körner Chair of Materials Science and Engineering for Metals, University of Erlangen-Nuernberg, Martensstr. 5, D-91058 Erlangen, Germany

a r t i c l e

i n f o

Article history: Received 25 April 2016 Received in revised form 1 July 2016 Accepted 5 July 2016

Keywords: Lattice Boltzmann method Cellular Automata Additive manufacturing Powder bed Grain structure

a b s t r a c t A 2D coupled Cellular Automaton (CA)–Lattice Boltzmann (LB) model has been developed to simulate grain structure evolution during powder–bed–based, layer by layer, additive manufacturing (AM). The presented model includes algorithms for random powder layer generation, electron beam energy absorption, evaporation, capillarity and wetting, meltpool dynamics, its temperature evolution and face centered cubic grain solidification. The model is first validated against the experimental findings of single track electron beam melting and solidification of a baseplate and then applied to simulate the grain structure produced during AM. Influence of the hatching strategy on grain structure is presented as well as stray grain formation resulting from partially molten powder particles. Ó 2016 Elsevier B.V. All rights reserved.

1. Introduction Powder-bed-based additive manufacturing (AM) is a state-ofthe-art technique of layer by layer building of components within a powder layer. Powder bed is applied on baseplate and a heat source, e.g., a laser or an electron beam, is used to selectively melt the powder. Individual meltpools develop and solidify to form first layer of the final part. After complete solidification of the first layer the baseplate is lowered vertically and the next powder layer is spread over the previous layer. This process is repeated until the final part is built. Depending on the heat source this process is referred to as selective laser melting (SLM) [1,2] or selective electron beam melting (SEBM) [3]. The final grain structure depends on the local thermal history and varies within the component [4,5]. Defects like residual porosity, partially molten powder particles and not well connected layers have been observed, which affect the grain structure evolution during solidification [6–8]. Partially molten powder particles act as nuclei and lead to stray grain formation in the bulk of the sample [4]. A surface roughness much larger than expected from the mean powder particle size is also observed [9] which originates from meltpool dynamics. The surface roughness, texture and mean grain width of the final component are complex functions of the process and material parameters [10]. The most influential process parameters are the beam power, speed, line offset (spacing between two ⇑ Corresponding author. E-mail address: [email protected] (A. Rai). http://dx.doi.org/10.1016/j.commatsci.2016.07.005 0927-0256/Ó 2016 Elsevier B.V. All rights reserved.

consecutive beam tracks), beam movement contours and the return time of the beam with respect to a fixed point [11]. The powder consolidation mechanism depends on the choice of material [12]. In case of metals it is mainly determined from the sequence of melting–solidification–remelting. For a given material and geometry of the component, choice of the optimum process parameters is still based on trial and error method. Therefore, a numerical model to simulate grain structure and texture during selective beam melting is needed. The numerical model should include physical phenomena like, beam energy absorption and losses in the powder layer, melting and solidification, meltpool dynamics, wetting, capillary effects, gravity, heat transfer (within the powder layer and meltpool) and competitive grain growth during solidification. Additionally, the stochastic effects arising due to the powder particle size distribution and variations of the local powder particle density should be taken into account. It must be emphasized here that the inclusion of the individual powder particle related effects play a very significant role in the prediction capabilities of the SEBM simulations. The net energy transfer to the powder bed is determined by the actual thermal conduction between individual powder particles which is proportional to their contact area. The wetting process of individual powder particles depends on their geometry and plays a very significant role during the consolidation process of the powder bed. Moreover, influence of the gas porosity from the powder bed and grain nucleation due to partially molten powder particles on the final built part can only be incorporated into the model when individual powder particles are considered. All these

38

A. Rai et al. / Computational Materials Science 124 (2016) 37–48

effects cannot be simulated using a homogenized approach of the powder bed. In the context of SLM numerous 3D models based on the finite element (FE) or finite difference (FD) methods can be found in literature to solve the global temperature fields and meltpool dynamics [13–17]. LB models have been presented by Körner et al. [12,18,19] and Ammer et al. [20] to simulate SEBM process in 2D and 3D, respectively. 3D CAFE models have been applied to study grain structure evolution during directional solidification [21,22], welding [23] and laser metal deposition [24,25]. The existing models either take into account the powder related stochastic aspects of the process (without the inclusion grain structure evolution) [12,18–20], or they predict the grain structure evolution based on homogenized approximation of the powder bed [21,22,24,25] but none of them include both, i.e., powder related effects and grain structure evolution into a single model. Aim of this paper is to simulate grain structure evolution during SEBM considering the individual powder related effects. For this purpose a CA based grain growth model will be coupled to the existing LB based meltpool dynamics model. Due to the computational times needed to simulate the meltpool dynamics with individual powder particle effects based model the multi-layer hatching simulations are limited to a few layers only [20]. Since the goal of the present work is to predict grain structure evolution for a multi-layer hatching with typical experimental build processes of at least 40–50 layers, the 2D LB model of Körner et al. [18] model is used in the present work to simulate meltpool dynamics. The resulting temperature field is transferred to a CA based grain structure evolution model [26]. The main advantage of the present model is its capability to simulate a multi-layer hatching process (
5 mm) with inclusion of individual powder related effects on the final grain structure. It will be shown in the results sections that 2D grain orientation histograms obtained from the presented model give a reasonable qualitative description of the texture evolution. The resulting tool is able to predict the asbuilt grain structure for a given material and beam parameters. For testing the model, powder particles and baseplate made from nickel based superalloy IN718 has been chosen which is well suited for applications under extreme temperature condition, e.g., in aerospace industry [27,28]. 2. Numerical model Fig. 1 shows the main physical processes present during SEBM of a powder layer. At the beginning, temperature of the material is lower than solidus temperature T solidus and the whole material is in solid phase. As the beam approaches the absorbed energy

from the beam causes a rise in temperature of the material, when it exceeds T solidus , the solid–liquid phase transformation sets in. For temperatures higher than the liquidus temperature T liquidus only liquid phase exists. The simulation domain consists of square cells and each cell can be categorized as, gas, liquid, solid, solid–liquid interface or gas interface cell. In the following part of this section main methods used to calculate the heat transfer, meltpool dynamics and grain growth are presented. Since the first two models have already been published a very short overview is provided. The emphasis is given to explain the solidification microstructure evolution model and its coupling to the meltpool dynamics model. 2.1. Meltpool dynamics The thermodynamic and hydrodynamic evolution of meltpool is numerically calculated using free surface, multi-distribution LB method. The main effort is focused to compute the evolution of fluid flow velocity v and thermal energy density E,

Z

E¼

T

q cp ðTÞ dT þ q DH;

ð1Þ

0

where q; T and cp denote mass density, temperature and specific heat at constant pressure, respectively. The first term represents the change in energy density due to the temperature change and DH is the latent enthalpy of a cell undergoing a phase change. The conservation equations are

rv ¼0

ð2Þ

@v 1 þ ðv  rÞv ¼  rp þ mr2 v þ F @t q

ð3Þ

@E þ r  ðv EÞ ¼ r  ðkðTÞrEÞ þ U: @t

ð4Þ

Here t; p; m; k and F denote time, pressure, kinematic viscosity, thermal diffusivity of the material and external forces like gravity, respectively. U is the source term for energy density conservation equation (Eq. (4)), e.g., the energy deposited by the electron beam. Since no grain movement is considered in the presented model, the solid is assumed to have a zero velocity. Temperature of a cell is derived from the calculated energy density E using the following relations

8 E=qcps ; E 6 qcps T s > > <   Eqcps T S T ¼ T S þ ðT L  T S Þ  qcps ðT L T S ÞþL ; qcps T S < E < qcps T L þ L > > : T L þ ðE  L  qcps T L Þ=qcpl ; E P qcps T L þ L:

Fig. 1. Schematic of the physical phenomena taking place during SEBM of a powder layer.

ð5Þ

39

A. Rai et al. / Computational Materials Science 124 (2016) 37–48

here L is the latent heat of melting, cps and cpl are the specific heat of the solid and liquid phases, respectively. When a cell solidifies the released latent heat modifies the energy density of the cell and the temperature field is updates as given above. For further details about the used LB model the reader is referred to [12,18]. When the beam melts the powder particles, wetting describes how the resulting liquid spreads out on the underlying solid surface. Capillary forces develop if the surface curvature of a liquid– gas interface does not vanish. The phenomena of capillarity and wetting are governed by surface and interface energies. These are strongly correlated with each other and play a crucial role during SEBM. The model uses the approach followed by Attar et al. to incorporates surface tension related effects of wetting and capillarity [29]. Wetting is treated as a kind of additional capillary force which vanishes when the dynamics wetting angle equals the equilibrium wetting angle. The liquid–gas interface is modeled using a volume of fluids method. Due to fluid motion, the fluid fraction of an interface cell changes. When an interface cell is fully filled or empty it is converted to a fluid or gas cell, respectively and the fluid–gas interface is also shifted accordingly. The effect of surface tension related forces is incorporated in Eq. (3) while applying the boundary conditions at the interface cells as a local modification of the gas pressure. Depending on the local curvature j at the interface between gas and liquid, the actual gas pressure p0G acting on the gas-interface differs from the ambient pressure pG and is given by

p0G ¼ pG  j  rLG

ð6Þ

here rLG represents the surface energy at the interface between liquid and gas. The curvature of the liquid–gas interface cells is determined using a discrete template-sphere method which uses an extended neighborhood of 25 cells. A special treatment is needed at triple points where fluid, gas and solid intersect. For further details about the implementation of the dynamic wetting and curvature calculations the reader is referred to [29]. In the present work gas porosity is not incorporated into the model and the equilibrium wetting angle between the fluid and solid is chosen to be 0, assuming perfect homologous wetting, therefore melting of the power bed results into a continuous meltpool. The value of surface tension and temperature dependent surface tension coefficient for IN718 can be found in [30,31]. An energy absorption model based on semi-empirical equations and theoretical considerations, is reported by Klassen et al. [19]. The beam is assumed to have a 2D Gaussian profile with standard deviation of r (beam width w ¼ 4 r). For an electron beam with beam current I and electron acceleration voltage U, the total energy deposited per unit time step Dt is given by, Etotal ¼ U  I  Dt. Fig. 2 shows the beam movement, it moves in xz plane, either along or across the grain structure simulation plane (xy plane). The amount of energy deposited in each cell as a function of time is given by:

Ecell ðx; y; z; tÞ ¼ Etotal  gðx; z; tÞ  hðyÞ

ð7Þ

hðyÞ and gðx; z; tÞ are the energy distribution functions in the vertical and horizontal directions, respectively. When the beam moves outof-plane (in z direction) with a deflection speed v z , the beam intensity profile in xz plane is given by:

# ðx  x0 Þ2 þ ðzðtÞ  z0 Þ2 ; gðx; z; tÞ ¼  exp  2pr2 2r 2 1

"

ð8Þ

here x0 and z0 are the co-ordinates of the beam’s focal spot and zðtÞ ¼ v z t. The value of hðyÞ depends on beam penetration which is determined from the combination of beam, material and powder properties. When electrons with incident energy E0 impinge on a target at an angle / with respect to the surface normal, the fraction of energy absorbed by the nth layer in y-direction is given by



   1 DE A 1 EA ynþ1 ; Z; E0 ; / þ EA yn ; Z; E0 ; / hðyn Þ ¼  ¼  ; f v  Dy E0 f v  Dy E0

ð9Þ

here Dy is the spatial resolution in y-direction, f v is the volume fraction of cell occupied by the material with atomic number Z; EA =E0 is the energy fraction absorbed within the target yn ¼ yn cos / and ynþ1 ¼ ðyn þ f v  DyÞ cos /. The presented energy absorption model has been validated against the experiments for SEBM and the details can be found in [19]. The calculation of evaporation effects and recoil pressure are performed based on an earlier work from our group [32]. 2.2. Solidification grain growth The underlying model for grain growth is based on the work of Gandin et al. [26]. Ni-based superalloys solidify with a face-centered cubic (fcc) lattice structure. Growth of a single fcc dendrite in an undercooled melt exhibits an octahedral grain morphology as observed in experiments [33] and simulations [34–36]. The main half diagonals of the octahedron are defined by <100> directions. The trunks and arms of the dendrite are also aligned along these directions. Under uniform thermal conditions all the diagonals of the octahedron have the same length (Fig. 3a). For mapping this 3D grain morphology to 2D it is assumed that one of the <100> direction of the grain is perpendicular to the 2D simulation plane (e. g. ½0 0 1 direction). Under uniform thermal conditions the grain envelop resembles a square in 2D with trunks and arms of the dendrites aligned along the <10> directions (Fig. 3b). For non-uniform thermal field e.g., a thermal gradient directed along the y-axis, the grain envelop looks as shown in Fig. 3c. For describing the basic algorithm, the growth of a single undercooled nucleus in a meltpool with uniform temperature is presented. The same physical and thermal parameters are used as in [26]. The simulation domain of size Lx  Ly is divided into nx and ny squares cells of size Dx in x and y directions, respectively. The temperature of a cell is denoted as T i;j , where i 2 f0; . . . ; nx  1g and j 2 f0; . . . ; ny  1g. The undercooling is defined with respect to the liquidus temperature, DT i;j ¼ T liquidus  T i;j . All cells are initialized with zero undercooling. The incubation time is ignored and the initial shape of the growing grain at the nucleation cell, denoted as g in Fig. 4a, is assumed to be a square with its ½1 0 direction aligned at an angle h ¼ 30 with respect to the x-axis. The square encloses the trunks and arms of the dendrite aligned along the crystallographic directions <10>. Depending on the actual local temperature, the dendritic trunks and arms grow. The overall grain envelop also evolves accordingly. Based on the model proposed by Kurz et al. [37,38] the relationship between the dendritic tip velocity (v tip ) and the local undercooling (DT) takes the form:

v tip ¼ AðDTÞ2 , where A is a material dependent constant calculated

either by analytical or more accurate phase field models [39]. An analytical model to determine the relationship between v tip and DT for directional solidification of IN718 is presented in [40,41]. The nucleation cell g has four first and second nearest neighbors denoted as l1 ; . . . ; l4 and l5 ; . . . ; l8 , respectively. Each CA cell is characterized by three parameters, orientation, grain identification number and length of the square. The first two parameters are defined when the cell begins to solidify and the third parameter increases with time depending on the local undercooling. Half size of a nuclei placed on cell g at the time of nucleation (t ¼ tN ¼ 1) is denoted as LtgN and the corresponding square is denoted as S1g in Fig. 4a. The system is then cooled down with a cooling rate of 0.1 K s1. Depending on local undercooling the square belonging to cell g grows with time and Ltg is updated in accordance with Eq. (10).

40

A. Rai et al. / Computational Materials Science 124 (2016) 37–48

Fig. 2. The Gaussian distributed beam moves in xz plane, either in-plane or out-of-plane with respect to the grain structure simulation domain (meshed 2D plane). Each color represents a specific grain orientation. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 3. The grain envelop under: uniform thermal conditions in 3D (a) and 2D (b), non-uniform thermal conditions in 2D (c).

Fig. 4. Schematic of the 2D decentered square algorithm [26]. (a) Growth of the square centered on nucleation cell g at time step 1, 5, 10, 15, 20, 22, 25 (S1g  S25 g ). At time step 25 22 all the neighbors of cell g are captured, therefore, S22 g and Sg completely overlap each other. (b) Capture of first neighbor cells l1 ; . . . l4 by the nucleation cell g. The distances Lt½face1 and Lt½face2 for cell l1 are also depicted. (c) Squares placed on all the second nearest neighbors of cell g at the time of capture. The cells l5 ; . . . ; l8 , their centers C l5 ; . . . ; C l8 and the overlapping corners K l5 ; . . . ; K l8 are marked.

t 1 X Ltg ¼ pffiffiffi v tip ½DT g ðt0 Þ  Dt0 : 2 t0 ¼tN

ð10Þ

g

The constant A relating the dendritic tip velocity to the under4

Fig. 4b, Lt½face1 and Lt½face2 are the distances of the cell center l1 from  1  faces of the capturing square, S17 , respectively. The ½11 and ½1

2

cooling is chosen as 10 m/(s K ). Fig. 4b depicts the capture of the neighboring cell labeled l1 by the nucleation cell g. As soon t

17

as center of the neighbor cell l1 lies within Sg (Sg in Fig. 4b), cell l1 is considered to be captured by cell g. The grain identification number and the orientation of the capturing cell are assigned to the captured cell. The state of the newly captured cell is changed from solid to liquid and a new square is placed on it. The procedure to place the new square is as follows. At the time of capture Lt½face1 and Lt½face2 are the distances along the dendrite front faces which has captured the cell on both sides of its center. For example, when cell l1 is captured by cell g in

length of the new square placed on the recently captured cell is given by

  pffiffiffi  pffiffiffi  Ltl1 ¼ min Lt½face1 ; 2Dx þ min Lt½face2 ; 2Dx :

l1

ð11Þ

The truncation of the growing square according to Eq. (11) ascertains that the size of squares is limited to the maximum size necessary to capture the eight nearest neighbors. This condition is essential for non uniform thermal conditions where the local temperature has to be taken into account. The next step is to find the corner of the capturing square closest to the captured cell center. Then a square of size Ltl1 and orientation h is placed on cell l1 such that one of its corner coincides with the closest corner of the

A. Rai et al. / Computational Materials Science 124 (2016) 37–48

capturing square. At least two out of four faces of the squares belonging to the capturing and captured cell overlap each other. In the present example, the grain is growing under a uniform thermal field. In this case all the first four neighbor cells l1 ; . . . ; l4 are captured at the same time (Fig. 4b, at time step 17). Both, Lt½face1 and pffiffiffi Lt½face2 are smaller than 2Dx and no truncation of the squares is needed. The squares placed on the newly captured cells

l1 ; . . . ; l4 have the same size as the one belonging to cell g and they fully overlap each other. At time step 22, the next nearest neighbors l5 ; . . . ; l8 are also captured simultaneously. Fig. 4c shows the four truncated squares placed on the next nearest neighbors, their centers C l5 ; . . . ; C l8 and the overlapping corners K l5 ; . . . ; K l8 . It can also be noticed that the centers of the squares (also growth centers) belonging to cells l5 ; . . . ; l8 no longer coincide with that of the capturing cell g giving the name decentered square growth algorithm. From now onwards, all newly captured cells grow according to the growth kinetics determined by their local undercooling. Growth of Stg continues until all eight nearest neighbors of g are captured (S22 g in Fig. 4a). In few cases it might happen that more than one cell try to capture a particular cell at the same time. For such scenarios, one of the following pre-defined conditions can be used to decide which capture is allowed (1) the capturing cell with highest undercooling, (2) the capturing cell belonging to the bigger grain or (3)

41

depending on the iteration order through the capturing cells. Choosing one of the first two criteria is physically more intuitive but implies higher memory usage. Numerical tests show that choosing any of these options doesn’t substantially influence the final outcome of the grain growth. This is due to the grain selection process embedded into the model. All the simulations presented in this paper are performed using the last criteria, i.e., the cell which comes first in the iteration order is allowed to capture. Fig. 5a–d shows how the nucleation cell situated at the center of the simulation domain and oriented at an angle of 30 with respect to the x-axis captures the neighboring cells. Each captured/solidified cell, its center and the corresponding square is depicted with the same color. In some cases, at the time of capture, the squares belonging to the captured and capturing cells completely overlap each other (e.g. in Fig. 5a when the first neighbor cells are captured by the nucleation cell). Since the growth velocity of the squares is determined from the local undercooling of the associated cell the squares grow with different velocity and do not overlap anymore at further time steps. Fig. 5d shows the grain envelop given by the superposition of all the squares. 2.3. Baseplate and powder layer generation The CA grain growth model is used to generate baseplates with equiaxed, columnar or single crystalline grain structures. The grain

Fig. 5. Cells captured by the nucleation cell (at the center of domain) at time-step (a) 17, (b) 22, (c) 30 and (d) 34. Each captured cell, its center and the corresponding square is depicted with the same color. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

42

A. Rai et al. / Computational Materials Science 124 (2016) 37–48

size distribution and texture of the initial baseplate is determined experimentally. This input is used to calculate the number of nuclei to be distributed in the baseplate and the range of their orientations. All the cells belonging to the baseplate are initialized as liquid cells and rest as gas cells. Guided by the experimental analysis of the baseplate, a given number of nuclei are stochastically placed in liquid cells. The minimum distance among two nuclei is about 10 lm. For generating an equiaxed baseplate, nuclei with grain orientations ranging from ½0 ; 90  are distributed randomly within the baseplate domain and allowed to grow under uniform thermal field with a constant external cooling. Fig. 6 shows an equiaxed baseplate with grain size between 10 lm and 80 lm which lies within the same ballpark as measured by experiments (20–60 lm). For generating a columnar baseplate nuclei with grain misorientations ranging from ½10 ; 10  are distributed at the bottom of the simulation domain and a thermal gradient is applied in the y-direction. This narrow range of orientations is chosen to reduce the computational domain due to the fact that during the grain selection process only grains with a misorientation of
±5° survive. For generating a single crystal baseplate all the cells within it are initialized with the same crystallographic orientation. An algorithm based on the rain drop model for random packing proposed by Meankin et al. [42] is used to generate a random powder layer consisting of circular powder particles with a pre-defined particle size distribution. The basic idea is to allow the particles to roll down with the help of gravity until it finds the state of minimum energy. All the powder particles are place randomly at a specific height along the build direction (much higher than the solid/liquid–gas interface cells) and are allowed to fall downwards until they come into contact with a stationary power particle. When no stationary powder particle is found then the falling powder particle is deposited at the basal line. In order to adjust the packing densities close to their experimental counterparts, some of the powder particles are removed from the originally produced dense powder layer. In the last step, the circular powder particles are mapped into the numerical grid. The details of the powder generation algorithm can be found in [18]. In the present work each powder particle is assumed to contain a single grain, i.e., all the cells within a powder particle are assigned a specific random

crystallographic orientation. Fig. 6 shows the cross-section of a generated powder layer applied on a baseplate. Due to the fact that a circular powder particle is mapped on a square grid, at the contact point of two powder particles an overlap of the order of one cell dimension is inevitable. The baseplate consists of equiaxed grains. This is also reflected in the grain orientation distribution plot in Fig. 6. It can be seen that the grain size distribution in the generated baseplate is in the same range as the baseplates used in the experiments. The temperature field obtained from the LB model is used to simulate the solidification grain structure using the CA model presented in Section 2.2, the resulting coupled model (CALB model) is presented in next section. The resulting phase change information is then transferred back to the LB model. The coupling back of the latent heat released during the solidification of a liquid cell to the meltpool dynamics using LB model is performed while calculating the energy density in the next step and consequently the temperature field using Eq. (5). 2.4. CALB model All the processes presented in Section 2 are coupled to form CALB model and the information flow among the modules is presented in Fig. 7. The main steps are described in the sequential order below. (1) Initialization: The domain is initialized into nx  ny cells along x and y directions, respectively. The build direction is the y direction. For generation of a nyBP cells high baseplate, all the cells within nx  nyBP domain are initialized with liquid cells. After baseplate generation, the temperature of all the solid cells within the simulation domain is set to the pre-heating temperature. Depending on the geometry to be manufactured, beam velocity and line off-set, the beam contour is determined for each layer. The time counter is initialized (t ¼ 0) and the time step for the LB calculations dt LB is defined. (2) Powder bed generation: A powder layer is generated upon the solidified material.

Fig. 6. Cross-section of a 50 lm thick powder layer on a baseplate. The color bar maps the crystallographic orientation with respect to the x-axis (top). Grain size distribution (bottom left) and grain orientation distribution (bottom right) in the baseplate at y ¼ 0:2 mm. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

A. Rai et al. / Computational Materials Science 124 (2016) 37–48

43

Fig. 7. Diagram showing the coupling and flow of information between CA and LB models. The dotted block outlines the grain structure solidification module.

(3) Energy absorption and LB calculations: Depending on the beam contour for the given layer the beam passes over the material. The energy absorption model is used to determine the energy deposited in each cell. The LB equations are solved to give the updated temperature and velocity fields (T i;j and v i;j ). The time counter is updated by dt LB (t ¼ t þ dt LB ). (4) Temperature interpolation: During SEBM the meltpool is very dynamic, e.g., the temperature variation within time step dtLB might be very large leading to extremely high solidification velocities. The following criteria is used to determine the time step for the CA dtCA and it is done in an adaptive way:

dt CA ¼ 0:1  Dx= maxfv ðDT i;j Þg:

ð12Þ

Here v ðDT i;j Þ is the growth velocity of the square belonging to cell ði; jÞ. As a rule of thumb, in order to keep the CA model for solidification stable it is ensured that the squares belonging to the cells do not grow more than 1/10th of the cell size. A linear interpolation is performed to calculate the temperature change within a dt CA time step. All solid cells which have temperatures higher than T liquidus are turned to liquid cells. (5) Grain growth: All the solid cells which have at least one of its neighbor cells as liquid are allowed to grow and the respective square sizes are updated according to the growth algorithm presented in Section 2.2. The growing cells eventually capture their undercooled liquid neighbors. Appropriate

squares are placed on the newly captured cells and their state is changed from liquid to solid. The steps 4–5 are repeated until dt LB P dt CA . Afterwards, LB calculations are performed again and steps 3–5 are repeated until all the cells are solidified. (6) New powder layer: As soon as all the liquid cells are solidified a new powder layer is generated and the steps 2–6 are repeated until the desired number of powder layers are simulated. 3. Validation of the grain growth algorithm The CA model for grain growth is validated against the analytical model proposed by Gandin et al. [43]. They developed an analytical model for grain growth which was validated against the experiments by Ovsienko et al. [33]. Fig. 8 shows the growth of a single nucleus placed at the center of a 10 mm  10 mm domain divided into 200 cells in each direction. All of the cells are initialized with zero undercooling except the cell containing the nucleus, which has an initial undercooling of 2 K and orientation 30 with respect to x-axis. The constant A relating the dendritic tip velocity to the undercooling is chosen as 104 ms1 K2 and the time step Dt ¼ 0:01 s. The nucleus grows under different thermal conditions. In Fig. 8a, it grows under externally applied spatially isothermal cooling. Since the underlying temperature field is uniform, dendrites grow with the same velocity in all directions and the grain envelop is a perfect square all the time and aligned at the initial

44

A. Rai et al. / Computational Materials Science 124 (2016) 37–48

Fig. 8. Grain growth predicted by CA model for a single nucleus growing under (a) spatially isothermal cooling (G ¼ 0:0 K m1 and T_ ¼ 0:1 K s1), (b) static thermal field (G ¼ 250 K m1 and T_ ¼ 0:0 K s1), and (c) Bridgman condition with nx ¼ ny ¼ 200 (G ¼ 250 K m1 and T_ ¼ 0:1 K s1). The grain envelop is plotted at an interval of 1 s.

orientation of the nucleus. Fig. 8b shows the grain growth under static thermal fields. In Fig. 8c the nucleus grows under a Bridgman thermal condition with an external thermal gradient G ¼ 250 K m1 in the y direction and cooling rate T_ ¼ 0:1 K s1. The undercooling increases with decreasing y, leading to higher growth velocities and the grain envelop is not a square anymore, but preserves its initial orientation. The grain envelop calculated using the CA model matches very well with the analytical predictions by Gandin et al. [43]. Favorable oriented grains (FOG, grain B in Fig. 9) are the ones where the <10> direction is better aligned to the heat flow/thermal gradient direction compared to the unfavorable oriented grains (UOG, grain A and C in Fig. 9). According to the classical rule of grain selection at the grain boundaries [44,38], in case of converging dendrites, i.e., the grain boundary between grain A and B, the dendritic tips of UOG (grain A) impinge upon the sides of FOG (grain B) dendrites and cannot grow further. Whereas in case of diverging dendrites, i.e., the grain boundary between grain B and C, secondary and tertiary arms develop from FOG. In both, converging and diverging grains, FOG overgrow UOG. The following simulation is performed to validate the embedded selection rules at the grain boundaries. A simulation domain of 10 mm  50 mm with cell size 50 lm is initialized with liquid cells and then 100 nuclei with random orientations are placed at the bottom of the simulation domain. A thermal gradient of 250 K m1 is applied in the y direction. The orientation and misorientation of a nucleus is defined as the angle of the [10] direction of the square with respect to the x-axis and the direction of thermal gradient, respectively. Fig. 10 a shows the solidification grain structure during the initial phase of grain selection. Misorientation of each grain with respect to the thermal gradient (y-axis) is depicted in the picture. Grains with large misorientations (dark blue and red) cease to grow as soon as the solidification front advances by <0.5 mm. It is clearly seen that the classical selection rules are followed in

the CA model. Fig. 10b shows the final grain structure when the whole domain is solidified. At the end of the solidification only the grains with a misorientation of
±5° survive the grain selection process. This characteristic feature of grain growth during directional solidification matches qualitatively very well with the experimental observations presented in [45,46].

4. Results The proposed CALB model has been applied to study powder based multi-layer hatching. The beam contour is pre-defined for each powder layer. In all the results presented in this work the beam hatches a rectangular cross section. In Fig. 11 the top view of the four variations of beam contours is shown. If the beam is deflected by 90 after every new powder layer then the first four layers traverse the contour presented through (a)–(d) and the same sequence is repeated for higher layers. In all the multi-layer hatching simulations presented in this work the beam is deflected by 180 after every new powder layer implying that it always moves out of the grain structure simulations plane (beam contours presented in Fig. 11 (a) or (c)). A 4.0 mm  0.25 mm baseplate consisting of randomly oriented grains, pre-heated at 900 C is constructed The beam has a power P ¼ 594 W, and line offset of 150 lm, velocity v ¼ 2:2 ms1 and moves out of the simulation plane. The beam movement takes place from [0.75 mm, 3.25 mm] and [1.0 mm, +1.0 mm] in the x and z directions, respectively. The spatial and temporal resolution is chosen as Dx ¼ 0:5 lm and dtLB ¼ 3  109 s, respectively. Pottlacher et al. have reported various temperature dependent thermo-physical material properties of IN718 in [47], these have been used in the present work and some of them are listed in Table 1. The growth kinetics coefficients, i.e., the coefficients relating the local undercooling to the corresponding dendritic tip velocity are derived from the results based on phase field method (PFM). PFM give very accurate description of the local undercooling dependent dendritic tip velocity for rapid solidification, e.g., for pure Ni by Bragard [39]. Since the data for complex alloy systems like IN718 are still not available, in the present work the values for pure Ni are used as the starting point. A polynomial function given by,

v ¼ A  DT þ B  DT 2 þ C  DT 3 ,

is fitted to the reported simulation

data from PFM for pure Ni, where A ¼ 2:2196  102 ms1 K1,

Fig. 9. Schematic of converging (A and B) and diverging (B and C) dendrites at grain boundaries.

B ¼ 1:0533  103 ms1 K2 and C ¼ 1:0302  106 ms1 K3. For IN718 alloy the solidification velocity of v IN718 ¼ 0:1  v pureNi is used in all the simulations because it gives a good qualitative agreement between simulated and experimental grain structure. Fig. 12 shows the layer by layer hatching of the baseplate with powder layers to simulate additive manufacturing. A tool has been developed to calculate the texture (in 2D histograms) and mean

A. Rai et al. / Computational Materials Science 124 (2016) 37–48

45

Fig. 10. Solidification grain structure for directional solidification. (a) Initial grain selection process at the grain boundaries. (b) Final solidification grain structure. The color bar shows the misorientation of each grain with respect to the thermal gradient. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 11. Top view of the possible beam contours for each layer. Line at z ¼ 0 is top view of the grain structure simulation plane, SP and EP are the starting and end points of the beam, respectively.

grain width at a given height. The temperature field calculated using the LB model is transferred to the CA model for grain structure solidification simulation. The resulting phase change information is then transferred back to the LB model. Fig. 13 shows the final grain structure for multi-layer hatching of the baseplate with powder layers to simulate additive manufacturing. The simulated grain structure is compared with the experimental measurements by Helmer et al. [10] for a layer by layer hatching of a 10 mm  10 mm sample produces with same beam

parameters. In experiment the beam is deflected by 90 and in simulation by 180 at each consecutive layer. Columnar grains growing epitaxially and parallel to the build direction are observed both in experiments and simulations. In experimental micrograph stray grains as a result of nucleation are also visible whereas in simulation no new nucleation takes place and the only source of nucleation is partially molten powder particle. A texture parameter s is defined to qualitatively compare the grain texture evolution. If NðyÞ is the number of grains in the

46

A. Rai et al. / Computational Materials Science 124 (2016) 37–48

Table 1 Thermo-physical properties of IN718. Physical property

Value

Unit

Density (liquid) Density (solid) Solidus temperature Liquidus temperature Viscosity

7580 8190 1528 1610

kg m3 kg m3 K K Pa s

Surface tension Specific heat (liquid) Specific heat (solid) Thermal conductivity (solid) at 1000 K Thermal conductivity (solid) at 1500 K Thermal conductivity (liquid) at 1700 K Thermal conductivity (liquid) at 2100 K Latent heat of fusion

4:9  103 1.73 778 652 22.4 30.7 28.1 33.5 227

N m1 J kg1 K1 J kg1 K1 W m1 K1 W m1 K1 W m1 K1 W m1 K1 kJ kg1

bulk at a given height y along the build direction, the corresponding texture parameter is given by

sðyÞ ¼

NðyÞ X li  dhi ;

ð13Þ

i¼1

where li and dhi are the size and absolute misorientation of grain i, respectively. Fig. 14 shows the texture evolution along the build height for the grain structure presented in Fig. 13. It is clearly seen that starting from an equiaxed baseplate (no texture) a strong texture develops due to grain selection with increasing build height. Fig. 15 shows the grain structure and texture for different melting strategies. In Fig. 15a the beam starts at x ¼ 1:5 mm and moves

Fig. 13. Simulated grain structure for multi layer hatching additive manufacturing using CALB. The color bar maps the grain mis-orientation with respect to the build direction (y-axis). Beam parameters: P ¼ 594 W, and loffset ¼ 150 lm and v ¼ 2:2 ms1. The beam movement takes place within [0.75 mm, 3.25 mm] and [1.0 mm, +1.0 mm] in the x and z directions, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

leftwards with a line offset of 150 lm till x ¼ 0:5 mm. Then the new powder layer is spread over the solidified sample and then the beam moves rightwards. This beam sequence is repeated with each new powder layer. The grain structure looks symmetric and the grains from the partially molten powder particles penetrate the sample from both the side (
1 mm). The corresponding texture measure at the build height of 1.7 mm also reflects the symmetric grain structure.

Fig. 12. Snapshot of the temperatue field and the resulting grain structure calculated using CALB model. Color bar: temperature in K (top) and grain misorientation w.r.t the build direction (bottom). The axis denote the length in lattice units (Dx ¼ Dy ¼ 5 lm). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

A. Rai et al. / Computational Materials Science 124 (2016) 37–48

47

leftwards and rightwards when the beam moves from right to left and left to right, respectively. A symmetric grain structure is observed when the scanning directions changes every consecutive layer. 5. Conclusions and outlook

Fig. 14. Texture evolution along the build height for the grain structure presented in Fig. 13.

In Fig. 15b the beam always moves leftwards and the last meltpool in each layer is situated on the left side which has higher temperature than the right side. This temperature difference leads to the surface tension driven convection and the meltpool moves towards the right. The final texture has a peak around
36° (grains on the right side pointing leftwards) The same trend is observed in Fig. 15c where the beam always moves rightwards. This is also reflected in the texture histogram where the grains with misorientation of
28° (grains on the left side pointing rightwards) are dominant. Similar observations have been made in the experiments also [11]. A multi-layer hatching is performed where the scanning direction is kept constant for 10 consecutive layers. The grains are tilted towards the beam scanning directions, i.e.,

The 2D CALB model has been applied to simulate the electron beam melting, powder consolidation and subsequent grain structure solidification. Simulation of the realistic grain texture evolution is a 3D phenomenon, but due to the large computational times required to simulate AM process with individual powder particles in 3D, the 2D CALB model is a valuable tool to study many grain growth related phenomenon, e.g., grain boundary wiggling, grain selection, penetration of new grains from the side walls and layer discontinuities due to partially molten powder particles. The simulated results for single line melting of baseplate have been compared with the experimental data. It has been applied to simulate powder based layer by layer additive manufacturing of IN718 and for the given process parameter columnar grain structure has been observed both in experiments and simulations. The model can be applied to carry out an extensive study to understand the influence of process parameters on the resulting grain structure. The developed tool has been applied to investigate the influence of beam power on stray grain formation due to partially molten powder particles. The solidification velocities used in the simulations is based on the PFM results for pure Ni followed by sensitivity studies to find a good agreement between simulations and experiment for single line melting. Concrete efforts are underway to obtain the solidification velocities for IN718 using PFM and will be used in the future work. In the present work each powder

Fig. 15. Solidified grain structure (top row) and the texture (bottom row) of the sample in the base plate and at the build height of 1.7 mm for a beam power of 594 W, line offset 150 lm and velocity 2.2 ms1. The color bar maps the grain mis-orientation with respect to the build direction (y-axis). The label of each column describes the beam movement. The deviation, dh shows the misorientation of the grain <1 0> axis with respect to the build direction. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

48

A. Rai et al. / Computational Materials Science 124 (2016) 37–48

particle is assumed to contain a single grain in it, whereas, a recent experimental investigation shows the presence of many small grains (
10.6 lm) in a single powder particle. The grain growth algorithm presented in Section 2.2 can be used to simulate the growth of several nuclei (depending on size of the powder particle) distributed within a powder particles under a uniform thermal field. This will enable the inclusion of sub-grain structure within the powder particles. In future this aspect should be taken into account. Since the data for complex alloy systems like IN718 are still not available, in the present work v IN718 ¼ 0:1  v pureNi is used in all the simulations because it gives a good qualitative agreement between simulated and experimental grain structure. In future detailed sensitivity studies will be performed to analyze the influence of this parameter on the final grain structure. Acknowledgments The authors would like to thank the Clean Sky Joint Undertaken under the Grant No. 326020 for the financial support. The authors also gratefully acknowledge funding by the Deutsche Forschungsgemeinschaft (DFG) through the Collaborative Research Center SFB 814, Project B4. References [1] D.C. Hofmann, S. Roberts, R. Otis, J. Kolodziejska, R.P. Dillon, J. Suh, A.A. Shapiro, Z.K. Liu, J.P. Borgonia, Sci. Rep. 4 (2014) 5357. [2] F. Abe, K. Osakada, M. Shiomi, K. Uematsu, M. Matsumoto, J. Mater. Proc. Tech. 111 (1–3) (2001) 210–213. [3] P. Heinl, A. Rottmair, C. Körner, R.F. Singer, Adv. Eng. Mater. 9 (5) (2007) 360– 364. [4] A.A. Antonsamy, J. Meyer, P.B. Prangnell, Mater. Charact. 84 (2013) 153–168. [5] R.R. Dehoff, M.M. Kirka, W.J. Sames, H. Bilheux, A.S. Tremsin, L.E. Lowe, S.S. Babu, Mater. Sci. Technol. 31 (8) (2015) 931–938. [6] Z.H. Liu, D.Q. Zhang, S.L. Sing, C.K. Chua, L.E. Loh, Mater. Character. 94 (2014) 116–125. [7] Q. Jia, D. Gu, J. Alloys and Comp. 585 (2014) 713–721. [8] A. Bauereiß, T. Scharowsky, C. Körner, J. Mater. Proc. Technol. 214 (2014) 2522–2528. [9] D.D. Gu, W. Meiners, K. Wissenbach, R. Poprawe, Int. Mater. Rev. 57 (3) (2012) 133. [10] H.E. Helmer, C. Körner, R.F. Singer, J. Mater. Res. 29 (17) (2014) 1987–1996. [11] C. Körner, H. Helmer, A. Bauereiß, R.F. Singer, MATEC Web of Conferences 14 (08001) (2014) 1–6. [12] C. Körner, A. Bauereiß, E. Attar, Modelling Simul. Mater. Sci. Eng. 21 (2013) 085011. [13] J.D. Williams, C.R. Deckard, Rapid Prototyping J. 4 (1998) 90–100. [14] M. Shiomi, A. Yoshidome, K. Osakada, F. Abe, Int. J. Mach. Tools Manuf. 39 (1999) 237–252.

[15] N.K. Tolochko, M.K. Arshinov, A.V. Gusarov, V.I. Titov, T. Laoui, L. Froyen, Rapid Prototyping J. 9 (2003) 314. [16] K. Dai, L. Shaw, Acta Mater. 52 (2004) 69–80. [17] M.F. Zäh, S. Lutzmann, Prod. Eng. Res. Devel. 4 (2010) 15–23. [18] C. Körner, E. Attar, P. Heinl, J. Mater. Proc. Technol. 211 (2011) 978–987. [19] A. Klassen, A. Bauereiß, C. Körner, J. Phys. D: Appl. Phys. 47 (2014) 065307. [20] R. Ammer, U. Rüde, M. Markl, V. Jüchter, C. Körner, Int. J. Mod. Phys. C 25 (12) (2014) 1441009. [21] T. Carozzani, H. Digonnet, Ch.-A. Gandin, Modelling Simul. Mater. Sci. Eng. 20 (2012) 015010. [22] T. Carozzani, Ch.-A. Gandin, H. Digonnet, Modelling Simul. Mater. Sci. Eng. 22 (2014) 015012. [23] G. Guillemot, O. Desmaison, S. Chen, M. Bellet, Ch.-A. Gandin, A multi-physic CAFE approach for the simulation of grain structure development in GMAW processes, in: C. Sommitsch, N. Enzinger, P. Mayr, H. Cerjak (Eds.), 11th Int. Seminar on Numerical Analysis of Weldability, Seggau, Austria, 27–30 September 2015, Proc. in Mathematical Modelling of Weld Phenomena 11, The Institute for Materials Science and Welding (IWS) at Graz University of Technology (TU Graz), 2015, p. 16. [24] Z. Fan, T. Sparks, E. Todd, F.W. Liou, A.B. Jambunathan, Y. Bao, J. Ruan, J.W. Newkirk, Faculty Research and Creative Works (2007) 11133. [25] J. Zhang, F.W. Liou, W. Seufzer, J.W. Newkirk, Z. Fan, H. Lui, T.E. Sparks, Faculty Research and Creative Works (2013) 11226. [26] Ch.-A. Gandin, M. Rappaz, Acta mater. 45 (5) (1997) 2187–2195. [27] T.M. Pollock, S. Tin, J. Propul. Power 22 (2006) 361–374. [28] Z. Hu, L. Liu, T. Jin, X. Sun, Aviation Engine 31 (2005) 1–7. [29] E. Attar, C. Körner, J. Colloid and Interf. Sci. 335 (2009) 84–93. [30] P.D. Lee, P.N. Quested, M. Mclean, Phil. Trans. R. Soc. Lond. A 356 (1998) 1027– 1043. [31] K.C. Mills, Recommended Values of Thermophysical Properties for Selected Commercial Alloys, Woodhead Publishing Limited, Cambridge CB1 6AH, England, 2002. [32] A. Klassen, T. Scharowsky, C. Körner, J. Phys. D: Appl. Phys. 47 (2014) 275303. [33] D.E. Ovsienko, G.A. Alfintsev, V.V. Maslov, J. Crystal Growth 26 (1974) 233– 238. [34] T.P. Schulze, Phys. Rev. E 78 (2008) 020601. [35] A.F. Ferreira, L.O. Ferreira, A.C. Assis, J. Braz. Soc. Mech. Sci. Eng. XXXIII (2) (2010) 125–130. [36] G.I. Tóth, G. Tegze, T. Pusztai, G. Tóth, L. Gránásy, J. Phys. Condens. Matter 22 (2010) 364101. [37] W. Kurz, B. Giovanola, R. Trivedi, Acta Metall. 34 (5) (1986) 823–830. [38] M. Rappaz, Ch.-A. Gandin, Acta Metall. Mater. 41 (2) (1993) 345–360. [39] J. Bragard, A. Karma, Y.H. Lee, Interface Sci. 10 (2002) 121–136. [40] Ch.-A. Gandin, M. Rappaz, R. Tintillier, Metall. Trans. A 24A (1993) 467–479. [41] Ch.-A. Gandin, J.-L. Desbiolles, M. Rappaz, Ph. Thévoz, Metall. Mater. Trans. A 30A (1999) 3153–3165. [42] P. Meakin, R. Juillien, J. Phys. 48 (1987) 1651–1662. [43] Ch.-A. Gandin, R.J. Schaefer, M. Rappaz, Acta Mater. 44 (8) (1996) 3339–3347. [44] D. Walton, B. Chalmers, Trans. Metall. Soc. AIME 215 (1959) 447. [45] D. Szeliga, K. Kubiak, A. Burbelko, M. Motyka, J. Sieniawski, J. Mater. Engg. Perfor. 23 (3) (2014) 1088–1095. [46] C.-A. Gandin, J.-L. Desbiolles, M. Rappaz, M. Swierkosz, P. Thévoz, Supercomput. Rev. (8) (1996) 11–15. [47] G. Pottlacher, H. Hosaeus, E. Kaschnitz, A. Seifter, Scandinavian J. Metall. 31 (2002) 161–168.