Nilanjan Chakraborty
email: n.chakraborty@liverpool.ac.uk

The effects of turbulence on molten pool transport during melting and solidification processes in continuous conduction mode laser welding of Copper-Nickel dissimilar couple

The effects of turbulence on molten pool transport during melting and solidification processes in continuous conduction mode laser welding of Copper-Nickel dissimilar couple

INTRODUCTION

Modern Laser welding applications often require the joining of two dissimilar metal parts, the physical properties of which are significantly different. The thermo-solutal history of the molten pool transport determines the weld-pool morphology, microstructure, and segregation pattern [START_REF] Sun | Laser welding of dissimilar metal combinations[END_REF].

Dissimilar metal welding has been extensively studied from the point of view of fusion characteristics and often these analyses are materials-specific [START_REF] Sun | Laser welding of dissimilar metal combinations[END_REF][START_REF] Metzger | Electron beam welding of dissimilar metals[END_REF]. Computational modelling of heat and mass transfer characteristics in the molten metal pool in welding applications has contributed significantly to the physical understanding of molten pool transport in similar metal welding applications [START_REF] Oreper | Heat and Fluid-Flow Phenomena in Weld Pools[END_REF][START_REF] Chan | Three Dimensional Axisymmetric Model for Convection in Laser Melted Pools[END_REF][START_REF] Kou | Three Dimensional Convection in Laser Melted Pools[END_REF][START_REF] Zacharia | Modelling of autogenous welding[END_REF][START_REF] Kanouff | The Unsteady Development of a GTA Weld Pool[END_REF][START_REF] Dutta | Determination of Gas Tungsten Arc Welding Efficiencies[END_REF][START_REF] Debroy | Physical processes in fusion welding[END_REF][START_REF] Dutta | Thermal Modelling of GTAW process with Non axisymmetric Boundary Conditions[END_REF][START_REF] Joshi | Non-axisymmetric Convection in Stationary Gas Tungsten Arc Weld Pools[END_REF][START_REF] Ehlen | Macroscopic modelling of Marangoni flow and Solute Redistribution during Laser Welding of Steel[END_REF]. Relatively few studies [START_REF] Chung | Mass, momentum and energy transport in molten pool when welding dissimilar metals[END_REF][START_REF] Phanikumar | Modelling of transport phenomena in laser welding of dissimilar metals[END_REF][START_REF] Phanikumar | Computational Modelling of Laser Welding of Cu-Ni Dissimilar couple[END_REF] addressed the modelling of dissimilar metal welding, despite its several industrial applications.

The scaling analysis and computational results [START_REF] Dutta | Thermal Modelling of GTAW process with Non axisymmetric Boundary Conditions[END_REF][START_REF] Joshi | Non-axisymmetric Convection in Stationary Gas Tungsten Arc Weld Pools[END_REF][START_REF] Chakraborty | Scaling analysis of momentum and heat transport in gas tungsten arc welding pools[END_REF] revealed that the maximum velocity obtained in the weld pool is in the order of 1.0m/s. For typical values of density and viscosity of . Atthey [START_REF] Atthey | A Mathematical Model for Fluid Flow in a Weld pool at high Currents[END_REF] indicated that in molten metal pools, laminar to turbulence transition sets in when a Reynolds number Re 2 Re = c exceeds 600, which is supported by subsequent studies [START_REF] Aboutalebi | Numerical Study of Coupled Turbulent Flow and Solidification for Steel Slab Clusters[END_REF][START_REF] Shyy | Modelling of Turbulent Transport and Solidification during Continuous Ingot Casting[END_REF] where surface tension driven convection (i.e.

Marangoni convection) is taken to be turbulent for Re in the order of 100. In addition to this Joshi et al. [START_REF] Joshi | Non-axisymmetric Convection in Stationary Gas Tungsten Arc Weld Pools[END_REF] reported perceptible surface oscillations in the molten pool for high power Gas Tungsten Arc Welding (GTAW) applications and indicated that such oscillations could be due to surface instabilities, which might eventually lead to turbulence. Moreover, the molten pool morphology obtained from experiments was found to be different from the prediction of numerical simulation using laminar fluid flow model [START_REF] Joshi | Non-axisymmetric Convection in Stationary Gas Tungsten Arc Weld Pools[END_REF]. It is therefore evident that the molten pool transport in weld pools is often likely to be turbulent for high power welding applications.

Mundra et al. [START_REF] Mundra | Role of Thermophysical Properties in Weld Pool Modeling[END_REF] attempted to mimic the turbulent momentum transport in similar metal welding by artificially increasing the molten metal viscosity. Choo and Szekely [START_REF] Choo | The Possible Role of Turbulence in GTA Weld Pool Behaviour[END_REF] and Hong et al. [START_REF] Hong | The Predicted Influence of Turbulence in Stationary Gas Tungsten Arc welds[END_REF] used two-equation models (i.e.

ε k

or w k -) to address molten pool convection in stationary GTAW applications. However, ε k or w k -equations are solved in a molten pool obtained from a conduction-based solution and were not solved in a coupled manner with other transport equations from the inception of molten pool. Moreover, the irregular morphology of the solid/liquid interface during welding would require wall treatment at phase boundaries. This demands treatment of solid-liquid interface with regard to the implementation of wall conditions for ε k equations. Some of the above issues were addressed in the computational methodology proposed by the present author and his co-workers, where modelled ε k equations were solved in a coupled manner with other transport equations to address molten pool transport in arc welding and laser surface alloying operations [START_REF] Chakraborty | Modelling of Turbulent Transport in Arc Welding Pools[END_REF][START_REF] Chakraborty | Three dimensional Modelling of Turbulent Weld Pool Convection in GTAW process[END_REF][START_REF] Chakraborty | Modelling of Turbulent Transport in Laser Surface Alloying[END_REF].

Recently, Chakraborty and Chakraborty [START_REF] Chakraborty | Modelling of Turbulent Molten Pool Convection in Laser Welding Copper-Nickel Dissimilar couple[END_REF] extended the methodology proposed by Chakraborty et al. [START_REF] Chakraborty | Three dimensional Modelling of Turbulent Weld Pool Convection in GTAW process[END_REF] to model turbulent transport in the molten pool during the melting stage in a typical dissimilar metal conduction mode laser welding. The thermo-solutal transport during the melting stage is likely to have significant influences on subsequent cooling characteristics and the species distribution in the cooled joint. In the present study, both heating and cooling cycles are simulated for a typical high power conduction mode laser welding of dissimilar metal couple where the effects of turbulent transport can be prominently realised. A binary metal couple based on Cu-Ni is taken for the study, as the physical properties of Copper and Nickel are significantly different and the phase diagram of Copper-Nickel does not show immiscible gaps and formation of intermetallic compounds [START_REF] Brandes | Smithells Metals Reference Book[END_REF].

In order to demonstrate the effects of turbulence in molten pool transport unsteady Reynolds Averaged Navier Stokes (RANS) simulations have been carried out for a continuous conduction mode laser welding of copper-nickel dissimilar couple. To demonstrate the contrasting features of turbulent convection, a simulation has been carried out corresponding to the process parameters used for the unsteady RANS simulation but without activating turbulence models. This simulation and the RANS simulation will henceforth be referred to as laminar and turbulent simulations respectively. The main objectives of this study are as follows:

(a) To model turbulent heat, mass and momentum transport in a copper-nickel dissimilar couple molten pool during both heating and cooling stages of the conduction mode laser welding process.

(b) To demonstrate the effects of turbulence on thermal history and species distribution of the welded joint.

The rest of the paper is organised as follows. In the next section, the necessary background for mathematical modelling will be provided along with the discussion on boundary conditions. The numerical implementation will be discussed in Section 3, followed by the presentation of the results and subsequent discussion. The main findings are summarised and conclusions are drawn in the final section of the paper.

MATHEMATICAL BACKGROUND

A schematic diagram of a typical continuous conduction mode laser welding of copper-nickel dissimilar couple is shown in Fig. 1a. Two blocks of copper and nickel are kept side by side in a butt joint arrangement. According to Fig. 1a the laser torch moves along the interface separating copper and nickel blocks. Only a small part of laser power acts to increase temperature of the metal blocks, which is specified by efficiency η following Phanikumar et al. [14,15] as:

Q q η = (1)
For the present study it is assumed that the laser torch is moving in the negative x direction with a scanning speed scan U . It is convenient to study the present problem in a moving co-ordinate system, which is carried out by the following co-ordinate transformation:

t U x x scan ) (- - ′ = ; y y ′ = and z z ′ = (2) where ) , , ( z y x ′ ′ ′
is the stationary co-ordinate system and ) , , ( z y x is the moving co-ordinate system translating with the laser torch. In Eq. ( 2) the negative sign in front of scan U indicates that the laser torch is moving in the negative x direction.

In order to model the melting and intermixing of two different metals, a locally homogeneous mixture model has been employed, where the mixture properties are determined by the relative proportion of the components. The following assumptions are invoked in the present mathematical model, which simplifies the analysis without losing the essential physics of the problem:

1. The welding is assumed to take place in conduction mode and no keyhole formation is considered.

2. The molten metal flows are assumed to incompressible and Newtonian in nature and thus the density change due to the temperature variation is accounted for by Boussinesq's approximation.

3. The efficiency of absorption of heat from laser beam η is assumed to be same for both Cu and Ni.

4. Laser power is assumed to be distributed in a Gaussian manner at the top surface.

5. The top surface of the welded joint is assumed to statistically flat [START_REF] Dutta | Thermal Modelling of GTAW process with Non axisymmetric Boundary Conditions[END_REF][START_REF] Phanikumar | Modelling of transport phenomena in laser welding of dissimilar metals[END_REF][START_REF] Phanikumar | Computational Modelling of Laser Welding of Cu-Ni Dissimilar couple[END_REF].

The generalised form of the governing equations for turbulent molten pool transport in tensorial notation (

z x y x x x = = = 3 2 1 , ,
) can be written as:

( ) ( ) φ φ ρ φ φ ρ φ ρ S x u f x x f u x f t j j j j j j + ∂ ′ ′ ∂ - ∂ ∂ Γ ∂ ∂ = ∂ ∂ + ∂ ∂ ) . . ( (3)
where φ represents the mean value and φ ′ represents the fluctuating component of a primitive variable φ , Γ is the diffusion coefficient associated with φ , and φ S is the source/sink term associated with the particular governing equation in question. The corresponding values ofφ , , f Γ and S φ , for various conservation equations, are given in Table 1.

For the momentum conservation equations, the reference temperature T ref is taken to be equal to the melting temperature m T , and C ref is the reference concentration, taken to be equal to the concentration of copper at that reference temperature. The source term i S in the momentum conservation equations originates from the consideration that the morphology of the phase change domain can be treated as an equivalent porous medium offering a frictional resistance towards fluid flow in that region. In a single-domain fixed-grid enthalpy-porosity formulation, this resistance can be conveniently formulated using Darcy' s model in association with the Cozeny-Karman relationship as [START_REF] Brent | Enthalpy-Porosity Technique for Modelling Convection-Diffusion Phase Change: Application to the Melting of a Pure Metal[END_REF]:

( ) i l l m i u b f f K S + - - = 3 2 1 (4)
where l f is the liquid fraction, which is given as:

L H f l / ∆ =
, with H is the latent enthalpy content of a control volume and L is the latent heat of fusion. The melting temperature m T and the latent heat of fusion L are evaluated as follows according to locally homogeneous mixing model [START_REF] Phanikumar | Modelling of transport phenomena in laser welding of dissimilar metals[END_REF][START_REF] Phanikumar | Computational Modelling of Laser Welding of Cu-Ni Dissimilar couple[END_REF][START_REF] Chakraborty | Modelling of Turbulent Molten Pool Convection in Laser Welding Copper-Nickel Dissimilar couple[END_REF]:

) 1 ( ) ( ) ( C T C T T Ni melt Cu melt m - + = and ) 1 ( ) ( ) ( C H C H L Ni Cu - + = (5)
where C is the mass fraction of copper in a given control volume, melt T

is the melting point temperature, H is the latent heat of fusion of pure material, and the subscripts Cu and Ni are used to refer to the properties of pure copper and pure nickel respectively. The latent enthalpy content of the computational cell H under consideration given by the following expression

L T f H = = ∆ ) ( : l T T > L f l = : l s T T T < ≤ (6) 0 = : s T T <
with T s and T l being the solidus and liquidus temperatures respectively. In Eq. ( 4), m K is a large number (~10 8 ) and b is a small number (~10 -30 ) to avoid division by zero. It is important to note that the momentum conservation equations are solved for the whole domain. In the solid phase the liquid fraction l f becomes identically zero which leads to a large magnitude for the sink term i S , which numerically ensures that fluid velocity is zero in the solid phase. In the fully molten state l f is equal to unity and as a result of this the contribution of i S becomes identically zero in the liquid phase (see Eq. ( 4)) and the momentum conservation equations take the usual form for pure liquid phase. The details of this formulation may be found in Brent et al. [START_REF] Brent | Enthalpy-Porosity Technique for Modelling Convection-Diffusion Phase Change: Application to the Melting of a Pure Metal[END_REF]. This technique has been used in several solid-liquid phase change related problems in the past [10, 11, 14-16, 18, 23-26].

Modelling of Reynolds stress terms

) ( j i u u ′ ′ -ρ
In the present analysis, the Reynolds stress terms

j i u u ′ ′ -ρ
appearing in momentum conservation are modelled by a gradient hypothesis of the form:

k x u x u u u ij i j j i t j i ρ δ µ ρ 3 2 ) ( - ∂ ∂ + ∂ ∂ = ′ ′ - (7)
where the eddy viscosity t µ is given by:

ε ρ µ µ / 2 k C f l t = (8)
where k and ε are turbulent kinetic energy and the dissipation rate of turbulent kinetic energy respectively. It is important to note that in a standard high Reynolds number According to Eq. ( 8), the eddy viscosity goes to zero in the solid phase (since liquid fraction is identically zero in solid), and assumes the standard expression for single-phase flow problems in a fully liquid phase (since l f is identically unity in liquid phase).

Modelling of Turbulent heat and species fluxes

The turbulent heat and species fluxes (Reynolds heat and mass fluxes) appearing in the energy and species conservation equations respectively are given by the following expressions:

i t i x T T u ∂ ∂ = ′ ′ - α and i c i x C C u ∂ ∂ = ′ ′ - α (9)
where t α is the eddy thermal diffusivity. From the analogy of laminar flow, t α can be expressed as: In the present case t σ and c σ are taken to be 0.9.

t t t ρσ µ α =

Transport equations of k and

The modelled governing equations for k and are given in Table 1. The model constants are given by:

Boundary conditions

In this study, the entire work-piece is taken to be the computational domain. The heat flux due to laser beam q ′ ′ is assumed to follow a Gaussian profile in the following manner

) / exp( ) / ( ) ( 2 2 2 q q r r r Q r q - = ′ ′ π η
. All the above conservation equations are solved in a coupled manner, in conjunction with a set of appropriate boundary conditions. Boundary conditions pertinent to various governing equations are summarized in Table 2.

k-transport at near wall regions

For k and transport, the phase changing interface is treated as a solid wall. For a high Reynolds number k-model, eddy diffusivities are specified at grid points next to the solid-liquid interface with the help of wall functions [START_REF] Durbin | Statistical Theory and Modelling for Turbulent Flows[END_REF]. Mean velocities gradient, turbulent kinetic energy and turbulent kinetic energy dissipation rate at the grid points next to the wall are described by:

+ - + = + + + κ σ µ σ µ ρ (12ii)
where Pr and Sc represent the molecular Prandtl and Schmidt numbers respectively. Values of k and are identically zero in the solid phase.

NUMERICAL IMPLEMENTATION

In this study, the simulations were carried out based on finite volume approach using the wellknown SIMPLER algorithm [START_REF] Patankar | Numerical Heat Transfer and Fluid Flow[END_REF]. The convective terms were discretised using the power law scheme proposed by Patankar [START_REF] Patankar | Numerical Heat Transfer and Fluid Flow[END_REF]. The time advancement has been carried out in an implicit manner. The phase change aspect is accounted for by the enthalpy-porosity technique in which the latent heat enthalpy content of every computational cell is updated using the following expression [START_REF] Brent | Enthalpy-Porosity Technique for Modelling Convection-Diffusion Phase Change: Application to the Melting of a Pure Metal[END_REF][START_REF] Voller | The Modeling of Heat, Mass and Solute Transport in Solidification Systems[END_REF]:

{ } melt n P P P m m T T c a a H H - + ∆ = ∆ + ] [ ) ( ] [ ] [ 0 1 λ (13)
where P a is the nodal point coefficient and 0 P a is the coefficient related with the transient part of the discretised energy conservation equation.

The physical properties used in the present study for both Ni and Cu are listed in Table 3 [START_REF] Brandes | Smithells Metals Reference Book[END_REF].

For a mixture of Ni and Cu, properties such as specific heat c, thermal conductivity K and viscosity are calculated with the use of the following expression:

) 1 ( C P C P P Ni Cu - + = (14)
where P represents the property concerned and subscripts Cu and Ni represent copper and nickel respectively. The mixture properties are calculated at each time step after solving the species conservation equation. The results presented in this paper correspond to process parameters typical of dissimilar metal laser welding (U scan = 0.008 m/s, Q =3.5 kW, η = 0.12, r q = 0.5 mm).

Grid resolution

The grid resolution in the present study is determined by two criteria which are as follows:

• Resolution of boundary layer formed beneath the top surface.

• Resolution of the flow physics in the vicinity of the phase changing interface acting as a wall to the fluid motion inside the molten pool.

In order to meet the first requirement, the boundary layer thickness v δ is estimated based on scaling arguments following Dutta et al. [START_REF] Dutta | Thermal Modelling of GTAW process with Non axisymmetric Boundary Conditions[END_REF] as:

v s U δ µ ~R T T sur ∆ ∂ ∂σ (15)
where U s represents the surface velocity scale, T/R represents the thermal gradient scale, R stands for maximum pool width and represents the temperature difference between the area just beneath the laser source and the phase changing boundary. About five grid points are kept within the boundary layer thickness.

In the context of high Reynolds number

ε - k model, the quantity µ ρ µ / 2 / 1 4 / 1 x k C x ∆ = + has to be greater than 40 (i.e. 40 > + x
) where x ∆ is the grid spacing. However, in order to have an accurate estimation of wall heat flux, the grid spacing has to be small. This is taken care of by ensuring that +

x remains within the following range 50 20 ≤ < + x , which ensures that the near wall points fall either within the buffer layer or within the logarithmic layer.

From the above indicated criteria it becomes clear that the grid size spacing along the x 2 axis is governed by criterion 1, and the grid size spacing along the 1

x and 3

x directions are principally governed by criterion 2. For the topmost five grids a spacing of 2.42 x 10 -6 m is chosen, 4.8 x 10 -6

for the three grids immediately beneath, followed by a grid spacing of 7.2 x 10 -6 m and subsequently a grid spacing of 9.6 x 10 -6 m over six consecutive grids after which grid spacings of 1.2 x 10 -5 m and 1.45 x 10 -5 respectively follow in the negative 2

x direction. For the remaining part of the pool, a uniform grid spacing of 2.42 x 10 -5 m is employed. In regions outside the molten pool, a non uniform less dense grid is chosen. For 1

x and 3

x directions near the solidliquid interface, the grid spacing assigned is to be 6.43 x 10 -5 m after which the grid spacing gradually increases to 1.34 x 10 -4 m away from the molten zone. In order to discretise the domain of

mm mm mm 8 5 8
× × a 69x48x69 Cartesian grid is used. The above choice of grid size is an optimal compromise between numerical accuracy and computational economy. The present study assumes that the laser torch is located at x 1 =4mm, x 2 =5mm and x 3 =4mm with respect to the moving co-ordinate system. The computational domain in the region

Time advancement

For the present problem parameters, the initiation of melting takes about 0.6 ms. Before melting, heat transfer takes place in conduction mode, where a relatively larger time step (about 0.2 ms) is chosen so that melting temperature is reached within three time steps. However, with the inception of melting, the high temperature gradient in the pool sets up a strong surface tensiondriven flow, leading to a convection-dominated transport. Since, over this time span, convection in the molten pool develops rapidly, a time step size only as small as about 0.05 ms can lead towards monotonic convergence of the solution during this period of initial transience. Typically after about 5 ms, the molten pool reaches a state when changes in dependent variables between the consecutive time steps are relatively small. At this stage, slightly larger time steps (typically about 0.07 ms) can safely be used to reduce computational time. After 10 ms, the pool becomes sufficiently developed, which allows time steps as high as 0.0001s without any oscillation. After 25 ms the laser power is switched off and the molten pool is allowed to cool down for 7.0ms as done previously by Phanikumar et al. [START_REF] Phanikumar | Computational Modelling of Laser Welding of Cu-Ni Dissimilar couple[END_REF]. By this time, the molten pool resolidifies and the final composition pattern in the weld pool is obtained.

Convergence criteria

The convergence for all the independent variable ( ) , , , , , ,

3 2 1 ε k C T u u u
values are checked after each iteration within a given time step. The convergence is assessed at each grid point based on the method of relative error given by the following criterion:

4 max 10 - ≤ - φ φ φ old (16)
The thermo-solutal dependence of surface tension and simultaneous convergence of all the primitive variables make the simulation procedure computationally intensive. During the pure conduction stage the convergence was achieved in about 50 iterations with a given time step and each time step took about 5 minutes CPU time on a single dual core Intel Xeon processor with 8GB RAM and Intel 5000X motherboard. For both the laminar and turbulent simulations each time step took about 500 iterations to converge during the inception of surface tension driven convection and an each time step typically took 45 minutes to complete. Once the pool is reasonably well developed it takes about 100 and 300 iterations to converge for each time step which required about 9 minutes and 30 minutes of CPU time for laminar and turbulent simulations respectively.

RESULTS & DISCUSSION

The Prandtl numbers of molten copper and nickel are about 0.008 and 0.04 respectively. This indicates that, for a given Reynolds number, the Peclet number assumes a much smaller value in the copper molten pool than in the nickel molten pool. This suggests that thermal conduction primarily governs energy transport in the copper molten pool, whereas this effect is likely to be weaker in nickel. For the present values of process parameters, Prandtl number and pool dimensions ( q r d / , where d is the penetration of the molten pool), the molten pool convection in the case of nickel falls within the turbulent flow regime (where turbulence has significant effects on momentum transport and has marginal influence on thermal transport), but also lies near the boundary with the fully turbulent regimes (where turbulence significantly affects both momentum and thermal transport) in the regime diagram as shown in Fig. 1b. By contrast, the molten pool convection in the case of copper falls in the boundary of conduction dominated (where thermal conduction plays the major role in molten pool transport and convection has marginal role in determining molten pool morphology) and turbulent flow regimes (see Fig. 1b). In Fig. 1b, the non-dimensional parameter Ch is given by [START_REF] Chakraborty | Thermal transport regimes and Effects of Prandtl number in molten pool transport in Laser Surface Melting processes[END_REF]:

1 ∂ ∂ = K Q T Ch sur π η µ σ ρ µ (17i)
The non-dimensional parameter Ch is closely related to Marangoni number Ma in the following manner:

µ ρ q s q r U Ma Ch r d = ) / ( 2 3 / 2 (17ii)
where Marangoni number can be taken as a surface-tension driven Reynolds number. This is consistent with earlier observations [START_REF] Phanikumar | Modelling of transport phenomena in laser welding of dissimilar metals[END_REF][START_REF] Phanikumar | Computational Modelling of Laser Welding of Cu-Ni Dissimilar couple[END_REF][START_REF] Chakraborty | Modelling of Turbulent Molten Pool Convection in Laser Welding Copper-Nickel Dissimilar couple[END_REF]. The temperature fields obtained from laminar and turbulent transport models are qualitatively similar to each other, but the maximum mean temperature obtained from the turbulent simulation is found to be smaller than the maximum temperature attained in the laminar simulation. The enhanced level of thermal diffusion in the turbulent pool acts to reduce the maximum mean temperature in the turbulent transport model (e.g. are the maximum temperatures obtained from laminar and turbulent simulations respectively), which is found to be in agreement with previous studies [START_REF] Choo | The Possible Role of Turbulence in GTA Weld Pool Behaviour[END_REF][START_REF] Chakraborty | Three dimensional Modelling of Turbulent Weld Pool Convection in GTAW process[END_REF][START_REF] Chakraborty | Modelling of Turbulent Transport in Laser Surface Alloying[END_REF][START_REF] Chakraborty | Modelling of Turbulent Molten Pool Convection in Laser Welding Copper-Nickel Dissimilar couple[END_REF].

The temperature field obtained from the melting stage in the welding process provides quantitative information about the solidification parameters at the solidification front. The temperature gradient to solidification rate under quasi-steady condition is given by ) cos /(

β scan U G
where G is the characteristic temperature gradient, β is the angle between the welding direction and normal on the solid-liquid interface on the solidification [START_REF] Flemmings | Solidification Processing[END_REF] ) is not satisfied in this particular case and the solidification front is likely to be non-planar. Some of the above nonplanar features of solidification in Cu-Ni dissimilar were observed in the experimental data of Phanikumar et al. [START_REF] Phanikumar | Continuous welding of Cu-Ni dissimilar couple[END_REF].

