Topology Optimization for Steady-State Navier-Stokes Flow Based on Parameterized Level Set Based Method

Peng Wei,Zirun Jiang,Weipeng Xu,Zhenyu Liu,Yongbo Deng,3 and Minqiang Pan

1State Key Laboratory of Subtropical Building Science,School of Civil Engineering and Transportation,South China University of Technology,Guangzhou,510641,China

2State Key Laboratory of Applied Optics(SKLAO),Changchun Institute of Optics,Fine Mechanics and Physics(CIOMP),Chinese Academy of Sciences,Changchun,130033,China

3Institute of Microstructure Technology(IMT),Karlsruhe Institute of Technology(KIT),Eggenstein-Leopoldshafen,76344,Germany

4School of Mechanical and Automotive Engineering,South China University of Technology,Guangzhou,510641,China

ABSTRACT In this paper,we consider solving the topology optimization for steady-state incompressible Navier-Stokes problems via a new topology optimization method called parameterized level set method,which can maintain a relatively smooth level set function with a local optimality condition.The objective of topology optimization is to find an optimal configuration of the fluid and solid materials that minimizes power dissipation under a prescribed fluid volume fraction constraint.An artificial friction force is added to the Navier-Stokes equations to apply the no-slip boundary condition.Although a great deal of work has been carried out for topology optimization of fluid flow in recent years, there are few researches on the topology optimization of fluid flow with physical body forces.To simulate the fluid flow in reality, the constant body force (e.g., gravity) is considered in this paper. Several 2D numerical examples are presented to discuss the relationships between the proposed method with Reynolds number and initial design,and demonstrate the feasibility and superiority of the proposed method in dealing with unstructured mesh problems.Three 3D numerical examples demonstrate the proposed method is feasible in threedimensional.

KEYWORDS Topology optimization;parameterized level set method;power dissipation;no-slip boundary condition;Navier-Stokes equations;gravity

1 Introduction

The topology optimization method originating from the field of structural mechanics,has received much attention in recent decades.As an advanced design means,the topology optimization method has been widely implemented in the structure design of aerospace, automotive, and construction engineering [1]. Compared with size optimization and shape optimization, topology optimization can not only change the structural boundaries but also change the layout of materials to obtain the best structures or components during the optimization procedure. Up to the present, the topology optimization method mainly covers the homogenization method [2,3], the density method, such as Solid Isotropic Material with Penalization(SIMP)[4,5],evolutionary structural optimization(ESO)/bi-directional evolutionary structural optimization(BESO)[6-8],the level set method[9-12],moving morphable component/void method [13,14], and so on. In recent years, the topology optimization method has been introduced into the fluid field to guide the design of fluid flow and heat conduction paths[15-18].

A lot of attention has been paid to solving the optimal design problem of fluid flow. In the 1970s,Pironnea[19,20]carried out shape optimization for incompressible viscous flow,a pioneering work that has received much attention. In 2003, Borrvall et al. [15] proposed a relaxed material distribution approach based on the density method to solve the topology optimization problem for Stokes flow, which anticipates minimizing the power dissipation with a prescribed fluid volume fraction.To distinguish the solid and fluid materials,an artificial friction force,which is proportional to the fluid velocity,was added to the Stokes equations.This approach was enlightening and worked well. Then, Gersborg-Hansen [16], Sigmund et al. [21], and Gersborg-Hansen et al. [17] extended this approach to the Navier-Stokes flow. In 2006, Olesen et al. [18] presented an efficient method for topology optimization for Navier-Stokes flow, in which topology optimization code was edited in the FEMLAB package to perform topology optimization for steady-state incompressible Navier-Stokes flow,and obtained optimal fluid flow paths that minimize power dissipation.Guest et al.[22]applied a new approach for the topology optimization of creeping fluid flows, in which the solid region is modeled as with Darcy flow of low permeability using an interpolated Darcy-Stokes finite element. Liu et al. [23] studied the topology optimization for Navier-Stokes flow with flow rate equality constraints, in which the lumped Lagrange multiplier method was used to implement the equality constraints on the specified boundaries.Deng et al.[24]studied the effects of dynamic inflow,Reynolds number, and target flux on specified boundaries for the optimal topology for unsteady Navier-Stokes flows. Deng et al. [25] used the density method to perform topology optimization of steady and unsteady incompressible Navier-Stokes flows driven by body forces,in which the artificial friction force and physical body forces are considered simultaneously. Pereira et al. [26] applied an educational code PolyTop [27] written in MATLAB to topology optimization problems for Stokes flow,using polygonal elements and the density method.Shen et al.[28]studied the three-phase(i.e.,solid,fluid,and porous materials)interpolation scheme using SIMP for the Navier-Stokes flow,aiming to minimize the pressure attenuation in multiple phases interpolation models. The density method is the most common and popular approach whether in structural or fluid topology optimization.However, the optimal configuration obtained by the density method contains a gray region, which makes it difficult to accurately describe the boundary between the fluid and solid materials.

