Rubén Ibáñez

Adrien Scheuer
email: adrien.scheuer@uclouvain.be

Emmanuelle Abisset-Chavanne
email: emmanuelle.abisset-chavanne@ec-nantes.fr

Francisco Chinesta
email: francisco.chinesta@ensam.eu

Antonio Huerta
email: antonio.huerta@upc.edu

Rub Én Ib Á Ñez

Roland Keunings
email: roland.keunings@uclouvain.be

Emmanuelle Abisset

A simple microstructural viscoelastic model for flowing foams

Keywords: Flowing foams, Viscoelasticity, Conformation, Microstructural description Francisco Chinesta

Introduction

Aqueous foams are concentrated dispersions of gas bubbles in a surfactant solution. Their structures are organized over a large range of length scales and complex flows take place at different scales [START_REF] Cohen-Addad | Flow in foams and flowing foams[END_REF].

The proposal of macroscopic constitutive equations allows for the efficient modelling and simulation of industrial processes involving the flow of foams [START_REF] Benito | An elasto-visco-plastic model for immortal foams or emulsions[END_REF][START_REF] Cheddadi | Numerical modelling of foam Couette flows[END_REF]. Usually, such descriptions remain however too phenomenological, and even though they predict accurately the flow kinematics, microstructure information remains often unaccessible. On the opposite side, fully microscopic simulations allow for very detailed descriptions of the foam microstructural evolution [START_REF] Bikard | Numerical prediction of the foam structure of polymeric materials by direct 3D simulation of their expansion by chemical reaction based on a multidomain method[END_REF]. However such approaches fail to address scenarios of industrial interest that usually involve the flow of foams in very large and complex 3D geometries.

The macroscopic flow model is expected to depend on the cellular structure: cell size, shape and orientation, as well as on the fluid rheology and the surface tension. Moreover, cell shape and orientation are induced by the flow. This microscopic information could be introduced into a macroscopic flow model by using standard upscaling and homogenization techniques. Thus inspired by [START_REF] Feyel | A multilevel finite element method (FE2) To describe the response of highly non-linear structures using generalized continua[END_REF], at some locations in the domain in which an effective homogeneous fluid flows, we could attach a representative volume containing several cells, whose size, shape and orientation depend on the considered location. Now, a detailed microscopic calculation could be carried out in order to determine the effective fluid rheology. However, such a route, widely and successfully considered in a variety of fields, remains expensive from a computational point of view despite some attempts at combining it with advanced model reduction techniques [START_REF] Lamari | Routes for efficient computational homogenization of nonlinear materials using the Proper Generalized Decomposition[END_REF].

The most appealing description consists of a macroscopic flow model making use of some conformational variables describing the main microstructural features, as widely considered in the field of multiscale polymer modelling [START_REF] Binetruy | F Flows in Polymers, Reinforced Polymers and Composites. A multiscale approach[END_REF][START_REF] Keunings | Micro-macro methods for the multiscale simulation viscoelastic flow using molecular models of kinetic theory[END_REF]. Thus, in [START_REF] Karimi | PUFOam: a novel open-source CFD solver for the simulation of polyurethane foams[END_REF] the authors study flows during foaming, considering cell evolutions but without addressing the shape orientation and then, the induced anisotropy. Richer microstructure descriptions can be obtained by using a set of configurational coordinates, from which a conformation tensor can be derived and a macroscopic constitutive equation established. This route was successfully considered in [START_REF] Tlili | Mechanical formalisms for tissue dynamics[END_REF]. In the present work, we propose an alternative simpler microstructure description and its coupling with the macroscopic flow.

Remark

Cell conformation

A very simple description of a 3D cell consists of a deformable ellipsoid (with constant volume) ranging from the spherical shape to the infinite aspect ratio ellipsoid (rod). Such an ellipsoid could be represented by means of three orthogonal extensible springs, with reference length 2L 0 and stiffness K. In the sequel, we restrict our analysis to 2D scenarios (by considering ellipses instead of ellipsoids) but the derived models and their numerical solution procedures can be straightforwardly extended to 3D. First, we consider the kinematics of a single linear elastic dumbbell as starting point for elaborating the cell conformation.

Kinematics of an extensible rod

The extensible rod, of reference length 2L 0 and assumed aligned in direction p (p having a unit norm), is represented by an elastic spring of length 2L (in the deformed state) and stiffness K equipped with two beads at its extremities where hydrodynamic forces act. In the sequel, the word hydrodynamic refers to the viscous drag force and not to the one considered in other works to describe the effects of a bead kinematics on the others from the use of the Oseen tensor. These forces scale with the fluid -bead relative velocity, the former given by v 0 + ∇v • pL and the latter by v G + ṗL + p L, where v 0 is the unperturbed fluid velocity at the rod center of gravity and v G the velocity of the rod centre of gravity. A sketch of the rod and the forces acting on it is depicted in Fig. 1.

The system is assumed inertialess, that implies the equilibrium of forces and torques. The first implies F H (pL) + F H (-pL) = 0, leading to v 0 = v G , that is, the rod centre of gravity moves with the fluid. Now, to prevent a resultant torque, force F H (pL) must align with p, i.e. F H (pL) = λp, λ ∈ R. Thus, we have

F H (pL) = ξ(∇v • pL -ṗL -p L) = λp, (1)
that multiplying by p and taking into account that p • p = 1 and consequently p • ṗ = 0, yields

ξ(∇v : (p ⊗ p)L -L) = λ, (2)
expression that introduced into Eq. 1 reads

ξ(∇v • pL -ṗL -p L) = ξ(∇v : (p ⊗ p)L -L)p, (3)
that leads to the rotary velocity ṗ ṗ = ∇v • p -∇v :

(p ⊗ p)p, (4)
that is nothing else than the standard Jeffery expression for ellipsoids of infinite aspect ratio (rods) [START_REF] Jeffery | The motion of ellipsoidal particles immersed in a viscous fluid[END_REF]. Now, by equating the force acting on the beads λ with the one within the spring, we have

2K(L -L 0 ) = ξ(∇v : (p ⊗ p)L -L), (5)
or

L = - 2K ξ (L -L 0 ) + ∇v : (p ⊗ p)L). ( 6 )
Thus, the kinematics of an elastic dumbbell of reference length 2L 0 with conformation at time t given by its orientation p and length 2L, read

