F Dell'isola

M Guarascio

K Hutter

A variational approach for the deformation of a saturated porous solid. A second-gradient theory extending Terzaghi's effective stress principle

Keywords: Principle of virtual power, second-gradient theory, saturated soil-water mixture, salt rock, porosity, static permeability

The principle of virtual power is used to derive the equilibrium ®eld equations of a porous solid saturated with a ¯uid, including second density-gradient effects; the intention is the elucidation and extension of the effective stress principle of Terzaghi and Fillunger. In the context of a ®rst density-gradient theory for a saturated solid we interpret the porewater pressure as a Lagrange multiplier in the expression for the deformation energy, assuring that the saturation constraint is veri®ed. We prove that this saturation pressure is distributed among the constituents according to their respective volume fraction (Delesse law) only if they are both true density-preserving. We generalize the Delesse law to the case of compressible constituents. If a material-dependent effective stress contribution is to arise, it is, in general, nonvanishing simultaneously in both the solid and ¯uid constituents. Moreover, saturation pressure, effective stresses and compressibility constitutive equations determine the exchange volume forces. In a theoretical formulation without non-isotropic strain measures, second density-gradient effects must be incorporated, not only to accommodate for the equilibriumsolid-shear stress and the interaction among neighboring solid-matrix pores, but also to describe internal capillarity effects. The earlier are accounted for by a dependence of the thermodynamic energy upon the density-gradient of the solid, while the latter derives from a mixed density-gradient dependence. Examples illustrate the necessity of these higher gradient effects for properly posed boundary value problems describing the mechanical behaviour of the disturbed rock zone surrounding salt caverns. In particular, we show that solid second-gradient strains allow for the de®nition of the concept of static permeability, which is distinct from the dynamic Darcy permeability.

Introduction

A saturated porous material is a mixture of a solid matrix with pores which are completely ®lled with a ¯uid, called the porewater ¯uid. The state of stress in a material point of this mixture consists of two contributions; one contribution is the constraint pressure p due to saturation, called in modern porous media theories the saturation pressure, in the soil science literature, however, better known as the porewater pressure. This pressure is distributed between the two constituents: in older theories via a postulate (Delesse law, s. e.g. [START_REF] De Boer | The development of the concept of effective stresses[END_REF]), according to volume fraction occupied by the respective constituents, p a pm a , m a is the volume fraction of constituent a; in newer theories by derivation. The other contribution to the stresses is constitutive, and thus responsible for the strength of the material. It is customary in the solid mechanics literature to refer to these second types of stresses as the effective stresses (note plural) in the solid and ¯uid, respectively. In anticipation that the ¯uid is true density-preserving, the solid carries all constitutive properties, and the term effective stress is only used in the singular form referring to the solid matrix.

This understanding was that of Terzaghi and Fillunger, and the latter was the ®rst soil engineer who clearly stated that a constitutive equation must be formulated only for the excess over the weighted porewater pressure and not for the total stress, s. [START_REF] De Boer | The development of the concept of effective stresses[END_REF] p. 81±83. However, also Terzaghi was of the opinion that the excess over the porewater stress ``has its seat exclusively in the solid phase of the soil'' (according to [START_REF] De Boer | The development of the concept of effective stresses[END_REF] p. 83, quoting Skempton 1960). Moreover, both Fillunger and Terzaghi postulate that solid volume fraction kinematically describes the strain of the solid matrix.

According to the above statements three questions are raised and will be answered in this paper: [START_REF] De Boer | The development of the concept of effective stresses[END_REF] Is the division of the porewater pressure between the solid and ¯uid constituents truly according to ``pressure equilibrium'', p a pm a ; or can it be different? [START_REF] Germain | La me Âthode des puissances virtuelles en me Âcanique des milieux continus. Premie Á partie: The Âorie du second gradient[END_REF] Is the concept of effective stress one that must be restricted to the solid matrix as expected by Terzaghi or does it apply for the solid and the ¯uid? (3) Is the solid volume fraction suf®cient to describe the solid volume pore strain?

To answer these questions we will use the principle of virtual power applied to a binary mixture of a solid saturated with a ¯uid. In order to simplify the required mathematical formalism we will limit ourselves to consider equilibrium con®gurations of the considered mixture. The ¯exibility in the underlying postulates will be assumed to be suf®ciently broad to cover the conjectured properties. We will assume, therefore, that the two constituents may be either compressible or true density-preserving, i.e. both situations can be imposed by specialization. Second, we shall assume the thermodynamic energy to depend both on the partial densities of the constituents as well as on their gradients. Third, by postulating the thermodynamic energy to depend on the density gradients an incompleteness of the solid matrix strain description will be removed. The latter dependence will give rise to double gradients of the densities in the constitutive relations for the pressure, which is why dependence on the density gradients in the energy is called a second-gradient theory. Gradients higher than the ®rst are e.g. needed to transmit nonspherical equilibrium stresses. In the application developed in the last Section we show that second-gradient strains of the solid matrix allow the de®nition of a static permeability, i.e. of a physical quantity measuring how under equilibrium conditions (in the absence of ¯uid motion) a porous solid matrix offers resistance to the variation of the saturating ¯uid volume fraction. Other dependencies could be incorporated but the intention is to keep the theory as simple as possible.

The ®rst two questions raised above and modifying the Terzaghi-Fillunger conjectures can be answered with a ®rst-gradient theory. It will be shown that, within the context of this variational formulation, the porewater pressure is, in general, not distributed between the constituents according to pressure equilibrium, as postulated by the Delesse law, and effective stresses arise for both the solid and the ¯uid constituents. Moreover, in order that the equilibrium equations derived from the variational principle be consistent mixture balances, the constitutive relations for the thermodynamic energy, the effective stresses, the constituent compressibilities and exchange volume force between constituents must be related by imposing objectivity of the power expended by internal contact actions [START_REF] Germain | La me Âthode des puissances virtuelles en me Âcanique des milieux continus. Premie Á partie: The Âorie du second gradient[END_REF][START_REF] Germain | The method of virtual power in continuum mechanics Part 2: Microstructure[END_REF]. The resultant relations (s. the following Eqs. (20)±( 22)) can be used in the framework of the present variational formulation to prove the following statements:

When both constituents are true density-preserving, the porewater pressure is distributed among the constituents according to the ``pressure equilibrium'' postulate (Delesse law). On the other hand, when the constituents are compressible, the Delesse law is no longer valid, and the porewater pressure is distributed according to a prescribed formula. Fillunger's [START_REF] Fillunger | Erdbaumechanik Wien: Selbstverlag des Verfassers[END_REF] statement that compressibility leads to a constitutive effective stress is corroborated.

The interaction force has a contribution proportional to the gradient of the solid volume fraction. The coef®cient is, in general, not equal to the porewater pressure, but reduces to the latter where the effective pressures vanish.

Effective stresses, in general, arise simultaneously in both constituents, and are related to the mixture thermodynamic pressure. This statement is also in agreement with other formulations of granular solid-¯uid mixtures such as [START_REF] Svendsen | A thermodynamic model for volume-fraction-based mixtures[END_REF] 1 or [START_REF] Krishnaswamy | A thermomechanical theory of solid-¯uid mixtures[END_REF]. The second gradient effects must be incorporated in a theory ignoring a dependence of the thermodynamic quantities on nonhydrostatic strain measures or rate-independent plastic effects, in order to allow for equilibrium shear stresses. These effects are also needed to describe the effects on stress in a material point related to the solid-matrix-pore deformation in the neighbouring material points.

These results seem to be new. One advantage of the variational formulation over direct approaches to higher gradient theories is the fact that the minimization principle automatically generates the natural boundary conditions. In other words, not only the dynamic equations within the body relating the introduced ®elds are generated, but also the conditions which these ®elds must ful®ll at the boundaries of the body. In ®rst-gradient theories these boundary conditions lend themselves easily through arguments of physics, for higher-gradient theories they must fall out from the formulation, because this physical intuition is generally not available. We adapt the arguments developed by [2, 3, 7±9], in order to physically interpret for the considered system the deduced boundary conditions. Furthermore, we demonstrate their appropriateness by solving a simple one-dimensional boundary value problem. It catches, at least qualitatively, some physical properties, [START_REF] Berest | A tentative classi®cation of salts according to their creep properties SMRI Proceedings Spring Meeting[END_REF], of deformable salt rocks fully saturated by a ¯uid, which cannot be described without the introduction of second density gradients.

