A stochastic topology optimization algorithm for improved fluid dynamics systems

Abstract

The use of topology optimization in the design of fluid dynamics systems is still in its infancy. With the decreasing cost of additive manufacture, the application of topology optimization in the design of structural components has begun to increase. This paper provides a method for using topology optimization to reduce the power dissipation of fluid dynamics systems, with the novelty of it being the first application of stochastic mechanisms in the design of 3D fluid–solid geometrical interfaces. The optimization algorithm uses the continuous adjoint method for sensitivity analysis and is optimized against an objective function for fluid power dissipation. The paper details the methodology behind a vanilla gradient descent approach before introducing stochastic behavior through a minibatch-based system. Both algorithms are then applied to a novel case study for an internal combustion engine's piston cooling gallery before the performance of each algorithm's resulting geometry is analyzed and compared. The vanilla gradient descent algorithm achieves an 8.9% improvement in pressure loss through the case study, and this is surpassed by the stochastic descent algorithm which achieved a 9.9% improvement, however this improvement came with a large time cost. Both approaches produced similarly unintuitive geometry solutions to successfully improve the performance of the cooling gallery.

Introduction

The application of algorithmic methods in the design of fluid–solid components aims to realize more complex and efficient designs, for internal flow components, that minimize undesirable energy losses. The algorithm should solve a problem in which a fluid–solid geometrical interface is optimized while one or more governing equation for fluid flow is solved (Alexandersen and Andreasen, 2020).

Topology optimization is an optimization method that has been commonly applied to reduce the weight of structural components, without compromising mechanical strength (Bendsoe and Sigmund, 2004). It is a mathematical method that spatially optimizes the distribution of material within a defined domain, for a given set of boundary conditions and constraints (Rosinha et al., 2015). The method acts to reduce a predefined cost function, aiming to find the minima of this objective. Topology optimization has been shown to be beneficial for solving multiphysics problems in which different phenomena interact, due to its ability to realize complex and unintuitive solutions (Gersborg-Hansen, 2003). This ability makes it a competitive approach to employ in the design of fluid–solid components, as the best way to design these is often poorly understood due to difficulty in visualizing the sensitivities of different designs. The importance of tools that aid the design of fluid–solid interacting surfaces is easy to understand, as even the simplest incompressible case will spontaneously give rise to a large variety of complex coherent structures and phenomena (Pierrehumbert, 2022).

Additive manufacturing technology is now reaching a point in which it is a viable solution to produce durable complex metal geometries, with motorsport manufacturers utilizing this technology to produce highly stressed components such as additive manufactured exhausts, or even laser metal fused piston heads as done in 2020 by MAHLE GmbH for Porsche (MAHLE, 2020). The development of a robust and accurate design process for optimizing fluid dynamic systems could allow for even greater efficiencies to be reached by these additive manufactured components, with little or no cost to the timeline of projects. Engineers using such computational methods as a tool could even speed up the development timeline of products, as optimum design concepts can quickly be defined without relying on a product being passed between design and simulation iterations.

Topology optimization of fluids

The topology optimization approach was first introduced in 1988 for stiffness design of mechanical structures (Bendsøe and Kikuchi, 1988). It took 15 years for the approach to be applied to fluid dynamics problems, for a 2D Stokes flow-based minimization problem (Borrvall and Petersson, 2002). The errors incurred by Stokes flow were rectified with a 2D Navier–Stokes based lubrication approximation approach, however, this method would naturally face a road block for application to 3D problems (Gersborg-Hansen, 2003).

To optimize 3D fluid–solid interfaces, the fundamental mathematical methods driving research are the level-set method (Duan et al., 2008) and the adjoint method (Othmer et al., 2007). The level-set method is widely used in structural optimization as it gives autonomous and flexible handling of topological changes (Wang et al., 2003). The main drawback of the level-set method is that all sensitivities are located at the interface, and thus, all design changes must propagate from this interface, meaning that no internal structures or holes can form automatically (Alexandersen and Andreasen, 2020), this limitation can be overcome by the adjoint method.

