F Laouafa

P Royis

An itérative algorithm for finite élément analysis

Keywords: finite élément method, mixed élément, GMRES algorithm, strain algorithm. 1

In this paper we state in a new form the algebraic problem arising from the one-field displacement finite élément method. The displacement approach, in this discrète form, can be considered as the dual approach (force or equilibrium) with subsidiary constraints. This approach dissociâtes the non-linear operator to the linear ones and their sizes are linear functions of intégration rule which is of interest in the case of reduced intégration. This new form of the problem leads to an inexpensive improvement of F.E.M. computations, which acts at local, elementary and global levels. We demonstrate the numerical performances of this approach which is independent of the mesh structure. By using the GMRES algorithm, we build for nonsymmetric problems, a new algorithm based upon the discretized field of strain. The new algorithms proposed are more closer to the mechanical problem than the classical ones because ail fields appear during the resolution process. The sizes of the différent operators arising in thèse new forms are linear functions of intégration rule, which is of great interest in the case of reduced intégration.

INTRODUCTION

In structural and solid mechanics analysis, the most frequently used finite élément method is probably the one based on the one field variational formulation involving the displacement field u(x,t) 10 > 17 > [START_REF] Zienkiewicz | The finite element method. Basic formulation and linear problems[END_REF] which leads classically to the well known algebraic problem : find U G M N such that :

K-U = F KeE M , FeJR N
where t U = {«i, u^, • • •, UN} is the vector of the unknown nodal displacements, in which N dénotes the total number of degrees of freedom. The vector F in the right-hand side of (1) is the vector of the nodal forces, whereas K dénotes the NxN (linear or not, symmetric or not) stiffness matrix. Classically, both K and F are built from the assembly of elementary contributions, where n e denotes the element number of degrees of freedom. In rocks mechanics or in soil mechanics as well as for material involving concrete the non-normality rule of the constitutive law leads to a non-symmetric stiffness matrices. The range of such mechanical problems involving non-symmetries is rather large and needs special attention when solving in boundary values problems. The building of the elementary stiffness matrices K e and therefore that of the global stiffness matrix K leads to an increasing computation cost and to the loss of the local mechanical informations, in particular those relating to the constitutive law, and then requires post-analysis to recover these local features.

Within the same mathematical framework we propose a new method for building and solving the algebraic problems arising from the one-field (displacement) FEM. This approach, called UDA for Unstructured Displacement Approach, can be considered as an IBI (Integration point By Integration point) method. In this method we state the problem (1) by using three operators B, H and I. The first one is the global strain-displacement operator, the second H is the 'rheological' operator the components of which are the local (i.e. by integration point) constitutive relationships. The third I is the one associated with the weak equilibrium and numerical integration rule. These three operators remain distinct and uncoupled during the whole resolution of the problem. In order to decrease the computational time linked with the updating and assembling process of the Hessian matrix, Axelsson et al. proposed an interesting and efficient factorization of the stiffness matrix (see ref. 1 and 2 for details). Our paper focus only on non-symmetric problems of solid continua and on the contribution of UDA to some iterative methods.

These methods take a great advantage of the numerical integration rules which makes possible to dissociate the kinetic, rheological (constitutive law) and static local features of the disctretized variational problem.

The structure of these operators remains mesh (h and p) independent. This approach takes advantage of the new form of the problem and of the features of the GMRES algorithm 24 ' [START_REF] Taylor | The patch test -a condition for assessing f.e.m. convergence[END_REF] , both of them providing then substancial improvements to the one-field displacement FEM. The memory storage required by the UDA method is less than that required by the classical sky line [START_REF] Golub | Matrix computations[END_REF] method, the sparse matrix method or by other methods like those based on an element by element approach [START_REF] Golub | Matrix computations[END_REF] . Thus, the computational times is greatly reduced. The above characteristics increase greatly when reduced integration is used. This new important feature is due to the fact that the parameters linked with the integration rule are dissociated during the whole resolution of the algebraic problem. Many mechanical problems require adaptive mesh refinement. In this framework the new structure of the problem increase the efficieny of the method.

In the UDA method the algebraic problem coming from the continuous mechanical problem is stated in a new form similar to the one arising from the stress-displacement two-fields mixed finite element method. In its algebraic form the displacement approach can be then considered as the dual approach (force or equilibrium) with subsidiary constraints. A unification of both one-field and two-fields methods is also proposed. From an algebraic point of view it completes the unification proposed by Malkus and Hughes [START_REF] Malkus | Mixed finite element methods-reduced and selective integration techniques: a unification of concepts[END_REF] and Zienkiewicz and Nakazawa [START_REF] Zienkiewicz | On variational formulation and its modification for numerical solution[END_REF] .

We mention that the algorithms presented in this paper are also available for a real stressdisplacement two-fields mixed finite element analysis.

2.

CLASSICAL APPROACHS

The general mechanical framework of our study lies in problems characterized by a solid continuum il subject to external body forces f(x,t) and surface forces g(x,t) prescribed on a part F 2 of its boundary F. The part Fj = F -F 2 of the boundary of the solid is fixed. Let us denote as u(x,t) the displacement field and as a(x,t) the second order tensor of the Cauchy stresses. We assume that Q is a materially simple continuum. We suppose also small transformations, so that the strain tensor e(u(x,t)) (symmetric part of the gradient of the displacement) is linear.

Eventually we restrict our analysis to quasi-static problems.

The constitutive law can be described by incrementally non-linear equations 3 ' 20 of the form :

71 " 1 xlR 71 "" 1 , J(d, x ) G ïïT™ Xn ™ ( 2
)
where ê is the strain rate vector (i.e. the vector built from the components of the strain rate tensor e), & the stress rate vector of size n an , and where x denotes the set of memory parameters. More generally the constitutive law can be an associated or non-associated elasto-plastic one written as follows:

a = H(x) • e, (e, a) G HT™ X HT™, H(x) € IR n ™ Xri ™ (3)
In this paper one can be consider equally well both of relationships ( 2) and (3). We only focus our attention on non-symmetric problem such as those encountered in soils, in a rocks mechanics or in physical problems involving non-symmetric constitutive law.

The well known one-field (displacement) variational formulation of this mechanical problem is classically stated as,

{

Find ÏÊV= {v e (H[ a) ) [START_REF] Brezzi | Mixed and hybrid finite element methods[END_REF] , v | T x -0} such that :

r f f (4)
/ (H : e(u)) : e(v) dQ = / / • v dtt + g • v dT V u G V and leads, after a space discretization, to the algebraic problem : J find U e M N such that : {S) \

The space discretization leads to the following finite sum involving the contributions of the finite elements U e 10 ' 29 nel K = ^Ke with K e e=l The numerical evaluation of K e is usually made by using a Gaussian quadrature method, so that the elementary stiffness matrices can be constructed as :

npi e K e = J2 'BefaiOHfcjOBeCarjkJwUfc (7)
fc=i with :

npi e : number of integration points used for the finite element e, 0Jk : the weight of the transformation at the integration point k, jk '• the value of the Jacobian of the transformation at the integration point k.

and the local displacement-strain relationship :

e k e =B e (x k )-u e , € k e en n ™, u e em ne B e (x k )eJR n ™ xn < (8)

The actual iterative methods 5 ' 9 ' 22 for solving (5) work generaly with the global stiffness matrix K or with its elementary contributions K e . The method described below takes advantage of the numerical evaluation of the element stiffness matrices.

= *B e HB e dfi (6)
3.

NEW FORM OF THE FINITE ELEMENT PROBLEM

The principal features of the new form of the finite element problem are linked to the numerical integration of the elementary stiffness matrices. Using a Gaussian quadrature method one can rewrite 12 ' 13 problem (1) in the form :

Find U eJR N such that :

We can also formulate the one-field problem in an equivalent stress-displacement two-fields mixed form 13;

Find (a, U) G JR N " X M N such that :

H 1 -B *BÎ 0 ( 10
)
where B is the The diagonal matrix I € IR N " xN " = J2 e ^-e contains all informations in relation with the integration rules : the weights (to) and the jacobian values (j) for all the integration points.

N c X N Global displacement-Global
where Id is the n an x n an identity matrix.

The matrix B is a sparse matrix, the components of which are the partial derivatives of the displacement interpolation functions. This operator is linear with respect to its argument U. It is of interest to remark that for non-linear problems arising, for instance, from elastoplasticity, this operator remains the same at each step of loading. So it is constant during all the iterations of the iterative process of resolution. Moreover, this operator remains constant for the whole finite element analysis of the problem if one does not proceed to the update of the node coordinates. This can be done without increasing error of evaluation under the assumption of small transformations. So only one computation is to be done to form B, whatever the number of loading steps, of non-linear iterations and the of iterations linked with the solving iterative process.

Lets us now recapitulate the sizes of the different matrices encountered in this section and the corresponding storage required.

1-Matrix B : a. size : N a X N, b. storage : fii J2 e n P^e • (nan X n e ). 2-Matrix H : a. size : N a x N a , b. storage : /3 2 E e n Vie • («<™ X n an ). 3-Matrix Î : a. size : N a X N a , b. storage : J2
The factors (3\ and /3 2 are reduction factors depending respectively of the kind of problem (plane stress or plane strain problem, axisymmetric problem, etc..) and of the symmetry properties of the constitutive relations. Their possible values are as follows :

3

For axisymmetric problems without torsion

2 n otherwise 1 1 -I if H is symmetric /3 2 = { 2 2n an 1 otherwise
When solving a boundary value problem, one have to check the quality of the solution by some energy criteria or stress equilibrium [START_REF] Zienkiewicz | Mixed methods for f.e.m. and patch test, some recent developments[END_REF] . When adaptive (h and p) refinement is performed, the modification of the mesh lead, in our method, only to an addition and/or modification of some finite elements operator. That do not destruct the structure of the algebraic problem and do not require the computation and modification of the other operators, and moreover the implementation in finite element code is very easy. This is a interesting feature when adaptive refinement is used intensively.

ONE AND EQUIVALENT TWO-FIELDS ALGEBRAIC PROBLEM

To obtain problem (9) we start from the discretized weak form of equilibrium relations which can be stated in term of stress unknowns and written as follows [START_REF] Laouafa | Analyse et développement de méthode d'éléments finis : unification et application aux lois incrémentales de type interpolation[END_REF] ' 13 :

*BÎ -a = F ( 13
)
By taking into account relation (12), we can express the weak form (13) in terms of strain :

By using relation (11) we can now write the algebraic problem (1) in its two equivalent forms :

Find U eM N such that :

HE) U = F ( 15
)
Form I (S) uda

Find (a, U) G 1R JV<7 X IR/ V such that :

Form II (S) uda H" 1 -B *BÎ 0 G U 4.1. Two-fields considerations.
Let us consider the algebraic problem written in the form (II).

Find (a, U) G JR."* X M iV such that :

\uda ) /-rV-i H -B 0 (16)
and let us multiply the first equation of the system (16) by I . We obtain

ET 1 -a-B-U = 0 o-î H" 1 -a -î B • C/= 0
Let us now consider the second equation in (16). Since I is a diagonal matrix this relation can be written as :

JD 1 • <7 = F v?
Then the system (16) takes the following form :

' Find (a, U) G B, N<r x M N such that : ÎH" 1 -(IB (sy (IB) 0 a U ( 17
)
The algebraic problem given by the system (17) arises from the displacement one-field approach. However, we recognize an algebraic problem coming from the discretization of a stressdisplacement two-fields mixed variational formulation. This last formulation can be stated as follows [START_REF] Laouafa | Analyse et développement de méthode d'éléments finis : unification et application aux lois incrémentales de type interpolation[END_REF] -13 -21 :

(VV) [START_REF] Fung | Foundations of solids mechanics[END_REF] Find (cr,u) G S X V such that :

H where S is the stresses trial functions space ; S = {sij G Imj Vi, j G {1, 2, 3} and s,j = Sji}.

The algebraic problem (17) can be also obtained from (18) under some conditions relating to the finite element spaces of displacements VjCV and of stresses S/j C S. These conditions, which are independent of the position of the stress nodes inside the elements, generalizes the ones done by Malkus and Hughes 14 and Zienkiewicz 28 ' [START_REF] Zienkiewicz | The finite element method. Basic formulation and linear problems[END_REF] . They can be stated as follows :

nel ÎH" 1 = ( 19
)
nel nel " = V / IB = V / Ve B e
where ifc e is the elementary interpolation matrix of the discontinuous stress field in the case of the two-fields mixed discretization.

<?(x) = i\> e {x) • 0~e, Va; G fi e , <7 e : elementary vector of stress components

We have then established, under algebraic considerations, that a one field algebraic problem can always be expressed in its equivalent two-fields algebraic form. From a certain point of view the displacement one-field algebraic problem without subsidiary constraint is equivalent to its stress dual algebraic problem with subsidiary constraints. These constraints are the discretized weak form of equilibrium. We remark that they are not formulated in a local form, but in a weak one.

Theorem : The stationary condition without subsidiary constraint of the discrete potential energy functional U$(U) (U G JR N such that u = 0 on Ti),

