Isogeometric Configuration Design Optimization of Three-dimensional Curved Beam Structures for Maximal Fundamental Frequency
NNoname manuscript No. (will be inserted by the editor)
Isogeometric Configuration Design Optimization ofThree-dimensional Curved Beam Structures for MaximalFundamental Frequency
Myung-Jin Choi · Jae-Hyun Kim · Bonyong Koo · Seonho Cho Published in
Structural and Multidisciplinary Optimization , DOI: https://doi.org/10.1007/s00158-020-02803-0Received: / Accepted: date
Abstract
This paper presents a configuration designoptimization method for three-dimensional curved beambuilt-up structures having maximized fundamental eigen-frequency. We develop the method of computation ofdesign velocity field and optimal design of beam struc-tures constrained on a curved surface, where both de-signs of the embedded beams and the curved surface aresimultaneously varied during the optimal design pro-cess. A shear-deformable beam model is used in theresponse analyses of structural vibrations within an iso-geometric framework using the NURBS basis functions.An analytical design sensitivity expression of repeatedeigenvalues is derived. The developed method is demon-strated through several illustrative examples.
Keywords
Configuration design · Beam structure · Fundamental frequency · Repeated Eigenvalues · Design velocity field · Isogeometric analysis
A structural design to have eigenfrequencies as far awayas possible from excitation frequencies has been re-garded as a significant process to avoid resonance and
Corresponding author (S. Cho)Tel.: +82 2 880 7322Fax: +82 2 888 9298E-mail: [email protected] Department of Naval Architecture and Ocean Engineering,Seoul National University, Seoul, Korea School of Mechanical Convergence System Engineering,Kunsan National University, Gunsan, KoreaThis document is the personal version of an article whosefinal publication is available at https://doi.org/10.1007/s00158-020-02803-0. assure structural stability. There have been lots of re-searches to develop effective and efficient design meth-ods to increase the structural fundamental frequency.Shi et al. [22] presented a parameter-free approach fora design optimization of orthotropic shells that maxi-mizes fundamental frequencies, where a shape variationwas generated by a displacement field due to an exter-nal force equivalent to negative shape gradient function.Hu and Tsai [11] found optimal fiber orientation anglesto maximize fundamental frequencies of fiber-reinforcedlaminated cylindrical shells. Blom et al. [2] optimizedthe fiber orientation on the surface of a conical shellfor the maximization of the fundamental frequency un-der manufacturing constraints on the curvature of fiberpaths. Silva and Nicoletti [23] sought an optimal shapedesign of a beam for desired natural frequencies by alinear combination of mode shapes whose coefficientsare determined by an optimization algorithm. Wang etal. [27] derived an analytic solution of minimum sup-port stiffness for maximizing a natural frequency of abeam. Zhu and Zhang [13] proposed an optimal sup-port layout for maximizing natural frequency using atopology optimization where stiffness values of supportsprings are considered as design variables. Allaire andJouve [1] utilized a level-set method to design two-and three-dimensional structures having maximal fun-damental frequencies. Picelli et al. [19] applied an evolu-tionary topology optimization method to maximize thefirst natural frequency in free-vibration problems con-sidering acoustic-structure interaction. Zuo et al. [28]employed a reanalysis method for eigenvalue problemto reduce computation cost during a design optimiza-tion by a genetic algorithm for minimizing the weight ofstructures with frequency constraints. In the design sen-sitivity analysis (DSA) of eigenvalue problems, a tech- a r X i v : . [ phy s i c s . a pp - ph ] J a n M.-J. Choi et al. nical difficulty may arise due to the lack of the usualdifferentiability in repeated eigenvalues. Seyranian etal. [21] presented a perturbation method to derive de-sign sensitivities of repeated eigenvalues, where it wasshown that sensitivities of repeated eigenvalues corre-spond to solutions of a sub-eigenvalue problem. Thismethod was used by Du and Olhoff [7] for topology op-timizations in order to maximize specific eigenfrequen-cies and gaps between two consecutive eigenfrequen-cies. It was shown, in Olhoff et al. [18], that the sep-aration of adjacent eigenfrequencies is closely relatedwith the generation of band gap property.
In this pa-per, we present a gradient-based configuration designoptimization method for three-dimensional beam built-up structures having maximal fundamental frequencies,using an analytical configuration DSA for repeated aswell as simple eigenvalues.
Isogeometric analysis (IGA) suggested by Hugheset al. [12] fundamentally pursues a seamless integrationof computer-aided design (CAD) and analysis. Espe-cially, within the IGA framework, shape design, andsizing design parameters like beam cross-section thick-ness and composite fiber orientation angle have higherregularity due to NURBS basis functions than the con-ventional finite element analysis (FEA) model. Taheriand Hassani [26] combined shape and fiber orienta-tion design optimizations for functionally graded struc-tures to have optimal eigenfrequencies. Nagy et al. [17]performed a sizing and configuration design optimiza-tion for a fundamental frequency maximization of pla-nar Euler beams. They also presented several types ofmanufacturing constraints to avoid the abrupt changeof cross-section thickness and neutral axis curvature.Liu et al. [15] presented a sizing optimal design of pla-nar Timoshenko beams, where the maximum rate ofchange of cross-section thickness was constrained, andthe constraints were aggregated by the K-S constraintscheme. These previous works on the isogeometric opti-mal structural design for maximizing the fundamentalfrequency using a beam element were limited to thedesign of two-dimensional structure and single patchmodel.
In this paper, we deal with configuration designof more practical three-dimensional built-up structures.
In this work, we particularly focus on the compu-tation of design velocity field and optimal design ofbeam structures constrained in curved surfaces. Struc-tures geometrically constrained on curved domains havebeen widely utilized in many engineering applicationsincluding substructures for offshore oil and wind tur-bine platforms. This kind of design problem has equal-ity constraints restricting their position on a specifiedcurved domain, and several researches regarding theconstrained optimal design problem have been reported recently. Liu et al. [14] performed topology optimiza-tion by combining the moving morphable component(MMC) method and coordinate transformations to ar-chitect structures constrained in planar curved domains.Choi and Cho [3] presented a constrained isogeomet-ric design optimization method for beam structures lo-cated on a specified curved surface, where the free-formdeformation (FFD) method was used to locate curvesexactly on the surface, and then the curves are inter-polated to determine control point positions. Choi andCho [4] presented the construction of initial orthonor-mal base vectors along a spatial curve embedded in acurved domain by using surface convected frames asreference frames in the smallest rotation (SR) method.However, the previous works [3] and [4] assumed thatthe surfaces where the curves are embedded do nothave design dependence, that is, beam structures aredesigned on a fixed curved surface.
In this paper, ex-tending those previous works, we additionally considerconfiguration design of the curved surfaces where beamstructures are embedded, so that more design degrees-of-freedom are furnished during the constrained optimaldesign process.
This article is organized as follows. In section 2,we explain an isogeometric analysis of structural vi-brations in shear-deformable beam structures. In sec-tion 3, we explain a configuration DSA of simple andrepeated eigenvalues for the beam structures, wherea design velocity computation method for constrainedbeam structures is detailed. In section 4, the developedoptimization method is demonstrated through severalnumerical examples for maximizations of fundamentalfrequencies. ξ = { ξ , ξ , · · · , ξ n + p +1 } , (1)where p and n are respectively the order of basis func-tion, and the number of control points. The B-splinebasis functions are obtained, in a recursive manner, as N I ( ξ ) = (cid:26) ξ I ≤ ξ < ξ I +1 , ( p = 0) , (2) onfiguration Design Optimization for Maximal Fundamental Frequency 3 and N pI ( ξ ) = ξ − ξ I ξ I + p − ξ I N p − I ( ξ ) + ξ I + p +1 − ξξ I + p +1 − ξ I +1 N p − I +1 ( ξ ) , ( p = 1 , , , . . . ) . (3)From the B-spline basis N pI ( ξ ) and the weight w I , theNURBS basis function W I ( ξ ) is defined by W I ( ξ ) ≡ N pI ( ξ ) w In (cid:80) J =1 N pJ ( ξ ) w J . (4)The NURBS basis function of order p has C p − k con-tinuity at a knot multiplicity k . A curve geometry isdescribed by a linear combination of NURBS basis func-tions and control point positions as C ( ξ ) = n (cid:88) I =1 W I ( ξ ) B I . (5)More details of the NURBS geometry can be found in[20].2.2 Linear kinematics and equilibrium equationsConsider an initial beam neutral axis ϕ ( ξ ) ∈ R char-acterized by a NURBS curve of Eq. (5). We define anarc-length parameter s ∈ Ω ≡ [0 , L ] that parameter-izes the beam neutral axis with the initial (undeformed)length L in a way that [6] s ( ϕ ( ξ )) = (cid:90) ξ = ξ (cid:13)(cid:13)(cid:13)(cid:13) ∂ ϕ ( ξ ) ∂ξ (cid:13)(cid:13)(cid:13)(cid:13) dξ ≡ s ( ξ ) , (6)which relates parametric domains of the arc-length pa-rameter ( Ω ) and NURBS parametric coordinate ( Ξ )through the physical domain Ω , as depicted in Fig. 1.The Jacobian of the mapping is defined as [6] J c ≡ ∂s∂ξ = (cid:13)(cid:13)(cid:13)(cid:13) ∂ ϕ ( ξ ) ∂ξ (cid:13)(cid:13)(cid:13)(cid:13) . (7)Hereafter, we often denote ϕ ≡ ϕ ( s ) for a brief ex-pression. A cross-section orientation is described by anorthonormal frame with basis j ( s ) , j ( s ) , j ( s ) ∈ R ,where the unit tangent vector is obtained by j ( s ) ≡ ϕ ,s due to the arc-length parameterization, and as-sumed to be orthogonal to the cross-section (see Fig.2 for an illustration). ( • ) ,s denotes the partial differen-tiation with respect to the arc-length parameter. Theglobal Cartesian coordinate system is denoted by X − Y − Z with base vectors e , e , e ∈ R . In this pa-per, we consider a structural dynamics with small (in-finitesimal) amplitudes, where a linearized kinematics Fig. 1: Parameterization of a beam neutral axis [4,6]Fig. 2: Problem definition of a spatial rod: neutral axisdefined by general curve ϕ ( s )consistently derived from a general nonlinear formula-tion of [24] is utilized. The (material form) axial-shearand bending-torsional strain measures derived in [25]are evaluated at an undeformed state, which respec-tively result in [5] Γ ( s, t ) ≡ Λ ( s ) T (cid:26) ∂∂s z t + j ( s ) × θ t (cid:27) , (8)and Ω ( s, t ) ≡ Λ ( s ) T ∂ θ t ∂s . (9) t ≥ z t ≡ z ( s, t ) denotesthe displacement of the neutral axis position, and θ t ≡ θ ( s, t ) denotes the (infinitesimal) rotation vector thatdescribes the rigid rotation of cross-section. Λ ( s ) ≡ [ j ( s ) , j ( s ) , j ( s )] defines an orthogonal transformationmatrix such that j I ( s ) = Λ ( s ) e I , I = 1 , , . (10)We define an elastic strain energy by [5] U ≡ (cid:90) Ω { Γ ( s, t ) · N ( s, t ) + Ω ( s, t ) · M ( s, t ) } ds, (11)where N ( s, t ) and M ( s, t ) are the material form resul-tant force and moment vectors over the cross-section at s ∈ Ω , which are related to corresponding spatial formsthrough n ( s, t ) = Λ ( s ) N ( s, t ) m ( s, t ) = Λ ( s ) M ( s, t ) (cid:27) . (12) M.-J. Choi et al.
As shown in Fig. 2, the (spatial form) resultant forceand moment can be respectively resolved in the or-thonormal basis as n ( s, t ) = n I ( s, t ) j I ( s ) and m ( s, t ) = m I ( s, t ) j I ( s ) where the repeated index I implies thesummation over the components 1, 2, and 3. We em-ploy a linear elastic constitutive relation as N ( s, t ) = C F Γ ( s, t ) M ( s, t ) = C M Ω ( s, t ) (cid:27) , (13)and the constitutive tensors C F and C M are definedby C F = diag [ EA, GA , GA ] { e I } C M = diag [ GI p , EI , EI ] { e I } (cid:27) , (14)where diag [ a, b, c ] represents a diagonal matrix of com-ponents a , b , and c . ( • ) { e I } represents the tensor is rep-resented in the global basis { e , e , e } . A representsthe cross-sectional area, and A ≡ k A , A ≡ k A where k and k denote the shear correction factors. I p denotes the polar moment of inertia, and I and I mean the second moments of inertia. E and G denoteYoung’s modulus and shear modulus, respectively. Akinetic energy is defined by [5] T ≡ (cid:90) Ω { ρA z t,t · z t,t + θ t,t · ( I ρ θ t,t ) } ds, (15)where ( • ) ,t denotes the partial differentiation with re-spect to the time, and ρ denotes the mass density. It isassumed that the initial local orthonormal basis { j , j , j } is directed along the principal axes of inertia of thecross-section; thus, the inertia tensor can be expressedby a diagonal matrix in the local basis, and representedin the global basis through the orthogonal transforma-tion of Eq. (10), as I ρ ≡ diag [ ρI p , ρI , ρI ] { j I } = Λ diag [ ρI p , ρI , ρI ] { e I } Λ T , (16)where ( • ) { j I } represents the corresponding tensor is rep-resented in the bases of { j , j , j } . A work done by ex-ternal loads is expressed as [5] W ≡ (cid:90) Ω { n ext ( s ) · z ( s, t ) + m ext ( s ) · θ ( s, t ) } ds, (17)where n ext ( s ) and m ext ( s ) denote the distributed ex-ternal force and moment, respectively. In order to de-rive equilibrium equations from time t to time t , theHamilton’s principle is employed as [8] δ (cid:90) t t ( U − T − W ) dt = 0 , (18)where δ ( • ) defines the first variation. Substituting Eqs.(11), (15), and (17) into Eq. (18), and applying the integration by parts lead to the following linear andangular momentum balance equations [5] n t,s + n ext = ρA z t,tt m t,s + j × n t + m ext = I ρ θ t,tt (cid:27) , (19)where n t ≡ n ( s, t ), m t ≡ m ( s, t ), and ( • ) ,tt denotesthe second-order partial derivative with respect to thetime.2.3 Variational formulationWe assume a time-harmonic solution such as z ( s, t ) = z ( s ) e − iωt θ ( s, t ) = θ ( s ) e − iωt (cid:27) , (20)where i ≡ √−
1, and ω denotes an angular frequency.From the governing equations of Eq. (19), using theprinciple of virtual work and substituting Eq. (20), weobtain the following generalized eigenvalue problem. a ( η , ¯ η ) = ζd ( η , ¯ η ) , ∀ ¯ η ∈ ¯ Z, (21)where ζ ≡ ω denotes the eigenvalue associated withthe eigenfunction η ≡ ( z ( s ) , θ ( s )). ( • ) ≡ δ ( • ) denotesthe first variation of ( • ) or the corresponding virtualquantity, and ¯ η ≡ ( ¯z ( s ) , ¯ θ ( s )). The strain energy andkinetic energy bilinear forms are respectively defined as[5] a ( η , ¯ η ) ≡ (cid:90) Ω (cid:26) ¯Γ¯Ω (cid:27) · (cid:26) NM (cid:27) ds, (22)and d ( η , ¯ η ) ≡ (cid:90) Ω (cid:26) ¯z¯ θ (cid:27) T (cid:20) ρA I 0 × × I ρ (cid:21) (cid:26) z θ (cid:27) ds, (23)where × ∈ R × and I ∈ R × represent the nulland identity matrices, respectively. It is assumed through-out this paper that the eigenfunctions are orthonormal-ized with respect to the bilinear form d ( • , • ) such that d ( η , η ) = 1 . (24)The solution space ¯ Z is expressed, considering a homo-geneous boundary condition, as¯ Z = (cid:8) ¯ η ∈ H × H : ¯z (0) = ¯ θ (0) = ¯z ( L ) = ¯ θ ( L ) = (cid:9) , (25)where H defines the Sobolev space of order 1. Thestatic (material form) axial-shear and bending-torsional onfiguration Design Optimization for Maximal Fundamental Frequency 5 strain measures are respectively obtained from Eqs. (8)and (9), as Γ ( η ) ≡ Λ T ( z ,s + j × θ ) Ω ( θ ) ≡ Λ T θ ,s (cid:27) , (26)and ¯Γ ≡ Γ ( ¯ η ) and ¯Ω ≡ Ω ( ¯ θ ) represent the virtualstrain measures. Eq. (26) can be rewritten in a compactform as (cid:26) Γ ( η ) Ω ( θ ) (cid:27) = Π T Ξ T (cid:26) z θ (cid:27) , (27)using the following matrix operators Π ≡ (cid:20) Λ × × Λ (cid:21) and Ξ ≡ (cid:20) dds I 0 × − [ j × ] dds I (cid:21) , (28)where [( • ) × ] denotes a skew-symmetric matrix associ-ated with the dual vector ( • ) ∈ R .2.4 Isogeometric discretization using the NURBS basisIn the IGA framework, the response field is approxi-mated by the same basis functions of NURBS as em-ployed to represent the geometry in CAD. The displace-ment and rotation amplitudes are approximated by theNURBS basis function, as z h = n (cid:88) N =1 W N ( ξ ) y N and θ h = n (cid:88) N =1 W N ( ξ ) θ N , (29)where y N and θ N are the coefficients associated withthe N -th control point. In the same manner, the vari-ation of the displacement and rotation amplitudes areapproximated by ¯z h = n (cid:88) N =1 W N ( ξ ) ¯y N and ¯ θ h = n (cid:88) N =1 W N ( ξ ) ¯ θ N , (30)where ¯y N and ¯ θ N are the N -th coefficients. Substitut-ing Eq. (29) into Eq. (27), we have (cid:26) Γ ( η h ) Ω ( θ h ) (cid:27) = Π T Ξ h N d N (31)where d N ≡ [ y N T , θ N T ] T , and we use Ξ h N ≡ (cid:20) W N,s
I 0 × − W N [ j × ] W N,s I (cid:21) . (32) W N,s denotes the differentiation of NURBS basis func-tion with respect to the arc-length coordinate, derivedas W N,s ≡ W I,ξ /J c [6]. ( • ) ,ξ denotes the partial dif-ferentiation with respect to the parametric coordinate ξ . Then the strain energy bilinear form of Eq. (22) isdiscretized by a ( η h , ¯ η h )= A ne e =1 (cid:20) ¯d TN (cid:26)(cid:90) Ξ e (cid:16) Ξ h N Π CΠ T Ξ h M T J c (cid:17) dξ (cid:27) d M (cid:21) ≡ ¯d T Kd , (33)where d N ≡ [ y N T , θ N T ] T and ¯d N ≡ [ ¯y TN , ¯ θ N T ] T . Therepeated indices N and M imply the summations overthe range from 1 to n , and n denotes the number ofcontrol points having local supports in the knot span Ξ e . ne denotes the number of knot spans. K , d , and ¯d are the assembled global stiffness matrix, the globalvectors of response coefficients, and virtual response co-efficients, respectively. Also, the kinetic energy form ofEq. (23) is discretized by d ( η h , ¯ η h )= A ne e =1 (cid:20) ¯d TN (cid:26)(cid:90) Ξ e (cid:20) ρA I 0 × × I ρ (cid:21) W N W M J c dξ (cid:27) d M (cid:21) ≡ ¯d T Md . (34)Thus, using Eqs. (33) and (34), we have the followingdiscretized forms of the generalized eigenvalue problemof Eq. (21). Kd = ζ Md , (35)where it is assumed that the eigenvectors are M -orthonormalized,from Eq. (24), as d T Md = 1 . (36)It is noted that a rotation continuity condition at theinterface between NURBS curve patches can be auto-matically satisfied by a typical assembly process, sincewe employ the rotation vectors as the control coeffi-cients in Eqs. (29) and (30).2.5 Geometric modeling of embedded structures:reviewIn the previous work of [3], geometric modeling of con-strained beam structure on a curved surface was con-ducted by combining a free-form deformation (FFD)and a global curve interpolation. Here, we briefly ex-plain the overall procedure. For example, we considerconstructing a beam structure on a cylindrical surface,as illustrated in Fig. 3. We first construct a structure, asshown in Fig. 3a, in a rectangular domain with width b and height b , so-called the “ initial ” parent domain ˜ Ω . M.-J. Choi et al.
The embedded structure in the initial parent domain isexpressed by a NURBS curve of Eq. (5), as X ( ξ ) ≡ n (cid:88) I =1 W I ( ξ ) B I ∈ ˜ Ω . (37) B I and n define a control point position and the numberof control points, respectively. We also define a “ target ”parent domain ( ˜ Ω ) like the cylindrical surface havingradius R and length L shown in Fig. 3b whose geometryis described by a NURBS surface as ˜X ( ˜ ξ , ˜ ξ ) = m (cid:88) J =1 ˜ W J ( ˜ ξ , ˜ ξ ) ˜B J , (38)where ˜ W J ( ˜ ξ , ˜ ξ ) and ˜ ξ i ( i = 1 ,
2) are the bivariateNURBS basis function and surface parametric coordi-nates, respectively, and ˜B J defines the control pointposition, and m denotes the number of control points.An FFD states that the embedded object in an initialparent domain is mapped onto a target parent domain,as shown in Fig. 3b in a way that [3] X ( ξ ) ≡ ˜X ( ˜ ξ ( ξ ) , ˜ ξ ( ξ ))= m (cid:88) J =1 ˜ W J ( ˜ ξ ( ξ ) , ˜ ξ ( ξ )) ˜B J , (39)where the surface parametric position ( ˜ ξ , ˜ ξ ) corre-sponding to the curve parametric position ξ is deter-mined by a simple algebra as˜ ξ i = X i ( ξ ) /b i , i = 1 , , (40)and X i denotes i -th component of X . It is noted thatthe embedded geometry X is exactly on the target par-ent domain; however, in order to determine a corre-sponding control net, a global curve interpolation needsto be additionally performed, which solves the followingsystem of linear equations [3]. N g B g = X g , (41)where an invertible matrix N g and the right-hand sidevector X g are respectively constructed by assemblingthe evaluations of the NURBS basis function W I ( ξ cL )and embedded curve position X ( ξ cL ) by Eq. (39) at se-lected parametric positions ξ cL ( L = 1 , ..., n ). From thesolution B g of Eq. (41), we extract a set of control pointpositions { B I } for each patch in the embedded curveson target parent domain (Fig. 3c), then a reconstructedgeometry, that is, a beam neutral axis geometry can beexpressed by X ∗ ( ξ ) = n (cid:88) I =1 W I ( ξ ) B I . (42) In [3], it was shown that a deviation of the reconstructedgeometry X ∗ ( ξ ) from the target one X ( ξ ), so-called a“ modeling error ”, is able to be reduced by increasingthe degrees of freedom (DOFs) of the reconstructedcurve through a knot insertion or a degree-elevationprocess. It was also investigated that the local non-smooth region of target parent domain can be effec-tively resolved by the multi-patch modeling of beamstructures. For more detailed discussions regarding theway to reduce the modeling error, one can refer to thereference [3].The configuration design DOFs of the reconstructedcurve, shown in Fig. 3c, can be divided into two parts:first, the geometry of embedded lattice structure in ini-tial parent domain (e.g., Fig. 3a), and second, the ge-ometry of target parent domain, for example, the cylin-drical surface shown in Fig. 3b may change its radius( R ) and length ( L ). In the previous work of [3], onlythe configuration design of embedded structures in theinitial parent domain was altered, that is, a structurewas constrained on a given specific surface. However,in this paper, we extend the previous work to considera design dependence of target parent domain, so that wehave additional design DOFs.
In the following section,we explain the design velocity computation consideringthe design dependence of target parent domain. V ( s ) as [6] ϕ τ ( s τ ) ≡ ϕ ( s ) + τ V ( s ) , (43)where τ ≥ τ represents a quantityevaluated at a perturbed design here and hereafter.Taking the material derivative of Eq. (43) gives ˙ ϕ ( s ) ≡ ddτ ϕ τ ( s τ ) | τ =0 = V ( s ) , (44)and the upper dot ( ˙ • ) denotes the material derivative,which is sometimes substituted by a superscript dot( • ) (cid:113) .We separate a configuration design variables intotwo sets. First , configuration changes of an embedded onfiguration Design Optimization for Maximal Fundamental Frequency 7(a) Initial parent domain (b) Modeling by FFD (c) Curve interpolation
Fig. 3: Geometric modeling of an embedded structurestructure in the initial parent domain are described bya set of design variables, and each variable has a givendesign velocity field, expressed by [3] V ( ξ ) ≡ ddτ X τ (cid:12)(cid:12) τ =0 = n (cid:88) I =1 W I ( ξ ) ˙B I = n (cid:88) I =1 W I ( ξ ) V I , (45)where V I ≡ ˙B I represents a given design velocity coef-ficient. For a given design velocity field of Eq. (45), thematerial derivative of surface parametric coordinates ofEq. (40) is obtained by [3]˙˜ ξ i = V i ( ξ ) /b i , i = 1 , , (46)where V i denotes the i -th component of V . Second ,the target parent domain has a set of design variables,and each variable has a given design velocity field, ex-pressed as ˜V ( ˜ ξ , ˜ ξ ) ≡ ddτ ˜X τ (cid:12)(cid:12)(cid:12) τ =0 = m (cid:88) J =1 ˜ W J ( ˜ ξ , ˜ ξ ) ˙˜B J = m (cid:88) J =1 ˜ W J ( ˜ ξ , ˜ ξ ) ˜V J , (47)where ˜V J denotes a given design velocity coefficient.Taking the material derivative of the surface positionof Eq. (39) yields a design velocity vector on the targetparent domain as ˙X ( ξ ) = m (cid:88) J =1 (cid:16) ˜ W J, ˜ ξ i ˙˜ ξ i ˜B J + ˜ W J ˜V J (cid:17) ≡ V ( ξ ) (48)where the repeated index i implies summation over 1and 2. ( • ) , ˜ ξ i denotes the differentiation with respect to the surface parametric coordinates ˜ ξ i ( i = 1 , N g V g = ˙X g , (49)where the same system matrix N g with that of Eq. (41)is used, and the right-hand side is constructed by an as-sembly of evaluations ˙X ( ξ cL ) of Eq. (48) at selected n discrete points ξ = ξ cL ( L = 1 , ..., n ). From the solution V g of Eq. (49), we extract a set of design velocity co-efficients { V I } for each patch in the embedded curveson the target parent domain, then the design velocityfield is finally expressed as V ( s ( ξ )) ≡ ddτ X ∗ τ | τ =0 = n (cid:88) I =1 W I ( ξ ) V I , (50)where V I denotes a design velocity coefficient at eachcontrol point.3.2 Material derivative of initial triads in embeddedstructuresFor a C -continuous spatial curve, a unit tangent vec-tor j ( ξ ) can be uniquely determined along the curve.However, the other two base vectors j ( ξ ) and j ( ξ ) ofthe orthonormal basis are not uniquely determined dueto the rotational DOF around the curve. In order to ex-plicitly parameterize local orthonormal frames, the SRmethod was employed in [16], as j k ( ξ ) = j ref k − j ref k T j ( ξ )1 + j ref1 T j ( ξ ) (cid:8) j ( ξ ) + j ref1 (cid:9) ≡ s k (cid:0) j ( ξ ) , j ref123 (cid:1) , k = 2 , , (51) M.-J. Choi et al. where j ref123 ≡ (cid:8) j ref1 , j ref2 , j ref3 (cid:9) denotes the reference basis,and s k ( • , ∗ ) denotes the SR operation which is inter-preted as rotating a given reference base vectors (*)in a way that the first reference base vector ( j ref1 ) isrotated to a given vector ( • ) with a minimal rotationangle. The SR method was also employed in [4] withinthe IGA framework, and a surface convected basis wasshown to be effectively utilized as the reference basis forcurves embedded in a smooth surface. In this section,we explain an extension of the previous work to derivematerial derivatives of initial orthonormal base vectors,considering the design dependence of target parent do-main.If a curve and its perturbed design are C -continuous,the material derivative j (cid:113) ( ξ ) can be obtained by takingthe material derivative of j ≡ ϕ ,s as [4] j (cid:113) ( ξ ) = V ,s − j ∇ s · V , (52)where ∇ s · ( • ) ≡ ( • ) ,s · j . The other two base vectors j ( ξ ) and j ( ξ ) of the orthonormal basis can be deter-mined by [4] j k ( ξ ) = s k (cid:0) j ( ξ ) , j ref123 ( ξ ) (cid:1) , k = 2 , , (53)where j ref123 ( ξ ) ≡ (cid:8) j ref1 ( ξ ) , j ref2 ( ξ ) , j ref3 ( ξ ) (cid:9) represents thereference orthonormal basis, defined by a surface con-vected basis in the following. The SR reference tangentvector, which is an exact tangent to the target parentdomain along the embedded curve, is defined as [4] j ref1 ( ξ ) ≡ ∂ ˜X /∂ξ (cid:13)(cid:13)(cid:13) ∂ ˜X /∂ξ (cid:13)(cid:13)(cid:13) , (54)where ∂∂ξ ˜X ≡ ˜X , ˜ ξ ∂ ˜ ξ ∂ξ + ˜X , ˜ ξ ∂ ˜ ξ ∂ξ = 1 b X ,ξ ˜X , ˜ ξ + 1 b X ,ξ ˜X , ˜ ξ . (55)The SR reference unit binormal and normal vectors arerespectively determined by [4] j ref3 ( ξ ) ≡ ˜X , ˜ ξ × ˜X , ˜ ξ (cid:13)(cid:13)(cid:13) ˜X , ˜ ξ × ˜X , ˜ ξ (cid:13)(cid:13)(cid:13) , (56)and j ref2 ( ξ ) ≡ j ref3 ( ξ ) × j ref1 ( ξ ) . (57)Taking the material derivative to Eq. (54) gives j ref (cid:113) ( ξ ) = 1 (cid:13)(cid:13)(cid:13) ∂ ˜X /∂ξ (cid:13)(cid:13)(cid:13) P j ref1 ∂ ˙˜X ∂ξ , (58) where P ( • ) ≡ I − ( • ) ⊗ ( • ), and I denotes the identitymatrix, and ∂ ˙˜X /∂ξ can be explicitly expressed in termsof the given design velocity fields of Eqs. (45) and (47)from Eq. (55) as ∂∂ξ ˙˜X ( ˜ ξ , ˜ ξ ) = 1 b (cid:110) V ,ξ ˜X , ˜ ξ + X ,ξ (cid:16) ˜X , ˜ ξ (cid:17) (cid:113) (cid:111) + 1 b (cid:110) V ,ξ ˜X , ˜ ξ + X ,ξ (cid:16) ˜X , ˜ ξ (cid:17) (cid:113) (cid:111) . (59)We also take the material derivative of Eq. (56) to have j ref (cid:113) = 1 (cid:13)(cid:13)(cid:13) ˜X , ˜ ξ × ˜X , ˜ ξ (cid:13)(cid:13)(cid:13) × P j ref3 (cid:110)(cid:16) ˜X , ˜ ξ (cid:17) (cid:113) × ˜X , ˜ ξ + ˜X , ˜ ξ × (cid:16) ˜X , ˜ ξ (cid:17) (cid:113) (cid:111) . (60)In Eqs. (59) and (60), (cid:16) ˜X , ˜ ξ i (cid:17) (cid:113) ( i = 1 ,
2) can be calcu-lated by (cid:16) ˜X , ˜ ξ i (cid:17) (cid:113) = ˜V , ˜ ξ i + V b ˜X , ˜ ξ i ˜ ξ + V b ˜X , ˜ ξ i ˜ ξ . (61)The material derivative of the second reference vectoris obtained from Eq. (57), as [4] j ref (cid:113) ( ξ ) ≡ j ref (cid:113) ( ξ ) × j ref1 ( ξ ) + j ref3 ( ξ ) × j ref (cid:113) ( ξ ) . (62)By taking the material derivative of the SR operationexpression of Eq. (53), the following was obtained [4]. j (cid:113) k = j ref (cid:113) k −
11 + j ref1 T j (cid:0) j ref (cid:113) k T j + j ref k T j (cid:113) (cid:1) (cid:0) j + j ref1 (cid:1) − (cid:0) j ref (cid:113) T j + j ref1 T j (cid:113) (cid:1) (cid:0) j ref k − j k (cid:1) + j ref k T j (cid:0) j (cid:113) + j ref (cid:113) (cid:1) = s (cid:113) k (cid:0) j , j ref123 ; j (cid:113) , j ref (cid:113) (cid:1) , k = 2 , . (63)where j ref (cid:113) ≡ (cid:8) j ref (cid:113) , j ref (cid:113) , j ref (cid:113) (cid:9) is obtained by Eqs. (58),(60), and (62). Then we can calculate ˙Λ = [ j (cid:113) , j (cid:113) , j (cid:113) ]using Eqs. (52) and (63). The material derivatives j ref (cid:113) I ( I = 1 , ,
3) reflect a configuration design variation oftarget parent domain by the terms regarding a givendesign velocity field of Eq. (47), in contrast to the cor-responding expressions in [4] where the designs of targetparent domains were not varied.3.3 Material derivatives of strain measuresFor future use, we recall the expression of the materialderivative of the infinitesimal arclength in [6], as( ds ) (cid:113) = ∇ s · V ds. (64)We also recall the material derivatives of the linearizedstrain measures from [3]. { Γ ( η ) } (cid:113) = Λ T ( ˙z ,s − ˙ θ × j )+ Λ T ( − z ,s ∇ s · V + j (cid:113) × θ ) + ˙Λ T Λ Γ ≡ Γ ( ˙ η ) + Γ (cid:48) V ( η ) , (65) onfiguration Design Optimization for Maximal Fundamental Frequency 9 and { Ω ( η ) } (cid:113) = Λ T ˙ θ ,s + (cid:110) Λ T ( − θ ,s ∇ s · V ) + ˙Λ T Λ Ω (cid:111) ≡ Ω ( ˙ η ) + Ω (cid:48) V ( η ) , (66)where ˙ η ≡ ( ˙z ( s ) , ˙ θ ( s )). Then, taking the material deriva-tive of Eq. (22), using Eqs. (64)-(66), yields [3][ a ( η , ¯ η )] (cid:113) = (cid:90) Ω (cid:110) Γ ( ˙¯ η ) T C F Γ ( η ) + Ω ( ˙¯ η ) T C M Ω ( η ) (cid:111) ds + (cid:90) Ω (cid:110) Γ ( ¯ η ) T C F Γ ( ˙ η ) + Ω ( ¯ η ) T C M Ω ( ˙ η ) (cid:111) ds + (cid:90) Ω Γ (cid:48) V ( ¯ η ) T C F Γ ( η )+ Ω (cid:48) V ( ¯ η ) T C M Ω ( η )+ Γ ( ¯ η ) T C F Γ (cid:48) V ( η )+ Ω ( ¯ η ) T C M Ω (cid:48) V ( η )+ (cid:26) Γ ( ¯ η ) T C F Γ ( η )+ Ω ( ¯ η ) T C M Ω ( η ) (cid:27) ∇ s · V ds ≡ a ( η , ˙¯ η ) + a ( ˙ η , ¯ η ) + a (cid:48) V ( η , ¯ η ) , (67)where ˙¯ η ≡ ( ˙¯z ( s ) , ˙¯ θ ( s )). Also, taking the material deriva-tive of Eq. (23) leads to[ d ( η , ¯ η )] (cid:113) = (cid:90) Ω (cid:18) ˙¯z T ρA z + ˙¯ θ T I ρ θ (cid:19) ds + (cid:90) Ω (cid:16) ¯z T ρA ˙z + ¯ θ T I ρ ˙ θ (cid:17) ds + (cid:90) Ω (cid:110) ¯ θ T ˙I ρ θ + (cid:16) ¯z T ρA z + ¯ θ T I ρ θ (cid:17) ∇ s · V (cid:111) ds ≡ d ( η , ˙¯ η ) + d ( ˙ η , ¯ η ) + d (cid:48) V ( η , ¯ η ) , (68)where I ρ of Eq. (16) has a configuration design depen-dence through Λ , so that ˙I ρ is obtained by ˙I ρ ≡ ˙Λ diag [ ρI p , ρI , ρI ] { e I } Λ T + Λ diag [ ρI p , ρI , ρI ] { e I } ˙Λ T . (69)It is noted that the explicit design dependence terms a (cid:48) V ( η , ¯ η ) and d (cid:48) V ( η , ¯ η ) are linear with respect to eigen-functions and the design velocity field of Eq. (50).3.4 DSA of simple eigenvaluesIn case of a simple eigenvalue, the eigenvalue and thecorresponding eigenfunction are Fr´echet differentiable[10]. Evaluating Eq. (21) at ¯ η = ˙ η since they belong tothe same function space ¯ Z , and using the symmetriesof strain and kinetic energy bilinear forms, we obtainthe following identity a ( ˙ η , η ) = ζd ( ˙ η , η ) . (70) Taking the material derivative of both sides of Eq. (21)and rearranging terms gives the following. a ( ˙ η , ¯ η ) + a (cid:48) V ( η , ¯ η ) = ζ { d ( ˙ η , ¯ η ) + d (cid:48) V ( η , ¯ η ) } + ˙ ζd ( η , ¯ η ) , ∀ ¯ η ∈ ¯ Z. (71)Evaluating Eq. (71) at ¯ η = η since they belong tothe same function space ¯ Z , and using the identity ofEq. (70) and the normalization condition of Eq. (24),we have the following configuration design sensitivityexpression for a simple eigenvalue [9].˙ ζ = a (cid:48) V ( η , η ) − ζd (cid:48) V ( η , η ) , (72)which is linear in the given design velocity field.3.5 DSA of repeated eigenvaluesConsider that the generalized eigenvalue problem of Eq.(21) yields S -fold repeated eigenvalue˜ ζ = ζ i , i = 1 , ..., S, (73)which is associated with its S linearly independent eigen-functions η i ≡ ( z i , θ i ), and any linear combinations ofthose eigenfunctions will also satisfy the Eq. (21). It wasshown that the repeated eigenvalue is not Fr´echet dif-ferentiable, but only directionally differentiable [10]. Ata perturbed design, the repeated eigenvalue may branchinto S eigenvalues as ζ iτ = ˜ ζ + τ µ i + o ( τ ) , i = 1 , ..., S, (74)where µ i denote directional derivatives of the repeatedeigenvalue, and the term o ( τ ) represents a quantitywhich approaches zero faster than τ , i.e., lim τ → o ( τ ) /τ =0 [10]. The directional derivatives µ i are derived aseigenvalues of a S × S matrix whose ( i, j ) component isexpressed by [21,10] H ij = a (cid:48) V ( η i , η j ) − ˜ ζd (cid:48) V ( η i , η j ) , i, j = 1 , ..., S. (75)It is noted that the matrix components of Eq. (75) arelinear in the design velocity field; however, the direc-tional derivatives, i.e., the eigenvalues may not be linearin the design velocity field. If the eigenvalue is simple( S = 1), the directional derivative reduces to the sensi-tivity expression of Eq. (72). D = 0 . m , and the material properties are the Young’smodulus E = 200GPa, the Poisson’s ratio ν = 0 . ρ = 7 , kg/m .Two design variables d and d are considered, as il-lustrated in Fig. 4b; the first one simultaneously changesthe horizontal and vertical space between ribs, whichmaintains the symmetrical design. The second one changesthe positions of the lowermost horizontal rib and theleftmost vertical ribs simultaneously, which asymmetri-cally varies the design. The design variables are initially d = d = 5 m . Fig. 5 shows the eigenmode shapes forthe lowest three eigenfrequencies. Hereafter, we denotethe i -th eigenfrequency by f i which is related with theangular frequency by f i = ω i / π . The first eigenfre-quency f = 16 . Hz is simple, but the second andthird eigenfrequencies are identical due to the symme-try in the structure and f = f = 30 . Hz . Resultsfrom three different DSA methods are compared; first,we calculate the directional derivatives of eigenvalues( µ i ) which are eigenvalues of the matrix of Eq. (75).Second, we calculate the sensitivities ( ˙ ζ i ) of eigenval-ues regarding them as simple using Eq. (72) in order tosee the results of this errorneous assumption. Third, inorder to check the accuracy of the aforementioned twoanalytical DSA methods, we calculate finite differencesensitivities using a forward difference method with aperturbation amount 10 − . Table 1 shows the DSA re-sults of the first three eigenvalues ζ , ζ , and ζ withrespect to two design variables d and d . As the firsteigenvalue is simple, both of the analytical DSA meth-ods give the same results which agree very well with thefinite difference one. Even though the second and thirdeigenvalues are repeated, the sensitivities ˙ ζ and ˙ ζ withrespect to the design variable d , under the erroneousassumption of simple eigenvalues, give nearly the samevalues respectively with the directional derivatives µ and µ . This is due to the symmetric design variationby changing d , which results in very weak off-diagonalterms in the matrix of Eq. (75) [21]. However, for theasymmetric design variation due to the change of d ,the assumption of simple eigenvalue gives erroneous re-sults while the directional derivatives agree very wellwith the finite difference ones. 4.2 Hexagonal lattice on torus-shaped domainWe consider an optimization problem for maximizingthe fundamental frequency under a volume constraint.The optimization problem can be stated asMaximize ψ ≡ ω j , (76)subject to g ≡ V /V − ≤ , (77)where ω j denotes the angular frequency of j -th eigen-mode, and the structural volume is calculated by V ≡ (cid:90) Ω Ads = ne (cid:88) i =1 (cid:90) Ξ i AJ c dξ. (78) A denotes the cross-sectional area. Let V denote thevolume of the original design, then an equality con-straint V = V can be implemented by employing anadditional inequality constraint g ≡ − V /V + 1 ≤ . (79)We consider a hexagonal lattice structure shown in Fig.6a, where the grey region represents an initial parentdomain, and the red-colored cell parameterized by lin-ear B-splines is repeated by a translation. 32 designvariables are selected as position changes of controlpoints in the planar domain. A toroidal surface of smallradius r = 4 m and large radius R = 10 m is constructedby the revolution of a half circle with a radius r aroundthe Z -axis, as illustrated in Fig. 6b, which is selected asa target parent domain, and consists of four quadraticNURBS patches. In this example, we assume that thetarget parent domain does not have a design depen-dence. The initial lattice embedded in the initial par-ent domain (grey-region) of Fig. 6a is mapped into eachof the four NURBS patches in the target parent do-main, for example, a green-colored region of Fig. 6b.Fig. 7 shows the constructed hexagonal lattice structurein three different views. We impose clamped bound-ary conditions at the whole bottom material points. Acircular cross-section with radius 10 cm is chosen. Thematerial properties are selected as Young’s modulus E = 210GPa, Poisson’s ratio ν = 0 .
29, mass density ρ = 7 , kg/m .We find an optimal design for maximizing the fun-damental frequency while maintaining the same volumeby solving the optimization problem of Eqs. (76), (77),and (79), where an SQP (Sequential Quadratic Pro-gramming) algorithm is utilized. The optimal design ofthe lattice structure is shown in Fig. 8. Fig. 9 shows theoptimization history where the fundamental frequencyincreases from 14.27Hz to 24.53Hz, while the volumeremains nearly the same within a specified toleranceof constraint violation. Fig. 10 compares the deformed onfiguration Design Optimization for Maximal Fundamental Frequency 11(a) Boundary condition (b) Design variables Fig. 4: Model description (a) Mode f = 16 . f = 30 . f = 30 . Fig. 5: Mode shapes (scaling factor 100 used for deformed configurations)Table 1: DSA results of the lowest three eigenvalues ( ζ i , i = 1 , , µ i (a) ˙ ζ i (b) Finite difference (c) Agreement (%)variable For a constrained structure on a planar surface, straightgeometries of structural members are often preferred fortheir manufacturability. However, even though a targetparent domain is a planar surface, the FFD processsometimes does not provide a straight curve geometry.To illustrate this issue and a solution, we consider thegeometric modeling of the curves embedded in a trape- zoidal area. An initial parent domain has a square shapewith a side length b = b = 1, and an embedded struc-ture is modeled by four single element linear B-splinepatches, as shown in Fig. 11a. A target parent domainis chosen as a trapezoidal domain shown in Fig. 11b.For each of the four linear B-spline curves in the initialparent domain, 11 equidistant points are selected, asindicated by triangular dots in Fig. 12a. Those pointsare mapped onto the target parent domain, and theyare not linearly aligned; thus, the reconstruction of eachpatch by a global curve interpolation of the 11 discretepoints does not provide a straight geometry, that is, thereconstructed curves deviate from the desired ones in-dicated by dotted lines, as illustrated in Fig. 12b. Formore detailed discussion of a nonlinear mapping be- Fig. 6: Geometric modeling of hexagonal honeycomb lattice structure on torus-shaped domain.Fig. 7: Three different views of the original design of lattice structure.tween the initial and target parent domains, see Ap-pendix A. In order to have a straight curve, it is re-quired to use a linear interpolation within each patch.For each of the four curves in the initial parent domain,we select two end points indicated by triangular dotsin Fig. 13a. Those points are mapped onto the targetparent domain, and are linearly interpolated so that astraight reconstructed geometry is obtained, as shownin Fig. 13b. These straight B-spline curves can be morerefined by a degree elevation or a knot insertion pro- cess for more accurate response analysis and DSA, andit should be noted that those refinement processes donot alter the curves either geometrically or parametri-cally.
In this example, we deal with a jacket tower struc-ture, which is modeled in a way that a lattice structureshown in Fig. 14a embedded in a rectangular domain onfiguration Design Optimization for Maximal Fundamental Frequency 13(a) Initial parent domain (b) Lattice structure
Fig. 8: Optimal design of lattice structure (a) Fundamental frequency (b) Volume
Fig. 9: Optimization historyhaving width b = 5 and height b = 10 is mappedinto each of the four surfaces comprising the prismatictarget parent domain shown in Fig. 14b, for example,the green-colored one. The target parent domain has aprismatic shape with height H = 58 . m and width B = 20 m . As illustrated in Fig. 14a, 7 configurationdesign variables are selected as control point positionsof curves embedded in the initial parent domain. Theother two design variables, indicated in Fig. 14b, changethe positions of top and bottom vertices in the pris-matic domain while maintaining the height and the flatgeometry. We seek an optimal design having straightmembers; thus, as investigated in section 4.3.1, a globalcurve interpolation process is conducted using a singlelinear B-spline element for each patch of the embed-ded curves. Fig. 14c shows the obtained jacket towerstructure, and we consider clamped boundary condi-tions at the bottom ends. The material properties areselected as Young’s modulus E = 210GPa, Poisson’sratio ν = 0 .
29, and mass density ρ = 7850 kg/m . Thecross-section has a circular shape with radius 40 cm . Us- ing the k -refinement strategy in IGA the simulationmodel is refined into 20 quartic B-spline elements ineach patch. MMFD (Modified Method of Feasible Di-rection) optimization algorithm is used. In the opti-mal design shown in Fig. 15b, it is noticeable that, thebottom dimension increases while the upper dimensiondecreases, which results in prismatic designs typicallyfound in many jacket tower structures. Fig. 16a showsthat fundamental frequency increases from 2.68Hz inthe original design to 3.88Hz in the optimal design.Meanwhile, the volume decreases, as shown in Fig. 16b.Fig. 17 compares mode shapes for the original and op-timal designs.4.4 Twisted tower modelA rectangular domain (grey region) having width b = 4and height b = 10 is considered as an initial parentdomain, where a rhombic lattice structure shown inFig. 18a is modeled by linear B-spline curves, and ithas 14 design variables d ∼ d as changes of con- f = 14 . f = 24 . Fig. 10: Comparison of the first eigenmode shape (scaling factor 200 used) (a) Initial parent domain (b) Target parent domain
Fig. 11: A modeling of curve geometry on a trapezoidal domaintrol point positions. As illustrated in Fig. 18b, a targetparent domain is a cylindrical surface which consistsof four NURBS patches, and is constructed by a sur-face of revolution of a straight line parameterized by 6control points with quadratic B-spline basis functions.The design of target parent domain is parameterizedby 6 design variables d ∼ d ; thus, total 20 con-figuration design variables are employed in this exam-ple. We map the rhombic lattice to each of the fourNURBS patches of the cylindrical surface, for example,the green-colored one. Then, we finally have the latticestructure shown in Fig. 19, where two kinds of boundary conditions are considered; first, bottom members arefixed (case E = 210GPa, Poisson’s ra-tio ν = 0 .
29, mass density ρ = 7 , kg/m . The cross-section has a circular shape with radius 40cm. We seeka structural configuration having maximal fundamentalfrequency under a volume constraint, that is, the opti-mization problem of Eqs. (76) and (77) is solved, wherea SQP algorithm is used. Figs. 20 and 21 illustrate theoptimal designs, where the changes of overall struc-tural shapes into tapered designs are noticeable. Fig. onfiguration Design Optimization for Maximal Fundamental Frequency 15(a) Initial parent domain (b) Target parent domain Fig. 12: Curved geometry of embedded object due to a nonlinear mapping between parent domains (a) Initial parent domain (b) Target parent domain
Fig. 13: A modeling of straight curve geometry on a trapezoidal domain (a) (b) (c)
Fig. 14: Geometric modeling and design variables of jacket tower structure. (a) Initial parent domain, (b) Targetparent domain, (c) Jacket tower
Fig. 15: Optimal design of jacket tower structure (a) Fundamental frequency (b) Volume
Fig. 16: Optimization history22a shows the change of fundamental frequency duringthe optimization process. At the original design, thecase
In this paper, we present a gradient-based optimal con-figuration design of three-dimensional built-up beamstructures for maximizing fundamental frequencies. Ashear-deformable beam model is employed for the anal-yses of structural vibrations within an isogeometric anal-ysis framework using the NURBS basis functions. Weparticularly extend the previous work [3] in the per-spective of considering a configuration design depen-dence of a curved surface where beam structures areembedded. This gives more design degrees-of-freedom,and provides an effective scheme for design parameter-ization and design velocity field computation in three- onfiguration Design Optimization for Maximal Fundamental Frequency 17(a) Original design ( f = 2 . f = 3 . Fig. 17: Comparison of the first mode shapes of the original and optimal designs (scaling factor 5,000 used) (a) Initial parent domain (b) Target parent domain
Fig. 18: Geometric modeling and design variables of cylindrical tower structuredimensional beam structures constrained on a curvedsurface. We also derive an analytical configuration DSAexpressions for repeated eigenvalues. The developed de-sign sensitivity analysis and optimization method isdemonstrated through various numerical examples.
Acknowledgements
M.-J. Choi and B. Koo were supported by the NationalResearch Foundation of Korea (NRF) grant funded bythe Korea government (Ministry of Science, ICT, andFuture Planning) (No. NRF-2018R1D1A1B07050370)
Conflict of interest
The authors declare that they have no conflict of inter-est.
Replication of results
All the expressions are implemented using FORTRANand the resulting data are processed by Tecplot. Theprocessed data used to support the findings of this studyare available from the corresponding author upon re-quest. If the information provided in the paper is notsufficient, interested readers are welcome to contact theauthors for further explanations.
Fig. 19: Lattice structure with two cases of boundary conditions (a) Initial parent domain (b) Twisted tower
Fig. 20: Optimal design for the boundary condition case
Appendix
In this appendix, we analytically illustrate the FFDprocess may yield a nonlinear mapping of embeddedcurves even if the target parent domain is a linear B-spline surface. We assume both of the square and trape-zoidal domains respectively for initial and target parentdomains, shown in Fig. 11, are composed of a single lin-ear B-spline surface element, and the basis functions are given by˜ W = (1 − ˜ ξ )(1 − ˜ ξ )˜ W = ˜ ξ (1 − ˜ ξ )˜ W = (1 − ˜ ξ ) ˜ ξ ˜ W = ˜ ξ ˜ ξ , (A.1)where the corresponding control point positions are ˜B (0 , ˜B (1 , ˜B (0 , ˜B (1 ,
1) for the initial parent do-main, and ˜B (0 , ˜B (2 , ˜B (0 . , ˜B (1 . , ˜B i onfiguration Design Optimization for Maximal Fundamental Frequency 19(a) Initial parent domain (b) Twisted tower Fig. 21: Optimal design for the boundary condition case (a) Fundamental frequency (b) Volume
Fig. 22: Optimization history( i = 1 , ...,
4) through Eq. (38), as ˜X ( ˜ ξ , ˜ ξ ) = (cid:20) ξ + 0 . ξ + 0 . ξ ˜ ξ ξ (cid:21) . (A.2)Using Eq. (40) and b = b = 1, the surface paramet-ric coordinates corresponding to the curve parametricposition is obtained by˜ ξ = X ˜ ξ = X (cid:27) . (A.3)Consequently, substituting Eq. (A.3) into Eq. (A.2) yields ˜X ( X , X ) = (cid:20) X + 0 . X + 0 . X X X (cid:21) , (A.4)which shows that the linear curve in the initial parentdomain is mapped into a quadratic curve in the targetparent domain. References
1. Allaire, G., Jouve, F.: A level-set method for vibrationand multiple loads structural optimization. Computermethods in applied mechanics and engineering (30-33), 3269–3290 (2005)2. Blom, A.W., Setoodeh, S., Hol, J.M., G¨urdal, Z.: Designof variable-stiffness conical shells for maximum funda-mental eigenfrequency. Computers & Structures (9),870–878 (2008)3. Choi, M.J., Cho, S.: Constrained isogeometric design op-timization of lattice structures on curved surfaces: com-putation of design velocity field. Structural and Multi-disciplinary Optimization (1), 17–34 (2018)4. Choi, M.J., Cho, S.: Isogeometric configuration de-sign sensitivity analysis of geometrically exact shear-deformable beam structures. Computer Methods in Ap-plied Mechanics and Engineering , 153–183 (2019)5. Choi, M.J., Oh, M.H., Koo, B., Cho, S.: Optimal designof lattice structures for controllable extremal band gaps.Scientific reports (1), 9976 (2019)0 M.-J. Choi et al.(a) Original design ( f = 1 . Hz ) (b) Optimal design ( f = 3 . Hz ) Fig. 23: Comparison of the first eigenmode shapes for case (a) Original design ( f = 4 . Hz ) (b) Optimal design ( f = 10 . Hz ) Fig. 24: Comparison of the first eigenmode shapes for case
6. Choi, M.J., Yoon, M., Cho, S.: Isogeometric configurationdesign sensitivity analysis of finite deformation curvedbeam structures using jaumann strain formulation. Com-puter Methods in Applied Mechanics and Engineering , 41–73 (2016)7. Du, J., Olhoff, N.: Topological design of freely vibratingcontinuum structures for maximum values of simple andmultiple eigenfrequencies and frequency gaps. Structuraland Multidisciplinary Optimization (2), 91–110 (2007)8. Goldstein, H., Poole, C., Safko, J.: Classical mechanics(2002)9. Ha, Y., Cho, S.: Design sensitivity analysis and topol-ogy optimization of eigenvalue problems for piezoelectricresonators. Smart Materials and Structures (6), 1513(2006)10. Haug, E.J., Choi, K.K., Komkov, V.: Design sensitivityanalysis of structural systems, vol. 177. Academic press (1986)11. Hu, H.T., Tsai, J.Y.: Maximization of the fundamentalfrequencies of laminated cylindrical shells with respect tofiber orientations. Journal of sound and vibration (4),723–740 (1999)12. Hughes, T.J., Cottrell, J.A., Bazilevs, Y.: Isogeometricanalysis: Cad, finite elements, nurbs, exact geometry andmesh refinement. Computer methods in applied mechan-ics and engineering (39-41), 4135–4195 (2005)13. Jihong, Z., Weihong, Z.: Maximization of structural nat-ural frequency with optimal support layout. Struc-tural and Multidisciplinary Optimization (6), 462–469(2006)14. Liu, C., Du, Z., Zhang, W., Zhu, Y., Guo, X.: Additivemanufacturing-oriented design of graded lattice struc-tures through explicit topology optimization. Journal ofApplied Mechanics (8), 081008 (2017)onfiguration Design Optimization for Maximal Fundamental Frequency 2115. Liu, H., Yang, D., Wang, X., Wang, Y., Liu, C., Wang,Z.P.: Smooth size design for the natural frequenciesof curved timoshenko beams using isogeometric analy-sis. Structural and Multidisciplinary Optimization (4),1143–1162 (2019)16. Meier, C., Popp, A., Wall, W.A.: An objective 3d largedeformation finite element formulation for geometricallyexact curved kirchhoff rods. Computer Methods in Ap-plied Mechanics and Engineering , 445–478 (2014)17. Nagy, A.P., Abdalla, M.M., G¨urdal, Z.: Isogeometricdesign of elastic arches for maximum fundamental fre-quency. Structural and Multidisciplinary Optimization (1), 135–149 (2011)18. Olhoff, N., Niu, B., Cheng, G.: Optimum design of band-gap beam structures. International Journal of Solids andStructures (22), 3158–3169 (2012)19. Picelli, R., Vicente, W., Pavanello, R., Xie, Y.: Evolu-tionary topology optimization for natural frequency max-imization problems considering acoustic–structure inter-action. Finite Elements in Analysis and Design , 56–64 (2015)20. Piegl, L., Tiller, W.: The NURBS book. Springer Science& Business Media (2012)21. Seyranian, A.P., Lund, E., Olhoff, N.: Multiple eigenval-ues in structural optimization problems. Structural opti-mization (4), 207–227 (1994)22. Shi, J.X., Nagano, T., Shimoda, M.: Fundamental fre-quency maximization of orthotropic shells using a free-form optimization method. Composite Structures ,135–145 (2017)23. da Silva, G.A.L., Nicoletti, R.: Optimization of naturalfrequencies of a slender beam shaped in a linear combina-tion of its mode shapes. Journal of Sound and Vibration , 92–107 (2017)24. Simo, J.C.: A finite strain beam formulation. the three-dimensional dynamic problem. part i. Computer methodsin applied mechanics and engineering (1), 55–70 (1985)25. Simo, J.C., Vu-Quoc, L.: A three-dimensional finite-strain rod model. part ii: Computational aspects. Com-puter methods in applied mechanics and engineering (1), 79–116 (1986)26. Taheri, A.H., Hassani, B.: Simultaneous isogeometricalshape and material design of functionally graded struc-tures for optimal eigenfrequencies. Computer Methods inApplied Mechanics and Engineering , 46–80 (2014)27. Wang, D., Friswell, M., Lei, Y.: Maximizing the naturalfrequency of a beam with an intermediate elastic sup-port. Journal of sound and vibration (3-5), 1229–1238(2006)28. Zuo, W., Xu, T., Zhang, H., Xu, T.: Fast structural opti-mization with frequency constraints by genetic algorithmusing adaptive eigenvalue reanalysis methods. Struc-tural and Multidisciplinary Optimization43