The adjoint method was used to bridge the gap to more refined results (Othmer et al., 2007). The adjoint method works to compute the sensitivities of a defined cost function depending on the porosity value of each cell, with the sensitivities allowing the optimizer to rank which cells have the largest influence on the cost function (Othmer et al., 2007). This method allows the geometry changes to be described by a parametrized surface, which has a lower computational cost compared with volume mesh methods. This parametrized surface is effectively an immersed boundary produced due to the porosity variations from cell to cell.

The adjoint method that forms the fundament of the algorithm developed in this work is based on a fictitious porous media approach. In a fictitious porous media approach, the design domain is treated as a porous medium with each cell's porosity being varied to be solid or to be fluid-like to determine the geometry of the domain. Each cell in the domain gets assigned an individual porosity value, which is modeled by Darcy's law. Olesen et al. (2005) first developed this method to optimize 2D designs with coarse meshes, however, their method struggled with finding the true global minimum of these problems.

Application of the algorithm

A useful application of topology optimization is to optimize cooling passages, as increased flow efficiency will produce improved cooling due to a steeper temperature gradient. This problem has been approached by Kim and Son who applied it in the design of a U-bend turbine serpentine cooling passage (Kim and Son, 2019). Their optimization produced a reduction in pressure loss; however, the process took a significant period of time to run, 6 h, and produced results that required a large degree of post-processing to remove inappropriate features and smooth rough surfaces. The requirement of significant post-processing indicates a failure to find the global minima.

This paper will apply a topology optimization algorithm to the design of a piston cooling gallery. A piston cooling gallery is a good test for the algorithm as it represents the expected use of the method to improve fluid dynamics systems with one inlet and one outlet. The algorithm is applicable to any fluid channel in which there is only one fluid and there are inlet and outlet boundary conditions. This piston cooling gallery also introduces some more complex flow interactions through the zones of fluid impingement derived from the fluid parting around either side of the gallery, and areas of separation at the inlet and outlet. This will demonstrate the algorithms applicability to solve problems with complex fluid interactions.

By improving the cooling of a piston, an internal combustion engine is able to run at a higher thermal load without piston temperature being so high as to reduce engine durability (Deng et al., 2018), be that through expansion-generated ring scuffing or the onset of surface ignition. With piston cooling galleries being the primary cooling strategy for high loaded pistons (Thiel et al., 2007), a small increase in cooling efficiency could produce a sizable difference in total possible engine power output.

Optimization method

It was noted that a frequent issue with research into the topology optimization of fluid systems was the failure of the optimizer to find the global minima of the problem. This is due to the majority of these algorithms relying on a vanilla gradient descent-based optimization process (Alexandersen and Andreasen, 2020), which is widely known to be prone to convergence inefficiencies. Taking recommendation from the further developed field of structural topology optimization, gradient-free optimization methods are not a solution to this problem due to large inefficiencies for problems with many variables, such as topology optimization, however, gradient-based random-process methods may be beneficial (Sigmund, 2011).

The stochastic version of gradient descent (Robbins and Monro, 1951) has become a fundament for optimizing deep neural networks (Lu, 2022), and offers the potential to benefit the topology optimization of fluid systems by avoiding bad local minima or saddle points in the objective function (Kleinberg et al., 2018) at decreased computational cost (De et al., 2019). The application of stochastic mechanism into topology optimization simply refers to the introduction of randomness into the form-finding algorithm. Stochastic optimization has achieved some attention for structural optimization (De et al., 2019) but none for fluid–solid interface optimization, in which one can argue that it is more beneficial, as fluid-based topology optimization problems are likely to lead to functions with many sharp local minima due to the complexity of the interaction between the two states of matter. One way of achieving stochastic gradient descent is to recompute gradients for similar topology examples before each parameter update, avoiding vanilla gradient descent's tendency to converge to the nearest local minimum and instead jump to new and potentially better local minima (Ruder, 2017).

