A Combined Probabilistic and Optimization Approach for Improved Chemical Mixing Systems Design ()
1. Introduction
Computational fluid dynamics (CFD) analyses, with software packages like Fluent™ (ANSYS, Inc.), have become accepted techniques for solving complex fluid flows using numerical implementation of fluid mechanics principles. CFD-driven optimization has been explored by several researchers [1-3]. Peigin et al. has proposed a design tool that implements multi-constrained optimization of shape design driven by Genetic Algorithms (GA) coupled with CFD. The benefit to GA is that they can handle a large number (20+) of design variables; however, with a large number of constraints, GA are computationally expensive. Peigin et al. reports 15-18 hours for a single point optimization and on average 8-12 optimization steps [1]. A similar multi-constrained optimization strategy has been implemented on a ship’s hull. Again the downside is the time constraints, where they report approximately 10 days for 20 shape generations on a PC-based cluster [2]. Many researchers are using these GA’s in order to optimize a system with numerous variables.
GA’s have also been used by Carroll (1996) for the modeling of complex chemical mixing systems. The GA was coupled with the Blaze II code [4]. The Blaze II code can address up to 500 chemical reactions and 40 species. It also contains 1-D fluid dynamic equations, with mixing terms derived from 2-D equations, that can be used to model axis symmetric and 2-D flow fields [4]. The GA used 5 parameters with either 32 or 16 discrete possibilities resulting in approximately 2 million permutations. Even with a 2-D model and only 5 variable parameters the GA required 8 days of continuous runtime [3].
The objective of the current project is to develop and implement a design analysis of a fluid mixing nozzle using a coupled optimization and probabilistic approach. Due to the potentially large computational times to run GA optimization, our approach is to reduce the number of constraints though probabilistic analysis in order to use single objective optimization techniques. This technique was applied to mixing nozzles for a generic chemical mixing system. Based on the insight into the features that most affect performance provided by the probabilistic analysis, optimization has the potential to improve design performance while also reducing the weight of the system. Additionally, a probabilistic analysis of the optimized design can confirm the original optimization parameters and provide insight into the effects of manufacturing tolerances of specific geometric variables. This will in turn, reduce manufacturing costs by determining which dimensions are important to the mixing performance of the nozzle and which geometric tolerances can be loosened. A clear benefit to this design approach is that it allows for fast optimization by reducing the number of constraints through probabilistic analysis, as well as, assessing the impact of manufacturing tolerances.
2. Methodology
Currently, there is no commercially available tool in industry to perform fluid mechanics analyses, optimization and probabilistic analysis in the proposed integrated manner. There are however, software packages that can execute individual components of the process. Our approach is to use commercially available CFD, probabilistic and optimization software and interface them together with custom scripting. The probabilistic software, or optimization routine, manages the CFD code by importing any number of variables within ranges and/or distributions set by the user. Figure 1 and 2 show a diagram linking the processes mentioned above.
2.1. CFD Model
The CFD software used for the simulations was Fluent™ Version 6.3.26 [5]. It is a widely used computational software package for modeling fluid flow and heat transfer in complex geometries. Easy mesh generation and ability to refine or coarsen the mesh autonomously based on the flow solution are just some of the features that make this CFD package extremely versatile and ideal for automation. Gambit ® was the pre-processor used for the solid modeling and mesh generation.
The algorithm employed was the pressure-based Navier-Stokes solution algorithm. Typically, this algorithm is used for low velocity incompressible flows, which fits this case where the flow velocity is about 83 m/s. The momentum equation provides the velocity field and the density is calculated from the equation of state. Other governing equations include energy and species conservation. The turbulence model chosen was the standard k-epsilon turbulence model. This model is robust and suitable for initial iterations, initial alternative design screenings, and parametric studies. The k-epsilon model will be ideal for the automated analysis where many different shapes will be analyzed. A no-slip boundary condition was placed at the walls. The primary inlet utilized a user-defined function (UDF) that modeled a fully developed fluid flow.
Definition of the boundary conditions is critical to both accuracy of the model and computational efficiency, the importance of which should not be under estimated when venturing into optimization and probabilistic analysis. Considerable attention was paid to identifying the best algorithmic parameters in order to achieve convergence in the minimum number of iterations, reliably, and with physically correct solutions regardless of geometry. It is also noteworthy to mention that this process was
Figure 1. Diagram of the probabilistic/optimization and CFD interface. Each program is linked by custom scripts.
Figure 2. Example flow chart of the probabilistic interface using the AMV method.
carried out in parallel on 8 processors in order to reduce computational time.
2.2. Probabilistic Model
As uncertainty is inherent in physical systems, a probabilistic analysis model input variables as distributions and then predicts a distribution of performance. Based on the distribution of performance, sensitivity and importance factors are used to identify critical parameters. The probabilistic software, Nessus®, implements a variety of probabilistic methods that vary in efficiency and accuracy of the solution. The most commonly used probabilistic method is Monte Carlo [6]. The Monte Carlo method generates random values for each variable according to its distribution and then predicts the distribution of performance through repeated trials. As the accuracy of the prediction is dependent on the number of trials performed, the Monte Carlo method is computationally expensive. The mean-value (MV) family of methods are approximate, but considerably more efficient than the Monte Carlo method. They create a mean-based response function and compute the most probable point (MPP), which is the shortest distance from the origin to the limit state surface and represents the combination of stochastic variables resulting in a specific level of performance [6]. In this study, the Advanced Mean Value (AMV) method was applied, which uses a higher order approximation to determine the MPP [7]. The number of trials used for the AMV method is 1 + the number of variables + the number of probability levels desired [7]. While the AMV method is a discrete and approximate solution, it has shown excellent agreement with Monte Carlo analysis for monotonic system. This has been shown in other applications which utilize FE analysis with geometric perturbations and realistic loading conditions [8,9]. The major advantage of the AMV is that it requires a small number of trials, which saves significant computational time when the CFD analysis requires multiple iterations to reach convergence. Both of these probabilistic methods provide sensitivity factors identifying and quantifying the contributions of each variable on performance of the system.
The AMV method was utilized early in the design process to identify the variables that contributed significantly to the behavior of the system reducing the number of parameters required for performance optimization. By reducing the number of variables in the optimization, the computation time and required resources are dramatically reduced. Instead of using a sensitivity study perturbing a single variable at a time, a probabilistic approach was used to allow for the interaction affects between the various input parameters. During the early stages of the design process, the emphasis was on identifying the important and unimportant variables. Following optimization, AMV methods were employed to evaluate the effects of manufacturing tolerances on system perform-ance.
2.3. Optimization Model
Numerical optimization techniques are designed to minimize an objective function subject to constraints, with many algorithms developed over the past several decades [10]. In general, the algorithms require a starting point, x0, and then iterate or step until there is no more progression, or the approximate solution falls within a user-defined tolerance. Typically, algorithms follow one of two types of strategies, line search or trust region. This study implemented a trust region [11] strategy in order to account for geometric changes that may result in the fluid domain acting non-linearly. A common problem in line searches is the fixed step size can causes them to miss a local minimum, where as the step size in the trust region search is not fixed and, therefore, has a better opportunity to find a minimum that is close to the current point.
The success and efficiency of an optimization is contingent on selection of an appropriate algorithm and an accuracy characterization of the problem. It was not known whether the variables would behave linear, nonlinear, or convex, only that they could hold any value between the bounds. The optimization algorithm had to be suitable for a continuous objective function with variables that are constrained by simple bounds and can solve for linear, non linear and convex variables.
A trust-region algorithm was implemented, which utilizes an active-set algorithm, for the optimiz-ation analysis. An active-set algorithm will employ linear techniques to estimate the active-set at each iteration and then solve an equality constrained quadratic program to generate a step [12]. This method was used because it tends to yield more exact solutions and is less sensitive to the initial starting point than interior point methods. Another benefit, in our case, of the active-set algorithm is that it uses a gradient projection method when only bounds are applied to the constraints [13]. The gradient projection method attempts to speed up the solution process within the active-set, but is only utilized when the variables are bounded. It consists of two different stages. First, the search direction will be along the path of steepest decent from the current point. The second stage investigates the face of the feasible region using the active-set constraints [12]. This can significantly reduce the optimization time.
2.4. Interfacing Model
To facilitate communication between the all of the software packages, custom interfacing was developed to build CFD models with perturbed parameters and calculate performance parameters from the analysis outputs. Interfacing was performed with components written in Matlab®, Dos and C. In addition, checks were performed to ensure mesh quality to prevent analyses that would fail or highly skewed elements, which may lead to convergence issues. This is noteworthy because the automated process can potentially take days and even weeks to run and computational efficiency will be a driving factor, especially as more complex flows are examined. An imported UDF specifies the physical boundary conditions for the specific system. Once the CFD simulation converges to a solution for the flow field; it calculates the values of a user defined fluid property (i.e. mole fraction of interested species) via a UDF and prints the results to a file. A script is also utilized to calculate the performance parameter and print the results for analysis by either the probabilistic or optimization routines. The performance parameter used is Mixedness, which is defined by:
(1)
The degree of mixing is measured by the ratio of the integral value for species mole fraction (Mf) across an exit plane divided by the homogeneous mole fraction (Mf_Homogeneous), where n is the number of nodes within the exit plane. Increased mixing of species in chemical systems should result in greater chemical efficiencies and better performance. This interfacing routine is continued until all of the probabilistic or optimization trials are completed.
3. Problem Description
3.1. CFD Model
The combined probabilistic and optimization approach is demonstrated for a low flow, subsonic Hydrogen-Iodide Chlorine (HICl) Laser. As opposed to other chemical lasers the HICl Laser has yet to reach its expected performance potentially due in part to a lack of homogeneous distribution of its excited chemical species.
Additionally:
1) the subsonic HICl Laser has simple mixing geometry, employing a repetitive nozzle array that allows the computational domain to be reduced based on simple lines of symmetry;
2) the cross-flow injection geometry would be easy to perturb within the code;
3) validation data exists for injection cross-flow at speeds similar to the HICl.
The geometry for the subsonic HICl Laser consists of a rectangular flow channel (Figure 3(a)). There are four rows of secondary flow injection nozzles on the top and bottom plates of the cavity (Figure 3(b)). The primary inlet flow for the subsonic HICl Laser is a mixture of helium (He) and hydrogen (H2). The secondary flow through the nozzle plates has two different mixtures that are being injected into the primary flow. The first three rows inject a mixture of helium and hydrogen-iodide (HI), where the fourth injects a mixture of helium and nitrogen-trichloride (NCl3). Table 1 shows the boundary conditions applied during this design optimization run. A constant cross-sectional area was maintained on each inlet, but the aspect ratio was allowed to vary.
It would be computationally expensive to model the entire flow channel for the CFD analysis; con-sequently, planes of symmetry are used again to reduce the size of
Figure 3. (a) A section view of the laser cavity showing the secondary inlet nozzle plates, and (b) detailed view of one of the nozzle plates. Note: the system has two identical sets of nozzle plates.
Table 1. Flow conditions for subsonic HICl laser.
the modeled fluid domain (Figure 4). The fluid domain was modeled and meshed using text commands within a journal file. This allowed for easy modification during the automated process.
The algorithms, boundary conditions, and meshing strategies for the CFD simulation were validated using experimental data from subsonic (yet compressible flows) jets in cross flow [14]. Downstream velocity and temperature profiles obtained from the CFD simulations showed good comparison (Figure 5) to the experimental data of Dizene et al. (2000). This effort highlighted the need to accurately describe the velocity profile at the
Figure 4. Schematic of the variables and geometry for the first probabi-listic analysis.
Figure 5. Velocity and Temperature profiles from experimental data compared to CFD simulations of jet in crossflow at X/D = 4.
entrance boundaries and confirmed that the dynamic adaptive grid technique and algorithms selected worked well in these types of fluid flow conditions.
3.2. Probabilistic Model
The initial probabilistic CFD analysis modeled eight geometric variables as distributions. The parameters were the size of the holes in two directions, the offset between the holes, and the hole placement, as illustrated in Figure 4. All of the parameters were modeled as normal distributions defined by mean and standard deviation (Table 2). Standard deviations were determined by the maximum range permissible by geometric limits.
3.3. Optimization Model
Optimization was then applied to determine the values of the variables within the constraints that maximized Mixedness. The four important variables, LHRx, LHRz, SHRx, and LHRz can be combined into two by keeping the cross sectional areas of the orifices constant. Therefore, only two variables, LHRz and SHRz, need to be perturbed for the HICl Laser optimization. The objective function and constraints for the optimization is stated below:
Maximize: M = f(LHRz, SHRz) (2)
Subject to: 0.0055 cm ≤ LHRz ≤ 0.023 cm
0.010 cm ≤ SHRz ≤ 0.042 cm (3)
The optimization routine was used to determine the values of LHRz and SHRz in each analysis. The initial guess in the optimization algorithm were the mean values from the probabilistic analysis. The automated process updated the geometry, generated the mesh, performed the CFD analysis and extracted the results to compute the Mixedness parameter. Since the maximum Mixedness is desired, the minimization of the ratio of HI concen-
Table 2. Probabilistic variables with respective mean, standard deviation, and distribution.
trations to the maximum concentration (the second term in Equation (1) was calculated. The optimization process continues until the Mixedness value has converged to within a 0.01% tolerance. In addition, initial guesses near the lower and upper bounds of the design space were also evaluated to check for local minimums that may be on either side of the mean orifice dimensions.
4. Results and Discussion
4.1. Identification of Important Parameters
Figure 6 shows the importance levels associated with each variable that was perturbed in the modeled geometry. The importance factors represent the design vector or value of each parameter that defines the MPP, which is proportionate to the output measure at the specified probability. As they are reported in the standard normal space, importance factors are relative measures and the sum of the squares for each measure will equal 1. The importance levels identify the variables that contribute the most to the reliability of the design; therefore, it is deduced from Figure 5 that the radii of the inlet orifices contribute the most to either increase or decrease of the response variable, Mixedness. The other important observation is that PRCNT1, PRCNT2, OFFSETH and L_2NDSET had little to no effect on the Mixedness of the system. Recall, they had the largest standard deviation of all the variables. Now, the elimination of four of the eight variables within the geometry allows for the focus on only the radii of the orifices.