The level set method has made great progress in many fields, such as image processing [29],computational geometry [30], structural optimization [31], and fluid dynamics [32,33]. In 1988, the level set method was first proposed by Osher et al.[34],and introduced into structural optimization by Sethian et al. [31], which has been greatly developed in the structural topology optimization field. The level set method can describe topology and shape changes by merging and breaking boundaries. Compared with SIMP and ESO/BESO, the main advantage of the level set method is that it always provides clear geometric boundaries during the optimization procedure, and thus has natural advantages in handling complex shape and topology changes[35].

The level set method has recently gained popularity in fluid topology optimization.Zhou et al.[36]presented a level set method for topology optimization of steady-state Navier-Stokes flow and constructed a variational form of the cost function based on the adjoint variable and Lagrangian multiplier technique.Duan et al.[37-39]proposed a new algorithm based on the variational level set method to investigate topology optimization problems for Stokes and Navier-Stokes flow,in which a relatively smooth evolution can be maintained without re-initialization.Challis et al.[40]used the level set method to solve the topology optimization problem for Stokes flow.To solve the Stokes equations only in the fluid portion of the design domain and achieve no-slip boundary conditions,the velocity DOF(degree of freedom)of the nodes touched by solid elements and the pressure DOF of the nodes surrounded by solid elements are removed from the matrix system, respectively. Pingen et al. [41]employed a parametric level-set description to study whether the convergence and the versatility of topology optimization methods for fluidic systems can be improved. Kreissl et al. [42] presented an explicit level set method that can avoid solving the Hamilton-Jacobi equation and allows using standard nonlinear programming methods.Kreissl et al.[43]studied the incompressible Navier-Stokes flow at low Reynolds numbers and the flow field was discretized by the extended finite element method.Deng et al. [44] used the level set method to solve the topology optimization problem of Navier-Stokes flow considering several kinds of body forces,the fluid and solid materials are distinguished by constraining the fluid velocity to be zero in solid material.Zhang et al.[45]employed the level set method to study the topology optimization problem of an arterial bypass in non-Newtonian flow governed by the Navier-Stokes equations. Deng et al. [46] used the level set method to study the topology optimization problem for two-phase flow. Koch et al. [47] proposed a novel transitional procedure,in which a parameterized solution has been created by post-processes 2D adjoint topology solutions which can be used as a CAD-compatible representation of the interface and a source for grid generation from which a shape optimization loop can be initialized.Dai et al.[48]proposed a piecewise constant level set method for topology optimization of both two-dimensional and three-dimensional steady-state incompressible Navier-Stokes flow. Duan et al. [49] presented a novel algorithm that belongs to the level set method but doesn’t capture the boundary explicitly as in the classical level set method and can achieve high efficiency.Nguyen et al.[50]attempted to improve the energy dissipation devices with vibration control and proposed a scheme for the topology optimization of transient flow using the level set method incorporating the lattice Boltzmann method. Kubo et al. [51] used a level-set boundary expression and immersed boundary method to study the total pressure drop minimization problems considering two-dimensional turbulent flows, under the frozen turbulence assumption.The topological derivative is used as design sensitivity and the level set function is updated by solving a reaction-diffusion equation. Cai et al. [52] proposed a derivative-free level-set based topology optimization method and use it for flow channel structure design, in which the optimized channel boundary was smoothed using the B-spline method.

However,few researches that had applied the parameterized level set method to study the topology optimization problems for fluid flow,which can solve the regular design domain problems discretized by structured mesh and complex irregular design domain problems discretized by unstructured mesh,respectively.Besides,the previous researches do not consider physical body forces which exist widely in reality and this is important for the optimized configuration.