Proving the applicability of stochastic descent for fluid dynamic topology optimization could unlock a path to improved optimization performance of fluid–solid interaction surfaces. The performance of stochastic descent can be fairly evaluated by comparing it to a vanilla gradient descent optimization method for a case with the same boundary conditions. Aside from the optimization method, the performance difference will be demonstrated while keeping the underlying methodology between algorithms the same.

Methodology

The design optimization presented in this study will be developed using the continuous adjoint method, as this method is able to provide insightful design improvements with low computational cost (Othmer, 2008).

The general idea of an optimization problem is to minimize a cost function while subject to different constraints. Optimization of fluid dynamics problems is challenging as each geometry change produces continuously changing variables that interact with each other in complicated ways. The optimization algorithm explored in this work follows the process presented in Figure 1.

Fig. 1. Topology optimization algorithm process, from the initial domain to the final CAD output.

The adjoint method solves fluid dynamics optimization problems through using only two applications of finite volume numerics per iteration, to solve the system of partial differential equations for the 3D domain, called the primal and the adjoint equations. These determine which areas of the domain should be solid in order to minimize a cost function (J). Thus, it requires relatively low computational cost due to calling on only two solving operations.

The adjoint method is commonly expressed as a minimization problem involving the flow variables velocity (v) and pressure (p), and a design variable alpha (α), which are constrained by the incompressible steady state Navier–Stokes equation (Othmer et al., 2007), depicted in Eq. (1).

(1)$$\eqalign{& {\rm Minimize}\;J = J( {v, \;p, \;\alpha } ), \cr & {\rm Such}\,{\rm that}\;( v\nabla ) v + \nabla p-\nabla ( {2vD( v ) } ) + \alpha v = 0, \cr & \nabla v = 0,} $$

where v is the kinematic viscosity, the sum of molecular and turbulent viscosity, and D(v) is the strain rate tensor represented by $D( v ) = 0.5( {\nabla v + {( {\nabla v} ) }^T} ).$ The Navier–Stokes constraints are referred to as the primal equations and these are solved by a Semi-Implicit Method for Pressure Linked Equations (SIMPLE) algorithm (Patankar, 1980). The primal equations use a k-epsilon viscous model for turbulent modeling (Launder and Spalding, 1974).

During optimization, the design domain is treated as a porous medium, and the design variable α determines cells in the domain that should be fluid-like (low porosity values) or solid (high porosity values) depending on sensitivities. This porosity field controls the geometry and is added as a sink term α ⋅ v to the Navier–Stokes equations using Darcy's law (Darcy, 1856).

This is a constrained optimization problem, which can be tackled by introducing adjoint velocity (u) and adjoint pressure (q) as Lagrange multipliers to a Lagrangian relaxation function L (Othmer et al., 2007).

Minimize:

(2)$$L = J + \int_{\rm \Omega } {( {u, \;q} ) R{\rm d\Omega }}. $$

In Eq. (2), Ω represents the flow domain and R is the incompressible steady state Navier–Stokes equations and continuity equation, as in Eq. (1). To get the desired sensitivity for the Lagrange relaxation with respect to the design variables, the total variation of L is introduced, with the adjoint velocity and pressure being chosen in such a way that the variation of L in respect to velocity and pressure becomes negligible.

This is done through:

$$\delta L = \delta _\alpha L + \delta _vL + \delta _pL,$$

(3)$$\delta _vL + \delta _pL = 0.$$

These combine to give an expression for the sensitivity of the cell i as below, where V _i is the volume of cell i:

(4)$$\displaystyle{{\partial L\;} \over {\partial \alpha _i}} = u_iv_iV_i.$$

Adjoint equations

To derive adjoint equations, put Eq. (2) into Eq. (3):