Conceptual prerequisites

Consider a binary mixture of a solid matrix with connected pores which are ®lled with a liquid. This arrangement can be thought of as being a saturated soil or rock. We have in mind as an application, however, a dome of salt which is saturated with a brine of given concentration. Solidi®cation of salt from the brine to the solid matrix is excluded. Let the two components be referred to as the solid and the porewater ¯uid and be indicated by the suf®ces s and f . Let, moreover, q s ; q f and v s ; v f be the solid and ¯uid densities and velocities, respectively, in the mixture. The mixture density and the barycentric velocity are then given by q q s q f qs m s qf m f qs m s qf 1 À m s ; 1

v q s q v s q f q v f : n s v s n f v f ; 2
in which m s is the solid volume fraction, m f the porosity; the saturation condition has been used stating that the pore-¯uid ®lls the entire pore space, and qs , qf are the true densities of the solid and the ¯uid. We conceive this mixture to be nonreactive so that the balances of mass for the constituents reduce to oq a ot r Á q a v a 0; a s; f : 3

In the ensuing analysis we shall restrict ourselves to purely mechanical processes; temperature will play no role, and so the constituent momentum equations are the only additional balance laws to be added to [START_REF] Germain | The method of virtual power in continuum mechanics Part 2: Microstructure[END_REF]. Instead of a direct application of these laws, we shall use to derive them the principle of virtual power, applied to the appropriate energy functional. Let w be this functional, and assume it to depend on the ®elds q a and rq a , a s; f , w wq a ; rq a : 4

Other dependences could also be introduced, e.g. a temperature dependence and dependences on tensorial strain measures; however, our interest is in the derivation of the most simple Terzaghi-type effective stress theory including second gradient effects. Imposing objectivity, Eq. ( 4) takes the form w wq s ; q f ; jrq s j 2 ; jrq f j 2 ; jrq s Á rq f j :

This form of the energy functional supposes that both the solid matrix as well as the ¯uid exhibit second-gradient effects, and that the two also give rise to an interaction energy of the two effects through the last variable in [START_REF] Svendsen | On the thermodynamics of a mixture of isotropic materials with constraints[END_REF]. The principle of virtual power states that the variation of the total energy in the body related to its motion equals the power of the external forces, i.e.

