Invariant isogeometric formulation for the geometric stiffness matrix of spatial curved Kirchhoff rods

doi:10.1016/j.cma.2021.113692

Computer Methods in Applied Mechanics and Engineering

Volume 377, 15 April 2021, 113692

https://doi.org/10.1016/j.cma.2021.113692 Get rights and content

Abstract

This paper presents an invariant isogeometric formulation for the geometric stiffness matrix of spatial curved Kirchhoff rods considering various end moments, i.e., the internal (member) moments and applied (conservative) moments. There are two levels of rigid-body qualification, one is on the buckling theory of the rod itself and the other on the isogeometric formulation for discretization. Both will be illustrated. Based on the updated Lagrangian formulation of three-dimensional continua, the rotational effect of end moments is naturally included in the external virtual work done by end tractions without introducing any definition of finite rotations. Both the geometric torsion and curvatures of the rod are considered closely for the centroidal axis, except with the omission of higher order terms. The geometric stiffness matrix for internal moments is consistent with that of the geometrically exact rod model with its rigid-body quality demonstrated. For structures rigorously defined for the deformed state, the geometric stiffness matrix after global assembly is always symmetric, for both the internal and external moments. By adopting the invariant isogeometric discretization following our previous work, a series of numerical examples, including the cases of external conservative moments, angled joint and complicated spatial geometry, were solved for buckling analysis, by which the reliability of the geometric stiffness matrix derived is verified via comparison with the analytical or straight beam solutions.

Introduction

General spatial Kirchhoff rod models have received considerable attention in the literature, due to their numerical advantages of, but not limited to, possibility of rotation-free formulations with fewer degrees of freedom, avoidance of shear-locking, and smoother approximation of the rod geometry. Previously, the finite element formulations of geometrically exact Kirchhoff rods were presented by Weiss [1], [2] and Boyer and Primault [3], where only initially straight rods with circular cross-section were considered. Researches that have been carried out achieve better numerical qualities including the applicability to complicated initial geometry [4], [5], [6], [7], [8], objectivity and path-independence [6], [9], avoidance of membrane-locking [10], [11], rotation-free discretization [4], [7], [12], [13], and singularity-free formulation [14]. For nonlinear analysis, few of these works have examined the quality of the geometric stiffness matrix with regard to rigid-body rotations. It was noted that the use of geometric stiffness matrix not qualified for rigid-body rotations may lead to failure of convergence in nonlinear analysis or incorrect solutions in buckling analysis.

Among the above nonlinear Kirchhoff rod models, despite the discretization strategies adopted, the local geometric stiffness matrices derived are generally different when using the equation of virtual work or the weak form [2], [6], [9], [13], or the so-called internal energy of the element considered [3], [7], [14]. In the latter approach, the internal energy of an element is directly obtained from the equation of virtual work, by which a symmetric local geometric stiffness matrix was obtained, not aware of the non-conjugateness in energy between the end moments and rotations expressed as the derivatives of displacements. A detailed review of the works conducted mainly by Yang and co-workers on the aforementioned problem for straight and curved beams is available in [15].

On the other hand, linearization of the equation of virtual work or the weak form with the finite rotation represented by an orthogonal transformation, e.g., the Rodrigues’ formula in the geometrically exact rod model, as was done by Simo and Vu-Quoc [16], leads to an asymmetric local geometric stiffness matrix with the rotational effect of end moments automatically taken into account. It was shown that the global matrix is symmetric since the anti-symmetric parts of the local matrices will cancel out in the assembly of the global matrix for conservative structures [16], [17], [18]. However, the geometric stiffness matrix of this kind is only valid for the internal (member) moments interpreted as stress resultants [17], [18], as special treatment is required for external (applied) moments generated by conservative force mechanisms.

