Model order reduction for hyperelastic materials

Siamak Niroomandi, Icíar Alfaro, Elías Cueto, Francisco Chinesta

► To cite this version:

Siamak Niroomandi, Icíar Alfaro, Elías Cueto, Francisco Chinesta. Model order reduction for hyper-elastic materials. International Journal for Numerical Methods in Engineering, 2010, 81 (9), pp.1180-1206. 10.1002/nme.2733 . hal-01007059

HAL Id: hal-01007059

https://hal.science/hal-01007059

Submitted on 2 Apr 2017

**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.

Public Domain
Model order reduction for hyperelastic materials

Siamak Niroomandi¹, Icıar Alfaro¹, Elías Cueto¹ and Francisco Chinesta²

¹Group of Structural Mechanics and Material Modelling, Aragón Institute of Engineering Research (I3A),
University of Zaragoza, María de Luna, 5. Campus Rio Ebro, E-50018 Zaragoza, Spain

²EADS Corporate International Chair, École Centrale de Nantes. 1, Rue de la Noë, 44300 Nantes, France

In this paper, we develop a novel algorithm for the dimensional reduction of the models of hyperelastic solids undergoing large strains. Unlike standard proper orthogonal decomposition methods, the proposed algorithm minimizes the use of the Newton algorithms in the search of non-linear equilibrium paths of elastic bodies.

The proposed technique is based upon two main ingredients. On one side, the use of classic proper orthogonal decomposition techniques, that extract the most valuable information from pre-computed, complete models. This information is used to build global shape functions in a Ritz-like framework.

On the other hand, to reduce the use of Newton procedures, an asymptotic expansion is made for some variables of interest. This expansion shows the interesting feature of possessing one unique tangent operator for all the terms of the expansion, thus minimizing the updating of the tangent stiffness matrix of the problem.

The paper is completed with some numerical examples in order to show the performance of the technique in the framework of hyperelastic (Kirchhoff–Saint Venant and neo-Hookean) solids.

KEY WORDS: model order reduction; proper orthogonal decomposition; asymptotic numerical method; Kirchhoff–Saint Venant; neo-Hookean material

# 1. INTRODUCTION

Model order reduction in computational mechanics has a strong tradition, starting from the early years in which the performance of computers implied the need for models with very few degrees of freedom, see [1, 2] and references therein, to name a few. But similar techniques are employed
in many other branches of Engineering and Applied Sciences: real-time simulation [3], Chemical Engineering [4], turbulence [5], weather forecast [6] and many others, linking back to the seminal works of Karhunen [7] and Loève [8].

Recent applications of model reduction techniques for the analysis of non-linear solids and structures include the works by Marsden and colleagues [2] and the works of Ryckelyck *et al*. for ‘*a priori*’ model reduction (see [9–11]). Also, recently, new extensions of the method have been developed, giving rise to the so-called *non-linear dimensionality reduction* [12].

In general, all these techniques share some common characteristics. It is necessary to have some previous results coming from previously computed, detailed models (in this frameworks, the work by Ryckelyck [11] is an exception). From this information, the most relevant structures of these results are computed by minimizing a functional that measures the distance from the data to the sought structures. This constitutes the so-called Karhunen–Loève (K–L) method, also known as proper orthogonal decomposition (POD) or principal component analysis (PCA). Finally, the reduced models are constructed by employing a reduced set of global, ‘good quality’ (in statistical terms) basis functions, as opposed to the FEM that employs compactly supported, general-purpose, piecewise polynomial shape functions.

Models that differ slightly from the original ones can thus be computed with a very important, sometimes impressive, computational saving, employing only a very limited number of degrees of freedom. This is, for instance, the approach employed by the authors in previous works, in which real-time simulation speed was sought [13].

For highly non-linear applications, however, these methods posses some limitations. The radius of convergence of the before-mentioned techniques is often low, unless some treatment is performed. For instance, it is often necessary to update the complete stiffness matrix in order to follow non-linear equilibrium paths, and this is sometimes burdensome or even impossible (as in real-time applications, for instance, in which results must be obtained at a frequency of 500 Hz [13]). Otherwise, the resulting simulation will be actually linear, even if it employs the best (in statistical terms) basis functions available. In the work by Ryckelycnk [11], on the contrary, Krylov subspaces are added to the reduced set of basis functions once the simulation provides an unsatisfactory norm of the residual. But this implies to come back to the last converged time step and to continue with the just enriched basis.

The approach followed in this paper is somewhat different. In order to avoid a frequent update of the stiffness matrix, an asymptotic expansion is made at the last converged time step of some variables of interest. When this expansion is introduced into the discrete weak form of the problem, a series of problems is encountered that share the same stiffness matrix for all the terms of the expansion. Of course, the expansion has a finite radius of convergence, but it is frequently enough to perform only a few updates of the stiffness matrix for complex equilibrium paths.

This technique thus combines two main ingredients: the proper orthogonal decomposition of the existing data coming from complete models, and an asymptotic expansion of variables of interest in the neighbourhood of equilibrium points [14, 15]. This technique was developed in [16] for the buckling analysis of foam structures, considered linear elastic (geometrically non-linear), although is generalized here for more complex, hyperelastic, constitutive equations including material non-linearities. This extension is not straightforward and, up to the best of our knowledge, has not been developed so far.

The structure of the paper is as follows. In Section 2 we review the basics of the technique proposed in [16] for linear elastic solids undergoing finite strains (Kirchhoff–Saint Venant models). In Section 3 we extend this technique to account for material non-linearities also. The proposed
technique is described in detail and, finally, in Section 4 we include some examples of application for Kirchhoff-Saint Venant and neo-Hookean models that show the performance of the technique and compare the computational cost with a standard FE solution.

## 2. STATE-OF-THE-ART: MODEL REDUCTION AND ASYMPTOTIC EXPANSIONS

### 2.1. The Karhunen-Loève decomposition

We assume that the evolution of a certain field $u(\mathbf{x}, t)$, governed by a PDE, is known. In practical applications, this field is expressed in a discrete form, that is, it is known typically at the nodes of a spatial mesh and for some time steps of existing simulations $u(\mathbf{x}_i, t^n) \equiv u_i^n$. The same can be wrote introducing a time discretization $u^n(\mathbf{x}) \equiv u(\mathbf{x}, t=n\Delta t)$; $\forall n \in [1, \dots, P]$. The main idea of the Karhunen-Loève (K-L) decomposition is how to obtain the most typical or characteristic structure $\phi(\mathbf{x})$ among these $u^n(\mathbf{x}) \forall n$. This is equivalent to obtaining a function $\phi(\mathbf{x})$ that maximizes $\alpha$ defined by

$$ \alpha = \frac{\sum_{n=1}^{n=P} \left[ \sum_{i=1}^{i=N} \phi(\mathbf{x}_i) u^n(\mathbf{x}_i) \right]^2}{\sum_{i=1}^{i=N} (\phi(\mathbf{x}_i))^2} \quad (1) $$

The maximization ($\delta\alpha=0$) leads to

$$ \sum_{n=1}^{n=P} \left[ \left( \sum_{i=1}^{i=N} \tilde{\phi}(\mathbf{x}_i) u^n(\mathbf{x}_i) \right) \left( \sum_{j=1}^{j=N} \phi(\mathbf{x}_j) u^n(\mathbf{x}_j) \right) \right] = \alpha \sum_{i=1}^{i=N} \tilde{\phi}(\mathbf{x}_i) \phi(\mathbf{x}_i), \quad \forall \tilde{\phi} \quad (2) $$

