The Cooling Rate and Residual Stresses in an AISI 310 Laser Weld: A Meso-Scale Approach

Abstract

A three-dimensional coupled temperature-displacement finite element model was developed to generate values of temperature distribution, cooling rate, and residual stresses at the meso-scale level in a thick sheet AISI 310 laser welding test sample. High cooling rates (cooling time from liquid-to-solid temperatures) ranging from 960 °C/s to 2400 °C/s were observed when the calculations were made at the meso-scale level. These high cooling rates that arise during the formation of the weld pool originate the highest observed residual stresses that evolve throughout the weld during the entire heating and cooling cycles. An ABAQUS CAE meso model with dimensions of 10 × 5 × 1 mm (element size 100 µ) constructed from a global macro model of 40 × 10 × 75 mm (element size 1 mm) via the submodeling technique is presented in the present paper. In both analyses, macro and meso, the C3D8T thermally coupled brick, trilinear displacement and temperature elements were used. To mesh the entire plate with elements of regular size 100 × 100 × 100 µ, a total of 30 million elements are necessary. With the present approach, 1 macro mesh of 30 thousand elements (1 × 1 × 1 mm) and a meso mesh of 50 thousand elements (100 × 100 × 100 µ) were enough to simulate the weld problem at the meso-scale level.

1. Introduction

Epitaxial solidification can be achieved if the welding process is adequately performed. It requires the solidification theory to control the columnar-to-equiaxed transition caused for the macroscopic heat input on the microstructure development. The analytical heat-flux solution methods presented in early approaches [1,2,3,4,5,6,7,8] use semi-empirically based computer programs with constant material properties. Consequently, several mathematical simplifications of the integral solution of the non-linear transient governing heat equation were made. As demonstrated in refs. [2,3,4,5], the heat input has a profound influence on the dendrite growth velocity, growth pattern, and melt pool geometry. The high peak temperatures and cooling rates (localized thermal cycling), grain structure, distortion, and reduced strength of a weld joint are produced by the intense heat input from fusion welding [9,10]. In the present work, the temperature-dependent material properties and forced convection due to Marangoni effects were incorporated into a three-dimensional coupled temperature-displacement finite element model.

The microstructure development and the evolution of the residual stresses are important aspects, both in the repair of gas turbine components and in relation to safety in the nuclear industry. The parameters that describe the heat input from the heat source are the most critical input data in the welding thermal analysis [9,10,11]. In regard to the residual stresses and distortions, the experimental optimization of the welding technique requires measurements and prototyping, which are expensive and time consuming. However, with the virtual prediction of residual stress evolution and microstructure development in the melt pool, the process can be optimized in the early stages of prototyping. Although these phenomena are difficult to simulate due to a lack of exact knowledge of welding conditions and material behavior, their prediction can provide substantial assistance in the accurate design and fabrication of welded structures. Accurately calculating the transient temperature field, that causes non-equilibrium phase formation in and around the welded joint, is the critical first step in creating a science base, not only for the design and analysis of welds, but also for the repair applications. The residual stresses that develop both in the heat-affected zone (HAZ) and the fusion zone (FZ) are harmful to the integrity and service-life of the welded part. In ref. [12], it is shown that residual stresses can cause problems, such as hydrogen-induced cracking, stress corrosion cracking, distortion, as well as initiating fracture and degrading the corrosion resistance of welded structures. To reduce them, pre- and post-heating is often used.

The accurate calculation and measurement of residual stresses at different scale levels still remains a major issue. The numerical prediction of thermal cycles is the prerequisite for subsequent calculations of residual stress distribution and microstructure development in welds [12,13,14,15]. Several key problems and issues remain to be addressed, especially at small scales. The main difficulties in the quantitative analysis of fusion welding are the scarcity of relevant data of properties at very high temperatures, i.e., much higher than the melting point, and the complexity of the physical processes. That is, the effects of fluid flow in the weld pool produced by the complex coupled transport phenomena are difficult to be addressed. Moreover, because of the extreme cooling conditions encountered during laser welding, the events that follow welding are far from equilibrium and, as a result, non-equilibrium phase formation in and around the welded joint occur. Therefore, the simulation of phase transformations (and the resulting microstructures), distortion and residual stress evolution in weldments remains a great challenge.