A proper treatment of the 3D rotational degrees of freedom of a rod is to adopt work-conjugate definitions for end rotations and moments. Argyris et al. [19] treated the end moments and rotations both to be of the semi-tangential type. Through this definition, a symmetric geometric stiffness matrix was obtained for straight beams for the nonlinear analysis of rigid-jointed frames. The same definition was also adopted in the nonlinear formulation of straight I-beams including the warping deformation [20]. Yang and Kuo [17], [18] derived the local geometric stiffness matrix by treating the end rotations as the derivatives of displacements, a definition commonly used by engineers. From the result of derivation, the bending moments and torque are identified as quasi- and semi-tangential, respectively, and the local geometric stiffness matrix is asymmetric. They also showed that the geometric stiffness matrix derived for internal moments can pass the rigid-body test [21]. Further, by enforcing the joint equilibrium at the deformed state, they proved that the global geometric stiffness matrix is symmetric for the structure [17], [18]. Other work-conjugate definitions of rotations [22], [23], [24] have been proposed in the literature, and the global geometric stiffness matrices for internal moments and external conservative moments are always symmetric for the structure rigorously defined for the deformed state. In fact, the definition of rotations can be artificially prescribed; it is known that identical results can be obtained if energy conjugateness is obeyed by each approach for the end quantities adopted. Overall, the aforementioned approaches are limited to beam elements that are straight or circular in geometry. Thus the consistency of the geometric stiffness matrix for the finite elements derived by the geometrically exact rod model (e.g., [16]) and that by the straight-beam approach (e.g., [17], [18], [19]) was not fully assessed in the literature.

Recently, the isogeometric analysis (IGA) proposed by Hughes et al. [25] provides an alternative technique for simulating spatial Kirchhoff rods with several advantages over the conventional finite element analysis (FEA), including the simplification in mesh generation process, exact representation of geometry and high-order continuity of the shape functions, etc. Although considerable researches (e.g., [4], [5], [12], [26]) have been conducted for the linear analysis of spatial Kirchhoff rods, only few were focused on the nonlinear analysis. Raknes et al. [27] proposed a nonlinear isogeometric formulation for bending-stabilized cables, where the torsional deformation was neglected. Bauer et al. [7] extended the natural reference frame (or the smallest-rotation reference frame) and rotation angle representation in [4] for the rod’s configuration to nonlinear analysis, including all the deformations of the rod. However, the element formulation in [7] is not fully objective as stated in Meier et al. [28], which may cause problems for certain rigid-body modes. This is not due to non-collinearly jointed rod structures but the application of the smallest-rotation triads.

Motivated by the above, the objective of this paper is to extend the invariant isogeometric formulation for linear analysis [26] to deriving the geometric stiffness matrix of general spatial Kirchhoff rods considering internal (member) moments and external (applied) conservative moments at ends. In the pre-buckling analysis, the displacements are assumed to be such small that the classical angle of twist and rotations about the principal axes can be adopted. The rod model is based on the updated Lagrangian formulation of 3D continua [29], through which the effect of rotational degrees of freedom can be captured by the external virtual work due to end tractions. The symmetry of the global geometric stiffness matrices for various end moments will be discussed. In this paper, the rigid-body qualification will be assessed for two levels, one is on the buckling theory of the rod itself and the other on the isogeometric formulation for discretization of the rod.

The outline of the paper is as follows. After a review of the geometry and kinematics of spatial Kirchhoff rods in Section 2, we present the equations of equilibrium of stress resultants in Section 3. In Section 4, the principle of virtual work for the rod is established. In Section 5, the geometric stiffness matrices for internal and external conservative moments acting at ends are derived, with discussion on their symmetry and qualification by the rigid-body test. Section 6 summarizes the isogeometric discretization that is invariant for spatial Kirchhoff rods. Finally, a series of representative numerical examples are implemented to examine the quality of the geometric stiffness matrices presented.

Section snippets

Geometry and kinematics description of the rod

