Thomas Milcent
email: thomas.milcent@u-bordeaux.fr

Emmanuel Maitre
email: manuel.maitre@imag.fr

EULERIAN MODEL OF IMMERSED ELASTIC SURFACES WITH FULL MEMBRANE ELASTICITY. *

Keywords: level set, fluid-structure, eulerian, immersed interface, Navier-Stokes, membrane elasticity subject classifications. 74F10, 76D05, 65M06

We introduce an Eulerian model for the coupling of a fluid governed by the Navier-Stokes equations, with an immersed interface endowed with full membrane elasticity (i.e. including shear effects). We show numerical evidences of its ability to account for large displacements/shear in a relatively simple way, avoiding some drawbacks of Lagrangian representation.

1. Introduction In biomechanics of immersed membranes, the ability for a numerical model and method to handle shear efficiently is a challenging issue. Human red blood cells (RBCs) are among the simplest animal cells, since they have no nucleus nor organelle. They are made, like vesicles, of a phospholipid bilayer, plus a protein network (of spectrin), the cytoskeleton. For simple vesicles, the involved energy is a curvature energy, which is minimized with an area and enclosed volume constraint. In that case, Eulerian and Lagrangian models were developed in the last decade and proved to successfully reproduce the dynamics of such objects in flow [START_REF] Beaucourt | Steady to unsteady dynamics of a vesicle in a flow[END_REF][START_REF] Cottet | A level-set formulation of immersed boundary methods for fluid-structure interaction problems[END_REF][START_REF] Thierry Biben | Phase-field approach to three-dimensional vesicle dynamics[END_REF][START_REF] Cottet | A level set method for fluid-structure interactions with immersed surfaces[END_REF][START_REF] Thierry Biben | Three-dimensional vesicles under shear flow: Numerical study of dynamics and phase diagram[END_REF][START_REF] Doyeux | Simulation of two-fluid flows using a finite element/level set method. application to bubbles and vesicle dynamics[END_REF]. For the red blood cells the cytoskeleton provides an extra resistance to membrane shear. Capsules are usually constituted of a liquid drop protected by a thin elastic membrane. They are used for applications in cosmetic, food, and pharmaceutical industries. They are more rigid objects than vesicles, so that the curvature energy is small compared to membrane elasticity, including shear. In the modeling of vesicles or capsules, taking into account shear efficiently is of paramount interest.

On one hand, the definition of local shear variation on a surface is apparently easy in Lagrangian coordinates, and in the context of capsules, has been formalized decades ago [START_REF] Skalak | Strain energy function of red blood cell membranes[END_REF][START_REF] Barthes-Biesel | Motion of a spherical microcapsule freely suspended in a linear shear flow[END_REF][START_REF] Barthes | The time-dependent deformation of a capsule freely suspended in a linear shear flow[END_REF]. However, when it comes to practical computations, the unavoidable singularity of Lagrangian parametrization of closed surfaces brings high complications [START_REF] Walter | Coupling of finite element and boundary integral methods for a capsule in a stokes flow[END_REF][START_REF] Walter | Ellipsoidal capsules in simple shear flow: prolate versus oblate initial shapes[END_REF]. Moreover, large deformation and volume conservation in this Lagrangian setting also rise numerous difficulties, which were studied by several authors in the framework of the immersed boundary method [START_REF] Charles S Peskin | The immersed boundary method[END_REF][START_REF] Lee | An immersed interface method for incompressible navierstokes equations[END_REF]. On the other hand, interface capturing methods using level-set, while more intrinsic and dimension-independent regarding surface localization, are usually reported as unable to capture tangential motions (see however somehow connected studies on surface flows and multicomponent vesicles [START_REF] Xu | A level-set method for interfacial flows with surfactant[END_REF][START_REF] John S Lowengrub | Phase-field modeling of the dynamics of multicomponent vesicles: Spinodal decomposition, coarsening, budding, and fission[END_REF][START_REF] Knut | A diffuse-interface method for two-phase flows with soluble surfactants[END_REF]). This is due to the fact that the transport equation for the level-set function, namely ∂ t φ + u • ∇φ = 0, while recording some information on the surface area variations (enabling a complex fluid formulation of fluid-interface coupling [START_REF] Cottet | A level-set formulation of immersed boundary methods for fluid-structure interaction problems[END_REF][START_REF] Cottet | A level set method for fluid-structure interactions with immersed surfaces[END_REF][START_REF] Beale | Locally corrected semi-lagrangian methods for stokes flow with moving elastic interfaces[END_REF]) just ignore any tangential component of u, since ∇φ is normal to the interface.

However, full 3D Eulerian elasticity and its coupling with fluid mechanics has already successfully been studied and implemented numerically by several teams both in the incompressible [START_REF] Cottet | Eulerian formulation and level set models for incompressible fluid-structure interaction[END_REF][START_REF] Maitre | Applications of level set methods in computational biophysics[END_REF][START_REF] Zhang | An eulerian numerical method for fluidsolid interaction[END_REF][START_REF] Richter | Finite elements for fluid-structure interaction in ale and fully eulerian coordinates[END_REF][START_REF] Richter | A fully eulerian formulation for fluid-structure-interaction problems[END_REF][START_REF] Valkov | Eulerian method for fluid-structure interaction and submerged solid-solid contact problems[END_REF] and compressible [START_REF] Kamrin | Reference map technique for finitestrain elasticity and fluid-solid interaction[END_REF][START_REF] Gorsse | A simple cartesian scheme for compressible multimaterials[END_REF] cases. The method relies on backward characteristics (also called reference map by some authors) which map back material points to their initial positions. This concept was also introduced early for point correspondence in level-set applied to image processing [START_REF] Pons | Maintaining the point correspondence in the level set framework[END_REF] and more recently to improve numerical accuracy of advection equations (e.g. level-set motions) [START_REF] Kohno | A new method for the level set equation using a hierarchical-gradient truncation and remapping technique[END_REF][START_REF] Mercier | The characteristic mapping method for the linear advection of arbitrary sets[END_REF].

This article is a first attempt, to the best of our knowledge, to write a full Eulerian model of an immersed interface endowed with full membrane elasticity (i.e. including shear effects) and to provide numerical evidences of its ability to account for large displacements/shear in a relatively simple way.

