Giulio Sciarra

Francesco Dell'isola

Kolumban Hutter

Dilatational and Compacting Behavior around a Cylindrical Cavern Leached Out in a Solid-Fluid Elastic Rock Salt

Keywords: Salt domes, Dilatancy, Compaction, Mixtures, Prestressing, Elastic analysis, Pressures

A fluid-filled cylindrical cavern of circular cross section in a homogeneous infinite fluid-saturated polycristalline ͑salt͒ formation subjected to isotropic stress is set under internal pressure that differs from the confining pressure. The fluid in the cavern and in the mixture is treated as ideal and the solid as elastic. The state of stress that is established as a consequence of an outside pressure and a cavern pressure serves as the reference state. Perturbing the cavern pressure induces small changes in the solid and fluid densities and in the solid displacements. We compute these and other fields as functions of the radial distance from the cavern center and show that, depending on the relative stress levels, the ͑salt͒ formation experiences either a dilatation or a compaction that is highly concentrated in a thin boundary layer near the cavern wall and tapers off as one moves away from it. The amount of dilatation/compaction of the cylindrical wall and the thickness of the boundary layer grow with an increase in the difference between the referential confining pressure and the pressure in the cavern.

Introduction

In geotechnique, caverns excavated in salt domes are used as storage volumes for the disposal of waste and hydrocarbons. These caverns, leached out in salt formations, are characterized by a damaged zone in the vicinity of their surface. In fact, permeability measurements ͑Cosenza et [START_REF] Cosenza | In situ rock salt permeability measurement for long term safety assessment of storage[END_REF][START_REF] Stormont | In situ gas permeability measurements to delineate damage in rock salt[END_REF] in the immediate vicinity of the cavern walls disclose a dilatation of the pores within a relatively small boundary layer of the polycrystalline salt. The origin of this disturbed rock zone is still somewhat debatable; however, a microstructured theory of solid-fluid mixtures was recently proposed ͑Sciarra 2002; [START_REF] Sciarra | A solid-fluid mixture model allowing for solid dilatation under external pressure[END_REF][START_REF] Sciarra | A variational approach to a micro-structured theory of solid-fluid mixtures[END_REF]dell'Isola et al. 2000͒ that offered a means of explaining, in the framework of an elastic model, the occurrence of the apparently counterintuitive expansion of the pore-space by the increase in internal cavern pressure. In our initial publication on this subject ͑Sciarra et al. 2001͒, second-gradient effects for the deformed porous salt matrix, which was filled with an ideal fluid, were claimed to be essential. Our present understanding tells us that second-gradient effects yield dominant qualitative behavior of the dilatancy deformation of the salt-fluid mixture only for specific simple load configurations, while in general it affects the deformation field only quantitatively, for instance by altering the dimension of dilatancy boundary layers.

To understand this, consider a cylindrical, very long cavern surrounded by an infinite salt-brine mixture; let us model it as a first-gradient elastic material and assume that the cylinder is subject to a steady internal pressure p 01 , while the outside boundary, ideally at infinity, is carrying the isotropic pressure p 02 ͑Fig. 1͒. Because both the cavern pressure and the outer boundary pressure, integrated over their individual closed boundaries, are in static equilibrium, p 01 and p 02 can individually be chosen; they thus define a multitude of isotropic stresses states. In particular, fluid can be pressed into or sucked from the cylinder, and the effect of such actions may be studied. Therefore, replacing p 01 by p 01 + ⌬p 01 ͑and applying the zero displacement condition at infin-ity͒ will, for a given p 02 , change the deformation state in the solid-fluid mixture around the fluid-filled cavern, in particular inducing some dilatancy boundary layers inside the solid matrix. This deformation state will depend upon the prestresses p 01 and p 02 . Of concern is how the apparent solid density ⌬ s changes due to the imposed additional ⌬p 01 as a function of the radial distance from the cavern wall.

For the first-gradient solid-fluid mixture model used in the present paper, the results are displayed in Fig. 2. If p 01 Ͻ p 02 , a linear increase in the cavern pressure results in a dilatancy, ⌬ s Ͻ 0, for all r ͓a , ϱ͒, where aϭradius of the cylinder; alternatively, when p 01 Ͼ p 02 , an increase in the cavern pressure will result in a compaction, ⌬ s Ͼ 0, for all r ͓a , ϱ͒. These results imply almost by default that, for p 01 = p 02 , ⌬ s = 0 for all r ͓a , ϱ͒. This is indeed the result that motivated us to improve the first-gradient model by a second-gradient model in our earlier paper ͑Sciarra et al. 2001͒ to recover dilatancy effects also in the wall-near boundary layer when p 01 = p 02 .

The monotonic behavior exhibited by the graphs of Fig. 2 prevails for all prestress conditions p 01 and p 02 . However, as p 01 increases beyond p 02 , the maximum value of ⌬ s at r = a increases and the boundary layer becomes thicker ͑Fig. 3͒.

1 Doctor, Dipartimento di Ingegneria Chimica, dei Materiali, delle Materie Prime e Metallurgia, Univ. of Rome "La Sapienza," via Eudossiana 18 00184 Rome, Italy. E-mail: giulio.sciarra@uniroma1.it 2 Professor, Dipartimento di Ingegneria Strutturale e Geotecnica, Univ. of Rome "La Sapienza," via Eudossiana 18 00184 Rome, Italy. E-mail: francesco.dellisola@uniroma1.it 3 Professor, Dept. of Mechanics, Darmstadt Univ. of Technology, D-64289 Darmstadt, Germany. E-mail: hutter@mechanik.tu-darmstadt.de

Note. Discussion open until February 1, 2006. Separate discussions must be submitted for individual papers. To extend the closing date by one month, a written request must be filed with the ASCE Managing Editor. The manuscript for this paper was submitted for review and possible publication on March 15, 2004;approved on December 9, 2004

Formulation of Problem

Consider a solid-fluid mixture; identify the solid and fluid constituents by the labels s and f, respectively. Material particles of the fluid and the solid are identified, respectively, by their position vectors X f and X s in fixed reference configurations ⍀ f 0 and ⍀ s 0 . As is usually done in the theory of mixtures, we presume that, at any time t, particles of both constituents occupy the same position x in the present configuration ⍀. The velocity v a ͑a = f , s͒ of the material particle X a is defined by

v a = d a u a ͑X a ,t͒ dt ͑1͒
where d a / dtϭmaterial time derivative following the motion of X a ; and u a ϭdisplacement of the ath constituent from its reference configuration. Let f and s denote, respectively, the apparent mass densities of the fluid and the solid; then the mass density, , of the mixture equals f + s . The mean barycentric velocity v of the mixture is defined by

v = f v f + s v s ͑2͒
We assume the reader to be familiar with the concepts of the theory of mixtures. Details can be taken from the literature, e.g., Truesdell ͑1957͒; Müller ͑1985͒; and Hutter and Jöhnk ͑2003͒.

Balance Laws

We assume a classical mixture of first-gradient formulation and use the principle of virtual power to derive the balance of linear momentum and the boundary conditions for each constituent; that is

͵ ⍀ ͑m s • v s + m f • v f + T s • ٌv s -p f div v f ͒dV = ͵ ⍀ ͑b s • v s + b f • v f ͒dV + ͵ ץ⍀ ͑t s • v s + t f • v f ͒dA ͑3͒
Here, m a ϭbulk solid-fluid interaction force; T s ϭpartial Cauchy stress in the solid; p f ϭhydrostatic pressure in the fluid ͑we assume that the partial Cauchy stress in the fluid is spherical͒; ٌϭgradient operator with respect to coordinates in the present configuration; b a ϭdensity of partial body forces; and t a ϭpartial surface tractions. The objectivity of the left-hand side of Eq. ͑3͒ implies that the sum of the two internal supplies m s and m f of linear momentum equals 0 and T s is symmetric.

By using the divergence theorem and exploiting the fact that Eq. ͑3͒ must hold for all virtual velocities, we obtain

div T s -m s + b s = 0, in ⍀ ͑4͒ -ٌp f -m f + b f = 0, in ⍀ ͑5͒ m s + m f = 0, in ⍀ ͑6͒ T s n = t s , on ץ ⍀ ͑7͒ -p f n = t f , on ץ ⍀ ͑8͒
Eqs. ͑4͒ and ͑5͒ are the momentum or force balance equations for the solid and fluid constituents, Eq. ͑6͒ expresses Newton's third law, and, Eqs. ͑7͒ and ͑8͒ are the boundary conditions of traction of the two constituents. The principle of virtual power can also allow us to take into account boundary conditions on the displace- ment field simply by restricting the space of virtual velocities to those fields that vanish on the boundary itself.

Constitutive Relations

The balance laws, Eqs. ͑4͒ and ͑5͒, are to be supplemented by constitutive relations; we express these in terms of the internal energy. We presume that the mixture is at a uniform temperature, the constituents are deforming quasi-statically so that their kinetic energy can be neglected, and no energy is dissipated. The internal energy density is assumed to exhibit only elastic properties of the solid and fluid constituents ⑀ = ⑀͑F s , f , X s ͒. Thus, the balance of energy takes the form

d v dt ͵ ⍀ ⑀͑F s , f ,X s ͒dV = ͵ ⍀ ͑b s • v s + b f • v f ͒dV + ͵ ץ⍀ ͑t s • v s + t f • v f ͒dA ͑9͒
in which F s ϭdeformation gradient of the solid; and d v / dtϭmaterial time derivative following the barycentric motion of the mixture. According to the Reynolds transport theorem, the following constitutive equations must hold:

T s = ⑀ץ ץF s F s T ͑10͒ p f = f ⑀ץ ץ f ͑11͒ m s = -ͫ f ٌ͑F s ͒ T ⑀ץ ץF s + f F s -T ⑀ץ ץX s -s ⑀ץ ץ f ٌ fͬ ͑12͒ in which f = f , s = s , = f + s ͑13͒
are the fluid and solid mass fractions and, by definition, s + f =1.

In this paper we limit attention to external actions for which b s = b f = 0, i.e., only external surface tractions are applied. When applying the interior pressure of the interface between the single fluid in the cavern, a suitable procedure must be known for how to divide this traction between the two constituents. Our earlier studies ͑Sciarra et al. 2001͒ have shown that, to avoid violation of the second law of thermodynamics, the existence of a potential function must be assumed such that the working W ˙ext of the external surface tractions is given by Computations were done with the values of Table 1 andp 02 = 18 MPa

. W ˙ext = d v dt ͵ ⍀ ext ͑x, s , f ͒dV ͑14͒
where W ˙ext ϭsurface integral on the right-hand side of Eq. ͑9͒.

The external surface tractions for which Eq. ͑14͒ holds are conservative.

Requiring that Eq. ͑14͒ holds for all choices of the velocity field, we obtain

ץ ext ץ s = C s , ץ ext ץ f = C f , in ⍀ ͑15͒ t s = - ץ ext ץ s s + s ext , on ץ ⍀ ͑16͒ t f = - ץ ext ץ f f + f ext , on ץ ⍀ ͑17͒
Eq. ͑15͒ with C s and C f as constants are necessary conditions for the existence of a ext for which b s = b f = 0.

