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 deﬁned 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 speciﬁcally 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 difﬁcult 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 ﬁnal 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 ﬁrst 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 speciﬁcally, 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 ﬁeld 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 ﬁnite 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 coefﬁcient 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 coefﬁcient 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 densiﬁcation 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 ﬁelds 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 ﬁelds (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, speciﬁc 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 veriﬁed with the use of data from Ref. [9]. In this study, the modeling of the ﬁrst 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, ﬁne 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 signiﬁcance 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 deﬁned 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 ﬁrst 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 ﬂowchart 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 ﬂowchart 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, speciﬁc 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 inﬁnite) 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 proﬁle 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 speciﬁed 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 coefﬁcient. 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 coefﬁcient, 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 speciﬁed 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 deﬁned 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 deﬁned 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 ﬂux 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: coefﬁcient, 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 deﬁned 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 speciﬁcally, 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 deﬁned. After the ns timesteps have passed, the heating boundary condition is applied to the ﬁrst 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 speciﬁcally, 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 ﬁrst 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 speciﬁed 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 speciﬁcally, 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 coefﬁcient 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 ﬁeld 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 ﬂowchart 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 solidiﬁed. 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, speciﬁc 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 speciﬁc heat of the porous medium: cpr ¼ ð1 f Þcm þ f ca ð20Þ However, the term ca cannot be ignored here, as the speciﬁc heat of air is in the same order of magnitude as that of aluminum. In Fig. 7, the function of the speciﬁc heat capacity of air versus temperature, which has been used for the calculation of the speciﬁc 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 deﬁned 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/solidiﬁcation phase change is modeled using the apparent heat capacity method [22], according to which, speciﬁc heat is artiﬁcially increased in a small temperature range around the melting point, so as to encompass the latent heat of the melting/solidiﬁcation 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 artiﬁcially 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 speciﬁcally, the value of the apparent heat capacity is deﬁned 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. Speciﬁc 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. Speciﬁc 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 speciﬁc heat capacity over temperature, in which both the standard heat capacity of the material and the apparent heat capacity are included. The artiﬁcially 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 ﬁrst, in z-axis, node in each layer) of the part will be used for the test at hand. With reference to the ﬁrst 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. Speciﬁc 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 Reﬂectivity 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 reﬂectivity 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 ﬁrst, 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 speciﬁc 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 ﬁrst 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 veriﬁcation (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 speciﬁcally, 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 speciﬁed constant maximum that is user speciﬁed (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 ﬁrst 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 ﬁrst 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 ﬁrst 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 ﬂux Moving heat source over time Temperature dependent material properties Porosity taken into account Temperature dependent convection coefﬁcient Dimensions of the simulation space Numerical method Gaussian distribution Yes Thermal conductivity, speciﬁc heat capacity and density Yes No 2D (can be extended) Finite Differences Gaussian distribution Yes Thermal conductivity and speciﬁc 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 proﬁle 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 ﬁrst 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 efﬁciency, 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 speciﬁcally, 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 proﬁle 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, speciﬁc 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 efﬁciency, 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 ﬁrst 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 Reﬂectivity 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 speciﬁcally, 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 speciﬁcally: 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.ﬁnel.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 ﬁnite 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 densiﬁcation 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 ﬂuid 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 solidiﬁcation 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 solidiﬁcation 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

1/--страниц