In section 2, we recall the basic definitions of interface capturing via level-set and backward characteristics. In section 3 we provide simple examples illustrating the possible drawbacks of a Lagrangian approach to record shear variations. In section 4, we build two invariants that will be used to record the full membrane elasticity. The first one records area change (in both compressible and incompressible settings), while the second measures shear. In section 5, we consider the fluid-structure coupling problem: an immersed membrane into an incompressible fluid. The whole model is reduced, following [START_REF] Cottet | Eulerian formulation and level set models for incompressible fluid-structure interaction[END_REF] to a Navier-Stokes equation with source terms expressed as a combination of partial derivatives of the backward characteristics, coupled with the vectorial transport equation giving these characteristics. We provide numerical illustrations of the ability of the numerical implementation of our model to simulate the dynamic relaxation of an initially sheared sphere. As the notion of shear is best understood through simple tests, we provide in the Annex A analytical examples of velocity fields under which the deformed continuous medium actually experiences a shear, and show that our invariants behave as expected. The Annex B deals with some technical lemmas.

Preliminaries

Forward and backward characteristics

Let Ω 0 ⊂ R 3 be the reference configuration of a continuous medium and assume that this medium is deformed by a smooth map X : R 3 × R + -→ R 3 (the forward characteristics) to Ω t = X(Ω 0 ,t). A velocity field u : R 3 × R + -→ R 3 is naturally associated with X:

∂ t X(ξ,t) = u(X(ξ,t),t), X(ξ,0) = ξ, ξ ∈ Ω 0 . (2.1)
We introduce the backward characteristics Y : R 3 × R + -→ R 3 by the formula Y (X(ξ,t),t) = ξ (see Fig 2.1). The physical interpretation of Y (x,t) is the position at time 0 of a material particle lying in x at time t and moving at speed u. The derivative of this relation with respect to t and ξ in turn gives with (2.1)

∂ t Y + (u • ∇ x )Y = 0 Y (x,0) = x, x ∈ Ω t , (2.2)
and

[∇ ξ X(ξ,t)] = [∇ x Y (x,t)] -1 for x = X(ξ,t). (2.3)
In the following, we use the notation ∇ ξ for the gradient with respect to ξ and we denote by ∇ the gradient with respect to x.

Eulerian representation of surfaces

We consider a surface Γ t captured by a level set function φ :

R 3 × R + -→ R Γ t = {x ∈ R 3 / φ(x,t) = 0}.
The initial position of the interface Γ 0 is associated to a given function φ 0 (thus we set φ(•,0) = φ 0 ). We assume that the interface is advected by the velocity field u, so that φ verifies:

∂ t φ + u • ∇φ = 0.
(2.4)

We have with (2.2) that the solution of (2.4) is given by φ(x,t) = φ 0 (Y (x,t)).

(2.5)

Let the normal in the reference configuration be denoted by n 0 (ξ) for ξ ∈ Γ 0 , and the normal in the deformed configuration by n(x,t) for x ∈ Γ t (See Fig 2 .1). These normals are expressed in terms of normalized gradient of the associated level set:

n(x,t) = ∇φ(x,t) |∇φ(x,t)| , n 0 (ξ) = ∇ ξ φ 0 (ξ) |∇ ξ φ 0 (ξ)| . (2.6)
The relation (2.5) gives ∇φ(x,t) = [∇Y (x,t)] T ∇φ 0 (Y (x,t)), therefore with the definitions (2.6),

n 0 (Y (x,t)) = [∇Y (x,t)] -T n(x,t) |[∇Y (x,t)] -T n(x,t)| . (2.7) With (2.
3), we get the Lagrangian equivalent of (2.7):

n(X(ξ,t),t) = [∇ ξ X(ξ,t)] -T n 0 (ξ) |[∇ ξ X(ξ,t)] -T n 0 (ξ)| . (2.8) Y (x,t) = ξ x = X(ξ,t) n 0 (ξ) n(x,t) Initial configuration Ω 0 Deformed configuration Ω t Γ 0 Γ t

Fig. 2.1: Initial and deformed configurations

3. Can we measure surface shear variation with a Lagrangian description of the surface? We present in this section two intuitive methods to measure shear variation with a parametrization of the surface. The first one is based on the angle between two normalized vectors and the second one uses the invariants of the metric tensor. We will show that even in the simplest case of a plane deformed in a shear velocity field, the two methods fail because the shear variation depends on the choice of initial parametrization. For general closed surfaces more difficulties arise from the singularity of Lagrangian parametrizations. The first method can not be applied because it is impossible to construct the normalized vectors (Hairy Ball theorem). We will show that for a sphere deformed in a shear velocity field the invariant of the second method is discontinuous at the poles. In any case the definition of shear depends completely on the choice of the parametrization and hence there is no canonical reference state. The limitations of these methods advocate for an immersion of the surface in R 3 using a Eulerian description to capture shear.

Two Lagrangian methods to capture shear

Let Γ 0 be the initial surface described by the parametrization Ψ 0 : U -→ R 3 where U is an open subset of R 2 . A parametrization of the deformed surface Γ t is given by Ψ : U × R + -→ R 3 where Ψ(θ 1 ,θ 2 ,t) = X(Ψ 0 (θ 1 ,θ 2 ),t) and the deformed vectors associated to the parametrization are given by a i (t) = ∂ θi Ψ(θ 1 ,θ 2 ,t). The first intuitive method to capture shear is to consider the angle between the two normalized deformed vectors

angle(t) = arccos a 1 (t) |a 1 (t)| • a 2 (t) |a 2 (t)| .
The second one is to consider the metric tensor

M(t) := a 1 (t) • a 1 (t) a 1 (t) • a 2 (t) a 1 (t) • a 2 (t) a 2 (t) • a 2 (t)
which is classically used to compute area, angles, and length of curves on the surface.

The trace of M(t) captures shear deformation but contains also information on local area. Dividing the trace by the square root of the determinant keeps only shear effects, as will be illustrated below. We thus introduce the following invariants

Z 1 (t) = det(M(t)), Z 2 (t) = Tr(M(t)) 2 det(M(t)) , (3.1)
where Z 1 (t) correspond to the local area, and we will show that Z 2 (t) is a good Lagrangian candidate to capture the local shear. In the following we are interested in the local variation of area and shear between time 0 and t. Therefore we consider the quantity angle(t)/angle(0) and the invariants Z i (t)/Z i (0). On the following examples we consider the shear velocity field (see Fig 7 .10) with the associated forward characteristics given by:

u(x,y,z,t) =   0 x 0   , X(ξ 1 ,ξ 2 ,ξ 3 ,t) =   ξ 1 ξ 2 + tξ 1 ξ 3   .
3.2. Case of the plane: orthogonal parametrization Let Γ 0 be the plane {(x,y,z) ∈ R 3 /z = 0} which is parametrized by Ψ 0 : R 2 -→ R 3 given by:

Ψ 0 : (θ 1 ,θ 2 ) →   cos(α)θ 1 -sin(α)θ 2 sin(α)θ 1 + cos(α)θ 2 0   .
where α ∈ [0,2π] is a parameter. The deformed surface Γ t is still geometrically Γ 0 in this case and the deformed vectors are given by

a 1 (t) =   cos(α) sin(α) + tcos(α) 0   , a 2 (t) =   -sin(α) cos(α) -tsin(α) 0   .
At t = 0 these vectors are orthonormal and have an α angle with the canonical basis (see Fig 3.1). We can show that for this example we have angle(t) angle(0) = 2 π arccos t(2cos(2α) -tsin(2α))

2 (1 + tsin(2α) + t 2 cos(α) 2 )(1 -tsin(2α) + t 2 sin(α) 2 ) . (3.2)
This quantity depends on the angle parameter α (even if all pair of initial vectors are orthonormal). Hence this first method is not appropriate to measure the local shear variation. For the second method, the invariants (3.1) in the case of the plane are given by

Z 1 (t) Z 1 (0) = 1, Z 2 (t) Z 2 (0) = 1 + t 2 2 . (3.3)
Hence there is no area variation, the shear is constant in space and increase in time which correspond intuitively to what we can expect of a shear deformation. The second invariant is independent of α so this second method seems appropriate to capture shear.

• α a 1 (0)

a 2 (0) • β a 1 (0) a 2 (0)
Fig. 3.1: Orthogonal (left) and non orthogonal (right) parametrization of the plane 3.3. Case of the plane: non orthogonal parametrization Let Γ 0 be the same plane {(x,y,z) ∈ R 3 /z = 0} described with another parametrization Ψ 0 : R 2 -→ R 3 given by:

Ψ 0 : (θ 1 ,θ 2 ) →   θ 1 + cos(β)θ 2 sin(β)θ 2 0   .
where β ∈ [0,2π] is a parameter. The deformed vectors are given by

a 1 (t) =   1 t 0   , a 2 (t) =   cos(β) sin(β) + tcos(β) 0   .
At t = 0 these vectors are not orthogonal and have an β angle (see Fig 3.1). For the second method, the invariants (3.1) are given by

Z 1 (t) Z 1 (0) = 1, Z 2 (t) Z 2 (0) = 1 2 1 + t 2 + cos(β) 2 + (tcos(β) + sin(β)) 2 (3.4)
There is no area variation and the shear depends on the parameter β. Therefore the second method fails to capture local shear variation because it depends on the choice of the parametrization.

3.4. Case of the sphere: orthogonal parametrization Let Γ 0 be the sphere

{(x,y,z) ∈ R 3 /x 2 + y 2 + z 2 = 1} which is parametrized by Ψ 0 : [0,π] × [0,2π] -→ R 3
given by:

Ψ 0 : (θ 1 ,θ 2 ) →   sin(θ 1 )cos(θ 2 ) sin(θ 1 )sin(θ 2 ) cos(θ 1 )   .
The deformed vectors are given by

a 1 (t) =   cos(θ 1 )cos(θ 2 ) cos(θ 1 )sin(θ 2 ) + cos(θ 1 )cos(θ 2 )t -sin(θ 1 )   , a 2 (t) =   -sin(θ 1 )sin(θ 2 ) sin(θ 1 )cos(θ 2 ) -sin(θ 1 )sin(θ 2 )t 0   .
At t = 0 these vectors are orthogonal (see Fig 3 .2). The first method with angle between vectors will also failed in this more complex case. Indeed it is impossible to construct a continuous tangent vector field on a sphere (Hairy Ball theorem). For the second method, the invariants (3.1) in the case of the sphere are given by

Z 1 (t) Z 1 (0) = 1 -tsin(θ 1 ) 2 sin(2θ 2 ) + t 2 sin(θ 1 ) 2 sin(θ 2 ) 2 , (3.5) Z 2 (t) Z 2 (0) = 1 + sin(θ 1 ) 2 + tcos(2θ 1 )sin(2θ 2 ) + t 2 2 (1 + cos(2θ 1 )cos(2θ 2 )) (1 + sin(θ 1 ) 2 ) 1 -tsin(θ 1 ) 2 sin(2θ 2 ) + t 2 sin(θ 1 ) 2 sin(θ 2 ) 2 . (3.6)
The first invariant (3.5) is well defined and represents the local area variation. The limit at the poles (θ 1 = 0 or θ 1 = π) of the second invariant (3.6) is cos(θ 2 ) 2 + (sin(θ 2 ) + tcos(θ 2 )) 2 and depends of θ 2 . Therefore there is a discontinuity of the second invariant at the poles. We can reparametrize periodically the surface but the algorithms are more complex. Moreover the results will always depends on the choice of the parametrisation, even in zones where there is no singularity.

To tackle this problem we propose to immerse the surface in R 3 and project on it the 3D deformations. In this case we use the flat metric of R 3 instead of metric of the surface and thus avoid the singularities of parametrizations. Also the immersion in R 3 will allow to get a reference state which not depends of the parametrization. Furthermore we will introduce Z i the Eulerian equivalent of the Lagrangian invariants Z i . We will show that Z 1 is the classical and well known local area variation. We will demonstrate with a variety of analytical illustrations in the Annex A that Z 2 is a "good candidate" to record the local shear variation. The surface tensor A We want to measure the deformations on the surface Γ t . Following the Lagrangian description of [START_REF] Skalak | Strain energy function of red blood cell membranes[END_REF][START_REF] Barthes-Biesel | Motion of a spherical microcapsule freely suspended in a linear shear flow[END_REF][START_REF] Barthes | The time-dependent deformation of a capsule freely suspended in a linear shear flow[END_REF] we introduce the tensor

• a2(0) a1(0)
M (X(ξ,t),t) := [∇ ξ X(ξ,t)][I -n 0 (ξ) ⊗ n 0 (ξ)]. (4.1)
Let v(ξ) be a vector defined at the point ξ ∈ Γ 0 . This vector is first projected with

[I -n 0 (ξ) ⊗ n 0 (ξ)] to v τ (ξ) ∈ T ξ Γ 0 , the tangent plane of Γ 0 at ξ. Then the vector v τ (ξ) is deformed with X in the vector [∇ ξ X(ξ,t)]v τ (ξ) at X(ξ,t). A first property is that M v is already in T X(ξ,t) Γ t , the tangent plane of Γ t at X(ξ,t). Indeed, using (2.8) and v τ (ξ) • n 0 (ξ) = 0, we have (M (X(ξ,t),t)v(ξ)) • n(X(ξ,t),t) = ([∇ ξ X] T n(X(ξ,t),t)) • v τ (ξ)) = 0
The tensor (4.1) is written in its Eulerian form with (2.3)

