Multi-Patch Black-White Topology Optimization in Isogeometric Analysis

Qingyuan Hu,Yuan Liang,Menghao Liu,Manfeng Hu and Yawen Mao

1School of Science,Jiangnan University,Wuxi,214122,China

2Department of Engineering Mechanics,Dalian University of Technology,Dalian,116024,China

ABSTRACT Topological optimization plays a guiding role in the conceptual design process.This paper conducts research on structural topology optimization algorithm within the framework of isogeometric analysis.For multi-component structures,the Nitsche’s method is used to glue different meshes to perform isogeometric multi-patch analysis.The discrete variable topology optimization algorithm based on integer programming is adopted in order to obtain clear boundaries for topology optimization.The sensitivity filtering method based on the Helmholtz equation is employed for averaging of curved elements’sensitivities.In addition,a simple averaging method along coupling interfaces is proposed in order to ensure the material distribution across coupling areas is reasonably smooth.Finally, the performance of the algorithm is demonstrated by numerical examples, and the effectiveness of the algorithm is verified by comparing it with the results obtained by single-patch and ABAQUS cases.

KEYWORDS Isogeometric;topology optimization;nitsche;multi-patch;integer programming

1 Introduction

Topology optimization of structures aims to find reasonable material distribution under certain objectives.Combined with Finite Element Method (FEM), by utilizing powerful computational resources,nowadays topology optimization is widely applied for engineering problems.As a new generation of FEM,Isogeometric Analysis(IGA)[1]adopts Non-Uniform Rational B-Spline(NURBS)functions as shape functions.Owing to some conceptual differences between IGA and classical FEM,since then researchers have been studying topology optimization methods in the framework of IGA in many angles.

In 2008,Wall et al.developed a shape optimization strategy by moving control points[2].In 2010,Youn and his collaborators proposed a topology optimization method using trimming techniques[3,4].In 2017,Wang et al.focused on lattice materials optimization using asymptotic homogenization schemes[5],then in 2020 Wang and his collaborators proposed a high-efficiency topology optimization method by using multilevel mesh, multi-grid conjugate gradient method and local-update strategy[6].By transplanting and extending the Moving Morphable Components (MMC) [7], and Moving Morphable Voids(MMV)methods[8],corresponding IGA-based methods were proposed by Hou et al.[9],Xie et al.[10],and Gai et al.[11].In 2019,Gao et al.conducted the research about isogeometric topology optimization for continuum structures by a constructed density distribution function [12].For a recent code implementation of topology optimization for structures using IGA in MATLAB,interested readers are referred to[13].

Since in engineer practice a structure model is usually made upon NURBS patches of parts,in this research we mainly study topology optimization of multi-patch IGA, especially in 2D.In FEM framework, Jeong et al.used a condensed mortar method to glue dissimilar interfaces for topology optimization and achieved good performance[14].In mortar method,the newly introduced corresponding Lagrange multipliers defined along the new interfaces act like forces to enforce the interfaces of both sides to stay tight.The mortar method adds new degrees of freedom(DoFs)into the system,however the newly introduced DoFs can be condensed as in[15].The mortar method is widely adopted to glue patches in IGA,such as in[15,16].Usually such a domain decomposition technique can also be used for contact problems,see[17,18],thus the mortar method can also be employed for contact topology optimization[19].

In this paper, we adopt the Nitsche’s method for NURBS patch coupling, by adding external stiffness by multiplying the averaging force and displacement gap along coupling interfaces.The Nitsche’s method introduces no additional DoFs,thus no special treatment is needed for the system,which makes the method suitable for contact problems[20].Moreover the Nitsche’s method is tested out to be a steady numerical method for patches coupling[21].In regard to the choice of the topology optimization method,we adopt the sequential integer programming method and Canonical relaxation algorithm proposed by Liang et al.[22].This method utilizes the sensitivity analysis to construct separable integer programming sub-problems.Besides,Liang et al.proposed the Canonical relaxation algorithm to efficiently solve the integer programming sub-problems.Thanks to the discrete variables,the post-processing treatment to eliminate the gray densities can be avoided.In addition,this method is especially suitable for multi-objective optimization.Combing with the IGA-based methods, the optimized results can be easily parameterized.For further description and code implementation,please refer to[23].

