J J Del Coz Díaz

P Menéndez Rodríguez

P J García Nieto

D Castro-Fresno

Comparative analysis of TIG welding distortions between Austenitic and Duplex Stainless Steels by FEM

Keywords: Welding process, distortion, modelling, FEM, temperature-dependent properties

In this study, thermal stress analyses were performed in the tungsten inert gas (TIG) welding process of two different stainless steel specimens in order to compare their distortion mode and magnitude. The growing presence of non-conventional stainless steel species like duplex family generates uncertainty about how their material properties could be affected under the welding process. To develop suitable welding numerical models, authors must consider the welding process parameters, geometrical constraints, material nonlinearities and all physical phenomena involved in welding, both thermal and structural. In this sense, four different premises are taken into account. Firstly, all finite elements corresponding to the deposition welding are deactivated and, next, they are reactivated according to the torch's movement to simulate mass addition from the filler metal into the weld pool. Secondly, the movement of the TIG torch was modelled in a discontinuous way assuming a constant welding speed. Thirdly, the arc heat input was applied to the weld zone using volumetric heat flux distribution functions. Fourthly, the evolution of the structural response has been tackled through a stepwise non-linear coupled analysis. The numerical simulations are validated by means of full-scale experimental welding tests on stainless steel plates. Finally, the results and conclusions of this research work are exposed.

Nomenclature

Q & thermal heat source rate, W/m 2 s L q surface heat loss, W/m 2 c L q convection surface heat loss, W/m 2 r L q radiation surface heat loss, W/m 2 ij s deviatoric stress tensor 1 3 Duplex stainless steels are finding increased application in the chemical, oil and gas industries, petrochemical process plants, the pulp and paper industry, pollution control equipment, transportation and for general engineering thanks to their outstanding corrosion resistance and mechanical properties [START_REF] Azuma | Duplex stainless steel excellent in corrosion resistance[END_REF].

ij kk ij σ σ δ     = -        , Pa t time, s U x nodal displacement in X direction, m U y nodal displacement in Y direction, m U z nodal displacement in Z direction
The high corrosion resistance and the excellent mechanical properties combination of duplex stainless steels can be explained by their chemical composition and balanced "duplex" microstructure of approximately equivalent volume fractions of ferrite and austenite. Firstly, the chemical composition based on high contents of Cr and Mo, improves intergranular and pitting corrosion resistance, respectively. Moreover, additions of nitrogen can promote structural hardening by an interstitial solid solution mechanism, which raises the yield strength and ultimate strength values without impairing toughness. Secondly, the two-phase microstructure guarantees higher resistance to pitting and stress corrosion cracking in comparison with conventional stainless steels [START_REF] Álvarez-Armas | Duplex Stainless Steel[END_REF]. Duplex stainless steels are typically twice as strong as common austenitic stainless steels. The thermal expansion in the duplex grades is intermediate to that of carbon steel and the austenitic stainless steels. The thermal conductivity in duplex stainless steels is also intermediate to that of carbon steels and the austenitic stainless steels. Duplex stainless steels, with higher thermal conductivity and lower coefficient of thermal expansion, do not create the same high intensity of local thermal stresses at the welds of austenitic stainless steels [START_REF] Cunn | DDS Microstructure, properties and applications[END_REF].

Aim

A complex state of thermal and residual stresses is developed in a welded structure as a direct consequence of the uniform heat supplied and the subsequent cooling process, these stresses being able to reduce the load-carrying capacity of the structure. The high level of traction stresses in the neighborhood of the seam can increase the tendency to brittle fracture and can originate serious changes in stress corrosion and fatigue behavior. The residual compression stresses in connection with the associated distortions can originate local or global instability. All of this justifies designer interest in knowledge of the residual stress state, both qualitatively and quantitatively. This knowledge would also permit a more rational and precise use of the stress-relief techniques, especially those based on a compensation of the existing residual stresses [START_REF] Cañas | A simplified numerical analysis of Residual stresses in aluminum welded plates[END_REF]. It is generally understood that welding involves highly complex phenomena, arising from interaction between thermal transients and the elastoplastic and structural responses of the assembly, together with the metallurgical transformation of the material being welded and the effects of turbulence on molten pool transport [START_REF] Camilleri | Alternative numerical techniques for distortion of thin plate due to fillet-welded stiffeners[END_REF][START_REF] Chakraborty | The effects of turbulence on molten pool transport during melting and solidification processes in continuous conduction mode laser welding of copper-nickel dissimilar couple[END_REF]. Following this last research trend, the main goal of this paper is to establish a set of approximations and simplifications in the numerical modeling of stainless steels welding, in order to reduce the huge computation time due to the coupled and non-linear formulations. Moreover, after reaching the above mentioned target, it is possible to tackle the different welding distortions for duplex and austenitic materials and compare them. Finally, the numerical results from the finite element model are compared with the experimental values measured in the specimens and conclusions are exposed.

Stainless Steel material properties and welding process parameters

The material properties of weld metal, base metal and heat affected zone are all temperature dependent. Due to the lack of information on material properties, both thermal and mechanical properties of weld metal and heat affected zone were assumed to be the same than those of the base metal in this analysis.

Fig. 1. Thermo-physical material properties for duplex and austenitic stainless steels: specific heat (left) and thermal conductivity (right). On the one hand, at about 1,773 K the duplex steel is melted and over this value there is liquid. The LDX-2101 (EN-1.4462 according to British Standard rule [START_REF]Design of steel structures. Part 1-2: General rules-Structural fire design[END_REF]) thermal material properties are reported in Fig. 1. On the other hand, the austenitic stainless steel EN-1.4404 [START_REF]Design of steel structures. Part 1-2: General rules-Structural fire design[END_REF] is melted about 1,573 K. From the above database it is not possible to know the material properties in the liquid phase, and the unique way at this stage is to keep the physical properties constant outside the solid range. The temperaturedependent mechanical properties for duplex and austenitic stainless steels used in the simulation are shown in Fig. 2. Fig. 2. Mechanical material properties for duplex and austenitic stainless steels: (a) yield strength, (b) ultimate tensile strength, and bilinear piece wise tensile-elongation relationship for different temperatures: (c) duplex steel and (d) austenitic steel. The amount of heat input for both stainless steels were found as the product of arc efficiency, voltage, current and welding speed, as it is shown in Table 1. The total net heat input was assigned to the weld metal as uniform body loads. Table 1 Welding parameters

Mathematical model

In order to build accurate numerical models of welding, we are going to consider all the parameters of the process (number and sequence of steps, welding speed, etc.), material nonlinearities, geometrical constraints, and all physical phenomena involved in welding following the previous work by Capriccioli and Frosi [8]. In this sense, it is a great daring to take into account all the important factors of the welding process at the same time. Therefore, this research work assumes the following main hypotheses and characteristics about the thermal model [8][START_REF] Zhu | Effects of temperature-dependent material properties on welding simulation[END_REF][START_REF] Zhu | Numerical simulation of transient temperature and residual stresses in friction stir welding of 304L stainless steel[END_REF]:

• Radiation and convection effects have been considered here.

• All the material properties are described till the liquid metal phase.

• During the welding process, the displacement of the parts does not affect the thermal distribution of the parts themselves. • The birth and death technique has been used [8] to simulate the material deposition, as it is explained in subsection 3.3. The weldment thermo-mechanical behaviour during the welding process is simulated using uncoupled formulation, because the dimensional changes in welding are negligible and the mechanical work done is insignificant compared to the thermal energy from the welding arc. Therefore, an uncoupled thermal and thermo-mechanical analysis is adopted in this work [START_REF] Zhu | Effects of temperature-dependent material properties on welding simulation[END_REF][START_REF] Zhu | Numerical simulation of transient temperature and residual stresses in friction stir welding of 304L stainless steel[END_REF][START_REF] Murthy | Numerical simulation of welding and quenching processes using transient thermal and thermo-elasto-plastic formulations[END_REF]. In the first place, a thermal analysis was carried out in the first step and during this phase the temperatures distributions are calculated and saved for every load step. In the second place, the transient temperature outputs from the previous analysis are used for a thermo-mechanical subsequent analysis. It is assumed that the thermal calculation at a given time is independent from the structural results obtained at a previous time according to the first hypothesis. Thus, the thermal and mechanical analyses can be uncoupled [8][START_REF] Zhu | Effects of temperature-dependent material properties on welding simulation[END_REF][START_REF] Zhu | Numerical simulation of transient temperature and residual stresses in friction stir welding of 304L stainless steel[END_REF][START_REF] Murthy | Numerical simulation of welding and quenching processes using transient thermal and thermo-elasto-plastic formulations[END_REF]. The heat conduction problem is solved independently from the stress problem to obtain the temperature history. However, the formulation considers the contributions of the transient temperature field to the stress analysis through thermal expansion, as well as temperature-dependent thermo-physical and mechanical properties [8]. The solution procedure consists of two steps. First, the temperature distribution and its history in the welding model are computed by the heat conduction analysis. Then, the temperature history is employed as a thermal load in the subsequent mechanical elastoplastic calculation of the residual stress field.

Thermal analysis

In the thermal analysis, the transient temperature field θ in the molten weld pool is a function of time t and the spatial coordinates ( )

, , x y z and is assumed to be governed by the same equation as the applied to the solid metal region. The governing partial differential equation for transient heat conduction, with internal heat generation, is given by the thermal equilibrium equation; that is to say, it is determined by the following three-dimensional nonlinear heat transfer equation [START_REF] Murthy | Numerical simulation of welding and quenching processes using transient thermal and thermo-elasto-plastic formulations[END_REF][START_REF] Bathe | Finite Element Procedures[END_REF][START_REF] Cook | Concepts and Applications of Finite Element Analysis[END_REF][START_REF] Chandrupatla | Introduction to Finite Element in Engineering[END_REF][START_REF] Mills | Heat Transfer[END_REF][START_REF] Del Coz Diaz | Nonlinear thermal optimization of external light concrete multi-holed brick walls by the finite element method[END_REF]:

