# Formulation of a 4-DoF torsion/bending element for the formfinding of elastic gridshells

Lionel Du Peloux, Frédéric Tayeb, Baptiste Lefevre, Olivier Baverel, Jean-François Caron

► To cite this version:

Lionel Du Peloux, Frédéric Tayeb, Baptiste Lefevre, Olivier Baverel, Jean-François Caron. Formulation of a 4-DoF torsion/bending element for the formfinding of elastic gridshells. IASS Symposium 2015: "Future Visions", Aug 2015, Amsterdam, Netherlands. hal-01199049

HAL Id: hal-01199049

https://hal.science/hal-01199049

Submitted on 15 Sep 2015

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
# Formulation of a 4-DoF torsion/bending element for the formfinding of elastic gridshells

Lionel DU PELOUX*, Frédéric TAYEBª, Baptiste LEFEVREª,
Olivier BAVERELª, Jean-François CARONª

* UMR Navier Ecole des Ponts ParisTech
6-8 av. Blaise Pascal, Cité Descartes, 77455, Marne-la-Vallée, France
lionel.dupeloux@gmail.com, www.thinkshell.fr/en

ª UMR Navier Ecole des Ponts ParisTech

## Abstract

The paper presents a completely novel approach to model elastic gridshells with a 4-DoF element based on recent advances in the field of hair modelling [1]. This element, based on Kirchhoff's beam theory, can account for both bending and torsion behaviours. The reduction from 6 to 4 degrees of freedom is achieved with an appropriate curve framing introduced by R. Bishop [2]. The resulting model is fast and efficient. Its accuracy has been validated on test cases. It opens new outlooks for the design of elastic gridshells.

**Keywords:** gridshell, torsion, bending, twist, Kirchhoff, elastic rod, formfinding, composite, GFRP

## 1. Introduction

Elastic gridshells offer an efficient way to build freeform shapes from an initially flat structure. Indeed, it is extremely relevant to materialize curved shapes using planar and regular objects that will be deformed to obtain the desired volume. Freeform is always complex and costly to build whereas standard planar and orthogonal construction is much more conventional and economic. This capacity to "form the form" efficiently is of peculiar importance in the current context where morphology is a predominant component of modern architecture (F. Gehry, Z. Hadid, etc.) and envelopes appear to be the neuralgic point for building environmental performance.

### 1.1. A brief history on elastic gridshells

The "proof of concept" was made by Frei Otto in 1974 when he built the emblematic Multihall of Mannheim – a wooden elastic gridshell of 7500m² – thanks to the expert knowledge he had developed on physical models and photogrammetry [3]. More recently, the development of numerical models has
<figure><img src="image_2.png" /><figcaption>Figure 1: Inside view of the Cathedral</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2: Grid erection with large deflections</figcaption></figure>

encouraged the emergence of few projects made of cardboard (Hanovre 2000), timber (Downland 2004, Savill 2006) or composite material (Solidays 2011, Créteil 2013).

**1.2. Recent developments**

For the last decade the Navier laboratory has developed a research program on elastic gridshells [4] focusing on both the use of composite materials [5] and the development of efficient design and structural analysis methods. These developments have been validated by the construction of two prototypes of 150 m² and the elaboration of specific design software.

More recently, its expert knowledge has led the laboratory to contribute to 2 large-scale prototypes of about 400m² each: the *Forum Café of Solidays* in 2011, which was the first gridshell in composite material to house public [6], [7]; and the *Ephemeral Cathedral of Créteil* in 2013 (Figure 1), a religious edifice still in service [8], [9]. Those large-scale prototypes represent a first in the building industry, which still shows excessive apprehension for the use of non-traditional materials such as composites, especially when it comes to structural applications, apart from the field of reinforcement [10].

**1.3. Recent developments**

The major difficulty in designing such structures is to compute the deformed geometry of the initially flat grid (Figure 2) and the corresponding stresses due to beams bending and twisting. 3-DoF elements have yet been formulated to achieve such computations using dynamic relaxation integration scheme [11], [12] but they don't account for torsional effects, anisotropic sections or complex joint cinematic, which could be very limiting. 6-DoF elements have also been formulated but they suffer from convergence instabilities [13]–[15].

**2. Elastic gridshells**

Literally, the word « gridshell » refers to grids behaving like shells: from a mechanical point of view that means stresses acting on the structure are mainly transmitted through compression and traction. These structures can cross large-span with very little material.
**Structural Typology.** Their mechanical behaviour is very similar to the one of real shells even if the material is discrete and located in a grid more or less open. In spite of that, gridshells benefit from the same advantages as the ones showed by an eggshell: they can cross large span using a low amount of material. Their stiffness is mainly linked to their double-curved shape.

**Material Flexibility for Structural Rigidity.** In this field of application, composite materials like glass fibre reinforced polymer (GFRP) could favourably replace wood, where both resistance and bending ability of the material is sought. The stiffness of the structure does not derive from the intrinsic material rigidity but principally from its geometric curvature. Ideally, the composite profiles are produced by pultrusion, an economic continuous moulded process. The standardization of the process guaranties very stable material and mechanical properties. It frees designers from the painful problematic of wood joining and wood durability.

**Erection Process.** Usually, the grid morphology is not trivial and leads to design numerous costly and complex joints. To overcome this issue, an original and innovative erection process was developed. It takes advantage of the flexibility inherent to slender elements. A regular planar grid made of long continuous linear members is built on the ground. The elements are pinned together so the grid has no in-plane shear stiffness and can accommodate large-scale deformations during erection. Then, the grid is bent elastically to its final shape (Figure 3.3). Finally, the grid is frozen in the desired shape with a third layer of bracing members and the structure becomes a shell (Figure 3.5).

<figure><img src="image_4.png" /><figcaption>Figure 3: Building stages</figcaption></figure>

## 3. A novel 4-DoF torsion/bending element

In this section a novel element with 4 degrees of freedom accounting for torsion and bending behaviours is presented. The beam is considered in Kirchhoff's theory framework, so that it is supposed to be inextensible and its sections are supposed to remain orthogonal to the centreline during deformation. The reduction from the classic 6-DoF model to this 4-DoF model is achieved by an appropriate curve-angle representation based on a relevant curve framing. Energies are then formulated and leads to internal forces and moments acting on the beam. The static equilibrium is deduced from a damped fictitious dynamic with an adapted dynamic relaxation algorithm.
## 3.1. Kirchhoff rod

