Automated shape and thickness optimization for non-matching isogeometric shells using free-form deformation

Isogeometric analysis (IGA) has emerged as a promising approach in the field of structural optimization, benefiting from the seamless integration between the computer-aided design (CAD) geometry and the analysis model by employing non-uniform rational B-splines (NURBS) as basis functions. However, structural optimization for real-world CAD geometries consisting of multiple non-matching NURBS patches remains a challenging task. In this work, we propose a unified formulation for shape and thickness optimization of separately-parametrized shell structures by adopting the free-form deformation (FFD) technique, so that continuity with respect to design variables is preserved at patch intersections during optimization. Shell patches are modeled with iso-geometric Kirchho ff –Love theory and coupled using a penalty-based method in the analysis. We use Lagrange extraction to link the control points associated with the B-spline FFD block and shell patches, and we perform IGA using the same extraction matrices by taking advantage of existing finite element assembly procedures in the FEniCS partial di ff erential equation (PDE) solution library. Moreover, we enable automated analytical derivative computation by leveraging advanced code generation in FEniCS, thereby facilitating e ffi cient gradient-based optimization algorithms. The framework is validated using a collection of benchmark problems, demonstrating its applications to shape and thickness optimization of aircraft wings with complex shell layouts.


Introduction
Shell structures exhibit exceptional stiffness and strength-to-self-weight ratios, and are extensively employed in various engineering fields, such as aerospace, automotive, and marine engineering [1].The performance of such structures is greatly influenced by geometric and material properties.Thus, structural optimization plays a vital role in obtaining superior designs for shell structures.In this paper, we present an optimization approach based on free-form deformation (FFD) [2] to achieve the optimal shape and thickness distribution for isogeometric shell structures.
Structural analysis is involved in the optimization process to evaluate the response of the current design and guide the subsequent iterations.The finite element (FE) method [3] is a well-established approach used to approximate PDEs with Lagrange polynomial basis functions.However, the discretization of the computational domain through interconnected simple elements, known as meshing, for complicated geometries is the primary challenge in FE analysis.The FE mesh generation and related process can account for up to 80% of total analysis time [4].Alternatively, isogeometric analysis (IGA) [5,6] offers the possibility to bypass FE mesh generation by approximating the solution using the smooth non-uniform rational B-spline (NURBS) [7] basis functions.NURBS is the industrial standard widely used to represent computer-aided design (CAD) models, making IGA an ideal method to streamline the design-through-analysis process.
IGA has gained growing interest since its introduction not only due to the unified description between design and analysis models but also the regularity provided by NURBS basis functions.The smoothness in splines allows for direct discretization of the Kirchhoff-Love shell model [8], a fourth-order PDE that requires at least C 1 continuity with the Galerkin method.Various applications including wind turbines [9][10][11], bioprosthetic heart valves [12][13][14], and car hoods [15] have employed the isogeometric Kirchhoff-Love shell formulations [16][17][18][19][20] and demonstrated exceptional results.However, practical CAD geometries are often too complex to be represented by a single tensor product NURBS surface.To make the CAD models with multiple NURBS patches directly available for analysis, [21] introduces a fictitious strip to add bending stiffness at patch interfaces with conforming discretizations.Additionally, various methods, such as Nitsche's method [22,23], penalty method [24,25], and super-penalty method [26,27] have been applied to CAD models with non-matching NURBS surfaces, further expanding the applicability of IGA in dealing with complex geometries.
The seamless integration between CAD and analysis models in IGA makes it a natural choice for design optimization.The updated design in the optimization process can be precisely captured in the analysis, which in turn ensures accurate responses due to the exact geometry representation and excellent approximation capabilities of spline basis functions [28].Shape optimization using IGA has been investigated extensively [29][30][31][32][33][34].Many applications, such as beam structures [35], vibrating membranes [36], shell structures [37,38], and complex photonic crystals [39] show superior design by employing IGA in optimization.Moreover, topology optimization [40,41] also benefits from the same spline basis in design models and analysis.Shape optimization for shell structures with stiffeners has been explored in [42,43] using the FFD concept, a B-spline solid is extruded from a "master" part, which is stiffened with several "slave" stiffeners, to modify the shape of the whole shell structure.
In this work, we employ the open-source framework PENGoLINS [44] for automated IGA of non-matching shell structures using the penalty method in [24].The shape of the shell structure is updated through a trivariate B-spline FFD block, which encompasses the entire shell structure, without differentiating the "master" and "slave" parts.The FFD block modifies the Lagrange control points of all shell patches concurrently to preserve the surface-surface intersections.Subsequently, we obtain the resulting NURBS surfaces of shells using the Lagrange extraction technique [45], which is also implemented in the IGA using FE subroutines.Moreover, this approach is also applicable to shell thickness optimization where the thickness distribution is continuous at patch intersections.By integrating these two design variables, simultaneous shape and thickness optimization for non-matching shells can be effectively achieved.This combined optimization approach enables the exploration of complex design spaces while preserving the geometric integrity of the non-matching shell structure.To demonstrate the capability of the proposed method, we apply it to the design optimization of aircraft wings, effectively navigating the unconventional design space.
The structure of this paper is outlined as follows, we introduce commonly used notations and terminologies in Section 2 for reference.Section 3 reviews the penaltybased formulation for coupling of non-matching isogeometric Kirchhoff-Love shells and the computational algorithm for automated IGA of non-matching shell structures.
We present the FFD-based shape and thickness optimization approach and formulate the sensitivities using the Lagrange extraction technique in Section 4. The optimization approach is validated using a suite of benchmark problems in Section 5 and is applied to aircraft wing design optimization in Section 6, where superior design solutions are demonstrated.Lastly, we conclude the effectiveness of the proposed method and discuss potential future directions in Section 7.

