The asymptotic concentration approach combined with isogeometric analysis for topology optimization of two-dimensional linear elasticity structures

1.   Introduction

As an advanced method, structural topology optimization (TO) aims to search and seek the optimal layout of materials in the design domain under the conditions of objective functions and various constraints, and generally involves multidisciplinary knowledge. Meanwhile, as a powerful tool of product design in the initial stage, structural TO can greatly reduce the labor of the designers and shorten the design cycle of the product. In addition to making the product design more efficient, structural TO can make the product meet various design requirements at a lower cost through various simulation experiments. Since Bendsøe and Kikuchi put forward the homogenization approach in 1988 ^[1], structural TO has attracted widespread attention and its related research has been rapidly developed. Many categories of theories and methodologies on TO have been successively proposed, such as the density-based TO method/algorithm ^[2,3,4,5,6], the evolutionary structural optimization (ESO) method ^[7,8,9], the level set method (LSM) ^{[10,11,12,13,14,15]}, the parameterized level set method (PLSM) ^[16,17,18] (it should be noted that parametric level sets have been presented to solve the inverse problem of detecting inclusion interfaces in a piezoelectric structure ^[16] 5–6 years prior to the approaches ^[17,18]), the meshless method ^[19], the phase field method ^[20], the modified guide-weight method ^[21], and so on. Zhou and Rozvany ^[2] and Mlejnek ^[3] proposed a function interpolation model by the solid isotropic material with penalisation (SIMP) scheme within the framework of variable density principle. In the SIMP scheme, a penalty factor is used to polarize the element densities during the iterative process of TO. As the optimization iteration progresses, the intermediate density variable is polarized to either 0 or 1 gradually, which prevents the optimization results from containing as much of the intermediate density as possible, so as to obtain an optimal topological structure with a clear boundary. Because the mathematical model of the SIMP scheme is relatively intuitive and possesses the advantages of high computational efficiency, it has attracted much attention and has been widely applied. Therefore, we carry out the relevant structural TO research under the framework of SIMP scheme in this work.

Structural optimization combined with finite element method is one of the commonly used technical approaches and means in engineering design. In the finite element method, interpolation polynomials are used to approximate various geometric models. To reduce the geometric description error as much as possible, it is necessary to refine the finite element meshes. This process can improve the accuracy of the geometric description. However, the refinement of meshes usually leads to a significant increase in the number of structural elements, which is bound to generate high computational costs ^[22]. Moreover, when the finite element method is used for structural TO, it is usually difficult to obtain smooth boundaries in the optimized results. These problems can be solved by some smoothing treatments. However, it is easy to produce large errors in the process of the smoothing treatment, which may lead to inconsistencies between the calculation model and the description model.

Accordingly, Hughes et al. ^[22] and Cottrell et al. ^[23] put forward the isogeometric analysis (IGA) method based on the integration of computer-aided design (CAD) and computer-aided engineering (CAE) frameworks. Based on the isoparametric element idea, the IGA approach uses spline basis functions to describe geometric models, and the same basis functions are also used in the analysis and design processes. IGA can ensure the accuracy of geometric descriptions, and the structural geometry is independent of the discrete degree of structure. Thus, even in the case of coarsely meshing, the structure geometry can be accurately represented, and the analysis error caused by geometric approximation can be eliminated by the IGA approach ^[24].