ṗ = ∇v • p -∇v : (p ⊗ p)p L = -2K ξ (L -L 0 ) + ∇v : (p ⊗ p)L) . ( 7
)
From rigid ellipses to orthogonal elastic bi-dumbells

In [START_REF] Chinesta | From single-scale to two-scales kinetic theory descriptions of rods suspensions[END_REF] it was proven that in order to represent a rigid ellipse whose kinematics are given by the Jeffery equation it suffices to consider a rigid system composed of two mutually orthogonal rods whose lengths correspond with the length of the ellipse axes. In this case, if p refers to the direction of the ellipse largest axis, and F = r 2 -1 r 2 +1 , with r the ellipse aspect ratio, we have

ṗ = • p + F (D • p -D : (p ⊗ p)p) , (8)
where and D are respectively the vorticity and the rate of strain tensors, 2 = ∇v -(∇v) T and 2D = ∇v + (∇v) T . Now, we address the more general case in which both rigid rods are replaced by two extensible and mutually perpendicular springs of reference lengths L 0 1 and L 0 2 . In the sequel, the same reference lengths are assumed for both dumbbells, i.e. 2L 0 1 = 2L 0 2 = 2L 0 . In the present configuration, and considering that as proven in our former works the centre of gravity moves with the fluid, the hydrodynamic forces applying at beads L 1 p 1 and L 2 p 2 , F H 1 and F H 2 read respectively

F H 1 = ξ(∇v • p 1 L 1 -ṗ1 L 1 -p 1 L1 ), (9)
and

F H 2 = ξ(∇v • p 2 L 2 -ṗ2 L 2 -p 2 L2 ), (10)
with p 1 ⊥ p 2 , and with their orientation rates of change expressed from

ṗ1 = ω × p 1 ṗ2 = ω × p 2 . ( 11
)
The angular momentum balance implies now

L 2 1 p 1 × (∇v • p 1 -ṗ1 ) + L 2 2 p 2 × (∇v • p 2 -ṗ2 ) = 0, (12
) which coincides with the expression obtained in the case of rigid rods [START_REF] Chinesta | From single-scale to two-scales kinetic theory descriptions of rods suspensions[END_REF], proving the validity of the Jeffery equation in the case of orthogonal elastic bi-dumbells.

Introducing the Jeffery Eq. 8 with r

= L 1 L 2 that implies F = L 2 1 -L 2 2 L 2 1 +L 2 2 , i.e. ṗ1 = • p 1 + F D • p 1 -p T 1 • D • p 1 p 1 , (13)
into the expression of the hydrodynamic force acting on bead p 1 L 1 , we have

