Matthieu Nicolas
email: matthieu.nicolas@inria.fr

Efficient renaming in Conflict-free Replicated Data Types (CRDTs)

Sequence Conflict-free Replicated Data Types (CRDTs) allow to replicate and edit, without any kind of coordination, sequences in distributed systems. To ensure convergence, existing works from the literature add metadata to each element but they do not bound its footprint, which impedes their adoption. Several approaches were proposed to address this issue but they do not fit a fully distributed setting. In this paper, we present our ongoing work on the design and validation of a fully distributed renaming mechanism, setting a bound to the metadata's footprint. Addressing this issue opens new perspectives of adoption of these CRDTs in distributed applications.

Research Statement

In order to serve an ever-growing number of users and provide an increasing volume of data, large-scale systems such as data stores or collaborative editing tools adopt a distributed and replicated architecture. However, as stated by the CAP theorem [START_REF] Brewer | Towards Robust Distributed Systems[END_REF], such systems cannot ensure both strong consistency and high availability in presence of network partitions. As a result, literature and companies increasingly adopt an optimistic replication model [START_REF] Saito | Optimistic Replication[END_REF] paired with the eventual consistency model to replicate data among nodes. This consistency model allows replicas to temporarily diverge to ensure high availability, even in case of network partition. Each node owns a copy of the data and can modify it before propagating updates to others. A conflict resolution mechanism is however required to handle updates generated concurrently by different nodes.

An approach, which gains in popularity since a few years, proposes to define Conflict-free Replicated Data Types (CRDTs) [START_REF] Shapiro | Conflict-Free Replicated Data Types[END_REF]. These data types behave as traditional ones, such as the Set or Sequence data type, but are designed for a distributed usage. Their specification ensures that concurrent updates are resolved deterministically, without requiring any kind of agreement, and that replicas eventually converge immediately after observing some set of updates, thus achieving Strong Eventual Consistency [START_REF] Shapiro | Conflict-Free Replicated Data Types[END_REF].

To achieve convergence, CRDTs proposed in the literature mostly rely on unique identifiers to reference updated elements. To generate such element identifiers, nodes often use their own identifier combined with a logical clock. Thus, regarding how node identifiers are generated, the size of element identifiers usually increases with the number of nodes. Element identifiers also have to comply to additional constraints according to the CRDT, for instance forming a dense set in the case of a Sequence data type. Consequently, the size of element identifiers also increases according to the number of updates performed on the data structure. Therefore, the size of element identifiers is usually not bounded.

Since the size of identifiers attached to each element is not bounded, the overhead of the replicated data structure increases over time. Since nodes have to store and broadcast these identifiers, the application's performances and efficiency decrease over time. This severely impedes the adoption of concerned CRDTs.

In this thesis, we aim to address this issue of the growth of identifiers for a specific category of CRDTs suffering particularly from this issue: Sequence CRDTs.

Related Work

The Sequence data type represents an ordered collection elements.

Two operations are usually defined to update the sequence: insert(index, element) adds the element at the given index whereas delete(index) removes the element at the given index.

This data type is often used in applications and, over the years, several CRDTs [START_REF] Shapiro | A comprehensive study of Convergent and Commutative Replicated Data Types[END_REF] were proposed to represent replicated sequences, which can be divided into two approaches.

The first one [START_REF] Briot | High Responsiveness for Group Editing CRDTs[END_REF] proposes to add to each element a timestamp and a reference to its predecessor. Using these data, a total order between elements is defined. However this approach requires to keep tombstones for deleted elements to handle concurrent insertions.

The second approach [START_REF] André | Supporting Adaptable Granularity of Changes for Massive-Scale Collaborative Editing[END_REF], on which we focus from this point, proposes to attach an identifier to each inserted element. These identifiers are used to achieve commutative updates by uniquely identifying each element and ordering them relatively to each other.

The downside of this approach is the increasing size of identifiers. Since identifiers are used to order elements, they have to form a dense set so that nodes are always able to insert a new element between two others. However, two identifiers of the same size can be contiguous. When inserting a new element between two such identifiers, there is no other choice than to increase the size of the generated identifier to be able to generate one respecting the intended order. This leads to a specification of identifiers which does not bound their size.

Several approaches were proposed in the literature in order to address the issue of ever-growing identifiers. Nédelec et al. [START_REF] Nédelec | LSEQ: an adaptive structure for sequences in distributed collaborative editing[END_REF] propose several strategies to generate identifiers and design a mechanism to switch between these strategies to limit the speed of the identifiers growth. Let , ia et al. [START_REF] Let | CRDTs: Consistency without concurrency control[END_REF] propose a renaming mechanism to reduce the size of currently used identifiers, relying on a consensus algorithm. To prevent a system-wide consensus, authors divide the system into two tiers: the Core, a small set of controlled and stable nodes, and the Nebula, an uncontrolled set of nodes. Both can perform updates but only the members of the Core participate in the consensus leading to a rename. A catch-up mechanism is provided for the Nebula to transform the updates performed concurrently.

Nonetheless, using consensuses is expensive and not suited to all kinds of applications. In a fully distributed application, there is no central authority to provide a set of stable nodes acting as the Core. A efficient distributed renaming mechanism is thus needed.

Proposed Approach

We propose a fully distributed renaming mechanism, allowing any node to perform a renaming without any kind of coordination with others. The goal of such mechanism is to reassign shorter identifiers to each current element of the data structure to reduce its footprint.

We define a new operation : rename. This operation generates a new, equivalent state of the current sequence but with arbitrary identifiers of minimal size. A map, called renamingMap, is also generated to link each identifier from the former sequence to the corresponding identifier in the new one. By propagating this map, other nodes are able to apply the same renaming to their own local copy. The renamingMap is only used for renaming, it can thus be garbage collected once the corresponding rename operation has been observed by every node.

However, being a distributed system, nodes can perform updates and broadcast them while another node is performing concurrently the process. Upon the reception of these concurrent updates, the node which triggered the renaming process cannot apply them to its copy. Indeed, since the identifiers have been modified, performing the updates would not insert the elements at the correct positions nor delete the correct elements. To address this issue, we designed rewriting rules for concurrent operations to a rename. By using the original identifier and the renamingMap, we are either able to find the new corresponding identifier, if it has been renamed, either able to deterministically generate a new identifier preserving the existing order, if it has been concurrently inserted. The same set of rules can be used by nodes to compute their locale state upon the reception of rename operations.

However, one rename operation can be concurrent to another one. To solve this scenario, we give priority to one operation over the others. Nodes which observed and applied one of the "losing" renaming operations have to undo its effect before applying the "winning" one. Thus we designed additional rewriting rules which reverse the effect of a previously applied rename while also generating new identifiers preserving the order for identifiers generated after the renaming. To determine which rename operation to proceed with, we define and use a total order between rename operations.

Based on this approach, we propose a new version of the Lo-gootSplit [START_REF] André | Supporting Adaptable Granularity of Changes for Massive-Scale Collaborative Editing[END_REF] algorithm implementing this renaming mechanism. LogootSplit is the state of art of Identifier-based Sequence CRDTs. It reduces the footprint of the data structure by aggregating elements into blocks. The renaming mechanism, paired to LogootSplit, allows us to reassign shorter identifiers to elements, but also to merge existing blocks into one, reducing even more the metadata used. 1 We also introduce an optimisation to reduce the bandwidth consumption of the rename operation by compressing the renam-ingMap.

Ongoing Validation

The proposed renaming mechanism has been implemented in MUTE [START_REF] Nicolas | MUTE: A Peer-to-Peer Web-based Real-time Collaborative Editor[END_REF], the real time peer-to-peer collaborative text editor developed by our research team. It is now in the process of being tested.

We are also taking steps to provide a formal specification of the algorithm and prove its correctness, either using a proof assistant either using an automatic theorem prover.

Although the renaming mechanism is bound to reduce the memory usage of the data structure over time, we need to assess that the additional computations do not degrade the overall performances. We are thus designing a benchmarker for collaborative editing tools. The goal is to replay a collaborative editing session, from its logs, using different conflicts resolution mechanisms to compare their respective performances. We intend to measure common metrics like the memory and CPU usage, bandwidth consumption but also more specific ones like the time to propagate and integrate an update and the global time needed to converge.

Conclusion and Future Work

The proposed renaming mechanism allows us to set a bound to the size of the identifiers. By bounding the size of identifiers, we open new perspectives for the adoption of concerned CRDTs in systems like distributed databases.

Nonetheless, we are still working on several aspects of the proposed renaming mechanism. The strategy used to trigger a renaming has yet to be designed and evaluated. It should keep the size of identifiers low while preventing many rename operations to be generated concurrently to achieve good performances. We could also rework the strategy used to pick a "winning" operation from a set of concurrent rename operations. The current one, defining and using a total order between these operations, can lead to cases in which nodes have to undo many rename operations because of only one concurrent operation. This could provoke performance issues. New strategies could be proposed to reduce the likelihood of these scenarios and improve the overall performances of the system.

Middleware'18, Rennes, France 2018. 978-x-xxxx-xxxx-x/YY/MM. . . $15.00 DOI: 10.1145/nnnnnnn.nnnnnnn

For further details on the proposed approach, please read https://bit.ly/2yHkdmA


M Otténio

M Destrade

R W Ogden

Incremental magnetoelastic deformations, with application to surface instability

In this paper the equations governing the deformations of infinitesimal (incremental) disturbances superimposed on finite static deformation fields involving magnetic and elastic interactions are presented. The coupling between the equations of mechanical equilibrium and Maxwell's equations complicates the incremental formulation and particular attention is therefore paid to the derivation of the incremental equations, of the tensors of magnetoelastic moduli and of the incremental boundary conditions at a magnetoelastic/vacuum interface.

The problem of surface stability for a solid half-space under plane strain with a magnetic field normal to its surface is used to illustrate the general results. The analysis involved leads to the simultaneous resolution of a bicubic and vanishing of a 7 × 7 determinant. In order to provide specific demonstration of the effect of the magnetic field, the material model is specialized to that of a "magnetoelastic Mooney-Rivlin solid". Depending on the magnitudes of the magnetic field and the magnetoelastic coupling parameters, this shows that the half-space may become either more stable or less stable than in the absence of a magnetic field.

Introduction

One of the main reasons for industrial interest in rubber-like materials resides in their ability to dampen vibrations and to absorb shocks. This paper is concerned with an extension of the nonlinear elasticity theory adopted for describing the properties of these materials to incorporate nonlinear magnetoelastic effects so as to embrace a class of solids referred to as magnetosensitive (MS) elastomers. These "smart" elastomers typically consist of an elastomeric matrix (rubber, silicon, for example) with a distribution of ferrous particles (with a diameter of the order of 1-5 micrometers) within their bulk. They are sensitive to magnetic fields in that they can deform significantly under the action of magnetic fields alone without mechanical loading, a phenomenon known as magnetostriction. As a result, their mechanical damping abilities can be controlled by applying suitable magnetic fields. This coupling between elasticity and magnetism was probably first observed by Joule in 1847 when he noticed that a sample of iron changed its length when magnetized.

In general, the physical properties of magnetoelastic materials depend on factors such as the choice of magnetizable particles, their volume fraction within the bulk, the choice of the matrix material, the chemical processes of curing, etc.; see [START_REF] Bellan | Field dependence of viscoelastic properties of MR elastomers[END_REF] for details, and also [START_REF] Rigbi | The response of an elastomer filled with soft ferrite to mechanical and magnetic influences[END_REF] for an experimental study on a magneto-sensitive elastomer.