The introduction of the IGA approach provides a new tool of structural TO. Due to the excellent characteristics of the IGA approach, researchers have gradually combined the IGA approach with various optimization methods to tackle corresponding TO problems. In structural TO, some scholars combine the SIMP model with the IGA approach because of its intuitive mathematical model, convenient and efficient calculation, and a wide range of applications. Hassani et al. ^[25] updated the design variables in structural TO using the optimization criterion on the basis of the SIMP framework and the isogeometric analysis method. Through a series of optimization examples, the optimization method was verified to be able to suppress the checkerboard phenomenon frequently occurring in the SIMP scheme. In fact, Kumar and Parthasarathy ^[26] used basis spline (B-spline) functions to represent structural geometry and density distribution functions in TO problems prior to the work. It can be observed from the numerical examples that the TO method based on B-spline basis functions can overcome the defects of the checkerboard phenomenon. However, the method cannot eliminate the disadvantages of grid correlation. Subsequently, Qian ^[27] proposed a TO method on the basis of B-spline space, where the optimized structure of arbitrary shape is embedded into a rectangular B-spline space, the density value corresponding to the B-spline control point is used as the design variable, and the density distribution function is represented by the B-spline basis functions. Moreover, it is shown by the method that this density distribution function can play a filtering role in the optimization iteration process. Lieu and Lee ^[28] put forward a multi-resolution TO method on the basis of the IGA approach and the SIMP framework, where high resolution design can be obtained by only adjusting the refinement level of the analysis unit without increasing the amount of calculation. Later, they combined the method with an active phase algorithm to address the TO of multi-material structures ^[29]. Taheri and Suresh ^[30] proposed an isogeometric method for TO with the perimeter constraints and utilized the method to carry on the TO of multi-material structures. Montemurro et al. ^{[31,32,33,34,35]} proposed the CAD-compatible topology optimization methods/algorithms based on non-uniform rational basis spline (NURBS) hyper-surfaces to design structures at multiple scales. As for the TO of composites, Ghasemi et al. ^[36] proposed a design methodology based on a combination of IGA, level set, and point wise density mapping techniques for the TO of piezoelectric/flexoelectric materials. The fourth order partial differential equations (PDEs) of flexoelectricity, which require at least C¹ continuous approximations, were discretized using NURBS. The point wise density mapping technique with consistent derivatives was directly used in the weak form of the governing equations. The boundary of the design domain was implicitly represented by a level set function. The accuracy of the IGA model was confirmed through numerical examples, including a cantilever beam under a point load and a truncated pyramid under compression with different electrical boundary conditions. Additionally, Ghasemi et al. ^[37] proposed a computational design methodology for the TO of multi-material-based flexoelectric composites. In the methodology, they adopted the multiphase vector level set (LS) model, which easily copes with various numbers of phases, efficiently satisfies multiple constraints, and intrinsically avoids either overlap or vacuum among different phases. They extended the point wise density mapping technique for multi-material design and used the B-spline elements to discretize the PDEs of flexoelectricity. The dependence of the objective function on the design variables was incorporated using the adjoint technique. The obtained design sensitivities were used in the Hamilton-Jacobi (H-J) equation to update the LS function. They provided numerical examples for two, three, and four phase flexoelectric composites to demonstrate the flexibility of the model, as well as the significant enhancement in electromechanical coupling coefficients that can be obtained using multi-material TO for flexoelectric composites. Wang et al. ^[38] proposed a multi-scale isogeometric method for the TO with high computational efficiency by combining the IGA approach and the evolutionary homogenization method. Gai et al. ^[39] proposed an explicit isogeometric TO approach based on moving morphable voids (MMVs) with closed B-spline boundary curves. They modeled the design domain by a NURBS patch, and then employed the NURBS-based IGA for a structural response and the well-established adjoint approach for sensitivity analysis. The IGA method based on the SIMP framework can overcome all shortcomings brought by the finite element meshes. Nevertheless, gray-scale interval problems still exist on boundaries of TO structures obtained by the method.

In this work, an asymptotic concentration method is proposed for structural TO with the B-spline entities and the SIMP framework. We employ the asymptotic concentration method to eliminate the gray-scale interval phenomenon within the SIMP framework for TO, while taking full advantage of the IGA approach. Here, only the pseudo-density values at the control points are considered as design variables, and the weights are not optimized. The optimization criterion method is used to update the design variable, and the asymptotic concentration method is used to update the element density variable to solve the TO problem of minimizing structural compliance under a volume constraint. Finally, several typical numerical examples are given to verify the optimization efficiency and feasibility of the proposed method.

The rest of this paper is organized as follows. The IGA formula of the plane elasticity problem is described in Section 2. The mathematical model of the asymptotic concentration approach combined with IGA for the TO of two-dimensional (2D) linear elasticity structures under the B-spline-SIMP framework is constructed in Section 3. Different values of parameters are used for optimization to study the influence of parameter values on TO in Section 4. The effectiveness of the proposed method is illustrated by typical numerical examples in Section 5. Finally, the summary is provided in Section 6.

2.   Formulation of isogeometric analysis for plane elasticity problems

According to the isoparametric element notion in finite element analysis, the IGA approach uses NURBS basis functions to describe structural geometric models instead of interpolation polynomials. There is no description error when describing the geometric model by the IGA approach. The basic principle of the IGA approach is that basis functions for geometric modeling are the same as those used for structural analysis. That is to say, the information representing geometry in CAD is used for mechanical calculations, so that the structural analysis with high-precision geometric models can be realized. NURBS is the most widely used spline in computing and modeling technology of IGA. NURBS can accurately express the conic curve, thereby making the free-form surface modeling convenient. In addition, NURBS can be combined with many efficient and stable algorithms. This study is limited to B-spline entities as opposed to general approaches using NURBS surfaces and hyper surfaces as in ^{[31,32,33,34,35]} (i.e., we are only considering B-spline entities and not NURBS entities to describe the pseudo-density field over the design domain).

2.1.   Description of the basis functions

The nodal vector of the B-spline curve in one-dimensional space, expressed by $[I] = \left\{{\xi }_{\text{1}}\text{, }{\xi }_{\text{2}}, \mathrm{...}, {\xi }_{n+\text{p}+1}\right\}$, is a monotone non-decreasing sequence. In the expression, ${\xi }_{\text{i}}$ represents the i-th node; i (i = 1, 2, …, n + p + 1) indicates the index coordinates of the nodes; p indicates the spline order; and n represents the total number of basis functions. If any two adjacent nodes are equidistant, the corresponding node vectors are entitled uniform node vectors. If the values of the first node and the last node are repeated (p + 1) times, the node vectors are named as open node vectors. In this work, we use the uniform open node vectors, which are expressed as follows:

In Eq (1), the values of a and b are equal to 0 and 1, respectively, and the values of a and b are repeated (p + 1) times.