The molecular diffusion driven thermal transport can be substantiated by examining the ratio of the turbulent thermal diffusivity to the molecular thermal diffusivity in the molten pool after ms 25 of laser aided heating as shown in Figs. 3a-c. Figure 3 in the region where surface tension differential due to temperature is in opposition to surface tension differential due to composition. As a result of this, in the present case,

1 1 / / / x C x T ∂ ∂ ∂ ∂ and 3 3 / / / x C x T ∂ ∂ ∂ ∂
remain close to the limiting value (in the order of 1000), and the relative strength of thermo-solutal variation determines the surface tension gradient in the top surface.

From the velocity vector plots in Figs. 4a and4b it is evident that, for both laminar and turbulent simulations, surface tension differential due to temperature prevails over the surface tension differential due to composition and the velocity is directed towards the pool boundary from the maximum temperature location. However, it is worth noting from the velocity vector plots (see Figs. 4a and4b) that the vectors are not exactly radially divergent from the hot temperature location, which should have been the case if the convection were determined solely by the temperature induced surface tension differential. In case of turbulent transport, the mean thermal and concentration gradient decreases due to the increase in diffusion strength. This eventually decreases the strength of mean convection in the turbulent pool. Moreover, the maximum value of velocity magnitude in the turbulent pool decreases in comparison to that obtained in case of laminar transport (e.g. the maximum surface velocities obtained from laminar and turbulent simulations are 3.16m/s and 1.90m/s respectively), caused by the greater degree of momentum transport due to turbulent dispersion. As the velocity vectors are directed towards the pool boundary from the maximum temperature location, the heat is transported from the maximum temperature location to pool boundaries by means of advection, in addition to thermal diffusion.