The preceding relations must be complemented by explicit expressions for the internal energy ⑀ and the external potential ext , and both must account for the prestress acting in the reference configuration. Denoting these stresses by T s 0 and p f 0 , respectively, they must satisfy the equilibrium conditions div T s 0 = 0, ٌp f 0 = 0 ͑18͒

We assume an internal energy ⑀ that depends on T s 0 , H s = ٌu s , where u s ϭsolid displacement field; ␥ f and ⌬ f = ff 0 , quadratic in H s and ⌬ f . This argument yields

⑀ = 1 0 ͫ T s 0 • H s + ␥ f ⌬ f + 1 2 C͓H s ͔ • H s - 1 4 ͑T s 0 H s T -H s T s 0 ͒ • H s + 1 8 ͑H s T s 0 -T s 0 H s + T s 0 H s T -H s T T s 0 ͒ • H s + 1 2 ␥ f f ⌬ f 2 + K sf ⌬ f I • H sͬ ͑19͒
in which 0 = ͑ s + f ͒ 0 ϭmixture density and Cϭelasticity tensor:

C͓H s ͔ = 2 symH s + ͑tr H s ͒I ͑20͒
where and ϭLamé constants. Moreover, ␥ f , ␥ f f , and K sf ϭmaterial parameters. The spherical tensor K sf I accounts for the interaction between the solid and the fluid phases because of the deformation of the pores, while ␥ f and ␥ f f mimic fluid apparent compressibilities due to the initial pressure p f 0 and the additional load ⌬p f . A constant term has been omitted in Eq. ͑19͒. The first two members of the right-hand side of Eq. ͑19͒ are linear in H s and ⌬ f and thus account for the primary effect of the nontrivial stress state in the reference configuration. All other terms are second order and yield a linear contribution of H s and ⌬ f to the solid and the fluid stress. They also contain prestress contributions that may contribute to an anisotropic response behavior also if C͓•͔ is isotropic, as suggested by Eq. ͑20͒.

The presence of a prestress exerted in the solid as well as in the fluid phase accounts for the initial state of stress due to the lithostatic confining traction and the internal pressure of the cavern. This contribution, which is typically omitted, will allow the determination of a prediction for the disturbance induced by a change of the cavern pressure in terms of changes in the solid and fluid densities close to the cavern walls with respect to their initial values. We underline that, in the present paper, we do not explicitly discuss the range of values of T s 0 , ␥ f , and K sf which guarantees the second-order polynomial in Eq. ͑19͒ to be positive definite; in other words, no stability analysis has been developed. Apparently, further steps should be devoted to investigate the buckling limit of the mixture.

It has to be noticed that Eq. ͑19͒ may be obtained by means of a Taylor expansion of the free-energy truncated at the second order. Its physical interpretation must be based on the remark that, in the reference configuration, the mixture elementary volume element is not stress-free, and moreover its state of stress is not spherical. In particular, as it is very well known in the literature, the deformation energy depends on the whole deformation gradient and not only on its symmetric part.

The constitutive choice for ext will depend upon the form of the traction applied at the cylindrical boundary and at infinity. Thus, the subsequent analysis will be specialized for the assumed state of stress that we are analyzing in this paper. It will be supposed that the deformations of the mixture are axisymmetric; i.e., in cylindrical coordinates the two in-plane components of the displacement u r and u are functions of the radial coordinate r only. This implies that Eq. ͑18͒ possesses the solutions

T s 0 = ͩa 0 -b 0 r 2 ͪe r e r + ͩa 0 + b 0 r 2 ͪe e ͑21͒ p f 0 = ␥ f f 0 = const ͑22͒
in which e r and e ϭradial and azimuthal unit vectors. Therefore, the solid stress T s 0 only has a radial and an azimuthal component. Its isotropic part a 0 defines a homogeneous pressure, while its deviatoric part exhibits the normal stress effects. These have the largest contribution at the cylindrical wall and decay quadratically as r -2 as one moves to infinity. The radial stress is decreased from the ambient pressure a 0 by b 0 r -2 ; the hoop stress, on the other hand, is increased by this same amount. The fluid pressure Eq. ͑22͒ is uniformly distributed in the entire exterior of the cylinder and linearly related to the initial apparent fluid density f 0 . The conceivable most simple representation of ext is of the form

ext = 0 ext + ⌬ ext = ͩC s s 0 + C f f 0 + ¯0 + ¯1 r ͪ+͑C s ⌬ s + C f ⌬ f + ⌬ 0 ͒
͑23͒ in which C s ,C f ϭphenomenological constants; and ¯0, ¯1, and ⌬ 0 are related to the pressure at the boundaries in the reference and in the perturbed configuration. To find these connections, Eq. ͑23͒ is substituted into Eqs. ͑16͒ and ͑17͒. A routine, but somewhat lengthy, computation yields the solid-fluid traction vectors

t s 0 = ͫ-f 0 s 0 0 ͑C s -C f ͒ + s 0 ͩ ¯0 + ¯1 r ͪͬn, at r = a,r → ϱ t f 0 = ͫ f 0 s 0 0 ͑C s -C f ͒ + f 0 ͩ ¯0 + ¯1 r ͪͬn, at r = a ͑24͒
Because t 0 = t s 0 + t f 0 , we deduce from Eq. ͑24͒ that

t 0 ͑a͒ = ͩ ¯0 + ¯1 a ͪn=-p 01 n, at r = a lim r→ϱ t 0 ͑r͒ = ͑-f 0 s 0 0 ͑C s -C f ͒ + s 0 ¯0͒n = -p 02 n, as r → ϱ ͑25͒
Finally, there remains the determination of a 0 and b 0 in Eq. ͑21͒.

With the aid of Eq. ͑24͒ we may write

ͩa 0 -b 0 a 2 ͪn=ͫ-f 0 s 0 0 ͑C s -C f ͒ + s 0 ͩ ¯0 + ¯1 a ͪͬn, at r = a a 0 n = ͓-f 0 s 0 0 ͑C s -C f ͒ + s 0 ¯0͔n, as r → ϱ ͑26͒ from which a 0 = ͓-f 0 s 0 0 ͑C s -C f ͒ + s 0 ¯0͔, b 0 = -s 0 a ¯1 ͑27͒
may be deduced, once the following expressions for 0 and 1 are taken into account from Eq. ͑25͒:

¯0 = - p 02 s 0 + f 0 0 ͑C s -C f ͒, ¯1 = a p 02 s 0 -ap 01 -a f 0 0 ͑C s -C f ͒ ͑28͒
For C s = C f , which means that the external potential simply depends on the mixture density, Eqs. ͑24͒ and ͑25͒ simplify considerably.

Eq. ͑25͒ defines the boundary conditions for the zeroth order problem setting. In other words, for a given internal pressure p 01 of the cylindrical boundary and a prescribed outer pressure p 02 , the density distributions f 0 and s 0 can be computed in terms of p 01 and p 02 , once suitable constitutive relations have been taken into account for solid and fluid reference stresses. In the subsequent sections, this solution is regarded as a basic reference state from which deviations will now be determined.

Perturbed Linear Boundary Value Problem

Consider now the cylindrical cavity embedded in the solid-fluid medium subjected to the boundary pressures p 01 and p 02 as described previously. Assume that, by the application of some external action, the pressure p 01 on the boundary of the cavity is quasi-statically changed to p 01 + ⌬p 01 , while on the external boundary the pressure p 02 is kept. This can for instance be achieved by pressing liquid into the cavern or by removing it. This will change the state of stress in the surrounding solid-fluid medium. For ⌬p 01 Ͼ 0, the fluid will be pressed into the pore space, the pores will tend to open, and the solid component will experience a density change and a radial displacement field u r s that will attenuate to zero as r tends to infinity.

It will now be assumed that all fields depend only on the radial coordinate and that the solid displacement field ͑u r s , u s ͒ as well as the change in densities ⌬ f and ⌬ s and pressure ⌬ 0 are small, such that nonlinear terms can be dropped. The linearized field equations, derivable from Eqs. ͑4͒, ͑5͒, ͑10͒, ͑11͒, and ͑12͒, then take the forms

ͫ 2 f 0 ͩa 0 -b 0 r 2 ͪ++2ͬu r Љ + 1 r ͫ 2 f 0 ͩa 0 + b 0 r 2 ͪ+ 2b 0 r 2 + + 2 ͬ u r Ј + 1 r 2 ͫ -2 f 0 ͩa 0 + b 0 r 2 ͪ-2b 0 r 2 -͑ + 2͒ ͬ u r + ͫ 1 0 ͩa 0 -b 0 r 2 ͪ+K sf -s 0 ␥ fͬ ⌬ f Ј = 0 ͑29͒ 1 r 2 ͫͩ a 0 2 + ͪr 2 -b 0ͬ u Љ + 1 r 3 ͫͩ a 0 2 + ͪr 2 -b 0ͬ u Ј - 1 r 4 ͫͩ a 0 2 + ͪr 2 -b 0ͬ u = 0 ͑30͒ ⌬ f Ј = f 0 ͓ s 0 ␥ f -0 K sf ͔ͩu r Ј + 1 r u r ͪ Ј -f 0 ͫͩ a 0 -b 0 r 2 ͪu r Ј ͬ Ј -ͫ f 0 ͩa 0 + b 0 r 2 ͪ 1 r u rͬ Ј 2 f 0 ␥ f + 0 f 0 ␥ f f ͑31͒
valid in a Ͻ r Ͻϱ, and the first-order boundary conditions follow from the ansatz Eq. ͑23͒ and the localizations Eqs. ͑16͒ and ͑17͒:

⌬t s = ͭͫ -͑ f 0 ͒ 2 ͑C s -C f ͒ + f 0 0 ͩ ¯0 + ¯1 r ͪͬ⌬ s + ͫ-͑ s 0 ͒ 2 ͑C s -C f ͒ + s 0 0 ͩ ¯0 + ¯1 r ͪͬ⌬ f + s 0 ⌬ 0ͮ n, at r = a ͑32͒ ⌬t f = ͭͫ ͑ s 0 ͒ 2 ͑C s -C f ͒ + s 0 0 ͩ ¯0 + ¯1 r ͪͬ⌬ f + ͫ͑ f 0 ͒ 2 ͑C s -C f ͒ - f 0 0 ͩ ¯0 + ¯1 r ͪͬ⌬ s + f 0 ⌬ 0ͮ n, at r = a ͑33͒
where n =-e r at r = a. Adding Eqs. ͑32͒ and ͑33͒ yields ⌬p 01 = -⌬ 0 , at r = a ͑34͒ thus identifying the perturbation ⌬ 0 with the perturbation pressure ⌬p 01 applied at the cylindrical wall and defining the solid and fluid boundary pressures according to A detailed derivation of Eqs. ͑29͒-͑33͒ is given in the Appendix. Thus, the proposed steady-state deformation problem appears mathematically as a standard linear two-point boundary value problem comprising the ordinary differential Eqs. ͑29͒-͑31͒ and the boundary conditions Eqs. ͑32͒, ͑33͒, ͑37͒, and ͑38͒. The differential equations are homogeneous and so are the boundary conditions Eqs. ͑37͒ and ͑38͒, but the traction boundary conditions are inhomogeneous via the prescription of ⌬p 01 .

⌬p s ͑a͒ = -ͭͫ -͑ f 0 ͒ 2 ͑C s -C f ͒ + f 0 0 ͩ ¯0 + ¯1 r ͪͬ⌬ s + ͫ-͑ s 0 ͒ 2 ͑C s -C f ͒ + s 0 0 ͩ ¯0 + ¯1 r ͪͬ⌬ f + s 0 ⌬ 0ͮ ͑35͒ ⌬p f ͑a͒ = -ͭͫ ͑ s 0 ͒ 2 ͑C s -C f ͒ + s 0 0 ͩ ¯0 + ¯1 r ͪͬ⌬ f + ͫ͑ f 0 ͒ 2 ͑C s -C f ͒ - f 0 0 ͩ ¯0 + ¯1 r ͪͬ⌬ s + f 0 ⌬ 0ͮ ͑36͒ In
Inspection of Eqs. ͑30͒ and ͑37͒, and the second part of Eq. ͑38͒ shows that the solution for u is decoupled from those of Eqs. ͑29͒, ͑31͒-͑33͒, and the first part of Eq. ͑38͒. Now, because the boundary value problem for u is homogeneous, u = 0 is a solution; in fact, this trivial result follows simply by the symmetry of the boundary value problem. What remains are the differential Eqs. ͑29͒ and ͑31͒ subjected to the boundary conditions Eqs. ͑32͒ and ͑33͒, and the first part of Eq. ͑38͒. The construction of their solution will be our next task.

Solution of Boundary Value Problem

The coupled system Eqs. ͑29͒ and ͑31͒ of ordinary differential equations can be reduced to a simple second-order differential equation for u r if Eq. ͑31͒ is substituted into Eq. ͑29͒. This new differential equations has the form a͑r͒u r Љ + b͑r͒u r Ј + c͑r͒u r = 0 ͑39͒ with coefficients that are singular at r = 0 and are such that Eq. ͑39͒ is not of the Fuchs type. It follows that no convergent series solutions can be found. Thus, it is hopeless to search for analytical solutions. Solutions must be constructed by numerical techniques, and these are not likely to offer any difficulties, because the only singularities of a͑r͒, b͑r͒, and c͑r͒ lie outside the range of integration, ͑a , ϱ͒, namely, at r =0.

For the numerical integration of the differential Eqs. ͑29͒ and ͑31͒, it is advantageous to reduce it to first-order standard form and simultaneously to nondimensionalize the equations. To this end, we select

r = ax, u r = ay 1 , du r dr = y 2 , ⌬ f = f 0 y 3 ͑40͒
With this choice, it is straightforward to show that Eqs. ͑29͒ and ͑31͒ reduce to

A͑x͒YЈ = B͑x͒Y ͑41͒
͑29͒ where Y T = ͑y 1 , y 2 , y 3 ͒; and

A͑x͒ = 1 0 0 0 A͑x͒ 0 B͑x͒ = 0 1 0 B͑x͒ 0 0 ͑42͒
in which A͑x͒ and B͑x͒ can easily be given in terms of the coefficients of Eqs. ͑29͒ and ͑31͒. Inverting A͑x͒, Eq. ͑41͒ may be written as

YЈ = A͑x͒ -1 B͑x͒Y = C͑x͒Y = 0 1 0 C͑x͒ 0 0 ͑43͒
This is now the standard form of a vector-valued first-order ordinary differential equation, which must be integrated between x = 1 and x → ϱ. The boundary conditions to which it is subjected can be deduced from Eqs. ͑32͒ and ͑33͒, and the first part of Eq. ͑38͒, and they imply that

a 1 y 1 ͑1͒ + a 2 y 2 ͑1͒ = P, y 3 ͑1͒ = y ¯3, at r = a lim x→ϱ y 1 ͑x͒ = 0, at r → ϱ ͑44͒
where P and y ¯3, as well as a 1 and a 2 , are prescribed constants that can be expressed in terms of p 01 , p 02 , ⌬p 01 , and all the constitutive parameters we have introduced.

The approach to solve Eqs. ͑43͒ and ͑44͒ was a fourth-order Runge-Kutta method combined with the shooting method. To this end the initial value problem

YЈ = C͑x͒Y Y͑1͒ = Y ¯͑45͒
was solved with arbitrarily selected y 1 ͑1͒, and this value was systematically altered until lim x→ϱ y 1 ͑x͒ = 0 was obtained with a preselected error ͉͉ Ӷ 1 ln. The actual computation x → ϱ was replaced by a large finite number x ϱ = 1 , 000 and ͉y 1 ͑x ϱ ͉͒ Ͻ =10 -15 was chosen. Furthermore, the shooting method in a linear problem that does not exhibit stiff properties only requires three integrations to achieve the final solution. The integrations were conducted with the use of the software MATHEMATICA.

Results

Numerical solutions of the two-point boundary value problem Eqs. ͑43͒ and ͑44͒ were sought for fixed assigned values of the outside reference pressure p 02 and several values of p 01 starting at p 01 Ͻ p 02 and ending at p 01 Ͼ p 02 ; a numerical value for p 02 is the overburden pressure of a salt cavern underground, and p 01 can vary from atmospheric pressure to approximately 2p 02 . Numerical values for all the phenomenological coefficients are collected in Table 1. Inspection of all the formulas in the previous section reveals that the constants C s and C f of the external potential only arise as ͑C s -C f ͒, so results are insensitive to the numerical values of these, if we choose C s = C f . The ensuing discussion will be restricted to this case. 2͑b͒ shows the results for p 01 = p 02 , and it states that ⌬ s = 0 for all x ͓1,ϱ͒, which is a trivial result. The two load cases p 01 Ͻ p 02 and p 01 Ͼ p 02 demonstrate qualitatively different physical behavior. Indeed, for the first case, Fig. 2͑a͒, the numerical values show ⌬ s Ͻ 0. This means that the pore space becomes bigger, ⌬ f Ͼ 0, where f is the porosity, corresponding to a dilatation. For the second case, ⌬ s Ͼ 0 for all x ͓1,ϱ͒, implying that the pore space becomes smaller, ⌬ f Ͻ 0, corresponding to a compaction. As a function of x, this behavior is monotonous; i.e., the dilatation and compaction, respectively, decrease with increasing x and seem to approach zero as x → ϱ.

When varying p 01 parametrically from p 01 = 1.8 MPaϽ p 02 in steps ␦p 01 = 5.4 MPa until p 01 = 34.2 MPaϾ p 02 , the graphs analogous to those of Fig. 2 are as shown in Fig. 3. Evidently, the maximum value of the dilatation/compaction occurs at the cylinder wall, and it increases until increasing the difference ͉p 01p 02 ͉. Moreover, as ͉p 01p 02 ͉ increases, the boundary layer with appreciable values of the dilatation/compaction thickens. Thus, the larger the difference ͉p 01p 02 ͉ are, the less localized the reactions to this difference will be. This behavior is evidently consistent with the identically vanishing solution of the apparent density, which is reached at p 01 = p 02 . An estimate of the boundary layer thickness ͑the thickness ␦ of the boundary layer is the distance from the cylinder wall at which the derivative of ⌬ s becomes vanishing͒ has been determined computationally and is summarized in Fig. 4.

The preceding results may equally be explained somewhat differently. To this end, consider a representative elementary volume ͑REV͒ comprising idealistically of four solid quadrants and a fluid sphere in the middle with connections to the REV sides ͑Fig. 5͒. When p 01 Ͻ p 02 ͓Fig. 5͑a͔͒, the dilatation is associated with an increase of the fluid volume that is filled by fluid injection from the neighboring elements or the cavern fluid. However, the true solid density and therefore the solid volume in the REV remain the same. The pressure of the solid-essentially the negative trace of the solid stress-is positive, a compression. Thus, the change of the apparent density of the solid is negative ⌬ s Ͻ 0. The effect of the fluid via its seepage flow is to enlarge the porosity, ⌬ f Ͼ 0, because fluid is injected.

When p 01 Ͼ p 02 ͓ Fig. 5͑b͔͒, the compaction is achieved by a reduction of the pore space, ⌬ f Ͻ 0, in the REV, and fluid is squeezed out. The solid stress now is under tension, but the true solid density does not change so that ⌬ s Ͼ 0.

In both cases this behavior is the result of seepage flow and effects of the fluid rather than deformation effects in the solid component within the REV.

Accompanied with the changes in ⌬ s are also changes in the fluid pressure, ⌬p f . Fig. 6 shows graphs for ⌬p f / p 01 against x for the same conditions as Fig. 2. For p 01 Ͻ p 02 ͓Fig. 6͑a͔͒, the perturbation in the fluid pressure is largest at the cylinder wall, and it decreases monotonically to a nonzero value as x → ϱ. Alternatively, for p 01 Ͼ p 02 ͓Fig. 6͑b͔͒, the perturbation pressure takes its minimum value at the cylinder wall and approaches a large asymptotic value as x → ϱ. In both cases typical boundary layer effects are seen.

This behavior, too, can be envisaged by sketches of typical REVs. We start from a reference configuration of the solid-fluid mixture, of which the pore structure is characterized by a positive or a negative porosity gradient when the confining state of stress satisfies the condition p 01 Ͻ p 02 ͓Fig. 7͑b͔͒ or p 01 Ͼ p 02 ͓Fig. 7͑d͔͒, respectively. This kind of equilibrium configuration for the fluid-saturated rock salt can be achieved, for instance, as the ultimate fate of a transient creep process tending to stabilize after a suitable length of time ͑Cristescu and Hunsche 1998͒ ͓see, e.g., the transition from Figs. 7͑a and c͒ to Figs. 7͑b andd͔͒; local creep deformations enhance the convergence of the solid within the REV and consequently the closure of the pores. The smaller the confining pressure, the more relevant this effect is; for this reason, we can appreciate a positive or negative gradient of porosity in the reference configuration when p 01 Ͻ p 02 or p 01 Ͼ p 02 , respectively. This hypothesis, endowed with the assumption of uniform apparent densities of the solid and fluid constituent in the initial equilibrium configuration, evidently implies the true densities ˆa 0 = a 0 / a 0 ͑a = s , f͒ not to be uniform with r. The deformation of the considered reference configuration is indeed consistent with the aforementioned results regarding the fluid pressure ⌬p f due to a positive ⌬p 01 . As we are dealing with a linear elastic model, the positive pressure ⌬p 01 , applied on the boundary of the cavity, is likely to imply an injection of fluid in the connected pores of the salt matrix.

Bearing in mind the previous results on the solid density profiles, let us try to discuss what happens to porosity when incompressible small deformations of the solid constituent ͑⌬ ˆs =0͒ take place; apparently, the following conditions have to be fulfilled:

⌬ s = ⌬͑ ˆs s ͒ = ⌬ ˆs s 0 + ˆs 0 ⌬ s = -ˆs 0 ⌬ f ͑46͒ ٌ⌬ s = ٌ⌬͑ ˆs s ͒ = -ٌ͑ ˆs 0 ⌬ f ͒ = -ٌ͑ ˆs 0 ͒⌬ f -ˆs 0 ٌ⌬ f ͑47͒
Because the porosity in each point X of the mixture corresponds to the ratio between the fluid volume in the REV ͑V f ͒ at X and the volume of the REV ͑V͒ itself, the following allusive chain of equalities can be established:

⌬ f = ⌬ͩ V f V ͪ= ⌬V f V 0 - V f 0 V 0 2 ⌬V = ⌬V f V 0 -f 0 ⌬V V 0 ͑48͒ ٌ͑⌬ f ͒ = ٌͩ ⌬V f V 0 ͪ-ٌͩ f 0 ⌬V V 0 ͪ = ٌͩ ⌬V f V 0 ͪ-f 0 ٌͩ ⌬V V 0 ͪ-ٌ f 0 ⌬V V 0 ͑49͒
where ⌬ f = ff 0 , i.e., the difference between the current and the

= ⌬ ˆs = - ˆs 0 V s 0 ⌬V s ⇒ ⌬V = ⌬V f ͑50͒
which implies, because of Eqs. ͑48͒ and ͑49͒

⌬ f = ͑1 -f 0 ͒ ⌬V V 0 ͑51͒ ٌ͑⌬ f ͒ = ͑1 -f 0 ٌ͒ͩ ⌬V V 0 ͪ-ٌ f 0 ⌬V V 0 ͑52͒
For p 01 Ͻ p 02 , the reference porosity gradient is positive ٌ͑ f 0 Ͼ 0͒; moreover, the dilatation of the REV ͑⌬ s Ͻ 0, say, ⌬ f Ͼ 0͒ progressively decreases ٌ͑⌬ s Ͼ 0͒ when r → ϱ. Incompressible small deformations of the solid constituent imply, according to Eq. ͑47͒, the variation of porosity to decrease with r, say, ٌ͑⌬ f ͒ Ͻ 0, as ٌ ˆs 0 = s 0 / ͑ s 0 ͒ 2 ٌ f 0 Ͼ 0 ͓see Fig. 7͑a and c͔͒. When the internal confining pressure is smaller than the external one, dilatancy occurs, and this effect progressively reduces going away from the cavern wall. As the behavior of the apparent density of the solid is essentially triggered by the change in volume of the REVs, the inequalities ⌬V / V 0 Ͼ 0 and ٌ͑⌬V / V 0 ͒ Ͻ 0 must hold true, and therefore ⌬V f / V 0 Ͼ 0 and ٌ͑⌬V f / V 0 ͒ Ͻ 0, because of Eq. ͑50͒. Thus, the macroscopic behavior described in Fig. 2͑a͒ is consistent with the microinterpretation given by Eq. ͑52͒.

On the other hand, when p 01 Ͼ p 02 , the reference porosity gradient is negative ٌ͑ f 0 Ͻ 0͒; moreover, compaction of the REV ͑⌬ s Ͼ 0, say, ⌬ f Ͻ 0͒ progressively decreases ٌ͑⌬ s Ͻ 0͒ as r tends to infinity. Incompressible small deformations of the solid constituent ͑⌬ ˆs =0͒ require, according to Eq. ͑47͒, the following condition to hold true:

ٌ⌬ f ͉⌬ f ͉ Ͼ ٌ f 0 1 -f 0 , ٌ f 0 /͑1 -f 0 ͒ Ͻ 0 ͑53͒
which does not univocally determine the sign of ٌ⌬ f . As the behavior of the apparent density of the solid is essentially triggered by the change in volume of the REVs, the inequalities

⌬V / V 0 Ͻ 0 ͑⌬V f / V 0 Ͻ 0͒ and ٌ͑⌬V / V 0 ͒ Ͼ 0ٌ͑͑⌬V f / V 0 ͒ Ͼ 0͒ must
hold true; these conditions are really consistent with the previous remarks on the sign of the gradient of porosity, according to Eq. ͑52͒.

As the mean pressure of the fluid p f can be essentially characterized, via a "naive" micro-macro identification, by the product p f , where p indicates the porewater pressure, the variation of the fluid pressure ⌬p f will be

⌬p f = p 0 ⌬ f + ⌬p f 0 ͑54͒
The pressure of the fluid at the micro level is essentially constant ٌ͑p = ٌ͑p 0 + ⌬p͒ =0͒. When p 01 Ͻ p 02 , the injection of a small amount of fluid mass into the cavern implies ٌ⌬p f Ͻ 0 ͓see Fig. 6͑a͔͒, and consequently

ٌ⌬ f Ͻ - ⌬p p 0 ٌ f 0 Ͻ 0 ͑55͒
which is consistent with Eq. ͑52͒; on the other hand, when p 01 Ͼ p 02 , the injection of the same amount of fluid into the cavern implies ٌ⌬p f Ͼ 0 ͓see Fig. 6͑b͔͒, and consequently

ٌ⌬ f Ͼ - ⌬p p 0 ٌ f 0 Ͼ 0 ͑56͒
which is consistent with Eq. ͑53͒.

Conclusion

In this paper, an analysis of the deformational behavior of the solid-fluid surrounding ground of a fluid filled cavern was pre- The results show that the solid density change, ⌬ s , associated with ⌬p 01 = 1, as a function of the radial distance from the cavern is a dilatation ͑compaction͒ if p 01 Ͻ p 02 ͑p 01 Ͼ p 02 ͒ ͑Fig. 2͒, with the maximum modulus arising at the cavern wall and quickly tapering off, thus giving rise to a boundary layer of pore space opening ͑or closure͒ close to the cavern wall. This nontrivial dilatancy is solely due to the state of prestress, and it collapses when p 01 = p 02 and when ⌬ s = 0 for r Ͼ a.

We note that the parametric analysis which has been developed represents a first step toward the study of the behavior of salt caverns, in particular when the internal pressure is greater than the lithostatic one. We noticed that the dilatancy phenomenon, induced by the injection of a liquid mass into the cavern, can occur, in the elastic regime, if the internal pressure is smaller than the external one. On the other hand, when the initial pressure of the cavern is greater than the external one, compaction is detected in the elastic regime. Apparently, dilatancy phenomena can also occur when reaching these pressures inside the cavern once the yield limit of the salt had been surpassed, but this is a wellknown issue. We prove that, in order to avoid dilatancy, it is not sufficient to keep the internal pressure below the lithostatic limit but also the conditions to assure an elastic regime have to be investigated.

Our analysis has obviously been restrictive in several respects and ought to be extended. First, the state of prestress does not contain shear stresses and may thus overlook a significant state of stress arising in practical applications. For instance, at the far-end boundaries, the vertical and horizontal confining pressures would differ from one another and, thus, induce a combination of an isotropic pressure plus a pure shear. Its analysis will be rather complicated. Second, the solutions of the linearly perturbed problem have not been analyzed with regard to their stability. More specifically, depending upon the initial state of stress, the perturbed solution may be linearly stable or unstable. Delineating the ranges of these stabilities ͑instabilities͒ is practically significant.

After the generalizations and improvements just delineated, it will also become possible to cope with the problem of formulating and applying a pore-plasticity model to describe more precisely the behavior of the disturbed rock zone around cavities.

Acknowledgment

This paper was written while the third writer was on sabbatical leave at the University of Rome, La Sapienza, the support of which is hereby acknowledged.

Appendix.

The constitutive Eqs. ͑10͒-͑12͒, together with the internal energy constitutive relation Eq. ͑19͒, are specialized bearing in mind the axisymmetric assumption of the deformation of the mixture; as in cylindrical coordinates, if the two in-plane components of the displacement u r and u depend only on r, we have

T s = ⌬ f 0 T s 0 -s 0 ͑tr H s ͒T s 0 + 1 2 ͑T s 0 H s T + H s T s 0T ͒ + 1 2 ͑W s T s 0 -T s 0 W s ͒ + C͓H s ͔ + ⌬ f K sf I ͑57͒ p f = ͓␥ f ͑1 + f 0 ͒ + f 0 ␥ f f ͔⌬ f + f 0 ͑K sf -s 0 ␥ f ͒tr H s ͑58͒ m s = -f 0 ٌ͑T s 0 • H s ͒ + s 0 ␥ f ٌ͑⌬ f ͒ ͑ 59͒
These last yield the following formulas for the divergence of the solid stress tensor and the gradient of the fluid pressure:

div

T s = T s 0 ͫ 1 0 ٌ͑⌬ f ͒ -s 0 ٌ͑tr H s ͒ ͬ + 1 2 ٌ͑T s 0 ͒H s T + 1 2 T s 0 div H s T + 1 2 ٌ͑H s ͒T s 0 + 1 2 ٌ͑W s ͒T s 0 + 1 2 T s 0 div W s + div C͓H s ͔ + ٌ͑⌬ f ͒K sf ͑60͒ ٌp f = ͓␥ f ͑1 + f 0 ͒ + f 0 ␥ f f ٌ͔͑⌬ f ͒ + f 0 ͑K sf -s 0 ␥ f ٌ͒͑tr H s ͒ ͑61͒
Eqs. ͑60͒ and ͑61͒, endowed with the following representation of the displacement gradient H s , assumed to depend just on r:

H s = u r,r e r e r + u ,r e r e -1 r u e e r + 1 r u r e e ͑62͒

imply the field Eqs. ͑29͒-͑31͒, once the radial and the azimuthal components of Eq. ͑4͒ as well as Eq. ͑5͒ are taken into account. Note that we adopt the following rule for the tensor product ͑a b͒u = ͑a • u͒b. Boundary conditions Eqs. ͑32͒ and ͑33͒ are then achieved, starting from the external potential ext ͓Eq. ͑23͔͒ and its partial derivatives ͓Eq. ͑15͔͒ with respect to the solid and fluid apparent densities, by replacing them into Eqs. ͑16͒ and ͑17͒.

Fig. 1 .

 1 Fig. 1. Cylindrical cavern of radius a filled with an inviscid fluid, in a solid-fluid mixture of polycrystalline salt and brine. The porous salt has a very low permebility. Plane deformation is assumed.

Fig. 2 .

 2 Fig. 2. Variation of the solid normalized apparent density change ⌬ ¯s = ⌬ s / s 0 plotted against x to a change in cavern pressure ⌬p 01 = 1 MPa. Computations were done for phenomenological constants of Table 1: ͑a͒ p 02 = 18 MPa and p 01 = 7.2 MPa; ͑b͒ p 01 = p 02 ; ͑c͒ p 01 = 28.8 MPa..