The coupling between magnetism and nonlinear elasticity has generated much interest over the last 50 years or so, as illustrated by the works of Truesdell and Toupin [START_REF] Truesdell | The classical field theories[END_REF], Brown [START_REF] Brown | Magnetoelastic Interactions[END_REF], Yu and Tang [START_REF] Yu | Magneto-elastic waves in initially stressed conductors[END_REF], Maugin [START_REF] Maugin | Continuum Mechanics of Electromagnetic Solids[END_REF], Eringen and Maugin [START_REF] Eringen | Electrodynamics of Continua[END_REF], Kovetz [START_REF] Kovetz | Electromagnetic Theory[END_REF], and others. The corresponding engineering applications are more recent (see Jolly et al. [START_REF] Jolly | A model of the behaviour of magnetorheological materials[END_REF], or Dapino [START_REF] Dapino | On magnetostrictive materials and their use in adaptive structures[END_REF], for instance) and have generated renewed impetus in theoretical modelling (see, for example, Dorfmann and Brigadnov [START_REF] Dorfmann | Constitutive modelling of magnetosensitive Cauchy-elastic solids[END_REF]; Dorfmann and Ogden [START_REF] Dorfmann | Nonlinear magnetoelastic deformations[END_REF]); Kankanala and Triantafyllidis [START_REF] Kankanala | On finitely strained magnetorheological elastomers[END_REF]. Here, we derive the (linearized) equations governing incremental effects in a magnetoelastic solid subject to finite deformation in the presence of a magnetic field. These equations are then used to examine the problem of surface stability of a homogeneously pre-strained half-space subject to a magnetic field normal to its (plane) boundary. Related works on this subject include the studies of McCarthy [START_REF] Mccarthy | Wave propagation in nonlinear magnetothermoelasticity[END_REF], van de Ven [START_REF] Van De Ven | Interaction of Electromagnetic and Elastic Fields in Solid[END_REF], Boulanger [START_REF] Boulanger | Influence d'un champ magnétique statique sur la réflexion d'une onde électromagnétique à la surface d'un conducteur parfait élastique[END_REF][START_REF] Boulanger | Contribution à l'électrodynamique des continus élastiques et viscoélastiques[END_REF], Maugin [START_REF] Maugin | Wave motion in magnetizable deformable solids[END_REF], Carroll and McCarthy [START_REF] Carroll | Finite amplitude wave propagation in magnetized perfectly electrically conducting elastic materials[END_REF] and Das et al. [START_REF] Das | Magneto-visco-elastic surface waves in stressed conducting media[END_REF].

We adopt the formulation of Dorfmann and Ogden [START_REF] Dorfmann | Nonlinear magnetoelastic deformations[END_REF] as the starting point for the derivation of the incremental equations. This involves a total stress tensor and a modified strain energy function or total energy function, which enable the constitutive law for the stress to be written in a form very similar to that in standard nonlinear elasticity theory. The coupled governing equations then have a simple structure. We summarize these equations in Section 2. For incompressible isotropic magnetoelastic materials the energy density is a function of five invariants, which we denote here by I 1 and I 2 , the first two principal invariants of the Cauchy-Green deformation tensors, and I 4 , I 5 , I 6 , three invariants involving a Cauchy-Green tensor and the magnetic induction vector. This formulation is similar in structure to that associated with transversely isotropic elastic solids (see Spencer [START_REF] Spencer | Theory of invariants[END_REF]). The general incremental equations of nonlinear magnetoelasticity are then derived in Section 3. Therein we define the various magnetoelastic 'moduli' tensors and provide general incremental boundary conditions. Care is needed in deriving the boundary equations since the Lagrangian fields in the solid and the Eulerian fields in the vacuum must be reconciled.

Section 4 provides a brief summary of the basic equations associated with the pure homogeneous plane strain of a half-space of magnetoelastic material with a magnetic field normal to its boundary. In Section 5, the general incremental equations are applied to the analysis of surface stability. Not surprisingly, the resulting bifurcation criterion is a complicated equation, even when the pre-stress corresponds to plane strain and the magnetic induction vector is aligned with a principal direction of strain, as is the case here. The bifurcation equation comes from the vanishing of the determinant of a 7 × 7 matrix, which must be solved simultaneously with a bicubic equation. To present a tractable example, we therefore focus on a "Mooney-Rivlin magnetoelastic solid" for which the total energy function is linear in the invariants I 1 , I 2 , I 4 , and I 5 . Of course, these invariants are nonlinear in the deformation and the theory remains highly nonlinear. The bicubic then factorizes and a complete analytical resolution follows. In addition to the two elastic Mooney-Rivlin parameters (material constants), the material model involves two magnetoelastic coupling parameters. The stability behaviour of the halfspace depends crucially on the values of these coupling parameters and also on the magnitude of the magnetic field. In particular, a judicious choice of parameters can stabilize the half-space relative to the situation in the absence of a magnetic field. Equally, the half-space can become de-stabilized for different choices of the parameters. Thus, even this very simple model illustrates the possible complicated nature of the magnetoelastic coupling in the nonlinear regime.

The equations of nonlinear magnetoelasticity

In this section the equations for nonlinear magnetoelastic deformations, as developed by Dorfmann and Ogden [START_REF] Dorfmann | Nonlinear magnetoelastic deformations[END_REF][START_REF] Dorfmann | Magnetoelastic modelling of elastomers[END_REF][START_REF] Dorfmann | Nonlinear magnetoelastic deformations of elastomers[END_REF][START_REF] Ogden | Magnetomechanical interactions in magneto-sensitive elastomers[END_REF], are summarized for subsequent use in the derivation of the incremental equations. We consider a magnetoelastic body in an undeformed configuration B 0 , with boundary ∂B 0 . A material point within the body in that configuration is identified by its position vector X. By the combined action of applied mechanical loads and magnetic fields, the material is then deformed from B 0 to the configuration B, with boundary ∂B, so that the particle located at X in B 0 now occupies the position x = χ(X) in the deformed configuration B.

The function χ describes the static deformation of the body and is a oneto-one, orientation-preserving mapping with suitable regularity properties. The deformation gradient tensor F relative to B 0 is defined by F = Gradχ, F iα = ∂x i /∂X α , Grad being the gradient operator in B 0 . The magnetic field vector in B is denoted H, the associated magnetic induction vector by B and the magnetization vector by M .

To avoid a conflict of standard notations, the Cauchy-Green tensors are represented here by lower case characters; thus, the left and right Cauchy-Green tensors are b = F F t and c = F t F , respectively, where t denotes the transpose. The Jacobian of the deformation gradient is J = det F , and the usual convention J > 0 is adopted.

Mechanical equilibrium

Conservation of the mass for the material is here expressed as