In the conventional level set method, the level set function is updated generally by solving a Hamilton-Jacobi partial differential equation (PDE) using the finite difference method. However,the time step size of the Hamilton-Jacobi PDE should be sufficiently small because of the Courant-Friedrichs-Lewy(CFL)condition[53,54].To ensure a smooth evolution of the level set method,a reinitialization scheme must be implemented frequently in the conventional level set method[53,54].In this paper,we would like to parameterize the level set function using the basis functions to improve the conventional level set method.Wang et al.[55,56],and Wang et al.[57]firstly introduced radial basis functions(RBFs)into the level set method and proposed an effective parameterized level set method using RBFs for shape and topology optimization.Later,Wei et al.[35]proposed a further improved parameterized level set method using RBFs,which is proven to be effective in complex design domains[58],and chosen here to study the topology optimization of steady-state incompressible Navier-Stokes problems.Recently it has been extended to parallel computation versions in both structured mesh[59]and unstructured mesh[60].Combined with parallel computing methods,the parameterized level set method can deal with larger-scale and more complex topology optimization problems.

Theoretically,the essential idea of the parameterized level set method using RBFs is that the level set function is replaced with a linear combination of a set of RBFs and corresponding expansion coefficients,thus converting the evolution of the level set function into updating expansion coefficients of RBFs.A natural advantage of the parameterized level set method is that the RBFs can refactor any possible structure combined with the level set topology optimization method,but it still preserves the advantages of an implicit description and avoids the problems caused by grid discretization at the same time[58].

In the present work, we focus on the steady-state Navier-Stokes flow, which is an important physical model in channel and reactor design. An artificial friction force, which is proportional to the fluid velocity, is added to Navier-Stokes equations to distinguish the solid and fluid regions.The physical body forces exist widely in reality so the physical body forces should be considered to simulate real fluid flow.The constant body force(e.g.,gravity)is considered in this paper.Combining the parameterized level set method,the artificial friction force term,and the gravity, we construct a topology optimization problem for steady-state incompressible Navier-Stokes flow.

The remainder of this paper is organized as follows. In Section 2, a brief introduction to the parameterized level set method is given.In Section 3,the governing equations in the fluid are given,the topology optimization problem is proposed and sensitivity information is derived. In Section 4,the numerical implementation of the topology optimization algorithm is discussed in detail. In Section 5, several numerical examples are presented to illustrate the method proposed in this paper has effectiveness.In Section 6,some conclusions of our work are obtained.

2 The Parameterized Level Set Method

A brief introduction to the parameterized level set method is given in this part. In the level set method,the dynamic boundary of the design is implicitly described by the zero contour or isosurface of the higher one-dimensional level set function φ(x,t),which is usually defined as follows:

where φ(x,t)and x = (x,y)represent the level set function and the spatial coordinate of any point in the full design domain,respectively,and t is the pseudo time.Ωs,Ωland Γ denote the solid region,fluid region, and the boundary between the solid and fluid regions. Ωsand Ωlare open, and they satisfy Ωs∩Ωl=Γ.Fig.1 shows the level set function and the corresponding level set function values of the solid region,fluid region,and boundary for a 2-D case.

Figure 1:The level set function and the corresponding level set function values of the solid region,the fluid region,and the boundary for a 2-D case(Ωs∪Ωl∪Γ =Ω)

In the conventional level set method, the level set function is updated by solving the following Hamilton-Jacobi partial differential equation(PDE),which is written as follows:

where ?φ denotes the gradient of the level set function, andVnis the normal velocity of the design boundary. It is important to note that the term |?φ| in (2) is replaced by the term δ(φ) here in this paper[35,55],so(2)can be modified to:

The essence of the parameterized level set method is to interpolate the level set function using a set of basis functions,to fit the structure into a higher one-dimensional smooth surface.For a given basis function,the level set function can be parameterized to the following expression:

where αk(t)is the expansion coefficient of the basis function corresponding to thekth knot,gk(x)is the basis function value of x coordinate for thekth basis function.The basis function includes RBFs(such as Gaussians,Inverse Multiquadric splines,Multiquadric splines,Compactly Supported Radial Basis Functions),B-spline functions,and shape functions in the finite element method(FEM).Fig.2 shows a schematic diagram of the level set function interpolated by three radial basis functions R1,R2,and R3,the corresponding expansion coefficient of R1,R2,and R3 is 1,-2.5,and 2,respectively.

In the parameterized level set method, the level set function φ(x,t) is interpolated by basis functions,which transforms the Hamilton-Jacobi partial differential equation(PDE)into an ordinary differential equation(ODE)updating the expansion coefficients[35]:

Theoretically, the matrix G is invertible [61],tirepresents theith time step and Δtdenotes the artificial time step size.Then,the level set function can be updated by the following expression:

In the conventional level set method,a re-initialization scheme should be implemented frequently to ensure a smooth evolution.In the work of this paper,to prevent|?φ|from being too large or too small, Wei et al. [35] adopted a simple approximate re-initialization scheme to update the level set function.The approximate re-initialization scheme can be defined as follows:

Figure 2:The schematic diagram of level set function interpolated by three radial basis functions R1,R2,and R3

3 Topology Optimization Formulation

In this section, several key formulas for topology optimization problems are given, including governing equations, objective function, constraint conditions, and sensitivity analysis. Firstly, the relationship between the level set function and the fluid-solid materials should be established.Fig.3 shows the scheme of the regions distinguished based on the level set function values.The impermeability is a medium that associates the level set function and the fluid-solid materials,which is zero in the fluid region and an infinite positive constant in the solid region.In this way,the fluid-solid results can be obtained without intermediate materials.Wei et al.[64]proposed a new method,called the level set band method,combining the density method and the level set method,to improve the continuity of objective and constraint functions.This method can deal with the problem of permeability distribution optimization well and is a promising research direction.

Figure 3:The scheme of the regions distinguished based on the level set function values(Ωs∪Ωl∪Γ =Ω)

3.1 The Governing Equations

The steady-state incompressible flow in the 2-D case is considered in the current work,which is governed by the Navier-Stokes equation and the continuity equation:

where η and ρ are the dynamic viscosity and the density of the fluid, respectively. u = (ux, uy) is the fluid velocity, andpdenotes the fluid pressure. The first and second above equations represent the momentum equation and the continuity equation,respectively.To distinguish the solid and fluid regions,an artificial friction force fais added into the momentum equation to implement the no-slip boundary condition[15].fais proportional to the fluid velocity which is expressed as follows:

where α(φ)can be defined as the impermeability which is expressed by the level set function:

where αminand αmaxrepresent impermeability in fluid and solid regions which are equal to zero and an infinite constant,respectively.However,to ensure numerical stability,αmaxshould be selected as a large but finite positive constant.H(φ)is the Heaviside function which is expressed as follows:

wherehis the support size,ε and α are small positive constants,respectively.

Because gravity will drive the fluid to flow and increase the fluid velocity in the solid region,gravity cannot be added to the artificial friction force directly[25].To avoid the fluidic flow driven by the gravity in the solid region,the gravity is expressed as follows:

where fbgis the gravity.Then,fbapproaches zero in the solid region;and approaches fbgin the fluid region.Therefore,the body force can be expressed as follows:

As general boundary conditions,Dirichlet boundary condition and Neumann boundary condition are defined as ΓDand ΓN,respectively.Then,the governing equations can be rewritten as:

where u0represents the known fluid velocity,which is a constant value,I denotes the identity matrix,n is the outward normal to the boundary, the body force term f is expressed as (14), and g = 0 corresponds to the open boundary condition.

3.2 The Topology Optimization Problem

In this paper, the topology optimization problem we want to solve is to find an optimal configuration of the solid and fluid materials,which can minimize the objective of power dissipation with the constraint of a prescribed fluid volume fraction. Similar problems have been proposed by Borrvall et al.[15],Gersborg-Hansen et al.[17],and Olesen et al.[18],and so on.The objective function of the topology optimization problem can be expressed as follows:

The objective function is governed by the Navier-Stokes equation and the continuity equation,and we consider the Dirichlet boundary condition and Neumann boundary condition on the inlet and outlet,respectively,the non-slip boundary condition on the wall,and the fluid volume constraint.Thus,the topology optimization problem here can be expressed as a constrained optimization problem:

where the body force term f is expressed as (14), consisting of artificial frictional force and gravity,andV0is the prescribed fluid volume fraction.

3.3 Sensitivity Analysis

In the level set method, the solid and fluid regions are distinguished by the zero contour or isosurface of the level set function. To obtain the normal velocity field driving the evolution of the boundary of the solid-fluid region,we need to conduct sensitivity analysis.

We use the intermediate variable γ (φ)to establish the connection between the nodal value of level set function φ and the local properties of the design domain. The intermediate variable γ (φ) varies from 0 to 1 and its expression is shown in(19).γ = 0 represents the solid region and γ = 1 the fluid region.