d dt B wdV B b s Á v s b f Á v f À Á dV oB C s v s ; rv s C f v f ; rv f dS : 6
Here, b a are the speci®c body forces which perform work on their constituent motions, indicated by the velocity ®elds v s and v f ; respectively, while integrals of C a are scalar-valued bilinear functionals. They represent the power expended on the velocity ®elds and their gradients by speci®c contact actions; this is so because, as in every second-gradient theory (s. e.g. [START_REF] Germain | La me Âthode des puissances virtuelles en me Âcanique des milieux continus. Premie Á partie: The Âorie du second gradient[END_REF][START_REF] Germain | The method of virtual power in continuum mechanics Part 2: Microstructure[END_REF] or [START_REF] Dell'isola | Edge contact forces and quasi-balanced power[END_REF]), contact forces exert power also on gradients of velocity ®elds. The variation of the total energy on the left-hand side of ( 6) is (i) referred to a barycentric motion in the following Sec. 3, and to constituent motions when considering the second gradient terms added in Sec. 5; and

(ii) calculated by independently varying the independent variables of the energy function subject to the kinematic constraint that the body is saturated. The latter condition is now incorporated by de®ning the energy according to

w e p m s m f À 1 À Á e q a ; rq a p q s qs q f qf À 1 2 3
; 7

in which qa are the true peculiar densities, assumed to be constants for density-preserving constituents; note that q a m a qa , and p is a Lagrange multiplier which will be determined by imposing the saturation constraint; the quantity e is the thermodynamic free energy. In the case that entropy is the independent thermodynamic ®eld variable, e is the internal energy. If this independent ®eld variable is the temperature, then e is the Helmholtz free energy. Because we do not specify this, we are either concerned with isentropic or isothermal processes. With the interpretation [START_REF] Seppecher | Etude des conditions aux limites en the Âorie du second-gradient: cas de la capillarite Â[END_REF], the variation in ( 6) can be performed for unconstrained q a . Finally, let us state that we do not require ab initio that the components be incompressible. In fact, it will be supposed that the true densities qa may be affected by the composition of the grains. At high porosity, qa will essentially be constant, at low porosity near the closest packing, qa will itself increase. This dependence was clearly noted by Fillunger, and can be accounted for by postulating a constitutive relation qa qa q a : For the ®nal application, we will assume that qf const. and qs qs q s : 3 The classical model Consider, as a preliminary approach, a theory without second gradient effects. Then The ®rst term in the integrand of ( 11) can be transformed according to

w e q b p q s qs q f qf À 1 2
B Àp a r Á v a dV À B r Á p a v a dV B rp a Á v a dV À oB p a v a Á ndS B rp a Á v a dV : 12
With ( 11) and ( 12), the principle of virtual work [START_REF] Krishnaswamy | A thermomechanical theory of solid-¯uid mixtures[END_REF], takes the form

a oB Àp a n Á v a dS a B r q a ow oq a À ow oq a rq a Á v a dV a B b a Á v a dV a oB C a v a ; 0dS : 13
This must hold for all ®elds v a de®ned over the body B and on its boundary oB, thus leading to Àq a r ow oq a b a 0; a s; f in B ; 14

Àp a n C a v a ; 0 0; a s; f on oB : 15

In order to identify in the last equation the addends representing partial stresses and exchange bulk force between the constituents, we impose the same Galilean invariance argument as done in [START_REF] Germain | La me Âthode des puissances virtuelles en me Âcanique des milieux continus. Premie Á partie: The Âorie du second gradient[END_REF][START_REF] Germain | The method of virtual power in continuum mechanics Part 2: Microstructure[END_REF] and [START_REF] Seppecher | Etude des conditions aux limites en the Âorie du second-gradient: cas de la capillarite Â[END_REF]: it follows that Àrp a m a b a 0; a s; f in B; Àp a n t a 0; a s; f in oB ; 16

in which the exchange bulk forces m a are given by

m s Àm f À 1 À n s rq s w q À ow oq s 1 À n f À Á rq f w q À ow oq f 2 3
; 17 and t a is the linear part of C a v a ; 0 on v a and represents a stress vector. We remark explicitly that when the saturation constraint is not considered, and therefore p in Eq. ( 8) vanish, Eq. ( 17) coincides with Eq. (4.6) in [START_REF] Krishnaswamy | A thermomechanical theory of solid-¯uid mixtures[END_REF].

With q a m a qa , qa qa q a , and w as given by [START_REF] Dell'isola | The relationship between edge contact forces, double forces and interstitial working allowed by the principle of virtual power[END_REF]

P b m b 1 À n b 1 À m b dq b dq b À1 : 20
We cannot resist to point out here the close similarity between Eq. ( 19) and the well-known formula

r eff m r m 1 À k 0 k s P ; 21
of soil mechanics (s. [11±13]), with the correspondence

r eff m 6 ÀP a r m 6 Àp a P 6 pm a ;
where r eff m is the mean effective stress, r m the mean applied stress and P the pore pressure. The factor of the second term on the right-hand side of ( 21) is called the Biot coef®cient, and k 0 ; k s are the moduli of compressibility of the drained soil and the rock material. It is tempting to identify k 0 =k s with m a dq a =dq a ; but this is no more than a suggestion.

The expression [START_REF] Dell'isola | What are the dominant thermomechanical processes in the basal sediment layer of large ice sheets?[END_REF] holds modulo p a m a À 1, which is zero, and the quantity

P a : q a oe oq a À e q 22
may be identi®ed with the thermodynamic pressure of the constituent a. Therefore, the equilibrium equation ( 16) 1 takes the form

1 2 3 Àrp a Mrm a b a 0 : 23
Notice that Eq. ( 23) expresses a balance between the pressure gradient (1), the exchange force (2) and the external body force (3), which necessarily must satisfy the basic postulates of the mixture balance laws. We also mention that, for dq b =dq b 0, the above expression for m a coincides with that found with different methods in [START_REF] Svendsen | A thermodynamic model for volume-fraction-based mixtures[END_REF]. The quantity P b =m b was denoted by b b in [START_REF] Svendsen | On the thermodynamics of a mixture of isotropic materials with constraints[END_REF].

Relation [START_REF] Dell'isola | A qualitative Analysis of the dynamics of sheared and pressurized layer of saturated soil[END_REF] serves as a restriction among constitutive equations: for instance, if qf and qs are both constant, i.e. if the two constituents are true density-preserving, then [START_REF] Dell'isola | A qualitative Analysis of the dynamics of sheared and pressurized layer of saturated soil[END_REF] requires that, once the constitutive quantity e and the pressure p are assigned, so is the exchange force M: In particular, if e 0 and dq b =dq b 0, Eq. ( 23 The ®rst form shows that the partial pressures p a pm a are obtained from the saturation (porewater) pressure by multiplying it with the constituent volume fraction, a property usually referred to as ``pressure equilibrium''. On the other hand, the interaction force is simply given by Àprm a : Such choices are the basis of a large number of porous solid theories (s. e.g.

[14±16]), but this choice is obviously very restrictive.

The above results show that 1. Fillunger [START_REF] Fillunger | Erdbaumechanik Wien: Selbstverlag des Verfassers[END_REF] was essentially correct when conjecturing that an extension of the classical theory exhibited by [START_REF]Fall Meeting[END_REF] can be obtained via compressibility assumptions of the constituents solid and ¯uid, (s. [START_REF] De Boer | The development of the concept of effective stresses[END_REF], p. 82, footnote 8): indeed Eq. ( 22) de®ne the Filllunger±Terzaghi effective partial pressures, and show that they are related to the compressibility constitutive relations qb q b : 2. Both Fillunger and Terzaghi (s. again [START_REF] De Boer | The development of the concept of effective stresses[END_REF]) did not investigate the implications of the balance of energy and of the second principle of the thermodynamics in the theory which they were developing; in particular, they assumed that (using the notation of the present paper) P s T 0; P f 0; M p ; 26 which clearly contradicts via [START_REF] Dell'isola | A qualitative Analysis of the dynamics of sheared and pressurized layer of saturated soil[END_REF] the principle of virtual powers on which we have based our treatment.

Limits of the classical model

By adding over the constituents their force balances [START_REF] Berest | Salt caverns and the compressibility factors[END_REF] an expression for the porewater pressure can be found. Indeed, by assuming that the body forces are conservative and possess a potential where k may be a function of time, which we take to be constant. This result shows that for homogeneous ¯uid-saturated solid matrices and for a vanishing potential ®eld / (i.e. b a 0 for a s; f ), the pressure ®eld p is constant if and only if one of the ®elds q a or m a is constant. On the other hand, it is easy to verify that, under the same hypothesis on /, both Eq. ( 23) are solved by constant ®elds q a , as can be checked immediately with the aid of ( 15). These results hold, no matter what the geometry of the solid matrix may be, and irrespective of whether the body be in equilibrium or in a dynamic Stokesian motion. On the other hand, it is easy to conceive physical situations in which a solid matrix saturated by a ¯uid shows a variable solid pore size distribution even in an equilibrium con®guration, and, consequently, a spatially variable solid-volume fraction. An instance of a physical system in which this occurs is given by the micro-cracked and permeable salt rock in the disturbed zone surrounding a ¯uid ®lled cavity, [START_REF] Pfei¯e | Damage-Induced Permeability Enhancement of Natural Rock Salt with Implications for Cavern Storage[END_REF]. It seems, therefore, useful to formulate a model in which the micro-mechanical interactions among neighbouring pores are described in such a way that the existence of a spatially varying equilibrium solid volume fraction ®eld is permissible in the absence of body forces. This will be done in the next Section, by introducing a second gradient theory.

A second-gradient theory

We choose now the thermodynamic energy in the form e e q a ; f ab e q a ; f ab ; f ab : rq a Á rq b ; 30 and thus need to complement the variation of B wdV performed in Secs. which, through rearrangements of differentiations and the application of Gauss' theorem, may be written as (s. [START_REF] Seppecher | Etude d'une modelisation des zones capillaires ¯uides: interfaces et lignes de contact[END_REF], p. 19

) B r Á A b Àr Á C b Áv b dV À oB v b Á A b Àr Á C b r v b : C b Á ndS : 38
Adding the transformed Eq. (34) to the variation calculated in the previous Section we obtain

a B r q a ow oq a À ow oq a rq a À r Á A a Àr Á C a ! Á v a dV a oB v a Á p a n A a Àr Á C a rv a : C a dS a B v a Á b a dV a oB C a v a ; rv a dS : 39
This identity proves that the expression on the right-hand side is a bilinear functional of the velocity ®eld v a and its gradient. Thus, we write C a v a ; rv a t a Á v a d a : rv a ; 40 and call t a surface traction vector and d a surface double-force tensor of constituent a. Formula (40) is nothing else than Cauchy's theorem extended to the second-gradient mixture theory treated here. Since (39) must be valid for all v a and rv a , we ®nally deduce the ®eld equations r q a ow oq a À ow oq a rq a À r Á A a Àr Á C a b a 0; in B ; 41 and the boundary conditions

Àp a I À A a À r Á C a Á n t a C a Á n d a '
on oB : 42

The surface traction vectors t a and the surface contact double forces d a for the constituent stresses are to be externally prescribed at the boundary oB (with a unit normal n) of a mixture body. In a ®rst-gradient theory, both A a and C a vanish, and Eqs. (42) reduce to [START_REF] Drew | Application of general constitutive principles to the derivation of multidimensional two-phase ¯ow equations[END_REF] where a and l are constitutive parameters describing the dynamic behaviour of the solid-matrix ¯uid interface. Moreover, we assume that

d dp i 47
is a given function of the incumbent pressure, which in ensuing application will be taken as linear

d Dp i :
Finally, it is easy to show that when Eqs. (44) are added we can generalize (28), and obtain

r a 1 À m 2 a dq a dq a p P a ! / @ A Àr Á k s q s Rq s k s 2 rq s j j 2 I À k s rq s rq s ! 0 : 48 6 Linearized constitutive equations
In this Section we consider the constitutive equations for the energy e and true densities qa in the neighbourhood of an equilibrium con®guration q 0 s ; q 0 f which is stress free. Having in mind the application to the study of equilibrium con®gurations of cracked salt rocks saturated with a ¯uid, the constitutive postulate (49) on e can be interpreted as follows: the rock has been damaged and, having an apparent density q 0 s ; can be saturated by a ¯uid of a given true density q0 f and apparent density q 0 f without any change of deformation. Therefore, in what follows the densities q 0 s ; q 0 f are assumed to be constitutive quantities describing the state of damage of the considered rock. Therefore, we assume for e the following expansion:

eq s ; q f 1 2 o 2 e oq 2 s q 0 s ;q 0 f q s À q 0 s À Á 2 1 2 o 2 e oq 2 f q 0 s ;q 0 f q f À q 0 f 2 o 2 e oq f oq s q 0 s ;q 0 f q s À q 0 s q f À q 0 f O3 : 49
A constant and a linear term are omitted to avoid pre-stresses in the natutal reference con®guration. Moreover, the reference state is assumed to be saturated, so that

q 0 s q0 f q 0 f q0 s q0 f q0 s : 50
Note, that in a ¯uid ®lled cavern of a salt-rock formation the applied ¯uid pressure in the cavern generates openings of the grain boundaries between the salt crystallites in the immediate vicinity of the cavern wall. As a result, its permeability and drainage properties change.

For more details see Chapter V of SMRI 1998 Technical Class [START_REF]Fall Meeting. Technical Class Guidelines for Safety Assessment of Salt Caverns[END_REF]. We do not model here this damage process but deal only with an equilibrium situation. As a consequence, we have for the pressures P a the expansions

P s q s ; q f q 0 s o 2 e oq 2 s q 0 s ;q 0 f q s À q 0 s q 0 s o 2 e oq f oq s q 0 s ;q 0 f q f À q 0 f O2; P f q s ; q f q 0 f o 2 e oq 2 f q 0 s ;q 0 f q f À q 0 f q 0 f o 2 e oq f oq s q 0 s ;q 0 f q s À q 0 s À Á O2 ;
51 which we will write in the compact form

P a b A ab q b À q 0 b À Á : 52
For compressible constituents, we assume the linear relation

qb q0 b c b q b À q 0 b À Á ; 53
where the positive material constant c b represents the compressibility of the constituent b: It vanishes for a density-preserving constituent.

The constitutive part of volume exchange force

b P b m b 1 À n b 1 À m b dq b dq b À1 54
can now be linearized: we will introduce the denotation

M p M 0 b M b q b À q 0 b À Á ; 55
where the constants M 0 and M b are easily de®ned in terms of the constants A ab ; c b ; q0 b and q 0 b . Their expressions are rather long and will be omitted here.

A one-dimensional pressure-driven fluid penetration problem

In this Section, we consider a homogeneous solid matrix saturated by a ¯uid occupying the halfspace x > 0. We will assume that all considered ®elds depend only on the variable x; that the potential / vanishes, and the boundary of the solid matrix is located at x 0: Here, the saturated solid matrix is in contact with the saturating ¯uid, which is kept at a ®xed pressure: the properties of the interface between the saturated solid matrix and the ¯uid are constitutively modelled by the coef®cients a and l appearing in (46) and introduced in [START_REF] Dell'isola | A qualitative Analysis of the dynamics of sheared and pressurized layer of saturated soil[END_REF], and by the coef®cient D introduced in the following formula (56) 3 : These coef®cients, respectively, determine the part of the ¯uid pressure transmitted to solid matrix, to the saturating ¯uid and the contact double force exerted by the ¯uid on the solid matrix. The boundary conditions at x I are obtained by assuming that, at this boundary, the saturated solid matrix is in contact with an impervious solid.

The general case of compressible constituents

The boundary conditions (45) at x 0 and at x I; implied by ( 46) and (47), and valid when a 1 and l 1; are, respectively, at x 0 where for determining the integration constant the boundary conditions (56) 1 and (56) 2 were used. On the other hand, Eq. ( 15) for a f admits the following ®rst integral:

À P s pm s 1 À m s dq s dq s ! k s q s d 2 q s dx 2 À k s 2 dq s dx 2 Àm s 0p i ; P f pm f 1 À m f dq f dq f 2
oe oq f p qf 1 À q f qf dq f dq f 2 3 c ; 59
where c is an integration constant. Eliminating p from Eqs. ( 58) and (59) yields an ordinary differential equation, which posseses the normal form

k s q s d 2 q s dx 2 À k s 2 dq s dx 2 Fq s ; q f ; c p i ; 60
Fq s ; q f ; c a P a 1 À a m 2 a dq a dq a 2 3 1 À q f qf dq f dq f 2 3 À1 c À oe oq f 2 3 qf : 61
Differentiating (60) yields d dx k s q s d 2 q s dx 2 À k s 2 dq s dx 2 4 5 k s q s d dx d 2 q s dx 2 d dx Fq s ; q f ; c ; 62
The integration constants c and g 1 c; p i ; D, together with the further two integration constants of (63), are determined by using the boundedness of the solution (57) as x 3 I; and (56) 2;3 for x 0:

Incompressible constituents

To show some of the qualitative features implied by (63)±(65) we assume that dq a dq a 0; q s q0 f q f q0 s q0 f q0 s : 66

The boundary conditions at x 0 are

À A ss q s À q 0 s A sf q f À q 0 f p q s q0 s ! k s q s d 2 q s dx 2 À k s 2 dq s dx 2 À q s 0 q0 s p i ; A ff q f À q 0 f A fs q s À q 0 s p 1 À q s q0 s 1 À q s 0 q0 s p i ; k s q s dq s dx Dp i ; 67
and Eq. (58) reduces to

p À a P a k s q s d 2 q s dx 2 À k s 2 dq s dx 2 p i : 68
Replacing the expression for p in (67) 2 , we get at x 0

A ff q f À q 0 f A fs q s À q 0 s À a P a k s q s d 2 q s dx 2 À k s 2 dq s dx 2 4 5 @ A 1 À q s q0 s 0 : 69
On the other hand, Eq. (59) takes the form

A ff q 0 f q f À q 0 f A fs q 0 f q s À q 0 s À Á p q0 f c : 70
By eliminating p from (68) and (70), we obtain the ordinary differential equation for q

s À A ff q 0 f q f À q 0 f À A fs q 0 f q s À q 0 s c 4 5 q0 f ÀA ff q f À q 0 f À A ss q s À q 0 s À Á k s q s d 2 q s dx 2 À k s 2 dq s dx 2 4 5
p i ; 71 which has the normal form

k s q s d 2 q s dx 2 À k s 2 dq s dx 2 A ff 1 À q0 f q 0 f A sf 4 5 q f À q 0 f A ss À A fs q0 f q 0 f 4 5 q s À q 0 s À Á À p i cq 0 f : Fq s ; q f ; c : 72
Because of (62) 2 , it implies

k s d 2 q s dx 2 L ln q s q 0 s g 9 L q s À q 0 s q 0 s g L q s q 0 s À L À g ; 73
where

L A ss À A fs q0 f q 0 f À A ff 1 À q0 f q 0 f 2 3 A sf 4 5 q0 f q0 s : 74
It possesses the dimension of energy per unit mass, and may, according to its de®nition, be positive or negative, depending on the relative magnitude of the coef®cients A ab : This may give way to a branching solution. However, in anticipation that A ss is much larger than A ff , A sf or A fs ; L is likely positive, and the solution of the ordinary differential equation ( 73) is decaying in x: In this spirit, therefore,

q s x q 0 s L À g L C 1 e L q 0 s ks q x C 2 e À L q 0 s ks q x : 75
Now we impose that 1. At x 3 I, the boundary conditions (57) hold. In particular, (57) 3 implies that

C 1 0 A d 2 q s dx 2
xI 0 and dq s dx xI 0 ; 76 and, therefore, with simple algebra, (57) 1;2 become P a j xI 0 A q a I q 0 a a s; f A g 0 : 77

In this deduction one also must use pI p i , which is the result of global equilibrium. 2. At x 0, dq s 0 dx Dp i k s q s 0 ; 78 which implies

C 2 q 0 s C 2 À Á ÀDp i q 0 s k s L s : 79
Because of the linearization procedure, which we have used, the only meaningful solution for C 2 is given by

C 2 q 0 s À 1 2 1 2 1 À 4 Dp i q 0 s À Á 2 q 0 s k s L s 2 3 v u u t 9 À Dp i q 0 s À Á 2 q 0 s k s L s : 80
In conclusion we have q s x q 0 s À Dp i q 0 s q 0 s k s L s e À L q 0 s ks q x : 81 8 Static permeability of porous solid matrices and conclusions Formula (81) shows that the second-gradient theory, which has been introduced in this paper, features important differences when compared with the classical theory of Fillunger and Terzaghi. Indeed, when k s 0 the apparent solid density is constant and independent of the pressure p i : On the other hand, when the second-gradient strain has an in¯uence on the stress, and, hence, k s T 0; then (81) shows that the length

x 0 : q 0 s k s L r ; 82
plays the role of an attenuation length for the apparent solid mass density. We will, henceforth, call the quantity x 0 static permeability of the solid matrix as it is this characteristic length of the exponential decay (81) that measures the amount of saturating ¯uid which has penetrated into the solid porous matrix under equilibrium conditions. It is evident that the introduced static permeability is quite different from the dynamic Darcy permeability, which is de®ned in terms of relative velocities driven by pressure jumps or gradients.

The density variation Dq s : Dp i q 0 s q 0 s k s L s ; 83 which depends on the contact double-force coef®cient D; measures the intensity of the p iinduced drop of the apparent solid mass density, and can be called static permeability of the solid matrix boundary.

To conclude, we remark that:

1. The expression (58), which is valid under equilibrium conditions, shows that in the generalized model the saturation pressure can no longer be interpreted immediately as a porewater pressure; indeed, under equilibrium conditions, the ¯uid pressure in a connected region is constant but p is not. 2. Again (58) shows that the compressibility of the saturating ¯uid and solid matrix greatly in¯uence the values of the saturation pressure. 3. The expressions for static permeabilities, which are found in this Section, are valid only for the case of incompressible solid and ¯uid constituents. 4. Although found for a one-dimensional case and, therefore, apparently far from engineering applications, formula (81) gives a ®rst interesting applicable answer to the problems raised recently in [22±24].

) reduces to Àr m a p prm a b a 0 ; 24 or

 24 Àm a rp b a 0 : 25

