Fiber bundle topology optimization for surface flows

This paper presents a topology optimization approach for the surface flows on variable design domains. Via this approach, the matching between the pattern of a surface flow and the 2-manifold used to define the pattern can be optimized, where the 2-manifold is implicitly defined on another fixed 2-manifold named as the base manifold. The fiber bundle topology optimization approach is developed based on the description of the topological structure of the surface flow by using the differential geometry concept of the fiber bundle. The material distribution method is used to achieve the evolution of the pattern of the surface flow. The evolution of the implicit 2-manifold is realized via a homeomorphous map. The design variable of the pattern of the surface flow and that of the implicit 2-manifold are regularized by two sequentially implemented surface-PDE filters. The two surface-PDE filters are coupled, because they are defined on the implicit 2-manifold and base manifold, respectively. The surface Navier-Stokes equations, defined on the implicit 2-manifold, are used to describe the surface flow. The fiber bundle topology optimization problem is analyzed using the continuous adjoint method implemented on the first-order Sobolev space. Several numerical examples have been provided to demonstrate this approach, where the combination of the viscous dissipation and pressure drop is used as the design objective.


Introduction
Surface flows can greatly decrease the computational cost in the numerical design of the related fluidic structures.The fluid flows in the channels attached on the walls of equipments can be described as surface flows on the curved surfaces corresponding to the outer shapes of fluidic structures.The streamsurfaces corresponding to the outer shapes of fluidic structures with complete-slip boundaries can be described as surface flows separated from the bulk flows, where the complete-slip boundaries can be approximated and achieved by chemically coating or physically structuring solid surfaces to derive the extreme hydrophobicity [1], using the optimal control method to manipulate the boundary velocity of flows [2], and producing vapor layers between the solid and liquid phases based on the Leidenfrost phenomenon [3], etc.
The topological structure of a surface flow can be described as the fiber bundle demonstrated in Fig. 1.Fiber bundle is a concept of differential geometry [4].It is composed of the base manifold and the fiber defined on it, where the manifold represents the topological space locally homeomorphous to an Euclidean space.For the surface flow, the flow pattern together with its definition domain corresponds to the fiber of the fiber bundle.If there exists a 2-manifold homeomorphous to the fiber, it can be set as the base manifold of the fiber bundle.In computation, the existence of the base manifold can be ensured by presetting a fixed geometrical surface as the base manifold, then the fiber can be found on the preset base manifold.That means the definition domain of the pattern is an implicit 2-manifold defined on the preset base manifold, where the implicit 2-manifold is the fluid/solid interface corresponding to the outer-shape surface.The reason for this paper to use the concept of fiber bundle is to describe the topological structure of the surface flow as an ensemble instead of three separated components.Therefore, the task for the fiber bundle topology optimization of the surface flow is to find the optimized matching between the pattern and the implicit 2-manifold defined on the preset base manifold.: Sketch for the fiber bundle of a surface flow, where Σ is the base manifold, Γ is the implicit 2-manifold used to define the pattern of the surface flow, γ p : Γ → [0, 1] is the pattern of the surface flow, u is the fluid velocity of the surface flow, u lv,Γ is the known fluid velocity at the boundary of Γ, n Γ is the unitary normal vector of Γ, τ Γ is the unitary tangential vector at ∂Γ, n τ Γ = n Γ × τ Γ is the outward unitary normal at ∂Γ, x Σ denotes a point on Σ, and x Γ denotes a point on Γ.This paper focuses on the laminar surface flows with low and moderate Reynolds numbers to demonstrate the fiber bundle topology optimization approach, although the sketched surface flows can be turbulent with high Reynolds number.
It is natural for one to ask if it is possible to implement the topology optimization to match the pattern of a surface flow and the implicit 2-manifold on which the pattern is defined.If such topology optimization can be achieved, the design space and design freedom will be further extended for flow problems by including the design domain for the pattern of the surface flow into the design space, where the design domain is the implicit 2-manifold.Therefore, this paper presents the fiber bundle topology optimization approach for the surface flow.
For the fiber bundle topology optimization approach, the material distribution method pioneered by [6] is used to determine the pattern of the surface flow.The implicit 2-manifold used to define the surface flow is described on the base manifold.Then, two sets of design variables are required for the pattern of the surface flow and the implicit 2-manifold, respectively.For the material distribution method, a porous medium model has been developed for Stokes flows [10].This model was then extended to implement topology optimization for steady and unsteady Navier-Stokes flows [12][13][14].In this model, the porous medium was filled in the two/threedimensional design domains.Correspondingly, an artificial Darcy friction was introduced into the force terms of the Stokes equations and Navier-Stokes equations.The impermeability of the porous medium was evolved in the topology optimization procedure to derive the fluidic structures.Inspired by the porous medium model, topology optimization for surface flows has been implemented by filling the porous medium onto fixed 2-manifolds, where an artificial Darcy friction is added to the surface Navier-Stokes equations [8].This paper inherits this model with the porous medium to implement the fiber bundle topology optimization for surface flows.
The remained sections of this paper are organized as follows.In Section 2, a monolithic description of the fiber bundle topology optimization problem for a surface flow is presented.In Section 3, numerical implementation for the iterative solution of the fiber bundle topology optimization problem is introduced.In Section 4, numerical tests are provided to demonstrate the developed fiber bundle topology optimization approach.In Sections 5 and 6, the conclusion and acknowledgment of this paper are provided.In Section 7, details are provided for the adjoint analysis of the fiber bundle topology optimization problem.All the mathematical descriptions are implemented in a Cartesian system.

