On local finite element refinements in multiscale air quality modeling

EnvironmentalSo]iware9 (1994) 61

66 © 1994 Elsevier Science Limited Printed in Great Britain. All rights reserved 0266-9838/94/$7.00 ELSEVIER

On local finite element refinements in multiscale air quality modeling Mehmet T. Odman MCNC, North Carolina Supercomputing Center, Research Triangle Park, North Carolina 27709-2889 USA

& Armistead G. Russell Department of Mechanical Engineering, Carnegie Mellon University, Pittsburgh, Pennsylvania 15213, USA

(Received 5 June 1992; revised version received 30 April 1993: accepted 28 July 1993)

A b s t r a c t - Variable resolution is a highly desirable property in air quality models, especially in regional applications. Resolution can be increased in dense source regions by using finite element refinements. Here, the important principles that must be obeyed at refinement boundaries are discussed. Mass conservation is achieved by making the element basis functions continuous. Constraint relations that assure continuity for various refinement ratios are described. A second issue is to keep the refinement boundaries free of noise. Since coarse and fine elements act like different media, aliasing errors usually lead to noise waves. A non-linear filter is used to remove some of this noise. Tests are conducted with different refinement ratios to see the effect of increased resolution on accuracy. In general, refinements increase accuracy by reducing diffusion errors. The peak concentrations are overpredicted during the transition from the fine to the coarse grid. These overpredictions are smaller when the refinements are gradual.

Key Words - - Finite element refinements, multiscale modeling, variable resolution.

