Topology optimization of irregular flow domain by parametric level set method in unstructured mesh

In this study, we present a novel method for the topology optimization of the irregular flow domain using a parametric level set method (PLSM). Some improvement was applied on the CS-RBFs (radial basis functions with compact support)-based PLSM to make it suitable for nonuniform mesh, expanding the range field of engineering application of the PLSM. The optimization problem is solved by a gradient-based algorithm with Stokes equations as state constraints, and the objective is set to minimize the power dissipation subject to the volume constraint of flow channels. A PLSM is introduced to avoid the direct solving of the Hamilton–Jacobi partial differential equation, which can have the potential to break through the restriction of relying on structured meshes because no finite difference scheme is required. Then, a self-adaption support radius approach is presented to allow the parametric level set to be evolved on the nonuniformed mesh, which can expand the application of the PLSM to more complicated engineering problems with irregular geometric shapes. A volume integration scheme is applied during the design sensitivity analysis to calculate the shape derivatives, allowing the nucleation of new holes. Numerical examples in two and three dimensions are provided to demonstrate the effectiveness of the proposed method.


Introduction
Topology optimization method was commonly used to optimize the performances (i.e. stiffness, frequency, and so on) of a structure by rationally redistributing materials in the solid design domain. Since the seminal work in 2003 (Borrvall & Petersson, 2003), topology optimization of fluid domains has become an active research field. Because of its important practical significance in engineering, many new approaches have been presented. As for the fluid domains, the topology optimization problem is reshaped from "where should the material be placed?" to "where should the fluid be flow?" (Alexandersen & Andreasen, 2020), and the objective is usually set to reduce the power dissipation. Such problems are common in practical engineering cases. For example, in the design of engines and air conditioners, it always meets difficulties in design manifold to divide fluids with the minimized pressure drop.
In this paper, the finite element method (FEM)-parametric level set method (PLSM) [CS-RBFs (radial basis functions with compact support)-based] coupling method was proved to be adequate for Stokes problem. In the meantime, some improvement was applied on the CS-RBFs-based PLSM to make it suitable for nonuniform mesh, which has expanded the range field of engineering application of the PLSM. Meanwhile, the combination of PLSM and nonuniform mesh is discussed theoretically and the performance of the PLSM in fluid problems was analysed via a comparative study.
The implementation of the PLSM based on CS-RBFs for fluid topology optimization problems is described first and how to apply PLSM on unstructured grids is discussed then. Then, the formulation of the fluid topology optimization problem with the objective of minimizing power dissipation and the design sensitivity analysis using volume integration scheme is presented. After that, we show the numerical implementation with a flowchart. At last, we provide numerical examples using unstructured grids in both two and three dimensions to illustrate the validity and applicability of the proposed method.