In case of turbulent transport, the mean advection is weaker compared to that obtained in the case of laminar transport. This tends to decrease the width of the turbulent pool in comparison to the pool obtained in laminar simulation, which is countered by the increased thermal diffusion strength in the case of turbulent transport. The relative strength of these two mechanisms ultimately determines the pool width in the turbulent pool. It can be seen from comparing Figs. 5a and 5b that the pool width in the turbulent case is found to be marginally smaller than that obtained from the corresponding laminar case.

The fluid flow from the hottest part of the pool to the pool boundary (i.e. solid-liquid interface) transports heat in the longitudinal and sidewise directions. The fluid stream turns in the downward direction as the solid-liquid boundary is approached, as governed by the continuity relation. This gives rise to an asymmetrical toroidal circulation pattern in the molten pool, as evident from the velocity vector plots in longitudinal and cross-sectional mid-planes as shown in 3). This suggests that turbulence in the molten pool significantly affects the mass transport in the dissimilar metal weld pool. The eddy thermal diffusivity t α is equal to the eddy mass diffusivity c α in the present analysis because of the same values of turbulent Prandtl and Schmidt numbers (see Eqs. ( 8) and ( 9)). As t α remains of the same order of molecular thermal diffusivity in the molten pool (see Figs. 3a-c) the ratio of eddy mass diffusivity to molecular diffusivity (i.e. D D t / ) remains in the order of Pr / Sc (i.e. ). This suggests that turbulent diffusion will have significant effects on species distribution in the molten pool, which can be substantiated by comparing the Cu mass fraction distribution obtained from laminar simulation with the Cu concentration field obtained from the corresponding turbulent simulation.

It is apparent from Figs. 5a and 5b that the distribution of Cu concentration in the turbulent pool is significantly different from the corresponding laminar pool. In both laminar and turbulent cases fresh copper is added at the melting front in the copper rich side (the pool boundary in the direction of laser torch movement), which locally increases the copper concentration in the molten pool. The added copper further diffuses towards the centre and then carried towards the Ni side of the molten pool by advective transport owing to fluid flow towards the boundary. In the Ni side of the pool the copper concentration is diluted by the fresh addition of molten Ni inside the molten pool. On the solidification front (pool boundary opposite to the laser torch movement) preferential rejection of solute (Ni here) in the Cu rich side of the molten pool decreases the copper concentration. The preferential rejection of Cu takes place in the solidification front of the Ni rich side, which locally increases the Cu concentration. These effects are evident in both laminar and turbulent pools. However, in the case of laminar pool the iso-concentration contours are severely deformed as the species transport in the laminar pool is primarily advection-driven because of the negligible value of mass diffusivity. By contrast, in case of the turbulent pool, the species distribution pattern shows clear stratification of Cu concentration, which is a characteristic of diffusion process with decreasing Cu concentration from the Cu to Ni side of the molten pool. In the turbulent pool, the mean advection strength becomes weaker than that in the corresponding laminar pool due to greater amount of momentum diffusion. Moreover, mass diffusion strength increases in the turbulent pool because of turbulent mixing; as a result of which eddy mass diffusion plays a significant role in species transport in turbulent molten pools. This leads to the so-called "banded" structure in the Cu concentration field in the turbulent molten pool [START_REF] Phanikumar | Modelling of transport phenomena in laser welding of dissimilar metals[END_REF][START_REF] Phanikumar | Computational Modelling of Laser Welding of Cu-Ni Dissimilar couple[END_REF]. The fresh nickel added at the melting front on the Ni side ultimately diffuses towards the middle of the molten pool which causes the iso-concentration contours of Cu mass fraction to be concave towards the fresh Ni side of the molten pool. The extent of this behaviour is more prominent in the turbulent pool because of enhanced diffusion rate of Ni. The isoconcentration contours of Cu close to Ni side of the molten pool appear as micro-structural banding in the post-weld experimental analysis of the weld pool section [START_REF] Phanikumar | Computational Modelling of Laser Welding of Cu-Ni Dissimilar couple[END_REF]. The sparsely packed iso-concentration contours on the Ni side and densely packed contours on the Cu side are also in agreement with earlier experimental observations [START_REF] Phanikumar | Computational Modelling of Laser Welding of Cu-Ni Dissimilar couple[END_REF]. Moreover, the asymmetrical pool shape in Cu and Ni side of the pool is also consistent with earlier findings [START_REF] Chung | Mass, momentum and energy transport in molten pool when welding dissimilar metals[END_REF][START_REF] Phanikumar | Modelling of transport phenomena in laser welding of dissimilar metals[END_REF][START_REF] Phanikumar | Computational Modelling of Laser Welding of Cu-Ni Dissimilar couple[END_REF].

The laser power is switched off after s m 25 and the workpiece is allowed to cool down for s m 7 .

By this time, the workpiece is subjected to convective and radiative cooling with the same values of heat transfer coefficients used for different surfaces of the domain during the melting stage.

The temperature contours after s m 0 . 7 of cooling of molten pools obtained from laminar and turbulent transport models are shown in Fig. 6. Figure 6 clearly indicates that the maximum temperature is smaller than the melting temperature ( melt T T <

) for both cases. This suggests that the weld pool resolidifies within s m 7 , and that the temperature fields obtained after resolidification based on laminar melting simulation are qualitatively similar to those obtained based on turbulent melting simulations. The copper mass fraction fields after re-solidification, resulting from both laminar and turbulent transport simulations are shown in Fig. 7. Comparing Figs. 5 and7 it can be seen that the composition of the molten pool after re-solidification remains similar to that obtained after the melting stage. The differences between Cu mass fraction fields obtained using laminar and turbulent transport models discussed in the context of Fig. 5 are therefore also valid for Fig. 7. This essentially indicates that the final composition pattern in the Cu-Ni dissimilar weld pool is essentially determined at the melting stage, and the mass transport during the solidification stage does not have a significant effect on the solute distribution in the molten pool.

