Wojciech Aniszewski
email: aniszewski@dalembert.upmc.fr

Stéphane Zaleski

Antoine Llor
email: antoine.llor@cea.fr

Leon Malan

Numerical simulations of pore isolation and competition in idealized micro-spall process

Keywords: micro spall, cavity, Volume of Fluid (VOF), free surface, pore competition

Introduction

Cavitation and micro-spall (see for instance [START_REF] Signor | Investigation of fragments size resulting from dynamic fragmentation in melted state of laser shock-loaded tin[END_REF], and references therein) appear when a weakly compressible (or expansible) liquid is suddenly submitted to a large volume growth (expansion) resulting in negative pressures. This process, which appears in many practical applications of interest, occurs when an initially solid medium is subject to an unsustained impact or detonation and is entirely or partly melted in the process. Once the system starts expanding, the pressure drop causes the onset of cavitation and pores (or bubbles filled with vapor) appear. After the initial phases of uniform expansion and pore opening, a longer-lasting phase of pore growth and competition appears, which is especially difficult to investigate either experimentally or numerically [START_REF] De Rességuier | Dynamic fragmentation of laser shock-melted tin: experiment and modelling[END_REF][START_REF] Signor | Investigation of fragments size resulting from dynamic fragmentation in melted state of laser shock-loaded tin[END_REF].

We present here numerical simulations of this latter phase for idealized conditions relevant to micro-spall. This paper picks up on a previously published study of [START_REF] Malan | Direct numerical simulations of pore competition in idealized micro-spall using the VOF method[END_REF] which presented low-Weber number expansions of the system (with constant expansion rates). [START_REF] Malan | Direct numerical simulations of pore competition in idealized micro-spall using the VOF method[END_REF] focused on a competition phenomenon in which expansion of some of the pores caused their neighboring pores to collapse: in other words, there was a volume transfer between the pores without actual contact and merger between them. The pore competition effect is important as it is the main phenomenon driving the evolution in time of the statistical distribution of pore sizes [START_REF] Everitt | Competition and interaction of polydisperse bubbles in polymer foams[END_REF]. Small pores shrink and eventually disappear as their volume is transferred to large pores. Pore statistics and pressure evolution profiles can then be obtained for future modeling purposes. Hypothetically, continued competition accompanied by coalescence could lead to formation of larger pores-or even of a single pore in the case of a spatially limited system. In this work, we focus on the practically relevant case of constant speed expansions whose rate decrease with time: the corresponding Weber numbers typically start and end respectively at ranges above and below the competition threshold.

Thus, our idealized assumptions are incompressible inviscid fluid, vanishing vapor pressure in cavities, homogeneous uniaxial "ballistic" expansion, perturbed face-centered-cubic lattice arrangement of pores. Under these assumptions, the system is characterized by a single dimensionless group, the Weber number based on the number of pores per unit volume. The "ballistic" expansion signifies evolution in which We rapidly decays: in the initial stages of expansion the flow is entirely dominated by inertia and the pores expand as if they were isolated [START_REF] Ilinskii | Bubble interaction dynamics in Lagrangian and Hamiltonian mechanics[END_REF]. As We drops the system transitions into a regime in which the bubbles interact and capillary pressure becomes significant. We investigate this transition in more detail below. Results indicate that both the uniform growth in the initial expansion phase, and the later "competition" [START_REF] Everitt | Competition and interaction of polydisperse bubbles in polymer foams[END_REF] regime can be captured.

Elementary characterization of ideal micro-spall

Mean expansion field

We investigate a fluid element undergoing expansion following a shock and release history. As sketched in the volume and pressure profiles of Fig. 1, pressure drops due to the accelerated As qualitatively represented, the initial pressure produced by the shock and the final spall stress are here assumed respectively large and small enough that the eventual expansion is nearly ballistic. This ballistic evolution can then be backward extrapolated to vanishing volume which defines the (arbitrary) time origin. As indicated, the system evolution is here simulated after pores have opened and local fluid density has returned to a nearly uniform and constant value.

expansion and eventually becomes negative. At some point, the cohesive limit of the fluid is reached and pores open: in this micro-spall phase, surface tension around the pores produces a macroscopic negative stress. If this spall stress is weak enough compared to the stored momentum the system eventually evolves in a quasi-ballistic way. The micro-spall phenomenon can appear regardless of the structure of the strain tensor-which in general displays three different eigenvaluesbut we shall focus here on the most common situation of uniaxial expansion and more marginally on the theoretically useful situation of isotropic expansion. For a small enough material volume, the average strain can be considered as uniform -which corresponds e.g. to uniform laser energy distribution in laser shock-loading (Signor et al., 2010, sect. 2.2) -and the backward-extrapolated positions of all the fluid elements collapse to a single point: this is the virtual "Big Bang" that provides the time and space origins in all the following, t = 0 and x = 0. Depending on the strain dimension and depending on the selected coordinate system, Lagrangian X or Eulerian x, the mean velocity field in the material volume can always be written as

Lagrangian Eulerian Uniaxial u(X, t) = u(X, t 0 ) = X/t 0 , u(x, t) = x/t, (1a) Isotropic u(X, t) = u(X, t 0 ) = X/t 0 , u(x, t) = x/t, (1b)
where the Lagrangian coordinate X is the actual position at reference time t 0 , while x and X represent the first coordinates of the x and X position vectors. Equations (1) describe expansion, as the "Hubble time" t 0 is assumed positive. The uniform but non-constant divergence of the mean velocity field is

∇ • u = d/t, (2)
where dimensionality d is 1 or 3 for respectively uniaxial or isotropic expansions, and is related to mean density through

d d t ln ρ = -∇ • u. (3)
Combining ( 2) and (3) yields

ρ = ρ 0 (t 0 /t) d , (4)
where the reference state is conveniently taken here at time t = t 0 where cavitation nuclei appear [START_REF] Batchelor | An introduction to fluid dynamics[END_REF] and where the mean density coincides with the fluid density ρ 0 = ρ l .

Equations ( 1) are solutions of the Euler equation

∂ t u + (u • ∇)u = - 1 ρ ∇p, (5)
with zero mean pressure gradient and uniform but non-constant mean density ρ. The ultimate goal of modeling is to provide a closure relationship between ρ, u, and p so as to simulate the effect of micro-spall on large scale dynamics (possibly elliptically unstable).

Mean pore characteristics

At microscopic scales, micro-spall displays fluctuations on all its smooth large-scale fields. As sketched in Fig. 2, a first approximation convenient to capture surface tension effects consists in assuming all pores as spherical, with no internal mass and pressure (vanishing vapor pressure), and furthermore with identical radii R(t) and no translation velocity with respect to their surroundings [START_REF] Caflisch | Effective equations for wave propagation in bubbly liquids[END_REF]. Therefore, no coalescence effects can be present, and the number of pores is constant and equal to the number of cavitation nuclei at t = t 0 . The mean number density of pores N then scales as mass density (4)

N = N 0 (ρ/ρ l ) = N 0 (t 0 /t) d , (6)
where N 0 is an initial number of cavitation nuclei. In a more realistic situation where pores may collapse (depicted in this work) or merge, this relationship does not hold.

Under the simplifying assumption of local spherical symmetry the number density of pores provides their geometric characteristics. Each pore influences a mean spherical volume defined by the Wigner-Seitz radius (t) [START_REF] Girifalco | Statistical Mechanics of Solids[END_REF]:

4π 3 3 N = 1. (7)
Seeing as the volume fraction of liquid is (ρ/ρ l ), identification of volume fractions in the Wigner-Seitz sphere yields the mean pore radius R(t) from the equivalent relationships

4π 3 R 3 N = 1 -(ρ/ρ l ) = 1 -(t 0 /t) d , (8a) 4π 3 R 3 N 0 = (ρ l /ρ) -1 = (t/t 0 ) d -1. ( 8b
)
Combining ( 8) and ( 7) yields the aspect ratio of the equivalent mean spherical shell

(R/ ) 3 = 1 -(ρ/ρ l ) = 1 -(t 0 /t) d . (9)
Equations ( 8) to ( 9) hold even if, due to collapse and merger, (6) may not.

Kinetic and surface energies, Weber number

Under expansion, the kinetic and surface energies evolve and the overall behavior of the system is characterized by their ratio: the Weber number. Because the fluid is assumed incompressible and inviscid, these two energies fully control the system and only two regimes are expected: all situations can be collapsed according to the Weber number alone. The energies can be estimated within the framework of the mean spherical pore of Section 2.2.

The surface energy of a single mean pore is given by the surface tension coefficient σ

E s = 4πR 2 σ. ( 10
)
The kinetic energy is obtained by assuming the velocity field to be radial and divergence free around the mean pore, with a vanishing mass flux at the Wigner-Seitz boundary

E k = R 1 2 ρ l ˚ 2 r 4 4πr 2 d r = 2πρ l d/3 t 2 6 1 R - 1 , ( 11
)
where ˚ is the time derivative of which according to ( 7) and ( 6) verifies ˚ / = (d/3)/t. Using identities (6) to (9) we can now represent both mean radius R and Wigner-Seitz radius only in terms of t and state of the system at t 0 . Namely, from (6) we have

N 0 t 0 t d = 3 4π 3 , (12)
thus can be expressed as:

l = 3 3 4πN 0 t 0 t d . (13)
Using the aspect ratio (8) we find that Subsequently, we substitute both ( 13) and ( 14) in ( 11), so the Weber number We = E k /E s is eventually obtained as