F H 1 = ξ(∇v • p 1 L 1 -ṗ1 L 1 -p 1 L1 ) (14) = ξ (∇v • p 1 L 1 -• p 1 L 1 -F D • p 1 L 1 -p T 1 • D • p 1 p 1 L 1 -p 1 L1 (15) = ξL 1 (1-F)D • p 1 +F p T 1 •D•p 1 p 1 -ξ p 1 L1 .( 16
)

Modelling incompressible ellipses from orthogonal elastic bi-dumbells

The projection of force F H 1 in the direction p 1 is the one that causes the spring extension, i.e.

2K(L 1 -L 0 ) -F I 1 = ξL 1 p T 1 • D • p 1 -ξ L1 , ( 17
)
where the force F I 1 ensures the incompressibility constraint. Similar calculations lead to

2K(L 2 -L 0 ) -F I 2 = ξL 2 p T 2 • D • p 2 -ξ L2 , ( 18
)
where it was assumed that

K 1 = K 2 = K.
Thus, the extension velocities read

L1 = - 2K ξ (L 1 -L 0 ) + L 1 p T 1 • D • p 1 + 1 ξ F I 1 , (19)
and

L2 = - 2K ξ (L 2 -L 0 ) + L 2 p T 2 • D • p 2 + 1 ξ F I 2 . ( 20
)
Taking into account

p T 1 • D • p 1 = D : (p 1 ⊗ p 1 ) p T 2 • D • p 2 = D : (p 2 ⊗ p 2 ) , ( 21
)
and the fact that, since p 1 and p 2 are mutually orthogonal,

(p 1 ⊗ p 1 ) + (p 2 ⊗ p 2 ) = I (with I the identity tensor), the incompressibility constraint d dt (L 1 L 2 ) = 0 reads L1 L 2 + L 1 L2 = - 2K ξ (L 1 -L 0 )L 2 - 2K ξ (L 2 -L 0 )L 1 +F I 1 L 2 ξ + F I 2 L 1 ξ = 0, (22)
which simplifies to

L1 L 2 + L 1 L2 = - 4K ξ L 1 L 2 + 2K ξ L 0 (L 1 + L 2 ) +F I 1 L 2 ξ + F I 2 L 1 ξ = 0. ( 23
)
Forces related to the incompressibility constraint are expected to contribute isotropically to the resulting macroscopic stress, and then taken in the flow incompressibility constraint. Thus, if we compute the contribution of F I 1 and F I 2 to the stress by using the Kramers rule, we have

σ I = F I 1 p 1 ⊗ L 1 p 1 + F I 2 p 2 ⊗ L 2 p 2 = F I 1 L 1 (p 1 ⊗ p 1 ) + F I 2 L 2 (p 2 ⊗ p 2 ), (24)
which suggests considering

F I 1 L 1 = F I 2 L 2 , since p 1 ⊗ p 1 + p 2 ⊗ p 2 = I.
Thus, considering

F I 2 = F I 1 L 1
L 2 in Eq. 23, we have

F I 1 L 2 1 + L 2 2 ξL 2 = 4K ξ L 1 L 2 - 2K ξ L 0 (L 1 + L 2 ), (25)
or

F I 1 = 4K L 1 L 2 2 L 2 1 + L 2 2 -2K L 0 (L 1 + L 2 )L 2 L 2 1 + L 2 2 , ( 26
)
that vanishes for the relaxed case

L 1 = L 2 = L 0 , i.e. F I 1 (L 1 = L 2 = L 0 ) = 0 and consequently F I 2 (L 1 = L 2 = L 0 ) = 0.
Thus, finally the governing equations for the orthogonal elastic bi-dumbbell representing an incompressible ellipse read:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ṗ1 = • p 1 + L 2 1 -L 2 2 L 2 1 +L 2 2 D • p 1 - L 2 1 -L 2 2 L 2 1 +L 2 2 p T 1 • D • p 1 p 1 L1 = -2K ξ (L 1 -L 0 ) -1 ξ F I 1 + L 1 p T 1 • D • p 1 L2 = -L 2 L 1 L1 F I 1 = 4K L 1 L 2 2 L 2 1 +L 2 2 -2K L 0 (L 1 +L 2 )L 2 L 2 1 +L 2 2 F I 2 = F I 1 L 1 L 2 . ( 27
)

Conformation descriptor

When considering control volumes (small enough with respect to the flow but large enough with respect to the cell size) in a flowing foam, it can be observed that the cells in each control volume have similar shapes and orientations.

In that case, the cell population in each volume element can be described using p 1 and L 1 (when L 0 is assumed known).

The contribution of a cell to the stress using the Kramers rule can be obtained from the elastic forces. The contribution due to the incompressibility constraint being isotropic, it can be aggregated to the pressure term. Thus, the conformation contribution σ c results

σ c = 2K(L 1 -L 0 )p 1 ⊗ L 1 p 1 +2K(L 2 -L 0 )p 2 ⊗ L 2 p 2 (28) = 2K(L 1 -L 0 )L 1 (p 1 ⊗ p 1 )+2K(L 2 -L 0 )L 2 (p 2 ⊗ p 2 ) (29) = 2K ( L 1 L 1 (p 1 ⊗ p 1 )+ L 2 L 2 (p 2 ⊗ p 2 )) , ( 30
)
that vanishes in the relaxed configuration

L 1 = L 2 = L 0 , with L 1 = 0 and L 2 = 0.
In these circumstances, the simplest choice for the conformation tensor consists of the second-order symmetric tensor c defined from

c = L 1 L 1 (p 1 ⊗ p 1 ) + L 2 L 2 (p 2 ⊗ p 2 ).
(31)

It is important to note that the conformation does not involve an averaging process affecting the different cells inside the control volume because we assumed that all of them share almost the same conformation. However, as it is usual in the modelling of suspensions, this hypothesis could easily be relaxed and we would thus consider the average of different cell conformation, as addressed in the micromacro simulations by Keunings and coworkers using the Lagrangian Particle Method -LPM - [START_REF] Halin | The Lagrangian particle method for macroscopic and micro-macro viscoelastic flow computations[END_REF][START_REF] Keunings | Micro-macro methods for the multiscale simulation viscoelastic flow using molecular models of kinetic theory[END_REF].

Even though many rheological behaviours could be associated to flowing foams, this paper focuses on the consideration of induced anisotropy due to the cell deformation. For this reason, and without loss of generality, we consider in what follows the simplest rheological behaviour consisting of a Newtonian behaviour complemented with an elastic contribution related to the cell deformation. In any case, more complex bulk rheologies could be considered, as for example the one related to a viscoplastic behaviour.

Thus, the macroscopic viscoelastic constitutive equation reads

σ = -pI + 2ηD + μc, ( 32
)
where p is the pressure, that can be viewed as the Lagrange multiplier associated with the macroscopic flow incompressibility constraints, η the effective homogenized fluid viscosity, and μ the rheological parameter affecting the contribution of the microscopic conformation that scales with K and the volume concentration of cells. The origin of the elasticity introduced into the conformation field evolution is related to the surface tension that resists the cell deformation from its spherical reference conformation. Topological changes can operate at the elemental cell level. However, on average their net effect is to avoid too large cellular distorsions. Thus, the coefficient μ should describe both the surface tension as well as the microscopic topological changes. In the proposed model, it remains purely phenomenological and should be identified from appropriate rheological tests.

Even though the present model does not address relaxation mechanisms, they could easily be incorporated by including a viscous component at the spring-beads level.

Rheological behaviour

In order to study the rheological response of the proposed model, we carried out a numerical study to obtain the loss and storage moduli of a foam. To proceed, consider a periodic shear strain

γ (t) = γ 0 sin(ωt) (33)
whose time derivative provides the evolution of the shear rate

γ (t) = γ 0 ω cos(ωt). ( 34
)
In simple shear flow, the velocity gradient thus reads ∇v = 0 γ (t) 0 0 , and we could use the model Eq. 27 to obtain the time evolution of the shape and orientation of a cell subject to such an oscillatory flow. We can thus obtain the time evolution c(t) of the conformation tensor introduced in the previous section.

We can now decompose the off-diagonal component of tensor c(t), that is c 12 (t) in its in-phase and out-of-phase contributions according to

c 12 (t) = c 12 sin(ωt) + c 12 cos(ωt), ( 35
)
where c 12 and c 12 are analogous to the so-called storage (G ) and loss (G ) modulus respectively (they actually correspond to the storage and loss moduli associated with the stress contribution σ c ).

Figure 2 shows the value of c 12 and c 12 for an oscillation angular frequency ranging from 0.1 to 100 rad/s. The parameters are ξ = 0.1, K = 1, L 0 = 0.005 (metric system of units). We can recognize a Maxwell-like model, with the classical slopes of 1 (G ) and 2 (G ) at low frequency.

To obtain properly the storage (G ) and loss (G ) moduli, we should also consider the viscous component of the stress tensor. Thus the moduli are given by

G = σ 12 γ 0 = μ c 12 γ 0 (36) G = σ 12 γ 0 = 2ηω + μ c 12 γ 0 . ( 37
)
Finally, we also conducted another rheological test to observe the behaviour of the conformation tensor when a constant shear flow is brutally stopped. As expected, we obtain an exponential relaxation of the tensor components towards zero.

Macroscopic flow problem

Neglecting inertia, the macroscopic flow problem is defined as follows:

⎧ ⎨ ⎩ ∇ • σ = 0 σ = -pI + 2ηD + μc ∇ • v = 0 , ( 38
)
where σ is the Cauchy's stress tensor, I the unit tensor and v the macroscopic fluid velocity field. It can be noticed that the problem reduces to the standard Stokes flow model as soon as μ = 0. Incompressibility is justified by the fact that drainage is neglected.

The flow model must be complemented with the conformation expression whose evolution is governed by the microstructural model

c = L 1 L 1 (p 1 ⊗ p 1 ) + L 2 L 2 (p 2 ⊗ p 2 ), (39)
⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ṗ1 = • p 1 + L 2 1 -L 2 2 L 2 1 +L 2 2 D • p 1 - L 2 1 -L 2 2 L 2 1 +L 2 2 p T 1 • D • p 1 p 1 L1 = -2K ξ (L 1 -L 0 ) -1 ξ F I + L 1 p T 1 • D • p 1 L2 = -L 2 L 1 L1 F I 1 = 4K L 1 L 2 2 L 2 1 +L 2 2 -2K L 0 (L 1 +L 2 )L 2 L 2 1 +L 2 2 F I 2 = F I 1 L 1 L 2 . ( 40
)
In order to solve the resulting flow model defined in ⊂ R 2 , appropriate boundary conditions must be enforced at the domain boundary

≡ ∂ v = v g in D t = σ • n = t g in N , ( 41
)
with D ∪ N = and D ∩ N = ∅.

Numerical solution strategy

Cells are assumed represented at the initial time by a spherical conformation tensor c = 0, since the cells are initially in the relaxed state,

L 1 = L 2 = L 0 .
1. The associated Stokes problem is solved at time t n , from the microstructural term given at the previous time step μc(t n-1 ). The problem is solved using a standard mixed velocity-pressure formulation using any finite element satisfying the stability conditions, the so-called LBB conditions. In the example reported later, we considered Q9/Q4 finite elements; 2. Then, from the just computed velocity field v(t n )

the orientation and extension fields, p 1 (x, t n ) and L 1 (x, t n ) respectively, are updated by using a first-order discontinuous Galerkin formulation; 3. Finally, the conformation tensor is updated c(t n ).

These three steps are repeated until reaching the maximum simulation time or the steady state.

Updating the conformation field

In the numerical experiments described below, we consider 2D flows defined in 2D geometries ⊂ R 2 . In the 2D case, the unit vector p 1 can be expressed by p 1 = (cos θ, sin θ) T . By taking the time derivative, we get ṗ1 = θ(sin θ, cos θ) T . Using θ as orientation descriptor, the Jeffery equation can be expressed as θ(sin θ, cos θ) T = G(θ, ∇v, L 1 ).

(

) 42
Multiplying the previous expression by (sin θ, cos θ), we obtain the scalar equation

θ = G(θ, ∇v, L 1 ), ( 43
)
where the material derivative can be expressed in an Eulerian framework by introducing the orientation field θ(x, t) whose evolution is governed by

∂θ ∂t + v • ∇θ = G(θ, ∇v, L 1 ). ( 44
)
The scalar equation governing the evolution of the spring length can also be written as

∂L 1 ∂t + v • ∇L 1 = H(θ, ∇v, L 1 ). ( 45
)
Both Eqs. 44 and 45 are purely advective, and thus appropriate discretization taking into account their hyperbolic character must be used. We make here the simplest choice, a first-order discontinuous Galerkin scheme, that considers the generic variable P (θ or L 1 ) constant in each element e belonging to the mesh M of . The test function is assumed also constant in each element, vanishing outside. The balance in element e reads, taking into account the flow incompressibility (∇

• v = 0): e ∂P ∂t dx + e ∇ • (vP)dx = e J dx, ( 46
)
where the source term J represents G or H depending on the considered equation.

Using the divergence theorem, the second term of the left-hand side can be written from the boundary flux, i.e.

+ ∂ e Pv • n dx = e J dx, ( 47
)
where n is the unit outward vector normal to the element boundary. As P is not defined on ∂ e , we consider the element boundary decomposition ∂ e = ∂ + e ∪ ∂ -e , where ∂ -e and ∂ + e represents the inflow and outflow element boundaries, both defined from v•n < 0 and v•n > 0 respectively. Then we assume that the property on the inflow boundary is given by its value at the upstream element, i.e. P(x ∈ ∂ -e ) = P e -, whereas on the outflow element boundary, it is given by the property at element e , i.e. P(x ∈ ∂ + e ) = P e .

Discussion

-In the flow model just proposed there is no size effects, which implies that the characteristic length of cells must be small in relation to the one characterizing the spatial variation of the macroscopic velocity field. -The mesh considered for integrating the velocity field must be small enough to capture all the macroscopic velocity field details. -The mesh considered for calculating the conformation field must be small enough for assuming cells in each element of the mesh described by the conformation field inside the element, and it must be large enough for assuming that it represents a population of cells. However, as soon as the model is described in a continuous way, there is no contradiction with the fact of considering elements smaller that the characteristic size of the cell (in the case of polymer flows one can consider elements smaller that the size of a molecule).

Numerical results

Uncoupled microstructure-flow calculations

In this section, we address some simple flows in order to evaluate the response of the cells. The model is uncoupled in the sense that the flow induces cell deformation, but the cell conformation does not affect the flow kinematics. Thus, even though some of the flows addressed here could exhibit rich kinematics in practice, in what follows the uncoupled solution does not allow to capture such a rich kinematics, as for example shear banding. Moreover, the only interest of this section being the evolution of the cell deformation and its macroscopic description, rheology is not considered in the analyses carried out.

In the numerical examples considered in this section a unit square is considered as depicted in Fig. 3. The initial (relaxed) conformation is given by a zero conformation tensor c(x, t = 0) = 0.

Four different flows are considered: (i) a simple shear; (ii) a contraction flow; (iii) the driven cavity flow problem and finally (iv) the flow around a square obstacle. In all cases, the microstructure was computed on the basis of the associated Stokes kinematics, i.e. considering the uncoupled flow associated with μ = 0.

In order to quantify the way the flow behaviour (shear, rigid motion or elongation) affects the microstructure evolution we consider a flow criterion. For that purpose, first, we introduce the relative rate of rotation W from

W = ω -w, (48)
where ω is related to the flow vorticity according to

= • ω, ( 49
)
where is the third-order permutation tensor (also known as the Levi-Civita tensor). The vector ω can also be written in terms of the curl of the velocity, (52)

ω = - 1 2 ∇ × v. (50)
Figure 4 shows the evolution of the microstructure conformation for different K in a pure shear with ξ = 0.1 and L 0 = 0.005 (metric system of units). The qualitative analysis performed here does not require a precise determination of these coefficients. As it can be noticed, the more stiffness the spring has, the less deformation the conformation presents. In order to prove that the conformation is not only accommodating a deformation but also that it is rotating, we included a small cercle on the ellipsoid surface to appreciate the way in which that point is evolving in time.

Simple shear flow

We consider V = 2 m • s -1 on the N-boundary, zero velocity on the S-boundary and a linear velocity evolution on the E-boundary and W-boundaries. The microstructure effects were introduced by using K = 1.

Figure 5 depicts the velocity field, that evolves almost linearly through the domain thickness, as well as the flow criterion, that as expected corresponds to a perfect shear behaviour. The pressure field is constant in the whole flow domain.

Figure 6 shows the steady state conformation. Because of the shear, the conformation is expected to rotate clockwise. The microstructure shows a significant variation along the domain as a consequence of the different velocity and constant shear rate.

Extrusion-like flow

In the present case, we consider again K = 1. A Poiseuille velocity profile (parabolic) is enforced on the W-boundary (-2y(y -1), 0) T , whereas at the Eboundary the fluid leaves the square domain throughout an exit where tension-free boundary conditions are enforced.

Figure 7 depicts both components of the velocity field whereas Fig. 8 shows the steady state conformation and the flow criterion. As it can be seen, the largest conformation axis remains aligned with the streamlines. A pure extension is noticed along the symmetry axis in agreement with the flow criterion.

Driven cavity flow

The only difference with respect to the previously analyzed flow is that now a unit horizontal velocity is applied on the top wall (N-boundary) that induces the fluid flow within the cavity.

Figure 9 depicts both components of the velocity field. Figure 10 shows the conformation at a given time as

Flow around an obstacle

In the present case, a unit horizontal velocity is enforced on the W-boundary, with non-slip conditions on the top and bottom walls and a free boundary condition is specified on the E-boundary. The domain contains in its center an obstacle where non-slip boundary conditions are enforced.

Figure 11 depicts both components of the velocity field. The velocity field presents two stagnation points located at the intersection between y = 0.5 and the obstacle. These stagnation points correspond with the maximum and minimum pressures.

Figure 12 depicts the conformation and flow criteria. As it can be noticed, elongation is specially present upstream and downstream, the shear being located in the neighbourhood of the upper and lower sides of the obstacle, with the expected effects on the microstructure conformation.

Coupled simulations

This section adresses a coupled simulation in which flow kinematics induces microstructure evolution and the latter affects at its turn the flow kinematics.

Figure 13 shows two snapshots taken from films of flowing foams where the flow and microstructures were assumed almost at steady-state. The microstructure of these snapshots will be used first to identify the model parameters and then to test the agreement between the computed and observed microstructures. Of course, from the information of kinematics only, we cannot conclude on rheological aspects but at least this serves to evaluate the proposed model in terms of conformation evolution.

In order to extract from these images the conformation tensor, both images were segmented (every pixel is transformed into either black or white) and then cell boundaries were easily identified. From that, the centre of gravity of each cell can easily be obtained as well as the conformation tensor. If X e i , i = 1, . . . , Q e are the points defining the wall of cell C e , the centre of gravity and the inertia tensor J e are obtained from

X e G = 1 Q e Q e i=1
X e i (53)

and

J e = 1 Q e Q e i=1
(X e i -X e G ) ⊗ (X e i -X e G ) (54) respectively. The numerical inertia tensor can be calculated as soon as p 1 and L 1 are known (p 2 and L 2 derive from p 1 and L 1 ).

J = L 2 1 (p 1 ⊗ p 1 ) + L 2 2 (p 2 ⊗ p 2 ). ( 55
)
In Fig. 14, the experimental ellipses related to J e are superimposed to the cells, and also to the solution predicted by using the proposed model. The model parameters K, ξ , η and μ are, as previously indicated, chosen in order to obtain cell shapes close to those observed experimentally: (metric system of units) K = 1, ξ = 0.1, η = 1 and μ = 0.1. L 0 = 0.005 was identified from the average cell surface observed experimentally. In absence of rheological data, we noticed that many choices of those parameters lead to similar microstructures and that the impact of those on the kinematics was almost negligible. Thus, a proper rheological characterization seems compulsory to obtain an adequate flow-microstructure coupling.

Conclusions

This work proposes a simple model for flowing foams, where the microstructure is introduced from a conformation tensor that describes the structural elasticity. The numerical predictions have been compared with some experimental results proving the ability of the model to describe the effective kinematics as well as the flow induced microstructure evolution.

More applicative analyses, in particular in the study of industrial processes, require an appropriate rheological characterization and very probably the proposal of appropriate rheometric devices. All these aspects will be addressed in future works.

  In the sequel, we consider the following tensor products, where Einstein's summation convention is assumed: -if a and b are first-order tensors, the single contraction • reads (a • b) = a j b j ; -if a and b are first-order tensors, the dyadic product ⊗ reads (a ⊗ b) jk = a j b k ; -if a and b are respectively second and first-order tensors, the single contraction • reads (a•b) j = a jm b m ; -if a and b are second-order tensors, the single contraction • reads (a • b) jk = a jm b mk ; -if a and b are second-order tensors, the double contraction : reads (a : b) = a jk b kj .

Fig. 1

 1 Fig. 1 Extensible rod immersed in a flow

Fig. 2

 2 Fig. 2 Oscillatory shear -c 12 and c 12 components

Fig. 3 Fig. 4

 34 Fig. 3 2D domain

Fig. 5 Fig. 6

 56 Fig. 5 Velocity field (left) and flow criterion (right) in a simple shear flow

Fig. 7 Fig. 8 Fig. 9 Fig. 11

 78911 Fig. 7 Extrusion-like flow problem. Velocity field: (left) x-component and (right) y-component

Fig. 13

 13 Fig. 13 Experimental snapshots (courtesy of F. Graner)

Acknowledgements This project has received funding from the European Unions Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No. 675919.

Compliance with Ethical Standards

Conflict of interests

The authors declare that they have no conflict of interest.

Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


Ezzahra Fatma

Salem

Tijani Chahed

Zwi Altman

Azeddine Gati
email: azeddine.gati@orange.comtijani.chahed@telecom-sudparis.eu

Traffic-aware Advanced Sleep Modes management in 5G networks

Keywords: Advanced Sleep Modes, 5G networks, Qlearning, traffic load, latency

Advanced Sleep Modes (ASMs) are defined as a progressive shutdown of the Base Station (BS) depending on the activation and the deactivation times of the different components. This transition duration defines different levels of sleep modes that can be implemented in future 5G networks. We propose in this paper a management strategy based on Q-learning approach which will enable to find the best combination and durations of ASM levels depending on the traffic load and the network operator's policy regarding energy reduction versus latency. Our results show that even in delay-sensitive scenarios, high energy gains can be achieved in low and moderate traffic loads, respectively 55% and 10%, without inducing an extra latency. Starting from a certain traffic load (approximately 30%), ASM should not be implemented in case of stringent latency constraint.

I. INTRODUCTION

The multiplicity of mobile devices and the upsurge of new services have led to an exponential growth of traffic demand. This has induced a huge increase in the energy consumed by the network operators [START_REF] Andrews | What Will 5G Be?[END_REF] which has an impact on both the operator's Operational Expenditures (OPEX) and the environment due to the increasing CO 2 emissions [START_REF] Fehske | The global footprint of mobile communications: The ecological and economic perspective[END_REF]. To cope with these problems, considerable efforts have been invested to reduce the energy consumption of mobile networks and more precisely of the Base Stations (BSs) as they are the most energy consumers (around 80% of the overall mobile energy consumption) [START_REF] Frenger | Reducing Energy Consumption in LTE with Cell DTX[END_REF]. Several techniques and technological levers have been investigated in this area such as Massive Multiple Input Multiple Output (MIMO) systems [START_REF] Salem | Energy consumption optimization in 5G networks using multilevel beamforming and large scale antenna systems[END_REF], cell zooming [START_REF] Niu | Cell zooming for cost-efficient green cellular networks[END_REF] and sleep modes [START_REF] Wu | Energy-Efficient Base-Stations Sleep-Mode Techniques in Green Cellular Networks: A Survey[END_REF].

We focus in this work on the feature of Advanced Sleep Modes (ASMs) which corresponds to a progressive deactivation of the BS's components [START_REF] Debaillie | A Flexible and Future-Proof Power Model for Cellular Base Stations[END_REF]. Different levels are considered in ASM according to the time needed for the activation and the deactivation (transition time) of each component. Going from one level to another, more components are deactivated which reduces the BS's power consumption, but this is at the expense of a longer delay as more time is needed to reactivate the BS again.

The sleep modes are already present in LTE systems under the name of µDTX (≈ 71 µs) [START_REF] Frenger | Reducing Energy Consumption in LTE with Cell DTX[END_REF]. With the lean carrier design in 5G networks permitting to have less control signals, deeper sleep modes can be implemented. It has been shown in previous works [START_REF] Salem | Advanced Sleep Modes and Their Impact on Flow-Level Performance of 5G Networks[END_REF] that ASMs can reduce up to 90% of the energy consumption in very low load. However, the delay increase may be up to 5ms which is critical for some 5G use cases like Ultra-Reliable Low Latency Communications (URLLC) where the maximum acceptable latency is around 1ms [START_REF]5G White Paper[END_REF]. For this reason, a delay-sensitive orchestration of the different ASM levels is needed in order to know how much each of them is allowed to be used depending on the length of the idle periods of a BS and according to the network operator's needs with respect to latency. A first study was presented in [START_REF] Salem | Reinforcement learning approach for Advanced Sleep Modes management in 5G networks[END_REF] where a Q-learning based implementation strategy has been applied for the case of very low loads. Hence, a more comprehensive investigation is needed in order to derive the optimal policies for higher loads values.

We present in this paper a traffic-aware implementation strategy which enables to have a codebook mapping the traffic load to the possible actions to be performed in the Q-learning algorithm. According to the constraint of the network operator imposed on the latency and to the traffic conditions, especially regarding the average inter-arrival times, several policies can be deduced.

The remainder of this paper is as follows: Section II-A describes the concept of ASMs, the 5G new radio opportunities that can be exploited to have deeper sleep modes, the buffering strategy and the proposed solution for ASM implementation. Section III develops our ASM management solution based on Q-learning approach. Section IV presents the simulation results and Section V concludes the paper and gives some perspectives for future works.

II. SYSTEM DESCRIPTION A. Advanced Sleep Modes description

The different components and sub-components of a BS can be classified into four categories according to the time needed by each of them to deactivate and then wake up again (transition time) as presented in [START_REF] Debaillie | A Flexible and Future-Proof Power Model for Cellular Base Stations[END_REF]. These categories define the different possible levels of sleep modes that a BS can enter, and are denoted as "Advanced Sleep Modes":

• SM 1 : represents the shortest sleep mode which needs a transition time of 71µs (OFDM symbol). The power amplifier and some components of the digital baseband and the analog front-end (both in Rx and Tx) are deactivated.

• SM 2 : corresponds to a longer sleep mode which needs 1ms as a transition time (1 sub-frame or TTI) and in which more components of the analog front-end are disabled compared to the first level. • SM 3 : the power amplifier, all the components of the digital baseband and almost all the parts of the analog front-end (except the clock generator) are switched off. This subset of components needs 10 ms (a frame) to deactivate then wake up again. • SM 4 : corresponds to the standby mode in which the BS retains its wake-up functionalities but almost all its components are deactivated. It groups the components that need more than 1s for the transition latency (deep sleep mode).

The ASMs are defined only based on the transition times of the BS's components. However, the sleep length of each of them has not been defined. We assume that the minimal sleep duration in each sleep level is in the same time scale of the transition time corresponding to that level, and that both the activation and deactivation delays are half of this transition time. This is summarized in Table I. With the lean carrier design of 5G new radio, synchronization signals are grouped in blocks and are sent periodically with a configurable periodicity among this set [5, 10, 20, 40, 80, 160 ms] as agreed in 3GPP [START_REF]Radio Resource Control (RRC) protocol specification (Release 15)[END_REF]. This allows to enter deeper sleep mode levels i.e., SM2 and SM3. For this reason, we consider in this work the first three ASM levels.

C. Buffering strategy

We consider in this work only users in connected mode and ignore users' initial access. When a user arrives during a sleep mode, his arrival triggers the activation of the BS. In the meantime, the user request is kept in a buffer until the BS wakes up and serves him. The buffereing delay depends on the depth of the sleep mode and on its phase, i.e., whether the BS is deactivating, activating or it is in the sleep mode itself. This is illustrated in Figure 1.

The deeper the sleep mode, the more power reduction we can have as we deactivate more hardware. This hampers the responsiveness of the BS due to the buffering time for the users arriving during the sleep mode. This may be a problem in some scenarios where there is a high constraint imposed on the delay such as URLLC type of service. Hence, a delaysensitive implementation strategy is needed in order to be able to profit from the ASM without inducing any impact on the latency.

D. Delay-sensitive implementation strategy

We proposed in [START_REF] Salem | Reinforcement learning approach for Advanced Sleep Modes management in 5G networks[END_REF] an implementation strategy for ASM which consists of putting the BS in different levels of sleep modes as follows: each time the BS is idle, we start with the deepest level (SM3), then we wake up the BS gradually: we put it into SM2 and then SM1 then in the idle mode as shown in Figure 2. Each sleep level can be repeated many times before going to the next one. The target is to determine the duration of each sleep mode.

The idea is to have different possibilities of sleep durations in each sleep mode level according to the policy of the network operator. If this policy consists of prioritizing the energy consumption reduction, then the best strategy is to remain in the deepest sleep level until the arrival of the next user which will trigger the activation of the BS. However, if the network operator decides to prioritize the users' latency, then the best strategy to follow is to anticipate the arrival of a user by taking more conservative actions (lower durations of sleep modes) so that whenever a user arrives, he finds the BS ready to serve him. In the case where both energy reduction and latency are to be considered, we make use of a utility function which is equal to the weighted sum of both criteria. We try to maximize it following the same strategy i.e., deepest sleep level first. We make so based on a Q-learning algorithm that depends on the traffic load, as shown next. As in [START_REF] Salem | Reinforcement learning approach for Advanced Sleep Modes management in 5G networks[END_REF], we consider the Q-learning approach [START_REF] Watkins | Q-learning[END_REF] in order to find the optimal duration for each sleep mode level. The advantage of this method is that it is model-free which makes it well adapted for the ASM problem.

A. Q-learning algorithm

The Q-learning approach is based on the estimation of a state-value function denoted by Q which links the state s to the action a and is updated in each iteration.

Table II presents the pseudo-code summarizing the Qlearning algorithm:

Q-learning algorithm Initialize Q(s, a) arbitrarily Repeat (for each episode): Initialize s Repeat (for each step of episode): Choose a from s using policy derived from Q (e.g., -greedy) Take action a, observe R, s • an episode is defined according to the users' dynamics.

Q(s, a) ← Q(s, a) + α[R + γ max a Q(s , a ) -Q(s, a) ] until s is terminal
The departure of a user corresponds to the begining of an episode, if there is any other user being served. • s denotes the current state of the system and s the following state after performing action a. • R is the reward achieved by choosing action a.

• α is the learning rate.

• γ is the discount factor describing the weight given to future rewards. • The -greedy algorithm operates in such a way to test random actions with probability (exploration) and to take the actions which maximize the Q value with a probability 1 -(exploitation).

B. Q-learning characteristics 1) State space: It corresponds to the different states of the BS. It can be either active (serving a user), idle (not transmitting anything but it is still activated) or in sleep mode (SM 1 , SM 2 or SM 3 ). S = {active, idle, SM 1 , SM 2 , SM 3 }.

