# Model Order Reduction in Hyperelasticity: A Proper Generalized Decomposition Approach

Siamak Niroomandi, Iciar Alfaro, David Gonzalez, Elías Cueto, Francisco Chinesta

► To cite this version:

Siamak Niroomandi, Iciar Alfaro, David Gonzalez, Elías Cueto, Francisco Chinesta. Model Order Reduction in Hyperelasticity: A Proper Generalized Decomposition Approach. International Journal for Numerical Methods in Engineering, 2013, 96 (3), pp.129-149. 10.1002/nme.4531. hal-01007060

HAL Id: hal-01007060

https://hal.science/hal-01007060

Submitted on 11 Dec 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
Model order reduction in hyperelasticity: a proper generalized
decomposition approach

Siamak Niroomandi¹, Icíar Alfaro¹, David González¹,
Elías Cueto¹,*,† and Francisco Chinesta²

¹Aragón Institute of Engineering Research (13A), Universidad de Zaragoza, Edificio Betancourt,
María de Luna, s.n. 50018 Zaragoza, Spain

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

This paper deals with the extension of proper generalized decomposition methods to non-linear problems,
in particular, to hyperelasticity. Among the different approaches that can be considered for the linearization
of the doubly weak form of the problem, we have implemented a new one that uses asymptotic numerical
methods in conjunction with proper generalized decomposition to avoid complex consistent linearization
schemes necessary in Newton strategies. This approach results in an approximation of the problem solution
in the form of a series expansion. Each term of the series is expressed as a finite sum of separated functions.
The advantage of this approach is the presence of only one tangent operator, identical for every term in the
series. The resulting approach has proved to render very accurate results that can be stored in the form of a
meta-model in a very compact format. This opens the possibility to use these results in real-time, reaching
kHz feedback rates, or to be used in deployed, handheld devices such as smartphones and tablets.

KEY WORDS: proper generalized decomposition; separated representations; asymptotic numerical methods; model order reduction; real time simulation; hyperelasticity

# 1. INTRODUCTION

Proper orthogonal decomposition (POD) methods have had a tremendous importance in many
branches of applied sciences and engineering, where they are known under a variety of names.
Starting from the pioneer works by Karhunen [1] and Loève [2] (also Lorenz [3]), POD methods,
also known as principal component analysis and subspace tracking, among other names, have been
applied to a wide range of problems, from chemical engineering [4] to real time visualization,
possibly with haptic response [5], structural mechanics [6, 7], or virtual surgery [8–10].

Proper orthogonal decomposition rely upon the computation of the solution to surrogate prob-
lems, the so-called snapshots. The most typical structure of these snapshots, obtained as an eigen-
value problem through a minimization process, is then used as a basis for subsequent problem
solving. Among the problems related to this approach, one can cite the need for an important num-
ber of snapshots. These snapshots should be generated by solving problems somewhat similar to
the one at hand, but maybe for slight changes in geometry or boundary conditions. The so-called a
priori POD [11–13] constitutes one of the first attempts to avoid the costly generation of snapshots
by starting from scratch and adding Krylov subspaces to the basis of the reduced model.
Another difficulty on the use of POD-based model order reduction is related to the need of recom-
puting the original system of equations in non-linear problems (even if it will not be solved) to
update tangent stiffness matrices. This prevents the method to provide as competitive costs as one
would expect from the limited size of the resulting systems of equations. This problem has gener-
ated recently quite an active area of research. As notable contributions, one can cite the empirical
interpolation method [14] or its discrete counterpart [15]. A somewhat related approach can be
found in [16].

On the other hand, non-linear POD methods (those designed for data lying on a manifold) are
always difficult to develop. See [5, 17–19] to name a few. Interpolation among reduced models
continues to be an issue too [18, 20].

Proper generalized decomposition (PGD) methods, on the contrary, arose recently as a gener-
alization of POD techniques. PGD roots can be traced back to the pioneer work by P. Ladeveze
on the LATIN method [21] and, particularly, the so-called radial loading approximation scheme,
in which a separated space–time approximation of the displacement in structural mechanics prob-
lems was used. Independently, the method was re-invented in the framework of problems defined
in high dimensional state spaces [22, 23]. It was then soon realized that PGD methods can be con-
sidered as a generalization of POD, in which the basis are computed on the fly without no previous
snapshot. Instead, the essential field is approximated as a finite sum of separable functions, very
much similar to the space–time structure of the radial approximation within the LATIN method. For
recent surveys on PGD methods, the reader is referred to [24–26]. Some of the nice properties of
POD are lost, however. For instance, optimality of the POD basis is no more guaranteed, although
convergence properties have been demonstrated recently [27], and error estimators have also been
proposed [28, 29].

When compared with the existing model order reduction techniques, such as the ones based on
POD, for instance, PGD provides the modes a priori, that is, without the need to compute snapshots
of the complete model. These snapshots would be extremely costly to generate for a problem defined
in high dimensional phase spaces, where the sampling of the space is not obvious. With respect to
the technique introduced in [5], on the contrary, it is stated that ‘deformation basis generation is a
hard open problem in solid mechanics, and there exist no algorithms for automatic proven-quality
global deformation basis generation under general forcing’. This is precisely where PGD comes into
play. The method here proposed provides with a suitable set of reduced basis a priori, without the
need to compute costly snapshots. In [5], two different methods were chosen to generate the basis;
one based upon standard modal analysis that includes derivatives of the eigenmodes. The second
one is basically a POD procedure in which the modes are extracted from a linear simulation and
successively enriched off-line. However, although in this reference, it is stated that the method is
general and can be applied to constitutive models other than Kirchhoff–Saint Venant (which is of
low practical interest because of its well-known instabilities in compression due to lack of convex-
ity of the strain energy functional), we believe that it is difficult to generalize for state-of-the-art
constitutive hyperelastic models. Up to our knowledge, it has not been generalized yet.

One key aspect in the rapid development of PGD methods is related to the fact that PGD can be
seen as both an efficient technique for high dimensional problems and as an a priori model order
reduction technique. This opened the door to re-interpreting parametric problems as high dimen-
sional ones, just by considering parameters as new dimensions of the state space of the model
[30–32].