R = 3 1 - t 0 t d . ( 14
)
We(t) = We 0 d 2 (t 0 /t) 2-d 3 1 -(t 0 /t) d 1 - 3 1 -(t 0 /t) d , (15a)
We

0 = ρ l 8πσN 0 t 2 0 . (15b)
Notice here that the Weber scaling We 0 built from quantities at t = t 0 is not equal to We(t 0 ): it is actually found that We(t) = We 0 at t/t 0 ≈ 1.16 and 1.38 for respectively d = 1 and 3with volume expansion ratios of ρ l /ρ = (t/t 0 ) d ≈ 1.16 and 2.60.

As illustrated in Fig. 3, the Weber number We(t) diverges as (tt 0 ) -1 for t → t 0 and vanishes as (tt 0 ) -2 for t → ∞.

The expressions in (15) were scaled in such a way that these asymptotes actually intersect at We 0 for either t/t 0 = 4/3 or 2 for uniaxial and isotropic expansions as shown in Fig. 3. The fact that the Weber number We(t) decreases from virtually unbounded values down to zero shows that two regimes should appear during expansion. Kinetic energy dominates at early times and expansion is quasi-ballistic, whereas later, surface tension dominates, thus slowing expansion and inducing pore interactions through Laplace pressure. The transition will appear more or less early depending on the initial conditions defined by the Weber scaling We 0 and possibly in either of the t → t 0 or (tt 0 ) -2 regimes represented in Fig. 3. According to the N 0 dependence in (15b), these two regimes correspond to respectively dense or sparse distributions of cavitation nuclei.

Pressure Evolution

Mean macroscopic p evolution has been sketched in Figure 1, presenting its sharp drop after the shock and during accelerated expansion phase. At the moment of pore opening pressure changes sign, while its growth is restored in the micro-spall regime. We can trace this evolution in a more detailed manner using velocity potential φ(t) and mean radii R. We first formulate the Bernoulli equation

φ t + u 2 2 + p -p ∞ ρ = 0, (16)
where u = ṘR2 /r 2 and thus φ = -ṘR 2 /r. This leads to a Rayleigh-Plesset type expression for pressure "at infinity"

p ∞ = -ρ l RR + 3 2 Ṙ2 - 2σ R . ( 17
)
Having discussed the mean spherical cell aspect ratio in the previous subsection, we can now solve for R using (8) leading to

R = - 2l 0 9 3 √ t (t -t 0 ) -5/3 , (18)
which allows us to retrieve p ∞ by substituting back to (17):

p ∞ (t) = ρ l l 2 0 18 3 t 2 0 (t -t 0 ) 4 - 2σ ρ l 3 t 0 t -t 0 . ( 19
)
It is easy to verify that the above formula depicts a curve whose shape corresponds closely to the pressure curve sketched in Figure 1 as long as t 0 l 0 . The latter is normally the case since mean Wigner-Seitz radii will be significant at the moment of pore opening.

While the simulations described in this work are set up in such a way that initial time value is t 1 > t 0 , it is interesting to study the effect of the propagating pressure pulse on the system. This reasoning is not without substance, as in computational practice the flow will respond to imposed expansion boundary conditions used here with a negative pressure pulse (see e.g. Figure 7 in [START_REF] Malan | Direct numerical simulations of pore competition in idealized micro-spall using the VOF method[END_REF] and results presented in Section 5). We assume incompressibility, however because of the negative pressure created by surface tension on the bubbles/pores, the system could be seen as having negative compressibility, so that the speed of sound c, given below:

c 2 = ∂p ∂ρ (20)
is imaginary. Taking this into account, and linearizing the Euler equation ( 5) and the mass conservation equation ( 3) we obtain

|c| 2 ∂ 2 xx p + ∂ 2 tt p = 0 (21)
A quick estimate based on Laplace's law leads to

|c| 2 = σ 4πNρ l R 4 . (22)
A conclusion that may be drawn from the pressure equation ( 21) is that as soon as the pressure has been affected by growing pores in a region, the pressure pulse will propagate to the rest of the system at a speed |c|. Such a propagating solution, of the form p = atan(x/|c|t), is possible even in the elliptic system case. The speed |c| of the wave may be comparable to the initial speed L 0 /(2t 1 ) of the expansion wave; comparison will be dependent on surface tension and Weber number.

Computational Methods

The PARIS-Simulator

The simulations presented in this paper use the Parallel Robust Interface Simulator (PARIS-Simulator [START_REF] Ling | Multiscale simulation of atomization with small droplets represented by a Lagrangian point-particle model[END_REF][START_REF] Malan | Direct numerical simulations of pore competition in idealized micro-spall using the VOF method[END_REF][START_REF] Salvador | Analysis on the effects of turbulent inflow conditions on spray primary atomization in the near-field by direct numerical simulation[END_REF], which is an in-home CFD code developed jointly at Institute ∂'Alembert, University of Notre-Dame and University of Bologna. PARIS-Simulator is a classical, MAC-type solver using uniform, cuboidal meshes. The strengths of the code lie in a very efficient MPI parallelisation and broad variety of implemented computational methods, especially for interfacial, two-phase flows. These include e.g. Front Tracking1 [START_REF] Tryggvason | Direct numerical simulations of gas-liquid multiphase flows[END_REF], Volume of Fluid [START_REF] Hirth | Volume of Fluid (VOF) method for the dynamics of free boundaries[END_REF] and tracking Lagrangian particles. The code is GPL licensed 2 and publicly available3 .

A well known projection scheme [START_REF] Tryggvason | Direct numerical simulations of gas-liquid multiphase flows[END_REF]) is used to solve momentum conservation equations. We will explain it briefly by showing a way the algorithm progresses from n-th to n + 1-st step. Knowing the values of all fields at the end of n-th step, we are starting with the definition of the temporary velocity field u

* u * -u n ∆t = -u n • ∇ h u n . ( 23
)
It can be found easily since, at the end of n-th time-step it is the only unknown in above formula. The symbol ∇ h stands for the discrete differential operator. This is a projection step, since velocity is projected onto a space with zero pressure field.

We can now write a discrete version of ( 5), which involves u * , only this time changing the way the approximation of temporal derivative is formulated

u n+1 -u * ∆t = - ∇ h p n+1 + σκnδ s ρ n , ( 24
)
where superscript n + 1 stands for the value at the end of nth time-step. Notice surface tension contribution has been included in (24): σ stands for the surface tension coefficient, n is the vector normal to the interface, while κ is scalar curvaturealso the restriction to the interfacial surface S is ensured by the δ S operator. If now we apply a divergence operator to both sides of ( 24) and remember that we assume ∇ • u n+1 to vanish, we obtain the Poisson equation from which pressure p n+1 can be found

∇ h • ∆t ρ n ∇ h p n+1 = ∇ h • u * + ∇ h • ∆t ρ n σκnδ s . (25)
In most computational codes numerical solution of ( 25) is the stage to which most computational cost (up to 90 percent) is associated; especially in multi-phase flows with variable density. Once ( 25) is solved and p n+1 found, both it and u * are used to find the divergence-free velocity at the end of the time-step

u n+1 = u * - ∆t ρ n ∇ h p n+1 + σκnδ s (26)
thus ending the procedure. Many authors call (26) a "correction" step as solenoidal character of u n+1 is ensured thereby.

Interface Tracking

As mentioned, the normal vectors n and the curvatures κ in above equations 4 are found from the color/fraction function C, whose advection equation follows in discrete form

C n+1 -C n ∆t + ∇ h • (Cu n ) = 0. ( 27
)
Equation ( 27) cannot be solved directly, as C is a sharp jump function: the jump would be diffused by numerical errors [START_REF] Youngs | Numerical simulation of turbulent mixing by Rayleigh-Taylor instability[END_REF]. Therefore, in most applications ( 27) is solved using specially crafted geometrical reconstruction-advection schemes such as CIAM [START_REF] Li | Calcul d'interface affine par morceaux[END_REF] or PLIC [START_REF] Aniszewski | Volume of Fluid (VOF) type advection methods in two-phase flow: A comparative study[END_REF]. Normal vectors and curvature are calculated using Height Functions technique [START_REF] Cummins | Estimating curvature from volume fractions[END_REF][START_REF] Popinet | Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF], which has multiple provisions for the cases of insufficient grid resolutions and/or specific interface configurations.

In the PLIC (Piecewise Linear Interface Calculation [START_REF] Tryggvason | Direct numerical simulations of gas-liquid multiphase flows[END_REF]) scheme, the fluxes Cu are calculated geometrically as intersection volumes between cubic grid-cells and planes, the latter representing the interface. Interface is represented in planar form which is done by first retrieving the vector n normal to the interface (oriented along the ∇C); allowing us to the interface equation:

nx = α, (28)
where the α term (shift along the normal vector) has yet to be found. Since, however, for each cell the fraction function value C i jk is known at the beginning of the timestep, we want the interface representation to be consistent with it, in other words the volume V i jk of fluid contained on one side of the interface has to be consistent with the fraction function value. However,