As shown in Fig. 1, the initial undeformed, the last known deformed and the current deformed configurations are denoted by $C_{0}$ , $C_{1}$ and $C_{2}$ , respectively. Bathe’s tensor notation [29] is adopted such that left superscripts indicate the configuration where the quantity occurs and left subscripts the configuration of reference. In this study, we adopt the updated Lagrangian formulation with reference to $C_{1}$ configuration, and the left subscripts are omitted for conciseness.

Stress resultants of the rod

In terms of the second Piola–Kirchhoff stresses in the Cartesian coordinates ${^{t} d_{i}}$ , the axial force $^{1} N_{1}$ and shear forces $^{t} N_{α}$ of the rod can be defined as $^{t} N_{1} ≔ \int_{A}^{t} {\tilde{S}}^{11} d A,^{t} N_{2} ≔ \int_{A}^{t} {\tilde{S}}^{12} d A,^{t} N_{3} ≔ \int_{A}^{t} {\tilde{S}}^{13} d A$ and the torque $^{t} M_{1}$ and bending moments $^{1} M_{α}$ as $^{t} M_{1} ≔ \int_{A} (^{t} {\tilde{S}}^{13} ξ_{2} -^{t} {\tilde{S}}^{12} ξ_{3}) d A,^{t} M_{2} ≔ \int_{A}^{t} {\tilde{S}}^{11} ξ_{3} d A,^{t} M_{3} ≔ - \int_{A}^{t} {\tilde{S}}^{11} ξ_{2} d A$ They may also be expressed in terms of the force and moment vectors $^{t} n$ and $^{t} m$ , respectively, as $^{t} n ≔^{t} N_{i}^{t} d_{i},^{t} m ≔^{t} M_{i}^{t} d_{i}$ Moreover, for a bisymmetric cross section, we have the following condition: $\int_{A}^{t} {\tilde{S}}^{12} ξ_{2} d A =$

Principle of virtual work

Based on the updated Lagrangian formulation of three-dimensional continua, the principle of virtual work for spatial Kirchhoff rods will be established by rigorous derivations of the virtual strain energy due to linear strains $δ U$ , the virtual potential energy due to initial stresses along the rod $δ V$ and the external virtual work $δ W$ . We then show that the total virtual potential energy due to initial stresses $δ \bar{V}$ , from which the local geometric stiffness matrix is obtained, is different from $δ V$

Geometric stiffness matrices

In this section, we will derive the total virtual potential energies due to initial stresses $δ \bar{V}$ or the geometric stiffness matrices for both the internal (member) moments and external (applied) conservative moments at the ends.

Isogeometric formulation

In this section, we use the isogeometric discretization proposed in [26] for the buckling analysis of spatial Kirchhoff rods for the following advantages compared with the finite element approaches: the geometry of the rod that significantly affects the accuracy of buckling analysis is exactly represented; higher-order continuity of the rod’s kinematics can be achieved, since the NURBS shape functions guarantee $C^{p - 1}$ continuity (where $p$ is the polynomial order); the isogeometric discretization

Representative numerical examples

In this section, we focus on the buckling analysis of a series of curved rod structures, to examine the quality of the rod theory presented, including particularly the geometric stiffness matrices, and the invariant isogeometric formulation for discretization. After examining the classical circular arch problems under pure bending and uniform compression, we will study the lateral buckling of two planar arches considering external (applied) conservative moments and angled joints, respectively.

Conclusions

Based on the updated Lagrangian formulation of three-dimensional continua, the total virtual potential energy due to initial stresses or the geometric stiffness matrix of spatial curved Kirchhoff rods is newly derived by considering the internal moments (as stress resultants) and external moments (generated by three conservative mechanisms) acting at the ends. Both the geometric torsion and curvatures of the rod are considered. In the virtual work formulation, we observe that the virtual

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

