Simulating Convective-Radiative Heat Sinks Effect by means of FEA-based Gaussian Heat Sources and Its Approximate Analytical Solutions for Semi-Infinite Body

FEA-based Gaussian density heat source models were developed to study the effect of convective and radiative heat sinks on the transient temperature field predicted by the available approximate analytical solution of the purely conduction-based Goldak’s heat source. A new complex 3D Gaussian heat source model, incorporating all three modes of heat transfer, i.e., conduction, convection and radiation, has been developed as an extension of the Goldak heat source. Its approximate transient analytical solutions for this 3-D moving heat source were derived and numerically benchmarked with the available measured temperature & weld pool geometry data by Matlab programming with ~5 to 6 times faster than FEA-based simulation. The new complex 3D Gaussian heat source model and its approximate solution could significantly reduce the computing time in generating the transient temperature field and become an efficient alternative to extensive FEA-based simulations of heating sequences, where virtual optimisation of a melting heat source (i.e. used in welding, heating, cutting or other advanced manufacturing processes) is desirable for characterisation of material behaviour in microstructure evolution, melted pool, microhardness, residual stress and distortions.


Introduction
Current increased application of virtual heating and welding simulation by laser melting heat source, requires more precise and rapid transient solutions for moving heat sources. To model the melting process accurately, it is critical to acquire sufficient knowledge on transient temperature field during manufacturing process to control and optimise component thermal history and metallurgy, e.g., microstructure evolution, residual stress and distortion, that influencing the mechanical properties of the fabricated products.
In general, for complex virtual heating and welding simulations, the FEA-based models are employed to generate the transient temperature required. However, if extensive heating sequences of the welding process and more details of the melting zone simulation are required, the FEA model would become bigger due to larger number of elements and heating sequences, and the FEA computing time could become an issue even though the PC speed has been rapidly improved for the last two decades. Therefore, available approximate analytical solutions could still be a very useful tool, not only for benchmarking the FEA solutions but also become an efficient alternative for the nominated typical applications.
Very first analytical solution of transient temperature field for a point heat source was traced back in 1941 by Rosenthal [1]. This solution found its broad application in welding, cutting and heating processes despite its singularity limitation at the heat source centre. The Rosenthal very first point heat source solution was purely based on conduction only and did not taking any heat loss by heat convection or radiation. However, the effect of heat convection was reported to be included into the analytical solution for the unit instantaneous plane source for semi-infinite body, which was published by earlier authors such as Brian (1891), Carslaw (1902) & Hobson (1887) as reported in Carslaw and Jaeger textbook [2].
Few decades later, Eager and Tsai [3] developed an analytical solution for 2D-Gaussian distributed flux densify heat source, which better simulated the shape of the welding arc and partly reducing the singularity issue of the Rosenthal solution. Two years later Goldak and his co-authors introduced a 3D-Gaussian distribution flux density heat source [4] and used Finite Element Method (FEM) to calculated transient temperature field in a welded plate.
Analytical solution for the Goldak's heat source was not found until 1999, when Nguyen and his co-workers [5] developed a transient analytical solution for this 3D-Gaussian density heat source in semi-infinite body. The authors carried out experimental transient temperature measurements in vicinity of weld pool region and numerical modelling to prove the validity of their analytical solutions for the double ellipsoidal heat sources.
Nguyen et. al. original analytical solution for 3-D moving Gaussian density heat source in semi-infinite body was also based on heat conduction-only, assuming the effect of heat conduction was predominant. Nguyen and his co-authors [6] later developed an approximate analytical solution for double ellipsoidal heat source in finite thick plate, i.e., incorporating the effect of plate thickness into their transient solution by means of error functions.
As the need for more efficient and accurate 3D heat source modelling has arisen, the researchers continue to improve the current conduction-only analytical solutions to include the temperature dependent material properties [7][8] and effects of heat convection and radiation [9][10] as well. Guo and Dai [7] and Wenji et al. [8] included temperature-dependent properties, i.e., thermal conductivity, density and specific heat, in their attempt to improve current conduction-only analytical solutions for transient temperature field, that were derived for the constant material properties, due to moving Gaussian heat source.
Pathak et al. [10] proposed a new method to account for convection and radiation heat losses from the boundary surfaces during simulation of Multipass welding thermal cycle. The authors used the fluid dynamics Nusselt number (Nu) to derive effective heat transfer film coefficient for both heat radiation and convectionon heat loss. The authors compared their proposed effective heat transfer film coefficient with the Vinokurov's empirical model (heff = 24.1x10 -4 T 1.61 where  is emissivity and T is temperature). It was reported that their effective heat transfer film coefficient is higher than the Vinokurov's model for temperatures that below 600 o C, but comparable to the Vinokurov's model for temperatures higher than 600 o C. The authors used 1-D transient solution of a point heat source [2], with equivalent heat lost amount, as temperature reduction due to the convection and radiation heat losses.
Veldman et al. [11] used Green function method to solve 2D conduction problem of infinite thin plate with a finite thickness and constant material properties, incorporating the heat convection boundary condition for the heat loss at top and bottom of the thin plate. The authors used a fundamental solution of the homogeneous differential 2D heat conduction equation in thin plate, Duhamel's principle and Taylor series to obtain an approximate transient temperature solution for a top surface heat source in infinite thin plate. Veldman et al. [11] managed to include a heat convection term (exp(-h.t), where h is heat convection coefficient and t is time) into their approximate analytical solution for thin plate.
Numerical techniques such as the finite element (FE) method are commonly used, particularly for simulating of heating process in complex geometries, to obtain transient temperature field for various heat sources [10,[12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29]. The FEA method requires calibration of the heat source parameters to get a solution with acceptable accuracy for predicting the transient temperature field in the heated products and melted pool geometry. The authors included the effect of convection and radiation into their transient temperature field modelling via FE boundary conditions as these effects are deem significant [23][24]. Usually, the FEA boundary conditions, imposed on the external surfaces of the spatial domain, using FEA packages such as ANSYS [14][15][16][17][18][25][26][27] or ABAQUS [28, 29]. Dutta and Narendranatah [23] studied the heat loss due to convection and radiation. The authors proposed a method to estimate total heat loss and compared their results with Vinokurov's empirical model due to radiation heat loss. They concluded that the radiation heat loss is dominating at higher temperature range (400 o C to 1300 o C) whereas convective heat loss is dominating at lower temperatures (30°C to 350°C).
Lemkeddem et al. [24] developed a numerical model to determine transient temperature distribution for the laser welding of two identical TA6V thin plates, taking account the energy lost by radiation and convection. They concluded that it is necessary to include energy losses in the welding area to improve laser welding efficiency. The energy lost due to radiation and convection to the surrounding environment would affect the total welding heat input into the welding region and significantly influence on the temperature field.
The analytical models are continuously improved overtime with modifications of the heat source geometry and boundary conditions [1,[3][4][5][6]. However, the heat loss in the analytical models due to convection and radiation were not included in these analytical models due to complexity of the solutions involved. The existing analytical solutions for the Gaussian heat source are conduction-only, hence, to our best knowledge, there is no 3D heat source and associated analytical solution that includes all three modes of heat transfer (i.e. heat conduction, convection, and radiation), is available for use to date. This paper aims to investigate and include the effects of two convective and conductive heat sinks on the transient temperature field, subjected to 3-D Gaussian heat sources by means of FEA modelling and simulation.
A new complex 3D Gaussian heat source model, incorporating all three modes of heat transfer (conduction, convection, and radiation), will be developed as an extension of the Goldak heat source. Its approximate transient analytical solutions for this complex 3-D moving heat source will also be derived and numerically benchmarked with the available measured temperature and weld pool geometry data and FEA-based simulations by Matlab programming.
Firstly, the new complex 3D Gaussian heat source model was theoretically developed based on the heat balance in semiinfinite body that including all three modes of heat transfer. The convective and radiative terms in the new complex heat source were effective heat sinks that corresponding to equivalent amount of heat loss from the infinite-body subjected to a moving Gaussian heat source.
Subsequently, the FEA modelling was used to generate the transient temperature field due to the new complex 3D Gaussian heat source, to demonstrate that the new complex heat source could be able to provide comparable transient temperature field as generated by FEA with convection and radiation boundary conditions. The new complex Gaussian heat source parameters were established by utilising the available transient measured temperature data from Nguyen et al. experiment [5]. A FEA model simulating the Nguyen et al. experiment was developed in Abaqus for this purpose.
Secondly, based on the Green function approach, superposition principle and Taylor series approximation, the approximate transient analytical solutions for the new complex Gaussian moving heat sources were developed for semiinfinite body.
The transient temperature results were then generated from these approximate analytical solutions using Matlab codes. The solution of transient temperature field is then compared with the FEA-based generated transient temperature field and Nguyen et al. experimental measured data to validate the new complex Gaussian heat source solutions. Goldak, et al. [4] developed Gaussian distribution double ellipsoidal heat source as shown in Fig. 1. There are two halves of un-even ellipsoidal heat source with different heat density parameters along the travelling x-axis, c hf and c hr , and with the same heat source density parameters along two other axes y and z, i.e. a h and b h. The volumetric heat flux, i.e., flux density, Q(x,y,z) at a point (x,y,z) within the semi-ellipsoid is given by the following equations (1a) and (1b) for the front half and rear half of the double ellipsoid, respectively as:

New Complex Heat
where a h , b h , c hf and c hr are ellipsoidal heat source parameters (as shown in Fig. 1) x, y, z is moving coordinates of the heat source. Q(x,y,z) is heat density at a point (x,y,z) V, I and  are welding voltage, current and arc efficiency, respectively. Q = arc heat input (Q = IV); r f , r r = proportion coefficient representing heat apportionment in front and rear half of the heat source, respectively, where r f + r r = 2.
Due to the condition of continuity of the volumetric heat source, Q(x,y,z) given by Equations (1a) and (1b) must be equal at the x = 0 plane, hence, r f and r r must satisfies r f /c hf = r r /c hr . Subsequently, these two coefficients are determined as r f =2c hf /(c hf + c hr ) and r r =2c hr /(c hf + c hr ).
This double ellipsoidal distribution heat source is conduction-only heat source, and it is defined by five unknown parameters, i.e., the arc efficiency  and four other ellipsoidal geometric parameters: a h , b h , c hf and c hr , which are linearly proportional to the measured weld pool geometry parameters. This means the Goldak's parametric density heat source requires a thorough calibration with the measured temperature data and weld pool geometry before it could be used for various thermal heat transfer applications.

Past analytical solutions for single and double ellipsoidal heat sources
Analytical solutions for the conduction only Goldak's moving 3D heat source had been obtained by Nguyen et al. where a = thermal diffusivity (a = k/ (ρ c); c = isobaric specific heat; k = thermal conductivity; ρ = mass density; t, t' = time for moving heat source with uniform speed v from time 0 to t).
Noting that when chf = chr, the solution for double ellipsoidal heat source in semi-infinite body reduces to the solution for single ellipsoidal heat source in semi-infinite body as expected.

Heat convection due to analytical solutions of the Goldak's heat sources
If the amount of heat loss (heat sink) by heat convection to the environment (i.e., the second half of the semi-infinite body) were subtracted from the total heat input, from the ellipsoidal heat sources, into semi-infinite body, we would obtain the effective conduction-convection solution for the complex Gaussian conductive-convective heat sources.
As the effect of the heat losses due to convection, radiation and other heat sinks (e.g. evaporative cooling of the welding arc due to melted metal evaporation et al. ) was effectively modelled by heating/welding arc efficiency coefficient (), it is reasonable to assume that the heat sink surface 2D temperature distribution from conduction-only solution, with effective heat input (Q = IV) as used in the Goldak's heat source, could be used to estimate the heat loss due to heat convection and radiation. This assumption is used as part of our assumption to derive approximate analytical solution that included the heat convection and radiation losses from the conduction only Goldak's heat sources.

Heat convections sink via semi-infinite body surface
Let us assume h is temperature independent heat transfer convection film coefficient on entire semi-infinite surface of the semi-infinite body and (dx.dy) is elemental area on that surface. The corresponding elemental flux heat sink to the environment through this area is expressed as where dT k-h (x,y,0,t) = temperature difference between the body surface and bulk environment temperature due to both heat conduction and convection.
Similarly, it is assumed that the elemental temperature differential at a time instant due to both heat convection and conduction (dT k-h (x,y,0,t-t') ) is approximate the same as that due to the effect of heat conduction only dT k (x,y,0,t-t') when subjected to effective heat input used by means of heating arc efficiency coefficient () as described above. Subsequently, Eq. (3) can be rewritten for convective heat sink rate as where dT k (x,y,0,t)/dt can be obtained from Eq. ) . . exp (− where Q max,h is equivalent peaked heat sink flux due to heat convection via top surface of the semi-infinite body and  h is thickness of the thin film through which convective heat sink is effective.
Subsequently, the total amount of heat lost through top surface of the semi-infinite body due to heat convection, in front and back half of the double surface heat sink flux, is given by The infinite integrals of equations (8a) & (8b) were solved and these reduced to The total heat sink Q h is sum of the heat sink amount from both halves (eqs. (9a) & (9b)) as On the other hand, this amount of heat sink rate due to convection can be calculated as From Eqs. (10) and (11) Where , = 3√3 √ Further simplifying Eq. (13), noting r f + r r =2 gives . Substituting Eq. (14) into (7a) and (7b) give the heat sink equations due to 3D Goldak's heat source as

Heat radiations sink via top surface of semi-infinite body
Similarly, the elemental heat sink approach is again applied here for the contribution of the radiation heat sink. The heat sink flux due to radiation, ( , , ), via top surface of the semi-infinite body is expressed as where T k, rad (x,y,0,t) and T o are the top surface temperature due to effect of conduction and radiation at time t and constant initial temperature, respectively;  = Stephen-Boltzmann constant,  = emissivity.
It is also assumed here that the temperature differential at a time instant due to both heat radiation and conduction (where dT k, rad (x,y,0,t)) is approximate the same as that due to the effect of heat conduction-only (dT k (x,y,0,t)) ) when subjected to effective heat input used by means of heating arc efficiency coefficient () as described above. Subsequently, Eq. (14) can be rewritten as Recall equation (6a) for single ellipsoidal heat source and noting that , = 3√3 √ : ( , ,0, ) = .
The T k (x, y, 0, dt) term can be approximately expressed as (17) gives Based on the Taylor series expression for exponential function exp(-x 2 ) as where n = 1, 2, …., n is natural number, we could express two exponential terms in Eq. (21) by taking only first two terms in the Taylor series as can be rewritten as Let's also consider that the radiative heat sink density flux due to single ellipsoidal heat source has following form: where Q max,rad is equivalent peaked heat sink flux due to heat radiation sink via top surface of the semi-infinite body and  rad is thickness of the thin film through which radiative heat sink is effective.
Subsequently, the total amount of heat sink rate through top surface of the semi-infinite body due to heat radiation sink flux is given by Eq. (25a) can be reduced to On the other hand, this amount of heat sink flux rate due to radiation can be calculated as Let's us solve two infinite integrals in equation (26): Using Integration by part rule and noting that: Then Eq. (24) can be rewritten for radiation heat sink density flux rate as: Equation (32) can be extended for double ellipsoidal heat sink flux as

Complex heat source including two heat sinks in semi-infinite body
It's reasonable to assume that the film thickness of effective convective and radiative heat sinks is the same, , we could express the convective and radiative density heat sinks for double ellipsoidal density heat source, noting that , = 3√3 √ , as: where k rad , as expressed by Eq.(23b) above, is considered as the correction factor for approximate solution for radiative heat sink rate, subsequently, the complex 3D density heat source can be expressed as putting r f = r b = 1 and c hf = c hb , Eqs. (36a) and (36b) would be reduced to equations for complex single ellipsoidal heat source as or in a shorten form as

Approximate analytical solutions for the complex heat source-heat sink
Recall the Fourier-Biot's differential equation for heat conduction: where Q v = heat generation or heat sink density (W/m 3 ).
As the point heat source solution is fundamental solution to Fourier-Biot's differential equation when Q v =0, the solution for Eq. (39) could be derived from the fundamental solution by volume integral as: Noting that Eq. (40) has been solved by Nguyen et al. [5] to obtain the temperature solution for both single and double ellipsoidal heat source, as expressed by Eqs.
Substituting the values of Q max , Q h,max and Q rad,max into Eqs. (41b), (41c) and (41d) and further simplifying gives

.1 Bead-on-plate Finite Element model for validation
In this section, a numerical model developed in ABAQUS [30] is described and used to simulate the experimental tests performed by Nguyen et al. [5]. The Goldak's double ellipsoidal and the new complex double ellipsoidal heat sources were used in the FEA model to generate the transient temperature fields and weld pool geometry for the purposes of this work. A numerical-experimental correlation of the available transient temperature distribution data and the weld bead sizes was executed to determine the suitable calibrated heat source parameters. This task is aimed to validate the proposed new complex heat source and its approximate analytical solution that incorporating the effect of the convection and radiation.

Model description
The specimen size used for the FEA model was the same as that used for the transient temperature measurement with plate dimensions of 240 x 240 x 20 mm by Nguyen et al [5]. The test specimen consists two perpendicular linear segments of weld beads, run on top of the plate, as shown in Fig. 2. This figure shows the coordinates of the three points A, B and C, used to compare the transient temperature around the weld toe, weld root and at the corner of the weld paths connection where the two weld paths change their orientation, respectively. The welding parameters used in the Abaqus FEA model are the same as used in the experimental test, i.e., U = 26 V, I = 230 A, v = 5 mm/s and the efficiency η = 0.80.

Material properties
The same heat transfer material properties used in the analytical analyses and the experimental test were used in the FEA model. The base material and weld were assumed to have the same thermal properties. The following material properties were used, i.e., material isochoric specific heat c =c p = 600J/kg °C, thermal conductivity k = 29 W/m/°C and density ρ = 7820 kg/m 3 .

Boundary condition and loading
The boundary conditions considered in the FEA analyses were the convection and the radiation films. Figure 3a, 3b and 3c show three different boundary conditions used in various FEA analyses in this work. Figure 3a shows FEA boundary condition for a base load case, which was used to determine the effective Goldak's heat source parameters used for the bead-on-plate experiment. This specimen is subjected to typical laboratory convection in still air (h = 5 W/m 2 / o C) and emissivity = 0.67 for carbon steel plate. This condition was used for radiation boundary conditions in the developed FEA model for the base load case.  Therefore, all three modes of heat transfer (i.e., conduction, convection, and radiation), as described by the new complex heat source equations, was implemented using body heat flux option and DFLUX code with no additional boundary conditions, as shown in Fig. 3c.   The radiation effect on the transient temperature profile was analysed in this work with Stefan-Boltzmann constant used is 5.67E-8 W/m 2 /K 4 . The emissivity coefficient was considered in all surfaces was of 0.67 for grey body in the base load case and ambient temperature of 20 °C was used.
The arc moving at the constant travelling speed of 5 mm/s and the Goldak's parameters, determined by benchmarking with the transient measured temperature data [5], that define the volumetric heat flux used for bead-on-plate welding is simulated and implemented through a written user subroutine DFLUX. The DFLUX codes were developed with defined the coordinates of the start-stop positions, welding path trajectory, velocity of the heat source and welding arc heat input that were implemented for both conduction-only and for the new complex heat source solution in this work.
The parameters of the double ellipsoidal heat source used in the FEA model, were determined from calibration process of the measured transient temperatures at three selected points A, B, and C and the weld pool size data [5], these were c hf =7 mm, c hr =14 mm, a h =8.5 mm and b h =1.0 mm.

Mesh model
Mesh model for 3D FEA analyses was generated by ABAQUS with cubic elements of 1 mm for the weld bead areas to capture the high-temperature gradients in weld pool and HAZ regions. The element size was gradually increased towards the transversal directions of the weld bead. Linear heat transfer brick elements with eight nodes, DC3D8-type elements, from ABAQUS library program were used. A total of 77,030 elements were generated for the mesh model (as in Fig. 4).

Heat source parameter calibration using measured temperature data and weld pool geometry
The existing transient temperature data and weld pool geometry of the bead-on-plate experiment by Nguyen et al. [5] were used to calibrate the FEA model by matching the transient temperature results generated for three selected points A, B and C (see Fig. 2) to determine a suitable heat source parameters for the double ellipsoidal heat source for the bead-onplate FEA model use.
The FEA model with conduction only Goldak's heat source and boundary condition of natural convection film of h = 5 W/m 2 / o C and radiative film of 0.67 emissivity for all surrounding surface at room temperature 20 o C, was used for the experimental condition in laboratory environment.

Transient FEA temperature results vs transient temperature data [5]
Initially, parameters for double elliptical heat source were estimated as c hf =7 mm, c hr =14 mm, a h = 8.0 mm & b h =1.2 mm and used as the volumetric heat source in the FEA model. Convection heat transfer film of h = 5 W/m 2 / o C and radiative film of emissivity of 0.67 for carbon steel at 20 o C ambience temperature was used to generate the transient temperature field, simulating the two-path Nguyen's bead-on-plate experiment [5]. The available temperature data measured at three selected points (A, B & C in Fig. 2), was used to calibrate the FEA generated transient temperatures to determine the best suitable heat source parameters that give as small deviation as possible.
Subsequently, the heat source parameters of c hf = 7 mm, c hb =14 mm, a h = 8.5 mm and b h =1.0 mm were determined to generate the transient FEA temperatures that are within ~6% compared with the measured temperature data shown in Fig. 5. The deviation between the simulated and measured transient temperature data is attributed to spatter mode of the welding arc and numerical and experimental errors, which is acceptable for the purpose of this study, i.e. to include the convective and radiative effects in an approximate transient analytical solution for the conduction-only double ellipsoidal heat source by Nguyen et al. [5]. A sensitivity study was carried out on the effect of heat convection and radiation films on the maximum transient temperatures, induced at the three selected points A, B and C, when all three modes of heat transfer were included into a FEA-based simulated conduction-only heat source. Figure 6 shows there are significant reduction in the maximum temperatures generated, when both heat convection and radiation sinks were included, compared to the conduction-only temperature results. These FEA simulated results confirmed the significant of the heat convection and radiation, which should be included in the numerical modelling of any transient temperature field due to Gaussian heat source to improve the accuracy of the simulation. Figure 7 shows how the weld pool geometry shape (i.e., weld pool width, length and depth) is subjected to the effect of convection and radiation. It can be observed that for the conduction-only heat source, the weld pool geometry is the largest as there is no heat loss by convection and radiation. However, as both effect of convection and radiation become more pronounced, the weld pool shape become smaller in all three dimensions, i.e., weld pool width decreases from 13.2 to 12.7 and 12.5 mm when the heat convection and radiation was included as normal laboratory conditions. The weld pool shape was reduced further when more intensive heat convection (i.e., in windy condition when heat convection coefficient was increased to h = 100 W/m 2 / o C) and when steel plate was in oxidised condition with emissivity increased to be 0.79 (instead of 0.67 at normal condition). The weld depths and lengths for three conditions were also decreased accordingly as shown in Table 1. In Fig. 7, the melting zone is characterized by the grey area, which indicated temperatures higher than the assumed melting temperature of 1,450 °C.    Figure 8a and 8b show the FEA simulated top plan view and transversal cross section of the weld pool, respectively, compared with the available measured experimental data [5]. It could be seen from the Fig. 8a that there is a very good agreement between the top plan view of the simulated weld pool shape and the measured weld bead size data, corresponding to the front half of the heat source as well as front part of the rear heat source. However, the measured tail of the weld bead is deviated from the FEA simulated weld pool, i.e., by ~4 mm. Figure 8b shows the maximum depth and width of the weld pool is in good agreement with the measured weld bead data. However, the finger-like shape of the weld bead geometry measured by Nguyen's et al. experiment [5] (Fig. 8) is typical spray mode of molten metal transfer from electrode to the work piece in Gas Metal Arc Welding (GMAW). The welding process used for this experiment was 80% Argon and 20% CO 2 , the welding current and voltage was 230 A and 26V, respectively. These welding parameters used were just above the critical current for spray-arc mode for 80% Argon and 20% CO 2 process. This explained the finger-like shape of the weld bead measured.

FEA bead-on-plate weld pool geometry vs. experimental weld pool data [5]
Therefore, the deviation between the measured and simulated weld pool (as shown in Fig. 8) could be justified, because the Gaussian density heat source is only valid for the first two modes of metal transfer in GWAW, i.e., the short-circuit and globular transfer mode.

Complex heat source with two heat sinks included in the analytical solution
In this section we are going to demonstrate the two points below, by using the FEA modelling:

Numerical implementation description
The new complex density heat source equations, as expressed by Eqs. (15a) and (15b) for convective effect and by Eqs.
(33a) and (33b) for the radiative effect, were implemented using Abaqus user subroutine DFLUX as body heat density, normally implemented for the conduction-only Gaussian heat source. The convective and radiative heat sink density were included in the DFLUX user-subroutine.
By changing the heat sink parameters (i.e. namely h, emisivity,  h ,  rad and k rad ) we would be able to include the convective, radiative effects or both into the complex double ellipsoidal heat sources, which would in turn produce similar transient temperature results as if these effects were implemented by the FEA-based heat convective and radiative boundary conditions, available in FEA packages such as in Abaqus [30] or ANSYS [31].

FEA simulated temperatures against those generated by the complex heat source equations
A combined effect of two heat sinks, i.e., heat convection film of h = 100 W/m 2 / o C and radiative film with emissivity of 0.67 at ambience temperature of 20 o C and with the same film thickness delta = delta2 = 0.001 mm, radiative correction factor k rad = 0.30, was used to model the new complex density heat source model. The DFLUX was used in the developed bead-on-plate FEA model to generate the simulated transient temperature results and compared with the temperatures that generated by the same bead-on-plate FEA model, with conduction-only Goldak's heat source subjected to the same two FEA-based heat sink conditions.
The transient results generated by both methods for three selected points A, B and C are in very good agreement as shown in Fig. 9. Thus, demonstrated that the two methods used to generating the transient temperature fields are equivalent as the transient temperature results almost identical to each other.
Thus confirmed the newly derived effective complex density heat source equations could be used to model the effect of heat sinks in the Gaussian heat source such as double ellipsoidal heat source as described above.  It is clearly seen from Fig. 9b that the effect of heat convection and radiation significantly reduced the peak temperatures generated at all three selected points A, B & C as the heat sinks levels were increased.    Figure 11a shows that the effect of heat radiation on the transient temperature distribution is more pronounced than the effect of the heat convection shown in Fig. 10a. The corresponding transient temperature reductions due to heat radiation effect are shown in Fig. 11b. The heat radiation effect caused the peak temperature reductions of ~143 o C, 117 o C and 82 o C at points A, B and C, respectively.  The maximum temperature reduction results showed that combined effect of both heat convection and radiation on the transient temperature distributions is significant. Therefore, these effects must be included in Gaussian distribution flux density heat source modelling to generate more accurate transient temperature field and metal molten pool geometry. Figure 12a shows the combined effect of heat convection and radiation on the transient temperature distributions of the three points A, B and C. The transient temperature generated are very similar to those generated for the effect of radiation of carbon steel of 0.67 emissivity at ambience temperature of 20 o C as shown in Fig. 11a. This is due to the fact that the effect of heat convection in windy air with h=100 W/m 2 / o C was small. The combined heat sink effects caused the peak temperature reduction of 143 o C, 117 o C and 82 o C at points A, B and C, respectively, as shown in Fig. 12b.

Simulate weld pool geometry due to the combined effects of heat sinks in the new complex density heat source
The combined effect of two heat sinks, convection (h = 100 W/m 2 / o C) and radiation (emissivity = 0.67) at 20 o C, also reflected in the weld pool geometric dimension reduction (i.e., in width, length and depth), compared to those generated by the analytical solution for the conduction-only Gaussian distribution heat source, as in Fig. 13. The width and length of the weld pool geometry was reduced by 1 and 2 mm, respectively, due to the combined effect of heat convection and radiation considered. However, the weld pool depth seems to remain unchanged as the weld pool depth is less likely to be affected by the effect of the two surface heat sinks considered. A FEA bead-on-plate model has been developed, simulating Nguyen's GMAW bead-on-plate experiment, to determine the parameters for the new complex heat source using available measured transient temperature data [5].
The GMAW process was simulated by using the Goldak's double ellipsoidal heat source model as was calibrated by using the measured temperature and weld pool geometry data by Nguyen et al. as measured from the same weld sample geometry and welding conditions [5]. The best-fitted heat source parameter matrix, determined for the numerical model, was [c hf c hr a h b h ] = [7 14 8.5 1.0]. These parameters were kept unchanged through-out all the FEA simulations to verify the validity of the new complex heat source. However, these heat source parameters were required to be adjusted for the corresponding approximate analytical solutions of this complex heat source, due to approximate nature of the analytical solutions.
The implementation of the new complex flux density heat source equations was carried out by using ABAQUS Fortranbased user-subroutine DFLUX codes, developed for body heat flux which is quite handy for Gaussian heat equations. The FORTRAN codes were written by the authors to simulate the two paths of the GMWA process on square bead-onplate sample, 240 x 240 x 20 mm, as shown in Fig. 2.
It has been shown that with suitable calibration coefficients (k rad ), both effects of convective and radiative heat sink films could be generated, either by implementing the complex flux density heat source model equations using DFLUX as body heat flux density only, or by using the conduction-only Goldak's heat source by using DFLUX in combination with two heat sink film FEA-based boundary conditions in the developed ABAQUS FEA model.
The transient temperature field generated from the calibrated Abaqus FEA model was in very good agreement with those calculated by the approximate analytical solution of the new complex density heat source, by using Matlab codes. That means that once the new complex heat source-sink was calibrated with the known FEA temperature results to determine the suitable complex heat source parameters, one could use the approximate analytical solutions of the complex density heat source to generate the transient temperature filed much more efficiently than by means of the transient FEA solvers. That would result in significant saving in CPU time for extensive and complex heat source modelling.
Lastly, the Matlab codes was used for generating weld pool geometry which was also agreed very well with the weld pool geometry data measured by Nguyen et al. [5] and by FEA-based simulated results. Thus, confirmed the credibility of the analytical solution of the new complex density heat source developed in this study.
However, the computing time of the Matlab code generated transient temperature field was 5 to 6 times faster than that of the FEA based modelling, using the same Intel® Dual Core™ i7-5700HQ CPU 2.70 GHz processor and 16GB RAM, i.e. it took ~9 to 10 min. per Matlab-based simulation, but ~53 to 60 min. per FEA-based simulation.
Thus, suggesting that the Matlab-based approximate analytical solution for the complex 3D heat source, developed in this work, could be an efficient alternative to the FEA-based heat source modelling for extensive simulation requirement. More details on the Matlab codes, developed in this research were briefly described in the Appendix A.

Comparison with the Vinokurov's heat convection film coefficient model
It has been reported in the literature that the effect of heat sinks, as a result of both heat convection and radiation, could be modelled as a sum of the heat transfer film coefficients due to convection and Vinokurov's effective to radiation loss model. The Vinokurov's model was tested by Guo Fig. 2) are compared with those obtained by FEA (conduction-only heat source with two heat sink films) and those generated by the new complex density heat source as shown in Fig. 14a. Figure 14a shows the Vinokurov's model consistently overestimated the temperature results at all three points (A, B and C) by 11 to 16% compared with those results generated by the Abaqus model using two heat sink films. This conclusion is consistent with conclusions from other authors on accuracy of the Vinokurov's model [7].
However, Fig. 14a shows very good agreement between the FEA-based transient temperatures and those generated by the new complex density heat source for the same condition. The deviation of the maximum temperatures for three points A, B and C were within 3%. Subsequently, an attempt in correcting the Vinokurov's temperature results was tried by a multiplication correction factor of 0.88. The corrected Vinokurov's model transient temperatures are plotted in Fig. 14b shows a good agreement with those generated by the FEA model with two heat sink films.
However, it is noticed that the transient temperature distribution of point C, generated by the corrected Vinokurov's model (with multiplication correction factor of 0.88), was still higher than those generated by the FEA-model with two heat sink films, particularly for the lower temperature regions. Based on extensive mathematical-analytical derivations and numerical simulations carried in this study, it is concluded that • The effect of heat convection and radiation could significantly reduce the peak temperature and weld pool size, therefore, must be included in all heat source simulation to improve the accuracy of the numerical modelling. • New 3D complex Gaussian density heat source model, developed in this study, has successfully included all three modes of heat transfer, i.e., heat conduction, convection and radiation. • The new complex Gaussian heat source has been successfully implemented by FEA-modelling by using usersubroutine DFLUX of Abaqus package. • It has been shown that the new complex Gaussian heat source was capable of generating the transient temperature field comparable with that generated by the finite element analysis, but with 5-6 times faster using the same PC configuration.  Potential future work would involve several case studies of the new complex density heat source and its approximate analytical solution applications to laser melting or 3D printing processes.