Given a node vector, when the spline order p equals 0, the basis functions of B-spline are the Bernstein's polynomials, defined recursively as follows:

Additionally, given a node vector, when the spline order p is greater than 0, the basis functions of B-spline are the Bernstein's polynomials, which can be defined as follows:

The non-zero second-order B-spline basis functions described on the prescribed node vector $[I] = \left\{\text{0}, \text{0}, \text{0}, \text{1}, \text{2}, \text{3}, \text{4}, \text{4}, \text{5}, \text{5}, \text{5}\right\}$ are shown in Figure 1.

2.2.   Refinement of computational domain

IGA is a numerical calculating approach, and its results are approximate to those of the original problem. To improve the calculation accuracy, the model of IGA usually needs to be refined. In the IGA approach, there are three refining methods, namely, the h refining method, the p refining method, and the k refining method. The first two refining methods are similar to those used in finite element method (FEM), whereas the third refining method, namely the k refining method, is unique to IGA. Therefore, the k refining method is briefly described as follows.

At present, k refinement has become the most widely used refining method in IGA. It makes full use of the property of non-interference between order increment and node insertion. There are two different algorithmic flows for k refinement: one is to insert the nodes first and then increase the order; the other is to increase the order first and then insert the nodes. At the same parameter points, the continuity of the curve obtained by the latter algorithm is higher than that obtained by the former algorithm.

Figure 2 shows the two different algorithmic flows of k refinement: the flow on the left side (from Figure 2(a) to Figure 2(b) and then to Figure 2(d)) represents the algorithm of inserting the nodes first and then increasing the order; the flow on the right side (from Figure 2(a) to Figure 2(c) and then to Figure 2(e)) represents the algorithm of increasing the order first and then inserting the nodes. According to Figure 2, one can conclude that the continuity of the basis function obtained by the algorithm of node insertion after order increment is higher than that obtained by the algorithm of order increment after node insertion.

2.3.   Solving planar elasticity problem based on the basis functions

To solve the planar elasticity problem, the B-spline-based approach under the SIMP framework is adopted in the TO of two-dimensional structure. Wherein we set a control point density for each control point and take each control point density as the design variable, the density of element e can be represented by the relevant control points as follows:

In Eq (4), (ε, η) indicates the centroid coordinates of element e; ${N}_{I}(\varepsilon , \eta )$ indicates the values of the basis functions at the centroid of element e; ${\rho }_{I}$ represents the density of control points relating to the basis function value ${N}_{I}(\varepsilon , \eta )$; and (p + 1)(q + 1) indicates the total number of the basis functions and the control points that support element e.

The algebraic equations of the two-dimensional linear elastic problem under static load are expressed as follows:

In Eq (5), K indicates the global stiffness matrix; U and F represent the global displacement and external force vectors at the control point of the structure, respectively; n indicates the total number of discrete elements in the parameter space; and the element stiffness matrix ${{\bf{k}}}_{\text{e}}$ can be expressed as follows:

In Eq (6), ${\Omega }_{\text{e}}$ represents the design area of element e; t represents the thickness of element e; J represents Jacobian determinant of geometrical mapping between structural measures ^[10]; and B represents strain displacement matrix of structural element, which is expressed as:

where, ${N}_{\text{i, x}}$ indicates the first order variation of the ith B-spline basis function with respect to the x coordinate; ${\boldsymbol{N}}_{\text{j, y}}$ indicates the first order variation of the jth B-spline basis function with respect to the y coordinate; and D indicates the element elastic constant matrix related to planar stress problem, which can be described as:

where, ν represents Poisson's ratio of structural material and E represents Young's modulus of the structural material.

3.   Mathematical model of the proposed method

An asymptotic concentration approach combined with isogeometric analysis is proposed for structural TO under the B-spline-SIMP framework, and the mathematical model of the proposed method is described in this section.

3.1.   Description of density-based topology optimization

In the traditional SIMP scheme ^[40], the relationship between the whole elasticity matrix of the materials and the density design variable are expressed as follows:

where, ${{\bf{D}}}_{\text{0}}$ represents the whole elasticity matrix of solid material; ${\rho }_{\text{e}}$ signifies the density design variable; and s signifies the penalty factor, and its value generally equals 3 so as to force intermediate density design variables to either 0 or 1. However, when the value of the density design variable equals zero, the corresponding Young's modulus also equals zero, which will lead to the ill-conditioned singularity of the structural stiffness matrix during the optimization calculation process. Hence, we utilize the improved SIMP interpolation model, wherein the relationship between the density design variable and Young's modulus in structural TO can be expressed as ^[41]:

where, ${E}_{\text{0}}$ represents Young's modulus of solid material and ${E}_{\text{min}}$ is set as a value slightly greater than zero, indicating Young's modulus of structural material when density design variable equals zero in structural TO.

According to Eq (8), the elastic constant matrix of solid materials for the planar stress problem can be described by the following formula:

In general, structural TO aims to search for and seek the optimal layout of materials in the design domain under certain constraint conditions, such as volume and stress constraints, so as to obtain the optimal value of structural stiffness. When structural stiffness reaches its maximum value, the total structural strain energy will reach its optimal value accordingly. Correspondingly, the structural compliance reaches its minimum value. In this work, structural compliance minimization is taken as the objective function of TO. According to Eq (5), the calculation formula of structural compliance is given as:

It should be noted that the definition of the compliance provided in this section and used in the rest of the manuscript is not the most general one and can be used only when the boundary conditions of the Dirichlet type (i.e., on generalised displacements) are zero. As widely discussed in ^[31], in the most general case of inhomogeneous Neumann-Dirichlet boundary conditions, the compliance does not coincide with neither the strain energy nor with the work of internal forces (which is twice the strain energy), but it is related to the total potential energy of the continuum.

Therefore, the TO model in this work can be expressed by the following mathematical formula:

In Eq (13), ${\rho }_{\text{e}}$ refers to the pseudo density of element e, the value of which belongs to the interval [0, 1]; F represents the forces applied on the structure; n represents the total number of finite element meshes after discretization of the structural design domain; ${{\bf{u}}}_{\text{e}}$ denotes the displacements of element e; U represents the displacements of the structural control points; K represents the stiffness matrix of the whole structure; ${k}_{\text{e}}$ represents the stiffness matrix of element e; ${V}_{\text{e}}$ denotes the volume of element e; $V$ denotes the total volume of the structural design domain; f represents the given volume fraction during structural optimization; and g represents the constraint function to ensure that the volume of solid material meets prescribed design requirements. $V$ and ${V}_{\text{e}}$ are respectively expressed as follows:

where, $\Omega$ denotes the global design area of the structure and ${\Omega }_{\text{e}}$ denotes the local design area of element e.

3.2.   Asymptotic concentration scheme

When the SIMP scheme is used for TO, there exists a gray-scale interval represented by the intermediate density element in the obtained optimal result. Even though the IGA approach under the SIMP framework is used for structural TO, the gray-scale interval still exists in the optimization results. To eliminate the gray-scale interval in the obtained optimization structure, the idea of asymptotic concentration ^[42] is employed in this work. The core point of this method is to gradually concentrate the material density variable to the specified variable value through multiple iterations. By virtue of the asymptotic concentration approach, all of the element density values are gradually concentrated from the set intermediate values to 0 and 1, so that the intermediate density values in the optimization results can be eliminated. The direction of the asymptotic concentration is from the specified value of the asymptotic intermediate density to the initial density of the material at the ends of both left and right. Meanwhile, the speed of the asymptotic concentration is determined by adjusting the values of the asymptotic parameters er and rr.

When the asymptotic concentration approach is adopted for structural TO, the asymptotic concentration process of the element density variables is described in Figure 3. Wherein, the number 0 represents element density of the void material phase; the number 1 represents element density of the solid material phase; the symbol m represents the specified intermediate density to be asymptotically concentrated; the symbols a and b represent two coefficients of asymptotic concentration, and the ratio of a to b is equal to that of the length of interval [0, m] to the length of interval [m, 1]; and the symbol rr is a gradually increasing asymptotic parameter with iterations in the structural TO, and its initial value is selected as 0. The asymptotic parameter is equal to the increment of rr in each iteration, so the asymptotic parameter represents an increasing rate of the parameter rr during optimization iterations. The density value within a certain distance of the intermediate density m is asymptotically concentrated on two specified values in the two directions of m. The distance is determined by the intermediate density m, the asymptotic parameter rr, and the asymptotic concentration coefficients a and b. Specifically, the element density values within the interval [m-a^*rr, m] will be concentrated on the lower limit of the interval, and the element density values within the interval [m, m+b^*rr] will be concentrated on the upper limit of the interval by the asymptotic concentration approach. The asymptotic parameter rr gradually increases with the number of iterations. As a result, the length of the density interval (m-a^*rr, m+b^*rr) affected by the asymptotic concentration gradually increases. Since the density value ranges from 0 to 1, when the density interval affected by the asymptotic concentration increases to a certain degree, there will be no elements with a density value other than 0 or 1, thereby eliminating the gray interval generated by the intermediate density. Therefore, it is crucial to reasonably select the asymptotic parameters a, b, and rr. The criterion of the asymptotic concentration approach to update the element density variable can be expressed as follows:

where, ${\rho }_{\text{e}}$ and ${\rho }_{\text{e}}^{\text{new}}$ represent the element densities before and after updating with the asymptotic concentration approach, respectively; m refers to the intermediate density; rr represents asymptotic parameters; and a and b represent the asymptotic coefficients. During the asymptotic concentration, the relationship of the asymptotic coefficients can be described as:

where, $\delta$ is a prescribed constant with a value equal to one.

3.3.   Sensitivity analysis

For 2D stress problem, the global stiffness matrix of the structure is expressed by element density as follows ^{[31,32,33,34,35,43]}:

where ${{\bf{k}}}_{e}$ is the element stiffness matrix, which is transformed into:

Since the element densities are represented by the B-spline basis functions and the control points in Eq (4), the element stiffness matrix is converted into a formula expressed by the densities of the control points by substituting Eq (4) into Eq (19). The formula is described as follows ^[43]:

By using the same method, the volume formula based on the densities of the control points is transformed into:

When the gradient based method is applied to the structural TO, it is always necessary to calculate the sensitivity of the objective and constraint functions with respect to the design variables. In particular, in order to obtain the sensitivity of structural compliance, the sensitivity of the element stiffness matrix with respect to the density of control points is calculated, which is expressed as ^[43]:

According to the TO model of Eq (12), the sensitivity of the objective function c and the element volume V_e with respect to the density of control points are respectively described as ^[43]:

The sensitivity of g in Eq (12) with respect to the density of control points is expressed as follows:

When the density-based method is applied to the structural TO, the checkerboard phenomenon usually appears in the optimization results. This is caused by the situation where the element density values 0-1 alternately appear in some areas of the optimal structure. In the process of density-based continuous TO, the checkerboard phenomenon will inevitably occur, which has a great impact on the optimal structures. The distribution of materials needs to be selected in the area where the element density values of 0-1 alternate. Even so, it is difficult to guarantee that the final optimized structure is the optimal structure.

When the traditional SIMP approach is only used for TO, an external filter is usually required to suppress the appearance of the checkerboard phenomenon. It is necessary to select the filter type suitable for the optimization method, so the addition of density filtering will inevitably cause an increase in the calculation amount. However, as for the B-spline basis function, the continuous density domain represented by the B-spline can provide a built-in filter for TO due to its excellent characteristics ^[44]. The principle is that B-spline basis functions have characteristics supported by local elements, and the density distribution function expressed by the control points and B-spline basis functions has the same representation form as density filtering. The built-in filter radii of the B-spline basis functions are affected by B-spline order, as well as by node span. Meanwhile, this feature also applies to NURBS. As discussed in ^{[31,32,33,34,35]}, in the context of the density-based TO algorithm reformulated in the context of NURBS entities, there is no need to introduce an artificial filter because the local support property of the basis functions acts as an implicit filter, thus the checkerboard effect and the mesh dependency of the solution can be strongly reduced or even avoided. Moreover, as widely discussed in ^[45], the integer parameters of the NURBS (i.e., number of control points and degrees of the Bernstein's polynomials) can be properly tuned to control the size of the local support and, hence, the minimum length scale requirements (without introducing an explicit constraint in the problem formulation).

3.4.   Corresponding optimization algorithm

A lot of existing technical means can be used to solve topological optimization problems, including the optimization criterion (OC) approach ^[7], the method of moving asymptote (MMA) ^[46], the sequential linear programming (SLP) method ^[47,48], and so on. In the OC approach, optimization criterion are derived according to the Karush-Kuhn-Tucker (KKT) conditions and formed the OC algorithm, through which the optimal solution can be obtained via multiple iterations. When it is used to tackle the single-constraint optimization problems, the OC approach either has a high speed of convergence or requires a small number of iterations. By virtue of this advantage, this work uses the OC approach to update the design variables.

In this method, we take densities of the control points as design variables, so that design variables updated by the OC algorithm are densities of the control points rather than densities of elements. The specific update criteria are defined as follows:

where, the parameter mov denotes a preset positive constant with the value of 0.2 used to indicate the speed of density change. The parameter η represents the damping coefficient with the prescribed value of 0.5. The calculation formula of B_e value is defined as follows ^[42]:

In Eq (27), λ denotes Lagrangian multiplier, the value of which can usually be calculated by the bisection method.

3.5.   Numerical implementation

The numerical implementation of the asymptotic concentration approach combined with IGA for TO under the B-Spline-SIMP framework is described as follows:

1) Input the node vector and control point vector needed to construct the geometric structure, and use the B-spline functions to construct the initial geometric model;

2) Construct the IGA grid, and refine the grid with the unique k refinement method of IGA;

3) Prescribe the structural boundary and load conditions;

4) Set the value of initial control point density to volume fraction, that is ${\rho }_{i, j} = f$, and use the control point densities and the B-spline basis functions to describe element densities. In this work, the densities of all points in an element are set to the same value, and the element densities are expressed via the centroid densities of the elements;

5) Calculate the local stiffness matrix of each element by using the B-spline basis functions, and assemble the structural global stiffness matrix according to local stiffness matrix of each element;

6) Perform the finite element analysis: introduce Young's modulus expressed by element densities into the structural stiffness matrix, calculate the displacement matrix, and then solve the objective function and volume fraction, as well as their sensitivity information;

7) Update control point densities by the OC method, represent the element density variables by using the B-spline basis functions and control point densities, and then update the element density variables by using the asymptotic concentration approach proposed in this work;

8) Judge whether the optimized structure satisfies the convergence condition: if optimization iteration meets convergence requirement, terminate the loop of optimization iteration and meanwhile output optimization result; otherwise, return to Step 6) to continue optimization iteration.

The flow chart of TO using the proposed method is shown in Figure 4.

4.   Selection of parameter value in topology optimization

4.1.   Influence of asymptotic parameters on optimization results

Here, we provide a numerical example involving the TO of planar bracket structure to analyze the influence of different asymptotic parameters on optimization results.