This new point of view offered new and interesting insights in many branches of applied science
and engineering. It is now envisaged how by looking at a parametric problem from a multidici-
mensional framework gives rise to important savings in simulation time. Furthermore, it has been
observed how important advantages could be obtained by formulating standard problems as high
dimensional ones. For instance, formulating plate and shell problems as 3D ones but with 2D com-
putational complexity [33] has revealed very interesting insights on a general finite element (FE)
formulation beyond Kirchhoff and Reissner–Mindlin plates.

In this work, we explore the possibilities of PGD methods applied to fast (real time) simula-
tion of hyperelastic solids. For a related reference, the interested reader can also consult [34]. It
will be shown how a multidimensional formulation of the problem, in which the displacement is
expressed as to be dependent on both the physical coordinates and the position and orientation of the applied loads (thus, defined in $\mathbb{R}^9$), opens the door to simulations with feedback rates on the order of 1 kHz. The developed strategy comprises an off-line computation strategy, in which a high-dimensional problem is solved. This solution provides in fact a sort of metamodel that can be stored in a very compact form. Then, an online simulation strategy is developed that solves the metamodel at impressive feedback rates, even on handheld devices.

But the price to pay in this case is the complex non-linear structure of the weak form of the problem, as will be shown in Section 2. Consistent linearization of PGD weak forms continues to be an issue in our community. Previous works include explicit linearizations [35]. In this work, we introduce a new approach by combining a PGD formulation of the problem with an asymptotic expansion of the displacement. The use of asymptotic expansions in the search of complex non-linear equilibrium paths of structures is due to Potier-Ferry and coworkers [36–39]. This gives rise to a formulation in which only one tangent operator is present, the same for all the terms of the expansion, and where the non-linearities are translated to the right-hand-side (force) term of the matrix form of the problem. The resulting technique provides results, as will be shown, at kHz feedback rates for general hyperelastic models. This constitutes, in our opinion, a major advancement in the state of the art of reduced order models. Up to our knowledge, only techniques working with Kirchhoff–Saint Venant model were able to obtain such performance rates.

In Section 3, we review the basics of asymptotic numerical methods (ANM) in their standard form, whereas in Section 4, we develop the combined PGD–ANM formalism and study how it provides a very convenient form for non-linear path following of hyperelastic solids and structures. In Section 5, we study the performance of the proposed technique by employing some classical benchmark tests and show how the resulting formalism can be employed advantageously to perform real-time simulation.

## 2. FORMULATION OF THE PROBLEM UNDER THE PGD POINT OF VIEW

As mentioned before, one of the main advantages of using PGD for model order reduction purposes relies in the possibility of rewriting the model in a multidimensional framework. Consider, for simplicity, the static equilibrium equations of a general solid under small strain assumptions,

$$ \nabla \cdot \sigma + b = 0 \text{ in } \Omega, \quad (1) $$

where $b$ represents the volumetric forces applied to the body, subjected to the following boundary conditions:

$$ u = \bar{u} \text{ on } \Gamma_u \quad (2) $$

$$ \sigma n = \bar{t} \text{ on } \Gamma_t. \quad (3) $$

The standard weak form of the problem is obtained after multiplying both sides of Equation (1) by an admissible variation of the displacement, $u^*$, and integrating over the domain $\Omega$. To fully exploit the characteristics of PGD methods, following [32], we convert the equilibrium problem given by Equations (1)–(3) into a parametric one, by considering, again for the sake of simplicity in the exposition, that the load $\bar{t}$ is punctual, unitary, vertical, and is applied at a position $s \in \Gamma_t$ that acts as a parameter in the formulation. Note that this alternative formulation of the problem gives rise to a new one defined in $\mathbb{R}^6$, because $u = u(x, s) \in \Omega \times \bar{\Gamma}$, where $\bar{\Gamma} \subseteq \Gamma_t$ represents the portion of the boundary where the load can be applied.

If inertia terms in Equation (1) were non-negligible, a possible approach that considers initial conditions as additional parameters has also been explored in [40].

In this spirit, an alternative (doubly) weak form of problem (1)–(3) consists in finding the displacement $u \in \mathcal{H}^1(\Omega) \times L_2(\bar{\Gamma})$ such that for all $u^* \in \mathcal{H}_0^1(\Omega) \times L_2(\bar{\Gamma})$ Siamak5,

$$ \int_{\bar{\Gamma}} \int_{\Omega} (\nabla_s u^*)^T \sigma d\Omega d\bar{\Gamma} = \int_{\bar{\Gamma}} \int_{\Gamma_{t2}} (u^*)^T t d\Gamma d\bar{\Gamma}, \quad (4) $$
where $\Gamma = \Gamma_u \cup \Gamma_t$ represents the boundary of the solid, divided into essential and natural regions, and where $\Gamma_t = \Gamma_{t1} \cup \Gamma_{t2}$, that is, regions of homogeneous and non-homogeneous, respectively, natural boundary conditions. In turn, we assume $t = e_z \delta(x - s)$, where $\delta$ represents the Dirac delta function and $e_z$ represents the unit vector along the $z$ direction, in this case. Note that a general-form load term does not include further complexity into this formulation. This Dirac delta term should be regularized for computation purposes and approximated by

$$t_j \approx \sum_{i=1}^{n} f_j^i(x) g_j^i(s) \qquad (5)$$

by simply performing a singular value decomposition of the load, for instance.

As explained before, PGD methods assume a separated representation of the unknown field (here, the displacement). This is the key ingredient of the method that allows solving efficiently high dimensional problems. The so-called *curse of dimensionality* associated to mesh-based solution of high dimensional problems is thus avoided by solving a sequence of low-dimensional problems in separated form. This was in fact the key ingredient in the space-time separated representation of the displacement in the LATIN method [21]. The PGD approach to the problem is characterized by the construction, in an iterative way, of an approximation to the solution in the form of a finite sum of separable functions. Assume that we have converged to a solution, at iteration $n$, of this procedure,

$$u_j^n(x, s) = \sum_{k=1}^{n} X_j^k(x) \cdot Y_j^k(s), \qquad (6)$$

where the term $u_j$ refers to the $j$th component of the displacement vector, $j = 1, 2, 3$.