which can be rewritten in the form

$$ \sum_{i=1}^{i=N} \left\{ \sum_{j=1}^{j=N} \left[ \sum_{n=1}^{n=P} u^n(\mathbf{x}_i) u^n(\mathbf{x}_j) \phi(\mathbf{x}_j) \right] \tilde{\phi}(\mathbf{x}_i) \right\} = \alpha \sum_{i=1}^{i=N} \tilde{\phi}(\mathbf{x}_i) \phi(\mathbf{x}_i), \quad \forall \tilde{\phi} \quad (3) $$

Defining the vector $\mathbf{\phi}$ such that its $i$-component is $\phi(\mathbf{x}_i)$, Equation (3) takes the following matrix form:

$$ \tilde{\mathbf{\phi}}^{\mathrm{T}} \mathbf{c} \mathbf{\phi} = \alpha \tilde{\mathbf{\phi}}^{\mathrm{T}} \mathbf{\phi}, \quad \forall \tilde{\mathbf{\phi}} \Rightarrow \mathbf{c} \mathbf{\phi} = \alpha \mathbf{\phi} \quad (4) $$

where the two-point correlation matrix is given by

$$ c_{ij} = \sum_{n=1}^{n=P} u^n(\mathbf{x}_i) u^n(\mathbf{x}_j) \Leftrightarrow \mathbf{c} = \sum_{n=1}^{n=P} \mathbf{u}^n(\mathbf{u}^n)^{\mathrm{T}} \quad (5) $$

which is symmetric and positive definite. If we define the matrix $\mathbf{Q}$ containing the discrete field history:

$$ \mathbf{Q} = \begin{pmatrix} u_1^1 & u_1^2 & \cdots & u_1^P \\ u_2^1 & u_2^2 & \cdots & u_2^P \\ \vdots & \vdots & \ddots & \vdots \\ u_N^1 & u_N^2 & \cdots & u_N^P \end{pmatrix} \quad (6) $$
it is straightforward to verify that the matrix **c** in Equation (4) results

$$ \mathbf{c} = \mathbf{Q}\mathbf{Q}^{\mathrm{T}} \qquad (7) $$

where the diagonal components are given by

$$ c_{ii} = (\mathbf{Q}\mathbf{Q}^{\mathrm{T}})_{ii} = \sum_{j=1}^{j=P} (u_i^j)^2 \qquad (8) $$

Thus, the functions defining the most characteristic structure of $u^n(\mathbf{x})$ are the eigenfunctions $\phi_k(\mathbf{x}) \equiv \phi_k$ associated with the highest eigenvalues.

## 2.2. A posteriori reduced-order modelling

If some direct simulations have been carried out, we can determine $u(\mathbf{x}_i, t^n) \equiv u_i^n$, $\forall i \in [1, \dots, N]$, $\forall n \in [1, \dots, P]$, and from these the $r$ eigenvectors related to the $r$-highest eigenvalues $\phi_k = \phi_k(\mathbf{x}_i)$, $\forall i \in [1, \dots, N]$, $\forall k \in [1, \dots, r]$ (with $r \ll N$). Now, we can try to use these $r$ eigenfunctions for approximating the solution of a problem slightly different to the one that has served to define $u(\mathbf{x}_i, t^n)$. For this purpose we need to define the matrix **A**

$$ \mathbf{A} = \begin{pmatrix} \phi_1(\mathbf{x}_1) & \phi_2(\mathbf{x}_1) & \cdots & \phi_r(\mathbf{x}_1) \\ \phi_1(\mathbf{x}_2) & \phi_2(\mathbf{x}_2) & \cdots & \phi_r(\mathbf{x}_2) \\ \vdots & \vdots & \ddots & \vdots \\ \phi_1(\mathbf{x}_N) & \phi_2(\mathbf{x}_N) & \cdots & \phi_r(\mathbf{x}_N) \end{pmatrix} \qquad (9) $$

Now, we consider the linear system of equations resulting from the discretization of a partial differential equation (PDE) in the form

$$ \mathbf{KU} = \mathbf{F} \qquad (10) $$

Obviously, in the case of evolution problems **F** contains the contribution of the solution at the previous time step.

Then, assuming that the unknown vector contains the nodal degrees of freedom, it can be expressed as

$$ \mathbf{U} = \sum_{i=1}^{i=r} \zeta_i \phi_i = \mathbf{A} \zeta \qquad (11) $$

it results

$$ \mathbf{KU} = \mathbf{F} \Rightarrow \mathbf{K} \mathbf{A} \zeta = \mathbf{F} \qquad (12) $$

and multiplying both terms by $\mathbf{A}^T$ it results

$$ \mathbf{A}^{\mathrm{T}} \mathbf{K} \mathbf{A} \zeta = \mathbf{A}^{\mathrm{T}} \mathbf{F} \qquad (13) $$

which proves that the final system of equations is of low order, that is, the dimensions of $\mathbf{A}^T \mathbf{K} \mathbf{A}$ are $r \times r$, with $r \ll N$, and the dimensions of both $\zeta$ and $\mathbf{A}^T \mathbf{F}$ are $r \times 1$.
### Remark 1

Equation (13) can also be derived by introducing the approximation (11) into the PDE Galerkin form.

As can be noticed from the paragraphs above, the main idea of model order reduction techniques is to employ, in a Ritz framework, a set of global basis that are, in a statistical sense, the best suited to reproduce the just computed, complete models. This is in sharp contrast with the finite element method that employs general purpose, piecewise polynomials to approximate the solution in a Galerkin framework.

## 2.3. Limitations of standard model reduction techniques

Classical model reduction techniques, when applied to the simulation of non-linear solids and structures, imply the need for frequent updating of the stiffness tangent matrix $K$ (or, equivalently, the reduced stiffness matrix $A^T K A$, but the matrix $A$ is assumed constant throughout the simulation). Otherwise, the reduced model will obviously be linear. This implies that once the residual of the discrete, algebraic form of the problems is estimated unacceptable, the only way of searching the equilibrium is by updating the stiffness matrix of the complete problem. Modified Newton-Raphson methods could equally be employed, but in the context of reduced models their convergence is often judged too slow. Even if this strategy can be employed for reducing non-linear models the computing cost is considerable, and of course this technique is unviable if one is looking for quasi-real-time simulations.

A different question arises in the so-called hyper-reduction methods [17] in which only a few rows of the stiffness matrix are integrated, based on the fact of the limited number of degrees of freedom that POD techniques select for a particular problem. The problem is precisely to select appropriately in which part of the model the integration is performed.

This stiffness matrix updating is usually an expensive procedure in terms of computational cost (and noteworthy in the context of reduced models). But sometimes is simply unaffordable. This is the case in the framework of real-time or near real-time simulations, where a frequency of 500-1000 Hz in the response of the simulation is needed, for instance, for haptic realism (30 Hz would be enough for video feedback, which is too much a limitation indeed).

## 2.4. Asymptotic numerical methods

