CFD Based Shape Optimization of IC Engine

Intense competition and global regulations in the automotive industry has placed unprecedented demands on the performance, efficiency, and emissions of today's IC engines. The success or failure of a new engine design to meet these often-conflicting requirements is primarily dictated by its capability to provide minimal restriction for the inducted and exhausted flow and by its capability to generate strong large-scale in-cylinder motion. The first criterion is directly linked to power performance of the engine, while the latter has been shown to control the burn rate in IC engines. Enhanced burn rates are favorable to engine efficiency and partial load performance. CFD based numerical simulations have recently made it possible to study the development of such engine flows in great details. However, they offer little guidance for modifying the ports and chamber geometry controlling the flow to meet the desired performance. This paper presents a methodology which combines 3D, steady state CFD techniques with robust numerical optimization tools to design, rather than just evaluate the performance, of IC engine ports and chambers.


INTRODUCTION
The ability to control the flow in engine components has long been recognized to be the key controlling mechanism to enhance IC engine performance. In particular, control of the flow through the ports and the combustion chamber has long become critical for meeting the more stringent emission regulations and CAFE fuel economy requirements.
Since engine flows are strictly controlled by ports and chamber geometry [1][2][3][4], an analytically based shape optimization scheme can be a very powerful tool in the design and development of the next generation IC engines.
Such an approach can eliminate the current ad hoc processes used in the design of engines and may prove very useful when targets are difficult to meet or design changes are required.
With the advent of high performance computers and more efficient numerical algorithms, computational fluids dynamics (CFD) lends itself well for the support of such an optimization process. The concept of coupling CFD analysis with numerical optimization procedures was introduced in the early seventies by Hicks and Vanderplaats [5] who developed a numerical optimization technique for the design of lowspeed airfoils. More recently, [6], renewed interest in the field has emerged but mostly dealing with aerospace type design problems where inviscid flow equations are considered [7,8].
To date, only a limited number of optimization studies, dealing with the full viscous and turbulent Navier-Stokes equations [9,10], has been reported.
Since most fluid flow problems in automotive powertrain applications are viscous, the optimization of these systems requires the solution of the full Navier-Stokes equations to evaluate the performance. Choi and Kim [11] faced this challenge by applying the CFD based numerical optimization to the design of a low-noise automotive cooling fan.
In this paper, we apply a similar approach to another class of problems in automotive applications, namely internal flows in engine intake systems. In particular, we demonstrate the feasibility and describe the procedure involved in coupling CFD automated analysis with numerical optimization techniques to optimize the shape of an intake port for a typical 2-valve IC engine.
There are three key challenges for CFD based design optimization. First, for the methodology to be possible, the entire process of building the model, performing CFD computations and post-processing to extract the relevant metrics must be fully automated. In addition, the automated grid generation schemes have to be very robust and must preserve a high mesh quality for each new design geometry. This is indeed what we have developed through a combination of internal and commercial pre-processing codes. The description of the automated meshing techniques is beyond the scope of this paper and will be described in a separate publication.
The second challenge is the efficiency of the optimization algorithm. The computational expense of real-world flow simulations demands that the optimal design be computed with relatively few CFD solutions. To accomplish this, various optimization algorithms [12] may be used to perform the engine airflow optimizations. The most common are gradient-based numerical optimization methods. First-order gradient-based algorithms typically require far fewer function evaluations than zeroth order methods, i.e. genetic algorithms [13], and thus are well suited for CFD-based design optimization. Indeed, Choi and Kim [14] used the conjugate gradient method [12] in their cooling fan optimization.
The third challenge is the quick and accurate evaluation of the design sensitivities that are required by the gradient-based search algorithms [15]. Significant advantages may be gained through the use of analytical design sensitivity analysis; i. e. Wang et al. [16] used a direct differentiation method to avoid costly finite difference perturbation analyses during the optimization, while Reuther and Jamison [17] developed an adjoint variable method for the evaluation of the design sensitivities for their potential flow solver. The development of such techniques, however, is very tedious and must be performed for each performance measure and each type of analysis. To the knowledge of the authors, such formulation is not yet available for full Navier-Stokes type analysis. For our optimizations, the finite difference method is used for the evaluation of the design sensitivities. While this method is computationally more expensive, it allows more flexibility in dealing with various performance measures and design parameters.
In this paper, we first present our CFD-based methodology for port and chamber optimization, followed by a discussion on the objective function formulation and geometric design parameterization. A brief discussion of the CFD solution method used in our optimization studies is then given. Finally, results for constrained and unconstrained optimization problems to maximize the discharge coefficient while maintaining nominal levels of total angular momentum flux levels are presented to illustrate the methodology.