( ) ( ) ( ) 2 2 2 2 2 2 Q C x y z t θ θ θ θ λ θ ρ θ θ   ∂ ∂ ∂ ∂ + + + =   ∂ ∂ ∂ ∂   & (1)
It is important to take into account that the heat flux to the system is input by a moving source on the boundary. Therefore, the materials are melted over 1,773 K and 1,573 K (see section two/2). In order to simulate the phase change we have used the enthalpy formulation ( ( ) ( )

H C d ρ θ θ θ ∆ = ∫
). This function depends on the previous values of density ρ and specific heat C at the same temperature [8]. The integration arbitrary constant is assumed equal to the product C ρ θ at 273 K (0ºC). In this way, at the solidus temperature, it is necessary to add the latent heat to obtain the enthalpy value at the liquidus temperature. Therefore, it is possible to model with only one material parameter the whole thermal history including heat exchange during solidification or melting. In the liquid phase, the parameters C and ρ are assumed constant, so that the enthalpy in that interval has a constant slope given by the last values defined for C and ρ .

Surface heat flow, that can take place into or from the model during the processes cited above, is specified by the following boundary conditions [START_REF] Murthy | Numerical simulation of welding and quenching processes using transient thermal and thermo-elasto-plastic formulations[END_REF]:

• Surface heat input condition through rate of heat input for unit surface area I q (applicable for welding to model the energy input from weld electrode).

• Surface heat loss due to convection and radiation:

( ) ( ) ( )( ) 1 2 1 2 c r L L L a a a q q q α θ θ α θ θ α α θ θ = + = - + - = + - (2)
Heat loss from the surfaces at temperatures near the melting point to ambient is predominantly due to radiation, and is given by [START_REF] Zhu | Effects of temperature-dependent material properties on welding simulation[END_REF][START_REF] Zhu | Numerical simulation of transient temperature and residual stresses in friction stir welding of 304L stainless steel[END_REF][START_REF] Murthy | Numerical simulation of welding and quenching processes using transient thermal and thermo-elasto-plastic formulations[END_REF]:

( ) ( )( ) ( ) 4 4 2 2 2 2 2 r L sb a sb a a a q A A εσ θ θ εσ θ θ θ θ α θ θ   = - = + - = -   (3)
Therefore, in order to consider heat convection and radiation on the plate surfaces, the heat flux loss is evaluated by the following expression [START_REF] Murthy | Numerical simulation of welding and quenching processes using transient thermal and thermo-elasto-plastic formulations[END_REF][START_REF] Siegel | Thermal Radiation Heat Transfer[END_REF]:

( ) ( ) ( )( ) 4 4 1 L a sb a a q A h α θ θ εσ θ θ θ θ θ = - + - = - (4)
where ( ) h θ is given by:

( ) ( )(
)

2 2 1 2 1 sb a a h A θ α α α εσ θ θ θ θ = + = + + + (5)
In this calculation, for the duplex stainless steel plates, ( ) h θ W/m 2 K is obtained through our experimental tests. In addition, a fusion latent heat of 500 J/kg K is taken in order to consider the phase transformation of this steel. This last value was obtained from the Outokumpu product specification data and other previous research works [START_REF] Pilipenko | Computer simulation of residual stress and distortion of thick plates in multi-electrode submerged arc welding[END_REF]. With respect to the austenitic stainless steel, a fusion latent heat of 330 J/kg K is used in order to consider its phase transformation [START_REF] Lindgren | Deformations and stresses in welding of shell structures[END_REF]. The variational formulation of this problem permits that the surface boundary conditions stated above may be treated implicitly together with the equilibrium equation. In this sense, the solution is obtained applying the principle of virtual temperatures [START_REF] Murthy | Numerical simulation of welding and quenching processes using transient thermal and thermo-elasto-plastic formulations[END_REF][START_REF] Bathe | Finite Element Procedures[END_REF][START_REF] Cook | Concepts and Applications of Finite Element Analysis[END_REF][START_REF] Chandrupatla | Introduction to Finite Element in Engineering[END_REF][START_REF] Mills | Heat Transfer[END_REF][START_REF] Del Coz Diaz | Nonlinear thermal optimization of external light concrete multi-holed brick walls by the finite element method[END_REF][START_REF] Siegel | Thermal Radiation Heat Transfer[END_REF].

Mechanical analysis 3.2.1. Plasticity

On the one hand, the plastic strain in both types of stainless steels fulfil the Von Mises yield criterion with isotropic strain hardening as well as its corresponding associated flow rule. Further, we have chosen the bilinear isotropic plasticity to model both types of steels. The rate relationship of different components between thermal stresses, σ , and strains, ε , is given by the following expression [9-10]:

( ) ≥ for plastic deformation ( 0 γ > ). On the other hand, in our structural thermo-mechanical analysis, a cut-off temperature around 2/3 of the melting temperature (see sub-section 4.3, entitled Mechanical analysis) is used in the numerical calculations to reduce unnecessary computational time [START_REF] Camilleri | Alternative numerical techniques for distortion of thin plate due to fillet-welded stiffeners[END_REF]. To save computational time without losing accuracy of results, different time steps are used for the two analyses. The initial time steps used in this study are 0.15s and 30s for the heat transfer analysis (heating and cooling phases respectively), and 0.2s and 15s for the stress-strain analysis (heating and cooling phases respectively). As such, the CPU time in a typical welding simulation is about 18min for the temperature computation and 256min for stress-strain computation in a welded seam between two plates whose size is 300 200 3 × × mm (length, width and thickness, respectively). The numerical models have been solved on a workstation computer with a CPU Intel Core2 Quad Q6600, 8 GB RAM memory and 4.5 TB hard disk. The total elapsed CPU average time for each model was 274min.