V i jk = V i jk (α) = C i jk (29)
which is solved allowing α to be found. Equation ( 29) can be solved analytically even in three dimensions [START_REF] Scardovelli | Analytical relations connecting linear interfaces and volume fractions in rectangular grids[END_REF][START_REF] Aniszewski | Large eddy simulation of turbulent two-phase flow[END_REF] however, approximate solutions have been published as well (e.g. Newton-Rhapson used in [START_REF] Ménard | Coupling level set/ volume of fluid/ ghost fluid methods, validation and application to 3d simulation of the primary breakup of a liquid jet[END_REF]). Moreover, various "simplified" versions of geometric Volume of Fluid method exist that either use an "incomplete" interface representation (e.g. without the α term [START_REF] Malan | Numerical simulation of bubble competition during micro-spalling[END_REF]) or don't employ planar representation at all [START_REF] Xiao | Revisit to the THINC scheme: A simple algebraic vof algorithm[END_REF]. In our implementation of PLIC/CIAM, once the full planar interface reconstruction is finished and ( 28) is known, the cube-chopping algorithm [START_REF] Scardovelli | Analytical relations connecting linear interfaces and volume fractions in rectangular grids[END_REF][START_REF] Tryggvason | Direct numerical simulations of gas-liquid multiphase flows[END_REF] is used to find fluxes Cu and enable solution of (27). In our implementation, fluxes are found in a "split" manner, e.g. along each of the axes separately -the order being permuted at each solver iteration [START_REF] Pilliod | Second order accurate volume of fluid algorithms for tracking material interfaces[END_REF]). 4 The same applies for the approximation of Dirac delta δ s .

More details about the interface tracking and advection methods used here can be found e.g. in comparative work [START_REF] Aniszewski | Volume of Fluid (VOF) type advection methods in two-phase flow: A comparative study[END_REF], comprehensive textbook [START_REF] Tryggvason | Direct numerical simulations of gas-liquid multiphase flows[END_REF] or, in context of PARIS-Simulator [START_REF] Ling | Multiscale simulation of atomization with small droplets represented by a Lagrangian point-particle model[END_REF][START_REF] Salvador | Analysis on the effects of turbulent inflow conditions on spray primary atomization in the near-field by direct numerical simulation[END_REF].

Free Surface Solver

To implement the boundary conditions on the interface, a free surface method has been programmed in PARIS-Simulator [START_REF] Malan | Numerical simulation of bubble competition during micro-spalling[END_REF]. The flow within the pores is not explicitly solved for (except up to two interior cell layers as explained below), additionally in this paper we initially set the p p pressure value to zero. With that assumption, the flow is in fact a quasi-single phase flow, but surface tension is accounted for on the pore surfaces. Still, in this case the pressure p p inside the pores is defined, although it is assumed spatially constant inside them and is calculated using the polytropic equation of state [START_REF] Malan | Direct numerical simulations of pore competition in idealized micro-spall using the VOF method[END_REF].

1. The pressure field is extrapolated onto the interface which is necessary for (25). The boundary value of the gradient is computed as

p * = p p -σκ, (30)
where p p is constant pore pressure, and σκ is Laplace pressure, thereby assuring that surface tension effect is accounted for. The p * value is used for ∇ h operators in cells neighboring the interface; 2. The velocity field is extrapolated to within the pores; these values are necessary for higher-order gradient operators for cells neighboring the interface. Extrapolation is based on liquid velocity, which is extended to two grid-cell layers (or 'levels') within pore/bubble interior by geometrical fitting and least-squares minimization; 3. Finally, the extrapolated velocities in level 1 and 2 cells are corrected to ensure the new field is divergence free.

Computational tests have shown that the above procedure amounts to approximately 10% of CPU cost in a massively parallel simulation compared with Poisson solver cost standing at 80%.

A limitation introduced by current implementation of the above algorithm is-apart from aforementioned necessity of pre-seeding the pores-that the voids cannot coalesce, as it is not yet made compatible with the point 2 above (velocity extrapolation to within the pores). of the above list, as topology change of the pores cannot be accounted for. The model however allows for the pore collapse, to which the system responds with a pressure pulse. It is also interesting to notice that the pores can undergo displacement (i.e. move) for example due to non-balanced surface tension force distributions. More details on the Free Surface sub-solver implementation can be found in [START_REF] Malan | Direct numerical simulations of pore competition in idealized micro-spall using the VOF method[END_REF].

Implementation

In implementation, we apply the Hypre package's [START_REF] Falgout | hypre: a library of high performance preconditioners[END_REF]) SMG (semi-coarsening multigrid solver with 3D plane smoothing) for numerical solution of the Poisson problems (25) on the structured, cubic meshes. Advection terms in ( 23) are calculated using central differencing with WENO [START_REF] Shu | Essentially Non-Oscillatory and Weighted Essentially Non-Oscillatory schemes for hyperbolic conservation laws[END_REF] and Superbee-type [START_REF] Sweby | High resolution schemes using flux limiters for hyperbolic conservation laws[END_REF] slope limiters. A cavity tagging and Lagrangian particle-tracking algorithm [START_REF] Ling | Multiscale simulation of atomization with small droplets represented by a Lagrangian point-particle model[END_REF] are used to retrieve statistics of cavity sizes. Interface curvatures for (24) are calculated with a combination of techniques providing for various stencil shapes and built within the framework of Height Function method [START_REF] Cummins | Estimating curvature from volume fractions[END_REF] and with the mixed-Youngs-central [START_REF] Youngs | Numerical simulation of turbulent mixing by Rayleigh-Taylor instability[END_REF] scheme for normal vector n computation. Temporal discretization is performed using the second-order Crank-Nicholson scheme [START_REF] Carnahan | Applied Numerical Methods[END_REF]; the code can also be switched to first order Euler scheme for testing purposes.

The method was validated comparing numerical solutions to solutions of the Rayleigh-Plesset equation for oscillating bubbles [START_REF] Malan | Numerical simulation of bubble competition during micro-spalling[END_REF]Malan et al., 2018, Fig. 5 in the latter). In the test, a perfect match was found between the Rayleigh-Plesset prediction for the cavity volume temporal evolution and a computed volume, proving that the chosen combinations of computational methods neither adds nor looses traced volume. Additional test is provided in section 3.5

Free Surface Advection Test: 2D Expanding Cavity

Here, we briefly present a simulation set up to test the accuracy of the VOF-PLIC advection scheme coupled with the Free Surface solver (as described in Sect. 3.3). Unlike synthetic advection tests performed traditionally in context of VOF using passive velocity fields -such as circle/sphere advections [START_REF] Aniszewski | Volume of Fluid (VOF) type advection methods in two-phase flow: A comparative study[END_REF]) -the "expanding cavity" test case involves solution of the Euler equation ( 5) i.e. u and p are solved for at each iteration. In this sense, it resembles the high-density droplet advection presented by Vaudor et al. [START_REF] Vaudor | A consistent mass and momentum flux computation method for two phase flows. application to atomization process[END_REF] which also solves for all terms of Navier-Stokes equations in a flow with simple initial conditions and a trivial analytic solution . Here, too, any errors in the advection procedure or Free Surface coupling would manifest itself.

The premise of the test is simple. A two-dimensional, square domain of size 1 is defined as filled with 'liquid' of density ρ l = 1, with a single circular pore (cavity) initialized in the center. Isotropic expansion is imposed by setting a constant normal outflow velocity u n such that ∀t : u n (t) = 0.05. This way liquid flux through each of the walls is equal to q = 0.05×1, thus the rate of cavity volume (surface) growth must be described by:

R c (t) = πr 2 0 + 0.05 • 4t, ( 31
)
where R C (t) is the cavity surface at time t, and r 0 stands for initial cavity radius. In the simulation, by convention the grid cells with C i j = 1 signify the cavity interior, while those with C i j = 0 are in the liquid (intermediary values 0 < C i j < 1 represent interface cells). Thus, total sum of fraction function may be used to represent the computed cavity surface V c (t) at any time:

V c (t) = N i, j C i j (t) • (∆x) 2 , ( 32
)
where N is the number of grid points in each direction, and ∆x is grid-cell size. Combining (31) and (32) yields a definition of the L 1 advection error as a function of time:

L 1 (t) = |πr 2 0 + 0.05 • 4t - N i, j C i j (t) • (∆x) 2 |. ( 33
)
For the test, a relatively large value of initial cavity radius was chosen at r 0 = 0.2, to ensure that the pore is properly resolved at 16 2 grid (3 grid-cells in radius). Cavity was expanded over 2.5 • 10 4 time-steps fixed at dt = 1 • 10 -4 resulting in t end = 2.5, and final radius of r end = 0.45. Surface tension was set to zero. We have investigated the error inherent in initial PLIC representation of the circle at each grid resolution -caused by the fact that the interface is represented by linear segments. This reconstruction error (noted E r ) is easily calculated as difference between πr 2 0 and V c (0). It is summarized in Table 1. We can see that second order convergence of the interface reconstruction scheme is obtained, as expected for PLIC type method (see for instance [START_REF] López | An improved plic-vof method for tracking thin fluid structures in incompressible two-phase flows[END_REF], Table 1). Representation errors have been accounted for in the final test result. 2. Proper convergence is observed at order roughly one, which only improves when the finest grid is used. This is consistent with order of the procedure we have chosen in this test to perform extrapolation of u to within the pores (see also Malan et al., 2018, Sect 3.3 therein) and with first order Euler discretization in time used in this test5 . Although not shown here, we have repeated the test in three dimensions finding a similar convergence. We thus conclude the method properly conserves tracked cavity volume, as the latter remains in good agreement with analytic prediction. Having confirmed the volume tracking capability of our coupled Free Surface-VOF code, we now focus our attention on its application in context of micro-spallation.

Micro-spall simulation setup

Initial Conditions

To simulate the idealized micro-spall phenomena specified in Section 2, we have set up the simulations as follows. Cubical computational domain is used (see Fig. 2), containing a given number of pre-seeded pores in a FCC (Face-Centered Cubic) lattice. It is "expanded" in x direction as assumed in Section 2. More specifically, the proper temporal evolution of the velocity field described by ( 1) is ensured by initializing:

• the x-velocity component to u = x with proper translation and scaling, so that for t = t 1 we have u(0, t 1 ) = -u n , u(L/2, t 1 ) = 0 and u(L, t 1 ) = u n where u n is the outflow velocity;

• the outflow boundary condition for the velocity field u = (u, v, w) in such manner, that

u n = 1/t 1 , (34)
where t 1 > 0; for subsequent moments of simulated time t (34) is held so that |u n (t)| = 1/t at all times;

• the v and w components equal to zero.

This implies usage of Dirichlet boundary condition for u in the x direction, which is time-dependent as was said above. Neumann boundary condition for p is thus applied on the x+ and x-walls. Imposing Dirichlet boundary conditions for pressure -instead of for velocity -could be a valid choice for this type of flow, however -unlike the cases in which pores are wholly contained in the domain [START_REF] Malan | Direct numerical simulations of pore competition in idealized micro-spall using the VOF method[END_REF] -uni-axial expansions discussed here involve pores crossing outer domains along x axis. This in turn would make Dirichlet condition for pressure inconsistent. Periodic boundary conditions are imposed on YZ walls. For the condition to be compatible with the pore cluster geometry, it is re-shaped by adjusting the pore-free buffer surrounding it. More precisely, due to periodic condition in YZ directions, the buffer is present only in the x direction (visible as empty region in Fig. 2), i.e. for x < 0.12 and x > 0.88 the domain contains no pores at t = t 1 . Cubic domain of size L = 1 is assumed with 365 pores. Thus, we arrive at ≈ 0.14.

Physical parameters and computational grids

We describe simulation results concerning the uniaxial expansion using example simulations whose parameters are given in Table 3.

Consider now the simulation domain (-L 0 /2, L 0 /2) 3 within which spherical pores are initialized with non-zero radii (which is necessary due to limitations of PARIS-Simulator code) i.e. numerical simulations are started at t 1 > t 0 . With this in mind, the mean radii distribution right after pore opening is

R(t 1 ) = L 0 3 4πN 0 t 1 t 0 -1 1 /3 , ( 35
)
which is how pores are defined in the initial simulation condition. In this computational configuration, it is desirable for R(t 1 ) to be possibly small, in order to offer a large range of scalesin other words, headroom for V(t) growth before the pores coalesce. A favourable relation would be

∆x R(t 1 ) L 0 , ( 36
)
where ∆x is the grid size. The leftmost and rightmost inequalities of ( 36) are however slightly relaxed for the results presented in this paper due to grid resolutions used, as will be detailed below.

Example simulations have been performed using 256 3 grid points. For the first simulation ("I"), value of Weber number (15a) at t 1 is 512.83. The pore lattice is configured as specified above, with liquid characteristics found in Table 3. In the Table, "I" and "C" are labels designating the "Isolation" and "Competition" regimes; We and Ma are dimensionless Weber and Mach numbers, ρ and σ are respectively liquid density and surface tension applied at the pore surfaces, while var(r) stands for the variation in pore radii applied to the initial condition. Mach number values presented in Table 3, are calculated using the speed of sound c as defined by ( 22).

As mentioned above (see ( 35)), due to the applied spatial resolution, certain restrictions on the initial pore radii R(t 1 ) are imposed since the pores must be properly represented by the interface tracking method at t 1 . This representation requires e.g. that gradients of the fraction function C can be resolved near the interface (to calculate interface normals) with finite difference operators. Having the size of domain L 0 and initial pore number N 0 fixed, and imposing additional restriction R > 3∆x we e.g. choose t 1 = 0.08 for the presented simulation of isolation regime (see Table 3). This amounts to R(t 1 )/∆x ≈ 3.44 which is considered resolved [START_REF] Tryggvason | Direct numerical simulations of gas-liquid multiphase flows[END_REF] in that proper values of curvature can be computed for pores of that radius6 . Due to the dependence (8) between simulation initial time t 1 and bubble/pore radii, any decrease of t 0 , and t 1 /t 0 would reduce R, which in turn would imply an increased grid resolution.

Results

Isolation regime

Figure 4a presents cluster geometry at t = t 1 + 1.929 • 10 -2 . Uniform expansion of the bubbles/pores is clearly visible with pore layers closest to the walls along the x axis visibly elongated due to the fact that convection is strongest there. At t ≥ t 1 + 1.92 • 10 -2 one observes pore radii of order . Nearly all pores are ellipsoidal in shape, and have expanded beyond the region of the simulated volume, including periodic yz walls. Moreover, the outer layer of pores has now completely disappeared (have been convected out) from the simulation domain, and only remnant interfacial cells are visible.

Temporal evolution of the pore volume fraction V g for "I" regime flow is seen in Fig. 5. Seeing as all volume change in the liquid results from growth of the pores, we can supplement (4) by a following expression for pore volumes:

V g = L 3 0 (1 - t 0 t ). ( 37
)
The initial pore volume fraction at t 1 , i.e. V g (t 1 = 0.08) should amount to 3.59 • 10 -3 , however it is slightly higher at simulation init, due to the initial variance in pore radii (Table 3). Proper scaling has thus been used for the analytic formula (37) to make up for that in Fig. 5.

Generally, one observes a difference between analytic prediction and simulated pore value in Figure 5. In our opinion, there are two main reasons for that. One is -as mentioned in context of Fig. 4a -that the pores leave the computational domain while it is expanded, which in the "I" simulation takes place at t + t 1 > 0.05 so very near to the beginning of the x range (this moment is marked in Fig. 5 as t ex ). Naturally, this means that vacuum volume (areas with C = 1) associated with these pores is no longer represented in the sum shown in Fig. 5. Moreover, once the pores are sufficiently large the case is no longer dilute, thus departing from the model assumptions as the R(t) requirement of (36) doesn't hold, also the pores are no longer spherical.

For the isolation regime simulation, we additionally present in Fig. 6 the profiles for p (Fig. 6a) and u x (Fig. 6b) measured along the x axis. These are taken from the early simulation stage (t = t 0 + 6.732 • 10 -3 ) and confirm that the evolution of the system is indeed in the ballistic regime, while the velocity profile conforms to linear as predicted in (1). The 'wavy' character of the plot in Fig. 6b, is caused by inclusion of zero velocities (inside the pores) to the average. Similarly, we observe no variance in p (Fig. 6a) when focusing our attention on the lattice-occupied part: for x < 0.12 and x > 0.88 there is an empty "buffer zone" in the x direction, as mentioned above.

Competition regime

We continue the description of the results with the second simulation ("C") which, as seen in Table 3 is characterized by a much lower Weber number. As a result of this, the evolution of the pore lattice is no more dominated by expansion and the pores are not isolated, and the elongated shapes of the pores visible in Fig. 4a give way to a "pore competition" phenomenon [START_REF] Everitt | Competition and interaction of polydisperse bubbles in polymer foams[END_REF]. Figure 4b, shows that number of pores has either shrunk or is at the verge of imploding, their volume being overtaken by a group of large pores. Distribution of radii in the latter group is rather isotropic, with all pores roughly spherical, while ellipsoidal forms are absent. For situation presented in Fig. 4b average pressure p has a negative value. This is due to the fact that positive pressure field is associated with expanding liquid, while for low Weber number negative Laplace pressures (capillary force) dominate the pressure distribution. We will revisit this subject in Fig. 10.

A comparison between (37) and calculated pore volume (analogous to Figure 5) is presented for the "competition" regime in Figure 7. For this low Weber number case, we note a virtually perfect fit of calculated total pore volume (obtained using a sum of fraction function C) up until t ex = 0.133 which is the moment at which first pores exit the calculation domain. The overall better fit might be due to the fact that in the absence of dominating uni-axial expansion, the pores remain generally spherical, as assumed in Sect. 2. It is interesting to note that in this case V g (t) minimally surpasses predicted values just after t ex , before visibly dropping below it for tt 1 > 0.25. Initial slight increase in V g (t) might be due e.g. to interaction with periodic walls, which the pore cluster contacts first. The subsequent decrease in V g (t) is coincident with first pore implosions and should be associated with them.

It is worth noting both Figures 5 and 7 may serve as an additional, indirect evidence of traced mass/volume conservation by our method. This is due to the fact that pore expansion described by ( 37) is the only source of the increase of the total C function, and a good agreement with it is reached at least as long as pores remain spherical and don't leave the domain.

Increased radii variation and implosion events that characterize the competition regime, are visible in Fig. 8a which displays individual pore volume histories for the larger (maroon lines) and smaller (turquoise lines) Weber numbers. To the accuracy allowed by post-processing software, each line displays volume of an individual pore, with up to 100 of the randomly chosen pores tracked for each simulation. Time is normalized by capillary timescale [START_REF] Malan | Direct numerical simulations of pore competition in idealized micro-spall using the VOF method[END_REF])

