J Bartel

B Nerlo-Pomorska

K Pomorski

A Dobrowolski

Transport coefficients in the Fourier shape parametrization

Keywords: liquid drop model, nuclear inertia, wall friction PACS: 21.10.Dr, 25.85.-w, 21.60.Ev

An innovative parametrization of nuclear shapes, based on a Fourier expansion of the square distance from the surface of the nucleus to the symmetry axis is introduced. Surface, curvature and Coulomb energy coefficients of a charged liquid drop, that determine the semiclassical nuclear energy, are evaluated within this shape parametrization, together with the wall-friction and the irrotational-flow mass tensors. These transport coefficients are important ingredients of many nuclear models describing nuclear structure and dynamics. A numerical code allowing for the determination of all these quantities for a huge variety of nuclear shapes is made available to the interested user.

PROGRAM SUMMARY

Program Title: inerfric Program Files doi: http://dx.doi.org/10.17632/gtvbc8b7sw.1 Licensing provisions: GPLv3 Programming language: Fortran Nature of problem: The inertia tensor evaluated in the hydrodynamical model using the Werner-Wheeler approximation generalized to nonaxial shapes is evaluated, together with the friction tensor using the wall formula, and the shape functions that define the liquid-drop energy in terms of surface, curvature, Coulomb

Introduction

An accurate description of nuclear shapes which involves as few collective variables as ever possible is a very demanding task, especially in connection with the fission process, which envolves an impressive variety of nuclear shapes. Several powerful shape parametrizations have been developed in the past [START_REF] Hasse | Geometrical Relationships of Macroscopic Nuclear Physics[END_REF]. We will use here a recently developped parametrization which allows to cover a rich variety of shapes with four collective deformation parameters only [START_REF] Schmitt | [END_REF]. This parametrization based on a Fourier expansion of the square distance of any point on the nuclear surface to the symmetry axis (chosen here as the z-axis) is very rapidly converging and easy to handle. In addition, it is able to overcome some of the limitations encountered with previous prescriptions. A detailed analysis of the deformation-energy landscapes, limitted to only 4 defomation parameters is thus presented in this work.

The investigations carried out in Refs. [START_REF] Schmitt | [END_REF]3] demonstrate the high performance of the here proposed approach and enables us to study large amplitude collective phenomena in a rather accurate way, often better as compared to those carried out with more advanced models involving a much larger number of deformation parameters, which, to our understanding, seems to indicate that our approach contains the essential physical ingredients. Our approach can be implemented very easily into any extended dynamical calculations of nuclear potential-energy surfaces as well as into the dynamical calculation related to the nuclear fission or fusion process. In particular, one could have in mind the Langevin approach of dissipative dynamics, coupled with the evaporation of light particles (see e.g. [4,6,7] and references therein).

Section 2 presents a brief summary of the Langevin approach of dissipative dynamics and of all transport functions (frequently called coefficients) entering the Langevin equation. In the following section 3 we give detailed expressions for the liquid-drop shape-dependent coefficients and some global macroscopic quantities characterizing the nuclear mater distribution, like the density multipole moments and the moments of inertia. Section 4 gives the expressions of the nuclear friction tensor and the irrotational-flow mass tensor in our Fourier shape parametrization. At the end we present in section 5 the description of the computer code that allows to evaluate the above transport coefficients and we give a few typical examples with their results. A detailed derivation of the velocity field of an irrotational flow is described in the Appendix.

Langevin dissipative dynamics

Nuclear reactions like the fusion or the fission process are very complex phenomena which are related to a large transfer of nuclear mass as well as a substantial transfer of nuclear energy from the collective to the singleparticle degrees of freedom, which can be considered, in such an approach, as a heat reservoir. The derivative of the collective potential V ( q) with respect to the collective coordinate (here the deformation parameter) q i plays the role of the driving force F i ( q). The collective kinetic energy is a quadratic form built from the inertia tensor M ij and the collective velocities ( qi ) which are the time derivatives of the collective coordinates describing the nuclear shape. The dissipation of collective energy into intrinsic (single-particle) excitation is described by the friction tensor γ ij , and the fluctuations of the collective coordinates which appear in such a statistical approach by the diffusion tensor D ij .

In many applications, it is assumed that the distribution probability w( q, p, t) of finding a nucleus at a given point in the collective phase-space built from the deformation parameters q i and conjugate momenta p i is described by the Fokker-Planck equation [8]:

∂w ∂t = i ∂V ∂q i ∂w ∂p i - i,j M -1 ij p i ∂w ∂q j + i,j ∂ ∂p i k γ ik M -1 kj p j w + i,j D ij ∂ 2 w ∂p i ∂p j . (1)
It has been shown (see e.g. [8]) that such a transport equation, derived on the basis of statistical (stochastic) mechanics, is equivalent to a Langevin equation with a normally distributed random force. One can then determine the transport coefficients that enter the Langevin equation governing the time evolution of the nuclear system

                   dq i dt = j M -1 i,j p j , dp i dt = - dV dq i - 1 2 j,k dM -1 jk dq i p j p k - j,k γ ij M -1 jk p k + F (L)
i ( q, t) .

(

) 2
Here

F (L) i
is the Langevin random force defined as

F (L) i ( q, t) = j g ij ( q ) Γ j (t) , (3)
where

k g ik g kj = D ij (4)
and Γ j (t) are Gaussian-distributed random numbers with vanishing mean value < Γ j (t) > and delta-correlated variance

< Γ i (t) Γ j (t ) >= 2 δ ij δ(t -t ) . (5)
In the following section we are going to collect the expressions that allow to evaluate, in a macroscopic approximation, all the transport functions entering the Fokker-Plank (1) or Langevin (2) equations.

The transport function in the Fourier shape parametrization

The deformation dependent transport functions, like the nuclear collective potential-energy surface, collective inertia, friction and diffusion tensors, which enter the Langevin or Focker-Planck equations should be evaluated using a given shape parametrization. In our present study we use the Fourier parametrization of deformed shapes written in cylindrical coordinates (ρ, ϕ, z), a parametrization that is not only rapidly converging as was shown in Ref. [START_REF] Schmitt | [END_REF], but it is also simple, analytical and not closed, which means that its convergence can be tested (contrary e.g. to the famous Funny-Hills shape parametrization [9]). To our opinion this is one of the best parametrizations which can be used both in nuclear-structure and nuclear-dynamics calculations.

The surface of an axially symmetric nucleus, or more precisely the square distance of a surface point to the symmetry z-axis, can be expanded in a Fourier series in the following way [START_REF] Schmitt | [END_REF]:

ρ2 s (u) R 2 0 = ∞ n=1 a 2n cos (2n -1)π 2 u + a 2n+1 sin 2nπ 2 u ( 6
)
where R 0 is the radius of a spherical nucleus having the same volume as the deformed one and u = (z -z sh )/z 0 , where 2z 0 is the length of the nucleus and z sh is a parameter that guarantees that the nuclear centre of mass is located at the origin of the coordinate system. More general nonaxial shapes can be obtained by assuming that any crosssection of the nucleus perpendicular to the z-axis has the form of an ellipsoid as shown in Fig. 1:

ρ 2 s (z, ϕ) = ρ2 s (z) f (ϕ) = ρ2 s [u(z)] 1 -η 2 1 + η 2 + 2η cos(2ϕ) . (7)
Here the non-axiality parameter η is a function of the coordinate z and defined by the ellipse half-axis a(z) and b(z):

η(z) = b(z) -a(z) b(z) + a(z) where a(z) • b(z) = ρ 2 s (z) . (8)
For the sake of simplicity, we will assume in the following that η is a constant. The nucleus being considered as an essentially incompressible fluid, the volume-conservation condition requires that

z 0 +z sh -z 0 +z sh dz 2π 0 dϕ ρ(z) 0 ρ dρ = 4π 3 R 3 0 , (9)
what leads to the relation:

∞ n=1 (-1) n-1 a 2n 2n -1 = πR 0 3z 0 . (10)
The above equation allows to evaluate the length of the deformed nucleus.

The condition that the centre of mass should be located at the origin of the coordinate system leads to the following expression for z sh :

z sh = 3 z 2 0 2πR 0 ∞ n=1 (-1) n a 2n+1 n . (11)
Using Eq. (7) and the above relations one can evaluate all transport functions entering the Langevin or Fokker-Plank equations. One usually assumes that the Einstein relation

D ij = γ ij T , (12)
between the diffusion and the friction tensor holds, where T is the nuclear temperature. At larger temperatures the friction tensor is well approximated by the wall-formula [10]

γ ij = ρ 0 2 v zmax z min dz 2π 0 dϕ ∂ρ 2 s ∂q i ∂ρ 2 s ∂q j 4ρ 2 s + 1 ρ 2 s ∂ρ 2 s ∂ϕ 2 + ∂ρ 2 s ∂z 2 (13)
frequently used in dissipative-dynamics calculations related to the fission of hot nuclei. The irrotational flow inertia-tensor in the Werner-Wheeler approximation [11] can be written as

M ij = ρ 0 V A ρ i A ρ j + A ϕ i A ϕ j ρ 2 s (z, ϕ) + A z i A z j d 3 r , (14)
where the A ν j are the expansion coefficients of the velocity field in ν direction (ν = z, ρ, ϕ) which are derived in the Appendix.

The potential energy of a deformed (deformation parameters q ), hot (temperature T ) and rotating (angular momentum L) nucleus is defined by the difference of the Helmholtz free energies F between deformed and spherical nucleus:

V ( q , T, L) = F ( q, T, L) -F ( q = 0, T, L) (15)
where the free energy is given by

F ( q, T, L) = E( q, L) -a( q )T (16)
with a( q ) the deformation dependent level-density parameter (see e.g. [13]) .

In the approximation where the nucleus can be considered as a charged liquid drop, the nuclear energy is given as the sum of volume, surface, curvature, Coulomb, and rotational energies:

E( q, L) = E vol + E surf ( q ) + E cur ( q ) + E Coul ( q ) + E rot ( q , L) , (17)
The first term, the volume term, in eq. ( 17) is deformation independent (the nucleus is considered to be incompressible), while the variation with deformation of all the other enery contributions can be described by the following dimensionless functions:

B surf ( q ) = E surf ( q )/E surf (sph) , B cur ( q ) = E cur ( q )/E cur (sph) , B Coul ( q ) = E Coul ( q )/E Coul (sph) . ( 18
)
Taking derivatives of the profile function (7), one is able to determine the shape functions corresponding to the surface energy

B surf = 1 4πR 2 0 z 0 -z 0 dz 2π 0 dϕ ρ 2 s (z, ϕ) + ∂ρ s ∂ϕ 2 + ρ s ∂ρ s ∂z 2 , (19)
to the curvature energy

B cur = 1 8πR 0 zmax z min dz 2π 0 dϕ ρ 2 s + ρ 2 s ∂ρ s ∂z 2 + ∂ρ s ∂ϕ 2 -1 • ρ 2 s ∂ρ s ∂z 2 -ρ s ∂ρ s ∂z 2 ∂ 2 ρ s ∂ϕ 2 + ρ 2 s -ρ s ∂ 2 ρ s ∂ϕ 2 + 2 ∂ρ s ∂ϕ 2 +2 ∂ρ s ∂z ∂ρ s ∂ϕ ρ s ∂ 2 ρ s ∂z∂ϕ -ρ 3 s ∂ 2 ρ s ∂z 2 -ρ s ∂ 2 ρ s ∂z 2 ∂ρ s ∂ϕ 2 , (20)
and to the Coulomb energy

B Coul = 5 64 π 2 R 5 0 zmax z min dz 2π 0 dϕ zmax z min dz 2π 0 dϕ • ρ 2 s + ρ 2 s -2 ρ s ρ s cos(ϕ -ϕ ) + (z -z ) 2 -1/2 • ρ 2 s -ρ s ρ s cos(ϕ-ϕ ) -ρ s sin(ϕ -ϕ ) ∂ρ s ∂ϕ -ρ s (z -z ) ∂ρ s ∂z • ρ 2 s -ρ s ρ s cos(ϕ -ϕ ) + ρ s sin(ϕ -ϕ ) ∂ρ s ∂ϕ + ρ s (z -z ) ∂ρ s ∂z , (21)
where ρ s is to be understood as a function of z and ϕ .

One finally needs to take the rotational degree of freedom into account. In our semiclassical approach, the rotational energy is given by

E rot ( q ) = L 2 2J ( q ) , ( 22
)
where L is the angular momentum and J the moment of inertia with respect to the rotation axis. The relative change of the rotational energy with deformation is then simply

B rot = E rot ( q ) E rot (sph) = J (sph) J ( q ) , (23)
where J (sph) is the moment of inertia of the spherical nucleus. For a rotation around the x-axis, the parameter B rot is given, in the rigid-body approximation, by the following expression:

B rot = 32 πR 5 0 15   zmax z min dz 2π 0 dϕ [ 2 z 2 ρ2 s + ρ4 s sin 2 (ϕ) ]   -1 . ( 24
)
All above functions entering the Langevin equation ( 2) and the macroscopic part of the collective potential (15) can be evaluated using the numerical program described in the next section which extensively uses the Fourier parametrization (7) of nuclear shapes.

Description of the numerical program

The computer code inerfric.f attached to the present paper is written according to the FORTRAN 90 and FORTRAN 95 standards. The module "constants " contains all physical and mathematical parameters used in the computer code. The parameters 'ng1', 'ng2', 'ng3' denote the order of the Gauss-Legendre integration quadratures (with ng1 = ng3 to avoid divergencies in the calculation of the Coulomb energy), while 'no' gives the number of mesh points used in the rectangular (Romberg) integrations. All these parameters can, of course, be freely adjusted by the user. The values proposed in the code have, however, been tested and turned out to be sufficient to obtain the numerical results with the chosen accuracy (we have e.g. required the LDM energies to be determined with an accuracy better than 0.1 MeV). If desired by the user, he can easily chose to improve that accuracy. The dimensionality 'ndim' of the {q i } deformation space is originally set to 'ndim = 6', but can be extended by the user up to 'ndim = 9'. Single-precision words are used in the program, but it is strongly recommended to compile it with the autodouble option in order to obtain the necessary precision of the transport coefficients needed in nuclear-physics calculations. Within the linux system that would be:

gfortran -fdefault-real-8 -o inerfric inerfric.f The command ./inerfric < inerfric.inp > inerfric.out& would then execute the program taking the input data from the inerfric.inp file and storing the results into the inerfric.out file.

The code allows to evaluate the deformations dependent liquid-drop surface (BS), (first order) curvature (BK), Coulomb (BC) and rotational energy coefficients (Bx,Bz) for rotation around the x-and z-axes respectively. In addition it gives the deformation dependence of the Wigner (congruence) energy (BW) according to Ref. [14], fission fragment mass ratio (BF) as well as the relative distance (R12) between the fragments, the length of nucleus (c) and the neck radius (Rn) in units of the spherical radius (R 0 ) and the quadrupole moment of the nucleus in R 2 0 units times the total mass or total charge of the nucleus. All above deformation dependent functions are evaluated in the subroutine "fcs ".

The inertia tensor in the irrotational flow approximation and the wall friction tensor are evaluated as well in the subroutine "tcoef ". These tensors T ij being symmetric, only the upper part (j ≥ i) needs to be evaluated in the four dimensional deformation-parameter space composed of non-axiality (q 1 = η), elongation (q 2 ), left-right asymmetry (q 3 ) and neck (q 4 ) degrees of freedom. One defines new collective coordinates

q 2 = a (0) 2 /a 2 -a 2 /a (0) 2 ,
q 3 = a 3 , q 4 = a 4 + ((q 2 /9) 2 + (a (0) 4 ) 2 ) , q 5 = a 5 -(q 2 -2)a 3 /10 , q 6 = a 6 -((q 2 /100) 2 + (a (0) 6 ) 2 ) , q 8 = a 8 + ((q 2 /300) 2 + (a

(0) 8 ) 2 ) , (25)
where the a (0) 2n = (-1) n-1 32 π 3 (2n -1) 3 and a

(0) 2n-1 = 0 (26)
are the Fourier components of the spherical shape.

The new parameters q i are chosen in such a way that the situations where all q i vanish corresponds to the spherical shape of the nucleus. Oblate shapes are obtained for q 2 < 0 and prolate ones for q 2 > 0. The trajectory S(q 2 ) for which q 1 = q 3 = q 4 = • • • = 0 corresponds, to a good approximation, to the liquid-drop path to fission as demonstated in Fig. 2. The deformation parameter q 5 is chosen in such a way that for a given elongation (triggered by q 2 ) the stiffness of the liquid-drop potential-energy surface with respect the left-right asymmetry parameter q 3 is as small as possible. The transforma-tion between the q i deformation parameters and the Fourier decomposition coefficients a i is made in the subroutine "qtoa ".

Figure 2: (Colour online) LD potential energy surface in the (q 2 , q 4 ) (top left), (q 2 , q 3 ) (top right), (q 2 , q 6 ) (bottom left) and (q 2 , q 5 ) (bottom right) deformation planes of a nucleus with fissility x=0.8. In each (q i , q j ) plot, the collective coordinates q k with k = i, j are set to zero.

To obtain the friction and the inertia tensor in the space of the q i parameters, we first calculate these quantities in the Fourier-parameter space to obtain the γ a i ,a j and M a i ,a j tensors and then transform these to the {q i } space according to

M q k q l = i,j M a i a j ∂a i ∂q k ∂a j ∂q l , γ q k q l = i,j γ a i a j ∂a i ∂q k ∂a j ∂q l , (27)
a calculation that is carried out at the end of the subroutine "tcoef ".

The components of the inertia and friction tensors obtained in this way are shown Fig. 3 along the LD path to fission as function of the relative distance R 12 between the fragments. One notices that for large values of R 12 the inertia approaches (as it should) to the reduced mass µ. We also show on the same figure the approximation for the inertia made in Ref. [15]. ) tensors as functions of the relative distance between the mass centers of the nascent fragments. The approximation for the inertia presented in Ref. [15] is shown as the dotted line.

Summary

The main purpose of the present work, and of the computer code that comes along with it, is the calculation of the transport functions necessary to perform studies of the nuclear fission process within the dissipative dynamics model. In addition the program evaluates the deformation dependent liquid-drop coefficients which are important ingredients of all macroscopicmicroscopic models. The rapidly converging Fourier type parametrization of nuclear shapes used in the program offers a very effective four dimensional (4D) description of the potential-energy surfaces of fissioning nuclei and their fission dynamics.

ρ = ρ B(z, s, ˙ s ) = ρ ρ s (z, ϕ) i B i (z, s ) ṡi . (A.7)
A natural generalisation of this approximation to the case of non-axial shapes, as given e.g. by Eq. ( 7), is the following

v x = (α + β) x = (α + β) ρ cos ϕ (A.8) v y = (α -β) y = (α -β) ρ sin ϕ
where the parameters α and β will generally depend on z and the collective coordinates s i . In cylindrical coordinates (x = ρ cos ϕ and y = ρ sin ϕ) the x and y components of the velocity field are given by v x = ẋ = ρ cos ϕ -ρ φ sin ϕ (A.9) v y = ẏ = ρ sin ϕ + ρ φ cos ϕ .

Using the above relations (A.8) and (A.9) we get the following expressions:

ρ = ρ (α + β cos 2ϕ) (A.10) φ = -β sin 2ϕ .
The velocity field of an ideal fluid should be both source-less and irrotational, i.e.

∇ • v = 0 and ∇ × v = 0 . (A.11)
The parameters α and β can now be evaluated from the divergence-free character of the velocity field, Eqs. (A.11), and the condition that the velocities of the wall and of the liquid should be equal at the surface of the body defined by ρ s (z, ϕ). For the velocity components on the x (ϕ = 0) and on the y axis (ϕ = π/2) one obtains more specifically

v x (ρ s , ϕ = 0, z) = ρs (z, ϕ = 0) = ρ s (z, ϕ = 0) [α(z) + β(z)] (A.12) v y (ρ s , ϕ = π 2 , z) = ρs (z, ϕ = π 2 ) = ρ s (z, ϕ = π 2 ) [α(z) -β(z)]
where, in the last step, Eq. (A.10) has been used and which directly leads to the following identities for α(z) and β(z)

α(z) = 1 2 ρs (z, ϕ = 0) ρ s (z, ϕ = 0) + ρs (z, ϕ = π 2 ) ρ s (z, ϕ = π 2 ) (A.13) β(z) = 1 2 ρs (z, ϕ = 0) ρ s (z, ϕ = 0) - ρs (z, ϕ = π 2 ) ρ s (z, ϕ = π 2 )
.

Let us mention here that, in the general case, the velocity components v x and v y along the x and y axis as given by Eq. (A.12) are, of course, not the only non-vanishing velocity components, but that, except at the zcoordinates where ∂ρ 2 s (z, ϕ)/∂z = 0, the z component v z of the velocity field is also present. But that is of no concern to us, since the coefficients α(z) and β(z) do only enter the definition of the x and y components of the velocity field through Eq. (A.8) and are unambiguously defined through Eq. (A.12).

Using the product ansatz of ρ 2 s (z, ϕ) as given in Eq. ( 7) one has

ρs (z, ϕ) ρ s (z, ϕ) = 1 ρ2 s (z) f (ϕ) d dt ρ2 s (z) f (ϕ) = 1 2 ˙ρ 2 s (z) ρ2 s (z) + ḟ (ϕ) f (ϕ) (A.14)
and with the particular definition for f (ϕ) of Eq. ( 7) one obtains ḟ (ϕ)

f (ϕ) = - 2 η 1 -η 2 2η + (1 + η 2 ) cos(2ϕ) 1 + η 2 + 2η cos(2ϕ) (A.15)
which then, together with (A.13), yields

α(z) = 1 2 ˙ρ 2 s (z) ρ2 s (z) = ˙ρ s (z) ρs (z) (A.16) and β(z) = η η 2 -1 . (A.17)
an expression that is only related to the non-axiality parameter η and its variation with time, but independent of z, something that we shall assume to be generally true in what follows, while the expression for α(z) will be determined below through Eq. (A.11).

Inserting α and β as found above into Eq. (A.10) one obtains the following expression for the velocities in a plane of constant z

ρ = ρ α(z) - η cos(2ϕ) (1 -η 2 ) and φ = η sin(2ϕ) (1 -η 2 ) . (A.18)
From the condition div v = 0, Eq. (A.11), one obtains together with (A.8)

∂v z ∂z = - ∂v x ∂x - ∂v y ∂y = -2α , (A.19)
which corresponds, in the case of axial symmetry, to Eq. (A7) of [11]. Proceeding as in Ref. [11], we suppose that, for an incompressible fluid, the convective time derivative of any fluid element of volume

V * (z; s) = 1 2 z z min 2π 0 ρ 2 s (z , ϕ) dz = π z z min ρ2 s (z ) dz must vanish: d dt V * (z; s) = ż ∂ ∂z V * (z; s) + i ṡi ∂ ∂s i V * (z; s) = 0 ,
which can be used to evaluate the z-component of the fluid velocity The kinetic energy of the liquid can then be written in the form:

v z (z) = - 1 ρ2 s (z) i ṡi z z min ∂ ρ2 s (z ) ∂s i dz = - 1 ρ2 s (z)
T = 1 2 ρ 0 V v 2 d 3 r = 1 2 i,j
B ij ṡi ṡj .

Figure 1 :

 1 Figure 1: (Colour online) Schematic visualization, in cylindrical coordinates, of the parameters entering the definition of the profile function defined by Eq. (7).

Figure 3 :

 3 Figure3: Relative distance component of the irrotational flow inertia (l.h.s) and the wallfriction (r.h.s.) tensors as functions of the relative distance between the mass centers of the nascent fragments. The approximation for the inertia presented in Ref.[15] is shown as the dotted line.

  v z is linear in the collective velocities ṡiv z = ż = i A zi (z; s) ṡi , Eq. (A.19) we can evaluate the function α(z) velocity components that appear in Eq. (A.5) are also linear in the collective velocitiesv ρ = ρ = ρ i A ρ i ( r; s) ṡi (A.24) v ϕ = ρ φ = ρ i A ϕ i ( r; s) ṡiwhere our collective coordinates s i are the Fourier deformation parameters a ν , and η. The expansion coefficients of the velocity field are then given by

ρ 2 2 s 4 s

 224 (A.26) Inserting Eqs. (A.21) and (A.24) for the velocity components one obtains the following expression for the inertia tensor B:B ij = ρ 0 A ρ i (z)A ρ j (z) + A ϕ i (z)A ϕ j (z) + A z i (z)A z j (z) ρ dρ . (z, ϕ) cos(2ϕ)dϕ = -2πη , (z, ϕ) cos(2ϕ)dϕ = -2π ρ4 s (z) 2η 1 -η 2 .(A.28) one obtains the final expressions for the components of the inertia tensor:

Acknowledgements

This work has been partly supported by the Polish-French COPIN-IN2P3 collaboration agreement under project number 08-131 and by the Polish National Science Centre, grant No. 2016/21/B/ST2/01227.

Appendix A. Expansion coefficients of the velocity field

The kinetic energy of an ideal incompressible fluid of constant matter density m = 3M/(4πR 3 0 ) and velocity field v( r ) can be written as

with the integration extending over the volume V of the matter distribution.

We will show that this kinetic energy can be expressed as a quadratic form in the collective velocities ṡi associated with the collective coordinates s i which could be identified in our case as the Fourier shape parameters a i entering Eq. ( 6), the non-axiality parameter η, Eq. ( 8), or the new collective coordinates defined by (25)

where M ij are the components of the collective mass tensor.

In order for Eq. (A.1) to imply Eq. (A.2) it is necessary that the position vector r of a fluid element will not depend on the time explicitly, but only on the shape of the system, thus implying that its velocity can be written as

which directly leads to Eq. (A.2) with

In the general case where no axial symmetry is imposed, the velocity field in cylindrical coordinates {ρ, ϕ, z} is written as

where e ρ , e ϕ and e z denote the unit vectors in ρ, ϕ and z directions respectively.

In the axially symmetric case, the Werner-Wheeler approximation consists in supposing that ż is independent of ρ and that ρ depends linearly on ρ, i.e. ż = A(z, s, ˙ s ) = i A i (z, s ) ṡi (A.6) and


Lan Fei

Jean-Yves Yves Dantan

Cyrille Baudouin

Shichang Du

Calibration and uncertainty estimation of non-contact coordinate measurement systems based on Kriging models

Keywords: coordinate measurement system, non-contact measurement, measurement uncertainty, calibration, Kriging model

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

Introduction

Since three dimensional (3D) coordinate measurement systems (CMSs) with optical imaging devices offer the option to output high density data for parts at high speed, they have been widely used in lab and industry [START_REF] Kolen | Self-calibration/compensation technique for microcontroller-based sensor arrays[END_REF][START_REF] Bacakoglu | A three-step camera calibration method[END_REF][START_REF] Sansoni | Calibration and performance evaluation of a 3-D imaging sensor based on the projection of structured light[END_REF][START_REF] Shen | Automatic camera calibration for a multiple-sensor integrated coordinate measurement system[END_REF]. In recent years, non-contact optical measurement systems have not only been seen with significant advances in technical developments, but also undergone a boost in applications, such as surface quality improvement [START_REF] Du | A selective multiclass support vector machine ensemble classifier for engineering surface classification using high definition metrology[END_REF][START_REF] Du | An adaptive support vector machine-based workpiece surface classification system using high definition metrology[END_REF][START_REF] Nguyen | Characterization of Cutting Force Induced Surface Shape Variation in Face Milling Using High-Definition Metrology[END_REF][START_REF] Nguyen | High-Definition Metrology Enabled Surface Variation Control by Cutting Load Balancing[END_REF], surface error evaluation [START_REF] Wang | 3D Surface form Error Evaluation using High Definition Metrology[END_REF], surface separation [START_REF] Du | A Shearlet-Based Separation Method of 3D Engineering Surface Using High Definition Metrology[END_REF][START_REF] Wang | A diffusion filter for discontinuous surface measured by high definition metrology[END_REF], and tool wear monitoring [START_REF] Wang | Tool Wear Monitoring of Wiper Inserts in Multi-insert Face Milling Using 3D Surface Form Indicators[END_REF]. However, for actual camera lenses, such as a fixed length lens, a zoom lens or even an expensive high-quality telecentric lens, image distortions unavoidably exist due to lens aberrations, misalignment of optical elements and non-parallelism between image plane and sensor plane [START_REF] Pan | Systematic errors in two-dimensional digital image correlation due to lens distortion[END_REF][START_REF] Cui | A generalized reference-plane-based calibration method in optical triangular profilometry[END_REF][START_REF] Pan | Systematic errors in two-dimensional digital image correlation due to lens distortion[END_REF]. Due to the non-uniform characteristics of camera lens distortion, actually the coordinates on object surface may cause different values in the sensor plane. Therefore, measurement results are usually affected by measurement uncertainty, which could lead to technical and economic risks in industrial applications.

Over the last 30 years, two different methods have been adopted to achieve high accuracy, namely, uncertainty avoidance and calibration. Uncertainty avoidance is the traditional method focusing on the reduction of possible uncertainty sources in the design and manufacturing processes [START_REF] Huang | On-line error compensation of coordinate measuring machines[END_REF]. Careful design and precise construction can reduce the uncertainty, but every subsequent micrometer/nanometer of uncertainty reduction will cause an exponentially increasing cost. So it is practical to calibrate CMSs between the measurement coordinate and the actual object.

The calibration and uncertainty estimation method has been explored as a cost-effective method of improving the accuracy of CMSs [START_REF] Satoa | Calibration and uncertainty evaluation of single pitch deviation by multiple-measurement technique[END_REF][START_REF] Barini | Uncertainty analysis of point-by-point sampling complex surfaces using touch probe CMMs_ DOE for complex surfaces verification with CMM[END_REF][START_REF] Hong | Non-contact R-test with laser displacement sensors for error calibration of five-axis machine tools[END_REF][START_REF] Kajima | Calibration of linear encoders with sub-nanometer uncertainty using an optical-zooming laser interferometer[END_REF]. Early developments in calibration are well described by Evans [START_REF] Evans | Precision Engineering: An Evolutionary View[END_REF]. Some effective methods are reported in the literature to model and calibrate CMSs. These methods include polynomial models and autoregressive moving average models [START_REF] Chen | Dynamic calibration and compensation of a 3D laser radar scanning system[END_REF], homogeneous coordinate transformation (HTM) methods [START_REF] Caja | Optical measuring equipments. Part I_ Calibration model and uncertainty estimation[END_REF], neural network based methods [START_REF] Shen | Error compensation of coordinate measurements in computer-Integrated manufacturing using neural networks[END_REF][START_REF] Yang | Error compensation of touch trigger probes[END_REF][START_REF] Dong | Prediction and compensation of dynamic errors for coordinate measuring machines[END_REF][START_REF] Tan | Geometrical error compensation of precision motion systems using radial basis function[END_REF][START_REF] Tan | Dynamic S-function for geometrical error compensation based on neural network approximations[END_REF], D-H models [START_REF] Barakat | Kinematic and geometric error compensation of a coordinate measuring machine[END_REF][START_REF] Li | The calibration and error compensation techniques for an articulated arm CMM with two parallel rotational axes[END_REF], fuzzy error interpolation techniques [START_REF] Bai | On the comparison of bilinear, cubic spline, and fuzzy interpolation techniques for robotic position measurements[END_REF], iterative learning and decoupling methods [START_REF] Teo | Dynamic geometric compensation for gantry stage using iterative learning control[END_REF][START_REF] Yang | Numerical Derivation-based serial iterative dynamic decoupling-compensation method for multiaxis force sensors[END_REF], recursive least squares identification techniques [START_REF] Wei | The geometric dynamic errors of CMMs in fast scanning-probing[END_REF], Monte Carlo method [START_REF] Wen | Adaptive Monte Carlo and GUM methods for the evaluation of measurement uncertainty of cylindricity error[END_REF] and other analytical methods [START_REF] Feng | Analysis of digitizing errors of a laser scanning system[END_REF][START_REF] Isheil | Systematic error correction of a 3D laser scanning measurement device[END_REF][START_REF] Johnson | Dynamic error characteristics of touch trigger probes fitted to coordinate measuring machines[END_REF][START_REF] Pereira | Characterization and compensation of dynamic errors of a scanning coordinate measuring machine[END_REF][START_REF] Schwenke | Geometric error measurement and compensation of machines-an update[END_REF][START_REF] Zhang | Error compensation of cylindrical coordinate measuring machines[END_REF].

So far, limited research on the spatial statistics-based methods to model and calibrate CMSs considering the spatial correlation of the measurement data has been conducted. The spatial correlation arises frequently in data measured within certain intervals of space, and the data indeed exhibits a significant amount of positive autocorrelation [START_REF] Suriano | Progressive measurement and monitoring for multi-resolution data in surface manufacturing considering spatial and cross correlations[END_REF]. Measurements are often spatially correlated because they are obtained in similar manufacturing conditions and related to similar properties of the machined material [START_REF] Poniatowska | Research on spatial interrelations of geometric deviations determined coordinate measurements of free-form surfaces[END_REF]. Spatial correlation is different from temporal correlation, which is usually represented via time series models. In fact, spatial correlation models allow one to represent contiguity in space rather than in time. Spatial statistics, e.g. Kriging method, is one of most important meta-models for describe spatial correlations in random field [START_REF] Chen | Enhancing stochastic Kriging metamodels with gradient estimators[END_REF]. Detailed descriptions of existing research on Kriging meta-models are provided in a review [START_REF] Kleijnen | Kriging metamodeling in simulation: A review[END_REF].

The literature on spatial statistics-based calibration methods is sparse. A method based on Co-Kriging models to estimate the errors of the surface form not only concerning spatial correlation but also concerning the influence of machining conditions is developed in [START_REF] Du | Co-Kriging method for form error estimation incorporating condition variable measurements[END_REF]. A spatial statistical method to design adaptive inspection plans for the geometric control of mechanical parts with CMSs is presented in [START_REF] Ascione | Adaptive inspection in coordinate metrology based on kriging models[END_REF]. The prediction uncertainty of geometric deviations is provided using Kriging models. A Kriging-based procedure to identify the minimum of measured points to check the conformity with a given confidence level in the inspection of large surfaces is presented in [START_REF] Dumas | AK-ILS: an active learning method based on Kriging for the inspection of large surfaces[END_REF]. An on-line inspection system with error calibration to obtain high machining accuracy for free-form surface components is proposed in [START_REF] Chen | Spatial statistical analysis and compensation of machining errors for complex surfaces[END_REF]. Through the spatial statistical analysis of the residual errors of a regression model, the errors are decomposed into systematic errors and random errors. The spatial error compensation methods for computer numerical control (CNC) machining center is proposed in [START_REF] Fan | Unified error model based spatial error compensation for four types of CNC machining center: Part II-unified model based spatial error compensation[END_REF][START_REF] Fan | Unified error model based spatial error compensation for[END_REF].

In these last applications, Kriging-based procedure to improve inspection has proved to be effective. This paper develops a Kriging-based method to achieve the calibration and uncertainty estimation of non-contact 3D CMSs. The advantages of this proposed solution are many in comparison to the classical approaches:

-It is not necessary to develop a complex phenomenological model of the sensor.

The systematic error of a simplified model of the sensor was corrected by the Kriging model.

-It is possible to adopt an adaptive approach for the calibration. Based on the uncertainty estimation of the kriging model, it is possible to identify the spatial area that should be recalibrated.

The proposed Kriging-based method for the calibration of non-contact 3D CMSs is presented in Section 2. An application of the approach is illustrated through a case study in Section 3.

The proposed calibration method

The aim of this proposed calibration method is to correct the systematic error of the simplified model of the sensor by a kriging model. To do so, a gauge is measured. It is important to know the geometrical characteristics of this gauge. In this application, the GOM gauges are used.

The calibration procedure

Fig. 1 shows the calibration procedure of a measurement system, and the main steps are described as follows.

[Insert Fig. 1 here

i i i i s x y z - - - = (1)
where measured measured measured , , ,  )

i i i i s x y z - - - = (2)
where theoretic theoretic theoretic , ,

i i i x y z - - -
are three dimensional coordinates of the theoretical reference points respectively.

Illustration: Case of the sphere (Figure 1) -To define the theoretical reference points, the gauge sphere is fitted to the measured points and each measured point is projected on the fitted gauge sphere.

Step 2. Calculating the measurement errors. A measurement error is the distance between the measured point and theoretical reference point at each direction of the axis X, Y, and Z,

measured measured measured theoretic theoretic theoretic ( , , ) ( , , ) ( , , )
i i i i i i i i i i s x y z x y z x y z - - - - - - ∆ = ∆ ∆ ∆ = - (3)
A measurement error is discrete in each measured point on the surface, which needs to be calibrated. It is necessary to use spatial interpolation method to estimate the whole surface of the new measured part surface to realize the calibration of 3D measurement systems. The measurement error estimation model is established and applied at any point of the whole standardized part surface.

Step 3. A Kriging model is used to estimate the calibration 0 0 0 0 ( , , ) C

x y z ∆ = ∆ ∆ ∆ on the specified untried site 0 s by the response. The measurement errors ( , , )

i i i i C x y z ∆ = ∆ ∆ ∆ at a set of measured points -measured measured measured measured ( , , )
i i i i s x y z - - - = (see Fig. 2) are represented as, 0 measured measured measured measured 0 measured measured measured measured 0 measured measured measured measured ( , , ) ( ) ( , , ) ( ) ( , , ) ( )
i i i i i i i i i i i i x u x y z u s y v x y z v s z w x y z w s - - - - - - - - - - - - ∆ = = ∆ = = ∆ = = (4)
where , , u v w are the Kriging models with different regression models and correlation models respectively.

[Insert Fig. 2 here]

Step 4. The calibration of any areas can be estimated by the Kriging spatial interpolation model. Therefore, this estimation could be used to calibrate the measurement error for a new measurement, which realizes the goal of calibration of 3D CMSs.

Kriging-based calibration model

( ) ( ), C D µ δ ∆ = + ∈ s s s s ( )
) 5
where s is the variable [Insert Fig. 3

here]

In order to develop the calibration model, the form of the deterministic structure ( ) µ ⋅ can be specified by three types of Kriging models:

-The simple Kriging: ( ) constant µ = s . The constant is known.

-The ordinary Kriging: 0 ( )

µ β = s
. The constant is unknown and estimated from the data.

-The universal Kriging: The measurement error

( ) C ∆ s for D ∈ s
is considered as a realization of a Gaussian random process:

( ) ( ) ( ) C f β δ ′ ∆ = + s s s (6)
where 1 2

( ) ( ( ), ( ), , ( ))

p ′ = L f s f s f s f s is a set of specified trend functions and 1 2 ( , , , ) p ′ = L β β β β is a set of coefficients.
The measurement error ( ) C ∆ s satisfies a Gaussian random process with zero mean and stationary covariance.

E( ( )) ( ) Cov( ( ), ( )) ( , ) C f C C d R d δ β σ θ ′ ∆ = ∆ ∆ + = s s s s (7)
where 2 δ σ is the process variance, and R is the stationary correlation function (SCF).

The SCF depends only on the displacement vector d between any pair of points in D and on a set of hyper-parameters θ .

, 1 ( , , ) ( , )

n i j j i j i j R R θ θ = = - ∏ s s s s (8)
Base on equation ( 2) and (3), the joint random variable, measurement error,

0 1 2 ( C( ), C( ), C( ), , C( )) n ∆ ∆ ∆ ∆ L s s s s satisfies a multivariate normal function, 2 0 (( , ) , ) N f F δ β σ ′ ′ Σ with 0 0 1 r r R ′   Σ =     (9)
where 0 r is the correlation vector

0 1 0 0 ( ) [ ( , , ) ( , , )] T m r R R θ θ = L s s s s s
, and R is the n n × correlation matrix whose ( , ) i j element is (

) ij i j R = - d s s .
According to the best unbiased estimate condition of Kriging method, the mean of estimation of the calibration value ( ) C ∆ s at the point 0

s is: ( ) ( ) ( ) ( ) ( ) n C C f r f r R C F µ β γ β β - ∆ = ∆ = ⋅ + ⋅ = + ∆ - s s s s (10)
where 0 f is the

1 p × vector of the trend functions in 0 s , F is the n p × matrix of 1, , , 1, , { ( )} j i i n j p f = = L L s of the trend functions calculated in 1 2 ( , , , ) n L s s s , β is unknown coefficient, and it is generalized by LS estimator * 1 1 1 ˆ( ) F R F F R C β - - - ′ ′ = ∆ .
Moreover, a confidence index about the prediction, called the Kriging variance, is the mean squared error between the prediction

ˆ( ) C

∆ s and the actual value ( ) C ∆ s . It is given by:

2 2 2 ' 1 ' ' 1 1 ˆ0 0 0 0 0 0 0 ( ) (( ( )) ) (1 ( ) ) C E C C r R r g F R F g δ σ σ - - - ∆ = ∆ - = - + s s (11)
where

' 1 0 0 0 g f F R r - = - , 0 0 0 0 1 1 : ( ), : ( ); :[ ( ) ( )], :[ ( , ) ( , )] n n f f r r F f f R r r L L s s s s s s s s and 2 δ σ is the variance of the process. 2 * ' * 1 ( ) ( ) C F C F m δ σ β β = ∆ - ∆ - (12)
Kriging is an exact interpolation method, so predictions made at any point s give

ˆ( ) ( ) C C ∆ = ∆ s s .
Thus, the Kriging variance is null at these points.

Cross-validation of calibration model

Cross-validation is a model validation technique for assessing how the results of a statistical analysis generalize to an independent data set. It is mainly used to estimate how accurately a predictive model performs in practice. In order to better realize the calibration of a measurement system, the cross-validation technique is applied to select the best regression model and correlation model. [Insert Fig. 4

here]

The main steps are described as follows.

Step 1: Generate measurement error with three dimensions ( , , )

i i i i C x y z ∆ = ∆ ∆ ∆
by the N measured points and nominal points which are projected on the fitting surface obtained by LS method.

Step 2: Select M points as measured points for the estimation of calibration and temporarily eliminate the measurement error of the remaining (N-M) points.

Step 3: Calculate the estimation of calibration and Kriging variance associated in the remaining (N-M) points.

Step 4: Repeat q times by steps (2) to (3) for each of M.

Step 5: Calculate the cross-validation criterion.

Step 6: Build the box plot (box-whisker plot) for evaluation value of cross-validation criterion on the M points.

Once the cross-validation is performed, several important criterions are used to evaluate the performances of the different Kriging methods. Here is a list of the most popular cross-validation criterions.

mean error:

* 1 1 ( ) ( ) n i i i ME C s C s n =   = ∆ -∆   ∑ (13)
mean squared error:

2 * 1 1 ( ) ( ) n i i i MSE C s C s n =   = ∆ -∆   ∑ (14)
root mean squared error:

2 * 1 1 ( ) ( ) n i i i RMSE C s C s n =   = ∆ -∆   ∑ (15)
average Kriging standard error:

n i i i i C s C s MSPE n s δ σ = ∆ -∆ = ∑ (17)
root mean square standardized prediction error: The model is chosen if the cross-validation criterions satisfy: 1) Mean estimation errors (ME) and mean standardized prediction error (MSPE) are close to 0. This criterion ensures the absence of bias.

* 2 1 ( ) ( ) 1 [ ] ( ) n i i i i C s C s RMSPE n s δ σ = ∆ -∆ = ∑ (18)
2) The variance of the estimation errors (MSE, RMSE) is lower. This criterion reflects the robustness of the estimation and it provides the information on the estimation accuracy.

3) The variance of the standard errors (RMSPE) is closer to 1. It indicates that the standard deviation of Kriging properly reflects the estimation accuracy.

Applications

In this section, several experiments were conducted to illustrate how the proposed method based on Kriging calibrates a 3D measurement system and realize the measurement more precise. The calibration evaluation is analyzed with the comparison for ordinary Kriging and universal Kriging.

Experimental setup

The stereoscopic sensor (ATOS II Triple Scan, GOM) is installed on the manipulator (Fig. 5), which can be used to scan complex surfaces. The measurement system GOM is calibrated through the proposed calibration method based on a simplified phenomenological model of the sensor.

[Insert Fig. 5 here]

The experiment measurements were executed on the planes of gauge (see Fig. 6) and spheres of the gauge (see Fig. 7). To eliminate the influence of plane position, different positions of planes were measured by the sensor GOM: lower plane, middle plane, upper plane (blue), and diagonal plane (green) shown in Fig. 6.

[Insert Fig. 6 here] [Insert Fig. 7 here]

Calculate measurement uncertainty

In order to reduce systematic errors and realize calibration of the GOM measurement system, it is firstly important to calculate measurement error on the surfaces.

Plane case

To model the calibration of the GOM sensor on the datum plane, ∆z, or ∆z and ∆y are represented as the measurement error corresponding to the different positions of the plane.

For lower, middle, and upper planes inclined with small angle (less than 10°), ∆z is used to represent the calibrations, and for the diagonal plane inclined with large angle (around 45°), ∆z and ∆y are all used to represent the calibrations (see Fig. 8).

[Insert Fig. 8 here]

In Fig. 8, large dots are known as the measure, and the gauge plane is represented by the line. For the inclined plane with small angle, ∆z are the differences between the measured points and projection points on the gauge plane. For inclined plane with large angle, ∆z and ∆y are the differences between the measured points (black) and projection points on the gauge plane.

The

i i i i s x y z - - - =
were obtained by GOM sensor on the four planes. The gauge planes are shown in Table 1.

[Insert Tab. 1 here]

According to equation ( 9) and Fig. 8, the measurement errors ( , , ) The results of measurement errors on different planes are shown in Fig. 9.

i i i x y z ∆ ∆ ∆
[Insert Fig. 9 here]

Sphere case

The calibrated radius of the gauge sphere is 15.0004 mm, and the coordinate of the sphere center is fitted by LS method using the measured points (see Fig. 10). The theoretical reference sphere (spherical datum surface curve) can be obtained. The measurement errors are the differences between the coordinates of actual measured points and the coordinates of their projection points on the spherical datum surface, expressed by ( , , )

i i i
x y z ∆ ∆ ∆ . The gauge sphere is shown in Table 2. The results of measurement errors on spherical surfaces are shown in Fig. 11. The measurement uncertainty of the calibrated radius of the gauge sphere has an impact on the errors assessment and therefore on the kriging model. GOM gauge was used in this application. The accuracy of the radius form is 0.0012mm (GOM calibration certificate).

[Insert Fig. 10

Calibration results and discussion

Plane case

According to the proposed calibration method based on Kriging for GOM sensor, the measurement errors calculated at different sites on the plane are considered as regional variables, and are used to estimate the calibrations in any other sites. For the lower plane, 1000 three dimensional data and the corresponding measurement errors i z ∆ are randomly sampled from 33678 measured points. Two experiments were conducted to provide calibrations in various areas: on the gauge plane and in a larger cubic space around the gauge plane. The measured points are collected and used to estimate the calibrations in Fig. 12.

[Insert Fig. 12 here]

The process of calibration calculation is explained for the unmeasured points in Fig. 12.

For sampling in a larger cubic space, 1,000 points are collected by rectangular grid method in the larger cubic space around the fitting mean plane (Fig. 12 (a)). Then the calibrations of these 1000 points are predicted using previously calculated measurement error data i z ∆ selected from 33678 by random sampling method. While sampling on the fitting plane, 2912

(52 * 56) points are collected by mesh grid method on the fitting mean plane (Fig. 12 (b)).

Then the calibrations of these 2912 points are estimated using the same known measurement error data.

To estimate calibrations for sampling on the fitting plane, two methods are defined in the calculation (see Fig. 13). The first one is that the calibration is directly based on the positions of points 0 measured measured measured ( , , )

i i i z w x y z - - - ∆ =
(called direct method). The second one is based on Kriging described as follow. The height value of Z has spatial correlation with the positions X and Y. Thus, the estimation of the value Z is calculated by Kriging Method [Insert Fig. 14

here]

The distribution of the estimated calibration depends greatly on the distribution of the measurement error. If the point is closer to the expected measured points, the variance of the estimate is lower. The variance of the estimation in the cubic space is symmetrical along the gauge plane and lower around the gauge plane.

It is important to choose the appropriate Kriging method including correlation model and regression model corresponding to the sampling strategy in order to reduce the estimation variance. Table 3 shows that the results of estimation and variance in process of calibration.

According to Table 3, when the exponential model is chosen as the correlation model and Ordinary Kriging is chosen as the regression model in Fig. 14 [Insert Tab. 3 here]

Sphere case

Two experiments were conducted to provide calibrations in various areas: in the cubic space and on the gauge sphere shown in Fig. 15. The estimation points on the spherical datum surface are shown in Fig. 16, and results of calibration estimation and estimation error variance on the spherical datum surface are shown in Fig. 17. The calibration of GOM sensor can be achieved at any sites of the full-field area of measurement by the proposed Kriging method. For a new measurement mission, the calibration of the measurement system can be automatically compensated to the results of the measures.

[Insert Fig. 15 here] [Insert Fig. 16 here]

[Insert Fig. 17 here]

Sensitivity analysis of different numbers of points to construct meta-model

The cross-validation method not only helps to choose the appropriate Kriging method to estimate the calibration of the GOM sensor, but also assists to decide the number of measurement error applied to construct meta-model of calibration. The sensitivity analysis on different initial number of points to construct meta-model is conducted according to the process of cross-validation in Fig. 4, and the main steps are described as follows.

1) Generate the measurement error i z ∆ by the N measured points (N=33678).

2) Randomly select M=10, 50, 100, 300, 500 points from the N measured points for the estimation of calibration and temporarily eliminate the measurement error of the remaining points (N-M).

3) Calculate the estimation and Kriging variance associated in the remaining points. [Insert Fig. 18 here]

According to Fig. 18, some conclusions can be obtained:

1) The average estimation errors (ME) and standard errors (MSE) is closer to 0 with the increasing number.

2) The variance of estimation errors (MSE, RMSE) is lower with greater number of points. This criterion reflects the robustness of the estimator and information on the accuracy of the estimation.

3) The values of MSPE with different numbers of points converge to 0. But RMSPE has a fluctuation with the increasing number of points, and at last converges to the specific value, 1.

Thus, the accuracy of the estimation of calibrations is increased by the number of measured points for the creation of meta-model. But the accuracy is priced by measurement cost and the calculation time. Generally, the more number of points used to estimate the calibration by Kriging method, the higher accuracy of estimation is achieved. So it is important to choose the appropriate number of points for the meta-model creation. The best indication is the variability RMSPE.

• RMSPE = 1 indicates that the prediction of the calibration is correct and effective.

• RMSPE> 1 indicates that the prediction is overestimated, which means to reduce the measurement information.

• RMSPE <1 indicates that the prediction is underestimated, which means to increase the measurement information.

The variability of RMSPE by Ordinary Kriging and by Universal Kriging are shown in Fig. 19. For Ordinary Kriging, 50 or 100 measured points lead to a less than 1 RMSPE, but for 500 points it is greater than 1. This indicates that 500 measured points are sufficient for the prediction of calibrations. If more points are used for the calculation, the variability of the prediction will be overestimated. For Universal Kriging, 500 measured points are not enough to perform the calibration. For an industrial application, a reference threshold of RMSPE could be defined to optimize the needed number of measured points.

[Insert Fig. 19 here]

Spatial correlation analysis for measurement uncertainty

A convenient choice for the correlation function is within the power exponential model and spherical model. Spatial correlation functions are [Insert Fig. 20 here]

Conclusion

In this paper, a novel Kriging models-based calibration and uncertainty estimation method is developed to improve the measurement accuracy of non-contact 3D CMSs. This calibration procedure has been tested on a measurement system with the GOM sensor, which served as a tool of measurement after calibration. The results of calibrations are used to compensate and correct the measurement system to perform a new measure. The

  coordinates of a measured point respectively. These values are obtained from the simplified model of the sensor. Based on the geometrical characteristics of the gauge, the coordinate values of theoretical reference points are determined by an orthogonal projection of the measured points on the gauge surface. Each measured point is mapped into the corresponding theoretical reference point,

  the expectation of deterministic structure and ( ) δ ⋅ is a stationary random function with zero mean and known dependence structure.

  . It is a polynomial function of order 1 (linear universal Kriging, LUK) and of order 2 (quadratic universal Kriging, QUK) at the location s with three dimension x, y and z. And ordinary Kriging can be considered as a special universal Kriging with 1 p = and 1 ( ) 1 f = s .

  2

Fig. 4

 4 Fig.4shows the procedure of cross-validation. Several points of the measurement errors (testing data) are temporarily eliminated and then the calibration values of the measurement systems are estimated by Kriging methods using the remaining data (training data). This operation is repeated for some or all points. Thus at any observation point, it is calculated by

  standard variance value and n is the number of measured points used for the Kriging estimation.

∆

  are calculated by the difference between the measured coordination data and the coordination of the projected points on the gauge plane. For lower, middle and upper planes, are calculated as the measurement errors.

  predicted) (point projected on the fitting plane) z z z ∆ = -. For the diagonal plane, measurement error along Y axis is calculated in the same way. [Insert Fig. 13 here] Calibrations in both areas are estimated by the Kriging method, the calibrations 0 z ∆ in the cubic space and on the fitting plane are shown in Fig. 14. The left image shows the estimated calibration, and the right image is the variance of estimation.

  (a) and Fig.14 (b), and Universal Kriging is chosen in Fig.14 (c) respectively, the estimation variances 2

4 ) 5 )

 45 Repeat 50 times steps (2) to (3) for each M points. Calculate the cross-validation criterion 6) Build the box plot (box-whisker plot) of cross-validation criterion for evaluation.

Fig. 18 shows

 18 Fig. 18 shows the results of different cross-validation criterion ME, MSE, RMSE, AKSE, MSPE and RMSPE with different numbers of measured points for the creation of meta-model Kriging.

Fig. 20

 20 Fig.20shows the spatial correlation influenced by distance and θ . From Fig.20, it is

KrigingFig. 2 .Fig. 3 .Fig. 4 .

 234 figure[START_REF] Satoa | Calibration and uncertainty evaluation of single pitch deviation by multiple-measurement technique[END_REF], the accuracy of the results with the kriging compensation is 10 times the accuracy of the results without the kriging compensation.

] Step 1. Obtaining

  measured points and theoretical reference points of the gauge. A

	measured point S obtained by 3D CMS is denoted as,
	-measured	(	measured	,	measured	,	measured	)

Table 1

 1 The fitting planes by LS method

	Plane	Number of points	Gauge plane
	Lower plane	33678	z = -122.9953 + 0.00010 * x -0.0516 * y
	Middle plane	41051	z = 5.5013 + 0.001953 * x -0.05085 * y
	Upper plane	31894	z = 78.7114 + 0.003069 * x -0.04998 * y
	Diagonal plane	57649	z = 20.2790 -0.0125 * x -0.9290 * y

Table 2 .

 2 The fitting sphere by LS method

		X	Y	Z
	Gauge center of sphere	-85.2707	65.6985	-55.4355
	Measured radius		15.0088 mm
	Gauge radius		15.0004 mm

Table 3 .

 3 Results of estimation and variance in process of calibration

		In the cubic space	On the fitting plane by	On the fitting plane by
				Direct method	Kriging method
		Estimation Variance Estimation Variance Estimation Variance
	Mean			2.08E-04	5.50E-05	2.30E-04	4.99E-05
	Min	-7.61E-02 2.40E-05 -7.94E-02 4.88E-06 -8.31E-02 1.39E-06
	Max	9.03E-02	7.10E-04	9.92E-02	1.90E-04	1.08E-01	1.71E-04
	Regression		Ordinary Kriging		Universal Kriging
	model
	Correlation			Exponential model
	model


A Jebrane
email: aissam.jebrane@enpc.fr

P Argoul

A Hakim

M El Rhabi

Estimating contact forces and pressure in a dense crowd: microscopic and macroscopic models

Keywords: Crowd modeling, contact forces, non-smooth contact dynamics, non-local interactions, simultaneous collisions

This paper deals with the estimation of pressure at collisions times during the movement of a dense crowd. Through the non-smooth contact dynamics approach for rigid and deformable solids, proposed by Frémond and his collaborators, the value of pressure and contact forces at collisions points, generated through congestion or panic situation are estimated. Firstly, we propose a second-order microscopic model, in which the crowd is treated as a system of rigid solids. Contact forces are rigorously defined by taking into account multiple simultaneous contacts and the non-overlapping condition between pedestrians.

We show that for a dense crowd, percussions can be seen as contact forces.

Secondly, in order to overcome the restrictive hypothesis related to the geometric form adapted to model the pedestrian, a continuous equivalent approach is proposed where the crowd is modeled as a deformable solid, the pressure is then defined by the divergence of the stress tensor and calculated according to volume and surface constraints. This approach makes it possible to retain an admissible right-velocity, including both the non-local interactions between nonneighbor pedestrians and the choice of displacement strategy of each pedestrian.

Introduction

The modeling of crowd dynamics has been deeply studied in the 25 past years. Particularly, approaches drawn from granular modeling [START_REF] Pecol | The non-smooth view for contact dynamics by Michel Frémond extended to the modeling of crowd movements[END_REF][START_REF] Kabalan | A crowd movement model for pedestrian flow through bottlenecks[END_REF][START_REF] Maury | A macroscopic crowd motion model of gradient flow type[END_REF] offer both a rich mathematical and physical framework to simulate and comprehend the dynamics of a dense crowd. Predicting and anticipating the different scenarios, that can occur when a crowd moves, is a challenge, which must be taken into account for infrastructure planning (airport, concert hall, train station ...), it becomes necessary to avoid past tragic dramas of a dense crowd in panic, and an essential tool to achieve this goal is the prediction of the pressure and the contact forces generated during collisions of pedestrians in a moving crowd. Two major existing approaches that have dealt with this issue can be distinguished:

• A microscopic approach proposed in [START_REF] Maury | Handling congestion in crowd motion modeling[END_REF][START_REF] Maury | A macroscopic crowd motion model of gradient flow type[END_REF], where pedestrians are identified with rigid disks that must obey the non-overlapping constraint. The rightvelocity after collision is the projection of the left-velocity, corresponding to the desired velocity, over the set of admissible velocities with regard to the non-overlapping constraint. This constrained optimization problem is rewritten in a saddle-point form, and the corresponding Lagrange multipliers are interpreted as the pressure between pedestrians in contact. In this first-order microscopic model, pedestrian mass is not taken into account and the dimension of the Lagrange multipliers is incompatible with that of a pressure. To overcome these problems, we propose a second-order model mainly developed by Frémond [START_REF] Frémond | Collisions Engineering: Theory and Applications[END_REF][START_REF] Frémond | Rigid bodies collisions[END_REF] for modeling granular media.

The application of this approach for a crowd has been the subject of our previous work [START_REF] Kabalan | A crowd movement model for pedestrian flow through bottlenecks[END_REF][START_REF] Pécol | Smooth/non-smooth contact modeling of human crowds movement: numerical aspects and application to emergency evacuations[END_REF]. By adopting this approach, we will rigorously define the contact forces from the notion of percussion. A comparison with the first-order microscopic model will be presented. We cite also the social force model [START_REF] Moussad | How simple rules determine pedestrian behavior and crowd disasters[END_REF][START_REF] Helbing | Simulating dynamical features of escape panic[END_REF][START_REF] Lhner | On the modeling of pedestrian motion[END_REF], in which repulsive forces can be considered as contact forces. This force is introduced primarily as a penalty term to keep a certain distance and, consequently, prevent inter-penetration between pedestrians.

• Macroscopic approaches: Macroscopic models have been developed, to describe the evolution of the crowd as a continuum medium, it can be mentioned among them, a classical fluid approach [START_REF] Helbing | A fluid dynamic model for the movement of pedestrians[END_REF] where the motion of a crowd and the generated pressure are described by Navier-Stokes equations. Important approaches are developed in [START_REF] Golas | Continuum modeling of crowd turbulence[END_REF], [START_REF] Twarogowska | Macroscopic modeling and simulations of room evacuation[END_REF] and [START_REF] Venuti | An interpretative model of the pedestrian fundamental relation[END_REF], which the pressure is described by constitutive laws, which involves the speed and density of the crowd. A first order macroscopic model is developed by authors in [START_REF] Maury | A macroscopic crowd motion model of gradient flow type[END_REF], this approach makes it possible to take into account the non-local character of the pressure, however, it has been shown [START_REF] Maury | A macroscopic crowd motion model of gradient flow type[END_REF] that this model in its current version, reproduces a contradictory effect to that observed experimentally; congestion accelerates the evacuation of a convex room through a single door. Note also that several recently developed approaches like [START_REF] Taneja | Network redesign for efficient crowd flow and evacuation[END_REF] are using kinetic model [START_REF] Bellomo | On the modeling of crowd dynamics: An overview and research perspectives[END_REF] did not yet deal with the pressure issue.

Contribution and paper content:

In this work, our interest is devoted to the estimation of contact forces and pressure, in a moving dense crowd, the theory is inspired from the non-smooth approach developed by Frémond and his collaborators both for the evolution of rigid and deformable bodies. To fix the ideas, we consider the configuration of a dense crowd at a given moment (i.e. a stationary problem), and we ask how to estimate the contact forces or the generated pressure. Therefore, two classes of models will be developed: In section 2, a microscopic second-order model, in which the crowd is identified with a system of rigid bodies is presented.

In previous works [START_REF] Pécol | Smooth/non-smooth contact modeling of human crowds movement: numerical aspects and application to emergency evacuations[END_REF][START_REF] Pecol | The non-smooth view for contact dynamics by Michel Frémond extended to the modeling of crowd movements[END_REF], we have evoked this question, and we have arrived at an estimation of the percussion which is the integral of the contact force during the collision. Now, we propose to define the contact forces rigorously and with regard to physical proprieties. We start with a brief overview on the mathematical and physical framework of the non-smooth dynamics approach, particularly, the collision equation, then, we show that in the case of a dense crowd where the contacts remain maintained that percussions (jumps of the momentum) qualify as contact forces, ie, they will have the dimension of a force. As a validation test, the analytical solution of the contact equation is developed in the case of a constrained pedestrian chain to study the shock wave effect. The contact forces obtained by the proposed approach are compared with the results gived by a microscopic model of the first order, and also with those obtained by models of the social force type. In order to overcome the restrictive hypothesis related to the choice of the geometric form of a pedestrian ( rigid disk), and as we are interested in the case of a dense crowd, we model the crowd as a continuous medium, which is the object of section 3. This approach makes it possible to take into account the deformability property by modeling a dense crowd as a deformable solid. Contact forces are deduced from the volume and surface stress tensor which is defined by a constitutive laws. The proposed approach is compared to three other macroscopic approaches.

2.

A second-order microscopic approach: Crowd is modeled as a system of rigid bodies

The contact equation

We consider a crowd as a system of rigid bodies, using the theory of nonsmooth contact dynamics by Frémond [START_REF] Frémond | Collisions Engineering: Theory and Applications[END_REF][START_REF] Frémond | Rigid bodies collisions[END_REF], we obtain the contact equation given by:

M u + -u -+ ∂φ D (u + + u -) 2 + N Uq (u + ) 0, (1)
where:

• M is the inertia matrix.

• u + (resp u -) Right-velocity or so-called after impact velocity (resp leftvelocity).

• φ is a pseudo-potential of dissipation, φ = φ i,j , φ i,j (x) = k (x) 2 , k is a coefficient of dissipation. This choice allows one to find the classical results when the coefficient of restitution is used.

• U q the set of admissible velocities defined by:

U q = u ∈ R 2np , D i,j (u).n i,j q ≤ 0, 1 ≤ i < j ≤ n p , D i,j (u)
is the velocity of deformation corresponding to contact of two pedestrians p i and p j , it is given by D i,j (u) = u i (A) -u j (A). In 2D representation, the velocity of deformation at the the point A is D i,j (u) = u i (q i ) + w i × G i A -u j (q j )w j × G j A, where w i is the rotation velocity of pedestrian p i . Finally,

D(u) = D i,j (u), 1 ≤ i < j ≤ n p is a column vector of size n c , where n c
is the number of contact.

• N Uq (u) the normal cone of U q at point u.

Non local interactions are introduced by percussion p i,j , 1 ≤ i < j ≤ n p , with:

p i,j reac ∈ ∂1 R + D i,j (u + ).n i,j q ,
where n i,j q is the normal oriented from p i to p j , and 1 R+ (x) = 0 if x ≥ 0, +∞ otherwise.

Total percussion of reaction p reac = p i,j reac , 1 ≤ i < j ≤ n p has the form:

p reac ∈ ∂1 U u + = N U (u + )
This approach makes it possible to manage multiple and simultaneous contact and to take into account the corresponding local and non-local (at a distance) interactions. Moreover, the capacity to use different constitutive laws make it possible to account for diverse experimental results.

After reformulation of the the contact equation in the saddle point form, the solution u + and the percussions of reaction are characterized by:

   M (u + -u -) + kD T D (u + + u -) + C T λ = 0 Cu + .λ = 0
where C the matrix defining the constraints on the normal deformation velocities, and λ i,j the corresponding Lagrange multiplier to the contact between pedestrians p i and p j . We note that p i,j reac = λ i,j n i,j . For a dense crowd, contacts are maintained,the constraint of non-interpenetration is saturated, i.e. : Cu + = 0, which means that contact is purely inelastic, and there is no dissipation (k = 0). Consequently, λ is a solution of:

Wλ = Cu -, W = CM -1 C T , (2)
W is the Delassus matrix.

According to the collision equation Eq.( 1), the percussion has the dimension of momentum, and hence it is perfectly for an instantaneous collision. For a dense crowd which is the case we are interested in, it does not remain true, indeed: the contact is maintained and the duration is important, because of these considerations, the motion of the crowd is considered smooth and Eq.( 1) can be rewritten for a purely inelastic contact as a second-order spewing process [START_REF] Moreau | An introduction to Unilateral Dynamics[END_REF][START_REF] Dzonou | A sweeping process approach to inelastic contact problems with general inertia operators[END_REF]:

-M du dt ∈ N Ũq (u),
where Ũq is the moving set of admissible velocities. The contact forces f contact correspond to a selector in the outward normal cone, i.e. :

-M du dt = f contact , wheref contact ∈ N V ad (u).
In this formalism the acceleration has a sens and consequently the contact force f contact which approached numerically by C T λ.

Shock wave and non-local interactions

We rest. Here, the pedestrians are modeled by rigid disks see Fig. 1 , generated contact force between the blue and the neighbored red disks is given by the formulate 3, this shows that the movement of each pedestrian is influenced by the movement of others.

p 1,2 reac = -m 1 m 2 (u 2,--u 1,-) + m 3 (u 3,--u 1,-) + m 4 (u 4,--u 1,-) m 1 + m 2 + m 3 + m 4 , (3)
As one can see the percussion p 1,2 reac involves u 3,-and u 4,-besides u 1,-and u 2,-, we obtain the same for p 2,3 reac and p 3,4 reac . The effect of the shock wave is responsible for the non-local interaction ( or the so-called indirect interactions ) between non-neighboring pedestrians [START_REF] Glocker | On frictionless impact models in rigid-body systems[END_REF].

The second example, see Fig. 2, it shows the displacement of 53 pedestrians of the same mass m = 75kg, they are supposed initially at rest in a corridor after being pushed by three moving pedestrians. Fig. 2c illustrates the transmission of contact force by local and non-local interaction in a crowd.

Estimation of contact forces: Comparison of the proposed discrete approach

to other microscopic models 2.3.1. The proposed discrete approach versus the first-order microscopic model

In [START_REF] Maury | Handling congestion in crowd motion modeling[END_REF], contact forces are characterized by the following discrete Poisson equation: where CC T is considered as the discrete Laplacian matrix. According to this equation Eq 4, we see clearly that the authors assumed that all pedestrians have the same mass m = 1kg, this leads to the equilibrium situations and produces undesirable scenarios like the complete congestion [START_REF] Maury | Handling congestion in crowd motion modeling[END_REF], which is specially related to the geometry of pedestrians (disks) and their size (radius of disks). By choosing the matrix inertia equals to the identity matrix, the equation Eq 4 is a special case of the proposed approach given by Eq..

CC T λ = Cu -, (4)
In order to show that the hypothesis "Pedestrians have the same mass m = 1kg" made by the authors is restrictive, we start by a one-dimensional case of a pedestrian chain. Eqs.5 and 6 represent the analytical solution in the trivial case of an aligned constrained pedestrian chain as illustrated in Fig. 3a, they

show the dependence on the mass and the Lagrange multipliers (and hence the contact forces).

   λ 1,2 = -m 1 (u -u 1,-), λ i,i+1 = λ i-1,i -m i (u -u i,-), 2 ≤ i < n p . (5)
and

u i,+ = u = np i=1 m i u i,- np i=1 m i , 1 ≤ i ≤ n p , (6)
In addition, Fig. 3b shows the difference in the percussion suffered by each pedestrian, in the cases where pedestrians have identical and different masses. In the 2D simulation case, we are interested in estimating the contact forces generated during the evacuation of a room, we choose the instant correspond to the formation of the ark around the exit, see Fig. 4a. The number of pedestrians considered is 594 with the same mass m = 75Kg and a maximum density of 9 pedestrians per square meter. The contacts are assumed purely inelastic (i.e. k = 0), and pedestrians have a desired velocity towards the exit of a norm equals to 2m.s -1 . Right-velocity and the contact forces are given by:

u + = u -+ C T λ, (7)
where Wλ = Cu -.

Obtained results are given in Fig. 4: We notice that the force is maximum in the neighborhood of the exit Fig. 4b, where the movement of the crowd is rallied Fig. 4c. As we mentioned, the proposed approach is general, in a way that we can take different masses, to produce a realistic and different behavior according to the calibration of the coefficient of dissipation k.

Fig. 5 shows the obtained results in the case where pedestrians have different masses, the left-velocity is chosen uniformly, and yet, we notice that neither the speed on the right nor the pressure are uniform with respect to the distance to the exit.

The proposed approach versus social force model

In [START_REF] Helbing | Simulating dynamical features of escape panic[END_REF], the authors have defined contact force f ij between two pedestrians i and j by:

f ij = {A i exp[(r ij -d ij )/B i ] + g(r ij -d ij )} n ij + ς, g(r ij -d ij )∆v t ij t ij (9)
where:

• A i and B i are constants, A i = 2.10 3 N and B i = 0.08m. [START_REF] Helbing | Simulating dynamical features of escape panic[END_REF].

• r ij := r i + r j the sum of radius r i and r j of pedestrians i and j.

• d ij denotes the distance between the pedestrians' mass center.

• n ij is the normalized vector pointing from pedestrian j to pedestrian i, and t ij the tangential direction.

• The function g(x) is zero if the pedestrians are not in contact (

d ij ≥ r ij ),
and is otherwise equal to the argument x.

• The parameters = 1.2 × 10 5 kg.s -2 and ς = 2.4 × 10 5 kg.m -1 .s -1 determine the obstruction effects in case of physical interactions.

Let us take again the one-dimensional example of the constrained pedestrian chain Fig. 3a, according to the formula (9), if pedestrians are considered as rigid disks / spheres, than the contact force between each pair of pedestrians i and j is constant throughout the chain of equal intensity ||f i,j || = A i (because

d i,j = r i,j
), which of course is not realistic. Else if the pedestrians are modeled as a deformable disk, the contact forces should be calculated over the entire contact surface and not at a single point of a contact, as in formula [START_REF] Helbing | Simulating dynamical features of escape panic[END_REF].

In the two-dimensional case, we consider the same evacuation simulation of a room made by the authors in [START_REF] Helbing | Simulating dynamical features of escape panic[END_REF]. Fig. 6 clearly shows that there is no overlap between the pedestrians, which is, to say, still according to the formula [START_REF] Helbing | Simulating dynamical features of escape panic[END_REF], that the contact forces are constant independently of the masses, positions, desired velocities and also the density of the crowd. Therefore, the main difference between the proposed approach and the social force model, at the level of the estimate of the contact force, is due to the physical basis of each approach.

In the social force model, which is inspired by the smooth dynamics theory, contact forces are introduced as a penalty term in order to avoid the overlapping of individuals, Conversely to the proposed approach based on the non-smooth theory of contact dynamics, and rely on the notion of percussion which is the integral of the contact force during the collision. When the duration of the contact becomes important, the percussion converges towards the contact force [START_REF] Moreau | An introduction to Unilateral Dynamics[END_REF]. An alternative contact force is proposed by the authors in [START_REF] Lhner | On the modeling of pedestrian motion[END_REF], where the Figure 6: Evacuation of a room by authors in [START_REF] Helbing | Simulating dynamical features of escape panic[END_REF], where there is no overlap between pedestrians, which gives a constant contact force between all pedestrians in contact contact force f i is applied to a pedestrian i by another pedestrian j, their model is given by:

f i =    -fmax 1+d 2 i if d i,j ≤ r i,j , -2fmax 1+d 2 i if d i,j > r i,j , (10)
where: d i,j the distance between pedestrian i and pedestrian j,

