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

HAL
open science

# B-splines modeling of free inflation of axisymmetric rubberlike membranes

Gilles Marckmann, Erwan Verron, Bernard Peseux

► To cite this version:

Gilles Marckmann, Erwan Verron, Bernard Peseux. B-splines modeling of free inflation of axisymmetric rubberlike membranes. European Congress on Computational Methods in Applied Sciences and Engineering (ECCOMAS 2000), Sep 2000, Barcelone, Spain. hal-01007931

HAL Id: hal-01007931

https://hal.science/hal-01007931

Submitted on 7 Oct 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
B-SPLINES MODELING OF FREE INFLATION OF
CYLINDRICAL RUBBERLIKE MEMBRANES

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

Laboratoire Mécanique et Matériaux

Division Structures

École Centrale de Nantes
BP 92101

44321 Nantes cedex 3, France

e-mail: erwan.verron@ec-nantes.fr

web page: http://www.ec-nantes.fr/Structures/

†Corresponding author

**Key words:** B-spline interpolation, Membrane inflation, Hyperelasticity, Non-linear numerical method.

**Abstract.** Free inflation of cylindrical rubberlike membranes is examined in this study. The material is assumed to obey the classical Mooney-Rivlin constitutive equation. The membrane profile is interpolated by cubic B-spline functions with natural end boundary conditions. The macro element formulation is detailed. The non-linear system obtained is solved by the Newton-Raphson algorithm associated with an arc-length continuation method. Numerical illustrations are given in order to demonstrate the method efficiency.
# 1 INTRODUCTION

Free inflation of hyperelastic cylindrical membranes has been extensively examined in the past. For infinite cylinders, some authors developed analytical solutions [1, 2]. In the case of finite cylinders, numerical procedures have to be used: governing equations reduce to a two-point boundary value problem of a system of ordinary differential equations [3]. Some authors used initial-value problem type procedures coupled with shooting algorithms [4]. Such methods reveal limitations and often fail. More recently, the finite element method with iterative solvers such as the Newton-Raphson algorithm were implemented [5]. Moreover, it is well-known that inflation of non-linear elastic membranes exhibits unstable behaviour due to the presence of both limit and bifurcation points [6, 7]. Thus special path-following techniques have to be developed to detect these points and to predict post-bifurcation behaviour [8, 9].

In the recent past, some structural mechanical studies were carried out with the help of spline functions. Due to their smoothness and continuity properties, spline functions offer great advantages compared with classical interpolation methods [10]. The two first works which report on the use of spline interpolation are due to Shik [11] and Cheung et al. [12]. Later, Gupta et al. [13], and Vermeulen and Heppler [14] analyzed the static behaviour of shells with cubic B-splines. More recently, dynamics characteristics of beams [15] and plates [16] are examined.

In the general context of spline functions in structural mechanics, the present paper reports a first attempt to model large deformations of hyperelastic membranes by B-splines interpolation. In the next section, the problem is presented. Governing equations of the inflation problem are briefly recalled. The rubberlike behaviour of the membrane is described by the classical Mooney-Rivlin model. Section 3 presents the B-splines definition and highlights the interpolation of the membrane. Some differentiation properties are proposed. The two next sections focus on the numerical procedure. In section 4 the interpolation is used to build the out of balance force vector and the numerical method which coupled the Newton-Raphson algorithm with the arc-length continuation method is examined in section 5. Last, the method efficiency is illustrated with some numerical examples.

# 2 GOVERNING EQUATIONS

Consider the axisymmetric deformations of a cylindrical membrane of non-uniform radius and uniform thickness in the undeformed state. The membrane is composed of homogeneous, isotropic, incompressible elastic material and undergoes large strains. By definition the membrane offers no moment or transverse shear resistance and the thickness is considered much less than any radius of curvature. In this context the geometry is described by the cylindrical coordinates systems $(r_0, \theta_0, z_0)$ and $(r, \theta, z)$ in the undeformed and deformed states respectively. Due to the symmetry there is no dependence on $\theta_0$, i.e.: $\theta = \theta_0$. Undeformed coordinates $r_0$ and $z_0$ are imposed functions of arc-length
coordinate s and a particle originally at $r_0(s), z_0(s)$ is positioned at $(r(s), z(s))$ in the deformed configuration. Respective thicknesses are denoted $h_0$ and $h(s)$. The membrane is subjected to an imposed internal inflation pressure $p$.

