A micro-structured continuum modelling compacting fluid-saturated grounds: the effects of pore-size scale parameter

. We show the dependence of field evolution equations on the micro-structure parameter l and compare our model with the homogenized asymptotic ones. The consideration of the pore size 1 allows us to forecast the onset of micro-displacement waves as a consequence of a ground settling and to suggest a possible description of the genesis of certain microearthquakes [5] [6].

Introduction

In [START_REF] Dell'isola | Dynamics of solids with micro-periodic non-connected fluid inclusions[END_REF] a new way is suggested to study the effect of the length scale parameter -characterizing the size of the pores -in the mechanical interaction of a fluid with a porous matrix.

The method used there (called in [START_REF] Wo~niak | Refined macro-dynamics of periodic structures[END_REF] macro-modelling procedure) is very powerful and relatively easy to apply. Here we want to use it to characterize the dynamical evolution of a system 5 P constituted by a porous solid visco-elastic matrix with non connected periodically distributed inclusions saturated by a viscous fluid. It is well known that the study of this system is very difficult from the mathematical point of view if it is modelled in a refined way: i.e. by means of two distinct 3-D Cauchy continua occupying, in every configuration, two disjoint but geometrically tightly nested regions of space.

Indeed the complex geometrical structure of the system leads to a coupled (by suitable boundary conditions) system of partial differential equations (PDE) for the displacement field of the solid matrix and for the velocity field of the included fluid. To overcome these difficulties different methods are available in the literature. Here we will mention three of them:

(i) In [START_REF] Bensoussan | Asymptotic analysis for periodic structures[END_REF] an asymptotic analysis is proposed which leads to the introduction of homogenized continua as models for considered periodically inhomogeneous systems. These continua are characterized in terms of effective coefficients appearing in their constitutive equations. These effective coefficients are found by means of a complex mathematical procedure involving the solution of a boundary value problem for a prototype PDE in the periodicity cell. (ii) In [START_REF] Capriz | Continui con microstruttura. ETS Editrice Pisa. Continua with microstructures[END_REF] -instead -a direct approach is considered. The kinematics of the substantial point of the continuum is described in terms of some descriptors of the micro state of the system. In order to find the evolution equations for these micro-state descriptors some balance equations have to be postulated in every considered instance.

(iii) In [START_REF] Dell'isola | Dynamics of solids with micro-periodic non-connected fluid inclusions[END_REF], [START_REF] Woiniak | Microdynamics: continuum modelling the simple composite material[END_REF], [START_REF] Germain | Cours de m6canique des milieux continus (tome 1)[END_REF], [START_REF] Di Carlo | Continuum modelling of beam like latticed truss: Identification of the constitutive functions for the contact and inertial actions[END_REF], although using different mathematical tools, a kind of compromise between the first and second approach is suggested. Refraining from a detailed description of the displacement fields at the micro-level, the evolution equations for the micro-state descriptors are deduced by assuming that they determine with a certain approximation the micro-fields introduced in the refined 3-1) description. This deduction is possible by the following two steps: a) First, suitable micro shape functions are chosen (respecting the criteria found in [START_REF] Wo~niak | Refined macro-dynamics of periodic structures[END_REF], [START_REF] Woiniak | Microdynamics: continuum modelling the simple composite material[END_REF]). b) Second, an expression for the micro-fields in terms of micro shape functions, the micro-descriptors and the macro-displacements is used in the weak form of force balance.

As an important by-product macro constitutive equations in terms of micro shape functions and micro constitutive equations are determined.

The first two methods show some drawbacks:

(i) In the method mentioned in (i) the solution of the prototype PDE in most cases is very difficult to construct because of the geometry of the periodicity cell and is solved, in general, only by means of numerical methods. Finally, as homogenized asymptotic models are obtained with a limiting procedure for vanishing periodicity cell size, they can account for the effects neither of length scale parameters nor of any structural change at the micro-level induced by deformation. (ii) In the method mentioned in (ii) one needs to postulate suitable balance equations from which deducing the evolution equations for micro-descriptors. No general axiomatic format seems presently available to supply valid criteria to select the appropriate form for these balance equations. For a more detailed discussion of this point we refer to [START_REF] Williams | Review of Continua with microstructure by G. Capriz[END_REF], [START_REF] Di Carlo | A non standard format for continuum mechanics[END_REF]. (iii) On the other hand the procedure proposed in [START_REF] Wo~niak | Refined macro-dynamics of periodic structures[END_REF] supplies a standard method for simplifying the refined description for ~ by means of the introduction of a micro-structured continuum, modelling it at a macroscopic level. In particular the aforementioned procedure allows for an easy determination of both 9 the evolution equations for the macro displacement and the relevant micro descriptors introduced for describing the kinematics of the micro structured continuum characterizing the coarse behavior of ~, 9 and the constitutive equations for the micro structured continuum in terms of the geometry of ~ at the micro-level and the constitutive equations of its micro-component.

The method used here is based on the deduction of evolution equations for both macro displacements and micro descriptors from the classical theory of Cauchy continua, assumed to be valid in the refined description at the micro level. The micro-displacement field for these continua is found -in principle -by using force balance law, when the geometry of the solid matrix, its interaction modalities with the field phase, and the constitutive equations for both solid and fluid phases are known.

Although their determination can be very difficult, in many instances it is possible to guess some features of the micro-displacements in order to obtain a simplified system of PDE governing the motion. In [START_REF] Wo~niak | Refined macro-dynamics of periodic structures[END_REF] this guessing procedure finds a formally correct frame of reference by means of the introduction of so called '"micro-shape functions". Their properties are listed in the following Sections 2, 3 or in [START_REF] Dell'isola | Dynamics of solids with micro-periodic non-connected fluid inclusions[END_REF], [START_REF] Wo~niak | Refined macro-dynamics of periodic structures[END_REF] to which we refer. The suitably chosen micro-shape functions allow for the representation of micro-displacement as the sum of a macro-displacement and a micro-disturbance whose evolution is controlled by the micro descriptors and whose form is fixed by micro-shape functions. In this way it is possible to account for those features of the deformation of the micro-structure relevant for the evolution of macro displacement. Indeed the weak form of the force balance equation -which we assume to hold at the micro-level -once the aforementioned representation for micro-displacement is accepted, allows for the deduction of the evolution equations for macro displacements and micro descriptors. This deduction -which makes possible the change from the refined to a coarse description of 5: -is based on 9 the determination of an accuracy parameter e and of the micro-periodicity cells; 9 the assumption that macro displacements and micro descriptors are slowly-varying fields.

More precisely macro displacements and micro descriptors are assumed to be e-macro functions, i.e. functions whose variation in every periodicity cell is e-negligible; 9 the possibility, in the macro description of 5 p, of neglecting -in the weak form of force balance -every quantity of order equal or greater than e.

While the equations for macro displacements are PDE, those obtained for micro descriptors are ODE. Therefore -in the sense used by [START_REF] Capriz | Continui con microstruttura. ETS Editrice Pisa. Continua with microstructures[END_REF] -the micro descriptors are internal state variables.

In [START_REF] Wo~niak | Refined macro-dynamics of periodic structures[END_REF] also a criterion is given to establish whether the chosen micro shape functions are reasonably guessed: for a discussion of the reliability of the micro shape functions we use in this paper we refer to [START_REF] Dell'isola | Dynamics of solids with micro-periodic non-connected fluid inclusions[END_REF]- [START_REF] Woiniak | Microdynamics: continuum modelling the simple composite material[END_REF].

In the present paper we:

9 formulate a new mathematical model for compacting fluid saturated grounds 9 obtain, in the framework of the aforementioned model, the description of the genesis of the microearthquakes related to the bradyseismic phenomena in the Phlegrean Fields in South Italy.

Geometrical structure of the solid matrix and fluid inclusions

The considered body ~ is assumed to be constructed -in a reference configuration -by translating a typical cell V. Let

,1,

We choose I as the micro-structure length scale parameter (cf. [START_REF] Wo~niak | Refined macro-dynamics of periodic structures[END_REF] and references there quoted).

In order to characterize the structure of a solid matrix with non connected fluid inclusions we consider the following partition of V: v=v~uv~, v~ c-, v~ = 0 where Vv_c V such that OVv~Vs=O.

We assume that V, Vs and VF are regular regions as defined by Truesdell [START_REF] Truesdell | A first course in rational continuum mechanics[END_REF]. Vs is occupied by a solid material while VF by a fluid one. For an arbitrary point x in the physical 3-space E we will use the following notation:

(i) V(x),= x + V, Vs(x) .'= x + Vs, Vv(x) ,= x + Vv.
Let f2 be the region occupied by the body ~ in the reference configuration. We assume that there exists a lattice A of places in E such that (where int A means the topological interior of A) A detailed (refined following [START_REF] Woiniak | Microdynamics: continuum modelling the simple composite material[END_REF]) dynamical description of the body ~ at a micro-level is obtained by introducing displacement fields whose variations in every cell V(x) (whose diameter is l) can be relevant. Instead a coarse description, i.e., a description at a macro-level (in which only some overall "averaged" properties of the body are to be considered) will be obtained by introducing "macro descriptors" of displacement which are nearly constant inside every V(x). (v) For a generic functionf defined respectively in V(x), Vs(x) and Vv(x) we denote the mean value as follows [START_REF] Woiniak | Microdynamics: continuum modelling the simple composite material[END_REF], [START_REF] Dell'isola | Dynamics of solids with micro-periodic non-connected fluid inclusions[END_REF]:

(ii) (V(x> x2)e A x A) [xl # x2 ~ int V(xl)c~ int V(x2) = 0]; O=int U lY(x),
,.,1 f ]V~I f v 1 f,=~ (f),,:= ~ f dV, <f)xS== f dV, <f ),, f dV. V(x) Vs(x) VF(x) If f is a V-periodic function, then, Vx ~ E we obtain <f),, = <f), (f)xs = (f)s, (f)xu = (f)F
where <f), <f>s, <f)P are constants.

(vi) Following [START_REF] Wo~niak | Refined macro-dynamics of periodic structures[END_REF] we call micro-shape functions a system of sufficiently regular V-periodic linear-independent functions h~A('), A = 1, ..., N, representing a functional basis for a finite dimensional subspace of the space of continuous V-periodic functions such that htA(x) ~ 0(l),

Vhla(x) e O(1), <h,A) s = (h,a) v = O, X Z E, ( 2
)
where the index I displays the dependence of the micro shape functions on the length scale parameter; O(1) represents a set of functions of I which are infinitesimal of equal or higher order than l, and O(1) is a set of functions of l whose maximum value remains finite when 1-0 0. (vii) For some e > 0, we call a real-valued function F defined on f2 an e-macro-function if