The structure of this paper is organized as follows.In Section 2, the isogeometric multi-patch analysis based on the Nitsche’s method is introduced.In Section 3,we outline the calculation process of the sequential integer programming and discuss the methods for sensitivity filtering and multi-patch treatment.In Section 4, several numerical examples are carried out to show the performance of the proposed method.Conclusions are given in Section 5.

2 Isogeometric Multi-Patch Analysis

2.1 Brief Introduction of IGA

FEM analysis needs two sets of basis functions, one is used for geometry model building, the other one is used for displacement field construction.The main idea of IGA consists in using NURBS as basis functions for both geometry models and the displacement field.The present research mainly focuses on 2D structures,thus only 2D surfaces are introduced,however 3D volumes can be naturally constructed by adding one more dimension.

A 2D surface, denoted by S(ξ,η), can be built by utilization of the NURBS basis functions as Eq.(1)

whereNi(ξ,η)are NURBS basis functions, generated by given sets of knot vectors together with corresponding weights,and Pi(x,y)are called control points.

Following the IGA concept,we approximate the displacement field by Eq.(2)

where Ui=(ui,vi)Tare nodal displacement degrees of freedom attached to each control point.Note here we use the same set of NURBS basis functionsNias Eq.(1).

As illustrated in Fig.1,different from FEM,in IGA the mesh and elements are generated naturally by the given knot vectors, and the elements are usually curved.In addition, in a traditional FEM element,nodes are located on the edges or corners of the element,while in IGA some control points are located outside of the element.For constructing IGA-suitable parameterization for complex CAD models,it is recommended to read[24].

Figure 1:Illustration of the concepts of elements and nodes in FEM(a)and IGA(b)

2.2 The Nitshce’s Method for Multi-Patch Analysis

Consider an interface problem,in which the domainΩis decomposed into two sub-domainsΩm,where the subscriptm=1,2 is used to mark the partitioned domain and corresponding variables.The shared boundary betweenΩ1andΩ2is denoted byΓ, andnmis the unit normal along the interfaceΓ, pointing out ofΩm.Denote the displacement field to be computed as u:Ω→Rd(d≥1), and denote by v ∈V an arbitrary test function that can represent the virtual displacement field,where V is a functional space of admissible fields.Then starting from the original FEM and IGA weak form in Eq.(3)

by further defining the jump and average operators along the interfaceΓin Eq.(4)

we finally reach at the weak form of the Nitsche’s method for multi-patch analysis as Eq.(5)find u ∈V:

where theθis named as the Nitsche parameter, and theγis called the stabilization parameter.The newly added terms are identified as Nitsche’s contribution,in which the second and the third terms are symmetric and they play a key role to glue the two patches,and the forth term is a stabilization term.

The Nitsche parameterθallows us to select some variants of the Nitsche formulation, that yield different theoretical properties and different degrees of dependency w.r.t.the stabilization parameterγ:

· forθ= 1,the standard symmetric Nitsche’s method is obtained,and a suitable choice forγis necessary in order to recover well-posedness and optimal accuracy;

· forθ=0,some terms are canceled and we obtain a simple formulation;

· forθ=-1,the skew-symmetric Nitsche method is obtained,and it is highlighted that for linear boundary/interface conditions,stability and optimal convergence are ensured even forγ= 0,resulting in a parameter-free formulation.

In this research,we setθ=1 to use the standard Nitsche’s method in order to keep the modified system symmetric.For the stabilization parameterγ,we adopt the experienced formula in[25].For more details about the derivation and theory of the Nitsche’s method,please refer to[20].

From the matrix point of view,imagine we have the original discretized IGA formulation in matrix form,as in Eq.(6)

in which K is the global stiffness matrix,U is the vector of nodal displacement DoFs,and F stands for the nodal force vector.Then assuming we have two NURBS patches to glue,the above formulation can be decomposed as Eq.(7)

where the subscript stands for the patch index.By pre-defining or other identification techniques,the interface elements are now marked by superscriptsandmfor interface elements along the two patches respectively.After collecting the indexes of the interface elements,we can rewrite the discretized formulation as Eq.(8)

the tilde means the remaining part,for example=K1-Ks1.

As illustrated in Fig.2, the Nitsche’s method couples the two patches by adding stiffness to the interface elements as Eq.(9)