Due to axial symmetry, principal stretch directions remain constant and coincide with the meridian, the circumference and the normal to the surface of the deformed membrane [3]. The corresponding principal stretch ratios are given by:

$$ \lambda_1 = \sqrt{\frac{r'^2 + z'^2}{r_0'^2 + z_0'^2}} , \quad \lambda_2 = \frac{r}{r_0} , \quad \lambda_3 = \frac{h}{h_0} \qquad (1) $$

in which the prime denotes differentiation according to $s$. The outward normal to the surface of the deformed membrane $n$ is:

$$ \mathbf{n} = \left\{ \begin{array}{c} \frac{z'}{\sqrt{r'^2 + z'^2}} \\ -\frac{r'}{\sqrt{r'^2 + z'^2}} \end{array} \right\} \qquad (2) $$

Considering that the reference configuration is well-known, we formulate the equilibrium equation in a variational form [8] with reference to the undeformed configuration for internal forces (Lagrangian approach). Noting $\mathbf{u}$ the displacement vector and $\delta\mathbf{u}$ a virtual displacement vector, and using the Principle of Virtual Work the residual $R(\mathbf{u}, \delta\mathbf{u}, p)$ is written as follow:

$$ R(\mathbf{u}, \delta\mathbf{u}, p) = \int_{\mathcal{B}_0} \delta\mathbf{E} : \mathbf{S} \, dV - \int_{\partial\mathcal{B}} \delta\mathbf{u} \cdot p \, dS \qquad (3) $$

in which $\mathcal{B}_0$ is the undeformed volume, $\partial\mathcal{B}$ is the deformed surface, and $\mathbf{E}$ and $\mathbf{S}$ are respectively the Green-Lagrange strain tensor and the second Piola-Kirchhoff stress tensor. The Green-Lagrange principal strains can be written in terms of principal stretch ratios:

$$ E_i = \frac{\lambda_i^2 - 1}{2} \qquad i = 1, 3 \qquad (4) $$

The first term on the right-hand side of (3) is the virtual work of internal loads and the second term stands for the virtual work of deformation dependent external loads. It is defined on the current deformed configuration. The equilibrium equation yields:

$$ \forall \delta\mathbf{u} \qquad R(\mathbf{u}, \delta\mathbf{u}, p) = 0 \qquad (5) $$

In this paper we study the inflation of rubberlike membranes. The corresponding material behaviour is considered hyperelastic and we examine the highly used Mooney-Rivlin model [17]. The corresponding strain energy function $W$ is expressed as:

$$ W = C [(I_1 - 3) + \alpha (I_2 - 3)] \qquad (6) $$
where $C$ and $\alpha$ are material parameters, and $I_1$ and $I_2$ are the first and second strain invariants respectively. In terms of principal stretch ratios and using the incompressibility assumption $\lambda_1\lambda_2\lambda_3 = 1$, strain invariants can be written as:

$$I_1 = \lambda_1^2 + \lambda_2^2 + \frac{1}{\lambda_1^2 \lambda_2^2} \qquad (7)$$

$$I_2 = \frac{1}{\lambda_1^2} + \frac{1}{\lambda_2^2} + \lambda_1^2 \lambda_2^2 \qquad (8)$$

As the material is isotropic $S$ and $E$ are coaxial tensors: principal stress directions coincide with principal strain directions. Principal stresses are denoted $S_i$ with $i = 1, 3$ and are functions of principal stretch ratios [18]:

$$S_i = -p_h \frac{1}{\lambda_i^2} + 2C \left(1 - \alpha \frac{1}{\lambda_i^4}\right) \qquad (9)$$

where $p_h$ is the hydrostatic pressure due to incompressibility. As the membrane is in a plane stress state, the hydrostatic pressure can be eliminated by using:

$$S_3 = 0 \qquad (10)$$

and the two first principal stresses becomes:

$$S_1 = 2C \left[ \left( 1 - \frac{1}{\lambda_1^4 \lambda_2^2} \right) - \alpha \left( \frac{1}{\lambda_1^2} - \lambda_2^2 \right) \right] \qquad (11)$$

$$S_2 = 2C \left[ \left( 1 - \frac{1}{\lambda_1^2 \lambda_2^4} \right) - \alpha \left( \frac{1}{\lambda_2^2} - \lambda_1^2 \right) \right] \qquad (12)$$

Finally, using geometrical definitions (1), the virtual work difference (3) can be casted in the following form:

$$R(\mathbf{u}, \delta\mathbf{u}, p) = \int_{0}^{l_0} 2\pi (\delta E_1 S_1 + \delta E_2 S_2) r_0 h_0 ds - \int_{0}^{l} 2\pi p (\delta u_r z' - \delta u_z r') \frac{r}{\sqrt{r'^2 + z'^2}} ds \qquad (13)$$

in which $l_0$ an $l$ represent respectively the undeformed and deformed lengths of the cylinder, $\delta u_r$ and $\delta u_z$ are virtual displacements in radial and axial directions.

# 3 B-SPLINES INTERPOLATION

## 3.1 Definition

A spline of degree $N$ is a piecewise $C^{N-1}$ continuous function composed of polynomial segments defined on each interval of a set of $n + 1$ knots $\xi_0, \dots, \xi_n$. In the present paper cubic B(asic)-splines are considered ($N = 3$).
Thus a spline $S(\xi)$ is a linear combination of B-spline functions:

$$S(\xi) = \sum_{i=-1}^{n+1} B_i(\xi) V_i \quad (14)$$

in which $(V_i)_{i=-1,n+1}$ are the spline parameters and $B_i(\xi)$ are $n+2$ polynomial functions build on the set of knots and defined by:

$$B_i(\xi) = \left\{ \begin{array}{ll} \frac{(\xi - \xi_{i-2})^3}{(\xi_{i+1} - \xi_{i-2})(\xi_i - \xi_{i-2})(\xi_{i-1} - \xi_{i-2})} & \xi_{i-2} \le \xi < \xi_{i-1} \\[2ex] \frac{(\xi - \xi_{i-2})^2(\xi_i - \xi)}{(\xi_{i+1} - \xi_{i-2})(\xi_i - \xi_{i-2})(\xi_i - \xi_{i-1})} + \frac{(\xi_{i+1} - \xi)(\xi - \xi_{i-1})(\xi - \xi_{i-2})}{(\xi_{i+1} - \xi_{i-1})(\xi_i - \xi_{i-1})(\xi_{i+1} - \xi_{i-2})} & \xi_{i-1} \le \xi < \xi_i \\[2ex] + \frac{(\xi - \xi_{i-1})^2(\xi_{i+2} - \xi)}{(\xi_{i+1} - \xi_{i-2})(\xi_i - \xi_{i-1})(\xi_{i+2} - \xi_{i-1})} & \xi_i \le \xi < \xi_{i+1} \\[2ex] \frac{(\xi_{i+1} - \xi)^2(\xi - \xi_{i-2})}{(\xi_{i+1} - \xi_{i-2})(\xi_{i+1} - \xi_{i-1})(\xi_{i+1} - \xi_i)} + \frac{(\xi_{i+2} - \xi)(\xi_{i+1} - \xi)(\xi - \xi_{i-1})}{(\xi_{i+2} - \xi_{i-1})(\xi_{i+1} - \xi_{i-1})(\xi_{i+1} - \xi_i)} & \xi_i \le \xi < \xi_{i+1} \\[2ex] + \frac{(\xi_{i+2} - \xi)^2(\xi - \xi_i)}{(\xi_{i+2} - \xi_{i-1})(\xi_{i+2} - \xi_i)(\xi_{i+1} - \xi_i)} & \xi_{i+1} \le \xi < \xi_{i+2} \\[2ex] \frac{(\xi_{i+2} - \xi)^3}{(\xi_{i+2} - \xi_{i-1})(\xi_{i+2} - \xi_i)(\xi_{i+2} - \xi_{i+1})} & \xi_{i+1} \le \xi < \xi_{i+2} \end{array} \right. \quad (15)$$

In this definition, if $j \le 0$ then $\xi_j = \xi_0$ and if $j \ge n$ then $\xi_j = \xi_n$. In these special cases, some denominators can be equal to 0 and we adopt the convention $0/0 = 0$ to keep valid the definition. Figure 1 shows the basic functions (15) for a set of ten equidistant knots.

<figure><img src="image_3.png" /><figcaption>Figure 1: The twelve basic polynomial functions built on a set of ten knots</figcaption></figure>

## 3.2 Membrane interpolation

The membrane is discretized in $m$ subintervals. Boundaries of these intervals are called nodes (similar to the classical finite element method) and are noted $(N^j)_{j=0,m}$.
Coordinates of the node $N^j$ are noted $r_0^j, z_0^j$. The coordinates of the undeformed membrane material points are interpolated on cubic B-splines basis:

$$r_0(\xi) = \sum_{i=-1}^{n+1} B_i(\xi) \alpha_i \qquad (16)$$

$$z_0(\xi) = \sum_{i=-1}^{n+1} B_i(\xi) \beta_i \qquad (17)$$

where $(\alpha_i)_{i=-1,n+1}$ and $(\beta_i)_{i=-1,n+1}$ are the parameters (similar to $V_i$ in Eq. (14)) for radial and axial coordinates respectively, and $\xi$ is defined in the undeformed state as:

$$\xi = \frac{s}{l_0} \quad \text{for } 0 \le s \le l_0 \qquad (18)$$

In order to determine parameters in a unique manner, $n + 3$ linearly independent equations are required. Interpolation at each node provides $m + 1$ equations. Moreover, mechanical aspects should be taken into account: boundary conditions relative to angles or moments at membrane extremities give two additional equations. Then the number of knots must be set equal to the number of nodes: $m = n$. In this work, the coordinate of the $i^{\text{iest}}$ knot, $\xi_i$, is set equal to $s(N^i)/l_0$. In the case of cylindrical membrane, boundary conditions used to calculate parameters are of the natural end type [19]:

$$\frac{d^2 r_0}{d\xi^2}(0) = \frac{d^2 r_0}{d\xi^2}(1) = 0 \quad \text{and} \quad \frac{d^2 z_0}{d\xi^2}(0) = \frac{d^2 z_0}{d\xi^2}(1) = 0 \qquad (19)$$

Figure 2 presents an example of a B-spline interpolation for ten nodes. The spline boundary conditions are similar to those in Eq. (19). Note that the first and the second derivatives of the spline are continuous.

<figure><img src="image_4.png" /><figcaption>Figure 2: A B-spline interpolation $S(\xi)$, its first derivative $S'(\xi)$ and its second derivative $S''(\xi)$ with natural end boundary conditions</figcaption></figure>
Equations (19) can be directly injected in the splines definition (16), (17) and interpolations become:

$$r_0(\xi) = \sum_{i=0}^{n} \bar{B}_i(\xi) \alpha_i \qquad (20)$$

$$z_0(\xi) = \sum_{i=0}^{n} \bar{B}_i(\xi) \beta_i \qquad (21)$$

in which $(\bar{B}_i(\xi))_{i=0,n}$ are linear combinations of the previous $(B_i(\xi))_{i=-1,n+1}$.

As in the FEM, our spline finite element is considered isoparametrical:

$$u_r(\xi) = \sum_{i=0}^{n} \bar{B}_i(\xi) \gamma_i \qquad (22)$$

$$u_z(\xi) = \sum_{i=0}^{n} \bar{B}_i(\xi) \delta_i \qquad (23)$$

where $(\gamma_i)_{i=0,n}$ and $(\delta_i)_{i=0,n}$ are spline parameters associated with radial and axial displacements respectively. Therefore, the deformed membrane material points coordinates are interpolated by:

$$r(\xi) = r_0(\xi) + u_r(\xi) = \sum_{i=0}^{n} \bar{B}_i(\xi) (\alpha_i + \gamma_i) \qquad (24)$$

$$z(\xi) = z_0(\xi) + u_z(\xi) = \sum_{i=0}^{n} \bar{B}_i(\xi) (\beta_i + \delta_i) \qquad (25)$$

(26)

Finally, all parameters can be calculated using the following relations:

$$r_0^i = \sum_{j=0}^{n} A_{ij} \alpha_j \qquad (27)$$

$$z_0^i = \sum_{j=0}^{n} A_{ij} \beta_j \qquad (28)$$

$$u_r^i = \sum_{j=0}^{n} A_{ij} \gamma_j \qquad (29)$$

$$u_z^i = \sum_{j=0}^{n} A_{ij} \delta_j \qquad (30)$$

where $u_r^i$ and $u_z^i$ are radial and axial displacement of node i, and $A_{ij} = \bar{B}_j(\xi_i)$ for $i = 0, n$ and $j = 0, n$.
### 3.3 Properties

Noting A the matrix containing the terms $A_{ij}$ and $A^{-1}$ its inverse, useful derivative formulae can be established:

$$ \frac{\partial r}{\partial u_r^i}(\xi) = \sum_{j=0}^{n} \bar{B}_j(\xi) A_{ji}^{-1} = D_i(\xi) \qquad \frac{\partial r}{\partial u_z^i}(\xi) = 0 \qquad (31) $$

$$ \frac{\partial z}{\partial u_r^i}(\xi) = 0 \qquad \frac{\partial z}{\partial u_z^i}(\xi) = \sum_{j=0}^{n} \bar{B}_j(\xi) A_{ji}^{-1} = D_i(\xi) \qquad (32) $$

$$ \frac{\partial \dot{r}}{\partial u_r^i}(\xi) = \sum_{j=0}^{n} \dot{B}_j(\xi) A_{ji}^{-1} = \dot{D}_i(\xi) \qquad \frac{\partial \dot{r}}{\partial u_z^i}(\xi) = 0 \qquad (33) $$

$$ \frac{\partial \dot{z}}{\partial u_r^i}(\xi) = 0 \qquad \frac{\partial \dot{z}}{\partial u_z^i}(\xi) = \sum_{j=0}^{n} \dot{B}_j(\xi) A_{ji}^{-1} = \dot{D}_i(\xi) \qquad (34) $$

in which the dot stands for differentiation with respect to $\xi$, and $(D_i(\xi))_{i=0,n}$ and $(\dot{D}_i(\xi))_{i=0,n}$ are functions which have to be computed only once.

## 4 B-SPLINE ELEMENT FORMULATION

Using the previous interpolation, we define the out of balance force $\mathbf{R}(\mathbf{U}, p)$ by:

$$ R(\mathbf{u}, \delta\mathbf{u}, p) = \delta\mathbf{U}^T \mathbf{R}(\mathbf{U}, p) \qquad (35) $$

with:

$$ \mathbf{R}(\mathbf{U}, p) = \mathbf{F}^{\text{int}}(\mathbf{U}) - \mathbf{F}^{\text{ext}}(\mathbf{U}, p) \qquad (36) $$

in which $\mathbf{F}^{\text{int}}$ and $\mathbf{F}^{\text{ext}}$ are respectively vectors of internal and external forces, and $\mathbf{U}$ is the nodal displacements vector which contains $(u_r^i, u_z^i)_{i=0,n}$. The equilibrium equation (5) becomes:

$$ \mathbf{R}(\mathbf{U}, p) = \mathbf{0} \qquad (37) $$

The internal forces vector is given by:

$$ \mathbf{F}^{\text{int}}(\mathbf{U}) = \int_{0}^{l_0} 2\pi r_0 h_0 \left[ S_1 \left\{ \frac{\partial E_1}{\partial \mathbf{U}} \right\} + S_2 \left\{ \frac{\partial E_2}{\partial \mathbf{U}} \right\} \right] ds \qquad (38) $$

and the deformation dependent external forces vector can be written as:

$$ \mathbf{F}^{\text{ext}}(\mathbf{U}, p) = \int_{0}^{l} 2\pi p \frac{r}{\sqrt{r'^2 + z'^2}} \left[ z' \left\{ \frac{\partial u_r}{\partial \mathbf{U}} \right\} - r' \left\{ \frac{\partial u_z}{\partial \mathbf{U}} \right\} \right] ds \qquad (39) $$
As mentioned earlier, all terms in Eqs (38) and (39) are completely determined by the B-spline interpolation.

From the differentiation of Eq. (36) with respect to the nodal displacement vector, the tangent stiffness matrix is defined by:

$$ \mathbf{K}_t = \mathbf{K}_t^{\text{int}} - \mathbf{K}_t^{\text{ext}} = \frac{\partial \mathbf{F}^{\text{int}}}{\partial \mathbf{U}} - \frac{\partial \mathbf{F}^{\text{ext}}}{\partial \mathbf{U}} \quad (40) $$

where $\mathbf{K}_t^{\text{int}}$ and $\mathbf{K}_t^{\text{ext}}$ are the tangent stiffness matrices given by differentiation of the internal and external vectors respectively.

All the previous expressions are analytically evaluated by using the differentiation chain rule and the relations of paragraph 3.3.

## 5 NUMERICAL SOLUTION

Spatial integration of vectors and matrices is performed by a simple Simpson's method. Faster methods will be studied in further work to improve our program.

The previous problem (37), (38) and (39) is highly non-linear, because of large strains and non-linearity of the constitutive equation. External forces vector depends on the current configuration, therefore it has to be calculated in each step. Moreover, such problem exhibits numerical difficulties because of the presence of bifurcation or limit points [6]. Some authors used manual control of the loading variable [20] but it is not as efficient as an automatic continuation method. In this study, a combination of the classical Newton-Raphson method and the arc-length method [21] is used to solve the problem. A similar method was successfully implemented in the case of the classical finite element method [9].

Denote $\mathbf{U}_e$ and $p_e$ an equilibrium point obtained on a previous loading step. In order to find a new equilibrium configuration, the problem consists in finding both displacement and pressure increments $\Delta \mathbf{U}$ and $\Delta p$ that satisfy simultaneously the following equations:

$$ \mathbf{R}(\mathbf{U}_e + \Delta \mathbf{U}, p_e + \Delta p) = 0 \quad (41) $$

$$ A(\Delta \mathbf{U}, \Delta p) = 0 \quad (42) $$

where $A$ is the arc-length constraint defined by:

$$ A(\Delta \mathbf{U}, \Delta p) = (\|\Delta \mathbf{U}\|^2 + \psi^2 \|\Delta \mathbf{F}^{\text{ext}}\|^2) - da^2 \quad (43) $$

in which $\Delta \mathbf{F}^{\text{ext}}$ is the loading increment, $\psi$ is a scale factor between displacement and force, and $da$ an user-defined arc-length. Here $\psi$ is set to zero and the arc-length method reduces to the classical displacement control method [22]. Using the Newton-Raphson algorithm, equilibrium points are obtained by iterations using the tangent stiffness matrix defined earlier (40). Consider the algorithm at the $k^{\text{iest}}$ iteration. Next values of $\mathbf{R}$, $A$,
$\Delta U$ and $\Delta p$ are calculated by:

$$ \mathbf{R}_k = \mathbf{R}_{k-1} + \mathbf{K}_t \, d\mathbf{U} - dp \, \mathbf{f}^{\text{ext}} = 0 \quad (44) $$

$$ A_k = A_{k-1} + 2 \Delta \mathbf{U}_{k-1}^T \, d\mathbf{U} = 0 \quad (45) $$

$$ \Delta \mathbf{U}_k = \Delta \mathbf{U}_{k-1} + d\mathbf{U} \quad (46) $$

$$ \Delta p_k = \Delta p_{k-1} + dp \quad (47) $$

where $d\mathbf{U}$ and $dp$ are changes of displacement and pressure increments to apply in the next iteration, and $\mathbf{f}^{\text{ext}}$ is the reduced external forces vector defined by:

$$ \mathbf{f}^{\text{ext}}(\mathbf{U}) = \frac{1}{p} \mathbf{F}^{\text{ext}}(\mathbf{U}, p) \quad (48) $$

Due to the follower force the tangent stiffness matrix is non-symmetric. Thus the arc-length constraint equation (45) can be considered as a bordered equation [23] of Eq. (44) without changing the nature of the resolution algorithm. Therefore, $d\mathbf{U}$ and $dp$ are the solutions of the extended following system:

$$ \begin{bmatrix} \mathbf{K}_t & -\mathbf{f}^{\text{ext}} \\ 2\Delta\mathbf{U}_{k-1}^T & 0 \end{bmatrix} \begin{Bmatrix} d\mathbf{U} \\ dp \end{Bmatrix} = - \begin{Bmatrix} \mathbf{R}_{k-1} \\ A_{k-1} \end{Bmatrix} \quad (49) $$

The present scheme is initiated by the forward-Euler tangential predictor solution [24].

# 6 NUMERICAL EXAMPLES

In order to illustrate our method, some numerical examples are examined. In these examples, only uniform radius cylinders are considered. Non-uniform cases may have an influence on non-linear deformations and deserve further studies. In order to simplify the discussion the results are presented in terms of dimensionless quantities. The aspect ratio and the inflation parameter are respectively noted $A_S$ and $I_P$, and defined by:

$$ A_S = \frac{l_0}{2r_0} \qquad I_P = \frac{r_0 \Delta p}{2Ch_0} \quad (50) $$

The dimensionless maximum radius at the middle of the cylinder is given by:

$$ R_{\text{max}} = \frac{r(s = l/2)}{r_0(s = 0)} \quad (51) $$

## 6.1 Influence of the aspect ratio

In order to compare our results with those of Khayat et al. [25], this study is restricted to the neo-Hookean material i.e.: $\alpha = 0$ in Eq. (6). Figure 3 presents the inflation parameter as a function of the maximum radius for three values of the aspect ratio. Our results are similar to those previously obtained [25]. As in the classical spherical case, neo-Hookean cylinders exhibit a maximum pressure for all geometries. At this maximum value, $I_P$ starts decreasing as $R_{\text{max}}$ continues to increase. As it could be expected, it can be observed that long cylinders are easier to inflate than short ones.
<figure><img src="image_5.png" /><figcaption>Figure 3: Inflation parameter versus maximum radius for different aspect ratios (neo-Hookean material)</figcaption></figure>

## 6.2 Influence of the Mooney parameter

In this study the aspect ratio is set to 5 and the inflation of uniform radius cylinders is examined for different values of the Mooney parameter $\alpha$. Figure 4 shows the inflation parameter versus the maximum radius in these cases. Behaviours of cylindrical mem-

<figure><img src="image_6.png" /><figcaption>Figure 4: Inflation parameter versus maximum radius for different Mooney parameters ($A_s = 5$)</figcaption></figure>
branes are similar to those corresponding to the inflation of spherical membranes [26]. As mentioned earlier, for a neo-Hookean membrane the curve presents a maximum pressure. As $\alpha$ increases, this maximum value tends to disappear: curves corresponding to $\alpha = 0.1$ and $\alpha = 0.2$ have both a local maximum and a local minimum pressures, and curves with $\alpha = 0.35$ and $\alpha = 0.5$ grow monotonically. For values of $\alpha$ different from 0, our results are similar to those of Khayat et al.. However for $\alpha = 0$ the results slightly differ: our curve decreases faster than the curve they obtained. As the neo-Hookean inflation of long cylinders without pre-stretching is recognized as a critical case [6], this difference might be due to boundary conditions. In the paper cited above, the evolution of the mid-cylinder radius is imposed as one of the boundary condition of the system of differential equations. In the present work, the value of $R_{max}$ is not imposed as a boundary condition: with the arc-length constraint, both the maximum radius and the pressure are calculated to ensure equilibrium.

## 6.3 Inflation profiles

Finally, two inflation profile evolutions are shown in Figure 5 for $\alpha = 0$ and $\alpha = 0.2$. In both cases, the aspect ratio is equal to 5. In the neo-Hookean case ($\alpha = 0$), the profiles are

<figure><img src="image_7.png" /><figcaption>Figure 5: Inflation profiles of cylinders with $A_S = 5$ for different values of $I_P$</figcaption></figure>
regular until $I_P$ reaches the maximum pressure 0.84. In the second part of the inflation, $I_P$ monotonically decreases and simultaneously a bulge develops in the mid-cylinder. This bulge never disappears. Similar observations are made by Khayat and Derdouri [27] for the same cylinder, and by Verron et al. [28] for an ellispoidal plane membrane. In both references, authors use the classical finite element method. In the case of $\alpha = 0.2$, the membrane inflates regularly and profiles are dome-like shape for all values of the inflation parameter.

# 7 CONCLUDING REMARKS

In this paper, the B-splines interpolation method is successfully used to model the axisymmetrical inflation of cylindrical rubberlike membranes. The numerical implementation does not present major difficulties and it highly reduces the number of degrees of freedom needed to obtain a good description of the membrane, thanks to the continuity of both B-spline functions and their derivatives. Moreover the use of the Newton-Raphson algorithm coupled with an arc-length method considered as a bordered equation ensures convergence of the scheme even in the vicinity of limit points.

As a perspective, further works are in progress on the B-splines boundary conditions to model the case of circular membranes with ends on the symmetry axis. The main goal of this approach remains the developments of small size models which will be associated with optimization procedures to determine rubberlike material parameters using experimental bubble inflation results [18].

# ACKNOWLEDGEMENTS

The authors wish to thank Dr. E. Arnoult who introduced them to the efficiency of the B-splines interpolation.
## REFERENCES

[1] H. Alexander, The tensile instability of an inflated cylindrical membrane as affected by an axial load, *Int. J. Mech. Sci.*, **13**, 87-95 (1971).

[2] W. Johnson and P. D. Soden, The discharge characteristics of confined rubber cylinders, *Int. J. Mech. Sci.*, **8**, 213-225, (1966).

[3] A. D. Kydonieſs and J. M. Spencer, Finite axisymmetrie deformations of an initially cylindrical membrane, *Q. J. Mech. Appl. Math.*, **22**, 87-95, (1969).

[4] R. Benedict, A. Wineman, and W. H. Yang, The determination of limiting pressure in simultaneous elongation and inflation of nonlinear elastic tubes, *Int. J. Solids Structures*, **15**, 241-249, (1979).

[5] G. A. Duffet and B. D. Reddy, The solution of multi-parameter systems of equations with application to problems in nonlinear elasticity, *Comp. Meth. Appl. Mech. Engng*, **59**, 179-213, (1986).

[6] R. E. Khayat, A. Derdouri, and A. Garcia-Réjon, Inflation of an elastic cylindrical membrane: non-linear deformation and instability, *Int. J. Solids Structures*, **29**, 69-87, (1992).

[7] S. Kyriakides, Propagation instabilities in structures, *Adv. Appl. Mech.*, **30**, 67-189, (1994).

[8] S. Reese and P. Wriggers, A finite element method for stability problems in finite elasticity, *Int. J. Numer. Meth. Engng*, **38**, 1171-1200, (1995).

[9] J. Shi and G. F. Moita, The post-critical analysis of axisymmetrie hyper-elastic membranes by the finite element method, *Comput. Methods Appl. Mech. Engng*, **135**, 265-281, (1996).

[10] C. De Boor, *A practical guide to splines*, New-York: Springer, (1978).

[11] C. T. Shik, On spline finite element, *J. Comput. Math.*, **1**, 50-72, (1979).

[12] Y. K. Cheung and S. C. Fan, Static analysis of right box girder bridges by finite strip method, *Proc. Inst. Civ. Engineers. Part 2*, **75**, 311-323, (1983).

[13] A. Gupta, J. Kiusalaas, and M. Saraph, Cubic b-spline for finite element analysis of axisymmetrie shells, *Comput. Struct.*, **38**, 463-468, (1991).

[14] A. H. Vermeulen and G. R. Heppler, Structural analysis of shells by the b-spline field approximation method, *Comput. Struct.*, **68**, 167-179, (1998).
[15] B. P. Patel, M. Ganapathi, and J. Saravanan, Shear flexible field-consistent curved spline beam element for vibration analysis, *Int. J. Numer. Meth. Engng*, **46**, 387-407, (1999).

[16] S. W. Yuen and G. M. Van Erp, Transient analysis of thin-walled structures using macro spline finite elements, *Engng Struct.*, **21**, 255-266, (1999).

[17] L. R. G. Treloar, The mechanics of rubber elasticity, *Proc. R. Soc. Lond., A351*, 301-330, (1976).

[18] E. Verron, *Experimental and numerical contribution to the blow-moulding and thermoforming processes (in french)*, PhD thesis, Ecole Centrale de Nantes, Nantes, France, (1997).

[19] G. Hammerlin and K. H. Hoffmann, *Numerical mathematics. (Chapter 6)*, Undergraduate texts in mathematics, New-York: Springer Verlag, (1991).

[20] S. Kyriakides and Y. Chang, The initiation and propagation of a localized instability in an inflated elastic tube, *Int. J. Solids Struct.*, **27**, 1085-1111, (1991).

[21] M. A. Crisfield, *Non-linear finite element analysis of solids and structures. Volume 1: Essentials*. (Chapter 9), John Wiley and Sons, (1994).

[22] J. L. Batoz and G. Dhatt, Incremental displacement algorithms for nonlinear problems, *Int. J. Num. Meth. Engng*, **14**, 1262-1267, (1979).

[23] E. Riks and C. C. Rankin, An incremental approach to the solution of snapping and buckling problems, *Int. J. Solids Struct.*, **15**, 529-551, (1979).

[24] R. Kouhia and M. Mikkola, Some aspects on efficient path-following, *Comput. Struct.*, **72**, 509-524, (1999).

[25] R. E. Khayat, A. Derdouri, and A. Garcia-Réjon, Multiple contact and axisymmetric inflation of a hyperelastic cylindrical membrane, *J. Mech. Eng. Sci.*, **207**, 175-183, (1993).

[26] E. Verron, R. E. Khayat, A. Derdouri, and B. Peseux, Dynamic inflation of hyperelastic spherical membranes, *J. Rheol.*, **43**, 1083-1097, (1999).

[27] R. E. Khayat and A. Derdouri, Inflation of hyperelastic cylindrical membranes as applied to blow moulding. Part I. Axisymmetric case, *Int. J. Numer. Meth. Engng*, **37**, 3773-3791, (1994).

[28] E. Verron, G. Marckmann, and B. Peseux, Dynamic inflation of non-linear elastic and viscoelastic rubberlike membranes, *Int. J. Numer. Meth. Engng*, submitted.