Open access

Study Regarding Numerical Simulation of Counter Flow Plate Heat Exchanger

Written By

Grigore Roxana, Popa Sorin, Hazi Aneta and Hazi Gheorghe

Submitted: 02 December 2010 Published: 09 September 2011

DOI: 10.5772/24113

From the Edited Volume

Numerical Analysis - Theory and Application

Edited by Jan Awrejcewicz

Chapter metrics overview

4,222 Chapter Downloads

View Full Metrics

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.

Advertisement

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 Ho= 0,003m and the number of channels is Nc =4. Overall heat transfer surface is S=0,218m2.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.

Table 1.

Geometric dimensions of the plate

Figure 1.

Schematic presentation of counter flow plate heat exchanger with eight thermal plates

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.

Table 2.

Measured values

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

V˙i=ΔViτE1

, [m3/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:

Q1=V˙1×ρ1×cp1×(T1,inT1,out)E2

, [W]

Where ρ1,[kg/m3]- density, cp1, [J/kg*oC] – heat capacity at constant pressure, determinate at average temperature of hot water T1,m = (T1,in+T1,out)/2.

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

Q2=V˙2×ρ2×cp2×(T2,inT2,out)E3

,[W]

Where ρ2,[kg/m3]- density, cp2, [J/kg*oC] – heat capacity at constant pressure, determinate at average temperature of hot water T2,m = (T2,in+T2,out)/2.

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

ηt=Q2Q1E4
,

The determinate values are presented in table 3.

Table 3.

Determinate values from experimental results

Advertisement

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:

ρt+ρux+ρvy+ρwz=0E5
,

x-, y-, z- momentum equations:

Figure 2.

model of counterflow plate heat exchanger with 8 thermal plates

ρut+ρuux+ρvuy+ρwuz==ρgxpx+x[2ηux]+y[η(uy+vx)]+z[η(uz+wx)]+Sω+SDRE6
,
ρvt+ρuvx+ρvvy+ρwvz==ρgypy+y[2ηvy]+x[η(uy+vx)]+z[η(vz+wy)]+Sω+SDRE7
,
ρwt+ρuwx+ρvwy+ρwwz==ρgzpz+z[2ηwz]+x[η(uz+wx)]+y[η(vz+wy)]+Sω+SDRE8
.

The two source terms in the momentum equations, Sω and SDR, are for rotating coordinates and distributed resistances, respectively.

The distributed resistance term can be written in general as:

SDR=(Ki+fd)ρVi22CηViE9
,

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:

Sω=2ρωi×Viρωi×ωi×riE10
,

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):

ρcpTt+ρcpuTx+ρcpvTy+ρcpwTz=x[UTx]+y[UTy]+z[UTz]+qVE11

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

Table 4 presents the variable of the equations:

Table 4.

Variables of the governing 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 106÷108) 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:

ρt+ρux+ρvy+ρwz=0E12
,

Momentum equations:

ρUt+ρUUx+ρVUy+ρWUz=ρgxPx+x[2ηUxρuu]+y[η(Uy+Vx)ρuv]+z[η(Uz+Wxρuv)]+Sω+SDRE13
,
ρVt+ρUVx+ρVVy+ρWVz=ρgyPy+y[2ηVyρvv]+x[η(Uy+Vxρuv)]+z[η(Vz+Wyρvw)]+Sω+SDRE14
ρWt+ρUWx+ρVWy+ρWWz=ρgzPz+z[2ηWzρww]+x[η(Uz+Wxρuw)]+y[η(Vz+Wyρvw)]+Sω+SDRE15

Energy equation:

ρcpTt+ρcpUTx+ρcpVTy+ρcpWTz=x[kTxρcPuT']+y[kTyρcpuT']+z[kTzρcpwT']+qVE16

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:

ηt=ρuu2Ux=ρuvUy+Vx=ρvwVz+Wy=...E17
,
kt=ρcputTx=ρcpvtTy=ρcpwtTzE18
.

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

ρt+ρux+ρvy+ρwz=0E19
ρUt+ρUUx+ρVUy+ρWUz=ρgxPx+x[2(η+ηt)Ux]+y[(η+ηt)(Uy+Vx)]+z[(η+ηt)(Uz+Wx)]+Sω+SDRE20
ρVt+ρUVx+ρVVy+ρWVz=ρgyPy+y[2(η+ηt)Vy]+x[(η+ηt)(Uy+Vx)]+z[(η+ηt)(Vz+Wy)]+Sω+SDRE21
ρWt+ρUWx+ρVWy+ρWWz=ρgzPz+z[2(η+ηt)Wz]+x[(η+ηt)(Uz+Wx)]+y[(η+ηt)(Vz+Wy)]+Sω+SDRE22
ρcpTt+ρcpUTx+ρcpVTy+ρcpWTz=x[(k+kt)Tx]+y[(k+kt)Ty]+z[(k+kt)Tz]+qVE23

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

ηt=CηρK2εE24
,
kt=ηtcpPrtE25
,

Where Prt - turbulent Prandtl number and Cη – empirical constant. The transport equations for K and ε are derived using momentum equations:

ρKt+ρUKx+ρVKy+ρWKz=z[(ηTσk+η)Kz]+x[(η+ηTσk)(Kx)]+y[(η+ηTσk)(Ky)]ρε+ηT[2(Ux)2+2(Vy)2+2(Wz)2+(Uy+Vx)2+(Uz+Wx)2+(Vz+Wy)2]E26
ρεt+ρUεx+ρVεy+ρWεz=z[(ηTσε+η)εz]+x[(η+ηTσε)(εx)]+y[(η+ηTσε)(εy)]ρC2ε2K+C1ηtεK[2(Ux)2+2(Vy)2+2(Wz)2++(Uy+Vx)2+(Uz+Wx)2+(Vz+Wy)2]E27

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

Table 5.

Used constants in model

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:

ρUϕx+ρVϕy+ρWϕz=ρUSϕSE28
,

Where s – stream-wise coordinate, Us – 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

Figure 3.

Analysis scheme

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

Figure 4.

Mesh

Table 6.

Boundary conditions

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:

  1. Solve x momentum equation;

  2. Solve y momentum equation;

  3. Solve z momentum equation;

  4. Solve pressure equation and velocities;

  5. Solve energy equation;

  6. Solve turbulent kinetic energy equation;

  7. Solve turbulent energy dissipation equation;

  8. 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).