Jρ = ρ 0 , (2.1)
where ρ 0 and ρ are the mass densities in the configurations B 0 and B, respectively. For an incompressible material, J = 1 is enforced so that ρ = ρ 0 . The equilibrium equation in the absence of mechanical body forces, is given in Eulerian form by divτ = 0, (

where τ is the total Cauchy stress tensor, which is symmetric, and div is the divergence operator in B. The total nominal stress tensor T is then defined by

T = JF -1 τ , (2.3)
so that the Lagrangian counterpart of the equilibrium equation (2.2) is

DivT = 0, (2.4)
Div being the divergence operator in B 0 .

Let N denote the unit outward normal vector to ∂B 0 and n the corresponding unit normal to ∂B. These are related by Nanson's formula nda = JF -t N dA, where dA and da are the associated area elements. The traction on the area element in ∂B may be written τ nda or as T t N dA. A traction boundary condition might therefore be expressed in the form

T t N = t a , (2.5)
where t a is the applied traction per unit reference area. If this is independent of the deformation then the traction is said to be a dead load.

Magnetic balance laws

In the Eulerian description, Maxwell's equations in the absence of time dependence, free charges and free currents reduce to divB = 0, curlH = 0, (

which hold both inside and outside a magnetic material, where curl relates to B. Thus, B and H can be regarded as fundamental field variables. A third vector field, the magnetization, when required, can be defined by the standard relation

B = µ 0 (H + M ). (2.7)
We shall not need to make explicit use of the magnetization in this paper. Associated with the equations (2.6) are the boundary continuity condi-

tions (B -B ⋆ ) • n = 0, (H -H ⋆ ) × n = 0, (2.8)
wherein B and H are the fields in the material and B ⋆ and H ⋆ the corresponding fields exterior to the material, but in each case evaluated on the boundary ∂B.

Lagrangian counterparts of B and H, denoted B l and H l , respectively, are defined by

B l = JF -1 B, H l = F t H, (2.9)
and in terms of these quantities equations (2.6) become

DivB l = 0, CurlH l = 0, (2.10)
where Curl is the curl operator in B 0 . We note in passing that a Lagrangian counterpart of M may also be defined, one possibility being M l = F t M . The boundary conditions (2.8) can also be expressed in Lagrangian form, namely

(B l -JF -1 B ⋆ ) • N = 0, (H l -F t H ⋆ ) × N = 0, (2.11)
evaluated on the boundary ∂B 0 .

Constitutive equations

There are many possible ways to formulate constitutive laws for magnetoelastic materials based on different choices of the independent magnetic variable and the form of energy function. For present purposes it is convenient to use a formulation involving a 'total energy function', or 'modified free energy function', which is denoted here by Ω, following Dorfmann and Ogden [START_REF] Dorfmann | Nonlinear magnetoelastic deformations[END_REF]. This is defined per unit reference volume and is a function of F and B l : Ω(F , B l ). This leads to the very simple expressions

T = ∂Ω ∂F , H l = ∂Ω ∂B l (2.12)
for a magnetoelastic material without internal mechanical constraints, and

T = ∂Ω ∂F -pF -1 , H l = ∂Ω ∂B l (2.13)
for an incompressible material, where p is a Lagrange multiplier associated with the constraint det F = 1. Note that the expression for H l is unchanged except that now det F = 1 in Ω. The Eulerian counterparts of the above equations are

τ = J -1 F ∂Ω ∂F , H = F -t ∂Ω ∂B l (2.14)
for an unconstrained material, where F -t = (F -1 ) t , and

τ = F ∂Ω ∂F -pI, H = F -t ∂Ω ∂B l , (2.15)
where I is the identity tensor. We emphasize that the first equation in each of (2.12)-(2.15) has exactly the same form as for a purely elastic material in the absence of a magnetic field.

Isotropic magnetoelastic materials

In general the mechanical properties of magnetoelastic elastomers have features that are similar to those of transversely isotropic materials. During the curing process a preferred direction is 'frozen in' to the material if the curing is done in the presence of a magnetic field, which aligns the magnetic particles. If cured without a magnetic field then the distribution of particles is essentially random and the resulting magnetoelastic response is isotropic. We focus on the latter case here for simplicity, but the corresponding analysis for the more general case follows the same pattern, albeit more complicated algebraically. A general constitutive theory for the former situation has been developed by Bustamante and Ogden [START_REF] Bustamante | Transversely isotropic nonlinearly magnetoelastic solids[END_REF] and applied to some simple problems. For isotropic materials, the energy function Ω depends only on c and B l ⊗ B l , through the six invariants

I 1 = tr c, I 2 = 1 2 (tr c) 2 -(tr c 2 ) , I 3 = det c = J 2 , I 4 = B l • B l , I 5 = (c B l ) • B l , I 6 = (c 2 B l ) • B l . (2.16)
For incompressible materials, I 3 = 1 and only the five invariants I 1 , I 2 , I 4 , I 5 , and I 6 remain. The total stress tensor τ is then expressed as (2.17) where Ω i = ∂Ω/∂I i , and the total nominal stress tensor T as

τ = -pI +2Ω 1 b+2Ω 2 (I 1 b-b 2 )+2Ω 5 B⊗B+2Ω 6 (B⊗bB+bB⊗B),
T = -pF -1 + 2Ω 1 F t + 2Ω 2 (I 1 F t -F t b) + 2Ω 5 B l ⊗ B + 2Ω 6 (B l ⊗ bB + F t B ⊗ B). (2.18)
Finally, the magnetic field vector H is found from (2.15) 2 as

H = 2(Ω 4 b -1 B + Ω 5 B + Ω 6 bB), (2.19)
and its Lagrangian counterpart is

H l = 2(Ω 4 B l + Ω 5 cB l + Ω 6 c 2 B l ).
(2.20)

Outside the material

In vacuum, there is no magnetization and the standard relation (2.7) reduces to

B ⋆ = µ 0 H ⋆ , (2.21)
where the star is again used to denote a quantity exterior to the material. Also, the stress tensor τ is now the Maxwell stress τ ⋆ , given by

τ ⋆ = µ -1 0 [B ⋆ ⊗ B ⋆ -1 2 (B ⋆ • B ⋆ )I], (2.22)
which, since divB ⋆ = 0 and curlB ⋆ = 0, satisfies divτ * = 0.

3 Incremental equations

Increments within the material

Suppose now that both the magnetic field and, within the material, the deformation undergo incremental changes (which are denoted by superposed dots). Let Ḟ and Ḃl be the increments in the independent variables F and B l . It follows from (2.12) that the increment Ṫ in T and the increment Ḣl in H l are given in the form

Ṫ = A Ḟ + Γ Ḃl , Ḣl = Γ Ḟ + K Ḃl , (3.1)
where A, Γ and K are, respectively, fourth-, third-and second-order tensors, with components defined by

A αiβj = ∂ 2 Ω ∂F iα ∂F jβ , Γ αiβ = ∂ 2 Ω ∂F iα ∂B l β = ∂ 2 Ω ∂B l β ∂F iα , K αβ = ∂ 2 Ω ∂B lα ∂B l β .
(3.2) We refer to these tensors as magnetoelastic moduli tensors. We note the symmetries

A αiβj = A βjαi , K αβ = K βα , (3.3)
and observe that Γ has no such indicial symmetry. The products in (3.1) are defined so that, in component form, we have

Ṫαi = A αiβj Ḟjβ + Γ αiβ Ḃl β , Ḣlα = Γ βiα Ḟiβ + K αβ Ḃl β . (3.4)
For an unconstrained isotropic material, Ω is a function of the six invariants I 1 , I 2 , I 3 , I 4 , I 5 , I 6 , and the expressions (3.2) can be expanded in the forms

A αiβj = 6 m=1, m =4 6 n=1, n =4 Ω mn ∂I n ∂F iα ∂I m ∂F jβ + 6 n=1, n =4 Ω n ∂ 2 I n ∂F iα ∂F jβ , Γ αiβ = 6 m=4 6 n=1, n =4 Ω mn ∂I m ∂B l β ∂I n ∂F iα + 6 n=5 Ω n ∂ 2 I n ∂F iα ∂B l β , K αβ = 6 m=4 6 n=4 Ω mn ∂I m ∂B lα ∂I n ∂B l β + n=4 Ω n ∂ 2 I n ∂B lα ∂B l β , (3.5)
where Ω n = ∂Ω/∂I n , Ω mn = ∂ 2 Ω/∂I m ∂I n . Expressions for the first and second derivatives of I n , n = 1, . . . , 6, are given in the Appendix. For an incompressible material, T is given by (2.13) 1 and its increment is then

Ṫ = A Ḟ + Γ Ḃl -ṗF -1 + pF -1 Ḟ F -1 , (3.6)
which replaces (3.1) in this case. On the other hand, H l is still given by (2.12) 2 and its increment is unaffected by the constraint of incompressibility, except, of course, since Ω is now independent of I 3 = 1, the summations in equations (3.5) omit m = 3 and n = 3.

It is now a simple matter to obtain the incremental forms of the (Lagrangian) governing equations. We have Div Ṫ = 0, Div Ḃl = 0, Curl Ḣl = 0.

(3.7)

These equations can be transformed into their Eulerian counterparts (indicated by a zero subscript) by means of the transformations

Ṫ 0 = J -1 F Ṫ , Ḃl0 = J -1 F Ḃl , Ḣl0 = F -t Ḣl (3.8)
(with J = 1 for an incompressible material), leading to div Ṫ 0 = 0, div Ḃl0 = 0, curl Ḣl0 = 0.

(3.9)

Now let u denote the incremental displacement vector x -X. Then, Ḟ = Gradu = (gradu)F , where grad is the gradient operator with respect to x. We use the notation d for the displacement gradient gradu, in components d ij = ∂u i /∂x j . From (3.8) and (3.1) we then have

Ṫ 0 = A 0 d + Γ 0 Ḃl0 , Ḣl0 = Γ 0 d + K 0 Ḃl0 , (3.10)
where, in index notation, the tensors A 0 , Γ 0 , and K 0 are defined by

A 0jisk = J -1 F jα F sβ A αiβk , Γ 0jik = F jα F -1 βk Γ αiβ , K 0ij = JF -1 αi F -1 βj K αβ (3.
11) for an unconstrained material. For an incompressible material J = 1 in the above and (3.10) is replaced by

Ṫ 0 = A 0 d + Γ 0 Ḃl0 + pd -ṗI, Ḣl0 = Γ 0 d + K 0 Ḃl0 , (3.12)
and the incremental incompressibility condition is divu = 0. (3.13)

Notice that A 0 and K 0 inherit the symmetries of A and K, respectively, so that

A 0jisk = A 0skji , K 0ij = K 0ji . (3.14)
Finally, using the incremental form of the rotational balance condition F T = (F T ) t , we find that Γ 0 has the symmetry

Γ 0ijk = Γ 0jik , (3.15)
and we uncover the connections

A 0jisk -A 0ijsk = τ js δ ik -τ is δ jk (3.16)
between the components of the tensors A 0 and τ for an unconstrained material (see, for example, Ogden [START_REF] Chadwick | On the definition of elastic moduli[END_REF] for the specialization of these in the purely elastic case), and

A 0jisk -A 0ijsk = (τ js + pδ js )δ ik -(τ is + pδ is )δ jk (3.17)
for incompressible materials (see Chadwick [START_REF] Chadwick | The application of the Stroh formalism to prestressed elastic media[END_REF] for the elastic specialization).

Following Prikazchikov [START_REF] Prikazchikov | Surface and Edge Phenomena in Pre-Stressed Incompressible Elastic Solids[END_REF], we decompose the tensor A 0 into the sum

A 0 = A (0) 0 + A (5) 0 + A (6)
0 .

(3.18)

The first term A (0) 0 does not involve any derivatives with respect to I 4 , I 5 , and I 6 . Clearly, this term is very similar to the tensor of elastic moduli associated with isotropic elasticity in the absence of magnetic fields. In component form it is given by

JA (0) 0jisk = 4b ij b ks Ω 11 + 4N ij N ks Ω 22 + 4J 4 δ ij δ ks Ω 33 + 4(b ks N ij + b ij N ks )Ω 12 +4(b ks δ ij + b ij δ ks )Ω 13 + 4J 2 (N ks δ ij + N ij δ ks )Ω 23 +2δ ik b js Ω 1 + 2(2b ij b ks + δ ik N js -b jk b is -b ik b js )Ω 2 +2J 2 (2δ ij δ ks -δ is δ jk )Ω 3 , (3.19)
where

N ij = b kk b ij -b ik b kj (3.20)
and b ij are the components of b.

The terms A

0jisk and A

0jisk may be expressed in the forms

A (5) 0jisk = A 0(5) 0jisk Ω 5 + 6 m=1, m =4 A m(5) 0jisk Ω m5 , A (6)
0αiβj = A 0(6) 0jisk Ω 6 + 6 m=1, m =4 A m(6) 0jisk Ω m6 , (3.21)
where A 4(5) 0jisk = 0 and A

0(5) 0jisk = 2J -1 a j a s δ ik , A 1(5) 0jisk = 4J -1 (a k a s b ij + a i a j b ks ), A 2(5) 0jisk = 4J -1 (a k a s N ij + a i a j N ks ), A 3(5) 0jisk = 4J(a k a s δ ij + a i a j δ ks ), A 5(5) 0jisk = 4J -1 a j a i a s a k , A 6(5) 0jisk = 2J -1 (a j a i H ks + a k a s H ij ), (3.22)
with

H ij = a j a k b ik + a i a k b jk , a i = F iα B lα . (3.23) Similarly, A 4(6) 0jisk = 0 and A 0(6) 0jisk = 2J -1 (δ ik H js + a i a s b jk + a j a k b is + a j a s b ik + a i a k b js ), A 1(6) 0jisk = 4J -1 (b ks H ij + b ij H ks ), A 2(6) 0jisk = 4J -1 (H ij N ks + H ks N ij ), A 3(6) 0jisk = 4J(H ks δ ij + H ij δ ks ), A 5(6) 0jisk = 2J -1 (a i a j H ks + a k a s H ij ), A 6(6) 0jisk = 4J -1 H ij H ks . (3.24)
The tensor Γ 0 is decomposed as

Γ 0 = Γ (1) 0 + Γ (2) 0 + Γ (3) 0 + Γ (5) 0 + Γ (6) 0 , (3.25)
with components given by

Γ (1) 0jik = 4b ij M 1k , Γ (2)
0jik = 4N ij M 2k , Γ (3)
0jik = 4J 2 δ ij M 3k , Γ (5) 0jik = 4a j a i M 5k + 2(a j δ ik + a i δ jk )Ω 5 , Γ (6) 0jik = 4H ij M 6k + 2(δ ik a s b js + a i b jk + δ jk a s b is + a j b ik )Ω 6 , (3.26)
where

M ik = F -1 αk B lα Ω i4 + a k Ω i5 + a j b jk Ω i6 . (3.27)
Finally, we represent K 0 in the form

K 0 = K (4) 0 + K (5) 0 + K (6) 0 , (3.28)
with components

K (4) 0ij = 2JF -1 αi (2B lα M 4j + F -1 αj )Ω 4 , K (5) 0ij = 2J(2a i M 5j + δ ij Ω 5 ), K (6) 0ij = 2J(2a k b ik M 6j + b ij Ω 6 ).
(3.29)

For an incompressible material, the above expressions are unaltered except that J = 1 and all the terms Ω 3 and Ω n3 , n = 1, . . . , 6 in A 0ijkl , Γ 0ijk , K 0ij are omitted.

Outside the material

The standard relation B = µ 0 H in vacuum is incremented to

Ḃ⋆ = µ 0 Ḣ⋆ , (3.30)
where Ḃ⋆ and Ḣ⋆ are the increments of B ⋆ and H ⋆ , respectively. These fields satisfy Maxwell's equations div Ḃ⋆ = 0, curl Ḣ⋆ = 0.

(3.31)

Finally, we increment the Maxwell stress of (2.22) to

τ ⋆ = µ -1 0 [ Ḃ⋆ ⊗ B ⋆ + B ⋆ ⊗ Ḃ⋆ -(B ⋆ • Ḃ⋆ )I], (3.32)
noting that div τ ⋆ = 0.

Incremental boundary conditions

At the boundary of the material, in addition to any applied traction t a (defined per unit reference area), there will in general be a contribution from the Maxwell stress exterior to the material. This is a traction τ ⋆ n per unit current area and can be 'pulled back' to the reference configuration to give a traction Jτ ⋆ F -t N per unit reference area, in which case the boundary condition (2.5) is modified to

T t N = Jτ ⋆ F -t N + t a . (3.33)
On taking the increment of this equation, we obtain

Ṫ t N = J τ ⋆ F -t N -Jτ ⋆ F -t Ḟ t F -t + Jτ ⋆ F -t N + ṫa , (3.34)
and hence, on updating this from the reference configuration to the current configuration,

Ṫ t 0 n = τ ⋆ n -τ ⋆ d t n + (divu)τ ⋆ n + ṫa . (3.35)
Proceeding in a similar fashion for the other fields, we increment the magnetic boundary conditions (2.11) to give, again after updating,

( Ḃl0 + dB ⋆ -(divu)B ⋆ - Ḃ⋆ ) • n = 0 (3.36) and ( Ḣl0 -d t H ⋆ - Ḣ⋆ ) × n = 0. (3.37)
4 Pure homogeneous deformation of a halfspace

Here we summarize the basic equations for the pure homogeneous deformation of a half-space in the presence of a magnetic field normal to its boundary prior to considering a superimposed incremental deformation in Section 5.

The deformed half-space

Let X 1 , X 2 , X 3 be rectangular Cartesian coordinates in the undeformed half-space B 0 and take X 2 = 0 to be the boundary ∂B 0 , with the material occupying the domain X 2 ≥ 0. In order to minimize the number of parameters, we consider the material to be incompressible and subject to a plane strain in the (X 1 , X 2 ) plane. With respect to the Cartesian axes, the deformation is then defined by

x 1 = λX 1 , x 2 = λ -1 X 2 , x 3 = X 3 .
The components of the deformation gradient tensor F and the right Cauchy-Green tensor c are written F and c, respectively, and are given by

F =   λ 0 0 0 λ -1 0 0 0 1   , c =   λ 2 0 0 0 λ -2 0 0 0 1   , (4.1)
where λ is the principal stretch in the X 1 direction. The invariants I 1 and I 2 are therefore

I 1 = I 2 = 1 + λ 2 + λ -2 . (4.2)
We take the magnetic induction vector B to be in the x 2 direction and to be independent of x 1 and x 3 . It then follows from divB = 0 that its component B 2 is constant. Thus,

B 1 = 0, B 2 = 0, B 3 = 0. ( 4.3)
The associated Lagrangian field B l = F -1 B then has components

B l1 = 0, B l2 = λB 2 , B l3 = 0, (4.4)
and the invariants involving the magnetic field are

I 4 = B 2 l2 , I 5 = λ -2 I 4 , I 6 = λ -4 I 4 . (4.5)
We may now compute the stress field using (2.17), (4.1) and (4.4). The resulting non-zero components of τ are

τ 11 = 2Ω 1 λ 2 + 2Ω 2 (λ 2 + 1) -p, τ 22 = 2Ω 1 λ -2 + 2Ω 2 (1 + λ -2 ) -p + 2Ω 5 λ -2 I 4 + 4Ω 6 λ -4 I 4 , τ 33 = 2Ω 1 + 2Ω 2 (λ 2 + λ -2 ) -p. (4.6)
The magnetic field H has components given by (2. [START_REF] Carroll | Finite amplitude wave propagation in magnetized perfectly electrically conducting elastic materials[END_REF] as

H 1 = 0, H 2 = 2(Ω 4 + λ -2 Ω 5 + λ -4 Ω 6 )λB l2 , H 3 = 0. (4.7)
Since B l2 and λ are constant, all the fields are uniform and the equilibrium equations and Maxwell's equations are satisfied. In view of (4.2) and (4.5), there are only two independent variables, λ and I 4 . We thus introduce a specialization ω(λ, I 4 ) of the total energy Ω, by the definition

ω(λ, I 4 ) = Ω(1 + λ 2 + λ -2 , 1 + λ 2 + λ -2 , I 4 , λ -2 I 4 , λ -4 I 4 ), (4.8)
from which it follows that

ω λ = 2λ -1 [(λ 2 -λ -2 )(Ω 1 + Ω 2 ) -λ -2 I 4 Ω 5 -2λ -4 I 4 Ω 6 ]
,

ω 4 = Ω 4 + λ -2 Ω 5 + λ -4 Ω 6 , (4.9)
where ω λ = ∂ω/∂λ, ω 4 = ∂ω/∂I 4 . Hence,

τ 11 -τ 22 = λω λ , H 2 = 2λB l2 ω 4 . (4.10)

Outside the material

From the boundary conditions (2.8) applied at the interface

x 2 = X 2 = 0, we have B ⋆ 2 = B 2 and H ⋆ 1 = H ⋆ 3 = 0, while from (2.21) it follows that B ⋆ 1 = B ⋆ 3 = 0 and H ⋆ 2 = µ -1 0 B ⋆ 2 = µ -1 0 B 2 .
Outside the material we take the magnetic field to be uniform and equal to its interface value, Maxwell's equations are then satisfied identically, B ⋆ therefore has components

B ⋆ 1 = 0, B ⋆ 2 = B 2 = λ -1 B l2 , B ⋆ 3 = 0, (4.11)
and H ⋆ has components

H ⋆ 1 = 0, H ⋆ 2 = µ -1 0 B 2 = µ -1 0 λ -1 B l2 , H ⋆ 3 = 0. (4.12)
From these expressions, we deduce that the non-zero components of the Maxwell stress (2.22) are given by

τ ⋆ 11 = -τ ⋆ 22 = -1 2 µ -1 0 B 2 2 = -1 2 µ -1 0 λ -2 I 4 = τ ⋆ 33 . (4.13)
The applied mechanical traction on x 2 = 0 required to maintain the plane strain deformation has a single non-zero component τ 22τ ⋆ 22 .

Surface stability

We now address the question of surface stability for the deformed half-space by establishing a bifurcation criterion based on the incremental static solution of the boundary-value problem. Biot [START_REF] Biot | Mechanics of Incremental Deformations[END_REF] initiated this approach, which has since been successfully applied to a great variety of boundary-value problems; see Ogden [START_REF] Ogden | Elements of the theory of finite elasticity[END_REF] for pointers to the vast literature on the subject.

Magnetoelastic moduli

First we note that since F ij = 0 for i = j and B l1 = B l3 = 0 several simplifications occur in the expressions for the components of the magnetoelastic moduli tensors A 0 , Γ 0 , K 0 . In particular, we have

A 0iijk = 0, K 0ij = 0, for j = k, Γ 0ii3 = Γ 03ii = Γ 0ii1 = Γ 01ii = 0, Γ 0ijk = 0, for i = j = k = i. (5.1)
For subsequent use we compute the quantities

a = A 01212 , 2b = A 01111 + A 02222 -2A 01221 -2A 01122 , c = A 02121 , d = Γ 0211 , e = Γ 0222 -Γ 0112 , f = K 011 , g = K 022 . (5.2)
Explicitly, we obtain

a = 2λ 2 (Ω 1 + Ω 2 ) + 2I 4 Ω 6 , b = (λ 2 + λ -2 )(Ω 1 + Ω 2 ) + I 4 [λ -2 Ω 5 + (6λ -4 -2)Ω 6 ] + 2(λ 4 + λ -4 -2)(Ω 11 + 2Ω 12 + Ω 22 ) + 4I 4 (λ -4 -1)[Ω 15 + Ω 25 + 2λ -2 (Ω 16 + Ω 26 )] + 2I 2 4 λ -4 (Ω 55 + 4λ -2 Ω 56 + 4λ -4 Ω 66 ), c = 2λ -2 (Ω 1 + Ω 2 ) + 2I 4 [λ -2 Ω 5 + (2λ -4 + 1)Ω 6 ], d = 2B l2 λ[λ -2 Ω 5 + (λ -4 + 1)Ω 6 ], e = 4B l2 λ -1 [Ω 5 + 2λ -2 Ω 6 + (1 -λ 4 )(Ω 14 + Ω 24 ) + (λ -2 -λ 2 )(Ω 15 + Ω 25 ) + (λ -4 -1)(Ω 16 + Ω 26 ) + I 4 (Ω 54 + λ -2 Ω 55 + 2λ -2 Ω 46 + 3λ -4 Ω 56 + 2λ -6 Ω 66 )], f = 2(λ -2 Ω 4 + Ω 5 + λ 2 Ω 6 ) g = 2(λ 2 Ω 4 + Ω 5 + λ -2 Ω 6 ) + 4I 4 (λ 2 Ω 44 + 2Ω 45 + 2λ -2 Ω 46 + λ -2 Ω 55 + 2λ -4 Ω 56 + λ -6 Ω 66 ).
(5.3)

In terms of the energy density ω(λ, I 4 ) we have the connections

a -c = λω λ , 2(b + c) = λ 2 ω λλ , e = -2B l2 λ 2 ω λ4 , g = 2λ 2 (ω 4 + 2I 4 ω 44 ), ( 5

Incremental fields and equations

We seek incremental solutions depending only on the in-plane variables x 1 and x 2 such that u 3 = 0 and Ḃl03 = 0. Hence u i = u i (x 1 , x 2 ) and Ḃl0i =

From the six possible roots we select s 1 , s 2 , s 3 to be the three roots satisfying (5.27). We then construct the general solution for the solid as

ψ = 3 j=1
A j e -ks j x 2 e ikx 1 , φ = k

3 j=1 D j e -ks j x 2 e ikx 1 , (5.30)
where A j , D j , j = 1, 2, 3, are constants.

For the half-space x 2 ≤ 0 (vacuum) we take a solution φ ⋆ to (5.18) that is localized near the interface x 2 = 0. Specifically, we write this as

φ ⋆ = ikC ⋆ e kx 2 e ikx 1 , (5.31)
where C ⋆ is a constant. The constants A j and D j are related through either equation in (5.28). From the second equation, for instance, we obtain

s j (ds 2 j + d -e)A j + (f s 2 j -g)D j = 0, j = 1, 2 , 3; no summation. (5.32)
We also have the two traction boundary conditions (5.21) and (5.22), which read (5.33) and

(c -τ 22 -1 2 µ -1 0 λ -2 I 4 )(A 1 + A 2 + A 3 ) + c(s 2 1 A 1 + s 2 2 A 2 + s 2 3 A 3 ) -d(s 1 D 1 + s 2 D 2 + s 3 D 3 ) -λ -1 B l2 C ⋆ = 0,
(τ 22 -1 2 µ -1 0 λ -2 I 4 -2b -c)(s 1 A 1 + s 2 A 2 + s 3 A 3 ) + c(s 3 1 A 1 + s 3 2 A 2 + s 3 3 A 3 ) + (e -ds 2 1 )D 1 + (e -ds 2 2 )D 2 + (e -ds 2 3 )D 3 -λ -1 B l2 C ⋆ = 0. (5.34)
Finally, the two magnetic boundary conditions (5.24) and (5.25) become

λ -1 B l2 (s 1 A 1 + s 2 A 2 + s 3 A 3 ) -(D 1 + D 2 + D 3 ) + µ 0 C ⋆ = 0, (5.35)
and

(d -µ -1 0 λ -1 B l2 )(A 1 + A 2 + A 3 ) + d(s 2 1 A 1 + s 2 2 A 2 + s 2 3 A 3 ) -f (s 1 D 1 + s 2 D 2 + s 3 D 3 ) -C ⋆ = 0. (5.36)
In total, there are seven homogeneous linear equations for the seven unknowns A j , D j , j = 1, 2, 3, and C ⋆ . The resulting determinant of coefficients must vanish and this equation is rather formidable to solve, particularly since it must be solved in conjunction with the bicubic (5.29). It is in principle possible to express the determinant in terms of the sums and products s 1 +s 2 +s 3 , s 1 s 2 + s 2 s 3 + s 3 s 1 , s 1 s 2 s 3 , and to find these from the bicubic (5.29), similarly to the analysis conducted in the purely elastic case (see Destrade et al. [START_REF] Destrade | Nonprincipal surface waves in deformed incompressible materials[END_REF]). However, the resulting algebraic expressions rapidly become too cumbersome for this approach to be pursued.

Instead, we propose either (a) to turn directly to a numerical treatment once Ω has been determined by curve fitting from experimental data for a given magnetoelastic solid, or (b) to use a simple form for Ω that allows some progress to be made.

Regarding approach (a), we remark that, as emphasized by Dorfmann and Ogden [START_REF] Dorfmann | Nonlinear magnetoelastic deformations[END_REF][START_REF] Dorfmann | Magnetoelastic modelling of elastomers[END_REF][START_REF] Dorfmann | Nonlinear magnetoelastic deformations of elastomers[END_REF][START_REF] Ogden | Magnetomechanical interactions in magneto-sensitive elastomers[END_REF], there is a shortage of, and a pressing need for, suitable experimental data and for the derivation of functions Ω from such data. In the next section we focus primarily on the analytical approach (b).

Example: a "Mooney-Rivlin magnetoelastic solid"

As a prototype for the energy function Ω, we propose

Ω = 1 4 µ(0)[(1 + γ)(I 1 -3) + (1 -γ)(I 2 -3)] + µ -1 0 (αI 4 + βI 5 ), (5.37)
where µ(0) is the shear modulus of the material in the absence of magnetic fields and α, β, γ are dimensionless material constants, α and β being magnetoelastic coupling parameters. For α = β = 0, (5.37) reduces to the strain energy of the elastic Mooney-Rivlin material, a model often used for elastomers.

In respect of (5.37) the stress τ in (2.17) reduces to

τ = -pI + 1 2 µ 0 (1 + γ)b + 1 2 µ 0 (1 -γ)(I 1 b -b 2 ) + 2µ -1 0 βB ⊗ B, (5.38)
while H in (2. [START_REF] Carroll | Finite amplitude wave propagation in magnetized perfectly electrically conducting elastic materials[END_REF]) becomes

H = 2µ -1 0 (αb -1 B + βB).
(5.39)

Clearly, equation (5.38) shows that the parameter α does not affect the stress. By contrast β, if positive, stiffens the material in the direction of the magnetic field, i.e. a larger normal stress in this direction is required to achieve a given extension in this direction than would be the case without the magnetic field. On the other hand, by reference to (5.39), we see that α provides a measure of how the magnetic properties of the material are influenced by the deformation (through b). If β = 0 the stress is unaffected by the magnetic field. On the other hand, if α = 0 then the magnetic constitutive equation (5.39) is unaffected by the deformation. Thus, a two-way coupling requires inclusion of both constants. The quantities defined in (5.2) and (5.3) now reduce to

a = µ(0)λ 2 , 2b = µ(0)(λ 2 + λ -2 + βλ -2 Ī4 ), c = µ(0)(λ -2 + βλ -2 Ī4 ), d = µ -1 0 µ(0)βλ -1 Bl2 , e = 2 µ -1 0 µ(0)βλ -1 Bl2 , f = µ -1 0 (αλ -2 + β), g = µ -1 0 (αλ 2 + β), (5.40)
where Bl2 , a dimensionless measure of the magnetic induction vector amplitude, and Ī4 are defined by Bl2 = B l2 / µ 0 µ(0), Ī4 = B2 l2 .

( and it follows that the relevant roots are

s 1 = 1, s 2 = λ 2 , s 3 = λ αλ 2 + β α + βλ 2 + αβ Ī4 .
(5.44) Note that for s 3 to be real for all λ > 0 and all Bl2 , the inequalities α ≥ 0, β > 0 or α > 0, β ≥ 0 (5.45) must hold. (The case in which there is no magnetic field corresponds to α = β = 0.) It is assumed here that these inequalities are satisfied, so that s 3 is indeed a qualifying root satisfying (5.27).

The equation (5.32) becomes

s j (s 2 j -1)βλ -1 Bl2 Âj -[(αλ -2 + β)s 2 j -αλ 2 -β] Dj = 0, j = 1, 2, 3, (5.46)
where Âj = µ -1 0 µ(0)A j , Dj = µ -1 0 D j , (5.47)

and the s j are given by (5.44).

Next, consider the four remaining boundary conditions (5.33)- (5.36). In order to keep the number of parameters to a minimum (so far, we have λ, Bl2 , α, β), and to make a simple connection with known results for the surface stability of an elastic Mooney-Rivlin material, we assume that there is no applied mechanical traction on the boundary x 2 = 0, and hence

τ 22 = τ ⋆ 22 = 1 2 µ -1 0 λ -2 I 4 .
(5.48)

The boundary conditions (5.21)-(5.25) now read

[1 + (β -1) Ī4 ]( Â1 + Â2 + Â3 ) + (1 + β Ī4 )(s 2 1 Â1 + s 2 2 Â2 + s 2 3 Â3 ) -βλ Bl2 (s 1 D1 + s 2 D2 + s 3 D3 ) -λ Bl2 C ⋆ = 0, (λ 4 + 2 + 2β Ī4 )(s 1 Â1 + s 2 Â2 + s 3 Â3 ) -(1 + β Ī4 )(s 3 1 Â1 + s 3 2 Â2 + s 3 3 Â3 ) + βλ Bl2 [(s 2 1 -2) D1 + (s 2 2 -2) D2 + (s 2 3 -2) D3 ] + λ Bl2 C ⋆ = 0, Bl2 (s 1 Â1 + s 2 Â2 + s 3 Â3 ) -λ( D1 + D2 + D3 ) + λC ⋆ = 0, λ Bl2 (β -1)( Â1 + Â2 + Â3 ) + λβ Bl2 (s 2 1 Â1 + s 2 2 Â2 + s 2 3 Â3 ) -(α + βλ 2 )(s 1 D1 + s 2 D2 + s 3 D3 ) -λ 2 C ⋆ = 0.
(5.49)

From the seven equations (5.46) and (5.49), we have derived a bifurcation criterion (vanishing of the determinant of coefficients) using a computer algebra package, but it is too long to reproduce here. It is a complicated rational function of the four parameters λ, Bl2 , α, β. However, it is easy to solve numerically, and for the numerical examples we fix the material parameters α and β and find the critical stretch λ cr in compression as a function of Bl2 . For Bl2 = 0, we recovered the well-known critical compression stretch for surface instability of the elastic Mooney-Rivlin material in plane strain, namely λ cr = 0.5437 [START_REF] Biot | Mechanics of Incremental Deformations[END_REF], as expected. For Figure 1a (Figure 1b), we set α = 0.5 (α = 2.0) and curves for β = 0.0, 0.5, 1.0, 1.5, 2.0 are shown. We found that λ cr is an even function of Bl2 and we therefore restricted attention to positive Bl2 (within the range 0 ≤ Bl2 ≤ 3). The behaviour as Bl2 becomes larger and larger (not shown here) indicates that the half-space becomes more and more unstable in compression. Moreover, it can even become unstable in tension (λ cr > 1). The figures also clearly demonstrate that for some values of α, β, and Bl2 the critical stretch ratio is smaller than that for the purely elastic case (λ cr < 0.5437), in which cases the magnetic field has a stabilizing effect.

Figure 1 Turning back to a phenomenological approach, we remark that the energy function (5.37) has quite good curve-fitting qualities for moderate fields. There are four parameters at hand, namely µ(0), α, β, γ, two of which, µ(0) and β, may be determined from shear tests. Indeed Dorfmann and Ogden [START_REF] Ogden | Magnetomechanical interactions in magneto-sensitive elastomers[END_REF] show that in general the shear modulus for isotropic nonlinear magnetoelasticity is 2[Ω 1 + Ω 2 + I 4 Ω 5 + I 4 Ω 6 (3 + 2κ 2 )], where κ is the amount of shear in a simple shear test. Here the modulus is independent of κ and is given by µ(B l2 ) = µ(0) + 2µ -1 0 βI 4 .

(5.50)

This highlights the role of β in increasing the mechanical stiffness of the material -through the shear modulus. Jolly et al. [START_REF] Jolly | A model of the behaviour of magnetorheological materials[END_REF] conducted double lap shear tests on magneto-sensitive elastomers containing 10, 20, and 30% by volume of iron particles. From their Figure 7, we see that in the range 0 ≤ B l2 ≤ 0.5 Tesla, the variations of µ(B l2 ) resemble those of a parabolic profile such as the one suggested by (5.50). For the 10% iron by volume elastomer specimen, Table 1 in Jolly et al. [START_REF] Jolly | A model of the behaviour of magnetorheological materials[END_REF] gives µ(0) = 0.26 MPa, and at B l2 = 0.5 Tesla, we read off their Figure 7 that µ(0.5)µ(0) ≃ 0.07 MPa, indicating that β ≃ 0.18. Similarly, for the 20% and the 30% iron by volume elastomer specimens we find β ≃ 0.53 and β ≃ 0.72, respectively. Figure 2a (Figure 2b) illustrates the variation of the critical compression stretch with the amplitude of the dimensional magnetic induction vector, from 0 to 0.5 Tesla, for the 20% (30%) iron by volume elastomer, and for several values of α. We remark than the presence of the magnetic field makes the two specimens slightly more stable than in the purely elastic case because all the critical compression stretch values are smaller than 0.5437. It is also clear that increasing the value of α makes the half-space more stable. However, it is worth noting that the 30% iron by volume specimen is slightly less stable than the 20% iron by volume specimen for the same values of α. Figure 2 A Derivatives of the invariants with respect to F and B l

We derive the expressions for the first derivatives of the six invariants with respect to F , The second derivatives of the invariants are computed as follows: first, the second derivatives with respect to F ,

∂I 1 ∂F iα = 2F
∂ 2 I 1 ∂F iα ∂F jβ = 2δ ij δ αβ , ∂ 2 I 2 ∂F iα ∂F jβ = 2(2F iα F jβ -F iβ F jα + c γγ δ ij δ αβ -b ij δ αβ -c αβ δ ij ), ∂ 2 I 3 ∂F iα ∂F jβ = 4I 3 F -1 αi F -1 βj -2I 3 F -1 αj F -

. 4 )

 4 whereω λλ = ∂ 2 ω/∂λ 2 , ω λ4 = ∂ 2 ω/∂λ∂I 4 and ω 44 = ∂ 2 ω/∂I 2 4 .

  .41) Note the connections 2b = a + c, e = 2d. (5.42) Now we find that the bicubic (5.29) factorizes in the form (s 2 -1)(s 2λ 4 )[αλ 4 + βλ 2 -(α + βλ 2 + αβ Ī4 )s 2 ] = 0, (5.43)

Figure 1 .

 1 Figure1. Dependence of the critical stretch λ cr < 1 for instability in compression for a magnetoelastic Mooney-Rivlin solid in plane strain on the non-dimensional measure Bl2 of the magnetic field for several values of the magnetoelastic coupling parameters α and β.

Figure 2 .

 2 Figure2. Dependence of the critical stretch λ cr < 1 for instability in compression for a magnetoelastic Mooney-Rivlin solid in plane strain with the dimensional measure B l2 of the magnetic field, for several values of the magnetoelastic coupling parameters α and β.

Figure 1 :Figure 2 :

 12 Figures

  iα , ∂I 2 ∂F iα = 2(c γγ F iαc αγ F iγ ), iγ B lγ c αβ B lβ + F iγ c γβ B lβ B lα ), (A.1)and with respect to B l ,

	∂I 3 ∂F iα ∂I 6 ∂F iα ∂B lα = 0, = 2I 3 F -1 αi , ∂I 2 ∂B lα = 2(F ∂I 1 ∂I 4 ∂B lα = 2B lα , ∂I 5 ∂B lα	∂I 4 ∂F iα = 0, = 2c αβ B lβ , = 0,	∂I 5 ∂F iα ∂I 3 ∂B lα ∂I 6 ∂B lα	= 2B lα (F iγ B lγ ), = 0, = 2c αγ c γβ B lβ .	(A.2)

  1 βi , ∂ 2 I 4 ∂F iα ∂F jβ = 0, ∂ 2 I 5 ∂F iα ∂F jβ = 2δ ij B lα B lβ , ∂ 2 I 6 ∂F iα ∂F jβ = 2[δ ij (c αγ B lγ B lβ + c βγ B lγ B lα ) + δ αβ F iγ B lγ F jδ B lδ + F iγ B lγ F jα B lβ + F jγ B lγ F iβ B lα + b ij B lα B lβ ]; (A.3)next, the mixed derivatives with respect to F and B l ,∂ 2 I 1 ∂F iα ∂B l β = 0, ∂ 2 I 2 ∂F iα ∂B l β = 0, ∂ 2 I 3 ∂F iα ∂B l β = 0, ∂ 2 I 4 ∂F iα ∂B l β = 0, ∂ 2 I 5 ∂F iα ∂B l β = 2δ αβ F iγ B lγ + 2B lα F iβ , ∂ 2 I 6 ∂F iα ∂B l β = 2F iβ c αγ B lγ + 2F iγ B lγ c αβ + 2F iγ c γβ B lα + 2δ αβ F iγ cγδ B lδ ; (A.4) finally, the second derivatives with respect to B l , ∂B lα ∂B lβ = 2δ αβ , ∂ 2 I 5 ∂B lα ∂B lβ = 2c αβ , ∂ 2 I 6 ∂B lα ∂B lβ = 2c αγ c γβ . (A.5)

	∂ 2 I 1 ∂B lα ∂B lβ	= 0,	∂ 2 I 2 ∂B lα ∂B lβ	= 0,	∂ 2 I 3 ∂B lα ∂B lβ	= 0,
	∂ 2 I 4

Ḃl0i (x 1 , x 2 ) for i = 1, 2 and ṗ = ṗ(x 1 , x 2 ). In the following, a subscripted comma followed by an index i signifies partial differentiation with respect to x i , i = 1, 2.

The incremental version (3.13) of the incompressibility constraint reduces here to u 1,1 + u 2,2 = 0, (5.5) and hence there exists a function ψ = ψ(x 1 , x 2 ) such that

Similarly, equation (3.9) 3 reduces to Ḃl01,1 + Ḃl02,2 = 0, (5.7)

The incremental equations of equilibrium (3.9) 1 simplify to Ṫ011,1 + Ṫ021,2 = 0, Ṫ012,1 + Ṫ022,2 = 0.

(5.9)

From the identities (5.1), the only non-zero components of the incremental stress Ṫ 0 are found to be

(5.10) Also, equation (3.9) 2 reduces to Ḣl01,2 -Ḣl02,1 = 0, (5.11) wherein are the only non-zero components of Ḣl , which, from (5.1), are given by

(5.12) In terms of the functions ψ and φ equations (5.9) and (5.11) become for ψ and φ.

Outside the material

In vacuum, Maxwell's equations (3.31) hold for Ḃ and Ḣ. From the second equation, and the assumption that all fields depend only on x 1 and x 2 , we deduce the existence of a scalar function

(5.19)

Boundary conditions

We now specialize the general incremental boundary conditions of Section 3.3 to the present deformed semi-infinite solid. First, for ṫa = 0, the incremental traction boundary conditions (3.35) reduce to

on x 2 = 0. Putting together the results of this section, using (4.13), (5.2), (5.6), (5.8), (5.10), (5.16) and (5.19), we express the two equations (5.20) as

and

which apply on x 2 = 0. In obtaining the latter we have differentiated (5.20) 2 with respect to x 1 and made use of (5.13) 1 .

Next, the incremental magnetic boundary conditions (3.36) and (3.37) reduce to

on x 2 = 0. Using again the results of the preceding sections, we write these as

and (µ

on x 2 = 0.

Resolution

We are now in a position to solve the incremental boundary value problem.

We seek small-amplitude solutions, localized near the interface x 2 = 0. Hence we take solutions in the solid (x 2 ≥ 0) to be of the form

where k > 0 (2π/k is the wavelength of the perturbation) and s is such that ℜ(s) > 0 (5.27)

to ensure decay with increasing x 2 (> 0). Substituting (5.26) into the incremental equilibrium equations (5.


Bernard Nabzar

P L J Bourbiaux

F Zitha

Roggero

S M Hosseini-Nasab

; F Douarche

; M Simjoo

; L Nabzar

; B Bourbiaux

; P L J Zitha

F Roggero

Numerical Simulation of Foam Flooding in Porous Media in the Absence and Presence of Oleic Phase

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Gas injection for enhanced oil recovery (EOR) suffers from poor sweep efficiency. Three reasons are associated with this deficiency of gas flooding: 1) segregation and gravity override due to the lower density of gas compared to oil and/or water, 2) viscous fingering due to high mobility ratio between injected gas and oil and/or water, 3) channelling through high-permeability streaks or layers in heterogeneous and layered reservoir [START_REF] Lake | Fundamentals of Enhanced Oil Recovery[END_REF].

Although not all the gas injection methods essentially lead to the poor oil recovery efficiency.

The recent method of the Gas-Assisted Gravity Drainage (GAGD) process has a good sweep efficiency and higher recovery than Continuous Gas Injection (CGI) and Water-Alternative Gas (WAG) processes, in both immiscible and miscible modes [START_REF] Al-Mudhafar | Lessons Learned from the Field-Scale Simulation of the Gas-Assisted Gravity Drainage GAGD Process in Heterogeneous Sandstone Oil Reservoirs[END_REF].

Foam can improve the volumetric sweep efficiency by reducing gas mobility, providing a favourable mobility ratio and contacting a larger fraction of the reservoir to mitigate the effect of heterogeneity, gas segregation and viscous instability [START_REF] Rossen | Foams in Enhanced Oil Recovery[END_REF][START_REF] Farajzadeh | Investigation of immiscible and miscible foam for enhancing oil recovery[END_REF]Simjoo et al., 2011). Design of the foam EOR process for field-scale application requires accurate prediction and description of foam behaviour in porous media with and without the oleic phase. Modelling of foam flow in a porous media can be categorized into two major methods: the local equilibrium with implicit texture (LE-IT) foam model, and the population balance approach. The LE-IT foam model assumes that a local steady state of foam dynamics in terms of creation and destruction of foam is reached instantaneously wherever gas and surfactant (as a foaming agent) coexist in porous media [START_REF] Rossen | Modeling Foams for Acid Diversion[END_REF][START_REF] Boeije | Fitting Foam Simulation Model Parameters to Data[END_REF][START_REF] Ma | Modeling Techniques for Foam Flow in Porous Media[END_REF]). This model implicitly takes into account the presence of foam generation and coalescence through a mobility interpolation factor which depends itself on water saturation, gas velocity and other factors [START_REF] Cheng | Simulating Foam Processes at High and Low Foam Qualities[END_REF].

