Geometrical design modes of dynamic structures

In this work, a hybrid, that is, discrete in time and continuous in space, sensitivity analysis for dynamic structures using isogeometric analysis is presented. The main focus is placed on using a direct differentiation technique to derive sensitivity matrices for displacement, velocity and acceleration. To gain further understanding of the sensitivity information, a singular value decomposition is used to decompose these sensitivity matrices. The findings are exemplified on an academic example.

In this work, the implicit Newmark method [4] is utilized for time discretization, which approximates the actual value of displacements, velocity and acceleration, respectively.The approximations within an implicit time interval [  ,  +1 ] read with the time step length  =  +1 −   and the scaling factors ) The Newmark parameters  and  control the stability and accuracy of the integration scheme.The space discretization using IGA yields a nonlinear equation for the actual displacements  +1 , which can be solved utilizing Newton's method.In each iteration, a correction Δ = − −1    is added to the current iterate.The corresponding values of velocity and acceleration can be updated according to the Newmark update scheme, compare Equations (2)(3), which leads to the updated solution of dynamic state variables at discrete time step  +1 .The update formulae with Δ u =  4 Δ and Δ ü =  1 Δ within Newton's method for displacement, velocity and acceleration read This process is repeated until the desired accuracy is achieved.The effective stiffness matrix in each time step and solution iteration is   =  1  +  4  + , with the Rayleigh damping matrix of the form  =  1  +  2 , where the coefficients  1 and  2 correspond to mass and stiffness proportional damping.

DIRECT SENSITIVITY ANALYSIS
After determination of a dynamic solution state, DSA is performed to obtain gradient information of the current state with respect to chosen geometric design variables , see refs.[3,5].Importantly, a total design variation must not violate force equilibrium of the dynamic system as defined in the residual equation, that is, In a first step, the functional dependencies  = ( ü+1 , u+1 ,  +1 , ) are considered and the variation of the internal part yields   =  ü+1 +  u+1 +  +1 +   , where ,  and  denote the already mentioned mass, damping and stiffness matrices, and   =    is the so-called internal pseudo load matrix.The Newmark method introduces the functional dependencies for the acceleration ü+1 = ü+1 ( +1 ,   , u , ü ) and for the velocity u+1 = u+1 ( +1 ,   , u , ü ), respectively.Consequently, the total variation of acceleration and velocity can be determined as follows where  ∈ ℝ   ×  denotes the identity matrix with dimension   .Assuming the external part as independent of the dynamic state, the variation of the external part takes the form   =   +   +   +   =   , where   =

𝜕𝒔
. Rearranging Equation ( 6) leads to the definition of the total sensitivity matrix  +1 , connecting the displacements with geometric design variations in the actual time step where  =   −   has been used and  denotes the dynamic pseudo load matrix of the form where   and   are history velocity and acceleration sensitivity matrices, respectively.Within the chosen time discretization scheme, the velocity and acceleration history sensitivity matrices  +1 and  +1 , respectively, are updated at the end of each time step as follows Depending on the specific update formula of the chosen time discretization scheme, the partial derivatives of the actual values of dynamic state variables at time  +1 w.r.t.their counterparts from the prior time step   can be computed straight forward.Based on the results presented in Section 3, geometric design sensitivity information of any physical quantity of interest Φ() can be computed by directly inserting the total variations of the dynamic state in form of the sensitivity matrices  +1 ,  +1 and  +1 , which leads to the total variation expressed only in terms of design variations

SVD OF RESPONSE SENSITIVITIES
SVD is a factorization method commonly used in linear algebra.Given a real matrix  ∈ ℝ × , it can be factorized as  =    T .Here,  ∈ ℝ × and  ∈ ℝ × are matrices containing orthonormal left and right singular vectors, respectively.The matrix  ∈ ℝ × is a diagonal matrix with non-negative singular values arranged in decreasing order on the diagonal.Model reduction based on SVD of sensitivity information is applied for shell elements in ref. [7].Applying SVD on the sensitivity matrices  +1 ∈ ℝ × ,  +1 ∈ ℝ × , and  +1 ∈ ℝ × at time  +1 yields F I G U R E 1 DSE algorithm using VDSA method.DSE, Design space exploration; VDSA, variational design sensitivity analysis.
SVD provides an enhanced interpretation of these sensitivity matrices.In this context, the right singular vectors (inputs) can be understood as design modes or shape modes, while the left singular vectors correspond to the associated displacement, velocity, and acceleration responses (outputs).The algorithm that utilizes SVD of sensitivities using direct differentiation is illustrated in Figure 1.
In order to perform structural response sensitivity analysis in the current time step, the values of displacements, velocities and accelerations from the previous time step are required.They have to be initialized for the first time step.Based on these values, the nonlinear structural dynamic problem is solved utilizing Newton's method and the Newmark time integration scheme as presented in Section 2. The time step length Δ can be calculated utilizing any time adaptivity algorithm, compare for example, [8].All time dependent variables need to be updated at the end of each time step according to the chosen time integration scheme.
In a converged state, sensitivity analysis is performed to compute the gradient information of quantities of interest with respect to chosen design parameters.This process involves loading sensitivity matrices from the previous time step, while the effective stiffness matrix is already calculated from the structural analysis at the current time.After calculating the pseudo load matrix and dynamic pseudo load matrix, the sensitivity matrices at the current time step are updated.Using these updated sensitivity matrices, it is possible to calculate the gradient of quantities of interest.After determining the sensitivities, DSE of sensitivity information can be performed using SVD at time  +1 to identify the most significant design modes.

