A 4-degree-of-freedom Kirchhoff beam model for the modeling of bending–torsion couplings in active-bending structures

Abstract

Gridshells are lightweight structures made of interconnected slender beams. Due to large displacements, high interaction between the beams, and bending–torsion coupling, modeling gridshells requires specific non-linear numerical tools to reach convergence within a reasonable time. In this article, the development of such a tool is presented. It is based on the Kirchhoff beam theory and uses the dynamic relaxation method. First, from Kirchhoff’s equations, the internal forces and moments acting on a beam are obtained. Once this mathematical work is done, the dynamic relaxation method is used in order to get the static equilibrium configuration of the beam. This new approach is tested on several examples and validated for slender beams with arbitrary rest-state configuration and cross sections. In particular, results for ribbons with high bending–torsion coupling are presented. Finally, this process enables the fast and precise modeling of gridshells including bending–torsion coupling.

Keywords

active bending composite materials dynamic relaxation gridshell Kirchhoff beam torsion

Introduction

Context: structures in active bending

When designing a structure, one has to keep in mind the process to erect it. And when it comes to building structures made of curved elements, two solutions arise: (1) the first is to produce curved elements and to assemble them on site and (2) the second is to produce straight elements and to bend them on site. The second solution appears to be the most cost-effective, as it avoids the complex production and transportation of curved elements.

In return, the resulting structures are in so-called active bending,³ giving rise to more complex behavior, as we are working with elements in finite displacements and in long-lasting elastic deformations. One example of such active-bending structures is the elastic gridshell shown in Figure 1.⁴ Elastic gridshells are indeed the result of the bending and the assembly of initially straight slender beams.

Figure 1.

Solidays’ gridshell built in June 2011.

The need for numerical tools

The simplicity and the cost-effectiveness of the construction of such structures have a price at the design step: the modeling is made more difficult by the finite displacements and the induced prestress. There is therefore an ongoing research for adequate fast and easy-to-use modeling tools.^5
–8

This article gives deep insights on the modeling of slender beams accounting for both torsion and bending behaviors in the context of largely bent structures.

They are indeed the structural elements used for gridshells, which is the starting point of our work, but they are also extensively used in the field of active bending, for their lightweight and great bending ability in elastic deformations.

The challenge to be met was to propose a tool suited to the modeling of elastic slender beams in large displacement. Such a tool should not only be fast and precise; as we are dealing with simulations in large displacement and with multiple boundary conditions, the calculation should also be easy to prepare and minimize the precautions to be taken to ensure convergence. Our purpose has then been to go beyond in the research for numerical tools, by meeting these requirements.

In the field of active bending, the dynamic relaxation method is often the preferred one for the modeling of slender beams.^5
–10 This is due to its simplicity of implementation and its good convergence, compared with other methods. For example, it has been found in Bouhaya¹¹ that the implicit finite element method (FEM) shows very poor convergence capabilities for the modeling of gridshells. The multiple interactions between the beams tend to lock the calculation, and the required precautions to ensure convergence are extremely time-consuming for this application. On the contrary, the dynamic relaxation method is robust to changes in the initial configuration.

In addition, its explicit formulation allows complex interaction between beams, and even the modeling of contact, if necessary. This versatility drove our choice for this dynamic relaxation method, in line with the previous works cited.

Originality of our approach

In the literature, one can find three different models dealing with slender beams in large displacement where the static equilibrium is calculated through a dynamic relaxation method:

In 1980, Wakefield¹² created a model with 6 degrees of freedom (DoFs) per node, later developed by Ong¹³ in 1992 and recently used by D’Amico et al.⁵ in 2014. This model is suited for the modeling of all slender beams.

In 2001, Adriaenssens¹⁴ reduced the previous model to 3 DoFs per node. This model does not take the torsion into account but applies in all cases where the torsion moments can be neglected. It is used in previous studies.^4,6,8
–10

In 2013, Barnes et al.⁷ improved the previous model by taking the torsion into account with only these 3 translation DoFs per node. Yet, this model only applies in cases where the beam stays in the vicinity of its rest configuration.

The novelty of our work is to improve these models by keeping the best of each of them. On one hand, according to Adriaenssens,¹⁴ the model (1) with 6 DoFs suffers from conditioning problems that could arise due to the coupling between the dynamic of the rotational and the translational DoFs. On the other hand, the model with 3 DoFs is not able to model properly the torsion in every configuration. In fact, the first model suffers from an excess of DoF, while the latter model lacks DoF. The solution has been to build a beam model solved by dynamic relaxation with 4 DoFs per node. 4 DoFs are indeed enough to model completely a slender beam in any configuration, taking torsion into account. It is in fact the minimal number of DoFs possible that enables this kind of modeling.

It is here worth underlining the need for taking torsion into account when modeling gridshells, as without it, the model with 3 DoFs would be sufficient. A first argument is the crucial role of torsion in the shape of beams with non-symmetrical cross sections, as shown in Figure 2: depending on how the beam is bent, out-of-plan displacement can be important.

Figure 2.

The need for modeling torsion: stable equilibriums of a beam with rectangular cross section in a model (a) with torsion or (b) without torsion.

A second argument is the negative influence of torsion on the creep rupture strength of composite beams submitted to bending, as shown in Kotelnikova-Weiler.¹⁵

By taking torsion into account, gridshells with non-symmetrical cross sections can then be precisely modeled, and, for all gridshells, even those with symmetrical cross sections, the torsion in the beams can be predicted, in order to avoid rupture under combined bending–torsion permanent stresses.

Hence, we built a model for slender beams, taking torsion into account, with 4 DoFs per node, solved by dynamic relaxation. This model can be used for numerous applications where modeling with finite displacements and multiple boundary conditions are used, ranging from the modeling of gridshells to the prediction of cable routing for the automotive industry.

Modeling large displacements of slender beams

The dynamic relaxation method

Before going further into the details of our model, let us rapidly explain the dynamic relaxation method, for the sake of completeness.

The dynamic relaxation method first appeared in the 1960s with Day¹⁶ and Otter.¹⁷ The goal of this method is to find a static equilibrium. To do so, starting from a given configuration with zero speed, a fictitious dynamic is simulated, where the system under study is immobilized each time it reaches a maximum of kinetic energy and is then released again with zero speed. By immobilizing the system, the kinetic energy of the system is dissipated. We then simulate a dynamic with “kinetic” damping.

For a conservative system, the maximum of kinetic energy is a minimum of potential energy, so that when we reach a kinetic energy peak on a trajectory, it is the local minimum of potential energy on this trajectory. Immobilizing the system at kinetic energy peaks is therefore an efficient way to converge toward the static equilibrium, defined as the minimum of potential energy. This is why the kinetic damping is chosen for this damped dynamic.

An example of such a fictitious dynamic is illustrated in Figure 3, where we look for the static equilibrium of a ball in an ellipsoidal bowl. In this example, the ball is stopped and released with zero speed each time it reaches a maximum speed.

Figure 3.

Static equilibrium of a ball in an ellipsoidal bowl, found by dynamic relaxation.

Let us now explicit the algorithm of the dynamic relaxation method for a system with $N$ nodes $x_{0}^{t}, \cdot, x_{N - 1}^{t}$ , with $M_{i}$ as the mass affected to the node $x_{i}^{t}$ . The time interval is denoted by ∆ $t$ .

Algorithm 1. Dynamic relaxation algorithm
1: initialization
2: convergence parameter: $E_{c}^{l i m}$
3: geometry: $x_{i} = x_{0}$
4: reset
5: reset velocity: $v_{i}^{t - \frac{Δ t}{2}} = 0$
6: compute lumped mass: $M_{i}$
7: loop
8: compute forces: $f_{i}^{t} (x_{i}^{t})$
9: compute accelerations: $a_{i}^{t} = \frac{f_{i}^{t}}{M_{i}}$
10: compute speeds: $v_{i}^{t + \frac{Δ t}{2}} = v_{i}^{t - \frac{Δ t}{2}} + Δ t a_{i}^{t}$
11: compute kinetic energy: $E_{c}^{t + \frac{Δ t}{2}} = \frac{1}{2} \sum M_{i} v_{i}^{t^{2}}$
12: if $E_{c}^{t + \frac{Δ t}{2}} > E_{c}^{t - \frac{Δ t}{2}}$ then ▹ no peak
13: compute positions: $x_{i}^{t + Δ t} = x_{i}^{t - Δ t} + Δ t v_{i}^{t + \frac{Δ t}{2}}$
14: go to step 7
15: else if $E_{c}^{t + \frac{Δ t}{2}} < E_{c}^{t - \frac{Δ t}{2}}$ and $E_{c}^{t - \frac{Δ t}{2}} > E_{c}^{l i m}$ then ▹peak
16: go to step 4
17: else ▹convergence
18: return
19: end if
20: end loop