Fig. 3 .

 3 Fig. 3. Variation of the normalized solid perturbation density ⌬ ¯s = ⌬ s / s 0 for ⌬p 01 = 1 MPa and different values of the cavern pressure p 01 , as indicated, and for ͑a͒ p 01 Ͻ p 02 and ͑b͒ p 01 Ͼ p 02 . For increasing p 01 , the value of ⌬ ¯s at the cavern wall increases and the boundary layer thickens. Computations were done with the values of Table1 and p 02 = 18 MPa.

  addition, the vanishing shear stresses in r = a require u Ј͑a͒ = 0 ͑37͒ and the boundary conditions at r → ϱ require lim r→ϱ u r ͑r͒ = lim r→ϱ u ͑r͒ = 0 ͑38͒

  Fig. 2 displays the dimensionless solid density change ⌬ s / s 0 plotted against the distance x = r / a for the particular values p 02 = 18 MPa and p 01 = 7.2 MPaϽ p 02 ͓Fig. 2͑a͔͒, p 01 = 28.8 MPaϾ p 02 ͓Fig. 2͑c͔͒. Fig.

Table 1 .Fig. 4 .

 14 Fig. 4. Scaled boundary layer thickness ␦ / a plotted against the scaled initial stress difference for p 02 = 18 MPa and p 01 parametrically varying from p 01 = 1.8 MPaϽ p 02 in steps ␦p 01 = 5.4 MPa until p 01 = 34.2 MPaϾ p 02 . For ͉p 01p 02 ͉ → 0, the thickness, of the boundary layer tends to zero ␦ / a → 0.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Representative elementary volume ͑REV͒ with solid and fluid ͑interconnected with neighboring cells͒ in the reference configuration ͑left͒ and in the present configuration ͑right͒: ͑a͒ shows behavior when p 01 Ͻ p 02 : the REV expands by fluid injection rather than solid deformation and thus yields an increase of the pore space; and ͑b͒ shows behavior when p 01 Ͼ p 02 : the REV contracts and fluid is squeezed out

Fig. 7 .

 7 Fig. 7. Rows of representative elementary volumes from the cavern wall through to infinity. The first row shows them before any load is applied with uniform solids and fluids fraction. When p 01 and p 02 are applied from the left and right, respectively, the difference in p 01 and p 02 induces an inhomogeneous distribution of the pore space. For ͑b͒ p 01 Ͻ p 02 , the porosity increases with distance from the cavern wall; for ͑d͒ p 01 Ͼ p 02 , it decreases. This corresponds to a fluid pressure change ⌬p f that is decreasing ͓when starting from ͑a͒ the reference configuration toward ͑b͒ the current configuration panel͔ and increasing ͓when starting from ͑c͒ the reference configuration toward ͑d͒ the current configuration panel͔, respectively, with x.

  sented, when this material is subjected to a confining isotropic stress p 02 from infinity and a different pressure p 01 from the circular cavern wall. Both constituents are treated as elastic, but the linear law takes the prestress that is established by the external pressures p 01 and p 02 into account. The model is set up only to generate radial dependencies of the process variables when the cavern pressure is changed from p 01 to p 01 + ⌬p 01 , where, owing to linearity, it suffices to analyze ⌬p 01 =1.


Mario Alves

Santos Junior

Carlos Antonio

França Sartori

José Roberto Cardoso

Damien Voyer

Laurent Krähenbühl

Ronan Perrussel

Christian Vollaire

Djonny Weinzierl

M A Santos

C A F Sartori

C Vollaire Ampère

A Load Effect Evaluation of a Transmission Line Exciting Chamber

Keywords: reverberation chamber, transmission lines, random excitations and loads I

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

anonical chambers -like Reverberation Chambers (RC) and TEM chambers -are generally used for electromagnetic immunity testing despite their particular operational restrictions. RCs using mechanical paddles or frequency stirring provide a statistical E-field uniformity in all the directions inside the working volume [START_REF] Hill | Electronic mode stirring for reverberation chambers[END_REF]. Nevertheless, the frequency operation of RCs is inversely proportional to the chamber dimensions and it is a constraint for low frequency tests. The International Standards recommend the RC configuration for immunity tests over 80 MHz frequencies [START_REF]Electromagnetic Compatibility (EMC), Part 4, Testing and Measurement Techniques -Section 21: Reverberation Chamber Test Methods[END_REF]. On the contrary, TEM chambers operate in the low frequency range. Introducing a stripline inside the chamber, a deterministic E-field uniformity is reached over a working area parallel to the plate, but not for all the polarizations in the chamber volume. The performances can however be extended for all the polarizations using three orthogonal striplines [START_REF] Klingler | On the use of 3D TEM Cells for Total Radiated Power Measurements[END_REF]. Recently, an alternative concept called Transmission Line Excitated Chamber (TLEC) has been proposed, based on a phase shifting excitation of several transmission lines (TL). For the sake of illustration, a structure constituted of three conductors with a phase shifting excitation has been investigated in [START_REF] Perini | Extending the frequency of mode stir chambers to low frequencies[END_REF]. Basically, a chamber excited by several TLs presents several TEM modes inside the closed metallic cavity. The resulting standing waves depend on the position of the TLs, on the amplitude and phase of the excitations and on the loading at the end of the TLs. Those parameters are important in the search for a suitable chamber working volume; they can be modified electronically, resulting in a random standing wave profile. Based on this, the performances of the TLEC can be improved to satisfy the pre-defined uniformity criteria within a wide frequency range, even at frequencies lower than 80 MHz. In this work, we present semi-analytical and numerical approaches for evaluating the performance of a TLEC.

II. SEMI-ANALYTIC APPROACH

A. Analytic expression for a TL

Consider the TL geometry given by Fig. 1. When one analyzes the E-field of the single TEM mode, the dependence with z longitudinal direction can be decoupled from the ones with the transverse x and y directions. Indeed, the propagation, which exists for any frequency, can be described by a function e -jk 0 z when an incident wave is considered; k 0 is the wave number in the vacuum. The problem of the impact of possible reflections at the end of the TL will be addressed further in section II.C. The E-field of a TEM mode in the transverse plane can be evaluated by the following analytical expressions:

( ) ( ) ( ) ∑ ∞ + =        <         + >         - ×         = 1 m 2 m 2 1 m 1 0 X TEM 0 y l y a m sh 0 y l y a m sh x a m cos a 2 y , x E π α π α π η (1) C ( ) ( ) ( ) ∑ ∞ + =        <         + >         - ×         = 1 m 2 m 2 1 m 1 0 Y TEM 0 y l y a m ch 0 y l y a m ch x a m sin a 2 y , x E π α π α π η (2) with 2 1 1 1 2 m m th l a J m m ch l a m m th l th l a a π π α π π             =     +         , 1 2 2 1 2 m m th l a J m m ch l a m m th l th l a a π π α π π             = -     +         (3)
where η 0 is the vacuum wave impedance and J m the harmonic coefficients related to the current density on the central conductor. Assuming that the conductor is an infinitely thin wire along z axis, one finds:

        = 0 0 m x a m sin a 2 I J π (4)

B. Several lines and Phase Shifting

When several TLs are ended with the same load, the E-field can still be written using the separation of variables in the transverse plane and in the longitudinal direction. The problem can then be treated by superposition:

( ) ( ) ∑ = i TEM i i TOT y , x E I y , x E (5)
where E i TEM is the E-field due to the i th TL for a current of amplitude 1 and I i represents the applied excitation current. The phase shifting between the TLs has then an effect on the repartition of the E-field in the transverse plane.

C. Load Shifting

The load at the end of the TLs introduces a reflection coefficient Γ that affects the longitudinal repartition of the field:

( ) ( ) ( ) z jk z jk TOT 0 0 e e y , x E z , y , x E + - + × = Γ (6)
Concerning the stripline used in canonical TEM chambers, the design of the TL is such that there is no reflection. Then the magnitude of the field is uniform in z direction. Using a TL with a thin wire, there is a discontinuity at the end of the line and it is difficult to match the line with the suitable load; then |Γ | ≠ 0 and a stationary wave is expected.

However, the maximum of E-field can be moved changing the phase of Γ. Consider the worst case |Γ | = 1. Suppose one can change the phase φ of Γ = e jφ ; this can be achieved for example with a phase shifter ended by a short circuit. Then the complete E-field is given by:

( ) ( ) ( ) ( ) 2 / j 0 TOT z jk j z jk TOT e 2 z k cos y , x E 2 e e e y , x E z , y , x E 0 0 φ φ φ       + × = + × = + - (7)
the phase shifting φ/2 in the cosinus function shows that the maxima and minima of the E-field in z direction move with φ.

When the phase varies linearly between 0 and 2π, the mean effect can be calculated as follows:

( ) ( ) ( ) y , x E 4 d 2 1 z , y , x E z , y , x E TOT 2 0 average π φ π π = = ∫ (8)
the average E-field is no more dependent on z direction. Thus, it is possible to homogenize the field in z direction even if the TL is unmatched.

III. NUMERICAL APPROACH

Numerical evaluation was performed with the Finite Integration Technique by using the commercial software CST-MWS ® [5]. The interest is that this approach is more realistic than the semi-analytic one since it takes into account the connection of the TLs outside the chamber; those details introduce discontinuities that can have an impact on the distribution of the E-field. A 3D TLEC model is built considering PEC walls and PEC TLs. The loads are set using ports defined at the end of the TLs while the phase shifting excitation is implemented by a post-processing approach.

A procedure applying Matlab® Activex commands was performed in order to set up the calculation parameters at CST-MWS® and to perform the indexes of merit calculation discussed in the following section.

IV. INDEXES OF MERIT

The uniformity of the field is evaluated using the indexes of merit defined in [START_REF] Weinzierl | Numerical evaluation of non-canonical reverberation chamber configurations[END_REF]. Due to the phase and load shifting, the Efield considered at any point of the chamber is the average field E average in time. Considering a spatial point k in the chamber, one finds for the x polarization of the average field:

( ) ( ) [ ] ∫ ∏ + = + = 1 M 2 , 0 M 1 i i 1 M M 1 k x k x _ average d d 2 , , , E E π ϕ φ π ϕ ϕ φ … (9)
with φ the phase of Γ, ϕ i the phase of the current I i in the i th TL and M the number of TLs. Numerically, the integral is computed using a finite number of phases.

Then the mean field x E and the standard deviation σ x in the spatial domain can be computed from the average field given in (9):

∑ = = N 1 k k x _ average x E N 1 E ( ) ∑ = - - = N 1 k 2 x k x _ average x E E 1 N 1 σ (10)
with N the number of points inside the chamber where the field is evaluated.

The normalized standard deviation in dB ˆx σ in x direction is defined by: 10 ˆ20log

x x x x E E σ σ   + =     (11)
The indexes in y and z directions are calculated similarly.

V. APPLICATION AND RESULTS

A TLEC with dimensions of 0.6m×0.6m×1.2m (see Fig. 2) has been considered; the volume under evaluation (WV) is a parallelepiped of dimensions 0.3m×0.3m×0.6m centered at the middle of the chamber. Besides that, the area under evaluation (WA) is a square of dimensions 0.3m x 0.3m at the center of WV. All the simulations are performed for a frequency of 80 MHz.

A. Semi-analytical results

The semi-analytical approach has been applied in the transverse plane {x, y} since the variation of E-field in z direction can be canceled using a suitable load shifting. The study concerns the influence of the phase shifting when several TLs are considered.

