Siamak Niroomandi

Icíar Alfaro

David Gonzalez

Elías Cueto

Francisco Chinesta

D González

Real Time Simulation of Surgery by Reduced Order Modelling and X-FEM Techniques

Keywords: real time, surgery, model order reduction, X-FEM

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Development of surgery simulators for the risk-free training of surgeons has been a very active field of research in the last decades. Of course, this development is a very challenging task because of the very different complexities arising in this kind of simulators [START_REF] Delingette | Hepatic surgery simulation[END_REF]. One of the sources of complexity is due to the highly non-linear behavior of soft living tissues, which are frequently modeled under the fiber-reinforced hyperelasticity framework [START_REF] Alastrue | Biomechanical modeling of refractive corneal surgery[END_REF]. Other source of complexity comes from the highly restrictive feedback rates imposed by the simulators (25 Hz for visual feedback and some 500 Hz if we want to add haptic feedback to the system). The third source of complexity comes from the multi-physic nature of the phenomena occurring in the actual surgery procedure: non-linear elasticity, contact, cutting, temperature, and so on.

Such simulators should provide a physically more or less accurate response such that, with the use of haptic devices, a realistic feedback is transmitted to the surgeon in terms of both visual feedback and force feedback. By 'accurate response', we mean that an advanced user should not encounter 'unphysical' sensations when handling the simulator. We definitely do not pursue an accurate solution in engineering terms. Following [START_REF] Bro-Nielsen | Real-time volumetric deformable models for surgery simulation using finite elements and condensation[END_REF], '. . . the model may be physically correct if it looks right'.

Ayache and co-workers [START_REF] Delingette | Hepatic surgery simulation[END_REF] defined the three generations of surgical simulators as those able to, respectively, reproduce accurately the anatomy (geometry), the physics (constitutive equations of soft tissue, temperature), and finally, physiology (blood flow, breathing, etc.). Undoubtedly, no third-generation simulation has ever been developed, and only some rigorous attempts have been made at the second-generation level (see [START_REF] Taylor | On modelling of anisotropic viscoelasticity for soft tissue simulation: numerical solution and GPU execution[END_REF][START_REF] Courtecuisse | Gpu-based real-time soft tissue deformation with cutting and haptic feedback[END_REF], among others). Most of the existing simulators can be classified into the first-generation category, even if they provide haptic feedback, because they only consider linear elastic response, for instance [START_REF] Bro-Nielsen | Real-time volumetric deformable models for surgery simulation using finite elements and condensation[END_REF][START_REF] Wang | Virtual reality simulation of surgery with haptic feedback based on the boundary element method[END_REF][START_REF] Cotin | Real-time elastic deformations of soft tissues for surgery simulation[END_REF][START_REF] Basdogan | Real-time simulation dynamically deformable finite element models using modal analysis and spectral Lanczos decomposition methods[END_REF]. Other are based upon spring and mass systems, which do not even reproduce the equations of linear elasticity (see [START_REF] Delingette | Soft tissue modeling for surgery simulation[END_REF] and references therein).

Incorporating non-linear tissue constitutive equations into real-time frameworks is challenging in itself. Among the various techniques developed to reach this objective, we can distinguish techniques on the basis of the use of explicit finite elements and GPUs (see, for instance [START_REF] Taylor | On modelling of anisotropic viscoelasticity for soft tissue simulation: numerical solution and GPU execution[END_REF][START_REF] Courtecuisse | Gpu-based real-time soft tissue deformation with cutting and haptic feedback[END_REF]). Limitations of explicit finite elements with respect to stability of the result for large time steps are well known, however. Others have preferred to employ reduced-order modeling through proper orthogonal decomposition (POD), also known as principal component analysis (PCA) [START_REF] Barbič | Real-time subspace integration for St. Venant-Kirchhoff deformable models[END_REF]. In either of these approaches, non-linear must be understood as 'geometrically non-linear', as the aforementioned approaches include simple co-rotational formulations or Saint Venant-Kirchhoff constitutive equations. References [START_REF] Niroomandi | Real-time deformable models of non-linear tissues by model reduction techniques[END_REF][START_REF] Niroomandi | Model order reduction for hyperelastic materials[END_REF][START_REF] Niroomandi | Accounting for large deformations in real-time simulations of soft tissues based on reduced-order models[END_REF] seem to be among the earliest works that incorporate material (not only geometrical) non-linearity into the real-time framework.

The third source of difficulties, as mentioned earlier, comes from the multi-physic nature of the problem. Although the thermal dependence of the problem is often-if not always-neglected, contact detection and cutting simulation are of utmost importance for a convincing result in terms of both visual and haptic perception.

Numerical simulation of contact mechanics is a well-known field in the computational mechanics community (see, for instance the recent book of Wriggers [START_REF] Wriggers | Computational Contact Mechanics[END_REF] to acquire an overall impression of the difficulty of the topic). An accurate simulation of the process of contact between surgical tools and organs, and between organs themselves, seems to be at this moment out of reach under real-time requirements. Some simplified algorithms, however, provide very realistic results at high feedback rates (see for instance [START_REF] Barbič | Time-critical distributed contact for 6-DoF haptic rendering of adaptively sampled reduced deformable models[END_REF]). This algorithm supports distributed contact detection between a complex object (the surgical tool, for instance) and a deformable body reduced model, at 1 kHz rates.

As mentioned before, cutting simulation is another important source of troubles for real-time modeling of surgery. This is so as it is necessary to modify the geometry and/or the topology of the domain and its associated mesh, and this needs to be carried out without penalizing the computation times of the integration of the equations of motion. A vast corps of literature has been devoted to this end. See, for instance, [START_REF] Courtecuisse | Gpu-based real-time soft tissue deformation with cutting and haptic feedback[END_REF][START_REF] Zhang | An advanced hybrid cutting method with an improved state machine for surgical simulation[END_REF][START_REF] Lee | Topology modification for surgical simulation using precomputed finite element models based on linear elasticity[END_REF][START_REF] Meier | Real-time deformable models for surgery simulation: a survey[END_REF], to name a few. All these approaches share the same spirit. All of these works propose more or less sophisticated algorithms that allow to remove finite elements from an existing mesh to create a cut or even to remove a whole part of the organ, if it is being ablated. None of them seem to be aware of the existence of extended finite element method (X-FEM) (see for instance [START_REF] Sukumar | Extended finite element method for three-dimensional crack modelling[END_REF] and the subsequent enormous list of papers derived from the seminal work of Moës and colleagues).

Therefore, after a more or less exhaustive review of the existing literature, a conclusion can be drawn. Techniques based upon model order reduction constitute nowadays an appealing choice for an efficient, real-time, integration of the equations of motion [START_REF] Krysl | Dimensional model reduction in non-linear finite element dynamics of solids and structures[END_REF]. They seem to be unique to efficiently handle material non-linearities such as those arising from state-of-the-art models for living soft tissues [START_REF] Niroomandi | Model order reduction for hyperelastic materials[END_REF][START_REF] Niroomandi | Accounting for large deformations in real-time simulations of soft tissues based on reduced-order models[END_REF]. And, finally, they seem to be specially well suited for efficient contact detection algorithms in real time [START_REF] Barbič | Time-critical distributed contact for 6-DoF haptic rendering of adaptively sampled reduced deformable models[END_REF].