Notation and terminology
In this section, we provide a summary of commonly used notions and terminologies for reference, as the formulations in the following sections can become complex due to the use of the extraction concept in IGA, interval quadrature meshes for coupling separate spline patches, and FFD B-spline blocks in optimization.
• p sh : order of spline surfaces for shell structures.• V I,IGA : IGA function space for shell patch S I .Superscript I indicates shell patch index.For single patch formulations, index I is neglected.• n I,IGA : number of degrees of freedom (DoFs) in V I,IGA .
• N I,IGA : spline basis functions in V I,IGA .
• u I,IGA : displacement in IGA DoFs for shell patch S I .• P I,IGA : NURBS control points for shell patch S I .• V I,FE : FE function space for shell patch S I .
• n I,FE : number of DoFs in V I,FE .
• N I,FE : Lagrange polynomial basis functions in V I,FE .• u I,FE : displacement in FE DoFs for shell patch S I .• P I,FE : Lagrange control points for shell patch S I .• M I : extraction matrix for I-th shell patch that represents spline basis functions with Lagrange basis functions.• ξ MI : parametric coordinates of intersection L with respect to shell patch S I .• V γM : function space of the interval quadrature mesh.Superscript γ ∈ {0, 1} indicates the derivative order of the functions interpolated from V I,FE to V γM .• u γMI : interpolated displacement with γ-th derivative from V I,FE to V γM .• P γMI : interpolated control point functions with γ-th derivative from V I,FE to V γM .• T γMI : interpolation or transfer matrix that interpolates γ-th derivative of functions from V I,FE to V γM .• p FFD : order of B-spline basis functions for the FFD block.
• P FFD : B-spline control points of FFD block.
In Section 3.1, we review the basic Kirchhoff-Love shell and penalty formulations, and the IGA and FE function spaces are not differentiated in the subsection.

Coupling of non-matching isogeometric shells
In aerospace structural applications, many components such as aircraft wings, empennage, and fuselage, can be modeled using shell theory.The Kirchhoff-Love shell model [16] is employed in this work with IGA discretization.Typically, aircraft CAD geometries are composed of a collection of NURBS patches, and we adopt a penalty-based coupling method [24] to perform analysis for shells composed of multiple patches with arbitrary intersections and assume St. Venant-Kirchhoff material model.We first review the basic kinematics of Kirchhoff-Love shell theory and the coupling formulation for a pair of intersecting shell patches.We then elucidate the computational procedures for the analysis of non-matching shell structures.

Penalty-based shell coupling
The Kirchhoff-Love shell model neglects transverse shear strains, with straight lines normal to the midsurface remaining straight and normal to the midsurface before and after the deformation, and shell thickness kept unchanged.Thus the displacement field of the 3D shell can be uniquely represented by the displacement of its midsurface.In the reference configuration, the midsurface of a single-patch shell can be represented by X(ξ) , where ξ = {ξ 1 , ξ 2 }, are parametric coordinates of the midsurface.In the deformed configuration with displacement u(ξ), the shell midsurface is expressed as Taking derivatives of midsurface with respect to parametric coordinates, we can obtain the covariant basis vectors in reference and deformed configurations on the tangent plane as and unit vectors normal to the midsurface read as Membrane strain ε and change of curvature κ can be derived from ( 2) and (3).Applying the appropriate material model, we can obtain the corresponding normal forces n and bending moments m.The internal virtual work of Kirchhoff-Love shell is given by where S is the midsurface of the shell.External virtual work is defined as where ρ is the density, B is the body force, and T is the traction.Ω is the 3D domain of the shell volume with a thickness of t about the midsurface, and Γ is the 2D surface boundary where T is applied.The principle of virtual work states (6) is a nonlinear equation and can be solved using the Newton-Raphson method, Here, we only go through the fundamental equations of the Kirchhoff-Love shell model, readers are referred to [17, Section 3] for detailed derivation.Practical complex shell structures typically consist of multiple NURBS patches.For example, a pair of intersecting shells joined at a certain angle, which is prevalent when modeling the stiffeners of an aircraft wing, cannot be described by one NURBS surface.To make the CAD geometry directly available for IGA, separate NURBS patches need to be coupled together so that both displacement continuity and joint angle conservation on the surface-surface intersections are maintained during deformation.
In this project, we use the penalty-based shell coupling formulation proposed by Herrema et al. [24], which offers a good balance between accuracy and simplicity.Consider two shells modeled with spline patches S A and S B that share one intersection L, as depicted in Figure 1.The unit vector a A t is tangent to the intersection in the deformed configuration.In the reference and deformed configurations, conormal vectors are defined as Subsequently, shell patches S A and S B are coupled through the following penalty energy where the α d term penalizes the differences of the displacement on the intersection to maintain C 0 continuity, and the α r term penalizes the change of the angle between S A and S B to keep joint angle conservation.Displacement and angle conservation penalty parameters α d and α r are defined based on a problem-independent, dimensionless penalty coefficient α scaled with element size and material properties.For isotropic materials, the penalty parameters are given by where E is Young's modulus, ν is Poisson ratio, t is shell thickness, and h is element size.Details about the derivation of penalty parameters can be found in [24, Section 2].
% Fig. 1: Spline patches S A and S B with one intersection L (indicated with red curves), where u and a 3 are displacement and unit normal vector of midsurface, a t and a n are unit tangent vector and unit conormal vector on the intersection.

Computational procedures for non-matching shells
In this section, we go through the computational procedures for structural analysis of isogeometric Kirchhoff-Love shells.The stiffness matrix of non-matching shells is formulated and used as the sensitivity in the design optimization in Section 4.

IGA of Kirchhoff-Love shells using extraction
The concept of extraction [45][46][47][48] is utilized in the implementation of IGA, whose spline basis functions can be represented exactly by the linear combination of Lagrange basis functions.These Lagrange basis functions can be used in the classical FEM, allowing IGA to be implemented using finite element software with pre-defined extraction operators.An open-source IGA Python library named tIGAr is developed by Kamensky et al. [49] using the finite element software FEniCS [50].Implementation and technical details are discussed in Section 4.4.In this section, we illustrate the basic mathematical operations and workflow of IGA using extraction.
To perform IGA, an extraction matrix M is generated to represent functions defined in spline function space V IGA with FE basis functions in V FE .The relation between these two sets of basis functions is given by where N IGA are IGA basis functions and N FE are FE basis functions.Each column of M is the line combination of N FE giving an IGA basis function.In the analysis, we first create an extraction matrix M and assemble the stiffness matrix K FE and force vector F FE in V FE using existing finite element subroutines.Then the displacement in V IGA is solved as with problem-specific boundary conditions applied to M T K FE M and M T F FE .
Remark 1 For the purpose of clarity, we assume control points of spline surfaces have unit weights.Therefore, rational spline basis functions are the same as homogeneous spline basis functions, both denoted as N IGA .In practice, weights need to be taken into consideration for correct geometric mapping and analysis.
For single patch Kirchhoff-Love shell analysis, stiffness matrix K FE is the second derivative of total work (∂u FE ) 2 from V FE to V IGA , and the formulation of IGA stiffness matrix can also be expressed as The right-hand side (RHS) of ( 12) for Kirchhoff-Love shell is equivalent to Therefore, ( 12) is recovered as the linear system in V IGA to solve for displacements increments in IGA DoFs,

