(S)

CrossMark

Available online at www.sciencedirect.com

ScienceDirect

Energy Procedía 69 (2015) 148 - 157

International Conference on Concentrating Solar Power and Chemical Energy Systems,

SolarPACES 2014

Computational flow optimization of heliostat aspect ratio for wind

direction and elevation angle

M.D. Maraisa, K.J. Craigb*, J.P. Meyerc

aGraduate student, Department of Mechanical and Aeronautical Engineering, University of Pretoria, Pretoria 0002, South Africa bPrEng, PhD, Professor, Department of Mechanical and Aeronautical Engineering, University of Pretoria, Pretoria 0002, South Africa cPrEng, PhD, Professor, HOD, Department of Mechanical and Aeronautical Engineering, University of Pretoria, Pretoria 0002, South Africa

Abstract

The choice of heliostat aspect ratio has previously been shown to have an effect on the wind loading coefficients experienced by the structure due to an attacking atmospheric wind. Knowing the wind loads on heliostats is important in order to allow for adequate sizing of all components. The preferred method for determining such loadings is experimental wind tunnel studies, with Computational Fluid Dynamics (CFD) being less popular. A CFD model capable of predicting the mean wind loadings on a heliostat structure subjected to an attacking atmospheric wind and turbulence intensity profile was developed and validated using experimental results available in the open literature. The commercial software package ANSYS Fluent V15.0 was used together with the ANSYS Workbench Design Explorer toolset to study the influence of aspect ratio on the wind loading coefficients.

©2015The Authors.Publishedby ElsevierLtd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.Org/licenses/by-nc-nd/4.0/).

Peer review by the scientific conference committee of SolarPACES 2014 under responsibility of PSE AG

Keywords: heliostat; optimization; aspect ratio; concentrated solar power; solar tower; wind load; computational fluid dynamics

1. Introduction

A heliostat field makes up 40% to 50% of the total capital cost of a solar tower power plant [1]. It is therefore essential that methods be investigated to design the field in an optimal manner, considering both the layout of heliostats to improve optical efficiency and the design of individual structures to be as light as possible, in an

* Corresponding author. Tel.: +27-12-420-3515. E-mail address: ken.craig@up.ac.za

1876-6102 © 2015 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.Org/licenses/by-nc-nd/4.0/).

Peer review by the scientific conference committee of SolarPACES 2014 under responsibility of PSE AG doi:10.1016/j.egypro.2015.03.018

attempt to reduce the relatively high cost of such plants. Recent research efforts directed towards reducing the wind loadings on a light-weight heliostat structure are indicative of the trend towards designing ever lighter structures that are still able to resist atmospheric wind loadings [2]. The choice of heliostat aspect ratio has also been shown to have an effect on the wind loading coefficients experienced by a heliostat structure due to an attacking wind [3]. Strategies for determining suitable heliostat aspect ratios other than exhaustive experimental studies are not yet widely available.

Investigations into heliostat aerodynamics are usually done experimentally in wind tunnels [3, 4] with CFD investigations being less popular. A three-dimensional CFD analysis of flow over a heliostat structure is available in the open literature [5]. Although the results did not compare favorably to wind tunnel tests, the study came to the conclusion that CFD can be a valuable tool in heliostat design and optimization. The ability of CFD to accurately predict wind loadings over solar collectors was proven recently when commercial software based on the Lattice-Boltzmann method was used to predict the wind loadings on parabolic trough collectors to within 10% of those obtained experimentally [6].

Nomenclature

Aref reflector area (m2) Reh Reynolds number (Reh = pUr<fhl^i)

b reflector width (m) t time (s)

CF force coefficient U(z) mean stream-wise velocity component (mls)

CM moment coefficient u*ABL ABL friction velocity (mls)

Cs roughness constant u fluctuating velocity component (mls)

Cel, Ce2 turbulence model constants Üi time averaged velocity component (m/s)

C„ ak turbulence model constants Uref reference velocity at height zrf (mls)

F wind force (N) x, y, z Cartesian coordinates

H height of elevation axis (m) Zo aerodynamic roughness length (m)

h reflector height (m) zref reference height (m)

ij indication of x, y or z a heliostat elevation angle (°)

lu longitudinal turbulence intensity ß wind azimuth angle (°)