The design domain, boundary, and load conditions of this example are shown in Figure 5. As for the design domain with the aspect ratio 2:1, a downward force equal to 1 is applied at the midpoint on the right side, whereas the left side is fixed. The entire design domain is discretized to 50 × 100 quadrilateral elements by virtue of the k refining method. The values of parameters p and q that determine the order of the B-spline basis functions are both set as 1. The thickness of the bracket structure is set as $t = 1$. Young's modulus associated with the solid material is selected as ${E}_{0} = \text{1}$. Young's modulus associated with void element is selected as ${E}_{\text{min}} = 1e-8$. Poisson's ratio is set as $\gamma = \text{0}\text{.3}$. The volume fraction is selected as $f = \text{0}\text{.4}$. In this example, several different asymptotic parameter values are used to analyze the influence of asymptotic parameters on the TO result. Moreover, in the result of this and the following examples, the solid material phase and the void material phase are represented by black and white, respectively, and the initial densities of the solid material phase and the void material phase are set to 1 and 0, respectively.

Table 1 displays the results of the TO of the planar bracket structure with different values of asymptotic parameter er, by which we can compare the effects of different asymptotic parameter values on the optimal topology structure, compliance value curve, number of iterations, and optimal compliance value. It is observed from Table 1 that the optimal topological structure obtained by the proposed method is almost the same when the asymptotic parameter er takes different values of 0.005, 0.007, 0.009, 0.011, and 0.013, respectively. In addition, the compliance value is close to the optimal compliance value after just 10 iterations. It can also be found from Table 1 that the value selection of the asymptotic parameter er for the same structure has little influence on the TO results. Specifically, the value selection of the asymptotic parameter er has only a small impact on the optimal compliance, and almost no influence on the optimal topology structure. Moreover, it can be seen from Table 1 that when the values of er are 0.005 and 0.007, respectively, the optimal compliance values obtained by the proposed method are almost the same; when the values of er are greater than 0.007, the optimal compliance value obtained by the proposed method gradually increases with the value of er; and the number of optimization iterations gradually decreases when the value of er gradually increases. Therefore, when selecting the values of asymptotic parameters, we should consider the balance between the calculation amount and the optimal compliance value. That is to say, for the selected values of asymptotic parameters, a relatively large number of optimization iterations are required to obtain a relatively low optimal compliance value; if completing the optimization with a relatively small number of optimization iterations, one will obtain a relatively high optimal compliance value.

4.2.   Influence of spline order on optimization results

In the asymptotic concentration approach for the TO under the B-spline-SIMP framework, the size of the spline order affects the number of control points and the size of the support interval of the basis functions, which in turn affects the computational complexity. Here, a numerical example involving the TO of two-point fixed beam structure is used to analyze the influence of different sizes of spline order on optimization results.

The design domain, boundary, and load conditions of this example are shown in Figure 6. The aspect ratio of the design domain is 2:1. The midpoints on the left and right boundaries of the structure are fixed, and a downward force equal to 1 is acted at the center of the structure. The entire design domain is discretized to 90 × 45 quadrilateral elements by virtue of the k refining method. The value of the asymptotic parameter er is selected as 0.007. The thickness of the beam structure is set as $t = 1$. Young's modulus associated with the solid material is selected as ${E}_{0} = \text{1}$. Young's modulus associated with the void element is selected as ${E}_{\text{min}} = 1e-8$. Poisson's ratio is set as $\gamma = \text{0}\text{.3}$. The volume fraction is selected as $f = \text{0}\text{.4}$. In this example, three different orders of B-splines, namely the first-order B-spline, the second-order B-spline, and the third-order B-spline are used to analyze the influence of B-spline order on optimization results.

Figure 7 illustrates the optimal topology results of the two-point fixed beam with different orders of B-splines. Figure 8 shows the compliance change curves corresponding to the optimization process of the two-point fixed beam with three different orders of B-splines. It is observed from Figure 7 that the optimal topology structures obtained by the proposed method are similar with three different orders of B-splines. Regarding the change of the compliance value during the optimization process, it can be found from Figure 8 that the compliance almost reaches its optimal value within 10 iterations. In other words, the compliance value changes very little after 10 iterations with all three orders of B-splines. For comparison, the data involved in the TO of two-point fixed beam with three different orders of B-splines are listed in Table 2. By comparing the compliance values corresponding to different spline orders, it is found that the optimal compliance value gradually increases with the order of splines. Nevertheless, the total number of optimization iterations seems hardly affected by the order of splines.

It can be found from this example that the required number of optimization iterations seems almost the same and similar optimal topology results can be obtained with three different orders of B-splines by the proposed method. Nonetheless, the use of different orders of B-splines has a great impact on the optimal compliance value. Therefore, for the relatively simple structures, the relatively small compliance value can be obtained by choosing a relatively small spline order; for the complex structures with high requirements for continuity, a large spline order should be selected to meet the requirements. However, the compliance value also gradually increases with the increase of the spline order. Therefore, it is necessary to select an appropriate spline order for the TO according to the geometric complexity of structure.

5.   Numerical examples

Next, the feasibility and optimization efficiency of the proposed method will be verified by a couple of typical numerical examples. In all of these examples, the minimum compliance is taken as an objective function, and structural TO is performed under the constraint of volume fraction by the algorithm proposed in this work.

