Available online at www.sciencedirect.com
ScienceDirect Materials Today: Proceedings 8 (2019) 632–641
www.materialstoday.com/proceedings
ECT_2017
Description of a 3D transient model to predict the performance of an experimental Thermoelectric Generator under varying inlet gas conditions J.E. Jimenez Aispuroa,b*, J-P. Bédécarratsa, S. Gibouta, D. Champierb
a
Universite De Pau & Pays Adour, Laboratoire De Thermique Energétique Et Procédés –Ipra, Ea1932, Bâtiment D’alembert, Rue Jules Ferry, Bp 7511, 64075, Pau Cedex, France. B Universite De Pau & Pays Adour, Laboratoire Des Sciences De L'ingénieur Appliquées A La Mécanique Et Au Génie Electrique – Ipra, Ea4581, Avenue De L'université ,64012, Pau , France.
Abstract
The main challenge of working with thermoelectric generators lies in their low efficiency, so finding the optimal design to achieve the best thermoelectric generator (TEG) performance under the working conditions becomes a strong problematic. In this study, a numerical transient simulation tool was developed using different modelling techniques to achieve the performance improvement of a TEG. This TEG is composed of thermoelectric modules wrapped between two heat exchangers. A three-dimensional heat transfer transient model was first developed for the heat exchanger using the finite volume method and then coupled with a thermoelectric generation model. Materials properties vary as a function of temperature and the Thomson effect is taken into account to be able to get numerical results close to reality. Thermoelectric modules properties used in this model are those given by the manufacturer. The originality of this study is that it takes into account unsteady-state conditions, necessary for simulating a multitude of processes like the automobile (engine and exhaust systems), combustion processes, turbines working conditions and energy production. The use of this numerical tool can help maximize the TEG performance and possibly reduce material cost. In order to validate the model, the TEG has been built with the capacity of using different module quantities for different configurations. Experiments have been conducted at different inlet gas conditions (mass flow rate and inlet temperature). A comparison of experimental and predicted values of temperatures and output electrical power has been done in steady state regime. Using the three dimensional model for the heat exchanger, a good precision of the thermoelectric power generation model is obtained. © 2019 Elsevier Ltd. All rights reserved. Selection and/or Peer-review under responsibility of 15th European Conference on Thermoelectrics.
Keywords: Thermoelectricity ; Thermoelectric generator ; Thermoelectric modules ; TEG ; Transient model. *corresponding:
[email protected] 2214-7853 © 2019 Elsevier Ltd. All rights reserved. Selection and/or Peer-review under responsibility of 15th European Conference on Thermoelectrics.
J.E. Jiménez Aispuro et al./ Materials Today: Proceedings 8 (2019) 632–641
Nomenclature C f I ṁ N Nu Pe Pr R Re Relec S t T V
specific heat capacity, J.kg-1.K-1 Fanning friction coefficient electric current, A mass flow rate, kg.s-1 number of pairs TE, Nusselt number. electric power, W Prandtl number, thermal resistance, K.W-1 Reynolds number, overall electrical resistance, Ω surface, m2 time, s temperatures, K or ° C volume, m3
633
Greek symbols Seebeck coefficient, V.K-1 thermal conductivity, W.m-1.K-1 dynamic viscosity of the fluid, kg.m-1.s-1 density, kg .m-3 Thomson coefficient, V.K-1 thermal flux, W Subscripts C cold side (cold side fluid in the absence of precision) co couple TE H hot side (hot side gas in the absence of precision)
α λ ρ τ
Abbreviation FVM Finite Volume Method TE Thermoelectric TEG Thermoelectric Generator TSI Thomson Seebeck Integrated
1. Introduction Around 60% of primary energy is devoted to electricity production [1]. Nevertheless, this energy conversion is not perfect and heat losses can be high. These losses go up to 70% in certain cases [2] and decreases up to 35% in optimum scenarios [3]. Thermoelectric (TE) technology is viewed as an effective and environmentally friendly solution for exploiting waste heat in order to directly transform it into electrical energy [4]. Thermoelectric generators (TEG) have several notable advantages like having a long useful life: they do not count with moving parts which means less maintenance required and they have a great scalability, meaning that TE technology can be applied to any size of heat source. Nevertheless, TE materials have a slow technology progression and a low energy conversion efficiency rate. For this reason, engineers have to optimize the TEGs to produce the most quantity of electric power. To do so, numerical models are a great tool which allows the engineers to represent the behavior of the TEG, and then, to find the best working conditions, quantity and size of TE modules or geometries for their facilities. Energy-harvesting systems which recover waste heat energy to convert it into electricity with the use of TEGs are being studied and introduced. M. Remeli et al. [5] made an experimental set-up of a heat pipe thermoelectric generator (HP-TEG) with 8 bismuth telluride TE modules arranged in series to the direction of air flow. The maximum power obtained in this experiment was around 7 W when the temperature gradient across the TEG surfaces was at 368 K with an air velocity of 1.1 m.s-1. T.J. Kim et J. Kim [6] performed a study to assess the potential of energy recovery from the exhaust gas of a passenger car engine under various drive cycles using a TEG. Their TEG system consisted of 40 thermoelectric modules coupled to a diesel engine, where the corresponding exhaust gas temperature varied from 700 to 870 K. They found that TEG modules are more efficient at low-speeds and high-load engine operations should be selected to maximize energy recovery during the real-world driving conditions. B. Orr et al. [7] proposed a system design using a TEG with 8 TE modules in series and heat pipes to recover heat from the exhaust gases of a car engine. The maximum power obtained in their experiments was around 38 W with a TEG efficiency of 2.26%, resulting in a potential 1.57% reduction in fuel consumption. N. Kempf et Y. Zhang [8] studied a TEG prototype tested on a diesel engine with semi-Heusler nanostructured high-temperature thermoelectric modules developed for the recovery of heat from automotive waste. The results indicated that the
634
J.E. Jiménez Aispuro et al./ Materials Today: Proceedings 8 (2019) 632–641
material, the pressure drop and the heat transfer efficiency of the heat exchanger have a strong influence on the increase in fuel consumption. G-Y Huang et al. [9] showed that the position of the TE modules affects the generated electrical power and that to maximize it, the internal flow uniformity of the velocity profile must also be taken into account. Other research teams prefer to build a test bench in order to validate a numerical model, allowing them to use it with different working conditions out of their experimental ranges. S. Lan et al. [10] developed a dynamic TEG model for vehicle applications where temperature dependence of material properties and thermal inertia are considered. The simulation results showed that a TEG with 20 thermoelectric modules can produce an electric power between 170 and 224W. Other numerical simulations in transient and steady-state regimes were performed by Li et al. [11] to investigate the thermal and electrical performances of a TEG module placed between a hot block heated by an electrical heater, and a water-cooled block. In their calculations using ANSYS software, the thermal conductivity, resistivity and Seebeck coefficient were temperature dependent. Their results were in agreement with the experimental data and showed that radiation from the hot source does not influence the TEG performances and can be neglected. Favarel C. et al [12] developed a Thomson Seebeck Integrated (TSI) 2D steady-state model and designed a thermoelectric generator to validate the numerical code and verify certain results concerning the optimization of the generator. They sought to evaluate the influence of using different quantities of thermoelectric modules in different configurations and found an optimal occupancy rate (total surface occupied by TE modules on the heat exchanger) for their TEG by reducing the number of modules. Considering the thermal dependence of thermoelectric properties, F. Meng et al. [13] have established a new model of thermoelectric generators. The results indicate that the temperature dependence of the thermoelectric properties has a strong effect on TEG power, efficiency and optimum variables. The majority of developed models and prototypes have been tested under steady-state conditions. However, many systems operate under large variable loads, resulting in significant variation in TEG performance. This paper will describe a simulation tool for steady but also unsteady state conditions, where the heat transfer model has a 3 dimensional geometry and takes into account the temperature dependence of the materials properties and all the TE effects for electric power generation for the thermoelectric model. 2. Experimental setup A thermoelectric generator (figure 1) was built in order to compare and validate the numeric model. This thermoelectric generator is installed in a test bench where the inlet temperature and mass flow rate of the cold and hot sources are controlled. The TEG consists of a double plate-fin heat exchanger in aluminum for the hot part, eight made-to-size exchangers in aluminum for the cold part and 16 thermoelectric modules (TEHP1-12656-0.3, [14]) from Thermonamic made of Bismuth telluride (Bi2Te3).
Figure 1. Experimental setup
J.E. Jiménez Aispuro et al./ Materials Today: Proceedings 8 (2019) 632–641
635
Both the hot and cold heat exchangers help to wrap the TE modules and provide the temperature difference for each of its contact faces. As working fluid for the hot and cold side, hot air and cooled ethanol were used, respectively. The TEG is 0.125 m wide, 0.138 m high and 0.3 m long, giving a maximal capacity of 16 TE modules (8 on the upper side and 8 on the lower side). The test bench is presented in form of diagram (figure 2) where the conforming parts of each loop of this experimental setup can be seen. For the cold loop shown in blue, the cryothermostat is a refrigerating machine (including a pump) that regulates a constant cold side temperature through the cold exchangers, with a working range that varies from 0°C to 20°C. The mass flow rate is measured by an electromagnetic flowmeter and the valves regulate the same mass flow rate for each cold exchanger. For the hot loop shown in red, the air is supplied by a fan with a working range that varies from 0,0055 kg.s-1 to 0,0333 kg.s-1. The mass flow rate is roughly adjusted with a flap box and then finely by acting on the speed of the fan. Consequently, the air passes through the heating resistances and lastly enters the hot heat exchanger. The last part (in green) corresponds to the TE modules with a maximum power point tracking DC/DC converter that adapts and maximizes the electric power obtained in order to stock it in a battery. The main conforming parts of the TEG are highlighted with a dotted line.
Figure 2. Schematic diagram of the experimental setup
Temperatures are measured with an accuracy of +/- 1.5 K using K-type thermocouples connected to a data acquisition system (Agilent 34970A). Sixteen measurements of the hot wall temperature and sixteen cold wall temperature measurements (eight cold exchangers with two measuring points on each) are available (figure 3).
Figure 3. Location of the thermocouples in the exchangers
636
J.E. Jiménez Aispuro et al./ Materials Today: Proceedings 8 (2019) 632–641
Electrical power (voltage and current) is measured using the Agilent central unit (34970A) for each module. Regarding the experiments, tests were conducted for different inlet mass flow rates: 0.0117 kg.s-1, 0.0195 kg.s-1 and 0.0297 kg.s-1, respectively, with an inlet hot gas temperature of 207 °C. The temperature of the cooling fluid was fixed at 20°C. 16 thermoelectric modules were used applying a pressure of 600 kPa. The experimental results will be shown and compared in a following section. 3. Numerical Model The numerical model developed in this study takes into account the coupling between the heat transfer equations and the thermoelectric power generation ones. This model was developed in Matlab, where an implemented toolbox called Coder was used to transform into C some time-consuming calculation part, dynamically linking a subroutine executed in the MATLAB environment. The choice to develop our own code instead of using a commercial one allows us to finely control the equations and hypotheses used. This will also allow us to link this code to the optimization tools developed in the laboratory. The hypotheses taken into account and the main equations used in each model will be mentioned. Finally, a description of the global model including the coupling of these two models is done. 3.1 Heat transfer model To simplify the numerical heat transfer model, the following hypotheses were made: Fluid flow is unidirectional; Heat transfer fluid is Newtonian and incompressible (Mach number M<<0.3); The power of the cryothermostat is sufficient to maintain a constant cold side temperature of 20 °C; Due to a good horizontal symmetry (plane x;z), only the upper part of the TEG is simulated (in gray on figure 4); Exchanges by radiation are negligible and are not considered; Density is constant over temperature;
Figure 4. Upper part of the TEG simulated due to symmetry (dark grey) and slices for parallel TE modules (blue circle)
Finite volumes method (FVM) was chosen for the spatial discretization. The starting point of this method is the breakdown of the field in small volume of control (VCs) in which physical values are considered as homogeneous. In order to define all the sub-systems (volume of control), the heat exchanger is cut in slices for each of its axes (figure 4), being M slices in x, N slices in y and K slices for z, which gives (non-constant) thicknesses of dx, dy and dz called spatial discretization step. This mesh is constructed in a manner that virtual volume frontiers coincide to real frontiers between the different materials of the system (metal, air, modules).
J.E. Jiménez Aispuro et al./ Materials Today: Proceedings 8 (2019) 632–641
637
In the system cut into MxNxK subsystems, continuous variables were replaced by discrete quantities which are applied at the centroid of each node. As a consequence, the continuous temperature field T(x,y,z,t) is approximated by a finite number of temperature , , ( ). The conservation of energy principle applied to each control volume then gives:
.C.(Tm,n ,k ).Vm,n ,k
dTm,n ,k (t )
frontiers
dt
m,n,k
(1)
which is a set of × × ordinary differential equations that can be easily solved using a classical algorithm like the first-order Euler method. Each time-evolving heat exchange term ( ) is evaluated according to classical laws: conduction for solid to solid interaction (i.e. in fins for example); convection using a heat transfer coefficient for solid to fluid interaction (interface between flow and metal); advection in fluid flows; eventually coupling with TEG. Thermophysical properties affecting the heat transfer, such as the thermal conductivity λ, specific heat capacity C, density ρ and viscosity of the fluid μair are temperature dependent. The temperature T, the thermal flux and thus thermoelectric power produced Pe are the response variables. The range of values for the controlled variables such as the mass flow rate, the hot and cold side inlet temperatures are chosen according to the limitations of our facilities (ventilator, heating resistances and cryothermostat). The convective exchanges of the hot air into the heat exchanger are determined using a heat transfer coefficient. The experiments being done in fully developed heat transfer conditions where the flow regime varies between transition and turbulent flow, the correlation of Gnielinski [15] is used : (
=
=
.
.
∗ (
) .
.
(
)
.
)
(2)
(3)
where f is the fanning friction coefficient. Pr the Prandtl number can be between 0.5 ≤ Pr ≤ 106 and Re the Reynold number is between 2300 ≤ Re ≤ 5x106. For the control volumes in contact with the thermoelectric modules (upper and lower sides of the heat exchanger), supplementary terms are added to take into account the effect of the modules. The heat absorbed by the modules is detailed hereafter. The temporal term is discretized using the classical forward-Euler method. All relevant values (mainly temperature) are calculated at discrete locations on the meshed geometry (typically the centroid of each volume of control). These values evolve with time because of the heat/mass exchanges through the frontiers of the corresponding volume of control. After initializing the heat transfer model, it starts to calculate the new temperature for each volume of control. These temperatures are calculated in order to predict the behavior of the heat exchanger and to provide the contact temperature facing the TE modules. This information will be used in the thermoelectric generation model to predict the electrical power produce by TE modules. 3.2 Thermoelectric generation model The methodology of this model consists in solving the equations of thermoelectric generation to simulate the behavior of the thermoelectric modules used in the TEG. The flows passing through the module are described as follows:
638
J.E. Jiménez Aispuro et al./ Materials Today: Proceedings 8 (2019) 632–641
1 2
1 2
TC _ co TH _ co Rco
(4)
1 2
1 2
TC _ co TH _ co Rco
(5)
H N . H .I .TH _ co .Relec .I 2 . .(TH _ co TC _ co ).I
C N . C .I .TC _ co .Relec .I 2 . .(TH _ co TC _ co ).I
Where H and C are the hot and cold fluxes passing through the TE modules, N the number of TE couples, αH and αC the Seebeck coefficients for the hot and cold sides of the TE couples respectively, Relec the electric resistance, I the electric current, τ the thomson effect coefficient, TH,co and TC,co the hot and cold temperatures for the TE couples and Rco the thermal resistance of the TE couples. Two equations representing the contact between the exchangers and the TE modules are also necessary:
H C
TH TH _ co RH TC _ co TC RC
(6) (7)
where TH and TC are the hot and cold side temperatures of the exchangers surfaces in contact with the TE modules and RH and RC are the thermal resistances of the hot and cold side respectively. The coupling with the heat transfer model appears here with the contact temperature of the exchangers TH and TC and with the exchanged fluxes H and
C.
As αH, αC, Relec and Rco are dependent of the TE couples temperature, it is necessary to solve these equations iteratively with a fixed convergence condition of 0.01 K for TH,co and TC,co. The electrical power produced by each TE module, taking into account all the TE effects, is calculated by:
Pe H C N . H .TH _ co C .TC _ co .I TH _ co TC _ co .I Relec .I 2
(8)
The values for H, C and Pe are calculated individually for each TE module. The heat flux H through the module is then supplied to the heat transfer model and represents the heat flux extracted by the modules. This extracted flux is a function of the position of the corresponding module. 3.3 Global numerical model of the TEG In order to have the numerical transient model for the TEG, the thermoelectric generation model has to be coupled to the heat transfer model. Once the input parameters are set, the heat transfer model starts and defines the surfaces and volumes of the cells in the mesh. A flow chart showing the procedure for numerical model convergence is shown in Figure 5. The initial working condition of the heat exchanger corresponds to the ambient temperature. When the iteration loop starts, the contact temperature for each TE module is provided as input for the thermoelectric generation model and once it converges, it returns the value of the electric power produced as well as the heat flux absorbed by each module. This absorbed heat flux is provided to the heat transfer model as the interacting flux for the cells conforming the contact surface (surface in contact with the TE modules) and takes part in the calculations for the temperature of this cells which are evolving with time.
J.E. Jiménez Aispuro et al./ Materials Today: Proceedings 8 (2019) 632–641
639
Figure 5. Flow chart describing the process of convergence of the numerical code
Finally, when the last iteration is effectuated, the heat transfer model achieves global convergence and all the thermal and electrical pertinent values for the study can be recovered, such as temperature evolution through the heat exchanger, hot and cold side temperature of TE couples as well as the heat fluxes through the TE modules and the total electric power produced by the TEG. The experimental and numerical results and comparisons for three different inlet gas conditions will be addressed in the next section. 4. Comparison and Results The comparisons made between the experimental setup and the numerical model are for a total of 16 TE modules positioned in the upper and lower sides of the hot heat exchanger for three different mass flows rates. The values of the mass flow rate are limited by the experimental setup (ventilator working range) but are sufficient to have a large range of working conditions to validate the model. The first validation relates to the thermal part. To be sure that the heat transfer model is correct, a comparison between the evolutions of temperatures needs to be done. The experimental temperature values are obtained from the thermocouples placed in precise positions in the walls of the TEG (figure 3). As the cryothermostat power is strong enough to maintain a constant cold side temperature, only the temperatures of the wall in the hot side are shown. For the three mass flow rates, the comparison between the experimental and predicted values of the temperature is shown in Figure 6. The blue lines correspond to the experimental measures and the red lines correspond to the values given by the numeric model. 8 measurements (4 for the upper part and 4 for the lower part) are available for the inner temperatures of the wall. Because of the symmetry of the system between the upper part and the lower part, only 4 temperatures corresponding to the upper part are presented. The temperature evolution predicted by the code is correct even if there is a slight underestimation of the temperatures for the two less significant mass flow rates. This underestimation decreases as the mass flow rate increases. Therefore, it can be said that the numerical code correctly predicts the inner wall temperatures of the TEG for different mass flow rates higher than 0.0117 kg.s-1.
640
J.E. Jiménez Aispuro et al./ Materials Today: Proceedings 8 (2019) 632–641
Figure 6. Profile temperature for the hot side wall (upper part).
The comparison is found in Table 1 and shows the electric power produced for each slice and the total electric power for the three different mass flows rates. Slice 1 corresponds to the inlet of the heat exchanger and slice 4 corresponds to the outlet. The electric power measured in the experimental setup (Pe_exp) is presented in the white column and the electric power predicted by the numerical model (Pe_model) is in the green column. Table 1. Experimental measures vs Predicted values for electric Power (W)
ṁ = 0 .0117 kg.s-1 # Slice 1 2 3 4 Pe total
ṁ = 0.0195 kg.s-1
ṁ = 0.0297 kg.s-1
Pe_exp 4.9 3.5 2.3 1.8
Pe_model 5.1 3.6 2.5 2.1
Pe_exp 8.2 6.8 4.8 4.2
Pe_model 8.4 6.7 4.9 4,2
Pe_exp 11.7 10.2 7.8 7.3
Pe_model 12.5 10.0 8.4 7.4
12.5
13.3 5.7
23.9
24.2 1.3
37
% Error
% Error
38.3 3.4
% Error
A good agreement between the experimental and numerical results is observed with the highest error for the total electric power equals to 6% for the lowest mass flow. The electric power increases with the mass flow rate due to a better heat exchange with the hot air. The convective heat transfer coefficient increases from 34 W.m-2.K-1 to 61 W.m-2.K-1 when the mass the flow rate increases from 0.0117 kg.s-1 to 0.0297 kg.s-1. It can be noticed the decreasing values of the electrical power along the heat exchanger (when increasing the number of slice). This is due to the heat absorption of the first TE modules in the first slice, they have a bigger temperature difference and thus a bigger TE power generation, however this leads to a bigger heat absorption that decreases the temperature difference for the next slice. These results allow us to validate the numerical model in steady-state in accordance to the predicted temperature and the electric power produced for the experimental measurements.
J.E. Jiménez Aispuro et al./ Materials Today: Proceedings 8 (2019) 632–641
641
5. Conclusion and further works In this study, a numerical model for TE power generation was tested and compared to the results of an experimental TEG with 16 TE modules under steady-state conditions. The results of the comparisons give a maximal error of about 6 percent. The good agreements of the results show that the model is highly accurate to predict the behavior of a TEG. The comparison results allow us to validate the numerical model in steady-state. Further works will include comparison and validation of the numerical model in unsteady state conditions and also carry out optimization studies of the TEG in order to produce more electric power. Acknowledgements This research is supported by Conacyt-SENER (Mexico), the Regional Council of New Aquitaine (Conseil Régional Nouvelle Aquitaine - France) and the Conurbation Committee of Pau-Pyrénées (Communauté d'Agglomération Pau-Pyrénées - France). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15]
IEA. Electricity Information. IEA Stat 2017:1–708. doi:http://dx.doi.org.ezproxy.lib.ryerson.ca/10.1787/electricity-2011-en. Gómez A, Dopazo C, Fueyo N. The future of energy in Uzbekistan. Energy 2015;85:329–38. doi:10.1016/j.energy.2015.03.073. Ibrahim TK, Rahman MM, Abdalla AN. Optimum gas turbine configuration for improving the performance of combined cycle power plant. Procedia Eng 2011;15:4216–23. doi:10.1016/j.proeng.2011.08.791. Champier D. Thermoelectric generators: A review of applications. Energy Convers Manag 2017;140:167–81. doi:10.1016/j.enconman.2017.02.070. Remeli MF, Date A, Orr B, Ding LC, Singh B, Affandi NDN, et al. Experimental investigation of combined heat recovery and power generation using a heat pipe assisted thermoelectric generator system. Energy Convers Manag 2016;111:147–57. doi:10.1016/j.enconman.2015.12.032. Kim TY, Kim J. Assessment of the energy recovery potential of a thermoelectric generator system for passenger vehicles under various drive cycles. Energy 2018;143:363–71. doi:10.1016/j.energy.2017.10.137. Orr B, Akbarzadeh A, Lappas P. An exhaust heat recovery system utilising thermoelectric generators and heat pipes. Appl Therm Eng 2017;126:1185–90. doi:10.1016/j.applthermaleng.2016.11.019. Kempf N, Zhang Y. Design and optimization of automotive thermoelectric generators for maximum fuel efficiency improvement. Energy Convers Manag 2016;121:224–31. doi:10.1016/j.enconman.2016.05.035. Huang G-Y, Hsu C-T, Fang C-J, Yao D-J. Optimization of a waste heat recovery system with thermoelectric generators by threedimensional thermal resistance analysis. Energy Convers Manag 2016;126:581–94. doi:10.1016/j.enconman.2016.08.038. Lan S, Yang Z, Chen R, Stobart R. A dynamic model for thermoelectric generator applied to vehicle waste heat recovery. Appl Energy 2018;210:327–38. doi:10.1016/j.apenergy.2017.11.004. Li W, Paul MC, Montecucco A, Siviter J, Knox AR, Sweet T, et al. Multiphysics simulations of thermoelectric generator modules with cold and hot blocks and effects of some factors. Case Stud Therm Eng 2017;10:63–72. doi:10.1016/j.csite.2017.03.005. Favarel C, Bédécarrats JP, Kousksou T, Champier D. Numerical optimization of the occupancy rate of thermoelectric generators to produce the highest electrical power. Energy 2014;68:104–16. doi:10.1016/j.energy.2014.02.030. Chen L. Thermoelectric Properties on the Power and. Int J ENERGY Environ 2014:pp.137-150. Characteristics G. Specification of Thermoelectric Module (TEHP1-12656-0.3) n.d. Gnielinski V. New equations for heat and mass transfer in turbulent pipe and channel flow 1976;16:359–68.