2 and 3

 3 only by the contributions due to the additional term involving the variables f ab : Indicating this contribution by the index add , we may write which rq H a denote the time derivatives following the motion of the solid and ¯uid, respectively, viz.

ab rq a rq b f ab I 2 3 : rv b dV 2 ab q b I rq a 2 3 ;

 323 r a Á v a : 32 Using (32) in (31) and employing the rule Àrq b Á rq H a rq a rq b f ab I : rv b q b I rq a . . . rrv b ; 33 obtained by recalling the balances of mass for the constituents, where no summation over repeated indices is performed, we may derive the identity 36 the right-hand side of (34) can be regarded as the sum of terms possessing the structure B A b : rv b C b . . . rrv b dV ; 37

  , since rm f Àrm s straightforward arithmetics shows that qa rm a rq a 1 À m a

		m
					dq a dq a	; q a	ow oq a	m a qa	oe oq a	p 1 À m a	dq a dq a	!	;	18
	p a Àq a	w q	q a	ow oq a	P a pm a 1 À m a	dq a dq a	;	19

a M rm a ; M : p b

  As an illustration, we choose oe=of ab 0, except for

	oe of ss	:	k s 2	:							43
	Then the ®eld equations become
	Àrp s m s b s r Á Àrp f m f b f 0 ;	k s q s Rq s	k s 2	j	rq s	j 2	I À k s rq s rq s	!	0;	44
	and the boundary conditions take the form
	Àp s k s q s Rq s Àp f n t f 0; k s q s k s 2 j oq s rq s on	j 2	n À k s rq s rq s d 0 ;	Án t s 0;	45
	in which
	d

