# Numerical solution of hyperelastic membranes by energy minimization

Rabah Bouzidi, Anh Le Van

► To cite this version:

Rabah Bouzidi, Anh Le Van. Numerical solution of hyperelastic membranes by energy minimization. Computers & Structures, 2004, 82 (23-26), pp.1961 - 1969. 10.1016/j.compstruc.2004.03.057. hal-01007135

HAL Id: hal-01007135

https://hal.science/hal-01007135

Submitted on 18 Nov 2016

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>
# Numerical solution of hyperelastic membranes by energy minimization

Rabah Bouzidi, Anh Le van

Laboratory of Civil Engineering, Faculty of Sciences, University of Nantes, 2, rue de la Houssinière BP 92208,
44322 Nantes cedex 3, France

In this work, a numerical approach is presented for solving problems of finitely deformed membrane structures made of compressible hyperelastic material and subjected to external pressure loadings. Instead of following the usual finite element procedure that requires computing the material tangent stiffness and the geometric stiffness, here we solve the membrane structures by directly minimizing the total potential energy, which proves to be an attractive alternative for inflatable structures.

The numerical computations are performed over two simple geometries: the circular and the rectangular membrane and over a more complex structure: a parabolic antenna using the Saint Venant Kirchhoff and neo Hook element models. Whenever available, analytical or semi analytic solutions are used to validate the finite element results.

*Keywords:* Energy minimization; Triangular finite element; Circular membrane; Rectangular membrane; Parabolic antenna

## 1. Introduction

In this work, a numerical approach is presented for solving problems of membrane structures subjected to an external pressure loading. The problem is both geo metrically nonlinear due to finite deformations and materially nonlinear through a hyperelastic constitutive relationship. Here, the structures are assumed to be made of a quasi incompressible hyperelastic material de

scribed by either the Saint Venant Kirchhoff or a neo Hookean type models.

Usually, the classical finite element method is used in order to solve such nonlinear problems. It involves an iterative scheme to satisfy the equilibrium equations and requires computing the material tangent stiffness and the geometric stiffness. Here, the membrane struc tures are solved by directly minimizing the total potential energy. Whereas the proposed approach is theoretically equivalent to the traditional finite element method, it proves to be an attractive alternative which is particu larly efficient for complex inflatable structures.

All the numerical computations are performed using a three dimensional flexible triangular finite element, which is able to handle the above mentioned nonlinearities. The finite element has no bending stiffness and the
tangent matrix is singular for the rotation degrees of freedom. The numerical solution is carried out by means of an iterative method like the conjugate gradient or Newton one.

The proposed finite element is validated through the well known Hencky problem [1] related to the inflation of a circular membrane clamped at its rim. The obtained numerical results will be compared to the standard non linear finite element solution, using the total Lagrangian formulation and membrane elements, and a semi analyt ical solution for the Saint Venant Kirchhoff model.

The second validation is related to problem of an orthotropic rectangular membrane. A semi analytical solution will be presented in the case of orthotropic behavior. This solution will be used to validate the finite element in the case of isotropic behavior.

In order to prove the ability of the proposed ap proach to deal with large scale problems, the final appli cation is carried out on a parabolic antenna made of two separate parts subjected to different pressures. By vary ing the applied pressures, it is shown that different bifur cated folding modes can be obtained on the deformed configurations.

## 2. The approach of total potential energy minimization

### 2.1. The potential energy of the external loading

Consider a body which undergoes the deformation $\phi$ or the displacement $\vec{U}$, carrying positions $X$ in the refer ence configuration $\Omega_0$ to positions $x$ $\phi(X)$ $X+\vec{U}$ in the current configuration $\Omega$.

The reference boundary surface $\partial\Omega_0$ of $\Omega_0$ is decom posed into disjoint parts $\partial\Omega_{0U}$ and $\partial\Omega_{0T}$, over which the boundary conditions are written as

$$ \vec{U}(X) = \vec{U}^g(X) \quad \text{on} \quad \partial\Omega_{0U} \qquad (1) $$

$$ \vec{T} = \Pi \vec{N} = \vec{T}^g \quad \text{on} \quad \partial\Omega_{0T} \qquad (2) $$

where $\vec{U}^g$ denotes given (prescribed) displacement, $\vec{T}^g$ the given nominal traction vector, $\Pi$ the first Piola Kirchhoff stress tensor (not symmetric) related to the Cauchy stress tensor $\sigma$ by $\Pi \mathbf{F}^{-T} \mathbf{J}\sigma$, where $\mathbf{F} = \mathbf{I} + \nabla_X \vec{U}$ is the deformation gradient, $\mathbf{J} \det \mathbf{F}$. Use will also be made of the second Piola Kirchhoff stress tensor $\Sigma$ (symmetric) related to the first one by $\Pi \mathbf{F} \Sigma$.