Reduced-order modeling employs global (Ritz) shape functions to construct the approximation spaces in a Galerkin framework. These shape functions are usually constructed to be optimal under certain requirements. For instance, in POD techniques, these functions are extracted after a statistical treatment of the results of problems similar to the one at hand [START_REF] Lorenz | Empirical Orthogonal Functions and Statistical Weather Prediction[END_REF][START_REF] Park | The use of the Karhunen-Loève decomposition for the modeling of distributed parameter systems[END_REF][START_REF] Brooks | A Karhunen-Loeve Galerkin Technique with Shock Fitting for Optimization of a Blunt Body Geometry[END_REF][START_REF] Karhunen | Uber lineare methoden in der wahrscheinlichkeitsrechnung[END_REF]. These techniques are known in many branches of science and technology under a wide variety of names, including the aforementioned POD, PCA, or Karhunen-Loève decomposition. Other techniques are even able to construct these shape functions without any prior knowledge on the results of similar problems (see [START_REF] Ryckelynck | A priori hyperreduction method: an adaptive approach[END_REF][START_REF] Ryckelynck | On the a priori model reduction: overview and recent developments[END_REF][START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids[END_REF][START_REF] Ammar | 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[END_REF][START_REF] Chinesta | Recent advances in the use of the proper generalized decomposition for solving multidimensional models[END_REF]).

However, the use of global, Ritz, shape functions imposes additional limitations. For instance, the use of X-FEM techniques is not straightforward in this context. In this paper, a new method is presented, which combines the features of existing methods based upon model order reduction [START_REF] Niroomandi | Model order reduction for hyperelastic materials[END_REF][START_REF] Niroomandi | Accounting for large deformations in real-time simulations of soft tissues based on reduced-order models[END_REF] and the ease of creating cuts and discontinuities without remeshing of the X-FEM technique. As will be noticed, the proposed technique has its origins in the s-FEM techniques by Fish [START_REF] Fish | The s-version of the finite-element method[END_REF] and the multiscale FEM by Rank [START_REF] Rank | A multiscale finite-element method[END_REF].

The paper is organized as follows. In Section 2, the basics of POD techniques is reviewed in order to show how it can deal with non-linear constitutive equations under real-time requirements. Then, in Section 3, a new technique is developed, which can efficiently create cuts and discontinuities in the reduced model without remeshing and with minor modifications in the problem's stiffness matrix. Finally, in Section 4, an analysis of its performance is performed in an application to refractive surgery of the cornea.

MODEL ORDER REDUCTION BY PROPER ORTHOGONAL DECOMPOSITION METHODS

The Karhunen-Loève decomposition

The technique here considered can be seen as a particular instance of a posteriori model order reduction techniques. This means that the process to achieve real-time performance is composed by a sequence of off-line simulation, storage of the results, and subsequent analysis, followed by an on-line process of real-time simulation that employs the Ritz functions obtained in the off-line phase of the method. Following Cotin and Bro-Nielsen, '. . . We do not care about the time taken for one-time pre-calculation such as setting up equations, inverting matrices, etc.' [START_REF] Bro-Nielsen | Real-time volumetric deformable models for surgery simulation using finite elements and condensation[END_REF]. Therefore, we assume that the evolution of a certain (scalar or vectorial) field u.x, t/, governed by a PDE, is known, typically from existing previous simulations on problems similar to that we are interested in. In practical applications, this field is expressed in a discrete form, that is, it is known at the nodes of a spatial finite element mesh and for some time steps of existing simulations u.x i , t n / Á u n i . The same can be written by introducing a time discretization u n .x/ Á u.x, t D nt /I 8n 2 OE1, , 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/ 8n. This is equivalent to obtaining a function .x/ that maximizes ˛defined by

˛D P nDP nD1 h P i DN i D1 .x i /u n .x i / i 2 P i DN i D1 . .x i // 2 (1)
This maximization (ı˛D 0) leads to

nDP X nD1 2 4 iDN X iD1 Q .x i /u n .x i / ! 0 @ j DN X j D1 .x j /u n .x j / 1 A 3 5 D ˛iDN X i D1 Q .x i / .x i /I 8 Q (2) which is equivalent to i DN X i D1 8 < : j DN X j D1 " nDP X nD1 u n .x i /u n .x j / .x j / # Q .x i / 9 = ; D ˛iDN X i D1 Q .x i / .x i /I 8 Q (3)
Defining the vector such that its i th component is .x i /, Equation (3) takes the following matrix form

Q T c D ˛Q T I 8 Q ) c D ˛ (4)
where the two-point correlation matrix is given by

c ij D nDP X nD1 u n .x i /u n .x j / , c D nDP X nD1 u n .u n / T (5)
which is symmetric and positive definite. If we define the matrix Q containing the discrete field history

Q D 0 B B B @ 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 1 C C C A (6)
it is straightforward to verify that the matrix c in Equation ( 4) results

c D Q Q T (7)
where the diagonal components are given by

c i i D .Q Q T / i i D j DP X j D1 .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 modeling

As mentioned before, the main ingredient of the POD method is to perform off-line some initial simulations on the complete model(s) that allow to obtain u.

x i , t n / Á u n i , 8i 2 OE1, , N , 8n 2 OE1,
, P , and from these, the r eigenvectors related to the r-highest

eigenvalues k D k .x i /, 8i 2 OE1,
, N , 8k 2 OE1, , r (with r N ). From these simulations, the r eigenfunctions can be successfully employed for approximating the solution of a problem slightly different than the one that has served to define u.x i , t n /. For this purpose, it is necessary to define the matrix

A A D 0 B B B @ 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 / 1 C C C A . ( 9
)
The typical structure of the discretized form of the (non-linear) equations of elasticity is of the type

K U D F . ( 10
)
Obviously, in the case of time-dependent problems, F contains the contribution of the solution at the previous converged time step. Assuming that the unknown vector contains the nodal (usually, displacement) degrees of freedom, it can be expressed as

U D i Dr X i D1 i i D A , (11)
which gives, for the complete system of equations,

K U D F ) K A D F . ( 12
)
Multiplying both terms by A T gives

A T K A D A T F , ( 13
)
which proves that the final system of equations is of low order, that is, the dimensions of A T K A are r r, with r N , and the dimensions of both and A T F are r 1.

From the paragraphs above, we recover the initial idea exposed in Section 1. 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 FEM, which employs general purpose, piecewise polynomial shape functions to approximate the solution in a Galerkin framework.

Limitations of classical model reduction techniques

Although the model order reduction technique reviewed before gives excellent results in some, very useful, situations (see, for instance, [START_REF] Niroomandi | Real-time deformable models of non-linear tissues by model reduction techniques[END_REF][START_REF] Krysl | Dimensional model reduction in non-linear finite element dynamics of solids and structures[END_REF]) and has given very promising results at real-time feedback rates (see [START_REF] Barbič | Real-time subspace integration for St. Venant-Kirchhoff deformable models[END_REF]), it still has some limitations. The main limitation deals with the severe limitations imposed by real-time constraints in haptic environments (set to some 0.5-1 kHz). This means that in practical situations, it is not possible to update the stiffness matrix K of Equation [START_REF] Niroomandi | Accounting for large deformations in real-time simulations of soft tissues based on reduced-order models[END_REF] in the course of the simulation. For strongly non-linear tissue constitutive equations, this can give errors below the 20% of the solution obtained by classical non-linear finite elements [START_REF] Niroomandi | Real-time deformable models of non-linear tissues by model reduction techniques[END_REF], still valid in some situations.

