A Numerical Study of Response Surface Based Shape Optimization Using Neural Networks

Industrial pipes that are used for fluid transport generally have to undergo many changes of shape to accommodate interfacing equipment related to plant operation, which results in flow maldistribution zones and higher pressure drops, and in turn leads to higher power consumption. In an attempt to redress this problem, ANSYS, a commercial Computational Fluid Dynamics (CFD) software, is used to perform numerical simulations based on a deterministic computational model of the internal fluid flow using the Reynolds Averaged Navier Stokes equations (RANS), a multi objective optimization study employing Response surface methodology and artificial neural networks. This numerical analysis has been performed on a galvanized steel duct for water recirculation. The focus of the paper is the study of the effect of a chosen set of several geometrical dimensions on the pressure drop and flow distribution inside the duct. Subsequently, a new set of designs with different geometrical parameters has been obtained to minimize the pressure drop and achieve a more uniform flow distribution by using artificial neural networks to generate a response surface and further employing Screening (Shifted-Hammersley sampling) as the optimization method that was used to select the best designs from amongst those that have been generated from the response surface.


Introduction
The shape of a duct in fluid transport systems is not often designed from a fluid dynamic perspective. This is mostly due to the interacting equipment being given priority over the ducts, resulting in a constraint of restricted size. The constraints are often imposed when developing a system, as the ducts are fitted in between and around other components. As a result, the duct shape can often be very complex, leading to exit flow maldistribution and higher pressure drops, resulting in higher consumption of pumping power. However, in recent years, optimization based on flow analysis is becoming increasingly popular in the field of engineering design but at the same time, manufacturers are also increasingly looking for ways in which to reduce the time taken during product development. One such way is having the products developed in the form of virtual prototypes in which CAD software is used to model a prototype which can be used to predict the performance prior to constructing physical prototypes.
Manufacturers can then explore the performance of a myriad of design alternatives without having to invest the money and time required to build physical prototypes [1,2]. Response surfaces methodology, as such that are employed in ANSYS, offers an efficient way to get the variation of a given performance with respect to input parameters and provide a continuous variation of the performance over a given variation of the input as evidenced by Winne et al and Behin et al [3][4][5]. Firstly, an acceptable range of variation for each input parameter (thus defining the design space) is given, either manually or automatically. A response surface is then generated from the Design of Experiments (DOE). A few points within the design space are solved in order to provide the approximated values of the output parameters everywhere else in the analyzed design space, without having to perform a complete solution. A response surface is thus an approximated best fit surface computed from the DOE results for each output parameter in terms of the input parameter and this approximation can be carried out using artificial neural networks as seen in the works of Wolfram et al and Sun et al [6,7].
The idea of artificial neural networks is based on the known architecture of the human brain. The human brain is made up of neurons, which are information-processing units operating in parallel and connected by synapses, which relay the information between neurons. It is capable of processing input signals (e.g. taste) that are picked up by various sensors (e.g. tongue) from the environment and of providing appropriate output signals (e.g. determining whether it is sweet or sour). The advantages of the human information processing paradigm are complexity, nonlinearity, and parallelism and so an artificial neural network resembles the human brain in many respects. It is similar in that it is a multi-layered structure constituted by neurons which are connected by synapses, has the ability of mapping input signals onto output signals and to learn certain tasks during a training phase by adjusting the values of connections (weights) between the neurons.
The output produced by a neural network in the current case is the response surface. The idea of such a response surface is then to replace the deterministic computational model for CFD analysis with a neural network [8]. Consequently, the network first needs to learn the features of the underlying deterministic computational model. That is, the input signals comprise the geometrical parameters and the network output provides the associated response surface in the form of pressure drop and flow distribution. The principal aim of this paper, to minimize the pressure drop and achieve a more uniform flow distribution within the duct, was thus achieved by drawing the best possible designs from the response surface given by the neural network by using the screening optimization method.