NUMERICAL EXAMPLE
A numerical example is provided, utilizing the approach presented above, to demonstrate the application of SVD for dynamic structures by using the obtained sensitivities.The IGA method is employed for structural analysis, integrating computer-aided design (CAD) and numerical analysis.Its objective is to bridge the gap between CAD geometry representation and numerical simulation by employing the same basis functions for both geometry description and structural analysis.In the case of shape sensitivity analysis using IGA, the coordinates of control points can be defined as design variables.The higher-order continuity of Non-Uniform Rational B-Splines (NURBS) basis functions offers advantages in sensitivity analysis, avoiding sharp edges.The concept of IGA is briefly introduced, followed by the presentation of the numerical example and subsequent discussion of the results.

NURBS-based isogeometric analysis
To represent a structural multi-part geometry using NURBS, it is usually divided into smaller regions called patches.Each patch can be described using a NURBS volume, which is defined by coordinates of a parameter space   ∈ .The NURBS volume in physical space with coordinates   ∈  is represented mathematically by a polynomial of predefined order  in each direction of the parameter space.Continuity conditions are defined for each spatial direction between adjacent elements of the parameter space   ∈ . Figure 2 illustrates a single patch NURBS-based isogeometric mapping.
To represent the geometry   , a set of elements with different continuities can be defined using parameter coordinates   .General rational basis functions are used to map discrete points from the parameter space   to the corresponding points in physical space.This results in the following geometry description Here,   are the coordinates of the control geometry in physical space and   ∈ ℝ   ×  contains the general rational basis functions defined as In the equation above,   ∈ ℝ 1×  contains control point weights and  ∈ ℝ 1×  represents the multivariate basis functions.In the case of three-dimensional space, the multivariate basis functions can be obtained by taking the tensor product of the univariate basis functions (  ) =  1 ⊗  2 ⊗  3 , where each univariate basis function  1 ,  2 , and  3 are B-Spline basis functions.For numerical integration using the Gaussian quadrature, an element with coordinates   ∈  is defined in parent space and integration point coordinates   ∈  0 and weights are calculated utilizing Legendre polynomials.After defining a single patch in the parameter space, knot vectors for each dimension can be calculated using the defined continuity conditions and coordinates of the parameter space.From the polynomial order and continuity conditions the required number of control points in each spatial direction can be calculated.Once the parameter space   , weights and coordinates of the control points are defined, general rational basis functions can be used to map each discrete point in the parameter space (defined as   ) to a corresponding point in physical space.This mapping creates a geometry that accurately represents the NURBS volume in physical space.Note that within the IGA framework, the values of the dynamic state variables are computed at the control point coordinates of the NURBS geometry, that is,  ∶=   , and can be projected to the actual physical geometry   within a post-processing step.

Problem description
A comprehensive 3D example, as introduced in ref. [1], is presented here.The patch of the geometry in the parametric space of the NURBS model is illustrated in Figure 3, demonstrating 12 elements with predefined continuity conditions.The geometry employs a polynomial order of two in each direction with four, three, and nine control points in the respective spatial directions, The control geometry and geometry in physical space is depicted in Figure 4.The problem incorporates the Neo-Hooke model, representing a constitutive nonlinear elasticity material.For the time integration scheme, the Newmark parameters are chosen as  = 0.25 and  = 0.5, while the time step size is set to Δ = 0.01.The problem involves the application of equal and opposite surface loads on the top surfaces.The displacement of the top surfaces is fixed in both the -direction and -direction.The applied surface load linearly increases over a duration of 0.1 s.Following this period, the surface load is fully released, and the response decay is simulated for a total duration of 3 s.In the sensitivity analysis, all control points are considered as design variables, leading to a total of  = (4 ⋅ 3 ⋅ 9) ⋅ 3 = 324 design variables to be accounted for.

Results
By utilizing the algorithm outlined in the preceding sections, as depicted in Figure 1, the computation of response sensitivities is performed at each time step.These sensitivity values undergo SVD analysis.To illustrate this, time step 39 is considered.The distribution of singular values associated with the number of principal components is shown in Figure 5.By analyzing the singular value distribution, valuable insights can be gained regarding the dominant factors driving the response and the impact of various design modes on the response fields.It becomes evident that the first few design modes have a remarkable influence on the response fields.A better understanding of the dominant structural features and their influence on the system's response can be obtained by visualizing these design modes.As an example, Figure 6 F I G U R E 7 Design modes 250 for displacement, velocity and acceleration at time step 39.
depicts the first design modes.The remaining design modes exhibit distortion in the design modes that are considered insignificant.For instance, Figure 7 shows design modes 250 as examples of such modes.On the one hand, the significant influence of the most prominent design modes on the response modes is a crucial observation in our analysis.These dominant design modes capture the primary structural characteristics.On the other hand, as we explore the higher-numbered design modes, particularly the last ones, their impact becomes increasingly negligible.

CONCLUSIONS AND OUTLOOK
In conclusion, this work has explored the design space through the utilization of SVD.The analysis of singular values has provided valuable insights into the shape modes identified in the study.These singular values have proven to be useful in quantifying the significance and importance of each mode.Consequently, this methodology can be effectively employed in future research for model reduction in shape optimization procedures.By selecting a subset of design modes, it becomes possible to represent changes in the structure's design as a linear combination of these modes.The scaling factors associated with these modes can be utilized as design variables in an adapted optimization procedure.Through this approach, the computational complexity of the optimization process can be significantly reduced.

A C K N O W L E D G M E N T S
Open access funding enabled and organized by Projekt DEAL.

F I G U R E 2
Single patch isogeometric mapping.

F I G U R E 4 F I G U R E 5 F I G U R E 6
Geometry and boundary conditions.Singular values of  +1 ,  +1 and  +1 at time step 39.Design modes 1 for displacement, velocity and acceleration at time step 39.