K turbulent kinetic energy (m2/s2) Kronecker delta

ks equivalent sand-grain roughness height (m) s turbulence dissipation rate (m2ls3)

ks,ABL ABL sand-grain roughness height (m) KK-e von Karman constant

Lref reference length, h at ra = 1 (m) Ml molecular viscosity of air (Pa.s)

M moment caused by wind (Nm) Mt turbulent viscosity (Pa.s)

n velocity power-law exponent P density of air (kglm3)

V time averaged pressure (Pa) Xj stress (Pa)

ra aspect ratio (ra = b/h)

The limited amount of CFD studies and few available strategies for determining optimal aspect ratios led to the objectives of this study being twofold: 1) Developing a robust CFD model capable of predicting the mean wind loadings on a heliostat structure subjected to an attacking atmospheric wind and turbulence intensity profile and 2) using the CFD model, together with an optimization technique, to find an optimum aspect ratio. The coordinate system and characteristic lengths that define the heliostat model and aerodynamic loadings used by Peterka and Derickson [7] have also been used in a previous paper concerning wind tunnel tests of heliostat models with different aspect ratios [3]. In keeping with the convention, the same coordinate system is used here and is presented in Fig. 1.

Fig. 1. Coordinate system and characteristic lengths as used by Peterka and Derickson [7].

2. Problem formulation

The optimum aspect ratio, when considering wind loadings, depends on the objective function. A heliostat with altitude-azimuth elevation tracks the sun by rotating its reflector about two axes, namely the z- and y-axes with reference to Fig. 1. In order to control this movement the heliostat utilizes two drive mechanisms, one for each rotation. These drives have to overcome the static moments induced by gravity loading as well as the moment coefficients CMz and CMHy caused by wind loadings. The effects that wind loadings will have on drive mechanisms will become all the more important as the mass, and consequent gravity loadings, of heliostat structures are reduced by refined design procedures. The drives have been estimated to constitute 30% of the total heliostat cost and drive deformation under wind loadings have been shown to cause spillage losses as the mirror facet departs from its ideal, aligned position [8]. The two drive moment coefficients have therefore been selected as objective functions.

3. Method

3.1. Governing equations

The K-e turbulence model used in this paper is based on the Reynolds Averaged Navier Stokes (RANS) equations. For an atmospheric boundary layer flow it is reasonable to treat air as an incompressible gas and therefore the momentum equation may be written as in Eq. 1.

Two additional conservation equations are introduced for the standard K-e turbulence model, one for the turbulent kinetic energy K and one for its rate of dissipation e:

with the supplementary relationships in Eq. 4 providing closure.

3 pKSij and y.t = C^p —

3.2. The problem of horizontal homogeneity of the atmospheric boundary layer

Since a heliostat structure is exposed to an atmospheric wind, the correct modeling of the velocity and turbulence profiles are just as important as accurately modeling the geometry of the structure and associated boundary conditions. An accurate CFD model of an atmospheric boundary layer (ABL) flow over a structure requires appropriate boundary conditions. At the inlet, profiles of mean velocity and turbulence quantities that are fully developed and in equilibrium should be specified. In order to reduce computational cost it is not feasible to model the infinitely large terrain upstream of the structure, but rather to introduce the profiles a short distance upstream. To ensure that the profiles are representative of the terrain roughness, an appropriate aerodynamic roughness length z0 is used in the case of the log-law, or a suitable power-law exponent n.

Blocken, et al. [9] have identified three important regions in an atmospheric computational domain, as well as three important profiles, named according to their stream-wise position relative to the model of interest. An inlet flow is the boundary condition specified by the modeler. The approach flow profiles are those travelling towards the model whilst the incident flow profiles are those that would be obtained in a similar but empty domain at the position where the model will be located. In order for the ABL to be horizontally homogenous (HH), all three these flow profiles must have the same shape in an equivalent, but empty (without the structure or building models), computational domain.

Several authors working in the computational wind engineering field have reported results where the inlet profiles changed rapidly after entering the computational domain, effectively causing the incident profiles to be very different from those specified at the inlet and making it impossible to exactly control the type of flow attacking the structure. These problems are mainly characterized by the acceleration of the flow close to the ground plane as well as dissipation of near-wall turbulence as the flow develops. According to Blocken, et al. [9] the unintended differences between inlet profiles and incident profiles can severely hamper the success of CFD simulations given that minor changes to the incident flow profiles can cause significant changes in the flow field.