The macrographs of the top surface and cross-sectional view of a Cu-Ni welded joint for the same laser power and scanning speed obtained from experimental data [START_REF] Phanikumar | Computational Modelling of Laser Welding of Cu-Ni Dissimilar couple[END_REF] are shown in Figs. 8a andb respectively. The comparison between the macrographs presented in Figs. 8a andb and the computational Cu mass fraction variations presented in Figs. 5 and7 suggests good qualitative agreement between the computational and experimental results. In order to assess the quantitative agreement between experimental and computational data, the variations of Ni mass fraction distribution with the normalised distance along the pool width for both laminar and turbulent simulations are compared with the experimental results reported by Phanikumar et al. [START_REF] Phanikumar | Computational Modelling of Laser Welding of Cu-Ni Dissimilar couple[END_REF] in Fig. 8c for the same laser power and scanning speed. In the case of laminar transport, molten copper is advected to a comparatively larger distance in the molten pool than in the turbulent transport simulation, which decreases the concentration of Ni in comparison to that obtained using turbulent transport, as the major part of the molten pool is created in the Ni-rich side of the welded joint. This effect is much weaker in the turbulent pool due to enhanced mixing, and the decrease in mean advective transport. This is consistent with greater (smaller) mean Cu (Ni) concentration in laminar pool in comparison to turbulent pool as depicted in Fig. 8c. It should also be noted that the concentration profile in laminar simulation is found to be steeper than that obtained in the case of turbulent simulation. This is due to the smaller mean concentration gradient in turbulent pool than the corresponding laminar pool owing to the more efficient mixing in turbulent transport. It is clearly evident that although both laminar and turbulent simulations capture the qualitative variation of the Ni mass fraction, the concentration profile obtained from turbulent simulation provides more accurate agreement with the experimental profile. Comparing Ni mass fraction distributions obtained from melting phase and after complete solidification, it is evident that the concentration pattern does not get altered significantly in the cooling and solidification stages for both laminar and turbulent transport models. This behaviour can be explained in the following manner. With decreasing temperature, the temperature gradient at the top surface starts to decrease which also weakens the strength of surface-tension driven advection in the molten metal pool. Moreover, the velocity gradient in the molten pool is directly proportional to the temperature gradient in Marangoni convection (see Eq. 18), which also affects the turbulent kinetic energy generation rate in the case of turbulent transport. This suggests that the generation of turbulent kinetic energy decreases with decreasing temperature, which results in a decay in the value of turbulent diffusivity and the flow ultimately re-laminarises before the solidification process. Decreasing advection strength in both laminar and turbulent transport reduces the extent of species transport in the molten pool. In the case of turbulent transport, the decay of turbulent diffusivity decreases the extent of diffusive transport, which further limits the extent of mass transport in the molten pool. As a result of this, the solute distribution in the weld pool is principally determined by the advective and diffusive transport during the melting stage and it does not get altered significantly during the solidification stage of dissimilar metal welding.

CONCLUSIONS

Melting and solidification stages of a conduction mode laser dissimilar welding of Cu-Ni couple have been numerically simulated using a finite-volume based approach. The solid-liquid phase change is accounted for by a modified enthalpy porosity technique [START_REF] Brent | Enthalpy-Porosity Technique for Modelling Convection-Diffusion Phase Change: Application to the Melting of a Pure Metal[END_REF]. The momentum, heat and mass transports in the molten pool are studied using both laminar and turbulent transport models for identical values of laser power and scanning speed. In the laminar transport model, the governing equations of mass, momentum and energy are solved simultaneously in three dimensions, whereas unsteady Reynolds averaged Navier Stokes simulations are carried out for the turbulent simulation, in which turbulent transport is accounted for a high Reynolds number ε k model. It has been found that temperature and velocity fields obtained from laminar and turbulent simulations are qualitatively similar to each other. However, the maximum values of mean temperature and velocity magnitude obtained from turbulent simulation are found to be smaller than the corresponding values obtained from the laminar simulation due to the enhanced level of diffusive transport in the case of turbulent simulation. The species distribution gets significantly affected by the enhanced diffusive strength in the turbulent molten pool and the species distribution obtained from turbulent transport shows a clear structure of concentration stratification. The stratified structure of species distribution, as obtained from the turbulent simulation, where iso-concentration contours are almost parallel to the interface between Cu and Ni workpieces is found to be consistent with earlier experimental results for the same process parameters [START_REF] Phanikumar | Computational Modelling of Laser Welding of Cu-Ni Dissimilar couple[END_REF]. The species distribution in the final welded joint is found to be determined principally during the melting stage, and the species distribution does not get altered significantly during cooling and solidification stages of the process because of the weakening of convection strength. This effect is further aided by the decrease in eddy diffusion strength with decreasing temperature in turbulent pools.

TABLES Table 1: Diffusion coefficients and source terms for various conservation equations

Governing equations

