Open access peer-reviewed chapter

Optimal Control of Thermal Pollution Emitted by Power Plants

Written By

Lèye Babacar, Tine Léon Matar and Sy Mamadou

Submitted: November 19th, 2018 Reviewed: July 17th, 2019 Published: December 13th, 2019

DOI: 10.5772/intechopen.88646

Chapter metrics overview

881 Chapter Downloads

View Full Metrics


The coastal areas near thermal or nuclear plants are subject to hot water discharges produced by cooling processes. These activities induce an increase of the temperature near the outlet vicinity, which can extend for miles. The temperature variation affects the metabolic rate of organisms and the level of dissolved oxygen. Cooling by cold water from an additional discharge can be considered in order to limit this thermal pollution. This paper present a methodology based on the implementation of a two-dimensional numerical model to study the dynamic of the temperature originated from the industrial discharges. Moreover the optimal injection rate of cold water is sought to keep the water temperature as close as possible to the survival of the ecosystem. Numerical simulations are performed to illustrate the efficiency approach.


  • CFD
  • power plant
  • modeling
  • thermal pollution
  • optimal control

1. Introduction

Thermal pollution is defined as the degradation of water quality by any process that changes ambient water temperature [1]. Coastal areas are often subject to thermal effluents originating from the cooling processes in industrial plants (nuclear reactors, electric power plants, petroleum refineries, steel melting factories, etc.) [2]. The industries collect water from lakes, rivers, or ocean, for cooling purpose, and return it in the environment at a high temperature. The hot water affects aquatic life, causes the substitution of aquatic fauna and flora, increases the mortality of certain species, and has indirect effects including bacterial development. More precisely, increasing the water temperature often increases the susceptibility of organisms to toxic substances (which are undoubtedly present in contaminated water) [3, 4, 5, 6].

Studying the thermal effluents in receiving environments can contribute to efficiently manage the discharges, reducing environmental and economic impacts. Hence the reduction of thermal pollution must be included in the installation of cooling systems. Moreover distance between inlet and outlet must be carefully determined to avoid a decrease of the power plant efficiency.

By the middle of the 1960s, there were many research projects concentrating on thermal discharges, where major publications focus on the environmental impacts of power plant thermal discharges. Early mathematical models took place with works of [7]. The first treatments addressed the equilibrium iso-contour of elevated temperature within the receiving waters. Slightly later more advanced models allow the analysis of thermal plumes across extensive data in relation to seasonal and climate change fluctuations [3, 8, 9, 10, 11].

In this research, thermal pollution due to industrial activities was modeled by a system of partial differential equations, and optimal control is applied to reduce the associated thermal pollution. The location of the understudy area is illustrated by Figure 1.

Figure 1.

Water domainΩ, industrial plantsΩ1, control zoneΩ2, observation zoneΩOBS, impermeable boundaryΓN, inflow boundaryΓIN, and outflow boundaryΓOUT.

The paper is organized as follows. First, the thermal pollution is modeled by a coupling of Navier-Stokes and heat equations. The cold water injection rate is minimum of a cost function, in order to reduce the temperature variation and the energy required to refresh injected water. Afterward, the well-posedness of this problem is investigated. It follows a numerical resolution of the optimal control by means of a gradient descent algorithm. Finally, numerical simulations are performed to illustrate our approach.


2. CFD modeling of the thermal dispersion

2.1 Geometry representation

We are interested in the evolution of the system in space and time. Then, we denote xand t, respectively, as the space and time variables. Ωis the domain occupied by the water. Its boundary is denoted as ∂Ωand is divided into three disjoint subborders. It is written


where ΓINis the entering border, ΓOUTis the outflow boundary, and ΓNis the impermeable part. Ωcontains three subdomains Ω1, Ω2, and ΩOBS. Ω1stands for the industrial plants, where the source of pollution modeled by fxtis defined. In Ω2, cold water is injected at a rate Uin order to control the temperature in ΩOBS. The objective consists in finding the optimal rate Uoptso that the temperature in ΩOBSwill be as close as possible to a desired value Td. Tdcan be the temperature favorable to the survival of the ecosystem. The geometric domain is illustrated by Figure 1.