The difficult part of this algorithm is to compute the forces at step 4, and it will be the purpose of the following work. At this step, we indeed need to deduce the internal and external forces $f_{i}^{t}$ acting on the beam, knowing its configuration $x_{i}^{t}$ , which is not straightforward. To achieve this task, we will therefore make use of Kirchhoff’s beam theory^1,2, as described in the following sections.

Modeling slender beams

This last introductive section is dedicated to the review of the models for slender beams available in the literature.^7,18
–26 The choices for the dynamic relaxation method and the 4 DoFs per node have indeed not been made alone, but validated by a careful analysis of the possible ways of modeling slender beams. This analysis, together with the previous considerations, will give all the foundations of the model: the underlying beam theory, the kinematic description, and the numerical tool.

The underlying beam theories

The beams under study being slender, two fundamental assumptions can be made:

The cross sections remain normal to the centerline

The centerline is inextensible

This last assumption is less common than the first one, but nevertheless relevant. As explained in Audoly et al.²⁷ and Starostin and Van der Heijden,²⁸ based on a scaling argument, two cases arise for slender beams: (1) either the centerline stretches and the bending and twisting forces become negligible compared to the axial forces; (2) either the centerline remains inextensible. As we are not interested in the first case, dealing with beams used as cables in tension, we will keep the inextensibility assumption.

With these assumptions, the Kirchhoff beam theory applies.^1,2 As in previous studies,^19
–24 we chose this theory. Ong¹³ and Theetten et al.²⁵ preferred the Euler–Bernoulli beam theory that drops the assumption of inextensibility, but as explained above, it is an unnecessary precaution for our application. Barnes et al.⁷ and Choe et al.²² chose homemade theories, avoided here, for reliability concerns. It is here worth noting that Kirchhoff’s beam theory is built assuming that the cross sections remain normal to the centerline, but they do not necessarily remain planar: the cross sections may warp in case of torsion.²⁷

The kinematic description

In order to describe the kinematic of a Kirchhoff beam, a description of the position and orientation of each cross section is needed. Barnes et al.⁷ only describe the position of the centerline (3 DoFs) but this is not enough. Bertails-Descoubes et al.^20,21 and Pai²⁴ describe the position of the first cross section, plus the rotation of each cross section with respect to the previous one (3 DoFs + 1 boundary condition). Previous works^18,19,23,25 describe the position of the centerline, plus the angle of rotation of each cross section around this centerline (4 DoFs). Previous works^13,22,23 describe the position of the centerline, plus the orientation of each cross section (6 DoFs). These three last kinematic descriptions are equally able to describe the Kirchhoff beams under study. We chose the one with 4 DoFs for two reasons. First, because it offers an explicit description of each cross section of the beam, whereas the one with 3 DoFs plus 1 boundary condition for each section is defined with respect to the first one, making the definition of other boundary conditions harder, for example at the last cross section. Second, because 4 DoFs is a minimum for the explicit description of each cross section of a Kirchhoff beam: more DoFs appear unnecessary.

The numerical tool

Different numerical methods can be used to find the static equilibrium of a Kirchhoff beam. As explained above, we chose the dynamic relaxation method, which simulates a kinetically damped dynamic. But the minimization of the potential energy by gradient descent²¹ or by the Newton method²³ could have been used too, as well as the static implicit FEM.²⁹ Our choice was yet driven by the will to use an explicit method and not an implicit one, because we are dealing with beams in multiple interactions (in gridshells for example). A beam in complex interaction has indeed very constrained displacements and shows a very tortuous load–displacement path, so that the time step of the numerical method cannot be large. As limiting the time step of an implicit method kills its efficiency, using an explicit method appeared the most reasonable solution, hence the choice for the dynamic relaxation method.

Summary

Now that all choices have been justified, let us sum them up:

Our slender beam model makes use of the Kirchhoff beam theory, based on the following assumptions: (1) the cross sections remain normal to the centerline and (2) the centerline is inextensible.

The kinematic description relies on 4 DoFs: 3 DoFs for the position of the centerline and 1 DoF for the angle of rotation of each section around the centerline.

The chosen numerical tool is the dynamic relaxation method.

These choices being made, the next section is dedicated to the construction of the model.

Solving Kirchhoff beam static equilibrium with the dynamic relaxation method

As seen in the previous section, in order to apply the dynamic relaxation method, the resultants of the internal forces have to be derived at each step, knowing the current configuration of the beam. This will be the purpose of the following sections (“Derivation of the elastic potential energy” and “Computation of the forces”).

Let us first precise how the configuration of the beam is described. As a Kirchhoff beam is inextensible, let us denote $s \in [0, L]$ as the curvilinear abscissa along a beam of length $L$ . The centerline can then be described by the curve $s \mapsto x (s) \in ℝ^{3}$ and the rotation of the sections around the centerline by an angle $s \mapsto θ (s) \in ℝ$ . This angle will be measured with respect to the Bishop frame, introduced later. So, the configuration of the beam can be described by 4 DoFs: $x$ and $θ$ .

The chosen strategy to derive the resultant of the internal forces from the configuration $(x, θ)$ is the following:

Derive the elastic potential energy as a function of $x$ and $θ$ (section “Derivation of the elastic potential energy”)

Differentiate the elastic potential energy with respect to $x$ and $θ$ to obtain the resultant of the internal forces (section “Computation of the forces”)

This approach is similar to the one chosen in Bergou et al.,¹⁹ but rather than a discrete formulation, we chose a continuous one, in order to have simpler expressions for the forces. Keeping a continuous formulation implies much more complex theoretical development, with the use of the functional derivative, but simplifies the expressions found in Bergou et al.,¹⁹ and therefore saves some simulation time.

After these two steps, the dynamic relaxation method can be applied (section “Implementation”).

Derivation of the elastic potential energy

In the following, the elastic potential energy will be denoted as $E_{p}$ . As writing $E_{p}$ as a function of $x$ and $θ$ is not straightforward, we will first present the classical expression of $E_{p}$ as a function of the curvatures and the twist in section “The elastic potential energy,” and $x$ and $θ$ will then be introduced in section “Introduction of the Bishop frame.”

The elastic potential energy

In this paragraph, a classic and basic way to describe the orientation of the sections is used. Let us define an orthonormal material frame $s \mapsto (d_{1} (s), d_{2} (s), d_{3} (s))$ attached to the sections of the beam, with $d_{3} (s) = x' (s)$ . Because this frame is orthonormal, there are three scalars here denoted by $κ_{1}$ , $κ_{2}$ , and $τ$ such that

\begin{array}{l} d_{1}' (s) = Ω (s) \times d_{1} (s) \\ d_{2}' (s) = Ω (s) \times d_{2} (s) \\ d_{3}' (s) = Ω (s) \times d_{3} (s) \end{array}

(1)

where Ω is the Darboux vector defined as

Ω (s) = κ_{1} d_{1} (s) + κ_{2} d_{2} (s) + τ d_{3} (s)

(2)

By definition, $κ_{1}$ , $κ_{2}$ , and $τ$ are, respectively, the curvatures and the twist of the beam (Figure 4).

Figure 4.

Definition of the curvatures and the twist.

Using these notations, and denoting $E$ as Young’s modulus; $G$ the shear modulus; $I_{1}$ and $I_{2}$ the principal moments of inertia; and $J$ the geometrical coefficient for torsion, the Kirchhoff beam theory²⁷ yields