2 : f 0; d s dI :

Following the results found in

[19±21]

for the constituent tractions, denoting p i the incumbent pressure exerted on the saturated solid matrix at its boundaries, we choose t s am l s p i n; t f 1 À am l s À Á p i n ; 46

In the formulation of[START_REF] Svendsen | On the thermodynamics of a mixture of isotropic materials with constraints[END_REF], the ¯uid stresses are absent in the reduced formulation, corresponding to this one due to a slip in calculations.

Fq s ; q f q s ; cdf s 64 and the function q f Á is implicitly de®ned by the saturation constraint q s qf q f q f qs q s qf q f q s q s :


Tao Wang
email: tao.wang@uclouvain.be

Fellow, IEEE Luc Vandendorpe
email: luc.vandendorpe@uclouvain.be

Iterative Power and Subcarrier Allocation for Maximizing WSMR in Cellular OFDMA Systems

Keywords: Orthogonal frequency division multiplexing, resource allocation, cochannel interference mitigation, spectrum management

We consider the resource allocation (RA) problem of maximizing the weighted sum of the minimal user rates (WSMR) of coordinated cells subject to a total power constraint at each base station (BS) in the downlink of a cellular orthogonal frequency division multi-access (OFDMA) system. In particular, the solution of this problem corresponds to a RA that guarantees similar rates to users in each cell. We propose a coordinate ascent (CA) based algorithm which alternatively updates the power allocation by solving a successive set of convex optimization problems with a duality-based numerical algorithm, and the subcarrier allocation by solving a mixed integer linear program for each cell. The effectiveness of the algorithm is illustrated by numerical experiments.

I. INTRODUCTION

Next-generation cellular systems envision full spectrum reuse to meet higher data-rate requirements. To this end, joint resource allocation (RA) among coordinated base stations (BSs) is a promising way to mitigate cochannel interference [START_REF] Gesbert | Adaptation, coordination, and distributed resource allocation in interference-limited wireless networks[END_REF]. In this paper, we consider the joint RA based BS coordination in the downlink of a cellular orthogonal frequency division multi-access (OFDMA) system.

Lately, a few margin-adaptive joint RA algorithms have been proposed for the downlink of cellular OFDMA systems [START_REF] Pietrzyk | Radio resource allocation for cellular networks based on OFDMA with QoS guarantees[END_REF], [START_REF] Abrardo | Centralized radio resource allocation for OFDMA cellular systems[END_REF]. In case that the joint RA is optimized for each subcarrier individually, the rate-adaptive algorithms in [START_REF] Kiani | Maximizing multicell capacity using distributed power allocation and scheduling[END_REF]- [START_REF] Gjendemsjφ | Binary power control for sum rate maximization over multiple interfering links[END_REF] can be adopted. However, a joint RA over all subcarriers is preferred in order to better exploit the frequency and multiuser diversity inherent in OFDMA systems. To this end, two rateadaptive joint RA algorithms, which aim to maximize the weighted sum-rate of all users, have been proposed in [START_REF] Venturino | An improved iterative waterfilling algorithm for multi-cell interference mitigation in downlink OFDMA networks[END_REF] and [START_REF]A successive convex approximation algorithm for weighted sumrate maximization in downlink OFDMA networks[END_REF]. In particular, the weight of each user represents the priority of this user's rate. However, it is not clear how to determine the weight of each user so as to fulfill a prescribed fairness criterion among users.

Motivated by the applications requiring the maximum fairness among users in each cell, we consider a novel RA problem which aims to maximize the weighted sum of the minimal user rates (WSMR) of coordinated cells subject to a total power constraint at each BS. Most interestingly, the optimum RA guarantees the maximum fairness in each cell, i.e., the users in each cell have similar rates. In a previous contribution, we proposed a coordinate ascent (CA) based iterative algorithm which optimizes the subcarrier allocation and the power allocation alternatively, so that the WSMR keeps increasing until convergence [START_REF] Wang | Resource allocation for maximizing weighted sum min-rate in downlink cellular ofdma systems[END_REF]. At each iteration, the subcarrier allocation is updated by solving a mixed integer linear program (MILP) for each cell, then the power allocation is updated by solving a successive set of convex-optimization problems. To solve these problems, we proposed a heuristic iterative method based on Karush-Kuhn-Tucker (KKT) conditions in [START_REF] Wang | Resource allocation for maximizing weighted sum min-rate in downlink cellular ofdma systems[END_REF]. Though simple and effective in case of convergence, this KKT-based method is not guaranteed to converge.

In this paper, we make improvement on the power allocation algorithm proposed in [START_REF] Wang | Resource allocation for maximizing weighted sum min-rate in downlink cellular ofdma systems[END_REF]. Specifically, our contributions lie in the following points:

• We add extra constraints to each convex-optimization problem involved, so as to ensure that its solution set is nonempty and the successive optimization process of the power allocation algorithm can proceed successfully. • We propose a duality-based method to solve each convexoptimization problem involved. Specifically, this method particularizes the subgradient method to find the dual optimum, and its convergence is guaranteed provided that the stepsize is sufficiently small. The rest of this paper is organized as follows. In the next Section, we present the system model as well as a concave lower bound of user rate. In Section III, we formulate the RA problem and describe the general framework of the proposed algorithm. In Section IV and V, we develop the methods for initialization and subcarrier allocation, as well as the power allocation algorithm, respectively. In Section VI, the effectiveness of the proposed algorithm is illustrated by numerical experiments. Finally, some conclusions wrap up this paper in Section VII.

II. SYSTEM MODEL AND A LOWER BOUND OF USER RATE

A. System model

We consider the downlink of a cellular OFDMA system with K BSs coordinated by a central controller running the proposed RA algorithm. Using OFDMA with N subcarriers, BS k serves a group of users belonging to the user set U k for cell k. We denote the channel gain of subcarrier n from BS l to a user u in cell k by g l u,k,n , and we assume g l u,k,n > 0, ∀u ∈ U k , ∀k, l, ∀n. In addition, σ 2 u,k,n represents the noiseplus-uncoordinated-interference power for subcarrier n at user u in cell k.

We make the following assumptions for the considered system. First, the number of users in each cell is smaller than N , and each subcarrier is allocated to only one user exclusively in each cell. Second, we assume the controller knows {g l u,k,n , ∀u, ∀k, l, ∀n} and {σ 2 u,k,n , ∀u, ∀k, ∀n} precisely. Third, we assume the channel between each BS-user pair remains unchanged for a sufficiently long duration, so that the RA is carried out accordingly for that duration.

We denote the total power available to BS k by P k,max . P k,n represents the proportion of P k,max allocated to subcarrier n by BS k. We stack {P k,n } K k=1 into a K × 1 vector p n = [P 1,n , • • • , P K,n ] T , and then stack {p n } N n=1 into a K × N matrix P satisfying [P] k,n = [p n ] k,1 = P k,n . Note that every possible P belongs to the set P f = {P|∀k, ∀n, P k,n ∈ [0, 1]}. We define a binary variable A u,k,n , which indicates that subcarrier n is allocated to user u in cell

k if A u,k,n = 1. We stack {A u,k,n } N n=1 into the vector A u,k = [A u,k,1 , • • • , A u,k
,N ] T and then stack A u,k 's of all users in cell k to a matrix A k column by column. Finally, we define a matrix A = [A 1 , • • • , A K ] which indicates how subcarriers are allocated to all users.

The rate of user u in cell k, as a function of A u,k and P, is expressed by

R u,k (A u,k , P) = N n=1 A u,k,n • R u,k,n (p n ) (1)
where R u,k,n (p n ) stands for the rate of subcarrier n when allocated to user u in cell k. Note that R u,k,n (p n ), as a function of p n , is calculated by

R u,k,n (p n ) = ln 1 + γ u,k,n (p n ) Γ (2)
in the unit of nats per OFDM symbol, where γ u,k,n (p n ) is the signal-to-noise-plus-interference ratio, and Γ represents the signal-to-noise-ratio gap between the adopted modulation and coding scheme and the one achieving channel capacity. Particularly, γ u,k,n (p n ) is evaluated by

γ u,k,n (p n ) = P k,max P k,n g k u,k,n σ 2 u,k,n + K l=1\k P l,max P l,n g l u,k,n = P k,n I u,k,n (p n ) (3)
where

I u,k,n (p n ) = x u,k,n + K l=1\k P l,n y l u,k,n , x u,k,n = σ 2 u,k,n P k,max g k u,k,n
, and y l u,k,n = P l,max g l u,k,n P k,max g k u,k,n .

B. A concave lower bound of user rate

We first define

R u,k,n (e q , p * n ) = α * u,k,n • ln(γ u,k,n (e q )) + β * u,k,n (4)
α * u,k,n = γ u,k,n (p * n ) Γ + γ u,k,n (p * n ) (5)
β * u,k,n = R u,k,n (p * n ) -α * u,k,n • ln(γ u,k,n (p * n )) (6)
where q is a K ×1 vector, P * is a given power allocation with its n-th column denoted by p * n . In particular, e X is a matrix mapped from X entrywise through the function y = e x . Note that P * should belong to the set P + = {P|∀k, ∀n, P k,n ∈ (0, 1]}.

We now define a function

R u,k (A u,k , e Q , P * ) = N n=1 A u,k,n R u,k,n (e qn , p * n ) (7)
where Q is a K × N matrix with its n-th column denoted by q n . We have shown in [START_REF] Wang | Resource allocation for maximizing weighted sum min-rate in downlink cellular ofdma systems[END_REF]

that R u,k (A u,k , e Q , P * ) is a concave function of Q as well as a lower bound of R u,k (A u,k , e Q ).
Besides, this bound is tight when e Q = P * .

III. PROBLEM FORMULATION AND OVERALL ALGORITHM

We consider the RA problem that maximizes the WSMR subject to a total power constraint for each BS, namely max

A,P f (A, P) = K k=1 w k • min u∈U k R u,k (A u,k , P) s.t. P ∈ P f , N n=1 P k,n ≤ 1 ∀ k, u∈U k A u,k,n ≤ 1 ∀ k, ∀ n, A u,k,n ∈ {0, 1} ∀u ∈ U k , ∀ k, ∀ n. (8)
where w k > 0 represents the weight assigned to cell k's minimal user rate. In particular, increasing w k leads to a higher priority assigned to the users in cell k.

Algorithm 1 Overall RA optimization algorithm

Initialize: m = 0, compute P 0 and A 0 using the initialization algorithm in Subsection IV-A;

repeat Use the power allocation algorithm proposed in Subsection V-A to find a suboptimal solution for (8) with A = A m , and assign it to P m+1 ; Use the subcarrier allocation algorithm proposed in Subsection IV-B to solve (8) with P = P m+1 for the optimal A, and assign it to

A m+1 ; m = m + 1;
until ||vec(P m -P m-1 )|| ≤ Ω 1 or m = M Output: P m and A m are adopted for the considered system.

We propose a CA based iterative algorithm which optimizes the subcarrier allocation and the power allocation alternatively. This algorithm is depicted in Algorithm 1, where the superscript m indicates the associated variable is produced after the m-th iteration. Note that the algorithm yields non-decreasing WSMRs as the iterations continue. Since each BS's total power is upper bounded, the optimal WSMR is upper bounded as well, and therefore the WSMR yielded will not increase infinitely as the iterations proceed. We stop the iteration when ||vec(P m -P m-1 )|| is smaller than a prescribe small positive value Ω 1 , or when the number of iterations exceeds a prescribed value M . Here || • || and vec(•) represent the Euclidean norm and the vectorization operator, respectively. Finally, P m and A m produced by the last iteration are adopted as a suboptimal solution to [START_REF]A successive convex approximation algorithm for weighted sumrate maximization in downlink OFDMA networks[END_REF].

IV. INITIALIZATION AND SUBCARRIER ALLOCATION ALGORITHM

A. Initialization methods

We use the uniform power allocation (UPA) for initialization, i.e., [P 0 ] k,n = 1 N , ∀n, ∀k. We consider two different methods for computing A 0 . One is to let users in each cell to choose in turn the subcarrier that is unassigned and corresponds to the highest SNIR. In this way, each user will be allocated almost the same number of subcarriers, and we refer to this method as the even subcarrier allocation (ESA). The other one is to compute A 0 as the best subcarrier allocation (BSA) that maximizes the WSMR when P = P 0 . Note that this subcarrier allocation can be obtained by solving a MILP problem for each cell, which will be elaborated in Section IV-B.

B. Subcarrier allocation optimization

The problem of finding the optimal subcarrier allocation A m+1 when P = P m+1 can be decomposed into K subproblems, the k-th of which is to find the optimal subcarrier allocation for cell k: max