The following term of this approximation, the $(n+1)$th one, will be similar to

$$u_j^{n+1}(\mathbf{x}, s) = u_j^n(\mathbf{x}, s) + R_j(\mathbf{x}) \cdot S_j(s), \qquad (7)$$

where $R(\mathbf{x})$ and $S(s)$ are the sought functions that improve the approximation. In this same way, the admissible variation of the displacement will be given by

$$u_j^*(\mathbf{x}, s) = R_j^*(\mathbf{x}) \cdot S_j(s) + R_j(\mathbf{x}) \cdot S_j^*(s). \qquad (8)$$

At this point, several options are at hand so as to determine the new pair of functions $R_j$ and $S_j$. The most frequently used, due to both its ease of implementation and good convergence properties, in general, is a fixed-point alternating directions algorithm in which functions $R_j$ and $S_j$ are sought iteratively. We describe briefly the implementation of this algorithm.

## 2.1. Computation of $S(s)$ assuming $R(\mathbf{x})$ is known

In this case, following standard assumptions of variational calculus, we have

$$u_j^*(\mathbf{x}, s) = R_j(\mathbf{x}) \cdot S_j^*(s), \qquad (9)$$

or, equivalently, $u^*(\mathbf{x}, s) = R \circ S^*$, where the symbol ‘$\circ$’ denotes the so-called entry-wise Hadamard or Schur multiplication for vectors. Once substituted into Equation (4), gives

$$\begin{aligned} \int_{\bar{\Gamma}} \int_{\Omega} \nabla_s (R \circ S^*) : \mathbf{C} : \nabla_s \left( \sum_{k=1}^{n} X^k \circ Y^k + R \circ S \right) d\Omega d\bar{\Gamma} \\ = \int_{\bar{\Gamma}} \int_{\Gamma_{t2}} (R \circ S^*) \cdot \sum_{k=1}^{m} f^k \circ g^k \right) d\Gamma d\bar{\Gamma}, \end{aligned} \qquad (10)$$

or equivalently (we omit obvious functional dependencies)

$$\begin{aligned} \int_{\bar{\Gamma}} \int_{\Omega} \nabla_s (R \circ S^*) : \mathbf{C} : \nabla_s (R \circ S) d\Omega d\bar{\Gamma} \\ = \int_{\bar{\Gamma}} \int_{\Gamma_{t2}} (R \circ S^*) \cdot \sum_{k=1}^{m} f^k \circ g^k \right) d\Gamma d\bar{\Gamma} - \int_{\bar{\Gamma}} \int_{\Omega} \nabla_s (R \circ S^*) \cdot \mathcal{R}^n d\Omega d\bar{\Gamma}, \end{aligned} \qquad (11)$$
where $\mathcal{R}^n$ represents

$$ \mathcal{R}^n = \mathbf{C} : \nabla_s \mathbf{u}^n. \qquad (12) $$

All the terms depending on $\mathbf{x}$ are known, and hence we can compute all integrals over $\Omega$ and $\Gamma_{t2}$ (support of the regularization of the initially punctual load) to derive an equation to compute $\mathbf{S}(\mathbf{s})$.

## 2.2. Computation of $\mathcal{R}(\mathbf{x})$ assuming $\mathbf{S}(\mathbf{s})$ is known

Equivalently, in this case, we have

$$ u_j^*(\mathbf{x}, \mathbf{s}) = R_j^*(\mathbf{x}) \cdot S_j(\mathbf{s}), \qquad (13) $$

which, once substituted into Equation (4), gives

$$ \begin{aligned} & \int_{\bar{\Gamma}} \int_{\Omega} \nabla_s (\mathbf{R}^* \circ \mathbf{S}) : \mathbf{C} : \nabla_s \sum_{k=1}^{n} \mathbf{X}^k \circ \mathbf{Y}^k + \mathbf{R} \circ \mathbf{S} ) \, d\Omega d\bar{\Gamma} \\ & = \int_{\bar{\Gamma}} \int_{\Gamma_{t2}} (\mathbf{R}^* \circ \mathbf{S}) \cdot \sum_{k=1}^{m} \mathbf{f}^k \circ \mathbf{g}^k ) \, d\Gamma d\bar{\Gamma}. \end{aligned} \qquad (14) $$

In this case, all the terms depending on $\mathbf{s}$ (load position) can be integrated over $\bar{\Gamma}$, leading to a generalized elastic problem to compute function $\mathbf{R}(\mathbf{x})$.

This simple algorithm renders generally excellent convergence properties (see [25] and references therein).

While fixed-point approaches (see [25] and references therein) render generally excellent results, it is generally very difficult to perform consistent linearizations on Equations (11) and (14). In [35], some explicit algorithms have been tested to overcome this problem that generally work well. By this, we mean that they provide errors, for reasonable pseudo-time step sizes, of the same order of those obtained by existing POD-based approaches such as [17, 18]. This level of error is generally higher than that of common engineering practice but is nevertheless acceptable in the real-time simulation community [41].

In what follows, a method that combines asymptotic expansions of the displacement and the second Piola–Kirchhoff stress tensor to avoid complex stiffness matrix update procedures is reviewed. The method, due to Potier-Ferry and coworkers [36–39, 42], has been used in a variety of structural mechanics problems to follow complex equilibrium paths of solids and structures.

# 3. A BRIEF REVIEW OF THE ANM

For completeness and following closely [42], we review here the basics of the ANM applied to the simplest non-linear case, namely, Kirchhoff–Saint Venant solids (i.e., linear elastic solids undergoing large strains). Under a Lagrangian frame of reference, we consider the displacement as given by

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

