Topology optimization on two-dimensional manifolds

TL;DR

This paper extends topology optimization techniques to two-dimensional manifolds, enabling design on complex surfaces like spheres and tori for applications in fluid mechanics, heat transfer, and electromagnetics.

Contribution

It introduces a method for implementing topology optimization on 2D manifolds, including material interpolation and boundary condition penalization on complex surfaces.

Findings

Successful numerical tests on sphere, torus, Möbius strip, Klein bottle

Applications demonstrated in fluid mechanics, heat transfer, electromagnetics

Method accommodates complex 2D geometries in optimization

Abstract

This paper implements topology optimization on two-dimensional manifolds. In this paper, the material interpolation is implemented on a material parameter in the partial differential equation used to describe a physical field, when this physical field is defined on a two-dimensional manifold; the material density is used to formulate a mixed boundary condition of the physical field and implement the penalization between two different types of boundary conditions, when this physical field is defined on a three-dimensional domain with its boundary conditions defined on the two-dimensional manifold corresponding a surface or an interface of this three-dimensional domain. Based on the homeomorphic property of two-dimensional manifolds, typical two-dimensional manifolds, e.g., sphere, torus, M\"{o}bius strip and Klein bottle, are included in the numerical tests, which are provided for the problems on fluidic mechanics, heat transfer and electromagnetics.

Figures

Tables

Table 1. Table 1: Pseudocode for the iterative solution of topology optimization on 2-manifolds. In the loop, n s ​ u ​ b superscript 𝑛 𝑠 𝑢 𝑏 n^{sub} is the loop-index, n m ​ a ​ x s ​ u ​ b subscript superscript 𝑛 𝑠 𝑢 𝑏 𝑚 𝑎 𝑥 n^{sub}_{max} is the maximal value of n s ​ u ​ b superscript 𝑛 𝑠 𝑢 𝑏 n^{sub} , n u ​ p ​ t superscript 𝑛 𝑢 𝑝 𝑡 n^{upt} is the updating interval of the projection parameter β 𝛽 \beta , and J n s ​ u ​ b subscript 𝐽 superscript 𝑛 𝑠 𝑢 𝑏 J_{n^{sub}} is the objective value corresponding to the loop-index n s ​ u ​ b superscript 𝑛 𝑠 𝑢 𝑏 n^{sub} .

Choose $r$, $V_f$ and $n_{max}^{sub}$;

Set $i\leftarrow 1$, $\gamma \leftarrow V_f$, $n^{sub}\leftarrow 1$, $\xi \leftarrow 0.5$ and $\beta \leftarrow 1$;

loop

Derive $\overline{\gamma }$ by filtering and projecting $\gamma$, and evaluate $V$;

Solve $u$ from the physical PDE, evaluate $J$, and set $J_{n^{sub}}\leftarrow J$;

Solve $\mu$ and $\nu$ from the adjoint equations for $J$, and evaluate $J^{\prime }$;

Solve $\nu$ from the adjoint equation for $V$, and evaluate $V^{\prime }$;

Update $\gamma$ based on $J^{\prime }$ and $V^{\prime }$;

if$mod(n^{sub},n^{upt})==0$

$\beta \leftarrow 2\beta$;

end(if)

if $(n^{sub}==n_{max}^{sub})$ or

$(\beta ==2^{10},\frac{1}{5}{\sum }_{m=0}^4|J_{n^{sub}}-J_{n^{sub}-m}|/J_0\le {10}^{-3},\mathrm{and}|V-V_f|\le {10}^{-3})$

break;

end(if)

$n^{sub}\leftarrow n^{sub}+1$

end(loop)

Table 2. Table 2: Parameters used to implement topology optimization of the micro-textures for the wetting behaviors in the Cassie-Baxter mode.

$L$	${\overline{\sigma }}_l$	$z_0$	$V_f$	$n_{max}^{sub}$	$n^{upt}$
$10\mu m$	${10}^2$	$1\mu m$	$0.2$	$315$	$30$

Table 3. (a) Sphere

${\overline{\sigma }}_l={10}^2$	$2.97\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟔}}$	$1.27\times {10}^{-5}$
${\overline{\sigma }}_l={10}^3$	$9.43\times {10}^{-8}$	$4.55\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟖}}$

Table 4. (a) Sphere

${\overline{\sigma }}_l={10}^2$	$2.97\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟔}}$	$1.27\times {10}^{-5}$
${\overline{\sigma }}_l={10}^3$	$9.43\times {10}^{-8}$	$4.55\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟖}}$

Table 5. (b) Torus

${\overline{\sigma }}_l={10}^2$	$9.78\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟔}}$	$1.13\times {10}^{-4}$
${\overline{\sigma }}_l={10}^3$	$2.04\times {10}^{-6}$	$1.10\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟔}}$

Table 6. (a) Möbius strip

	Figure 10(a1)	Figure 10(b1)
$P_1$	$1.31$	$2.58$
$P_2$	$2.23$	$1.21$

Table 7. (a) Möbius strip

	Figure 10(a1)	Figure 10(b1)
$P_1$	$1.31$	$2.58$
$P_2$	$2.23$	$1.21$

Table 8. (b) Klein bottle

	Figure 11(a1)	Figure 11(b1)
$P_3$	$0.09$	$0.34$
$P_4$	$0.20$	$0.14$

Table 9. (a) Torus

	Figure 13(a)	Figure 13(b)	Figure 13(c)
$x$-polarization	$107.8$	$81.7$	$97.9$
$y$-polarization	$58.2$	$110.1$	$68.5$
Circular polarization	$54.3$	$102.1$	$125.5$

Table 10. (a) Torus

	Figure 13(a)	Figure 13(b)	Figure 13(c)
$x$-polarization	$107.8$	$81.7$	$97.9$
$y$-polarization	$58.2$	$110.1$	$68.5$
Circular polarization	$54.3$	$102.1$	$125.5$

Table 11. (b) Möbius strip

	Figure 14(a)	Figure 14(b)	Figure 14(c)
$x$-polarization	$22.8$	$12.9$	$10.1$
$y$-polarization	$18.3$	$19.9$	$8.4$
Circular polarization	$16.0$	$15.8$	$16.7$

Equations

$$\begin{array}{rl}& {\nabla }_s\cdot \left[p\left(\bar{\gamma }\right)\mathbf{g}\left({\nabla }_{s}u,u\right)\right]=c_s,\mathrm{in}{\Sigma }_S\\ & \left[p\left(\bar{\gamma }\right)\mathbf{g}\left({\nabla }_{s}u,u\right)\right]\cdot \mathbf{\tau }=c_b,\mathrm{on}\partial {\Sigma }_S\\ & u=u_0,\mathrm{at}\mathcal{P}\subset {\Sigma }_S\end{array}$$

$$\begin{array}{rl}& e\left(u;\bar{\gamma }\right):={\int }_{{\Sigma }_S}-p\left(\bar{\gamma }\right)\mathbf{g}\left({\nabla }_{s}u,u\right)\cdot {\nabla }_s\hat{u}-c_s\hat{u}\mathrm{d}s+{\int }_{\partial {\Sigma }_s}{c}_b\hat{u}\mathrm{d}l=0,\\ & \forall \hat{u}\in \mathcal{H}\left({\Sigma }_S\right)\end{array}$$

$$\begin{array}{rl}& \nabla \cdot \left[p\mathbf{g}\left(\nabla u,u\right)\right]=c_s,\mathrm{in}\Omega \\ & \left[p\mathbf{g}\left({\nabla }_{s}u,u\right)\right]\cdot \mathbf{n}=\alpha \left(\bar{\gamma }\right)\left(u-u_d\right),\mathrm{on}{\Sigma }_S=\partial \Omega \\ & u=u_0,\mathrm{at}\mathcal{P}\subset \Omega \end{array}$$

$$\begin{array}{r}e\left(u;\bar{\gamma }\right):={\int }_{\Omega }-p\mathbf{g}\left(\nabla u,u\right)\cdot \nabla \hat{u}\mathrm{d}\Omega +{\int }_{{\Sigma }_S}\alpha \left(\bar{\gamma }\right)\left(u-u_d\right)\hat{u}\mathrm{d}s=0,\forall \hat{u}\in \mathcal{H}\left(\Omega \right)\end{array}$$