2) Decision epochs: The decision points are defined at the end of each state: if the BS was in the active state, the departure of a user marks the end of this state and an action should be performed. This is the start of an episode. After peforming that action, the BS will be in the next state, first in SM3 if applicable, or else SM2, or else SM1. Once the time allowed for the corresponding sleep level elapses, another decision should be taken, it consists of going to a higher SM level. This is repeated until the arrival of a user which will end the current episode and impose the activation of the BS or after having entered all the possible sleep mode levels then activated the BS.

3) Action space: The action space corresponds to the decisions that can be taken at each desicion point. An action consists of finding n i the number of times the BS can repeat the following sleep mode level i. Let N i be the set of possible values that n i can take, for i ∈ {1, 2, 3}.

So, the action space is:

A s = {N 3 × N 2 × N 1 }.
The objective is to find the optimal triplet (n 3 , n 2 , n 1 ) which maximize a reward function R.

4) Reward function:

The reward R for a given episode can be defined as the weighted sum of the energy gain E gain and the added delay D, that are achieved following the sleep strategy during that episode.

Hence, R can be written as follows:

R = -D + (1 -)E gain (1)
where is a normalized weight ( ∈ [0, 1]) that denotes the importance given to the two factors D and E gain . A small means that the energy reduction is prioritized over the delay and vice versa.