Application of LE-IT foam model for the field scale requires dependency of the model parameters upon the variety of geological properties of the reservoir, in particular permeability, porosity or fracture geometry in each direction [START_REF] Farajzadeh | Effect of Permeability on Implicit-Texture Foam Model Parameters and the Limiting Capillary Pressure[END_REF]. The LE-IT empirical foam model is unable to model the correct foam density as a function of foam quality where the foam quality can vary in thick reservoirs based upon various factors such as permeability, pore size distribution etc.

The population balance approach describes the dynamics of foam generation and destruction where the foam mobility reduction is based on the bubble size and bubble density (foam texture) [START_REF] Kovscek | A mechanistic population balance model for transient and steady-state foam flow in Boise sandstone[END_REF]2010). Although population balance model provides a comprehensive mechanistic description of foam flow in porous media, it requires many physical parameters which are not easy to determine for the field application. On the other hand, the empirical LE-IT model is more pragmatic and simpler for the field-scale simulation as it requires fewer simulation parameters and numerical difficulties and lower computational cost compared to the population balance model [START_REF] Kapetas | Effect of Permeability on Foam-model parameters -An Integrated Approach from Coreflood Experiments through to Foam Diversion Calculations[END_REF]Rossen & Boeije, 2015).

Much work in the literature has focused on the simulation of foam flooding in the absence of an oleic phase. [START_REF] Simjoo | Modeling of Foam Flow Using Stochastic Bubble Population Model and Experimental Validation[END_REF] studied the transient foam flow in an oil-free porous media by using the stochastic bubble population model. The premise of this model is that foam flow in porous media is a complex fluid and bubble generation is a stochastic process. They obtained a good match between the numerically calculated fluid saturation and pressure data with those obtained from the experiments at which foam was generated by coinjecting nitrogen and alpha olefin sulfonate surfactant in Bentheimer sandstone. [START_REF] Boeije | Fitting Foam-Simulation-Model Parameters to Data : I . Coinjection of Gas and Liquid[END_REF] proposed a method to calculate the parameters of the water-saturation-dependent function and also shear-thinning function from the foam pressure gradient data at low and high-quality regimes at fixed total velocity. This method can provide the initial estimates for the foam model physical parameters to be used in the reservoir simulator for foam simulation at a large scale. [START_REF] Ma | Estimation of Parameters for the Simulation of Foam Flow through Porous Media[END_REF] estimated the parameters of the water-saturation-dependent function to describe the dry-out effect in the absence of the oil phase. In their approach, shearthinning behaviour at the low-quality regime of foam flooding was ignored.

