Sequentially coupled gradient-based topology and domain shape optimization

A coupled topology and domain shape optimization framework is presented that is based on incorporating the shape design variables of the design domain in the Solid Isotropic Material with Penalization topology optimization method. The shape and topology design variables are incrementally updated in a sequential fashion, using a staggered numerical update scheme. Non-Uniform Rational B-Splines are employed to parameterize the shape of the design domain. This not only guarantees a highly accurate description of the shape boundaries by means of smooth basis functions with compact support, but also enables an efficient control of the design domain with only a few control points. Furthermore, the optimization process is performed in a computationally efficient way by applying a gradient-based optimization algorithm, for which the sensitivities can be computed in closed form. The usefulness of the coupled optimization approach is demonstrated by analyzing several benchmark problems that are subjected to different types of initial conditions and domain bounds. The variation in simulation results denotes that a careful construction of the initial design domain is necessary and meaningful.


Introduction
Topology optimization is a mathematical method that optimally places the material within a given design domain for the specific loading and boundary conditions applied. It has been recognized as the most active research area in the field of structural optimization over the past two decades (Deaton and Grandhi 2014). Driven by a major interest from both academia and industry, various topology optimization approaches have been developed (Garreau et al. 2001;Sigmund and Maute 2013;van Dijk et al. 2013;Huang and Xie 2010;Zhang et al. 2017). The density-based methods, especially the so-called SIMP (Solid Isotropic Material with Penalization) method (Sigmund 2001;Bendsøe and Sigmund 2003), are the most widely used techniques. They start from a discretization of the prescribed design domain into a finite element model, whereby a design variable, which is the relative density, is assigned to each finite element. The relative density can take values between zero (a void) and unity (a solid), and is directly related to the element stiffness via a penalization factor. Accordingly, the optimized distribution of the relative density can be determined by adopting an optimization algorithm that computes the minimal structural compliance (which corresponds to the maximal structural stiffness) of the design domain, see also Andreassen et al. (2011), Liu and Tovar (2014), Stolpe (2015, 2016) and references therein.
Due to the nature of the SIMP method, however, the possible solutions of the density distribution are bounded by the specific size and shape of the design domain chosen. This limitation can be reduced by choosing the design domain as large as possible, although this may lead to relatively large computational times. Alternatively, the construction of the design domain may be formulated as an optimization problem (Kim and Kwak 2002). Such a strategy not only offers the possibility to keep computational times more manageable, it also facilitates the optimization of aeroelastic structures, such as wind turbine blades (Buckney et al. 2013) and aircraft wings , whereby the applied loads (i.e., the aerodynamic forces) are characterized by the actual size and shape of the design domain.
Previous research has demonstrated that the variation of the design domain within an optimization strategy can be accomplished by means of several approaches. A first possible approach is based on the so-called design space optimization method (Kim and Kwak 2002;Jang and Kwak 2008), which adjusts the design domain boundary during topology optimization by adding new design elements. This approach seems to be relatively straightforward to implement, but has the drawback that the boundaries of the domain are represented by a jagged geometry characterized by the element boundaries of the finite element discretization. A second approach that enables the variation of the design domain is the level-set method (Osher and Sethian 1988). As demonstrated in Luo et al. (2007), Wang et al. (2007), Zhang et al. (2012), Wang and Kang (2018), this method allows to elegantly combine shape and topology optimizations within a single framework. In comparison to density-based topology optimization methods, with the level-set method it is more difficult to mathematically describe the appearance of holes of arbitrary shape inside the design domain (Allaire and Jouve 2006). In addition, level-set functions may become too flat or too steep during the optimization procedure, such that they need to be periodically reinitialized in order to obtain a solution with sufficient numerical accuracy, which increases the computational cost. An approach somewhat similar to the level-set method is the phasefield method, which uses a partial differential equation for describing the evolution of topologies in a structural optimization problem (Takezawa et al. 2010;Gain and Paulino 2012). Although the topologies are solved over the complete design domain without the need of prior information on the location of internal boundaries, the phase-field method does not facilitate the generation of holes, which thus limits the range of possible topological solutions. The same limitation applies to the moving morphable components (MMC) method (Guo et al. 2014a;Zhang et al. 2016), in which the topology of a structure is optimized by adapting the layout of a set of morphable structural components. In contrast to the SIMP method, however, in the MMC method the description of the structural topology is independent of the finite element discretization, which has the advantage that the number of topology design variables reduces strongly, such that the computational efficiency improves significantly.
In order to avoid the jagged geometrical boundaries typically resulting from optimizing the topology of a discretized finite element model with the SIMP method, in Ramm (1995, 1997), Maute et al. (1998), Schwarz et al. (2001) the structural layout calculated with the SIMP method is smoothed using spline functions. Apart from improving the quality of the optimization result, this adaptation also increases the computational speed of the SIMP approach by reducing the number of active optimization variables. In Ansola et al. (2002) this so-called adaptive topology optimization technique has been subsequently extended towards a combined shape and topology approach for shell structures, whereby the design domain is defined by B-splines and their variability is accounted for by modifying the positions of the control points. The determination of the shape sensitivities is performed numerically using a finite difference method, which, although computationally demanding, has shown to lead to stable, converged results.
Along the lines of Ansola et al. (2002), in the present communication a coupled domain shape and topology optimization method is formulated by incorporating the shape design variables of the design domain in the SIMP topology optimization method. However, the formulation is generalized by considering design domains of arbitrary shape. The shape and topology design variables are incrementally updated in a sequential fashion, using a so-called staggered numerical update scheme. Non-Uniform Rational B-Splines (NURBS) are employed to parameterize the shape of the design domain (Piegl and Tiller 1995). This not only guarantees a highly accurate description of the shape boundaries by means of smooth basis functions with compact support, but also enables an efficient control of the design domain with only a few control points (Hsu 1994). Further, the optimization process is performed in a computationally efficient way by applying a gradient-based optimization algorithm, for which the sensitivities can be computed in closed form; this is another important difference compared to the approach advocated in Ansola et al. (2002).
The objective function adopted in the coupled optimization method refers to the structural compliance, which is the common objective used in topology optimization. However, if required, it is relatively straightforward to replace this objective function in the formulation by an alternative choice, or by a combination of various objective functions, in accordance with a multi-objective optimization approach. Minimizing the structural compliance under the constraint of a significant weight reduction by searching for optima regarding the shape as well as the topology certainly has significant practical value. As will be shown through various numerical examples, it may lead to substantially different structural designs compared to when this optimization is performed on the topology only, thus offering interesting, alternative structural solutions. The main features of the coupled optimization method are demonstrated by means of illustrative 2D benchmark problems; for the application of the method in advanced 3D problems the reader is referred to Wang et al. (2020aWang et al. ( , 2020b. The manuscript is organized as follows. Section 2 introduces the basic features of the coupled topology and domain shape optimization framework, by starting from the classical formulation for topology optimization. The minimization problem is formulated and the staggered numerical update scheme is described. Section 3 reviews the mathematical description of the geometry of the design domain by means of NURBS. In Sect. 4, the Finite Element Method (FEM) formulation for the coupled optimization formulation is summarized. The FEM formulation forms the basis for the computation of the sensitivities of the structural compliance to the shape and topology design variables, which is presented in Sect. 5 together with a description of the update of the design variables. Section 6 treats several numerical benchmark problems that demonstrate the basic features of the coupled optimization method. Finally, in Sect. 7 some concluding remarks are given.

Problem formulation and solution strategy
In accordance with the SIMP approach, topology optimization can be mathematically formulated by means of the following minimization problem (Sigmund 2001;Bendsøe and Sigmund 2003): where the objective function c typically is referred to as the structural compliance, is the vector of design variables that is composed of the relative densities e of the elements e, min is the prescribed minimum density (which is non-zero to avoid a singular stiffness matrix), N is the total number of elements, V( ) and V 0 represent the current (material) volume and the initial volume of the design domain, respectively,

3
Sequentially coupled gradient-based topology and domain… and f r is the volume fraction of material. Further, and are the global force and displacement vectors, respectively, which are related via with the global stiffness matrix. The formulation given by Eq.
(1) refers to a fixed outer shape of the design domain. To allow for a variable outer shape, Eq. (1) is extended as where a is the vector of shape-related design variables with the elements l m and u m representing the lower and upper bounds of the design variable a m , and M is the total number of shape-related design variables. Note that in Eq. (3) the material volume V is a function of both the shape-related design variables and the density variables , and that the initial volume of the design domain V 0 corresponds to the initial shape.
The formulation given by Eq.
(3) can be solved by updating the shape variables and density variables in parallel using a so-called monolithic scheme, or by performing the update of the shape and density variables in a sequential fashion using a staggered scheme. In this work a staggered update scheme is used, which at a specific incremental step keeps the shape variables momentarily fixed while solving for the density variables in accordance with the topology optimization formulation given by Eq. (1). Subsequently, the density variables are frozen, and the shape variables are optimized by solving the minimization problem: This sequence of optimization steps is repeated until both optimization formulations have reached their convergence criterion. Figure 1 shows the flowchart of the staggered solution scheme. The first step is to construct a geometry description of the design domain to be analyzed, which is done by decomposing the design domain into design elements and determining the coordinates of the control points of these elements. Subsequently, a finite element model is constructed by meshing the design elements. The relative densities of the finite elements are optimized by solving Eq.
(1), which is done by performing a structural analysis, a topology sensitivity analysis, and a density update procedure. Accordingly, the optimized relative density distribution is found for the design domain (for which the shape has been momentarily kept fixed). The next step is to optimize the shape of the design domain, based on the optimized densities of the finite elements just computed in the topology optimization step. This step is carried out by solving Eq. (4), in which a structural analysis, a shape sensitivity analysis, and a shape update procedure are performed. When the relative element densities computed in the last topology optimization step are not optimal for the new shape, a new topology optimization step needs to be conducted based on the new shape. The above sequence of steps is repeated until the stop criteria for both topology optimization and domain shape optimization are satisfied.
Note that the solution strategy in Fig. 1 starts with topology optimization and may therefore be referred to as Topology Shape Optimization (TSO). Alternatively, it can start with domain shape optimization, and will then be referred to as Shape Topology Optimization (STO). Although the element densities are kept fixed during an incremental shape optimization step, the overall density distribution, and thus the internal topology, changes as a result of a change in the nodal coordinates of the finite elements. However, the shape changes during an incremental shape optimization step are such that the corresponding alterations in the topology typically are relatively small and thus acceptable. Moreover, possible inaccuracies generated by these topology alterations may be expected to vanish in the subsequent incremental topology optimization step.

Geometry description
As visualized in Fig. 2, the description of the geometry of the design domain can be obtained from the parametric domain 0 using a geometry function , see also Vuong et al. (2010). Accordingly, for an arbitrary point ( , ) in the parametric domain 0 , the corresponding point (x, y) in the physical design domain is computed via the projection  Fig. 1 Flowchart of the staggered solution scheme for coupled topology and domain shape optimization with iterations g and h referring to topology and shape optimization subloops, respectively, and iteration k referring to the outer loop with the geometry function ( , ) defined by a set of basis functions and control points in the physical design domain. Observe from Fig. 2 that the geometry function projects the fixed auxiliary mesh defined in the parametric domain to an adequate finite element mesh for the actual design domain, see also Zhang et al. (2010). The geometry function is here constructed by NURBS, which have the advantage that they can describe a large variety of relatively complex geometries with high accuracy (Hughes et al. 2005). Since NURBS are based on B-splines, these will be reviewed first for reasons of clarity.

B-spline surface
For constructing a B-spline surface, the coordinates in the two-dimensional parametric domain are assembled in a non-decreasing fashion by using so-called knot vectors (Piegl and Tiller 1995) where i, j are the knot indices, s and t are the number of control points (and thus the number of basis functions used for constructing the B-spline surface) in the x-and y-directions, and p, q are integers defining the polynomial orders of the basis functions. The B-spline basis functions S i,p ( ) in the -direction can be obtained recursively as Parameterization of the physical design domain via a coordinate projection with a geometry function connected to the parametric domain 0 . The solid lines in the parametric domain 0 define a fixed auxiliary mesh, which is projected through to the actual finite element mesh that discretizes the design domain with the index i running from 1 to the number of control points s. The B-spline basis functions T j,q ( ) in the -direction can be determined in a similar fashion as shown in Eq. (7). The control points are defined by vectors including their x-and y-coordinates, which construct the geometry in the physical design domain via a two-dimensional control net i,j , with i = 1, 2, … s, j = 1, 2, … t . The B-spline surface follows from combining the control net with the tensor product of the basis functions S i,p ( ) and T j,q ( ) as Accordingly, the B-spline surface ( , ) may represent the geometry function ( , ) indicated in Fig. 2.

NURBS surface
From the B-spline basis functions S i,p ( ) and T j,q ( ) , a rational surface R i,j ( , ) can be constructed as (Hughes et al. 2005) where i,j are the weights of the control points. The rational surface can be subsequently combined with the control net i,j , to obtain the NURBS surface: Using matrix-vector notation, Eq. (9) can be compactly written as with and (7) for

3
Sequentially coupled gradient-based topology and domain… and With these expressions, the NURBS surface given by Eq. (10) becomes Since the basis functions S i,p ( ) and T j,q ( ) form a partition of unity, the NURBS surface, Eq. (15), reduces to the B-spline surface, Eq. (8), when all weights i,j are set equal to unity (Hughes et al. 2005).
The control points with their corresponding weights can be conveniently stored in the matrices X and Y as and where X i,j , Y i,j (with i = 1, 2, … , s;j = 1, 2, … , t ) are the control point coordinates. Hence, the projection, Eq. (5), of coordinates from the parametric domain 0 to the physical design domain can be reformulated as With Eq. (18), a surface in the physical design domain thus is completely determined once the knot vectors and control points with their weights are defined.

Structural analysis
In order to solve an optimization problem with the finite element method, the design domain needs to be meshed by finite elements. This is done by creating a fixed auxiliary mesh through defining the finite element node locations in the parametric domain, see also Zhang et al. (2010), and subsequently mapping these node locations to the actual physical design domain using Eq. (18), see Fig. 2. This enables the construction of an explicit relationship between the shape design variables (represented by the locations of the NURBS control points) and the FEM mesh, so that the FEM mesh is automatically adapted for changes of the structural boundary caused by shape optimization.
The main equations of the finite element method are reviewed below. These equations are subsequently used for deriving a gradient-based optimization framework, by determining the sensitivities of the objective function with respect to the shape and topology design variables. Hence, for an individual finite element, the element stiffness matrix e is given by where e is the strain-displacement matrix, e is the constitutive matrix, | e | is the determinant of the Jacobian matrix e , and e and e are the parametric coordinates of the finite element shape functions (not to be confused with the parametric coordinates of the geometry function presented in the previous section). The stiffness matrix is calculated by means of numerical integration using Gauss quadrature, i.e., where i * , j * are the indices of the integration points, t * equals the number of integration points along each direction, and w i * and w j * are the weighting factors of the integration points. In the present study the two-dimensional optimization problems analyzed are meshed by 4-node quadrilateral plane-stress elements, for which the element strain-displacement matrix e present in Eq. (20) is composed as in which

3
Sequentially coupled gradient-based topology and domain… where x and y represent the spatial directions of the coordinate system in the physical design domain, the index n reflects the node number, and N n represents the corresponding displacement shape function. The 4-node plane-stress elements used in this study are equipped with standard quadratic shape functions: Note that in Eq. (22) the spatial derivatives of the shape functions are defined with respect to the spatial coordinates x and y in the physical domain, while in Eq. (23) the shape functions are defined by the spatial coordinates e and e in the parametric domain. In order to relate the derivatives of the shape functions in the two domains, the vectors e n and e n are introduced as The relation between the vectors e n and e n follows as with the Jacobian J e expressed by The Jacobian can be elaborated further from the definition of the displacement field e within an element e: where x n and y n are the coordinates of the element nodes. Writing the displacement at an arbitrary point within an element e as the difference between its actual coordinate and its initial coordinate, u e x = x − x 0 and u e y = y − y 0 , together with Eq. (27) the Jacobian, Eq. (26), becomes by which e in Eq. (21) can be expressed in terms of e and e using Eqs. (22), (24), (25) and (28).
For a 2D plane-stress state, the constitutive matrix e within an element reads where is Poisson's ratio (which, for simplicity, is taken as uniform across the whole design domain), and E e is the Young's modulus of an element, which, in accordance with the SIMP approach for topology optimization, may be formulated as Sigmund 1999, 2003;Sigmund 2001;Andreassen et al. 2011) with e the relative density within an element, p a penalization factor on the bulk stiffness (a typical value being 3, see Sigmund (2001)), and E 0 the initial Young's modulus of the material. With Eqs. (21)-(30), Eq. (20) can be used to calculate the stiffness matrix e for each element: with 0 representing the initial element stiffness matrix, computed from Eq. (20), but with the constitutive matrix e replaced by Sequentially coupled gradient-based topology and domain… Finally, the global stiffness matrix can be constructed by assembling the contributions of the element stiffness matrices e : where N is the total number of elements, and the sum operator represents the typical element assembly procedure used in the finite element method. Denoting bc as the global stiffness matrix obtained after incorporating the boundary conditions, the nodal displacements result from inverting Eq. (2), i.e., This result can be used to compute the structural compliance as In addition, the volume of each element e follows from which, in accordance with Eq. (1), can be employed to calculate the total material volume required for topology optimization as well as the volume of the initial design domain

Sensitivity analysis and update of design variables
In order to establish a gradient-based optimization framework, the sensitivity of the objective (i.e., the structural compliance) with respect to the shape design variables (net of control points) and the topology optimization variables (element relative densities) will be derived in this section. Subsequently, the update procedure of the design variables is discussed.

Shape sensitivity analysis
As already explained, the meshing of the design domain is performed by projecting the nodal coordinates ( n , n ) of an auxiliary mesh in the parametric domain to the nodal coordinates (x n , y n ) of the actual finite elements in the physical domain. Using the NURBS surface definition given by Eq. (18), with the dependencies on the nodal coordinates of the auxiliary mesh indicated, this projection reads Inserting Eq. (39) into Eq. (28), the Jacobian e at element level can be expressed in terms of NURBS surface characteristics: in which the dependency of the basis functions and on, respectively, n and n have been omitted for reasons of brevity. In the domain shape optimization process, a limited number of control points is typically used for varying the shape. These control points are stored in the vector of design variables a. With the use of Eq. (40), the derivative of the Jacobian e with respect to a single design variable a m can be expressed as Furthermore, from Eq. (40) the determinant of the Jacobian matrix becomes

3
Sequentially coupled gradient-based topology and domain… With this expression, the partial derivative of the determinant of the Jacobian | e | with respect to design variable a m results from applying the chain rule: Combining Eqs. (21), (22) and (24)     in which e ∕ a m and | e |∕ a m are provided by Eqs. (45) and (43), respectively. Subsequently, from Eq. (31) the partial derivative of the element stiffness matrix e with respect to the design variable a m is calculated via with 0 ∕ a m thus given by Eq. (47). Equation (48) is substituted into the partial derivative of the total stiffness matrix with respect to the design variable a m , which via Eq. (33) becomes: After taking into account the boundary conditions, → bc , Eq. (49) is used to construct the partial derivative of the displacement with respect to design variable a m , which, with Eq. (34), leads to When assuming that the external loads are independent of the design variables a m , Eq. (50) reduces to In order to perform the domain shape optimization procedure given by Eq. (4), from Eq. (35) the partial derivative of the structural compliance c with respect to design variable a m needs to be calculated, which, after substituting Eq. (51), results in The partial derivative of the material volume V with respect to design variable a m can be formulated based on Eqs. (36) and (37), leading to whereby | e |∕ a m is given by Eq. (43).

Topology sensitivity analysis
To perform the topology optimization procedure presented in Eq.
(1), similar to Eq. (52) the partial derivative of the structural compliance c with respect to the relative density e is computed as The partial derivative of the element stiffness matrix e with respect to e can be determined via Eq. (31): Hence, the dependency of the total stiffness matrix bc on e enters through the individual element stiffness matrices e . Accordingly, Eq. (55) is developed for a specific element e using Eq. (56), i.e., In order to increase the robustness of the solution procedure and to alleviate possible checkerboard patterns in the spatial distribution of the element relative densities, Eq. (57) is slightly modified by adopting the sensitivity filter suggested in Sigmund (2001), i.e., in which Ĥ f = max(0, r min − dist(e, f )) , where dist(e, f ) is the distance between the center of element e and the center of element f, and r min is the radius adopted for the sensitivity filter. The filter radius is given by r min = r 0 √ V e , where r 0 is a predefined scaling parameter. Although the dependency of the filter radius r min on the element volume V e compromises the mesh-independent aspect of the sensitivity filter, it avoids the emergence of checkerboard patterns in regions where shape enlargement causes enlarged elements to lose their filtering capacity in the case the filter radius is kept fixed. Finally, from Eq. (37) the partial derivative of the material volume V with respect to the relative density e is calculated as

Update of design variables
The domain shape optimization procedure formulated by Eq. (4) is carried out by applying a Sequential Quadratic Programming (SQP) method, using the fmincon solver in MATLAB (Version R2016b), see Nocedal and Wright (2006) for more details. The element densities within the topology optimization procedure are updated using the so-called Optimality Criteria (OC) method described in Sigmund (2001). The topology optimization algorithm used in the present work is based on the formulation presented in Sigmund (2001). In this algorithm, however, all finite element volumes are taken as constant, which in the present approach is not the case, due to the combination with the shape optimization algorithm. Therefore, some minor adaptations were made to this topology optimization algorithm in order to account for the variability in element volume.

Cantilever beam problem
A popular benchmark problem for evaluating optimization methods is the cantilever beam problem (Wall et al. 2008;Qian 2010;Sigmund 2001;Andreassen et al. 2011). The coupled gradient-based topology and domain shape optimization approach will be applied to this benchmark problem, whereby the goal is to determine the geometry with a minimal structural compliance c, subject to appropriate geometrical constraints and a target material volume. The initial configuration of an elastic cantilever beam with Young's modulus E and Poisson's ratio is shown in Fig. 3. The left side of the beam is fully clamped (zero displacements along the left edge in all directions), and a vertical load F is applied at point Q at the right edge. The beam has a fixed length L and an initial

Fig. 3
Elastic cantilever beam with the control points i,j indicated height H 0 , whereby during domain shape optimization the actual height H of the beam may vary between H min and H max . The material volume V obtained after the optimization procedure should be 40% of that of the initial configuration, V 0 , in correspondence with a material volume fraction f r = 0.4 . The first problem that will be analyzed is referred to as Case 1, for which the values of the parameters introduced above are summarized in Table 1. Subsequently, two other problems will be studied, referred to as Case 2 and Case 3. Here, Case 2 is similar to Case 1, but the simulation now starts from the maximal beam height, H 0 = H max , while in Case 3 the beam height is left unbounded by not prescribing a maximal value H max . The geometry of the cantilever beam is described by NURBS. Setting the polynomial order of the basis functions, Eq. (7), as p=2 and q=1, the knot vectors for the B-spline surface, Eq. (6), can be defined as which correspond to a total of 6 control points. Table 2 lists the initial coordinates and weights of these control points. For simplicity, all weights have been chosen equal to unity, by which the NURBS surface reduces to a B-spline surface, see Eqs.
(8) and (10)-(14). The y-coordinates of the control points 1,2 , 2,2 , 3,2 govern the height, and thus the shape of the cantilever beam. As such, the vector with shape design variables a can be expressed as The auxiliary mesh in the initial parametric domain 0 is constructed by 200 × 25 equal-sized, rectangular parts, of which the corners via Eq. (5) are projected on the finite element nodes in the physical domain. Accordingly, the geometry in the physical domain is meshed with 5000 4-node plane-stress elements, equipped with quadratic shape functions, see Eq. (23). A mesh sensitivity study presented further in this communication will demonstrate that this FEM discretization is sufficiently fine for (60) = [0, 0, 0, 1, 1, 1] T , = [0, 0, 1, 1] T , (61) a = Y 1,2 , Y 2,2 , Y 3,2 T .  Table 2 Coordinates i,j (in m) and weights i,j of control points i j = 1 j = 2 1 1,1 = (0, 0) 1,1 = 1 1,2 = (0, 7) 1,2 = 1 2 2,1 = (20, 0) 2,1 = 1 2,2 = (20, 7) 2,2 = 1 3 3,1 = (30, 0) 3,1 = 1 3,2 = (30, 7) 3,2 = 1 obtaining accurate numerical results. For the domain shape optimization procedure, the settings of the MATLAB solver fmincon are as follows: the "Algorithm" is set to "sqp" (sequential quadratic programming), the "GradObj" and "GradConstr" are set to "on", which means that the sensitivities of the objective function and constraints used in the gradient-based domain shape optimization procedure are provided by the user. The remaining options are default. The domain shape optimization procedure is considered to be finished if the structural compliance meets the condition: | c h+1 − c h |≤1e-6, with h the iteration number, or if the shape design variables satisfy the criterion: max(| h+1 − h |) ≤1e-6. For the topology optimization, the scaling parameter r 0 characterizing the filter radius r min used in Eq. (58) is set equal to r 0 = 1.5 . Furthermore, the move limit appearing in the OC method (Sigmund 2001) is set to = 0.4 . The topology optimization procedure is finished if the structural compliance satisfies the requirement: | c g+1 − c g |≤1e-6, with g the iteration number. Finally, the stop criterion used in the outer loop in Fig. 1 where k is the iteration number of the outer loop. The above value of 1e-6 chosen for the stop criteria of the shape and topology optimization inner subloops was selected based on an separate variation study, in which for the STO and the TSO update schemes described in Sect. 2 the value of these stop criteria was systematically varied in the Case 1 and Case 3 problems introduced above. The variation study clearly demonstrated that the value found for the optimized compliance c increases monotonically with a growing value of the stop criteria. In specific, increasing the value of the stop criteria from 1e-6 to 1e-5 typically results in a percentual increase of the optimized compliance between 0.03% and 0.44%, while increasing the value from 1e-6 to 1e-3 increases the optimized compliance between 1.2% and 8.2%. Based on this result, it was concluded that a value of 1e-6 of the stop criteria for the shape and topology inner subloops results in solutions with sufficient numerical accuracy.
First the coupled topology and domain shape optimization with a TSO sequence is analyzed for Case 1. The optimization starts with a topology optimization step Sequentially coupled gradient-based topology and domain… of the initial design domain with a fixed shape, as shown in Fig. 4. Subsequently, a domain shape optimization step is performed whereby the finite element densities found at the end of the last topology optimization step are used as initial values, see Fig. 5. Note that the height of the clamped support at the left side of the cantilever beam is allowed to vary during the process of domain shape optimization. According to the solution strategy illustrated in Fig. 1, the alternating topology and shape optimization steps will continue until the stop criterion of the outer loop is satisfied. Figure 6 presents the value c of the objective function and the values of the shape design variables, Eq. (61), as a function of the outer loop iterations k, while Fig. 7 shows the shapes and topologies of the cantilever beam at specific iterations k of the outer loop. It can be observed that during the iterative process the shape changes from convex to concave, and that the topological layout of the interior "trusses" changes significantly as well. During the optimization procedure the shape clearly has a strong influence on the topology and vice versa; hence, the optimized structural layout is a trade-off between the density distribution and the shape. In order to investigate whether it makes a difference if the optimization procedure starts with a shape or a topology optimization step, Case 1 is also performed using an STO sequence. Figure 8 presents the initial, domain shape optimization step, whereby the relative densities are taken in accordance with a homogeneous distribution. Using the shape obtained at the end of the first step, a topology optimization step is carried out in turn, see Fig. 9. The results obtained after successive, combined shape and topology optimization steps are shown in Figs. 10 and 11, illustrating, respectively, the convergence behavior of the structural compliance c and the shape design variables a , and sketching the corresponding structural configurations at various iterations k of the outer loop. Comparing the results for the present STO sequence with those for the TSO sequence depicted in Figs. 6 and 7 shows that the TSO and STO sequences lead to a similar final configuration, but that the TSO sequence requires more iterations ( k = 14 ) to reach this configuration than the STO sequence ( k = 10 ). For both the TSO and STO sequence the number of the inner iterations g (topology) and h (domain shape) decrease with an increase of the    Figs. 6 and 10, respectively. The appearance of different local minima under different algorithmic schemes and search paths is well-known within the structural optimization community, see, for example, Sigmund and Petersson (1998). The number of local minima in structural optimization problems typically is large, and is characterized by the "landscape" created by the objective function(s), the design variables and constraints, the characteristics of the boundary value problem and the FEM discretization.
Since in the present coupled topology and domain shape optimization method the finite element mesh is chosen to morph together with design domain, it is necessary to analyze the effect of the FEM discretization on the optimization result in more detail. Accordingly, a comparison study has been performed considering four different initial FEM meshes for Case 1 (STO sequence), which are characterized by uniformly discretizing the parametric domain 0 into (a) 160 × 20 , (b) 180 × 20 , (c) 120 × 40 and (d) 200 × 25 elements. Note that the effect of the number of elements on the computational result is analyzed here by choosing meshes (a) and (b) to be 30-40% coarser than meshes (c) and (d). Figure 12 depicts the optimized configurations computed for the different initial meshes, together with the final FEM discretizations (whereby specific details can be distinguished by zooming in on the figures). Clearly, the various FEM meshes lead to similar layouts. In addition, for the relatively fine mesh (d) the final structural compliance c is only 1% higher than the value obtained for the other fine mesh (c), from which it is concluded that the difference in the initial width-to-height element aspect ratio for the two meshes (i.e., about 3:2 for mesh (c) versus 1:2 for mesh (d)) only marginally influences the final numerical result. Note, however, that the significant reduction of the structural height at the right beam end resulting from the domain shape optimization procedure causes the element width-to-height aspect ratio of mesh (c) locally to become nearly two times larger -and thus less optimized -than that of mesh (d). For the relatively coarse meshes (a) and (b) the final structural compliance respectively is only 7% and 5% higher than for the fine mesh (c). Additionally, the number of outer iterations k is the lowest for the relatively coarse mesh (a) (i.e., k = 7 ) and the highest for the relatively fine mesh (d) (i.e., k = 10 ). From the differences in the optimization results for the relatively coarse and fine meshes, and the minor sensitivity to the element aspect ratio used in the discretization of the two fine meshes, it may be concluded that the fine mesh (d) adopted in the forthcoming comparison study of Cases 1, 2 and 3 will provide adequate numerical results.
In the comparison study of the three cases, for Case 2 the height of the initial domain is set equal to the maximal value, H 0 = H max = 10 m. The volume (per unit depth) of the initial structure ( V 0 = 300 m 2 ) is larger than for Case 1 ( V 0 = 210 m 2 ), which is compensated for by scaling down the required volume fraction f r in Sequentially coupled gradient-based topology and domain… accordance with the ratio of the volumes, such that an adequate comparison between the two cases can be made, i.e., f r = (210∕300) × 0.4 = 0.28.
In order to clarify the mutual influence of topology and shape optimization for this case, the cantilever beam is first subjected to topology optimization only. By starting from a homogeneous density distribution, the converged density distribution obtained after g = 70 iterations is shown in Fig. 13. The structural compliance at convergence, c = 166.325 Nmm, turns out to be lower than for the STO optimization procedure of Case 1, c = 183.313 Nmm. In other words, the search path followed for Case 2 with a pure topology optimization procedure has led to a different local minimum for the structural compliance. The question now remains whether this result will alter when the topology optimization procedure for Case 2 is combined with domain shape optimization. The outcome of this analysis is shown in Fig. 14, illustrating that the topologies and the corresponding values of the structural compliance c computed with the TSO (Fig. 14a) and STO (Fig. 14b) sequences are virtually identical to those calculated with the pure topology optimization procedure shown in Fig. 13. Accordingly, it may be concluded that domain shape optimization has a negligible effect on the optimization of Case 2.  The computational results for Case 3, in which the beam height is left unbounded, are shown in Figs. 15 and 16 for the coupled optimization procedure with the TSO sequence, and in Figs. 17 and 18 for the coupled optimization procedure with the STO sequence. A comparison of the topologies of the two sequences depicted in Figs. 16 and 18 illustrates that these are quite similar, both leading to a triangular "truss" structure with no material present inside. The possible appearance of unfavorable element distortions in the triangular structure is avoided by locally using a relatively fine FEM mesh. As illustrated in Fig. 19, the FEM discretization of the structure as obtained after the optimization procedure of Case 3 indeed preserves the mesh quality at sharp corners. An alternative approach that could have been applied for preventing element distortions is based on the use of so-called boundary-cut elements, see Liu et al. (2016), Noël et al. (2017 for more details. Observe that the final value of the structural compliance c is different for the two update schemes, i.e., c = 159.647 Nmm (TSO) versus c = 149.705 Nmm (STO), and that the outer loop of the TSO sequence requires more iterations ( k = 11 ) to converge than the STO sequence ( k = 8 ). Furthermore, from Figs. 15 and 17 it can be seen that the final values of the shape variables for the two optimization sequences also differ somewhat, i.e., a = [13.0, 5.  Table 3 the minimal compliance values are summarized together with the computational times of Cases 1, 2 and 3 for the coupled optimization procedures with TSO and STO update schemes. The STO sequence appears to be 19-46% faster than the TSO sequence, and results in a minimal compliance value that is up to 6% lower. Hence, it may be concluded that the performance of the coupled optimization procedure with the STO update scheme here is superior to that with the TSO update scheme.

Curved L-shape beam problem
For the cantilever beam problem analyzed above the initial design domain has a basic, rectangular shape. In this section the features of the coupled topology and domain shape optimization approach are further explored for the curved L-shape beam configuration depicted in Fig. 20, which is characterized by a more complicated initial design domain. The L-shape beam configuration has proven to be a useful benchmark problem for the validation of topology optimization techniques, see, e.g., Zhang et al. (2012), Guo et al. (2014b), Duarte et al. (2015), Lian et al. (2017), Lieu and Lee (2017), Picelli et al. (2018).
As illustrated in Fig. 20, the top side of the L-shape beam is fully clamped, and a vertical load F = 10 kN is applied at the top corner of the free edge. The polynomial orders of the NURBS basis functions expressed by Eq. (7) are p=3 and q=1. Accordingly, the knot vectors for the B-spline surface, Eq. (6), are = [0, 0, 0, 0, 1, 1, 1, 1] T

Fig. 19
Finite element discretization of the final optimized design domain for Case 3 (STO sequence), with the mesh details at the upper left corner and right side of the beam illustrated in the insets  Fig. 20 are listed in Table 4. The shape design variables a are selected as a = X 1,1 , X 2,1 , Y 2,1 , X 3,1 , Y 3,1 T , with their upper and lower bounds given by For clarity, the FEM mesh of the curved L-shape beam is illustrated in Fig. 21, which has been constructed by 5000 plane-stress elements equipped with quadratic shape functions, see Eq. (23). Figure 22 illustrates the final configurations calculated with the coupled topology and domain shape optimization approach using the TSO sequence, Fig. 22b, and the STO sequence, Fig. 22c, together with the final configuration computed with pure topology optimization, Fig. 22d. Although the final configurations appear as similar, a comparison of the calculated values of the minimal compliance c illustrates that the coupled optimization approach performed with the STO sequence leads to an optimum c = 146.232 Nmm that lies 10% below the value found for the topology optimization procedure, c = 161.592 Nmm. This is, because for the coupled optimization approach the design space during the  Table 4 Coordinates i,j (in m) and weights i,j of control points i j = 1 j = 2 1 1,1 = (−12, 0) 1,1 = 1 1,2 = (0, 0) 1,2 = 1 2 2,1 = (−12, −28) 2,1 = 1 2,2 = (0, −18) 2,2 = 1 3 3,1 = (3, −30) 3,1 = 1 3,2 = (3, −20) 3,2 = 1 4 3,1 = (18, −19) 3,1 = 1 3,2 = (18, −18) 3,2 = 1 computational procedure is slightly altered as a result of the shape optimization process, which allows to find a lower optimum than for pure topology optimization in which the design space is kept fixed. In other words, the solutions found for coupled optimization with TSO and STO sequences are not incorporated in  Final configurations for the a curved L-shape beam, calculated with coupled topology and domain shape optimization using the b TSO sequence, and the c STO sequence, and computed with d pure Topology Optimization (TO) using the initial design domain 1 3 Sequentially coupled gradient-based topology and domain… the design space for pure topology optimization. In addition, the value determined from the coupled optimization procedure with the TSO sequence, c = 149.389 Nmm, lies 2% above the value c = 146.232 Nmm calculated with the STO sequence, indicating that the STO update scheme results in a more optimized design than the TSO update scheme. Conversely, as shown in Fig. 23, the number of iterations used in the outer loop of the TSO update scheme, k = 8 , is considerably lower than for the STO update scheme, k = 20 . In summary, from the above comparison study it is concluded that the coupled topology and domain shape optimization approach may find a better optimum than a pure topology optimization procedure when the shape of the initial design domain is specified, whereby the dimensions of this shape are variable within certain bounds.

Concluding remarks
In the current contribution a modeling framework is presented that couples topology and domain shape optimization. The shape and topology design variables are incrementally updated in a sequential fashion, using a so-called staggered numerical update scheme. Here, the shape of the design domain is momentarily frozen when performing a topology optimization step, and the topology of the design domain is temporarily kept fixed when carrying out a shape optimization step. The sequence of alternating shape and topology optimization steps is repeated until the value of the structural compliance c meets the convergence criteria. The geometry for domain shape optimization is described by NURBS, in order to allow for the analysis of relatively complex structural forms. The domain shape optimization procedure is performed using an SQP method in combination with an efficient gradient-based solution strategy. The topology optimization procedure is carried out in accordance with the SIMP approach. The usefulness of the coupled optimization approach has been demonstrated by analyzing several benchmark problems. The simulation results confirm a strong dependency on the initial conditions and bounds selected, showing that a careful construction of the design domain generally is necessary and meaningful. Under the specific conditions whereby the initial design domain has a basic, rectangular shape, of which the size is prescribed by a relatively large value, the effect of domain shape optimization on the computational result may be relatively small, due to which the coupled topology and domain shape optimization framework leads to a similar outcome as that obtained by pure topology optimization. Conversely, the coupled optimization approach may find a better optimum than the pure topology optimization approach when the shape of the initial design domain is specified, but the dimensions of the shape are allowed to vary substantially during the coupled optimization procedure. The latter approach is suitable for the optimization of aeroelastic structures, whereby the applied loads are determined by the actual size and shape of the design domain, see Wang et al. (2020b, c).
The focus of the present work is on the formulation of a coupled topology and domain shape optimization method, and the demonstration of its main features by means of illustrative numerical benchmark problems. Although the benchmark problems address several numerical discretization effects, such as the mesh sensitivity of the numerical result and the possible appearance of unfavorable element distortions, more research on such aspects is necessary in order to apply the coupled optimization method for advanced three-dimensional structures with complex geometries. Along this line, an interesting topic for future research may be to integrate the current optimization framework with an isogeometric, higher-order discretization method that uses an automated mesh generation procedure.