Numerical analysis and simulation of the heat recovery from wastewater using heat exchanger

The problem of global warming and the reduction of energy consumption have led to an evolutionary progress of research directed towards finding as many solutions as possible to these environmental issues. Firstly, this paper presents the background information on the role of wastewater as a source of heat for the future. Next, the paper includes the analysis elements that define a system for recovering thermal energy from wastewater. The main objective was to identify the parameters that determine the heat transfer. It has started from a conceptual model of the technological system that involves inputs and outputs characterized by technological, physical-chemical, measurable or imposed properties. In the second part this paper presents a numerical model elaborated for the analysis and simulation of the main physical processes, the mass and heat transfer, which underlie the operation of the heat pipe heat exchangers (HPHE). The numerical simulation of heat and mass transfer in the HPHE is computed by using Delphi 7 solver program. This program contained a series of sub-programs for the meshing of the field occupied by the HPHE, another subprogram for solving the meshing equations and the third for post processing. The design of HPHE is the key to provide a heat exchanger system to work proficient as expected. Finally, the result is used to optimize and improving heat recovery systems of the increasing demand for energy efficiency in residential buildings or industry.


Introduction
One of the secondary energetic resources that on which a greater accent has been put worldwide are wastewaters (liquid waste). Wastewater from all sources except toilets, enters the sewage system at a relatively high temperature and with an exergy content [1,2]. Clearly, the heat content of wastewater is a major loss source, so it must be addressed as a solution to reduce energy consumption in both residential and commercial buildings [3]. The potential for domestic wastewater recovery is quite high, approximately 3.5 kWh of energy / person / day, which could be capitalized and used directly in several areas. [2] The recovery of secondary energetic resources and of liquid waste with high energetic content also led to the development and implementation of various types of heat exchangers and heat recuperators [4,5]. Heat exchangers are widely used in the industry, transportation, constructions, thermal power plants, heating and air conditioning systems, electronic equipment. In all these applications, improving the efficiency of heat recovery can lead to substantial reduction of costs of production and use. The research is oriented in this direction aim in particular at selecting the work fluid with high thermal conductivity and how to direct the flow towards the surfaces destined for the actual heat transfer [6,7].
The current paper presents results of research obtained from simulation and numerical analysis of heat recovery from wastewater using a HPHE.
The paper describes the physical mathematical model, establishing the system of governing equations, initial and limits conditions. It also presents the numerical model associated to the mathematical model and numerical methods of discretization when solving the differential system of equations underlying the mathematical model.
Simulation of thermal transfer is conducted for a liquid-liquid type heat exchanger, where heat pipes operate on the principle of phase change of the work fluid. [8,9] Worldwide there are a series of papers that deal with theoretical and experimental research of the physical processes associated with the various types of heat exchangers, and more recently a great emphasis is placed on the use of numerical models [10][11][12].
Numerical models used for the analysis and optimization of mass and heat transfer processes between work fluids in heat exchangers are less costly than experimental research [13]. But it is important to note that these numerical models are improved based on the results of experimental research [14][15][16][17].

Analysis and modelling 2.1 Description of the experimental equipment
The heat exchanger used for experimental research ( Fig.  1) was equipped with 40 heat pipes from S.C. EnergiQ Indco S.R.L Cluj-Napoca, with the following characteristics: outside diameter 10 mm; length 1390 mm; working fluid: refrigerant occupied approximately 20% of the total volume of the heat pipe at the vacuum at 10 -8 Torr. The heat exchanger was divided in two compartments (Fig. 2) -the lower one destined for the warm wastewater and the upper one destined for the cold water from the network, which is to be preheated. To these are added the outer wall and partition wall between the two compartments. The upper compartment is provided with a horizontal flow divider, thus increasing the time interval in which the main flow of cold water is in contact with the heat pipes. The heat exchanger occupies a range that is approximated by parallelipipedic cells having the dx, dy and dz sides, in number of im x jm x km, where im, jm and km are the number of cells in the direction of coordinate axes Ox, Oy and Oz.
In Fig. 3 is represented a horizontal section parallel with the xOy coordinate plane. The 40 heat pipes are represented into two parallel lines. The wall of the heat pipes is approximated also with parallelipipedic cells, with a volume Vcell = dxdydz. The flow of wastewater and water from the network is treated as a viscous flow, and the heat transfer at the interface between the heat pipes and water is achieved through conduction. Within the domain occupied by water, the heat transfer is achieved through both conduction and convection.
The structure of the water flow in cavities equipped with heat pipes is complex, mainly because of the occurrence of stationary vortexes above or laterally to the pipes, or of stagnation areas downstream of the pipes. These influence the heat transfer; vortexes contribute to the increase in efficiency of the local heat exchange, while the stagnation areas inhibit the heat exchange. [12,18]

