<figure><img src="image_1.png" /></figure>

HAL
open science

# Topology optimization for 3D fluid diode design considering wall-connected structures

Takamitsu Sasaki, Kozo Furuta, Naoyuki Ishida, Hao Li, Tsuguo Kondoh, Kazuhiro Izui, Shinji Nishiwaki

<figure><img src="image_2.png" /></figure>

## To cite this version:

Takamitsu Sasaki, Kozo Furuta, Naoyuki Ishida, Hao Li, Tsuguo Kondoh, et al., Topology optimization for 3D fluid diode design considering wall-connected structures. 2024. hal-04618128

HAL Id: hal-04618128

https://hal.science/hal-04618128

Preprint submitted on 20 Jun 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Topology optimization for 3D fluid diode design considering
wall-connected structures

Takamitsu Sasaki¹*, Kozo Furuta²†, Naoyuki Ishida¹†, Hao Li¹†,
Tsuguo Kondoh¹†, Kazuhiro Izui²†, Shinji Nishiwaki¹†

¹*Department of Mechanical Engineering and Science, Kyoto University, Kyoto Daigaku Katsura, Nishikyo-ku, Kyoto, 615-8530, Japan.

²Department of Micro Engineering, Kyoto University, Kyoto Daigaku Katsura, Nishikyo-ku, Kyoto, 615-8530, Japan.

*Corresponding author(s). E-mail(s): sasaki.takamitsu.36c@st.kyoto-u.ac.jp;
†These authors contributed equally to this work.

Abstract

This paper proposes a density-based topology optimization method for the three-dimensional design of fluid diodes considering wall-connected structures based on the fictitious physical modelling approach. The optimum design problem of fluid diodes is formulated as maximizing the energy dissipation in the reverse flow subject to the upper bound constraint of the energy dissipation in the forward flow. A fictitious physical model and a geometric constraint are constructed to detect and restrict the “floating” solid domains, which are not connected to the outer boundaries. The sensitivities of cost functions are derived and computed based on the continuous adjoint method. The finite volume method is employed to discretize the governing and adjoint equations to mitigate the huge computational costs of three-dimensional fluid analysis. Numerical investigations are presented to validate the fictitious physical model and the geometric constraint for excluding “floating” islands. Finally, topology optimization for fluid diodes with and without the geometric constraint is performed, and the result demonstrates that the proposed method is capable of generating fluid diodes with wall connectivity while maintaining a good functional performance.

**Keywords:** Topology optimization, Fluid diodes design, Wall-connected structures, Geometric constraint, Finite volume method

# 1 Introduction

Fixed-geometry fluid diodes or no-moving-part valves (e.g. Tesla valves) are check valves without moving parts or deformation whose resistance depends on flow directions. The geometry of fluid diodes is designed so that the fluid can flow with little resistance in the forward direction while fluidic inertial forces generate huge resistance in the reverse direction. Since fluid diodes without

any moving parts are easy to fabricate and resistant to wear and fatigue, they are widely used as components of micropumps [Forster et al., 1995, Morganti et al., 2005, Nabavi, 2009].

A number of studies have been conducted to improve the performance of fluid diodes. Most of the previous studies have employed a measure of diodicity, defined as the ratio of the pressure drop in the reverse flow to that in the forward flow,

001

002

003

004

005

006

007

008

009

010

011

012

013

014

015

016

017

018

019

020

021

022

023

024

025

026

027

028

029

030

031

032

033

034

035

036

037

038

039

040

041

042

043

044

045

046

047

048

049

050

051
# Topology optimization for 3D fluid diode design considering wall-connected structures

Takamitsu Sasaki¹*, Kozo Furuta²†, Naoyuki Ishida¹†, Hao Li¹†, Tsuguo Kondoh¹†, Kazuhiro Izui²†, Shinji Nishiwaki¹†

¹*Department of Mechanical Engineering and Science, Kyoto University, Kyoto Daigaku Katsura, Nishikyo-ku, Kyoto, 615-8530, Japan.

²Department of Micro Engineering, Kyoto University, Kyoto Daigaku Katsura, Nishikyo-ku, Kyoto, 615-8530, Japan.

*Corresponding author(s). E-mail(s): sasaki.takamitsu.36c@st.kyoto-u.ac.jp;

†These authors contributed equally to this work.

## Abstract

This paper proposes a density-based topology optimization method for the three-dimensional design of fluid diodes considering wall-connected structures based on the fictitious physical modelling approach. The optimum design problem of fluid diodes is formulated as maximizing the energy dissipation in the reverse flow subject to the upper bound constraint of the energy dissipation in the forward flow. A fictitious physical model and a geometric constraint are constructed to detect and restrict the “floating” solid domains, which are not connected to the outer boundaries. The sensitivities of cost functions are derived and computed based on the continuous adjoint method. The finite volume method is employed to discretize the governing and adjoint equations to mitigate the huge computational costs of three-dimensional fluid analysis. Numerical investigations are presented to validate the fictitious physical model and the geometric constraint for excluding “floating” islands. Finally, topology optimization for fluid diodes with and without the geometric constraint is performed, and the result demonstrates that the proposed method is capable of generating fluid diodes with wall connectivity while maintaining a good functional performance.

**Keywords:** Topology optimization, Fluid diodes design, Wall-connected structures, Geometric constraint, Finite volume method

## 1 Introduction

Fixed-geometry fluid diodes or no-moving-part valves (e.g. Tesla valves) are check valves without moving parts or deformation whose resistance depends on flow directions. The geometry of fluid diodes is designed so that the fluid can flow with little resistance in the forward direction while fluidic inertial forces generate huge resistance in the reverse direction. Since fluid diodes without

any moving parts are easy to fabricate and resistant to wear and fatigue, they are widely used as components of micropumps [Forster et al., 1995, Morganti et al., 2005, Nabavi, 2009].

A number of studies have been conducted to improve the performance of fluid diodes. Most of the previous studies have employed a measure of diodicity, defined as the ratio of the pressure drop in the reverse flow to that in the forward flow,

001

002

003

004

005

006

007

008

009

010

011

012

013

014

015

016

017

018

019

020

021

022

023

024

025

026

027

028

029

030

031

032

033

034

035

036

037

038

039

040

041

042

043

044

045

046

047

048

049

050

051
as a measure of the performance of fluid diodes. Truong and Nguyen [2003] investigated the design optimization parameters for fluid diodes such as the length and the cross angle of channels and conducted a parametric optimization. Gamboa et al. [2005] applied shape optimization and obtained the geometry with significantly improved diodicity over the Reynolds number $0 < Re \le 2000$. Although these studies have achieved improvements in fluid diode performance, the topology of target geometry does not change during the optimization process in these methods.

Topology optimization has advantages over those optimization methods in that the topology of the structure can be changed during the optimization process, leading to a significant performance improvement. In this method, the structural optimization problem is replaced with a material distribution problem, and optimal distribution is obtained based on mathematical and physical principles. Topology optimization originated in solid mechanics [Bendsøe and Kikuchi, 1988] and has since been extended to other physical phenomena. The pioneer work of topology optimization in fluid problems was done by Borrvall and Petersson [2003], in which Stokes equations were considered with a density-based method introducing Darcy law (also called the Brinkman model). Based on this approach, topology optimization has been extended to Navier-Stokes equations [Gersborg-Hansen et al., 2005, Olesen et al., 2006, Zhou and Li, 2008], non-Newtonian [Pingen and Maute, 2010], unsteady [Kreissl et al., 2011, Deng et al., 2011], and turbulent flows [Kontoleontos et al., 2013, Yoon, 2016, Dilgen et al., 2018], and applied to a variety of design problems, e.g., optimum design of rotor machines [Romero and Silva, 2014], aerodynamics [Li et al., 2022b, Ghasemi and Elham, 2022], heat sinks [Yaji et al., 2015, Li et al., 2022a] in three dimensions.

Topology optimization was also applied to design problems of fluid diodes [Deng et al., 2010, Lin et al., 2015]. Lin et al. [2015] formulated the optimal design problem for fluid diodes as maximizing diodicity subject to a volume constraint with an additional penalty on flow through intermediate density. Sato et al. [2017a] proposed a bi-objective topology optimization method for fluid diodes and investigated Pareto-optimal solutions. In work done by Guo et al.

[2020], length scale constraints were introduced into two-dimensional topology optimization for the design of fluid diodes to reduce the difference between the performance of the two-dimensional model and the three-dimensional model obtained from the extrusion.

Although these previous studies showed that topology optimization can generate novel designs of fluid diodes with improved performances, one limitation in the design problems is the two-dimensional assumption of fluid flow. This assumption is valid as long as the effect of upper and lower boundaries is negligible. In practice, however, the height of the flow channel is not infinite, and the discrepancy between the performances of the two-dimensional model and the extruded three-dimensional model should not be neglected. Lin et al. [2015] reported that the diodicity values of the extruded three-dimensional model of fluid diodes measured by both numerical simulations and experiments were less than those of the two-dimensional model, and they concluded that topology optimization in three dimensions should be performed when the height of the channel cannot be neglected. Although Guo et al. [2020] showed that length scale control can mitigate the discrepancy, difficulty in estimating the final performance still remains. Thus, this work sheds a light on the topology optimization of fluid diodes in full-scale three-dimensional setting.

However, topology optimization of fluid diodes in three dimensions poses another challenging issue. The optimized structures may feature some “floating” islands disconnected to the external walls, which is not ideal in practice, therefore, cannot be manufactured. The Brinkman model allows the nucleation of such structures owing to the nature of its mathematical modelling. Therefore, we need to construct a design methodology for obtaining wall-connected structures in the fluid topology optimization.

Previous works have investigated geometric constraints for ensuring desired geometrical features of the structure, especially in the fields of solid mechanics and additive manufacturing. In topology optimization, the control of minimum and maximum length scales of structures has long been a typical issue of geometric constraints [Poulsen, 2003, Guest et al., 2004, Zhou et al., 2015, Guest, 2009, Li et al., 2023], and the concept of geometric constraints have been
as a measure of the performance of fluid diodes. Truong and Nguyen [2003] investigated the design optimization parameters for fluid diodes such as the length and the cross angle of channels and conducted a parametric optimization. Gamboa et al. [2005] applied shape optimization and obtained the geometry with significantly improved diodicity over the Reynolds number $0 < Re \le 2000$. Although these studies have achieved improvements in fluid diode performance, the topology of target geometry does not change during the optimization process in these methods.

Topology optimization has advantages over those optimization methods in that the topology of the structure can be changed during the optimization process, leading to a significant performance improvement. In this method, the structural optimization problem is replaced with a material distribution problem, and optimal distribution is obtained based on mathematical and physical principles. Topology optimization originated in solid mechanics [Bendsøe and Kikuchi, 1988] and has since been extended to other physical phenomena. The pioneer work of topology optimization in fluid problems was done by Borrvall and Petersson [2003], in which Stokes equations were considered with a density-based method introducing Darcy law (also called the Brinkman model). Based on this approach, topology optimization has been extended to Navier-Stokes equations [Gersborg-Hansen et al., 2005, Olesen et al., 2006, Zhou and Li, 2008], non-Newtonian [Pingen and Maute, 2010], unsteady [Kreissl et al., 2011, Deng et al., 2011], and turbulent flows [Kontoleontos et al., 2013, Yoon, 2016, Dilgen et al., 2018], and applied to a variety of design problems, e.g., optimum design of rotor machines [Romero and Silva, 2014], aerodynamics [Li et al., 2022b, Ghasemi and Elham, 2022], heat sinks [Yaji et al., 2015, Li et al., 2022a] in three dimensions.

Topology optimization was also applied to design problems of fluid diodes [Deng et al., 2010, Lin et al., 2015]. Lin et al. [2015] formulated the optimal design problem for fluid diodes as maximizing diodicity subject to a volume constraint with an additional penalty on flow through intermediate density. Sato et al. [2017a] proposed a bi-objective topology optimization method for fluid diodes and investigated Pareto-optimal solutions. In work done by Guo et al.

[2020], length scale constraints were introduced into two-dimensional topology optimization for the design of fluid diodes to reduce the difference between the performance of the two-dimensional model and the three-dimensional model obtained from the extrusion.

Although these previous studies showed that topology optimization can generate novel designs of fluid diodes with improved performances, one limitation in the design problems is the two-dimensional assumption of fluid flow. This assumption is valid as long as the effect of upper and lower boundaries is negligible. In practice, however, the height of the flow channel is not infinite, and the discrepancy between the performances of the two-dimensional model and the extruded three-dimensional model should not be neglected. Lin et al. [2015] reported that the diodicity values of the extruded three-dimensional model of fluid diodes measured by both numerical simulations and experiments were less than those of the two-dimensional model, and they concluded that topology optimization in three dimensions should be performed when the height of the channel cannot be neglected. Although Guo et al. [2020] showed that length scale control can mitigate the discrepancy, difficulty in estimating the final performance still remains. Thus, this work sheds a light on the topology optimization of fluid diodes in full-scale three-dimensional setting.

However, topology optimization of fluid diodes in three dimensions poses another challenging issue. The optimized structures may feature some “floating” islands disconnected to the external walls, which is not ideal in practice, therefore, cannot be manufactured. The Brinkman model allows the nucleation of such structures owing to the nature of its mathematical modelling. Therefore, we need to construct a design methodology for obtaining wall-connected structures in the fluid topology optimization.

Previous works have investigated geometric constraints for ensuring desired geometrical features of the structure, especially in the fields of solid mechanics and additive manufacturing. In topology optimization, the control of minimum and maximum length scales of structures has long been a typical issue of geometric constraints [Poulsen, 2003, Guest et al., 2004, Zhou et al., 2015, Guest, 2009, Li et al., 2023], and the concept of geometric constraints have been
extended to molding constraints [Xia et al., 2010, Allaire et al., 2016, Sato et al., 2017b] and enclosed voids restriction constraints [Liu et al., 2015, Li et al., 2016, Yamada and Noguchi, 2022]. In the work of Liu et al. [2015], a Virtual Temperature Method (VTM) was constructed, in which a virtual boundary value problem is introduced so that the enclosed void regions can appear with a high value of virtual temperature. The fictitious physical model to constrain closed cavities proposed by Yamada and Noguchi [2022] is similar to VTM in concept, but the governing equations for the fictitious physical field are different. Both of these methods successfully obtained optimal structures in which enclosed void regions do not exist, and all void regions have connections to the outer boundaries. We can impose a geometric constraint for ensuring wall-connected structures of the fluid diodes by considering the enclosed void regions in these methods as floating solid regions.

This paper proposes a topology optimization of fluid diodes considering wall-connected structures. Since excessive energy dissipation in the forward flow impairs the original function as a flow path, the optimal design problem of fluid diodes is formulated as maximizing the energy dissipation in the reverse flow subject to a maximum allowed energy dissipation in the forward flow. An additional geometric constraint based on the fictitious physical approach is formulated to avoid the nucleation of “floating” islands in the optimum design problem. To the authors’ best knowledge, very limited works have been proposed to guarantee wall-connected structures in the three-dimensional fluid topology optimization. In numerical implementation, to overcome the problems of numerical costs of three-dimensional fluid flow analysis, every implementation is performed in OpenFOAM, an open-source finite volume software which is widely used in computational fluid dynamics and topology optimization related to fluid problems [Yu et al., 2020, Pan et al., 2022, Galanos et al., 2022].

The remainder of this paper is organized as follows. In Section 2, the governing equations and density-based topology optimization formulation of incompressible steady-state flow are described, and the fictitious physical field and the geometric constraint are introduced. Then, the optimal design problem of fluid diodes considering wall-connected structures is defined. In Section 3,

we describe the implementation details and the optimization algorithm. In Section 4, we show a numerical verification of the fictitious physical model and design examples of fluid diodes obtained from the proposed method. Finally, conclusions are provided in Section 5.

## 2 Formulation

### 2.1 Governing equations

The governing equations for an incompressible, steady-state flow in the non-dimensional form are given as:

$$ 
\left\{ 
\begin{array}{ll}
\nabla \cdot (u u^{\top}) = \nabla \cdot \bar{\sigma}(u, p) + f & \text{in } \Omega \\
\nabla \cdot u = 0 & \text{in } \Omega \\
u = u_{\text{D}} & \text{on } \Gamma_{\text{D}} \\
\bar{\sigma}(u, p) \cdot n = 0 & \text{on } \Gamma_{\text{N}},
\end{array}
\right.
\quad (1) 
$$

where $n$ is the unit outward normal vector on the boundary, and $u, p, f$, and $\nabla$ denote the non-dimensional velocity, pressure, body force, and gradient operator, respectively. These non-dimensional quantities are defined as follows:

$$ u = \frac{u^*}{U}, \quad p = \frac{p^*}{\rho U^2}, \quad f = \frac{L}{U^2} f^*, \quad \nabla = L \nabla^*, \quad (2) $$

where $U$ and $L$ are the characteristic velocity and length, $\rho$ and $\nu$ are the fluid density and dynamic viscosity, and $*$ denotes the dimensional quantities. The fluid stress tensor $\bar{\sigma}(u, p)$ is defined as follows:

$$ \bar{\sigma}(u, p) := -p\mathbf{I} + \frac{1}{\text{Re}}(\nabla u + \nabla u^{\top}), \quad \text{Re} = \frac{UL}{\nu}. \quad (3) $$

where $\mathbf{I}$ is the identity tensor, and $\text{Re}$ is the Reynolds number. The boundary $\Gamma := \partial\Omega$ is composed of the Dirichlet boundary $\Gamma_{\text{D}}$, where a velocity profile is given, and the Neumann boundary $\Gamma_{\text{N}}$, where the boundary stress vector $\sigma := \bar{\sigma} \cdot n$ is specified as zero.

### 2.2 Formulation of the density-based topology optimization

In this research, the density-based topology optimization method with filtering and projection is used [Lazarov and Sigmund, 2011, Wang et al., 2011, Kawamoto et al., 2011]. A pseudo-density
extended to molding constraints [Xia et al., 2010, Allaire et al., 2016, Sato et al., 2017b] and enclosed voids restriction constraints [Liu et al., 2015, Li et al., 2016, Yamada and Noguchi, 2022]. In the work of Liu et al. [2015], a Virtual Temperature Method (VTM) was constructed, in which a virtual boundary value problem is introduced so that the enclosed void regions can appear with a high value of virtual temperature. The fictitious physical model to constrain closed cavities proposed by Yamada and Noguchi [2022] is similar to VTM in concept, but the governing equations for the fictitious physical field are different. Both of these methods successfully obtained optimal structures in which enclosed void regions do not exist, and all void regions have connections to the outer boundaries. We can impose a geometric constraint for ensuring wall-connected structures of the fluid diodes by considering the enclosed void regions in these methods as floating solid regions.

This paper proposes a topology optimization of fluid diodes considering wall-connected structures. Since excessive energy dissipation in the forward flow impairs the original function as a flow path, the optimal design problem of fluid diodes is formulated as maximizing the energy dissipation in the reverse flow subject to a maximum allowed energy dissipation in the forward flow. An additional geometric constraint based on the fictitious physical approach is formulated to avoid the nucleation of “floating” islands in the optimum design problem. To the authors’ best knowledge, very limited works have been proposed to guarantee wall-connected structures in the three-dimensional fluid topology optimization. In numerical implementation, to overcome the problems of numerical costs of three-dimensional fluid flow analysis, every implementation is performed in OpenFOAM, an open-source finite volume software which is widely used in computational fluid dynamics and topology optimization related to fluid problems [Yu et al., 2020, Pan et al., 2022, Galanos et al., 2022].

The remainder of this paper is organized as follows. In Section 2, the governing equations and density-based topology optimization formulation of incompressible steady-state flow are described, and the fictitious physical field and the geometric constraint are introduced. Then, the optimal design problem of fluid diodes considering wall-connected structures is defined. In Section 3,

we describe the implementation details and the optimization algorithm. In Section 4, we show a numerical verification of the fictitious physical model and design examples of fluid diodes obtained from the proposed method. Finally, conclusions are provided in Section 5.

# 2 Formulation

## 2.1 Governing equations

The governing equations for an incompressible, steady-state flow in the non-dimensional form are given as:

$$ 
\left\{ 
\begin{array}{ll} 
\nabla \cdot (u u^{\top}) = \nabla \cdot \bar{\sigma}(u, p) + f & \text{in } \Omega \\ 
\nabla \cdot u = 0 & \text{in } \Omega \\ 
u = u_D & \text{on } \Gamma_D \\ 
\bar{\sigma}(u, p) \cdot n = 0 & \text{on } \Gamma_N, 
\end{array} 
\right. 
\quad (1) 
$$

where $n$ is the unit outward normal vector on the boundary, and $u, p, f$, and $\nabla$ denote the non-dimensional velocity, pressure, body force, and gradient operator, respectively. These non-dimensional quantities are defined as follows:

$$ u = \frac{u^*}{U}, \quad p = \frac{p^*}{\rho U^2}, \quad f = \frac{L}{U^2} f^*, \quad \nabla = L \nabla^*, \quad (2) $$

where $U$ and $L$ are the characteristic velocity and length, $\rho$ and $\nu$ are the fluid density and dynamic viscosity, and $*$ denotes the dimensional quantities. The fluid stress tensor $\bar{\sigma}(u, p)$ is defined as follows:

$$ \bar{\sigma}(u, p) := -p\mathbf{I} + \frac{1}{\text{Re}}(\nabla u + \nabla u^{\top}), \quad \text{Re} = \frac{UL}{\nu}. \quad (3) $$

where $\mathbf{I}$ is the identity tensor, and $\text{Re}$ is the Reynolds number. The boundary $\Gamma := \partial\Omega$ is composed of the Dirichlet boundary $\Gamma_D$, where a velocity profile is given, and the Neumann boundary $\Gamma_N$, where the boundary stress vector $\sigma := \bar{\sigma} \cdot n$ is specified as zero.

## 2.2 Formulation of the density-based topology optimization

In this research, the density-based topology optimization method with filtering and projection is used [Lazarov and Sigmund, 2011, Wang et al., 2011, Kawamoto et al., 2011]. A pseudo-density
field $\gamma(\boldsymbol{x})$, which is a continuous function defined in $[0, 1]$, is introduced to describe the material distribution. To perform topology optimization, we consider the computational domain $\Omega$ which consists of the solid sub-domain $\Omega_s$ and the fluid sub-domain $\Omega_f$. The solid domain $\Omega_s$ is defined as the domain where $\gamma = 0$, and the fluid domain $\Omega_f$ is defined as the domain where $\gamma = 1$.

$$ \gamma(\boldsymbol{x}) = \begin{cases} 0 & \text{in } \Omega_s \\ 1 & \text{in } \Omega_f. \end{cases} \qquad (4) $$

The density filtering is employed to prevent numerical instabilities in topology optimization. The filtered density field $\tilde{\gamma}$ is obtained as the solution of the following Helmholtz-type partial differential equation.

$$ \begin{cases} -r^2 \nabla^2 \tilde{\gamma} + \tilde{\gamma} = \gamma & \text{in } \Omega \\ \nabla \tilde{\gamma} \cdot \boldsymbol{n} = 0 & \text{on } \Gamma, \end{cases} \qquad (5) $$

where $r$ denotes the filter radius. Finally, the projected density field $\hat{\gamma}$ is obtained by applying the smoothed Heaviside function $H(\tilde{\gamma}; \eta, \beta)$ as follows:

$$ \hat{\gamma} := H(\tilde{\gamma}; \eta, \beta) = \frac{\tanh(\beta\eta) + \tanh(\beta(\tilde{\gamma} - \eta))}{\tanh(\beta\eta) + \tanh(\beta(1 - \eta))}, \qquad (6) $$

where $\eta$ and $\beta$ denote the threshold and the sharpness parameters, respectively. Fig. 1 shows a plot of the smoothed Heaviside function for different $\beta$.

<figure><img src="image_1.png" /><figcaption>Fig. 1: Smoothed Heaviside function $H(\tilde{\gamma}; \eta, \beta)$, $\eta = 0.5$</figcaption></figure>

In order to make the governing equations to be applicable in the entire fixed design domain, a fictitious force $\boldsymbol{f} = -\alpha \boldsymbol{u}$ is introduced based on the Darcy law (also called the Brinkman model) [Borrvall and Petersson, 2003].

$$ \left\{ \begin{array}{ll} \nabla \cdot (\boldsymbol{u}\boldsymbol{u}^{\top}) = \nabla \cdot \bar{\boldsymbol{\sigma}}(\boldsymbol{u}, p) - \alpha(\hat{\gamma})\boldsymbol{u} & \text{in } \Omega \\ \nabla \cdot \boldsymbol{u} = 0 & \text{in } \Omega \\ \boldsymbol{u} = \boldsymbol{u}_{\mathrm{D}} & \text{on } \Gamma_{\mathrm{D}} \\ \bar{\boldsymbol{\sigma}}(\boldsymbol{u}, p) \cdot \boldsymbol{n} = \mathbf{0} & \text{on } \Gamma_{\mathrm{N}}, \end{array} \right. \qquad (7) $$

The solid domain is represented as a porous material, where the inverse permeability $\alpha(\hat{\gamma})$ takes the value much greater than 1. In this work, the function $\alpha(\hat{\gamma})$ is defined as follows:

$$ \alpha(\hat{\gamma}) = \alpha_{\max} (1 - \hat{\gamma}), \qquad (8) $$

where $\alpha_{\max}$ is the maximum value of the inverse permeability.

## 2.3 Geometric constraint for excluding “floating” islands

As has been discussed in Section 1, “floating” solid regions entirely surrounded by the fluid may be nucleated in the design problem of fluid diodes. In order to obtain wall-connected structures in the fluid topology optimization, we must eliminate such “floating” solid regions or make them connected to some fixed parts, e.g., no-slip wall boundaries. Therefore, we introduce a geometric constraint for excluding “floating” islands based on the fictitious physical model [Yamada and Noguchi, 2022].

The key point is to set up a boundary value problem similar to a heat conduction problem where the floating solid regions are detected as the hot portion in the solution. The governing equations for the fictitious temperature field $T$ are defined as follows:

$$ \left\{ \begin{array}{ll} \nabla \cdot (d(\hat{\gamma}) \nabla T) + (1 - \hat{\gamma} - T) = 0 & \text{in } \Omega \\ T = 0 & \text{on } \Gamma_T \qquad (9) \\ d \nabla T \cdot \boldsymbol{n} = 0 & \text{on } \Gamma_q, \end{array} \right. $$

where $d$ is the diffusion coefficient defined as:

$$ d(\hat{\gamma}) = d_{\max} + (d_{\min} - d_{\max}) \hat{\gamma}. \qquad (10) $$

The diffusion coefficient $d$ is made to take a large value $d_{\max}$ in the solid domain and a small value
field $\gamma(\boldsymbol{x})$, which is a continuous function defined in $[0, 1]$, is introduced to describe the material distribution. To perform topology optimization, we consider the computational domain $\Omega$ which consists of the solid sub-domain $\Omega_s$ and the fluid sub-domain $\Omega_f$. The solid domain $\Omega_s$ is defined as the domain where $\gamma = 0$, and the fluid domain $\Omega_f$ is defined as the domain where $\gamma = 1$.

$$ \gamma(\boldsymbol{x}) = \begin{cases} 0 & \text{in } \Omega_s \\ 1 & \text{in } \Omega_f. \end{cases} \qquad (4) $$

The density filtering is employed to prevent numerical instabilities in topology optimization. The filtered density field $\tilde{\gamma}$ is obtained as the solution of the following Helmholtz-type partial differential equation.

$$ \begin{cases} -r^2 \nabla^2 \tilde{\gamma} + \tilde{\gamma} = \gamma & \text{in } \Omega \\ \nabla \tilde{\gamma} \cdot \boldsymbol{n} = 0 & \text{on } \Gamma, \end{cases} \qquad (5) $$

where $r$ denotes the filter radius. Finally, the projected density field $\hat{\gamma}$ is obtained by applying the smoothed Heaviside function $H(\tilde{\gamma}; \eta, \beta)$ as follows:

$$ \hat{\gamma} := H(\tilde{\gamma}; \eta, \beta) = \frac{\tanh(\beta\eta) + \tanh(\beta(\tilde{\gamma} - \eta))}{\tanh(\beta\eta) + \tanh(\beta(1 - \eta))}, \qquad (6) $$

where $\eta$ and $\beta$ denote the threshold and the sharpness parameters, respectively. Fig. 1 shows a plot of the smoothed Heaviside function for different $\beta$.

<figure><img src="image_4.png" /><figcaption>Fig. 1: Smoothed Heaviside function $H(\tilde{\gamma}; \eta, \beta)$, $\eta = 0.5$</figcaption></figure>

In order to make the governing equations to be applicable in the entire fixed design domain, a fictitious force $\boldsymbol{f} = -\alpha \boldsymbol{u}$ is introduced based on the Darcy law (also called the Brinkman model) [Borrvall and Petersson, 2003].

$$ \left\{ \begin{array}{ll} \nabla \cdot (\boldsymbol{u}\boldsymbol{u}^{\top}) = \nabla \cdot \bar{\boldsymbol{\sigma}}(\boldsymbol{u}, p) - \alpha(\hat{\gamma})\boldsymbol{u} & \text{in } \Omega \\ \nabla \cdot \boldsymbol{u} = 0 & \text{in } \Omega \\ \boldsymbol{u} = \boldsymbol{u}_{\mathrm{D}} & \text{on } \Gamma_{\mathrm{D}} \\ \bar{\boldsymbol{\sigma}}(\boldsymbol{u}, p) \cdot \boldsymbol{n} = \mathbf{0} & \text{on } \Gamma_{\mathrm{N}}, \end{array} \right. \qquad (7) $$

The solid domain is represented as a porous material, where the inverse permeability $\alpha(\hat{\gamma})$ takes the value much greater than 1. In this work, the function $\alpha(\hat{\gamma})$ is defined as follows:

$$ \alpha(\hat{\gamma}) = \alpha_{\max}(1 - \hat{\gamma}), \qquad (8) $$

where $\alpha_{\max}$ is the maximum value of the inverse permeability.

## 2.3 Geometric constraint for excluding “floating” islands

As has been discussed In Section 1, “floating” solid regions entirely surrounded by the fluid may be nucleated in the design problem of fluid diodes. In order to obtain wall-connected structures in the fluid topology optimization, we must eliminate such “floating” solid regions or make them connected to some fixed parts, e.g., no-slip wall boundaries. Therefore, we introduce a geometric constraint for excluding “floating” islands based on the fictitious physical model [Yamada and Noguchi, 2022].

The key point is to set up a boundary value problem similar to a heat conduction problem where the floating solid regions are detected as the hot portion in the solution. The governing equations for the fictitious temperature field $T$ are defined as follows:

$$ \left\{ \begin{array}{ll} \nabla \cdot (d(\hat{\gamma}) \nabla T) + (1 - \hat{\gamma} - T) = 0 & \text{in } \Omega \\ T = 0 & \text{on } \Gamma_T \quad (9) \\ d \nabla T \cdot \boldsymbol{n} = 0 & \text{on } \Gamma_q, \end{array} \right. $$

where $d$ is the diffusion coefficient defined as:

$$ d(\hat{\gamma}) = d_{\max} + (d_{\min} - d_{\max}) \hat{\gamma}. \qquad (10) $$

The diffusion coefficient $d$ is made to take a large value $d_{\max}$ in the solid domain and a small value
$d_{min}$ in the fluid domain. There is a heat source of $O(1)$ in the solid domain and no heat source in the fluid domain. Also, the fictitious temperature value is set to zero on the Dirichlet boundary $\Gamma_T$ to which all the solid regions are expected to be connected.

Under this problem setup, in the solid regions connected to a Dirichlet boundary $\Gamma_T$, heat is transferred instantaneously by the very high diffusion coefficient $d_{max}$, and the fictitious temperature becomes nearly zero. But in the solid regions not connected to the Dirichlet boundary $\Gamma_T$, heat cannot be transferred anywhere since it is surrounded by the fluid region with very low conductivity $d_{min}$, and the fictitious temperature takes the value around one. In this way, “floating” islands can be detected as regions with fictitious temperature of $O(1)$.

In addition, to restrict these floating solid regions, the following cost function is introduced:

$$J_h = \int_{\Omega} (1 - \hat{\gamma}) T^2 d\Omega. \quad (11)$$

It should be noted that, as mentioned in the work of Yamada and Noguchi [2022], the choice of the cost function $J_h$ is not unique. This function (11) is expressed as the integral of the square of the temperature in the solid domain $\Omega_s$ and is expected to have a large value when there is a floating solid region. Then, setting a sufficiently small upper bound for this cost function can constrain the nucleation of those “floating” solid regions.

## 2.4 Optimum design of fluid diodes with the geometric constraint for excluding “floating” islands

Let $u_f, p_f$ denote the state variables representing the fluid velocity and pressure in the forward flow, respectively. Similarly, let $u_r, p_r$ denote those in the reverse flow. The viscous energy dissipation can be evaluated as the boundary integral of total pressure flux as follows:

$$J_f = - \int_{\Gamma} \left( p_f + \frac{u_f \cdot u_f}{2} \right) (u_f \cdot n) d\Gamma, \quad (12)$$

$$J_r = - \int_{\Gamma} \left( p_r + \frac{u_r \cdot u_r}{2} \right) (u_r \cdot n) d\Gamma, \quad (13)$$

where $J_f$ and $J_r$ denote the energy dissipation of the forward and the reverse flow, respectively. It should be noted that this expression in boundary integral form is identical to that in domain integral form used in the previous works [Lin et al., 2015, Sato et al., 2017a], see Appendix A.

In the most previous works, the diodicity, the ratio of the pressure drop in the reverse flow to that in the forward flow, is used as a performance indicator for fluid diodes.

$$Di = \frac{\Delta p_r}{\Delta p_f} \quad (14)$$

According to the previous study [Lin et al., 2015], the diodicity can also be defined using the energy dissipation as follows:

$$Di' = \frac{J_r}{J_f} = \frac{-\int_{\Gamma} \left( p_r + \frac{u_r \cdot u_r}{2} \right) (u_r \cdot n) d\Gamma}{-\int_{\Gamma} \left( p_f + \frac{u_f \cdot u_f}{2} \right) (u_f \cdot n) d\Gamma} \approx Di. \quad (15)$$

A large value of the diodicity $Di'$ indicates a significant difference between the energy lost in the forward and the reverse flows, quantifying the directionality of the flow.

However, as mentioned in [Sato et al., 2017a], a fluid diode with high diodicity is not always an appropriate design. It is desirable to avoid excessive energy dissipation in the forward flow to mitigate the driving force. From this standpoint, we formulate the optimal design problem of maximizing the energy dissipation in the reverse flow subject to a maximum allowed energy dissipation in the forward flow. Assuming that this constraint becomes active in the process of maximizing the energy dissipation in the reverse flow, this problem is equivalent to the problem for maximizing diodicity. In addition, to restrict the floating solid regions, a geometric constraint is introduced using the fictitious physical model, cf. (9).

Finally, the optimization problem is described as follows:

$$\underset{\gamma}{\text{minimize}} \quad J = -J_r \quad (16)$$

$$\text{subject to:} \quad G_1 = J_f - J_{f,\text{max}} \le 0 \quad (17)$$

$$G_2 = J_h - J_{h,\text{max}} \le 0 \quad (18)$$

$$0 \le \gamma \le 1$$

Governing equations (7) and (9),
$d_{min}$ in the fluid domain. There is a heat source of $O(1)$ in the solid domain and no heat source in the fluid domain. Also, the fictitious temperature value is set to zero on the Dirichlet boundary $\Gamma_T$ to which all the solid regions are expected to be connected.

Under this problem setup, in the solid regions connected to a Dirichlet boundary $\Gamma_T$, heat is transferred instantaneously by the very high diffusion coefficient $d_{max}$, and the fictitious temperature becomes nearly zero. But in the solid regions not connected to the Dirichlet boundary $\Gamma_T$, heat cannot be transferred anywhere since it is surrounded by the fluid region with very low conductivity $d_{min}$, and the fictitious temperature takes the value around one. In this way, “floating” islands can be detected as regions with fictitious temperature of $O(1)$.

In addition, to restrict these floating solid regions, the following cost function is introduced:

$$J_h = \int_{\Omega} (1 - \hat{\gamma}) T^2 d\Omega. \quad (11)$$

It should be noted that, as mentioned in the work of Yamada and Noguchi [2022], the choice of the cost function $J_h$ is not unique. This function (11) is expressed as the integral of the square of the temperature in the solid domain $\Omega_s$ and is expected to have a large value when there is a floating solid region. Then, setting a sufficiently small upper bound for this cost function can constrain the nucleation of those “floating” solid regions.

## 2.4 Optimum design of fluid diodes with the geometric constraint for excluding “floating” islands

Let $u_f, p_f$ denote the state variables representing the fluid velocity and pressure in the forward flow, respectively. Similarly, let $u_r, p_r$ denote those in the reverse flow. The viscous energy dissipation can be evaluated as the boundary integral of total pressure flux as follows:

$$J_f = - \int_{\Gamma} \left( p_f + \frac{u_f \cdot u_f}{2} \right) (u_f \cdot n) d\Gamma, \quad (12)$$

$$J_r = - \int_{\Gamma} \left( p_r + \frac{u_r \cdot u_r}{2} \right) (u_r \cdot n) d\Gamma, \quad (13)$$

where $J_f$ and $J_r$ denote the energy dissipation of the forward and the reverse flow, respectively. It should be noted that this expression in boundary integral form is identical to that in domain integral form used in the previous works [Lin et al., 2015, Sato et al., 2017a], see Appendix A.

In the most previous works, the diodicity, the ratio of the pressure drop in the reverse flow to that in the forward flow, is used as a performance indicator for fluid diodes.

$$Di = \frac{\Delta p_r}{\Delta p_f} \quad (14)$$

According to the previous study [Lin et al., 2015], the diodicity can also be defined using the energy dissipation as follows:

$$Di' = \frac{J_r}{J_f} = \frac{-\int_{\Gamma} \left( p_r + \frac{u_r \cdot u_r}{2} \right) (u_r \cdot n) d\Gamma}{-\int_{\Gamma} \left( p_f + \frac{u_f \cdot u_f}{2} \right) (u_f \cdot n) d\Gamma} \approx Di. \quad (15)$$

A large value of the diodicity $Di'$ indicates a significant difference between the energy lost in the forward and the reverse flows, quantifying the directionality of the flow.

However, as mentioned in [Sato et al., 2017a], a fluid diode with high diodicity is not always an appropriate design. It is desirable to avoid excessive energy dissipation in the forward flow to mitigate the driving force. From this standpoint, we formulate the optimal design problem of maximizing the energy dissipation in the reverse flow subject to a maximum allowed energy dissipation in the forward flow. Assuming that this constraint becomes active in the process of maximizing the energy dissipation in the reverse flow, this problem is equivalent to the problem for maximizing diodicity. In addition, to restrict the floating solid regions, a geometric constraint is introduced using the fictitious physical model, cf. (9).

Finally, the optimization problem is described as follows:

$$\underset{\gamma}{\text{minimize}} \quad J = -J_r \quad (16)$$

$$\text{subject to:} \quad G_1 = J_f - J_{f,max} \le 0 \quad (17)$$

$$G_2 = J_h - J_{h,max} \le 0 \quad (18)$$

$$0 \le \gamma \le 1$$

Governing equations (7) and (9),
where $J_{f,max}$ and $J_{h,max}$ are the upper bound values of $J_f$ and $J_h$, respectively. The boundary conditions for the forward flow $\mathbf{u}_f$, $p_f$, the reverse flow $\mathbf{u}_r$, $p_r$, and the fictitious temperature field $T$ are described as follows:

$$ \left\{ \begin{array}{ll} \mathbf{u}_f = \mathbf{u}_{\text{in,f}} & \text{on } \Gamma_{\text{in,f}} \\ \mathbf{u}_f = \mathbf{0} & \text{on } \Gamma_{\text{wall}} \\ \bar{\sigma}(\mathbf{u}_f, p_f) \cdot \mathbf{n} = \mathbf{0} & \text{on } \Gamma_{\text{out,f}} \\ \mathbf{u}_r = \mathbf{u}_{\text{in,r}} & \text{on } \Gamma_{\text{in,r}} \\ \mathbf{u}_r = \mathbf{0} & \text{on } \Gamma_{\text{wall}} \\ \bar{\sigma}(\mathbf{u}_r, p_r) \cdot \mathbf{n} = \mathbf{0} & \text{on } \Gamma_{\text{out,r}} \\ T = 0 & \text{on } \Gamma_T \\ d\nabla T \cdot \mathbf{n} = 0 & \text{on } \Gamma_q, \end{array} \right. \quad (19) $$

where $\Gamma_{\text{in,f}}$ and $\Gamma_{\text{out,f}}$ represent the inlet and the outlet for the forward flow, and $\Gamma_{\text{in,r}}$ and $\Gamma_{\text{out,r}}$ represent those for the reverse flow, respectively. The inlet of the forward flow $\Gamma_{\text{in,f}}$ corresponds to the outlet of the reverse flow $\Gamma_{\text{out,r}}$, and the outlet of the forward flow $\Gamma_{\text{out,f}}$ corresponds to the inlet of the reverse flow $\Gamma_{\text{in,r}}$. For the fictitious physical field, the Dirichlet boundary $\Gamma_T$ corresponds to the portion of the wall boundary $\Gamma_{\text{wall}}$ that is adjacent to the design domain, while the other boundaries $\Gamma \setminus \Gamma_T$ are considered as the Neumann boundary $\Gamma_q$.

The sensitivities of $J$, $G_1$ and $G_2$ with respect to the design variable $\gamma$ are derived based on the continuous adjoint approach. Adjoint equations should be solved in each optimization step, and the sensitivities are calculated using both state and adjoint variables. See Appendix B for details on the adjoint systems and the sensitivities.

# 3 Implementation details

## 3.1 Implementation in OpenFOAM

As mentioned above, topology optimization of fluid diode design requires forward and reverse flow simulations and the solutions of their adjoint equations. Furthermore, the computational costs of three-dimensional fluid analysis are considerably higher than that of two-dimensional analysis. To mitigate such computational burdens, efficient fluid analysis methods must be employed.

In this work, we employed an open-source finite volume software, OpenFOAM v5, to discretize the PDEs. A Semi-Implicit Method for Pressure Linked Equation (SIMPLE) algorithm was used to solve Navier-Stokes equations. Adjoint Navier-Stokes equations were also solved using a SIMPLE-type algorithm based on the work done by Othmer [2008]. The fictitious heat equation was solved based on the implementation of laplacianFoam, a default solver for heat equations in OpenFOAM.

## 3.2 Optimization algorithm

The optimization procedure is described in Algorithm 1. First, we prepare a mesh for computational domain and optimization parameters, such as upper bound values of the constraint functions. Then, we decompose the computational domain for parallel computing, declare unknowns, initialize the psuedo density field $\gamma$, and compute $\tilde{\gamma}$ from (5), and $\hat{\gamma}$ from (6). Material properties such as the inverse permeability $\alpha$ and the fictitious thermal conductivity $d$ are interpolated using the projected density field $\hat{\gamma}$ from (8) and (10).

Next, the optimization loop starts by solving governing equations of primal problems, i.e., Navier-Stokes equations for forward and reverse flows (7) and the fictitious heat equation (9). Then, objective function (16) and constraint functions (17) and (18) are evaluated. If the convergence criterion is satisfied, the optimization is terminated. Otherwise, the adjoint governing equations of forward and reverse flows and fictitious heat problems are solved. Sensitivities with respect to the projected density $\hat{\gamma}$ are computed using primal and adjoint variables. We apply the chain rule to obtain the sensitivities with respect to the density $\gamma$. The density field $\gamma$ is updated using the Method of Moving Asymptotes (MMA) [Svanberg, 1987]. Then, the filtered density $\tilde{\gamma}$ and the projected density $\hat{\gamma}$ are obtained from the filtering (5) and projection (6). Material properties $\alpha$ and $d$ are interpolated using the updated projected density $\hat{\gamma}$ from (8) and (10).

In this work, the following continuation method is employed for sharpness parameter $\beta$ of Heaviside function (6).

$$ \beta_t = \begin{cases} 8.0 & (1 \le t < 100) \\ \min(40, \beta_{t-1} + 0.1) & (t \ge 100), \end{cases} \quad (20) $$
where $J_{f,max}$ and $J_{h,max}$ are the upper bound values of $J_f$ and $J_h$, respectively. The boundary conditions for the forward flow $u_f$, $p_f$, the reverse flow $u_r$, $p_r$, and the fictitious temperature field $T$ are described as follows:

$$ \left\{ \begin{array}{ll} \mathbf{u}_f = \mathbf{u}_{\text{in,f}} & \text{on } \Gamma_{\text{in,f}} \\ \mathbf{u}_f = \mathbf{0} & \text{on } \Gamma_{\text{wall}} \\ \bar{\sigma}(\mathbf{u}_f, p_f) \cdot \mathbf{n} = \mathbf{0} & \text{on } \Gamma_{\text{out,f}} \\ \mathbf{u}_r = \mathbf{u}_{\text{in,r}} & \text{on } \Gamma_{\text{in,r}} \\ \mathbf{u}_r = \mathbf{0} & \text{on } \Gamma_{\text{wall}} \\ \bar{\sigma}(\mathbf{u}_r, p_r) \cdot \mathbf{n} = \mathbf{0} & \text{on } \Gamma_{\text{out,r}} \\ T = 0 & \text{on } \Gamma_T \\ d\nabla T \cdot \mathbf{n} = 0 & \text{on } \Gamma_q, \end{array} \right. \quad (19) $$

where $\Gamma_{\text{in,f}}$ and $\Gamma_{\text{out,f}}$ represent the inlet and the outlet for the forward flow, and $\Gamma_{\text{in,r}}$ and $\Gamma_{\text{out,r}}$ represent those for the reverse flow, respectively. The inlet of the forward flow $\Gamma_{\text{in,f}}$ corresponds to the outlet of the reverse flow $\Gamma_{\text{out,r}}$, and the outlet of the forward flow $\Gamma_{\text{out,f}}$ corresponds to the inlet of the reverse flow $\Gamma_{\text{in,r}}$. For the fictitious physical field, the Dirichlet boundary $\Gamma_T$ corresponds to the portion of the wall boundary $\Gamma_{\text{wall}}$ that is adjacent to the design domain, while the other boundaries $\Gamma \setminus \Gamma_T$ are considered as the Neumann boundary $\Gamma_q$.

The sensitivities of $J$, $G_1$ and $G_2$ with respect to the design variable $\gamma$ are derived based on the continuous adjoint approach. Adjoint equations should be solved in each optimization step, and the sensitivities are calculated using both state and adjoint variables. See Appendix B for details on the adjoint systems and the sensitivities.

# 3 Implementation details

## 3.1 Implementation in OpenFOAM

As mentioned above, topology optimization of fluid diode design requires forward and reverse flow simulations and the solutions of their adjoint equations. Furthermore, the computational costs of three-dimensional fluid analysis are considerably higher than that of two-dimensional analysis. To mitigate such computational burdens, efficient fluid analysis methods must be employed.

In this work, we employed an open-source finite volume software, OpenFOAM v5, to discretize the PDEs. A Semi-Implicit Method for Pressure Linked Equation (SIMPLE) algorithm was used to solve Navier-Stokes equations. Adjoint Navier-Stokes equations were also solved using a SIMPLE-type algorithm based on the work done by Othmer [2008]. The fictitious heat equation was solved based on the implementation of laplacianFoam, a default solver for heat equations in OpenFOAM.

## 3.2 Optimization algorithm

The optimization procedure is described in Algorithm 1. First, we prepare a mesh for computational domain and optimization parameters, such as upper bound values of the constraint functions. Then, we decompose the computational domain for parallel computing, declare unknowns, initialize the psuedo density field $\gamma$, and compute $\tilde{\gamma}$ from (5), and $\hat{\gamma}$ from (6). Material properties such as the inverse permeability $\alpha$ and the fictitious thermal conductivity $d$ are interpolated using the projected density field $\hat{\gamma}$ from (8) and (10).

Next, the optimization loop starts by solving governing equations of primal problems, i.e., Navier-Stokes equations for forward and reverse flows (7) and the fictitious heat equation (9). Then, objective function (16) and constraint functions (17) and (18) are evaluated. If the convergence criterion is satisfied, the optimization is terminated. Otherwise, the adjoint governing equations of forward and reverse flows and fictitious heat problems are solved. Sensitivities with respect to the projected density $\hat{\gamma}$ are computed using primal and adjoint variables. We apply the chain rule to obtain the sensitivities with respect to the density $\gamma$. The density field $\gamma$ is updated using the Method of Moving Asymptotes (MMA) [Svanberg, 1987]. Then, the filtered density $\tilde{\gamma}$ and the projected density $\hat{\gamma}$ are obtained from the filtering (5) and projection (6). Material properties $\alpha$ and $d$ are interpolated using the updated projected density $\hat{\gamma}$ from (8) and (10).

In this work, the following continuation method is employed for sharpness parameter $\beta$ of Heaviside function (6).

$$ \beta_t = \begin{cases} 8.0 & (1 \le t < 100) \\ \min(40, \beta_{t-1} + 0.1) & (t \ge 100), \end{cases} \quad (20) $$
where the subscript *t* denotes the iteration number. Also, the threshold parameter *η* of Heaviside function (6) is chosen as 0.5.

The above procedure is repeated until the following criterion is satisfied: after the continuation method is finished, the relative difference between the objective function value of current and previous time steps becomes smaller than *ε*<sub>c</sub> while the constraints are satisfied.

Numerical computations in this paper are performed on a Linux computer with 20 Intel Xeon Gold 6242R CPU cores (3.10GHz).

# 4 Numerical examples

## 4.1 Validation of the fictitious physical model

In this section, we performed two-dimensional numerical analysis on the proposed fictitious physical model. The target geometry is shown in Fig. 2, whose boundary is set as Dirichlet boundary $Γ_T$ where $T = 0$. The size of the domain is set as $1.5 \times 2.0$.

<figure><img src="image_2.png" /><figcaption>**Fig. 2:** Target geometry for validation of the fictitious physical field. Black region represents solid $γ = 0$ and white region represents fluid $γ = 1$. Every boundary is considered as $Γ = Γ_T$ where $T = 0$.</figcaption></figure>

Fig. 3 shows fictitious temperature fields for different sets of maximum and minimum diffusion coefficients, and Fig. 4 shows the corresponding sensitivity distributions of $J_h$. The aim of the fictitious physical model is to detect the “floating” islands as regions with high value of $T$. The solid domain connected to the outer boundaries and the fluid domain should have a value of $T$ close to 0, so that only the structures without wall

connectivity will be penalized by the geometric constraint. Taking this into account, the maximum value of diffusion coefficient $d_{max}$ should be sufficiently large and the minimum value of diffusion coefficient $d_{min}$ should be sufficiently small.

However, in the case with $d_{max} = 1 \times 10^3$ such as Figs. 3o and t, the value of $T$ in the “floating” islands becomes smaller compared to the result with $d_{max} = 1 \times 10^2$, because of the excessively large maximum value of diffusion coefficient. Another characteristic of the fictitious physical field is its dependence on the size of the floating solid regions. In the fictitious physical model proposed in [Yamada and Noguchi, 2022], it was reported that the void region takes $T = 1$ regardless of its size. Although this feature is desirable in detection, our formulation may have an advantage in the context of capturing the boundaries of “floating” islands, which is discussed in Appendix C. From the above observations, it can be concluded that the proposed fictitious physical model is effective in detecting “floating” islands, and that the maximum value of diffusion coefficient should be $1 \times 10^1 \le d_{max} \le 1 \times 10^2$ and the minimum value of diffusion coefficient should be $d_{min} \le 1 \times 10^{-5}$.

Next, to verify the effectiveness of the cost function (11), we performed a non-constrained optimization problem of minimizing $J_h$. For this numerical example, the steepest decent method with a constant step size was employed as the optimizer, and the sensitivity of $J_h$ was normalized by its absolute maximum value at each time step. Here, we chose the diffusion coefficients as $(d_{max}, d_{min}) = (1 \times 10^1, 1 \times 10^{-5})$. Figs. 5 and 6 show the history of density distribution and the history of evaluation function $J_h$, respectively. In the first 20 iterations, most solid regions either vanished or were connected to the outer boundary $Γ_T$, then the solid regions gradually decreased, and after approximately 80 iterations, the solid regions vanished completely. The cost function $J_h$ also decreased monotonically during the optimization. Consequently, we confirm that the proposed fictitious physical model and the cost function $J_h$ can penalize the “floating” islands. It is worth noting that those $d_{min}$ and $d_{max}$ should be chosen carefully depending on the problem-at-hand, because they affect how the floating structures will be modified by the geometric constraint, i.e.,
where the subscript *t* denotes the iteration number. Also, the threshold parameter *η* of Heaviside function (6) is chosen as 0.5.

The above procedure is repeated until the following criterion is satisfied: after the continuation method is finished, the relative difference between the objective function value of current and previous time steps becomes smaller than *ε*<sub>c</sub> while the constraints are satisfied.

Numerical computations in this paper are performed on a Linux computer with 20 Intel Xeon Gold 6242R CPU cores (3.10GHz).

# 4 Numerical examples

## 4.1 Validation of the fictitious physical model

In this section, we performed two-dimensional numerical analysis on the proposed fictitious physical model. The target geometry is shown in Fig. 2, whose boundary is set as Dirichlet boundary $Γ_T$ where $T = 0$. The size of the domain is set as $1.5 \times 2.0$.

<figure><img src="image_5.png" /><figcaption>**Fig. 2:** Target geometry for validation of the fictitious physical field. Black region represents solid $γ = 0$ and white region represents fluid $γ = 1$. Every boundary is considered as $Γ = Γ_T$ where $T = 0$.</figcaption></figure>

Fig. 3 shows fictitious temperature fields for different sets of maximum and minimum diffusion coefficients, and Fig. 4 shows the corresponding sensitivity distributions of $J_h$. The aim of the fictitious physical model is to detect the “floating” islands as regions with high value of $T$. The solid domain connected to the outer boundaries and the fluid domain should have a value of $T$ close to 0, so that only the structures without wall

connectivity will be penalized by the geometric constraint. Taking this into account, the maximum value of diffusion coefficient $d_{max}$ should be sufficiently large and the minimum value of diffusion coefficient $d_{min}$ should be sufficiently small.

However, in the case with $d_{max} = 1 \times 10^3$ such as Figs. 3o and t, the value of $T$ in the “floating” islands becomes smaller compared to the result with $d_{max} = 1 \times 10^2$, because of the excessively large maximum value of diffusion coefficient. Another characteristic of the fictitious physical field is its dependence on the size of the floating solid regions. In the fictitious physical model proposed in [Yamada and Noguchi, 2022], it was reported that the void region takes $T = 1$ regardless of its size. Although this feature is desirable in detection, our formulation may have an advantage in the context of capturing the boundaries of “floating” islands, which is discussed in Appendix C. From the above observations, it can be concluded that the proposed fictitious physical model is effective in detecting “floating” islands, and that the maximum value of diffusion coefficient should be $1 \times 10^1 \le d_{max} \le 1 \times 10^2$ and the minimum value of diffusion coefficient should be $d_{min} \le 1 \times 10^{-5}$.

Next, to verify the effectiveness of the cost function (11), we performed a non-constrained optimization problem of minimizing $J_h$. For this numerical example, the steepest decent method with a constant step size was employed as the optimizer, and the sensitivity of $J_h$ was normalized by its absolute maximum value at each time step. Here, we chose the diffusion coefficients as $(d_{max}, d_{min}) = (1 \times 10^1, 1 \times 10^{-5})$. Figs. 5 and 6 show the history of density distribution and the history of evaluation function $J_h$, respectively. In the first 20 iterations, most solid regions either vanished or were connected to the outer boundary $Γ_T$, then the solid regions gradually decreased, and after approximately 80 iterations, the solid regions vanished completely. The cost function $J_h$ also decreased monotonically during the optimization. Consequently, we confirm that the proposed fictitious physical model and the cost function $J_h$ can penalize the “floating” islands. It is worth noting that those $d_{min}$ and $d_{max}$ should be chosen carefully depending on the problem-at-hand, because they affect how the floating structures will be modified by the geometric constraint, i.e.,
<table>
  <thead>
    <tr>
      <th>Step</th>
      <th>Description</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>358</td>
      <td><b>Algorithm 1:</b> Optimization problem of fluid diodes with the geometric constraint for excluding "floating" islands</td>
    </tr>
    <tr>
      <td>359</td>
      <td>input: mesh, optimization parameters</td>
    </tr>
    <tr>
      <td>360</td>
      <td>initialization:</td>
    </tr>
    <tr>
      <td>361</td>
      <td>decompose the computational domain for parallel computing</td>
    </tr>
    <tr>
      <td>362</td>
      <td>declare unknowns</td>
    </tr>
    <tr>
      <td>363</td>
      <td>initialize the density fields &gamma;, &gamma;&#771;, and &gamma;&#770;</td>
    </tr>
    <tr>
      <td>364</td>
      <td>interpolate material properties &alpha; and d</td>
    </tr>
    <tr>
      <td>365</td>
      <td>begin</td>
    </tr>
    <tr>
      <td>366</td>
      <td>for it = 0; it &lt; Maxloop; it = it + 1 do</td>
    </tr>
    <tr>
      <td>367</td>
      <td>solve Navier-Stokes equations for forward flow using the SIMPLE algorithm</td>
    </tr>
    <tr>
      <td>368</td>
      <td>solve Navier-Stokes equations for reverse flow using the SIMPLE algorithm</td>
    </tr>
    <tr>
      <td>369</td>
      <td>solve fictitious heat conduction equation</td>
    </tr>
    <tr>
      <td>370</td>
      <td>evaluate objective and constraint functions</td>
    </tr>
    <tr>
      <td>371</td>
      <td>if ||<sup>J<sub>it</sub> - J<sub>it-1</sub></sup>&frasl;<sub>J<sub>it</sub></sub>|| &lt; &epsilon;<sub>c</sub>, G<sub>1</sub> &lt; 0, G<sub>2</sub> &lt; 0 then</td>
    </tr>
    <tr>
      <td>372</td>
      <td>break;</td>
    </tr>
    <tr>
      <td>373</td>
      <td>else</td>
    </tr>
    <tr>
      <td>374</td>
      <td>solve adjoint Navier-Stokes equations for forward flow using the SIMPLE algorithm</td>
    </tr>
    <tr>
      <td>375</td>
      <td>solve adjoint Navier-Stokes equations for reverse flow using the SIMPLE algorithm</td>
    </tr>
    <tr>
      <td>376</td>
      <td>solve adjoint fictitious heat conduction equation</td>
    </tr>
    <tr>
      <td>377</td>
      <td>calculate sensitivities</td>
    </tr>
    <tr>
      <td>378</td>
      <td>update the density field &gamma; using MMA</td>
    </tr>
    <tr>
      <td>379</td>
      <td>solve PDE-filter to get &gamma;&#771;</td>
    </tr>
    <tr>
      <td>380</td>
      <td>adopt Heaviside projection to get &gamma;&#770;</td>
    </tr>
    <tr>
      <td>381</td>
      <td>interpolate material properties &alpha; and d</td>
    </tr>
    <tr>
      <td>382</td>
      <td>end</td>
    </tr>
    <tr>
      <td>383</td>
      <td>end</td>
    </tr>
    <tr>
      <td>384</td>
      <td>end</td>
    </tr>
    <tr>
      <td>385</td>
      <td>end</td>
    </tr>
    <tr>
      <td>386</td>
      <td>end</td>
    </tr>
    <tr>
      <td>387</td>
      <td>end</td>
    </tr>
    <tr>
      <td>388</td>
      <td></td>
    </tr>
    <tr>
      <td>389</td>
      <td></td>
    </tr>
    <tr>
      <td>390</td>
      <td>whether the "floating" islands will disappear or be connected to the outer boundaries.</td>
    </tr>
    <tr>
      <td>391</td>
      <td></td>
    </tr>
    <tr>
      <td>392</td>
      <td></td>
    </tr>
    <tr>
      <td>393</td>
      <td><b>4.2 Three-dimensional design of fluid diodes with symmetric condition</b></td>
    </tr>
    <tr>
      <td>394</td>
      <td></td>
    </tr>
    <tr>
      <td>395</td>
      <td></td>
    </tr>
    <tr>
      <td>396</td>
      <td></td>
    </tr>
    <tr>
      <td>397</td>
      <td>In this section, optimizations of fluid diodes with and without the geometric constraint (18) were performed. Fig. 7a shows the overall view of the channel to be considered in this numerical example. Fig. 7b shows the computational domain and the boundary conditions for the fluid problem. To avoid excessive computational costs, we handle a quarter model of the entire domain and assume symmetry of the structure in this case. The computational domain was discretized with 906, 250 structured hexahedral cells. For parallel computing, the computational domain was decomposed</td>
    </tr>
    <tr>
      <td>398</td>
      <td>into 12 subdomains. In this work, the characteristic length L and the characteristic velocity U represent the inlet width and the average velocity at the inlet, respectively, and an incompressible flow with Reynolds number Re of 100 is assumed. Before optimization, the energy dissipation of the channel was evaluated and found to be approximately 0.25 assuming pure fluid phase in the computational domain. Based on the simulation result, the upper bound value of the energy dissipation in the forward flow J<sub>f,max</sub> was set to 1.0. The maximum value of inverse permeability &alpha;<sub>max</sub> was set to 1.0 &times; 10<sup>3</sup>. For parameters of density filtering, the filter radius r was set to 1.7 &times; 10<sup>-2</sup>, which corresponds to the length of 3 cells. For convergence criteria, &epsilon;<sub>c</sub> = 1 &times; 10<sup>-5</sup> was used. The initial guess was given as an uniform distribution &gamma; = 0.5 in the design domain.</td>
    </tr>
    <tr>
      <td>399</td>
      <td></td>
    </tr>
    <tr>
      <td>400</td>
      <td></td>
    </tr>
    <tr>
      <td>401</td>
      <td></td>
    </tr>
    <tr>
      <td>402</td>
      <td></td>
    </tr>
    <tr>
      <td>403</td>
      <td></td>
    </tr>
    <tr>
      <td>404</td>
      <td></td>
    </tr>
    <tr>
      <td>405</td>
      <td></td>
    </tr>
    <tr>
      <td>406</td>
      <td></td>
    </tr>
    <tr>
      <td>407</td>
      <td></td>
    </tr>
    <tr>
      <td>408</td>
      <td></td>
    </tr>
  </tbody>
</table>

whether the “floating” islands will disappear or be
connected to the outer boundaries.

## 4.2 Three-dimensional design of fluid diodes with symmetric condition

In this section, optimizations of fluid diodes with and without the geometric constraint (18) were performed. Fig. 7a shows the overall view of the channel to be considered in this numerical example. Fig. 7b shows the computational domain and the boundary conditions for the fluid problem. To avoid excessive computational costs, we handle a quarter model of the entire domain and assume symmetry of the structure in this case. The computational domain was discretized with 906, 250 structured hexahedral cells. For parallel computing, the computational domain was decomposed

into 12 subdomains. In this work, the character-
istic length *L* and the characteristic velocity *U*
represent the inlet width and the average velocity
at the inlet, respectively, and an incompressible
flow with Reynolds number *Re* of 100 is assumed.
Before optimization, the energy dissipation of the
channel was evaluated and found to be approx-
imately 0.25 assuming pure fluid phase in the
computational domain. Based on the simulation
result, the upper bound value of the energy dissi-
pation in the forward flow *J*<sub>*f,max*</sub> was set to 1.0.
The maximum value of inverse permeability *α*<sub>*max*</sub>
was set to 1.0 × 10<sup>3</sup>. For parameters of density fil-
tering, the filter radius *r* was set to 1.7 × 10<sup>-2</sup>,
which corresponds to the length of 3 cells. For
convergence criteria, *ε*<sub>*c*</sub> = 1 × 10<sup>-5</sup> was used. The
initial guess was given as an uniform distribution
γ = 0.5 in the design domain.
<table><caption>**Algorithm 1:** Optimization problem of fluid diodes with the geometric constraint for excluding "floating" islands</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>input:</td><td>mesh, optimization parameters</td></tr><tr><td>initialization:</td><td>decompose the computational domain for parallel computing<br>declare unknowns<br>initialize the density fields γ, γ̃, and ŷ<br>interpolate material properties α and d</td></tr><tr><td>begin</td><td></td></tr><tr><td>for it = 0; it &lt; Maxloop; it = it + 1 do</td><td>solve Navier-Stokes equations for forward flow using the SIMPLE algorithm<br>solve Navier-Stokes equations for reverse flow using the SIMPLE algorithm<br>solve fictitious heat conduction equation<br>evaluate objective and constraint functions</td></tr><tr><td>if</td><td>||<div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">J<sub>it</sub> - J<sub>it-1</sub></div><div style="padding-top: 1px;">J<sub>it</sub></div></div>|| &lt; ε<sub>c</sub>, G<sub>1</sub> &lt; 0, G<sub>2</sub> &lt; 0 then</td></tr><tr><td>break;</td><td></td></tr><tr><td>else</td><td>solve adjoint Navier-Stokes equations for forward flow using the SIMPLE algorithm<br>solve adjoint Navier-Stokes equations for reverse flow using the SIMPLE algorithm<br>solve adjoint fictitious heat conduction equation<br>calculate sensitivities<br>update the density field γ using MMA<br>solve PDE-filter to get γ̃<br>adopt Heaviside projection to get ŷ<br>interpolate material properties α and d</td></tr><tr><td>end</td><td></td></tr><tr><td>end</td><td></td></tr><tr><td>end</td><td></td></tr></tbody></table>

whether the "floating" islands will disappear or be connected to the outer boundaries.

## 4.2 Three-dimensional design of fluid diodes with symmetric condition

In this section, optimizations of fluid diodes with and without the geometric constraint (18) were performed. Fig. 7a shows the overall view of the channel to be considered in this numerical example. Fig. 7b shows the computational domain and the boundary conditions for the fluid problem. To avoid excessive computational costs, we handle a quarter model of the entire domain and assume symmetry of the structure in this case. The computational domain was discretized with 906, 250 structured hexahedral cells. For parallel computing, the computational domain was decomposed

into 12 subdomains. In this work, the characteristic length L and the characteristic velocity U represent the inlet width and the average velocity at the inlet, respectively, and an incompressible flow with Reynolds number Re of 100 is assumed. Before optimization, the energy dissipation of the channel was evaluated and found to be approximately 0.25 assuming pure fluid phase in the computational domain. Based on the simulation result, the upper bound value of the energy dissipation in the forward flow $J_{f,max}$ was set to 1.0. The maximum value of inverse permeability $α_{max}$ was set to $1.0 \times 10^3$. For parameters of density filtering, the filter radius $r$ was set to $1.7 \times 10^{-2}$, which corresponds to the length of 3 cells. For convergence criteria, $ε_c = 1 \times 10^{-5}$ was used. The initial guess was given as an uniform distribution $γ = 0.5$ in the design domain.
<figure><img src="image_3.png" /><figcaption>Fig. 3: Fictitious physical field $T$ for various diffusion coefficients ($d_{\text{max}}$, $d_{\text{min}}$)</figcaption></figure>

Boundary conditions for the fictitious physical field are shown in Fig. 8, where the wall adjacent to the design domain is considered as the Dirichlet boundary $\Gamma_T$. The parameters for the fictitious physical field were set as $d_{\text{max}} = 1 \times 10^2$ and $d_{\text{min}} = 1 \times 10^{-5}$ based on the result of the previous Section 4.1. Although the geometric constraint (18) was not imposed in Section 4.2.1, the function value (11) was evaluated for the comparative studies in later sections.

### 4.2.1 Optimization results without the geometric constraint

Fig. 9 shows the obtained configuration of a fluid diode plotted by the isovolume of $\hat{\gamma} \in [0, 0.5]$,

where Figs. 9a and b show the quarter and half model of the entire design, and Figs. 9c and d show streamlines of the forward and reverse flows, all respectively. Fig. 10 shows the optimization history of $J_f$, $J_r$, and $J_h$. The convergence criteria was satisfied in 458 optimization iterations. From Figs. 9c and d, it can be seen that the flow has generally straight paths from the inlet to the outlet in the forward flow, while the flow is bent by inertial forces and the structure, containing many vortices. The solid domains have holes, and exhibit a wall-like structure that separates the design domain into the main channel and the outside. These holes and outer spaces seem to be effectively utilized for changing the flow directions

409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
<figure><img src="image_6.png" /><figcaption>Fig. 3: Fictitious physical field $T$ for various diffusion coefficients ($d_{\text{max}}$, $d_{\text{min}}$)</figcaption></figure>

Boundary conditions for the fictitious physical field are shown in Fig. 8, where the wall adjacent to the design domain is considered as the Dirichlet boundary $\Gamma_T$. The parameters for the fictitious physical field were set as $d_{\text{max}} = 1 \times 10^2$ and $d_{\text{min}} = 1 \times 10^{-5}$ based on the result of the previous Section 4.1. Although the geometric constraint (18) was not imposed in Section 4.2.1, the function value (11) was evaluated for the comparative studies in later sections.

### 4.2.1 Optimization results without the geometric constraint

Fig. 9 shows the obtained configuration of a fluid diode plotted by the isovolume of $\hat{\gamma} \in [0, 0.5]$,

where Figs. 9a and b show the quarter and half model of the entire design, and Figs. 9c and d show streamlines of the forward and reverse flows, all respectively. Fig. 10 shows the optimization history of $J_f$, $J_r$, and $J_h$. The convergence criteria was satisfied in 458 optimization iterations. From Figs. 9c and d, it can be seen that the flow has generally straight paths from the inlet to the outlet in the forward flow, while the flow is bent by inertial forces and the structure, containing many vortices. The solid domains have holes, and exhibit a wall-like structure that separates the design domain into the main channel and the outside. These holes and outer spaces seem to be effectively utilized for changing the flow directions

409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
<figure><img src="image_4.png" /><figcaption>**Fig. 4:** Sensitivity distributions of the geometric constraint $J_h$ for various diffusion coefficients ($d_{max}$, $d_{min}$)</figcaption></figure>

in the reverse flow, maximizing the energy dissipation. The energy dissipation in the reverse flow $J_r$ increased during the optimization and reached approximately 1.52 at the end of the optimization process. Since the upper bound constraint for the energy dissipation of forward flow $J_f \le 1.0$ was satisfied with $J_f \approx 1.0$, the diodicity of the obtained structure was $Di \approx 1.52$.

Although the optimized result shows a satisfactory diodicity, the structure features floating solid islands, namely, the solid domains without wall connectivity. To manufacture the obtained structure, the design must be modified so that every solid domain can be supported by outer

boundaries. One possibility is to remove all “floating” islands, but it is clear that the “floating” regions act an essential role in the enhancement of diodicity in this example. However, it is also challenging for designers to modify it to a freestanding structure while maintaining diodicity, e.g., by adding supports. Therefore, we need to employ a geometric constraint to ensure wall-connected structures for the three-dimensional design of fluid diodes.

The history of function $J_h$ is plotted in Fig. 10, and the fictitious physical field over solid regions of the obtained structure is depicted in Fig. 11. It is seen that “floating” islands are detected as
<figure><img src="image_7.png" /><figcaption>**Fig. 4:** Sensitivity distributions of the geometric constraint $J_h$ for various diffusion coefficients ($d_{max}$, $d_{min}$)</figcaption></figure>

in the reverse flow, maximizing the energy dissipation. The energy dissipation in the reverse flow $J_r$ increased during the optimization and reached approximately 1.52 at the end of the optimization process. Since the upper bound constraint for the energy dissipation of forward flow $J_f \le 1.0$ was satisfied with $J_f \approx 1.0$, the diodicity of the obtained structure was $Di \approx 1.52$.

Although the optimized result shows a satisfactory diodicity, the structure features floating solid islands, namely, the solid domains without wall connectivity. To manufacture the obtained structure, the design must be modified so that every solid domain can be supported by outer

boundaries. One possibility is to remove all “floating” islands, but it is clear that the “floating” regions act an essential role in the enhancement of diodicity in this example. However, it is also challenging for designers to modify it to a freestanding structure while maintaining diodicity, e.g., by adding supports. Therefore, we need to employ a geometric constraint to ensure wall-connected structures for the three-dimensional design of fluid diodes.

The history of function $J_h$ is plotted in Fig. 10, and the fictitious physical field over solid regions of the obtained structure is depicted in Fig. 11. It is seen that “floating” islands are detected as
<figure><img src="image_5.png" /><figcaption>**Fig. 5:** Density distribution during minimizing $J_h$ with ($d_{max}$, $d_{min}$) = (1 × 10¹, 1 × 10⁻⁵)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 6:** Optimization history of minimizing $J_h$</figcaption></figure>

regions with high value of fictitious physical field $T$. The value of $J_h$ for the final configuration was $J_h \approx 3.9 \times 10^{-2}$. We call this value as $J_{h,ref}$, and we impose a sufficiently smaller upper bound value $J_{h,max}$ than $J_{h,ref}$ to constrain “floating” islands in the following section.

### 4.2.2 Optimization results with the geometric constraint

The problem settings and parameters were the same as in Section 4.2.1 except for introducing the geometric constraint (18) into the optimization problem. For the upper bound value $J_{h,max}$, 1 × 10⁻⁵ was used.

Fig. 12 shows the obtained configuration of a fluid diode with the geometric constraint plotted by the isovolume of $\hat{\gamma} \in [0, 0.5]$, where Figs. 12a and b show the quarter and half model of the

<figure><img src="image_7.png" /><figcaption>**Fig. 7:** Computational domain and boundary conditions for the fluid flow</figcaption></figure>

entire design, and Figs. 12c and d show streamlines of the forward and reverse flows, respectively. Fig. 13 shows the optimization history of $J_f$, $J_r$, and $J_h$. In this case, the convergence criteria was satisfied in 420 optimization iterations. From Figs. 12c and d, the flow in the obtained structure has features similar to the solution in the previous example, i.e., the fluid flows straight in the forward direction, while the many vortices exist in the reverse flow. Although the number of

511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
<figure><img src="image_8.png" /><figcaption>**Fig. 5:** Density distribution during minimizing $J_h$ with ($d_{max}$, $d_{min}$) = (1 × 10¹, 1 × 10⁻⁵)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 6:** Optimization history of minimizing $J_h$</figcaption></figure>

regions with high value of fictitious physical field $T$. The value of $J_h$ for the final configuration was $J_h \approx 3.9 \times 10^{-2}$. We call this value as $J_{h,ref}$, and we impose a sufficiently smaller upper bound value $J_{h,max}$ than $J_{h,ref}$ to constrain “floating” islands in the following section.

### 4.2.2 Optimization results with the geometric constraint

The problem settings and parameters were the same as in Section 4.2.1 except for introducing the geometric constraint (18) into the optimization problem. For the upper bound value $J_{h,max}$, 1 × 10⁻⁵ was used.

Fig. 12 shows the obtained configuration of a fluid diode with the geometric constraint plotted by the isovolume of $\hat{\gamma} \in [0, 0.5]$, where Figs. 12a and b show the quarter and half model of the

<figure><img src="image_10.png" /><figcaption>**Fig. 7:** Computational domain and boundary conditions for the fluid flow</figcaption></figure>

entire design, and Figs. 12c and d show streamlines of the forward and reverse flows, respectively. Fig. 13 shows the optimization history of $J_f$, $J_r$, and $J_h$. In this case, the convergence criteria was satisfied in 420 optimization iterations. From Figs. 12c and d, the flow in the obtained structure has features similar to the solution in the previous example, i.e., the fluid flows straight in the forward direction, while the many vortices exist in the reverse flow. Although the number of

511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
<figure><img src="image_8.png" /><figcaption>**Fig. 8:** Boundary conditions for the fictitious physical field</figcaption></figure>

holes appears to be less than that of the structure in Fig. 9, the resulting structure still has a unique topology. The energy dissipation in the reverse flow $J_r$ reached approximately 1.63, and the diodicity was $Di \approx 1.63$.

The difference between the obtained structures of Figs. 9 and 12 is the existence of floating solid regions. In the solution of Fig. 12, every solid part connects to the outer boundaries, and therefore the structure has wall connectivity. As shown in Fig. 13, the geometric constraint (18) was satisfied while the optimization process of fluid diodes proceeded in the same manner as the case without the geometric constraint. It can be seen that the floating solid regions in Fig. 9 disappeared in this solution since such structure is penalized by the geometric constraint, and instead a structure connecting to the outer boundaries emerged near the inlet of the reverse flow, which significantly changes the direction of the reverse flow. The value of the evaluation function $J_h$ of the obtained structure was $J_h \approx 9.1 \times 10^{-6}$. Fig. 14 shows the fictitious temperature field over solid regions. The fictitious temperature field of this result has a quite smaller value than that of the case without the constraint for excluding “floating” islands.

While the geometric constraint worked well in this example, we remark that it is non-trivial to

know how to set the upper bound value of the geometric constraint. If the upper bound value is not sufficiently small, the “floating” islands might not be penalized and remain in a solution. However, an overly strict geometric constraint value can lead to numerical instability or significantly affect the original optimization problem. Also, the result considering wall-connected structures has better performance in terms of diodicity, but this is not always the case. Both obtained structures in Figs. 9 and 12 are local optima of the problem, and in the latter case, the geometric constraint may have caused the structure to converge to another local optimum. It can also be observed that the solution contains elongated structures of several cell sizes. This structure might be attributed to the effect of the geometric constraint rather than maximizing diodicity. To prevent such small structures, one can employ the minimum length control technique. We hope to extend our proposed method to a robust methodology that can create three-dimensional manufacturable fluid diodes with high diodicity.

## 4.3 Full-scale three-dimensional design of fluid diodes

In this section, we present optimization results of full-scale three-dimensional design of fluid diodes with and without the geometric constraint. Fig. 15a shows the overall view of the channel to be considered in this numerical example. Fig. 15b shows the computational domain and the boundary conditions for the fluid problem. As can be seen, the target channel has a three-dimensional torsional structure, and unlike the previous section, symmetry boundary conditions cannot be used. The computational domain was discretized with $3.63 \cdot 10^5$ structured hexahedral cells. For parallel computing, the computational domain was decomposed into 20 subdomains. As in the previous section, a flow with $Re = 100$ is considered. Before optimization, the energy dissipation of the channel was evaluated and found to be approximately 2.08, assuming that the computational domain is composed of pure fluid phase. For this case, the upper bound value of the energy dissipation in the forward flow $J_{f,max}$ was set to 10. The maximum value of inverse permeability $\alpha_{max}$ was set to $1.0 \times 10^3$. For parameters of density filtering, the filter radius $r$ was set to $1.7 \times 10^{-2}$.
<figure><img src="image_11.png" /><figcaption>**Fig. 8:** Boundary conditions for the fictitious physical field</figcaption></figure>

holes appears to be less than that of the structure in Fig. 9, the resulting structure still has a unique topology. The energy dissipation in the reverse flow $J_r$ reached approximately 1.63, and the diodicity was $Di \approx 1.63$.

The difference between the obtained structures of Figs. 9 and 12 is the existence of floating solid regions. In the solution of Fig. 12, every solid part connects to the outer boundaries, and therefore the structure has wall connectivity. As shown in Fig. 13, the geometric constraint (18) was satisfied while the optimization process of fluid diodes proceeded in the same manner as the case without the geometric constraint. It can be seen that the floating solid regions in Fig. 9 disappeared in this solution since such structure is penalized by the geometric constraint, and instead a structure connecting to the outer boundaries emerged near the inlet of the reverse flow, which significantly changes the direction of the reverse flow. The value of the evaluation function $J_h$ of the obtained structure was $J_h \approx 9.1 \times 10^{-6}$. Fig. 14 shows the fictitious temperature field over solid regions. The fictitious temperature field of this result has a quite smaller value than that of the case without the constraint for excluding “floating” islands.

While the geometric constraint worked well in this example, we remark that it is non-trivial to

know how to set the upper bound value of the geometric constraint. If the upper bound value is not sufficiently small, the “floating” islands might not be penalized and remain in a solution. However, an overly strict geometric constraint value can lead to numerical instability or significantly affect the original optimization problem. Also, the result considering wall-connected structures has better performance in terms of diodicity, but this is not always the case. Both obtained structures in Figs. 9 and 12 are local optima of the problem, and in the latter case, the geometric constraint may have caused the structure to converge to another local optimum. It can also be observed that the solution contains elongated structures of several cell sizes. This structure might be attributed to the effect of the geometric constraint rather than maximizing diodicity. To prevent such small structures, one can employ the minimum length control technique. We hope to extend our proposed method to a robust methodology that can create three-dimensional manufacturable fluid diodes with high diodicity.

## 4.3 Full-scale three-dimensional design of fluid diodes

In this section, we present optimization results of full-scale three-dimensional design of fluid diodes with and without the geometric constraint. Fig. 15a shows the overall view of the channel to be considered in this numerical example. Fig. 15b shows the computational domain and the boundary conditions for the fluid problem. As can be seen, the target channel has a three-dimensional torsional structure, and unlike the previous section, symmetry boundary conditions cannot be used. The computational domain was discretized with $3.63 \times 10^5$ structured hexahedral cells. For parallel computing, the computational domain was decomposed into 20 subdomains. As in the previous section, a flow with $Re = 100$ is considered. Before optimization, the energy dissipation of the channel was evaluated and found to be approximately 2.08, assuming that the computational domain is composed of pure fluid phase. For this case, the upper bound value of the energy dissipation in the forward flow $J_{f,max}$ was set to 10. The maximum value of inverse permeability $\alpha_{max}$ was set to $1.0 \times 10^3$. For parameters of density filtering, the filter radius $r$ was set to $1.7 \times 10^{-2}$.
<figure><img src="image_9.png" /><figcaption>Fig. 9: Obtained configuration without the geometric constraint</figcaption></figure>

The initial guess was uniformly given as $\gamma = 0.5$ in the design domain.

Boundary conditions for the fictitious physical field are shown in Fig. 16. The parameters for the fictitious physical field were set as $d_{\text{max}} = 1 \times 10^2$ and $d_{\text{min}} = 1 \times 10^{-5}$.

### 4.3.1 Optimization results without the geometric constraint

Fig. 17 shows the obtained configuration of a fluid diode represented by the isovolume of $\hat{\gamma} \in [0, 0.5]$, where Figs. 17a and b show the front and side view

613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
<figure><img src="image_12.png" /><figcaption>**Fig. 9:** Obtained configuration without the geometric constraint</figcaption></figure>

The initial guess was uniformly given as $\gamma = 0.5$ in the design domain.

Boundary conditions for the fictitious physical field are shown in Fig. 16. The parameters for the fictitious physical field were set as $d_{\text{max}} = 1 \times 10^2$ and $d_{\text{min}} = 1 \times 10^{-5}$.

### 4.3.1 Optimization results without the geometric constraint

Fig. 17 shows the obtained configuration of a fluid diode represented by the isovolume of $\hat{\gamma} \in [0, 0.5]$, where Figs. 17a and b show the front and side view

613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
<figure><img src="image_10.png" /><figcaption>Fig. 10: Optimization history of the result without the geometric constraint</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 11: Fictitious physical field of obtained solution without the geometric constraint. The floating solid domains have high value of the fictitious temperature</figcaption></figure>

of obtained structure, and Figs. 17c and d show streamlines of the forward and reverse flows, all respectively. Fig. 18 shows the optimization history of $J_f$, $J_r$, and $J_h$. The convergence criteria was satisfied in 432 optimization iterations. As in the previous example, it can be seen that the reverse flow contains a relatively large number of vortices. In this example, the obtained structure had the diodicity with Di ≈ 3.53.

However, the obtained result has “floating” islands similar to what has been observed in the previous example in Section 4.2.1. Fig. 19

shows the fictitious physical field on the solid domains. The “floating” solids are detected with large value of T. The history of function $J_h$ is plotted in Fig. 18, and the value of $J_h$ for the final configuration was $J_h \approx 8.8 \times 10^{-3}$.

### 4.3.2 Optimization results with the geometric constraint

In order to obtain a solution without “floating” islands, we introduced the geometric constraint into this example. In this case, the geometric constraint was used as a post-processing step. That is, using the solution obtained without the geometric constraint in Section 4.3.1 as the initial structure, the optimization problem was solved by adding the geometric constraint. For the upper bound value $J_{h,max}$, $1 \times 10^{-5}$ was used. For convergence criteria, $\varepsilon_c = 1 \times 10^{-4}$ was used. The sharpness parameter $\beta$ was set to a constant value of 40.

In this example, it seems to be rather difficult to obtain the solution which fully satisfies the convergence criteria. Given that the objective value was somewhat stabilized after approximately 200 iterations with a slight fluctuation, we therefore terminated the optimization procedure at the 300th iteration. The final structure is shown in Fig. 20 represented by the isovolume of $\hat{\gamma} \in [0, 0.5]$, and the optimization history is plotted in Fig. 21. It shows that the obtained structure contained a relatively large amount of regions with greyscale density. Therefore, the performance of the solution was evaluated by binarizing the density with a threshold value of 0.5. The energy dissipation in the forward flow was $J_f \approx 9.9$, the energy dissipation in the reverse flow was $J_r \approx 34.3$, and the evaluation function of the fictitious field was $J_h \approx 2.9 \times 10^{-5}$. Although the constraint for the fictitious physical field $G_2$ was violated after binarization, the energy dissipation in the forward flow remained within the constraint. The value of diodicity was Di ≈ 3.47, a slight decrease from the initial structure, cf. Fig. 17.

Fig. 22 shows the fictitious temperature field over solid regions. Compared with the structure in Fig. 19, most of the “floating” islands were successfully eliminated. Nevertheless, there still remain several tiny solid domains without wall connectivity. This can be attributed to the nature of the constraint function Eq. (18). The geometric constraint limits the mean square of the
<figure><img src="image_13.png" /><figcaption>Fig. 10: Optimization history of the result without the geometric constraint</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 11: Fictitious physical field of obtained solution without the geometric constraint. The floating solid domains have high value of the fictitious temperature</figcaption></figure>

of obtained structure, and Figs. 17c and d show streamlines of the forward and reverse flows, all respectively. Fig. 18 shows the optimization history of $J_f$, $J_r$, and $J_h$. The convergence criteria was satisfied in 432 optimization iterations. As in the previous example, it can be seen that the reverse flow contains a relatively large number of vortices. In this example, the obtained structure had the diodicity with $Di \approx 3.53$.

However, the obtained result has “floating” islands similar to what has been observed in the previous example in Section 4.2.1. Fig. 19

shows the fictitious physical field on the solid domains. The “floating” solids are detected with large value of $T$. The history of function $J_h$ is plotted in Fig. 18, and the value of $J_h$ for the final configuration was $J_h \approx 8.8 \times 10^{-3}$.

### 4.3.2 Optimization results with the geometric constraint

In order to obtain a solution without “floating” islands, we introduced the geometric constraint into this example. In this case, the geometric constraint was used as a post-processing step. That is, using the solution obtained without the geometric constraint in Section 4.3.1 as the initial structure, the optimization problem was solved by adding the geometric constraint. For the upper bound value $J_{h,max}$, $1 \times 10^{-5}$ was used. For convergence criteria, $\varepsilon_c = 1 \times 10^{-4}$ was used. The sharpness parameter $\beta$ was set to a constant value of 40.

In this example, it seems to be rather difficult to obtain the solution which fully satisfies the convergence criteria. Given that the objective value was somewhat stabilized after approximately 200 iterations with a slight fluctuation, we therefore terminated the optimization procedure at the 300th iteration. The final structure is shown in Fig. 20 represented by the isovolume of $\hat{\gamma} \in [0, 0.5]$, and the optimization history is plotted in Fig. 21. It shows that the obtained structure contained a relatively large amount of regions with greyscale density. Therefore, the performance of the solution was evaluated by binarizing the density with a threshold value of 0.5. The energy dissipation in the forward flow was $J_f \approx 9.9$, the energy dissipation in the reverse flow was $J_r \approx 34.3$, and the evaluation function of the fictitious field was $J_h \approx 2.9 \times 10^{-5}$. Although the constraint for the fictitious physical field $G_2$ was violated after binarization, the energy dissipation in the forward flow remained within the constraint. The value of diodicity was $Di \approx 3.47$, a slight decrease from the initial structure, cf. Fig. 17.

Fig. 22 shows the fictitious temperature field over solid regions. Compared with the structure in Fig. 19, most of the “floating” islands were successfully eliminated. Nevertheless, there still remain several tiny solid domains without wall connectivity. This can be attributed to the nature of the constraint function Eq. (18). The geometric constraint limits the mean square of the
<figure><img src="image_12.png" /><figcaption>(a) Quarter model.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(b) Half model.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(c) Streamlines of forward flow.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(d) Streamlines of reverse flow.</figcaption></figure>

<caption>Fig. 12: Obtained configuration with the geometric constraint $J_{h,max} = 1 \times 10^{-5}$</caption>

fictitious temperature over the solid domain. On the other hand, although small “floating” islands have high fictitious temperatures, their contribution to the mean square temperature is rather insignificant due to their small volumes. As an alternative, an upper constraint on the maximum

temperature could be considered, which might effectively exclude “floating” islands by setting an appropriate upper bound temperature.

Furthermore, mesh-size structures are observed in Fig. 19, and this may also be due to the geometric constraint. To mitigate the

715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
<figure><img src="image_15.png" /><figcaption>(a) Quarter model.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(b) Half model.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(c) Streamlines of forward flow.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(d) Streamlines of reverse flow.</figcaption></figure>

<caption>**Fig. 12:** Obtained configuration with the geometric constraint $J_{h,max} = 1 \times 10^{-5}$</caption>

fictitious temperature over the solid domain. On the other hand, although small “floating” islands have high fictitious temperatures, their contribution to the mean square temperature is rather insignificant due to their small volumes. As an alternative, an upper constraint on the maximum

temperature could be considered, which might effectively exclude “floating” islands by setting an appropriate upper bound temperature.

Furthermore, mesh-size structures are observed in Fig. 19, and this may also be due to the geometric constraint. To mitigate the

715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
<figure><img src="image_16.png" /><figcaption>Fig. 13: Optimization history of the result with the geometric constraint $J_{h,max} = 1 \times 10^{-5}$</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 14: Fictitious physical field of obtained solution with the geometric constraint $J_{h,max} = 1 \times 10^{-5}$</figcaption></figure>

mesh size structures, an additional constraint or further regularization would be necessary. As mentioned in Section 4.2.2, the minimum length control technique will work effectively. Employing a larger filter radius or introducing a volume constraint could also be successful. A volume constraint can usually regulate the optimization problem, making it easier to obtain a fluid diode without “floating” islands.

<figure><img src="image_18.png" /><figcaption>(a) Entire domain</figcaption></figure>

$\Gamma_{\text{out,f}}$: Outlet of forward flow

$\Gamma_{\text{in,r}}$: Inlet of reverse flow

<figure><img src="image_19.png" /></figure>

$\Gamma_{\text{in,f}}$: Inlet of forward flow

$\Gamma_{\text{out,r}}$: Outlet of reverse flow

(b) Design domain and boundary conditions for fluid flow

<caption>Fig. 15: Computational domain and boundary conditions for the fluid flow of full-scale three-dimensional case</caption>
<figure><img src="image_19.png" /><figcaption>Fig. 13: Optimization history of the result with the geometric constraint $J_{h,max} = 1 \times 10^{-5}$</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 14: Fictitious physical field of obtained solution with the geometric constraint $J_{h,max} = 1 \times 10^{-5}$</figcaption></figure>

mesh size structures, an additional constraint or further regularization would be necessary. As mentioned in Section 4.2.2, the minimum length control technique will work effectively. Employing a larger filter radius or introducing a volume constraint could also be successful. A volume constraint can usually regulate the optimization problem, making it easier to obtain a fluid diode without “floating” islands.

<figure><img src="image_21.png" /><figcaption>(a) Entire domain</figcaption></figure>

$\Gamma_{\text{out,f}}$: Outlet of forward flow

$\Gamma_{\text{in,r}}$: Inlet of reverse flow

<figure><img src="image_22.png" /></figure>

$\Gamma_{\text{in,f}}$: Inlet of forward flow

$\Gamma_{\text{out,r}}$: Outlet of reverse flow

(b) Design domain and boundary conditions for fluid flow

<caption>Fig. 15: Computational domain and boundary conditions for the fluid flow of full-scale three-dimensional case</caption>
<figure><img src="image_20.png" /><figcaption>Fig. 16: Boundary conditions for the fictitious physical field</figcaption></figure>

## 5 Conclusions

In this paper, we proposed a density-based topology optimization method for the three-dimensional design of fluid diodes considering wall-connected structures based on the fictitious physical model approach. One possible problem of three-dimensional design of fluid diodes was the emergence of the “floating” solid domains, the structures without wall connectivity. To deal with this issue, we formulated a fictitious physical model and a geometric constraint to detect and restrict the “floating” solid domains. The main findings are summarized as follows:

1. The optimum design problem of fluid diodes was formulated as maximizing the energy dissipation in the reverse flow subject to the upper bound constraint of the energy dissipation in the forward flow. To detect “floating” islands, the fictitious physical problem was introduced so that the solution has value of $O(1)$ in the “floating” islands and a near-zero value in the other regions. To restrict the “floating” islands, the upper bound constraint was introduced for the integral of the square of the fictitious temperature in the solid domains. The sensitivities were derived based on continuous adjoint

method and computed by solving primal and adjoint equations. All the numerical simulations were implemented using OpenFOAM.

2. The proposed fictitious physical model was investigated through the two-dimensional example. As a result, our model has proved to be able to detect the “floating” solids by adjusting the maximum and the minimum value of the diffusion coefficient to an appropriate value. Unlike the fictitious model in previous work, the fictitious temperature has a larger value in large solid domains than in smaller ones. Also, our proposed constraint function for excluding floating solids could successfully eliminate them during the optimization.

3. We first performed topology optimization for designing three-dimensional fluid diodes without the geometric constraints in the two types of channels. In both cases, the obtained structures showed diodicity of $Di > 1$ within the constraint for the energy dissipation of forward flow. Nevertheless, our numerical examples demonstrated that the solutions of three-dimensional fluid diodes can contain “floating” solid islands and therefore they cannot be manufactured. Meanwhile, those “floating” islands were detected clearly through the proposed fictitious physical model.

4. We then applied our formulation to design fluid diodes considering wall-connected structures. In the symmetric test case, the optimization procedure was initialized by a uniform density distribution, and in the non-symmetric test case, we performed optimization from the obtained result without the geometric constraint. As a result, our method obtained the solutions without “floating” islands while the relatively high diodicity was maintained. Although we observed that mesh-size structures emerged as a side effect of the geometric constraint, the solutions had similar performance to those without wall connectivity and had no or much smaller “floating” islands.

817

818

819

820

821

822

823

824

825

826

827

828

829

830

831

832

833

834

835

836

837

838

839

840

841

842

843

844

845

846

847

848

849

850

851

852

853

854

855

856

857

858

859

860

861

862

863

864

865

866

867
<figure><img src="image_23.png" /><figcaption>**Fig. 16:** Boundary conditions for the fictitious physical field</figcaption></figure>

## 5 Conclusions

In this paper, we proposed a density-based topology optimization method for the three-dimensional design of fluid diodes considering wall-connected structures based on the fictitious physical model approach. One possible problem of three-dimensional design of fluid diodes was the emergence of the “floating” solid domains, the structures without wall connectivity. To deal with this issue, we formulated a fictitious physical model and a geometric constraint to detect and restrict the “floating” solid domains. The main findings are summarized as follows:

1. The optimum design problem of fluid diodes was formulated as maximizing the energy dissipation in the reverse flow subject to the upper bound constraint of the energy dissipation in the forward flow. To detect “floating” islands, the fictitious physical problem was introduced so that the solution has value of $O(1)$ in the “floating” islands and a near-zero value in the other regions. To restrict the “floating” islands, the upper bound constraint was introduced for the integral of the square of the fictitious temperature in the solid domains. The sensitivities were derived based on continuous adjoint

method and computed by solving primal and adjoint equations. All the numerical simulations were implemented using OpenFOAM.

2. The proposed fictitious physical model was investigated through the two-dimensional example. As a result, our model has proved to be able to detect the “floating” solids by adjusting the maximum and the minimum value of the diffusion coefficient to an appropriate value. Unlike the fictitious model in previous work, the fictitious temperature has a larger value in large solid domains than in smaller ones. Also, our proposed constraint function for excluding floating solids could successfully eliminate them during the optimization.

3. We first performed topology optimization for designing three-dimensional fluid diodes without the geometric constraints in the two types of channels. In both cases, the obtained structures showed diodicity of $Di > 1$ within the constraint for the energy dissipation of forward flow. Nevertheless, our numerical examples demonstrated that the solutions of three-dimensional fluid diodes can contain “floating” solid islands and therefore they cannot be manufactured. Meanwhile, those “floating” islands were detected clearly through the proposed fictitious physical model.

4. We then applied our formulation to design fluid diodes considering wall-connected structures. In the symmetric test case, the optimization procedure was initialized by a uniform density distribution, and in the non-symmetric test case, we performed optimization from the obtained result without the geometric constraint. As a result, our method obtained the solutions without “floating” islands while the relatively high diodicity was maintained. Although we observed that mesh-size structures emerged as a side effect of the geometric constraint, the solutions had similar performance to those without wall connectivity and had no or much smaller “floating” islands.

817

818

819

820

821

822

823

824

825

826

827

828

829

830

831

832

833

834

835

836

837

838

839

840

841

842

843

844

845

846

847

848

849

850

851

852

853

854

855

856

857

858

859

860

861

862

863

864

865

866

867
<figure><img src="image_21.png" /><figcaption>Fig. 17: Obtained configuration without the geometric constraint in full-scale 3D case</figcaption></figure>

868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
<figure><img src="image_24.png" /><figcaption>Fig. 17: Obtained configuration without the geometric constraint in full-scale 3D case</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Fig. 18:** Optimization history of the result without the geometric constraint in full-scale 3D case</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Fig. 19:** Fictitious physical field of obtained solution without the geometric constraint</figcaption></figure>

919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
<figure><img src="image_25.png" /><figcaption>**Fig. 18:** Optimization history of the result without the geometric constraint in full-scale 3D case</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>**Fig. 19:** Fictitious physical field of obtained solution without the geometric constraint</figcaption></figure>

919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
<figure><img src="image_24.png" /><figcaption>(a) Front view</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>(c) Streamlines of forward flow</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>(d) Streamlines of reverse flow</figcaption></figure>

<caption>**Fig. 20:** Obtained configuration with the geometric constraint in full-scale 3D case</caption>

970

971

972

973

974

975

976

977

978

979

980

981

982

983

984

985

986

987

988

989

990

991

992

993

994

995

996

997

998

999

1000

1001

1002

1003

1004

1005

1006

1007

1008

1009

1010

1011

1012

1013

1014

1015

1016

1017

1018

1019

1020
<figure><img src="image_27.png" /><figcaption>(a) Front view</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>(c) Streamlines of forward flow</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>(d) Streamlines of reverse flow</figcaption></figure>

<caption>**Fig. 20:** Obtained configuration with the geometric constraint in full-scale 3D case</caption>

970

971

972

973

974

975

976

977

978

979

980

981

982

983

984

985

986

987

988

989

990

991

992

993

994

995

996

997

998

999

1000

1001

1002

1003

1004

1005

1006

1007

1008

1009

1010

1011

1012

1013

1014

1015

1016

1017

1018

1019

1020
<figure><img src="image_27.png" /><figcaption>Fig. 21: Optimization history of the result with the geometric constraint in full-scale 3D case</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Fig. 22: Fictitious physical field of obtained solution with the geometric constraint $J_{h,max} = 1 \times 10^{-5}$</figcaption></figure>

1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
<figure><img src="image_30.png" /><figcaption>Fig. 21: Optimization history of the result with the geometric constraint in full-scale 3D case</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>Fig. 22: Fictitious physical field of obtained solution with the geometric constraint $J_{h,max} = 1 \times 10^{-5}$</figcaption></figure>

1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
## Appendix A Expressions of the energy dissipation

In this section, we describe the conservation of the mechanical energy of the Newtonian fluid in the control volume $\Omega$ shown in Fig. A.1. Similar explanations can be seen in [Lin et al., 2015]. The governing equations of a steady-state incompressible flow with Brinkman penalization are described as follows:

$$ \begin{cases} -\nabla \cdot (\boldsymbol{u}\boldsymbol{u}^{\top}) - \nabla p + \nabla \cdot \bar{\boldsymbol{\tau}} - \alpha \boldsymbol{u} = \mathbf{0} & \text{in } \Omega \\ \nabla \cdot \boldsymbol{u} = 0 & \text{in } \Omega, \end{cases} \quad (A.1) $$

where $\bar{\boldsymbol{\tau}}$ denotes the viscous stress tensor defined as,

$$ \bar{\boldsymbol{\tau}} := \frac{1}{\text{Re}} (\nabla \boldsymbol{u} + \nabla \boldsymbol{u}^{\top}). \quad (A.2) $$

Multiplying the momentum equation by the velocity $\boldsymbol{u}$ as inner product yields, after some manipulations,

$$ \nabla \cdot \left\{ -\left(p + \frac{u^2}{2}\right) \boldsymbol{u} \right\} + \nabla \cdot \boldsymbol{\tau} = \Phi + \alpha \boldsymbol{u}^2, \quad (A.3) $$

where $\boldsymbol{\tau} := \bar{\boldsymbol{\tau}} \cdot \boldsymbol{n}$, and $\Phi$ is the dissipation function defined as,

$$ \Phi := \bar{\boldsymbol{\tau}} : \nabla \boldsymbol{u} = \frac{1}{2\text{Re}} (\nabla \boldsymbol{u} + \nabla \boldsymbol{u}^{\top})^2 . \quad (A.4) $$

<figure><img src="image_29.png" /><figcaption>Fig. A.1: Control volume $\Omega$ for considering energy conservation</figcaption></figure>

Taking the domain integral of (A.3) and applying the divergence theorem yields the following relation.

$$ \begin{aligned} & - \int_{\Gamma} \left( p + \frac{u^2}{2} \right) (\boldsymbol{u} \cdot \boldsymbol{n}) \, d\Gamma + \int_{\Gamma} \boldsymbol{\tau} \cdot \boldsymbol{u} \, d\Gamma \\ &= \int_{\Omega} (\Phi + \alpha \boldsymbol{u}^2) \, d\Omega. \end{aligned} \quad (A.5) $$

The boundary $\Gamma = \partial\Omega$ of the control volume $\Omega$ consists of inlet $\Gamma_{in}$, outlet $\Gamma_{out}$ and wall $\Gamma_{wall}$. On the wall $\Gamma_{wall}$, the velocity is equal to zero.

$$ \boldsymbol{u} = \mathbf{0} \text{ on } \Gamma_{\text{wall}}. \quad (A.6) $$

Therefore, the equation (A.5) can be transformed to:

$$ \begin{aligned} & - \int_{\Gamma_{\text{in}} \cup \Gamma_{\text{out}}} \left( p + \frac{u^2}{2} \right) (\boldsymbol{u} \cdot \boldsymbol{n}) \, d\Gamma + \int_{\Gamma_{\text{in}} \cup \Gamma_{\text{out}}} \boldsymbol{\tau} \cdot \boldsymbol{u} \, d\Gamma \\ &= \int_{\Omega} (\Phi + \alpha \boldsymbol{u}^2) \, d\Omega. \end{aligned} \quad (A.7) $$

The second term in the left-hand side of (A.7) implies work done by viscous stress on the inlet and outlet boundaries. Although it can be considered as negligible in general, this term can be also theoretically neglected if fully developed flow is assumed. Then, the following relation is obtained.

$$ - \int_{\Gamma_{\text{in}} \cup \Gamma_{\text{out}}} \left( p + \frac{u^2}{2} \right) (\boldsymbol{u} \cdot \boldsymbol{n}) \, d\Gamma = \int_{\Omega} (\Phi + \alpha \boldsymbol{u}^2) \, d\Omega. \quad (A.8) $$

Here, we have two equivalent expressions of energy dissipation. Therefore, the definition of diodicity Di' in the previous works [Lin et al., 2015, Sato et al., 2017a] and the definition in this work (15) are equivalent. Additionally, if the shapes of the inlet and the outlet are identical, the fully developed velocity profiles on the inlet and the outlet are the same, and thus the total flux of the kinematic energies is equal to zero:

$$ \int_{\Gamma_{\text{in}}} \frac{u^2}{2} (\boldsymbol{u} \cdot \boldsymbol{n}) \, d\Gamma + \int_{\Gamma_{\text{out}}} \frac{u^2}{2} (\boldsymbol{u} \cdot \boldsymbol{n}) \, d\Gamma = 0 \quad (A.9) $$
## Appendix A Expressions of the energy dissipation

In this section, we describe the conservation of the mechanical energy of the Newtonian fluid in the control volume $\Omega$ shown in Fig. A.1. Similar explanations can be seen in [Lin et al., 2015]. The governing equations of a steady-state incompressible flow with Brinkman penalization are described as follows:

$$ \begin{cases} -\nabla \cdot (\boldsymbol{u}\boldsymbol{u}^{\top}) - \nabla p + \nabla \cdot \bar{\boldsymbol{\tau}} - \alpha \boldsymbol{u} = \mathbf{0} & \text{in } \Omega \\ \nabla \cdot \boldsymbol{u} = 0 & \text{in } \Omega, \end{cases} \quad \text{(A.1)} $$

where $\bar{\boldsymbol{\tau}}$ denotes the viscous stress tensor defined as,

$$ \bar{\boldsymbol{\tau}} := \frac{1}{\text{Re}} (\nabla \boldsymbol{u} + \nabla \boldsymbol{u}^{\top}). \quad \text{(A.2)} $$

Multiplying the momentum equation by the velocity $\boldsymbol{u}$ as inner product yields, after some manipulations,

$$ \nabla \cdot \left\{ -\left(p + \frac{u^2}{2}\right) \boldsymbol{u} \right\} + \nabla \cdot \boldsymbol{\tau} = \Phi + \alpha \boldsymbol{u}^2, \quad \text{(A.3)} $$

where $\boldsymbol{\tau} := \bar{\boldsymbol{\tau}} \cdot \boldsymbol{n}$, and $\Phi$ is the dissipation function defined as,

$$ \Phi := \bar{\boldsymbol{\tau}} : \nabla \boldsymbol{u} = \frac{1}{2\text{Re}} (\nabla \boldsymbol{u} + \nabla \boldsymbol{u}^{\top})^2 . \quad \text{(A.4)} $$

<figure><img src="image_32.png" /><figcaption>Fig. A.1: Control volume $\Omega$ for considering energy conservation</figcaption></figure>

Taking the domain integral of (A.3) and applying the divergence theorem yields the following relation.

$$ \begin{aligned} & - \int_{\Gamma} \left( p + \frac{u^2}{2} \right) (\boldsymbol{u} \cdot \boldsymbol{n}) \, d\Gamma + \int_{\Gamma} \boldsymbol{\tau} \cdot \boldsymbol{u} \, d\Gamma \\ &= \int_{\Omega} (\Phi + \alpha \boldsymbol{u}^2) \, d\Omega. \end{aligned} \quad \text{(A.5)} $$

The boundary $\Gamma = \partial\Omega$ of the control volume $\Omega$ consists of inlet $\Gamma_{\text{in}}$, outlet $\Gamma_{\text{out}}$ and wall $\Gamma_{\text{wall}}$. On the wall $\Gamma_{\text{wall}}$, the velocity is equal to zero.

$$ \boldsymbol{u} = \mathbf{0} \text{ on } \Gamma_{\text{wall}}. \quad \text{(A.6)} $$

Therefore, the equation (A.5) can be transformed to:

$$ \begin{aligned} & - \int_{\Gamma_{\text{in}} \cup \Gamma_{\text{out}}} \left( p + \frac{u^2}{2} \right) (\boldsymbol{u} \cdot \boldsymbol{n}) \, d\Gamma + \int_{\Gamma_{\text{in}} \cup \Gamma_{\text{out}}} \boldsymbol{\tau} \cdot \boldsymbol{u} \, d\Gamma \\ &= \int_{\Omega} (\Phi + \alpha \boldsymbol{u}^2) \, d\Omega. \end{aligned} \quad \text{(A.7)} $$

The second term in the left-hand side of (A.7) implies work done by viscous stress on the inlet and outlet boundaries. Although it can be considered as negligible in general, this term can be also theoretically neglected if fully developed flow is assumed. Then, the following relation is obtained.

$$ - \int_{\Gamma_{\text{in}} \cup \Gamma_{\text{out}}} \left( p + \frac{u^2}{2} \right) (\boldsymbol{u} \cdot \boldsymbol{n}) \, d\Gamma = \int_{\Omega} (\Phi + \alpha \boldsymbol{u}^2) \, d\Omega. \quad \text{(A.8)} $$

Here, we have two equivalent expressions of energy dissipation. Therefore, the definition of diodicity Di' in the previous works [Lin et al., 2015, Sato et al., 2017a] and the definition in this work (15) are equivalent. Additionally, if the shapes of the inlet and the outlet are identical, the fully developed velocity profiles on the inlet and the outlet are the same, and thus the total flux of the kinematic energies is equal to zero:

$$ \int_{\Gamma_{\text{in}}} \frac{u^2}{2} (\boldsymbol{u} \cdot \boldsymbol{n}) \, d\Gamma + \int_{\Gamma_{\text{out}}} \frac{u^2}{2} (\boldsymbol{u} \cdot \boldsymbol{n}) \, d\Gamma = 0 \quad \text{(A.9)} $$
Therefore, the left-hand side of (A.8) can be rewritten as,

$$
\begin{align}
& - \int_{\Gamma_{\text{in}} \cup \Gamma_{\text{out}}} \left( p + \frac{u^2}{2} \right) (\mathbf{u} \cdot \mathbf{n}) \, d\Gamma \nonumber \\
&= - \int_{\Gamma_{\text{in}} \cup \Gamma_{\text{out}}} p (\mathbf{u} \cdot \mathbf{n}) \, d\Gamma \tag{A.10} \\
&= (p_{\text{in}} - p_{\text{out}}) \cdot Q = \Delta p \cdot Q, \nonumber
\end{align}
$$

where $\Delta p$ is the static pressure drop and $Q := \int_{\Gamma_{\text{in}}} (-\mathbf{u} \cdot \mathbf{n}) \, d\Gamma = \int_{\Gamma_{\text{out}}} \mathbf{u} \cdot \mathbf{n} \, d\Gamma$ is the flow rate. Thus, from equations (14) and (15), the two expressions of diodicity $Di$ and $Di'$ are equivalent.

Appendix B Sensitivity analysis

In this section, the derivation of adjoint systems and sensitivities are briefly described. For simplicity, we use the following expressions of the governing equations in the residual form.

$$
\left\{
\begin{aligned}
\mathcal{R}_u(\mathbf{u}, p, \gamma) &= \nabla \cdot (\mathbf{u}\mathbf{u}^\top) - \nabla \cdot \bar{\sigma}(\mathbf{u}, p) + \alpha(\gamma)\mathbf{u} = \mathbf{0} \\
R_p(\mathbf{u}, \gamma) &= -\nabla \cdot \mathbf{u} = \mathbf{0} \\
R_T(T, \gamma) &= \nabla \cdot (d(\gamma)\nabla T) + (1 - \hat{\gamma} - T) = \mathbf{0}.
\end{aligned}
\right.
\tag{B.1}
$$

The extended objective function $\mathcal{J}(\mathbf{u}_r, p_r, \gamma)$ and constraint functions $\mathcal{G}_1(\mathbf{u}_f, p_f, \gamma)$ and $\mathcal{G}_2(\mathbf{u}_f, p_f, \gamma)$ are defined as follows:

$$
\begin{align}
\mathcal{J}(\mathbf{u}_r, p_r, \gamma) &= J(\mathbf{u}_r, p_r) \nonumber \\
&\quad + \int_{\Omega} (-\mathbf{R}_u(\mathbf{u}_r, p_r, \gamma)) \cdot \mathbf{u}_r' \, d\Omega \nonumber \\
&\quad + \int_{\Omega} (-R_p(\mathbf{u}_r, \gamma)) p_r' \, d\Omega \nonumber \\
\mathcal{G}_1(\mathbf{u}_f, p_f, \gamma) &= G_1(\mathbf{u}_f, p_f) \nonumber \\
&\quad + \int_{\Omega} \mathbf{R}_u(\mathbf{u}_f, p_f, \gamma) \cdot \mathbf{u}_f' \, d\Omega \nonumber \\
&\quad + \int_{\Omega} R_p(\mathbf{u}_f, \gamma) p_f' \, d\Omega \nonumber \\
\mathcal{G}_2(T, \gamma) &= G_2(T, \gamma) \nonumber \\
&\quad + \int_{\Omega} R_T(T, \gamma) T' \, d\Omega \tag{B.2}
\end{align}
$$

where $\mathbf{u}_f'$ and $p_f'$ are adjoint variables for the forward flow, $\mathbf{u}_r'$ and $p_r'$ are adjoint variables for the

reverse flow, and $T'$ is the adjoint variable for the
fictitious physical field. The residuals are multi-
plied by -1 in the definition of $\mathcal{J}$ for simplicity of
the following derivation.

The adjoint equations are derived from the
conditions that partial variations of the extended
objective and constraint functions with respect to
the state variables are zero. The resultant adjoint
systems are described as follows:

$$
\left\{
\begin{array}{ll}
- (\nabla \mathbf{u}_r' + \nabla \mathbf{u}_r'^{\top}) \cdot \mathbf{u}_r & \quad \text{in } \Omega \\
- \nabla \cdot \bar{\sigma} (\mathbf{u}_r', p_r') + \alpha(\hat{\gamma}) \mathbf{u}_r' = \mathbf{0} & \quad \text{in } \Omega \\
\nabla \cdot \mathbf{u}_r' = 0 & \quad \text{in } \Omega \\
\mathbf{u}_r' = \mathbf{u}_D & \quad \text{on } \Gamma_D \\
\bar{\sigma} (\mathbf{u}_r', p_r') \cdot \mathbf{n} & \quad \text{in } \Gamma_N, \\
\frac{\mathbf{u}_r^2}{2} + (\mathbf{u}_r \cdot \mathbf{n}) \mathbf{u}_r & \quad \text{on } \Gamma_N, \\
- (\mathbf{u}_r \cdot \mathbf{n}) \mathbf{u}_r' - (\mathbf{u}_r \cdot \mathbf{u}_r') \mathbf{n} & \quad \text{in } \Omega
\end{array}
\right.
\tag{B.3}
$$

$$
\left\{
\begin{array}{ll}
- (\nabla \mathbf{u}_f' + \nabla \mathbf{u}_f'^{\top}) \cdot \mathbf{u}_f & \quad \text{in } \Omega \\
- \nabla \cdot \bar{\sigma} (\mathbf{u}_f', p_f') + \alpha(\hat{\gamma}) \mathbf{u}_f' = \mathbf{0} & \quad \text{in } \Omega \\
\nabla \cdot \mathbf{u}_f' = 0 & \quad \text{in } \Omega \\
\mathbf{u}_f' = \mathbf{u}_D & \quad \text{on } \Gamma_D \\
\bar{\sigma} (\mathbf{u}_f', p_f') \cdot \mathbf{n} & \quad \text{in } \Gamma_N, \\
\frac{\mathbf{u}_f^2}{2} + (\mathbf{u}_f \cdot \mathbf{n}) \mathbf{u}_f & \quad \text{on } \Gamma_N, \\
- (\mathbf{u}_f \cdot \mathbf{n}) \mathbf{u}_f' - (\mathbf{u}_f \cdot \mathbf{u}_f') \mathbf{n} & \quad \text{in } \Omega
\end{array}
\right.
\tag{B.4}
$$

$$
\left\{
\begin{array}{ll}
\nabla \cdot (d(\hat{\gamma}) \nabla T') - T' + 2(1-\hat{\gamma})T = 0 & \text{in } \Omega \\
T' = 0 & \text{on } \Gamma_T \\
d\nabla T' \cdot \mathbf{n} = 0 & \text{on } \Gamma_q.
\end{array}
\right.
\tag{B.5}
$$

The sensitivity of the objective function $J$ with
respect to the design variable $\gamma$ can be described
as follows:

$$
\frac{\delta \mathcal{J}}{\delta \hat{\gamma}} = \frac{\delta \mathcal{J} \partial \hat{\gamma} \delta \tilde{\gamma}}{\delta \tilde{\gamma} \partial \tilde{\gamma} \delta \gamma}, \quad (B.6)
$$

where the functional derivative of $\mathcal{J}$ with respect
to the projected density $\hat{\gamma}$ is derived as:

$$
\frac{\delta \mathcal{J}}{\delta \hat{\gamma}} = - \frac{\partial \alpha}{\partial \hat{\gamma}} (\mathbf{u}_r \cdot \mathbf{u}_r') = \alpha_{\max} (\mathbf{u}_r \cdot \mathbf{u}_r') \text{ in } \Omega, \quad (B.7)
$$

and the term $\frac{\partial \hat{\gamma}}{\partial \tilde{\gamma}}$ is the derivative of the smoothed
Heaviside function (6) with respect to the filtered
density $\tilde{\gamma}$:

$$
\frac{\partial \hat{\gamma}}{\partial \tilde{\gamma}} = \frac{\beta (1 - \tanh^2(\beta(\tilde{\gamma} - \eta)))}{\tanh(\beta\eta) + \tanh(\beta(1 - \eta))}. \quad (B.8)
$$
Therefore, the left-hand side of (A.8) can be rewritten as,

$$ 
\begin{aligned}
& - \int_{\Gamma_{\text{in}} \cup \Gamma_{\text{out}}} \left( p + \frac{u^2}{2} \right) (\mathbf{u} \cdot \mathbf{n}) \, d\Gamma \\
&= - \int_{\Gamma_{\text{in}} \cup \Gamma_{\text{out}}} p (\mathbf{u} \cdot \mathbf{n}) \, d\Gamma \\
&= (p_{\text{in}} - p_{\text{out}}) \cdot Q = \Delta p \cdot Q,
\end{aligned} 
\quad (A.10) $$

where $\Delta p$ is the static pressure drop and $Q := \int_{\Gamma_{\text{in}}} (-\mathbf{u} \cdot \mathbf{n}) \, d\Gamma = \int_{\Gamma_{\text{out}}} \mathbf{u} \cdot \mathbf{n} \, d\Gamma$ is the flow rate. Thus, from equations (14) and (15), the two expressions of diodicity Di and Di' are equivalent.

## Appendix B Sensitivity analysis

In this section, the derivation of adjoint systems and sensitivities are briefly described. For simplicity, we use the following expressions of the governing equations in the residual form.

$$ 
\left\{
\begin{aligned}
& \mathbf{R}_u(\mathbf{u}, p, \gamma) \\
& = \nabla \cdot (\mathbf{u}\mathbf{u}^\top) - \nabla \cdot \bar{\sigma}(\mathbf{u}, p) + \alpha(\gamma)\mathbf{u} = \mathbf{0} \\
& R_p(\mathbf{u}, \gamma) = -\nabla \cdot \mathbf{u} = \mathbf{0} \\
& R_T(T, \gamma) \\
& = \nabla \cdot (d(\gamma)\nabla T) + (1 - \hat{\gamma} - T) = \mathbf{0}.
\end{aligned}
\right.
\quad (B.1) $$

The extended objective function $\mathcal{J}(\mathbf{u}_r, p_r, \gamma)$ and constraint functions $\mathcal{G}_1(\mathbf{u}_f, p_f, \gamma)$ and $\mathcal{G}_2(\mathbf{u}_f, p_f, \gamma)$ are defined as follows:

$$ 
\begin{aligned}
\mathcal{J}(\mathbf{u}_r, p_r, \gamma) &= J(\mathbf{u}_r, p_r) \\
&\quad + \int_{\Omega} (-\mathbf{R}_u(\mathbf{u}_r, p_r, \gamma)) \cdot \mathbf{u}_r' \, d\Omega \\
&\quad + \int_{\Omega} (-R_p(\mathbf{u}_r, \gamma)) p_r' \, d\Omega \\
\mathcal{G}_1(\mathbf{u}_f, p_f, \gamma) &= G_1(\mathbf{u}_f, p_f) \\
&\quad + \int_{\Omega} \mathbf{R}_u(\mathbf{u}_f, p_f, \gamma) \cdot \mathbf{u}_f' \, d\Omega \\
&\quad + \int_{\Omega} R_p(\mathbf{u}_f, \gamma) p_f' \, d\Omega \\
\mathcal{G}_2(T, \gamma) &= G_2(T, \gamma) \\
&\quad + \int_{\Omega} R_T(T, \gamma) T' \, d\Omega
\end{aligned}
\quad (B.2) $$

where $\mathbf{u}_f'$ and $p_f'$ are adjoint variables for the forward flow, $\mathbf{u}_r'$ and $p_r'$ are adjoint variables for the

reverse flow, and $T'$ is the adjoint variable for the fictitious physical field. The residuals are multiplied by $-1$ in the definition of $\mathcal{J}$ for simplicity of the following derivation.

The adjoint equations are derived from the conditions that partial variations of the extended objective and constraint functions with respect to the state variables are zero. The resultant adjoint systems are described as follows:

$$ 
\left\{
\begin{aligned}
& - (\nabla \mathbf{u}_r' + \nabla \mathbf{u}_r'^{\top}) \cdot \mathbf{u}_r && \text{in } \Omega \\
& - \nabla \cdot \bar{\sigma} (\mathbf{u}_r', p_r') + \alpha(\hat{\gamma}) \mathbf{u}_r' = \mathbf{0} && \text{in } \Omega \\
& \nabla \cdot \mathbf{u}_r' = 0 && \text{in } \Omega \\
& \mathbf{u}_r' = \mathbf{u}_D && \text{on } \Gamma_D \\
& \bar{\sigma} (\mathbf{u}_r', p_r') \cdot \mathbf{n} && \text{on } \Gamma_N, \\
& \frac{\mathbf{u}_r'^2}{2} + (\mathbf{u}_r \cdot \mathbf{n}) \mathbf{u}_r && \text{on } \Gamma_N, \\
& - (\mathbf{u}_r \cdot \mathbf{n}) \mathbf{u}_r' - (\mathbf{u}_r \cdot \mathbf{u}_r') \mathbf{n} && \text{on } \Gamma_N, \\
\end{aligned}
\right.
\quad (B.3) $$

$$ 
\left\{
\begin{aligned}
& - (\nabla \mathbf{u}_f' + \nabla \mathbf{u}_f'^{\top}) \cdot \mathbf{u}_f && \text{in } \Omega \\
& - \nabla \cdot \bar{\sigma} (\mathbf{u}_f', p_f') + \alpha(\hat{\gamma}) \mathbf{u}_f' = \mathbf{0} && \text{in } \Omega \\
& \nabla \cdot \mathbf{u}_f' = 0 && \text{in } \Omega \\
& \mathbf{u}_f' = \mathbf{u}_D && \text{on } \Gamma_D \\
& \bar{\sigma} (\mathbf{u}_f', p_f') \cdot \mathbf{n} && \text{on } \Gamma_N, \\
& \frac{\mathbf{u}_f'^2}{2} + (\mathbf{u}_f \cdot \mathbf{n}) \mathbf{u}_f && \text{on } \Gamma_N, \\
& - (\mathbf{u}_f \cdot \mathbf{n}) \mathbf{u}_f' - (\mathbf{u}_f \cdot \mathbf{u}_f') \mathbf{n} && \text{on } \Gamma_N, \\
\end{aligned}
\right.
\quad (B.4) $$

$$ 
\left\{
\begin{aligned}
& \nabla \cdot (d(\hat{\gamma}) \nabla T') - T' + 2(1 - \hat{\gamma}) T = 0 && \text{in } \Omega \\
& T' = 0 && \text{on } \Gamma_T \\
& d \nabla T' \cdot \mathbf{n} = 0 && \text{on } \Gamma_q. \\
\end{aligned}
\right.
\quad (B.5) $$

The sensitivity of the objective function $J$ with respect to the design variable $\gamma$ can be described as follows:

$$ \frac{\delta \mathcal{J}}{\delta \hat{\gamma}} = \frac{\delta \mathcal{J}}{\delta \hat{\gamma}} \frac{\partial \hat{\gamma}}{\partial \hat{\gamma}} \frac{\delta \tilde{\gamma}}{\delta \gamma}, \quad (B.6) $$

where the functional derivative of $\mathcal{J}$ with respect to the projected density $\hat{\gamma}$ is derived as:

$$ \frac{\delta \mathcal{J}}{\delta \hat{\gamma}} = - \frac{\partial \alpha}{\partial \hat{\gamma}} (\mathbf{u}_r \cdot \mathbf{u}_r') = \alpha_{\max} (\mathbf{u}_r \cdot \mathbf{u}_r') \text{ in } \Omega, \quad (B.7) $$

and the term $\frac{\partial \hat{\gamma}}{\partial \tilde{\gamma}}$ is the derivative of the smoothed Heaviside function (6) with respect to the filtered density $\tilde{\gamma}$:

$$ \frac{\partial \hat{\gamma}}{\partial \tilde{\gamma}} = \frac{\beta (1 - \tanh^2(\beta(\tilde{\gamma} - \eta)))}{\tanh(\beta\eta) + \tanh(\beta(1 - \eta))}. \quad (B.8) $$
1174 Then, the functional derivative of $\mathcal{J}$ with respect
to the design variable $\gamma$ is obtained as the solution
1175
1176 of the adjoint Helmholtz-type partial differential
equation, that is:

$$
\begin{equation}
\begin{aligned}
& 1178 \\
& 1179 \quad \left\{ -r^2 \nabla^2 S_{\mathcal{J}} + S_{\mathcal{J}} = S_{\mathcal{J}, \hat{\gamma}} \quad \begin{array}{@{}l@{}} \text{in } \Omega \\ \text{on } \Gamma, \end{array} \right. \\
& 1180 \\
& 1181 \quad \text{where,} \\
& 1182
\end{aligned}
\tag{B.9}
\end{equation}
$$

$$
\begin{align}
& 1183 \nonumber \\
& 1184 \quad S_{\mathcal{J}} := \frac{\delta \mathcal{J}}{\delta \hat{\gamma}}, \quad S_{\mathcal{J}, \hat{\gamma}} := \frac{\delta \mathcal{J}}{\delta \hat{\gamma}} \frac{\partial \hat{\gamma}}{\partial \hat{\gamma}}. \tag{B.10} \\
& 1185 \nonumber \\
& 1186 \nonumber \\
& 1187 \nonumber
\end{align}
$$

In the same manner, the sensitivities of con-
straint functions G₁ and G₂ with respect to the
design variable γ can be derived by replacing J in
the above process with G₁ and G₂. The resultant
functional derivatives of G₁ and G₂ with respect
to the projected density ŷ are described as:

$$
\begin{align}
& 1193 \quad \frac{\delta \mathcal{G}_1}{\delta \hat{\gamma}} = \frac{\partial \alpha}{\partial \hat{\gamma}} (\mathbf{u}_f \cdot \mathbf{u}_f') = -\alpha_{\max} (\mathbf{u}_f \cdot \mathbf{u}_f') \text{ in } \Omega, \notag \\
& 1194 \notag \\
& 1195 \quad \frac{\delta \mathcal{G}_2}{\delta \hat{\gamma}} = -\frac{\partial d}{\partial \hat{\gamma}} (\nabla T \cdot \nabla T') - T' - T^2 \notag \\
& 1196 \notag \\
& 1197 \quad = (d_{\max} - d_{\min}) (\nabla T \cdot \nabla T') - T' - T^2 \text{ in } \Omega. \tag{B.11} \\
& 1198 \notag \\
& 1199 \notag \\
& 1200 \notag
\end{align}
$$

1200
1201 Appendix C Fictitious
1202
1203 physi-
1204
1205 cal model
1206
1207 formulation

In the previous work [Yamada and Noguchi, 2022],
the governing equation for the fictitious physical
field is given as follows:

$$
\begin{equation}
\begin{split}
& 1211 \quad \left\{ \begin{array}{@{}l@{}} \nabla \cdot (d(\hat{\gamma}) \nabla T) + (1-\hat{\gamma})(1-T) = 0 \\ T=0 \\ d\nabla T \cdot \mathbf{n}=0 \end{array} \right. \quad \begin{array}{@{}l@{}} \text{in } \Omega \\ \text{on } \Gamma_T \\ \text{on } \Gamma_q, \end{array} \\
& 1212 \quad \text{(C.1)}
\end{split}
\end{equation}
$$

The difference is in the definition of the heat
source,

$$
\begin{equation}
\begin{split}
& 1218 \\
& 1219 \quad Q(\hat{\gamma}, T) = \begin{cases} (1 - \hat{\gamma} - T) := Q_1(\hat{\gamma}, T) \\ (1 - \hat{\gamma})(1 - T) := Q_2(\hat{\gamma}, T), \end{cases} \tag{C.2}
\end{split}
\end{equation}
$$

1221
1222 where $Q_1(\hat{\gamma}, T)$ is used in this work and $Q_2(\hat{\gamma}, T)$
1223 is used in the previous work.
1224

Fig. C.1 shows the fictitious physical field when the governing equation (C.1) is used for the same problem setting in Section 4.1. In this case, it seems that $d_{\max} = 1 \times 10^1$ and $d_{\min} \le 1 \times 10^{-5}$ are the appropriate values for the diffusion coefficients. As mentioned in [Yamada and Noguchi, 2022], $T \approx 1$ in the “floating” islands (or, closed cavities in [Yamada and Noguchi, 2022]) is achieved regardless of their size.

However, in the fluid domain (or the mate-
rial domain) ŷ = 1, the fictitious temperature T
has non-zero value especially near the “floating”
islands, which may lead to undesirable effects on
the optimization process. This behavior may have
not become an important issue in the previous
work because the target of penalization was the
void cavities in the mean compliance minimiza-
tion problem and the thermal diffusivity problem,
where the material domains ŷ = 1 are usually sur-
rounded by the void domains ŷ = 0. In the fluid
problems, however, it is often the case that the
fluid domains occupy most of the design domain,
especially in the case without volume constraint.
Therefore, the fictitious physical field governed
by (C.1) might not work as expected, and we
employed another expression Q₁(ŷ, T) in (C.2) for
the heat source Q(ŷ, T) in this paper.

**Supplementary Information**

No supplementary information is used in this
paper.

Acknowledgement

No funding was received for conducting this study.

Declarations

Conflict of interest

On behalf of all authors, the corresponding author
states that there is no conflict of interest.

Replication of results

The presented results were obtained using Open-
FOAM and we can provide them upon reasonable
request.
1174 Then, the functional derivative of $\mathcal{J}$ with respect to the design variable $\gamma$ is obtained as the solution of the adjoint Helmholtz-type partial differential equation, that is:

$$ 
\begin{array}{ll}
1178 & -r^2 \nabla^2 S_{\mathcal{J}} + S_{\mathcal{J}} = S_{\mathcal{J}, \hat{\gamma}} \quad \text{in } \Omega \\
1179 & \nabla S_{\mathcal{J}} \cdot \mathbf{n} = 0 \quad \text{on } \Gamma, \\
1180 & \\
1181 & \text{where,} \\
1182 & \\
1183 & S_{\mathcal{J}} := \frac{\delta \mathcal{J}}{\delta \hat{\gamma}}, \quad S_{\mathcal{J}, \hat{\gamma}} := \frac{\delta \mathcal{J}}{\delta \hat{\gamma}} \frac{\partial \hat{\gamma}}{\partial \hat{\gamma}}. \qquad (\text{B.10})
\end{array}
$$

In the same manner, the sensitivities of constraint functions $G_1$ and $G_2$ with respect to the design variable $\gamma$ can be derived by replacing $\mathcal{J}$ in the above process with $G_1$ and $G_2$. The resultant functional derivatives of $G_1$ and $G_2$ with respect to the projected density $\hat{\gamma}$ are described as:

$$ 
\begin{array}{ll}
1193 & \frac{\delta G_1}{\delta \hat{\gamma}} = \frac{\partial \alpha}{\partial \hat{\gamma}} (\mathbf{u}_f \cdot \mathbf{u}_{f'}) = -\alpha_{\max} (\mathbf{u}_f \cdot \mathbf{u}_{f'}) \quad \text{in } \Omega, \\
1194 & \\
1195 & \frac{\delta G_2}{\delta \hat{\gamma}} = -\frac{\partial d}{\partial \hat{\gamma}} (\nabla T \cdot \nabla T') - T' - T^2 \\
1196 & \\
1197 & = (d_{\max} - d_{\min}) (\nabla T \cdot \nabla T') - T' - T^2 \quad \text{in } \Omega. \\
1198 & \\
1199 & \\
1200 & \\
1201 & \text{Appendix C Fictitious} \\
1202 & \text{physical model} \\
1203 & \text{formulation}
\end{array}
$$

## Appendix C Fictitious physical model formulation

In the previous work [Yamada and Noguchi, 2022], the governing equation for the fictitious physical field is given as follows:

$$ 
\begin{array}{ll}
1211 & \nabla \cdot (d(\hat{\gamma}) \nabla T) + (1-\hat{\gamma})(1-T) = 0 \quad \text{in } \Omega \\
1212 & T = 0 \quad \text{on } \Gamma_T \\
1213 & d\nabla T \cdot \mathbf{n} = 0 \quad \text{on } \Gamma_q, \\
1214 & \\
1215 & (\text{C.1})
\end{array}
$$

The difference is in the definition of the heat source,

$$ 
\begin{array}{ll}
1218 & Q(\hat{\gamma}, T) = \left\{ \begin{array}{ll} (1-\hat{\gamma}-T) := Q_1(\hat{\gamma}, T) & \text{on } \Gamma_T \\ (1-\hat{\gamma})(1-T) := Q_2(\hat{\gamma}, T), & \text{on } \Gamma_q \end{array} \right. \\
1219 & \\
1220 & \\
1221 & \\
1222 & \text{where } Q_1(\hat{\gamma}, T) \text{ is used in this work and } Q_2(\hat{\gamma}, T) \text{ is used in the previous work.}
\end{array}
$$

1223
1224

Fig. C.1 shows the fictitious physical field when the governing equation (C.1) is used for the same problem setting in Section 4.1. In this case, it seems that $d_{\max} = 1 \times 10^1$ and $d_{\min} \le 1 \times 10^{-5}$ are the appropriate values for the diffusion coefficients. As mentioned in [Yamada and Noguchi, 2022], $T \approx 1$ in the “floating” islands (or, closed cavities in [Yamada and Noguchi, 2022]) is achieved regardless of their size.

However, in the fluid domain (or the material domain) $\hat{\gamma} = 1$, the fictitious temperature $T$ has non-zero value especially near the “floating” islands, which may lead to undesirable effects on the optimization process. This behavior may have not become an important issue in the previous work because the target of penalization was the void cavities in the mean compliance minimization problem and the thermal diffusivity problem, where the material domains $\hat{\gamma} = 1$ are usually surrounded by the void domains $\hat{\gamma} = 0$. In the fluid problems, however, it is often the case that the fluid domains occupy most of the design domain, especially in the case without volume constraint. Therefore, the fictitious physical field governed by (C.1) might not work as expected, and we employed another expression $Q_1(\hat{\gamma}, T)$ in (C.2) for the heat source $Q(\hat{\gamma}, T)$ in this paper.

## Supplementary Information

No supplementary information is used in this paper.

## Acknowledgement

No funding was received for conducting this study.

## Declarations

### Conflict of interest

On behalf of all authors, the corresponding author states that there is no conflict of interest.

### Replication of results

The presented results were obtained using OpenFOAM and we can provide them upon reasonable request.
<figure><img src="image_30.png" /><figcaption>**Fig. C.1:** Fictitious physical field *T* in the work by Yamada and Noguchi [2022] for various diffusion coefficients (*d*<sub>max</sub>, *d*<sub>min</sub>)</figcaption></figure>

References

G. Allaire, F. Jouve, and G. Michailidis. Molding direction constraints in structural optimization via a level-set method. *Variational Analysis and Aerospace Engineering: Mathematical Challenges for the Aerospace of the Future*, pages 1–39, 2016.

M. P. Bendsøe and N. Kikuchi. Generating optimal topologies in structural design using a homogenization method. *Computer Methods in Applied Mechanics and Engineering*, 71(2): 197–224, 1988. ISSN 0045-7825.

T. Borrvall and J. Petersson. Topology optimization of fluids in stokes flow. *International Journal for Numerical Methods in Fluids*, 41(1): 77–107, 2003.

Y. Deng, Z. Liu, P. Zhang, Y. Wu, and J. G. Korvink. Optimization of no-moving part fluidic resistance microvalves with low reynolds number. In *2010 IEEE 23rd International conference on micro electro mechanical systems (MEMS)*, pages 67–70. IEEE, 2010.

Y. Deng, Z. Liu, P. Zhang, Y. Liu, and Y. Wu. Topology optimization of unsteady incompressible navier–stokes flows. *Journal of Computational Physics*, 230(17):6688–6708, 2011.

C. B. Dilgen, S. B. Dilgen, D. R. Fuhrman, O. Sigmund, and B. S. Lazarov. Topology optimization of turbulent flows. *Computer Methods in Fluids*

1225

1226

1227

1228

1229

1230

1231

1232

1233

1234

1235

1236

1237

1238

1239

1240

1241

1242

1243

1244

1245

1246

1247

1248

1249

1250

1251

1252

1253

1254

1255

1256

1257

1258

1259

1260

1261

1262

1263

1264

1265

1266

1267

1268

1269

1270

1271

1272

1273

1274

1275
<figure><img src="image_33.png" /><figcaption>**Fig. C.1:** Fictitious physical field *T* in the work by Yamada and Noguchi [2022] for various diffusion coefficients (*d*<sub>max</sub>, *d*<sub>min</sub>)</figcaption></figure>

References

G. Allaire, F. Jouve, and G. Michailidis. Molding direction constraints in structural optimization via a level-set method. *Variational Analysis and Aerospace Engineering: Mathematical Challenges for the Aerospace of the Future*, pages 1–39, 2016.

M. P. Bendsøe and N. Kikuchi. Generating optimal topologies in structural design using a homogenization method. *Computer Methods in Applied Mechanics and Engineering*, 71(2): 197–224, 1988. ISSN 0045-7825.

T. Borrvall and J. Petersson. Topology optimization of fluids in stokes flow. *International Journal for Numerical Methods in Fluids*, 41(1): 77–107, 2003.

Y. Deng, Z. Liu, P. Zhang, Y. Wu, and J. G. Korvink. Optimization of no-moving part fluidic resistance microvalves with low reynolds number. In *2010 IEEE 23rd International conference on micro electro mechanical systems (MEMS)*, pages 67–70. IEEE, 2010.

Y. Deng, Z. Liu, P. Zhang, Y. Liu, and Y. Wu. Topology optimization of unsteady incompressible navier–stokes flows. *Journal of Computational Physics*, 230(17):6688–6708, 2011.

C. B. Dilgen, S. B. Dilgen, D. R. Fuhrman, O. Sigmund, and B. S. Lazarov. Topology optimization of turbulent flows. *Computer Methods in Fluids*

1225

1226

1227

1228

1229

1230

1231

1232

1233

1234

1235

1236

1237

1238

1239

1240

1241

1242

1243

1244

1245

1246

1247

1248

1249

1250

1251

1252

1253

1254

1255

1256

1257

1258

1259

1260

1261

1262

1263

1264

1265

1266

1267

1268

1269

1270

1271

1272

1273

1274

1275
1276 *in Applied Mechanics and Engineering*, 331: 363-393, 2018.

1277 F. K. Forster, R. L. Bardell, M. A. Afromowitz,

1279 N. R. Sharma, and A. Blanchard. Design, fabrication and testing of fixed-valve micro-pumps.

1280 *Asme-Publications-Fed*, 234:39-44, 1995.

1282 N. Galanos, E. M. Papoutsis-Kiachagias, K. C.

1283 Giannakoglou, Y. Kondo, and K. Tanimoto.

1284 Synergistic use of adjoint-based topology and shape optimization for the design of bi-fluid heat exchangers. *Structural and Multidisciplinary Optimization*, 65(9):245, 2022.

1288 A. R. Gamboa, C. J. Morris, and F. K. Forster.

1289 Improvements in fixed-valve micropump performance through shape optimization of valves. *J. Fluids Eng.*, 127(2):339-346, 2005.

1291

1292 A. Gersborg-Hansen, O. Sigmund, and R. B.

1293 Haber. Topology optimization of channel flow problems. *Structural and Multidisciplinary Optimization*, 30(3):181-192, 2005.

1295 *Structural and Multidisciplinary Optimization*, 30(3):181-192, 2005.

1296 A. Ghasemi and A. Elham. Efficient multi-stage aerodynamic topology optimization using an operator-based analytical differentiation. *Structural and Multidisciplinary Optimization*, 65(4): 130-147, 2022.

1297 aerodynamic topology optimization using an

1298 operator-based analytical differentiation. *Structural and Multidisciplinary Optimization*, 65(8):1-15, 2022b.

1299 *Structural and Multidisciplinary Optimization*, 65(4): 130-147, 2022.

1300 130, 2022.

1301 J. K. Guest. Imposing maximum length scale in topology optimization. *Structural and Multidisciplinary Optimization*, 37:463-473, 2009.

1302 *Structural and Multidisciplinary Optimization*, 37:463-473, 2009.

1303 *Structural and Multidisciplinary Optimization*, 37:463-473, 2009.

1304 J. K. Guest, J. H. Prévost, and T. Belytschko.

1305 Achieving minimum length scale in topology optimization using nodal design variables and projection functions. *International Journal for Numerical Methods in Engineering*, 61(2): 238-254, 2004.

1306 projection functions. *International Journal for Numerical Methods in Engineering*, 61(2): 238-254, 2004.

1307 *International Journal for Numerical Methods in Engineering*, 61(2): 238-254, 2004.

1308 *International Journal for Numerical Methods in Engineering*, 61(2): 238-254, 2004.

1309 238-254, 2004.

1310 Y. Guo, H. Pan, E. Wadbro, and Z. Liu. Design applicable 3d microfluidic functional units using 2d topology optimization with length scale constraints. *Micromachines*, 11(6):613, 2020.

1311 applicable 3d microfluidic functional units using

1312 2d topology optimization with length scale con-

1313 constraints. *Micromachines*, 11(6):613, 2020.

1314 A. Kawamoto, T. Matsumori, S. Yamasaki, T. Nomura, T. Kondoh, and S. Nishiwaki.

1315 Heaviside projection based topology optimization by a pde-filtered scalar function. *Structural and Multidisciplinary Optimization*, 44(1): 19-24, 2011.

1316 Heaviside projection based topology optimization by a pde-filtered scalar function. *Structural and Multidisciplinary Optimization*, 44(1): 19-24, 2011.

1317 *Structural and Multidisciplinary Optimization*, 44(1): 19-24, 2011.

1318 *Structural and Multidisciplinary Optimization*, 44(1): 19-24, 2011.

1319 19-24, 2011.

1320 E. Kontoleontos, E. Papoutsis-Kiachagias, A. Zymaris, D. Papadimitriou, and K. Giannakoglou. Adjoint-based constrained topology optimization for viscous flows, including heat transfer. *Engineering Optimization*, 45(8): 941-961, 2013.

1321 A. Zymaris, D. Papadimitriou, and K. Giannakoglou. Adjoint-based constrained topology optimization for viscous flows, including heat transfer. *Engineering Optimization*, 45(8): 941-961, 2013.

1322 optimization for viscous flows, including heat transfer. *Engineering Optimization*, 45(8): 941-961, 2013.

1323 transfer. *Engineering Optimization*, 45(8): 941-961, 2013.

1324 transfer. *Engineering Optimization*, 45(8): 941-961, 2013.

1325 941-961, 2013.

1326

S. Kreissl, G. Pingen, and K. Maute. Topology optimization for unsteady flow. *International Journal for Numerical Methods in Engineering*, 87(13):1229-1253, 2011.

B. S. Lazarov and O. Sigmund. Filters in topology optimization based on helmholtz-type differential equations. *International Journal for Numerical Methods in Engineering*, 86(6): 765-781, 2011.

H. Li, T. Kondoh, P. Jolivet, K. Furuta, T. Yamada, B. Zhu, H. Zhang, K. Izui, and S. Nishiwaki. Optimum design and thermal modeling for 2d and 3d natural convection problems incorporating level set-based topology optimization with body-fitted mesh. *International Journal for Numerical Methods in Engineering*, 123(9):1954-1990, 2022a.

H. Li, T. Kondoh, P. Jolivet, N. Nakayama, K. Furuta, H. Zhang, B. Zhu, K. Izui, and S. Nishiwaki. Topology optimization for lift-drag problems incorporated with distributed unstructured mesh adaptation. *Structural and Multidisciplinary Optimization*, 65(8):1-15, 2022b.

H. Li, M. Yu, P. Jolivet, J. Alexandersen, T. Kondoh, T. Hu, K. Furuta, K. Izui, and S. Nishiwaki. Reaction-diffusion equation driven topology optimization of high-resolution and feature-rich structures using unstructured meshes. *Advances in Engineering Software*, 180:103457, 2023.

Q. Li, W. Chen, S. Liu, and L. Tong. Structural topology optimization considering connectivity constraint. *Structural and Multidisciplinary Optimization*, 54:971-984, 2016.

S. Lin, L. Zhao, J. K. Guest, T. P. Weihs, and Z. Liu. Topology optimization of fixed-geometry fluid diodes. *Journal of Mechanical Design*, 137(8):081402, 2015.

S. Liu, Q. Li, W. Chen, L. Tong, and G. Cheng. An identification method for enclosed voids restriction in manufacturability design for additive manufacturing structures. *Frontiers of Mechanical Engineering*, 10:126-137, 2015.

E. Morganti, G. Pignatel, et al. Microfluidics for the treatment of the hydrocephalus. In *Proceedings of International Conference on Sensing Technology, Palmerston North*, pages 483-487, 2005.
1276 *in Applied Mechanics and Engineering*, 331: 363-393, 2018.

1277 F. K. Forster, R. L. Bardell, M. A. Afromowitz,

1279 N. R. Sharma, and A. Blanchard. Design, fabrication and testing of fixed-valve micro-pumps.

1280 *Asme-Publications-Fed*, 234:39-44, 1995.

1282 N. Galanos, E. M. Papoutsis-Kiachagias, K. C.

1283 Giannakoglou, Y. Kondo, and K. Tanimoto.

1284 Synergistic use of adjoint-based topology and shape optimization for the design of bi-fluid heat exchangers. *Structural and Multidisciplinary Optimization*, 65(9):245, 2022.

1288 A. R. Gamboa, C. J. Morris, and F. K. Forster.

1289 Improvements in fixed-valve micropump performance through shape optimization of valves. *J. Fluids Eng.*, 127(2):339-346, 2005.

1291

1292 A. Gersborg-Hansen, O. Sigmund, and R. B.

1293 Haber. Topology optimization of channel flow problems. *Structural and Multidisciplinary Optimization*, 30(3):181-192, 2005.

1294

1295 *Structural and Multidisciplinary Optimization*, 30(3):181-192, 2005.

1296 A. Ghasemi and A. Elham. Efficient multi-stage aerodynamic topology optimization using an operator-based analytical differentiation. *Structural and Multidisciplinary Optimization*, 65(4): 130-147, 2022.

1297

1298 operator-based analytical differentiation. *Structural and Multidisciplinary Optimization*, 65(8):1-15, 2022b.

1299 *Structural and Multidisciplinary Optimization*, 65(4): 130-147, 2022.

1300

1301 J. K. Guest. Imposing maximum length scale in topology optimization. *Structural and Multidisciplinary Optimization*, 37:463-473, 2009.

1302

1303 *Structural and Multidisciplinary Optimization*, 37:463-473, 2009.

1304 J. K. Guest, J. H. Prévost, and T. Belytschko.

1305 Achieving minimum length scale in topology optimization using nodal design variables and projection functions. *International Journal for Numerical Methods in Engineering*, 61(2): 238-254, 2004.

1306

1307 projection functions. *International Journal for Numerical Methods in Engineering*, 61(2): 238-254, 2004.

1308

1309

1310 Y. Guo, H. Pan, E. Wadbro, and Z. Liu. Design applicable 3d microfluidic functional units using 2d topology optimization with length scale constraints. *Micromachines*, 11(6):613, 2020.

1311

1312 2d topology optimization with length scale constraints. *Micromachines*, 11(6):613, 2020.

1314 A. Kawamoto, T. Matsumori, S. Yamasaki, T. Nomura, T. Kondoh, and S. Nishiwaki.

1315

1316 Heaviside projection based topology optimization by a pde-filtered scalar function. *Structural and Multidisciplinary Optimization*, 44(1): 19-24, 2011.

1317

1318 *Structural and Multidisciplinary Optimization*, 44(1): 19-24, 2011.

1319

1320 E. Kontoleontos, E. Papoutsis-Kiachagias, A. Zymaris, D. Papadimitriou, and K. Giannakoglou. Adjoint-based constrained topology optimization for viscous flows, including heat transfer. *Engineering Optimization*, 45(8): 941-961, 2013.

1321

1322 optimization for viscous flows, including heat transfer. *Engineering Optimization*, 45(8): 941-961, 2013.

1323

1324 transfer. *Engineering Optimization*, 45(8): 941-961, 2013.

1325

1326

S. Kreissl, G. Pingen, and K. Maute. Topology optimization for unsteady flow. *International Journal for Numerical Methods in Engineering*, 87(13):1229-1253, 2011.

B. S. Lazarov and O. Sigmund. Filters in topology optimization based on helmholtz-type differential equations. *International Journal for Numerical Methods in Engineering*, 86(6): 765-781, 2011.

H. Li, T. Kondoh, P. Jolivet, K. Furuta, T. Yamada, B. Zhu, H. Zhang, K. Izui, and S. Nishiwaki. Optimum design and thermal modeling for 2d and 3d natural convection problems incorporating level set-based topology optimization with body-fitted mesh. *International Journal for Numerical Methods in Engineering*, 123(9):1954-1990, 2022a.

H. Li, T. Kondoh, P. Jolivet, N. Nakayama, K. Furuta, H. Zhang, B. Zhu, K. Izui, and S. Nishiwaki. Topology optimization for lift-drag problems incorporated with distributed unstructured mesh adaptation. *Structural and Multidisciplinary Optimization*, 65(8):1-15, 2022b.

H. Li, M. Yu, P. Jolivet, J. Alexandersen, T. Kondoh, T. Hu, K. Furuta, K. Izui, and S. Nishiwaki. Reaction-diffusion equation driven topology optimization of high-resolution and feature-rich structures using unstructured meshes. *Advances in Engineering Software*, 180:103457, 2023.

Q. Li, W. Chen, S. Liu, and L. Tong. Structural topology optimization considering connectivity constraint. *Structural and Multidisciplinary Optimization*, 54:971-984, 2016.

S. Lin, L. Zhao, J. K. Guest, T. P. Weihs, and Z. Liu. Topology optimization of fixed-geometry fluid diodes. *Journal of Mechanical Design*, 137(8):081402, 2015.

S. Liu, Q. Li, W. Chen, L. Tong, and G. Cheng. An identification method for enclosed voids restriction in manufacturability design for additive manufacturing structures. *Frontiers of Mechanical Engineering*, 10:126-137, 2015.

E. Morganti, G. Pignatel, et al. Microfluidics for the treatment of the hydrocephalus. In *Proceedings of International Conference on Sensing Technology, Palmerston North*, pages 483-487, 2005.
<table>
  <tbody>
    <tr>
      <td>M. Nabavi. Steady and unsteady flow analysis in microdiffusers and micropumps: a critical review. <i>Microfluidics and Nanofluidics</i>, 7: 599–619, 2009.</td>
      <td>1327</td>
    </tr>
    <tr>
      <td>L. H. Olesen, F. Okkels, and H. Bruus. A high-level programming-language implementation of topology optimization applied to steady-state Navier-Stokes flow. <i>International Journal for Numerical Methods in Engineering</i>, 65(7): 975–1001, 2006.</td>
      <td>1328<br>1329<br>1330<br>1331<br>1332</td>
    </tr>
    <tr>
      <td>C. Othmer. A continuous adjoint formulation for the computation of topological and surface sensitivities of ducted flows. <i>International Journal for Numerical Methods in Fluids</i>, 58(8): 861–877, 2008.</td>
      <td>1333<br>1334<br>1335<br>1336</td>
    </tr>
    <tr>
      <td>S. Pan, M. Yu, H. Li, Z. Li, M. Ren, J. Gu, and C. Shen. An integrated two-step strategy for an optimal design of liquid-cooled channel layout based on the mmc-density approach. <i>Structural and Multidisciplinary Optimization</i>, 65(8):221, 2022.</td>
      <td>1337<br>1338<br>1339<br>1340<br>1341</td>
    </tr>
    <tr>
      <td>G. Pingen and K. Maute. Optimal design for non-newtonian flows using a topology optimization approach. <i>Computers &amp; Mathematics with Applications</i>, 59(7):2340–2350, 2010.</td>
      <td>1342<br>1343<br>1344<br>1345<br>1346</td>
    </tr>
    <tr>
      <td>T. A. Poulsen. A new scheme for imposing a minimum length scale in topology optimization. <i>International Journal for Numerical Methods in Engineering</i>, 57(6):741–760, 2003.</td>
      <td>1347<br>1348<br>1349<br>1350</td>
    </tr>
    <tr>
      <td>J. Romero and E. C. N. Silva. A topology optimization approach applied to laminar flow machine rotor design. <i>Computer Methods in Applied Mechanics and Engineering</i>, 279: 268–300, 2014.</td>
      <td>1351<br>1352<br>1353<br>1354<br>1355</td>
    </tr>
    <tr>
      <td>Y. Sato, K. Yaji, K. Izui, T. Yamada, and S. Nishiwaki. Topology optimization of a non-moving-part valve incorporating pareto frontier exploration. <i>Structural and Multidisciplinary Optimization</i>, 56:839–851, 2017a.</td>
      <td>1356<br>1357<br>1358<br>1359<br>1360</td>
    </tr>
    <tr>
      <td>Y. Sato, T. Yamada, K. Izui, and S. Nishiwaki. Manufacturability evaluation for molded parts using fictitious physical models, and its application in topology optimization. <i>The International Journal of Advanced Manufacturing Technology</i>, 92:1391–1409, 2017b.</td>
      <td>1361<br>1362<br>1363<br>1364<br>1365</td>
    </tr>
    <tr>
      <td>K. Svanberg. The method of moving asymptotes—a new method for structural optimization. <i>International Journal for Numerical Methods in Engineering</i>, 24(2):359–373, 1987.</td>
      <td>1366<br>1367<br>1368<br>1369<br>1370<br>1371</td>
    </tr>
    <tr>
      <td>T. Truong and N. Nguyen. Simulation and optimization of tesla valves. <i>Nanotech</i>, 1:178–181, 2003.</td>
      <td>1372<br>1373<br>1374<br>1375<br>1376<br>1377</td>
    </tr>
  </tbody>
</table>
<table>
  <tbody>
    <tr>
      <td>M. Nabavi. Steady and unsteady flow analysis in microdiffusers and micropumps: a critical review. <i>Microfluidics and Nanofluidics</i>, 7: 599–619, 2009.</td>
      <td>1327</td>
    </tr>
    <tr>
      <td>L. H. Olesen, F. Okkels, and H. Bruus. A high-level programming-language implementation of topology optimization applied to steady-state Navier-Stokes flow. <i>International Journal for Numerical Methods in Engineering</i>, 65(7): 975–1001, 2006.</td>
      <td>1328<br>1329<br>1330<br>1331<br>1332</td>
    </tr>
    <tr>
      <td>C. Othmer. A continuous adjoint formulation for the computation of topological and surface sensitivities of ducted flows. <i>International Journal for Numerical Methods in Fluids</i>, 58(8): 861–877, 2008.</td>
      <td>1333<br>1334<br>1335<br>1336</td>
    </tr>
    <tr>
      <td>S. Pan, M. Yu, H. Li, Z. Li, M. Ren, J. Gu, and C. Shen. An integrated two-step strategy for an optimal design of liquid-cooled channel layout based on the mmc-density approach. <i>Structural and Multidisciplinary Optimization</i>, 65(8):221, 2022.</td>
      <td>1337<br>1338<br>1339<br>1340<br>1341</td>
    </tr>
    <tr>
      <td>G. Pingen and K. Maute. Optimal design for non-newtonian flows using a topology optimization approach. <i>Computers &amp; Mathematics with Applications</i>, 59(7):2340–2350, 2010.</td>
      <td>1342<br>1343<br>1344<br>1345<br>1346</td>
    </tr>
    <tr>
      <td>T. A. Poulsen. A new scheme for imposing a minimum length scale in topology optimization. <i>International Journal for Numerical Methods in Engineering</i>, 57(6):741–760, 2003.</td>
      <td>1347<br>1348<br>1349<br>1350</td>
    </tr>
    <tr>
      <td>J. Romero and E. C. N. Silva. A topology optimization approach applied to laminar flow machine rotor design. <i>Computer Methods in Applied Mechanics and Engineering</i>, 279: 268–300, 2014.</td>
      <td>1351<br>1352<br>1353<br>1354<br>1355</td>
    </tr>
    <tr>
      <td>Y. Sato, K. Yaji, K. Izui, T. Yamada, and S. Nishiwaki. Topology optimization of a non-moving-part valve incorporating pareto frontier exploration. <i>Structural and Multidisciplinary Optimization</i>, 56:839–851, 2017a.</td>
      <td>1356<br>1357<br>1358<br>1359<br>1360</td>
    </tr>
    <tr>
      <td>Y. Sato, T. Yamada, K. Izui, and S. Nishiwaki. Manufacturability evaluation for molded parts using fictitious physical models, and its application in topology optimization. <i>The International Journal of Advanced Manufacturing Technology</i>, 92:1391–1409, 2017b.</td>
      <td>1361<br>1362<br>1363<br>1364<br>1365</td>
    </tr>
    <tr>
      <td>K. Svanberg. The method of moving asymptotes—a new method for structural optimization. <i>International Journal for Numerical Methods in Engineering</i>, 24(2):359–373, 1987.</td>
      <td>1366<br>1367<br>1368<br>1369<br>1370<br>1371</td>
    </tr>
    <tr>
      <td>T. Truong and N. Nguyen. Simulation and optimization of tesla valves. <i>Nanotech</i>, 1:178–181, 2003.</td>
      <td>1372<br>1373<br>1374<br>1375<br>1376<br>1377</td>
    </tr>
  </tbody>
</table>