# Dynamic inflation of non-linear elastic and viscoelastic
rubberlike membranes

Erwan Verron, Gilles Marckmann, Bernard Peseux

► To cite this version:

Erwan Verron, Gilles Marckmann, Bernard Peseux. Dynamic inflation of non-linear elastic and viscoelastic rubberlike membranes. International Journal for Numerical Methods in Engineering, 2001, 50 (5), pp.1233-1251. 10.1002/1097-0207(20010220)50:53.0.CO;2-W. hal-01006759v2

HAL Id: hal-01006759

https://hal.science/hal-01006759v2

Submitted on 4 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>

Distributed under a Creative Commons Attribution 4.0 International License
Dynamic inflation of non-linear elastic and viscoelastic
rubber-like membranes

E. Verron*,†, G. Marckmann and B. Peseux

Laboratoire Mécanique et Matériaux, École Centrale de Nantes, BP 92101, 44321 Nantes cedex 3, France

SUMMARY

The present paper deals with the dynamic inflation of rubber-like membranes. The material is assumed to obey the hyperelastic Mooney's model or the non-linear viscoelastic Christensen's model. The governing equations of free inflation are solved by a total Lagrangian finite element method for the spatial discretization and an explicit finite-difference algorithm for the time-integration scheme. The numerical implementation of constitutive equations is highlighted and the special case of integral viscoelastic models is examined in detail. The external force consists in a gas flow rate, which is more realistic than a pressure time history. Then, an original method is used to calculate the pressure evolution inside the bubble depending on the deformation state. Our numerical procedure is illustrated through different examples and compared with both analytical and experimental results. These comparisons yield good agreement and show the ability of our approach to simulate both stable and unstable large strain inflations of rubber-like membranes.

KEY WORDS: dynamic inflation; membrane; large strains; hyperelasticity; viscoelasticity; finite element method

# 1. INTRODUCTION

Some of the earlier works on rubber-like membranes inflation dealt with static equilibrium of these inflated thin structures. The general non-linear theory for membrane deformation was firstly derived by Green and Adkins [1]. Later in the 1960s, Hart-Smith and Crisp [2] and Klingbeil and Shield [3] examined the special case of axisymmetrical hyperelastic membranes inflation. Authors proposed analytical solutions for the circular plane membrane inflation problems using different hyperelastic non-linear constitutive equations. The results obtained are compared with experimental data for rubber. At the same time, Oden and Sato [4] used the finite element method in this context. They set the foundation of a Galerkin finite element approach to solve non-linear elastic membrane inflation problems. They considered three-nodes triangular elements and the corresponding non-linear governing equations are solved using the Newton-Raphson algorithm.

<footnote>*Correspondence to: E. Verron, Laboratoire Mécanique et Matériaux, École Centrale de Nantes, BP 92101, 44321 Nantes cedex 3, France</footnote>

<footnote>†E-mail: erwan.verron@ec-nantes.fr</footnote>
Later, some authors extensively analysed viscoelastic membranes inflation cases. Wineman [5] and Feng [6] used semi-analytical time-discretization schemes to solve the time-dependent inflation of an initially plane circular non-linear viscoelastic membrane. In these papers, rubber-like materials are assumed to obey non-linear integral viscoelastic constitutive equations such as K-BKZ and Christensen's models. Shrivastava and Tang [7] developed a three-dimensional finite element method for such problems: owing to the material large strain time-dependent behaviour, described by the Christensen's model, and geometric non-linearities, the assembled non-linear equations of equilibrium are solved by the way of an iterative method.

More recently, there has also been considerable interest in the dynamic inflation, which takes into account membranes inertia. Such investigations are reviewed in details by Jenkins and Leonard [8], and Jenkins [9]. These papers show that dynamic approaches are very few in the membrane inflation field. But it is evident that such works present advantages: fast processes such as thermoforming or blow-moulding of plastics can be better simulated and numerical solution procedures are more efficient.

In this context, the present paper reports a numerical investigation on the dynamic inflation of rubber-like membranes in 3D shapes. Membranes are deformed by internal pressure under large strains and the corresponding governing equations are solved using a finite element procedure. Materials are assumed to be isotropic and incompressible, and we use two types of constitutive equations: hyperelastic and non-linear viscoelastic. The dynamic formulation is derived: the implementation of the finite element method and the temporal integration procedure are reported. The problem is solved using three-nodes triangular finite elements and the central-difference explicit time-integration scheme. The numerical implementation of the constitutive equations is highlighted and the special case of integral non-linear viscoelastic Christensen's model is carried out in detail. Moreover, we develop an original method for the determination of the pressure evolution based on a simple thermodynamic state equation of the inflating gas. Finally, different examples demonstrate the accuracy of numerical solutions by comparison with both analytical and experimental results.

## 2. RUBBER-LIKE MATERIALS

### 2.1. Preliminaries

The general constitutive equation of an isotropic, homogeneous, incompressible simple material is a tensor function of the strain history [10]:

$$ S(t) = -pC^{-1}(t) + \mathcal{S}\{C(\tau)\} \quad (1) $$

in which $C(t)$ is the right Cauchy-Green deformation tensor at time $t$, $p$ is an indeterminate hydrostatic pressure due to the incompressibility assumption, $S(t)$ stands for the second Piola-Kirchhoff stress tensor and $\mathcal{S}$ is the response functional.

Noting $F(t)$ the deformation gradient and $J(t)$ the Jacobian of the transformation (determinant of $F(t)$), the stress tensor $S(t)$ is related to the true stress tensor $\sigma(t)$ by

$$ S(t) = J(t)F^{-1}(t)\sigma(t)F^{-T}(t) \quad (2) $$

In our case, the material incompressibility implies

$$ J(t) = \det(F(t)) = 1 \quad \forall t \quad (3) $$
In the two next sections, we present two classical types of behaviour used to describe rubber-like materials such as natural rubber, elastomeric materials or high-temperature thermoplastics: hyperelastic and non-linear viscoelastic constitutive relations.

## 2.2. Hyperelastic model