Integration of penalty energy on patch intersections
For shell geometries that comprise multiple NURBS surfaces with non-matching intersections, the penalty energy is introduced to couple separate shell patches.In this section, we demonstrate the process to integrate the penalty energy (9) by using quadrature meshes.Again, taking two intersecting shells with one intersection as an illustrative example, the schematic configuration in parametric and physical spaces is shown in Figure 2. We generate a geometrically 2D, topologically 1D interval mesh Ω M , which is named quadrature mesh, in the parametric space to represent the intersection curve for the penalty energy integration.The parametric coordinates of the surface-surface intersection with respect to two spline patches are denoted with ξ MA and ξ MB .The overall algorithm to compute W AB pen on Ω M are outlined as follows: 1. Define function spaces V 0M and V 1M on Ω M with linear FE basis functions to create displacements, control point functions, and their first derivatives.2. Move quadrature mesh Ω M to parametric coordinate ξ MA to generate transfer matrices between function spaces of shell patch S A and quadrature mesh Ω M , where the configuration is shown in the upper left part of Figure 2. The following transfer matrices are created: (a) T 0MA : V A,FE → V 0M , which interpolates functions from V A,FE to V 0M .Each entry is defined as where N A,FE i is the i-th basis function in V A,FE , and ξ MA j is the j-th coordinate of Ω M .(b) T 1MA : V A,FE → V 1M , which interpolates the first derivative of functions from V A,FE to V 1M .T 1MA is defined as where N A,FE i, ξ is the parametric derivative of i-th basis function in V A,FE .3. Compute displacement u 0MA , geometric mapping P 0MA , and their first derivatives u 1MA , P 1MA from parametric domain of S A to Ω M using transfer matrices generated from last step, 4. Substitute interpolated quantities from previous step into (1)-( 3) to obtain basis vectors on Ω M before and after deformation: Remark 2 In this section, we assume S A and S B have the same thickness t M = t A = t B in penalty parameters, same for Young's modulus and Poisson ratio.For the thickness optimization problem discussed in Section 4.3, shell patches may have different thicknesses, or each shell patch has variable thickness distribution.In these scenarios, we use an identical method to the element size calculation to interpolate thickness from shell patches to quadrature mesh and take the average, i.e., t M = 1 2 (t MA + t MB ).
10. Based on (9), we can integrate the penalty energy on Ω M with quantities computed from previous steps, where J M is the associated line Jacobian mapping Ω M to L.

Assembly of non-matching system
With the coupled formulations for intersecting shells, we can perform IGA on the non-matching shells directly.For a pair of shells with one intersection illustrated in Figure 2, the total virtual energy in equilibrium is given by which can be solved by the Newton-Raphson method.We can solve the linearized problem of ( 21) to obtain the full displacement Parametric space Physical space Fig. 2: A pair of shell patches share one surface-surface intersection in the physical space, associated parametric surfaces are illustrated.A topologically 1D, geometrically 2D quadrature mesh (red interval mesh) is created in the parametric space and moved to parametric locations of the intersection with respect to two shells to create transfer matrices.Geometric and displacement quantities of shell patches are interpolated to the quadrature mesh, where the penalty energy is integrated.
Formulations for the derivative of W t can be obtained by means of the chain rule.Using S A for demonstration, the derivative of W t with respect to u A,IGA is given by The second derivatives of W in the left-hand side (LHS) of ( 22) can be computed with the following formulations.Taking S A for illustration, the diagonal block reads as (24) and the off-diagonal block can be expressed as In ( 23)-( 25), expressions of derivatives of W A can be found in [16,Section 3], and [24, Section 2.2] spells out derivatives of W AB pen .These derivatives are computed automatically by the implementation discussed in Section 4.4 using computer algebra and code generation capabilities in FEniCS.Substituting ( 23)-( 25) into (22), one can solve the displacements in IGA DoFs for S A and S B .
Remark 3 For clarity, we only demonstrate the non-matching system with two shell patches.However, ( 22) can be readily extended to shell geometries consisting of more than two surfaces.We refer interested readers to [44,Section 3.2] for the assembly of shell structures with m patches and arbitrary intersections.

FFD-based design optimization
CAD geometries of Kirchhoff-Love shells can be used for analysis without finite element mesh generation by employing formulations from Section 3.This provides attractive features for shape optimization of shell structures, where the discretization of shell patches stays unaltered during shape evolution.The updated geometry in optimization iterations stays consistent with the analysis model.As a result, this approach minimizes the effort required for geometry processing while simultaneously enhancing accuracy.
In the context of shape optimization for non-matching shell structures, it is crucial to ensure that updated shell patches remain properly connected.Failure to maintain this connectivity can result in separation or self-penetration of shell patches during the optimization iteration.Such issues would lead to false analysis and yield unrealistic optimal shapes.To tackle this challenge, we adopt the FFD-based [2] technique combined with Lagrange extraction to update shell geometry and demonstrate the workflow in Section 4.1.A comparable concept can be applied to thickness optimization to ensure continuous thickness distribution at the surface-surface intersections if needed.Sensitivities for both shape and thickness optimization are given in the subsequent sections.

Non-matching shells update through FFD block
We use a cylindrical roof consisting of four non-matching shell patches that are shown in the upper-left part of Figure 3 as an example to demonstrate the FFDbased shape optimization approach.Note that this approach can be applied to shell structures with an arbitrary number of patches.