1 ij ij kk ij ij a s E E t ν ν α θ ε σ σ δ γ α θ θ θ + ∂ ∂   = - + + + -   ∂ ∂   & & & (6

Large displacements

In many problems it has been implicitly assumed that both displacements and strains developed in the structure are small. In practical terms this means that the geometry of the elements remains basically unchanged during the loading process and that firstorder, infinitesimal, linear strain approximations can be used. In practice such assumptions fail frequently even though actual strains may be small and elastic limits of ordinary structural materials are not exceed. If accurate determination of the displacements is needed, geometric non-linearity may have to be considered in some structures [START_REF] Zienkiewicz | The Finite Element Method[END_REF][START_REF] Del Coz Diaz | Non-linear analysis of the tubular 'heart' joint by FEM and experimental validation[END_REF]. Conversely, it may be found that a load is reached where deflections increase more rapidly that predicted by a linear solution and indeed a state may be attained where load-carrying capacity decreases with continuing deformation. This classic problem is that of structural stability and obviously has many practical implications. The applications of such an analysis are clearly of considerable importance in aerospace engineering, design of radio telescopes, cooling towers, box girders bridges and other relatively slender structures. In many cases very large displacements may occur without causing large strains. Typical in this context is the classical problem of the 'elastica' of which an example is a watch spring. Whether the displacements (or strains) are large or small, equilibrium conditions between internal and external 'forces' have to be satisfied. Thus, if the displacements are prescribed in the usual manner by a finite number of nodal parameters a , we can obtain the necessary equilibrium equations using the virtual work principle [START_REF] Bathe | Finite Element Procedures[END_REF][START_REF] Cook | Concepts and Applications of Finite Element Analysis[END_REF][START_REF] Chandrupatla | Introduction to Finite Element in Engineering[END_REF][START_REF] Zienkiewicz | The Finite Element Method[END_REF]:

( ) 0 T B d Ω Ψ = Ω -= ∫ a σ f (7)
in which B is defined from the strain definition ε as:

d B d = ε a (8)
The bar suffix has now been added for, if displacements are large, the strains depend non-linearly on displacement, and the matrix B is now dependent on a . We see that it can be conveniently written [START_REF] Zienkiewicz | The Finite Element Method[END_REF]:

( ) 0 L B B B = + a (9)
in which 0 B is the same matrix as in linear infinitesimal strain analysis and only L B depends on the displacement. In general, L B will be found to be a linear function of such displacements. Clearly, the solution of Eq. ( 7) will have to be approached iteratively. If, for instance, the Newton-Raphson process is to be adopted we have to find the relation between da and Ψ d . Thus taking appropriate variations of Eq. ( 7) with respect to da we have [START_REF] Zienkiewicz | The Finite Element Method[END_REF]:

T T T d dB d B d d K d Ω Ω Ψ = Ω + Ω = ∫ ∫ σ σ a (10)
and using equation d D d = σ ε and Eq. ( 8) it is obtained:

d D d D B d = = σ ε a (11)
and taking into account the Eq. ( 9), it is verified that

L dB B d =
. Therefore,

∫ + = + = Ψ V T L a d K a d K a d K dV dB d r r r r σ σ (12)
where

∫ + = = V L T K K dV B D B K 0 ( 13
)
in which 0 K and L K are given by [START_REF] Bathe | Finite Element Procedures[END_REF][START_REF] Zienkiewicz | The Finite Element Method[END_REF]:

∫ = V T dV B D B K 0 0 0 , ( )
0 0 T T T L L L L L V K B D B B D B B D B dV = + + ∫ (14)
To summarize, Eq. ( 12) can be expressed globally as [START_REF] Zienkiewicz | The Finite Element Method[END_REF]:

( ) 0 L T d K K K d K d σ Ψ = + + = a a (15)
Newton-type iteration can be once more applied precisely in order to solve the final non-linear problem.

Birth and death technique

If material is added to a system, certain elements in the FEM model may become existent. In such cases, it is possible to use element birth and death technique to deactivate or reactivate some finite elements. In this paper, to simulate the filling deposition of TIG welding procedure, the previous technique is chosen [8]. These elements are collected in groups according to the torch's movement, as it is shown in Fig. 3. To achieve the "element death" effect, the numerical procedure does not actually remove "killed" elements. Instead, it deactivates them by multiplying their stiffness by a reduction factor of 6 1.0 10 - × [START_REF] Madenci | The Finite Element Method and Applications in Engineering using ANSYS[END_REF]. In the first place, the elements describing the molten material are selected and deactivated. Next, they are reactivated in the thermal analysis according to the actual welding speed of the TIG manufacturing process. Fig. 3. Element birth and death sequence in the TIG welding process.

Elements loads associated with deactivated elements are zeroed out of the load vector.

In case of the TIG welding procedure, the deactivated elements mass and energy are not considered. Furthermore, the element's strain is also set to zero as soon as that element is killed. When the elements are activated, the heat input generation is applied as a load in the thermal FEM model, and the thermal strains are computed for newly-activated elements based on the current load step temperature and the reference temperature.

Finite element analysis

Welding processes are very difficult to model due to the complicated physics involved, which includes phenomena such as fluid mechanics with phase changes, shrinkage and porosity, macrosegregation among others. Each of the phenomena is a subject of intense research [START_REF] Camilleri | Alternative numerical techniques for distortion of thin plate due to fillet-welded stiffeners[END_REF]8,[START_REF] Sergey | Numerical simulation of filling and solidification of permanent mold castings[END_REF][START_REF] Pathak | Effects of mould filling on evolution of the solid-liquid interface during solidification[END_REF]. Taking into consideration even a few of these phenomena simultaneously can make welding analysis an extremely difficult problem. To make numerical simulation of welding, a number of simplification assumptions have to be made.

The present research paper deals with a finite element analysis that was carried out in two steps. A non-linear transient thermal analysis was conducted first to obtain the global temperature history generated during the welding and cooling processes. Then, a structural analysis was developed with the temperatures obtained from the thermal analysis used as loading to the stress model. The advantage of a decoupled approach is computational simplification and efficiency. Because of the thermo-mechanical nature of welding process theoretical considerations can be assessed either using a thermal or a mechanical analysis. Welding thermal stresses are calculated from the temperature distributions determined by the thermal model. The residual stresses from each temperature increment are added to the nodal point location to determine the updated behaviour of the model before the next temperature increment.

A simplified flowchart corresponding to the numerical analysis is shown in Fig. 4. Fig. 4. Flowchart corresponding to the numerical analysis. Welding in large unrestrained structures can cause significant distortions in modes that simplified two dimensional analyses of the weld region may not capture. Specifically, the parts may move relatively to each other, and the weld may be placed in a different location than that of the undeformed configuration, for these reasons, in this simulation, a three dimensional finite element analysis was used to simulate the welding process.

Finite element meshing

The mesh used in the stress analysis was identical to that in the thermal analysis. Fig. 5 shows the entire FEM model of the butt weld. To increment the computer efficiency, three zones with different mesh densities are used. The finest mesh corresponds with the 'fusion zone', whereas the coarsest corresponds with the 'base metal'. Between them, a transition mesh was imposed with a properly element size which matches the 'heat affected zone'. Fig. 5. Entire FEM model of the butt weld. Eight-node brick elements with linear shape functions are used in the model meshing. For thermal analysis, the element type is SOLID70 which has single degree of freedom, temperature, on its each node. For structural analysis the element type is SOLID185 with three translational degrees of freedom at each node [START_REF] Madenci | The Finite Element Method and Applications in Engineering using ANSYS[END_REF].

Thermal analysis

The heat source movement has been simulated with a thermal load applied to the elements that represent the fused material at a given time. The birth and death technique has been used. Figs. 3 and5 show the finite elements that simulate the filling material, these elements are removed at the beginning of the thermal analysis and then, for every load step, a cross sectional strip of elements are activated and loaded to simulate the passage of the heat source and the filling material deposition. The arc heat input was applied in the welding zone using a piecewise volumetric heat flux distribution function. Thermo-fluids phenomena are neglected in this direct approach, and conduction, convection and radiation are simulated and the thermo-physical characteristics are considered as temperature dependent. On the one hand, the effect of latent heat of fusion is taken into account by use of an enthalpy method. On the other hand, it is assumed that the thermal conductivity is constant after melting, whose value is three times that corresponding to this melting temperature (see Fig. 1 lower). In the thermal analysis, a total of 33 load steps with a maximum of 23 sub-steps and a minimum of 5 sub-steps were required to complete the heating cycle. Automatic time stepping and full Newton-Raphson method were used in each time step for the heat balance iteration.

Mechanical analysis

The material elastic-plastic behaviour has been implemented in the finite element code using a piece-wise linear fit to the true stress-plastic strain curves. The simplified model is presented in Figs. 2c and2d for duplex and austenitic stainless steels, respectively. Rate-independent plasticity is modelled using the incremental theory with a von Mises yield surface, associated flow rule and isotropic hardening. The convergence´s problems coming from the small tangent modulus at high temperature can be overcome using a meshing parameter of 2mm value and activating seven axial elements for a single drop in the welding zone. The material is assumed to present zero stiffness for temperatures exceeding a specified temperature, named "cut-off" temperature. Above this temperature any thermal strain is not considered [START_REF] Camilleri | Alternative numerical techniques for distortion of thin plate due to fillet-welded stiffeners[END_REF]. Moreover, the effect of specifying different cut-off temperatures has been investigated, as it is indicated in Table 2. Fig. 6 shows the effect of the cut-off temperatures variation in the maximum vertical displacement as well as the total CPU time for both stainless steels.

On the one hand, a numerical convergence problem can occur when the cut-off temperature increases, as it is possible to see in the total CPU time graph. The origin of this numerical problem can be explained by the appearance of excessive distortions due to a huge reduction of the material stiffness for such high temperatures. On the other hand, when the cut-off temperature is varied in the range considered in this study (see Table 2), a maximum vertical displacement is obtained. Consequently, the cut-off temperature values shown in Table 2 are used in the subsequent numerical calculations.

Fig. 6. Cut-off temperature effects for both stainless steels: total CPU time (left) and maximum vertical displacement (right). As it is shown in Fig. 7, the data from the temperature field is applied via a loop comparing the nodal temperature in all nodes for each time step. If the nodal temperature is greater than the cut-off temperature, then the temperature in this node is fitted to this value. Once all nodal temperatures are defined, the non-linear, stepwise static structural analysis is solved for current time step until all elements of the model have been activated. Fig. 7. Flowchart corresponding to the structural analysis

Results and discussion

In order to check the numerical results obtained from the finite element model, several experimental measurements have been made. Simple criteria have been followed to develop the testing methodologies. These methodologies are very useful to validate the simulation tool from a practical point of view. All of them and the main results obtained with the corresponding comparisons are explained as follows.

Experimental and numerical thermal results

To validate the numerical results of the thermal analysis, a set of thermographies have been taken at different times. In this way, a sequence of the different thermal states of the piece along the welding process can be observed. The thermographic principle is based on capturing infrared emission from the subject under investigation. For this purpose, a FLIR ThermaCAM™ E2 thermographic camera with a focal plane detector type (array of 160×120 pixels) has been employed. This device captures infrared radiation at wavelengths in the range of 7.5 to 13µm and temperatures ranging from 253 to 523 K. The accuracy of the measured temperature is ±2K and the value adopted for the material emissivity was 0.85 [START_REF] Wolfe | The Infrared Handbook[END_REF]. The thermographies were taken at a focus distance of 1m. It is possible to directly compare the pictures obtained from the thermographic camera with the thermal fields obtained from the numerical analysis at each step of the process. Fig. 8. shows the results at the following time steps: 274s, 780s, and 1,300s, respectively. Fig. 8. Comparison between experimental thermographies (left, a, c, e) and numerical results (right, b, d, f). Fig. 8 also shows the comparison between numerical and experimental results for different time steps. A good agreement is observed for the temperature fields in several locations. The small differences are due to material emissivity changes, caused by the variation of the temperature field and surface texture during the welding process. However, the maximum and minimum temperature values are a very close match.

Experimental measurements and structural analysis results

The results of the structural analysis are the practical and main purpose of this work. From this analysis it is possible to find out the stress, strains and displacements in both plates and weld pool when they are subject to a suddenly thermal process as welding. In order to validate the numerical FEM displacements, a DEA Global coordinate measurement machine (CMM) type 09. [START_REF] Mills | Heat Transfer[END_REF].08 has been used (see Fig. 9). The main technical characteristics of this device are the following:

• Measuring range: 900 × 1,500 × 800mm in X, Y and Z directions.

• Accuracy: from 1.4 to maximum length divided by 333µm (1µm in this case). Fig. 9. Experimental setup of the coordinate measuring machine (CMM): overall view (left) and a detail (right). Fig. 10 shows a set of pictures corresponding to the deformed shape of the plates after being welded compared with the numerical results explained here. Fig. 10. Comparison between the experimental deformed shape (left, a, c, e, g) and the numerical results (right, b, d, f, h). In order to validate the structural analysis, vertical displacements at nine different locations have been measured over a plate. These locations correspond to three transversal sections of the plate (two of them in the vertex and the others in the middle) and, over each of them, at three distances from the weld pool (see Fig. 11). Fig. 11. Geometrical description of the locations considered for results comparison. The numerical method developed in the current work helps to study the evolution of the welded plates' displacements. As in a common numerical simulation, it is possible to use two different ways of post-processing:

• Post-I: Results at specific points in the model over all time steps.

• Post-II: Overall model results at a specific time. Fig. 12. Vertical displacement evolution at nine different locations for both specimens: duplex stainless steel plates (left) and austenitic stainless steel (right).

With the basis of Post-I procedure, Fig. 12 shows the vertical displacement evolution at the previously described locations during the whole welding process for both materials considered: duplex and austenitic stainless steels. It is possible to observe that the absolute value for the maximum vertical displacement is similar for both steels but the sign is opposite. This means that the overall deformations of the welded plates are opposite for these two materials, as it was well observed in the experimental tests. Fig. 13. Numerical FEM displacement results for duplex (left) and austenitic (right) stainless steels in three different time steps. Through a combination of both types of post-processing methods it is possible to compare experimental and numerical results. On the one hand, with the Post-II procedure, vertical displacements at the selected locations at the final instant of the cooling stage are obtained. On the other hand, the displacement values at these locations were experimentally obtained by means of the CMM previously indicated. The numerical and experimental comparison is shown in Fig. 14. A very good agreement between them is obtained. Fig. 14. Comparison between vertical displacement values at selected locations for duplex (left) and austenitic (right) stainless steel. Numerical simulation vs experimental results.

Conclusions

Firstly, a computer code is developed for transient thermal and thermo-elasto-plastic stress coupled analyses considering various important nonlinearities due to temperaturedependent material properties. Secondly, the methodology was verified through comparison with experimental tests for different stainless steels. Thirdly, from this research work, it is possible to conclude that:

• The finite element method including the birth and death technique is an efficient procedure to analyze the vertical displacement evolution at different plate locations in the welding process of stainless steels, as well as their angular deformation and longitudinal bending. • The use of thermographic imaging in the context of TIG welding has been successfully demonstrated on a practical scale and helps to compare the thermal results in the subsequent numerical simulation at different locations and a specific time. The maximum and minimum temperature values in the numerical simulation and the experimental setup are a very close match. • The experimental setup of the coordinate measuring machine is appropriate to quantify vertical displacements in the welding plates at different locations. The numerical and experimental comparison shows a very good agreement between them. • The effect of specifying different cut-off temperatures has been investigated and the effect in the total CPU time and the maximum vertical displacement was calculated. As final conclusion, the optimum value in the calculations for this parameter was about 2/3 of the fusion temperature, specifically 1,173 K (for duplex steel) and 1,073 K (for austenitic steel). • This paper shows the important influence of the material properties in the final deformed shape of two different stainless steels, noting that opposite angular deformation and longitudinal bending occur when these steels are welded with a similar TIG procedure. • Finally, the methodology developed here can be applied to other types of welding processes and other different materials with success.

B 1 C

 1 A surface area of the body through which heat flows, m 2 a nodal displacements vector, m B global strain shape function, dependent on a, m -1 0 infinitesimal strain shape function, m -1 L B large strain shape function, dependent on a, m -specific heat, temperature-dependent, J/kg K D elastic constant matrix, temperature-dependent, Pa E Young's modulus temperature-dependent, Pa f external forces vector, N h equivalent film coefficient, W/m 2 K matrix, N m -1 K σ initial stress or geometric stiffness matrix, N m -1

, m 1 α 2 α

 12 convective heat transfer, W/m 2 K radiative heat transfer, W/m 2 K α temperature-dependent thermal expansion coefficient, K -1 ij δ Kronecker's delta, dimensionless ε emissivity of the body Ω (dimensionless) ij ε& strain rate tensor, s -1 γ plastic flow factor, Pa -1 s -1 λ material thermal conductivity, temperature-dependent, W/m K ν Poisson's ratio, dimensionless and temperature-dependent θ temperature of the studied body Ω , K

Fig. 1 .Fig. 2 .

 12 Fig. 1. Thermo-physical material properties for duplex and austenitic stainless steels: specific heat (upper) and thermal conductivity (lower).

Fig. 3 .

 3 Fig. 3. Element birth and death sequence in the TIG welding process.

Fig. 4 .

 4 Fig. 4. Flowchart corresponding to the numerical analysis.

Fig. 6 .

 6 Fig. 6. Cut-off temperature effects for both stainless steels: total CPU time (left) and maximum vertical displacement (right).

Fig. 7 .

 7 Fig. 7. Flowchart corresponding to the structural analysis.

Fig. 9 .

 9 Fig. 9. Experimental setup of the coordinate measuring machine (CMM): overall view (left) and a detail (right).

Fig. 10 .

 10 Comparison between the experimental deformed shape (left, a, c, e, g) and the numerical results (right, b, d, f, h).

Fig. 11 .

 11 Fig. 11. Geometrical description of the locations considered for results comparison.

Fig. 12 .

 12 Fig. 12. Vertical displacement evolution at nine different locations for both specimens: duplex stainless steel plates (left) and austenitic stainless steel (right).

Fig. 14 .

 14 Fig. 14. Comparison between vertical displacement values at selected locations for duplex (left) and austenitic (right) stainless steel. Numerical simulation vs experimental results.

Table 1

 1

	Welding parameters
	Item	Units	Duplex steel	Austenitic steel
			EN-1.4462	EN-1.4404
	Welding process efficiency	%	75		85
	Voltage	V	23		13
	Current intensity	A	180		115
	Welding speed	m/s	1.67 10 -×	3	3 4.8 10 -×

Table 2

 2 Cut-off temperatures determination Item

	f [K] θ	θ [K] cut	Ranging cut-off temperatures [K]


F Illán

M Alarcón

Numerical analysis of combustion and transient heat transfer processes in a two stroke SI engine

Keywords: Two zone combustion model, Spark ignition engine, Heat transfer, Network simulation model

A zero-dimensional model is presented to simulate the transient processes occurring within a two stroke SI engine. A two zone combustion model, with a spherically expanding flame front originating from the spark location, is applied. The model is numerically solved using the network simulation model which allows coupling the combustion model with a heat transfer model where both radiant and convective heat contributions have been taken into account for the in-cylinder gases. The boundary conditions for this model are the convective heat transferred to the cooling medium. A gas mixture model has been used to obtain the influence of working fluid properties on combustion development.

Nomenclature a

Wiebe's efficiency parameter (-) Combustion modelling in spark ignition (SI) engines is an established technique for improving engine performance. There are two approaches to combustion modelling:

A
thermodynamic and multidimensional models. The simplest type of thermodynamic model is termed zero-dimensional model and is ideal for general parametric studies. It can be used if the burn rate is known or by the use of semi-empirical relationships to express the burn rate; other thermodynamic models, termed quasi-dimensional models, attempt to simulate actual flame propagation and rely on the prediction of fundamental quantities such as eddy size, turbulence intensity and entrainment velocity. Multidimensional models attempt to compute fluid motion by numerically solving the partialdifferential equations which represent the basic principles of conservation of mass, momentum, energy and chemical species; these models, owing to their complexity, are impractical for basic engine design considerations [START_REF] Reid | A closed cycle simulation model with particular reference to two-stroke cycle engines[END_REF].

Zero-dimensional models have been used in a large number of studies in order to study different engine parameters: Descieux and Feidt [START_REF] Descieux | One zone thermodynamic model simulation of an ignition compression engine[END_REF] used a one-zone model to obtain the influence of several engine parameters on the power and efficiency of an air standard diesel engine; Guezennec and Hamama [START_REF] Guezennec | Two-zone heat release analysis of combustion data and calibration of heat transfer correlation in an I[END_REF] used a two-zone model to extract the heat release rate from pressure data for a SI 6-cylinder engine; Hajireza et al. [START_REF] Hajireza | A three-zone model for investigation of gas behavior in the combustion chamber of S.I. engines in relation to knock[END_REF] used a three-zone model in order to study the gas thermodynamic characteristics and their relation to knock in SI engines. A comprehensive assessment of various quasidimensional models was presented by Agarwal et al. [START_REF] Agarwal | Assessment of single-and two-zone turbulence formulations for quasi-dimensional modeling of spark-ignition engine combustion[END_REF] and Verhelst and Sheppard [START_REF] Verhelst | Multi-zone thermodynamic modelling of sparkignition engine combustion -An overview[END_REF].

The heat release process is, together with combustion, the process most studied in internal combustion (IC) engines. For this reason there exist a large number of studies which treat jointly both processes. Nevertheless, in most cases, to obtain the heat release, wall temperature or heat transfer rate were imposed [START_REF] Guezennec | Two-zone heat release analysis of combustion data and calibration of heat transfer correlation in an I[END_REF][START_REF] Hajireza | A three-zone model for investigation of gas behavior in the combustion chamber of S.I. engines in relation to knock[END_REF][START_REF] Jennings | A Computational Study of Wall Temperature Effects on Engine Heat Transfer[END_REF][START_REF] Lu | Studies on the cyclic variations of single cylinder twostroke engines -cycle analysis[END_REF], or wall temperature was obtained as a function of engine operation parameters [START_REF] Wu | Heat transfer model for small-scale aircooled spark-ignition four-stroke engines[END_REF][START_REF] Wu | Heat transfer model for smallscale spark-ignition engines[END_REF]. Other works [START_REF] Esfahanian | Thermal analysis of an SI engine piston using different combustion boundary condition treatments[END_REF][START_REF] Morel | Detailed analysis of heat flow pattern in a piston[END_REF] are mainly focussed on the piston and do not analyse other engine parts.

This work presents a numerical study of combustion and heat transfer processes in a geometrically simplified two-stroke crank case scavenged spark-ignition engine.

Assuming as known the burning law and the ambient temperature, the model allows obtaining the heat transfer characteristics of a 49 cm 3 air-cooled SI engine, including the temperature distribution in the whole engine, from inner cylinder wall to fins wall. The temperature maps obtained allow finding those parts of the engine subjected to high temperatures, which is useful to establish the number and distribution of fins in crankcase and block, in an air refrigerated engine, as well as to design surface coating of inner engine surfaces, including piston, or lubrication system.

Analysis has been made for the transient period comprised between the engine start and the moment when a steady temperature distribution is reached. As a result of this analysis, the engine p-θ diagram has been obtained, which contains relevant information about engine run, such as the evolution along the start of the engine (including transient period), which is very important on engine warm up and pollutant emissions, or maximum pressure in cylinder, which alerts in engine knock prevention. The numerical technique used to solve the problem is the Network Simulation Method (NSM) whose efficiency has been established for different non-lineal problems, both in the heat transfer domain [START_REF] Alarcón | Timedependent heat transfer in a fin-wall assembly. New performance coefficient: thermal reverse admittance[END_REF][START_REF] Zueco | Thermophoretic hydromagnetic dissipative heat and mass transfer with lateral mass flux, heat source, Ohmic heating and thermal conductivity effects: Network simulation numerical study[END_REF], and other fields [START_REF] Bég | Numerical study of magnetohydrodynamic viscous plasma flow in rotating porous media with Hall currents and inclined magnetic field influence[END_REF].

-Problem statement and governing equations.

The engine simulated in this work is a 49 cm 3 crankcase scavenged two-stroke SI engine. Simulation has been made at full load. Table 1 summarizes the engine's main features and operation conditions.

The problem in this work has been approached from a predictive point of view which requires the consideration of all relevant processes in the form of a combustion model [START_REF] Tinaut | Utilization of a quasi-dimensional model for predicting pollutant emissions in SI engines[END_REF]. Wiebe function has been used to obtain mass fraction burned (MFB):

(

where a (efficiency parameter) and n (form factor) are constants, θ is the crank angle during the combustion, θ 0 is the crank angle of ignition (25º bTDC) and ∆θ is the duration of the combustion (60º). According to Heywood [START_REF] Heywood | The Two-Stroke Cycle Engine. Its Development, Operation and Design[END_REF], the curve of MFB has been fitted with the values recommended by Sher for small two-stroke cycle SI engines (a = 3 and n = 5). Changing either a or n will alter the shape of the curve significantly.

In the case of comparison to the experimental results these parameters must be adjusted to the engine.

Only two zones, burned and unburned gases, are considered in this work, neglecting the influence of the thermal boundary layer and the flame front thickness. As the flame front is propagating from burned to unburned gas region, the volume of the burned gas increases and consequently the volume of the unburned gas decreases. According to Heywood [START_REF] Heywood | Internal Combustion Engine Fundamentals[END_REF], the density ratio is close to 4 for most SI engine operation conditions:

(

The gases in both zones are assumed to follow the ideal gas law [2-6, 8-10, 18-19].

Combustion is assumed complete (no unburned gases are present in burned zone) and there is neither exhaust gas recirculation (EGR) nor residual gases. The influence of residual gases has been taken into account by assuming a scavenging efficiency of 0.5, calculated applying the perfect mixing model, with a delivery ratio of 0.6 (charging efficiency equal to 0.45). Since the lubrication oil is mixed with the gasoline at a very low mixing ratio, its influence in the combustion process has been also neglected.

The unknown values for the model proposed are temperature, pressure, mass and volume for both the unburned and the burned zones and the mass flow from unburned to burned zone. Following the traditional two-zone engine modelling approach the analysis is based on applying the mass and energy conservation principles for the two zones, supplemented with the state equation in each zone, the mass fraction burned and the volume and area information (based on the engine geometry and kinematics). Since these equations are applicable at each crank angle, it is very easy to formulate them in a differential form as a function of crank angle:

(

where the instantaneous cylinder total volume depends on the piston position and is calculated by the following equation: [START_REF] Wu | Heat transfer model for smallscale spark-ignition engines[END_REF] The areas for heat transfer to the walls from the burned and unburned zones are obtained from equations ( 1), ( 2), ( 9) and [START_REF] Wu | Heat transfer model for smallscale spark-ignition engines[END_REF] assuming that the flame front is perfectly spherical and centred in the plug. The main engine parameters are summarized in Table 1.

Assuming that any thermodynamic property of the system (ϕ) can be expressed as a function of temperature, pressure, and relative air fuel ratio: [START_REF] Esfahanian | Thermal analysis of an SI engine piston using different combustion boundary condition treatments[END_REF] and operating with equations (1) to (9) it can be obtained:

In order to simplify equation [START_REF] Morel | Detailed analysis of heat flow pattern in a piston[END_REF] there are some assumptions that can be made for the system analysed:

• Relative air fuel ratio remains constant for the fresh gases during the entire cycle.

• Molecular mass is different for burned and unburned gases, but remains constant in each zone during the entire cycle.

• Since the gases are assumed to follow the ideal gas law, their enthalpy can be obtained as a function of temperature only.

• For both the unburned and burned zones it can be established:

(13) (14)
Additionally, the heat addition due to combustion is considered in the burned zone. Its value can be obtained as: [START_REF] Bég | Numerical study of magnetohydrodynamic viscous plasma flow in rotating porous media with Hall currents and inclined magnetic field influence[END_REF] where Q LHV is the lower heating value of the fuel (gasoline, assumed to be C 8.26 H 15.5 ).

Taking into account all these assumptions, equation ( 12) can be considerably simplified, although it must be expressed in different terms for unburned and burned zones. For the unburned gases zone, it can be established: [START_REF] Tinaut | Utilization of a quasi-dimensional model for predicting pollutant emissions in SI engines[END_REF] whereas for the burned gases zone it can be expressed as:

The enthalpy variation was obtained assuming that the fresh gases are a homogeneous mixture of air (0.79 N 2 + 0.21 O 2 ) and fuel (C 8.26 H 15.5 ) in stoichiometric ratio and the burned gases are a mixture of N 2 , CO 2 and H 2 O (complete combustion).

The standard state enthalpy for species j is assumed to be temperature dependent and it takes the following form [START_REF] Heywood | Internal Combustion Engine Fundamentals[END_REF]: [START_REF] Heywood | Internal Combustion Engine Fundamentals[END_REF] where the constants a jk for N 2 , O 2 , CO 2 and H 2 O are given in [START_REF] Heywood | Internal Combustion Engine Fundamentals[END_REF] in two different temperature ranges, appropriate for unburned and burned mixture property calculations.

Enthalpy for the fuel is also fitted to a polynomial form: [START_REF] Tillock | Measurement and Modeling of Thermal Flows in an Air-Cooled Engine[END_REF] where θ = T(K)/1000, A f6 is the constant for the datum of zero enthalpy for C 2 , H 2 , O 2 and N 2 , at 298.15 K (for a 0 K datum, A f8 is added to A f6 ) and the values for coefficients A fj are given in [START_REF] Heywood | Internal Combustion Engine Fundamentals[END_REF].

The convective heat transfer coefficient in both zones is obtained by applying the model proposed by Rao and Bardon [START_REF] Rao | Convective heat transfer in reciprocating engines[END_REF]: [START_REF] Rao | Convective heat transfer in reciprocating engines[END_REF] where the characteristic velocity adopted to obtain Reynolds number is the turbulence intensity, u'.

To obtain the turbulence intensity, the following model [START_REF] Heywood | Internal Combustion Engine Fundamentals[END_REF] is adopted:

(21)
where the subscript 0 refers to the start of the combustion process. The turbulence intensity at this moment has been assumed equal to about half the mean piston speed.

In the burned gases zone, a radiation term proposed by Annand [START_REF] Heywood | Internal Combustion Engine Fundamentals[END_REF] is added to the convective heat transfer: [START_REF]Pspice 6.0[END_REF] where σ is the Stefan-Boltzman constant, is the mean gas temperature, T wi is the inner wall temperature and β is an adjust parameter for which a value of 0.6 has been adopted.

Therefore, the heat transfer rate in both unburned and burned zone is obtained as:

(23) (24)
To obtain the inner wall temperature, heat conduction through the cylinder walls is considered:

(25)

Finally, to obtain the outer wall temperature, convection between cylinder outer wall and cooling medium is considered:

(26)
where the subscript j refers to the cooling medium, T j is the temperature for cooling medium j, A wo,j and T wo,j are the area and temperature of the engine outer wall in contact with cooling medium j and h j is the convective heat transfer coefficient for cooling medium j. Three different fluids have been considered:

• Air at a temperature T a = 298.16 K; α a value depends on air velocity and external surface geometry. A constant value of 66 W•m -2 •K -1 obtained applying the Colburn correlation [START_REF] Chapman | Heat transfer[END_REF] for forced convection with an air velocity of 12 m/s was assumed in the entire external engine finned surface. The engine axis has been considered in vertical position (Fig. 2a) and the air flows perpendicular to this vertical axis.

• Crankcase gases at a temperature T cg = 333.16 K; α cg is equal to 415 W•m -2 •K -1 and was calculated applying the correlation proposed by Rao and Bardon [START_REF] Rao | Convective heat transfer in reciprocating engines[END_REF].

• Lubricating oil, which forms a separating film between piston and cylinder, has been assumed to be at a constant temperature T lub = 393.16 K; α lub is equal to 1668 W•m - 2 •K -1 and was calculated applying the correlation proposed by Sieder and Tate [START_REF] Chapman | Heat transfer[END_REF].

-Numerical approach.

The differential equations ( 16) and ( 17) under boundary conditions (26) have been solved using the network simulation method (NSM), a powerful numerical method extensively used in a diverse variety of thermal and fluid problems [START_REF] Alarcón | Timedependent heat transfer in a fin-wall assembly. New performance coefficient: thermal reverse admittance[END_REF][START_REF] Zueco | Thermophoretic hydromagnetic dissipative heat and mass transfer with lateral mass flux, heat source, Ohmic heating and thermal conductivity effects: Network simulation numerical study[END_REF][START_REF] Bég | Numerical study of magnetohydrodynamic viscous plasma flow in rotating porous media with Hall currents and inclined magnetic field influence[END_REF].

NSM is founded on the classical thermoelectric analogy between thermal and electrical variables. The starting point is a network circuit design, whose equations are formally equivalent to equations ( 16) and [START_REF] Heywood | The Two-Stroke Cycle Engine. Its Development, Operation and Design[END_REF]. The engine walls have been discretized and a sufficient number of networks are parallel and serially connected to form the engine walls. Heat transfer conditions ( 23) to (25) and boundary conditions (26) are added by means of special electrical devices. Once the network model is designed, its simulation is carried out through adequate circuit resolution software; PSpice® [START_REF]Pspice 6.0[END_REF] is the software selected for this study.

The electrical analogy relates the electrical potential to the temperature, while the electrical current is equivalent to the energy flow. According to this analogy, the term for unburned and burned gases in equations ( 16) and ( 17) must be implemented by means of a capacitor. Since depends on unburned (or burned) gases temperature and PSpice does not allow the defining of variable capacitors, equations ( 16) and ( 17) must be rewritten in the following form:

(27) (28)
In order to avoid convergence problems due to discontinuities during burning zone creation and unburned zone consumption, equations ( 27) and ( 28) are solved in an integrated way. Figure 1a schematically shows the network model corresponding to the integrated model for equations ( 27) and (28). A capacitor C of value 1 represents the temporary variation in gases temperature, , where is an averaged temperature which depends on burned and unburned gases temperature and mass proportion. All the terms on the right of equations ( 27) and ( 28 27) and ( 28) is implemented by means of a variable voltage source whose output voltage varies according to a polynomial function of in agreement with equations ( 18) and ( 19). This source is included in the definition of all the previously described current generators as an additional controlling voltage source.

Finally, to obtain , the network models defined for the combustion chamber must be connected to the engine walls. Figure 1b schematically shows how this connection is implemented. As is schematically represented in Figure 2b, the engine inner walls are considered to be divided into eleven different zones, namely Cells 1 to Cells 11. Each "j" wall zone is connected to two voltage controlled current generators and that represent the heat transferred from gases to "j" zone. In order to assess the piston motion, in Cells 6 to Cells 8 zones, the current generators vary not only with gases temperature but with piston position.

In order to simplify the definition of the model's geometry and reduce the computational cost, a 2-D discretization of the solid body has been built: cylindrical symmetry is assumed for the piston and cylinder geometry, whereas spherical symmetry is assumed for the cylinder head. Figure 2a shows a view of the simplified engine geometry. Annular geometry has been chosen for cylinder fins, and an innovative spherical fin has been implemented for the cylinder head. Extended surfaces of the actual engine do not have such geometry, but they have been adjusted to obtain the same exchange surface area as in the real engine, which assures in practical terms the validity of the heat transfer model.

As can be seen in Figure 2b, to discretize the entire engine, a total of 234 unitary networks are serially connected. Figure 3 shows the four types of unitary network used: a) middle symmetrical cell, b) symmetrical cell with convection, c) middle fin cell, and d) end fin cell (with convection on the end surface).