Methodology
In this section, the fiber bundle topology optimization problem is described to match the pattern of the surface flow and the implicit 2-manifold on which the surface flow is defined.The implicit 2-manifold is defined on the base manifold.The incompressible surface fluid is considered.

Physical model and material interpolation
In the fiber bundle topology optimization for the surface flow, the porous medium model is utilized.In this model, the porous medium is filled onto the implicit 2-manifold.Correspondingly, the artificial Darcy friction is added to the surface Navier-Stokes equations.The artificial Darcy friction is derived based on the constitutive law of the porous medium.It is assumed to be proportional to the fluid velocity [10,12]: where α is the impermeability; Γ is the implicit 2-manifold; and xΓ denotes a point on Γ.
When the porosity of the porous medium is zero, it corresponds to a solid material with infinite impermeability and zero fluid velocity caused by the infinite friction force.When the porosity is infinite, it corresponds to the structural void for the transport of the fluid with zero impermeability.Therefore, the impermeability can be described as where γp ∈ {0, 1} is a binary distribution defined on Γ, with 0 and 1 representing the solid and fluid phases, respectively; ΓD is the design domain for the pattern of the surface flow, ΓF is the fluid domain with the material density enforced to be γp = 1, respectively, where ΓD and ΓF satisfy ΓD ∪ ΓF = Γ and ΓD ∩ ΓF = ∅.Especially, Γ is the design domain, when there is no enforced fluid domain, i.e., ΓF = ∅ and Γ = ΓD.
To avoid the numerical difficulty on solving a binary optimization problem, the binary variable γp in the design domain is relaxed to vary continuously in [0, 1].The relaxed binary variable is referred to as the material density of the impermeability.Based on the description of the impermeability in Eq. 2, the material interpolation of the impermeability can be implemented by using the convex and q-parameterized scheme [10]: where αs and α f are the impermeability of the solid and fluid phases, respectively; q is the parameter used to tune the convexity of this interpolation.For the fluid phase, the impermeability is zero, i.e., α f = 0.For the solid phase, αs should be infinite theoretically; numerically, a finite value much larger than the fluid density ρ is chosen for αs, to ensure the stability of the numerical implementation and approximate the solid phase with enough accuracy.Based on numerical tests, q is valued as 1 and αs is chosen as 10 4 ρ to satisfy αs ρ in this paper.

Design variables
In the fiber bundle topology optimization for the surface flow, two sets of design variables are required to be defined for the implicit 2-manifold and the pattern of the surface flow, respectively.

Design variable for implicit 2-manifold
To describe the implicit 2-manifold, the design variable that takes continuous values in [0, 1] is defined on the base manifold.This design variable is used to describe the distribution of the normal displacement of the implicit 2-manifold relative to the base manifold.Equivalently, the pattern of the surface flow is defined on a variable design domain.Then, the result of the fiber bundle topology optimization can be regarded to be a two-order hierarchical structure composed of the base and secondary structures corresponding to the implicit 2-manifold and the pattern of the surface flow, respectively.
To control the smoothness of the implicit 2-manifold and ensure the well-poseness of the solution, a surface-PDE filter sketched in Fig. 2 is imposed on the design variable of the implicit 2-manifold [50]: where dm is the design variable for the implicit 2-manifold; d f is the filtered design variable; rm is the filter radius, and it is constant; Σ is the base manifold used to define the implicit 2-manifold; xΣ denotes a point on Σ; ∇Σ and divΣ are the tangential gradient operator and tangential divergence operator defined on Σ, respectively; nτ Σ = nΣ ×τ Σ is the outward unitary conormal vector normal to ∂Σ and tangent to Σ at ∂Σ, with nΣ and τ Σ representing the unitary normal vector on Σ and the unitary tangential vector at ∂Σ, respectively; A d is a parameter used to specify the amplitude of the normal displacement of the implicit 2-manifold relative to the base manifold, and it is nonnegative After the filter operation, the implicit 2-manifold can be described by the filtered design variable as where Γ is the implicit 2-manifold; xΓ denotes a point on Γ.From Eq. 5, a differential homeomorphism can be determined corresponding to the bijection The Jacobian matrix of the homeomorphism (Eq.5) for the implicit 2-manifold in the curvilinear coordinate system of the base manifold can be transformed into the following formulation: with ∂x Γ ∂x Σ representing its determinant.The variational formulation of the surface-PDE filter in Eq. 4 is considered in the first order Sobolev space defined on Σ.It can be derived based on the Galerkin method as where df is the test function of d f ; H (Σ) represents the first order Sobolev space defined on Σ; L 2 (Σ) represents the second order Lebesque space defined on Σ.
Surface-PDE filter in Eq. 4 Figure 2: Sketch for the surface-PDE filter for the design variable of the implicit 2-manifold Γ defined on the base manifold Σ.

