Mechanical models in computational form finding of bending-active structures

Abstract

This article reviews the different aspects involved in computational form finding of bending-active structures based on the dynamic relaxation technique. Dynamic relaxation has been applied to form-finding problems of bending-active structures in a number of references. Due to the complex nature of large spatial deformations of flexible beams, the implementation of suitable mechanical beam models in the dynamic relaxation algorithm is a non-trivial task. Type of discretization and underlying beam theory have been identified as key aspects for numerical implementations. References can be classified into two groups depending on the selected discretization: finite-difference-like and finite-element-like. The first group includes 3- and 4-degree-of-freedom implementations based on increasingly complex beam models. The second gathers 6-degree-of-freedom discretizations based on co-rotational three-dimensional Kirchhoff–Love beam elements and geometrically exact Reissner–Simo beam elements. After reviewing and comparing implementation details, the advantages and drawbacks of each group have been discussed, and open aspects for future work have been pointed out.

Keywords

bending-active structures computational form finding dynamic relaxation Euler–Bernoulli beam theorem finite-difference-like discretization finite-element-like discretization geometrically exact beam model Kirchhoff–Love beam theory

Introduction

Bending-active structures are a kind of structures in which some members are pre-bent and then stabilized to achieve a certain desired configuration. The resulting system may reach considerable stiffness in relation to its weight due to the combination of the shape and slenderness of bent members and the post-bending stabilization by additional means as cables or membranes. The active bending principle has been used in vernacular architecture, for example, in the construction of tent-like dwellings (yurts or gers) by Asian nomadic peoples, where flexible timber slats are joined to form a lattice and bent to a cylindrical-shaped wall structure and a dome-shaped roof.

Since the construction of the Mannheim Multihalle¹ designed by Frei Otto, which is a pioneering modern application of the active bending principle in architecture, a number of dome-shaped grid-shells have been designed and built. Many of them have been designed as temporary or experimental structures as reported by Douthe et al.,² Nicholas et al.,³ Pone et al.⁴ or Harding et al.⁵ The basic principle of grid-shells is the use of very flexible members which are joined and bent into a target shape and stabilized by fixing the ground supports in a convenient way. The result is a lightweight shell-like structure which is stiff enough due to its shape, but whose members are light and easy to manipulate during construction. The arrangement of joints and the length and stiffness of members allows for different overall shapes of the dome.

Bending-active structures are not limited to grid-shells; other realizations include small-scale grandstand roofs,⁶ umbrellas and sculptural or ornamental applications. Lienhard et al.⁷ have prepared a comprehensive review of this structural type.

The design process is non-trivial, and involves three stages: (1) determination of the initial shape (form finding), (2) determination of the initial stress state, and (3) modelling of the post-stabilization behaviour. Ideally, steps (1) and (2) are carried on simultaneously, as in the case of computational form finding of cable nets. Lienhard⁸ compared governing variables in form-finding problems of form-active structures with variables in form finding of bending-active structures; he points out that in the former, the mechanical properties of the material are not governing the result, because shape is solely the result of equilibrium. As stated in this reference

… the form finding of bending-active structures is largely influenced by the length of a beam […] that is bent as a result of the constraining boundary conditions as well as the mechanical behaviour of the beam or shell elements.

Therefore, the complex nature of the mechanical problem has led in some cases to design processes in which the shape is initially searched using physical models and afterwards the initial state is simulated by means of finite element models, in a so-called ‘integral approach’ to the design process.⁸

Mainstream computational form-finding techniques rely on dynamic relaxation (DR) with different underlying mechanical models. From a numerical point of view, DR can be classified as an explicit method to find an equilibrium configuration. As suggested by D’Amico et al.,⁹ explicit methods are well suited to form-finding problems – for which tentative initial configurations may be far from equilibrium – in contrast with implicit solution methods, widely used in physical simulation of flexible structures under prescribed forces or displacements.

In addition to the pure form-finding problem, the need to solve the mechanical problem appears in the aforementioned steps (2) and (3). The second step may require tracking the full deformation path of the structural elements from an initial unstressed configuration. This usually involves traversing critical points in the equilibrium path, because the target configuration is a post-buckling state of the initial system. After the target shape has been reached and stabilized, the behaviour may not involve large displacements any more (unless the loads are again very large), but it is influenced by the inherited stress state. In this context, four basic features are required for a reliable simulation: (1) a sound mechanical model for flexible members undergoing large displacements of the reference lines and large rotations of the cross sections, and a numerical implementation capable of (2) traversing critical points in the equilibrium path, (3) allowing the addition and removal (in each stage) of certain members required to bring the system to the desired shape and/or stabilize it and (4) inheriting the stress state corresponding to a previous stage/configuration. Lienhard⁸ has analysed two commercial finite element software packages and shown some drawbacks regarding features (2), (3) and (4) that lead to disregarding one of them and adapting and completing the other to meet these needs. In this reference, Lienhard proposes a form-finding procedure based on finite element simulations that start from undeformed/unstressed members and drive them to a target configuration by shortening notional elastic cables. This process has been recently improved using an isogeometric finite element implementation (Bauer et al.¹⁰) that allows integration into a computer-aided design (CAD) environment (Bauer et al.¹¹ and Längst et al.¹²).

Regarding the first mentioned requirement, the mechanical model, numerous references can be found in the context of finite element techniques. However, in the applications to bending-active structures, most attention has been directed to developing the geometry-related aspects of form finding, and to a lesser extent to the mechanics-related aspects. The following specific issues need be taken into consideration in the mechanical modelling of structures with bending-active members:

Large displacements and large rotations of member cross sections with respect to an initial configuration (generally not in equilibrium in a form-finding process) are to be expected.

Because active structures need to behave elastically in the target (bent) configuration – or equivalently, they should be far from yielding – active members are designed using materials with a high strength-to-Young’s modulus ratio and low-depth cross sections to reach high flexibility and resilience.

The previous observation implies that, in spite of displacements and rotations being large, strains will remain small in the target configuration.