-Results and discussion.

In order to analyze the transient period, the simulation runs for a period of 300 seconds which, as can be seen in figures 4 to 6, is enough to assure that a stable temperature distribution is reached for the entire engine. Figures 4 to 6 show the temperature evolution at twelve different points as well as the thermal map obtained for the cylinder body, the cylinder head and the piston six minutes after engine start. As can be seen, for the piston the stable situation is reached in less than one minute, whereas in the case of the cylinder body and the cylinder head this time increases up to nearly six minutes due to the presence of fins.

The temperature distribution shown in Figure 4 for the cylinder body is in good agreement with that obtained by Šelih et al. [START_REF] Šelih | A novel cylinder cooling system of air-cooled engines[END_REF] or Wu and Chiu [START_REF] Wu | A study on the integration model of heat transfer in an engine fin system[END_REF]. As can be seen in Figure 5, due to its lower thermal conductivity, the temperature in the spark plug is higher than in the rest of the cylinder head. As expected, temperature decreases from the inner wall, in contact with combustion gases, to the fins walls in contact with ambient air. Finally, the temperature distribution shown in Figure 6 for the piston is in good agreement with that presented by Heywood for SI engines [START_REF] Heywood | Internal Combustion Engine Fundamentals[END_REF] and that obtained by previous researchers such as Esfahanian et al. [11] or Morel and Keribar [START_REF] Morel | Detailed analysis of heat flow pattern in a piston[END_REF].