Governing equations for the system
In the xOyz Cartesian Coordinate System we have the equations of the numerical model, based on fundamental physical principles, expressed by the relationship of conservations laws for movement quantity -Navier Stokes equations (1-3), mass -continuity equation (4) and heat -Fourier equation (5) [19][20][21]: where: u, v and w are components of the velocity vector; μ is dynamic viscosity; x, y and z are independent variables in the xOyz system; T(x,y,z,t) is temperature; k is the thermal conductivity coefficient; T0 is the reference temperature, in this case the wastewater or cold water input temperature; β is the thermal expansion coefficient; c is the specific heat at constant pressure.
For the control the advence of the free surface of the liquid used the equation of the volume fraction (6).
where: f(x,y,z,t) is the volume fraction function, if in point (x,y,z) there is liquid it is equal to one and if there is no liquid it is equal to zero. [22,23]

Numerical model
During the discretization of the numeric modeling, the differential equations from VOF were used [18,22,23]. A horizontal section of the discretization cells applied to the heat pipes is sketched in Fig. 4. Inside the discretization network, cells colored in blue indicate the pipe's wall.  A parallel section to plane yOz is presented in Fig. 6. It passes through the center of the current P cell and its neighbors, representing the notations and associated grid points; the notations are those used by Patankar [20]. Components of the speed vector are defined in central points of the faces of the current cell, from axis Oy in w and e, and from axis Oz in t and b. It can be observed a cross grid formed by the central points and face centers. The grid points for the v component, by axis Oy, are on the faces parallel with plane yOz. An advantage of the cross grid is the fact that mass flows passing through the faces of the cells can be calculated without interpolation of the speed vector components (Fig. 6).
The parallelepiped centered in point e having the rectangle tcdb as vertical section parallel to yOz, will be the control volume of the v components (i, j, k).

Fig. 6.
Vertical section, parallel with plane yOz, through the center of current cell P and its neighbors.
A horizontal section through the center of the current cell P and its neighbors shows the arrangement of the speed components, u and v, in the cross grid, being represented in Fig. 7. Due to the fact that network cells can contain various materials, it is required the control of cell content, thus an additional variable was introduced in the form of a tridimensional matrix. Next, there was conducted the temporal discretization using an Euler scheme assumed with the time step dt. The time step is variable and is determined in such a way that conditions of numeric stability to be satisfied. It integrates the basic equation system on cells in the network in a way that differs by the usual developments in Taylor series that are applied in the classic finite differences methods [21,24]. The integration of the differential equation with a dependent variable Φ ϵ {u, v, w, T, p} leads to an equation in the form: where: a0=ρ0V0/dt, and V0 and ρ0 are the volume and density of the cell in point P at the previous time step, and dt is the dimension of the time step. Coefficients aE, aW, aT, aB, aS and aN contain the sum of advective and diffusive contributions of the neighboring cells (Fig. 5). The diffusive component that was used has the general form Γf A/Dn where A is the area of the common face between the current cell and its neighbor, and dn is the distance between the centers of the two cells. The exchange coefficient at the interface between cells, Γf, is calculated as an arithmetic mean for density and viscosity, and as a harmonic mean for temperature values in the neighboring central nodes. When calculating the advective component it's used an "upwind" type interpolation which leads to a coefficient in the form of ρA max (0, v) where v is the component of the speed vector that is perpendicular on the interface between cells [19,20].
The system of the mesh equations for the movement quantity conservation (1-3), continuity equation (4) and heat equation (5) are solved in an iterative manner.
The Fourier equation (5) can be solved independently, but the continuity equation (4) and movement quantity equations (1-3) must be solved simultaneously. For the latter, is used the SIMPLE (Semi-Implicit Method for Pressure -Linked Equations) method, which simultaneously solves equations (1-4) using pressure and speed corrections.
Discretization the continuity equation (4) and movement quantity equations (1-3) we obtained (8) equation what is the standard discretization form of the Eq. (2) of movement quantity conservation in the Oy axis direction, similarly obtained for Ox (9) and Oz (10) axis.
Following the discretization of the heat conservation equation is resulted eq. (11):

