Physical parameters of the river.
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.
- power plant
- thermal pollution
- optimal control
Thermal pollution is defined as the degradation of water quality by any process that changes ambient water temperature . 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.) . 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 . 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 and , 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 is the entering border, is the outflow boundary, and is the impermeable part. contains three subdomains , , and . stands for the industrial plants, where the source of pollution modeled by is defined. In , cold water is injected at a rate in order to control the temperature in . The objective consists in finding the optimal rate so that the temperature in will be as close as possible to a desired value . 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.
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 , the temperature dynamic in is described by the equation
represents the fluid temperature at position and time . stands for the thermal diffusion coefficient. is the fluid velocity inducing the advection process. The velocity is obtained by solving the Navier-Stokes system, described below. and are, respectively, and characteristic functions. They allow to localize the source term and control , respectively, in the subdomains and . The source term is given, while the control must be computed as a solution of an optimal control problem, described in the sequel. is the temperature distribution in the inlet border:
On the impermeable boundary, no heat flux boundary condition is considered:
where the vector 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 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 represents the distribution of the temperature at the initial time:
2.2.2 Velocity and pressure
The fluid velocity is obtained by solving, in , the incompressible Navier-Stokes system:
where is the water pressure; is the kinematic viscosity; and are, respectively, the velocity sources in and . At the inlet, the velocity is known and given by a function . 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:
2.2.3 Cost functional
In order to reduce the pollution in an arbitrary area , a freshwater is introduced in the subdomain . We are seeking the optimal rate at which the freshwater is introduced, such that the temperature in must be as closed as possible to a prescribed threshold denoted . This optimal control must be the minimum of the cost function:
where is the ideal control rate and is the cost-efficiency ratio. More is great; more energy must be provided to maintain the temperature in close to .
3. Optimal control
3.1 Cost functional
The aim is to find an optimal control minimizing the cost function:
where is the admissible function space. By considering the symmetric, continuous, coercive bilinear form
for all , and the linear bounded functional
the cost function is written
We are in the framework of Theorem 16.1 in  that establishes the existence and uniqueness of solution to the minimization problem.
3.2 Directional derivative
First, for a fixed , two function sequences
are considered, for all .
3.2.1 Sequences convergence
The difference between equations satisfied by and results to the following one:
If , this previous system admits a unique weak solution satisfying
with . The functional spaces are defined by 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 is solution of the equation:
By passing to the limit , 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 . 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 and integrated over . The result is
Integrations by parts lead to
By using the initial and boundary conditions of , the terms and become
From the condition on , it becomes
Hence, we assume that is solution of the adjoint problem:
Consequently, it becomes , , and
Using this above equality in relation (26), we obtain
A change of variables is made where 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.
3.4 Iterative algorithm
First, the Navier-Stokes system is solved on to obtain the fluid velocity. Secondly, the optimal control , , 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:
This algorithm is summarized by Figure 4.
4. Numerical scheme
4.1 State equation
The state equation is solved by using a method of 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  is introduced in the scheme; hence the solved state equation is as follows:
being the maximal mesh element diameter.
4.2 Adjoint equation
As the state equation, the adjoint problem is solved by using a method of 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:
4.3 Navier-Stokes system
The Navier-Stokes system is solved by means of a Lagrange finite element method for the velocity and the pressure. The following algorithm proposed by Chorin  is used for the time discretization:
Computation of an intermediate solution
Computation of the pressure
Computation of the velocity
where , , and , , are, respectively, the approximated velocity, pressure, and source term at the 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.
5. 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 . The initial velocity is given by . The temperature at the inlet boundary is set to , while the velocity profile is described by the parabolic function:
is the maximal value of the velocity. At the outflow boundary, mixed boundary conditions are used with and . The velocity source at the discharge is given by
with . For the optimal control, the target temperature in the observation area is equal to , and the target control is of . The cost-efficiency ratio of the objective functional is defined by . The time step is set to . The stopping criteria tolerance of the iterative algorithm is given by , and the step of the descent gradient algorithm by .
5.1 Thermal pollution simulation
We assume that hot water is discharged in 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.
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 is applied on . The control on is initialized to . The velocity sources on and are, respectively:
The cost function value according to the optimization iteration is represented by Figure 8. At the initial step, we assume that , thus obtaining a cost function (Figure 9). The optimal solution is obtained after 10 iterations, for an optimal control rate of order −10illustrated by Figure 12.
In Figure 10, the stopping criteria in terms of the number of iterations are reported. At initial step, the stopping criteria is equal to 1 and then decreases to reach after 10 iterations. From this observation, it can be deduced that the control sequence converges to the optimal control rate.
Figure 11 compares the temperature evolution in for the thermal plume dispersion simulation , the initial step , 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. For the initial step , the temperature is lower than 30and reaches a minimum of 29.87. For the optimal solution, the temperature is >30but does not exceed 30.08. 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 at a value close to the desired threshold 30.
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 .
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.