Open access peer-reviewed chapter - ONLINE FIRST

# Optimal Control of Thermal Pollution Emitted by Power Plants

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

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

DOI: 10.5772/intechopen.88646

## Abstract

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.

### Keywords

• 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.

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

∂Ω=ΓNΓINΓOUT,E1

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

TtvariationkΔTdiffusion+u.Tconvection=fψ1discharges+Uψ2control.E2

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:

T=TinonΓIN×]0,T[.E3

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

kTn=0onΓN×]0,T[,E4

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:

kTn+αTu.nT=0onΓOUT×]0,T[,E5

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.

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:

utνΔu+u.u+p=f1ψ1+f2ψ2,divu=0,E7

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.

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:

JTU=120TfΩOBSTTd2dxdt+β0TfΩ2UUd2dxdt,E12

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

πsv=<TsT0,TvT0>+β<sUd,vUd>E14

for all s,vUad, and the linear bounded functional

Fv=<TdT0,TvT0>,E15

the cost function is written

J˜v=12πvv+TdT02Fv.E16

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

Tλ=TU+λh,wλ=TλTE17

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

V=vH1Ωsuchthatv=0onΓIN.E21

It can be deduced from the preceding inequality that

limλ0wλ=0inL20,TfVC0TfH.E22

#### 3.2.2 Directional derivative computation

A direct computation gives us

J˜U+λhJ˜U=120Tf<Tλ+T2TdψOBSwλ>+<λβ2UUd+λhψ2h>dt.E23

By dividing this last equality by λ, it becomes

J˜U+λhJ˜Uλ=120Tf<Tλ+T2TdψOBSwu>+<β2UUd+λhψ2h>dt,E24

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

J˜'U·h=0Tf<TTdψOBSwu>+<βUUdψ2h>dt.E26

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

0Tf<TTdψOBS,wu>dtE27

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.

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

0Tf<hψΩ2,p˜>dt=0Tf<wutkΔwu+u.wu,p˜>dt.E28

0Tf<hψΩ2,p˜>dt=0Tf<p˜tkΔp˜u.p˜,wu>dt+I1+I2,E29

where

• I1=ΩwuTfp˜Tfwu0p˜0dx,

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

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

• I1=ΩwuTfp˜Tfdx,

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

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

I2=0TfΓINkwunp˜+ΓNΓOUTkp˜nwudt.

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

0Tf<hψΩ2,p˜>dt=0Tf<TTdψOBS,wu>dt.E31

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

J˜U·h=0Tf<p˜+βUUdψΩ2,h>dt,E32

hence

J˜U=p˜+βUUdψΩ2.E33

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

J˜Ut=p¯Tft+βUUdtψΩ2.E35

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.

## 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.

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

uunΔt=un.un+νΔun+f˜n,

1. Computation of the pressure pn+1

Δpn+1=1Δt.u,

1. Computation of the velocity un+1

un+1=uΔtpn+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 athttp://www.freefem.org/ff++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.

ParametersNotationValueUnit
Viscosityν8.84104m2.s1
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:

hxyt=umaxy2ym.s1.E39

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

f1=0vmaxm.s1,E40

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.

### 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.

### 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:

f1=010m.s1,f2=010m.s1.E41

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.

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 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 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.

## 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.

## How to cite and reference

### Cite this chapter Copy to clipboard

Lèye Babacar, Tine Léon Matar and Sy Mamadou (December 13th 2019). Optimal Control of Thermal Pollution Emitted by Power Plants [Online First], IntechOpen, DOI: 10.5772/intechopen.88646. Available from: