Emanuela Abbate
email: emanuela.abbate@inria.fr

Angelo Iollo
email: angelo.iollo@inria.fr

Gabriella Puppo
email: gabirella.puppo@uninsubria.it

An implicit scheme for moving walls and multi-material interfaces in weakly compressible materials

Keywords: AMS subject classifications: 65M08, 65M99, 76T15 multi-material interfaces, moving boundaries, low Mach regime, Eulerian elasticity

We propose a numerical method for the simulation of flows from weakly compressible to low Mach regimes in domains with moving boundaries. Non-miscible weakly compressible materials separated by an interface are included as well. The scheme is fully implicit and it exploits the relaxation all-speed scheme introduced in [1]. We consider media with significantly different physical properties and constitutive laws, as fluids and hyperelastic solids. The proposed numerical scheme is fully Eulerian and it is the same for all materials. We present numerical validations by simulating weakly compressible fluid/fluid, solid/solid and solid/fluid interactions.

Introduction

Physical phenomena involving different materials arise in several applications such as multiphase flows, fluid-structure interaction and impacts. A possible simulation approach for these phenomena is the modelisation of the material discontinuity with a diffused interface [START_REF] Abgrall | How to prevent pressure oscillations in multicomponent flow calculations: a quasi conservative approach[END_REF][START_REF] Karni | Multicomponent flow calculations by a consistent primitive algorithm[END_REF]. This method is robust and prevents the formation of oscillations, thanks to an interface that corresponds to an artificial mixture of the two fluids. However, the interface amplitude is due to the numerical diffusion, which increases with time evolution.

The interface can also be considered as a contact discontinuity, especially in presence of interactions among different materials such as solid/solid and solid/fluid, or also between two non-miscible fluids. In this framework, the ghost fluid method has been introduced in [START_REF] Fedkiw | A non-oscillatory eulerian approach to interfaces in multimaterial flows (the ghost fluid method)[END_REF] to model multiphase flows with sharp interfaces on Cartesian meshes. It consists in considering in every cell of the domain both the real fluid and the "ghost fluid", extrapolating from this latter some quantities at the interface. Other variants and extensions of the ghost fluid method have been proposed in several works [START_REF] Abgrall | Computations of compressible multifluids[END_REF][START_REF] Liu | Ghost fluid method for strong shock impacting on material interface[END_REF][START_REF] Barton | An Eulerian method for multi-component problems in non-linear elasticity with sliding interfaces[END_REF][START_REF] Barton | A conservative level-set based method for compressible solid/fluid problems on fixed grids[END_REF][START_REF] Zeng | A systematic approach for constructing higher-order immersed boundary and ghost fluid methods for fluid-structure interaction problems[END_REF]. Immersed boundary methods [START_REF] Hu | On the HLLC Riemann solver for interface interaction in compressible multi-fluid flow[END_REF][START_REF] Luo | An immersedboundary method for flow-structure interaction in biological systems with application to phonation[END_REF][START_REF] Gorsse | A simple Cartesian scheme for compressible multimaterials[END_REF] are another option to keep the material interface sharp. With these methods, the contact discontinuity representing the interface can arbitrarily cross the grid and the transmission conditions are applied via interpolation. Several strategies have been proposed in literature to solve both moving boundaries and multi-material flows in the fully compressible regime, including Lagrangian models [START_REF] Scovazzi | Stabilized shock hydrodynamics: I. A Lagrangian method[END_REF][START_REF] Fazio | Central schemes and second order boundary conditions for 1D interface and piston problems in Lagrangian coordinates[END_REF] and Eulerian models [START_REF] Barton | Eulerian adaptive finite-difference method for high-velocity impact and penetration problems[END_REF][START_REF] Favrie | Solid-fluid diffuse interface model in cases of extreme deformations[END_REF][START_REF] Miller | A conservative three-dimensional Eulerian method for coupled solidfluid shock capturing[END_REF][START_REF] Chertock | A simple eulerian finite-volume method for compressible fluids in domains with moving boundaries[END_REF].

In the present work we are concerned with the numerical simulation of multi-material flows in the weakly compressible regime and especially when the Mach number tends to zero. We address the solution of fluid-dynamics problems and also compressible solids deformations. The propagation of waves in heterogeneous compressible media can be affected by drastic changes in the speed of sound or in the speed of elastic waves. This is related to the different local stiffness of the considered materials. Thus, it is important to have a scheme that is able to deal with different regimes. Standard multi-material schemes are usually based on the explicit-upwind framework for the simulation of compressible flows. Such schemes may fail when approaching the low Mach number limit, since upwind discretizations provide an excessive numerical viscosity on the slow waves in this regime [START_REF] Guillard | On the behaviour of upwind schemes in the low Mach number limit[END_REF][START_REF] Dellacherie | Analysis of Godunov type schemes applied to the compressible Euler system at low Mach number[END_REF]. The other drawback in adopting explicit methods for low Mach flows is related to the need of extremely small time steps, because in such schemes the CFL condition has to be imposed on the fastest wave speed to ensure stability. Due to this, the time step of compressible codes becomes extremely small as the incompressible regime is approached, increasing the computational time. In literature, several different schemes have been devised to deal with low Mach regimes, including implicit preconditioning methods [START_REF] Li | An all-speed Roe-type scheme and its asymptotic analysis of low Mach number behaviour[END_REF][START_REF] Viozat | Implicit upwind schemes for low Mach number compressible flows[END_REF][START_REF] Turkel | Preconditioned methods for solving the incompressible and low speed compressible equations[END_REF] and semi-implicit asymptotic preserving schemes [START_REF] Cordier | An asymptotic-preserving all-speed scheme for the Euler and Navier-Stokes equations[END_REF][START_REF] Noelle | An asymptotic preserving all Mach number scheme for the Euler equations of gas dynamics[END_REF][START_REF] Haack | An all-speed asymptotic-preserving method for the isentropic Euler and Navier-Stokes equations[END_REF]. All these methods are specifically build to solve fluid flows.

We are able to overcome the problems mentioned above by coupling an implicit treatment of the interface with the all-speed scheme proposed in [START_REF] Abbate | An all-speed relaxation scheme for gases and compressible materials[END_REF] for the simulation of compressible materials, including elastic solids. The scheme is general enough to deal with different materials, since it is based on a monolithic Eulerian model, describing each material (gas, liquid or solid) with the same system of conservation laws and an appropriate general formulation of the constitutive law [START_REF] Godunov | Elements of continuum mechanics and conservation laws[END_REF][START_REF] Plohr | A conservative Eulerian formulation of the equations for elastic flow[END_REF][START_REF] Plohr | A conservative formulation for plasticity[END_REF]. The scheme is based on the Jin-Xin relaxation and it is fully implicit. In the present paper, we propose a method to impose boundary conditions on moving walls (piston problem) and a technique to solve flows across multi-material interfaces, even in the case when the interface velocity is unknown and determined by the flow. In particular, we are interested in an all-speed treatment in order of accurately approximating moving boundaries and multi-material interfaces in the weakly compressible and low Mach regimes. For this reason the treatment of the boundary and of the interface must be fully implicit, as in the all-speed scheme [START_REF] Abbate | An all-speed relaxation scheme for gases and compressible materials[END_REF] used in the bulk of the flow. To our knowledge, this is the first completely implicit scheme to deal with moving interfaces in compressible and low Mach non-viscous flows.

Firstly, we show how to implement conditions for moving boundaries when the ve-locity of the wall is externally imposed, hence known at all times. Some examples can be piston problems or beam elongations/compressions. In literature, several explicitupwind schemes have been extended to the simulation of moving boundaries in the fully compressible regime. In the present work, our aim is to approach the approximation of moving walls also when the Mach number tends to zero. Thus, also the boundary treatment has to be implicit and has to preserve the asymptotic behaviour of the solution.

Then, the multi-material model is introduced. The interface between different materials is treated as a moving wall and suitable equilibrium conditions are imposed via extrapolations, as in immersed boundary methods. This method accurately approximates multi-material flows with Mach number ranging from the incompressible limit to compressible regimes with M O(1). The evolution of the material discontinuity is sharp by construction and it stays sharp even after long times.

In this paper, we focus on one dimensional simulations, where most difficulties already occur. The numerical results show convergence of the schemes and a perfect resolution of the interfaces. A comparison of different discretizations of the advective operator (upwind, centered and all-speed) is carried out. It is observed that the proposed wall/interface conditions can be accurately imposed no matter how the spatial derivatives are discretized. The ability of the multi-material scheme in approximating the physical interface is exhibited in both compressible and low Mach regimes.

The outline of the paper is as follows. In Section 2, we describe the Eulerian model for simulating compressible materials and in Section 3 we revise the relaxation all-speed scheme introduced in [START_REF] Abbate | An all-speed relaxation scheme for gases and compressible materials[END_REF]. Section 4 is devoted to the proposed numerical model for walls and multi-material interfaces, which is included in the implicit all-speed scheme. The numerical validation of the model is performed in Section 5, where its applicability to different regimes and materials is shown. Finally, conclusions are drawn in Section 6.