Endowment of the Fengtay Chair Professorship to the senior author is greatly appreciated. Agencies that have sponsored this research include: National Natural Science Foundation of China (Grant No. 51678091), Special Program for Scientific and Technological Talents (Grant No. cstc2018jcyj-yszxX0013), Chongqing Municipal Natural Science Foundation (Grant No. cstc2019yszx-jcyjX0001), Chongqing Municipal Natural Science Foundation (Grant No. cstc2017zdcy-yszxX0006), and Science, Technology

References (35)

GrecoL. et al.
B-spline interpolation of kirchhoff–love space rods
Comput. Methods Appl. Mech. Engrg.
(2013)
GrecoL. et al.
An implicit G1 multi patch B-spline interpolation for kirchhoff–love space rod
Comput. Methods Appl. Mech. Engrg.
(2014)
MeierC. et al.
An objective 3D large deformation finite element formulation for geometrically exact curved kirchhoff rods
Comput. Methods Appl. Mech. Engrg.
(2014)
BauerA.M. et al.
Nonlinear isogeometric spatial Bernoulli beam
Comput. Methods Appl. Mech. Engrg.
(2016)
MeierC. et al.
A locking-free finite element formulation and reduced models for geometrically exact kirchhoff rods
Comput. Methods Appl. Mech. Engrg.
(2015)
GrecoL. et al.
An efficient blended mixed B-spline formulation for removing membrane locking in plane curved kirchhoff rods
Comput. Methods Appl. Mech. Engrg.
(2017)
RadenkovićG. et al.
Linear static isogeometric analysis of an arbitrarily curved spatial Bernoulli–Euler beam
Comput. Methods Appl. Mech. Engrg.
(2018)
YangY.B. et al.
Rational nonlinear analysis of framed structures and curved beams considering joint equilibrium in deformed state
Int. J. NonLinear. Mech.
(2020)
SimoJ.C. et al.
A three-dimensional finite-strain rod model. part II: Computational aspects
Comput. Methods Appl. Mech. Engrg.
(1986)
ArgyrisJ.H. et al.
On large displacement-small strain analysis of structures with rotational degrees of freedom
Comput. Methods Appl. Mech. Engrg.
(1978)

BesselingJ.F.

Derivatives of deformation parameters for bar elements and their use in buckling and postbuckling analysis

Comput. Methods Appl. Mech. Engrg.

(1977)

IzzuddinB.A.

Conceptual issues in geometrically nonlinear analysis of 3D framed structures

Comput. Methods Appl. Mech. Engrg.

(2001)

HughesT.J.R. et al.

Isogeometric analysis: CAD finite elements NURBS exact geometry and mesh refinement

Comput. Methods Appl. Mech. Engrg.

(2005)

RaknesS.B. et al.

Isogeometric rotation-free bending-stabilized cables: Statics, dynamics, bending strips and coupling with shells

Comput. Methods Appl. Mech. Eng.

(2013)

ArmeroF. et al.

Invariant hermitian finite elements for thin kirchhoff rods. II: The linear three-dimensional case

Comput. Methods Appl. Mech. Engrg.

(2012)

WeissH.

Dynamics of geometrically nonlinear rods: I. Mechanical models and equations of motion

Nonlinear Dynam.

(2002)

WeissH.

Dynamics of geometrically nonlinear rods: II. Numerical methods and computational examples

Nonlinear Dynam.

(2002)

Cited by (15)