τ R = ρR 3 σ 1 /2 , (38)
which results (due to parameters presented in Table 3) in factor of 5 between the two simulations presented in Fig. 8a, accounted for in the figure. Unlike the high-We regime, in which expansion is rapid and nearly uniform, the low-We competition regime exhibits higher variance in pore volumes, with a number of implosion events visible towards the end of recorded time. This proves that the presented numerical method is capable of capturing the transition between the two regimes, provided that simulation covers a sufficient temporal range7 . Additionally in Fig. 8a we observe a sub-population of bubbles exiting the domain through the side-walls towards τ D = 0.05. This flow stage is visualized again in Figure 8b which contains a magnification for τ D ∈ [0.02, 0.05]. Even at this time scale we observe that a distinct group of pores loses volume at identical rate: those lines represent the group of approximately 50 pores that are carried out of the simulation domain. It is a visibly different evolution than that presented by the maroon lines in Fig. 8a. The latter is characterized by individual collapses/implosions taking place for τ D > 5.5, and absence of distinct pore sub-groups. The same character of volume evolution in the competition regime has previously been reported in [START_REF] Malan | Direct numerical simulations of pore competition in idealized micro-spall using the VOF method[END_REF] (see Figure 7 therein, note their plots show all individual pores while Figure 8 presents only popula- tion subset). Most importantly, their result -even if obtained for isometric expansions -shows a trend of pore population to create discrete groups characterized by specific growth rates, which is more visible the higher We is. Our result visible in Fig. 8a and 8b presents an extreme case of that trend: for "ballistic" expansion only one group exists at We = 512. Both pore populations visible in Figure 8a reach the coalescence region: for the We = 512 population (maroon) the volumes recorded reach 2 • 10 -3 , while for competition regime (turquoise) they are slightly higher at 3.5 • 10 -3 . This translates to final equivalent sphere radii of r ≈ 8 • 10 -2 and r ≈ 9.42 • 10 -2 for isolation and competition regimes respectively. Since these radii are in both cases larger than /2, coalescence will appear. The probability of pore contact occurring is of course higher in the isolation regime with all pores present. Ultimately, coalescence is the reason why simulation of either regime cannot be continued beyond the τ D values presented in Figure 8a due to implementation limitations mentioned in Section 3.3.

Evolution of the pressure field

We continue our comparison of the presented uniaxial expansion regimes with Figs 9a and 9b, first of which displays the evolution of averaged pressure field

p = 1 V V pdV (39)
for both analyzed flow regimes. As the high-We regime is characterized by a more violent decrease in the average pressure in the initial parts of expansion, it was necessary to plot p in logarithmic scale in Fig. 9a. However due to the function being negative, we instead plot in Fig. 9a the

log(| p |) • sgn( p ). ( 40
)
It is clearly visible that the rapid decrease from e 10 ≈ 2.2 • 10 5 to zero is followed by a period for which p < 0. Figure 9b presents a plot comparable to Fig. 9a, but obtained by discarding Laplace pressures (i.e. capillary contributions). Here, we are using only input from grid cells that do not coincide with pore surfaces (this can be found from values of C function)-in other words, grid-cells positioned away from the pores and the interface. This can be roughly approximated by

p = 1 V V p(1 -C)dV (41)
where C stands for the fraction function (C is associated with pore volume, meaning that a grid cell with C = 1 is inside a pore, C = 0 corresponds to the liquid, while 0 < C < 1 signifies an interface cell). Accuracy provided by ( 41) is moderate, as it will still include contributions from direct neighbours of the interface grid-cells; those neighbouring cells are generally still influenced by Laplace pressures. While better accuracy could be provided by a distance function-based methods [START_REF] Aniszewski | Volume of Fluid (VOF) type advection methods in two-phase flow: A comparative study[END_REF], the above formula is enough to yield a much different result than depicted in Fig. 9a. Indeed, while a similar pressure evolution as in Fig. 9a is observed in Fig. 9b, we notice

∀t : p (t) > 0 ∧ ∂ p ∂t < 0 . (42)
This is to say average pressure is always positive and decreasing. Additionally, slight oscillations visible in Fig. 9a for 0.025 < t < 0.035, are not found in Fig. 9b, thereby reassuring us that they were attributed to the capillary forces. A different p evolution is observed for the "competition" regime. In Fig. 10 for low We, which (as we mentioned in context of Fig. 4) is dominated by capillary pressures nearly from the onset, we notice p ∈ [-50, -20] within the computed time interval. For t > the individual peaks are visible (as triangles above other curves in Fig. 10) which should be associated with implosion (collapse) events that result in shocks to the liquid.

Pore implosions, encountered only in competition regime (low We), are difficult to treat numerically in codes that use uniform grid size cells. Once a pore decreases in size heading into collapse, so does the number of grid-cells coinciding with its surface (i.e. the interface). It eventually becomes so small that interfacial curvature term σκnδ s of (24) can not be reliably calculated. In our implementation [START_REF] Malan | Direct numerical simulations of pore competition in idealized micro-spall using the VOF method[END_REF] when encountering this situation PARIS-Simulator removes the vapor cells from the solution, replacing them with liquid cells, supplying a volume sink in form of a source term added to the r-h-s of (25). A consequence of that is a shock wave propagating to the system as described in Section 2.4. This explains pressure peaks visible in Figure 10.

It is easy to note that unlike in the isolation regime, in the competition regime we find for t > 0.02

∂ p ∂t > 0, (43)
i.e. average pressure is steadily growing, which in fact is predicted by ( 19). We would associate this with the initially rapid expansion slowing down, and competition process taking over. The first expansion phase (or the isolation regime) corresponded to rapid expansion when the pores are still small compared to , while the latter is due to Laplace (capillary) pressures which-in spite of being inversely proportional to pore radii-become dominant, and competition phenomenon appears [START_REF] Malan | Direct numerical simulations of pore competition in idealized micro-spall using the VOF method[END_REF]. Note that measured in absolute time values, evolution presented in Fig. 10 involves a temporal interval nearly ten times longer than that of Fig. 9.

To further validate the predictions of the numerical simulation result presented in Figure 10 we have included a curve displaying temporal evolution of p ∞ (t) calculated using (19). The sharp pressure drop resulting from initial shock, as well as the beginning of pressure growth after minimum is reached (mentioned in Section 2), are well visible in Fig. 10 for both curves, and there is overall agreement in their shape. Obviously, mean pressure p is not equal to p ∞ and is strongly influenced by the capillary pressures originating from surface tension on pore surfaces. For this reason, minimum p value is predicted too low. Interestingly, if instead of using σ = 1 in (19) we calculated p ∞ using σ = 1.3, we get min t (p ∞ (t)) = min t ( p (t)) and a better agreement of the p ∞ curve to p overall which strongly suggests that indeed the contributions of Laplace pressures substitute the difference between the analytic prediction and the simulation result. Figure 11 presents an example of an additional simulation 8 with We(t) = 9.56 • 10 -2 at t = t 0 + 0.157. This third presented simulation corresponds to a transition from isolation to competition regimes predicted by (15a). We can illustrate this by directing our attention to values of pressure inside the pore cluster. In Fig. 11, three isosurfaces are presented. First, a gray (semi-transparent) isosurface of fraction function C is visible. This marks the position of individual pores, although at first glance it is hard to conclude if the competition phenomenon is present. Similarly, regions with p ≈ -42 (green, semitransparent isosurface in Fig. 11) exhibit little spatial variability: all pores are surrounded by them. However, once an isosurface for p ≈ -63 is drawn (solid blue surface in Fig. 11) we notice that homogeneity is gone close to certain pores' surface. Upon closer inspection it is evident that these pores have smaller radii (hence larger mean curvatures); they are being acted upon by other, growing pores to eventually succumb to the competition [START_REF] Malan | Direct numerical simulations of pore competition in idealized micro-spall using the VOF method[END_REF]. 8 The simulation uses the same boundary conditions as those mentioned above, resolution is 512 3 grid points.

Conclusions

We have presented a numerical simulation setup that permits prediction of the behavior of a pore lattice within rapidly expanding medium. The expansion corresponds to a uniaxial (one-dimensional) ballistic phenomenon characterized by large initial outward velocity magnitudes which then immediately diminish. On one hand, this setup facilitates recognition of momentum-dominated phenomena such as rapid pore elongation in "isolation" regime. On the other, the decline in expansion rate promotes the onset of pore "competition" and more anisotropic radii distribution in surviving pores. Moreover, simulations allow us to observe transition between these two regimes. We have presented that instantaneous velocity and pressure profiles match analytic predictions, and so does the temporal pressure evolution in the simulated pore systems.

It is relatively easy to point out numerous simplifications of the presented model. One example is the uni-axial mean expansion assumption, another: lack of solution inside pore interiors, yet another is using an incompressible medium. Still, the work presented herein provides us with a rare insight into the micro-spall process at the FCC lattice scale which is otherwise inaccessible experimentally, forcing the researchers to apply modelling and/or indirect analysis. In this context, the present study-especially with its choice of (partly periodic) boundary conditions-may easily be used as a departure point for larger-scale modeling e.g. in the context of metals undergoing spallation under load.

Results presented here can be employed to improve predictions of micro-spallation type damage in an indirect manner. Published computations of micro-spall utilize larger domains i.e. simulate larger material samples -up to 180nm nickel sample in [START_REF] Qiu | Non-equilibrium molecular dynamics simulations of the spallation in ni: Effect of vacancies[END_REF] -over which stress is imposed in similar one-dimensional fashion i.e. [1 0 0] direction of the FCC lattice (consistent with the X axis direction in this work). Individual Ni atoms are presented as Lagrangian particles using the NEMD (NonEquilibrium Molecular Dynamics) approach [START_REF] Wang | An atomic study on the shock-induced plasticity and phase transition for iron-based single crystals[END_REF], and cavities are created by randomly removing atoms. Meanwhile our simulations include domains orders of magnitude smaller -only a few dozen basic 12-atom FCC clusters versus e.g. 5 • 10 6 atoms in NEMD [START_REF] Wang | An atomic study on the shock-induced plasticity and phase transition for iron-based single crystals[END_REF] -but account for some phenomena forcibly neglected at larger scales, such as non-spherical pore shapes and competition.