Eulerian model for compressible materials

To fix notations, we briefly revise the monolithic Eulerian model discussed in several works [START_REF] Gorsse | A simple Cartesian scheme for compressible multimaterials[END_REF][START_REF] Godunov | Elements of continuum mechanics and conservation laws[END_REF][START_REF] Plohr | A conservative Eulerian formulation of the equations for elastic flow[END_REF][START_REF] Plohr | A conservative formulation for plasticity[END_REF][START_REF] Brauer | A Cartesian scheme for compressible multimaterial models in 3d[END_REF]. Different compressible materials are described with the same system of conservation laws. The one dimensional version and the low Mach limits have been analized in [START_REF] Abbate | An all-speed relaxation scheme for gases and compressible materials[END_REF].

General 2D model

When describing a continuous medium, let Ω 0 ∈ R 2 be the reference or initial configuration and Ω t ∈ R 2 the deformed configuration at time t. We define the forward characteristics X (ξ,t) as the image at time t in the deformed configuration of a material point ξ belonging to the initial configuration, i.e., X :

Ω 0 ×[0,T] → Ω t , (ξ,t) → X (ξ,t). The corre- sponding Eulerian velocity field u : Ω t ×[0,T] → R 2 , (x,t) → u(x,t) corresponds to ∂ t X (ξ,t) = u(X (ξ,t),t) X (ξ,0) = ξ, ξ ∈ Ω 0 .
If one aims at describing the continuum in the Eulerian framework,the backward characteristics Y(x,t) associate to a material point in position x in the deformed configuration at time t, the corresponding position ξ in the initial configuration of the same material point, i.e., Y : Ω t ×[0,T] → Ω 0 , (x,t) → Y(x,t). When differentiating relation Y(X (ξ,t),t) = ξ with respect to time and space we get:

∂ t Y+u•∇ x Y = 0 Y(x,0) = x, x ∈ Ω t . (2.1)
In general, the stress tensor has a direct dependence on the gradient of the deformation, which can be computed in the Eulerian framework via Y. Thus, the gradient of (2.1) is taken as a governing equation. The other governing equations are given by the conservation of mass, momentum and energy in the Eulerian framework (for details, see [START_REF] Gorsse | A simple Cartesian scheme for compressible multimaterials[END_REF][START_REF] Brauer | A Cartesian scheme for compressible multimaterial models in 3d[END_REF]). In two dimensions, the conservative form of the equations of a general medium in the deformed configuration reads

           ∂ t ρ+div x (ρu) = 0 ∂ t (ρu)+div x (ρu⊗u-σ) = 0 ∂ t ([∇ x Y])+∇ x (u•[∇ x Y]) = 0 ∂ t (ρe)+div x ρeu-σ T u = 0. (2.2)
Here ρ is the density, u is the Eulerian velocity field, [∇ x Y] is the gradient of the backward characteristics and σ is the Cauchy stress tensor. e is the total energy per unit mass and it is given by the sum of the kinetic energy and the internal energy per unit mass :

e = 1 2 |u| 2 + .
The adopted constitutive law has a general formulation that includes different behaviours and describes gases, fluids and elastic solids at the same time. The state law reads [START_REF] Gorsse | A simple Cartesian scheme for compressible multimaterials[END_REF][START_REF] Brauer | A Cartesian scheme for compressible multimaterial models in 3d[END_REF]

(ρ,s,[∇ x Y]) = κ (s) γ-1 ρ γ-1 + p ∞ ρ gas + χ ρ trB-2 neohookean solid , (2.3)
where κ (s) = exp(s/c v ) (s being the entropy), χ, p ∞ , γ are positive constants that char- acterize a given material and B is the 2D right Cauchy-Green tensor, which is defined as follows The Cauchy stress tensor σ has the following formulation given by (see [START_REF] Holzapfel | Nonlinear solid mechanics: a continuum approach for engineering science[END_REF] for the derivation):