Under the framework of level set method,the level set function is updated by the normal velocityVn.Here we use the derivative of the objective Φ function with respect to the intermediate variable γ as the sensitivity information to update the level set function. This is quite similar to the sensitivity analysis procedure in density-based topology optimization design. Besides, a physically meaningful and efficient scheme called the velocity extension method[55]is adopted to extend the velocity to the entire design domain.The evolution of the level set function with the velocity extended method can lead to the creation of new holes inside the design domain,which can reduce dependency on the initial design.Therefore,the normal velocityVnis replaced by the extended velocityVEas follows:

where dΦ/dγ is the total derivative of the objective function with respect to the intermediate variable and λ is the Lagrange multiplier to deal with the constraint of volume fraction. When the objective function does not include the contribution of the boundary part as in this paper,the form of(20)is the same as the form of shape derivative in[44].

Based on the extended velocityVE,the level set function can be updated according to Eq.(6).And the updated nodal value can be used to calculate the latest state variable distribution in the design domain.

Because of the large number of nodal design variables,the adjoint method is adopted in this paper to calculate dΦ/dγ,just as Pingen et al.[41],Kreissl et al.[42,43],Duan et al.[65]did,and the same adjoint sensitivity algorithm can also be found in the COMSOL Multiphysics reference guide [66].The detailed derivation process of dΦ/dγ is given in the Appendix A.

4 Numerical Implementations

In this section,several details for the numerical application of topology optimization are discussed.The optimization procedure is shown in Fig.4,which includes the following steps:

1)The level set function φ is initialized;2)The intermediate variable γ is calculated;3)The fluid velocity u and the fluid pressurepare obtained by solving the governing equations; 4) Calculate the objective function and the sensitivity analysis is performed to calculate dΦ/dγ;5)The extended velocity fieldVEof the design boundary is obtained by sensitivity analysis;6)The level set function φ is updated;7)Check the convergence conditions:the optimization procedure is terminated when the difference between the current fluid volume and the prescribed fluid volume is less than 0.1%, and the difference between the current objective function value and the previous nine values are all less than 0.1%,respectively.The topology optimization procedure ends if the convergence conditions are reached, otherwise, the second step to the seventh step is repeated until the convergence conditions are achieved.

The proposed method in this paper is implemented based on numerical software MATLAB(Version 2008b) and commercial finite element software COMSOL Multiphysics (Version 3.5).COMSOL Multiphysics has significant advantages in dealing with multi-physical field problems and users can easily define desirable partial differential equations to solve different problems[67,68].

Figure 4:The flow chart of the optimization procedure

We edit codes in MATLAB, and then call the commands in COMSOL Multiphysics to solve the Navier-Stokes equation and adjoint equations. In the current work, the essential idea of the 88-line MATLAB code presented by Wei et al. [35] for the compliance minimization problem and the local optimality condition[62,63]are adopted here and extended to the fluid topology optimization problem. In addition, the design domain is discretized by the linear rectangle elements and triangle elements in 2D cases and tetrahedron elements in 3D cases.

5 Numerical Examples

In this section,several numerical examples are presented to demonstrate the effectiveness of the proposed method.The numerical examples in Section 5.1 do not consider gravity,and the body force f only includes the artificial friction-αu.In Section 5.2,artificial friction and gravity are considered to study the influence of gravity on the flow channels.For all numerical examples,the dynamic viscosity η and density ρ of the fluid are both set to 1. Dirichlet boundaries with the prescribed flow have parabolic flow profiles normal to the boundary with magnitude,which is given as follows:

whereDdenotes a characteristic length scale which is set to 1 in this paper.

5.1 Numerical Examples without Gravity

5.1.1 Double Pipe

The first numerical example is to study the topology optimization of a double pipe problem, in which the design domain is discretized by structured mesh and unstructured mesh,respectively.The design domain and boundary conditions information are shown in Fig.5.The prescribed fluid volume fraction of the optimization problem is 0.35, and the Reynolds number is set to Re = 10. One of the advantages of the parameterized level set method is to refactor any possible structure by basis functions combined with the level set topology optimization method. Therefore, the parameterized level set method has superiority in dealing with irregular design domain problems. The structured mesh is discretized by 7100 regular quadrilateral elements while the unstructured mesh is discretized by 6216 irregular triangular elements.Fig.6 displays the optimization history of the double pipe problem whose design domain is discretized by structured mesh.

Figure 5:Design domain and boundary conditions for the double pipe problem

Figure 6:Optimization history for the double pipe problem