in which the Nitsche related terms can be derived from Eq.(5)as shown in Eq.(10)

Figure 2:Illustration of the nitsche’s method

In cases of 2D linear elasticity,we have Eq.(11)

as displacement matrix,and Eq.(12)

as strain-displacement matrix.In addition,D1and D2are material constitutive matrix,however in this study we adopt same materials for both sides.Note that the integrations are performed along the slave elementsΓs,and n1is the direction matrix for the strain-displacement adaptation,as in Eq.(13)

wherenxandnyare components of vectorn1along x-and y-directions respectively at each quadrature point of the slave interface.

3 Multi-patch IGA Topology Optimization

3.1 Algorithm Description of Sequential Integer Programming

Consider a minimization problem of structural compliancec(p)under material quantity constraint,in the context of pixel(element-wise)description,the optimization formulation is give by Eq.(14)

where the structure is discretized intoNelements,each element area is represented byvi.The discrete design variableρstands for the element density by taking a value of 0 or 1.

In the algorithm proposed by Liang et al.[22,23], the above problem is first approximated by a series of separable integer programming sub-problems as Eq.(15)show

where the compliance is defined in Eq.(16)as

and the sensitivity of design variable is defined in Eq.(17)as

It should be noted that the density design variableρstill remains discrete.The relationship between the Young’s modulus and the design variables is assumed to beE=E0ρp,wherepis the penalty factor of design variables and it is usually taken asp=3 in common practices.In addition,in order to avoid the singularity of the stiffness matrix due to material reduction,ρmin=10-3is taken.

After filtering the sensitivities of the design variables frombto ˜busing the method described in Subsection 3.2, the series of separable integer programming sub-problems Eq.(15) are solved using the regularized relaxation algorithm.Firstly, the constraints for the discrete design variablesρiare relaxed into the following equality constraints as Eq.(18)

For Eqs.(15)and(18),define the Lagrangian function as Eq.(19)

where the dual variableλis introduced for the material quantity constraints,and the dual variableσis introduced for the equality constraints.Since the Lagrangian function above is differentiable,its first order critical condition in terms of the density variables is given in Eq.(20)

Therefore the original design variableρiand the two dual design variablesλi,σisatisfy the formula

which can be simply referred as the original duality relationship.Substitute Eq.(21)into Eq.(19),we get the dual problem of explicit differentiability with any order as the following Eq.(22)

In order to ensure the above dual problem has critical dual variables that are strictly larger than zero,and thereby to ensure the positive definiteness of the above optimization problem,the objective function of Eq.(22)is perturbed to

whereβ=2×104is the perturbation parameter.

The presented algorithm is composed by Eqs.(21),(24)and(25).Specifically,the iterative solving process is as follows.In the first step,we initialize the dual variableλ0,and solve for the dual variableσ0according to the formula Eq.(25).In the next step,the design variableρi0is obtained by using the Eq.(21).If a certain convergence criterion is satisfied,it is assumed that approximate integer solutions of the integer programming sub-problems are obtained.Otherwise,we takeσ0into Eq.(24)to get a new dual variableλ1,and the iteration continues.

Due to the integer programming sub-problems are constructed by sensitivity information, the accuracy of these sub-problems is only valid within a local range of current design variables.In order to ensure the approximation accuracy of the integer programming sub-problems, and to make the solutions of the sub-problems converge slowly to the solutions of the original problem,a moving limit strategy should be adopted.Specifically,we gradually reduce the material consumption constraint by a volume fraction reduction factorχ∈[0.95,1).Only when the relative change of the objective function in two adjacent iterations is less than a given value (for example 10-3), the volume fraction can be reduced.Otherwise,it is necessary to continue to iterate at the current volume fraction until the target volume fraction is reached.This strategy is expressed as Eq.(26)

whererkdenotes the number of volume ratio reductions.In this paper,χ=0.98 is taken,which means only 2%of the current volume change is allowed for each iteration loop.

For the present topology optimization technique by the sequential integer programming method,by setting the element density as 1 or 0 in the parameter space,we could achieve a black-white topology design accordingly in the physical space with clear and curved boundaries,as shown in Fig.3.

Figure 3:Black-white topology optimization in IGA

3.2 Curved Elements Sensitivity filtering by the Helmholtz Equation