I[ d p (U) = \ l U -K-U -f U -F
is equivalent to the stationary condition with subsidiary constraints of the discrete complementary energy functional Ilf (S),

nf (S) = i *E • î H" 1 • S, S G K K = {S G JR N ° such that *BÎ • S =
So, the equivalent discrete two-fields formulation is obtain by imposing the stationary condition of the discrete Lagrangian £ d (E, V), obtained by introducing the subsidiary constraint E G K

C d (E, V) = nf (S) -V • (*BÎ • E -F
The algebraic form (II) of the displacement finite element formulation is interesting for several reasons. The first of them is that for a certain class of non-linear constitutive equations 4 ' 20 written as :

e(x) = J(x) • cr(x), x G Q.
the constitutive operator J coming from J and built like H does not need to be explicitly inverted. This constitutes a great difference and provide substantial advantages in comparison to the classical approach, because of the decreasing of calculations and of the decreasing of errors due to that inversion in the case of ill conditioned constitutive relations. The second reason, which holds also for the form (I) of the algebraic problem, is that the performances of the algorithms described below increase if one adopts reduced integration rules.

Considering the form (II) of the algebraic problem, we can easily give a condition for integration rule. Two levels of condition are required.

1-The global level, which ensures the injectivity of B :

nel e n u -bel e=i n<yn
where bel is the number of components of U prescribed.

2-The element level, which ensures the injectivity of B e :

. n Ue -bcl e npi e >

where bcl e is the number of components of U e prescribed. As established by Zienkiewicz [START_REF] Zienkiewicz | Mixed methods for f.e.m. and patch test, some recent developments[END_REF] ' 28 these conditions give some restrictions for reduced integration rule. As we can see by considering the algebraic forms (II) the reduced integration rule can be interpreted, from a two-fields mixed point of view, by the relaxation of the conditions of equilibrium. A physical interpretation can be done by considering the stress-displacement two-fields approach, and one can linked the LBB 2 condition to the integration rule.

5.

FIRST COMPUTATIONAL ADVANTAGES

In this section we describe some features of the new approach proposed, which are independent of the method used to solve the linear or linearized problem. We give some numerical examples showing comparisons with the EBE method. In section 6 we give some real applications comparing our method and the sparse matrix method which used only nonzero entries.

Storage allowed.

The total storage required for H, B and I is : In the frame of an Element By Element approach, the saving of storage increases with h and p refinement. Under the same assumption that all the needed operators are stored, The ratio p s between the storages needed by EBE and UDA approach is : neln\ In the case of a same p refinement (npi e = npi, n e = n Ve) we have : n 2 01 npi n an n + 02 npi n [START_REF] Axelsson | An efficient finite element for nonlinear diffusion problems[END_REF] an + npi By using an exact integration rule for all the finite elements we have, for plane (n an = 3) and 3-D problems (n an = 6) : For EBE or for the classical approach using K the cost of storage does not depend on the number of integration points. So this cost remains the same if a reduced integration is used and thus whatever the methods employed. As we can see by considering equation (20), St is a function of the size of the various operators and of parameters 0i and 02. Therefore, a reduced integration leads to an important reduction of storage. The magnification m of p s , which is a function of the number of reduced integration point npi r is equal to :

S t =
m = npi npi r (21)
The reduced integration rule is based on the order of the integrated polynomial. If npi is the quadrature rule used to integrate a polynomial function of degree n then npi r is the one which integrates a polynomial function of order n -1 As we can see in tables 3 the magnification m of p s is significant. This parameter is also equal to the ratio between the stress (or strain) size of the problems obtained by using an exact integration and that of problems corresponding to a reduced integration. This reduced integration leads then to an increasing speed of computation of the resolution process. This is not a feature of the classical approach.

Computation time.

The matrix-vector products in a EBE method is worth n 2 e for one element. In UDA method, the cost matrix-vector products for one element is equal to the sum of : This subsection dealt with the reduction of the computation time induced by the consideration of operators H, B and I instead of K or K e , but not with the computation time required by the solver. These advantages relating to this last point are to be added to those described above, and to the post-processing step.

B e • U e =2
For a non-linear problem where the non-linearities are induced by the material behaviour, the only non-linear operator is the constitutive one H, since we assume small transformations. This leads to a unique computation of operators B and I for the global analysis of the mechanical problem.

Thus, with classical EBE or with K a non-linear analysis leads to update all of them at ecah iteration. Theses methods do not benefit of the fact that the matrices B and I remain constant, because the main operators K or K e have to be formed for each iteration of the linearization process.

6.

THE ITERATIVE ALGORITHMS

The various algorithms presented below can be used equally well for displacement one-field finite element analyses or for real stress-displacement two-fields mixed ones. For this last problems, we have only to substitute A to IH" 1 and B to ÎB, where A and B are defined by the following relations.

nel " Lets us remind the reader that the problem modelled together with its equivalent form obtained after linearization at a given loading step are as follows:

^--v

\ classical

Find UeJR N such that :

K-U = F
Find U £ M n such that :

(22)
6.1 GMRES method.

Let ||.|| and (.,.) denote respectively the Euclidian norm and inner product. The iterative method used for the resolution of the non-symmetric algebraic problem ( 22) is the GMRES one [START_REF] Saad | GMRES: a Generalized Minimal Residual algorithm for solving nonsymmetric linear systems[END_REF] . This numerical method, which is a particular case of the Petrov-Galerkin method, is based upon a minimization of the Euclidian norm of the residual vectors r m , in a Krylov space K m of dimension TO,

K m = /C m (K) = span {r 0 , K • r 0 , K 2 • r 0 , • • •,
• r 0 }.

Let us consider an approximate solution of the form XQ + z m , where XQ is an initial guess and z m belongs to the space K, m . Then the above minimization problem can be expressed as follows [START_REF] Saad | GMRES: a Generalized Minimal Residual algorithm for solving nonsymmetric linear systems[END_REF] :

min \\F -K • (x 0 + z m )\\ = min [|r 0 -K • z v £fc zÇK ( 23
)
This algorithm uses an orthonormal basis

V m = [vi,V2, • • • , v m ] of K, m
which is built with the classical Arnoldi [START_REF] Saad | GMRES: a Generalized Minimal Residual algorithm for solving nonsymmetric linear systems[END_REF] procedure. If we take as the operator the stiffness matrix K , the restarted GMRES (TO) flowchart is given, in a synthetic form, by : Initialization :

VQ initial guess [f] r 0 = F -K • «o Building of the orthonormal basis V m of /C h i,j = (K-Vj,vi) i = {l,2,---,j} } j \vj j = {1,2, i+i|l "i+ij
Resolution of the reduced minimization problem :

v m = vo + V m • y, with y G IR m minimizes ||e -iî r
Convergence check :

||r m = F -K • V m \\ < £ (userspecified tolerance)
If convergence check is not satisfied then affect v m to vo and go to f.

In the algorithm above H m is the (TO + 1) x (m) upper Hessenberg matrix :

i,m 0 INI L o o \\Vm 0
and e = *(lko||, 0, • • -, 0) G IR m+1 is a basis vector. Let us denote this algorithm as GMRES-K. We shall now proceed to the substitution of K by the operators described in the previous sections.

Substitution.

The algorithm given above is written for an approach based on the stiffness matrix. If we now apply this algorithm to the dissociated form (I), it will be called GMRES-UDA and is as follows

Flowchart of GMRES-UDA 1 VQ initial guess V "' " «roll L/-i_L1 1 * J v m = v o + V m • Convergence check If not satisfied : VQ IB-Vj^i) HB-Vj- y, with y llr -F - II' m -± = '• V m g°0 3 -J2 hi i-l -*BÎ to [f] {1,2,--.J} jV{ minimizes ||e -H m -y\\

HB -u TO || < e

In the algorithm above, we have not to form the stiffness matrix, so that all the operators remain distinct at each step of the GMRES iterations. The algorithm GMRES-UDA is more efficient than those based on the stiffness matrix such as GMRES-K. Let us note that when the GMRES-UDA has converged in U = v m , we also get the strain vector e = B • U and the stress vector çr = H • e, since they are intermediate computed values required by the construction of V m . We only need two auxiliary arrays to store the stress and strain vector, however their size a rather small compared to that of other arrays.

Numerical results.

In order to demonstrate the performance of the method proposed we apply the algorithm described above on two non-linear geomechanical problems. To check the quality of the proposed method we calculate the parameter p defined by:

t[ N 2 f 2 -
where N\ is the number of linear or linearized algebraic problems, and iV 2 the number of one-GMRES(m = 1) iterations for the whole problem. The quantities t\ and t 2 are respectively the times allowed to one formation of K and GMRES(TO = 1) iteration, whereas t[ and t' 2 have the same meaning in the case of UDA.

The first problem modelled is a loading of a rigid circular fondation on an elastoplastic soil (fig. 1). From left to right, figure 1 showns the geometry of the problem, the field of vertical displacement and the horizontal one. The soil behaviour is modelled by a non-linear constitutive model [START_REF] Darve | Incrementally non-linear constitutive relationships[END_REF] .

fig. 1 : Problem of the rigid circular fondation.

The number of elements is equal to 720 quadratic (Te) or cubic (Tio) triangular elements. In the table 6 we recapitulate the principal features of the problem : the number of elements ne/,the number of degrees of freedom N and the number of nonzero entries (NZK) in K used in a sparse matrix method. The second problem modelled is a soil slope loaded on its head (fig. 5). From left to right on figure 4 we have shown the geometry of the problem, the field of vertical displacement and the horizontal one. The number of elements is again equal to 720 quadratic (T 6 ) or cubic (Tio) triangular elements.

In the table 8 we have recapitulated the main features of the problem.

Element T 6 Tio

JV 3078 6776 NZK 128029 429497 Table 7. Features of the algebraic plane strain problem, and the results obtained with different integration rules are shown of figures 5 and 6. As we can see the values of ratio p are significant. They increase with the number of elements and the space-dimension of the problem. We shall now focus on the way to modify the original displacement algorithm GMRES-UDA in terms of strain or strain rate variables.

6.4. The strain approach.

In the previous section we have shown that the dissociated form (I) of the displacement problem is more interesting than the classical one using K. In this section we shall show how to increase these numerical performances. It is well known that the displacement field does not belong to the set of memory parameters of the constitutive equations. This feature is very interesting in the non-linear case. And indeed, we don't need to compute the discrete displacement field of the linearized problems arising from some methods like the Newton's one.

If one is not interested in the displacement fields of these mechanical problems, one can do without computing them. If one is interested in these fields for some time values (not necessarily for all the steps of loading) one can easily recover them. Eventually, if one wants to build the displacement fields for all the steps of loading, one can do it easily without performing useless computations at each iteration of the GMRES iterative process. These features are offered by the algorithm described below.

Considering the calculation of the coefficient inner product in the algorithm :

= ||*B î H B • VJ|| 2 + % • f Vj Vj • hj -2 t hj-t Vj t BÏti.B-v 3 = ||*B î H B • VJ || 2 + % • h 3 -2 \j • V|| 2 %•£
We then obtain the relation :

ll^+ijll 2 = Pi+iH 2 = || f BÎHB -vj\\ 2 -J2 h h j
and the vector ê^+ 1 is defined by :

S +1 = B • v j+1 = - B • ù i+1
After all theses algebraic transformations, the new problem does not consist in the research of the displacement vector U, but in that of the associated strain vector e. In other terms

v m -VQ + V m • y, with y G H m minimizes ||e -H m • y\\ is replaced by êvm = tvo + [ê]m • V-, with y G H m minimizes ||e -H m • y||,
The new algorithm expressed in term of strain can be written as follows :

Flowchart of GMRES-e hi hi

hj.

am -

initial guess IMI .1 = 3*3 n,3 = f II'B îo + [ê] m • y, Convergence check | If not satisfied: CQ = 1, :) ÎH- ÎH j i ê 3 • ê 3 with r m = -F go = {1, 3 ||2 _ ' % 2,-- i =i yeiR m -*B to[f] ÎH 0" •i = {l,2,-,m} minimizes ||e -iï m • y|| • Ul < E with [ê] m = [êi, ê 2 , • • •, ê m ] and ê; = ê Vs
The difference between GMRES-UDA and GMRES-e lies in the number of B-products which constitute the most expensive operations. In GMRES algorithm the most important computation time is that allowed to the construction of the orthonormal basis V m of the Krylov space K m .

A substantial reduction of this computation time leads naturally to a significant improvement of this algorithm. GMRES-UDA is more paying than GMRES-K. However, GMRES-e is more paying than GMRES-UDA. And indeed the number N of B-products avoided is equal to : A^ = m g x n is (25) where nj g is the number of GMRES iterations, tin the number of linearized problems (in the non-linear case) and n ts the number of loading steps of the problem. For example, the table 8 give the contribution of B product in percent in each K X X product. Eventually let us note that when the GMRES-e have converged in e = e m we obtain also the stress vector g_ = H • e, since it is an intermediate computed variable for the building of matrix V m .

If one wants to recover the displacement vector U at some given loading step, one has to solve, out of the GMRES loop and of course out of the non-linear loop, the IR linear problem :

Find U e M N such that : (26) BU = e

Preconditioning purpose

This short section is devoted to some recommendations which can be taken into account when preconditioning of the linear or linearized problem is performed. Since the problem of an optimal preconditioning remains open, we give only some indications.

A strong restriction of the preconditioned methods is that they require to maintain the special dissociated or factorized (following Axelsson et al. [START_REF] Axelsson | On the updating and assembly of the hessian matrix in finite element methods[END_REF] ' 2 ) form, at each inner (for the linearized algebraic problem) and outer (for the linearization procedure) iteration of the solving process.

The question of preconditioner is natural within iterative solving framework. They are performed in order to enhance the robustness of the iterative solvers.

The simplest preconditioned method is the diagonal one : An other direction which have to be investigated, lie in the interesting spectral equilavence theorem established by Axelsson et al [START_REF] Axelsson | On the updating and assembly of the hessian matrix in finite element methods[END_REF] . Following Axelsson 1 ' 2 two matrices M and N of dimension TV are called spectrally equivalent if :

3c x > 0, c 2 > 0 such that Cl *X • M • X < f X -N • X < c 2 l X • M -X VX G IRÂ
xelsson and Gustafsson 2 used successfully this method in a non-linear diffusion problems.

In Krylov subspace framework, more specially in GMRES subspace iterative methods, one could used the so call FGMRES referred as flexibles GMRES 23 . In this framework we dont used a unique preconditioner for the linear algebraic problem, but it perform variation with subspace iteration in order to improve the robustness of the preconditioner. Saad 23 propose to incorporate the changes in the precondtioner into the classical GMRES algorithm.

All these methods will be developped and investigated.

CONCLUSION

In this paper we have shown how to increase substantially the performances of the Gmres iterative algorithm when used for the numerical resolution of algebraic systems coming from the f.e.m. The theorem established shows that from an algebraic point of view the displacement one-field approach can always be considered as the dual one with subsidiary constraints. This allows us to use algorithms which can not be employed for classical algebraic systems based on the global or elementary stiffness matrices. The numerical performances of two iterative algorithms were described and proved. One has to note that these algorithms can be easily implemented in a finite element code without requiring substantial modifications of its structure. The B-bar method [START_REF] Golub | Résolution numérique des grands systèmes linéaires[END_REF] can be also improved thanks to our two algorithms GMRES-UDA and GMRES-e. A spectral analysis can be performed by considering the dissociated form (I) of the problem, since the structure of the corresponding algebraic systems is more interesting than that of systems built with the stiffness matrix. More developments will be done, especially as concerns the coupling between the displacement one-field approach and the stress-displacement two-fields mixed one as well as the extension of the previous algorithms made to include large transformations. This last extension can be done easily without adding geometrical non-linearities to the problems modelled if one adopts a rate-type point of view involving both objective stress-rate and velocity fields. [START_REF] Royis | Formulation mathématique de loi de comportement-Modélisation numérique de problèmes aux limites en mécanique des solides déformables[END_REF] ' [START_REF] Royis | Modélisation par éléments finis des géomatériaux[END_REF] The method will be also applied in adaptive refinement framework.

  net nel K = J2 K e, K e G M n « xn « and F = ^Fe , F e G E" £ e=l ^ e=l

  " 1 :a) : sdÇî -[s: e(v) dQ = 0 Vs G S 1 Jo, / a : e(w) dQ, = I f • w dO. + I g • w dT Vu; G V (18)

  0 2 y^ npi e n 2 an + 0\ Y^ npi e n an n e + V^ npi e e Lets us call p s the ratio between the storage required by K and that required by B, H and I. storage required by K storage required by B , H and I

If iA

  = > / -0. H TO. dÇl _ Ja e nel " B = V / V e B e dû ^ In e in which ip e denotes the interpolation functions of stress : <T(X) = 1p e (x) • <7 e , CT e : elementary vector of stress components.

Table 6 .fig. 2 :

 62 fig. 2 : Values of p for the quadratic element, npi = 3(a) npi =

1 fig. 4 :

 14 fig.4: Soil slope loaded on its head.

fig. 5 :

 5 fig. 5 : Values of p for the quadratic element, npi = 3(a) npi = 1(6).

  hij = (K • Vj, Vi) ij of the Hessenberg matrix we can replace the by a relation which shows off the vector of strain components hi.i -€" • e,. € Vi -B • According to the orthonormality of vectors V{ G M N i G {1, 2, • • -,j] we can rewrite the expression of the components of the Hessenberg matrix /ÏJ+IJ. j \\h j+ i tj \\ 2 = \\v j+ i\\ 2 = I^BÎHB -VJ -Yj h i,j v i\? (24) Let us define t hj = (hij, h,2j, • • • , hjj) and the basis Vj = [v-y,V2, • • • , Vj]. Taking account ofthis notation the expression (24) can be written as :

D

  = diag [*BÎHÉ]which tansform the problem,(t BÎÊ.B) -U = F => B 1 / 2 ( t BÎÈ.B]D 1/2 -Û = F with Û = D 1 / 2 • U and F -D 1 / 2• F. This simple method do not need to form explicitly the stiffness matrices but require only the computing of the iV diagonal entries Da.

Table 3 .

 3 Magnification m of p s for some 2-D and 3-D elements.

	. For instance we have, for the finite element of

Table 4 .

 4 npi e n eThus equal to npi e (4 n e + 9) for two dimensional problems and to npi e (6 n e + 36) for 3-D problems. Let us call p t this cost obtain with an exact integration rule and p tr the one obtained with reduced integration rule. We give in tables 4 and 5 the values of n 2 e) p t and pt r for the most used finite elements. Values of nl, p t and p ir for some 2-dimensional elements.

	He€ R	=	°=

Table 5 .

 5 Values of n%, p t and p tr for some 3-dimensional elements.

Table 8 .

 8 Contribution of B product in each K X X product.

	2D-Finite element	B x X (%)	3D-Finite element	B X X (%)
	Bilinear Lagrange	39.0	Trilinear Lagrange	40.0
	Linear triangle	36.4	Linear tetrahedron	33.3
	Biquadratic Lagrange	44.4	Triquadratic Lagrange	46.6
	Quadratic triangle	42.1	Quadratic tetrahedron	44.4
	Bicubic Lagrange	46.7	Tricubic Lagrange	48.5
	Cubic triangle	44.9	Cubic tetrahedron	44.4


O Deck

M Al Heib

F Homand

The take into account of soil-structure interaction phenomena to study mining subsidence consequences on buildings. Prise en compte des phénomènes d'interaction sol-structure pour l'étude des conséquences des affaissements miniers sur le bâti

The authors used a numerical finite element software for several modellings to study soil-structure interaction phenomena which occur during mining subsidences. We broke down the ground movements into several basic movements to highlight the impact and relative importance of each other (ground curvature and horizontal strain). We elaborate synthetic schemes of the ground and structure behaviour which improve the traditional apprehension of phenomena.

RESUME -Les auteurs ont utilisé un code de calcul aux éléments finis pour étudier les phénomènes d'interaction sol-structure qui ont lieu durant les affaissements miniers. Nous avons distingué les différents mouvements des terrains (déformation horizontale et courbure) afin de mettre en évidence leur impacte respective. Nous avons élaboré des schémas synthétiques, du comportement des terrains et des structures, qui permettent d'améliorer la conception traditionnelle des phénomènes en place. Compressive strain : e=-5.10" 3 Extensive strain : E=-5.10' 3 -a/4 Convex :R =1OOOm Concave: R =-1000m =-Failure line Confined area with no failure Horizontal displacements along vertical sections © Horizontal stresses in structure, traction ® (+) or compression (-) Vertical stresses along horizontal sections

Introduction

The underground mining of raw materials may be the cause of movements on the ground surface. Whether planned or accidental, such movements can cause considerable damage to structures located within the area of influence of underground mining work. This is borne out by the recent mining subsidences that took place at the end of the 1990s in the Lorraine iron-ore field. A better understanding of how ground surface movements can be imparted to the supporting structure and damage it, has now become necessary. Indeed, it is too often considered that damage depends only on ground strain with no account being taken of soil-structure interaction phenomena which may affect structural behaviour considerably.

We will first describe the ground movements caused by a mining cave-in. Then we will investigate the relative effect of each ground movement on a rigid structure, which behaviour is complex because of strong interaction between the ground and the structure.

Description of the phenomenon

Mining subsidences produce considerable horizontal and vertical movements on the surface of soils (figure 1). It is usual to consider, on the one hand, the vertical subsidence, the maximum value "Sm" of which is considered as a characteristic of the trough, and, on the other, the horizontal strain of the soil, its curvature and its slope which are the three movements loading structure and causing structural damage.

Methods for predicting structural damage generally use the definition of threshold values of strain and curvature from which such damage is expected. Some authors offer some behaviour schemes providing an overall view of the structure behaviour under the effect of each ground movement. A few examples of behaviour in a curvature area and in a strain area are shown in the figure 2 and figure 3. The outcome of these few examples is that the fundamental question is whether soil movements are imparted integrally, partially or not at all. "Integrally" means that the structure is loaded by displacements, "not at all" or "partially" indicates the existence of complex soil-structure interaction phenomena.

Investigation of basic loading

The aim of this investigation is to highlight the overall behaviour of a solid-wall type structure, i.e. very rigid compared with the ground. For this purpose, we carried out a set of numerical modellings in order to investigate the behaviour of a structure with regard to each basic loading: horizontal strain and curvature. The effect of the slope was not investigated since the slope, as a rigid body movement of the structure, produces very few stresses in the structure.

We used the finite element software called CESAR-LCPC. Elements are six-node triangular, eight node quadrangular and six-node interface elements. The latter one make it possible to take into account a friction type behaviour with an interface separation possibility between the ground and the structure (zero cohesion and tensile strength, friction angle 20° and dilation angle 0°). In all the modellings, the ground behaviour was considered to be elastoplastic with a Mohr-Coulomb failure criteria (zero cohesion, friction angle 30° and dilation angle 10°). The structure was considered to be elastic. Thus, results which we are going to present here correspond to the structure and soil behaviours before structure damage. This choice is entirely legitimate given the very high mechanical strengths of the structure compared with those of the ground.

Curvature

Vertical section in the structure Horizontal section in the ground Figure 4 : Description of the models used to investigate the effect of the ground curvature and horizontal strain.

Two sets of models were produced, and are shown in figure 4. The first model is used to investigate the effect of a curvature of the ground. The second model is used to investigate the effect of a horizontal strain of the ground. The two models briefly represent a ground massive measuring 15 m by 15 m with a solid structure on its surface. The left border represents a symmetry axis. The calculations were carried out for plane strains. To model the effect of the ground curvature, we modelled the ground feature with a surface curvature (concave or convex curve) on which we modelled a perfectly horizontal structure. Between the structure and the ground, we placed interface elements which are initially detached except at the ends of the structure for the concave area and in the middle for the convex area. The ground was consolidated under its own weight, then we loaded it by the weight of the structure. The results were exploited by analysing stresses in the ground, in the structure and along privileged sections shown in figure 4. In addition, particular attention was paid to the localization of the failure points in the model. To model the effect of the horizontal strain of the ground (tensile or compression), the ground feature, whose surface is horizontal, was consolidated under its own weight and then loaded by the weight of the structure and finally deformed by imposing a uniform horizontal displacement on the right border of the model. Various curvatures (+/-1000 m to +/-7500m) and strains (+/-0.5 10" 3 to +/-5 10" 3 ) were investigated as well as various Young modulus of the ground (30 to 100 MPa and Poisson's ratio 0.25) and of the structure (1500 to 15000 MPa and Poisson's ratio 0.2).

Curvature investigation

The effect of a curve of the ground on a structure is easy to understand. Figure 2 shows different schemes of the behaviour in a curvature area. In a concave area, we expect to observe a concentration of vertical stresses in the ground at the edges of the structure and in the convex area under the middle. We can compare these schemes with results of numerical modellings in figure 5 which show results in the whole model. Large changes can be observed of the plastified area, which differs from the area obtained in the reference investigation. The plastic points are concentrated perpendicularly to the areas on which the structure rests preferentially. They are thus deduced by the elevation of the vertical stresses.

The representation of the field of the main stresses highlights remarkably the redistribution of the stresses in the ground. The field observed in the structure revealed bending behaviour consistent with the ground curvature. However, the structure curvature is under no circumstances equal to the ground curvature. The structure stiffness is such that it retains a horizontal geometry. Tensile stresses appear in the lower fibre (0.1 MPa) in convex area and in the upper fibre in concave area.

The analysis of the vertical stresses in the ground under the structure shows that the disturbances are greater for a concave curvature than for a convex curvature. The vertical stresses in the ground, where the structure rests preferentially, increase by 60% in convex area and by about 100% in concave area. In addition, a structure is more likely to be found partially detached from the ground in the concave area than in the convex area. We can thus confirm the intuitive approach adopted by Kratzsch (1983), who suggested that the structure remains horizontal, compared to that of Boscardin and Cording (1989), who suggested that the curvature of the ground is transmitted to the structure. The study showed that the possible separation between the ground and the structure is not related to the bending flexibility of the structure but instead to the vertical flexibility of the ground.

Investigation of strains

In the same way as for curvature, the effect of horizontal strain of the ground on a structure has often been grasped qualitatively (figure 3). The main idea adopted has been to assume that the horizontal strain of the ground only produces a uniform horizontal strain in the structure. The investigation presented here (figure 6) provides a visualization of the behaviour of a structure under such a movement. The failure area in the model vary differently according to the strain direction. A generalised failure in the model can be observed in the tensile area, which can be readily understood given that the tensile strain tends to decrease horizontal minor stresses whereas vertical major stresses remain the same. In compression area, however, the presence of a confined area under the structure can be observed, similar to that obtained with the reference investigation (structure resting on a horizontal ground without strain) and which is similar to the cone of ground which remains integral under a vertically loaded foundation.

The stress field in the model shows that a compression strain generates a concentration of stresses in the ground under the middle of the structure as well as bending stresses in the structure. In a tensile area, however, the stresses seem to be relatively little modified in the ground or in the structure.

The variation in the vertical stresses under the structure provides a fuller explanation of the overall behaviour observed. The compression strain leads the structure to rest preferentially on its middle section. The increase in the stresses is then about 40%, i.e. an amplitude comparable with that observed for the curvature. In tensile area, however, a re-homogenization of the vertical stresses under the structure can be observed.

The horizontal stresses in the ground increase only slightly in tensile area. In compression area, however, a penetration of the horizontal stresses under the structure can be observed. No increase of stresses is observed towards the edge of the structure. Indeed, the presence of the free surface of the ground, does not allow a significant increase in the horizontal stresses because the vertical stresses are very low. It is of interest to compare this result against one of the remedial actions proposed to minimise the effect of the strains: the digging of trenches around the periphery of the structure is often cited as a technical solution. But horizontal stresses observed on the ground surface are low because of the failure criteria adopted (Mohr-Coulomb). So, such a solution does not completely cancel the effect of the horizontal strain of the ground. A shearing phenomenon appears naturally under the structure. Contrary to what is often assumed (Kratzsch, 1983), the effect of shear is not constant. In a compression area, the variation in the shear stresses is more complex. This phenomenon is associated with the distribution of the vertical stresses which vary significantly along the structure.

The structure is affected by horizontal strains of the ground much less than if such strains had been integrally transmitted. Indeed, a strain of +/-5.10' 3 would then have the effect of a horizontal stress of 25 MPa in the structure. Instead of such a value, the compression area causes a compression stress of -0.2 MPa in the lower fibre and a tensile stress of 0.1 MPa in the upper fibre. The structure is thus bent. To a less extent, the same phenomenon can be observed in a tensile area since the latter generates tensile stresses of 0.04 MPa in the lower fibre and -0.02 MPa in the upper fibre. Compression stresses distribution in the soil and the direction of shear stresses are both responsible for the bending observed.

Synthesis

Using the previous set of results obtained, added with those of parametrical study, it is possible to draw behaviour schemes of the structure in a curvature area or in a strain area. These models are shown in figure 7. These figures show the location of slippage lines and the distribution of vertical stresses in the ground and under the structure at various depths. These diagrams must be compared in figures 2 and 3. A number of elements are common. In particular, the slippage lines proposed by Speck and Bruhn (1995) or Kratzsch (1983) are fairly similar, which thus consolidates our results since, for the latter, the slippage lines had been observed experimentally on scale models in sand. The horizontal stresses observed along the middle section of the structure were also highlighted. The stresses are shown as a function of a parameter "a" for the single purpose of being able to make comparisons between them. These values correspond to the investigation presented in detail (E grO und=100 MPa, E struc , ure =5000 MPa).

Figure 7: Behaviour of a structure with regard to horizontal ground strain and ground curvature.

Conclusions

The purpose of the investigation presented was to study the soil-structure interaction phenomena at work during a mining subsidence. Structures concerned are loaded, among other things, by a ground curvature and horizontal strain. Because of complex phenomena, we decided to dissociate the ground curvature from strain. The analysis provide a better identification of the various behaviours of grounds and structures under the effect of each ground movement. It remains for us to investigate the combination of these loadings. The main results of our investigation are outlined below: Soil-structure phenomena may be very important during subsidence so that buildings movements may be real lower and different from ground displacements. The investigation of the ground curvature is consistent with the intuitive behaviour schemes generally adopted. The strain investigation shows that the behaviour is more complex than assumed until now. In particular, strain causes a bending moment in the structure which is not negligible compared with the one induced in the curvature area. The effect of tensile strain is also to produce immediate failures in the ground which prevent stresses from being imparted to the structure. The numerical investigation made it possible to compare the effect of the various loadings with respect to each other and the effect of the ground mechanical properties. We propose overall behaviour schemes of a structure in a curvature area and in a strain area. The latter provides a better understanding of the phenomenon and a better assessment of the intrinsic effectiveness of the preventive measures to protect buildings. The overall behaviour schemes proposed in a tensile area would suggest that the structure stability cannot be guaranteed since the ground is no longer secured. Such instability cannot, however, be understood by means of the numerical software used.

Boscardin M.D., Cording E.J. (1989) Building response to excavation-induced settlement. Journal of Geotechnical Engineering, Vol. 115, No, 1. Kratzsch 1-1.(1983) Building standards in West Germany. Mining Subsidence Engineering, Springer-Verlag. Speck R.C. et Bruhn R.W. (1995) Non-Uniform Subsidence Ground Movement and Resulting Surface-Structure Damage. Environmental and Engineering Geoscience,Vol. 1,No. 1,

Figure 3 :

 3 Figure 1 : Description of surface ground movements as a result of a mining subsidence.

Figure 5 :

 5 Figure 5: Stress fields and failure areas in a curvature area (+/-1000m)

Figure 6 :

 6 Figure 6: Stress fields and failure areas in a strain area (+/-5.10" 3 )


S Brohez

C Delvosalle

G Marlair

A Carrau

A zone model for the prediction of thermal and chemical effects of a pool fire in a forced ventilated enclosure

Zone models have been developed since the early 60s for the prediction of fire parameters such as smoke temperature, smoke filling and movement in multi-compartment buildings. Unfortunately, one major difficulty in current zone models is that heat and chemical species release rates in relation with a given fire source term are usually to be provided as input data. A new zone model for prediction of thermal and chemical effects of pool fire in a forced ventilated enclosure has been developed. The novelty of the approach relies in particular on the provision of three sub-models that are used for reducing the number of input data needed for a given simulation. The burning rate history of liquid pool fire is calculated from an original submodel of pyrolysis. A solid flame sub-model is used for predicting radiative properties of flame. Yields of chemical species are estimated from a dedicated sub-model of combustion. This zone model has received some validation for use in forced ventilated enclosures only. Promising results have been obtained.

INTRODUCTION

A large fire in a compartment is often considered as one of the most hazardous accidental event which may affect safety in industries. The fire damages may be thermal or non thermal. For examples, intense radiation produced by large fires may cause serious burn injuries to the staff of the industrial premises and the fire fighters. Moreover, the fire plume may transport a variety of toxic pollutants which can be very dangerous for people, and polluted extinction waters while unconfined may greatly affect the aquatic environment.

Zone models have been developed since the early 60s for predicting fire parameters such as smoke temperature, smoke filling and movement in multi-compartment buildings [START_REF] Janssens | Room Fires Models. General[END_REF]. In this zone modelling approach, the space within each compartment is generally divided into one or two control volumes [START_REF] Cox | Compartment Fire Modelling[END_REF][START_REF] Jones | A Technical Reference for CFAST: An Engineering Tool for Estimating Fire and Smoke Transport[END_REF][START_REF] Quintiere | Fundamentals of Enclosure Fire Zone Models[END_REF]). Following experimental observations of thermal stratification of gases inside the compartment (pre-flashover fires), two-zone models divide the room gas volume into two distinct and uniform layers: an upper layer of hot gases and a lower layer of air. Inside each of those layers, gas-phase physical parameters such as temperature and species concentrations are assumed to be uniform. On the other hand (post-flashover fires), one-zone models consider a single layer in the compartment, gases inside the compartment are supposed well-mixed. The physical parameters (temperature and composition of gas) of the layers are predicted by solving mass and energy conservation equations on each layer.

Unfortunately, zone models still have important limitations. One major difficulty in current zone models is that heat and chemical species release rates in relation with a given fire source term have usually to be provided as input data. Consequently, there is no feedback interaction between ventilation conditions and these input parameters, which could be assumed unrealistic in many practical cases since ventilation conditions are generally unknown prior to running the simulation. In order to avoid the mentioned drawbacks, a new zone model has been developed for simulating the behaviour of liquid pool fires in enclosures. The novelty of the approach relies in particular on the provision of three sub-models that are used for reducing the number of input data needed for a given simulation. The burning rate history of liquid pool fire is calculated from an original sub-model of pyrolysis. A solid flame sub-model is used for predicting radiative properties of flame. Yields of chemical species are estimated from a dedicated submodel of combustion. At present time, this zone model has received some validation for use in forced ventilated enclosures only.

CONTROL VOLUMES OF THE ZONE MODEL

In the conditions of the experiments conducted for validation purposes of the proposed model (in an enclosure with forced ventilation), fires carried out with low ventilation rates were observed to produce single layer environment in terms of chemical species concentrations (Backovsky et al 1988[START_REF] Peatross | Ventilation Effects on Compartment Fire Characterization[END_REF]. Accordingly, in the zone model, the gas phase inside the enclosure is thus supposed well stirred.

Moreover, in classical zone modelling approach, the point source thermal model is considered for the flame, the energy radiated by the flame is specified by the user and is a fraction of the heat release rate (typically 30%). In the development of this zone model, the classical solid flame radiation model initially developed for fire in open space [START_REF] Mudan | Thermal Radiation Hazards from Hydrocarbon Pool Fires[END_REF]) is considered for modelling the flame in the compartment. Accordingly, there are two control volumes: one for the flame and one for the fumes (see fig. 1).

Figure 1. Control volumes selected in zone modelling 3 PYROLYSIS SUB-MODEL One major difficulty in classical zone models is that mass burning rate (or heat release rate) in relation with a given fire source term is usually to be provided as input data. As this zone fire model is dedicated to simulate liquid pool fires in enclosures, the mass burning rate is predicted from correlation of Zabetakis & Burgess (1961) valid in open space:

m" f =10 -3 AH,

0)

vap where m" f = mass burning rate (kg.m" 2 .s"'); AH comb = heat of combustion ; and AH vap = heat of vaporisation at the boiling point of the liquid fuel. However, in the configuration of the one-zone model, the air entrained into the fire is partially depleted in oxygen. The reduced oxygen concentrations at the base of flame leads to a reduced mass burning rate [START_REF] Peatross | Ventilation Effects on Compartment Fire Characterization[END_REF][START_REF] Nikitin | Variations of Mass Combustion Rate with Oxygen Concentrations and Gas Pressure of a Milieu[END_REF] . The influence of air partly depleted in Q 2 on the mass burning rate is introduced additionally by use of the correlation recently proposed by [START_REF] Peatross | Ventilation Effects on Compartment Fire Characterization[END_REF]:

m. m f,21%O 2 (2) 4 COMBUSTION (CHEMISTRY) SUB-MODEL
In classical zone models, the yields of chemical species in relation with a given fire source term are usually to be provided as input data. In this zone model, the combustion sub-model is based on the concept of global equivalence ratio Q [START_REF] Pitts | The Global Equivalence Ratio Concept and the Formation Mechanisms of Carbon Monoxide in Enclosures Fires[END_REF].

In order to develop the combustion sub-model, experiments were carried out at lab-scale on the Fire Propagation Apparatus operated by INEPJS (apparatus in agreement with ASTM E2058) with different mass flow rates of incoming air. For several chemical substances (pyridine, adiponitrile, 1-chlorobutane, thiophène,...), the yields of CO 2 , CO, total unburned hydrocarbons, soot, NO, HCN, SO 2 , and O 2 were measured as function of the global equivalence ratio [START_REF] Marlair | Use of ASTM E2058 fire propagation apparatus for the evaluation of underventilated fires[END_REF]. As examples, yields of CO 2 and CO are presented in figure 2, as a function of the equivalence ratio for the solvent pyridine. Following the data processing procedure proposed by [START_REF] Tewarson | Ventilation-Controlled Combustion of Polymers[END_REF], the ventilation controlled fire properties are related to those applying for well-ventilated conditions using the following generalised relationship:

fijb.v. = 1 + a (2)
where,,/; = fire property valid for any ventilation rate determined by the equivalence ratio value; f i:W . v . = the same fire property applying for well ventilated fires; and a, p and £ = experimental correlation parameter.

The user has to choose a fuel from a chemical substances database file. This file contains empirical correlation parameters (a, P and Ç) of 8 chemicals species (CO 2 , CO, total unburned hydrocarbons, soot, NO, HCN, SO 2 , and O 2 ) as a function of the global equivalence ratio. The yields of N 2 is obtained from the conservation of nitrogen by considering the chemical formulae of the fuel (

C a H b O c N d X e S f
) as well as the NO and the HCN yields:

JN 3 _ d /NO fncN Q\ M N 2 M comb M NO M HCN
The yields of H 2 O is obtained in a similar way as [START_REF] Beyler | Major Species Production by Diffusion Flames in a Two-Layer Compartment Fire Environment[END_REF].

It can be pointed out that the database file also contains fuel properties such as heat of vaporisation, specific heat and boiling point. The chemical substances database file today contains fuels such as heptane, adiponitrile, pyridine, TDI, isoproturon, 1-chlorobutane and thiophene.

SOLED FLAMME SUB-MODEL

In the solid flame sub-model, the flame is modelled as an opaque body. The user only specifies the flame diameter. The flame height is given by the well-known Thomas' correlation (1963). Air entrainment into the burning region of pool fire is estimated from Delichatsios' work and related correlation [START_REF] Delichatsios | Air Entrainment into Buoyant Jet Flames and Pool Fires[END_REF] established in open space.

MASS AND ENERGY CONSERVATION EQUATIONS

Fumes control volume

Mass flow rate, chemical species concentrations and temperature of smoke are calculated respectively from mass, chemical species and energy balances on a control volume for gas (smoke) inside the enclosure. The conservation of mass in the fumes control volume leads to the following equation:

, 'fum-Pfim M fum m air + m pan ~m fum ~ m ent = ~7 W
where m k = molar flow rate of flux k; V fum , p fum and M fum are respectively the volume, the density an the molecular weight of fumes.

It is convenient to considered all gaseous species as ideal gases. An ideal gas is any gas that obeys the following ideal gas equation:

(5) where P = the pressure; and R = the universal gas constant.

As the pressure inside the compartment can be supposed constant [START_REF] Zukoski | Development of a stratified Ceiling Layer in the Early Stages of a Closed-Room Fires[END_REF], then substituting p fum in equation 4 by equation ( 5) leads to the following equation ( 6) :

fum ffum m air + m pan ~ m fum ~ m ent R (6)
The conservation of species i in the fumes leads to the following equation:

Kir yi,air + ' + ,) = fum • yi. l fum R dt (7)
where y t k = molar fraction of species i in flux k.

The conservation of energy in the fumes leads to the following equation:

fh f ) ™air Kir ~ (m fum + th ent )h fum + q fum =- ' l fum ffum R dt ( 8
)
where h k = molar enthalpy of flux k ; and q fum = net heat transfer to the fumes (W)

Flame control volume

Mass flow rate and chemical species concentrations of plume are calculated respectively from mass and chemical species balances on a control volume for flame. It is assumed that time constant of flame is small in comparison with that of fumes. Hence, a quasi-steady state is assumed for flame (the volume, composition and temperature of flame are assumed to change slowly over the duration of the time step). The conservation of mass in flame control volume leads to the following equation:

N+lf M ™ent + ™fuel S fi i=l \ fuel M, =0 (9)
where N+l = number of species ; and,/* = chemical yields of species / (g of z.g" 1 of fuel)

The conservation of species i in the flame leads to the following equation:

m ent s m fuel m fuelfi- m pan ( 10
)
The conservation of energy in the flame leads to the following equation:

™ ent + m fuel hfuel ~ ™pan ^pan ~ Ïflam = (ID
where qji am = net heat transfer to the flame (W)

Finite-difference approximation

The differential equations on the fumes control volume are discretised by a classical finite difference scheme. The second-order accurate central difference formula is used to discretise the first derivatives.

7 HEAT TRANSFER SUB-MODEL

Radiation heat transfer

Governing equations

The ceiling, the floor and the vertical walls of the compartment are assumed to be made in the same material and to have the same behaviour (same thermophysical properties and same surface temperature). Consequently, the inside surface of the enclosure forms an unique opaque surface. In addition to the wall surface, the flame is modelled as an opaque body (a surface enveloping the flame has to be considered). The wall surface and the flame surface are modelled as diffuse-grey bodies. The layer of fumes filling the enclosure is an absorbing-emitting medium. Radiative heat transfer between N surfaces in a participating medium is described by the following equation [START_REF] Siegel | The Engineering Treatment of Gas Radiation in Enclosures[END_REF]:

F 7 ~% •* IT-r * i-lc K-/ /-ft q] = zfo* ~Fk _j Tj_ k )aTf -F^a^aT^) ( 12
)
where The net radiation equation ( 12) applied to the wall surface (S w ) and the flame surface ] ) leads the following equations:

J_ JJ *• & w . " j-, " flam . " "w-w ^w-w rlrad.wall ~ ^w-flam, tfrad.flam, £ w s w J s flam 'V ~ ^w-w^w-w)^ ^-wall ~rw-flam, T flam, -w cr -' flam ~~£ fum°'• * fum F 1 ~£w • " 1 • " ' • * * flamj-w T w-flam] Qrad.wall + Qrad.flam] £ w E flam
• ~Pflami -w T w-flam, a ^wall + a ? flam ~ E fum a ^fum A solvable system of two equations (14 and 15) with two unknowns {q ra d, wa ii an d q'rad flam ) * s obtained. The net radiation equation ( 12) applied to the fuel surface (Sf ue i) and flame surface ( Sjj am ) leads the following equations:

•" aT flam -aT fuel 9rod, flam 2 -~1rad, fuel ~ j j 1

| 1 I e flam £ fuel
The net radiative heat flux leaving flame surface q rad -is given by the following equation:

Qrad,flam = Sflaml Qrad.flaml + ^'flam 2 Qrad,flam2

View factors

The view factor Fy is defined as the fraction of radiant energy leaving surface k and that is incident to surface j. The view factorF^a mi _ w from surfaceSfl ami to wall surface is equal to unity.

From the reciprocity relation, we can obtained F w _j lam :

$ walls Fw-flam, = ^flamj Fflamj-w (18)

Finally, as the wall surface and the flame surface Sjj am form an enclosure, we have the following equation for F w _ w (all the fractions of energy leaving wall surfaces and reaching the surfaces of the enclosure must total to unity):

Gas emissivity

Radiation properties depend on several parameters such as gas properties (temperature, pressure, nature of constituents) wavelength and geometry of the gas volume [START_REF] Siegel | The Engineering Treatment of Gas Radiation in Enclosures[END_REF][START_REF] Tien | Radiation Heat Transfer[END_REF]. For engineering applications, the diffuse-grey assumption is used for flame surface and fumes; hence, mean radiation properties (integrated over the whole spectrum of wavelengths) are used. Moreover, the notion of equivalent mean beam length is used for the flame and the fumes in order to simplify the analysis:

J-6 Vflame J-Wfumes eq, flame <-, eq, fumes <-, Aflame ^ fumes
Emissivities of smoke and flame are calculated by taking into account main contributing emissivities of CO 2 , H 2 O and soot [START_REF] Tien | Flame Radiation[END_REF]): e = e g +e s -s g s s s g = e CO2 + s HlO -As (21)

The emissivities of CO 2 and H 2 O are calculated using Modak correlations [START_REF] Modak | Radiation from Products of Combustion[END_REF]). The emissivity of soot is calculated according to [START_REF] Yuen | A Simple Calculation Scheme for the Luminous-Flame Emissivity[END_REF].

Convection heat transfer

Convection heat transfer is described by the following equation:

4conv,walls ~ ^ ( T walls ~ T gaz )

where, h = Mean convection coefficient to walls (W.m^.K" 1 )

The convective heat transfer coefficient can be obtained from empirical correlations for natural convection based on Nusselt, Grashof and Prandtl numbers [START_REF] Atreya | Convection Heat Transfer[END_REF]:

Nu = F(Gr,Pr) (23)
The convective heat transfer coefficients are approximately the same for vertical and horizontal walls. In order to simplify the procedure, the same weighted average value is taken for the vertical and horizontal walls.

Conduction heat transfer

In radiation and convection heat transfer equations, the wall temperature T waU is still an unknown. This temperature is calculated from the one-dimensional heat conduction model (Ôzisik 1994):

dT(t,x) , Ô 2 T(t,x) pCp ~ET^-^- ( 25
)
where, p, c p and X = respectively the density, the heat capacity and the conductivity of the wall material ; x = space coordinate ; T(t, x) = temperature profile in the wall at time t.

The user may choose the nature of the compartment surfaces from a thermophysical properties database file (containing values of density, heat capacity, conductivity and emissivity of walls).

In order to be coherent with the discretisation of differential equations obtained on the fumes control volume, the Crank-Nicolson finite difference method is applied for discretising the onedimensional heat equation. Hence, equation 25 at any of the interior nodes (i=l ...M-l) is discretised as follows (the walls are subdivided into M equal parts of mesh size):

-rTflY +(2 + 2r)T i n+1 -rTfi 1 = rT/L } +(2-2r)T l n + rTfa (26) .., a At , X
with, r and a = (Ax) 2 pc p where, T" = temperature of node i at the time level n ; a = material thermal diffusivity

Equation 26 provides M-l algebraic equations but contains M+l unknown node temperatures. We need also to consider the two boundary conditions (flux boundary condition for node x=0 and convection boundary condition for node x=L):

x=0, (2 + 2r)T 0 n+1 -2rTf +1 =(2-2r)T 0 n + 2rTj n + 2 r f* (q' o (tj + q o (t n+1 )) (27) x=L, -2rTÛ + Ji +(2 + 2rp L )T^+ 1 = 2rT&_i +(2-2r/3 L )T^ + 4ry L (28)
with,

PL=1 + T^ and y L = ~^T air

A A

where, q" 0 (t n+1 ) = net heat flux (W.m~2) applied to the boundary surface x=0 at time level t n+1 ;

h M = convection heat fransfer coefficient at boundary surface x=L.

The net heat flux received by the inside surface is composed of convection and radiation heat transfers to the inside walls:

"

" "

HO ~ ~~lray,walls ~~ "conv.walls

Equations 26 to 28 lead to a tridiagonal system of M+l algebraic equations with M+l unknown node temperatures T" +I . In order to solve this system by the Thomas' algorithm (which is valuable for linear systems), the source term q 0 in the boundary equation 27 was linearised by the Newton-Raphson iterative method.

As air temperature inside the compartment can be different of ambient air temperature T air , conduction equations must be solved at time level 0 in order to obtain the initial conditions T(x,0). At time level 0, we can point out that we have a convection boundary condition for node x=0.

VALIDATION

Experiments were carried out in a forced ventilated enclosure of 80 m 3 (at INERIS) with pyridine pool fires. The test compartment is approximately 4.9 m wide by 4 m deep by 4 m high (figure 4). The vertical walls, floor and ceiling consist of concrete of 20 mm thickness. The centre of the inlet opening (0.2 m in diameter) is 0.5 m above the floor. The centre of the exit opening (0.5 m x 0.5 m) is 3.25 m above the floor. Pyridine was placed in a steel pan of 0.564 m diameter. This pan was placed on a load cell in the centre of the floor. A total number of five tests with different ventilation conditions were performed (ventilation flow rates decreased successively from 1300 to 280 Nmlh" 1 ).

The ventilation rate of air was measured using a Pitot probe and a thermocouple located inside the inlet duct. Humidity of ambient air was also measured. The ventilation rate of fumes was measured using a bi-directional probe and a thermocouple located in the exhaust duct. Molar fractions of O 2 , CO 2 , CO and NO in the exhaust duct were measured. Molar fraction of THC and generation of HCN in the exhaust duct were measured using respectively a flame ionisation detector and an automatic titrimeter (a heated sampling line was used). Generation of soot was measured using a gravimetric soot measuring device. Molar fraction of O 2 inside the compartment (50 cm from floor -North tree) was also measured.

There were four thermocouple trees to measure gas temperatures inside the compartment. Trees were located in the median plans of the compartment at one meter from the vertical walls. Gas temperatures were measured with Type K thermocouples of 1 mm in diameter. East and west trees contain 2 thermocouples located at 1 m et 2.5 m from floor. South tree contain 15 thermocouples spaced 25 cm apart beginning 25 cm from the floor. North tree contain 8 thermocouples spaced 50 cm apart beginning 25 cm from the floor. hi terms of input data, the model basically requires the knowledge of the mass of fuel, the pool fire (pan) area, the mass flow rate, the temperature and humidity of incoming air and the compartment geometry in relation with a given scenario. In a final step of the initial procedure, the user also chooses the fuel and the nature of compartment surfaces from two database files.

The flame temperature predicted by the model is strongly underestimated. Indeed, the conservation of energy in flame leads to temperature of about 300°C, which is typical of plume temperature. In order to solve this problem, flame temperature must be introduced as an input data. In a forthcoming development, this parameter will be calculated from the conservation of energy on a new control volume, the persistent region of flame. Figure 4 presents the comparisons for fuel burning rate (for the five tests, in steady state conditions) between the experimental results and the model predictions. The smaller is the ventilation rate, the smaller is the fuel burning rate because of vitiation of air entrained into the base of the flame. The model leads to fairly good predictions of fuel burning rate. hi order to analyse more finely the combustion sub-model of the zone model (comparisons of chemical effects between the experiments and the model predictions), the same input parameters were imposed, only that the actual fuel burning rate was also fixed as input data (the yields of combustion products depend on the value of the equivalence ratio which depends on the fuel burning rate). The paper now discusses the results obtained for the experiment carried out with an air ventilation rate of 280 Nm 3 .h"' (the least ventilated test configuration), hi this configuration, we can point out that the equivalence ratio was about 0.55.

For concentrations of O 2 and CO 2 (see fig. 5 and6), there is a good agreement between the predicted and measured values. Prediction of concentrations of O 2 and CO 2 are respectively overestimated and underestimated by 10%. This difference could be the consequence of an overestimation of the measured ventilation rate of the compartment. CO concentrations (ventilation rate =280 Nm 3 .h"')

For concentration of CO (see fig. 7), there is a strong difference between calculated and measured values (zone model version 1). The concentration of CO is underestimated by a factor of about 2. We can point out that molar fraction of oxygen at the flame base has an effect on the CO yield [START_REF] Pitts | The Global Equivalence Ratio Concept and the Formation Mechanisms of Carbon Monoxide in Enclosures Fires[END_REF][START_REF] Mulholland | The Effect of Oxygen Concentration on CO and Smoke Produced by Flames[END_REF]. In our case, molar fraction of oxygen at the flame base is equal to about 13%. Hence, in the sub-model of combustion, we have introduced the correlation of [START_REF] Mulholland | The Effect of Oxygen Concentration on CO and Smoke Produced by Flames[END_REF] in order to take into account the influence of oxygen concentration on the CO yield. As can be seen, there is now a good prediction of CO concentration (zone model version 2).

CONCLUSION

A new zone model for prediction of thermal and chemical effects of pool fire in a forced ventilated enclosure has been developed. The main differences between this model and classical zone models rely on the development of three sub-models that are used for reducing the number of input data needed for a given simulation. The burning rate history of liquid pool fire is calculated from an original sub-model of pyrolysis. A solid flame sub-model has been developed for predicting radiative properties of flame. Yields of chemical species are estimated from a dedicated sub-model of combustion which is based on experimental results obtained at labscale and introduced in a fuel database. This zone model has received some validation for use in forced ventilated enclosures. Promising results have been obtained. For a more general point of view, this is very encouraging as fire safety engineering techniques based on this type of approach are clearly needed for addressing some important safety issues in enclosures such as warehouses or tunnels.

Figure 2 .

 2 Figure 2. Yields of CO 2 and CO (g/g) as a function of equivalence ratio for pyridine

  ç = net radiative heat flux leaving surface j ; Su = the Kronecker symbol ; e , • = emissivity of surface j ; F k j = view factor from surface k to surface j ; T.-_ k = fraction of the radiation from surface j to surface k transmitted through the gas ; ctj_ k = fraction of the radiation from surface j to surface k absorbed by the fumes ; Tj -temperature of surface j ; and Tfy m = temperature of fumes As the flame is opaque, the fuel surface only see the surface of the flame base. The flame surface is divided into two surfaces (of same temperature and same emissivity) Sfl am and

Figure 3 .

 3 Figure 3. Schematic of the INERIS forced ventilated fire compartment

Figure 4 .Figure 6 .

 46 Figure 4 .Comparison of measured and predicted Figure 5 . Comparison of measured and predicted fuel burning rate in steady state conditions O 2 concentrations (ventilation rate =280 Nmlh' 1 )


C Rollin

T Heidi

J Laversanne

Natural attenuation of PAHs under anaerobic conditions

Introduction

Risk assessment for groundwater at contaminated sites uses flow and transport modelling. Taking into account the low mobility of PAH's, simulations are carried out for large time scale and natural attenuation, even if low, cannot be discarded. Sensitivity analysis have shown that the integration of biodégradation in the models has a strong effect on the results and on the conclusions of the risk assessment. At the moment, the position of the French Ministry of Environment is close to the position advocated by the Environmental Protection Agency (EPA) in the USA, based on the principle of precaution : "the phenomena of degradation of chemicals can be considered within the risk assessment only if irrefutable evidences of its existence can be brought ".

Laboratory experiments

Studies have been carried out to check whether natural attenuation of PAH's (mainly sorption and biodégradation) could occur at a disused coke plant site. Two aquifers separated by a clay layer are present at the site location. We performed batch tests in order to measure the partition coefficient of naphthalene on the clay. The kinetics of adsorption are notably fast and the partition coefficient was found to have a value of 7,64 cm 3 /g. We also performed column and microcosm tests in order to check whether biodégradation could occur in the groundwater. A column experiment was performed with soil and groundwater taken from the site. The experimental soil was supplemented with PAH's in the same ratios as found at the source, resulting in a total addition of 200 mg/kg. The soil was then mixed with filter gravel (50 % by volume). The filter gravel increases the hydraulic conductivity and shows nearly no sorption of PAH. The soil was filled into a glass column (10 cm in diameter, 50 cm in length) connected with a reservoir filled with the original groundwater. The water circulated through the soil column from the bottom to the top with the help of a peristaltic pump at a rate of 45 ml/h. The glass reservoirs were wrapped in aluminium foil to avoid photo catalytic degradation. At different time intervals water samples were taken at the sampling port located behind the water reservoir. To ensure anaerobic conditions, the volume of the samples taken was replaced by argon supplied from an argon filled compressible gas bag. After an initial rapid decrease of the concentrations (1 st to 60 th day : phase 1), no degradation occurred from the 60 th to the 100 th day (phase 2). It is assumed that the first period (1 st to 60 th day) is characterised by an adjustment of equilibrium of the PAH concentrations between the soil and the water phase. The decrease of the PAH concentrations initially observed might be due to a slow migration of the pollutants to binding sites on the soil particles (surface and pores). The adaptation of the microflora to the degradation of the PAH under sulfate reducing conditions might have taken about 100 days. Hence, during the 2" d phase (60 th to 100 th day) no decrease of the PAH concentrations occurred. During the last period (100* to 248 th day : phase 3) degradation of all water-bound PAH took place.

In a first, orientating microcosm experiment, the hypothesis that PAH degradation is linked to sulfate reduction was tested. Three microcosms all containing soil were constructed (abiotic control was achieved by poisoning with HgCl 2 ). After an incubation period of 2,5 months, the solid and the liquid phases were separated by centrifugation and the liquid phase only was analysed for its PAH content. The inventory of PAH in the liquid phase was calculated. While the poisoned control exhibited a high residual PAH inventory (1900 u.g), only a small amount of PAH (1,1 u.g) could be determined in the biotic assay. The inhibition of sulfate reduction by addition of molybdic acid resulted in a significantly higher residual PAH inventory (1180 u.g), indicating that the PAH degradation was due to sulfate reduction.

Modelling

The eventual aim of the research is to achieve best knowledge of sorption and biodegration rates in order to fit the constants in the modelling. We have already conducted sensitivity analyses that will be included in the paper.