Trivariate Isogeometric Analysis for Flexible Multibody Dynamics

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, Ω = {ξ₁, ξ₂, …, ξ_{n + p + 1}}, is constituted by nondecreasing real numbers, where p is the order of the B-Spline and n 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 p + 1 times. In this paper, the open knot vectors are adopted. In addition, ξ₁ and ξ_{n + p + 1} are set to 0 and 1, respectively.

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

$N_{i,\mathrm{0}}\left(\xi \right)=\{\begin{array}{cc}\mathrm{1},\hfill & \xi ∊\left[{\xi }_i,{\xi }_{i+\mathrm{1}}\right),\hfill \\ \mathrm{0},\hfill & \xi \notin \left[{\xi }_i,{\xi }_{i+\mathrm{1}}\right).\hfill \end{array}$ (1)

For p > 1,

$N_{i,p}\left(\xi \right)=\frac{\xi -{\xi }_i}{{\xi }_{i+p}-{\xi }_i}{N}_{i,p-\mathrm{1}}\left(\xi \right)+\frac{{\xi }_{i+p+\mathrm{1}}-\xi }{{\xi }_{i+p+\mathrm{1}}-{\xi }_{i+\mathrm{1}}}{N}_{i+\mathrm{1},p-\mathrm{1}}\left(\xi \right).$ (2)

If the internal knots are not repeated, B-Spline basis functions are C^p−1-continuous. If a knot has multiplicity k, the basis is C^p–k-continuous at that knot. In particular, when a knot has multiplicity p, the basis is C⁰-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 C is expressed as follows:

$\mathbf{C}\left({\xi }_X,{\xi }_Y,{\xi }_Z\right)=\sum _{i=\mathrm{1}}^n\hspace{0.17em}\hspace{0.17em}\sum _{j=\mathrm{1}}^m\hspace{0.17em}\hspace{0.17em}\sum _{k=\mathrm{1}}^l\hspace{0.17em}{S}_{i,j,k}\left({\xi }_X,{\xi }_Y,{\xi }_Z\right){\mathbf{e}}_{i,j,k},$ (3)

where ξ _X , ξ _Y , and ξ _Z are the three parameters located in the knot vectors Ω _X , Ω _Y , and Ω _Z , respectively; e_{i, j, k} are the control points; S_{i, j, k} (ξ _X , ξ _Y , ξ _Z ) are the NURBS basis functions as follows:

$S_{i,j,k}\left({\xi }_X,{\xi }_Y,{\xi }_Z\right)=\frac{N_{i,u}\left({\xi }_X\right)N_{j,v}\left({\xi }_Y\right)N_{k,w}\left({\xi }_Z\right){\omega }_{i,j,k}}{{\sum }_{\widehat{i}=\mathrm{1}}^n\hspace{0.17em}{\sum }_{\widehat{j}=\mathrm{1}}^m\hspace{0.17em}{\sum }_{\widehat{k}=\mathrm{1}}^l{N}_{\widehat{i},u}\left({\xi }_X\right)N_{\widehat{j},v}\left({\xi }_Y\right)N_{\widehat{k},w}\left({\xi }_Z\right){\omega }_{\widehat{i},\widehat{j},\widehat{k}}},$ (4)

where ω_{i, j, k} are weights; {N_{i, u} (ξ _X )}, N_{j, v} (ξ _Y ), and N_{k, w} (ξ _Z ) are the B-Spline basis functions defined on the knot vectors Ω _X , Ω _Y , and Ω _Z , respectively. It can be easily deduced that the number of NURBS basis functions is N _e = n × m × l.

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 ${\widehat{\Omega }}_X$ , ${\widehat{\Omega }}_Y$ , and ${\widehat{\Omega }}_Z$ are generated by deleting the repeated knot values in knot vectors Ω _X , Ω _Y , and Ω _Z , separately. Taking ${\widehat{\Omega }}_X$ as an example, every knot span $[{\widehat{\xi }}_i,{\widehat{\xi }}_{i+\mathrm{1}})$ is regarded as an element in the parametric space. Then, the product of knot vectors ${\widehat{\Omega }}_X\times {\widehat{\Omega }}_Y\times {\widehat{\Omega }}_Z$ 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 r of an arbitrary point on a flexible body can be expressed in the global coordinate system as

$\mathbf{r}=\mathbf{S}\mathbf{e},$ (5)

where S is the shape function given as follows:

$\mathbf{S}={\left[\begin{array}{ccc}{\mathbf{S}}_{\mathrm{3}D}& \mathit{0}& \mathit{0}\\ \mathit{0}& {\mathbf{S}}_{\mathrm{3}D}& \mathit{0}\\ \mathit{0}& \mathit{0}& {\mathbf{S}}_{\mathrm{3}D}\end{array}\right]}^T,$ (6)

where ${\mathbf{S}}_{\mathrm{3}D}={[S_{\mathrm{1}},S_{\mathrm{2}},\dots ,S_{N_e}]}^T$ , an arrangement of all of the NURBS bases in (4). e is the corresponding arrangement of all of the control points e_{i, j, k}:

$\mathbf{e}={\left[\begin{array}{ccc}{\mathbf{x}}_{\mathbf{e}}^T& {\mathbf{y}}_{\mathbf{e}}^T& {\mathbf{z}}_{\mathbf{e}}^T\end{array}\right]}^T,$ (7)

where x _e , y _e , and z _e are the coordinates of the control points in three directions.

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

$\mathbf{J}=\frac{\partial \mathbf{r}}{\partial {\mathbf{r}}_o}=\frac{\partial \mathbf{r}}{\partial \mathbf{\xi }}\frac{\partial \mathbf{\xi }}{\partial {\mathbf{r}}_o}=\frac{\partial \mathbf{r}}{\partial \mathbf{\xi }}{\left(\frac{\partial {\mathbf{r}}_o}{\partial \mathbf{\xi }}\right)}^{-\mathrm{1}}.$ (8)

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

$\mathbf{∊}=\frac{\mathrm{1}}{\mathrm{2}}\left({\mathbf{J}}^T\mathbf{J}-{\mathbf{I}}_{\mathrm{3}\times \mathrm{3}}\right).$ (9)

With this strain tensor, the weak form of dynamic equilibrium equations of the continuous body in the reference configuration can be written as

${\int }_V\hspace{0.17em}\rho \mathbf{a}·\delta \mathbf{r}dV+{\int }_V\hspace{0.17em}\left|\mathbf{J}\right|{\mathbf{J}}^{-\mathrm{1}}\mathit{\sigma }{\mathbf{J}}^{-T}:\delta \mathbf{∊}dV-{\int }_A\hspace{0.17em}\frac{\left|\mathbf{J}\right|}{\sqrt{{\mathbf{n}}^T\mathbf{J}{\mathbf{J}}^T\mathbf{n}}}{\mathit{\sigma }}_n·\delta \mathbf{r}dA-{\int }_V\hspace{0.17em}{\mathbf{f}}_b·\delta \mathbf{r}dV=\mathrm{0}.$ (10)

The first integral in (10) represents the virtual work of inertia forces, where ρ is the density of material, a is the acceleration vector, dV 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 |J| 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 dA is the area element, n is a unit vector directed along the outward normal to the surface, and σ _n is the mean surface traction. The fourth integral is the virtual work of the body forces, where f _b is the vector of body forces.

3.1. Inertia Forces

The virtual work of inertia forces can be written as

${\int }_V\hspace{0.17em}\rho \mathbf{a}·\delta \mathbf{r}dV={\int }_V\hspace{0.17em}\rho {\mathbf{a}}^T\delta \mathbf{r}dV={\stackrel{‥}{\mathbf{e}}}^T{\int }_V\hspace{0.17em}\rho {\mathbf{S}}^T\mathbf{S}dV\delta \mathbf{e}={\left({\int }_V\hspace{0.17em}\rho {\mathbf{S}}^T\mathbf{S}dV\stackrel{‥}{\mathbf{e}}\right)}^T\delta \mathbf{e}.$ (11)

It can be inferred from (11) that the vector of inertia forces is

${\mathbf{Q}}_i={\int }_V\hspace{0.17em}\rho {\mathbf{S}}^T\mathbf{S}dV\stackrel{‥}{\mathbf{e}},$ (12)

and the mass matrix is

$\mathbf{M}={\int }_V\hspace{0.17em}\rho {\mathbf{S}}^T\mathbf{S}dV.$ (13)

It can be deduced that the mass matrix M 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 S_{i, ξ1}, S_{i, ξ2}, and S_{i, ξ3}, i = 1, 2, 3, are the partial derivatives of the ith row of shape functions S with respect to parameters ξ₁, ξ₂, and ξ₃, respectively, and

${\left({\mathbf{S}}_{,j}\right)}_i\triangleq {\left(\frac{\partial {\mathbf{r}}_o}{\partial \mathbf{\xi }}\right)}_j^{-T}{\left[\begin{array}{ccc}{\mathbf{S}}_{i,{\xi }_{\mathrm{1}}}^T& {\mathbf{S}}_{i,{\xi }_{\mathrm{2}}}^T& {\mathbf{S}}_{i,{\xi }_{\mathrm{3}}}^T\end{array}\right]}^T,\hspace{1em}i,j=\mathrm{1,2},\mathrm{3}.$ (14)

Substituting (14) into (8), the position vector gradients can be written as

$\mathbf{J}=\left[\begin{array}{ccc}{\mathbf{S}}_{,\mathrm{1}}\mathbf{e}& {\mathbf{S}}_{,\mathrm{2}}\mathbf{e}& {\mathbf{S}}_{,\mathrm{3}}\mathbf{e}\end{array}\right].$ (15)

Substituting (15) into (9), the elements of the Green-Lagrange strain tensor can be written as

$\begin{array}{c}{\mathbf{∊}}_{\mathrm{11}}=\frac{\mathrm{1}}{\mathrm{2}}\left({\mathbf{e}}^T{\mathbf{S}}_{,\mathrm{1}}^T{\mathbf{S}}_{,\mathrm{1}}\mathbf{e}-\mathrm{1}\right),\hspace{1em}\hspace{1em}{\mathbf{∊}}_{\mathrm{22}}=\frac{\mathrm{1}}{\mathrm{2}}\left({\mathbf{e}}^T{\mathbf{S}}_{,\mathrm{2}}^T{\mathbf{S}}_{,\mathrm{2}}\mathbf{e}-\mathrm{1}\right),\\ {\mathbf{∊}}_{\mathrm{33}}=\frac{\mathrm{1}}{\mathrm{2}}\left({\mathbf{e}}^T{\mathbf{S}}_{,\mathrm{3}}^T{\mathbf{S}}_{,\mathrm{3}}\mathbf{e}-\mathrm{1}\right),\hspace{1em}\hspace{1em}{\mathbf{∊}}_{\mathrm{12}}={\mathbf{∊}}_{\mathrm{21}}=\frac{\mathrm{1}}{\mathrm{2}}\left({\mathbf{e}}^T{\mathbf{S}}_{,\mathrm{1}}^T{\mathbf{S}}_{,\mathrm{2}}\mathbf{e}\right),\\ {\mathbf{∊}}_{\mathrm{13}}={\mathbf{∊}}_{\mathrm{31}}=\frac{\mathrm{1}}{\mathrm{2}}\left({\mathbf{e}}^T{\mathbf{S}}_{,\mathrm{1}}^T{\mathbf{S}}_{,\mathrm{3}}\mathbf{e}\right),\hspace{1em}\hspace{1em}{\mathbf{∊}}_{\mathrm{23}}={\mathbf{∊}}_{\mathrm{32}}=\frac{\mathrm{1}}{\mathrm{2}}\left({\mathbf{e}}^T{\mathbf{S}}_{,\mathrm{2}}^T{\mathbf{S}}_{,\mathrm{3}}\mathbf{e}\right).\end{array}$ (16)

For isotropic elastic materials, the virtual work of elastic forces, U _s , can be written as

$\begin{array}{c}{U}_s={\int }_V\hspace{0.17em}\frac{\lambda +\mathrm{2}G}{\mathrm{2}}\left({\mathbf{∊}}_{\mathrm{11}}^{\mathrm{2}}+{\mathbf{∊}}_{\mathrm{22}}^{\mathrm{2}}+{\mathbf{∊}}_{\mathrm{33}}^{\mathrm{2}}\right)+\lambda \left({\mathbf{∊}}_{\mathrm{11}}{\mathbf{∊}}_{\mathrm{22}}+{\mathbf{∊}}_{\mathrm{11}}{\mathbf{∊}}_{\mathrm{33}}+{\mathbf{∊}}_{\mathrm{22}}{\mathbf{∊}}_{\mathrm{33}}\right)+\mathrm{2}G\left({\mathbf{∊}}_{\mathrm{12}}^{\mathrm{2}}+{\mathbf{∊}}_{\mathrm{13}}^{\mathrm{2}}+{\mathbf{∊}}_{\mathrm{23}}^{\mathrm{2}}\right)dV={\int }_{\Omega }\hspace{0.17em}\left|\frac{\partial {\mathbf{r}}_o}{\partial \mathbf{\xi }}\right|[\frac{\lambda +\mathrm{2}G}{\mathrm{2}}\left({\mathbf{∊}}_{\mathrm{11}}^{\mathrm{2}}+{\mathbf{∊}}_{\mathrm{22}}^{\mathrm{2}}+{\mathbf{∊}}_{\mathrm{33}}^{\mathrm{2}}\right)\phantom{\text{hhhhhhhh}}+\lambda \left({\mathbf{∊}}_{\mathrm{11}}{\mathbf{∊}}_{\mathrm{22}}+{\mathbf{∊}}_{\mathrm{11}}{\mathbf{∊}}_{\mathrm{33}}+{\mathbf{∊}}_{\mathrm{22}}{\mathbf{∊}}_{\mathrm{33}}\right)\phantom{\text{hhhhhhhh}}+\mathrm{2}G\left({\mathbf{∊}}_{\mathrm{12}}^{\mathrm{2}}+{\mathbf{∊}}_{\mathrm{13}}^{\mathrm{2}}+{\mathbf{∊}}_{\mathrm{23}}^{\mathrm{2}}\right)\phantom{\frac{\lambda +\mathrm{2}G}{\mathrm{2}}}]d\Omega ,\end{array}$ (17)

where λ is Lame's constant and G is the shear modulus. The vector of elastic forces is defined as the derivative of U _s with respect to generalized coordinates e:

${\mathbf{Q}}_s=\frac{\partial U_s}{\partial \mathbf{e}}={\int }_V\hspace{0.17em}\frac{\lambda +\mathrm{2}G}{\mathrm{2}}\left(\mathrm{2}{\mathbf{∊}}_{\mathrm{11}}\frac{\partial {\mathbf{∊}}_{\mathrm{11}}}{\partial \mathbf{e}}+\mathrm{2}{\mathbf{∊}}_{\mathrm{22}}\frac{\partial {\mathbf{∊}}_{\mathrm{22}}}{\partial \mathbf{e}}+\mathrm{2}{\mathbf{∊}}_{\mathrm{33}}\frac{\partial {\mathbf{∊}}_{\mathrm{33}}}{\partial \mathbf{e}}\right)\hspace{1em}+\lambda ({\mathbf{∊}}_{\mathrm{11}}\frac{\partial {\mathbf{∊}}_{\mathrm{22}}}{\partial \mathbf{e}}+\frac{\partial {\mathbf{∊}}_{\mathrm{11}}}{\partial \mathbf{e}}{\mathbf{∊}}_{\mathrm{22}}+{\mathbf{∊}}_{\mathrm{11}}\frac{\partial {\mathbf{∊}}_{\mathrm{33}}}{\partial \mathbf{e}}+\frac{\partial {\mathbf{∊}}_{\mathrm{11}}}{\partial \mathbf{e}}{\mathbf{∊}}_{\mathrm{33}}\phantom{\text{hhhhh}}+{\mathbf{∊}}_{\mathrm{22}}\frac{\partial {\mathbf{∊}}_{\mathrm{33}}}{\partial \mathbf{e}}+\frac{\partial {\mathbf{∊}}_{\mathrm{22}}}{\partial \mathbf{e}}{\mathbf{∊}}_{\mathrm{33}})\hspace{1em}+\mathrm{2}G\left(\mathrm{2}{\mathbf{∊}}_{\mathrm{12}}\frac{\partial {\mathbf{∊}}_{\mathrm{12}}}{\partial \mathbf{e}}+\mathrm{2}{\mathbf{∊}}_{\mathrm{13}}\frac{\partial {\mathbf{∊}}_{\mathrm{13}}}{\partial \mathbf{e}}+\mathrm{2}{\mathbf{∊}}_{\mathrm{23}}\frac{\partial {\mathbf{∊}}_{\mathrm{23}}}{\partial \mathbf{e}}\right)dV={\int }_{\Omega }\hspace{0.17em}\left|\frac{\partial {\mathbf{r}}_o}{\partial \mathbf{\xi }}\right|[\frac{\lambda +\mathrm{2}G}{\mathrm{2}}\left(\mathrm{2}{\mathbf{∊}}_{\mathrm{11}}\frac{\partial {\mathbf{∊}}_{\mathrm{11}}}{\partial \mathbf{e}}+\mathrm{2}{\mathbf{∊}}_{\mathrm{22}}\frac{\partial {\mathbf{∊}}_{\mathrm{22}}}{\partial \mathbf{e}}+\mathrm{2}{\mathbf{∊}}_{\mathrm{33}}\frac{\partial {\mathbf{∊}}_{\mathrm{33}}}{\partial \mathbf{e}}\right)\phantom{\text{hhhhhhhh}}+\lambda ({\mathbf{∊}}_{\mathrm{11}}\frac{\partial {\mathbf{∊}}_{\mathrm{22}}}{\partial \mathbf{e}}+\frac{\partial {\mathbf{∊}}_{\mathrm{11}}}{\partial \mathbf{e}}{\mathbf{∊}}_{\mathrm{22}}+{\mathbf{∊}}_{\mathrm{11}}\frac{\partial {\mathbf{∊}}_{\mathrm{33}}}{\partial \mathbf{e}}\phantom{\text{hhhhhhhhhhhh}}+\frac{\partial {\mathbf{∊}}_{\mathrm{11}}}{\partial \mathbf{e}}{\mathbf{∊}}_{\mathrm{33}}+{\mathbf{∊}}_{\mathrm{22}}\frac{\partial {\mathbf{∊}}_{\mathrm{33}}}{\partial \mathbf{e}}+\frac{\partial {\mathbf{∊}}_{\mathrm{22}}}{\partial \mathbf{e}}{\mathbf{∊}}_{\mathrm{33}})\phantom{\text{hhhhhhhh}}+\mathrm{2}G(\mathrm{2}{\mathbf{∊}}_{\mathrm{12}}\frac{\partial {\mathbf{∊}}_{\mathrm{12}}}{\partial \mathbf{e}}+\mathrm{2}{\mathbf{∊}}_{\mathrm{13}}\frac{\partial {\mathbf{∊}}_{\mathrm{13}}}{\partial \mathbf{e}}\phantom{\text{hhhhhhhhhhhhh}}+\mathrm{2}{\mathbf{∊}}_{\mathrm{23}}\frac{\partial {\mathbf{∊}}_{\mathrm{23}}}{\partial \mathbf{e}})]d\Omega ,$ (18)

where

$\frac{\partial {∊}_{ij}}{\partial \mathbf{e}}=\frac{\mathrm{1}}{\mathrm{2}}\left({\mathbf{S}}_{,i}^T{\mathbf{S}}_{,j}+{\mathbf{S}}_{,j}^T{\mathbf{S}}_{,i}\right)\mathbf{e}\mathbf{,}\hspace{1em}i,j=\mathrm{1,2},\mathrm{3},\hspace{0.17em}\hspace{0.17em}j\ge i.$ (19)

It can be seen from (14), (16), (18), and (19) that S_{i, ξ j} is the key to evaluate the elastic forces. S_{i, ξ j} 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:

$\frac{\partial {\mathbf{Q}}_s}{\partial \mathbf{e}}=\sum _{\alpha =\mathrm{1}}^{\mathrm{3}}\hspace{0.17em}\frac{\lambda +\mathrm{2}G}{\mathrm{2}}{\int }_V\hspace{0.17em}(\left({\mathbf{e}}^T{\mathbf{S}}_{,\alpha }^T{\mathbf{S}}_{,\alpha }\mathbf{e}\right){\mathbf{S}}_{,\alpha }^T{\mathbf{S}}_{,\alpha }\phantom{\text{hhhhhhhhhhhhhh}}+\mathrm{2}{\mathbf{S}}_{,\alpha }^T{\mathbf{S}}_{,\alpha }\mathbf{e}{\mathbf{e}}^T{\mathbf{S}}_{,\alpha }^T{\mathbf{S}}_{,\alpha }-{\mathbf{S}}_{,\alpha }^T{\mathbf{S}}_{,\alpha })dV\hspace{1em}+\sum _{\alpha =\mathrm{1}}^{\mathrm{3}}\hspace{0.17em}\sum _{\begin{array}{c}\beta =\mathrm{1}\\ \beta \ne \alpha \end{array}}^{\mathrm{3}}\frac{\lambda }{\mathrm{2}}{\int }_V(\left({\mathbf{e}}^T{\mathbf{S}}_{,\beta }^T{\mathbf{S}}_{,\beta }\mathbf{e}\right){\mathbf{S}}_{,\alpha }^T{\mathbf{S}}_{,\alpha }\hspace{0.17em}\phantom{hhhhhhhhhhhhhh}+\mathrm{2}{\mathbf{S}}_{,\alpha }^T{\mathbf{S}}_{,\alpha }\mathbf{e}{\mathbf{e}}^T{\mathbf{S}}_{,\beta }^T{\mathbf{S}}_{,\beta }-{\mathbf{S}}_{,\alpha }^T{\mathbf{S}}_{,\alpha })dV\hspace{1em}+\sum _{\alpha =\mathrm{1}}^{\mathrm{3}}\hspace{0.17em}\sum _{\begin{array}{c}\beta =\mathrm{1}\\ \beta \ne \alpha \end{array}}^{\mathrm{3}}G{\int }_V\hspace{0.17em}(\left({\mathbf{e}}^T{\mathbf{S}}_{,\alpha }^T{\mathbf{S}}_{,\beta }\mathbf{e}\right){\mathbf{S}}_{,\alpha }^T{\mathbf{S}}_{,\beta }\phantom{hhhhhhhhhhhhhh}+{\mathbf{S}}_{,\alpha }^T{\mathbf{S}}_{,\beta }\mathbf{e}{\mathbf{e}}^T\left({\mathbf{S}}_{,\alpha }^T{\mathbf{S}}_{,\beta }+{\mathbf{S}}_{,\beta }^T{\mathbf{S}}_{,\alpha }\right))dV=\sum _{\alpha =\mathrm{1}}^{\mathrm{3}}\hspace{0.17em}\frac{\lambda +\mathrm{2}G}{\mathrm{2}}\left|\frac{\partial {\mathbf{r}}_o}{\partial \mathbf{\xi }}\right|{\int }_{\Omega }(\left({\mathbf{e}}^T{\mathbf{S}}_{,\alpha }^T{\mathbf{S}}_{,\alpha }\mathbf{e}\right){\mathbf{S}}_{,\alpha }^T{\mathbf{S}}_{,\alpha }\hspace{0.17em}\phantom{hhhhhhhhhhhhhh}\phantom{hhhh}+\mathrm{2}{\mathbf{S}}_{,\alpha }^T{\mathbf{S}}_{,\alpha }\mathbf{e}{\mathbf{e}}^T{\mathbf{S}}_{,\alpha }^T{\mathbf{S}}_{,\alpha }-{\mathbf{S}}_{,\alpha }^T{\mathbf{S}}_{,\alpha })d\Omega \hspace{1em}+\sum _{\alpha =\mathrm{1}}^{\mathrm{3}}\hspace{0.17em}\sum _{\begin{array}{c}\beta =\mathrm{1}\\ \beta \ne \alpha \end{array}}^{\mathrm{3}}\frac{\lambda }{\mathrm{2}}\left|\frac{\partial {\mathbf{r}}_o}{\partial \mathbf{\xi }}\right|{\int }_{\Omega }(\left({\mathbf{e}}^T{\mathbf{S}}_{,\beta }^T{\mathbf{S}}_{,\beta }\mathbf{e}\right){\mathbf{S}}_{,\alpha }^T{\mathbf{S}}_{,\alpha }\hspace{0.17em}\phantom{hhhhhhhhhhhhhhhhhh}+\mathrm{2}{\mathbf{S}}_{,\alpha }^T{\mathbf{S}}_{,\alpha }\mathbf{e}{\mathbf{e}}^T{\mathbf{S}}_{,\beta }^T{\mathbf{S}}_{,\beta }-{\mathbf{S}}_{,\alpha }^T{\mathbf{S}}_{,\alpha })d\Omega \hspace{1em}+\sum _{\alpha =\mathrm{1}}^{\mathrm{3}}\hspace{0.17em}\sum _{\begin{array}{c}\beta =\mathrm{1}\\ \beta \ne \alpha \end{array}}^{\mathrm{3}}G\left|\frac{\partial {\mathbf{r}}_o}{\partial \mathbf{\xi }}\right|{\int }_{\Omega }(\left({\mathbf{e}}^T{\mathbf{S}}_{,\alpha }^T{\mathbf{S}}_{,\beta }\mathbf{e}\right){\mathbf{S}}_{,\alpha }^T{\mathbf{S}}_{,\beta }\hspace{0.17em}\phantom{hhhhhhhhhhhhhhhhhh}+{\mathbf{S}}_{,\alpha }^T{\mathbf{S}}_{,\beta }\mathbf{e}{\mathbf{e}}^T\left({\mathbf{S}}_{,\alpha }^T{\mathbf{S}}_{,\beta }+{\mathbf{S}}_{,\beta }^T{\mathbf{S}}_{,\alpha }\right))d\Omega .$ (20)

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, E _Y , and Poisson's ratio, μ. The relationships between these parameters and Lame's constants are the following:

$\begin{array}{c}\lambda =\frac{\mu E_Y}{\left(\mathrm{1}+\mu \right)\left(\mathrm{1}-\mathrm{2}\mu \right)},\\ G=\frac{E_Y}{\mathrm{2}\left(\mathrm{1}+\mu \right)}.\end{array}$ (21)

3.3. Body Forces

Taking gravity forces as an example, the virtual work of body forces can be written as

${\int }_V\hspace{0.17em}{\mathbf{f}}_b·\delta \mathbf{r}dV={\int }_V\hspace{0.17em}\left[\begin{array}{ccc}0& 0& -\rho g\end{array}\right]\mathbf{S}dV\delta \mathbf{e}=\left[\begin{array}{ccc}0& 0& -\rho g\end{array}\right]\left|\frac{\partial {\mathbf{r}}_o}{\partial \mathbf{\xi }}\right|{\int }_{\Omega }\hspace{0.17em}\mathbf{S}d\Omega \delta \mathbf{e}.$ (22)

It can be inferred that the vector of gravity forces is

${\mathbf{Q}}_g=\left|\frac{\partial {\mathbf{r}}_o}{\partial \mathbf{\xi }}\right|{\int }_{\Omega }\hspace{0.17em}{\mathbf{S}}^{T}d\Omega \left[\begin{array}{c}0\\ 0\\ -\rho g\end{array}\right].$ (23)

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

$\begin{array}{c}\mathbf{M}\left(\mathbf{e}\right)\stackrel{‥}{\mathbf{e}}+{\mathrm{\Phi }}_{\mathbf{e}}^T\left(\mathbf{e},t\right)\mathbf{\lambda }-\mathbf{Q}\left(t,\mathbf{e},\dot{\mathbf{e}}\right)=\mathbf{0},\\ \mathrm{\Phi }\left(\mathbf{e},t\right)=\mathbf{0},\end{array}$ (24)

where M is the system mass matrix, e, $\dot{\mathbf{e}}$ , and $\stackrel{‥}{\mathbf{e}}$ are the position coordinates, velocities, and accelerations of control points in the global coordinate system, t is the time variable, Φ are the system constraint equations, Φ _e is the Jacobian of constraints Φ, λ is the vector of the Lagrange multipliers, and Q 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 n + 1 are discretized as

$\begin{array}{c}\mathbf{M}\left({\mathbf{e}}_{n+\mathrm{1}}\right){\stackrel{‥}{\mathbf{e}}}_{n+\mathrm{1}}+{\left({\mathrm{\Phi }}_{\mathbf{e}}^T\mathbf{\lambda }\right)}_{n+\mathrm{1}}-{\mathbf{Q}}_{n+\mathrm{1}}=\mathbf{0},\\ \frac{\mathrm{1}-{\alpha }_m}{\left(\mathrm{1}-{\alpha }_f\right)\beta h^{\mathrm{2}}}\mathrm{\Phi }\left({\mathbf{e}}_{n+\mathrm{1}},t_{n+\mathrm{1}}\right)=\mathbf{0},\end{array}$ (25)

where h is the time step, α _m = (2ρ_∞ – 1)/(ρ_∞ + 1), α _f = (ρ_∞)/(ρ_∞ + 1), ρ_∞ ∊ [0, 1], and β = (1 + α _f – α _m )²/4. ρ_∞ is a parameter controlling numerical dissipation (ρ_∞ = 0 for maximal dissipation). The position and velocity variables of time step n + 1 are discretized as

$\begin{array}{c}{\mathbf{e}}_{n+\mathrm{1}}={\mathbf{e}}_n+h{\dot{\mathbf{e}}}_n+\frac{h^{\mathrm{2}}}{\mathrm{2}}\left[\left(\mathrm{1}-\mathrm{2}\beta \right){\mathbf{q}}_n+\mathrm{2}\beta {\mathbf{q}}_{n+\mathrm{1}}\right],\\ {\dot{\mathbf{e}}}_{n+\mathrm{1}}={\dot{\mathbf{e}}}_n+h\left[\left(\mathrm{1}-\gamma \right){\mathbf{q}}_n+\gamma {\mathbf{q}}_{n+\mathrm{1}}\right],\end{array}$ (26)

where γ = 1/2 + α _f – α _m ; q _n is constituted by the acceleration-like auxiliary variables, which require that

$\begin{array}{c}\left(\mathrm{1}-{\alpha }_m\right){\mathbf{q}}_{n+\mathrm{1}}+{\alpha }_m{\mathbf{q}}_n=\left(\mathrm{1}-{\alpha }_f\right){\stackrel{‥}{\mathbf{e}}}_{n+\mathrm{1}}+{\alpha }_f{\stackrel{‥}{\mathbf{e}}}_n,\\ {\mathbf{q}}_{\mathrm{0}}={\stackrel{‥}{\mathbf{e}}}_{\mathrm{0}}.\end{array}$ (27)

Algorithm 1 sums the process of generalized-α scheme to solve dynamic equilibrium at time t_{n + 1}.

Algorithm 1:

Generalized-α scheme.

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

${\mathbf{J}}_{\text{Gen-}\alpha }\hspace{1em}=\left[\begin{array}{cc}\mathbf{M}+\beta h^{\mathrm{2}}\frac{\mathrm{1}-{\alpha }_f}{\mathrm{1}-{\alpha }_m}{\left(\mathbf{M}\left(\mathbf{e}\right)\stackrel{‥}{\mathbf{e}}+{\mathrm{\Phi }}_{\mathbf{e}}^T\mathbf{\lambda }-\mathbf{Q}\right)}_{\mathbf{e}}-\gamma h\frac{\mathrm{1}-{\alpha }_f}{\mathrm{1}-{\alpha }_m}{\mathbf{Q}}_{\dot{\mathbf{e}}}& {\mathrm{\Phi }}_{\mathbf{e}}^T\\ {\mathrm{\Phi }}_{\mathbf{e}}& \mathit{0}\end{array}\right].$ (28)

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 3 m, and its thickness is chosen to be 0.1 m. The density, Young's modulus, and Poisson ratio of the material are assumed to be 5000 kg/cm³, 1. 0 × 10⁷ N/m², and 0.3, respectively.

Figure 1:

The configuration of the thin plate at different moments (36 elements).

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

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, Ω _x = Ω _y = {0, 0, 0, 0, 1, 1, 1, 1}, and Ω _z = {0, 0, 1, 1}. 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.

Figure 2:

Relationship between the element and the corresponding control points.

To see the influence of refinement on the plate's deformation, other 36-elements (Ω _x = Ω _y = {0, 0, 0, 0, 1/6, 1/3, 1/2, 2/3, 5/6, 1, 1, 1, 1}, Ω _z = {0, 0, 1, 1}) and 100-element (Ω _x = Ω _y = {0, 0, 0, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 1, 1, 1}, Ω _z = {0, 0, 1, 1}) meshing plans are performed separately. A 2 s simulation with 2⁻¹² 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 36-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, 36-element and 100-element meshing plans at different moments. It is obvious that the plate configuration changes significantly when the number of elements increases from 1 to 36, while it differs not that much from 36 to 100. 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.

Figure 3:

The configuration of the plate and related control points at different moments (1 element).

Figure 4:

Comparison of plate deformation of different meshing plans (1 element: yellow; 36 elements: cyan; 100 elements: red).

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.