CFD BASED SHAPE OPTIMIZATION
The methodology employed in our CFD-based Port and chamber shape optimization is shown in Figure 1. Starting from a given set of design variables, in-house codes (Pro-Duct and Chamgen) are used to generate the 3D surface description of the port and chamber. A third in-house code (FloMod) uses both the parametric and the surface data generated to automatically grid the entire geometry and to set up the solver parameters and boundary conditions. A commercial CFD solver (Star-CD) is then used to conduct the analysis and to return the results back to FloMod for the   Yes No X l+1 X l evaluation of the objective and constrain functions, in this study discharge coefficient, Cd, and total angular momentum flux, Sw, respectively.
While the methodology is quite general, we choose to modify the geometry of the engine intake port only in this study. The performance measures are based on steady state type flow analysis at a fixed valve lift, i.e. discharge coefficient and total angular momentum flux. The problem as well as the metrics is described in more details in the following sections.

Optimization Problem: Formulation and Solution
The most general single objective optimization problem is one which minimizes an objective function F defined over the n design parameters X i , i = 1, 2, ... n while satisfying both equality and inequality constraints. In mathematical terms we write [18]: find to minimize or maximize , subject to: Constraint functions g and h are defined to mark the boundaries between which designs are allowed and which are infeasible. These constraints are divided into two groups, m inequality constraints g j , j = 1, 2, ..., m, and l equality constraints h k , k = 1, 2, ... The design parameters are assembled in a vector X and initial values are chosen for each component. The design space of all allowable designs is defined through side constraints that bound the upper and lower limits for each design variable X i as X i U and X i L , respectively.
For this study, as described above, gradient based algorithms are used since they require far fewer function evaluations than zeroth order methods. The unconstrained and constrained minimizations in this study are solved using the Broydon-Fletcher-Goldfarb-Shanno (BFGS) and the modified method of feasible directions algorithms, respectively (DOT) [19].

Design Sensitivity Analysis
Descent algorithms use design sensitivity information to compute the search direction and thus the optimal design. These design sensitivities (or design derivatives) quantify the relationships between the design variables and the performance measures. In this study, the design sensitivity of a performance measure F, with respect to a design variable X i is evaluated by the forward finite difference approximation as: where is a zero vector with the exception of the i-th location that contains . The main advantage of using the finite difference method is its ease of implementation. However, disadvantages of the finite difference method must be addressed when used with CFD analysis tools, namely, computational expense and degree of accuracy. Performance measures and thus the system response must be evaluated n + 1 times to compute the design sensitivity with respect to each of the n design variables. The computational burden of CFD-type analyses, therefore, present a limitation on the number of design variables that may be allowed to vary independently.
Additionally, when F is nonlinear in the design variable X i , the finite difference approximation depends on the perturbation size . When is too large, discretization errors occur and when is too small, round-off error due to limits in machine precision and numerical noise in the CFD solution corrupts the sensitivity calculation. Small but non-zero sensitivities are particularly susceptible to round-off error. Several numerical experiments were conducted to choose the proper settings for the size of the perturbation of each design variable.