In some previous works of the authors, see [START_REF] Niroomandi | Model order reduction for hyperelastic materials[END_REF][START_REF] Niroomandi | Accounting for large deformations in real-time simulations of soft tissues based on reduced-order models[END_REF], a method has been developed and validated that allows for a fully non-linear simulation (in some neighborhood of previous, equilibrated, solutions) without the need of tangent stiffness matrix updating. The method is based on the combination of model order reduction techniques and asymptotic expansions of the solution (the so-called asymptotic numerical method [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][START_REF] Cao | An improved iterative method for large strain viscoplastic problems[END_REF][START_REF] Cochelin | The asymptotic numerical method: an efficient perturbation technique for nonlinear structural mechanics[END_REF]).

The second, very important, limitation of classical model reduction techniques (also of more recent versions such as those in [START_REF] Niroomandi | Real-time deformable models of non-linear tissues by model reduction techniques[END_REF][START_REF] Niroomandi | Model order reduction for hyperelastic materials[END_REF]) is that the employ of Ritz functions very much complicates the issue of simulating surgical cutting. The appearance of a cut implies a change of topology in the geometry of the domain and also an important change in the stiffness matrix, element connectivities, and so on. X-FEM techniques [START_REF] Sukumar | Extended finite element method for three-dimensional crack modelling[END_REF] very much limit these inconveniences for standard finite element simulations, with only minimal changes in the number of degrees of freedom and element connectivity, but these are not readily applicable to reduced-order models employing Ritz basis functions.

In the next section, we briefly review the X-FEM in the framework of standard finite elements, prior to its introduction in conjunction with reduced basis.

INTRODUCING CUTS IN THE REDUCED MODEL BY X-FEM TECHNIQUES

A brief review of the X-FEM

In order to make this work self-contained, a brief review of the well-known X-FEM is included here. The interested reader is referred to [START_REF] Sukumar | Extended finite element method for three-dimensional crack modelling[END_REF], for instance, among other classical references in the field.

The basic ingredient of X-FEM methods is to consider a cut or crack d as a discontinuity in the displacement field. Therefore, by simply enriching those nodes whose shape function's support intersects, the crack with a discontinuous functions will suffice to obtain a conforming discretization without the need of remeshing. Thus, the new approximation of the displacement will be, in the context of X-FEM, for a single crack or cut,

u h .x/ D X i2I u i N i .x/ C X j 2J b j N j .x/H.x/, (14)
where J D fj 2 I W ! j \D ¤ ;g, I represents the set of all nodes in the FE mesh, N i .x/ represents the i th node shape function, whose support is ! i , evaluated at x, u i represents node i's displacement vector, and H.x/ is discontinuous across the crack d , whose geometry is represented by D.

Finally, b j represents a new set of enriched degrees of freedom that control the magnitude of the displacement discontinuity across the crack. Classical works of X-FEM in the context of fracture mechanics such as in [START_REF] Sukumar | Extended finite element method for three-dimensional crack modelling[END_REF] still add some more degrees of freedom in order to reproduce the linear elastic fracture mechanics asymptotic solution at the crack front. These are not considered here for several reasons. Firstly, for living soft tissues, the form of the asymptotic solution (if any) at the crack front is in general not known. In any case, the gain in accuracy would be, in general, negligible, if compared with all the simplifying assumptions taken up to now in order to reach a real-time performance.

Coupling reduced models and X-FEM descriptions of surgical cuts

The approach here suggested is composed by the coupling between the (Ritz) reduced basis and a superimposed patch of finite elements where the force exerted by the scalpel exceeds some limit value, thus appearing a cut in the organ. This approach can be seen as a generalization of some existing techniques in the field of finite elements, notably the so-called s-version of the FEM by Fish [START_REF] Fish | The s-version of the finite-element method[END_REF] and the multiscale FEM by Rank [START_REF] Rank | A multiscale finite-element method[END_REF], among others [START_REF] Ammar | Coupling finite elements and reduced approximation bases[END_REF][START_REF] Ammar | Coupling finite elements and proper generalized decompositions[END_REF]. The basic of the method is easily understood from Figure 1.

The domain is discretized by finite elements, and a reduced-order model is constructed on top of it that employs Ritz basis functions, as explained before. Once the force exerted by the scalpel reaches a prescribed threshold (see the sections that follow), a cut is supposed to appear at the contact location. Around that region, a superimposed patch L of finite elements is placed around the scalpel tip. In L , the displacement field is approximated as

u D u 0 C u L in L ( 15
)
where, to guarantee the compatibility of the displacement field, we enforce

u L D 0 on 0L . ( 16
)
The reduced-order model approximates the displacement everywhere, in the spirit of Equation ( 11),

u 0 .x/ D i Dr X i D1 i i .x/, (17)
whereas in L a traditional FEM mesh, including discontinuous enrichment, is added:

u L .x/ D X j 2 L d j N j .x/ C X k2J b k N k .x/H.x/. ( 18
)
Figure 1. Basics of the method for reduced basis-X-FEM coupling.

The discrete equilibrium equations can be obtained after the weak form of the problem, namely,

find u 2 U D fu W u 2 H 1 , u D u on u g such that Z " W d D Z t u td C Z u bd , 8u 2 U 0 D fu W u 2 H 1 0 , u D 0 on u g, (19)
where, as usual, u represents the Dirichlet (essential) part of the boundary of the domain, u represents the value of the prescribed displacement at that location, t is the prescribed vector of traction at the boundary, and finally, t is the Neumann (natural) part of the boundary . Admissible variations of the displacement are calculated as

u D u 0 C u L . ( 20
)
By substituting the displacement field approximations ( 17) and ( 18) into the weak form [START_REF] Sukumar | Extended finite element method for three-dimensional crack modelling[END_REF], and after invoking the arbitrariness of admissible displacements u , we arrive at a discrete system of equations of the type 0

@ K K d K b K d K dd K db K b K bd K bb 1 A 0 @ d b 1 A D 0 @ f f d f b 1 A , ( 21
)
where the stiffness matrices K ˛ˇare given by

K ˛Ǐ J D Z i .B Į / T DB J d (22)
where matrix B Į , ˛D , d , b represent the standard shape function derivative matrix for any of the three types of approximation functions here considered, namely, the Ritz functions .x/, the standard finite element shape functions N.x/, or the discontinuous enrichment shape functions, N.x/H.x/. i represents the domain of integration, either or L , respectively. Finally, D represents the consistently linearized constitutive tensor. As mentioned before, because of the severe restrictions placed by real-time requirements in haptic environments, it is not possible, in general, to update this tangent stiffness tensor. In this work, for the sake of simplicity in the exposition of the method, an approach similar to that in [START_REF] Niroomandi | Real-time deformable models of non-linear tissues by model reduction techniques[END_REF] is employed, in which no updating is accounted for. This leads, of course, to higher errors in the results, still acceptable in virtual surgery environments.

If more accurate results are needed, an alternative formulation can be used [START_REF] Niroomandi | Accounting for large deformations in real-time simulations of soft tissues based on reduced-order models[END_REF], which employs asymptotic expansions of the displacement field and, without any updating, allows to closely follow complex non-linear force-displacement paths. The extension of the technique developed in some of our former works [START_REF] Niroomandi | Accounting for large deformations in real-time simulations of soft tissues based on reduced-order models[END_REF] to this framework is a work in progress.