The geometric configuration of the rod is described by its centerline $\mathbf{x}(s)$ and its cross sections. The centerline is parameterized by its arc-length. Cross sections orientations are followed along the centerline by their material frame $\{\mathbf{d}_1(s), \mathbf{d}_2(s), \mathbf{d}_3(s)\}$ which is an adapted orthonormal moving frame aligned to section's principal axes of inertia. Here, "adapted" means $\mathbf{d}_3(s) = \mathbf{x}'(s) = \mathbf{t}(s)$ is aligned to the centerline's tangent. In the literature, this description is also known as a Cosserat Curve.

### 3.1.1. Inextensibility.

Note the previous description is only valid for inextensible rods in order to follow material points by their arc-length indifferently in their rest or deformed configuration. As explained in [16], this hypothesis is usually relevant for slender beams. Indeed, in practice, if a slender member faces substantial axial strain the bending behaviour would become negligible due to the important difference between axial and bending stiffness. The length of the rod will be denoted L and the arc-length s will vary (with no loss of generality) in [0, L].

### 3.1.2. Euler-Bernoulli

Strains are supposed to remain small so that material frame remains orthogonal to the centerline in the deformed configuration. Thus, differentiating the conditions of orthonormality leads to the following differential equations governing the evolution of $\{\mathbf{d}_1(s), \mathbf{d}_2(s), \mathbf{d}_3(s)\}$ along the centerline:

$$ 
\begin{aligned}
\mathbf{d}_i \cdot \mathbf{d}_j &= 0 \Rightarrow \mathbf{d}_i', \mathbf{d}_j + \mathbf{d}_i \cdot \mathbf{d}_j' = 0 & \mathbf{d}_1'(s) &= \tau(s)\mathbf{d}_2(s) - \kappa_2(s)\mathbf{d}_3(s) \\
\mathbf{d}_i \cdot \mathbf{d}_i &= 1 \Rightarrow \mathbf{d}_i', \mathbf{d}_i = 0 & \mathbf{d}_2'(s) &= \kappa_1(s)\mathbf{d}_3(s) - \tau(s)\mathbf{d}_1(s) \\
& & \mathbf{d}_3'(s) &= \kappa_2(s)\mathbf{d}_1(s) - \kappa_1(s)\mathbf{d}_2(s)
\end{aligned}
$$

### 3.1.3. Darboux vector.

Those equations can be formulated with the Darboux vector of the chosen material frame, which represents the rotational velocity of the frame along $\mathbf{x}(s)$:

$$ 
\begin{aligned}
\mathbf{d}_1'(s) &= \mathbf{\Omega}_m(s) \times \mathbf{d}_1(s) \\
\mathbf{d}_2'(s) &= \mathbf{\Omega}_m(s) \times \mathbf{d}_2(s) \\
\mathbf{d}_3'(s) &= \mathbf{\Omega}_m(s) \times \mathbf{d}_3(s)
\end{aligned}
\quad , \quad \mathbf{\Omega}_m(s) = \begin{bmatrix} \kappa_1 \\ \kappa_2 \\ \tau \end{bmatrix}(s)
$$

$\kappa_1(s), \kappa_2(s), \tau(s)$ represent respectively the rate of rotation of the material frame around the axis $\mathbf{d}_1(s), \mathbf{d}_2(s), \mathbf{d}_3(s)$.

### 3.1.4. Curvatures & Twist

The material curvatures are denoted $\kappa_1(s)$ and $\kappa_2(s)$ and represent the rod's flexion in the principal planes respectively normal to $\mathbf{d}_1(s)$ and $\mathbf{d}_2(s)$. The material twist is denoted $\tau(s)$ and represents the section's rate of rotation around $\mathbf{d}_3(s)$. Those scalar functions measure directly the strain as defined in Kirchhoff's theory (Figure 4). Recall that the Frenet frame $\{\mathbf{t}(s), \mathbf{n}(s), \mathbf{b}(s)\}$ defines the osculating plane and the total curvature ($\kappa$) of a spatial curve:

$$ \mathbf{t}' = \kappa \mathbf{n} \qquad \kappa = \|\mathbf{t}'\| \qquad \mathbf{b} = \mathbf{t} \times \mathbf{n} = \frac{1}{\kappa} \mathbf{t} \times \mathbf{t}' $$
<figure><img src="image_5.png" /><figcaption>Figure 4: Curvatures and twist</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: Curve-angle representation</figcaption></figure>

To describe the osculating plane in which lies the bending part of the deformation, let's introduce the curvature binormal $\kappa\mathbf{b} = \mathbf{t} \times \mathbf{t}'$, a vector of direction $\mathbf{b}$ and norm $\kappa$. The osculating plane is normal to $\kappa\mathbf{b}$.

### 3.1.5. Elastic energy

Kirchhoff's theory assigns an elastic energy to beams according to their strain [16]. In this theory, a beam is supposed to be inextensible. Thus the elastic energy ($E_p$) only accounts for torsion and bending behaviours and is given by:

$$E_p = E_b + E_t = \frac{1}{2} \int_0^L EI_1(\kappa_1 - \overline{\kappa_1})^2 + EI_2(\kappa_2 - \overline{\kappa_2})^2 + \frac{1}{2} \int_0^L GJ(\tau - \overline{\tau})^2$$

Here, $\overline{\kappa_1}$, $\overline{\kappa_2}$, $\overline{\tau}$ denote the natural curvature and twist of the rod in the rest position (no stress).

## 3.2. Curve-Angle representation

The previous paragraph has shown how the elastic potential energy of a rod can be computed following both its centerline and its cross sections orientations, which represents a model with 6-DoF: 3 for centerline positions and 3 for cross section orientations. Following [1], let's introduce a reduced coordinate formulation of the rod that accounts for only 4-DoF. This reduction of DoF relies on the concept of zero-twisting frame which gives a reference frame with zero twist along a given centerline. Thus, cross section orientations {$d_1(s), d_2(s), d_3(s)$} can be tracked only by the measure of an angle $\theta$ from this reference frame denoted {$u(s), v(s), d_3(s)$} (Figure 5). An alternative could be to parameterize the global rotations of local material frame and to compute the rotation needed to align two successive frames along the curve's tangent.

### 3.2.1. Zero-twisting frame

Zero-twisting frame, also known as Bishop frame, was introduced by Bishop in 1964. Bishop remarked that there was more than one way to frame a curve [2]. Indeed, for a given curve, any orthonormal moving frame would satisfy the following differential equations, where $k_1(s), k_2(s), \tau(s)$ are scalar functions that define completely the moving frame:

$$\begin{bmatrix} \mathbf{e}_1' \\ \mathbf{e}_2' \\ \mathbf{e}_3' \end{bmatrix} (s) = \begin{bmatrix} 0 & \tau(s) & -k_2(s) \\ -\tau(s) & 0 & k_1(s) \\ k_2(s) & -k_1(s) & 0 \end{bmatrix} \begin{bmatrix} \mathbf{e}_1 \\ \mathbf{e}_2 \\ \mathbf{e}_3 \end{bmatrix} (s)$$
For instance, a Frenet frame {$t(s), n(s), b(s)$} is a frame which satisfies $k_1 = 0$. Note that this frame suffers from major disadvantages: it is undefined where the curvature vanishes and it flips at inflexion points.

A Bishop frame {$t(s), u(s), v(s)$} is a frame which satisfies $\tau = 0$. By construction, this frame has no angular velocity (i.e. no twist) around the curve's tangent: $u \cdot v' = u' \cdot v = 0$. Its evolution along the curve is described by the corresponding Darboux vector: $\Omega_b(s) = \kappa b = t \times t'$. Remark that $\Omega_b(s)$ only depends on the centerline and is well defined even when the curvature vanishes. Thus, by the help of $\Omega_b(s)$, it's possible to transport a given vector $e$ along the centerline with no twist: $e' = \kappa b \times e$. This is called *parallel transport*.

### 3.2.2. Bending strains

Let's compute the bending strains $\kappa_1(s)$ and $\kappa_2(s)$ regarding the geometric configuration of the rod:

$$ \begin{aligned} \kappa b \cdot \mathbf{d}_1 &= (\mathbf{d}_3 \times \mathbf{d}_3'). \mathbf{d}_1 = (\mathbf{d}_1 \times \mathbf{d}_3). \mathbf{d}_3' = -\mathbf{d}_2 \cdot \mathbf{d}_3' = \kappa_1 \\ \kappa b \cdot \mathbf{d}_2 &= (\mathbf{d}_3 \times \mathbf{d}_3'). \mathbf{d}_2 = (\mathbf{d}_2 \times \mathbf{d}_3). \mathbf{d}_3' = \mathbf{d}_1 \cdot \mathbf{d}_3' = \kappa_2 \end{aligned} \quad \omega = \begin{bmatrix} \kappa_1 \\ \kappa_2 \end{bmatrix} = \begin{bmatrix} \kappa b \cdot \mathbf{d}_1 \\ \kappa b \cdot \mathbf{d}_2 \end{bmatrix} = \begin{bmatrix} -\mathbf{x}'' \cdot \mathbf{d}_2 \\ \mathbf{x}'' \cdot \mathbf{d}_1 \end{bmatrix} $$

Where $\omega(s)$ is the vector of material curvatures expressed on material frame axes {$\mathbf{d}_1(s), \mathbf{d}_2(s)$}.

### 3.2.3. Torsion strain

Let's compute the twist or torsion strain $\tau(s)$ regarding the geometric configuration of the rod. Decomposing the material frame on the bishop frame gives:

$$ \begin{bmatrix} \mathbf{d}_1 \\ \mathbf{d}_2 \end{bmatrix} = \begin{bmatrix} \cos \theta \mathbf{u} + \sin \theta \mathbf{v} \\ -\sin \theta \mathbf{u} + \cos \theta \mathbf{v} \end{bmatrix} \quad \tau = \mathbf{d}'_1 \cdot \mathbf{d}_2 = (\theta' \mathbf{d}_2 + \kappa b \times \mathbf{d}_1) \cdot \mathbf{d}_2 = (\theta' \mathbf{d}_2 - \kappa_2 \mathbf{d}_3) \cdot \mathbf{d}_2 = \theta' $$

Finally, remark that the twist is given directly by the variation of $\theta$. Note that the Frenet frame does not lead to a correct evaluation of the twist.

## 3.3. Elastic energy

Introducing theta, the elastic energy can be written as follow:

$$ E_p = E_b + E_t = \frac{1}{2} \int_0^L (\omega - \bar{\omega})^T B (\omega - \bar{\omega}) + \frac{1}{2} \int_0^L \beta (\theta' - \bar{\theta}')^2 \quad B = \begin{bmatrix} EI_1 & 0 \\ 0 & EI_2 \end{bmatrix}, \quad \beta = GJ $$

## 3.4. Inextensibility

Recall that the rod is supposed to be inextensible in the Kirchhoff's theory. Thus, there is no stretching energy associated with an axial strain. However, this constraint will be enforced via a penalty energy, which in practice is somehow very similar as considering an axial stiffness in the beam...

The assumption that $s$ is an arc length parameterization of the centreline leads to consider the rod as

## 3.5. Time-scale assumption

Following [1], it is relevant to assume that the propagation of twist waves is instantaneous compared to the one of bending waves. Thus, internal forces $f^{\text{int}}$ and moment $m^{\text{int}}$ acts on two different timescales in the rod dynamic. Thus on the timescale of action of the force $f^{\text{int}}$ on the center line, driving the
bending waves, the twist waves propagate instantaneously, so that $\forall s \in [0, L]$, $\delta E_p / \delta \theta = 0$ for the computation of $f^{\text{int}}$. This assumption may not be enforced, as in [17], but leads to simpler and faster computations.

### 3.6. Gradients

Internal torsional moments and forces acting on the rod are classically obtained by differentiating the potential energy of the system with respect to $\theta$ and $\mathbf{x}$. Here, the calculus is a bit tricky as far as the differentiation takes place in function spaces. After a brief reminder on *functional derivative*, the main results of the calculations of the energy derivatives are given.

#### 3.6.1. Derivative in function spaces

The usual notion of derivative can be extended in function spaces [18]. It is known as the Fréchet or Gâteaux derivative and is defined and computed as explained in the followings:

Let V and B be two Banach spaces and $U \subset V$ an open subset of V. Let $f: U \to f(u)$ be a function of $U^W$. $f$ is said to be *Gâteaux differentiable* at $u_0 \in U$ if there exists a continuous linear operator $Df(u_0) \in L(V,W)$ such that:

$$ \forall h \in U, \quad \lim_{h \to 0} \frac{f(u_0 + \lambda h) - f(u_0)}{\lambda} = \frac{d}{d\lambda} f(u_0 + \lambda h) \Big|_{\lambda=0} = Df(u_0) \cdot h $$

Or, equivalently:

$$ \forall h \in U, \quad f(u_0 + \lambda h) = f(u_0) + \lambda Df(u_0) \cdot h + o(h) \quad , \quad \lim_{h \to 0} \frac{o(h)}{\|h\|} = 0 $$

The dot denotes the evaluation of $Df(u_0)$ at $h$ and must not be confused with a product. The gradient is associated to the derivative in the particular case of *scalar functions*, also known as *functionals*.

Let H be a Hilbert space with its inner product denoted $\langle \cdot ; \cdot \rangle$. Let $U \subset H$ be an open subset of H. Let $F: U \to F(u)$ be a scalar function of $U^R$. The *gradient* of $F$ is the map: $x \to (\text{grad } F)(x)$ of $U^H$ such that:

$$ \forall h \in U, \quad \langle (\text{grad } F)(x); h \rangle = Df(x) \cdot h $$

On $L_2$ function space of square-integrable functions, this definition leads to the following notion of functional derivative or gradient:

$$ \forall h \in U, \quad Df(x) \cdot h = \langle (\text{grad } F)(x); h \rangle = \int (\text{grad } F)h = \int \frac{\delta f}{\delta x} h $$

Note that the chain rule is still valid for this notion of derivative and will be important to do the calculations:

$$ D(f \circ g)(x) = Df(g(x)) \circ Dg(x) $$
### 3.6.2. Derivative of theta with respect to x

It's important to notice that the centreline **x** is independent of θ. However, the conversion does not hold. In fact, a variation of the centerline would cause the cross section orientations to be perturbed and thus the potential energy to be modified.

The derivative of θ with respect to **x** is evaluated by the change of write [19] in the curve:

$$D_x\theta(s).h_x = \int_0^L \frac{\delta\theta}{\delta x}(s)h_x = \int_0^s \frac{-x \times (x + h_x)'}{1 + x' \cdot (x + h_x)'} (x'' + (x + h_x)')' = \int_0^L (\kappa b(\delta_s - \delta_0) - (1 - H_s)\kappa b').h_x$$