Objective and Constraint Definition and Evaluation
It is well recognized that WOT (wide open throtle) performance measures, such as torque and power, are mainly controlled by the engine breathing capability. In contrast, at idle and part load conditions, where the engine characteristics of interests are idle stability, fuel consumption, and emissions, the performance is mainly controlled by the combustion rates. In turn, the combustion is controlled by the engine ability to generate enough turbulence at the time of spark and by its capability to mix the fuel and air [20,21]. A key mechanism to create and sustain a strong turbulent flow field at time of spark is the creation of large-scale in-cylinder motion, such as tumble and swirl, during the intake process. Such large scale flow structures act as storage devices for kinetic energy that is converted during the compression stroke to generate and sustain a higher levels of turbulence. Thus, the optimization of an intake port should allow for easy breathing of the engine and should generate a high level of in-cylinder motion. The ability of an engine to breathe is measured by the pressure losses incurred when flow travels through the intake port and chamber. This pressure drop is traditionally measured in steady state flow rig and the metric used is discharge coefficient, Cd. The in-cylinder motion has also traditionally been measured in steady state experimental facilities. In such experiments, the angular momentum fluxes are measured and then used to reconstruct the transient flow conditions [22,23]. Since the current optimization is based on steady state analyses, we define similar variables, namely discharge coefficient, Cd, and total angular momentum flux, Sw, to evaluate the performance of each design. The latter variable is an extension of the commonly used angular momentum flux coefficient for swirl and tumble. It defines the magnitude of the angular momentum flux vector, including all three components; namely the swirling, the tumbling and the cross tumbling. The three components of the momentum flux vector are evaluated by computing the moment of momentum that the inducted flow can generate around a specific reference point. The computational cells used for this analysis are located at the intake valve curtain area, while the reference point is at the center of the chamber at a distance of a half-stroke from the engine deck interface. Discharge coefficient, Cd, in this case is defined in the usual fashion, as ratio of actual mass flow rate to the isentropic mass flow produced by a similar pressure drop across the port and chamber.
As will be discussed in the result section, the discharge coefficient is used as the design objective function, while the total angular momentum flux defines the constraint function. The exception is only in case 2 described below, where the port shape is optimized solely to maximize Sw.

Design Parameterization
The intake port geometry significantly influences both the discharge coefficient and charge motion in the combustion chamber. In this study, we alter the shape of the intake port to generate better flow conditions and thus improved engine performance. In numerical shape optimization, there are two common methods used to generate new simulation models from design data: -geometry-based mesh parameterization [24] -design basis vectors [18]. The design basis approach is often employed when a geometry parameterization is inconvenient or not possible. For our problem, the geometry-based method is used. In this approach, all grid locations are related to higher level geometry data such as centerline radii or cross-sectional diameters so that the optimization output is generated in a form that is easily interpreted by design engineers. The parametric approach also allows the integration of casting rules and limitations, thus guaranteeing the manufacturing feasibility of the final optimum port shape. In our designs, the core split line with proper draft angles and uniform wall thickness are guaranteed by the parameterization.
Shape design parameterizations are often complicated since the relationship between each grid point and the shape design variable must be defined prior to the optimization. Each of the grid locations must also be specified prior to each function evaluation. In addition, it is very important, particularly for these CFD type applications, that mesh quality be preserved over the entire design space. Changes in the design variables, for example, should not distort the mesh such that the numerical solution diverges or reaches unacceptable numerical results. In our case the parameterization is carefully planned so that it satisfies the above criteria as well as it fully spans the intended design space with a small number of variables. As mentioned above, the casting rules imposed on the design, present a drastic reduction in the number of independent variables, and thus make this type of optimization practical.
In this work, only four independent design variables are used to parametrically represent the intake port surface. The design variables are geometric parameters that define a parametric surface geometry of the intake port. Here, we construct the shape of the port by sweeping between control cross sections along a path line. Each control cross section has sixteen parameters defining its perimeter. With the casting rules applied, the sixteen section control variables are not all independent. The path is defined by fitting a Bspline curve through three control points, which have prescribed first and second derivatives. It should be mentioned that the plane of each control section is always kept normal to the path.
The four design variables we chose for this study are defined as follows: X1, as shown in Figure 2a, is a path control variable that defines the inner radius of the arc defining a part of the path that originates at the interface with the valve seat. X2, as shown in Figure 2a, is a second path control variable that gives the angle of the arc with radius X1. X3, as shown in Figure 2b, is a section control variable that defines the distance along the topside of section perimeter measured before filleting the section corners. X4, as shown in Figure 2b, is a section control variable that defines the distance along the bottom side of section perimeter measured before filleting the section corners. The lower and upper limits used in the optimization for each design variable described above are shown in Table 1.
As shown in Figure 2a, X1 and X2 are used to control the shape of the path along which the port sections are swept. Here, the path originates at the base of the valve seat in a direction parallel to the valve axis. It starts with an arc defined by X1 and X2. The remaining part of the path is construct by a B-spline curve connecting the end point of the arc and the path point located at the port interface with the intake manifold.  The final two design variables are used to control the port cross sections. A control section, as shown in Figure 2b, is parametrized by X3 and X4 which denote respectively the lateral dimension of the top and bottom corner points before filleting. The filleting is automatically chosen to the highest value possible between the corner control points and the central points defined at the split line of the port. The port geometry, in this case, is symmetric and the cross sectional shape defined by X3 and X4 is maintained constant for a third of the length of the port, before it is allowed to smoothly blend into a circle at the base.
The symmetry constraint as well as the start of blending control point could have been allowed to vary at the same time as the four selected variables, but initial analysis indicated that they were not as sensitive. These sensitivity analyses will be the subject of another paper. Figure 3 presents the different shapes generated by the upper and lower bounds of the design variables described above. Figure 3a presents the extreme port path shapes considered, while Figure 3b presents the extreme area cross section shapes considered.

