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

HAL
open science

# Model reduction method: an application to the buckling analysis of laminated rubber bearings

Stéphane Lejeunes, Adnane Boukamel, Bruno Cochelin

► To cite this version:

Stéphane Lejeunes, Adnane Boukamel, Bruno Cochelin. Model reduction method: an application to the buckling analysis of laminated rubber bearings. Revue Européenne des Éléments Finis, 2006, 15 (1-3), pp.281-292. 10.3166/remn.15.281-292. hal-00089047

HAL Id: hal-00089047

https://hal.science/hal-00089047v1

Submitted on 9 Aug 2006

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 CC BY 4.0 - Attribution - International License
# Model reduction method:
an application to the buckling analysis
of laminated rubber bearings

Stéphane Lejeunes — Adnane Boukamel — Bruno Cochelin

Ecole Généraliste d'Ingénieurs de Marseille
IMT/Technopôle Chateau-Gombert, F-13451 Marseille Cedex 20
Laboratoire de Mécanique et d'Acoustique de Marseille
31 chemin Joseph-Aiguier, F-13402 Marseille cedex 20
lejeunes@esim.fr, adnane.boukamel@egim-mrs.fr, bruno.cochelin@egim-mrs.fr

**ABSTRACT.** In this paper, we apply a model reduction method to find the equilibrium state at finite strain of geometrically complex structures which have periodic properties in one direction and exhibit a non-linear material behavior. This method, based on a finite-element approach, consists in projecting the unknowns fields onto a polynomial basis in order to reduce the size of the problem. This method was combined with a continuation resolution scheme to find the instabilities of a laminated rubber bearing subjected to compression loading. Comparisons with standard finite-element models show the reliability of the present method.

**RÉSUMÉ.** On présente, dans ce papier, une application d'une technique de réduction de modèles pour résoudre les problèmes d'équilibre de structures géométriquement complexes possédant des propriétés de périodicité selon une direction et dont le comportement présente des non-linéarités géométriques et matérielles. Basée sur une formulation en éléments finis, cette méthode consiste à projeter sur une base polynômiale, bien choisie vis-à-vis des propriétés d'invariance, les champs inconnus. Elle permet ainsi de réduire la dimension du problème à résoudre. Le couplage d'une telle méthode avec une méthode de continuation, nous permet d'analyser les instabilités d'un lamifié élastomère-métal en compression. Des comparaisons avec des modèles éléments-finis standard sont présentées.

**KEYWORDS:** model reduction, hyperelasticity, finite element, stability, bifurcation, continuation.

**MOTS-CLÉS** : réduction de modèle, hyperélasticité, éléments-finis, stabilité, bifurcation, continuation.
# 1. Introduction

Rubber bearings are used in many fields of application such as seismic base isolation and aeronautical construction. In most cases, these parts are designed to damp or transmit complex loads in which compression is the main component. A fine numerical method is often required for simulating this kind of parts, especially for performing buckling and post-buckling analysis. However, classical numerical tools (such as the standard finite-element method) result in some cases in prohibitively large models. Various approaches which decrease the model size have been proposed.

The first of these models is based on kinematically enriched beam theory or on behavioral assumptions (macro-models), where the bearing is taken to be an elastic (or viscoelastic) column. These models are useful for predicting the behavior in single loading cases or when it is required to determine the buckling load (Kelly 2003, Iizuka 2000). Sub-structuring and homogenization techniques have been proposed to avoid the problems associated with the prohibitively large size of finite element models. We do not intend here to review the whole literature on this subject. Readers interested can refer to (Dumontet 1990, Devries 1998, Léné *et al.* 2001). The last group of models are based on specific elements. Boukamel, for instance has used a decomposition method into Fourier series to deal with the case of a cylindrical rubber bearing (Boukamel 1988) subjected to overall loading.

The model reduction method used in this study is based on the semi-analytical finite element method (which is usually called the finite-strip method) developed by Cheung in the mid 70's for the analysis of shell elastic structures. The finite-strip method is now applied in many fields (Cheung *et al.* 1995, Zhong *et al.* 1998). Our aim here is to extend this concept to the case of non linear behavior, where the non linearity is due to both the large strains and the quasi-incompressibility constraint.