Approximating shell patches with extraction
For the initial CAD geometry consisting of m Kirchhoff-Love shell patches, define a set of NURBS surface functions {S A (ξ), S B (ξ), . . ., S I m (ξ)}, and the I-th shell patch S I (ξ) is expressed as where N I,IGA (ξ) are the spline basis functions of degree p sh in V I,IGA .We omit degree p sh in the notation for clarity.P I,IGA are the NURBS control points for surface function S I .Using the extraction concept discussed in Section 3.2.1,NURBS surface function S I (ξ) can be represented with Lagrange polynomials as well, where N I,FE (ξ) are basis functions in the finite element function space V I,FE s with nodal interpolatory property, and P I,FE are Lagrange control points, or nodal values of S I .Plugging nodal coordinate ξ I,FE of V I,FE s into (27), coordinate of the NURBS surface S I is represented with nodal values in the discrete setting, S I (ξ I,FE ) = N I,FE (ξ I,FE )P I,FE = P I,FE . ( Based on (11), Lagrange control points can be obtained through the extraction matrix and NURBS control points.We have the following relation, S I (ξ I,FE ) = P I,FE = M I P I,IGA .(29)

Relating FFD block control points and shell geometry
The first step of Figure 3 illustrates the initial configuration of a collection of intersecting non-matching shell patches S , where red control nets are displayed.To enforce connectivity of the intersections during optimization, we immerse S in a trivariate B-spline block, which is refered to as an FFD block, and use control points of the FFD block as design variables.A schematic demonstration is shown in the second step of Figure 3.The FFD B-spline block is defined as where θ is the parametric coordinate of the FFD block, N FFD (θ) are B-spline solid basis functions of degree p FFD with knots vector, and P FFD are B-spline block control points.
To simplify formulation and implementation, we use an identity mapping for the FFD block B-spline block, so that the parametric coordinate coincides with the physical coordinate, Substituting ( 29) into (31), NURBS surfaces of the non-matching shells can be expressed using the FFD block basis functions and control points, In the continuous context of (32), shell patches will not separate in the final configuration as long as they are interconnected in the initial geometry.As the shape update of the FFD block is continuous, there is no relative movement between patches within the FFD block.In the discrete space, we can relate the NURBS control points of the shell patches to the control points of the FFD block, The control points of the NURBS surface P I,IGA of shell patches can be updated through the control points of the FFD block P FFD .Let N FFD ( PI,FE ) A I FFD , where PI,FE denotes Lagrange control points of spline patch I in the baseline configuration.Then P I,IGA can be computed as It is noted that we need to solve the system using Moore-Penrose pseudo inverse due to the non-square nature of the extraction matrix M I , which has dimensions of n I,FE × n I,IGA .For the extraction matrix, we have n I,FE > n I,IGA , which means that we are solving an overdetermined system.Therefore, P I,IGA is considered as a least square fit in (34) rather than an exact solution.The shape update strategy using FFD block is illustrated in the third step in Figure 3, and the resulting shell patches with control net are depicted in the fourth step.A comparison between the initial nonmatching cylindrical roof and updated NURBS surfaces is shown in Figure 4, where the surface-surface intersections keep overlapping within tolerance in the updated configuration.The procedures to update control points of non-matching shells with m patches are summarized as follows: 1.In the preprocessing step, generate sparse matrices of evaluation of FFD block B-spline basis functions at shells' Lagrange control points in the initial configuration A I FFD and Lagrange extraction matrices M I , for I ∈ {A, B, . . ., I m }. 2. At optimization iteration step i opt , obtain updated control points of the FFD block (P FFD ) i opt .Compute updated Lagrange control points P I,FE i opt for all shell patches, 3. Solve NURBS control points P I,IGA i opt at step i opt through Moore-Penrose pseudo inverse, 4. Perform IGA with updated shell geometry, evaluate objective function and derivatives if needed, then proceed with optimization iteration.
Though the control points of the shell patches are computed in the least square fit sense, the updated geometry can still retain the intersection with sufficient discretization.A sliced view of the intersection, in accompaniment with NURBS and Lagrange control points, between the right two spline patches S C and S D in the first step of Figure 3 is shown in Figure 5, where we use coarser discretizations to make the comparison clearer.In the updated configuration, the two cubic intersecting edges are still overlapping even with only 5 and 6 NURBS control points.Remark 4 Since there is no relative movement between intersecting spline patches within the FFD block, which can be achieved with adequate control points in the discrete context, parametric coordinates of surface-surface intersections remain unchanged during shape updates.Therefore, transfer matrices introduced in Section 3.2.2can be reused to interpolate data from spline patches to quadrature meshes when integrating penalty energies in the optimization iteration.These matrices only need to be generated once at the preprocessing stage.

Sensitivities for shape optimization
By utilizing the capabilities of direct analysis for non-matching isogeometric shells and incorporating FFD-based shape updates, we are able to conduct shape optimization for the shell structures in a seamless manner.The problem that optimizes the shape of non-matching shells can be formulated as follows, minimize f obj (P FFD ) subject to g con i g (P FFD ) ≤ 0 , for i g ∈ {1, 2, . . ., n g } h con i h (P FFD ) = 0 , for i h ∈ {1, 2, . . ., n h } P FFDl ≤ P FFD ≤ P FFDu , (37) where control points of the FFD block P FFD are design variables, f obj is the objective function, g con i g are inequality constraints, and h con i h are equality constraints.P FFDl and P FFDu are lower and upper limits for the design variables.
To perform gradient-based design optimization, we formulate total derivatives of the objective function with respect to design variables as where P IGA = P A,IGA P B,IGA . . .P I m ,IGA T is the full vector of NURBS control points, and similarly, u IGA = u A,IGA u B,IGA . . .u I m ,IGA T is the full vector of shell displacements in IGA DoFs.
Partial derivatives ∂ f obj ∂P IGA and ∂ f obj ∂u IGA in (38) can be computed and depend on the objective function combined with extraction matrices, where P FE = P A,FE P B,FE . . .
is a diagonal block matrix for global extraction.Calculation of partial derivatives in (39) is automated using FEniCS.Formulation for total derivative dP IGA dP FFD is introduced in (34).As for total derivative du IGA dP IGA , we have the implicit relation between P IGA and u IGA , where W t is the total energy of the non-matching shells defined in (21).Once an updated P IGA is obtained, the shell displacements need to be solved using (21) until the residual vector r t reaches a tolerance.Thus, r t is supposed to remain as 0 despite the change of P IGA , and we have the following derivative and the total derivative du IGA dP IGA in (38) is given by Partial derivative ∂R t ∂u IGA is equivalent to ∂ 2 W t ∂(u IGA ) 2 and is the stiffness matrix defined in (22).Analogously, we use a pair of shell patches to illustrate the formulation of partial derivative ∂R t ∂P IGA , Partial derivatives in (43) have identical expressions to ( 24) and ( 25).Extend partial derivatives in (43) to shell structures with an arbitrary number of patches, and substitute du IGA dP IGA in (38) with ( 42), we can obtain the total derivative of the shape optimization

Sensitivities for thickness optimization
The idea of FFD-based shape update can be applied to shell thickness optimization, where the shell thickness is treated as an extra field of the NURBS control points.We can use (34) to build the relation of the thickness between shell patches and FFD block, where t I,IGA is the thickness for shell S I in IGA DoFs, and t FFD is the corresponding thickness field of the FFD block.Subscript s in M I s and A I FFDs denotes matrices for scalar fields.Note that t FFD is not the actual thickness of the B-spline solid geometry, but an extra set of the control points on the FFD block to update the thickness of the non-matching shells.Accordingly, the identical shape update strategy Section 4.1.2is applicable to thickness update.FFD-based thickness optimization also offers the benefit that shell thickness remains continuous on the surface-surface intersections.
Replacing control points of the FFD block in (37) with t FFD , one can have the problem description of thickness optimization.Since both Kirchhoff-Love shell total work W A and W B , and penalty energy W AB pen involve shell thickness, the total derivative and associated partial derivatives of the thickness optimization problem can be acquired by replacing P FFD , P FE , P IGA with t FFD , t FE and t IGA in Eqs. ( 42)- (44), respectively.The total derivative of FFD-based thickness optimization reads as In some applications, one may choose to have a constant thickness for each shell patch.This can be easily achieved within the current framework by relating the shell thickness in IGA DoFs to one scalar value t I,const as where T is a unit column vector contains n I,IGA entries.The total derivative for piecewise constant thickness optimization can be obtained by replacing (46) with the following derivative, The FFD block is not needed in piecewise constant only thickness optimization.These two approaches can be combined together to achieve a more realistic design, where specific sections of the structure necessitate a continuous thickness distribution while constant thickness is better suited for other patches.In the implementation, shell patches can be separated into various groups.One group comprises the shell patches immersed within an FFD block allowing for a continuous thickness distribution.On the other hand, the shell patches not contained in an FFD block are assumed to have a constant thickness.Moreover, shell patches originating from different FFD blocks would exhibit discontinuous thickness at their intersections.Therefore, the combined thickness optimization approach provides more flexibility.

Implementation
This section introduces the open-source implementation for the proposed optimization method in Sections 4.1-4.3.We employ the open-source Python library PENGoLINS [44] for IGA of non-matching shell structures.PENGoLINS utilizes functionalities in pythonOCC [51] to approximate patch intersections in CAD geometries.Additionally, it leverages advanced code generation in FEniCS [50] and extraction operators in tIGAr [49] to perform IGA.CAD models in IGES or STEP formats can be imported into PENGoLINS, making them directly available for IGA.The code generation capabilities and computer algebra in FEniCS allow us to obtain partial derivatives in ( 23)-( 25), as well as (39) automatically.Integrating computed analytical derivatives with matrices M = diag(M A , M B , . . ., M I m ) and A FFD = diag(A A FFD , A B FFD , . . ., A I m FFD ), we can calculate the total derivatives in ( 44) and (46) for shape and thickness optimization.
With the availability of automated derivatives, we use the open-source framework CSDL [52] for conducting gradient-based large-scale optimization.Each partial derivative computed from FEniCS, in combination with predefined matrices, is passed to a CSDL sub-model, which contains explicit or implicit operations.A system-level model is created to connect all sub-models to enable efficient and modularized design optimization.For the FFD-based shape optimization problem, which involves minimizing the internal energy of a non-matching shell structure, four essential sub-models are required and listed in Table 1.All the derivatives required for these models are defined in Sections 3. 2  The models listed in Table 1 and other constraint models are modularized in a Python library named GOLDFISH (Gradient-based Optimization, Large-scale Design Framework for Isogeometric SHells).This library provides users with the option to use predefined models or create custom models to build the system-level model for their specific optimization needs.The library also includes an open-source optimizer, SLSQP [53], which can be utilized for shape or thickness optimization of non-matching shells.For optimization problems involving a large number of design variables, the SNOPT [54] optimizer is available for faster convergence.Furthermore, we have integrated the analysis framework with another optimization library OpenMDAO [55], which is developed by NASA.This integration further enhances the accessibility of the GOLDFISH library, enabling users to benefit from the collective expertise of the OpenMDAO community.The GOLDFISH library is hosted on a GitHub repository [56] and open to the public.

Benchmark Problems
A series of benchmark problems are considered to verify the effectiveness of the optimization method.Sections 5.1-5.3illustrate that baseline non-matching shell structures with arbitrary intersections are able to accurately converge to the analytical optimum.Section 5.4 studies the capability and flexibility of the framework for thickness optimization.

Arch shape optimization
An arch fixed at two ends and subjected to a constant downward load per unit horizontal length is modeled by a Kirchhoff-Love shell theory.Detailed problem definitions can be found in [38,Section 8].To test the effectiveness of FFD-based shape optimization for the non-matching shells approach, we model the arch using four NURBS patches with three intersections, where the arch geometry in the baseline configuration is shown in Figure 6a.We immerse the arch geometry in a trivariate B-spline block in the initial configuration, as is illustrated in Figure 6b.The analytical optimal solution is given by a quadratic parabola, where the ratio between the height of the arch and the horizontal distance of two fixed edges is 0.54779.We use quadratic B-spline for the FFD block in all three directions, p FFD = 2.The arch shell patches are described by cubic NURBS surfaces, p sh = 3, with 1086 DoFs in total.This benchmark problem minimizes the internal energy of the shell structure, with vertical coordinates of the control points of the FFD block being the design variables.From the control net in Figure 6b, it can be observed that there are 54 design variables.Two constraints are applied to this problem.The first constraint ensures that the lines in the FFD control net are parallel to the axial direction of the arch, keeping the arch devoid of tilting or twisting during the optimization process.The second constraint fixes the bottom layer of FFD control points so that the two edges of the arch remain in the initial position.We use the SLSQP optimizer with a tolerance of 10 −12 to perform the optimization, snapshots of the optimization iteration are demonstrated in Figure 7.
The arch converges to the analytical optimum shape after 40 iterations.The shape of the FFD block in the final configuration is shown in Figure 7.As anticipated, the optimized arch geometry is still contained in the FFD block.The height to base span ratio in this problem is measured as 0.54748, exhibiting a relative error of 0.057% compared to the exact value.Considering the coarse discretization of the arch geometry, the results are encouraging.Fig. 7: Snapshots of non-matching arch shape optimization.

Tube shape optimization
A square tube in the baseline configuration is subjected to an internal constant pressure.The optimal shape is given by a cylindrical tube [38,Section 7].We use four cubic NURBS surfaces to model one-quarter of the square tube, where the initial geometry is illustrated in Figure 8a.The square tube geometry contains four nonmatching intersections with 1035 DoFs in total, and symmetric boundary conditions are applied in the analysis.The geometry is immersed in a cubic B-spline block to perform FFD-based shape optimization, as shown in Figure 8b, where the crosssection of the optimal shape is indicated by the red curve.In this example, the horizontal and vertical coordinates of the FFD control points are design variables, totaling 200 design variables.Similar constraints are employed in this problem as those in the arch shape optimization.Control points on the left and bottom layers are fixed, where the lines in the FFD control net that are parallel to the tube axis maintain their orientations.
The optimization problem converges successfully after 42 iterations using SLSQP optimizer with a tolerance of 10 −12 , and snapshots are depicted in Figure 9.As expected, the initial square tube converges to the cylindrical tube.

T-beam shape optimization
In this section, a T-beam is considered to test the method for shell structures with intersections in the middle.We model a T-beam using two NURBS patches which are shown in Figure 10a.In the baseline design, the vertical patch in the T-beam is   In this benchmark problem, we aim to minimize the internal energy of the T-beam by updating the horizontal coordinates of shell patches' control points.Subsequently, the T-beam is placed in a cubic FFD block, where the horizontal coordinates of the FFD block's control points serve as design variables.The rear end of the T-beam is fixed.Given a downward distributed load, the optimal shape of the T-beam is expected to have a junction at the center of the horizontal patch under a constant volume constraint.The optimal position of the vertical patch is depicted in Figure 10b with red lines.The right and left layers of the FFD block's control points are fixed by employing equality constraints, and control points stay collinear in the axial direction.Using a tolerance of 10 −15 , the SLSQP optimizer converges successfully after 16 iterations, and the optimization process is shown in Figure 11.
In Figure 11, the T-beam converges to the optimal solution, indicating that the proposed method is effective for non-matching shell structures with intersections.A more complicated demonstration is presented in Section 6.2, where the geometry contains curved T-junctions.

Thickness optimization of a clamped plate
As stated in Section 4.3, the FFD-based optimization methodology can also be applied to thickness optimization.In the following section, we first demonstrate a piecewise constant thickness optimization for a clamped non-matching plate, in which the FFD block is not needed.Subsequently, we proceed to perform the variable thickness optimization for the same geometry.

Piecewise constant thickness optimization
For the thickness optimization example, a unit square plate composed of six cubic non-matching NURBS surfaces is considered.The geometry, which is shown in Figure 12a, exhibits 5 intersections with a total of 1449 DoFs.We apply a clamped boundary condition on the left side and with a line force applied to the right side in the normal direction of the plate.All patches of the plate have an initial thickness of 0.01 m.Using the strategy introduced in Section 4.3, we perform piecewise constant thickness optimization for the clamped plate to minimize the internal energy under the constant volume constraint.This problem only has 6 design variables.In this and the following sections, the SNOPT optimizer is used for faster convergence.The optimal thickness is plotted in Figure 12b.
The observed optimal piecewise constant thickness in Figure 12b shows material redistributes toward the clamped side, which provides enhanced support to the plate.The internal energy in the final configuration is 37.17% less than the baseline configuration.

Variable thickness optimization
We now perform variable thickness optimization for the non-matching plate.The plate is placed in a cubic B-spline block, which is shown in Figure 13a.Besides the constant volume constraint, we only optimize the plate thickness in one direction that is perpendicular to the intersections.
With the FFD-based approach, the continuity of shell thickness is maintained at all intersections.Figure 13b depicts the converged solution, where a smooth thickness distribution is observed.The smooth thickness profile offers an improved design compared to the piecewise constant thickness approach, resulting in a 40.20% reduction of the initial internal energy.A comparison of optimization iteration of normalized internal energy between the two methods is illustrated in Figure 14a.The FFD-based thickness optimization approach converges to a smaller internal energy.
To validate the proposed method, we compare the continuous thickness profile to an optimal thickness configuration of a cantilever beam [57].The cantilever beam is modeled using the Euler-Bernoulli beam theory, where a point load is applied to the free end.Since the Kirchhoff-Love shell is an extension of the Euler-Bernoulli beam, both models are expected to yield identical thickness distributions.The normalized thickness profiles of these two models, along with the piecewise constant thickness profile, are plotted in Figure 14b.A good agreement is observed between the variable thickness of the Kirchhoff-Love shell at the center line and the Euler-Bernoulli beam.Meanwhile, the cross-sectional view of the piecewise constant thickness shows a similar trend to the Euler-Bernoulli beam, albeit with discontinuities at the intersections.
We then investigate the effect of basis function order of continuity in the FFD block.Using the same knots vectors as illustrated in Figure 13a, we increase the order of the B-spline basis functions from linear (C 0 ) to quartic (C 3 ) and compare the amounts of reduced internal energy relative to the baseline configuration.These data points are summarized in Table 2.The results presented in Table 2 indicate that an FFD block with quadratic B-spline basis functions can achieve a better optimal thickness distribution for the clamped plate.The internal energy with quadratic FFD block only exhibits 0.27% of relative difference compared to the quartic FFD block.Table 2 also suggests that elevating the order of continuity of the FFD block leads to better designs with lower internal energy, particularly when transitioning from linear to quadratic B-spline basis functions.The plate example demonstrates that both piecewise constant and variable thickness optimization can be conducted in the proposed framework.One can select desired thickness distribution, or a mixed approach of these two, demonstrated in Section 6, based on the physical conditions and problem requirements to obtain an optimal design.

Application to aircraft wings
Two aircraft wing design optimization problems are performed.Section 6.1 considers a Parallel Electric-Gas Architecture with Synergistic Utilization Scheme (PEGASUS) wing thickness optimization problem, and Section 6.2 demonstrates a simultaneous shape and thickness optimization for an Electric Vertical Takeoff and Landing (eVTOL) aircraft wing.Both applications illustrate that a design-analysisoptimization workflow can be achieved with the proposed framework.

PEGASUS wing thickness optimization
For the PEGASUS wing problem, we first verify the accuracy of the structural analysis using PENGoLINS for a shell structure with a large number of patches and intersections.Then two types of thickness design optimizations are performed in the following section.

Structural analysis of the PEGASUS wing
The PEGASUS wing CAD model is created using a customized geometry engine, an exploded view of the wing with IGA discretization is shown in Figure 15  The PEGASUS wing is made of material with Young's modulus 69 GPa and Poisson's ratio 0.35, and the wing span is 12.22 m.At the wing root, the chord is 1.52 m and the airfoil thickness is 0.37 m.A uniform initial thickness is 5 mm for all patches.Considering an aircraft take-off weight of 9000 kg, a distributed upward pressure with a magnitude of 132.5 N/m 2 is determined by dividing half of the take-off weight by the surface area of the wing.Clamped boundary conditions are imposed at the wing root.Importing the PEGASUS wing geometry in IGES format to PENGoLINS, we perform structural analysis directly without finite element mesh generation.Given an analysis-suitable CAD file, the only required geometry preprocessing is to approximate surface-surface intersections, which is a much easier effort than finite element mesh generation and is automated in PENGoLINS using the functionality provided by pythonOCC.Figure 16a shows all the intersections presented in the PEGASUS wing, while the displacements computed by PENGoLINS are visualized in Figure 16b.To validate the calculated displacements, we conduct FE analysis for the PEGA-SUS wing using COMSOL.Figure 16c displays an extensively refined finite element mesh for the COMSOL FE analysis.Displacements solved in COMSOL, utilizing the Reissner-Mindlin shell theory, are depicted in Figure 16d.In this analysis, quadratic triangular elements with 118644 DoFs are used.Figure 16 indicates that the displacements obtained from PENGoLINS closely match the results from COMSOL.The maximum displacement magnitude in PENGoLINS is 0.06662 m, which has a relative difference of 0.15% compared to the corresponding value of 0.06672 m in COMSOL.This aligns well with the findings of a numerical experiment presented in [44,Section 5.3], where a relative difference of 0.76% is observed for the converged vertical displacement at the wingtip of an eVTOL wing between PENGoLINS and an open-source Reissner-Mindlin shell solver [58] using formulations from [59].The simulation results for the PEGASUS wing indicate that PENGoLINS provides good accuracy for complex shell structures, which is crucial for the subsequent design optimization.

Thickness optimization of the PEGASUS wing
Similar to the thickness optimization of the clamped plate discussed in Section 5.4, we apply the same methodology to the PEGASUS wing for piecewise constant thickness optimization.The same boundary conditions are employed throughout the optimization.In the piecewise constant thickness optimization case, a total of 90 design variables are included with lower and upper bounds of 1 mm and 100 mm, respectively.The initial thickness for all patches is taken as 5 mm.A constant volume constraint is employed in the optimization.The resulting shell thickness with minimum internal energy is depicted in Figure 17.The shell patch with the maximum thickness is observed at the wing root in the outer skins, while the thickness decreases consistently along the span direction for both the outer skins and spars, following a pattern similar to that of the clamped plate.The thicknesses of the internal ribs and the wingtip are close to the lower bound since the bending moments are mainly carried by the lower and upper skins given the distributed upward load.Therefore, the majority of material is redistributed towards the clamped root of the outer skins.The optimized design in Figure 17 gives an internal energy 38.17% less than that of the baseline configuration.To achieve an improved design, we consider variable thickness in the outer skins and spars of the PEGASUS wing, while keeping the internal ribs with piecewise constant thickness.The configuration of the FFD blocks for variable thickness optimization is illustrated in Figure 18a.Four FFD blocks with quadratic bases are created to allow for variation in thickness within a single spline patch while ensuring continuity at patch intersections, with a total of 402 design variables used for this problem.By minimizing the internal energy again, the optimal thickness distribution is obtained as shown in Figure 18b.Both applications in this section use the SNOPT optimizer with a tolerance of 10 −4 .In the optimized design, the outer skins of the PEGASUS wing at the clamped edge still have the largest thickness, where the thickness distribution is smooth at the surface-surface intersections within each FFD block.The decreasing thickness trend in outer skins and spars remains consistent with the optimal piecewise constant thickness case.Comparing the combined optimization strategy with the piecewise constant method, the maximum thickness in the former is higher while maintaining the same volume.Additionally, the internal energy is reduced by 44.71% compared to the baseline design.These observations indicate that the combined thickness optimization method demonstrates a more efficient utilization of material than the piecewise constant method.

Simultaneous optimization for eVTOL wing
With the continuous advancements in aviation battery technology [60], eVTOL aircraft have emerged as a promising solution for cost-effective urban mobility [61].
In this section, we use a more advanced eVTOL wing to demonstrate the capabilities of the FFD-based optimization approach, where both the thickness and shape control points are considered as design variables simultaneously.By incorporating both thickness and shape coordinates in the design of shell structures, we can utilize the material more efficiently than considering thickness optimization only.The CAD geometry of the eVTOL wing, material parameters, and the corresponding structural analysis can be found in [44, Section 5].For the optimization problem, we implement the same clamped boundary conditions and distributed upward pressure as those in the previous application.The magnitude of pressure is set to 120 N/m 2 , approximated by dividing the take-off weight by the surface area of the eVTOL wing.The baseline configuration of the eVTOL wing, which is composed of 27 cubic NURBS patches with 87 intersections, is illustrated in Figure 19.We note that we perform the shape optimization for the eVTOL wing without including an aerodynamic model (which would be needed for a well-posed wing design problem) purely to provide a demonstration of the FFD-based method for complex aerospace structures.To achieve a meaningful design for the eVTOL wing, we use two quadratic Bspline FFD blocks for thickness optimization.This configuration allows for varying thicknesses in the lower and upper skins of the eVTOL, while using piecewise constant thicknesses for the internal stiffeners and wingtips.The arrangement of the thickness FFD blocks is illustrated in Figure 19a.Furthermore, a cubic B-spline FFD block is created for shape optimization, as depicted in Figure 19b.Only the vertical coordinates of control points for the shape FFD block, denoted as P FFD3 , are updated.In total, there are 642 design variables involved in this optimization process, and a constant volume constraint is applied as well.Regarding the thickness design variables, the lower and upper limits are selected as 1 mm and 50 mm, respectively.All shell components have initial thicknesses of 3 mm.
One challenge encountered during shape optimization of complex shell structures, such as eVTOL wings, is the potential occurrence of oscillatory or highly distorted geometries in the updated designs.These distorted shapes can lead to poor element quality, resulting in spurious energy and affecting the accuracy of the analysis.To mitigate oscillation or radical change of shell components, an additional term is introduced in the objective function to regularize the gradient of the shape.The objective is formulated as where I is the index of shell patches, λ is a dimensionless regularization coefficient, h I A is the element area of shell patch I in the physical space, P I 3 is the vertical component of the shape variable for shell patch I, and S I is the midsurface of shell patch I. ( •) denotes quantities in the baseline configuration.The regularization term in (49) serves as an additional artificial bending energy associated with the curvature of the eVTOL wing.Therefore, the shape oscillation can be eliminated by adjusting the regularization coefficient λ.The SNOPT optimizer is employed with a tolerance of 10 −3 .Figure 20 demonstrates the optimized eVTOL wing designs achieved with varying λ values, providing insights into the influence of regularization on the final shape and thickness distribution.It can be seen that the patch intersections of the optimized shapes in Figure 20 are still connected using the FFD-based approach.The reductions of internal energy for different λ are listed in Table 3.For this eVTOL wing optimization problem, different regularization coefficients yield distinct outcomes.A regularization coefficient of 1 and 0.1 provide optimal designs dominated by thickness update, resembling the patterns observed in the combined thickness optimization of the PEGASUS wing, refer to Figure 18b.The optimized shapes in these two cases are similar to the baseline configuration since the artificial bending energies are the major contributor to the objective function due to the large regularization coefficients.Even slight variations in the shape variables result in substantial increases in the objective function.Conversely, small regularization coefficients, i.e., 10 −4 and 10 −5 , lead to considerable changes in the shape of the eVTOL wing and reduction of internal energy, amounting to 92.82% and 95.09%, respectively.However, these cases exhibit noticeable oscillations in the wingtip area, which can lead to ill-conditioned systems.On the other hand, employing regularization coefficients with values of 10 −2 and 10 −3 yields balanced solutions, where both the thickness redistribution and shape updates contribute to the optimal design.The material moves towards the clamped area, accompanied by a widening of the crosssection to provide increased wing support.An exploded view of the optimal design with λ = 10 −3 is shown in Figure 21.This optimal design achieves an internal energy reduction of 83.23% compared to the baseline configuration.

Conclusion
In this paper, we have introduced an FFD-based shape and thickness optimization approach for shell structures composed of separately-parametrized NURBS surfaces.The integration of this method with the Lagrange extraction technique enables IGA with existing FE toolkits and provides a connection between the FFD block and nonmatching shell patches.By employing the FFD block approach, the updated shell geometry and thickness remain properly connected at patch intersections throughout the optimization process.This feature prevents undesired shape discontinuities in shell structures.We made use of the penalty-based coupling formulation for Kirchhoff-Love shells to perform IGA in the optimization.The automation of analytical derivative computations is achieved through code generation in FEniCS, enabling gradientbased multidisciplinary design optimization.This automation streamlines the optimization process and allows for efficient exploration of design spaces.The unified NURBS representation shared by both the design geometry and analysis model enhances accuracy per DoF in the analysis and precise design updates.Moreover, the proposed framework circumvents FE mesh generation and streamlines design-analysis-optimization workflow for complex shell structures.Consequently, the automated workflow is expected to accelerate the conceptual design of novel eVTOL aircraft with minimal manual effort.
A suite of benchmark problems is adopted to verify the effectiveness of the FFD-based optimization approach proposed in this paper.Both the shape optimization and thickness optimization results agree well with analytical solutions or other established references.Furthermore, we have applied the framework to two different aircraft wings.This demonstration highlights the potential of the proposed method in exploring complex design spaces and obtaining superior designs for innovative aircraft structures.Future works can extend the shape optimization for aircraft wings with more realistic loads by coupling the structural solver with an aerodynamic solver [62].Additionally, practical constraints such as von Mises stress can be incorporated into the optimization process to ensure that the obtained designs meet the requirements of real-world applications.To promote code transparency and potential contributions to the shell optimization community, we make the Python library GOLDFISH openly accessible in a GitHub repository [56].

Fig. 3 :
Fig.3: Workflow of FFD-based shape optimization for non-matching shell structures.A cylindrical roof consisting of four non-matching NURBS patches is first immersed in a trivariate B-spline block.We update the control points of the FFD block to deform the shape of the non-matching cylindrical roof.Control nets of the NURBS surfaces are indicated with red color, and black is used for the control net of the FFD block.

Fig. 4 :
Fig. 4: (a) Initial configuration of the cylindrical roof geometry consisting of four non-matching NURBS patches.(b) Updated NURBS surfaces using FFD block.

Fig. 5 :
Fig. 5: Sliced view of the intersecting edges between shell patches S C and S D of the cylindrical roof.The two edges remain overlapping in the updated configuration.

Fig. 6 :
Fig. 6: (a) Baseline geometry of a non-matching arch consisting of four NRUBS patches, three surface-surface intersections are indicated with red lines.(b) Initial configuration of the arch immersed in an FFD block, where black lines and dots denote the control net.The analytical optimal shape is plotted with a red curve.

Fig. 8 :
Fig. 8: (a) Baseline geometry of the square tube, a quarter of the tube is modeled using four non-matching B-spline patches with four intersections.(b) Initial configuration of the FFD block with control net.The optimal cross-section is depicted by a red circle.

Fig. 10 :
Fig. 10: (a) Baseline configuration of a T-beam whose vertical patch is at the threequarter position of the horizontal patch.(b) The T-beam is placed in an FFD B-spline block.The optimal position of the vertical patch is depicted with red lines.

Fig. 12 :Fig. 13 :
Fig. 12: (a) A unit square plate consisting of six non-matching patches, intersections are indicated with red lines.(b) Final plate thickness for piecewise constant thickness optimization.

Fig. 14 :
Fig. 14: (a) Optimization process of normalized internal energy for two approaches.(b) Cross-sectional view of piecewise constant thickness and variable thickness, and comparison with Euler-Bernoulli beam thickness optimization.
. The CAD model consists of 90 NURBS patches, two outer skins (one lower skin and one upper skin) and two spars (one front spar and one rear spar) connecting two adjacent ribs.The NURBS surfaces in the PEGASUS wing are represented using cubic basis functions with maximal continuity, resulting in 19572 DoFs in total.Additionally, 280 patch intersections are formed in the wing structure.

Fig. 15 :
Fig. 15: CAD geometry of the PEGASUS wing which is composed of 90 NURBS patches with 280 intersections, totaling 19572 DoFs.
(a) Illustration of surface-surface intersections in the PEGASUS wing geometry.(b) Visualization of displacement magnitude solved by PENGoLINS in baseline design.(c) Finite element mesh of PEGASUS wing generated in COMSOL.(d) Analysis result obtained from COMSOL using Reissner-Mindlin shell theory.

Fig. 16 :
Fig. 16: Structural analysis of the PEGASUS wing using PENGoLINS, and the resulting displacement magnitude is compared with the corresponding output from COMOSL.

Fig. 17 :
Fig. 17: Optimization result of the PEGASUS wing with piecewise constant thickness.

Fig. 18 :
Fig. 18: (a) Configuration of the combined thickness optimization.Each group of outer skins and spars is placed in one FFD block, and the remaining internal ribs have a piecewise constant thickness.(b) Optimal thickness distribution of PEGASUS wing.

Fig. 19 :
Fig. 19: (a) FFD blocks for eVTOL wing thickness optimization, where the lower and upper skins have a variable thickness.Wingtips and internal ribs and spars have a piecewise constant thickness.(b) FFD block for shape optimization.

Fig. 21 :
Fig. 21: Exploded view of optimal design of eVTOL wing with regularization coefficient λ = 10 −3 , which results in 83.23% reduction of internal energy compared to the baseline design.
Move quadrature mesh Ω M to parametric coordinate ξ MB , which is depicted in the lower left configuration in Figure2.Then create transfer matrices T 0MB and T 1MB. 8. Repeat steps 3, 4 and 6 to obtain corresponding displacement u 0MB , normal vectors A BM 3 and a BM 3 , and element size h MB X from S B . 9. Calculate penalty parameters based on (10) using averaged physical element size h

Table 2 :
Reduction of internal energy of the clamped plate for different degrees of the FFD block.

Table 3 :
Reduction of internal energy of the eVTOL wing after simultaneous optimization with varying regularization coefficients.