The in-cylinder pressure obtained is shown in Figure 7 as a function of the crank angle.

As can be seen, only slight differences between cold engine and hot engine conditions are found. Once the temperature distribution in the entire engine reaches a stable situation, the peak in-cylinder pressure rises to near 26 bars, which represents only 2 % of increase with respect to the value obtained for cold engine conditions. Similarly, the maximum gases temperature obtained for hot engine conditions is 2101 K whereas for cold engine conditions the maximum temperature obtained is 2061 K.

Finally, Figure 8 shows the overall energy balance in the cylinder as a function of crank angle. According to this figure, around 30.9 % of total energy released during fuel combustion is transformed into useful work (the friction effect has not been taken into account) whereas around 38.6 % of energy added is lost in the exhaust gases enthalpy and the other 30.5 % represents the heat transfer rate to the cooling medium. The values obtained are in good agreement with typical values given in the specialized bibliography [START_REF] Heywood | Internal Combustion Engine Fundamentals[END_REF] and the energy balance (between energy input and energy output) is in near equilibrium, with only a slight deviation, lower than 1 %, attributable to the use of a relatively high data sampling rate. A closer energy balance could be reached by reducing the data sampling rate but, on the other hand, the computational cost would increase. In order to achieve a balance between accuracy and computational cost, the data sampling rate has been fixed at 1•10 -5 seconds.

