## 1. Introduction

Heat exchangers are equipments commonly found in industrial applications. Virtually not exist almost industrial area which has not a heat exchanger. This is used to exchange heat between two fluids, cooling and heating processes, heat recovery. The most performant, in terms of heat transfer, are the plate heat exchangers. These types of heat exchangers have a lot of advantages, including a high heat exchange area per unit volume and good heat transfer performance.

An important number of numerical studies applying finite element method have been made to research fluid flow and heat transfer into heat exchangers (Gut&Pinto,2003), (Saber&Mazaher Ashtiani,2010), (Awrejcewicz et al., 2007).

The paper presents a theoretical and experimental study on plate heat exchanger. It is performed a numerical simulation of a counter flow plate heat exchanger using finite element method. A 3D model was developed to analyze thermal transfer and fluid flow along the plate heat exchanger, using COSMOS/Flow program. The results are presented graphically and numerically. In parallel, starting from the same input data, it makes thermal calculations for the studied plate heat exchanger. The basic equations are the equation of heat balance for thermal agents and plate heat transfer equation. The calculation is iterative and has certain features related to channel geometry.

Validation of the models presented is made by comparing the measured values obtained on experimental study.

## 2. Presentation of studied plate heat exchangers - experimental results.

The studied heat exchanger is a pack of 8 stainless steel thermal plates with gaskets. These plates are assembled together in cast iron frames and there are chevron type plates. The hot water flows are in one direction in alternating chambers while the cold water flows are in counter flow in the other alternating chambers like in figure1. The number of passes is 1 and the thermal agents are directed into their proper chambers either by a suitable gasket made from ethylene propylene rubber (EPDM). The width of channel between plates is *H*_{o}= 0,003m and the number of channels is *N*_{c} =4. Overall heat transfer surface is *S*=0,218m^{2}.The geometric dimensions of the thermal plate are represented in table 1.

Such a heat exchanger can be used to warm the cold water, considered a secondary thermal agent, with hot water, a primary thermal agent. Figure 1 shows the simplified presentation of a plate heat exchanger with eight plates in counter flow arrangements.

The plate heat exchanger is equipped with instrumentation for measuring pressure and temperature at the entrance and exit of thermal agents. Also it is measured the volume of hot water and cold water volume in adequate testing stand. The measured values are presented in table 2.

The volume flow rate of hot water and the volume flow rate of cold water are calculated with the next formula:

, [m^{3}/s]

Where *i=1* for hot water and *i=2* for cold water.

The heat transfer rate from hot water is calculated with the next equation, for steady state conditions:

, [W]

Where ρ_{1,}[kg/m^{3}]- density, *c*_{p1}*,* [J/kg*^{o}C] – heat capacity at constant pressure, determinate at average temperature of hot water *T*_{1,m}
*= (T*_{1,in}*+T*_{1,out}*)/2*.

The thermal flux received by cold water is given by equation (3), for steady state conditions:

,[W]

Where ρ_{2,}[kg/m^{3}]- density, *c*_{p2}*,* [J/kg*^{o}C] – heat capacity at constant pressure, determinate at average temperature of hot water *T*_{2,m}
*= (T*_{2,in}*+T*_{2,out}*)/2*.

The coefficient of heat retention for studied counter flow plate heat exchanger is defined with next relation:

,The determinate values are presented in table 3.

## 3. Numerical simulation of the counter flow plate heat exchanger using finite element method

3D geometric model of the heat exchangers is created using SolidWorks program. Figure 2 shows the model of plate heat exchanger with eight thermal plates.

Mathematical modelling includes assignation of governing equations. The partial differential equations (pdes) governing fluid flow and heat transfer include the continuity equation, the Navier-Stokes equations and the energy equation. These equations are intimately coupled and non-linear making a general analytic solution almost impossible.

The governing equations for fluid flow and heat can be written as (Grigore&Popa, 2009):

Continuity equation:

,x-, y-, z- momentum equations:

,,.The two source terms in the momentum equations, *S*_{ω} and *S*_{DR}, are for rotating coordinates and distributed resistances, respectively.

The distributed resistance term can be written in general as:

,Where *i* refer to the global coordinate direction (*u, v, w* momentum equation), *f-* friction factor, *d*- hydraulic diameter, *C* – permeability and the other factors are descript in table 4. Note that the *K*-factor term can operate on a single momentum equation at a time because each direction has its own unique *K*-factor. The other two resistance types operate equally on each momentum equation (Grigore&Popa, 2009), (Cosmos/Flow,2001).

The other source term is for rotating flow. This term can be written in general as:

,Where *i* refer to the global coordinate direction, *ω* is the rotational speed and r is the distance from the axis of rotation.

For incompressible and subsonic compressible flow, the energy equation is written in terms of static temperature (Grigore et al., 2010):

The volumetric heat source term from equation (10) is considered zero for this model.

Table 4 presents the variable of the equations:

The equations describe the fluid flow and heat transfer under steady-state conditions for Cartesian geometries. For the turbulent flow, the solution of these equations would require a great deal of finite elements (on the order of 10^{6}÷10^{8}) even for a simple geometry as well as near infinitesimal time steps. In this paper is used COSMOS/Flow program which solves the *time-averaged* governing equations.

The time-averaged equations are obtained by assuming that the dependent variables can be represented as a superposition of a mean value and a fluctuating value, where the fluctuation is about the mean value and a fluctuating value, where the fluctuation is about the mean. For example, the velocity component in y-direction can be written (Cosmos/Flow, 2001), (Grigore et al., 2010):

*V = V + v’*, [m/s]

where *V,* [m/s] – the mean velocity, *v*’, [m/s] – the fluctuation about the mean. This representation is introduced into the governing equations and the equations themselves are averaged over time. If it uses the notation that the uppercase letters represent the mean values and lowercase letters represents fluctuating values, it can be written the governing equations (Cosmos/Flow, 2001), (Grigore et al., 2010):

Continuity equation:

,Momentum equations:

(13) |

(14) |

(15) |

Energy equation:

The averaging process has produced extra terms in the momentum and energy equations. For turbulent flow, equations of continuity, of momentum and energy is a system of 5 equations with 14 unknowns (Cosmos/Flow,2001). To solve, it is used Boussinesq approximation which defines an eddy viscosity and eddy conductivity:

, .These terms imply that the effect of turbulence is isotropic. With these approximations the governing equations become:

(20) |

(21) |

(22) |

In these conditions should be determined in addition to the 5 unknowns only *k*_{t} and *η*_{t}. The program COSMOS Flow uses a model with two equations for their determination (Cosmos/Flow,2001).

Where Pr_{t} - turbulent Prandtl number and C_{η} – empirical constant. The transport equations for K and ε are derived using momentum equations:

(26) |

(27) |

Where σ_{K} and σ_{ε} – turbulent Schmidt numbers, C_{1}, C_{2} – empirical constants. With these 2 equations, there are 9 equations in 9 unknowns: U,V,W, P,T, μ_{t}, k_{t}, K,ε. Table 5 presents the constants associated to model (Cosmos/Flow,2001).

Finite element method is used to *discretize* the flow domain, thereby transforming the governing partial differential equations into a set of algebraic equations whose solution represents an approximation to the exact analytical solution (Awrejcewicz &Krysko, 2003), (Grigore et al., 2010), (Andrianov et al., 2004). A set of simplified hypothesis are introduced:

Hot water and cold water are Newtonian fluids;

No phase change occurs, the fluids are unmixed;

Turbulent flow is fully developed;

Working fluids are incompressible,

Steady state conditions;

Coefficient of heat retention equal with 1.

It is applied the Streamline Upwind Petrov Galerkin (SUPG) method. The method is used directly on the diffusion and source terms and for the advection terms, the streamline upwind method is used with the weighted integral method. These terms are transformed to stream-wise coordinates, like in next expression:

,Where *s* – stream-wise coordinate, *U*_{s} – the velocity component in the stream-wise coordinate direction, *Φ*- transported quantity.

In figure 3 is shown the used analysis scheme.

The disparagement mode (mesh), is very important for the final results. The models are divided in a multitude of little parts with simple geometrical forms, defined as finite parts, and connected in common points called nods, like in figure 4.

The quality of mesh is high (10-node tetrahedral), mesh type is solid mesh, element size – 5,3588 mm, 365377 nodes and 244363 elements.

Boundary conditions from table 6 are proposed.

For incompressible flows, the most robust condition for the pressure equation is to specify a value at the outlet. Since only relative pressures are calculated by COSMOS/Flow, a value of 0 is recommended (Cosmos/Flow,2001). The numerical simulation of turbulent flow is modeled by *k-*
*ε* turbulence model. *ε* represents the turbulent energy dissipation. To

calculate the boundary layer, either “wall functions” are used, overriding the calculation of *k* and *ε* in the wall adjacent nodes, or integration is performed to the surface, using a “low turbulent Reynolds (*low-Re*) *k-ε*” model (Grigore et. al., 2010).

After the analysis was processed it can be visualized the results, under graphical form or numerical value. Because the governing equations are non-linear, they must be solved iteratively. A Picard or successive substitution is used. In this method estimates of the solution variables are substituted in the governing equations. The equations are solved for new values which are the used as the estimates for the next pass. The convergence criterion

is the level at which the specified variable’s residual norm must reach. With each pass, the residuals should become smaller if the solution is converging. The global iterations is shown below:

Solve x momentum equation;

Solve y momentum equation;

Solve z momentum equation;

Solve pressure equation and velocities;

Solve energy equation;

Solve turbulent kinetic energy equation;

Solve turbulent energy dissipation equation;

Check convergenge (GOTO1)

