# A numerical and experimental study of the quasi-static deployment of membrane tubes

Rabah Bouzidi, Sandrine Buytet, Anh Le van

LUNAM University, University of Nantes, GeM, Institute for Research in Civil and Mechanical Engineering, CNRS UMR 6183, 2 rue de la Houssiniere, BP 92208 44322 Nantes Cedex 3, France

The paper is concerned with the deployment of inflatable structures in quasi static conditions. The problem is formulated as a minimization problem under constraints with emphasis on the self contact aspect and it is solved by a descent algorithm of the first order which allows to overcome the singularity or ill conditioning difficulty specific to membrane structures. Numerical computations are carried out on a simply folded membrane tube and compared with experimental results in order to prove the ability of the proposed numerical approach to correctly predict the deployment process.

## 1. Introduction

Inflatable structures are membrane like structures which acquire their stiffness and bearing capacities when inflated by air or another gas. They possess highly interesting mechanical properties which account for their increasing development in space industry: compared to standard structures, they are light, easy to be folded, rapidly deployable, their storage volumes are small, their manufacturing is less expensive, and moreover they have better thermal properties. The successful inflatable antenna experiment (IAE) carried out by Freeland and Veal (1998) showed that inflatable structures can be used in space applications as an alternative to conventional electromechanical deployable systems. Many types of inflatable structures have then been developed especially in building large scale systems such as solar arrays, solar sails, radar antennas, thermal or light shields, and a lot of studies have been devoted to them. Yet not everything is easy with the deployment of inflatable structures. The main difficulty is to control the change of geometry all along the deployment stage and further experimental or numerical investigations are still necessary in order to better understand the deployment mechanism.

As far as the numerical computations are concerned, one has to take into account different types of nonlinearities due to the large

deformations and finite strains, the material nonlinearities, the inflating pressure which is a follower loading, the contact of the membranes with external obstacles and particularly the self contact. Bifurcation may also be involved as wrinkles and folds can appear in the membranes. Such nonlinear problems are usually dealt with by means of the standard finite element method based on the discretization of principle of virtual work, where the equation of motion is discretized and solved at each increment by an iterative Newton like scheme. The specific difficulty encountered is due to the fact that the in plane stiffness is much higher than the flexural stiffness which is quite negligible, and, as a consequence, the stiffness matrix is almost singular as long as the inflatable structure is not fully deployed and taut. This singularity can be avoided by considering the deployment in the dynamic context, as done by Salama et al. (1999), Clem et al. (2000) and Wang and Johnson (2002) using an explicit scheme. In the first computational steps when the stiffness matrix is singular, the inertial terms allow the global tangent matrix to be regularized. In general, the deployment scenario depends on the rate of the inflating pressure, as shown in Wang and Johnson (2002).

On the other hand, the dynamical solution may be difficult especially when the deployment takes place in so low speed that the inertial forces are not significant. In this case, one had better abandon the dynamic point of view and solve the deployment problem in the quasi static framework. Few works have been devoted to the self contact of the inflatable membranes in large deformations. The aim of this paper is to propose a numerical
approach to solve the quasi static deployment problem, which
cope with the above mentioned singularity difficulty in a simple
and efficient way. The self contact is assumed to be frictionless
so that all forces are conservative and the solution is obtained by
minimizing a convenient total potential energy. A descent algo
rithm of first order will be chosen to achieve the minimization,
as done in Bouzidi and Levan (2004) and Gil (2006). Although first
order methods have only a linear rate of convergence, they have a
great advantage over the Newton iterative scheme of requiring the
gradient of the energy only not the stiffness matrix thus allow
ing one to compute all the equilibrium configurations, even when
the inflatable structure is not fully stretched. The proposed ap
proach is particularly interesting in the case of weak inflating pres
sures, as those used in space applications.

The paper is outlined as follows. First, the basic equations of the contact problem will be set out and then it will be shown how to recast the resulting strong formulation into a minimization problem under constraints. Next, the algorithm of the projected gradient method used to solve the minimization problem and the finite implementation will be described. Numerical computations of the deployment will be carried out on a simply folded membrane tube in Section 6. Experimental results obtained on the same tube will be presented in Section 7, which prove the ability of the proposed numerical approach to correctly predict the deployment process.

## 2. Strong formulation of the self-contact problem

### 2.1. *The membrane structure*

Let us consider a membrane structure with no bending stiffness,
which undergoes large deformations so that it may come into self
contact over a part of its surface and during some time interval. We
shall describe the basic equations of the problem using a Lagrang
ian description. The reference position of the membrane is defined
by the reference middle surface S₀ and the reference thickness H.
The current position S of the middle surface at any time t in the
time interval [O, T] of interest is given by the deformation φ:

$$
\forall t \in [O, T], \quad \phi(\cdot, t) : S_0 \ni \mathbf{X} \to \mathbf{X} \quad \phi(\mathbf{X}, t) \quad (1)
$$

where **X** and **x** are the reference and current positions of any parti-
cle on the middle surface. The positions **X** and **x** are parameterized
by the curvilinear coordinates ξ (ξ¹, ξ²) via the following chain
(Fig. 1):

$$
\begin{align}
& \text{A region of } \mathbb{R}^2 \xrightarrow{\theta} S_0 \xrightarrow{\phi} S \notag \\
& \phantom{\text{A region of }} \xi (\xi^1, \xi^2) \mapsto \mathbf{X} \quad \theta(\xi) \mapsto \mathbf{X} \quad \phi(\mathbf{X}, t) \quad \kappa(\xi, t) \quad (\kappa \quad \phi \circ \theta) \tag{2}
\end{align}
$$

Assuming that all the mappings introduced are smooth enough,
the covariant local base on the reference middle surface is defined
as

$$
\mathbf{A}_{\alpha} \quad \mathbf{X}_{,\alpha} \quad \frac{\partial \theta}{\partial \xi^{\alpha}} \quad \mathbf{A}_{3} \quad \mathbf{N} \quad \frac{\mathbf{A}_{1} \times \mathbf{A}_{2}}{\|\mathbf{A}_{1} \times \mathbf{A}_{2}\|} \quad (3)
$$

<figure><img src="image_1.png" /><figcaption>Fig. 1. Parametrization of the membrane middle surface.</figcaption></figure>

where Greek indices are equal to 1 or 2. **A**_3 **N** is the unit normal
vector at point **X** ∈ S₀. Similarly, the covariant local base on the
current middle surface is defined as

$$
\mathbf{a}_{\alpha} \quad \mathbf{x}_{,\alpha} \quad \frac{\partial \mathbf{k}}{\partial \xi^{\alpha}} \quad \mathbf{n} \quad \mathbf{a}_{3} \quad \frac{\mathbf{a}_{1} \times \mathbf{a}_{2}}{\| \mathbf{a}_{1} \times \mathbf{a}_{2} \|} \quad (4)
$$

The covariant components of the Green strain tensor **E** on the
middle surface are

$$
E_{\alpha\beta} \frac{1}{2} (a_{\alpha\beta} A_{\beta\alpha}) \quad (5)
$$

where $A_{\alpha\beta}$ $A_\alpha A_\beta$ and $a_{\alpha\beta}$ $a_\alpha$ $a_\beta$ are the covariant components of the metric tensor in the reference and current configurations, respectively.

The membrane is assumed to be made of a hyperelastic material defined by a strain energy density w(E) per unit volume. The second Piola Kirchhoff stress tensor Σ is related to strain tensor E via the constitutive relationship

$$
\Sigma \frac{\partial w(\mathbf{E})}{\partial \mathbf{E}} \qquad (6)
$$

The plane stress assumption in the membranes is satisfied by zeroing the normal stress, $\Sigma^{33} 0$, in the energy density evaluation. The reference surface $S_0$ is partitioned into three non overlapping and pairwise disjoint parts denoted $S_{0U}$, $S_{0T}$ and $S_{0c}$, where $S_{0U}$ and $S_{0T}$ are the parts where the displacement $\mathbf{U}$ $\mathbf{x}$ $\mathbf{X}$ and the tractions are prescribed, respectively, and $S_{0c}$ is the part where self contact potentially takes place during the time interval $[O, T]$. The boundary conditions on $S_{0U}$ and $S_{0T}$ are standard:

$$
\begin{align}
\forall t \in [O, T], \quad \forall \mathbf{X} \in S_{0U}, \mathbf{U}(\mathbf{X}, t) \quad \bar{\mathbf{U}}(\mathbf{X}, t), \notag \\
\forall \mathbf{X} \in S_{0T}, \mathbf{\Pi}(\mathbf{X}, t) \cdot \mathbf{N}(\mathbf{X}) \quad \bar{\mathbf{T}}(\mathbf{X}, t) \tag{7}
\end{align}
$$

where $\bar{\mathbf{U}}$, $\bar{\mathbf{T}}$ are the prescribed displacement and nominal traction, respectively; $\mathbf{\Pi}$ is the first Piola Kirchhoff (nominal) stress tensor. The boundary conditions on the contact surface are deferred to Section 2.4 where the contact law is specified.

2.2. Contact kinematics

Contact, or more precisely self contact, is one of the chief agents
of the deployment and has to be taken into account in the analysis.
It is assumed that one portion S<sub>0c</sub><sup>(1)</sup> of the reference membrane sur
face can potentially come into contact with another portion S<sub>0c</sub><sup>(2)</sup> (so
that the above mentioned contact surface S<sub>0c</sub> is S<sub>0c</sub><sup>(1)</sup> &cup; S<sub>0c</sub><sup>(2)</sup>). In order
to describe the relative motion of surfaces S<sub>0c</sub><sup>(1)</sup> and S<sub>0c</sub><sup>(2)</sup>, we arbi
trarily choose one, say S<sub>0c</sub><sup>(2)</sup>, as the reference (the target or master)
and evaluate the proximity of the other, thus S<sub>0c</sub><sup>(1)</sup> (the contactor

<figure><img src="image_2.png" /><figcaption>Fig. 2. Parametrizations of the contact surfaces.</figcaption></figure>
or slave), with respect to it. In the sequel, the spatial counterpart of $S_{0c}^{(i)}$, $i \in \{1, 2\}$, will be denoted $S_c^{(i)}$, e.g., $S_c^{(i)} \phi(S_{0c}^{(i)}, t)$.

For contact purposes, let us define the following parametrizations for the reference surface $S_{0c}^{(i)}$ and the current surface $S_c^{(i)}$ (Fig. 2):

$$ \text{A region of } \mathbb{R}^2 \to S_{0c}^{(i)} \to S_c^{(i)} \\ \xi^{(i)} \quad (\xi^{(i)1}, \xi^{(i)2}) \mapsto \mathbf{X}^{(i)} \quad \theta^{(i)}(\xi^{(i)}) \mapsto \mathbf{X}^{(i)} \quad \phi(\mathbf{X}^{(i)}, t) \quad \mathbf{k}^{(i)}(\xi^{(i)}, t) \quad (8) $$

The mappings $\theta^{(i)}$ and $\mathbf{k}^{(i)}(\cdot, t)$ $\phi(\cdot, t) \circ \theta^{(i)}$ are specific to the contact surfaces, they are not necessarily the same as the membrane parametrization (2). The target surface $S_c^{(2)}$ should be conveniently parametrized so that the normal at point $\mathbf{x}^{(2)}$ to $S_c^{(2)}$

$$ \mathbf{n}^{(2)} \frac{\partial \mathbf{k}^{(2)}}{\partial \xi^{(2)1}} \times \frac{\partial \mathbf{k}^{(2)}}{\partial \xi^{(2)2}} / \left\| \frac{\partial \mathbf{k}^{(2)}}{\partial \xi^{(2)1}} \times \frac{\partial \mathbf{k}^{(2)}}{\partial \xi^{(2)2}} \right\| \quad (9) $$

is directed toward the contactor $S_c^{(1)}$.

Now let $\mathbf{x}$ be a particular point on contactor surface $S_c^{(1)}$, one defines the opposite point $\mathbf{y} \in S_c^{(2)}$ as the closest point to $\mathbf{x}$ in the Euclidean sense (Fig. 3):

$$ \mathbf{y} \quad \arg \min_{\mathbf{x}^{(2)} \in S_c^{(2)}} \|\mathbf{x} \mathbf{x}^{(2)}\| \quad (10) $$

In writing (10), we assume that point $\mathbf{y}$ is unique without dis cussing the case of non uniqueness, which is not of major conse quence in numerical computations. The interested readers may consult comprehensive comments on this issue in the literature, e.g., Heegaard and Curnier (1993), Wriggers (2002, p. 48), Laursen (2002, p. 115), and references quoted therein. If the target surface $S_c^{(2)}$ is smooth at $\mathbf{y}$, then $\mathbf{y}$ is the projection of $\mathbf{x}$ on $S_c^{(2)}$. As may be seen, the notation $\mathbf{x}^{(i)}$, $i \in \{1, 2\}$, in (8) designates any point on $S_c^{(i)}$, whereas the notation $\mathbf{x}$ designates a particular point of interest on $S_c^{(1)}$ and $\mathbf{y}$ the opposite point of $\mathbf{x}$ on $S_c^{(2)}$ defined by (10).

The point $\mathbf{x} \in S_c^{(1)}$ being associated to (at least) one point $\mathbf{y} \in S_c^{(2)}$, the proximity $g$ is defined by

$$ g \quad \mathbf{v}(\mathbf{x}, \mathbf{y}) \quad (11) $$

where $\mathbf{v}$ is the normal at point $\mathbf{y}$. The sign convention chosen in (11) entails that the proximity $g$ is positive when there is interpenetration and negative when there is a gap between the two bodies. Thus, symbol $g$ is the *opposite* of what is referred to as the gap in the literature. In order to check that $g$ is an *objective* quantity, let us consider a change of reference frames defined by

$$ \mathbf{x}^{*} \quad \mathbf{Q}(t) \cdot \mathbf{x} + \mathbf{c}(t) \quad (12) $$

where vector $\mathbf{c}(t)$ represents the relative rigid body translation between the original reference frame and the new one (denoted

by superscript $*$) and $\mathbf{Q}(t)$ is an orthogonal tensor representing the relative rotation of the reference frames. From relations $\mathbf{x}^* \mathbf{y}^* \mathbf{Q}(t) (\mathbf{x} \mathbf{y})$ and $\mathbf{v}^* \mathbf{Q}(t) \mathbf{v}$, one derives

$$ \begin{aligned} g^* \quad & \mathbf{v}^* \cdot (\mathbf{x}^* \mathbf{y}^*) \quad (\mathbf{Q}(t)\mathbf{v}) \cdot (\mathbf{Q}(t)(\mathbf{x} \mathbf{y})) \\ & \mathbf{v} \cdot \mathbf{Q}^T(t)\mathbf{Q}(t) \cdot (\mathbf{x} \mathbf{y}) \quad g \end{aligned} \quad (13) $$

using equality $\mathbf{Q}^T(t)\mathbf{Q}(t) \mathbf{I}$ (identity tensor). Relation (13) shows that the proximity $g$ is objective.

Let us now convert all the spatial quantities defined above into material ones defined on the reference configuration, by means of the deformation $\phi$ in (1). Considering the point $\mathbf{X} \in S_{0c}^{(1)}$ related to the given point $\mathbf{x}$ by $\mathbf{x} \cdot \phi(\mathbf{X}, t)$, and the point $\mathbf{Y} \in S_{0c}^{(2)}$ related to point $\mathbf{y}$ by $\mathbf{y} \cdot \phi(\mathbf{Y}, t)$, the inverse image of $\mathbf{X}$ (resp. $\mathbf{Y}$) under the parametrization $\theta^{(1)}$ (resp. $\theta^{(2)}$) will be denoted $\xi$ (resp. $\eta$). Points $\xi$, $\mathbf{X}$ and $\mathbf{y}$ are linked through the following chain:

$$ \text{A region of } \mathbb{R}^2 \to S_{0c}^{(1)} \to S_c^{(1)} \\ \xi \mapsto \mathbf{X} \quad \theta^{(1)}(\xi) \mapsto \mathbf{X} \cdot \phi(\mathbf{X}, t) \quad \mathbf{k}^{(1)}(\xi, t) \quad (14) $$

whereas points $\eta$, $\mathbf{Y}$ and $\mathbf{y}$ are linked by

$$ \text{A region of } \mathbb{R}^2 \to S_{0c}^{(2)} \to S_c^{(2)} \\ \eta \mapsto \mathbf{Y} \quad \theta^{(2)}(\eta) \mapsto \mathbf{Y} \cdot \phi(\mathbf{Y}, t) \quad \mathbf{k}^{(2)}(\eta, t) \quad (15) $$

By making $\mathbf{x}^{(2)} \quad \mathbf{k}^{(2)}(\xi^{(2)}, t)$ in (10), one finds that $\eta$ is characterized by

$$ \eta \quad \arg\min_{\xi^{(2)}} \|\mathbf{x} \quad \mathbf{k}^{(2)}(\xi^{(2)}, t)\| \quad (16) $$

For notational conveniences, the dependencies on $\phi$ and $\mathbf{k}^{(2)}$ will be omitted in the material kinematics variables and only the dependency on $(\mathbf{X}, t)$ will be shown. Thus, for instance, one will write $\eta \quad \eta(\mathbf{X}, t)$ and $g \quad g(\mathbf{X}, t)$.

## 2.3. Contact tractions

We assume frictionless contact during the deployment process and resolve the contact Cauchy traction vector $\mathbf{t}(\mathbf{x}, t)$ at a typical point $\mathbf{x} \in S_c^{(1)}$ in terms of the normal vector $\mathbf{v}$ at the projection point $\mathbf{y} \in S_c^{(2)}$ as follows

$$ \mathbf{t} \quad t_N \mathbf{v} \quad (17) $$

The normal component $t_N$ of traction vector $\mathbf{t}$ in the direction of normal $\mathbf{v}$ is the contact pressure, positive when the contact occurs. Consider the change of reference frames defined by (12) again. Since the Cauchy traction vector $\mathbf{t}$ itself is objective, $\mathbf{t}^* \mathbf{Q}(t) \mathbf{t}$, one can recast the contact pressure in the new frame as

$$ t_N^* \quad \mathbf{t}^* \cdot \mathbf{v}^* \quad (\mathbf{Q}(t) \cdot \mathbf{t}) \cdot (\mathbf{Q}(t) \cdot \mathbf{v}) \quad \mathbf{t} \cdot \mathbf{Q}^T(t)\mathbf{Q}(t) \\ \cdot \mathbf{v} \cdot \mathbf{t} \cdot \mathbf{v} \quad t_N \quad (18) $$

which means that the contact pressure $t_N$ is an *objective* variable.

Let $\mathbf{T}(\mathbf{X}, t) \quad \mathbf{\Pi} \quad \mathbf{N}$ be the nominal Piola Kirchhoff traction vector at point $\mathbf{X} \in S_{0c}^{(1)}$ ($\mathbf{N}$ is the normal vector at $\mathbf{X}$ to $S_{0c}^{(1)}$). According to relation $\mathbf{t} \mathbf{d} S^{(1)} \quad \mathbf{T} \mathbf{d} S_o^{(1)}$ ($\mathbf{d} S_o^{(1)}$ is a differential reference area, $\mathbf{d} S^{(1)}$ its spatial counterpart), the nominal traction $\mathbf{T}$ can be resolved in terms of normal $\mathbf{v}$ at $\mathbf{y} \in S_c^{(2)}$ as follows

$$ \mathbf{T}(\mathbf{X}, t) \quad T_N(\mathbf{X}, t)\mathbf{v} \quad (19) $$

where the material and spatial normal traction components are related by $t_N \mathbf{d} S^{(1)} \quad T_N \mathbf{d} S_o^{(1)}$.

## 2.4. Contact law

The contact law describing the interfacial response at current time $t$ is naturally expressed with spatial variables. Here, it is recast

<figure><img src="image_3.png" /><figcaption>Fig. 3. A contactor point $\mathbf{x}$ and its projection $\mathbf{y}$ on the target surface. In comparison with Fig. 2, $\mathbf{x}$ is a particular point $\mathbf{x}^{(1)}$ and $\mathbf{y}$ a particular point $\mathbf{x}^{(2)}$.</figcaption></figure>
by means of relation $tdS^{(1)}$ $TdS_{o}^{(1)}$ into the mixed form involving the spatial proximity and the nominal contact tractions, which is more suitable for the numerical implementation within the Lagrangian description framework. The contact law in frictionless contact then reads

$$ \begin{cases} \forall t \in [0, T], \quad \forall \mathbf{X} \in S_{0c}^{(1)}, \quad g(\mathbf{X}, t) \le 0 \\ \cdot \text{ if } g(\mathbf{X}, t) < 0, \text{ then } T_N(\mathbf{X}, t) = 0 \\ \cdot \text{ if } g(\mathbf{X}, t) = 0, \text{ then } T_N(\mathbf{X}, t) \ge 0 \end{cases} \quad (20) $$

In this work, use will be made of the penalty method and the above contact law will be replaced by a smoother one

$$ T_N \quad \epsilon_N \langle g \rangle \quad (21) $$

where $\epsilon_N > 0$ is the penalty parameter and the Macaulay brackets $\langle g \rangle$ are defined as

$$ \langle g \rangle = \begin{cases} g & \text{if } g \ge 0 \\ 0 & \text{if } g < 0 \end{cases} \quad (22) $$

One can check that the contact laws (20) and (21) are objective by combining Relations (13) and (18).

### 3. The constrained minimization problem for the deployment

#### 3.1. Potentials

The internal potential energy of elastic strains $\Pi_{int}$ is defined from the strain energy density $w(E)$ via

$$ \Pi_{int}(\mathbf{U}) \quad \int_{S_0} H w(\mathbf{E}) dS_0 \quad (23) $$

The external loadings are of two types: (i) a follower force which is the internal pressure $p(t)$ inflating the membrane struc- ture during its deployment and (ii) a dead load such as the weight of the membrane itself. Since the pressure is applied on a closed surface, it is a conservative load like the dead load. Thus, all the external loadings give rise to the potential

$$ \Pi_{ext}(\mathbf{U}, t) \quad \int_{S_{0T}} \bar{\mathbf{T}} \cdot \mathbf{U} dS_0 + p(t) V(\mathbf{U}) \quad (24) $$

where $\bar{\mathbf{T}}$ is the dead load per unit surface prescribed over $S_{0T}$, $V(\mathbf{U})$ is the volume enclosed by the membrane surface, which depends on the current position of the membrane and can be expressed as a function of the displacement field $\mathbf{U}$ by

$$ V(\mathbf{U}) \quad \frac{1}{3} \int_S \mathbf{x} \cdot \mathbf{n} dS \quad \frac{1}{3} \int_S (\mathbf{X} + \mathbf{U}) \cdot \mathbf{n} dS \quad (25) $$

where $\mathbf{n}$ is the outward unit normal vector to the membrane sur- face, as shown in Fig. 1.

In fact, the membrane considered in this work is very thin so that its weight is almost negligible. One then has to add some other dead loads which will be specified later in the numerical and experimental parts in order to stabilize the quasi static deploy- ment process.

From the contact law (21), the potential of the contact forces $\Pi_c$ is a quadratic functional of $\langle g \rangle$:

$$ \Pi_c(\mathbf{U}) \quad \int_{S_{0c}^{(1)}} \frac{1}{2} \epsilon_N \langle g \rangle^2 dS_0 \quad (26) $$

#### 3.2. The constrained minimization problem

To determine the solution of the problem, one may use the well known principle of minimum total potential energy which states that the actual displacement field $\mathbf{U}$ minimizes the total

potential energy $\Pi(\mathbf{U}, t)$ $\Pi_{int}(\mathbf{U})$ $\Pi_{ext}(\mathbf{U}, t)$ + $\Pi_c(\mathbf{U})$ among all kinematically admissible displacements. The energy minimization method was successfully applied by Bonet et al. (2000) to air supported membrane structures and Gil (2006) to prestressed membranes undergoing large strains.

Here, the minimization problem will be transformed into an other one which proves to be more convenient for the deployment computations. Indeed, as we shall see in the experimental works (Section 7), the inflating pressure $p(t)$ is not a monotonous function during the deployment: it first increases up to a maximum value and then decreases as soon as the deployment is in a sufficiently advanced stage. There is not a one to one relationship between the internal pressure and the maximum displacement for instance. Meanwhile, the volume $V(\mathbf{U})$ enclosed by the membrane surface keeps on increasing till the end of the deployment.

As a consequence, one had better take volume $V(\mathbf{U})$ as the con- trol parameter, instead of pressure $p(t)$. To do this, let us rewrite the total potential energy as

$$ \Pi(\mathbf{U}) \quad \int_{S_0} H w(\mathbf{E}) dS_0 \quad \int_{S_{0T}} \bar{\mathbf{T}} \cdot \mathbf{U} dS_0 + \int_{S_{0c}^{(1)}} \frac{1}{2} \epsilon_N \langle g \rangle^2 dS_0 \quad (27) $$

where only the potential due to external dead loads is retained in $\Pi_{ext}$ and there is no more potential due to the pressure. Solving the problem then reduces to the following minimization problem of $\Pi(\mathbf{U})$ under constraint: given a prescribed volume $\bar{V}$,

$$ \begin{aligned} \text{find } \mathbf{U} \text{ minimizing } \Pi(\mathbf{U}) \quad & \int_{S_0} H w(\mathbf{E}) dS_0 \quad \int_{S_{0T}} \bar{\mathbf{T}} \cdot \mathbf{U} dS_0 \\ & + \int_{S_{0c}^{(1)}} \frac{1}{2} \epsilon_N \langle g \rangle^2 dS_0 \text{ subject to } h(\mathbf{U}) \\ & \frac{1}{3} \int_{S} \mathbf{x} \cdot \mathbf{n} dS \quad \bar{V} = 0 \end{aligned} \quad (28) $$

Besides the volume constraint $h(\mathbf{U})$ $V(\mathbf{U})$ $\bar{V}$ = 0, there exist other constraints to be taken into account as well, such as the boundary conditions and the symmetry conditions (only half the tube is modeled in finite element computations). In Relation (28), only the volume constraint is mentioned since it is less standard than others and deserves particular attention in the deployment problem. The other constraints are not mentioned explicitly for the sake of clarity.

#### 3.3. Algorithm of the projected gradient method

Minimization problems are usually solved by means of either first order methods (Gil, 2006) or second order methods (Bonet et al., 2000; Gil, 2006). Second order methods necessitate the tan- gent stiffness matrix which is often singular when the membrane is not yet taut. Although first order methods have only a linear rate of convergence, they prove to be particularly efficient for inflatable structures as they require computing the first gradient of the total potential energy only, not the Hessian matrix. Thus, it is possible to compute equilibrium configurations which exist when the membranes are not taut. Here, use will be made of the projected gradient method formulated by Goldstein (1964) and Levitin and Polyak (1966) in order to solve the minimization problem under volume control (28). This is a first order method of the same type as those used in Bouzidi and Levan (2004) and Gil (2006).

Let $\Pi(\mathbf{U})$ be the functional to be minimized under the con- straints $\mathbf{h}(\mathbf{U})$ ($h_1(\mathbf{U}), h_2(\mathbf{U}), \dots, h_n(\mathbf{U})$) = 0. The simple update relation $\mathbf{U}^{i+1} = \mathbf{U}^i + \alpha \mathbf{s}^i$, where $\mathbf{s}^i$ is the descent direction at itera- tion $i$ and $\alpha$ the descent step, needs to be corrected since it does not fulfill the constraints $\mathbf{h}(\mathbf{U}) = 0$. The key idea of the algorithm lies in using the projected point $\bar{\mathbf{U}}^i$ of $\mathbf{U}^i$ onto the constraints sub space $\mathbf{h}(\mathbf{U}) = 0$ and searching the solution increment $\mathbf{U}^{i+1} = \mathbf{U}^i$ in
<figure><img src="image_4.png" /><figcaption>Fig. 4. Projected gradient algorithm.</figcaption></figure>

the tangent hyperplane of $\mathbf{h}(\mathbf{U})$ **0**, see Fig. 4. Thus, the update relation writes

$$ \mathbf{U}^{i+1} \quad \mathbf{U}^i + \alpha \mathbf{s}^i + (\bar{\mathbf{U}}^i \quad \mathbf{U}^i) \quad (29) $$

where the descent direction $\mathbf{s}^i$ is chosen to be a unit vector tangent to the constraints subspace.

Therefore, problem (28) with unknown **U** is transformed into another one whose unknown is the descent direction **s**:

$$ \text{Find } \mathbf{s} \text{ minimizing } \mathbf{s}^T \nabla_{\mathbf{U}} \Pi \quad \text{subject to} \quad \mathbf{B}^T \mathbf{s} \quad \mathbf{0} \text{ and } \mathbf{s}^T \mathbf{s} \quad \mathbf{1} \quad (30) $$

where **B** **V**_U**h** is the gradient matrix of **h**. Here, the functional to minimize represents the first variation of energy Π defined in (27). The Lagrangian $\mathcal{L}$ associated to problem (30) is

$$ \mathcal{L}(\mathbf{s}, \lambda, \mu) \quad \mathbf{s}^T \nabla_{\mathbf{U}} \Pi \quad \mathbf{s}^T \mathbf{B} \lambda \quad \mu(\mathbf{s}^T \mathbf{s}) \quad \mathbf{1} \quad (31) $$

whence the stationarity condition of $\mathcal{L}$

$$ \frac{\partial \mathcal{L}}{\partial \mathbf{s}} \quad \mathbf{V}_U \Pi \quad \mathbf{B} \lambda \quad 2\mu \mathbf{s} \quad \mathbf{0} \quad (32) $$