Physically based models used to simulate the entire thermal cycle of the welding process can provide a new understanding of the distortion evolution and microstructure development [14]. However, the physically based modeling of the entire thick plate and the melt pool as a whole, is a computationally formidable task. The lack of accurate stress–strain constitutive relations and the use of coarse grids create differences in the calculated and measured residual stresses. The uncertainties in the calculation procedures include the inaccuracies in the calculation of the thermal cycles as well as the approximations in the general frame for the mechanical model, especially when important solid-state transformations take place [12,13,14,15]. In previous works [10,16], it was well established that the cooling time from liquid-to-solid temperatures is longer at the fusion line and shorter at the weld centerline. As such, the cooling rate through the solidification temperature range increases and the dendrite arm spacing decreases from the fusion line to the centerline. In agreement with these previous numerical results, in ref. [16], it is also shown through weld micrographs that the solidification microstructure becomes finer from the fusion line to the centerline.

The modeling of forming processes involving the melting and solidification of metals, such as additive manufacturing (AM), casting or welding, all share some common features associated with those thermomechanical aspects that play an important role in the final material properties at the macroscale. The multiscale and Multiphysics phenomena involved, at different levels of significance, are imperative to be addressed when it comes to characterizing important aspects that perform a role in the final properties of the resulting parts. When running simulation studies with large assemblies, the submodeling feature allows us to refine the results for critical components without having to re-run the analysis for the entire assembly. In technical terms, submodeling helps transfer complex global loads from the entire structure to local sub-regions to obtain much more accurate results for the unknown field variables, such as temperature, stresses, and strains. Submodeling allows us to start the analysis process on a large assembly by running a “relatively coarse” parent analysis. This analysis may use coarser mesh than is really needed, might make simplified assumptions about contact and boundary conditions, or might neglect connection details.

The present research aims not only to verify that, at the microscale level, materials display strong size effects when the characteristic length-scale associated with non-uniform plastic deformation is on the order of microns, but also to clarify the link between cooling rates and solidification grain substructures. The prediction of crack initiation and propagation in weld samples is very important to identify direct energy deposition (DED) additive manufacturing (AM) process parameters for preserving the single-crystal nature of nickel-base superalloy gas turbine blades during repair procedures. The accurate prediction of the melt-pool shape and solidification parameters as a function of the processing parameters (including filler metal additions) is also important in the development of microstructural models for the design of functionally graded materials, using either laser additive manufacturing, or Wire + Arc additive manufacturing (WAAM).

2. The Model

The temperature T(x; y; z; t) at any location (x; y; z) and time (t) with respect to the moving heat source is calculated by solving the 3D transient nonlinear heat equation:

$$\frac{\partial }{\partial x}\left(k\frac{\partial T}{\partial x}\right)+\frac{\partial }{\partial y}\left(k\frac{\partial T}{\partial y}\right)+\frac{\partial }{\partial z}\left(k\frac{\partial T}{\partial z}\right)+\dot{Q}=\rho c_p\frac{\partial T}{\partial t}$$

(1)

Here, k, ρ and c_p are the temperature-dependent thermal conductivity, density and specific heat, respectively. T is temperature, t is time, and $\dot{Q }$ the internal heat-source term. The solution depends on the physical conditions existing at the boundaries of the workpiece and the conditions existing in the workpiece at some initial time. In this work, $\dot{Q }$ is 0, the initial temperature T_o is 20 °C, and the latent heat was ignored. Convective and radiative boundary conditions for the heat exchange between the top surface of the workpiece and the surroundings beyond the laser heat source are expressed in the following equation:

$$-k\frac{\partial T}{\partial y}{|}_{top}+q\left(x,z,t\right)=h_t\left(T-T_s\right)+\sigma \epsilon \left(T^4-T_s^4\right)$$

(2)

Here, h_t is the convection heat-transfer coefficient, T_s is the surrounding temperature, σ is the Stefan–Boltzman constant and ε is the surface emissivity. Values of ε = 0.7 and h_t= 242 Wm⁻² K⁻¹ to include forced convection in the area directly beneath the nozzle of the laser heat source were used. Convection heat loss at the workpiece top surface and all the other weldment surfaces were also accounted for using a value of h_t = 10 Wm⁻²K⁻¹. The heat input q(x,z,t) from the heat source to the workpiece is represented by a moving Gaussian power density distribution (Figure 1a and Equation (3)) applied over the top surface of the specimen during a period of time that depends on the welding speed (v). The characteristic heat distribution parameter, C, was selected based on previous experimental measurements of the fusion zone and heat-affected zone widths [17]. The mentioned parameter is very important because the shape and distribution of the heat input from the heat source greatly depend on its value. A distribution parameter C = 4.5 mm was used in all the simulations. The Gaussian distribution model is usually used in simulations of welding processes with high power densities, i.e., arc, laser, or electron beam welding. With the help of this heat source model, in which energy is distributed according to a Gaussian profile, it is possible to map deep penetration, while maintaining a small width, characteristic of “keyhole” welding techniques [18].

The Gaussian power density distribution is described as follows:

$$q\left(x,z,t\right)=\frac{3Q}{\pi C^2}ex{p}^{\left\{-3\left[{\left(z-\mathrm{v}t\right)}^2+x^2\right]/C^2\right\}}$$

(3)

Q = η_thVi, where η_th is the thermal efficiency, V is the voltage and i is the electric current. The thermal efficiency n_th is defined as the ratio of the power required for melting the volume of metal in the fusion zone P_th and the emitted laser beam power P_E [19], and is described by:

$$n_{th}=\frac{P_{th}}{P_E}=\frac{A .\rho . \mathrm{v} . \left(c_{p }. \left(T_M-T_A\right)+h_M\right)}{P_E}$$

(4)

where A is the area of the cross-section of the molten pool, $\rho$ the density, v the welding speed, c_p the specific heat capacity, T_M the melting temperature, T_A the ambient temperature, and h_M the latent heat for melting. Using Equation (4), it should be considered that the energy for heating the material outside of the molten pool as well as the energy for overheating the molten pool are assumed to be energy losses. Equation (4) can also be written as:

$$Heat input=\frac{n_{th }.P_E }{\mathrm{v}} \left[\frac{J}{mm}\right]$$

(5)

A thermal efficiency of η_th = 0.5 and a heat input of 120 J/mm were used in all the simulations. The ABAQUS [20] user subroutines FILM and DFLUX were written to account for the convection, radiation and heat-input distribution. Of a similar manner, the dependent elastic modulus and yield strength (0.2% offset) for AISI 310 austenitic steel obtained from ref. [21] and plotted in Figure 2d,e, respectively, were included into the ABAQUS material module. It is important to note that the AISI 310 steel can be welded relatively easily to other weldable steels, and also to the HAYNES^® 230^® alloy. Material property values at very high temperatures (i.e., above the liquidus temperature) are difficult to obtain, and most of the information found in the literature is related to fits, assumptions, simplifications, and effective properties used to simulate phase transformations and stirring of the weld pool. Therefore, we prefer to feed the model with the available data documented in Figure 2, and allow the code’s material module adjust the values itself based on the most appropriate trend.

The mechanical model framework used to represent the physical problem is presented in refs. [22,23] as follows:

The components of the total strain-rate tensor ${\dot{\epsilon }}_{ij}$ are presented by the sum of the elastic ${\dot{\epsilon }}_{ij}^e$ and plastic ${\dot{\epsilon }}_{ij}^p$ components:

$${\dot{\epsilon }}_{ij}={\dot{\epsilon }}_{ij}^e+{\dot{\epsilon }}_{ij }^p$$

(6)

The elastic part obeys Hooke’s law and reads as follows:

$${\dot{\epsilon }}_{ij}^e=C_{ijkl}{\dot{\sigma }}_{kl} , C_{ijkl}=\frac{1+\nu }{E} \left({\delta }_{ik}{\delta }_{jl}-\frac{\nu }{1+\nu }{\delta }_{ij}{\delta }_{kl}\right) ,$$

(7)

where $C_{ijkl}$ is the elastic stiffness matrix, ${\sigma }_{ij}$ is the stress tensor and a dot denotes the differentiation with respect to time t, E is Young’s modulus, ν is the Poisson’s ratio, and ${\delta }_{ij}$ is Kronecker’s delta. The Levy–Mises equation accounts for the plastic part of the strain-rate tensor as follows:

$${\dot{\epsilon }}_{ij}^p=\frac{3}{2} \frac{\dot{p}}{q_{vm}}{\sigma }_{ij}^{\prime }$$

(8)

Here, ${\sigma }_{ij}^{\prime }$ is the deviatoric stress tensor; the quantities

$$\dot{p}={\left(\frac{2}{3}{\dot{\epsilon }}_{ij}^p{\dot{\epsilon }}_{ij}^p\right)}^{1/2}$$

(9)

and

$$q_{vm}=\sqrt{\frac{3}{2}{\sigma }_{ij}^{\prime }{\sigma }_{ij}^{\prime }}$$

(10)

are the equivalent von Mises plastic strain rate and stress, respectively. The substitution of the above equations into the equilibrium Equation (11)

$$\frac{\partial {\sigma }_{ji}}{\partial x_j}=0$$

(11)

leads to the governing equation to be solved by means of the finite element method. The elastic-plastic problem is a boundary-value problem that requires essential boundary conditions at each boundary point. Burnet [24] has shown that, for a well-posed problem, that is, to prevent the rigid-body motion of the entire structure, enough displacements (BCs) must be specified. In this work, to represent the fixation of the specimen to the base plate, appropriate mechanical boundary conditions were applied using the ABAQUS-BCs module. The eight corners of the thick sheet (the macro model) were firmly fixed with the ABAQUS option *ENCASTRE, and the ABAQUS option* ZASYMM was used to prevent rigid-body motion in the sub-model. Of this manner, loading from the master study was transferred as temperature or/and displacement values to the boundary between the surfaces included in the sub-model, and those neglected.

In a coupled temperature-displacement finite element model, the total strain increment at a material point from t → t + Δt can be expressed as:

$$\mathsf{∆}{\epsilon }_{ij}=\mathsf{∆}{\epsilon }_{ij}^e+\mathsf{∆}{\epsilon }_{ij}^p+\mathsf{∆}{\epsilon }_{ij}^{th}$$

(12)

Here, the total strain increment derived from the classical incremental plasticity theory can be accounted for an additive decomposition of the elastic, plastic and thermal strain component increments. On the other hand, a fairly general formulation adopted by most of the commercial finite element codes is the von Mises yield function defined by

$$f=q_{vm}-{\sigma }_y=\sqrt{\frac{3}{2}{\sigma }_{ij}^{\prime }{\sigma }_{ij}^{\prime }}-{\sigma }_y$$

(13)

Here, ${\sigma }_y$ is the yield strength at temperature T, and q_vm is the von Mises stress. The yield criterion is given by

f < 0: Elastic deformation

f = 0: Plastic deformation

The use of the von Mises equivalent quantities implies the plastic isotropy (ideal plasticity) of the material. The C3D8T ABAQUS element type and appropriate built-in-type boundary conditions (thermal and mechanical) were used in the macro and meso models. From there, a transient nonlinear coupled temperature-displacement finite element model was performed to determine the evolution of temperature, cooling rate and residual stresses at the meso-scale level in a thick sheet AISI 310 laser welding test sample. The calculations of the distortions, stresses, and strains are computationally intensive, and they take significantly more computer time than the calculation of the transient temperature field.