$$\begin{array}{rl}& \nabla \cdot \left[p\mathbf{g}\left(\nabla u,u\right)\right]=c_s,\mathrm{in}\Omega \\ & [[p\mathbf{g}\left({\nabla }_{s}u,u\right)]]\cdot \mathbf{n}=\alpha \left(\bar{\gamma }\right)\left(u-u_d\right),\mathrm{on}{\Sigma }_S\hookrightarrow \Omega \\ & u=u_d,\mathrm{at}\mathcal{P}\subset {\Sigma }_S\\ & u=u_0,\mathrm{on}\partial \Omega \end{array}$$

$$\begin{array}{r}e\left(u;\bar{\gamma }\right):={\int }_{\Omega }-p\mathbf{g}\left(\nabla u,u\right)\cdot \nabla \hat{u}\mathrm{d}\Omega +{\int }_{{\Sigma }_S}\alpha \left(\bar{\gamma }\right)\left(u-u_d\right)\hat{u}\mathrm{d}s=0,\forall \hat{u}\in \mathcal{H}\left(\Omega \right)\end{array}$$

$$\begin{array}{rl}& {\nabla }_s\cdot \left(-r^2{\nabla }_s\tilde{\gamma }\right)+\tilde{\gamma }=\gamma ,\mathrm{on}{\Sigma }_S\\ & -r^2{\nabla }_s\tilde{\gamma }\cdot \mathbf{\tau }=0,\mathrm{at}\partial {\Sigma }_S\end{array}$$

$$\begin{array}{r}f\left(\tilde{\gamma };\gamma \right):={\int }_{{\Sigma }_S}{r}^2{\nabla }_s\tilde{\gamma }\cdot {\nabla }_s\hat{\tilde{\gamma }}+\tilde{\gamma }\hat{\tilde{\gamma }}-\gamma \hat{\tilde{\gamma }}\mathrm{d}s=0,\forall \hat{\tilde{\gamma }}\in \mathcal{H}\left({\Sigma }_S\right)\end{array}$$

$$\bar{\gamma }=\frac{\tanh \left(\beta \xi \right)+\tanh \left(\beta \left(\tilde{\gamma }-\xi \right)\right)}{\tanh \left(\beta \xi \right)+\tanh \left(\beta \left(1-\xi \right)\right)}$$