Effect of permeability variation of porous media on the behaviour of foam flow in the absence of oil was studied experimentally and theoretically by [START_REF] Kapetas | Effect of Permeability on Foam-model parameters -An Integrated Approach from Coreflood Experiments through to Foam Diversion Calculations[END_REF]. They showed permeability can have a significant impact on the critical foam saturation such that the higher permeability layer exhibits lower critical water saturation ( * ). However, they did not come up with a robust correlation between the permeability and transition abruptness characteristic of the LE-IT foam model for the foam quality-scan experiments at different foam-flow regimes. [START_REF] Jones | Effect of surfactant concentration on foam: From coreflood experiments to implicittexture foam-model parameters[END_REF] fitted the core-flood results of foam flooding without oil, for different surfactant concentration by the LE-IT foam model. To be able to predict the effect of the concentration on the foam apparent viscosity, they extend the model such that five foam parameters vary with surfactant concentration. [START_REF] Lotfollahi | A Four-Phase Flow Model to Simulate Chemical EOR with Gas[END_REF] presented a numerical model to simulate foam flooding in the presence of micro-emulsion phase. However, in this work no validation was given with respect to the experimental data. Similarly, [START_REF] Lashgari | Modeling of Low-Tension Surfactant-Gas Flooding Process in a Four-Phase Flow Simulator[END_REF] applied the black-oil model system coupled with the micro-emulsion phase behaviour model for simulation of low-tension gas flooding. They used the interfacial tension (IFT) reduction as the main factor to control the incremental oil recovery even though this mechanism in foam flooding does not always function. Therefore, most of the works in the literature have modelled foam flow in porous media either in the absence or the presence of oleic phase. On the other hand, many experimental data of foam flooding for EOR purposes have been reported.