The numerical simulations presented in the next section correspond to an autogeneous laser weld test (Figure 3) with a heat input of 120 J/mm (emitted laser beam power P_E = 1500 Watts, heat source speed v = 6.25 mm/s and thermal efficiency n_th = 0.5) deposited on an AISI 310 austenitic steel macro mesh.

3. Results and Discussions

The multiscale simulation of microstructure and distortion evolution into the melt-pool is an important issue. It is, at the same time, a heavy and costly task, giving way to the use of different approaches traditionally used in multidomain approaches, such as level set or phase field. This may include, for example, the influence of cooling rates on the microstructure (formation of different phases, and grain texture), and resulting macro-scale properties. This also encompasses the modeling of heat sources combining fluid flow, solid mechanics, gas or plasma interaction, in order to predict the occurrence of various defects (for instance, porosity, segregation, and cracks) and evaluate their influence on local properties and part behavior. Since differences in the measured and calculated residual stresses are attributed to the absence of accurate stress–strain constitutive relations, uncertainties in the calculations of thermal cycles and the use of coarse grids, in this work, an efficient model order reduction approach is addressed. The numerical results obtained along three different paths located in the fusion zone were calculated in coarse and fine meshes. Figure 4 shows an autogenous laser weld deposited on a macro mesh of 40 × 10 × 75 mm, composed of 30,000 C3D8T elements of size 1 × 1 × 1 mm, and 3 thermal cycles calculated at points N1, N2 and N3 located ouside the fusion zone. It is observed that the cooling rates decrease as the analyzed locations are further away from the fusion zone.

Figure 5 shows the rough Gaussian profiles of temperature vs. true distance along the documented paths.

To obtain smoother Gaussian profiles, a fine meso grid was constructed from the global macro model via the submodeling technique. The domain of the meso model is located symmetrically at the top in the middle of the plate (see Figure 6a). In this figure, a description of the meso model with dimensions of 10 × 5 × 1 mm constructed from a global macro model of dimensions 40 × 10 × 75 mm (element size = 1 mm) via the submodeling technique using the Shape → Cut → Extrude options available in ABAQUS CAE is presented. A meso mesh of 50,000 C3D8T elements (element size = 100 µ) was generated. As expected, smoother Gaussian profiles of the temperature distribution along Paths 1 and 2 were obtained and are plotted in Figure 6b,c, respectively. The temperature–distance curves were obtained at the step time of 4.756 s.

Figure 7 shows a meso weld-pool cross-section, thermal cycles and cooling rates at six reported locations. The liquidus (T_L) and solidus (T_S) temperatures are also shown. For AISI 310 steel, the liquidus and the solidus temperatures are 1402 °C and 1354 °C, respectively [25]. Figure 7c reveals high cooling rates (through the solidification temperature range) in the points located at the meso fusion zone ranging from 960 °C/s at point D6 to 2400 °C/s at point D1.

Figure 8 shows the meso-level residual stresses (normal, shear and mises) at a cooling time of 60 s along Paths 1 and 2 (see the paths in Figure 5 and Figure 6). The maximum observed value of 520 MPa corresponds to the normal stress S33 that follows the direction of the heat-source speed (longitudinal axes Z or 3). These high values explain the evolution of the distortion of the fusion zone due to the welding process. Macro residual stresses are commonly calculated with the finite element method using macro grids. However, rough residual stress profiles are commonly plotted due to the scarcity of numerical points. In the proposed meso model, as a fine mesh is used, smoothed profiles of residual stresses can be better graphed, as shown in Figure 8.

