## 2 Methodology
### Governing Equations  
We seek the optimum body shapes in the flow. To achieve that, we first introduce \(\gamma\in\{0,1\}\) as the design variable, which is a binary distribution defined in the design domain. Respectively, \(0\) and \(1\) represent the solid phase and liquid phase, corresponding to the presence or absence of material at an arbitrary point within the design domain. In this problem, the distribution of the solid material is the optimal body structure we expect to obtain within the flow field, so the essence of the topology optimization is an optimization problem with \(0\)-\(1\) discrete variables. In order to adopt an optimization algorithm based on continuous variables and avoid the difficulties of solving binary optimization problems, we can relax \(\gamma\) to [0, 1] to achieve continuous variation. Figure 1 is a schematic representation of the fluid topology optimization problem. The shape optimization problem is defined in a bounded domain \(\Omega\). \(\Omega\) is the fluid domain and \(\Omega_{\text{S}}\) is the design domain; both of them are separated by the boundary \(\Gamma_{\text{S}}\). The above model in a steady incompressible laminar flow (ignoring gravity effect) can be described by the stationary incompressible Navier-Stokes equation with boundary conditions as:  
\[\left\{ \begin{aligned} \rho\mathbf{u}\cdot\
abla\mathbf{u}+ \
abla\cdot\left[-\mu(\
abla\mathbf{u}+\
abla^{T}\mathbf{u})+p\mathbf{I} \right]&=\mathbf{f}\quad\text{in}\;\Omega\\ -\
abla\cdot\mathbf{u}&=0\quad\text{in}\;\Omega\\ \mathbf{u}&=u_{0}\quad\text{on}\;\partial\Omega^{in} \\ \left[-\mu(\
abla\mathbf{u}+\
abla^{T}\mathbf{u})+p\mathbf{I} \right]\cdot\mathbf{n}&=0\quad\text{on}\;\partial\Omega^{out}\\ \mathbf{u}&=0\quad\text{on}\;\Gamma_{s}\\ \mathbf{u}\cdot\mathbf{n}&=0\quad\text{on}\;\Gamma \end{aligned} \right. \tag{1}\]  
where \(\rho\) is the fluid density, \(\mu\) is the viscosity, \(\mathbf{I}\) is the unit tensor, \(\mathbf{f}\) is the body force vector (source term), and \(\mathbf{n}\) is the unit outward normal vector on the boundary of flow domain. The fluid velocity vector \(\mathbf{u}\) and pressure \(p\), which characterize the fundamental motion of the particles in the flow field, are the only solutions to the system (1). In order to solve the Navier-Stokes equation, boundary conditions are necessary to set in Equation (1), where the normal stress at the outlet boundary \(\partial\Omega^{out}\) is zero; initial velocity at the inlet boundary \(\partial\Omega^{in}\) is a constant \(\mathfrak{u}_{0}\); the velocity at the fluid-solid interface \(\Gamma_{s}\) is zero; and \(\Gamma=\partial\Omega\backslash\left(\partial\Omega^{in}\cup\partial\Omega^{out}\right)\) denotes the free-slip boundaries. To facilitate the calculation, we introduce the reference length L and the reference velocity U to make all variables in a dimensionless form as:  
\[\mathbf{u}=\frac{\mathbf{u}}{U},\;p=\frac{pL}{\mu U},\;\mathbf{f}=\frac{ \mathbf{f}L^{2}}{\mu U},\;\
abla=L\
abla\;\;\text{Re}\equiv\frac{\rho LU}{\mu} \tag{2}\]  
where Re is the Reynolds number defined as the ratio of the inertial force to the viscous force. The dimensionless momentum equation is as follows:  
\[\text{Re}\mathbf{u}\cdot\
abla\mathbf{u}+\
abla\cdot\left[-(\
abla\mathbf{u} +\
abla^{T}\mathbf{u})+p\mathbf{I}\right]=\mathbf{f}\;\text{in}\;\Omega \tag{3}\]  
For using adjoint analysis later, we give a variational description of Navier-Stokes equations as follows:  
Figure 1: Schematic of the fluid flow problem. \[\text{Find}\left\{