A k min u∈U k R u,k (A u,k , P m+1 ) s.t. u∈U k A u,k,n ≤ 1, ∀ n, A u,k,n ∈ {0, 1}, ∀ u ∈ U k , ∀ n. (9)
Let's denote the optimal A k for (9) by A m+1 k . Then the optimal subcarrier allocation is constructed by

A m+1 = [A m+1 1 , • • • , A m+1 K ].
In fact, (9) can be reformulated into an equivalent form:

max A k ,S k S k s.t. N n=1 A u,k,n • R u,k,n (p m+1 n ) ≥ S k , ∀u ∈ U k , u∈U k A u,k,n ≤ 1, ∀ n, A u,k,n ∈ {0, 1}, ∀ u ∈ U k , ∀ n. (10
) where p m+1 n is the n-th column of P m+1 , and S k is the extra optimization variable introduced to guarantee the equivalence between ( 9) and [START_REF]Iterative resource allocation for maximizing weighted sum minrate in downlink cellular ofdma systems[END_REF]. Note that (10) is a MILP problem involving both integer and continuous variables. In practice, we can solve it by using a commercial optimization software such as CPLEX.

V. POWER ALLOCATION ALGORITHM

A. Successive convex optimization algorithm for power allocation

When A = A m , we use R m u,k (P) and R m u,k (e Q , P * ) to represent respectively the rate of user u in cell k and R u,k (A m u,k , e Q , P * ) for notation simplicity. Here, A m u,k represents the A u,k contained in A m . In the following, we denote the user in cell k scheduled for subcarrier n by u m k,n , and we define the sets P = {P|∀k, ∀n, P k,n ∈ [e , 1]} and Q = {Q|∀k, ∀n, q k,n ∈ [, 0]}, where q k,n = [Q] k,n and is a prescribed negative value.

To develop a method for optimizing the power allocation, let's consider a parameterized convex optimization problem

min Q,{C k } K k=1 - N k=1 w k • C k s.t. Q ∈ Q , C k -R m u,k (e Q , P * ) ≤ 0, ∀ u ∈ U k , ∀ k, N n=1 e q k,n -1 ≤ 0, ∀ k (11
) According to Theorem 2 in [START_REF]Iterative resource allocation for maximizing weighted sum minrate in downlink cellular ofdma systems[END_REF], the solution set to (11) is nonempty. Let's denote a solution to [START_REF] Bertsekas | Nonlinear programming[END_REF] by Q o and suppose P * is feasible for (8) and P * ∈ P , P o = e Qo is feasible for (8) and f (A m , P o )≥f (A m , P * ) holds. Capitalizing on these facts, we propose an iterative power-allocation optimization algorithm depicted in Algorithm 2, which introduces an inner loop nested in the overall algorithm, to find a power allocation no worse than P m when A = A m . Note that T and Ω 2 are a prescribed integer and a small positive value, respectively.

Algorithm 2 Power allocation optimization algorithm

Initialize: t = 0, P m,0 = P m ; repeat Solve [START_REF] Bertsekas | Nonlinear programming[END_REF] with P * = P m,t for Q o using Algorithm 3

proposed in Section V-B; P m,t+1 = e Qo ; t = t + 1;

until ||vec(Q o -Q * )|| ≤ Ω 2 or t = T Output: P m+1 = P m,T .
We can see that Algorithm 2 produces a set of power allocations yielding nondecreasing WSMRs. When the WSMR stops increasing, f (A m , P o ) = f (A m , P * ) holds. In this case, P * might be a local minimum of (8) with A = A m and P ∈ P , since the KKT conditions are satisfied by P * according to Theorem 3 in [START_REF]Iterative resource allocation for maximizing weighted sum minrate in downlink cellular ofdma systems[END_REF].

In fact, Algorithm 2 iteratively looks for a suboptimal power allocation in P , while a possible P actually lies in a bigger set P f ⊃ P . However, we can not replace P by P f , because we can not construct (11) for a P * ∈ P f with zero entries as explained in Section II-B. We can not replace P by P + either, since when Q + is used the set of feasible Q's is unbounded, and the existence of an optimal solution to (11) is not guaranteed [START_REF] Bertsekas | Nonlinear programming[END_REF]. On one hand, e should not be too small to be rounded to zero by the central controller. On the other hand, should be much smaller than 0, so that P is a good approximation to P + and Algorithm 1 is able to start with P 0 ∈ P . B. A duality-based algorithm to solve [START_REF] Bertsekas | Nonlinear programming[END_REF] We denote the Lagrange multipliers related to the rate constraint of user u in cell k and the power constraint of BS k by λ u,k and µ k , respectively. We stack the λ u,k 's related to all users in cell k into a row vector λ k , and stack the µ k 's related to all BSs into a row vector µ = [µ 1 , • • • , µ K ]. We then stack all λ k 's and µ into a row vector ρ

= [λ 1 , • • • , λ K , µ].
We define the Lagrangian of [START_REF] Bertsekas | Nonlinear programming[END_REF] as [START_REF] Bertsekas | Nonlinear programming[END_REF]. Since ∀k C k can take any real value for the minimization operator in (12), u∈U k λ u,k = w k , ∀k must be fulfilled so as to guarantee that D(ρ) > -∞. As a result, D(ρ) can be reduced to

D(ρ) = min Q∈Q ,C h(Q, ρ) + K k=1 ( u∈U k λ u,k -w k )C k (12) where C = [C 1 , • • • , C K ] and h(Q, ρ) = K k=1 µ k ( N n=1 e q k,n -1) - K k=1 u∈U k λ u,k R m u,k (e Q , P * ) (13) Note that the domain of D(ρ) is D ρ = {ρ|D(ρ) > -∞; ∀k, µ k ≥ 0; ∀u ∈ U k , ∀k, λ u,k ≥ 0}
D(ρ) = min Q∈Q h(Q, ρ) (14)
and

D ρ is reduced to {ρ|∀k, λ k ∈ D k ; µ ∈ D µ }, where D µ = {µ|∀k, µ k ≥ 0}, and D k = {λ k | u∈U k λ u,k = w k ; ∀u ∈ U k , λ u,k ≥ 0}.
Apparently, D µ and D k are both convex, and thus D ρ is also a convex set.

Let's define Q ρ = arg min Q∈Q h(Q, ρ). Specifically, the duality-based algorithm to solve [START_REF] Bertsekas | Nonlinear programming[END_REF] is to first find the optimal dual variable ρ = arg max ρ∈Dρ D(ρ), and then compute Q ρ as the optimal solution to [START_REF] Bertsekas | Nonlinear programming[END_REF]. We propose Algorithm 3 which particularizes the subgradient method in [START_REF] Bertsekas | Nonlinear programming[END_REF] to find ρ . Here, ρ i represents the dual variable produced after the i-th iteration, η i denotes the subgradient of D(ρ) when ρ = ρ i , and δ and ∆ 1 are both prescribed small positive values.

Algorithm 3

The duality-based algorithm to solve [START_REF] Bertsekas | Nonlinear programming[END_REF]. [START_REF] Bertsekas | Nonlinear programming[END_REF].

Initialize: i = 0; compute ρ 0 in which ∀u ∈ U k , ∀k, λ u,k 0 = w k |U k | and ∀k, µ k 0 = 1. repeat Compute Q ρ i using the algorithm proposed in Section V-C; Compute η i and ρ i+1 = arg min ρ∈Dρ ||ρ i + δ • η i -ρ|| 2 . i = i + 1; until ||ρ i -ρ i-1 || ≤ ∆ 1 Output: Q ρ i is the optimal solution to
Note that |U k | is the number of users in cell k, η i has the same size as ρ, and the entries in η i associated with λ u,k and µ k are equal to -R m u,k (e Q ρ i , P * ) and

N n=1 e [Q ρ i ] k,n -1, respectively.
Using the associated KKT conditions, we can compute ρ i+1 by