$$\begin{array}{rl}& \mathrm{find}\gamma \left({\Sigma }_S\right)\in \left[0,1\right]\mathrm{for}\left({\Sigma }_S\times \gamma \left({\Sigma }_S\right),{\Sigma }_S,pro{j}_1,\gamma \left({\Sigma }_S\right)\right),\\ & \mathrm{to}\mathrm{minimize}J\left(u;\bar{\gamma }\right)\mathrm{constrained}\mathrm{by}\\ & \{\begin{array}{rl}& e\left(u;\bar{\gamma }\right)=0(\mathrm{Physical}\mathrm{PDE})\\ & f\left(\tilde{\gamma };\gamma \right)=0(\mathrm{PDE}\mathrm{filter})\\ & \bar{\gamma }=\frac{\tanh \left(\beta \xi \right)+\tanh \left(\beta \left(\tilde{\gamma }-\xi \right)\right)}{\tanh \left(\beta \xi \right)+\tanh \left(\beta \left(1-\xi \right)\right)}(\mathrm{Threshold}\mathrm{projection})\\ & \left|\frac{1}{\left|{\Sigma }_S\right|}{\int }_{{\Sigma }_S}\bar{\gamma }\mathrm{d}s-V_f\right|\le 10^{-3}(\mathrm{Area}\mathrm{constraint})\end{array}\end{array}$$

$$u_{\bar{\gamma }}\left(\bar{\gamma }\right)=-e_u^{-1}\left(u;\bar{\gamma }\right)e_{\bar{\gamma }}\left(u;\bar{\gamma }\right).$$

$${\tilde{\gamma }}_{\gamma }\left(\gamma \right)=-f_{\tilde{\gamma }}^{-1}\left(\tilde{\gamma };\gamma \right)f_{\gamma }\left(\tilde{\gamma };\gamma \right).$$

$$\begin{array}{rl}& {⟨J^{\prime },t⟩}_{{\mathcal{L}}^2\left({\Sigma }_S\right),{\mathcal{L}}^2\left({\Sigma }_S\right)}\\ =& {⟨J_u,u_{\bar{\gamma }}{\bar{\gamma }}_{\tilde{\gamma }}{\tilde{\gamma }}_{\gamma }t⟩}_{{\mathcal{U}}^{\ast }\left(\Omega \right),\mathcal{U}\left(\Omega \right)}+{⟨J_{\bar{\gamma }},{\bar{\gamma }}_{\tilde{\gamma }}{\tilde{\gamma }}_{\gamma }t⟩}_{\mathcal{H}\left({\Sigma }_S\right),\mathcal{H}\left({\Sigma }_S\right)}\\ =& {⟨{\tilde{\gamma }}_{\gamma }^{\ast }{\bar{\gamma }}_{\tilde{\gamma }}\left(u_{\bar{\gamma }}^{\ast }{J}_u+J_{\bar{\gamma }}\right),t⟩}_{{\mathcal{L}}^2\left({\Sigma }_S\right),{\mathcal{L}}^2\left({\Sigma }_S\right)}\\ =& {⟨-f_{\gamma }^{\ast }{\left(f_{\tilde{\gamma }}^{-1}\right)}^{\ast }{\bar{\gamma }}_{\tilde{\gamma }}\left(-e_{\bar{\gamma }}^{\ast }{\left(e_u^{-1}\right)}^{\ast }{J}_u+J_{\bar{\gamma }}\right),t⟩}_{{\mathcal{L}}^2\left({\Sigma }_S\right),{\mathcal{L}}^2\left({\Sigma }_S\right)},\forall t\in {\mathcal{L}}^2\left({\Sigma }_S\right)\end{array}$$

$$e\left(u;\bar{\gamma }\right)=e^{{\Sigma }_S}\left(u;\bar{\gamma }\right)+e^{\partial {\Sigma }_S}\left(u\right).$$

$$\begin{array}{rl}\mu & :=-{\left(e_u^{-1}\right)}^{\ast }{J}_u={\left(e_u^{{\Sigma }_S\ast }+e_u^{\partial {\Sigma }_S\ast }\right)}^{-1}{J}_u,\\ \nu & :=-{\left(f_{\tilde{\gamma }}^{-1}\right)}^{\ast }{\bar{\gamma }}_{\tilde{\gamma }}{J}_{\bar{\gamma }}={\left(f_{\tilde{\gamma }}^{\ast }\right)}^{-1}{\bar{\gamma }}_{\tilde{\gamma }}{J}_{\bar{\gamma }},\end{array}$$

$$\begin{array}{rl}& {⟨\left(e_u^{{\Sigma }_S\ast }+e_u^{\partial {\Sigma }_S\ast }\right)\mu ,v⟩}_{\mathcal{H}\left({\Sigma }_S\right),\mathcal{H}\left({\Sigma }_S\right)}={⟨-J_u,v⟩}_{\mathcal{H}\left({\Sigma }_S\right),\mathcal{H}\left({\Sigma }_S\right)},\forall v\in \mathcal{H}\left({\Sigma }_S\right)\\ & {⟨f_{\tilde{\gamma }}^{\ast }\nu ,y⟩}_{\mathcal{H}\left({\Sigma }_S\right),\mathcal{H}\left({\Sigma }_S\right)}={⟨-J_{\bar{\gamma }}{\bar{\gamma }}_{\tilde{\gamma }},y⟩}_{\mathcal{H}\left({\Sigma }_S\right),\mathcal{H}\left({\Sigma }_S\right)},\forall y\in \mathcal{H}\left({\Sigma }_S\right)\end{array}$$

$$\begin{array}{r}{⟨J^{\prime },t⟩}_{{\mathcal{L}}^2\left({\Sigma }_S\right),{\mathcal{L}}^2\left({\Sigma }_S\right)}={⟨f_{\gamma }^{\ast }\nu J_{\bar{\gamma }}^{-1}{\left(e_{\bar{\gamma }}^{{\Sigma }_S}\right)}^{\ast }\mu +f_{\gamma }^{\ast }\nu ,t⟩}_{{\mathcal{L}}^2\left({\Sigma }_S\right),{\mathcal{L}}^2\left({\Sigma }_S\right)},\forall t\in {\mathcal{L}}^2\left({\Sigma }_S\right)\end{array}$$

$$e\left(u;\bar{\gamma }\right)=e^{\Omega }\left(u\right)+e^{{\Sigma }_S}\left(u;\bar{\gamma }\right).$$

$$\begin{array}{rl}& {⟨\left(e_u^{\Omega \ast }+e_u^{{\Sigma }_S\ast }\right)\mu ,v⟩}_{\mathcal{H}\left(\Omega \right),\mathcal{H}\left(\Omega \right)}={⟨-J_u,v⟩}_{\mathcal{H}\left(\Omega \right),\mathcal{H}\left(\Omega \right)},\forall v\in \mathcal{H}\left(\Omega \right)\\ & {⟨f_{\tilde{\gamma }}^{\ast }\nu ,y⟩}_{\mathcal{H}\left({\Sigma }_S\right),\mathcal{H}\left({\Sigma }_S\right)}={⟨-J_{\bar{\gamma }}{\bar{\gamma }}_{\tilde{\gamma }},y⟩}_{\mathcal{H}\left({\Sigma }_S\right),\mathcal{H}\left({\Sigma }_S\right)},\forall y\in \mathcal{H}\left({\Sigma }_S\right)\end{array}$$

$$\begin{array}{r}{⟨V^{\prime },t⟩}_{{\mathcal{L}}^2\left({\Sigma }_S\right),{\mathcal{L}}^2\left({\Sigma }_S\right)}={⟨f_{\gamma }^{\ast }\nu ,t⟩}_{{\mathcal{L}}^2\left({\Sigma }_S\right),{\mathcal{L}}^2\left({\Sigma }_S\right)},\forall t\in {\mathcal{L}}^2\left({\Sigma }_S\right)\end{array}$$

$$\begin{array}{r}{⟨f_{\tilde{\gamma }}^{\ast }\nu ,y⟩}_{\mathcal{H}\left({\Sigma }_S\right),\mathcal{H}\left({\Sigma }_S\right)}={⟨-V_{\bar{\gamma }}{\bar{\gamma }}_{\tilde{\gamma }},y⟩}_{\mathcal{H}\left({\Sigma }_S\right),\mathcal{H}\left({\Sigma }_S\right)},\forall y\in \mathcal{H}\left({\Sigma }_S\right).\end{array}$$

$$\begin{array}{rl}& {\nabla }_s\cdot \left(\bar{\sigma }\frac{{\nabla }_s\bar{z}}{\sqrt{{\left(L/z_0\right)}^2+{\left|{\nabla }_s\bar{z}\right|}^2}}\right)=1-2{\bar{\sigma }}_l\kappa ,\mathrm{on}{\Sigma }_S\\ & \bar{\sigma }\frac{{\nabla }_s\bar{z}}{\sqrt{{\left(L/z_0\right)}^2+{\left|{\nabla }_s\bar{z}\right|}^2}}\cdot \mathbf{\tau }=0,\mathrm{at}\partial {\Sigma }_S\\ & \bar{z}=0,\mathrm{at}\mathcal{P}\end{array}$$

$$J=\frac{1}{{\left|{\Sigma }_S\right|}^2}{\int }_{{\Sigma }_S}{\bar{z}}^2\mathrm{d}s$$

$$\bar{\sigma }={\bar{\sigma }}_l+\left({\bar{\sigma }}_s-{\bar{\sigma }}_l\right)\frac{1-\bar{\gamma }}{1+10^4\bar{\gamma }}$$

$${⟨J^{\prime },t⟩}_{{\mathcal{L}}^2\left({\Sigma }_S\right),{\mathcal{L}}^2\left({\Sigma }_S\right)}=-{\int }_{{\Sigma }_S}{\tilde{\gamma }}_{a}t\mathrm{d}s,\forall t\in {\mathcal{L}}^2\left({\Sigma }_S\right).$$

$$\begin{array}{rl}& {\int }_{{\Sigma }_S}\frac{2}{{\left|{\Sigma }_S\right|}^2}\bar{z}{\hat{\bar{z}}}_a-\bar{\sigma }\frac{{\nabla }_s{\bar{z}}_a\cdot {\nabla }_s{\hat{\bar{z}}}_a}{\sqrt{{\left(L/z_0\right)}^2+{\left|{\nabla }_s\bar{z}\right|}^2}}+\bar{\sigma }\frac{\left({\nabla }_s\bar{z}\cdot {\nabla }_s{\bar{z}}_a\right)\left({\nabla }_s\bar{z}\cdot {\nabla }_s{\hat{\bar{z}}}_a\right)}{{\left(\sqrt{{\left(L/z_0\right)}^2+{\left|{\nabla }_s\bar{z}\right|}^2}\right)}^3}\mathrm{d}s=0,\\ & \forall {\hat{\bar{z}}}_a\in \mathcal{H}\left({\Sigma }_S\right)\end{array}$$

$$\begin{array}{rl}& {\int }_{{\Sigma }_S}-r^2{\nabla }_s{\tilde{\gamma }}_a\cdot {\nabla }_s{\hat{\tilde{\gamma }}}_a+{\tilde{\gamma }}_a{\hat{\tilde{\gamma }}}_a-\frac{\partial \bar{\sigma }}{\partial \bar{\gamma }}\frac{\partial \bar{\gamma }}{\partial \tilde{\gamma }}\frac{{\nabla }_s\bar{z}\cdot {\nabla }_s{\bar{z}}_a}{\sqrt{{\left(L/z_0\right)}^2+{\left|{\nabla }_s\bar{z}\right|}^2}}{\hat{\tilde{\gamma }}}_a\mathrm{d}s=0,\\ & \forall {\hat{\tilde{\gamma }}}_a\in \mathcal{H}\left({\Sigma }_S\right)\end{array}$$

$$\begin{array}{r}{⟨V^{\prime },t⟩}_{{\mathcal{L}}^2\left({\Sigma }_S\right),{\mathcal{L}}^2\left({\Sigma }_S\right)}=-{\int }_{{\Sigma }_S}{\tilde{\gamma }}_{a}t\mathrm{d}s,\forall t\in {\mathcal{L}}^2\left({\Sigma }_S\right)\end{array}$$

$$\begin{array}{r}{\int }_{{\Sigma }_S}\frac{1}{\left|{\Sigma }_S\right|}\frac{\partial \bar{\gamma }}{\partial \tilde{\gamma }}{\hat{\tilde{\gamma }}}_a\mathrm{d}s+r^2{\nabla }_s{\tilde{\gamma }}_a\cdot {\nabla }_s{\hat{\tilde{\gamma }}}_a+{\tilde{\gamma }}_a{\hat{\tilde{\gamma }}}_a\mathrm{d}s=0,\forall {\hat{\tilde{\gamma }}}_a\in \mathcal{H}\left({\Sigma }_S\right).\end{array}$$

$$\begin{array}{rl}& -\nabla \cdot \left(k\nabla T\right)=Q,\mathrm{in}\Omega \\ & T=T_d,\mathrm{on}{\Gamma }_d\\ & -k\nabla T\cdot \mathbf{n}=0,\mathrm{on}{\Gamma }_n\\ & T=T_0,\mathrm{at}\mathcal{P}\end{array}$$

Full text

Keywords: Topology optimization; two-dimensional manifold (2-manifold); material distribution method; tangential gradient operator; mixed boundary condition.

1 Introduction

Topology optimization is a robust method used to determine the structural configuration, which corresponds to the material distribution in a structure [1]. In contrast to designing devices by tuning a handful of structural parameters in size and shape optimization, topology optimization utilizes the full-parameter space to design structures based on the user-desired performance, and it is more flexible and robust, because of its low dependence on initial guess and implicit expression of the structures.

Optimization of structural topology was investigated as early as 1904 for trusses by Michell [2]. Material distribution method for topology optimization was pioneered by Bendsøe and Kikuchi for elasticity [3], and this method was extended to a variety of areas, e.g. acoustics, electromagnetics, fluid dynamics and thermodynamics [4, 5, 6, 19, 7, 20, 21, 24, 22, 23, 25, 26, 29, 18, 27, 28, 10, 8, 12, 11, 13, 9, 14, 15, 16, 17]. Several other methods have also been proposed and developed for the implementation of topology optimization, e.g. the level set method [30, 31, 32, 33], the evolutionary techniques [34, 35, 36], the moving morphable components [37, 38] and the phase field method [39].

In topology optimization, structures were usually defined on a three-dimensional (3D) domain in $\mathbb{R}^{3}$ or reduced two-dimensional (2D) plane on $\mathbb{R}^{2}$. The current development of additive manufacturing, e.g., 3D-printing, has effectively enlarged the structural-design space. Implementing topology optimization on 2-manifolds can be helpful to enhance the freedom of structural design, where a 2-manifold is a topological space with its arbitrary interior point having a neighborhood homeomorphic to a sub-region of $\mathbb{R}^{2}$ and it can be used to describe a structural surface or a material interface. Related researches have been implemented for the structural design on 2-manifolds based on the conformal geometry theory [40, 41], layouts on surfaces for shell structures [51, 52, 53, 54, 55], and material interfaces for stiffness and multi-material structures [42, 43, 44, 45, 46, 47, 48, 49, 50], fluid-structure interaction [56, 57, 58], energy absorption [59], cohesion [60] and actuation [61].

In this paper, topology optimization on 2-manifolds is implemented by using the material distribution method and this implementation includes the following two cases. When a physical field is defined on a 2-manifold, the material interpolation of topology optimization is implemented on a material parameter in the partial differential equation (PDE) with the tangential gradient operator used to describe the physical field. When the physical field is defined on a 3D domain and its boundary conditions are defined on a 2-manifold of codimension one in a 3D Euclidian space, the material density representing the structural configuration is used to implement the penalization between two different types of boundary conditions and formulate a mixed boundary condition of the physical field. In the second case, the 2-manifold can be an exterior surface or an interface of the 3D domain. When the 2-manifold is an exterior surface, it can mix the Dirichlet and Neumann types of boundary conditions; when the 2-manifold is an interface, it can mix the no-jump and Dirichlet types of boundary conditions.

For the two cases of topology optimization on 2-manifolds, a monolithic description of the corresponding optimization problems is presented in Section 2, including the adjoint analysis and numerical implementation; test problems on fluidic mechanics, heat transfer and electromagnetics are implemented in Section 3. In the following, all the mathematical descriptions are performed in the Cartesian system.

2 Methodology

Details for topology optimization implemented on 2-manifolds are provided as follows.

2.1 2-manifold

According to the classification theorem [62], a 2-manifold without boundary is compact if every open cover of it has a finite subcover; its family can be exhausted by the two infinite families \{$$\mathbb{S}^{2} (sphere), $\mathbb{T}^{2}$ (torus), $\mathbb{T}^{2}\#\mathbb{T}^{2}$ (double torus), \cdots$$\} and \{$$\mathbb{P}^{2} (projective plane), $\mathbb{P}^{2}\#\mathbb{P}^{2}$ (Klein bottle), \cdots$$\}, where $\#$ denotes the connected sum of two manifolds. A 2-manifold with boundary can be derived by removing an open disk from a 2-manifold without boundary. All 2-manifolds without boundary can be derived by gluing the basic 2-manifolds with boundaries, e.g., the disk, cylinder, Möbius strip. Structural surfaces can be described as the orientable 2-manifolds, e.g., sphere, on which the normal vector can be defined globally. Structural interfaces can include both the orientable and non-orientable 2-manifolds, e.g. Möbius strip is a non-orientable 2-manifold, on which the normal vector can be defined locally instead of globally. Therefore, without loosing generality, topology optimization in this paper is implemented on the 2-manifolds homeomorphic to sphere, torus, Möbius strip, and their connected sum or glued manifolds (Figure 1).

2.2 Physical PDEs

To implement topology optimization on a 2-manifold, a design variable, which is a relaxed binary distribution, is defined on this 2-manifold to represent the structural configuration; and it is bounded in the typical interval $\left[0,1\right]$, with [math] and $1$ representing two different material phases, respectively. An optimization problem can be formulated by minimizing or maximizing a design objective used to evaluate a desired performance of the structure implicitly expressed on the 2-manifold. The physical field related with this structure can be described by using a corresponding PDE. Then, this optimization problem is a PDE constrained optimization problem. It is nonlinear and challenging to be solved directly. Therefore, the iterative solution procedure is widely utilized. To ensure the monolithic convergence of the iterative procedure, regularization based on a PDE filter and threshold projection is imposed on the design variable, with the projected design variable referred to as the material density. This iterative procedure including the PDE filter and the threshold projection can control the feature size of the structure and remove the gray regions from the derived structural pattern.

For the case with the physical field defined on a 2-manifold, the physical PDE, used to describe the distribution of the physical field, can be expressed in a typical abstract form with the material parameter interpolated by the material density:

[TABLE]

where $u$ is the physical variable; $\Sigma_{S}$ is the 2-manifold; $c_{s}$ and $c_{b}$ are known distributions defined on $\Sigma_{S}$ and $\partial\Sigma_{S}$, respectively; if $\Sigma_{S}$ is treated as the 2-manifold of codimension one in an 3D Euclidian space, $\nabla_{s}=\nabla-\left(\mathbf{n}\cdot\nabla\right)\mathbf{n}$ is the tangential gradient operator defined on $\Sigma_{S}$, with $\mathbf{n}$ denoting the unitary normal vector of $\Sigma_{S}$; $\boldsymbol{\tau}$ perpendicular to $\partial\Sigma_{S}$ is the unitary outward tangential vector of $\Sigma_{S}$; $\mathcal{P}\subset\Sigma_{S}$ is a finite point set with known physical variable $u_{0}$; $p\left(\bar{\gamma}\right)$ is the material parameter interpolated by the material density $\bar{\gamma}\in\left[0,1\right]$, derived from sequential filtering and projecting the design variable defined on $\Sigma_{S}$ [63, 64, 65]; $\mathbf{g}\left(\nabla_{s}u,u\right)$ is a vector functional of $\nabla_{s}u$ and $u$. Based on the Green’s formula including the tangential gradient operator [66], the physical PDE in equation 1 can be transformed into the following variational formulation:

• •

find $u\in\mathcal{H}\left(\Sigma_{S}\right)$ with $u=u_{0}$ at $\mathcal{P}$, satisfying

[TABLE]

where $\hat{u}$ is the test function of $u$; \mathcal{H}\left(\Sigma_{S}\right)=\big{\{}u\in\mathcal{L}^{2}\left(\Sigma_{S}\right)\big{|}\nabla_{s}u\in\left(\mathcal{L}^{2}\left(\Sigma_{S}\right)\right)^{3}\big{\}} is the first-order Sobolev space on $\Sigma_{S}$, with $\mathcal{L}^{2}\left(\Sigma_{S}\right)$ denoting the second-order Lebesgue integrable functional space on $\Sigma_{S}$; $\mathrm{d}s$ and $\mathrm{d}l$ are the corresponding Riemann metrics on $\Sigma_{S}$ and $\partial\Sigma_{S}$, respectively.

For the case with the physical field defined on a volume domain and design variable defined on the 2-manifold corresponding to a boundary surface of this volume domain, the physical PDE can be expressed in a typical abstract form as

[TABLE]

where $u_{d}$ is the known physical variable on the structure defined on $\Sigma_{S}$; $\Omega$ is the open and bounded volume domain with the boundary of Lipschitz type; $\alpha\left(\bar{\gamma}\right)$ is the material interpolation used to implement the penalization between the Neumann and Dirichlet types of boundary conditions. Based on this material interpolation, the mixed boundary condition is formulated as that in equation 3. When the material density $\bar{\gamma}$ takes the value of [math], $\alpha$ is valued to be large enough to ensure the dominance of the Dirichlet term $\left(u-u_{d}\right)$. The mixed boundary condition degenerates into the Neumann type, when $\alpha$ is valued as [math] with $\bar{\gamma}$ taking the value of $1$. Sequentially, the structural pattern with the known physical variable $u_{d}$ can be determined. The physical PDE in equation 3 can be transformed into the variational formulation as

• •

find $u\in\mathcal{H}\left(\Omega\right)$ with $u=u_{0}$ at $\mathcal{P}$, satisfying

[TABLE]

where \mathcal{H}\left(\Omega\right)=\big{\{}u\in\mathcal{L}^{2}\left(\Omega\right)\big{|}\nabla u\in\left(\mathcal{L}^{2}\left(\Omega\right)\right)^{3}\big{\}} is the first-order Sobolev space defined on $\Omega$. In this case, when the design variable is defined on an interface of the volume domain, the typical abstract form of the physical PDE can be

[TABLE]

where $\llbracket\cdot\rrbracket$ denotes the local jump of a variable across $\Sigma_{S}$ embedded in $\Omega$. Correspondingly, the variational formulation of this physical PDE is

• •

find $u\in\mathcal{H}\left(\Omega\right)$ with $u=u_{d}$ at $\mathcal{P}$ and $u=u_{0}$ on $\partial\Omega$, satisfying

[TABLE]

The above abstract forms of physical PDEs are expressed with divergence operator. They can be directly transformed into the forms with curl operator in the following numerical examples for electromagnetics.

2.3 Regularization

The material density is derived by sequentially implementing the filter and projection operations on the design variable. The filter operation of the design variable is implemented by solving the PDE defined on the 2-manifold [63]:

[TABLE]

where $r$ is the filter radius; $\gamma\in\mathcal{L}^{2}\left(\Sigma_{S}\right)$ is the design variable; $\tilde{\gamma}$ is the filtered design variable; $\hat{\tilde{\gamma}}$ is the test function for $\tilde{\gamma}$. The variational formulation of the PDE filter in equation 7 is

• •

find $\tilde{\gamma}\in\mathcal{H}\left(\Omega\right)$ satisfying

[TABLE]

The filtered design variable is projected by the threshold projection[64, 65], to derive the material density:

[TABLE]

where $\beta$ and $\xi$ are the projection parameters and their values are chosen based on numerical experiments [65].

2.4 Topology optimization problem

The 2-manifold and the material density together composes a fiber bundle [67], where $\Sigma_{S}$ is its base manifold and $\bar{\gamma}:\Sigma_{S}\rightarrow\left[0,1\right]$ is its fibers. This fiber bundle can be expressed as $\left(\Sigma_{S}\times\bar{\gamma}\left(\Sigma_{S}\right),\Sigma_{S},proj_{1},\bar{\gamma}\left(\Sigma_{S}\right)\right)$, with the natural projection $proj_{1}:\Sigma_{S}\times\bar{\gamma}\left(\Sigma_{S}\right)\rightarrow\Sigma_{S}$ satisfying $proj_{1}\left(\mathbf{x},\bar{\gamma}\left(\mathbf{x}\right)\right)=\mathbf{x}$ for $\forall\mathbf{x}\in\Sigma_{S}$. Then the topology optimization problem in this paper can be formulated as an optimization problem constrained by the PDEs defined on the 2-manifold or a volume domain with this 2-manifold as the boundary or interface. The topology optimization problem can be described in the following form:

[TABLE]

where $J:\mathcal{H}\left(\Sigma_{S}\right)\times\mathcal{H}\left(\Sigma_{S}\right)\rightarrow\mathbb{R}$ or $\mathcal{H}\left(\Omega\right)\times\mathcal{H}\left(\Sigma_{S}\right)\rightarrow\mathbb{R}$ is a bounded continuous mapping operator; $V_{f}\in\left(0,1\right)$ is the area fraction of the structure on the 2-manifold with an admissible tolerance of $10^{-3}$; $\left|\Sigma_{S}\right|$ is the area of the 2-manifold $\Sigma_{S}$.

2.5 Adjoint analysis

The adjoint analysis can be implemented to derive the adjoint derivative of the topology optimization problem and define the descent direction of the design objective for the iterative solution procedure. For the PDE used to describe a physical field, the usual situation is that $e$ is continuously Fréchet-differentiable and $e_{u}\left(u;\bar{\gamma}\right)$ is a linear operator with bounded inverse. According to the implicit function theorem [68], $e\left(u;\bar{\gamma}\right)=0$ locally defines a continuously Fréchet-differentiable $\bar{\gamma}\mapsto u\left(\bar{\gamma}\right)$ with the following Fréchet derivative:

[TABLE]

The PDE filter defines a continuously Fréchet-differentiable $\gamma\mapsto\tilde{\gamma}\left(\gamma\right)$ with the following Fréchet derivative:

[TABLE]

Then, the Gâteaux derivative of $J$ is

[TABLE]

where $*$ denotes the adjoint of an operator; $\left\langle\cdot,\cdot\right\rangle$ denotes the pairing between a functional space and its dual. Under the precondition that $J$ is Fréchet-differentiable, the adjoint equations and the adjoint derivative can be derived from the Gâteaux derivative, according to the Kurash-Kuhn-Tucker condition [69].

When the physical field is defined on the 2-manifold $\Sigma_{S}$, $e\left(u;\bar{\gamma}\right)$ can include the terms of surface-integrals and curve-integral; then, it is expressed to be the sum of the corresponding terms as

[TABLE]

By setting

[TABLE]

the adjoint equations can be obtained as

• •

find $\mu\in\mathcal{H}\left(\Sigma_{S}\right)$ and $\nu\in\mathcal{H}\left(\Sigma_{S}\right)$, satisfying

[TABLE]

where $\mu$ and $\nu$ are the adjoint variables of $u$ and $\tilde{\gamma}$, respectively. The adjoint derivative of $J$ can be derived as

[TABLE]

where the adjoint variables $\mu$ and $\nu$ are derived by solving equation 16.

When the physical field is defined on the volume domain $\Omega$, $e\left(u;\bar{\gamma}\right)$ includes the terms of volume-integral and surface-integrals; then, it is expressed to be the sum of the corresponding terms as

[TABLE]

By setting $\mu=-\left(e_{u}^{-1}\right)^{*}J_{u}=\left(e_{u}^{\Omega*}+e_{u}^{\Sigma_{S}*}\right)^{-1}J_{u}$ and $\nu=-\left(f_{\tilde{\gamma}}^{-1}\right)^{*}\bar{\gamma}_{\tilde{\gamma}}J_{\bar{\gamma}}=\left(f_{\tilde{\gamma}}^{*}\right)^{-1}\bar{\gamma}_{\tilde{\gamma}}J_{\bar{\gamma}}$, the adjoint equations can be obtained as

• •

find $\mu\in\mathcal{H}\left(\Omega\right)$ and $\nu\in\mathcal{H}\left(\Sigma_{S}\right)$, satisfying

[TABLE]

where $\mu$ and $\nu$ are the adjoint variables of $u$ and $\tilde{\gamma}$, respectively. The adjoint derivative of $J$ can be derived in the same form as that in equation 17, with the adjoint variables $\mu$ and $\nu$ derived by solving equation 19.

Similarly, for $V={1\over\left|\Sigma_{S}\right|}\int_{\Sigma_{S}}\bar{\gamma}\,\mathrm{d}s$ in the area constraint, the adjoint derivative of $V$ can be derived as

[TABLE]

where the adjoint variable $\nu$ is derived by solving the following adjoint equation:

• •

find $\nu\in\mathcal{H}\left(\Sigma_{S}\right)$ satisfying

[TABLE]

Based on the derived adjoint derivatives with adjoint variables solved from the corresponding adjoint equations, the design variable can be iteratively updated.

2.6 Numerical implementation

By using an iterative procedure, the design variable is penalized and evolved into a binary distribution on the 2-manifold. The iterative procedure is implemented as outlined by the pseudocode in Table 1. The surface finite element method is utilized to solve the relevant PDEs and the adjoint equations defined on the 2-manifold [70]. And it can be implemented by choosing a software-package including a surface finite element solver. In this iterative procedure, the projection parameter $\beta$ with initial value $1$ is doubled after every $30$ iterations and $\xi$ is set to be $0.5$; the loop for solving the optimization problem in equation 10 is stopped when the maximal iteration number is reached, or the averaged variation of the design objective over continuous 5 iterations and the residual of the area constraint are simultaneously less than the specified tolerance of $10^{-3}$; and the design variable is updated using the method of moving asymptotes [71], which has the merits on dealing with multiple integral constraints and bound constraint of the design variable.

3 Test problems

In this section, the first case of topology optimization implemented on 2-manifolds is demonstrated by the optimization of the micro-textures for wetting behaviors on the solid surfaces in the form of 2-manifolds; and the second case is demonstrated by topology optimization of patterns of heat sink for heat transfer and patterns of the perfect conductor for electromagnetics.

3.1 Microtextures for wetting behaviors in Cassie-Baxter mode

Wettability is an important aspect of fluidic mechanics [72]. Constructing artificial micro-textures on a solid surface with complicated geometries can be interesting and attractive for stable hydrophobic wettability. Recently, topology optimization has been implemented to design micro-textures on flat solid surfaces [73]. On solid surfaces, the micro-textures can support two modes of hydrophobic performance, i.e. the Wenzel mode with the liquid completely filling the micro-textures and the Cassie-Baxter mode with vapour pockets trapped in the micro-textures. The Cassie-Baxter mode can be transitioned into the Wenzel mode, when the liquid is pressurized [74, 75]. Therefore, it is desired to use reasonable micro-textures to keep the Cassie-Baxter mode from transition and enhance the stability of hydrophobic wettability. On a solid surface with complicated geometry, the pattern of the micro-textures is defined on a 2-manifold used to describe the solid surface. In this case, the mean curvature of the lquid/vapor interface supported by the micro-textures can be regarded to be a physical field defined on the 2-manifold.

The mean curvature of the lquid/vapor interface supported by the micro-textures on the 2-manifold can be described by the dimensionless Young-Laplace equation [76, 77]:

[TABLE]

where $\bar{z}=z/z_{0}$ is the normalized displacement of the liquid/vapor interface relative to the 2-manifold $\Sigma_{S}$, with $z_{0}$ denoting the magnitude of the original displacement $z$; $L$ is the characteristic size of the solid surface; $\bar{\sigma}=\sigma/\left(LP\right)$ is the dimensionless surface tension, with $\sigma$ denoting the surface tension and $P$ denoting the static pressure at the liquid/vapor interface; $\bar{\sigma}_{l}$ is the dimensionless surface tension of the liquid; $\kappa$ is the mean curvature of the 2-manifold $\Sigma_{S}$. To ensure the uniqueness of the solution, the liquid/vapor interface is constrained at a specified point set $\mathcal{P}$ localized on the 2-manifold.

For the wetting behavior in the Cassie-Baxter mode, the performance of the micro-textures can be measured by the volume of the liquid bulges enclosed by the convex liquid/vapor interface and untextured solid surface. Therefore, the topology optimization of the micro-textures is implemented by minimizing the normalized volume of the liquid bulges expressed as

[TABLE]

which is constrained by an area constraint with the specified area fraction of $V_{f}$. In topology optimization of the micro-textures, the material density is used to interpolate the dimensionless surface tension as

[TABLE]

where $\bar{\sigma}_{s}$ is set as $10^{5}\sup_{\mathbf{x}\in\Sigma_{S}}\left|1-2\bar{\sigma}_{l}\kappa\left(\mathbf{x}\right)\right|$ to approximate the liquid/solid interface.

Based on the adjoint analysis introduced in Section 2.5, the adjoint derivative of $J$ is derived as

[TABLE]

The adjoint variable $\tilde{\gamma}_{a}$ in equation 35 can be derived by solving the following adjoint equations:

• •

find $\bar{z}_{a}\in\mathcal{H}\left(\Sigma_{S}\right)$ with $\bar{z}_{a}=0$ at $\mathcal{P}$, satisfying

[TABLE]

• •

find $\tilde{\gamma}_{a}\in\mathcal{H}\left(\Sigma_{S}\right)$, satisfying

[TABLE]

where $\bar{z}_{a}$ is the adjoint variable of $\bar{z}$. For $V={1\over\left|\Sigma_{S}\right|}\int_{\Sigma_{S}}\bar{\gamma}\,\mathrm{d}s$ in the area constraint, the adjoint derivative is

[TABLE]

where the adjoint variable $\tilde{\gamma}_{a}$ is derived by solving the following adjoint equation:

• •

find $\tilde{\gamma}_{a}\in\mathcal{H}\left(\Sigma_{S}\right)$, satisfying

[TABLE]

After adjoint analysis, the topology optimization problem is solved by the numerical implementation procedure introduced in Section 2.6.

Topology optimization of the micro-textures is implemented for the sphere- and torus-shaped surfaces. By setting the parameters as listed in Table 2, the patterns of the micro-textures are derived for different choices of the points in the point set $\mathcal{P}$. For the sphere, the surface is divided into regular spherical-triangles (Figure 2a) and spherical-quadrangles (Figure 3a), respectively. The point set $\mathcal{P}$ is set as the central points of these spherical-polygons. For the torus, the surface is divided by two sets of circles (Figure 4a and 5a), where $\mathcal{P}$ is set as the intersection points of the circles. Then, the patterns of the micro-textures are derived as shown in Figure 2b and 5b, with the corresponding liquid/vapor interfaces shown in Figure 2c and 5c. In Figure 2d and 5d, 2e and 5e, the derived patterns and the corresponding liquid/vapor interfaces are shown in the deformation views. The convergent histories of the normalized optimization objective and the area constraints are shown in Figure 2f and 4f, including the snapshots for the evolution of the material density. From the convergent histories and the snapshots, the convergent performance of the topology optimization procedures can be confirmed. The topology optimization of the micro-textures is further implemented on a Möbius ring derived by gluing three Möbius strips shown in Figure 6a, where the point set $\mathcal{P}$ is set as the six points sketched in the cutaway view shown in Figure 6b. The micro-textures are derived as shown in Figure 6c with the corresponding liquid/vapor interface shown in Figure 6d.

In Table 3, the volumes of the liquid-bulges are cross compared for the derived micro-textures on the sphere and torus for the cases with two different dimensionless surface tension of $10^{2}$ and $10^{3}$, respectively. The point set $\mathcal{P}$ is set to be the same as that in Figure 2 and 4. From a cross comparison of the values in every row of the sub-tables in Table 3, the improved performance of the micro-textures corresponding to the derived patterns can be confirmed.

When the solid objects are textured using the derived patterns, the three-phase contact lines of the liquid/vapour interfaces can be anchored at the geometrically singular corners formed by the top and side walls of the micro-textures. This anchoring effect is caused by the hysteresis of contact angle. The wetting behaviors are thus in the Cassie-Baxter mode.

If the difference between the static pressure imposed on the liquid is large enough to make the contact angle between the liquid/vapour interface and side wall of the micro-textures reach the critical advancing angle, the liquid/vapour interface can burst and transition can occur from the Cassie-Baxter mode to the Wenzel mode. As the contact angle evolves towards the critical advancing angle, the robustness of the Cassie-Baxter mode decreases; simultaneously, the liquid/vapour interface supported on the micro-textures becomes more convex, and the volume of the liquid bulges suspended at the liquid/vapour interface increases. Therefore, reasonable micro-textures on a solid surface can keep the Cassie-Baxter mode from transition by keeping the contact angle aloof of the crucial advancing angle.

3.2 Heat sink for heat transfer problem

The heat transfer problem has been investigated by using topology optimization to find the layout of the heat conductive materials [78, 79, 80] and optimal match of the structural topology and heat source [81]. This section implements topology optimization on 2-manifolds for heat transfer problems to determine the patterns of heat sinks corresponding to the minimized thermal compliance. The temperature distribution in a 3D domain can be described as

[TABLE]

where $k$ is the thermal conductivity; $T$ is temperature; $Q$ is the heat source; $\Gamma_{d}$ is heat sink with known temperature $T_{d}$ and it is attributed to the Dirichlet type; $\Gamma_{n}$ is the insulation boundary, attributed to the Neumann type; $\Gamma_{d}$ and $\Gamma_{n}$ satisfy $\Gamma_{d}\cup\Gamma_{n}=\Sigma_{M}$; $\Omega$ is the domain enclosed by $\Sigma_{M}$; $T_{0}$ is the known temperature at the points in the point set $\mathcal{P}$. By setting the design domain as $\Sigma_{S}=\Sigma_{M}$, topology optimization can be implemented to determine the crosswise distribution of the heat sink $\Gamma_{d}\left(\bar{\gamma}=0\right)$ and the insulation $\Gamma_{n}\left(\bar{\gamma}=1\right)$ on $\Sigma_{S}$. Then, the material interpolation is implemented on $\Sigma_{S}$ to define the mixed boundary condition:

[TABLE]

where $\alpha\left(\bar{\gamma}\right)$ is the penalization function of the material density. The penalization function $\alpha\left(\bar{\gamma}\right)$ is expressed as

[TABLE]

where $\alpha_{\max}$ and $q$ are the parameters used to implement the penalization and to tune the convexity of the penalization function, respectively. The value of $\alpha_{\max}$ should be chosen to be large enough to ensure the domination of the term $\left(T-T_{d}\right)$ in equation 31, when the material density takes on the value of [math]; equation 31 degenerates into the insulation boundary condition, when the material density takes on the value of $1$. Based on numerical experiments, $\alpha_{\max}$ and $q$ are set as $10^{4}$ and $10^{-6}$, respectively. The variational formulation of equation 10 is

• •

find $T\in\mathcal{H}\left(\Omega\right)$ with $T=T_{0}$ at $\mathcal{P}$, satisfying

[TABLE]

where $\hat{T}$ is the test function of $T$; $\mathcal{H}\left(\Omega\right)$ is the first order Sobolev space on $\Omega$.

The distribution of the heat sink on $\Sigma_{S}$ is determined to minimize the thermal compliance:

[TABLE]

Based on the adjoint analysis introduced in Section 2.5, the adjoint derivative of $J$ is derived as

[TABLE]

The adjoint variable $\tilde{\gamma}_{a}$ in equation 35 can be derived by sequentially solving the following adjoint equations:

• •

find $T_{a}\in\mathcal{H}\left(\Omega\right)$ with $T_{a}=0$ at $\mathcal{P}$, satisfying

[TABLE]

• •

find $\tilde{\gamma}_{a}\in\mathcal{H}\left(\Sigma_{S}\right)$ satisfying

[TABLE]

where $T_{a}$ is the adjoint variables of $T$. For the area constraint, the adjoint derivative and the adjoint equation are the same as those in equations 28 and 29.

Based on the numerical implementation introduced in Section 2.6, the computational domain $\Omega$ is set to be the volume domains with the genus 0 and 1, respectively. The corresponding typical 2-manifolds are sphere and torus (Figure 1a and 1b). The coordinate origin is set as the centers of the sphere and torus; the radius of the sphere is $1$; the inner radius and outer radii of the torus are $3/4$ and $7/4$, respectively. For these two 2-manifolds, the thermal conductivity is set as $1$; the temperature of the heat sink $T_{d}$ and known point-temperature $T_{0}$ are set as [math]; the heat source is set as $Q=1/\left(1+\mathbf{x}^{2}\right)$; the area fractions are set as $0.75$ and $0.7$, respectively. For the sphere, the point set $\mathcal{P}$ is set as the one composed of the vertexes of a cube with the center and one of the vertexes localized at the sphere center and $\left(0.5,0.5,0.5\right)$, respectively. For the torus, $\mathcal{P}$ is set as $\left\{\left(-1.65,0,0\right),\left(0,1.25,0.4\right),\left(0.85,0,0\right),\left(0,-1.25,-0.4\right)\right\}$. The maximal iteration number and the updating interval are set as $n^{sub}_{max}=200$ and $n^{upt}=40$, respectively. The patterns of the heat sinks are derived as shown in Figure 7 and 8, including the corresponding convergent histories and snapshots for the evolution of the material density. From the convergent histories, the convergent performance of the topology optimization procedure can be confirmed for this heat transfer problem. The temperature distributions in the cross-sections of the volume domains are shown in Figure 7b and 8b. The results show that the heat insulations are localized at the parts of the manifolds nearest the zero-temperature points and the heat sinks distribute around the heat insulations. Such distributions of the heat sinks and insulations can preserve the thermal energy in the regions around the zero-temperature points and reduce the temperature gradient to minimize the thermal compliance.

The thermal compliance on 2-manifolds can also be minimized by using topology optimization, where the 2-manifolds are imbedded in a block-shaped 3D domain $\Omega$ enclosed by the insulation boundaries and the center of $\Omega$ is localized at the coordinate origin. In this case, the optimization objective is

[TABLE]

where the 2-manifold $\Sigma_{S}$ is an interface of $\Omega$. The material interpolation implemented on $\Sigma_{S}$ is the mixture of the Dirichlet ($T=T_{d}$) and no-jump (-k\left\llbracket\nabla T\right\rrbracket\cdot\mathbf{n}=0) boundary conditions:

[TABLE]

The adjoint equation in equation 36 is changed as

• •

find $T_{a}\in\mathcal{H}\left(\Omega\right)$ with $T_{a}=0$ at $\mathcal{P}$, satisfying

[TABLE]

By setting the area fraction as $0.6$ and keeping the other parameters without change, topology optimization for the patterns of the heat sinks is implemented on the Möbius strip and Klein bottle which can be derived by gluing two Möbius strips (Figure 9a). Because these two 2-manifolds are non-orientable, the untary normal vector $\mathbf{n}$ is defined locally. The patterns of the heat sinks corresponding to different choices of the point set $\mathcal{P}$ (Figure 9b) are derived as shown in Figure 10 and 11, where the temperature distributions on these 2-manifolds are included. The derived patterns of the heat sinks change the temperature distribution in the volume domain $\Omega$ with the tendency to reduce or eliminate the temperature gradient on the remained part of the interface $\Sigma_{S}$. From the temperature distributions, the performance of the derived patterns of the heat sinks on minimizing the thermal compliance can be confirmed.

To check the optimality, the thermal compliance is cross-compared as listed in Table 4 for the derived patterns of the heat sinks in Figure 10a1, 10b1, 11a1 and 11b1. From a cross comparison of the values in every row of the sub-tables in Table 4, the optimized performance of the derived patterns can be confirmed.

3.3 Perfect conductor for electromagnetics

In electromagnetics, perfect conductor condition is widely used to approximate a completely conductive and thin metal layer [82]. Topology optimization on 2-manifolds can be used to determine the patterns of the perfect conductor by maximizing the scattering energy of an incident field. The electric field scattered by the perfect conductor can be described by the wave equation, where the tangential component of the electric field on the 2-manifold is zero. The scattering field of the perfect conductor is described as

[TABLE]

where $\mathbf{E}_{s}$ and $\mathbf{E}_{i}$ are the scattering and incident fields, respectively; the electric field $\mathbf{E}=\mathbf{E}_{s}+\mathbf{E}_{i}$ is the total field; the second equation is the divergence-free condition of the electric displacement; $\Omega$ is a cuboid-shaped domain. The infinite computational space is truncated by the perfectly matched layers (PMLs). In the PMLs, the wave equations with the complex-valued coordinate transformation is described as [83, 84]

[TABLE]

where $\mathbf{x}^{\prime}$ is the complex-valued coordinate transformed from the original Cartesian coordinate in $\Omega_{P}$; $\mathbf{E}_{s}$ in the PMLs satisfy the divergence-free condition described in the original Cartesian coordinate system, because the fields are source-free in the PMLs; $\nabla_{\mathbf{x^{\prime}}}$ is the gradient operator in the PMLs with transformed coordinates; $\Omega_{P}$ is the union of the PML domains. The transformed coordinates and the original Cartesian coordinates satisfy the following transformation:

[TABLE]

where $\mathbf{T}$ and $\mathbf{x}$ are the transformation matrix and the original Cartesian coordinate, respectively. The transformation matrix is [83]

[TABLE]

where $\lambda$ the wavelength of the incident wave; $d$ the thickness of the PMLs; $\Omega_{fb}$, $\Omega_{lr}$ and $\Omega_{td}$ are the PMLs attached on the surfaces of $\Omega$ with normal vector parallel to $x$, $y$ and $z$ axes, respectively; $\Omega_{fbe}$, $\Omega_{lre}$ and $\Omega_{tde}$ are the PMLs attached on the edges of $\Omega$ with tangential vectors perpendicular to $yOz$, $zOx$ and $xOy$ planes, respectively; $\Omega_{c}$ are the PMLs attached on the vertexes of $\Omega$. The no-jump boundary condition for the scattering field is imposed on the interface $\partial\Omega$ between $\Omega_{P}$ and $\Omega$:

[TABLE]

The perfect electric conductor condition $\mathbf{n}\times\mathbf{E}_{s}=\mathbf{0}$ is imposed on the external boundaries $\Gamma_{D}=\partial\left(\Omega\cup\Omega_{P}\right)$ of the PMLs.

The perfect conductor layer is attached on the 2-manifold $\Sigma_{S}$ immersed in $\Omega$. The design variable is used to indicate the no-jump boundary (\mu_{r}^{-1}\left\llbracket\nabla\times\mathbf{E}\right\rrbracket\times\mathbf{n}=\mathbf{0}) and perfect conductor ($\mathbf{n}\times\mathbf{E}=\mathbf{0}$) parts of $\Sigma_{S}$. The corresponding material interpolation is implemented as

[TABLE]

where $\alpha\left(\bar{\gamma}\right)$ is the penalization function in equation 32. In this penalization function, $\alpha_{\max}$ is chosen as a large but finite value to ensure the domination of the term $\mathbf{n}\times\left(\mathbf{E}_{s}+\mathbf{E}_{i}\right)$ corresponding to perfect conductor boudanry, when the material density takes the value of [math]; equation 46 degenerates into the no-jump boundary condition, when the material density takes the value of $1$. Based on numerical tests, $\alpha_{\max}=1\times 10^{4}$ and $q=1\times 10^{0}$ are chosen in this section. The perfect conductor layer is optimized to maximize the energy of the scattering field:

[TABLE]

Based on the adjoint analysis introduced in [85], the adjoint derivative is derived with the same form as that in equation 35. The adjoint variables are derived by solving the adjoint equations in the following variational formulations:

• •

find $\mathbf{E}_{sa}$ with $Re\left(\mathbf{E}_{sa}\right)\in\mathcal{V}_{\mathbf{E}}$, $Im\left(\mathbf{E}_{sa}\right)\in\mathcal{V}_{\mathbf{E}}$ and $\mathbf{n}\times\mathbf{E}_{sa}=\mathbf{0}$ on $\Gamma_{D}$, satisfying

[TABLE]

• •

find $\tilde{\gamma}_{a}\in\mathcal{H}\left(\Sigma_{s}\right)$ satisfying

[TABLE]

where $Re$ and $Im$ are operators used to extract the real and imaginary parts of a complex; $\mathbf{E}_{sa}$ is the adjoint variables of $\mathbf{E}_{s}$; $\mathcal{V}_{\mathbf{E}}$ is the functional space \big{\{}\mathbf{u}\in\mathcal{H}\left(\mathrm{curl};\Omega\cup\Omega_{P}\right)\big{|}\nabla\cdot\mathbf{u}=0\big{\}} with \mathcal{H}\left(\mathrm{curl};\Omega\cup\Omega_{P}\right)=\big{\{}\mathbf{u}\in\left(\mathcal{L}^{2}\left(\Omega\cup\Omega_{P}\right)\right)^{3}\big{|}\nabla\times\mathbf{u}\in\left(\mathcal{L}^{2}\left(\Omega\cup\Omega_{P}\right)\right)^{3}\big{\}} and $\mathcal{L}^{2}\left(\Omega\cup\Omega_{P}\right)$ denoting the second order Lebesque space for the real functionals defined on $\Omega\cup\Omega_{P}$.

Based on the numerical implementation in Section 2.6, topology optimization for the patterns of the perfect conductor is implemented on the torus and Möbius strip (Figure 12a and 12b). These 2-manifolds are immersed into the volume domain $\Omega$, a brick enclosed by the PMLs. The size of $\Omega$ is $2.4\times 2.4\times 0.96$ in the unit of $\mathrm{m}$. This domain is discretized by the cubic elements with a size of $0.08$. The coordinate origin is set as the center of $\Omega$. For these two 2-manifolds, the area fraction in the area constraint is set as $0.5$. The point set $\mathcal{P}$ is set as shown in Figure 12c and 12d. The maximal iteration number and the updating interval are set as $n^{sub}_{max}=315$ and $n^{upt}=30$, respectively. The incident waves is set to propagate in the $+z$ with a wavelength of $0.8$m. For three different polarizations of the incident waves, the optimized patterns of the perfect conductor are derived as shown in Figure 13 and 14. The evolution snapshots for the material density are shown in Figure 15 and 16 for the incident waves with circular polarization. From the evolution snapshots, the performance of the topology optimization procedures can be confirmed for this electromagnetic problem.

The performance of the derived patterns of the perfect conductor on maximizing the energy of the scattering fields can be confirmed from the vector distributions of the scattering field compared with those of the incident field (Figure 17 and 18). The comparison of the vector distributions shows that the perfect layers can enhance the energy of the scattering field by effectively reflecting and disturbing the incident waves.

To check the optimality, the optimized patterns of the perfect conductor in Figure 13 and 14 are used to scatter the incident waves with different polarizations and the energy of the scattering fields is cross-compared, where the other parameters are kept without change. From the cross comparison of the values in every row of the sub-tables in Table 5, the scattering performance of the derived patterns of the perfect conductor can be confirmed.

4 Conclusions

This paper discussed the topology optimization implemented on 2-manifolds. When a physical field is defined on a 2-manifold, the topology optimization is implemented by interpolating a material parameter in the PDE used to describe this physical field; this case has been demonstrated by the topology optimization of the micro-textures for the wetting behaviors in the Cassie-Baxter mode. When the physical field is defined on a 3D domain and its boundary conditions are defined on a 2-manifold, the material density is used to formulate a mixed boundary condition of the physical field and implement the penalization between two different types of boundary conditions; this case has been demonstrated by the topology optimization of the patterns of the heat sinks for heat transfer and the patterns of the perfect conductor for electromagnetics. Typical 2-manifolds, e.g., sphere, torus, Möbius strip and Klein bottle, have been included in the numerical examples. Based on the homeomorphic property of 2-manifolds, it can be concluded that this topology optimization can be implemented on any compact 2-manifolds.

Acknowledgements

Y. Deng acknowledges the support from a Humboldt Research Fellowship for Experienced Researchers (Humboldt-ID: 1197305), the National Natural Science Foundation of China (No. 51875545), the Youth Innovation Promotion Association of the Chinese Academy of Sciences (No. 2018253), and the Open Fund of SKLAO; Z.  Liu acknowledges the support from the National Natural Science Foundation of China (No. 51675506); J.G. Korvink acknowledges support from an EU2020 FET grant (TiSuMR, 737043), the DFG under grant KO 1883/20-1 Metacoils, in the framework of the German Excellence Intitiative under grant EXC 2082 ”3D Matter Made to Order”, and by the VirtMat initiative ”Virtual Materials Design”. The authors are grateful to Prof. K. Svanberg for supplying the MMA codes.