Fig. 3 gives the E-field computed inside the TLEC when two kinds of TLs are considered: a TL with a thin wire conductor and a stripline with a plate conductor as it is used in canonical TEM chambers. It appears that none of the TL induces a uniform field; however the variation of the E-field magnitude inside the working area is greater in the case of the TL with a thin wire conductor (30 dB) compared to the stripline (15 dB). This result should prove that the stripline is a better candidate.

However when one takes into account some practical considerations in the realization of the lines, the TL with wire conductor remains an attractive solution. Then, to improve the E-field uniformity, a solution consists in using several TLs. Fig. 4 illustrates when currents with I 1 =I 2 =1 and I 3 =I 4 =-1 flows along 4 TLs inside the TLEC; as one can see, it still remains minima and maxima of the E-field. A better uniformity of the E-field can then be achieved by introducing a random variation of the excitations. The excitation shifting can be applied to the amplitude or to the phase of the different currents. Fig. 5 shows that the effect on the average E-field is the same whenever the phase or the amplitude is supposed to be random. However, due to practical considerations, the phase shifting technique should be retained here. The variation of the E-field magnitude inside the working area is less important in the case of the 4 TLs with phase shifting (10 dB) compared to the stripline (15 dB) studied previously (see Fig. 3b). I. It appears that the number of TLs is an important parameter: the standard deviation decreases of 3 dB when 4 TLs are introduced instead of 2 TLs. Moreover, the phase shifting improves of 0.5 dB to 3 dB the performances of the chamber.

B. Numerical results

Table II presents the indexes of merit calculated using CST MWS ® software when the TLs are terminated by a 50 Ω load (the characteristic impedance of the TLs is 250 Ω). As one can see, the results in the WA are closed to the ones given in table I: for example the difference is about 0.5 dB in the case of 4 TLs with random phase shifting. However, there is a difference with the results when the WV is considered: this is due to the standing wave that appears because the TLs are unmatched. However, the difference never exceeds 1 dB because the dimension of the WV in z direction (0.6 m) is not large compared to the wavelength (3.75 m at 80 MHz). The performances in the WV can be improved by a load shifting. The effect of changing the load is presented in Fig. 7.

It appears that the maximum of E-field moves between a load of -j500 Ω and +j500 Ω. Considering that the characteric impedance of the TLs is 250 Ω, one finds that the reflection coefficient Γ is respectively e -j0.93 and e +j0.93 ; going back to equation (7), one expects that the minimum of E-field moves of 0.93÷ k 0 = 0.55 m between both loads. This result is illustrated by the simulations reported Fig. 7 (remind that the dimension of the chamber in z direction is 1.2 m, that is to say approximatively 2×0.55 m). Table III presents the indexes of merit calculated when the load shifting is applied to the TLs. As one can see, the performances of the TLEC is sligthly improved compared to the results obtained without the load shifting. Moreover, the results between WA and WV are very close: for example the difference is 0.1 to 0.2 dB in the case of a single TL.

CONCLUSION

It has been shown that a random profile of the standing waves in a TLEC can be obtained by an appropriate phase and load shifting of the TLs, which can be performed electronically. Then, when a TLEC is well designed and operated, it can reach and overcome the uniformity standards stated to canonical chambers.

Figure 1 .

 1 Figure 1. Geometry of a TL in the transverse plane

Figure 2 .

 2 Figure 2. WV and WA definitions in the TLEC

Figure 3 .

 3 Figure 3. Distribution of the normalized E-field in the transverse plane inside the chamber (on the left) and inside the WA (on the right) in dB: (a) TL with a thin wire conductor (b) stripline used in canonical TEM chambers

Figure 4 .Figure 5 .

 45 Figure 4. Distribution of the normalized E-field inside the chamber with 4TL excited by I 1 =I 2 =1 and I 3 =I 4 =-1 Basically, it is possible to move the minima and maxima of the total E-field inside the working area by changing the excitations I 1 , I 2 , I 3 and I 4 of the 4 TLs. For a sake of illustration, two examples are given Fig.5.

Figure 6 .

 6 Figure 6. Distribution of the mean normalized E-field inside the working area with an amplitude shifting (on the left) or a phase shifting (on the right)

Figure 7

 7 Figure 7. . Distribution of the E-field modulus inside the chamber for differents loads (a) -j500 Ω load (b) +j500 Ω load

TABLE I .

 I STANDARD DEVIATIONS FOR SEMI-ANALYTICAL APPROACH

	Chamber Configuration	Standard Deviation (dB)
	Number of TL	Phase shifting	ˆx σ	ˆy σ	, ˆx y σ
	1	-------	5.7	6.0	5.9
	2	No	6.5	3.3	5.3
	2	Yes Random	5.2	4.3	4.8
	4	No	6.4	2.0	5.0
	4	Yes Random	3.3	3.3	3.4
	TEM Chamber	-------	5.1	3.1	4.6

TABLE II

 II

	.	STANDARD DEVIATIONS FOR NUMERICAL APPROACH (LOAD
			50Ω)
	Chamber Configuration	Standard Deviation (dB)
	Number of TL	Phase shifting	Working region	ˆx σ	ˆy σ	, ˆx y σ
	1	----	WA WV	6.1 7.3	6.1 7.3	6.1 7.3
	2	No	WA WV	7.7 7.8	3.6 3.8	4.9 5.1
	2	Yes Random	WA WV	5.9 6.0	4.0 4.2	4.8 5.0
	4	No	WA WV	6.5 6.6	2.7 2.9	3.8 4.0
	4	Yes Random	WA WV	3.9 4.0	3.9 4.1	3.9 4.1

TABLE III .

 III STANDARD DEVIATIONS FOR NUMERICAL APPROACH (LOAD SHIFTING)

		Chamber Configuration	Standard Deviation (dB)
	Number of TL	Load and Phase shifting	Working region	ˆx σ	ˆy σ	, ˆx y σ
	1	Yes, random	WA WV	6.1 6.2	6.0 6.2	6.0 6.2
	2	Yes, random	WA WV	5.5 5.7	4.0 4.2	4.8 4.4
	4	Yes, random	WA WV	3.5 3.7	3.4 3.7	3.5 3.7

This work was partially supported by Capes-Cofecub (07/0568), FAPESP (2007/51192-6) and Brazilian Navy.


Dionysia Triantafyllopoulou

Nikos Passas
email: passas@di.uoa.gr

Evangelos Zervas
email: zervas@ee.teiath.gr

Lazaros Merakos
email: merakos@di.uoa.gr

Optimized Cross-Layer Adaptation for Improved Resource Management in Wireless Networks

The theoretical analysis of a cross-layer mechanism for improving the quality of service of real-time applications in wireless networks is presented. The mechanism coordinates adaptations of the modulation order at the Physical layer and the media encoding mode at the Application layer, to improve packet loss rate, throughput and mean delay. With the use of Continuous Flow Modeling, the system is considered as a "fluid" queue with inflow and outflow rates representing its traffic generation and service rates, respectively. Each data source is modeled as a Markov chain, from the steady-state of which the optimal adaptation thresholds of the cross-layer mechanism are derived. Performance evaluation results show that the optimized operation of the mechanism attains a significant performance improvement compared to both the sub-optimal mechanism, and a legacy system.

I. INTRODUCTION

Recently, cross-layer designs for the performance improvement of wireless communication systems have attracted an increasing research interest. Most proposals in this area coordinate the operation of various adaptation mechanisms residing in different layers of the wireless protocol stack to attain efficient resource utilization and better Quality of Service (QoS) provision. This requires proper adaptation decisions based on information regarding the overall system conditions. In the recent bibliography, many cross-layer design proposals deal with the appropriate adaptation of various operational parameters based on the values of dynamic [START_REF] Choi | Collision-aware design of rate adaptation for multi-rate 802.11 WLANs[END_REF] or static [START_REF] Hoang | Cross-layer adaptive transmission with incomplete system state information[END_REF], [START_REF] Triantafyllopoulou | A Heuristic Cross-Layer Mechanism for Real-Time Traffic over IEEE 802.16 Networks[END_REF] thresholds.

In this paper, we focus on a cross-layer mechanism, introduced in [START_REF] Triantafyllopoulou | A Heuristic Cross-Layer Mechanism for Real-Time Traffic over IEEE 802.16 Networks[END_REF], for the performance improvement of real-time applications in mobile WiMAX networks [START_REF]IEEE Standard for Local and Metropolitan Area Networks Part 16: Air Interface for Fixed and Mobile Broadband Wireless Access Systems. Amendment 2: Physical and Medium Access Control Layers for Combined Fixed and Mobile Operation in Licensed Bands and Corrigendum 1[END_REF]. Although this mechanism was introduced for mobile WiMAX, it can be applied to any wireless network that provides similar adaptation and signaling capabilities. The said mechanism, referred to as Cross-Layer Encoding and Modulation Adaptation (CLEMA), coordinates the adaptive modulation capability at the Physical layer and the multi-encoding mode feature at the Application layer to avoid inefficiencies caused by their independent operation. However, in [START_REF] Triantafyllopoulou | A Heuristic Cross-Layer Mechanism for Real-Time Traffic over IEEE 802.16 Networks[END_REF], the values of the thresholds used by the CLEMA decision algorithm for activating modulation order and encoding mode adaptations are statically determined after extensive simulation experiments. To overcome this inefficiency and derive optimal threshold values, this paper proposes a Continuous Flow Model (CFM) [START_REF] Wardi | Variational Bounds and Sensitivity Analysis of Traffic Processes in Continuous Flow Models[END_REF], [START_REF] Wardi | Loss Volume in Continuous Flow Models: Fast Simulation and Sensitivity Analysis via IPA[END_REF] consisting of a "fluid" queue whose inflow and outflow processes represent the traffic generation and service processes, respectively. The optimal values of the thresholds that activate adaptations of the modulation order or the encoding mode are the ones that lead to optimal system performance in terms of packet loss rate and mean delay. Performance evaluation results show that using optimal values, significant improvement can be attained compared to both a legacy system that adapts the modulation order and encoding mode separately and independently of each other, and CLEMA using static, sub-optimal values.

The rest of this paper is organized as follows. Section II provides an overview of the CLEMA mechanism operation. Section III describes the theoretical analysis of the proposed system using CFM. In Section IV, the system model is described in detail, leading to the derivation of optimal threshold values. Section V presents comparative results that reveal the performance improvement attained with optimal CLEMA. Finally, Section VI contains conclusions and plans for future work.

II. CLEMA: CROSS-LAYER ENCODING AND MODULATION ADAPTATION

The system under consideration is a cellular wireless network, where each cell consists of a Base Station (BS) and a number of Mobile Stations (MSs) within the BS area of coverage. CLEMA is split into two parts, namely the BS part and the MS part, residing at the BS and each MS, respectively. The operation of CLEMA (Figure 1) for both the uplink and the downlink directions can be divided in three main phases. information regarding the performance status of each of its active connections. This information includes channel state conditions in the uplink and downlink directions, packet timeout rate and mean delay. In Phase 2, the BS part uses the collected information to run a decision algorithm (described below) and decide between modifying the encoding mode or using a different burst profile, which is defined as the modulation and coding scheme used. In Phase 3, the various system components are informed of the BS part decisions and perform the required adjustments.