Nanson's formula [5,8] for the transformation of sur face elements shows that the nominal traction vector $\vec{T}$ is in general a deformation sensitive loading depending on displacement $\vec{U}$ and $\mathbf{F}$. In the case of an inflated mem brane, one part of $\partial\Omega_{0T}$ (the outer side of the membrane) is traction free while the remaining part (the inner side) denoted by $S_0$ is subjected to a pressure $p$, so that traction vector $\vec{T}$ there writes

$$ \vec{T} = \Pi \vec{N} = \vec{T}^g = p \mathbf{J} \mathbf{F}^{-T} \vec{N} \quad \text{on} \quad S_0 \qquad (3) $$

### 2.1.1. Conservative surface loading [7]

A prescribed nominal surface loading $\vec{T}^g(\vec{X}, \vec{U}, \mathbf{F})$ de fined over the reference surface $S_0$ derives from external potential energy $V_{\text{ext}}(\vec{U})$ if $\forall \vec{v}^*$ virtual velocity field satis fying $\vec{v}^* = \vec{0}$ over $S_{0U}$, the virtual power of the surface loading can be expressed as:

$$ \wp^*(\vec{T}^g) = \int_{S_0} \vec{T}^g(\vec{X}, \vec{U}, \mathbf{F}) \vec{v}^*(X) \, dS_0 = \frac{\partial V_{\text{ext}}}{\partial \vec{U}}[\vec{v}^*] \qquad (4) $$

where $\frac{\partial V_{\text{ext}}}{\partial \vec{U}}$ is the Gateaux derivative of $V_{\text{ext}}$, which is related to the first variation of $V_{\text{ext}}$ by: $\forall \delta \vec{U}, \delta V_{\text{ext}} = \frac{\partial V_{\text{ext}}}{\partial \vec{U}}[\delta \vec{U}]$.

When (4) is satisfied, $\vec{T}^g$ is referred to as a conservative surface loading.

### 2.1.2. Definition of the pressurized volume $\theta$

From (1), the boundary line $\partial S_0$ of $S_0$ is subjected to the prescribed displacements

$$ \forall \vec{X} \in \partial S_0, \quad \vec{U}(X) = \vec{U}^g(X) \qquad (5) $$

Let us denote by $S_{\text{lat}}$ the surface swept in the three dimensional space, up to current time, by line $\partial S$. The pressurized volume $\theta$ is then defined as that bounded by surfaces $S_0$, $S$ and $S_{\text{lat}}$ at current time (Fig. 1a).

More often than not the prescribed displacements $\vec{U}^g$ are zero, so that $S_{\text{lat}}$ is the empty set and $\theta$ is merely the volume between $S_0$ and $S$ (Fig. 1b).

### 2.1.3. Properties of volume $\theta$

* The pressurized volume $\theta$ is a functional of the dis placement field $\vec{U}$. In order to express this, we write $\theta = \theta(\vec{U})$.

* In general, the explicit expression for $\theta(\vec{U})$ as a func tional of $\vec{U}$ is unknown. Nevertheless, its first varia tion $\delta\theta$ corresponding to a variation $\delta\vec{U}$ of $\vec{U}$ can be determined by the following relation:

$$ \forall \delta \vec{U}, \quad \delta\theta = \int_{S_0} \delta \vec{U} \mathbf{J} \mathbf{F}^{-T} \vec{N} \, dS_0 \qquad (6) $$

<figure><img src="image_2.png" /><figcaption>Fig. 1. (a) Definition of pressurized volume $\theta$. The reference surface $S_0$ maps into the current surface $S$. (b) Case where the prescribed displacements are zero.</figcaption></figure>
### 2.1.4. Pressure potential energy

Assuming that the whole line $\partial S$ is subjected to pre scribed displacements (5) and that pressure $p$ is uniform over surface $S$, then the pressure loading derives from the potential [6]:

$$V_{\text{ext}}(\vec{U}) = \int p(\theta) d\theta \quad (7)$$

where the symbol $\int$ stands for a primitive. In particular, if the pressure is so controlled that it is independent of volume $\theta$, then the pressure potential energy is $V_{\text{ext}}(\vec{U}) = p\theta(\vec{U})$. If the pressure varies as a function of volume $\theta$ according to the law $p\theta K$ constant, then the potential energy is $V_{\text{ext}}(\vec{U}) = K \log \theta(\vec{U})$.

