- About this Journal ·
- Abstracting and Indexing ·
- Aims and Scope ·
- Annual Issues ·
- Article Processing Charges ·
- Articles in Press ·
- Author Guidelines ·
- Bibliographic Information ·
- Citations to this Journal ·
- Contact Information ·
- Editorial Board ·
- Editorial Workflow ·
- Free eTOC Alerts ·
- Publication Ethics ·
- Reviewers Acknowledgment ·
- Submit a Manuscript ·
- Subscription Information ·
- Table of Contents

Advances in Mechanical Engineering

Volume 2013 (2013), Article ID 139498, 8 pages

http://dx.doi.org/10.1155/2013/139498

## Trivariate Isogeometric Analysis for Flexible Multibody Dynamics

Wuhan Second Ship Design & Research Institute, Wuhan 430064, China

Received 25 July 2013; Accepted 5 October 2013

Academic Editor: Xiaoting Rui

Copyright © 2013 Ting Pi. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

Isogeometric analysis (IGA) has been a fundamental step forward in the computational mechanics for the past few years, which maintains the accuracy of the description of computational domain geometry throughout the analysis process. However, the research on IGA in the area of flexible multibody dynamics is little and mainly concentrates on the univariate or bivariate NURBS geometry. This paper applies the trivariate IGA to the flexible multibody dynamics and proposes a continuum mechanics-based method to construct the system dynamic equations within the framework of IGA. A significant feature of this method is that it only employs the position coordinates of the control points as the system variables. To solve the large rotation and deformation coupled problems without introducing any rotation angles, the Green-Lagrange strain tensor is adopted. The evaluation of the elastic force and its Jacobian is easy and accurate by exploiting the appropriate mathematical transformation. In addition, the mass matrix and the generalized gravity force remain constant, and the centrifugal and Coriolis inertia forces equal zero. A numerical experiment is conducted using a thin-plate pendulum, which proves the feasibility and effectiveness of this method.

#### 1. Introduction

The kinematic description of flexible bodies that undergo large rotation and deformation is a research hotspot in the area of flexible multibody dynamics. Various methods have been proposed over the years [1], among which the floating frame of reference (FFRF) is the most widely used way. FFRF uses two sets of coordinates to describe the configuration of the deformable bodies; one set describes the position and orientation of a body-fixed coordinate system, and the other describes the deformation of the body with respect to its body-fixed coordinate system. As a consequence, the stiffness matrix used to obtain elastic forces remains the same. However, the mass matrix, centrifugal and Coriolis inertia forces, and generalized gravity forces are highly nonlinear. Moreover, the small deformation assumption limits the use of FFRF in large deformation problems.

In order to overcome the limitation of FFRF, the incremental finite-element approach introduces infinitesimal rotation angles as nodal variables. However, this approach cannot exactly describe the rigid-body motion, which is a fundamental requirement in flexible multibody dynamics.

To solve the drawback of the incremental finite-element approach, the large rotation vector formulation replaces the infinitesimal rotation angles with finite ones. However, this method has not been widely accepted due to the redundancy of representing the large rotation of the cross-section which may lead to singularity problems and unrealistic shear forces.

Absolute nodal coordinate formulation (ANCF) introduces absolute displacements and global slopes as nodal coordinates, which prevents the system equations from being highly nonlinear because the mass matrix and generalized force remain constant and centrifugal and Coriolis inertia forces equal zero. ANCF can exactly describe the rigid body modes and solve the large rotation and deformation problems, while FFRF is just effective in small deformation situations.

Different from the above methods, isogeometric analysis (IGA) [2, 3], which has been considered as a fundamental step forward in computational mechanics, aims to maintain the same exact description of the computational domain geometry throughout the analysis process. The basic idea is the use of the same NURBS-based mathematical model, which is the predominant technology used by CAD, for both modeling and FEA analysis. This makes it possible to exactly represent certain geometries that can only be approximated by polynomial functions, such as conic and circular sections. Moreover, it uses the shape functions with higher regularity than the traditional FEM. IGA exhibits many advantages and has been successfully applied in many problems, for example, contact problems [4], optimization problems [5], shell and plate problems [6–10], fluid-structure interaction problems [11], and so on. In flexible multibody dynamics, [12, 13] compare IGA and ANCF. Reference [14] analyzes the Kirchhoff-Love shells within the framework of IGA. However, most of the prior research uses the univariate or bivariate NURBS geometries to describe the deformable bodies. This paper applies the trivariate IGA to flexible multibody dynamics and proposes a continuum mechanics-based systematic method to solve the large rotation and deformation coupled problems.

The rest of the paper is organized as follows. Section 2 makes a short review of the geometric modeling and computational domain discretization techniques used in IGA. Section 3 proposes a systematic method to construct the dynamic models within the framework of IGA. Section 4 describes the computational strategies used to solve the system equations. To verify the feasibility and effect of this method, numerical tests are performed and the results are presented in Section 5. A conclusion is made at last in Section 6.

#### 2. Geometric Modeling and Space Discretization

IGA forms the computational domain in a way much the same as that of the traditional FEMs. In these methods, the position of any point in an element is evaluated by interpolation of element nodes, and the interpolation weights are determined by the element shape functions. IGA employs NURBS as the geometric description method which is based on the knot vectors, and control points. Meshing and element shape functions are formed by knot vectors and the control points play the same role as that of the element nodes.

A knot vector, , is constituted by nondecreasing real numbers, where is the order of the B-Spline and is the number of basis functions (and control points) necessary to describe it. A knot vector is said to be uniform if its knots are uniformly spaced and nonuniform otherwise. Moreover, a knot vector is said to be open if its first and last knots are repeated times. In this paper, the open knot vectors are adopted. In addition, and are set to and , respectively.

Given a knot vector, , B-Spline basis functions are defined recursively starting with (piecewise constants) as follows: For ,

If the internal knots are not repeated, B-Spline basis functions are -continuous. If a knot has multiplicity , the basis is -continuous at that knot. In particular, when a knot has multiplicity , the basis is -continuous and interpolatory at that location.

With the knot vectors and the control points, NURBS can form curves, surfaces, and solids exactly. In this paper, we concentrate on the trivariate NURBS solid, on which any point is expressed as follows: where , , and are the three parameters located in the knot vectors , , and , respectively; are the control points; are the NURBS basis functions as follows: where are weights; , , and are the B-Spline basis functions defined on the knot vectors , , and , respectively. It can be easily deduced that the number of NURBS basis functions is .

The most distinguished feature of IGA is the discretization of the computational domain. To directly perform the meshing on the same exact NURBS geometry, knot vectors , , and are generated by deleting the repeated knot values in knot vectors , , and , separately. Taking as an example, every knot span is regarded as an element in the parametric space. Then, the product of knot vectors forms a mesh of the three-dimensional parametric space. Finally, the corresponding mesh of the NURBS solids is formed by projecting these parametric elements into the physical space. Because of this “inherent” meshing scheme, IGA can not only perform FEM typical refinement strategy by means of knot insertion and degree elevation, but it can also offer an effective high-regularity refinement strategy [2].

#### 3. Dynamic Equations

The large rotation and deformation problems of the flexible bodies belong to the highly nonlinear problems. To extend IGA to these problems, a continuum mechanics-based method is proposed to describe the elastic forces.

Suppose that the position vector of an arbitrary point on a flexible body can be expressed in the global coordinate system as where is the shape function given as follows: where , an arrangement of all of the NURBS bases in (4). is the corresponding arrangement of all of the control points : where , , and are the coordinates of the control points in three directions.

Given the initial configuration and current configuration , the position vector gradients are defined as

To solve the large rotation and deformation problems within the framework of IGA, the Green-Lagrange strain tensor and the corresponding second Piola-Kirchhoff stress tensor are employed to guarantee that the work of elastic forces and strain energy are not affected when the continuum undergoes pure rigid-body motion.

The Green-Lagrange strain tensor is defined as With this strain tensor, the weak form of dynamic equilibrium equations of the continuous body in the reference configuration can be written as The first integral in (10) represents the virtual work of inertia forces, where is the density of material, is the acceleration vector, is the volume element, and the symbol means the scalar product of two vectors. The second integral is the virtual work of the internal elastic forces, where is the determinant of position vector gradients, is the Cauchy stress tensor, and the symbol : means the scalar product of two second-order tensors. The third integral is the virtual work of surface forces, where is the area element, is a unit vector directed along the outward normal to the surface, and is the mean surface traction. The fourth integral is the virtual work of the body forces, where is the vector of body forces.

##### 3.1. Inertia Forces

The virtual work of inertia forces can be written as It can be inferred from (11) that the vector of inertia forces is and the mass matrix is It can be deduced that the mass matrix is a constant symmetric matrix, which leads to zero centrifugal and Coriolis forces when the body experiences an arbitrary deformation and rotation.

##### 3.2. Elastic Forces

As shown in (10), the Green-Lagrange strain tensor is used to measure the degree of deformation in this paper. In order to deduce the expression of the elastic forces, the detailed formula of the Green-Lagrange strain tensor is needed. Suppose , , and , , are the partial derivatives of the th row of shape functions with respect to parameters , , and , respectively, and Substituting (14) into (8), the position vector gradients can be written as Substituting (15) into (9), the elements of the Green-Lagrange strain tensor can be written as

For isotropic elastic materials, the virtual work of elastic forces, , can be written as where is Lame’s constant and is the shear modulus. The vector of elastic forces is defined as the derivative of with respect to generalized coordinates : where

It can be seen from (14), (16), (18), and (19) that is the key to evaluate the elastic forces. is the collection of all of the derivatives of NURBS bases with respect to their parameters, which has been fully studied in computer graphics [15]. In the numerical integration methods, the derivatives of elastic forces with respect to generalized coordinates are necessary for the calculation of the Jacobian. Reference [16] gives an elaborately deduced form of the Jacobian:With (18) and (20), the elastic forces and the Jacobian can be effectively evaluated. In engineering problems, mechanical properties of isotropic elastic materials are usually given in terms of Young’s modulus, , and Poisson’s ratio, . The relationships between these parameters and Lame’s constants are the following:

##### 3.3. Body Forces

Taking gravity forces as an example, the virtual work of body forces can be written as It can be inferred that the vector of gravity forces is

#### 4. Computational Strategies

Multibody systems usually contain various forms of constraints. Different kinds of system equations for the constrained multibody systems have been proposed over the years [17, 18]. This paper adopts the frequently used DAE system of index-3 with holonomic constraints which is usually written as where is the system mass matrix, , , and are the position coordinates, velocities, and accelerations of control points in the global coordinate system, is the time variable, are the system constraint equations, is the Jacobian of constraints , is the vector of the Lagrange multipliers, and is the vector of generalized forces, including body forces and elastic forces.

In general, DAE-solvers can be classified into two categories [17, 18]: algebraic elimination and constraint violation stabilization. The latter is paid continuously increasing attention due to the simplicity and efficiency in implementation, among which the commonly used methods include Newmark [19], BDF [20], HHT [21], generalized- [22, 23]. The Newmark method has only order-1 convergence of errors, while the others possess order-2 convergence. In addition, the generalized- method has adjustable control of numerical dissipation; which leads to an optimal combination of high-frequency and low-frequency dissipation, that is, for a desired level of high-frequency dissipation, the low-frequency dissipation is minimized. Based on these considerations, the generalized- method is employed to solve the system equation (24). With this method, system equations of time step are discretized as where is the time step, , , , and . is a parameter controlling numerical dissipation ( for maximal dissipation). The position and velocity variables of time step are discretized as where ; is constituted by the acceleration-like auxiliary variables, which require that

Algorithm 1 sums the process of generalized- scheme to solve dynamic equilibrium at time .

The Jacobian, of system equations (25) in Algorithm 1 is

#### 5. Numerical Tests

In this section, a numerical experiment with a thin-plate pendulum is performed to verify the proposed continuum mechanics-based trivariate isogeometric analysis for the large rotation and deformation coupled problems.

As shown in Figure 1, the plate is fixed in one corner without any deformation at the initial time. Then, it falls down freely under the action of gravity. Both the length and width of the plate are set to m, and its thickness is chosen to be m. The density, Young’s modulus, and Poisson ratio of the material are assumed to be kg/cm^{3}, N/m^{2}, and , respectively.

The employed plate element adopts NURBS bases of order in the direction of both length and width. In order not to make a distraction, these NURBS bases are assumed to be -continuous at all of the internal knot values, which are uniformly distributed between and . The polynomial order in the thickness direction is chosen to be . Therefore, the number of control points in a single element is .

To better illustrate the NURBS element, Figure 2 presents the relationship between the element and its related control points when the whole plate is regarded as a single element, that is, , and . As shown in the figure, control points are marked with circles and connected with dashed line to form the control net. It should be noted that this is a special case that NURBS basis functions are interpolatory at all of the corners of the element. In most cases, control points do not need to lie on the geometry.

To see the influence of refinement on the plate’s deformation, other -elements ( = = , ,, = ) and -element ( = = , ) meshing plans are performed separately. A s simulation with s step size is performed for each of the above three meshing plans. Figure 1 shows the configurations of the moving thin plate at different moments for the -element meshing plan. Figure 3 shows the position relationship between the control points and the thin plate configuration for the single-element meshing plan at different moments. Figure 4 compares the plate deformation of single-element, -element and -element meshing plans at different moments. It is obvious that the plate configuration changes significantly when the number of elements increases from to , while it differs not that much from to . But the ratio of CPU time used in the three cases is 1 : 22 : 66. Therefore, it can be concluded that, using this method, a balance between the accuracy and efficiency can be achieved with a relatively coarse meshing.

#### 6. Conclusion

This paper applies the trivariate IGA to the flexible multibody system dynamics and proposes a continuum mechanics-based method to construct the system dynamic equations within the framework of IGA. A significant feature of this method is that it only employs the position coordinates of the control points as the system variables. To solve the large rotation and deformation coupled problems without introducing any rotation angles, the Green-Lagrange strain tensor is employed. Elastic force and its Jacobian are easily and accurately evaluated by using appropriate mathematical transformations. In addition, the mass matrix and generalized gravity force remain constant, and the centrifugal and Coriolis inertia forces equal zero. The feasibility and effectiveness of this method are verified through a thin-plate pendulum experiment. In addition, the trivariate IGA has a good symmetric form, which facilitates the computer implementation. The element types and the comparison with other methods will be studied in the future.

#### References

- A. A. Shabana, “Flexible multibody dynamics: review of past and recent developments,”
*Multibody System Dynamics*, vol. 1, no. 2, pp. 189–222, 1997. View at Scopus - T. J. R. Hughes, J. A. Cottrell, and Y. Bazilevs, “Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement,”
*Computer Methods in Applied Mechanics and Engineering*, vol. 194, no. 39–41, pp. 4135–4195, 2005. View at Publisher · View at Google Scholar · View at Scopus - J. A. Cottrell, T. J. Hughes, and Y. Bazilevs,
*Isogeometric Analysis: toward Integration of CAD and FEA*, John Wiley & Sons, New York, NY, USA, 2009. - T. Temizer, P. Wriggers, and T. J. R. Hughes, “Contact treatment in isogeometric analysis with NURBS,”
*Computer Methods in Applied Mechanics and Engineering*, vol. 200, no. 9–12, pp. 1100–1112, 2011. View at Publisher · View at Google Scholar · View at Scopus - W. A. Wall, M. A. Frenzel, and C. Cyron, “Isogeometric structural shape optimization,”
*Computer Methods in Applied Mechanics and Engineering*, vol. 197, no. 33–40, pp. 2976–2988, 2008. View at Publisher · View at Google Scholar · View at Scopus - D. J. Benson, Y. Bazilevs, M.-C. Hsu, and T. J. Hughes, “A large deformation, rotation-free, isogeometric shell,”
*Computer Methods in Applied Mechanics and Engineering*, vol. 200, no. 13–16, pp. 1367–1378, 2011. View at Publisher · View at Google Scholar · View at Scopus - D. J. Benson, Y. Bazilevs, M. C. Hsu, and T. J. R. Hughes, “Isogeometric shell analysis: the Reissner-Mindlin shell,”
*Computer Methods in Applied Mechanics and Engineering*, vol. 199, no. 5–8, pp. 276–289, 2010. View at Publisher · View at Google Scholar · View at Scopus - J. Kiendl, K.-U. Bletzinger, J. Linhard, and R. Wüchner, “Isogeometric shell analysis with Kirchhoff-Love elements,”
*Computer Methods in Applied Mechanics and Engineering*, vol. 198, no. 49–52, pp. 3902–3914, 2009. View at Publisher · View at Google Scholar · View at Scopus - L. Beirão da Veiga, A. Buffa, C. Lovadina, M. Martinelli, and G. Sangalli, “An isogeometric method for the Reissner-Mindlin plate bending problem,”
*Computer Methods in Applied Mechanics and Engineering*, vol. 209-212, pp. 45–53, 2012. View at Publisher · View at Google Scholar · View at Scopus - R. Echter, B. Oesterle, and M. Bischoff, “A hierarchic family of isogeometric shell finite elements,”
*Computer Methods in Applied Mechanics and Engineering*, vol. 254, pp. 170–180, 2013. - Y. Bazilevs, V. M. Calo, T. J. R. Hughes, and Y. Zhang, “Isogeometric fluid-structure interaction: theory, algorithms, and computations,”
*Computational Mechanics*, vol. 43, no. 1, pp. 3–37, 2008. View at Publisher · View at Google Scholar · View at Scopus - G. G. Sanborn and A. A. Shabana, “On the integration of computer aided design and analysis using the finite element absolute nodal coordinate formulation,”
*Multibody System Dynamics*, vol. 22, no. 2, pp. 181–197, 2009. View at Publisher · View at Google Scholar · View at Scopus - P. Lan and A. A. Shabana, “Integration of B-spline geometry and ANCF finite element analysis,”
*Nonlinear Dynamics*, vol. 61, no. 1-2, pp. 193–206, 2010. View at Publisher · View at Google Scholar · View at Scopus - A. Goyal, M. R. Dörfel, B. Simeon, and A. V. Vuong, “Isogeometric shell discretizations for flexible multibody dynamics,”
*Multibody System Dynamics*, vol. 30, no. 2, pp. 139–151, 2013. - L. A. Piegl and W. Tiller,
*The NURBS Book*, Springer, New York, NY, USA, 1997. - D. García-Vallejo, J. Mayo, J. L. Escalona, and J. Domínguez, “Efficient evaluation of the elastic forces and the jacobian in the absolute nodal coordinate formulation,”
*Nonlinear Dynamics*, vol. 35, no. 4, pp. 313–329, 2004. View at Publisher · View at Google Scholar · View at Scopus - A. Laulusa and O. A. Bauchau, “Review of classical approaches for constraint enforcement in multibody systems,”
*Journal of Computational and Nonlinear Dynamics*, vol. 3, no. 1, Article ID 011004, 2008. View at Publisher · View at Google Scholar · View at Scopus - O. A. Bauchau and A. Laulusa, “Review of contemporary approaches for constraint enforcement in multibody systems,”
*Journal of Computational and Nonlinear Dynamics*, vol. 3, no. 1, Article ID 011 005, 2008. - N. M. Newmark, “A method of computation for structural dynamics,”
*Journal of the Engineering Mechanics Division*, vol. 85, pp. 67–94, 1959. - C. W. Gear,
*Numerical Initial Value Problems in Ordinary Differential Equations*, Prentice Hall; PTR, New York, NY, USA, 1971. - H. M. Hilber, T. J. Hughes, and R. L. Taylor, “Improved numerical dissipation for time integration algorithms in structural dynamics,”
*Earthquake Engineering & Structural Dynamics*, vol. 5, no. 3, pp. 283–292, 1977. - J. Chung and G. M. Hulbert, “Time integration algorithm for structural dynamics with improved numerical dissipation: the generalized-
*α*method,”*Journal of Applied Mechanics, Transactions ASME*, vol. 60, no. 2, pp. 371–375, 1993. View at Scopus - O. L. Jay and D. Negrut, “A second order extension of the generalized-$\alpha $ method for constrained systems in mechanics,”
*Multibody Dynamics*, pp. 143–158, 2008.