-Conclusion.

A zero-dimensional combustion model consisting of two zones, coupled with a heat release model has been numerically solved in this paper. A gas mixture model has been used to obtain the influence of working fluid properties on engine performance. The model allows obtaining the in-cylinder temperature and pressure profile, the engine energy balance or the engine temperature map as well as other derived results. The results obtained are in good agreement with other results obtained by previous researchers, proving the ability of the model proposed in this work.

Approximately six minutes have been found as the duration of the transient period from start up in the whole engine, although the piston reaches the regime conditions in less than one minute. Referring to temperature, the most critical point of the engine has been found in the spark plug and its vicinity, more even than in piston head.

Using the model proposed in this work it is relatively simple to analyze, in further work, the influence of different operation parameters like cooling medium conditions, engine speed, combustion characteristics (start of combustion, total combustion duration, or Wiebe function adjust parameters). Finally, using a mesh model more adjusted to the engine geometry, the computational cost will increase, but it will be possible to analyze the influence of fin geometry on the engine temperature distribution.

  rod length to crank radius (r=l/a=2•l/L) R specific gas constant (J kg -1 K -1 ) wall in contact with burned gases w,u wall in contact with unburned gases wi inner engine wall wo outer engine wall 1. -Introduction.

  ), except the heat addition due to combustion, are implemented by means of voltage controlled current generators. In burned gases zone, the heat is transferred through the cylinder walls by means of two different mechanisms and therefore two different voltage controlled current generators are necessary to implement the first term on the right of equations (27) and (28): current generator simulates the convection heat transfer and therefore must be the controlling voltage for this source; current generator simulates the radiation heat transfer and therefore must be the controlling voltage for this other source. Additionally, the definition of these current generators takes into account a flame propagation model and the variations in the heat exchange coefficients. is the controlling voltage for the second term on the right of equations (27) and (28), whereas the variations in the value of are included in the source definition. The heat addition due to combustion only depends on the MFB and the lower heating value of the fuel; therefore the last term on the right of equation (28) is implemented by means of a source which varies according to the MFB. The denominator in all the terms on the right of equations (

Figure 1 .

 1 Figure 1. Schematic representation of the engine network model: 1.a. integrated combustion model corresponding to equations (27) and (28); 1.b. connection between the combustion model and the engine walls.

Figure 2 .

 2 Figure 2. Engine solid model (a) and mesh (b). Eleven different wall zones have been considered in order to assess the burned gases zone motion. Three different convective boundary conditions exist: ambient air at 298.16 K, crankcase gases at 333.16 K, and oil film at 398.16 K.

Figure 3 .

 3 Figure 3. Unitary networks used in the engine network model; a) middle symmetrical cell, b) symmetrical cell with convection, c) middle fin cell, and d) end fin cell (with convection on the end surface).

Figure 4 .Figure 5 .

 45 Figure 4. Cylinder body temperature map at t = 300 sec and temperature evolution at fin end (6), fin base[START_REF] Jennings | A Computational Study of Wall Temperature Effects on Engine Heat Transfer[END_REF] and cylinder body inner wall[START_REF] Lu | Studies on the cyclic variations of single cylinder twostroke engines -cycle analysis[END_REF]

Figure 6 .Figure 7 .

 67 Figure 6. Piston temperature map at t = 300 sec and temperature evolution at piston base[START_REF] Wu | Heat transfer model for small-scale aircooled spark-ignition four-stroke engines[END_REF], piston middle point[START_REF] Wu | Heat transfer model for smallscale spark-ignition engines[END_REF], piston crown outer edge[START_REF] Esfahanian | Thermal analysis of an SI engine piston using different combustion boundary condition treatments[END_REF] and piston crown center[START_REF] Morel | Detailed analysis of heat flow pattern in a piston[END_REF]

Figure 8 .

 8 Figure 8. Overall energy balance in the cylinder as a function of crank angle.

Table 1 .

 1 Engine parameters.

	Compression ratio (ε)		9.167
	Connecting rod length (l)		6.248 cm
	Ratio of connecting rod length
	to crank radius (r=l/a=2•l/L)	3.704
	Delivery ratio ( )		0.6
	Scavenging efficiency (	) 0.5
	Charging efficiency (	)	0.45


Yupeng Wu
email: yupeng.wu@nottingham.ac.uk

Guohui Gan

Anne Verhoef

Pier Luigi Vidale

Raquel Garcia Gonzalez

Experimental measurement and numerical simulation of horizontal-coupled Slinky Ground Source Heat Exchangers

Keywords: ATE Ground source heat pump, horizontal-coupled, slinky, specific heat extraction

Introduction

To reduce dependence on fossil fuel energy resources and environmental degradation resulting from their combustion, renewable energy sources must play an increasingly significant role. Geothermal energy, one of the renewable energy resources, can be used to provide electricity, heating and cooling for buildings. Shallow Ground Source Heat Pump systems (GSHP) which make use of the relatively constant temperature of the earth are mainly used for space conditioning. However, the installation of the ground Heat Exchanger (HE) is costly, which is one of the main causes preventing the wide adoption of the GSHPs in the UK. Therefore, the horizontal-coupled ground source heat pump has become increasingly important. It can reduce installation cost compared to that of other GSHP systems as no drilling is necessary and only a trench 1 to 2 m in depth is required. Nevertheless, A single tier ground loop requires a large area of ground to lay the pipe network. This problem can be alleviated to some extent by employing double tier pipe arrangement or slinky ground loop.

A slinky ground loop requires horizontal trench lengths of 20-30% of those for a horizontal coupled single pipe configuration, however the thermal performance of the slinky ground loop may be doubled. A first recorded slinky heat exchanger coupled GSHP system was developed by Bose and Smith [START_REF] Bose | Performance of new ground heat exchanger configurations for heat pump[END_REF] at Oklahoma State University. They indicated that the slinky heat exchanger enhanced the heat transfer area within a limited space when compared to a horizontal buried heat exchanger.

The performance of a heat pump is usually expressed in terms of its Coefficient Of Performance (COP) which is the ratio of energy output to supplied energy (electricity for the compressor, pump etc) of GSHP. COPs in the heating and cooling mode are indentified by COP h and COP c , respectively. A typical heat pump has a COP of around 4 which indicates that the heat pump produces four units of heating energy for every unit of electrical energy input. The COP of the GSHP depends on the soil type at the installation, etc. Research has been carried out to determine the performance of horizontal-coupled GSHP. Metz [START_REF] Metz | Ground coupled heat pump system experimental results[END_REF] used a horizontal-coupled GSHP to provide heating and cooling for a 104m 2 house in Long Island, New York. During a period of one year monitoring, the obtained average COP h and COP c of the GSHP system were 2.46 and 1.91, respectively. İnallı and Esen [START_REF] İnallı | Experimental thermal performance evaluation of a horizontal ground source heat pump system[END_REF] carried out experimental measurements of a GSHP system with horizontal-coupled heat exchanger in selected depths of 1 and 2m in Elaziğ, Turkey. After running the GSHP for heating from November 2002 to April 2003 for a room with 16.24m 2 floor area, the obtained average COP h of the GSHP system was 2.66 and 2.81 for the heat exchanger buried at depths of 1m and 2m, respectively. Doherty et al [START_REF] Doherty | Ground source heat pumpdescription and preliminary results of the Eco House system[END_REF] installed and undertook experimental measurements of a horizontal slinky GSHP system in the University of Nottingham, UK. An average COP h of approximately 2.7 was reported. Coşkun et al [START_REF] Coşkun | Experimental performance investigation of a horizontal ground source compression refrigeration machine[END_REF] studied the performance of a horizontal ground source compression refrigeration system in Bursa city, Turkey. The horizontal coupled heat exchanger was buried at a depth of 2.0m. The COP c of the overall system varied between 2 and 2.5.