Premultiplying Relation (32) by $\mathbf{B}^T$ gives the Lagrange multiplier $\lambda$:

$$ \lambda \quad [\mathbf{B}^T \mathbf{B}]^{-1} \mathbf{B}^T \nabla_{\mathbf{U}} \Pi \quad (33) $$

Re inserting $\lambda$ in (32) then leads to **s**:

$$ \mathbf{s} \quad \frac{1}{2\mu} \left[ \mathbf{I} \quad \mathbf{B} \quad [\mathbf{B}^T \mathbf{B}]^{-1} \mathbf{B}^T \right] \nabla_{\mathbf{U}} \Pi \quad (34) $$

where the multiplier $\mu$ has to be adjusted so that $||\mathbf{s}||$ 1. Further more, it can be shown from the relation $\mathbf{h}(\mathbf{U}^i) \quad \mathbf{B}^T(\mathbf{U}^i \quad \bar{\mathbf{U}}^i)$ that

$$ \bar{\mathbf{U}}^i \quad \mathbf{U}^i \quad [\mathbf{B}\mathbf{B}^T]^{-1} \mathbf{B}\mathbf{h}(\mathbf{U}^i) \quad (35) $$

Finally, the update scheme (29) becomes

$$ \mathbf{U}^{i+1} \quad \mathbf{U}^i + \frac{\alpha}{2\mu} \left[ \mathbf{I} \quad \mathbf{B} \quad [\mathbf{B}^T \mathbf{B}]^{-1} \mathbf{B}^T \right] \nabla_{\mathbf{U}} \Pi \quad [\mathbf{B}\mathbf{B}^T]^{-1} \mathbf{B}\mathbf{h}(\mathbf{U}^i) \quad (36) $$

The unknown $\alpha$ is obtained by minimizing the energy $\Pi(\mathbf{U}^{i+1}) \quad \Pi(\bar{\mathbf{U}}^i + \alpha \mathbf{s}^i)$ with respect to scalar variable $\alpha$ by the so called line search method. The iterative process is stopped when the displacement increment $||\mathbf{U}^{i+1} \quad \mathbf{U}^i||$ is small enough.

## 4. Finite element implementation

The projected gradient method is implemented by means of the finite element method. We shall focus on the discretization of the

contact potential (26) and the volume (25), as other terms in (28) can be dealt with in a standard way. In the sequel, the following notational convention is adopted for the matrix representations: curly brackets {} designate a column vector and the square brackets [] a general matrix.

Let the contactor surface $S_{0c}^{(1)}$ be subdivided into a number of isoparametric finite elements $e^{(1)}$ (the superscript '(1)' is used throughout to remind of the contactor surface $S_{0c}^{(1)}$), then the contact potential (26) can be recast as the sum of the element terms related to elements $e^{(1)}$ and evaluated using the Gauss quadrature rule:

$$ \Pi_c \sum_{e^{(1)}} \frac{\epsilon_N}{2} \int_{e_\xi} \langle g \rangle^2 j(\xi) d\xi^1 d\xi^2 \quad \sum_{e^{(1)}} \frac{\epsilon_N}{2} \sum_{\text{Gauss points}} \langle g \rangle^2 j(\xi) w(\xi) \quad (37) $$

where $e_\xi$ is the reference element of $e^{(1)}, \xi$ ($\xi_1, \xi_2) \in \mathbb{R}^2$ the parent coordinates of point $\mathbf{X} \in e^{(1)}$ (see Fig. 3) and $j(\xi) \quad ||\mathbf{X}_{e_\xi}^T(\xi) \times \mathbf{X}_{e_\xi}^T(\xi)||$ is the area Jacobian. In the last sum in Relation (37) all terms are computed at Gauss point $\xi$ and $w(\xi)$ is the associated weight.

In any element $e^{(1)} \subset S_{0c}^{(1)}$, the reference and current coordinates are interpolated by

$$ \{\mathbf{X}\} \quad [\mathcal{N}^{(1)}(\xi)]\{\mathbf{X}\}e^{(1)} \quad \{\mathbf{X}\} \quad [\mathcal{N}^{(1)}(\xi)]\{\mathbf{X}\}e^{(1)} \quad (38) $$

In (38)$_b$ for instance, $\{\mathbf{X}\}$ is the column vector containing the three coordinates of point $\mathbf{x}$ in a fixed Cartesian orthonormal basis ($\mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3$). The column vector $\{\mathbf{x}\}^{e^{(1)}}$ contains the nodal coordinates of element $e^{(1)}$, and the matrix $[\mathcal{N}^{(1)}(\xi)]$ contains the shape functions of element $e^{(1)}$.

In a similar way, the target surface $S_{0c}^{(2)}$ is subdivided into a number of isoparametric finite elements. Given a quadrature point $\mathbf{X}$ in element $e^{(1)}$, i.e. $\mathbf{x} \quad \phi(\mathbf{X}, t)$, let $e^{(2)}$ the (or an) element of $S_{0c}^{(2)}$ containing the reference point $\mathbf{Y}$ corresponding to the projection $\mathbf{y} \quad \text{proj}_{S_{0c}^{(2)}}(\mathbf{x}) \in S_{0c}^{(2)}$. In element $e^{(2)}$, the reference and current coordinates are interpolated by

$$ \{\mathbf{Y}\} \quad [\mathcal{N}^{(2)}(\eta)]\{\mathbf{Y}\}e^{(2)} \quad \{\mathbf{Y}\} \quad [\mathcal{N}^{(2)}(\eta)]\{\mathbf{Y}\}e^{(2)} \quad (39) $$

where $\eta$ ($\eta_1, \eta_2) \in \mathbb{R}^2$ is the parent coordinate of $\mathbf{Y}$ and the super script '(2)' is used to remind of element $e^{(2)} \subset S_{0c}^{(2)}$. The column vector $\{\mathbf{y}\}^{e^{(2)}}$ contains the nodal coordinates of element $e^{(2)}$ and the matrix $[\mathcal{N}^{(2)}(\eta)]$ contains the shape functions of element $e^{(2)}$.

Given a quadrature point $\mathbf{X}$ in element $e^{(1)} \subset S_{0c}^{(1)}$, there is at least one element $e^{(2)} \subset S_{0c}^{(2)}$ containing the reference point $\mathbf{Y}$ of $\mathbf{y} \quad \text{proj}_{S_{0c}^{(2)}}(\mathbf{x})$. Both element $e^{(2)}$ and its node number, which depend on the projection $\mathbf{y}$, vary as functions of the quadrature point $\mathbf{X}$. Moreover, one element $e^{(1)}$ can be in contact with several elements $e^{(2)}$; consequently, the whole set of quadrature points of a given element $e^{(1)}$ in (37) generally involves several elements $e^{(2)}$ in $S_{0c}^{(2)}$. It results that the approach presented herein is of the node to segment rather than the segment to segment type (for the mentioned types, see, e.g., Kikuchi and Song, 1981; Harnau et al., 2005 and Puso and Laursen, 2004).

By using Relations (38) and (39), the potential of the contact forces (26) is computed from

$$ \Pi_c \sum_{e^{(1)}} \frac{\epsilon_N}{2} \sum_{\text{Gauss points}} \left\langle \left( [\mathcal{N}^{(2)}(\eta)]\{\mathbf{y}\}e^{(2)} - [\mathcal{N}^{(1)}(\xi)]\{\mathbf{x}\}e^{(1)} \right) \mathbf{v} \right\rangle^2 j(\xi)w(\xi) \quad (40) $$