For numerical purposes, we assume pressure $p$ to be constant. The membrane surface is discretized into trian gular finite elements (as shown in Fig. 2). In order to compute the pressurized volume $\theta$ we transform it into surface integrals by means of the Gauss theorem:

$$\theta = \frac{1}{3} \int_{S} \vec{x} \cdot \vec{n} \, dS \quad (8)$$

where $\vec{x}$ is the current vector position of a particle on the surface $S$, $\vec{n}$ the outward normal to $S$. The surface integral (8) becomes a discrete summation over the finite elements:

$$V_{\text{ext}} = \frac{1}{3} \sum_{\text{elements}} p_i S_i \sum_{i=1}^{3} \left[ \frac{x_i^a + x_i^b + x_i^c}{3} n_i \right] \quad (9)$$

## 2.2. The internal potential energy (strain energy)

We consider a hyperelastic material described by a strain energy $\psi$ per unit reference volume (rather than per unit mass), function of either the right stretch tensor $C = g_{ij} \vec{G}^i \otimes \vec{G}^j$ or the Green strain tensor $E = \frac{1}{2}(C \cdot I)$. Coefficients $g_{ij}$ are the covariant components of metric tensor $g$ defined in the current configuration and vectors $\vec{G}^i$ form the dual of the natural basis in the reference configuration.

The first and second Piola Kirchhoff stress tensors $\Pi$ and $\Sigma$, respectively, are given by

$$\Sigma = \frac{\partial \psi(\mathbf{E})}{\partial \mathbf{E}} = 2 \frac{\partial \psi(\mathbf{C})}{\partial \mathbf{C}} \quad \Pi^T = \frac{\partial \bar{\psi}(\mathbf{F})}{\partial \mathbf{F}} \quad (10)$$

<figure><img src="image_3.png" /><figcaption>Fig. 2. Deformation of a triangle element.</figcaption></figure>

where $\tilde{\psi}(\mathbf{F}) \psi(C \mathbf{F}^T \mathbf{F})$. By introducing the strain energy

$$V_{\text{int}}(\vec{U}) = \int_{\Omega_0} \bar{\psi}(\mathbf{F}(\vec{U})) d\Omega_0 \quad (11)$$

the virtual stress power along a virtual velocity field $\vec{v}^*$ can be expressed as a function of the Gateaux derivative of $V_{\text{int}}$

$$ \begin{aligned} \mathcal{G}_{\text{int}}^* &\equiv \int_{\Omega} \boldsymbol{\sigma} : \nabla_x \vec{v}^* d\Omega = \int_{\Omega_0} \mathbf{\Pi}^T : \nabla_x \vec{v}^* d\Omega_0 \\ &= \int_{\Omega_0} \frac{\partial \bar{\psi}(\mathbf{F})}{\partial \mathbf{F}} : \nabla_x \vec{v}^* d\Omega_0 = \frac{\partial V_{\text{int}}}{\partial \vec{U}}[\vec{v}^*] \end{aligned} \quad (12) $$

From relation (12), one says that the Piola Kirchhoff stress $\mathbf{\Pi}$ derives from internal potential energy $V_{\text{int}}(\vec{U})$.

For an isotropic material without internal con straints, the strain energy per unit volume $\psi$ can be expressed in terms of the invariants of the right Cauchy Green tensor $\mathbf{C}$

$$\psi = \psi(I_1(\mathbf{C}), I_2(\mathbf{C}), I_3(\mathbf{C})) \quad (13)$$

with $I_1(\mathbf{C}) \text{ tr } \mathbf{C}$, $I_2(\mathbf{C}) = \frac{1}{2}((\text{tr } \mathbf{C})^2 \text{ tr } \mathbf{C})^2$, $I_3(\mathbf{C}) \text{ det } \mathbf{C}$ $J^2$.

Two constitutive models are used for numerical purposes:

(i) the Saint Venant Kirchhoff model, for which the strain energy is

$$\psi = \frac{\lambda}{2} (\text{tr } \mathbf{E})^2 + \mu \mathbf{E} : \mathbf{E} \quad (14)$$

where $\lambda$ and $\mu$ are material parameters;

(ii) and the compressible isotropic model described by a strain energy of neo Hookean type

$$\psi = \frac{\mu}{2} (\text{tr } \mathbf{C})^2 + \frac{\lambda}{2} (\text{tr } \mathbf{J})^2 \quad (15)$$

which yields from the constitutive law (10)