Flow Analysis
In this paper, the boundary conditions were a fixed velocity at the inlet and a fixed pressure condition at the outlet. The CFD analyses were conducted using the regular k-ε turbulent model, the SIMPLE algorithm, and a first order upwind difference scheme. These solver settings guaranteed existence of a fully converged solution for each set of parameters. The optimization technique required high CFD accuracy for the evaluation of objective and constrain functions to compute the sensitivity for each variable. In order to keep the analysis at a manageable cost; the CFD mesh (Fig. 4) was kept relatively coarse, namely 37 000 cells defining the flow domain. In addition, to further improve the robustness of the grid generation tools, the CFD model was simplified by excluding the valve guide and stem in the port. An entry section with an area fifty times that of the port opening at the base of the valve seat was added to the CFD model to provide smooth inlet conditions. The shape of this inlet expansion section was set to resemble a bell-mouth shape in order to eliminate any additional entry losses.
In order to reduce total time of computation, the analysis were initiated from previous design iteration results whenever available. This capability was possible since the same number of grid cell distributed on exactly the same topology was used in every function evaluation.

RESULTS AND DISCUSSION
This study involves both constrained and unconstrained design optimization problems. First unconstrained methods are used to maximize the discharge coefficient Cd and total angular momentum flux Sw, defined here as cases 1 and 2 respectively. These represent ideal situations where an uncompromised design can be reached. They, therefore, do not particularly represent the real engine design problem where a compromise between the engine WOT performance, idle stability, dilution tolerance, etc must be reached. These results are included in this paper as an indication of the validity of the subsequent analysis and as limiting cases for the problem of interest.
The actual optimization problem of interest is addressed in the subsequent constrained optimizations, cases 3, 4, and 5. These design optimizations are conducted to maximize discharge coefficient while maintaining preset minimum values of total angular momentum flux.

