Sequentially coupled shape and topology optimization for 2.5D and 3D beam models

A sequentially coupled shape and topology optimization framework is presented in which the outer geometry and the internal topological layout of beam-type structures are optimized simultaneously. The outer geometry of the beam-type structures is parametrically described by non-uniform rational B-splines (NURBS), which guarantees a highly accurate description of the structural shape and enable an efficient control of the design domain with only a few control points. The computational efficiency of the coupled optimization approach is assured by applying a gradient-based optimization algorithm, for which the sensitivities are derived in closed form. The formulation of the coupled optimization approach is tailored toward 2.5D and full 3D representations of beam structures used in engineering applications. The 2.5D beam model, which has been taken from the literature, uses standard beam elements to simulate the beam response in the longitudinal direction, whereby the cross-sectional properties of the beam elements are calculated from additional 2D finite element method (FEM) analyses. A comparison study of a cantilever beam problem subjected to pure shape optimization and pure topology optimization illustrates that the 2.5D and 3D beam models lead to similar shape and topology designs, but that the 2.5D beam model has a significantly higher computational efficiency. Specifically, the computational times for the 2.5D model are about a factor 70 (shape optimization) and 1.4 (topology optimization) lower than for the 3D model, which indicates that in the coupled optimization approach the optimization of the shape provides the largest contribution to the higher computational efficiency of the 2.5D model. The coupled shape and topology optimization analysis subsequently performed on the 2.5D cantilever beam model demonstrates that the specific order at which the alternating shape and topology optimization increments are performed in the staggered update procedure turns out to have some influence on the computational speed and the value of the minimal compliance computed. Despite these differences, the final beam structures following from the different staggered update procedures illustrate how shape and topology can be efficiently optimized in an integrated, coupled fashion.


NURBS curve
A NURBS curve is a piecewise polynomial curve composed of B-spline basis functions. The definition of such a basis function is based on a set of non-decreasing real numbers in a parametric space, defined by a so-called knot vector Z [ζ 1 , ζ 2 , ..., ζ l+r +1 ] T with where each real number ζ k reflects a knot with k being the knot index, l is the number of basis functions that comprises the B-spline, and r is the polynomial order of the B-spline. With the knot index k running from 1 to l, the B-spline basis functions L k,r are defined recursively for an arbitrary value of ζ , i.e., for r 0 : L k,0 (ζ ) 1 if ζ k ≤ ζ < ζ k+1 , 0 otherwise, and can be conveniently stored in a vector as The so-called control points define the geometry of the structure, and their locations during the optimization procedure may change. Together with the corresponding weights, the control points are stored in a matrix as where C k ∈ {X k , Y k , Z k }, in which X k , Y k , Z k with k 1, 2, . . . , l, are the x-, y-and z-coordinates of the control points and ω k is the corresponding weight. Accordingly, the coordinates of an arbitrary point p [x, y, z] T of the NURBS curve are determined by with the weights of the control points assembled as Ω 1D [ω 1 , ω 2 , . . . , ω l ] T .

NURBS surface
In order to construct a NURBS surface in a 2D parametric domain, the coordinates characterizing the basis functions are stored in two knot vectors, i.e., where i, j are the knot indexes, n, m are the number of basis functions and p, q are the polynomial orders. The basis functions N i, p and M j,q in the ξ -and η-directions can be obtained in a similar fashion as shown in Eq.
(2). These basis functions are stored in the vectors N b and M b as For generality, the 2D net of control points defining a NURBS surface is obtained here from the 3D grid of control points describing a NURBS solid. This is done by selecting a single layer of control points from the grid of control points describing the NURBS solid. Accordingly, by considering the directions of the knot vectors Ξ , H and Z in the parametric domain to correspond to the x-, y-and z-directions in the physical domain, the net of control points becomes P 2D k,C ⎡ ⎢ ⎢ ⎣ ω 1,1,k C 1,1,k ω 1,2,k C 1,2,k · · · ω 1,m,k C 1,m,k ω 2,1,k C 2,1,k ω 2,2,k C 2,2,k · · · ω 2,m,k C 2,m,k . . . . . . . . . . . . ω n,1,k C n,1,k ω n,2,k C n,2,k · · · ω n,m,k C n,m,k with C i, j,k ∈ {X i, j,k , Y i, j,k , Z i, j,k }, where X i, j,k , Y i, j,k and Z i, j,k represent the x-, y-and z-coordinates of the control points, respectively. Here, the first two indices vary as i 1, 2, . . . , n, j 1, 2, . . . , m, while k is kept fixed as it refers to a specific layer of the 3D grid of control points. The resulting NURBS surface is defined by the tensor product of the 2D net of control points and the basis functions N b and M b . As such, an arbitrary point p of the NURBS surface can be calculated as with the weights of the control points as Ω 2D k ⎡ ⎢ ⎢ ⎣ ω 1,1,k ω 1,2,k · · · ω 1,m,k ω 2,1,k ω 2,2,k · · · ω 2,m,k . . . . . . . . . . . . ω n,1,k ω n,2,k · · · ω n,m,k where k ∈ {1, 2, . . . , l}, and P 2D k,C follows from Eq. (10).

NURBS solid
Although a NURBS solid can be defined in a similar fashion as the NURBS curve and surface, it is relatively straightforward to formulate the NURBS solid by applying projections of NURBS surfaces. Using the 2D control net, Eq. (10), and the NURBS surface basis functions N b and M b given by Eqs. (8) and (9), a matrix R 3D C can be constructed as where C ∈ {X, Y, Z } represents the x-, y-and z-coordinates of the control points. With Eq. (13), an arbitrary point p of the NURBS solid is described by with where Ω 2D k (with k 1, 2, . . . , l) is defined by Eq. (12), and N T b , M T b and L T b are the basis functions in the three directions of the NURBS solid.

