Скачать 46.67 Kb.

Прикладне програмне забезпечення UDC 004.75Preservation of the qualitative properties in the vertical advectiondiffusion DEM submodel of the air pollution process Maryam A.T. Elshebli Department of Applied Analysis, Eötvös Loránd University H1117 Budapest, Pázmány P. s. 1/c, Hungary The analysis of the transport of air pollutants is one of the most widely investigated problems and it has great practical importance. The operator splitting method is a widely used tool for the solution of such problems. Because the numerical model has to preserve the main qualitative properties of the original model, we need such methods which preserve them as well. In this paper we investigate this question by use of the operator splitting approach. Keywords: air pollution modelling, operator splitting, numerical solution. Introduction The transport of air pollutants is one of the most widely investigated phenomena and it is of great practical importance (Prussov and Dorosenko 2003, Zlatev 1995). The numerical handling of this problem needs many efforts due to the extremely big size of the unknown in the discretized models. The physical model describes of the simultaneous effect of several different subprocesses. The mathematical equivalents of these processes are operators describing the subprocesses. They are as a rule simpler than the whole spatial differential operator. Their separation serves as the basic idea of the operator splitting method. The preservation of the main qualitative properties of the continuous model is a natural requirement for the global numerical models. In order to get such a discretization method we guarantee the preservation in each submodel. Clearly, this serves as a sufficient condition for the global preservation. The application of operator splitting for the air pollution model results in an advectiondiffusion submodel, as a separate task. For the solution of such a problem we need numerical methods. There are several robust numerical methods and some of them are very powerful (c.f. Prussov and Dorosenko 2006). However, usually the convergence of these methods are investigated only. In this paper we investigate the qualitative properties, namely, the preservation of the maximum principle in the discrete model. This paper is organised as follows. In Section 2 we describe the mathematical model of the air pollution process. In the Section 3 we formulate the operator splitting methods. We give the algorithm of the most widely used operator splitting methods like the sequetial, MarchukStrang splittings and the weighted sequential splitting. We also mention some further splittings like the additive splitting and the iterated splitting. In Section 4 we analyse the problem of how to choose the suboperators in the operator splitting process. We show that this choice can be performed differently in the air pollution model given in Section 2. In Section 5 we formulate the requirement of the qualitative property preservation for the global discretization method and we show how the operator splitting idea can be used to this aim. We analyse in detail one submodel in the DEM (Danish Eulerian Model), namely, the vertical advectiondiffusion decomposition part. We examine the preservation property of the maximumminimum principle (the socalled discrete maximumminimum principle) and we give the condition of the choice of the discretization parameter for the bilinear finite element space and method time discretization method. Finally we conclude how nonnegativity preservation is related to this property.
In this section we present the mathematical model of the air pollution process. Let denote the concentration of the jth air pollutant, and c the vector function of these functions. Then the time evolution of the vector c can be described mathematically by the system of partial differential equations (Zlatev, 1995) (1) where is a vector function describing the wind velocity, is the diffusion coefficient function, is the function of emission,describes the deposition and R defines the chemical reactions of the pollutants. The initial function c_{0}(x) is given. Using these notations, the terms in equation (10) have the following physical meaning. The first term on the righthand side describes the transportation due to the velocity field, which is called advection. The second term expresses the turbulent diffusion, the third term the emission, the fourth term describes the deposition and the last term defines the chemistry of the pollutants. We can see that (1) is a system of partial differential equations which is nonlinear due to the chemistry term. The analytical solution of such a system is impossible. In order to simplify the task, we use the idea of operator splitting. The number of chemical species involved in an air pollution model sometimes reaches 200, or even more, which results in a huge system of partial differential equations. The analytical solution of such a problem is obviously impossible to find. Hence we have to treat it numerically. We note that in case of semidiscretization usually the number of spatial grid points equals many millions. This means that the system of ordinary differential equations obtained after spatial discretization is extremely big, hence the use of any numerical method developed for systems of ODE’s is rather complicated. Moreover, the model equations contain terms that have different physical meanings and so different mathematical properties (e.g. linear, nonlinear, stiff and nonstiff). Therefore, it is impossible to find such a universal numerical method which would perform well when applied directly to the original system. The application of operator splitting allows us to treat the different physical terms separately. Operator splitting method (OSM) is a kind of problem decomposition: we divide the spatial differential operator of the global system into a few simpler operators and solve the corresponding problems one after the other, by connecting them through their initial conditions. The simpler systems which are obtained in this manner, and are sometimes called subsystems, might have some special properties that can be exploited in the numerical solution. The subsystems are usually easier to treat numerically than the whole system. Splitting can be performed in several ways. We expect the method to be accurate as well as efficient enough. The latter property depends on the number of computations and the possibility of performing the computations in parallel. Taking into account the latter requirement, we made attempts to construct a new splitting scheme which does not require a lot of computational work, and is parallelizable on the operator level. 2. Description of the OSM In the sequel the frequently used splitting methods are described and compared. We describe the methods only for two operators (i.e., n = 2), however, the generalisation for n operators is straightforward. (For more details, see Hundsdorfer and Verwer, 2003, Zlatev, 1995, and Dimov at al., 2006.) Hence, we consider the ACP (2) where now A and B are given operators. In fact, the OSM yields a timediscretization method where instead of solving the continuous (in time) problem (2) directly, we seek the split (discretized) solution on the grid points of the mesh _{}_{ } {t_{j }= j, j = 0, 1,…M}, where M = T. (3) Here τ > 0 denotes the splitting time step. In the following we summarise some widely used splitting methods. 2.1 Sequential splitting. The scheme of this method is the following. As a first step, we solve the system with operator A using the initial condition of the original problem, and then, applying the obtained solution at time τ as an initial condition, we solve the system with operator B. The solution obtained in this way is considered as the splitting solution in τ. This procedure is performed cyclically in the following way: (4) (5) for k = 1, 2,…, m, where . The split solution at time t = kτ is defined as for k = 1, 2,…, M. (6) It can be shown that the sequential splitting is a firstorder method (see e.g. Havasi et al., 2001). 2.2 StrangMarchuk splitting. Using this method, at each time step we begin and end the computation with operator A (we apply it over a distance τ/2 twice) and put B to the middle (we apply it over a distance τ once) as follows: (7) (8) (9) for k = 1, 2,…, M, where . The splitting solution at t = kτ is defined as for k = 1, 2,…, M. (10) This method has secondorder accuracy. 2.3 Weighted sequential splitting. This method can be obtained by symmetrizing the sequential splitting in the following way: in each time step we apply sequential splitting both in the order A → B as follows: (11) (12) and B → A as follows: (13) (14) The splitting solution at t = kτ is defined as , for k = 1, 2,…, M, (15) where , and is some fixed weight parameter. The weighted sequential splitting is a secondorder method for , otherwise it has first order. An important property of the weighted sequential splitting is that it can be parallelized in a natural way (on the operator level). The theoretical investigation of the method can be found in Csomós et al. (2003), and an application of the method in a onecolumn transportchemistry model is described in Botchev et al. 2003. 2.4 Other kinds of splittings. There are also other OSM’s, which has been recently developed. We want to mentioned two methods of them, namely, the additive splitting and the iterated splitting. The additive splitting is similar to the sequential splitting with the following only difference: in both subproblems (4) and (5) we use the same initial value i.e., the split solution at the previous time level. If we denote by and the corresponding solutions, then the split solution at the new time level is defined as for k = 1, 2,…, M. (16) The method has first order accuracy. The main advantage of this OSM is its easy parallelization. (For more details we refer to Gnandt, 2005.) The iterated splitting suggests the following algorithm: on the interval [t_{k1},t_{k}] we solve the following subproblems consecutively, for i = 1,3,5, … 2m+1. (17) (18) where the initial function u_{0}(t) is any fixed function for each iteration, and . The split solution at t = kτ is defined as for k = 1, 2,…, M. (19) The order of the accuracy of this OSM depends on the number of the inner iteration, and equals 2m+1. The main advantages of the constructed OSM are the following. Each of the subproblems result in a consistent approximation of the exact solution. Moreover, using the builtin iteration, we can achieve (at least, theoretically) arbitrary high order for the local splitting error. 3. Choice of the suboperators For some given complex mathematical model the suboperators can be chosen in different ways. In the particular case of the air pollution modelling, we sketch the following two basic possibilities. 3.1 Physical decomposition. It is natural to define the different suboperators on the base of the separate physical processes, namely, we can define the following operators
3.2 DEM decomposition. This kind of decomposition is used in the Danish Eulerian Model (DEM) and is called DEM decomposition.
The main advantage of the DEM decomposition is its high flexibility for 2D problems because only the last operator contains the vertical part. We note that in the DEM splitting the operators A_{1} and A_{2 }are twodimensional operators and which is crucial in our later analysis the operator A_{5} is a onedimensional operator. Hence the numerical investigation of the latter is simple. 4. Qualitative analysis in the vertical advectiondiffusion DEM submodel Using the OSM, the split subproblem for the vertical suboperators means the following onedimensional partial differential equation of parabolic type (20) where denotes the partial derivative w.r.t. the vertical variable z. Here, for simplicity, we assume that the advection and diffusion coefficients are constant. For this case the maximum principle can be proved, which states the following: on some bounded domain the unknown function c(x,t) attains its maximum and minimum on the parabolic boundary, i.e., either for t = 0 or for where denotes the boundary of the bounded space domain . It is worth preserving this property for the numerical solution, i.e., to guarantee the following property for the fully discretized numerical model: whenever we introduce the mesh _{}_{z,}_{}_{t} on the solution domain _{T }= (0,T), the discretization of the problem (20) should take its maximal and minimal value at the corresponding discrete parabolic boundary points of the mesh _{}_{z,}_{}_{t}. If denotes the approximation at the meshpoint (z_{i},t_{n}) (where z_{i} is the ith space meshpoint and t_{n} is the nth time level, respectively) then the maximal and minimal value of the numerical solution are achieved either for n = 0, or where the indices correspond to the space boundary points. If we use either the finite difference approximation or the linear finite element method to the semidiscretization for the space variable and the usual method for the time discretization (which is, in fact the generalized trapezoidal formula), then we can show the following sufficient condition which guarantees the discrete maximumminimum principle. Theorem. Assume that the space discretization parameter is sufficiently small. Then for the bilinear finite element approximation combined with the method the discrete maximum principle holds if the conditions (21) and (22) are fulfilled. The proof is complicated and it can be found in Elshebli, 2006. We note that it usually happens in the qualitative analysis of finite element approximations that there are both upper and lower bounds for the timestep, which means that t can be chosen neither too small nor too large. Such a kind of error bounds can be found in some works, e.g. Faragó et al. 2004 and 2005. The conditions (21) and (22) imply that this criterion can be applied only in the case where the lower bound is less than the upper bound. This requirement implies that for the choice of the requirement is . (23) This shows that for negative u the condition (24) should be satisfied. It is worth mentioning that many other important qualitative properties can be shown from the discrete maximum principle. One of the most crucial requirements is the nonnegativity preservation property, which means the following: if the output functions of the continuous problem are nonnegative, then the solution of the problem is also nonnegative. The notion of the discrete nonnegativity preservation is clear: if the vectors which correspond to the initial function and the boundary condition are nonnegative then the numerical solution should also be nonnegative. It can be shown that in both the continuous and discrete cases the discrete nonnegativity preservation property and the discrete maximum principle are equivalent. Hence, the conditions (21) and (22) together guarantee the nonnegativity property in the discretized model, too.
3. Elshebli, M., 2006, Discrete maximum principle for bilinear finite elements on uniform square meshes (submitted) 4. Faragó.I., Horváth,,R. and Korotov, S., 2004, Discrete maximum principle for Galerkin finite element solutions to parabolic problems on rectangular meshes, Numerical Mathematics and Advanced Applications, Springer Verlag, Berlin, 298307. 5. Faragó.I., Horváth,,R. and Korotov, S., 2005, Discrete maximum principle for linear parabolic problems solved on hybrid meshes, Appl. Numer. Math. 43, 249264. 6. Gnandt, B.2005, A new operator splitting method and its numerical investigation, Advances in Air Pollution Modeling for Environmental Security, 54, Kluwer, Amsterdam, 229241. 7. Havasi, Á., Bartholy, J., Faragó, I., 2001, Splitting method and its application in air pollution modeling. Időjárás 105, 39–58. 8. Hundsdorfer, W., Verwer, J. G., 2003, Numerical solution of timedependent advectiondiffusionreaction equations. Springer, Berlin 9. Prusov, V., Dorosenko, A., 2003, Modelling and forecasting atmospheric pollution over region, Annales Univ. Sci. ELTE 46, 2745. 10. Prusov, V., Dorosenko, A., 2006, On the numerical solution of the onedimensional convectiondiffusion equation, Int. J. of Environmental. Pollution (to appear) 11. Zlatev, Z., 1995, Computer treatment of large air pollution models. Kluwer, Amsterdam. © K. Georgiev, E. Donev, 2006 ISSN 17274907. Проблеми програмування. 2006. №23. Спеціальний випуск 