The optimal solutions for different meshes of the double pipe problem are shown in Table 1.The optimal configurations in the first row of Table 1 display the coherence between the structured quadrilateral meshes and the unstructured triangular meshes.Due to the use of unstructured mesh,a small portion of the boundary is not sufficiently smooth,which may be further improved by using the polygonal elements [58,69]. However, the result for unstructured meshed has only slight differences from that of structured meshes. Further, Fluid velocity fields, fluid pressure fields, and the power dissipation values show negligible differences.By comparing the results of two different types of mesh,the method proposed in this paper is proven to apply to unstructured grids and works well.

Table 1: Comparison of solutions with structured and unstructured meshes for the double pipe problem

5.1.2 Four-Terminal Device with Different Reynolds Numbers

In the second numerical example, a four-terminal device is considered to minimize power dissipation with a prescribed fluid volume fraction which is similar to the problems proposed by Borrvall et al. [15], Gersborg-Hansen et al. [17], Olesen et al. [18] and Dai et al. [48]. The design domain and boundary conditions information are shown in Fig.7.The design domain is discretized by 100×70 rectangular elements and the fluid volume fraction is 0.4.To study the effect of the Reynolds number on the optimal configuration, the Reynolds number is set as Re = 20 and Re = 200,respectively.

Figure 7:Design domain and boundary conditions for the four-terminal device problem

Fig.8 displays the optimal configurations for different Reynolds numbers,in which the optimal four-terminal device has two bending channels for the flow with a small Reynolds number,while the optimal four-terminal device with a large Reynolds number has two parallel straight channels.With the increase of the Reynolds number,the fluid power dissipation caused by the bending channels increases.When the inertia effect dominates,a large velocity gradient appears in the bending channels,and this increases power dissipation compared to the low Reynolds number case.Therefore,the optimal fourterminal device has two parallel straight channels for large Reynolds number flow. Fig.9 displays the convergence curves of power dissipation and volume fraction when the Reynolds number is equal to 20.

Figure 8:Optimal configurations corresponding to different Reynolds numbers for the four-terminal device problem

Figure 9:Convergence curves of power dissipation and volume fraction when the Reynolds number is equal to 20 for the four-terminal device problem

In addition,Yaji et al.[70]have proven that the value of the Reynolds number affects the channel’s configuration,and the curvature radius of the channel decrease as the Reynolds number is increased.Furthermore,we study the double pipe problem in Section 5.1.1 at different Reynolds numbers,which is started from an initial region that is completely fluid.The prescribed fluid volume fraction of the optimization problem is 1/3,and the Reynolds numbers are set to Re = 20,Re = 200,and Re =360,respectively.The design domain and boundary conditions of the double pipe problem are shown in Fig.5, and the optimal double pipes for different Reynolds numbers are shown in Fig.10. The four-terminal device and the double pipe problem demonstrate that the Reynolds number affects the configuration of the flow channel,and the larger the Reynolds number,the more inclined to obtain straight flow channels.

Figure 10:Optimal configurations corresponding to different Reynolds numbers for the double pipe problem

Theoretically, when the Reynolds number is equal to 360, the optimal configuration should be two parallel straight channels, which minimizes power dissipation. However, when the initial design as shown in Fig.11a is adopted,two parallel straight channels can be obtained as shown in Fig.11b.The power dissipation in Figs.10c and 11b are 8.9648×106and 7.9960×106,respectively.Therefore,the proposed method is sensitive to the initial design for some problems.

5.1.3 Rectangular Splitter with Different Initial Designs

The third example is the topology optimization of a rectangular splitter problem,considering three different initial designs to examine the dependency of the method proposed in this paper on the initial design.Fig.12 shows the design domain and boundary conditions.The objective is to minimize the power dissipation in optimal configuration between the inlet and the outlets. The design domain of the rectangular splitter is discretized by 60×100 rectangular elements. The prescribed fluid volume fraction of the optimization problem is 0.5,and the Reynolds number is set to Re = 10 in all cases.Three cases start from completely fluid,3×3 solid holes and 4×5 solid holes,respectively.

Figure 11:Another initial design and corresponding optimal configuration when Reynolds number is equal to 360

Figure 12:Design domain and boundary conditions for the rectangular splitter problem

Table 2 exhibits three different initial designs and several iterative steps in their corresponding optimization procedures.Although the three initial designs are different,the optimal configurations obtained are roughly the same, and any two values of the objective function differ by less than 5%.Therefore, the method proposed in this paper is not sensitive to the initial design. In addition, the evolution histories of contours corresponding to three initial designs are shown in Table 2.So,here’s what we can see is optimal configurations have clear and smooth boundaries in the optimization procedure, which demonstrates that the method proposed in this paper can deal with geometric boundary and topology change problems well. Fig.13 gives the convergence curves of the objective function and the fluid volume fraction for case 3.