Structural analysis model
Consider the 3D beam-type structure shown in Fig. 1a, which will be used for demonstrating specific features of the shape and topology optimization procedure. The beam is modeled by means of the FEM, whereby two types of element discretization are adopted, namely (i) a discretization by 3D continuum elements, see Fig. 1b, and (ii) a discretization by 1D beam elements in the longitudinal direction of the beam, combined with additional 2D continuum elements for a detailed modeling of the cross-sectional properties of a beam element, see Fig. 1c. For brevity, in the following these two types of discretization will be denoted as the "3D beam model" and the "2.5D beam model", respectively.
The finite element discretizations for the 3D and 2.5D beam models are obtained by projecting a fixed, auxiliary mesh, defined in the parametric domain, onto the actual mesh in the physical domain, using Eqs. (5), (11), and (14) for 1D, 2D and 3D finite elements. This projection is schematically illustrated in Fig. 2, and enables obtaining an explicit relationship between the structural response of the beam structure and the parameters that control the geometry, i.e., the NURBS control points. Accordingly, analytical expressions can be derived for the parameter sensitivities governing the structural shape optimization, which will be presented in Sect. 5.

3D beam model
Given the knot vectors, the control points with their weights, and a fixed auxiliary mesh in the parametric domain, the discretization of the 3D finite element model is carried out by projecting the nodal coordinates (ξ n , η n , ζ n ) of the fixed, auxiliary mesh onto the nodal coordinates (x n , y n , z n ) of the actual mesh in the physical domain, in correspondence with Eq. (14). Hence, for each isoparametric element e, its nodal coordinates can be cast into a vector:p where the coordinates p n of node n, with n 1, 2, . . . , n e , are obtained from Eq. (14), with n e the total number of nodes of element e. The location of an arbitrary pointp e [x, y, z] T within the finite element e can be obtained from an interpolation of the coordinates of the element nodesp e as: with where N n , with n 1, 2, . . . , n e , are the polynomial interpolation functions of the element, and I d represents the d × d identity matrix, with d the number of degrees of freedom of the element node. By introducing the gradient operator ∂ [∂/∂ξ e , ∂/∂η e , ∂/∂ζ e ] T , the Jacobian matrix of element e follows as Inserting Eq. (17) into Eq. (19), the Jacobian matrix J e can be directly expressed in terms of the nodal coordinatesp e , i.e., wherep e is given by Eq. (16). When denoting the displacement vector for an arbitrary point of the structure as s [s x , s y , s z ] T , the corresponding strain vector [ x x , yy , 2 xy , 2 xz , 2 yz , zz ] T can be calculated as Bs, with the matrix B given by: At the element level, the matrix B e can be expressed as [5] with where the index c refers to the column of the inverse of the Jacobian matrix, Eq. (20). Using Eqs. (16) to (24), the stiffness matrix k e of the isoparametric 3D element is computed as with |J e | the determinant of the Jacobian matrix J e . Here, the 6 × 6 matrix Q e incorporates the constitutive properties of the element. In accordance with the SIMP approach used for topology optimization, the stiffness matrix is scaled by the relative density of the element via a power-law expression [2,18,23] Q e (ρ e ) p Q, (26) where ρ e is the relative density of element e, p is a penalization factor with a typical value of 3 [2] that is also adopted in the present work, and Q represents the constitutive behavior in a material point-taken here as linear-elastic-in accordance with σ Q , with the stress given by σ [σ x x , σ yy , σ xy , σ xz , σ yz , σ zz ] T . Finally, the global stiffness matrix K is found by assembling the element stiffness matrices: with N the total number of elements and the sum operator representing the assembly procedure typically used in the finite element method. The global displacement vector u is obtained by solving the system of equilibrium equations where f is the global force vector.

2.5D beam model
For the 2.5D beam model, standard 1D beam elements are used to discretize the beam in the longitudinal direction, while the cross-sectional properties of the beam elements are calculated from additional 2D FEM analyses. In this section, the 2.5D FEM approach originally presented in [12] is reviewed, see also [3][4][5][6]11] for more details; these expressions serve as input for the computation of the shape and topology sensitivities used in the gradient-based optimization algorithm, see Sect. 5.

Beam finite element model
In the beam FEM model, the deformations are described by translations χ [χ x , χ y , χ z ] T and rotations ϕ [ϕ x , ϕ y , ϕ z ] T of a cross-sectional reference point. These deformations are grouped into a single vector as r [χ T , ϕ T ] T . In addition, the shear strains in the x-and y-directions and the normal strain in the z-direction are assembled as τ [τ x , τ y , τ z ] T , and the curvatures about the three directions are given by κ [κ x , κ y , κ z ] T . The above two deformation measures can be assembled into a single vector as ψ [τ T , κ T ] T . Accordingly, the strain-displacement relation becomes ψ B r, with the operatorB given bȳ with Further, I 6 is the 6 × 6 identity matrix, 0 3 is the 3 × 3 matrix null matrix, and the z-direction corresponds to the longitudinal direction of the beam, see Fig. 1.
From the strains and curvatures, the normal force and shear forces T [T x , T y , T z ] T and bending moment and torsional moments M [M x , M y , M z ] T can be determined as follows: with θ [T T , M T ] T and K s representing the stiffness matrix of the beam cross section.
Analogous to Eq.(28) for a 3D beam model, the equilibrium equations for a 2.5D beam model are defined byKū where the stiffness k b e , of the 1D beam element is computed as Here,K is the global stiffness matrix,ū contains the nodal displacements and rotations of the beam, N e includes the interpolation functions for the beam displacements and rotations, the operator B e for each element e is given by Eq. (29), and |J e | represents the determinant of the Jacobian matrix, withJ e expressed as where ∂z/∂ζ e is calculated from Eq. (17) using the beam interpolation functions. The derivation of the crosssectional stiffness matrix K s e of the beam element, which appears in Eq. (34), is discussed in Sect. 3.2.2 below.