M (x,t) := [∇Y (x,t)] -1 [I -n 0 (Y (x,t)) ⊗ n 0 (Y (x,t))].
The associate Cauchy-Green tensor is defined as (I -n 0 ⊗ n 0 is a projector hence involutive)

A := M M T = [∇Y ] -1 (I -n 0 (Y ) ⊗ n 0 (Y ))[∇Y ] -T .
After introducing the right Cauchy-Green tensor ). As A is real and symmetric, there exists an orthonormal basis of eigenvectors. Moreover A is positive since Ax • x = |M T x| 2 ≥ 0, so we denote its eigenvalues by 0, λ 2 1 and λ 2 2 . Therefore Tr(A) = λ 2 1 + λ 2 2 and Tr(Cof(A)) = (λ 1 λ 2 ) 2 . We introduce the following invariants:

B = [∇ ξ X][∇ ξ X] T = [∇Y ] -1 [∇Y ] -T , (4.2) using (2.7), the relation A(v ⊗ v)A T = (Av) ⊗ (Av) and |[∇Y ] -T n| 2 = (Bn) • n we get A = B - (Bn) ⊗ (Bn) (Bn) • n . ( 4
Z 1 = Tr(Cof(A)) = |λ 1 λ 2 |, (4.5)
Z 2 = Tr(A) 2 Tr(Cof(A)) = 1 2
λ 1 λ 2 + λ 2 λ 1 . (4.6)
In the reference configuration (usually at initial time in our case), we have that A(0) = I -n 0 ⊗ n 0 so Tr(A(0)) = 2 and Tr(Cof(A(0))) = 1. Under the smoothness assumption on u, there exists a time T that for t < T , Tr(A) and Tr(Cof(A)) do not vanish. Therefore the invariant Z 2 is well defined. We have that Z 1 = Z 2 = 1 at initial time and the inequalities Z 1 ≥ 0, Z 2 ≥ 1. Note that Z 1 and Z 2 are both Eulerian fields, and we will derive in 5.2 equations to compute them. We will prove in the Annex B that Z 1 is exactly equal to the local area variation of the surface (even if incompressibility is not assumed on u). We will demonstrate in Annex A that Z 2 is a "good candidate" to record the local shear variation.

Immersed membrane in a fluid

In this section we apply the above considerations to build an Eulerian fluid-structure model of an immersed elastic interface, which accounts for the full membrane elasticity through an energy involving invariants Z 1 and Z 2 . We derive Eulerian equations verified by the two invariants, and the forces induced by the energy variation during the structure motion. The final model is recast as a Navier-Stokes system with a source term, the latter depending of space derivatives of Z 1 and Z 2 , coupled with their Eulerian equations. This is a generalization to full membrane elasticity of a simple level-set formulation of the immersed boundary method that was introduced previously, and that we recall now.

Mathematical model

We consider an elastic membrane Γ t immersed into a fluid governed by e.g. incompressible Navier Stokes equations, lying in a bounded domain Ω of R 3 , for a time interval [0,T ]. The interface Γ t is captured by a level-set function φ : Ω × [0,T ] → R. Based on previous works on Eulerian description of immersed interface [START_REF] Cottet | A level-set formulation of immersed boundary methods for fluid-structure interaction problems[END_REF][START_REF] Cottet | A level set method for fluid-structure interactions with immersed surfaces[END_REF] we write the coupling between the immersed surface and the surrounding fluid under the formulation:

     ∂ t u + u • ∇u -1 Re ∆u + ∇p = f (φ) on Ω × [0,T ] divu = 0 on Ω × [0,T ] ∂ t φ + u • ∇φ = 0 on Ω × [0,T ]
where f (φ) accounts for the elastic force imparted on fluid by the immersed surface. Dependence on φ means here dependence on φ and its derivatives. We have shown that in the incompressible case, |∇φ| captures the local area variation and that the force associated to the regularized energy

E = Ω E(|∇φ|) 1 ε ζ φ ε dx (5.1)
can be expressed by the formula

f (φ) = div E (|∇φ|)|∇φ| I - ∇φ ⊗ ∇φ |∇φ| 2 1 ε ζ φ ε (5.2)
where r → E(r) is the elastic constitutive law, ε the width of the interface and ζ is a cut-off function used to spread the interface near {φ = 0}. However in order to capture the full membrane energy, the function φ is not sufficient. Indeed a velocity field u tangential to Γ t would verify u • ∇φ = 0 and therefore would not change the value of φ in its evolution equation. In order to capture some information on how points are moving on the surface, one way is to use the backward characteristics Y of the velocity field u. We propose the new model:

     ∂ t u + u • ∇u -1 Re ∆u + ∇p = F 1 (φ,Y ) + F 2 (φ,Y ) on Ω × [0,T ] divu = 0 on Ω × [0,T ] ∂ t Y + u • ∇Y = 0 on Ω × [0,T ] (5.3)
where F i is the force associated to the invariant Z i , which expressions will be derived in the next section. In most cases, equations (5.3) will be supplemented with homogeneous initial and boundary conditions on the velocity, and with initial condition on Y given by identity:

u = 0 on ∂Ω × [0,T ], u = 0 on Ω × {0}, Y (0,x) = x on Ω.
Under such initial and boundary conditions, characteristics remain for all time inside Ω, thus the initial and deformed whole continuous medium (fluid and structure) always occupies the same domain of R 3 (thus there is no need to introduce an initial and deformed configuration as in 2.1). Moreover with such initial and boundary conditions, there is no need to solve an equation on φ, since the knowledge of Y and φ 0 is sufficient. Indeed the solution to the transport equation on φ is given by (see (2.5))

φ(x,t) = φ 0 (Y (x,t)).
However if we consider an initially deformed interface (see (5.20) in the numerical tests), we could choose to set Y (0,•) to some maps not equal to identity. In that case, we need to solve the transport equation on φ

∂ t φ + u • ∇φ = 0.
(5.4)

Remark 5.1. One could wonder about the optimality of our representation, regarding the number of fields to capture full membrane deformation. Note that to capture the interface position itself, we need one field (a level-set function φ). In the general case, to record the area change we need to compute Z 1 (or |∇φ| and J, see (8.1) in the Annex B), and therefore an extra field to get the shear. Thus as we solve for Y , which has three components, and obtain φ, Z 1 and Z 2 from it, we are optimal. In the incompressible case, area change is recorded by |∇φ| alone. We were not able, for the moment, to express shear by only introducing another field that we could compute without involving Y .

Eulerian equations on Z i

In order to build an Eulerian model we need to find evolution equations associated with Eulerian invariants. Proposition 5.1. Under the smoothness assumption made on u, the invariants verify:

∂ t Z 1 + u • ∇Z 1 = Z 1 [∇u] : C 1 , C 1 = I -n ⊗ n, (5.5)
∂ t Z 2 + u • ∇Z 2 = Z 2 [∇u] : C 2 , C 2 = 2A Tr(A) -(I -n ⊗ n). (5.6)
Proof. We first focus on the equations related to A. Taking the gradient of (2.2), we get

∂ t [∇Y ] + u • ∇([∇Y ]) = -[∇Y ][∇u]
and its inverse verifies:

∂ t ([∇Y ] -1 ) + u • ∇([∇Y ] -1 ) = [∇u][∇Y ] -1 .
(

The equation (4.2) on B is then given by

∂ t B + u • ∇B = [∇u]B + B[∇u] T . (5.8)
Using (2.2) we get that

∂ t (I -n 0 (Y ) ⊗ n 0 (Y )) + u • ∇(I -n 0 (Y ) ⊗ n 0 (Y )) = 0.
This relation together with (5.7) gives

∂ t M + u • ∇M = [∇u]M, ∂ t M T + u • ∇(M T ) = M T [∇u] T .
(5.9)

Using the relations (5.9) we get

∂ t A + u • ∇A = [∇u]A + A[∇u] T (5.10)
where the initial condition is given by A(0) = I -n 0 ⊗ n 0 . Note that B verifies the same equation (5.8) with a different initial condition B(0) = I. Following (5.10), we get

∂ t Tr(A) + u • ∇Tr(A) = 2[∇u] : A (5.11)
and

∂ t Tr(Cof(A)) + u • ∇Tr(Cof(A)) = 2[ATr(A) -A 2 ] : [∇u].
(5.12) Using (5.12) we get

∂ t Z 1 + u • ∇Z 1 = Z 1 [ATr(A) -A 2 ] Tr(Cof(A)) : [∇u].
We show in Lemma 8.1 of the Annex B that ATr(A) -

A 2 = Tr(Cof(A))[I -n ⊗ n].
Using (5.11) and the equation on Z 1 we get the equation for Z 2 . We have that C i n = 0 because An = 0 (see 4.4). We will prove in proposition 8.3 of Annex B that Z 1 captures the local area variation. In Annex A we provided a lot of illustrations to explain why Z 2 is a relevant measure of shear variation.

Computation of the elastic forces F i We introduce the regularized energy

E i = Ω E i (Z i ) 1 ε ζ φ ε dx. (5.13)
Here E i is the constitutive law associated to the invariant Z i . Proposition 5.2. The time variation of E i , using the principle of virtual power

∂ t E i = - Ω F i • u dx (5.14)
corresponds to the following force:

F i = ∇ E i (Z i ) 1 ε ζ φ ε + div E i (Z i )Z i C i 1 ε ζ φ ε .
(5.15)

Proof. The derivative with respect to t gives

∂ t E i = Ω E i (Z i )(Z i ) t 1 ε ζ φ ε dx + Ω E i (Z i ) 1 ε 2 ζ φ ε φ t dx.
Using the transport equation on φ and (5.5, 5.6) we get

∂ t E i = Ω E i (Z i )(-u • ∇Z i + [∇u] : Z i C i ) 1 ε ζ φ ε dx + Ω E i (Z i ) 1 ε 2 ζ φ ε (-u • ∇φ) dx.
Integrating the second term by parts gives (the integral on ∂Ω vanishes since ζ( φ ε ) = 0 on ∂Ω)

∂ t E i = - Ω u • ∇(E i (Z i )) 1 ε ζ φ ε + div E i (Z i )Z i C i 1 ε ζ φ ε • u + E i (Z i )u • ∇ 1 ε ζ φ ε dx.
Grouping the first and last term and using (5.14) leads to the expression (5.15). In the incompressible case, the gradient can be forgotten up to a redefinition of pressure and Z 1 = |∇φ| thanks to (8.1). Therefore we found the previous result (5.2). We have that

div E i (Z i )Z i C i 1 ε ζ φ ε = div(E i (Z i )Z i C i ) 1 ε ζ φ ε + E i (Z i )Z i C i n |∇φ| ε 2 ζ φ ε .
The last term vanishes because C i n = 0. With the identity div(Av) = div(A T ) • v + A T : ∇v and the symmetry of

C i we get div(E i (Z i )Z i C i n) = 0 = div(E i (Z i )Z i C i ) • n + E i (Z i )Z i C i : [∇n].
With the definition of C 1 (5.5) and (5.6) we find the normal part of forces F i (without the pressure term)

F 1 • n = -Tr([∇n])E 1 (Z 1 )Z 1 1 ε ζ φ ε , (5.16)
and with the relation 2Z 1 Z 2 = Tr(A) (see (4.6)), we get

F 2 • n = -(2Tr(A[∇n]) -Tr(A)Tr([∇n])) E 2 (Z 2 ) 2Z 1 1 ε ζ φ ε .
(5.17)

In the previous formulas we used the relation [I -n ⊗ n] : [∇n] = Tr([∇n]) ([∇n] T n = 0 since |n| 2 = 1) which is by definition the surface mean curvature. Hence the normal force associated to the area variation (5.16) only depends on the geometry through the surface curvature. The normal force associated to the shear variation (5.17) not only depends of the geometry but also of the deformations on the surface through A.

In the general case, the normal part of the shear force does not vanish.

Numerical method

The fluid-structure equations (5.3) and (5.4) are discretized with finite difference schemes on a staggered grid (see Fig 5 .1 for a 2D configuration) with the classical projection method.

• • • • • • • • • p,φ,Y × × u 1 × × u 2 ∆x 2 ∆x 1

Fig. 5.1: Staggered grid with position of unknowns

Let ∆t be the time step and u n ,p n ,φ n ,Y n the time discretization of the variables at t n = n∆t. The semi-discretization in time is given by

Step 1:

u -u n ∆t + u n • ∇u n - 1 R e ∆u n = F 1 (φ n ,Y n ) + F 2 (φ n ,Y n ) Step 2: ∆p n+1 = div(u ) ∆t Step 3: u n+1 = u -∆t∇p n+1
Step 4:

φ n+1 -φ n ∆t + u n+1 • ∇φ n = 0 Y n+1 -Y n ∆t + u n+1 • ∇Y n = 0
In step 1 a prediction of the velocity is computed with an explicit Euler scheme in time and classical central staggered velocity-pressure schemes of order two for the convection, diffusion and the source terms. An explicit scheme in time is used for the diffusion because the Reynolds number is large enough. In step 2, the Poisson equation for the pressure is performed with the Fishpack library [START_REF] Adams | Fishpack -efficient fortran subprograms for the solution of separable elliptic partial differential equations[END_REF] (we use homogeneous Neumann boundary conditions). In step 3 the the velocity is corrected to enforce the incompressibility condition. In step 4 the transport equations are discretiezd with an explicit Euler scheme in time and a WENO5 scheme in space [START_REF] Jiang | Efficient implementation of weighted eno schemes[END_REF].

We do not perform the redistancing. Instead, we use the renormalization φ |∇φ| to measure the distance to interface. Thus,

|∇φ| 1 ε ζ φ ε is replaced by 1 ε ζ φ |∇φ|ε .
This approach was proved in [START_REF] Cottet | A level set method for fluid-structure interactions with immersed surfaces[END_REF] to be efficient from the point of view of both volume conservation and interface force calculations. For the cut-off function, we considered the following expression ζ(r) = 1 2 (1 + cos(πr)) on [-1,1] and ζ(r) = 0 elsewhere. We use in our simulations the linear elastic laws

E 1 (r) = λ(r -1) E 2 (r) = µ(r -1) (5.18)
where λ and µ are the elastic stretching and shear modulus (however the model is still strongly nonlinear due to the geometric nonlinearities and the coupling with Navier-Stokes equations).

Numerical tests

In this section we provide evidences of numerical convergence of our method on the illustrative test case of a sheared elastic sphere. The domain Q = [-1,1] 3 is discretized on a Cartesian mesh with 128 points in each direction. We choose in our simulation a Reynolds number R e = 100, an elastic stretching modulus λ = 1 and a shear modulus µ = 0.1. The parameter ε is equal to 3.5∆x in the simulations where ∆x is the grid size. We take the time step ∆t = 1.310 -3 s. We impose zero velocity for the initial and boundary conditions. The immersed initial surface is a sphere, so that φ 0 (x,y,z) = x 2 + y 2 + z 2 -0.5, (5.19) but which is not in its rest state, since material points have been moved so that Y (x,y,z,0) = (xcos(t 0 z) + y sin(t 0 z),-xsin(t 0 z) + y cos(t 0 z),z) (5.20) This corresponds to a deformation of the sphere when a 3D circular shear (see expression of Y for TC4 in Table 7.1 and Fig 7 .8) is applied on opposite directions on the north and south poles (with respect to the z axis) until t = t 0 . Here we take t 0 = π. While this initial deformation was (artificially) imposed with no area variation (actually the sheared surface is still geometrically a sphere), the area will change when it starts to relax, so the force F 1 is also involved. The motion is however driven by F 2 .

The numerical results at different time steps are presented in Fig 5 .2. For plotting purposes, we represented on the deformed surface a check board pattern which was tracked with markers, to visualize the shear relaxation. In Fig 5 .3 we plot the velocity magnitude at t = 0.5s, while the relaxation takes place. Due to the large imposed shear, the surface undergoes a complex deformation involving some small folds. These kind of ripples were also observed in [START_REF] Walter | Coupling of finite element and boundary integral methods for a capsule in a stokes flow[END_REF] in the simulation of a capsule in a simple shear flow. A nice feature of our numerical code is its ability to relax toward a stable solution without any curvature energy. Some grid effects are however present for large time (see last picture of Fig 5 .2). The use of an unstructured mesh and finite-element solver to remove this effect is presently under study. computation were performed with N points in each direction for three resolutions with N = 64,128,256. Accordingly, we chose ε = 1.75∆x,3.5∆x,7∆x where ∆x is the corresponding grid size, so that the numerical interface width remains constant. We observe numerical convergence. Likewise, Fig 5 .6 depicts the variation of vertical radius during time. As the flow is incompressible, the map Y (•,t) : Ω → Ω conserves volume, that is, det∇Y = 1 in the continuous case. However, after time and space discretization, and due to numerical errors introduce while solving for the transport equations on Y , such nonlinear combination of derivatives will not fully respect this constraint. On Fig 5 .7 and Fig 5 .8, we depict the RMS (Root Mean Square) of det∇Y -1 as a function of time, both for the whole domain, and on the interface. Thus the plotted quantities are respectively:

RM S Ω (t) = 1 |Ω| Ω |det∇Y (x,t) -1| 2 dx 1 2
,

and RM S Γt (t) = 1 |Γ t | Γt |det∇Y (x,t) -1| 2 ds 1 2
.

(5.21)

While N = 64 is clearly under-resolved, the figures show numerical convergence of those RMS as N increases.

Conclusion

In this paper we introduced a way to account for membrane shear in a completely Eulerian framework. The idea was to build two invariants from algebraic combinations of the space derivatives of the backward characteristics of the continuous medium velocity field. The first invariant accounts for changes of local area of the membrane, while the second records the local membrane shear. We showed how to compute those invariants in an Eulerian way, and produced some evidences that the second invariant indeed captures the membrane shear. The fluid-structure coupling is therefore rephrased as a Navier-Stokes equation with source term involving the two invariants, coupled to a vectorial transport equation for the backward characteristics. The formulation is therefore very versatile since it could be implemented as an addon on a pre-existing fluid solver. In this first work, we presented outputs of a 3D numerical code using finite differences and a projection method on a staggered mesh, which simulate the relaxation of a spherical membrane which has been sheared from its poles. Our model and code is able to capture ripples due to the bucking of membrane, in a stable way without any extra smoothing curvature energy. Future work will focus on the application of this method to model red blood cells and capsules in 3D flows.

7. Annex A: Analytical illustrations for Z 2 We present some analytical illustrations to show that Z 2 is intuitively a "good" measure of the local shear variation of a surface. In all the test cases, we define an initial surface Γ 0 = {φ 0 = 0} and a velocity field u which will move material points of this surface. Then we compute the backward characteristics Y and the invariants Z 1 and Z 2 to see how these Eulerian quantities could record this motion. The test cases and the results are described in Table 7.1.

Deformations where Γ t = Γ 0

In this first serie of test cases (TC1 to TC4), while material points could have moved, globally the initial and deformed surface are the same (Γ 0 = Γ t ). Moreover, the deformations are uniform in space in the sense that the invariants Z i on the surface do not depend on the spatial variables (except for TC4).

What we refer to the α and β deformations (TC1 and TC2) are 2D and the initial surface is the plane z = 0. The velocity field of each deformation is plotted in the figures below with the corresponding values of Z 1 and Z 2 . The β = -1 deformation is a rotation and as expected (slipping from hyperbolic to circular functions with √ β identified to i ∈ C), there is no area variation and shear variation (see Fig 7.1). The α = 1 deformation is a pure dilatation and as expected, there is only area variation (see Fig 7 .2). The β = 0, β = 1, α = -1 deformations correspond to different shear transformations and as expected, there is only shear variation (see Fig 7.3,7.4,7.5). Note that for the β = 0 test case we take the limit β -→ 0 the results for Y and Z i . Moreover, for the β = 0 deformation we found the same invariants as in the Lagrangian case (see (3.3)). The α = 0 deformation is an uniaxial deformation and there is area and shear variation (see Fig 7.6). It can be quite surprising at first sight, but when a surface is stretched in two directions with different magnitude (α = 1) we do have shear and therefore as expected Z 2 = ch(t(1 -α)) = 1. In the "3D circular shear" test cases the initial surfaces are a cylinder (TC3 see Fig 7 .7) and a sphere (TC4 see Fig 7.8). In each plane {z = α} the velocity is a rotation of magnitude α. In these test cases there is no area variation, but pure shear variation. For TC3, Z 2 = 1 + t 2 2 is constant on the surface (x 2 + y 2 = 1 on the cylinder). This test case is clearly a 3D generalization on a cylinder of the 2D β = 0 deformation and this is why we found the same invariants. For TC4, Z 2 = 1 + t 2 2 (1 -z 2 ) 2 on the surface and depends only on the height z which is intuitive (x 2 + y 2 + z 2 = 1 on the sphere).

TC φ 0 (x,y,z) u Y (x,y,z,t) Z 1 Z 2 1 z   x αy 0     xe -t ye -αt z   e t(1+α) ch(t(1 -α)) 2 z   βy x 0     xch(t √ β) -y √ βsh(t √ β) -x √ β sh(t √ β) + ych(t √ β) z   1 1 + (1+β) 2 2β sh 2 (t √ β) 3 x 2 + y 2 -1   -yz xz 0     xcos(tz) + ysin(tz) -xsin(tz) + ycos(tz) z   1 1 + t 2 2 (x 2 + y 2 ) 4 x 2 + y 2 + z 2 -1   -yz xz 0     xcos(tz) + ysin(tz) -xsin(tz) + ycos(tz) z   1 1 + t 2 (x 2 +y 2 ) 2 2(x 2 +y 2 +z 2 ) 5 x 2 + y 2 + z 2 -1   x y z     e -t x e -t y e -t z   e 2t 1
Z 1 = 1 Z 2 = 1 + t 2 2 (x 2 + y 2 )
Fig. 7.7: Velocity field (-yz,xz,0) of the 3D circular shear and initial shape of the cylinder (TC3) The calculations are done independently on each plane that composes the cube. We have the following results

Z 1 = 1 Z 2 = 1 + t 2 (x 2 +y 2 ) 2 2

On the faces {x

= ±1} Z 1 = 1, Z 2 = 1 (7.1)
On the faces {y -tx = ±1}

Z 1 = 1 + t 2 Z 2 = 2 + t 2 2 √ 1 + t 2 (7.2) On the faces {z = ±1} Z 1 = 1, Z 2 = 1 + t 2 2 (7.3)
On the faces corresponding to {x = ±1} there is no area and shear variation (the faces are just translated). For the faces of equations {y -tx = ±1} the faces are stretched in only one direction and there is area and shear variation as in the 2D α = 0 deformation. For the faces {z = ±1} there is only shear variation as in the 2D β = 0 deformation. For the test case TC7 the deformed shape of the sphere Γ t is given by the zero level set of φ(x,y,z,t) = x 2 + (y -tx) 2 + z 2 -1.

We have that

Z 1 = 1 + t 2 + t 2 (x 2 -z 2 ) -2txy x 2 + (y -tx) 2 + z 2 (7.4) Z 2 = 1 + t 2 2 + t 2 x 2 -2txy 2(x 2 + (y -tx) 2 + z 2 ) 1 Z 1 (7.5)
In Fig 7 .12, iso-contours of Z 1 and Z 2 are plotted on the deformed surface Γ t . The results are intuitive: the area variation reaches its maximum on points that are orthogonal to the y = tx plane while the shear variation is larger along the z axis. The deformation and the initial surface are the same as in the Lagrangian test case in section 3.3. In the formula (7.5) there is no discontinuities at the poles (x 2 + (y -tx) 2 + z 2 = 1 on the surface Γ t ) unlike in the Lagrangian formula (3.6) Fig. 7.12: Plot of iso-contours of Z 1 and Z 2 on the deformed surface at t = 0.5s (TC7) Proof. It is equivalent to prove the lemma in any basis of R 3 . We consider the orthonormal basis of R 3 B = (e 1 ,e 2 ,e 3 ) = (τ 1 ,τ 2 ,n) where (τ 1 ,τ 2 ) is an orthonormal basis of the tangent plane orthogonal to n. Let A be the matrix of the tensor A in the basis B and we denote by A ij its coefficients. We have with (4.4) that An = 0, hence we have A i3 = 0. Since A is a symmetric tensor and B is an orthonormal basis, the matrix A is a symmetric and its structure is given by: It is now easy to check that A Tr(A ) -(A ) 2 = Tr(Cof(A ))[I -n ⊗ n] and the lemma is proved.

8.2. Equivalent expression for Z 1 In this subsection we give an equivalent expression of Z 1 which links the present formulation with our previous work [START_REF] Cottet | A level set method for fluid-structure interactions with immersed surfaces[END_REF]. Proposition 8.2. Let J = det(∇ ξ X) = det(∇Y ) -1 be the volume ratio, then We have the relation det(B) = J 2 so the expression of Z 1 (8.2) reduces to (8.1). This recovers the result of [START_REF] Cottet | A level set method for fluid-structure interactions with immersed surfaces[END_REF] where the area variation captured by |∇φ| in the incompressible case (J = 1).

Z 1 = J |∇φ| |∇φ 0 (Y )| . ( 8

8.3.

Proof that Z 1 is the local area variation Proposition 8.3. Let Γ 0 a smooth surface, deformed with smooth forward characteristics in Γ t = X(Γ 0 ,t). We denote by u the associated smooth velocity field. Let a parametrization of the surface Γ t be given by (θ 1 ,θ 2 ) → Ψ(θ 1 ,θ 2 ,t) where Ψ : This result holds for all function g. Therefore (8.3) holds and Z 1 measures the local area variation (in the compressible and incompressible cases).

U × R + -→ R 3 is

Fig. 3 . 2 :

 32 Fig. 3.2: Orthogonal parametrization of the sphere

. 3 ) 4 . 2 .

 342 Invariants of A According to (4.3) we have An = 0. (4.4) Thus 0 is an eigenvalue and det(A) = 0. The other invariants are Tr(A) and Tr(Cof(A)) = 1 2 (Tr(A) 2 -Tr(A 2 )

Fig. 5 . 2 :

 52 Fig. 5.2: Numerical simulation of the relaxation of a sphere due to shear force at times t = 0, 0.5s, 1s, 1.5s, 2s, 2.5s, 3s, 3.5, 4s, 4.5s, 5s, 9s. Lagrangian markers are used for plotting a check board pattern coming back to undeformed state when shear decreases.

Fig. 5 . 4 :

 54 Fig. 5.4: Pressure along x axis at time t = 0.1 (Top) and t = 1.2 (Bottom) for N = 64,128,256.

Fig. 5 . 5 :

 55 Fig. 5.5: Pressure along z axis at time t = 0.1 (Top) and t = 1.2 (Bottom) for N = 64,128,256.

Fig. 5 . 6 :

 56 Fig. 5.6: Vertical radius up to t = 30 (Top). Zoom on t ∈ [0,5] for N = 64,128,256.(Bottom)

Fig. 5 . 7 :

 57 Fig. 5.7: Time evolution, for N = 64,128,256 of t → RM S Ω (t) defined in (5.21).

Fig. 5 . 8 :

 58 Fig. 5.8: Time evolution, for N = 64,128,256 of t → RM S Γt (t) defined in (5.21).

5 Z 1 = 1 Z 2 = 1 Fig. 7 . 1 : 5 Z 1 = 2t Z 2 = 1 Fig. 7 . 2 : 5 Z 1 = Z 2 = 1 + t 2 2 Fig. 7 . 3 :Fig. 7 . 4 :Fig. 7 . 5 : 5 Z 1 Fig. 7 . 6 :

 5112171512172512127374755176 Fig. 7.1: Velocity field (-y,x,0) of the β = -1 deformation (TC2)

(x 2 +y 2 +z 2 )Fig. 7 . 8 := 1 Fig. 7 . 9 :Fig. 7 . 10 :Fig. 7 .

 2781797107 Fig. 7.8: Velocity field (-yz,xz,0) of the 3D circular shear and initial shape of the sphere (TC4)

1 .

 1 The following identity holds:ATr(A) -A 2 = Tr(Cof(A))[I -n ⊗ n].

A

  ) 2 + (A 12 ) 2 A 11 A 12 + A 12 A 22 0 A 11 A 12 + A 12 A 22 (A 12 ) 2 + (A 22 )In the basis B we haveI -n ⊗ n =and A have the same invariants so Tr(A) = Tr(A ) = A 11 + A 22 Tr(Cof(A)) = Tr(Cof(A )) = A 11 A 22 -(A 12 ) 2

. 1 ) 2 and 2 2 =

 1222 Proof. The invariant Z 1 is defined by (4.5)Z 1 = Tr(Cof(A)) (8.2)With the definition (4.3), we getTr(A) 2 = Tr(B) -(B 2 n) • n (Bn) • n Tr(A 2 ) = Tr(B 2 ) -2 (B 3 n) • n (Bn) • n + (B 2 n) • n(Bn) • n Using the Cayley-Hamilton theorem B 3 -Tr(B)B 2 + Tr(Cof(B))B -det(B)I = 0 we get Tr(Cof(A)) = 1 (Bn) • n (B 3 n) • n -Tr(B)(B 2 n) • n + Tr(Cof(B))(Bn) • n = det(B) (Bn) • nThe gradient of (2.5) gives the relation ∇φ = [∇Y ] T ∇φ 0 (Y ) so (Bn) • n = [∇Y ] -T ∇φ |∇φ| |∇φ 0 (Y )| 2 |∇φ| 2

Table 7 .

 7

1: Table of test cases

  smooth on an open set U of R 2 . The local area variation verifies|Ψ ,1 (θ 1 ,θ 2 ,t) ∧ Ψ ,2 (θ 1 ,θ 2 ,t)| |Ψ ,1 (θ 1 ,θ 2 ,0) ∧ Ψ ,2 (θ 1 ,θ 2 ,0)| = Z 1 (Ψ(θ 1 ,θ 2 ,t),t) Z 1 (Ψ(θ 1 ,θ 2 ,0),0) = Z 1 (x,t) Z 1 (Y (x,t),0) (8.3)for x = Ψ(θ 1 ,θ 2 ,t) and with the notationΨ ,i (θ 1 ,θ 2 ,t) := ∂ θi Ψ(θ 1 ,θ 2 ,t).Proof. Let f : R 3 × R + -→ R be a smooth function. The Reynolds formula for surfaces readsΓt ∂ t f + u • ∇f + f [∇u] : [I -n ⊗ n] ds.Let g : R 3 -→ R be smooth function and let f (x,t) = g(Y (x,t)) Z1(x,t) . The expression in the previous integral becomes1 Z 1 (∂ t (g(Y )) + u • ∇(g(Y ))) -g(Y ) (Z 1 ) 2 (∂ t Z 1 + u • ∇Z 1 -Z 1 [∇u] : [I -n ⊗ n])(8.4)The first term vanishes with (2.2). The second term is also zero with(5.5). ThereforeSince Ψ(θ 1 ,θ 2 ,0) = Y (Ψ(θ 1 ,θ 2 ,t),t) the equation (8.5) becomes |Ψ ,1 (θ 1 ,θ 2 ,t) ∧ Ψ ,2 (θ 1 ,θ 2 ,t)| dθ 1 dθ 2 = U g(Ψ(θ 1 ,θ 2 ,0)) Z 1 (Ψ(θ 1 ,θ 2 ,0),0) |Ψ ,1 (θ 1 ,θ 2 ,0) ∧ Ψ ,2 (θ 1 ,θ 2 ,0)| dθ 1 dθ 2

	d dt Γt	f ds =
		d dt Γt	g(Y (x,t)) Z 1 (x,t)	ds = 0	(8.5)

U g(Ψ(θ 1 ,θ 2 ,0)) Z 1 (Ψ(θ 1 ,θ 2 ,t),t)