Thus, this study aimed to investigate the application of LE-IT foam model for numerical modelling of foam flow in the sandstone rock for both in the absence and presence of an oleic phase. The foam model of Puma-Flow reservoir simulator (IFP Energies nouvelles) was used, which is similar to the foam model in CMG-STARS simulator (Computer Modelling Group Ltd., 2007). To this end, first the parameters of the dry-out and gas velocity functions in the foam model were determined by a least-square matching of the model to the experimental data obtained at low and high-quality regimes at a constant superficial velocity. Then, numerical simulations were conducted to investigate the effect of surfactant concentration on the pressure drop across the core sample for the core-flood laboratory data of foam flooding in the absence of oil. Thereafter, the effect of oil on the modelling of foam flooding was elaborated by fitting foam parameters to experimental data of the foam flooding in the presence of oil. The structure of this paper is as follows. First, we describe the main features of the foam model used in this study. Next we present an overview of the experimental study on foam flooding in sandstone porous media with and without oleic phase. The paper proceeds with the simulation results obtained from numerical modelling of the experimental data, and finally the main conclusions are drawn.

Theoretical Description of LE-IT Foam Model

Features of local equilibrium and implicit-texture (LE-IT) foam model, in the PumaFlow simulator (in-house reservoir simulator of IFPEN), are described as follows. Reduction of gas mobility due to presence of a foaming agent is assigned to the relative permeability function, while gas-phase viscosity is assumed unchanged no matter whether a foaming agent is present or not. The relative permeability reduction factor is interpolated between a (maximum) reference value known or measured in optimal foaming conditions and a unity value in the absence of the foaming agent. Gas relative permeability in the presence of foam, , is rescaled to gas relative permeability in the absence of foam, , by multiplying by a dimensionless interpolation factor, FM.

(1)

The interpolation function (FM), also known as the mobility reduction factor, is a product of various functions to capture the different physical effects on the foam stability in porous media. The FM function generally includes the contributions of four variables, namely surfactant concentration, water saturation, gas phase velocity and oil saturation.

1 1 1 (2)
The term M ref is the (reference) maximum foam mobility reduction factor obtained without the impact of other factors (F 1 , F 2 …). The functions F 1 , F 2 , F 3 and F 4 incorporate the effects of foaming agent concentration, water saturation, oil saturation and capillary number (gas phase velocity) respectively. The F i values are within the range [0; 1]: the closer they are to 1 the more efficient the foam will be. Through the indices one can decide to neglect the dependency on some parameters by simply inputting e i = 0.

, , , , (5)