C. Traffic-aware action space

The ASMs are intended to be used during the idle periods of the BS. They are most profitable when the traffic load is very low as we have enough time between the users to put the BS in deep sleep mode for large periods. It has been shown that in almost zero-load condition, we can achieve up to 90% of energy reduction [START_REF] Salem | Advanced Sleep Modes and Their Impact on Flow-Level Performance of 5G Networks[END_REF]. When the load increases, the users' inter-arrival time decreases which will restrict the use of sleep modes. The action space denoting the number of times to repeat each sleep mode level is a finite set of some possible values that should be well adapted to the traffic load.

An action n i can be written as follows:

n i = x i n (2)
where:

• x i is the number of times to repeat the sleep mode i in order to cover all the idle period only by this level. It can be defined as x i = I Ti with I is the average inter-arrival time of users and T i the duration of sleep mode i.

• n is a multiplicative factor taken in such a way to have different possibilities of sleep duration for each sleep mode level while being around the optimal value x i . We consider n ∈ [0, 0.1, 0.2, 0.5, 1, 2, 5, 7, 10] in our numerical experiments shown next.

IV. NUMERICAL RESULTS

A. Simulator description

We consider a 2x2 MIMO BS with a single sector and we simulate an FTP service where the users request to download a file of exponential size and once the download is complete they quit the network. The simulator is event-based where an event corresponds to the dynamics of users and sleep modes (activation, deactivation, change of sleep level). The different characteristics of the simulator are given in Table III. where Ns,a is the number of visits to the action-state pair (s, a). Discount factor γ = 0.1 Exploration factor exp = 0.1