Thus, assumptions we have taken allow re-using our flow domain as a control volume in a finite volume approach for modelling bigger systems. Slices (of thickness L 0 or L 0 /2) can be created by repeating the domain in periodic YZ directions. In X, either by using entire domain (x ∈ [-L 0 /2, L 0 /2], u ∈ [-u n , u n ]) or only half (x ∈ [0, L 0 /2], u ∈ [0, u n ])) one could construct neighbouring slices by using u n at the right-hand-side of one volume as a condition on the left-hand-side of the neighbouring volume -in this scenario expansion would be uni-directional along X with increasing u n . Pressure could be treated as locally constant in each control volume (see Fig. 6a). Having thus constructed and analysed the larger system, predictions could be either compared with NEMD-type simulations or used as an initial condition for the latter concerning, for example, initial vacancy concentration calculated from the V g (t) evolution presented in this work. This concentration is a critical parameter in spall strength prediction (see e.g. Table 1 in Qiu et al., 2017, and references therein), and, as evident from our results, dependent strongly on the expansion regime.

Figure 1 :

 1 Figure1: Sketched evolution of mean macroscopic volume and pressure of a material element undergoing mico-spall. As qualitatively represented, the initial pressure produced by the shock and the final spall stress are here assumed respectively large and small enough that the eventual expansion is nearly ballistic. This ballistic evolution can then be backward extrapolated to vanishing volume which defines the (arbitrary) time origin. As indicated, the system evolution is here simulated after pores have opened and local fluid density has returned to a nearly uniform and constant value.

Figure 2 :

 2 Figure 2: Schematic representation of idealized micro-spall as a pore cluster undergoing expansion along the x direction.

Figure 4 :

 4 Figure 4: (a) Simulation of the flow characterized by We = 512. (isolation regime) at t = t 1 +1.92•10 -2 . (b) Simulation of the pore lattice at We = 3.3•10 -2("competition" regime) for t = t 1 + 0.237.

Figure 5 :

 5 Figure 5: Volume fraction temporal evolution (isolation regime). Continuous line: simulation; dashed line: using (37). The moment at which first pores touch domain boundary is marked with vertical line and the t ex label.

Figure 6 :

 6 Figure 6: Isolation regime: Averaged (in YZ) profiles of p (a) and u (b) for tt 0 ≈ 6 • 10 -3 . (The x range applies to both sub-plots.)

Figure 7 :

 7 Figure 7: Volume fraction temporal evolution (competition regime). Continuous line: simulation; dashed line: using (37). The moment at which first pores touch domain boundary is marked with vertical line and the t ex label.

Figure 8 :

 8 Figure 8: (a): Pores volume history for We = 512 (maroon), and We = 0.03 (turquoise) flows. (b) Enlargement for τ R ∈ [0.02, 0.05], We = 512.

Figure 9 :

 9 Figure 9: Isolation regime: (a) Evolution of the domain-averaged pressure for the flow at We = 512; (b) The same quantity prepared by including only grid cells in the liquid.

  that

Figure 10 :

 10 Figure 10: Evolution of the mean pressure field for the competition regime. Numerical simulation (black triangles) and analytical p ∞ (t) using (19) (red line).

Figure 11 :

 11 Figure 11: The bubble cluster at We ≈ 0.1, corresponding to a transition from isolation to competition regime.

Table 1 :

 1 Reconstruction error of the circular pore E r at varying grid levels.

	Grid points 16 2 32 2 64 2 128 2	E r 2.97 • 10 -5 8.29 • 10 -6 3.59 Ratio Order --1.79 1.99 • 10 -6 4.16 2.08 4.93 • 10 -7 4.03 2.01

Table 2 :

 2 Final L 1 error values for expanding cavity test at varying grid levels.

	Grid points L 1 (t end )(×10 -5 ) Ratio Order 16 2 7.59 --32 2 4.721 1.608 0.8 64 2 2.801 1.685 0.84 128 2 0.669 4.18 2.09
	Final L 1 (t

end ) values are presented in Table

Table 3 :

 3 Parameters for two examples of simulations of the uniaxial expansion, as explained in text.

	Label	ρ	σ	t 1	t 0	var (r)	We(t 1 )	Ma(t 1 )
	I C	1000 0.1 8 • 10 -2 7.797 • 10 -2 444 1 2.1 2.092	50% 50%	512.83 3.3 • 10 -2	15 0.12

Not used in this work.

https://www.gnu.org/licenses/gpl.html

Available at http://www.lmm.jussieu.fr/∼zaleski/paris/index.html

Other results presented in the paper employ the Crank-Nicholson scheme.

Not accounting for about 5 percent variance in R in initialization.

Which in turns depends directly on numerical grid resolution. The smaller R(t 1 ) the more timespan will be included in the simulation before the pores coalesce.

Acknowledgements

Authors acknowledge partial funding from the CEA (Commissariat à l'énergie atomique et aux énergies alternatives); the computations have been partly performed using the CEA's TGCC facility. Visualizations have been performed using Gnuplot (Williams et al., 2010) and Visit (Childs et al., 2012). The authors want to thank D. Fuster, S. Mercier and the anonymous reviewers for helpful discussions, corrections and suggestions that have helped us improve the paper. W. Aniszewski thanks A.A.S. for proofreading the manuscript.


D Levesque

J P Hansen

The origin of computational Statistical Mechanics in France

The two main methodologies of computational Statistical Mechanics, namely the stochastic Monte Carlo and the deterministic Molecular Dynamics methods, were developed in the USA in the mid nineteen fifties. In the present paper we show how these "computer experiments" migrated to Europe in the sixties, and first bloomed at the Orsay Science Faculty, before spreading throughout Europe. Collaborations between the Orsay group, led by Loup Verlet, and pioneering groups in the USA and Europe are pointed out. Finally it is shown how the celebrated Verlet algorithm for the integration of classical equations of motion can be traced back to Isaac Newton.

Introduction

The basic principles of Statistical Mechanics were established during the second half of the nineteenth and the beginning of the 20th century by Maxwell, Boltzmann and Gibbs. Exact results for the thermodynamic properties of many-particle systems could be derived analytically for some simple models, culminating in the solution of the twodimensional Ising model by Onsager [START_REF] Onsager | Crystal Statistics. I. A Two-Dimensional Model with an Order-Disorder Transition[END_REF]]. However very few models have analytic solutions, so that approximations, like mean-field, must be made to predict thermodynamic, structural or transport properties of physically relevant models.

This situation changed radically with the appearance, during and after World-War II, of the first electronic computers which were mostly based, for defence reasons, at the US National Laboratories at Los Alamos and at Livermore. This led to the development of numerical simulation methods for classical many-particle systems, namely the Metropolis Monte-Carlo (MC) method [START_REF] Metropolis | Equation of State Calculations by Fast Computing Machines[END_REF][START_REF] Rosenbluth | Monte Carlo Calculation of the Average Extension of Molecular Chains[END_REF], Wood 1957a] and the Molecular Dynamics (MD) method [START_REF] Alder | Phase Transition for a Hard Sphere System[END_REF][START_REF] Alder | Studies in Molecular Dynamics. I. General Method[END_REF]]. The initial MD code of Berni Alder and Tom Wainwright, assisted by Mary Ann Mansigh, was designed for hard disks (in 2d) or hard spheres undergoing elastic binary collisions. The extension to systems interacting via "soft", continuous pair potentials was later put forward by Anees Rahman at the Argonne National Laboratory [START_REF] Rahman | Correlations in the Motion of Atoms in Liquid Argon[END_REF]]. The early MC and MD simulations led immediately to a major, unexpected discovery, namely the purely entropy-driven freezing of hard disks and spheres at sufficiently high packing fractions [START_REF] Alder | Phase Transition for a Hard Sphere System[END_REF], Wood 1957b].

The pioneering period of computational Statistical Mechanics is well described in a recent paper by Michel Mareschal [START_REF] Mareschal | Early Years of Computational Statistical Mechanics[END_REF]]. In the present paper we report how computational Statistical Mechanics was imported to Europe, and more specifically to France, where the first group developing "computer experiments" worked at the Orsay Faculty of Science. In Section 2 we follow the emergence of a computational Statistical Mechanics group within the newly founded Laboratory of Theoretical Physics and High Energies at Orsay in the late 1950's. Section 3 describes the context of the development of MD and MC at Orsay throughout the sixties and beyond. Prolongations and the diversification of these computational activities are analysed in Section 4, while in the Conclusion we report the link between the most commonly used MD algorithm for the integration of the equations of motion and a Theorem due to Isaac Newton!

The Orsay Statistical Physics group

In 1960, the Orsay Science Campus (located 20 km south of Paris), still part of the Faculty of Science of the University of Paris (Sorbonne), is in a phase of rapid development. Several buildings are under construction, to host new laboratories, in particular the "Laboratoire de Physique Théorique et Hautes Énergies" (LPTHE). From the onset, its founder and director, Professor Maurice Lévy, wished the LPTHE to be multi-disciplinary, bringing together theoreticians of elementary particles physics, nuclear physics, statistical physics and mathematical physics. In these early days, LPTHE numbered about 20-30 scientists, including PhD students, and occupied half a floor of the freshly built "211" building, which also hosted solid-state and plasma physics laboratories.

