Isogeometric sizing and shape optimization of 3D beams and lattice structures at large deformations

A computational method for optimizing the shape of the centerline curve and the spatial variation of geometric and material sizing parameters of the cross-sections of elastic, 3-dimensional beams and beam structures subject to large deformations is presented in this work. The approach is based on the concept of isogeometric analysis, i.e., the representation of geometry and the discretization of the numerical solution using spline functions. Here, mixed isogeometric collocation methods are used to discretize the geometrically exact 3D beam model. These spline representations are extended to the parameterization of the design variables, which are the initial centerline curves of the beams, as well as cross-sectional sizing properties, which may be varying along the beam axis and can be functionally graded through the cross-sections. To tailor the mechanical deformation behavior of a beam or beam structure, a nonlinear optimization problem is formulated and solved using gradient-based methods. For this purpose, all required gradients and sensitivities are derived analytically. The potential of this holistic design optimization approach is demonstrated in application to tailoring of elastic metamaterials and beam lattice structures, as well as 4D printing of multi-material laminate beams.


Introduction
In recent years, the development of advanced and additive manufacturing (AM) technologies has substantially expedited the development of computational design and optimization methods. Besides the ability to realize structures with complex topologies and shapes, AM also provides unprecedented design freedom to integrate multiple materials and graded properties into monolithically fabricated parts (Gibson et al. 2015). So far, many AM technologies and applications, as well as the computational tools being developed therefore focus on stiff structures in the linear elastic regime, e.g., for load-carrying parts and lightweighting. However, in recent years tremendous progress has been made in particular in the fabrication of soft, flexible, and multi-functional structures, as well as architected metamaterials (Kadic et al. 2019;Surjadi et al. 2019). This includes, for instance, flexible and highly stretchable (Jiang and Wang 2016;, buckling-dominated (Liu et al. 2016;Janbaz et al. 2019;Jamshidian et al. 2020) or 4D printed, i.e., actively deforming, self-assembling, or reconfigurable (Ding et al. 2018;Boley et al. 2019) beam lattices and (micro-)structures. To exploit these capabilities and rationally design structures and architected metamaterials with, e.g., functionally graded material distributions (Kuang et al. 2019), shape-matching deformation behavior (Mirzaali et al. 2018), or anisotropic and curved lattice members (Vangelatos et al. 2020), computational design optimization methods and tools are required that can efficiently model free-form beam structures fabricated from soft materials and subject to large (elastic) deformations.
Even though spline-based and computer-aided design (CAD)-integrated optimization methods were already proposed much earlier (Braibant and Fleury 1984;Olhoff et al. 1991;Schramm and Pilkey 1993), many advancements in structural simulation and optimization have been made in recent years through the adoption of the isogeometric analysis (IGA) concept . Meanwhile, this idea of a unified representation of geometry and numerical discretization using basis splines (B-splines), non-uniform rational basis splines (NURBS), or variants thereof have already been applied in virtually all areas of computational mechanics and engineering (Cottrell et al. 2006). Already soon after its inception, the advantages of IGA for shape optimization problems have been exploited (Wall et al. 2008;Cho and Ha 2009), since IGA allows to use the control points of the spline geometry directly as design variables for optimization. Subsequently, IGA methods have also been developed for shape optimization of linear (Nagy et al. 2010(Nagy et al. , 2011Hosseini et al. 2018) and nonlinear beams (Radaelli and Herder 2014;. In the context of shape optimization, the derivation of sensitivities and gradients using analytical or semi-analytical approaches has received great attention (Cho and Ha 2009;Kiendl et al. 2014;Qian 2010;Fußeder et al. 2015). Moreover, spline representations were also extended to material parameterizations in the context of material distribution (Nagy et al. 2013;Weeger et al. 2016;Izzi et al. 2021) and topology optimization formulations (Seo et al. 2010;Kang and Youn 2016;Costa et al. 2019). In general, the advantage of the parameterization of the optimization variables using splines is that these parameterizations are smooth, have local support, and can be chosen independently from the resolution of the discretization.
Focusing on the optimal design of beams and lattice structures, most previous works are restricted to stiff structures in the linear elastic, infinitesimal deformation regime, e.g., for load-carrying parts and light-weighting applications (Haftka and Grandhi 1986). For instance, sizing, shape, and topology optimization methods were developed for the cross-sectional design of thin-walled or reinforced beams (Schramm et al. 1995;Zhang et al. 2009;Amir and Shakour 2018). For the topological design of truss structures, the ground structure method was introduced (Hagishita and Ohsaki 2009;Zegard and Paulino 2015) and also extended to nonlinear 2D beams (Changizi and Warn 2020). In (Mergel et al. 2014), nonlinear beam formulations were employed for the optimization of adhesive microstructures. A sheardeformable beam model with linearized kinematics and an isogeometric finite element discretization was employed for the optimal design of band gaps of lattice structures in . There, also the beam centerline shapes and cross-sectional sizing parameters were introduced as design variables.
In this work, we essentially combine the ideas of our previous works on isogeometric analysis for the discretization of a geometrically exact and nonlinear 3D beam model (Weeger et al. 2017) with the spline parameterization of cross-sectional properties (Weeger et al. 2016 and of the centerline curve ) for a holistic approach to the computational design optimization of beams and beam structures. Thus, by shape we refer here to the initial centerline curve of a beam and by sizing we refer to the geometric and material parameters of the cross-section (radius, Young's modulus, etc.), which can be varying along the centerline and also be functionally graded through the cross-section, i.e., be composed of different materials. To the best of our knowledge, so far no method has been presented yet that would facilitate such optimal design of beam structures or metamaterials with curved, free-form shapes, and varying graded cross-section properties at large deformations. For the shape sensitivity analysis and optimization of geometrically exact, nonlinear beams, only  introduced another approach. It is based on an isogeometric finite element discretization, a multiplicative update of rotations, and the smallest rotation method for determining the initial cross-sectional orientations.
The further outline of this manuscript is as follows: The 3D beam model with the general constitutive model for functionally graded cross-sections is introduced in Section Geometrically exact 3D beam model with functionally graded cross-sections. Then, the isogeometric analysis concept using B-splines or NURBS for the parameterization of the initial shape and cross-sectional sizing properties, as well as for the discretization of the beam model with mixed collocation methods is described in Section Isogeometric parameterization, discretization, and mixed collocation. In Section Isogeometric design and shape optimization, the combined design and shape optimization problem is formulated and the required gradients and sensitivities are derived. In Section Numerical results, the numerical implementation of the method is verified and it is demonstrated in several applications. Finally, the paper concludes with a summary of methods and results presented in Section Summary and conclusions.

Geometrically exact 3D beam model with functionally graded cross-sections
The mechanical modeling of the elastic deformation of 3-dimensional beams and beam structures in this work is based on the shear-deformable, geometrically exact beam model, which is usually referred to as the Cosserat rod, Timoshenko, or Simo-Reissner beam theory (Antman 2005;Simo 1985). It accounts for large deformations and rotational changes of moderately thick beams. However, due to the use of the linear elastic constitutive model and rigid cross-sections, it is only valid for small strains and stresses.

Kinematics
Kinematically, 3-dimensional beams are described as framed curves. The initial, undeformed configuration of a beam is given by a 3D curve r ∶ [0, L] → ℝ 3 , which describes its centerline, and an orthonormal field of rotation matrices R ∶ [0, L] → SO(3) , which describes the spatial orientation of the cross-section at each point along the centerline. L is the length of the initial centerline curve r , which must be arclength parameterized, i.e., ‖r � (s)‖ = ‖ dr ds ‖ = 1 ∀s ∈ [0, L]. Although the initial orientation of the cross-section can be defined arbitrarily, it is usually chosen such that the cross-section is orthogonal to the curve, which here defines g 3 (s) =r � (s) ∀s ∈ [0, L] , where R = (g 1 ,g 2 ,g 3 ) . This convention still yields g 1 ,g 2 as arbitrary up to a rotation around g 3 . Here, we always use the so-called Bishop frame for initially curved beams, which minimizes the torsion in the initial configuration, see ) for further details.
The current deformed configuration of a beam is described by the centerline curve r ∶ [0, L] → ℝ 3 and the orthonormal field R ∶ [0, L] → SO(3) , see Fig. 1 for illustration. The kinematic model then includes two types of strain vectors, the normal and shear strains ∈ ℝ 3 , as well as the bending and torsion strains ∈ ℝ 3 : Here k is the curvature vector: where axl refers to the axial or cross-product vector of a skew-symmetric matrix. (1) (2) k(s) =

Constitutive model for functionally graded beams
Assuming that strains and stresses are small and that the Bernoulli hypothesis holds, i.e., cross-sections remain plain and undeformed, a linear elastic constitutive model is used to relate the strain measures to stress resultants, which reads in its most general form as: where ∈ ℝ 3 is the vector of the normal and shear stresses and ∈ ℝ 3 is the vector of bending and torsion stresses. The constitutive matrices A, B, C ∈ ℝ 3×3 take the following forms: The topology and material constitution of the beam crosssection, which may be functionally graded or layered, define the coefficients of these matrices, see (Bîrsan et al. 2012) for a derivation. Without loss of generality, here we only consider circular cross-sections with either homogeneous or bilayer laminate material distributions, see Fig. 2. For homogeneous, isotropic cross-sections the coefficients are as follows: where E > 0 is the Young's modulus of the material, G = E 2(1+ ) its shear modulus, and ∈ [0, 0.5] its Poisson's ratio. For circular cross-section, the geometric properties depend only on the radius r > 0 . Then, A = r 2 is the crosssection area, I 1 = I 2 = r 4 12 the second moment of area, J = I 1 + I 2 the polar moment of area, and k 1 = k 2 = 5 6 are the shear correction factors.

Equilibrium equations
The vectors of stress resultants, which are defined in the material configuration of the beam through the constitutive model (3), are then transformed to the spatial configuration, which yields the vectors of internal forces n ∈ ℝ 3 and moments m ∈ ℝ 3 : The balance equations of linear and angular momentum are then formulated in the current spatial configuration: Here, n ∶ [0, L] → ℝ 3 and m ∶ [0, L] → ℝ 3 denote externally applied distributed forces and moments, respectively. To formulate a boundary value problem with a uniquely defined solution in terms of r and R , this coupled system of nonlinear ordinary differential equations must be completed with boundary conditions for s = 0 and s = L . Essential or Dirichlet boundary conditions are formulated in terms of prescribed positions or rotations as r(s) =r s , R(s) =R s at s = 0, L and natural or Neumann boundary conditions in terms of prescribed forces or moments n(s) =n s , m(s) =m s at s = 0, L , where the quantities ⋅ denote the prescribed given values.

Isogeometric parameterization, discretization, and mixed collocation
To approximate the solution of the boundary value problem of the beam, we employ the concept of isogeometric analysis , i.e., the kinematic unknowns are discretized using spline shape functions, such as B-splines and NURBS. In the same way, the geometry, here the initial configuration of the beam, is represented using spline curves, which provide a direct link to CAD and also serves as a parameterization for shape optimization, see . Furthermore, the parameters that define the cross-sections are also parameterized as spline curves, which enable the design of tapered beams with axially varying cross-sections, see , and thus also design optimization.

Spline basics
B-splines and NURBS are commonly used in CAD to mathematically represent curves and surfaces. Here, we briefly review the most relevant definitions, terminology and properties of B-splines and NURBS, which will be used in the following sections. A detailed introduction into NURBS, including algorithms for their numerical implementation, can be found in (Piegl and Tiller 1997). B-splines are piece-wise polynomial functions of degree p ∈ ℕ and order p + 1 . n ∈ ℕ B-spline basis functions B i ∶Ω → [0, 1] are defined on a knot vector Ξ = 1 , … , m with m = n + p + 1 , e.g., using the Cox-De Boor recursion formula, see (Piegl and Tiller 1997). The knot vector is a non-decreasing sequence of knots i ∈ ℝ (i = 1, … , m), i ≤ i+1 (i = 1, … , m − 1) and the parameter domain Ω = [ 1 , m ] ⊂ ℝ is, without loss of generality, often chosen as Ω = [0, 1] . For two distinct knots i ≠ i+1 the half-open interval [ i , i+1 ) is denoted as the i-th knot span, which may also be called an element. The total number of nonzero knot spans or elements in Ξ is denoted by ∈ ℕ . Here, we employ only open knot vectors that are interpolatory at 1 and m . This means that the first and last knot have multiplicity p + 1 , while inner knots have multiplicity k with 1 ≤ k < p.
Important properties of B-splines, especially in the context of numerical methods, are that they are non-negative, i.e., B i ( ) ≥ 0 ∀ ∈Ω , form a partition of unity, i.e., ∑ n i=1 B i ( ) = 1 ∀ ∈Ω , and that they are linearly independent, i.e., ∑ n i=1 B i ( ) u i = 0 ∀ ∈Ω ⇔ u i = 0 ∀i = 1, … , n . Furthermore, B-splines are p-times continuously differentiable within the elements and (p − k)-times at knots of multiplicity k, i.e., B i ∈ C p−k (Ω) . Moreover, the basis functions have compact support, meaning that B i is nonzero only in [ i , i+p+1 ] , i.e., on a maximum of p + 1 elements.
NURBS are piece-wise rational functions, which can be built from n B-spline basis functions B i and associated weights w i > 0: Since the NURBS basis functions N i inherit the above-mentioned properties of B-splines and B-splines can be considered as NURBS with all weights being equal to 1, we use the symbols N i for both B-spline and NURBS basis functions interchangeably.
A spline curve c ∶Ω → ℝ d can then be defined as a linear combination of n spline basis functions N i with control points c i ∈ ℝ d : where d ≥ 1 is the dimension of the curve. The differentiability properties of the spline basis functions then also apply to spline curves.
B-splines and NURBS offer the concepts of h-, p-, and k-refinement in order to adjust the resolution of a spline space and associated geometries, such as curves without changing their shape, see Piegl and Tiller 1997) for details.

Spline parameterization of the initial configuration
Now the reference configuration of a 3D beam is parameterized using spline curves. As in (Weeger et al. 2017), we use unit quaternions q = (q 1 , q 2 , q 3 , q 4 ) ⊤ ∈ ℝ 4 , ‖q‖ = 1 for the parameterization of the frame: Thus, the spline parameterizations of the initial configuration are as follows: with n 0 spline basis functions N 0 i ∶Ω → ℝ , which are defined on a knot vector Ξ 0 with degree p 0 and 0 elements, as well as n 0 control points r i ∈ ℝ 3 ,q i ∈ ℝ 4 and, for NURBS, weights ẘ i > 0.
For the initial strains and curvatures in (1) to be welldefined, the initial centerline curve r and quaternion field q must be continuously differentiable w.r.t. the arc-length variable, i.e., r,q ∈ C 1 [0, L] . However, both spline curves are

R(s) ≡ R(q(s)).
(11) defined on an arbitrary parameter domain Ω and thus they are in general not arc-length parameterized. The arc-length derivatives of any vector field → y( ) ∶Ω → ℝ d depending on the spline parameterization, such as r,q , but also strains , , can be computed as: where J is the Jacobian of the parameterization: and ̇y = dy d denote the parametric derivative w.r.t. the spline parameter . Evaluating the strong form of the balance equations (7) would even require second-order arc-length derivatives; however, this is avoided in the mixed methods presented below in Section Mixed discretizations and collocation methods.

Parameterization of cross-section properties
The concept of isogeometric parameterization is now also extended to the geometric and material parameters that define the cross-sections, which allows to introduce functional grading and variation along the beam axis, see .
By u ∈ ℝ d u we denote the vector of parameters of the cross-section, which can be both material or geometric specifications. For instance, for a homogeneous circular crosssection, we consider the radius r and the Young's modulus E, i.e., u = (E, r) , and for a circular bilayer cross-section we could consider u = (E 0 , E 1 , r, h R , ) . To introduce axial variation, this parameter vector is now also expressed as a function of the arc-length parameter s of the beam -or rather the spline parameter -in terms of a spline curve: Here, n u spline basis functions N u i ∶Ω → ℝ are used, which are defined using a knot vector Ξ u with degree p u and u elements. n u control points of this curve are u i ∈ ℝ d u , where d u is the dimension of the parameter vector. Additionally, for NURBS n u weights w u i > 0 are required. Depending on the desired axially varying design of the beam or the design freedom wanted for optimization, this spline basis can be chosen arbitrarily. Thus, it does not necessarily have to be the same or a refined version of the one used to parameterize the initial configuration of the beam in (11).
This axial parameterization of the cross-section properties implies that the matrices in the linear elastic constitutive model (3) also become dependent on the position along the beam centerline: As detailed in , this additional dependency on the arc-length parameter also has to be considered when derivatives in the strong form of the balance equations (7) are evaluated. For the transformation from parametric to arc-length derivatives, (12) applies.

Mixed discretizations and collocation methods
To determine the deformed configuration, i.e., the solution of the boundary value problem of the balance equations of linear and angular momentum (7), the current centerline curve r and quaternions q are also discretized and approximated as spline curves r ≈ r h ∶Ω → ℝ 3 and q ≈ q h ∶Ω → ℝ 4 : Here, the n basis functions N i ∶Ω → ℝ typically refer to a p/h/k-refined version of the basis functions N i from (11), for NURBS with n refined weights w i > 0 , see . The corresponding n control points are Here, the unknown control points are obtained using isogeometric collocation approaches (Auricchio et al. 2010), which are based on the solution of the strong form of the balance equations as given in (7). Since shear locking affects the accuracy and convergence of numerical discretizations of shear-deformable beams with small to moderate thicknessto-length ratios, i.e., r∕L ≪ 1 , we employ two variants of mixed methods, in which also the stress resultants are discretized independently to avoid shear locking, see (Weeger et al. 2017: Mixed method. For homogeneous cross-sections, where the constitutive tensor B ≡ 0 vanishes, also the internal forces n(s) and moments m(s) are independently discretized as spline curves: Here, the same spline shape functions as above in (16) are used, see (Auricchio et al. 2013;Weeger et al. 2017) and

(15) A(s) ≡ A(u(s)), B(s) ≡ B(u(s)), C(s) ≡ C(u(s)).
(16) additional unknowns in terms of the coefficients n i , m i ∈ ℝ 3 are introduced. The discretizations from (16) and (17) are substituted into the strong form of the balance equations, the quaternion normalization condition, and the additional compatibility conditions of the stress resultants and are then collocated: As the n collocation points ̂k ∈Ω , we here use the Greville abscissae of the spline knot vector Ξ: which guarantees the stability of the method (Auricchio et al. 2013). To include the boundary conditions, f n (̂k) and f m (̂k) for k = 1, n are replaced with either the Dirichlet or Neumann boundary conditions. This mixed isogeometric collocation method yields a nonlinear system of N = 13n algebraic equations: For the solution of this nonlinear system we employ Newton's method, which uses the analytically derived tangent stiffness matrix K(x) = df ∕dx for linearization, see Appendix 1 for details.
Enhanced mixed method. For non-homogeneous crosssections, such as the circular bilayer, this mixed method may not fully resolve the locking problem due to the translational-rotational coupling with B ≠ 0 and the scale differences between the constitutive matrices A ∼ r 2 , C ∼ r 4 , and B ∼ r 3 . Thus, in ) an enhanced mixed method was introduced, in which the translational and rotational strain contributions to the internal stress resultants are separated: Subsequently, the respective internal force and moment contributions are discretized independently: This approach yields N = 19n unknowns, since in addition to r and q , 4 ⋅ 3n control points n i , n i , m i , m i ∈ ℝ 3 are introduced, which are gathered in the vectors n , n , m , m ∈ ℝ 3n . The collocated governing equations then read as: where again f n (̂k) and f m (̂k) for k = 1, n are replaced with either the Dirichlet or Neumann boundary conditions. This yields a system of 19n nonlinear equations and unknowns: Also here Newton's method is used to solve the nonlinear system of equations with an analytically derived tangent stiffness matrix K(x) , see Appendix 1 for details. In the following, we will use the notation f (x) = 0 interchangeably for both (20) and (24). As shown in , these mixed methods resolve the shear locking issue and yield optimal convergence behavior. Furthermore, they reduce the continuity requirements on the shape functions to N i ∈ C 1 [0, L] , compared to a primal collocation method, which would require N i ∈ C 2 [0, L] . While the computational effort for solving the linear systems with the tangent stiffness matrix in the Newton's method is, of course, larger for the mixed methods since they involve more degrees of freedom, the total computational effort is actually reduced, since evaluations of the nonlinear force residual vectors in (20) and (24) are faster, the methods require less Newton iterations until convergence and also allow larger load steps. The reduced order of differentiation is also beneficial for shape and design optimization since the analytical derivation of sensitivities is simplified and their numerical evaluation accelerated, as will be shown in Section Gradient and design sensitivities.
Beam structures. To model structures with interconnected beams, such as meshes and lattices, the nonlinear residual vectors f (x) of all beams are assembled into a global vector. Then, Dirichlet boundary conditions can be eliminated and interface conditions must be applied. Here, we enforce a rigid coupling, i.e., the current centerline positions r and changes of orientation ΔR = RR −1 must be equal for all n I beams that are interfacing with each other at a node. Furthermore, the equilibria of linear and angular momentum must be fulfilled at the interfaces. As in (Weeger et al. 2017), the equilibrium conditions are implemented by summing up the (signed) internal force contributions of each adjacent beam end-point into the collocated nonlinear force vector entry of the "first" beam of the interface: where where q * denotes the conjugate of a unit quaternion q. Iterative solution process. As already mentioned, Newton's method is used for the iterative solution of the nonlinear system f (x) = 0 . Now and in the following, f stands for the assembled residual vector of a beam structure, to which coupling conditions have already been applied and from which Dirichlet boundary conditions have been eliminated.
x denotes the combined vector of all degrees of freedom from all beams. Thus, also the tangent stiffness matrices of the beams are assembled into a global stiffness matrix K(x) = df ∕dx , which is also modified according to the interface conditions. Then, given a starting value x 0 , the solution vector is iteratively updated for i = 0, 1, 2, … as: until convergence is reached in terms of ‖Δx i ‖ < and ‖f (x i )‖ < . Typically, loads or prescribed Dirichlet boundary conditions are applied gradually over several load steps, for which the solution is computed using the converged solution from the previous load step as starting value x 0 .

Isogeometric design and shape optimization
The spline parameterizations of the centerline curve and the cross-sectional properties are now used as design variables for optimizing beams and beam structures. Thus, a nonlinear optimization problem is formulated and then solved using gradient-based optimization methods. The gradient of the objective function is computed using the adjoint method and for this purpose analytical sensitivities are derived. For conciseness, in the following we formulate the approach in terms of a single beam. Nevertheless, the extension to a beam structure is straightforward. Thus, quantities such as x, f , can be understood as referring to a single beam or to assembled vectors of a beam structure.

Nonlinear optimization problem
The discrete representation of the initial centerline curve r and the parameterization of the axial variation of the crosssectional properties u were introduced as spline curves in Section 3.2 and Section Parameterization of cross-section properties, see Eqs. (11) and (14), respectively. While we take the definitions of the spline spaces and shape functions N 0 i and N u i in terms of the respective knot vectors Ξ 0 and Ξ u and also the corresponding NURBS weights ẘ i , w u i as fixed, the control points r i and u i of both parameterizations serve as design variables for the mechanical optimization of a beam. Therefore, we summarize them into two vectors r = (r i ) i=1,…,n 0 ∈ ℝ N 0 with N 0 = 3n 0 and u = (u i ) i=1,…,n u ∈ ℝ N u with N u = d u n u . Note that including the weights of NURBS as design variables would offer even more design freedom, see (Costa et al. 2019), but also complicate the derivations and evaluations of the design sensitivities.
A nonlinear optimization problem is now mathematically formulated for the minimization of a scalar objective function g in terms of the design variables r and u: . The objective function g depends (often explicitly) on the state variables x in the equilibrium configuration of the beam, which has to be determined through solving f (x) = 0 , see (20) or (24), for the current design instance given by (r, u) . Thus, f (x;r, u) = 0 is here included in the formulation as an implicit constraint with a dependence also on r and u . Furthermore, the design variables are typically constrained by lower and upper bounds, r m ,r M ∈ ℝ N 0 and u m , u M ∈ ℝ N u , respectively. Additionally, other types of (nonlinear) inequality constraints, such as volume, parameterization, or directional constraints, can be defined in terms of a vector-valued function h ∶ ℝ N 0 × ℝ N u → ℝ N c . For the formulation of such constraints, we refer to (Weeger et al. 2019, Sect. 4.3). In the examples demonstrated here in Section Numerical results, no nonlinear inequality constraints are applied. However, they could be adopted if required by a concrete application of the method.
Since the evaluation of g requires the computation of the deformed configuration x for the current design (r, u) , i.e., the solution of the nonlinear equation system f (x) = 0 through iterative procedures that potentially require many load steps with several Newton iterations each, the objective function evaluation can generally be considered as computationally expensive, while the subsequent evaluation of its gradient dg∕du is relatively cheap. This motivates the use of gradient-based local optimization methods, in contrast to derivative-free global methods, which typically require many more function iterations (Haftka and Gürdal 1992).

Gradient and design sensitivities
Since the objective function g generally depends on the design variables also implicitly through the current configuration, its gradient must be calculated as the total derivative: While the explicit derivatives of g w.r.t. r, u and x can be derived in a straightforward manner and then be evaluated analytically, see Section Objective function, the implicit derivatives require the terms dx∕dr and dx∕du . Using the adjoint method, see (Haftka and Gürdal 1992;, it follows that Here, ∈ ℝ N is the solution of: where df ∕dx = K(x;r, u) is the tangent stiffness matrix evaluated for the deformed configuration x for the beam design given by (r, u) and df ∕dr = K̊r(x;r, u) and df ∕du = K u (x;r, u) are the design sensitivities of the nonlinear force vector f (x;r, u) . Note that the tangent stiffness matrix K is not symmetric for collocation methods and thus the numerical solution of (31) requires the transposition of K. Shape sensitivities. The analytical derivation of the matrix of shape sensitivities: with ◻ ∈ {n, m, q, , } was described in detail in ) for a primal collocation method. For the mixed methods used here, this evaluation greatly simplifies since only first-order arc-length derivatives appear in f . The required partial derivatives w.r.t. the control points r i are briefly summarized in the following. First, the partial derivatives of the arc-length derivatives of r and of any (vectorvalued) function y , that does not explicitly depend on r , are required: Then, it follows for the enhanced mixed method from (23) The required expressions for r i , r i can be found in (Weeger et al. 2019, Sect. 4.2.1). Note that these derivatives also include the partial derivatives of the initial cross-section frame R w.r.t. r . Since R and q depend on r implicitly through the Bishop frame, similar to (31), the computation of the corresponding design sensitivities requires another sparse linear system solve for each gradient computation, see (Weeger et al. 2019, Sect. 4.2.3) for details. However, although not as efficient as the explicit smallest rotation method used in , the computational cost is negligible especially for larger beam structures as this system is local to each beam. For the standard mixed method from (18), it must simply be used that Cross-sectional sizing sensitivities. The analytical derivation of the matrix of the sensitivities w.r.t. the design parameters of the cross-sections: is also simpler for the mixed methods compared to primal collocation methods, since the constitutive matrices only appear in the compatibility equations of the mixed variables: Instead of explicitly representing the third-order tensors, such as dA du ∈ ℝ 3×3×n u , or implementing their application to the strain vectors, such as dA du , we rather compute the partial derivatives w.r.t. the individual parameters gathered in u , which form the columns of the terms, such as For both cases of homogeneous circular cross-sections with u = (E, r) and circular bilayer laminate cross-sections with u = (E 0 , E 1 , r, h R , ) , these partial derivatives can be found in Appendices 1 and 2, respectively.

Remark 1
It is important to mention that B-spline and NURBS shape functions such as N i , N 0 i , N u i all have local support and only p + 1, p 0 + 1, p u + 1 shape functions are nonzero at each collocation point ̂k , see Section Spline basics. Thus, the influence of each design variable, i.e., each control point r i and u i , is only local and the expressions in (34) and (37) need only be evaluated if ̂k is within the support of N 0 i or N u i , respectively.

Objective function
In design, shape, and topology optimization of structures subject to small, linear elastic deformations, the typical objective is to maximize the stiffness of a structure, i.e., minimize its compliance, subject to constraints that limit the weight. However, at finite deformations maximizing the stiffness is typically not a desirable target. Instead, we formulate objective functions for matching desired deformed shapes or resulting forces and moments: The first term g r corresponds to the mismatch between the current deformed centerline curve r h and a desired target given as a spline curve r t ∶ [0, L] → ℝ 3 . The error is evaluated in terms of the squared vector norm and summed over the collocation points ̂k . Depending on the initial configuration given by (r, u) , the mismatch between r h and r t can initially be very large, which may result in undesirable local minima being obtained by gradient-based optimizers. Thus, another term g is introduced, which measures the difference between the curvature in the deformed state and the curvature of the target curve t ≡ (r t ) ∶ [0, L] → ℝ 3 . As the numerical examples in Section Numerical results will show, combining these two terms often yields better optimization results. Furthermore, we also introduce the term g p , which only considers matching of the end-points of the centerline curve. It is specified by the two control points r 1 and r n and their respective target points r t 1 , r t n ∈ ℝ 3 . The final terms g n and g m consider the matching of given values n t k , m t k ∈ ℝ 3 (k = 1, n) for the internal forces n and moments m at the beam ends, i.e., for ̂1 and ̂n . For the (enhanced) mixed formulations used here, it is simply n(̂1) = n 1 + n 1 ≡ n 1 , n(̂n) = n n + n n ≡ n n , etc. All of these terms can be enabled and weighted using coefficients C r , C , C 1 p , C n p , C 1 n , C n n , C 1 m , C n m ≥ 0 , depending on the desired formulation of the objective.
Except for , which requires an arc-length derivative, see (2), none of the terms and expressions in the formulation of the objective function in (38) depends explicitly on the design variables r and u . Thus, it is (38) g(r, u;x) = g r + g + g p + g n + g m with The gradient of g w.r.t. the control points of the discretization x can be analytically derived from (38) in a straightforward manner: with the nonzero partial derivatives being

Numerical results
In the following, we present the numerical application of the unified isogeometric sizing and shape optimization approach. For this purpose, the methods have been implemented in a C++ code, which is based on the isogeometric analysis library G+SMo (Jüttler et al. 2014) and the nonlinear optimization library NLopt (Johnson n.d.), from which we employ the gradient-based Method of Moving Asymptotes (MMA) and Sequential Least Squares Programming (SLSQP) algorithms.

2D cantilever beam
To demonstrate the capabilities of sizing and shape optimization of beams, we first discuss a rather simple 2D example, which is illustrated in Fig. 3. A cantilever beam of length L = 1 with a circular cross-section with radius r = 0.02 , Young's modulus E = 10 8 , and Poisson's ratio = 0.5 is subject to a point force n L = (0, 0, −20) ⊤ at its free end. Numerically, the beam is discretized using B-splines with . p = 6, = 16 . Now, the beam shall be optimized such that its deformed configuration shows a constant curvature, i.e., is a circular arc, as it would result from applying a moment m L = (0, 10, 0) ⊤ instead of the force, see Fig. 3a. Using the shape-matching objective function g = g r + g with C r = 10, C = 0.2 , see (38), the initial objective value is g = 0.76. For the following optimizations, the cross-sectional properties and the initial centerline curve are parameterized with p u = p 0 = 3, u = 0 = 2, n u = n 0 = 5 . Then, either only the Young's modulus, u ≡ (E) , only the radius, u ≡ (r) , only the shape of the centerline curve, r , both modulus and radius, u ≡ (E, r) , or all three variables combined, u ≡ (E, r),r , are optimized. As bounds for the design variables, 10 7 ≤ E i ≤ 10 9 , 0.1 ≤ r ≤ 0.4 , and (0, −0.05, −0.5) ⊤ ≤r i ≤ (1.5, 0.05, 0.5) ⊤ are chosen. Since the beam must be clamped at its left end, for shape optimization the first two control points of the centerline curve are excluded, i.e., r = (r 3 , … ,r n 0 ) . As can be seen in Fig. 3b-3f, in all five cases the deformed shapes of the optimized beam designs resulting from n L are visually indistinguishable from the target shape and overlap it almost perfectly. Similar results are also achieved when the homogeneous cross-section is replaced by a circular bilayer, as is  shown in Fig. 3g-3i. Here, the initial design is chosen with E 0 = E 1 = 10 8 , r = 0.02, h R = 0.5 , i.e., equivalent to the homogeneous beam, and then only the moduli, u = (E 0 , E 1 ) , all cross-section parameters, u = (E 0 , E 1 , r, h R ) , or cross-section and shape, u = (E 0 , E 1 , r, h R ),r , are optimized. For the verification of the analytical derivation and numerical implementation of the sensitivities and gradients as detailed in Section Gradient and design sensitivities, Table 1 shows a comparison of the gradients for optimization of (E, r,r) with forward finite differences with a spacing of 10 −5 . The values agree very well, suggesting that the implementation is correct. Note that the objective function value here is g = O(1) and that the order of magnitude of the gradients is inverse to the order of the respective design variables, where E = O(10 8 ) , r = O(10 −2 ) , r i,1 =r i,3 = O(1) , and r i,2 = O(10 −1 ) . Thus, to improve the effectiveness of the numerical optimizers, in the implementation of the method all design variables are normalized to the range [0, 1].
The convergence behavior of the gradient-based optimization procedures is shown in Fig. 4 in terms of the evolution of g over the iterations. In all cases, convergence is reached in less than 40 iterations, but the final values of g vary substantially. When optimizing E and/or r, the target deformations can be matched well, as Fig. 3 shows. However, although g r becomes very small, the numerical values of the resulting curvatures actually match not as well (even this is not visible) and thus the objective is only reduced to g ≈ g ≈ 0.1 . The combined optimization of E and r finds an optimum with a slightly lower g-value, but only when r is (also) optimized, also the curvature can be matched well, yielding g ≈ 10 −4 . Here, we have also included a shape optimization that uses a NURBS parameterization with nonuniform weights, which is obtained by refining an initially quadratic spline with one element and weights {1, 1∕ √ 2, 1} to p 0 = 3, 0 = 2, n 0 = 5 . The optimized shape coincides with Fig. 3d and also the optimization behavior shown in Fig. 4 is similar to the shape optimization with uniform weights. The fast convergence behaviors shown in Fig. 4 also indicate the correct derivation and implementation of the gradients and sensitivities. The spiky behaviors most likely stem from the optimizer first attempting a design update with a full gradient step and then reducing the step size. Note that this behavior highly depends on the scaling of the objective function and the design variables, as well as on the choice of optimizer and maximum step sizes.
Interestingly, the combined optimization of E, r,r yields a slightly higher value for g than the optimization of r alone. Similarly, also the bilayer result for E 0 , E 1 , r, h R ,r is worse than the shape optimization alone. In these cases, probably less optimal local minima are found due to the strong nonconvexity of the optimization problems. In particular, for shape optimization, too much design freedom increases the risk of gradient-based optimizers being trapped in suboptimal local minima. This is not shown for conciseness, but while refining p u or u would slightly improve the optimal solution, further refinement of p 0 and especially 0 rather lead to less optimal local minima. Thus, the degree and number of elements should be kept rather low in particular for shape optimization and stabilizing constraints such as length or "speed of parameterization" constraints may be applied, see .
Altogether, we conclude that cross-section parameter and centerline curve optimization and the combination thereof all delivered good results in this demonstration problem.

Lattice structure
Beam lattices and metamaterials have become ubiquitous through advances in additive manufacturing technologies. Besides stiffness-to-volume ratios, also the large deformation and buckling behavior of soft and flexible microstructures are gaining increasing attention.
As an example, we consider the optimization of the deformation behavior of a lattice structure consisting of 3 × 3 × 3 truncated octahedron (Kelvin foam) unit cells, see Fig. 5a. The cubic unit cells each have a cell size of 10 × 10 × 10 , such that the total 756 beams each have length L = 2.5 √ 2 ≈ 3.5355 . The radii of the homogeneous circular cross-sections are initially set to r = 0.25 and the material Fig. 4 Convergence behavior for the optimization of the 2D cantilever beam parameters to E = 10, = 0.45 . This lattice structure is clamped at the bottom ( z = 0 ) with 0-displacement boundary conditions and at each beam end-point at the top ( z = 30 ) a tensile force n = (0, 0, 0.1) is applied. The deformed configuration is simulated with p = 6, = 12, n = 18 (total N = 78, 300 DOFs) and also shown in Fig. 5a. As can be seen, the z-direction displacements of the top points are almost homogeneous with u z ≈ 15.5. Now, the objective for optimization is to achieve, given the same applied forces, an inhomogeneous displacement of the top points such that u z (x, y) = 6 ⋅ (1 + x∕60 + y∕60) , which is indicated by the green planes in Fig. 5. Thus, the objective function is formulated with the end-point position mismatch g p , see (38), but here only includes the top points and only considers the position in z-direction. First, a sizing optimization of only the beam radii with u = (r) and p u = 2, u = 1, n u = 3 is attempted with bounds 0.125 ≤ r ≤ 0.375 and a total of N u = 2268 design variables. Then, a pure shape optimization of the centerlines r with p 0 = 2, 0 = 1, n 0 = 3 is performed with a total of N 0 = 4536 design variables, since the curves at the bottom ( z = 0 ) and top ( z = 30 ) are not optimized. Due to the large number of design variables, the MMA optimizer is used. In both cases, it exhibits good convergence behavior and reduces the objective function from g ≈ 1 to g ≈ 10 −5 in 40 iterations, see Fig. 6. The resulting optimized undeformed and deformed configurations are shown in Fig. 5b and Fig. 5c, respectively. It can be seen that the optimized lattices exhibit the desired deformation behavior under the applied tensile forces since the deformations of the top points align with the green surface. For the sizing optimization, the local stiffening and softening of the structures in form of increased and decreased radii is clearly visible in Fig. 5b. The shape changes in Fig. 5c are more subtle, with only beams in the upper part becoming significantly curved. Since both sizing and shape optimization yield very good results, a combined optimization is not attempted here.

Metamaterial unit cell
Next, we want to demonstrate the application of our approach to tailoring the compressive behavior of a periodic lattice metamaterial that avoids instantaneous buckling. As initial design we use the unit cell of a body-centered cubic (BCC) lattice of cell size 10 × 10 × 10 , which consists of 8 straight beams connecting the center of the cell with its corners, see Fig. 7a. The beams all have a uniform cross-section with Young's modulus E = 10 , Poisson's ratio = 0.45 , and radius r = 0.2 . For the simulation, they are discretized with p = 6, = 12, n = 18 . The effective behavior of the microstructure is homogenized by applying a uniaxial deformation to the corner nodes of the cell and enforcing periodic rotations at all corners. The response force f in x 1 -direction, which is the sum of the end-point forces of the 4 beams on the right face of the unit cell, compare g n in (38), is plotted in Fig. 7b. As can be seen, when the applied deformation is equivalent to a compressive strain of ≈ −3% , an instability occurs. Due to strut buckling the force-displacement behavior instantaneously softens significantly. Now, the goal is to design a unit cell that has the same initial stiffness, but avoids this instability and exhibits a smooth softening behavior, as indicated by the target curve in Fig. 7b. To solve this force-matching problem, we take the shapes of the centerline curves r and the radii of the crosssections u = (r) as design variables, which are parameterized for all beams using p u = p 0 = 3, u = 0 = 2, n u = n 0 = 5 . In total, this yields 112 design variables, as N u = 8 ⋅ 5 = 40 and N 0 = 8 ⋅ 3 ⋅ 3 = 72 (since the first and last control point of each beam are excluded from shape optimization). Here, an optimal microstructure, which almost perfectly fits the desired target response force curve, see Fig. 7b, is found by the SLSQP optimizer in only 11 iterations. Figure 7c shows the optimized unit cell with curved beams and axially varying thickened radii.

Active direct 4D-printed beam structures
The direct 4D printing approach presented in (Ding et al. 2018) allows to fabricate structures that actively deform. Essentially, this is achieved by incorporating an eigenstrain into a material through the 3D printing process, which is activated only upon heating due to the temperature dependence of material properties. When a bilayer laminate is fabricated from this active expanding material and a non-active material, the bilayer bends upon activating the eigenstrain through heating. The incorporation of these effects into the beam model is detailed in Appendix 1. Similarly, also many other 4D printing methods can be modeled.
Buckyball. As a first example, we want to optimize an initially flat beam mesh such that it assumes the perfectly spherical shape of a buckyball upon temperature activation, compare (Ding et al. 2018). The initial geometry consists of 102 beams of length L = 31 mm laid out in hexagonal patterns in the xy-plane, see all have a circular bilayer cross-section with radius r = 1.5 mm. The initial layer height ratio is set to h R = 0.44 , since based on the analytical model from (Ding et al. 2018) this should yield a resulting curvature of ≈ 12.9 m -1 and thus deform the mesh onto a sphere of radius 0.0775 m. However, as observed both numerically and experimentally in (Ding et al. 2018), due to the axial extension and coupling of bending and twisting at the joints, the resulting curvature and deformation of this initial configuration are much smaller than expected, see Fig. 8a. Now, the layer heights and the layer orientations shall be optimized such that the deformed beam mesh coincides with the target buckyball shape. Therefore, the design variables u = (h R , ) are parameterized using p u = 2, u = 1, n u = 3 and the shape-matching objective function g r is chosen with C r = 10 −3 . For the simulation of the actively deformed shapes, the beams are discretized using p = 8, = 8, n = 16 and the temperature is increased from T 0 = 25 • C to T 1 = 65 • C in 20 steps.
Using only 18 iterations, the SLSQP optimizer then converges from an initial objective function value of g = 0.55 to g = 2.2 ⋅ 10 −4 . As can be seen in Fig. 8b, the resulting deformed shape almost perfectly matches the desired target buckyball. The optimized values of the height ratio h R range from 0.444 to 0.517, yielding higher curvatures than the initial value of 0.44, and the layer orientations vary only slightly with ∈ [−0.12 , 0.04 ].
Helix. As another example for design optimization of active beams, we now want to tailor the height ratio and layer orientation of a long, direct 4D-printed rod such that it coils into a helix. As shown in Fig. 9a, the desired target helix has two coils, a constant diameter of 100 mm, and a pitch/spacing of 100 mm. This target curve can be exactly parameterized by a NURBS curve and has a cord length of (a) (b) Fig. 8 Optimization of a direct 4D-printed beam mesh to form a spherical buckyball upon activation. The flat initial shape is shown in gray and the target shape in transparent green. The deformed shapes are colored by the height ratio h R in both, the 3D views on the left and xz-planar views on the right 43 Page 16 of 22 659 mm. Considering the axial expansion of the 4D-printed rod, its initial length is chosen as 651 mm. With a bilayer cross-section with radius r = 1.5 mm that is initially oriented along the y-axis, the straight beam simply bends around the y-axis, as can be seen in Fig. 9a for h R = 0.2 and h R = 0.4 , the latter forming a circle. Now, the design variables u = (h R , ) are to be optimized such that the deformed shape coincides with the helix. For the design parameterization, p u = 2, u = 4, n u = 6 and p u = 2, u = 8, n u = 10 are applied. Due to the very large deformation from straight to helical shape to be obtained here, it is crucial that the shape-matching objective function includes both g r and g (here with C r = 0.1 and C = 10 −4 ). For the simulation, the beam is discretized with p = 8, = 64, n = 72 and (again due to the large deformation) 80 temperature increments are used from T 0 to T 1 .
As shown in Fig. 9d, about 35-50 iterations of the SLSQP optimizer are required to converge from objective function values g ≈ 1 to g ≈ 0.01 . For p u = 2, u = 4 , i.e., less design freedom, the minimal objective function value is slightly higher than when p u = 2, u = 8 is, which can also be seen from the agreement of deformed configurations with the target shape shown in Fig. 9b and Fig. 9c. The initial height ratio has only minor impact on the optimization, as the results for p u = 2, u = 8 with h R = 0.2 and h R = 0.4 are basically indistinguishable. The resulting distributions of the design variables are plotted in Fig. 9e. As could be expected, the desired constant curvature and pitch can be achieved with an (almost) constant h R and a that varies (almost) linearly along the beam length.

Summary and conclusions
We have presented a numerical framework for the combined design optimization of the cross-sectional sizing parameters and the shapes of the centerline curves of 3D beams and beam structures subject to large elastic deformations. The approach is based on the concept of isogeometric analysis, whereby not only the numerical discretization of the nonlinear differential equation but also the initial centerline curves and the cross-sectional sizing variables are parameterized using splines. The discretization of the strong form of the boundary value problem of the geometrically exact, shear-deformable 3D beam model is carried out using mixed isogeometric collocation formulations. These approaches exploit the efficiency of collocation methods. Compared to the primal approach, they alleviate shear locking for thin and functionally graded beams and greatly simplify the analytical derivation of shape and sizing design sensitivities for gradient-based optimization. We have demonstrated the applicability of this unified isogeometric design optimization framework to various shape and force-displacement curve matching problems of different complexities. In particular, it could be shown that the mechanical behavior of 3D lattice structures and periodic metamaterial unit cells, as well as the very large deformations of 4D-printed structures, could be optimized successfully. While in some problems a pure design or sizing optimization alone was sufficient to achieve the optimization target (e.g., lattice structure in Section Lattice structure), in others the combined optimization was necessary (e.g., metamaterial unit cell in Section Metamaterial unit cell). This also demonstrated the flexibility of the framework. In future research, this numerical optimization framework can be directly applied for the design of beam structures, in particular in the context of additive manufacturing of lattice structures, metamaterials, and multi-material fabrication. Furthermore, to make it even more versatile, it could be combined with a ground structure topology optimization method. A current limitation of the employed beam model is that it assumes rigid cross-sections and a linear constitutive model, which only applies to elastic materials and small strains. Thus, it would be valuable to generalize the approach to beam models with deformable cross-sections, hyperelastic, inelastic, or multi-physical material behaviors, since soft polymer or hydrogel materials often exhibit very large strains, viscoelastic, and damage effects and can, for instance, be combined with magneto-or electro-active battery materials. For applications with only moderately soft materials, stress constraints should be added to the optimization formulation so that the material behavior remains restricted to the linear elastic regime and nonlinear or failure effects can be avoided.

Constitutive coefficients for circular bilayer laminate cross-sections
According to (Bîrsan et al. 2012), it is A 12 = B 13 = B 23 = 0 in (4) for the case of cross-sections made from two isotropic materials. Then, the nonzero constitutive coefficients can be expressed as: where S ⊂ ℝ 2 is the geometric domain of the cross-section, k 1 and k 2 are the shear correction factors, G = E∕(2 + 2 ) and E are the location-dependent shear and Young's moduli for ∈ S , and is the Poisson's ratio. The torsion coefficient C 33 can typically not be expressed in a closed-form way.
For circular bilayer laminate cross-sections, we assume that 1 -coordinate of the parameterization of S is aligned with the layer interface, see Fig. Fig. 2. The circular domain S = { ∶ ‖ ‖ < r} is defined by the radius r. Using the layer height h L ∈ (−r, r) or the layer ratio h R = (h L + r)∕2r ∈ (0, 1) , the lower layer with Young's modulus E 0 is defined by the subdomain S 0 = { ∈ S ∶ 2 ≤ h L } and the upper layer with Young's modulus E 1 by S 1 = { ∈ S ∶ 2 > h L } . For both layers we assume that Poisson's ratios = 0 = 1 and mass densities = 0 = 1 are equal, which is the case for different kinds of polymers. As shown in , the constitutive coefficients can then be found as: where k 1 = k 2 = 5 6 and (43) 1 22 , C 33 ≡ C 33 (S 0 , E 0 , S 1 , E 1 ), with = 2 arccos(2h R − 1) . For the torsion constant C 33 , we use the approximation: If the layer interface is not oriented along 1 -direction, which we assume to coincide with the local g 1 -direction of the cross-section, it is rotated in 1 2 -plane by an angle ∈ ℝ . Then, the constitutive matrices can be expressed as: where A * , B * , C * are the "non-rotated" matrices as defined above using (44) and as also shown in ).

Extension to modeling of direct 4D printing
In (Ding et al. 2018), the modeling of the deformation of direct 4D-printed structures is introduced in the context of the geometrically exact, shear-deformable beam model. The constitutive model (3) is enhanced by a thermal strain, which includes thermal expansion and the release of the eigenstrain over a one-time temperature activation through heating from T 0 to T 1 > T 0 : where ΔT = T − T 0 for T ∈ [T 0 , T 1 ] . For a circular bilayer cross-section, the nonzero coefficients of the additional constitutive matrices A , B ∈ ℝ 3×3 take the form: compare (44) and (45). Here, 0 = t 0 + p 0 ∕(T 1 − T 0 ) and 1 = t 1 are the modified thermal expansion coefficients of material 0 and material 1, respectively. For the two materials used in (Ding et al. 2018), it holds that T 0 = 25 • C and T 1 = 65 • C , the material parameters (at T 1 ) are E 0 = 0.6 MPa, E 1 = 6.0 MPa, and = 0.5 , the thermal expansions coefficients are t 0 = 2.3 ⋅ 10 −4 K −1 and t 1 = 1.7 ⋅ 10 −4 K −1 , and the printing-induced eigenstrain of material 0 is p 0 = 0.044 . Rotation of the layering in the cross-section plane affect A , B as shown above in (47).
To obtain the final activated shape of a direct 4D-printed beam, the eigenstrain is applied through a temperaturecontrolled simulation where T is increased from T 0 to T 1 . Note that although it is essential for this 4D printing concept that the Young's moduli are temperature-dependent, we here only consider the material parameters at T 1 . Thus, the deformations during the activation/heating phase may not be physically correct, but the final deformations at T = T 1 are.

Linearization of collocated equilibrium equations
The iterative solution of the nonlinear systems of equations stemming from the mixed isogeometric collocation approaches require their linearization, i.e., the derivation and evaluation of the tangent stiffness matrices.
Mixed method. For the mixed method from (20), the tangent stiffness matrix is derived as: The required partial derivatives can be obtained as:

Circular bilayer laminate
For circular bilayer laminate cross-sections made from two isotropic materials, we allow as design parameters the two Young's moduli, the radius, the height ratio, and the layer orientation angle, i.e., u = (E 0 , E 1 , r, h R , ) , see A.1. Thus, we derive the sensitivities w.r.t. these parameters in the following.
Since the dependency on the layer orientation angle is realized through the rotational transformations in (47), it holds that with and ◻ ∈ {E 0 , E 1 , r, h R } being any of the other design parameters. The derivatives of the coefficients of the (non-rotated) matrices w.r.t. the Young's moduli E i , i = 0, 1 , can be easily inferred from (44):

Conflict of interest
The authors declare they have no competing financial interests.

Replication of results
The methods are presented and described in all steps required for implementation. The examples contain all necessary information about geometric, material, and loading parameters required for replication. Program codes can be shared by the author upon reasonable request.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long (61) as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.