Computational modelling of the selective laser sintering process

Universidade do Minho Escola de Engenharia João da Costa Teixeira e Castro Computational Modelling of the Selective Laser Sintering Process October 2022 UMinho | 2022 João Castro Computational Modelling of the Selective Laser Sintering Process João da Costa Teixeira e Castro Computational Modelling of the Selective Laser Sintering Process Master's Dissertation Integrated Masters in Polymer Engineering Work accomplished under the supervision of: Professor Doctor João Miguel Nóbrega Doctor Ricardo Costa Universidade do Minho Escola de Engenharia October 2022 Direitos de autor e condições de utilização do trabalho por terceiros Este é um trabalho académico que pode ser utilizado por terceiros desde que respeitadas as regras e boas práticas internacionalmente aceites, no que concerne aos direitos de autor e direitos conexos. Assim, o presente trabalho pode ser utilizado nos termos previstos na licença abaixo indicada. Caso o utilizador necessite de permissão para poder fazer um uso do trabalho em condições não previstas no licenciamento indicado, deverá contactar o autor, através do RepositóriUM da Universidade do Minho. Licença concedida aos utilizadores deste trabalho: Atribuição CC BY https://creativecommons.org/licenses/by/4.0/ ii Acknowledgements The present work is a culmination of a five years journey from which I leave a better and wiser person, so, before advancing, I could not fail to thank those who contributed to mine and this work development. I would like to start by thanking my professor and supervisor João Miguel Nóbrega for introducing and teaching me computational mechanics, igniting a flame of curiosity and interest that illuminated this work. I would also like to extend my thanks to my second supervisor, Ricardo Costa, who, despite only meeting at the end of my academic path, greatly contributed to the development of my Master’s Dissertation. It is a certainty that both their insights and knowledge on the work’s subject steered me through it and, without any doubt, allowed its completion. Financially, I acknowledge the support of National Funds through FCT - Portuguese Foundation for Science and Technology, Reference UID/CTM/50025/2019 and UIDB/04436/2020, and project SIFA - Sistema Inteligente de Fabricação Aditiva (POCI 01-0247-FEDER-047108). I also acknowledge the computing facilities support by Search-ON2: Revitalization of HPC Infrastructure of UMinho (project no. NORTE-07-0162-FEDER-000086), co-funded by the North Portugal Regional Operational Programme (ON.2 -- O Novo Norte), under the National Strategic Reference Framework (NSRF), through the European Regional Development Fund (ERDF). Finally, a special thanks to my family and friends who have accompanied me for so long and, without knowing, have also contributed to this work, not with their knowledge, but with their fondness, which has, undoubtedly, supported me through this long journey. I would like to end by dedicating my dissertation to my father, who I am sure would have loved to see me finish it. To everyone, I am so grateful for your support, as I could have not done it without you. iii Statement of integrity I hereby declare having conducted this academic work with integrity. I confirm that I have not used plagiarism or any form of undue use of information or falsification of results along the process leading to its elaboration. I further declare that I have fully acknowledged the Code of Ethical Conduct of the University of Minho. University of Minho, 31st October 2022 Name: João da Costa Teixeira e Castro Signature: iv Resumo A Manufatura Aditiva (AM do inglês Additive Manufacturing) tem ganho popularidade em várias indústrias importantes e exigentes devido à sua capacidade de produzir peças com geometrias complexas e pouco desperdício. Como uma das suas mais populares técnicas, a Sinterização Seletiva a Laser (SLS do inglês Selective Laser Sintering) é muito procurada por diversas indústrias que pretendem substituir processos convencionais mais dispendiosos. No entanto, o processo de SLS é intrinsecamente complexo devido aos vários fenómenos multi-fisícos e são necessários mais estudos para obter uma melhor perceção dos mesmos. Isto tem originado um elevado interesse académico em otimizar o processo para que este cumpra os requisitos industriais. Grande parte destas otimizações são feitas através de métodos experimentais que são demorados, dispendiosos e nem sempre resultam nas configurações ótimas. Este enquadramento tem motivado investigadores a recorrer à modelação computacional com o objetivo de entender melhor o processo, de modo a antecipar e corrigir defeitos. O objetivo principal do presente trabalho foi desenvolver um modelo capaz de simular o processo de SLS para aplicações poliméricas, em código de distribuição livre, à escala do tamanho da partícula. Como são necessárias abordagens distintas para simular com precisão cada etapa do processo, diferentes métodos numéricos foram aplicados para desenvolver uma ferramenta capaz de estudar o impacto, numa secção representativa da cama de pó, dos parâmetros físicos que podem ser ajustados no processo. O trabalho desenvolvido integrou diversas etapas, começando por um estudo extenso dos aspetos teóricos do processo de SLS que visou a familiarização com os fenómenos envolvidos, o desenrolar do processo, os seus parâmetros e respetivas influências, assim como a avaliação das limitações e desafios existentes. Este estudo foi seguido por uma análise detalhada dos modelos mais populares empregues para representar os principais fenómenos associados ao processo e do nível de precisão das abordagens, com base nas simplificações consideradas. Um conjunto de ferramentas computacionais foi posteriormente apresentado e os seus respetivos modelos selecionados, quando possível, de acordo com a revisão bibliográfica efetuada. Por último, vários testes foram executados, visando uma validação experimental qualitativa dos códigos utilzados, para garantir que o modelo em uso era adequado para simular o processo, permitindo o estudo e observação da influência dos principais parâmetros do processo e a evolução da sinterização. Os desenvolvimentos obtidos representam um avanço significativo para a simulação do processo de SLS. Com o uso de software open-source (LIGGGHTS e OpenFOAM), vários estudos foram feitos numa geometria realista e, apesar da ausência de dados experimentais suficientes e mais detalhados, os resultados da simulação mostraram-se bem correlacionados com os usados para comparação. Em suma, o trabalho efetuado permitiu concluir que a ferramenta desenvolvida apresenta um elevado potencial para estudar, com detalhe, o processo de SLS e a influência dos seus parâmetros e, deste modo, contribuir para a sua otimização. Palavras-chave: Sinterização Seletiva a Laser (SLS), polímeros, modelação numérica, OpenFOAM, LIGGGHTS v Abstract Additive Manufacturing (AM) has increased in popularity in numerous important and demanding industries due to the capability of manufacturing parts with complex geometries with little wastage. As one of its most popular techniques, Selective Laser Sintering (SLS) is sought after by several industries that aim to replace conventional and more expensive processes. However, the SLS process is intrinsically complex due to the various underlying multi-physics phenomena and more studies are needed to obtain more insights about it. These has resulted in many academical interests to optimize the process and allow it to achieve industrial standards. Most of these optimization attempts are performed through experimental methods that are time consuming, expensive and do not always provide the optimal configurations. This has lead researchers to resort to computational modelling, aiming at better understanding the process to anticipate and fix the defects. The main objective of the present work was to develop a model capable of simulating the SLS process for polymeric applications, within an open-source framework, at particle length scale. Since distinct approaches are required for accurately simulating each step of the SLS process, different numerical methods were employed to develop a tool capable of studying the impact, in a representative section of the powder bed, of the physical parameters that can be adjusted in the process. The developed work comprised several steps, starting with an extensive study of the theoretical aspects of the SLS process that aimed at the acquaintance with the involved phenomena, process unwind, its parameters and their influence, as well as evaluating the existing limitations and challenges. This study was then followed by a detailed analysis of the most common employed models to represent the major phenomena and of the accuracy level of the approaches, based on the employed simplifications. A set of computational tools was then exhibited and their built in models were selected, when possible, according to the precedent literature review. Lastly, various tests were carried to obtain an experimental qualitative validation of the used code, to assure that the used model was adequate to simulate the process, allowing the study and observation of the principal parameter influence and sintering progression. The achieved developments represent a significant advance towards the SLS process simulation. With the use of open-source software (LIGGGHTS e OpenFOAM), several studies were performed on a realistic geometry and, despite the absence of enough and more detailed experimental data, the simulation results are in agreement with the ones used for comparison. Overall, the accomplished work allowed to conclude that the developed tool constitutes a great potential to study, in detail, the SLS process and its parameters influence and, therefore, contribute to its optimization. Keywords: Selective Laser Sintering (SLS), polymer, numerical modelling, OpenFOAM, LIGGGHTS vi Table of contents Chapter 1 -- Introduction ......................................................... 1 1.1 AdditiveManufacturing ....................................................... 1 1.2 PowderBedFusion .......................................................... 2 1.3 Numerical Modelling and Simulation State-of-the-Art . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.4 MotivationandObjectives ..................................................... 7 1.5 DissertationOutline .......................................................... 8 Chapter 2 -- Selective Laser Sintering ............................................ 10 2.1 ProcessDescription ......................................................... 10 2.2 BindingMechanism ......................................................... 12 2.2.1 SolidStateSintering .................................................... 13 2.2.2 ChemicallyInducedBinding .............................................. 13 2.2.3 Liquid Phase Sintering/Partial Melting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.2.4 FullMelting ........................................................... 15 2.3 PolymersinSLS ............................................................ 15 2.4 ProcessVariablesandtheirInfluence ........................................... 17 2.4.1 ProcessParameters .................................................... 18 2.4.1.1 LaserandScanParameters ......................................... 18 2.4.1.2 BuildParameters .................................................. 21 2.4.2 PowderProperties ...................................................... 22 2.4.2.1 PowderManufacture ............................................... 23 2.4.2.2 Particle Size and Shape Influence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.4.3 Material Thermal and Physical Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.4.3.1 Coalescence ...................................................... 26 2.4.3.2 SinteringWindow .................................................. 29 2.4.3.3 Crystallization ..................................................... 31 2.4.3.4 Reusability ....................................................... 32 Chapter 3 -- Modelling and Simulation ............................................ 35 3.1 ThermalModel ............................................................. 36 3.1.1 LatentHeat ........................................................... 36 3.1.2 LaserBeam ........................................................... 38 3.1.3 HeatLosses ........................................................... 41 vii 1.3. Numerical Modelling and Simulation State-of-the-Art Method (OTM) and Smoothed Particle Hydrodynamics (SPH), while mesh-based counterparts comprise the Finite Element Method (FEM), the Finite Difference Method (FDM) and the Finite Volume Method (FVM). The DEM is the most common mesh-free method employed for simulating the process, but mainly for the powder bed distribution [ 24 ]. It accounts for the discrete nature of the powder bed and considers the interactions between each particle at their contact point, being well suited for the study of granular media. Despite its elementary Newtonian mechanics, it has been extended to include other interactions between particles, such as heat transfer [ 25 ], that allowed the temperature distribution study within a realistic powder bed [ 26 ]. Unfortunately, this approach fails to consider the surrounding air or the fusion and coalescence of the particles. On the other end, mesh-based methods, such as the FEM, firstly introduced in 1956 [ 27 ], stands out for its extensive use for simulating the SLS process [ 28 -- 30 ] due to its simplicity, efficiency and mainly the vast amount of literature on the subject over the past decades. Even though it provides useful information, such as the temperature evolution, liquid fraction and residual stresses, this method makes it difficult to account for volume shrinkage, layer deposition and the particle nature of the bed [ 31 ], which motivates some authors to opt for the FVM [ 32 ]. Also, since meshed-based methods turn out to be more cumbersome to simulate the powder distribution and deposition, DEM is often coupled with those approaches to take care of this step [ 33 ]. Unlike conventional manufacturing processes, such as injection and extrusion, whose process have been extensively studied and optimized in the past, the optimization of the SLS process is more recent. It is usually performed through experimental methods [ 34 ], but such approaches are time-consuming and expensive. The primary purpose of most simulations on the SLS process is to assess the influence of the various process parameters into the final part mechanical and geometrical quality. Although the process being the same, distinct modelling approaches are taken for different materials, mainly polymers and metals, due to their different behaviours. For example, metals are exposed to a much wider temperature range, so considering their temperature dependent properties is imperative [ 35 ]. Also, due to their reduced viscosity, liquid metal flows are heavily influenced by density variations and many other effects that need to be accurately modelled [ 36 ]. The numerical modelling of temperature distribution, deformation and thermal stresses of a SLS process is similar to that of multi-pass welding [ 37 ], which is an older technology that has been developed and modelled for more than four decades. Unfortunately, for the SLS process, the model has to be more advanced due to all the additional complex physical phenomena involved, such as the laser interactions with the powder, namely absorption, reflection and scattering. Despite being a relatively new technology, a significant amount of attempts to simulate the process was already done. Some approaches aimed at developing 5 1.3. Numerical Modelling and Simulation State-of-the-Art a model and validating it by comparing simulation results with experimental ones [ 38 ], while others aimed at optimizing a process or study the influence of certain parameters. Considering the powder bed as a continuous media, Childs et al. [ 39 ] investigated the influence of laser power and scanning speed and Dai and Shaw [ 40 ] the influence of scanning pattern and speed on the temperature and stress fields. Bugeda et al. [ 41 ] also developed a model to simulate the 3D sintering process of a single track, obtaining data such as the temperature field, solid fraction and sintering depth. Similarly, Dong et al. [ 42 ] modelled the transient heat transfer during the sintering process while considering the phase transition, but also the temperature dependent material properties to evaluate the temperature and density distribution. A more accurate model, considering a more realistic volumetric heat source, was developed by Riedlbauer et al. [ 43 ] to predict melt pool dimensions and validate the results with experimental data. Still related to melt pool, Peyre et al. [ 44 ] and Foroozmehr et al. [ 45 ] tested the effects of the process parameters on its temperature distribution and dimensions, validating the results with experimental tests. Bierwisch et al. [ 46 ] developed a 2D model with a realistic powder bed for two different materials, where the various phases of the process, namely scan, heat dissipation and cooling, were considered. Osmanlic et al. [ 47 ] developed a similar model with an accurate laser source, where phenomena such as reflection, refraction and attenuation were considered, but for much shorter time scales and with the objective of assessing the bed density influence on the process. Some authors compiled and reviewed several studies [ 36 , 48 -- 51 ] and the interested reader is addressed to those publications. Besides isolated models, there are also some commercial softwares that offer the option to simulate the LS process, such as ANSYS [ 52 ], SIMULIA [ 53 ], COMSOL [ 54 ] and NETFABB [ 55 ]. Also, for DEM simulations, LIGGGHTS [56] and LAMMPS [57] are commonly used. Modelling SLS processes comes with several challenges. There are non-linearities to be considered, such as the temperature dependence of the material properties. Besides that, the problem must be discretized both spatially and temporally and, depending on the model scale, very small time scales and high mesh resolutions might be needed, which requires a significant amount of computational power. For the interaction of the laser with the material, there are several different physical behaviours to consider, such as absorption, reflection and radiation, along with the subsequent heat transfer, phase transformation, moving interface between the different phases, fluid flow and possible chemical reactions. Moreover, the powder distribution also affects most of the previously mentioned aspects, so an adequate geometry definition is also needed. In the end, for an accurate AM simulation, one must model the addition of the material, have a realistic model for the heat source, account for the thermal losses, consider the temperature dependence of some properties, model 6 1.4. Motivation and Objectives elasto-plastic stresses and strains, couple the thermo-mechanical behaviour in a manageable way and, finally, model phenomena that change material properties [ 58 ]. However, these requirements are computationally demanding, therefore it is important to verify if simpler models are also able to provide useful insights about the process. Because of all the complexity associated, one aims to start with a simple approach, while maintaining the main features taking place. 1.4 Motivation and Objectives Since their arrival, AM technologies have shown an immense potential to fabricate parts in a way that, if the quality is achieved, would make conventional and other recent alternative methods look outdated. They proved to have the ability to improve everyone’s life and be worth their attention, something that can be seen by the multiple industrial applications and the consequent exponential market growth [ 7 ]. With PBF being one of the most popular AM methods, SLS falls no short in potential when compared to the other AM processes. Indeed, the ability to easily manufacture parts with very complex geometries that would be impossible, or extremely hard, to achieve with other methods and would take multiple other technologies and waste a lot more energy and material, is the most prominent advantage of the SLS process and the main reason for its increasing attention in the industry. For these reasons, increasing efforts have been put into further improving the SLS process. Despite the ever-increasing interest of the industry in the SLS process, this recent technology is still not fully dominated and the current limitations are still substantial, attending to the demanding and specific thermo-mechanical and visual characteristics of the parts required for several applications. In that regard, the numerical simulation turns out to be a very useful approach for a better understanding of the process, since one can observe and analyse details that would otherwise not be accessible with the conventional experimental approaches. Moreover, with a good understanding of the process and its parameters, the numerical simulation allows to acquire useful insight on the process control, to prevent defects, and to guide process optimization, to improve specific desired properties. Nonetheless, such a task is challenging due to the complexity of predicting accurately the effects of the process parameters changes, motivated by the process’ complex nature. The common practice to perform numerical studies of AM processes is to acquire a software license, which is expensive. Even after that, it is nearly (or even totally) impossible to get access to the code and adapt it to specific needs, because they are proprietary. This framework results in the 7 1.5. Dissertation Outline financial impossibility to utilize the software by everyone but large companies, and the users cannot know if the modelling approach is rigorous enough or adequate for their specific needs. For these reasons, the use of open-source software, such as OpenFOAM [ 59 ], is gaining attention and becoming more widespread with growing communities of users and developers, where one can, free of any cost, use the libraries and contribute to the code development. Unfortunately, the literature on the use of open-source software for the computational modelling of AM processes is still very scarce. Most of the research works cited in the previous section do not make use of open-source software and, when they do, little to no detail is given to the reader, making it very hard to draw adequate conclusions. Moreover, a large number of researchers chooses to simulate the process at a macro scale, with a continuous medium powder bed, that provides no information on phenomena at the particle length scale, such as the porosity or sintering evolution. This work is mostly motivated by the absence of enough and detailed information on the topic of SLS simulation, especially within OpenFOAM, as well as the need to develop a model of easy access to conduct the envisaged studies. Based on the aforementioned information, the main objective for this dissertation is the development of a computational model capable of accurately simulating the SLS process, within an open-source framework. The developed tool should also allow studying the impact of the physical parameters that can be adjusted in the process. To achieve that, on the geometry side, one must represent the powder bed as realistic as possible, with adequate powder particle spatial and size distribution. For the computational model, a solver available in the OpenFOAM library, which was previously studied to assess its adequacy [ 60 ], will be further investigated and appropriately adapted to better represent the physical phenomena involved in the SLS process. Towards the end of the project, with the developed computational model, a representative section of the powder bed will be simulated, allowing the studies of the sintering progression and the influence of different underlying parameters. 1.5 Dissertation Outline The remaining contents of the dissertation are organised into several parts that allow the fulfilment of the outlined objectives for this work. After presenting the motivation and objectives, the SLS process is described in detail in Chapter 2, covering the main mechanisms, phenomena and process variables as well as the associated challenges and limitations. Then, Chapter 3 addresses the process modelling and simulation, based on the state-of-the-art mathematical models and numerical methods employed to describe and solve the associated phenomena, to build a base to critically select the most 8 1.5. Dissertation Outline adequate computational tools. Afterwards, in Chapter 4, the selected computational tools are described together with the careful characterization of the underlying models being solved. In Chapter 5, the computational model setup is illustrated, comprising the simulation of the representative geometry, mesh generation and material properties, followed by the testing and assessment of the underlying models to ensure that they can correctly represent the process. At the end of this chapter, the proposed case studies are introduced, whose results are presented and discussed in Chapter 6, to evaluate the solver’s capabilities and effectiveness to simulate the SLS process. The document ends with Chapter 7, where the main conclusions are drawn and future work is proposed to surpass the identified limitations. 9 CHAPTER 2Selective Laser Sintering 2.1 Process Description As mentioned in Section 1.2, SLS is a process derived from PBF, which implies that it inherits a set of basic principles, such as the three main steps that summarize the process: powder recoating, energy input and material coalescence. A general representation of the SLS process is illustrated in Figure 2.1. Feed chamber Feed chamberBuild chamber Laser Lenses Mirror Levelling roller Laser beam Powder feed supply Powder feed piston Powder bed Build piston Sintered parts Figure 2.1: General representation of the SLS process. SLS machines are divided into two main zones, the build chamber, where, as the name suggests, the parts are built, and the feed chamber, that stores the necessary powder for the process and provides it to the build chamber when needed. Two feed chambers are commonly employed so the roller can create a layer starting from any side, removing the need for a return phase. Also, both can collect the excess powder dragged by it, whereas if only one feed chamber is present, an additional zone is required at the opposite side to catch the excess powder. Although possible to use the laser to 10 2.1. Process Description melt the material from room temperature, in SLS, the powder is pre-heated to a temperature just a few degrees below its melting point. Besides minimizing the laser power requirements and allowing a faster scanning speed, this reduces temperature gradients between scanned and non-scanned particles, facilitating the fusion between layers and keeping a more uniform temperature on the bed, which helps to prevent defects, such as local shrinkage and thermal distortions caused by non-uniform thermal expansions and contractions [ 61 ]. The pre-heating is done in bulk, being achieved through a combination of heat lamps and conductive elements above and around the perimeter of the build and feed chambers, respectively, and is maintained throughout the entire build process. Considering the possibility of the material thermal degradation due to these high temperatures during both the pre-heating and building phases, the process takes place inside a closed chamber filled with an inert gas, generally nitrogen, in order to minimize oxidation. After the pre-heating phase is complete, the other steps of the process start taking place. Initially, the build piston starts at its highest position and a levelling roller spreads material from the feed supply to the build platform, creating the first powder layer. The amount of material dragged by the roller is always superior to the required for the layer, to assure that it is completely filled. Some machines use a blade instead of a roller to move the powder, however, it has been found that their dragging mechanism produces a less dense powder bed with higher layer surface roughness, when compared to the one obtained with the counter rotating roller [ 62 ]. When the layer is complete, a laser beam is directed to the powder bed through a set of lenses and is directed to the desired positions, at a pre-selected speed, with the use of mirrors. The laser type should be selected according to the material present in the bed. For example, for metals, the laser interacts with their electrons, while in polymers the interactions is with the resonant vibration of larger segments of the material molecule. It is known that for aliphatic polymers, the CH 2 segments resonates near 10.6 µ m [ 63 ], which justifies the wide use of CO 2 laser beams (capable of outputting infrared wavelength from 9 to 11 µ m) at that wavelength [64]. As the laser spot moves, it scans the material, providing it enough energy to melt and coalesce with adjacent granules, including the ones on the previous layer (if it exists). The unscanned powder stays in place and works as support for the next layer. The amount of energy provided to the material is crucial to determine the particle consolidation and, consequently, the final part porosity. After the scanning of the present layer is complete, the build piston lowers the powder bed by the thickness of one layer and a new powder layer is spread on top. These operations are repeated until the desired object is fully built and, after that, the heaters are turned off and the part is left to cool inside the machine chamber, so the cooling is as uniform as possible and within an inert ambient, in order to 11 2.2. Binding Mechanism minimize thermal defects. Although mentioning multiple times that one part is built, multiple objects can be manufactured inside the same bed, as long as the material is the same and there is enough space. A bed with multiple parts is often referred to as part cake. The more parts present in one bed, the more efficient the production process is, so this comes as an advantage. Once cooled, the parts can be removed, cleaned off and, if necessary, subjected to finishing operations. 2.2 Binding Mechanism SLS is a complex manufacturing technique with various process variables that can be adjusted according to specific needs. However, some adjustments, mainly related to the material and energy density, result in significant differences in the process approach and development, consequently affecting the final part properties. These differences come from the distinct binding mechanisms between the powdered material and, because of that, in order to facilitate communication, different terminologies were given to different technologies according to the predominant binding mechanism. These include four main groups, solid state sintering, chemically induced binding, liquid phase sintering and full melting [ 65 ]. It is important to note that, despite being called sintering, most of the approaches do actually melt some of the involved material, as it will be discussed bellow. Independently of the mechanism, because there are no additional forces applied to the powder bed, the coalescence between liquid particles is mainly driven by surface tension forces that pull them together. This phenomenon depends on the material and process parameters and a schematic representation is illustrated in Figure 2.2, where the formation and development of a "neck" between particles is visible. Increasing temperature and/or time Unsintered particles Neck formation Increased necking Decreased porosity Pore Figure 2.2: Formation of the "neck" between particles. 12 2.2. Binding Mechanism 2.2.1 Solid State Sintering Solid State Sintering (SSS) is a thermal consolidation process that occurs bellow the material’s melting temperature. It aims at achieving densification by reshaping of the powder in such a way that no liquid is formed. This reshaping is a consequence of reducing energy by eliminating a solid-gas interface and replacing it with a solid-solid one through the diffusion of atoms in the solid state. The tendency of the particles to move to their lowest energy state will "fuse" the adjacent ones together through the formation of a "neck", similar to that shown in Figure 2.2. The "neck" development should be accentuated to minimize porosity, however, the process is very slow because it depends only on diffusion on the solid state so, in order to achieve high sintering rates, it requires either higher temperatures or longer processing times, which comes as a disadvantage, since efficiency both in manufacturing time and energy translate to reduced costs. Because of the negative effects, mainly the time consumption, SSS is only applied to very specific materials, such as ceramics, or for post processing from other techniques [66]. 2.2.2 Chemically Induced Binding This technique is based on thermal activated chemical reactions between two types of powders, or between the powder and a gas, to bind together the particles. It is a fusion mechanism used mainly for ceramics, for example, processing SiC within an atmosphere that contains O 2 , forms SiO 2 that binds itself to the original material, forming a product of SiO 2 and SiC [ 66 ]. This technique originates parts with high porosity, making it rely on post-process techniques, such as infiltration or high temperature furnace sintering, to reduce it and obtain appropriate properties [67]. 2.2.3 Liquid Phase Sintering/Partial Melting Liquid Phase Sintering (LPS) is a very versatile mechanism that unifies different types of technologies. Different from other techniques, there are two types of materials with distinct functions, where one, the binder, becomes fully molten and acts as a glue to bind the other, the structural, that remains solid throughout the production process. This approach is useful for cases where the structural material, that contains the desired properties, is hard to process, so another material is used to bind it. The binder and structural materials can be combined in the powder bed in three different ways: 13 2.2. Binding Mechanism  Separate Particles - This approach is the simplest of the three. Creating a good mixture of both structural and binder particles is, in many cases, sufficient to obtain good results. There are, however, some aspects that need to be considered to assure a good quality mixture that will yield parts with good properties. For example, the binder particles must be smaller than the structural ones so the packing is more efficient, leading to less shrinkage, higher density and less porosity. It also increases the binding capability, because the smaller particles will easily fit between the larger ones and, due to their reduced size, melt faster. However, because the laser source moves quickly, there is not always sufficient heat or surface tension forces for the binder to fully flow, meaning that parts formed from separate particles present a higher porosity, thus requiring post processing in a furnace to reduce it [66].  Composite particles - These particles comprise both the structural and binder, which are obtained by mechanically alloying the two materials though grinding cast, extruded or moulded mixtures. This approach aims at reducing the porosity in parts resulting from separate particles and also offers a better surface finish [ 68 ]. However, composite particles are also used when the combination of properties from both materials is desired. The most common composite particle is glass-filled nylon, an atypical application where the binder does not act as a glue but instead is reinforced by the structural material, the glass [69].  Coated Particles - The third and last possibility to combine the two materials is to coat the structural particle with the binder material. In the other two cases, as the laser scans the powder bed, the radiation will be absorbed and distributed through all the particles, independent of the material, meaning that the energy is not being focused on the desired constituent, the binder. With coated particles this is not an issue because the binder is always at the particle surface, which is the zone that receives more radiation and, therefore, will absorb more energy and be forced to melt faster. Besides that, the fact that every structural particle is coated with the binder results in an overall better binding efficiency, not only because the surface is the best place for particles to fuse/sinter together, but also because there are no places where the binder is missing due to the particles random mixture, as is seen in the other cases [66]. There are, however, some special cases where there is no distinction between the particles, meaning that there is only one type of material present in the powder bed. For these cases, instead of distinguishing binder and structural particles, there is a distinction of the material physical state. In this situation, "partial melting" is a better suited nomenclature. During the process, when the heat supplied 14 2.4. Process Variables and their Influence complete review and discussion can be found in Jia et al., 2021 [83]. 2.4.1.2 Build Parameters Process build parameters comprise the layer thickness and build temperatures to be maintained throughout the process. Generally, they are less complicated to setup due to their more easily predictable influence. Starting with the layer thickness, it is a parameter that dictates how much the build piston should move downwards after the scan of each layer. It is defined before the machine virtually slices the digital model and, depending on its value, more or less layers are generated. Higher layer thickness values generate less layers for a determined part, which translates to faster production speeds, but sacrifices dimensional accuracy and surface roughness since the resolution is inferior. The maximum value for this parameter is determined by the depth of penetration of the laser energy that, itself, depends on the material absorbency. Each material absorbs a set of electromagnetic wavelengths differently, hence the importance of selecting an appropriate material-laser pair, as stated in Section 2.4.1.1. Usually, the absorbency is measured using thin films with a known thickness, however, in practical cases, the powder bed layer cannot be simplified to such configuration and the penetration depth is always greater due to the empty spaces between particles. Therefore, parameters such as particle size and bed density also influence the penetration depth. Researchers experimentally investigate this depth by scanning and gradually increasing the thickness of each layer and measuring the transmitted power. For PA12, it was reported a consistent absorption of over 90% for a layer thickness of 200 µ m [ 85 ]. Paired with the penetration depth, the fusion depth, that indicates the limit of laser penetration that also results in particle melting, is usually also evaluated [ 86 ]. The knowledge of these values is of great importance to adequately select a maximum layer thickness, which should be smaller than the fusion depth, to assure full coalescence within the layer and a strong connection with the previous one. On the other end, the minimum value for the layer thickness is restricted by the particle size or, in other words, each layer must contain at least one particle of thickness. Generally, the layer thickness is two to three times the average particle size (see Section 2.4.2.1). The build temperatures are selected according to the material and energy related parameters. All systems are capable of controlling the temperature on the build chamber for the powder bed and its surrounding atmosphere by using a combination of heating elements that surround the powder bed. Some machines have the capability of heating the feed chamber as well, which comes with advantages 21 2.4. Process Variables and their Influence and disadvantages. On one side, heating the powder before spreading it reduces the temperature difference between the bed and the fresh layer, which helps reducing temperature gradients and premature cooling. On the other side, the material is held at high temperatures for longer periods, which influences its properties and reusability (more information in Section 2.4.3.4). As explained in Section 2.1, controlling the temperature, as well as pre-heating the powder, is essential for a good process unwinding. Despite being counter-intuitive, most of the energy provided comes from the heating elements and the laser is simply used to overcome the material enthalpy of fusion and thus melting it. Each material and build will have an unique ideal build temperature, but, in general, it is beneficial to select the highest possible temperature that does not promote the consolidation of unscanned particles, to minimize laser power requirements, thermal gradients between scanned and unscanned particles and thermal expansion promoted by the laser [ 87 ]. For semi-crystalline polymers, this temperature is a few degrees below 𝑇𝑚 and equal, or very close, to 𝑇𝑔 , for the amorphous ones. If the build temperature is too high (assuming that it never reaches 𝑇𝑚 ), the energy required to melt the particles will be so small that non-scanned particles adjacent to scanned ones will still receive enough energy to melt and coalesce. This is often referred to as part growth and is a defect that reduces accuracy and surface finish. Moreover, keeping the powder at even higher temperatures throughout the process further hinders its reusability and also reduces flowability, as mentioned in Section 2.4.3.4 and Section 2.4.2, respectively. Conversely, if the selected build temperature is too low, premature cooling is more prone to happen, accompanied by subsequent contraction and its possible defects (further discussed in Section 2.4.3.3). It is important to keep in mind that, due to the larger size of the builds and machines, there is always a variation on the temperature despite being set to constant. This variation is attributed to some causes, such as losses of energy from the material to the surrounding elements (machine walls and the naturally generated convection currents in the atmosphere), energy gains from the laser scanned zones and the heating elements inability to provide an even heat source [ 61 ]. Although minimal, this variation can have a significant negative impact on the process, especially for polymers with a narrow processing window (see Section 2.4.3.2). 2.4.2 Powder Properties As suggested by the name of the parent process (PBF), the raw material is used in its powdered form. Although extrinsic, the powder properties interact significantly with intrinsic properties and dictate 22 2.4. Process Variables and their Influence the success of the manufacturing process. Selecting and successfully obtaining a powder with the optimal properties is essential to produce consistent, predictable and high quality parts. 2.4.2.1 Powder Manufacture The majority of polymers are not produced directly in powder form and must be converted to an appropriate particle size prior to sintering. That is usually done via milling/grinding procedures and solvent precipitation, although existing other less common methods that use immiscible blends and melt prilling. Melt-based and mechanical size reduction approaches can be employed for all the materials, while solution-based ones are exclusive for polymers. Each method results in particles with differences in shape and size distribution, however that is still not well defined or predictable and, unfortunately, it is an area that lacks detailed research studies. Nonetheless, the particles, in terms of shape, fall into three main categories, spherical, potato-shaped and irregular, as shown in Figure 2.6 [78]. 100 μm100 μm100 μm Spherical particles Potato-shaped particles Irregular particles Figure 2.6: Particles obtained by different production technologies [78]. Grinding Mechanical approaches, such as grinding, are challenging to use with polymers because of their viscoelastic behaviour. Usually, they require cryogenic temperatures to exhibit a more brittle response that allows fracturing, as it is typical in grinding, which increases output and produces a finer product [ 88 ]. Because the grinding process converts a large amount of energy into heat, which works against keeping the material cool, other grinding approaches are being evaluated [ 88 ]. As of now, despite the efforts, the powders produced by this technique are not recommended due to their non-spherical nature and high portion of small particles that, together, work against their free flow and difficult powder spreading [89]. 23 2.4. Process Variables and their Influence Solvent Precipitation The solvent precipitation method is inspired by the dispersion polymerization method and relies on the capability of dissolving together a polymer and a dispersant, creating a solution. The reactive monomer is soluble in the dispersant medium, but the reactant polymer is not, so, when a solvent is added to the solution and the polymer is formed, it immediately precipitates and self assembles into micro-spheres. Another variation is spray drying, where a solution is atomized into a hot drying gas, usually air, evaporating the liquid phase and leaving small particles of the solid phase. In both approaches, the resulting particles have a desirable shape and size distribution [ 90 , 91 ]. Unfortunately, the solvent precipitation method is not compatible with all polymers and only a fraction can be obtained through this technique. Immiscible Blends This approach aims at controlling the morphology of the powder during its manufacturing phase to obtain spherical particles. It works by mixing two distinct polymers in a controlled way to induce a strain rate that will spontaneous form the desired micro-spheres and then, a solvent dissolves the unwanted polymer and the spheres are obtained [ 92 ]. Although interesting, this approach is highly dependent on the material properties and imposed parameters, making it challenging to use. 2.4.2.2 Particle Size and Shape Influence The SLS process has high feedstock requirements, which further restricts the use of many materials. Indeed, one of the reasons that contributes to the short list of available materials in SLS is the need for very specific and challenging to obtain powder properties, such as particles shape, size and size distribution. Particle Size and Shape The selected particle size is, generally, as small as possible based on the benefits from denser powder beds and faster sintering. Larger particles can be used, but are less desirable. It is possible to obtain very fine powders, however, regardless of their shape, from sizes lower than 45 µ m, static Van der Waals forces start to create agglomerates that oppose their flow and prevents the recoating of a 24 2.4. Process Variables and their Influence smooth layer [ 79 ], heavily influencing their usability. Therefore, particle size is usually between 45 µ m and 90 µ m [ 93 ]. It is not possible to obtain uniform sized particles and a distribution of sized will always exist so, technically speaking, the powder will have a particle size distribution (PSD). A good PSD should be centred around the desired particle size and not too wide to avoid uneven melting, with smaller particles prematurely fused and larger particles still solid. To assure a good PSD and to diminish problems induced by smaller particles, it is common to filter particles outside the selected range [ 93 ]. Despite the mentioned problems, smaller particles can be used when proper treatments to disrupt Van der Waals forces are employed [ 94 ]. For particle shape, as expected, spherical (or as spherical as possible) particles are easier to spread and will naturally generate better packing. If the shape is irregular it is very likely to interact and form arches that lead to empty spaces and, consequently, higher porosity [95]. The different shapes are shown in Figure 2.6. Powder Flow and Bed Density Particle morphology has a direct impact in the powder deposition phase, also called recoating, that, itself, influences the ease of processing and final product properties [ 61 ]. The energy parameters may be more intuitively linked to creating a better product, but having a densely packed bed prior to sintering will minimize void space and, overall, increase desired properties, being a crucial factor for good quality printed parts [ 79 ]. Powder flow is a secondary property that depends on the particles size and shape and is generally used to characterize the powder in order to predict its performance. With an increasing flow, density will be higher, so this parameter should be as high as possible. The Hausner ratio is widely considered as a good predictor of flowability [ 96 ]. It numerically compares the powder in its highest packed state, often called tapped density, to the resulting packing from a standardized flow condition. Both packing types are illustrated in Figure 2.7. For PA12, this ratio ranges from 0.4 to 0.6, which means that the apparent density is 40% to 60% of the tapped density. This method is often considered insufficiently precise and more PBF related approaches are employed, using modern powder rheometers or other powder spreading machines that aim at mimicking the real process [ 97 ]. It is important to note that polymers and metals are submitted to different tests due to their different cohesion levels. Also, despite most of the tests being performed at ambient conditions [ 96 ], it should be imperative to consider the higher temperatures in the chamber, that are well above the polymer 𝑇𝑔 and promote changes in its properties, which reduce flowability. Differences in humidity should also be measured, as it can promote different levels of agglomeration for hygroscopic polymers [ 98 ]. Correctly 25 2.4. Process Variables and their Influence Tapped density Apparent density Figure 2.7: Illustration and comparison of tapped and apparent density. analysing the powder flow and consequent bed density is extremely important because the variations in these parameters are considered to be leading causes for uncertainty and lack of repeatability in the process [ 99 ]. Numerous tests have been preformed under various conditions for a variety of materials and, overall, it has been observed that flowability increases with the PSD decrease, is improved for larger particles and will decrease in moist environments [98]. 2.4.3 Material Thermal and Physical Properties The chemical structure of the polymer particles is, as expected, very important because it determines other parameters that dictate the usage of the material in the process. Altering it to obtain properties that favour the SLS process is far from trivial, forcing researchers to select materials that naturally fulfil all process needs, which greatly reduces the diversity of polymers that can be employed in the process, as mentioned in Section 2.3. 2.4.3.1 Coalescence The particles coalescence is one, if not the most, important characteristics of the process. Indeed, numerous process parameters are carefully selected to assure maximum coalescence between scanned particles. Even the employed materials are selected based on their capability to coalesce during the process. If the particles do not fuse together at a sufficiently high rate, the final consolidated part will have poor mechanical properties because of the insufficient adhesion between particles and consequent high porosity. Coalescence takes place after the material is scanned and melts, thus changing to a state with high enough mobility. The amount of mobility that a material has in liquid state determines the rate of coalescence, as shown by the Frenkel model in Equation (2.2) [100]: 𝑥 𝑟2 =3𝜎𝑡 2𝑟𝜂0 ,(2.2) 26 2.4. Process Variables and their Influence where 𝑥 is the growing neck between two particles, 𝑟 the original particle radius, 𝜎 the surface tension in the liquid state, 𝑡 is the time and 𝜂0 the zero-shear viscosity. This model fails to consider the varying particle radius that results from coalescence, making it only valid for the early stages of coalescence, while the particles spherical shape is preserved. This led other researchers to make adjustments to the model, further increasing its similarity to the SLS particle behaviour [ 101 ]. Nonetheless, despite this inaccuracy, the Frenkel model is adequate to assess parameters influence and proves that high coalescence rates are achieved when surface tension is high and viscosity is low. Surface Tension At small length scales such as the size of the particles present in the powder bed, surface tension has a prominent impact on the molten material movement, especially because the inertia forces are much smaller. Surface tension can be described as the molecules attractive force, present at the surface of a liquid, towards each other. It is also known as surface energy, although this nomenclature is more correct when referring to the attractive forces between solid molecules, instead of liquid. Any system will always tend to move to its lowest energy state, so, when given enough mobility, that is, when the particles melt, because surface tension pulls surface molecules to each other, the liquid will move towards an equilibrium state that can be represented according to the Young equation [102]: 𝜎sg =𝜎sl +𝜎lg cos𝜃, (2.3) where 𝜎sg , 𝜎sl and 𝜎lg are the interfacial surface tensions for solid-gas, solid-liquid and liquid-gas, respectively. An illustration of an arbitrary equilibrium state is shown in Figure 2.8. Based on the equation and the figure, one can predict the behaviour of the liquid if the interfacial surface tensions for the multiple phases are known. For example, assuming that 𝜎lg stays constant, if 𝜎sg increases, then 𝜃 will decrease. Contrarily, for a constant 𝜎sg , 𝜃 will increase with 𝜎lg . In summary, when 𝜎sg is higher than 𝜎lg , the adhesive forces are stronger than the cohesive forces and the liquid will wet the solid ( 𝜃 < 90◦ ). On the opposite side, for 𝜎sg < 𝜎lg the liquid will not wet the solid, but instead tend to form a spherical shape ( 𝜃 > 90◦ ). For SLS, the predominant surface tension is 𝜎lg because the molten particles will be mostly in contact with others in the same state, therefore, the interaction is only between liquid and gas. Since there are no external forces applied to the melt, the surface tension has to be high enough to overcome its high viscosity and promote coalescence. Also, for the external parts of the layer that are in contact with solid particles, it is not desirable that the liquid material totally wets 27 2.4. Process Variables and their Influence Solid Liquid Gas σsl σsg σlg θ Figure 2.8: Three-phase surface tension equilibrium and resulting contact angle. the solid particles, so 𝜎sg should not be much greater than 𝜎lg [ 102 ]. Researchers have analysed the 𝜎lg in fuction of temperature and have reported values, for PBF grade PA12 within its melt range, to be between 0.035-0.040 N/m [102]. Viscosity When in their liquid state, polymers, such as PA12, are very viscous. Since it is desirable to maximize the coalescence rate, according to the Frenkel model, viscosity should be as low as possible. One should note that, in the model, the used viscosity, 𝜂0 , stands for zero-shear viscosity that differs from other values of viscosity for the same material. That is because polymers, being non-Newtonian fluids, have a viscosity that depends not only on temperature [ 103 ], but also shear rate [ 104 ]. Due to these dependencies, tests performed to evaluate viscosity for a certain process must be done at process equivalent conditions. Polymers exhibit a shear thinning behaviour, typical of pseudoplastics, in other words, their viscosity decreases with increasing shear rate. Considering this, it is not adequate to select a material for SLS, a process where there are no forces applied to the melt, based on results from a capillary rheometer that mimics injection moulding/extrusion conditions. For SLS, cone and plate or parallel discs rheometers are more adequate due to the possibility of much lower, close to zero, shear rates [ 105 ]. A polymer viscosity is highly dependent on its molecular weight ( 𝑀𝑤 ) [ 106 ] because, contrarily to most materials, polymers are composed of very long molecular chains, that are a combination of many repeating units, and their length is what determines the molecular weight of the material. These chains have a spatial distribution and are able to interact with each other in a similar way that a bunch of wires would. Higher 𝑀𝑤values mean that there are longer chains and, consequently, their interactions, mainly entanglements, are harder to overcome, which, in practical terms, translates into difficulties from the molten material to flow or, in other words, an increase in viscosity. Also, for the same 𝑀𝑤 , depending on the entanglement intensity, the viscosity will vary [ 107 ]. This justifies the shear thinning behaviour, a consequence of increasing the shear rate, that forces the molecular chains to untangle, hence reducing viscosity. While 28 2.4. Process Variables and their Influence low viscosity is desirable for fast sintering, the necessary low 𝑀𝑤 to achieve it is not, because it impairs the mechanical properties of the final part. The process applies no external forces to the material so the viscosity is totally dependent on 𝑀𝑤 . To circumvent such inconvenience, the polymer is fabricated with an initial low molecular weight that can increase later, during the process. For example, for PA12, this is achieved during its synthesis, more specifically, the ring opening polyaddition of lauryl lactam, where the amino and carboxylic acid end groups are left purposely "unprotected", instead of being chemically terminated [ 108 ]. Later, when higher temperatures are achieved inside the build chamber, these groups combine through condensation, increasing 𝑀𝑤 [ 109 ]. This comes with an advantage of increasing mechanical properties through the later increase in 𝑀𝑤 , while maintaining an initial low 𝜂0 for high consolidation. Unfortunately, all the unscanned material in the chamber will also increase in 𝑀𝑤, severely hindering its reusability [110], as further discussed in Section 2.4.3.4. 2.4.3.2 Sintering Window Particles coalescence is important to assure good final part properties and, as mentioned before, it is not a fast process due to the typical high viscosity of polymers. Obtaining a material capable of sintering at high rates within process conditions is ideal, but it can easily be overthrown if one very important condition, keeping the material in liquid state, is not met. Independently of the selected material, coalescence is only possible in its liquid state, so it is imperative that the scanned particles stay liquid long enough to effectively coalesce. In polymers, for which the solidification temperature is known as crystallization temperature, it is common to have distinct melting and crystallization temperatures. This temperature gap between the onset of melting ( 𝑇OM ) and crystallization ( 𝑇OC ) (illustrated in Figure 2.9) is known as sintering/process/supercooling window. When in between these temperatures (assuming that the polymer is in liquid state and cooling down), solidification/crystallization is not possible and, despite the high viscosity, the material has enough mobility to allow coalescence. Avoiding early crystallization also helps prevent cooling associated defects, as further discussed in Section 2.4.3.3. One material limitation that prevents the use of many polymers in SLS are too narrow, or non-existent, sintering windows. It is challenging to fabricate a part when extremely precise laser parameters and build temperatures are required to avoid premature cooling. Materials that have larger sintering windows are the most desirable since they allow for wider processing temperature ranges, meaning that they are less sensitive to temperature fluctuations and do not require high process precision, allowing the process parameters to be adjusted to optimize other conditions. The lack of a sintering window is one of the reasons why amorphous polymers are less desirable [75]. 29 2.4. Process Variables and their Influence The compatibility of a new material with the SLS process is often evaluated through a well-known technique for polymer characterization, Differential Scanning Calorimetry (DSC). This technique measures the energy that has to be transferred to, or from, a sample to induce temperature change, thus allowing the identification of thermal transitions, along with their respective enthalpy/entropy. Figure 2.9 presents a schematic representation of a result from a DSC test for a PA12. It shows the various information obtained such as the 𝑇OM and 𝑇OC but, most importantly, it highlights the interval between them, the sintering window. This large interval, coupled with narrow melting range, for a more centred melting point, low melting temperature, so less energy is required, and a high melting enthalpy, to minimize unwanted sintering from conduction, is why, as mentioned in Section 2.3, PA12 is one of the most common used material for SLS. DSC tests, although also melting and solidifying the sample, are not a totally accurate representation of the SLS process because the heating and cooling rates, contrary to what is observed in the process, are fixed. The values for 𝑇OM and 𝑇OC are not the same for different heat rates and the material will crystallize sooner when the cooling rates are very low [ 74 ]. This results in the narrowing of the sintering window and is the reason why, despite the material still being at a temperature where it should be, theoretically, liquid, after being held at that temperature for long periods (very slow cooling rates), crystallization starts to happen [74]. TEC TOC TC TM SINTERING WINDOW TOM TEM Heating Cooling Temperature (K) Heat flow (W) Endo Exo PA12 TOM - Onset of melting TEM - End of melting TM - Peak of melting ΔH - Enthalpy of fusion TOC - Onset of crystallization TEC - End of crystallization TC - Peak of crystallization ΔH Figure 2.9: Schematic representation of a DSC test for a PA12 highlighting the gap between the melting and crystallization onset. 30 3.1. Thermal Model modified specific heat term as [125]: 𝐶𝑝=𝐶𝑝0+Δ𝐻𝑚 √︁𝜋(𝑇mf −𝑇ms)exp−(𝑇−𝑇ms)2 (𝑇mf −𝑇ms)2,(3.3) where 𝐶𝑝0 is the base specific heat of the material to which the enthalpy heat of fusion Δ𝐻𝑚 is added, distributed between the fusion start and finish temperatures, 𝑇ms and 𝑇mf, through a Gaussian Law. Another option would be to consider the latent heat of melting and the heat released during crystallization as a heat sink and heat source, respectively. With this approach, instead of changing the 𝐶𝑝 term, the contributions are added to the right-hand side of the energy conservation equation (Equation (3.1)). The latent heat of fusion (heat sink), can be modelled as [126]: ¤ 𝑄sink =𝛿ℎ𝛿 𝑓liq 𝛿𝑡 ,(3.4) where 𝑓liq is the volumetric fraction of liquid phase, calculated as [126]: 𝑓liq =(𝑇−𝑇sol) 𝑇liq −𝑇sol,(3.5) where, given the range of melting temperature Δ𝑇=𝑇liq −𝑇sol , its value varies from 0, when 𝑇≤𝑇sol to 1, when 𝑇≥𝑇liq , always satisfying 𝑓liq +𝑓sol =1 . 𝛿ℎ is the total absorbed latent heat over a melting temperature range for a two-phase medium given by 𝛿ℎ =𝐻liq −𝐻sol , with 𝐻 representing the enthalpy of both liquid and solid constituents [127]. Both terms are calculated as [127]: 𝐻liq =𝜌liq ∫𝑇 𝑇sol 𝐶𝑝liq 𝑑𝜃 +𝛿𝐻, (3.6) 𝐻sol =𝜌sol ∫𝑇 𝑇liq 𝐶𝑝sol 𝑑𝜃, (3.7) where 𝜌liq , 𝜌sol , 𝐶𝑝liq and 𝐶𝑝sol are the density and specific heat of the liquid and solid material, respectively, and 𝛿𝐻 is the reference enthalpy for the liquid. Oppositely, in the cooling phase, polymers crystallize and release energy, acting as a heat source that can be modelled by the following equation [32]: ¤ 𝑄source =𝜌sol 𝛿𝛼𝑐 𝛿𝑡 Δ𝐻𝑐,(3.8) 37 3.1. Thermal Model with 𝛼𝑐 representing the crystallinity rate and Δ𝐻𝑐 the crystallization enthalpy. The term 𝛿𝛼𝑐 𝛿𝑡 accounts for the evolution of the relative degree of crystallization as a function of time and temperature and can be described by the Nakamura’s model [128]: 𝛿𝛼𝑐 𝛿𝑡 =𝑛 𝐾 (𝑇)(1−𝛼𝑐)ln 1 1−𝛼𝑐𝑛−1 𝑛 ,(3.9) where 𝑛 is the Avrami index that depends on the growth geometry of crystallites and 𝐾(𝑇) is the non-isothermal crystallization rate. 𝐾(𝑇) can be related to 𝑛 and determined by the following expression [129]: 𝐾(𝑇)=ln2 1 𝑛1 𝑡1/2,(3.10) with 𝑡1/2 representing the half crystallization time for a defined isothermal temperature. The determination of that time can be expressed by the Hoffman-Lauritzen theory [130]: 1 𝑡1/2=𝐾0exp−𝑈 𝑅(𝑇−𝑇∞)exp−𝐾𝐺(𝑇+𝑇0) 2𝑇2(𝑇0−𝑇),(3.11) where 𝐾0 is a temperature independent constant, 𝑈 is the activation energy of the crystallization transport with the universal value of 6270 J/mol, 𝑅 the universal gas constant that is equal to 8.314 J/mol/K, 𝑇∞ is the temperature at which the crystallization transport finishes, that takes the value of 𝑇𝑔−30𝐾 , 𝐾𝐺 is a parameter related to the nucleation characteristics and 𝑇0 is the equilibrium melting point. All the unknown parameters are determined by means of DSC. 3.1.2 Laser Beam When modelling the laser beam it is important to formulate an equation that takes into account its shape and power distribution. From Section 2.4.1.1, it is known that the types of laser beams employed into SLS processes have a Gaussian intensity profile and a circular shape, hence only the radius being mentioned as a parameter. When treated as a surface heat source, following such conditions, the power intensity of the laser 𝑄LB can be expressed as [28]: 𝑄LB(𝑥, 𝑦)=𝑃𝜂 2𝜋𝜎2𝑒 𝑥2+𝑦2 −2𝜎2,(3.12) 38 3.1. Thermal Model where 𝑃 is the laser beam power and 𝜂 the material absorptivity, that can be calculated, if the material reflectivity ( 𝜆 ) is measured experimentally, with 𝜂=1−𝜆 [ 131 ]. The position on the surface is given by coordinates 𝑥 and 𝑦 and the point of highest intensity, the centre of the beam, corresponds to the origin of this local coordinate system. The standard deviation 𝜎 indicates the variation in intensity along the laser spot radius 𝑟LB . One should aim at obtaining an adequate distribution throughout the radius, where the maximum intensity is at the centre and decays to values close to zero as it approaches the periphery, so the selection of 𝜎 cannot be random. To assure a good distribution that fits inside the spot radius, researchers define the standard deviation proportionally to it, as to satisfy 𝜎=𝑟LB/2 . This slightly modifies Equation (3.12) into [51, 131]: 𝑄LB(𝑥, 𝑦)=2𝑃𝜂 𝜋𝑟2 LB 𝑒 −2𝑥2+𝑦2 𝑟2 LB .(3.13) This way, when the distance to the centre is equal to 𝑟LB , the intensity value is reduced by a factor of 𝑒2 of the maximum value. This approach also makes the equation easier to read, because it replaces the standard deviation of the intensity distribution with the laser radius, an easier parameter to define. However, from previous attained information, it is known that the laser has the ability to go through the material and, in fact, it is necessary that it does, to effectively melt the entire thickness of the scanned powder. Modelling the laser as a surface 2D heat source is simpler, but far from realistic, since, in reality, it exhibits a volumetric 3D behaviour. In this case, the energy attenuates as it penetrates the material because it is being absorbed until, eventually, it extinguishes. This behaviour follows the Beer-Lambert equation, which is an exponential equation that can be written as [44, 132]: 𝑄(𝑧)=𝑄0𝑒−𝛼𝑧,(3.14) where 𝑄 is the transmitted laser power in the thickness direction 𝑧 , 𝑄0 the surface heat source from Equation (3.13) and 𝛼 the extinction coefficient, with the term 𝑒−𝛼𝑧 representing the exponential decay along the powder bed thickness. From this equation, the energy distribution function can be written as [132]: 𝑞(𝑧)=𝑄0−𝑄0𝑒−𝛼𝑧,(3.15) 39 3.1. Thermal Model that, after a simple derivation, is transformed into [132]: 𝑞(𝑧)=𝛼𝑄0𝑒−𝛼𝑧.(3.16) If 𝑄0 is replaced with the surface heat source from Equation (3.13) , the equation of the volumetric heat source is obtained and it can be added to Equation (3.1), in the ¤ 𝑄term, as a heat source: 𝑄𝑣(𝑥, 𝑦, 𝑧)=𝛼2𝑃𝜂 𝜋𝑟2 LB 𝑒 −2𝑥2+𝑦2 𝑟2 LB 𝑒−𝛼𝑧.(3.17) Besides the energy distribution of the beam, given by a Gaussian distribution along the radius, and the energy absorption of the material modelled by the Beer-Lambert equation, more phenomena are present during the scanning phase. One popular modelling approach for the laser beam is to disctretize its domain into a finite number of rays so that each individual one can interact with intersecting media. Such consideration is important because the laser is an electromagnetic wave and it is well known that they suffer an alteration in their direction when they encounter a change in medium, as illustrated in Figure 3.1. θ n1 αα n2 L LR LT n Figure 3.1: Laser reflection and refraction when changing media. Looking at this illustration, the initial ray, represented by the vector L, when changing medium, from n 1 to n 2 , gets terminated and two new rays, L R and L T , are cast in different directions, representing the reflection and refraction suffered by L, respectively. The reflected ray, since its travelling in the same medium, has the same angle 𝛼 as its former ray relatively to the surface normal n. Because the initial ray gets split in two, each one will keep a portion of its total energy 𝐸𝐼 . By using the Fresnel equation for 40 3.1. Thermal Model unpolarized light and non-magnetic media, the portion of reflected energy 𝐸𝑅 can be calculated as [ 47 ]: 𝐸𝑅=1 2"𝑛1cos𝛼−𝑛2cos𝜃 𝑛1cos𝛼+𝑛2cos𝜃2 +𝑛1cos𝜃−𝑛2cos𝛼 𝑛1cos𝜃+𝑛2cos𝛼2#𝐸𝐼.(3.18) The remaining non reflected energy, given by 𝐸𝑇=𝐸𝐼−𝐸𝑅 is refracted through the new medium. Contrarily to the reflected ray, the refracted one travels in a different medium at which the wave will propagate at a different velocity, causing an alteration on the angle in relation to n. The newly formed refraction angle 𝜃can be determined using Snell’s law [47]: sin𝜃=𝑛1 𝑛2 sin𝛼, (3.19) where 𝑛1and 𝑛2are the refractive index of each medium. 3.1.3 Heat Losses During the building process, the energy losses from the powder bed to the environment are minimal when compared to the energy absorbed by the material. Nonetheless, they exist in the form of convection to the surrounding air, emissive radiation and, even less pronounced, conduction to the external walls and material ablation. The process takes place inside an insulated chamber where there are not many elements taking heat away from the system, which leads to some models neglecting their contribution to the energy equation [ 133 ]. However, despite being common practice to ignore energy losses when modelling small sections, because the errors scale with time and volume, it is recommended to consider them for larger builds. In that case, even only the predominant losses are considered, with conduction and ablation being frequently ignored, and the heat transfer coefficients for convection and radiation are grouped together as [134]: ℎ=(ℎconv +ℎrad).(3.20) The heat loss flux is then modelled according to Newton’s Law of cooling [44, 134]: −𝑘𝜕𝑇 𝜕𝑛 =ℎ(𝑇𝑠−𝑇∞),(3.21) where −𝑘𝜕𝑇 𝜕𝑛 is the convective heat flux, with 𝑛 being the surface’s unit normal vector, ℎ is the integrated heat transfer coefficient and 𝑇𝑠 and 𝑇∞ the surface and environment temperatures, respectively. The 41 3.1. Thermal Model values for ℎconv can be determined as a function of the Nusselt number (𝑁𝑢): ℎconv =𝑁𝑢 𝑘 𝐿,(3.22) where 𝑘/𝐿 represents the conductive heat transfer, with 𝐿 being the characteristic length of the specimen. If the scanned, molten zone is considered a horizontal plate, 𝑁𝑢 itself can be determined by [135]: 𝑁𝑢 =0.54𝑅𝑎0.25,(3.23) where 𝑅𝑎 is the Rayleigh number that is given by the product of the Prandtl (Pr) and Grashof (Gr) numbers (𝑅𝑎 =𝑃𝑟 𝐺𝑟), which can be calculated as [135]: 𝑃𝑟 =𝐶𝑝𝜇 𝑘,(3.24) 𝐺𝑟 =𝑔𝜌𝛽(𝑇𝑠−𝑇∞)𝐿3 𝜇2,(3.25) with 𝑔 being the gravitational acceleration, 𝜇 the viscosity and 𝛽 the volumetric expansion coefficient. All the variables present in these equations are for the fluid properties only. In a different manner, the heat flux due to radiation is described by the Stefan-Boltzmann Law [134]: 𝑞rad =𝜀 𝜎𝑇𝑠4−𝑇∞4,(3.26) with 𝜀 being the surface emissivity and 𝜎 the Stefan-Boltzmann constant with the value of 5.6710−8 W/m 2 K 4 . The radiation is then linearised and treated as an effective heat transfer coefficient so it can be utilized in Equation (3.20) and fit into Equation (3.21) format [134]: 𝑞rad =ℎrad(𝑇𝑠−𝑇∞),(3.27) ℎrad =𝜀 𝜎(𝑇𝑠+𝑇∞)𝑇𝑠2+𝑇∞2.(3.28) It is important to note that there is no energy loss to the air in the beginning of the scanning process, right after the pre-heating, because the air was also pre-heated, so there is a thermal equilibrium. The 42 3.3. Continuity Equation major aforementioned heat losses from the material only start to happen after the laser starts to scan the powder bed and the system thermal equilibrium is disturbed. 3.2 Linear Momentum Most works on the simulation of SLS processes only consider the temperature profiles, sintering depth and phase change, whereas mass transfer is neglected [ 28 ]. However, sometimes researchers want to study the sintering evolution between scanned particles or, in the case of SLM, motion effects caused by differences in density in the liquid phase. In these cases, it is imperative to consider the momentum, that is given by Newton’s Laws of motion, and assure its conservation. Compared to energy or mass, momentum is harder to deal with because it is represented as a vector, meaning that it has both magnitude and direction. Considering that, the governing equation for the momentum conservation is given as [124]: 𝜕(𝜌𝒖) 𝜕𝑡 + ∇ · (𝜌𝒖⊗𝒖)=∇ · (𝜇∇𝒖)− ∇𝑝+𝜌𝒈+𝐹, (3.29) with 𝜇 representing the viscosity and 𝑔 the gravitational acceleration. 𝐹 accounts for extra forces that are considered depending on the modelling approach. For metals, these additional forces generally comprise buoyancy, the Darcy’s term and surface forces including capillary forces (surface tension), thermo-capillary forces (Marangoni effect) and recoil pressure [ 136 , 137 ]. Again, for SLS with polymer applications, due to the high viscosity and consequent reduced flux, the model is generally simplified and 𝐹only accounts for surface tension, that is modelled as [138]: 𝐹𝜎=𝜎𝜅 ∇𝛼, (3.30) where 𝜎 is the surface tension constant value, 𝛼 the phase fraction and 𝜅 is the mean surface curvature, that can be calculated as [138]: 𝜅=−∇ · 𝒏=−∇ · ∇𝛼 ∥∇𝛼∥.(3.31) 43 3.4. Powder Deposition 3.3 Continuity Equation The continuity equation, also known as mass conservation equation, describes the conservation of mass in the problem. The change of mass over time, represented as 𝑑𝑚 𝑑𝑡 , refers to the mass that can flow into and out of the control volume. Physically, volume can also change overtime, but considering that the control volume is constant, mass is related to it by density. This denotes that the change in mass depends on the density and on the flux of mass through the volume as the following equation describes: 𝑑𝑚 𝑑𝑡 =𝜕𝜌 𝜕𝑡 + ∇ · (𝜌𝒖).(3.32) The models for the SLS process generally consider the involved fluids to be incompressible in the calculation process, therefore, because the density is constant, the term 𝜕𝜌 𝜕𝑡 is equal to 0 and the continuity equation only depends on the sum of the converging/diverging fluxes that, too assure mass conservation, must also be 0: ∇ · 𝒖=0.(3.33) 3.4 Powder Deposition The powder deposition is one of the main steps of the process. It occurs once for every layer and is the first action towards its fabrication. In Section 2.4.2, powder properties such as their shape, size and size distribution were associated with the quality of the powder bed, more specifically, its density and absence of voids. Since the powder bed is the foundation of the process, it is imperative to assure that it stays consistent and of good quality throughout various processes. Even with the other parameters optimized, if the base of the process is defective, so will be the final part. When simulating the SLS process, some approaches consider the powder layers as a continuum domain with mixed properties between material and air, according to the expected powder bed density [ 43 ]. Modelling the layer this way is much easier, but it only provides information on the macroscopic scale. Nonetheless, it is a much simpler way to determine temperature profiles and sintering/melting depth. However, this approach, due to its excessive geometry simplification, does not provide any information on phenomena that happens at particle length scale, such as porosity and 44 3.4. Powder Deposition sintering evolution. More useful data could be obtained if the powder bed was modelled as realistic as possible, but this requires the mesh to be much more refined, in order to accommodate many small particles, and the coupling of different simulation methods. The powder bed is considered a granular material, more specifically, a group of numerous discrete solid particles that interact with each other through collisions and/or other interaction types dispersed over a domain. From a macroscopic scale, a granular material cannot be assigned to a state type, such as solid or liquid, because none can represent the complete behaviour of the group of particles. Despite the particles being solid, because of their small size, when under certain conditions, such as being dragged by a roller over the powder bed, they behave similar to a fluid, hence researchers studying powder flowability (see Section 2.4.2.2). Moreover, under highly agitated systems, groups of particles can behave as a gas, as seen daily with dust particles in the air. That is why, when grouped together, small solid particles form a granular material. Modelling the interaction of the particles during their deposition over the powder bed cannot be done with mesh based methods, such as the ones mentioned in the previous sections, and, instead, a meshless approach is employed in the form of a suitable numerical method. The Discrete Element Method (DEM) is widely used for particle simulations. In fact, it was developed in 1979 by Cundall [ 139 ] to simulate the behaviour of granules such as soil. DEM discretizes granular media in a way where each particle is a single element and then calculates the interactions between them at a microscopic level, so one can observe the influence of such interactions at a macroscopic level. This method requires a contact detection strategy, interaction Laws and time discretized equations of motion that govern the particle displacements. The motion of individual particles is governed by the Newton-Euler Laws of motion, while the interactions between them are calculated using a variety of models suitable for a particular geometry and material behaviour, considering the necessary forces. As in any simulation, there are numerous models available [ 140 ] and their selection may vary according to process conditions, materials and/or the complexity desired. Generally, DEM approaches consider translational and rotational motion and the governing equations are given by [141]: 𝑚𝑑2𝑷 𝑑𝑡2=∑︁(𝑭𝒏+𝑭𝒕+𝑭𝒄) + 𝑚𝒈,(3.34) 𝐼𝑑𝜔 𝑑𝑡 =∑︁(𝑇−𝑅𝑓),(3.35) where 𝑚 is the mass of the particle, 𝑷 the position vector, 𝑭𝒏 and 𝑭𝒕 are the contact forces in 45 3.4. Powder Deposition the normal and tangential direction, respectively, 𝑭𝒄 the cohesion forces and 𝒈 is the gravitational acceleration. Moreover, 𝐼is the inertial moment, given by 𝐼=0.4𝑚𝑟 [141], 𝜔is the angular velocity, 𝑇 is the torque that the particles exert on each other due to contact forces and 𝑅𝑓 is the rolling friction. In Figure 3.2, a schematic representation of the force diagram, with the relevant forces for the powder bed simulation, for two contacting particles with radii 𝑟1 and 𝑟2 and masses 𝑚1 and 𝑚2 , is illustrated. The torque 𝑇from other contacting particles is given by [24]: 𝑇=𝑟𝐹𝑡,(3.36) and the rolling friction is calculated as [141]: 𝑅𝑓=𝜇𝑟𝑭𝒏𝜔, (3.37) where 𝑟 is the particle radius, 𝜇 is the friction coefficient for the sliding friction and 𝜔 is the angular velocity unit vector. The cohesion forces 𝑭𝒄 are based on the Johnson-Kendall-Roberts (JKR) cohesion theory. They represent attractive forces due to van-der-Waals effects and are essential for a correct powder bed simulation because of the particles small size, as discussed in Section 2.4.2.2. They are implemented as 𝑭𝑱𝑲𝑹 [142]: 𝐹𝐽𝐾𝑅 =−4√︄𝜋𝑎3𝐸𝜆 2(1−𝜉2) −→ 𝑒𝑛,(3.38) where 𝑎 is the contact radius, 𝐸 the Young’s modulus, 𝜆 is the surface energy density and 𝜉 represents the Poisson’s ratio. The contact forces 𝑭𝒏 and 𝑭𝒕 are calculated based on the Hertz-Mindlin contact theory as [143]: 𝑭𝒏=(𝜅𝑛𝛿3/2 𝑛+𝜁𝑛𝒗𝒏)−→ 𝑒𝑛,(3.39) 𝑭𝒕=min(𝜅𝑡𝛿𝑡+𝜁𝑡𝒗𝒕, 𝜇𝑭𝒏)−→ 𝑒𝑡,(3.40) with 𝛿𝑛 , given by 𝛿𝑛=𝑟1+𝑟2−𝑟12 , where the last term is the distance between the centre of the particles, representing the compression of colliding particles in the normal direction and 𝛿𝑡 the deformation between particles in the tangential direction at the contact points. The velocity is represented as 𝑣𝑛 and 𝑣𝑡 and the unit vectors as −→ 𝑒𝑛 and −→ 𝑒𝑡 , both in the normal and tangential direction, 46 3.5. Finite Volume Method ∫𝑡𝑛 𝑡𝑛−1 𝜙(𝑡)𝑑𝑡 =1 2(𝜙𝑛−1+𝜙𝑛)Δ𝑡, (3.58) where 𝜙𝑛−1 denotes the value from the previous time step, 𝑡𝑛−1 , and 𝜙𝑛 the value at the time step being solved, 𝑡𝑛 . By applying these notions to Equation (3.56) , it is finally fully discretized in the popular Crank-Nicholson form [147]: 𝜌𝑛 𝑃𝜙𝑛 𝑃−𝜌𝑛−1 𝑃𝜙𝑛−1 𝑃 Δ𝑡𝑉𝑃+1 2∑︁ 𝑓 𝐹𝜙𝑛−1 𝑓+𝜙𝑛 𝑓−1 2∑︁ 𝑓Γ𝜙𝑓𝑺𝒇·(∇𝜙)𝑛−1 𝑓+(∇𝜙)𝑛 𝑓 =𝑆𝐸𝑉𝑃+1 2𝑆𝐼𝑉𝑃𝜙𝑛−1 𝑃+𝜙𝑛 𝑃. (3.59) System of Equations From the analysis of the above equation, one can conclude that, in order to solve the problem and determine the present values of 𝜙 at the faces, the contributions from all the terms are required as well as the values of 𝜙 from the previous iteration. That leads to an algebraic equation with the following form for each CV [146, 147]: 𝑎𝑃𝜙𝑛 𝑃+∑︁ 𝑓 𝑎𝑁𝜙𝑛 𝑁=𝑆𝑃(3.60) with 𝑎𝑃 representing the coefficient for the CV of interest, 𝑎𝑁 the the coefficient for the neighbour cells and 𝑆𝑃 the source terms. By assembling the equations for all computational cells, a system of equations is obtained in the form of: [A]{𝜙}={b}(3.61) where {𝜙} contains the unknown values that will be calculated, [A] is a sparse matrix composed of the the 𝑎𝑃 coefficients on the diagonal, and the coefficients of the neighbour cells, 𝑎𝑁 , on the off-diagonal positions. Lastly, {b} represents the source terms. By solving this system of equations, the solution of the mass, momentum and energy conservation equations can be obtained and the main variables values are calculated. 53 CHAPTER 4Computational Tools After reviewing the modelling of the most important important phenomena in the SLS process, a set of computer based techniques, also known as computational tools, are necessary to analyse and solve the equations. For the problem in hands, multiple tools are required, since different methods need to be employed. The different tools for each method, along with their process compatible models, will be described in the following sections. 4.1 OpenFOAM Open-source Field Operation And Manipulation or, in short, OpenFOAM [ 148 ], is a free open-source software used to solve a multitude of complex CFD problems. It is composed of a vast library of models and equations written in C++ programming language, which has been used and greatly expanded throughout the years by a growing community of contributors across most areas of engineering and science. Since 2004, it has been mainly developed by OpenCFD, with two major releases every six months that include customer sponsored developments and contributions from the community. The software is compatible with most Linux distributions, mainly Ubuntu, but it can also be used on other major operative systems, such as Windows and macOS, through the Docker [ 149 ] or WSL [ 150 ], that provides a self-contained Linux subsystem with a compatible environment for it. Since OpenFOAM base version lacks a graphical user interface, executing commands is done through the command prompt and the visualization of the data is done through post-processing with external software, for example, ParaView [151]. 54 4.1. OpenFOAM 4.1.1 General Case Structure Each OpenFOAM case consists of a directory structure that contains specific files, stored in specific folders, that are required to run the simulation. Because OpenFOAM is a toolbox that contains many solvers, according to the application, particular files might be needed. Figure 4.1 represents the general structure of a case. Case 0 U p T alpha system controlDict fvSchemes fvSolution constant polyMesh thermophysical Properties turbulence Properties g Figure 4.1: General case directory structure. Time directories The time directories contain the solution fields for specific time-steps and are generated during the simulation, apart from the initial data (usually in the 0 folder), which contains the initial and boundary conditions for the simulation. The variables are specified using the International System of Units (SI) and are displayed as SI [kg m s K mol A cd]. Typically, this folder contains information regarding the velocity ( 𝑈 ), pressure ( 𝑝 ), temperature ( 𝑇 ) and phase (alpha). For multiphase problems, each phase requires an alpha file. The simulation can also be started from a subsequent time. 55 4.1. OpenFOAM Constant directory This directory contains files that allow characterization of the problem physics, as well as the computational mesh. The mesh is stored in the polyMesh folder, distributed between multiple files that contain information on the points and the faces that they form when orderly connected, as well as the boundary faces that are generally grouped into patches used when specifying boundary conditions. Most solvers require additional information on gravity ( 𝑔 ), the presence or absence of turbulence and the thermophysical properties of the involved phases, such as density and specific heat. System directory The system directory contains files that dictate how the problem should be solved. The three included in Figure 4.1 are mandatory to every case. The controlDict file is for the specification of more general settings, such as the simulation time-step, start and end time, how often a time-step solution should be written to the case directory, what format the data files have and, for transient solvers, the maximum Courant number. The other two files are for controlling simulation options related to the methods. fvSchemes contains the necessary information related to the discretization schemes for the governing equations differential terms, namely how each equation is integrated with respect to time, how the gradient of each field is calculated, the discretization of the divergence, Laplacian and surface-normal gradient terms and how the interpolation from cell-centred values to face-centred values is computed. fvSolutions is where the user specifies the controls for the solution and algorithm, such as residual control. Another optional files consist of mesh generation dictionaries, such as blockMeshDict, subdivision of the domain for parallel processing with the decomposeParDict and, very frequently used, setFieldsDict, which allows the user to set a region anywhere in the domain with a specific value. 4.1.2 icoReactingMultiphaseInterFoam Solver The icoReactingMultiphaseInterFoam solver was originally released with OpenFOAM version 1806 [ 152 ]. It is a subcategory of the multiphaseInterFoam solver that, itself, is based on the parent solver interFoam, which solves the continuity and momentum Navier-Stokes equations for two fluids with a phase fraction that is tracked resorting to the volume of fluid algorithm [ 153 ]. multiphaseInterFoam follows the same principles, but for multiple incompressible fluids, tracking phase fraction, while also 56 4.1. OpenFOAM considering surface tension and contact-angle effects for each phase [ 154 ]. Being a subcategory of these solvers, icoReactingMultiphaseInterFoam shares many features, but was complemented with additional functionalities. For instance, it allows the user to chose the thermodynamic model for each phase, but the main assumption of the shared fields among the phases from interFoam, namely velocity, pressure and temperature, remain. The solver supports mass and heat transfer between different phases of the same material, similar to the ones on which it is based, but also allows it between different phases of the different materials. Among the many solvers present in the OpenFOAM library, the icoReactingMultiphaseInterFoam was selected because it meets the modelling requirements for the SLS process presented previously, such as modelling the polymers phase change, considering fluid flow and taking into account interactions between multiple phases, including interface tension forces. The solver also contains a laser radiation model capable of modelling a volumetric heat source, following a Gaussian distribution, that attenuates as it passes through the material. Previous work was also done to explore and assure the solver capability to simulate the SLS process [ 60 ], which further reinforces its selection to model the problem of interest. 4.1.2.1 phaseProperties Being a multiphase solver that accounts for multiple interactions between the different phases, a dictionary for the specification of the models for certain interactions alongside their parameters is needed. For this solver, that dictionary is located in the constant folder and is named phaseProperties. The available models that can be indicated on it are described bellow. Phase models The solver comprises three basic phase models:  pureMovingPhaseModel - adequate for flowing phases (fluids).  multiComponentMovingPhaseModel - used for multi-component fluids (fluid mixtures)  pureStaticSolidPhaseModel - meant for solids. For the SLS process modelling, the adequate models should be pureStaticSolidPhaseModel for the powder particles in solid state and pureMovingPhaseModel for regions where the temperature 57 4.1. OpenFOAM surpasses a defined value, transforming the solid into a liquid. Selecting a model for a phase in this dictionary will later dictate what properties are necessary to provide for each phase, for example, when specifying a phase as a solid, the viscosity is not requested. Mass transfer models Mass transfer models should be selected when two previously specified phases can be transformed into each other due to phase change. The solver has only two options, the kineticGasEvaporation model, used for condensation and evaporation, and the Lee model, for melting and solidification applications. Since the polymer particles can be either in solid or liquid state, melting and solidifying as the process evolves, the relevant model is the latter. The Lee model [ 155 ] assumes that the phase change occurs at a constant pressure and at a quasi-thermal equilibrium state and that it mainly depends on the difference between the interfacial cell and saturation temperatures. Because the Lee model is used more often, in the literature, for condensation and evaporation, the following equations are slightly adapted so it is easier to understand them in this context. However, the original ones can be found in some recent studies [156--158]. For melting, when 𝑇 > 𝑇activate, the mass transfer rate from solid to liquid 𝑑𝑚𝑠𝑙 𝑑𝑡 is: 𝑑 𝑚𝑠𝑙 𝑑𝑡 =𝐶𝜌𝑙𝛼𝑙 𝑇−𝑇activate 𝑇activate (4.1) and for solidification, when 𝑇 < 𝑇activate, the mass transfer rate from liquid to solid 𝑑𝑚𝑙𝑠 𝑑𝑡 is: 𝑑 𝑚𝑙𝑠 𝑑𝑡 =𝐶𝜌𝑠𝛼𝑠 𝑇activate −𝑇 𝑇activate .(4.2) In both equations, 𝜌𝑙 , 𝜌𝑠 , 𝛼𝑙 and 𝛼𝑠 are the density and phase volume fraction for the liquid and solid material, respectively, 𝑇activate is the activation temperature, commonly referred to as saturation temperature ( 𝑇sat ) and 𝐶 is a constant value that represents an empirical coefficient, called the mass transfer intensity, which is used to adjust the phase change rate. Some researches have found that small values of 𝐶 result in a significant deviation of the interface temperature from 𝑇sat and that increasing this value helps reduce that. Used values range from 0.1, for flow boiling, to 7.5x 105 , for micro-channel condensation. Such information and more on this topic, with several references, can be found in Lee et al. [ 158 ]. In the solver adaptation of the model, 𝐶 also works as a switch, meaning that if it has a positive value, only melting can occur at that temperature, otherwise, there is only 58 4.1. OpenFOAM solidification. Both phase change directions are possible if both are individually specified. Inter-phase porosity models The inter-phase porosity models purpose is to add an artificial momentum source on or next to the interface between the two phases, which can influence the behaviour during phase change, mainly for solidifying or melting. The solver has only one model available for said function, named VollerPrakash after its creators [ 159 , 160 ]. This model aims to solve the major problem of assuring a null velocity as the liquid region turns into solid. The approach is to treat the cells in which phase change is occurring as pseudo porous media, where the porosity ( 𝜆 ) ranges between 1 and 0 for a fully liquid or solid cell, respectively, with an intermediate value during phase change. Most applications of convection-diffusion phase change numerical methods have been made to isothermal phase change problems, assuming that the material is pure and has a very well defined phage change temperature. Despite being a general method that can handle such situation, the model was developed to also consider the "mushy" region created during the phase change of materials, such as alloys and polymers, that, for example, melt over a temperature range, which implies that the evolution of latent heat has a functional relationship with temperature (Δ𝐻=𝑓(𝑇)) , as opposed to pure materials where there is a step transition. Modelling the "mushy" zone is more challenging since it represents a solid-liquid mixed state of the material. The VollerPrakash implements two source terms in the momentum balance equation, 𝑆𝒖 and 𝑆𝑏. The first term is given by [159]: 𝑆𝒖=−𝐴𝒖,(4.3) with 𝒖=(𝑢, 𝑣, 𝑤)representing the velocity vector in the "mushy" zone, governed by the Darcy law: 𝒖=−𝜅 𝜇∇𝑃, (4.4) where 𝜅 is the permeability, given as a function of porosity by 𝜆=1−𝐹𝑠 , with the latter term being the solid fraction, whereas 𝐴 works as an intensification factor for the velocity calculated in Equation (4.4) . For an isothermal problem, where the porosity approach is typically a numerical "fix", any increasing function for this term would be suitable, however, for a "mushy" zone, where a porous zone indeed exists, one should resort to physics to assure a proper form for the 𝐴 term. The well-known 59 4.1. OpenFOAM Carman-Koseny equation [ 161 ], derived from the Darcy Law, shows the relation between the cell gradient of pressure, porosity and velocity: ∇𝑃=−𝐶(1−𝜆)2 𝜆3𝒖.(4.5) This equation suggests that 𝐴should be represented as [159]: 𝐴=−𝐶(1−𝜆)2 (𝜆3+𝑞),(4.6) with 𝑞 being a constant, which takes very small values, introduced in the equation just to avoid division by 0. The value of 𝐶 , another constant, depends on the morphology of the porous media. 𝐴 will increase from 0 to large values as 𝐹𝑠 in the cell increases from its lowest liquid value of 0 to its highest fully solid value of 1. This defines the source term as negligible in the liquid region and the momentum equation is given in terms of the actual fluid velocities. However, as phase change stars to occur and the "mushy" zone is formed, the source term’s influence increases along with 𝐴 . The values of 𝐴 increase in such a way that the value of the sources overpass the transient, convective, and diffusive terms contributions and the momentum balance equation approximates the Darcy Law (Equation (4.4) ). Following this equation, as 𝐹𝑠 approaches 1, the velocity will tend to 0 and along high values of 𝐴 , this source term wil dominate all other terms in the momentum balance equation and force the predicted superficial velocities to approach zero. The second source term added by the model, 𝑆𝑏 , is a buoyancy term implemented to induce natural convection and it adds momentum only in the vertical direction. It assumes that the Boussinesq approximation is valid, which means that the density is constant in all terms except a gravity source term, and is described as [159]: 𝑆𝑏=𝜌𝑔𝛽(ℎ−ℎref) 𝐶𝑝 ,(4.7) where 𝛽is the thermal expansion coefficient, ℎis the sensible heat and ℎref its reference value. Surface Tension Models To include surface tension forces, the user must specify the value for a pair of phases. The solver is only capable of treating the surface tension as a constant continuum surface force. The equation that 60 4.1. OpenFOAM models this phenomena is based on the general approach for considering surface tension forces, which is illustrated in Equation (3.30). 4.1.2.2 Radiation model As mentioned previously, the availability of an adequate radiation model was one of the motivations to select this solver. It uses the Discrete Transfer Radiation Model (DTRM) to simulate a collimated radiation flux that allows for interaction between the laser and the participating media [ 152 ]. DTRM based methods work by discretizing the laser into multiple rays that will interact with the domain. Using such methods allows for the consideration of reflection and refraction of the individual rays, however, one should find a balance between the number of rays and their respective interactions, as ray tracing tends to be computational expensive. In the laser dictionary, radiationProperties, it is possible to adjust the available laser parameters, as well as some additional models. Within the available parameters, there is the laser power, expressed in Watts, direction, indicated by a vector, and position, which can be either constant or vary with time to simulate the laser target movement. The shape of the beam is a circular disk with a user defined radius and the distribution of rays inside the laser periphery is given by a polar division of that disk, with 𝑛𝑇 ℎ𝑒𝑡𝑎 number of particles in the theta direction, and 𝑛𝑟 particles along the radial direction. The power distribution of the laser can be either manually defined, uniform, or, similar to the reality, have a Gaussian distribution. It is also possible to define how the laser interacts with each phase, with the localDensityAbsorptionEmission model, allowing for different phases to have different emission and absorption coefficients. Both the power distribution and emission coefficients are treated similar to the common approaches in literature, as described in Section 3.1.2. The power intensity follows the same formula as Equation (3.12) and the user defined radius is not used directly into it, but instead to truncate the Gaussian distribution. In this dictionary, the user can also select a model to account for the reflection. One of them is named Fresnel and has similar but more complex approach than the commonly employed models in the literature, which were described in Section 3.1.2. This model calculates the reflectivity by decomposing the wave into two linear components that are vibrating within the plane of incidence 𝜌𝑃 and its normal 𝜌𝑁, represented by [162]: 𝜌𝑁=(𝑛1cos𝛼−𝑝)2+𝑞2 (𝑛1cos𝛼+𝑝)2+𝑞2,(4.8) 61 4.2. LIGGGHTS 𝜌𝑃=(𝑝−𝑛1sin𝛼tan𝛼)2+𝑞2 (𝑝+𝑛1sin𝛼tan𝛼)2+𝑞2𝜌𝑁,(4.9) where 𝑛1 is the refractive index of the medium of the incident ray and 𝛼 is the angle between the incident ray and the surface normal. The total reflectivity 𝜌 is given by (𝜌𝑁+𝜌𝑃)/2 . Lastly 𝑝 and 𝑞 are calculated as [162]: 𝑝2=1 2"√︂𝑛2 2−𝑘2 2−𝑛2 1sin2𝛼2 +4𝑛2 2𝑘2 2+𝑛2 2−𝑘2 2−𝑛2 1sin2𝛼#,(4.10) 𝑞2=1 2"√︂𝑛2 2−𝑘2 2−𝑛2 1sin2𝛼2 +4𝑛2 2𝑘2 2−𝑛2 2−𝑘2 2−𝑛2 1sin2𝛼#,(4.11) with 𝑛2 being the refractive index of the remaining medium and 𝑘2 is the absorptive index, also known as the imaginary part of the index of refraction. This implies that, instead of just mentioning the refractive index of each medium, for the absorption medium, one has to specify the complex index of refraction 𝑚in the form of 𝑚=𝑛2+𝑖𝑘2. 4.2 LIGGGHTS LIGGGHTS [ 56 ] is an open-source Discrete Element Method particle simulation software whose code is written in C++ programming language. The software’s name is an acronym that stands for LAMMPS Improved General Granular and Granular-Heat Transfer Simulations, where LAMMPS is a widely used simulator in the Molecular Dynamics (MD) field. LIGGGHTS uses LAMMPS as a platform, making improvements on it by moving from MD to DEM through the addition of characteristic features of DEM such as simulation of contact and cohesion forces, rolling friction and heat conduction between particles. LIGGGHTS is being developed with the goal of being used in industrial applications, hence the possibility to import complex CAD geometries, utilize moving meshes and include a variety of particle insertion options. 4.2.1 Input Script Contrary to OpenFOAM, LIGGGHTS reads all simulation parameters from a single input script, all written by the user. The input file must follow a specific structure because the software will read from 62 5.1. Geometry and position were stored and the visualisation was performed with Paraview. For the simulation, two geometries, a box and a blade, that are displayed, together with their dimensions, in Figure 5.1.a, were created using a CAD software. Both act as walls, which means that particles collide with them without getting through. The box has a deeper section that can be viewed as a smaller box under the larger one, where the particles will fall, simulating the beginning of the deposition of one layer. This smaller box will contain the representative section of the powder bed. The first step in the simulation is the particles insertion. To ensure that the smaller box is totally filled, the selected amount of particles was large enough so that when the blade moves past it, it still drags a significant amount of powder. In this case, the amount of particles is given by mass and the selected value that fulfils these requirements (without inserting an exaggerated large amount of particles) was 0.115 g. The insertion phase lasts for two seconds and takes place over a area in the far end of the larger box (Figure 5.1.b). It is then followed by the blade movement, that moves at a constant speed in the 𝑥 direction. Before the inserted powder reaches the smaller box, it has to be dragged over a long flat section that was purposely made in such way, not only to accommodate the particles insertion, but also to assure a more realistic powder distribution because, this way, the particles are packed by the blade movement and their "flow" is more developed before falling into the hole (Figure 5.1.c). Once past the smaller box, the blade will keep moving, dragging the remaining particles past the boundary, outside the simulation domain, deleting them (Figure 5.1.d). At the end, only the particles inside the smaller box will be used for the mesh generation. The particle interaction models for the powder bed simulation, which are described in Section 4.2, were selected by comparing the models available in the software with the literature on the topic. The software also requires the input of various material properties that are indicated on Table 5.1 [ 168 ]. The powder particle size distribution was also taken into account. In LIGGGHTS, each individual particle size has to be specified with its corresponding percentage in mass or quantity. Because of the impossibility to consider and infinite number of particle diameters, or the impracticality of specifying a large amount of them, nine different particles sizes, centred on the highest percentages, starting at 35 µ m of diameter and increasing by 5 µ m, until the maximum diameter of 75 µ m was reached, were defined. Each individual size was represented alongside its respective percentage, in order to approximate the real material distribution that is illustrated in Figure 5.2 [ 169 ]. Lastly, for simplification purposes, and also due to the impossibility of manually defining an unique shape for each particle, they were considered to be perfectly spherical, as commonly done in similar studies [47]. 69 5.1. Geometry 4000μm 1000μm 400μm 200μm 500μm (a) (b) (c) (d) Figure 5.1: Geometry dimensions and various stages of the simulation. Initial time (a), end of particle insertion (b), particles falling into the representative section (c) and end of the simulation (d). Table 5.1: Polymer powder properties for the simulation in LIGGGHTS. Parameter Value Material Density 1000 kg/m3 Young’s Modulus 2.3x107Pa Poisson’s Ratio 0.4 Coloumb’s Friction Coefficient 0.5 Surface Energy 0.1 mJ/m2 Volume (%) Diameter (μm) Particle Size Distribution PA12 - PA2200 12 10 8 6 4 2 020 40 60 80 100 Figure 5.2: Particle size distribution of the selected material for the simulation. 70 5.2. Computational Mesh Finally, to validate the results, the bulk density of the simulated representative section was calculated using Paraview. The obtained value was 0.43 g/cm 3 , which is satisfactorily comparable to the material data sheet value of 0.45 g/cm3[170]. 5.2 Computational Mesh The representative section of the powder bed simulated using LIGGGHTS needs to be converted into a computational mesh before it can be used by the CFD software. Multiple approaches can be employed for that purpose, but, considering the large number of particles and consequent complexity, some could reveal themselves too complicated or time-consuming, whereas the following approach is simple and practical. Using OpenFOAM tools, a meshed box with dimensions of 500x1000x230 µ m (slightly taller to accommodate particles that exceeded the top "limits" of the box) was generated using blockMesh [ 171 ]. Then, the particles in the last time step of the LIGGGHTS simulation needed to be converted into a STL file so they can be used for the mesh generation step. Since the particles are perfectly round, a spherical glyph, that scales with the radius, was created around each particle centre point in Paraview, allowing the simulation data to be converted into a geometry file, in STL format. Finally, with the setFields utility [ 172 ] from OpenFOAM, the STL file was inserted into the meshed box and the cells that intersected it were defined as solid. By using blockMesh, a fully orthogonal mesh with equal sized cells was obtained, while each particle and atmosphere regions are properly set. When creating a computational mesh, it is important to balance the accuracy with the simulation time, since both will increase with the cells number. For this specific case, the limiting factor are the smaller particles that require a higher resolution (smaller cells) to be adequately represented. Four different refinement levels were tested, starting with the coarsest mesh, that has one cell per 8 µ m, and then doubling the number of cells (when possible) in each direction for the next refinement levels, as shown in Table 5.2. The Coarse mesh (Figure 5.3.a) is clearly inadequate since the geometry is not properly represented. The resolution is so low that the particles barely resemble a sphere, with the smaller ones fitting inside four cell sided cubes. The Medium mesh (Figure 5.3.b), with double the refinement level, already provides a good geometry representation, with well defined spheres and an acceptable definition for the smaller particles. For preliminary simulations, where the main goal is to quickly 71 5.3. Material Properties Table 5.2: Number of cells for the different levels of mesh refinement. Name Cells x-direction Cells y-direction Cells z-direction Total Cells Coarse 63 29 125 228 375 Medium 125 58 250 1 812 500 Fine 250 115 500 14 375 000 Extra Fine 300 138 600 24 840 000 test some parameters, this mesh is sufficient. For more detailed studies, however, the Fine mesh (Figure 5.3.c) is significantly more adequate. Despite the large number of cells, which will increase the computational cost, the geometry is remarkably well defined, with the particles looking almost spherical. A good refinement level is always desirable since it greatly reduces the error associated with the simulation. Moreover, there will be movement in the liquid phase, caused by surface tension forces, that is expected to be substantially better observed with this mesh. Despite the satisfactory results, an additional refinement was tried. Increasing the refinement further than this is very challenging in terms of computational power since a large amount of memory is required. Due to memory limitations, this time, the increase could not be doubled and, instead, using the available hardware, an increase in 20% of cells in each direction was tested (Figure 5.3.d). As expected, the particles are even closer to their original shape in the STL file, however, considering that the refinement level obtained in the previous attempt was already sufficient to properly represent the geometry and that the Extra Fine mesh has ten million more cells, it was discarded for the sake of computational time. Among the four meshes, the Fine mesh, with one cell per two micrometers and a total of around fourteen million cells, was the selected to perform the studies. The approach for selecting the mesh is only based on the balance between their geometry representation quality and the number of cells, which is not totally correct, since various simulations with each mesh should have been performed to assess the refinement level influence on the results. However, due to time limitations, such approach was not possible and the next best option was employed. 5.3 Material Properties The selected OpenFOAM solver, contrary to other popular solvers, considers each state of the material as a separate phase and requires the specification of its properties separately. This translates 72 5.3. Material Properties (a) (c) (b) (d) Figure 5.3: Visualisation of the different mesh refinement levels. Coarse mesh (a), Medium mesh (b), Fine mesh (c) and Extra Fine mesh (d). into some properties, such as viscosity, only being required for one of the phases, in this case, the liquid one. For this specific case, a total of three phase files should be provided: alpha.solid, for the initial solid polymer, alpha.liquid, for the polymer in liquid state and alpha.air, for the surrounding air present in the chamber. 5.3.1 Polymer The material utilized in the process used for comparison purposes is a widely used polyamide for SLS applications, PA12 - PA 2200, produced by EOS [ 173 ]. According to the manufacturer, this PA12 is a high-performance equivalent to well regarded materials in injection moulding, such as ABS or PA6, capable of producing equally strong, flexible and durable parts. Accurately defining each property for polymeric materials is challenging because of their intrinsically complex behaviours. Most properties, in order to be accurately determined, require the replication of the SLS process conditions that are hard to achieve, or obtained through unconventional characterization methods, which leads to less adequate information available. Taking the viscosity as a first example, the test conditions must be as close as possible to null shear rate to capture the real value without the shear thinning influence, which requires specific rheometers. Another example is the melting point whose values are obtained through DSC tests, generally at a constant heating rate. 73 5.3. Material Properties However, in the real process, the material is only going to heat at such rates in the pre-heating phase, where no melting occurs. After that, it heats much faster by action of the laser. This is a less accurate way to determined the desired properties, because the rate at which polymers heat can change their melting point and consequently its enthalpy of fusion. Moreover, polymer properties heavily depend on the surrounding conditions, namely the temperature. As already stated, the increase in temperature during the process causes the material to age, a phenomena that will change its properties. Using the viscosity as an example again, even if the powder bed is considered to contain only virgin material, by the time it reaches the pre-heating temperature, due to ageing, the viscosity will no longer be the one initially defined. Similar to these examples, many other properties are challenging to obtain due to suchlike reasons, thus some simplifications and assumptions need to be considered. A common one is to assume temperature independent properties. Although the solver allows temperature dependent properties, no adequate information was found that allowed such definition for several parameters. This should not present itself as a large hindrance since that, although there is a large temperature gradient after the scan, the material will cool down to temperatures closer to the initial one relatively fast [ 174 ]. Between the list of parameters that the solver can work with, the useful ones for the simulation are the density, specific heat, thermal conductivity, viscosity, surface tension, melting point, complex index of refraction and coefficient of absorption for the laser. These parameters are indicated in Table 5.3. Starting with the density, the material datasheet only provides information on the bulk powder bed and sintered parts value, with no actual values for the density of the material itself. Therefore, a general value for PA12 of 1000 kg/m 3 was selected. The thermal conductivity was provided by the in-house process operators and the surface tension [ 102 ], absorption coefficient [ 47 ] and viscosity [ 110 ] values were taken from the literature, from articles that considered the SLS process conditions. The latter contains two values, the lower corresponds to the virgin material, and the higher, due to ageing, to the processed one. The two components of the complex index of refraction are the refractive index, with the value 1.6 [ 175 ], and the imaginary part, 𝑘 . Unfortunately, no information was found for the latter component, however, it can be related to the absorption coefficient as 𝑘=𝛼𝜆 4𝜋 [ 176 ], giving 𝑘 the value of 0.011. All the aforementioned properties were defined as constant, mainly due to the lack of information to specify them otherwise. The two remaining parameters, the specific heat and the melting point, require a different treatment for various reasons. One of them, the solver does not take into account the enthalpy of fusion, which would help into more accurately accounting the total energy required for phase change. Fortunately, in-house tests were performed on the material that allowed the determination of 74 5.3. Material Properties Table 5.3: Some of the required polymer properties for the computational studies. Property Value Units Density 1000 kg/m3 Thermal Conductivity 0.2 W/m K Viscosity (at 474k) 390/5095 Pa/s Surface Tension 0.035 N/m Absorption Coefficient 1.3x104m-1 Complex Index of Refraction 1.6+𝑖0.011 − the specific heat as function of the temperature, which is expected to greatly increase the accuracy of the simulation. Within the temperature range of interest, the values were introduced in the material properties dictionaries in the form of a table. Additionally, not considering the enthalpy of fusion brings an additional challenge for selecting the melting point. Polymers do not have a well defined melting point, but instead, melt over a range of temperatures. Lamentably, the solver does not consider such behaviour and the material will change phase when it reaches the transition temperature. Although the total energy is accounted for with the temperature dependent specific heat, it is still necessary to define a temperature for phase change. In reality, process operators tend to consider the beginning of the melting range, also known as onset of melting, as the melting point that, for this material, is 449.15 K. Doing that for the simulation is inaccurate because it would be largely over predicting the amount of liquid phase, since it would consider as fully liquid any solid cell above that temperature. Similarly, on the other end, only considering the melting point at the end of the melting range would, much likely, promote an under prediction, especially considering that it is around 469 K [ 177 ]. Due to these reasons, the selected phase transition temperature is somewhere in between those two, at the peak energy point during the DSC test, commonly referred to as simply the melting temperature. According to the literature, the melting point for fresh PA12 is at 459.75 K and it suffers a slight increase after being used, with a mid point of 460.65 K [ 175 ]. Another study, this time for the same PA2200 used in the process being simulated, exhibits a melting point of 461.75 K [ 177 ]. Despite the first study providing various temperatures depending on the material use, because there is a very small difference between all of them, the selected one was from the latter study. Although describing in Section 3.1 the different approaches for modelling the material crystallization, such phenomena was not considered. In fact, once the material liquefies, it will stay in that state throughout the entire simulation. That is because the selected domain is but a representative section that, considering the simulation time and the sintering window, will never reach temperatures 75 5.4. Model Parameter Studies low enough to solidify. Besides that, the main objective of the work is to evaluate the parameters influence on the particles melting and coalescence, which denotes the crystallization as irrelevant. 5.3.2 Chamber Atmosphere The chamber is filled with an inert gas, usually nitrogen, in order to avoid degradation reactions, such as oxidation. Even though a significant amount of information was found regarding the temperature dependency of the nitrogen properties, they were all defined as constant to simplify the model and alleviate the calculations. Besides that, the material at study is the polymer so, as long as the gas properties are somewhat correct, it should not drastically influence the results. The values for the density [ 178 ], thermal conductivity [ 179 ], specific heat [ 180 ] and viscosity [ 181 ], for temperatures within the process range, were all obtained from well-known databases. For the absorption coefficient, no information was found regarding its value for the employed laser, however, a study for other gases, with a CO 2 laser beam at a close wavelength value, showed that the absorption coefficient for gases is much smaller and almost negligible [ 182 ]. Since no study on the process mentions the interaction of the laser with the air, it was assumed that it was insignificant, therefore a small value was selected. The gas properties are displayed in Table 5.4. Table 5.4: Nitrogen properties employed in the computational studies. Property Value Units Density 0.676 kg/m3 Thermal Conductivity 3.887x10−2W/m K Specific Heat 1049 J/kg K Viscosity 2.6x10−7Pa/s Absorption Coefficient 1 m-1 5.4 Model Parameter Studies The proposed solver for the CFD simulation is extremely complex and contains a variety of parameters for tuning the models that are, most of the times, case dependant and, therefore, not available in the literature, especially numerical ones that have no physical meaning. It is necessary to evaluate these parameters influence so one can obtain results that are comparable to those of the real process to assure that there is physical meaning behind what is being simulated. Only after that is it 76 5.4. Model Parameter Studies possible to assess the process parameters influence. 5.4.1 Surface Tension Studies As already stated in Section 2.4.3.1, the surface tension is the driving force for the coalescence between particles, which makes it extremely important for the process dynamics. Such crucial mechanism has to be accounted for in the model and its behaviour needs to be studied before advancing with more complex simulations. To test the model, two spheres with 50 µ m of diameter (the average particle size) were positioned inside a domain where they are not touching anything but themselves, as shown in Figure 5.4. To isolate the studies, all optional models are disabled and the spheres already start in a liquid state at a uniform and constant temperature. Because surface tension is a force opposed by the viscosity of the material, while keeping it constant at a value of 0.035 N/m, two different viscosities were employed in distinct tests, 390 and 5095 Pa · s. The values are, respectively, from the virgin and used material. Figure 5.4: Base case geometry for the surface tension studies. Virgin Polymer The results from the tests with the virgin material, illustrated in Figure 5.5, show a complete coalescence between the two particles after five seconds, which, considering that the building process lasts much longer, is a fast consolidation. Unfortunately, no information was found regarding the sintering time between two particles that allowed for a comparison of the obtained results, however, a study on the sintering evolution between two particles of PA2200 by slowly increasing their temperature is illustrated on Figure 5.6 [ 183 ]. Despite the absence of a time scale, it is possible to analyse and compare the particles shape and the neck that they form. Figure 5.6.a shows the particles forming a very small neck, in a stage slightly 77 5.4. Model Parameter Studies 2 0 0.25s 0.50s 0.75s 1s 2s 3s 4s 5s Velocity magnitude (μm/s) 0 Figure 5.5: Surface tensions studies with virgin material. before what is visible after 0.25s on the simulation. At Figure 5.6.b, the neck is further developed and the shape similar to the results at 0.5s. Lastly, in Figure 5.6.c, the two particles are now a single one with an oval shape that resembles the results at 2s. (a) (b) (c) Figure 5.6: Hot stage microscopy of PA12 - PA2200 (adapted from [183]). Used Polymer The tests with the used material (Figure 5.7) differ significantly from the previous ones. The used material viscosity is more than ten times larger, which highly hinders the movement of the liquid phase. For these results, a different time scale had to be considered, since it took three seconds for the neck formation to be similar of just a quarter of a second with the virgin material. In this case, the much larger viscosity does not allow for a quick development of the neck and, instead, the particles coalesce partially by attracting themselves. Mixed Polymer As previously noted, the material will suffer alterations in its properties, namely an increase in viscosity, due to exposure to high temperatures. There is usually a mix between the used powder and 78 5.5. Initial and Boundary Conditions (a) (b) 446.15 498 Melting point Temperature (K) Figure 5.16: Comparison results between not using (a) and using (b) the Fresnel reflection model. 5.5 Initial and Boundary Conditions OpenFOAM, similar to other CFD software, requires the specification of the initial conditions in order to define a starting point for the simulation. Besides that, the solutions for the boundary problems are also necessary to solve the system of equations. For this specific simulation on the selected solver, the user must specify the initial and boundary conditions for the temperature, pressure, velocity and fraction of each phase. Temperature Since the performed studies aim at evaluating the process evolution during and after the laser scan, the temperature conditions are those of the bed after the pre-heating phase. At this stage, as mentioned previously, both the powder and the air are at the same temperature, usually a few degrees below the material melting point. According to the experimental reference data, the initial temperature for the entire domain is 446.15 K. Forcing all the boundaries to maintain this value throughout the simulation would be unrealistic, since it would impose a large difference between the boundary face and the cell value, especially for the regions where the laser passes through. Therefore, the lateral and bottom boundaries were defined with a null gradient condition, which should avoid such large temperature gradients, although considering these boundaries as adiabatic is still not representative of the real process. The top patch, because it is not in contact with the spheres, was kept constantly at the initial temperature, since natural convection is expected to redistribute energy to the large process 85 5.5. Initial and Boundary Conditions chamber and, therefore, the global temperature increment in the air is expected to be negligible. Velocity In the beginning of the simulation, the particles are at a solid state and stationary, therefore, a null velocity was imposed in the entire domain. As the process evolves and the liquid fraction is formed, due to surface tension and gravity forces, the liquid phase develops a velocity profile. To assure that no liquid fraction is lost through the boundaries, they can be treated as walls by imposing a null velocity, however, with this approach, the molten particles will stick to the boundaries, which is unrealistic since, in reality, there are no physical walls limiting the section in the powder bed. To prevent that, a null velocity gradient was imposed in the boundaries with the obvious exception of the bottom patch that was treated as a wall, otherwise the entire mass would fall outside the domain due to gravity. By employing these boundary conditions on the lateral and top boundaries, mass exchange with the exterior is possible. Although not entirely realistic, this approach, among the available options, is the most suitable because it does not promote infeasible results as the null velocity does. Moreover, in this specific case, the particles on the periphery of the domain will always be less representative of reality, due to the unavoidable unrealistic boundary conditions, therefore, it is more adequate to select the conditions that do not impact negatively the rest of the domain. Pressure The process takes place inside a closed air tight chamber to prevent non-inert gases from contaminating the controlled atmosphere. However, at such a small scale, the representative section is but a small amount of particles surrounded by an atmosphere in a large space. Considering that, and also the fact that the velocity boundary conditions were imposed as a null gradient, the pressure for all the boundaries was set with constant atmospheric value. The initial time is under the same conditions, therefore, the same value was imposed for the internal cells. Phase Fraction The absence of inlets in this problem simplifies the boundary conditions for the phase fraction of all phases to a null gradient. The initial conditions that represent the starting phase distribution were defined with the setFields utility. This way, the cells in the internal field are represented with 86 5.6. Assessment Studies a non-uniform list of scalars where all the cells containing the particles are defined as alpha.solid with values of 1 and the remaining cells with alpha.air values of 1. The alpha.liquid phase is also considered, however not at the initial state. 5.6 Assessment Studies After generating a computational mesh, defining the properties for each phase, studying the models for the major phenomena, selecting values for their numerical constants and providing the initial and boundary conditions, the first simulations can be performed. However, before directly assessing the parameters influence, one should first compare the process with experimental results to evaluate the solver behaviour and accuracy. Only if the results resemble those of the experimental cases, considering, obviously, that some simplification are being employed and some parameters are taken from external sources, can it be concluded that the developed methodology is adequate to study the process and its parameters influence. Conveniently, the experimental cases that are going to be used as a reference are from the same in-house process and are available in an article containing information on the process parameters as well as the respective results at a microscopic level [ 185 ]. In this article, a complex geometry with various details was produced for dimensional and geometric evaluation. Among those details, there are squared structures with ten millimetres sides and various heights. Because these structures were later analysed in the article, the representative section was assumed to be in the center one of those squares, and the resulting laser path is shown in Figure 5.17. As a note, the only relevant effect of the representative section location is the time between adjacent scans that, in practice, would translate to the energy from the first scan having more time to dissipate before the second scan hits. In these studies, the geometry will be scanned by two adjacent lines and the particles inside the scan area will be analysed. For the following studies, the material properties, model constants, mesh and simplifications will be the ones stated in the previous sections. The only differing variables between tests will be the process parameters that will be selected according to the study. Three different cases, from the article of reference, with low, medium and high energy parameters, will be used as a base for comparison. A vertical cross section of the scanned part for each case, which will be used later to compare with the simulation results, is illustrated in Figure 5.18. 87 5.6. Assessment Studies 0.5 mm 1 mm 4.75mm 4.75mm Figure 5.17: Representative section top view and laser path. Figure 5.18: Cross section of samples obtained with low (left), medium (middle) and high (right) energy parameters (adapted from [185]). 5.6.1 High Energy Case The experimental test case with high energy parameters lead to the degradation of material, as it is clearly visible in the right image in Figure 5.18, by the burned aspect of the material and large holes present in the cross section. Unfortunately, the current solver is not equipped with a degradation model, but validation can be performed by comparing the material temperature profile after the scan with the results from the in-house thermogravimetric analysis (TGA), shown in Figure 5.19. The laser scan parameters used for this test case are shown in Table 5.5. Table 5.5: Laser scan for the high energy case. Parameter Value Units Power 38.7 W Scan Speed 3000 mm/s Hatch Distance 0.3 mm 88 5.6. Assessment Studies Figure 5.19: TGA of PA2200. The simulation results, displayed in Figure 5.20, show that, with the selected parameters, there is, indeed, an excessive amount of energy. Considering that the geometry contains approximately two layers of material, the energy transferred to the material is enough to fully melt both immediately. Besides that, looking at the temperature scale, every zone in red, or a reddish tone, is above the temperature for the beginning of degradation of 574.15 K/ 300◦ C, that was taken from the TGA results in Figure 5.19, and is associated with the temperature where the degradation curve begins. Degradation seems to be more prominent at the powder bed surface, especially where the two subsequent scans overlap, as visible in Figure 5.20.b. The experimental results do not show major degradation, which indicates that the material did not reach the onset of degradation of 669.6 K and was, most likely, closer to values where the degradation rate was much less accentuated. Based on such comparison, it seems that the temperature profile obtained in the simulation is close to the experimental observations. 5.6.2 Medium Energy Case As previously discussed in Section 2.4.1.1, despite its inaccuracy, researchers commonly use Drexler’s energy density equation (Equation (2.1) ) to calculate the energy density provided by the laser. These values are then used, along with their respective experimental results, to determine whether or not more energy is needed. This information is being mentioned again because, unfortunately, the experimental results for the medium energy case do not come with their respective process parameters, 89 5.6. Assessment Studies 446.15 616 Beginning of degradation Melting point Temperature (K) (a) (b) Figure 5.20: Temperature distribution after the first (a) and second (b) laser scans for the high energy case. but the energy density calculated with the aforementioned equation is known. It is also known that the layer thickness and hatch distance were the same for all test cases so, assuming that the scan speed is also maintained, the laser power can be calculated. The resulting laser scan parameters are indicated in Table 5.6 . Table 5.6: Laser scan parameters for the medium energy case. Parameter Value Units Power 30 W Scan Speed 3000 mm/s Hatch Distance 0.3 mm Looking at the simulation results in Figure 5.21, one can conclude that the provided energy is more than enough to fully melt the scanned zone for one layer, with plenty still remaining to reheat the previous one for a good layer adhesion. Compared with the high energy case, by lowering the power by 22.5%, the maximum temperature lowered by 45 K, just a few degrees over the beginning of minor degradation. The experimental results do not display any signs of degradation, only an almost fully dense region with minor porosity from, most likely, powder bed formation defects. This implies that the material achieved high temperatures that promoted fluid flow, by lowering the viscosity, and allowed the particles to coalesce easier. Judging by the impossibility to detect layer separation, the energy provided must also have been enough to reheat the previous layers, which lead to a good layer adhesion. The numerical results, based on the temperature profile after scan, reasonably match the experimental ones. 90 5.6. Assessment Studies (a) (b) 446.15 571 Melting point Temperature (K) Figure 5.21: Temperature distribution after the first (a) and second (b) scan for the medium energy case. 5.6.3 Low Energy Case The last case, whose laser scan parameters are indicated in Table 5.7, has a significant decrease in laser power and, consequently, energy density. After analysing the experimental results, one can conclude that the diminishing energy clearly affected the built part since there is poor adhesion between some interlayer sections. Nonetheless, most of the cross section is fully dense and no unmelted particles are visible, which denotes that the lack of coalescence in the defective zones might be related to something else. Table 5.7: Laser scan parameters for the low energy case. Parameter Value Units Power 17.1 W Scan Speed 3000 mm/s Hatch Distance 0.3 mm The numerical results for the representative section temperature profile after the laser scan, visible in Figure 5.22, are in agreement with what was mentioned in the previous paragraph. There is, indeed, enough energy to effectively melt one layer and reheat a significant portion of the previous one. However, due to the notable power decrease, the melt depth and overall temperatures are now much lower. Since in this case, contrarily to the previous ones, the energy transferred to the material is insufficient to melt the entire section thickness right after the scan, the simulation was left to run until the remaining particles either melt through conduction, or the overall temperature reaches an 91 5.6. Assessment Studies equilibrium. (a) (b) 446.15 512 Melting point Temperature (K) Figure 5.22: Temperature distribution after the first (a) and second (b) laser scan for the low energy case. After just two seconds, the dissipated heat was enough to melt the bottom of the section, as observed in Figure 5.23. As a note, in reality, the average temperature of the domain would be lower due to natural convection currents that further cool the material beyond the heat dissipation promoted by conduction. These phenomena are generally considered for domains with a continuous media with a heat convection boundary condition at the surface [ 44 ], but, for these studies, that alone is insufficient, because the polymer and air are two distinct phases, which results in natural convection inside the domain as well. To consider natural convection within the current model, a temperature dependant density for the air would be necessary, whereas a Boussinesq approximation is also often employed provided that only small temperature gradients exist. However, for the former strategy, a compressible model needs to be considered for a variable density, whereas the current solver assumes the incompressibility condition. On the other side, the latter strategy cannot be considered since high temperature gradients arise in the SLS process. Therefore, the air was defined with a constant density and the only possible consideration for the natural convection was the temperature imposed on the top patch (see Section 5.5). One could argue that, if the entire section (that represents two layers) is fully melted, poor layer connectivity should not exist, after all, the particles are liquid and will remain in that state until the cooling phase. In the simulation, because the viscosity is set as a constant and has an a value between the virgin and used polymer, that would be the outcome. In reality, however, such simplification do not apply. Starting with the powder mixture, half of the material is guaranteed to be virgin, while the 92 5.7. Parameter Influence Studies 446.15 477 Melting point Temperature (K) Figure 5.23: Temperature distribution two seconds after the scan for the low energy case. remaining half is used, which is a very vague definition, since it does not specify how many hours or cycles each portion has. This entails that, among the used particles, different levels of viscosity are possible. Besides that, the surface tension, also set as a constant in the solver, decreases with material use [ 110 ], further contributing to a reduced coalescence. Additionally, and most likely the main reason for the poor interlayer coalescence in some regions, the viscosity levels decrease exponential with increasing temperature [ 109 ], which means that, although being liquid, particles that receive more energy will coalesce more and faster. In fact, by comparing only the experimental results with low and medium energy, it is safe to assume that both would contain formation defects and particles with high viscosity levels, still, in the medium energy case, the material, because it reaches much higher temperatures that greatly reduce its viscosity, is able to effective coalesce into an almost fully dense part. Taking that into account, when analysing the simulation results, the molten zone, after just three seconds, has a maximum temperature of 474K on the core of the first layer that, coincidently, corresponds to the temperature at which the viscosities for both the used and virgin material are known. Below and around that zone, the temperatures are considerably lower thus the viscosity will greatly increase. Despite the lack of information on the temperature dependent viscosity and its implementation in the model, the temperature profile on the simulation results, considering what was discussed, present the same behaviour as the experimental results. 93 5.7. Parameter Influence Studies 5.7 Parameter Influence Studies During the modelling phase, several simplifications had to be done due to some existing limitations, mainly in time, computational power and lack of appropriate data. Despite that, based on the comparison with the experimental results, the solver proved to be capable of reproducing the SLS process, at least for the temperature profile and temperature evolution. Knowing that the solver is not only behaving in a physically correct way, but also providing results close to the reality, allows the continuation to the next step, the study of parameters influence. With the comparison tests alone, it was visible what happens when more or less energy is given to the material. Performing more simulations with various energy density levels would not provide much more information on the solver capability or provide additional insights about the process, since the experimental results already provide the extreme energy levels for the employed material. The objective of the following studies is to change some parameters, one at a time, in the direction that would help improve the process by making it faster and/or cheaper, while securing the desired properties. During these analysis, other phenomena, such as dimensional accuracy and sintering evolution, can also be studied. For the base case, the used parameters will be those of the low energy case because, despite not being considered optimal, it was the only one where the entire scanned zone was not fully melted right after the scan. This will, hopefully, allow for a better analysis of the parameters influence, since the sintering depth can also be measured. This section will serve as an introduction to the proposed tests and the respective results and their discussion will be presented in Chapter 6. Fusion Depth From the experimental results, both extremes in terms of energy were assessed. That being said, in all the tested cases, the fusion depth was always superior to one layer as a consequence of the high energy needed to lower the material viscosity and avoid the defects seen in the low energy case. It would be interesting to find the lowest energy needed to effectively fuse one layer and reheat just a small portion of the previous one because, as mentioned in Section 2.4.3.4, not all PA12 behave in the same way PA2200 does. PA12 - Orgasol, for example, is much more resistant to ageing and, as a result, the difference in properties, such as viscosity, between virgin and used material, is much less accentuated. For a powder bed that contains material with similar properties, with reasonable viscosity levels, that does not require higher temperatures to lower them for increased coalescence, having a 94 6.2. Hatch Distance significant portion of unmelted material between the scans, which obviously is not desirable. 446.15 483 Melting point Temperature (K) (a) (b) Figure 6.6: Temperature distribution after 1 (a) and 2.5 (b) seconds for the highest hatch distance increase case. (a) (b) Figure 6.7: Phase distribution (melted as orange and solid as blue) after 1 (a) and 2.5 (b) seconds for the highest hatch distance increase case. As indicated earlier, since the current scan distance was excessive, one more test is going to be performed but, this time, with a hatch distance value corresponding to half of that considered previously. The results after the first scan are displayed in Figure 6.8. By decreasing the hatch distance, there is now an overlap of molten particles right after the scan. However, this only happens at the surface, where the temperature is higher, and might not be enough to melt the remaining of the layer. Similarly to the previous analysis, the case was run for a couple more seconds until the liquid fraction increase was negligible. Looking at the temperature profile in Figure 6.9, the maximum temperature is identical to the previous case. That was expected since there is no significant overlap that would greatly increase the temperature in one zone. Moreover, despite being closer, the hot spots are very similar for both cases. The liquid fraction, displayed in orange in Figure 6.10, shows that, with this lower hatch distance, a liquid bridge, with the thickness of almost one layer, is formed between the two scanned zones. By looking at the temperature profile and comparing with the base experimental 101 6.3. Coalescence case (Figure 5.22), the zone between scans, despite being liquid, has a lower temperature distribution, which will, most likely, cause insufficient coalescence. 446.15 512 Melting point Temperature (K) (a) (b) Figure 6.8: Temperature distribution after the first (a) and second (b) laser scans for the lowest hatch distance increase case. 446.15 483 Melting point Temperature (K) (a) (b) Figure 6.9: Temperature distribution after 1 (a) and 2.5 (b) seconds for the lowest hatch distance increase case. The hatch distance increase tests prove that, for the selected energy parameters, the minimum value for a good adhesion between two consecutive scans is the starting value used in the experimental tests. This implies that, despite not being mentioned in the article of reference, some tests should have been performed to establish this parameter and 0.3 mm was found to be an appropriate, which is in agreement with the simulation results. 102 6.3. Coalescence (b) (a) Figure 6.10: phase distribution (melted as orange and solid as blue) after 1 (a) and 2.5 (b) seconds for the lowest hatch distance increase case. 6.3 Coalescence Two identical cases, differing only on the viscosity of the liquid phase, were extended in time in order to visualize the coalescence development. Mixed Material The first test (Figure 6.11) is for the mixed material viscosity (2742.5 Pa · s). This test purpose is to be a base case for comparison with the next one, where the viscosity will be increased. Unfortunately, it is not possible to evaluate the accuracy of the results in terms of coalescence since no experimental data that allowed any type of comparison was found. Nonetheless, the observed behaviour was the expected, with the coalescence increasing with time. Used Material The results for the increase in viscosity are illustrated in Figure 6.12. As expected, by increasing the viscosity to that of the used material, the coalesce was significantly slower. Moreover, in agreement with the tests performed in Section 5.4.1, the particles, due to the increased viscosity, retain their shape for longer and seem to be more pulled towards each other instead of mainly diffusing. Because the surface tension is modelled as a constant and the values are the same for both cases, the results are almost identical, but slower and, as seen in Section 5.4.1, eventually, the final form will be the same. This is a limitation that could be surpassed if the surface tension values for the used and virgin material were known and if the model allowed for a temperature dependent value, which would require experimental tests to determine that dependency. Assuming that these results are a correct representation, considering the simplifications employed, besides the coalesce, one can also observe 103 6.4. Dimensional Accuracy (c) (d) (a) (b) Figure 6.11: Coalescence development after 3 (a), 5 (b), 7 (c) and 10 (d) seconds for the mixed material, with the respective phase distribution (melted as orange and solid as blue). layer densification due to the elimination of empty spaces between particles. Moreover, the partial melting of particles is also visible. The particle unmelted region will stick to the molten counterpart and that is why, in order to remove them, post processing operations, such as sand blasting [ 185 ], are often required. 6.4 Dimensional Accuracy In preceding sections, after just a few seconds, the melted fraction reached its peak. For the other experimental cases, since the energy provided to the material was higher, a longer time for stabilization is expected. The remaining cases will also be extended past their after scan time and an analysis to the melted region will be performed to assess the part dimensional accuracy. 104 6.4. Dimensional Accuracy (c) (d) (a) (b) Figure 6.12: Coalescence development after 3 (a), 5 (b), 7 (c) and 10 (d) seconds for the used material, with the respective phase distribution (melted as orange and solid as blue). Low Energy Case The first test case, with the low energy parameters, is shown in Figure 6.13. Beyond around two thirds of the laser radius, the material remains solid. Due to the low energy density, a significant portion of the scanned zone is left unmelted, especially along the thickness. Medium Energy Case The medium energy case results, displayed in Figure 6.14, surprisingly, despite the great increase in laser power, do not possess a much larger liquid fraction. In fact, the maximum lateral growth is only around three quarters of the laser radius and most of the increase in fused particles is visible in the layer thickness direction. This is, most likely, motivated by the typical polymer thermal insulating behaviour. 105 6.4. Dimensional Accuracy 1st scan 2nd scan Figure 6.13: Phase distribution (melted as orange and solid as blue) in relation to the scanned zone for the low energy case. 1st scan 2nd scan Figure 6.14: Phase distribution (melted as orange and solid as blue) in relation to the scanned zone for the medium energy case. High Energy Case Figure 6.15 shows the results for the high energy case. As anticipated, by increasing the laser power, more of the scanned area was able to phase change and the molten area is now closer to the scanned one. Despite the energy, for these process parameters, being so high that the material degrades, the phase fraction on the periphery of the scanned area does not increase significantly. This seems to be a consequence of the Gaussian intensity distribution that concentrates the energy at the 106 6.4. Dimensional Accuracy center of the scan in a way that, despite the laser power increase, not much of that energy reaches the periphery, and less notable changes are visible in this region, as evidenced by the temperature profiles from the tests in Section 5.6. Such behaviour indicates that the increase in width of the molten zone is mainly due to the accumulated energy at the high intensity zones that diffuses through heat conduction and not so much because the laser provides more energy at the periphery. However, the low polymer thermal conductivity, coupled with the low contact area between particles, hinders the energy diffusion. The contact area increases when the material coalesces, but that is a slow process and by the time a larger contact area is formed, most of the energy is lost through convection, as one can observe, for example, in Figure 6.9. All these reasons lead to a small lateral growth difference between the test cases, despite the distinct laser power values employed. The part growth is, most likely, a much larger concern for metal applications, since they have a much higher thermal conductivity and lower viscosity that allows for a greater contact area. 1st scan 2nd scan Figure 6.15: Phase distribution (melted as orange and solid as blue) in relation to the scanned zone for the high energy case. 107 CHAPTER 7Conclusions and Future Work SLS is one of the most popular AM techniques due to its capability of processing parts with complex geometries and good mechanical properties with almost no material wastage. However, the process is extremely complex due to its multi-physics nature and, as a consequence, it is yet not fully understood. There is still a significant amount of defects associated with the process and solving them experimentally has proved to be challenging, expensive and time-consuming. Nowadays, there is commercial software that simulates the SLS process, however, they present some major drawbacks, such as the high license cost, limiting the usage to only large companies, and the low flexibility to adapt the code to the user specific needs. For those who cannot afford such costs or limitations, open-source software is viable and an interesting alternative. Unfortunately, the literature on the numerical modelling of the process is not developed enough, especially for the use of open-source software. The present work focused on a global view of the SLS process and its simulations, with various aspects related to the process description, challenges and modelling discussed in more detail. Such analysis was necessary for the correct development of a computational tool to simulate the process at a particle length scale within an open-source framework. By using LIGGGHTS and OpenFOAM, a realistic powder bed was created and several studies were performed to assess the process evolution. Among the various studies, despite some necessary simplifications, the results were very satisfactory and in agreement with the available experimental data used for comparison. Starting with the thermal model, the laser source is behaving as expected and, based exclusively on the temperature profile after scan, the temperatures achieved by the material are a good representation of the real behaviour. The reflection model did not result in significant changes to the results, which was expected, according to the literature. Nonetheless, more efforts are needed to evaluate its adequacy and possible implementation in future simulations. Regarding the coalescence, arguably the most important phenomena on the process, several simplifications had to be employed. Starting with the viscosity of the material, it was 108 7.0. Dimensional Accuracy not possible to consider a portion of the particles with distinct viscosity levels to fully replicate the experimental practice, due to the incorrect interaction, predicted by the solver, between distinct liquid phases, thus, an average value had to be selected. Moreover, the viscosity was defined as constant, when in reality it varies significantly with temperature. By adapting the code and obtaining information on the temperature dependency of the viscosity for both the used and virgin material, a mixed powder bed could be used and the insufficient coalescence between particles, due to their high viscosity in low temperature zones, would be much more prominent and representative of the reality. Similarly, the surface tension is also treated as a constant in the current model, but if a temperature dependency was implemented, the results would likely also be more realistic. Sadly, no experimental data exists for that dependency as well. Lastly, as a closing statement, the developed strategy, using exclusively open-source software, manifests an immense potential. The powder bed is a realistic representation of reality and the particle size distribution influence is clearly visible in the test cases. The solver is also behaving as expected and the parameters influence is evident for all the performed studies. For a future work, if efforts are put into accurately determining the material properties that could not be correctly defined and the model is adapted to accept some dependencies and behave correctly between two distinct fluids, the current methodology could become a powerful tool to further assess, with higher detail, the process parameters influence. Moreover, for further validation of the solver capability, experimental cases with specific conditions, that would allow a more rigorous comparison with the simulation results, could also be performed. 109 References References [1] David Bourell, Joseph Beaman, Ming Leu, and David Rosen. A brief history of additive manufacturing and the 2009 roadmap for additive manufacturing: Looking back and looking ahead. Proceedings of RapidTech, pages 24--25, 2009. [2] Otto Munz. Photo-glyph recording, U.S. Patent 2 775 758, December 1956. [3] Wyn Swainson. Method, medium and apparatusfor producing three-dimensional figure product, U.S. Patent 4 041 476, July 1971. [4] Pierre Ciraud. Method and device for manufacturing any articles from any meltable material. Gernamy Patent DE2263777A1, December 1972. [5] Ross Housholder. Molding process, US Patent 4 247 508, January 1981. [6] Terry Wohlers and Tim Gornet. History of additive manufacturing. Wohlers report, (2016):38, 2016. [7] Gideon N. Levy, Ralf Schindel, and J.P. Kruth. Rapid manufacturing and rapid tooling with layer manufacturing (LM) technologies, state of the art and future perspectives. CIRP Annals, 52(2), 2003. doi: 10.1016/S0007-8506(07)60206-6. [8] Ian Gibson, David Rosen, and Brent Stucker. Powder Bed Fusion Processes. In Additive Manufacturing Technologies: 3D Printing, Rapid Prototyping, and Direct Digital Manufacturing, pages 107--145. Springer, New York, NY, 2015. ISBN 978-1-4939-2113-3. doi: 10.1007/978-1-4939-2113-3_5. [9] Samuel Clark Ligon, Robert Liska, Jürgen Stampfl, Matthias Gurr, and Rolf Mülhaupt. Polymers for 3D Printing and Customized Additive Manufacturing. Chemical Reviews, 117(15):10212--10290, 2017. ISSN 0009-2665. doi: 10.1021/acs.chemrev.7b00074. [10] Zhen Jiang, Broden Diggle, Ming Li Tan, Jekaterina Viktorova, Christopher W Bennett, and Luke A. Connal. Extrusion 3D Printing of Polymeric Materials with Advanced Properties. Advanced Science, 7(17):2001379, 2020. ISSN 2198-3844. doi: 10.1002/advs.202001379. [11] F42 Committee. Terminology for Additive Manufacturing Technologies. ASTM International. doi: 10.1520/F2792-12A. [12] Zhangwei Chen, Ziyong Li, Junjie Li, Chengbo Liu, Changshi Lao, Yuelong Fu, Changyong Liu, Yang Li, Pei Wang, and Yi He. 3D printing of ceramics: A review. Journal of the European Ceramic Society, 39(4):661--687, 2019. ISSN 0955-2219. doi: 10.1016/j.jeurceramsoc.2018.11.013. [13] Fuse 1+ 30W: Compact Selective Laser Sintering (SLS) 3D Printer, Available online: https://formlabs.com/3d-printers/fuse-1/ (Accessed on 28 October 2022). [14] 3D Printing Technology Comparison: FDM vs. SLA vs. SLS, Available online: https://formlabs.com/blog/fdm-vs-sla-vs-sls-how-to-choose-the-right-3d-printing-technology/ (Accessed on 28 October 2022). 110 References [90] Gexia Wang, Pingli Wang, Zhichao Zhen, Wei Zhang, and Junhui Ji. Preparation of PA12 microspheres with tunable morphology and size for use in SLS processing. Materials & Design, 87: 656--662, 2015. ISSN 0264-1275. doi: 10.1016/j.matdes.2015.08.083. [91] Alejandro Sosnik and Katia P. Seremeta. Advantages and challenges of the spray-drying technology for the production of pure drug particles and drug-loaded polymeric carriers. Advances in Colloid and Interface Science, 223:40--54, 2015. ISSN 0001-8686. doi: 10.1016/j.cis.2015.05.003. [92] Dietmar Drummer, Martha Medina-Hernández, Maximilian Drexler, and Katrin Wudy. Polymer Powder Production for Laser Melting Through Immiscible Blends. Procedia Engineering, 102: 1918--1925, 2015. ISSN 1877-7058. doi: 10.1016/j.proeng.2015.01.332. [93] R. D. Goodridge, K. W. Dalgarno, and D. J. Wood. Indirect selective laser sintering of an apatite-mullite glass-ceramic for potential use in bone replacement applications. Proceedings of the Institution of Mechanical Engineers, Part H: Journal of Engineering in Medicine, 220(1):57--68, 2006. ISSN 0954-4119. doi: 10.1243/095441105X69051. [94] Jochen Schmidt, Marius Sachs, Stephanie Fanselow, Meng Zhao, Stefan Romeis, Dietmar Drummer, Karl-Ernst Wirth, and Wolfgang Peukert. Optimized polybutylene terephthalate powders for selective laser beam melting. Chemical Engineering Science, 156:1--10, 2016. ISSN 0009-2509. doi: 10.1016/j.ces.2016.09.009. [95] G. Lumay, F. Boschini, K. Traina, S. Bontempi, J. C. Remy, R. Cloots, and N. Vandewalle. Measuring the flowing properties of powders and grains. Powder Technology, 224:19--27, 2012. ISSN 0032-5910. doi: 10.1016/j.powtec.2012.02.015. [96] Michael Van den Eynde, Leander Verbelen, and Peter Van Puyvelde. Assessing polymer powder flow for the application of laser sintering. Powder Technology, 286:151--155, 2015. ISSN 0032-5910. doi: 10.1016/j.powtec.2015.08.004. [97] Camden A. Chatham, Timothy E. Long, and Christopher B. Williams. A review of the process physics and material screening methods for polymer powder bed fusion additive manufacturing. Progress in Polymer Science, 93:68--95, 2019. ISSN 0079-6700. doi: 10.1016/j.progpolymsci.2019.03.003. [98] Silvia Vock, Burghardt Klöden, Alexander Kirchner, Thomas Weißgärber, and Bernd Kieback. Powders for powder bed fusion: A review. Progress in Additive Manufacturing, 4(4):383--397, 2019. ISSN 2363-9520. doi: 10.1007/s40964-019-00078-6. [99] Li Ma, Jeffrey T. Fong, Brandon Lane, Shawn P. Moylan, James J. Filliben, N. Alan Heckert, and Lyle E. Levine. Using Design of Experiments in Finite Element Modeling to Identify Critical Variables in Laser Powder Bed Fusion. NIST, 2016. [100] J. Frenkel. Viscous flow of crystalline bodies under the action of surface tension. J. phys., 9:385, 1945. [101] Ondej Pokluda, Céline T. Bellehumeur, and John Vlachopoulos. Modification of Frenkel’s model for sintering. AIChE Journal, 43(12):3253--3256, 1997. ISSN 1547-5905. doi: 10.1002/aic.690431213. [102] K. Wudy, D. Drummer, and M. Drexler. Characterization of polymer materials and powders for selective laser melting. AIP Conference Proceedings, 1593(1):702--707, 2014. ISSN 0094-243X. doi: 10.1063/1.4873875. 117 References [103] Jeou-shyong Wang and Roger S. Porter. On the viscosity-temperature behavior of polymer melts. Rheologica Acta, 34(5):496--503, 1995. ISSN 1435-1528. doi: 10.1007/BF00396562. [104] C. A. Hieber and H. H. Chiang. Shear-rate-dependence modeling of polymer melt viscosity. Polymer Engineering & Science, 32(14):931--938, 1992. ISSN 1548-2634. doi: 10.1002/pen.760321404. [105] Costas G. Gogos and Zehev Tadmor. Principles of polymer processing. John Wiley & Sons, 2013. [106] Ralph H. Colby, Lewis J. Fetters, and William W. Graessley. The melt viscosity-molecular weight relationship for linear polymers. Macromolecules, 20(9):2226--2237, 1987. ISSN 0024-9297. doi: 10.1021/ma00175a030. [107] Timothy P. Lodge. Reconciliation of the Molecular Weight Dependence of Diffusion and Viscosity in Entangled Polymers. Physical Review Letters, 83(16):3218--3221, 1999. doi: 10.1103/PhysRevLett.83.3218. [108] Manfred Schmid, Rob Kleijnen, Marc Vetterli, and Konrad Wegener. Influence of the Origin of Polyamide 12 Powder on the Laser Sintering Process and Laser Sintered Parts. Applied Sciences, 7(5):462, 2017. ISSN 2076-3417. doi: 10.3390/app7050462. [109] Leander Verbelen, Sasan Dadbakhsh, Michael Van den Eynde, Jean-Pierre Kruth, Bart Goderis, and Peter Van Puyvelde. Characterization of polyamide powders for determination of laser sintering processability. European Polymer Journal, 75:163--174, 2016. ISSN 0014-3057. doi: 10.1016/j.eurpolymj.2015.12.014. [110] Barry Haworth, Neil Hopkinson, David Hitt, and Xiaotao Zhong. Shear viscosity measurements on Polyamide-12 polymers for laser sintering. Rapid Prototyping Journal, 19(1):28--36, 2013. ISSN 1355-2546. doi: 10.1108/13552541311292709. [111] Yoshitomo Furushima, Christoph Schick, and Akihiko Toda. Crystallization, recrystallization, and melting of polymer crystals on heating and cooling examined with fast scanning calorimetry. POLYMER CRYSTALLIZATION, 1(2):e10005, 2018. ISSN 2573-7619. doi: 10.1002/pcr2.10005. [112] Y. Shi, Z. Li, H. Sun, S. Huang, and F. Zeng. Effect of the properties of the polymer materials on the quality of selective laser sintering parts. Proceedings of the Institution of Mechanical Engineers, Part L: Journal of Materials: Design and Applications, 218(3):247--252, 2004. ISSN 1464-4207. doi: 10.1177/146442070421800308. [113] Janaina L. Leite, Gean Salmoria, R. A. Paggi, Carlos Ahrens, and A. S. Pouzada. A study on morphological properties of laser sintered functionally graded blends of amorphous thermoplastics. 2010. ISSN 0268-1900. doi: 10.1504/IJMPT.2010.034272. [114] J. Choren, V. Gervasi, T. Herman, S. Kamara, and J. Mitchell. SLS Powder Life Study. pages 39--45, 2001. doi: 10.26153/tsw/3234. [115] Stefan Josupeit, Johannes Lohn, Eduard Hermann, Monika Gessler, Stephan Tenbrink, and Hans-Joachim Schmid. Material Properties of Laser Sintered Polyamide 12 as Function of Build Cycles Using Low Refresh Rates. pages 540--548. University of Texas at Austin, 2015. [116] T. J. Gornet, K. R. Davis, T. L. Starr, and K. M. Mulloy. Characterization of Selective Laser Sintering Materials to Determine Process Stability. pages 546--553, 2002. doi: 10.26153/tsw/4531. 118 References [117] Krassimir Dotchev and Wan Yusoff. Recycling of polyamide 12 based powders in the laser sintering process. Rapid Prototyping Journal, 15(3):192--203, 2009. ISSN 1355-2546. doi: 10.1108/13552540910960299. [118] Sasan Dadbakhsh, Leander Verbelen, Olivier Verkinderen, Dieter Strobbe, Peter Van Puyvelde, and Jean-Pierre Kruth. Effect of PA12 powder reuse on coalescence behaviour and microstructure of SLS parts. European Polymer Journal, 92:250--262, 2017. ISSN 0014-3057. doi: 10.1016/j.eurpolymj.2017.05.014. [119] D. T. Pham, K. D. Dotchev, and W. A. Y. Yusoff. Deterioration of polyamide powder properties in the laser sintering process. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, 222(11):2163--2176, 2008. ISSN 0954-4062. doi: 10.1243/09544062JMES839. [120] Bastian E. Rapp. Finite Element Method. In Microfluidics: Modelling, Mechanics and Mathematics, pages 655--678. Elsevier, 2017. ISBN 978-1-4557-3141-1. doi: 10.1016/B978-1-4557-3141-1. 50032-0. [121] Bastian E. Rapp. Finite Volume Method. In Microfluidics: Modelling, Mechanics and Mathematics, pages 633--654. Elsevier, 2017. ISBN 978-1-4557-3141-1. doi: 10.1016/B978-1-4557-3141-1. 50031-9. [122] Bastian E. Rapp. Finite Difference Method. In Microfluidics: Modelling, Mechanics and Mathematics, pages 623--631. Elsevier, 2017. ISBN 978-1-4557-3141-1. doi: 10.1016/ B978-1-4557-3141-1.50030-7. [123] Loong-Ee Loh, Chee-Kai Chua, Wai-Yee Yeong, Jie Song, Mahta Mapar, Swee-Leong Sing, ZhongHong Liu, and Dan-Qing Zhang. Numerical investigation and an effective modelling on the Selective Laser Melting (SLM) process with aluminium alloy 6061. International Journal of Heat and Mass Transfer, 80:288--300, 2015. ISSN 0017-9310. doi: 10.1016/j.ijheatmasstransfer.2014.09.014. [124] Liu Cao and Xuefeng Yuan. Study on the Numerical Simulation of the SLM Molten Pool Dynamic Behavior of a Nickel-Based Superalloy on the Workpiece Scale. Materials, 12(14):2272, 2019. ISSN 1996-1944. doi: 10.3390/ma12142272. [125] Liu Xin, Mhamed Boutaous, Shihe Xin, and Dennis A. Siginer. Numerical modeling of the heating phase of the selective laser sintering process. International Journal of Thermal Sciences, 120: 50--62, 2017. ISSN 1290-0729. doi: 10.1016/j.ijthermalsci.2017.05.017. [126] V. R. Voller and C. R. Swaminathan. Eral Source-Based Method for Solidification Phase Change. Numerical Heat Transfer, Part B: Fundamentals, 19(2):175--189, 1991. ISSN 1040-7790. doi: 10.1080/10407799108944962. [127] A. Laouadi, M. Lacroix, and N. Galanis. A numerical method for the treatment of discontinuous thermal conductivity in phase change problems. International Journal of Numerical Methods for Heat & Fluid Flow, 8(3):265--287, 1998. ISSN 0961-5539. doi: 10.1108/09615539810206348. [128] K. Nakamura, T. Watanabe, K. Katayama, and T. Amano. Some aspects of nonisothermal crystallization of polymers. I. Relationship between crystallization temperature, crystallinity, and cooling conditions. Journal of Applied Polymer Science, 16(5):1077--1091, 1972. ISSN 1097-4628. doi: 10.1002/app.1972.070160503. 119 References [129] Antonio Amado, Manfred Schmid, and Konrad Wegener. Simulation of warpage induced by non-isothermal crystallization of co-polypropylene during the SLS process. AIP Conference Proceedings, 1664(1):160002, 2015. ISSN 0094-243X. doi: 10.1063/1.4918509. [130] Rajen M. Patel. Crystallization kinetics modeling of high density and linear low density polyethylene resins. Journal of Applied Polymer Science, 124(2):1542--1552, 2012. ISSN 1097-4628. doi: 10.1002/app.35177. [131] Ahmed Hussein, Liang Hao, Chunze Yan, and Richard Everson. Finite element simulation of the temperature and stress fields in single layers built without-support in selective laser melting. Materials & Design (1980-2015), 52:638--647, 2013. ISSN 0261-3069. doi: 10.1016/j.matdes.2013.05.070. [132] Xiaoyong Tian, Gang Peng, Mengxue Yan, Shunwen He, and Ruijuan Yao. Process prediction of selective laser sintering based on heat transfer analysis for polyamide composite powders. International Journal of Heat and Mass Transfer, 120:379--386, 2018. ISSN 0017-9310. doi: 10.1016/j.ijheatmasstransfer.2017.12.045. [133] Jingjing Yang, Hanchen Yu, Huihui Yang, Fanzhi Li, Zemin Wang, and Xiaoyan Zeng. Prediction of microstructure in selective laser melted Ti6Al4V alloy by cellular automaton. Journal of Alloys and Compounds, 748:281--290, 2018. ISSN 0925-8388. doi: 10.1016/j.jallcom.2018.03.116. [134] Michael Gouge, Pan Michaleris, Erik Denlinger, and Jeff Irwin. Chapter 2 - The Finite Element Method for the Thermo-Mechanical Modeling of Additive Manufacturing Processes. In Michael Gouge and Pan Michaleris, editors, Thermo-Mechanical Modeling of Additive Manufacturing, pages 19--38. Butterworth-Heinemann, 2018. ISBN 978-0-12-811820-7. doi: 10.1016/B978-0-12-811820-7.00003-3. [135] Jyun-Rong Zhuang, Yee-Ting Lee, Wen-Hsin Hsieh, and An-Shik Yang. Determination of melt pool dimensions using DOE-FEM and RSM with process window during SLM of Ti6Al4V powder. Optics & Laser Technology, 103:59--76, 2018. ISSN 0030-3992. doi: 10.1016/j.optlastec.2018.01.013. [136] Chinnapat Panwisawas, Chunlei Qiu, Magnus J. Anderson, Yogesh Sovani, Richard P. Turner, Moataz M. Attallah, Jeffery W. Brooks, and Hector C. Basoalto. Mesoscale modelling of selective laser melting: Thermal fluid dynamics and microstructural evolution. Computational Materials Science, 126:479--490, 2017. ISSN 0927-0256. doi: 10.1016/j.commatsci.2016.10.011. [137] Chunlei Qiu, Chinnapat Panwisawas, Mark Ward, Hector C. Basoalto, Jeffery W. Brooks, and Moataz M. Attallah. On the role of melt flow into the surface structure and porosity development during selective laser melting. Acta Materialia, 96:72--79, 2015. ISSN 1359-6454. doi: 10.1016/j.actamat.2015.06.004. [138] Jennifer Lundkvist. CFD Simulation of Fluid Flow During Laser Metal Wire Deposition Using OpenFOAM : 3D Printing. PhD thesis, 2019. [139] P. A. Cundall and O. D. L. Strack. A discrete numerical model for granular assemblies. Géotechnique, 29(1):47--65, 1979. ISSN 0016-8505. doi: 10.1680/geot.1979.29.1.47. [140] H. Kruggel-Emden, S. Rickelt, S. Wirtz, and V. Scherer. A study on the validity of the multi-sphere Discrete Element Method. Powder Technology, 188(2):153--165, 2008. ISSN 0032-5910. doi: 10.1016/j.powtec.2008.04.037. 120 References [141] Zhaowei Xiang, Ming Yin, Zhenbo Deng, Xiaoqin Mei, and Guofu Yin. Simulation of Forming Process of Powder Bed for Additive Manufacturing. Journal of Manufacturing Science and Engineering, 138(8), 2016. ISSN 1087-1357. doi: 10.1115/1.4032970. [142] K. L. Johnson and I. Sridhar. Adhesion between a spherical indenter and an elastic solid with a compliant elastic coating. Journal of Physics D: Applied Physics, 34(5):683, 2001. ISSN 0022-3727. doi: 10.1088/0022-3727/34/5/304. [143] Alberto Di Renzo and Francesco Paolo Di Maio. Comparison of contact-force models for the simulation of collisions in DEM-based granular flow codes. Chemical Engineering Science, 59(3): 525--541, 2004. ISSN 0009-2509. doi: 10.1016/j.ces.2003.09.037. [144] Wentao Yan, Ya Qian, Wenjun Ge, Stephen Lin, Wing Kam Liu, Feng Lin, and Gregory J. Wagner. Meso-scale modeling of multiple-layer fabrication process in Selective Electron Beam Melting: Inter-layer/track voids formation. Materials & Design, 141:210--219, 2018. ISSN 0264-1275. doi: 10.1016/j.matdes.2017.12.031. [145] João Luís Oliveira Pedro. Numerical modelling of the filling stage of the injection moulding process. Master Dissertation, Integrated Master in Polymer Engineering, University of Minho, 2020. [146] Henk Kaarle Versteeg and Weeratunge Malalasekera. An Introduction to Computational Fluid Dynamics: The Finite Volume Method. Pearson Education Ltd, Harlow, England ; New York, 2nd ed edition, 2007. ISBN 978-0-13-127498-3. [147] F. Moukalled, L. Mangani, and M. Darwish. The Finite Volume Method in Computational Fluid Dynamics: An Advanced Introduction with OpenFOAM and Matlab, volume 113 of Fluid Mechanics and Its Applications. Springer International Publishing, Cham, 2016. ISBN 978-3-319-16873-9 978-3-319-16874-6. doi: 10.1007/978-3-319-16874-6. [148] OpenFOAM, Available online: https://www.openfoam.com/ (Accessed on 28 October 2022). [149] Docker: Accelerated, Containerized Application Development, Available online: https://www.docker.com/ (Accessed on 28 October 2022). [150] Windows Subsystem for Linux Documentation, Available online: https://learn.microsoft.com/enus/windows/wsl/ (Accessed on 28 October 2022). [151] ParaView, Available online: https://www.paraview.org/ (Accessed on 28 October 2022). [152] OpenFOAM v1806: New and updated solvers and physics, Available online: https://www.openfoam.com/news/main-news/openfoam-v1806/solver-and-physics (Accessed on 28 October 2022). [153] InterFoam - OpenFOAMWiki, Available online: https://openfoamwiki.net/index.php/InterFoam (Accessed on 28 October 2022). [154] OpenFOAM v6 User Guide: 3.5 Standard solvers, Available online: https://cfd.direct/openfoam/user-guide/v6-standard-solvers/ (Accessed on 28 October 2022). [155] W. H. Lee. A pressure iteration scheme for two-phase modeling. Technical Report LA-UR, 79-975, 1979. doi: 10.1016/j.ijheatmasstransfer.2015.02.037. 121 References [156] Guang Chen, Taotao Nie, and Xiaohong Yan. An explicit expression of the empirical factor in a widely used phase change model. International Journal of Heat and Mass Transfer, 150:119279, 2020. ISSN 0017-9310. doi: 10.1016/j.ijheatmasstransfer.2019.119279. [157] Dong-Liang Sun, Jin-Liang Xu, and Li Wang. Development of a vapor–liquid phase change model for volume-of-fluid method in FLUENT. International Communications in Heat and Mass Transfer, 39(8):1101--1106, 2012. ISSN 0735-1933. doi: 10.1016/j.icheatmasstransfer.2012.07.020. [158] Hyoungsoon Lee, Chirag R. Kharangate, Nikhin Mascarenhas, Ilchung Park, and Issam Mudawar. Experimental and computational investigation of vertical downflow condensation. International Journal of Heat and Mass Transfer, 85:865--879, 2015. ISSN 0017-9310. doi: 10.1016/j.ijheatmasstransfer.2015.02.037. [159] V. R. Voller and C. Prakash. A fixed grid numerical modelling methodology for convection-diffusion mushy region phase-change problems. International Journal of Heat and Mass Transfer, 30(8): 1709--1719, 1987. ISSN 0017-9310. doi: 10.1016/0017-9310(87)90317-6. [160] C. R. Swaminathan and V. R. Voller. A general enthalpy method for modeling solidification processes. Metallurgical Transactions B, 23(5):651--664, 1992. ISSN 1543-1916. doi: 10.1007/BF02649725. [161] P. C. Carman. Fluid flow through granular beds. Chemical Engineering Research and Design, 75:S32--S48, 1997. ISSN 0263-8762. doi: 10.1016/S0263-8762(97)80003-2. [162] M. F. Modest. Radiative Heat Transfer. Academic Press, New York, third edition edition, 2013. ISBN 978-0-12-386944-9. [163] gran model hertz model — LIGGGHTS v3.X documentation, Available online: https://www.cfdem.com/media/DEM/docu/gran_model_hertz.html (Accessed on 28 October 2022). [164] gran cohesion sjkr model — LIGGGHTS v3.X documentation, Available online: https://www.cfdem.com/media/DEM/docu/gran_cohesion_sjkr.html (Accessed on 28 October 2022). . [165] gran cohesion sjkr2 model — LIGGGHTS v3.X documentation, Available online: https://www.cfdem.com/media/DEM/docu/gran_cohesion_sjkr2.html (Accessed on 28 October 2022). . [166] Jun Ai, Jian-Fei Chen, J. Michael Rotter, and Jin Y. Ooi. Assessment of rolling resistance models in discrete element simulations. Powder Technology, 206(3):269--282, 2011. ISSN 0032-5910. doi: 10.1016/j.powtec.2010.09.030. [167] gran rolling_friction cdt model — LIGGGHTS v3.X documentation, Available online: https://www.cfdem.com/media/DEM/docu/gran_rolling_friction_cdt.html (Accessed on 28 October 2022). [168] Eric J. R. Parteli and Thorsten Pöschel. Particle-based simulation of powder application in additive manufacturing. Powder Technology, 288:96--102, 2016. ISSN 0032-5910. doi: 10.1016/j.powtec.2015.10.035. 122 References [169] Manfred Schmid and Konrad Wegener. Additive Manufacturing: Polymers Applicable for Laser Sintering (LS). Procedia Engineering, 149:457--464, 2016. ISSN 1877-7058. doi: 10.1016/j.proeng.2016.06.692. [170] PA2200 material data sheet 12-08 en, available online: https://www.shapeways.com/materials/pa11/attachment/material-data-sheet-nylon-12 (accessed on 28 october 2022). [171] 4.3 Mesh generation with the blockMesh utility, Available online: https://www.openfoam.com/documentation/user-guide/4-mesh-generation-and-conversion/4.3mesh-generation-with-the-blockmesh-utility (Accessed on 28 October 2022). [172] OpenFOAM: Manual Pages: setFields(1), Available online: https://www.openfoam.com/documentation/guides/latest/man/setFields.html (Accessed on 28 October 2022). [173] PA 12 - PA2200: Nylon for Industrial 3D Printing | EOS GmbH, available online: https://www.eos.info/en/additive-manufacturing/3d-printing-plastic/sls-polymermaterials/polyamide-pa-12-alumide (accessed on 28 october 2022). [174] M. C. Mielicki. Prediction of PA12 melt viscosity in Laser Sintering by a Time and Temperature dependent rheological model. 2012(9), 2012. [175] Gražyna Simha Martynková, Aleš Slíva, Gabriela Kratošová, Karla Čech Barabaszová, Soa Študentová, Jan Klusák, Silvie Brožová, Tomáš Dokoupil, and Sylva Holešová. Polyamide 12 Materials Study of Morpho-Structural Changes during Laser Sintering of 3D Printing. Polymers, 13 (5):810, 2021. ISSN 2073-4360. doi: 10.3390/polym13050810. [176] David R. Lide, Grace Baysinger, Swain Chemistry, Lev I. Berger, Robert N. Goldberg, and Henry V. Kehiaian. CRC Handbook of Chemistry and Physics. 2005. [177] P. Amend, C. Pscherer, T. Rechtenwald, T. Frick, and M. Schmidt. A fast and flexible method for manufacturing 3D molded interconnect devices by the use of a rapid prototyping technology. Physics Procedia, 5:561--572, 2010. ISSN 1875-3892. doi: 10.1016/j.phpro.2010.08.084. [178] Nitrogen - Density and Specific Weight vs. Temperature and Pressure, available online: https://www.engineeringtoolbox.com/nitrogen-n2-density-specific-weight-temperaturepressure-d_2039.html (accessed on 28 october 2022). . [179] Nitrogen - Thermal Conductivity vs. Temperature and Pressure, available online: https://www.engineeringtoolbox.com/nitrogen-n2-thermal-conductivity-temperature-pressured_2084.html (accessed on 28 october 2022). . [180] Nitrogen Gas - Specific Heat vs. Temperature, available online: https://www.engineeringtoolbox.com/nitrogen-d_977.html (accessed on 28 october 2022). . [181] Nitrogen - Dynamic and Kinematic Viscosity vs. Temperature and Pressure, available online: https://www.engineeringtoolbox.com/nitrogen-n2-dynamic-kinematic-viscosity-temperaturepressure-d_2067.html (accessed on 28 october 2022). . 123 References [182] Walter Schnell and Gaston Fischer. Carbon dioxide laser absorption coefficients of various air pollutants. Applied Optics, 14(9):2058--2059, 1975. ISSN 2155-3165. doi: 10.1364/AO.14. 002058. [183] G. M. Vasquez, C. E. Majewski, B. Haworth, and N. Hopkinson. A targeted material selection process for polymers in laser sintering. Additive Manufacturing, 1--4:127--138, 2014. ISSN 2214-8604. doi: 10.1016/j.addma.2014.09.003. [184] Tobias Laumer, Thomas Stichel, Marius Sachs, Philipp Amend, and Michael Schmidt. Qualification and modification of new polymer powders for laser beam melting using Ulbricht spheres. pages 255--260. 2013. ISBN 978-1-138-00137-4. doi: 10.1201/b15961-48. [185] A. C. Lopes, A. M. Sampaio, and A. J. Pontes. The influence of the energy density on dimensional, geometric, mechanical and morphological properties of SLS parts produced with single and multiple exposure types. Progress in Additive Manufacturing, 7(4):683--698, 2022. ISSN 2363-9520. doi: 10.1007/s40964-021-00254-7. [186] Christian Nelson, Kevin McAlea, and Damien Gray. Improvements in SLS Part Accuracy. pages 159--169, 1995. doi: 10.15781/T2222RR22. 124