Francisco Chinesta

Ph Lorong

Elias Cueto

NEW TRENDS ON THE NUMERICAL MODELING OF FORMING PROCESSES

Keywords: Numerical modeling, Meshless methods, Forming Processes, Multiscale

In this work, we explore some recent numerical developments able to perform fast and accurate simulations of complex systems, as the ones encountered in forming processes. In particular, we will focus on the use of meshless techniques and a discussion on the perspectives introduced by multiscale modeling.

INTRODUCTION

Numerical techniques for solving the complex models encountered in forming processes have experienced a remarkable progress in recent years. In convectional and nonconventional forming processes different difficulties are encountered from the modeling and its associated numerical treatment: (i) large thermomechanical transformation that usually involve large geometrical transformations; (ii) strong multi-physics coupling; (iii) localization in space and time with the consequent effect on the discrete resolution (in both space and time); (iv) multiscale behavior (that could be defined from the quantum level to the macroscopic one); (v) the large resulting discrete system, usually non-linear and that involves millions of degree of freedom; (vi) existence of evolving discontinuities; among many others. Obviously, today does not exist a technique able to alleviate all the just referred difficulties, but some progresses have been made for developing new strategies able to circumvent some of them, and that potentially could be integrated into the more experienced techniques in order to expand the domain of applicability of numerical simulation. A novel technique that seems an appealing choice for treating this kind of models is the one known as Natural Element Method that will be revisited in the next section, before illustrating its application in the forming process framework. Further considerations in non-conventional machining are addressed in the last section.

NATURAL ELEMENT METHOD

The NEM [START_REF] Elias | Overview and recent advances in Natural Neighbour Galerkin methods[END_REF] is a Galerkin procedure based on the natural neighbor interpolation scheme. This interpolation scheme is the key ingredient of the method here proposed. It relies on the concepts of Voronoi diagrams and Delaunay triangulations (see Figure 1), to build Galerkin trial and test functions. These are defined as the Natural Neighbor coordinates of the point under consideration, that is, with respect to Figure 2, the value in the point x of the shape function associated with the node 1 is:

abcd abfe A A x = ) ( 1 φ (1)
In addition, the NEM has other interesting properties such as linear consistency and smooth shape functions. These functions are dependent on the position and density of nodes, leading to standard FE constant strain triangle shape functions, bilinear shape functions or rational quartic functions in different situations (see Fig. 3 for a typical shape function). These properties permit an exact reproduction of linear displacement fields on the boundary of convex domains. In the non-convex case other possibilities exist [START_REF] Elias | Overview and recent advances in Natural Neighbour Galerkin methods[END_REF][START_REF] Julien | A new extension of the Natural Element Method for non convex and discontinuous problems: the Constrained Natural Element Method[END_REF].

SIMULATION PROCEDURE

This strategy allows to proceed with the same cloud of nodes from the beginning of the simulation to its end, because even is some Delaunay triangles become very distorted during the simulation, the solution accuracy does not depend significantly on the diagrams quality.

In this way remeshing are no more required, and consequently field projection between the old and new meshes avoided. It is well known that this remeshing procedure involves a non negligible amount of numerical diffusion. The solution strategy is very simple. Let

( ) ( ) ( , u t f t ℑ = x ) , x (2)
be the model (system of partial differential equations) to be solved. One can proceed, as in the finite element framework, to consider the associated variational formulation

( ) ( ) ( ) 0, t u f d Ω Ψ ℑ - Ω = ∀Ψ ∫ (3)
The initial domain is denoted by in which we consider the Delaunay triangulation related to a cloud of N nodes with known initial positions. The unknown field can be approximated in a standard way:

0 Ω ( ) 1 1 N N u u u φ = + + x x L ( ) φ (4)
using the Sibson's shape functions, and then introduced in the variational formulation (3) leading to a discrete linear system whose solution results in the nodal values. Then, using (4) the solution can be computed everywhere. Now, the domain is updated, by advancing all the nodes from their velocities:

(5)

i i i t ← + Δ x x v
the Delaunay and Voronoi diagrams recomputed, and the solution related to the new time step computed. This procedure continues until the end of the process simulation.

All the internal variables, describing the thermo-mechanical history, are associated with the nodes and then their evolutions are computed at the nodal positions by performing a local calculation.

From the last analysis it seems natural that this strategy could be applied either in the solid or in the fluid mechanics frameworks. Moreover, its use simplifies significantly the solid-fluid coupling because both are defined in an updated Lagrangian framework. Its generalization to the 3D case is quite simple.

In figure 4 we depict some snapshots concerning the elastoviscoplastic deformation of a plate. This computation has been performed by using the same cloud of nodes during the entire transformation. In some cases, the accuracy of the solution needs for higher nodal densities in some regions, as is the case of models involving localization. Figure 5 depicts the appearance of an adiabatic shear band after impact, by assuming again an elastoviscoplastic constitutive equation in the context of an explicit dynamics, and by using an appropriate automatic procedure to adapt the nodal density. Other specificity of the proposed numerical technique is its ability to label the nodes, that is, if we have two materials that are mixing we could label the nodes belonging at the initial time to the first one by assigning a color (red for example) being the others labeled with the blue color. As we can proceed with the same cloud of nodes in the entire time interval, each node keeps its color allowing to quantify the mixing at the end of the simulation. Figure 6 and 7 shows the initial and final snapshots related to the numerical simulation of the friction stir welding process of two plates, whose initial nodes are labeled with the red and blue colors. We can observe that a nice material mixing is obtained in the interface neighborhood as the welding process (pin rotation) progress. Other appealing property of this technique is its ability to capture free interfaces (in particular the alpha-NEM version widely described in [START_REF] Elias | Overview and recent advances in Natural Neighbour Galerkin methods[END_REF]). Thus, the combination of an the alpha shape constructor and the natural neighbor interpolation allows to define a robust technique for simulating free surface problems as usually encountered in fluid mechanics. Some welding processes need supply an amount of material that usually consists in a drop that impacts a solid or a liquid phase. Figure 8 depicts the impact of a drop on a solid and rigid surface. This technique has been successfully applied in the simulation of laser coating processes, arising, for instance, in superconductor texturing processes or ceramic tiles processing, illustrated in Fig. 9.

Fig. 8. Drop impact on a rigid

Fig. 9. Typical apparatus to texture a superconductor

The resulting Stefan problem (thermal model including change of phase) can be accurately solved by using the natural element method as described in [START_REF] David | Natural neighbour strategies for the simulation of laser surface coating process[END_REF]. The evolution of the temperature field when the heat source is moving with uniform velocity on the plate surface is illustrated in Fig. 10.

Fig. 10. Tempeture field evolution

The melted zones at two different times are depicted in figure 11 and compared in figure 12 with the experimental results.

MICRO AND NONOMETRIC MODELING

Fine descriptions of the structure and mechanics of materials at the micro, nano and subnanometric scales involve numerous computational challenges.

Accurate descriptions of such materials need for a multiscale approach and the definition of pertinent bridges between these different scales.

Fig. 11. Melted zones

Fig. 12. Transversal cut of a textured Bi

The finest description concerns the atomic level where quantum mechanics leads to molecular dynamics simulations. Even if the resolution of the Schrödinger equation could lead to the establishment of accurate interatomic potentials, that could be used in molecular dynamics simulations, the curse of dimensionality that this equation induces represents even today a computational challenge.

The large number of atoms to be included in MD simulations, the extremely small time step imposed by the numerical stability and the approximated (at best only quantum inspired) interaction potentials, are the main drawbacks related to its use.

The next description scale introduces some simplifications leading to a coarse grained molecular dynamics, the DPD (dissipative particle dynamics) being one of such approaches. In that approach atoms are substituted by aggregates (or clusters) of atoms that interact from a variety of "empirical" interaction potentials.

Coarser descriptions focus on some entities by substituting the other entities by their effects usually modeled by random distributions.

This is the scale of Brownian dynamics simulations where the Langevin equation is a prototype.

Other descriptions come back to the continuous medium.

In those approaches the physical entities are substituted by a distribution defined in the physical space, in time and in the conformation coordinates: e.g. velocities in the case of models related to systems composed by charged of non-charged particles (Boltzmann or Vlasof-Poisson-Boltzmann equations, the last one modeling quantum relativistic plasma); or the molecular conformation in the case of macromolecular systems (multi-bead-springs or multi-beadrods models) leading to the Fokker-Planck formalism.

This approach, despite its mathematical simplicity, introduces a density function that is defined in a multidimensional space, and then the associated partial differential equations must be solved in a multidimensional domain.

To circumvent the curse of dimensionality that these high-dimensional partial differential equations induce, stochastic techniques have been applied intensively in the last decade. Some improvements have been proposed, the Brownian Configurations Fields being one of such approaches.

Recently some incipient techniques based on sparse grids or those based on separated representations have allowed solving models defined in highly multidimensional spaces [START_REF] Amine | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids[END_REF].

The coarsest scale in the description of such material is the one related to the macroscopic flow.

Micro-macro Lagrangian/Eulerian approaches have been intensively applied. Today meshless approaches, as the one described in the first part of this work, are opening interesting possibilities because its natural Lagrangian character that makes possible simulate flows using the same cloud of nodes that are advected with the material velocity.

The key point in that multiscale approaches is the definition of appropriate bridges between these scales. In the case of macromolecular systems, where the inertia effects can be (and are) usually neglected, the bridge is simple (Kremer's rule). For other systems the establishment of such bridges is at present an open problem.

Fig. 2 .Fig. 2 .Fig. 3 .

 223 Fig. 2. Delaunay triangulation and Voronoi diagram

Fig. 4 .Fig. 7 .

 47 Fig. 4. Large elastoviscoplastic transformation


Siamak Niroomandi

Icıar Alfaro

Elıas Cueto

Francisco Chinesta

Model order reduction for hyperelastic materials

Keywords: model order reduction, proper orthogonal decomposition, asymptotic numerical method, Kirchhoff-Saint Venant, neo-Hookean material

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.

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 [START_REF] Idelsohn | A reduction method for nonlinear structural dynamics analysis[END_REF][START_REF] Krysl | Dimensional model reduction in non-linear finite element dynamics of solids and structures[END_REF] and references therein, to name a few. But similar techniques are employed in many other branches of Engineering and Applied Sciences: real-time simulation [START_REF] Jernej | Real-time subspace integration for St. Venant-Kirchhoff deformable models[END_REF], Chemical Engineering [START_REF] Park | The use of the Karhunen-Loève decomposition for the modeling of distributed parameter systems[END_REF], turbulence [START_REF] Sirovich | Turbulence and the dynamics of coherent structures coherent structures part I: coherent structures[END_REF], weather forecast [START_REF] Lorenz | Empirical Orthogonal Functions and Statistical Weather Prediction[END_REF] and many others, linking back to the seminal works of Karhunen [START_REF] Karhunen | Uber lineare methoden in der wahrscheinlichkeitsrechnung[END_REF] and Loève [START_REF] Loève | Probability Theory (3rd edn). The University Series in Higher Mathematics[END_REF].

Recent applications of model reduction techniques for the analysis of non-linear solids and structures include the works by Marsden and colleagues [START_REF] Krysl | Dimensional model reduction in non-linear finite element dynamics of solids and structures[END_REF] and the works of Ryckelynck et al. for 'a priori' model reduction (see [START_REF] Ryckelynck | On the a priori model reduction: overview and recent developments[END_REF][START_REF] Ryckelynck | A priori model reduction method for the optimization of complex problems[END_REF][START_REF] Ryckelynck | A priori hyperreduction method: an adaptive approach[END_REF]). Also, recently, new extensions of the method have been developed, giving rise to the so-called non-linear dimensionality reduction [START_REF] Tenenbaum | A global framework for nonlinear dimensionality reduction[END_REF].

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 [START_REF] Ryckelynck | A priori hyperreduction method: an adaptive approach[END_REF] 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 [START_REF] Niroomandi | Real-time deformable models of non-linear tissues by model reduction techniques[END_REF].

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 realtime applications, for instance, in which results must be obtained at a frequency of 500 Hz [START_REF] Niroomandi | Real-time deformable models of non-linear tissues by model reduction techniques[END_REF]). 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 [START_REF] Ryckelynck | A priori hyperreduction method: an adaptive approach[END_REF], 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 [START_REF] Cochelin | Asymptotic-numerical methods and Padé approximants for non-linear elastic structures[END_REF][START_REF] Abichou | Asymptotic numerical method for problems coupling several nonlinearities[END_REF]. This technique was developed in [START_REF] Yvonnet | A model reduction method for the post-buckling analysis of cellular microstructures[END_REF] 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 nonlinearities. 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 [START_REF] Yvonnet | A model reduction method for the post-buckling analysis of cellular microstructures[END_REF] 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.

STATE-OF-THE-ART: MODEL REDUCTION AND ASYMPTOTIC EXPANSIONS

The Karhunen-Loève decomposition

We assume that the evolution of a certain field u(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(x i , t n ) ≡ u n i . The same can be wrote introducing a time discretization u n (x) ≡ u(x, t = n t); ∀n ∈[1, . . . , P]. The main idea of the Karhunen-Loève (K-L) decomposition is how to obtain the most typical or characteristic structure (x) among these u n (x)∀n. This is equivalent to obtaining a function (x) that maximizes defined by

= n=P n=1 i=N i=1 (x i )u n (x i ) 2 i=N i=1 ((x i )) 2 (1)
The maximization ( = 0) leads to

n=P n=1 i=N i=1 (x i )u n (x i ) j=N j=1 (x j )u n (x j ) = i=N i=1 (x i )(x i ), ∀ (2)
which can be rewritten in the form

i=N i=1 j=N j=1 n=P n=1 u n (x i )u n (x j )(x j ) (x i ) = i=N i=1 (x i )(x i ), ∀ (3)
Defining the vector / such that its i-component is (x i ), Equation (3) takes the following matrix form:

/T c/ = /T /, ∀ / ⇒ c/ = / (4)
where the two-point correlation matrix is given by

c ij = n=P n=1 u n (x i )u n (x j ) ⇔ c = n=P n=1 u n (u n ) T (5)
which is symmetric and positive definite. If we define the matrix Q containing the discrete field history:

Q = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ u 1 1 u 2 1 ••• u P 1 u 1 2 u 2 2 ••• u P 2 . . . . . . . . . . . . u 1 N u 2 N ••• u P N ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (6)
it is straightforward to verify that the matrix c in Equation ( 4) results

c = QQ T (7)
where the diagonal components are given by

c ii = (QQ T ) ii = j=P j=1 (u j i ) 2 (8)
Thus, the functions defining the most characteristic structure of u n (x) are the eigenfunctions k (x) ≡ / k associated with the highest eigenvalues.

A posteriori reduced-order modelling

If some direct simulations have been carried out, we can determine u(x i , t n ) ≡ u n i , ∀i ∈[1, . . . , N ], ∀n ∈[1, . . . , P], and from these the r eigenvectors related to the r -highest eigenvalues / k = k (x i ), ∀i ∈[1, . . . , N ], ∀k ∈[1, . . . ,r ] (with r 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(x i , t n ). For this purpose we need to define the matrix

A A = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 1 (x 1 ) 2 (x 1 ) ••• r (x 1 ) 1 (x 2 ) 2 (x 2 ) ••• r (x 2 ) . . . . . . . . . . . . 1 (x N ) 2 (x N ) ••• r (x N ) ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (9)
Now, we consider the linear system of equations resulting from the discretization of a partial differential equation (PDE) in the form KU = F [START_REF] Ryckelynck | A priori model reduction method for the optimization of complex problems[END_REF] 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

U = i=r i=1 i / i = Af (11)
it results

KU = F ⇒ KAf = F (12)
and multiplying both terms by A T it results

A T KAf = A T F (13)
which proves that the final system of equations is of low order, that is, the dimensions of A T KA are r ×r , with r N , and the dimensions of both f and A T F are r ×1.

Remark 1

Equation ( 13) can also be derived by introducing the approximation [START_REF] Ryckelynck | A priori hyperreduction method: an adaptive approach[END_REF] 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.

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 KA, 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 [START_REF] Ryckelynck | Hyper-reduction of mechanical models involving internal variables[END_REF] 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).

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 [START_REF] Cochelin | Asymptotic-numerical methods and Padé approximants for non-linear elastic structures[END_REF][START_REF] Abichou | Asymptotic numerical method for problems coupling several nonlinearities[END_REF]. 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 -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.

Problem formulation for Kirchhoff-Saint venant models

Here, we review the development made in [START_REF] Yvonnet | A model reduction method for the post-buckling analysis of cellular microstructures[END_REF] 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 0 whose boundary is denoted by . Essential and natural boundary conditions are applied to the non-overlapping portions of the boundary denoted by u and t . The deformed configuration of each point is given by the vector x, such that

x = X+u (14)
where the displacement field is denoted by u.

Following [START_REF] Cochelin | The asymptotic numerical method: an efficient perturbation technique for nonlinear structural mechanics[END_REF], we consider a linear and a non-linear terms for the Green-Lagrange strain tensor, E, in the form

E = 1 2 (F T F-1) = l (u)+ nl (u, u) (15)
where F = ∇u+I is the gradient of the deformation tensor and using the notation in [START_REF] Cochelin | The asymptotic numerical method: an efficient perturbation technique for nonlinear structural mechanics[END_REF],

l (u) = 1 2 (∇(u T )+∇(u)) (16a) nl (u, u) = 1 2 ∇(u T )∇(u) (16b)
Hyperelastic materials are based on the assumption of a particular strain-energy function, . Then the second Piola-Kirchhoff stress tensor S can thus be obtained by

S = * *E (17)
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 looks like ∇P+B = 0 in 0 [START_REF] Cochelin | The asymptotic numerical method: an efficient perturbation technique for nonlinear structural mechanics[END_REF] in which B is the body force. The boundary conditions of the body are defined by

u(X) = ū on u PN = t on t ( 19
)
where N is the unit vector normal to , t is an applied traction and 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 0

S : E d = t t• u d ∀ u ∈ H 1 ( ) ( 20
)
where in the above equation E is expressed by

E = 1 2 [F T ∇(du)+∇(du) T F] = l (du)+ nl S (u, du) (21)
where, in turn, nl S (u, u) is defined by

nl S (u, u) = nl (u, u)+ nl ( u, u) (22)
The technique is best viewed by considering the Saint-Venant Kirchhoff model, as done in [START_REF] Yvonnet | A model reduction method for the post-buckling analysis of cellular microstructures[END_REF][START_REF] Cochelin | The asymptotic numerical method: an efficient perturbation technique for nonlinear structural mechanics[END_REF], 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

= 2 (tr(E)) 2 + E : E (23)
where and are Lame's constants. The second Piola-Kirchhoff stress tensor can be obtained by

S = * (E) *E = C : E (24)
in which 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 [START_REF] Jernej | Real-time subspace integration for St. Venant-Kirchhoff deformable models[END_REF] 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).

Asymptotic numerical method for geometrically non-linear problems

In the ANM [START_REF] Cochelin | Asymptotic-numerical methods and Padé approximants for non-linear elastic structures[END_REF][START_REF] Abichou | Asymptotic numerical method for problems coupling several nonlinearities[END_REF] 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 (u n ; S n ; 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 are expanded in series prior to their introduction in the weak form of the problem

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ u n+1 (a) S n+1 (a) n+1 (a) ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ u n (a) S n (a) n (a) ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ + N p=1 a p ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ u p S p p ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ (25)
where (u p , S p , p ) are unknowns. Above, (u n+1 (a), S n+1 (a), 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 E(u) gives

E n+1 (a) = l (du)+ nl S (du, u n )+ n p=1 a p nl S (du, u p ) (26)
The series expansions of S gives in turn

S n+1 (a) = C : E n+1 (a) = C: nl (u n , u n )+ l (u n )+ n p=1 a p l (u p )+ nl S (u n , u p )+ p-1 i=1 nl (u i , u p-i ) (27)
and at order p we obtain

S p = C : l (u p )+ nl S (u n , u p )+ p-1 i=1 nl (u i , u p-i ) (28)
with W ext ( u) = t t• u d . 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, . . . , N ) takes the form

L(du, u n ) = p W ext ( u)+ F nl p ( u) (30)
with

L(du, u n ) = {S n : nl S (u n , u)+[ l ( u)+ nl S (u p , u)] : C : [ l (u p )+ nl S (u n , u p )]} d (31)
and where F nl p ( u) is equal to zero at order one and at order p it can be calculated as

F nl p ( u)=- p-1 i=1 S i : nl S (u p-i , u)+ p-1 i=1 [ nl (u i , u p-i )]:C:[ l ( u)+ nl S (u n , u)] d (32)
Discretization of Equation (30) by using finite elements leads to a sequence of linear problems in the form [START_REF] Cochelin | The asymptotic numerical method: an efficient perturbation technique for nonlinear structural mechanics[END_REF] Order 1

⎧ ⎨ ⎩ K t u 1 = 1 f u T 1 u 1 + 2 1 = 1 (33) Order p ⎧ ⎨ ⎩ K t u p = p f+f nl p (u i ), i< p u T p u 1 + p 1 = 0 ( 34
)
where 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, u p is the discretized form of the displacement field at order p, f is the loading vector and f nl p represents the discretized form associated with F nl p ( u) in Equation ( 32), which at order p only depends on the values of u i , i< p. The solution of these problems can be obtained as follows: at order 1

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ û = {K t } -1 f 1 = 1 √ ûT û+1 u 1 = 1 û (35)
at order p

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ u nl p = {K t } -1 f nl p p = -1 {u nl p } T u 1 u p = p 1 u 1 +u nl p (36)

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 [START_REF] Yvonnet | A model reduction method for the post-buckling analysis of cellular microstructures[END_REF] 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

u p = M m=1 / m m p = Af p ( 37
)
where f p are unknowns. So the new asymptotic expansion of u is expressed by

u n+1 (a) = A f n + N p=1 a p f p (38)
where (f n , n ) represents the previous converged solution. Introducing Equation (38) into Equation(30) results in

L(Af p , A f) = p W ext (A f)+ F nl p (A f) (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:

At order 1

⎧ ⎨ ⎩ A T K t Af 1 = 1 A T f f T 1 f 1 + 2 1 = 1 (40) At order p A T K t Af p = p A T f+f nl p f T p f 1 + p 1 = 0 (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:

At order 1

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ f = {A T K t A} -1 A T f 1 = 1 fT f+1 f 1 = 1 f û = A f, u 1 = 1 f 1 (42) At order p ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ f nl p = {A T K t A} -1 A T f nl p p = -1 {f nl p } T f 1 f p = p 1 f 1 +f nl p u p = Af p (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.

NON-LINEAR MODEL REDUCTION OF HYPERELASTIC MATERIALS

The extension of the technique explained before to other hyperelastic materials, in which nonlinearities 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 [START_REF] Bonet | Nonlinear Continuum Mechanics for Finite Element Analysis[END_REF].

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

= 2 (tr(C)-3)-ln J + 2 (ln J ) 2 (44)
where and 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 = * (E) *E = (I-C -1 )+ (ln J )C -1 (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 [START_REF] Cao | An improved iterative method for large strain viscoplastic problems[END_REF], 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 + N p=1 U p a p ( 46
)
where U 0 = U(0) and

U p = 1 p! d p U da p a=0 ( 47
)
In this case we have selected the following variables to perform the expansion:

U(a) = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ u(a)
S(a)

J 2 C -1 (a) ln J 2 (a) 1 J 2 (a) (a) ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (48)
By performing the substitution of the before-mentioned variables into the weak form of the problem, see Equation [START_REF] Cao | An improved iterative method for large strain viscoplastic problems[END_REF], we arrive to a problem entirely similar to that in Equations ( 33) and (34), in which the tangent stiffness matrix takes the form

K t = 0 (B T DB+G T S0 G) d (49)
where

D = 1 J 2 0 C -1 0 M T 0 +2( -ln J 0 ) 1 J 2 0 (C -1 0 M T 0 )-C0 (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, C0 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 S0 . B represents the usual strain-displacement matrix, G relates the nodal displacements u and the gradient of displacements vector, and, finally, S0 represents a matrix that contains the initial stresses (we have chosen the same notation as in [START_REF] Cochelin | The asymptotic numerical method: an efficient perturbation technique for nonlinear structural mechanics[END_REF]).

In the right hand side of Equation ( 34), the non-linear load vector f nl p 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 nl p = 0 (B T (S nlmat p +S nlgeom p )+G T S * p ) d (51)
As in the stiffness matrix, S nlgeom p and S * p represent the standard matrices found in the literature when ANM is used to solve geometrical non-linear problems with linear materials. S nlmat p takes into account the material behaviour:

S nlmat p = ( ln J 0 -) CC 0 R Z p - R J p J 4 0 + RCC p J 2 0 +RC -1 p + CC 0 J 2 0 RY p + R J p 2J 2 0 +RS p (52)
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 = C0 C p +RCC p ( 53
)
RCC p = p-1 r =1 ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ C r 22 C p-r 33 -C r 23 C p-r 23 C r 11 C p-r 33 -C r 13 C p-r 13 C r 11 C p-r 22 -C r 12 C p-r 12 C r 13 C p-r 23 -C r 12 C p-r 33 C r 13 C p-r 12 -C r 11 C p-r 23 C r 12 C p-r 23 -C r 13 C p-r 22 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (54)
Here, R J 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 T 0 C p + R J p ( 55
)
RS p collects the terms concerning the expansion of Y = ln J and C -1 :

RS p = p-1 r =1 Y r C -1 p-r (56)
RC -1 p collects the terms concerning Z = J -2 and cofactor matrix of C expansions:

RC -1 p = p-1 r =1 Z r CC p-r (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 = 1 2J 2 0 (J 2 ) p + RY p and Z p = -1 J 4 0 (J 2 ) p + R Z p ( 58
)
where

RY 1 = 0 RY 2 = -1 4J 4 0 (J 2 ) 2 1 ) RY 3 = 1 6J 6 0 (J 2 ) 3 1 +2 -1 4J 4 0 (J 2 ) 1 (J 2 ) 2 R Z 1 = 0 R Z 2 = 1 J 6 0 (J 2 ) 2 1 R Z 3 = -1 J 8 0 (J 2 ) 3 1 +2 1 J 6 0 (J 2 ) 1 (J 2 ) 2 . . .
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.

NUMERICAL EXAMPLES

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

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.

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.

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

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 hexaedral elements. Kirchhoff-Saint Venant and neo-Hookean behaviours are assumed in this example.

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, . . . , 6 have been obtained and depicted in Figure 7. 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 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 ≈ -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 = 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, ) curve shows great similarity for p = 2, . . . , 4 and a wide convergence radius.

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 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 [START_REF] Alastrué | Biomechanical modeling of refractive corneal surgery[END_REF]. Forces of 0.014 N 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 = 0.48, although a realistic model of corneal tissue is far more complicated [START_REF] Niroomandi | Real-time deformable models of non-linear tissues by model reduction techniques[END_REF].

In this example nine modes were applied that provide decent approximation. The solution has been obtained using ANM-POD for p = 1, . . . , 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 ( ) 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 = 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 proof of robustness). The use of truly incompressible (mixed) formulations is currently under investigation.

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.

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

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 [START_REF] Yvonnet | A model reduction method for the post-buckling analysis of cellular microstructures[END_REF] 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.

  p : l ( u)+ nl S (u n ,

Figure 1 .

 1 Figure 1. Geometry of the beam in traction.

Figure 2 .

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

Figure 3 .

 3 Figure 3. Beam in traction. Load-displacement results (in terms of ) for different orders of approximation, p, compared with the analytical solution. Kirchhoff-Saint Venant behaviour.

Figure 4 .

 4 Figure 4. Beam in traction. Load-displacement results (in terms of ) for different orders of approximation, p, compared with the analytical solution. Neo-Hookean behaviour.

Figure 5 .

 5 Figure 5. Six most important eigenmodes for the simulation of the neo-Hookean 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.

Figure 6 .Figure 7 .

 67 Figure 6. Hemisphere mesh.

Figure 8 . 5 Figure 9 .

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

Figure 10 .

 10 Figure 10. Deformed hemisphere computed by the full FE Newton-Raphson equilibrium iteration.

Figure 11 .Figure 12 .Figure 13 .Figure 14 .Figure 15 .Figure 16 .

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

Figure 17 .Figure 18 .

 1718 Figure 17. u y -contour of the pinched cornea obtained by a full FE model with Newton-Raphson iterations.

Figure 19 .

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

Table I .

 I Time (CPU seconds) employed for the solution of the three examples. was used) AMD Quad Opteron running at 2.2 GHz and with 16 Gb RAM, under Scientific Linux and running Matlab v. 2007.

	Example	FEM+Newton-Raphson	POD+ANM
	Beam	1.48	0.11
	Hemisphere	225.1	3.92
	Cornea	8.64E04	79.62
	no parallel computing


J Yvonnet

D Umbrello

F Chinesta

F Micari

A simple inverse procedure to determine heat flux on the tool in orthogonal cutting

Keywords: Heat flux, Split-tool, FEM, Inverse approach

The applications of numerical simulation to machining processes have been more and more increasing in the last decade: today, a quite effective predictive capability has been reached, at least as far as global cutting variables (for instance cutting forces) are concerned. On the other hand, the capability to predict local cutting variables (i.e. stresses acting on the tool, temperature distribution, residual stresses in the machined surface) has to be furtherly improved, as well as effective experimental procedures to validate numerical results have to be developed.

The aim of this paper is the proposition of an innovative approach, based on an simple inverse procedure, in order to identify both the heat flux flowing into the tool through the rake face and the heat transfer coefficient between the tool and the environment during a typical orthogonal cutting process. It is worth pointing out that the effective determination of such quantities is necessary in order to carry out a reliable prediction of the temperature distribution in the tool during the process.

The procedure is based on the integration of numerical simulations, inverse approach algorithms and experimental tests.

Introduction

In the last two decades, finite element techniques have achieved a very satisfactory level of development and reliability; in particular, the capability of the utilized codes to model cutting processes has supplied new powerful tools for machining processes analysis, although orthogonal cutting has been mainly taken into account, due the well known geometrical simplifications.

Nowadays, a relevant know-how about some relevant simulation issues can be considered as fully assessed: several topics, such as material behaviour under high strain, strain rate and temperature conditions, friction on the rake face of the tool and chip generation and development have been investigated by several researchers all over the world [START_REF] Merchant | Basic mechanics of the metal-cutting process[END_REF][START_REF] Shaw | Shear angle relationship in metal cutting[END_REF][START_REF] Ceretti | Basic Aspects and Modeling of Friction in Cutting[END_REF][START_REF] Ceretti | Prediction of Chip Morphology in Orthogonal Cutting by Means of a Customized Finite Element Code[END_REF][START_REF] Bacaria | Transient numerical models of discontinous chip formation based on damage effect[END_REF][START_REF] Bariani | Analysis of material behaviour at high strain rates for modelling machine processes[END_REF][START_REF] Zhao | Material constitutive laws at high strain rate for numerical simulations: testing and modeling in[END_REF][START_REF] Changeaux | High speed shear tests for the identification of the Johnos-Cook law[END_REF][START_REF] Tonshoff | Amor cutting of hardened steel[END_REF]. As well as a satisfactory prediction of global variables, for instance the cutting forces, is nowadays possible, even utilizing simple simulation tools.

On the other hand, several other issues related to machining processes require a deeper research effort, many of which are strictly related to the industrial application of machined parts: problems such as residual stresses prediction and control on the machined surfaces, simulation of complex fully 3D operations or, finally, cutting temperatures and tool wear prediction during machining operations are just few examples of new emerging needs.

During machining processes cutting tools are subjected to mechanical, thermal and tribo-chemical loads. In order to improve the scientific knowledge on machining technology and, moreover, in order to develop new powerful and reliable cutting tools, the simple knowledge of the global mechanical and thermal loads acting on the tool is not sufficient. In turn, a sound prediction of pressure and heat flux distributions on the rake face is necessary, i.e. the analysis must be developed at a local level.

In the past, many attempts to determine the local distribution of normal and shear stresses on the tool surface were carried out using photoelastic materials for either the tool or the workpiece, cutting tool materials with a well defined limit of elasticity and tools characterized by a reduced contact length including split-tools [START_REF] Coker | Report on the action of cutting tools[END_REF][START_REF] Dietrich | Untersuchungen des Spannungszustandes vor einer Drehstahlschneide im polarisierten Licht[END_REF][START_REF] Amini | Photoelastic analysis of stresses and force in steady cutting[END_REF][13][14].

More recently, some researchers proposed different approaches to quantify the stress distribution at the toolchip interface [START_REF] Kato | Stress Distribution at the Interface between Tool and Chip in Machining[END_REF][16][START_REF] Childs | On the stress distribution between the chip and the tool during metal turning[END_REF], mainly based on the use of split-tool dynamometers. Several experiments based on the same techniques were also performed by Buryta et al. [START_REF] Buryta | Stress distribution on the rake face during orthogonal machining[END_REF] and by Tonshoff et al. [START_REF] Tonshoff | Analysis of the Effect of Thermal and Mechanical Stress on Tool Wear in Continuous Cutting Report of the Institute for Production Engineering and Machine Tool[END_REF].

Unfortunately, the use of experimental approaches, such as the split dynamometer, is very difficult and time consuming; on the other hand the development of a reliable numerical simulation, able to determine the stress distributions, would offer relevant advantages. A finite element analysis of the normal and the shear stresses distribution on the rake face in orthogonal cutting was proposed by McClain et al. [START_REF] Mcclain | A numeric investigation of the rake face stress distribution in orthogonal machining[END_REF]. Actually, the predicted normal stresses on the rake face only qualitatively agreed with the experimental distribution determined by other researchers.

Recently, a new approach, based on the integration of numerical and experimental techniques was proposed by Filice et al. [START_REF] Filice | Pressure distribution on the tool in cutting: Prediction and measurement[END_REF]. They utilized a set of properly designed tools, characterized by proper cuts at well defined distances from the tool tip, to measure the pressure distribution in orthogonal cutting; as well they carried out extensive numerical simulations, achieving a satisfactory agreement between the experimental and the predicted results.

The basic idea of this paper is to apply the same philosophy of the previous paper (Filice et al. [START_REF] Filice | Pressure distribution on the tool in cutting: Prediction and measurement[END_REF]), in order to identify the heat flux distribution flowing into the tool through the rake face by means of an inverse approach. The method here addressed is based in fact on the partitioning concept utilized by Filice et al. [START_REF] Filice | Pressure distribution on the tool in cutting: Prediction and measurement[END_REF] to measure the pressure distribution on the tool in orthogonal cutting.

A set of different tools is utilized: each of them is characterized by a tight slot cut at a proper distance from the tool tip. In this way, the tool is partitioned into two regions fully independent from the thermal point of view: it means that, during machining, a portion of the total heat flux flows on each of the regions, determining in each of them a temperature distribution which only depends on the portion itself. More in detail, the temperature distribution in each region is not affected by the heat flux portion flowing into the other region of the tool.

In this way, some temperature values in the tool region closer to the tool tip are accurately measured by means of a thermocouple directly mounted on the tool flank. Then an inverse procedure is carried out to determine both the convection coefficient on the external surface of the cutting tool and the average heat flux flowing into the region of the tool closer to the tool tip. Repeating the same approach as above for each tool, i.e. at the varying of the distance of the tight slot from the tool tip, the whole heat flux distribution can be reconstructed.

In the following paragraphs of the paper, the subsequent steps of the research will be presented, namely the preliminary investigations to verify if the partitioning concept holds for thermal phenomena, the design of the utilized tools, the experimental tests and finally the inverse approach.

Tool design

As mentioned the inverse procedure is based on temperature data supplied by a thermocouple mounted on the flank face of the tool. Furthermore, the tight slot cut into the tool must ensure the effectiveness of the thermal partitioning described above. Therefore, a proper design of the instrumented tool is necessary.

Two constrains are involved in this design: (i) the thermal field in the former region of the tool, i.e. the one closer to the tool tip, must be independent on the heat flux flowing into the remaining portion of the chip-tool contact surface, at least for a cutting time t max long enough for the application of the inverse procedure; (ii) the position of the thermocouple along the flank of the tool must permit to avoid local effects, i.e. an unique relationship should exist between the heat flux flowing into the former region of the tool and the temperature history measured by a thermocouple. A calibrated plate made of a material with a very low thermal conductivity was introduced in the slot. The material was Silica-Aluminous based ceramic fibres, characterized by service temperature up to 1260 8C, high stability, low thermal conductivity, low thermal storage capacity, excellent resistance to thermal shocks. The application of the calibrated plate was necessary both in order to ensure the thermal effects separation and to avoid chip introduction into the slot.

Nomenclature

a tip-slot distance, mm A thermocouple position E error C heat capacity, J/(m 3 K) h convection coefficient, W/(m 2 K) k thermal conductivity W/(m K) K tangent matrix in the Newton strategy Q heat flux, J/(m 2 s) R residual
Table 1 reports the thermal parameters for the workpiece (Aluminum Alloy 6060-T5), the tool (HSS) and the calibrated plate (Silica-Alumina), as used in the FEM analysis.

In addition, the same authors have shown in a recent work [START_REF] Yvonnet | An inverse procedure to determine heat flux on the tool in orthogonal cutting[END_REF], the validation of the proposed technique taking into account the cutting temperature distributions and the thermal separation effects by a comparison between the experimental evidences and the numerical results.

Design of the optimum slot geometry

A set of numerical simulations at the varying of the slot geometry was carried out in order to determine the above mentioned t max value, i.e. the time length during which the thermal field in the former region of the tool, is independent on the heat flux flowing into the remaining portion of the chip-tool contact surface. The optimum tight slot design is reported in Fig. 1. It is characterized by a depth of 8 mm and an inclination angle of 458.

Such design permits to obtain a cutting time t max equal to 17 s as confirmed by the numerical results reported in Figs. 2 and3. Fig. 2 shows that applying a thermal flux only on the tool portion far away from the tool tip, after 17 s only the 10% of such flux flows into the former portion of the tool.

In turn, Fig. 3 shows that applying a thermal flux only in the portion of the tool closer to the tool closer to the tool tip, heat remains confined in such portion for more than 17 s.

Design of the optimum thermocouple position

The thermal model is assumed isotropic and both the heat capacity and the thermal conductivity are assumed known. If, as mentioned, the cutting tool is characterized by a slot filled with an non-conductive material (whose conductivity k isol can be neglected, i.e. k isol z0), and the intensity of the convection effects on the tool surface are much smaller that the heat conduction mechanism, the following thermal result holds: it exists a position x (see Fig. 4) that must be determined such that the temperature evolution at this position within the time interval [0,t max ] depends on the average heat flux

Q Z Ð a 0 QðyÞdy a (1)
and in consequence it does not depend on the local details of Q(y).

If the thermocouple is too near to the tool tip, the measured temperature evolution will depend on the local details of the incoming heat flux Q(y) and not only on its average value Q. On the other hand, if the thermocouple is too far from the tool tip, the measured temperatures will not depend only on the average heat flux incoming in the former portion of the tool, but also they will be influenced by the heat flux flowing through the second portion of the tool.

According to the above reported reasons, it is possible to determine the best position x within the t max time interval minimizing the error E(x) defined by EðxÞ Z Tðx; QÞKTðx; QðyÞÞ Tðx; QðyÞÞ

being Tðx; QÞ the temperature at position x due to average heat flux Q and T(x,Q(y)) the one due to the actual heat distribution.

For the sake of simplicity, the error minimization was carried out by assuming the heat flux distribution incoming on the former portion of the tool as reported in Fig. 5. If Q a and Q b are assumed equal to 1.5 Q and 0.5 Q, respectively (their average value being Q), the validity of the results will be guaranteed for variations of the thermal fluxes with respect to the average heat flux lower than 50%.

From the calculated error for different x values within the [0,t max ] time interval we can conclude that for x greater than 2a the error is negligible. According to these results and taking into account several constrains due to the dimensions of the thermocouple (1 mm diameter), the optimum distance from the tool tip was detected as equal to 2.25 mm, as reported in Fig. 6.

Experimental tests

The aim of the experimental tests was to measure the temperature profile related to each tool by using a Chromel w /Alumel w thermocouple (type K) mounted on the former portion of the tool (Fig. 7).

Four different tools were utilized (Fig. 8): three of them were manufactured by electrical discharge machine (EDM) cutting the tight slot with different distances between the tool tip and the slot, namely 0.35, 0.5 and 0.60 mm; the last one was instrumented only with the thermocouple in order to measure the whole average heat flux flowing into the tool through the rake face.

The thermocouple with a diameter of 1 mm was embedded in the tool applying a force of about 50 N in order to ensure a sufficient adhesion [START_REF] Astakhov | Metal Cutting Mechanics[END_REF]. Both the contact length and the thermocouple positions were measured utilizing a proper optical microscope. The former by depositing on the rake face a special paint and observing the rake face, the latter by observing a posteriori the tool section (along the hole axes) obtained by EDM.

The cutting speed and the feed rate was fixed to 125 m/min and 0.4 mm/rev, respectively, to guarantee a tool-chip contact surface long enough to include all the slots (the measured contact length was about 2 mm).

The experimental equipment was completed by a Personal Computer based on Intel Pentium IV technology equipped with a data acquisition board (AT-MIO-16E-10 by National Instruments). The software used for the acquisition and the analysis of the data supplied by the thermocouple measuring devices was SG411 by Sint Technology. Fig. 9 shows the average measured temperature histories during the experimental tests. All the tests were executed in the same environment conditions with almost three repetitions for each of them.

Depending on the slot position, i.e. depending on the position of the non-conductive plate, different fractions of the total heat flux flow into the former portion of the tool, determining different temperature curves; the more the non-conductive plane is far away from the tool tip, the higher is the heat amount flowing into the tool and the higher is the measured temperature. Of course the difference among the curves reduces after a certain cutting time that is when the partitioning effect disappears.

Fig. 9 also demonstrates that the partitioning effects is confirmed until 15 s, while for a longer time the different tools reach the same temperature value of the tool without the slot, since the heat flux flowing into the second portion of the tool affects the whole tool.

Heat flux identification using an inverse approach

Using the divided tool technique described previously, the lower part of the tool was thermically isolated. The evolution of the temperature in time during cutting was measured at point A (see Fig. 10) using a thermocouple, which was positioned according to the discussion in Section 3 (Fig. 6).

For each tool, i an average flux Qi on the isolated part of the rake face was determined through the inverse procedure described below, which in fact is associated with the real flux distribution Q(y) on the considered part by the

expression Qi Z Ð a i 0 Q i ðyÞdy.
Given different distances of the isolated band a i and the associated average fluxes Qi , the real heat flux distribution Q(y) can be approximated by a piecewise constant function (see Fig. 11) by means of the following scheme

Qi Z Q 1 ; Q2 Z Q 1 a 1 C Q 2 ða 2 Ka 1 Þ a 2 ; Qi Z Q 1 a 1 C P i jZ2 Q j ða j Ka jK1 Þ a i (3)
from which For a given a i value and the corresponding experimental temperature history T exp (t) at point A, the average flux Qi was determined using the finite element model depicted in Fig. 12.

Q 1 Z Q1 ; Q 2 Z Q2 a 2 KQ 1 a 1 a 2 Ka 1 ; Q i Z Qi a i KQ 1 a 1 K P iK1 jZ2 Q j ða j Ka jK1 Þ a i Ka iK1 (4)
Identification of both unknowns, i.e. the average heat flux Qi and the convection coefficient h (see Fig.

i ; h k g, f Qk i C 3 Q ; h k g, f Qk i K3 Q ; h k g, f Qk i ; h k C 3 h g
vE k v Qi Z E k 3Q KE k K3Q 23 Q and vE k vh Z E k 3h KE k K3h 23 h
If kRk!TOL then f Qi Z Qk i ; hZ h k g; else continue

Construct the tangent matrix. For this purpose compute four more temperature histories at position A are calculated, namely, the ones associated to the following conditions:

f Qk i C 3 Q ; h k C 3 h g, f Qk i K3 Q ; h k C 3 h g, f Qk i C 3 Q ; h k K3 h g and f Qk i K3 Q ; h k K3 h g. They permit to compute the associated errors, denoted by E k 3Q;3h , E k K3Q;3h E k 3Q;
K3h and E k K3Q; K3h , and therefore, the tangent matrix

K Z v 2 E k v Q2 i v 2 E k v Qi vh v 2 E k v Qi vh v 2 E k vh 2 0 B B B B @ 1 C C C C A
where the different components of the tangent matrix are the following

v 2 E k v Q2 i Z E k 3Q K2E k C E k K3Q 3 2 Q ; v 2 E k vh 2 Z E k 3h K2E k C E k K3h 3 2 h ; v 2 E k v Qi vh Z E k 3Q;3h KE k K3Q;3h KE k 3Q; K3h C E k K3Q; K3h 43 Q 3 h
i ; Dh k g QkC1 i Z Qk i C D Qk i ; h kC1 Z h k C Dh k
Return to step 2.

Results And discussion

The described approach was applied to the complete set of tools reported in Section 3. For each of them, the experimental temperature histories were shown in Fig. 9.

Fig. 13 reports the calculated temperature distribution for one of the cases taken into account, namely the tool characterized by a distance of the slot equal to 0.65 mm. In this case, only four Newton iterations were necessary until convergence.

Finally, Fig. 14 shows the obtained piecewise constant heat flux distribution.

Such distribution is quite consistent with the already available scientific knowledge. In fact, while for the pressure distribution, there is a monotonic exponential distribution with the maximum at the tool tip, for the heat flux distribution the maximum is located close to the crater region, i.e. where the maximum energy consumption due to friction occurs. The convection coefficient, h, was found equal to 0.5 W/m 2 K.

Conclusions

According to results presented in this paper, the following conclusions can be addressed: an innovative approach based on the integration of experimental tests, numerical simulations and inverse procedures was developed in order to determine the heat flux distribution flowing into the tool in orthogonal cutting; the basic idea behind this research is the application of the so called partitioning concept: in this way, utilizing a set of properly designed tools, it was possible to effectively calculate the heat flux distribution flowing into the tool; such result represents the most relevant outcome of the proposed approach: according to the available knowledge on the subject area, only an average heat flux value is generally calculated and heat flux distribution is assumed constant on the contact area; on the contrary a detailed prediction of the heat flux distribution was obtained, therefore, achieving a deeper insight on the physical phenomena occurring during machining; the utilized approach required the development of a set of instrumented split-tools, characterized by embedded thermo-couples and properly designed tight slots filled with non-conductive material; as well simple and well suitable inverse procedures were implemented in order to achieve the best matching between numerical and experimental results.

Fig. 2 .

 2 Fig. 2. After 17 s only the 10% of the applied heat flux flows into the former portion of the tool.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Heat remains confined in the former portion of the tool for more than 17 s.

Fig. 5 .

 5 Fig. 5. Used model for error minimization.

Fig. 6 .Fig. 7 .

 67 Fig. 6. Sketch of the tool utilized in the tests.

  12) was possible by sing an iterative Newton-Raphson procedure to minimize the error defined by difference between the calculated and the experimental temperature history in point A: Eð Qi ; hÞ Z ð t max 0 ½T num ðA; tÞKT exp ðtÞ 2 dt The algorithm was the following: Assume the initial values Q0 i and h 0 , or consider known their values at iteration k: Qk i and h k Compute the residual: RZ fvE k =v Qi ; vE k =vhg. For this purpose, five temperature histories at position A associated with the following conditions: f Qk

Fig. 8 .Fig. 9 .

 89 Fig. 8. The tools utilized in the experimental tests.

Fig. 10 .

 10 Fig. 10. Schematic view of the experimental setup.

Fig. 11 .Fig. 12 .

 1112 Fig. 11. Heat flux approximation.

Table 1

 1

	Thermal parameters
		Heat capacity	Thermal conduc-	Mass density
		J/(m 3 K)	tivity W/(m K)	(Kg/m 3 )
	Tool	3.79!10 6	50	7870
	Workpiece	2.42!10 6	160	2700
	Calibrated plate	2.83!10 6	0.09	2700

Fig. 1. The optimum slot design.

Acknowledgements

The authors wish to thank Mr Franco Pulice of the Department of Mechanical Engineering of the University of Calabria for his continuous and qualified support in the preparation and development of the experimental tests.

Authors also thank the reviewers for their inputs and comments.