Research has also been undertaken for the horizontal-coupled GSHP systems by other investigators including the experimental characterisations of the thermal performance of the horizontal coupled heat exchanger and different approaches to enhance the heat transfer between the heat exchanger and the soil [START_REF] Jones | Thermal performance of horizontal closed-loop ground coupled heat pump systems using flowable-fill[END_REF][START_REF] Piechowski | Heat and mass transfer model of a ground heat exchanger: validation and sensitivity analysis[END_REF][START_REF] Piechowski | Heat and mass transfer model of a ground heat exchanger: theoretical development[END_REF][START_REF] Esen | Numerical and experimental analysis of a horizontal ground-coupled heat pump system[END_REF][START_REF] Demir | Heat transfer of horizontal parallel pipe ground heat exchanger and experimental verification[END_REF][START_REF] Koyun | Experimental study of heat transfer of buried finned pip for ground source heat pump applications[END_REF]. Florides and Kalogirou [START_REF] Florides | Ground heat exchangers-A review of systems, models and applications[END_REF] reviewed the performance and numerical model of ground-coupled heat exchangers. They indicated that the ambient climatic conditions would affect the temperature profile below the ground surface and this should be considered when designing a horizontal coupled heat exchanger. In the UK, most of the recent studies of the horizontal coupled GSHP were focused on the heat pump output to maintain daily indoor air temperature [START_REF] Rawlings | Ground Source Heat Pump: A technology review[END_REF][START_REF] Singh | Factors influencing the uptake of heat pump technology by the UK domestic sector[END_REF][START_REF]Heat pumps in the UK -a monitoring report[END_REF][START_REF]Domestic ground source heat pumps: design and installation of closed-loop systems[END_REF]. Few studied the long term effect of horizontal-coupled heat exchanger on the ground temperature, which is vital for correct sizing of heat pumps to achieve an accurate prediction of their long term performance over their useful working life. In this work, the performance of a horizontal coupled slinky GSHP installed in the UK was monitored and measured results were used to validate a CFD model. The validated model was then used to predict the thermal performance of ground coupled slinky heat exchangers for different slinky diameters and slinky interval distances.

Experimental measurements

Experimental apparatus installation

Experimental measurements were carried out at Talbot Cottage, Drayton St Leonard site, Oxfordshire, UK where 4 parallel horizontal-coupled slinky Heat Exchanger loops were installed in a 80m long by 20m wide paddock area at a depth of around 1.2m below the ground surface. Figure 1 shows an image for one of the loops. The fluid in the heat exchanger was water-ethylene glycol (30% by weight) mixture. The ground surface was overgrown with wild flowers. Soft sand was used below and on top of the slinky heat exchanger.

To monitor the soil temperature distribution and also the performance of the horizontalcoupled heat exchanger, two holes were dug in the ground to install thermistors at various depths. One hole was dug above a portion of the ground coupled heat exchanger, the second hole served as a reference hole dug at a distance of around 2m from the other hole and with no heat exchanger installed around it. After installing the sensors, all the holes were refilled with the original soil. Monitoring started 1 month after switching on the GSHP system. A sketch of the positions of the sensors used to monitor the soil and the heat exchanger thermal performance with detailed sensor positions is shown in Figure 2. T-type thermocouples were installed at the inlet and outlet of the ground coupled heat exchangers to monitor the thermal energy extracted by the ground coupled slinky heat exchanger. Simultaneously, a mini weather station was installed to monitor the ambient air temperature, wind speed, global solar radiation, relative humidity and rain fall. All the sensors were calibrated before installation.

Analysis of experimental results

The sensor readings were taken every 30mins. Ambient air temperature and soil temperatures above a portion of heat exchanger and in the reference hole for a period of around two soil temperature at a depth of 0.25m had a similar characteristic to ambient air temperature. There were short term strong and irregular fluctuations of temperature, caused by the daily change of weather conditions. However, the temperature fluctuation ceased at a larger depth, below the depth of 0.5m, the temperature fluctuations were hardly discernible. It can be seen that soil with and without heat exchanger below had a similar temperature variation pattern at the same depth. However, the soil temperatures with a heat exchanger installed below were lower than those of the soil in the reference hole without heat exchanger under it. With the soil depth increase, the soil temperature increased without heat exchanger below it, whereas the soil temperature decreased where the heat exchanger was installed below it. This is due to the thermal energy being extracted by the heat exchanger. At a depth of 0.02m below the ground surface, there were no significant temperature differences for the soil with and without heat exchanger below them. Around 20th Dec 2009, the soil without heat exchanger installed below had a higher temperature of around 2°C than that with heat exchanger installed below. This deviation may be caused by the different soil physical properties or the heat exchanger.

At a larger depth, the temperature difference between the soil with and without heat exchanger installed was larger. The temperature difference at a depth of 0.5m was around 1.5°C, and 3°C at a depth of 1m. it can be seen that the soil with and without heat exchanger installed, their surface temperature had a similar characteristic to the ambient air temperature on that day, however, the ground surface temperature was slightly higher with no heat exchanger installed below.

This might be caused by different soil physical properties. They all decreased in the night from 00:00 to 8:00 due to convective and long wave radiative heat losses, and reached the lowest value of around 3°C of the day at 8:00 in the morning, and then increased to around 13.5°C (no HE) and 11.5°C(HE installed), respectively at noon. After 13:30 they all began to decrease due to the reduced solar radiation intensity. Because the thermal energy in the ground was extracted by the heat exchanger, the soil temperature varied between 7 and 8°C at a depth of around 1.14m below the ground. The environmental conditions influenced the soil temperature up to a depth of around 0.5m. It can also be seen that at the same depth of around 1m from the ground surface, the soil with no heat exchanger installed below had a constant temperature around 11°C for the whole day, which was over 3°C higher than the soil temperature with heat exchanger installed below. At a depth of 0.25m from the ground surface, the soil without heat exchanger installed below had a temperature around 1°C higher than that with heat exchanger installed below. It can be concluded that the extraction of thermal energy by the heat exchanger affected the soil temperature at a distance of around 0.9m from the heat exchanger during the operation period. and 6, it can be seen that the soil at a depth of 0.02m had a temperature around 0°C for the whole day, because the ground surface was covered with approximately 0.05m thick snow, which prevented the radiative and convective heat loss, and also the chance to absorb the diffused solar radiation at the ground surface. Its temperature was higher than the ambient air temperature at the same time. At a depth of around 1.0 m from the ground surface, the soil with no heat exchanger below had a constant temperature around 5°C for the whole day, it was approximately 2.5°C higher than that of the soil with heat exchanger installed below. At a depth of 0.1m from the ground surface, the soil with no heat exchanger below had a temperature around 0.5°C higher than that with heat exchanger installed below, similar to characteristics of the soil temperature distribution found on 7th Nov 2009.

Figure 7 illustrates the measured inlet and outlet fluid temperatures of the ground coupled heat exchanger. The temperatures of the inlet and outlet of the heat exchanger significantly decreased with time. However, there were no significant change in temperature difference between the inlet and outlet. The GSHP system had been continuously running for nearly two month by the time (nearly 24 hours a day). Only for a few hours (less than 1 hour a week), the temperature difference dropped to approximately 1°C from 2°C during a two month operation, because the GSHP system stopped running, and the inlet temperature for the heat exchanger began to increase for a short period. The calculated COP of the ground source heat pump system is shown in Figure 8. It can be seen that the system COP decreased slightly with the time from 6th Nov 2009 to 31st Dec 2009. It was around 2.7 on 6th Nov 2009, and decreased to approximately 2.5 on 31st Dec 2009.

Measurement of site soil properties

Three soil samples at depths of around 0.2, 0.8 and 1m below the ground surface at the test site were taken to the lab at the Department of Built Environment, University of Nottingham using 50mm deep by 50mm inner diameter aluminium cylinder tins. The thermal conductivity, thermal capacity, density and volumetric moisture content of the in situ undisturbed soil were measured using the following facilities KD2 pro used to measure the thermal conductivity and capacity, A thermostatically controlled oven, capable of maintaining a temperature between 105 and 110°C, A balance readable and accurate of 0.01g

The measured thermophysical properties of the soil samples are shown in Table 1. The measured thermophysical properties varied at different depths. The thermal conductivity and density were higher for the soil near the ground surface and heat exchanger than in the middle but the soil thermal capacity and the volumetric moisture content decreased with the depth.

Average thermophysical properties were used in the numerical prediction.

Numerical prediction

A commercial CFD software package FLUENT was used to predict the thermal performance of a portion of the horizontal coupled slinky and straight heat exchangers. The transient 3 dimensional sensible heat transfer model for the ground coupled heat exchanger is as follows:

c q T c k T V t T ) ( ) (  (1)
where is the density (kg/m 3 ), t is the time (s), T is the temperature (K), V is the fluid flow velocity (m/s), k is the thermal conductivity (W/mK), c is the specific heat (J/kgK) and q is the heat source (W/m 3 ).

The model was validated using the experimentally measured results using the method described below. The measured environmental conditions at 00:00 7th Nov 2009 were used as the initial conditions for the model. The predicted and the measured soil temperature distributions at selected times are shown in Figure 9. It can be seen that the predictions generally agreed with the experimental measurements. The maximum difference between the predicted and measured soil temperatures was less than 1°C, which was found near the ground surface. This might be the result of using constant and uniform soil thermophysical properties in the model, whereas the soil thermophysical properties might not be uniform particularly near the surface.

The thermal performance of a portion of a horizontal coupled slinky and straight heat exchangers was predicted using the validated CFD model. Comparisons of thermal performances of the slinky heat exchanger at different slinky (coil) diameters and different slinky (coil) interval distance were also conducted. The measured average thermophysical properties of the in situ soil detailed in section 2.3 were used in the numerical model. In the predictions, it was assumed that the horizontal-coupled straight and slinky heat exchangers were buried at a depth of 1.2m in the ground as illustrated in Further assumptions used to simplify the model were as follows,

For a short length of the heat exchanger, the pipe outer surface temperature was assumed constant at 1°C for heating mode. Prediction of the thermal performance of the heat exchanger with similar dimension and environmental conditions with 1°C

refrigerant fluid flow at the inlet of the heat exchanger was also undertaken, and the results showed that after running the system for more than 140 hours, there was no significant difference in the temperature distribution in the soil or the specific heat extraction between the conditions with fixed pipe temperature and fluid flow inside the heat exchanger, A constant temperature of 10°C was utilised at a depth of 4m from the ground surface for a short period of operating the GSHP,

The soil thermal properties (soil thermal conductivity and diffusivity) were constant and uniform,

The wind speed was 3m/s on the ground surface with an ambient air temperature of 5°C.

Predicted thermal performance of the horizontal-coupled straight and slinky heat exchanger

Isotherms generated from prediction for the horizontal-coupled straight and slinky heat exchanger at the elapsed time of 1 and 50hours are shown in Figures 11 and12, respectively.

For the horizontal-coupled straight heat exchanger, after running the system for 1 hour, the temperature of the soil around the heat exchanger decreased to 1°C. The temperature of the soil decreased at a distance within 0.1m from the central long axis of the heat exchanger, due to the thermal energy being extracted by the heat exchanger. At 50 hours, the heat exchanger affected the temperature of the soil within a distance of around 0.6m from the central long axis of the heat exchanger. With further running of the system, the soil temperature continuously decreased, due to the convective heat loss from the soil to the ambient air and heat extraction by the heat exchanger. The slinky heat exchanger and its surrounding soil had similar thermal characteristics as that for the straight heat exchanger under similar conditions.

However, the temperature of the soil within a distance of around 0.8m from the central long axis of the heat exchanger was influenced by the heat exchanger, it was larger than that of the single straight heat exchanger system. This is because the slinky system comprised of a straight and coil heat exchangers with a small portion of the heat exchanger being in parallel located.

A parameter used to calculate the required length of borehole heat exchangers is the specific heat extraction, expressed in Watt per meter length. Typical values are in the range from 40 to 70W/m for borehole heat exchangers, dependent on soil thermal conductivity, heat pump annual operation hours, number of neighbouring boreholes, etc [START_REF]European Geothermal Energy Council, Ground source heat pump[END_REF]. The variation of the predicted specific heat extraction with time for a portion of the horizontal-coupled straight and slinky heat exchangers, and the difference of the specific heat extraction between the two systems are illustrated in Figure 13. It can be seen that the specific heat extraction for both the horizontal-coupled straight and slinky heat exchangers was approximately at 46W/m initially, decreased to approximately 33 and 30W/m respectively after continuously running the GSHP system for 10hours, and further decreased to around 18 and 15W/m after running the GSHP system for 140hours. Thus, it can be concluded that the difference of the specific heat extraction between the straight and the slinky heat exchangers increased with running time. It was around 1.5W/m for the first hour, and increased to around 3W/m at 140hour. The predicted specific heat extraction was lower than that commonly used for designing GSHPs.

If the total length of the heat exchanger was estimated from the over optimistic value of the specific heat extraction, the output and COP of the heat pump would be lower than expected.

The variation of the heat extraction per meter length of the soil occupied by heat exchangers is illustrated in Figure 14. Although the specific heat extraction of the slinky heat exchanger was lower than that of the straight one, the heat extraction per meter length of the soil for the slinky heat exchanger was significantly higher than that of the straight heat exchanger. At the beginning, the difference between the slinky and straight heat exchanger was about 45W per meter length of soil, and decreased to approximately 12W per meter length of soil after running the simulation for 140hours.

Effect of the coil diameter on the thermal performance of the slinky heat exchangers

The thermal characteristics of the ground coupled slinky heat exchangers at different coil diameters of 0.6, 0.8 and 1.0m were determined by running the transient state simulation for the similar environmental conditions used in section 3.1. The dimensions of the domain and the heat exchangers are shown in Table 2.

From the prediction, it was found that there was no significant difference of the specific heat extraction of the slinky heat exchanger at different coil diameters of 0.6, 0.8 and 1.0m. This might be because the horizontal-coupled straight part of the heat exchanger occupied over 50% of the total length in the simulated slinky heat exchanger. Thus by changing the diameters of the coil, there was no obvious differences of the specific heat extraction. However, the larger the diameter of the coil, the higher the heat extraction per meter length of soil. The variation of the heat extraction per meter length of the soil occupied by heat exchangers at different coil diameters are illustrated Figure 15. At the beginning, the heat extraction per meter length of soil was around 91W/m for the coil diameter of 1m, approximately 82W/m for the coil diameter of 0.8m and 72W/m for the coil diameter of 0.6m. It decreased to approximately 30W/m with a coil diameter of 1m, around 27W/m with a coil diameter of 0.8m and 24.6W/m with a coil diameter of 0.6m at 140 hours.

Effect of coil central interval distance on the thermal performance of the slinky heat exchangers

The thermal characteristics of the ground coupled slinky heat exchangers at different coil central interval distance of 1.2, 1.6, 2.0 and 3.0m were predicted using the similar environmental conditions used in sections 3.1 and 3.

Conclusions

The thermal performance of horizontal-coupled slinky GSHP was determined both experimentally and numerically in a UK climate. During a two month period of monitoring the performance of the GSHP system, the COP of the GSHP decreased with the running time.

The average COP of the horizontal coupled GSHP was 2.5. To increase the COP of the system, it required longer heat exchanger and a larger land area for heat exchanger installation.

The thermal performance of a portion of slinky GSHP at different coil diameters and different coil central interval distances was predicted using a validated 3D model. The performance of slinky heat exchangers was also compared with that of straight heat exchangers. The specific heat extraction for both the straight and slinky heat exchangers initially at 46W/m decreased with increasing system running time but at different rates. After running the systems for 140hours, the specific heat extraction of the straight pipe would be 3.5W/m higher than that of the slinky pipe. However, the heat extraction per meter length of soil for the slinky heat exchanger was significantly higher than of the straight system. There was no significant difference in the specific heat extraction of the slinky heat exchanger at different coil diameters of 0.6, 0.8 and 1.0m. However, the larger the diameter of coil, the higher the heat extraction per meter length of soil. The specific heat extraction also increased with the increase in the coil central interval distances. On the other hand, the heat extraction per meter length of soil for a slinky decreased with the increase of coil central interval distance.

  months from 6th Nov 2009 to 31st Dec 2009 are shown in Figures 3. The

  measured

  The ground temperature distribution was analysed in more detail for two different days-7th Nov and 28th Dec 2009. It was sunny on 7th Nov 2009, and a partly cloudy day and the ground was covered by snow on 28th Dec 2009. The measured solar radiation intensity, ambient air temperature and wind speed for the two days are shown in Figure 4(a). The measured solar radiation intensity on 7th Nov. 2009 was significantly higher than that on 28th Dec. 2009. The highest solar radiation intensity was around 350W/m 2 on 7th Nov 2009 and approximately 150W/m 2 on 28th Dec 2009. The measured ambient temperature was also higher on 7th Nov 2009. The measured average wind speed was around 1 and 2m/s on 7th Nov and 28th Dec 2009, respectively. The various ground surface temperatures with ambient air temperature are shown in Figure 4(b), and the ground temperature distribution above the heat exchanger and without heat exchanger installed below on 7th Nov 2009 is illustrated in Figure 5. From Figures 4(b) and 5,

Figure 6

 6 Figure 6 illustrates ground temperature distribution above the heat exchanger and ground temperature distribution in the reference hole on 28th Dec 2009. From Figures 4(b) and 6, it

  )

Figure 10 .

 10 Due to a large number of meshes required to model the 3D transient heat transfer, only half of the slinky heat exchanger was applied in the model and a symmetric surface was defined at the vertical central surface as shown in Figure 10(b).

2 .

 2 The variation of the specific heat extraction for the slinky heat exchangers and the variation of the heat extraction per meter length of the soil occupied by heat exchangers at different coil central interval distances are illustrated Figures16 and 17, respectively. From the prediction, it was observed that the specific heat extraction increased with the increase in the coil central interval distance. At the beginning, the specific heat extraction was 44.6W/m for a slinky with coil central interval distance of 3m, 44.2W/m for that at 2m, 44.0W/m at 1.6m and 43.8W/m at 1.2m. At 140hours, the specific heat extraction decreased to 14.9W/m for a slinky with coil central interval distance at 3m, 14.0W/m for that at 2m, 13.3W/m at 1.6m and 11.8W/m at 1.2m. However, the heat extraction per meter length of soil for a slinky decreased with the increase of coil central interval distance. At the beginning, the heat extraction per meter length of soil was 158.4W/m for a slinky with coil central interval distance of 1.2m, 130.5W/m for that at 1.6m, 113.7W/m at 2.0m and 91.3W/m at 3.0m. At 140hours, the heat extraction per meter length of soil decreased to 42.5W/m for a slinky with coil central interval distance of 1.2m, 39.3W/m for that at 1.6m, 36.1W/m for that at 2.0m and 30.5W/m at 3.0m

Figure 1

 1 Figure 1 Image of the horizontal-coupled ground source heat exchanger at Talbot Cottage, Drayton St Leonard site, Oxfordshire, UK

Figure(s)Figure 2 Figure 3

 23 Figure(s)

Figure 9 - 1 Figure 10 Figure 13 Figure 14 Figure 15 Figure 16 Figure 17

 91101314151617 Figure 9 Predicted and measured soil temperatures at selected times on 07th Nov. 2009

Table 1

 1 Measured thermophysical properties of the undisturbed soil samples from the site

	Sample 1 (Depth at 1m)	Sample 2 (Depth at 0.6m)	Sample 3 (Depth at 0.2m)	Average

Table(s)

Acknowledgements:

This work has been supported by the Natural Environment Research Council through the funding of the project. Also, the authors would like to thank Mr Alan Johns for allowing us to access the Ground Source Heat Exchanger used in this work.