Cross-sectional stiffness
The displacements s of an arbitrary point of the beam cross section can be expressed in terms of the rigidbody displacements and rotations r of the cross-sectional reference point, complemented by the in-plane and out-of-plane warping displacements s w associated with the deformation of the cross section, i.e., where x, y, and z are the coordinates of the specific location in the cross section, with x and y defined with respect to the location of the cross-sectional reference point. Using a FEM model for the cross section, the warping displacements s w are discretized as where N contains the interpolation functions of the cross-sectional finite element andû represents the vector with the nodal warping displacements of the element. As will be shown below, it is convenient to isolate the partial derivative terms ∂/∂z in the strain-displacement relation, which is done by combining Bs with Eq. (22), resulting in which uses the relation [5]B with T r given by Eq. (30). Invoking the kinematic relation ψ B r withB given by Eqs. (29), (42) finally turns into As a next step, the equilibrium equations for the beam cross section need to be formulated. In accordance with the variational framework presented in [3][4][5][6], equilibrium is described in terms ofû, ψ and θ by the following set of partial differential equations: Note that the generalized loads θ given by Eq. (31) are related to the stress components by θ A Z T σ b d A, with Z given by Eq. (37) and the stress σ b presented by σ b [σ xz , σ yz , σ zz ] T . Further, the matrices A, R, E, C, L and M presented by Eq. (46) are determined from the FEM model at cross-sectional level. In Eq. (46), the matrix Z e is obtained from Eq. (37), whereby the x-and y-coordinates within a cross-sectional element are calculated from combining the element interpolation functions with the nodal coordinates of the element. The nodal coordinates of an arbitrary cross-sectional finite element e are computed via Eq. (16), where p n represents the nodal coordinates determined from Eq. (11), with the z-coordinate of the control points Z i, j,k in Eq. (10) set to zero. Furthermore, the matrix S e is computed from Eq. (41), and Q e follows from Eq. (26), in which ρ e refers to the relative material density in finite element e, and |Ĵ e | is the determinant of the Jacobian matrix of the element. The Jacobian matrix has the usual form which, by invoking Eqs. (10), (11), (16), (17), (18) and (19), can be expressed as a function of the element nodal coordinates (and thus as function of the NURBS control points): with the derivative operator∂ [∂/∂ξ e , ∂/∂η e , 0] T and the matrixP e given by Eq. (21). After computing the inverse of the Jacobian matrix Eq. (48), the matrixB e appearing in Eq. (46) can be developed from Eq. (40) aŝ whereB In order to eliminate the rigid body motions and to ensure that the cross-sectional warping displacements satisfy the orthogonality conditions formulated in Ladevèze and Simmonds [16], El Fatmi [8] and Genoese et al. [10], the following constraint equations need to be satisfied: with where n n is the total number of element nodes in the FEM mesh and Z n is obtained by inserting the nodal coordinates of node n in Eq. (37). Formulating Eq. (45) in matrix-vector notation and incorporating the constraint relations, Eq. (51), results in the system of equilibrium equations [6] Kw f ⇐⇒ with where For computational convenience, a 6 × 6 matrix W is constructed by column-wise assembling the solutions of Eq. (53) for 6 different right-hand sides, which are distinguished by successively setting one specific entry of θ to unity and the remaining 5 entries to zero. Accordingly, for an arbitrary force vector θ , with the definition of θ given below Eq. (31), the solution vector w [w 1 , w 2 ] T can be straightforwardly computed as In addition, by equating the external complementary virtual work to the internal complementary virtual work, the compliance matrix F s of the cross section follows as [6] Fig. 3 Flowchart of the incremental-iterative solution strategy for coupled shape and topology optimization, with iterations g and h referring to topology and shape optimization sub-loops, respectively, and iteration w referring to the outer loop with where the matrices A, R, E, M are provided by Eq. (46), and the matrix H is given by Eq. (54). Finally, the stiffness matrix of the element cross section used in Eq. (31) can be calculated via K s F −1 s , with F s in accordance with Eq. (56).