Table 2: Optimization procedures and several iterative steps for three different initial designs

5.1.4 A Z-Pipe

This numerical example is a Z-pipe problem that has an irregular design domain. The design domain and boundary conditions are shown in Fig.14a. The design domain is discretized by 4216 unstructured triangular elements and the Reynolds number is set to Re = 10.The objective function for this Z-pipe problem is still to minimize the power dissipation, with the prescribed fluid volume fraction being 0.3.As shown in Fig.14b,the optimal configuration is a curved channel between the inlet and outlet, with clear boundaries, although a small portion of the boundary is not sufficiently smooth. The example shows that the method proposed in this paper can be applied to the complex irregular design domain and a reasonable optimal configuration can be obtained. It also shows the unique advantage of the parameterized level set method in the topology optimization of irregular structures.

Figure 13:Convergence curves of power dissipation and volume fraction of case 3 for the rectangular splitter problem

Figure 14:Design domains,boundary conditions and optimal configuration for the Z-pipe problem

5.1.5 Bend Pipe in 3D

The first numerical example in 3D is a bend pipe problem which is used to demonstrate that the proposed method can be successfully applied to three-dimensional problems.The design domain of the bend pipe problem is shown in Fig.15, which is discretized by 161,851 tetrahedron elements.Note the types of boundary conditions for the inlets and outlets are consistent with the 2D numerical examples, but the initial velocities of the inlet are uniform flows rather than parabolic flows. The Reynolds number is set to Re = 10 and the fluid volume fraction is 10%.

Fig.16 shows the optimal configurations in different directions for the bend pipe problem.The optimal configuration is a bent pipe connecting the inlet and outlet. However, the surfaces of the optimal configurations are not smooth due to the unstructured mesh,i.e.,the tetrahedron element.A refinement method[71]is adopted here to obtain smooth results.Refined mesh schemes with 874,387 tetrahedron elements are adopted to smooth the surfaces of optimal configuration,the smooth results are shown in Fig.17. Therefore, this numerical example demonstrates that the proposed method is feasible in three-dimensional and smooth results can be obtained with a refinement method[71].

Figure 15:Design domain for the bend pipe problem

Figure 16:Optimal configurations in different directions for the bend pipe problem

Figure 17:Smooth optimal configurations in different directions for the bend pipe problem

5.1.6 Multi-Outlet Terminal in 3D

The second numerical example in 3D is a multi-outlet and its design domain is shown in Fig.18,which is discretized by 178,562 tetrahedron elements. The initial velocities of the inlet are uniform flows rather than parabolic flows. The Reynolds number is set to Re = 10, and the fluid volume fractions are 20%,15%,and 10%,respectively.

Figure 18:Design domain for the multi-outlet problem

Consistent with the bend pipe problem,the refinement method[71]is adopted to obtain smooth results.Refined mesh schemes with 2,192,849 tetrahedron elements are adopted to smooth the surfaces of optimal configuration,the smooth results corresponding to different volume fractions are shown in Table 3.With the decrease of volume fraction,the pipes become thinner.In addition,it can be seen that when the volume fraction is 10%, the pipe on the right side is thinner than those on the other sides,the reason may be the finite element mesh dividing the design domain is unstructured.

Table 3: Smooth optimal configurations corresponding to different volume fractions for the multioutlet problem

5.2 Numerical Examples with Gravity

5.2.1 Horizontal Channel

In this numerical example,the simplest body force,e.g.,gravity,is considered in the Navier-Stokes equations. Comparing the optimal fluid channels with and without gravity, the effectiveness of the proposed method in considering gravity is studied.The design domain and boundary conditions are shown in Fig.19.Note that the inlet is the middle third of the length on the left side while the outlet is the entire length on the right side.The design domain is discretized by 3600 rectangle elements and the Reynolds number is set to Re = 1.The objective function is still to minimize the power dissipation,with the prescribed fluid volume fraction being 1/3.The gravity fbgis set to be(0,-10),which is the downward direction.Topology optimization starts from completely fluid.

Figure 19:Design domain and boundary conditions for the horizontal channel in gravity problem