Figure 5.

Distribution of nodal temperature

Table 7.

Obtained values for average temperatures at the outlet

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.

Figure 6.

Convergence of the nodal temperature after 90 iterations

Advertisement

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:

Q=U×S×LMTDE29

, [W]

Where U, [W/m2K] - 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):

LMTD=ΔTmaxΔTminln(ΔTmaxΔTmin)E30

,[K]

Where ΔTmax = max(ΔT1, ΔT2), ΔTmin =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 Tm1 and Tm2 and approximation of plates temperatures Tp =(Tm1 +Tm2/2.

Figure 7.

The hot and cold fluid temperature distributions in the counter flow heat exchanger

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

dh=4×l×Ho2×(l+Ho)E31

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

wi=V˙iNc×Ho×lE32

,[m/s]

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

4. Calculation of Reynolds numbers.

Rei=wi×dhνiE33
,

Where i=1 for hot water and i=2 for cold water, υi [ m/s2] – 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):

αi=63×Rei0,61×Pri0,4×λi×(ηiηp)0,14E34

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

6. It calculates temperature of the plates:

Tpi=Tm1×α1+Tm2×α2α1+α2E35

, [oC]

And the error ε=|TpTpiTp×100|.

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:

U=11α1+δpλp+1α2E36

, [W/m2K]

8. The maximum number of transfer units:

NTUmax=U×SWminE37
,

Where Wmin = min(W1, W2),[W/K] – minimum heat capacity rate.

W1=V˙1×cp,1×ρ1E38

, [W/K]

W2=V˙2×cp,2×ρ2E39

, [W/K]

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

ε=1eNTUmax×(1WminWmax)1WminWmax×eNTUmax×(1WminWmax)E40
,

Where Wmax =max(Wmin,Wmin).

10. The heat exchanger duty is:

Q=ε×Wmin×(T1,inT2,in)E41

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

T1,out=T1,inQW1E42

,[oC]

T2,out=T2,in+QW2E43

,[oC]

Obtained results are: T1,out = 22,36 oC and T2,out = 19,38 oC.

Advertisement

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.

Table 8.

Results for outlet temperatures

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.

References

  1. 1. AndrianovI. V.AwrejcewiczJ.ManevitchL. I.2004Asymptotical Mechanics of Thin Walled Structures.A Handbook, Springer-Verlag, 354087602Germany
  2. 2. AwrejcewiczJ.KriskoV. A.2004Nonclassical Thermoelastic Problems in Nonlinear Dynamics of Shells, Springer_Verlag, 3-54043-880-7Germany
  3. 3. AwrejcewiczJ.KryskoV. A.KriskoA. V.2007Thermodynamics of Plates and Shells, Springer-Verlag, 978-3-54034-261-8Berlin,Germany
  4. 4. A. Badea, H. Necula, M. Stan, L. Ionescu, P. Blaga, G. Darie2003 Echipamente şi instalaţii termice, Editura Tehnică, ISBN: 973-31-2183-5, Bucureşti, Romania
  5. 5. COSMOS/Flow.2001Technical Reference
  6. 6. GrigoreR.PopaS. 2009 Modeling a Counter-flow Plate Heat Exchanger, Proceedings of 4TH International Conference on Energy and Environment, CIEM 2009, ISSN: 1454-23xx,București, Romania
  7. 7. GrigoreR.PopaS.HaziA.HaziG.2010Study Regarding the Influence of the Plate Heat Exchanger Configuration on Its Performance, WSEAS Transactions on heat and Mass Transfer, 351331421790-5079
  8. 8. GutJ. A. W.PintoJ. M.2003Modeling of plate heat exchangers with generalized configurations, International Journal of Heat and Mass Transfer,, 257125850017-9310
  9. 9. Facultatea de Energetica, Indrumar de schimbatoare de caldura, Bucuresti.2010Available from http://insttermind.3x.ro/Index_files/Iti_indrumar.pdf
  10. 10. SaberM. H.MazaherA. H.2010Simulation and CFD Analysis of heat pipe exchanger using Fluent to increase of the thermal efficiency, Proceedings of 5th IASME/WSEAS International Conference on Continum Mechanics, 978-9-60474-158-8University of Cambridge, UK,February 23-25, 184189

Written By

Grigore Roxana, Popa Sorin, Hazi Aneta and Hazi Gheorghe

Submitted: 02 December 2010 Published: 09 September 2011