Stabilization of the speed correction equation
The momentum equations are discretized as follows: Resolution of discretized equations is possible only if the pressure field is given or estimated. If this pressure field is incorrect, the continuity equation cannot be obtained by the speed field given by the momentum equations.
Components of an imperfect speed field (momentum equations with an estimated p * pressure field) are noted with u * , v * and w * [20]. The following discretized equations had as a solution this "*" speed field, as follows:  (13) obtained from the discretization equations of the moment, in which the terms containing the pressure forces, and thus the free terms bs, be and bt, were highlighted, and became  V The field of stellar speeds must satisfy the continuity equation, in this case p * must be modified until it does not satisfy this condition. The correct pressure p is obtained with relation [20]: where p' is the pressure correction. Modification of the pressure value induces a change in the speed components in the form of: where u, v and w are the correct speed components, and u', v' and w' are the speed corrections. Subtracting in turn equations (15) from equations (12) we obtained:  A  p  p  w  a  w  a  w  a  w  a  w  a  w  a  w  a   A  p  p  v  a  v  a  v  a  v  a  v  a  v  a  v  a   A  p  p  u  a  u  a  u  a  u  a  u  a  u  a u a (16) where relations (14) and (15) where taken into account, and: As=dydz, Ae=dxdz, At=dxdy For the calculation of speed corrections obtained with equations (16) we used the SIMPLE procedure. The method is called semi-default because the terms containing speed corrections from the neighboring grid points are omitted, thus the correction speeds became: (17) or and regrouping the terms we have: It should be noted that, although the notations of the coefficients resemble the notations used in the discretization of the heat conservation equation, here they have a completely different meaning (pressure).
In this context we have seen a lot of interests in the literature about the experimental and theoretical investigation of the heat pipes. [25][26][27][28]