CFD Approach
ANSYS Fluent was used to perform the simulations. The materials used were, steel with a thick-ness of 2mm and water at 80 degrees Celsius with an inlet speed of 0.5m/s. The Navier-Stokes equations are solved using the SIMPLE scheme and the second order upwind discretization methods of the finite volume method. The realizable k-model was applied to model the flow turbulence and standard wall functions are used to manage the near wall turbulence layer.

Governing Equations
To account for turbulence, each variable is decomposed into an average component and a fluctuating component (for example the speed, pressure, temperature) [9,10]. If we replace the instantaneous values of the transport equation by their expression in the conservation equation, it yields the Reynolds stress and equations called the RANS (Reynolds Averaged Navier-Stokes). The governing equations for the CFD models are presented below: Where (1) is the continuity equation Where is the mass density and u is the flow velocity. The equation of motion describing an incompressible Newtonian fluid is shown in the equation (2).
Here, p represents the pressure, g is the gravity, and is the viscous stress tensor. Equation (3) is the general differential energy equation.
Where, E is the total energy, k is the fluid thermal conductivity, T is the absolute temperature and Φ is a function representing the dissipation of energy due to viscous effects.

Grid Independence
The dimensions and boundary layer profiles vary widely depending on the properties of the fluid used and the type of flow. In order to have a broader approach, and one that is independent of these values, we introduce the concept of dimensionless distance. The distance y + is a normal dimensionless distance from the wall and is obtained by the following formula: With * the shear rate, + shear stress at the wall, y is the distance to the wall and 1 the kinematic viscosity of the fluid. The value of y + can be calculated anywhere in the area, though it has a different meaning when calculated in the CFD software. In this case, the value of y + actually corresponds to the value of y + at the first node of the mesh. The practical value of this value is to define the different areas of the boundary layer as well as verify that the fineness of the mesh is a good fit with the chosen turbulence model.

Goal Driven Optimization
There are various optimization tools available in Ansys such as the adjoint method or Goal driven optimization (GDO), which will be used in this paper, which is a set of multi-objective techniques in which the "best" possible designs are obtained from a sample, given the goals that have been chosen for the parameters [11]. There are two types of GDO systems, namely: Direct optimization, wherein the "best" possible designs are selected from actual analysis results which have been solved and response surface optimization wherein the "best" possible designs are drawn from the generated response surface which provides a continuous variation of a chosen output/(s) with respect a given range of the inputs without having to solve for the entire design space. There are different techniques for generating response surfaces such as the kriging meta-model but this paper will be focusing on using neural networks [12].

Artificial Neural Networks
Artificial neural networks are a type of supervised learning system. The system is trained to perform a particular task by adjusting the connections (weights) between the neurons of the network. Below, in Figure 2 is a typical neural network with an input layer, an output layer, and 3 hidden layers. Each synapse has a value, w, which can be adjusted.
As shown in Figure 2, the values of each neuron within a specific layer, are dependent on the values of neurons in the preceding layer and the value of the weights connecting them. Each neuron has a function, g, which collates the incoming signals into a summation function, with each input signal weighted by a synaptic weight, w. The result is then squashed into a range of 0 to 1 by multiplying it with an activation function, 2, and thus producing the value of the next neuron.

The Back Propagation Learning Rule
A learning rule, or training algorithm, is an iterative technique for adjusting the weights of a network in order to train the network to perform some task [13,14]. This is done by first introducing a cost function, G H 68 9,3 ; >, based on the root mean square error, e rms . It is a function of the known values (sample of solved design points, d m ) and the output of the neural network (the response surface), y m . This is a function that is used to evaluate the performance of the neural network at each iteration. G H 68 9,3 ; > = IJ K;L M (7) The goal is to minimize this function by using the method of gradient descent and the algorithm for computing this function is called back propagation. The negative gradient of the cost function is used to adjust the weights of the neural network thus minimizing the error whilst checking the response surface at each iteration against the solved design points. The gradient descent algorithm for the root mean square error is written as follows, for the i th iteration, Where O is the learning rate, denoting the magnitude of the step in the direction of steepest descent:

Optimization Method
There are various optimization algorithms that can be used, such as the Genetic algorithms [15,16]. For this paper however, Screening (Shifted-Hammersley sampling) was used. It is a goal-driven optimization method that draws a number of design points (samples) from the response surface (population), with the goal of finding the "best" designs within the response surface, without having to sample the entire design space. It is a non-iterative direct sampling method by a quasi-random number generator which has a low discrepancy. Quasi-random numbers act like random numbers but unlike true random numbers, are uniformly distributed in some domain, generally the unit hypercube [17]. Whereas true random numbers have no correlation whatsoever, quasi-random numbers are quite correlated and chosen in such a way that new numbers do also uniformly fill the gaps between the old ones. A quantitative measure for this property is the discrepancy of the point set, which should be as small as possible. The discrepancy is therefore computed by comparing the actual number of sample points in a given volume of multidimensional space with the number of sample points that should be there assuming a uniform distribution. Hammersley Sequence Sampling (HSS) uses Hammersley points to uniformly sample a unit hypercube and inverts these points over the joint cumulative probability distribution (0, 1) to provide a sample set for the variables of interest. The design of Hammersley points is given below. Any integer n can be represented as a sequence of digits ) 0 , ) = , ) C , ) U , V , ) ; by the following equation: The integer n can then be written in radix-R notation (R is an integer) as follows: A unique fraction between 0 and 1 called the inverse radix number can be constructed by reversing the order of the digits of n around the decimal point as follows: X Y ) ) 0 W <= ) = " W <C V ) ; " W < ;<= (12) Thus, for a k-dimensional search space, the Hammersley points are given by the following expression: Where W = , W C , V , W 3<= are the first k-1 prime numbers and i=0,…, N indicates the sample points. This algorithm generates a set of N points in the n dimensional design space I0, 1M N . Now, from the plot of these points, it is seen that the first row (corresponding to the first sample point) of the Hammersley matrix is zero and the last row is not 1. This implies that, for the k-dimensional hypercube, the Hammersley sampler generates a block of points that are skewed more toward the origin of the cube and away from the far edges and faces. To compensate for this bias, a pointshifting process is proposed that shifts all Hammersley points by the following amount: ∆ = C`.

Near-Wall Mesh Quality
Standard wall functions, which are the most widely used in the industrial flows, were applied for capturing the near-wall flow features. It is recommended in ANSYS-Inc. (2009) that for using the standard wall function, the value of y+ in the first cell should be within the logarithmic law layer (30<y+<300). During the meshing of the duct, an inflation method was employed with a first layer thickness of 0,0025 m and thus maintaining a constant first cell height. Figure 3 shows the value of y+ near the wall of the duct for different meshes with an increasing number of elements. The study of y+ shows that successful computation of turbulent flow close to the walls of the duct can be achieved by using any one of the different mesh qualities.

Mesh Independence Test
A grid convergence study was performed, referencing the work of Karimi et al by running simulations on the eight different meshes from the y + study to predict the drag coefficient and determine how the mesh quality affects the results [18].    Table 1 displays the range of the inputs used to generate the response surface. The flow uniformity is accounted for by the standard deviation of the outlet velocity. Standard deviation is a statistic that measures the dispersion of a dataset (outlet velocity flow field) to its mean. This means that uniform flow corresponds to a minimized standard deviation of velocity. 24 design points were generated using the Custom design of experiment type in conjunction with the central composite design. The 24 design points then served as the source of the training data used by the neural network to generate a response surface. Figure 5 shows a Predicted versus Observed chart, which displays the goodness of fit data for both the output parameters. It is used to check how accurately the response surface generated by the neural network, can predict the design points from the DOE. Initially, only 18 design points had been generated by central composite design type but the goodness of fit metric had not been satisfactory. The custom design type was then employed to add 6 more design points to the training data to improve the goodness of fit of the response surface. Additional design points, which have been solved, were also added to further verify the accuracy of the response surface and these are known as verification points. Some of these verification points were then further used as refinement points and added to the training data for improving the accuracy of the response surface. As observed from Figure 5, most of the points fall either on or near the line which means that the response surface generated by the neural network was able to sufficiently predict most of the values of the design points within the given range, including the verification points.

Output Sensitivities
Local sensitivity data is also available to outline how the output parameters change with respect to the individual input parameters, at a specific response point, by varying a single input parameter whilst keeping the other input parameters constant. The local sensitivity bar plot in Figure 6 highlights how the Incident angle (P3) is inversely correlated to both the pressure drop (P6) and the standard deviation (P7) of the output velocity whilst also having the highest sensitivity.
Conversely, the outlet length (P1), middle length (P4) and Inlet length (P5) have a positive correlation with the pressure drop but have a much smaller impact as compared to the incident angle. Relative to the standard deviation, the outlet length (P1) has a positive correlation and stronger impact compared to the middle length (P2) and Inlet length which both have a negative correlation and weaker impact on the standard deviation of the outlet velocity.

Response Surfaces
With a three-dimensional response surface, the combined effect of two different input parameters on an output parameter can be visualized. Figure 7, highlights how the outlet length (P1) and the Incident angle (P3) affect the pressure drop (P6). As observed from the local sensitivities, the Incident angle (P3) has an overwhelmingly stronger influence on the pressure drop (P6) than the outlet length (P1). The variation of the outlet length (P1) along its axis, has very little effect on the pressure drop (P6), on the vertical axis. Unlike the outlet length (P1) though, and as also observed in the previous section on local sensitivities, variations in the Incident angle (P3), greatly affect the pressure drop (P6).

Optimization
The screening optimization method, also known as Shifted-Hammersley sampling, is then used to sample the response surface in search of the optimum design dependent on the set objectives. An optimized design, in this case, would have a minimized pressure drop (P6) and also a minimized standard deviation of the output velocity (P7). This is, therefore, a case of multi-objective optimization where the aim is to simultaneously treat a number of conflicting objectives. A trade-off surface can thus be used to view how the different output parameters conflict.

Tradeoff Chart
A Tradeoff chart is a scatter chart representing the generated samples with the output parameters set on different axes. This allows for the visualizing of how an output goal can be achieved in relation to other output goals. In essence, it highlights whether achieving one goal might be detrimental to achieving the other objectives. The groups of samples generated by the optimization algorithm are then colored according to the Pareto front they belong to with blue representing the best set of samples (first Pareto front) and red representing the worst (worst Pareto front) as displayed in Figure 9.

Screening
After running the optimization algorithm, it produces the 3 best candidates as per the set objectives. Table 2 shows the input parameters of the 3 candidates drawn from the response surface along with the predicted outputs. The "real solve" of the candidate points is also shown along-side the predicted outputs.
Candidate point 1 had the lowest pressure drop and the lowest standard deviation of the outlet velocity and hence was the best candidate point produced by the algorithm. It is also apparent how the input parameters of the produced candidate points are as expected based on the local sensitivity data and the upper and lower bounds of the response surface. A maximized Incident Angle (P3), for example, was commensurate with the set objectives and the optimization algorithm produced the highest possible value of 120 degrees which was observed on the generated response surface in Figure 8. The available global sensitivity data also backs up the fidelity of the results. Global sensitivities are based on a correlation analysis using the generated sample points, which are drawn from all the possible values for the input parameters on the design space.     Table 3 shows a side by side comparison of input and output parameters of the original design against the optimized design. The optimized design managed to achieve a 7,89% decrease in the pressure drop and a 24,27% improvement in the flow uniformity.

Comparison Between the Original Design and the Optimized Design
The attained improvements are visualized in Figure 11, with the top images showing the flow uniformity on the outlet velocity profile and the second images showing the physical differences of the designs.

Conclusions
This study has demonstrated how competent ANSYS is at solving shape optimization for a fluid domain using parametric models to successfully decrease the pressure drop and simultaneously make the flow more uniform. However in future works, experimental data, collected using apparatus such as a hydraulic bench and Particle Image Velocimetry, should be used in order to validate the results along with a comparative analysis against other optimization methods such as Bayesian optimization.