First, it is helpful to consider that rubber-like materials stresses only depend on the current deformations:

$$S = -pC^{-1} + \mathcal{S}\{C} \quad (4)$$

The time $t$ can be omitted because of the elastic behaviour.

We also assume that the material is characterized by the existence of a strain-energy function, $W$. Such constitutive equations are called hyperelastic. The general stress-strain relationship is written as

$$S = -pC^{-1} + 2 \frac{\partial W}{\partial C} \quad (5)$$

Taking into account isotropic hypothesis, the strain energy can be expressed as a function of the strain invariants $I_1$, $I_2$ and $I_3$. Noting that $I_3 = J = 1$ (incompressibility), $W$ takes the following form:

$$W = W^*(I_1, I_2) \quad (6)$$

where

$$I_1 = \text{tr}(C) \quad \text{and} \quad I_2 = \frac{1}{2}[I_1^2 - \text{tr}(C^2)] \quad (7)$$

Here, we consider the special Mooney's form of the strain energy function. This model is one of the earlier and most employed hyperelastic models because of its mathematical simplicity and relatively good accuracy for reasonably large strains (≤150 per cent) [11]. The form of the strain energy function is

$$W = c_1(I_1 - 3) + c_2(I_2 - 3) \quad (8)$$

where $c_1$ and $c_2$ are the two material parameters. Note that the particular case with $c_2 = 0$ corresponds to the neo-Hookean model [12].

## 2.3. Non-linear viscoelastic model

We previously assumed that rubber-like materials are hyperelastic. In fact, the behaviour of such polymers is slightly viscoelastic [13]. Stresses at time $t$ does not only depend on the present strain state but also on the strain history. It is commonly said that the material has a memory. Consequently, the general constitutive relation of non-linear incompressible viscoelastic materials is given by Equation (1).

Concerning the non-linear viscoelasticity, the choice of models which give realistic description of the material behaviour are mathematically simple enough for numerical implementation is far from conclusion. In this paper, we consider the Christensen's theory [14], which is the simplest non-linear theory for the viscoelastic behaviour of elastomers. The author presented his theory as
‘the viscoelastic counterpart of the kinetic theory of rubber elasticity’ because, for a sufficiently slow process, the material is assumed to be neo-Hookean.

The stress–strain relationship is obtained using a first-order Green–Rivlin expansion of the response functional ($\mathcal{S}$ in Equation (1)). The second Piola–Kirchhoff stress tensor is given by

$$S(t) = -pC^{-1}(t) + g_0I + \frac{1}{2} \int_0^t g_1(t - \tau) \frac{\partial C(\tau)}{\partial \tau} d\tau \quad (9)$$

where $g_0$ is an elastic constant and $g_1(\tau)$ a relaxation function which satisfies

$$\lim_{\tau \to \infty} g_1(\tau) = 0 \quad (10)$$

### 3. FINITE ELEMENT PROCEDURE

In this section, we develop a finite element formulation for dynamic inflation of non-linear membrane problems. By definition, a membrane offers no moment or transverse shear resistance. The mid-surface of the sheet is a two-dimensional continuum and the thickness is assumed to be a function of the position of the mid-surface. In the present formulation, no difference is made between tension and compression (wrinkling phenomena are not taken into account). In reality, a membrane is unable to afford compressive stresses, therefore our analysis is valid only until compressive stresses take place, i.e. mid-plane principal stresses must always be positive.

In this general context, we present our spatial finite element procedure and the temporal integration scheme. Then, implementation procedures for the previous chosen constitutive equations are highlighted and the treatment of the viscoelasticity is extensively investigated. Last, the relationship between membrane deformation and blowing pressure is studied.

#### 3.1. Formulation and temporal scheme

The most general form of the principle of virtual work in the Lagrangian description is given by

$$\int_{\mathcal{B}_0} \delta \mathbf{u}(t) \rho_0 \ddot{\mathbf{u}}(t) dV_0 = \int_{\partial \mathcal{B}_0} \delta \mathbf{u}(t) \mathbf{T}_0(t) dS_0 - \int_{\mathcal{B}_0} \delta \mathbf{E}(t) : \mathbf{S}(t) dV_0 \quad \forall \delta \mathbf{u}(t) \quad (11)$$

where $\mathcal{B}_0$ and $\partial \mathcal{B}_0$ are, respectively, the volume and the boundary surface of the undeformed membrane, $\rho_0$ is the constant mass density, $\ddot{\mathbf{u}}(t)$ is the acceleration vector, $\mathbf{T}_0(t)$ is the external surface force with reference to the undeformed state, $\mathbf{E}(t)$ is the Green–Lagrange strain tensor, $\mathbf{S}(t)$ is the second Piola–Kirchhoff stress tensor and finally, $\delta \mathbf{u}(t)$ is a compatible virtual displacement vector.

##### 3.1.1. Spatial discretization.
First, the previous equation is discretized in space by the finite element method. In order to perform this semi-discretization, we consider the equilibrium of the membrane at time $t$ including the acceleration dependent inertia force (here, there is no damping force). Then, the problem reduces to a system of ordinary differential equations:

$$\tilde{\mathbf{M}}\ddot{\mathbf{U}}(t) = \mathbf{F}^{\text{ext}}(t) - \mathbf{F}^{\text{int}}(t) \quad (12)$$

in which $\tilde{\mathbf{M}}$ is the mass matrix constant in time, $\ddot{\mathbf{U}}(t)$ is the nodal acceleration vector ($\mathbf{U}(t)$ stands for the vector of nodal displacements) and, $\mathbf{F}^{\text{ext}}(t)$ and $\mathbf{F}^{\text{int}}(t)$ are, respectively, the external and
the internal nodal forces. In this equation, the term on the left-hand side stands for inertia forces and fits to the term on the left-hand side of Equation (11), and the terms on the right-hand side correspond to the same position terms in Equation (11) and depend on the strain state, i.e. on the nodal displacement vector $\mathbf{U}(t)$.