Related Work
Topology optimization methods aiming at obtaining the minimum power dissipation of fluid flow were first presented in 2003 (Borrvall & Petersson, 2003), and then similar objective functions were used in many works (Gersborg-Hansen et al., 2005;Olesen et al., 2006;Wiker et al., 2007;Abdelwahed et al., 2009;Challis & Guest, 2009;Pingen et al., 2010;Deng et al., 2011;Koga et al., 2013;Dilgen et al., 2018). M. P. Bendsøe and O. Sigmund have pointed out that this objective function corresponds to maximizing the flow compliance (Bendsoe & Sigmund, 2004). This conclusion has revealed some relevance between minimizing power dissipation of Stokes flow and minimizing compliance of solid structures. After then, some topology optimization methods, including the bi-directional evolutionary structural optimization (Munk et al., 2017), phase-field (Garcke & Hecht, 2016), and level set method (LSM; Challis & Guest, 2009;Duan et al., 2008a, b, c;Zhou & Li, 2008), have also been used to solve such power dissipation problems. In particular, LSM is able to directly describe boundaries of structures and provide the pure 0-1 topology, making the fluid simulation more precise.
The LSM (Osher & Sethian, 1988;Sethian & Wiegmann, 2000;Wang et al., 2003;Allaire et al., 2004;Oliveira & Leonel, 2019;Oliveira et al., 2020) has become a popular approach for topology optimization. The LSM implicitly represents the boundary of a structure as the zero-level set of a higher dimensional function, which makes it easier to get an exact description of the structural boundaries. However, in the traditional LSM, boundaries are changed by solving the Hamilton-Jacobi partial differential equation (H-J PDE). This procedure may lead to some unfavorable numerical issues (Luo et al., 2008), and new holes cannot be conveniently created during the optimization process (van Dijk et al., 2013). Therefore, many alternative approaches have been developed to overcome the defects of the traditional LSM. For the fluid topology optimization problem, a variational LSM was presented by Duan et al. (2008a, b, c). In this modified LSM, a relatively smooth evolution can be maintained without any re-initialization scheme, and the drastic topology change can be handled easily (Duan et al., 2008a). After that, some studies based on a similar level-set model have also been presented to expand its applications in the field of fluid topology optimization (Zhou & Li, 2008;Challis & Guest, 2009). However, the variational formulation introduced in the LSM may cause some other computational complexity. Concretely speaking, a functional integral was built as a metric to characterize how close the function is to a signed distance function, which involved a new term into the H-J PDE (Duan et al., 2020).
Another way to avoid solving the H-J PDE directly is to use the PLSMs, which approximates the level-set function by the interpolation schemes. The early attempt to parametrize the levelset function (Wang & Wang, 2006a, b) includes the interpolation scheme via using the multiquadric and inverse multiquadric splines (Belytschko et al., 2003;Norato et al., 2004). Aimed at the numerical difficulties of the above method, a full parametrization LSM using CS-RBFs was published by Luo et al. (2008), which was named the PLSM. PLSM has been applied to different topology optimization problems (Li et al., 2016(Li et al., , 2018Wang et al., 2014Wang et al., , 2016. To solve fluid topology optimization problems, Pingen et al. (2010) formulated a parametric level-set approach using the Lattice Boltzmann method (LBM). Although there are no distinct advantages or disadvantages of the FEM over the LBM scheme in solving pure fluid problems (Pingen et al., 2010), the FEM has been integrated inside many commercial and open-source simulation software, such as ANSYS, Abaqus, Comsol, FEniCS, etc. Using FEM makes our method easier to be integrated inside those software and thus gives our method more practical value. At the same time, the universality of the FEM scheme makes it suitable to solve multiphysics problems.
In many real-world engineering problems, a design domain with complicated geometry is often encountered. The uniformly structured grids are usually used to discretize regular computing domains, and thus they have limitations in solving complex engineering problems. For instance, most of the studies based on LSM only considered uniformly structured grids. Since the H-J PDE is no longer needed to be solved directly as is mentioned above, an explicit finite difference scheme can thus be avoided, which makes the structured grid needed by the finitedifference calculation not necessary. Because of that, there is no fatal hindrance to applying PLSM on unstructured grids other than some possible numerical problems caused by CS-RBFs (Luo et al., 2009). By applying a self-adaption support radius approach to the PLSM algorithm, the numerical problems can be avoided, leading to no significant differences between the final optimal design using the structured and unstructured grids. In terms of the problem that new holes cannot be nucleated in the optimization process, a hole nucleation mechanism (Burger et al., 2004;He et al., 2007;Abdelwahed et al., 2009) was introduced into the LSM. However, in PLSM, as the velocity field has been extended naturally over the design domain, the hole generation problem may be caused by the integration scheme. Luo et al. (2009) have approved that a volume integration scheme can be used to create new holes; thus, the additional nucleation mechanism is no longer required, which gives the PLSM more engineering generality. However, the application of the PLSM was mainly on structural related problems, so the applicability on multiphysics problems such as the Stokes problem was rarely discussed. Meanwhile, the method was always discussed only in uniform structured mesh, which has limited its potential in solving practical engineering problems.

Level set parametrization scheme
We introduced the PLSM at the irregular fluid domain. The traditional LSM is using an implicit level-set model (H-J PDE) to describe structural boundary as the zero-level set [ (x, t) = 0] of a higher dimensional scalar function (x, t): where D is the design domain and is the admissible shape ( ⊂ D). The parameters x and t, respectively, represent the coordinate and the pseudo-time. In the design domain D ⊂ ℝ d (d = 2 Downloaded from https://academic.oup.com/jcde/article/9/1/100/6490315 by guest on 06 January 2022 or 3), the material filled or not is decided by positive or negative of the scalar value; then, we can get a clear "0-1" topology design naturally. For the topology optimization problem of the flow domain, the structural boundary represents the fluid-solid interface. By introducing the pseudo-time t, the motion of the fluid-solid interface is equivalent to solving the first-order H-J PDE: where v n is the normal velocity field (Wang et al., 2003;Allaire et al., 2004). The PLSM can be regarded as a discretion of the classical HJlevel-set function by RBFs; in other words, a radial based function (RBF) interpolation of the HJ-level-set function was implemented. The interpolation with RBFs is a well-established technique in the area of multivariate approximations (Wendland, 1995(Wendland, , 2006. The level-set function was interpolated as follows: where ϕ = [ϕ 1 , ϕ 2 , ϕ 3 , ..., ϕ N ] are the RBFs and s = [s 1 , s 2 , s 3 , ..., s N ] is the vector of the weights or expansion coefficients of the RBFs. N is the number of RBF knots. By substituting equation (3) into equation (2), we can transform the H-J PDE into a system of first-order ODEs (Wang et al., 2003): where RBFs introduced into shape and topology optimization problems can be approximately classified into two types, one is CS-RBFs and the other is GS-RBFs (RBFs with global support). To implement the interpolation leads to a collection matrix, whose size is the square of the number of RBFs' center knots. The GS-RBFs always lead to a full collection matrix (Belytschko et al., 2003); however, when the number of the center knots becomes large, the matrix will take too much memory; thus, the sparsity of collection matrices becomes important. Because of that, in this study, we choose the CS-RBFs with 2k continuity presented by Wendland (1995), which may be more preferable due to their simple forms with minimal degree (Luo et al., 2008). We adopt the CS-RBF with C 2 smoothness: whose gradients are where r defined in Euclidean space can be expressed by where x I is the coordinate of knot I, d ml is the support size that can be calculated by using d mI = d max · C I , where d max decides the basis function's influence range at knot I, and C I is used to guarantee that the neighborhood knots can be included in the influence range.

The self-adaption support radius approach
To facilitate the calculation and improve the accuracy, the knots are usually set to coincide with the finite element node when using the FEM, as in Fig. 1a. In the uniform grid, the C I is a fixed value. However, in the unstructured grid, the size of each element may have a significant difference; for example, in Fig. 1b, a remarkable nonuniform mesh was formulated. In this figure, it can be noticed that with a fixed value of support radius, when the element size gets small, there are too many knots that are included, which will increase the computational cost. On the opposite, few knots included inside the influence ranges because of large element sizes may lead to meet the nonsingularity of the interpolation (Luo et al., 2009). To summarize, the support radius is an important concept in CS-RBFs. To avoid meeting the above problems, the C I should be set differently in each knot. We presented a self-adaption support radius approach based on the average of the distances from knot I to all its neighborhood knots. For example, the C I of the CS-RBF on knot N 0 shown in Fig. 2a should be .
(9)  This approach makes support radius of CS-RBFs on each knot changes along with the sizes of its surrounding elements, guarantees that a similar number of knots can be included in each influence range, as is shown in Fig. 2b.
However, this approach also has certain limitations, one of which is that the low element Jacobi will lead to an inaccurate value of C I . However, in the meantime, the low mesh Jacobi will also make the finite element calculation inaccurate; such elements will thus be avoided during the mesh generation.

The "ersatz material" model
When an element is crossing and divided by the zero-level set boundary, the material parameters of those elements are relevant to the evaluation of the design sensitivities. The material parameters represent some material proprieties to model the "filled" or "void" status of one region. For example, in the stiffness optimization problem, the material parameters represent Young's module, and in the fluid optimization problem they represent the permeability. In the numerical implementation, we chose the "ersatz material" model (Allaire et al., 2004). The "ersatz material" model amounts to fill the elements crossing the boundary by a weak phase, and the elements' material parameters are handled under the assumption that it is proportional to the area fraction of the filling material. Because that the dimensions of one element are small comparing to the whole grid, we can locally assume the level-set function to be flat. It means that in an FEM element, the level-set function is approximated to be linearly varying. The material parameters of such elements are thus calculated by linear interpolation, as is shown in Fig. 3a. In Fig. 3b, we have > 0 in the red area, the material parameters of the elements inside the red area are one and outside are zero, and the material parameters of the elements crossing the boundary are between one and zero as so-called "ersatz material."

Formulation of the optimization problem
The first 2D fluid topology optimization approach (Borrvall & Petersson, 2003) was based on the parallel-plates model; for threedimensional (3D) problems, the model of the parallel plates loses its physical meaning, yet the Brinkman-type model that consid-Downloaded from https://academic.oup.com/jcde/article/9/1/100/6490315 by guest on 06 January 2022 ers the material parameter as frictional resistance has no such problem. The Brinkman equations can be expressed as a frictional resistance term added to a standard-steady-Stokes equations, which is given by where u and p are the velocity and pressure fields, respectively, μ is the fluid kinematic viscosity, and κ is the porosity that provides the frictional resistance. Combined with the PLSM approach above, the porosity coefficient κ can be modeled as a linear function of the material parameter of the "ersatz material": where κ ∈ [κ min , κ max ] and H( ) is the material parameter in the "ersatz material" model related to the level-set function. For the elements inside of the "unfilled" area, the κ becomes very large, which forces near-zero flow, consistent with the solid region, and the "filled" elements corresponding to low κ indicate the fluid region accordingly. The elements crossing the boundary can be regarded as porous medium regions. The cost function of the problem has two major forms, one is the linear combination of nodal field variables and the other is the dissipated energy inside the fluid flow (Gersborg-Hansen, 2003). In this study, we chose to use the dissipated energy that determines the most energy-efficient design. The cost function was defined as half the energy loss due to viscosity, which is rewritten as where u is the velocity solution of the Brinkman equations (10) and the two integrated terms are, respectively, associated with the viscous force and the frictional resistance. Using standard FEM discrete matrix form, equation (10) can be obtained as references (Zienkiewicz et al., 2013): and N is the number of the elements in the design domain, G is the divergent operator matrix for the continuity equation, and G T is the gradient operator matrix. K is the velocity stiffness matrix consisting of symmetrical matrices K d and K κ that come from the diffusive term and the frictional resistance term separately. The objective functions described in equation (12) are rewritten in a discrete form: where U denotes the solution vector of velocity coefficients for each node, i.e. U = [u 1 , u 2 , u 3 , ..., u n ] T , and n is the number of nodes. Finally, the topology optimization problem in discrete form can be established by where J is the objective function, V max is the volume limitation of the fluid region. s max and s min are the upper and lower bounds of the design variables, respectively.

Sensitivity analysis
Here, the shape derivative (Sokolowski & Zolesio, 1992;Choi & Kim, 2005) is introduced to compute the design sensitivity of the fluid regions represented by the level sets (Wang et al., 2003;Allaire et al., 2004). The weak form of Stokes equations can be expressed as Downloaded from https://academic.oup.com/jcde/article/9/1/100/6490315 by guest on 06 January 2022  where g represents the boundary loads of the fluid domain and q and v are the test functions of the weak form (Borrvall & Petersson, 2003). Differentiating the first equation in the equations above with respect to t yields Since boundary loads of the fluid domain do not get any influence from the evolution of the level-set function, the shape derivatives of both sides of equation (18) can be respectively expressed as follows: where δ( ) is the derivative of the Heaviside function H( ) anḋ x = dx/dt. Considering the terms containingv in equation (19), we can establish the following equation (Li et al., 2018): Substituting equation (20) into equations (19) and (18) yields Downloaded from https://academic.oup.com/jcde/article/9/1/100/6490315 by guest on 06 January 2022 Since the velocity stiffness matrix is symmetrical, equation (21) can be rewritten as Meanwhile, according to the incompressibility of the fluid, we have Then, equation (22) can be reduced as The shape derivatives of the objective function J can be expressed as By substituting equation (24) into equation (25), we have According to the chain rule, the derivative of objective function J with respect to t can be calculated as By comparing the equations (26) and (27), we can find that the derivative of objective function J concerning to the design variables s becomes Similarly, the derivative of the volume constraint in equation (16) with respect to the design variables can be obtained as In equations (28) and (29), a volume integration scheme is used to calculate the design sensitivity. The advantage of it is that new holes can be added since the design sensitivity is evaluated over the entire design domain and the velocity field is naturally extended to the whole domain (Luo et al., 2009). More specifically, the traditional LSM only calculates the sensitivities on the boundary of structure, which causes that the velocity field only exists on the boundary and other regions' evolution is nearly stopped. However, in the present method, the evolution velocity exists at every point in the design domain during the whole iteration process just like the variable density method does, which means that the holes can naturally be created during the iterations.

Optimization Algorithm Implementation
In this section, we introduce the updating scheme of design variables and then show a flowchart for the numerical implementation. After obtaining the sensitivities of the objective function, an OC-based algorithm is used to update the design variables (Sigmund, 2001). Based on the optimization problem defined in equation (16), one can form a heuristic scheme as Downloaded from https://academic.oup.com/jcde/article/9/1/100/6490315 by guest on 06 January 2022  where In order to stabilize the iterative process, some artificial parameters, the damping factor ζ (0 < ζ < 1), the move limitm = 0.01(0 < m < 1), and a small positive number μ, are introduced (Li et al., 2018). The Lagrangian multiplier can be updated by a bi-sectioning algorithm during the numerical calculation.
The flowchart for the numerical implementation is shown in Fig. 4. Here, values of the level-set function at each knot can be calculated by multiplying the expansion coefficients s with a constant interpolation matrix that needs to be calculated only once before starting the optimization loop. At each step of the optimization iteration, the finite element analysis is solved first and the derivatives are then obtained using the solutions. The value of the objective functional J and its change |J (k) − J (k−1) | are also calculated in every step to determine whether the iterative process converges. Generally speaking, the present method contains the following steps: Step 1: Calculate C I at each knot (FE node) and set up the interpolation matrix and the initial level-set function; Step 2: Solve the Stokes problem with FEM and get the velocity field u; Step 3: Calculate the derivatives of the objective functional with u; Figure 11: Design domain and boundary conditions of 2D manifold.    Step 4: Update the level-set function using OC and check the convergence. If not converged, go to step 2.

Numerical Examples
In the numerical implementation, we solve systems of coupled PDEs by the FEM. As mentioned before, we used the "ersatz material" model to handle the material parameters that represented the frictional resistance for fluid problems. To achieve that, the coefficients of the frictional resistance term were set differently in each element, more specifically, at each quadrature point. For the design sensitivity analysis, we respectively calculated integrations over each element of the design region.
In this section, we first introduce a 2D classical numerical case to prove the validity of our method and then change the element type of the grid and initial structure to discuss the influence of different factors. After that, we respectively give the results of 2D and 3D cases with irregular design domains to imitate divided manifold design problems using unstructured grids. In each example, the upper and lower bounds of the porosity coefficient were, respectively, set to be 1e4 and 1e-4.

2D examples
For the first example, a rectangular design domain was considered. The inflow and outflow pipes were placed, respectively, in the up left and bottom right corners of the design domain. The inflow boundary condition was set on the left boundary of the inflow pipes with a constant inflow velocity u 0 = 0.1 along the boundary's normal. The outflow boundary condition on the bottom was implemented as p = 0. On other boundaries, no-slip boundary conditions were implemented, as is shown in Fig. 5.
The fluid kinematic viscosity μ was set to 1 and the volume fraction was set to 0.3.
To prove the mesh adaptability of our method, we prepared two different mesh generation forms, one is pure rectangular gird and the other is a pure triangular gird. We used HyperMesh for mesh generation and the element size of both meshes was set to be 0.2, as shown in Fig. 6.
Figures 7 and 8 display the topologies at different stages using two different meshes. The topology convergence processes and the final design of the rectangular and triangular girds were similar to each other. Also, in the triangular girds, the convergence seems to be faster. The reason is that under similar element sizes, the influence range of one CS-RBF can include more knots in the triangular girds, which has strengthened the contribution of a local evolution of level-set function to that of the entire domain. In Fig. 8e, we can see the initial sign distance functions that are common in LSM-based topology optimization. In these two examples, the results only contain a single path, which is similar to the results of shape optimization. However, in Fig. 8e-h, the LSM shows a strong capability of changing topologies, which gives the method better freedom of design than shape optimization. Besides, in Fig. 8e-h we can see that the evolutionary velocity field exists everywhere in the design domain, which guarantees the capability of our method to create new topologies.
The cases mentioned above are started with a lot of initial holes inside the design domain, which makes it difficult to show the ability to create new holes of our method because the holes are simply merged and faded away. So, we choose an initial design with no holes. The topology convergence process was shown in Fig. 9a-d; as we can see, although the final design has no holes, two holes temporarily appeared during the Downloaded from https://academic.oup.com/jcde/article/9/1/100/6490315 by guest on 06 January 2022  optimization process in Fig. 9b and c. This has confirmed that our method has the ability to create new holes. Besides, the final design is similar to that using an initial design with holes inside. Besides, Fig. 9e-h shows that the size and positions of the initial holes have no influence on the final design.
Finally, we give the result with pressure and velocity fields. The optimization histories of objective function value and volume fraction are shown in Fig. 10. The maximum number of iteration steps was artificially set to 100. However, the optimization iteration had nearly converged at the 30th step, and the design topology had been established at that step. After that, the optimization process was acting more like a shape optimization. In the result's plot, the pressure and velocity fields and the final design accord with that of the classical cases in references (Borrvall & Petersson, 2003;Gersborg-Hansen, 2003).
To prove the adaptability of irregular flow domains, we provided a 2D manifold example. The manifold has one inflow pipe joint on the left and four outflow pipe joints on the bottom. The shape of the design domain was chosen to be an irregular polygon to simulate the interference by other parts when assembly, as is shown in Fig. 11. The inflow velocity u 0 was set to 0.2 and the inflow velocity u 1 was set to 0.1; the volume fraction was set to 0.75. The grid was generated by the Delaunay triangulation algorithm that is shown in Fig. 12.
The initial design was started from sign distance functions mentioned before. In the first 100 steps, the change of topology was obvious, but after step 100 the iterative process was more like shape optimization, as is shown in Fig. 13. In the result's plot in Fig. 14, the iterative process was stable and no numerical oscillation was produced.
Finally, in order to test the performance of our method in irregular meshes, the manifold case with irregular mesh is given. The mesh in the design domain was generated by gradient-sized elements with mesh refinement on the bounding layer, representing common fluid meshes, as is shown in Fig. 15a. The result is shown in Fig. 15b, which is similar to that shown in Figs 13 and 14, confirming the robustness of the method in problems with irregular meshes.

3D example
The 3D example is a divided manifold design problem. It can represent most of complex engineering problems in which the unstructured mesh is needed. In the case, we set five circular  pipe joints at both bases of the cylindrical design domain, one for the liquid to flow in and four in the opposite base for flow out. The normal outflow velocity was set to 0.1, and the normal outflow velocity was set to make the inflow flux equal to the outflow flux to satisfy the incompressibility. The volume fraction was set to 0.3. Figure 16 shows the design domain and boundary conditions of the problem. Figure 17 shows the tetrahedral mesh generation result of the above geometry. The hexagonal mesh was also tested in this example, yet either the final design or the convergence shows no obvious difference compared to those using the tetrahedral mesh. Because of that, only the result using the tetrahedral mesh is shown in this paper.
As the iterative process and final topology shown in Figs 18 and 19, identical channels connecting with each outflow pipe joint were created after the optimization process owing to the domain's symmetry, which was consistent with engineering experience. The fluid channel structure has a trumpet-like outline that can reduce resistance during the liquid diversion, and the shape of which was close to the classical "diffuser" case presented in 2D Stokes flow topology optimization (Borrvall & Petersson, 2003). Besides, the circular holes of the initial design in the 2D design domain were replaced by spherical holes in the 3D domain.

Computational time consuming
With regard to the computational time consuming, the CPU times of both 2D and 3D cases (in Figs 10 and 19) were tested. The problems are solved on AMD Ryzen5 3500 × 6C6T @4.1GHz. During the process of FEM, MUMPS was used as the sparse di-rect solver. As a control, a density-based method with sensitivity filtering (Sigmund, 2001) and density projection method (Wang et al., 2011) using the same mesh and boundaries is also tested. The results are shown in Table 1.
In the table, we can see that compared to the traditional density method, the PLSM takes more time in a single step. However, the density method takes more steps to reduce the number of the elements with intermediate density and achieve a 0-1 topology; thus, by using the PLSM, the total time was saved.

Conclusions
We proposed a new methodology of performing topology optimization for fluid flow by a PLSM using the FEM. Comparing to the traditional LSM, PLSM has overcome the numerical difficulties in solving H-J PDEs. Based on the advantages of PLSM that can avoid solving H-J PDEs in the finite difference scheme, the method was applied to unstructured grids. A self-adaption support radius approach was also presented to reduce precision loss relevant to the feature of CS-RBFs. During the optimization process, the PLSM can directly get clear structure boundaries, which profits computational accuracy for CFD and simplifies postprocessing. For the design sensitivity analysis, a volume integration scheme has been applied to give our algorithm the ability to create new holes in the design domain. Benefit from this, the additional nucleation mechanism is not required, which has reduced the computation cost. Compatibility with unstructured mesh has made our method easier to be integrated inside common simulation software, and using FEM also gives it adaptation to multiphysics problems. The numerical examples show that the proposed method can be applied to 2D and 3D irregular Downloaded from https://academic.oup.com/jcde/article/9/1/100/6490315 by guest on 06 January 2022 domains and create a smooth shape of the flow path. However, the present method is still limited in solving nonlinear problems. To overcome that, the sensitivity analysis format needs to be further improved. Meanwhile, the elements cut by the level set boundary are actually anisotropic, but in the "ersatz material" method they are modeled to be isotropic, causing some precision losses. Using the adaptive mesh generation method to refine mesh near the level set boundary might help. Because of the maturity of automatic generation methods of finite element meshes, unstructured gird has been widely applied to the simulation and optimization of mechanical parts with irregular geometric shapes, which has given the present approach great application value, and it can be easily applied in the engineering area and extended to solve other fluid-related design problems.