Non-linear structural problems are generally solved using iterative methods, such as Newton-Raphson or modified Newton schemes. Such algorithms are successful for solving the non-linear equilibrium equations of the model. However, the computing time is usually large as compared with a linear solution, due to both the number of iterations usually needed and the computation of tangent stiffness matrices. A family of asymptotic numerical methods (ANM) based on perturbation techniques and finite element methods have been proposed and intensively used by Potier-Ferry and co-workers for computing perturbed bifurcations, and applied in computing the post-buckling behaviour of elastic plates and shells. Next, they have extended the method to many non-linear elastic solutions, plastic deformations, etc. For a complete review the interested reader can refer to [14, 15]. In contrast to predictor-corrector algorithms, the non-linear equilibrium paths are determined by means of asymptotic expansions: the unknown nodal vector $U$ and the load parameter $\bar{\lambda}$—that represents a pseudo-time, the overline is used to avoid confusion with the Lamé coefficient—are represented by power series expansions with respect to a control parameter $a$.
By introducing the expansions into the equilibrium equation, the non-linear problem is transformed into a sequence of linear problems in a recurrent manner and is solved by the finite element method, for instance. Noteworthy, because all the linear problems have the same stiffness matrix, the method requires only one matrix inversion. Moreover, one gets a continuous analytic representation of the load-displacement (response) curve that differs from the point by point representation of standard algorithms.

It is precisely this feature (the existence of only one tangent operator per equilibrium point) that makes the ANM an appealing choice to combine with model reduction techniques in order to get a very efficient technique.

## 2.5. Problem formulation for Kirchhoff–Saint venant models

Here, we review the development made in [16] for the buckling analysis of foam structures that we will also considered in our numerical examples, where they will be compared with the constitutive non-linear models. We consider, as usual, a Lagrangian description of the movement whose material coordinates are given by the vector **X**. The solid occupies a volume $\Omega_0$ whose boundary is denoted by $\Gamma_u$. Essential and natural boundary conditions are applied to the non-overlapping portions of the boundary denoted by $\Gamma_u$ and $\Gamma_t$. The deformed configuration of each point is given by the vector **x**, such that

$$ \mathbf{x} = \mathbf{X} + \mathbf{u} \qquad (14) $$

where the displacement field is denoted by **u**.

Following [18], we consider a linear and a non-linear terms for the Green–Lagrange strain tensor, **E**, in the form

$$ \mathbf{E} = \frac{1}{2}(\mathbf{F}^T\mathbf{F} - 1) = \gamma_l(\mathbf{u}) + \gamma_{nl}(\mathbf{u}, \mathbf{u}) \qquad (15) $$

where $\mathbf{F} = \nabla\mathbf{u} + \mathbf{I}$ is the gradient of the deformation tensor and using the notation in [18],

$$ \gamma_l(\mathbf{u}) = \frac{1}{2}(\nabla(\mathbf{u}^T) + \nabla(\mathbf{u})) \qquad (16a) $$

$$ \gamma_{nl}(\mathbf{u}, \mathbf{u}) = \frac{1}{2} \nabla(\mathbf{u}^T) \nabla(\mathbf{u}) \qquad (16b) $$

Hyperelastic materials are based on the assumption of a particular strain-energy function, $\Psi$. Then the second Piola–Kirchhoff stress tensor **S** can thus be obtained by

$$ \mathbf{S} = \frac{\partial \Psi}{\partial \mathbf{E}} \qquad (17) $$

that is a symmetric tensor and is related to the first Piola–Kirchhoff stress tensor, **P**, by $\mathbf{P} = \mathbf{F}\mathbf{S}$.

The equilibrium equation stated in the reference configuration looks like

$$ \nabla \mathbf{P} + \mathbf{B} = \mathbf{0} \quad \text{in } \Omega_0 \qquad (18) $$

in which **B** is the body force. The boundary conditions of the body are defined by

$$ \begin{aligned} \mathbf{u}(\mathbf{X}) &= \bar{\mathbf{u}} & \text{on } \Gamma_u \\ \mathbf{P}\mathbf{N} &= \bar{\lambda}\bar{\mathbf{t}} & \text{on } \Gamma_t \end{aligned} \qquad (19) $$
where N is the unit vector normal to Γ, $\bar{\mathbf{t}}$ is an applied traction and $\bar{\lambda}$ is a loading parameter, equivalent to a pseudo-time, and ranging from 0 to 1. The weak form of the problem is then given by

$$ \int_{\Omega_0} \mathbf{S} : \delta \mathbf{E} \, d\Omega = \bar{\lambda} \int_{\Gamma_t} \bar{\mathbf{t}} \cdot \delta \mathbf{u} \, d\Gamma \quad \forall \delta \mathbf{u} \in H^1(\Omega) \qquad (20) $$

where in the above equation $\delta \mathbf{E}$ is expressed by

$$ \delta \mathbf{E} = \frac{1}{2} [\mathbf{F}^T \nabla(\delta \mathbf{u}) + \nabla(\delta \mathbf{u})^T \mathbf{F}] = \gamma_l(\delta \mathbf{u}) + \gamma_{nl_S}(\mathbf{u}, \delta \mathbf{u}) \qquad (21) $$

where, in turn, $\gamma_{nl_S}(\mathbf{u}, \delta \mathbf{u})$ is defined by

$$ \gamma_{nl_S}(\mathbf{u}, \delta \mathbf{u}) = \gamma_{nl}(\mathbf{u}, \delta \mathbf{u}) + \gamma_{nl}(\delta \mathbf{u}, \mathbf{u}) \qquad (22) $$

The technique is best viewed by considering the Saint-Venant Kirchhoff model, as done in [16, 18], since only geometric non-linearities are present. Later on, in Section 3, we will generalize the method for other hyperelastic models. The Kirchhoff–Saint Venant model is characterized by the energy function given by

$$ \Psi = \frac{\lambda}{2} (\mathrm{tr}(\mathbf{E}))^2 + \mu \mathbf{E} : \mathbf{E} \qquad (23) $$

where $\lambda$ and $\mu$ are Lame's constants. The second Piola–Kirchhoff stress tensor can be obtained by

$$ \mathbf{S} = \frac{\partial \Psi(\mathbf{E})}{\partial \mathbf{E}} = \mathbf{C} : \mathbf{E} \qquad (24) $$

in which $\mathbf{C}$ is the fourth-order constitutive (elastic) tensor.

The Saint Venant–Kirchhoff model possesses well-known limitations, particularly some instabilities when subjected to pure compression. Nevertheless, it remains to be interesting for some applications. Noteworthy, in real-time simulation environments (see [3] and references therein, for instance) it is among the state-of-the-art models that can be computed under the severe limitations that real-time frameworks impose (30 Hz for video feedback and 500 Hz for haptic feedback).

## 2.6. Asymptotic numerical method for geometrically non-linear problems

In the ANM [14, 15] the displacement of each material point is expanded asymptotically in terms of a control parameter 'a'. This expansion is developed in the neighbourhood of a known equilibrium point $(\mathbf{u}^n; \mathbf{S}^n; \bar{\lambda}^n)$ at step n and the series is truncated at order N. To simplify the resulting expressions, also the second Piola–Kirchhoff stress tensor and the load parameter $\lambda$ are expanded in series prior to their introduction in the weak form of the problem

$$ \left\{ \begin{matrix} \mathbf{u}^{n+1}(a) \\ \mathbf{S}^{n+1}(a) \\ \bar{\lambda}^{n+1}(a) \end{matrix} \right\} = \left\{ \begin{matrix} \mathbf{u}^n(a) \\ \mathbf{S}^n(a) \\ \bar{\lambda}^n(a) \end{matrix} \right\} + \sum_{p=1}^{N} a^p \left\{ \begin{matrix} \mathbf{u}_p \\ \mathbf{S}_p \\ \bar{\lambda}_p \end{matrix} \right\} \qquad (25) $$

