a new finite element model for welding heat sources

A New Finite Element Model for Welding Heat Sources JOHN GOLDAK, ADITYA CHAKRAVARTI, and MALCOLM BIBBY A mathematical model for weld heat sources based on a Gaussian distribution of power density in space is presented. In particular a double ellipsoidal geometry is proposed so that the size and shape of the heat source can be easily changed to model both the shallow penetration arc welding processes and the deeper penetration laser and electron beam processes. In addition, it has the versatility and flexibility to handle non-axisymmetric cases such as strip electrodes or dissimilar metal joining. Previous models assumed circular or spherical symmetry. The computations are performed with ASGARD, a nonlinear transient finite element (FEM) heat flow program developed for the thermal stress analysis of welds.* Computed temperature distributions for submerged arc welds in thick workpieces are compared to the measured values reported by Christensen ~ and the FEM calculated values (surface heat source model) of Krutz and Segerlind. 2 In addition the computed thermal history of deep penetration electron beam welds are compared to measured values reported by Chong. 3 The agreement between the computed and measured values is shown to be excellent. I. INTRODUCTION THE problems of distortion, residual stresses, and reduced strength of a structure in and around a welded joint are a major concern of the welding industry. These problems re- sult directly from the thermal cycle caused by the localized intense heat input of fusion welding. The high temperatures developed by the heat source cause significant metallurgical changes around the weld area of low carbon structural steels. The thermal history, particularly soaking time at high temperatures and cooling time from 800 to 500 ~ deter- mines the microstructure and mechanical properties for a given composition. The cooling time from 400 to 150 ~ is a controlling factor in the diffusion of hydrogen and the cold cracking of welds. Accurate predictions of residual stress, distortion, and strength of welded structures require an accu- rate analysis of the thermal cycle. The importance of a good model for the weld heat source in the analysis of the thermal cycle has been emphasized by several investigators. 12'4-9 After examining the performance of several models, a new weld heat source model is proposed that is not only more accurate than those now available but is the first one capable of handling cases that lack radial symmetry. In addition, the model smoothes the load vector which reduces the error and the computing costs of FEM analysis. The basic theory of heat flow developed by Fourier and applied to moving heat sources by RosenthaP ~ in the late 1930s is still the most popular analytical method for calcu- lating the thermal history of welds. As many researchers have shown, Rosenthal's analysis (which assumes either a point, line, or plane source of heat) is subject to serious error for temperatures in or near the fusion and heat-affected zones. In regions of the workpiece where the temperature is less than about 20 pct of the melting point, Rosenthal's solution can give quite accurate results. However, the in- finite temperature at the heat source assumed in this model and the temperature sensitivity of the material thermal prop- erties (a temperature independent mean value is assumed) *Developed in the Department of Mechanical and Aeronautical Engi- neering, Carleton University, Canada. JOHN GOLDAK and MALCOLM BIBBY are Professors, Department of Mechanical and Aeronautical Engineering, Carleton University, Ottawa, Canada, K1S 5B6. ADITYA CHAKRAVARTI is Research Engineer, AMCA Corporation, Ottawa, Canada. Manuscript submitted February 14, 1983. METALLURGICAL TRANSACTIONS B increases the error as the heat source is approached. The effect of these assumptions and others on the accuracy of temperature distributions from the Rosenthal analysis has been discussed in detail by Myers et al. u To overcome most of these limitations several authors have used FEM to analyze heat flow in welds. Since Rosenthal's point or line models assume that the flux and temperature is infinite at the source, the temperature distri- bution has many similarities to the stress distribution around the crack tip in linear elastic fracture mechanics. Therefore many of the FEM techniques developed for fracture me- chanics can be adapted to the Rosenthal model. Certainly it would be possible to use singular FEM elements to analyze Rosenthal's formulation for arbitrary geometries. This would retain most of the limitations of Rosenthal's analysis but would permit complex geometries to be analyzed easily. However, since it would not account for the actual distribu- tion of heat in the arc and hence would not accurately predict temperatures near the arc, this approach is not pursued here. Pavelic et al. 8 first suggested that the heat source should be distributed. He proposed a Gaussian distribution of flux (W/m:) deposited on the surface of the workpiece. The subsequent works of Andersson, 5 Krutz and Segerlind,: and Friedman 7 are particularly notable. In References 2 and 7 Pavelic's disc model is combined with FEM analysis to achieve significantly better temperature distributions in the fusion and heat-affected zones than those computed with the Rosenthal model. While Pavelic's 'disc' model is certainly a significant step forward, some authors have suggested that the heat should be distributed throughout the molten zone to reflect more accurately the digging action of the arc. This approach was followed by Paley 6 and Westby 4 who used a constant power density distribution in the fusion zone (FZ) with a finite difference analysis, but no criteria for estimating the length of the molten pool was offered. In addition, it is difficult to accommodate the complex geometry of real weld pools with the finite difference method. A nonaxisymmetric three-dimensional heat source model is proposed in this investigation. It is argued on the basis of molten zone observations that this is a more realistic model and more flexible than any other model yet proposed for weld heat sources. Both shallow and deep penetration welds can be accommodated as well as asymmetrical situations. VOLUME 15B, JUNE 1984--299 The advantages of the model are demonstrated by compar- ing it with the Rosenthal calculations, other FEM models, and experimental results. II. THEORETICAL FORMULATIONS The model proposed in this investigation is a 'double ellipsoid' configuration. It is shown that the 'disc' of Pavelic et al. s and the volume source of Paley and Hibbert 6 and Westby 4 are special cases of this model. In order to present and justify the double ellipsoid model, a brief description of the Pavelic 'disc' and of the Friedman 7 modification for FEM analysis is necessary. In addition, the mathematics of the disc is extended to spherical, ellipsoidal, and finally to the double ellipsoidal configuration. In this way both the physics and mathematics can be presented and discussed in a coherent manner. A. Model Considerations The interaction of a heat source (arc, electron beam, laser) with a weld pool is a complex physical phenomenon that still cannot be modeled rigorously. At this time little is known about the distribution of pressure from the arc source, the precise effects of surface tension, buoyancy forces, and molten metal viscosity. However, it is known that these factors combine to cause weld puddle distortion and considerable stirring. Because of these complexities, modeling of the fluid flow phenomena directly is not at- tempted in this presentation (or elsewhere for that matter). However, because of the arc "digging" and stirring, it is clear that the heat input is effectively distributed through- out a volume in the workpiece. The 'disc' model is more realistic than the point source because it distributes the heat input over a source area. In fact, for a preheat torch that causes no melting this may be a very accurate model indeed. However, in the absence of modeling the weld pool free boundary position, the ap- plied tractions, and convective and radiative conditions between the weld pool and the arc, some form of ideali- zation of the heat source is necessary to achieve a solution. The disc model does not account for the rapid transfer of heat throughout the FZ. In particular, it is not possible to predict the deep penetration FZ of an EB or laser weld with the surface disc model. A comparison of calculated thermal history data (disc model) with measured values during this investigation underscored the need for an 'effective volume source' such as the one suggested by Paley and Hibbert. 6 In addition, it was found necessary to generate a volume source with considerable flexibility, i.e., the double ellip- soid model. With less general shapes such as a hemisphere or a single ellipsoid significant discrepancies between the computed and measured temperature distributions could not be resolved. The size and shape of the "double ellipsoid" can be fixed, i.e., the semi-axes lengths, by recognizing that the solid- liquid interface is the melting point isotherm (assuming two phase effects are negligible). At the same time weld pool temperature measurements have shown that there is little superheating in the molten zone. l The accuracy with which the heat source model predicts the size and shape of the FZ and the peak temperatures is probably the most stringent test of the performance of the model. In this investigation it was found that the most accuracy was obtained when the ellip- soid size and shape were equal to that of the weld pool. The nondimensional system suggested by Christensen 1 can be used to estimate the ellipsoid parameters. In the Paley and Hibbert 6 "effective volume heat source" the power density is constant throughout the molten zone. This is unrealistic physically because the stirring velocity must decay to zero at the FZ boundary and rise to a maxi- mum at the arc-weld interface. It is undesirable mathema- tically because the step in the power density requires a fine mesh in a FEM analysis to obtain accurate results which is computationally unacceptable. In this investigation a Gaussian distribution is assumed centered at the origin of the heat source. Intuitively this is preferable both mathema- tically and physically. The results support this contention. B. Gaussian Surface Flux Distribution In the 'disc' model proposed by Pavelic et al.,8 the ther- mal flux has a Gaussian or normal distribution in the z-z plane (Figure 1): q(r) = q(O)e -cr2 [1] where: q(r) = surface flux at radius r (W/m 2) q(0) = maximum flux at the center of the heat source (W/m 2) C = concentration coefficient (m-2). r = radial distance from the center of the heat source (m) A simple physical meaning can be associated with C. If a uniform flux of mag_nitude q(0) is distributed in a circle of diameter d = 2~X/C, the rate of energy input would be rllV, i.e., the circle would receive exactly the energy from the arc. Therefore the coefficient, C, is related to the source width; a more concentrated source would have a smaller diameter d and a larger value of C (Figure 1). Experiments have shown that a significant amount of heat is transferred by radiation and convection from the arc directly to the solid metal without passing through the molten pool. Based on this observation, Pavelic et al. s de- veloped a correlation showing the amount and the distri- bution of this heat over the solid material. In their study, I Arc Flame Spread ~ , ~ / j ' q L dH _ ~. . End of Flame !:::::::." I f Heot Distribution e I :=- c 2 ~ c 3 d H --- Hot Spot (Dia.) c = Concentration Coeff=cient Fig. 1 --Circular disc heat source (Pavelic et al.S). 300--VOLUME 15B, JUNE 1984 METALLURGICAL TRANSACTIONS B provisions were made for convective and radiative losses from the heated plate to the surroundings as well as variable material properties. Friedman 7 and Krutz and Segerlind z suggested an alter- native form for the Pavelic 'disc'. Expressed in a coordinate system that moves with the heat source as shown in Fig- ure 2, Eq. [1] takes the form: q(x, ~) = 3--~Q e-3X2/C2e-3~2/c2 [2] ,.B.C 2 where: Q = energy input rate (W) c = is the characteristic radius of flux distribution (m) It is convenient to introduce an (x, y, z) coordinate sys- tem fixed in the workpiece. In addition, a lag factor ~" is needed to define the position of the source at time t = 0. The transformation relating the fixed and moving coordinate systems is: s t= z + v ( r - t) [3] where v = the welding speed (m/s). In the (x, y, z) coordi- nate system Eq. [2] takes the form: q(x, z, t) = 3Q e_3Xz/C2e_3tz+v(r_t)]%cZ [4] 7"gC 2 fo rx 2+ ~:2c 2,q(x,~,t ) =0. To avoid the cost of a full three-dimensional FEM analysis some authors assume negligible heat flow in the longitudinal direction; i.e., OT/Oz = 0. Hence, heat flow is restricted to an x-y plane, usually positioned at z = 0. This has been shown to cause little error except for low speed high heat input welds. 5 The disc moves along the surface of the work- piece in the z direction and deposits heat on the reference plane as it crosses. The heat then diffuses outward (x-y direction) until the weld cools. C. Hemispherical Power Density Distribution For welding situations, where the effective depth of pene- tration is small, the surface heat source model of Pavelic, Friedman, and Krutz has been quite successful. However, for high power density sources such as the laser or electron beam, it ignores the digging action of the arc that trans- ports heat well below the surface. In such cases a hemi- spherical Gaussian distribution of power density (W/m 3) / Fig. 2--Coordinate system used for the FEM analysis of the disc model according to Krutz and Segerlind. 2 would be a step toward a more realistic model. The power density distribution for a hemispherical volume source can be written as: 6 "~/ 3Q _312/c2 _ 3;.2/c2 _ 3~.2/c2 q(x,y,~) = c3~V,--~e e e [5] where q(x, y, ~) is the power density (W/m3). Eq. [5] is a special case of the more general ellipsoidal formulation de- veloped in the next section. Though the hemispherical heat source is expected to model an arc weld better than a disc source, it, too, has limitations. The molten pool in many welds is often far from spherical. Also, a hemispherical source is not appropriate for welds that are not spherically symmetric such as strip electrode, deep penetration electron beam, or laser beam welds. In order to remove these constraints, and make the formulation more accurate, an ellipsoidal volume source is now proposed. D. Ellipsoidal Power Density Distribution The Gaussian distribution of the power density in an ellipsoid with center at (0, 0, 0) and semi-axes a, b, c paral- lel to coordinate axes x, y, s c can be written as: q(x, y, st) = q(O)e-~2e-B~2e-C~2 [6] where q(0) is the maximum value of the power density at the center of the ellipsoid. Conservation of energy requires that: 2Q = 2~1V1 = 8 q(O)e-'W-e-8~2e-C~2 dx dy dst [7] 000 where: 77 = heat source efficiency V = voltage I = current Evaluating Eq. [7] produces the following: q(0)TrV'~ [81 2Q= 2 Q ~/-A-ffC q(0) - 7 r~ [9] To evaluate the constants, A,B, C, the semi-axes of the ellipsoid a, b, c in the directions x, y, s t are defined such that the power density falls to 0.05q(0) at the surface of the ellipsoid. In the x direction: q(a, O, O) = q(O)e -A~2 -= 0.05q(0) [10] Hence Similarly In 20 3 A = a2 a2 [111 3 B -~ - - [121 b 2 3 C -~ -- [131 C 2 Substituting A,B ,C from Eqs. [11] to [13] and q(0) from METALLURGICAL TRANSACTIONS B VOLUME 15B, JUNE 1984-- 301 Eq. [9] into Eq. [6]: q(x, y, ~) = 6X/'3Q e_3X2/a2e_3y2/b2e_3~%c2 abcTrV'-~ [ 14] The coordinate transformation (Eq. [3], Figure 2) can be substituted into Eq. [14] to provide an expression for the ellipsoid in the fixed coordinate system. q(x, y, z, t) = 6X/3Q e_3X2/a2e_3yZ/b2e_3[z+v(r_t)]2/c 2 [15] abcTr~,~ If heat flow in the z direction is neglected, an analysis can be performed on the z-y plane located at z = 0 which is similar to the 'disc' source (Figure 2). Where the ellipsoidal source intersects this plane the power density is calculated for each time increment. E. Double Ellipsoidal Power Density Distribution Calculation experience with the ellipsoidal heat source model revealed that the temperature gradient in front of the heat source was not as steep as expected and the gentler gradient at the trailing edge of the molten pool was steeper than experimental experience. To overcome this limitation, two ellipsoidal sources are combined as shown in Figure 3. The front half of the source is the quadrant of one ellipsoidal source, and the rear half is the quadrant of another ellipsoid. The power density distribution along the ~ axis is shown in Figure 3. In this model, the fractions fy andfr of the heat deposited in the front and rear quadrants are needed, where fl + fr = 2. The power density distribution inside the front quadrant becomes: 6 %/-3f/2 e_ 3x2/a2e_ 3y2/b2e_ 3[z+v(~,_t)]2/c2 q(x, y, z, t) - abcTr%/--~ [ 16] Similarly, for the rear quadrant of the source the power density distribution inside the ellipsoid becomes: 6 X/3frQ e_3x2/a2e_ 3y2/b2e_ 3[x +v(r_t)]%c2 q(x, y, z, t) - abcTrX ~ ~ [17] ul ms Fig. 3--Double ellipsoid heat source configuration together with the power distribution function along the ~ axis. In Eqs. [16] and [17], the parameters a ,b ,c can have different values in the front and rear quadrants since they are independent. Indeed, in welding dissimilar metals, it may be necessary to use four octants, each with independent values of a, b, and c. FEM CALCULATIONS III. EVALUATION OF THE DOUBLE ELLIPSOID MODEL In order to minimize the computing cost the initial analy- sis was done in the plane normal to the welding direction as shown in Figures 4 and 5. Thus, heat flow in the welding direction was neglected. The above simplification is accu- rate in situations where comparatively little heat flows from the arc in the welding direction. This is reasonable when the arc speed is high. An estimate of the effect of this approxi- mation has been given by Andersson 5 who argues that the errors introduced by neglecting heat flow in the direction of the moving electrode are not large, except in the immediate vicinity of the electrode. A. Verification of the Model In order to demonstrate the flexibility and assess the va- lidity of the double ellipsoidal heat source model two quite different welding situations were considered. The first case analyzed was a thick section (10 cm (4 inches)) submerged arc bead on plate (low carbon structural steel - 0.23 pct C) weld shown schematically in Figure 4. The welding condi- tions are contained in the figure. Christensen ''4 reported a I Ib y / / / / / / FZ HAZ 1480 =C 4' 7230C- / Experimental I x , Voltage = 32.9% Current = lI70A Welding speed = 0.005 ra/s Ef f ic iency n = 0.95 / El l ipso ida l Parameters (Equat ions 16,17 ) Semi -axes (cm) a = 2.0 b = 2.0 l c,- 1.5 c2= 3.0 Heat Input F rac t ions ff= 0.6 IOcm -- fr 1.4 Boundaries FZ (xl) 1.6 1.4 HAZ (x2) 2.1 2.1 Calculated (cm) Ex)erimental (cm) Fig. 4--Experimental arrangement and FEM mesh for the thick section bead on plate weld reported by Christensen. 302--VOLUME 15B, JUNE 1984 METALLURGICAL TRANSACTIONS B / / / 1 / J Weld Pro f