Following the notation in [39], we consider a linear and a non-linear term for the Green-Lagrange strain tensor, $\mathbf{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 (16) $$

where $\mathbf{F} = \nabla \mathbf{u} + \mathbf{I}$ is the gradient of deformation tensor and

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

$$ \gamma_{nl}(\mathbf{u}, \mathbf{u}) = \frac{1}{2} \nabla (\mathbf{u}^T) \nabla(\mathbf{u}). \qquad (17b) $$
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

$$S = \frac{\partial \Psi}{\partial E}, \qquad (18)$$

that is a symmetric tensor and is related to the first Piola–Kirchhoff stress tensor, $P$, by $P = FS$. The equilibrium equation stated in the reference configuration is

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

in which $\mathbf{B}$ is the body force. The boundary conditions of the body are defined by (we do not consider the case of follower loads for simplicity)

$$\begin{aligned} \mathbf{u}(X) &= \bar{\mathbf{u}} \quad \text{on } \Gamma_u, \\ PN &= \bar{\lambda}\bar{\mathbf{t}} \quad \text{on } \Gamma_t, \end{aligned} \qquad (20)$$

where $N$ is the unit vector normal to $\Gamma = \partial\Omega_0$, $\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} S : E^* d\Omega = \bar{\lambda} \int_{\Gamma_t} \bar{\mathbf{t}} \cdot \mathbf{u}^* d\Gamma \quad \forall \mathbf{u}^* \in H^1(\Omega), \qquad (21)$$

where in the aforementioned equation, $E^*$ is expressed by

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

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

$$\gamma_{nl_S}(\mathbf{u}, \mathbf{u}^*) = \gamma_{nl}(\mathbf{u}, \mathbf{u}^*) + \gamma_{nl}(\mathbf{u}^*, \mathbf{u}). \qquad (23)$$

The Kirchhoff–Saint Venant model is characterized by the energy function given by

$$\Psi = \frac{\lambda}{2}(\text{tr}(E))^2 + \mu E : E, \qquad (24)$$

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

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

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

Although it is well known that the Kirchhoff–Saint Venant model is unstable under compression, and thus of limited importance in engineering applications, we have considered it here for the sake of simplicity in the formulation and because it is still in use in fields such as real-time simulation of surgery and rendering in general, where it is among the state of the art models [5, 17, 18].

