Numerical solution of hyperelastic membranes by energy minimization

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 membranes-and over a more complex structure-a parabolic antenna-using the Saint-Venant Kirchhoff and neo-Hookean models. Whenever available, analytical or semi-analytic solutions are used to validate the finite element results.


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 nonlineari ties. 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 / or the displacementŨ , carrying positions X in the refer ence configuration X 0 to positions x /(X) X +Ũ in the current configuration X.
The reference boundary surface oX 0 of X 0 is decom posed into disjoint parts oX 0U and oX 0T , over which the boundary conditions are written as whereŨ g denotes given (prescribed) displacement,T g the given nominal traction vector, P the first Piola Kirchhoff stress tensor (not symmetric) related to the Cauchy stress tensor r by PF ÀT Jr, where F ¼ I þ r XŨ is the deformation gradient, J detF. Use will also be made of the second Piola Kirchhoff stress tensor R (symmetric) related to the first one by P FR. NansonÕs formula [5,8] for the transformation of sur face elements shows that the nominal traction vectorT is in general a deformation sensitive loading depending on displacementŨ and F. In the case of an inflated mem brane, one part of oX 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 vectorT there writes 2.1.1. Conservative surface loading [7] A prescribed nominal surface loadingT g ðX ;Ũ ; FÞ de fined over the reference surface S 0 derives from external potential energy V ext ðŨ Þ if 8ṽ Ã virtual velocity field satis fyingṽ Ã ¼0 over S 0U , the virtual power of the surface loading can be expressed as: where oV ext oŨ is the Gateaux derivative of V ext , which is related to the first variation of V ext by: 8dŨ ; dV ext ¼ oV ext oŨ ½dŨ .
When (4) is satisfied,T g is referred to as a conserva tive surface loading.

Definition of the pressurized volume h
From (1), the boundary line oS 0 of S 0 is subjected to the prescribed displacements Let us denote by S lat the surface swept in the three dimensional space, up to current time, by line oS. The pressurized volume h is then defined as that bounded by surfaces S 0 , S and S lat at current time (Fig. 1a).
More often than not the prescribed displacementsŨ g are zero, so that S lat is the empty set and h is merely the volume between S 0 and S (Fig. 1b).

Properties of volume h
• The pressurized volume h is a functional of the dis placement fieldŨ . In order to express this, we write h ¼ hðŨ Þ. • In general, the explicit expression for hðŨ Þ as a func tional ofŨ is unknown. Nevertheless, its first varia tion dh corresponding to a variation dŨ ofŨ can be determined by the following relation:

Pressure potential energy
Assuming that the whole line oS 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]: where the symbol ò stands for a primitive. In particular, if the pressure is so controlled that it is independent of volume h, then the pressure potential energy is V ext ðŨ Þ ¼ phðŨ Þ. If the pressure varies as a function of volume h according to the law ph K constant, then the potential energy is V ext ðŨ Þ ¼ K log hðŨ Þ.
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 h we transform it into surface integrals by means of the Gauss theorem: where x ! is the current vector position of a particle on the surface S, n ! the outward normal to S. The surface integral (8) becomes a discrete summation over the finite elements: 2.2. The internal potential energy (strain energy) We consider a hyperelastic material described by a strain energy w per unit reference volume (rather than per unit mass), function of either the right stretch tensor or the Green strain tensor E ¼ 1 2 ðC IÞ. Coefficients g ij are the covariant components of metric tensor g defined in the current configuration and vectors G i ! form the dual of the natural basis in the reference configuration.
The first and second Piola Kirchhoff stress tensors P and R, respectively, are given by wherew(F) w(C F T F). By introducing the strain energy the virtual stress power along a virtual velocity fieldṽ Ã can be expressed as a function of the Gateaux derivative From relation (12), one says that the Piola Kirchhoff stress P derives from internal potential energy V int ðŨ Þ.
For an isotropic material without internal con straints, the strain energy per unit volume w can be expressed in terms of the invariants of the right Cauchy Green tensor C w ¼ wðI 1 ðCÞ; I 2 ðCÞ; I 3 ðCÞÞ ð13Þ with I 1 (C) tr C, Two constitutive models are used for numerical purposes: (i) the Saint Venant Kirchhoff model, for which the strain energy is where k and l are material parameters; (ii) and the compressible isotropic model described by a strain energy of neo Hookean type which yields from the constitutive law (10) If the body undergoes a rigid body motion, the defor mation gradient F is equal to an orthogonal tensor, ten sor 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 m related to parameters (k, l) as in small strains: When using the triangular finite elements (Fig. 2), the mixed components of C are defined by where D and d are the metric tensors in the reference and current configuration with their covariant components defined as:

½D ¼Ã
Since the strain tensor is constant inside each finite element, the strain energy is simply obtained as the prod uct of the strain energy by the volume of the element where S 0 is the reference area of the membrane element and h the reference thickness.

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 poten tial energy V ext (e.g. relation (9)), (ii) and the constitutive law is hyperelastic, so that the stresses derive from the internal potential energy V int (11).
Then, the total potential energy VV ext + V int is stationary at the solution displacement fieldŨ .
Proof. The principle of virtual power gives : 8ṽ Ã CAH, Relation (20) shows that the first variation of V is dV ¼ oV oŨ ½dŨ ¼ 0. h Furthermore it is shown that a stable solutionŨ cor responds 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 suf ficient accuracy and to correctly handle some phenom ena like membrane folding.

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 mem brane 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 solu tion 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 consid ered the pressure as a dead load. These approximations led to a simplified solution which is chosen here for com parison in moderate rotations.

Fichter's semi analytical solution
Let us summarize FichterÕs semi analytical solution given in [2]. For brevityÕs sake, the dimensionless radial co ordinate r/a will be denoted q. The deflection w is searched in the form of a power series: By replacing this expression in the equilibrium equa tion, one obtains the circumferential stress R hh 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] It should be noted that all the coefficients depend on b 0 . By using the boundary condition w(a) 0, Fichter showed that b 0 is related to the PoissonÕs ratio by the fol lowing equation: It follows that all the other coefficients also depend on the PoissonÕs ratio only.

Numerical results
The geometry and mechanical properties of the circu lar membrane are: radius a 0.1375 m, equivalent Young modulus E * E AE h 600 kPa m, Poisson ratio m 0.3. The membrane inflation is modelled by triangu lar finite elements as described in Section 2. The mesh contains 1024 triangular elements.
The initial and deformed shapes of the membrane are shown in Fig. 3a and b. Fig. 4 shows the central deflection (at r 0) versus the pressure. Fig. 4b shows that for small deflections the numerical results are in very good agreement with the semi analytical solution given by Fichter (Relation (21)) which is not valid for higher pressures, as explained above.
Moreover, the proposed numerical scheme gives ex actly the same values as the standard nonlinear finite ele ment method, over the whole load range, for both the Saint Venant Kirchhoff and the neo Hookean models, Eqs. (14) and (15). It should be noted that the neo Hookean model leads to a limit point, at the pressure p 3,006,180 Pa.

Inflation of an orthotropic rectangular membrane
In this section, we deal with the problem of a rectan gular membrane inflation. We first give a semi analytical solution in the case of orthotropic behavior and then its modelling by the triangular finite element based on the minimization of the energy. For simplicity, a comparison between these two solutions is presented for isotropic behavior case only.

Semi analytical solution
Let us consider a rectangular membrane with sides 2a and 2b, clamped at its rim, subjected to a uniform pres sure 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 n X/a and g Y/b.
We build the semi analytical solution in the large deflection case: the in plane displacements of the mem brane u(n,g), v(n,g) and the deflection w(n,g). In order to automatically satisfy the kinematic boundary condi tions, we assume the following polynomial expressions for the displacement components: g(n,g) is chosen so that the boundary kinematic condi tions are satisfied.
gðn; gÞ ¼ ð1 n 2 Þð1 g 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 rec tangular membrane. For this we need to evaluate the fi nite strain components: 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. Conse quently, 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: The D ij terms depend on the material constants E X , E Y and m XY 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 m XY the Poisson ratio (m XY /E X m YX /E Y ). Then, the strain en ergy density for an orthotropic material takes the form: The strain energy is obtained after integration of the strain energy density over the membrane: A symbolic calculus of V 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 approx imation 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 int by means of the actual symbolic calculation systems.

Pressure potential energy
The pressure potential energy is the product of the pressure by the integral of the deflection on the mem brane. 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 approxi mation 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: 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 displace ment is computed by minimizing the total potential en ergy V tot V int + V ext relatively to the coefficients of the power series: The analytic expression of V tot upto n 5 (Eq. (26)) is obtained by means of symbolic calculus tools. The min imization is achieved numerically by a Newton like iter ative 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.

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 coeffi cients of the displacement series by a numerical minimi zation of the total potential energy. We used some data issued from simple tension experiments for measurement of YoungÕs modulus listed in Table 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.

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 con vergence 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 cri terion of convergence relates to total energy. The varia tion of the energy between successive iteration must be less than 10 À6 N m.

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 trian gular finite element noted w fe . In order to underline the difference between the two solutions, lets define the rel ative difference between them as: (w fe w a )/max(w fe ). We have represented on Fig. 7 the distribution of the rel ative 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 mem brane is about 2.5%. The restriction on the order of the polynomial series can explain this difference.

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 parabo loids. One of them is transparent to the electromagnetic waves and the second one acts as a concave reflector. Un der 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 * E AE h 350,000 Pa m and Poisson ration m 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 incompatibil ity, the torus warps out of its plane and displays more vis ible folds regularly distributed along the circumference.    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.

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.