Mathematical model of compact type evaporator

In this paper, development of the mathematical model for evaporator used in heat pump circuits is covered, with focus on air dehumidification application. Main target of this ad-hoc numerical model is to simulate heat and mass transfer in evaporator for prescribed inlet conditions and different geometrical parameters. Simplified 2D mathematical model is developed in MATLAB SW. Solvers for multiple heat and mass transfer problems plate surface temperature, condensate film temperature, local heat and mass transfer coefficients, refrigerant temperature distribution, humid air enthalpy change are included as subprocedures of this model. An automatic procedure of data transfer is developed in order to use results of MATLAB model in more complex simulation within commercial CFD code. In the end, Proper Orthogonal Decomposition (POD) method is introduced and implemented into MATLAB model.


Introduction
An increasing effort to decrease the cost of any device requires more advanced methods for precise simulation and modeling of physical phenomena. Heat and mass transfer in heat exchangers/evaporator with complex geometry is still a complicated problem, because an analytical solution is not often available, so numerical simulation along with measurement is the only option that can be used for heat duty estimation and geometry optimization. In this paper, focus is on modeling heat and mass transfer inside evaporators for air dehumidification. Selected evaporator for this application is classified as a cross flow compact type heat exchanger (or finned tube heat exchanger) -tube bundle (single tube with multiple passes) with common fins (every rectangular fin is in contact with all passes of the tube), using R134a as a refrigerant. An evaporator ("Inside Heat exchanger" in Fig. 1.) is placed in heat pump circuit that can be often found in many air dehumidification devices. Other parts of heat pump circuits are: condenser, compressor and capillary tube/reduction valve.

Mathematical model
The mathematical model is developed in MATLAB SW, because of user-friendly environment, simple implementation of matrix operations and ability to use advanced optimization tools for model validation and development. The model contains simple menu interface for user convenience and tool for data transfer into commercial CFD code. Optimization of code has been performed so model can be used for optimization of heat exchanger geometry.

Model properties
The mathematical model simulates heat transfer between humid air and coolant/refrigerant flowing through tubes. Next, simulates mass transfer on outer side -water vapor condensation from air into water film on fins and tubes, and refrigerant evaporation on internal (tube) side. The model predicts locations, where condensation occurs, based on local flow parameters and local fin/tube temperature. Model is parametric in inlet conditions (air/refrigerant temperature and mass flow, humidity, pressure) and also in heat exchanger geometrical parameters. A semi-automatic procedure, that will manage data transfer of results from the MATLAB model into the simulation in commercial CFD package, using scripts and macros is also developed and implemented. A very important attribute of the model is its dimension -in how many axis parameter/property distribution is solved. Choice of model dimension and its complexity affects precision, computation time and computational resources. To obtain a balance between model precision, computation time and computation requirements, 2D model is used. 2 dimensions (coordinate x and y) where flow parameters are solved correspond to 2 edges of a fin as shown in Fig. 2.

Coordinate system and flow orientation
Heat exchanger with arbitrary dimensions and its coordinate system is in Fig. 3. (x -length, y -height, zwidth)

Fig. 3. Coordinate systemorientation
Air flow is always orientated in x direction (from left to right), and refrigerant flow is oriented in z direction -"inlet tube section" is placed at upper right corner of fin and "outlet tube section" at upper left corner (only even number of tube columns with staggered pattern is supported in this version).

Simplifying assumptions
Modeling heat and mass transfer is complicated geometry is a complex problem, so couple simplifications are assumed: [2] 1) The fin is flat. (corrugated fins are usually present, the effect of corrugation in included in AreaFx, AreaFy parameters -area increase factors) 2) For the first guess, refrigerant temperature is constant (equal to boiling temperature) through whole evaporator. During simulation, refrigerant temperature distribution is developing.
3) Condensate from tubes is dropping down, not flowing down the fin.

4)
Temperature of fin is constant through its thickness.