* The reference middle surface $S_0$ is discretized into isoparametric elements $e$. The volume (25) of the structure is

$$ V(\mathbf{U}) \quad \frac{1}{3} \sum_e \int_{e_\xi} \mathbf{x} \cdot (\mathbf{a}_1 \times \mathbf{a}_2) d\xi^1 d\xi^2 \quad (41) $$
where $e_ζ$ is the reference element of $e$, the current position $\mathbf{x}$ is interpolated as in (38) by $\{\mathbf{x}\}$ [$N(\xi)\{\mathbf{x}\}^e$ and $\mathbf{a}_\alpha$ [$N_\alpha(\xi)\{\mathbf{x}\}^e$.

## 5. Contact detection

The main task in evaluating the contact potential (40) and its gradient is to determine the projection $\mathbf{y}$ of an integration point $\mathbf{x}$ onto the target surface $S_c^{(2)}$. The contact detection is a nonlinear procedure which may be highly CPU time consuming and requires an iterative process as explained in Benson and Hallquist (1990). In order to reduce the computational cost, most algorithms involve a preliminary vicinity test, which consists in identifying the closest elements likely to interpenetrate so that the contact computation can be restricted to those elements only, for more details see Kane et al. (1999).

It should be mentioned here another difficulty specific to the deployment of membrane structures. Usually, the contactor sur face $S_c^{(1)}$ may come into contact with the target surface $S_c^{(2)}$ from one side of the target body, namely the outward side as defined by the outward normal vector $\mathbf{n}^{(2)}$ to $S_c^{(2)}$, see the convention about normal $\mathbf{n}^{(2)}$ in definition (9). This situation is shown in Fig. 5a for the points $\mathbf{x}, \mathbf{y}$ and the normal $\mathbf{v}$ in Fig. 3. On the other hand, when the contact involves hollow bodies as is the case with (closed) membranes, the target surface $S_c^{(2)}$ may encounter the contactor surface $S_c^{(1)}$ from either side (Wang and Nakamachi, 1997). As can be seen in Fig. 5b, point $\mathbf{y}$ may potentially come into contact with point $\mathbf{x}$ or $\mathbf{x}'$, located in both sides of $S_c^{(2)}$ with respect to normal $\mathbf{v}$. Point $\mathbf{x}$ is in the same situation as in Fig. 5a and is treated as usual. Point $\mathbf{x}'$ needs a special operation: in order to prevent it from cross ing point $\mathbf{y}$, one has to replace normal $\mathbf{v}$ by the opposite normal $\mathbf{v}$.

### 5.1. Validation of the contact algorithm

The contact algorithm is tested on the following simple example which admits of an analytical expression for the contact sur face. Consider two identical spherical balloons of initial radius $R$ inflated by an internal pressure $p$ and compressed by two opposite forces of same magnitude $F$ passing through the centers of the spheres, as shown in Fig. 6a. The problem is symmetrical with respect to the mid plane perpendicular to forces $F$ and the contact surface is a planar disk. The balloons are modeled as membranes (thus having no bending stiffness) and the contact is assumed to be frictionless. By writing the equilibrium balance for one of the two balloons with its air and without the surface of contact, see Fig. 6b, one finds that the contact area $S_c$ is simply the ratio of the compressive load to the internal pressure: $S_c F/p$, thus inde pendent of the geometry and the constitutive material.

The numerical computations are done with the radius $R$ 50 cm, the membrane thickness $H$ 50 µm and by assuming

<figure><img src="image_5.png" /><figcaption>Fig. 6. Contact between two spherical balloons.</figcaption></figure>

the Saint Venant Kirchhoff potential $w(E)$ $\mu E : E + \frac{\lambda}{2}(trE)^2$, so that it comes from (6):

$$ \Sigma = 2\mu E + \lambda \mathrm{tr}(E)I \iff E = \frac{1+\nu}{E} \Sigma = \frac{\nu}{E} \mathrm{tr}(\Sigma)I \quad (42) $$

where $I$ is the unit second order tensor, $\lambda$ and $\mu$ are the Lamé con stants, $E$ is the Young's modulus and $\nu$ the Poisson's ratio. Here, we take $E = 5 \cdot 10^7$ Pa, $\nu = 0.3$. For symmetry reasons, only a quarter of the system is discretized, the mesh is made of three node triangle elements as shown in Fig. 7. The penalty parameter is $\epsilon_N = 10^6$ N/m$^3$ and the convergence criterion in the iterative loop (36) is set to $10^{-8}$. Numerical tests show that the penalty parameter $\epsilon_N$ may be given any value within a large range, say $10^5$ N/m$^3$ to $10^8$ N/m$^3$, without much changing the final results. However, if $\epsilon_N$ is greater than $10^6$ N/m$^3$, the iteration number required for conver gence rises notably. As a consequence, the penalty parameter should be chosen as a compromise between the accuracy of the re sults and the speed of convergence.

Two series of finite element computations are done: (i) the first with constant internal pressure $p$ 1000 Pa and force $F$ ranging from 10 N to 70 N, (ii) the second with constant force $F$ 10 N

<figure><img src="image_6.png" /><figcaption>Fig. 5. (a) Usual contact condition; (b) contact from either sides in deployment.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 7. Contact between two spherical balloons. Initial and deformed meshes.</figcaption></figure>
<table><caption>**Table 1**
Analytical and numerical contact areas versus the applied load (pressure $p$ 1000 Pa).</caption><thead><tr><th rowspan="2">Force<br>$F$ (N)</th><th colspan="3">Contact area $S_c(m^2)$</th><th rowspan="2">Top<br>displacement<br>$u$ (m)</th></tr><tr><th>Analytical<br>values</th><th>Numerical<br>values</th><th>Difference<br>(%)</th></tr></thead><tbody><tr><td>20</td><td>0.020</td><td>0.0203</td><td>1.5</td><td>0.1041</td></tr><tr><td>30</td><td>0.030</td><td>0.0310</td><td>3.6</td><td>0.1753</td></tr><tr><td>60</td><td>0.060</td><td>0.0613</td><td>2.2</td><td>0.3814</td></tr><tr><td>70</td><td>0.070</td><td>0.0693</td><td>1.0</td><td>0.4571</td></tr></tbody></table>

<table><caption>**Table 2**
Analytical and numerical contact areas versus the internal pressure (force $F$ 10 N).</caption><thead><tr><th rowspan="2">Pressure $p$ (Pa)</th><th colspan="3">Contact area $S_c(m^2)$</th><th rowspan="2">Top displacement $u$ (m)</th></tr><tr><th>Analytical values</th><th>Numerical values</th><th>Difference (%)</th></tr></thead><tbody><tr><td>350</td><td>0.0286</td><td>0.0290</td><td>1.6</td><td>0.171</td></tr><tr><td>400</td><td>0.0250</td><td>0.0256</td><td>2.5</td><td>0.151</td></tr><tr><td>450</td><td>0.0220</td><td>0.0224</td><td>2.0</td><td>0.134</td></tr><tr><td>500</td><td>0.0200</td><td>0.0202</td><td>1.3</td><td>0.120</td></tr><tr><td>550</td><td>0.0182</td><td>0.0182</td><td>0.1</td><td>0.107</td></tr><tr><td>650</td><td>0.0154</td><td>0.0158</td><td>2.7</td><td>0.085</td></tr><tr><td>750</td><td>0.0133</td><td>0.0137</td><td>2.7</td><td>0.066</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Fig. 8.** Deployment of a folded membrane tube.</figcaption></figure>

and internal pressure $p$ increasing from 350 Pa to 850 Pa. At each computational step, we search the whole surface of the balloons and as soon as an element is found to be in contact, its entire area is added to the contact area. Thus, the numerical value of the contact area can only vary in a discontinuous way, however fine the mesh is and however small the computational step is. It is difficult to precisely evaluate the contact area and the mesh chosen here is fine enough for the discontinuity to be acceptable. The results are displayed in Tables 1 and 2 for the numerical contact area in comparison with the theoretical values, together with the axial displacement $u$ of the point of application of force $F$. In the second calculation, the contact area decreases as pressure $p$ increases, it is the volume of the two balloons that increases with $p$.

Because of the above mentioned discontinuity, the error in creases when few finite elements are in contact, i.e. when the load is small or the pressure is high. The numerical values are in good agreement with the analytical ones: the maximum error is less than 4% in the first calculation and 3% in the second one. This example shows that the numerical developments for deployment by contact and self contact work well.

## 6. Numerical computations of the deployment

The numerical computations of the quasi static deployment are performed on a folded cylindrical membrane tube (length = 80 cm,

diameter = 15 cm or 20 cm, thickness = 50 µm), made of a hyperelastic St Venant Kirchhoff material with Young's modulus $E$ 5.10⁹ Pa and Poisson's ratio $\nu$ 0.3.

Initially, the tube is completely flattened and then folded in two equal parts as shown in Fig. 8. The bottom edge is fixed while the top edge is constrained to move freely along the $z$ axis perpendicular to the flattened position of the tube. At each step of the deployment, the volume of the tube is set to a prescribed value as stated in (28), which amounts to inflate the tube with a certain internal pressure. The computations are carried out in quasi statics and a dead load $Fz$ is uniformly distributed on the top edge in order to stabilize the deployment process. Here, the dead load is due to a mass of 691.7 g, 1383.9 g, 2050.7 g or 2742.9 g, i.e. $F$ 6.78 N, 13.62 N, 20.11 N or 26.9 N.

Because of the symmetry of the problem with respect to the middle vertical plane, only half of the tube is modeled. A very small initial gap of 0.01 cm is prescribed between the top and the bottom edges (equivalent to an opening angle of about 0.01°), so that the contact surfaces are very close yet not in contact at the first computational step. Fig. 8 shows the mesh used containing 529 nodes and 1024 three node triangular elements. The penalty parameter for the contact is $\epsilon_N$ 10⁶ N/m³ in all the numerical computations.

As noticed by Salama et al. (1999) in their study of the dynamic deployment of inflatable structures, high values of the Young's modulus may cause a numerical instability. This kind of difficulty

<figure><img src="image_9.png" /><figcaption>**Fig. 9.** Deployment configurations of the 20cm-diameter tube with applied dead load $F$ 13.62 N.</figcaption></figure>
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Internal pressure versus top displacement.</figcaption></figure>

is commonly observed in the deployment computations. In the present work, one finds a similar numerical locking phenomenon, which can be accounted for by two main reasons: (i) the facettization of the membrane surface due to the finite element mesh, and (ii) the significant difference between the in plane and out of plane stiffnesses inherent to the membrane model. With a high value of the Young's modulus *E*, the strain energy is found to be prevalent compared to other energies and generates a numerical noise, leading to a premature convergence towards local minima of the total potential energy. This noise may decrease with the mesh refinement but it remains a major hindrance to the

convergence. To remedy the situation, a simple solution is to artificially and temporarily diminish the membrane stiffness of the tube, by reducing *E* in the first iterations of each computational step and resetting it to its real value in the final iterations of the step. The choice of the artificial value for *E* has little influence on the deployment kinematics since the deployment involves large displacements and rotations, rather than large strains. For our example considered herein, we set $E \ 5.10^7$ Pa in the first iterations a small value conveniently chosen so as to avoid the above mentioned numerical difficulties and yet maintain small strains and reset the Young's modulus to its real value $E \ 5.10^9$ Pa in the final iterations of the step.

Fig. 9 shows the 20 cm diameter tube at different stages of deployment with applied dead load *F* 13.62 N. The numerical results show that the deployment is possible only if the diameter of the tube is large enough or load *F* is small enough: for the heaviest load *F* 26.9 N (2742.9 g), convergence is achieved for the 20 cm diameter tube, not for the smaller tube. The inflating pressure *p* is

<figure><img src="image_12.png" /><figcaption>**Fig. 12.** Deployment stages of the Z-folded tube with applied dead load *F* 1 N.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Finite element model of a Z-folded tube.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Fig. 13. Experimental schema.</figcaption></figure>

plotted versus the top displacement along direction **z** in Fig. 10 for the prescribed *F* values. As can be seen, for a given dead load *F*, the pressure reaches a maximum value at about halfway and then decreases until the end of the deployment. The higher the load *F*, the higher the maximum pressure is. On the other hand, the volume of the tube is monotonically increasing throughout the deployment process.

### 6.1. Deployment of a Z folded tube

Another example is briefly described here to show how the numerical computations work on more complex geometries with multiple self contact surfaces. We consider the same tube as in the previous section and assume now that it is folded in four equal parts. As before, only half of the tube with respect to the middle vertical plane is modeled. The boundary conditions are prescribed so that the top and medium edges of the tube slide freely along the vertical direction and a dead load force *F* 1 N is applied on the top edge. The finite element model is shown in Fig. 11.

As may be seen in Fig. 12, the upper and lower parts of the tube are deployed in a symmetrical way in the first stages of the deployment. However, the symmetry is broken as soon as one crosses the limit point corresponding to the maximum pressure on the deployment curve. The displacement field of the tube then bifurcates from the symmetrical branch and the lower half of the tube is completely deployed before the upper one. Further investigations could be made to control the symmetrical and asymmetrical branches in the bifurcation context but they fall out of the scope of this work.

## 7. Experimental validation

This section describes the experimental works conducted in order to validate the numerical results obtained in the previous Section. Experiments on deployment have been investigated by a large number of authors in the dynamical context, see, e.g., Free land and Veal (1998), and Miyasaki and Uchiki (2002). The study presented here differs from these works in that the deployment is conducted in a quasi static way, as done in Buytet et al. (2006). We consider the membrane tube folded in two, having the same geometries and material properties as in the numerical part. The tube has a length of 80 cm, a diameter of 15 cm or 20 cm, it is made of a flexible and isotropic thin film (50 µm thick) with Young's modulus *E* 5.10⁹ Pa and Poisson's ratio *v* 0.3.

At the beginning, the tube is completely flattened and folded in two. During the deployment, the bottom of the tube is kept in fixed position by a pivot pin and the top is attached to an airborne carriage mounted on a fixed vertical axis, see Figs. 13 and 14. The air borne carriage is designed so as to allow the top of the tube to move freely with as little friction as possible. A vertical dead weight *F* is applied on the top of the tube in order to stabilize the deployment in quasi statics.

The deployment takes place as air is injected in the tube and the tube increases in volume. For each given quantity of air, the position of the top is measured when the internal pressure inside the tube is stable. The quantity of the inflating air is increased until the tube is almost completely deployed, namely when the top displacement is about 65 cm. From this moment on, the tube is gradually deflated until it returns to the initial flattened and folded configuration. Fig. 15 displays the pressure versus the top displacement for the 15 cm and 20 cm diameter tubes, under four dead loads *F* = 6.78 N, 13.62 N, 20.11 N and 26.9 N. For a given geometry and dead load, several tests are carried out, giving rise to slightly different response curves in Fig. 15.

It is noteworthy that wrinkles may appear virtually randomly all over the membranes throughout the deployment process. For a given geometry and dead load, the distribution of the wrinkles changes from one test to the next; however, this has insignificant influence on the global response of the deployment.

As can be seen in Fig. 15, the inflation and deflation curves do not coincide, giving rise to a hysteresis, contrary to what can be expected for a conservative system. For a given position of the tube top, the pressure in the inflation stage is always higher than that in the deflation stage. The same phenomenon is observed in all tubes, for all applied dead loads and for various materials used.

One may think that the hysteresis is due to the fact that some friction can exist between the airborne carriage and the vertical

<figure><img src="image_15.png" /><figcaption>Fig. 14. Experimental device.</figcaption></figure>
<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 15. Experimental deployment curves.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 17. Correction of the vertical dead load.</figcaption></figure>

axis, or over the surfaces where the membrane comes into self contact. However, that should not be the cause since a lot of pre cautions have been taken in the experimental settings in order to diminish the friction: in addition to the above mentioned air cushion system, the surfaces of the membrane are lubricated with oil. Also, the vertical dead load applied to the top of the tube is precisely measured by means of a force sensor inserted between the dead load and the top of the tube (see the load gauge in Fig. 13). The measured load is plotted in Fig. 16 in the case of the 20 cm diameter tube submitted to F = 6.78 N and 26.9 N. It is found that the measured load is equal to the theoretical load F in the inflation stage, but it is slightly lowered in the deflation stage, about 0.4 N for F = 26.9 N. Additional masses are then used to compensate for this difference so that the same dead load is actually applied during the inflation and the deflation, as shown in Fig. 17. In spite of this correction, the hysteresis observed in Fig. 15 remains unchanged. The conclusion is that the friction introduced by the carriage, if any, cannot be accountable for this hysteresis.

<figure><img src="image_19.png" /><figcaption>Fig. 16. Actual vertical dead load during the inflation and the deflation stages in the 20 cm-diameter tube. (a) F 6.78 N; (b) F 26.9 N.</figcaption></figure>

<figure><img src="image_20.png" /></figure>
<figure><img src="image_21.png" /><figcaption>**Fig. 18.** Numerical versus experimental deployment curves.</figcaption></figure>

A close look at the membrane surface reveals that the wrinkles are pronounced in the beginning of the inflation stage and most of them disappear when pressure *p* gets across its peak value. In the subsequent deflation stage, the membrane is always smoother with fewer wrinkles. The wrinkles appearing in the first stage are so marked that the membrane could be plastified at these locations and the mechanical energy used to smooth out the wrinkles may be the cause of the observed hysteresis.

Fig. 18 shows that the numerical results obtained in Section 6 agree very well with the experimental results. In particular, the limit point in the inflating pressure *p* encountered at halfway of the full deployment is predicted with a good precision. The difference between the numerical and the experimental pressure in creases with the dead load; its maximal value is about 10%.

## 8. Conclusions

The deployment of inflatable structures in quasi static conditions is a difficult problem as it involves different types of nonlinearities such as the large deformations, the inflating pressure which is a follower loading, the contact and self contact phenomenon; and

most importantly the difficulty due to the singularity of the structural stiffness during the deployment stage. The minimization formulation combined with the method of projected gradient has been proposed as a simple tool to solve the problems of this type. The efficiency of the approach has been proved through the case of a simply folded tube as the numerical results obtained have been found to be in good accordance with the experimental results.

## Acknowledgement

The National Space Center (Centre National d'Etudes Spatiales, CNES) in Toulouse, France, is gratefully acknowledged for having partly supported this work. The authors would also like to thank K. A. Brakke for his Surface Evolver freeware (Brakke, 1992) which has been very helpful throughout all numerical developments and implementations in this paper.

## References

- Benson, D.J., Hallquist, J.O., 1990. A single surface contact algorithm for the post-buckling analysis of shell structures. Computer Methods in Applied Mechanics and Engineering 78 (2), 141-163.

- Bonet, J., Wood, R., Mahaney, J., Heywood, P., 2000. Finite element analysis of air supported membrane structures. Computer Methods in Applied Mechanics and Engineering 190 (5-7), 579-595.

- Bouzidi, R., Levan, A., 2004. Numerical solution of hyperelastic membranes by energy minimization. Computers & Structures 82 (23-26), 1961-1969 (Computational Structures Technology).

- Brakke, K., 1992. The surface evolver. Experimental Mathematics 1 (2), 141-165.

- Buytet, S., Bouzidi, R., Dupuy, C., 2006. Deployment of membranous tubes by air inflation at low pressure. In: Topping, B., Montero, G., Montenegro, R. (Eds.), Proceedings of the Eighth International Conference on Computational Structures Technology. Civil-Comp Press, Las Palmas de Gran Canaria, Spain.

- Clem, A., Smith, S., Main, J., 2000. A pressurized deployment model for inflatable space structures. In: Proceeding of the 41st AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference, Atlanta.

- Freeland, R., Veal, G., 1998. Significance of the inflatable antenna experiment technology. In: Proceeding of the 39th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference, Long Beach.

- Gil, A.J., 2006. Structural analysis of prestressed Saint Venant-Kirchhoff hyperelastic membranes subjected to moderate strains. Computers & Structures 84 (15-16), 1012-1028.

- Goldstein, A., 1964. Convex programming in Hilbert space. Bulletin of the American Mathematical Society 70, 709-710.

- Harnau, M., Konyukhov, A., Schweizerhof, K., 2005. Algorithmic aspects in large deformation contact analysis using 'solid-shell' elements. Computers & Structures 83 (21-22), 1804-1823.

- Heegaard, J.-H., Curnier, A., 1993. An augmented Lagrangian method for discrete large-slip contact problems. International Journal for Numerical Methods in Engineering 36 (4), 569-593.

- Kane, C., Repetto, E., Ortiz, M., Marsden, J., 1999. Finite element analysis of nonsmooth contact. Computer Methods in Applied Mechanics and Engineering 180 (1-2), 1-26.

- Kikuchi, N., Song, Y., 1981. Penalty finite element approximation of a class of unilateral problems in linear elasticity. Quarterly of Applied Mathematics 39, 1-22.

- Laursen, T., 2002. Computational Contact and Impact Mechanics: Fundamentals of Modeling Interfacial Phenomena in Nonlinear Finite Element Analysis. Springer (Engineering Online Library).

- Levitin, E., Polyak, B., 1966. Constrained minimization methods. U.S.S.R. Computational Mathematics and Mathematical Physics 6 (1), 1-50.

- Miyasaki, Y., Uchiki, M., 2002. Deployment dynamics of inflatable tube. In: Proceeding of the 43rd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference, Denver.

- Puso, M.A., Laursen, T.A., 2004. A mortar segment-to-segment frictional contact method for large deformations. Computer Methods in Applied Mechanics and Engineering 193 (45-47), 4891-4913.

- Salama, M., Kuo, C., Lou, M., 1999. Simulation of the deployment dynamics of inflatable structures. In: Proceeding of the 40th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference, Reston.

- Wang, J., Johnson, A., 2002. Deployment simulation of ultra-lightweight inflatable structures. In: Proceeding of the 43rd AIAA/ASME/ASCE/AHS/ASC structures, Structural Dynamics and Materials Conference, Denver.

- Wang, S., Nakamachi, E., 1997. The inside-outside contact search algorithm for finite element analysis. International Journal for Numerical Methods in Engineering 40, 3665-3685.

- Wriggers, P., 2002. Computational Contact Mechanics. J. Wiley & Sons, Ltd.