Unconstrained Port Shape Optimization: Maximizing Cd and Sw Independently
In this section, the results of the unconstrained optimization problem for the two first cases are presented. In both cases, the optimization started from the initial design shown in Table 2. This initial shape was chosen to represent a port design that is within the design domain conforming to well established design guidelines; i.e. area ratios, maximum bends, etc. The same design is also shown in graphical form in Figure 5a. In this case, the end section shape is perfectly circular having the same area as the valve throat at the base of the seat insert. The path is a smooth, monatomic B-spline function generated by a large radius of 40 mm extended over a angle of 18 degrees. The optimization for these two cases used the BFGS method [19] and converged in three design iterations. Each design iteration required approximately seven CFD function evaluations. Also shown in Table 2 are the results of both unconstrained optimizations. In the first case, where only Cd is maximized, the optimization converged to the design that provides the largest possible inlet area (largest possible values for both X3 and X4), as shown by the shape in Figure 5b. The path in this case converged to a shape having the smallest values of radius Figure 5 Initial and final port shapes.
X1 and angle X2. A close inspection of Figure 3a indicates that in contrast to all other X1 and X2 extreme combinations, the shape of the path in this case does not show any inflection point. Such a shape represents the ideal conditions for eliminating separation or stagnation in the port, reducing skin friction at the walls, and minimizing expansion losses at the inlet of the valve. The velocity and pressure fields for this design, illustrated in Figure 6 confirm this conclusion. Both the velocity vector field (Fig. 6a) and the pressure field (Fig. 6b) present a gradual decrease in their respective magnitudes and show no indication of separation. In particular, the pressure field in the port appears to decrease linearly along the path of the port in a similar fashion to fully developed pipe flow, thus indicating a lack of local losses in the port. In addition, due to the large cross section of the port  itself, the velocities in the port are relatively low for the set pressure drop across the engine (only reaching a maximum value of about 60 m/s). These lower velocities result in a reduction of the friction losses in the port. It is worth noting that the low velocities in this case also contribute to the reduced values of the total angular momentum flux, Sw, shown in Table 2.
In contrast, in the case where only the total angular momentum flux is maximized, the optimized shape of the port converged to the one where the cross sectional area is reduced to its lowest possible value, i.e. minimum values for both X3 and X4. The path, in this case, also reaches limiting values for both of its controlling variables, namely the lower limit for the radius of curvature X1 and the upper limit for the angle of curvature X2. A close inspection of this geometry (Fig. 5c) indicates that this geometry represents the best condition for targeting the flow to the far side of the valve. Therefore, this geometry biases the inlet of the flow in the chamber thus creating potential for high level of in-cylinder large-scale motion, in this case potential to develop tumble. The large increase in the value of total angular flux (Sw) also shown in Table 2, is the results of two phenomena; i.e. the increase in the velocities in the port due to the narrowing of the port cross section and the severe biasing of the flow caused by the Figure 7a, shows the resulting flow field for this case. The flow is highly biased to the far side of the intake valve and the velocity magnitudes are much higher than the earlier case (Fig. 6a), reaching levels of up to 95 m/s in the port. Note that this increase in total angular momentum flux was obtained at the expense of the port breathing capabilities. As can be seen in Table 2, a 20% reduction in flow capability was incurred with this design. The reasons for this reduction in discharge coefficient (or increase in pressure drop across the port and chamber) is apparent from the pressure field presented in Figure 7b. Velocity and pressure fields (maximizing Sw) sharp lowering of the port path close to the valve seat.
In contrast with case 1, the pressure in the port reaches a maximum value close to the valve seat. This increased pressure is due to the slowing down of the flow near the valve imposed by the sharp bend in the geometry. Actually, the bend is so severe in this case that the flow almost stagnate at the upper part of the port, Figure 7a.
Not surprisingly, the optimum design for an intake port capable to generate a high level of in-cylinder motion is achieved at the expense of the breathing capability of the port and vice versa for an optimum high flow port. Therefore, in practical applications a compromise always exists because of the different requirements on the engine flow characteristics to satisfy the combustion needs at the each operating conditions, i.e. high Cd for increased torque and power at WOT and high in-cylinder motion for enhanced burn rates at idle and part load conditions. The higher in-cylinder motion is in turn responsible for enhanced burn rates, which are more favorable to more stable idle conditions and higher dilution limits for emission control and fuel economy. The design decisions in this case should be based on a trade-off analysis, which is the subject of the next section.