B = [∇ x Y] -1 [∇ x Y] -T , J = det[∇ x Y] -1 , B = B J . ( 2
     σ(ρ,s,[∇ x Y]) = -p(ρ,s) I +2χJ -1 B- trB 2 I p(ρ,s) = -p ∞ +k(s)ρ γ .
(2.5)

The model in one dimension

We reduce the 2D model by considering only variations that occur along direction x 1 .

Deformations in both the two directions are taken into account , but the derivatives along x 2 are set equal to zero. This way, tensor [∇ x Y] reduces to

[∇ x Y] = Y 1 ,1 0 Y 2 ,1 1
.

The equation on Y 1 ,1 is actually equivalent to conservation of mass, due to the fact that ρ(x,t) = det([∇ x Y](x,t))ρ 0 (x), where ρ 0 is the initial density. Thanks to this we have Y 1

,1 = ρ/ρ 0 and one equation is redundant. Moreover, the Cauchy stress tensor (2.5) reduces to only two non-zero components, which are the normal and the tangential stress respectively. They take the following formulation:

σ 11 =-p(ρ,s)+2χJ -1 B 11 - trB 2 = -p(ρ,s)+χ 1-Y 2 ,1 2 -(ρ/ρ 0 ) 2 σ 21 =2χJ -1 B 21 = -2χY 2 ,1 .
We write a compact formulation of the 1D Eulerian model as

∂ t ψ+∂ x 1 F(ψ) = 0, (2.6)
where the vector of conservative variables ψ and the vector of fluxes F(ψ) are both in R 5 and read

ψ =       ρ ρu 1 ρu 2 Y 2 ,1 ρe       , F(ψ) =       ρu 1 ρu 2 1 -σ 11 ρu 1 u 2 -σ 21 u 1 Y 2 ,1 +u 2 ρe-σ 11 u 1 -σ 21 u 2       .
The Mach number in fluid-dynamics is calculated as M=u 1 /c, where the sound speed is defined as c = ∂p/∂ρ| s=const = γ(p+ p ∞ )/ρ. When the Mach number tends to zero, system (2.6) becomes stiff due to the acoustic part of the stress tensor gradient, namely to the pressure part. Moreover, if the material properties are such that O (χ) O (p ∞ ), the stiffness in the low Mach limit is also due to the elastic part. In this case, also the "elastic Mach number" M χ = ρu 2 1 /(2χ) tends to zero. The analysis of the different low Mach limits that can occur in this monolithic Eulerian framework has been carried out in [START_REF] Abbate | An all-speed relaxation scheme for gases and compressible materials[END_REF].

Implicit all-speed scheme

We present the all-speed relaxation scheme recently introduced in [START_REF] Abbate | An all-speed relaxation scheme for gases and compressible materials[END_REF][START_REF] Abbate | A relaxation scheme for the simulation of low Mach number flows[END_REF]. The scheme is the same for all the considered compressible materials, since it does not have a direct dependence on the state law. It is based on the Jin-Xin relaxation [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimensions[END_REF] in order to have a linear advective operator. We consider the hyperbolic system of conservation laws (2.6) with the vector of conservative variables ψ ∈ R n , n = 5. By introducing the vector of relaxation variables v ∈ R n , the relaxation system of dimension 2n reads as follows:

   ∂ t ψ+∂ x 1 v = 0 ∂ t v+A∂ x 1 ψ = 1 η (F(ψ)-v). (3.1)
The parameter η > 0 is the relaxation rate, a small parameter that in the present work is always taken equal to 10 -8 . The relaxation matrix A is chosen by enforcing the subcharacteristic condition [START_REF] Whitham | Linear and nonlinear waves[END_REF]:

A-F (ψ) 2 ≥ 0. (3.2)
More details on the construction of the relaxation matrix are given in [START_REF] Abbate | An all-speed relaxation scheme for gases and compressible materials[END_REF].

Implicit time discretization

The acoustic CFL constraints become computationally demanding when the low Mach limit approaches. The scheme is fully implicit, in order to get rid of the stability constraint. One of the advantages of the relaxation procedure consists in the linearity of the advective operator in (3.1). Thanks to this, the introduction of Riemann solvers may be avoided and this allows for a straightforward use of implicit time discretizations. At first order, a simple backward Euler scheme is introduced:

       ψ n+1 -ψ n ∆t +∂ x 1 v n+1 = 0 v n+1 -v n ∆t +A∂ x 1 ψ n+1 = 1 η F ψ n+1 -v n+1 (3.3)
The non-linear fluxes F(ψ) in the stiff relaxation term are solved with a Newton method. Actually computational experiments have shown that one iteration is always enough. Thus, the resulting approximation is the following truncated Taylor expansion:

F ψ n+1 = F(ψ n )+F (ψ n ) ψ n+1 -ψ n , (3.4)
where F (ψ n ) is the Jacobian of the flux and is computed analytically. Since the fully implicit time discretization is unconditionally stable, the time step has to be chosen only for accuracy reasons. We define a Courant number in general as ν = ∆t/∆x•λ, where λ is a wave speed of the problem. If we need accuracy on the acoustic waves, we enforce an acoustic constraint

∆t ≤ ν ac ∆x λ max , (3.5)
where λ max is the fastest characteristic speed, namely the speed of the acoustic waves.

Choosing ∆t as in (3.5), the time step is small enough to follow all the propagating waves: specifically, the acoustic waves are fully resolved. If instead we only want to focus on the accurate solution of material waves, the time step can be chosen with a material condition

∆t ≤ ν mat ∆x λ mat , (3.6)
where λ mat is the material wave speed. This way, the computational time is reduced when approaching the low Mach limit. With this constraint, the material waves are accurately resolved, while the fast longitudinal and shear waves are under resolved. This may not represent an issue in low Mach regimes, since usually these waves carry a small amount of energy and the main interest is the flow velocity approximation (or the deformation inside solids). For more details on the accuracy on the different waves, we refer to [START_REF] Abbate | An all-speed relaxation scheme for gases and compressible materials[END_REF].

Spatial discretization

System (3.1) is discretized with finite volumes on a Cartesian mesh. A 1D computational domain

[0,L] of length L is divided into N cells C i = [x i-1/2 ,x i+1/2 ] i = 1,.., N, letting ∆x = x i+1/2 -x i-1/2
be the grid spacing. w i denotes the approximate cell average of a quantity w in the cell C i and w i+1/2 denotes the approximate point value of w at the cell interface x = x i+1/2 . The semi-discretization in space for system (3.1) reads

     ∂ t ψ i + v i+1/2 -v i-1/2 ∆x = 0 ∂ t v i +A ψ i+1/2 -ψ i-1/2 ∆x = 1 η (F(ψ i )-v i ). (3.7)
The variables at the interfaces x i+1/2 inside (3.7) have to be computed with an "ad hoc" stable scheme.

It has been proved in previous works [START_REF] Guillard | On the behaviour of upwind schemes in the low Mach number limit[END_REF][START_REF] Dellacherie | Analysis of Godunov type schemes applied to the compressible Euler system at low Mach number[END_REF] that upwind spatial discretizations present an excess of numerical viscosity in the case of low Mach flows. Instead, centered spatial discretizations of the pressure gradient in the Euler equations limit the viscosity on any mesh [START_REF] Dellacherie | Analysis of Godunov type schemes applied to the compressible Euler system at low Mach number[END_REF][START_REF] Dellacherie | Construction of modified Godunov-type schemes accurate at any Mach number for the compressible Euler system[END_REF] and the low Mach regime is accurately solved. For the adopted Eulerian model, this corresponds to centering the stiff parts of the stress tensor σ, as it has been proposed in [START_REF] Abbate | An all-speed relaxation scheme for gases and compressible materials[END_REF].

A hybrid spatial discretization is introduced with the use of a convex combination of upwind and centered schemes. An upwind spatial discretization is desired when M O (1) in order to prevent spurious oscillations. On the other hand, in the limit M → 0, a centered scheme is required for accuracy. The combination of these two schemes is based on the local Mach number of the specific flow. This allows to build an all-speed scheme. For a general variable w, the interface values are computed as follows:

w i+1/2 = f (M loc )(w i+1/2 ) upw +(1-f (M loc ))(w i+1/2 ) cent , (3.8)
where M loc is the local Mach number at the previous time step. The function f (M loc ) has to satisfy the criterion 0 ≤ f (M loc ) ≤ 1. We simply choose a minimum function f (M loc ) = min{1, M loc }. The interface values with a centered scheme are computed in the following way:

(w i+1/2 ) cent = 1 2
(w i+1 +w i ).

(3.9)

It has been shown in [START_REF] Abbate | An all-speed relaxation scheme for gases and compressible materials[END_REF] that in the low Mach limit this discretization is equivalent to centering the whole fluxes. Therefore, in this case also the gradient of σ is centered and consequently the accuracy in space does not depend on the Mach number. The upwind scheme is built as in [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimensions[END_REF] along the two characteristic variables v±A 1/2 ψ of the relaxation system (3.1). The interface values obtained with the upwind discretization are the following:

ψ i+1/2 upw = 1 2 ψ i+1 +ψ i -A -1/2 2 (v i+1 -v i ) (v i+1/2 ) upw = 1 2 (v i+1 +v i )-A 1/2 2 ψ i+1 -ψ i .
(3.10) Discretizations (3.9) and (3.10) can be plugged into the convex combination (3.8) along both directions, obtaining the following "all-speed" interface values:

ψ i+1/2 all = 1 2 ψ i+1 +ψ i -f (M loc ) 2 A -1/2 (v i+1 -v i ) (v i+1/2 ) all = 1 2 (v i+1 +v i )- f (M loc ) 2 A 1/2 ψ i+1 -ψ i . (3.11)
Central differencing has potentially no numerical diffusion in space, unlike upwinding. However, a centered scheme is not stable in an explicit time-integration. Here we are allowed to adopt this scheme since we have introduced the fully implicit framework (3.3). Moreover, the convex combination moderates the numerical viscosity as M loc tends to zero. This is clear by comparing the all-speed discretization (3.11) with the upwind one (3.10): in both schemes a diffusive term is present, but in the all-speed discretization it is reduced via the multiplication by f (M loc ) ≤ 1 (for a detailed analysis see [START_REF] Abbate | An all-speed relaxation scheme for gases and compressible materials[END_REF]).

Then, the all-speed implicit scheme has the following formulation:

               ψ n+1 i -ψ n i ∆t + 1 2∆x v n+1 i+1 -v n+1 i-1 - f (M loc )A 1/2 2∆x ψ n+1 i+1 -2ψ n+1 i +ψ n+1 i-1 = 0 v n+1 i -v n i ∆t + A 2∆x ψ n+1 i+1 -ψ n+1 i-1 - f (M loc )A 1/2 2∆x v n+1 i+1 -2v n+1 i +v n+1 i-1 = 1 η F ψ n+1 i -v n+1 i (3.12)

Numerical method for moving walls and multi-material interfaces

A numerical method for the simulation of moving boundaries and multi-material interfaces is derived. With the aim of solving weakly compressible and low Mach multimaterial flows, we propose an implicit treatment of the wall/interface, involving a layer of cells close to the interface. This is coupled with the implicit relaxation scheme with general spatial discretization (3.7), which is solved in the bulk of the flow, where upwind (3.10), centered (3.9) and all-speed (3.8) discretizations of the transport operator may be used.

General model for walls and interfaces

We introduce a general framework to model both walls moving at prescribed velocity and physical interfaces where the velocity is dictated by the flow. The 1D domain [0,L] is discretized as described in Sec. 3.2. Let the wall/interface be in position x k ≤ x B (t) < x k+1 at time t with material 1 on the left side of the domain [0,x B ] and material 2 on the right side [x B ,L]. We distinguish among internal cells (of each material) and interface cells. The internal cells are fully occupied by one specific material, the interface cells form a thin layer between the two materials. The interface cells have to be introduced since the moving wall/interface cannot in general be forced to coincide with the cell edges. As a result, the interface cells are partially filled with one material and partially with the other one. This can be inconvenient in a finite volume logic, because numerical solutions are represented in terms of the cell averages. We point out that in 1D we have only one interface cell for each considered wall or physical interface.

In the present work, we assign to material 1 cell C k and to material 2 cell C k+1 , as it is depicted in Fig. 1(a). This means that the wall/interface is "artificially" set to coincide with the closest numerical interface, which is x k+1/2 in Fig. 1(a). This produces a numerical error of order O (∆x). Only when the wall/interface overcomes a cell center, this cell changes the material to which it is assigned. This approximation is consistent with the first order scheme (3.12).

Since we are dealing with two different materials, we need to introduce two different interface values at x k+1/2 , one from the left and one from the right. For a generic variable w, the interface value at the left of the interface (from material 1) will be w - k+1/2 and the value at the right of the interface (from material 2) will be w + k+1/2 , as depicted in Fig. 1(b). We have then

ψ k+1/2 -= ψ k+1/2 + and (v k+1/2 ) -= (v k+1/2 ) + . (4.1)
The implicit time discretization of (3.7) is therefore modified for cells k and k+1 by distinguishing among left and right interface values as follows:

             ψ n+1 k -ψ n k ∆t + v n+1 k+1/2 - -v n+1 k-1/2 ∆x = 0 v n+1 k -v n k ∆t +A ψ n+1 k+1/2 - -ψ n+1 k-1/2 ∆x = 1 η F ψ n+1 k -v n+1 k , (4.2)
             ψ n+1 k+1 -ψ n k+1 ∆t + v n+1 k+3/2 -v n+1 k+1/2 + ∆x = 0 v n+1 k+1 -v n k+1 ∆t +A ψ n+1 k+3/2 -ψ n+1 k+1/2 + ∆x = 1 η F ψ n+1 k+1 -v n+1 k+1 .
(4.3)

v n+1 k-1/2 , v n+1
k+3/2 and the corresponding conservative variables are computed with the chosen scheme, which can be upwind (3.10), centered (3.9) or all-speed (3.11). The interface material 1 material 2 variables (v k+1/2 ) -, (v k+1/2 ) + , ψ k+1/2 -and ψ k+1/2 + need to be reconstructed with ad hoc conditions. In all other cells C i , i = k,k+1, scheme (3.12) is solved. Since we are using a fully implicit scheme, it is more practical to avoid the introduction of ghost cells to impose wall and multi-material interface conditions. Ghost cells would consistently increase the dimension of the matrix and the memory use, especially in multi-dimensions. This explains why we work directly on the interface x k+1/2 and we introduce left and right reconstructions to be used directly inside discretizations (4.2)-(4.3).

x k x k+1 x k+1/2 interface region C k C k+1 x B (a) Wall/Interface position material 1 material 2 x k x k+1 w -k+1/2 w + k+1/2 x k+1/2 (b) Interface values
Specific conditions for reconstructing the interface values have to be derived. In this perspective, we firstly approach the numerical modeling of walls moving at an externally imposed velocity. In this case, the wall velocity is known at all times, whereas the interface values of the other quantities need to be reconstructed. As a second step, we derive equilibrium conditions for the simulation of interfaces between different materials. In this latter case, the velocity of the interface is unknown, being dictated by the interaction between the two materials.

Conditions for moving walls with imposed velocity

We consider walls separating two media that do not interact. Within this framework, pistons in fluids that expand or compress according to the piston motion and elongations or compressions of elastic beams can be simulated.

We derive the numerical conditions that have to be imposed to reconstruct the interface values (v k+1/2 ) -, (v k+1/2 ) + , (ψ k+1/2 ) -and (ψ k+1/2 ) + for the case of a wall that moves with a prescribed velocity. Let u * 1 be the velocity imposed to the wall and known at all times. The position of the wall at time t n+1 is computed as x B (t n+1 )=x B (t n )+u * 1 •∆t. In the case of Euler gas-dynamics equations, the density and the pressure are extrapolated from the left and from the right, having:

           ρ n+1 k+1/2 - = 3 2 ρ n+1 k -1 2 ρ n+1 k-1 u n+1 k+1/2 - = u * 1 t n+1 p n+1 k+1/2 - = 3 2 p n+1 k -1 2 p n+1 k-1 and            ρ n+1 k+1/2 + = 3 2 ρ n+1 k+1 -1 2 ρ n+1 k+2 u n+1 k+1/2 + = u * 1 t n+1 p n+1 k+1/2 + = 3 2 p n+1 k+1 -1 2 p n+1 k+2 . (4.4)
By considering the full Eulerian framework (2.2), we also need to extrapolate the backward characteristic Y 2 ,1 from each material and to impose the transverse velocity of the wall u * 2 . Thus, ((Y 2 ,1 ) n+1 k+1/2 ) ± share the same discretization of (ρ n+1 k+1/2 ) ± in (4.4). Then, these conditions have to be formulated on conservative and relaxation variables, in order to be included in (4.2)-(4.3) with the correct formulation. Conditions on the momentum are easily imposed with u * 1 and the density extrapolations. For the energy, we use the pressure extrapolations and the state law. Due to the fully implicit time discretization, a linearization with respect to the conservative variables is applied to the non-linear terms. The relaxation variables are written as

   v n+1 k+1/2 - = F ψ n+1 k+1/2 - v n+1 k+1/2 + = F ψ n+1 k+1/2 + , (4.5)
where again a linearization has to be applied to F(ψ).

The model and the conditions derived above solve the case of internal moving walls, separating two non-interacting materials. Of course, the model can be easily adapted to the simulation of a moving boundary wall, where the focus is on the piston or the beam only, and the surrounding medium is not considered. In this case, the right part of the domain is considered empty.

Linear system for moving walls

By defining the vector containing the grid point values of conservative and relaxation variables as w = [ψ,v] ∈R N•2n (N is the number of cells and n the number of conservative variables), the full linear problem can be written in the following way:

M(w n ,u * 1 ,u * 2 )w n+1 = r(w n ,u * 1 ,u * 2 ) (4.6)
where M ∈ R (N•2n)×(N•2n) and r ∈ R N•2n . The matrix structure comes from the spatial discretization introduced in Sec. 3 and from the linearized conditions. The matrix and the right hand side are both functions of the wall velocity u * 1 (and u * 2 when it is non zero). The two materials are not interacting with each other, hence the full system can be split into two sub-systems, each one associated to the corresponding material. The size of the two sub-systems varies with k every time the wall overcomes a cell center.

Conditions for multi-material interfaces

We extend the model to the case of physical interfaces separating interacting materials. We propose a "sharp interface" model, since the interface is considered as a moving wall and thus treated with the framework introduced above. The velocity of the wall is now dictated by the flow, thus becoming one of the unknowns of the problem. Therefore, additional equilibrium conditions have to be introduced to account for the interactions between the two materials and to recover the correct velocity of the interface. Due to this, the linear system cannot be split into two different sub-systems associated to the two materials.

It is physically required that the forces are balanced at the multi-material interface. For simplicity, we start by reasoning on gas/gas interactions. In this case, a force balancing simply corresponds to a pressure balancing, namely p - k+1/2 = p + k+1/2 . With the perfect gas state law, one gets:

γ --1 (ρe) - k+1/2 - 1 2 ρ - k+1/2 (u * 1 ) 2 = γ + -1 (ρe) + k+1/2 - 1 2 ρ + k+1/2 (u * 1 ) 2 , (4.7)
where γ -and γ + are the heat ratios for material 1 and material 2 respectively and u * 1 is the unknown velocity of the interface. The balancing relation (4.7) may be used to find u * 1 : this way, system (4.6) for the moving wall can be closed. However, the matrix of the full moving wall system depends on u * 1 itself. Moreover, ρ ± k+1/2 and (ρe) ± k+1/2 in (4.7) are unknowns, because the scheme is implicit. This procedure would then require a fixed point iteration to find the correct interface velocity at every time step.

In order to avoid this problem, trasmission conditions at interface are imposed in terms of pressure and velocity via interpolations, as in immersed boundary methods [START_REF] Peskin | The immersed boundary method[END_REF]. This means that, for fluid-dynamics problems, we prescribe continuity of velocity and pressure across the interface in the numerical scheme. Specifically, this is done by computing a mean of the values on the left and on the right of the interface for continuous quantities:

   u n+1 k+1/2 - = u n+1 k+1/2 + = 1 2 u n+1 k +u n+1 k+1 p n+1 k+1/2 - = p n+1 k+1/2 + = 1 2 p n+1 k + p n+1 k+1 (4.8)
This way, the force balancing (4.7) is automatically taken into account by the numerical model. Other quantities such as density and energy may be discontinuous across the interface, hence we compute them with an extrapolation from each material (left and right). The density extrapolation is done as in (4.4) for the case of the moving wall.

We extend the interface treatment also to the simulation of solid/solid and fluid/solid interactions. The multi-material conditions are enforced on the full Eulerian system (2.6), where the stress tensor σ takes the role that pressure has in gas-dynamics. When dealing with an interface between two elastic solids, normal and tangential stress are continuous across the contact discontinuity. Moreover, also both velocity components are continuous. Instead, in the case where at least one of the materials is a fluid, namely χ =0 on one side, the tangential stress σ 21 vanishes at the interface and thus the transverse velocity can be discontinuous.

Interface conditions to be imposed for every kind of interaction (fluid/fluid, solid/solid and solid/fluid) consist in the left/right extrapolation of the density (ρ n+1 k+1/2 ) ± and of the backward characteristics ((Y 2 ,1 ) n+1 k+1/2 ) ± as in (4.4). Then, conditions ((u 1 ) n+1 k+1/2 ) -= ((u 1 ) n+1 k+1/2 ) + and (σ 11 k+1/2 ) -= (σ 11 k+1/2 ) + are imposed with the mean computation (4.8). To these relations, we need to add conditions on the remaining variables, distinguishing two cases:

1. solid/solid interface: (σ 21 k+1/2 ) -= (σ 21 k+1/2 ) + and ((u 2 ) n+1 k+1/2 ) -= ((u 2 ) n+1 k+1/2 ) + are computed with the mean as in (4.8); 2. solid/fluid interface: we impose that (σ 21 k+1/2 ) -= (σ 21 k+1/2 ) + = 0 and ((u 2 ) n+1 k+1/2 ) ± is computed with the left/right extrapolation as in (4.4).

We remark that the interface conditions are written on primitive variables, to impose the continuity of u and σ. These relations need to be transferred on the conservative variables ψ, in order to include them inside the implicit relaxation scheme (4.2)-(4.3). The relation between primitive and conservative variables is non-linear, therefore a Newton method could be used to deal with the non-linearities. The linearization of the reported multi-material conditions with respect to the conservative variables corresponds to one Newton iteration, as in (3.4). The linearization itself is usually enough, however, for cases when the Mach number gets of order one or for wave patterns with shocks that are close to the interface, running the Newton sub-iteration is required to recover the correct speed and position of the interface.

In what follows, we will adopt the interface method (4.2)-(4.3) with the linearized multi-material conditions, which is our "multi-material scheme". This can be coupled with different spatial discretization away from interfaces, such as upwind (3.10), centered (3.9) and all-speed (3.12) schemes. As for ghost-fluid methods, the scheme is locally nonconservative. In the numerical tests section the scheme proves to be consistent, correctly predicting shock speeds and positions. Indeed, the number of cell interfaces for which a non-conservative numerical flux is employed is always negligible compared to the total number of mesh cells.

Multi-material algorithm

We summarize the algorithm for the multi-material scheme, describing one time step. Let the interface be in position x B (t n ) at time t n , having x k ≤ x B (t n ) ≤ x k+1 . A material CFL condition on ∆t is always enforced, thus the physical interface can cross at most one cell interface at every time step. The following procedure is implemented from time t n to time t n+1 : 1. the physical interface is advected, computing its position with the velocity from the previous time step as follows:

x B t n+1 = x B (t n )+ (u 1 ) n k +(u 1 ) n k+1 2 ∆t.
2. the new position is located on the grid by finding the cell index k such that x k ≤ x B (t n+1 ) ≤ x k+1 . We distinguish two cases:

• k = k, i.e. the multi-material interface has not overcome a cell center in the advection from time t n to time t n+1 . In this case the cells occupied by material 1 and the ones occupied by material 2 are the same as they were at time t n .

• k = k, i.e. the multi-material interface at time t n+1 has overcome one cell center with respect to time t n , either having k = k-1 or k = k+1. We examine the case C k = C k+1 (forward-moving interface): this cell at time t n+1 is assigned to a different material with respect to time t n . The state of cell C k is needed at time t n to assemble the matrix of the implicit scheme. In this cell, the state at time t n is known for material 2, but is now needed as initial condition for material 1. We "correct the initial condition" at the interface cell by creating a fictitious state at time t n for the cell that changes material at time t n+1 . To do this, we simply extrapolate the quantities from the material in expansion.

3. with all the computed quantities (interface position, assignment of the cells to the corresponding material, fictitious old state), we assemble the linear system associated to scheme (3.12)-(4.2)-(4.3), with the derived reconstructions for the interface conditions. Then the linear system is solved to obtain all variables at time t n+1 .

We point out that the accuracy of the proposed method largely depends on the computation of the interface velocity, via the conditions derived in Sec. 4.1.2. These conditions are built with the aim of solving low Mach and weakly compressible multi-material flows.

For high speed impacting problems (e. g. Mach numbers of O(10) or higher), the use of suitable multi-material Riemann solvers might be required to recover the correct interface velocity. As shown below with the numerical experiments, we analyze multimaterial flows ranging from the low Mach and weakly compressible regimes, to the fully compressible one with M O(1).

Numerical results

Moving walls with imposed velocity

The moving wall model is here validated with simulations of gas pistons and beams elongations. In all the simulations we apply at the left boundary of the domain homogeneous Neumann boundary conditions and the final times are taken sufficiently small so that no waves reach the left end. We present four gas piston problems in different regimes, the piston being a moving boundary. The initial state, the initial position of the wall x B (0) and the velocity u * 1 imposed at the wall are reported in Table 2. A biatomic gas with γ = 1.4 fills the entire domain. In tests 1, 2 and 3 the imposed velocity is constant in time, whereas in test 4 we study a piston with a sinusoidal trajectory to see expansion and compression. Test 3 and 4 present a discontinuity in the gas chamber at x disc (0) = 0.5. For test 1 and 2, the initial gas density and pressure are constant in space. The pistons move with a constant velocity u * 1 > 0, consequently the gas expands with a rarefaction wave moving to the left. We compare the results obtained with the three discretizations modified at the wall with ∆x = 2•10 -3 and ν mat = 0.2. Two regimes obtained with two different values of u * 1 are analyzed. In test 1, setting u * 1 = 0.5m/s produces a subsonic but fully compressible flow with M 0.4 on the piston. The piston moves from x B (0) = 0.5 to x B (t end ) = 0.6. In Fig. 2(a) we observe that the upwind discretization suffers from the overheating problem on the density profile (entropy errors occuring at the interface may cause overshoots/undershoots on the density and on the temperature: see [START_REF] Fedkiw | An isobaric fix for the overheating problem in multimaterial compressible flows[END_REF]). As expected, the overheating disappears when reducing numerical diffusion: this is why the centered discretization produces the smallest spike among the three. However, the centered scheme is oscillatory on the acoustic waves. The all-speed scheme produces the most accurate results, preventing oscillations and reducing the overheating spikes. It also reduces the diffusion on the head and tail of the rarefaction with respect to the upwind scheme. We get ∆t = 8•10 -4 with the chosen CFL, which corresponds to ν ac = 0.9 for this fully compressible simulation. In fact, to resolve the rarefaction wave an acoustic CFL constraint has to be enforced.

Gas pistons

Test L x B (0) x disc (0) t end ρ 1 ρ 2 (u 1 ) 1 (u 1 ) 2 p 1 p 2 u * 1 (m) (m) (s) (Kg/m 3 ) (Kg/m 3 ) (m/s) (m/s) (Pa) (Pa) (m/s) 1 1 0.5 - 0.2 1 - 0 - 1 - 0.5 2
In test 2, setting u * 1 = 0.01m/s gives a Mach number M 0.8•10 -3 on the moving boundary. The piston moves from x B (0) = 1.75 to x B (t end ) = 1.76, namely it crosses 5 grid cells. In Fig. 2(b), we observe that the rarefaction is smeared because the time step is too large to resolve acoustic waves in this regime (we have ∆t = 4•10 -2 ). The three schemes behave in a similar manner, only the centered scheme develops some small oscillations on the density profile at the wall. The mass conservation errors obtained with the all speed scheme for tests 1 and 2 are reported in Table 3. As expected, the mass conservation error decreases as the mesh is refined coherently with the truncation error, exhibiting 1.03 the numerical convergence of order one for test 1 and of almost order one for test 2. Moreover, these errors are about 10 times smaller with respect to the L 1 -norm error of the numerical solution, reported in Table 4. The mass loss introduced by the method may be controlled through a local mesh refinement, especially for 2D test cases. The results of the first two tests have shown that the implicit scheme with the upwind and with the all-speed scheme behave in a similar manner on acoustic waves and at the wall. However, they provide different results when approximating contact discontinuities in low Mach flows. This is clearly observed with test 3: here we solve a Riemann problem with a discontinuity initially in x disc (0) =0.5 and a moving boundary initially in x B (0) = 1. The chosen initial state produces M 6•10 -3 on the contact discontinuity. In Fig. 3 it is evident the all-speed scheme is able to keep the contact discontinuity sharper than the upwind scheme, because it provides the correct numerical viscosity on the material wave. This test illustrates the need of introducing a numerical model for moving walls inside an implicit scheme, which is the main point of the present work. At the end of the simulation, the wall position is x B (t end ) = 2.0048 (crossing 4 cells on the chosen grid) and the contact discontinuity position is x disc (t end ) = 1.0027 (crossing 2 cells). With test 4 we study a gas tube with a moving boundary wall oscillating around x B (0) = 0.9 with a time-dependent velocity u * 1 (t) = 0.5•πcos(10πt). The initial data correspond to a right moving wave, initially positioned in x disc = 0.5, with a small pressure ratio applied. The moving boundary generates shocks and rarefaction waves that interact with the incoming wave, creating a complicated solution structure. This test shows that a sinusoidal piston trajectory is well resolved in both expansion and compression phases. In Fig. 4 we compare the density profiles obtained with the upwind and all-speed spatial discretizations at different times: at time t = 0.05 seconds the wall position is in x B = 0.95 (maximum expansion) and at time t =0.15 seconds in x B =0.85 (maximum compression). The all-speed scheme keeps the waves structure slightly sharper with respect to the up-wind scheme, because the convex combination (3.8) is able to moderate the numerical viscosity. For this test we employ an acoustic CFL condition ν ac = 0.9 in order to resolve all the waves occurring during the simulation.

Two materials tests

Test L

x B (0 Table 5: Parameters for the two materials tests (copper/gas): initial state (subscripts: 1 for copper on the left and 2 for a biatomic gas on the right), initial position of the wall x B (0) and imposed velocity u * 1 . Tests 5 and 6 simulate the elongation of two copper beams surrounded by a perfect gas. The initial state is reported in Table 5 for the two tests. An external velocity is applied at the right boundary of the beam:

) t end ρ 1 ρ 2 (u 1 ) 1 (u 1 ) 2 p 1 p 2 u * 1 (m) (m) (s) (Kg/m 3 ) (Kg/m 3 ) (m/s) (m/s) (Pa) (Pa) (m/s) 5
• u * 1 =500m/s for test 5, leading to a subsonic but fully compressible regime in copper (M -0.14) and a Mach number of order 1 in the gas side. The right boundary of the beam moves from x B (0) = 0.6 to x B (t end ) = 0.65;

• u * 1 = 25m/s for test 6, with a Mach number around 5.5•10 -3 on the copper side of the interface and of order 1 in the gas. The right boundary of the beam moves from x B (0) = 1.75 to x B (t end ) = 1.756, namely it crosses 3 cells on the chosen grid.

Copper is deformed in the normal direction (elongation), with a rarefaction wave moving to the left. The gas is compressed and a shock moving to the right occurs. In Fig. 5, the correct solution is recovered for both simulations. The moving boundary is accurately modeled and corresponds to the discontinuity on the density profile. As expected, the centered scheme is oscillatory on the acoustic waves, especially on the shocks. For both simulations we use ∆x =2•10 -3 and we impose ν mat =0.2. This gives a time step ∆t = 8•10 -7 for test 5 (which is equivalent to impose ν ac = 0.9) and ∆t = 1.6•10 -5 , for test 6. In this latter case, the rarefaction wave is under-resolved because it is consistently faster than the wall, whose velocity dictates the large time step. The mass conservation errors of the all-speed scheme for the two copper beams are reported in Table 6. The mass conservation is guaranteed, with a convergence of roughly order one.

Multi-material interfaces

We validate the multi-material all-speed scheme introduced in Sec. 4.1.2 by presenting fluid/fluid, solid/solid and solid/fluid Sod tube problems. The tests are conducted for different regimes: the Mach numbers of the multi-material interface are reported in Tables 7 and 9. In some cases we need to distinguish between the Mach number on the left of the interface, which we call M -, and the Mach number on the right M + . The initial condition and the physical parameters of the test cases are described in Tables 8 and10, where L is the length of the tube and x B (0) is the initial position of the interface. For all test cases we use homogeneous Neumann boundary conditions.

Fluid/fluid

We perform numerical validations on fluid/fluid interaction problems. The first tests are conducted with a progressive refinement of the grid, in order to validate the spatial convergence of the all-speed multi-material scheme. Then we carry out a comparison with the mono-material all-speed scheme of [START_REF] Abbate | An all-speed relaxation scheme for gases and compressible materials[END_REF], for low Mach interfaces separating chambers filled with the same material (gas-gas and water-water). The three different discretizations of the advective operator are also compared for water/gas problems.

Grid refinement

We perform two tests with different grid spacing, with the aim of showing that the multimaterial all-speed scheme is convergent in different regimes. Test 7 is a gas-gas Sod shock tube, filled by a biatomic perfect gas. This is a monomaterial test, thus the interface corresponds to the material wave. The Mach number on the material wave is around 0.9, hence we are dealing with a fully compressible regime. Test 8 is a low Mach gas-gas Sod tube, with two different adiabatic constants γ. In the tube, a small pressure ratio and a small velocity on the right are imposed. Both gases expand in the two directions with two small rarefaction waves. The interface between the two gases moves very slowly if compared to the acoustic waves (the Mach number on the interface is M 6•10 -3 ). After 0.25s, the interface has only moved from x B (0)=0.5 to x B (t end )=0.5012, namely it has crossed 1 cell for a grid spacing ∆x=10 -3 . Fig. 6 shows that the scheme is convergent to the exact solution and it is stable. Test 7 suffers from the overheating effect, due to the fact that the interface (material wave) is treated as a moving wall. This effect slowly decreases when refining the grid.

Comparison with mono-material schemes

We compare the results of the multi-material scheme (in its upwind and all-speed versions) with the mono-material all-speed scheme of [START_REF] Abbate | An all-speed relaxation scheme for gases and compressible materials[END_REF]. Of course, this can be done only if the interface separates chambers filled with the same material, i.e. the interface corresponds to a material wave. We focus on low Mach material waves, propagating at slow velocity.

In Fig. 7 we plot a zoom on the contact wave for tests 8.1 and 9. These tests are low Mach tubes filled with perfect gas and water respectively. For test 8.1 the material wave is initially in x B (0) =200 Comparison of the multi-material schemes (upwind and all-speed) with the "mono-material" all-speed scheme and a standard explicit-upwind scheme.

spacing ∆x = 10 -1 it crosses 7 cells. For test 9, after 0.095 seconds the material wave is in x B (t end ) = 200.76, namely for a grid spacing ∆x = 10 -1 it crosses 8 cells. The different curves are obtained with a standard mono-material explicit-upwind scheme (yellow), with the mono-material all-speed scheme (green), and with the implicit multi-material schemes (blue and red). As already discussed in [START_REF] Abbate | An all-speed relaxation scheme for gases and compressible materials[END_REF], the mono-material all-speed scheme is superior in accurately capturing material waves in the low Mach regime with respect to a standard explicit-upwind scheme. With the multi-material model (no matter how the spatial derivatives are discretized) the material wave is kept sharper with respect to the mono-material scheme, since it is treated as an interface. This is due to the fact that we are introducing a "material wave reconstruction", whereas a mono-material scheme produces a "material wave approximation". A small overheating effect is present.

Comparison of the spatial discretizations

We compare the results obtained using the implicit relaxation scheme with the three discretizations of the advective part: the upwind scheme (3.10), the centered scheme (3.9) and the all-speed scheme (3.11). Each one of the three discretizations is modified to model multi-material interfaces as described in Sec. [START_REF] Fedkiw | A non-oscillatory eulerian approach to interfaces in multimaterial flows (the ghost fluid method)[END_REF] The results of tests 10 and 11 are shown in Fig. 8. We can observe that for both water/air shock tubes the interface is kept sharp and its correct position and velocity are recovered. For test 10, the interface moves from x B (0) = 0.7 to x B (t end ) = 0.712. It means that for a grid spacing ∆x = 2•10 -3 it overcomes 12 cell centers. We are dealing with a "multi-regime" test, with M -0.035 in water and M + 0.13 in the gas. The water is expanding with a rarefaction, which moves to the left very fast if compared to the velocity of the interface. A shock occurs in the gas due to compression. Test 11 deals with a low Mach regime, with M -6.5•10 -3 and M + 0.027. The interface moves very slowly, hence we consider a long integration in time. For this simulation we enforce ν mat = 0.2, giving a time step ∆t = 8•10 -4 on a grid spacing ∆x = 10 -1 . The interface moves from x B (0) = 200 to x B (t end ) = 201.103, namely it crosses 11 cells. The two acoustic waves are smoothed due to the use of a material constraint on the time step and the velocity profiles are continuous at the interface, as expected.

Solid/solid and solid/fluid

As described in [START_REF] Abbate | An all-speed relaxation scheme for gases and compressible materials[END_REF], for solids an elastic Mach number M χ = ρu 2 1 /(2χ) can be introduced. We validate the scheme on a copper/copper interface with shear and then on a copper/gas interface (see Tables 9 and10 for regimes and initial conditions). Due to copper intrinsic properties O (χ)=O (p ∞ ), we have that M χ M. This means that the speeds of shear and of longitudinal waves are of the same order of magnitude. In test 12, an interface separates a copper high pressure chamber at rest and the same material on the right at lower pressure, where a tangential velocity is imposed. With these initial state, the Mach number at the interface is around 0.14 and we expect five waves to appear. The fastest waves are relative to the normal stress, the middle one is the interface, and the two intermediate waves are relative to the tangential stress. Fig. 9 shows that the two components of the velocity field and of the stress tensor are continuous at the interface, whereas the density and the pressure are discontinuous. The interface moves from x B (0) = 0.75 to x B (t end ) = 0.794, namely for a grid spacing ∆x = 2•10 -3 it crosses 22 cells. These results are obtained for a grid spacing ∆x = 2•10 -3 and for a time step ∆t = 5•10 -7 , given by a material constraint ν mat = 0.2.

Test L x B (0) t end ρ 1 ρ 2 (u 1 ) 1 (u 1 ) 2 (u 2 ) 1 (u 2 ) 2 p 1 p 2 (m) (m) (s)
Test 13 is a shock tube containing copper at high pressure and air at atmospheric pressure. This test case is stiff because at the initial time, the pressure and density ratios are very large, since the copper is compressed. The results for an increasing number of grid points are plotted in Fig. 10, where numerical convergence is observed. The normal velocity and the normal stress are continuous at the interface. The pressure and the density are discontinuous and we can see a shock wave transmitted in air. For these results, a Newton sub-iteration to deal with the linearizations on the interface conditions is needed. As explained in Sec. 4.1.2, this is due to highly non-linear wave pattern, with

Conclusions and future developments

A completely implicit numerical method for moving boundaries and multi-material interfaces is proposed. The main purpose of the present paper is the simulation of weakly compressible and low Mach multi-material flows. To our knowledge this is the first fully implicit treatment of moving interfaces. Moving walls are accurately solved and multimaterial interfaces are kept sharp, accurately capturing their position and velocity. The scheme is able to deal with fluid/fluid, fluid/solid and solid/solid interfaces, as illustrated by the numerical tests. The multi-material model is validated on compressible and weakly compressible flows, adopting the all-speed scheme proposed in [START_REF] Abbate | An all-speed relaxation scheme for gases and compressible materials[END_REF] away from boundaries.

Several tests, including gas/gas, gas/liquid, solid/solid and solid/gas interfaces, show the excellent resolution of the multi-material interface. It is also observed that numerical oscillations almost never occur. All tests are performed by enforcing a CFL condition that depends only on the material velocity and not on the acoustic waves speed.

In the future, the scheme will be extended to solve two and three dimensional multimaterial problems. The interface will be described with a level set function, coherently with the adopted fully Eulerian approach.

Figure 1 :

 1 Figure 1: Wall/interface position in the discretized domain and interface values at the numerical interface.

Figure 2 :

 2 Figure 2: Comparison of the three spatial discretizations for two gas pistons. Test 1 : ∆x = 2•10 -3 and ∆t = 8•10 -4 . Test 2 : ∆x = 2•10 -3 and ∆t = 4•10 -2 given by ν mat = 0.2.

  -3 1.6824•10 -4

Figure 3 :

 3 Figure 3: Test 3 : Low Mach gas piston with contact discontinuity (∆x=2•10 -3 , ∆t=5•10 -2 given by ν mat =0.2).

9 Figure 4 :

 94 Figure 4: Test 4 : density profiles obtained with the upwind and all-speed spatial discretizations (∆x = 2•10 -3 and ∆t = 4.9•10 -4 ).

6 Figure 5 :

 65 Figure 5: Density profiles for copper-gas beam elongations: tests 5 and 6 (∆x = 2•10 -3 and ∆t = 8•10 -7 for test 5 and ∆t = 1.6•10 -5 for test 6).

Figure 7 :

 7 Figure 7: Test 8.1 in panel (a) and test 9 in panel (b): low Mach gas-gas and water-water tubes (zoom on the material wave in the density profile).Comparison of the multi-material schemes (upwind and all-speed) with the "mono-material" all-speed scheme and a standard explicit-upwind scheme.

Figure 8 :

 8 Figure 8: Tests 10 and 11 : water-gas tubes, results obtained with the three spatial discretizations modified with the multi-material interface conditions.

Table 10 :

 10 (Kg/m 3 ) (Kg/m 3 ) (m/s) (m/s) (m/s) (m/s) (Pa) Parameters for the solid/solid and solid/fluid test cases: initial state (material 1 on the left and material 2 on the right) and initial position of the interface x B (0).

Figure 9 : 1 .

 91 Figure 9: Test 12 : results obtained with the three spatial discretizations modified with the multi-material conditions.

Figure 10 :

 10 Figure10: Test 13 : copper-gas tube, results obtained with the multi-material all-speed scheme for an increasing number of points.

  .[START_REF] Fedkiw | A non-oscillatory eulerian approach to interfaces in multimaterial flows (the ghost fluid method)[END_REF] The first term in(2.3) is the general law of polytropic gases and the p ∞ term describes the the intermolecular forces that are present in liquids and solids. The last term is a standard neohookean law for elastic materials, where χ represents the shear elastic modulus. With a specific choices of the coefficients, we obtain standard models of different materials (see some examples in Table1

  ).

	Material	γ	p ∞	χ
			[Pa]	[Pa]
	Perfect gas	1.4	0	0
	Stiffened gas (water) 4.4	6.8•10 8	0
	Elastic solid (copper) 4.22 3.42•10 10 5•10 10

Table 1 :

 1 Examples of typical parameters for specific materials.

Table 2 :

 2 Parameters for the gas pistons: initial state (subscripts: 1 for left chamber and 2 right chamber), initial position of the wall x B (0), of the initial discontinuity x disc (0) and imposed velocity u * 1 . A biatomic gas with γ = 1.4 fills the entire domain.

		2	1.75	-	0.9	1	-	0	-	1	-	0.01
	3	2.2	2	0.5	0.6	1	1	0	0.008 0.4	0.399	0.008
	4	1	0.9	0.5	0.2	4/3	1	35/99 0	1.5	1.499	0.5•πcos(10πt)

Table 3 :

 3 Mass conservation errors for tests 1 and 2 at final time.

		(a) Test 1			(b) Test 2
	∆x	Numerical error Rate	∆x	Numerical error Rate
	8•10 -3 7.601•10 -3	-	8•10 -3 3.016•10 -3	-
	4•10 -3 4.937•10 -3	0.62	4•10 -3 2.033•10 -3	0.57
	2•10 -3 3.115•10 -3	0.66	2•10 -3 1.301•10 -3	0.64
	1•10 -3 1.892•10 -3	0.72	1•10 -3 8.309•10 -4	0.65

Table 4 :

 4 L 1 -norm errors of the density ρ and convergence rates for tests 1 and 2 at final time.

Table 6 :

 6 Mass conservation errors for tests 5 and 6 at final time.

Table 7 :

 7 Parameters for the fluid/fluid test cases: materials and regime on the interface.

	Test Materials	Regime		γ 1	γ 2	p ∞,1	p ∞,2
								(Pa)	(Pa)
	7	gas-gas	M 0.9		1.4 1.4 0	0
	8	gas-gas	M 6•10 -3		1.4 1.6 0	0
	8.1	gas-gas	M 6•10 -3		1.4 1.4 0	0
	9	water-water M 2.5•10 -3	4.4 4.4 6.8•10 8 6.8•10 8
	10	water-gas	M -0.035; M + 0.13	4.4 1.4 6.8•10 8 0
	11	water-gas	M -6.5•10 -3 ; M + 0.027 4.4 1.4 6.8•10 8 0
	Test	L (m) (m) x B (0) t end (s)	ρ 1 (Kg/m 3 ) (Kg/m 3 ) (m/s) (m/s) (Pa) ρ 2 (u 1 ) 1 (u 1 ) 2 p 1	p 2 (Pa)
	7: gas-gas	1	0.5	0.1644	1	0.125	0	0	1	0.1
	8: gas-gas	1	0.5	0.25	1	1	0	0.008	0.4	0.399
	8.1: gas-gas	400 200	150	1	1	0	0.008	0.4	0.399
	9: water-water 400 200	0.095	1000	1000	0	15	10 8	0.98•10 8
	10: water-gas	1	0.7	2.4•10 -4 1000	100	0	0	10 8	10 7
	11: water-gas	400 200	0.095	1000	50	0	25	5•10 7 2.5•10 7

Table 8 :

 8 Parameters for the fluid/fluid test cases: initial state (material 1 on the left and material 2 on the right) and initial position of the interface x B (0).

  and after 150 seconds it is in x B (t end ) =200.71, namely for a grid Figure 6: Density profiles for tests 7 and 8 : results obtained with the all-speed multi-material scheme for an increasing number of points with ν mat = 0.3.

								0.999 1.000			250 pts 500 pts 1000 pts
							density [Kg/m 3 ]	0.995 0.996 0.997 0.998
								0.994
								0.0	0.2	0.4	0.6	0.8	1.0
				(a) Test 7					(b) Test 8
	density [Kg/m 3 ]	0.9946 0.9956		exact explicit-upwind all-speed mono upwind multi all-speed multi			density [Kg/m 3 ]	995.8 995.9 996.0 996.1 996.2 996.3	exact explicit-upwind all-speed mono upwind multi all-speed multi
								995.7
		0.9936	196	198	200	202	204	995.6	196	198	200	202	204
				(a) Test 8.1 (gas-gas)

  .1.2.

		10 3									50
	density [Kg/m 3 ]	10 2	upwind centered all-speed						normal velocity [m/s]	0 10 20 30 40				upwind centered all-speed
		0.0		0.2	0.4		0.6	0.8		1.0	0.0	0.2		0.4		0.6	0.8		1.0
				(a) Density (test 10)						(b) Velocity (test 10)
		10 3									25
	density [Kg/m 3 ]	10 2	upwind centered all-speed						normal velocity [m/s]	0 5 10 15 20				upwind centered all-speed
		0	50	100	150	200	250	300	350	400	0	50	100	150	200	250	300	350	400
				(c) Density (test 11)						(d) Velocity (test 11)

Table 9 :

 9 Parameters for the solid/solid and solid/fluid test cases: materials and regime on the interface.

	Test Materials	Regime	γ 1	γ 2	p ∞,1	p ∞,2	χ 1	χ 2
						(Pa)	(Pa)	(Pa)	(Pa)
	12	copper-copper M 0.14	4.22 4.22 3.42•10 10 3.42•10 10 5•10 10 5•10 10
	13	copper-gas	M -2.5•10 -2 ; M + 0.9 4.22 1.4	3.42•10 10 0	5•10 10 0

Acknowledgments

E. Abbate acknowledges the support of the Università Italo-Francese/Université Franco-Italienne ("Bando Vinci 2016"). G. Puppo and E. Abbate were partially supported by the "National Group for Scientific Computation (GNCS-INDAM), 2018". Some experiments were carried out using the PLAFRIM experimental testbed, being developed under the Inria PlaFRIM development action with support from Bordeaux INP, LABRI and IMB and other entities: Conseil Régional dAquitaine, Université de Bordeaux and CNRS (and ANR in accordance to the programme dinvestissements dAvenir http://www.plafrim.fr/).


Jean-Claude Bermond

Graphes

Graphes orientés fortement k-connexes et graphes "k arc hamiltoniens"


Seonho Seok
email: seonho.seok@u-psud.fr

Elie Lefeuvre

Cong Wang

Jungyul Park

Electrical Characterization of Textile-Based Enzymatic Biofuel Cell for Energy Harvesting Interface Circuit

Keywords: Textile-based glucose fuel cell, energy harvesting, DC-DC converter, power management

This paper presents electrical characterization and power management circuit of textile-based enzymatic biofuel cell. Firstly, static electrical characteristics (open circuit voltage, short circuit current, and maximum power) of the enzymatic biofuel cell have been investigated by imposing a voltage across the biofuel cell and measuring corresponding output current. Output power of the biofuel cell has been characterized as a function of load resistance, while a glucose solution of 100 mM is fed from one end of the textile-based biofuel cell. Two-step approach based on supercapacitor and DC-DC converter has been proposed and tested as power management of the proposed biofuel cell.

I. INTRODUCTION

Wearable devices have gained significant interest from academia and industry owing to their immense applications such as biomedical, health and entertainment etc [START_REF] Jia | Wearable textile biofuel cells for powering electronics[END_REF][START_REF] Yi | Wearable energy sources based on 2D materials[END_REF]. To realize the internet of things (IoT) for human body data or human-machine interactions, wearable devices need to be capable of undertaking multiple complex tasks and requires long term and sustainable energy sources.

To address the above issues, the wearable power supply needs to be developed, which enables to generate electric power from body fluid for sustainable working [START_REF] Chang | Paper-based energy harvesting from salinity gradients[END_REF][START_REF] Pang | Flexible and stretchable biobatteries: monolithic integration of membrane-free microbial fuel cells in a single textile layer[END_REF]. Here, a wearable textile-based glucose fuel cell is proposed using moisture management fabrics (MMF) for improved and long-term power generation as shown in Fig. 1. MMF is widely used as basic materials for a sportswear and shows very fast sweat absorption and water evaporation. It is composed of polyester with modified cross-sectional shapes for quick water absorbing and wicking, and has the series of closely spaced channels for high evaporation. The balance between fast water evaporation and absorption enables high flow rate and continuous flow within MMF. Owing to low-cost and scalable fabrication process, it can be easily integrated into clothes, socks or underwear to utilize the human body fluid like sweat, tear or urine as fuel resource to generate energy for wearable devices. The role of power management circuit for energy harvesters is critical for delivering maximum power from energy harvesters to energy storage element, for example, battery. Basically, it is achieved by matching load resistance with the internal resistance (or impedance) of the energy harvesters in many energy harvesters. Different approaches have been reported to achieve high-efficiency energy harvesters [START_REF] Ottman | Adaptive Piezoelectric Energy Harvesting Circuit for Wireless Remote Power Supply[END_REF][START_REF] Ramadass | An Efficient Piezoelectric Energy Harvesting Interface Circuit Using a Bias-Flip Rectifier and Shared Inductor[END_REF][START_REF] Lefeuvre | Power and frequency bandwidth improvement of piezoelectric energy harvesting devices using phase-shifted SECE interface circuit[END_REF]. Direct powering of low power electronics with a biofuel cell may be undesirable due to its low output voltage and load dependent characteristics [START_REF] Zhang | A High-Efficiency DC-DC Boost Converter for a Miniaturized Microbial Fuel Cell[END_REF]. Thus, a DC-DC converter is need as an interface circuit between biofuel cell and low power electronics. Due to low energy of biofuel cell energy harvester, multiple fuel cell element interconnected, for example, microbial fuel cell (MFC), is frequently used to increase output voltage as well as power capacity [START_REF] Park | Efficient Energy Harvester for Microbial Fuel Cells using DC/DC Converters[END_REF]. In this case, two-staged DC-DC converter could be used to achieve maximal energy harvesting; the first converter harvests the energy from an MFC and charges the storage capacitor, and the second layer converter boosts the voltage to an appropriate level for the connected load [START_REF] Park | Efficient Energy Harvester for Microbial Fuel Cells using DC/DC Converters[END_REF]. In this paper, preliminary test results of a two-staged power management of a textile-based biofuel cell is presented.

II. TEXTILE-BASED ENZAYMATIC BIOFUEL CELL

The textile-based glucose fuel cell consists of MMF, carbon cloth-based cathode and anode, and conductive thread as shown Fig. 1 (a). As depicted in Fig. 1 (b), the cathode and anode were fabricated by coating carbon cloth with Prussian blue (PB) nanoparticles and glucose oxidase (GOD), respectively. Especially, the carbon cloth for the cathode is treated with 30 % PTFE to enable the air exchange within the cathode for longterm use. MMF in the middle is designed as channels with hydrophobic boundary for fuel resource absorption and transport. The conductive thread is used as a current collector to connect with the external loads.

When 100 mM glucose in PBS solution was introduced into the chamber shown in Fig. 3 (a), the glucose fuel cell started to absorb the solution through MMF layer and guided it into the reaction region where PB and GOD electrodes existed. Subsequently, the voltage signal and current signal were recorded by a computer connected to picoammeter/voltage source through GPIB card using Labview program. To analyze the I-V curve, the electrical characteristics of the glucose fuel cell including the open circuit voltage, the short circuit current and the maximum power were obtained. Fig. 2 shows static electrical characteristics of the enzymatic biofuel cell. Open circuit voltage (Voc), short circuit current (Isc) and maximum output power have been extracted from the curves; Voc = 0.3 V, Isc = 200 μA and Pmax = 12 μW/cm 2 for the best case, Voc = 0.24 V, Isc = 60 μA and Pmax = 3 μW/cm 2 for the worst case. First of all, output power of the biofuel cell has been characterized as a function of load resistance, while a glucose solution of 100 mM is fed from one end of the textile-based biofuel cell as depicted in Fig. 3. The output power has been measured up to 30 minutes at a step of 15 minutes and the optimal load resistance value is found around 100 kΩ. Given with the optimal load resistance, maximum output power is about 0.3 μW/cm 2 . Note that the output power has increased until 30 minutes after starting the measurement and the electrode (black color in the figure) size is 1 cm 2 .

(a)

(b) Fig. 3 Electrical characterization of the biofuel cell As the harvested power is not large, another sample with lager electrode has been tested, while a wind of 0.8 m/s is applied onto the biofuel cell. The wind plays a role of accelerating of liquid evaporation and thus preventing the channel from being saturated by the glucose solution. This 0.8 m/s wind speed corresponds to slow walking speed of human (the average is about 1.4 m/s [START_REF] Browning | Effects of obesity and sex on the energetic cost and preferred speed of walking[END_REF]). In the same way as the previous one, the output power has been measured as shown in Fig. 4. Interestingly, the biofuel cell with larger electrode shows higher power output than the previous one, 1.4 µW/cm 2 with an optimal resistance of 2 kΩ. After that, dynamic electrical characteristic of the biofuel cell is also characterized as presented in Fig. 5. The output voltage needs stabilization time of about 100 seconds when the load resistance is changed from 200 kΩ to 100 kΩ. This dynamic response will be useful to determine capacitance value of electrical circuit model of the biofuel cell. Due to small output power of the biofuel cell, a supercapacitor has been used as an intermediate element between the biofuel cell and DC-DC converter. The approach converting the energy generated by the biofuel cell is that it is first stored in a supercapacitor and then it is fed into DC-DC converter for power management. Hence, a supercapacitor of 220 mF is directly connected with the biofuel cell to find charging characteristic as a function of charging time as shown in Fig. 6. Charging time to reach Vsupercap = 100 mV is 40 minutes after which DC-DC conversion can be started. IV. CONCLUSION Textile-based enzymatic biofuel cell has been electrically characterized and combined with power management circuit as a wearable energy source. The biofuel cell shows electrode-size dependent internal impedance, which is a critical parameter for efficient power transfer. A two-step power management circuit, composed of a DC-DC converter combined with an intermediate supercapacitor, has been successfully demonstrated. This circuit enables to increase the output voltage, making thus possible to power various electronic devices with the harvested energy.

Fig. 1

 1 Fig. 1 Textile-based biofuel cell

Fig. 4 Fig. 5

 45 Fig. 4 Output power vs. Rload with larger electrode

Fig. 6

 6 Fig. 6 Evolution of output voltage and energy stored at 0.22 F supercapacitor A Buck-boost DC-DC converter is then connected with the supercapacitor after it is charged up to 100 mV as shown in Fig.

7 .

 7 The DC-DC converter is used to step up the supercapacitor voltage suitable for conventional battery. Energy stored at the supercapacitor is converted through inductor L1 and diode D before it is charged into output capacitor C1 by adjusting duty cycle of clock control signal of switch. Evolution of the voltage on supercapacitor and C1 at different conditions is presented in Fig.8. Output voltage is ranged from 2.5 V to 4 V and conversion time (i.e. time needed to extract more than 90% of the energy stored in the supercapacitor) is around one minute for all of cases.

Fig. 7

 7 Fig. 7 Overall energy conversion circuit of the biofuel cell