The present paper deals with the buckling analysis of laminated structures. This mechanical analysis is carried out with an arc-length continuation method (Crisfield 1997, Riks 2004). The singular points are detected and analyzed to identify the presence of a bifurcation. In the case of a bifurcation, a simple branch switching strategy is adopted to follow the post-buckling path.

# 2. Modeling isotropic nearly-incompressible hyperelastic materials

Methods of modeling isotropic nearly-incompressible hyperelastic bodies are by now readily available. Several books contain a detailed review of the basic and more advanced concepts developed in the framework of hyperelasticity (Fu *et al.* 2001, Holzapfel 2004). The aim of this section is mainly to describe the theoretical framework required to write the variational form of the equilibrium equations.
## 2.1. Constitutive equations

Taking **F** to denote the deformation gradient and $J = \det \mathbf{F}$ to denote the volume variation, the deformation gradient can be split into a volumetric part $J^{1/3}\mathbf{I}$ and an isochoric part $\overline{\mathbf{F}}$ such that: $\mathbf{F} = (J^{1/3})\overline{\mathbf{F}}$. A similar decomposition can be obtained with the right Cauchy-Green tensor: $\mathbf{C} = \mathbf{F}^T \mathbf{F} = (J^{2/3})\overline{\mathbf{C}}$.

In line with Miehe (1994) (and other authors), we consider an additive split of the isotropic Helmholtz free energy into a volumetric part and an isochoric part:

$$ \psi(\mathbf{F}) = \psi_{iso}(\overline{\mathbf{C}}) + \psi_{vol}(J) \quad [1] $$

Taking the first Piola-Kirchoff stress tensor $\mathbf{\Pi}$, the following constitutive equation is obtained:

$$ \mathbf{\Pi} = \overline{\mathbf{\Pi}} : \mathbb{P} + p \mathbf{J} \mathbf{F}^{-T} \quad [2] $$

where $\overline{\mathbf{\Pi}} = 2\overline{\mathbf{F}}\partial\psi_{iso}/\partial\overline{\mathbf{C}}$. The fourth-order tensor $\mathbb{P}$ is a deviatoric operator for the mixed configuration:

$$ \mathbb{P} = \frac{\partial \overline{\mathbf{F}}}{\partial \mathbf{F}} = J^{-1/3} \left( \mathbb{I} - \frac{1}{3} \mathbf{F} \otimes \mathbf{F}^{-T} \right) \quad [3] $$

where $\mathbb{I}$ is the fourth order identity tensor and $\otimes$ denotes the tensor product (with e.g. $\mathbf{A} = \mathbf{B} \otimes \mathbf{C}$ corresponding to $A_{ijkl} = B_{ij}C_{kl}$ in indicial notation). The hydrostatic pressure $p$ is defined from a compressibility law:

$$ p = \frac{\partial \tilde{\psi}^{vol}(J)}{\partial J} = \tilde{\psi}_{J}^{vol} \quad [4] $$

One can find many forms of the strain energy density based on statistical physics or phenomenological assumptions; see (Holzapfel 2004) for a review. In this paper, we use the so-called generalized Mooney-Rivlin model, extended from (Mooney 1940), and a compressibility law proposed in (Hartmann et al. 2003):

$$ \psi_{iso}(\overline{I}_1, \overline{I}_2) = a_{10}(\overline{I}_1 - 3) + a_{01}(\overline{I}_2 - 3) \quad [5] $$

$$ \psi_{vol}(J) = \frac{k}{2} G(J)^2 = \frac{k}{2} \left( \frac{1}{10} (J^4 - J^{-6}) \right)^2 \quad [6] $$

where $\overline{I}_1$ and $\overline{I}_2$ are the first principal invariants of $\overline{\mathbf{C}}$. The material characteristics are $a_{10}$ and $a_{01}$, which are the two Mooney-Rivlin parameters, and $k$ is the modulus of compressibility.

## 2.2. Variational form

