Optimal Control of Thermal Pollution Emitted by Power Plants

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


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.

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 f x; t ð Þ 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 U opt so that the temperature in Ω OBS will be as close as possible to a desired value T d . T d can be the temperature favorable to the survival of the ecosystem. The geometric domain is illustrated by Figure 1.

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.

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 T f . 0, the temperature dynamic in ΩÂ0, T f ½ is described by the equation T x; t ð Þ represents the fluid temperature at position x ∈ Ω and time t ∈ 0, T f ½. k stands for the thermal diffusion coefficient. u x; t ð Þ is the fluid velocity inducing the advection process. The velocity is obtained by solving the Navier-Stokes system, described below. ψ 1 x ð Þ and ψ 2 x ð Þ are, respectively, Ω 1 and Ω 2 characteristic functions. They allow to localize the source term f t ð Þ and control U t ð Þ, respectively, in the subdomains Ω 1 and Ω 2 . The source term f t ð Þ is given, while the control U t ð Þ must be computed as a solution of an optimal control problem, described in the sequel. T in x; t ð Þ 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. Function T 0 x ð Þ represents the distribution of the temperature at the initial time:

Velocity and pressure
The fluid velocity u is obtained by solving, in ΩÂ0, T f ½, the incompressible Navier-Stokes system: where p x; t ð Þ is the water pressure; ν . 0 is the kinematic viscosity; f 1 t ð Þ and f 2 t ð Þ are, respectively, the velocity sources in Ω 1 and Ω 2 . At the inlet, the velocity is known and given by a function u in x; t ð Þ. It is written On the impermeable boundary, it is assumed that the velocity is equal to zeros due to the viscosity: On the outflow boundary, the pressure is equal to zeros: 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:

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 T d . This optimal control must be the minimum of the cost function: where U d 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 T d .

Cost functional
The aim is to find an optimal control U minimizing the cost function: where U ad ¼ L 2 Ω 2 ð Þ is the admissible function space. By considering the symmetric, continuous, coercive bilinear form for all s, v ∈ U ad , and the linear bounded functional the cost function is writteñ We are in the framework of Theorem 16.1 in [12] that establishes the existence and uniqueness of solution to the minimization problem.

Directional derivative
First, for a fixed h ∈ L 2 ð0, T f ; U ad ½ Þ, two function sequences are considered, for all λ . 0.

Sequences convergence
The difference between equations satisfied by T λ and T results to the following one: If u ∈ L ∞ ΩÂ ð 0, T f ½Þ, this previous system admits a unique weak solution satisfying with C . 0 [13]. The functional spaces are defined by H ¼ L 2 Ω ð Þ and It can be deduced from the preceding inequality that

Directional derivative computation
A direct computation gives us By dividing this last equality by λ, it becomes where w u ¼ w λ λ is solution of the equation: By passing to the limit λ ! 0, the directional derivative is writteñ 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.

Explicit gradient
Equation (25) is multiplied by the adjoint functionp and integrated over Integrations by parts lead to where By using the initial and boundary conditions of w u , the terms I 1 and I 2 become From the condition u ¼ 0 on Γ N , it becomes Hence, we assume thatp is solution of the adjoint problem: Consequently, it becomes I 1 ¼ 0, Using this above equality in relation (26), we obtaiñ A change of variables p x; t ð Þ ¼p x; T f À t À Á is made where p is the solution of The gradient becomes This gradient allows to solve the minimization problem (12). The gradient descent algorithm is used to compute the optimal control.

Iterative algorithm
First, the Navier-Stokes system is solved on Ω Â 0; T f Â Ã to obtain the fluid velocity. Secondly, the optimal control U t ð Þ, t ∈ 0; T f Â Ã , 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: U 0 t ð Þ, Maximal number of iterations: m max , Tolerance: tol. This algorithm is summarized by Figure 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: H being the maximal mesh element diameter.

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:

Navier-Stokes system
The Navier-Stokes system is solved by means of a P 1 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 * u * À u n Δt ¼ À u n :∇ ð Þu n þ νΔu n þf n,

Computation of the pressure p nþ1
3. Computation of the velocity u nþ1 where u n , p n , andf n, n ∈ N * , are, respectively, the approximated velocity, pressure, and source term at the n th 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 http://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.

Numerical results
This section presents numerical tests to illustrate the validity of our approach. The river parameters are listed in Table 1.
The initial temperature is always constant and equal to T 0 ¼ 30°C. The initial velocity is given by u 0 ¼ 0 m:s À1 . The temperature at the inlet boundary is set to T in ¼ 30°C, while the velocity profile is described by the parabolic function: u max is the maximal value of the velocity. At the outflow boundary, mixed boundary conditions are used with α T ¼ 10 À8 and α u ¼ 10 À8 . The velocity source at the discharge is given by with v max . 0. For the optimal control, the target temperature in the observation area is equal to T d ¼ 30°C, and the target control is of U d ¼ 0°C:s À1 . The costefficiency ratio of the objective functional is defined by β ¼ 1. The time step is set to Δt ¼ 0:1 s. 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.

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.

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.

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:s À1 . 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:s À1 , thus obtaining a cost function J T; U ð Þ¼1:37274 * 10 À2 (Figure 9). The optimal solution is obtained after 10 iterations, for an optimal control rate U t ð Þ of order À10 À2 illustrated by Figure 12.
In Figure 10, the stopping criteria |∇J U m ð Þ|=|∇J U 0 À Á | 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:56428 * 10 À3 after 10 iterations. From this observation, it can be deduced that the control sequence U m converges to the optimal control rate. Figure 11 compares the temperature evolution in Ω OBS for the thermal plume dispersion simulation U ¼ 0°C:s À1 , the initial step U ¼ À1°C:s À1 , 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:s À1 , the temperature is lower than 30°C and reaches a minimum of 29.87. For the optimal   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:s À1 . In red, temperature evolution is corresponding to the optimal control. 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 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 .

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.