Corotational isogeometric shear deformable geometrically exact spatial form beam element for general large deformation analysis of flexible thin-walled beam structures
2024, Thin-Walled Structures
In this paper, a new 3D generic order corotational geometrically exact beam element undergoing arbitrarily large deformations is presented. The novelty of the formulation lies in the use of the standard corotational framework, that is the decomposition into rigid body motion and pure deformation, to derive not only the objective, path-independent and singularity-free spatial exact strain measure but also the internal force and the consistent tangent stiffness matrix. The new formulation is derived by firstly conducting the consistent linearization of the weak formulation for the equilibrium equations of Reissner–Simo theory in the spatial form. By constructing a co-rotated frame that continuously rotates and translates with the element, the nodal strain-producing deformations can be extracted with the rigid-body translations and rotations removed from the element total motion. Subsequently, the corotational deformation variables are exploited to establish the distribution of deformation field within the beam element through B-spline basis functions, and derive the invariant and path independent exact strain measure in the spatial form by Lie derivative formalism. Through performing the variational operations of the compact transformation between corotational and global variables, the consistent tangent stiffness matrix are derived in closed form of the discrete variables for deformation field of any order. The proposed formulation is applicable to not only low-order two-node geometrically exact elements but also high-order multi-node straight or curved ones. The accuracy and robustness of the new formulation are demonstrated through the solutions of various benchmark problems.
Geometrically exact isogeometric Bernoulli–Euler beam based on the Frenet–Serret frame
2023, Computer Methods in Applied Mechanics and Engineering
A novel geometrically exact model of the spatially curved Bernoulli–Euler beam is developed. The formulation utilizes the Frenet–Serret frame as the reference for updating the orientation of a cross section. The weak form is consistently derived and linearized, including the contributions from kinematic constraints and configuration-dependent load. The nonlinear terms with respect to the cross-sectional coordinates are strictly considered, and the obtained constitutive model is scrutinized. The main features of the formulation are invariance with respect to the rigid-body motion, path-independence, and improved accuracy for strongly curved beams. A new reduced beam model is conceived as a special case, by omitting the rotational DOF. Although rotation-free, the reduced model includes the torsion of the beam axis, which allows simulations of spatial beams that are predominantly bent with respect to the binormal. The applicability of the obtained isogeometric finite element is verified via a set of standard academic benchmark examples. The formulation is able to accurately model strongly curved Bernoulli–Euler beams that have well-defined Frenet–Serret frames.
A novel four-field mixed FE approximation for Kirchhoff rods using Cartan's moving frames
2022, Computer Methods in Applied Mechanics and Engineering
A four-field mixed variational principle is proposed for large deformation analysis of Kirchhoff rods with the lowest order $C^{0}$ mixed FE approximations. The core idea behind the approach is to introduce a one-parameter family of points (the centerline) and a separate one-parameter family of orthonormal frames (the Cartan moving frame) that are specified independently. The curvature and torsion of the curve are then related to the relative rotation of neighboring frames. The relationship between the frame and the centerline is then enforced at the solution step using a Lagrange multiplier (which plays the role of section force). In case of Kirchhoff rods, the cross sectional orientation can be described using frames like the Frenet–Serret, which are defined only using the centerline, thereby demanding higher-order smoothness for the centerline approximation. Decoupling the frame from the position vector of the base curve leads to a description of torsion and curvature that is independent of the position information, thus allowing for simpler interpolations. The four-field mixed variational principle we propose has the frame, section force, extension strain and position vector as input arguments. We discretize the position vector using linear Lagrange shape functions, while the frames are interpolated as piecewise geodesics on the rotation group. Finite element approximations for extensional strain and section force are constructed using constant shape functions. Using these discrete approximations, a discrete mixed variational principle is laid out which is then numerically extremized. The discrete approximation is then applied to a few benchmark problems. Vis-á-vis most available approaches, our numerical studies reveal an impressive performance of the proposed method without numerical instabilities or locking.
Geometrically exact static isogeometric analysis of an arbitrarily curved spatial Bernoulli–Euler beam
2022, Computer Methods in Applied Mechanics and Engineering
Citation Excerpt :
Another recent contribution to the analysis of assemblies of BE beams is developed in [45], where the authors have obtained the symmetric tangent stiffness matrix by the second covariant derivative of the internal energy. An invariant geometric stiffness matrix is derived in [46] considering various end moments. However, the formulation is applied solely to the buckling analysis.
The objective of this research is the development of a geometrically exact model for the analysis of arbitrarily curved spatial Bernoulli–Euler beams. The complete metric of the beam is utilized in order to include the effect of curviness on the nonlinear distribution of axial strain over the cross section. The exact constitutive relation between energetically conjugated pairs is employed, along with four reduced relations. The isogeometric approach, which allows smooth connections between finite elements, is used for the spatial discretization of the weak form. Two methods for updating the local vector basis are applied and discussed in the context of finite rotations. All the requirements of geometrically exact beam theory are satisfied, such as objectivity and path-independence. The accuracy of the formulation is verified by a thorough numerical analysis. The influence of the curviness on the structural response is scrutinized for two classic examples. If the exact response of the structure is sought, the curviness must be considered when choosing the appropriate beam model.
A non-linear symmetric G<sup>1</sup>-conforming Bézier finite element formulation for the analysis of Kirchhoff beam assemblies
2021, Computer Methods in Applied Mechanics and Engineering
Citation Excerpt :
Isogeometric formulations based on the Timoshenko beam model can be found in [11–14] and [15,16] respectively for the collocation and the Galerkin schemes. Recently in [23,24] an invariant isogeometric formulation for the Kirchhoff beam model, based on a suitable interpolation for the twist angle, has been presented for the linear case. An alternative Bernoulli beam model based on the Frenet–Serret triad is proposed in [25].
We present a formulation for slender space curved rods and rod assemblies that implicitly accounts for the Kirchhoff constraints and for the G $^{1}$ -continuity conditions (i.e. continuity of the geometric tangent) between elements. The whole formulation is developed in tensorial coordinate free form, apt to any numerical interpolation to be implemented. A symmetric tangent stiffness operator is obtained performing the second covariant derivative of the internal energy functional, for which the Levi-Civita connection of the configurations manifold of the rod is needed. The G $^{1}$ -continuity conditions are fulfilled by means of a change of basis, from the original configuration parameters (position of the beam axis and rotation around the axis tangent) to a new set of configuration parameters, whose relation to the original set is non-linear. For this reason an additional geometric term, specific for the G $^{1}$ -formulation, appears in the tangent stiffness matrix. The robustness and accuracy of the obtained Kirchhoff model is demonstrated with numerical examples that employ Bézier interpolation for the position and for the rotation angle.
Nonlinear statics of three-dimensional curved geometrically exact beams by a hierarchal quadrature element method
2024, International Journal for Numerical Methods in Engineering

