Topology Optimization
Introduction
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

$$M(u,\gamma )$$
 Subject to

$$D(u,\gamma )\ge 0$$
where $u=({u}_{1},{u}_{2},{u}_{3},p,\mathrm{...})$ is the solution vector of flow variables ${u}_{1}$ for velocity vector component; $p$ for pressure, and $\gamma $ is the design topology field which determines the presence of solid material. The function $M(u,\gamma )$ is the objective, $D(u,\gamma )$ and ${G}_{k}(u,\gamma )$ are constraint functions. The flow equations are represented by $R(u,\gamma )=0$ acting as an additional constraint.
Here, $M(u,\gamma )$ represents the mechanical losses as a single objective, integrating over the inlet and outlet surfaces.
where $\rho $ is the density, $S$ is the surface.
The function $D(u,\gamma )=D(\gamma )$ 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 $\Omega $ is the volume.
As an alternative, the constraint can also be given as an upper bound constraint $D(u,\gamma )\le 0$ .
Further, the functions ${G}_{k}(u,\gamma )$ are used to constrain mass flow rates at outlets 𝑘 for the 𝑁 outlets,
where ${m}_{k}$ is the mass flow rate at exit 𝑘 and ${\overline{m}}_{k}$ 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.
 Minimize

$$M(u,\gamma )+{c}_{s}{\epsilon}_{n}D\left(\gamma \right)$$
 Subject to

$$R(u,\gamma )=0$$
Flow Equations
The state equation residuals of the incompressible equations are defined as, ${R}_{ui}=0$ , ${R}_{p}=0$ where
where $p$ , ${u}_{i}$ denote the static pressure and velocity components, respectively and 𝜌 denotes density and ${\widehat{b}}_{i}$ are the specific body force components. The stresses ${\tau}_{ij}$ are given by
where $\mu $ 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.
The Darcy porosity model for the momentum flow equations is used; the residual is expressed as
where ${f}_{i}$ are the components of the porous media force vector
and ${C}_{\text{Darcy}}$ is the Darcy coefficient,
The value of ${C}_{\text{Darcy}}^{s}$ 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 ${f}_{i}$ to avoid poor convergence.
The adjoint equations to the linearized incompressible NavierStokes equations are given by,
where the adjoint variables ${u}^{a}$ and ${p}^{a}$ measures the sensitivities of varying the field $\gamma $ around the fluid solution $u$ and $p$ of the forward problem.
Optimization Algorithm
The variation of the mechanical losses objective function (Equation 5) with respect to the design field variables is given by the expression
where ${u}_{i}^{a}$ is the adjoint velocity vector. The coefficient ${\tilde{C}}_{\text{Darcy}}\left(\gamma \right)$ is the variation of the darcy factor, which is constant in the design domain for the SIMPmethod 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 $\tilde{\gamma}$ .
The function ${G}_{k}\left(u,\gamma \right)$ given by Equation 7, used to account for the constraint on mass flow at outlet $k$ , is specified as an equality constraint and actually treated as a double inequality constraint defined below as
or
where ${\delta}_{k}^{\text{Tol}}$ is a userdefined 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.
The design variable 𝛾 is continuous, but a welldefined 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 welldefined solution of either fluid or solid, a quality index can be used,
where 0 ≤ 𝐼𝑄 ≤ 1, and 𝐼𝑄 = 1 for a welldefined geometry. Here, 𝑉 is the volume of the design domain.
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.
Examples
 Minimize

$$M(u,\gamma )$$
 Subject to

$$R(u,\gamma )=0$$
 Minimize

$$M(u,\gamma )+{c}_{s}{\epsilon}_{n}D(\gamma )$$
 Subject to

$$R(u,\gamma )=0$$
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%.
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 ×10^{5} 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.
A body fitted mesh with 6.5×10^{5}nodes 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 spacefilled and all the parts of the design are taken into consideration while performing the shape optimization step.
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.
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.20211892.
[2] Svanberg, K., “MMA and GCMMAtwo methods for nonlinear optimization,” vol, Vol. 1, 2007, pp. 1–15.