In order to evaluate the previous integrals leading to the stiffness matrix of the problem, no special integration procedure is employed, in sharp contrast with the original works of Fish [START_REF] Fish | The s-version of the finite-element method[END_REF], as we assume that the superimposed mesh conforms with the existing one. It is necessary, however, to perform some form of tailored numerical integration in those elements enriched by discontinuous displacements, following the standard procedures of X-FEM.

Simplified physics of the cutting procedure

The last ingredient in the method is related to the placement of the patch L during the surgery. Once contact between the scalpel and the organ has been detected by a suitable contact algorithm (see [START_REF] Barbič | Time-critical distributed contact for 6-DoF haptic rendering of adaptively sampled reduced deformable models[END_REF] for instance, for a valid contact criterion in reduced model settings), a criterion must be set in order to determine if cutting appears, thus generating a new boundary in the domain, or not. We follow closely the criteria set in [START_REF] Bielser | Interactive simulation of surgical cuts[END_REF]. Although greatly simplified, these criteria have demonstrated to provide realistic enough results in haptic environments. A scalpel cuts along its blade, so a decomposition of the acting force as in Figure 2 is employed:

F ext D F ? C F k D F ? C F a C F n ( 23
)
Figure 2. Force decomposition at the scalpel point of contact.

A threshold value of the force F cut is considered such that lower modulus of the force F k produces static friction, but no cut. Once kF k k exceeds F cut , the finite element patch L is added to the scalpel tip. In this work, a value F cut D 3N has been considered.

In order to simplify the process and to make it simpler and (notably) faster, once the threshold value F cut is reached, a whole finite element is then cut. No cut of length smaller than the typical element size h is considered. If the finite element mesh is dense enough, this limitation does not very much affect the results. Remember that the size of the global finite element mesh, N , does not affect the size of the reduced model, r (see Equation [START_REF] Niroomandi | Accounting for large deformations in real-time simulations of soft tissues based on reduced-order models[END_REF].

Speeding up the X-FEM

The standard X-FEM technique, as presented before, provides excellent results but is, in general, somewhat heavy for real-time simulation. In this section, several assumptions are taken into account in order to alleviate these limitations.

The first one is related to the support of enriched degrees of freedom. The so-called shifted enrichment functions [START_REF] Zi | New crack-tip elements for XFEM and applications to cohesive cracks[END_REF] allow to minimize the support of enriched approximation functions and therefore the width of the superimposed finite element mesh. These functions have the simple form

i .x/ D H.x/ H i 2 (24)
where H i represents the value of H.x/ at the ith node.

The second big problem of the use of X-FEM for real-time simulation is that of the integration of the resulting weak form. The discontinuity forces to explicitly take into account the two resulting elements from a cut element and to apply standard Gauss integration on each volume. Location of the integration points in these resulting volumes is often a very heavy task. For that purpose, in this work, an alternative approach has been employed. It is based on the use of linear tetrahedra (in spite of their well-known numerical stiffness), for which exact, closed-form expressions exist for the stiffness matrix. Hexahedral meshes in the vicinity of the cut are then split into tetrahedra. Thus, the resulting stiffness matrix of an element traversed by the cut will be

K D  K uu K ua K au K aa à (25)
where K uu corresponds to the standard, unenriched element stiffness matrix, and K ua , K au , and K aa correspond to the enriched degrees of freedom. As the linear tetrahedra have constant strain (and therefore constant derivatives matrix B), we arrive at [40]

K uu ij D K ij (26)
K ua ij D  V a V ‰ aj C V b V ‰ bj à K ij (27)
K au ij D  V a V ‰ ai C V b V ‰ bi à K ij (28)
K aa ij D  V a V ‰ ai ‰ aj C V b V ‰ bi ‰ bj à K ij (29)
where K uu is identical to that of the non-enriched tetrahedron, and ‰ a.b/ represents the enrichment functions above (respectively, below) the cut plane. In the same spirit, V , V a , and V b represent the volume of the tetrahedron, or the volume above or below the cut, respectively. The formulation presented before is apt for linear FEM. In the case of large displacements or strains, more sophisticate approaches should be used. This would be imperative in the case of resections, with large displacements of the resected part, for instance. However, for the application pursued in this paper, and for the sake of speed of simulation, it is maintained as presented before. As will be shown, the results are good in spite of this formal contradiction.

APPLICATION TO THE SIMULATION OF CORNEAL REFRACTIVE SURGERY

It must be highlighted at this point that the purpose of this work is not to validate the use of reducedorder models in the context of real-time applications. This has been performed in many previous works (see for instance [START_REF] Barbič | Real-time subspace integration for St. Venant-Kirchhoff deformable models[END_REF][START_REF] Niroomandi | Real-time deformable models of non-linear tissues by model reduction techniques[END_REF][START_REF] Niroomandi | Model order reduction for hyperelastic materials[END_REF][START_REF] Niroomandi | Accounting for large deformations in real-time simulations of soft tissues based on reduced-order models[END_REF], among others). Another purpose of this work is also not to validate X-FEM technique to accurately simulate moving discontinuities in solid mechanics problems (for instance, see, among others, [START_REF] Sukumar | Extended finite element method for three-dimensional crack modelling[END_REF][START_REF] Combescure | X-FEM a good candidate for energy conservation in simulation of brittle dynamic crack propagation[END_REF][START_REF] Meschke | Energy-based modeling of cohesive and cohesionless cracks via X-FEM[END_REF][START_REF] Daux | Arbitrary branched and intersecting cracks with the extended finite element method[END_REF][START_REF] Stolarska | Modelling crack growth by level sets in the extended finite element method[END_REF]). However, the combination of both provides a unique set of features that makes the resulting technique an appealing choice. Among these features, the feasibility of real-time simulation of complex non-linear tissues (including, for instance, fiber-reinforced hyperelastic tissue, see [START_REF] Niroomandi | Real-time deformable models of non-linear tissues by model reduction techniques[END_REF]) combined with the possibility of simulating surgical cutting without remeshing (at haptic feedback rates, it is obvious that for realistic rendering, the cut must be represented accurately, but only at some 30-Hz feedback rates and not at 500 Hz). To this end, the case of corneal refractive surgery is analyzed here under the framework of the proposed method.

Astigmatism is a refractive error due to the non-spherical shape of the cornea, that is, the refractive power is not uniform in all meridians. Refractive surgery techniques are used to modify the curvature in order to repair the refractive error of the eye [START_REF] Alastrue | Biomechanical modeling of refractive corneal surgery[END_REF]. This defect may be corrected by making the cornea as spherical as possible, through the application of cuts.

In addition, properties of the cornea can be quite different between patients with the same level of pathology; therefore, the technique presented before can be seen as an efficient means to plan a patient-specific surgery that minimizes uncertainty in the results, providing the surgeon with the sensations and results he will obtain later in the true surgery.

Finite element model of the cornea

The model of the cornea employed here is based upon that presented in [START_REF] Alastrue | Biomechanical modeling of refractive corneal surgery[END_REF]. The cornea was meshed using trilinear hexahedral elements. It consisted of 8514 nodes and 7182 elements. The mesh is shown in Figure 3 in two views. The cornea is clamped at its base, giving a dome-like problem that showed buckling under some types of loads.

The cornea was assumed to be, without loss of generality, hyperelastic, with a Kirchhoff-Saint Venant behavior. More sophisticated material behaviors can also be efficiently tackled with this technique, as in [START_REF] Niroomandi | Real-time deformable models of non-linear tissues by model reduction techniques[END_REF], where a two-families of fibers reinforced hyperelasticity model was successfully employed.

The Kirchhoff-Saint Venant model is characterized by the energy function given by where and are Lame's constants. The Green-Lagrange strain tensor, E , has the form

‰ D 2 .tr.E // 2 C E W E (30)
E D 1 2 .F T F 1/ D l .u/ C .u, u/ (31)
where F D r u C I is the gradient of deformation tensor. The second Piola-Kirchhoff stress tensor can be obtained by

S D @‰.E / @E D C W E (32)
in which C is the fourth-order constitutive (elastic) tensor. The Kirchhoff-Saint Venant hyperelastic 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] Barbič | 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. The material properties of the cornea were assumed to be E D 2 MPa and D 0.48 [START_REF] Alastrue | Biomechanical modeling of refractive corneal surgery[END_REF].