Note that external forces are limited to the action of blowing pressure $P(t)$ on the deformed membrane surface. Thus, the Eulerian pressure force, $\mathbf{T}(t)$, acting on an infinitesimal deformed membrane surface $dS$, depends on the deformed geometry and is

$$ \mathbf{T}(t) = P(t) \, dS(t) \mathbf{n}(t) \qquad (13) $$

where $\mathbf{n}(t)$ is a normal vector to the deformed membrane. It implies that external forces vector $\mathbf{F}^{\text{ext}}(t)$ must be computed in the deformed configuration by replacing the previous Lagrangian external virtual work, $\int_{\partial \mathcal{B}_0} \delta \mathbf{u}(t) \mathbf{T}_0(t) \, dS_0$, in (11) by its Eulerian counterpart $\int_{\partial \mathcal{B}} \delta \mathbf{u}(t) \mathbf{T}(t) \, dS$, where $\partial \mathcal{B}$ stands for the surface of the deformed body.

In membrane context, we typically use 3-nodes triangular elements which only deform in their plane, remain flat and triangular with straight edges [7]. Computations are made by assemblage of elementary contributions. As shown later, each element contribution is calculated by exact integration in a local co-ordinates system. Therefore, no spatial numerical integration procedure (as Gauss points integration) is needed.

### 3.1.2. Time-integration scheme.
As mentioned earlier, the spatial discretization provides a system of ordinary differential equations in time. In order to solve this system, we use the explicit second-order central difference method [15]. The nodal velocity and acceleration vectors are given by

$$ \dot{\mathbf{U}}(t) = \frac{\mathbf{U}(t + \Delta t) - \mathbf{U}(t - \Delta t)}{2\Delta t} \qquad (14) $$

and

$$ \ddot{\mathbf{U}}(t) = \frac{\mathbf{U}(t + \Delta t) - 2\mathbf{U}(t) + \mathbf{U}(t - \Delta t)}{\Delta t^2} \qquad (15) $$

Note that this scheme is conditionally stable.

Using the previous temporal discretization equations, problem (12) becomes

$$ \frac{\tilde{\mathbf{M}}}{\Delta t^2} \mathbf{U}(t + \Delta t) = \mathbf{F}^{\text{ext}}(t) - \mathbf{F}^{\text{int}}(t) + \frac{\tilde{\mathbf{M}}}{\Delta t^2} [2\mathbf{U}(t) - \mathbf{U}(t - \Delta t)] \qquad (16) $$

Consequently, the nodal displacement vector at time $t + \Delta t$ can be explicitly computed using its previous values at times $t$ and $t - \Delta t$. Usually, initial displacements and velocities are assumed equal to zero.

The effectiveness of this temporal scheme is improved by the use of a diagonal mass matrix obtained by the special lumping technique [16]. The total mass of a finite element is defined by the summation of its mass matrix terms. This total mass is proportionally splitted into the diagonal terms of a new elementary mass matrix. Last, elementary matrices are assembled in $\tilde{\mathbf{M}}'$, the new diagonal mass matrix of the membrane. In that case, each term of the nodal displacement vector at time $t + \Delta t$, $U_i(t + \Delta t)$, can be obtained by

