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: 19 November 2018 Reviewed: 17 July 2019 Published: 13 December 2019

DOI: 10.5772/intechopen.88646

From the Edited Volume

Numerical Modeling and Computer Simulation

Edited by Dragan M. Cvetković and Gunvant A. Birajdar

Chapter metrics overview

1,060 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 x and 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 ΓIN is the entering border, ΓOUT is the outflow boundary, and ΓN is the impermeable part. Ω contains three subdomains Ω1, Ω2, and ΩOBS. Ω1 stands for the industrial plants, where the source of pollution modeled by fxt is defined. In Ω2, cold water is injected at a rate U in order to control the temperature in ΩOBS. The objective consists in finding the optimal rate Uopt so that the temperature in ΩOBS will be as close as possible to a desired value Td. Td can 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


Txt represents the fluid temperature at position xΩ and time t]0,Tf[. k stands for the thermal diffusion coefficient. uxt is the fluid velocity inducing the advection process. The velocity is obtained by solving the Navier-Stokes system, described below. ψ1x and ψ2x are, respectively, Ω1 and Ω2 characteristic functions. They allow to localize the source term ft and control Ut, respectively, in the subdomains Ω1 and Ω2. The source term ft is given, while the control Ut must be computed as a solution of an optimal control problem, described in the sequel. Tinxt is the temperature distribution in the inlet border:


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


where the vector n defined 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>0 is 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 T0x represents the distribution of the temperature at the initial time:

T0=T0 on Ω¯.E6

2.2.2 Velocity and pressure

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


where pxt is the water pressure; ν>0 is the kinematic viscosity; f1t and f2t are, respectively, the velocity sources in Ω1 and Ω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 U at which the freshwater is introduced, such that the temperature in ΩOBS must be as closed as possible to a prescribed threshold denoted Td. This optimal control must be the minimum of the cost function:


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


3. Optimal control

3.1 Cost functional

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


where Uad=L2Ω2 is 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 T results 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 I1 and I2 become

  • I1=ΩwuTfp˜Tfdx,

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

    From the condition u=0 on Γ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˜xTft is 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 Ω¯×0Tf to 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 1 discontinuous 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

H being the maximal mesh element diameter.

4.2 Adjoint equation

As the state equation, the adjoint problem is solved by using a method of 1 discontinuous 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 P1 Lagrange 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 nth time 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:


umax is the maximal value of the velocity. At the outflow boundary, mixed boundary conditions are used with αT=108 and α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 Ω1 by 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 f is applied on Ω1. The control on Ω2 is initialized to U=1°C.s1. The velocity sources on Ω1 and Ω2 are, 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 Ut of order −102 illustrated 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˜U0 in terms of the number of iterations are reported. At initial step, the stopping criteria is equal to 1 and then decreases to reach 9.56428103 after 10 iterations. From this observation, it can be deduced that the control sequence Um converges to the optimal control rate.

Figure 10.

Stopping criteria according to the number of iterations.

Figure 11 compares the temperature evolution in ΩOBS for 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°C and reaches a minimum of 29.87. For the optimal solution, the temperature is >30°C but 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 ΩOBS at 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 Ω1 is considered. In green, simulations are carried out by using the initial control U=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: 19 November 2018 Reviewed: 17 July 2019 Published: 13 December 2019