Numerical analysis of liquid jet impingement cooling of a thermoelectric generator

Designs of heat exchangers are quite often disconnected to the performance of thermoelectric generators (TEG). In this work, the TEG and the heat exchanger are numerical modelled simultaneously in a computational fluid dynamics (CFD) environment (OpenFOAM) to maximize the output power of the system while minimize the hydraulic power required. A preliminary work was done where the modelling of the heat exchanger, a single laminar slot jet, and the modelling of a 16 element TEG are validated. The considered heat exchanger is a laminar slot jet consists of a linear array of discrete heat sources which accord with the geometry of a thermoelectric generator. The considered 16 element TEG is modelled using the temperature dependent material properties which require a solution of a system of nonlinear differential equations, namely the conservation of energy and the conservation of electric current. The conjugate heat transfer OpenFOAM solver chtMultiRegionFoam is extended by an additional differential equation for the solid region to model the conservation of current. The conservation of energy is expanded by additional source terms based on Peltier/Thomson effect and Joule heat. To simplify the calculation, interface and 1D resistor load boundary conditions are developed and implemented. The heat exchanger and the TEG model, both, are validated by comparisons with measurements, where a good agreement is observed.


Nomenclature
Electric field intensity vector, V/m² Electric current density vector, A/m² R Resistance, Ω Greek Symbols α Seebeck coefficient, V/K σ Electric conductivity, S/m φ Electric scalar potential, V Subscripts C cold H hot n n-type material p p-type material

