Topology optimization for fluid flow in AcuSolve is based
on the solution of the flow equations combined with the solution of associated
linearized adjoint equations. The objective is to minimize the mechanical losses of
the fluid flow system in combination with constraints limiting the used material. To
represent the solid geometry a field variable is used to indicate the presence of
material which adds a resistance modeled by a porous media force with a positive
Darcy coefficient. Refer to Sandboge et al. (2021) for details.
This document shows the equations used to solve the topology optimization problem for
flow.
Optimization
Minimize mechanical
losses
The considered optimization problems take the form:
Minimize
Subject to
where is the solution vector of flow
variables for velocity vector
component; for pressure, and
is the design topology field which
determines the presence of solid material. The function
is the objective,
and
are constraint functions. The flow
equations are represented by acting as
an additional constraint.
Here, represents the mechanical losses as a single
objective, integrating over the inlet and outlet surfaces.
where is the density,
is the surface.
The function is here the function of design fraction occupying
the solid volume and is used as an optional constraint to limit the amount of solid
material.
where is the volume.
As an alternative, the constraint can also be given as an upper bound constraint
.
Further, the functions are used to constrain mass flow rates at outlets 𝑘
for the 𝑁 outlets,
where is the mass flow rate at exit 𝑘
and is the target mass flow rate at the
same exit.
The constraints are naturally optional. It is possible to have an optimization
problem without any constraints, except for the flow constraint Equation 4 which is always active.
Minimizing the mechanical losses and fluid
volume
If represents the mechanical losses in Equation 1, regions of zero velocity will not change the objective
value. Thus, zero velocity "pockets" in the solid region can occur in the final
solution. To avoid such unclean solid geometry, it is possible to define as a multi-objective adding the size of the fluid
volume in the design domain to the mechanical losses. The strength of minimization
of the fluid volume size is defined by a user-defined coefficient multiplied by a normalization coefficient ,
Minimize
Subject to
Flow Equations
Computational domain
The flow equations are solved in a computational domain which is allowed to vary
during the optimization process by modifying the material properties. At the first
iteration, the design domain may use material properties representing a fluid, such
as air, see Figure 1, and at convergence some parts of the domain have built
up regions with material properties representing a solid material, see Figures Figure 2–Figure 3. The computational domain is fixed, but the material properties are functions
of a design field 𝛾 as is explained in the continuous adjoint for topology
optimization section.
Incompressible Navier-Stokes
equations
The state equation residuals of the incompressible equations are defined as, , where
where , denote the static pressure and velocity components,
respectively and 𝜌 denotes density and are the specific body force components. The stresses are given by
where is the viscosity.
Continuous Adjoint for Topology Optimization
The optimization problem (Equation 1) is solved using the adjoint methods. The continuous
adjoint approach for topology optimization is presented in this section.
Momentum equations for topology
optimization
The Darcy porosity model for the momentum flow equations is used; the residual is
expressed as
where are the components of the porous media force
vector
and is the Darcy coefficient,
where the Solid Isotropic Material with Penalization (SIMP) defined by , a constant whose value is selected to . Thus, the design field determines the level of resistance in the design
domain; zero resistance for fluid regions and a high resistance to minimize the flow
through the solid regions.
The value of is automatically computed
by the optimization algorithm, aiming to minimize the flow through the solid region
but at the same time limiting the size of the force
to avoid poor convergence.
Adjoint equations
The adjoint equations to the linearized incompressible Navier-Stokes equations are
given by,
where the adjoint variables and
measures the sensitivities of
varying the field around the fluid
solution and
of the forward problem.
Optimization Algorithm
Derivative of responses
The variation of the mechanical losses objective function (Equation 5) with respect to the design field variables is given by the
expression
where is the adjoint velocity vector. The coefficient is the variation of the darcy factor, which is
constant in the design domain for the SIMP-method of order 1. In essence, the
derivative of the design material is determined from the scalar product between the
velocity and the adjoint velocity vectors.
A volume constraint based on the topology design fraction in Equation 6 gives the variation
which gives a constant value in the design domain for the same volume variation .
The function given by Equation 7, used to account for the constraint on mass flow at outlet , is specified as an equality constraint and actually
treated as a double inequality constraint defined below as
or
where is a user-defined tolerance value for the exit 𝑘.
The variation of this function is equal to
where 𝑆𝑂𝑘 is the surface for outlet 𝑘.
Given the values of the objectives and constraints, and the field derivatives with
respect to 𝛾 for the objective, (Equation 18 for mechanical losses), and constraints, (constant
for the design fraction constraint), the optimizer can update the design field 𝛾.
The default optimizer is the methods of moving asymptotes (MMA). Refer to Svanberg
(2007) for details of the MMA.
Metrics
The design variable 𝛾 is continuous, but a well-defined geometry should have the
value either 𝛾 = 0 for a fluid or 𝛾 = 1 for a solid with small regions of the
design domain where 0 < 𝛾 < 1. To measure the amount a well-defined solution
of either fluid or solid, a quality index can be used,
where 0 ≤ 𝐼𝑄 ≤ 1, and 𝐼𝑄 = 1 for a well-defined geometry. Here, 𝑉 is the volume
of the design domain.
Illustration of the optimization
algorithm
Equation 14 provides information on how the design material should be
altered in the design space at each design iteration. To determine if solid material
should be added or removed, it is possible to view the velocity 𝑢 and adjoint
velocity 𝑢𝑎 vectors in conjunction; if the angle is less than 90
degrees, material should be removed and if the angle is greater than 90 degrees
material should be added.
Figure 4 shows (i) the initial velocity and adjoint velocity
fields where the vectors have an angle greater than 90 degrees in some regions in
front and rear of the obstacle; and (ii) the final velocity and adjoint velocity
fields where all angles are less than 90 degrees in the entire fluid design
domain.
A close-up of the region behind the obstacle is shown in Figure 5, where the velocity shows some circulation while the
adjoint velocity is attached without circulation.
Examples
Duct with sharp turn
As a simple example consider the duct domain (Figure 6), where the following optimization problem is solved,
Minimize
Subject to
where 𝑀(𝑢,𝛾) is the function of Equation 5, representing the mechanical losses of the fluid
flow.
To solve the multi objective problem where you also factor in the size of the
occupied fluid volume you have,
Minimize
Subject to
By using a larger coefficient 𝑐𝑠, the duct cross-section will tend to be
smaller due to a higher weight minimizing the occupied fluid volume
().
Design-optimization of an HVAC Ducting
System
The optimization technology is illustrated in this section in the design of the HVAC
system of the Altair CX1 car. The cabin geometry of the
Altair car is illustrated in Figure 8.
In this figure the ten different inlets to the cabin are shown. The six of them are
aiming at the front and rear upper part of the cabin and the remaining four are
providing the air to the front and rear lower part of the cabin. The supply of the
air from the blower has also been separated into two parts. The first supply has to
be connected with the upper inlets and the second supply has to be connected to the
lower inlets to the cabin.
The first part of the algorithmic procedure is based on topology optimization for
minimum power dissipation or mechanical energy losses to account for minimum fuel
consumption and minimum volume to account for minimum manufacturing cost constrained
by specific mass flow rates at each exit. In this study all mass flow rates selected
as constraints were equal to 0.001kg/s and the tolerance of the constraints were
selected equal to 5%.
Given the available space from the cabin and generally the car components, the first
step is to define the initial box within which the duct must be confined. Apart from
the availability in space the only other restriction to apply for topology
optimization is to define the inlet to the box as the air supply location from the
blower(s) and the outlet from the box as the inlets to the cabin. The CAD geometry
can be in general automatically created and can be extracted from the available
space. In this case a set of simple boxes were created to define the initial domain
(Figure 9).
The next step is to generate the computational mesh at the initial domain. This mesh
can consist of either hexaedral or tetrehedral elements. For this study, a
hexahedral mesh with 2.5 ×105 nodes is used, Figure 10, left, which has sufficient resolution for topology
optimization.
The topology optimization of the two ducts was performed in a sequential manner.
First, the duct connecting the one inlet with the six upper outlets (cabin inlets)
was optimized. For that reason, a part of the domain was excluded from the solution
at this step to allow the evolution of the design of the first duct in a way that
would not block the evolution of the other duct. This exclusion was made possible by
applying nodal boundary conditions on the design variable at the excluded part of
the domain and not through separate meshing.
After, the first optimal duct geometry is computed, the space it requires is excluded
from the initial box, using again nodal boundary conditions on the design variable.
To avoid intersections and very close proximity between the two duct geometries a
user selected area that surrounds the first duct is also excluded. The second duct
is optimized, and the optimal geometry uses the same computational mesh. Thus, a
single initial mesh should only be created and no remeshing is required while
performing the optimization of the separate ducts. This procedure can be continued
to a next round of optimization of the two ducts; the first duct can be optimized
again, leaving out not the initial random space but the space covered by the second
duct (plus an additional layer), then the second duct can be optimized again based
on the new shape of the first duct and so forth. This approach is expected to result
in designs with even better performance but was omitted in this study.
The final design of the two ducts using topology optimization is afterwards extracted
and smoothed. The final shape is shown in Figure 10, right.
A body fitted mesh with 6.5×105nodes is generated at this geometry in
order to further optimize the shape using the shape optimization methodology with
morph shapes. Shape optimization is used to further improve the value of the power
dissipation objective function and ensure that the constraints on mass flow rates at
the exits of the domain are as accurately as possible satisfied. It has been
observed that the values of the mass flow rates at the exits are very sensitive to
the geometrical design and mesh and, thus, the smoothing and remeshing of the
optimized geometry affects them quite a lot. The body fitted mesh is also expected
to give more accurate estimation of the quantities of interest and, thus shape
optimization is expected to give the completely final design of the ducts. It should
also be mentioned that the two ducts are handled as one during this step to avoid
intersections and high proximity in the final design.
The morph shapes, required for the shape optimization described in the theory, are
generated automatically using the flow field computed in the body fitted mesh. This
way the two methodologies, topology and shape optimization described above are not
disconnected but rather constitute the ingredients of a unified methodology. The
flow velocity is the metric used to define optimally the morph shapes so that
several types of morph shape can be formed as well as to ensure that the whole
domain is appropriately space-filled and all the parts of the design are taken into
consideration while performing the shape optimization step.
The automatically computed centers of the morph shapes are shown in Figure 11, left. In this case, the length scale of morph shape
generator was equal to 0.2m which resulted in the generation of 22 control points in
total. In Figure 11, right, the morph shape vectors are shown.
Two different types of morph shapes are generated; the first one is depicted by the
red vectors, and it allows for the bending of the duct during optimization. Two
design variables per node are needed for this type of morph shape to cover all the
design space. The second type of morph shapes is illustrated by the blue arrows and
is used to increase or decrease the cross section of the duct geometry. One design
variable per node is needed for this morph shape which includes four or eight
control points with radial vectors that allow for the change in the cross section.
The total number of morph shapes and thus the number of design variables for shape
optimization is 66.
Based on the automatically computed morph shapes, shape optimization is performed
further minimizing the power dissipation, with the same constraints on the mass flow
rates at the exits and constraint on the volume of the design to be constant with a
small tolerance. The initial and optimal geometries during this stage of
optimization are shown in Figure 12. The distribution of the magnitude of the grid
displacement that corresponds to the optimal geometry is also shown. The mass
balance constraint at the exits is satisfied at this stage with a tolerance of
1%.
The optimal design of the HVAC with streamlines and velocity vectors is shown in
Figure 13, together with the continuing velocity vectors inside the
cabin.
References
[1] Sandboge, R., Papadimitriou, D., and Reddy, M., “Shape and Topology
Optimization in Computational Fluid Dynamics Including Heat Transfer Using
Gaussian Processes and Adjoint Methods,” 2021.
https://doi.org/10.2514/6.2021-1892.
[2] Svanberg, K., “MMA and GCMMA-two methods for nonlinear optimization,” vol,
Vol. 1, 2007, pp. 1–15.