In order to quantify the outcome of the proposed solution, we focus on two performance metrics: energy consumption and latency.

The energy consumption can be defined as follows:

EC = Ne i=1 τ i-1,i P i-1,i (3)
where:

• N e is the number of events.

• τ i-1,i is the duration between two events i -1 and i.

• P i-1,i is the power consumed by the BS during τ i-1,i and is computed using IMEC power model tool [START_REF]IMEC Power Model Tool[END_REF]. This power value depends on the state of the BS during that period as mentioned in Table IV. The latency of a user can be computed as follows:

∆ Latency = T serve -T req (4)
where T serve is the time when the BS starts serving the user and T req is the time when the request occurs. We focus here on the latency induced by the sleep modes due to the buffering approach.

We consider different values of arrival rates λ ∈{1, 3, 6, 10, 20, 35} (users/s/Km 2 ) which are translated into loads of 1%, 4%, 9%, 15%, 30% and 55%, respectively. The load is defined as the fraction of time during which the BS was active:

ρ = 1 T Ne i=1 τ i-1,i 1 {xe(τi-1,i)=1} (5)
where :

• T = Ne i=1 τ i-1,i , is the total observation time.

• x e (τ i-1,i ) = 1 if the BS was active during τ i-1,i . 0 otherwise.

B. Convergence analysis