where $(\mathbf{u}_p, \mathbf{S}_p, \bar{\lambda}_p)$ are unknowns. Above, $(\mathbf{u}^{n+1}(a), \mathbf{S}^{n+1}(a), \bar{\lambda}^{n+1}(a))$ represents the solution along a portion of the loading curve. Noteworthy, the behaviour of the solid is described continuously
with respect to 'a'. The introduction of Equation (25) into Equations (20) and (24) leads to a series of linear problems with the same tangent operator, thus avoiding the burden associated with stiffness matrix updating in the Newton-Raphson scheme.

A general procedure for constructing quadratic forms of equations consists in developing auxiliary variables into finite series, as will be detailed in Section 3 for neo-Hookean materials. The series expansion of $\delta E(\mathbf{u})$ gives

$$ \delta E^{n+1}(a) = \gamma_l(\delta\mathbf{u}) + \gamma_{nl_S}(\delta\mathbf{u}, \mathbf{u}^n) + \sum_{p=1}^{n} a^p \gamma_{nl_S}(\delta\mathbf{u}, \mathbf{u}_p) \quad (26) $$

The series expansions of $S$ gives in turn

$$ \begin{aligned} S^{n+1}(a) &= C : E^{n+1}(a) \\ &= C : \left[ \gamma_{nl}(\mathbf{u}^n, \mathbf{u}^n) + \gamma_l(\mathbf{u}^n) + \sum_{p=1}^{n} a^p \left( \gamma_l(\mathbf{u}_p) + \gamma_{nl_S}(\mathbf{u}^n, \mathbf{u}_p) + \sum_{i=1}^{p-1} \gamma_{nl}(\mathbf{u}_i, \mathbf{u}_{p-i}) \right) \right] \end{aligned} \quad (27) $$

and at order $p$ we obtain

$$ S_p = C : \left\{ \gamma_l(\mathbf{u}_p) + \gamma_{nl_S}(\mathbf{u}^n, \mathbf{u}_p) + \sum_{i=1}^{p-1} \gamma_{nl}(\mathbf{u}_i, \mathbf{u}_{p-i}) \right\} \quad (28) $$

Introducing the asymptotic expansion into Equation (24) results in

$$ \begin{aligned} & \int_{\Omega_0} \left\{ \left( S^n + \sum_{p=1}^{N} a^p S_p \right) : \left( \gamma_l(\delta \mathbf{u}) + \gamma_{nl_S}(\mathbf{u}^n, \delta \mathbf{u}) + \sum_{p=1}^{N} a^p \gamma_{nl_S}(\mathbf{u}_p, \delta \mathbf{u}) \right) \right\} d\Omega \\ &= \left( \bar{\lambda}^n + \sum_{p=1}^{N} a^p \bar{\lambda}_p \right) \Psi_{\text{ext}}(\delta \mathbf{u}) \end{aligned} \quad (29) $$

with $\Psi_{\text{ext}}(\delta\mathbf{u}) = \int_{\Gamma_l} \mathbf{t} \cdot \delta\mathbf{u} \, d\Gamma$. Introducing Equation (28) into Equation (29) and identifying terms with the same power of 'a' results in a successive series of linear problems which at order $p$, ($p=1, \dots, N$) takes the form

$$ \mathcal{L}(\delta\mathbf{u}, \mathbf{u}^n) = \bar{\lambda}_p \Psi_{\text{ext}}(\delta\mathbf{u}) + F_p^{nl}(\delta\mathbf{u}) \quad (30) $$

with

$$ \mathcal{L}(\delta\mathbf{u}, \mathbf{u}^n) = \int_{\Omega} \{ S^n : \gamma_{nl_S}(\mathbf{u}^n, \delta\mathbf{u}) + [\gamma_l(\delta\mathbf{u}) + \gamma_{nl_S}(\mathbf{u}_p, \delta\mathbf{u})] : C : [\gamma_l(\mathbf{u}_p) + \gamma_{nl_S}(\mathbf{u}^n, \mathbf{u}_p)] \} d\Omega \quad (31) $$

and where $F_p^{nl}(\delta\mathbf{u})$ is equal to zero at order one and at order $p$ it can be calculated as

$$ F_p^{nl}(\delta\mathbf{u}) = - \int_{\Omega} \left\{ \sum_{i=1}^{p-1} S_i : \gamma_{nl_S}(\mathbf{u}_{p-i}, \delta\mathbf{u}) + \sum_{i=1}^{p-1} [\gamma_{nl}(\mathbf{u}_i, \mathbf{u}_{p-i})] : C : [\gamma_l(\delta\mathbf{u}) + \gamma_{nl_S}(\mathbf{u}^n, \delta\mathbf{u})] \right\} d\Omega \quad (32) $$
Discretization of Equation (30) by using finite elements leads to a sequence of linear problems in the form [18]