In IGA, although the meshes in the parameter space are structured, however the meshes in the physical space are usually curved after mapping.The classical sensitivity filtering technique needs to search adjacent elements and then perform averaging methods,this is not suitable for complex curved meshes in IGA.In this study,the implicit filtering method based on Helmholtz equation[26]is used,which not only saves the amount of calculation effort for searching of adjacent elements, but also ensures the filtering effect.The implicit filtering equation is written as Eq.(27)

where bEis the element-wise sensitivity array before filtering,is the node-wise sensitivity array after filtering,and in Eq.(28)

we have the isogeometric stiffness matrix for the scalar problem.Reand ∇Redenote the basis function interpolation matrix and the basis function gradient matrix respectively, symbol ∑stands for the element stiffness matrix assembly operator.The filter radiusrmincontrols the range of the filtering circle by Eq.(29)

The transformation matrix is given as Eq.(30)

and it is used for the conversion between element-wise array and node-wise array, for instance in Eq.(31)

Finally,the filtering formula is obtained as Eq.(32)

Note that KFis only related to the original mesh since its generation and remains unchanged in the following topology optimization process,thuscan be calculated after the mesh initialization,then be utilized repeatedly to save computation resources.

3.3 Multi-Patch Treatment

The flow chart of the discrete variable topology optimization algorithm is drawn in Fig.4.After the process of IGA analysis with the Nitsche’s method, we obtain a single stiffness matrix K which contains the stiffness contribution of all patches as well as the Nitsche’s coupling.By solving the matrix equation Eq.(6)we reach at U in which the displacements of all the control points of all the patches are listed.According to the sequential integer programming algorithm,the calculation of the compliancecein Eq.(16) and the following variables are defined and performed element-wisely, thus all the optimization process can be done without considering the multi-patch issue.

Figure 4:Flow chart of the discrete variable topology optimization algorithm

However,for multi-patch sensitivities,since the mesh resolutions of different patches are usually different,the magnitudes of the sensitivity values correspond to meshes of different densities are also different.If these sensitivities are directly substituted into the topology optimization algorithm, it could lead to inappropriate optimization results,especially along the coupling interfaces.In this paper,a simple sensitivity filtering method for the coupling interface is proposed.Suppose that the numbers of elements along both sides of the coupling interface isM1andM2,and the corresponding sensitivities areb1andb2,respectively.Refine the meshes on both sides to the same densityM,that is,to find the least common multiple ofM1andM2as Eq.(33)

Correspondingly,the sensitivities of both sides is scaled as Eq.(34)

Since the meshes on both sides are in the same level at this stage, sensitivity filtering can be performed by using a simple averaging method,and the filtered sensitivitiesandcan be calculated by,like Eq.(35)for instance

Finally, the sensitivity is restored to the original meshes by adding all the values within each element

Note that although the relationship between the sensitivities and the mesh resolutions is not linear,we still scale the sensitivities linearly according to mesh refinement by Eq.(34),which does not harm much because in the last step(Eq.(36))an opposite treatment is performed.

4 Numerical Examples

In the following examples,the material of Young’s modulusE=106Mpa as well as Poisson’s ratioν= 0.3 are adopted for IGA analysis,and the concentrated force isF= 1000 N.For the topology optimization process,it is requited that current volume fractioncannot be further reduced until the relative change of the adjacent two objective functions reaches 10-3or even less.In addition,the sensitivity filter radiusrminis set to be 4 times of the minimum cell length.

4.1 Double-Beam Structure

In order to show the performance of the method in terms of coupling and optimization,and to compare with the existing result of the Mortar method [14], we consider a double-beam structure with overlapped interface as shown in Fig.5.The beam size is 20 mm × 2 mm, and the two beams are partially overlapped,resulting a shared interface of length 15 mm.The upper beam is meshed into 256×32 elements,and the lower beam is meshed into 512×64 elements.The 60%volume constraint is taken,i.e.,=0.6.

Figure 5:The double-beam model