Fig.20 shows the Optimal configurations for the horizontal channel with and without gravity.It can be concluded that gravity causes the bending horizontal channel by comparing the optimal configuration without gravity.Extra work is applied to the fluid by gravity,causing the fluid to bend in the direction of gravity,just as water flows down from a high place,which is consistent with reality.And it also shows the effectiveness of the proposed method in the fluid with gravity.In addition,we find a phenomenon where two small curved angles appear at the fluid inlets,as shown in Fig.20.It can be attributed that the artificial friction force is added to the Navier-Stokes equation to distinguish the solid and fluid regions,the impermeability in the solid region is a large but finite positive constant and the velocity in the solid region tend to zero but not equal to zero;the initial velocity at the inlet has a parabolic shape,and the velocities at the upper and lower parts of the inlet are close to zero,so they can flow through the small curved angles.

Figure 20:Optimal configurations for the horizontal channel with and without gravity

5.2.2 Splitter in 3D

This numerical example extends the proposed method with gravity to the 3D case. The design domain of the splitter is shown in Fig.21,which is discretized by 169,510 tetrahedron elements.The Reynolds number is set to Re = 10, and the fluid volume fraction is 10%. Gravity is a constant force, so the sign of gravity represents its direction. The gravity fbgis set to be (0,-10) and (0,10),respectively.

Figure 21:Design domain for the splitter problem

Refined mesh schemes with 1,592,844 tetrahedron elements are adopted to smooth the surfaces of optimal configuration,the smooth results corresponding to different directions are shown in Table 4.When the direction of gravity is downward,the fluid flows downward and the lower outlet is chosen,otherwise,the fluid flows upward and the upper outlet is chosen.It can be seen that when the direction of gravity is changed,the bending direction of the flow channel will change accordingly.It also further verifies the influence of gravity on the optimal configuration in a 3D case.

Table 4:Smooth optimal configurations corresponding to a different direction for the splitter problem

6 Conclusions

In this paper, we employ the parameterized level set method using basis functions to study the topology optimization problems of steady-state incompressible Navier-Stokes flow.The intermediate variable γ (φ)is defined to perform sensitivity analysis,which is used to drive the evolution of the level set function.Several numerical examples are shown to illustrate the validity of the proposed method applied in channel optimal design.

In Section 5.1, six numerical examples demonstrate the proposed method has flexibility and effectiveness in different kinds of design domains and mesh categories. The four-terminal device examples in Section 5.1.2 give similar observations as the research of others. That is the Reynolds number affects the shape of the optimal Configuration,and the larger the Reynolds number,the more inclined to obtain straight flow channels. In Section 5.2, two numerical examples of the horizontal channel in gravity illustrate the effectiveness of the proposed method in channel optimal design problems with gravity by comparing the optimal configuration with and without gravity.In addition,the proposed method is successfully applied to three-dimensional design problems and relatively smooth results can be obtained through the refinement method, so the proposed method has the potential in solving practical engineering problems.

However, the surface of the 3D results shown in this paper is still a little rough, which can be further improved by refining mesh and adopting other post-processing methods.Besides,the current work only considers constant body force(e.g.,gravity),and it is necessary to be extended to fluid flow problems with other kinds of physical body forces in the future,such as centrifugal force and Coriolis force.

Funding Statement:This research was supported by the National Natural Science Foundation of China (Grant No. 12072114), the National Key Research and Development Plan (Grant No.2020YFB1709401), and the Guangdong Provincial Key Laboratory of Modern Civil Engineering Technology(2021B1212040003).

Conflicts of Interest:The authors declare that they have no conflicts of interest to report regarding the present study.

Appendix A

The derivation of dΦ/dγ is shown in this appendix.The objective function Φ can be augmented as follows:

where U(γ) is the state variables of the governing equations holding the velocity u and the pressurep,Λ is the Lagrange multiplier and R(U(γ), γ)is the discretized Navier-Stokes equation in residual form.More detailed derivation of the form of residual can be found in[43].Using the chain rule,the total derivative d/dγ can be written as follows:

Note that,the sensitivity of the state variables with respect to the design variables,?U/?γ,is still an unknown quantity in the above expression. Therefore, we can introduce the adjoint Eq.(25) to eliminate ?U/?γ:

Therefore,the total derivative of the augmented objective function with respect to the intermediate variable can be written as follows:

Since the state solution at each loop will always make the residual vector R(U(γ), γ) to zero,the total derivative d/dγ become:

Therefore, the total derivative of the objective function with respect to the intermediate variable is:

It should be noted that the above procedures are automatically calculated by the commands in COMSOL Multiphysics(Version 3.5)called by the codes edited in MATLAB(Version 2008b).Then we only need to extract the calculation results of sensitivity information and use it to update the level set function.