A Layer-by-Layer Approach for Simulating Residual Stresses in AM

11th European LS-DYNA Conference 2017, Salzburg, Austria A Layer-by-Layer Approach for Simulating Residual Stresses in AM Mikael Schill1 , Niclas Strömberg2 1 Dynamore Nordic, [email protected] 2 Örebro University, [email protected] 1 Abstract Distortions and residual stresses developed during additive manufacturing (AM) might be so severe that the design does not meet tolerance requirements or even cracks and fail. Predictions of such distortions and residual stresses by finite element analysis are therefore preferable performed early during the design development in order to reduce such problems. One approach for this kind of analysis is to apply techniques from welding simulations letting the Goldak heat source move over each layer of the material. However, this is computationally most costly and one might also argue how well this actually represent the AM process. We suggest to simplify this boundary condition to instead apply the heat layer-by-layer. We mean that this will represent the AM process as well as an approach of using the Goldak heat source and it is much more computational efficient than applying Goldaks heat source. The model is simply obtained by slicing the geometry in several layers using an in-house script. Then, a volume heat source is applied layer-by-layer, where each layer is activated starting from the build plate until the final build layer. Sequentially thermomechanical analysis using LS-Dyna is applied, where the heat capacity, conductivity, Youngs modulus, Poissons ratio, expansion coefficient, yield stress, and the hardening modulus are given as functions of temperature. The classical heat equation and J2-plasticity with kinematic or isotropic hardening are adopted by using *MAT THERMAL CWM and *MAT CWM, which makes the activation of the layers straight-forward using the ghost feature implemented in these material models. After all layers are activated and heated, a restart of cooling analysis is performed and, finally, the part is cut from the building plate using a third restart analysis. The approach is demonstrated for a benchmark of an open cylinder printed in Inconel 718 showing the development of distortions and residual stresses in all stages of the AM process from heating during building until the final cut from the build plate. We suggest that this benchmark can serve as an experimental setup in order to validate the suggested approach. This is discussed in the paper and is a topic for future work. c 2017 Copyright by DYNAmore GmbH 11th European LS-DYNA Conference 2017, Salzburg, Austria 2 Introduction Components subjected to changes in temperature will expand or contract. The corresponding thermal strains will result in thermal stresses if the deformations are prevented by external boundary conditions or internally by different regions of the component. These thermal stresses might cause severe malfunction. A well-known example shown in every basic course on solid mechanics is railway buckling in hot sunny days [1]. Frictional heating of disc brakes when retarding vehicles is another example where corresponding thermal stresses might lead to severe malfunction [2]. In manufacturing processes, such as e.g. casting, thermal stresses develop during solidification which eventually result in distortions or residual stresses which might influence the lifetime of the component [3]. A standard approach to simulate these kinds of residual stresses is to set up a J2-plasticity model with temperature dependent properties, i.e. letting Young’s modulus, Poisson’s ratio, thermal expansion coefficient, yield stress and hardening modulus depend on the temperature [4]. In addition, one also need the temperature history of the process, which can be found by solving the classical heat equation with temperature dependent specific heat and thermal conduction for proper boundary conditions. In Gustafsson et al. [5], these two kind of simulations were done sequentially in order to find residual stresses from solidification of grey iron castings. The residual stresses in a stress lattice were simulated and then validated with experiments. The results from simulations correlated well with measured values. The same J2-plasticity model was used in Rashid and Strömberg [6] to simulate residual stresses in disc brakes as a result from severe braking downhill of a heavy truck. Of course, the temperature history was now obtained in a different manner. Instead of cooling the component from an initial state of melt temperature, a new thermomechanical contact approach was utilized where the disc is formulated in an Eulerian framework treating the disc as a fluid [7]. Recently, this framework was used to se up a thermo-flexible multi-body model of a brake dynamometer [8]. The temperature history from these simulations was then exported to LS-Dyna to find the corresponding residual stresses. The J2-plasticity model was set up using *MAT CWM, which originally was implemented to simulate welding [9]. Recent developments of the models for welding in LS-Dyna were recently presented by Schill et al. [10]. These models are equipped with a ghost functionality which makes it possible to activate the material when it reaches a certain temperature interval. In this work, we use *MAT CWM to simulate residual stresses and thermal distortions developed in the AM process. Instead of using Goldak’s heat source [11], we apply the heat layer-by-layer in the build direction as an internal volume heat production by using *LOAD HEAT GENERATION. The model is set up automatically with an in-house script. The component to be build and the build c 2017 Copyright by DYNAmore GmbH 11th European LS-DYNA Conference 2017, Salzburg, Austria plate is imported as one part, which automatically is sliced into several layers and then exported to three LS-dyna jobs. In the first job, the component is built by applying the internal heat source layer-by-layer. Then, the component and build plate are cooled in the chamber and outside in room temperature and, finally, the component is cut from the build plate. The suggested layer-by-layer approach seems promising and this is demonstrated at the end of the paper by studying on open cylinder and a Michell structure obtained from topology optimization. 35 600 30 k [W/m] c [J/kgK] 550 500 25 20 450 15 400 10 0 200 400 600 800 1000 1200 1400 1600 0 200 400 600 T (a) 800 1000 1200 1400 1600 T (b) Specific heat Conductivity Figure 1: Specific heat and thermal conductivity as function of temperature. 3 Governing equations In this section a brief overview of the governing equations for the thermomechanical analysis performed in this work using LS-Dyna is presented. The analysis is sequential coupled using *MAT CWM THERMAL and *MAT CWM. The balance in linear momentum and energy are governed by divσ = 0, ρcṪ = k∆T + σ : ǫ̇p + r, (1) where σ is a stress, ρ is the mass density, c = c(T ) is the heat capacity, k = k(T ) is the thermal conductivity, σ : ǫ̇p represents dissipation due to plastic work and r is internal heat production. The latter term r is discussed in the next section. The stress is governed by temperature dependent J2-plasticity with isotropic and/or kinematic hardening. More details are presented below. The heat capacity and the conductivity are given by where c(T ) = cm (T )γ + cg (1 − γ), k(T ) = km (T )γ + kg (1 − γ), (2)  T max − T start γ = min 1, max(0, end ) , T − T start (3)  c 2017 Copyright by DYNAmore GmbH 11th European LS-DYNA Conference 2017, Salzburg, Austria cm = cm (T ) is heat capacity as function of temperature T for active material, km = km (T ) is thermal conductivity for active material, cg is a constant ghost parameter used as artificial specific heat capacity before cm is fully activated and kg is a ghost parameter used as artificial conductivity before km is activated. The activation is adjusted by γ in (3), where T start and T end define a temperature interval for activation. Furthermore, T max is the highest temperature obtained during the simulation. Thus, if T max reaches T end , then our thermal properties are completely activated meaning that c(T ) = cm (T ) and k(T ) = km (T ). The temperature dependent specific heat and conductivity used in this work are presented in Figure 1. These values are taken from Deshpande et al. [12] and represent properties for Inconel 718. Other material properties presented in this section are also taken from that work. 2 ×10 5 2.2 ×10 -5 2 1.8 1.6 α [1/K] E [MPa] 1.8 1.6 1.4 1.4 1.2 1 secant tangent 1.2 0.8 1 0.6 0 200 400 600 800 1000 1200 1400 1600 0 200 400 600 T (a) 800 1000 1200 1400 1600 T Young’s modulus (b) Thermal expansion Figure 2: Young’s modulus and thermal expansion as function of temperature. Stress σ and strain ǫ are related by temperature dependent J2-plasticity with isotropic and/or kinematic hardening. The strain is decomposed into one elastic part ǫe , a plastic strain ǫp and a thermal strain ǫt , i.e. ǫ = ǫe + ǫp + ǫt . (4) Isotropic elasticity is set up for the stress and the elastic strain, where Young’s modulus E = E(T ) and Poisson’s ratio ν = ν(T ) are governed by E(T ) = Em (T )γ + Eg (1 − γ), ν(T ) = νm (T )γ + νg (1 − γ), (5) where γ follows (3), Em = Em (T ) and νm = νm (T ) are temperature dependent elastic properties of our material, and Eg and νg are artificial elastic properties for non-activated ghost material. Consequently, E = Em and ν = νm when T max reaches T end . Young’s modulus Em (T ) as function of time used in this work is presented in Figure 2. c 2017 Copyright by DYNAmore GmbH 11th European LS-DYNA Conference 2017, Salzburg, Austria 300 2500 2000 H [MPa] σy [MPa] 250 200 150 100 1500 1000 500 50 0 0 0 200 400 600 800 1000 1200 1400 1600 0 200 400 600 T (a) 800 1000 1200 1400 1600 T Yield stress (b) Hardening modulus Figure 3: Yield stress and hardening modulus as function of temperature. The plasticity is governed by temperature dependent J2-plasticity with isotropic and/or kinematic hardening. This is formulated by introducing the deviatoric stress 1 s = σ − tr(σ)I (6) 3 and the back-stress η which is used to translate the yield surface. The yield surface is represented by r 3 f= (s − η) : (s − η) − σy (T ) − βH(T )ǫp , (7) 2 where σy = σy (T ) is the yield stress as function of temperature and H = H(T ) is the hardening modulus as function of temperature. The yield stress and hardening modulus used in this work is presented in Figure 3. Furthermore, 0 ≤ β ≤ 1 apperaing in (7) is a parameter defining the level of isotropic and kinematic hardening. Thus, for β = 1 pure isotropic hardening is obtained and β = 0 corresponds to pure kinematic hardening, see also (11). Otherwise, a mixture of isotropic and kinematic hardening is obtained. The effective plastic strain ǫp in (7) is defined by r 2 p p p ǫ̇ = ǫ̇ : ǫ̇ . (8) 3 In addition, the plastic strain evolves according to ǫ̇p = λ̇ 3 s−η ∂f , = λ̇ r ∂σ 2 3 (s − η) : (s − η) 2 (9) where ǫ̇p = λ̇ ≥ 0, f ≤ 0, λ̇f = 0. c 2017 Copyright by DYNAmore GmbH (10) 11th European LS-DYNA Conference 2017, Salzburg, Austria The back-stress is updated using η̇ = λ̇(1 − β)H(T ) r s−η 3 (s − η) : (s − η) 2 . (11) The plastic flow rules in (9) and (11) are equipped with annealing that is activated in a temperature interval Tas ≤ T ≤ Tae such that   T − Tae p p ǫn+1 = ǫn max 0, min(1, s ) , Ta − Tae   (12) T − Tae p p η n+1 = η n max 0, min(1, s ) Ta − Tae are performed in the stress update after convergence in the radial return step. Figure 4: The in-house script setting up the model automatically. Finally, the thermal strains are governed by ǫ̇t = α(T )Ṫ I, (13) α = α(T ) = αm (T )γ + αg (1 − γ) (14) where is the tangent thermal expansion coefficient using αm = αm (T ) for the active material and αg for the ghost material following the same procedure as presented previously. It is most important to understand that the tangent thermal expansion coefficient is used and not the secant coefficient α̂ = α̂(T ) which appears more frequently. For instance, the secant expansion coefficient is used c 2017 Copyright by DYNAmore GmbH 11th European LS-DYNA Conference 2017, Salzburg, Austria in [12]. By using the secant coefficient, the thermal strains can be expressed as ǫt = α̂(T )(T − T0 )I, (15) where T0 is a reference temperature. By taking the time derivative of (15), one obtains dα̂(T ) Ṫ (T − T0 )I + α̂(T )Ṫ I. (16) ǫ̇t = dT From this expression one can see that the tangent thermal expansion coefficient dǫt α(T ) = (17) dT is related to the secant coefficient by α(T ) = dα̂(T ) (T − T0 )I + α̂(T )I. dT (18) Figure 5: A sequential layer-by-layer by approach for simulating residual stresses in AM. 4 The layer-by-layer approach The acronym CWM stands for Computational Welding Mechanics. Thus, the materials *MAT CWM THERMAL and *MAT CWM were originally developed for simulating thermal stresses and distortions developed during welding followed by cooling. In addition, Goldak’s moving weld heat source [11] is implemented in *BOUNDARY THERMAL CWM. In this work, we use *MAT CWM THERMAL and *MAT CWM in order to simulate residual stresses and distortions develop during additive manufacturing followed by cooling. Instead of applying the heat power using Goldak’s approach, we simply suggest to slice the part into several layers c 2017 Copyright by DYNAmore GmbH 11th European LS-DYNA Conference 2017, Salzburg, Austria and then apply the heat power as an internal heat source for each layer using *LOAD HEAT GENERATION, see r in (1). This is done sequentially in the build direction starting from the build plate towards the top of the component to be build. In order to set up this model, an automatic preprocessor is developed using Matlab with a simple GUI. By setting different parameters, such as e.g. number of layers, the GUI imports the component to be build together with the build plate as one part, then a script automatically slice the component into several layers. Each layer gets unique material properties and internal heat sources are provided for each layer which are activated sequentially in the build direction. The model is exported to three jobs: one for building the component (BUILD), a second one for cooling the component and build plate (COOL), and, finally, a job for cutting the component from the build plate (CUT). In Figure 4, this pre-processing is illustrated for an open cylinder which we think could serve as a good geometry for validating the simulations with experiments. In Figure 5, we show results from a simulation of the open cylinder using the proposed layer-by-layer approach. The sequential analysis starts by building the cylinder layer-by-layer. The building of the cylinder is almost completed in the plot to the left in Figure 5, where the color in red represents melted material and blue color above the red is ghost material. Then cooling of the component and the build plate follows in sequence by a dynain restart analysis and finally the cylinder is cut from the build plate. The pink geometry presented here in the right plot of Figure 5 represents the reference geometry and the green geometry is the actual geometry obtained due to the thermal distortions from the AM process. Finally, in Figure 6, we illustrate how our suggested approach can be used in a CAE-driven process using topology optimization. A block fixed at the corners and a force acting at the center is optimized using topology optimization [13] and the optimal geometry is then imported into the in-house script that automatically slice the geometry and set up the AM simulations by exporting the following three jobs: BUILD, COOL and CUT. The last picture in the loop is showing the heat profile when only a few layers remain before the component is completely build. 5 Concluding remarks In this work we suggest a layer-by-layer approach for simulating residual stresses and thermal distortions developed in additive manufacturing. A model of a build plate with the component to be build is set up automatically using an in-house script. The approach is promising and the next step is to validate the simulations with experiments. The open cylinder studied in this work could server as a proper benchmark for such validations. c 2017 Copyright by DYNAmore GmbH 11th European LS-DYNA Conference 2017, Salzburg, Austria Figure 6: CAE-driven process of a Michell structure using topology optimization and our proposed layer-by-layer approach for AM simulation. Acknowledgement - We are grateful for the financial support from Vinnova (www.vinnova.se) in the project Digi3D (www.digi3d.org). 6 Literature [1] G. Yang & M.A. Bradford, Thermal-induced Buckling and Postbuckling Analysis of Continuous Railway Tracks, International Journal of Solids and Structures, 97-98, 637-649, 2016. [2] P. Dufrénoy & D. Weichert, A Thermomechanical Model for the Analysis of Disc Brake Fracture Mechanisms, Journal of Thermal Stresses, 26, 815– 828, 2003. [3] R.W. Lewis & K. Ravindren, Finite Element Simulation of Metal Casting, International Journal for Numerical Methods in Engineering, 47, 29-59, 2000. [4] L.E. Lindgren, Finite Element Modelling and Simulation of Welding, Part 1 - Part 3, Journal of Thermal Stresses, 24, 141–334, 2001. [5] E. Gustafsson, M. Hofwing & N. Strömberg, Residual Stresses in a Stress Lattice - Experiments and Finite Element Simulations Journal of Materials Processing Technology, 209, 4320–4328, 2009. c 2017 Copyright by DYNAmore GmbH 11th European LS-DYNA Conference 2017, Salzburg, Austria [6] A. Rashid & N. Strömberg, Sequential Simulation of Thermal Stresses in Disc Brakes for Repeated Braking, Journal of Engineering Tribology, 227, 919–929, 2013. [7] N. Strömberg, An Eulerian Approach for Simulating Frictional Heating in Disc-Pad Systems, European Journal of Mechanics, A/Solids, 30, 673–683, 2011. [8] N. Strömberg, A Virtual Test Rig for Disc Brake Systems by Adopting a Thermo-Flexible Multi-Body Approach, in the proceedings of the 34th annual SAE 2016 Brake Colloquium & Exhibition, Arizona, USA, September 25–28, 2016. [9] P. Lindström, Improved CWM Platform for Modelling Welding Procedures and Their Effects on Structural Behaviour, PhD Thesis, Production Technology, No. 6, University West, 2015. [10] M. Schill, A. Jernberg & T. Klöppel, Recent Developments for Welding Simulations in LS-DYNA and LS-PrePost, in the proceedings of the 14th International LS-Dyna Users Conference, Michigan, USA, June 12–14, 2016. [11] J. Goldak, A. Chakravarti & M. Bibby, A New Finite Element Model for Welding Heat Sources, Metallurgical Transactions B, 15, 299-305,1984. [12] A.A. Deshpande, D.W.J. Tanner, W. Sun, T.H. Hyde & G. McCartney, Combined Butt Joint Welding and Post Weld Heat Treatment Simulation using SYSWELD and ABAQUS, in the proceedings of the Institution of Mechanical Engineers Part L Journal of Materials Design and Applications, 225, 1–10, 2011. [13] N. Strömberg, Topology Optimization of Structures with Manufacturing and Unilateral Contact Constraints by Minimizing an Adjustable ComplianceVolume Product, Structural and Multidisciplinary Optimization, 42, 341-350, 2010. c 2017 Copyright by DYNAmore GmbH