λ i+1 u,k = [λ i u,k -δ • R m u,k (e Q ρ i , P * ) + ξ i k ] + , ∀u ∈ U k , ∀k (15)
µ i+1 k = [µ i k + δ • ( N n=1 e [Q ρ i ] k,n -1)] + , ∀k (16)
where the superscript i indicates the associated dual variable is an entry in ρ i , [x] + = max{0, x}, and ξ i k is a real value that satisfies u∈U k λ i+1 u,k = w k . We find the value of ξ i k by the bisection method, since u∈U k λ i+1 u,k increases monotonically with respect to ξ i k . Provided that δ is sufficiently small, the convergence of Algorithm 3 is guaranteed [START_REF] Bertsekas | Nonlinear programming[END_REF]. We terminate the iterations when ||ρ iρ i-1 || is smaller than a prescribed small positive value ∆ 1 , then Q ρ i produced in the last iteration is taken as the optimal solution to [START_REF] Bertsekas | Nonlinear programming[END_REF].

C. An algorithm to compute Q ρ i

Thanks to the convexity of h(Q, ρ i ), we resort to a gradientprojection based iterative algorithm to find Q ρ i . This method begins with initializing Q by

Q * if i = 0, otherwise by Q ρ i-1 . Then, Q is iteratively updated by Q = [Q-τ •B] Q ,
where B is a matrix containing the gradients of h(Q, ρ i ) with respect to every entry of Q, [•] Q is the operator of projection into Q , and τ is a prescribed small positive value that guarantees convergence. We terminate the iterations when ||vec(

[Q -τ • B] Q -Q)|| is smaller than a prescribed small value ∆ 2 . Then Q produced by the last iteration is taken as Q ρ i .
In fact, every entry of Q is updated by

∀k, ∀n, q k,n =    if q k,n -τ • [B] k,n ≤ 0 if q k,n -τ • [B] k,n ≥ 0 q k,n -τ • [B] k,n otherwise (17) where [B] k,n = ∂h(Q,ρ i ) ∂q k,n is derived as [B] k,n = -λ i u m k,n ,k α * u m k,n ,k,n + e q k,n µ i k + (18) l =k λ i u m l,n ,l α * u m l,n ,l,n y k u m l,n ,l,n I u m l,n ,l,n (e qn ) . (19)

VI. NUMERICAL EXPERIMENTS

For illustration purposes, we consider the downlink of a cellular OFDMA system with K = 3 coordinated cells and N = 8 subcarriers, as shown in Figure 1. Note that

U 1 = {u 1 , u 2 }, U 2 = {u 3 , u 4 }, and U 3 = {u 5 , u 6 }. We set w 1 = w 2 = w 3 = 1, Γ = 1, σ 2 u,k,n = 1, ∀ u ∈ U k , ∀ k, ∀ n, P k,max = P M , ∀ k, = -10 which corresponds to e = 4.54 × 10 -5 , δ = τ = 0.01, Ω 1 = Ω 2 = ∆ 1 = ∆ 2 = 0.001, M =
4, and T = 5 in numerical experiments. We conducted the following numerical experiments with Matlab v7.1 on a laptop equipped with an AMD Turion CPU of speed 2.2 GHz and a memory of 2 GBytes. In particular, we used the TOMLAB/CPLEX toolbox v12.1 to solve [START_REF]Iterative resource allocation for maximizing weighted sum minrate in downlink cellular ofdma systems[END_REF].

We generate the channel between each BS-user pair according to the following assumptions. First, it is modeled by a 6-tap delay line, and the average received power decays exponentially with distance. In particular, the propagation exponent is equal to 3. Second, we assume the amplitude of the i-th tap is a circularly symmetric complex Gaussian random variable with zero mean and variance as σ 2 i . Finally, we assume

σ 2 i+1 σ 2 i
= e 3 , i.e., the tap power decays exponentially with the coefficient equal to 3.

We have tested the proposed algorithms on a random realization of the BS-user channel gains shown in Figure 2 when P M = 90 dBW. To facilitate the following explanation, we define a 3 × 8 user-scheduling matrix U m which contains u m k,n at the k-th row and the n-column. Apparently, there exists an one-to-one mapping between U m and A m . When the UPA as well as the BSA and the ESA are used for initialization, the corresponding U 0 's are evaluated respectively as After Algorithm 1 terminates, the computed WSMR's with respect to mT + t, as well as user rates at initialization and after optimization are shown in Figure 3. In particular, each line segment with the same mark shows the WSMR corresponding to iterations of Algorithm 2 with the subcarrier allocation fixed. After the algorithm terminates, P 4 and U 4 corresponding to the BSA and ESA initialization are feasible for [START_REF]A successive convex approximation algorithm for weighted sumrate maximization in downlink OFDMA networks[END_REF], and they are computed respectively as

U 0 =   2 1 2 2 2
P 4 =            
0.0066 0.0008 0.3626 0.1075 0.0063 0.0553 0.0005 0.4094 0.0053 0.0152 0.0380 0.0905 0.0018 0.0003 0.4443 0.1497 0.0701 0.0003 0.6557 0.0115 e 0.0039 0.4678 0.0004

            , (22)
and

P 4 =             0
            , ( 24
) U 4 =   2
We can see that both initialization methods lead to the similar WSMR after optimization, even though the BSA initialization results in a higher WSMR at the beginning. Besides, the power and subcarrier allocation after optimization results in a WSMR about 10 nats/OFDM-symbol higher than the initial UPA and BSA/ESA. Despite the initialization method, the following points can be observed. First, the WSMR is increased significantly by the power allocation optimization at the beginning, while the subcarrier allocation optimization leads to a slight increase of the WSMR. Second, the power allocation after optimization shows that each BS should allocate the dominant part of its power to only a few subcarriers. Once a BS allocates a significant part of its power to a subcarrier, the other BSs should allocate only a small part of its power to this subcarrier, so as to reduce CCI. Third, the minimal user rate in each cell is enhanced, and the users in each cell have similar rates after optimization. These observations show that the power and subcarrier allocation yielded by the proposed algorithm leads to a significantly better WSMR performance than the UPA and the BSA/ESA for the considered random realization of channel gains.

We have also tested the proposed algorithm on 100 random realizations of the BS-user channel gains, when P M varies from 60 to 100 dBW. The average value of the WSMR at initialization and after optimization are shown in Figure [START_REF] Kiani | Maximizing multicell capacity using distributed power allocation and scheduling[END_REF]. We can see that two initialization methods result in the similar average WSMR after optimization, although the BSA initialization leads to a higher average initial WSMR. Regardless of the initialization method, the average WSMR after optimization is slightly higher than the average initial WSMR when P M is below 80 dB. This is because the CCI is relatively low due to the relatively low total power of each BS, and therefore the UPA is nearly optimal as explained in Section IV-A. However, the average initial WSMR increases slightly, while the average WSMR after optimization increases significantly, when P M increases from 80 dB to 100 dB. Besides, the average WSMR after optimization is significantly higher than the average initial WSMR when P M is greater than 80 dB. These observations show that the proposed algorithm can effectively produce the power and subcarrier allocations that lead to the significantly better WSMR performance in average than the UPA and the BSA/ESA, especially when the total power of each BS is high.

VII. CONCLUSION

We have considered the downlink of a cellular OFDMA system, in which multiple BSs are coordinated by a centralized resource allocation algorithm. We have proposed an iterative algorithm to maximize the WSMR subject to a total power constraint at each BS, in terms of jointly optimizing coordinated BSs' subcarrier and power allocation. The effectiveness of the algorithm has been illustrated by numerical experiments.

6 Fig. 1 .

 61 Fig. 1. An exemplary system setup.

Fig. 2 .

 2 Fig.2. A random realization of channel gains, where the blue, green, and red bars represent the channel gains from BS 1, 2, and 3, respectively.

Fig. 3 .

 3 Fig. 3. The computed results for the random realization of channel gains.

Fig. 4 .

 4 Fig. 4. Average WSMR over 100 random realizations of channel gains.

ACKNOWLEDGMENT

The authors would like to thank the Walloon Region for funding the project MIMOCOM, the ARC SCOOP and the FP7 Network NEWCOM++, as well as Tomlab Optimization Inc. for offering a demonstration license.