, ( 6
)
F 1 is a power-law relationship to describe the effect of foaming agent concentration. It is controlled by parameters ̅ , and e S , where ̅ , is the critical surfactant concentration above which gas mobility is independent of surfactant concentration (C s ). ̅ , depends on surfactant type and is larger than the critical micelle concentration (CMC). When surfactant concentration is larger than the CMC value, foam becomes more stable due to the increase of disjoining pressure [START_REF] Kovscek | A mechanistic population balance model for transient and steady-state foam flow in Boise sandstone[END_REF][START_REF] Schramm | Foams: Fundamentals and Applications in the Petroleum Industry[END_REF]. Note that foam coalescence rate is lower than foam film creation rate, as long as the disjoining pressure is larger than the local capillary pressure [START_REF] Buchavzov | A disjoining pressure study of foam films stabilized by mixtures of nonionic and ionic surfactants[END_REF][START_REF] Ekserova | Foam and foam films : theory, experiment, application[END_REF][START_REF] Schulze-Schlarmann | A disjoining pressure study of foam films stabilized by tetradecyl trimethyl ammonium bromide C14TAB[END_REF].

F 2 is namely foam dry-out function which describes the dependence of foam strength on water saturation S w . * in the function F 2 is the critical water saturation at which the maximum foam strength is reached in a high-quality regime ('foam quality' referring to gas fractional flow), and below which foam weakens in a low-quality regime [START_REF] Kam | A Model for Foam Generation in Homogeneous Media[END_REF]. For the water saturation lower than * as foam quality increases, foam (in a highquality regime) begins to significantly dry out and coarsen by the 'limiting capillary pressure, * [START_REF] Khatib | Effects of Capillary Pressure on Coalescence and Phase Mobilities in Foams Flowing Through Porous Media[END_REF]. During a foam quality-scan test with a fixed superficial velocity, the high-quality regime is the range of foam quality where the pressure gradient decreases with increasing foam quality, while in the low-quality regime the pressure gradient increases with increasing foam quality [START_REF] Tanzil | Conditions for Foam Generation in Homogeneous Porous Media[END_REF]. In the F 2 function, f w is the dryout coefficient that controls the collapse rate of the foam in the high-quality regime as a function of water saturation [START_REF] Kapetas | Effect of Permeability on Foam-model parameters -An Integrated Approach from Coreflood Experiments through to Foam Diversion Calculations[END_REF]. This coefficient (f w ) has a significant effect on the predictions of pressure gradient particularly near the transition between the low/high-quality regimes.

Larger values of f w yield a sharper transition between the high and low-quality regimes, which means that foam dries out faster, while small values cause a gradual transition between the foam flow regimes, which means that the foam gradually dries out. If the transition between high and low-quality regimes is abrupt, then the * parameter corresponds to the water saturation at the limiting capillary pressure P c *, i.e. the water saturation at which foam strength collapses [START_REF] Gauglitz | Foam generation in homogeneous porous media[END_REF]. F 3 represents the destabilization effect of oil on the foam stability. The * parameter is the oil saturation threshold above which foam is destroyed and gas mobility is not reduced at all. The e o is the exponent of the function of F 3 , which controls the rate of disappearance of the foam when oil is present. By setting a zero value for this exponent, one can neglect any detrimental effect of the presence of oil on the generated foam and thus, on the gas mobility control properties.

F 4 is the shear-thinning velocity effect that represents the influence of non-Newtonian power-law shear thinning of foam rheology on the foam strength. It has been seen that the apparent foam viscosity decreases as the velocity increases, which makes the foam look like a non-Newtonian shear thinning fluid [START_REF] Hirasaki | Mechanisms of Foam Flow in Porous Media: Apparent Viscosity in Smooth Capillaries[END_REF][START_REF] Xu | Effective viscosity of foam in periodically constricted tubes[END_REF]. In F 4 function, N c is the capillary number, and 0< <1 is a range of reference capillary number values. N c is a dimensionless number, which is defined by the ratio of viscous forces to capillary forces. It is defined by the following equation:

/ ⁄ (7)
where is the total Darcy velocity (gas+liquid), is the apparent viscosity of the displacing fluid in the porous media, is porosity of rock, k is absolute permeability of rock, is the pressure gradient, and is the surface tension between gas and foaming agent solution, which is a function of the concentration of foaming agent in the aqueous phase.

Other definitions of the capillary number have been used in the literature [START_REF] Lake | Fundamentals of Enhanced Oil Recovery[END_REF], but Eq. ( 7) is the most commonly used for foam flooding [START_REF] Simjoo | Novel Insight into Foam Mobility Control[END_REF]. The capillary pressure has an important impact on the foam stability, and in this study, the capillary number is used to model this effect.

Since the maximum value of function F 4 is equal to 1, therefore the value of should be set equal to the lowest capillary number expected in the simulations. Below this value of capillary number (i.e.

), the shear-thinning behaviour is not represented. The parameter e C controls the significance of the shear-thinning effect; in which this parameter the larger it is, the stronger the shear-thinning behaviour becomes. A value of e C = 0 represents Newtonian behaviour. It is an advantage for foam to be a shear thinning fluid for improved sweep efficiency purposes in reservoir applications. This is because near the injection wellbore where the velocity is high, the apparent foam viscosity will be low and thus, the injectivity will be high [START_REF] Ashoori | Dynamic foam behavior in the entrance region of a porous medium[END_REF]. However, deep in the reservoir, the apparent foam viscosity will be high, and therefore foam will act as an efficient mobility control agent.

Overview of the Core-flood Experiments

Materials and Methods

Brine was prepared by adding sodium chloride (NaCl, Merck) at a fixed concentration of 3 wt% in de-ionized water (pH=6.8±0.1). Surfactant solutions with different surfactant concentrations were prepared with the same salinity as the brine. Normal hexadecane (n-C 16 )

with a density of 0.78±0.01 g/cm 3 was used as model oil. The surface and interfacial tensions were measured using a KSV Sigma tensiometer by the DuNouy ring method. Surface tension of 1.0 wt%AOS surfactant solution was 28.0 ±0.1 mN/m, and IFT between the surfactant solution and model oil was found to be 1.9 ±0.1 mN/m (both measured at 20 o C). The critical micelle concentration of AOS solution in the presence of brine was 4.0×10 -3 wt%. The properties of the chemicals used in this work are summarized in Table 1, and the physical properties of the core samples are presented in Table 2. Details of the experiment conditions and results can be found elsewhere [START_REF] Simjoo | Effects of Oil on Foam Generation and Propagation in Porous Media[END_REF].

Core-flooding Set-up and Procedure

The set-up used to perform the core-floods is shown schematically in Figure 1. The sequences and conditions used to conduct the core-flooding experiments are summarized in Table 3. First, air from the core sample was removed by flushing it with CO 2 at 5 bar injection pressure. Then the dry core sample was saturated by injecting approximately 10 pore volumes (PV) of brine while increasing back-pressure up to 30 bar. Back pressure was increased to ensure that any CO 2 present in the core was dissolved into the water phase (100% core saturation with brine). For each individual experiment after core was fully saturated with brine, absolute brine permeability was measured by the standard method [START_REF] Christiansen | Literature Review and Recommendation of Methods for Measuring Relative Permeability of Anhydrite from the Salado Formation at the Waste Isolation Pilot Plant[END_REF]. For the experiments in absence of oil, core was pre-flushed with a surfactant prior to foam flooding. Surfactant pre-flush was done to satisfy the adsorption capacity of the rock surface and thus to shorten the delay in foam generation. A different injection sequence was undertaken for the experiments mimicking foam EOR: after full core saturation with brine, first oil was injected into the core till connate water saturation was reached (primary drainage). Then core was water-flooded till residual oil was reached (imbibition). Next similarly to the experiments without oil, core was pre-flushed with a surfactant solution and then subjected to foam flooding. Nitrogen gas (N 2 ) and surfactant solution were co-injected at a fixed superficial velocity of 4.58 ft/day and foam quality 91%, to generate foam in-situ at a back-pressure of 20 bar at ambient temperature (21 ± 1°C). Foam quality , f g , is given by f g = q g /(q g + q l ) × 100, where q l is flow rate of the liquid phase (surfactant solution), and q g is flow rates of the gas phase (nitrogen).

Gas superficial velocity was determined from its nominal value by using two corrections:

(1) with respect to the calibration factor for the mass flow controller (MFC) towards an adjustment of the pressure before the MFC and injection pressure after the MFC (2) with respect to the adjustment of the effect of gas compressibility of N 2 due to 30 bar backpressure. For the foam-quality scan experiment with the low and high-quality regimes, the foam quality was controlled by varying the relative rates of injection of N 2 gas and AOS solution at a constant superficial velocity.

Results and Discussion

Numerical Modelling of Foam Flooding in Absence of Oil

The objective of this section is to model the dynamic of foam generation and propagation in porous media in the absence of oil at the transient and steady-state flow conditions. We then model the effect of surfactant concentration on the pressure-gradient build-up during the foam flow. We illustrated the dynamics of foam generation and propagation in the porous media by the mobility reduction of full-strength foam, surfactant concentration effect, the limiting water saturation at which foam collapses, and the parameters governing foam destruction at limiting capillary pressure. To this end, we demonstrated the determination of dedicated parameters to the aforementioned characteristics of foam flow for the LE-IT foam model. First, parameters of the water saturation dependent function (F 1 ), and also the parameters of capillary number dependent function (F 4 ) were obtained by modelling the pressure gradient core-flood data of one scan of N 2 -foam quality at a fixed superficial velocity. In the absence of oil at a fixed surfactant concentration and total velocity, the LE-IT foam model related the foam mobility reduction factor, FM presented in Eq. ( 2), only to two functions of water saturation and capillary number as follows:

1

1 1 0.5 * max , (8)
In the absence of foam, Corey-type relative permeability for water and gas phases were used as follows:

1 (9) 1 1 (10)
where and are the endpoint relative permeabilities for water and gas, respectively, and are the corresponding exponents, is connate water saturation and is residual gas saturation. The LE-IT foam model captures the value of each function to a maximum of unity. Thus, the value of the reference capillary number ( ) was set equal to the lowest capillary number expected in the simulations, where below this value ( ), shear-thinning behaviour in the low-quality regime does not exist (Boeije & Rossen, 2015).

This implies that was not considered a foam matching parameter while other four independent parameters , * , , and were considered fitting parameter in the low and high-quality regimes of foam flooding in porous media. To determine these parameters, a non-linear least-square optimization was used in which all the four foam parameters were computed simultaneously. For this, an initial guess and an allowed range were considered for each parameter. Table 4 presents the input parameters of rock-fluid properties, relative permeability parameters of gas and liquid phases (Eqs. ( 9) and ( 10)), and conditions of the foam flooding experiment that were implemented in the numerical simulation. In order to compute the capillary number, a function describing the dependence of the gaswater surface tension versus surfactant concentration was implemented. Then, the stabilizing effect of the surfactant concentration on the generated foam was modelled using the dependence of capillary number due to the variation of gas-water surface tension. The values of the foam matching parameters obtained by least-square optimization method are summarized in Table 5. Numerical simulation of the foam flooding pressure drop of foam-quality scan experiment without oil is shown in Figure 2. It demonstrates a very good fit of the simulations to the experimental data was obtained. 5. Numerically simulated behaviour of foam-quality scans at a fixed superficial velocity that are higher and lower superficial velocities than the experimental condition are also shown in this figure.

Calculated pressure drops from numerical simulation for higher and lower total superficial velocity than the experimental results are also presented in Figure 2. It illustrates that the rheological characteristics and foam-flow regimes in a foam quality-scan experiment with the different velocities can also be captured by the numerical results, where the pressure gradient first increased with increasing foam quality, and then the pressure gradient decreased due to the foam dry-out. The constant of driving foam evolution (f w ) of foam dry-out function was not considered a large value for the set of data examined in this work, because the transition from low-quality to high-quality regime is not abrupt and the generated foam does not collapse at a constant water saturation as shown in Figure 2 [START_REF] Kapetas | Effect of Permeability on Foam-model parameters -An Integrated Approach from Coreflood Experiments through to Foam Diversion Calculations[END_REF]Rossen & Boeije, 2015).

In order to test the capability of the LE-IT foam model to describe the effect of the surfactant concentration, numerical simulation of the pressure drop for foam flooding at two different concentrations was performed. Five foam model parameters of dry-out function (F 2 )