Coupled optimization model
This section presents a coupled optimization model that is able to simultaneously optimize the outer shape and the topological layout of a beam structure. For this purpose, the formulation recently proposed in [30] is followed, which incorporates the shape-related design variables in the SIMP topology optimization method. Accordingly, the minimization of the compliance of a structure under a predefined material volume constraint is formulated as in which the structural compliance c depends on both the shape design variables a (i.e., the spatial control points of the NURBS) and the topology design variables ρ (i.e., the relative densities). Further, the variables l s and u s are the lower and upper bounds of the shape design variable a s , S is the total number of shape design variables, ρ min is the minimum density, e is the element number and N is the total number of elements. The material volume V is a function of both a and ρ, and the parameter V 0 represents the initial volume of the design domain. The variable f r is the prescribed volume fraction, and f is the global force vector. The displacement vector u of the modeled system is obtained by solving the corresponding equilibrium equations, i.e., Eq. (28) for a 3D beam model and Eq. (32) for a 2.5D beam model. The coupled shape and topology optimization problem for the beam structure is solved in an incrementaliterative fashion using a staggered (sequential) solution strategy, which is illustrated in Fig. 3: 1. The geometry of the design domain is parametrically described by NURBS, as introduced in Sect. 2.
Subsequently, for performing a structural analysis, the geometry data are discretized into a FEM model, as described in Sect. 3.
2. By keeping the shape design variables a momentarily fixed, the problem formulated in Eq. (58) can be regarded as a classical topology optimization problem, formulated as Accordingly, the relative densities ρ are optimized by solving Eq. (59), which is carried out by following the typical topology optimization procedure, consisting of a structural analysis, a topology sensitivity analysis, and the update of the densities. 3. By momentarily freezing the updated densities ρ, Eq. (58) reduces to a shape optimization problem, formulated as The shape design variables a are calculated by solving Eq. (60) based on the densities ρ found in step 2.
Step 3 is carried out by performing a standard gradient-based shape optimization procedure, consisting of a structural analysis, a shape sensitivity analysis, and the updating of the shape. 4. When the element densities ρ obtained from step 2 are not optimal for the shape subsequently computed in step 3, in a sense that they do not lead to an immediate satisfaction of the shape convergence criteria (i.e., after the first iteration), a new topology optimization increment (step 2) is performed, based on the updated shape design variables a. After this step has converged and the topology design variables ρ have been updated, a new shape optimization increment (step 3) is performed, and the above check on immediate convergence is repeated. The above process is continued until the incremental-iterative solution procedure has converged. Note that updating the topology and shape configurations of the beam after convergence of the corresponding optimization increment causes the coupling effects between topology and shape to be automatically accounted for. The convergence (or stop) criteria for the individual topology and shape optimization increments and for the overall procedure are, respectively, evaluated during every iteration g, h and w, see Fig. 3, and will be specified when discussing the numerical examples in Sect. 6.
The staggered solution strategy depicted in Fig. 3 starts with an incremental topology optimization step, and, as such, is referred to as coupled topology shape optimization (CTSO). Alternatively, it can start with an incremental shape optimization step, and will then be referred to as coupled shape topology optimization (CSTO), see also [30]. In shape and topology optimization problems, the landscape of solutions is characterized by numerous local minima, whereby the convergence speed and the specific optimized solution calculated may be expected to be sensitive to the algorithmic features of the numerical update scheme applied. Despite that the main focus of the present work is on numerical discretization aspects and on comparing the outer and internal geometries of the 3D and 2.5D beam models under shape and topology optimization, possible differences in the results from the CTSO and CSTO update schemes will be also analyzed and discussed. Note that the choice of a staggered update scheme in this study is arbitrary; the coupled optimization problem, Eq. (58), could also be solved concurrently by using a monolithic scheme. A systematic comparison of the influence of staggered and monolithic update schemes on the computational efficiency and outcome of the coupled optimization approach, however, falls outside the scope of this work; for more details on this aspect, the reader is referred to [29].

Sensitivity analysis
For the coupled, gradient-based optimization approach, the analytical sensitivities of the structural compliance c need to be determined with respect to the shape design variables a and relative element densities ρ. These analytical sensitivities are elaborated below for the 2.5D and 3D beam models.

2.5D beam model
For simplicity of notation, the variable d is introduced, which represents either the shape design variable a s or the element density ρ e . In correspondence with Eq. (56), the derivative of the cross-sectional compliance matrix F s with respect to d can be expressed as Since the matrix W contains the solutions of the equilibrium equations given by Eq. (53), whereby inf one component of θ equals unity and the remaining components are zero, the derivative ∂W/∂d in Eq. (61) follows from Eq. (53) as Inserting Eq. (62) into Eq. (61), the derivative of the cross-sectional compliance matrix becomes where the matrixK is obtained from Eqs. (53) and (54), and G follows from Eqs. (56) and (57). Since the stiffness matrix K s of the cross section is the inverse of the compliance matrix F s , the derivative of the stiffness matrix with respect to d may be expressed as The global stiffness matrixK of the beam model is based on the cross-sectional stiffness matrix K s e of the beam element via Eqs. (33) and (34). Accordingly, the derivative ofK with respect to the design variable d becomes Because the loads f are independent of the design variable d, from Eqs. (58), (32) and the symmetry relation K K T , the sensitivity of the compliance c with respect to d follows as The specific forms of the derivatives ∂K/∂d and ∂G/∂d appearing in the above equations depend on the type of design variable. The derivation of these terms is presented below.

Shape sensitivity analysis
The shape sensitivity analysis starts with the computation of the partial derivative of the Jacobian matrix, Eq.
(48), with respect to the shape design variable a s , leading to in which, for reasons of brevity, the notation ∂ s (·) ∂(·)/∂a s has been used. The derivative ∂ s (P e ) appearing in Eq. (68) is calculated by combining Eqs. (21) and (16), in which the derivative of the coordinates of an arbitrary point on the NURBS surface with respect to the shape design variable follows from Eq. (11) as where ∂ s (P 2D k,C ), with C ∈ {X, Y, Z }, can be easily determined from Eq. (10), in which the z-coordinate of the control points Z i, j,k is set to zero. Further, the matrices N b , M b and Ω 2D k are obtained from Eqs. (8), (9) and (12), respectively. Eq. (68) can be used for the calculation of the derivative of the determinant of the Jacobian with respect to the shape design variable a s via Further, Eq. (68) can be employed in the expression for the derivative of the inverse of the Jacobian matrix with respect to the shape design variable a s , which reads: whereB c e , with c ∈ {1, 2}, is obtained from Eq. (50). Accordingly, the derivative ofB e N e with respect to the shape design variable a s can be formulated as with From the last two expressions in Eq. (75), the derivatives with respect to the shape design variable of the matrices A, R and L in Eq. (46) are calculated as with D given by Eq. (52). Finally, the terms ∂K/∂a s and ∂G/∂a s in Eq. (63) are calculated from Eqs. (53)-(57), thereby inserting the derivatives of the coefficient matrices given by Eqs. (75) and (76).