2.2 Mathematical model

We present the system of partial differential equations representing the evolution of the river parameters (temperature, velocity, and pressure). Then, the cost function to be minimized in order to reduce the thermal pollution is described.

2.2.1 Temperature

By hypothesis, three processes influence the temperature evolution: the thermal conduction, the convection, and the internal reactions. The thermal conduction translates the fact that the heat flux is proportional to the temperature gradient. The convection expresses the temperature transfer by the fluid velocity. The internal reactions are represented by the different sources of temperature and the industrial plant discharges in this situation. By taking into account these processes, for a time Tf>0, the temperature dynamic in Ω×]0,Tf[is described by the equation


Txtrepresents the fluid temperature at position xΩand time t]0,Tf[. kstands for the thermal diffusion coefficient. uxtis the fluid velocity inducing the advection process. The velocity is obtained by solving the Navier-Stokes system, described below. ψ1xand ψ2xare, respectively, Ω1and Ω2characteristic functions. They allow to localize the source term ftand control Ut, respectively, in the subdomains Ω1and Ω2. The source term ftis given, while the control Utmust be computed as a solution of an optimal control problem, described in the sequel. Tinxtis the temperature distribution in the inlet border:


On the impermeable boundary, no heat flux boundary condition is considered:


where the vector ndefined on the boundary constitutes the outward unit normal vector. On the outflux boundary, the heat flux is proportional to the velocity and the temperature:


where αT>0is a constant. The boundary condition allows us, as we will see in the sequel in Subsection 3.2, to obtain an explicit formula for the cost function gradient. These boundary conditions for the temperature are summarized in Figure 2.

Figure 2.

Boundary conditions for temperature.

Function T0xrepresents the distribution of the temperature at the initial time:

T0=T0 on Ω¯.E6

2.2.2 Velocity and pressure

The fluid velocity uis obtained by solving, in Ω×]0,Tf[, the incompressible Navier-Stokes system:


where pxtis the water pressure; ν>0is the kinematic viscosity; f1tand f2tare, respectively, the velocity sources in Ω1and Ω2. At the inlet, the velocity is known and given by a function uinxt. It is written

u=uin on ΓIN×]0,Tf[.E8

On the impermeable boundary, it is assumed that the velocity is equal to zeros due to the viscosity:

u=0 on ΓN×]0,Tf[.E9

On the outflow boundary, the pressure is equal to zeros:

p=0 on ΓOUT×]0,Tf[,E10

The boundary conditions applied to the velocity and the pressure are summarized in Figure 3.

Figure 3.

Boundary conditions for the velocity and pressure.

The system is also completed by the initial condition for the velocity:

u0=u0 on Ω¯.E11

2.2.3 Cost functional

In order to reduce the pollution in an arbitrary area ΩOBS, a freshwater is introduced in the subdomain Ω2. We are seeking the optimal rate Uat which the freshwater is introduced, such that the temperature in ΩOBSmust be as closed as possible to a prescribed threshold denoted Td. This optimal control must be the minimum of the cost function:


where Udis the ideal control rate and β>0is the cost-efficiency ratio. More βis great; more energy must be provided to maintain the temperature in ΩOBSclose to Td.


3. Optimal control

3.1 Cost functional

The aim is to find an optimal control Uminimizing the cost function:


where Uad=L2Ω2is the admissible function space. By considering the symmetric, continuous, coercive bilinear form


for all s,vUad, and the linear bounded functional


the cost function is written


We are in the framework of Theorem 16.1 in [12] that establishes the existence and uniqueness of solution to the minimization problem.

3.2 Directional derivative

First, for a fixed hL20,TfUad, two function sequences


are considered, for all λ>0.

3.2.1 Sequences convergence

The difference between equations satisfied by Tλand Tresults to the following one:

wλtkΔwλ+u.wλ=λhψ2 in Ω×]0,Tf[,wλ=0 on ΓIN×]0,Tf[,kwλn=0 on ΓN×]0,Tf[,E18
kwλn+αTu.nwλ=0 on ΓOUT×]0,Tf[,wλ0=0 on Ω¯.E19

If uLΩ×0,Tf, this previous system admits a unique weak solution satisfying


with C>0[13]. The functional spaces are defined by H=L2Ωand


It can be deduced from the preceding inequality that


3.2.2 Directional derivative computation

A direct computation gives us


By dividing this last equality by λ, it becomes


where wu=wλλis solution of the equation:

wutkΔwu+u.wu=hψ2 in Ω×]0,Tf[,wu=0 on ΓIN×]0,Tf[,kwun=0 on ΓN×]0,Tf[,kwun+αTu.nwu=0 on ΓOUT×]0,Tf[,wu0=0 on Ω¯.E25

By passing to the limit λ0, the directional derivative is written


Unfortunately, this directional derivative does not provide an explicit expression of the gradient. To achieve this, the term


must be written as a scalar product on h. In this scope, the Lagrangian approach is used and consists of solving an adjoint system, stated below.

3.3 Explicit gradient

Equation (25) is multiplied by the adjoint function p˜and integrated over Ω×]0,Tf[. The result is


Integrations by parts lead to



  • I1=ΩwuTfp˜Tfwu0p˜0dx,

  • I2=0Tf∂Ωkp˜n+α˜wukwunp˜dσdt.

    By using the initial and boundary conditions of wu, the terms I1and I2become

  • I1=ΩwuTfp˜Tfdx,

  • I2=0TfΓINkwunp˜+ΓNkp˜n+α˜wu+ΓOUTkp˜nwudt.

    From the condition u=0on ΓN, it becomes


Hence, we assume that p˜is solution of the adjoint problem:

p˜tkΔp˜u.p˜=TTdψOBS in Ω×]0,Tf[,p˜=0 on ΓIN×]0,Tf[,kp˜n=0 on ΓNΓOUT×]0,Tf[,p˜T=0 on Ω¯.E30

Consequently, it becomes I1=0, I2=0, and


Using this above equality in relation (26), we obtain




A change of variables p¯xt=p˜xTftis made where p¯is the solution of

p¯tkΔp¯u.p¯t=ψOBSTTdTft in Ω×]0,Tf[,p¯=0 on ΓIN×]0,Tf[,kp¯n=0 on ΓNΓOUT×]0,Tf[,p¯0=0 on Ω¯.E34

The gradient becomes


This gradient allows to solve the minimization problem (12). The gradient descent algorithm is used to compute the optimal control.

3.4 Iterative algorithm

First, the Navier-Stokes system is solved on Ω¯×0Tfto obtain the fluid velocity. Secondly, the optimal control Ut, t0Tf, is computed by means of a descent algorithm with a fixed step. And finally, this optimal control is used in the state equation to simulate the fluid temperature propagation. The optimal control is the limit of the sequence:


τbeing the step. The algorithm used is described as follows:

Input: Initial control: U0t, Maximal number of iterations: mmax, Tolerance: tol.

This algorithm is summarized by Figure 4.

Figure 4.

Flow chart of the iterative algorithm of the solution to the optimal control problem.

Algorithm 1.

Optimal control algorithm.


4. Numerical scheme

4.1 State equation

The state equation is solved by using a method of 1discontinuous Galerkin in space and implicit finite difference in time. The fluid velocity is very high in relation to its thermal conductivity. To stabilize the induced oscillations, streamline diffusion [12] is introduced in the scheme; hence the solved state equation is as follows:

TtkΔT+u.T+HuuΔuTstreamline diffusion=fψ1+Uψ2,E37

Hbeing the maximal mesh element diameter.

4.2 Adjoint equation

As the state equation, the adjoint problem is solved by using a method of 1discontinuous Galerkin in space and implicit finite difference in time. Streamline diffusion is introduced in the scheme; hence the solved adjoint state is as follows:

p¯tkΔp¯u.p¯+HuuΔup¯streamline diffusiont=ψOBSTTdTft.E38

4.3 Navier-Stokes system

The Navier-Stokes system is solved by means of a P1Lagrange finite element method for the velocity and the pressure. The following algorithm proposed by Chorin [14] is used for the time discretization:

  1. Computation of an intermediate solution u


  2. Computation of the pressure pn+1


  3. Computation of the velocity un+1


where un, pn, and f˜n, nN, are, respectively, the approximated velocity, pressure, and source term at the nthtime step. The mesh is frequently adapted to improve the solution efficiency. For the numerical implementation, the solver of partial differential equations FreeFem++ downloadable at allows to do the domain meshing, the computation, and the post-processing of the solution. The numerical code is run on a computer of characteristics ProBook 250 G2, processor Intel(R) Core(TM) I5-5200U CPU @ 2.20 GHz 2.20 GHz, and RAM memory 8.00 Go.


5. Numerical results

This section presents numerical tests to illustrate the validity of our approach. The river parameters are listed in Table 1.

Thermal diffusionk1.5107m2.s1

Table 1

Physical parameters of the river.

The initial temperature is always constant and equal to T0=30°C. The initial velocity is given by u0=0m.s1. The temperature at the inlet boundary is set to Tin=30°C, while the velocity profile is described by the parabolic function:


umaxis the maximal value of the velocity. At the outflow boundary, mixed boundary conditions are used with αT=108and αu=108. The velocity source at the discharge is given by


with vmax>0. For the optimal control, the target temperature in the observation area is equal to Td=30°C, and the target control is of Ud=0°C.s1. The cost-efficiency ratio of the objective functional is defined by β=1. The time step is set to Δt=0.1s. The stopping criteria tolerance of the iterative algorithm is given by tol=0.02, and the step of the descent gradient algorithm by τ=0.5.

5.1 Thermal pollution simulation

We assume that hot water is discharged in Ω1by power plants. The distribution of the water temperature at different time steps is shown in Figure 5. The flow velocity is presented in Figure 6. It can be observed that flow displaces the thermal plume from the power plants to the right hand side.

Figure 5.

Temperature contours. (a) Temperature contour after 0 second, (b) Temperature contour after 10 second, (c) Temperature contour after 20 second, (d) Temperature contour after 30 second, (e) Temperature contour after 40 second, and (f) Temperature contour after 50 second.

Figure 6.

Velocity fields. (a) Velocity field after 0 second, (b) Velocity field after 10 second, (c) Velocity field after 20 second, (d) Velocity field after 30 second, (e) Velocity field after 40 second, and (f) Velocity field after 50 second.

5.2 Influence of the discharge rates

In Figure 7, we present the influence of different discharge rates on the thermal plume area. For low rates, we observe a high temperature far from the discharge area. However, for high rates, the temperature seems to be high near the emission zone and small far from the discharge zone. In this last case, the polluted area is more extended.

Figure 7.

Temperature contours for different discharge rates. (a) vmax = 20, (b) vmax = 30, (c) vmax = 40, and (d) vmax = 50.

5.3 Heat emission optimal control

Similarly to the first case, a source term fis applied on Ω1. The control on Ω2is initialized to U=1°C.s1. The velocity sources on Ω1and Ω2are, respectively:


The cost function value according to the optimization iteration is represented by Figure 8. At the initial step, we assume that U=1°C.s1, thus obtaining a cost function JTU=1.37274102(Figure 9). The optimal solution is obtained after 10 iterations, for an optimal control rate Utof order −102illustrated by Figure 12.

Figure 8.

Cost functional after each iteration of the steepest descent algorithm.

Figure 9.

Control rate after iterations 0, 4, and 10.

In Figure 10, the stopping criteria J˜Um/J˜U0in terms of the number of iterations are reported. At initial step, the stopping criteria is equal to 1 and then decreases to reach 9.56428103after 10 iterations. From this observation, it can be deduced that the control sequence Umconverges to the optimal control rate.

Figure 10.

Stopping criteria according to the number of iterations.

Figure 11 compares the temperature evolution in ΩOBSfor the thermal plume dispersion simulation U=0°C.s1, the initial step U=1°C.s1, and the optimal case. For the thermal dispersion, it can be noticed that after 8 s the temperature increases to reach a maximum of 30.15°C. For the initial step U=1°C.s1, the temperature is lower than 30°Cand reaches a minimum of 29.87. For the optimal solution, the temperature is >30°Cbut does not exceed 30.08°C. Moreover for all the time, it is closer to the optimal value than for the thermal dispersion case. According to these remarks, it can be concluded that the computed optimal rate allows to maintain the temperature in ΩOBSat a value close to the desired threshold 30°C.

Figure 11.

Time series of water mean temperature in the observation zoneΩOBS. In blue, only hot water discharge inΩ1is considered. In green, simulations are carried out by using the initial controlU=1°C.s1. In red, temperature evolution is corresponding to the optimal control.

Figure 12 illustrates the temperature distributions at times 1, 5, 10, 30, 45, and 59.9 s. A reduction of the thermal pollution is observed, due to the cold water source in Ω2.

Figure 12.

Optimal control: temperature field (in °C) at time intervals of 1 (a), 5 (b), 10 (c), 30 (d), 45 (e), and 59.9 s (f).


6. Conclusion

Numerical models are essential to predict the thermal effluent impacts on natural systems. This work is of particular relevance for the coastal area managements, by contributing to a better understanding of these consequences on coastal dynamics. The cooling water discharge causes an increase of the water temperature. Model simulations show that water dynamic plays a significant role on the temperature dispersion. The optimal control of the model allows to define a strategy to limit this pollution. The simulations show that an injection of freshwater, at an appropriate rate, allows to reduce this pollution and keeps water temperature favorable for ecosystem survival.


  1. 1. Precht H, Christophersen J, Hensel H, Larcher W. Heat Exchange with the Environment. Berlin, Heidelberg: Springer; 1973. pp. 545-564
  2. 2. Issakhov A, Zhandaulet Y. Numerical study of technogenic thermal pollution zones’ formations in the water environment from the activities of the power plant. Environmental Modeling and Assessment. 2019;24:1-16
  3. 3. Abbaspour M, Javid AH, Moghimi P, Kayhan K. Modeling of thermal pollution in coastal area and its economical and environmental assessment. International journal of Environmental Science and Technology. 2005;2(1):13-26
  4. 4. Coulter DP, Sepúlveda MS, Troy CD, Höök TO. Thermal habitat quality of aquatic organisms near power plant discharges: Potential exacerbating effects of climate warming. Fisheries Management and Ecology. 2014;21(3):196-210
  5. 5. Hester ET, Doyle MW. Human impacts to river temperature and their effects on biological processes: A quantitative synthesis 1. JAWRA Journal of the American Water Resources Association. 2011;47(3):571-587
  6. 6. Kelso JRM, Milburn GS. Entrainment and impingement of fish by power plants in the great lakes which use the once-through cooling process. Journal of Great Lakes Research. 1979;5(2):182-194
  7. 7. Edinger JE, Brady DK, Geyer JC. Heat Exchange and Transport in the Environment. Report No. 14. Technical Report. NTIS Issue No. 198314. 1974, p. 137
  8. 8. El-Ghorab EAS. Physical model to investigate the effect of the thermal discharge on the mixing zone (case study: North Giza Power Plant, Egypt). Alexandria Engineering Journal. 2013;52(2):175-185
  9. 9. Hunt CD, Mansfield AD, Mickelson MJ, Albro CS, Geyer WR, Roberts PJW. Plume tracking and dilution of effluent from the Boston sewage outfall. Marine Environmental Research. 2010;70(2):150-161
  10. 10. Kalinowska MB, Rowiński PM. Thermal Pollution in Rivers—Modelling of the Spread of Thermal Plumes. Cham: Springer International Publishing; 2015. pp. 591-613
  11. 11. Issakhov A, Zhandaulet Y. Numerical simulation of thermal pollution zones’ formations in the water environment from the activities of the power plant. Engineering Applications of Computational Fluid Mechanics. 2019;13(01):279-299
  12. 12. Quarteroni A. Numerical Models for Differential Problems. New York, USA: Springer; 1988
  13. 13. Brézis H. Analyse fonctionelle: Collection Mathématiques appliquées pour la matrise. New York, USA: Masson; 1983
  14. 14. Chorin A. Numerical solution of the Navier-Stokes equations. Mathematics of Computation. 1968;22:10

Written By

Lèye Babacar, Tine Léon Matar and Sy Mamadou

Submitted: November 19th, 2018 Reviewed: July 17th, 2019 Published: December 13th, 2019