(5)$$\delta _{v}\left({J + \int_{\rm \Omega } {( {u, \;q} ) R{\rm d\Omega }} } \right) + \delta _p\left({J + \int_{\rm \Omega } {( {u, \;q} ) R{\rm d\Omega }} } \right) = 0.$$

For the adjoint equation calculations, turbulence is assumed to be frozen in order to allow for the primal turbulent viscosity to be re-used for the adjoint diffusion term (Othmer et al., 2007). This is a commonly used approximation, called Taylor's hypothesis (Taylor, 1938), that introduces a small source of error in the optimization problem that holds for situations where relative turbulence intensity is small, such as ducted flows.

Further expansion of Eq. (5) produces:

(6)$$\eqalign{ & \delta _vJ + \delta _pJ + \int_\Omega {u[(\delta v\nabla )v + \nabla p-\nabla (2vD(\delta v))]{\rm d}\Omega } -\int_\Omega {q\nabla \delta v{\rm d}\Omega } \cr & + \int_\Omega {u\nabla \delta p{\rm d}\Omega = 0} .} $$

After splitting the cost function J into boundary Γ and domain Ω parts (Lindberg, 2015):

$$J = \int_{\rm \Gamma } {J_{\rm \Gamma }{\rm d\Gamma } + \int_{\rm \Omega } {J_{\rm \Omega }{\rm d\Omega }} }.$$

Equation (6) can be reformulated into:

$$\eqalign{& \int_{\rm \Gamma } {\left({uv + \displaystyle{{\partial J_{\rm \Gamma }} \over {\partial p}}} \right)\delta p{\rm d\Gamma } + } \int_{\rm \Omega } {\left({-\nabla u + \displaystyle{{\partial J_{\rm \Omega }} \over {\partial p}}} \right)\delta p{\rm d\Omega } }\cr &+ \int_{\rm \Gamma } {\left({n( {uv} ) + u( {vn} ) + 2vnD\left({n + \displaystyle{{\partial J_{\rm \Gamma }} \over {\partial v}}} \right)} \right)\delta v{\rm d\Gamma }} -\int_{\rm \Gamma } {2vnD( {\delta {\rm v}} ) u{\rm d\Gamma } } \cr &+\int_{\rm \Omega } {\left({-\nabla uv-vu-\nabla \left({2vD( u ) + \nabla q + \displaystyle{{\partial J_{\rm \Omega }} \over {\partial p}}} \right)} \right)\delta v{\rm d\Omega } = 0.} }$$

To resolve this equation, each function to be integrated has to be equal to 0 individually for any δv and δp value. The parts of the equation that are integrated across the domain Ω give the following adjoint equations.

The adjoint momentum equation, with the sink term included as with the primal equation:

$$-2D( u ) v = {-}\nabla q + \nabla ( {2vD( u ) } ) -\displaystyle{{\partial J_{\rm \Omega }} \over {\partial v}}-\alpha u.$$

The adjoint pressure equation:

$$\nabla u = \displaystyle{{\partial J_{\rm \Omega }} \over {\partial p}}.$$

The adjoint boundary conditions for velocity and pressure:

$$\eqalign{&\int_{\rm \Gamma } {\left({n( {uv} ) + u( {vn} ) + 2vnD( u ) -qn + \displaystyle{{\partial J_{\rm \Gamma }} \over {\partial v}}} \right)\delta v{\rm d\Gamma }}\cr & - \int_{\rm \Gamma } {2vnD( {\delta v} ) u{\rm d\Gamma } = 0,} }$$

$$\int_{\rm \Gamma } {\left({un + \displaystyle{{\partial J_{\rm \Gamma }} \over {\partial p}}} \right)\delta p{\rm d\Gamma } = 0.}$$

Primal boundary conditions

Boundary conditions (BCs) for the primal equations are implemented as presented in Table 1.

Table 1. Primal BCs for the inlet, outlet, and wall for ducted flow problems.

Cost function and adjoint boundary conditions