The statistical physics group was led by Loup Verlet, a "transfuge" from elementary particle and nuclear physics [START_REF] Drell | Multiple Scattering Corrections in π ± Deuteron Scattering[END_REF]]. In 1959, he derived the hyper-nettedchain (HNC) equation [START_REF] Verlet | On the theory of classical fluids[END_REF]], simultaneously and independently of other authors, using a diagrammatic re-summation method, a highly "fashionable" exercise at that time. The HNC equation, as well as the Percus-Yevick (PY) equation derived around the same time [START_REF] Percus | Analysis of Classical Statistical Mechanics by Means of Collective Coordinates[END_REF]] are approximate closure relations which supplement the Ornstein-Zernike relation between total and direct pair correlation functions [START_REF] Ornstein | Accidental deviations of density and opalescence at the critical point of a single substance[END_REF]]. These integral equations were at the basis of theoretical developments in statistical physics of dense fluids and liquids throughout the sixties and beyond.

While the PY equation admits an analytic solution for the case of a hard sphere fluid [START_REF] Thiele | Equation of State for Hard Spheres[END_REF][START_REF] Wertheim | Exact Solution of the Percus-Yevick Integral Equation for Hard Spheres[END_REF][START_REF] Baxter | Ornstein-Zernike Relation for a Disordered Fluid[END_REF]], no such solutions exist for the HNC equation, nor for more "realistic" models, like the Lennard-Jones (LJ) pair potential, of real fluids, within either PY or HNC. It is hence clear that one had to resort to numerical solutions of these integral equations in order to allow a quantitative comparison of the predictions of these approximate theories with experimental data.

In that perspective, the statistical physics group at LPTHE was to benefit from the creation by Maurice Lévy and Loup Verlet, with the support of other laboratories, of the computing centre of the Orsay Faculty. The computer which was installed in building 211, was one of the first commercially available machines in the world, namely the IBM 650. The central processor unit still used vacuum tubes, like radio sets of that period! The central memory of 2000 words of ten decimal digits was located on a magnetic drum which rotated at the frequency of several thousand revolutions per minute and contained both instructions and data. The execution time of instructions and arithmetic operations was of the order of a millisecond. Programming was preferentially in machine language, each elementary operation being coded on punched cards. Shortly thereafter a second machine was installed at Orsay Faculty, namely a CAB 500 computer built by SEA (Société d'Électronique et d'Automatisme) company. This machine had a magnetic drum memory of 16000 words of 10 digits; the typical execution time of arithmetic operations was of the order of 40 milliseconds.

Loup Verlet and his PhD student Dominique Levesque used mainly the IBM 650 to get acquainted with calculations on a computer. They developed programs to solve the HNC and PY equations by an iterative method carried out in real space since the Fourier transformation did not provide any advantage within the IBM 650 architecture and its Fast Fourier version [START_REF] Cooley | An algorithm for the machine calculation of complex Fourier series[END_REF] clearly not adapted to decimal arithmetics.

Meanwhile the seminal papers of the MC and MD methods by Metropolis et al., Rosenbluth and Rosenbluth, Wood and Parker and Alder and Wainwright already quoted, had not escaped the attention of the statistical physics group of LPTHE and were studied in detail by the Verlet group. The results of these papers, in particular those of Alder and Wainwright for hard spheres, allowed an appreciation of the quantitative validity of the PY and HNC equations for that system. A similar comparison with experimental data for rare gases also provided such an appreciation, which is however less reliable because of the uncertainty of the interatomic interaction modeled by the LJ pair potential. In view of these comparisons, improvements of the PY and HNC equations were put forward by the Verlet group. Functional expansions were used to derive corrections to the PY equation, while an attempt was made to improve the HNC closure by including some so called "bridge" diagrams neglected within HNC [START_REF] Verlet | On the theory of classical fluids-III[END_REF][START_REF] Verlet | On the theory of classical fluids-IV[END_REF]. It is for computing the simplest bridge diagram that the Verlet group used for the first time the MC method [START_REF] Verlet | On the theory of classical fluids II[END_REF]].

In 1965-1966, Loup Verlet was invited by Joël Lebowitz to spend one year at the Belfer Graduate School in New-York. This was the year after the publication of the seminal paper by Anees Rahman on the MD simulation of a system of particles interacting via the LJ pair potential [START_REF] Rahman | Correlations in the Motion of Atoms in Liquid Argon[END_REF]]. Loup Verlet who had noticed this article before his departure for the US, decided to extend Rahman's pioneering work along several lines. First, Lebowitz, Percus and Verlet established that, since a time average taken within the course of a MD simulation, at constant total energy and momentum, is equivalent to a micro-canonical average, all thermodynamics quantities can be computed from the time evolution of the dynamical variables themselves or their fluctuations [START_REF] Lebowitz | Ensemble Dependence of Fluctuations with Application to Machine Computations[END_REF]]. For instance, the time average of the kinetic energy fluctuations allows the calculation of the system's specific heat. Next, on an up-to-date CDC computer newly accessible at the Courant Institute of New-York University, Loup Verlet implemented a MD code innovating along two lines, namely an algorithm for the integration of the Newtonian equations of motion and the awareness of the range of the molecular and atomic interactions [START_REF] Verlet | Computer "Experiments" on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules[END_REF]]. These fundamental innovations are now universally acknowledged and used.

Concerning the second point, the work on the HNC and PY equation had shown that the asymptotic decay of the LJ potential could be accounted for by assuming that at large distances the pair distribution function is equal to its asymptotic value g(r → ∞) = 1. Verlet exploited this observation in his MD code by using appropriate near neighbour tables around each particle, optimally up-dated to account for the relative motion of the N particles. This procedure is easily validated by comparing the results obtained with these tables to simulations taking into account all N (N -1)/2 pairs of particles. The tables allow the simulation of much larger systems, since the computation time becomes of the order N ln N rather N 2 . In Verlet's initial implementation [START_REF] Verlet | Computer "Experiments" on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules[END_REF]] for a system of 864 particles, the simulation time, for comparable statistical errors, was reduced by more than a factor of 4.

Concerning the integration algorithm, Verlet replaced the somewhat cumbersome Runge-Kutta algorithm used by Anees Rahman by the appealingly simple algorithm:

- → r (t + dt) = 2 - → r (t) -- → r (t -dt) + dt 2 - → f (t)/m (1)
where -→ r (t) is the position of a particle with mass m at time t and -→ f (t) is the total force acting on the particle at time t. The algorithm is obviously reversible and valid to order dt 4 .

The original MD program was written by Verlet in Fortran II CDC. On his return from New-York, it was immediately transcribed on the new computer installed in Orsay: the Univac 1107, followed by the Univac 1108. The latter was, during the period 1967-1970, one of the most powerful commercially available computers. The execution times of arithmetic operations were of the order of few microseconds, while the operating system was among the first multi-tasking. Engineers from IBM-France came to Orsay to observe de visu this functionality which is now standard!

The early MD runs sparked a number of questions concerning the size of the simulated systems and the rigour with which the conservation laws had to be satisfied. The non-conservation of the momentum was the simplest case, since it results exclusively from the choice of initial conditions and the round-off errors affecting the calculation of the forces. Regarding the latter, the Univac 1108 had the advantage of coding floating point numbers with 36 bits (rather 32 on IBM). Moreover the Verlet algorithm Eq. ( 1) is independent of particle velocities, the latter being only required for the calculation of the total energy and the temperature via:

- → v (t) = - → r (t + dt) -- → r (t -dt) 2 dt (2)
The integration time step was chosen empirically to ensure energy conservation within a relative error 10 -3 -10 -4 over runs of several thousand time-steps after equilibration. Simulations were generally run at night or during week-ends by members of the Verlet team so as to benefit from favourable financial conditions, and to avoid delaying the calculations of other users. While the team had full confidence in the validity of the calculated thermodynamic and structural quantities due to the equivalence between time average, computed by MD simulation at constant energy, and micro-canonical average, questions arose on the precision to be achieved on the particle trajectories for a realistic estimate of the time-dependent correlation functions, like the velocity autocorrelation functions. In the late sixties, the importance of accounting for the conservation of the phase space volume had not yet been clearly recognized. In particular the fact that Anees Rahman's algorithm is not symplectic, while the Verlet algorithm Eq. ( 1) is, had not yet been noted.

It was hence on the basis of the equivalence, anticipated by Boltzmann, between MD and the micro-canonical ensemble, combined with a heuristic implementation of practical details that the use of MD flourished at the LPTHE in Orsay. It was, of course, recognized both at Orsay and elsewhere, that a MD program is readily transformed into a MC program at the cost of a few minor modification and the inclusion of a subroutine providing a truly reliable random number generator.

In particular, independently of the seminal Quantum Monte Carlo (QMC) paper by W.E. [START_REF] Mcmillan | Ground State of Liquid He 4[END_REF]], Verlet and his students used the formal analogy between a Jastrow trial wavefunction [START_REF] Jastrow | Many-Body Problem with Strong Forces[END_REF]] and the Boltzmann factor for a system of N classical pair-wise interacting particles, to calculate the ground-state energy of liquid He 4 [START_REF] Levesque | On the Ground State of Liquid and Solid He 4 at Zero Temperature[END_REF]]; the preprint remained unpublished because McMillan's paper appeared in print just as the Orsay preprint was completed. The preprint was later supplemented by a paper on the ground state of He 4 and He 3 , using a Jastrow-Slater wave-function for the fermion case, written by D. Schiff and L. Verlet while both were in New-York [START_REF] Schiff | Ground State of Liquid Helium-4 and Helium-3[END_REF]]. At about the same time, back in Orsay, J.P. Hansen and D. Levesque used a Jastrow wave function, multiplied by a product of Gaussians centered on the crystal lattice positions, within a variational QMC calculation, to determine the ground-state energies of solid He 4 and He 3 [START_REF] Hansen | Ground State of Solid Helium-4 and -3[END_REF]]. These investigations of the ground states of liquid and solid Helium, based on variational approaches, were extended in 1973 -74 within the framework of a collaboration with Mal Kalos. Kalos had initiated and promoted a ground-breaking Monte-Carlo method allowing the exact sampling of the ground-state wavefunction of many-boson systems, referred to as Green's Function Monte-Carlo. The collaboration led to the first application of this exact algorithm to an N-body system, namely a quantum system of hard spheres [START_REF] Kalos | Helium at Zero Temperature with Hard Sphere and Other Forces[END_REF]].