(Vx ~ ~2o)(Vzr V(x)) f(z) ~-f(x) where a ~ b means I[ a -b II --< e. ( 3
)
For each F the parameter ~ will be considered as an error implied by the choice of giving up the microscopically accurate description of the variation in V(x) of F. Although in general different F have different e (i.e. e = ee), in the following we assume that one and the same e can be chosen for all macro-functions necessary to the description of the system.

Evolution equations for macro-displacements and macro-descriptors

We assume that ~ consists of 9 a linear elastic solid matrix 9 saturated inclusions of a Stokesian incompressible viscous fluid [START_REF] Truesdell | The mechanical foundations of elasticity and fluid dynamics[END_REF]. Therefore, at the micro-level we assume that the following constitutive equations hold:

) 5
where u is the micro-displacement vector field defined in ~?s for every substantial point belonging to the matrix; v is the micro-velocity field in the fluid inclusions. Since H Vu [I ~ 1 we may set = sym (grad (v)); C is the linear elasticity fourth order tensor; D is the Navier-Stokes linear viscosity tensor; p, ts, te, fir are the hydrostatic pressure field, the Cauchy tensor for the solid and the fluid media, the Stokes viscosity in the fluid inclusions, and I is the identity tensor, respectively. Let Qs, 0r represent mass density fields of the solid and fluid components, and let b be the external body force.

We assume the following weak-form of force balance equations (see for instance [START_REF] Capriz | Continui con microstruttura. ETS Editrice Pisa. Continua with microstructures[END_REF], [START_REF] Mason | Methods of functional analysis for application in solid mechanics[END_REF]): for every ~u continuous in fl and of class C 1 in intOsuintOr, and such that ~u[0a = 0, ~e = Sym (grad (~u)).

The dynamical problem for ~ is given by Eqs. ( 4), [START_REF] Corrado | Inflation and micro earthquake activity of Phlegraean Fields[END_REF], and [START_REF] Vesuviano | Rapporto di Sorveglianza Bradisismo Flegreo[END_REF]. It describes completely the micro-behavior of the considered continuum. No qualitative study of the resulting PDE has been performed. In particular up to now it has not been possible to describe the influence of the geometry of the fluid inclusions on the average displacement of the ceils constituting the solid matrix. Only numerical solutions could be determined for particular geometrical shapes of the solid matrix.

In this paper we try the aforementioned qualitative study by refraining from a detailed description of the micro-displacement field inside the typical cell V(x). Indeed, we guess a representation for the micro-displacement field in terms of a macro-displacement field and some disturbances obtained by means of some micro descriptors and some micro shape functions. The micro shape functions account for the geometry of the inclusions and for some aspects of the mechanical behavior of the solid matrix, while the micro-descriptors account for those aspects of the micro-displacements which are expected to influence the macro-displacements. Formally this can be done by introducing the following hypotheses (for more details we refer to [START_REF] Dell'isola | Dynamics of solids with micro-periodic non-connected fluid inclusions[END_REF]- [START_REF] Woiniak | Microdynamics: continuum modelling the simple composite material[END_REF], [START_REF] Wo~niak | Micro-macro dynamics of periodic materials structures[END_REF]):

(j) The micro-displacement field u is represented as (summation over A holds):

u(x, t) = U(x, t) + hff(x) QA(x, t), x ~ s t E (to, tfl,

where hff(-), (A = 1 ..... N) are the micro shape functions describing the kinematics of the system at the microscopic level. They are to be chosen a priori in order to account for those features of the solutions at the micro-level which are relevant at the macro-level in the considered phenomena. The more suitable the guess represented by the micro shape functions is the more accurate will be the qualitative analysis they allow. For instance, in [START_REF] Dell'isola | Dynamics of solids with micro-periodic non-connected fluid inclusions[END_REF] it is shown that for empty inclusions the most appropriate micro shape functions differ from those which have to be chosen in the case of filled inclusions. More examples are available in [START_REF] Woiniak | Microdynamics: continuum modelling the simple composite material[END_REF], [START_REF] Wo~niak | Micro-macro dynamics of periodic materials structures[END_REF], [START_REF] Wo~niak | Refined macro-dynamics of periodic structures[END_REF].

Remark 1. The functions U and QA are sufficiently regular macro-functions which will be treated as new kinematical, independent variables. They will be called, respectively, the macro-displacement and the micro-descriptors of micro-displacement fields.

(jj) In modelling the macro-behavior of the system we dispence with the description of some of the aspects that can be considered micro. Consequently, we neglect all terms appearing in the equations that are of the same (or higher) order of the accuracy parameter e that characterizes the precision to which we can measure the macro-functions U(-, t), Qa(., t) and all their derivatives. In the sequel this kind of approximation will be denoted by ~. It is better to underline at this point that in this approach the micro-structure length scale parameter l is a known physical constant (independent of e) and will not be neglected. therefore the macro-motion of the matrix is the same as the macro-motion of the fluid inclusions. The hypotheses formulated for the micro shape functions are therefore suitable for describing the behavior of a solid matrix with non-connected fluid inclusions. This fact allows to use a Iagrangian description also for the fluid component of the considered system.

In this way we see that at micro-level the motion, described by u, can be regarded as a superimposition of "micro-disturbances" (hffQ a) over the "macro-motion" U. Thus the QA(., t) are those macro-quantities describing those disturbances at the micro-level which are relevant at the macro-level. Choosing different micro shape functions allows for the study of different aspects of the considered process, whose more detailed description can also be obtained by increasing the number N of micro shape functions. where 6U and 6Q A are e-macro-functions. Using (j), (jj), and (jjj) we obtain, by manipulating each term in Eq. ( 6), the weak form of the evolution equations for U and Qa. To this end observe that (for a detailed derivation of these expressions we refer to Eq. ( 61 Remark 3. The last two equations imply that both S and H A are macro-functions in f2. Moreover, we underline that the tensor S because of its definition accounts also for the dissipative component of the micro stress tensor in the fluid inclusions.

In the same way we find (see Eqs. ( 64) and (65) in the Appendix):

6h" ~u dV+ ~ 66" 3u dV ~-~ ((6) U + (6hl A) QA). (~U dV f~F ~s + f ((6 h'A) (-J + (6 htAhlB) QB). ~QA dV

(6h ahte) Qe + (Ohza) ~ + H A + pa _ (6bhlA) = 0. ( 19
)
Equations ( 18), [START_REF] Jong | Three dimensional consolidation[END_REF] together with ( 11), ( 12) and ( 16) represent the system of equations for the fields U and QA, provided that the terms in which the fluid pressure in an inclusions appears E delrIsola et al.

are known. In the following we will limit ourselves to consider incompressible fluids. This implies that

(div (u)} = 0 ~ (1} e div U + (Vh~A} F" QA = 0. ( 20
)
Remark 4. The coefficients of the system ( 18), [START_REF] Jong | Three dimensional consolidation[END_REF] and [START_REF] Verruijt | Generation and dissipation of pore water pressure, finite elements in geomechanics[END_REF] are constant. This is a consequence of the hypothesis of V-periodicity we have made on the solid matrix. Note, also, that some of the aforementioned coefficients depend on l because of the terms ~ohiAht B) ~ 12 , (~ht A) ~ 1.

Remark 5. Equation ( 19) represents a set of ODE. Therefore there is no boundary condition to assign for the QA : in the sense of [START_REF] Capriz | Continui con microstruttura. ETS Editrice Pisa. Continua with microstructures[END_REF] they must be regarded as internal state variables. Remark 6. Simple inspection of Eqs. ( 12) and ( 19) shows that it is not possible to use them to express S -as given by ( 11) -as functions of the quantities E and E only. Indeed if the length scale parameter I is not negligible, [START_REF] Jong | Three dimensional consolidation[END_REF] (after substitution of ( 12)) becomes an ODE of the second order for the variables QA. If l is negligible Eq. ( 19) reduces to an ODE of the first order for the QA's. In both cases the possibility to represent S -by means of suitable Effective Moduli (as in the mogenized theories) -as a function of macro-strain and its time derivative is lost.

A compaction problem: the effect of fluid inclusions on vibrations induced by ground settling

Let us apply the model to the case in which the solid matrix undergoes a compaction process.

In particular we wish to study the effect of the geometry of the microscopic fluid inclusions on the micro-displacement waves induced by ground settling. To our knowledge this problem cannot be handled with other methods. We assume that C and D act on E and l); as follows: We will distinguish two regimes, separate in time: the first (t e [0, to]) in which the solid compacts in the xl direction. Because of compression some micro-motions arise which are able to influence the behavior of the body at the macro level also. We study these effects in an interval following the compaction: t c [to, ts], (eventually (t z --+ oo)), analyzing the influence of the micro structure size on the macro-behavior and compare our results with those obtained with the homogenized asymptotic model. The initial conditions for t = to are determined as the final conditions of the compaction regime so that we must necessarily analyze the compaction process in some detail. In general in order to model these processes (see [START_REF] Biot | General theory of three dimensional consolidation[END_REF]- [START_REF] Verruijt | Generation and dissipation of pore water pressure, finite elements in geomechanics[END_REF]) PDE in three dimensional domains need to be solved. Here we use a semi-inverse method for solving the system of PDE (for MD) coupled with the ODE governing the evolution of micro-descriptors. We use the following micro-shape functions (see Fig. 1):

C'E=~ 0, if x~g2r (21
fhl-2 hi } --xl, if 0 < xl < hi, hi(xl) = gl xl --gi hi + 2hi T ~-, if gi < xi < ll, Ii g~ i= 1, 2, 3. (22)
This choice has been fully justified from the mathematical point of view in [START_REF] Dell'isola | Dynamics of solids with micro-periodic non-connected fluid inclusions[END_REF] - [START_REF] Woiniak | Microdynamics: continuum modelling the simple composite material[END_REF], [START_REF] Wo~niak | Micro-macro dynamics of periodic materials structures[END_REF]. Here we want to stress the physical meaning: it states that the unit V-cell and its fluid inclusion VF undergo deformations mapping them into parallelepipeds only. This assumption is well-grounded only ifa solid matrix is in contact with an incompressible or nearly incompressible fluid inclusion. For instance if we apply our modelling procedure to a solid with voids then we need as micro shape functions higher order polynomials in the x variables.

The initial phase: ground settling

We use the following boundary conditions (BC) for the macro displacement:

t UI(0 , x2, x3, t) = 0, UI(L1, x2, x3, t) = -qL1 ~o' t t UI(xl, O, x3, t) = -rlXl--to' Ul(xl, L2, x3, t) = -fix1 to (23) t t Ul(xl, x2, O, t) = -rlxl to' Uz(xl, x2, L3, t) = -rlx~ ~, U2= U3=0 on c~f2
where r/is a measure of the amount of the ground settling. Indeed the total variation of the ground level 3 at x3 = 0 is given by 3 = r/L1. To these boundary conditions we must add the following initial conditions (IC):

UI(X1, X2, X3, O) = O, D~I(X1, X2, X3, O) = --~0 Xj. to Ll' U2 (x1, x2, x3, 0) = 0, [)2(x1, x2, x3, 0) = 0, ~)2(xl, x2, x3, 0) = 0, ~23(xl, x2, x3, 0) = 0, (24) QiA(xl, x2, X3, 0) = 0, (A = 1, 2, 3; i = 1, 2, 3)
, QiA(xI, X2, X3, 0) = 0, (A = 1, 2, 3; i = 1, 2, 3).

It is very easy to verify that

U1= -3o ~ t U2 = U3 = O, P = P(t) (25) L 1 t o ' is a solution of div S -(~o) [3 + VP = 0 (26)
satisfying BC and IC.

In this way Eqs. ( 18), ( 19) and ( 20) simplify and can be written

(~o(h~) 2) Qll + Hi ~ + P~ = 0, ( 27
) (0(h2) 2) Q2 2 -}-H22 + P2 =-0, (28)
(o(h3) 2) Q33 -I-H33 + P3 = 0, (29)

The subsequent phase." micro-vibrations

We assume that in the interval [to, tl] the IC for the macro displacement and the microdescriptors are given by their corresponding final value in the settling phase. Therefore, when the

Towards an application to the genesis of micro-earthquake induced by compaction

We note that, because of the dependence of~o and/~o on I and g, we can prove that there are two ranges in the set of geometric and material constants. Bodies characterized by constants belonging to the first set are such that all initial strain deformations are damped in an exponentially decreasing way in a short "transient" time interval. On the other hand media with material constants belonging to the second class experience a longer period of damped displacement oscillations.

The quoted results open interesting possibilities towards the description of some among the phenomena described in [START_REF] Casertano | Hydrodynamics and geodynamics in the Phlegraean Fields area of Italy[END_REF], [START_REF] Corrado | Inflation and micro earthquake activity of Phlegraean Fields[END_REF] and [START_REF] Vesuviano | Rapporto di Sorveglianza Bradisismo Flegreo[END_REF]. These papers show how peculiar the seismic and volcanic activity is in the region of Phlegraean Fields in southern Italy. In particular in [START_REF] Corrado | Inflation and micro earthquake activity of Phlegraean Fields[END_REF] microearthquakes are decribed which are typical of that area.

On the other hand in [START_REF] Casertano | Hydrodynamics and geodynamics in the Phlegraean Fields area of Italy[END_REF] -in a very intuitive and somehow tentative way -some ideas are presented about the genesis of those earthquakes and the whole bradyseismi& phenomenon. The ideas are:

(i) the vertical movement of the crust in the Phlegraean Fields could be caused by the pore pressure variation in the porous fluid saturated media filling their volcanic basin. (ii) the micro-earthquakes accompanying both the positive (crustal uplift) and the negative (crustal settling) phases of the bradyseism could be a consequence of a sudden micro-strain release following by ground settling.

It seems to us that the results of the present paper allow us to perceive the formulation of a mathematical model suitable for the qualitative and quantitative substantiation of the idea mentioned in (ii). Indeed we are already able to prove that micro-displacement and macropressure damped waves can arise -under particular conditions -as a consequence of ground settling. On the other hand, some improvements of the present treatment need to be performed: in particular, a model capable to account for the phenomena occurring in solid matrices with micro-connected fluid inclusions need to be formulated. This generalization will be presented in a future paper.

We conclude by underlining that the ideas mentioned in (ii) can be found -in a more intuitive, but also more suggestive form -already in a very ancient tractate, by Lucretius:

Nunc age quae ratio terrai motibus exstet percipe, et in primis terram fac ut esse rearis subter item ut supera ventosis undique plenam speluncis multosque et rupis deruptaque saxa multaque sub tergo terrai flumina tecta volvere vi fluctus summersaque saxa putandumst, undique enim similem esse sui res postulat ipsa. his igitur rebus subiunctis suppositisque terra superne tremit magnis concussa ruinis, subter ubi ingentis speluncas subruit aetas ; quippe cadunt toti montes magnoque repente concussu late disserpunt inde tremores.

Lucretius, De Rerum Natura Liber VI, 535-547

We believe to be useful to the reader quoting here the English translation of this extract by Rouse and Ferguson Smith [START_REF] Lucretius | De rerum natura. With an english translation[END_REF]: Now attend and learn what is the reason for earthquakes. And in the first place, be sure to consider the earth below as above to be everywhere full of windy caverns, bearing many lakes and many pools in her bosom with rocks and steep cliffs; and we must suppose that many a hidden stream beneath A bradyseism is a slow quiet upward or downward movement of the earth's crust first observed in the Phlegraean Fields in southern Italy. the earth's back violently rolls its waves and submerged boulders; for the facts themselves demand that she be everywhere like herself. Since therefore she has these things attached beneath her and ranged beneath, the upper earth trembles under the shock of some great collapse when time undermines those huge caverns beneath; for whole mountains fall, and with the great shock the tremblings in an instant creep abroad from the place far and wide.

Appendix: Deduction of Eqs. ( 18) and (19) from Eq. ( 6)

We start from the weak-form of the force balance equations:

., (ts(X, t) VX ~ f2s S(X, t) ,= ~t~ --pI VX E YaF.

(60)

We obtain, by substituting for each term the expressions derived in ( 7) and ( 8), and by neglecting all terms of order equal to or higher than z,

H A = <C.VhlA> : E + <C : (Vht A @ VhlS)) .QB + <D.Vhla> : 1~ + <D : (ght a @ VhzB)> .QB, (63)

Vte(to, ty) and Vxe~2 ~ :S(x, t),= <s>, s, HA(x, t):= <sVhlA>~ s.

In the same way we have (69)

~. oil. au dv+ [. o~" au dV

  f2s == int U Vs(x), x~A is a regular region. OF,= int U 9"v(X), f2 ~ ".= {x s O" V(x) c f2}, xEA will represent, respectively, the placement of solid matrix, of fluid inclusions and of what is called the macro interior of O. (iv) The smallest characteristic length dimension L of g2 satisfies 1 --41. Le

(

  we use the conventions C.v ..= Cijkzvt, C:w ,= Cimwkz with summation over doubly repeated indices) ts = C : e + D : ~; e = Sym (grad (u)); for the solid, (4) tr = pI + 2fir Sym (grad (v)) -~ div (v) for the fluid, (

Remark 2 .

 2 Because of the definition of micro shape functions given in Section 2 we obtain

(

  jjj) In what follows we will assume, coherently with (7+ Sym (V(hzAQA)) ~ E + Sym (Vhff | Qa) de ~ dE + Sym (Vht a @ 6Q a)

  ) in the Appendix) (to -pI): 6e dV+ ~ ts : 6e dV ~ -~ (S : 6E + H A. 6QA) dV (to, t f) and Vx ~ ~2 ~ we define the following system of macro-internal forces: S(x, t):= (ts)x s + (tr --pI)x r, HA(x, t):= (s.VhtA)x s. (10) Because of the last equation we obtain the following macro constitutive equations: S ~ (C) :E + (C.Vh~ A) .QA + (D) :1~ + (D.Vhz A) .QA (11) H A ~-(C.Vh~ A) : E + (C : (Vh~ A | Vh~B)) .Q8 + (D.VhzA) : E + (D : (Vhz A | Vh~B)) .QB. (12)

  au dV+ 5 6b" au dV ~ -5 [(6b) 9 aU + (6bh,a) 9 ,~Qa] dV, (14) f2F ~s s p div (t~u) dV ~-~ (-VP. ~U + pa. ~QA) dV (15) s where p = (p)x F, pA = (pl7hta)xF, VX E O ~ (16) Therefore we obtain the following weak form for the PDE and ODE governing, respectively, the evolution of U and QA: S (S: ~E + H A. ~QA) dV+ ~ ((0) U + (oh, A) OA) " (~U dV f~ I2 + I ((6 hla) ~ + (6 hlAhIB) ~e). 3QA dV = ~ [((6b) 9 ~U + (6bh~A) 9 ~QA) + (VP. ~U + pA. ~QA)] dV (

) 2 /

 2 ~sE + 2s Tr (E) I, if x s fas I( 1 ) D:I~=, 2fiF l~-~Tr(l~)I , if x~2v 12~s( l~-l-Tr(l~)I)3 , if x~g2s.

Fig. 1 .

 1 Fig. 1. Shape and size of a typical cell, plots of used micro-shape-functions

  dV+ tr : Je dV+ ~ Osfi Ju dV+ OF~ f2s ~F DF -f Qsb . ,~u dV-f qFb . au dV. DS ~F Let us define gx e K2 -Ju dV)

  S s:JedV~ Z I s:JedV~ 2 f [ s:JE+s:(Vh'A| g~ xeA V(x) xea V(x) = ~ [<s>x : JE + <s.VhlA>x .jQa] IVI ~ 5 (S:JE + nt A. aQ A) dE, " E + (C.VhtA> .QA + <D> : E + <D.Vhta> .QA,

  ( (J + h'AQA)'( Ju + ht"JO") dV+ Z f 0((5 + hffQa) 9 (JV + h,"JO") dV xeA VF(X) xeA Vs(x) I (<0> (I + <Qh,a> QA). JU dV+ S (<~ h,A> (f + <0 h~Ah~B> QB). jQA dV, -diU + ob'htA(~0 A]dV+ Z S [0 b'~U + eb'htAb0 A]dV xeA VF(X) xeA Vs(x) ]" [<~ob). ~U + <obhza> QB). ~QA] dV, (65) A micro-structured continuum I P div (au) dV= 181 ~, S pdiv(Ou) dV= • I pdiv(OU+hAOQ A) dV xaA VF(x) x~a VF(X) Y~ I {p div (aU) dr+ pVh a. aQ A dr+ ph ~ div (aQ A) dr} x~A VF(X) = ~ {(p div (OU))x ~ + (pVhA)~ r" r176 + (pha)~ v div (~QA)} Iv[ x~A I {P div (OU) + pA. ~Qa + (ph a) div (OQA)} dV = p6Uloe -~ VP" 6U dV+ I pa. OQA dV+ (ph A) OQAlo ~ -I (V(P ha) 6QA dg = I (-VP' 6U + pA. 6QA) dV (66) g2 where we have defined p = (p) F, pa = @VhlA)J, Vxef2o. Equation (65) holds because (ph A) = 0, Eqs. (60), (63), (64) and (65) hold VOU, 65Q a : 6U[o~ = 0OE + H a. ~$QA) dV+ S ((O) ~J + (oht A) QA). 6U dV + S ((O h'a) (J + (0 h'ahtB) QB). 6QA dV D -S [((0b) 9 6U + (obh~ A) 9 6Q A) + (VP. 6U + pA. 6Qa)] dV, (0) U -(oh~ a) QA + VP + (ob) = 0 (68) <oh ahtS) Q8 + <oh A) ~] + H a + pA _ <0bh n) = 0.

H22 =f2 81U1 + g2 810a + C22Q22 + C2aQ~ a + C23Q33 + D220~2 z + D210a 1 -t-023033, (33) H33 =f3 8aUa + g3 01[-)1 -Jr-C33Q33 + C3aQa 1 + C32Q22 + 033~)33 + D31~)a 1 + 032{)22, (34) where (Q(hl) 2) = (hl)2/2 (i, j = 1, 2, 3; i r j).

To give an example of the possible application of our model we limit ourselves to the case in which

(50)

where Oi A and Qi A are, respectively, the values ofQi A and ~)A at the time t = to as obtained in the previous Section. The solution of Eqs. ( 38), (40

where

Remark 1. Qt(t) is obtained by substituting Eq. ( 51) in Eq. ( 40).

70

Remark 2. Observe that lim Qo(t) = -q # Qo(to). This means that, at the micro level, after the first phase of ground settling a second settling phase -charaterized by micro vibrations -occurs.

Remark 3. If eo 2 -4flo < 0 the solid matrix undergoes -in general -damped oscillations at the micro-level, along all three spatial directions. Therefore -in the framework o four model -we can prove that -under some particular geometrical and mechanical conditions, of the type represented by Eq. ( 47) with the reversed inequality sign -some damped displacement waves arise in the solid matrix -filled with a viscous fluid -as a consequence of a compaction process.

Remark 4. 7he pressure P(t) can be obtained by one of the two Eqs. ( 38), (39), once the expressions of Qo(t) and Ql(t) are obtained. P(t) undergoes damped oscillations if Qo(t) and Ql(t)

do.

E delHsola et al.

Homogenized asymptotic model

In this case the equations are readily obtained from ( 18) by letting l ~ 0 in [START_REF] Jong | Three dimensional consolidation[END_REF]. Equations ( 38), ( 39) and (40) then become

and with the same calculations as in the previous cases we obtain, in the phase of ground settling,

where go = ~o,/7o = rio, % = 7o, ~z = 71; therefore ~o t

)

with A3 =/7 /Joto ' On the other hand, in the subsequent phase, we have

Remark 5. Because of the dissipative terms appearing in the constitutive equations ( 4) and [START_REF] Corrado | Inflation and micro earthquake activity of Phlegraean Fields[END_REF] and consequently in Eq. ( 19) also in the case of the homogenized asymptotic model the possibility of introducing effective elasticity moduli is lost. Remark 6. Once the "pore-size" length scale parameter is neglected, i.e. once an homogenized asymptotic model is accepted, it is not possible to account for those cases in which damped oscillations arise after compaction.

6 Conclusions

Comparison with the homogenized asymptotic models

We start by comparing the results of our refined "length scale parameter" model with the asymptotic model. Calculating the limit when 1 -~ 0 we find from (42)-(45) that Thus, in this limit the solutions found in Section 4 turn into the asymptotic one found in Section 5. We observe that within the homogenized asymptotic model that part of the general solution depending on A1 cannot be obtained.


Sylvie Gibet
email: sylvie.gibet@univ-ubs.fr

Alexandre Bouenard Valoria

Sound synthesis by simulation of percussive gesture of a virtual animated agent

This paper presents a system which simulates a physically-based sound synthesis algorithm, controlled by a virtual animated agent. This system provides a global framework for studying the interaction between virtual autonomous and realistic humanoids and audionumeric processes. We are more specifically concerned by physically based simulation, both for the virtual character and the virtual instruments, since we aim at producing natural gestures thus leading to realistic musical interaction and credible synthesized sounds.

Introduction

Most interactive sound synthesis applications controlled by gestures try to express some relationship between the gestural actions and the effects on the produced sounds. Different methods are employed. For the approaches which directly manipulate acoustical parameters of sound, the main difficulty concerns the proper definition of mapping strategies between the gesture and the sound parameters. With such approaches, we generally observe a lack of realism in the produced sounds and in the interaction process. Other approaches use recorded sound clips for providing sounds corresponding to object interaction in a scene. Although this approach has the advantage that the sounds are realistic, there are many physical effects which cannot be taken into account with this technique. In particular, when interacting with percussive instruments, the location and the magnitude of the force impulse play a major role in the produced sound timbre. In recent years, there has been more and more emphasis placed on physical modeling in order to restore the causal link as the natural basis for sound synthesis. A basic hypothesis which supports this causal approach relies on the motor theory of perception which states that the objects of the perception emerge out of the image of physical processes that produce it. Following this idea, we believe that the organization of acoustical structures implicitly contains gesture patterns. Moreover, this approach makes possible the fine modeling of the interaction process between sound and gesture.

In this paper we propose a general physical framework, both for gesture modeling and sound synthesis. The sound-producing instrument (mechanical or acoustic) is represented as a collection of vibrating structures responding to external demands. The model consists in simulating the propagation of the vibrations along these structures. This approach requires the development of a numerical algorithm which simulates this wave equation. Furthermore, producing natural sounds necessitates the modeling of interaction mechanisms between the excitatory device and the vibrating structure, as well as the high-level control of such mechanisms. This physical modeling approach is placed in a broader simulation and animation framework, in which a virtual humanoid dynamically interacts with a physical simulated instrument. The idea here is to show the influence of the gesture control in the context of sound simulation. This influence can be evaluated as the subtle effects of the interaction processes when objects are colliding but also can be measured by the quality of the preparatory gesture. This simulation framework makes also possible the analysis of the spatio-temporal coordination of motion responsible for musical expressivity. In this paper we present the foundation of this system and the first simulation modules that are being developed. The paper is organized as follows. Related works are reviewed in section 2. In section 3 we present a general overview of the system. The environment for the dynamical modeling of virtual characters is described in section 4. The physical model of a specific percussion instrument is proposed in section 5. The interaction methodology is presented in section 6. Finally, we conclude with further perspectives.

Related works

Previous work concerns both the sound synthesis using physical models and dynamic animation of virtual characters.

Modeling the wave equation for vibrating objects was first proposed by Hiller and Ruiz in 1971 [START_REF] Hiller | Synthesizing Musical sounds by solving the wave equation for vibrating objects[END_REF]. Different methods have been proposed in the last two decades. Finite element methods can be used for boundary value differential equations. These methods were proposed for numerical integration of vibrating systems [START_REF] Chaigne | Numerical simulations of xylophones[END_REF][START_REF] O'brien | Synthesizing sounds from physically based motion[END_REF]. Some approaches use an assembling of masses and damped springs to model physically-based sound synthesis [START_REF] Florens | The physical model: modeling and simulating the instrumental universe[END_REF]. The main drawback with this approach is the loss of accuracy in the simulation results due to the propagation of numerical errors. Another physical modeling technique, called modal synthesis computes analytically the different modes of vibration instead of numerical integration [START_REF] Adrien | The missing link : modal synthesis synthesis[END_REF]. A software for modal synthesis called Modalys has been developed by IRCAM (France). More recently, several approaches proposed to link the vibration modes to physical models and to fit the sound synthesis with real time graphical synthesis. O'Brien and al [START_REF] O'brien | Synthesizing sounds from rigid-body simulations[END_REF] found a correspondence between arbitrary shapes and tetrahedral volume elements modeled by finite element equations. A similar approach is proposed in [START_REF] Raghuvanshi | Interactive Sound Synthesis for Large Scale Environments[END_REF], where the surface vibrations are modeled by a system of pointmasses and damped springs, which fit an analytical solution. A simpler method, using finite difference technique has been proposed in [START_REF] Gibet | Gestural control of sound synthesis[END_REF] to simulate a vibrating circular membrane with forces applied on various impact points.

Concerning human simulation, several physics environments are now available to design complex graphical scenes including physical entities. A specific physics engine, called PhysX (AGEIA) helps to construct rigid articulated bodies, to define collision shapes of the bodies, materials that define the surface properties and constraints between objects and their environment [9]. Numerous methods have been proposed to control dynamical characters. Physical models of humanoids were first designed by Hodgins et al. [START_REF] Hodgins | Animating human athletics[END_REF]. More recently, controllers based on automatics theory have been proposed to control dynamical characters [START_REF] Faloutsos | Composable controllers for physics-based character animation[END_REF][START_REF] Shapiro | Hybrid control for interactive character animation[END_REF]. In [START_REF] Marteau | Models with Biological Relevance to Control Anthropomorphic Limbs : A Survey[END_REF][START_REF] Gibet | Expressive Gesture Animation Based on Non Parametric Learning of Sensory-Motor Models[END_REF], methods based on learning algorithms have been developed, and are used to control cinematic or dynamical hand-arm systems.

In this study, we use the AGEIA PhysX framework, which aims at modeling and animating in real time several entities in a scene [9]. The humanoid is constructed by the attachment of joints between segments.

Overview of the system

An experimental platform has been built to test our humanoid animation and sound synthesis, using the physical engine PhysX. As shown in Figure 1, the functional architecture is divided into three main parts. The first one deals with the motion database. A task specification, directly extracted from captured motion controls the second module, called Animation Engine. This one is composed of a set of specific controllers which drive the mechanical model of the humanoid through the specification of forces and torques applied on specific joints. The third module, called Instrument Simulation is dedicated to the simulation of numerical vibrating objects. Both the humanoid and the instrument are displayed in real time, thanks to a rendering engine.

Γ ∈ ∀ = x t x 0 ) , ( φ
For numerical simulation, we discretize the surface with a set of uniform squares of equal size and keep only those contained in Γ. We use the finite difference method. This method approximates the Laplacian and the second derivative with respect to time with Taylor formula. This method does not make any assumption concerning the form of the solution. It has also the advantage to be very easily implemented, and to behave well numerically.

The only substantial difficulty which arises with such an approach is the real-time requirement necessary if we want to interact with the instrument using forces provided by the virtual agent.

Gesture interaction with the drum

Different modes of interaction are considered, depending on whether both systems -the stick manipulated by the virtual agent and the drum -are or not physically coupled. When these systems are separated, each mechanical system independently evolves. Their respective evolution can be described by the wave equation for the drum (the forces being null) and the motion equation for the virtual agent. When they are in contact, these equations are coupled through interaction forces equal in magnitude and acting in opposite sense (Action-Reaction mechanical principle). The condition of the separation of the two systems depends on the nature of the motion. In the analysis of simple percussive gestures we consider that the rupture occurs when the interaction forces are nearly zero.

In our approach, we make the assumption that the quality of the produced sounds is directly linked to the control of the percussionist gesture, taking into account both the planning of the gesture and the fine interaction process. The proposed methodology used to study these different levels of influence between gesture parameters and sound parameters can be stated as follows. First, we analyze the influence of the physical properties of the drum, in particular its elasticity and its shape. Second, we analyze two categories of gestures and the effect of their variation on the synthesized sounds: hitting gestures with a stick and sliding gestures with the finger. For both classes, the influence of the gesture parameters are qualitatively analyzed through the observation of impact forces and velocity profiles. The corresponding acoustical spectrum around the instant of the impact is also studied.

Hitting gestures

For hitting gestures we want to emphasize the importance of the force magnitude and the force profile, and also the effect of the point

Animation Engine

Instrument Simulation impact of the stick on the membrane. In a first step we only consider vertical forces applied on the horizontal membrane. With our gesture control model, the planning gesture is defined by two target points: the initial target represents the beginning of the gesture; the final target the planned objective. This objective can be located directly on the membrane, but it can also be located under the surface, leading to a gesture with a greater velocity at the impact point. Therefore, by programming a fictive target more or less close to the membrane, we can act on the force of the strike. The total duration to reach the target is not a priori known. It depends on the task and on the motion controller.

Furthermore, we want to analyze how the specification of different set of mechanical parameters for the arm-stick system (stiffness of the joints and inertia of the segments) as well as the physical properties of the membrane (mostly the elasticity) may influence the dynamics of the interaction. If for instance the end of the arm is flexible during the preparatory motion, a narrow peak appears in the force profile.

Sliding gestures

In order to produce sliding movements on the membrane, we specify a slightly moving target.

Here again, we can modify the trajectory of the gesture, spatially and temporally, and analyze the influence on the produced sounds.

Conclusion and perspectives

We have presented in this paper a general physically-based framework for analyzing gesture interaction between a virtual percussionist and an artificial drum. This approach requires both the design of a virtual actor and the computing of a virtual 2D membrane. The actor is modeled using the AGEIA PhysX dynamic engine, and dedicated controllers are designed and applied on the actors' joints. A methodology for studying the influence of the physical parameters of both the gesture and the instrument, as well as the influence of the preparatory gesture is presented. This is a work in progress. First results have already been obtained for the humanoid control and a simple 2D membrane. The synthesis of sound in real-time according to various hand-instrument interaction, and with little loss in perceived sound quality will require the development of acceleration techniques.

Fig. 1 :

 1 Fig. 1: Architecture of the system


Ugo Andreaus

Francesco Dell'isola

ON THERMOKINEMATIC ANALYSIS OF PIPE SHAPING IN CAST INGOTS: A NUMERICAL SIMULATION VIA FDM

Almr~--A model is proposed to describe the influence of solidification shrinkage on the final shape of cast ingots. We regard it as a starting point towards the modelling of the phenomena producing thermal stresses. The model is applied to 'nearly' purely thermal solidification processes occurring in 'long' cylindrical molds. In this instance evolution equations are greatly simplified: their solution in a closed form is found in some technologically interesting instances by means of a semi-inverse method. This solution is used to validate the numerical procedure conceived for tracking the free moving boundary and for determining the final shape of the ingot in a wider range of temperature and heat flux boundary conditions. The above mentioned procedure is based on a nonconventional space-time finite difference discretization of the governing field and boundary equations. A locally uniaxial discretization in space is worked out along either direction of an orthogonal two-dimensional grid. It is refined by iteratively solving the sequence of two uniaxial linear problems, each of which exhibits a narrowly banded coefficient matrix: thus their solutions are found via a factorization algorithm.

INTRODUCTION

In [START_REF] Chalmers | Principles of Solidification[END_REF] the phenomenology of solidification processes is discussed: melting and solidification problems occur in numerous important areas of science, engineering and industry. For example, freezing and thawing of foods, production of ice, ice formation on a pipe surface, solidification of metals and chemical reactions all involve either a melting or solidification process. Mathematically, melting/solidification problems are special cases of free moving boundary problems (FMBPs). One of the most interesting classes of phenomena considered in [START_REF] Chalmers | Principles of Solidification[END_REF] concerns the effects of solidification modalities upon the final shape and mechanical properties of cast ingots. In [START_REF] Szekely | Some mathematical, physical and engineering aspects of melting and solidifications problems[END_REF] a list of unsolved or partially solved problems arising in engineering applications is discussed, the aim of that paper being twofold: (i) to draw the reader's attention to some problems which are really relevant for their importance in applied sciences and technology, and (ii) to show how--if engineering answers to practical problems are required--the actual state-of-art makes necessary a compromise in both mathematical rigor and accuracy in the representation of physical reality.

The aim of the present paper is to use the pragmatic approach defended in [START_REF] Szekely | Some mathematical, physical and engineering aspects of melting and solidifications problems[END_REF] to face one technologically important problem (see [START_REF] Chalmers | Principles of Solidification[END_REF]): finding a model suitable for the description of slow solidification processes of cast ingots.

This model is obtained using the results found in [START_REF] Dell'isola | [END_REF] and [4], where the methods of continuum theory [START_REF] Truesdell | Rational Continuum Mechanics[END_REF][START_REF] Noll | The Foundation of Mechanics and Thermodynamics[END_REF] are extended to describe multi-phase systems with interfaces. However, the complete solution of mathematical problems arising when many thermomechanical phenomena are considered, is far from being attained. This circumstance is related to the 'impossibility' results found in [START_REF] Friedman | [END_REF]: the scope of classical analysis hardly includes FMBPs in more than 1D. Therefore we start a sequence of more and more particularizing hypotheses restricting the scope of applicability of the proposed model but trying to preserve the main features of the physical phenomenon to be described.

We finally consider 2D 'nearly purely thermal' or 'thermokinetic' processes of solidification. Using a semi-inverse method we obtain a quantitative description of those thermokinetic 1349 solidification processes occurring under some well-defined externally imposed conditions. This description becomes qualitative in nature when quoted conditions are not strictly controlled. Indeed the features of the particular class of solution which we have determined are most probably common to the solutions pertaining to more general cases. Hence some engineering answers can now be guessed using a somehow more rational argument than those available before. Such an approach (nowadays also implemented by means of numerical methods) was frequently and successfully exploited in engineering sciences: the best known case is represented by the intensive exploitation in the theory of structures of the Saint-Venant semi-inverse method. In particular our analysis allows the determination of (i) the scope of some classical semi-empirical rules (e.g. Chvonorinov rule about solidification time), and (ii) some features of the shape of cast ingot. Indeed the analytical solution we have found-assuming that (i) the mold has rectangular cross-section, and (ii) latent heat is only extracted from its vertical sides--shows that the pipes formed in the completely solidified ingot are exactly as deep as the mold, even in the presence of very small differences in mass density between the phases. This circumstance leads us to deduce that, in general, one should expect the formation of deep pipes if latent heat is extracted mainly from lateral sides. Similar considerations have guided the development of the numerical code we use to determine the final shapes of cast ingots in a wide range of heat extraction modalities. It seems that it could be easily improved, in order to 'handle' more physically accurate evolution equations, and is tested upon determined analytical solution. Such a test is not always possible as there are only very limited exact solutions to FMBPs which are concerned with allowable initial conditions and boundary conditions.

As far as front-tracking methods are concerned, finite-difference methods and the finite element technique have been used extensively for the numerical solution of FMB problems. Crank [START_REF] Crank | Free and Moving Boundary Problems[END_REF] examines a selection of methods for computational efficiency and accuracy, programming complexity, and ease of generalization to more complicated problems. These methods compute, at each step in time, the position of the moving boundary. If the solution is computed on a fixed grid in the space-time domain, the boundary will in general be between two grid points at any given time. Hence special formulae are needed to cope with the PDE and their boundary conditions, in the neighborhood of the free moving boundary (FMB). These formulae must allow for unequal space intervals. Implicit finite difference formulae are used [START_REF] Mochnacki | Solidification problem. Numerical model on the basis of generalized FDM[END_REF] and, since the position of the FMB is not known at the new time, an iterative procedure is inevitable. Alternatively, the grid itself is somehow deformed, or some transformation of variables adopted, so that the FMB is always on a grid line or is fixed in the transformed domain; special averaging procedures are to be used in the element containing the moving boundary, if the problem is discretized via adaptive meshes. Various ways of modifying the grid were proposed, all aiming to avoid the complication and loss of accuracy associated with unequal space intervals near the moving boundary. Adopting a variable time step approach, rather than using a fixed time step and searching for the boundary, one can decide to determine, as a part of the solution, a variable time step such that the FMB coincides with a grid line in space at each time level. Variable space grid methods keep the number of space intervals between a fixed and a moving boundary constant for all time; thus, the space interval varies in each time step, while the FMB is always on the same grid line. The complications due to the unequal grid size near FMB are avoided by moving the whole uniform grid system with the velocity of FMBs and so transfering the unequal interval to the end of the domain, which is a more tractable region. The space grid is adapted at each time step to construct finite elements in space and time for a non-rectangular grid [START_REF] Tran | Multi-adaptative finite element remeshing with specific reference to metal casting[END_REF]; a weak form of the heat-flow equation is solved, while elements containing the FMB have special features incorporating discontinuous interpolation [START_REF] Mitchell | The Finite Element Method in Partial Differential Equations[END_REF]. Usually, the position of the FMB is assumed to be linear in a special element, and fixed temperature distributions are assumed on its sides. Methods which combine unconditional stability (at least up to 2D) with calculation simplicity are the so-called alternating-direction schemes; these are closely related to methods known as splitting or fractional steps methods [START_REF] Richtmyer | Difference Methods for Initial-Value Problems[END_REF][START_REF] Shih | Numerical Heat Transfer. Series in Computational Methods in Mechanics and Thermal Sciences[END_REF].

In the present work the above mentioned alternating-direction method for multidimensional heat-flow problem and implicit finite-difference approximations on a fixed grid are used to solve two-phase solidification problems in 2D. The developed procedure enables the computation of boundary movement in the directions of both coordinate axes. Standard finite-difference approximations to the heat-flow equation are used at grid points far enough from the solidification boundary. Near the boundary, formulae for unequal intervals are incorporated into the equations. Interpolation formulae of Lagrangian type can be used [START_REF] Crank | Proc. of the Conf. "Moving boundary problems in heat flow and diffusion[END_REF]. To avoid loss of accuracy associated with singularities arising if the FMB is too near a grid point, localized quadratic temperature profiles are used [START_REF] Lazaridls | [END_REF]. The derivation of these equations is outlined for the two-dimensional case, numerical experiments are reported for solidification problems when the casting is a long bar of rectangular cross-section with dimensions LIL2, for various slenderness ratios Lz/LI and fixed boundary conditions. From the thermokinetic point of view the problem described in this paper is a purely two-dimensional task (2D); however the procedure can be quite easily extended to the case of time-dependent boundary conditions, temperature-dependent thermophysical properties, mixed initial-boundary-value problems. The achieved results are in agreement with those obtained by an analytical method (Section 4.3). The advantage of the proposed method is its ability to track the FMB which moves along two coordinate axes even when any coordinate line is not assumed to cut the FMB at most once. We use a simple front-tracking method, i.e. a longitudinal-transversal scheme or implicit method with variable direction, determine the temperature field of the casting and the solidification progress, and calculate the total solidification time [START_REF] Hlousek | A coherence between a castings shape and its total solidification time[END_REF]. However, in numerical applications considered here, the FMB moves essentially along one coordinate axis, and hence only the associated component of its velocity is assumed nonvanishing.

PHYSICAL AND ENGINEERING ASPECTS OF SOLIDIFICATION PHENOMENA

In [START_REF] Chalmers | Principles of Solidification[END_REF] is described how the solidification process influences the properties of solidified ingots. The importance of considered phenomena resides in the need to optimize the ingot production process. We consider solidification of pure metals cast in molds whose surfaces are solidification catalyzers, so that solidification begins exactly close to them. We aim to describe the effects of mass density variations determined by thermodynamic phenomena. They occur (1) in phase transition processes, being more relevant, and (2) in both liquid and solid phases, being not completely negligible.

In case [START_REF] Chalmers | Principles of Solidification[END_REF] the volume variation between solid and liquid phases leads to the formation of empty regions (cavities or pipes) in the mold.

In case (2) in liquid phases convective motions arise, while in solid ones a thermal strain is induced, leading either to plasticity phenomena or to residual stresses or both (cf. p. 230 of [START_REF] Chalmers | Principles of Solidification[END_REF]). In fact, if the shape of the mold does not allow uniform shrinkage of solidifying or solidified metal the resulting configuration of the ingot is not stress free.

While their volume is determined by the difference of liquid and solid mass density, the shape of cavities and pipes depend at least on (in decreasing order of importance, el. [START_REF] Chalmers | Principles of Solidification[END_REF]): (i) relative rate of heat extracted from the bottom and the sides of the mold, (ii) the thermoconvective motions arising in the liquid phase during the solidification, (iii) the field of residual stress arising in the solid phase (see [START_REF] Boley | Theory of Thermal Stresses[END_REF]), and (iv) the interracial phenomena occurring in the considered physical system.

The a priori determination of the effects of these phenomena is important in engineering applications: to control the quality of the final product one must avoid formation of cavities, determine the placement of pipes, minimize residual and thermal stresses, and limit the volume of regions where plasticity occurs.

We do not try to model all listed phenomena: we use the results found in [START_REF] Dell'isola | [END_REF] and [4] in view of the applications which, for the time being, seem possible. We model all the interfaces present in the system simply by means of discontinuity surfaces to which no physical properties are attached. In particular we will not be able to forecast the effects of surface tension on the mechanical properties (cf. pp. 291 ft. of [START_REF] Chalmers | Principles of Solidification[END_REF]) of the solidified ingot. Instead (see Section 4), we can describe 'nearly purely thermal' or 'thermokinematic' 2D solidification.

More particularly, we consider molds whose shape is very 'long' (i.e. a cylinder whose axis is much longer than the diameter of the section), so that the solidification phenomena is described by specifing what occurs in a 2D section of the molds. The importance of this case in the applications is proved by the interest paid in [START_REF]Free Moving Boundary Problems: Theory and Applications[END_REF] to 2D problems in the modelization of continuous casting.

FORMULATION OF THE MODEL

In the first column of the following table we name the compact surfaces which model the placement of the interfaces specified in the second column: AA liquid phase--surrounding ambient ZA solid phase--surrounding ambient AK liquid phase--catalyzing inner surface of the mold XK solid phase---catalyzing inner surface of the mold All quoted compact surfaces can move during solidification. Therefore enormous mathematical difficulties arise when a complete set of evolution equations (EE) for all the fields describing the state of considered systems is formulated. Indeed a FMBP for EE is to be formulated where five free moving boundaries appear. Let K be the whole catalyzing inner surface of the mold. We have

AK t3 Y K _ K. (3.1)
We assume that: (H1) K is a fixed known surface; (H2) there exists a fixed compact domain D (in the sense of TruesdeU [START_REF] Truesdell | Rational Continuum Mechanics[END_REF]) in which the solidification process occurs whose boundary OD includes (and eventually coincides with) K. In processes occurring in the absence of gravity K is closed, in order to bound the domain D, while in the presence of gravity the 'upper' part of 0D includes the interface solidifying metal-surrounding ambient (e.g. air or supplementary refrigerating system).

In the absence of an external supply of liquid metal, XA is not empty and does not move anymore due to solidification. Indeed, became of mass density change occurring in phase transition, (i) the liquid phase in general is not covering all the boundary of S not in contact with K, and (ii) a flow arises in the liquid phase at least to assure mass balance.

We call S and L the domains occupied by liquid and solid phases. These sets are time dependent. We assume that the state of the considered system is determined if the following fields are known:

O: temperature field defined in S and in L p: mass density field defined in S and in L v: velocity field defined in S and in L u: displacement field in S (3.2)

We remark that a problem arises when defining the displacement field in the solid phase, as no global reference configuration can be introduced for a solid which before the process did not exist. Following what is usually assumed in the literature (for instance see [START_REF] Boley | Theory of Thermal Stresses[END_REF]) we regard u as the displacement from the place occupied by the substantial solid particle at the instant in which it became solid. We can formulate for the fields (3.2) the EE local consequences of the corresponding global balance laws. We denote with an upper dot the material time derivative, with n the unit normal vector field to considered surfaces whose normal speed is c, and with the symbol [FL,~ ~p = s,1 the limit value of the field evaluated on the considered surface respectively from the solid and liquid side.

Balance of mass iJ + pdivv=OinSO L;fsp + fLp= M;

(3.3) {[p(v.n -c.)]~ -[p(v.n -c.)],} = 0 on ZA (3.4)
Equations (3.3) and (3.4) ensure that there is no mass production inside bulk phases and at the IgA interface; the total mass M of the solidifying specimen will be assumed constant. stating that (i) solid-ambient, solid-catalyzer, liquid-catalyzer interfaces do not move, and (ii) external ambient exerts an assigned contact force on the liquid phase. Some initial conditions have to be added to the listed equations, together with all constitutive equations, in order to well-pose the mathematical problem representing the searched model.

Balance of force

Balance of energy and entropy

To our knowledge this problem is not solved: no general theorem of existence and uniqueness seems to be available and no known numerical method (see [START_REF]Computational Modelling of Free and Moving Boundaries I![END_REF]) seems to be comprehensive enough to be applicable in the considered instance.

THERMOKINEMATIC 2D-SOLIDIFICATION: EXACT SOLUTIONS BY MEANS OF A SEMI-INVERSE METHOD

We call a process thermokinematic or nearly purely thermal 2D-solidification if the following hypotheses are verified: (H1) Solid and liquid phases are incompressible and their densities are uniform and constant. We call PL and Ps respectively the liquid and solid phase mass densities. (H2.1) Some symmetry considerations reduce the problem to a two-dimensional one: the considered mold has the shape of a long beam with constant section. Then the spatial coordinates reduce to (xl,x2) and the surfaces modelling the interfaces can be determined by curves in the plane of the section. (H2.2) The section of the mold is rectangular and its sides are parallel to the (x~,x2) axes.

Let 3' be the curve representing the union EA U YA. Latent heat is extracted mainly from the lateral (i.e. parallel to the x2-axis) side of the mold in such a way that 3' can be determined always during the solidification process as the graph of a function f of the variable x2. (H3) When phase transition occurs the total volume of the present matter varies because of the mass density difference between the phases. In solidification processes the volume decreases and because of the mass conservation law a velocity field has to arise in the liquid phase in order to allow such a volume change. This velocity field has negligible effects on all thermomechanical phenomena: it only assures the global balance of mass. Moreover the solidification process is slow enough to allow the liquid phase to flow in an equilibrium configuration once a new solid part is formed, and before the solidification process is able to subtract to the liquid phase again a relevant quantity of mass. This means that AA interface is fiat and parallel to xl-axis (i.e. horizontal): its level h has to be determined by the global balance of mass. We assume that convective velocity fields are negligible. (H4) The external ambient exerts a pressure Pe on the solidifying system and Pe >> (Ps + PL)(ll + 12)g SO that stress tensor inside the liquid and solid phase is equal to p¢l. then for every profile f there exists a unique level h verifying (4.1). For a solidification process, a family of solidified profiles has to be considered. We then introduce a family of functions f depending on the time variable t describing the profile at instant t. A particular solution for EE is obtained assuming that:

Relation between the level h of liquid phase and the solidified phase profile f. The case of flat Y.A

(1-I5) the part of 3' representing the interface Y A is, at every instant t, a vertical straight line. If hypothesis (H5) is verified and if T is the solidification time, then f(x2, T) is easily obtained using equation (4.1) and the condition f(12,T)= 0. Indeed we have:

f(x2,T)=I,(I-(~Y~ \12] ] (4.3)
with a = p2/Ap.

The proof of (4.3) is easily obtained as follows. Consider for every a E [0,12] the function go(x2) defined in the interval [0,12] as follows:

~f(x2,T) if x2 > a (4.4) go(X2): = tf(a, T) if x2 < a.
REMARK. When (H5) is verified, the function go represents the solidified profile when the level of liquid phase is equal to a. Equation (4.4)~ means that once one part of the solidified profile is no more in contact with the liquid phase it will not evolve.

Then, using M = pLlll2 and (4.1) we have: Deriving equation (4.5) with respect to the variable a we obtain the ODE whose solution, verifying the condition f(12,T)= 0, is given by (4.3).

Modalities of extraction of latent heat

The balance of energy and hypotheses (H1), (H3) and (H4) imply that latent heat is extracted only because of heat conduction and that the derived FMBP is the classical Stefan problem. We report here the Neumann solution for one-dimensional solidification (cf. [START_REF] Crank | Free and Moving Boundary Problems[END_REF]). Indeed in Section 5 we find numerical solutions for a wide class of 2D thermokinetic solidification processes to be tested on analytical solutions. Such a solution is found in Section 4.3 by means of a semi-inverse method in which the Neumann solution is needed.

Let x and t be the space and time variables. We consider at t = 0 the liquid phase occupying the [0,~] space interval with uniform temperature Or (melting temperature). Moreover we assume that for every t > 0 the temperature at x = 0 is OR < Or. Then a solidification process takes place and (i) the solid-liquid interface moves following the law

~(t) = aV~ (4.6)
where a is the solution of the following transcendent equation:

(OR -OF)Klexp( -a2/4al)

L aV~ = (4.7)
X/~alerf(a/2X/~al) in which L is volume latent heat, K~ is heat conductivity of the solid phase, and al: = K~/c~ with c~ denoting volume specific heat; 2. (ii) the temperature field O~ in solid phase is given by

(OR --Op) O,=OR erf-~a ~a ~) erf(x/2V~a ~ t).
(4.8)

An exact solution by means of a semi-inverse method

We determine a class of thermokinematic processes assigning some of the necessary boundary conditions and assuming some particular properties of the temperature fields. Analytical expressions for these fields and for the position of the free moving boundaries can be thus determined. We assume that (cf. Fig. 3):

1. at the initial instant to the interface ZA is formed on the lateral walls (LWs: i.e. two parallel sides of the rectangular section of the mold) and that the liquid temperature is equal to Or; moreover heat extraction occurs only through LWs, while the remaining contact interface of the liquid with the mold (on the basis) does not allow any heat flux; 2. at every instant t, during solidification at the ZA and AA interfaces the temperature has the constant value Or; the AA interface--which is flat and parallel to basis of the section of the mold----does not allow any heat flux; 3. also the ZA interface is flat and parallel to the LWs of the mold and its distance ~ from the LWs is given by equation (4.6); and 4. at every instant t, the solid phase in contact with one of the LWs occupies the region whose intersection with the plane of the section of the mold is bounded by the LWs and the graph of the function defined by equation (4.4) [with a = ~(t) as given by (4.6)]; moreover the temperature field O~ in the solid phase does not depend on the coordinate x2 along the direction of the LWs but only on the coordinate x t along the direction parallel to the mold basis. Thus we obtain equation (4.8) when replacing the variable x with the variable x~.

REMARKS. (i) AS a result of (4.8), the Stefan condition is verified at the Y.A interface.

(ii) the boundary condition at the Z,A interface can no longer be assigned, as it is specified by the expression for the temperature field determined in point (4) of the list above. In particular we observe that the temperature field on the interface solid-external ambient is equal to OR on LWs and to Or at its intersection with the YA interface.

In conclusion we have found a particular solution to the solidification problem formulated in the previous sections. The most striking feature of this solution is implied by the shape of the function f defined by equation (4.3): also a very small mass density difference between the phases--when heat extraction occurs only through the LWs---implies the solidification of two different ingots inside the mold. This can be regarded as a property common to all slow thermokinetic solidification processes: in these processes the shrinkage pipes are very deep. -k~o-~n r.~

In equations (5.1b), F~, and F~ denote the outer boundary regions where temperature 0 r and heat flux ~ are respectively assigned; nr,~ denotes the outer normal to the contour line a=l. F,q. For simplicity's sake, in the following we will assume a, On the FMB • between liquid and solid phases the temperature ag(~:;t)= O*,~r~, is constant and equal to the melting temperature O*; the heat fluxes are discontinuous and their difference is c*~,,,, where c* is the latent heat at melting temperature and ~,. is the FMB speed in the direction normal n. to the FMB itself; the normal pointing from solid to liquid phase is taken positively oriented. On the FMB 4) between liquid and solid phases the following boundary conditions do hold:

Discretization in the space-time domain

Let the space-time mesh llTh,h2 r (closure of mh,h2~) be generated in ¢3 = G + F and denoted mh,h2~ = mh, × mh2 × mr = {(X~.i~ = i~h~, tj =jr); i~ = 0,1,N~, j = 0,1 ..... No}.

This mesh is uniform in each of the variables (xl,x2;t). The mesh mh,h2~ is comprised, obviously, of the meet point of the planes xo=x~.~.(is=O,1 ..... N,,) and of the planes t = tj(j =0,1 ..... No). The sets of all the inner and boundary nodes of the mesh mh,h2~ are respectively denoted 0Jhth2 r = ~0h, X ~Ok2 X 0) r = {(Xl .i, ,X2,il;tj), 0 < i~ < N~, 0 <j <-No} Yh,h2r = {(X1.,,, X2.,2; tj), i = 0, 0 --<j --No; is = Ns, 0 <-j <-No; j = O, 0 <-is <-N~}.

The index a is omitted for the sake of brevity; dividing the range [0,L] into N steps by points xl <x2 < ... <xN-1 < L yields the mesh mh = {xi,i = 0,1,2 ..... N, Xo = 0, xN = L}, having step hi = x~ -x,._~ which depends upon the index i of the node x~. If h~ # ha+~ for at least one index i, then the mesh mh = m~' is called non-uniform. If h~ = constant = h = L/N for all i = 1,2 ..... N, we have a uniform mesh. In the irregular nodes of mh,h~ the difference operator applies to a non-uniform mesh [START_REF] Crank | Proc. of the Conf. "Moving boundary problems in heat flow and diffusion[END_REF]. The following differential operators should be approximated for the u. ANDREAUS and F. DELL'ISOLA Od problem at hand: the temperature time rate -~-, the temperature space Laplacian operator AO, the temperature space gradient Vd, the FMB time rate ~ the normal to the FMB n~. So far, a numerical procedure has been developed in order to solve the governing equations with both continuous and discontinuous coefficients. The total number of required calculations, i.e. the total number of arithmetical operations needed to solve the problem at the desired degree of accuracy, turns out to be a crucial point in the multi-dimensional mesh methods. The discretization schemes which combine the better features of the explicit algorithm (total amount of operations proportional to the total number of mesh nodes (N~ -1) * (N 2 -1)) with those of the implicit algorithm (unconditional stability, i.e. for any ~ and h~,h2) are called economical schemes. This basic idea is embodied in the so-called longitudinal-transversal scheme or implicit method with variable direction. This locally unidimensional method is suitable for use in a broad range of circumstances, namely thermal conduction equation with constant, variable and discontinuous coetficients in an arbitary domain G having any dimension p.

PROBLEM DISCRETIZATION

Discretization of the thermal conduction equation OO

The temperature time rate --and the temperature space Laplacian operator Ad are to be Ot approximated in the difference scheme discretizing the thermal conduction equation.

An approximated solution v/+t for t = tj+l is generated for equation (5.1a) by successively involving (for ot = 1, 2) the unidimensional equations of thermal conduction with the space-time continuity conditions between discretization directions: A,,~,v{+j + a,,_,,,_,v{,,_,,l \pc/~ a=1,2 (a=1~a-l=2) = 1,2.

(6.2)

In equation (6.2), ~p = 1 and ~ = 2 denote solid and liquid phase respectively. The operator A~ is substituted by the difference operator L~ in the neighborhood of the moving boundary [START_REF] Crank | Proc. of the Conf. "Moving boundary problems in heat flow and diffusion[END_REF].

Discretization of the evolutional equation of FMB

&P and the outward normal n to the The space gradient operator VO, the FMB time rate at' FMB are to be approximated in the finite difference scheme discretizing the evolutional equation of the FMB. We denote by v, sty) and nt~ ) the approximated solutions of temperature O, FMB ~, and FMB normal n¢, respectively. Solving a unidimensional problem on each row (column) yields the following form of the boundary conditions on the FMB (Fig. 1).

In more detail, the position and velocity vectors of a moving boundary point, and the vector normal to the FMB and pointing from solid to liquid phase have components

a,8 = 1,2, ~st~)~ = s<~)8 = , a 1 t x~-i J' =a 1 ~n¢,)~= -Of,,/Ox,_t ;(a=l~a-l=2). (6.4)
The transition temperature and the evolutional equation of the moving boundary are approximated by a = 1,2: vt~)(st,,);t ) = 0% a = 1 =), a -1 = 2; (

[/.0v,.,,

,o,o-, ---,o,o. 0£ /k---x-~/ -/k //n~)~ + -k t~-l~ = c* (6.6) \ 0xa/, \ Ox,,/i-I L\ Ox,,-t /s dx,,-i /if Ot [(k ( ]1-
The space-time continuity equations between discretization directions read as follows:

• j+l

Cj _ ¢1+1 ~ " " " v~) = vt~_o, jt~-Yt~-I) (a = 2); v~l)= v',fl m =fL (6.7)
Finally, knowledge of fo = 0 yields if+'. In equations (6.3)-(6.7) and in the following, lower indices within round parentheses denote differentiation direction, whereas Cartesian components are without parentheses.

In each equation (6.5)-(6.7) having index a the differential operator a is substituted by the difference operators L~ and As [START_REF] Crank | Proc. of the Conf. "Moving boundary problems in heat flow and diffusion[END_REF][START_REF] Tichonov | Equazioni della Fisica-matematica[END_REF]: The normal to the moving boundary is explicitly approximated by means of the weighted central difference (weighted average of left and right differences, the weight being the layer thickness): 3hi -hi+l hi -hi+t 3hi+l --hi AW~'Yi = -2hi(hi + hi+t) yl-t + 2hihi+"""~ yi +2h,+,(hi + hi+l) yi+t' instead of the difference operator A°. [START_REF] Tichonov | Equazioni della Fisica-matematica[END_REF] at the 'comer' point where the slope of the moving boundary is abruptly varying near the liquid head.

U. ANDREAUS and F. DELL'ISOLA

Discretization of the initial and boundary conditions

In order to formulate the boundary conditions for v(~), we consider an inner point x E ¢bh,h2 and we trace through it the straight line Ca parallel to the axis Ox~. We study the simpler case where Co meets F through two points P~, and P~+. We denote by Yh. (a = 1,2) the set of all points P~ and P+. If G = Go = 0 -< x~ <-L~, a = 1,2 is a rectangle, hence Yh° is comprised of nodes (ilhl,i2h2) lying on the sides x~=0(i~=0) and x~=L~,(i~,=N~), a = 1,2. The alternative boundary conditions for v(~) are obviously assigned only on Yha: v~,,~ = O(x;t) for x E yh,, a = 1,2 (6.9a)

{-kgrad[v¢,o]}..[ ~:e(.,]= [k O--v-~] Ox. J~, = gt(x;t)'e¢.) = 0~,(x;t) forx E Yh., a = 1,2 (6.9b)
where e¢,.~ denotes the unit vector in the a-coordinate direction, and the minus/plus sign refers respectively to side x. = 0 or x. = L.. At the initial time t = 0 the following condition is given:

v¢~)(x;0) = g°(x) forx E tOh,h~. (6.10)
Equations (6.9) and (6.10) uniquely define v / for any i~ = 1,2,...,N~ and x E tOh,h~. More explicitly, equations (6.9b) take the form

+ [kA+Lo = ~(x;t) forx ~ Th°, a = 1,2 (6.11)
where the differential operator 0~ has been approximated by means of the difference operator A 2 [201.

ITERATIVE METHOD OF SOLUTION

Time discretization

The jump from time j to time j + 1 is performed in two steps. First, equations (6.2) and (6.8) are solved for ot = 1 implicitly in the xl-direction and explicitly in the x2-direction, and then the same equations are solved for a = 2 explicitly in the x~-direction and implicitly in the x2-direction. We thus obtain a sequence of locally unidimensional schemes. Typically for the FDM approach this leads to a linear system of equations, which has a general band form; an effective algorithm based on the Gauss elimination which uses certain properties of general band matrices and computer data type architecture can be found in [START_REF] Dongarra | LINPACK Users' Guide[END_REF][START_REF] Cline | [END_REF]. In the locally unidimensional scheme not all the discretization directions exhibit the same 'weight'. Numerical experiments allowed to find out the following optimal ratios:

1 ' h , -~ ~ L , , h 2 -~ Z ~-\-~l / .
We write down the locally unidimensional (longitudinal-transversal) scheme ('variable direction' implicit scheme) at the (j + 1)th time and for the rth iteration in the following iterative form. Time discretization is given by j = 0,1 ..... No.

The initial conditions at t = 0 and at the beginning of the current time step (r = 0 iteration) are respectively expressed by vto),,i2o _-O(ilh~,i2h2;O), ij = 0,1 ..... N,, i2 = 0,1 ..... N2 (7.1) j+l,0 j ,ej+ 1,0 _ /rj it/+ 1,0

(2)ili 2 : V (2)ili2, Jl.i 2 --J I.i2, J2.i, = fi2.i,. m2h2 <s 2<( m2----+ 1)h2, vJ+~"t2)i.s2--0". (7.6) The solution at the rth iteration and the convergence condition of the longitudinaltransversal scheme for a = 1,2 respectively read as follows:

~/~ rr~+"°' :/÷,.~,-,: When convergence is attained (roth iteration) we have j+l . j+l __ j+l,rll fj+.l __~j+l,ro Viii2 = Y(2)ili2-V(2)ili2' J ~a,,-I --Ja,io-I"

The current height of liquid head H{~+~ '" has been evaluated by discretizing equation (3.3), which expresses mass conservation. The trapezoidal rule has been applied in order to calculate the current volume of solid phase. A two-dimensional cast of rectangular shape has been considered which is initially filled with liquid aluminium (pure metal); thermophysical properties of the material have been shown in Table 1, while phase change data are given in Table 2. A uniform mesh has been adopted with N1 = 5, and N2 = 60. The initial temperature and the boundary condition at the FMB (below the liquid head) are respectively v ° = ~(ilhl,izh2;O) = 0", il = 0,1 .... ,Nl, iz = 0,1 ..... N2

(~t)ill 2 l/J+l ~,,) = 0", a = 1,2.

Adiabatic boundary conditions have been assumed on the symmetry axis, liquid head and bottom respectively:

[} ] _ ",' z" A-v~+,1 (I)N, i~J' r = q Li,izhz;(j + 1)r "nr(~,=L,r2) = O, nr~,=t.,a) = {1,0} _ktAi-s+,.,_ [. /+1:-1. • ] ={0,1} v(2)i,h-q qhl,Ht2) ,(I + 1)z "nr(~-h~ = 0, nr(~:h~ // A+ j+t,r -- ~-,,oz'l V(2)i,0 : q[ilhl,0;(j + 1)z]'nr(~:o) = 0, nr(~:o) = {0, -1}.
Here, H = ,,(i)uJ+l"--height of liquid head at ts÷t = (j + 1)z at the (r -1)th iteration. Two cases have been investigated which differ for the boundary conditions assumed along the side: (1) given temperature uniform in the fit = 2 direction, and (2) assigned heat flux. In both cases the FMB velocity has been restricted to have only one nonzero component, namely in the fit = 1 direction. Such an assumption is exactly fulfilled in the first of the following study cases, on the basis of the boundary conditions and of the analytical solution given in Section 4.3, wheras it seems to be reasonable in the second example, and it has been validated by means of a suitable numerical test. Introducing data of Table 1 and equation (8.1) in equation (4.7) and solving in respect to the solidification coefficient a gives the relevant value 0.0146623. The total solidification time is equal to 41.8 sec, whereas its analytical forecast, obtained by equation (4.6) at ~ = LI, is 46.5 see. Comparison between analytical, equation (4.6), and numerical solidification progresses is made in Fig. 2 It is of value to note that the FMB does remain parallel to the x2-axis throughout the solidification process (see Fig. 2(b)), as predicted by the analytical solution (see Section 4.3). Furthermore, temperature of the liquid phase is kept constant at melting value and temperature distribution in x2-direction in solid phase does remain uniform.

Given temperature at mold sides

Given heat flux at mold sides

In the second sample, the rectangle sizes are L] = 0.0125 m and L2 = 0.1 m, L2/2L1 = 4.

The boundary conditions at the mold sides are (Fig. 4 A simple test has been successfully performed in order to validate the assumption of null 2-component of FMB velocity: the maximum value of the ratio between heat flux components through the FMB has been numerically evaluated at each time according to

  pt) = div T + b in S U L (3.5) where T is the Cauchy stress tensor and b is the volume density of external force; {[pv(v.n -c,) -Tn]s -[pv(v.n -c,) -Tn],} = 0 on ~A. (3.6) Equation (3.6) holds because of the absence of surface tension and related phenomena. We then add the following boundary condition: c, = 0 on ~A U ~:K; Tn assigned on AA (3.7)

  If • is the internal energy per unit mass and h is heat flux we have that:p~ = T:grad v + div h + r in S U L {[p(-~ + ,)(v.n-c,)-vTn-h.n]-[p(-~ + ,)(v.n-c,)-vTn-h.n],}(3.8)= 0 on YA. (3.9)Moreover we assume that 0 or h are assigned at every point on ~A U YK U AK. These u ~mH conditions determine how the solidification process is controlled. Finally we assume that on Y.A no entropy dissipation occurs:

  If (H2) and (H3) are verified the global balance of mass implies: F(h;f(.)): = Apf(y)dy + psf(y)dy + hl, pL--M : 0 (4.1)where Ap: = Ps -PL, l~ and 12 are the lengths of the sides of the mold parallel to the x~-and x2-axes and F is a function of h E [0,/2] and a functional of the profile fof the solidified profile. As F(0~f) < 0; F(/2;f) > 0; OF/Oh > 0 (4.2)

  p,{ J~2f(x2,T)dx2 + af(a,T) } = pL{af(a,T) + l,(12--a)}.

( 4 . 5 )

 45 REMARK. Equation (4.5) can be interpreted as follows: the mass of solidified phase is equal to the mass of the 'disappeared' liquid.

  Let G '~UG'~ be a two-dimensional domaint in the (x~,x2)-plane, F ~ UI",~ be its outer ~0 ~o boundary, and • be the FMB between two homogeneous phases: a solid phase (phase s or tWe adopt a usual index notation in which Greek minuscules have the values 1,2 and repeated indices do not mean summation. cp = 1) and a liquid phase (phase ! or ~p = 2); these phases are characterized by specific heat c,, thermal conductivity k,, and mass density p, constant in either phase. In either of the above mentioned phases the temperature O(x;t) is governed by the thermal conduction equation together with the initial and boundary conditions respectively given in the form () (;) ,,,a, x={xl,x2I~G.,O<t<-T,O.=O.(x;t);OO-~-a~ap u. a~. = . Otlap(tg.)=AO.=~O2tg*-~_~Dstg.; c)x~ -O,(x), x E ¢~0, t = 0, tg, lr,. ' = O,(x;t), x E F~,, 0 < t -< T; O~(x;0) = o = = (i,(x;t)'nr,~, x E r.~, 0 < t <--T. (5.1b) -[kgrad(O)l,o.nr,.

~

  ~,, ~. = ~(t), ~ a~o, = at O(~:;t)= O*, (5.2a) {[kgrad(O)]. -[kgrad(O)],}.n. = c*d~n. = c*~.®, dt [grad(O)lo = [V(O)]~ = oO, O(x;t) < O*,x ~ G~; O(x;t ) > O*,x E C,,. Oxa (5.2b)

(a = 1 , 2 )

 12 , ti<-t<-d ÷1 2 Ov(~)= ~ DaY(m, (6.1a) 8t ~=l via+l-,) = vi~ ) (a = 2) v{1 ) = v j, (6.1b) and the boundary and initial conditions (5.1b). The solution of such a problem is v/+1= "/+~ v(2) • Assigning v °= d°(x)gives v/÷1. In either equation (6.1) having index a we substitute the differential operators ~, and D with the difference operators A. and A,. [20]: A: "j+tv.,) = (k)[

Fig. 1 .

 1 Fig. 1. Approximation of the normal to the moving boundary.

  vt,~_t)A~_tf~ =c*A;-f~ +l. (6.8) ~,(~) O*;k,[Lov(~) -A~-tv~-lA~-tf~] ~ntn..,~v(~)