Software Availability Name o(the software: Urban and Regional Multiscale (URM) model Brief dega-iption; URM is a photochemical air quality model for use in urban and regional modeling. It allows variable resolution both in the horizontal and vertical domains. The LCC mechanism is used to describe the chemistry. Dgvelotx~Contact Addresses; Armistead G. Russell, Department of Mechanical Engineering, Carnegie-Mellon University, Pittsburgh, Pennsylvania, 15213 U.S.A. Phone: (412) 268-3614 Fax: (412) 268-3348 E-Mail: russeU @pollution.me.cmu.edu Melimet T. Odman, North Carolina Supercomputing Center, MCNC, Research Triangle Park, North Carolina, 277092889, U.S.A. Fax: (919) 248-9245 Phone: (919) 248-9235 E-Mail: odman @ncse.org Year first available: 1992 32 Mbytes ROM, 512 Mbytes I/O Disk I:Iallaalr,,glaa~gk Ikaglam.L~amlagg FORTRAN and C 512 Kbytes 61

M. 7". Odman, A. G. Russell

62 1

INTRODUCTION

Regional applications of air quality models include acid deposition and ozone simulations. It is well known that urban emissions play an important role in both of these pollution problems. To study the potential impact of emissions, modeling domains must be sized to include both the source and the receptor areas. Consequently, typical domains cover large geographic regions (~1000 kin). On the other hand, many processes involved in the formation of air pollution occur on fine scales. To resolve the high concentration gradients around source areas, scales as small as 2 km are often desired. Since the horizontal resolution in most current models is fixed, the entire domain may need to be modeled with the finest scale. This approach is computationally inefficient. The gradients are not large over rural areas that usually constitute a Large portion of modeling domains. Therefore, f'me resolution is not really needed in rural areas. A better approach is to change the structure of current models by allowing variable resolution. Recently, Odman and Russell[ l] proposed a multiscale modeling approach. In this approach, various scales are incorporated into a single model by using finite elements. Coarse elements are used in rural areas for efficiency. Resolution is increased in source areas by refining the size of the elements. This is different from one-way nesting [23], because coarse and fine elements fully interact with each other. Mass conservation usually becomes a problem in one-way nesting where a regional model with a coarse grid provides boundary conditions to an urban model with a fine grid. In multiscale modeling, continuity of finite element basis functions assures mass conservation at the coarse-fine element boundaries. Another important issue is the transition of resolvable waves from one scale to another, ideally, without any noise at the interface. This is usually more difficult to achieve than mass conservation. In this paper, the effect of refinements on the smoothness of resolvable waves as they travel through the boundary is investigated. Several refinement ratios are considered. In each case, the so called aliasing effect is observed.

2

LOCAL FINITE ELEMENT REFINEMENTS

While developing the multiscale air quality model, Odman and Russell [I] proposed two refinement techniques. These were the constraint element and the transition element techniques. Later, the transition element technique was abandoned due to significant amplification of mass conservation errors in the presence of non-linear c h e m i s t r y [ 4 ] . In this section, the theory behind the constraint element will be described. This element will be generalized for use in various refinement ratios. Also, the weaknesses of the transition element in air quality modeling will be discussed in more detail. There are essentially two ways of increasing the local resolution in finite element methods. The first one, locally increasing the order of element basis functions may result in large matrices and will not be discussed here. The

CG

3

FG

2

Fig. 1 Finite element refinement with a ratio of 1:2: each coarse-grid (CG) element faces two fine-grid (FG) elements.

second one, refining the element size while keeping the order of basis functions the same, is more practical in air quality modeling. Assume that a domain has been subdivided into quadrilateral finite elements with bilinear basis functions. If the element basis functions are continuous everywhere, this would result in mass conservation. When elements are locally refined, continuity can be maintained by introducing constraint relations[5"6]. Assuming that node 3 is bisecting the side with corner nodes 1 and 2 (Figure 1), the following constraint is obtained: c 3 = l ( c I + Cz).

(1)

Thus, to achieve continuity of basis functions, the concentration at node 3 must be equal to the average of the concentrations at parent nodes 1 and 2. Since the concentration at node 3 is not a true degree of freedom (unknown), it is eliminated from the analysis in the finite element code. However, this is not done until the assembly stage to maintain uniformity and vectorization of the code. Element stiffness matrices are still computed as if node 3 was a degree of freedom, however, in the global stiffness matrix, the residuals are distributed equally to the parent nodes. The constraint in Equation (1) is valid only for a refinement ratio of 1:2. If other refinement ratios are used, there will be more than one node that needs to be constrained. There would be N - 1 such nodes for a 1 :N refinement ratio (i.e., each coarse element faces N fine elements). A more general constraint relation can be written as:

¢i+2 = Wi+2,1Cl + Wi+2,2C2, i = 1, 2 ..... N - 1,

(2)

where wi+2,1 and wi+2,2 are the weighting factors of the parent node concentrations. These weights may be obtained easily by linear interpolations using the nodal coordinates. Note that the interpolations should be linear to assure continuity of the bilinear basis functions. For

On/oea/finite elenwnt refinements in multiseale air quality modeling example, in a 1:4 refinement where the constrained nodes are numbered as 3, 4 and 5, with 3 closer to parent node 1, the weight w5, 2 would be equal to 0.75. The transition element technique proposed earlier used quadratic basis functions along the refinement boundary. For the situation shown in Figure 1, the coarse element would be a transition element that connects a bilinear coarse element to two bilinear fine elements. Node 3 is connected to the transition element by quadratic basis functions along the boundary 1-3-2. By making node 3 a true degree of freedom it was hoped that smoother transitions could be obtained. However, the continuity requirement is obviously violated with these transition elements. There are several cases in finite element applications where the continuity requirement is violated by using nonconforming elements[7]. It is still possible to obtain convergent solutions as long as these violations are local. The transition element yielded convergent solutions and the aliasing errors were slightly smaller than the constraint element. However, the mass conservation errors were significantly amplified in the air quality model due to high nonlinearity of the chemistry. It became apparent that any non-conservative technique must be avoided in air quality modeling. Conserving mass is not the only requirement for refinement techniques. It is also important that resolvable waves pass through the refinement boundaries with no significant noise[8]. Since coarse and fine grids act like different media to traveling waves, noise will be generated at such boundaries. This phenomenon, known as aliasing, can be filtered by adding artificial diffusion[9]. A review of filtering techniques at refined grid boundaries for finite difference methods can be found in Koch and McQueen[10]. Recently, Odman and Russell[ 11] developed a nonlinear filter that eliminates the spurious oscillations observed near sharp gradients in linear finite element solutions. The technique fu'st determines the location of potential noise waves and then applies an artificial diffusion operator to the immediate vicinity of this location. This way, the global high-order accuracy of the solution is maintained. Here, the same filter is used to eliminate the noise generated by aliasing errors at refinement boundaries.

3

NUMERICAL TESTS AND RESULTS

The rotating puff problem is used intensively to test the robustness of advective transport algorithms. It will be used here to evaluate the effectiveness of finite element refinements. Firsk a square domain is divided into 32x32 bilinear finite elements. With the origin at the center of the domain, the lower left and the upper fight comers have the coordinates (-16, -16) and (+16, +16) respectively. Then, the 13x26 elements in a rectangular subdomain are replaced by finer elements. The lower left and the upper right coordinates of this subdomaln are (-13, -13) and (-0, +13) respectively (Figure 2). For a 1:2 refinement ratio, the fine elements are 1/2 the size of the coarse elements. Thus, two fine elements are facing one coarse element.

63

CG FG

(0, 13) Y

~

In.

5

j

x

U

(-13, -13)

(a)

,

,

,

IIIImlmmmlllmmlu mlmuulnnmmmmmlnnn Ilnmmlllnlnmmmnl ::::::::::::::::::::::::::::: ::::::::::::::::::::::::::::: ::::::::::::::::::::::::::::: ::::::::::::::::::::::::::::: ::::::::::::::::::::::::::::: ::::::::::::::::::::::::::::: ::::::::::::::::::::::::::::: ::::::::::::::::::::::::::::: ::::::::::::::::::::::::::::: ::::::::::::::::::::::::::::: ::::::::::::::::::::::::::::: ::::::::::::::::::::::::::::: ::::::::::::::::::::::::::::: ::::::::::::::::::::::::::::: ::::::::::::::::::::::::::::: ::::::::::::::::::::::::::::: ::::::::::::::::::::::::::::: ::::::::::::::::::::::::::::: ::::::::::::::::::::::::::::: ::::::::::::::::::::::::::::: ::::::::::::::::::::::::::::: ::::::::::::::::::::::::::::: ::::::::::::::::::::::::::::: :::::::::::::::::::::::::::::

::::::::::::::::::::::::::::: ', i : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : imimmimiinninimnn

Immmmmmmnmmmmmmu I l i l l l l l l l i l l l i l

-

lmnmnmmmm immnmmmmm imnmmmmmm immin|nnn lmmmm|mn| Imnnnmmmm immmnimmn Inmmnnmmm ilmimlmii lmiimlmio umnnnmmmm llliillll lUliilllU ilaamanui lllllllii mmimilimm Illllllll Illllllll IIIllllll Illllllll Illllllll IIIIIIlll IIIllllll Illllllll IIIIIIIll Illllllll Illllllll

]

, : : - ~- ~

_

- -

linnmnnnnnnnnnn limmmmmim limmiimiim

Immmmmnmm IlllllllU

Fig. 2 Definition of the rotating puff test with element refinements: a) 13x26 elements of the coarse-grid (CG) are replaced by finer elements in the subdomain marked by FG. The initial location and shape of the puff is shown. The contour lines represent concentrations of: (1) 1.00, (2) 0.85, (3) 0.50, (4) 0.15, (5) 0.00. b) Element refinements is illustrated further with a 1:2 refinement ratio. There are 26x52 fine elements. For refinement ratios of 1:4 and 1:8, the subdomain is divided into 52x104 and 104x208 bilinear fine elements respectively. The puff is initialized as a cosine-hill:

c(x,y)=I2(l+cos---~),

L

0,

R<4 R>4

(3)

M. T. Odman, A. G. Russell

64

where, R is the distance from the peak of the puff located at the point (-8, 0). The velocity field is a rigid-body rotation field: U ---- --(O'y

(4)

Y = +(/J0f

The angular velocity, t~, is adjusted such that one full 360 ° rotation of the puff is performed in 240 time steps. This corresponds to a Courant number of nil5, when the peak of the puff is in the coarse grid. With this

experimental design, the puff remains in the refined subdomain for one-half of its revolution. Waves of wavelength less than or equal to 4h are assumed to be noise, where h is the typical coarse grid length (unity). The nonlinear filter operates on these waves to dissipate them with the same artificial diffusivity everywhere (0.16), regardless of the refinement ratio. The solutions obtained after two complete rotations are compared to the exact solution in Figure 3. When the same test was conducted with no refinement (not shown), the peak height of the puff dropped to 75.6% of its

I : I 720 Peok Ntn. -

Peok Mtn.

7.560E-01 -2.000E-05

1:2 720 8. OBSE-O 1 -" -3.200E-05

i~0"00 ~oo"o s0. °

p.0 0 .0 Qo

.~0"°

-

1

(aI

Peok NLn.

Q

" •

(bt

1:4 720 - 8.112C-01 - -3.007E-03

Peok Ntn.

.

Q

1:8 720 - 8.303E-01 - -4.581C-03

"

tcl Id) Fig. 3 The shape of the puff after two complete rotations as predicted by the a) Exact solution, and numerical solutions using element refinement ratios orb) 1:2, c) 1:4, and d) 1:8.

On local finite eh,ment refinements in multiscale air quali O, modeling original value. The puff is followed more accurately in the ref'med grids and a higher percentage of the original peak height is retained at the end of the tests. The peak height is 80.9%, 81.8% and 83.9% of the original when 1:2, 1:4 and 1:8 refinement ratios are used, respectively. The final shape of the puff is in good agreement with the exact solution in all cases (Figure 3). Though aliasing errors increased with decreasing refinement ratio, the filter was able to eliminate noise waves that might lead to negative concentrations. This is illustrated by the smooth background in all the cases of Figure 3. The time history of the peak heights is shown in Figure 4. The maximum concentrations in the field are sampled at each 10° degrees of rotation for two complete revolutions. The exact location of the peak may not always coincide with a nodal point, therefore, maximum concentrations displayed may be lower than the peak concentration. This leads to the oscillations in the curves of Figure 4. It is important to note that the amplitude of the oscillations decrease in the refined grid, because there is a better chance that the peak will coincide with a grid point in the refined grid. Here, the deterministic nature of the curves will be used in the analysis. When no refinement is used (1:1) there is a steady drop in the peak height due to the diffusion errors. When element refinements are used, the rates of reduction of the peak heights decrease noticeably between 0 ° and 90 °, between 270 ° and 450 °, and between 630 ° and 720 ° of rotation. During these periods, the peak of the puff is within the fine grid where the resolution is increased. Diffusion errors are much smaller in the fine grid, therefore the peak heights are better retained. An overshoot is observed when the peak is entering the coarse grid mesh, and there is a slight drop when the peak is entering the fine grid. Obviously, the fine grid acts as a more compliant media to traveling waves. Though the filter is very effective in eliminating aliasing errors at the base of sharp gradients, it offers no mechanism that will stop the peak from overshooting. 1.1