3 The context of the development of MD and MC at Orsay: the late sixties and beyond

The development of MD and MC simulations at LPTHE in the sixties and early seventies was based on the close collaboration of Loup Verlet with his students: Dominique Levesque, Daniel Schiff, Jacques Vieillard-Baron, Jean-Pierre Hansen, and Jean-Jacques Weis and the mutual collaboration of group members. Verlet's group also benefited from a very active international collaboration with distinguished visitors including: Berni Alder, Anees Rahman, John Valleau, Ian Mc Donald, Konrad Singer, Mal Kalos, George Stell, Mark Nelkin and post-docs Juhani Kurkijarvi, Roy Pollock among others.

In the academic research environnement of France during the sixties and early seventies, the simulation work of the LPTHE group was favourably regarded as witnessed by the Langevin prize of the French Physical Society awarded to Loup Verlet in 1971. At this time, it was clear that MD and MC simulations applied to simple classical or quantum systems were able to establish the limits of validity of approximate theoretical approaches, and even to predict unexpected physical phenomena like the freezing of hard spheres or the algebraic (rather than exponential) decay of time dependent correlation functions [START_REF] Alder | Velocity Autocorrelation for Hard Spheres[END_REF], Levesque 1974]. But it was much less clear that these simulation methods would eventually be able to study genuinely complex systems, just as it was difficult to foresee the near exponential growth of computer power, at exponentially decreasing costs, which we have been witnessing in the last decades.

Meanwhile, strengthening the development of the Orsay Faculty, the Centre National de la Recherche Scientifique (CNRS) established on the Orsay campus, at the end of the sixties, the Centre Inter-Régional de Calcul Électronique (CIRCÉ), a national computational facility equipped with an IBM computer. However at the LPTHE, the coding of the MD and MC programs which was partly in Univac machine language, in order to reduce computational times, prevented an immediate portability of the codes to other computers. Moreover, since the installation of the Univac computers at Orsay had hardly been welcomed by IBM, the LPTHE group had developed a kind of "allergy" towards IBM, and used the CIRCÉ facilities only very marginally. Within the CIRCÉ building, Carl Moser, a quantum chemist, created the Centre Européen de Calcul Atomique et Moléculaire (CECAM) in october 1969. Members of the Verlet group did not participate, due to initial misunderstandings and the aforementioned "allergy" to IBM, in the workshops organized by CECAM which contributed strongly to the development of MD and MC methods and their applications throughout Europe. But this initial "boycott" did not affect scientific collaborations between LPTHE group and participants of the workshops who regularly visited the LPTHE and spent sabbatical leaves with the Verlet group.

While the Orsay group focused mostly on improving and exploiting Molecular Dynamics methodologies, the Monte Carlo method in Statistical Mechanics spread simultaneously to other European countries. At Royal Holloway College near London, Ian McDonald and Konrad Singer introduced the very efficient "histogram reweighting" method [START_REF] Mcdonald | Calculation of thermodynamic properties of liquid argon from Lennard-Jones parameters by a Monte Carlo method[END_REF]], and applied it to classical fluids. In Austria, at the Technical University of Vienna, Kurt Binder applied MC sampling to discrete spin systems, in an effort to model magnetic materials and analyse neutron scattering data [START_REF] Binder | Calculation of spin-correlation functions in a ferromagnet with a Monte Carlo method[END_REF]]. At the Munich Technical University, Dietrich Stauffer was the first to use MC methods to study particle percolation on a lattice, which later led to a well-known monograph [START_REF] Stauffer | Introduction to Percolation Theory[END_REF]]. More generally throughout the seventies, participants in CECAM workshops started new computational statistical mechanics groups throughout Europe, including Giovanni Ciccotti in Rome and Herman Berendsen in Groningen among others.

Some prolongations

In the late sixties, and early seventies, the MD and MC simulation projects at LPTHE diversified along several lines, of which we give a few examples. Daniel Schiff used "computer experiments" to investigate simple liquid metals based on effective ion-ion potentials [START_REF] Schiff | Computer "Experiments" on Liquid Metals[END_REF]]. Jacques Vieillard-Baron was the first to simulate assemblies of hard ellipses as a simple model for two dimensional nematic ordering [START_REF] Vieillard-Baron | Phase Transitions of the Classical Hard-Ellipse System[END_REF]; his pioneering work inspired the seminal papers by Daan Frenkel and collaborators on the phase diagram of lyotropic liquid crystals [START_REF] Frenkel | The Hard Ellipsoid of Revolution Fluid: I. Monte-Carlo Simulations[END_REF]]. Jean-Jacques Weis joined forces with Berni Alder and Berkeley experimentalist Herbert Strauss to investigate depolarized light scattering via MD simulations [START_REF] Alder | Depolarization of Light in Atomic Fluids[END_REF]]. Dominique Levesque and collaborators were the first to extend MD to simple models of rigid diatomic molecules [START_REF] Barojas | Simulation of Diatomic Homonuclear Liquids[END_REF]]. Jean-Pierre Hansen joined forces with Ian Mc Donald and Roy Pollock to investigate the static and dynamical properties of the "one component plasma" (OCP) a model for highly compressed plasmas generated in inertial confinement fusion experiments, both by MC [START_REF] Hansen | Statistical Mechanics of Dense Ionized Matter. I. Equilibrium Properties of Classical One-component Plasma[END_REF]] and MD [START_REF] Hansen | Statistical Mechanics of Dense Ionized Matter. III. Dynamical Properties of Classical One-component Plasma[END_REF]] simulations. The Verlet group also engaged in a fruitful collaboration with LPTHE theoretical physicists, in particular with Bernard Jancovici, a former classmate of Loup Verlet at École Normale Supérieure [START_REF] Hansen | Phase Diagram of a Charged Bose Gas[END_REF]].

Around 1974, Loup Verlet and Daniel Schiff moved to other areas of research and human endeavour. In 1974 they published an opinion column in the November issue of the journal La Recherche [Verlet 1974], putting forward their questioning of the objectives of fundamental research in theoretical physics and their motivations to contribute to other areas. Their concerns, in particular those related to ecological issues are still of crucial relevance today. Loup Verlet's practice of theoretical physics also led him to an epistemological analysis of the creation and development of scientific theories published in his two books "La Malle de Newton" [START_REF] Verlet | La Malle de Newton[END_REF]] and "Chimères et Paradoxes" [START_REF] Verlet | Chimères et Paradoxes[END_REF]].

Conclusion

To conclude this brief historical review of the development of the MD and MC simulations at Orsay Faculty, we quote an extract from Loup Verlet's book "Chimères et Paradoxes", on pages 174-175, where Loup Verlet gives his personal account of this period: « Immediately after World War II, the increase in power of computers allowed to reconsider the basis of the question raised by the atomistic hypothesis and addressed by Maxwell and Boltzmann in the 19th century: how to determine the macroscopic properties of gases and liquids, knowing that they are composed by large numbers of atoms obeying Newton's equations of motion and assuming that they interact via a given force law? This question became all the more pressing towards the end of the century, in light of growing evidence that atoms and molecules really exist. Until the end of the World War II, this question remained essentially unsolved except in the case of low density gases. The situation changed in the fifties due to the advancement of powerful computers. It became possible to follow, over a sufficiently long time span, the trajectories of several tens or even hundreds of particles behaving like billiard balls. A few years later it became feasible to tackle a more realistic model, namely about one thousand particles interacting via pair-wise Lennard Jones forces. Within a collective research effort, I suggested, among other procedures to speed up the calculations, a new algorithm for the integration of Newton's equation of motion of the dynamics of the molecules. Simple, stable and accurate, this algorithm can be used to simulate the motion of very complex molecules. At that time, I did not realise that the algorithm I had suggested was nothing but the generalization of the geometric algorithm of Newton ... ».

In Fig. 1, we show this geometric algorithm extracted from Newton's "Principia" [Newton ed. 1883]: The particle's discrete trajectory -→ r (t) is the broken line A, B, C, D, E, F, ... and the source of the force acting on the particle is at S. As Loup Verlet notes on page 124 of "Chimères et Paradoxes", even if « the translation of proportionality ratios into vector equalities is anachronistic, it reflects Newton's indications », (cf. [Principia Translation 1934]). In particular the AB segment can be interpreted as -→ v (t -dt) dt or -→ r (t) --→ r (t -dt), the segment BV as -→ f (t) dt 2 /m, etc. This translation into vectorial quantities, and noting that ABCV is a parallelogram, leads to the relation:

Fig. 1 .

 1 Fig. 1. From a facsimile edition of Newton's Philosophiae naturalis principia mathematica, schema of a particle trajectory submitted to a central force [Newton ed., 1883], cf. text below.

Acknowledgements

This paper is dedicated to Loup Verlet, the central actor in our presentation, who was our mentor and thesis supervisor in the sixties. We are grateful to Giovanni Ciccotti, Gianni Battimelli and Rodolphe Vuilleumier, who organised a very fruitful workshop in Paris, bringing together computational scientists and historians of science in May 2018.