Introduction
One of the key answers to increase the overall efficiency of industrial processes is the waste heat recovery. The energy emission by industrial processes is recognized as low-grade waste heat in many circumstances. This implies that the temperature is often lower than 600K. It is nearly impossible to recover this kind of waste heat with conventional heat to electricity conversion systems, like steam or gas turbines. A conversion process which can deal with low-grade waste heat is the Organic Rankine Cycle (ORC) [1][2][3], but it is rather complex and the involved turbine contain rotary parts [4] which, in general, needs a high-maintenance effort. The thermoelectric generator (TEG) is a waste heat recovery system, which doesn't contain rotary parts and, thus, needs a low maintenance effort. The today's efficiency of around 5% for the commercial TEGs is a drawback of this technology. Nevertheless, there is an ongoing process in the material research to increase the efficiencies of TEGs and cost optimization of the TEG systems to implement them in various processes optimally with a minimum amortization time. The heat management system around the TEG is a further necessary development, e.g. heat source or cooling, which would decrease the TEG efficiency significantly, if not working in an optimal manner. Even though the TEG cooling system should provide an adequate cooling to achieve a large temperature difference across the TEG and to allow the TEG to work efficiently, the cooling costs should be low to achieve high overall system efficiency. There are diverse approaches to cool a TEG. One approach is the use of the natural convection, where the cooling does not need any power, and, thus, does not influence the system efficiency negatively [5][6][7][8]. The alternative approach is the use of forced convection. It is important to minimize the cooling costs for this approach, in order to preserve high system efficiency. For forced convection, there are mostly two working fluids, which are normally considered for practical applications, i.e. air [9] and water [6,[10][11][12][13]. Though an amount of electric energy needs to be invested to cool by the forced convection, it is generally observed [9][10][11][12][13] that a higher overall TEG system efficiency can still be realized in comparison to the natural convection systems.
A constant temperature boundary condition is mostly assumed for the hot and cold side surfaces of the TEG in the practical development of TEG systems, which represents an ideal condition for an optimal behaviour of the TEG. A one-dimensional model for the prediction the performance and efficiency of a TEG can be derived with this assumption, under the consideration of temperature independent [14] and temperature dependent [15] material properties. However, the constant boundary temperature assumption of these simplified, onedimensional models cannot manage the behaviour of most practical fluid systems that reveal multidimensional, irregular geometries, spatial variations in the flow parameters. Hence, multi-dimensional methods [16,17] are required for being able to perform realistic simulations of practical systems. This is also the goal of the present study, where a prediction instrument for TEG performance with ability to three-dimensional applications is developed.
One of the early investigations on the computational modelling of thermo-electricity is done by Antonova and Looman [16], who modelled the thermo-electric phenomena by using the Finite Element Method (FEM), within the ANSYS software package. Although, the suggested procedure is capable to model in three dimensions, the presented application on the TEG was only two-dimensional. A three-dimensional investigation of a TEG system was presented by Chen et al. [17], where the TEG model was implemented in the commercial computational fluid dynamics software ANSYS Fluent based on the Finite Volume Method (FVM). Nevertheless, the validation of the model in comparison to measurements and other predictions were not very suitable [17]. Baskaya et. al. [18] studied the application of TEG in a rather complex practical case, i.e. to a condensing combi-boiler using the commercial FVM based software ANSYS Fluent. However, their study was indirectly related with the modelling of thermo-electric processes. In the study, six TEGs were considered in the boiler, and a simulation of the hot and cold side of each TEG was performed, without considering thermoelectric effects. Thus, TEGs were modelled simply as heat-conducting parts. However, this work [18] demonstrates the importance of multidimensional effects in TEG applications and is therefore a significant contribution. Rabari et al. [15] presented a one-dimensional, analytical model of a TEG, as well as a two-dimensional computational study of a TEG considering temperature dependent material properties based on an in-house code using the FEM. However, the enhancement of the model to include three-dimensional effects and, beyond this, its coupling to the fluid domain do not seem to be straightforwardly feasible, and were not put as a possible future scope [15], which limits the perspectives for the application of the model in complex practical situations. Hu et al. [19] presented a further, recent study on TEG modelling. They presented a study of different operation conditions for a TEG using the commercial software package COMSOL Multiphysics®. The used TEG model, which is based on the FEM, is a part of the general-purpose software COMSOL Multiphysics®. Thus, three-dimensional applications of the model and its coupling with different physical phenomena are possible.
There are two main goals targeted by the present work. The first is the development of a fully threedimensional TEG model, by properly accounting for all relevant thermo-electric effects, which can readily be coupled with thermo-fluid dynamics of the heating and cooling mediums as well as with the structural mechanics of the solid parts. The second part is a validation study of a laminar slot jet as a preliminary work for the coupling of the thermo-fluid dynamics with the TEG model. The coupling is part of the future work.
A distinguishing feature of the present work is that the model is implemented in a general-purpose, open source C++ toolbox for field operation and manipulation (OpenFOAM) [20], which is principally based on an unstructured FVM that can inherently cope with complex three-dimensional geometries. Coupled with the fluid dynamics and structural mechanics modules of OpenFOAM, a complete multi-physics model of the TEG system can be obtained. Furthermore, the open source platform allows a great amount of freedom in implementing different modelling aspects.
Based on the implemented models for thermoelectricity in OpenFOAM, the before developed solver entitled tegFoam [21] is extended by the ability to handle temperature dependent material properties. The extended solver is again validated by comparisons with the experiments and numerical predictions of Hu et al. [19], who used a commercial package, for a three-dimensional problem. The three-dimensionality results from the spatial distribution of the combined TEG array, whereas each TEG element has a one-dimensional behaviour for its own, in accordance with Hu et al. [19]. Independent from this validation test case, the present model can inherently consider multi-dimensional field distributions within each TEG module.
Furthermore the steady state solver for buoyant, turbulent fluid flow and solid heat conduction with conjugate heat transfer between solid and fluid regions (chtMultiRegionSimpleFoam) is used to model the confined liquid laminar slot jet investigated experimental and numerical by Schafer et. al. [22,23].

Validation of extended solver tegFoam
As already mentioned above, the present study deals with the modelling of the thermo-electric processes in the solid phase of the TEG module. The corresponding governing equations are presented in this section. For a steady state conservation of energy and charge of a solid thermoelectric material with temperature independent material properties, the governing equations can be expressed as [24] (1) (2) with the electric field intensity vector ( ) and electric current density vector ( ) given as Here, the variables T and φ are the resulting fields of temperature and the electric scalar potential, respectively. The coefficients α, k, σ represent the material properties, i.e. the Seebeck coefficient, thermal conductivity, and electric conductivity, respectively.
The differential governing equations are discretized by the FVM by a 2nd order accurate central differencing scheme. The discretized equations (Eq. 1, 2) are solved subsequently, within an iterative procedure based on successive substitution. In the present applications, equidistant, rectangular grids are used. For problems, where the material properties exhibit substantial local variations, the details of the interpolation procedure used for the calculation of diffusion coefficients on cell faces play an important role for the accuracy for a given grid. In the tegFoam solver, the diffusion coefficients are interpolated to the cell face centers by linear interpolation. Higher accuracy could be obtained, if the harmonic mean, which was suggested by Patankar [25], were used. This would enable, in comparison to linear interpolation, the same level of accuracy with a coarser mesh, or a better accuracy with the same mesh. This scheme is not implemented in the tegFoam solver yet. However, grid independence studies have always been performed to ensure a minimal discretization error, while using the linear interpolation scheme.
The applied boundary conditions shall also be discussed in this section. The formulation and implementation of the boundary conditions of the energy conservation equation (Eq. 1) in terms of temperature are quite straightforward to formulate and apply. Principally, boundary conditions of all three kinds can be used. In a coupled treatment of the TEG with the heat exchangers, the temperature distributions on the TEG boundary surfaces will directly result from the solution of the conjugate heat transfer problem, which may also indicate a spatially non-uniform distribution of temperature on these surfaces. Note that the p-type and n-type legs of the TEG are not thermally connected and are subject to hot-side and cold-side temperatures as boundary conditions of their own, which may also differ from one leg to the other.
The formulations of the boundary and interface conditions for the conservation of charge (Eq. 2) are not very straightforward and deserve special attention. They formulate the electric connection between the materials, which ensures that the electric circuit is closed. A special situation arises here, since the legs are connected electrically, which needs to be described by an adequate formulation. The applied formulations for the boundary and interface conditions are described by Pfeiffelmann et al. [21] directly in terms of the discretized equations.

Preliminary studyslot jet
In the second part of the present study a laminar slot jet will be considered, which is a preliminary work for coupling the tegFoam solver with thermo-fluid dynamics. The continuity equation, the Navier-Stokes equations and the energy equation are solved for laminar, incompressible flow of a Newtonian fluid with constant material properties. For the prevailing small dimensions, turbulent flow is unlikely and assumption of laminar flow is realistic [26][27][28], in contrast to the majority of applications in forced convection [29,36]. Buoyancy is neglected. The radiative heat transfer [37] is also omitted currently, as the temperatures are rather low.
For the computational modelling, the steady state solver chtMultiRegionSimpleFoam and the transient solver chtMultiRegionFoam, is used out of the basic software platform OpenFoam [20]. For the velocitypressure coupling, the SIMPLEC [38] and PISO [39] schemes are used for steady-state and unsteady formulations, respectively. For unsteady calculations, a first-order accurate backward differencing scheme is used for the time discretization, whereas the time stepsize is selected to assure the cell Courant number to be smaller than unity. For the discretization of the convective terms, a second-order accurate upwind scheme is used. In all cases grid independence studies are performed and grid independence is ensured.

Validation of extended solver tegFoam
As a validation test case, a 16 element Bi2Te3 TEG module is considered, which was experimentally and computationally investigated by Hu et al. [19]. This module consists of 8 p-type and 8 n-type semiconductor legs as illustrated in Figure 1. This module was measured and modelled for three cases, where temperatures of 403 K, 453 K and 503 K are applied to the hot surface, while keeping the cold surface temperature constantly at 303K. In [19], temperature dependent material properties (λ, α, σ) are considered by means of the following cubic polynomial fits. These polynomials are also adopted in the present study.
p-type Bi2Te3 n-type Bi2Te3 The presently predicted electrical power and efficiency as a function of current are compared to the measurements [19] and the COMSOL Multiphysics® predictions [19] in Figure 2 and Figure 3, respectively. It can be observed that the present predictions and the numerical results of Hu et al. [19] for the power show a quite good agreement (Fig. 2). For efficiency (Fig. 3), the present predictions and the numerical results of Hu et al. [19] show a very good agreement. However, both predictions show deviations from the measurements. As Hu et al. [19] argued, this might be caused by errors in predicting the heat flux by neglecting the parasitic radiative heat flow and imperfect soldering.

Preliminary study -Slot jet
The confined laminar slot jet investigated experimentally and numerically by Schafer et al. [22,23] is considered in the present study as a preliminary work for coupling the tegFoam solver with thermo-fluid dynamics.
In the symmetric slot jet, water flows through a 12.7 mm x 3.18 mm (Lh x B) slot nozzle and impinges on a bottom plate containing five square heat sources (Lh x Lh = 12.7 mm x 12.7 mm) in a 15.87 mm wide channel.
The first heat source was centered directly below the nozzle exit, and the remaining heat sources are placed downstream  The Nusselt number NuL=hLh/λ is averaged over each discrete heat source, where h=q/(Th-To) is the heat transfer coefficient, q the local heat flux, Te the temperature of the fluid at the nozzle exit and Th the local temperature on the heater surface. The other walls are assumed as adiabatic. For the outlet a constant pressure boundary condition and for the inlet a fully developed velocity profile of a duct flow is assumed.
The presently predicted average Nusselt number for the central heater with H/B of 1.5 and 1 is shown in

Conclusions
The new solver for OpenFOAM, tegFoam, was extended to model thermoelectric effects with temperature dependent material properties. The solver was validated for 16 element TEG modules. The results obtained by the developed solver is observed to show a good agreement for current and voltage and a fairly good agreement for heat flow with the experiments [19] and an overall good agreement with the predictions of the commercial finite element solver COMSOL Multiphysics® [19]. It is expected that the accuracy will be increased by the implementation of radiation [37], as well as the harmonic interpolation scheme [25] in the new solver, which will be performed in the future work.
Additionally, numerical simulations have been performed to determine flow and heat transfer of a confined laminar slot jet driven by liquids. The prediction of OpenFoam solver chtMultiregionFoam shows a good agreement with the experiments performed by Schafer et al. [22].
The further modelling step will be the combined modelling of the cold side laminar slot jet together with the TEG module.