The present contribution is in line with the multiscale framework that covers the thermal and mechanical failure analysis of extended microstructural regions. The submodeling technique (a flexible strategy used to calculate plastic strains, temperature gradients, and residual stresses to any number of levels) was used to mesh a local portion of the part with a refined mesh, based on the interpolation of the solution from the initial coarse macro global model. In a future work, the presented meso-level residual stresses will be validated by a comparison to the experimental results performed under the same welding and boundary conditions. In addition, more numerical and experimental work will be performed to compare the best welding parameter combination to minimize residual stresses and the formation of stray and equiaxed grains. Depending on the experimental technique and magnification used, it is observed that dendrite secondary arm spacing belongs to the micro scale, cooling rate belongs to the macro scale, and residual stresses belong to the macro, meso or micro scale. Therefore, to find new expressions that link, for example, the grain substructure with the cooling rate that is imposed on the welding during solidification, the experimental and numerical work should be conducted at the same scale. At the present time, the connection between the scales is yet unclear [26]. An interesting overview of the state of the art in laser welding simulation can be found in refs. [27,28,29].

Figure 9a describes Path 3 (see also Figure 7a), which is used to calculate the stresses vs. true distance along Path 3 shown in Figure 9b. The stress evolution at location D1 during the first 10 s of the welding process is shown in Figure 9c. The point D1 located at the top surface is the starting point of Path 3 (coordinate 0.00), and the point D6 located at 0.5 mm from D1 in the Y (2) direction is the end point of Path 3.

After the process has started, it takes 4.756 s for the heat source to be positioned just above the meso-model domain located at the central part of the plate. At that time, maximum negative values of stresses appear (see Figure 9c). Figure 9c also shows the evolution of the stresses that change signs (from negative-to-positive values) as the cooling time increases.

The evolution of the meso-normal stresses S33 is shown in contour graphs in Figure 10. These stress contours evolve from 100% negative values at step time 3 s, to almost 100% positive values at step time 60 s. At this step time, the plate temperature has already cooled down to approximately 70 °C. It therefore means that only a few more seconds are needed to finish the cooling cycle, and thus the complete conversion to positive values. It is important to note that, in order to mesh the 40 × 10 × 75 mm plate with elements of a 0.1 × 0.1 × 0.1 mm size, a total of 30 million elements is necessary. With the present approach, one macro mesh of 30,000 elements (1 × 1 × 1 mm) and a meso mesh of 50,000 (0.1 × 0.1 × 0.1 mm) elements were enough to simulate the weld problem at the meso-scale level. Multiscale modeling is based on fundamental physical principles and experimental data. The goal is to predict the behavior and performance of complex materials across all relevant time-and-length scales. The challenge is tremendous due to the fact that, at the macro scale (centimeters), stresses arising from temperature gradients may be the controlling elements of the materials’ performance. At the microscale (tens of micrometers), defects, such as dislocations controlling the mechanical behavior, occur, while large collections of such defects, including grain boundaries and other microstructural elements, govern the mesoscopic properties (hundreds of micrometers). The net outcome of these interactions governs the continuum behavior that can be described as a constitutive law. To analyze the thermal and mechanical phenomena on the micro scale, the proposed meso-model based on the interpolation of the solution from an initial relatively coarse macro global model has to be considered as the global model for a subsequent micro submodel. In this manner, the cooling rate–microstructure relationship can be calculated at the same scale level.

4. Conclusions

• A new three-dimensional coupled temperature-displacement finite element meso model constructed using appropriate ABAQUS CAE submodeling options is proposed.
• High cooling rates (through the solidification temperature range) ranging from 960 °C/s to 2400 °C/s were observed when the calculations were made at the mesoscale level. These high cooling rates that arise during the formation of the weld pool originate the highest observed residual stresses that evolve throughout the weld during the entire heating and cooling cycles.
• To mesh the entire plate with elements of a regular size of 100 × 100 × 100 µ, a total of 30 million elements are necessary. With the present approach, one macro mesh of 30 thousand elements (1 × 1 × 1 mm) and a meso mesh of 50 thousand elements (100 × 100 × 100 µ) were enough to simulate the weld problem at the meso-scale level.