Sensitivity computations using first and second orders perturbative methods for the advection-diffusion-reaction model of pollutant transport

Sensitivity computations using first and second orders perturbative methods for the advection-diffusion-reaction model of pollutant transport

A. Fraidenraich^I; P. M. Jacovkis^II; F. R. de A. Lima^III

^IDepartamento de Matemática, Facultad de Ingeniería, Universidad de Buenos Aires, Av. Paseo Colón 850, 1063 Buenos Aires. Argentina. E-mail: [e-mail] ^IIDepartamento de Matemática, Facultad de Ingeniería, Universidad de Buenos Aires, Departamento de Computación, Facultad de Ciencias Exactas y Naturales, Universidad de Buenos Aires, Ciudad Universitaria 1428, Buenos Aires. Argentina. E-mail: [e-mail] ^IIIDepartamento de Energia Nuclear, Centro de Tecnologia, Universidade Federal de Pernambuco, Av. Prof. Luiz Freire, 1000, Cidade Universitária, 50740-540 Recife, PE. Brazil. E-mail: [e-mail]

ABSTRACT

This work aims to apply the disturbance theory to accomplish sensitivity computations in problems of pollutant transported in liquid media modeled through the advection-diffusion-reaction equation. The numerical solution of the differential equation that describes the behavior of the system was found via the SUPG ("Streamline Unwinding Petrov Galerkin") finite element technique. Simulations were done for different Péclet numbers. Then, the adjoint equation of the advection-diffusion-reaction equation was derived for the one-dimensional case and the expression of the coefficient of sensitivity of a generic functional related to a generic parameter was obtained. The sensitivities of the mean and instantaneous pollutant rates were analyzed with relation to the following parameters: drag speed of the flowing current and Péclet number. Results of the sensitivity coefficient obtained with first and second order perturbation methodology satisfactorily matched the same values calculated by the direct method, that is, by means of the direct solution of the advection-diffusion-reaction equation by changing the values of input data parameters.

Keywords: Sensitivity analysis, SUPG method, advection-diffusion-reaction

Introduction

Problems of dispersion and transport of substances immersed in a fluid medium, such as water, can be modeled by means of the advection-diffusion-reaction partial differential equation. The numerical solution can be obtained, for instance, via the "Streamline Upwinding Petrov-Galerkin (SUPG)" technique (Brooks and Huges, 1982; Codina, 1993).

With the purpose of quantifying the amount of dispersed substance, the variations of two functionals are analyzed in this work: the mean value, in the time-space domain of the transported substance, and its instantaneous value, related to variations of the following parameters of system definition: drag speed of the fluid, Péclet mesh number and isotropic diffusivity coefficient.

These variations can be quantified by determining the so-called sensitivity coefficients, which can be obtained directly through a new solution of the advection-diffusion-reaction differential equation, by changing the parameter in the input data. This is the technique known as Construction of Response Surfaces. Another alternative is the broadly used methodology of Sensitivity Analysis via Perturbation Methods widely used in physics and thermal-hydraulics of reactors (Gandini, 1981 and 1987; Oblow, 1976 and 1978; Cacuci et al., 1980; Lima and Alvim, 1986; Lima, 1990; Lima et al., 1993). Recently, the advances in perturbation methods applied to nuclear engineering problems were compiled as a review paper (Lima et al., 1998) and sensitivity analysis for waterhammer problems in hydraulic networks was performed using the differential perturbative method (Baliño et al., 2002). References were not found of applications of perturbations methods to perform sensitivity analysis to pollutant transport.

Sensitivity analysis through perturbation methods presents some advantages because it can be accomplished without previously choosing the parameter to be studied, and the calculations are faster and more efficient, as the equation system that describes the physical behavior of the problems is solved only once. However, the main disadvantage of this methodology is that the answer is linearized around a point of interest. Such disadvantage can be minimized through the application of perturbation methods of higher order (Gandini, 1987). We include this approach in this paper.

Nomenclature

Latin Letters

Greek Symbols

Subscripts

Superscript

Theoretical Model

We study now the transport of a substance in a fluid medium characterized by a field of known scalar speed u (x, t) that satisfies the incompressibility equation.

The transport of a substance and the discharge of the pollutant in a shallow water problem can be described by the advection-diffusion-reaction equation, i. e.,

where f is the depth averaged concentration of the pollutant transported, u is the drag speed of the substance in the fluid (m/s), k is the constant isotropic diffusivity coefficient (m²/s). It is taken as a constant parameter because it is difficult to determine it experimentally in a shallow water problem (Baptista, E. M., 1984). k is also called the dispersion coefficient, as may be seen in (Hunt, B., 1999). g is the absorption coefficient (_s^-1) denominated a first order decay mechanism, ƒ is the source term (_s^-1), h_o is the averaged flow depth (m). W is a domain in R^d, d = 1, 2, or 3 and I is a time interval [0, T].

We simplify the problem taking d = 1 (one-dimensional flow), g = 0 (null absorption) and ƒ = 0 (absence of source), h_o is taken constant. Then Eq. (1) can be rewritten as

where W = [0, L] in R¹.

The initial and boundary conditions of the problems are:

Numeric Solution

Equation (2) with the initial and boundary conditions given by Eq (3) was numerically solved using the "Streamline Upwinding Petrov-Galerkin (SUPG)" technique (Brooks and Huges, 1982; Codina, 1993; Heinrich et al, 1996). The numerical solution can be written in the following matrix form:

where is the matrix of the global mass of the system, is the global capacitive matrix, is the global rigidity matrix, is the vector of the unknown values of the global nodes and is the vector of the speeds of the unknown values of the global nodes. Integrating this equation over time using an implicit finite difference approach we obtain

where and are the values corresponding to times tⁿ⁺¹ and tⁿ respectively, and 0 < q < 1 is the "implicitness" parameter (for Crank-Nicholson method q = 0.5).

Sensitivity Analysis

Eq. (2) can be rewritten as

where

By deriving Eq. (6) with relation to a generic parameter p_i , according to the differential formalism (Oblow, 1978; Cacuci et al, 1980; Lima, 1990), we obtain

where

The adjoint operator of H is given by

Therefore, the adjoint equation of Eq. (8) can be written as

where S⁺ represents the source term of the adjoint equation, related to the functional to be studied. The conditions for the solution of the adjunct equation are chosen in order to cancel the unknown terms of the concomitant bilinear P between f^*and f/ i given by the expression

Examples

To show the use of the perturbative technique in the sensitivity study, we analyzed two results: the mean concentration f_m in the phase space and the instantaneous value at t = 0.5 s and x = 0.975 m. This point has been chosen because it belongs to the numerical boundary layer and presents strong gradients of the variable f.

Adjoint function

i) mean velocity f_m in the phase space

Sensitivity Study of the Mean Concentration f_m in the Phase Space

The mean value concentration f_m in the phase space is defined by

where T represents the final time and L represents the length of the channel as indicated in the definition of W x I.

From (Cacucci et al, 1980) the following can be concluded:

Then, Eq. (11) takes the form:

with the following final and boundary conditions chosen to nullify the terms of the concomitant bilinear:

Note that taking the adjoint equation is a "reversal operation" that changes the sign of the time derivative.

ii) Instantaneous velocity at t = 0.5 s and x = 0.95 m

The instantaneous velocity, , at t = 0.5 s and x = 0.95 m can be defined by

where d (x) is the Dirac delta.

From Eq. (17) we get

Then, using Eq.(18), Eq. (11) takes the form

with the following final and boundary conditions, chosen to nullify the terms of the concomitant bilinear (Gandini,1987):

with e arbitrarily small. The numeric treatment of Eq (20) is made using integration (i. e., as a weak solution), as can be seen for instance in (Jacovkis and Rosales, 1993).

Sensitivity Coefficient

The general functional variation dR respect to the independent parameter dp_i, as we can see for instance in (Gandini, A., 1987), is given by

where M/_j is the derivative with respect to p_j of the stationary linear operator M given by

Considering d_pi = d_pj, Eq.(21) can be simplified to

and the sensitivity coefficient is given by

The first order approximation for the sensitivity coefficient can be written as

From Eq.(24), we can obtain the sensitivity coefficients with respect to the drag speed of the fluid (u), diffusivity coefficient (k) and Péclet number (Pe), i. e.,

where Dt is the time step, and h the element length.

From Eq.(24) we can obtain the second order sensitivity coefficients, which are

Results and Discussion

To exemplify the study developed herein, consideration was given to the case of the one-dimensional transport of a pollutant substance in a fluid medium (water), moving in a channel of length L, at speed u.

The initial and boundary conditions of the problem are

Figure 1 shows the results of Eq. (2) that represent the evolution of the pollutant concentration along the channel for several times. These results are very important for obtaining the sensitivity coefficients in equations (21), (22) and (23). The flux analyzed corresponds to the high Péclet number (12.5) that is significant in the upwinding functions adopted. A migration of the pollutant to the exterior of the channel is due to the principle of mass conservation written in Eq. (2). The initial and boundary conditions of the pollutant concentration are taken as in (Codina, 1993) with the objective of comparing results for both methods. It represents the evolution of the pollutant concentration in the fluid medium when this substance is introduced against the stream and is maintained constant along the time.

Figure 2 shows the graph of the adjoint function (as can be seen in Eqs.(15) and (16)) which corresponds to the mean concentration functional of the pollutant in the fluid medium. represents the physical contribution (of the increment of the pollutant situated at point r_o) to the functional. In this case the contribution affects the mean concentration of the pollutant given by Eq.(13) in the phase space.

Figure 2

This justifies the fact that the adjoint function is called "importance function" by some authors (Gandini, 1987; Lima, 1990).

Figure 3 shows the graph of the adjoint function which corresponds to the instantaneous functional at t = 0.5s and at x = 0.975m due to the instantaneous pollutant concentration in the fluid medium given by Eq. (19) and Eq. (20). For a time interval the contribution of a pollutant increment will be larger at points closer to the exit of the channel and to the final time. On the other hand, as expected, for fixed points inside the channel, the contribution is more significant when the instant considered is closer to the final time.

Figs. 4 and 5 show the graphs of the mean and instantaneous concentration of the pollutant f_m and f_r0 which corresponds to Eqs. (13) and (17) as a function of the diffusivity coefficient k, obtained via the disturbance formalism (perturbation method) and by the direct computations (via SUPG method). The difference between the results obtained by the perturbation and direct methods is smaller because the values of Pe gradually decline in order to guarantee the stability of the solution of Eq. (2).

Figs. 6 and 7 show the graphs of the mean and instantaneous concentrations of the pollutant f_m and f_r0 as a function of the drag speed of the fluid u, which corresponds to Eqs. (13) and (17), obtained via the disturbance formalism and via the direct computation. The variation of f_m as a function of parameter u is practically symmetrical with regard to the point of reference. The reduction in the value of f_m as u increases is an expected result due to the increase in Péclet number.

Figs. 8 and 9 show the graphs of the mean and instantaneous concentrations of the contaminant f_m and f_r0 as a function of the Péclet number, which corresponds to Eqs. (13) and (17), obtained via the disturbance formalism and via the direct computation . The variation of f_m as a function of parameter Pe is also practically symmetrical with regard to the point of reference. The decrease in the value of f_m as Pe increases is an expected result due to the increase in Péclet number.

Table 2 shows the results of the first order sensitivity coefficients of the functional analyzed with regard to the three selected parameters. The values considered in the reference solution are also shown in this table. For comparison purposes, the dimensionless sensitivity coefficients were also included. Positive values indicate that the functional increases with the parameter, and negative values indicate a decrease in the functional as the parameter value increases. It should be observed that the Péclet mesh number is the parameter for which both functionals studied show larger sensitivity, and the diffusivity coefficient is the parameter that has the least influence on the functionals. On the other hand, it is also seen that the instantaneous concentration functional shows less sensitivity to the parameters than the mean concentration functional.

Table 3 shows the errors found in the first order sensitivity coefficients of the functionals studied with relation to the chosen parameters. It can be seen that there are major errors of the order of 3.8%, what suggests the use of higher orders for the sensitivity study via perturbative methods.

Table 4 shows the results of the second order sensitivity coefficients of the functional analyzed with relation to the three selected parameters. For comparison purposes, the dimensionless sensitivity coefficients were also included. Positive values indicate that the functional increases with the parameter, and negative values indicate a decrease in the functional as the parameter value increases The values considered in the reference solution are also shown in this table. We consider only a mean value functional because the sensitivity coefficient, which corresponds to the instantaneous functional value, does not need to be corrected with a theory of second order.

In a first approximation, the perturbation method is efficient to compute sensitivity coefficients of the functional f_m and f₀, in relation to the three parameters: drag velocity u, Péclet number Pe and the diffusivity k. The second order perturbation theory corrects the sensitivity coefficients related to the mean functional. It can be observed that this correction is small because the variation of the unknown variable respects to the parameters is small too.

Conclusions

The methodology proposed in this paper for the sensitivity analysis, via the differential formalism, allowed us to compute the sensitivity coefficients of the mean and instantaneous concentrations of the pollutant in a fluid medium (water) without having to build the response surfaces imposed by the direct methodology, that is, the solution of the diffusion-convection-reaction equation repeated several times, that describes the physical behavior of the system. Computational saving is an evident consequence of this approach.

The values that were found showed that the functionals studied present higher sensitivity to variations of the Péclet mesh number. On the other hand, the sensitivity coefficients found indicate correctly the direction of the variations of the answer functionals. The results obtained with second order perturbations show a small difference respect to first order errors. It can be explained due to the small nonlinearity of the Burgers equation. The results obtained with the second order computations do not justify the use of this more sophisticated analysis in this problem.