E_{p} = \int_{0}^{L} \frac{E I_{1}}{2} κ_{1}^{2} + \frac{E I_{2}}{2} κ_{2}^{2} + \frac{G J}{2} τ_{1}^{2}

(3)

Note that this equation is derived under the assumption of small deformations, that is, slowly varying curvature and twist, excluding sharp-edged beams. This assumption is indeed verified for elastically bent rods.

Introduction of the Bishop frame

Based on the expression of $E_{p}$ in equation (3), the expression $E_{p} (x, θ)$ can be derived. To do so, a precise definition of $θ$ is first needed, coming along with the definition of the Bishop frame.³⁰ As mentioned above, $θ (s)$ is defined as the angle of rotation of the section attached to the frame ${d_{1} (s), d_{2} (s)}$ with respect to the Bishop frame ${u (s), v (s)}$ (Figure 5).

Figure 5.

Definition of θ $θ$ .

The Bishop frame ${u (s), v (s)}$ is a frame defined on the curve $s \mapsto x (s)$ by parallel transport of an initial frame ${u_{0}, v_{0}}$ chosen in $s = 0$ . A frame ${u_{0}, v_{0}}$ being chosen in $s = 0$ , parallel transport along a curve is defined as follows

\begin{array}{l} u' (s) = κ b (s) \times v (s) \\ v' (s) = κ b (s) \times u (s) \end{array}

(4)

where $κ b$ stands for the curvature binormal