Acknowledgement
Contribution of Mohammad Nasiri in earlier stage of this project for initial Matlab codes implementing the heat convection effect, was greatly appreciated.
This project was purely volunteer research project, based on the tremendous efforts and endless time of the authors who have strong common interest and passion for advancing research in Gaussian heat sources, FEA-based simulation and its transient analytical solutions for potential industrial applications.
The second author gratefully acknowledges the financial support from the National Agency for Petroleum, Natural Gas and Biofuels (ANP)-PRH 18 and FINEP for his current fellowship, not directly related to this project as he worked on this project voluntarily in addition to his work's project.

The Declarations
a. Funding: No funding was received for this project.
b. Conflict of interest: There is no conflict of interest that the authors of this work aware of. h. Authors' contributions: The first author has derived the new complex Gaussian heat source by extending the Goldak's heat source to include all three modes of heat transfer (i.e., conduction, convection, and radiation); derivation of approximate analytical solutions and give direction to Abaqus-based finite element modelling; carried out Matlab implementation of the analytical solutions for transient temperature field for bead-on-plate and manuscript.
The second author has carried out the Abaqus based FEA modelling of the complex Gaussian heat source, heat source calibration, sensitivity study and manuscript.  Figure 15 shows a captured screen of the Matlab code run display at the end of the path 2 of the bead-on-plate GMAW simulation. The temperature contour surrounding the heat source is shown on the steel plate cross section, cutting through the centre planes of the weld paths.
There are four small windows designated for the real time display of the weld pool shapes, i.e., showing the real time sizes of the weld pool length, with and depth with constant iso-boundary temperature depicted at 1,370 o C.
The top two windows (Left & Right) showing the weld pool boundary due to heat conduction-only transient temperature solution for the Goldak's heat source by Nguyen et al. [5], using calibrated heat source parameters. The bottom two windows (Left & Right) showing the effective weld pool shape geometry (i.e., pool length, width and depth), determined from the approximate analytical transient solution for the complex heat source which incorporated all three modes of heat transfer, i.e., heat conduction, convection and radiation.
The approximate transient analytical solution of the complex heat source model was implemented by written Matlab codes, with calibrated Goldak's heat source parameter [c hf , c hb a h b h ] = [10 20 10 1.5], typical heat convection film in quiet air h = 5 W/m 2 / o C, grey body of steel plate with emissivity of 0.67 and radiative correction factor of k rad = 0.02 for temperature solution reduction, due to radiative heat sink as per the approximate analytical solution for the complex heat source equations, developed in this work.
The transient temperature depicted at three selected point A, B and C and its corresponding temperature reductions due to convective and radiative heat sinks are shown in Fig. 16 (a) and 16(b), respectively.