Topology sensitivity analysis
According to Eqs. (26) and (46), the coefficient matrices A, R, E, C, L and M of the cross section can be explicitly expressed as a function of the element density ρ e . Correspondingly, by referring to an arbitrary coefficient matrix as M c , its derivative with respect to the element density is calculated as where the variable M c 0 refers to the corresponding coefficient matrices obtained from Eq. (46) by substituting the density-independent constitutive matrix Q (instead of the density-dependent constitutive matrix Q e , see Eq. (26)). With Eq. (77), the terms ∂K/∂ρ e and ∂G/∂ρ e in Eq. (63) can be determined from Eqs. (53) to (57). To alleviate mesh dependency and possible checkerboard patterns in the solution computed by topology optimization, a sensitivity filter [23] is introduced that smoothens the sensitivity of the structural compliance as follows: where ∂c/∂ρ f is calculated employing Eqs. (63)-(67). Further,Ĥ f max(0, r min − dist(e, f )), in which dist(e, f ) is the distance between the center of element e and the center of element f, and r min is the radius of the circle within which smoothing takes place.

Shape sensitivity analysis
In the three-dimensional FEM model for the beam structure, the Jacobian matrix of an element e is presented by Eq. (20). Correspondingly, its derivative with respect to the shape design variable a s follows from The term ∂ s (P e ) in Eq. (79) is calculated by combining Eqs. (21) and (16) with the derivative of the nodal coordinates, Eq. (14), with respect to a s : with where ∂ s (P 2D k,C ), with C ∈ {X, Y, Z } and k 1, 2, . . . , l, is determined based on Eq. (10), and the matrices L b , N b , M b and Ω 3D k are obtained from Eqs. (3), (8), (9) and (15) while the derivative of the inverse of the Jacobian matrix becomes with the Jacobian matrix J e given by Eq. (20). In addition, the derivative of the matrix product B e N e used in Eq. (25) is calculated with the use of Eqs. (23) and (24) as with Following the assumption that the load vector f is independent of the design variables, the derivative of the structural compliance with respect to the shape design variable a s becomes where ∂ s (K) is given by Eq. (86) and u is the global displacement vector obtained from solving the system of equilibrium equations, Eq. (28).

Topology sensitivity analysis
In correspondence with Eqs. (28) and (58), the derivative of the structural compliance c with respect to the relative density ρ e equals ∂c ∂ρ e −u T ∂K ∂ρ e u.
Since the three-dimensional element stiffness matrix is an explicit function of the element density ρ e , see Eqs. (25) and (26), the derivative of the element stiffness matrix k e with respect to ρ e follows as ∂k e ∂ρ e p(ρ e ) p−1 k 0 , Fig. 4 Cantilever beam with a uniform rectangular cross section Table 1 Parameter values of the cantilever beam problem where k 0 is obtained from Eq. (25), by substituting the density-independent constitutive matrix Q (instead of the density-dependent constitutive matrix Q e , see Eq. (26)). Correspondingly, the derivative ∂c/∂ρ e provided by Eq. (88) is reformulated for a specific element e using Eqs. (89) and (27), i.e., The filtering of Eq. (90) is performed in accordance with Eq. (78), where r min is the radius of a sphere (instead of a circle) within which smoothing occurs.

Validation case studies
In order to validate the accuracy and efficiency of the coupled shape and topology optimization procedure for 2.5D beam models, an elastic cantilever beam is considered whereby the solution obtained for the 2.5D beam model is compared to that for the 3D beam model. The initial configuration of the cantilever beam is illustrated in Fig. 4, showing a slender beam with length L 0 and a uniform rectangular cross section with dimensions W 0 × H 0 . The left end of the beam is fully clamped (i.e., zero displacements and rotations in all directions), and the right end is subjected to a vertical load F. The values of the model parameters are listed in Table 1, whereby E and ν are the Young's modulus and Poisson's ratio of the beam material, respectively.