Compute this example using IGA with the Nitsche’s method to glue the overlapped interface,we plot the displacement distribution in Fig.6,proving the analysis correctness for this example.We reach at the optimization result as shown in Fig.7.It is noted that in paper[14],two kinds of interface conditions are composed,i.e.,tie condition and small sliding condition,since in this paper our objective is to glue both beams by the Nitsche’s method,our results is similar to the optimum results obtained by the interface tie condition, while the small sliding condition is not tested.The stress distribution of the optimized configuration is plotted in Fig.8,and the iteration history of topology optimization is shown in Fig.9.This example is also computed by CPS4R elements (4-node bilinear plane stress quadrilateral, reduced integration, hourglass control) using the SIMP optimization method with a penalty coefficientp= 3 for design variables.Please note that the result from ABAQUS seems smoother than the proposed result,which is mainly due to the post-processing treatment in ABAQUS,and this could be simply described as drawing black solid lines along the material boundaries.In addition, in order to get a relatively good result as reference, in ABAQUS a well refined mesh is adopted.The objective function values of the configuration by the proposed method shown in Fig.7 and by ABAQUS shown in Fig.10 are 382 and 396,respectively,meaning that these two configurations are similar in terms of mechanical characteristics, which once again proves the correctness of the proposed method.

Figure 6:Contour plots of displacement distribution of the double-beam example

Figure 7:Topology optimization result of the double-beam example

Moreover, the curve about the compliance in Fig.9 appears violent oscillation, which can be explained as follows.In order to ensure the approximation accuracy of the integer programming subproblems,a volume fraction reduction factor is adopted to gradually reduce the material consumption as a moving limit strategy, as indicated in Eq.(26).Thus the structural compliance is temporarily improved as the material usage is reduced, and the structural compliance is then effectively reduced owing to the proposed optimization algorithm.Until the convergence criterion is reached,the material usage will be once again reduced, and this process is repeated as explained in the flow chart in Fig.4.Therefore,the value of the object function,i.e.,the compliance,appears oscillation unavoidably throughout the iteration history.Although not adopted in this research, this problem could be effectively addressed as reported in paper [27], in which Liang and Cheng proposed a sequential approximate integer programming method with trust region,in order to directly control the amount of the change of the design variables and thus to stabilize the iterative history significantly.

Figure 8:Contour plot of Von-Mises stress distribution of the optimized double-beam example

Figure 9:Iteration history of the double-beam example

Figure 10:Topology optimization result by ABAQUS of the double-beam example

4.2 Cantilever Beam

As the second example,consider the cantilever beam problem as shown in Fig.11,then verify the results by comparing with the single-patch result and ABAQUS result.The beam size is 100 mm ×50 mm.The cantilever beam is divided into one upper and one lower patches according to the neutral layer, wherein the upper patch is meshed by 128×32 elements, and the lower patch is meshed into 256×64 elements.For better comparison of the effectiveness of the proposed method,the 60%volume constraint,i.e.,=0.6 and 50%volume constraint,i.e.,=0.5 are both examined.

Figure 11:Cantilever beam model made by two patches

The contour plot of displacement distribution by multi-patch IGA is shown in Fig.12,indicating that the Nitsche’s method can effectively glue the two patches.The final optimized material distribution is shown in Fig.13,it can be seen that the material transition at the coupling interface is smooth.The Mises stress distribution is shown in Fig.14, which shows that the optimized structure is uniformly stressed.The iterative curve of the objective function and material usage is shown in Fig.15.As the material usage decreases,the structural stiffness decreases,and when the material usage converges to the target volume fraction,the objective function also converges to a stable value.

Figure 12:Contour plots of displacement distribution of the cantilever beam example

Figure 13:Topology optimization results of the cantilever beam example

Figure 14:Contour plots of Von-Mises stress distribution of the optimized cantilever beam example

In order to illustrate the effectiveness of the proposed method, the result by our method is compared with the results by single patch (Fig.16), and by ABAQUS (Fig.17).It can be seen that these two configurations are similar as Fig.13 by our method.Once again,the result from ABAQUS seems smoother than the proposed result, which is mainly because of its post-processing treatment by drawing solid lines along boundaries as well as the adoption a more refined mesh resolution.The topology optimization results of continuum structures usually depend on mesh generation, thus the density distribution details of these three results are not strictly the same.Furthermore,the values of the objective function correspond to the configurations in Figs.13–17 are respectively 32.9,32.1 and 33.9 for= 0.6, and 37.6, 35.8 and 38.6 for= 0.5.The objective function values of the three configurations are similar, indicating that the mechanical characteristics of the three configurations are similar as well.