\begin{array}{l} κ b (s) = x^{'} (s) \times x^{″} (s) \end{array}

(5)

Figure 6 shows two examples of Bishop frames built by parallel transport. We can see that the Bishop frame “rotates around the curvature binormal,” as equation (4) suggests.

Figure 6.

Examples of Bishop frames (left: curve made of arcs, right: helix).

It is worth underlining that the definition of the Bishop frame only depends on the position of the centerline $x$ , and on a choice of ${u_{0}, v_{0}}$ in $s = 0$ . So, knowing $x$ , $θ$ , and ${u_{0}, v_{0}}$ , it is possible to draw the beam, first by drawing the centerline x , then by reconstructing the Bishop frame ${u, v}$ by parallel transport of ${u_{0}, v_{0}}$ , and finally, by placing the sections with respect to the Bishop frame, thanks to the knowledge of $θ$ .

The Bishop frame is not a randomly chosen reference frame. The Darboux vector associated with the parallel transport $Ω = κ b$ is indeed normal to the centerline, meaning that $τ = 0$ for this frame: the Bishop frame is then a “zero twisting frame.” The use of this frame for the study of Kirchhoff beams has been initiated in Langer and Singer³¹ and continued in Bergou et al.,¹⁹ because it offers a convenient reference for the study of torsion. The Bishop frame not only simplifies the equations but also avoids the numerical instability encountered by Bertails-Descoubes et al.²¹ with the better known Frenet frame. It is therefore our best choice as a reference for the definition of $θ$ . $x$ and $θ$ being precisely defined, $E_{p} (x, θ)$ can now be derived from equation (3), by remarking that

\begin{array}{l} κ_{1} = - x^{″} \cdot d_{2} = - x^{″} \cdot (c o s θ v - s i n θ u) \\ κ_{2} = x^{″} \cdot d_{1} = x^{″} \cdot (c o s θ u + s i n θ v) \\ τ = θ^{'} \end{array}

(6)

Let us denote by $ω$ the curvature binormal projected in the material frame, $B$ the bending stiffness matrix, and $β$ the torsional stiffness ( $B$ and $β$ can vary along the beam)

\begin{array}{l} ω = (\begin{matrix} x^{″} \cdot d_{1} \\ x^{″} \cdot d_{2} \end{matrix}) = (\begin{matrix} x^{″} \cdot (c o s θ u + s i n θ v) \\ x^{″} \cdot (c o s θ v - s i n θ u) \end{matrix}) \\ B = (\begin{matrix} E I_{1} & 0 \\ 0 & E I_{2} \end{matrix}), β = G J \end{array}

(7)

Then, $E_{p}$ written as a function of $x$ and $θ$ is simply

E_{p} = \frac{1}{2} \int_{0}^{L} ω^{T} B ω + β {(θ^{'})}^{2}

(8)

This expression can be generalized for non-straight undeformed configurations. Denoting by $\bar{ω}$ and $\bar{θ}$ the properties of the beam in the undeformed configuration, $E_{p}$ can be rewritten as follows

E_{p} (x, θ) = \frac{1}{2} \int_{0}^{L} {(ω - \bar{ω})}^{T} B (ω - \bar{ω}) + β {(θ^{'} - \bar{θ^{'}})}^{2}

(9)

Derivation of the internal forces an moments

In this section, the quasi-static out-of-balance internal forces and moments acting on the beam are calculated by differentiating bending and torsion elastic energies relatively to the DoFs of the system $(x, θ)$ . This forces and moments will be used later in a damped explicit dynamic procedure to solve the equilibrium of the system.

Differentiating $E_{p} (x, θ)$ with respect to $x$ and $θ$ yields the linear resultants of, respectively, the flexural and torsional internal forces and moments. As $x$ and $θ$ are functions (of $s$ ), the differentiation of $E_{p} (x, θ)$ requires the use of a functional derivative. In a functional derivative, a function is differentiated with respect to a function, instead of differentiating a function with respect to a variable as in a usual derivative. In section “The functional derivative,” the functional derivative will be defined. The linear resultants of the torsional (section “Computation of the moments”) and flexural (section “Computation of the forces”) internal forces are then calculated.

The functional derivative

For $h \in ℝ$ and $ϕ$ any arbitrary function, the functional derivative $s \mapsto \frac{d E_{p}}{d θ} (s)$ of $E_{p}$ with respect to $θ$ is such that

E_{p} (x, θ + h ϕ) - E_{p} (x, θ) = h \int_{0}^{L} \frac{d E_{p}}{d θ} ϕ + o (h^{2})

(10)

For $h \in ℝ$ and $ϕ$ any arbitrary function, the functional derivative $s \mapsto \frac{d E_{p}}{d x} (s)$ of $E_{p}$ with respect to $x$ is such that

E_{p} (x + h ϕ, θ) - E_{p} (x, θ) = h \int_{0}^{L} \frac{d E_{p}}{d x} ϕ + o (h^{2})

(11)

With these definitions, and denoting by $f$ the linear resultant of the flexural internal forces and $m$ the linear resultant of the torsional moments, we have

\begin{array}{l} \begin{array}{l} f = - \frac{d E_{p}}{d x} = - \frac{\partial E_{p}}{\partial ω} \frac{\partial ω}{\partial x} \\ m = - \frac{d E_{p}}{d θ} = - \frac{\partial E_{p}}{\partial θ} - \frac{\partial E_{p}}{\partial ω} \frac{\partial ω}{\partial θ} \end{array} \end{array}

(12)

Note the dependence of $ω$ in both $x$ and $θ$ . Consequently, the total derivative of $E_{p}$ with respect to the variables $x$ and $θ$ are given by the chain rule.

In practice, $f$ is a linear force acting on the centerline and $m$ is a linear moment acting on the orientation of the cross sections.

The timescale separation assumption

In order to simplify the computations and to save some simulation time, Bergou et al.¹⁹ noted that the speed of the twist waves increases when the typical size $l$ of the cross sections decreases and, on the contrary, that the speed of the bending waves increases with the ratio $l / λ$ , $λ$ being the wavelength. For the problems studied here, we have both vanishing cross sections and large wavelength-bending modes, so that the following timescale separation assumption can be made: the propagation of the twist waves is instantaneous compared to the one of the bending waves. So, the forces $f$ and the moment $m$ act on two different timescales.

Consequently, on the timescale of action of the force $f$ on the centerline which drives the bending waves, the twist waves propagate instantaneously. Thus, the twisting kinematics regarding the DoF $θ$ will be considered as a quasi-static phenomena in the computation of $f$ . In other words, when calculating $f$ , the following assumption will be made

\forall s \in [0, L], m (s) = - \frac{d E_{p}}{d θ} (s) = 0

(13)

This assumption will later lead to simplifications in the computation of the internal forces.

Computation of the moments

The quasi-static out-of-balance internal linear moment acting on the beam is calculated by differentiating bending and torsion elastic energies relatively to the twisting DoF of the system ( $θ$ )

\begin{matrix} m = - \frac{d E_{p}}{d θ} \\ = - \frac{d}{d θ} (\frac{1}{2} \int_{0}^{L} {(ω - \bar{ω})}^{T} B (ω - \bar{ω}) + β {(θ^{'} - {\bar{θ}}^{'})}^{2}) \end{matrix}

(14)

In order to compute the functional derivative as in equation (10), an expression for $ω (x, θ + h ϕ)$ is needed. Starting from equation (7) and using a finite expansion of $s i n$ and $c o s$ , we find

ω (x, θ + h ϕ) = ω (x, θ) - R ω (x, θ) h ϕ + o (h^{2})

(15)

where $R = (\begin{matrix} 0 & - 1 \\ 1 & 0 \end{matrix})$ is the 90° rotation matrix.

Then, $E_{p} (x, θ + h ϕ)$ can be computed

\begin{matrix} E_{p} (x, θ + h ϕ) = E_{p} (x, θ) + h \int_{0}^{L} β (θ - \bar{θ})^{'} ϕ^{'} \\ + h \int_{0}^{L} (ω - \bar{ω}) T B R ω ϕ \\ + o (h^{2}) \end{matrix}

(16)

By integration by parts, and denoting by $δ_{s}$ the Dirac delta function centered in $s$

\begin{matrix} \int_{0}^{L} β (θ - \bar{θ})^{'} ϕ^{'} = - \int_{0}^{L} (β (θ - \bar{θ})^{'})^{'} ϕ \\ + \int_{0}^{L} β (θ - \bar{θ})^{'} (δ_{L} - δ_{0}) ϕ \end{matrix}

(17)

From equations (10), (14), (16), and (17), we can conclude on the expression of $s \mapsto m (s)$

\begin{matrix} m (s) = (β (θ - \bar{θ})^{'})^{'} - β (θ - \bar{θ})^{'} (δ_{L} - δ_{0}) \\ + {(ω - \bar{ω})}^{T} B R ω \end{matrix}

(18)

Note that the second term, involving the Dirac delta functions, vanishes everywhere but at beam’s ends.

Computation of the forces

The quasi-static out-of-balance internal linear force acting on the beam is calculated by differentiating bending and torsion elastic energies relatively to the displacement DoFs of the system ( $x$ )

f = - \frac{d E_{p}}{d x} = - \frac{d}{d x} (\frac{1}{2} \int_{0}^{L} {(ω - \bar{ω})}^{T} B (ω - \bar{ω}))

(19)

In order to compute this partial derivative as in equation (11), an expression for $ω (x + h ϕ)$ is required

ω (x + h ϕ) = (\begin{matrix} (x + h ϕ)^{″} \cdot d_{1} (x + h ϕ) \\ (x + h ϕ)^{″} \cdot d_{2} (x + h ϕ) \end{matrix})

(20)

Now, an expression for $d_{1} (x + h ϕ)$ is required too. As $h ϕ$ is a small variation of the centerline, $d_{1} (x + h ϕ)$ can be written as a small variation of $d_{1} (x)$ , in the frame ${d_{1}, d_{2}, x'}$

d_{1} (x + h ϕ) = d_{1} (x) + Δ Ψ d_{2} (x) + Δ Φ x' + o (h^{2})

(21)

where $Δ Ψ$ and $Δ Φ$ are the first-order terms, respectively, equal to the small angles of rotation around $x'$ and $d_{2}$ . Then

\begin{matrix} (x + h ϕ)^{″} \cdot d_{1} (x + h ϕ) = x^{″} \cdot d_{1} (x) + h ϕ^{″} \cdot d_{1} (x) \\ + Δ Φ x^{″} \cdot d_{2} (x) + o (h^{2}) \end{matrix}

(22)

A similar relation is found for $(x + h ϕ)^{″} \cdot d_{2} (x + h ϕ)$ , so that equation (20) can be rewritten as follows

\begin{matrix} ω (x + h ϕ) = ω (x) + h (\begin{matrix} d_{1}^{T} \\ d_{2}^{T} \end{matrix}) \cdot ϕ^{″} - Δ Ψ R ω (x) + o (h^{2}) \end{matrix}

(23)

Let us now focus on $t \in [0, L], t \mapsto Δ Ψ (t)$ . $Δ Ψ$ is the angle of rotation of $d_{1}$ around $x'$ when the centerline undergoes a variation $h ϕ$ . As $θ$ remains constant in this partial derivative with respect to $x$ , the angle between $d_{1}$ and the Bishop vector $u$ remains constant too. Then, the angle of rotation of $d_{1}$ around $x'$ is equal to the angle of rotation of $u$ around $x'$ . So, $Δ Ψ$ is the angle of rotation of the Bishop frame around the centerline when the centerline undergoes a variation $h ϕ$ .

Parallel transport being a propagation of frame from $s = 0$ , $Δ Ψ (t)$ is the cumulated angle of rotation of $u (x + h ϕ)$ around $x'$ during its parallel transport along the portion $[0, t]$ of the curve $x + h ϕ$ . Looking at the definition of parallel transport (4), we note that the Bishop vector $u$ rotates with a rotation rate $κ b$ during its parallel transport along the curve $x$ .^19,32 Then, $u (x + h ϕ)$ rotates with a rotation rate vector $κ b (x + h ϕ)$ during its parallel transport along $x + h ϕ$ . Hence, the cumulated angle of rotation of $u (x + h ϕ)$ around $x'$ is

Δ Ψ = \int_{0}^{t} κ b (x + h ϕ) \cdot x'

(24)

Remembering that $κ b = x^{'} \times x^{″}$ gives

κ b (x + h ϕ) = κ b (x) + h (ϕ^{'} \times x^{″} + ϕ^{″} \times x^{'}) + o (h^{2})

(25)

Substituting equation (25) in (24) yields

Δ Ψ = h \int_{0}^{t} (ϕ^{'} \times x^{″} + ϕ^{″} \times x^{'}) \cdot x^{'} + o (h^{2})

(26)

Using the properties of the triple product, we have $ϕ^{'} \times x^{″} \cdot x^{'} = - x^{'} \times x^{″} \cdot ϕ^{'} = - κ b \cdot ϕ^{'}$ . Then, choosing $L$ as upper integration limit, we get

Δ Ψ (t) = h \int_{0}^{L} κ b ϕ^{'} 1_{[0, t [} + o (h^{2})

(27)

where $1_{[0, t [}$ denotes the indicator function of $[0, t [$ . Now that all terms appearing in the expression of $ω (x + h ϕ)$ given in equation (23) have been made explicit, and using the symmetry of $B$ , we can derive the first equation for $E_{p} (x + h ϕ)$

\begin{matrix} E_{p} (x + h ϕ) = E_{p} (x) + h \int_{0}^{L} {(ω - \bar{ω})}^{T} B (\begin{matrix} d_{1}^{T} \\ d_{2}^{T} \end{matrix}) \cdot ϕ^{″} \\ - \int_{0}^{L} {(ω - \bar{ω})}^{T} B R ω Δ Ψ + o (h^{2}) \end{matrix}

(28)

In this last equation, $d E_{p} / d x$ cannot be explicitly identified as presented in equation (11). In equations (29)–(32), some rearrangements will therefore be done, in order to make an inner product with $ϕ$ appearing explicitly. For a synthetic presentation of the next three equations, let us denote by $I = \int_{0}^{L} {(ω - ω)}^{T} B R ω Δ Ψ$ the last integral of equation (28) and $a = (ω - \bar{ω})^{T} B R ω$ . First, using equation (27), we have

\begin{array}{l} I = h \int_{0}^{L} a (t) \int_{0}^{L} κ b (s) \cdot ϕ^{'} (s) 1_{[0, t [} (s) d s d t + o (h^{2}) \end{array}

(29)

Second, using Fubini’s theorem and remarking that $1_{[0, t [} (s) = 1_{] s, L]} (t)$ yields

\begin{array}{l} I = h \int_{0}^{L} κ b (s) \cdot ϕ^{'} (s) \int_{0}^{L} a (t) 1_{] s, L]} (t) d t d s + o (h^{2}) \end{array}

(30)

Third, noting that $a 1_{] s, L]} (t) = a$ and integrating by parts gives

\begin{array}{l} I = h {\int_{0}^{L} (κ b (s) \int_{s}^{L} a)}^{'} \cdot ϕ (s) d s + o (h^{2}) \end{array}

(31)

Substituting equation (32) in (28), and integrating by parts the first integral of equation (28) two times yields the following convenient expression for $E_{p} (x + h ϕ)$

\begin{matrix} E_{p} (x + h ϕ) = E_{p} (x) + h {\int_{0}^{L} ({(ω - \bar{ω})}^{T} B (\begin{matrix} d_{1}^{T} \\ d_{2}^{T} \end{matrix}))}^{″} \cdot ϕ \\ - h {\int_{0}^{L} (κ b \int_{s}^{L} {(ω - \bar{ω})}^{T} B R ω)}^{'} \cdot ϕ \\ + o (h^{2}) \end{matrix}

(32)

With equation (32), we can conclude on the value of $f$

\begin{matrix} f = - \frac{d E_{p}}{d x} = - {({(ω - ω)}^{T} B (\begin{matrix} d_{1}^{T} \\ d_{2}^{T} \end{matrix}))}^{″} \\ + {(κ b \int_{s}^{L} {(ω - \bar{ω})}^{T} B R ω)}^{'} \end{matrix}

(33)

A last simplification can be made on this equation, thanks to the timescale separation assumption. We have indeed $\forall s \in [0, L], d E_{p} / d θ = 0$ . Using the value of $d E_{p} / d θ$ found in equation (18) then yields

\begin{array}{l} {(ω - \bar{ω})}^{T} B R ω = - (β (θ - \bar{θ})^{'})^{'} + β (θ - \bar{θ})^{'} (δ_{L} - δ_{0}) \end{array}

(34)

Substituting equation (34) in (33) finally gives

\begin{array}{l} f = - ({(ω - \bar{ω})}^{T} B (\begin{matrix} d_{1}^{T} \\ d_{2}^{T} \end{matrix})) ″ + (β (θ^{'} - {\bar{θ}}^{'}) κ b)^{'} \end{array}

(35)

Note that, unlike in equation (17), the integration by parts for equations (31) and (32) does not involve $δ_{0}$ or $δ_{L}$ . We indeed preferred to nullify these terms because they would involve $κ b (0)$ or $κ b (L)$ which are poorly defined at the ends of a curve. As a consequence, equation (35) is valid for the open interval $] 0, L [$ , but not for $s = 0$ or $s = L$ , where the nullified terms involving the Dirac delta functions would be non-zero. The boundaries will be treated in section “Boundary conditions.”

Implementation

Knowing the resultants of the internal forces, the dynamic relaxation method can be applied to the Kirchhoff beam, as described in section “The dynamic relaxation method.” In this section, the geometry and relevant quantities will first be discretized (section “Discretization”). Then, solutions to impose the inextensibility (section “Inextensibility”), the boundary conditions (section “Boundary conditions”), and the external forces (section “External forces”) will be presented. Finally, the dynamic relaxation method adapted to Kirchhoff beams will be made explicit (section “Dynamic relaxation method”).

Discretization

The discrete beam will be described with $N$ nodes $x_{0}, \dots, x_{N - 1}$ and $N$ cross-section orientations $θ_{0}, \dots, θ_{N - 1}$ at these nodes. The segments will be denoted by $e_{i} = x_{i + 1} - x_{i}$ and their lengths $l_{i} = ∥ e_{i} ∥$ . These notations are illustrated in Figure 7. The curvature binormal $κ b_{i}$ at node $x_{i}$ is defined by considering that the curvature $κ_{i} = ∥ κ b ∥_{i}$ is equal to the inverse of the radius of the circle passing through $x_{i - 1}$ , $x_{i}$ , and $x_{i + 1}$

\begin{array}{l} κ b_{i} = \frac{2}{l_{i - 1} l_{i}} \frac{e_{i - 1} \times e_{i}}{∥ e_{i - 1} + e_{i} ∥} \end{array}

(36)

Figure 7.

Portion of the discrete beam.

Let us denote by ${u_{0}, v_{0}}$ the Bishop frame chosen in $x_{0}$ . As parallel transport is invariant along a straight line, the Bishop frame is constant along each segment of the discrete beam. Let us denote by ${u_{i / i + 1}, v_{/ i + 1}}$ the constant Bishop frame on the segment $e_{i}$ , with $u_{0 / 1} = u_{0}$ and $v_{0 / 1} = v_{0}$ , and by ${u_{i}, v}$ the Bishop frames defined at the nodes $x_{i}$ . Denoting by $R_{a}$ the rotation matrix of axis $a$ and angle $∥ a ∥$ , the parallel transport is defined as follows

(\begin{matrix} u_{i} = R_{α_{i - 1}^{i} b_{i}} u_{i - 1 / i} \\ u_{i / i + 1} = R_{α_{i}^{i} b_{i}} u_{i} \end{matrix}, (\begin{matrix} v_{i} = R_{α_{i - 1}^{i} b_{i}} v_{i - 1 / i} \\ v_{i / i + 1} = R_{α_{i}^{i} b_{i}} v_{i} \end{matrix}

(37)

where

\begin{array}{l} b_{i} = κ b_{i} / κ_{i} \\ α_{i - 1}^{i} = \arcsin (l_{i - 1} κ_{i} / 2) \\ α_{i}^{i} = \arcsin (l_{i} κ_{i} / 2) \end{array}

(38)

These notations are illustrated in Figure 8. This definition of ${u_{i}, v_{i}}$ is such that ${u_{i}, v_{i}}$ is orthogonal to the circle passing through $x_{i - 1}$ , $x_{i}$ , and $x_{i + 1}$ . ${u_{i}, v_{i}}$ and $θ_{i}$ being defined, ${d_{1 i}, d_{2 i}}$ is given by

\begin{array}{l} d_{1 i} = \cos θ_{i} u_{i} + \sin θ_{i} v_{i} \\ d_{2 i} = - \sin θ_{i} u_{i} + \cos θ_{i} v_{i} \end{array}

(39)

and

ω = (\begin{matrix} κ b_{i} \cdot d_{2 i} \\ - κ b_{i} \cdot d_{1 i} \end{matrix})

(40)

Figure 8.

Discrete parallel transport.

Finally, denoting the following terms by

\begin{array}{l} Δ θ'_{i + \frac{1}{2}} = ((θ_{i + 1} - {\bar{θ}}_{i + 1}) - (θ_{i} - {\bar{θ}}_{i})) / l_{i} \\ M_{i} = - (ω_{i} - {\bar{ω}}_{i}) B_{i} (\begin{matrix} d_{1 i}^{T} \\ d_{2 i}^{T} \end{matrix}) \\ r_{i} = (ω_{i} - {\bar{ω}}_{i})^{T} B_{i} R ω_{i} \end{array}

(41)

We can write the internal discretized forces and moments

\begin{matrix} f_{i}^{i n t} = \frac{1}{l_{i - 1}} M_{i - 1} - (\frac{1}{l_{i}} + \frac{1}{l_{i - 1}}) M_{i} + \frac{1}{l_{i}} M_{i + 1} \\ + \frac{β_{i} κ b_{i} + β_{i + 1} κ b_{i + 1}}{2} Δ θ'_{i + \frac{1}{2}} \\ - \frac{β_{i - 1} κ b_{i - 1} + β_{i} κ b_{i}}{2} Δ θ'_{i - \frac{1}{2}} \end{matrix}

(42)

\begin{matrix} m_{i}^{i n t} = \frac{l_{i - 1} + l_{i}}{2} r_{i} + \frac{β_{i} + β_{i + 1}}{2} Δ θ'_{i + \frac{1}{2}} \\ - \frac{β_{i - 1} + β_{i}}{2} Δ θ'_{i - \frac{1}{2}} \end{matrix}

(43)

Inextensibility

In order to enforce the inextensibility, we chose to use penalty forces. Denoting by $L_{i}$ the length of $e_{i}$ in the rest configuration and by $A_{i}$ its axial stiffness, the penalty force in $x_{i}$ is given by

\begin{array}{l} f_{i}^{p e n} = A_{i} (\frac{l_{i}}{L_{i}} - 1) \frac{e_{i}}{l_{i}} - A_{i - 1} (\frac{l_{i - 1}}{L_{i - 1}} - 1) \frac{e_{i - 1}}{l_{i - 1}} \end{array}

(44)

The maximal time step ∆ $t$ ensuring stability for the dynamic relaxation method decreases when the stiffness of the system increases, so that the axial stiffness $A_{i}$ affects the efficiency of the numerical method. The value of $A_{i}$ must be sufficient to enforce the inextensibility, but low enough to guarantee an efficient simulation. Starting with $A_{i}$ equal to the real axial stiffness of the beam, this value can then be adjusted. Note that if $A_{i}$ is equal to the real axial stiffness, $f_{i}^{p e n}$ is equal to the axial forces.

Boundary conditions

As observed when deriving the internal forces, some quantities ( $x^{″}$ and $κ b$ for example) are poorly defined at the extremities of the beams. To cope with this problem, the discretized beam can be artificially extended. The fictitious extension of the beam will vary depending on the boundary conditions:

For a ball-joint connection, the extension is chosen aligned with the last segment so that $x^{″} = 0$

For a clamped beam, the extension is chosen symmetrical to the last segment with respect to the clamping plane. In this case, $x^{″}$ can be non-zero.

These situations are illustrated in Figure 9. Using these fictitious extensions, the internal forces can be defined at the real ends of the beam, and balanced with the external forces induced by the connections. This trick avoids to define the internal forces at the ends of the considered curve, as the beam is strictly included in this curve.

Figure 9.

Extensions for the boundary conditions.

External forces

In order to be added to the internal forces, the external forces have to be expressed in the same way as the internal forces, namely, forces $f_{i}^{e x t}$ acting on the centerline and moments $m_{i}^{e x t}$ acting on the orientation of the cross sections. Generally, external forces are expressed in the form of a couple force/moment ${({\tilde{f}}_{i}^{e x t}, {\tilde{M}}_{i}^{e x t})}_{i \in [1, N - 1]}$ . The goal of this paragraph is to project this couple ${({\tilde{f}}_{i}^{e x t}, {\tilde{M}}_{i}^{e x t})}_{i \in [1, N - 1]}$ in the space of the couples ${(f_{i}^{e x t}, m_{i}^{e x t})}_{i \in [1, N - 1]}$ accepted by our model. To do so, the moment ${\tilde{M}}_{i}^{e x t}$ will be decomposed into two components: a tangential component ${({\tilde{M}}_{i}^{e x t})}_{∥}$ and an orthogonal component ${({\tilde{M}}_{i}^{e x t})}_{⊥}$ , as illustrated in Figure 10. Denoting by $t_{i} = u_{i} \times v_{i}$ the tangent at $x_{i}$ , they are defined as

\begin{array}{l} {({\tilde{M}}_{i}^{e x t})}_{∥} = {\tilde{M}}_{i}^{e x t} \cdot t_{i} \\ {({\tilde{M}}_{i}^{e x t})}_{⊥} = {\tilde{M}}_{i}^{e x t} - {({\tilde{M}}_{i}^{e x t})}_{∥} \end{array}

(45)

Figure 10.

Projection of the external moment.

The tangential component tends to rotate the cross sections, so that we simply have

m_{i}^{e x t} = ({\tilde{M}}_{i}^{e x t})_{∥}

(46)

The orthogonal component locally rotates the centerline, so that it can be transformed into couples of forces acting on the adjacent points. The forces are then given by

\begin{matrix} f_{i}^{e x t} = {\tilde{f}}_{i}^{e x t} + \frac{1}{l_{i - 2} + l_{i - 1}} \frac{1}{l_{i - 1}} {({\tilde{M}}_{i - 1}^{e x t})}_{⊥} \times e_{i - 1} \\ + \frac{1}{l_{i - 1} + l_{i}} \frac{1}{l_{i - 1}} {({\tilde{M}}_{i}^{e x t})}_{⊥} \times e_{i - 1} \\ + \frac{1}{l_{i - 1} + l_{i}} \frac{1}{l_{i}} {({\tilde{M}}_{i}^{e x t})}_{⊥} \times e_{i} \\ + \frac{1}{l_{i} + l_{i + 1}} \frac{1}{l_{i}} {({\tilde{M}}_{i + 1}^{e x t})}_{⊥} \times e_{i} \end{matrix}

(47)

With this choice of external forces, the moment induced in $x_{i}$ by the forces $f_{i - 1}^{e x t}$ , $f_{i}^{e x t}$ , and $f_{i + 1}^{e x t}$ is equal to ${({\tilde{M}}_{i}^{e x t})}_{⊥}$ , so that ${(f_{i}^{e x t}, m_{i}^{e x t})}_{i \in [1, N - 1]}$ is indeed equivalent to ${({\tilde{f}}_{i}^{e x t}, {\tilde{M}}_{i}^{e x t})}_{i \in [1, N - 1]}$ .

Dynamic relaxation method

In section “The dynamic relaxation method,” we proposed a simple description of the dynamic relaxation method, with only the translational DoFs $x$ . The dynamic of the cross sections has now to be added through $θ$ . As suggested by the timescale separation assumption, the dynamics of $x$ and $θ$ will be computed in parallel, based on two different kinetic energies:

1. Translational kinetic energy

E_{c}^{x} = \frac{1}{2} \sum_{i = 0}^{N - 1} M_{i} {\dot{x}}_{i}^{2}

(48)

where $M_{i}$ denotes the mass assigned to the node $x_{i}$ .

2. Rotational kinetic energy

E_{c}^{θ} = \frac{1}{2} \sum_{i = 0}^{N - 1} I_{i} {\dot{θ}}_{i}^{2}

(49)

where $I_{i}$ denotes the inertia assigned to the section located in $x_{i}$ .

Then, the algorithm can be modified as follows (see Algorithm 2): With this algorithm for the dynamic relaxation method, the two dynamics of $x$ and $θ$ evolve in parallel, and the convergence is reached when the translational kinetic energy is less than the specified threshold $E_{c, l i m}^{x}$ . This threshold is usually chosen as a fraction of the first peak of kinetic energy, so that we consider that the equilibrium is reached when a sufficient amount of kinetic energy has been dissipated. Note that due to the timescale separation assumption, if the equilibrium is reached for the translational dynamic, it is also the case for the rotational dynamic.

Algorithm 2. Dynamic relaxation algorithm
1: initialization
2: define convergence parameter: $E_{c}^{l i m}$
3: define rest configuration: $({\bar{x}}_{i}^{0}, {\bar{θ}}_{i}^{0})$
4: define initial configuration: $(x_{i}^{0}, θ_{i}^{0})$
5: reset
6: reset translational velocity: ${\dot{x}}_{i}^{t - \frac{Δ t}{2}} = 0$
7: reset rotational velocity: ${\dot{θ}}_{i}^{t - \frac{Δ t}{2}} = 0$
8: compute lumped mass: $M_{i}$
9: loop
10: compute forces: $f_{i}^{t} = f_{i}^{i n t, t} + f_{i}^{p e n, t} + f_{i}^{e x t, t}$
11: compute translational accelerations: $a_{i}^{t} = f_{i}^{t} / M_{i}$
12: compute translational speeds: ${\dot{x}}_{i}^{t + \frac{Δ t}{2}} = {\dot{x}}_{i}^{t - \frac{Δ t}{2}} + Δ t a_{i}^{t}$
13: procedure detect translational kinetic energy peak
14: compute kinetic energy: $E_{c}^{x, t + \frac{Δ t}{2}} = \frac{1}{2} \sum M_{i} {({\dot{x}}_{i}^{t + \frac{Δ t}{2}})}^{2}$
15: if $E_{c}^{x, t + \frac{Δ t}{2}} > E_{c}^{x, t - \frac{Δ t}{2}}$ then ▹no peak
16: compute positions: $x_{i}^{t + Δ t} = x_{i}^{t - Δ t} + Δ t {\dot{x}}_{i}^{t + \frac{Δ t}{2}}$
17: else if $E_{c}^{t + \frac{Δ t}{2}} > E_{c}^{l i m}$ then ▹peak
18: go to step 6
19: else ▹convergence
20: return
21: end if
22: end procedure
23: compute moments: $m_{i}^{t} = m_{i}^{i n t, t} + m_{i}^{e x t, t}$
24: compute angular accelerations $α_{i}^{t} = m_{i}^{t} / I_{i}$
25: compute rotational speeds: ${\dot{θ}}_{i}^{t + \frac{Δ t}{2}} = {\dot{θ}}_{i}^{t - \frac{Δ t}{2}} + Δ t α_{i}^{t}$
26: procedure detect rotational kinetic energy peak
27: compute kinetic energy: $E_{c}^{θ, t + \frac{Δ t}{2}} = \frac{1}{2} \sum I_{i} {({\dot{θ}}_{i}^{t + \frac{Δ t}{2}})}^{2}$
28: if $E_{c}^{θ, t + \frac{Δ t}{2}} > E_{c}^{θ, t - \frac{Δ t}{2}}$ then ▹no peak
29: compute angles: $θ_{i}^{t + Δ t} = θ_{i}^{t - Δ t} + Δ t {\dot{θ}}_{i}^{t + \frac{Δ t}{2}}$
30:
31: else
32: go to step 7
33: end if
34: end procedure
35: end loop

Let us finish this paragraph by the choice of the masses and inertia $M_{i}$ and $I_{i}$ . As we are looking for a static equilibrium, they have no influence on the equilibrium found. But as explained in Lewis,³³ in order to ensure stability, the following Courant–Friedrichs–Lewy (CFL) condition holds

\frac{k δ t^{2}}{M} < 2

(50)

where $k$ is the stiffness of the system and $M$ the mass. This relation is true for a system composed of one element, and can be generalized to more complex systems. Based on equation (50), we indeed used the general formula $M = k δ t^{2}$ , with $k$ equal to the sum of the stiffnesses of the segments converging in $x_{i}$ . Noting that the bending stiffness is negligible compared to the axial stiffness for the evaluation of $M_{i}$ , we chose

\begin{array}{l} M_{i} = δ t^{2} (\frac{A_{i - 1}}{l_{i - 1}} + \frac{A_{i}}{l_{i}}) \\ M_{i} = δ t^{2} β_{i} (\frac{1}{l_{i - 1}} + \frac{1}{l_{i}}) \end{array}

(51)

These choices were validated by our simulations: we encountered instabilities for lower values, and the use of greater values unnecessarily slowed the simulations.

Given these values, a typical evolution of the two kinetic energies during a simulation is shown in Figure 11. It can be observed that the dynamic in $θ$ is much faster than the one in $x$ , so that the timescale separation assumption indeed applies. Now that all modeling details have been given, the model can be validated. In the next section, we will show the competitive precision and speed of our algorithm.

Figure 11.

Evolution of the kinetic energies.

Validation of the model

The validation of the model is done through a comparison with the finite element software Abaqus. Many situations have been studied. First, the model has been validated in simple cases, that is to say when only one elastic energy is at stake. These situations offered satisfying results, but they are not developed here. Second, the model has been validated, thanks to the simulation of situations with bending, torsion, and bending–torsion coupling. For these studies, the beams have an anisotropic section. More precisely, they have a rectilinear rest-state configuration, and a uniform rectangular cross section. The study is presented here through two examples. For these two examples, the initial configuration is very different from the equilibrium one, which gives an indicator of the robustness of the model. The studied beam has the following properties:

Length L = 10 m

Rectangular section: 4 × 6 cm

Young’s modulus: E = 25 GPa

Shear modulus: G = 10 GPa

The computations performed with Abaqus use B31 Timoshenko beam elements. The geometric non-linearities are obviously allowed. The simulation is explicit because implicit simulations do not accept large displacements.

Case of “thwarted” flexion

This kind of flexion happens, for example, when an anisotropic beam is bent along the direction of highest inertia. For a rectangular cross section, the bending starts classically in-plane. When increasing the curvature, the geometry evolves suddenly out-of-plane. This is a bifurcation. Numerically, the in-plane configuration can be stable even if it is physically unstable. This situation is illustrated in Figure 12. The position in-plane is in black; the equilibrium configuration is in gray. In the case of study, the beam is clamped at both ends. The tangents at the ends are set to do an π $/ 3$ angle from the straight line passing through both ends. The distance between the two ends is equal to $d = 2 m$ . The accuracy of the model is monitored through two quantities, one geometric characterized by the $h$ parameter, and one mechanic—the maximum flexion stress along the beam in direction 1, $M_{1}^{f l e x i o n}$ . The results are presented in Figure 13. They are normalized with the converged values obtained through FEM simulations $(h = 1.46 m; M_{1}^{f l e x i o n} = 2763 Nm)$ . The results are satisfying in terms of accuracy and in terms of mechanics: the computed values for the bending moment in the other direction and the torsion moment are close to the Abaqus ones. Moreover, the computation time is good. For 100 nodes, a Python implementation of the algorithm converges in 1 min with a convergence criterion of 10^-10 times the value of energy of the first translation energy peak. Abaqus, with a standard setting, also converges in 1 min. A recent C++ implementation of our algorithm ran about two orders of magnitude faster, demonstrating the potential of this numeric tool to perform quasi-real-time structural analysis studies on models with a few hundreds of nodes. To sum up this study, in this situation, where flexion and torsion are coupled, the accuracy of the model is good, in terms of geometry and mechanics.

Figure 12.

Out-of-plane bending of a rectangular beam.

Figure 13.

Evolutions of out-of-plane displacement and flexion moment.

Case of “buckled” torsion

To complete the study of the beam model, another configuration is studied. This configuration is the one of a rectilinear anisotropic beam, put in torsion and then clamped at both ends. Finally, this twisted beam is buckled by bringing its ends closer, without modification of the clamped ends. This situation is illustrated in Figure 14. In this figure, the results of the FEM anlaysis are presented on the left side, whereas the results of the developed model are presented on the right side. Once again, both the geometry and the stresses are monitored, through the $u_{r}^{\max}$ (radial maximal displacement) and ${(G J τ)}^{\max}$ (maximal torsion moment) parameters. The results are presented below, for a 100-elements model.

	Abaqus	Dynamic relaxation
$u_{r}^{\max}$	$2.34 m$	$2.33 m$
${(G J τ)}^{\max}$	$4269 Nm$	$4266 Nm$

Figure 14.

Comparison of the equilibrium of “buckled” torsion beam obtained with Abaqus (left) and the model (right).

Once again, a good concordance is reported between the results of the model and the results obtained through FEM analysis.

Results

The 4-DoF model taking torsion and flexion into account and developed to obtain the equilibrium configuration of slender inextensible beams has been compared to a reference FEM. The comparison has been performed in several situations. Some of them with a high torsion–flexion coupling are presented in this article. It appears that our model offers a good accuracy in terms of geometry and mechanic stresses. With a C++ implementation, the computation time is in the order of magnitude of 1 s for 100 nodes, with an expected linear increase, thanks to the explicit formulation, opening for quasi-real-time applications. It is important to note that the implemented algorithm is stable and robust as it can deal with large displacements. Finally, it has been implemented in an efficient way to do form-finding, which is a large advantage in comparison with Abaqus.

Conclusion and outlooks

The beam model developed in this article relies on a 4-DoF kinematics that captures both bending and torsion behaviors of slender beams. These 4 DoFs (three for the nodal position plus one for the orientation of the cross section around the centerline) offer the minimal way to represent a Kirchhoff beam. According to Bergou et al.,¹⁹ cross-section orientations are given relatively to a reference frame (the Bishop frame), thanks to the angular field $θ$ .

Bending and torsion elastic energies are then calculated: the quasi-static out-of-balance internal forces and moments acting on the centerline are obtained by differentiating these energies relatively to the DoFs ( $x$ , $θ$ ). The internal axial force is taken into account trough a penalty force to enforce the inextensibility assumption during the simulation process. Finally, internal forces and moments are combined with external applied loads and support reactions to compute the overall quasi-static out-of-balance internal forces (acting on the centerline) and moments (acting on the cross sections).

The equilibrium configuration of the system is obtained through a damped dynamic analysis with an explicit time integration scheme, where the system is set in motion as long as the overall quasi-static internal forces and moments are out-of-balance. In this article, an adapted dynamic relaxation procedure with kinetic damping is presented to deal with both $x$ and $θ$ DoFs. This procedure relies on the timescales separation assumption, which leads to simplifications in the equations and is verified since twisting waves commonly travel much faster than bending waves in slender beams.

The discretization of the continuous model is presented and implemented. A benchmark is conducted on a test case, comparing this new model with a reference model embedded in a well-known commercial finite element software (Abaqus). Computation times are very similar to the ones obtained with Abaqus. The results are satisfying in terms of accuracy. Displacements and stresses converge for situations presenting high bending–torsion coupling at about 1% of the reference solution. A major advantage of the model developed here is that it can be implemented in a dedicated environment to the study of slender beams. This environment offers very quick pre-processing which matters for form-finding processes in design stages that require fast back-and-forth exchanges between architects and structural engineers. For instance, a gridshell, which is a structure made of interconnected slender beams,^4,34–36 has been modeled as illustrated in Figure 15.

Figure 15.

Structural skeleton of the gridshell built in Paris in 2013.

Indeed, a gridshell can be seen as a structure made of independent beams subjected to external forces and moments at each connection. The internal forces in the beams can also be computed. These forces and moments can be implemented in the model to simulate the whole structure. Then, it becomes possible to use the modeling for different purposes as modeling a full-scale structure or compare the influence of the choice of beams. For example, Figure 16 presents the equilibrium configuration of a gridshell with square cross sections (left—no bending–torsion coupling because the beams have square cross section and are rectilinear in their rest-state configuration) and the equilibrium of a gridshell having the same geometry at the flat stage, but made of rectangular cross sections (right—in this case, the coupling between bending and torsion exists). The two gridshells are very different from a geometric point of view, and also in terms of stress. The torsion moment has been illustrated. The characteristic values are significantly different for the two types of beams. Finally, the model implemented will offer large developments in modeling structures made of interconnected slender beams.

Figure 16.

Computation of whole gridshells, without (left) and with (right) flexion–torsion coupling.

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) received no financial support for the research, authorship, and/or publication of this article.

References

Love

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

Dill

EH.

Kirchhoff’s theory of rods. Arch Hist Exact Sci 1992; 23.

Lienhard

Alpermann

Gengnagel

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

Baverel

Caron

Tayeb

. Gridshells in composite materials: construction of a 300 m² forum for the Solidays’ festival in Paris. Struct Eng Int 2012; 22(3): 408–414.

D’Amico

Kermani

Zhang

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

Adriaenssens

Barnes

Williams

A new analytic and numerical basis for the form-finding and analysis of spline and gridshell structures. In: Kumar

Topping

BHV

(eds) Computing developments in civil and structural engineering. Edinburgh: Civil-Comp Press, 1999, pp. 83–91, http://opus.bath.ac.uk/13096/

Barnes

Adriaenssens

Krupka

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

Douthe

Baverel

Design of nexorades or reciprocal frame systems with the dynamic relaxation method. Comput Struct 2009; 87(21–22): 1296–1307.

Douthe

Baverel

Caron

JF.

Formfinding of a grid shell in composite materials. J IASS 2006; 47: 53–62.

10.

Douthe

Caron

Baverel

Gridshell in composite materials: towards large span shelters. Eur J Environ Civil Eng 2009; 13(9): 1083–1093.

11.

Bouhaya

Optimisation structurelle des gridshells. PhD Thesis, Université Paris-Est, 2010, https://tel.archives-ouvertes.fr/tel-00583409/

12.

Wakefield

Dynamic relaxation analysis of pretensioned networks supported by compression arches. PhD Thesis, City, University London, London, 1980.

13.

Ong

CF.

Computer-aided design of tension structures. PhD Thesis, City, University London, London, 1992.

14.

Adriaenssens

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

15.

Kotelnikova-Weiler

Optimisation mécanique et énergétique d’enveloppes en matériaux composites pour les bâtiments. PhD Thesis, Université Paris-Est, 2012, https://pastel.archives-ouvertes.fr/pastel-00807510

16.

Day

An introduction to dynamic relaxation. Engineer 1965; 219: 218–221.

17.

Otter

Computations for prestressed concrete reactor pressure vessels using dynamic relaxation. Nucl Struct Eng 1965; 1(1): 61–75.

18.

Bergou

Audoly

Vouga

. Discrete viscous threads. In: Proceedings of the ACM transactions on graphics (TOG), 2010, pp. 1–10, http://dl.acm.org/citation.cfm?id=1778853

19.

Bergou

Wardetzky

Robinson

. Discrete elastic rods. In: Proceedings of the ACM SIGGRAPH, 2008, pp. 1–12, http://portal.acm.org/citation.cfm?doid=1508044

20.

Bertails-Descoubes

Audoly

Cani

. Super-helices for predicting the dynamics of natural hair. ACM T Graphics 2006; 25(3): 1180.

21.

Bertails-Descoubes

Audoly

Querleux

. Predicting natural hair shapes by solving the statics of flexible rods. Eurographics 2005; 3–6.

22.

Choe

Choi

HS.

Simulating complex hair with robust collision handling. In: Proceedings of the 2005 ACM SIGGRAPH/Eurographics symposium, 2005, pp. 153–160, http://portal.acm.org/citation.cfm?doid=1073368.1073389

23.

Grégoire

Schömer

Interactive simulation of one-dimensional flexible parts. Comput Aided Design 2007; 39(8): 694–707.

24.

Pai

DK.

STRANDS: interactive simulation of thin solids using Cosserat models. Eurographics 2002; 21(3): 347–352.

25.

Theetten

Grisoni

Andriot

. Geometrically exact dynamic splines. Comput Aided Design 2008; 40: 35–48.

26.

Meier

Popp

Wall

WA.

An objective 3D large deformation finite element formulation for geometrically exact curved Kirchhoff rods. Comput Method Appl M 2014; 278: 445–478.

27.

Audoly

Amar

Pomeau

Elasticity and geometry: from hair curls to the nonlinear response of shells. Oxford University Press, 2010.

28.

Starostin

van der Heijden

The equilibrium shape of an elastic developable Möbius strip. Proc Appl Math Mech 2007; 7: 2020115–2020116.

29.

Dassault Systèmes. SIMULIA: Abaqus theory manual. Providence, RI: Dassault Systèmes, 2013.

30.

Bishop

There is more than one way to frame a curve. Washington, DC: Mathematical Association of America, 1975.

31.

Langer

Singer

Lagrangian aspects of the Kirchhoff elastic rod. Soc Ind Appl Math Rev 1996; 38(4):605–618.

32.

Alexander

Antman

The ambiguous twist of love. Q Appl Math 1982; 40(1):83–92.

33.

Lewis

Tension structures: form and behaviour. London: Telford Thomas, 2003.

34.

Burkhardt

Otto

IL 13: Multihalle Mannheim. Stuttgart: Universität Stuttgart, 1976.

35.

Harris

Romer

Kelly

. Design and construction of the Downland Gridshell. Build Res Inform 2003; 31(6):427–454.

36.

Du Peloux

Tayeb

Baverel

. Construction of a large composite gridshell structure: a lightweight structure made with pultruded glass fibre reinforced polymer tubes. Struct Eng Int 2016; 26(2): 160–167.