5.1.   Two-corner-fixed bridge structure

Here, an example involving the TO of two-corner-fixed bridge structure is given to analyze the influence of mesh refinement on the optimization results obtained by the proposed method. The design domain, boundary, and load conditions of this example are shown in Figure 9. The left and the right lower corners of the bridge structure are both fixed, whereas a downward force equal to 1 is applied at the midpoint on the lower side of the bridge structure. As for the design domain, the aspect ratio is 8:3. To analyze the influence of mesh refinement on the optimization results, the entire design domain is discretized by the k refining method to three types of quadrilateral meshes, namely 80 × 30, 96 × 36, and 120 × 45. The value of the asymptotic parameter er is set to 0.007. The thickness of the bridge structure is set to $t = 1$. Young's modulus associated with the solid material is selected as ${E}_{0} = \text{1}$. Young's modulus associated with the void element is selected as ${E}_{\text{min}} = 1e-8$. Poisson's ratio is set as $\gamma = \text{0}\text{.3}$. The volume fraction is selected as $f = \text{0}\text{.4}$.

In this example, three different degrees of mesh refinement are used for the TO to analyze the influence of mesh refinement on the optimal topological structures obtained by the proposed method. Figures 10(a)–(c) show the optimal topological structures of the two-corner-fixed bridge with different types of quadrilateral meshes obtained by the proposed method, respectively. It is observed from Figure 10 that similar optimal structures can be obtained by using meshes with different degrees of refinement, and the obtained optimized structures contain no gray-scale intervals. Nonetheless, different degrees of mesh refinement lead to different degrees of smoothness of the optimized structure boundary. The finer the meshes are, the weaker the jagged structure boundary is, and correspondingly the smoother the structure boundary is; on the contrary, when the meshes are coarse, the jagged boundary is relatively obvious, and the smoothness of the structural boundary is relatively weak.

The clear advantages and characteristics of the proposed method are verified by this example: even at the level of coarse mesh refinement, the correct and reasonable optimal structure can be obtained, and the obtained optimal structure does not contain a gray-scale interval.

5.2.   Cantilever beam structure

Here, we provide another numerical example involving the TO of a cantilever beam structure to substantiate the effectiveness of the method proposed in this work. The design domain, boundary, and load conditions of this example are shown in Figure 11. As for the design domain, the aspect ratio is 2:1. The cantilever beam is fixed on its left side, whereas a downward force equal to 1 is acted at the midpoint on the right side of the cantilever beam structure. The entire design domain is discretized to 100 × 50 quadrilateral meshes by virtue of the k refining method. The value of the asymptotic parameter er is selected as 0.007. The thickness of the beam structure is set to $t = 1$. The B-spline orders are set as $p = 1$ and $q = 1$. Young's modulus associated with the solid material is selected as ${E}_{0} = \text{1}$. T Young's modulus associated with the void element is selected as ${E}_{\text{min}} = 1e-8$. Poisson's ratio is set as $\gamma = \text{0}\text{.3}$. The volume fraction is selected as $f = \text{0}\text{.4}$. This example is composed of two numerical cases as follows.

In the first case, the proposed method, namely the B-spline-SIMP based asymptotic concentration method, is used to accomplish the TO of a cantilever beam structure. When the convergence condition is satisfied after 102 iterations, the volume fraction reaches 0.400. At the same time, the value of the corresponding structural compliance reaches 75.0760. In this case, Figure 12 shows the obtained optimal topological structures. It is observed from Figure 12 that the obtained optimal topological structures are in line with those obtained by other classical methods. After about 25 iterations, the TO result has already been similar to the optimal structure in the frame. As the optimization iteration continues, the structure boundary becomes gradually clear, and the objective function moves towards the stable direction until the convergence condition is satisfied. In Figures 12(a)–(c), the TO results contain gray-scale intervals. In other words, there is still intermediate density in the optimization results. The reason is that the influence area of the asymptotic concentration cannot fully cover the density interval. However, as the TO iteration progresses, the area of influence in the asymptotic concentration gradually expands. Consequently, the intermediate density will be concentrated, thereby causing the gray-scale interval to gradually decrease, until the area of influence in the asymptotic concentration covers the entire density interval, and the intermediate density is completely polarized to either 0 or 1. The final TO result shown in Figure 12(d) does not contain the gray-scale interval any more.

In the second case, the TO of a cantilever beam structure is respectively accomplished by using the proposed method and the existing SIMP scheme based on B-spline for comparison.

In this case, Figure 13 shows the obtained optimal topological structures. It is observed from Figure 13 that the obtained optimal topological structures by using the two above-mentioned methods are basically the same. Nevertheless, the obtained optimal topological structures by using the two methods have obvious differences in terms of clarity. The reasons for this difference are described as follows. The obtained optimal topological structure by using the proposed method has no intermediate-density element, so its boundary is completely clear. However, the obtained optimal topological structure by using the existing SIMP scheme based on B-spline contains a certain number of intermediate-density elements, so its boundary is fuzzy.