$$ U_i(t + \Delta t) = \frac{\Delta t^2}{M_{ii}'} [F_i^{\text{ext}}(t) - F_i^{\text{int}}(t)] + 2U_i(t) - U_i(t - \Delta t) \qquad (17) $$
<figure><img src="image_2.png" /><figcaption>Figure 1. Deformation of a finite element in its local co-ordinates system (R): (ABC) undeformed element, (abc) deformed element.</figcaption></figure>

where subscript ·i stands for the ith component of a vector and M'_ii is the corresponding diagonal mass term. We can clearly see that numerical operations needed to solve the dynamic problem (12) are reduced very significantly by the use of this method.

## 3.2. Numerical implementation of constitutive relations

### 3.2.1. Preliminaries.
We have seen that dynamic equations at time t are used to compute the displacement vector at time t + Δt. In order to compute the internal force vector at time t, F'int(t), we have to determine the Lagrangian deformation and stress tensors of each element, respectively, C(t) and S(t) in Equation (11). The right Cauchy-Green deformation tensor is obtained using the nodal displacement vector U(t) and the initial position of element nodes, and the second Piola-Kirchhoff stress tensor is next obtained by taking into account the constitutive relation of the material.

Then, we study each element deformation by placing the undeformed triangular element in the local co-ordinates system of the deformed element so that the two sets of local axes and their origins coincide. The origin of the local co-ordinate system is the first node of the deformed triangular element, its two first axes correspond to the element plane and the third one is normal to this plane. This situation is shown in Figure 1.

Considering the plane stress state due to membrane hypothesis, the deformation and stress matrices, $\tilde{C}(t)$ and $\tilde{S}(t)$, respectively, have the following forms:

$$ \tilde{C}(t) = \begin{bmatrix} C_{11}(t) & C_{12}(t) & 0 \\ C_{12}(t) & C_{22}(t) & 0 \\ 0 & 0 & C_{33}(t) \end{bmatrix} \qquad (18) $$

and

$$ \tilde{S}(t) = \begin{bmatrix} S_{11}(t) & S_{12}(t) & 0 \\ S_{12}(t) & S_{22}(t) & 0 \\ 0 & 0 & 0 \end{bmatrix} \qquad (19) $$

Upon use of the incompressibility constraint, the term $C_{33}(t)$ in Equation (18) can be directly computed from the other components of the deformation tensor:

$$ C_{33}(t) = \frac{1}{C_{11}(t)C_{22}(t) - C_{12}(t)^2} \qquad (20) $$
$C_{33}(t)$ is the square of the principal stretch ratio in the thickness direction, said $\lambda_3(t)$, defined by

$$ \lambda_3(t) = \sqrt{C_{33}(t)} = \frac{h(t)}{h_0} \qquad (21) $$

where $h(t)$ and $h_0$ are membrane thicknesses (uniform in each element), respectively, in the deformed and undeformed configurations. The two other principal stretch ratios, $\lambda_1(t)$ and $\lambda_2(t)$ (we assume that $\lambda_1(t) \ge \lambda_2(t)$), and the respective angles of inclination of the principal directions, $\beta_1(t)$ and $\beta_2(t)$, are calculated by diagonalizing $\tilde{C}(t)$:

$$ \lambda_1^2(t) = \frac{1}{2} [C_{11}(t) + C_{22}(t) + \sqrt{[C_{11}(t) - C_{22}(t)]^2 + 4C_{12}(t)^2}] \qquad (22) $$

$$ \tan \beta_1(t) = \frac{C_{12}(t)}{\lambda_1^2(t) - C_{22}(t)} \qquad (23) $$

and

$$ \lambda_2^2(t) = \frac{1}{2} [C_{11}(t) + C_{22}(t) - \sqrt{[C_{11}(t) - C_{22}(t)]^2 + 4C_{12}(t)^2}] \qquad (24) $$

$$ \tan \beta_2(t) = \frac{C_{12}(t)}{\lambda_2^2(t) - C_{22}(t)} \qquad (25) $$

We next examine the numerical implementation of the constitutive relations described previously.

### 3.2.2. Hyperelastic case.
Assuming the elastic behaviour of the membrane, the time $t$ is omitted in this paragraph and we next examine the numerical implementation of the Mooney's model.

Using Equations (5) and (8), the stress-strain relationship is given by

$$ S = -pC^{-1} + 2[(c_1 + c_2I_1)I - c_2C] \qquad (26) $$

in which the hydrostatic pressure is determined by the use of $S_{33}=0$:

$$ p = 2(c_1 + c_1I_1 - c_2\lambda_3^2)\lambda_3^2 \qquad (27) $$

where $\lambda_3$ is given by Equations (20) and (21).

### 3.2.3. Viscoelastic case.
Since now the current stresses depend on deformation history, we must store all physical values from initial time to current time $t$ to compute the integral in Equation (9). Then, the required computing time becomes longer and longer. In order to save both computing time and data storage, we use the recurrence formula developed by Feng [17]. With this recurrence formula, the value of the hereditary integral at current time step $t$ depends only on different variables values at the previous time step $t - \Delta t$. Assuming that the relaxation function in Equation (9) is written in terms of the exponential law:

$$ g_1(t) = g_1 e^{-t/\tau} \qquad (28) $$
in which $g_1$ is a constant and $t_R$ is the relaxation time, this formula is given by

$$ S(t) = -pC^{-1}(t) + g_0I + \frac{1}{2}e^{-\Delta t/t_R} \int_0^{t-\Delta t} g_1(t-\Delta t-\tau) \frac{\partial C(\tau)}{\partial \tau} d\tau \\ + g_1e^{-\Delta t/2t_R} [C(t) - C(t-\Delta t)] \quad (29) $$

where the hereditary integral from 0 to $t - \Delta t$ is a function of the previous time step $t - \Delta t$.

Another difficulty of such models implementation is that the kinematical variables (elementary local axes for example) change during the time step. As we cannot evaluate this change, we have to approximate it. Therefore, we assume that the principal directions of the deformation tensor remain constant during the time step $\Delta t$ between discrete times $t - \Delta t$ and $t$, and that these directions are equal to the principal directions at current time $t$. Because of large rotations involved in the present problem, this stress update procedure must be used with very small load steps. Here, because of the time-integration procedure, the time steps must be less than the critical time step to ensure convergence. This critical time-step size is a function of elements size and of the constitutive equation. It is ever very small [15]. Thus, the corresponding load steps are sufficiently small to adopt the stress update method. Rachik et al. [18] use a similar procedure for the implementation of a differential viscoelastic constitutive law in the context of blow-moulding and thermoforming processes simulation. The previous recurrence formula (29) remains valid only for the principal directions

$$ S_i(t) = -p \frac{1}{\lambda_i^2(t)} + g_0 + \frac{1}{2} e^{-\Delta t/t_R} J_i(t - \Delta t) + g_1 e^{-\Delta t/2t_R} [\lambda_i^2(t) - \lambda_i^2(t - \Delta t)], \quad i=1,3 \quad (30) $$

where $S_3(t)=0$ and where $J_i(t - \Delta t)$ are the convolution integrals

$$ J_i(t - \Delta t) = \int_0^{t - \Delta t} g_1(t - \Delta t) \frac{\partial \lambda_i^2(\tau)}{\partial \tau} d\tau \quad i = 1,3 \quad (31) $$

and the hydrostatic pressure $p$ is given by

$$ p = \lambda_3^2(t) \left\{ g_0 + \frac{1}{2} e^{-\Delta t/t_R} J_3(t - \Delta t) + g_1 e^{-\Delta t/2t_R} [\lambda_3^2(t) - \lambda_3^2(t - \Delta t)] \right\} \quad (32) $$

in which $\lambda_3(t)$ is deduced from Equations (20) and (21).

For each finite element of the mesh, we must store five principal values relative to the previous time $t - \Delta t$: two principal components of $C$, $\lambda_1^2$ and $\lambda_2^2$ ($\lambda_3^2$ is computed using the incompressibility property), and the three principal viscous hereditary integrals, $J_i(t - \Delta t)$ for $i = 1,3$.

Now, the numerical implementation of Christensen's model can be summarized as follows:

1. as $U(t)$ is known from the previous deformed state, $C(t)$ is computed in element local axes,

2. by diagonalizing $C(t)$, we obtain its plane principal stretch ratios $\lambda_1(t)$, $\lambda_2(t)$ and the corresponding principal directions, $\beta_1(t)$ and $\beta_2(t)$,

3. using incompressibility and plane stress state assumptions, data stored at $t - \Delta t$ and Equations (30), (31) and (32), principal stresses, $S_1(t)$ and $S_2(t)$, are computed,

4. this allows us to determine the stress tensor $S(t)$ in the local axes by a simple change of basis.
<figure><img src="image_3.png" /><figcaption>Figure 2. Notations for the pressure evolution calculus: (a) initially non-plane membrane; (b) initially plane membrane.</figcaption></figure>

### 3.3. Pressure evolution calculus

As mentioned earlier, the external forces are only due to the inflating pressure acting on the deformed membrane. In most of the cases, numerical attempts to simulate membrane inflations assume that the pressure-time history is an external loading data (i.e. is independent of the membrane deformation state).

In most practical applications, only gas flow rate can be mastered. Thus, the pressure inside the volume limited by the membrane highly depends on the geometry of the membrane: as the membrane inflates, the pressure changes following the thermodynamical law of the inflating gas. Therefore, it is more realistic to impose a gas flow rate function and to compute the corresponding pressure evolution.

In this paper, we assume that the gas used to inflate the membrane obey the Perfect Gas Law. Denoting $P_0$, $V_0$ and $m_0$ the initial pressure around the membrane, the volume delimited by the undeformed membrane and the initial gas mass contained inside this membrane, and $m^{mol}$ the gas molecular mass, $R$ the perfect gas universal constant, $T^{gas}$ the gas temperature assumed to be constant during inflation and $V_{init}$ an arbitrary reservoir volume equal to zero in the case of an initial non-plane membrane ($V_0 \neq 0$) and non-zero (arbitrary or experimentally defined) in the case of a plane membrane ($V_0 = 0$), the initial thermodynamic state of the gas is written as

$$P_0(V_0 + V_{init}) = \frac{m_0}{m^{gas}} RT^{gas} \quad (33)$$

This equation can be used to determine the gas mass initially contained in the membrane (i.e. $m_0$). Figure 2 shows the previous notations for the two different cases.

At time $t$, the current internal pressure $P(t)$ is given by

$$P(t) [V^{init} + V(t)] = \frac{m_0 + m(t)}{m^{mol}} RT^{gas} \quad (34)$$

where $V(t)$ is the volume limited by the deformed membrane, $m(t)$ is the extra gas mass added between initial time and $t$. Thus, the pressure difference $P(t) - P_0$, denoted by $P^*(t)$, takes the following form:

$$P^*(t) = \frac{m(t)RT^{gas} + P_0[V(t) - V_0]}{V_{init} + V(t)} \quad (35)$$

In this equation, the gas flow rate $m(t)$ is imposed.
The numerical implementation of the pressure evolution calculus is very simple. Initially, we assume that the external equilibrium pressure $P_0$ and the reservoir volume $V_{init}$ are well known and that $V_0$ can be computed from the undeformed configuration. At the current time $t$, the internal volume delimited by the deformed membrane $V(t)$ is computed using displacements. Equation (35) allows us to directly obtain the uniform loading pressure $P^*(t)$ acting on the finite elements mesh. This pressure value is used to solve the problem at time $t$, i.e. to compute the next displacement vector $U(t + \Delta t)$.

## 4. NUMERICAL EXAMPLES

In order to illustrate the capability of our work, three different examples are treated and results are compared with analytical or experimental results.

### 4.1. Instability of spherical membranes under a pressure step

The relationship between the inflating pressure and the geometry of inflated rubber-like membranes was extensively analysed in the past, experimentally [19] and theoretically [20]. The instability phenomena associated with these inflations were examined by Khayat et al. [21] in the case of the static inflation of non-linear elastic cylinders.

More recently, Verron et al. [22] showed the existence of instabilities during the dynamic inflation of membranes. They extensively investigated the special case of spherical Mooney membranes under constant pressure steps.

We study this case using our numerical procedure. Both hyperelastic and non-linear viscoelastic membrane inflations are next examined. The corresponding finite element results obtained on a complete spherical mesh (without any symmetry condition) of 2000 finite elements are compared with semi-analytical results of the one-dimensional reduced problem (obtained by assuming the spherical symmetry hypothesis) solved by a fifth- or sixth-order Runge-Kutta method (subroutine IVPRK from IMSL library [23]). In these examples, the membrane is supposed to be inflated by a constant pressure step. Numerically, the corresponding pressure function is an exponential one:

$$P^*(t) = P_{\infty}^*(1 - e^{-100t}) \qquad (36)$$

Therefore, the prescribed pressure step value $P_{\infty}^*$ is reached very early and in the following examples the dependence of $P^*$ on $t$ is omitted. The time step is set to $10^{-6}$ s due to elements size. Therefore, the loading steps remain very small during computations and the use of the stress update procedure is justified.

#### 4.1.1. Mooney's membrane.
The resolution of the corresponding semi-analytical problem is detailed in Reference [22]. Consider a spherical hyperelastic membrane. The initial radius, thickness and mass density are noted $R_0$, $h_0$ and $\rho_0$, respectively. The material is assumed to obey the Mooney's model and the constitutive parameters are $c_1$ and $c_2$. Moreover, we define the classical non-dimensional Mooney's constant $\alpha$ by

$$\alpha = \frac{c_2}{c_1} \qquad (37)$$
Denoting the time $t$ and the loading pressure step $P^*$, we define the corresponding reduced time, $\tau$, and pressure step, $p^*$, by

$$ \tau = t \frac{2}{R_0} \sqrt{\frac{c_1}{h_0}} \qquad (38) $$

and

$$ p^* = P^* \frac{R_0}{4c_1h_0} \qquad (39) $$

The non-dimensional governing equation of the membrane inflation is a second-order differential equation:

$$ \ddot{\lambda} = p^* \dot{\lambda}^2 + \left( \frac{1}{\dot{\lambda}^5} - \dot{\lambda} \right) (1 + \alpha \dot{\lambda}^2) \qquad (40) $$

where $\lambda$ is the circumferential principal extension defined as the ratio of the deformed radius to the undeformed radius. We only consider the case of the membrane initially in equilibrium, i.e. subjected to the following initial conditions:

$$ \lambda(\tau=0) = 1 \quad \text{and} \quad \dot{\lambda}(\tau=0) = 0 \qquad (41) $$

Depending on the values of material parameter $\alpha$ and the pressure step $p^*$, three different inflation modes can take place:

* the membrane may oscillate around equilibrium points,

* it may tend to inflate indefinitely,

* it may tend to reach an unstable equilibrium point, then it neither oscillates nor inflates any further.

The comparison between one-dimensional Runge-Kutta results and finite element results for each previous case are now examined. Figure 3 shows the results for three different values of $\alpha$. Figure 3(a) is relative to the case where $\alpha=0$ and exhibits two different dynamic behaviours. The numerical results closely reproduce both oscillatory behaviour ($p^*=0.5$) and unstable growing stage ($p^*=0.6$). Now, consider the case where $\alpha=0.1$. The corresponding dynamic behaviours are shown in Figure 3(b). For this value of the material parameter, there is an unstable equilibrium point, i.e. a saddle point, which corresponds to $p^* \approx 0.687$. Far from this point, numerical results are highly similar to the analytical results ($p^*=0.8$): the loss of linearity of oscillations is well reproduced. But for values of $p^*$ near this saddle point ($p^*=0.687$ and $0.7$), the finite element results can only approach qualitatively the real behaviour: the temporal evolution of the deformed radius obtained by the finite element method differs from the analytical solution. In such cases, the membrane tends to stabilize at the saddle point and thus, the period of oscillations tends asymptotically to infinity. That is the reason why the numerical results are highly sensitive to the pressure step value. For $\alpha=0.25$, the membrane oscillates whatever the value of the pressure step and the 3D calculations closely simulate the behaviour (see Figure 3(c) for $p^*=0.7$ and $1.2$).

Moreover, for the three previous values of $\alpha$, we analytically obtain the curves of oscillations period versus pressure step and, for some values of the pressure step, 3D simulations are made. The corresponding analytical curves and numerical points are shown in Figure 4. In this figure, the two vertical lines correspond to the limit pressure values for $\alpha=0$ and $0.1$. In the former case,
<figure><img src="image_4.png" /><figcaption>Figure 3. Dynamic inflation of a Mooney spherical membrane, (—) one-dimensional Runge-Kutta, (○) finite elements: (a) $\alpha = 0$, (b) $\alpha = 0.1$, (c) $\alpha = 0.25$.</figcaption></figure>

the membrane cannot withstand pressure greater than $p^*=0.556$, and in the latter $p^*=0.687$ corresponds to the saddle point and the membrane can afford greater inflating pressure values.

These results show that our finite element procedure is able to simulate the dynamic behaviour of inflated membranes even in the vicinity of the saddle points (near the vertical line in the figure) and to reproduce the physical unstable phenomena, which take place during the dynamic inflation of hyperelastic membranes.

### 4.1.2. Christensen's membrane.
In this section, we study the dynamic inflation of a non-linear viscoelastic spherical membrane, which obey Christensen's model. The one-dimensional analytical equation of this problem is obtained by the same manner as in the previous hyperelastic case [24].
<figure><img src="image_5.png" /><figcaption>Figure 4. Influence of the pressure step on the period of motion for three values of $\alpha$: (—) one-dimensional Runge–Kutta, (⊙) finite elements.</figcaption></figure>

The initial spherical bubble radius, thickness and mass density are noted $R_0$, $h_0$ and $\rho_0$, respectively. Assuming that the constitutive parameters are $g_0$, $g_1$ and $t_R$ as defined in Section 3.2.3, we consider again reduced time and pressure defined by

$$ \tau = t \frac{\sqrt{2}}{R_0} \sqrt{\frac{g_0}{\rho_0}} \qquad (42) $$

and

$$ p^* = P^* \frac{R_0}{2g_0h_0} \qquad (43) $$

Moreover, we consider a reduced relaxation time $\tau_R$ and a reduced relaxation function $\alpha(\tau)$:

$$ \tau_R = t_R \frac{\sqrt{2}}{R_0} \sqrt{\frac{g_0}{\rho_0}} \qquad (44) $$

and

$$ \alpha(\tau) = \frac{g_1}{g_0} e^{-\tau/\tau_R} \qquad (45) $$

Assuming the spherical symmetry condition, the one-dimensional governing differential equation of the viscoelastic problem is written in terms of the principal circumferential extension, said $\lambda$:

$$ \begin{aligned} \ddot{\lambda}(\tau) = & p^* \lambda^2(\tau) + \left[ \frac{1}{\lambda^5(\tau)} - \lambda(\tau) \right] - \frac{\lambda(\tau)}{2} \int_0^\tau \alpha(\tau - u) \frac{d\lambda^2(u)}{du} du \\ & + \frac{1}{2\lambda^5(\tau)} \int_0^\tau \alpha(\tau - u) \frac{d\lambda^{-4}(u)}{du} du \end{aligned} \qquad (46) $$
<figure><img src="image_6.png" /><figcaption>Figure 5. Dynamic inflation of a Christensen spherical membrane: (—) one-dimensional Runge-Kutta, (○) finite elements: (a) $p^*=0.5$; (b) $p^*=1.0$.</figcaption></figure>

The convolution integrals are computed using the Feng's recurrence formula (31) and this differential equation is solved by a fifth- or sixth-order Runge-Kutta method.

Taking into account the complexity of a parametric study, we choose to fix two parameter values:

$$ \alpha = 0.2 \quad \text{and} \quad \tau_R = 1 \qquad (47) $$

Then, two different behaviours are highlighted depending on the value of the pressure step:

* for sufficiently small values of $p^*$, the membrane oscillates and oscillations tend to decrease with time in amplitude,

* for greater values of $p^*$, the membrane exhibits an unstable behaviour for which its radius grows to infinity with small oscillations.

We next examine the comparison between one-dimensional analytical and finite element results. The two previous cases are illustrated in Figure 5: Figures 5(a) and 5(b) show, respectively, a stable oscillatory motion ($p^* = 0.5$) and an unstable continuously growing motion ($p^* = 1$). In this figure, the continuous curves correspond to the semi-analytical results and the circles represent numerical calculations.

For these two cases, the numerical simulations are in good agreement with analytical results. However, we note that both amplitude and period of motion in the oscillatory case and only amplitude in the unstable case, are underestimated. These small differences between numerical and analytical results are due to the approximations made for the implementation of the viscoelastic constitutive relation in the program (see Section 3.2.3).
<figure><img src="image_7.png" /><figcaption>Figure 6. Dynamic inflation of a Mooney spherical membrane by a constant gas flow rate for (a) $\alpha=0$ and (b) $\alpha=0.25$: (—) one-dimensional Runge-Kutta, (◦) finite elements, (---) static solution.</figcaption></figure>

## 4.2. Inflation of a Mooney spherical membrane by a constant gas flow rate

The simulation of the inflation of a hyperelastic spherical membrane by a constant gas flow rate is carried out in this example. The geometry of the initial membrane and the definition of the reduced parameters are the same as in the first simulations. The time step size is also $10^{-6}$ s.

The corresponding one-dimensional governing equation is given by [24]

$$ \ddot{\lambda} = p_0 \left[ \frac{1}{\lambda} (1 + q\tau) \right] + \left( \frac{1}{\lambda^5} - \lambda \right) (1 + \alpha\lambda^2) \quad (48) $$

where $p_0$ is the initial reduced equilibrium pressure inside the bubble and $q$ a reduced constant gas flow rate.

In this problem, initial pressure and gas flow rate are both equal to 1 and numerical simulations are performed for two Mooney's reduced parameter $\alpha$. The pressure is computed using the method developed in Section 3.3. The results provide the pressure inside the bubble versus the circumferential extension. Comparisons between analytical and numerical results are shown in Figures 6(a) and 6(b), respectively, for $\alpha=0$ and $0.25$.

Firstly, we note that numerical solutions closely fit to analytical results in the two cases. Secondly, the dynamic curves oscillate around the corresponding static curves (dotted-dash lines in Figures 6) obtained by setting $\ddot{\lambda}=0$ in Equation (48). The dynamic numerical curve relative to the case where $\alpha=0$ (Figure 6(a)) oscillates around the static curve which exhibits a maximum pressure point beyond which the membrane grows as the pressure decreases (called snap-through behaviour). Classically, it is difficult to capture equilibrium beyond the maximum pressure value and special numerical algorithms have to be developed [25, 26]. In this case, we show that if we consider the thermodynamical behaviour of the gas, we are able to directly obtain the real motion and pressure evolution of the membrane without any convergence problem.
<table><caption>Table I. Initial geometrical data for elliptical membranes.</caption><thead><tr><th></th><th>Membrane 1</th><th>Membrane 2</th></tr></thead><tbody><tr><td>a</td><td>70.7 mm</td><td>50 mm</td></tr><tr><td>b</td><td>141.4 mm</td><td>200 mm</td></tr><tr><td>e</td><td>2</td><td>4</td></tr><tr><td>h<sub>0</sub></td><td>0.23 mm</td><td>0.23 mm</td></tr></tbody></table>

### 4.3. Inflation of elliptical membranes

To test the ability of our finite element model in 3D problems, the dynamic inflation of the initially plane elliptical neo-Hookean membranes has been analysed. This problem is based on the work of Charrier et al. [27] which performed experiments and quasi-static calculations on such latex rubber membranes.

We consider two membranes defined by the lengths of their semi-minor axis, a, semi-major axis, b, and by a reduced aspect ratio, $e=b/a$. The undeformed thickness is assumed to be uniform and is denoted by $h_0$. These geometrical data are presented in Table I for the two membranes. The material is a natural latex rubber, which is described by a neo-Hookean model with $c_1 = 141.5$ kPa. The mass density is 1200 kg/m³. The constant gas flow rate and the time-step size are, respectively, set to $15 \times 10^4$ mm³/s and $6.6 \times 10^{-5}$ s. In order to compare our results with Charrier's experiments, we define the adimensional pressure inside the bubble, $\bar{P}$, and the adimensional apex height, $\bar{z}_a$, by

$$ \bar{P} = \frac{Pa}{2c_1h_0} \qquad (49) $$

and

$$ \bar{z}_a = \frac{z_a}{a} \qquad (50) $$

where $P$ is the pressure inside the bubble and $z_a$ is the apex height of the inflated membrane.

The symmetric quadrants of the two membranes are, respectively, meshed with 331 nodes and 584 elements for the $e=2$ membrane, and 343 nodes and 591 elements for the $e=4$ membrane. They are submitted to a constant gas flow rate, then the pressure inside the bubble is calculated by the program. Figures 7 and 8 show four inflation stages for the $e=2$ and 4 membranes, respectively. It is to note that in the latter case, the membrane undergoes an instability shown in Figure 8(d). This instability is called a bulge and is observed in experiments [27]. It might be due to the presence of compressive stresses in the membrane.

Curves of Figure 9 presents the pressure evolution versus the membrane apex height for the two membranes. The pressure reaches a maximum value after which it decreases as the inflation goes on. The numerical curves match experimental points reasonably with a more important difference for the $e=4$ membrane, due to the occurrence of the bulge.

Moreover, we compare our numerical results with the measurements of the principal stretch ratios at the apex. Figure 10 presents these curves for $e=2$ (curve (a)) and 4 (curve (b)). In this figure, we can see that our predictions of principal stretch ratios are better in the case $e=4$, i.e. the bulging case. Taking into account the previous comparisons on the pressure evolution, we should
<figure><img src="image_8.png" /><figcaption>Figure 7. Four inflation stages of an elliptical membrane with $e=2$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8. Four inflation stages of an elliptical membrane with $e=4$.</figcaption></figure>

have a better fit for the $e=2$ membrane, which inflates regularly. Following Charrier *et al.*, these differences are due to experimental difficulties in measuring the principal stretch ratios at the apex. The use of a neo-Hookean constitutive equation to describe latex behaviour is also questionable.

## 5. CONCLUSION

A three-dimensional finite element formulation of the dynamic inflation of rubber-like membranes has been presented in this paper. The material sheets are assumed to be hyperelastic or non-linear
<figure><img src="image_10.png" /><figcaption>Figure 9. Reduced inflation pressure inside the bubble versus reduced apex height of elliptical membranes defined by: (a) $e=2$ and (b) $e=4$: (—) finite elements, (○) experiments [27].</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 10. Apex principal stretch ratios versus reduced apex height of elliptical membranes defined by: (a) $e=2$ and (b) $e=4$: (—) finite elements for $\lambda_1$, (○) experiments for $\lambda_1$ [27], (· · ·) finite elements for $\lambda_2$ and (●) experiments for $\lambda_2$ [27].</figcaption></figure>

viscoelastic. Then, the hyperelastic Mooney's model and the integral viscoelastic Christensen's model have been described and implemented. The numerical implementation of the viscoelastic behaviour is made by the use of the recurrence formula of Feng [17] to calculate the convolution integrals with minimum data storage. A membrane formulation is adopted, thus the incompressibility constraint is imposed in a trivial manner. The explicit finite difference scheme for time integration was used. Moreover, in order to be more realistic, we consider an external loading by a gas flow rate rather than by a pressure-time history: as the pressure is highly dependent on the deformation state geometry (volume inside the bubble), the real external independent loading is the gas flow rate and the pressure has to be recalculated inside the membrane, taking into account the deformed geometry.

Finally, different comparisons between numerical results, and analytical and experimental data show that our formulation has the capability to reproduce non-linear dynamic large strains of
inflated membrane and to predict instability phenomena as unstable growing behaviour or bulge occurrence.

The present study leaves some issues of fundamental importance. The formulation does not handle wrinkle phenomenon which arises as the membrane is compressed. In order to predict the post-bulging behaviour, it will be necessary to take into account this phenomenon in the future.

REFERENCES

1. Green AE, Adkins JE. *Large Elastic Deformations*. The Clarendon Press: Oxford, 1960.

2. Hart-Smith LJ, Crisp JDC. Large elastic deformations of thin rubber membranes. *International Journal of Engineering Science* 1967; **5**:1–24.

3. Klingbeil WW, Shield RT. Some numerical investigations on empirical strain energy functions in the large axisymmetric extensions of rubber membranes. *Zeitschrift für Angewandte Mathemathik und Physik* 1964; **15**:608–629.

4. Oden JT, Sato T. Finite strains and displacements of elastic membranes by the finite element method. *International Journal of Solids and Structures* 1967; **3**:471–488.

5. Wineman A. On axisymmetric deformations of nonlinear viscoelastic membranes. *Journal of Non-Newtonian Fluid Mechanics* 1978; **4**:249–260.

6. Feng WW. Viscoelastic behavior of elastomeric membranes. *Journal of Applied Mechanics* 1992; **59**:S29–S34.

7. Shrivastava S, Tang J. Large deformation finite element analysis of non-linear viscoelastic membranes with reference to thermoforming. *Journal of Strain Analysis* 1993; **28**(1):31–51.

8. Jenkins CH, Leonard JW. Nonlinear dynamic response of membranes: state of the art. *Applied Mechanic Review* 1991; **44**(7):319–328.

9. Jenkins CH. Nonlinear dynamic response of membranes: state of the art—update. *Applied Mechanic Review* 1996; **49**(10):S41–S48.

10. Truesdell C. *The Elements of Continuum Mechanics*. Springer: New York, 1966.

11. Mooney M. A theory of large elastic deformation. *Journal of Applied Physics* 1940; **11**:582–592.

12. Treloar LRG. The mechanics of rubber elasticity. *Proceedings of the Royal Society of London* 1976; **A351**:301–330.

13. Ward IM, Hadley DW. *An Introduction to the Mechanical Properties of Solid Polymers*. Wiley: New York, 1993.

14. Christensen RM. *Theory of Viscoelasticity. An Introduction* (2nd edn). Academic Press: New York, 1982.

15. Dokainish MA, Subbaraj K. A survey of direct time-integration methods in computational structural dynamics. I. Explicit methods. *Computers and Structures* 1989; **32**(6):1371–1386.

16. Zienkiewicz OC, Taylor RL. *The Finite Element Method. Basic Formulation and Linear Problems*, vol. I. MacGraw-Hill: New York, 1989.

17. Feng WW. A recurrence formula for viscoelastic constitutive equations. *Proceedings of the International Conference on Computational Mechanics*, Tokyo, vol. IV, May 25–29, 1986; 77–81.

18. Rachik M, Roelandt JM, Batoz JL. Numerical simulation of blow-moulding and thermoforming of plastics. *Revue Européenne des Éléments Finis* 1994; **3**(2):187–210 (in French).

19. Alexander H. The tensile instability of initially spherical balloons. *International Journal of Engineering Science* 1971; **9**:151–162.

20. Beatty MF. Topics in finite elasticity: hyperelasticity of rubber, elastomers, and biological tissues—with examples. *Applied Mechanic Review* 1987; **40**(12):1699–1734.

21. Khayat RE, Derdouri A, Garcia-Réjon A. Inflation of an elastic cylindrical membrane: non-linear deformation and instability. *International Journal of Solids and Structures* 1992; **29**(1):69–87.

22. Verron E, Khayat RE, Derdouri A, Peseux B. Dynamic inflation of hyperelastic spherical membranes. *Journal of Rheology* 1999; **43**:1083–1097.

23. IMSL Inc. *IMSL Math Library. FORTRAN Subroutines for Mathematical Applications*. User's Manual Version 1.0, 1987.

24. Verron E. Experimental and numerical contribution to the blow-moulding and thermoforming processes. *Ph.D. Thesis*, Ecole Centrale de Nantes, Nantes, France, 1997 (in French).

25. De Souza Neto EA, Peric D, Owen DRJ. Finite elasticity in spatial description: linearization aspects with 3-D membrane applications. *International Journal for Numerical Methods in Engineering* 1995; **38**:3365–3381.

26. Shi J, Moita GF. The post-critical analysis of axisymmetric hyper-elastic membranes by the finite element method. *Computer Methods in Applied Mechanics and Engineering* 1996; **135**:265–281.

27. Charrier JM, Shrivastava S, Wu R. Free and constrained inflation of elastic membranes in relation to thermo-forming—non-axisymmetric problems. *Journal of Strain Analysis* 1989; **24**(2):55–73.