Assuming N active connections in a cell, the QoS parameters used by the CLEMA decision algorithm are the following: i. R err,i is the packet error rate of the i th connection, i ∈ {1, 2, . . . , N}, i.e., the percentage of packets lost due to channel errors. ii. R timeout,i is the packet timeout rate of the i th connection, i.e., the percentage of packets that were lost due to deadline expiration (packet delay has exceeded maximum acceptable delay). iii. R loss,i = R timeout,i +R err,i •(1 -R timeout,i ) is the total packet loss rate of the i th connection. iv. S i is the mean delay of the i th connection. v. μ i is the encoding mode of the i th connection. μ i ∈ {1, 2, . . . , D}, depending on the corresponding multi-rate application. vi. ε maxi and S maxi are the maximum tolerable loss rate and maximum acceptable delay of the i th connection, respectively, while ε mini is a parameter indicating very low loss rate (low loss rate indicator). vii. δ i = Rerr,i R loss,i is the fraction of lost packets (whose loss was due to channel errors) for the i th connection. δ i gives information on the nature of the packet losses experienced by a connection, which is valuable for the adaptation decision process. viii. δ low , δ med and β are the thresholds upon which the algorithm decides on the appropriate adaptation actions. The algorithm is initiated at regular time instants and takes proper adaptation decisions based on the packet loss rate of each connection i (see (D-1) in Figure 2). The actions to be taken depend on the nature of these losses, as follows: 1) In case δ i > δ med (D-2), most of the losses are due to packet errors caused by unfavorable channel conditions, so the most appropriate action is to switch to a more robust burst profile. 2) In case δ i < δ med , most of the losses are the result of delays that cause packet timeouts. The action to be performed depends on the contribution of these timeouts to the overall packet losses: i. If δ i < δ low (D-3), the overall loss rate is caused almost exclusively by packet timeouts and the BS part instructs for an increase of the burst profile. ii. If δ low < δ i < δ med , a significant percentage of packet losses is caused by errors due to poor channel conditions, and the BS part instructs the MS part for a media encoding rate reduction. To achieve an efficient performance under all possible conditions, the algorithm must make adaptation decisions also when the QoS for a specific connection is improved. Thus, when the loss rate decreases significantly, R loss,i < ε mini (D-4), 1) If Si Smax i > β i (D-5), the mean delay is close to the connection delay bound. The algorithm instructs for a burst profile upgrade, which will immediately increase the data rate, and reduce the mean delay.

2) Otherwise, if the mean delay is relatively low compared to the delay bound ( Si Smax i < β i ), the algorithm instructs for a media encoding rate increase to improve the QoS provided to the user.

III. ANALYSIS USING A CONTINUOUS FLOW MODEL

A. General Continuous Flow Model

The operation of the system that employs the CLEMA mechanism can be modeled as a CFM [START_REF] Wardi | Variational Bounds and Sensitivity Analysis of Traffic Processes in Continuous Flow Models[END_REF], [START_REF] Wardi | Loss Volume in Continuous Flow Models: Fast Simulation and Sensitivity Analysis via IPA[END_REF], such as the one depicted in Figure 3. This model consists of a "fluid" queue whose inflow and outflow processes are characterized by flow rates, while its content is defined by the volume of the stored fluid. The size of the queue is finite. Thus, in case the queue is full, the excess flow cannot be absorbed, leading to overflow. The basic storage unit of this model is an Orthogonal Frequency Division Multiple Access (OFDMA) symbol. In the rest of this paper, for presentation purposes, an OFDMA symbol will be referred to as "symbol". Each connection is fed by a data source producing a continuous data flow with maximum tolerable delay S max .

The basic parameters of the CFM are the following: i. α (t): Inflow rate in symbols/s. ii. ϕ (t): Constant service rate in symbols/s (the MAC layer time frame is considered to have a duration of T f seconds and serve c symbols). iii. C = c • Smax T f : Queue size in symbols. Its value is such that overflow occurs when the delay of a queued symbol exceeds the maximum tolerable value S max .

i (t) = μ i,d (t) bi,m(t) while the total inflow rate is α (t) = N i=1 α i (t). The buffer load x (t) is described by the differential equation dx (t) dt + = ⎧ ⎨ ⎩ 0, if x (t) = 0 and α (t) ≤ ϕ (t) 0, if x (t) = C and α (t) -ϕ (t) ≥ 0 α (t) -ϕ (t) , else (1)
The outflow rate is defined as

η (t) = ϕ (t) , if x (t) > 0 α (t) , if x (t) = 0 .
The total overflow rate λ (t) is

λ (t) = α (t) -ϕ (t) , if x (t) = C and α (t) -ϕ (t) ≥ 0 0, else ,
while the overflow rate λ i (t) of each source is proportional to the percentage of the source inflow rate with respect to the total system inflow rate. The total overflow volume during the time interval [0, T ] is L = T 0 λ (τ ) dτ . Thus, the total overflow rate during the time interval [0, T] is R overf low = L T 0 α(τ )dτ . This corresponds to the connection packet timeout rate R timeout defined in Section II.

B. Transmission over the wireless channel

The data that are successfully served by the queue are transmitted over the wireless medium. The data loss rate owing to errors over the wireless channel of the i th source, assuming independent and uniformly distributed bit errors, is (t) where η i (t) is the outflow rate, and BER i (t) is the Bit Error Rate (BER) of the i th source. BER i (t) is a function of the Signal to Noise Ratio (SNR) γ i (t) and the modulation order b i,m (t):

ζ i (t) = η i (t) • 1 -(1 -BER i (t)) bi,m
BER i (t) = f (γ i (t) , b i,m (t)).
The total error rate in symbols/s is defined as

ζ (t) = N i=1 ζ i (t)
. The total loss due to wireless channel errors during

[0, T] is E = T 0 ζ (τ ) dτ . The total error rate is R err = E T 0 η(τ )dτ .
Finally, the total loss rate of the described

system is R loss = R overf low + R err • (1 -R overf low
). The decision algorithm, executed at the BS part of CLEMA, has to determine the appropriate values of μ i,d (t) and b i,m (t) in order to minimize R loss .

C. Calculation of loss and mean delay

In this section, the main QoS parameters based on which the CLEMA decision algorithm adjusts the data source modulation order and encoding mode are calculated using the previously described CFM.

Let us first consider a system that consists of a single data source, the fluid queue and the wireless medium. The source generates data with an encoding rate μ d (t) bits/s, d ∈ {1, 2, ..., D} and modulation order b m (t) bits/symbol, m ∈ {1, 2, ..., M }. The data source mean delay is denoted as S. It is assumed that during the observation interval [0, T] the functions α (t) and ϕ (t) are piecewise constant with a finite number of "jumps". In this case, the function x (t) is piecewise linear, while the functions η (t) and λ (t) are piecewise constant as well. The decision algorithm has to determine the appropriate values of μ d (t) and b m (t), thus to adapt the value of α (t) every ΔT seconds (where ΔT is a constant). If the quality of the wireless channel is assumed to be constant during a time frame, then ΔT = T f . Thus, the

[0, T] interval is divided in K = T T f sub-intervals of T f seconds: [T k , T k+1 ), k ∈ {0, 1, ..., K -1},
x k+1 = min {max {x k + [α k -ϕ k ] • T f , 0} , C} , (2) while the overflow volume during [T k , T k+1 ) is L k = [α k -ϕ k ] • T f + x k -C, if x k+1 = C 0, else .
Thus, in [T k , T k+1 ) the overflow rate is R overf low,k = L k α k •T f . The error rate during [T k , T k+1 ), assuming constant wireless channel quality during a time frame, is R err,k = 1 - (1 -BER k )
bm k where BER k is the connection BER during the time frame. Finally, the total loss rate during [T k , T k+1 ), according to which the decision algorithm will instruct for proper adaptations of the modulation order or the encoding mode, is as follows:

R loss,k = R overf low,k + R err,k • (1 -R overf low,k ) = L k α k • T f + 1 -(1 -BER k ) bm k • 1 - L k α k • T f . (3)
2) Calculation of the mean delay: The value of the data source mean delay, Sk , during [T k , T k+1 ) is derived from the mean buffer load (given the constant service rate of the queue). More specifically, in the time interval [T k , T k+1 ) the mean buffer load is xk = 1

T f T k+1 T k
x k (τ ) dτ . Thus, for the mean

delay we have Sk = xk • T f c .

D. Wireless Channel Model

To model how the BER affects the data source loss rate we need a more specific wireless channel model. The connection is assumed to have a constant mean SNR value, denoted by γ. According to [START_REF] Liu | Queuing With Adaptive Modulation and Coding Over Wireless Links: Cross-Layer Analysis and Design[END_REF], the instantaneous SNR, referred to as γ, is statistically described by the Nakagamim channel model. According to this model, γ is a random variable following the Gamma distribution with probability density function (pdf) t) dt is the Gamma function, and m is the fading parameter of the Nakagami-m channel model, m ≥ 1/2. In this paper, we use m = 1, i.e., we assume a Rayleigh fading wireless channel.

p γ (γ) = m m γ m-1 γm Γ(m) exp -mγ γ where Γ (m) = ∞ 0 t m-1 exp (-

IV. SYSTEM MODEL

Consider a single source CFM, with the source generating data with encoding rate μ d (t) bits/s and modulation order b m (t) bits/symbol. The maximum tolerable loss rate is ε max while the maximum tolerable waiting time in the buffer is S max . The δ parameter, which is used by the decision algorithm and is the basic criterion for the determination of the most appropriate modulation order and encoding mode, is defined as the ratio of the error rate over the total loss rate: δ = Rerr R loss . This single source system can be modeled as a Markov chain whose state is the triplet (d, m, x), where d is the encoding mode, m is the modulation order, and x is the discretized buffer load. For the state space E of this chain we have

E = {(d, m, x) : d ∈ {1, 2, ..., D}, m ∈ {1, 2, ..., M }, x ∈ {0, 1, ..., C} }.
This model can be extended to a symmetric multiple source system as follows. Consider N identical data sources of the same traffic class. If the size of the queue equals C symbols, it can be assumed that C N symbols are dedicated to each source. Each source is then modeled as an independent Markov chain of D × M × ( C N + 1) states.

A. Calculation of the Markov Chain transition probabilities

Since adaptations of the modulation order and the encoding mode are only allowed between successive levels, we have

P ((d, m, x) , (d , m , x )) = 0, if |d-d | > 1 or |m-m | > 1.
According to the CLEMA decision algorithm, the state transitions occur as follows:

1) Transition to a state with lower modulation order:

P ((d, m, x) , (d, m -1, x )) = P (R loss > ε max , δ > δ med ) , m > 1,
where x is the buffer load after the transition to the new state.

2) Transition to a state with higher modulation order:

P ((d, m, x) , (d, m + 1, x )) = P (R loss > ε max , δ < δ low ∪ R loss < ε min , S S max > β , m < M.
3) Transition to a state with lower encoding mode:

P ((d, m, x) , (d -1, m, x )) = P (R loss > ε max , δ low < δ < δ med ) , d > 1.
4) Transition to a state with higher encoding mode:

P ((d, m, x) , (d + 1, m, x )) = P R loss < ε min , S S max < β , d < D.
In any other case, only transition from state (d, m, x) to state (d, m, x ) occurs.

Using the previously described CFM, the probability of transition to a state with lower modulation order, i.e., (

d i , m i , x i ) → (d i , m i -1, x j ), m i > 1, i, j ∈ E can be calculated. According to (2), x j = min max x i + μ d i bm i -1 -ϕ • T f , 0 , C and x i = min x i • bm i bm i -1 , C
. The probability of this transition is

P (i,j)=P (R loss >εmax, δ>δ med ) = P L i α i •T f +(1-(1-BERi) bm i )• 1-L i α i •T f >εmax , 1-(1-BER i ) bm i L i α i •T f + ( 1-(1-BER i ) bm i ) • 1-L i α i •T f >δ med ⎞ ⎠ = P ⎛ ⎜ ⎝BERi>1- ⎛ ⎝ 1-εmax 1- L i α i •T f ⎞ ⎠ 1 bm i , BERi>1- ⎛ ⎝ 1- L i α i •T f •δ med 1-δ med • 1- L i α i •T f ⎞ ⎠ 1 bm i ⎞ ⎟ ⎠ = P ⎛ ⎜ ⎝BERi>max ⎧ ⎪ ⎨ ⎪ ⎩ 1- ⎛ ⎝ 1-εmax 1- L i α i •T f ⎞ ⎠ 1 bm i , 1- ⎛ ⎝ 1- L i α i •T f •δ med 1-δ med • 1- L i α i •T f ⎞ ⎠ 1 bm i ⎫ ⎪ ⎬ ⎪ ⎭ ⎞ ⎟ ⎠.
The probabilities of all other possible transitions are derived in a similar way.

B. Calculation of the mean system loss

To evaluate the performance of the system that employs the CLEMA mechanism, the mean overflow and the mean error rate have to be calculated. Let the values of the thresholds δ low , δ med and β be such that the Markov chain does not contain any absorbing states, and let E c be its minimum closed set. Since E c has a finite number of states, the system πP = π, π1 = 1, where P is the transition matrix of E c , has a strictly unique solution and π (j) = lim n→∞ P n (i, j), ∀ i, j ∈ E c [START_REF] Cinlar | Introduction to Stochastic Processes[END_REF].

For the system mean loss rate we have

R loss = i∈Ec π (i) • R loss,i where R loss,i = R overf low,i + R err,i • (1 -R overf low,i
) is the loss rate in state i, according to (3).

Let the system be in state i ∈ E c , with inflow rate α i symbols/s, modulation order b mi bits/symbol, overflow volume L i symbols, buffer load x i symbols, and mean bit error rate BER i = f ( γi , b mi ).

Apart from the loss that is a result of overflow due to excess inflow rate, additional overflow can occur in case the system performs a transition to a state with lower modulation order. More specifically, due to the fact that the decrease of the modulation order results in an increase of the buffer load, it is possible that overflow occurs due to such a transition. This loss is calculated for transitions to state i ∈ E c from each state j ∈ E c such that j = i, b mj > b mi , and is equal to L transition (j, i) = max x j bm j bm i -C, 0 . Thus, the total data loss (in bits) in state i due to overflow, which is the result of modulation order decreases, is L transition (i) = j∈Ec, bm j >bm i P (j, i)•L transition (j, i)•b mi . Thus, the system mean loss volume due to overflow caused by modulation order decreases is

L transition = i∈Ec L transition (i) • π (i).

C. The Legacy System

The performance of the system employing CLEMA is compared against a typical system, referred to as "legacy" system, which performs the operations of adaptive modulation and encoding mode adjustment separately and independently of each other. More specifically, in the legacy system a decrease of the modulation order is performed if the current value of the BER is higher than the maximum tolerable BER, BER max , while an increase of the modulation order is performed if the BER is lower than the minimum BER, BER min . Similarly, a transition to a state with lower encoding mode is performed if the total loss rate exceeds the maximum tolerable loss rate ε max , while a transition to a state with higher encoding mode is performed if the loss rate is lower than ε min . The modeling of the legacy system is analogous to that of the system that employs CLEMA using a CFM, while a data source is modeled as a Markov chain of D × M × (C + 1) states. Using a similar methodology, the transition probabilities between the various states of the Markov chain that models the legacy system can be calculated, and the mean delay as well as the additional losses caused by modulation order adaptations can be derived.

D. Derivation of the optimal adaptation thresholds

The optimal values δ lowopt , δ medopt , and β opt of the thresholds δ low , δ med and β, used by the CLEMA decision algorithm are the ones that minimize the system loss rate R loss , while maintaining the system's mean error rate R errCLEMA lower than or at most equal to the error rate of the respective legacy system R errLeg .

V. PERFORMANCE EVALUATION

With the use of Matlab, the performance of the system employing CLEMA with optimal values of the thresholds δ low , δ med and β (referred to as "optimal CLEMA") was compared against i) the legacy system, and ii) a system that employs CLEMA with constant sub-optimal values of the adaptation thresholds (referred to as "sub-optimal CLEMA"). The performance metrics used for the comparison of the three systems were the mean loss rate R loss , and the mean delay S.

The evaluation scenario considers a data source that produces Adaptive Multi-Rate (AMR) [START_REF]3GPP; Technical Specification Group Services and System Aspects; Mandatory speech CODEC speech processing functions[END_REF] voice traffic with five discrete encoding modes. The modulation orders used are QPSK, 16-QAM and 64-QAM, which correspond to 2, 4 and 6 bits/symbol, respectively. The duration of the MAC layer time frame, which is the period used for the adaptation decisions, is 5 ms. The maximum tolerable waiting time in the queue S max is 40 ms. The maximum tolerable loss rate ε max is equal to 7 • 10 -3 , while the value of ε min is 10 -9 . The values of the thresholds based on which the legacy system performs the adaptations of the modulation order are: BER max = 10 -6 and BER min = 10 -9 . Figure 4 depicts the mean overflow rate, R overf low , and the mean error rate R err of the three systems (optimal CLEMA, sub-optimal CLEMA and legacy) versus the mean SNR (γ) of the data source. The service rate per traffic frame equals c = 10 symbols/5 ms, which corresponds to a queue size of C = 80 symbols. All three systems manage to maintain a relatively low mean error rate, which generally decreases as the mean SNR increases. This is a result of the fact that, as the conditions of the wireless medium improve, the systems have the ability to increase the system data rate by instructing more frequent transitions to states with higher modulation orders. Nevertheless, the legacy system is characterized by a relatively high overflow rate, as a result of the independent operation of the adaptive modulation at the Physical layer and the encoding mode adaptations at the Application layer. On the contrary, optimal CLEMA results in a significantly lower overflow rate, especially in cases of increased mean SNR, taking advantage of the coordinated operation and the optimal values.

The importance of the dynamic adaptation of the values of the thresholds δ low , δ med and β is depicted from the performance of sub-optimal CLEMA, whose threshold values are constant and optimal only for the case of γ=30 dB, as shown in Figure 4. Although this system results in improved performance compared to the legacy system, it does not manage to obtain the minimum values of the mean overflow rate of optimal CLEMA for all values of the mean SNR. Figure 5 depicts the data loss L transition of the legacy, optimal CLEMA and sub-optimal CLEMA systems caused by transitions in states with lower modulation order that, depending on the queue load, may result in data overflow. In the three systems under consideration, L transition follows a declining course as the mean SNR increases. This is a result of the fact that, as the quality of the wireless medium improves, transitions to states with lower modulation orders are instructed less frequently. Thus, in such cases, the volume of the data loss is lower. As shown in this figure, the operation of adaptive modulation in the legacy system, which is performed based exclusively on the mean SNR, leads to more frequent reductions of the modulation order compared to the systems that employ the optimal and suboptimal CLEMA mechanisms, resulting in significantly higher data loss. Both optimal and sub-optimal CLEMA have similar and significantly improved performance in terms of additional data losses compared to the legacy system, as this depends mainly on the operation of the CLEMA decision algorithm.

Figure 6 plots the mean delay versus the mean SNR for the three systems under comparison. Clearly, optimal CLEMA has significantly improved performance in terms of mean delay compared to the legacy and sub-optimal CLEMA systems, especially in cases of favorable wireless channel conditions.

Finally, Figure 7 depicts the total loss rate of the legacy system and the system employing the CLEMA mechanism as a function of the thresholds δ low and δ med , assuming a mean SNR equal to 20 dB, and a service rate per traffic frame c=20 symbols/5ms, which corresponds to a queue size of C = 160 symbols. Additionally, for the CLEMA mechanism the value of the threshold β equals 0.3. As expected, the loss rate of the legacy system has a constant value of 1.21 • 10 -2 , as its performance is not affected by the values of the thresholds. On the contrary, the loss rate of CLEMA decreases as the value of δ low decreases and the value of δ med increases, with the optimal performance obtained when δ low =0 and δ med =1. This means that, in this case, optimal CLEMA instructs an increase of the modulation order only when no errors occur at the wireless channel, and a decrease of the modulation order only when the queue experiences no overflow.

VI. CONCLUSIONS

In this paper, we have introduced a Continuous Flow Model for the analysis and optimization of the CLEMA mechanism proposed in [START_REF] Triantafyllopoulou | A Heuristic Cross-Layer Mechanism for Real-Time Traffic over IEEE 802.16 Networks[END_REF]. The model consists of a "fluid" queue whose inflow and outflow rates represent the system's traffic generation and service rates, respectively. Each data source is modeled as a discrete-time Markov chain, with each state representing its current encoding mode, modulation order and respective buffer load. The Markov chain's transition probabilities are derived from the CLEMA decision algorithm. The system performance is assessed in terms of overflow rate, error rate and mean delay. Performance evaluation results show that CLEMA with the optimal values for the adaptation thresholds, derived by the proposed CFM, significantly outperforms both a legacy system that performs the adaptations of the modulation order and encoding mode separately and independently of each other, and a system employing sub-optimal CLEMA using static values of its adaptation thresholds.

Fig. 1 .

 1 Fig. 1. The CLEMA mechanism for the downlink

Fig. 2 .

 2 Fig. 2. CLEMA decision algorithm flow chart

Fig. 3 .

 3 Fig. 3. Multiple source CFM

  where T 0 = 0 and T K = T , during each of which the value of the function α (t) is constant. The value of ϕ (t) is by definition constant during all these time intervals. 1) Calculation of the loss rate: Let a k , ϕ k , η k and b m k be the values of a (t), ϕ (t), η (t) and b m (t), respectively, during the time interval [T k , T k+1 ). Similarly, x k = x (T k ) is the buffer load at T k and L k = T k+1 T k λ (τ ) dτ is the overflow volume during [T k , T k+1 ). The buffer load at T k+1 is

Fig. 4 .

 4 Fig. 4. Rerr and R overf low vs. γ

Fig. 5 .

 5 Fig. 5. L transition vs. γ (service rate per traffic frame c=10 symbols/5ms)

Fig. 7 .

 7 Fig. 6. Mean delay vs. γ