Design variable for pattern of surface flow
The pattern of the surface flow is represented by the material density defined on the implicit 2-manifold.The material density in Eqs. 2 and 3 is obtained by sequentially implementing the surface-PDE filter and the threshold projection on the design variable for the material density, as sketched in Fig. 3.This design variable is also valued continuously in [0, 1].Here, the threshold projection is used to remove the gray regions and control the minimum length scale in the derived pattern.The surface-PDE filter for the design variable of the pattern is implemented by solving the following surface-PDE [50]: where γ is the design variable; γ f is the filtered design variable; r f is the filter radius, and it is constant; ∇Γ and divΓ are the tangential gradient operator and tangential divergence operator defined on the implicit 2-manifold Γ, respectively; nτ Γ = nΓ × τ Γ is the outward unitary conormal vector normal to ∂Γ and tangent to Γ at ∂Γ, with nΓ and τ Γ representing the unitary normal vector on Γ and the unitary tangential vector at ∂Γ, respectively.The threshold projection of the filtered design variable is implemented as [51,52] where β and ξ are the parameters for the threshold projection, with values chosen based on numerical experiments [52].
Projection in Eq. 10 Surface-PDE filter in Eq. 8 The variational formulation of the surface-PDE filter is considered in the first order Sobolev space defined on Γ.It can be derived based on the Galerkin method as where γf is the test function of γ f ; H (Γ) represents the first order Sobolev space defined on Γ; L 2 (Γ) represents the second order Lebesque space defined on Γ.

Coupling of design variables
The design variable introduced in Section 2.2.2 for the pattern of the surface flow is defined on the implicit 2-manifold introduced in Section 2.2.1.Their coupling relation can be derived by transforming the tangential gradient operator ∇Γ, the tangential divergence operator divΓ and the unitary normal nΓ into the forms defined on the base manifold Σ.The transformation of the tangential gradient operator ∇Γ is implemented based on the following relation: where PΓ is the normal projector on the tangential space of Γ.The unitary normal vector nΓ is transformed as n (df ) where • 2 is the 2-norm of a vector.In Eq. 12, the transformed unitary normal vector is distinguished from the original form by using the filtered design variable d f as the superscript, and this identification method is used in the following for the other transformed operators and variables.The normal projector PΓ is sequentially transformed as where I is the two-dimensional unitary tensor; and the superscript T represents the transposition operation of a vector or tensor.The tangential gradient operator ∇Γ can then be transformed as Based on the transformed tangential gradient operator, the tangential divergence operator divΓ can be transformed as div (df ) where tr is the trace operator used to extract the trace of a tensor.Because the tangential gradient operator ∇Γ depends on d f , its first-order variational to d f can be derived as Similarly, the first-order variational of divΓ to d f can be derived as Because d f is differential homeomorphism, it can induce the Riemannian metric.Then, the differential on the base manifold and implicit 2-manifold satisfies where dl ∂Γ and dl ∂Σ are the differential of the boundary curves of Γ and Σ, respectively.Based on the transformed tangential gradient operator in Eq. 14 and the homeomorphism between H (Γ) and H (Σ) described in Eq. 5, the coupling relation between the two sets of design variables can be derived by instituting Eq. 14 into Eq.10: where the tangential gradient operator ∇Γ on Γ is replaced to be its transformed form ∇ (df ) Γ in Eq. 14.

Surface Navier-Stokes equations defined on implicit 2-manifold
The governing equations for the motion of a Newtonian surface fluid can be formulated intrinsically on a 2-manifold of codimension one in an Euclidian space.Based on the conservation laws of momentum and mass, the surface Navier-Stokes equations can be derived to describe the incompressible surface flows [53][54][55]: where u is the fluid velocity; p is the fluid pressure; ρ is the fluid density; η is the dynamic viscosity; u • nΓ = 0 is the tangential constraint of the fluid velocity.The tangential constraint is imposed, because the fluid spatially flows on the 2-manifold Γ and the fluid velocity is a vector in the tangential space of Γ.
To solve the surface Navier-Stokes equations, the fluid velocity and pressure are required to be specified at some boundaries, interfaces or points of the 2-manifold Γ: where u l v,Γ is a known distribution of the fluid velocity depending on the specified fluid velocity u l v,Σ at lv,Σ representing a boundary or interface curve of Σ; lv,Γ satisfies lv,Γ ⊂ ∂Γ when lv,Γ is a boundary curve of Γ, and it satisfies lv,Γ ⊂ Γ when lv,Γ is an interface curve of Γ; ls,Γ is the boundary curve with open boundary condition, and it satisfies ls,Γ ⊂ ∂Γ; p0,Γ is the known fluid pressure depending on the specified fluid pressure p0,Σ at PΣ representing a finite point set on Σ; PΓ is a finite point set on Γ.In Eq.21, when the known fluid velocity is 0, the inlet or interfacial boundary condition degenerates into the no-slip boundary condition: where u l v,Γ is equal to 0 on lv0,Γ ⊂ lv,Γ, and lv0,Γ is the no-slip part of the boundary curve.
The variational formulation of the surface Navier-Stokes equations is considered in the functional spaces without containing the tangential constraint of the fluid velocity.The tangential constraint of the fluid velocity is imposed by using the Lagrangian multiplier [56,57].Based on the Galerkin method, the variational formulation of the surface Navier-Stokes equations can be derived as where λ is the Lagrange multiplier used to impose the tangential constraint of the fluid velocity; ũ, p and λ are the test functions of u, p and λ, respectively.The Lagrangian multiplier in Eq. 23 is used to impose the tangential constraint of the fluid velocity and acts as a distributed force in the normal direction of Γ.Such distributed force cancels out the centrifugal, Coriolis and Euler forces exerted on the fluid particles in the normal direction of Γ to satisfy the tangential constraint.
Because lv,Γ is homeomorphous to lv,Σ, PΓ is homeomorphous to PΣ and lv,Σ and PΣ are fixed, u l v,Γ and p0,Γ are homeomorphous to u l v,Σ and p0,Σ, respectively.Then, based on the coupling relations in Section 2.2.3, the variational formulation in Eq. 23 can be transformed into the form defined on the base manifold Σ:

Design objective in general form
The design objective of the fiber bundle topology optimization problem for the surface flow is considered in the following general form: where A and B are the integrands of the design objective.Based on the coupling relations in Section 2.2.3, the design objective in Eq. 25 can be transformed into the form defined on the base manifold Σ: In Eq. 26, the unitary tangential vector τ Γ at ∂Γ satisfies the relation sketched in Fig. 4: Eq. 26 can then be transformed into Based on the transformed design objective, the adjoint analysis of the fiber bundle topology optimization problem can then be implemented on the functional spaces defined on the base manifold.

Fiber bundle topology optimization problem
The fiber bundle of the surface flow is composed of the base manifold together with the implicit 2-manifold and the pattern, where Σ is the base manifold and Γ × [0, 1] is the fiber, respectively.It can be expressed as ) with the diagram shown in Fig. 5, where proj1 is the natural projection proj1 Based on the above introduction, the fiber bundle topology optimization problem of the surface flow can be constructed as Find where J0 is the value of the design objective corresponding to the initial distribution of the design variables; to regularize this optimization problem, area and volume constraints are imposed on the pattern of the surface flow and implicit 2-manifold, respectively; s is the area fraction of the pattern of the surface flow; v is volume fraction of spacial domain enclosed by the implicit 2-manifold and the base manifold; s0 ∈ (0, 1) and v0 ∈ [0, 1) are the specified area and volume fractions, respectively.The coupling relations among the variables, functions, tangential divergence operator and tangential gradient operator in the fiber bundle topology optimization problem are illustrated by the following arrow chart: Eq. 9

−−−→ γp
Eq. 29 where the design variables dm and γ, marked in blue, are the inputs; the design objective J, the area fraction s and the volume fraction v, marked in red, are the outputs.

Adjoint analysis
The fiber bundle topology optimization problem in Eq. 29 can be solved by using a gradient information-based iterative procedure, where the adjoint sensitivities are used to determine the relevant gradient information.The adjoint analysis is implemented for the design objective and the area and volume constraints to derive the adjoint sensitivities.The details for the adjoint analysis have been provided in the appendix in Section 7.
Based on the continuous adjoint analysis method [58], the adjoint sensitivity of the design objective J is derived as where γ f a and d f a are the adjoint variables of the filtered design variables γ f and d f , respectively; δ is the first-order variational operator.The adjoint variables can be derived from the adjoint equations in the variational formulations.The variational formulation for the adjoint equations of the surface Naiver-Stokes equations is derived as where ua, pa and λa are the adjoint variables of u, p and λ, respectively; ũa, pa and λa are the test functions of ua, pa and λa, respectively.The variational formulations for the adjoint equations of the surface-PDE filters for γ and dm are derived as and where γfa and dfa are the test functions of γ f a and d f a , respectively.For the area constraint, the adjoint sensitivity of the area s is derived as In Eq. 34, the adjoint sensitivity δ (s |Γ|) can be derived based on the adjoint analysis of s |Γ| = Γ γp dΓ: In Eq. 35, the adjoint variables γ f a and d f a are derived by solving the variational formulations for the adjoint equations of the surface-PDE filters for γ and dm, respectively: and The adjoint sensitivity δ |Γ| in Eq. 34 can be derived based on the adjoint analysis of |Γ| = Γ 1 dΓ: In Eq. 38, the adjoint variables γ f a and d f a are derived by solving the variational formulations for the adjoint equations of the surface-PDE filters for γ and dm, respectively: and For the volume constraint, the adjoint sensitivity of the volume v is derived as In Eq. 41, the adjoint variable d f a is derived by solving the variational formulation for the adjoint equation of the surface-PDE filter for dm: After the derivation of the adjoint sensitivities in Eqs. 30, 34 and 41, the design variables γ and dm can be evolved iteratively to determine the fiber bundle of the surface flow.