$$\Sigma = \mu(\mathbf{I} \cdot \mathbf{C}^{-1}) + \frac{\lambda}{2} \ln \det \mathbf{C} \cdot \mathbf{C}^{-1} \quad (16)$$

If the body undergoes a rigid body motion, the defor mation gradient $\mathbf{F}$ is equal to an orthogonal tensor, ten sor $\mathbf{C}$ is equal to the identity tensor and the stress (16) is zero, as expected.

Use will be also made of Young's modulus $E$ and Poisson's ratio $\nu$ related to parameters $(\lambda, \mu)$ as in small strains: $\lambda E\nu/(1+\nu)(1-2\nu)$, $\mu E/(2(1+\nu))$.

When using the triangular finite elements (Fig. 2), the mixed components of $\mathbf{C}$ are defined by

$$C_i^j = d_{ik} D^{kj} \quad (17)$$

where $D$ and $d$ are the metric tensors in the reference and current configuration with their covariant components defined as:
$$[D] = \begin{bmatrix} \vec{A}B \cdot \vec{A}B & \vec{A}B \cdot \vec{A}C \\ \vec{A}B \cdot \vec{A}C & \vec{A}C \cdot \vec{A}C \end{bmatrix} \quad [d] = \begin{bmatrix} \vec{a}b \cdot \vec{a}b & \vec{a}b \cdot \vec{a}c \\ \vec{a}b \cdot \vec{a}c & \vec{a}c \cdot \vec{a}c \end{bmatrix} \quad (18)$$

Since the strain tensor is constant inside each finite element, the strain energy is simply obtained as the product of the strain energy by the volume of the element

$$V_{\text{element}}^{\text{int}} = \int_{\text{element}} \psi \, d\Omega_0 = \psi S_0 h \quad (19)$$

where $S_0$ is the reference area of the membrane element and $h$ the reference thickness.

## 2.3. Minimization of the total potential energy

The minimization approach is based on the following result.

**Proposition. If:**

(i) the external loading derives from an external potential energy $V_{\text{ext}}$ (e.g. relation (9)),

(ii) and the constitutive law is hyperelastic, so that the stresses derive from the internal potential energy $V_{\text{int}}$ (11).

Then, the total potential energy $V \equiv V_{\text{ext}} + V_{\text{int}}$ is stationary at the solution displacement field $\vec{U}$.

**Proof.** The principle of virtual power gives: $\forall \vec{u}^* \text{ CAH}$,

$$\begin{align*} 0 &= \varphi_{\text{ext}}^* + \varphi_{\text{int}}^* \\ &\equiv \int_{S_{0r}} \vec{T}^d(\vec{X}, \vec{U}, \mathbf{F}) \cdot \vec{v}^*(X) \, dS_0 \\ &\qquad \int_{\Omega_0} \mathbf{\Pi}^T : \nabla_X \vec{v}^* \, d\Omega_0 = \frac{\partial(V_{\text{ext}} + V_{\text{int}})}{\partial \vec{U}} [\vec{v}^*] \quad (20) \end{align*}$$

Relation (20) shows that the first variation of $V$ is $\delta V = \frac{\partial V}{\partial \vec{U}} [\delta \vec{U}] = 0$. $\square$

Furthermore it is shown that a stable solution $\vec{U}$ corresponds to a minimum of $V$.

The minimization algorithm is rather classical. It is based on the descent method like the conjugate gradient or the Newton method. Nevertheless, we insist on the fact that the gradient of the total potential, which gives the direction of descent, must be computed in exact way. Thus, it is necessary to provide the analytical expression of the gradient of the potential in order to have the sufficient accuracy and to correctly handle some phenomenon like membrane folding.

## 3. Inflation of an isotropic circular membrane

In order to validate the proposed numerical model, we consider the Hencky's problem [1, 4], which consists

in a circular membrane of (initial) radius $a$ and thickness $h$, clamped on its rim and subjected to lateral pressure $p$. Here it is assumed that the material obeys either the Saint Venant Kirchhoff or neo Hookean models (Eqs. (14) and (15)).

The obtained numerical results will be compared to:

(i) The standard nonlinear finite element solution, using the total Lagrangian formulation and membrane elements (six node triangles and eight node quadrilaterals).

(ii) Fichter's semi analytical solution [2] for the Saint Venant Kirchhoff model.

Analytical solutions for circular membranes made in incompressible isotropic materials can be found in [10,11]. For compressible materials, analytical solutions are rather few since the absence of isochoric constraints leads to more complicated kinematics. A review of solution strategies for compressible isotropic materials was presented by Horgan in Chapter 4 of [8]. When solving a circular membrane, Fichter [2] dropped some second order terms in the Green strain components and considered the pressure as a dead load. These approximations led to a simplified solution which is chosen here for comparison in moderate rotations.