Several variational principles exist in the case of rubber-like behavior, and particular attention has been paid to the multi-field methods based on Lagrange multipliers and Hellinger-Reisner variational theorem since the 70's, see for example the work
of Hermann, Hughes and Malkus, Simo and Taylor, Glowinski and Le Tallec and the references therein (Malkus et al. 1978, Glowinski et al. 1982, Herrmann 1963). In the present study, we use an approach adapted to modeling nearly-incompressible behavior, consisting of introducing the following two field functional of the total potential energy (Rüter et al. 2000):

$$P(\mathbf{u}, J) = \int_{\Omega_0} \left[ \tilde{\psi}_{iso}(\overline{\mathbf{C}}(\mathbf{u})) + \frac{1}{2}kG(J)^2 \right] dV - P_{ext}(\mathbf{u}) \quad [7]$$

where the external load potential is defined by

$$P_{ext}(\mathbf{u}) = \int_{\Omega_0} \mathbf{f} \cdot \delta\mathbf{u} dV - \int_{\partial\Omega_F} \mathbf{T} \cdot \delta\mathbf{u} dS \quad [8]$$

where $\mathbf{f}$ and $\mathbf{T}$ are the prescribed body or surface forces acting on either the domain $\Omega_0$ or the boundary $\partial\Omega_F$. Partial Legendre transformation of [7] leads to the same expression for the potential energy as with the perturbed Lagrange-multiplier method:

$$P(\mathbf{u}, p) = \int_{\Omega_0} \left[ \tilde{\psi}_{iso}(\overline{\mathbf{C}}(\mathbf{u})) + pG(J) - \frac{1}{2k}p^2 \right] dV - P_{ext}(\mathbf{u}) \quad [9]$$

The aim of this formulation is to solve the saddle-point problem:

$$\min_{\mathbf{u} \in \mathcal{V}} \max_{p \in \mathcal{P}} P(\mathbf{u}, p) \quad [10]$$

where $\mathcal{V} = (H_0^1(\Omega_0))^2$ and $\mathcal{P} = L_2(\Omega_0)$ are the space of smooth kinematic and pressure fields such that $\mathbf{u} = 0$ on the boundary $\partial\Omega_u$. The stationary of the functional [9] gives the equilibrium and quasi-compressibility equations