Finally, it is important to emphasize that the problem under discussion in this paper is an assay for the application of sensitivity studies via perturbation methods in more complex situations.

Acknowledgements

The authors would like to acknowledge the Conselho Nacional de Desenvolvimento Científico e Tecnológico-CNPq of Brazil and the University of Buenos Aires (grant I050) for financial support.

Paper accepted December, 2002. Technical Editor: Atila P.Silva Freire

• Baliño, J. L., Larreteguy, A. E., Lorenzo, A. C., Padilla, A. G. and Lima, F. R. A., 2002, "The Differential Perturbative Method Applied to Sensitivity Analysis for Waterhammer Problems in Hydraulic Networks", Applied Mathematical Modelling, to be published.
• Baptista, E. M., 1984, "Eulerian-Lagrangian Analysis of Pollutant Transport in Shallow Water", Master of science thesis in Civil Engineering, Massachusetts Institute of Technology, Cambridge, Mass.
• Brooks, A.N. and Huges, T.R.J., 1982, "Streamline Upwinding/Petrov-Galerkin Formulations for Convection Dominated Flows with Particular Emphasis on the Incompressible Navier-Stokes Equations", Comput. Meths. Appl. Mech. Engrg., vol. 32, pp. 199-259.
• Cacuci, D.G., Weber, C.F., Oblow, E.M. and Marable, J.H., 1980, "Sensitivity Theory for General Systems of Nonlinear Equations", Nuclear Science and Engineering, vol. 75, pp. 88-110.
• Codina, R., 1993, "A Finite Element Formulation for the Numerical Solution of the Convection-Diffusion Equation", PhD Thesis International Center for Numerical Methods in Engineering, Barcelona.
• Gandini, A., 1981, "Generalized Perturbation Theory for Nonlinear Systems for the Importance Conservation Principle", Nucl. Sci. Eng., vol. 77, pp. 316.
• Gandini, A., 1987, "Generalized Perturbation Theory (GPT) Methods. A Heuristic Approach", Advances in Nuclear Science and Technology, vol. 19, pp. 205-380.
• Heinrich, J.C., Idelsohn, S.R. and Vionnet, C.A , 1996, "Boundary Conditions for Finite Element Simulations of Convective Flows with Artificial Boundaries", International Journal For Numerical Methods In Engineering, vol. 39, pp. 1053-1071.
• Hunt, B., 1999, "Dispersion Model for Mountain Streams", Journal of Hydraulic Engineering, Vol. 125, pp. 99-105.
• Jacovkis, P.M. and Rosales, R.R., 1993, "Análisis Numérico de Ondas de Detonación Retardada", Revista Internacional de Métodos Numéricos para Cálculo y Diseño en Ingeniería, Vol. 9, pp. 335-356.
• Lima, F.R.A., 1990, "Aplicações de Métodos Perturbativos ao Modelo Multicanal COBRA IV_1 para Cálculos de Sensibilidade em Núcleos de Reatores Nucleares", Tese de doutorado, COPPE/UFRJ, Brasil.
• Lima, F.R.A. and Alvim, A.C.M., 1986, "Tempera-V2: Um Programa para Análise de Sensibilidade num Canal Refrigerante de Reatores Nucleares", PEN-139, Programa de Engenharia Nuclear COPPE/UFRJ.
• Lima, F. R. A., Lira, C.A.B.O. and Gandini, A., 1993, "Sensitivity Analysis of Thermohydraulic Systems via Heuristic Generalized Perturbation Theory (HGPT) Methods", Ann. Nucl. Energy, vol. 20, No 10, pp. 679-690.
• Lima, F.R.A., Gandini, A., Blanco, A., Lira, C.A.B.O., Maciel, E. S. G., Alvim, A.C.M, Silva, F.C., Melo, P.F.F., França, W.F.L, Baliño, J.L., Larreteguy, A.E., Lorenzo, A., 1998, "Recent Advances in Perturbative Methods Applied to Nuclear Engineering Problems", Progress in Nuclear Energy, Vol. 33, N 1/2, pp. 23-97.
• Oblow, E. M.,1976, "Sensitivity Theory from a Differential Viewpoint", Nucl. Sci. Eng., vol. 59, pp. 187-189.
• Oblow, E. M.,1978, "Sensitivity Theory for Reactor Thermal-Hydraulics Problems, Nucl. Sci. Eng., vol. 68, pp. 322-337.

Publication Dates

• Publication in this collection
  18 Mar 2004
• Date of issue
  Mar 2003

History

• Received
  Dec 2002