φ f Γ φ S Continuity 1 1 0 0 Momentum in i th direction i u 1 µ - ) ( ) ( 1 scan 2 2 x u U S C C g T T g x P i i ref C i ref T i i ∂ ∂ + - + - + ∂ ∂ - ρ β δ ρ β δ ρ Energy T c K ( ) 1 ) ( ) ( x T c U H u x t H scan i i ∂ ∂ - ∆ ∂ ∂ - ∂ ∆ ∂ - ρ ρ ρ Species (Cu mass fraction) C 1 D ρ 1 x C U scan ∂ ∂ -ρ Turbulent kinetic energy k 1 k t σ µ µ + 1 2 2 x k U x C g x T g x u x u x u scan C c t T t t j i i j j i t ∂ ∂ - - ∂ ∂ - ∂ ∂ - ∂ ∂ ∂ ∂ + ∂ ∂ ρ ρε β σ µ β σ µ µ Dissipation rate of turbulent kinetic energy ε 1 ε σ µ µ t + 1 2 2 2 1 2 1 1 x U k C k x C g C k x T g C k x u x u x u C scan c c t T k t j i i j j i t ∂ ∂ - - ∂ ∂ - ∂ ∂ - ∂ ∂ ∂ ∂ + ∂ ∂ ε ρ ε ρ ε β σ µ ε β σ µ ε µ ε ε ε ε
0 2 = u i sur i sur i i x C C x T T x u ∂ ∂ ∂ ∂ + ∂ ∂ ∂ ∂ = ∂ ∂ = σ σ µ τ 2 2
(where = i 1 and 3)

0 1 = u 0 2 = u 0 3 = u 0 1 = u 0 2 = u 0 3 = u Not required (follows from enthalpy-porosity formulation) Energy conservation equation ) ( ) ( ) ( 4 4 2 α α σ ε T T T T h r q x T K rad r - + - + ′ ′ - = ∂ ∂ - 0 2 = ∂ ∂ x T ) ( / ∝ - = ∂ ∂ - T T h n T K wall Not required Species conservation equation 0 / 2 = ∂ ∂ x C 0 / 2 = ∂ ∂ x C 0 C n ∂ ∂ = Freezing front : C k v n C D P n ) 1 ( - = ∂ ∂ - Fusion front : C v n C D n = ∂ ∂ - Turbulent kinetic energy equation 0 2 = ∂ ∂ x k k = 0 k = 0 k = 0 Dissipation rate of turbulent kinetic energy 0 2 = ∂ ∂ x ε = 0 = 0 0 = ∂ ∂ n ε
The value of ε at near wall points are specified through wall condition

Table 3: List of physical properties

Physical properties

Values

β β β β T (Compressibility) 0.45×10 -4 K -1 β β β β C (Compressibility) 11.0×10 -2 T m (Melting point) 1085 C 0 (Cu) 1453 C 0 (Ni) L (Latent heat) 1.88×10 5 J/kg (Cu) 2.90×10 5 J/kg (Ni) ρ ρ ρ ρ (

f

  does not appear in the eddy viscosity expression. In the present problem, one encounters both solid and liquid phases in the computational domain. As a result, the formulation should ensure a reduction of eddy viscosity, eddy thermal conductivity and eddy mass diffusion coefficient (all of which is referred to as eddy diffusivity) in such a manner that the effective diffusivity (molecular + eddy) values approach their respective molecular values along the solidliquid interface, and gradually become equal to the single-phase turbulent flow conditions in the fully liquid region. In order to satisfy this requirement, Shyy et al.[START_REF] Shyy | Modelling of Turbulent Transport and Solidification during Continuous Ingot Casting[END_REF] used a modified eddy viscosity expression in the context of a low Reynolds number ε k model where the damping factor l f was used to address the solid-liquid phase change problem. The present work is based on a high Reynolds number ε k model, and accordingly, the wall-damping factor used by Shyy et al. [19] in the context of a low Reynolds number ε k model is not employed.

σ

  are the turbulent Prandtl number and turbulent Schmidt number respectively.

C

  are standard model constants for high-Reynolds number ε k model and these model constants were evaluated by calibrating the computational simulation results with experimental data under decaying homogeneous isotropic turbulence and homogeneous shear flow configurations. For more information on these model constants, interested readers are directed to Durbin and Pettersson Reif [29, see pp.120-122].

u

  represents the mean velocity component tangential to the solidliquid interface at the near wall grid point. The quantity P y represents the normal distance of the grid-point next to the interface, w represents the wall shear stress, and is the Von Karman's constant which will be taken as =0.4. The thermal and solutal diffusivities at near wall points are given as[START_REF] Kader | Temperature and concentration profiles in fully turbulent boundary layer[END_REF]:[ ]

mm x 4 3 ≤ 3 ≥ 2 plane at mm x 4 3 =

 3323 is taken to be the Cu work piece and the rest (i.e.mm x 4) is considered to be the Ni work piece. The x 1 -x corresponds to the initial contact plane between Cu and Ni work pieces.

Figure

  Figure1bsuggests that thermal transport in dissimilar welding of the copper-nickel couple is primarily driven by molecular diffusion (thermal conduction) in the molten pool. As a result, turbulence in the molten pool significantly affects the momentum transport, but the temperature field distribution only gets marginally affected.

10 -

 10 Figs. 5a and 5b it is evident that ) 10 ( ~2 O T ∆ , whereas ) 10 ( ~1 -∆ O C

Figs

  Figs. 4c-f for both laminar and turbulent simulations.

1

  -1.3 ×10 -4 (T -T m ) N/m (Cu) 1.778-3.4 ×10 -4 (T -T m ) N/m (Ni) Based on private communication with Prof. H. K. Bhadeshia, Cambridge University

Fig. 1 :Fig. 2 :.

 12 Fig. 1: (a) Schematic diagram of continuous laser welding of a Cu-Ni dissimilar couple; (b) The location of molten pool convection for copper and nickel on the regime diagram for thermal transport proposed by Chakraborty [36] for the present processing parameters.Fig. 2: Isotherms after 25ms of laser heating at the top surface: (a) laminar simulation, (b) turbulent simulation; at the longitudinal mid-plane: (c) laminar simulation, (d) turbulent simulation; at the cross-sectional mid-plane: (e) laminar simulation, (f) turbulent simulation. All the temperatures are shown in C 0 . In moving co-ordinate system the co-ordinate of the laser torch is x = 4mm, and z = 4mm. At time t = 0s z= 4.0mm indicates the separation line between Cu and Ni workpieces.

Fig. 3 :

 3 Fig. 3: The ratio of eddy thermal diffusivity to molecular thermal diffusivity α α / t

Fig. 4 :

 4 Fig. 4: Velocity vectors after 25ms of laser heating at the top surface: (a) laminar simulation, (b) turbulent simulation; longitudinal mid-plane: (c) laminar simulation, (d) turbulent simulation; cross-sectional mid-plane: (e) laminar simulation, (f) turbulent simulation. The isotherms corresponding to

Fig. 5 :

 5 Fig. 5: Contours of Cu mass fraction after 25ms of laser heating at the top surface: (a) laminar simulation, (b) turbulent simulation; at the longitudinal mid-plane: (c) laminar simulation, (d) turbulent simulation; at the cross-sectional mid-plane: (e) laminar simulation, (f) turbulent simulation. In moving co-ordinate system the co-ordinate of the laser torch is x = 4mm, and z = 4mm. At time t = 0s z= 4mm indicates the separation line between Cu and Ni workpieces.

Fig. 6 :Fig. 7 :

 67 Fig. 6: Isotherms after 7ms of cooling at the top surface: (a) laminar simulation, (b) turbulent simulation; at the longitudinal mid-plane: (c) laminar simulation, (d) turbulent simulation; at the cross-sectional mid-plane: (e) laminar simulation, (f) turbulent simulation. All the temperatures are shown in C 0 .In moving co-ordinate system the co-ordinate of the laser torch is x =

Fig. 8 :

 8 Fig. 8: Experimental macrographs of the Cu-Ni welded joint for KW Q 5 . 3 = and

Fig. 1 :Fig. 2 .

 12 Fig. 1: (a) Schematic diagram of continuous laser welding of a Cu-Ni dissimilar couple; (b) The location of molten pool convection for copper and nickel on the regime diagram for thermal transport proposed by Chakraborty [36] for the present processing parameters.

Fig. 4 :Fig. 5 :

 45 Fig. 4: Velocity vectors after 25ms of laser heating at the top surface: (a) laminar simulation, (b) turbulent simulation; longitudinal mid-plane: (c) laminar simulation, (d) turbulent simulation; crosssectional mid-plane: (e) laminar simulation, (f) turbulent simulation. The isotherms corresponding to C 0 1000 , C 0 1500 , C 0 2000 and C 0 2200 are shown in Figs. 4a and b by dotted lines. The temperature value increases from circumference to the centre of the plot. In Figs. 4c-f isotherms corresponding to C 0 1200 , C 0 1500 , C 0 2000 and C 0 2200 are shown by dotted lines. In moving co-ordinate system laser torch located at x = 4mm, and z = 4mm. At time t = 0s z= 4mm indicates the separation line between Cu and Ni workpieces.

Fig. 6 :.Fig. 7 :

 67 Fig. 7: Contours of Cu mass fraction after 7ms of cooling at the top surface: (a) laminar simulation, (b) turbulent simulation; at the longitudinal mid-plane: (c) laminar simulation, (d) turbulent simulation; at the cross-sectional mid-plane: (e) laminar simulation, (f) turbulent simulation. In moving co-ordinate system the co-ordinate of the laser torch is x = 4mm, and z = 4mm. At time t = 0s z= 4.0mm indicates the separation line between Cu and Ni workpieces.

Fig. 8 :

 8 Fig. 8: Experimental macrographs of the Cu-Ni welded joint for KW Q 5 . 3 = and

Table 2 : Table of boundary conditions Boundary conditions Governing Equation Top surface Bottom surface Side faces Solid/liquid interface

 2

	Mass and
	momentum
	conservation
	equations

4mm, and z = 4mm. At time t = 0s z= 4.0mm indicates the separation line between Cu and Ni workpieces.

ACKNOWLEDGEMENTS

The author is grateful to Drs. S. Chakraborty and G. Phanikumar for many useful discussions.

The author also acknowledges the help of Dr. M. Home while preparing the manuscript. Figures 8a andb are taken from Figs. 9a and 12a of Ref. [15] with the kind permission of Springer Science and Business media.


Miguel Torres García

Francisco José

Jiménez-Espadafor Aguilar

; Tomás

Sánchez Lencero

José Antonio Becerra

A new Heat release rate (HRR) law for Homogeneous Charge Compression Ignition (HCCI) combustion mode

Homogeneous charge compression ignition (HCCI) engines are drawing attracting attention as the next-generation's internal combustion engine, mainly because of its very low NOx and soot emissions and also for improvement in engine efficiency. Much research has been carried out in order to go deeper in this combustion process using

multizone models or CFD codes. These simulation tools, although they can give a detailed view of the combustion process, are very time consuming and the results depend a lot on the initial conditions. A previous step to be considered in the simulation of the HCCI process is a heat release law evaluated from results of the experiment and a zero dimensional model. This paper focuses on the development of a new heat release rate (HRR) law that models the HCCI process when the combustion chamber is considered as a homogeneous volume. The parameters of this law have been adjusted through an optimization process that has allowed to fit the combustion chamber pressure. All the engine operative conditions from low to full load have been successfully simulated with this HRR law, with the maximum error in the estimation of combustion chamber pressure less than 2%.

Keywords: HRR (heat release rate); EGR (exhaust gas recirculation); HCCI (homogeneous charge compression ignition), in-cylinder NO x . Woschnic [START_REF] Woschni | Equation for the Instantaneous Heat Transfer Coefficient in the Internal Combustion Engine[END_REF] or Hohenberg [4], which have been elaborated in engine conditions that vary significantly from those in HCCI combustion.

Nomenclature

The dominant heat transfer mechanism in HCCI combustion mode is forced convection from hot gases to the combustion chamber walls. Due to very low soot formation and relatively low bulk gas temperature, the radiation effect is very small, the opposite being the case in the conventional diesel combustion mode. For this reason, any HRR analysis for HCCI combustion mode should be developed through a dedicated heat transfer model. In this paper, that of Junseok Chang has been used [5].

Regarding the rate of heat release in HCCI combustion mode, it is not controlled by the rate of fuel injection as in DI engines, nor by finite turbulent flame propagation as in SI engines. The absence of an ignition control mechanism has led researchers to explore a range of control strategies. Performing these explorations solely in the laboratory would be inefficient, expensive, and impractical since there are many variables that exhibit complex interaction. Fundamental tools, such as computational fluid dynamics codes (CFD) with detailed chemistry, need to be applied in order to provide insight into the combustion process. Codes of this nature are very computationally intensive and usually require some simplifications to expedite the solution while attempting to maintain accuracy. Fully coupled CFD/kinetic models have employed a number of different combustion methodologies. Kusaka [6] had no explicit modification of the reaction rate due to mixing; the only mixing that occurred during combustion was that normally calculated by the CFD code. Results with a detailed kinetic scheme showed a more rapid combustion event than the experiment. Because of the relatively coarse grid that was used, there is the possibility that the mass in colder zones was underestimated. Kong et al. (2001) [START_REF] Kong | Use of Detailed Chemical Kinetics to Study HCCI Engine Combustion with Consideration of Turbulent Mixing Effects[END_REF] coupled the KIVA-3V CFD code with kinetics modified for turbulence effects and found a similar effect; i.e., when no turbulence effect was used, they predicted a higher energy release rate than was seen in the experiment. All these models present a severe drawback: the huge computational load that precludes efficient and quick analysis of the HCCI process. In this paper, a theory-experimental model is developed which allows evaluating the HRR in HCCI combustion mode in a modified diesel engine (DEUTZ-DITER FL1 906) from an analytical law of HRR. In all experimental studies carried out on HCCI combustion mode, it is observed that the HRR is quite different from the one in diesel combustion mode. The HRR of the HCCI process presents only a premixed combustion shape while that of the diesel mode shows premixed and diffusive appearance combustion. The knowledge of an analytical HRR law validated in a wide range of operating conditions in HCCI combustion mode will be a great advance in knowledge of HCCI combustion and could allow developing a predictive computational model for this combustion process.

This paper develops a new heat release rate (HRR) law that models the HCCI process.

The parameters of this law have been adjusted through an optimization process that has allowed to fit the combustion chamber pressure. The paper is organised as follows:

 Engine experimental setup, where the test rig and experimental procedure is described.



Model for the evaluation of the HRR. Here the model developed for HRR from experimental measurement is explained.



Analytical law for HRR in HCCI combustion mode, where the new HRR is presented.

 Numerical solution procedure.

 Results and conclusions.

ENGINE AND EXPERIMENTAL SETUP

The experimental part of this work is based on the tests carried out on the modified DEUTZ FL1 906 engine. The original diesel engine was modified to adapt it to HCCI combustion. The main characteristics of the engine are shown in Table 1 and the experimental installation is represented in Figure 1. The modified engine systems were:

1º Injection System: the injection pump was extracted from the engine block, and its design modified to control the injection point. The gear of the crankshaft coupling to the injection pump was substituted by a transmission belt.

2º Exhaust gas recirculation system (EGR): an EGR system was designed which included an exhaust gases refrigeration system. With this system, the EGR by-passed mass fraction and the temperature of the recirculated gases can be controlled.

3º Measurement equipment: different sensors were located in the engine to monitor the pressure in the combustion chamber, the injection pressure, needle rising, engine speed, exhaust temperature, fuel mass flow, the intake system, pressure evolution, temperature and air mass flow.

Start of combustion is highly dependent on intake temperature. To carry out a full analysis of the intake temperature effect in different operating conditions, a heating system for the intake air was designed (figure 1). Also, a heat exchanger cooler for the EGR gases was installed to control the recirculated gases temperature and so to control the temperature of the cylinder intake mixture. The installation incorporates stagnation boxes to eliminate the strong pulsating waves of the mono-cylinder engine.

Table 1 Tests were done with a commercial diesel fuel compatible with European Standard EN590. The recirculated gases (EGR) were cooled in order to control the intake temperature that never surpassed the ambient temperature. The external EGR rate is evaluated measuring the EGR mass recirculated and the total mass inlet. The formula is as follows:

100 ) (%    Air EGR EGR m m m mass EGR    (1)
where EGR m  and Air m  are measured with a hot wire flow meter

MODEL FOR THE EVALUATION OF THE HRR

In this work, a thermodynamic zero-dimensional model has been used to carry out the analysis and evaluation of the parameters of an analytical law of HRR, following an iterative optimization process. The combustion chamber is considered as a perfect mixture reactor with variable volume, with even pressure distribution, temperature and concentration of chemical species (see Figure 2) with heat losses. Equations and simplifications that govern the mathematical model are needed due to complex interaction between physics and chemical phenomenon during combustion and are described in the following equations [8]:  Heat losses  Blow by has not been considered.

The heat loss model is written as follows:

    1 W g c c W T T h A d dQ   (2)
The heat transfer coefficient h c (W/m 2 K) is based on the correlation formula given by Junseok Chang [5]. W Q is the wall heat loss, A c is the area in contact with the gases, T g (K) is the gas temperature, T w (K) is the cylinder wall temperature and ω (rad/s) is the average angular speed. The global heat transfer coefficient can be written as:

8 . 0 73 . 0 8 . 0 2 . 0 ) ( ) ( ) ( ) ( ) ( t v t T t p t L t h scaling c     (3)
A scaling factor scaling  is used for tuning of the coefficient to match specific engine geometry. Combustion-induced gas velocity is a function of the difference between motoring and firing pressure [START_REF] Xingcai | The effects of external exhaust Gas recirculation and Cetane number improver on the gasoline homogeneous charge compression ignition engines[END_REF].

  motoring r r r d p p p V p T V C S C t v    6 ) ( 2 1 (4)
The main idea for using this equation is to keep the velocity constant during the unfiring period of the cycle, and to then impose a steep velocity rise once combustion pressure departs from motoring pressure. The subscript r denotes a reference crank angle, such as intake valve closing.

Figure 3 shows the net cumulative heat release obtained from a fuel consumption of 0.015 g/cycle, 1800 RPM and 0% EGR. Combustion efficiency comes from measured exhaust gas composition, 93-96%, so the energy released by the fuel per cycle (LHV=42.5 MJ/kg) is between 593-612J. Crevice losses are considered very low [5], about 2%. Parameter scaling  from equation (3) has been tuned so that the cumulative heat released by fuel plus heat losses, plus 2% from crevices, is equal to the energy released by the fuel. This analysis was repeated over the whole range of operating conditions tested.

Figure 3 In the experimental results there is no kinetic information, therefore the start of combustion is defined from the Heat Release Rate (HRR) obtained from the registered in-cylinder pressure combined with a zero-dimensional model of the combustion system: the combustion starts when the cumulative energy released is 5% of the total [START_REF] Xingcai | The effects of external exhaust Gas recirculation and Cetane number improver on the gasoline homogeneous charge compression ignition engines[END_REF][START_REF] Huber | Investigations on Heat Transfer in Internal Combustion Engines under Low load and motoring conditions[END_REF]. The combustion duration has been established as the angle between the start of combustion and that corresponding to 10% of the magnitude of the peak of HRR on the falling side of the curve [START_REF] Xingcai | Heat Release Analysis on Combustion and Parametric Study on Emissions of HCCI Engines Fuelled with 2-Propanol/n-Heptane Blend Fuels[END_REF].

HRR ANALYSIS

In general, cold EGR plays a significant role on HCCI combustion characteristics. As shown in Figure 4, the maximum in-cylinder pressure decreases, the ignition delay increases, the peak value of HRR diminishes and the combustion process is prolonged with the introduction of cold external EGR. The increment of ignition delay diminishes the amount of fuel burned before TDC, which reduces the maximum combustion pressure. Also, the temperature of the gases is lower, which decelerates the combustion process. Although combustion is prolonged with any increment of the EGR rate, which makes that combustion process separate from constant volume combustion, the specific fuel consumption improves at the end. In Figure 4, any increment of EGR produces an increment of engine torque due to a reduction of the fuel burned and, therefore, of the pressure BTDC up to a threshold level that the torque reduces. This is associated with lack of oxygen, that is, more than about 20% of EGR to 2100 RPM reduces the inlet air and so there is not enough oxygen to burn the fuel. With an increase of cold EGR rate, the mixture temperature rising rate during the compression process decreases because of the increase of heat capacity. Also, the oxidation and decomposition reaction rate were depressed for the final products such as CO 2 , NO x , and H 2 O.

ANALITICAL LAW FOR HRR IN HCCI COMBUSTION MODE

The HRR curves shown in Figure 4 reveal the very high heat release rate that causes a rapid increase in pressure in the engine cylinder. From this, combustion is characterized by the sudden steep pressure rise on the cylinder pressure curve and the peak on the HRR. Once all the already formed flammable mixture is burned, the HRR decreases until the end of combustion. For all fuel consumption rates, engine speeds and EGR rates tested in HCCI mode, the HRR has shown the same development.

Therefore, combustion is completely controlled by chemical kinetics. Because gas cylinder temperature is not sufficiently high, the fuel does not reach pyrolysis conditions almost avoiding soot formation and diminishing NO x formation drastically [START_REF] Peng | Characteristic of Homogeneous charge Compression Ignition (HCCI) combustion and Emissions of n-Heptane[END_REF].

A model for the HRR that can be adapted efficiently to a combustion process controlled only by chemical kinetics is that of Wiebe. In this work, it has been found that modifying the Wiebe function has allowed the formulation of a new HRR law that fits quite well with the HCCI combustion process [START_REF] Yasar | Double-Wiebe function: An approach for single-zone HCCI engine modeling[END_REF].

                           1 2 1 exp ) 1 ( pp p M p M p p p p a M Q a d dQ       ( 5
)
Where a 1 , a 2 , M p , M pp are shape factors. In order to reproduce the HRR the parameter M p has had to be considered different from the parameter M pp , see (5). M p has been considered as a different parameter from M pp (the opposite of the original Wiebe function), θ p is the duration of the energy release, and Q p characterizes the heat release in HCCI combustion mode.

Defining the terms:

  1 1 1 1    p M p p M Q a K p  1 2 2 1   pp M p a K  (6)
It follows that the HRR is:

  1 2 1 exp    pp p M M K K d dQ    (7)

NUMERICAL SOLUTION PROCEDURE

Once the HRR law has been chosen, the problem is the identification of values of the parameters of this law, K 1 , K 2 , M p , M pp for each operative engine condition.

Finding the values of the sought parameters is an optimization problem that can be formulated as:

 Minimize L          n i n i pp p i simulated i erimental i M M K K P P L 1 1 2 2 1 exp 2 ) , , , , ( ) (    ( 8
)
Where i  is a measure of the error, P experimental and P simulated are the measured and simulated combustion chamber pressure respectively, and θ i is every crank angle where pressure is measured and simulated.  Heat losses [5] The equations of the thermodynamic model presented are solved numerically step-bystep using a simple time marching technique, and they give the combustion chamber pressure, P simulated .

The kinetic chemical model is based on a replacement fuel (PRF) fuel to provide the beginning of combustion. The start of combustion, or ignition point, is controlled by physical processes, such as the atomization of the fuel, fuel evaporation and fuel -air mix, and by the chemical kinetic reactions among the fuel, air and combustion products [START_REF] Torres Garcia | Analysis of the Start of Combustion of a Diesel Fuel in a HCCI Process through an Integral Chemical Kinetic Model and Experimentation[END_REF]. The direct simulation of the oxidation reactions of commercial diesel fuel is unavailable due to its composition, formed by thousands of chemical species, and the uncertainty of some mechanisms of oxidation. So commercial fuel is usually modelled on a PRF of a well-known composition and kinetic mechanisms [START_REF] Richter | The Influence of charge Inhomogeneity on the HCCI Combustion Process[END_REF].

OPTIMIZATION PROCESS

The scheme of the numerical solution procedure is represented in Figure 6. The theoryexperimental model is a direct type, the start is the pressure chamber curve and the end is the HRR curve. The time of computation can be raised because, during the optimization process, a direct model has to be resolved in every iteration and the rate of convergence depends on the initial conditions used and the accuracy. One of the most decisive factors in computation time is the initial condition. With the new HRR law, the optimized parameters have allowed good concordance between experimental and combustion pressure chamber simulations. Figure 7 compares simulated cylinder pressures, obtained as discussed in previous sections, with measured ones, derived from the experimental study. The modelled cylinder pressure shown has been obtained from the new HRR law for different mass flow and EGR; the angular velocity is 1800 RPM and the fuel consumption is constant (1.3e-5 kg/cycle).

In terms of the area under the pressure curve between IVC and EVO, the error is about

1%. An interesting aspect to highlight is the good reproduction of pressure during the combustion process, that is, when an abrupt increase of pressure takes place.

Table 2 shows the four parameters of the new HRR law for quite different operating engine conditions.

Table 2 Figures 8 and9 show the measured and predicted cylinder pressures and the heat release rate obtained from the new HRR law. The results show good concordance between experimental and numerical results; both the angular velocity and fuel consumption are constant. As can be observed, the new HRR law adapts perfectly to any load condition.

From the minimum EGR rate (5%), each increment of the EGR rate has three effects on the HRR a) A delay in start of combustion.

b) A diminution of the maximum heat release rate.