Figure 14 shows compliance changes in this case. It is observed from Figure 14 that the compliance value 75.0760 is obtained by the proposed method when the corresponding requirements for convergence are met after 102 iterations, and the compliance value 79.3368 is obtained by the existing SIMP scheme based on B-spline when the corresponding requirements for convergence are met after 118 iterations.

From the aforementioned TO results in this case, it is observed that the proposed method can obtain a smaller compliance value and a clearer structure boundary than the existing SIMP scheme based on B-spline under the same conditions. Furthermore, it is also observed that the required iterations by the proposed method are fewer than those by the existing SIMP scheme based on the B-spline for the same TO problem. Thereupon, the advantages of the proposed method are verified through this example.

5.3.   Multi-load bridge structure

Next, we use another example involving the TO of a multi-load bridge structure to affirm the superior performance of the proposed method. The design domain, boundary, and load conditions of this example are shown in Figure 15. As for the design domain, the aspect ratio is 8:3. The bridge structure is only supported at the lower left and right corners. The former is completely fixed, whereas the latter is simply supported. Three downward external forces with values of 2, 1, and 1 are respectively applied to the midpoint, quarter point, and three-quarter point of the lower edge of the structure.

In this example, Young's modulus associated with the solid material is selected as ${E}_{0} = \text{1}$. Young's modulus associated with void element is selected as ${E}_{\text{min}} = 1e-8$. Poisson's ratio is set as $\gamma = \text{0}\text{.3}$. The volume fraction is selected as $f = \text{0}\text{.4}$. The entire design domain is discretized to 120 × 45 quadrilateral elements by virtue of the k refining method. The B-spline orders are set as $p = 1$ and $q = 1$. The thickness of the bridge structure is set to $t = 1$. The asymptotic parameter is selected as er = 0.007.

Figure 16 shows the optimized topological structures of multi-load bridge by the existing SIMP scheme based on B-spline. All the topological structures shown in Figure 16 have fuzzy boundaries, which are caused by the intermediate-density elements in the TO process. As the process of the TO proceeds, the intermediate density will be gradually polarized, the number of intermediate-density elements will decrease, and the structure boundary will become relatively clear accordingly. Nevertheless, the existing SIMP scheme based on B-spline cannot completely polarize the intermediate density to either 0 or 1. Therefore, there are still some intermediate-density elements in the final topological structure, so the structure boundary contains a certain degree of gray-scale interval. Figure 17 shows the optimized topological structures of the multi-load bridge by the proposed method. The topology structures shown in Figures 17(a)–(c) have fuzzy boundaries, which is due to the fact that the intermediate-density element has not yet been fully polarized during the iteration process. The topology structure shown in Figure 17(d) has a clear structural boundary, which signifies that the proposed method can completely polarize the intermediate density value to 0 or 1. That is to say, when the requirements for convergence are met, the optimization result does not contain the intermediate-density element anymore. Hence the structure boundary is completely represented by the elements with a density value of either 0 or 1.

Figure 18 shows the compliance changes for the TO of the multi-load bridge structure obtained by using the existing SIMP scheme based on B-spline and the proposed method, respectively. It is observed from Figure 18 that the compliance value 254.9795 is obtained by the proposed method when the corresponding requirements for convergence are met after 104 iterations, and the compliance value 267.4641 is obtained by the existing SIMP scheme based on B-spline when the corresponding requirements for convergence are met after 150 iterations.

From the aforementioned TO results obtained by the two methods, it is observed that the proposed method can obtain a smaller compliance value and a clearer structure boundary than the existing SIMP scheme based on B-spline under the same conditions. Furthermore, it can also be found that the required iterations by the proposed method are fewer than those by the existing SIMP scheme based on the B-spline for the same TO problem. Thereupon, the advantages of the proposed method are verified once more through this example.

6.   Conclusions

We propose an asymptotic concentration approach combined with IGA for the structural TO of 2D linear elasticity structures under the B-spline-SIMP framework. In this approach, the control point density is first updated by using the optimization criterion method, and then the element density represented by the B-spline basis function and control point density is polarized by using the asymptotic concentration idea, so as to gradually concentrate all the element density values on either 0 or 1. In the process of TO, with the increase of iteration number, the element density values are gradually polarized to the solid material phase density or the void material phase density, thereby resulting in a final structural element density value of either 0 or 1. Thus, the phenomenon of the gray-scale interval, which leads to fuzzy structure boundary, is completely avoided. The good performance of the proposed method is illustrated by several numerical examples involving the structural TO. Compared with the traditional SIMP scheme, this method can accurately describe the structural geometry by virtue of IGA in the case of a coarse mesh, and the optimal topology structure obtained by this method does not contain the gray-scale intervals. Furthermore, compared with the existing SIMP scheme based on B-spline, this method can obtain a relatively clear structural boundary and a relatively small compliance with relatively few numbers of iterations in the same structural TO problem.

Acknowledgments

This work was supported by Natural Science Foundation of Shaanxi Province (2020JM-200) and China Scholarship Council (201506965015). The authors would like to thank the editors and the reviewers for their valuable comments, which have improved the quality of the manuscript.

Conflict of interest

The authors declare there is no conflict of interest.