Numerical implementation
The fiber bundle topology optimization problem in Eq. 29 is solved by using an iterative procedure described as the pseudocode in Tab. 1, where a loop is included for the iterative solution.The surface finite element method is utilized to solve the variational formulations of the relevant PDEs and adjoint equations.On the details for the surface finite element solution, one can refer to [59].Especially, when the surface finite element method is used to solve the surface flow problems on the implicit 2-manifold filled with the porous medium, the Lagrange multiplier method is used to enforce the tangential constraints of the fluid velocity [56,57].To avoid the numerical singularity caused by the null value of the denominator, the 2-norm of a vector function is approximated in the numerical implementation as , where f is the vector function and 0 is the value of floating point precision.
To ensure the well-posedness of the numerical solution of the variational formulations of the surface Navier-Stokes equations and their adjoint equations (Eqs.23 and 31), the Taylor-Hood elements satisfying the inf-sup condition are used [60].The linear elements are used to interpolate the design variable of the pattern of the surface flow and solve the variational formulation of the surface-PDE filter for this design variable and the corresponding adjoint equation.The quadratic elements are used to interpolate the design variable for the implicit 2-manifold and solve the variational formulation of the surface-PDE filter for this design variable and the corresponding adjoint equation.The meshes of the Taylor-Hood, linear and quadratic elements of the quadrangular-element based discretization of the base manifold have been sketched in Fig. 6, including the mapping meshes on the implicit 2-manifold.
Algorithm: iterative solution of Eq. 29 Set u lv p 0 , ρ, η, A d , v 0 , and s 0 ; Set Table 1: Pseudocode used to solve the fiber bundle topology optimization problem for the surface flow.In the iterative solution loop, n i is the loop-index; n max is the maximal value of n i ; J ni−m and J ni−(m+1) are the values of J in the (n i − m)-th and (n i − (m + 1))-th iterations; and mod is the operator used to take the remainder.
In the iterative procedure, the projection parameter β with the initial value of 1 is doubled after every 30 iterations; the loop is stopped when the maximal iteration number is reached, or if the averaged variation of the design objective in continuous 5 iterations and the residuals of the area and volume constraints are simultaneously satisfied.The design variable is updated by using the method of moving asymptotes [61].

Results and discussion
In this section, the fiber bundle topology optimization is carried out for the surface flows defined on several different base manifolds, including the flat surfaces for the bending channel and the four-terminal device, the curved surfaces deformed from a square to a sphere and the ones deformed from a cylinder to a Möbius.
The design objective is set to be the combination of the power of the viscous dissipation and pressure drop between the inlet and outlet: where ω is the weight of the viscous dissipation and it is valued to be 9/10 and the weight of the pressure drop is hence 1/10.The density and dynamic viscosity of the fluid are assigned to be unitary.The surface flows are driven by the boundary velocity at the inlets, in the forms of parabolic distribution as the functions of arc length and with the magnitude set to be U0 with U0 = sup ∀x∈l v,Γ u l v,Γ 2 .The outlets are set to be open boundaries.The remained boundaries are in the type of no slip.

Bending channel
For the bending channel, the fiber bundle topology optimization is implemented on the flat surface Σ composed of the design domain ΣD and the fluid domain ΣF as shown in Fig. 7.For different values of the magnitude parameter A d , the optimized fiber bundles and their components are derived as shown in Fig. 8(a∼h) including the distribution of the velocity vectors, where the area and volume fractions are set to be s0 = 0.3 and v0 = 0, respectively.Especially, the fiber bundle topology optimization problem degenerates into the topology optimization problem for the bending flow on the flat surface as shown in Fig. 8(a), when the magnitude parameter A d is set to be 0. By setting the volume fraction to be 0, the implicit 2-manifold Γ is derived with the same absolute values of the positive part and negative part of the enclosed volumes at the two sides of Σ, respectively.This can be confirmed from the distribution of the filtered design variable d f shown in Fig. 8(a1, b1, c1, d1, e1, f1, g1 and h1).The objective values for the optimized results derived in Fig. 8 have been listed in Tab. 2. From Tab. 2, it can be concluded that higher value of A d is helpful to decrease the viscous dissipation and pressure drop of the bending flow, because the design space of the fiber bundle topology optimization problem in Eq. 29 can be enlarged by increasing the value of the magnitude parameter in Eq. 4.
The convergent histories of the objective values and area and volume constraints have been plotted in Fig. 9 for the results in Fig. 8c with the magnitude parameter A d = 2, including the snapshots for the evolution of the fiber bundles.From the convergent histories, the robust convergence of the numerical solution of the fiber bundle topology optimization problem can be confirmed for the bending flow.In the convergent histories, there are jumps of the objective values and area constraints, and those phenomena are caused by updating the projection parameter β in Eq. 9.Meanwhile, the convergent histories of the volume constraints are smooth, without jumping phenomenon.This is because that no projection operation is imposed to regularize the design variable of the implicit 2-manifolds.By choosing the magnitude parameter to be A d = 2, the fiber bundle topology optimization problem is further investigated for different values of the velocity magnitude at the inlet of the bending flow.The optimized results are derived as shown in Fig. 10 by setting U0 to be the elements of 1 × 10 0 , 2 × 10 2 , 5 × 10 2 , 8 × 10 2 , sequentially.Because larger value of U0 corresponds to stronger Reynolds effect of the surface flow, increasing the velocity magnitude at the inlet can strengthen the convection of the surface flow.Therefore, different fiber bundles are derived as shown in Fig. 10.Table 3: Values of the design objective in Eq. 43 for the fiber bundles in Fig. 10.The optimized entries have been noted in bold.

