Topology optimization on two-dimensional manifolds
TTopology optimization on two-dimensional manifolds
Yongbo Deng , ∗ , Zhenyu Liu , Jan G. Korvink † Abstract
This paper implements topology optimization on two-dimensional manifolds. In thispaper, the material interpolation is implemented on a material parameter in the partialdifferential equation used to describe a physical field, when this physical field is de-fined on a two-dimensional manifold; the material density is used to formulate a mixedboundary condition of the physical field and implement the penalization between twodifferent types of boundary conditions, when this physical field is defined on a three-dimensional domain with its boundary conditions defined on the two-dimensional man-ifold corresponding a surface or an interface of this three-dimensional domain. Basedon the homeomorphic property of two-dimensional manifolds, typical two-dimensionalmanifolds, e.g., sphere, torus, M¨obius strip and Klein bottle, are included in the numer-ical tests, which are provided for the problems on fluidic mechanics, heat transfer andelectromagnetics.
Keywords:
Topology optimization; two-dimensional manifold (2-manifold); mate-rial distribution method; tangential gradient operator; mixed boundary condition.
Topology optimization is a robust method used to determine the structural configura-tion, which corresponds to the material distribution in a structure [1]. In contrast todesigning devices by tuning a handful of structural parameters in size and shape opti-mization, topology optimization utilizes the full-parameter space to design structuresbased on the user-desired performance, and it is more flexible and robust, because ofits low dependence on initial guess and implicit expression of the structures.Optimization of structural topology was investigated as early as 1904 for trussesby Michell [2]. Material distribution method for topology optimization was pioneeredby Bendsøe and Kikuchi for elasticity [3], and this method was extended to a varietyof areas, e.g. acoustics, electromagnetics, fluid dynamics and thermodynamics [4–29].Several other methods have also been proposed and developed for the implementationof topology optimization, e.g. the level set method [30–33], the evolutionary techniques[34–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 R or reduced two-dimensional (2D) plane on R . The current de-velopment of additive manufacturing, e.g., 3D-printing, has effectively enlarged thestructural-design space. Implementing topology optimization on 2-manifolds can be ∗ [email protected] (Y. Deng) † [email protected] (J. G. Korvink) a r X i v : . [ phy s i c s . c o m p - ph ] N ov elpful to enhance the freedom of structural design, where a 2-manifold is a topo-logical space with its arbitrary interior point having a neighborhood homeomorphicto a sub-region of R and it can be used to describe a structural surface or a mate-rial interface. Related researches have been implemented for the structural design on2-manifolds based on the conformal geometry theory [40, 41], layouts on surfaces forshell structures [51–55], and material interfaces for stiffness and multi-material struc-tures [42–50], fluid-structure interaction [56–58], energy absorption [59], cohesion [60]and actuation [61].In this paper, topology optimization on 2-manifolds is implemented by using thematerial distribution method and this implementation includes the following two cases.When a physical field is defined on a 2-manifold, the material interpolation of topologyoptimization is implemented on a material parameter in the partial differential equation(PDE) with the tangential gradient operator used to describe the physical field. Whenthe physical field is defined on a 3D domain and its boundary conditions are defined on a2-manifold of codimension one in a 3D Euclidian space, the material density representingthe structural configuration is used to implement the penalization between two differenttypes of boundary conditions and formulate a mixed boundary condition of the physicalfield. In the second case, the 2-manifold can be an exterior surface or an interface ofthe 3D domain. When the 2-manifold is an exterior surface, it can mix the Dirichletand Neumann types of boundary conditions; when the 2-manifold is an interface, it canmix the no-jump and Dirichlet types of boundary conditions.For the two cases of topology optimization on 2-manifolds, a monolithic descriptionof the corresponding optimization problems is presented in Section 2, including theadjoint analysis and numerical implementation; test problems on fluidic mechanics,heat transfer and electromagnetics are implemented in Section 3. In the following, allthe mathematical descriptions are performed in the Cartesian system. Details for topology optimization implemented on 2-manifolds are provided as follows.
According to the classification theorem [62], a 2-manifold without boundary is com-pact if every open cover of it has a finite subcover; its family can be exhausted bythe two infinite families { S (sphere), T (torus), T T (double torus), · · · } and { P (projective plane), P P (Klein bottle), · · · } , where To implement topology optimization on a 2-manifold, a design variable, which is arelaxed binary distribution, is defined on this 2-manifold to represent the structuralconfiguration; and it is bounded in the typical interval [0 , a) Sphere (b) Torus (c) M¨obius strip (d) Klein bottle Figure 1: Typical examples of 2-manifolds. on a PDE filter and threshold projection is imposed on the design variable, with theprojected design variable referred to as the material density. This iterative procedureincluding the PDE filter and the threshold projection can control the feature size of thestructure 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, usedto describe the distribution of the physical field, can be expressed in a typical abstractform with the material parameter interpolated by the material density: ∇ s · [ p (¯ γ ) g ( ∇ s u, u )] = c s , in Σ S [ p (¯ γ ) g ( ∇ s u, u )] · τ = c b , on ∂ Σ S u = u , at P ⊂ Σ S (1)where u is the physical variable; Σ S is the 2-manifold; c s and c b are known distributionsdefined on Σ S and ∂ Σ S , respectively; if Σ S is treated as the 2-manifold of codimensionone in an 3D Euclidian space, ∇ s = ∇ − ( n · ∇ ) n is the tangential gradient operatordefined on Σ S , with n denoting the unitary normal vector of Σ S ; τ perpendicular to ∂ Σ S is the unitary outward tangential vector of Σ S ; P ⊂ Σ S is a finite point set withknown physical variable u ; p (¯ γ ) is the material parameter interpolated by the materialdensity ¯ γ ∈ [0 , S [63–65]; g ( ∇ s u, u ) is a vector functional of ∇ s u and u . Based on theGreen’s formula including the tangential gradient operator [66], the physical PDE inequation 1 can be transformed into the following variational formulation: • find u ∈ H (Σ S ) with u = u at P , satisfying e ( u ; ¯ γ ) := (cid:90) Σ S − p (¯ γ ) g ( ∇ s u, u ) · ∇ s ˆ u − c s ˆ u d s + (cid:90) ∂ Σ s c b ˆ u d l = 0 , ∀ ˆ u ∈ H (Σ S ) (2)where ˆ u is the test function of u ; H (Σ S ) = (cid:8) u ∈ L (Σ S ) (cid:12)(cid:12) ∇ s u ∈ (cid:0) L (Σ S ) (cid:1) (cid:9) isthe first-order Sobolev space on Σ S , with L (Σ S ) denoting the second-order Lebesgueintegrable functional space on Σ S ; d s and d l are the corresponding Riemann metrics onΣ S and ∂ Σ S , respectively.For the case with the physical field defined on a volume domain and design variabledefined 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 ∇ · [ p g ( ∇ u, u )] = c s , in Ω[ p g ( ∇ s u, u )] · n = α (¯ γ ) ( u − u d ) , on Σ S = ∂ Ω u = u , at P ⊂
Ω (3)where u d is the known physical variable on the structure defined on Σ S ; Ω is the openand bounded volume domain with the boundary of Lipschitz type; α (¯ γ ) is the materialinterpolation used to implement the penalization between the Neumann and Dirichlettypes of boundary conditions. Based on this material interpolation, the mixed boundarycondition is formulated as that in equation 3. When the material density ¯ γ takes the alue of 0, α is valued to be large enough to ensure the dominance of the Dirichlet term( u − u d ). The mixed boundary condition degenerates into the Neumann type, when α is valued as 0 with ¯ γ taking the value of 1. Sequentially, the structural pattern with theknown physical variable u d can be determined. The physical PDE in equation 3 can betransformed into the variational formulation as • find u ∈ H (Ω) with u = u at P , satisfying e ( u ; ¯ γ ) := (cid:90) Ω − p g ( ∇ u, u ) · ∇ ˆ u dΩ + (cid:90) Σ S α (¯ γ ) ( u − u d ) ˆ u d s = 0 , ∀ ˆ u ∈ H (Ω) (4)where H (Ω) = (cid:8) u ∈ L (Ω) (cid:12)(cid:12) ∇ u ∈ (cid:0) L (Ω) (cid:1) (cid:9) is the first-order Sobolev space definedon Ω. In this case, when the design variable is defined on an interface of the volumedomain, the typical abstract form of the physical PDE can be ∇ · [ p g ( ∇ u, u )] = c s , in Ω (cid:74) p g ( ∇ s u, u ) (cid:75) · n = α (¯ γ ) ( u − u d ) , on Σ S (cid:44) → Ω u = u d , at P ⊂ Σ S u = u , on ∂ Ω (5)where (cid:74) · (cid:75) denotes the local jump of a variable across Σ S embedded in Ω. Correspond-ingly, the variational formulation of this physical PDE is • find u ∈ H (Ω) with u = u d at P and u = u on ∂ Ω, satisfying e ( u ; ¯ γ ) := (cid:90) Ω − p g ( ∇ u, u ) · ∇ ˆ u dΩ + (cid:90) Σ S α (¯ γ ) ( u − u d ) ˆ u d s = 0 , ∀ ˆ u ∈ H (Ω) (6)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 followingnumerical examples for electromagnetics. The material density is derived by sequentially implementing the filter and projectionoperations on the design variable. The filter operation of the design variable is imple-mented by solving the PDE defined on the 2-manifold [63]: ∇ s · (cid:0) − r ∇ s ˜ γ (cid:1) + ˜ γ = γ, on Σ S − r ∇ s ˜ γ · τ = 0 , at ∂ Σ S (7)where r is the filter radius; γ ∈ L (Σ S ) is the design variable; ˜ γ is the filtered designvariable; ˆ˜ γ is the test function for ˜ γ . The variational formulation of the PDE filter inequation 7 is • find ˜ γ ∈ H (Ω) satisfying f (˜ γ ; γ ) := (cid:90) Σ S r ∇ s ˜ γ · ∇ s ˆ˜ γ + ˜ γ ˆ˜ γ − γ ˆ˜ γ d s = 0 , ∀ ˆ˜ γ ∈ H (Σ S ) (8)The filtered design variable is projected by the threshold projection [64, 65], to derivethe material density: ¯ γ = tanh ( βξ ) + tanh ( β (˜ γ − ξ ))tanh ( βξ ) + tanh ( β (1 − ξ )) (9)where β and ξ are the projection parameters and their values are chosen based onnumerical experiments [65]. The 2-manifold and the material density together composes a fiber bundle [67], where Σ S is its base manifold and ¯ γ : Σ S → [0 ,
1] is its fibers. This fiber bundle can be expressed as(Σ S × ¯ γ (Σ S ) , Σ S , proj , ¯ γ (Σ S )), with the natural projection proj : Σ S × ¯ γ (Σ S ) → Σ S atisfying proj ( x , ¯ γ ( x )) = x for ∀ x ∈ Σ S . Then the topology optimization problemin this paper can be formulated as an optimization problem constrained by the PDEsdefined on the 2-manifold or a volume domain with this 2-manifold as the boundary orinterface. The topology optimization problem can be described in the following form:find γ (Σ S ) ∈ [0 ,
1] for (Σ S × γ (Σ S ) , Σ S , proj , γ (Σ S )) , to minimize J ( u ; ¯ γ ) constrained by e ( u ; ¯ γ ) = 0 (Physical PDE) f (˜ γ ; γ ) = 0 (PDE filter)¯ γ = tanh ( βξ ) + tanh ( β (˜ γ − ξ ))tanh ( βξ ) + tanh ( β (1 − ξ )) (Threshold projection) (cid:12)(cid:12)(cid:12)(cid:12) | Σ S | (cid:90) Σ S ¯ γ d s − V f (cid:12)(cid:12)(cid:12)(cid:12) ≤ − (Area constraint) (10)where J : H (Σ S ) × H (Σ S ) → R or H (Ω) × H (Σ S ) → R is a bounded continuousmapping operator; V f ∈ (0 ,
1) is the area fraction of the structure on the 2-manifoldwith an admissible tolerance of 10 − ; | Σ S | is the area of the 2-manifold Σ S . The adjoint analysis can be implemented to derive the adjoint derivative of the topologyoptimization problem and define the descent direction of the design objective for theiterative solution procedure. For the PDE used to describe a physical field, the usualsituation is that e is continuously Fr´echet-differentiable and e u ( u ; ¯ γ ) is a linear operatorwith bounded inverse. According to the implicit function theorem [68], e ( u ; ¯ γ ) = 0locally defines a continuously Fr´echet-differentiable ¯ γ (cid:55)→ u (¯ γ ) with the following Fr´echetderivative: u ¯ γ (¯ γ ) = − e − u ( u ; ¯ γ ) e ¯ γ ( u ; ¯ γ ) . (11)The PDE filter defines a continuously Fr´echet-differentiable γ (cid:55)→ ˜ γ ( γ ) with the followingFr´echet derivative: ˜ γ γ ( γ ) = − f − γ (˜ γ ; γ ) f γ (˜ γ ; γ ) . (12)Then, the Gˆateaux derivative of J is (cid:104) J (cid:48) , t (cid:105) L (Σ S ) , L (Σ S ) = (cid:104) J u , u ¯ γ ¯ γ ˜ γ ˜ γ γ t (cid:105) U ∗ (Ω) , U (Ω) + (cid:104) J ¯ γ , ¯ γ ˜ γ ˜ γ γ t (cid:105) H (Σ S ) , H (Σ S ) = (cid:10) ˜ γ ∗ γ ¯ γ ˜ γ (cid:0) u ∗ ¯ γ J u + J ¯ γ (cid:1) , t (cid:11) L (Σ S ) , L (Σ S ) = ¨ − f ∗ γ Ä f − γ ä ∗ ¯ γ ˜ γ Ä − e ∗ ¯ γ (cid:0) e − u (cid:1) ∗ J u + J ¯ γ ä , t ∂ L (Σ S ) , L (Σ S ) , ∀ t ∈ L (Σ S ) (13)where ∗ denotes the adjoint of an operator; (cid:104)· , ·(cid:105) denotes the pairing between a func-tional space and its dual. Under the precondition that J is Fr´echet-differentiable, theadjoint equations and the adjoint derivative can be derived from the Gˆateaux derivative,according to the Kurash-Kuhn-Tucker condition [69].When the physical field is defined on the 2-manifold Σ S , e ( u ; ¯ γ ) can include theterms of surface-integrals and curve-integral; then, it is expressed to be the sum of thecorresponding terms as e ( u ; ¯ γ ) = e Σ S ( u ; ¯ γ ) + e ∂ Σ S ( u ) . (14)By setting µ := − (cid:0) e − u (cid:1) ∗ J u = (cid:0) e Σ S ∗ u + e ∂ Σ S ∗ u (cid:1) − J u ,ν := − Ä f − γ ä ∗ ¯ γ ˜ γ J ¯ γ = (cid:0) f ∗ ˜ γ (cid:1) − ¯ γ ˜ γ J ¯ γ , (15)the adjoint equations can be obtained as • find µ ∈ H (Σ S ) and ν ∈ H (Σ S ), satisfying (cid:10)(cid:0) e Σ S ∗ u + e ∂ Σ S ∗ u (cid:1) µ, v (cid:11) H (Σ S ) , H (Σ S ) = (cid:104)− J u , v (cid:105) H (Σ S ) , H (Σ S ) , ∀ v ∈ H (Σ S ) (cid:10) f ∗ ˜ γ ν, y (cid:11) H (Σ S ) , H (Σ S ) = (cid:104)− J ¯ γ ¯ γ ˜ γ , y (cid:105) H (Σ S ) , H (Σ S ) , ∀ y ∈ H (Σ S ) (16) here µ and ν are the adjoint variables of u and ˜ γ , respectively. The adjoint derivativeof J can be derived as (cid:104) J (cid:48) , t (cid:105) L (Σ S ) , L (Σ S ) = ¨ f ∗ γ νJ − γ Ä e Σ S ¯ γ ä ∗ µ + f ∗ γ ν, t ∂ L (Σ S ) , L (Σ S ) , ∀ t ∈ L (Σ S ) (17)where the adjoint variables µ and ν are derived by solving equation 16.When the physical field is defined on the volume domain Ω, e ( u ; ¯ γ ) includes theterms of volume-integral and surface-integrals; then, it is expressed to be the sum ofthe corresponding terms as e ( u ; ¯ γ ) = e Ω ( u ) + e Σ S ( u ; ¯ γ ) . (18)By setting µ = − (cid:0) e − u (cid:1) ∗ J u = (cid:0) e Ω ∗ u + e Σ S ∗ u (cid:1) − J u and ν = − Ä f − γ ä ∗ ¯ γ ˜ γ J ¯ γ = Ä f ∗ ˜ γ ä − ¯ γ ˜ γ J ¯ γ ,the adjoint equations can be obtained as • find µ ∈ H (Ω) and ν ∈ H (Σ S ), satisfying (cid:10)(cid:0) e Ω ∗ u + e Σ S ∗ u (cid:1) µ, v (cid:11) H (Ω) , H (Ω) = (cid:104)− J u , v (cid:105) H (Ω) , H (Ω) , ∀ v ∈ H (Ω) (cid:10) f ∗ ˜ γ ν, y (cid:11) H (Σ S ) , H (Σ S ) = (cid:104)− J ¯ γ ¯ γ ˜ γ , y (cid:105) H (Σ S ) , H (Σ S ) , ∀ y ∈ H (Σ S ) (19)where µ and ν are the adjoint variables of u and ˜ γ , respectively. The adjoint derivativeof J can be derived in the same form as that in equation 17, with the adjoint variables µ and ν derived by solving equation 19.Similarly, for V = | Σ S | (cid:82) Σ S ¯ γ d s in the area constraint, the adjoint derivative of V can be derived as (cid:104) V (cid:48) , t (cid:105) L (Σ S ) , L (Σ S ) = (cid:10) f ∗ γ ν, t (cid:11) L (Σ S ) , L (Σ S ) , ∀ t ∈ L (Σ S ) (20)where the adjoint variable ν is derived by solving the following adjoint equation: • find ν ∈ H (Σ S ) satisfying (cid:10) f ∗ ˜ γ ν, y (cid:11) H (Σ S ) , H (Σ S ) = (cid:104)− V ¯ γ ¯ γ ˜ γ , y (cid:105) H (Σ S ) , H (Σ S ) , ∀ y ∈ H (Σ S ) . (21)Based on the derived adjoint derivatives with adjoint variables solved from the cor-responding adjoint equations, the design variable can be iteratively updated. By using an iterative procedure, the design variable is penalized and evolved into a bi-nary distribution on the 2-manifold. The iterative procedure is implemented as outlinedby the pseudocode in Table 1. The surface finite element method is utilized to solve therelevant PDEs and the adjoint equations defined on the 2-manifold [70]. And it can beimplemented by choosing a software-package including a surface finite element solver.In this iterative procedure, the projection parameter β with initial value 1 is doubledafter every 30 iterations and ξ is set to be 0 .
5; the loop for solving the optimizationproblem in equation 10 is stopped when the maximal iteration number is reached, or theaveraged variation of the design objective over continuous 5 iterations and the residualof the area constraint are simultaneously less than the specified tolerance of 10 − ; andthe design variable is updated using the method of moving asymptotes [71], which hasthe merits on dealing with multiple integral constraints and bound constraint of thedesign variable. In this section, the first case of topology optimization implemented on 2-manifoldsis demonstrated by the optimization of the micro-textures for wetting behaviors onthe solid surfaces in the form of 2-manifolds; and the second case is demonstrated bytopology optimization of patterns of heat sink for heat transfer and patterns of theperfect conductor for electromagnetics. r , V f and n submax ;Set i ← γ ← V f , n sub ← ξ ← . β ← loop Derive ¯ γ by filtering and projecting γ , and evaluate V ;Solve u from the physical PDE, evaluate J , and set J n sub ← J ;Solve µ and ν from the adjoint equations for J , and evaluate J (cid:48) ;Solve ν from the adjoint equation for V , and evaluate V (cid:48) ;Update γ based on J (cid:48) and V (cid:48) ; if mod Ä n sub , n upt ä == 0 β ← β ; end ( if ) if Ä n sub == n submax ä or (cid:16) β == 2 , (cid:80) m =0 (cid:12)(cid:12) J n sub − J n sub − m (cid:12)(cid:12) (cid:46) J ≤ − , and | V − V f | ≤ − (cid:17) break; end ( if ) n sub ← n sub + 1 end ( loop )Table 1: Pseudocode for the iterative solution of topology optimization on 2-manifolds. Inthe loop, n sub is the loop-index, n submax is the maximal value of n sub , n upt is the updatinginterval of the projection parameter β , and J n sub is the objective value corresponding to theloop-index n sub . Wettability is an important aspect of fluidic mechanics [72]. Constructing artificialmicro-textures on a solid surface with complicated geometries can be interesting andattractive for stable hydrophobic wettability. Recently, topology optimization has beenimplemented to design micro-textures on flat solid surfaces [73]. On solid surfaces, themicro-textures can support two modes of hydrophobic performance, i.e. the Wenzelmode with the liquid completely filling the micro-textures and the Cassie-Baxter modewith vapour pockets trapped in the micro-textures. The Cassie-Baxter mode can betransitioned 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 fromtransition and enhance the stability of hydrophobic wettability. On a solid surfacewith complicated geometry, the pattern of the micro-textures is defined on a 2-manifoldused to describe the solid surface. In this case, the mean curvature of the lquid/vaporinterface supported by the micro-textures can be regarded to be a physical field definedon the 2-manifold.The mean curvature of the lquid/vapor interface supported by the micro-textures onthe 2-manifold can be described by the dimensionless Young-Laplace equation [76, 77]: ∇ s · Ñ ¯ σ ∇ s ¯ z » ( L/z ) + |∇ s ¯ z | é = 1 − σ l κ, on Σ S ¯ σ ∇ s ¯ z » ( L/z ) + |∇ s ¯ z | · τ = 0 , at ∂ Σ S ¯ z = 0 , at P (22)where ¯ z = z/z is the normalized displacement of the liquid/vapor interface relativeto the 2-manifold Σ S , with z denoting the magnitude of the original displacement z ; L is the characteristic size of the solid surface; ¯ σ = σ/ ( LP ) is the dimensionlesssurface tension, with σ denoting the surface tension and P denoting the static pressureat the liquid/vapor interface; ¯ σ l is the dimensionless surface tension of the liquid; κ isthe mean curvature of the 2-manifold Σ S . To ensure the uniqueness of the solution,the liquid/vapor interface is constrained at a specified point set P localized on the -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 convexliquid/vapor interface and untextured solid surface. Therefore, the topology optimiza-tion of the micro-textures is implemented by minimizing the normalized volume of theliquid bulges expressed as J = 1 | Σ S | (cid:90) Σ S ¯ z d s (23)which is constrained by an area constraint with the specified area fraction of V f . Intopology optimization of the micro-textures, the material density is used to interpolatethe dimensionless surface tension as¯ σ = ¯ σ l + (¯ σ s − ¯ σ l ) 1 − ¯ γ ¯ γ (24)where ¯ σ s is set as 10 sup x ∈ Σ S | − σ l κ ( x ) | to approximate the liquid/solid interface.Based on the adjoint analysis introduced in Section 2.5, the adjoint derivative of J is derived as (cid:104) J (cid:48) , t (cid:105) L (Σ S ) , L (Σ S ) = − (cid:90) Σ S ˜ γ a t d s, ∀ t ∈ L (Σ S ) . (25)The adjoint variable ˜ γ a in equation 35 can be derived by solving the following adjointequations: • find ¯ z a ∈ H (Σ S ) with ¯ z a = 0 at P , satisfying (cid:90) Σ S | Σ S | ¯ z ˆ¯ z a − ¯ σ ∇ s ¯ z a · ∇ s ˆ¯ z a » ( L/z ) + |∇ s ¯ z | + ¯ σ ( ∇ s ¯ z · ∇ s ¯ z a ) (cid:0) ∇ s ¯ z · ∇ s ˆ¯ z a (cid:1)(cid:16) » ( L/z ) + |∇ s ¯ z | (cid:17) d s = 0 , ∀ ˆ¯ z a ∈ H (Σ S ) (26) • find ˜ γ a ∈ H (Σ S ), satisfying (cid:90) Σ S − r ∇ s ˜ γ a · ∇ s ˆ˜ γ a + ˜ γ a ˆ˜ γ a − ∂ ¯ σ∂ ¯ γ ∂ ¯ γ∂ ˜ γ ∇ s ¯ z · ∇ s ¯ z a » ( L/z ) + |∇ s ¯ z | ˆ˜ γ a d s = 0 , ∀ ˆ˜ γ a ∈ H (Σ S ) (27)where ¯ z a is the adjoint variable of ¯ z . For V = | Σ S | (cid:82) Σ S ¯ γ d s in the area constraint, theadjoint derivative is (cid:104) V (cid:48) , t (cid:105) L (Σ S ) , L (Σ S ) = − (cid:90) Σ S ˜ γ a t d s, ∀ t ∈ L (Σ S ) (28)where the adjoint variable ˜ γ a is derived by solving the following adjoint equation: • find ˜ γ a ∈ H (Σ S ), satisfying (cid:90) Σ S | Σ S | ∂ ¯ γ∂ ˜ γ ˆ˜ γ a d s + r ∇ s ˜ γ a · ∇ s ˆ˜ γ a + ˜ γ a ˆ˜ γ a d s = 0 , ∀ ˆ˜ γ a ∈ H (Σ S ) . (29)After adjoint analysis, the topology optimization problem is solved by the numericalimplementation procedure introduced in Section 2.6.Topology optimization of the micro-textures is implemented for the sphere- andtorus-shaped surfaces. By setting the parameters as listed in Table 2, the patterns of themicro-textures are derived for different choices of the points in the point set P . For thesphere, the surface is divided into regular spherical-triangles (Figure 2a) and spherical-quadrangles (Figure 3a), respectively. The point set P is set as the central points of thesespherical-polygons. For the torus, the surface is divided by two sets of circles (Figure4a and 5a), where P is set as the intersection points of the circles. Then, the patternsof the micro-textures are derived as shown in Figure 2b and 5b, with the correspondingliquid/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 ¯ σ l z V f n submax n upt µ m 10 µ m 0 . deformation views. The convergent histories of the normalized optimization objectiveand the area constraints are shown in Figure 2f and 4f, including the snapshots for theevolution of the material density. From the convergent histories and the snapshots, theconvergent performance of the topology optimization procedures can be confirmed. Thetopology optimization of the micro-textures is further implemented on a M¨obius ringderived by gluing three M¨obius strips shown in Figure 6a, where the point set P is setas the six points sketched in the cutaway view shown in Figure 6b. The micro-texturesare derived as shown in Figure 6c with the corresponding liquid/vapor interface shownin Figure 6d. Figure 2: (a-c) Stereo and vertical views of the sphere surface divided by regular spherical-triangles, the derived pattern of the micro-textures, and the normalized displacement ofthe liquid/vapor interface supported on the micro-textures with the derived pattern; (d-e) deformation views of the micro-textures and the liquid/vapor interfaces; (f) convergenthistories of the optimization objective and the area constraint, including the snapshots forthe evolution of the material density.
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 surfacetension of 10 and 10 , respectively. The point set P is set to be the same as that inFigure 2 and 4. From a cross comparison of the values in every row of the sub-tables inTable 3, the improved performance of the micro-textures corresponding to the derivedpatterns can be confirmed.When the solid objects are textured using the derived patterns, the three-phase con-tact lines of the liquid/vapour interfaces can be anchored at the geometrically singularcorners formed by the top and side walls of the micro-textures. This anchoring effect P ; (c-d) deformationviews of the micro-textures and the liquid/vapor interfaces. (a) Sphere ¯ σ l = 10 . × − . × − ¯ σ l = 10 . × − . × − (b) Torus ¯ σ l = 10 . × − . × − ¯ σ l = 10 . × − . × − Table 3: Volume of the liquid bulges supported by the micro-textures with the derivedpatterns on the sphere and the torus for the dimensionless surface tension 10 and 10 ,respectively. The optimized entries are noted in bold.11 s caused by the hysteresis of contact angle. The wetting behaviors are thus in theCassie-Baxter mode.If the difference between the static pressure imposed on the liquid is large enoughto make the contact angle between the liquid/vapour interface and side wall of themicro-textures reach the critical advancing angle, the liquid/vapour interface can burstand transition can occur from the Cassie-Baxter mode to the Wenzel mode. As thecontact angle evolves towards the critical advancing angle, the robustness of the Cassie-Baxter mode decreases; simultaneously, the liquid/vapour interface supported on themicro-textures becomes more convex, and the volume of the liquid bulges suspended atthe liquid/vapour interface increases. Therefore, reasonable micro-textures on a solidsurface can keep the Cassie-Baxter mode from transition by keeping the contact anglealoof of the crucial advancing angle. The heat transfer problem has been investigated by using topology optimization to findthe layout of the heat conductive materials [78–80] and optimal match of the struc-tural topology and heat source [81]. This section implements topology optimization on2-manifolds for heat transfer problems to determine the patterns of heat sinks corre-sponding to the minimized thermal compliance. The temperature distribution in a 3Ddomain can be described as − ∇ · ( k ∇ T ) = Q, in Ω T = T d , on Γ d − k ∇ T · n = 0 , on Γ n T = T , at P (30)where k is the thermal conductivity; T is temperature; Q is the heat source; Γ d is heatsink with known temperature T d and it is attributed to the Dirichlet type; Γ n is theinsulation boundary, attributed to the Neumann type; Γ d and Γ n satisfy Γ d ∪ Γ n = Σ M ;Ω is the domain enclosed by Σ M ; T is the known temperature at the points in thepoint set P . By setting the design domain as Σ S = Σ M , topology optimization can beimplemented to determine the crosswise distribution of the heat sink Γ d (¯ γ = 0) and theinsulation Γ n (¯ γ = 1) on Σ S . Then, the material interpolation is implemented on Σ S todefine the mixed boundary condition: − k ∇ T · n = α (¯ γ ) ( T − T d ) , on Σ S (31)where α (¯ γ ) is the penalization function of the material density. The penalization func-tion α (¯ γ ) is expressed as α (¯ γ ) = α max q (1 − ¯ γ ) q + ¯ γ (32)where α max and q are the parameters used to implement the penalization and to tunethe convexity of the penalization function, respectively. The value of α max should bechosen to be large enough to ensure the domination of the term ( T − T d ) in equation31, when the material density takes on the value of 0; equation 31 degenerates intothe insulation boundary condition, when the material density takes on the value of 1.Based on numerical experiments, α max and q are set as 10 and 10 − , respectively. Thevariational formulation of equation 10 is • find T ∈ H (Ω) with T = T at P , satisfying (cid:90) Ω k ∇ T · ∇ ˆ T − Q ˆ T d v + (cid:90) Σ S α ( T − T d ) ˆ T d s = 0 , ∀ ˆ T ∈ H (Ω) (33)where ˆ T is the test function of T ; H (Ω) is the first order Sobolev space on Ω.The distribution of the heat sink on Σ S is determined to minimize the thermalcompliance: J = (cid:90) Ω k ∇ T · ∇ T dΩ . (34) ased on the adjoint analysis introduced in Section 2.5, the adjoint derivative of J isderived as (cid:104) J (cid:48) , t (cid:105) L (Σ S ) , L (Σ S ) = − (cid:90) Σ S ˜ γ a t d s, ∀ t ∈ L (Σ S ) . (35)The adjoint variable ˜ γ a in equation 35 can be derived by sequentially solving the fol-lowing adjoint equations: • find T a ∈ H (Ω) with T a = 0 at P , satisfying (cid:90) Ω k ∇ T · ∇ ˆ T a + k ∇ T a · ∇ ˆ T a d v + (cid:90) Σ S αT a ˆ T a d s = 0 , ∀ ˆ T a ∈ H (Ω) (36) • find ˜ γ a ∈ H (Σ S ) satisfying (cid:90) Σ S ∂α∂ ¯ γ ∂ ¯ γ∂ ˜ γ ( T − T d ) T a ˆ˜ γ a d s + r ∇ s ˜ γ a · ∇ s ˆ˜ γ a + ˜ γ a ˆ˜ γ a d s = 0 , ∀ ˆ˜ γ a ∈ H (Σ S ) (37)where T a is the adjoint variables of T . For the area constraint, the adjoint derivativeand the adjoint equation are the same as those in equations 28 and 29.Based on the numerical implementation introduced in Section 2.6, the computa-tional domain Ω 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). Thecoordinate origin is set as the centers of the sphere and torus; the radius of the sphereis 1; the inner radius and outer radii of the torus are 3 / /
4, respectively. Forthese two 2-manifolds, the thermal conductivity is set as 1; the temperature of theheat sink T d and known point-temperature T are set as 0; the heat source is set as Q = 1 / (cid:0) x (cid:1) ; the area fractions are set as 0 .
75 and 0 .
7, respectively. For the sphere,the point set P is set as the one composed of the vertexes of a cube with the center andone of the vertexes localized at the sphere center and (0 . , . , . P is set as { ( − . , , , (0 , . , . , (0 . , , , (0 , − . , − . } . The maxi-mal iteration number and the updating interval are set as n submax = 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 thematerial density. From the convergent histories, the convergent performance of thetopology optimization procedure can be confirmed for this heat transfer problem. Thetemperature distributions in the cross-sections of the volume domains are shown in Fig-ure 7b and 8b. The results show that the heat insulations are localized at the parts ofthe manifolds nearest the zero-temperature points and the heat sinks distribute aroundthe heat insulations. Such distributions of the heat sinks and insulations can preservethe thermal energy in the regions around the zero-temperature points and reduce thetemperature gradient to minimize the thermal compliance.The thermal compliance on 2-manifolds can also be minimized by using topologyoptimization, where the 2-manifolds are imbedded in a block-shaped 3D domain Ωenclosed by the insulation boundaries and the center of Ω is localized at the coordinateorigin. In this case, the optimization objective is J = (cid:90) Σ S k ∇ s T · ∇ s T d s, (38)where the 2-manifold Σ S is an interface of Ω. The material interpolation implementedon Σ S is the mixture of the Dirichlet ( T = T d ) and no-jump ( − k (cid:74) ∇ T (cid:75) · n = 0) boundaryconditions: − k (cid:74) ∇ T (cid:75) · n = α (¯ γ ) ( T − T d ) , on Σ S . (39)The adjoint equation in equation 36 is changed as • find T a ∈ H (Ω) with T a = 0 at P , satisfying (cid:90) Ω k ∇ T a · ∇ ˆ T a d v + (cid:90) Σ S k ∇ s T · ∇ s ˆ T a + αT a ˆ T a d s = 0 , ∀ ˆ T a ∈ H (Ω) . (40)By setting the area fraction as 0 . trip and Klein bottle which can be derived by gluing two M¨obius strips (Figure 9a).Because these two 2-manifolds are non-orientable, the untary normal vector n is definedlocally. The patterns of the heat sinks corresponding to different choices of the pointset P (Figure 9b) are derived as shown in Figure 10 and 11, where the temperaturedistributions on these 2-manifolds are included. The derived patterns of the heat sinkschange the temperature distribution in the volume domain Ω with the tendency toreduce or eliminate the temperature gradient on the remained part of the interface Σ S .From the temperature distributions, the performance of the derived patterns of the heatsinks on minimizing the thermal compliance can be confirmed. (a)Möbius strip 1 Möbius strip 2 Klein bottle xy zO xy zO (b) xy O ∙ ∙ ∙ ∙ P P P P xzO xzO xz O xz O xz O xz O xyO ∙ ∙ xy O x z O ∙ ∙ P P P P Gluing
Figure 9: (a) Sketch for the Klein bottle, derived by gluing two M¨obius strips; (b) sketchesfor the points of the point set P .Figure 10: Derived patterns of the heat sinks and the corresponding temperature distribu-tions on the M¨obius strip: (a1-a2) for the case with P = { P } ; (b1-b2) for the case with P = { P } . P and P has been sketched in Figure 9b. To check the optimality, the thermal compliance is cross-compared as listed in Table4 for the derived patterns of the heat sinks in Figure 10a1, 10b1, 11a1 and 11b1. Froma cross comparison of the values in every row of the sub-tables in Table 4, the optimizedperformance of the derived patterns can be confirmed. P = { P } ;(b1-b2) the case with P = { P } . P and P has been sketched in Figure 9b. (a) M¨obius strip Figure 10(a1) Figure 10(b1) P . . P . . (b) Klein bottle Figure 11(a1) Figure 11(b1) P . . P . . Table 4: Values of the thermal compliance corresponding to the patterns of the heat sinkson the M¨obius strip and the Klein bottle for the different choices of the point set P . Theoptimized entries are noted in bold. 17 .3 Perfect conductor for electromagnetics In electromagnetics, perfect conductor condition is widely used to approximate a com-pletely conductive and thin metal layer [82]. Topology optimization on 2-manifolds canbe used to determine the patterns of the perfect conductor by maximizing the scatteringenergy of an incident field. The electric field scattered by the perfect conductor can bedescribed by the wave equation, where the tangential component of the electric field onthe 2-manifold is zero. The scattering field of the perfect conductor is described as ∇ × (cid:2) µ − r ∇ × ( E s + E i ) (cid:3) − k (cid:15) r ( E s + E i ) = , in Ω ∇ · [ (cid:15) r ( E s + E i )] = 0 , in Ω (41)where E s and E i are the scattering and incident fields, respectively; the electric field E = E s + E i is the total field; the second equation is the divergence-free condition of theelectric displacement; Ω is a cuboid-shaped domain. The infinite computational spaceis truncated by the perfectly matched layers (PMLs). In the PMLs, the wave equationswith the complex-valued coordinate transformation is described as [83, 84] ∇ x (cid:48) × (cid:0) µ − r ∇ x (cid:48) × E s (cid:1) − k (cid:15) r E s = , in Ω P ∇ · E s = 0 , in Ω P (42)where x (cid:48) is the complex-valued coordinate transformed from the original Cartesiancoordinate in Ω P ; E s in the PMLs satisfy the divergence-free condition described in theoriginal Cartesian coordinate system, because the fields are source-free in the PMLs; ∇ x (cid:48) is the gradient operator in the PMLs with transformed coordinates; Ω P is theunion of the PML domains. The transformed coordinates and the original Cartesiancoordinates satisfy the following transformation: x (cid:48) = Tx , ∀ x ∈ Ω P (43)where T and x are the transformation matrix and the original Cartesian coordinate,respectively. The transformation matrix is [83] T = diag ((1 − j ) λ/d, , , in Ω fb diag (1 , (1 − j ) λ/d, , in Ω lr diag (1 , , (1 − j ) λ/d ) , in Ω td diag (1 , (1 − j ) λ/d, (1 − j ) λ/d ) , in Ω fbe diag ((1 − j ) λ/d, , (1 − j ) λ/d ) , in Ω lre diag ((1 − j ) λ/d, (1 − j ) λ/d, , in Ω tde diag ((1 − j ) λ/d, (1 − j ) λ/d, (1 − j ) λ/d ) , in Ω c (44)where λ the wavelength of the incident wave; d the thickness of the PMLs; Ω fb , Ω lr and Ω td are the PMLs attached on the surfaces of Ω with normal vector parallel to x , y and z axes, respectively; Ω fbe , Ω lre and Ω tde are the PMLs attached on the edges ofΩ with tangential vectors perpendicular to yOz , zOx and xOy planes, respectively; Ω c are the PMLs attached on the vertexes of Ω. The no-jump boundary condition for thescattering field is imposed on the interface ∂ Ω between Ω P and Ω: µ − r ( ∇ × E s − ∇ x (cid:48) × E s ) × n = , on ∂ Ω . (45)The perfect electric conductor condition n × E s = is imposed on the external bound-aries Γ D = ∂ (Ω ∪ Ω P ) of the PMLs.The perfect conductor layer is attached on the 2-manifold Σ S immersed in Ω. Thedesign variable is used to indicate the no-jump boundary ( µ − r (cid:74) ∇ × E (cid:75) × n = ) andperfect conductor ( n × E = ) parts of Σ S . The corresponding material interpolationis implemented as µ − r (cid:74) ∇ × ( E s + E i ) (cid:75) × n = α (¯ γ ) [ n × ( E s + E i )] , on Σ S (46)where α (¯ γ ) is the penalization function in equation 32. In this penalization function, α max is chosen as a large but finite value to ensure the domination of the term n × E s + E i ) corresponding to perfect conductor boudanry, when the material density takesthe value of 0; equation 46 degenerates into the no-jump boundary condition, when thematerial density takes the value of 1. Based on numerical tests, α max = 1 × and q = 1 × are chosen in this section. The perfect conductor layer is optimized tomaximize the energy of the scattering field: J = (cid:90) Ω E s · E ∗ s dΩ . (47)Based on the adjoint analysis introduced in [85], the adjoint derivative is derived withthe same form as that in equation 35. The adjoint variables are derived by solving theadjoint equations in the following variational formulations: • find E sa with Re ( E sa ) ∈ V E , Im ( E sa ) ∈ V E and n × E sa = on Γ D , satisfying (cid:90) Ω E ∗ s · ˆ E sa + µ − r ( ∇ × E ∗ sa ) · Ä ∇ × ˆ E sa ä − k (cid:15) r E ∗ sa · ˆ E sa dΩ+ (cid:90) Ω P µ − r ( T ∇ × E ∗ sa ) · Ä T ∇ × ˆ E sa ä | T | − − k (cid:15) r E ∗ sa · ˆ E sa | T | dΩ = 0 , ∀ ˆ E sa ∈ V E (48) • find ˜ γ a ∈ H (Σ s ) satisfying (cid:90) Σ s r ∇ ˜ γ a · ∇ ˆ˜ γ a + ˜ γ a ˆ˜ γ a + ∂α∂ ¯ γ ∂ ¯ γ∂ ˜ γ (cid:2) Re ( n × ( E s + E i )) · Re ( n × E ∗ sa ) + Im ( n × ( E s + E i )) · Im ( n × E ∗ sa ) (cid:3) ˆ˜ γ a dΣ = 0 , ∀ ˆ˜ γ a ∈ H (Σ s ) (49)where Re and Im are operators used to extract the real and imaginary parts of a com-plex; E sa is the adjoint variables of E s ; V E is the functional space (cid:8) u ∈ H (curl; Ω ∪ Ω P ) (cid:12)(cid:12) ∇· u = 0 (cid:9) with H (curl; Ω ∪ Ω P ) = (cid:8) u ∈ (cid:0) L (Ω ∪ Ω P ) (cid:1) (cid:12)(cid:12) ∇ × u ∈ (cid:0) L (Ω ∪ Ω P ) (cid:1) (cid:9) and L (Ω ∪ Ω P ) denoting the second order Lebesque space for the real functionals definedon Ω ∪ Ω P .Based on the numerical implementation in Section 2.6, topology optimization forthe patterns of the perfect conductor is implemented on the torus and M¨obius strip(Figure 12a and 12b). These 2-manifolds are immersed into the volume domain Ω, abrick enclosed by the PMLs. The size of Ω is 2 . × . × .
96 in the unit of m. Thisdomain is discretized by the cubic elements with a size of 0 .
08. The coordinate origin isset as the center of Ω. For these two 2-manifolds, the area fraction in the area constraintis set as 0 .
5. The point set P is set as shown in Figure 12c and 12d. The maximaliteration number and the updating interval are set as n submax = 315 and n upt = 30,respectively. The incident waves is set to propagate in the + z with a wavelength of0 . This paper discussed the topology optimization implemented on 2-manifolds. When aphysical field is defined on a 2-manifold, the topology optimization is implemented by a) Torus (b) M¨obius strip (a) Möbius strip 1 Möbius strip 2 Klein bottle xy zO xy zO
Möbius strip 1 Möbius strip 2 Klein bottle xy zO (b) xy O ∙ ∙ ∙ ∙ P P P P xzO xzO xz O xz O xz O xz O xyO ∙ ∙ xy O x z O ∙ ∙ P P P P (c) P = { P , P } ontorus (a) Möbius strip 1 Möbius strip 2 Klein bottle xy zO xy zO
Möbius strip 1 Möbius strip 2 Klein bottle xy zO (b) xy O ∙ ∙ ∙ ∙ P P P P xzO xzO xz O xz O xz O xz O xyO ∙ ∙ xy O x z O ∙ ∙ P P P P (d) P = { P , P } on M¨obius strip Figure 12: (a-b) Sketches of the torus and M¨obius strip immersed in the volume domain Ω;(c-d) sketches of the point set P on the torus and the M¨obius strip. yz O z x O xy O xy z O xy z O Perfect conductor yz O z x O xy O xy z O Perfect conductor y z O z x O x y O xy z O xy z O Perfect conductor y z O z x O x y O xy z O Perfect conductor yz O z x O x y O xy z O xy z O Perfect conductor yz O z x O x y O xy z O Perfect conductor y z O z x O x y O xy O xy z O xy z O Perfect conductor y z O z x O xy O xy z O xy z O y z O z x O x y O x y O xy z O xy z O Perfect conductor y z O z x O x y O x y O xy z O Perfect conductor x y O Perfect conductor (a) x -polarization yz O z x O xy O xy z O xy z O Perfect conductor yz O z x O xy O xy z O Perfect conductor y z O z x O x y O xy z O xy z O Perfect conductor y z O z x O x y O xy z O Perfect conductor yz O z x O x y O xy z O xy z O Perfect conductor yz O z x O x y O xy z O Perfect conductor y z O z x O x y O xy O xy z O xy z O Perfect conductor y z O z x O xy O xy z O xy z O y z O z x O x y O x y O xy z O xy z O Perfect conductor y z O z x O x y O x y O xy z O Perfect conductor x y O Perfect conductor (b) y -polarization yz O z x O xy O xy z O xy z O Perfect conductor yz O z x O xy O xy z O Perfect conductor y z O z x O x y O xy z O xy z O Perfect conductor y z O z x O x y O xy z O Perfect conductor yz O z x O x y O xy z O xy z O Perfect conductor yz O z x O x y O xy z O Perfect conductor y z O z x O x y O xy O xy z O xy z O Perfect conductor y z O z x O xy O xy z O xy z O y z O z x O x y O x y O xy z O xy z O Perfect conductor y z O z x O x y O x y O xy z O Perfect conductor x y O Perfect conductor (c) Circular polarization
Figure 13: Derived patterns of the perfect conductor on the torus: (a) for the incident wavepolarized in the x -axis; (b) for the incident wave polarized in the y -axis; (c) for the incidentwave right-circularly polarized in the xOy -plane. yz O zx O xy O xy z O xy z O Perfect conductor yz O zx O xy O xy z O Perfect conductor yz O zx O x y O xy z O xy z O Perfect conductor yz O zx O x y O xy z O Perfect conductor yz O zx O x y O xy z O xy z O Perfect conductor yz O zx O x y O xy z O Perfect conductor (a) x -polarization yz O zx O xy O xy z O xy z O Perfect conductor yz O zx O xy O xy z O Perfect conductor yz O zx O x y O xy z O xy z O Perfect conductor yz O zx O x y O xy z O Perfect conductor yz O zx O x y O xy z O xy z O Perfect conductor yz O zx O x y O xy z O Perfect conductor (b) y -polarization yz O zx O xy O xy z O xy z O Perfect conductor yz O zx O xy O xy z O Perfect conductor yz O zx O x y O xy z O xy z O Perfect conductor yz O zx O x y O xy z O Perfect conductor yz O zx O x y O xy z O xy z O Perfect conductor yz O zx O x y O xy z O Perfect conductor (c) Circular polarization
Figure 14: Derived patterns of the perfect conductor on the M¨obius strip: (a) for theincident wave polarized in the x -axis; (b) for the incident wave polarized in the y -axis; (c)for the incident wave right-circularly polarized in the xOy -plane.20 a) Iteration 1 (b) Iteration 5 (c) Iteration 20 (d) Iteration 70 (e) Iteration 100(f) Iteration 110 (g) Iteration 150 (h) Iteration 200 (i) Iteration 260 (j) Iteration 310 Figure 15: Snapshots for the evolution of the material density in the topology optimizationof the patterns of the perfect conductor on the torus, where the incident wave is right-circularly polarized in the xOy -plane. (a) Iteration 1 (b) Iteration 5 (c) Iteration 20 (d) Iteration 40 (e) Iteration 70(f) Iteration 110 (g) Iteration 150 (h) Iteration 200 (i) Iteration 260 (j) Iteration 310 Figure 16: Snapshots for the evolution of the material density in the topology optimizationof the patterns of the perfect conductor on the M¨obius strip, where the incident wave isright-circularly polarized in the xOy -plane. (a) x -polarization (b) y -polarization (c) Circular polarization Figure 17: Distributions of the vectors of the scattering field (red arrows) and the incidentfield (blue arrows) correspond to the three cases in Figure 13, respectively.21 a) x -polarization (b) y -polarization (c) Circular polarization Figure 18: Distributions of the vectors of the scattering field (red arrows) and the incidentfield (blue arrows) correspond to the three cases in Figure 14, respectively. (a) Torus
Figure 13(a) Figure 13(b) Figure 13(c) x -polarization . . . y -polarization 58 . . . . . . (b) M¨obius strip Figure 14(a) Figure 14(b) Figure 14(c) x -polarization . . . y -polarization 18 . . . . . . Table 5: Energy of the scattering field corresponds to the patterns of the perfect conductoroptimized on the torus and the M¨obius strip for the incident waves with different polariza-tions. The optimized entries are noted in bold.22 nterpolating a material parameter in the PDE used to describe this physical field; thiscase has been demonstrated by the topology optimization of the micro-textures for thewetting behaviors in the Cassie-Baxter mode. When the physical field is defined on a 3Ddomain and its boundary conditions are defined on a 2-manifold, the material densityis used to formulate a mixed boundary condition of the physical field and implementthe penalization between two different types of boundary conditions; this case has beendemonstrated by the topology optimization of the patterns of the heat sinks for heattransfer and the patterns of the perfect conductor for electromagnetics. Typical 2-manifolds, e.g., sphere, torus, M¨obius strip and Klein bottle, have been included inthe numerical examples. Based on the homeomorphic property of 2-manifolds, it canbe 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 Experi-enced Researchers (Humboldt-ID: 1197305), the National Natural Science Foundationof China (No. 51875545), the Youth Innovation Promotion Association of the Chi-nese Academy of Sciences (No. 2018253), and the Open Fund of SKLAO; Z. Liuacknowledges 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 Ger-man Excellence Intitiative under grant EXC 2082 ”3D Matter Made to Order”, and bythe VirtMat initiative ”Virtual Materials Design”. The authors are grateful to Prof. K.Svanberg for supplying the MMA codes.
References [1] M. P. Bendsøe, O. Sigmund, Topology optimization—theory methods and applica-tions, Springer, Berlin, 2003.[2] A. G. M. Michell, The limit of economy of material in frame-structures,
Phil. Mag. , 589-597.[3] M. Bendsøe, N. Kikuchi, Generating optimal topologies in optimal design using ahomogenization method, Comput. Methods Appl. Mech. Eng. , 197-224.[4] O. Sigmund, A 99-line topology optimization code written in Matlab, Struct. Mul-tidisc. Optim. 21 (2001) 120-127.[5] O. Sigmund, On the design of compliant mechanisms using topology optimization,Mech. Struct. Mach. 25 (1997) 495-526.[6] G. I. N. Rozvany, Aims scope methods history and unified terminology of computer-aided optimization in structural mechanics, Struct. Multidisc. Optim. , 2001, , 90-108.[7] M. P. Bendsøe, O. Sigmund, Material interpolations in topology optimization, Arch.Appl. Mech. 69 (1999) 635-654.[8] A. Gersborg-Hansen, M.P. Bendsøe, O. Sigmund, Topology optimization of heatconduction problems using the finite volume method, Struct. Multidisc. Optim. , 251-259.[9] W. Akl, A. El-Sabbagh, K. Al-Mitani, A. Baz, Topology optimization of a platecoupled with acoustic cavity, Int. J. Solids Struct. , 2060-2074.[10] T. Borrvall, J. Petersson, Topology optimization of fluid in Stokes flow, Int. J.Numer. Methods Fluids , 77-107.[11] O. Sigmund, K. G. Hougaard, Geometric properties of optimal photonic crystals, Phys. Rev. Lett. , 153904.[12] T. Nomura, K. Sato, K. Taguchi, T. Kashiwa, S. Nishiwaki, Structural topologyoptimization for the design of broadband dielectric resonator antennas using the finitedifference time domain technique,
Int. J. Numer. Methods Eng. , 1261-1296.
13] M. B. Duhring, J. S. Jensen, O. Sigmund, Acoustic design by topology optimiza-tion,
J. Sound Vibr. , 557-575.[14] S. Kreissl, G. Pingen, K. Maute, An explicit level-set approach for generalizedshape optimization of fluids with the lattice Boltzmann method,
Int. J. Numer.Meth. Fluids , 2011, ,496-519.[15] S. Zhou, Q. Li, A variational level set method for the topology optimization ofsteady-state Navier-Stokes flow, J. Comput. Phys. , 2008, , 10178-10195.[16] J. K. Guest, J. H. Prevost, Topology optimization of creeping fluid flows using aDarcy-Stokes finite element,
Int. J. Numer. Methods Eng. , 2006, , 461-484.[17] A. Takezawa, M. Haraguchi. T. Okamoto, M. Kitamura, Cross-sectional optimiza-tion of whispering-gallery mode sensor with high electric field intensity in the detec-tion domain, IEEE. J. Sel. Top. Quant. Electron. (6), 1-10.[18] Y. Deng, Z. Liu, P. Zhang, Y. Liu, Y. Wu, Topology optimization of unsteadyincompressible Navier-Stokes flows, J. Comput. Phys., 230 (2011), 6688-6708.[19] J. Andkjær, O. Sigmund, Topology optimized low-contrast all-dielectric opticalcloak, Appl. Phys. Lett. , 021112.[20] J. Andkjær, N. A. Mortensen, O. Sigmund, Towards all-dielectric, polarization-independent optical cloaks, Appl. Phys. Lett. , 101106.[21] G. Fujii, H. Watanabe, T. Yamada, T. Ueta, M. Mizuno, Level set based topologyoptimization for optical cloaks,
Appl. Phys. Lett. , 251106.[22] A. R. Diaz, O. Sigmund, A topology optimization method for design of negativepermeability metamaterials,
Struct Multidisc Optim , 163-177.[23] S. Zhou, W. Li, G. Sun, Q. Li, A level-set procedure for the design of electromag-netic metamaterials, Optics Express , 6693-6702.[24] S. Zhou, W. Li, Y. Chen, G. Sun, Q. Li, Topology optimization for negative perme-ability metamaterials using level-set algorithm, Acta Materialia , 2624-2636.[25] J. Andkjær, S. Nishiwaki, T. Nomura, O. Sigmund, Topology optimization of grat-ing couplers for the efficient excitation of surface plasmons, JOSA B , 1828-1832.[26] S. Zhou, W. Li, Q. Li, Level-set based topology optimization for electromagneticdipole antenna design, Journal of Computational Physics , 6915-6930.[27] E. Hassan, E. Wadbro, M. Berggren, Topology optimization of metallic antennas,
IEEE Transactions on Antennas And Propagations , 2488-2500.[28] H. Shim, V. T. T. Ho, S. Wang, D. A. Tortorelli, Level set-based topology op-timization for electromagnetic systems, IEEE Transactions on Magnetics ,1582-1585.[29] T. Feichtner, O. Selig, M. Kiunke, B. Hecht, Evolutionary optimization of opticalantennas, Phys. Rev. Lett. , 127701.[30] M. Y. Wang, X. Wang, D. Guo, A level set method for structural optimization,Comput. Methods Appl. Mech. Eng. 192 (2003) 227-246.[31] G. Allaire, F. Jouve, A. Toader, Structural optimization using sensitivity analysisand a level-set method, J. Comput. Phys. 194 (2004) 363-393.[32] Z. Liu, J.G. Korvink, Adaptive moving mesh level set method for structure opti-mization, Eng. Optim. 40 (2008) 529-558.[33] X. Xing, P. Wei, M.Y. Wang, A finite element-based level set method for structuraloptimization, Int. J. Numer. Methods Eng. 82 (2010) 805-842.[34] Y. M. Xie, G. P. Steven, Evolutionary structural optimization, Springer, 1997.[35] G. P. Steven, Q. Li, Y. M. Xie, Evolutionary topology and shape design for physicalfield problems, Comput. Mech., 26 (2000), 129-139.[36] P. Tanskanen, The evolutionary structural optimization method: theoretical as-pects, Comput. Methods Appl. Mech. Engrg., 191 (2002), 47-48.[37] Guo X, Zhang W, Zhong W (2014) Doing topology optimization explicitly andgeometrically — a new moving morphable components based framework. J ApplMech 81:081009.
38] Guo X, Zhang W, Zhang J, Yuan J (2016) Explicit structural topology optimiza-tionbased onmoving morphable components (MMC) with curved skeletons. ComputMethods Appl Mech Eng, 310:711-748.[39] Akihiro Takezawa, ShinjiNishiwaki, MitsuruKitamura (2010) Shape and topologyoptimization based on the phase field method and sensitivity analysis. Journal ofComputational Physics, 229:2697-2718.[40] Qian Ye, Yang Guo, Shikui Chen, Na Lei, Xianfeng Gu, Topology optimizationof conformal structures on manifolds using extended level set methods (X-LSM) andconformal geometry theory, Computer Methods in Applied Mechanics and Engineer-ing, 344 (2019): 164-185.[41] Panagiotis Vogiatzis, Ming Ma, Shikui Chen and Xianfeng Gu, Computationaldesign and additive manufacturing of periodic conformal metasurfaces by synthesiz-ing topology optimization with conformal mapping, Computer Methods in AppliedMechanics and Engineering, 328 (2018): 477-497.[42] G. Allaire, C. Dapogny, G. Delgado, G. Michailidis. Multi-phase structural op-timization via a level set method,
ESAIM: Control, Optimisation and Calculus ofVariations , 2014, 20, 576-611.[43] N. Vermaak, G. Michailidis, G. Parry, R. Estevez, G. Allaire, Y. Br´echet. Materialinterface effects on the topology optimization of multi-phase structures using a levelset method,
Struct. Multidisc. Optim.
J. Mech. Phys. Solids
J. Mech. Phys. Solids
Int. J. Numer. Meth. Engng.
AIAA J.
Struct. Mul-tidisc. Optim.
Comput. Methods Appl. Mech. Eng.
Struct. Multidisc. Optim.
Computersand Structures
Computers and Structures
Struct. Multidisc. Optim.
Shell structures for architecture-formfinding and optimization , Routledge, New York, 2014.[55] A. Clausen, E. Andreassen, O. Sigmund, Topology optimization of 3D shell struc-tures with porous infill,
Acta Mechanica Sinica
Int. J. Numer. Meth. Engng
57] C. Lundgaard, J. Alexandersen, M. Zhou, C. S. Andreasen, O. Sigmund, Revisitingdensity-based topology optimization for fluid-structure-interaction problems,
Struct.Multidisc. Optim.
Struct. Multidisc. Optim.
A topology optimization interface for LS-DYNA . In: 11.LS-DYNA Forum, Ulm, 2012.[60] R. Behrou, M. Lawry, K. Maute, Level set topology optimization of structuralproblems with interface cohesion,
Int. J. Numer. Meth. Engng.
Struct. Multidisc. Optim.
Int. J. Numer. Methods Eng.
Struct. Multidiscip. Optim.
Int. J. Numer.Methods Eng.
Elliptic partial differential equations of second order ,Springer, 1988.[67] S. S. Chern, W. H. Chen, K. S. Lam,
Lectures on differential geometry , WorldScientific Publishing, 1999.[68] E. Zeidler, Nonlinear Functional Analysis and Its Applications. I, Fixed-Point The-orems, Springer, Berlin, 1986.[69] M. Hinze, R. Pinnau, M. Ulbrich, S. Ulbrich, Optimization with PDE Constraints,Springer, Berlin, 2009.[70] G. Dziuk, C. M. Elliott, Finite element methods for surface PDEs, Acta Numerica.(2013) 289-396.[71] K. Svanberg, The method of moving asymptotes: a new method for structuraloptimization, Int. J. Numer. Methods Eng. 24 (1987) 359-373.[72] X. Feng, L. Jiang, Design and creation of superwetting/antiwetting surfaces.
Adv.Mater.
Comput. Methods Appl. Mech. Engrg.
Europhys. Lett.
Nat. Mater.
Phil. Trans.
81] Y. Deng, Z. Liu, Y. Liu, Y. Wu, Combination of topology optimization and optimalcontrol method, Journal of Computational Physics 257 (2014) 374-399.[82] J. D. Kraus, K. R. Carver, Electromagnetics, 1973, McGraw-Hill edition, in En-glish, 2nd ed.[83] J. Jin, The Finite Element Method in Electromagnetics, 2nd edition, John Wiley& Sons, New York, 2002.[84] J. Berenger, A perfectly matched layer for the absorption of electromagnetic waves,J. Comput. Phys. 114 (1994) 185-200.[85] Y. Deng, J. G. Korvink, Self-consistent adjoint analysis for topology optimizationof electromagnetic waves, J. Comput. Phys. 361 (2018) 353-376.81] Y. Deng, Z. Liu, Y. Liu, Y. Wu, Combination of topology optimization and optimalcontrol method, Journal of Computational Physics 257 (2014) 374-399.[82] J. D. Kraus, K. R. Carver, Electromagnetics, 1973, McGraw-Hill edition, in En-glish, 2nd ed.[83] J. Jin, The Finite Element Method in Electromagnetics, 2nd edition, John Wiley& Sons, New York, 2002.[84] J. Berenger, A perfectly matched layer for the absorption of electromagnetic waves,J. Comput. Phys. 114 (1994) 185-200.[85] Y. Deng, J. G. Korvink, Self-consistent adjoint analysis for topology optimizationof electromagnetic waves, J. Comput. Phys. 361 (2018) 353-376.