Shape optimization
As a first step, the comparison study focuses only on the optimization of the shape of the beam sketched in Fig. 4, and the topology optimization is temporarily left out of consideration. In the 3D model, the shape of the cantilever beam is parametrically described by a NURBS solid. As explained in Sect. 2.3, the NURBS solid is constructed using a one-dimensional NURBS curve in the axial direction of the beam and a two-dimensional NURBS surface for the description of the beam cross section. According to the definition of a NURBS surface introduced in Sect. 2.2, the polynomial orders in Eq.
Correspondingly, the NURBS surface is characterized by a net of 9 control points, with the arrangement of the control points at the kth location along the axial (z-)direction of the beam illustrated in Fig. 5, where C i, j,k (X i, j,k , Y i, j,k , Z i, j,k ), with i 1, 2, 3, j 1, 2, 3 and k 1, 2, . . . , l, see also Eq. (10). The initial coordinates and weights of these control points are listed in Table 2, in which Z k with k 1, 2, . . . , l is used to replace Z i, j,k since all control points of a specific surface have the same z-coordinate. The value of y x C 3,3,k C 2,3,k C 1,3,k C 1,2,k C 2,2,k C 3,2,k C 1,1,k C 2,1,k C 3,1,k  Table 2 Coordinates (in m) and corresponding weights of the control points of the kth (with k 1, 2, . . . , l) NURBS surface the z-coordinate Z k of each NURBS surface is determined based on the NURBS curve, see Sect. 2.1, which describes the geometry along the beam length. The polynomial order of the NURBS curve is r 1, which, in correspondence with Eq. (1), represents a linear curve. In order to construct the complete beam geometry, l 81 NURBS surfaces are equally spaced along the beam length L 0 . Accordingly, the knot vector Z in Eq.
(1) is defined as and the z-coordinate Z k of the control points is with the beam length L 0 given in Table 1.
With the above characterization of the beam geometry, the nodal positions in parametric ξ -, ηand ζcoordinates are obtained after dividing the parametric domain W 0 × H 0 × L 0 into equal 6 × 10 × 80 hexahedral parts. Subsequently, the nodal positions of the FEM model are found by projecting the auxiliary mesh in the parametric domain to the physical domain using Eq. (14), see also Fig. 2. These nodal positions are finally utilized to construct the 3D and 2.5D models of the cantilever beam, as illustrated in Fig. 1b and c, respectively. In accordance with the shape optimization procedure formulated in Eq. (60), the objective-the structural compliance c-and constraints-the volume V -are evaluated using the relations for the 3D FEM model introduced in Sect. 3.1. The relative element densities ρ are uniformly kept fixed at a value of 1. In order to limit the computational demand of the 3D FEM model in the comparison study on shape optimization of 3D and 2.5D beam models, only the height H of the beam is allowed to change along the beam length, whereby the cross section is prescribed to remain rectangular. Note that this restriction in the shape design variables does not affect the nature and objectivity of the comparison study; in fact, it could have been introduced equally well as a geometrical constraint following from structural design considerations. It is further emphasized that it is straightforward to incorporate the width of beam cross sections into the set of shape design variables, as is done for the coupled shape and topology optimization analysis of the 2.5D beam model considered in Sect. 6.2. In correspondence with the above arguments, the shape design variables a are formulated as where Y 1,3,k is the y-coordinate of control point C 1,3,k illustrated in Fig. 5. The upper and lower bounds of this y-coordinate are set as 0.5 m ≤ Y 1,3,k ≤ 1.5m, which also apply to the y-coordinate of neighboring control points via the equality Y 1,3,k Y 2,3,k Y 3,3,k . Further, the coordinates of the control points of the first NURBS surface (k 1) are taken the same as those of the adjacent, second NURBS surface (k 2); in fact, the first cross section is only used for post-processing of the analysis results, and does not contribute to the optimization procedure, see also [6]. The values of the shape design variables a are efficiently updated by incorporating the analytical shape sensitivities presented in Sect. 5.2.1 in a sequential quadratic programming (SQP) procedure that was implemented in the MATLAB solver fmincon. The shape optimization is assumed to be converged when the structural compliance meets the stop criterion | c h+1 − c h |≤1e-5, with h the iteration number, or when the shape design variables satisfy the criterion max(| a h+1 − a h |) ≤1e-5, whereby | . | refers to the absolute value of a scalar or the absolute value of each specific component of a vector.
Alternatively, the optimization procedure can be carried out by modeling the cantilever beam as a 2.5D model. The 2.5D model is composed of 80 two-node beam elements along the beam (z-)axis (shown as a dotted black line in Fig. 1c). The z-coordinate of the beam nodes is given by Eq. (93). The stiffness of each beam finite element is obtained from a 2D FEM analysis of the beam cross section (evaluated at one of the beam element nodes). The meshes of the 81 cross sections have a similar mesh density as in the 3D FEM model. In Table 2, the z-coordinate Z k (with k 1, 2, . . . , 81) of the control points depicted in Fig. 5 essentially is irrelevant as a result of the two-dimensional character of the NURBS surface, and therefore is arbitrarily set to zero. The shape design variables a relate to the height of the cross sections and are presented by Eq. (94). The shape optimization procedure is performed in accordance with Eq. (60), for which the analytical shape sensitivities are provided in Sect. 5.1.1.
The results of the shape optimization procedures with the 3D and 2.5D beam models are shown in Figs. 6 and 7. It can be observed that both models provide almost the same optimized shape, whereby the beam height near the clamped end, at which the bending moment is relatively large, is maximal and equal to 1.5m. Along the middle part of the beam the height gradually decreases, whereby it reaches the minimum value of 0.5 m close to the end at which the external load is applied. The convergence behavior of the shape optimization procedure for the two beam-type structures is illustrated in Fig. 8. The convergence rate of the structural compliance c for the 2.5D model is somewhat lower than for the 3D model, although the two models eventually provide the same minimum value of c. The computational time of the 2.5D model is about a factor of 70 lower than that of the 3D model, which illustrates that the 2.5D model is computationally much more efficient in this case.