5) No heat flow from fin edges.
6) All fields of flow parameters (temperature, velocity, humidity …) are constant in z direction in every slit. 7) Steady-state flow.
8) Perfect contact between outer tube surface and fin.

Governing equations
Differential equations describing heat and mass transfer problems are listed below. Equations (1) and (2) describe the general balance of mass for air and film flow. Ad hoc equation (3), (4) and (5) describe conservation of enthalpy in 1 slot (space between 2 neighboring fins). Discretization of the equation (6) is used for a surface temperature analysis.
▪ Continuity equation for water film is ▪ Enthalpy conservation equation for humid air is ▪ Enthalpy conservation equation for refrigerant/coolant can be written as

Meshing
Automatic mesh generation (and tube bundle generation) on fin is included in the model. (only parameter of the mesh is Telem -#of elements per tube cross-section). Number of elements in x direction -N and y direction -M is solved automatically for given Telem parameter by sub-optimization (minimization) of tube cross section area error -.  For numerical simulation, it is necessary to divide tube into elements/segments in which temperature/enthalpy increments are evaluated. In this model, tube element is defined by tube segment between 2 neighboring fins.

Solver sub-procedures
MATLAB model itself is composed of multiple subprocedures, which solve various heat and mass transfer problems. [2] ▪ Surface temperature solver ▪ Refrigerant temperature solver ▪ "i + j" cycles to solve air and film local parameters (Air/film enthalpy solver) ▪ Film thickness solver ▪ Film temperature solver

Fig. 6. Solver schema
The iterative solution procedure is obvious from the diagram above. First, the initial distribution of air temperature and condensate flow field is solved from guessed refrigerant temperature distribution. Based on new air temperature, a new refrigerant temperature distribution is solved, from which new fin surface temperature is evaluated. New refrigerant and surface temperature are then used to solve for new air temperature and condensate flow and iteration loop is completed. A number of iterations required to obtain stable solution depends on inlet conditions, but varies from 50 to 150.

Film thickness solver
Because the focus is on evaporators for air dehumidification, low condensate flow is expected, thus, small thickness of a film is assumed. Velocity gradient for a simple 1D case of fluid film on vertical wall can be then derived as [3] [4] Ordinary differential equation above is integrated and boundary conditions applied, velocity profile is obtained.
Where range for z coordinate is 0 ≤ ≤ , , condensate film density, -mixture density and is shear stress at interface between mixture and film. Local mean velocity in condensate film is obtained by integration over local film thickness Local mass flow of condensate film is given by ̇, , = ̅ , , , ∆ . (11)

Fig. 7. Viscosity effect in condensate film
These equations are used to solve local film thickness -, , that is necessary for solving heat transfer through film.

Film Temperature solver
Local temperature of film on tubes and fins is very important parameter for heat transfer and condensation and it is determined from heat balance between humid air, film and fin/tube surface. [3] [6] Heat balance between air, film and surface can be written as ′ ( − ) = ̇′ ′ ℎ 23 + ( − ).
Increment of film mass flow Combining equations above we obtain Equation (15) is solved iteratively From equation (16), mean film temperature is solved after couple iterations. This equation is also used as a part of solution stability monitoring -the first indication of solution divergence is unstable film temperature, if this fact is detected, simulation is terminated.

Fin surface temperature solver
After considering simplifications, problem of fin surface temperature can be then solved separately (based on informations from air and refrigerant side) as 2D heat conduction problem in rectangular domain. To ensure unique solution, boundary conditions are requiredrefrigerant temperatures at fin/tube intersections are used. As a volumetric heat loads, local heat flow from gas (or film) into fin is used. [7] [2] Heat flow (flux) vector is (17) Fourier's law in matrix form can be written as = − .
(21) Equations from weak formulation can be then rewritten as a simple matrix equation = + = .
(22) This method is then applied on heat exchanger geometry for arbitrary inlet conditions, resulting surface temperature is shown in Fig. 8.

Fig. 8. Fin surface temperature
Because plate/fin surface temperature is solved every iteration, code optimization is required to decrease overall computation time even for fine meshes. Temperature in mesh nodes can solved from equation (22) by = − (23) Matrix inverse operation (equation (23)) is very demanding and imprecise for large matrices. The dimension of stiffness matrix for a coarse mesh is typically 7300x7300, which takes 22 seconds on average (4 core CPU, 3,3GHz). This is not suitable for simulations that usually take 150 iterations to converge. Matrix inverse can be replaced by mldivide operation implemented in MATLAB. Another improvement is using algorithms for sparse matrices (CRS formcompressed row storage) and LU factorization.
[ , ] = ( ) Equation (23) Plate/fin surface also needs to be divided into sub-regions (areas), from where all the heat flow coming from the air is assigned to specific tube cross section. This allows evaluating enthalpy/temperature increment of tube pass. In general, these sub-regions can be traced by finding locations, where local heat flux inside the fin is 0 or very small. A boundary of these sub-regions are very complicated (Fig. 9); so an additional simplification is required.

Fig. 9. Complex fin sub-regions
Simplified sub-regions of fin are shown in Fig. 10.

Fig. 10. Simplified fin sub-regions
This simplification also brings new problem -for a case with big temperature gradients between 2 neighboring tubes, some of the heat that in reality "flows" into tube 1 is within this simplification associated with tube 2. This causes an error in heat balance (difference of enthalpy on air side is not equal to enthalpy on the refrigerant side).
To prevent this situation, Redistribution Correction Solver is added. RCS eliminates this error by transferring additional heat from 1 sub-region to the other, based on estimation from heat flux integral on its boundary.

Superposition of heat and mass transfer
Theoretical heat transfer coefficient on fin and tube is known, but effect of tube presence on heat transfer coefficient on fins (and fin presence to heat transfer coefficient on tubes), needs to be included. If we assume that evaporator is only composed of tubes, total heat duty without condensation (mass transfer) can be written symbolically as [2] [8] where α tube is heat transfer coefficient on tubes.
Similarly, if evaporator (heat exchanger) is made only from fins, total heat duty is where ̃ is heat transfer coefficient on fins/plates. Superposing fins and tubes, total heat duty without condensation (mass transfer) is symbolically: = Where superposition (effectiveness) coefficients , are derived from experimental data by optimizationfitting the model results to results from experiment. These coefficients describe influence of temperature and velocity field of one case to the other. Similar hypothesis can be made with condensate mass flow, so another 2 effectiveness coefficients are introduced. Local heat and mass transfer coefficients are obtained from correlations for the case of flow in small slot (as a limit of concentric annular duct) and flow over staggered tube bundle. These coefficients are then used to solve overall heat transfer coefficient according to Fig. 11.   Fig. 11. Overall HTC

POD and RBF method
Another improvement of data transfer process is done by incorporating the discrete version of POD (Proper Orthogonal Decomposition) method -Method of Snapshots. This method allows to simplify a large amount of data into a smaller package and also to reconstruct solution for given (directly unsolved) inlet parameters very quickly. Input for POD method is set of direct MATLAB solutions for different parameters -inlet conditions.
Where is arbitrary points in n-D parameter space and is a point for which direct simulation is performed. Coefficients and are trained with data from direct solutions by MATLAB optimization method fminsearch. Target (fitness) function is defined as where sample error function is evaluated as ∆( ) = (ℎ( ) − ) 2 (33) Total approximation error function can be then written as (34) Necessary number of modes required to "precise" approximation is determined from quantity called level of correlation and it is defined as Where N is a number of direct simulations and n is a number of modes, choose so > 0.99. Implementation of POD and RBF method into the mathematical model is covered more deeply in [2].

Conclusion
The main subject of this paper is the heat and mass transfer modeling inside heat exchangers (evaporator) used in air dehumidification applications. Theoretical and