Where $H_s$ and $\delta_s$ are respectively the Heaviside step function and the Dirac function positioned at arc-length s. In other words, this means that a variation $h_x$ of **x** causes a variation $\int_0^L \frac{\delta\theta}{\delta x}(s)h_x$ of θ at s. Note that this is an integrated quantity, all along the centerline.

### 3.6.3. Derivative of the material curvatures vector with respect to x

Let's express the variation of the material frame after a variation $\lambda h_x$ of **x**:

$$\begin{aligned} d_1[\mathbf{x} + \lambda \mathbf{h}_x] &= d_1[\mathbf{x}] + \Delta\psi d_2[\mathbf{x}] + \Delta\phi d_3[\mathbf{x}] \\ d_2[\mathbf{x} + \lambda \mathbf{h}_x] &= d_2[\mathbf{x}] - \Delta\psi[\mathbf{x}] d_1[\mathbf{x}] - \Delta\phi d_3[\mathbf{x}] \end{aligned}$$

Firstly one can remark that the variation $\Delta\phi$ is of first order in $\lambda$ (recall that $d_3 = x'$):

$$\begin{aligned} d_1[\mathbf{x} + \lambda \mathbf{h}_x].(\mathbf{x} + \lambda \mathbf{h}_x)' &= 0 \Leftrightarrow d_1[\mathbf{x} + \lambda \mathbf{h}_x].\mathbf{x}' = \Delta\phi = -\lambda d_1[\mathbf{x} + \lambda \mathbf{h}_x].\mathbf{h}_x' = -\lambda d_1[\mathbf{x}].\mathbf{h}_x' + o(\lambda) \\ d_2[\mathbf{x} + \lambda \mathbf{h}_x].(\mathbf{x} + \lambda \mathbf{h}_x)' &= 0 \Leftrightarrow d_2[\mathbf{x} + \lambda \mathbf{h}_x].\mathbf{x}' = -\Delta\phi = -\lambda d_2[\mathbf{x} + \lambda \mathbf{h}_x].\mathbf{h}_x' = -\lambda d_2[\mathbf{x}].\mathbf{h}_x' + o(\lambda) \end{aligned}$$

The variation $\Delta\psi$ was computed before and is nothing but:

$$\Delta\psi = -\lambda \int_0^L \frac{\delta\theta}{\delta x}(s)h_x$$

Thus, the variation of the material curvatures vector after a variation $\lambda h_x$ of **x** can be calculated as (recall that $d_3 \cdot x'' = 0$):

$$\begin{aligned} \omega[\mathbf{x} + \lambda \mathbf{h}_x] &= \begin{bmatrix} -(\mathbf{x} + \lambda \mathbf{h}_x)'' \cdot d_2[\mathbf{x} + \lambda \mathbf{h}_x] \\ (\mathbf{x} + \lambda \mathbf{h}_x)'' \cdot d_1[\mathbf{x} + \lambda \mathbf{h}_x] \end{bmatrix} \\ &= \begin{bmatrix} -(\mathbf{x} + \lambda \mathbf{h}_x)'' \cdot (d_2[\mathbf{x}] + \lambda \int_0^L \frac{\delta\theta}{\delta x}(s)h_x) d_1[\mathbf{x}] + \lambda(d_2[\mathbf{x}].h_x)d_3[\mathbf{x}] \\ (\mathbf{x} + \lambda \mathbf{h}_x)'' \cdot (d_1[\mathbf{x}] - \lambda \int_0^L \frac{\delta\theta}{\delta x}(s)h_x) d_2[\mathbf{x}] - \lambda(d_1[\mathbf{x}].h_x)d_3[\mathbf{x}] \end{bmatrix} \\ &= \omega[\mathbf{x}] + \lambda \begin{bmatrix} -d_2[\mathbf{x}]^T \cdot h_x'' - \lambda \int_0^L \frac{\delta\theta}{\delta x}(s)h_x \begin{bmatrix} x'' \cdot d_1[\mathbf{x}] \\ x'' \cdot d_2[\mathbf{x}] \end{bmatrix} - \lambda \begin{bmatrix} d_2[\mathbf{x}].h_x' \\ d_1[\mathbf{x}].h_x' \end{bmatrix} (x'' \cdot d_3[\mathbf{x}]) + o(\lambda) \end{aligned}$$

Which finally leads to the derivative of the material curvatures vector, denoting $R_{\pi/2}$ the counter-clockwise rotation matrix of angle $\frac{\pi}{2}$:
$$ \mathbf{D}_x \omega(\mathbf{x}). \mathbf{h}_x = \begin{bmatrix} -\mathbf{d}_2[\mathbf{x}]^T \\ \mathbf{d}_1[\mathbf{x}]^T \end{bmatrix} \cdot \mathbf{h}_x'' + \left( \int_0^L \frac{\delta\theta}{\delta x}(s) \mathbf{h}_x \right) R_{\pi/2} \omega[\theta] \qquad R_{\pi/2} = \begin{bmatrix} 0 & -1 \\ 1 & 0 \end{bmatrix} $$

Note that because $\omega$ is a vector function, its derivative is no more a scalar function. By linearity of derivation, product and integration, the operator $\mathbf{D}_x\omega(\mathbf{x})$ is linear regarding $\mathbf{h}_x$.

### 3.6.4. Derivative of the material curvatures vector with respect to $\theta$

Let's express the variation of the material frame after a variation $h_\theta$ of $\theta$:

$$ \mathbf{d}_1[\theta + h_\theta] = \mathbf{d}_1[\theta] + \sin(h_\theta)\mathbf{d}_2[\theta] - (1 - \cos h_\theta)\mathbf{d}_1[\theta] = \mathbf{d}_1[\theta] + h_\theta\mathbf{d}_2[\theta] + o(h_\theta) $$

$$ \mathbf{d}_2[\theta + h_\theta] = \mathbf{d}_2[\theta] - \sin(h_\theta)\mathbf{d}_1[\theta] - (1 - \cos h_\theta)\mathbf{d}_2[\theta] = \mathbf{d}_2[\theta] - h_\theta\mathbf{d}_1[\theta] + o(h_\theta) $$

Thus, it's possible to calculate the variation of the material curvatures vector after a variation $h_\theta$ of $\theta$:

$$ \omega[\theta + h_\theta] = \begin{bmatrix} -\mathbf{x}'' \cdot \mathbf{d}_2[\theta + h_\theta] \\ \mathbf{x}'' \cdot \mathbf{d}_1[\theta + h_\theta] \end{bmatrix} = \omega[\theta] + \begin{bmatrix} \mathbf{x}'' \cdot \mathbf{d}_1[\theta] \\ \mathbf{x}'' \cdot \mathbf{d}_2[\theta] \end{bmatrix} \cdot h_\theta + o(h_\theta) $$

Which finally leads to the derivative of the material curvatures vector, denoting $R_{\pi/2}$ the counter-clockwise rotation matrix of angle $\frac{\pi}{2}$:

$$ \mathbf{D}_{\theta}\omega(\theta).h_{\theta} = -R_{\pi/2}\omega[\theta].h_{\theta} \qquad , \qquad R_{\pi/2} = \begin{bmatrix} 0 & -1 \\ 1 & 0 \end{bmatrix} $$

### 3.6.5. Moment of torsion

The moment of torsion is the functional derivative of the potential energy with respect to $\theta$:

$$ \langle -m^{\text{int}}(s); h_{\theta} \rangle = \int \frac{\delta E_p}{\delta \theta} (s) h_{\theta} = \mathbf{D}_{\theta} E_p(\theta) \cdot h_{\theta} = \mathbf{D}_{\theta} E_b(\omega(\theta)) \cdot h_{\theta} + \mathbf{D}_{\theta} E_t(\theta) \cdot h_{\theta} = \mathbf{D}_{\theta} E_t(\theta) \cdot h_{\theta} $$

Decomposing the previous calculus gives:

$$ \mathbf{D}_{\theta} E_t(\theta) \cdot h_{\theta} = \left. \frac{d}{d\lambda} E_t(\theta + \lambda h_{\theta}) \right|_{\lambda=0} = \int_0^L (\beta(\theta' - \bar{\theta}')(\delta_L - \delta_0) - (\beta(\theta' - \bar{\theta}'))') \cdot h_{\theta} $$

$$ \mathbf{D}_{\omega} E_b(\omega) \cdot \mathbf{h}_{\omega} = \left. \frac{d}{d\lambda} E_b(\omega + \lambda \mathbf{h}_{\omega}) \right|_{\lambda=0} = \int_0^L (\omega - \bar{\omega})^{\text{T}} \mathbf{B} \cdot \mathbf{h}_{\omega} $$

$$ \mathbf{D}_{\theta} E_b(\omega(\theta)) \cdot h_{\theta} = \mathbf{D}_{\omega} E_b(\omega(\theta)) \cdot (\mathbf{D}_{\theta} \omega(\theta) \cdot h_{\theta}) = \int_{0}^{L} -(\omega - \bar{\omega})^{\text{T}} \mathbf{B} R_{\pi/2} \omega \cdot h_{\theta} $$

Finally, the moment of torsion is given by:

$$ m^{\text{int}}(s) = (\beta(\theta' - \bar{\theta}'))' + (\omega - \bar{\omega})^{\text{T}} \mathbf{B} R_{\pi/2} \omega - \beta(\theta' - \bar{\theta}')(\delta_L - \delta_0) \qquad , \quad m^{\text{int}}(s) = m^{\text{int}}(s) \mathbf{d}_3 $$

### 3.6.6. Forces

The internal force acting on the centerline is the functional derivative of the potential energy with respect to $\mathbf{x}$:
$$ \begin{aligned} \langle -\mathbf{f}^{\text{int}}(s); \mathbf{h}_x \rangle &= \int \frac{\delta E_p}{\delta \mathbf{x}}(s) \mathbf{h}_x = D_x E_p(\mathbf{x}). \mathbf{h}_x = D_x E_b(\omega(\mathbf{x})). \mathbf{h}_x + D_x E_t(\theta(\mathbf{x})). \mathbf{h}_x \\ &= D_\omega E_b(\omega(\mathbf{x})). (D_x \omega(\mathbf{x}). \mathbf{h}_x) + D_\theta E_t(\theta(\mathbf{x})). (D_x \theta(\mathbf{x}). \mathbf{h}_x) \\ &= D_\omega E_b(\omega(\mathbf{x})). (D_x \omega(\mathbf{x}). \mathbf{h}_x) \end{aligned} $$

Recall that with the quasi-static assumption, $D_x \theta(\mathbf{x}) = 0$ in the calculus of $\mathbf{f}^{\text{int}}$. By several integration by parts, using Fubini's theorem once and supposing that the terms in 0 and L vanish:

$$ \begin{aligned} D_x E_b(\omega(\mathbf{x})). \mathbf{h}_x &= D_\omega E_b(\omega(\mathbf{x})). (D_x \omega(\mathbf{x}). \mathbf{h}_x) \\ &= \int_0^L ((\omega - \bar{\omega})^\mathrm{T} B \cdot \begin{bmatrix} -\mathbf{d}_2^\mathrm{T} \\ \mathbf{d}_1^\mathrm{T} \end{bmatrix}) \cdot \mathbf{h}_x'' + \left( \int_0^L \frac{\delta \theta}{\delta \mathbf{x}}(s) \mathbf{h}_x dt \right) R \omega[\theta] \mathrm{d}s \\ &= \int_0^L \left( \begin{aligned}[t] & (\omega - \bar{\omega})^\mathrm{T} B \begin{bmatrix} -\mathbf{d}_2^\mathrm{T} \\ \mathbf{d}_1^\mathrm{T} \end{bmatrix} \end{aligned} \right)'' \cdot \mathbf{h}_x \\ &\quad + (\omega - \bar{\omega})^\mathrm{T} B \left( \int_0^L (\kappa \mathbf{b}(\delta_s - \delta_0) - (1 - H_s)\kappa \mathbf{b}') \cdot \mathbf{h}_x dt \right) R \omega \mathrm{d}s \\ &= \int_0^L \left( \begin{aligned}[t] & (\omega - \bar{\omega})^\mathrm{T} B \begin{bmatrix} -\mathbf{d}_2^\mathrm{T} \\ \mathbf{d}_1^\mathrm{T} \end{bmatrix} \end{aligned} \right)'' + (\omega - \bar{\omega})^\mathrm{T} B R \omega \kappa \mathbf{b} - \left( \int_s^L (\omega - \bar{\omega})^\mathrm{T} B R \omega \right) \kappa \mathbf{b}' \cdot \mathbf{h}_x \mathrm{d}s \\ &= \int_0^L \left( \begin{aligned}[t] & (\omega - \bar{\omega})^\mathrm{T} B \begin{bmatrix} -\mathbf{d}_2^\mathrm{T} \\ \mathbf{d}_1^\mathrm{T} \end{bmatrix} \end{aligned} \right)'' - \left( \left( \int_s^L (\omega - \bar{\omega})^\mathrm{T} B R \omega \right) \kappa \mathbf{b} \right)' \cdot \mathbf{h}_x \mathrm{d}s \end{aligned} $$

Using again the quasi-static hypothesis, this expression can be rewritten to give:

$$ \mathbf{f}^{\text{int}}(s) = - \left( (\omega - \bar{\omega})^{\mathrm{T}} B \begin{bmatrix} -\mathbf{d}_2^{\mathrm{T}} \\ \mathbf{d}_1^{\mathrm{T}} \end{bmatrix} \right)'' + (\beta(\theta' - \bar{\theta}') \kappa \mathbf{b})' $$

This expression is only valid for the open interval $]0, L[$. The boundary conditions will be treated directly within the discretized model, by conditions on $\kappa \mathbf{b}$.

## 4. Numerical model

The discrete beam will be described with n nodes $\mathbf{x}_0, \dots, \mathbf{x}_{n-1}$ and n cross-section orientations $\theta_0, \dots, \theta_{n-1}$ at these nodes (Figure 6). The segments will be denoted by $\mathbf{e}_i = \mathbf{x}_{i+1} - \mathbf{x}_i$, and their lengths $\|\mathbf{e}_i\|$.

### 4.1. Discrete curvature

The curvature binormal $\kappa \mathbf{b}_i$ at node $\mathbf{x}_i$ is defined by considering that the curvature $\kappa_i = \|\kappa \mathbf{b}_i\|$ is equal to the inverse of the radius of the circle passing through $\mathbf{x}_{i-1}, \mathbf{x}_i$ and $\mathbf{x}_{i+1}$ (Figure 7):

$$ \kappa \mathbf{b}_i = \frac{2}{\|\mathbf{e}_i\|\|\mathbf{e}_{i-1}\|\|\mathbf{e}_{i+1}\|} \frac{\mathbf{e}_{i-1} \times \mathbf{e}_i}{\|\mathbf{e}_{i-1} + \mathbf{e}_i\|} $$
<figure><img src="image_7.png" /><figcaption>Figure 6: Discretization of the beam</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: Discrete curvature and parallel transport</figcaption></figure>

## 4.2. Discrete Bishop frame and parallel transport

Let us denote by {$\mathbf{u}_0, \mathbf{v}_0$} the Bishop frame chosen in $\mathbf{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 {$\mathbf{u}_{i/i+1}, \mathbf{v}_{i/i+1}$} the constant Bishop frame on the segment $\mathbf{e}_i$, with $\mathbf{u}_{0/1} = \mathbf{u}_0$ and $\mathbf{v}_{0/1} = \mathbf{v}_0$, and by {$\mathbf{u}_i, \mathbf{v}_i$} the Bishop frames defined at the nodes $\mathbf{x}_i$. Denoting by $R_a$ the rotation matrix of axis $\mathbf{a}$ and angle $\|\mathbf{a}\|$, the parallel transport is defined as follows:

$$ 
\begin{aligned}
\mathbf{u}_i &= R_{\alpha_{i-1}^i, \mathbf{b}_i} \mathbf{u}_{i-1/i} \quad \mathbf{v}_i = R_{\alpha_{i-1}^i, \mathbf{b}_i} \mathbf{v}_{i-1/i} \\
\mathbf{u}_{i/i+1} &= R_{\alpha_i^i, \mathbf{b}_i} \mathbf{u}_i \quad \mathbf{v}_{i/i+1} = R_{\alpha_i^i, \mathbf{b}_i} \mathbf{v}_i \\
\alpha_{i-1}^i &= \sin^{-1}(\|\mathbf{e}_{i-1}\|\kappa_i/2) \\
\alpha_i^i &= \sin^{-1}(\|\mathbf{e}_i\|\kappa_i/2)
\end{aligned}
$$

These notations are illustrated in Figure 7. This definition of {$\mathbf{u}_i, \mathbf{v}_i$} is such that {$\mathbf{u}_i, \mathbf{v}_i$} is orthogonal to the circle passing through $\mathbf{x}_{i-1}, \mathbf{x}_i$ and $\mathbf{x}_{i+1}, \mathbf{v}_i$ and $\theta_i$ being defined, {$(\mathbf{d}_1)_i, (\mathbf{d}_2)_i$} is given by:

$$ 
\begin{aligned}
(\mathbf{d}_1)_i &= \cos \theta_i \mathbf{u}_i + \sin \theta_i \mathbf{v}_i \\
(\mathbf{d}_2)_i &= -\sin \theta_i \mathbf{u}_i + \cos \theta_i \mathbf{v}_i \\
\mathbf{\omega}_i &= \begin{bmatrix} \kappa \mathbf{b}_i \cdot (\mathbf{d}_1)_i \\ \kappa \mathbf{b}_i \cdot (\mathbf{d}_2)_i \end{bmatrix}
\end{aligned}
$$

## 4.3. Internal forces and moment of torsion

Finally, denoting the following terms by:

$$ \Delta\theta'_{i+\frac{1}{2}} = \frac{(\theta_{i+1} - \bar{\theta}_{i+1} - \theta_i + \bar{\theta}_i)}{\|\mathbf{e}_i\|} \qquad \mathbf{M}_i = -(\omega_i - \bar{\omega}_i)B_i \begin{bmatrix} -(\mathbf{d}_2)_i^T \\ (\mathbf{d}_1)_i^T \end{bmatrix} \qquad r_i = (\omega_i - \bar{\omega}_i)^T B_i R \omega_i $$

We can write the internal discretized forces and moments at node i:

$$ 
\begin{aligned}
\mathbf{f}_i^{\text{int}} &= \frac{1}{\|\mathbf{e}_{i-1}\|} \mathbf{M}_{i-1} - \left( \frac{1}{\|\mathbf{e}_i\|} + \frac{1}{\|\mathbf{e}_{i-1}\|} \right) \mathbf{M}_i + \frac{1}{\|\mathbf{e}_i\|} \mathbf{M}_{i+1} + \frac{\beta_i \kappa \mathbf{b}_i + \beta_{i+1} \kappa \mathbf{b}_{i+1}}{2} \Delta\theta'_{i+\frac{1}{2}} - \frac{\beta_{i-1} \kappa \mathbf{b}_{i-1} + \beta_i \kappa \mathbf{b}_i}{2} \Delta\theta'_{i-\frac{1}{2}} \\
\mathbf{m}_i^{\text{int}} &= \frac{\|\mathbf{e}_{i-1}\| + \|\mathbf{e}_i\|}{2} \mathbf{r}_i + \frac{\beta_i + \beta_{i+1}}{2} \Delta\theta'_{i+\frac{1}{2}} - \frac{\beta_{i-1} + \beta_i}{2} \Delta\theta'_{i-\frac{1}{2}}
\end{aligned}
$$

## 4.4. Dynamic Relaxation algorithm

Dynamic relaxation is a well-known explicit integration scheme [20], [21]. Here, a combined algorithm leads to the coupled static equilibrium regarding $\mathbf{x}$ and $\theta$:
<caption>**Algorithm 1** Coupled bending and twisting dynamic relaxation algorithm</caption>

1. Define the rest configuration $(\bar{\mathbf{x}}_i^0, \bar{\theta}_i^0)$
2. Define the initial configuration $(\mathbf{x}_i^0, \theta_i^0)$
3. Choose the convergence parameter $E_{c,\text{lim}}^x$
4. Nullify the translational speeds : $\dot{\mathbf{x}}_i^{t-\Delta t/2} = 0$
5. Nullify the angular speeds : $\dot{\theta}_i^{t-\Delta t/2} = 0$
6. Compute the forces $\mathbf{f}_i^t = \mathbf{f}_i^{\text{int},t} + \mathbf{f}_i^{\text{pen},t} + \mathbf{f}_i^{\text{ext},t}$
7. Compute the acceleration $\mathbf{a}_i^t = \mathbf{f}_i^t/M_i$
8. Speed increment : $\dot{\mathbf{x}}_i^{t+\Delta t/2} = \dot{\mathbf{x}}_i^{t-\Delta t/2} + \Delta t \mathbf{a}_i^t$

Compute the translational kinetic energy : $E_c^{x,t+\Delta t/2} = \frac{1}{2} \sum_{i=0}^{N-1} M_i (\dot{\mathbf{x}}_i^{t+\Delta t/2})^2$

10a. **If** $E_c^{x,t+\Delta t/2} > E_c^{x,t-\Delta t/2}$ : no translational kinetic energy peak
    position increment : $\mathbf{x}_i^{t+\Delta t} = \mathbf{x}_i^t + \Delta t \dot{\mathbf{x}}_i^{t+\Delta t/2}$

10b. **Else if** $E_c^{x,t+\Delta t/2} > E_{c,\text{lim}}^x$ : translational kinetic energy peak but no convergence
    **goto 4**

10c. **Else** : **end**

11. Compute the moments $\mathbf{m}_i^t = \mathbf{m}_i^{\text{int},t} + \mathbf{m}_i^{\text{ext},t}$
12. Compute the angular acceleration $\alpha_i^t = \mathbf{m}_i^t/I_i$
13. Speed increment: $\dot{\theta}_i^{t+\Delta t/2} = \dot{\theta}_i^{t-\Delta t/2} + \Delta t \alpha_i^t$
14. Compute the rotational kinetic energy : $E_c^{\theta,t+\Delta t/2} = \frac{1}{2} \sum_{i=0}^{N-1} I_i (\dot{\theta}_i^{t+\Delta t/2})^2$

15a. **If** $E_c^{\theta,t+\Delta t/2} > E_c^{\theta,t-\Delta t/2}$ : no rotational kinetic energy peak
    angle increment: $\theta_i^{t+\Delta t} = \theta_i^t + \Delta t \dot{\theta}_i^{t+\Delta t/2}$
    **goto 6**

15b. **Else** **goto 5**

**4.5. Boundary conditions**

To take account for boundary conditions, fictitious extensions are considered to compute $\kappa \mathbf{b}$ at beam's extremities depending if it is clamped (continuity of $\kappa$) or free ($\kappa = 0$).

**4.6. Inextensibility**

The inextensibility hypothesis is enforced by a penalty energy, which leads to axial forces. Denoting by $||\bar{\mathbf{e}}_i||$ the length of $\mathbf{e}_i$ in the rest configuration and by $A_i$ its axial stiffness, the penalty force at node $\mathbf{x}_i$ is given by:
$$f_i^{\text{pen}} = A_i \left( \frac{\|e_i\|}{\|e_i\|} - 1 \right) \frac{e_i}{\|e_i\|} - A_{i-1} \left( \frac{\|e_{i-1}\|}{\|e_{i-1}\|} - 1 \right) \frac{e_{i-1}}{\|e_{i-1}\|}$$

The value of $A_i$ must be sufficient so the inextensibility is satisfied. In practice, true axial stiffness of members will give good results and have the advantage to give directly the true axial strength in the considered beams, otherwise the axial strength must be computed with the ratio $ES/A_i$.

### 4.7. Fictitious lumped mass

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 [22], in order to ensure stability the following CFL condition holds:

$$\frac{k\Delta t^2}{M} < 2$$

Where $k$ is the stiffness of the system and $M$ the mass. The general relationship $M = k\Delta t^2$, with $k$ equal to the sum of the stiffnesses of the segments converging at $x_i$, holds. Noting that the bending stiffness is negligible compared to the axial stiffness for the evaluation of $M_i$, we chose:

$$M_i = \Delta t^2 \left( \frac{A_{i-1}}{\|e_{i-1}\|} + \frac{A_i}{\|e_i\|} \right)$$

$$I_i = \Delta t^2 \beta_i \left( \frac{1}{\|e_{i-1}\|} + \frac{1}{\|e_i\|} \right)$$

## 5. Test case

The model has been implemented in a python plugin for the software Rhino & Grasshopper. In this last section a benchmark with ABAQUS (a standard commercial FEM solver) is presented. A test case evaluating the influence of anisotropic cross sections over the relaxed shape of an initially spherical grid is also presented.

### 5.1. Benchmark with ABAQUS

In this benchmark, an initially rectilinear beam with a rectangular cross section is clamped at its two ends. The distance between the supports is then reduced while at one end the cross section is rotated by an angle of $\pi/2$. The beam is discretized in 100 elements. The results show an excellent accuracy even with a very large bending/torsion coupling.

<figure><img src="image_9.png" /><figcaption>Figure 8: deformed shape from custom plugin</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9: deformed shape from ABAQUS</figcaption></figure>
## 5.2. Influence of the cross section anisotropy over the relaxed shape of a grid

The following test case compare two beam grids with the same initial geometry but with different cross sections: square versus rectangular. Both grids are constructed with the compass method on a sphere of radius 5m, with a step length of 4m, and orthogonal directrix. The young modulus is set to 25GPa and the shear modulus to 10GPa. In the case of a square cross section, the torsion is negligible and the relaxed shape is closed to the sphere. In the case of a rectangular cross section the torsional stiffness induce important local bending at connections. Thus, the relaxed shape is quite far from the initial sphere.

<figure><img src="image_11.png" /><figcaption>Figure 10: comparing relaxed grids with different cross sections (square vs. rectangular)</figcaption></figure>

## 6. Conclusion

The paper presents a novel Kirchhoff beam element that accounts for both bending and torsion behaviours. This element has a minimal number of DoF. This reduction from classical 6 DoF elements is achieved thanks to a proper curve framing called zero-twisting frame. This element is implemented in a fast and robust extended dynamic relaxation algorithm which allow to find the static equilibrium of interconnected beam networks in non-linear large scale deformations.

This new tool, inspired from recent advanced works in the field of hair modeling, opens up new possibilities regarding formal explorations and structural design of actively bent structures.

## References

[1] M. Bergou, M. Wardetzky, S. Robinson, B. Audoly, and E. Grinspun, "Discrete elastic rods," *ACM SIGGRAPH*, pp. 1–12, 2008.

[2] R. L. Bishop, "There is more than one way to frame a curve," *Math. Assoc. Am.*, 1975.

[3] F. Otto, E. Schauer, J. Hennicke, and T. Hasegawa, *IL10 Grid Shells*. 1974.

[4] C. Douthe, J.-F. Caron, and O. Baverel, "Gridshell structures in glass fibre reinforced polymers," *Constr. Build. Mater.*, vol. 24, no. 9, pp. 1580–1589, 2010.

[5] N. Kotelnikova-Weiler, "Optimisation mécanique et énergétique d'enveloppes en matériaux composites pour les bâtiments," Université Paris-Est, 2012.
[6] O. Baverel, J.-F. Caron, F. Tayeb, and L. du Peloux, “Gridshells in composite materials : construction of a 300 m2 forum for the solidays’ festival in Paris,” *Struct. Eng. Int.*, vol. 22, no. 3, pp. 408–414, 2012.

[7] L. du Peloux and F. Tayeb, “Construction of the Forum Solidays,” 2011. [Online]. Available: https://www.youtube.com/watch?v=24LLfcVIZWw.

[8] L. du Peloux, F. Tayeb, O. Baverel, and J.-F. Caron, “Faith can also move composite gridshells,” in *IASS symposium 2013*, 2013, pp. 1–7.

[9] L. du Peloux, “Construction of the Ephemeral Cathedral of Créteil,” 2013. [Online]. Available: http://vimeo.com/59726330.

[10] G. Foret and O. Limam, “Experimental and numerical analysis of RC two-way slabs strengthened with NSM CFRP rods,” *Constr. Build. Mater.*, vol. 22, no. 10, pp. 2025–2030, 2008.

[11] S. Adriaenssens, M. Barnes, and C. Williams, “A new analytic and numerical basis for the form-finding and analysis of spline and gridshell structures,” in *Computing Developments in Civil and Structural Engineering*, B. Kumar and B. H. V Topping, Eds. Edinburgh: Civil-Comp Press, 1999, pp. 83–91.

[12] M. Barnes, S. Adriaenssens, and M. Krupka, “A novel torsion/bending element for dynamic relaxation modeling,” *Comput. Struct.*, vol. 119, pp. 60–67, Apr. 2013.

[13] S. Adriaenssens, “Stressed spline structures,” University of Bath, 2000.

[14] J. Li and J. Knippers, “Rotation formulations for dynamic relaxation with application in 3D framed structures with large displacements and rotations,” in *International Symposium of the IASS-APCS*, 2012.

[15] B. D'Amico, A. Kermani, and H. Zhang, “Form finding and structural analysis of actively bent timber grid shells,” *Eng. Struct.*, vol. 81, pp. 195–207, 2014.

[16] B. Audoly, M. Amar, and Y. Pomeau, *Elasticity and geometry*. 2010.

[17] S. Nabei, “Mechanical form-finding of timber fabric structures,” EPFL, 2014.

[18] R. Abraham, J. E. Marsde, and T. Ratiu, *Manifolds, Tensor Analysis, and Applications (Ralph Abraham, Jerrold E. Marsden and Tudor Ratiu)*. 2002.

[19] R. de Vries, “Evaluating changes of writhe in computer simulations of supercoiled DNA,” *J. Chem. Phys.*, vol. 122, no. 6, 2005.

[20] A. Day, “An Introduction to dynamic relaxation,” *Eng.*, 1965.

[21] M. Papadrakakis, “A method for the automatic evaluation of the dynamic relaxation parameters,” 1981.

[22] W. Lewis, *Tension structures: form and behaviour*. Telford, Thomas, 2003.