Analysis runs for 100 iterations, in turbulence conditions, for all eight cases. Profiles are obtained for the following parameters: *u, v, w, T, k*, *ε.* In figure 5 is shown the distribution of the nodal temperature, after 50 iterations.

The hot water temperature and the cold water temperature vary along their flow path, even in the case of constant thermal resistance, because of the flow distribution and temperature gradient variations across the plates.

Convergence control of a solution variable is accomplished by reduction the solution progression rate so that the change of divergence is minimized. COSMOS/Flow has the Graphical Convergence Monitor, where are presented the numerical data like the average, the average, minimum, maximum values for each degree of freedom over the completed range of iterations(Grigore et al.,2010).

The convergence of the nodal temperature is shown in figure 6.

The values obtained for average temperature at the outlet of hot water and outlet of cold water are presented in table 7.

## 4. Theoretical analysis of studied plate heat exchanger

The theoretic analysis is based on the equation of heat balance for thermal agents and plate heat transfer equation. The total rate of heat transfer between the hot and cold fluids passing through a plate heat exchanger may be expressed as:

, [W]

Where *U,* [W/m^{2}K] - the overall heat transfer coefficient, *LMTD*, [K] - the log mean temperature difference in K. *U* is dependent upon the heat transfer coefficients in the hot and cold streams. *LMTD* is computed under assumption of counter flow condition with next relation (Badea et al., 2003):

,[K]

Where *ΔT*_{max} = max(ΔT1, ΔT2), *ΔT*_{min} =min (ΔT1, ΔT2); ΔT1, ΔT2 from figure 7. Figure 7 shows the hot and cold fluid temperature distributions in the counter flow heat exchanger. The heat transfer surface aria *S* is represented along the x-axis and the fluid stream temperature along the y-axis.

The boundary conditions are the same like in experimental case and analyze with finite element method. The following simplifications were considered:

steady-state conditions;

no leakage flow;

no phase change;

physical proprieties are constant in the plate heat exchanger;

uniform temperature and uniform fluid distribution;

no heat losses to the surrounding;

efficiency of counter flow plate heat exchanger is considered 1.

The calculation is iterative and has certain features related to channel geometry and in the same time depends on the flow regimes and criterial relations for convection heat transfer coefficients. Next are presented the principal steps:

1. Approximation of the average temperatures of the thermal agents *Tm*_{1} and *Tm*_{2} and approximation of plates temperatures *Tp* =(*Tm*_{1}
*+Tm*_{2}*/2.*

2. Determination of hydraulic diameter characteristic for studied plate heat exchanger configuration, with next formula:

3. Determination of average velocity under channels for both thermal agents:

,[m/s]

Where *i=1* for hot water and *i=2* for cold water.

4. Calculation of Reynolds numbers.

,Where *i=1* for hot water and *i=2* for cold water, υ_{i} [ m/s^{2}] – cinematic viscosity.

5. For turbulent regimes, with next formula is calculated convection heat transfer coefficient, for each thermal agents( hot water and cold water)(Facultatea de Energetica, 2010):

Where *i=1* for hot water and *i=2* for cold water, ɳ_{i} [ m/s^{2}] – dynamic viscosity, λ_{i} [W/mK] – thermal conductivity, Pr- Prandtl number.

6. It calculates temperature of the plates:

, [^{o}C]

And the error

If ε >2%, it goes to point 1 and the calculus begin again. If ε <2% is still calculating.

7. The overall heat transfer coefficient is determinated:

, [W/m^{2}K]

8. The maximum number of transfer units:

,Where W_{min} = min(W_{1}, W_{2}),[W/K] – minimum heat capacity rate.

^{,} [W/K]

^{,} [W/K]

9. Heat transfer effectiveness is definite like actual heat transfer divided by the maximum possible heat transfer.

,Where Wmax =max(W_{min},W_{min}).

10. The heat exchanger duty is:

11. The outlet temperature of the hot water and the outlet temperature of the cold water:

,[^{o}C]

,[^{o}C]

Obtained results are: T_{1,out} = 22,36 ^{o}C and T_{2,out} = 19,38 ^{o}C.

## 4. Conclusions

The paper presents a simplified model for a plate heat exchanger in a counter flow arrangement. It is realized a numerical simulation and it is observed that the model is in concordance with the experimental results and with the results from theoretical analysis.

Numerical simulation of plate heat exchanger using finite element method is very representative, although it is very laborious and consume more resources from a computer (the geometrical model is much complex, the simulation is more difficult), the results are well presented visual, graphic and numerical.

There are small differences between results. The differences appear due to the simplifying assumptions considered and due to presence of fouling on the surface of the plates. Also a relative degree of uncertainty is introduced by the criterial relations used to calculate convection heat transfer coefficients. In the same time, the plate heat exchanger has corrugated plates patern. Numerical simulation cannot reflect the influence of the corrugation angle and corrugation height, but it offers a good understanding of the temperature distribution and fluid flow under turbulent motion.