Results and discussion
For the analysis of the flow and heat transfer phenomena in the upper compartment was used a spatial grid with steps dx=dy=0.002 m and dz=0.004 m, resulting a number of im x jm x km = 383 x 46 x 78 = 1374204 computational cells. But observing that the vertical median plane -perpendicular on the wall that contains input and output orifices from the upper compartmentparallel with lateral walls, is a symmetry plane, number of computational cells was reduced in half by replacing jm=46 cu jm=23, resulting a number of 687102 cells.
Required temperatures: lower wall of upper compartment tinf=29°C, wall of heat pipes at constant temperature of tpipe=32°C. To take account of heat losses at the interface with the environment, an additional layer of air was added near the lateral faces and the upper face of the compartment, with a constant temperature of tair = 22°C.
Thermophysical characteristics of materials that are part of the upper compartment are shown in Table 1. Values of the parameters used for the initial conditions are shown in Table 2. Two runs were performed and for each of them the initial condition varied: The runs were conducted with a variable time step, determined based on the maximum flow speed at the beginning of each cycle.
The flow divider in the upper compartment has the purpose of increasing the heat exchanger's efficiency, because the main flow passes twice between the condensing sides of the heat pipes [11], thus the cold water will accumulate a larger amount of heat -in a given interval of time. Figure 8 a) and b) presents the position of the free surface in vertical transversal sections close to the input orifice, at 10 mm from the wall with the input orifice, and in the section that contains the first two heat pipes. We can see the partially filled cells and the formation of swirling currents near the walls, centered in points A and B. Also, there are swirling currents formed between pipes and lateral walls. Results obtained from running the program, in which the heat exchange is taken into account, are presented below.
In Figure 9 is presented the distribution of temperature in the central vertical section at t=168 s from the beginning of operation of the heat exchanger and at t'1= 35°C, t'2 = 12°C, tair = 22°C, vw.i. = 0.083 m/s. Near the lateral walls, the upper wall and in the flow divider plate, the temperature is lower, with values between 13 -15°C. Near the inferior wall, in contact with the lower compartment and heat pipes, there is a warm layer with temperatures of 25 -27°C. Inside, there are warmer traces next to the heat pipes. At the left end of the divider plate there is a deformation of these traces, due to a swirling current formed in this phase of operation of the heat exchanger. In Figure 10 is shown the intensity of speeds, in horizontal section made a distance of 2 cm from the inferior wall. We can see the main current between the 2 rows of pipes, with maximum speeds in the center of the section, slowing of the flow towards the lateral wall opposed to the input. Speed intensity is lower between two consecutive pipes, but also near the central area of the lateral walls. Temperature of the preheated water outputting the upper compartment was obtained by a mean of the temperature values in the computational cells, which represent the input sections, namely: where: Toutput is the temperature in the evacuation section; Devac is the domain formed by cells that are part of the evacuation section; t is time that passed since the beginning of operation of the heat exchanger; N is the number of computational cells in the evacuation section.
Ties(t) values were saved every step of time and based on these values were constructed graphs of temperature evolution in the evacuation section. In analogy, was determined the mean pressure in the input sections pinput and in the evacuation sections poutput, in order to determine the pressure drop between the input orifice and the evacuation orifice: In Figure 11 are presented the variation of temperatures of preheated water, at the output orifice from the upper compartment. We can observe that time required to obtain a state of quasi-stationary state depends of the input speed, initial temperature of the cold water and temperature of the environment. In this figure we can observe how the cold water temperature varies depending on time at the output orifice of the upper compartment in the case of the first run.
Due to the heat insulating paint applied on the lateral and upper walls, the heat exchange coefficient between the walls and the environment have a relatively low value of 2 W/m 2 °C. At lower input speeds, the quasistationary phase is obtained later. Similarly, at lower environment temperatures, the time required to obtain the quasi-stationary state increases. a) b) Fig. 11. Variation of temperature of heated water at the output orifice: a) first run; b) second run.
A centralization of the values obtained for temperature during the two runs, is shown in Table 3. Results obtained for the estimated efficiency are presented in this table based on relation [29].  We can observe that a maximum efficiency was obtained for run number one, and the lowest efficiency for run number two. The cause of low efficiency for run number two is high input speed, which reduces the contact time of cold water with the heat pipes.

Conclusions
The program used in the heat and mass transfer simulation was done in the Delphi 7 programming language, which is a variant of the Turbo Pascal language. This program contained a series of subprograms for discretization of the domain occupied by the heat exchanger with heat pipes, another subprogram for solving the discretization equations and a third one for post-processing.
Parameters for the flow and heat transfer analysis were chosen to be as close as possible to the real existing conditions, due to the fact that at an experimental level on the laboratory stand we could not use a high enough flow and due to weather conditions (cold water temperature was higher).
The results, obtained by computerized simulation, revealed that the recovered temperature depends greatly on the flow of wastewater and cold water, which validates the experimental results.
Also, the obtained results revealed that the distribution of speed field inside the upper compartment, between the heat pipes, leads to the formation of swirls, where the heat transfer is much lower, due to the way in which the placement of pipes that are in line.
In the context of the above, we can conclude that heat recovery from wastewater in quasi-stationary regime depends on: the speed of the input cold water, the initial temperature of the cold water and the temperature of the environment, which was confirmed by the conducted experimental research.