The convergence criteria is based on the computation of the maximum variation of the quality-value Q(s, a) in each iteration for all the state-action pairs (s, a). Figure 3 shows the convergence behaviour in a low load scenario (1%) with = 0.75 during the learning phase. The maximum variation tends to zero for a sufficiently high number of iterations (in the order of 10 4 ).

C. Performances results

We consider five different values of : 0.001, 0.25, 0.5, 0.75 and 0.999. For ≈ 0, the energy consumption reduction is prioritized, whereas when ≈ 1, it is the delay that counts the most. In the cases in between, different weights are given to the two factors. 1) Output policies: Table V presents the differents policies, i.e., the triplet (n 3 , n 2 , n 1 ), that are deduced depending on the value of for the different considered arrival rates. The more we increase the more we reduce the time spent in SM 3 giving more opportunities for the next levels.

2) Energy consumption reduction with ASM: Figure 4 shows the energy consumption reduction achieved for different loads and different values of . The highest energy gain (90%) is achieved in low load with the lowest . The more we increase the load, the more we reduce the energy gain as we use less sleep modes between the users' arrivals. For a given value of load, the energy gains decrease with , as when we increase , we give more weight to the delay so that we use the sleep modes for shorter periods to prevent the arrival of users during the sleep state. Even with this constraint, we can achieve up to 55% of energy reduction in low loads. For average loads, the energy gain is also significant (10% of reduction for 10% of load). Starting from 30% of the load, the best policy to take in order to satisfy the delay constraint is to not use any sleep mode level. 3) Influence on users' latency: Figure 5 shows the variation of the added delay caused by the sleep strategy for different values of load and . We observe that if we increase , we can reduce this impact on latency until having no (or very little) influence (lowest curve).