$$ \text{Order 1} \left\{ \begin{array}{l} \mathbf{K}_t \mathbf{u}_1 = \bar{\lambda}_1 \mathbf{f} \\ \mathbf{u}_1^T \mathbf{u}_1 + \bar{\lambda}_1^2 = 1 \end{array} \right. \qquad (33) $$

$$ \text{Order } p \left\{ \begin{array}{l} \mathbf{K}_t \mathbf{u}_p = \bar{\lambda}_p \mathbf{f} + \mathbf{f}_p^{nl}(\mathbf{u}_i), \quad i<p \\ \mathbf{u}_p^T \mathbf{u}_1 + \bar{\lambda}_p \lambda_1 = 0 \end{array} \right. \qquad (34) $$

where $\mathbf{K}_t$ denotes the tangent stiffness matrix associated with Equation (31), common to the problems at different orders $p$. It is the same as the one applied in a classical iterative algorithm like Newton-Raphson (in the first iteration). In the above, $\mathbf{u}_p$ is the discretized form of the displacement field at order $p$, $\mathbf{f}$ is the loading vector and $\mathbf{f}_p^{nl}$ represents the discretized form associated with $\mathbf{F}_p^{nl}(\delta\mathbf{u})$ in Equation (32), which at order $p$ only depends on the values of $\mathbf{u}_i$, $i<p$. The solution of these problems can be obtained as follows:

$$ \text{at order } 1 \left\{ \begin{array}{l} \hat{\mathbf{u}} = \{\mathbf{K}_t\}^{-1} \mathbf{f} \\ \bar{\lambda}_1 = \frac{1}{\sqrt{\hat{\mathbf{u}}^T \hat{\mathbf{u}} + 1}} \\ \mathbf{u}_1 = \bar{\lambda}_1 \hat{\mathbf{u}} \end{array} \right. \qquad (35) $$

$$ \text{at order } p \left\{ \begin{array}{l} \mathbf{u}_p^{nl} = \{\mathbf{K}_t\}^{-1} \mathbf{f}_p^{nl} \\ \bar{\lambda}_p = -\bar{\lambda}_1 \{\mathbf{u}_p^{nl}\}^T \mathbf{u}_1 \\ \mathbf{u}_p = \frac{\bar{\lambda}_p}{\bar{\lambda}_1} \mathbf{u}_1 + \mathbf{u}_p^{nl} \end{array} \right. \qquad (36) $$

## 2.7. Combined POD-ANM method

As explained before the use of the reduced-order basis within the standard strategy, without updating the stiffness matrix, leads to a very fast calculation of the system of equations but with some error. On the other hand, by using ANM we can obtain the solution of non-linear problems accurately in a neighbourhood of an equilibrium point, and without the need for iterative procedures. In [16] a combination of these two methods to solve material and geometrically non-linear problems was employed. In the next section, we extend this procedure for solving non-linear constitutive models. Here, we assume that the POD basis has been calculated as explained in Section 2.1. The terms of the asymptotic expansions associated with the displacements are in turn expressed as functions of POD basis as

$$ \mathbf{u}_p = \sum_{m=1}^{M} \Phi^m \zeta_p^m = \mathbf{A} \zeta_p \qquad (37) $$
where $\zeta_p$ are unknowns. So the new asymptotic expansion of $\mathbf{u}$ is expressed by

$$ \mathbf{u}^{n+1}(a) = \mathbf{A} \left( \zeta^n + \sum_{p=1}^{N} a^p \zeta_p \right) \qquad (38) $$

where $(\zeta^n, \bar{\lambda}^n)$ represents the previous converged solution. Introducing Equation (38) into Equation(30) results in

$$ \mathcal{L}(\mathbf{A}\zeta_p, \mathbf{A}\delta\zeta) = \bar{\lambda}_p \delta\Psi_{\text{ext}}(\mathbf{A}\delta\zeta) + F_p^{nl}(\mathbf{A}\delta\zeta) \qquad (39) $$

After finite element discretization, and making use of the fact that the POD basis is orthonormal, we have the following sequence of linear systems of equations:

$$ \text{At order } 1 \left\{ \begin{array}{l} \mathbf{A}^{\mathrm{T}} \mathbf{K}_t \mathbf{A} \zeta_1 = \bar{\lambda}_1 \mathbf{A}^{\mathrm{T}} \mathbf{f} \\ \zeta_1^{\mathrm{T}} \zeta_1 + \bar{\lambda}_1^2 = 1 \end{array} \right. \qquad (40) $$

$$ \text{At order } p \left\{ \begin{array}{l} \mathbf{A}^{\mathrm{T}} \mathbf{K}_t \mathbf{A} \zeta_p = \bar{\lambda}_p \mathbf{A}^{\mathrm{T}} \mathbf{f} + \mathbf{f}_p^{nl} \\ \zeta_p^{\mathrm{T}} \zeta_1 + \bar{\lambda}_p \bar{\lambda}_1 = 0 \end{array} \right. \qquad (41) $$

The size of the above equations depends on the number of the POD basis functions, but it is very low, as explained before (typically less than ten), so they can be computed in a very short time.

The solution of these equations can be obtained as follows:

$$ \text{At order } 1 \left\{ \begin{array}{l} \hat{\zeta} = \{\mathbf{A}^{\mathrm{T}} \mathbf{K}_t \mathbf{A}\}^{-1} \mathbf{A}^{\mathrm{T}} \mathbf{f} \\ \bar{\lambda}_1 = \frac{1}{\sqrt{\hat{\zeta}^{\mathrm{T}} \hat{\zeta} + 1}} \\ \zeta_1 = \bar{\lambda}_1 \hat{\zeta} \\ \hat{u} = \mathbf{A} \hat{\zeta}, \mathbf{u}_1 = \bar{\lambda}_1 \zeta_1 \end{array} \right. \qquad (42) $$

$$ \text{At order } p \left\{ \begin{array}{l} \zeta_p^{nl} = \{\mathbf{A}^{\mathrm{T}} \mathbf{K}_t \mathbf{A}\}^{-1} \mathbf{A}^{\mathrm{T}} \mathbf{f}_p^{nl} \\ \bar{\lambda}_p = -\bar{\lambda}_1 \{\zeta_p^{nl}\}^{\mathrm{T}} \zeta_1 \\ \zeta_p = \frac{\bar{\lambda}_p}{\bar{\lambda}_1} \zeta_1 + \zeta_p^{nl} \\ \mathbf{u}_p = \mathbf{A} \zeta_p \end{array} \right. \qquad (43) $$

Note that, as mentioned before, in Equations (42)–(43) there exists only one tangent stiffness matrix per time step, shared by all the terms for any order $p$. The advantage, of course, comes from the fact that the radius of convergence of the method, as will be seen, and consequently the size of time steps, is considerably bigger than in traditional Newton procedures.
### 3. NON-LINEAR MODEL REDUCTION OF HYPERELASTIC MATERIALS

The extension of the technique explained before to other hyperelastic materials, in which non-linearities other than geometrical ones are present, is not straightforward and, up to our knowledge it has not been done before. In this section we extend the technique to neo-Hookean materials [19]. Extension to other hyperelastic materials would follow the same guidelines, the major difference with the Kirchhoff–Saint Venant model being the presence of material non-linearities, in addition to the geometrical ones.

The compressible Neo-Hookean model is characterized by a strain energy function given by

$$ \Psi = \frac{\mu}{2}(\mathrm{tr}(\mathbf{C}) - 3) - \mu \ln J + \frac{\lambda}{2}(\ln J)^2 \quad (44) $$

where $\lambda$ and $\mu$ are Lame's constants and $\mathbf{C}=\mathbf{I}+2\mathbf{E}$ is the right Cauchy–Green strain tensor. The second Piola–Kirchhoff stress tensor can be obtained by

$$ \mathbf{S} = \frac{\partial \Psi(\mathbf{E})}{\partial \mathbf{E}} = \mu(\mathbf{I} - \mathbf{C}^{-1}) + \lambda(\ln J)\mathbf{C}^{-1} \quad (45) $$

In this case, an expansion similar to that in Equation (25) is done, but in this case the intricate expansion procedure becomes more clear if we identify, as in [20], the asymptotic expansion with a Taylor series of the variables of interest, denoted by $\mathbf{U}(a)$, in the vicinity of $a=0$. Truncating at order $N$:

$$ \mathbf{U}(a) = \mathbf{U}_0 + \sum_{p=1}^{N} \mathbf{U}_p a^p \quad (46) $$

where $\mathbf{U}_0 = \mathbf{U}(0)$ and

$$ \mathbf{U}_p = \left. \frac{1}{p!} \frac{d^p \mathbf{U}}{da^p} \right|_{a=0} \quad (47) $$

In this case we have selected the following variables to perform the expansion:

$$ \mathbf{U}(a) = \begin{pmatrix} \mathbf{u}(a) \\ \mathbf{S}(a) \\ J^2 \mathbf{C}^{-1}(a) \\ \ln \sqrt{J^2}(a) \\ \frac{1}{J^2}(a) \\ \bar{\lambda}(a) \end{pmatrix} \quad (48) $$

By performing the substitution of the before-mentioned variables into the weak form of the problem, see Equation (20), we arrive to a problem entirely similar to that in Equations (33) and (34), in which the tangent stiffness matrix takes the form

$$ \mathbf{K}_t = \int_{\Omega_0} (\mathbf{B}^\mathrm{T} \mathbf{D} \mathbf{B} + \mathbf{G}^\mathrm{T} \tilde{\mathbf{S}}_0 \mathbf{G}) \, d\Omega \quad (49) $$
where

$$ \mathbf{D} = \lambda \left( \frac{1}{J_0^2} \mathbf{C}_0^{-1} \mathbf{M}_0^\mathrm{T} \right) + 2(\mu - \lambda \ln J_0) \left( \frac{1}{J_0^2} (\mathbf{C}_0^{-1} \mathbf{M}_0^\mathrm{T}) - \tilde{\mathbf{C}}_0 \right) \quad (50) $$

now takes into account the material non-linearity and has a somewhat similar appearance to the Lagrangian elastic tensor at the initial state. $J_0$ and $C_0$ represent the Jacobian and right Cauchy-Green strain tensor of the initial solution. $M_0$ is obtained from the series expansion of the Jacobian, and contains minors of $C_0$. Finally, $\tilde{C}_0$ is obtained from the series expansion of $C^{-1}$ and contains components of $C_0$, arranged in a particular way.

The geometrical non-linearities are included in the matrices **B**, **G** and $\tilde{\mathbf{S}}_0$. **B** represents the usual strain-displacement matrix, **G** relates the nodal displacements **u** and the gradient of displacements vector, and, finally, $\tilde{\mathbf{S}}_0$ represents a matrix that contains the initial stresses (we have chosen the same notation as in [18]).

In the right hand side of Equation (34), the non-linear load vector $\mathbf{f}_p^{nl}$ is a vector containing information of material and geometrical non-linearities of all order problems ranging from order one to $p-1$. It can be written as:

$$ \mathbf{f}_p^{nl} = \int_{\Omega_0} (\mathbf{B}^{\mathrm{T}}(\mathbf{S}_p^{\mathrm{nlmat}} + \mathbf{S}_p^{\mathrm{nlgeom}}) + \mathbf{G}^{\mathrm{T}}\mathbf{S}_p^{*}) \, d\Omega \quad (51) $$

As in the stiffness matrix, $\mathbf{S}_p^{\mathrm{nlgeom}}$ and $\mathbf{S}_p^*$ represent the standard matrices found in the literature when ANM is used to solve geometrical non-linear problems with linear materials. $\mathbf{S}_p^{\mathrm{nlmat}}$ takes into account the material behaviour:

$$ \mathbf{S}_p^{\mathrm{nlmat}} = (\lambda \ln J_0 - \mu) \left( \mathbf{C}\mathbf{C}_0 \left( RZ_p - \frac{RJ_p}{J_0^4} \right) + \frac{\mathbf{R}\mathbf{C}\mathbf{C}_p}{J_0^2} + \mathbf{R}\mathbf{C}_p^{-1} \right) \\ + \lambda \left( \frac{\mathbf{C}\mathbf{C}_0}{J_0^2} \left( RY_p + \frac{RJ_p}{2J_0^2} \right) + \mathbf{R}\mathbf{S}_p \right) \quad (52) $$

In this equation, $\mathbf{C}\mathbf{C}_0$ represents the cofactor matrix of $C_0$ and $\mathbf{R}\mathbf{C}\mathbf{C}_p$ is a vector containing values of $C_{ij}$ of all problems from order one to $p-1$, obtained when the cofactor matrix of $C$ is expanded in Taylor series:

$$ \mathbf{C}\mathbf{C}_p = \tilde{\mathbf{C}}_0\mathbf{C}_p + \mathbf{R}\mathbf{C}\mathbf{C}_p \quad (53) $$

$$ \mathbf{R}\mathbf{C}\mathbf{C}_p = \sum_{r=1}^{p-1} \begin{pmatrix} C_{22}^r C_{33}^{p-r} - C_{23}^r C_{23}^{p-r} \\ C_{11}^r C_{33}^{p-r} - C_{13}^r C_{13}^{p-r} \\ C_{11}^r C_{22}^{p-r} - C_{12}^r C_{12}^{p-r} \\ C_{13}^r C_{23}^{p-r} - C_{12}^r C_{33}^{p-r} \\ C_{13}^r C_{12}^{p-r} - C_{11}^r C_{23}^{p-r} \\ C_{12}^r C_{23}^{p-r} - C_{13}^r C_{22}^{p-r} \end{pmatrix} \quad (54) $$
Here, $RJ_p$ is a summation of products of different components of $C_p$ and is obtained when the squared Jacobian is expanded in Taylor series:

$$ (J^2)_p = M_0^T C_p + R J_p \quad (55) $$

$RS_p$ collects the terms concerning the expansion of $Y=\ln J$ and $C^{-1}$:

$$ \mathbf{RS}_p = \sum_{r=1}^{p-1} Y_r C_{p-r}^{-1} \quad (56) $$

$RC_p^{-1}$ collects the terms concerning $Z=J^{-2}$ and cofactor matrix of C expansions:

$$ \mathbf{RC}_p^{-1} = \sum_{r=1}^{p-1} Z_r \mathbf{CC}_{p-r} \quad (57) $$

Finally, it is necessary to expand $Y=\ln J$ and $Z=J^{-2}$ by using the Taylor series and the chain rule generalized to higher derivatives:

$$ Y_p = \frac{1}{2J_0^2} (J^2)_p + RY_p \quad \text{and} \quad Z_p = \frac{-1}{J_0^4} (J^2)_p + RZ_p \quad (58) $$

where

$$ RY_1 = 0 $$

$$ RY_2 = \frac{-1}{4J_0^4} (J^2)_1^2 $$

$$ RY_3 = \frac{1}{6J_0^6} (J^2)_1^3 + 2 \frac{-1}{4J_0^4} (J^2)_1 (J^2)_2 $$

$$ RZ_1 = 0 $$

$$ RZ_2 = \frac{1}{J_0^6} (J^2)_1^2 $$

$$ RZ_3 = \frac{-1}{J_0^8} (J^2)_1^3 + 2 \frac{1}{J_0^6} (J^2)_1 (J^2)_2 $$

$\vdots$

Once the expansion of the ANM has been performed, the next step of the method consists of the use in this framework of a reduced model, as explained in the previous sections for a general case, regardless of the constitutive model chosen.
4. NUMERICAL EXAMPLES

In order to show the performance of this method we have applied it to some hyperelastic, academic examples.

4.1. *Beam in traction*

In this example we have applied four concentrated 100 N forces to the left face of a beam in the axial direction. The nodes on the right face are constrained in the axial direction. The length of the beam is 400 mm and the cross-section is a square with sides of length 40 mm. The beam is shown in Figure 1. The material properties of the beam are E = 1 MPa and ν = 0.25.

4.1.1. *Kirchhoff–Saint Venant material.* In this example, the dimension of the basis of the reduced model is six, which is the number of modes necessary to capture 99.9% of the energy of the system. They are depicted in Figure 2. The load–displacement curve has been obtained for a node on the right face and is shown in Figure 3 for *p* = 1, . . . , 6. The linear solution obtained with *p* = 1 is the same as the one that can be obtained using standard POD only, without updating the tangent stiffness matrix. In this case, for the higher *p* the load–displacement curves are clearly non-linear and in this example they coincide with the analytic solution and the one calculated using the full FE Newton–Raphson equilibrium iterations (the curves are indistinguishable for *p*≥2).

Notice how, for *p*≥2, the radius of convergence of the proposed technique is very high (on the order of 60 mm of tip displacement) for a beam with a total initial length of 400 mm. If we need to follow the equilibrium path beyond this point, the method can be restarted at any equilibrium point, by a new computation of the stiffness matrix at that point.

<figure><img src="image_1.png" /><figcaption>Figure 1. Geometry of the beam in traction.</figcaption></figure>
<figure><img src="image_2.png" /><figcaption>Figure 2. Six most important eigenmodes for the simulation of the Kirchhoff–Saint Venant beam under traction. The corresponding eigenvalues are: (a) 38 488.48; (b) 0.04; (c) 1.04E−11; (d) 3.44E−12; (e) 1.34E−12; and (f) 1.28E−12.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 3. Beam in traction. Load–displacement results (in terms of $\bar{\lambda}$) for different orders of approximation, $p$, compared with the analytical solution. Kirchhoff–Saint Venant behaviour.</figcaption></figure>

4.1.2. Neo-Hookean material. For the neo-Hookean material the behaviour of the proposed technique is very similar to that of the Kirchhoff–Saint Venant. The load–displacement curve is depicted in Figure 4, where an excellent agreement between the FEM results and the $p=4$ reduced model has been noticed.

In Figure 5 the six most important modes of the complete solution, obtained by Newton–Raphson methods, are depicted.
<figure><img src="image_4.png" /><figcaption>Figure 4. Beam in traction. Load-displacement results (in terms of $\bar{\lambda}$) for different orders of approximation, $p$, compared with the analytical solution. Neo-Hookean behaviour.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 5. Six most important eigenmodes for the simulation of the neo-Hookean beam under traction. The corresponding eigenvalues are: (a) 38.488; (b) 0.04; (c) 1.04E-11; (d) 3.44E-12; (e) 1.34E-12; and (f) 1.28E-12.</figcaption></figure>

The excellent agreement of the results for very large beam tip displacements also deserve some comments. Note that we have plotted the solution of the FE, complete model until a tip displacement of around 120 mm. Until this very large level of strain (the beam is 400 mm long), the reported solution by the reduced model with an order $p=4$ is remarkable, while it seems to reproduce the 'expected' solution far beyond.
## 4.2. Pinched hemisphere

In this example, a concentrated force of 30 N is applied to the pole of a hemisphere and directed towards its centre. The bottom of the hemisphere is fixed. Its radius is 25 mm and its Young's modulus and the Poisson ratio are 2 MPa and 0.48, respectively. The mesh is shown in Figure 6 and it is composed of 614 nodes and 448 linear hexagonal elements. Kirchhoff-Saint Venant and neo-Hookean behaviours are assumed in this example.

### 4.2.1. Kirchhoff-saint venant material.
In this example the empirical basis has 10 modes, depicted in Figure 8. The load-displacement curves for the node on the pole of the hemisphere for $p=1, \dots, 6$ have been obtained and depicted in Figure 7.

<figure><img src="image_6.png" /><figcaption>Figure 6. Hemisphere mesh.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 7. Load-displacement curves for the hemisphere problem.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 8. Ten most important eigenmodes employed for the simulation of the pinched hemisphere. The corresponding eigenvalues are: (a) 25.42; (b) 0.19; (c) 0.02; (d) 4.80E-05; (e) 6.04E-07; (f) 1.15E-08; (g) 2.85E-10; (h) 4.12E-12; (i) 6.52E-14; and (j) 2.32E-16. Kirchhoff-Saint Venant behaviour.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 9. Deformed hemisphere obtained using the combined technique POD-ANM.</figcaption></figure>

Again it is worth noting that the solution with $p=1$ that is linear is the same as the solution that one obtains using standard model reduction without stiffness matrix updating. For $p \le 5$ it can be seen that the ANM-POD solution has a good agreement with the solution computed using full Newton-Raphson equilibrium iteration up to $u \approx -2.5$ mm, which can be considered as the radius of convergence (note that this is in the order of 10% of the sphere radius). But for $p=6$ the size of the radius of convergence is remarkable, going beyond 4 mm (more than 20% of the radius of the sphere). Again, if the load is to be extended beyond this point, a new tangent stiffness matrix should be computed at this point, performing again the algorithm given by Equations (42)–(43).

The deformed hemisphere for $p=6$ and $\lambda=1$ is shown in Figure 9 and the displacement field in the $z$ direction, $u_z$, obtained by Newton-Raphson procedures is shown in Figure 10.

4.2.2. *Neo-Hookean behaviour.* As in the previous example, the results were tested also for the neo-Hookean behaviour, showing similar levels of accuracy as the Kirchhoff-Saint Venant examples. The modes employed as global basis in the example are depicted in Figure 11.

Results for the complete model, solved by means of the full Newton-Raphson iterations, and the reduced model, for different approximation orders, are compared in Figure 12. The great similitude for both the models is noteworthy. Again, the load-displacement (in terms of the load factor, $\bar{\lambda}$) curve shows great similarity for $p=2, \dots, 4$ and a wide convergence radius.

4.3. *Pinching the human cornea*

As mentioned before, real-time simulation in surgical environments is one of the fields in which very fast simulations are needed. This field arises as a natural potential application of the technique
<figure><img src="image_10.png" /><figcaption>Figure 10. Deformed hemisphere computed by the full FE Newton-Raphson equilibrium iteration.</figcaption></figure>

proposed here, which can solve with moderate accuracy simplified models of non-linear solids
without the need for tangent stiffness matrix updating nor iterative procedures.

In this example we consider a human cornea, whose geometry has been taken from [21]. Forces of 0.014N each have been applied to nine neighbour nodes located at the centre of the cornea, thus simulating the contact of a rounded tool. The cornea was meshed using trilinear hexahedral elements. It consisted of 8514 nodes and 7182 elements. The mesh is shown in Figure 13 in two views. The cornea is clamped at its base, giving a dome-like problem that showed buckling under some types of loads.

4.3.1. *Kirchhoff-Saint Venant model.* The material properties of the cornea are assumed to be $E=2$ MPa and $v=0.48$, although a realistic model of corneal tissue is far more complicated [13].

In this example nine modes were applied that provide decent approximation. The solution has been obtained using ANM-POD for $p=1, \dots, 6$, plotted in Figure 14. Note how the first mode resembles very closely the final solution. Higher modes concentrate in capturing the solution near the region where the load is applied.

In order to verify the results we have computed the solution by the full FE Newton-Raphson method. The loading factor ($\bar{\lambda}$) has been plotted versus the minimum displacement in Figure 15.

As it can be noticed the results have good accuracy with Newton-Raphson solution. The deformed cornea obtained using ANM-POD for $\bar{\lambda}=1$ is shown in Figure 16 and the one obtained using the full FE Newton-Raphson is depicted in Figure 17.

Note that, despite the high Poisson's ratio employed, close to 0.5, the method shows no tendency to lock in the examples studied (although it is well known that this does not constitute a valid
<figure><img src="image_11.png" /><figcaption>Figure 11. Six most important eigenmodes for the simulation of the neo-Hookean pinched hemisphere. The corresponding eigenvalues are: (a) 30.0; (b) 0.013; (c) 2.95E-05; (d) 9.42E-08; (e) 3.12E-10; (f) and 4.14E-12.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 12. Load–displacement curves for the hemisphere problem with neo-Hookean behaviour.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 13. Geometry of the finite element model for the human cornea.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 14. Six most important eigenmodes for the simulation of the cornea. The corresponding eigenvalues are: (a) 67.5; (b) 0.07; (c) 5.95−4; (d) 7.54E−7; (e) 1.26E−8; and (f) 3.48E−10.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 15. The loading factor vs minimum displacement for the Pinched cornea. Kirchhoff–Saint Venant model.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 16. $u_y$-contour of the pinched cornea obtained by ANM-POD.</figcaption></figure>

proof of robustness). The use of truly incompressible (mixed) formulations is currently under investigation.

4.3.2. *Neo-Hookean model.* Again, results for the neo-Hookean cornea follow the same guidelines of accuracy as those for the Kirchhoff–Saint Venant model. With an expansion of only three or four terms (see Figure 18) a great accuracy is obtained along a big radius of convergence. In this case, only six modes were enough to simulate the cornea, thus making a stiffness matrix of size 6 × 6. These modes are depicted in Figure 19.
<figure><img src="image_17.png" /><figcaption>Figure 17. $u_y$-contour of the pinched cornea obtained by a full FE model with Newton-Raphson iterations.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 18. The loading factor vs minimum displacement for the Pinched cornea. Neo-Hookean model.</figcaption></figure>

### 4.4. *Timing*

An essential feature of the proposed technique is the reduction of computing time that it allows. In order to have a final, overall, impression of the savings achieved with this technique, we reproduce here the total amount of time (in terms of CPU seconds) employed for the three examples, in the neo-Hookean version (timings for the Kirchhoff-Saint Venant problems showed similar values). All the examples were performed on a PC equipped with two processors (only one was employed,
<figure><img src="image_19.png" /><figcaption>Figure 19. Six most important eigenmodes for the simulation of the neo-Hookean cornea. The corresponding eigenvalues are: (a) 53.5; (b) 2.44E-02; (c) 2.09E-04; (d) 7.60E-08; (e) 3.78E-10; and (f) 7.19E-13.</figcaption></figure>

<table><caption>Table I. Time (CPU seconds) employed for the solution of the three examples.</caption><thead><tr><th>Example</th><th>FEM+Newton-Raphson</th><th>POD+ANM</th></tr></thead><tbody><tr><td>Beam</td><td>1.48</td><td>0.11</td></tr><tr><td>Hemisphere</td><td>225.1</td><td>3.92</td></tr><tr><td>Cornea</td><td>8.64E04</td><td>79.62</td></tr></tbody></table>

no parallel computing was used) AMD Quad Opteron running at 2.2 GHz and with 16 Gb RAM, under Scientific Linux and running Matlab v. 2007.

As can be noticed, savings are more impressive the larger the model is. Notice also that in these examples no tangent stiffness matrix updating was performed, which eventually leads to some loss of accuracy for very large strains. But even if such an update would have been accomplished, the presented technique can be considered as very competitive for some applications (Table I).

## 5. CONCLUSIONS

In this paper, we study a technique for the construction of reduced models of hyperelastic solids. It is based on the use of model reduction techniques (based upon proper orthogonal decomposition
or Karhunen-Loève methods) and an asymptotic expansion of the solution in the neighbourhood
of the last converged equilibrium point. This technique was previously introduced in [16] for linear
(finite strain) elastic solids and has been generalized here for hyperelastic materials involving
material non-linearities.

These two ingredients render a method with some very attractive features, such as the use of
very few degrees of freedom (up to 10 in the examples showed in this paper) and the need of
only one tangent stiffness matrix within the region covered by the radius of convergence of the
asymptotic expansion. This radius of convergence is very often relatively high (on the order of
10% of the total dimensions of the solid for the examples showed herein). Within this radius of
convergence no iterative procedures (typical of modified Newton methods) nor the update of the
stiffness matrix (typical in standard Newton-Raphson methods) are needed.

The examples developed in this work showed that 6–10 eigenmodes are enough for the vast
majority of the examples considered, and usually an expansion of order 4–6 renders very acceptable
results, with excellent accuracy within a vast radius of convergence.

The final simplicity of the resulting method suggests its usage in real-time environments, for
instance, where feedback frequencies of 30 Hz for realistic video streaming are needed. Even
haptic simulations could be envisaged as a potential application of this kind of techniques, where
500–1000 Hz are needed for realistic feedback of forces. The evaluation of the true performance
of these algorithms in the true prototype codes is currently under investigation by the authors.

REFERENCES

1. Idelsohn SR, Cardona R. A reduction method for nonlinear structural dynamics analysis. *Computer Methods in Applied Mechanics and Engineering* 1985; **49**:253-279.

2. Krysl P, Lall S, Marsden JE. Dimensional model reduction in non-linear finite element dynamics of solids and structures. *International Journal for Numerical Methods in Engineering* 2001; **51**:479–504.

3. Jernej B, Doug James L. Real-time subspace integration for St. Venant–Kirchhoff deformable models. *ACM Transactions on Graphics (SIGGRAPH 2005)*, Los Angeles, U.S.A., vol. 24(3), August 2005; 982–990.

4. Park HM, Cho DH. The use of the Karhunen-Loève decomposition for the modeling of distributed parameter systems. *Chemical Engineering Science* 1996; **51**(1):81-98.

5. Sirovich L. Turbulence and the dynamics of coherent structures coherent structures part I: coherent structures. *Quaterly of Applied Mathematics* 1987; *XLV*:561–571.

6. Lorenz EN. *Empirical Orthogonal Functions and Statistical Weather Prediction*. MIT, Department of Meteorology, Scientific Report Number 1, Statistical Forecasting Project, 1956.

7. Karhunen K. Uber lineare methoden in der wahrscheinlichkeitsrechnung. *Annales Academiae Scientiarum Fennicae*, Series Al. Mathematical Physics 1946; **37**:1–79.

8. Loève MM. *Probability Theory* (3rd edn). The University Series in Higher Mathematics. Van Nostrand: Princeton, NJ, 1963.

9. Ryckelynck D, Chinesta F, Cueto E, Ammar A. On the a priori model reduction: overview and recent developments. *Archives of Computational Methods in Engineering* 2006; **12**(1):91–128.

10. Ryckelynck D. A priori model reduction method for the optimization of complex problems. *Workshop on Optimal Design of Materials and Structures*, Ecole Polytechnique, Palaiseau, Paris (France), 2003.

11. Ryckelynck D. A priori hyperreduction method: an adaptive approach. *Journal of Computational Physics* 2005; **202**(1):346–366.

12. Tenenbaum JB, de Silva V, Langford JC. A global framework for nonlinear dimensionality reduction. *Science* 2000; **290**:2319–2323.

13. Niroomandi S, Alfaro I, Cueto E, Chinesta F. Real-time deformable models of non-linear tissues by model reduction techniques. *Computer Methods and Programs in Biomedicine* 2008; **91**(3):223–231.

14. Cochelin B, Damil N, Potier-Ferry M. Asymptotic-numerical methods and Padé approximants for non-linear elastic structures. *International Journal for Numerical Methods in Engineering* 1994; **37**:1187–1213.
15. Abichou H, Zahrouni H, Potier-Ferry M. Asymptotic numerical method for problems coupling several nonlinearities. *Computer Methods in Applied Mechanics and Engineering* 2002; **191**(51–52):5795–5810.

16. Yvonnet J, Zahrouni H, Potier-Ferry M. A model reduction method for the post-buckling analysis of cellular microstructures. *Computer Methods in Applied Mechanics and Engineering* 2007; **197**:265–280.

17. Ryckelynck D. Hyper-reduction of mechanical models involving internal variables. *International Journal for Numerical Methods in Engineering* 2008; **77**(1):75–89.

18. Cochelin B, Damil N, Potier-Ferry M. The asymptotic numerical method: an efficient perturbation technique for nonlinear structural mechanics. *Revue Europeenne des Elements Finis* 1994; **3**:281–297.

19. Bonet J, Wood RD. *Nonlinear Continuum Mechanics for Finite Element Analysis*. Cambridge University Press: Cambridge, 2008.

20. Cao H-L, Potier-Ferry M. An improved iterative method for large strain viscoplastic problems. *International Journal for Numerical Methods in Engineering* 1999; **44**:155–176.

21. Alastrué V, Calvo B, Peña E, Doblaré M. Biomechanical modeling of refractive corneal surgery. *Journal of Biomechanical Engineering—Transactions of the ASME* 2006; **128**:150–160.