c) An increase of combustion duration.

From Table 3 and for engine conditions in Figure 9 it can be observed that every EGR increment also increases engine torque. The increase of torque is due to the diminution of combustion pressure during the compression stroke because of the delay in start of combustion, effect a. Because of this power increase at the same engine speed and fuel consumption, there is also an improvement in the engine's specific fuel consumption to weight of effect c that reduces thermodynamic efficiency [START_REF] Nishijima | Premixed Lean Diesel Combustion (PREDIC) using Impingement Spray System[END_REF]. Effect b produces a diminution of the combustion chamber temperature and so there will be a diminution of NO x emissions [START_REF] Murata | Emissions suppression mechanism of premixed diesel combustion with variable valve timing[END_REF].

Figure 8 Figure 9 Table 3 quantifies the error between measurement and simulation built from the new HRR law. The difference between maximum pressure measured and predicted is less than 0.8%, the error of the area under cylinder pressure is less than 2% but is reduced to 1.4% when cylinder pressure is between inlet valve opening (IVO) and exhaust valve closing (EVC).

Figure 2 

 2 Figure 2

Figure 5

 5 Figure 5 shows the combustion ignition angle and combustion duration versus EGR rate

Figure 5 According

 5 Figure 5

5. 1 

 1 EVALUATION OF THE COMBUSTION CHAMBER PRESSURE FROM HRR LAW. Combustion pressure, P simulated, is evaluated through a zero zone model which includes the following submodels: Thermal state equation  Combustion chamber volume  Heat release rate law  First principle of thermodynamic  Start of combustion model[START_REF] Torres Garcia | Analysis of the Start of Combustion of a Diesel Fuel in a HCCI Process through an Integral Chemical Kinetic Model and Experimentation[END_REF]