Richards and Hoxey [10] addressed the appropriate boundary conditions for computational wind engineering models using the K-e turbulence model and state that the existence of a homogenous flow in an empty computational domain has three main implications: 1) The vertical velocity is zero, 2) pressure is constant (i.e. no pressure gradients exist) and 3) shear stress is constant throughout the domain. From this it follows that a constant friction velocity exists throughout the domain. It is also necessary to ensure the equilibrium of the inlet profiles, ground shear stress and turbulence model. In a more recent paper Richard and Norris [11] show that the appropriate profiles for mean velocity and turbulence properties may be directly derived from Eq. 1 to 4 for the K-e turbulence model and also follow the same process for other models. The profiles appear in Eq. 5 and it must be stressed that in an attempt to model a HH boundary layer, no other form for the profiles are acceptable, even the frequently used power-law velocity profile is not recommended [11].

where the friction velocity u*ABL may be calculated using a reference velocity at some reference height by replacing the U and z with reference values. By replacing u*ABL in Eq. 5 with its new expression in terms of reference values, the velocity profile may be expressed as in Eq. 6.

ln (*)

U(Z) = ^/"T^A (6)

ln ("f)

The recommendation [11] continues by stating that the von Karman constant is effectively determined by the turbulence model and is given by Eq. 7.

KK-E =

(C£l - Ce2)ae |c„ = 0.433 (7)

A last requirement is that the top boundary condition of the domain must include a driving shear stress, zero-flux condition for K, and the flux of e across the boundary. In this paper the top boundary condition, as well as the two side boundary conditions, were slip-walls, since a wind tunnel study was being replicated instead of true external ABL flow.

When using the standard K-e turbulence model with standard wall functions, the near-wall flow is not resolved but rather wall functions are used to solve for the unknowns. Blocken, et al. [9] showed that the inlet profile, based on some aerodynamic roughness length z0, can be partly maintained by applying an equivalent sand grain roughness ks to the bottom wall. The relationship between z0 and ks is specific to the CFD package and its implementation of wall functions. For ANSYS Fluent the following relationship is given [9], where Cs is a roughness constant with a value of 0.5:

_ 9.793z0

ks,ABL _ "p (8)

It may be argued that the problem of an appropriate incident profile can be mitigated by simply placing the inlet plane very close to the model. This can be regarded as one remedy [9], but experimentation will be required to ensure that the inlet is not so close that it incorrectly prevents the flow from reacting naturally to the presence of the model. Ensuring that the ABL being modeled is homogenous gives the modeler the freedom to allow sufficient upstream distance whilst being assured that the inlet and incident profiles are the same.

3.3. Computational domain and boundary conditions

Two different heliostat models were used for this study, one that exactly matched the dimensions of the 1:60 scale model used by Peterka, et al. [4], including the three facets, pylon and torque tube, to validate the CFD model. Dimensions were obtained from an engineering diagram in the stated reference. Dimensions not explicitly given like the diameters of the torque tube and pylon, as well as the reflector thickness, were measured from the drawing.

Previous numerical studies have shown that the gap between mirror facets have an effect on wind loadings [12] and since it was unclear how the gaps would scale with aspect ratio ra, it was decided to use a heliostat model with a single facet to observe the effects of varying the ra in the current study. The pylon and torque tube were kept the same diameter but their lengths did scale with aspect ratio. The varying aspect ratio heliostat had a fixed reflector area of Aref = 0.0117 m2. In a previous study [3], the gap between the facet and the ground for a = 90° was kept constant. For this study a gap size of 0.015 m was used, the same as that of a heliostat with ra = 1 [4]. This has the effect of causing the elevation axis height H to be an implicit function of aspect ratio through H = h/2 + 0.015 where h = (Are/r/'5.

The mesh is a hybrid quad/tri mesh that automates the re-meshing process as the heliostat changes angle and aspect ratio, consisting of 6.8 million cells. ANSYSFluent V15.0 was used to solve the steady-state flow field using the K-e turbulence model with standard wall functions. The outlet boundary condition was a zero gauge pressure outlet and convergence was monitored through scaled residuals, monitoring the two objective functions and checking for mass balance between the inlet and outlet of the domains. The coupled pressure-based segregated solver was used. Fig. 2 shows the computational domain and mesh for the validation case.

Fig. 2. (a) Computational domain divided into blocks; (b) unstructured block with validation case heliostat model.

Obtaining the inlet profiles required experimental data so that the CFD results could be validated and for this purpose the velocity and turbulence intensity profiles of Peterka, et al. [4] were chosen. Eq. 6 was fitted to the experimental data by solving a least squares minimization problem to find an optimum aerodynamic roughness length.

It is therefore noted that for a chosen reference height and reference velocity, the only variable that can be tuned to match a certain velocity profile is z0. For this study, the reference height (zref = 0.168 m) was chosen to be the same as that of Peterka, et al. [4]. Reference velocity (Uref = 7.82 m/s) was chosen to ensure a Reynolds number of above 30 000, high enough to ensure Reynolds number independency of the aerodynamic coefficients [4]. The second profile to be matched was that of turbulence intensity, but we have to specify inlets for K and e using Eq. 5. From the definition of K and the assumption of isotropic turbulence, we may relate the turbulence intensity Iu to K:

U(z) U(z)

From Eq. 9 we therefore conclude that, short of modifying turbulence model constants like CM, we have no real control over the profile of Iu after the velocity profile has been defined since K is defined by Eq. 5. The turbulence intensity profile is therefore a direct consequence of the velocity profile and for each individual velocity profile only one matching Iu profile exists if we want to ensure horizontal homogeneity. This seems to be one of the major shortcomings of using RANS turbulence models to model ABL flow and the subsequent wind loadings of structures.

Fig. 3.a shows the profiles used for this study (the profile locations are indicated in Fig. 2. a) and how they compare with those of Peterka, et al. [4]. Incident CFD profiles are shown to ensure that data compared to experimental values are indeed those observed at the model location. It can be seen that Iu is under-predicted in the CFD model when compared to the experiment. The incident profiles in Fig. 3.a are as a result of applying a sand-grain roughness to the bottom wall of the domain. Fig. 3.b shows the improvement obtained by applying an equivalent sand-grain roughness to the bottom wall of an empty computational domain using Eq. 8.

Fig. 3. (a) CFD incident profiles compared to those of [4]; (b) improvement in horizontal homogeneity obtained by roughening the bottom wall.

3.4. Optimization procedure

The optimization procedure makes use of both a multi-objective (MO) maximization step to obtain a worst-case orientation and an MO minimization step to find an optimal aspect ratio at said orientation. This nested optimization approach has been used previously to minimize the effect of automotive pollution in an urban environment [13], albeit with a single objective function. The design variables are the two orientation angles (x1 _ [a fi]T) and the aspect ratio of the reflector (x2 _ [ra]) where the aspect ratio is given by ra _ b/h. The optimization may then be described as minimizewrtx2{maximum.wrtxif(x1, x2)} with fix1, x2) _/(/iC*1, x2) _ |cMHy|, f2(x1, x2) _

|CMz|j as the multi-objective function. Response surfaces are created from the solver results using the ANSYS Workbench Design Explorer toolset and Kriging method. The built-in Multi-Objective Genetic Algorithm (MOGA) is used to optimize the response surfaces and obtain candidate optima.

4. Results

4.1. Validation

In order to judge the accuracy of the CFD simulations, results were compared to those of Peterka, et al. [4] for a _ 30° and a p sweep from 0° to 180° in increments of 45°. A comparison is made in Fig. 4 for two force coefficients and two moment coefficients. It may be seen that the turbulence model tends to underpredict the magnitude of the forces on the heliostat model but does predict the correct trend for the aerodynamic coefficients. The hinge moment coefficient is also underpredicted, but the CFD model closely follows the trend of the experimental values for all values of p, an important qualitative aspect that should ensure reasonable optimization results that rely on the shapes of response surfaces and the consequent locations of global and local minima and maxima. The moment about the z-axis is reasonably predicted up until p _ 90°, after which the values do not agree well. It is however hard to explain why a symmetrical structure such as a heliostat would have a non-zero CMz as measured and reported in the experimental study [4] when the wind is approaching from directly behind (¡3 _ 180°) and also a zero moment when the wind is approaching from some angle incident to the reflector, as was the case with the experiment for /? _ 135°. A possible reason for the general under prediction might be the fact that the turbulence intensity is lower in the CFD simulation than in the experiment (refer to Fig. 3.a). The force and moment coefficients are defined in Eq. 10.

Cm —

015 0.1 oos o

-O.ÛS -m

OÉXp.CMI üE*p,CMHy • CFD, CM; ACFO, CMIIy

Fig. 4. Comparison of CFD results with experimental values from Peterka, et al. [4] for a = 30° and varying {1.

or CP — ■

0.5pU'refArefLref 0.5pUrefAref

The results indicate that the magnitudes of the force and moment coefficients are close to symmetrical about P = 90°, indicating that the effect of the reflector dominates the small contribution due to the exposed torque tube and pylon when the wind approaches from the back. For this reason the following investigations regarding aspect ratio were done for 0° < P < 90° to limit the amount of simulations and hence the computational cost.

4.2. Optimization results

For the MO maximization step, combinations of a and y3 were used to create design points. Angles were varied from 0° to 90° in increments of 22.5° for a heliostat with ra = 1, requiring a total of 25 CFD simulations. The objective functions were then maximized to obtain candidate worst-case orientations. Contour plots of the resulting response surfaces are shown in Fig. 5 together with the candidate maxima.

Fig. 5. Contour plots of the moment coefficient Kriging responses for a heliostat with ra = 1, candidate maxima indicated by black dots.

After choosing a worst-case orientation (a = 43.9°, ft = 40.5°) from the maximization results, ra was varied between 0.5 and 3 to investigate the effect of aspect ratio at that orientation. CMz increases with an increase in aspect ratio whilst CMHy decreases (see Fig. 6.a). This is true for the most part because the moments are directly proportional to their respective moment arms as measured between the centre of the reflector and the point where the net aerodynamic force vector acts. The arm increases in the y-direction and decreases in the z-direction as ra increases. There is however a plateau near ra = 1 where CMz remains almost constant as aspect ratio is increased (see Fig. 6.b). This is due the existence of the atmospheric boundary layer and although the moment arm increases, the

overall height of the heliostat decreases, causing it to experience lower velocities and thus lower loadings near the ground plane. A local minimum is reached in the region of ra = 1.1 and the MO minimization step predicts this by suggesting a new aspect ratio of 1.1.

b u ï

E Ï С

§ 0.05 £

— - CM Ну Kriging line ■ CM H Y

-СМ/ Kriging line

A CM г

Fig. 6. (a) Moment coefficient responses for the range of ra tested at a = 43.9°, p = 40.5°; (b) Kriging responses near ra = 1

The new aspect ratio of 1.1 has the effect of decreasing the maximum CMHy by 7.5%, although increasing the maximum CMz by 5.8%. The z-moment is therefore more sensitive to an increase in moment arm at its maximum orientation than to a reduction in wind load by being in closer proximity to the ground. For the angles considered it seems as if reducing the one moment coefficient by varying the aspect ratio will increase the other, requiring a tradeoff. It could be argued that a = 90°, y3 = 67.5° is not an operational angle that a heliostat will typically achieve during its tracking of the sun and future work could include limiting the range of angles to operational values for a specific solar tower site. Larger reductions are however obtainable for CMHy for larger aspect ratios in the region of ra = 1 than the resulting increase in Cmj. Lift and drag coefficients also favour larger aspect ratios due to lower velocities near the ground and will influence the optimization accordingly if included in the objective functions.

4.3. CFD results

Fig. 7 shows contour plots of static pressure on and around the heliostat. The largest aerodynamic forces at the considered orientation are pressure forces with shear forces making a smaller contribution. Fig. 7.a shows the pressure difference between the up- and downstream faces of the reflector, with the largest difference occurring at the leading edge, resulting in a moment about the ;y-axis that wants to tilt the heliostat in the direction shown. In Fig. 7.b the reason for the existence of a moment about the vertical axis is also clear since one side of the reflector shows higher static pressure values than the other, causing a moment that wants to pivot the heliostat about its pylon.

Fig. 7. (a) CFD gauge pressure contours (Pa) on a vertical plane for ra = 1.1 at a = 43.9°, ß = 40.5°; (b) pressure contours on heliostat surface.

5. Conclusion and recommendations

In this paper a method for optimizing the aspect ratio of a heliostat has been presented, based on minimizing the magnitudes of both the drive moment coefficients. Both objective functions were weighted equally and in reality the objective functions may require more complex considerations like, for instance, cost and manufacturability of components. In such cases different weighting factors may be used and additional force and moment coefficients could be included in the analysis.

It is also of major importance that a method be developed for predicting peak loadings on the structure instead of only mean values. For this purpose it would be necessary to run transient simulations using more advanced turbulence models such as Large Eddy Simulation (LES) or Detached Eddy Simulation (DES). The optimization procedure developed here could then be employed by using peak loadings as the objective functions.

A major drawback of Reynolds Averaged Navier Stokes (RANS) simulations concerning atmospheric boundary flow is that the velocity and turbulence intensity profiles cannot be decoupled if horizontal homogeneity is to be achieved. Turbulence intensities encountered in wind tunnel studies and at full scale sites seem to be higher than can be obtained by using Eq. 9. If CFD models using the RANS approach are to be properly validated, then velocity and turbulence profiles generated in experimental studies should attempt to replicate those required to model horizontally homogenous boundary layers.

Acknowledgements

The authors would like to acknowledge the support from the University of Pretoria (South Africa) and the South African National Research Foundation (DST-NRF Solar Spoke).

References

[1] Kolb, G.J., Jones, S.A., Donnelly, M.W., Gorman, D., Thomas, R., Davenport, R., and Lumia, R. Heliostat Cost Reduction Study, SAND2007-3293, Sandia, New Mexico, 2007.

[2] Pfahl, A., Brucks, A., and Holze, C. Wind Load Reduction for Light-Weight Heliostats. Proc. SolarPACES 2013 conference. Las Vegas, Nevada, USA. Published by Elsevier: Energy Procedia, 2014. 49: p. 193-200.

[3] Pfahl, A., Buselmeier, M., and Zaschke, M., Wind loads on heliostats and photovoltaic trackers of various aspect ratios. Solar Energy, 2011. 85(9): p. 2185-2201.

[4] Peterka, J.A., Hosoya, N., Bienkiewicz, B., and Cermak, J.E. Wind load reduction for heliostats, Report SERI/STR-253-2859, Solar Energy Research Institute, Golden, Colorado, USA, 1986.

[5] Wu, Z. and Wang, Z., Numerical study of wind load on heliostat. Progress in Computational Fluid Dynamics, an International Journal, 2008. 8(7): p. 503-509.

[6] Mier-Torrecilla, M., Herrera, E., and Doblaré, M. Numerical Calculation of Wind Loads over Solar Collectors. Proc. SolarPACES 2013 conference. Las Vegas, Nevada, USA. Published by Elsevier: Energy Procedia, 2014. 49: p. 163-173.

[7] Peterka, J.A. and Derickson, R.G. Wind load design methods for ground-based heliostats and parabolic dish collectors, Report SAND92-7009, Sandia National Laboratories, Springfield, Alabama, USA, 1992.

[8] Teufel, E., Buck, R., Pfahl, A., Boing, G., and Kunert, J. Dimensioning of heliostat components under wind andgraviy load: The map approach. Proc. SolarPACES 2008 conference. Las Vegas, Nevada, USA.

[9] Blocken, B., Stathopoulos, T., and Carmeliet, J., CFD simulation of the atmospheric boundary layer: wall function problems. Atmospheric environment, 2007. 41(2): p. 238-252.

[10] Richards, P.J. and Hoxey, R.P., Appropriate boundary conditions for computational wind engineering models using the k-e turbulence model. Journal of Wind Engineering and Industrial Aerodynamics, 1993. 46: p. 145-153.

[11] Richards, P.J. and Norris, S.E., Appropriate boundary conditions for computational wind engineering models revisited. Journal of Wind Engineering and Industrial Aerodynamics, 2011. 99(4): p. 257-266.

[12] Pfahl, A., Buselmeier, M., and Zaschke, M. Determination of wind loads on heliostats. Proc. SolarPACES 2011 conference. Granada, Spain.

[13] Craig, K.J., De Kock, D.J., and Snyman, J.A., Minimizing the effect of automotive pollution in urban geometry using mathematical optimization. Atmospheric environment, 2001. 35(3): p. 579-587.