d i = di,j
ri , f max is a maximum force such that fmax mi = O(max(α i , α j ), with m i the mass of the pedestrian i, and α i its acceleration. The formula shows that the force depends on the distance between the pedestrians, and that the contact force can be nonzero even if no contact takes place, contrary to the proposed approach where the contact force is active only at the instant of contact.

To conclude, the proposed approach based on the theory of rigid solid dynamics allows to manage the simultaneous contact in a crowd of pedestrians with different masses, and to estimate the corresponding contact forces which is rigorously defined from both a mathematical and physical point of view. It is also possible to show the effect of the shock wave responsible for non-local interactions.

Managing and treating multiple and simultaneous collisions (including the non interpenetration condition) in a dense crowd, where panic situations are very common, is the key point of the approach developed above. However, to model a pedestrian as a rigid disk stays a restrictive hypothesis. It is because of this problem, that we are focusing to treat in the following section by identifying the entire crowd as a deformable solid.

A macroscopic approach: Crowd is modeled as a deformable solid

The proposed continuous approach is also inspired by the non-smooth view of contact dynamics for deformable solids developed by Frémond [START_REF] Frémond | Rigid bodies collisions[END_REF]. As previously pointed out, we are interested in the problem of contact in a dense crowd, i.e. calculating the right-velocity and the generated pressure. Generally, macroscopic models consist on transporting the density, which represents the evolution of the crowd, by an admissible velocity, this latter can depend on several quantities such as the desired direction, overpopulated zones, flow direction, etc. We assume that all these considerations are taken into account by the left-velocity (velocity before the impact). We take inspiration from the basic approach established by Frémond to model the collision of a deformable solid with a rigid plane, the contact forces between pedestrians are given by the stress tensor, and pedestrians-obstacles contact forces are deduced from the constraints of non interpenetration.

The model: Mathematical and physical framework

We consider a dense moving crowd, whose density is denoted at instant t and point x by ρ(t, x), the left-velocity is u -(t c , x), where t c is a collision instant.

Applying the principle of virtual work on the crowd, we obtain at collision instant the motion equations given by:

ρ(t c , x) u + (t c , x) -u -(t c , x) = divσ(t c , x), x ∈ Ω, σ(t c , x)n(t c , x) = -r int (t c , x), x ∈ Γ w , σ(t c , x)n(t c , x) = r ext (t c , x), x ∈ Γ f , σ(t c , x)n(t c , x) = 0, x ∈ Γ e ,
where Ω is the area occupied by the crowd, whose border ∂Ω is defined by:

∂Ω = Γ w ∪ Γ f ∪ Γ e ,
where Γ w are the contact points with obstacles (walls, ...), Γ e the exits, and Γ f the free border that has no contact with any obstacle, see Fig. 7. Constitutive laws are needed to define the contact percussions r int , and the percussion stresses σ. They are defined respectively by a surface pseudo-potential of dissipation φ s ? that takes into account the impenetrability condition, and a volume pseudo-potential of dissipation φ v that defines the volumetric constraints:

φ s D u + (t c , x) + u -(t c , x) 2 = k s u + (t c , x) + u -(t c , x) 2 2
...

+I -u + (t c , x).n(t c , x) , φ v D u + (t c , x) + u -(t c , x) 2 = k v D u + (t c , x) + u -(t c , x) 2 2 , with D(u) = 1 2 ∂ui ∂xj + ∂uj ∂xi . r int (respectively σ) is the sub-differential of φ s (respectively of φ v ) at point u + (tc,x)+u -(tc,x) 2 : r int (t c , x) = k s u + (t c , x) + u -(t c , x) ... +∂I -u + (t c , x).n(t c , x) .n(t c , x), σ(t c , x) = k v D u + (t c , x) + u -(t c , x) .
Finally, the problem describing the crowd at the instant of collision is:

                     ρ(t c , x) [u(t c , x)] = k v div (D (u + (t c , x) + u -(t c , x))) , x ∈ Ω, σ(t c , x)n(t c , x) + k s (u + (t c , x) + u -(t c , x)) ... +∂I -(u + (t c , x).n(t c , x)) .n(t c , x) 0, x ∈ Γ w , σ(t c , x)n(t c , x) = r ext (t c , x), x ∈ Γ f , σ(t c , x)n(t c , x) = 0, x ∈ Γ e . . ( 11
) [u(t c , x)] := (u + (t c , x) -u -(t c , x
)) the jump of velocity discontinuity. For the existence and uniqueness of the solution of the problem described by Eqs. [START_REF] Helbing | A fluid dynamic model for the movement of pedestrians[END_REF], it is shown in [START_REF] Frémond | Collisions Engineering: Theory and Applications[END_REF] that at the collision instant, if the domain Ω is regular, with mes(Γ w ) > 0, and

ρ(t c , x) ≥ ρ 0 > 0, if u -(t c , .) ∈ H 1 (Ω), r ext (t c , .) ∈ L 2 (Ω)
then there is an unique right-velocity u + (t c , .) ∈ H 1 (Ω) solution of the problem (11).

Numerical simulations and comparisons

3.2.1. Validation tests: Jostling in a pedestrians chain.

We take the one-dimensional example of the pedestrians chain introduced previously. For a constrained pedestrian chain, i.e. a fixed obstacle is placed in front of pedestrians at the point x = l, see Fig. 3a. The problem [START_REF] Helbing | A fluid dynamic model for the movement of pedestrians[END_REF] becomes :

         ρ(t c , x) (u + (t c , x) -u -(t c , x)) = k v ∂ 2 (u + (tc,x)+u -(tc,x) ∂x 2 , x ∈]0, l[, -k v ∂(u + (tc,l)+u -(tc,l)) ∂x + k s (u + (t c , l) + u -(t c , l)) + ∂I -(u + (l)) 0. ∂(u + (tc,0)+u -(tc,0)) ∂x = 0. (12)
For a constant left-velocity and an uniform density, under the following assumption:

k ≥ e -2l/h -1 e -2l/h + 1 , k = k s k v , 1 h 2 = ρ k v , (13)
the solution of the unidimensional problem ( 12) is:

u + (t c , x) = u -(t c , x)(1 - e x/h + e -x/h e l/h + e -l/h , with h = k ρ ,
and internal percussion is given by ( 14), it becomes maximal as x goes to l.

r int (t c , x) = u -(t c , x) e x/h -e -x/h h e l/h + e -l/h (14)
The condition 13 implies that u + (l) = 0, i.e., the right pedestrian remains crushed against the obstacle. However, the other pedestrians stay glued to each other after the impact . If the condition ( 13) is not verified, rebounds will appear and the solution must be changed. Note that for an unconstrained pedestrians chain with a constant left-velocity, the right-velocity after collision is constant and the pressure is equal to zero. However, for an unconstrained pedestrian chain with a decreasing left-velocity, we get the same result as in the rigid case.

By analogy to the constrained pedestrian chain in the second-order microscopic model, we evaluate the pressure for a uniform density and a variable density.

For this, it is quite enough to mark that the percussion given by equation 14 is an increasing function, with respect to the density ρ. This is exactly the same result as in microscopic case.

The proposed macroscopic approach versus continuous approach for modeling crowd turbulence

In [START_REF] Lhner | On the modeling of pedestrian motion[END_REF], inter-pedestrians stress is given by a symmetric 2 × 2 tensor field representing the internal forces in the crowd:

σ = pI + S (15)
where the normal component p represents the normal or repulsive forces, while the deviatoric stress S represents the frictional forces. The pressure p is defined by an increasing density function g :

g(ρ) =          0 if ρ ≤ ρ c c 1 ρ √ ρy-ρ+ -c 2 if ρ c < ρ < ρ y c 3 ρ -c 2 if ρ ≥ ρ y . (16)
where the values of the parameters ρ c , ρ y , c 1 , c 2 , c 3 and are given in the table 1.

The value of frictional stress S is determined by minimizing the kinetic energy

E: E = 1 ρ y ||ρ(v -∇p) + δt∇.S|| 2 dV ( 17
)
subjected to the following constraint:

||S|| F ≤ √ 3βg(ρ) (18)
with ||.|| F is the Frobenius norm, and α the coefficient of friction.

The contact forces are given by:

∇.σ = ∇p + ∇.S (19)
To compare this approach with our proposed macroscopic model, we consider the example of an arch formation during an emergency evacuation of a room, previously introduced in the discrete case. The room is a square with exit side of 1m. Desired velocity u d is defined by: u d = u d e d , where u d is the free speed and e d the desired direction, e d (x) =

x-PΓ e (x)

||x-PΓ e (x)|| , where P Γe (x) is the projection of point x on the exit. Fig. 8 shows the error corresponding to the difference between the contact force intensity, found by the authors in [START_REF] Lhner | On the modeling of pedestrian motion[END_REF], and our proposed approach. We note the following remarks:

• We obtained similar results for k v >= 4 and a time step δt = 0.01s (error of the order of 10 -2). As in the discrete case, the dissipation coefficient reflects the aggressive nature of pedestrians [START_REF] Kabalan | A crowd movement model for pedestrian flow through bottlenecks[END_REF]. • In both approaches, the estimation of contact forces is inspired by the modeling of granular media. Authors in [START_REF] Lhner | On the modeling of pedestrian motion[END_REF] proposed a constitutive law allowing to define the density-pressure dependence, and frictional forces minimizing kinetic energy over a short period. In the proposed approach, all these considerations have been taken into account by the pseudo-potential of dissipation that defines the percussions tensor.

• The contact forces must be calculated independently to the time discretization step, this is not possible within the framework of the smooth dynamics that is the basis of the approach in [START_REF] Lhner | On the modeling of pedestrian motion[END_REF].

• The dependence of the kinetic energy by the step of discretization in time gives different values of frictional force according to the time step.

The proposed approach versus a first order macroscopic model

In [START_REF] Maury | A macroscopic crowd motion model of gradient flow type[END_REF], the authors proposed a first order continuous model based on the mass conservation, in which the pressure is characterized by the following place equation:

         -∆p = -∇.u d in Ω ∂p ∂n = 0 on Γ w p = 0 on Γ e ∪ Γ f (20)
Fig. 9 shows the results obtained by the proposed approach and the first order model for a uniform desired velocity u d = 1.7m/s and a uniform density ρ = 7

pedestrians per square meter. The proposed approach gives a Gaussian form of pressure concentrated around the exit, this result reflects the accumulation of forces near the exit. However, in the case of the first order macroscopic model, the pressure is minimal around the exit. Consequently, as the authors in [START_REF] Maury | A macroscopic crowd motion model of gradient flow type[END_REF] confirm, the velocity is maximum at these zones, which is contradictory with reality.

The proposed approach versus constitutive laws

In this part, we will compare the variation of pressure and the right-velocity with respect to the density in the case of our approach, with constitutive laws The complete congestion is undesirable, as it can be seen in actual emergency evacuation situations, this is consistent with the result obtained by the proposed approach Fig. 10b (blue line), unlike the two constitutive laws Fig. 10b (black and red line).

Figs.11a and 11b represent respectively the right-velocity and the corresponding pressure for a given uniform left-velocity and initial density. Figs. 11c and11d correspond respectively, to the right-velocity and the generated pressure for a given random left-velocity.

According to Fig11, the right-velocity is canceled at the boundary of the domain which reflects the condition of non-interpenetration with the wall, and the pressure is maximal around the exit, as for the microscopic case. Particularly, the general form of the right-velocity depends on the left-velocity, for (a) and (b) where the left-velocity is chosen uniformly with respect to the distance to the exit, the right-velocity keeps the same form. For a random left-velocity, we find a random right-velocity . We recall here for this example, the density is chosen uniformly, the speed also depends on the density.

Conclusion

In this paper, we are interested in estimating contact forces generated in a dense crowd from both discrete and continuous point of view. The mathematical and physical framework is inspired by the non-smooth contact dynamics theory developed by Frémond and his collaborators for rigid and deformable bodies assemblies. Starting with the microscopic approach, it has been shown the ability of this approach to manage multiple and simultaneous contacts, and to estimate the corresponding contact force via the notion of percussion. The case of arch formation is studied, the obtained results and the comparative study with other contact force approximation techniques show the superiority of the proposed approach. The pedestrian choice form remains the main criticism of all discrete models, different values can be obtained depending on the pedestrian representation and the type of contact. To overcome this challenge, a continuous version of the discrete approach is introduced. It make it possible to keep highlights of the discrete case and to estimate rigorously the contact forces by the stress tensor. It has been shown also that, for a good choice of the values of modeling parameters, we find results more significant and realistic than those found in the literature. As perspectives, the performance of the continuous approach is considered. In particular, taking into account the condition of noninterpenetration between pedestrians as in the discrete case and the study of the convergence of the microscopic approach towards the macroscopic approach.

Exploring the case where each pedestrian is modeled as a deformable solid also holds our attention in future work.

Acknowledgment:

The first author would like to thank Anouar BEN-LOGHFYRY, for his help editing the English in this paper.

Figure 1 :

 1 Figure 1: Simultaneous contact of four aligned disks. (a) The initial state, (b) The state after impact for an inelastic contact (k=0), (c) the state after impact for an elastic contact (k=10)

Figure 2 :

 2 Figure 2: Effect of the shock wave for a group of 56 pedestrians in a corridor, blue discs: pedestrians supposed to rest at the initial state while the three pedestrians in red push with an initial velocity v = (2, 2). (a) initial stat, (b) final state, (c) schematization of local and non-local interactions in the discrete model.

Figure 3 :

 3 Figure 3: (a) Constrained pedestrian chain, (b) the generated percussion, red line: pedestrians have the same masses, blue line: pedestrians have different masses.

Figure 4 :

 4 Figure 4: Estimation of contact forces in the situation of an arch formation with pedestrians of the same masses. (a) Initial configuration as an arch, (b) Right-velocity, (c) Norm of the corresponding contact force applied on each pedestrian, (d)Influence of the dissipation coefficient k on the contact force.

Figure 5 :

 5 Figure 5: Estimation of contact forces in the situation of an arch formation with pedestrians of different masses. (a) Initial configuration as an arch, (b) Right-velocity, (c) Norm of the corresponding contact force applied on each pedestrian.

Figure 7 :

 7 Figure 7: 2D domain of evolution.

Figure 8 :

 8 Figure 8: The norm of the difference between the contact force modulus calculated by the authors [10] and that by the proposed approach.(a) Error corresponding to a random desired velocity of an average 1.7m/s, standard deviation 0.5m/s and an uniform density equals to 7 pedestrians per square meter.(b) Error for an uniform desired velocity equals to 1.23 m : s and density equals to 7 pedestrians per square meter

Figure 9 :

 9 Figure 9: pressure estimation in the case of arch formation, (a) the proposed approach, (b) the first order macroscopic model, (c) the relative error.

Figure 10 :

 10 Figure 10: (a) Pressure versus density, Red line (M1), black line (M2) blue line (proposed approach), (b) Density versus velocity, red line (L1), black line (L2) blue line (proposed approach).

Figure 11 :

 11 Figure 11: Right-velocity (a) and the correspond pressure (b) for an uniform leftvelocity. (c) right-velocity and (d) pressure for a random left-velocity of average 1.34m/s and standard deviation of 0.26m/s.

Table 1 :

 1 Parameter values used in the numerical simulations

	Parameter name	Symbol Value	Units
	volume dissipation coefficient	k v	9	N.s 2 .m -2
	surface dissipation coefficient	k s	1	N.s 2 .m -3
	desired speed	u max	2	m.s -1
	maximal density	ρ max	8.5	ped.m -2
	density-speed coefficient	α	7.5	-
	adiabatic exponent	γ	2	-
	pressure coefficient	p 0	0.005 ped 1-γ .m 2+γ .s -2
	measurement parameter	R	0.7	m
	yield density	ρ y	5.5	ped/m 2
	critical density	ρ c	3.5	ped/m 2
	epsilon		0.01	-
	coefficient 1	c 1	13.33	-
	coefficient 2	c 2	1.5	-
	coefficient 3	c 3	c1 √	-
	frictional coefficient	β	0.75	-
	step time	δt	0.01	s

describing the density-pressure and the density-speed dependence, Fig. 10. Fig. 10a shows pressure as a function of the right-velocity.

We are interested in tracing the behavior of the density-pressure and densityright-velocity, then compare it with others models. For the dependence densitypressure, we compare our result with two models. The first one is proposed in [START_REF] Lhner | On the modeling of pedestrian motion[END_REF] and it is also restrained by the authors in [START_REF] Golas | Continuum modeling of crowd turbulence[END_REF]. Pressure-density relation is given by:

where ρ(x, t) indicates the density at location x, Var(V) is the current velocity variance, and

is the distance of location x and pedestrian i, and f

, where R is a measurement parameter. In the second model [START_REF] Twarogowska | Macroscopic modeling and simulations of room evacuation[END_REF], the density-pressure relation is defined by the following constitutive law:

where p 0 is the pressure coefficient and γ the adiabatic exponent.

Concerning the variation of the speed as a function of the density, we compare our result with two constitutive laws, the first is adapted by authors in [START_REF] Twarogowska | Macroscopic modeling and simulations of room evacuation[END_REF],

which is given by:

, the second is the Kladek formula introduced by Weidmann [START_REF] Weidmann | Transporttechnik der fussgänger[END_REF], it is adopted in the revisited version proposed in [START_REF] Venuti | An interpretative model of the pedestrian fundamental relation[END_REF]:

where ρ max is the congestion density, u max is the maximum walking speed in unconstrained conditions (also called free speed), and γ is an exponent.

The simulation data are given in the following table 1, note that we kept the same values of the parameters appearing in the formulas (M 1 ), (M 2 ), (L 1 ) and (L 2 ), that ae given in [START_REF] Golas | Continuum modeling of crowd turbulence[END_REF], [START_REF] Twarogowska | Macroscopic modeling and simulations of room evacuation[END_REF] and [START_REF] Venuti | An interpretative model of the pedestrian fundamental relation[END_REF].