ISOGEOMETRIC ANALYSIS APPLIED TO NONLINEAR ELASTODYNAMICS

A numerical investigation is carried out in this work in order to evaluate the computational performance of a numerical model based on isogeometric analysis for applications on geometrically nonlinear elastodynamics. In the FEM (finite element method) practice, it is well known that the Newmark’s method is unconditionally stable for linear structural systems, but this characteristic is frequently lost when problems with nonlinear behavior are analyzed. Since numerical instability is induced owing to the lack of energy balance within every time step of the time integration process, numerical schemes with energy-conserving and controllable numerical dissipation properties are usually adopted to stabilize the time integration procedure. Therefore, the main objective of the present work is to determine if isogeometric analysis can handle dynamic problems more efficiently and more accurately than finite element models, especially for applications where the nonlinear response presents some drawbacks. In this sense, geometrically nonlinear dynamic problems are analyzed employing an isogeometric model based on NURBS (non-uniform rational B-splines), which is obtained considering the Bubnov-Galerkin weighted residual method. The kinematical description is performed using the corotational approach formulated in the context of isogeometric analysis and a transformation matrix given by the classical polar decomposition. The constitutive equation is written in terms of corotational variables according to the hypoelastic theory, where the small strain hypothesis is adopted. Large displacements and large rotations are also considered in the present scheme. Some numerical examples are analyzed to determine the behavior of the present formulation for nonlinear dynamic problems.


INTRODUCTION
Isogeometric analysis is dedicated to unifying numerical techniques adopted in geometrical design and analysis, which is accomplished by using a single parameterization framework where the same basis functions are utilized in both procedures.Geometrical de-sign and analysis have been performed independently with pre-processing programs based on computer aided design (CAD) technologies and numerical solvers based on the finite element method (FEM).However, it is frequently observed that the finite element model obtained after mesh generation does not match the geometric shape of complex models reproduced with CAD tools, since insufficient approximations may be utilized depending on the basis functions adopted in the finite element formulation.In order to eliminate this drawback, isogeometric analysis was proposed, where B-splines and NURBS, which are employed in CAD to build the geometric models, are also adopted to approximate the solution field.The basic concepts on isogeometric analysis were introduced by Hughes et al. [18] and many investigations have been performed in order to extend the applicability of the proposed formulation to different fields of computational applied mechanics (see, for instance, Bazilevs et al. [3]; Benson et al. [5]; Kiendl et al. [20]).In the field of elastodynamics, the development of numerical algorithms to simulate the dynamic response of linear and nonlinear elastic bodies is a major topic.Traditional time-stepping schemes, which present excellent stability properties in the linear range, are frequently subjected to numerical instabilities when they are applied to nonlinear problems using numerical models based on the FEM.In this sense, investigations on elastodynamics using isogeometric analysis must be performed in order to study aspects related to numerical stability in the time integration process.
Finite element analysis based on CAD formulation is referred to as isogeometric analysis, where finite element procedures are integrated with CAD numerical techniques such as B-splines and NURBS parameterizations.B-splines may be seen as a generalization of the Bézier formulation (Bézier [6], [7], [8]), with basis functions presenting minimal support with respect to degree, smoothness and domain partition.On the other hand, NURBS utilizes weighting of the control points to obtain rational geometries.Geometrical entities are then represented using rational polynomials, where conic sections such as circles and ellipses are exactly constructed taking into account projective transformations of piecewise quadratic curves.NURBS objects are also able to be refined by inserting knots into the knot vector, considering that a sequence of parametric values are specified in order to define the influence of the control points over the geometrical entity.Moreover, they possess C p-k continuity for basis functions of degree p and knot multiplicity k, variation diminishing and convex hull properties.Another important aspect related to isogeometric analysis is the isoparametric concept, since the solution space is represented with the same basis functions utilized to represent the geometry.Additional information about B-splines and NURBS parameterizations are obtained in Piegl and Tiller [26] and a comprehensive work on isogeometric analysis may be found in Cotrell et al. [11].
The dynamic equilibrium equation can be discretized in the time domain using implicit or explicit algorithms.The main advantages of explicit methods are related to the economy of computational memory and short processing times observed in each time step of the numerical analysis, which are specially noticed in nonlinear problems.However, the time steps adopted by explicit methods are strongly limited owing to stability restrictions imposed by the speed of propagation of the elastic wave in solid media, which is much higher than that observed in fluids, for example.On the other hand, implicit methods such as the Newmark method are unconditionally stable and the time steps are usually determined considering lower frequencies of the energy spectrum.
It is well known that the Newmark method is unconditionally stable only for linear systems.This unconditional stability is frequently lost when the method is applied to nonlinear problems, since numerical oscillations may arise owing to the lack of energy conservation and dissipation properties of the numerical algorithm.Stabilization schemes based on numerical dissipation have been adopted since the work presented by Newmark [24], where some dissipation was introduced to control higher vibration modes in the linear range by employing a user parameter instead of the time step.Other algorithms based on numerical dissipation were proposed, such as the the Hilber-α method (Hilber et al. [16]) and the Bossak-α method (Wood et al. [30]), where unconditional stability, second order precision and numerical dissipation of higher modes are observed.Chung and Hulbert [10] presented the generalized-α method in order to obtain a numerical model with second order accuracy where minimal dissipation is observed for lower modes and maximal dissipation is verified for higher modes when linear applications are analyzed.Unfortunately, energy conservation properties of the dynamic system are not verified.
One of the first authors to analyze numerical instabilities in time integration procedures for nonlinear dynamics is due to Hughes et al. [17], where the trapezoidal rule introduced by Newmark [24] is extended by using the Lagrange multiplier method to enforce energy conservation.Optimized parameters for the α methods were obtained by Khul and Ramm [22] to extend the formulation to nonlinear problems.The optimized parameters maintained an integration process with less numerical dissipation for lower frequencies and more dissipation on higher frequencies of the energy spectrum.In addition, the authors proposed an algorithmic combination of the internal forces evaluated at the beginning and the end of every time step of the time integration process.Khul and Crisfield [21] presented a general model where the concepts suggested by Simo and Tarnow [29] are introduced into the framework of the generalized-α method (energy conservation is verified algorithmically).Currently, a sufficient condition for numerical stability of time integration schemes is given by energy conservation or reduction during each time interval of the time integration process (see Belytschko and Schoeberle [4]; Khul and Crisfield [21]).
In the present work, a numerical model based on isogeometric analysis is developed for applications on nonlinear elastodynamics, which is obtained by extending the formulation proposed by Espath et al. [15] for applications on nonlinear elastostatics.The kinematic description of the continuum is performed using the corotational approach in the context of isogeometric analysis.A hypoelastic constitutive model is adopted utilizing corotational stress and strain tensors, where the small strain hypothesis and large displacements and rotations are considered.The numerical model is obtained by applying the Bubnov-Galerkin weighted residual method over the Cauchy's equation of motion and a Newton-Raphson scheme is adopted for linearization of the residual vector in the nonlinear range.Geometry and solution fields are approximated using NURBS basis functions according to the isoparametric concept.The generalized-α method is implemented into the isogeometric formulation in order to obtain a stable scheme for time integration.The influence of aspects related to the isogeometric discre-tization is investigated for numerical applications where numerical instabilities are expected when standard finite element models are utilized.

B-splines and NURBS
B-splines are piecewise polynomials which are built from a linear combination of basis functions with local support and controlled continuity.The coefficients of this linear combination are associated with points in space, referred to as control points, and the basis functions are associated with piecewise polynomials with the pieces being joined along the knot lines.Unlike finite elements, where each element presents its own parametric space, the parametric space for B-splines is related to the patch, which is characterized by knots spans defined according to the knot vector.
A knot vector is composed of a non-decreasing set of coordinates defined in the parametric space, which also defines the extension of the patch.The knots split the parameter space into knot spans or elements where the basis functions are smooth (C ∞ ), observing that element boundaries are represented as points, lines or surfaces in the physical space depending on the geometric topology associated with the problem analyzed.After the knot coordinates are established, the knot spans are then defined and the extent of control of the control points over the geometry is also determined.By changing the knot span lengths, more sample points can be used in regions where geometric nonlinearities are identified.If the level of refinement is insufficient, the basis may be also refined, considering that the refinement procedure maintains both the geometry and the parameterization unchanged.
Knots spans are always bounded by two consecutive knots, which constitute the basic entities for isogeometric analysis in the same manner as elements are basic entities for finite element analysis.When consecutive knots present the same value in the knot vector, knot spans of zero length are defined and knot multiplicity is obtained.Knot vectors with repeated values implicate reduction of continuity of the corresponding basis functions and knot multiplicity, which is limited by the degree of the basis function, when the basis function becomes interpolatory.The maximum level of continuity across an element boundary is determined by continuity of the basis across the corresponding knot span, since the basis functions are C p-k continuous across these knots.In addition, it is also important to notice that the support of each basis function over the knot spans is given by p+1, where p is the polynomial degree.
Knot vectors are classified as uniform if the knots are equally spaced in the parametric space.Otherwise, they are defined as non-uniform.A knot vector is also classified as open if the first and last knots appear p+1 times in the knot vector, when basis functions are interpolatory at the ends of the parametric space.At interior knots, the basis is not, in general, interpolatory, unless repeated knots are considered.The use of open vectors is very common in CAD, which also permits that patches can be assembled in isogeometric analysis in the same manner as that utilized in the FEM.A one-dimensional knot vector may be written as follows: where p is the polynomial degree and n is the number of B-spline basis functions, which is also associated with the number of control points.The Cox-de Boor recursive formulation (Cox [12]; De Boor [13]) is usually adopted to evaluate B-spline basis functions, which are obtained considering a given knot vector ξ defined over the parametric space, the number of control points and the polynomial degree.According to the Cox-de Boor formulation, the B-spline basis functions may be expressed as: where p is the polynomial degree, which is valid for p ≥ 1, and i is the knot index.The recursion is performed over the polynomial degree until p = 0 is obtained, when the following equation is utilized: Derivatives of the B-spline basis functions with respect to the parametric directions are usually required by numerical models based on isogeometric analysis, which are represented in terms of B-spline lower order bases owing to the recursive definition of the basis functions (see Espath et al. [15]).Algorithms for numerical evaluation of derivatives of B-spline basis functions may be found in Piegl and Tiller [26].
The most important properties related to B-splines may be summarized as follows: 1) In the absence of repeated knots or control points, continuous derivatives of order p-1 are maintained.
2) The basis functions are non-negative over the entire parametric domain.
3) The basis functions constitute a partition of unity for open vectors, that is: The number of continuous derivatives is decreased by k when a knot or a control point is repeated k times.
5) The affine covariance property guarantees that any transformation on a B-spline curve is obtained applying the transformation directly to the control points.

Representation of geometric entities
Given a set of n basis functions N i,p of degree p and the corresponding vector of control points P i , a B-spline curve may be obtained as: where the knot vector is specified according to Eq. ( 1).
B-spline surfaces are obtained considering a tensor product of univariate B-spline basis functions and a two-dimensional net of control points P i,j , with i ∈ [1,n] and j ∈ [1,m], that is: ) where N i,p and M j,q are B-spline basis functions of degree p and q, respectively, and the corresponding knot vectors are defined by ξ = {ξ 1 , ξ 2 ,…, ξ n+p+1 } and η = {η 1 , η 2 ,…, η m+q+1 }.B-spline solids are obtained analogously to B-spline surfaces.Given a threedimensional net of control points ) where p, q and r denote the degree of the basis functions N i,p , M j,q and L k,r , respectively, while n, m and l determine the number of control points of the corresponding basis functions.Following the terminology adopted in the finite element practice, the polynomial degrees 0, 1 and 2, for instance, refer to constant, linear and quadratic polynomial functions.Geometrical entities can be also represented using rational polynomials, where conic sections such as circles and ellipses are exactly constructed taking into account projective transformations of piecewise quadratic curves (see Roberts [28]; Riesenfeld [27]; Patterson [25]).A rational function is defined as any function that can be written as the ratio of two polynomial functions.Homogeneous coordinates may be utilized to represent rational polynomials in n-dimensional space by using polynomials in n+1-dimensional homogeneous space.When the rational concept is applied to non-uniform B-splines, non-uniform rational B-splines (NURBS) are obtained, which represent a significant improvement over standard Bsplines, since complex objects cannot be exactly represented using simple polynomials.
A NURBS curve of degree p may be defined with the following expression: 1 ( ) ( ) where the control points P i are obtained through the point mapping w → P P, that is: The rational basis functions for NURBS curves are given by: , , 1 ( ) ( ) ( ) NURBS surfaces and solids are analogously defined considering that: ( , , ) ( , , ) with the rational basis functions given by: , , , where subscripts have the same meaning as those previously defined.If the weights are specified with unit values, then p i R = N i,p .Derivatives of NURBS basis functions are evaluated here considering the formulation presented by Espath et al. [15], which is based on the numerical algorithms proposed by Piegl and Tiller [26].

THE COROTATIONAL APPROACH FOR NONLINEAR ISOGEOMETRIC ANALYSIS
Corotational formulations adopted in finite element models are usually defined considering local coordinate systems positioned at the quadrature points of the elements constituting the finite element mesh.The same procedure is employed in isogeometric analysis, taking into account that the motion at element level is now described in terms of displacements defined at the control points.Assuming that all kinematical variables at the previous configuration t n of the body are known, the displacement field at the end of the current load step can be obtained from integration of the strain rate tensor over the time interval defining the present load increment [t n , t n+1 ].In addition, this integration to obtain the strain increment must be performed in the corotational coordinate system, where only the deformational part of the incremental displacement field is considered.The strain rate tensor in the corotational system is defined as: where def v represents the velocity field associated with the deformation part of the motion in the corotational system x .In order to obtain strain increments, some methodology must be adopted to integrate the strain tensor over the time interval [t n , t n+1 ].According to the midpoint integration (see Hughes and Winget [19]), the strain increment may be obtained from: where def ∆u is the deformation part of the displacement increment in the corotational system and is the intermediate configuration of the body defined in the corotational system, which can be determined according to the following expression: ( ) x +x (16) where R n+1/2 is the orthogonal transformation matrix performing rotation from the global system to the corotational system defined locally at the intermediate configuration t n+1/2 .The displacement increment referring to the present time interval [t n , t n+1 ] can be decomposed as follows: where ∆u def and ∆u rot are, respectively, the deformation and rotation parts of the displacement increment defined in the global coordinate system.It is important to notice that the decomposition described in Eq. ( 17) is locally performed at element level.The deformation displacement increment in the corotational system can be obtained from the following expression: where the transformation matrix is evaluated at the intermediate configuration t n+1/2 of the current time interval [t n , t n+1 ], since the strain rate tensor must be referred to the body configuration at t n+1/2 .Coordinates corresponding to the previous and current configurations of the body in the corotational system are obtained with following transformations: where R n and R n+1 are orthogonal transformation matrices performing rotations from the global system to the corotational system defined locally at t n and t n+1 , respectively.A hypoelastic constitutive formulation is very effective for corotational descriptions, since the nonlinear problem can be posed in rate form by considering the small strain hypothesis and an objective rate of the Cauchy stress tensor.Consequently, after determining the strain increment in the corotational system, strain and stress updates can be performed with the following equations: where n and n+1 denote the previous and current configurations in the corotational system, respectively.In order to obtain an incrementally objective constitutive formulation, the Truesdell rate tensor is adopted in this work, which may be described as follows: where ˆˆ= + L ε ω is the spatial velocity gradient tensor defined in the corotational system.The corotational spin tensor ˆ ω must be also integrated over the time interval [t n , t n+1 ] considering the same mid-point rule adopted in Eq. ( 15).In the present work, a classical polar decomposition is utilized to obtain the orthogonal transformation matrix R, where spectral decomposition of the right Cauchy-Green deformation tensor C is adopted to obtain the right stretch tensor U.The right Cauchy-Green deformation tensor is defined as function of the deformation gradient tensor F as follows: T = C F .F (22) where the polar decomposition theorem F = Q.U is invoked in order to obtain the relation between C and U, that is: By using spectral decomposition of C, the following expression is obtained: where λ i and N i are, respectively, the eigenvalues and the eigenvectors of C. Consequently, the rotation tensor Q can be evaluated from: ( ) The transformation matrix utilized in the corotational formulation is obtained considering that R = Q T .

THE NUMERICAL MODEL
Problems on elastodynamics may be formulated considering the Cauchy's equation of motion, where mass and energy conservation must be also enforced over the volume enclosing the body (see Malvern [23]).Considering a classical Lagrangian kinematical description in the Cartesian coordinate system and in the absence of temperature changes, the conservation equations are reduced to the following expressions: where X and x are vectors containing components of the material (X i ) and spatial (x i ) coordinates in the Cartesian coordinate system, respectively, t represents time, ρ is the specific mass, ∇ is the differential operator, b is the body force vector, σ contains components of the Cauchy stress tensor and v is the velocity vector.The Cauchy's equation of motion is derived considering the current deformed configuration of the body (Ω).
In the present model, geometrically nonlinear problems are analyzed taking into account the corotational formulation presented in the previous section, where stress and strain are described according to a coordinate system locally attached to every element of the physi-cal mesh.A linear constitutive model restricted to small strains is adopted in order to relate strain and stress measures, which may be written as: where σ and ε are the Cauchy stress tensor and the small strain tensor, both defined in the corotational system, C is the fourth-order elastic tensor, which may be described in terms of the Lamé constants λ and µ.When infinitesimal displacements and rotations are observed, the geometrical linear approach can be utilized, where the undeformed configuration of the body (Ω 0 ) is taken as reference throughout the analysis.
Applying the Bubnov-Galerkin weighted residual method in conjunction with the Green-Gauss theorem over the equation of motion given in Eq. ( 26), the following expression is obtained: where t is the traction vector and Ω and Γ are, respectively, volume and boundary surface referring to the physical space where the problem takes place.In order to define the element concept in the context of isogeometric analysis, geometry, velocities, displacements and displacement variations are represented with the following expressions: where R a is the NURBS basis function related to control point a, which is defined as function of the parametric coordinates (ξ,η,ζ), and n cp is the number of global control points.Knot vectors corresponding to the different directions in the parametric space must be specified defining the non-zero knot spans where elements are then identified.Consequently, the equation of motion given by Eq. ( 28) can be rewritten as: where Ω e and Γ e are, respectively, volume and boundary surface corresponding to element e in the physical mesh.Considering n, m and l as the number of basis functions related to the parametric directions ξ, η and ζ, respectively, and their respective polynomial degrees denoted by p, q and r, element e is defined by determining the indices at which the corresponding non-zero knot span begins in the index space, that is: where p+1 ≤ i ≤ n, q+1 ≤ j ≤ m and r+1 ≤ k ≤ l.The total number of elements in which the spatial field is discretized in the parametric domain is defined as: By substituting the NURBS approximation related to the displacement field (see Eq. 29) into the constitutive equation (see Eq. 27), an element level approximation of the stressstrain relation is obtained, where the strain components in the corotational system are given by: where B and û are the gradient matrix and the displacements field, which are evaluated re- ferring to the current configuration of the body in the corotational coordinate system.When infinitesimal displacements and rotations are observed, Eq. ( 33) is described in terms of the undeformed configuration of the body (Ω 0 ).Introducing the expansions of velocity and displacements components and their corresponding variations (see Eq. 29) and the relationship given by Eq. (33) into Eq.( 30), a matrix equation representing a system of algebraic equations is obtained for the equation of motion, which may be expressed as: where M e and K e are the element mass and element stiffness matrices, respectively, and f e is the force vector at element level.The matrix and vector dimensions associated to M e and K e , and f e , are specified as (n eq .neq .neq ) and (n eq ), respectively, where n eq = n en .ndf , with n df denoting the number of degrees of freedom at the control point level.The summation symbol indicates the assembling procedure to evaluate the global system of equations, considering the element contributions given according to connectivity relations established among the control points.The global stiffness matrix is always sparse because the support of each basis function is highly localized.
In the geometrically nonlinear regime, the system of equations represented by Eq. 34 must be iteratively satisfied using the incremental approach (see Bathe [1]), since internal forces are given now as functions of the current configuration of the body.The nonlinear equation of motion is obtained employing a linearization procedure given by the Newton-Raphson method, where the residual vector is submitted to a Taylor series expansion within the increment interval [t n , t n+1 ].Consequently, Eq. (34) must be rewritten as follows: where K tan is the tangent stiffness matrix.At each iterative step, the tangent stiffness matrix and the internal force vector are initially evaluated in the corotational coordinate system with the following expressions: where ˆe Ω is referred to the current configuration of element e in the corotational coordinate system, D and σ are stress tensors related to the Truesdell rate tensor and the corotational Cauchy stress tensor, respectively, with both evaluated in the corotational coordinate system.See Duarte Filho and Awruch [14] and Braun and Awruch [9] for further details.
In order to solve the system of nonlinear equations, the tangent stiffness matrix and the internal force vector must be obtained in the global coordinate system through an objective transformation from the corotational system, that is: where R is the transformation matrix defined in the previous section.When a numerical model based on isogeometric analysis is formulated with the Galerkin method and NURBS basis functions, homogeneous boundary conditions are exactly enforced by setting the corresponding control variables as zero.A trivial procedure for imposition of essential boundary conditions is then obtained, which is similar to that utilized by finite element models.In the present model, the Kroenecker delta property of the NURBS basis functions can be applied on the displacement field as follows: where vector x B specifies Cartesian coordinates of control points with parametric coordinates defined by ξ B , which are located at boundary knots with essential boundary conditions.
In the field of elastodynamics, the Newmark method is usually adopted for the solution of the equation of motion in the time domain.The Newmark scheme is an extension of the linear acceleration method in which a linear variation of acceleration is assumed in the time interval [t n , t n+1 ].Although this scheme is unconditionally stable for linear problems, it may be unstable for nonlinear problems.In the generalized-α method, the equilibrium of the equation of motion is verified at a general mid-point instead of the end-point used by the classical Newmark scheme.The modified equation of motion, including damping effects, is presented in the following form (see, for instance, Khul and Ramm [22]): ) with subscripts denoting time positions within the time interval [t n , t n+1 ] as functions of the time integration parameters α m and α f , where n and n+1 correspond to initial and end points of the time interval ∆t, respectively.Introducing Newmark approximations for n+1 u and n+1 u (see, for instance, Bathe [1]) into the expressions given by Eq. ( 40), the number of unknowns of the modified equation of motion is reduced to the displacements vector n+1 u , as it is demonstrated below: Substituting Eqs. ( 41) and (42) into the modified equation of motion (Eq.39), the effective dynamic equation can be obtained in the following form: where: The time integration parameters α, δ, α m and α f are defined as functions of the spectral radius r ∝ ( 0 r 1 ∞ ≤ ≤ ) according to: ( ) An algorithm referring to the numerical model proposed in this work is presented in Table 1, where TOL is identified as a tolerance criterion and a 0 , a 1 , a 2 , a 3 , a 4 and a 5 are constants of the generalized-α method, which are obtained from the time integration parameters α, δ, α m and α f .Table 1.Numerical algorithm for the model proposed in this work.
2. Compute the time integration parameters α, δ, α m and α f 3. Compute the time integration constants a 0 , a 1 , a 2 , a 3 , a 4 and a 5 :

Cantilever beam
A two-dimensional cantilever beam subject to pressure loading and undergoing large displacements is analyzed in this example, where plane strain state is also considered.Geometrical and load description for the present simulation are shown in Fig. 1 and material properties of the structure as well as the time step adopted in the time integration procedure are found in Table 2. Information on computational parameters utilized in the numerical analyses carried out here, which are referring to the isogeometric formulation and the generalizedα method, are summarized in Table 3.  Predictions obtained with the numerical model proposed in this work are shown in Fig. 2, where the structural response of the beam is presented taking into account displacements evaluated at its tip.The response presented here is related to a specific simulation performed considering basis functions with C 4 continuity and a time integration process characterized by a spectral radius r ∞ = 0.95.One can verify that the present results reproduced important aspects of the motion related to the beam investigated in this example.The displacement amplitude is not decreasing with time, which is frequently observed when excessive numerical dissipation is employed, and the period of vibration is maintained.It is important to notice that a stable solution with the same characteristics was also observed when a spectral radius r ∞ = 0.90 was adopted.However, numerical instabilities were identified when the spectral radius was greater than r ∞ = 0.95, which is probably associated with the lack of sufficient dissipation to maintain numerical stability in this case.Studies performed with respect to the continuity class of the basis functions demonstrated that this aspect is not significant for the dynamic response of the cantilever beam analyzed here.When compared to numerical predictions presented by Braun and Awruch [9], where a finite element formulation with eight-node hexahedral elements and one-point quadrature was utilized, the results shown here were obtained with less amount of numerical dissipation and a smaller time step.In addition, the range of stable spectral radii obtained with the isogeometric model was significantly wider than that presented by the finite element model.The structural response presented in Fig. 2 is very similar to that presented by Bathe and Baig [2].Time histories of the energy variables utilized in this work are plotted in Fig. 3, which also correspond to a specific simulation performed considering a spectral radius r ∞ = 0.95 and C 4 continuity.A similar response is obtained for r ∞ = 0.90 and the different continuity classes investigated.Therefore, stable solutions can be obtained with the generalized-α method when adequate numerical damping is employed.By comparing the present model with other dissipation-based algorithms (see, for instance, Braun and Awruch [9]), the scheme presented in this work is not disturbed by excessive numerical damping, which leads to excessive decrease of the energy and motion variables (total energy, displacement, velocity and acceleration).Other important aspect related to isogeometric analysis is the mathematical characteristic of the NURBS basis functions, which lead to much better approximations for the modes of vibration when compared with finite element formulations with linear interpolation functions.The motion of the cantilever beam during the period of oscillation [0.07s, 0.15s] can be visualized in Fig. 4, where the high nonlinear behavior is evidenced.
where E pot , E kin and E tot are the strain, kinetic and total energy, respectively, W ext is the work done by external forces, J x is angular momentum around x axis, J y is angular momentum around y axis and J z is angular momentum around z axis.The symbols (:), (.) and ( ⊗ ) denote tensor, scalar and vector products, respectively.The vectors x and x are the position vector in global coordinates and its respective time derivative, nel is the total number of elements, Ω e is the element volume and σ and ε are stress and strain energetically conjugated tensors.

Toss rule
A numerical investigation of the plane movement of a toss rule is performed in this example, where a geometrically nonlinear dynamic analysis is carried out.Geometry and load information for the present application are described in Fig. 5 and material properties of the structure as well as the time step adopted in the time integration procedure are presented in Table 4.It is important to notice that distributed loads are applied to the structure to produce the plane motion of the rule, which is free to fly in the absence of displacement restrictions and gravity action.Computational parameters regarding the numerical analyses performed here are presented in Table 5.In Fig. 6 the evolution of the different energy variables associated with the vibration time history of the toss rule is presented taking into account the time interval covered by the numerical analysis.The results shown in Fig. 6 refer to only one of many cases investigated here, which was carried out considering basis functions with C 4 continuity and spectral radius r ∝ = 0.95.It was observed that the generalized-α method with r ∝ = 1.0 led to numerical failure in the time integration procedure, since that value corresponds to the case of no algorithmic damping.On the other hand, a stable solution was obtained even with small numerical dissipation, i.e. r ∝ = 0.99, where the total energy is perfectly maintained during the time interval of the numerical analysis.The same behavior was observed for the components of the angular and linear momenta, which demonstrates that the formulation proposed in this work presents excellent performance with respect to numerical stability associated with the time integration process.Studies were also performed considering the continuity class of the basis functions, where an improvement can be identified when C 4 continuity is adopted instead of lower values.Results with C 1 , C 2 and C 3 continuity show no difference.Unlike the numerical predictions obtained by Braun and Awruch [9], where a finite element formulation with linear elements was utilized, the present results show no fluctuations in the time history of the total energy.In addition, the range of stable spectral radii obtained with the isogeometric model is slightly wider that that presented by the finite element model.The range of energy values obtained with the present algorithm is in agreement with those obtained by Kuhl and Ramm [21].The motion referring to the toss rule can be visualized in Fig. 7, where a sequence of deformed configurations obtained with the algorithm proposed in this paper is shown.One can observe that the inertial motion is developed after the initial load is removed and structural displacements take place on the plane x-z in accordance with the initial load configuration.
Figure 7. Deformed configurations of the rule during numerical analysis.

CONCLUSIONS
In the present work, a numerical model based on isogeometric analysis was proposed for investigations on geometrically nonlinear elastodynamics.In order to stabilize the time integration procedure in the nonlinear range, the generalized-α method was adopted in conjunction with an algorithmic combination of the internal forces at the beginning and the end of each time step.The simulations performed with the generalized-α method obtained stable solutions when appropriate numerical dissipation was employed.Conservation of energy and angular momentum were shown through the examples.The algorithm keeps also the great advantage of the dissipation-based models, which refers to the low computational costs de-manded by one-step schemes (similar to the classical Newmark's scheme).Results demonstrated that isogeometric analysis can provide numerical predictions accurately with much less numerical dissipation than that required by a finite element model previously developed by the present authors.Moreover, the range of stable spectral radii is wider when isogeometric analysis is carried out.Other important aspect related to isogeometric analysis is the mathematical characteristic of the NURBS basis functions, which leads to better approximations for the modes of vibration observed in dynamic analyses.A significant improvement for the numerical model proposed here, which refers to a more extended applicability of the present formulation, refers to the treatment of general material models.Therefore, the present scheme will be reformulated in future works in order to take into account hyperelastic and elastoplastic constitutive models.

Figure 1 .
Figure 1.Geometrical and load characteristics for the cantilever beam analysis.

Figure 2 .
Figure 2. Displacement response for the cantilever beam analysis using r ∞ = 0.95 and C 4 .

Figure 3 .
Figure 3. Energy response for the cantilever beam analysis using r ∞ = 0.95 and C 4 .

Figure 4 .
Figure 4. Deformed configurations of the beam during the time interval [0.07s, 0.15s].The energy variables are evaluated here utilizing the following expressions: { }

Figure 6 .
Figure 6.Energy response for the motion analysis of the toss rule using r ∞ = 0.95 and C 4 .

Table 2 .
Material properties and time step utilized in the cantilever beam analysis.

Table 3 .
Computational parameters employed in the cantilever beam analysis.

Table 5 .
Computational parameters employed in the toss rule analysis.