Fluid power dissipation represents the power consumption experienced by the flow of a fluid through an orifice as a function of the pressure drop and the flow rate. Power dissipation includes the effect of velocity, rather than the option of focusing solely on pressure reduction. It is computed as the net inward flux of total pressure through the device boundaries. The cost function, and therefore the objective, of our optimizer is:

$$J = {-}\int_{\rm \Gamma } {( {\,p + 0.5v^2} ) vn{\rm \;d\Gamma }}.$$

For this cost function, the volumetric contribution $J_{\rm \Omega } = 0$. The adjoint velocity boundary conditions therefore reduce to (Othmer et al., 2007):

$$u_t = 0,$$

$$u_n = \left\{{\matrix{ 0, \hfill & {{\rm at}\,{\rm wall}}, \hfill \cr {v_n}, \hfill & {{\rm at}\,{\rm inlet}}. \hfill \cr } } \right.$$

The outlet adjoint boundary conditions with $\partial J_{\rm \Gamma }/\partial v_n$ fulfilled become:

$$uv + u_nv_n + v( {n\nabla } ) u_n-\displaystyle{1 \over 2}v^2-v_n^2 = q,$$

$$u_tv_n-v_nv_t + v( {n\nabla } ) u_t = 0.$$

Optimization algorithm

Initially, a steepest gradient descent algorithm (Cauchy, 1847) is implemented to update the porosity, the design variable α, using the calculated sensitivities. The steepest descent is a gradient-based algorithm which can be used to update porosity as shown:

$$\alpha _{n + 1} = \alpha _n-\displaystyle{{\partial L\;} \over {\partial \alpha _i}} = \alpha _n-u_iv_iV_i\delta.$$

For the implementation of the steepest descent algorithm, commonly called vanilla gradient descent, δ is the step-length (lambda) and V _i is the cell volume. Modifications to α will employ an underrelaxation factor (γ), to ensure each steps solution is not too different from the last, improving stability.

$$\alpha _{n + 1} = \alpha _n( {1-\gamma } ) + \gamma \min ( {\max ( {( {\alpha_n-u_iv_iV_i\delta } ) , \;0} ) , \;\alpha_{{\rm max}}} ),$$

where α _max represents the maximum porosity value for each cell.

The descent algorithm is then adapted to introduce a stochastic mechanism to compare gradient descent of a fluid topology problem with a custom stochastic gradient descent method. The stochastic mechanism of the optimization takes inspiration from He et al.'s work in which they penalize the sensitivity of each element with a random coefficient (He et al., 2020). For each element i:

$$\hat{\alpha }_i = \alpha _i\cdot \varphi _i, \quad \varphi _{i}\in [ {1.0-\omega , \; 1.0 + \omega } ],$$

where ω is the penalty coefficient, which is randomly chosen for each element before the form-finding process is carried out.

This method is adapted by applying elementwise random penalties to the porosity field α before the sensitivities are calculated, however, the penalty is scaled down to a hundredth in order to promote convergence. At each iteration of the loop, 10 instances of penalized porosity fields are created, to produce a minibatch of size 10 that is computed against the objective function, with the porosity field of the smallest cost being progressed to the next stage. The advantage is that this method will evaluate many different non-optimal geometries allowing it to study more local minima compared with vanilla gradient descent. The effect of this stochastic adaption on the performance of the resulting topology will be compared in the hopes that this method guides the algorithm away from false minima.

Case study

The topology optimization algorithm was applied to the novel case study of an internal combustion engine's piston cooling gallery. This case study is of interest as although it only has one inlet and one outlet, which is the most common use case for topology optimization of fluid systems, the torus region in which flow is split to either side of the gallery induces dynamic flow interactions, such as a large amount of flow impingement at the outlet, and separation zones near the inlet and outlet, which is a flow problem that has not yet been tackled in the literature.

By optimizing the cooling channel of a piston, a greater quantity of cooling medium can be passed through the piston, improving heat transfer. A small increase in fluid flow would be a step toward realizing internal combustion engine designs that can achieve higher thermal efficiencies, as this allows for an increase in engine operating temperature without piston failure. For high-performance engines, such as the Porsche 911 GT2 RS from MAHLE's study (MAHLE, 2020), reduced piston temperature will reduce the onset threshold of knock and pre-ignition, which allows greater performance to be achieved.

Optimization setup

A piston cooling gallery represents a relatively simple yet interesting case for investigation into the efficacy of the approach. The passage consists of one inlet and one outlet, linked by a torus shaped cooling path, meaning that the fluid can flow in two directions, clockwise or anticlockwise, from inlet to outlet. This flow path is of interest as it leads to a large amount of flow impingement at the outlet, which is a challenge to efficiency. The geometry created to represent a simple piston cooling channel can be seen in Figure 2, and this acted as the domain for the optimization, meaning that any optimized design could not exceed these defined boundaries, and all changes must be internal.

Fig. 2. Starting domain of cooling passage within piston.

The cooling passage to be simulated has circular inlets and outlets of 10 mm with the torus passage being 12 mm in diameter. The domain included a small chamfer between the inlet and outlet passages and the torus.

The optimizer was run using an inlet velocity of 15 m/s which represents a Reynolds number of 300,000 for 10 W oil at 100°C (SAE, 1999). The inlet velocity was chosen to simulate a 911 GT2 RS engine with a 77.5 mm engine stroke running at 6000 rpm (Porsche, 2020), to replicate the boundary conditions MAHLE's additive piston would be under.

The domain was discretized using a mesh element size of 0.25 mm. This meant that the optimizer was faced with a domain of 321,724 cell zones. Figure 3 presents a mesh convergence study that indicated a mesh resolution of 0.25 mm offered an equilibrium between good result quality and computational cost. The time increase from elements of 0.25 to 0.1 mm was 7-fold.

Fig. 3. Mesh convergence study for objective function.

As required by the optimizer, the flow is assumed to be incompressible and isothermal. Frozen turbulence was assumed for the adjoint equations, with the primal equations using a k-epsilon turbulence model whilst calculating the Navier–Stokes equations. The simulation's objective function converged after 220 iterations which took 500 s to run on a 4.4 GHz processor.

Optimized geometry

For ease of comparison, Figure 4 presents the initial domain and the optimized geometry side by side. Although the optimizer had limited domain space to alter, significant design changes were realized.

Fig. 4. Initial domain (left) and vanilla gradient descent optimized geometry (right).

An obvious geometry change to improve the performance of the gallery would be to have a feature that directs the flow 90° from the inlet, around the torus, and 90° again to the outlet. This transition geometry was produced for the inlet, as shown in the optimized geometry in Figure 4, however, it was not produced for the outlet. The optimizer did fulfill the expectation that a separating geometry at the inlet would reduce power dissipation, however, the outputted geometry could benefit from a smoother and more arcuate shape which could be added in post-processing.

The most unexpected geometry changed occurred toward the outlet, this being the algorithms implementation of large but thin internal vanes that lead toward the outlet, as displayed in the semi-transparent view of the geometry shown in Figure 6. Figure 5 shows the outlet and the base of the internal section from the outside and Figure 7 shows a side on view of the internal section. This addition goes against general convention that introducing geometry that reduces the cross-sectional area of a tube will decrease and obstruct fluid flow, however, the algorithm has found this geometry to reduce power dissipation. This geometry appears to direct flow toward the outlet, reducing fluid flow impingement at the outlet, and this flow direction change will have overcome the negative effect associated with decreased passage cross-sectional area in order for the solver to have minimized the optimizer's objective function.

Fig. 5. Optimized geometry, internal section to exit interface.

Fig. 6. Semi-transparent view of optimized geometry with internal section.

Fig. 7. Semi-transparent side view of optimized geometry.

The overhead view of the topology, Figure 6, presents a lack of symmetry of the internal vanes. This gives evidence to the gradient descent-based optimizer having failed to find the true global minimum for the problem.

Stochastic addition

The optimization algorithm was adapted to incorporate a stochastic mechanism as described in the methodology section. These changes were made in order to explore the applicability of stochastic gradient descent to improve the algorithm's ability to avoid false minima for the objective function compared with the vanilla gradient descent method.

Figure 8 presents a comparison of the two approaches ability to reduce the cost function for the piston cooling case study. The plot compares gradient descent to stochastic gradient descent highlighting that stochastic gradient descent achieved a lower objective function and therefore should lead to the best performing geometry. Due to the addition of the minibatch system to the stochastic method, the time cost of the stochastic gradient descent algorithm increased to 4500 s, which is a significant 9-fold increase.

Fig. 8. Comparison of objective function between stochastic and vanilla gradient descent.

Figure 9 presents the geometry produced by the stochastic gradient descent. Compared with the vanilla gradient descent's topology output, the most notable differences are the introduction of a separating geometry at the outlet, increased symmetry and shorter vanes leading toward the outlet. The addition of the flow directing geometry at outlet and the improved symmetry hint at the discovery of an improved minima, but performance analysis will give quantitative evidence for this claim.

Fig. 9. Stochastic gradient descent geometry.

Performance

The performance of the vanilla gradient descent and stochastic gradient descent was analyzed and compared with the starting domain. This was done through commercial CFD package ANSYS Fluent.

The boundary conditions for these simulations were chosen to represent 10 W 40 oil, which is a commonly used artificial engine oil with a density of 865 kg/m³ and a viscosity of 0.0143 kg/ms at 100°C (SAE, 1999). The inlet velocity was 15 m/s representing a Reynolds number of 300,000 as used in the optimization. The SIMPLE method was used for the solution, as used in the optimization, along with a k-epsilon viscous model, with wall enhancement, as used in the optimization. The oil was assumed incompressible. The effect of crankcase pressure on the efficiency of the piston cooling channel was ignored during simulation.

The velocity streamlines for the starting domain and the optimized geometry can be seen through Figures 10–12.

Fig. 10. Velocity streamlines for the initial starting domain.

Fig. 11. Velocity streamlines for the vanilla gradient descent topology.

Fig. 12. Velocity streamlines for the stochastic gradient descent topology.

It is clear that the two optimized geometries have reduced flow impingment over the starting domain. The gradient descent displays less impingment at inlet but more toward the outlet compared with the stochastic result. The geometry change at the inlet has resulted in improved fluid velocity throughout the stochastic and vanilla gradient descent topologies. Both optimized geometries exhibit fluid particle pathlines that take a larger radius route around the torus curve, compared with the original domain, which represents more efficient flow, with less energy being lost due to sharp direction changes.

Table 2 presents the resulting area weighted average pressure drop from inlet to outlet. The pressure drop through the cooling passage was 8.9% less in the optimized geometry compared with the initial domain. This is a sizeable performance improvement, especially given the constraints the starting domain would have put on design options, as it has been shown that piston cooling velocity has a linear relation to the average heat transfer coefficient at the bottom of a piston (Najafabadi et al., 2014). The stochastic approach generated a further improvement of 1.1% on the vanilla gradient descent geometry which coincides with the improvement it saw in its objective function. A performance increase greater than 1% can be considered a useful step to improving efficiency and reducing emissions, for example Luff et al. found that just a 3°C drop in piston temperature decreases CO output by 1% (Luff et al., 2012).

Table 2. Pressure drop through initial domain and optimized geometries.

Figure 13 compares the outlet velocities for both geometries, showing an increase in outlet velocity. The stochastic method can be seen to have produced a more uniform outlet distribution, accountable to a more effective exit vane structure. The reduced pressure loss and increased velocity mean that the optimized geometry will more efficiently allow for the passage of fluid, and thus improve the cooling capability of the piston head.

Fig. 13. Outlet velocity contour for initial domain (left), vanilla gradient descent geometry (center), and stochastic gradient descent geometry (right).

The most unusual feature of minima that the algorithm found was the internal vanes that lead toward the outlet, Figure 14 compares the stochastic and vanilla gradient descent topologies and the effect these sections have on the flow path.

Fig. 14. Outlet internal vane geometry streamlines for the vanilla gradient descent method (top) and the stochastic gradient descent method (bottom).

The streamline plots show that the stochastic gradient descent topology has further reduced flow impingement at the outlet, which is the main driver in performance improvement. The shorter vane has produced a more coherent flow structure reducing inefficiencies and thus proving the algorithm has better avoided false minima. For both geometries, the primary sensitivity reducing aspect of this unusual vaned section appears to be to lead the flow by a larger radius toward the outlet and reduce the impingement with the opposing side of the torus.

Further discussion

The resultant geometries produced by both algorithms exhibit creative and unintuitive solutions. The application of either methodology could be employed by engineers in the early stages of product development and, by combining machine and human, produce a better product than could be created by each counterpart alone. The internal fins are testament to this, as it goes against intuition to design a geometry, with the goal of improving ducted flow efficiency, that decreases the products cross-sectional area. The algorithm has proven merit to exploring an unusual design, with the internal fins proving to be a useful mechanism to improve flow transition into the outlet.

The addition of this specific stochastic mechanism to the optimization problem improved results, however this came at a significant time cost. With the optimization of a fluid–solid interface representing a nonconvex problem, the application of vanilla gradient descent is very prone to missing global minima. Although there is a time penalty associated with this stochastic method, it represents a straightforward addition to increase the probability of the algorithm finding the global minima. From the case study, it is clear that the stochastic method tended closer to the global minima due to improved performance. Both approaches have produced a geometry that is closer to optimum than the base case however neither have truly found the global minima due to both displaying an amount of asymmetry, as the true global minima would be fully symmetric due to the starting domain being fully symmetric. To further quantify the capability of the stochastic algorithm, further work should compare the algorithm with a design that is known to be optimum, or close to optimum, given the same boundary conditions.

The time cost increase associated with the stochastic mechanism are largely due to the minibatch addition; it is possible that reducing the relaxation on the random penalties could lead to faster convergence, however, this should be handled in an adaptive way, increasing relaxation with iterations. Compared with the vanilla gradient descent method, the stochastic additions greatly reduce the agility of the tool to produce fast design improvements, however the increased time has led to quantifiable improvements in results. A lot of post-processing would be required if the designer were to make either algorithm's resulting geometry manufacturable, for example, the internal vanes would need shortening and thickening to meet minimum wall thickness requirements, and radii on several geometry transitions would likely need to be increased. This is a drawback to the algorithm, as the changes required for manufacture are likely to reduce the performance of the design, and this represents a necessary recommendation for further algorithm adaptations. Further work into the best approach to manufacture topology optimized fluid systems could look at using the optimizer to define an envelope of variation for different regions of the geometry, as some areas of the geometry will have reduced sensitivity regarding the cost function and thus could be changed to aid manufacture.

Conclusion

We have introduced a novel stochastic gradient descent optimization method to the optimization of fluid dynamics systems. The case study has shown the benefits derived from both the vanilla and stochastic gradient descent algorithms for improving designs. From the starting domain, we achieve a 10% reduction in pressure loss through the stochastic method, which is a useful step toward improved efficiency of an internal combustion engine. A post-processed version of the geometry produced by the optimizer could easily be added to the pistons MAHLE has produced through additive manufacture, furthering the performance of this product, representing a clear example of the use of the optimizer in the real world.

This work has highlighted the effectiveness of applying stochastic methods to the nonconvex problem of fluid–solid surface optimization, presenting a method that should act as a platform for further research to develop on. The creative results presented in the case study have demonstrated the power of this approach as a tool in the engineer's toolbox. There is a clear position for the algorithm in the development process of new products, and it has shown to be able to improve the performance of a piston's cooling channel through design changes that are less obvious.