### 3.1. Fichter's semi analytical solution

Let us summarize Fichter's semi analytical solution given in [2]. For brevity's sake, the dimensionless radial coordinate $r/a$ will be denoted $\rho$. The deflection $w$ is searched in the form of a power series:

$$w(\rho) = a \left(\frac{pa}{Eh}\right)^{1/3} \sum_{0}^{\infty} a_{2n} (1 - \rho^{2n+2}) \quad (21)$$

By replacing this expression in the equilibrium equation, one obtains the circumferential stress $\Sigma_{00}$

$$\Sigma_{00}(\rho) = \frac{E}{4} \left(\frac{pa}{Eh}\right)^{2/3} \sum_{0}^{\infty} (2n + 1) b_{2n} \rho^{2n} \quad (22)$$

Taking into account the boundary condition $w(a) = 0$ enables one to determine the coefficients $a_{2n}$ and $b_{2n}$. The values up to $n=10$ are given in [2]

$$\begin{align} b_2 &= \frac{1}{b_0^3} & b_4 &= \frac{2}{3b_0^3} & b_6 &= \frac{13}{18b_0^8} \notag \\ b_8 &= \frac{17}{18b_0^{11}} & b_{10} &= \frac{37}{27b_0^{14}} & b_{12} &= \frac{1205}{567b_0^{17}} \notag \\ b_{14} &= \frac{219241}{63504b_0^{20}} & b_{16} &= \frac{6634069}{1143072b_0^{23}} & b_{18} &= \frac{51523763}{5143824b_0^{26}} \notag \\ b_{20} &= \frac{998796305}{56582064b_0^{29}} & & & & \tag{23} \end{align}$$
comparison between these two solutions is presented for isotropic behavior case only.

## 4.1. Semi analytical solution

Let us consider a rectangular membrane with sides 2a and 2b, clamped at its rim, subjected to a uniform pressure p and made of an orthotropic elastic material. In order to simplify the subsequent equations, Cartesian co ordinates X and Y will be replaced by dimensionless quantities $\xi$ X/a and $\eta$ Y/b.

We build the semi analytical solution in the large deflection case: the in plane displacements of the membrane $u(\xi, \eta)$, $v(\xi, \eta)$ and the deflection $w(\xi, \eta)$. In order to automatically satisfy the kinematic boundary conditions, we assume the following polynomial expressions for the displacement components:

$$ 
\begin{align*}
u(\xi, \eta) &= g(\xi, \eta) \sum_{i=1}^{n} a_{2i-1} \xi^{2i-1} \\
v(\xi, \eta) &= g(\xi, \eta) \sum_{i=1}^{n} b_{2i-1} \eta^{2i-1} \\
w(\xi, \eta) &= g(\xi, \eta) \left[ w_0 + \sum_{i=1}^{n} c_{2i} \xi^{2i} + d_{2i} \eta^{2i} \right]
\end{align*}
\qquad (26) $$

$g(\xi, \eta)$ is chosen so that the boundary kinematic conditions are satisfied.

$$ g(\xi, \eta) = (1 \quad \xi^2)(1 \quad \eta^2) $$

The orthotropy is taken into account with distinct coefficients for u and v displacement components. We also introduced different coefficients in the expression of the deflection w for X and Y directions. Minimizing the total potential energy provides the coefficients of the polynomial series and the displacement field in the rectangular membrane. For this we need to evaluate the finite strain components:

$$ 
\begin{align}
E_{XX} &= \frac{1}{a} \frac{\partial u}{\partial \xi} + \frac{1}{2a^2} \left( \frac{\partial w}{\partial \xi} \right)^2 \nonumber \\
E_{XX} &= \frac{1}{b} \frac{\partial v}{\partial \eta} + \frac{1}{2b^2} \left( \frac{\partial w}{\partial \eta} \right)^2 \tag{27} \\
E_{XY} &= \frac{1}{2b} \frac{\partial u}{\partial \eta} + \frac{1}{2a} \frac{\partial v}{\partial \xi} + \frac{1}{2ab} \frac{\partial w}{\partial \xi} \frac{\partial w}{\partial \eta} \nonumber
\end{align}
$$

As a matter of fact, the large strain terms in the above definition have been omitted in order to facilitate the obtention of a semi analytical solution. Consequently, while the proposed semi analytical solution is able to account for large rotations, the comparison is only valid for moderate strains. On the other hand, the finite element procedure described below includes all the nonlinear terms and should be more accurate than the semi analytical solution. Recall that the elasticity matrix for an elastic orthotropic material writes:

$$ 
\left\{ 
\begin{array}{c} 
\Sigma_{XX} \\ 
\Sigma_{YY} \\ 
\Sigma_{XY} 
\end{array} 
\right\} = \begin{bmatrix} D_{XX} & D_{XY} & 0 \\ D_{XY} & D_{YY} & 0 \\ 0 & 0 & 2G_{XY} \end{bmatrix} \left\{ 
\begin{array}{c} 
E_{XX} \\ 
E_{YY} \\ 
E_{XY} 
\end{array} 
\right\} \qquad (28) $$

The $D_{ij}$ terms depend on the material constants $E_X$, $E_Y$ and $v_{XY}$

$$ 
\begin{align}
D_{XX} &= \frac{E_X h}{v_{XY} v_{YX}}, & D_{YY} &= \frac{E_Y h}{v_{YX} v_{XY}}, \nonumber \\
D_{XY} &= \frac{v_{XY} E_Y h}{v_{YX} v_{XY}} \tag{29}
\end{align}
$$

where $E_X$ is the Young modulus in the length direction, $E_Y$ the Young modulus in the width direction, $G_{XY}$ the shearing modulus, $h$ the membrane thickness and $v_{XY}$ the Poisson ratio ($v_{XY}/E_X$, $v_{YX}/E_Y$). Then, the strain energy density for an orthotropic material takes the form:

$$ \psi = \frac{1}{2} (D_{XX} E_{XX}^2 + 2D_{XY} E_{XX} E_{YY} + D_{YY} E_{YY}^2 + 4G_{XY} E_{XY}^2) \qquad (30) $$

The strain energy is obtained after integration of the strain energy density over the membrane:

$$ V_{\text{int}} = hab \int_{-1}^{1} \int_{-1}^{1} \psi \, dX \, dY \qquad (31) $$

A symbolic calculus of $V_{\text{int}}$ is carried out up to the fifth order terms of the displacement series. The solution is not presented here because of its complexity. However at this stage, we have evaluated the strain energy in the case of a large deflection without any numerical approximation like linearization or other simplification usually used in the classical finite elements method. The only terms which were neglected are those of the large strains (Eq. (27)). Taking into account these terms makes it very difficult to obtain the analytical expression of $V_{\text{int}}$ by means of the actual symbolic calculation systems.

## 4.2. Pressure potential energy

The pressure potential energy is the product of the pressure by the integral of the deflection on the membrane. Rigorously, we must evaluate the exact volume of the membrane by taking into account the horizontal displacements. It means that it the pressure work only on the deflection component. In this case, the pressure vector is assumed to still parallel to z during the inflation and do not follow the membrane normal. This approximation affects slightly the shape of the membrane and gives less rounded surface. Taking into account the dimensionless coordinates, the pressure potential energy can be written in the following form:

$$ 
\begin{align}
V_{\text{ext}} &= P 2a2b \int_{-1}^{1} \int_{-1}^{1} w(X, Y) \, dX \, dY \nonumber \\
&= 16pab \sum_{i=1}^{n} \frac{c_{2i} + d_{2i}}{4i^2} \frac{1}{1} \tag{32}
\end{align}
$$
where *n* is the chosen degree of the polynomial series. This expression of the pressure potential is obtained with the help of a symbolic calculus tools. The displacement is computed by minimizing the total potential energy $V_{tot} = V_{int} + V_{ext}$ relatively to the coefficients of the power series:

$$ \frac{\partial V}{\partial a_i} = 0 \quad \frac{\partial V}{\partial b_j} = 0 \quad \frac{\partial V}{\partial c_k} = 0 \quad \frac{\partial V}{\partial d_l} = 0 \quad \text{for all } i, j, k, l \qquad (33) $$

The analytic expression of $V_{tot}$ upto *n* 5 (Eq. (26)) is obtained by means of symbolic calculus tools. The minimization is achieved numerically by a Newton like iterative scheme or the conjugate gradient method. All the coefficients depend on the following parameters: the pressure *p*, the dimensions 2a and 2b of the membrane and the material properties.

In fact, the initial purpose of this analytic solution is to use it in an inverse analysis. The symbolic expressions of the displacement's polynomial coefficients have been obtained up to *n* 1. However, this analytic solution is not presented here until more complete developments.

### 4.3. Numerical results

