**3. Numerical Method**

The numerical analysis is conducted on the physical domain of the thermoelectric generator system for waste heat recovery. The physical domain of the thermoelectric generator system for waste heat recovery is comprised of the heat exchanger, cold fluid channels, thermoelectric modules and fluid domains of the hot gas and water. The coupled numerical approach comprises the CFX and the thermal electric solvers of ANSYS 19.1 commercial software is used to predict the current, power and thermal efficiency of the thermoelectric generator system for waste heat recovery [4]. The continuity, momentum, energy and thermoelectric coupling equations are solved using the coupled numerical approach. The boundary conditions for solving governing equations are shown in Table 1. The experimental voltage conditions with time are used as the high potential voltage conditions for the numerical analysis. The low potential voltage is set to 0 V to ensure the flow of current from the high potential to the low potential of the thermoelectric module. The tetrahedrons mesh structure and fluid domains are used to solve the governing equations because of the complex geometrical structures of the thermoelectric generator system for waste heat recovery with the presence of fins and sharp corners. To verify the convergence of the predicted results, the grid dependency test is carried out for five different grid element numbers. Figure 2 shows the grid dependency test for the simulated power and thermal efficiency of the thermoelectric generator system for waste heat recovery. As shown in Figure 2, the simulated power and thermal efficiency values are converged within ±1% above a grid element number of 11,431,310. Therefore, considering the computational time and the computational cost, the grid element number of 11,431,310 is selected as the final mesh configuration for the thermoelectric generator system for waste heat recovery to predict its current, power and thermal efficiency numerically. In numerical analysis, the density, specific heat, thermal conductivity and dynamic viscosity of the hot gas are set as 1.19 kg/m<sup>3</sup> , 1005 J/k· ◦C, 0.026 W/m·K and 1.8 × 10−<sup>4</sup> kg/m·s, respectively. The density, specific heat, thermal conductivity and dynamic viscosity of the water are set as 997 kg/m<sup>3</sup> , 4182 J/kg· ◦C, 0.607 W/m·K and 8.9 × 10−<sup>4</sup> kg/m·s, respectively. For skutterudite, the density, specific heat and thermal conductivity are used as 7598 kg/m<sup>3</sup> , 350 J/kg· ◦C and 3.4 W/m·K, respectively. In addition, the Seebeck coefficient is set as 142.8 µV/K for p-leg and −183.5 µV/K for n-leg [22]. The continuity, momentum and energy equations are expressed with Equations (2) to (5) [23].

Continuity Equation

$$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathcal{U}) = 0 \tag{2}$$

Momentum Equation

$$\frac{\partial(\rho \mathcal{U})}{\partial t} + \nabla \cdot (\rho \mathcal{U} \times \mathcal{U}) = -\nabla p + \nabla \tau + \mathcal{S}\_M \tag{3}$$

Stress tensor τ is expressed in terms of strain rate as follows:

$$
\tau = \mu \left(\nabla \mathcal{U} + \left(\nabla \mathcal{U}\right)^{T} - \frac{2}{3} \delta \nabla \cdot \mathcal{U}\right) \tag{4}
$$

Energy equation Energy equation

$$\frac{\partial(\rho h)}{\partial t} + \nabla \cdot (\rho \mathcal{U} h) = \nabla \cdot (\lambda \nabla T) + \tau : \nabla \mathcal{U} + \mathcal{S}\_E \tag{5}$$

where ρ is the density (kg/m<sup>3</sup> ), *U* is the average velocity (m/s), ∇ is the nabla operator, *p* is the static pressure (Pa), τ is the stress tensor, *S<sup>M</sup>* is the momentum source, µ is the dynamic viscosity (Pa·s), *h* is the enthalpy (J), λ is the thermal conductivity (W/m·K) and *S<sup>E</sup>* is the energy source. where is the density (kg/m3), is the average velocity (m/s), ∇ is the nabla operator, is the static pressure (Pa), is the stress tensor, ெ is the momentum source, is the dynamic viscosity (Pa∙s), ℎ is the enthalpy (J), is the thermal conductivity (W/m∙K) and ா is the energy source. To deal with the turbulence of the hot gas and cold water, the k-ɛ turbulence model [3] is used

*Symmetry* **2020**, *12*, x FOR PEER REVIEW 5 of 30

To deal with the turbulence of the hot gas and cold water, the k-ε turbulence model [3] is used with Equations (6) and (7): with Equations (6) and (7):

$$\frac{\partial(\rho u\_i k)}{\partial \mathbf{x}\_i} = \frac{\partial}{\partial \mathbf{x}\_i} \left[ \left( v + \frac{v\_t}{\sigma\_k} \right) \frac{\partial k}{\partial \mathbf{x}\_i} \right] + P\_k - \rho \varepsilon \tag{6}$$

$$\frac{\partial(\rho u\_i \varepsilon)}{\partial \mathbf{x}\_i} = \frac{\partial}{\partial \mathbf{x}\_i} \left[ (v + \frac{v\_t}{\sigma\_\varepsilon}) \frac{\partial \varepsilon}{\partial \mathbf{x}\_i} \right] + \rho \mathbf{C}\_1 \mathbf{S} \varepsilon - \rho \mathbf{C}\_2 \frac{\varepsilon^2}{k + \sqrt{v\varepsilon}} \tag{7}$$

where *u<sup>i</sup>* is the velocity component, *k* is the turbulence kinetic energy, *x<sup>i</sup>* is the cartesian co-ordinates, *v<sup>t</sup>* is the turbulent eddy viscosity, σ*<sup>k</sup>* is 1 based on reference [3], σ<sup>ε</sup> is 1.2 based on reference [3], *P<sup>k</sup>* is the shear production of the turbulent kinetic energy, ε is the dissipation rate of the turbulence energy and *S* is the modulus of the mean rate of the strain tensor. where is the velocity component, is the turbulence kinetic energy, is the cartesian coordinates, ௧ is the turbulent eddy viscosity, is 1 based on reference [3], ɛ is 1.2 based on reference [3], is the shear production of the turbulent kinetic energy, ɛ is the dissipation rate of the turbulence energy and is the modulus of the mean rate of the strain tensor. The thermoelectric coupling equations [24] of Equations (8) and (9) are calculated:

The thermoelectric coupling equations [24] of Equations (8) and (9) are calculated:

→

$$
\nabla \cdot (\mathbf{p} \cdot \overrightarrow{J}) - \nabla \cdot (\mathbf{k} \cdot \nabla T) = \overrightarrow{J} \cdot \overrightarrow{E} \tag{8}
$$

$$
\nabla \cdot (\sigma \cdot \mathbf{a} \cdot \nabla T) + \nabla \cdot (\sigma \cdot \nabla \mathbf{Z}) = 0 \tag{9}
$$

where *p* is the Peltier coefficient (V), *J* is the electric current intensity (A/m<sup>2</sup> ), *k* is the thermal conductivity (W/m·K), ∇*T* is the temperature gradient, → *E* is the electric field intensity (V/m), σ is the electrical conductivity (Ω−<sup>1</sup> ·m−<sup>1</sup> ), α is the Seebeck coefficient (V/K), and ∇∅ is the electric potential (J/C). conductivity (W/m·K), ∇ is the temperature gradient, ሬ⃗ is the electric field intensity (V/m), is the electrical conductivity (Ω−1·m−1), is the Seebeck coefficient (V/K), and ∇∅ is the electric potential (J/C).

**Figure 2.** Grid dependency test for the simulated power and thermal efficiency. **Figure 2.** Grid dependency test for the simulated power and thermal efficiency.


**Table 1.** Boundary conditions.