Therefore, beam theories considering linear elastic material will be adequate to reproduce the behaviour of bending-active structures. The objectives of this work are to provide an overview of the mechanical models that have been used in computational form finding of bending-active structures and to review and classify the most relevant references with respect to this point of view. Only references in which explicit methods as DR are used for form finding will be considered, although as it has been mentioned before, other methods have also been used.

Key aspects in the specialization of flexible beam models for DR

It is assumed that the reader is familiar with the DR method to determine the equilibrium configuration of a mechanical system (refer, for example, to Barnes¹³). The system to be analysed is discretized in a set of nodes. Fictitious masses are associated with model nodes. Stiffness relationships between nodes are defined. The analysis is started from an initial out-of-equilibrium configuration. Internal link forces corresponding to nodal positions at a given instant are calculated using stiffness relationships; external forces acting at nodes can be also introduced. With the resulting residual forces acting on every node, accelerations, velocities and displacements are evaluated. Artificial damping (or suppression of the kinetic energy) is used to bring the system to rest in static equilibrium.

The DR framework has been specialized to form finding of bending-active structures in different manners. The following aspects are tightly related and define the model specialization:

Type of discretization;

Underlying beam theory.

There are two types of discretizations in DR implementations reported in the literature: finite-difference-like discretizations and finite-element-like discretizations. In both cases, a key question is the way in which force–configuration relationships between nodes are established. Those relationships are determined in each case by the underlying beam theory, defined by the assumed kinematics, the selected elastic energy terms and the corresponding constitutive relations.

Beam theories for form finding of bending-active structures

Computational models for flexible rods rely on a suitable beam theory that models the three-dimensional (3D) deformation of the elastic body using a reduced set of kinematic and static variables. All examined references make use of the classical assumption of cross sections remaining plane and perpendicular to the deformed centreline of the rod, although with some differences between them. This assumption is a natural choice because strains induced by shear forces, which cause loss of planarity or perpendicularity of cross sections, are negligible in very slender members. The next sections review relevant aspects of the main beam theories in the context of form finding of bending-active structures. They are presented in order of increasing complexity; each theory may be considered a subset of the following one.

Euler–Bernoulli theory

The first reported solution for planar large deformations of a rod was developed by Euler in 1744, applying calculus of variations to the expression of the bending energy of the rod and assuming that (1) curvatures are proportional to bending moments – as suggested to him by Bernoulli – and (2) the rod length remains constant; this is equivalent to assuming inextensibility of the rod centreline and neglecting shear deformability. Conjugate variables in Euler–Bernoulli theory are centreline curvatures $κ$ (strains) and bending moments $M$ (section forces). The assumed constitutive relation is

M = E I κ

(1)

and the deformation energy is given by

U_{E - B} = \frac{1}{2} \int_{L} E I κ^{2} d S

(2)