View all citing articles on Scopus

View full text

Invariant isogeometric formulation for the geometric stiffness matrix of spatial curved Kirchhoff rods

Abstract

Introduction

Section snippets

Geometry and kinematics description of the rod

Stress resultants of the rod

Principle of virtual work

Geometric stiffness matrices

Isogeometric formulation

Representative numerical examples

Conclusions

Declaration of Competing Interest

Acknowledgments

Comput. Methods Appl. Mech. Engrg.

Comput. Methods Appl. Mech. Engrg.

Comput. Methods Appl. Mech. Engrg.

Comput. Methods Appl. Mech. Engrg.

Comput. Methods Appl. Mech. Engrg.

Comput. Methods Appl. Mech. Engrg.

Comput. Methods Appl. Mech. Engrg.

Int. J. NonLinear. Mech.

Comput. Methods Appl. Mech. Engrg.

Comput. Methods Appl. Mech. Engrg.

Comput. Methods Appl. Mech. Engrg.

Comput. Methods Appl. Mech. Engrg.

Comput. Methods Appl. Mech. Engrg.

Comput. Methods Appl. Mech. Eng.

Comput. Methods Appl. Mech. Engrg.

Dynamics of geometrically nonlinear rods: I. Mechanical models and equations of motion

Nonlinear Dynam.

Dynamics of geometrically nonlinear rods: II. Numerical methods and computational examples

Nonlinear Dynam.