and gas velocity effect (F 4 ) determined from the previous simulation of foam flood at low/high-quality regimes were used. Then, the reference surfactant concentration and the exponent of the surfactant-effect function (F 1 ) of the foam model, by fitting to the experimental data of pressure drop of foam flow at two different concentrations were estimated as earlier presented in Table 5. The foam existence was considered to be a function of the surfactant (foaming agent) concentration and water saturation. Figure 3 shows numerically calculated pressure drop versus experimental data at transient and steady-state regimes for 0.5 wt% and 1.0 wt% surfactant concentrations. state conditions. However, initiation of foam generation and propagation by the model are slightly higher than those obtained experimentally. Simulated pressure drop profile in the case of the lower surfactant concentration (0.5 wt%) at the earlier time of the transient state mimics less features of the experimental results due to the delayed foam generation.

Simulation of Foam Flooding in Presence of Oil

The objective of this section is the numerical simulation of the core-flood experimental data of immiscible foam flooding during the displacement of the water-flooded remaining oil.

Modelling foam propagation through the porous media in the presence of oil is essential to predict the performance of the immiscible foam EOR process. In this section, we assumed that the destabilizing effect of oil on the foam is a function of oil and water saturations. The total foam coalescence rate is the summation of the coalescence rate due to the water dry-out near the limiting water saturation (Eq. ( 4)), and also due to the presence of oil (Eq. ( 5)).

Table 6 summarizes the input parameters, and other conditions of the foam flooding experiment that were implemented in the numerical simulation. The reservoir model used for the vertical one-dimensional simulation of the core-flow system had dimensions 1×1×100 grid blocks. No-flow boundries were imposed on the lateral sides of the core sample along with in-and out-flow conditions. Reservoir model was isothermal with the temperature constant at 22°C, the initial pressure of 20 bar, the porosity of 23% and the permeability was 2.1 Darcy. Immediately before foam flooding, the reservoir model was at the water-flooded residual oil saturation (S or =0.44), with no gas initially present in the core. The injector and producer were placed at opposite ends of the reservoir model, and injection rate was performed at a rate of 4.2 ft/day. The maximum bottom-hole pressure (BHP) in the injector and producer was set 50 bar and 30 bar, respectively. N 2 gas and surfactant solution were co-injected at a fixed gas fraction into the flow system of the reservoir model with a uniform water fraction of 0.09 (foam quality of 91%). The binary relative permeability curves between gas/oil and water/oil were used for the 1D displacement simulation. A Corey-type equation was used to describe the relative permeabilities. The binary relative permeability curves were combined into a ternary relative permeability function using a three-phase model of Stone I (Stone, 1973). The three-phase oil relative permeability was computed with the Stone I model and applied to calculate the residual oil saturation parameter. The effects of oil saturation on foam was incorporated by implementing the parameters of F 3 function (Eq. ( 5)) as presented in Table 7.

It is essential to note that the F i in the FM function were multiplied together. Then, if several parameters are not favorable, then the foam stability will be strongly impacted. Figure 4 shows a good match of the simulated and measured oil recovery by immiscible foam flooding. However, there is a small discrepancy between the numerically calculated and measured oil recovery data between 1.5 and 3.0 PV, most likely due to the lower oil relative permeability during the three phase flow of oil displacement experiment by foam flooding [START_REF] Heins | Oil Relative Permeability During Enhanced Oil Recovery by Foam Flooding Coreflood Experiments[END_REF][START_REF] Simjoo | Modeling of Foam Flow Using Stochastic Bubble Population Model and Experimental Validation[END_REF]. Thus, this plot (Figure 4) demonstrates an acceptable fit to the series of pressure drops and oil recover factor at both the transient and steady-state conditions of foam generation and propagation through the porous media. Figure 5 compares the oil-cut (normalized oil production rate) profile obtained from numerical simulation with those obtained from the experiments. The numerically predicted oil-cut profile is in a good agreement with the one obtained from the experiment. In the first 1.2 pore volumes, the oil-cut increases until it reaches a peak value of 54 % and then decreases towards 5.8% over more than 15 pore volumes. Visual inspection of the effluents showed that in the first 1.8 PV clean oil was produced while after oil was produced as an oilin-water emulsion [START_REF] Simjoo | Effects of Oil on Foam Generation and Propagation in Porous Media[END_REF]. This means that foam flooding first induces the formation of an oil bank and then, after the breakthrough of the trailing edge of the oil bank, the diffuse oil is obtained as a dispersed phase. The above oil production features due to foam flooding, i.e. oil bank followed by a long tailing oil recovery, were captured remarkably well by the numerical simulations. The tailing oil production occurred during the highest capillary number (or pressure build-up, see Figure 5), and produced oil was attributed to the emulsified oil. The desaturation of the oleic phase from the core can also be elaborated from the series of CT scan images obtained during foam flooding as presented in Figure 6. The red color corresponds to the liquid phase consisting of residual oil plus surfactant solution. As gas and surfactant solution were co-injected from the top of the core, the intensity of the orange colour

Conclusion

The study presented in this work demonstrated an integrated approach for the numerical simulation of foam transport in porous media with and without the oleic phase. First, the LE- Two distinct regimes of incremental oil production were captured by the numerical simulation, first by the formation of an oil bank and then followed by a long tail production.

The results of this study showed that the developed numerical model reproduced the main feature of the transient and steady-state foam flow regimes in presence of oleic phase. This was supported by a good match between pressure and oil recovery profiles obtained numerically with those obtained from the foam flooding EOR experiment in the Bentheimer sandstone core samples. The small discrepancy between the numerically calculated and measured oil recovery data, most likely was due to the lower oil relative permeability during the three phase flow of oil displacement experiment by foam flooding.

Nomenclature

SI units are assumed for all parameters used in calculations.

,

reference

Figure 1 :

 1 Figure 1: Schematic of the experimental set-up used to perform the core-flooding experiments. The core-holder was held vertically on the bench of the CT scanner.

Foam

Figure 2 :

 2 Figure 2: History match of the pressure drop of foam-quality scan experiment without oil at fixed total superficial velocity (U g + U w ) (at 4.58 ft/day and 91% foam quality) in Bentheimer core by the LE-IT foam model and the simulation parameters presented in Table5. Numerically simulated behaviour of foam-quality scans at a fixed superficial velocity that are higher and lower superficial velocities than the experimental condition are also shown in this figure.

Figure 3 :

 3 Figure 3: Comparison of numerically calculated and measurement pressure drops of foam flooding in the absence of oil at two different surfactant concentrations. Total superficial velocity and foam quality at injection side of the core was fixed at 4.58 ft/day and 91%, respectively.

Figure 4

 4 Figure 4 compares the pressure drop and oil recovery factor (with respect to OIIP) obtained from numerical simulation and the experiments. The results show a good qualitative agreement between the simulated and measured pressure drop data. The simulated pressure drop increased until reaching a minimum-pressure gradient to generate strong foam in the transient state, and after pressure drop raised, then levelled off to a plateau value during the steady-state condition of foam flow regime. One can see small fluctuations in the pressure drop after the chemical breakthrough for the experimental data and modelling results. Also

Figure 4 :

 4 Figure 4: Comparison of numerically calculated and measured pressure drops and oil recovery factor (OIIP) during oil displacement by foam flooding.

Figure 5 :

 5 Figure 5: Comparison of numerically calculated and measured oil-cut by foam flooding.

  diminishes from the left to right progressively in favour of more blue, corresponding threephase flow. This gives a qualitative impression of the change in fluid saturations including oil saturation in the core.

Figure 6 :

 6 Figure 6: CT images obtained during foam flooding in a core containing water-flooded residual oil[START_REF] Simjoo | Effects of Oil on Foam Generation and Propagation in Porous Media[END_REF]. The red colour stands for the water-flooded core including residual oil saturation, and the blue colour indicates the presence of the foam phase. Foam was injected from top to bottom to ensure gravity stable conditions. Foam breakthrough occurred at 0.63±0.02 PV.

  IT foam model was matched to the foam core-flood experiments in the absence of the oleic phase by taking into account the foam dry-out phenomena, rheological characteristics along with shear thinning properties and foam-quality regimes. The parameters of the watersaturation and the shear-rate dependent functions (F 2 and F 4 ) were determined by historymatching the foam quality-scan data at the steady-state condition, in the low/high quality regimes. The corresponding parameters of F 2 and F 4 functions were then used to model the effect of surfactant concentration on the pressure build-up induced by foam. The parameters of the surfactant-concentration dependent function (F 1 ) were obtained by matching the numerical and experimental pressure drop data for the transient and steady-state conditions.The numerical simulation of foam flooding in the presence of water-flooding residual oil was performed by estimating the oil-saturation dependent function of the foam mobility interpolation factor. The numerically calculated pressure drop fitted well the experimental lab data of the transient and steady-state foam generation and propagation in the presence of oil.

Table 1 :

 1 Properties of the materials used to conduct the experiments at the ambient Temperature 21°C and atmospheric pressure

	Materials	Formulation	Molecular weight (g/mole)	Viscosity (cP)	Active (wt%) content	Supplier
	Salt	NaCl	58.50	solid	99.98	Merck
	Oil	n-C16	226.00	3.3 ± 0.01	99.99	Sigma
	Surfactant	AOS	315.00	1.08 ± 0.01	40.00	Stepan
	Gas	N 2	28.01	0.017	99.98%	3M
	Rock	mainly SiO 2	--	--	93.00%	Bentheim mine

Table 2 :

 2 Physical properties of the core samples used in the

	core-flooding experiment
	Core sample	Bentheimer
	Length (cm)	17.0±0.1
	Diameter (cm)	3.8 ±0.1
	Porosity (%)	21.0 ±0.1
	Pore volume (cm 3 )	42.5 ±0.5
	Core density (g/ cm 3 )	2.65 ±0.02
	Brine permeability (Darcy)	2.5 ±0.1

Table 3 :

 3 Overview of the experimental procedure used for the foam flooding experiments with and without oil

	Injection step sequence	Flowrate (cm 3 /min)	Back pressure (bar)	Injection direction

Table 4 :

 4 Input values of rock, fluids and relative permeability parameters in the simulation

	Input parameters	Values
	Number of grid blocks	1 1 100
	Grid block size (m m m)	0.0367 0.0367 0.0017
	Brine permeability	2.30 (Darcy)
	Porosity	0.25
	Temperature	21°C
	Brine viscosity	1.0 (cP)
	Brine density	1.0 (gr/cm 3 )
	Water compressibility factor	5.0E -5
	Gas viscosity	0.01744-0.01832 (cP)
	Surfactant concentration	1.0, 0.5, wt%
	Salinity (water flooding)	3.0 wt% (NaCl)
	Injection flowrate	1.1 (cm 3 /min)
	Producer bottom hole pressure	20 bar
	Foam injection time	20 (PV)
	Connate (irreducible) water saturation	0.15
	Residual gas saturation	0
	Water relative permeability endpoint	0.25
	Gas relative permeability endpoint	0.76
	Exponent of water relative permeability	2
	Exponent of gas relative permeability	1.3

Table 5 :

 5 Foam model parameters for the numerical modelling of foam flooding pressure drop at low/high-quality regimes in the absence of an oil phase

	Parameters

Table 6 :

 6 Input values of the parameters of the reservoir model in the simulation

	Input parameters

Table 7 :

 7 Foam model parameters for numerical simulation of the pressure drop and oil recovery by foam flooding in the presence of oil

	Parameters	Values
	Reference mobility reduction factor, M Ref	45000
	Critical oil saturation, *	0.5
	Exponent of oil effect function, e o	0.2
	Reference surfactant concentration, ,	0.35 wt%
	Exponent of surfactant-dependent function, e s	0.5
	Reference capillary number,	9.78E -6

Acknowledgement

Shell Global Solutions B.V. is gratefully acknowledged for partially funding of the laboratory experiments presented in this paper through its Gamechanger programme. Our especial thanks to Descartes French Excellence Scholarships funded by TOTAL company, for the financial support to conduct this research at the IFP Energy Nouvel (IFPEN) in France.

We would like to thank M. Friebl, E, Meivogel and H. Asten for technical supports at Dietz Laboratory of the Geoscience and Engineering department of TU Delft.