where $d s = d S$ is the centreline length differential in the actual and the reference configuration, and they are equal due to the inextensibility assumption. Axial and shear forces $N and V$ are related to $M$ through equilibrium equations and are not associated with the deformation energy. Configuration variables can be reduced to $x$ and $y$ , positions of the centreline points referred to an external fixed frame. Then, $d s = \sqrt{1 + {(y^{'})}^{2}} d x$ , with the prime $(^{'})$ denoting derivatives with respect to $x$ . The exact expression of the curvature in terms of centreline positions is

κ = \frac{1}{R} = \frac{y^{″}}{{(1 + {(y^{'})}^{2})}^{3 / 2}}

(3)

The configuration of a flexible rod under compressive loads applied at both ends deduced from this theory is called elastica. The linearized version of the curvature, $κ = y^{″} = v^{″}$ , with $v$ denoting the beam transverse displacement, is used to formulate the linear planar bending theory of beams.

A refinement of Euler–Bernoulli theory can be obtained relaxing the inextensibility assumption. It involves the elongation energy

\frac{1}{2} \int_{L} E A ε^{2} d S

(4)

together with the constitutive relation between axial force and centreline strain

N = E A ε

(5)

and the definition of the axial strain

ε = \frac{d s}{d S} - 1

(6)

Non-linear solutions derived with this theory are called extensible elasticas. For the linear beam theory, equation (6) reduces to $ε = u^{'}$ .

Kirchhoff–Love theory

In 1859, Kirchhoff¹⁴ proposed a 3D theory for very slender rods with extensible centreline that was developed by Love.¹⁵ Kirchhoff–Love theory leads to negligible shear deformations, keeps Bernoulli’s moment-to-curvature proportionality in both principal axes and models torsion according to Saint–Venant’s theory (Dill¹⁶); it has been extensively used for developing finite elements for beams in the realm of small or moderately large displacements and rotations. Conjugate variables can be defined with respect to a set of orthogonal axes ${1, 2, 3}$ located at each point of the deformed centreline, with axis $1$ being tangent to the centreline, and $2,3$ coincident with the cross-sectional principal axes: centreline axial strains $ε$ , torsional strains $κ_{1}$ and curvatures $κ_{2}, κ_{3}$ are conjugate to axial forces $N$ , torsional moments $M_{1}$ and bending moments $M_{2}, M_{3}$ , respectively. Shear forces $V_{2}, V_{3}$ do not deform the rod and are related to previous variables through equilibrium relations. Constitutive equations are

N = E A ε M_{1} = G J κ_{1} M_{2} = E I_{2} κ_{2} M_{3} = E I_{3} κ_{3}

(7)

and the deformation energy is

U_{K - L} = \frac{1}{2} \int_{L} (E A ε^{2} + G J {κ_{1}}^{2} + E I_{2} {κ_{2}}^{2} + E I_{3} {κ_{3}}^{2}) d S

(8)

Configuration variables are positions of centreline points referred to a fixed frame, $x, y, z$ , and rotational parameters needed to describe cross-sectional rotations. Several choices can be made in this respect, and an adequate selection of sectional reference frames (using torsion-free Bishop¹⁷ frames) allows the use of only one rotational parameter describing the section torsional rotation $θ_{1}$ , because cross sections remain perpendicular to centreline tangents. Thus, $κ_{1} = d θ_{1} / d S$ , and the curvatures $κ_{2}$ and $κ_{3}$ can be expressed through differential-geometric relations in terms of centreline positions and their derivatives. The derivation of linear theories based on Kirchhoff’s kinematics is much simpler because the curvatures reduce to $κ_{2} = v^{″}$ and $κ_{3} = w^{″}$ , with $v$ and $w$ denoting transverse displacements in the beam. In finite element implementations, it is usually convenient to keep the (small) rotational degrees of freedom (DoFs) $θ_{2}$ and $θ_{3}$ and use the linearized constraints $θ_{2} = v^{'}$ and $θ_{3} = w'$ . Then, curvatures are simply the first derivatives of the flexural rotations.

If the rod is considered inextensible, axial forces become a consequence of equilibrium, and Kirchhoff–Love theory specializes into the 3D version of Euler–Bernoulli theory.

Other theories

Reissner¹⁸ proposed a 3D theory in which cross-sectional rotations are independent from centreline tangents. In this theory, kinematics is based on Timoshenko’s assumption for shear-deformable beams – originally introduced for small displacements and deformations – according to which cross sections remain plane but not necessarily perpendicular to the centreline. In Reissner’s theory, as in both previously described ones, displacements and rotations may be arbitrarily large. In addition to the previous constitutive equations (7), equations relating shear forces to shear strains $γ_{2} and γ_{3}$ have to be introduced

V_{2} = G A_{v 2} γ_{2} V_{3} = G A_{v 3} γ_{3}

(9)

and additional terms for shear deformation energies need also be added to equation (8)

\frac{1}{2} \int_{L} (G A_{v 2} {γ_{2}}^{2} + G A_{v 3} {γ_{3}}^{2}) d S

(10)

In this case, centreline positions and sectional (finite) rotations are used as configuration variables. This theory has been successfully developed by Simo¹⁹ and other authors for the simulation of flexible beams in the context of finite elements, using finite rotations to model the behaviour of cross sections. The adjective geometrically exact has been used to refer to these models.

To put these theoretical developments into perspective, it is necessary to mention the more general Cosserat rod theory that models the rod as a space curve with three vectors attached to each point. The first is the tangent vector to the centreline curve; the other two are called directors and characterize the configuration of material fibres of the cross section. Cosserat rods were developed by Green et al.²⁰ on the basis of Cosserats’²¹ theory for directed continua.

Both Reissner–Simo and Kirchhoff–Love theories can be classified as special cases of Cosserat rods: in the first case, directors are chosen orthonormal in the undeformed state and are constrained to remain orthonormal during the deformation; in the second (Kirchhoff–Love), orthogonality to the tangent vector is also enforced. More complex theories considering out-of-plane deformation modes of cross sections have been developed by other authors (as Simo and Vu-Quoc²² or Hodges²³) although their increased number of variables makes them more difficult to apply to shape finding of bending-active structures.

Types of discretization in DR procedures for flexible rods

Two groups of discretizations can be found in the literature: (1) finite-difference-like discretizations using 3 or 4 DoFs per node and (2) finite-element-like methods with 6 DoFs per node. The main distinction is found in the way strain measures and internal forces exerted to model nodes are evaluated.

In finite-difference discretizations, discrete strain measures at nodes are deduced from difference schemes between nodal degrees of freedom using groups of two, three or four nodes. Internal forces acting on nodes may be calculated either (1) using the discrete strain measures (assuming linear constitutive behaviour) and ad hoc equilibrium relations or (2) deriving expressions of the internal forces from energy methods and discretizing these expressions by means of finite-difference schemes.

In the case of finite-element-like discretizations, the link between nodes is considered as a beam finite element. Relations between element end forces and relative displacements/rotations of end-nodes are derived making use of beam theory and energy methods. In the following sections, we review the main contributions in each group.

Finite-difference discretizations

In this group, we are placing in chronological order the contributions by Adriaenssens and Barnes,²⁴ Barneset al.,²⁵ Du Peloux et al.²⁶ and D’Amico et al.⁹

Adriaenssens and Barnes²⁴ is based on the Euler–Bernoulli theory with the extensibility assumption. In a certain configuration, axial strains $ε$ are calculated from nodal coordinate differences, and link forces in the direction of lines joining nodes are obtained using the constitutive equation (5) for axial forces. In order to account for bending effects, groups of three nodes are considered: the bending moment at node $j$ is obtained using a discrete measure of the curvature of each three-node group instead of equation (3)

κ_{j} = 2 \frac{\sin α_{j}}{d_{i k}}

(11)

and the constitutive equation (1) $M_{j} = E I κ_{j}$ (see Figure 1). The contribution of the three-node group to in-plane shear forces acting on each node is deduced from equilibrium, so that end forces produce a moment $M_{j}$ on the middle $j$ node

S_{i j} = \frac{M_{j}}{d_{i j}} S_{j k} = \frac{M_{j}}{d_{j k}}

(12)

Figure 1.

Discretization in Adriaenssens and Barnes.²⁴

Adriaenssens and Barnes prove that, for isotropic cross sections (sections with uniform second moment of area), torsional moments are related to bending moments through an equilibrium equation and no additional forces due to torsion are needed for the DR procedure. Thus, shear forces are vector-added to link (axial) forces in order to compute residual forces acting on each node, and the DR procedure can be performed. In this way, only 3 DoFs per node are needed, and the simplest beam theory is used in the process, with the counterpart of the limitation imposed to cross sections, which need be square or circular shaped (full or hollow).

Barnes et al.²⁵ account for the case of non-uniform cross sections: in addition to considering the effect of bending through three-node in-plane shear forces as in the previous reference, out-of-plane shear forces acting on four-node groups are calculated using the constitutive equation for torsion (equation (7)) $M_{j k 1} = G J κ_{j k 1}$ and a discrete measure of the torsional strain $κ_{1}$ obtained from the nodal coordinates of the four-node group (see Figure 2)

κ_{j k 1} = \frac{ϕ - ϕ_{0}}{d_{j k}}

(13)

Figure 2.

Discretization in Barnes et al.²⁵

where $ϕ_{0}$ is the angle formed by planes $i j k$ and $j k l$ at the start of the process. Out-of-plane forces $P_{i j}$ and $P_{k l}$ are determined to be in equilibrium with $M_{j k 1}$ . With this setup, in-plane as well as out-of-plane forces are added to axial forces to calculate residuals for DR. The underlying theory is a simplified version of Kirchhoff–Love theory, in which (apart from the constitutive equations for extension and torsion) only one bending constitutive equation is entering the formulation (in-plane bending). This simplifying assumption has two consequences: on one side, the orientation of one principal bending axis at a node is forcefully contained into the plane formed by the node itself and the two adjacent ones; on the other, out-of-plane bending stiffness has to be indirectly considered through a factor that modifies the torsional stiffness $G J$ . This discretization keeps the advantage of working with only 3 DoFs in the DR process.

The method proposed by Du Peloux et al.²⁶ is close to recent developments in the field of computer animation (Bertails et al.²⁷ and Bergou et al.²⁸) which relate the geometry of discrete curves (an ordered set of nodes) to their mechanical behaviour using a finite-difference scheme (Figure 3). In this reference, the minimal set of degrees of freedom per node that keeps full consistency with Kirchhoff–Love theory is used: 3 DoFs for nodal positions and one rotational DoF to keep track of torsion. The latter measures the angular difference between cross-sectional material (principal) frames and torsion-free (Bishop) frames. A discrete measure of oriented curvatures using the coordinates of three consecutive nodes is employed. Torsional strains are calculated as finite differences between nodal torsional rotations.

Figure 3.

Discretization in Du Peloux et al.²⁶

In order to calculate nodal forces and twisting moments, Du Peloux et al. discretize the expressions of the internal forces and torsional moment in the rod, obtained as derivatives of the elastic deformation energy with respect to positions and torsional rotation. A simplifying assumption for the evaluation of the derivatives is that torsional waves propagate in an instantaneous manner compared to bending waves. The discretization of internal forces and twisting moments to obtain forces and moments acting on nodes follows a finite-difference scheme. This computation depends on the torsional stiffness $G J$ , and the bending stiffness values $E I_{2}$ and $E I_{3}$ . The rod is considered, therefore, inextensible. Force and moment residuals are calculated at each node for performing DR. Inextensibility is enforced by means of penalty forces at nodes.

D’Amico et al.⁹ make also use of a discrete approach to 3D curved rods as Du Peloux et al. It can be considered an improvement of Barnes et al. because it keeps track of the angular difference between material frames and Bishop frames in case that torsional constraints are imposed at the ends of the rod. The orientation of tangent vectors at nodes is found using Catmull–Rom interpolation; the other two vectors of the Bishop frame are found by parallel transport and then rotated to obtain the material frame. However, in contrast to Du Peloux et al., these rotation angles are not considered DoFs of the system, but initially given data. Curvatures and torsional strains are calculated operating with coordinates and frame orientations at a given node and the immediately adjacent ones (three-node groups). The total moment at a node is the result of vector-adding the components given by constitutive equation (7). Nodal shears acting on every node of a group of three are calculated imposing equilibrium with the total moment at the mid node; they are vector-added to forces from subsequent groups and to the axial link forces at nodes in a similar fashion as in Barnes et al. to calculate residuals for DR. This setup makes, therefore, use of Kirchhoff–Love theory but forces the orientation of the principal axes at a node to keep a fixed angle to the corresponding three-node plane due to the restricted number of DoFs.

Table 1 includes a summary of the main features of finite-difference-like parametrizations. Moreover, all references use equilibrium relations to obtain nodal forces for the DR process, excepting the one by Du Peloux et al. that evaluates internal forces as derivatives of the deformation energy with respect to kinematic variables.

Table 1.

Comparison between references using finite-difference-like discretizations.

Reference	No. of DoFs	Beam theory	Stiffness parameters
Adriaenssens et al.²⁴	3	Extensible E-B	$E A$ and $E I$
Barnes et al.²⁵	3	Modified extensible K-L	$E A$ , $E I$ and $G J$
Du Peloux et al.²⁶	4	Inextensible K-L	$G J$ , $E I_{2}$ and $E I_{3}$
D’Amico et al.⁹	3	Extensible K-L	$E A$ , $G J$ , $E I_{2}$ and $E I_{3}$

DoF: degrees of freedom; E-B: Euler–Bernoulli; K-L: Kirchhoff–Love.

Finite-element discretizations

The following references have been included in this category: Li and Knippers,²⁹ D’Amico et al.,³⁰ D’Amico et al.³¹ and Senatore and Piker.³² They base on previous work by Williams thoroughly explained in Adriaenssens.³³ All of them share a common implementation feature that allows their classification into a well-established formulation for flexible rod finite elements: the co-rotational formulation. The rod is discretized into elements that can undergo large rotations, but at local level, use small-displacement/rotation relationships. We have also included Bessini et al.³⁴ that implements Reissner–Simo model, which is capable to deal with large rotations of cross sections with no limitation to their relative magnitude at the element level.

The proper numerical treatment of large rotations is a main concern in the computational solution of this type of problems. This fact was readily recognized by Argyris³⁵ to which the reader is referred for a detailed account of the following ideas. From a mathematical point of view, a rotation is an element of the special orthogonal group SO(3). The elements of this non-additive and non-commutative group can be numerically represented in several manners. Before reviewing the co-rotational formulation, the following section summarizes (in a non-exhaustive manner) the most important alternative parametrizations for rotations.

Parametrization of finite rotations

Rotations, as members of the special orthogonal group, can be represented by means of 3 × 3 orthogonal matrices $Λ$ of determinant 1. The columns of the rotation matrix are the components of the rotated frame. The main drawback of this representation is the high number of required parameters: nine non-independent scalars related by six conditions.

Euler’s theorem establishes the equivalence between the matrix representation of rotations and the pseudo-vectorial representation $θ$ : the modulus $θ$ of $θ$ represents the rotation angle and its orientation determines the rotation axis. The relationship between $Λ$ and $θ$ is given by Rodrigues formula³⁶

Λ = \cos θ 1 + \frac{\sin θ}{θ} \hat{θ} + \frac{1 - \cos θ}{θ^{2}} θ \otimes θ

(14)

where $\hat{θ}$ is the skew-symmetric matrix with axial vector $θ$ . This representation uses a minimum number of parameters and the correspondence with $Λ$ is one-to-one up to $θ = π$ , but it has the drawback that addition of pseudo-vectors does not represent composition of rotations (hence the pseudo- prefix). However, it has been successfully employed to formulate finite elements for large rotation problems.

An alternative pseudo-vectorial three-parametrization consists in using $ψ = 2 \tan (θ / 2) e$ , with $θ$ equal to the rotation angle and $e$ being a unit vector in the direction of the rotation axis. The relation between both pseudo-vectors is $ψ = (2 \tan (θ / 2) / θ) θ$ . The components of $ψ$ are sometimes called Rodrigues parameters (Crisfield³⁷). Rodrigues formula transforms through substitution into an expression with no trigonometric functions

Λ = 1 + \frac{1}{1 + \frac{1}{4} ψ^{T} ψ} (\hat{ψ} + \frac{1}{2} {\hat{ψ}}^{2})

(15)

This alternative pseudo-vector has been partially used by Williams (Adriaenssens³³) and by D’Amico et al.^30,31 (we will expand on this in the next section).

Other classical three-parameter representations, as Euler’s angles or Cardano’s angles are not advantageous compared to pseudo-vector representations and therefore haven’t found use in mechanical modelling of slender beams.

An alternative system to represent rotations is the use of unit quaternions. They are four-tuples of real numbers $(q_{0}, q_{1}, q_{2}, and q_{3})$ related by the equation $q_{0}^{2} + q_{1}^{2} + q_{2}^{2} + q_{3}^{2} = 1$ . The three later quantities can be interpreted as a vector $q$ . The relationship to the rotation matrix is

Λ = (q_{0}^{2} - 1) 1 + 2 q_{0} \hat{q} + 2 q \otimes q

(16)

Unit quaternions and rotation pseudo-vectors are related as follows

(q_{0}, q) = (\cos \frac{θ}{2}, \sin \frac{θ}{2} n) with n = \frac{θ}{θ}

(17)

Quaternions are endowed with an algebra that provides an efficient tool to operate with rotations; this fact and the limited number of parameters compared to rotation matrices make unit quaternions a preferred choice to store and keep track of rotations.

Argyris³⁵ showed that the exponential operator acting on a skew-symmetric matrix, $\exp (\hat{a}) = \sum_{1}^{\infty} (1 / n!) {\hat{a}}^{n}$ , results in an orthogonal matrix that can be calculated using Rodrigues formula (14). Considering the group character of rotations, this observation leads to the introduction of a local parametrization through spins $Δ θ$ , also called local rotational parameters: a new rotation $Λ^{(t + 1)}$ close to the original one $Λ^{(t)}$ is got through the formula

Λ^{(t + 1)} = \exp (\hat{Δ θ}) Λ^{(t)}

(18)

Spins have been used as rotational DoFs for finite element implementations based on Reissner’s beam theory (Simo¹⁹). They have been also used to update nodal frames in the finite-element-based discretizations that will be reviewed in the next section.

Discretizations based on the co-rotational formulation

The co-rotational technique was initially proposed in the 1970s by Wempner and by Belytschko et al. (Crisfield,³⁸ which includes a review of the history and the essentials of this method.)

The discretization associated with the co-rotational formulation splits the rod into two-node elements; at a given instant, the configuration is defined by nodal positions and nodal frames. The large-displacement/rotation mechanical problem is divided into two sub-problems: (1) modelling the mechanical behaviour of finite elements between nodes in terms of local displacements and rotations, and relating the latter to the changes in positions and orientations of each node; (2) keeping track of changes in nodal positions and nodal frames, with no restriction in their magnitude.

The first sub-problem requires to define and keep track of element frames, in order to quantify angular differences with nodal frames (Figure 4). If the discretization is sufficiently refined, these differences will be small, and the element behaviour can be even modelled with a linear beam theory, or considering moderately large displacements and small rotations.

Figure 4.

Co-rotational setup.

Nodal unit base vectors $t_{i α}$ , $t_{j α}$ and element base vectors $e_{α}$ allow to approximate the (small) nodal rotations at some instant as (Crisfield³⁷)

φ_{i 1} = \frac{1}{2} (t_{i 3}^{T} e_{2} - t_{i 2}^{T} e_{3})

(19a)

φ_{i 2} = \frac{1}{2} (t_{i 1}^{T} e_{3} - t_{i 3}^{T} e_{1})

(19b)

φ_{i 3} = \frac{1}{2} (t_{i 2}^{T} e_{1} - t_{i 1}^{T} e_{2})

(19c)

φ_{j 1} = \frac{1}{2} (t_{j 3}^{T} e_{2} - t_{j 2}^{T} e_{3})

(19d)

φ_{j 2} = \frac{1}{2} (t_{j 1}^{T} e_{3} - t_{j 3}^{T} e_{1})

(19e)

φ_{j 3} = \frac{1}{2} (t_{j 2}^{T} e_{1} - t_{j 1}^{T} e_{2})

(19f)

For explicit problems, once nodal rotations are known, linear beam relations between forces, moments, displacements and rotations in a beam can be used to get the applied moments and forces at the ends of each element. Alternatively, axial force–dependent non-linear terms could also be used to obtain end forces and moments. Force and moment residuals $f_{i}$ , $m_{i}$ acting on node $i$ are subsequently evaluated, and from them, linear and angular accelerations $a_{i}$ , $α_{i}$ follow from D’Alembert’s principle

f_{i} = m_{i} a_{i}

(20a)

m_{i} = I_{i} α_{i} + ω_{i} \times I_{i} ω_{i}

(20b)

Because there is no need to exactly reproduce the dynamics in the DR process, an isotropic inertia tensor may be chosen in order to simplify equation (20b)

m_{i} = I_{i} α_{i}

(21)

With the selected time-integration technique, nodal displacements ∆x _i and local rotation increments (spins) $Δ θ_{i}$ corresponding to the time step $t \to t + 1$ are calculated.

The second sub-problem involves updating large rotations of nodal and element frames from a given configuration to the next. In the case of nodal frames, the multiplicative update from equation (18) together with Rodrigues formula (14) to evaluate the exponential are used

t_{i α}^{(t + 1)} = \exp ({\hat{Δ θ}}_{i}) t_{i α}^{(t)}

(22)

for node $i$ and $α = 1, 2, 3$ . For element frames, the unit vector $e_{1}^{(t + 1)}$ is computed in a straightforward way from the updated element nodal positions. However, updating the other two vectors requires an interpolation between previously evaluated element nodal frames. Two approaches to this problem are described in Crisfield.³⁷

All reviewed references in this group use the co-rotational framework, with some differences in the rotation updating schemes and the mechanical model for the beam element.

A simplification, common to all references, is to avoid the introduction of full element frames: only $e_{1}$ is computed (as the normalized difference of updated nodal positions). Then, equations (19) reduce to (Figure 5)

Δ φ_{i j 1} = \frac{1}{2} (t_{i 2}^{T} t_{j 3} - t_{i 3}^{T} t_{i 2})

(23a)

φ_{i 2} = - t_{i 3}^{T} e_{1}

(23b)

φ_{i 3} = t_{i 2}^{T} e_{1}

(23c)

φ_{j 2} = - t_{j 3}^{T} e_{1}

(23d)

φ_{j 3} = t_{j 2}^{T} e_{1}

(23e)

Figure 5.

Simplified co-rotational setup.

Considering that equations (19a)–(19f) were already an approximation, the simplified version requires a fine discretization in order to neglect terms that involve element frame vectors $e_{2}$ and $e_{3}$ .

There are some differences in the way that spins $Δ θ$ are calculated: Li and Knippers,²⁹ as well as Senatore and Piker,³² follow a direct process at each node $i$ : (1) calculate the angular acceleration from the non-balanced moment at the node $α_{i} = I^{- 1} m_{i}$ , (2) calculate the angular velocity using the time step suited to the selected time-integration scheme $ω_{i}^{(t + 1)} = ω_{i}^{(t)} + Δ t α_{i}$ and (3) approximate the spin for the step $t \to t + 1$ by $Δ θ_{i} = Δ t ω_{i}^{(t + 1)}$ . However, D’Amico et al.,^30,31 following Williams (Adriaenssens),³³ include an additional procedure to transform the incremental rotations $Δ t ω_{i}^{(t + 1)}$ into the spin $Δ θ_{i}$ . For that purpose, they calculate the rotation axis through successive rotations of the nodal local axes using the alternative rotational pseudo-vectors $ψ_{α} = 2 \tan (Δ θ_{i α} / 2) t_{i α}$ presented in the previous section and calculated with the components $Δ θ_{i α}$ of $Δ t ω_{i}^{(t + 1)}$ ; in this process, they miss the advantage of using the alternative version of Rodrigues formula (15).

With regard to the mechanical modelling of the beam element, Li and Knippers use linear stiffness relations corresponding to a (linear) Kirchhoff beam between nodal local displacements and rotations and local end forces.

However, Adriaenssens³³ describes the way followed by Williams to compute element end forces and moments. First, non-linear strain measures (elongation $ε$ , twist $κ_{1}$ and curvatures $κ_{2}$ , $κ_{3}$ ) are computed in terms of nodal distance (initial and actual) and nodal local rotations (equation (23)); second, Kirchhoff’s strain energy (equation (8)) is derived with respect to the kinematic variables in order to get axial force, torsional moment and end bending moments at nodes $i$ and $j$

N = \frac{\partial U}{\partial ε}

(24a)

M_{i 1} = \frac{\partial U}{\partial κ_{1}} M_{j 1} = - M_{i 1}

(24b)

M_{i β} = \frac{\partial U}{\partial φ_{i β}} M_{j β} = \frac{\partial U}{\partial φ_{j β}} β = 2, 3

(24c)

D’Amico et al.^30,31 adopt the same setup. The first reference has a summary of the resulting expressions. It is remarkable that the axial force has non-linear terms as well as rotation dependent terms due to the structure of the adopted elongation strain, and the end bending moments have geometric terms depending on the axial force, added to the standard linear terms as used by Li and Knippers.

In order to transform the scalar values into oriented force and moment terms, all references in this group need to make up for the lack of element frame vectors $e_{2}$ and $e_{3}$ . This is done in two steps: (1) assuming that the components of moment vectors at each end of the element are equal to the moments got in equations (24), despite the orientation of each nodal frame is different

M_{i} = \sum_{α = 1}^{3} M_{i α} t_{i} M_{j} = \sum_{α = 1}^{3} M_{j α} t_{j}

(25)

and (2) calculating end shear forces such that the element is in equilibrium of moments, and vector-adding them to the oriented axial forces to get force vectors acting at element ends.

The reviewed parametrizations based on the co-rotational formulation have very similar features: all of them use 6 DoFs, rely on Kirchhoff–Love beam theory and therefore use the corresponding stiffness parameters $E A$ , $G J$ , $E I_{2}$ and $E I_{3}$ . Only Li and Knippers use linear stiffness relations between kinematic variables and element end forces/moments. There is also a difference in the way the references by D’Amico et al. and the others perform updating of nodal frames.

Discretization based on the geometrically exact model

Simo¹⁹ and Simo and Vu-Quoc³⁹ carried out the first finite element implementation of Reissner’s model. The implementation is termed geometrically exact in the literature, because of the mathematically exact treatment of rotations. The element kinematics considers that centreline positions $x$ and cross-sectional rotations $Λ$ vary along the element. For an initially straight element with initial tangent vector $t_{1}^{(0)}$ , strains are expressed by means of two vectors

γ = x^{'} - t_{1}^{(0)}

(26a)

κ = Λ^{'} Λ^{T}

(26b)

The first vector gathers the elongation and the shear strain components; the second one is a measure of the rate of rotation of cross sections along the centreline. Note that, in general, cross-sectional normals won’t be coincident with centreline tangents due to shear deformation. Linear elastic constitutive relations are implemented in a straightforward manner considering six stiffness parameters: $E A, G A_{v 2}, G A_{v 3}$ and $G J, E I_{2}, E I_{3}$ which relate axial and shear forces, and twisting and bending moments to the strain components expressed in the intrinsic reference frame.

The simplest implementation discretizes the structure into two-node elements, selects displacements ∆x and spins $Δ θ$ as discrete variables (6 DoFs) and uses linear Lagrangian interpolating functions and reduced integration to evaluate integrals. The translation of the geometrically exact model to the DR algorithm has been developed and is described in Bessini et al.³⁴ At a certain step in the DR process, forces and moments acting at element ends result from the discretized internal virtual work. With them, force and moment residuals at nodes are evaluated and linear and angular accelerations yield from D’Alembert’s principle – equations (20a) and (21). A multiplicative update is used to update rotations – equation (18) – and rotations are stored using quaternions.

Discussion

The main argument for developing 3-DoF finite-difference-like discretizations for form finding of bending-active structures is that coupling of rotational DoFs with translational DoFs and axial stiffness is the cause of ill-conditioning of the DR process.²⁵ Adriaenssens³³ compared the results of the 3-DoF model with those obtained using the 6-DoF discretization proposed by Williams for the simulation of the planar elastica problem with different load magnitudes and found that while both methods provided similar accuracy in the geometry of the solution, the 3-DoF model was remarkably faster. This was an expected conclusion because the latter model is substantially simpler. It is interesting to follow how the 3-DoF model has been gradually refined, from the one first developed by Adriaenssens,^24,33 only valid for isotropic cross sections, through the model by Barnes et al.²⁵ that considers non-isotropic cross sections but models out-of-plane bending stiffness through the introduction of an artificial torsion constant affecting the torsional stiffness, to the model by D’Amico et al.⁹ that calculates Bishop frames to refer the orientation of the bending axes and therefore uses all strain measures and stiffness constants of the Kirchhoff–Love beam at a higher computational cost than the other two models. The fact is that using only translational DoFs imposes an internal constraint to the orientation of the section axes: in the model by Barnes, the reference bending axis at a node is forced to stay in the plane formed by the node and the two adjacent ones; in the model by D’Amico, it is forced to keep a pre-determined orientation with respect to the Bishop frame defined by the discrete curve in each step.

The 4-DoF discretization by Du Peloux et al. offers a trade-off between computational cost and model accuracy, as it uses the minimal set of parameters that is fully consistent with Kirchhoff–Love theory. Although the theoretical foundation presented by the authors considers the inextensible Kirchhoff model, the computational implementation imposes the inextensibility constraint through node-to-node directed penalty forces that can be adjusted to be actual axial forces by using the real cross-sectional axial stiffness.

All 6-DoF finite-element-like discretizations share the disadvantage of requiring the highest number of parameters. The comparative computational cost with respect to 4-DoF discretizations remains to be analysed, because in the latter, the cost of keeping track of Bishop frames and material frames during the process might be comparable to the cost of updating nodal frames in the former. In spite of the conditioning problems of 6-DoF discretizations mentioned by Adriaenssens and Barnes, all reviewed references report an adequate performance of their models. A further advantage of these discretizations is the ability to directly integrate the result of form finding into a finite element tool to further analyse the structural behaviour under additional loads, or to trace the construction sequence back to an undeformed configuration. Integration with an analysis tool may be much more difficult in the case of finite-difference discretizations (with perhaps the exception of the 4-DoF implementation). A particular feature of the implementation by Bessini et al. is the ability to deal with shear-deformable structures. In the realm of active bending, this can be applied to the design of structures with members consisting of two or more layers of laths when considered as structural elements with large shear deformations, such as the Mannheim Multihalle.

Final remarks

Computational form finding of bending-active structures is a challenging problem. DR has been used in many references as the basis for form finding; however, under this common umbrella, approaches show significant differences in both the underlying beam theory and the type of discretization.

Discretizations have been classified into two groups: finite-difference-like and finite-element-like. Both rely on a distribution of nodes along the centreline of a flexible beam and consider masses (and inertias) to be lumped at nodes. The first group evaluates strains and internal forces acting at nodes using the relative positions of groups of nodes and the stiffness parameters of cross sections, whereas the second group treats the links between elements as finite elements and evaluates strains and internal forces from the configuration of adjacent nodes.

Finite-difference discretizations use 3 or 4 DoFs per node, and range from simple implementations based on planar Euler–Bernoulli theory, suited for isotropic cross sections, to more complex models based on 3D Kirchhoff–Love theory that account for torsion and spatial bending in an accurate way. References report advantages in terms of computational cost and robustness of the DR algorithm.

Finite-element discretizations employ 6 DoFs per node and are built on Kirchhoff–Love or Reissner–Simo beam theory. Reviewed references implement simplified versions of co-rotational elements or geometrically exact elements. Their main advantage is that results from the form finding process can be directly exported to finite element models.

The reviewed references show that there is a range of available models/discretizations to be applied in DR algorithms for very flexible beams. However, there are still unexplored questions that are worth studying, as carrying out performance comparisons between models, or even implementing alternative discretizations for flexible beams that have been tested in the finite element realm.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The financial support from the Spanish Ministry of Economy and Competitiveness through grant BIA2015-69330-P (MINECO) is gratefully acknowledged.

ORCID iD

Carlos Lázaro

References

Happold

Liddell

Timber lattice roof for the Mannheim Bundesgartenschau. Struct Eng 1975; 53(3): 99–135.

Douthe

Baverel

Caron

JF.

Gridshell in composite materials: towards wide span shelters. J Int Ass Shell Spatial Struct 2007; 48(3): 175–180.

Nicholas

Lafuente-Hern á ndez

Gengnagel

. The Faraday pavilion: activating bending in the design and analysis of an elastic gridshell. In: Proceedings of the symposium on simulation for architecture & urban design (SimAUD 2013), San Diego, CA, 7–10 April 2013, pp. 154–161. San Diego, CA: Society for Modeling & Simulation International.

Pone

Colabella

Parenti

et al . Construction and form-finding of a post-formed timber grid-shell. In: Proceedings of the 2nd international conference on structures and architecture (ICSA 2013), Guimaraes, 24–26 July 2013, pp. 245–252. Boca Raton, FL: CRC Press.

Harding

Pearson

Lewis

et al . The ongreening Pavilion. In: Adriaenssens

Gramazio

Kohler

(eds) Advances in architectural geometry. New York: Springer International Publishing, 2016, pp. 295–308.

Gengnagel

Mobile Membrankonstruktionen. PhD Thesis, Technische Universität München, München, 2005.

Lienhard

Alpermann

Gengnagel

et al . Active bending, a review on structures where bending is used as a self-formation Process. Int J Space Struct 2013; 28(3): 187–196.

Lienhard

Bending-active Structures: form-finding Strategies using elastic deformation in static and kinetic systems and the structural potentials therein. Number 36 in itke Forschungsberichte, Universität Stuttgart – Institut für Tragkonstruktionen und Konstruktives Entwerfen, 2014, https://elib.uni-stuttgart.de/handle/11682/124

D’Amico

Zhang

Kermani

A finite-difference formulation of elastic rod for the design of actively bent structures. Eng Struct 2016; 117: 518–527.

10.

Bauer

Breitenberger

Philipp

et al . Nonlinear isogeometric spatial Bernoulli beam. Comput Method Appl Mech Eng 2016; 303: 101–127.

11.

Bauer

Längst

Wüchner

et al . Isogeometric analysis for modelling and simulation of building processes. In: Proceedings of the IASS annual symposium 2017–Interfaces: architecture engineering science, Hamburg, 25 September 2017.

12.

Längst

Bauer

Michalski

et al . The potentials of isogeometric analysis methods in integrated design processes. In: Proceedings of the IASS annual symposium 2017–Interfaces: architecture engineering science, Hamburg, 25 September 2017.

13.

Barnes

Form finding and analysis of tension structures by dynamic relaxation. Int J Space Struct 1999; 14: 89–104.

14.

Kirchhoff

Ueber das Gleichgewicht und die Bewegung eines unendlich dünnen elastischen Stabes. J Pure Appl Math 1859; 56: 285–313.

15.

Love

AEH

. A treatise on the mathematical theory of elasticity. Cambridge: Cambridge University Press, 1927.

16.

Dill

EH.

Kirchhoff’s theory of rods. Arch Hist Exact Sci 1992; 44(1): 1–23.

17.

Bishop

RL.

There is more than one way to frame a curve. Am Math Mon 1975; 82(3): 246–251.

18.

Reissner

On one-dimensional large-displacement finite-strain beam theory. Stud Appl Math 1973; 52(2): 87–95.

19.

Simo

JC.

A finite strain beam formulation: the three-dimensional dynamic problem – part I. Comput Method Appl Mech Eng 1985; 49: 55–70.

20.

Green

Naghdi

Wenner

ML.

On the theory of rods – part II: developments by direct approach. P R Soc Lond A 1974; 337: 485–507.

21.

Cosserat

Théorie Des Corps Déformables. Paris: Librairie Scientifique A. Hermann et Fils, 1909.

22.

Simo

Vu-Quoc

A geometrically exact rod model incorporating shear and torsion-warping deformation. Int J Solids Struct 1991; 27(3): 371–393.

23.

Hodges

DH.

Nonlinear composite beam theory. Reston, VA: American Institute of Aeronautics and Astronautics, 2006.

24.

Adriaenssens

Barnes

Tensegrity spline beam and grid shell structures. Eng Struct 2001; 23: 29–36.

25.

Barnes

Adriaenssens

Krupka

A novel torsion/bending element for dynamic relaxation modeling. Eng Struct 2013; 119: 60–67.

26.

Du Peloux

Tayeb

Lefevre

et al . Formulation of a 4-dof torsion/bending element for the formfinding of elastic gridshells. In: Proceedings of the IASS symposium, Amsterdam, 17–20 August 2015.

27.

Bertails

Audoly

Cani

et al . Super-helices for predicting the dynamics of natural hair. In: Proceedings of ACM SIGGRAPH, vol. 25, Boston, MA, 30 July–3 August 2006, pp. 1180–1187. New York: ACM Transactions on Graphics.

28.

Bergou

Wardetzky

Robinson

et al . Discrete elastic rods. In: Proceedings of ACM SIGGRAPH, vol. 27, Los Angeles, CA, 11–15 August 2008, Article no. 63. New York: ACM Transactions on Graphics.

29.

Knippers

. Rotation formulations for dynamic relaxation – with application in 3D framed structures with large displacements and rotations. In: Proceedings of the IASS-APCS Symposium, Seoul, 21–24 May 2012.

30.

D’Amico

Kermani

Zhang

Form finding and structural analysis of actively bent timber grid shells. Eng Struct 2014; 81: 195–207.

31.

D’Amico

Kermani

Zhang

et al . Optimization of cross-section of actively bent grid shells with strength and geometric compatibility constraints. Comput Struct 2015; 154: 163–176.

32.

Senatore

Piker

Interactive real-time physics: an intuitive approach to form-finding and structural analysis for design and education. Comput Aided Design 2015; 61: 32–41.

33.

Adriaenssens

Stressed spline structures. PhD Thesis, University of Bath, Bath, 2000.

34.

Bessini

Lázaro

Monleón

A form-finding method based on the geometrically exact rod model for bending-active structures. Eng Struct 2017; 152: 549–558.

35.

Argyris

An excursion into large rotations. Comput Method Appl Mech Eng 1982; 32(1–3): 85–155.

36.

Rodrigues

Des lois géométriques qui régissent les déplacements d’un système solide dans l’espace et de la variation des coordonnées provenant de ces déplacements considérés ind’ependamment des causes qui peuvent les produire. J Math Pures Appl 1840; 5: 380–440.

37.

Crisfield

MA.

Non-linear finite element analysis of solids and structures, vol. 2. Hoboken, NJ: John Wiley & Sons, 1197.

38.

Crisfield

MA.

Non-linear finite element analysis of solids and structures, vol. 1. Hoboken, NJ: John Wiley & Sons, 1991.

39.

Simo

Vu-Quoc

A three-dimensional finite-strain rod model – part II: computational aspects. Comput Method Appl Mech Eng 1986; 58: 79–116.