Topology optimization
As a next step, the cantilever beam problem illustrated in Fig. 4 is subjected to pure topology optimization, whereby the results computed with the 2.5D and 3D beam models are mutually compared. The geometry and FEM models of the problem are in accordance with the description in Sect. 6.1.1. In the topology optimization procedure, the volume fraction f r appearing in Eq. (59) is set equal to 0.4. Further, the interior-point (IP) algorithm implemented in the MATLAB solver fmincon [20] is employed to update the element densities ρ. The sensitivities required for the gradient-based topology optimization procedure are presented in Sects. 5.1.2 and 5.2.2 for the 2.5D and 3D beam models, respectively. The topology optimization procedure is assumed to be converged when the structural compliance meets the stop criterion | c g+1 − c g |≤1e-5, with g the iteration number, or when the topology design variables satisfy the criterion max(| ρ g+1 − ρ g |) ≤1e-5. In order to warrant an objective comparison of the results calculated for the 2.5D and 3D beam models, the filtering of possible checkerboard patterns, which is defined by Eq. (78) and may work out differently for the 2.5D and 3D beam models, is omitted here. Figure 9 shows the convergence behavior of the structural compliance c of the topology optimization process for the 2.5D (circles) and 3D (triangles) beam models. The overall convergence behavior and the corresponding compliance values clearly are similar for the two beam models. As a result of three-dimensional stress effects, some minor differences emerge in the relative density distributions at the clamped end and at the point of load application, see Figs. 10 and 11. The density distributions across the cross sections illustrate that most of the material is located at the top and bottom of the beam height, which is logical for a beam structure loaded in bending. The optimized configurations presented in Figs. 10 and 11 further show the presence of intermediate element densities. Note that in topology optimization analyses the appearance of intermediate densities commonly is considered as undesirable, as it fades the distinction between a solid material (material density is unity) and a void (material density is zero). An efficient way to reduce intermediate densities in topology optimization results is by increasing the penalization factor p in Eq. (26), see [25], which in the present work has been selected in accordance with the common value p 3 [2]. Alternative, more elaborate strategies that mitigate the appearance of intermediate densities have been presented in [13,24,27]. It should be further mentioned that the canonical duality method recently proposed by [9] can produce O(1) optimal Structural compliance c (N·m) g Fig. 9 Convergence behavior of the compliance c for pure topology optimization of 3D (triangles) and 2.5D (circles) beam models. The iteration number g refers to the loop for topology optimization, see Fig. 3 density distributions without the use of a filtering step. The main purpose of Figs. 10 and 11, however, is to objectively compare the optimized configurations computed for 3D and 2.5D beam models, for which the appearance of intermediate densities is insignificant. In addition, in the analyses performed with the coupled shape and topology optimization approach as presented in Sect. 6.2, the appearance of intermediate densities turns out to be negligible. The computational time required for solving the 3D beam model is 1.4 times larger than the time needed for solving the 2.5D beam model, indicating that the reduction in spatial dimension here only provides a relatively small gain in computational efficiency as compared to the reduction of a factor of 70 found in Sect. 6.1.1 for pure shape optimization. This can be explained by considering the computational time of the topology optimization procedure in more detail. In particular, the time required for the structural analysis of the beam for the 2.5D model is one order of magnitude smaller than for the 3D model. However, this gain in time is counterbalanced by the fact that the time required for the sensitivity analysis of the 2.5D model is one order of magnitude larger than that of the 3D model. The latter result is caused by two reasons. Firstly, in 3D topology optimization, the sensitivity of the global stiffness matrix can be efficiently obtained from the element stiffness matrices stored during the structural analysis, see Eqs. (88) and (89). In contrast, in 2.5D topology optimization, the calculation of the sensitivity of the global stiffness matrix requires the performance of a 1D FEM beam analysis, whereby the sensitivity of the stiffness matrix of each cross section needs to be determined a priori, see Eqs. (61)-(66) and (77), which obviously is computationally demanding. Secondly, the 1D beam FEM analysis must be repeated for the calculation of the sensitivity with respect to each element density. These two aspects thus make the topology sensitivity analysis for the 2.5D model computationally considerable more expensive than for the 3D model, and therefore result in only a relatively small decrease of the overall computational time by a factor of 1.4.
It should be noted that in the 2.5D beam model the number of cross sections has been chosen in accordance with the element density of the 3D beam model in longitudinal direction, in order to perform a consistent comparison of the simulation results at the cross-sectional level, see Fig. 10. From the viewpoint of computational efficiency, however, the number of cross sections in the 2.5D beam model could have been taken substantially lower without significantly compromising on the accuracy of the simulation results. It can be further noticed that the initial compliances plotted in Fig. 9 are different from those plotted in Fig. 8, although the same initial beam configuration is used in the shape and topology optimization analyses. This is, because in the shape optimization analysis the relative element densities are kept fixed at a value of 1, while the relative element densities in the topology optimization analysis are initially uniformly set equal to 0.4, in correspondence with the constraint on the volume fraction, f r 0.4.