In this section, we present a comparison between the orthotropic semi analytical solution presented in the previous section and the triangular finite element on an isotropic rectangular membrane. Given the material properties of the membrane, we can evaluate the coefficients of the displacement series by a numerical minimization of the total potential energy. We used some data issued from simple tension experiments for measurement of Young's modulus listed in Table 1.

#### 4.3.1. Semi analytical solution results

The minimization of the total potential energy (33) with the semi analytical solution gives the coefficients of the polynomial series (26). The obtained values are listed in the Table 2.

In spite of the isotropic behaviour that we retained for this study, one can note that the coefficients obtained are not isotropic by permutation of axes *X* and *Y*. This orthotropism of the shape is due only to the dimension of the plate which are not the same according to two dimensions.

#### 4.3.2. Triangular finite element solution results

We have modelled the same rectangular membrane with the triangular finite element. Fig. 5 shows one of the meshes used to discretize the membrane.

For the results analysis, we have considered the convergence of the results obtained by the proposed finite element. Fig. 6 shows the dependence of the central deflection and the total potential energy on the mesh refinement. Convergence on displacement is faster than convergence on total energy. It is suitable to not penalize the process of convergence by an energy criterion with a great restriction. For the results presented here, the criterion of convergence relates to total energy. The variation of the energy between successive iteration must be less than $10^{-6}$ Nm.

<figure><img src="image_4.png" /><figcaption>Fig. 5. Rectangular membrane simulation.</figcaption></figure>

<table><caption>Table 1
Isotropic membrane data</caption><thead><tr><th>p (Pa)</th><th>a (m)</th><th>b (m)</th><th>E<sub>Xh</sub> (Pa)</th><th>E<sub>Yh</sub> (Pa)</th><th>v<sub>XY</sub></th><th>G<sub>XY</sub></th></tr></thead><tbody><tr><td>60,000</td><td>0.15</td><td>0.20</td><td>80,000</td><td>80,000</td><td>0.3</td><td>5000</td></tr></tbody></table>

<table><caption>Table 2
Polynomial series coefficients obtained for the semi analytical solution</caption><tr><td>w<sub>0</sub></td><td>0.05716</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>a<sub>i</sub></td><td>0.09261</td><td>1.73120</td><td>18.11561</td><td>0.58344</td><td>0.01450</td><td></td></tr><tr><td>b<sub>i</sub></td><td>0.04561</td><td>1.58834</td><td>10.31184</td><td>0.22604</td><td>0.00187</td><td></td></tr><tr><td>c<sub>i</sub></td><td>0.55763</td><td>9.40701</td><td>0.44526</td><td>0.02019</td><td>0.00060</td><td></td></tr><tr><td>d<sub>i</sub></td><td>0.95364</td><td>3.95415</td><td>2.27647</td><td>0.13073</td><td>0.00547</td><td></td></tr></table>
<figure><img src="image_5.png" /><figcaption>Fig. 6. Results convergence of the triangular finite element.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 7. Relative difference between the semi analytical and triangular finite element solutions.</figcaption></figure>

### 4.3.3. Semi analytical and finite elements results comparison

We have compared the deflection $w$ obtained from the two models: semi analytical noted $w_a$ and the triangular finite element noted $w_{fe}$. In order to underline the difference between the two solutions, let's define the relative difference between them as: $(w_{fe} - w_a)/\max(w_{fe})$. We have represented on Fig. 7 the distribution of the relative difference over the membrane surface. At the center of the membrane, we obtained the following deflection: $w_{fe}$ 0.0611 m and $w_a$ 0.0572 m which gives 6.4% of relative difference.

The difference between these two solutions is due essentially to the difference of their shape but not to their amplitude. The average relative difference over the membrane is about 2.5%. The restriction on the order of the polynomial series can explain this difference.

## 5. Inflation of a parabolic antenna

As an application to large scale problems, we consider the inflation of a parabolic antenna used in space devices.

The antenna is made of two separate closed membranes, sewed together along a circumferential rim and subjected to different internal pressures $p_1$ and $p_2$. In the unstressed state, the whole structure lies in a plane. Under pressure $p_1$, the inner circular membranes become two paraboloids. One of them is transparent to the electromagnetic waves and the second one acts as a concave reflector. Under pressure $p_2$, the outer membrane becomes a torus, and the final shape of the antenna can be controlled by adjusting the two pressures. It should be noted that this study only aims to show typical results of such structures due to geometrical incompatibilities. We take the inner radius equal to 1 m, the radius of the torus is equal to 0.1 m. The material verifies the Saint Venant Kirchhoff model with the equivalent Young modulus product of the Young modulus $E$ by the thickness of the membrane $h$, $E^* \cdot E \cdot h$ 350,000 Pa and Poisson ratio $\nu$ 0.3. Fig. 8 shows the structure under pressures $p_1$ 1000 Pa and $p_2$ 10,000 Pa. The paraboloids remain quite plane and contain some shallow folds regularly distributed in the radial direction. Due to the geometric incompatibility, the torus warps out of its plane and displays more visible folds regularly distributed along the circumference.