The mathematical nature of the topology optimization problem with discrete variables consists in nonlinear integer programming with partial differential equation constraints, and different discretizations will lead to different local optimal solutions.The method proposed in this paper provides a feasible method for mechanical analysis and topology optimization under different meshes.In practical engineering,structures are usually made upon parts assembling,the advantage of using multipatches of NURBS models is that it allows engineers to use independent meshes of different parts for analysis and design,which could greatly reduce the mesh requirements.

Figure 15:Iteration histories of the cantilever beam example

Figure 16:Topology optimization results of the cantilever beam example(single patch)

Figure 17:Topology optimization results by ABAQUS of the cantilever beam example

4.3 One-Quarter Annulus

Finally consider the one-quarter annulus problem as shown in Fig.18.The model is divided into two patches,and the mesh division is 128×32 and 256×64,respectively.In order to show the advantage of IGA in terms of curved meshes,the two patches are connected by curved edges.Both=0.6 and=0.5 are taken in order to verify the effectiveness of the proposed method.

Figure 18:One-quarter annulus model made by two patches

The displacement distribution of the initial configuration is drawn in Fig.19, and the final optimized material distribution is shown in Fig.20.Although the material distribution under the coarse grid is relatively rough due to the grid resolution of the upper patch,it can be seen that the material transition along the coupling interface is still relatively smooth.It can be inferred that the analysis and optimization method proposed in this paper can effectively deal with the connection of curved interfaces.

Figure 19:Contour plots of displacement distribution of the one-quarter annulus example

The contour plot of the Mises stress is shown in Fig.21, indicating that the stress distribution of the structure is generally uniform.The iteration curves of the objective function and the material consumption is plotted in Fig.22.As the material consumption is gradually reduced, the objective function converges to a stable value of 28.2 for= 0.6, and 32.2= 0.5.The optimized results by single patch and by ABAQUS are shown in Figs.23 and 24,the corresponding objective function values are respectively 28.4 and 30.5 for= 0.6, and 31.5 and 34.4 for= 0.5.By comparing of the optimization results and the objective function values,it is concluded that the obtained optimized configuration in Fig.20 is effective.

Figure 20:Topology optimization results of the one-quarter annulus example

Figure 21:Contour plots of von-mises stress distribution of the optimized one-quarter annulus example

Figure 22:Iteration histories of the one-quarter annulus example

Figure 23:Topology optimization results of the one-quarter annulus example(single patch)

Figure 24:Topology optimization results by ABAQUS of the one-quarter annulus example

5 Conclusions

For multi-patch models in engineer practices, it is necessary to analyze and optimize from the perspective of the entire structure.In this paper, the Nitsche’s method is used to perform the isogeometric multi-patch analysis,which allows the patches of independent meshes to be glued with each other, thus greatly reducing the quality requirements for mesh generation.Furthermore, the discrete topology optimization method based on integer programming is introduced to obtain blackwhite boundaries of the conceptual design.Taking advantages of curved-edge elements in IGA,topology optimization in the framework of IGA can obtain locally smoother material boundaries than traditional FEM under the same level of mesh resolution.In addition,a sensitivity filtering method based on the Helmholtz equation is used for complex curved-edge elements,which not only ensures the filtering effect,but also saves the computational resources for searching of adjacent elements.In addition, for the coupling interfaces, an averaging method of interface sensitivities is proposed to ensure the transition smoothness of materials.Finally the performance of the proposed method is tested and verified by numerical examples.

It should be noted that although the topology optimization results in this research is totally black and white, the resulting design is still unavoidably jagged due to the nature of the IGA meshes and the nature of the adopted pixel-based element-wise optimization method.However, we believe that the research results have good potential to be adopted as initial configurations for some optimization methods,for instance to provide initial configuration guesses for MMC[7],MMV[8]and the trimmed method[14],and this could be treated as our further research topics.

Acknowledgement:The authors wish to express their appreciation to the reviewers for their helpful suggestions which greatly improved the presentation of this paper.

Funding Statement:Qingyuan HU is supported by the Fundamental Research Funds for the Central Universities (No.JUSRP12038), the Natural Science Foundation of Jiangsu Province (No.BK20200611),and the National Natural Science Foundation of China(No.12102146).

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