V. CONCLUSION

We presented in this paper a traffic-aware, energy and delaysensitive implementation strategy for ASM in 5G networks. It takes into account the traffic variation and the constraint that may be imposed by the network operator on energy and latency in order to derive the policy that captures the projected trade-off between them. Different policies can be achieved based on different traffic loads as well as the network operator's policy regarding delay and energy reduction. Our results show that significant energy gains can be achieved by the ASM, especially in low load scenarios and when there is no constraint imposed on the delay: up to 90% of energy reduction. Even with a high constraint on the delay, This study was performed without considering the control signals. As an extension of this work, we aim to add the periodic signaling and find out the possible changes in the policies that we can deduce according to the traffic conditions. It is also interesting to extend the work to a network with multiple BSs and consider cooperation between them in order to have the highest possible energy gains with the minimal impact on the delays. A centralized orchestration of the different BS sleeping cycles is thus needed.

Fig. 1 :

 1 Fig. 1: Buffering delays in ASMs

𝒏Fig. 2 :

 2 Fig. 2: Delay-sensitive implementation strategy of ASMs

Fig. 3 :

 3 Fig. 3: Convergence evolution

Fig. 4 :

 4 Fig. 4: Energy consumption reduction

Fig. 5 :

 5 Fig. 5: Impact on latency

TABLE I :

 I Advanced Sleep Modes characteristics PBCH) which is transmitted in the first sub-frame of the LTE radio frame. With such frequent signaling, only the first level of sleep modes (SM1) is applicable.

	Sleep Deactivation	Minimum	Activation
	level	duration	sleep duration	duration
	SM1	35.5 µs	71 µs	35.5 µs
	SM2	0.5 ms	1 ms	0.5 ms
	SM3	5 ms	10 ms	5 ms
	SM4	0.5 s	1 s	0.5 s
	B. Signaling in 5G new radio
	In LTE systems, the BS sends periodically control signals
	such as Cell-specific Reference Signals (CRS) broadcasted
	during 4 OFDM symbols out of 14, Primary and Secondary
	Synchronization Signals (PSS and SSS) transmitted every 5ms
	and Physical Broadcast CHannel (

TABLE II

 II

	: Q-learning algorithm
	where :

TABLE III :

 III Simulator characteristics

	Network parameters
	Antenna height	30 m
	Bandwidth	20MHz
	Scheduling type	Round Robin
	Channel characteristics
	Thermal noise	-174 dBm/Hz
	Path loss (d in km)	128.1 + 37.6 log 10 (d) dB
	Shadowing	Log-normal (6dB)
	Traffic characteristics
	Users' arrivals	Log-normal with mean λ
		and variance v = λ 10
	Service type	FTP
	Average file size	4 Mbits
	Q-learning parameters
	Learning rate	α = 1 Ns,a

TABLE IV

 IV

	: Power consumption of 2x2 MIMO BS (1 sector)
		in different states.
	Radiated Power: 46 dBm, Bandwidth: 20 MHz
	Active	Idle	SM 1	SM 2	SM 3
	250 W	109 W	52.3 W	14.3 W	9.51 W

TABLE V :

 V Policies

	Arrival rate		≈ 0	= 0.25	= 0.5	= 0.75	≈ 1
	(users/s/Km 2 )
		n 3	6944	1389	694	694	694
	λ 1 = 1	n 2	69444	69444	13889	6944	0
		n 1 1369327 391236	978091	391236 39124
					1369327 978091
					1956182 1369327
						1956182
		n 3	4630	463	231	231	46
	λ 2 = 3	n 2	9259	46296	4630	2315	0
			23148
			32407
			46296
		n 1	130412	652061	130412	326030 13041
			326030		326030
			456442		456442
			652061		652061
		n 3	1157	139	116	116	0
	λ 3 = 6	n 2	4630	13889	4630	1157	0
			11574		16204
		n 1	65206	136933	163015	65206	6521
			163015		228221
			228221		326030
			326030
		n 3	1389	69	69	69	0
	λ 4 = 10	n 2	13889	6944	2778	694	0
		n 1	195618	68466	97809	136933	3912
				97809	136933	195618
					195618
		n 3	694	40	35	14	0
	λ 5 = 20	n 2	6944	1984	3472	694	0
		n 1	97809	39124	68466	97809	0
				55891	97809
		n 3	397	40	8	4	0
	λ 6 = 35	n 2	3968	1984	1984	394	0
		n 1	55891	39124	39124	55891	0
				55891