<figure><img src="image_7.png" /><figcaption>Fig. 8. Deformed shape of the parabolic antenna. $p_1$=1000 Pa and $p_2$=10,000 Pa.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 9. Deformed shapes under higher internal pressure. $p_1$=10,000 Pa and $p_2$=100,000 Pa.</figcaption></figure>
When the inner pressure $p_1$ is increased up to 10,000 Pa and the outer pressure $p_2$ to 100,000 Pa, the parabo loids became more convex and the torus breaks down into rounded sectors separated by more marked folds (as shown in Fig. 9). This example shows the ability of the proposed approach to deal with the ill conditioned stiffness matrix when passing bifurcations from axisym metric states. Although the solution does not handle the self contact of the membranes for instance, the exter nal torus may overlap the parabolic part of the antenna at some locations the obtained deformed shapes are quite realistic. Computations including contact aspects are in progress in order to obtain more precise numerical results.

## 6. Conclusions

In this work, we have proposed a numerical ap proach to deal with hyperelastic membrane structures undergoing large deformations. The membrane surface has been divided into triangular finite elements and the solution achieved by directly minimizing the total poten tial energy, instead of satisfying the equilibrium equa tions as done with the usual finite element method.

The finite element procedure has been validated on the Hencky's problem of a circular membrane. Two compressible isotropic hyperelastic potentials have been used: the Saint Venant Kirchhoff and the neo Hookean ones. The obtained numerical results have been found to be in exact agreement with those given by the standard finite element method using the total Lagrangian formu lation and membrane elements. They also agreed very well with the semi analytical solution developed by Fichter [2], using a polynomial series developed up to order 10.

Contrary to the circular membrane case, for which expressions for the displacement coefficients are found explicitly, it is not so easy to analytically solve the elast ostatic problem of the rectangular membrane and we have to develop a semi analytical solution for the rectan gular membrane. Again, the displacement field is ex pressed in polynomial series forms. However, the order of the polynomial series is limited here to 5 because the analytic expression of the strain energy (31) becomes lengthy as the order increases. The comparison between semi analytical and finite element minimization shows a greater divergence than in the case of the circular mem

brane. Nevertheless, the maximum difference is about 8 % and the average difference is about 2.5%.

Application to the inflation of a parabolic antenna shows the ability of the proposed approach to deal with large scale problems with solutions bifurcated from axi symmetric states.

Further developments in progress prove that the con tact problem of inflatable structures encountering fric tionless obstacles can be easily dealt with by the proposed minimization technique as well.

## Acknowledgment

All the minimization computations and simulations have been done using the surface evolver code developed by Brakke [9]. We would like to thank him for this help ful tool.

## References

[1] Hencky H. On the stress state in circular plates with vanishing bending stiffness. Z Math Phys 1915;63:311 7.

[2] Fichter WB. Some solutions for the large deflections of uniformly loaded circular membranes. NASA Technical Paper 3658 NASA Langley Research Center; Hampton, VA; 1997.

[3] Campbell JD. On the theory of initially tensioned circular membranes subjected to uniform pressure. Quart J Mech Appl Math 1956;IC(Pt 1):84 93.

[4] Marker DK, Jenkins CH. Surface precision of optical membranes with curvature. Opt Exp 63 1997;1(11):311 7.

[5] Ogden RW. Non linear elastic deformations. Dover; 1997.

[6] Fischer D. Configuration dependent pressure potentials. J Elast 1998;19:77 84.

[7] Ciarlet PG. Mathematical elasticity, vol. 1: three dimen sional elasticity. Amsterdam: North Holland; 1988.

[8] Fu YB, Ogden RW, editors. Nonlinear elasticity Theory and applications. Cambridge University Press; 2001.

[9] Brakke KA. The surface evolver. Exp Math 1992;1(2): 141 65.

[10] Ogden RW. Large deformation isotropic elasticity: on the correlation of theory and experiment for incompressible rubberlike solids. Proc Royal Soc London A 1972;326: 565 84.

[11] de Souza Neto EA, Peric D, Owen DRJ. Finite elasticity in spatial description: linearization aspects with 3D mem brane applications. Int J Num Meth Engng 1995;38: 3365 81.