$$\left\{ \begin{array}{l} \displaystyle \int_{\Omega_0} (\overline{\mathbf{\Pi}} : \mathbb{P} + pJG'(J)\mathbf{F}^{-T}) : \delta\mathbf{F} dV - \int_{\Omega_0} \mathbf{f} \cdot \delta\mathbf{u} dV - \int_{\partial\Omega_F} \mathbf{T} \cdot \delta\mathbf{u} dS = 0 \\ \displaystyle \int_{\Omega_0} (G(J) - \frac{1}{k}p)\delta p dV = 0 \end{array} \right. \quad [11]$$

$$\forall \delta \mathbf{u} \in \mathcal{V}, \forall \delta p \in \mathcal{P}$$

### 3. Implementation of the model reduction method

The main idea underlying this method is to approximate the kinematics and the pressure fields in the transversal direction of a layer in terms of smooth polynomial shape functions. By projecting the unknown fields onto a sufficiently complete polynomial basis, we expect to obtain the boundary effects occurring near the free edges. This approximation leads to the condensation of a geometrical dimension before the discretization process is performed. A plane strain elastomeric layer will therefore be modeled by an assembly of 1 dimensional finite elements (figure 1).
<figure><img src="image_3.png" /><figcaption>Figure 1. 2D-1D reduction</figcaption></figure>

We take $(X, Y)$ to denote the position of a material point in the undeformed configuration. The polynomial expansion of the unknown fields is carried out using the mapping $\xi = 2X/L$:

$$ \mathbf{u}(X, Y) = \sum_{i=0}^{n_u} \mathbf{u}^i(Y) T_i(\xi) \quad [12] $$

$$ p(X, Y) = \sum_{i=0}^{n_p} p^i(Y) T_i(\xi) \quad [13] $$

where $n_u$ and $n_p$ are the polynomial orders of series expansion. The order must be chosen to obtain a compromise between the accuracy, the computing time and the numerical stability of the finite elements.

The polynomial basis consists of order one Lagrangian functions (to account for kinematic constraints such as translations and rotations) and bubble functions. Using the order $n$ Legendre polynomial $L_n(\xi)$, we have:

$$ T_0(\xi) = \frac{1 - \xi}{2}, \quad T_1(\xi) = \frac{1 + \xi}{2}, \quad T_n(\xi) = \frac{L_n(\xi) - L_{n-2}(\xi)}{\sqrt{2(2i-1)}} \quad [14] $$

Starting with equations [12] and [13], we can make a finite-element approximation in the vertical direction with Lagrangian shape functions $N(Y)$. The kinematic and pressure fields in an element are:

$$ \mathbf{u}(X, Y) = \sum_{i=0}^{n_u} \sum_{j=1}^{l_u} N_u^j(Y) \mathbf{u}_j^i T_i(\xi) \quad [15] $$

$$ p(X, Y) = \sum_{i=0}^{n_p} \sum_{j=1}^{l_p} N_p^j(Y) p_j^i T_i(\xi) \quad [16] $$

The degrees of freedom are the components of the polynomial expansion, denoted $\mathbf{u}_j^i$ for the kinematics and $p_j^i$ for the pressure. Linearizing the equilibrium equations [12], yields the following system:

$$ \begin{aligned} \sum_{e}^{\text{nel}} \langle \Delta u_e, \Delta p_e \rangle & \left( \begin{bmatrix} [k_t] & [g] \\ [g]^t & -[m_p] \end{bmatrix} \right) \left\{ \begin{array}{c} \{\Delta u_e\} \\ \{\Delta p_e\} \end{array} \right\} \\ & - \left\{ \begin{array}{c} \{f_{\text{ext}}\} - \{r\} \\ -\{i\} \end{array} \right\} = 0 \end{aligned} \quad [17] $$
where nel is the number of elements in the mesh. The elementary matrices are defined by:

$$[k_t] = \int_{\Omega_e} [B]^t \left[ \frac{\partial \mathbf{\Pi}}{\partial \mathbf{F}} \right] [B] dV \\ [g] = \int_{\Omega_e} \sum_{i=0}^{n_p} T_i \frac{\partial G(J)}{\partial J} \{[B]^t \mathbf{C} \mathbf{o} f \mathbf{F}\} \otimes \{N_p^e\}^T dV \quad [18]$$

$$[m_p] = \alpha \int_{\Omega_e} \sum_{i=0}^{n_p} \sum_{j=0}^{n_p} T_i T_j \{N_p^e\} \otimes \{N_p^e\} dV$$

The residual and external forces are:

$$\{r\} = \int_{\Omega_e} [B]^t \{\mathbf{\Pi}\} dV \\ \{f_{ext}\} = \int_{\Omega_e} \sum_{i=0}^{n_u} \sum_{j=0}^{n_u} T_i T_j \{N_u^e\} \{\mathbf{f}^j\} dV + \int_{\partial \Omega_e^F} \sum_{i=0}^{n_u} \sum_{j=0}^{n_u} T_i T_j \{N_u^e\} \{\mathbf{F}^j\} dS \\ \{i\} = \int_{\Omega_e} \sum_{i=0}^{n_p} T_i \{N_p^e\} (G(J) - \frac{1}{k} p_e) dV \quad [19]$$

where $[B]$ is the linear operator relating the variation of the deformation gradient to the degrees of freedom: $\delta \mathbf{F} = [B]\delta \mathbf{u}$. The pressure degrees of freedom of each element are eliminated at elementary level by performing a static condensation step. This condensation of the pressure, gives a global number of dof per node define by:

$$n_{dof} = \dim \times (n_u + 1) \quad [20]$$

where dim is the physical dimension of the initial problem.

These developments have been implemented using the finite element software ZerbuLoN (Foerch et al. 1996).

## 4. Validation of the reduced model

In this section, we present some examples tests carried out to check the validity of the present reduced models. Shear and flexion tests were performed on a plane strain layer. This layer, 1 mm thick and 50 mm long, is composed of a Mooney-Rivlin materials which has the following characteristics: $a_{10} = 0.31$ Mpa, $a_{01} = 0.11$ Mpa and $k = 2500$ Mpa. In all these tests, the bottom part of the layer was fixed and the upper part was constrained so as to obtain rigid body motions. In this way, only one degree of freedom of the upper part was activated in each test. The occurrence of coupled motions was reflected in the resulting loads. The standard 2-D model is composed of Q9P3 elements (quadratic in displacement, linear in pressure) and the
<figure><img src="image_4.png" /></figure>

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

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

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

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

<figure><img src="image_9.png" /><figcaption>Figure 2. Flexion and shear test on plane strain layer</figcaption></figure>

reduced one, of L3P2 elements (quadratic in displacement, linear in pressure). The same mechanical formulation and the same software is used to compare the models.

Figures 2 show that good agreement was obtained between the reduced and standard models in both the shear and the flexion tests. It can be seen from figure 2(b) that the coupling between transversal and normal motions required greater enrichment of the polynomial basis than the other types of coupling. The other tests carried out showed good agreement with the reduced model in terms of local behavior: for exam-
ple, the boundary effects occurring near the free edges were accurately modelled in terms of the stress values (Lejeunes et al. 2005).

## 5. Buckling and post-buckling analysis

In studies on the non-linear responses of structures, three main problems tend to arise. First, the primary equilibrium path has to be followed, secondly, the singular points have to be detected and calculated, and thirdly, in the case of bifurcation points, the secondary path has to be followed. Nowadays, several approaches can be used to solve these three main problems. Incremental-iterative strategies with arc-length control are the most popular, see (Crisfield 1997, Riks 2004). Alternative methods such as the asymptotic numerical method have also been proposed, see (Damil et al. 1990).

In this study, we used a pseudo arc-length method, in line with Riks, to follow an equilibrium path. By expressing the equilibrium equation [12] as:

$$ \mathbf{R}(\mathbf{a}) - \lambda \mathbf{P} = 0 \qquad [21] $$

where a denotes the vector of nodal displacements, $\mathbf{R}(\mathbf{a})$ the internal forces, $\mathbf{P}$ some constant load pattern and $\lambda$ a scalar load factor, the pseudo arc-length method consist in the resolution of the following augmented system:

$$ \left[ \begin{array}{c} \mathbf{R}(\mathbf{a}) - \lambda \mathbf{P} \\ f(\Delta\mathbf{a}, \Delta\lambda) = \Delta\mathbf{a}^T \Delta\mathbf{a} - dl^2 \end{array} \right] = 0 \qquad [22] $$

where $\Delta\mathbf{a}$ relates to the lastest converged equilibrium states, $dl$ is the given initial radius of the constraint. During the corrector process, this procedure is linearized (Riks 1984).

The pseudo arc-length control is suitable to follow equilibrium path that exhibit snap-through or limit points. To detect a singular points, we use the smallest pivot of the tangent matrix $\mathbf{K}_T$, for this purpose we define:

$$ \tau = \min \{\text{pivot } \mathbf{K}_T\} \qquad [23] $$

A singular point is crossed if $\tau$ makes a changes of sign.

To analyze the characteristics of the singular point, we need to compute it precisely. This was done with an iterative bracketing procedure, which consists in performing a linear interpolation of the load factor between two points straddling the bifurcation point, such that:

$$ \lambda_i = \lambda_{i-2} - \tau_{i-2} \frac{\lambda_{i-1} - \lambda_{i-2}}{\tau_{i-1} - \tau_{i-2}} \qquad [24] $$

where $\lambda_{i-2}, \lambda_{i-1}, \tau_{i-2}, \tau_{i-1}$ are the control parameter and the value of the smallest pivot at the latest solutions that straddle the bifurcation point. More elaborated procedures exist for the determination of singular point see for example (Magnusson et al. 1998).
When the singular point has been carefully obtained, we compute Z the lowest eigenmode (corresponding to the lowest eigenvalue) at the singular point. We are now able to determine whether this point is a bifurcation point or a limit point by computing:

$$ P^T Z = 0 \quad \text{bifurcation point} \qquad [25] $$

$$ P^T Z \neq 0 \quad \text{limit point} \qquad [26] $$

In the case of a bifurcation point, we could use the bifurcation equation (see Crisfield (1997), section 20.3) to determine the nature of the bifurcation point, and in the case of a simple bifurcation to exactly compute the two tangent. However, this procedure involves calculating higher order derivatives of the equilibrium equation. This is generally not an easy task and alternative simplified procedures have been proposed. For instance, the step at the bifurcation point is initiated by taking the predictor in the direction of the eigenvector Z:

$$ \Delta a = dlZ \quad \Delta \lambda = 0 \qquad [27] $$

This procedure is exact in the case of a pitchfork bifurcation, as it often occurs in structural analysis. For a transcritical bifurcation, it could still be used as an approximated predictor (Crisfield 1997).

## 6. Application

As an example of how the reduced model can be applied to stability analysis, we propose to determine the buckling load and the post buckling behavior of a plane strain laminated bearing. This structure is composed of 10 layers: the rubber layers are 1 mm thick and the steel reinforcements are 1 mm thick. All the layers are 10 mm wide and the steel material is modeled as an elastic medium with E=210000 MPa and ν=0.3. The rubber is modeled as a Mooney-Rivlin material, such that $a_{10} = 0.31$ Mpa, $a_{01} = 0.11$ Mpa and the compressibility modulus is $k = 2500$ Mpa. We study the response of the bearing in two loading cases. In the first case, the bottom layer is fixed and the upper layer is subjected to a compression load, this layer can only move in the vertical direction. In the second case, the bottom layer is still fixed and the compression load is applied to the upper layer which can move in both the vertical and transversal directions.

Figure 3 (on which, for simplicity, the beginning of the primary path is not shown) allows to compare the responses obtained with the reduced and standard models. Very good agreement is observed between the buckling loads. As far as the post-buckling behavior is concerned, a good agreement is obtained in the second boundary case, and a softening response of the reduced model is obtained in the first boundary case. In terms of buckling modes, the reduced and standard models both give the same response, as can be seen from figures 4.
<figure><img src="image_10.png" /><figcaption>**Figure 3.** Equilibrium paths in the two boundary case</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 4.** Buckling modes in the fixed/fixed case and the fixed/free in transversal motion case</figcaption></figure>

To illustrate the advantages of reduced models, table 1 compares the reduced and standard models in terms of model size and the computing time required for the stability analysis (in the first boundary case). The corresponding computations have been performed on an Athlon MP processor (1533 Mhz) with a 768-Mb memory.

<table><caption>**Table 1.** Model size and cpu time required by the reduced and standard models for the laminated bearing</caption><thead><tr><th></th><th>dof</th><th>cpu time</th><th>gain</th></tr></thead><tbody><tr><td>2-D standard</td><td>10282</td><td>154</td><td></td></tr><tr><td>2-D reduced n<sub>u</sub> = 9, n<sub>p</sub> = 8</td><td>1494</td><td>12</td><td>12.8</td></tr></tbody></table>
## 7. Conclusion

In this paper, we have applied a model reduction method accounting for non-linear behavior under finite strain conditions, which can be used to simulate the behavior of laminated bearing structures. This study shows the reliability of this method in the context of an incremental-iterative procedure for performing buckling and post-buckling analysis on rubber bearing structures.

The reduced method was compared with standard plane strain finite element models, and proved to be advantageous in terms of both the computing time (which was nearly 13 fold shorter) and the model size (which was nearly 7 fold smaller). This method is an useful alternative to the classical finite element method, which involves either prohibitively large numbers of dof, or elements with a high aspect ratio. Other kinds of reduced finite element approaches have been developed elsewhere based on these concepts, which are suitable adapted for the modeling of three-dimensional beam-like structures (3D-2D reduction) and 3-D rectangular shape laminated bearings (3D-1D reduction), see (Lejeunes et al. 2005).

## 8. References

Boukamel A., Etude théorique et expérimentale d'un stratifié caoutchouc-acier en grandes déformations, Thèse de doctorat, Université d'Aix-Marseille II, 1988.

Cheung Y., Kong J., “The application of a new finite strip to the free vibration of rectangular plates of varying complexity”, *Journal of Sound and Vibration*, vol. 181, p. 341-353, 1995.

Crisfield M., *Non-linear Finite Element Analysis of Solids and Structures*, vol. 2, Wiley, 1997.

Damil N., M. P.-F., “A new method to compute perturbed bifurcations: Application to the buckling of imperfect elastic structures”, *International Journal of Engineering and Sciences*, vol. 28, p. 943-957, 1990.

Devries F., “Homogenization of elastomer matrix composites: method and validation”, *Composites Part B: Engineering*, vol. 29, p. 753-762, 1998.

Dumontet H., Homogénéisation et effets de bords dans les matériaux composites, Thèse d'état, Université Pierre et Marie Curie Paris 6, 1990.

Foerch R., Besson J., Cailletaud G., Pivlin P., “Polymorphic constitutive equations in finite element codes”, *Comput. Methods Appl. Mech. Engrg.*, vol. 141, p. 355-372, 1996.

Fu Y., Ogden R. (eds), *Nonlinear Elasticity: Theory and Applications*, Cambridge university press, 2001. London Mathematical Society Lecture Note Series, 283.

Glowinski R., Le Tallec P., “Numerical solution of problems in incompressible finite elasticity by augmented Lagrangian methods I. two-dimensional and axisymmetric problems”, *SIAM J. Appl. Math.*, vol. 42, p. 400-429, 1982.

Hartmann S., Neff P., “Polyconvexity of generalized polynomial-type hyperelastic strain energy functions for near incompressibility”, *International Journal of Solids and Structures*, vol. 40, p. 2767-2791, 2003.

Herrmann L. R., “Elasticity equations for incompressible and nearly-incompressible materials by a variational theorem”, *AIAA J.*, vol. 3, p. 1896-1900, 1963.

Holzapfel G., *Nonlinear Solid Mechanics*, Wiley, 2004.
Iizuka M., “A macroscopic model for predicting large-deformation behaviors of laminated rubber bearings”, *Engineering Structures*, vol. 22, p. 323-334, 2000.

Kelly J., “Tension Buckling in Multilayer Elastomeric Bearings”, *16th Engineering Mechanics Conference*, American Society of Civil Engineers, Seattle, July, 2003.

Lejeunes S., Boukamel A., Cochelin B., “Model reduction method for composites structures with elastomeric matrix”, in Austrell, Kari (eds), *Constitutive Models for Rubber IV*, Taylor & Francis Group., London, p. 391-396, 2005.

Léné F., Rey C., “Some strategies to compute elastomeric lamified composite structures”, *Compositive Structures*, vol. 54, p. 231-241, 2001.

Magnusson A., Svensson I., “Numerical treatment of complete load-deflection curves”, *Int. J. Numer. Meth. Engng.*, vol. 41, p. 955-971, 1998.

Malkus D., Hughes T., “Mixed finite element methods - reduced and selective integration techniques: a unification of concepts”, *Comp. Meth. Appl. Mech. Eng.*, vol. 15, p. 63-81, 1978.

Miehe C., “Aspects of the formulation and finite element implementation of large strain isotropic elasticity”, *International Journal for Numerical Methods In Engineering*, vol. 37, p. 1981-2004, 1994.

Mooney M., “A theory of large elastic deformation”, *J. Appl. Phys.*, vol. 11, p. 582-592, 1940.

Riks E., “Some computational aspects of the stability analysis of nonlinear structures”, *Computers Methods in Applied Mechanics and Engineering*, vol. 47, p. 219-259, 1984.

Riks E., “Buckling”, in E. Stein, R. De Borst, T. J. R. Hughes (eds), *Encyclopedia of Computational Mechanics*, vol. 2, Wiley, p. 139-167, 2004.

Rüter M., Stein E., “Analysis, finite element computation and error estimation in transversely isotropic nearly incompressible finite elasticity”, *Comput. Methods Appl. Mech. Engrg.*, vol. 190, p. 519-541, 2000.

Zhong W., Cheung Y., Li Y., “The precise finite strip method”, *Computers and Structures*, vol. 69, p. 773-783, 1998.