For a punctual load of increasing value, the complete model gave the five modes that capture 99% of the energy of the model. These modes are shown in Figure 4. Therefore, the reduced-order model has a stiffness matrix of 15 15 (five modes for each component of the displacement vector).

With these five modes, standard reduced-order models provide an error under 20% for loads placed at different positions to the one considered in the evaluation of the modes in Figure 4 and with no updating of the stiffness matrix. This level of error is considered valid in many real-time applications [START_REF] Niroomandi | Real-time deformable models of non-linear tissues by model reduction techniques[END_REF]. If more accurate simulations are needed, the method proposed in [START_REF] Niroomandi | Accounting for large deformations in real-time simulations of soft tissues based on reduced-order models[END_REF] provides a nearly exact simulation for displacement under the scalpel tip up to 1 mm (possibly much more than that, but at the scale of the cornea, 1-mm displacement involves very large strains and no reference finite elements results could be obtained to compare with).

Simulating limbal relaxing incisions

Limbal relaxing incisions are one of the three main types of cuts made in corneal refractive surgery. These incisions are made near the outer edge of the iris and are used to correct minor astigmatism (typically less than 2 diopters). In this case, a cut is made roughly at this position. In Figure 5, a detailed view of the deformed cornea, together with the superimposed X-FEM mesh, is presented. For post-processing purposes, those elements cut by the scalpel are represented as two different finite elements in the figure, although no such elements exist in the simulation, as explained before.

The obtained displacement field is shown in Figure 6. A detailed view of the x-direction displacement field is shown in Figure 7, where the appearance of the cut can be noticed, in spite of the low magnitude of the displacement between crack lips.

In order to see the difference that a practitioner would see when dealing with the simulator, in Figure 8, a comparison is made between the solutions obtained by employing the complete model, explicitly meshing the cut by separate hexahedral finite elements, and that obtained by the combined POD-X-FEM model. The practically indistinguishable displacement field that is obtained (both are superimposed in Figure 8, in practice no difference can be noticed visually) is noticeable.

Timings

The results presented in this chapter have been obtained with a MacBook Pro laptop, running Matlab 2010a with an Intel Core 2 duo processor at 2.80 GHz and 4 Gb of DDR3 RAM memory (The MathWorks, Inc. Natick, Massachusetts, U.S.A.). Despite the use of rude Matlab code prototypes, the examples of the cornea ran at some 25 Hz, which is enough for visual real-time requirements. The use of more powerful computers, possibly running in parallel, because many of the procedures here described can easily be programmed to run in that way, or the use of GPUs seem to be a feasible alternative to improve these promising results, is still compatible with the state-of-the-art techniques in the field (see for instance [START_REF] Jerabkova | Stable cutting of deformable objects in virtual environments using XFEM[END_REF]).

CONCLUSIONS

In this paper, a novel method for the real-time simulation of surgery in haptic environments has been presented. The method is based on the use of reduced-order modeling. Reduced models are, to the author's knowledge, the only technique able to simulate at real-time feedback rates, highly complex constitutive models for living tissues (fiber-reinforced hyperelastic models, for instance). The problem with such an approach lies in the introduction of surgical cuts. The global character of the approximation functions (Ritz) precludes the possibility of employing standard methods in the literature that use very efficient algorithms for partitioning or eliminating the mesh in the neighboring of the cut. The approach here developed is based on the use of X-FEM techniques, coupled with the existing reduced model through a multi-scale-like method. Thus, the superimposed finite element mesh is capable of reproducing the displacement discontinuities produced by the scalpel, whereas the underlying reduced model is able to accurately reproduce the global behavior of the organ.

The proposed method runs at feedback rate, thus allowing to take part in a surgery simulator, together with some specialized contact detection algorithm. This opens the possibility to incorporate complex, state-of-the-art soft tissue constitutive laws into real-time simulation of surgery, thus leading the possibility of making a true second-generation simulator at hand.

Figure 3 .

 3 Figure 3. Geometry of the finite element model for the human cornea.

Figure 4 .

 4 Figure 4. Five most important eigenmodes for the simulation of the cornea. The corresponding eigenvalues are 10.77 (a), 0.014 (b), 2.58 10 4 (c), 6.53 10 7 (d), and 4.43 10 11 (e).

Figure 5 .

 5 Figure5. Superimposed X-FEM mesh at the position of limbal relaxing incision.

Figure 6 .

 6 Figure 6. Displacement field (y-direction) in the cutting procedure.

Figure 7 .Figure 8 .

 78 Figure 7. Cutting procedure. Displacement field (x-direction). Enrichment degrees of freedom are magnified 10 times to highlight the magnitude of the cut.


P Mounanga

G Bastian

G Bron

R Coué

Rapid determination of the thermophysical properties of various materials. Method of the hot plane film

A new technical approach of a traditional method for the rapid (and therefore, non interfering) determination of the thermophysical properties of porous materials likely, a priori, to undergo coupled heat and moisture transfers is described. The improvements discussed here consist in the use of a heater film suiting the "method of the hot plane film" especially well. A comparison with former works by one of the method's author is made. Finally, the updated method is applied to a variety of particularly significant materials.

Introduction

Most of the civil engineering materials (also called "geomaterials") are porous, and, therefore, can absorb, retain and release water under liquid or vapour phases. They can be subjected to coupled heat and moisture transfers. The transfer-involved independent variables are, if the materials are non-saturated, the temperature, on the one hand, and the moisture content, the suction [1] or the relative humidity, if the materials are hygroscopic [START_REF] Bastian | Drying 89[END_REF], on the other hand. Moreover, according to Moyne [START_REF] Moyne | [END_REF], a third variable (the pressure of the gaseous phase (air + vapour)) is necessary. However, at least as regards civil engineering researches, situations where this variable occurs are exceptional [START_REF] Kalifa | [END_REF]. The use of geomaterials, in particular for durability evaluations, requires to know their transfer properties. Several measurement methods for parameters characterizing heat and moisture transfers within hygroscopic materials have been developed at the laboratory [5][6][7][8]. Later, some methods suiting the case of close to saturation and weakly permeable materials have been perfected [9]. Bastian's former works on the thermophysical property rapid determination using the method of the hot plane film [10] are resumed and improved to be now used in recent [START_REF] Mounanga | Proceedings of the 16th ECTP[END_REF] or currently in progress studies conducted on very fresh cement pastes (less than 24 h). If the method is difficult to apply to a fresh cement paste because of the strong heat release due to the water-cement reaction of hydration, its relevance, on the other hand, is unquestionable for porous or not, heat stable materials with mean thermal conductivity (from 0.2 to 3 W m -1 K -1 ). The technical changes brought about in the initial method and their application to different materials are presented here. The problem of thermophysical property determination for porous materials having a certain moisture content is first evoked.