Therefore, these overshoots that are not present in the unrefined (1:1) case can be used as a measure for aliasing errors. The amount of the overshoot increases with decreasing refinement ratio. The fwst time the puff leaves the refined subdomain, the peak suddenly increases by 1.7%, 1.9%, and 2.4% for refinement ratios of 1:2, 1:4, and 1:8, respectively. Part of the increase in final peak heights compared to the unrefined case is due to these overshoots, but there is an apparent improvement due to the increased accuracy in the refined grid.

4

S U M M A R Y AND C O N C L U S I O N S

A finite element refinement technique is developed for locally increasing the resolution in air quality models. Mass conservation that is a primary requirement in air quality modeling is assured by maintaining the continuity of element basis functions. Noise waves that are generated at the refinement boundaries are partially eliminated by using a nonlinear filter. Here, several refinement ratios are tested in a rotating puff problem to evaluate the performance of element refinements in increasing the accuracy and the effectiveness of the filter in eliminating the aliasing errors. The accuracy of following the transport increases when element refinements are used. However, these refinements must be applied gradually (e.g., 1:2 or 1:4), because more abrupt refinements (e.g., 1:8) cause larger aliasing errors. The local nonlinear filter used here was effective in controlling the aliasing errors. The element refinement technique described here may be beneficial in regional air quality models where there is a need for increasing the model's resolution near emission sources.

Acknowledgements

The support for this project was provided by the U.S. Environmental Protection Agency, Division of University Grants; the Coordinating Research Council, the Pittsburgh Supercomputing Center and the California Air Resources Board. -

-

i

1:1 .............. 1:2 ........ 1:4 1:8

1.0

REFERENCES 1

2

°.

I=

65

0.9

3 0.8 4 0.7

' o

' I

' 2

5

Number of Rotations

Fig. 4 History of the peak heights in rotating puff tests with element refinement ratios of 1:1 (no refinement), 1:2, 1:4 and 1:8.

6

Odman, M.T. and Russell A.G. A multiscale finite element pollutant transport scheme for urban and regional

modeling. Atmos. Environ. 1991, 25A, 2385-2394 Rat, S.T., Sistla G., Ku J.Y., Schere K., Scheffe R. and Godowitch J. Nested grid modeling approach for assessing urban ozone air quality. Paper 89-42A.2 Air and Waste Management Association 1989, Pittsburgh, PA. Pleim, J.E., Chang J.S. and Zhang K. A nested grid mesoscale atmospheric chemistry model. Y. of Geophys. Res. 1991, 96, 3065-3084 Odman, M.T. and Russell A.G. Multiscale modeling of pollutant transport and chemistry. J. of Geophys. Res. 1991, 96, 7363-7370 Oden, J.T., Strouboulis T. and Devloo P. Adaptive finite element methods for the analysis of inviscid compressible flow: Part I. Fast reflnement/unreflnement and moving mesh methods for unstructured meshes. Coraput. Meth. Appl. Mech. Engng. 1986, 59, 327-362 Young, D.P., Melvin R.G. and Bieterman M.B. A locally

M. T. Odman, A. G. Russell

66

refined rectangular grid finite element method: Application to computational fluid dynamics and computational physics. .L of Comput. Phys. 1991, 92, 1-66

Zienkiewicz, O.C. and R.L. Taylor The finite element method, Vol. I, Basic formulations and linear problems. McGraw-Hill 1988, London, 290-318 Berger, M.J. On conservation at grid interfaces. SIAM J. Numer. Anal. 1987, 24, 967-984

9

Alpert, P. Implicit filtering in conjunction with explicit filtering. ,L of Comput. Phys. 1981, 44, 212-219 10 Koch, S.E and McQueen J.T. A survey of nested grid techniques and their potential for use within the MASS weather prediction model. NASA Technical Memorandum 87808, 1987 1 1 0 d m a n , M.T. and Russell A.G. A nonlinear filtering algorithm for multi-dimensional finite element pollutant advection schemes. Atmos. Environ. 1993, in press