Four-terminal device
For the fiber bundle topology optimization problem with the flat surface as its base manifold, the four-terminal device is further investigated by setting the magnitude parameter A d to be 0 and 2, respectively.The computational domain is set as the flat surface shown in In Tab. 4, the four-terminal device has the topology of double bending channels, when the surface flow has relatively weak Reynolds effect; the topology changes into double straight channels, as the Reynolds effect is strengthened.From the comparison of the results corresponding to A d = 0 and A d = 2, it can be concluded that the increase of the magnitude parameter can speed up the change of the optimized topology of the four-terminal device, when the Reynolds effect is strengthened.This is because that the increase of the magnitude parameter enlarges the design space of the four-terminal device, i.e., the characteristic size and area of the channels are increased and the averaged velocity is decreased at the inlets, then the gradient of the velocity decreases, and hence the viscous dissipation and pressure drop decrease.

Flows on deformed surfaces
For continuously deformed base manifolds with keeping the conservation of area, the fiber bundle topology optimization is implemented for the surface flows.By deforming a square to a sphere as shown in Fig. 12(a1∼c1), the optimized fiber bundles are derived as shown in Fig. 12(a2∼a4), 12(b2∼b4) and 12(c2∼c4) including the distribution of the velocity vectors, where the area and volume fractions and the magnitude parameter are set to be s0 = 0.4, v0 = 0 and A d = 2, respectively.
In Fig. 12, the fiber bundle for the surface flow on the square is composed of the base manifold together with the pattern of the flat diffuser and the implicit 2-manifold coinciding with the base manifold.The flat diffuser is consistent with the previously reported results derived by using topology optimization [10].When the square deforms into the shape of a semi-sphere, the pattern of the surface flow spits into two branches; and the implicit 2-manifold shrinks to straighten the channels corresponding to the pattern of the surface flow.When the semi-sphere further deforms into a sphere, the two branches merges to remove one part of the no-slip boundary, and the surface flow evolves into the enclosed mode with two vortexes.The underlying mechanism for the evolution of the fiber bundles along with the deformation of base ( manifolds is that the fluid is prone to moving in the short path and widening the channel, and detachment of the no-slip boundary can help to decrease the viscous dissipation and pressure drop.Additionally, the viscous dissipation and pressure drop of the surface flows with the patterns in the optimized fiber bundles decreases, along with the base manifold deforming from a square to a sphere.This can be conformed from the pressure distribution in Fig. 13.Further, the fiber bundle topology optimization is implemented on the base manifolds derived by deforming a cylinder firstly to a strip and then to a Möbius as shown in Fig. with the sizes marked on the strip, to minimize the viscous dissipation and pressure drop for the surface flows.The area and volume fractions and magnitude parameter are set without change.The patterns of the surface flows and the implicit 2-manifolds of the optimized fiber bundles are derived as shown in Figs.14(a2∼a4), 14(b2∼b4), 14(c2∼c4), 14(d2∼d4) and 14(e2∼e4), including the distribution of the velocity vectors represented by the red arrows.Especially, the derived implicit 2-manifold on the Möbius is broken by setting the inlet simultaneously to be the outlet with the same known velocity distribution.The destination of such setting is to make the derived implicit 2-manifold be orientable and remove the singularity of the normal direction on the non-orientable Möbius.
In Fig. 14, the derived pattern and implicit 2-manifold defined on a cylinder is composed of a circle channel and a curved surface with asymmetry.Then, the cylinder is opened and sequentially evolved into the shapes of semi-cylinder, strip, semi-Möbius until being enclosed again into the shape of Möbius, where the conservation of area is kept.The derived patterns and implicit 2-manifolds defined on the base manifolds corresponding to the sequential evolution of the opened cylinder are all composed of channel-shaped patterns and curved surfaces with asymmetry.The asymmetry assists the derived implicit 2-manifolds to satisfy the volume constraint with the volume fraction of 0, and it is caused by the asymmetrical property of the convection of the surface flows.The asymmetry is advantageous to shorten the path of the surface flows, then to decrease the viscous dissipation and pressure drop.
The pressure distribution in the derived fiber bundles has been provided in Fig. 15, which shows that the deformation of the base manifold from the cylinder to the Möbius is advantageous to decrease the viscous dissipation and pressure drop.This is because that the deformation can shorten the path of surface flows with the optimized patterns on the derived implicit 2-manifolds.Additionally, the case of Möbius has similar pressure drop to that of cylinder.Therefore, the orientability of 2-manifolds can provide similar effectivity on the design domain for fiber bundle topology optimization to minimize the viscous dissipation and pressure drop of a surface flow.

Conclusions
A fiber bundle topology optimization approach for the surface flow has been developed to match the implicit 2-manifold and the pattern defined on it, where the surface flow is described by the surface Navier-Stokes equations defined on the implicit 2-manifold.The material distribution method is used to implement the topology evolution of the pattern, where an artificial Darcy friction force of the porous model is added to the surface Navier-Stokes equations.The implicit 2-manifold is evolved based on the homeomorphous map between it and the base manifold.Continuous adjoint analysis method has been used to analyze the fiber bundle topology optimization problem.
Numerical tests have been presented to demonstrate this approach, including the fiber bundle topology optimization of bending channel, four-terminal device and fluid channels on continuously deformed base manifolds.The bending channel has been optimized to present the effect of the magnitude parameter used to determine the design space of the implicit 2-manifold, and the results show that the increase of the magnitude parameter can enlarge the design space of the fiber bundle for the surface flow.The Reynolds effect has been demonstrated by the  (  fiber bundle topology optimization of the bending channel and four-terminal device, by setting different velocity magnitude at the inlets.In the results of the bending channel, the valley and slope-shaped implicit 2-manifolds are derived to shorten the fluid path.In the results of the four-terminal device, the magnitude parameter of the implicit 2-manifold can speed up the topology change from double bending channels to double straight channels as the Reynolds effect is strengthened, compared with the results of topology optimization for the surface flow on the flat surface corresponding to the degenerated case with the null value of the magnitude parameter.The fiber bundle topology optimization has also been implemented on the base manifolds derived by deforming a square into a sphere and deforming a cylinder into a Möbius, where the area conservation is kept during the deformation.The derived results show that the non-orientable base manifolds have similar performance to the orientable ones on shortening the fluid path and minimizing the viscous dissipation and pressure drop.The presented fiber bundle topology optimization approach includes the design domain into the design space of fluidic structures.This approach achieves the topology optimization for fluid flows on the variable design domain.It provides a topology optimization method for the conformal design of fluidic channels, where the channel topology and the outer shapes of structural walls can be optimized simultaneously to achieve the matching optimization.Especially, the fiber bundle topology optimization problem will degenerate into the topology optimization problem for the fluid flow on a flat surface, if the null value is chosen for the magnitude parameter of the surface-PDE filter and the flat surface is set as the base manifold.Therefore, the presented fiber bundle topology optimization is the generalization of topology optimization for two-dimensional flow problems.This paper focuses on the laminar surface flows.In the future, it can be promoted for the turbulent surface flows.

Acknowledgements
The authors acknowledge the support of the National Natural Science Foundation of China (No. 51875545), the Innovation Grant of Changchun Institute of Optics, Fine Mechanics and Physics (CIOMP), the Youth Innovation Promotion Association of the Chinese Academy of Sciences (No. 2018253) and the Fund of State Key Laboratory of Applied Optics (SKLAO).They are also grateful to Prof. K. Svanberg of KTH for supplying the codes for the method of moving asymptotes.

Appendix
This section provides the details for the adjoint analysis of the fiber bundle topology optimization problem in Eq. 29.

Adjoint analysis for design objective
Based on the transformed design objective in Eq. 28, the variational formulations of the surface-PDE filters in Eqs.7 and 10 and the surface Navier-Stokes equations in Eq. 24, the augmented Lagrangian of the design objective in Eq. 29 can be derived as Based on the transformed operators in Eqs. 14 and 15 and their first order variationals in Eqs.16 and 17, together with the first order variational of the 2-norm of a vector function with f representing the vector function, the first order variational of the augmented Lagrangian in Eq. 44 can be derived as with the satisfication of the constraints in Eq. 45 and According to the Karush-Kuhn-Tucker conditions of the PDE constrained optimization problem [58], the first order variational of the augmented Lagrangian to the variables u, p and λ can be set to be zero as the first order variational of the augmented Lagrangian to the variable γ f can be set to be zero as and the first order variational of the augmented Lagrangian to the variable d f can be set to be zero as The constraints in Eqs.Without losing the arbitrariness of δu, δp, δλ, δγ f , δd f , δγ and δdm, one can set δu = ũa with ∀ũa ∈ (H (Σ)) 3 , δp = pa with ∀pa ∈ H (Σ), δλ = λa with ∀ λa ∈ L 2 (Σ), δγ f = γfa with ∀γ f a ∈ H (Σ), δd f = dfa with ∀ dfa ∈ H (Σ), δγ = γ with ∀γ ∈ L 2 (Σ) and δdm = dm with ∀ dm ∈ L 2 (Σ), to derive the adjoint system composed of Eqs. 30, 31, 32 and 33.

Adjoint analysis for area constraint
Based on the variational formulations of the surface-PDE filters in Eqs.7 and 10, the augmented Lagrangian of the pattern area s |Γ| can be derived as Based on the transformed operators in Eq. 14 and its first order variational in Eq. 16, the first order variational of ' s |Γ| can be derived as According to the Karush-Kuhn-Tucker conditions of the PDE constrained optimization problem [58], the first order variational of the augmented Lagrangian to the variable γ f can be set to be zero as and the first order variational of the augmented Lagrangian to the variable d f can be set to be zero as Based on the variational formulations of the surface-PDE filters in Eqs.7 and 10, the augmented Lagrangian of the area of the implicit 2-manifold can be derived as Based on the transformed operators in Eq. 14 and its first order variational in Eq. 16, the first order variational of |Γ| can be derived as According to the Karush-Kuhn-Tucker conditions of the PDE constrained optimization problem [58], the first order variational of the augmented Lagrangian to the variable γ f can be set to be zero as and the first order variational of the augmented Lagrangian to the variable d f can be set to be zero as Without losing the arbitrariness of δγ f , δd f , δγ and δdm, one can set δγ f = γfa with ∀γ f a ∈ H (Σ), δd f = dfa with ∀ dfa ∈ H (Σ), δγ = γ with ∀γ ∈ L 2 (Σ) and δdm = dm with ∀ dm ∈ L 2 (Σ), to derive the adjoint systems composed of Eqs.35,36,37 and Eqs. 38,39,40, respectively.Then, the adjoint sensitivity of the area fraction s can be derived from Eq. 34.

Adjoint analysis for volume constraint
Based on the variational formulations of the surface-PDE filter in Eq. 7, the augmented Lagrangian of the volume fraction v can be derived as The first order variational of v can be derived as According to the Karush-Kuhn-Tucker conditions of the PDE constrained optimization problem [58], the first order variational of the augmented Lagrangian to the variable d f can be set to be zero as Without losing the arbitrariness of δd f and δdm, one can set δd f = dfa with ∀γ f a ∈ H (Σ), and δdm = dm with ∀ dm ∈ L 2 (Σ), to derive the adjoint system composed of Eqs.41 and 42.

Figure 1
Figure1: Sketch for the fiber bundle of a surface flow, where Σ is the base manifold, Γ is the implicit 2-manifold used to define the pattern of the surface flow, γ p : Γ → [0, 1] is the pattern of the surface flow, u is the fluid velocity of the surface flow, u lv,Γ is the known fluid velocity at the boundary of Γ, n Γ is the unitary normal vector of Γ, τ Γ is the unitary tangential vector at ∂Γ, n τ Γ = n Γ × τ Γ is the outward unitary normal at ∂Γ, x Σ denotes a point on Σ, and x Γ denotes a point on Γ.This paper focuses on the laminar surface flows with low and moderate Reynolds numbers to demonstrate the fiber bundle topology optimization approach, although the sketched surface flows can be turbulent with high Reynolds number.

Figure 3 :
Figure 3: Sketch for the surface-PDE filter and projection operation for the design variable of the pattern of the surface flow.

Figure 4 :
Figure 4: Sketch for relation among the unitary tangential vector τ Γ at ∂Γ, the unitary normal vector n Σ on Σ and the tangential gradient ∇ Σ d f .

Figure 5 :
Figure 5: Diagram for the fiber bundle composed of the base manifold, the implicit 2-manifold and the pattern of the surface flow.

Figure 6 :
Figure 6: Sketch for the meshes of the Taylor-Hood, linear and quadratic elements of the quadrangular-element based discretization of the base manifold Σ and the mapping meshes on the implicit 2-manifold Γ.

Figure 7 :
Figure 7: Sketch for the base manifold for the fiber bundle topology optimization of the bending channel, where Σ D is the design domain and Σ F is the channel domains.

Figure 8 :
Figure 8: Fiber bundle topology optimization of the bending channel with different values of the magnitude parameter A d , including the distribution of the filtered design variables, the pattern of the bending flow projected on the base manifold and the fiber bundle composed of the base manifold together with the implicit 2-manifold and pattern of the surface flow, where the red arrows represent the distribution of the fluid velocity.(a1)∼(a3) are the results for A d = 0; (b1)∼(b3) are the results for A d = 1; (c1)∼(c3) are the results for A d = 2; (d1)∼(d3) are the results for A d = 3; (e1)∼(e3) are the results for A d = 4; (f1)∼(f3) are the results for A d = 5; (g1)∼(g3) are the results for A d = 6; (h1)∼(h3) are the results for A d = 7.

Figure 10 :
Figure 10: Fiber bundle topology optimization of the bending channel with different values of the velocity magnitude U 0 at the inlet, including the distribution of the filtered design variable d f , the pattern of the bending channel γ p projected on Σ and the fiber bundle composed of the pattern and 2-manifolds, where the red arrows represent the distribution of the fluid velocity.(a1)∼(a3) are the results for U 0 = 1 × 10 0 ; (b1)∼(b3) are the results for U 0 = 2 × 10 2 ; (c1)∼(c3) are the results for U 0 = 5 × 10 2 ; (d1)∼(d3) are the results for U 0 = 8 × 10 2 .

Figure 11 :
Figure 11: Sketch for the base manifold for the fiber bundle topology optimization of the fourterminal device, where Σ D is the design domain and Σ F is the channel domains.

U 0 = 5 . 50 × 10 2 U 2 A d = 0 A d = 2 Table 4 :
0 = 5.75 × 10 2 U 0 = 6.25 × 10 Fiber bundle topology optimization of the four-terminal device for different values of the velocity magnitude U 0 at the inlets, where the red arrows represent the distribution of the fluid velocity.

Figure 12 :
Figure 12: Fiber bundle topology optimization for the surface flows on the base manifolds deformed from a square to a sphere with keeping the conservation of area, where the base manifolds are sketched in Figs.12(a1), 12(b1) and 12(c1), the distribution of the filtered design variables for the implicit 2-manifolds are shown in Figs.12(a2), 12(b2) and 12(c2), the projected patterns on the base manifolds are shown in Figs.12(a3), 12(b3) and 12(c3), and the fiber bundles are derived as shown in Figs.12(a4), 12(b4) and 12(c4) with the red arrows representing the distribution of the fluid velocity.

Figure 13 :
Figure 13: Distribution of the fluid pressure for the surface flows on the fiber bundles derived as shown in Figs.12(a4), 12(b4) and 12(c4).

2 −
45 and 48 are imposed to Eq. 49.Further, the adjoint sensitivity of J A d d f a δdm dΣ.(52)