Problematics. Theory formulation 2.1 Transfer equations

By neglecting pressure variations of the gaseous phase and hysteresis phenomena, and by assuming the structural changes and the concomitant reaction of hydration to be slow and weak, respectively (which becomes true 24 h after casting for cement-based materials), the transfer equations, in the case of vertical transfers with an upward vertical Ox-axis only, can be written:

ρ 0 c ∂T ∂t = ∇ • (λ∇T + ρ l L v D θv ∇θ) ( 1
)
∂θ ∂t = ∇ • (D θ ∇θ + D T v ∇T ) + ∂K l ∂z (2)
θ, the volume moisture content and T , the temperature, are independent variables. ρ 0 is the density of the dry material (i.e., which has been incubated), ρ l is the density of water, c is the wet material mass heat capacity [J kg -1 K -1 ], and L v is the water latent heat of vaporization [J m -3 ]. D θ is the total hydrous diffusivity, D T v and D θv are vapour diffusion coefficients under temperature and moisture content gradients, and K l is the hydraulic conductivity.

λ(θ, T ) [W m -1 K -1 ], the "apparent" thermal conductivity [1,[START_REF] Moyne | [END_REF] results from the cumulative effect of pure conduction phenomena (conduction is considered as pure at 0 • C (λ 0 (θ): Fourier's law) and phase change heat releases due to vapour transfers. We then have:

λ = λ 0 + ρ l L v D T v .
As shown by Daïan [13], unlike suction (∇ • (D θ ∇θ)), ∂K l /∂z, which corresponds to gravity is generally neglected in the case of cement materials.

Moisture diffusion being much slower than the diffusion of heat within the cement-based materials studied, rapid dynamic measurement methods, which would not generate interferences in the moisture content (contrary to static methods like "the guarded hot plate"), can be used here. Equations (1, 2) then reduce to:

ρ 0 c ∂T ∂t = ∇ • (λ∇T ) . ( 3
)
If thermal interferences are weak, this equation can be linearized. We then obtain the classical heat diffusion equation:

ρ 0 c ∂T ∂t = λ∇ 2 T. ( 4
)
The determination of λ (apparent thermal conductivity), of ρ 0 c (volume heat capacity), of the ratios a = λ/ρ 0 c (thermal diffusivity) and b = (λρ 0 c) 1/2 (thermal effusivity) is carried out with a simplified experimental set-up (simple geometry and boundary conditions -thermal flux step, for instance) and by solving analytically the heat transfer equation, and, finally, by comparing the experimental and theoretical results.

Method of the hot plane film

One side of a l thick cement material plate is heated using a film generating a power P [W m -2 ] from the instant t = 0. The film/plate set is placed between two thick insulation plates considered as perfect. Two thermocouples are used to measure the temperatures on both sides of the plate (heated and non-heated faces). With regard to the configuration, the solution of equation ( 4) gives the temperature theoretical evolutions of both faces, T h and T nh , respectively. If t < 3l 2 /aπ 2 , then

T h = 2P b √ π √ t ( 5
)
and (6,7)

if t > 3l 2 /aπ 2 T h = P t ρ 0 cl + l 3λ (6)
T nh = P t ρ 0 cl - l 6λ . ( 7
)
Insofar as the temperature experimental evolutions satisfactorily agree with the theoretical predictions, the values of λ, ρ 0 c, a and b can be deduced (the "anamorphose" study of T h vs. t 1/2 , for instance, introduces the value of the effusivity, b [J s -1/2 m -2 K -1 ]). ∆ being some linear variations.

Equation (5) gives:

b = 2P ∆( √ t) √ π∆T (8)
equations (6,7) give:

λ = P l 2∆T (9)
and ρ 0 c = P ∆t l∆T (10) and consequently:

a = l 2 2∆t ( 11
)
and

b = P ∆T ∆t 2 . ( 12
)
3 Experimental setting

Heating film

A former experimental set-up has been improved for this research [10]. Data acquisition system changes are not described here. The objective here is to describe the use of higher performance heating films. Because thermal inertia of the previous films (TRESCO c and SEMME c ), traditionally used for domestic low temperature ceiling panel heating, is not insignificant, and since they suit small-scale set-ups with difficulty, the performances of an alternative paper coated with a semiconducting graphited material (TELEDELTOS paper) deserve some studies. This paper has been used for many years now to represent two-dimensional stationary thermal fields [START_REF] Alder | Représentations analogiques et homologiques dans les techniques de la chaleur[END_REF]. The paper, through which a weak current generated by silverlacquered electrodes run, and which simulates the boundary conditions of the thermal conduction problem, is used to determine equipotential lines. On conditions that the voltage, which is applied between both opposite edges, is high enough, the TELEDELTOS paper is precious to heat material samples.

Set-up

The objective of the set-up presented in Figures 1 and2 is:

to characterize the TELEDELTOS paper as heating film, i.e., to measure its electrical resistance according to its surface temperature (as a rule, the resistance also depends on the relative humidity of the surrounding air [START_REF] Alder | Représentations analogiques et homologiques dans les techniques de la chaleur[END_REF] but, since the laboratory is air-conditioned, it is considered constant here); to heat a glass plate used as reference material whose thermal properties have previously been determined using two independent methods [10]. The dimensions of the plate (as well as the samples studied later) are 141×141×20 mm (that is to say a surface of 200 cm 2 ). In fact, the plate is a combination of two 10-mm thick glass plates glued together using specific grease maintaining a good thermal contact.

Two type K thermocouples, placed on upper and bottom sides, respectively, and connected to a data acquisition station, are used to record the temperature evolution as a function of time. Data are exported in an EXCEL format. The set-up is thermally insulated using polyurethane foam, which can be cut precisely.

The wiring scheme of Figure 2 consists of a direct tension of 150 V, supplied by a regulated power system and applied to a set of three series resistors, at the terminals of which tensions are measured, of the film with the two lacquered electrodes on both opposite edges, of a shunt to determine the amount of current (approximately 30 mA) and of a rheostat regulated to a resistance value more or less equal to that of the film. Under these conditions (i.e., constant voltage), an elementary calculation reveals that although the film resistance varies a little during heating, the released power due to the Joule effect, on the other hand, remains approximately constant.

This consistency, verified by calculating the product of the voltage measured at the film terminals by the current intensity through the shunt, is essential. Flux density consistency, indeed, is one of the basic hypotheses of the thermal property measurement method.

The preliminary set-up, with which the reference glass plate is studied, is then modified so as to heat the reference sample and a sample from another material simultaneously (Fig. 3). The heating film surface being double, the voltage emitted by the power supply is regulated to 200 V (Fig. 4). Measurements of the either absolute (like for glass) or relative (as described below) values of the material thermal properties can then be carried out.

Experimental results

Heating film

Figure 5 presents the calibration curve showing the relationship between the resistance of a 141-mm TELEDEL-TOS paper square strip and its temperature.

Figure 6 confirms the consistency of the released power. Table 1 gives a review of the experimental results of the glass thermophysical properties and a comparison with former results achieved with the same sample.

Thermophysical properties of the reference glass plate

Former and new results satisfactorily agree. Furthermore, thanks to a more modern acquisition system, data processing is automated.

The comparison between Figures 8 and 9 reveals a reduction in the thermal inertia of the new heating film, which is, moreover, much easier to implement. Paper strips fitting the required dimensions are easily cut. Furthermore, considering the paper resistivity, the electrical power supply does not pose problems for small size sheets.

Application to different materials

The materials studied here are characterized by their mean conductivity (i.e., they are not insulating and not metallic). For some, which are not porous (glass and PVC), the measured conductivity is called "real". For others (wood, sand with various moisture contents, clay, fresh cement paste), which are porous, the measured conductivity is called "apparent". Wood, moreover, is anisotropic.

Among these materials, wood, dry sand and PVC are low heat conductors. In [10], Bastian has shown that, for these materials, the choice of the set-up insulator can have an impact on the thermophysical property measurements because a non-negligible portion of the heat flux released by the film can be absorbed by the insulator. The author, thus, demonstrates that, for a ratio of the insulator effusivity (b 2 ) to the material effusivity (b 1 ) lower than 0.02, the error introduced by the insulator is negligible. Beyond this value, the results must be corrected.

Determination of the error due to the imperfect insulation of both external faces

We set:

b * = b 2 b 1 ; T * = T λ P l ; τ = l 2 aπ 2 ; t * = t τ .
Let ∆T * be the temperature difference between both external faces of the studied sample. In the case of an imperfect insulation (b * = 0), for λ, ρ 0 c, a and b, we Table 1. Comparison between the experimental results obtained for glass using the new and the former set-up. [10] [10] Thermophysical units Current results obtain [10]:

Hot guarded plate Dynamic measurements

λ [W m -1 K -1 ] 1.00 1.00 0.99 ρ0c [J K -1 m -3 ] - 2 .30 × 10 6 2.21 × 10 6 a [m 2 s -1 ] - 4 .35 × 10 -7 4.53 × 10 -7 b = √ λρ0c [W s 1/2 m -2 K -1 ] - 1520 1479 b [W s 1/2 m -2 K -1 ] - 1490 1496 (direct measurement)
λ = P l 2∆T ∆T * b * =0 ∆T * b * =0 ( 13
)
ρ 0 c = P ∆t l∆T ∆T * b * =0 ∆T * b * =0 ∆t b * =0 ∆t b * =0 (14) b 1 = 2P √ t √ πT (0, t) -b 2 (15) a = l 2 2∆t . ( 16
)
Table 2. Polyurethane foam thermal properties.

Thermal conductivity

[W m -1 K -1 ] 0.030 Mass heat capacity [J K -1 kg -1 ] 840 Density [kg m -3 ] 3 5 Thermal effusivity [W s 1/2 m -2 K -1 ] 3 0
The properties of the insulator used for this research are presented in Table 2.

For instance, for the wood plate heated in the normal direction to the grain (Tab. 3), the correction on the experimental results is approximately 5%.

Application to PVC, wood, clay and sand

The corrected results obtained with these materials are plotted in Table 3.

Two different methods can be used for the determination of the various physical units: first, like for glass, by using equations (7,[START_REF] Mounanga | Proceedings of the 16th ECTP[END_REF] separately, an absolute determination is conducted, which implies that the applied flux density, P , is known; or, like in [10] where equations (8,10), within which time is eliminated, are applied to two samples of materials A and B with the same thickness and heated simultaneously. We then obtain:

-short periods of time: equation (8) gives

b B = b A (∆T hA /∆T hB ); (17)
-long periods of time: equation ( 10) gives

(ρ 0 c) B = (ρ 0 c) A (∆T hA /∆T hB ). ( 18
)
Figure 10 is an illustration of relative determination using PVC, glass being the reference material.

Relative determination results satisfactorily agree with the direct measurement results of Table 3. Globally, the thermophysical property values measured for all materials are in good agreement with data from the literature [START_REF] Farouki | Thermal Properties of Soils[END_REF][START_REF]Aide-mémoire du thermicien A 3 .E 2 .Th[END_REF][START_REF] Hazard | Mémotech -Bois et matériaux associés[END_REF][START_REF] Bastian | [END_REF][START_REF] Götz | Construire en bois[END_REF][START_REF] Lignum | L'Union Suisse en Faveur du Bois)[END_REF]. For dry sand, De Vries, mentioned by Farouki [START_REF] Farouki | Thermal Properties of Soils[END_REF], obtains a thermal conductivity of approximately 0.30 W m -1 K -1 , whereas the information on sea sand found in the Aide Mémoire du Thermicien [START_REF]Aide-mémoire du thermicien A 3 .E 2 .Th[END_REF] gives values within the range 0.35-0.43 W m -1 K -1 .

Mixing of water and sand increases the material thermal conductivity significantly because of the formation of water bridges between the grains, which further heat transfers inside the mixture [START_REF] Farouki | Thermal Properties of Soils[END_REF]. For a sand with a moisture content of 10%, the thermal conductivity obtained by De Vries is approximately 2 W m -1 K -1 . The composition and the particle size distribution, which can vary depending on the sand, accounts for the difference between this value and the value achieved here.

[J K -1 m -3 ]
The thermophysical property values achieved for PVC agree with the following values found in the literature [START_REF]Aide-mémoire du thermicien A 3 .E 2 .Th[END_REF], variations in values being due to the differences in the composition of the materials: at 20

• C, λ = 0.166 W m -1 K -1 , ρ 0 c = 1.38 × 10 6 J K -1 m -3 and a = 1.20 × 10 -7 m 2 s -1 .
Moreover, for a PVC with a density within the range 1300-1400 kg m -3 , the thermal conductivity obtained by Hazard et al. is 0.20 W m -1 K -1 [START_REF] Hazard | Mémotech -Bois et matériaux associés[END_REF].

The thermophysical properties of clay are highly dependent upon the material moisture content [START_REF] Farouki | Thermal Properties of Soils[END_REF]. Both thermal conductivity and volume heat increase when the moisture content increases. The moisture content of the saturated clay studied here is 123%. In a former study, conducted on the same material save for the moisture content, which was lower, the values obtained were λ = 0.70 W m -1 K -1 , ρ 0 c = 2.29 × 10 6 J K -1 m -3 [START_REF] Bastian | [END_REF], which is lower than the values obtained here. This difference is probably due to the difference in the moisture content.

The wooden material studied (fir) is highly anisotropic (Fig. 11). In order to examine the impact of this on the values of its heat transfer parameters, measurements are carried out on two plates with the same dimensions, one being cut with the grain and the other perpendicularly. Depending on the cutting direction, the thermal conductivity varies with a multiplying factor of three, whereas volume heat remains identical whatever the configuration. This result is not surprising since volume heat depends on the material composition only, whereas the thermal conductivity depends on both material composition and structure. The results for wood agree with those found in the literature. Thus, the thermal conductivity obtained by Hazard et al. [START_REF] Hazard | Mémotech -Bois et matériaux associés[END_REF] and Götz et al. [START_REF] Götz | Construire en bois[END_REF] for fir tree wood (with a density ranging between 300 and 450 kg m -3 for Hazard et al.) is 0.14 W m -1 K -1 and 0.12 W m -1 K -1 , respectively. Few pieces of information, on the other hand, are available in the literature regarding the thermal conductivity of wood cut with the grain. The LIGNUM (The Swiss Union for Wood) gives a value of λ ranging between 0.25 and 0.29 W m -1 K -1 [START_REF] Lignum | L'Union Suisse en Faveur du Bois)[END_REF].

Application to fresh cement paste

The tests conducted on fresh cement pastes with the setup presented in Figure 3 are not conclusive. The objective of this research was to examine the correlation between the progress of the hydration process and the possible variations of the material thermophysical properties. The exothermicity of the hydration reaction is at the origin of a sharp increase in the temperature of the sample (Fig. 12), which interferes with measurements. During the test, indeed, the heat liberated by the hydration process (variable with time and difficult to assess accurately) is added to the power released by the heating film, which makes the application of equations (5, 7) and, consequently, the determination of very fresh cement thermal properties impossible.

In order to free the measurements from the influence of hydration heat on cement paste, a variant of the initial set-up (Fig. 3) has been developed, which consists in transferring heat away through the upper surface of the sample using a cooling box. This box is made of two 141 × 1-mm square plates (Fig. 13). One is made of copper and in contact with the cement paste, the other is made of transparent PVC and can be seen in Figure 13. The sidewalls of the box and the inner baffles used to direct the direction of the water flow are also in PVC. Two copper tubes are used as inlet and outlet for the cooling water pumped from a thermostatic bath. Thanks to this device the temperature of the cement paste is satisfactorily regulated (Fig. 12).

So, thanks to the modification, one side of the cement paste sample with a thickness of l, is now heated, from the instant t = 0, by the film emitting a power P [W m -2 ], whereas the temperature of the upper surface remains constant. The whole combination is placed between two thick insulation plates. Some thermocouples measure the temperature on both sides of the plate, the heated face and the isotherm one. When applying the solution of equa- tion (3), we obtain the theoretical evolution of the temperature, T h , of the heated side, that is, according to [START_REF]Aide-mémoire du thermicien A 3 .E 2 .Th[END_REF], equations [START_REF] Götz | Construire en bois[END_REF][START_REF] Lignum | L'Union Suisse en Faveur du Bois)[END_REF]:

T h = 2P b √ π √ t, if t → 0 ( 19
)
T h = P l λ , if t → ∞. ( 20
)
Both equations are used to determine λ b, ρ 0 c and a: equation [START_REF] Götz | Construire en bois[END_REF] gives:

b = 2P ∆( √ t) √ π∆T (21)
equation [START_REF] Lignum | L'Union Suisse en Faveur du Bois)[END_REF] gives:

λ = P l ∆T ( 22
)
and consequently:

ρ 0 c = 4P (∆ √ t) 2 lπ∆T ( 23
)
and

a = lπ 4(∆ √ t) 2 . ( 24
)
In practice, both equations can be applied only when the temperature experimental evolutions agree with theoretical estimations. However, the tests conducted here reveal that the heating time (>30 min) required to reach the steady state (T h constant) is much too long for a precise follow-up of the evolution of the material thermophysical properties. The application of this method to fresh cement paste is, therefore, abandoned for a specific set-up based on the method of heated and non-heated wires [7].

Conclusion

The improvements of the experimental principle described here have made it easy to implement and capable of performing repetitive tests. This technique can be used for mean conductivity materials undergoing no sharp thermal changes (exothermic reactions, for instance).

The TELEDELTOS paper is a high quality and versatility heating film. Unfortunately, the sources of supply appear uncertain.

Results on thermophysical properties of materials, at last, are reliable and can be reproduced. A certain redundancy in the relationships presented here makes it possible to assess all properties through two distinctive approaches at least and, therefore, to lessen the possibility of error in the validation of the measurements.

The miniaturization of the set-up requires homogeneous materials. If not, like concrete containing big-sized aggregates, for instance, the experimental technique describes in reference [7] suits their study much better. In a more general way, a classical approach [START_REF] L'hermite | Méthodes générales d'essai et de contrôle en laboratoire[END_REF] consists in considering a sample size 10 times higher than that of largest of the components. For example, in the presence of aggregates with a maximum size of 20 mm, the lowest dimension of the sample should be higher than 200 mm. It is then necessary to increase both the size of the device (while keeping a reasonable thickness of sample) and the power delivered by the heating film, in order to obtain a sufficient thermal signal.

In addition, in case of strongly heterogeneous medium, the local temperature measurement may be improved by installing several thermocouples, placed in series, on the lower and the upper sides of the sample [START_REF] Bouzidi | Journée AFM, Identification des propriétés mécaniques de matériaux et de structures à partir de mesures de champs[END_REF]: it is thus possible to obtain a stronger thermal signal and to measure an average temperature value on both sides of the sample.

In spite of these practical solutions, the measurement of the thermophysical properties of strongly heterogeneous materials remains a complex problem, which has not been completely elucidated. The theoretical resolution of this problem requires, according to us, the application of "volume averaging" techniques [START_REF] Whitaker | [END_REF].

Fig. 1 .

 1 Fig. 1. Thermal scheme.

Fig. 2 .

 2 Fig. 2. Wiring scheme.

Figure 7

 7 Figure 7 presents the evolution of the temperatures, T h and T nh , according to time, from which the values of λ, ρ 0 c, then a and b are deduced. Similarly, T h vs. √ t and finally effusivity b again are obtained from Figure 8.Table1gives a review of the experimental results of the glass thermophysical properties and a comparison with former results achieved with the same sample.Former and new results satisfactorily agree. Furthermore, thanks to a more modern acquisition system, data processing is automated.The comparison between Figures8 and 9reveals a reduction in the thermal inertia of the new heating film, which is, moreover, much easier to implement. Paper strips fitting the required dimensions are easily cut. Furthermore, considering the paper resistivity, the electrical power supply does not pose problems for small size sheets.

Fig. 3 .

 3 Fig. 3. Thermal scheme with both heating films connected in series.

Fig. 4 .

 4 Fig. 4. Wiring scheme with both heating films.

Fig. 5 .

 5 Fig. 5. Evolution of the resistance of the TELEDELTOS film during heating.

Fig. 6 .

 6 Fig. 6. Variation of the electric power released within the heating film during test.

Fig. 7 .

 7 Fig. 7. Evolution of temperatures T h and T nh vs. t.

Fig. 8 .Fig. 9 .

 89 Fig. 8. Evolution of temperature T h vs. t.

Fig. 10 .

 10 Fig. 10. X and Y recordings of glass and PVC temperatures.

Fig. 11 .

 11 Fig. 11. Photographs of the wooden plates showing their anisotropic features.

Fig. 12 .

 12 Fig. 12. Temperature evolution of cement paste (a 141×141× 20-mm plate) kept isotherm thanks to cooling system, on the one hand, and in semi-adiabatic conditions, on the other hand, during hydration.

Fig. 13 .

 13 Fig. 13. Photograph of the cooling box added to the set-up for the study of fresh cement paste.

Table 3 .

 3 Experimental results. 10 -7 1.62 × 10 -7 4.33 × 10 -7 2.53 × 10 -7 2.31 × 10 -7 11.1 × 10 -7 [m 2 s -1 ]

	Wood