Under these assumptions, the ANM is based upon expanding the displacement in the neighborhood of each material point in terms of a control parameter 'a'. This expansion is developed in the neighborhood of a known equilibrium point $(\mathbf{u}^m; \mathbf{S}^m; \bar{\lambda}^m)$ at step $m$, 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}^{m+1}(a) \\ \mathbf{S}^{m+1}(a) \\ \bar{\lambda}^{m+1}(a) \end{matrix} \right\} = \left\{ \begin{matrix} \mathbf{u}^m(a) \\ \mathbf{S}^m(a) \\ \bar{\lambda}^m(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 (26)$$

where $(\mathbf{u}_p, \mathbf{S}_p, \bar{\lambda}_p)$ are unknowns. Previously, $(\mathbf{u}^{m+1}(a), \mathbf{S}^{m+1}(a), \bar{\lambda}^{m+1}(a))$ represents the solution along a portion of the loading curve. Noteworthy, the behavior of the solid is described continuously with respect to 'a'. The introduction of Equation (26) into Equations (21) and (25) 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.

In general, any variable can be expanded in terms of $a$, so, for instance, the series expansion of $E^*(u)$ gives

$$E^{*m+1}(a) = \gamma_l(u^*) + \gamma_{nl_S}(u^*, u^m) + \sum_{p=1}^{N} a^p \gamma_{nl_S}(u^*, u_p). \quad (27)$$

As can be noticed, $E^{*m+1}$ includes terms in both $u$ and $u^*$. This is not surprising, because we are computing a correction to $u^m$. The series expansions of $S$ gives in turn

$$\begin{aligned} S^{m+1}(a) &= \mathbf{C} : E^{m+1}(a) \\ &= \mathbf{C} : \left[ \gamma_{nl}(u^m, u^m) + \gamma_l(u^m) + \sum_{p=1}^{N} a^p \gamma_l(u_p) + \gamma_{nl_S}(u^m, u_p) + \sum_{i=1}^{p-1} \gamma_{nl}(u_i, u_{p-i}) \right], \end{aligned} \quad (28)$$

and at order $p$, we obtain

$$S_p = \mathbf{C} : \left\{ \gamma_l(u_p) + \gamma_{nl_S}(u^m, u_p) + \sum_{i=1}^{p-1} \gamma_{nl}(u_i, u_{p-i}) \right\}. \quad (29)$$

Introducing the asymptotic expansion into Equation (25) results in

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

with $\Psi_{\text{ext}}(u^*) = \int_{\Gamma_t} t \cdot u^* d\Gamma$. Introducing Equation (29) into Equation (30) 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}(u^*, u^{m+1}) = \bar{\lambda}_p \Psi_{\text{ext}}(u^*) + F_p^{nl}(u^*) \quad (31)$$

with

$$\mathcal{L}(u^*, u^{m+1}) = \int_{\Omega} \{ S^m : \gamma_{nl_S}(u^m, u^*) + [\gamma_l(u^*) + \gamma_{nl_S}(u_p, u^*)] : \mathbf{C} : [\gamma_l(u_p) + \gamma_{nl_S}(u^m, u_p)] \} d\Omega \quad (32)$$

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

$$F_p^{nl}(u^*) = - \int_{\Omega} \left\{ \sum_{i=1}^{p-1} S_i : \gamma_{nl_S}(u_{p-i}, u^*) + \sum_{i=1}^{p-1} [\gamma_{nl}(u_i, u_{p-i})] : \mathbf{C} : [\gamma_l(u^*) + \gamma_{nl_S}(u^m, u^*)] \right\} d\Omega \quad (33)$$

Discretization of Equation (31) by using FEs leads to a sequence of linear problems in the form [39]

$$\text{Order 1} \left\{ \begin{array}{l} K_t u_1 = \bar{\lambda}_1 f \\ u_1^T u_1 + \bar{\lambda}_1^2 = 1, \end{array} \right. \quad (34)$$

$$\text{Order } p \left\{ \begin{array}{l} K_t u_p = \bar{\lambda}_p f + f_p^{nl}(u_i) \quad i < p \\ u_p^T u_1 + \bar{\lambda}_p \lambda_1 = 0, \end{array} \right. \quad (35)$$
where $K_t$ denotes the tangent stiffness matrix associated with Equation (32), common to the prob-
lems at different orders $p$. It is the same as the one applied in a classical iterative algorithm such
as Newton-Raphson (in the first iteration). In the aforementioned equation, $u_p$ is the discretized
form of the displacement field at order $p$, $f$ is the loading vector, and $f_p^{nl}$ represents the dis-
cretized form associated with $F_p^{nl}(u^*)$ in Equation (33), which at order $p$ only depends on the
values of $u_i$, $i < p$.

4. A COMBINED PGD–ANM APPROACH TO HYPERELASTICITY

The combination of the two previously introduced tools, namely the PGD approach to parametrized
problems (in this case the position of the load is the parameter) and the ANM for a consistent lin-
earization of the weak form of the problem gives rise to a particularly useful formulation. In it,

<figure><img src="image_2.png" /><figcaption>Figure 1. Load-displacement curve (in terms of the load parameter $\lambda$) for one particular load position of the pinched cube problem. Different solutions for different orders of expansion ($p = 1, \dots, 5$) compared with the FE solution by employing Newton-Raphson algorithms.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2. Deformed configuration of the cube for one particular load position.</figcaption></figure>
the displacement field is approximated as a series expansion around the last equilibrium point, whereas each term of the series is considered to be approximated by a finite sum of separated functions,

$$ \mathbf{u}^{m+1} = \mathbf{u}^m + a \sum_{i=1}^{n_1} (\mathbf{F}_i^1 \circ \mathbf{G}_i^1) + a^2 \sum_{j=1}^{n_2} (\mathbf{F}_j^2 \circ \mathbf{G}_j^2) + \dots + a^p \sum_{l=1}^{n_p} (\mathbf{F}_l^p \circ \mathbf{G}_l^p). \quad (36) $$

This gives rise to a series of problems of the form (34)–(35), within which a traditional PGD problem is solved for functions $\mathbf{F}_i^j(\mathbf{x})$ and $\mathbf{G}_i^j(\mathbf{s})$ at each order of the expansion $j$. No further modification of the method is necessary, resulting in a series of standard PGD problems that can be solved by employing any of the available non-linear solvers. Here, as in many of our previous works, we have employed a fixed point algorithm, similar to the one sketched in Sections 2.1–2.2 before.

We study the behavior of the proposed technique by means of a series of benchmark problems in Section 5.

<figure><img src="image_4.png" /><figcaption>Figure 3. Deformed configuration of the beam for the load positioned at the beam tip. Note that no artificial gain in volume can be perceived.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4. Load-displacement curve (in terms of the load parameter $\lambda$) for one particular load position of the beam bending problem. Different solutions for different orders of expansion ($p = 1, \dots, 5$) compared with the FE solution by employing Newton-Raphson algorithms.</figcaption></figure>
5. NUMERICAL EXAMPLES

5.1. Kirchhoff–Saint Venant material

We will consider two sets of examples. The first one is composed by three different benchmark
tests over Kirchhoff–Saint Venant materials. Despite its simplicity and well-known limitations,
Kirchhoff–Saint Venant approaches are very useful in the field of real-time simulation, because they
provide a good compromise between realism in the deformation and computational cost [5, 10].

5.1.1. Pinched cube. A unit cube modeled by 3 × 3 × 3 nodes and a tetrahedral mesh is con- sidered. Young’s modulus of 1 MPa and Poisson’s coefficient of 0.25 are assumed. The cube is loaded by a vertical force (0.01 N) acting at any point of the top face. Results obtained with the before presented PGD–ANM approach are compared with traditional FE approaches, solved by standard Newton–Raphson linearization strategies. For one particular position of the load (one cor- ner), the load-displacement curve is shown in Figure 1, whereas its deformed configuration is shown in Figure 2.

In general, as was the case in previous references such as [18], expansions up to order six are
enough to obtain a reasonable accuracy, despite the fact that in the ANM community, much higher
order expansions are usually employed on the order of 15 [38]. The number of modes (separated

<figure><img src="image_6.png" /><figcaption>Figure 5. Geometry of the finite element model for the liver.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6. Load-displacement curve (in terms of the load parameter λ ) for one particular load position of the liver palpation problem. Different solutions for different orders of expansion ( p = 1 , . . . , 6) compared with the FE solution by employing Newton–Raphson algorithms. Noticeably, the Newton–Raphson algorithm did not converge when solving the case λ = 1 .</figcaption></figure>
functions) employed at each expansion term in this example was 17, 17, 11, 2, 1 and 1, respectively, for terms of orders one to six. The accuracy of the approach is noteworthy, despite the fact that terms 5 and 6 of the expansion (solution for the sixth order expansion is not depicted in Figure 1, because it is practically indistinguishable from that of order five) are obtained with only one couple $F(x) \circ G(s)$, that is, $n_5 = n_6 = 1$ in Equation (36).

<figure><img src="image_8.png" /><figcaption>Figure 7. Vertical displacement field of the Kirchhoff–Saint Venant liver for one particular position of the load.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8. Functions $F_i^1(x)$, for $i = 1, 2, 3$, and $20$, respectively, for the first order expansion, in the simulation of the Kirchhoff–Saint Venant liver.</figcaption></figure>
**5.1.2. Beam bending.** In this case, we consider the problem of beam bending. This simple test is among the most popular ones in the field of real-time simulation, because it readily shows a great divergence from physical results if a poor formulation is used [43]. In fact, if linear elasticity formulations are employed, a great gain of volume is observed, leading the observer to perceive a clearly non-physical result, see Figure 3.

In this case, we consider a squared cross-section beam, with 40 × 40 × 400 mm. Young's modulus was assumed to be 209,000 MPa, whereas Poisson coefficient was set to 0.3. A load of 10⁶ N can be applied at any point of the upper face of the beam. Under these conditions, a comparison was established between the tip displacement obtained for the load position at the beam rightmost side with the value obtained by using standard FE analysis and Newton-Raphson iterations to solve the non-linear equations.

Results for different expansion order are shown in Figure 4. For a fifth order expansion, almost no difference can be perceived between the FE result and the reduced PGD-AnM result. For this example, the number of separated functions necessary at each expansion order was 162, 102, 79, 161, 24, and 8, respectively. As stopping criterion, we employed a relative one, which stops the fixed point algorithm if the new functional pair contributes less than 10⁻⁴ times the initial pair of functions, and a general one in which the algorithm is stopped if the modulus of the new functional pair is less than 10⁻¹⁵.

In general, results are below 5% error with respect to the target FE solution of the problem. In any case, a higher number of modes for each term, or a higher number of terms in the expansion, has been shown to converge to the right solution. It is true, in general, that the number of off-line computed modes depends on the error estimator considered (see the works by Huerta [28] or Chamoin

<figure><img src="image_10.png" /><figcaption>Figure 9. Functions $F_i^4(x)$, for $i = 1, 2, 3$, and $20$, respectively, for the fourth order expansion, in the simulation of the Kirchhoff-Saint Venant liver.</figcaption></figure>
and Ladeveze [29]). At present, there are no efficient and robust error estimators applicable to general multiparametric non-linear models. For this reason, we computed many terms (remember that this calculation is carried out off-line and only once) and then we compared the solution for some particular choices of the parameters with the one obtained with FEs to check the convergence and solution accuracy.

5.1.3. Palpation of the liver. One of the most typical examples in real-time applications is that of liver palpation in a laparoscopic virtual surgery procedure. The liver is the biggest gland in the human body, after the skin. Liver geometry has been obtained from the SOFA project [43] and post-processed to obtain a mesh composed by 2853 nodes and 10,519 tetrahedra, whose geometry is shown in Figure 5. The anterior surface of the liver is considered free, whereas the posterior one was assumed to be supported over different organs (it is connected to the diaphragm by the coro- nary ligament, for instance). The inferior vena cava travels along the posterior surface, and the liver is frequently assumed clamped at that location. Although the assumed boundary conditions are not strictly correct from a physiological point of view, our main interest is to show that the model can be solved under real-time constraints with reasonable accuracy.

Here, the human liver is considered as a Kirchhoff–Saint Venant material with *E* = 0.17 MPa and *ν* = 0.48 [44]. A vertical load of 5 N is considered at any point of the anterior surface of the liver. Again, for comparison purposes, results are evaluated at a particular position of the load and compared with those obtained with a traditional FE analysis and standard Newton–Raphson iterative schemes. It can be noticed how the sixth order expansion gives very reasonable results if compared with that of the FE model (Figure 6).

<figure><img src="image_11.png" /><figcaption>Figure 10. Functions $G_i^{(1)}(s)$, for $i = 1, 2, 3$, and $20$, respectively, for the first order expansion, in the simulation of the Kirchhoff–Saint Venant liver.</figcaption></figure>
The vertical displacement field for that particular position of the load is depicted in Figure 7. To check the overall behavior of the technique, functions $F_i$ and $G_i$ (Equation (36)) for $i = 1, 2, 3$, and 20 (terms of first and fourth order expansions) are shown in Figures 8–11.

## 5.2. Neo-Hookean behavior

Extension of the before presented technique to neo-Hookean materials [45] is in principle straightforward, although lengthy (this is not important, in fact, because the parametric calculation is carried out off-line). The major difference with the Kirchhoff–Saint Venant model is the presence of material non-linearities, in addition to the geometrical ones. A POD–ANM approach has been suggested in a previous work of the authors, see [18], and is here extended to the before presented PGD framework.

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

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

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

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

<figure><img src="image_12.png" /><figcaption>Figure 11. Functions $G_i^4(s)$, for $i = 1, 2, 3$, and 20, respectively, for the fourth order expansion, in the simulation of the Kirchhoff–Saint Venant liver.</figcaption></figure>
In this case, the intricate expansion procedure becomes easier if we identify, as in [37], the asymptotic expansion with a Taylor series of the variables of interest, denoted by $U(a)$, in the vicinity of $a = 0$. Truncating at order $N$.

$$U(a) = U_0 + \sum_{p=1}^{N} U_p a^p, \qquad (39)$$

where $U_0 = U(0)$ and

$$U_p = \frac{1}{p!} \frac{d^p U}{da^p} \Big|_{a=0}. \qquad (40)$$

<figure><img src="image_13.png" /><figcaption>Figure 12. Solution for the neo-Hookean beam bending problem at different expansion orders.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 13. Solution for the neo-Hookean beam bending problem at different expansion orders.</figcaption></figure>
In this case, as in [17], we have selected the following variables to perform the expansion,

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

By performing the substitution of the aforementioned variables into the weak form of the problem (Equation (21)), we arrive to a problem entirely similar to that in Equations (34) and (35). The entire details are provided, for completeness, in Appendix 1.

### 5.2.1. Neo-Hookean beam under bending.
We reproduce here the problem in Section 5.1.2 but considering a neo-Hookean constitutive model. The deformed configuration of the beam at expansion orders one to four is depicted in Figure 12. Note how the first order expansion (linear approach) shows a tremendous gain in volume that renders the simulations clearly non-physical. Again, expansions up to orders four to six were judged sufficient to obtain a good approximation to the reference solution (Figure 13). Obviously, higher accuracy can be obtained by increasing even more the expansion order. Up to $p = 15$ is a typical value of the expansion order in the ANM literature.

### 5.2.2. Palpation of a neo-Hookean liver.
The same procedure has been applied to the problem in Section 5.1.3 but now considering neo-Hookean behavior. The neo-Hookean law in Equation (37) has been now particularized to $E = 0.17$ MPa and $\nu = 0.48$. As in previous examples, a PGD-ANM solution has been obtained and compared with a standard FE solution at a particular node. To this end, standard Newton-Raphson procedures have been employed for the solution of the resulting non-linear system of equations. The load-displacement curve for this particular node is shown in Figure 14.

In this case, the observed agreement between the PGD-ANM solution and that of the FE model is even higher than in the previous example. For the fourth order expansion, the agreement between the predicted load-displacement curves is almost exact.

Remember that our approach is based upon an off-line/online procedure such that, once the off-line computation has been carried out, its solution is stored in the form of a series of 1D vectors that are evaluated in real time very efficiently. The solution thus computed is a multidimensional one, that is particularized online. This procedure is sketched in Figure 15.

<figure><img src="image_15.png" /><figcaption>Figure 14. Comparison between the PGD-ANM solution, at different expansion orders, and that for a standard finite element solution for a load at a particularized position. Neo-Hookean behavior.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 15. The result of the off-line PGD simulation is a multidimensional solution that depends on the load position. Thus, evaluating the solution for a prescribed load is entirely similar to performing a cut on a hyper volume $\Omega \times \bar{\Gamma}$. In the figure, the red line on the blue liver represents approximately the path followed by the load over the surface of the liver. The displacement is magnified by a factor of 2.0.</figcaption></figure>

## 6. CONCLUSIONS

The problem of PGD approaches to non-linear problems has been addressed here, with a particu- lar emphasis on hyperelasticity. The development of suitable linearizations for complex non-linear problems formulated under a PGD framework has been a major focus of attention for our commu- nity in the last years. First, because consistent linearizations of the complex weak forms appearing in these formulations are far from being readily available. Previous works include explicit approaches to these problems [35] but are difficult to generalize to arbitrary hyperelastic models or constitutive laws. The presented technique could be applied to any hyperelastic constitutive law, and, in view of existing results in the ANM, possibly to other non-linearities such as plasticity.

Here, we have proposed a combined PGD–ANM approach. The main advantages of this approach is that it produces a series approach to the solution that involves the same tangent operator for all the terms in the expansion. Thus, costly stiffness matrix updating procedures, which constitute nowa- days to one of the main problems in the model order reduction community, are avoided. It has been shown how the proposed method provides with a very accurate solution to complex problems in the form of a metamodel that can be evaluated at very high feedback rates (even at real time rates on the order of 1 kHz) with very little computational cost. This is precisely one of the most attractive features of the just-proposed method. It combines an off-line/online procedure by which a truncated expansion of the multidimensional form of the displacement is obtained and stored in memory. After, these modes are evaluated in the online procedure of the method at extremely fast feedback rates. We believe that this philosophy opens new possibilities in the field of real-time simulation that deserve to be explored.
Other problems remain open, however. Notably, optimality of PGD approaches (i.e., under what circumstances a priori PGD modes are equivalent to a posteriori POD or SVD modes) is not well understood. It has been noticed in the examples throughout this work that, unlike those in previous explicit approaches, PGD-ANM modes seem not to be optimal. They are highly oscillating and numerous for a prescribed tolerance. The solution is simple, however. It suffices to obtain the SVD modes that show most of the energy of the system to obtain a very compact representation of the solution (the so-called projected PGD, see [46]). But in any case, it seems pertinent to work in discerning what are the key ingredients of optimality in a PGD approach. This constitutes one of our current efforts of research.

In any case, although not optimal, the proposed method provides with a very competitive solution to highly demanding problems in applied sciences and engineering such as dynamic data-driven problems, real-time response even under very restrictive scenarios (haptic peripherals, for instance), simulation-driven control of structures and processes, and many others, where non-linear simulations are nowadays standard in industrial practice.

## APPENDIX A: DERIVATION OF THE TANGENT STIFFNESS MATRIX FOR THE NEO-HOOKEAN CASE

For the neo-Hookean case explained in Section 5.2, the tangent stiffness matrix takes the following form:

$$ K_t = \int_{\Omega_0} (\mathbf{B}^T \mathbf{D} \mathbf{B} + \mathbf{G}^T \tilde{\mathbf{S}}_0 \mathbf{G}) \, d\Omega, \quad (A.1) $$

where

$$ \mathbf{D} = \lambda \left( \frac{1}{J_0^2} C_0^{-1} M_0^T \right) + 2(\mu - \lambda \ln J_0) \left( \frac{1}{J_0^2} (C_0^{-1} M_0^T) - \tilde{C}_0 \right) \quad (A.2) $$

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_0^{-1}$ and contains components of $C_0$, arranged in a particular way.

The geometrical non-linearities are included in the matrices $\mathbf{B}$, $\mathbf{G}$, and $\tilde{\mathbf{S}}_0$. $\mathbf{B}$ represents the usual strain-displacement matrix, $\mathbf{G}$ relates the nodal displacements $\mathbf{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 [39]).

In the right-hand side of Equation (35), the non linear load vector $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

$$ f_p^{nl} = \int_{\Omega_0} \left( \mathbf{B}^T (\mathbf{S}_p^{nlmat} + \mathbf{S}_p^{nlgeom}) + \mathbf{G}^T \mathbf{S}_p^* \right) d\Omega \quad (A.3) $$

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

$$ \begin{aligned} \mathbf{S}_p^{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). \end{aligned} \quad (A.4) $$
In this equation, $CC_0$ represents the cofactor matrix of $C_0$, and $RCC_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,

$$CC_p = \tilde{C}_0 C_p + RCC_p, \qquad (A.5)$$

$$RCC_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}. \qquad (A.6)$$

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 + RJ_p. \qquad (A.7)$$

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

$$RS_p = \sum_{r=1}^{p-1} Y_r C_{p-r}^{-1}. \qquad (A.8)$$

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

$$RC_p^{-1} = \sum_{r=1}^{p-1} Z_r CC_{p-r}. \qquad (A.9)$$

Finally, it is necessary to expand $Y = \ln J$ and $Z = J^{-2}$ by using 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, \qquad (A.10)$$

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,$$

...

At this point, a procedure entirely similar to that of Equation (36) is performed, leading to an equivalent expression of the displacement in terms of a power series in parameter $a$, where each term is composed by a finite sum of separable functions.
ACKNOWLEDGEMENT

The authors would wish to thank the Spanish Ministry of Economy and Innovation (grant number, CICYT-DPI2011-27778-C02-01) for their support.

REFERENCES

1. Karhunen K. Über lineare methoden in der wahrscheinlichkeitsrechnung. *Annales Academiae Scientiarum Fennicae Series A1 Mathematical Physics* 1946; **37**:1–79.

2. Loève MM. *Probability Theory*, 3rd ed., The University Series in Higher Mathematics. Van Nostrand: Princeton, NJ, 1963.

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

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. Barbič J, James DL. Real-time subspace integration for St. Venant–Kirchhoff deformable models. *ACM Transactions on Graphics* (SIGGRAPH 2005) 2005; **24**(3):982–990.

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

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

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

9. Taylor ZA, Crozier S, Ourselin S. A reduced order explicit dynamic finite element algorithm for surgical simulation. *IEEE Transactions on Medical Imaging* 2011; **30**(9):1713–1721.

10. Taylor ZA, Ourselin S, Crozier S. A reduced order finite element algorithm for surgical simulation. In *Engineering in Medicine and Biology Society (EMBC), 2010 Annual International Conference of the IEEE*, Buenos Aires, Argentina, September 4–31, 2010; 239–242.

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

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

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

14. Barrault M, Maday Y, Nguyen NC, Patera AT. An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations. *Comptes Rendus Mathematique* 2004; **339**(9):667–672.

15. Chaturantabut S, Sorensen DC. Nonlinear model reduction via discrete empirical interpolation. *SIAM Journal of Scientific Computing* 2010; **32**:2737–2764.

16. Nguyen NC, Patera AT, Peraire J. A 'best points' interpolation method for efficient approximation of parametrized functions. *International Journal for Numerical Methods in Engineering* 2008; **73**(4):521–543.

17. Niroomandi S, Alfaro I, Cueto E, Chinesta F. Accounting for large deformations in real-time simulations of soft tissues based on reduced-order models. *Computer Methods and Programs in Biomedicine* 2012; **105**(1):1–12.

18. Niroomandi S, Alfaro I, Cueto E, Chinesta F. Model order reduction for hyperelastic materials. *International Journal for Numerical Methods in Engineering* 2010; **81**(9):1180–1206.

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

20. Amsallem D, Farhat Ch. An interpolation method for adapting reduced-order models and application to aeroelasticity. *AIAA Journal* 2008; **46**:1803–1813.

21. Ladeveze P. *Nonlinear Computational Structural Mechanics*. Springer: N.Y., 1999.

22. Ammar A, Mokdad B, Chinesta F, Keunings R. A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids. Part II: transient simulation using space-time separated representations. *Journal of Non-Newtonian Fluid Mechanics* 2007; **144**:98–121.

23. Ammar A, Mokdad B, Chinesta F, Keunings R. A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids. *Journal of Non-Newtonian Fluid Mechanics* 2006; **139**:153–176.

24. Chinesta F, Ammar A, Cueto E. Recent advances in the use of the proper generalized decomposition for solving multidimensional models. *Archives of Computational Methods in Engineering* 2010; **17**(4):327–350.

25. Chinesta F, Ladeveze P, Cueto E. A short review on model order reduction based on proper generalized decomposition. *Archives of Computational Methods in Engineering* 2011; **18**:395–404.

26. Ladeveze P, Passieux J-C, Neron D. The Latin multiscale computational method and the proper generalized decomposition. *Computer Methods in Applied Mechanics and Engineering* 2010; **199**(21-22):1287–1296.

27. Le Bris C, Lelièvre T, Maday Y. Results and questions on a nonlinear approximation approach for solving high-dimensional partial differential equations. *Constructive Approximation* 2009; **30**:621–651. 10.1007/s00365-009-9071-1.
28. Ammar A, Chinesta F, Diez P, Huerta A. An error estimator for separated representations of highly multidimensional models. *Computer Methods in Applied Mechanics and Engineering* 2010; **199**(25-28):1872-1880.

29. Ladeveze P, Chamoin L. On the verification of model reduction methods based on the proper generalized decomposition. *Computer Methods in Applied Mechanics and Engineering* 2011; **200**(23-24):2032-2047.

30. Ghnatios Ch, Chinesta F, Cueto E, Leygue A, Poitou A, Breitkopf P, Villon P. Methodological approach to efficient modeling and optimization of thermal processes taking place in a die: application to pultrusion. *Composites Part A: Applied Science and Manufacturing* 2011; **42**(9):1169-1178.

31. Ghnatios Ch, Masson F, Huerta A, Leygue A, Cueto E, Chinesta F. Proper generalized decomposition-based dynamic data-driven control of thermal processes. *Computer Methods in Applied Mechanics and Engineering* 2012; **213-216**(0):29-41.

32. Pruliere E, Chinesta F, Ammar A. On the deterministic solution of multidimensional parametric models using the Proper generalized decomposition. *Mathematics and Computers in Simulation* 2010; **81**(4):791-810.

33. Bognet B, Bordeu F, Chinesta F, Leygue A, Poitou A. Advanced simulation of models defined in plate geometries: 3D solutions with 2D computational complexity. *Computer Methods in Applied Mechanics and Engineering* 2012; **201-204**(0):1-12.

34. Relun N, Neron D, Boucard P-A. A model reduction technique based on the PGD for elastic-viscoplastic computational analysis. *Computational Mechanics* 2013; **51**:83-92.

35. Niroomandi S, Gonzalez D, Alfaro I, Bordeu F, Leygue A, Cueto E, Chinesta F. Real time simulation of biological soft tissues: a PGD approach. *International Journal for Numerical Methods in Biomedical Engineering* 2013; **29**(5):586-600.

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

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

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

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

40. Gonzalez D, Masson F, Poulhaon F, Cueto E, Chinesta F. Proper generalized decomposition based dynamic data driven inverse identification. *Mathematics and Computers in Simulation* 2012; **82**:1677-1695.

41. Bro-Nielsen M, Cotin S. Real-time volumetric deformable models for surgery simulation using finite elements and condensation. *Computer Graphics Forum* 1996; **15**(3):57-66.

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

43. Allard J, Cotin S, Faure F, Bensoussan P-J, Poyer F, Duriez C, Delingette H, Grisoni L. SOFA – an open source framework for medical simulation. In *Medicine Meets Virtual Reality (MMVR’15)*, Long Beach, USA, February 2007; 13-18.

44. Delingette H, Ayache N. Soft tissue modeling for surgery simulation. In *Computational Models for the Human Body*, Ayache N (ed.), Handbook of Numerical Analysis (Ph. Ciarlet, Ed.) Elsevier: Juan-les-Pins, France, 2004; 453-550.

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

46. Nouy A. A priori model reduction through proper generalized decomposition for solving time-dependent partial differential equations. *Computer Methods in Applied Mechanics and Engineering* 2010; **199**(23-24):1603-1626.