Constrained Problem: Maximizing Cd while Maintaining Nominal Levels of Sw
Here, we show results of the constrained optimization problem, which deals with maximizing the discharge coefficient, Cd, while maintaining prescribed levels of total angular momentum flux, Sw. Three different cases are investigated dealing with different levels of Sw, namely minimum Sw values of 0.4, 0.35, and 0.3 respectively. These inequality constraint values are chosen in order to establish the trade-off relationship that exists between Cd and Sw. The optimization for these three cases used the modified method of feasible directions [19] and converged in five design iterations, which required more than fifty CFD function evaluations Figure 8, shows the optimization history for the case where the inequality of Sw > 0.3 was applied. Note that Cd and Sw converge in opposite directions. The discharge coefficient increases from 0.501 to 0.626 while Sw decreases from 0.42 to 0.3. Note also that the angular momentum flux constraint is active at the optimum design. Again, this suggests an inverse relationship between the engine capability to breathe and to generate strong in-cylinder motion. This point will be further clarified when the trade-off curve is presented later in this section.
The results for the three constrained problems are summarized in Table 3. Note that, as expected, the optimum designs requiring higher Sw are associated with larger values of design variable X3 and X4 which determine the size of the port inlet area. This is clearly seen in Figure 9 where the optimum shapes of the three different cases are illustrated. Note that the shape of the path appears to be moving monotonically from that of a straight directed port in the first case, Figure 9a, to the one of a low floor type port, Figure 9c.  Optimization history of Cd.

Figure 9
Optimum design shapes (maximize Cd with constraints). Table 3 indicates that the optimum value for Cd varies inversely with Sw confirming the existence of a trade-off between Cd and Sw. This trade-off is illustrated in Figure 10 where all the optimization results from the constrained and unconstrained analyses are combined. Actually for the purpose of such a plot, the unconstrained analyses are regarded as the limiting cases for the constrained cases. The trade-off curve appears to be quite flat up to an Sw values of about 0.34. Beyond this value, the tradeoff curve becomes very steep. This dual behavior suggests that the flow is controlled by two different regimes at either asymptotic part the trade-off curve. Below this value of 0.34 for Sw, the flow capability of the intake port design is relatively independent of the level of total angular momentum required. In this case, it is expected that the flow does not separate significantly in the port and that the pressure drop is mainly controlled by the skin friction in the port and the final expansion loss in the chamber. However, higher angular momentum flux levels are obtained at significant loss to the engine flow capability. For Sw > 0.34, it is expected that the pressure drop becomes dominated by the reduced effective area of the port resulting from flow separations. Trade-off curve. Note that the point separating the two regimes is not believed to be universal but primarily dictated by the combination of valve angle and intake manifold flange location and angle.
From a design point of view, the trade-off curve offers an idealized look at all possible optimum shapes. It may be used to identify the benefits and limitations of an engine and manifold layout. Actually, since each point on curve represents an optimum design, one can decide on the ideal design based on the requirements of the engine or vehicle line at hand; i.e. the optimum design of engines intended for a heavy truck application, where power and low end torque is important, should be selected from the left side of the curve (low constraint values on strength of cylinder flow motion), while the optimum design of an engine intended to be used in a luxury vehicle line, where combustion refinement is required, should be selected from the right hand side of the curve. In the latter case, the required level of in-cylinder motion can be predicted through the use of engine simulation analysis [20,21].

CONCLUSION
A methodology for design of ports and chambers of IC engines with improved performance, efficiency and emissions has been presented. This methodology combines three-dimensional, steady state CFD techniques with robust numerical optimization tools to design, rather than just evaluate the performance, of IC engine ports and chambers. Significant challenges in automating the CFD mesh generation, data flow control, and post processing of the resulting flow fields are addressed in order to make the methodology possible. In addition, mathematically based design parameterization of the ports and chamber are implemented to reduce the number of design variables and thus to enable the optimization procedures with reasonable computational cost.
Various optimization problems are solved to illustrate the design methodology for an intake port of a typical two-valve IC engine. Both constrained and unconstrained optimization techniques are employed to generate a trade-off relationship between the engine flow resistance and its capability to generate large-scale in-cylinder motion. These results may be used by engine cylinder head engineers to select the optimized geometry that best addresses their particular design needs. The trade-off curve clearly identifies the envelope of best possible designs, providing the design engineers with a significant advantage in today's engine development environment.