(7. 2 )__. 4 ) 7 . 2 1 $ 2 -

 247212 Space discretization in the coordinate direction a = 1 (longitudinal scheme) is given by i2 = 1,2 ..... N2-1, il =0,1 ..... N1. The equation of thermal conduction in either the solid (tp = 1) and liquid (tp = 2) phase and the evolution equation of moving boundary between solid and liquid phases are respectively A-,,j+I,, l.,j+l,r _ AO~,j+l,r-lAOpj+l.r-ll __ b [l+lvl +l.r AO4,j+l.r-lAOfj+l.r-I1 ,A-pj+l.r sLZ-,(l)¢(l~lia zx2¢(2)mli a zL2Jl,i2 j n, lLZ-~(1)sli 2 --zx2v(2)mli 2 zL2JI,i2 j = C A t f'l,i2 --(ml + 1)hi, vo)~,i 2j+l,r = 0". (7Space discretization Space discretization in the coordinate direction a = 2 (transversal scheme) is given by il = 1,2,...,N1 -1, i2 = 0,1,...,N2. The equation of thermal conduction in either the solid (,p = 1) or liquid (~p = 2) phase and the evolution equation of moving boundary between solid and liquid phases are respectively t V(2)ili2 --I V(I)ill 2 + -/~221/(2)ili2] k fAOvj+l,r A0tC/+l,r-I I-j+l.rl kl[ --AO j+l.r AOfj+l. r-I -1-I + /+l,rl ,A-pj+l,r sLZLI (I)ilm 2 --LL1J2,il "~ /"2 'V(2)ilszJ --LtlV(i)ilm 2 --ZXlj2,il --z-~ 2 ~"(2)ils2j ~-" C '/~t f2,i I __ ~j+ 1.r---f2.i.

1¢I2 = 1 ,

 1 /. [ fJa:ilj rl'-I ] v ~,~)~,i2, .,.,i. , fll~ -I~I= + 1 <-~'where E is a suitable tolerance, and ~/~,~/~ are defined as follows: Hi;)'"-h2 <-hj~: <_ Hi~'.'; 0 <-j2.~, <-h2, ,,~) -h2 <-:2:, <-rt<2) .

( 8 . 1 )Fig. 2 .

 812 Fig. 2. (a) Solidification progress: (b) intermediate shapes for given temperature at mold sides at t = 20 sec.

Fig. 3 .

 3 Fig. 3. (a) Final shapes for given temperature at mold sides; (b) time-histories of FMB average speed for given temperature at mold sides.

  (a). Figures 2(b) and 3(a) show the comparisons between (i) intermediate shapes (at t = 20 see) and (ii) final shapes obtained in an analytical [equation (4.3)] (solid line) and numerical (thin line) way.

Figure 3 (

 3 b) shows the comparison between the time-histories of FMB average speed obtained in analytical, [equation (4.6) divided by t], (thick line) and numerical (thin line) way.

Fig. 4 .

 4 (a) Heat flux given at mold sides; (b) time-history Of liquid head level for given heat flux at mold sides.

' 2 where t 7 '

 27 lV(I)o/2 = q[O, i2h2;(j + 1)v]'nrtx,.o) = q~2, nr(x,-o) = { -1,0}, nr(x,=2L,) = {1,0}. The boundary condition at the solidified frontier above the liquid head is eft(----z ] --tJ, z v~,)~,,= 01(g~,i2h2;(j+ 1)'r) = O* \2aJ sl -yl.,,, j>~,, : ,2h2 > H(~'; , z = ~t~t i , a = ~; Ingot half width [m] (xl-axis) Fig. 5. (a) Time-history of liquid head average speed for given heat flux at mold sides; (b) final shape for given heat flux at mold sides. where tJ'2 is the solidification time of the x2 = i2h2_ mesh line, ~1 is the xl-ordinate of the /4s'2'r-I is the height of liquid head at time solidified frontier at the x2 = i2h2 mesh line, and "'0) t z, =]~2v at the (r-1)th iteration; v~2},j2 = 02(ilhl,sz;(j + 1)l') , is solidification time of the x, --i~h~ mesh line, ~z is the x2-ordinate of the solidified ": -7 frontier at the xl = i~h~ mesh line, and .us'"'. (2) is the height of liquid head at time t s', =j~,v at the (r -1)th iteration. The total solidification time is 16.9 sec. Figures 4(b) and 5(a) illustrate time-histories of level and average speed of liquid head respectively; Fig. 5(b) shows the final shape obtained in the numerical approach.

Table 1 .

 1 Thermophysical properties of aluminium[START_REF] Brandes | Smithells Metals Reference Book[END_REF]

		Mass density	Specific heat	Thermal conductivity
		[Kgm 3]	[J°K" Kg "l ]	[W°K"m '' ]
	Liquid phase	2385	1080	100
	Solid phase	2700	1080	230

Table 2 .

 2 Change of phase[START_REF] Brandes | Smithells Metals Reference Book[END_REF]

	Melting temperature [*K]	933.5
	Formula weight (Wf)	26.9815 (see [24])
	Latent heat [JKg-atom "l]	10.47/W C I0 +t'
	Volume change during melting	6.5
	aV m = [(Vliquid --Vsolid)lVsolid] (°/o)

where the temperature space rate has been approximated by Li- [START_REF] Crank | Proc. of the Conf. "Moving boundary problems in heat flow and diffusion[END_REF] and A ° [START_REF] Tichonov | Equazioni della Fisica-matematica[END_REF], operators in the 1-and 2-directions respectively. This maximum ratio remained less than a few percent.