Figure 6 6.

 6 Figure 6

3 .

 3 Accuracy of numeric model.

Figure 1 .

 1 Figure 1. Experimental installation.

Figure 2 .Figure 3

 23 Figure 2. Control volume for zero-dimensional analysis

Figure 4

 4 Figure 4 Combustion chamber pressure versus crank angle with a compression ratio of

Figure 5

 5 Figure 5 Combustion ignition angle and combustion duration as a function of EGR rate

Figure 6 .

 6 Figure 6. Scheme of the optimization procedure.

Figure 7 .

 7 Figure 7. Measured and predicted cylinder pressure traces with a 1.3e-5kg/cycle of fuel

Figure 8 .

 8 Figure 8. Measured and predicted cylinder pressure and HRR simulated in HCCI combustion

Figure 9 .

 9 Figure 9. Measured and predicted cylinder pressure and HRR simulated in HCCI combustion

Table 3

 3 As has been shown, the HRR law depends on four parameters that can be related to any load condition. Research is in progress on the development of a predictive model of the engine in HCCI combustion mode. 9. Measured and predicted cylinder pressure and HRR simulated in HCCI combustion mode (angular velocity is 1500 RPM, constant fuel consumption 1.8e-5 kg/cycle, initial Ф=0.3 and different EGR rate).

	7. CONCLUSIONS
	In this work, an experimental and simulation study has been carried out to compare the
	performance of a new HRR law that defines a proportion of slower combustion for
	HCCI engine modelling. The new HRR law was implemented in an engine model to
	evaluate performance in comparison to the experimental data obtained in detailed tests.
	The study showed that by describing a proportion of slower combustion with the new
	HRR proposed, it was possible to achieve a very good match to experimental data. The
	new HRR law allows predicting the cylinder pressure curve perfectly with minimum
	error.

2. The Best-fit parameters from new HRR law for different tests.

Table 1 .

 1 Deutz FL1 906 Engine characteristics.

	Type	DI Monocylinder 4 stroke
	Cylinder displacement	708 cm 3
	Bore	95 mm
	Stroke	100 mm
	Compression ratio	19:1

Table 2 .

 2 The Best-fit parameters from new HRR law for different tests.

	RPM	Fuel consumption (10 -5 kg/cycle)	EGR % (mass)	K 1	K 2	M p	M pp
	1350	1.24	35	54.7	0.28	0.76	2.65
	1650	2.45	30	67.4	0.17	1.38	1.92
	1800	2.98	28	98.4	2.89	2.55	1.15
	2000	2.52	13	84.	0.26	1.98	1.61
	2400	3.02	21	43.4	0.13	3.01	1.48
	1200	1.74	38	75.3	0.037	0.51	1.40

Table 3 .

 3 Accuracy of numeric model.

								Error
	RPM	Torque (Nm)	m f (10 -5 kg/cycle)	% EGR (mass)	Error % Pressure maximum	Error % presusure area (total)*	Error % pressure area (parcial)**	ε rms (curve cylinder
								pressure)
	1200	5	1.63	2	0.142	0.15	0.35	0.046
	1200	10	1.63	47	-0.24	0.413	0.41	0.038
	1200	10	2.19	2	0.03	-1.12	-1.25	0.047
	1200	12	2.19	13	-0.11	-1.36	-1.65	0.274

* Between IVC until EVO. ** from 35º BTDC until 35º ATDC.

ACKNOWLEDMENTS

The work in this paper is a part of item CTQ2007-68026-CO2-02/PPQ within the I+D+i national plan in the period 2007-2009 and has been supported by the Spanish Government (Ministry of Science and Education). The authors are grateful to the Ministry of Science and Education of Spain for their financial support of this work.