close

Вход

Забыли?

вход по аккаунту

?

j.cirpj.2017.09.007

код для вставкиСкачать
G Model
CIRPJ 441 No. of Pages 18
CIRP Journal of Manufacturing Science and Technology xxx (2017) xxx–xxx
Contents lists available at ScienceDirect
CIRP Journal of Manufacturing Science and Technology
journal homepage: www.elsevier.com/locate/cirpj
On thermal modeling of Additive Manufacturing processes
Panagis Foteinopoulos, Alexios Papacharalampopoulos, Panagiotis Stavropoulos*
Laboratory for Manufacturing Systems and Automation, Department of Mechanical Engineering and Aeronautics, University of Patras, Patras 26504, Greece
A R T I C L E I N F O
Article history:
Available online xxx
Keywords:
Additive Manufacturing
Modeling
Simulation
Thermal
Finite Differences
Adaptive Mesh
A B S T R A C T
A two-dimensional Finite Difference (FD) model of the thermal history of parts manufactured in powder
bed fusion Additive Manufacturing (AM) processes is presented. The temperature of the part is calculated
in each time-step taking into account the moving laser heat source, the melting phase change and
functions of both temperature and porosity are used for the material thermal properties. Also, an
algorithm for node birth and distance adaptation over time is utilized, minimizing computational time
and memory. A validation of the results of the model is included.
© 2017 The Author(s). This is an open access article under the CC BY-NC-ND license (http://
creativecommons.org/licenses/by-nc-nd/4.0/).
Introduction
Additive Manufacturing (AM) is defined as “the process of
joining materials to build objects from 3D model data, usually layer
upon layer, as opposed to subtractive manufacturing processes,
such as traditional machining” [1]. The main difference between
Rapid Prototyping is that AM specifically aims to the manufacturing of end user parts rather than just prototypes [2]. AM processes
have more than 25 years of history [3] and the interest in them is
steadily increasing due to the design freedom [4], the potential of
producing near net shape structural components and the
environmental and ecological promise they offer.
In most of the AM processes, parts are manufactured layer by
layer, using a source of thermal energy to fuse the different layers
together. As a result, anisotropic material properties and residual
stresses are common, because of the non-homogenous thermal
and cooling phenomena that take place [5]. Except from the
uncertainty of the mechanical properties, other important issues
are the low productivity and poor surface quality, the optimization
of which is difficult because of limited modeling approaches to the
topic [6].
The residual stresses and distortions, which are caused by the
non-homogenous thermal phenomena (heating and cooling) [7]
that take place in AM, deteriorate the mechanical properties and
the dimensional accuracy of the parts. As a result, the thermal
history of a part’s manufacturing procedure is essential, because it
determines its microstructure, mechanical properties and final
dimensions. To this effect, the thermal modeling of the AM
processes can be utilized for the optimization of those important
Key Performance Indicators (KPIs), without the requirement for
time-consuming and costly experiments and it is the first step to
establishing relationships between the KPIs or Quality Performance Measures (QPM) of a part and the variables of the process
(Fig. 1).
There are many different approaches to the modeling of the
thermal history of parts, manufactured by AM processes. Most of
the existing studies utilize numerical methods, due to the
complexity of the phenomena that take place. More specifically,
the modeling of heat transfer of AM metal deposition, via Finite
Elements (FE), takes place in the work of Ref. [8], along with an
error minimization. A temperature field simulation of the Selective
Laser Melting (SLM) process, also by using the FE method, is
presented in Ref. [9]; the same numerical method is utilized by Ref.
[10] for the simulation of the temperature distribution and the
melt pool size, when the bulk of powder is heated by a laser source.
The FE method has also been used by Refs. [10–14]. Different
modeling methods have been followed by some studies, like that of
Ref. [15], in which a computational tool has been developed by
assembling models of many interacting particles in the small scale.
Also, the laser energy was correlated to the Total Area of Sintering
(TAS) via a convex hull based approach by Ref. [16]. Heat transfer
modeling for the SLM process has been carried out via discrete grid
models [17], which take porosity into account. In the works of Refs.
[18,19], the finite volume method has been used for the thermal
* Corresponding author.
E-mail address: pstavr@lms.mech.upatras.gr (P. Stavropoulos).
https://doi.org/10.1016/j.cirpj.2017.09.007
1755-5817/© 2017 The Author(s). This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
Please cite this article in press as: P. Foteinopoulos, et al., On thermal modeling of Additive Manufacturing processes, NULL (2017), https://doi.
org/10.1016/j.cirpj.2017.09.007
G Model
CIRPJ 441 No. of Pages 18
2
P. Foteinopoulos et al. / NULL xxx (2017) xxx–xxx
Nomenclature
T2
Tb
Dt
Dx
Dz
Dza
p
r
ra
rm
rpr
A
c
ca
cap
cl
cm
cn
cpr
cs
ct
Dl
d
f
h
I
I0
k
km
kpr
L
LL
Ln
Lo
ld
m
NL
Nl
Nt
Nx
Nz
n
nh
ns
P
r
T
T1
Time duration of a time-step
Length between two nodes in the x-axis
Length between two nodes in the z-axis (standard, nonadaptive meshing)
Array storing the distances between the adaptive nodes
of the mesh
Pi constant
Material density
Density of air
Density of fully dense solid material
Density of porous solid material
Apparent heat capacity coefficient
Heat capacity of the material
Heat capacity of air
Apparent heat capacity of the material (function of
temperature)
Heat capacity of the liquid phase of the material
(function of temperature)
Heat capacity of fully dense solid material
Heat capacity of the solid and liquid phases of the
material (function of temperature)
Heat capacity of porous solid material
Heat capacity of the solid phase (function of temperature)
Total heat capacity (solid, liquid phases and apparent
heat capacity) of the material (function of temperature)
Diameter of the laser spot
Adaptive mesh node distance non-uniformity exponent
Volumetric fraction of porosity
Convective heat transfer coefficient
Laser beam intensity
Laser beam intensity at the beam axis and at the focal
level
Thermal conductivity of the material
Thermal conductivity of fully dense solid material
Thermal conductivity of porous solid material
Latent heat of the material
Layer thickness
The length, taking into account the addition of a new
layer, in the z direction in which the meshing will be
created using the adaptive algorithm
The length, without taking into account the addition of
a new layer, in the z direction in which the meshing will
be created using the adaptive algorithm
Distance from the laser beam axis
The x coordinate of a node
Number of nodes in the z direction in the thickness of a
layer of a part
Number of nodes in the x direction in the length of the
diameter of the laser spot
Total number of time-steps
Total number nodes in the x-axis
Total number nodes in the z-axis
The z coordinate of a node
Number of time-steps a node is being heated
Number of time-steps needed for the addition of a new
layer of powder
Laser power
Radius of the laser beam spot
Nodal temperature
Lower temperature boundary of the mushy area of the
apparent heat capacity method
T env
Tm
T pre
t
th
ts
vl
x; z
zn
zo
Upper temperature boundary of the mushy area of the
apparent heat capacity method
Temperature of the building platform of the AM
machine
Environmental Temperature
Melting temperature of the material
Temperature of the new layers of powder that are
added over time
Time
Time a node is being heated
Time needed for the addition of a new layer of powder
Laser head scan speed
Cartesian coordinates
Final positions of the nodes of the adaptive mesh in the
z-axis
Previous positions of the nodes of the adaptive mesh in
the z-axis
modeling of the SLM process; in the latter, the densification of the
material (WC/CU composite powder) and the induced surface
tensions are also simulated. A different approach has been
followed by Ref. [20], in which the OpenFOAM software has been
utilized for the modeling of the process dynamics of the laser beam
melting AM process.
However, due to the speed of the process and the high
complexity of the spatial and temporal dynamic thermo-mechanical phenomena that take place, the computational cost, time and
memory needed for the numerical modeling of AM processes tends
to be very high, especially when combined with the need of the
simulation of the entire thermal history [7]. As a result, most of the
models simulate only a short time-span of the manufacturing of a
part and not the whole process. However, such approaches are
unable to provide the necessary information for the calculation of
the thermal induced stress fields and deformations, because the
entire thermal history, including the cooling down rates, is
necessary for this. It has to be pointed out that such information
is very important for the design and manufacturing engineers, in
order to take the necessary actions, like changes in the design that
will enable a more homogenous cooling, creation of supports that
will minimize the distortions and simultaneously offer force
cooling, or change the process parameters in a way that will
minimize or even prevent such unnecessary phenomena (thermal
distortions, non-homogenous mechanical properties).
Addressing the gap in the existing state of the art, this study
emphasizes in the creation of a practical and fast to run, yet
accurate in its predictions, model of the thermal history of a part
which is manufactured in a powder bed fusion AM process. This
model’s simulation was not created through a ready to use
software, but it was custom made instead, so as to be tailored to the
complex and dynamic problem at hand. This decision was made,
because such a solution provides better adaptability, easier
coupling with other fields (e.g. mechanical) and offers increased
connectivity with other modules, such as optimizers. The FD
method has been used in this study because of advantages such as
strict formulation and ease concerning user inputs [21]. In order to
keep the computational time and cost, as well as the accuracy loss,
to a minimum, a two-dimensional (2D) space combined with a
non-uniform mesh has been used, which is dense in the regions
where complex dynamic phenomena take place, while it becomes
coarser in places that less dynamic temporal and spatial changes
occur. A further increase of the accuracy is achieved by assuming
temperature dependent material thermal properties; namely
thermal conductivity, specific heat capacity and density. In
Please cite this article in press as: P. Foteinopoulos, et al., On thermal modeling of Additive Manufacturing processes, NULL (2017), https://doi.
org/10.1016/j.cirpj.2017.09.007
G Model
CIRPJ 441 No. of Pages 18
P. Foteinopoulos et al. / NULL xxx (2017) xxx–xxx
3
Fig. 1. Modeling in Additive Manufacturing: part quality, cost and production rate.
addition, different thermal properties are used for the unmelted
powder and the fully dense part (porosity of the powder is taken
into account). Moreover, the latent heat of the solid–liquid phase
change is also taken into consideration, utilizing the apparent heat
capacity method [22]. The combination of the aforementioned
actions lead to a highly integrable model, which can also be used
for machine control, because it keeps simulation time, memory
and cost to a minimum, while maintaining highly accurate results,
which have also been verified with the use of data from Ref. [9]. In
this study, the modeling of the first level KPIs, as presented in Fig. 1,
takes place, whereas, the remaining aspects will be the subject of
future work.
Power Bed Fusion (PBF) is one of the biggest groups of AM
processes, encompassing Selective Laser Sintering (SLS), Selective
Laser Melting (SLM) and Electron Beam Melting (EBM) [23]. In the
PBF processes, fine powder, which has been spread on a building
table, is heated by a laser beam (an electron beam in EBM) so as to
allow the grains to fuse together. Different powder consolidation
mechanisms are used in each process [24] and the unused material
powder may be cleaned away and recycled [25,26]. The heat
energy absorbed by the top layers is conducted to the rest of the
part’s mass, mainly perpendicularly to the heated surface, starting
from the laser heated top, the temperature of which is the highest,
to the bottom of the part that is in contact with the machine table.
Furthermore, conduction occurs in the horizontal axis to the
surface; however, the significance of this phenomenon diminishes
in the lower region of the part. Thermal losses occur in the form of
radiation or convection to the air and conduction to the machine
table. The model that has been developed in this study has been
calibrated for the full melting of the powder particles, which is a
characteristic of the SLM process [27], but it can be easily adapted
to other PBF processes as well.
In the following sections, the problem and the approach that is
followed are defined and the equations that describe the
phenomenon in hand are stated, both in their analytical and
discretized form, as well as the assumptions that were made. There
is a description of the meshing (uniform and non-uniform), as well
as the apparent heat capacity method (simulation of melting) and
the temperature/powder-density material functions of the thermal
properties. Finally, the discussion of the results and the conclusion
are presented and in Appendix sections a comparison is presented
between the adaptive meshing algorithm and the non-adaptive
conventional one concerning the accuracy of the results and the
required calculation time.
Modeling approach
The thermal model developed in this study calculates the
temperature of a part under construction over time. A 2D space has
been used for the modeling of the part. Moreover, the material
thermal properties (conductivity, heat capacity) and density are
functions of temperature and not constant values. In addition, the
difference in the material properties, due to the porosity of the
unheated powder, has been taken into account. The thermal
analysis starts at the point when the first layer of powder is laid
upon the machine table and the laser head begins to offer energy to
the powder. The analysis continues over time, by simulating the
movement of the laser scanner head. This is achieved by moving
the heating boundary condition across the top layer over time.
When the laser head has completed the scanning of the entire top
layer the addition of an extra layer of powder is simulated, as well
as the time required for the roller to spread it. Then, the reverse
movement of the moving heat source is modelled. In this way, the
creation of a wall of material (2D part) is simulated. A reduction of
all the used units has been made in order for them to correspond to
the 2D approach. The nodal temperature is calculated in each time
step and it is stored, calculating in this way the temperature history
of the part. The model uses as input the part dimensions and the
material, laser and machine properties, as well as certain
environmental conditions data. A flowchart of the modeling
approach can be seen in Fig. 2.
The differential equation of heat conduction through an
isotropic material, in 2D Cartesian coordinates, has been used:
r2 kT ¼ rc
@T
@t
@2 T
@2 T
@T
þ k 2 ¼ rc
2
@t
@x
@z
k
ð1Þ
ð2Þ
Please cite this article in press as: P. Foteinopoulos, et al., On thermal modeling of Additive Manufacturing processes, NULL (2017), https://doi.
org/10.1016/j.cirpj.2017.09.007
G Model
CIRPJ 441 No. of Pages 18
4
P. Foteinopoulos et al. / NULL xxx (2017) xxx–xxx
Fig. 2. General flowchart of the modeling approach.
The assumptions that have been made in this study are the
following:
The part created is situated in a 2D space and its shape is
rectangular and continuous.
Since this study focuses on the Additive Manufacturing of metals,
which are heat conductors, the heat transfer through radiation is
negligible in comparison to that due to conduction.
The thermal conductivity, specific heat capacity and density
properties that have been used in this study, are temperature
dependent. Any other material properties used are constant and
not temperature dependent.
The metal’s vaporization is not taken into account, as it is
negligible in the processes at hand.
The machine building platform is considered as a heat sink (its
heat capacity is infinite) and as a result, its temperature remains
constant (equal to Tb).
Each new layer of powder is spread in one iteration, and has a
temperature of Tpre; its spreading is simulated by introducing a
time delay to the heating source.
Changes in dimensions due to temperature induced differences
in density or the powder/solid state of the material, have been
neglected.
Please cite this article in press as: P. Foteinopoulos, et al., On thermal modeling of Additive Manufacturing processes, NULL (2017), https://doi.
org/10.1016/j.cirpj.2017.09.007
G Model
CIRPJ 441 No. of Pages 18
P. Foteinopoulos et al. / NULL xxx (2017) xxx–xxx
The boundary condition used for each category of boundary
nodes is presented here in detail, where m, n are the x and z
coordinates respectively of the corresponding nodes.
Top
i. Heating
The following equation is used for heating:
@T
¼I
j
@z x ¼ m
ð3Þ
k
z¼n
A Gaussian profile is assumed for the laser beam intensity I,
which is described as a function of distance from the laser beam
axis by:
l 2
d
ð4Þ
I ¼ I0 e2 r
where, ld is the distance, in the x-axis, of the node from the laser
beam axis, r is the radius of the laser beam and I0 is the intensity of
the laser beam at the beam axis and at the focal level [21] and it is
specified by the laser power P and the laser spot radius r as:
I0 ¼
2P
ð5Þ
pr 2
However, since the case studies in this work have been
conducted in a 2D Cartesian space, setting the quantity @@yT2 equal
to zero (as aforementioned in Section “Modeling approach”), a
correction factor for power has been introduced; t
2
ii. Convection
Thermal losses to the environment, from the top side of the part
due to convection are modeled with the following equation
(Neumann boundary condition):
@T
¼ h½T env T ð0Þ
j
@z x ¼ m
k
ð6Þ
z¼n
where, T ð0Þ is the temperature of the boundary and h is the
convective heat transfer coefficient. A conservative value of
W
h ¼ 10 =ðm2 KÞ has been used for the convective heat transfer
Bottom
Considering the bed as a heat sink of a constant temperature
equal to Tb, the thermal losses of this row of elements were
modeled using the convection equation. An experimentally
calculated heat transfer coefficient, with a higher value, has been
used for the modeling of the higher rate of energy exchange that
takes place between the part and the bed, in comparison to that
between the part and the air.
Eq. (2) has been discretized using central differences and it is
applied in the stencil nodes (nodes that are not situated on the
edges of the grid), when a uniform mesh is used.
The boundary conditions have also been discretized with the
use of either forward or backward differences, according to their
placement in the grid. The nodes, which the different boundary
conditions are applied to, can be seen in Fig. 3.
The implicit method has been used for the solution of the
system: rearrangement of the terms leads to a linear algebraic
system, the solution of which provides the nodal temperatures at
each time step. Then, the time step progresses by one and the same
procedure is repeated until the entire part is constructed or the
user specified maximum calculation time has been reached. The
graphical representation of the algorithm in two consecutive timesteps is depicted in Fig. 4.
As far as mesh creation is concerned, the selection of the
discretization parameters, for the constant part of the mesh, is
conducted by the following procedure: Each new layer of powder is
modeled using NL layers of nodes, with NL 3. In this way, at any
time at least one bottom, one stencil and one top node exist,
rendering the aforementioned calculations possible. As a result, Dz
is defined by:
Dz ¼
LL
NL 1
Sides (left–right)
Dx is calculated in order for a Nl number of nodes to be included
in the length of the spot diameter. The number of nodes Nl has been
defined so as to ensure accuracy of the results as well as
minimization of the computation time and cost. The following
equation provides Dx:
Dx ¼
Dl
Nl 1
Dl
vl
The convection equation that was previously described, has also
been used here, where m, n are the coordinates of the
corresponding nodes.
th ¼
Corners
nh ¼
The mesh of the part has four corners at any given moment. The
temperature of the corner nodes is affected by their adjacent
nodes, but it does not affect them in turn, since the boundary nodes
interact only in one dimension with other nodes. As a result, an
insulation boundary condition, with a vertical interaction, has been
chosen and after the solution has been obtained, the temperature
of those nodes is replaced by the mean value of the temperatures of
their two adjacent ones. The equation used for the corner nodes is
the Neumann boundary condition with a heat flux of zero:
As a result, Dt is:
@T
¼0
@x
ð8Þ
ð9Þ
The heating boundary condition is applied to Nl nodes that are
under the laser spot at any given time-step. The time, in seconds,
those nodes are heated is determined by:
coefficient, according to Ref. [28].
k
5
ð7Þ
ð10Þ
This time is converted to a number of time-steps by:
Dt ¼
th
Dt
Dl
ðnh 1Þvl
ð11Þ
ð12Þ
The values of NL, Nl and nh have been defined in order for the
convergence of the results to be ensured, while simultaneously
keeping the computing time and cost to a minimum. The
convergence analysis can be found in detail in a later chapter of
this study.
The mesh of the part is created using the above mentioned
discretization. However, neither the part’s dimensions nor the region
where the laser heating is applied to are constant over time.
Subsequently, both the sum of nodes and the nodes to which the
Please cite this article in press as: P. Foteinopoulos, et al., On thermal modeling of Additive Manufacturing processes, NULL (2017), https://doi.
org/10.1016/j.cirpj.2017.09.007
G Model
CIRPJ 441 No. of Pages 18
6
P. Foteinopoulos et al. / NULL xxx (2017) xxx–xxx
Fig. 3. Constant mesh of two layers. The boundary conditions can also be seen. Heat conduction takes place in the stencil nodes (big black circles).
boundary conditions are applied are not constant over time. More
specifically, the convection boundary condition is applied to all the
top nodes and in some of them the heating convection is also applied.
The group of the top nodes, to which the heating boundary
condition is applied, depends on the time-step of the simulation.
As a result, at any given time step, the top layer of nodes is divided
into two groups: (a) those that are heated (heating and convection
boundary condition) and (b) those the temperature of which drops
because of convective heat transfer to the environmental air
(convection boundary condition). The laser spot travels over time
with a speed equal to that of the laser scan head and this is
simulated by changing accordingly the group to which the top
nodes belong. This procedure is repeated until the laser scanner
head has completed traversing the full length of the part and as a
result, all the top layer nodes have been heated. Upon the
completion of the heating of all the top nodes, a new set of NL layers
of nodes are added above the previous top one, in order to simulate
the new layer of powder.
A delay of ns time-steps is introduced in order to simulate the
time needed for the spreading of the powder, for the duration of
which no nodes are heated and the convection boundary condition
is applied to all the top layer nodes. The number of time-steps, ns,
needed for the spreading process is calculated by:
ns ¼
ts
Dt
ð13Þ
where, ts is the amount of time required for the machine to lay out
the new layer of powder, which is user defined. After the ns timesteps have passed, the heating boundary condition is applied to the
first group of Nl nodes of the new layer. Those nodes are situated in
the left or right end of the top side (new layer), according to the
direction the laser scanner head had been moving during the
heating of the previous layer. More specifically, when simulation
starts, the laser scanner head starts moving from the left to the
right of the top side, heating the corresponding nodes as it moves.
Upon the completion of the heating of the whole length of the top
side, the laser scanner head will be situated on the right end of the
top side of the part. When the new layer is added, the first group of
heated nodes will be a number of Nl nodes situated on the right end
of the top side of the new layer of the part. As a result, the scanner’s
speed direction changes after each layer has been completed. The
above procedure is repeated as many times as it is required for the
entire part’s construction to be simulated or up to when the
maximum user specified time-step is reached.
Please cite this article in press as: P. Foteinopoulos, et al., On thermal modeling of Additive Manufacturing processes, NULL (2017), https://doi.
org/10.1016/j.cirpj.2017.09.007
G Model
CIRPJ 441 No. of Pages 18
P. Foteinopoulos et al. / NULL xxx (2017) xxx–xxx
7
Fig. 4. Graphical representation of the algorithm. The temperatures of the stencil nodes of the “t+1” time-step are calculated using as input the temperature of the highlighted
node of the “t” time-step.
The selection of Dx and Dz is made so as to ensure the needed
mesh density for (a) the accurate description of the dynamic heat
transfer phenomena that take place and (b) the convergence and
stability of the numerical method. However, as more layers are
created, the z-dimension of the part becomes larger and many of
the nodes are now situated far from the heat source. However, the
heat transfer phenomena that take place in the lower part of the zaxis are far less dynamic than those taking place close to the heated
surface. As a result, the mesh density, which was designed for its
capacity to describe the highly dynamic phenomena, which now
take place only in the upper part of the mesh, is no longer
necessary. In fact, the use of a mesh with sparser nodes in the lower
part, would be able to accurately describe the phenomenon, while
simultaneously reducing the computational time and cost.
This was the reason why an adaptive meshing algorithm has
been developed. The main idea of this algorithm is that the
addition of new layers of material is described by adding new
layers of nodes up to the point where a certain number of layers of
nodes has been reached. After that, any addition of new layers of
material is simulated by increasing the distance between the
already existing nodes instead of adding more layers of nodes.
More specifically, each time the addition of a new layer of material
is to be simulated, the distance in the z direction of the nodes that
are situated in the lower section of the part, where the thermal
phenomena are less dynamic in space and time, is increased,
according to the following equation:
Ln ¼ Lo þ LL
ð14Þ
In this way, the length of the mesh is increased in the z direction
proportionally to the length of a new layer. Consequently, the new
layers of the part are simulated without adding any more layers of
nodes. The resulting mesh is non-uniform: dense in the upper part
of z-axis (uniform part of the mesh, in which the original Dx and Dz
are used), while in the lower part of the z-axis it gets coarser (nonuniform part of the mesh). The lower, non-uniform part of the
mesh will be referred to as adaptive mesh in the rest of the study.
The next step is the distribution of the nodes of the adaptive
mesh in the z direction of the Ln length. The new positions zn in the
z direction of the nodes of the adaptive mesh are calculated by the
following equation:
zn ¼
d
zo
Ln
ð15Þ
where zo is the vector containing the old positions of the nodes and
the exponent d is a constant, the use of which leads to non-uniform
distances between the nodes of the adaptive mesh. In this way the
distances between the nodes are increased exponentially (the
closer to the bottom of the part, the coarser the mesh); however,
the d coefficient is chosen so as to ensure that the adaptive mesh
density decreases gradually, thus minimizing the possibility for
stability issues of the numerical method. In Fig. 5 a schematic of
the adaptive mesh is depicted, in which NL = 3 has been used. The
temperature values of the nodes of the adaptive mesh in their new
zn positions are calculated utilizing the known temperature nodal
field and a spline regression algorithm (more details in Appendix
section b).
The next step is the calculation of the distances in the z
direction between the nodes of the adaptive mesh. The values of
those distances, which are different between each pair of nodes,
are stored into the Dza array. Those values are used in replacement
of the Dz of the uniform mesh in the FD equation for the nodes of
the adaptive mesh. As a result, a new FD equation had to be derived
by utilizing the truncated Taylor series formula, adopted for two
different mesh sizes. The FD equation is the following:
tþ1
tþ1
T tþ1
mþ1;n 2T m;n þ T m1;n
k
2
Dx
tþ1
tþ1
Dz1 T tþ1
m;nþ1 2 Dz1 þ Dz2 T m;n þ Dz2 T m;n1
þ 2k
2
2
Dz 1 Dz 2 þ Dz 1 Dz 2
t
T tþ1
m;n T m;n
¼ rc
dt
ð16Þ
The equation is solved using the same implicit scheme that has
been previously described. The flowchart for the procedure of the
adaptive mesh’s creation can be seen in Fig. 6. The thermal
properties of a material depend both on its temperature and, as far
as AM is concerned, on whether the material is in powder form or
solidified. As it was previously stated, the material properties
Please cite this article in press as: P. Foteinopoulos, et al., On thermal modeling of Additive Manufacturing processes, NULL (2017), https://doi.
org/10.1016/j.cirpj.2017.09.007
G Model
CIRPJ 441 No. of Pages 18
8
P. Foteinopoulos et al. / NULL xxx (2017) xxx–xxx
Fig. 5. Adaptive meshing.
(namely thermal conductivity, specific heat capacity and density)
that have been used in this model are functions of temperature and
porosity. The thermal conductivity and density of the porous
material are calculated by the Maxwell model [29]:
kpr
km
;0 f < 1
¼
1 0:75f
ð17Þ
where, f is the volumetric fraction of porosity, with f = 0 referring to
a material with zero porosity. Concerning the change of density,
due to porosity, the following equation, from Ref. [30], is used:
rpr ¼ ð1 f Þrm þ f ra
ð18Þ
Due to a difference of three orders of magnitude between ra and
rm, the previous equation can be written as:
rpr ¼ ð1 f Þrm
ð19Þ
The above equations can be used for the calculation of the thermal
conductivity and the porosity of the powder by using as inputs the
volumetric fraction of porosity and the properties of the fully dense
material. The same type of equation is used for the calculation of
the specific heat of the porous medium:
cpr ¼ ð1 f Þcm þ f ca
ð20Þ
However, the term ca cannot be ignored here, as the specific heat of
air is in the same order of magnitude as that of aluminum. In Fig. 7,
the function of the specific heat capacity of air versus temperature,
which has been used for the calculation of the specific heat of the
porous medium, can be seen.
The thermal properties of the material also depend on its phase.
Experimental data [32,33,34] have been utilized for the creation of
analytical functions of property-temperature, using polynomials,
for both the solid and liquid phase of each property. Two different
polynomials have been created for each property: one for the solid
and one for the liquid phase of the material. As a result, there is a
loss in the continuity of the property-temperature function when
the temperature is equal to that of the phase change. A third
function has been created for each property, leading to the smooth
transition (no continuity loss) between the liquid and the solid
property-temperature functions, ensuring their numerical stability and solution convergence. This function is used in a temperature
Please cite this article in press as: P. Foteinopoulos, et al., On thermal modeling of Additive Manufacturing processes, NULL (2017), https://doi.
org/10.1016/j.cirpj.2017.09.007
G Model
CIRPJ 441 No. of Pages 18
P. Foteinopoulos et al. / NULL xxx (2017) xxx–xxx
9
Fig. 6. Flowchart of the adaptive-mesh algorithm.
range around the melting point and it is defined by a third grade
polynomial, which ensures the same value and slope in its two
boundaries. In Figs. 8–10 there is a depiction of the calculated
functions of the thermal properties of aluminum over the
temperature, which takes into account the phase changes as well.
In the case of heat capacity, without the addendum of the apparent
heat capacity method that is later described, the following function
has been used:
c cs
T Tm
1 þ tanh 8
ð21Þ
cn ¼ cs þ l
2
T2 T1
Fig. 8. Thermal conductivity of aluminum as a function of temperature and phase
(solid–liquid). Experimental data from Ref. [33].
where, cs and cl are the functions of heat capacity of the solid and
liquid phase that have been created using experimental data from.
This selection for the cn function ensures that it maintains its
continuity and as result, it will not disrupt the stability and
convergence of the numerical scheme. The “mushy area”
temperature range is between T1 and T2. In Fig. 10, the plot of cn
over temperature can be seen.
The melting/solidification phase change is modeled using the
apparent heat capacity method [22], according to which, specific
heat is artificially increased in a small temperature range around
the melting point, so as to encompass the latent heat of the
melting/solidification phase change [22]. In this way, the material’s
resistance to the changing of its temperature is increased and as a
result, more energy is required for a temperature change to occur,
encompassing the latent heat of melting. This artificially increased
heat capacity is called apparent heat capacity and the temperature
range that is used is technically referred to as “mushy area”. This
term is not used to indicate the existence of an actual mushy zone,
such as that of alloys, but as an assumption of the apparent heat
capacity method. The apparent heat capacity can be calculated to
either be a constant or a function of temperature. In either case, the
area encompassed below the graph of the apparent heat capacity
versus temperature in the “mushy area” temperature range, has to
be equal to the material’s latent heat. More specifically, the value of
the apparent heat capacity is defined by the necessary energy in
order for the phase change to take place and it is calculated by the
integration of the function that is used for the apparent heat
Fig. 7. Specific heat capacity of air as a function of temperature. Experimental data
from Ref. [31].
Fig. 9. Density of aluminum as a function of temperature and phase (solid–liquid).
Experimental data from Ref. [34].
Please cite this article in press as: P. Foteinopoulos, et al., On thermal modeling of Additive Manufacturing processes, NULL (2017), https://doi.
org/10.1016/j.cirpj.2017.09.007
G Model
CIRPJ 441 No. of Pages 18
10
P. Foteinopoulos et al. / NULL xxx (2017) xxx–xxx
Fig. 10. Specific heat capacity of aluminum as a function of temperature and phase.
Data from Ref. [32].
capacity over the temperature range of the “mushy area”. The
function used in this study for the apparent heat capacity is the
squared cosine with an argument of half a period:
T Tm
ð22Þ
cap ¼ Acos2 P
T2 T1
where, A is calculated by the following equation:
ZT 2
L¼
ZT 2
cap dT ¼
T1
ZT 2
¼A
T1
cos2 p
T1
T Tm
Acos2 p
dT
T2 T1
T Tm
AðT 2 T 1 Þ
dT ¼
2
T2 T1
ð23Þ
Solving for A:
A¼
2L
T2 T1
ð24Þ
Subsequently, the sum of cn and cap constitutes the heat
capacity function that is used in the entire temperature range:
ct ¼ cn þ cap
ð25Þ
In Fig. 11 can be seen the total specific heat capacity over
temperature, in which both the standard heat capacity of the
material and the apparent heat capacity are included. The
artificially enhanced value, which simulates the latent heat of
the phase change, can be seen in the “mushy area” temperature
range in red color.
Implementation
A convergence analysis has been conducted in order to
determine the values of NL, Nl and nh. The part dimensions used
for this analysis, as well as the laser and model parameters, can be
seen in Table 1. The central nodes of the x-axis of the top four part
layers (the first, in z-axis, node in each layer) of the part will be
used for the test at hand. With reference to the first convergence
analysis, for the determination of NL, three different simulations
have been conducted: one using NL = 6, another one using NL = 9
and one using NL = 12, the results of which can be seen in Fig. 12–14
respectively. It can be observed that the difference between the
results of the simulations that use NL = 6 and NL = 9 is small enough
Fig. 11. Specific heat capacity of aluminum as a function of temperature and phase.
Data used from Ref. [32]. Use of the apparent heat capacity method for the
simulation of the latent heat of the phase change. (For interpretation of the
references to color in the text, the reader is referred to the web version of this
article.)
in order to assume that convergence has been achieved when
NL = 6. However, NL = 9 is used in order to ensure the convergence of
the numerical scheme when the adaptive meshing strategy is used,
in which the distance between the nodes is increased to simulate
the addition of new layers. The high temperature drops of the
bottom right graph are completely normal and take place because
of the addition of new layers, the temperature of which is Tpre,
which are modeled using the same nodes. As a result, the
temperature of the top nodes which represent the current top layer
changes when a new layer is added (adaptive meshing strategy,
described in previous section).
The convergence analysis for the determination of Nl and nh has
been conducted via the same procedure. The difference in the
results between Nl = 5 and Nl = 10 was small enough to assume
convergence for Nl = 5; for the same reasons, nh = 5 has been
selected between nh = 5 and nh = 10.
The use of a 2D space instead of a 3D one brings about the
following issue: in the 3D space, the heat of the laser is offered to
an area of the part (laser spot area), while in the 2D one, it is offered
to a line segment (spot diameter). The procedure followed for the
calculation of a correction factor (as aforementioned in
Section “Modeling approach”) is the following: a simulation of
the thermal history of the same part, using common process
parameters and boundary conditions, has been conducted both in a
3D part, using FEM, and in the 2D model that has been developed in
this study. In this way, a power correction factor that allows the 2D
FD model to produce equivalent results to those of a 3D FEM one,
Table 1
Parameters used in all the convergence analysis tests.
Part dimensions in m
x-axis
z-axis
Process parameters
Spot size
Scanner head speed
Layer thickness
Material: aluminum
Thermal properties
Reflectivity
Model parameters
Material layers described by non-adaptive mesh
Material layers described by adaptive mesh
0.004 m
0.012 m
0.0005 m
0.004 m/s
0.001 m
Figs. 8,9 and 11
82% [35]
2
2
Please cite this article in press as: P. Foteinopoulos, et al., On thermal modeling of Additive Manufacturing processes, NULL (2017), https://doi.
org/10.1016/j.cirpj.2017.09.007
G Model
CIRPJ 441 No. of Pages 18
P. Foteinopoulos et al. / NULL xxx (2017) xxx–xxx
11
Fig. 12. Temperature over time of nodes of different layers for NL = 6.
Fig. 13. Temperature over time of nodes of different layers for NL = 9.
has been calculated. The parameters of the analysis can be seen in
Table 2. The results of the 3D FEM and the 2D FD model can be seen
in Fig. 15. Important to be mentioned is that the thickness of the
part is much smaller compared to its length. The convection
boundary condition has been used for all the surfaces. The material
properties that have been used are the ones presented in this
paper’s main body (Figs. 8–11). In addition, zero reflectivity has
been assumed for those tests only.
In order to prove that this calibration is not geometrydependent, a second simulation has been conducted on a part
of different dimensions. The results derived from the analysis of
the two different parts, using the 3D FEM and the 2D FD model, can
be seen in Fig. 16. This 2D power correction factor that has been
calculated and validated by the two aforementioned comparisons,
minimizes the accuracy loss caused by the 2D space, which now is
on a par with that of an equivalent 3D one, when referring to thinwalled parts. Moreover, the correction factor addresses shapedimension dependent properties and also constitutes a first,
numerical validation of the developed model.
Results and discussion
The thermal history results and the melt pool dimensions of the
2D FD model that has been developed in this paper have also been
validated with the help of the experimentally validated model of
Ref. [9], which has the same inputs and outputs to the model of this
Please cite this article in press as: P. Foteinopoulos, et al., On thermal modeling of Additive Manufacturing processes, NULL (2017), https://doi.
org/10.1016/j.cirpj.2017.09.007
G Model
CIRPJ 441 No. of Pages 18
12
P. Foteinopoulos et al. / NULL xxx (2017) xxx–xxx
Fig. 14. Temperature over time of nodes of different layers for NL = 12.
Fig. 15. Power Calibration factor calculation (left: FEM, right: developed model).
study, allowing for a direct comparison. Most of the models that
use experimental data for validation, utilize microstructure data
for indirect comparison, as melt pool dimensions have a direct
impact on the resulting microstructure [36]. However, such
Table 2
Parameters used in power correction analysis tests.
Part 1 dimensions
x-axis
z-axis
Part 2 dimensions
x-axis
z-axis
Process parameters
Spot size
Power
FD model parameters
Dx, Dz
Dt
FE model parameters
Element size
Dt
0.03 m
0.005 m
0.002 m
0.001 m
0.0005 m
1666.66 W
0.000125 m
0.0002 s
0.00006 m
0.000006 s
comparisons are more extensive in terms of space and a further
analysis is required in order to provide the necessary correlation
between the different KPIs (melt pool and thermal history
microstructure). The specific study focuses on the presentation
of a modeling approach that will render possible the calculation of
the entire thermal history of a part’s manufacturing, by minimizing
the computational time and memory without a noticeable
accuracy loss; a procedure that requires an extensive presentation.
Consequently, a lengthy correlation section of indirect experimental data with the thermal history and melt pool dimensions results
of the model would be out of the scope of this study and as a result,
the utilization of an easily comparable experimentally validated
model has been preferred instead. The input data that have been
used (common in both models) are presented in Tables 3 and 4, in
which the most important similarities and differences of the two
modeling approaches can be seen. In Fig. 17 there is a depiction of
the temperature variation over time of a node, at the center of the
first layer, as it has been calculated by both models, whereas in
Figs. 18 and 19 the simulation results of the melt pool dimensions
can be seen. A comparison between the results of the two models
Please cite this article in press as: P. Foteinopoulos, et al., On thermal modeling of Additive Manufacturing processes, NULL (2017), https://doi.
org/10.1016/j.cirpj.2017.09.007
G Model
CIRPJ 441 No. of Pages 18
P. Foteinopoulos et al. / NULL xxx (2017) xxx–xxx
13
Fig. 16. Power Calibration factor verification (left: FEM, right: developed model).
Table 3
Common analysis parameters used by the model of this study and that of Ref. [9].
Simulation analysis parameters
Absorptance
Laser spot size
Ambient temperature
Laser power
Scan speed
Part dimensions
Length
Width (only in the 3D analysis of Ref. [9])
Thickness
9% [37]
70 mm
20 C
200 W
200 mm/s
1.54 mm
0.7 mm
0.1 mm
can be seen in Table 5. It can be observed that the temperature
history and melt pool dimensions results of the developed model
are fairly accurate and the slight differences can be attributed to
the fact that a 2D space is used in the presented model, whereas in
the study of Ref. [9] a 3D one is used. As a result, thermal
conduction also takes place in the y direction, which is not taken
into account in the 2D space of this model (x, z-axis).
The use of the adaptive meshing strategy greatly saves
memory space and computational time. More specifically,
memory is minimized by the fact that even very big parts can
be modeled with only a few layers of nodes. As far as
computational time is concerned, the time needed for the LU
decomposition (which is used in this study) of a matrix increases
exponentially to the order of the matrix [38]. Using the adaptive
meshing, the matrix size is kept up to a specified constant
maximum that is user specified (maximum number of nodes in
the z direction). Subsequently, the time saved when the adaptive
Fig. 17. Temperature variation over time at the center of the first layer. Model
developed in this study and model of Ref. [9].
mesh is used increases exponentially (the time needed for the
calculation decreases exponentially) and depends on the size of
the part: the bigger the part, the more time is saved. In the first
part of Appendix section a comparison can be found, which proves
the difference between the computational time and memory
needed when the adaptive or the standard meshing strategy is
followed. The analysis is more than 25 times faster when the
adaptive meshing strategy is used. Furthermore, a comparison of
the accuracy of the two meshing strategies is presented and
minimal differences can be observed.
This model can be used for the calibration of the process
variables for the maximization of the KPIs, as well as for the first
module of a thermo-mechanical model that will use the
Table 4
Similarities and differences between the model of this study and that of Ref. [9].
Model feature
Model of this study
Model of Ref. [9]
Heat flux
Moving heat source over time
Temperature dependent material properties
Porosity taken into account
Temperature dependent convection coefficient
Dimensions of the simulation space
Numerical method
Gaussian distribution
Yes
Thermal conductivity, specific heat capacity and density
Yes
No
2D (can be extended)
Finite Differences
Gaussian distribution
Yes
Thermal conductivity and specific heat capacity
Yes
No
3D
Finite Elements
Please cite this article in press as: P. Foteinopoulos, et al., On thermal modeling of Additive Manufacturing processes, NULL (2017), https://doi.
org/10.1016/j.cirpj.2017.09.007
G Model
CIRPJ 441 No. of Pages 18
14
P. Foteinopoulos et al. / NULL xxx (2017) xxx–xxx
Fig. 18. Temperature profile on the cross-section of the molten pool (model of Ref.
[9]).
thermalhistory data of this model for the calculation of the
thermal distortions and residual stresses. In Fig. 20 there is a
depiction of two more uses to be made of this model. The first one
is that many different simulations can be conducted and their
results can be stored into a library, the utilization of which can
lead to enhanced production efficiency, as there will be no need
for an analysis to be run on the spot; instead the library could be
accessed and the required result could be retrieved. The second
part utilizes the fact that the model developed in this study is
capable of simulating and storing the temperature history of the
entire manufacturing time, of a part manufactured in the SLM
calculation of the entire thermal history is viable, from a
computational time and memory perspective, because of the
adaptive meshing strategy. As a result, the evolution of this
phenomenon can be simulated and the cool down rates of a part
(over space and time), which play a major role to its mechanical
properties, can be calculated. This can be utilized in the
calibration of the process parameters, as well as for the
optimization of the process itself, in order for the optimum
cooling down rates to be achieved, ensuring the desired material
properties and minimizing the warping and part oversizing
phenomena (utilization as input for a thermo-mechanical model).
In the case of a very long part, the simulation strategy shown in
Fig. 21 can be adopted. More specifically, when the length of the
part (x-axis) is greatly larger than its thickness (z-axis), the
temperature outside the orange highlighted area is almost
constant. As a result, in order for the computational time to be
further decreased, the temperature of the nodes outside the orange
highlighted zone can be assumed as constant with a low loss of
accuracy.
Finally, the fact that the model uses a 2D space instead of a 3D
one can be overcome for thin walled parts, through the calibration
of the model. This can be achieved by running the analysis of the
same simple part, both on the developed 2D model and on a 3D
one, using as the third dimension the desired wall thickness (it has
to be small enough to be considered thin-walled). The calculation
of a power input correction factor, which leads to the coincidence
of the results of the two models, can be achieved. In this way, the
calibrated model will be capable of providing accurate results for
thin walled parts, with a thickness of the same order of magnitude
with the part of the 3D analysis that has been used for its
calibration.
Fig. 19. Temperature profile on the cross-section of the molten pool (model developed in this study).
process and can easily be adapted to other PBF processes. The
Table 5
Comparison of results.
Result type
Model of this study
Model of Ref. [9]
Difference
Melt pool length (x direction)
Melt pool length (z direction)
Max temperature
88 mm
53 mm
1112 C
103.8 mm
50.2 mm
1254 C
15.2%
5.6%
11.3%
Please cite this article in press as: P. Foteinopoulos, et al., On thermal modeling of Additive Manufacturing processes, NULL (2017), https://doi.
org/10.1016/j.cirpj.2017.09.007
G Model
CIRPJ 441 No. of Pages 18
P. Foteinopoulos et al. / NULL xxx (2017) xxx–xxx
15
Fig. 20. Optimization framework.
Fig. 21. Simulation strategy for long parts.
Conclusions and future work
The idea behind this paper is the creation of a model capable of
accurately simulating and storing in memory the full temperature
history of a 2D component (or a thin walled 3D part) manufactured
in a powder bed AM process. This was made possible by the
utilization of an adaptive meshing strategy, which dramatically
decreases the computational time and the necessary memory for
the simulation. The thermal properties (namely thermal conductivity, specific heat and density) used in this model are a function of
temperature and porosity. The combination of these actions
constitutes the model not only a fast but also an accurate tool,
capable of simulating the temperature history of the total volume
of a printed part, over the whole production time, enabling the user
to optimize the process parameters for enhanced production
efficiency, from a time and energy perspective. It also provides the
evolution of the phenomenon, including the cooling rates, which
play a major role to the mechanical properties of manufactured
parts.
Subsequently, it enables the user to calibrate the process
parameters in such a way so as for the optimum cooling rates to be
ensured. In addition, the melt pool dimensions can be measured at
any time step, by utilizing the color graded temperature diagrams
and taking into account the material’s melting temperature.
Moreover, the combination of (i) the accuracy of results, (ii) the
very high computational speed, (iii) the minimal computational
memory, (iv) the high ease of application, (v) the easy connectivity
with other modules (use of self-developed algorithm), leads to the
conclusion that the developed model can also be used for machine
control.
Finally, as far as a future work is concerned, the study that has
been presented in this paper is the first step of a complete process
model, which will be capable of establishing relationships between
a part’s KPIs and the process variables (Fig. 1): the development of a
complete thermo-mechanical model will take place. Also, this
paper consists the theoretical background that will be used for a
future 3D extension. The next steps, required for the completion of
this work, will be published in the future.
Acknowledgement
This work was supported by the EU Project Borealis [grant
number 636992] of the European Union’s Horizon 2020 research
and innovation program.
Appendix.
a) Adaptive and non-adaptive mesh: accuracy of results and
calculation time comparison
For the validation of the adaptive mesh, a comparison between
the results of a part modeled utilizing the adaptive meshing and
Please cite this article in press as: P. Foteinopoulos, et al., On thermal modeling of Additive Manufacturing processes, NULL (2017), https://doi.
org/10.1016/j.cirpj.2017.09.007
G Model
CIRPJ 441 No. of Pages 18
16
P. Foteinopoulos et al. / NULL xxx (2017) xxx–xxx
Table 6
Parameters used in the non-adaptive/adaptive analysis tests.
Part dimensions in m
x-axis
z-axis
Process parameters
Spot size
Scanner head speed
Layer thickness
Material: aluminum
Thermal properties
Reflectivity
Adaptive model parameters
Layers of nodes used, build and cooling time of t = 19 s
Layers of nodes used, build and cooling time of t = 60 s
Time needed for analysis, build and cooling time of t = 19 s
Time needed for analysis, build and cooling time of t = 60 s
Non-adaptive model parameters
Layers of nodes used, build and cooling time of t = 19 s
Layers of nodes used, build and cooling time of t = 60 s
Time needed for analysis, build and cooling time of t = 19 s
Time needed for analysis, build and cooling time of t = 60 s
0.008 m
0.02 m
0.0005 m
0.004 m/s
0.001 m
Figs. 8, 9 and 11
82% [35]
36
36
2 min 8 s
9 min 41 s
63
180
4 min 52 s
4 h 11 min 11 s
one utilizing the constant one, has been conducted. The
parameters of the test can be seen in Table 6. Both of the analyses
have been conducted in the same system.
When the analysis reaches t = 19 s, seven material layers, with
thickness of 1 mm each, will have been simulated by both models.
Each layer of the material is modeled by nine nodes per layer in the
z direction, which leads to the use of 63 layers of nodes when using
a standard, non-adaptive, mesh. However, when the adaptive mesh
is used, only 36 layers of nodes are required for the modeling of the
same part. At the end of the analysis, after the simulation of 60 s of
build time, the number of nodes used in the adaptive mesh is
exactly the same (36 layers of nodes) as that used when t = 19 s.
When multiplied with the nodes of the x direction, it leads to 2,304
nodes for the whole part (x and z directions). However, the
standard meshing now needs 180 layers of nodes (z direction),
leading to 11,520 nodes for the whole part (x and z directions). As a
result, this difference leads to a much lower computational time
and memory when the adaptive meshing strategy is used. More
specifically, the time and resources required for the analysis to take
place in MatLAB [39], via the two different mesh types, is the
following: when the adaptive meshing strategy has been used,
0.2 GB of RAM and 581 s (9 min and 41 s) were necessary for the
simulation, whereas when the non-adaptive meshing strategy was
used, 4.25 GB of RAM and 15,071 s (4 h, 11 min and 11 s) were
required. This makes this model over 25 times faster when the
adaptive mesh is used instead of the conventional one, for the
given part. In addition, this difference of the computational time
and the memory will increase exponentially along with the
dimensions of the simulated part. As a result, the adaptive meshing
enables the conduction of simulations that would otherwise have
been very costly, time-consuming and subsequently impractical,
due to the resources that would be required in terms of memory
and computational time. In Figs. 22 and 23 the temperature over
time of the middle nodes of different layers can be seen, using the
adaptive and the non-adaptive meshing respectively. The difference in the bottom right plot of Fig. 23 (adaptive meshing) to that
of Fig. 22 (non-adaptive meshing) is completely normal and it is
attributed to the fact that when the non-adaptive meshing is used,
the creation of a new layer is simulated by the addition of a set of
nine nodes with temperature of Tpre. In the case of the adaptive
meshing, the temperature values of the set of the nine top nodes
are moved to the second from the top set of nodes and the value of
Tpre is given to the top set. The difference seen in the plot is
attributed to this change in the temperature of the top set of nodes.
Taking this fact into account, it can be observed that the plots of
temperature over the time of the two different types of mesh are
almost identical, thus, validating the accuracy of the adaptive
meshing. The same can be concluded by Figs. 24 and 25, in which
the color-map plots of the part under construction at the t = 19 s for
both meshing strategies are depicted.
b) Regression accuracy validation
A validation, concerning the accuracy of the regression used in
the adaptive meshing algorithm, for the calculation of the nodal
temperatures in the new positions of the nodes (change in their z
coordinate), has been performed in this section. This regression
Fig. 22. Temperature over time of nodes of different layers using the non-adaptive mesh.
Please cite this article in press as: P. Foteinopoulos, et al., On thermal modeling of Additive Manufacturing processes, NULL (2017), https://doi.
org/10.1016/j.cirpj.2017.09.007
G Model
CIRPJ 441 No. of Pages 18
P. Foteinopoulos et al. / NULL xxx (2017) xxx–xxx
17
Fig. 23. Temperature over time of nodes of different layers using the adaptive mesh.
calculates the temperature value each node should have when
moved to a new position, based on the existing solution of nodal
temperatures of the current time-step. A procedure like the one
used in the changing of the node position, in the z direction, will
take place and then its results will be validated, by conducting the
reverse procedure (Fig. 26). More specifically: starting from the
initial, green-star set of points and with the use of Eq. (15), which is
also used for the adaptive meshing, the x-values of a new set of
points are calculated. The z-values of this new set of points is
calculated by utilizing the spline regression, which is under
validation, having as input the x and z of the green-star set of points
and the new x values, which were previously calculated via
Eq. (15). The new set of points created is depicted in red-circles. At
this point, the reverse procedure will be followed: using as input
for the spline regression a x set of values, equal to those of the
green-star set of points and the x and z values of the red-circle set
of points, the black-X set of points is created. It can be observed in
Fig. 26 that those points coincide almost perfectly with the original
green-star set of points, thus validating the accuracy of this
regression. Taking into account the way that the temperature
function changes over time, the spline regression is suitable for the
task in hand and will maintain its accuracy.
Fig. 24. Colour-map plot of the part under construction at t = 19 s using the nonadaptive mesh.
Fig. 25. Colour-map plot of the part under construction at t = 19 s using the adaptive
mesh.
Please cite this article in press as: P. Foteinopoulos, et al., On thermal modeling of Additive Manufacturing processes, NULL (2017), https://doi.
org/10.1016/j.cirpj.2017.09.007
G Model
CIRPJ 441 No. of Pages 18
18
P. Foteinopoulos et al. / NULL xxx (2017) xxx–xxx
[15]
[16]
[17]
[18]
[19]
[20]
Fig. 26. Validation of the accuracy of the spline regression that is used for the
calculation of the temperature value of the nodes, when their position changes in
the z direction in the context of adaptive meshing. (For interpretation of the
references to color in the text, the reader is referred to the web version of this
article.)
[21]
[22]
[23]
References
[1] ISO/ASTM52921-13. 2013, Standard Terminology for Additive Manufacturingcoordinate Systems and Test Methodologies.ASTM International, West Conshohocken, PA. http://dx.doi.org/10.1520/ISOASTM52921-13.
[2] Hopkinson NN., Hague RR.J.M.JM, Dickens PP.M.M, (Eds.) (2006), Rapid
Manufacturing: An Industrial Revolution for the Digital age. John Wiley &
Sons.
[3] Levy, G.N., Schindel, R., Kruth, J.P., 2003, Rapid manufacturing and rapid tooling
with layer manufacturing (LM) technologies, state of the art and future
perspectives. CIRP Ann Manuf Technol, 52:589–609. http://dx.doi.org/
10.1016/S0007-8506(07)60206-6.
[4] Adam, G.A.O., Zimmer, D., 2014, Design for additive manufacturing—element
transitions and aggregated structures. CIRP J Manuf Sci Technol, 7:20–28.
http://dx.doi.org/10.1016/j.cirpj.2013.10.001.
[5] Childs, T.H.C., Berzins, M., Ryder, G.R., Tontowi, A., 1999, Selective laser sintering of an amorphous polymer—simulations and experiments. Proc Inst Mech
Eng B, 213:333–349. http://dx.doi.org/10.1243/0954405991516822.
[6] Bikas, H., Stavropoulos, P., Chryssolouris, G., 2016, Additive manufacturing
methods and modelling approaches: a critical review. Int J Adv Manuf Technol,
83/1–4: 389–405. http://dx.doi.org/10.1007/s00170-015-7576-2.
[7] Gibson, I., Rosen, D., Stucker, B., 2014, Additive Manufacturing Technologies:
3D Printing, Rapid Prototyping, and Direct Digital Manufacturing. 2nd ed.
Springer. http://dx.doi.org/10.1007/978-1-4939-2113-3.
[8] Michaleris, P., 2014, Modeling metal deposition in heat transfer analyses of
additive manufacturing processes. Finite Elem Anal Des, 86:51–60. http://dx.
doi.org/10.1016/j.finel.2014.04.003.
[9] Li, Y., Gu, D., 2014, Parametric analysis of thermal behavior during selective
laser melting additive manufacturing of aluminum alloy powder. Mater Des,
63:856–867. http://dx.doi.org/10.1016/j.matdes.2014.07.006.
[10] Riedlbauer, D., Drexler, M., Drummer, D., Steinmann, P., Mergheim, J., 2014,
Modelling, simulation and experimental validation of heat transfer in selective
laser melting of the polymeric material PA12. Comput Mater Sci, 93:239–248.
http://dx.doi.org/10.1016/j.commatsci.2014.06.046.
[11] Dong, L., Makradi, A., Ahzi, S., Remond, Y., 2009, Three-dimensional transient
finite element analysis of the selective laser sintering process. J Mater Process
Technol, 209:700–706. http://dx.doi.org/10.1016/j.jmatprotec.2008.02.040.
[12] Kolossov, S., Boillat, E., Glardon, R., Fischer, P., Locher, M., 2004, 3D FE
simulation for temperature evolution in the selective laser sintering process.
Int J Mach Tools Manuf, 44:117–123. http://dx.doi.org/10.1016/j.ijmachtools.2003.10.019.
[13] Liu, F.R., Zhang, Q., Zhou, W.P., Zhao, J.J., Chen, J.M., 2012, Micro scale 3D FEM
simulation on thermal evolution within the porous structure in selective laser
sintering. J Mater Process Technol, 212:2058–2065. http://dx.doi.org/10.1016/j.
jmatprotec.2012.05.010.
[14] Loh, L.-E., Chua, C.-K., Yeong, W.-Y., Song, J., Mapar, M., Sing, S.-L., Liu, Z.-H.,
Zhang, D.-Q., 2015, Numerical investigation and an effective modelling on the
[24]
[25]
[26]
[27]
[28]
[29]
[30]
[31]
[32]
[33]
[34]
[35]
[36]
[37]
[38]
[39]
selective laser melting (SLM) process with aluminium alloy 6061. Int J Heat
Mass
Transf, 80:288–300.
http://dx.doi.org/10.1016/j.ijheatmasstransfer.2014.09.014.
Ganeriwala, R., Zohdi, T.I., 2014, Multiphysics modeling and simulation of
selective laser sintering manufacturing processes. Procedia CIRP, 14:299–304.
http://dx.doi.org/10.1016/j.procir.2014.03.015.
Paul, R., Anand, S., 2012, Process energy analysis and optimization in selective
laser sintering. J Manuf Syst, 31:429–437. http://dx.doi.org/10.1016/j.
jmsy.2012.07.004.
Kovaleva, I., Kovalev, O., Smurov, I., 2014, Model of heat and mass transfer in
random packing layer of powder particles in selective laser melting. Phys
Procedia, 56:400–410. http://dx.doi.org/10.1016/j.phpro.2014.08.143.
Mohanty, S., Hattel, J.H., 2014, Numerical model based reliability estimation of
selective laser melting process. Phys Procedia, 56:379–389. http://dx.doi.org/
10.1016/j.phpro.2014.08.135.
Dai, D., Gu, D., 2014, Thermal behavior and densification mechanism during
selective laser melting of copper matrix composites: simulation and experiments.
Mater
Des,
55:482–491.
http://dx.doi.org/10.1016/j.
matdes.2013.10.006.
Gürtler, F.-J., Karg, M., Leitz, K.-H., Schmidt, M., 2013, Simulation of laser beam
melting of steel powders using the three-dimensional volume of fluid method.
Phys Procedia, 41:881–886. http://dx.doi.org/10.1016/j.phpro.2013.03.162.
Pastras, G., Fysikopoulos, A., Giannoulis, C., Chryssolouris, G., 2015, A numerical approach to modeling keyhole laser welding. Int J Adv Manuf Technol,
78:723–736. http://dx.doi.org/10.1007/s00170-014-6674-x.
Hu, H., Argyropoulos, S.A., 1996, Mathematical modelling of solidification and
melting: a review. Model Simul Mater Sci Eng, 4:371. http://dx.doi.org/
10.1088/0965-0393/4/4/004.
Monzón, M.D., Ortega, Z., Martínez, A., Ortega, F., 2015, Standardization in
additive manufacturing: activities carried out by international organizations
and projects. Int J Adv Manuf Technol, 76:1111–1121. http://dx.doi.org/
10.1007/s00170-014-6334-1.
Kruth, J.-P., Levy, G., Klocke, F., Childs, T.H.C., 2007, Consolidation phenomena
in laser and powder-bed based layered manufacturing. CIRP Ann Manuf
Technol, 56:730–759. http://dx.doi.org/10.1016/j.cirp.2007.10.004.
Chryssolouris, G., 2006, Manufacturing Systems: Theory and Practice. 2nd ed.
Springer, New York.
Schleifenbaum, H., Meiners, W., Wissenbach, K., Hinke, C., 2010, Individualized
production by means of high power selective laser melting. CIRP J Manuf Sci
Technol, 2:161–169. http://dx.doi.org/10.1016/j.cirpj.2010.03.005.
Kruth, J.P., Froyen, L., Van Vaerenbergh, J., Mercelis, P., Rombouts, M., Lauwers,
B., 2004, Selective laser melting of iron-based powder. J Mater Process Technol,
149:616–622. http://dx.doi.org/10.1016/j.jmatprotec.2003.11.051.
Sheikhi, M., Malek Ghaini, F., Assadi, H., 2015, Prediction of solidification
cracking in pulsed laser welding of 2024 aluminum alloy. Acta Mater,
82:491–502. http://dx.doi.org/10.1016/j.actamat.2014.09.002.
Cernuschi, F., Ahmaniemi, S., Vuoristo, P., Mäntylä, T., 2004, Modelling of
thermal conductivity of porous materials: application to thick thermal barrier
coatings. J Eur Ceram Soc, 24:2657–2667. http://dx.doi.org/10.1016/j.jeurceramsoc.2003.09.012.
German, R.M., 2009, Handbook of Mathematical Relations in Particulate
Materials Processing. John Wiley & Sons.
Lemmon, E.W., Jacobsen, R.T., Penoncello, S.G., Friend, D.G., 2000, Thermodynamic properties of air and mixtures of nitrogen argon, and oxygen from 60 to
2000 K at pressures to 2000 Mpa. J Phys Chem Ref Data, 29:331–385.
Chase, M.W., Davies, C.A., Downey, J.R., Frurip, D.J., McDonald, R.A., Syverud, A.
N., 1985, Janaf thermochemical tables-2. J Phys Chem Ref Data, 14:927–1856.
Ho, C.Y., Powell, R.W., Liley, P.E., 1974, Thermal Conductivity of the Elements: A
Comprehensive Review.the American Chemical Society and the American
Institute of Physics for the National Bureau of Standards, West Lafayette,
Indiana.
Jensen, J.E., Stewart, R.G., Tuttle, W.A., Brechna, H., 1980, Brookhaven National
Laboratory Selected Cryogenic Data Notebook: Sections I–IX. Brookhaven
National Laboratory.
Samsonov, G.V., 1968, Handbook of the Physicochemical Properties of the
Elements.Plenum Publishing Corporation, New York.
Gade, R., Moeslund, T.B., 2014, Thermal cameras and applications: a survey.
Mach Vis Appl, 25:245–262. http://dx.doi.org/10.1007/s00138-013-0570-5.
Louvis, E., Fox, P., Sutcliffe, C.J., 2011, Selective laser melting of aluminium
components. J Mater Process Technol, 211:275–284. http://dx.doi.org/10.1016/
j.jmatprotec.2010.09.019.
Bunch, J.R., Hopcroft, J.E., 1974, Triangular factorization and inversion by fast
matrix multiplication. Math Comput, 28:231–236. http://dx.doi.org/10.1090/
S0025-5718-1974-0331751-8.
The MathWorks, Inc. 2016, MATLAB Release R2016a.The MathWorks, Inc.,
Natick, Massachusetts, United States.
Please cite this article in press as: P. Foteinopoulos, et al., On thermal modeling of Additive Manufacturing processes, NULL (2017), https://doi.
org/10.1016/j.cirpj.2017.09.007
Документ
Категория
Без категории
Просмотров
4
Размер файла
5 985 Кб
Теги
2017, 007, cirpj
1/--страниц
Пожаловаться на содержимое документа