Coupled shape and topology optimization case studies
Now that the shape and topology implementations have been individually tested and validated, the cantilever beam problem depicted in Fig. 4 is subjected to a coupled shape and topology procedure in accordance with the sequentially-coupled computational scheme depicted in Fig. 3. The variation in shape and topology in the longitudinal beam direction is described by using 8 cross sections, which are distributed along the beam length as shown in Fig. 14a. As explained before, for post-processing of the analysis results, the first cross section   at the clamped end of the cantilever beam (with index k 1) is taken the same as the adjacent second cross section (with index k 2), and thus does not contribute to the optimization procedure [6]. The locations of the control points of a cross section and the corresponding coordinates are given in Fig. 5 and Table 2. By leaving out the cross section related to the clamped end, the shape design variables a of the 7 remaining cross sections are assembled as a (a CS 2 ) T , (a CS 3 ) T , . . . , (a CS k ) T T with k 2, 3, . . . , 8, in which where X 1,1,k and Y 1,1,k are the x-and y-coordinates of control points C 1,1,k , and X 1,2,k and Y 2,1,k are the xand y-coordinates of control points C 1,2,k and C 2,1,k , respectively. Note that the control points at the corners of a cross section can move in both the x-and y-directions, while the control points at the half-width and half-height of the cross section boundaries can only move in the y-direction and x-direction, respectively. The control points of a cross section are assumed to be line-symmetric with respect to the x-and y-axes, so that the locations of control points C 3,1,k , C 3,2,k , C 1,3,k , C 2,3,k and C 3,3,k are prescribed by the corresponding control points above. The lower and upper bounds of a CS k are, respectively, defined as lb [−1.5, −1.5, −2, −2] and ub [−0.2, −0.2, −0.2, −0.2]. The parametric domain of the NURBS surface that defines a cross section is divided into 12 × 20 equally sized, rectangular parts, which define the grid for the 2D auxiliary FEM mesh. In the longitudinal direction, the model uses a discretization by 120 two-node beam elements, for which the z-coordinates are given by Eq. (93) with l 121. The stiffness k b e of each beam element e is computed from the stiffness of the corresponding cross section K s e via Eq. (34). Accordingly, under a stepwise change in cross-sectional properties, the optimized solution corresponds to a shape that varies linearly in the longitudinal direction. The volume fraction in Eq. (58) is defined as f r 0.4. The structural compliance c and the beam volume V are calculated from the beam FEM model. The SQP and IP algorithms in the MATLAB solver fmincon are used to iteratively update the shape design variables a and element relative densities ρ appearing in Eq. (58). The convergence criteria adopted for the individual shape and topology optimization steps are the same as those defined in Sect. 6.1. In addition, the outer loop of the coupled optimization procedure, which has been illustrated in the flowchart in Fig. 3, is assumed to be converged when the structural compliance meets the condition | c w+1 − c w |≤1e−5, where w is the iteration number of the outer loop. Finally, the shape and topology sensitivities of the objective are computed with the expressions given in Sect. 5.1, whereby the 2D sensitivity filter formulated in Eq. (78) is applied in the topology optimization procedure. The filter size is set to r min 1.5 √ V e , where V e is the area of element e [30]. The convergence behavior for the structural compliance c is shown in Fig. 12. The CTSO solution procedure that starts with a topology optimization step takes more iterations (w 10) than the CSTO solution procedure that starts with a shape optimization step (w 7). Figure 13 shows the evolution of the first cross section at the clamped beam end during the convergence process for the CSTO and CTSO update schemes. The evolution of the cross section clearly visualizes that the CTSO and CSTO update sequences during the convergence process follow a different search path in the solution space of design variables, eventually leading to somewhat different local minima (or optimized compliances), i.e., c 10.467 Nm (CSTO) and c 8.516 Nm (CTSO). Figures 14 and 15 illustrate the optimized designs predicted by CSTO and CTSO sequences. The two solution procedures result in a comparable design concept, namely a tapered beam with a hollow cross section composed of relatively thin vertical webs and thick horizontal flanges. The cross section computed by the CSTO solution procedure has a smaller height-to-width ratio than the cross section following from the CTSO procedure. Further, the beam structure following from the CSTO update procedure contains one cross section (with index k 4) with an extra vertical web.
In summary, in structural optimization problems, the landscape of solutions is determined by the objective function(s), the design variables and constraints, the characteristics of the boundary value problem and the FEM discretization, which typically leads to a large number of local minima, see e.g., [25]. The above comparison of the analyses performed with the CTSO and CSTO update schemes illustrates that the convergence speed  Fig. 13 Evolution of the structural compliance c and the geometry of the first cross section at the clamped beam end during the convergence process for the CSTO and CTSO update schemes. The NURBS control points are indicated by the black dots, and w is the iteration of the outer loop of the update scheme, see Fig. 3 (b) Coupled shape topology o y x z (c) Coupled topology shape (a) Fig. 14 Optimized control nets of 2.5D beam model. Initial design (a) and solutions obtained by coupled shape and topology optimization using CSTO (b) and CTSO (c) update schemes and the specific local minimum computed (i.e., the optimized compliance) are to some extent sensitive to the algorithmic features of the numerical update scheme applied, but that the final, optimized beam configurations of the two update schemes nevertheless are qualitatively similar.

Conclusions
A coupled gradient-based optimization framework is presented that simultaneously optimizes the outer shape and the internal topology of beam-type structures by using a staggered update procedure. The objective function refers to the structural compliance, which is evaluated using 2.5D and 3D FEM models. In the 2.5D model, standard two-node beam elements are used along the longitudinal direction of the beam, whereby the crosssectional properties of the beam elements are calculated from additional 2D FEM analyses. Conversely, in the 3D model, the beam geometry is simulated using 3D continuum elements. The shape of the beam-type structures is parameterized using NURBS, with the shape design variables being represented by the NURBS control points. The topology design variables are reflected by the relative densities assigned to each finite element. The design variables are iteratively updated by applying NURBS-based shape and density-based topology optimization techniques, which use analytic shape and topology sensitivities that ensure an accurate and computationally efficient solution procedure. A comparison study of a cantilever beam problem subjected to pure shape optimization and pure topology optimization illustrates that the 2.5D and 3D beam models lead to similar shape and topology designs, but that the 2.5D beam model has a significantly higher computational efficiency. Specifically, the computational times for the 2.5D model are about a factor 70 (shape optimization) and 1.4 (topology optimization) lower than for the 3D model, which indicates that in the coupled optimization approach the optimization of the shape provides the largest contribution to the higher computational efficiency of the 2.5D model. The coupled shape and topology optimization analysis subsequently performed on the 2.5D cantilever beam model demonstrates that the specific order at which the alternating shape and topology optimization increments are performed in the staggered update procedure turns out to have some influence on the final computational result for the boundary value problem considered. Further, the convergence speed for the optimization procedure starting with an incremental topology optimization step appears to be somewhat lower than that resulting from starting with an incremental shape optimization step, although it may be reasonably expected that this feature is problem dependent. Despite these differences, the final beam structures following from the two staggered update schemes illustrate how shape and topology can be efficiently optimized in an integrated, coupled fashion.