J S Perkell

D H Klatt

Ward L Evel

in INVARIANCE AND VARIABILITY lN SPEECH PROCESSES edited

As is wcll known, the word rcprescllts thc llIaill poillt of illtcrnctioll bctwccn ail sources of knowlcdge shnrcd by both thc spcnker alld thc listcllcr (Klntt, ]977). Such SOllrces of kllowledge inclllde thc lexicoll which contai liS thc list or nvnilable possiblc words ill n uLllgllnge, the sYlltnx, thc sclnnlltics, nnd thc prnglllntics. Those factors constraillS scqllcnccs or words tü rorlll (ill tllc idenl CiLSC) grmlllllntically, scllmllticnJly, nlld prnglllnticnlly correct sè'lltcnccs. The ultilllnte pllrpose of specch COllllllllllicatioll is to collvey mcnnillg from spcakcr to listcncr. As thc basic imormatioll-cl\rrying 1I1lit, (~ohclI POilltS ollt, the word rnthcr thall thc phonellle might be collsidcrcd thc buildillg block in spcech rese..rch. As for whnt concerns illvnrinnce nt the word levcl, Cohen considers the listcner's point of view nlld concludes thl\t, "the only invariance we might clnim is the OllC bnscd on word types, to be seen as moulds in the shape of gestalts, stored in our mentnllexicon, detcrlllined by the phonotactic constrains of the langunge." Cohen emphasizes the raIe or other types of imormation than acoustic îmormation availnble to the listencr to decode the successive words in the message, and cites the LAFS llIodcl (Klatt, 1979a) and the Logogen model (Morton, 1964) as possible models for speech perception. From Cohen's mentalist point of view (see Section 23a), the problems of segmentation into phoneme-size units and invariance of phoneme and feature become srnall relative to the problern of scgrnentation into word-size units and their invariance. 1 agree with the views expresscd in Cohen's paper, rnost of which serve to point out the importance of the word as an unit in speech. However, 1 would like to try to fill in Borne of the dctails of these views.

WaRD DEFINITION

Thc paper does not define in enough detail what a "word" is. The definition of the word, either as a unit of meaning or as an acoustic unit, is Dot clcar cut. First, it may refer to the graphic form of the word, and blank spaces are considcred boundary markeTs between successive words. Sccond, thc word may also be defincd as the basic unit of mcaning, "the morphcme." A "graphic" word, however, may be composed by scveral morphemcs, and morphemes (and Dot thc graphic word) are considered by linguists as the basic unit of me..ning. Morcover some words such as the grammatical words (articles or auxiliarics) do not have a mcaning by themsclvcs. Third, the separation of the "graphic words" into two word typcs and the notion of thc prosodic word have bcen introduced by refcrence to the acoustic level for languages like French and English.

The words of the first type (type A), corresponding typically to function (grammatical) words, arc generally vcry short words in terrns of numbcr of syllables, have a rcduced duration, an avcrl\gc lowcr fundamental frequcncy and arc not prcciscly articulatcdj the words of the second type (typc B), typically lexical words, correspond to a local pcnk of promincnce on one syllnble (markcd by durational and fUlldillllcntm frcqllcllcy colltrastS), with a IclIgthcllillg of thc word lill/LI syllable, alld Il strellgthcllillg of thc worrl illitial phollcme. Ucpellrlillg 011 thc contcxt, Il Icxicm ward m/LY howcvcr behavc as ."\ word of typc A /\lld fllllCtioll woro as /\ worrl of I.ype l~. Words of typc A alld of typc il nl/\Y bc rcgrollpcrl into wh/\t is c/\llcd a sillgle prosorlic word (relatcrl to th(: notiOII of "rhythmic IInit," "syntagnl," "strcss group," or "hat-p."\ttCrll") prim/\rily ch."\r/\ct.(!rizerll>y two FI) movcmcllts ill opposit.e dircctioll. SlIch regrOIlpillg obscllrcs thc /\Collstic irlclltity of thc componcllt worrls. Whcre the prosodic word st./\rts /\110 wherc it cllds is rlifficult to dctcrmine, silice thcrc is no solid tllcorctical basis for sllch a dcm/\rcation.

For cx/\mplc, the two Icxic."\1 worrls "petit g/\rçon" /\re gcllcr/\lly rcgroupcrl illtO /\ sillglc four-syll."\blc prosodic word. The t-wo words 111/\Y tclld to scp1\r/\tc howevcr illtO two dilfercllt aconstic units in certaill circnmst1\lIccS, which depcnd on the speaker, the style, the r/\tc of speech, thc ICllgth 01: tllC words (shorter worrls may be regrouperl nlore oftcn), /\nd the frequency of occurrence of thc words in the Icxicon and in the discourse. Also Il long gr/\phic, lexic1\1 worrl bas a tendency to "sprc."\d" into two prosorlic units depcnding 011 iLs morphemic composition (ilS marked in French by a Icngthcning of the fin/\l syllablc of /\n intern/\l morphcme). Intcgration or complcte scparation of two "graphic" worrls into a single acoustic unit, or separ/\tion of a graphic word into two "morphemic" units may Dot be a binary decision, but may vary along a continuum. ln other words, the rclationship anlong graphic worrls, units of meaning, and prosodic words is Dot ncccssarily a one-to-one corrcspondence.

REPRE6ENTATION OF THE WaRD

Cohen does Dot express clearly what a gestalt at the lcvel of the ward may look like, since by definition, a gestalt must be invariant.

A gestalt point-of-view is Dot inconlpatible with the vie..v of the ward as composed of smaller constitucnts (phonemes, for example) at least at Borne abstract level: Such a theory may be equally applied to the syllablc, the demisyllable, the diphonc, the phone, and the fc."\ture. Thc detniled acoustic shape of a ward is known to be iIIflucnccd by a large number of factors. By painting out the importance of other sources of knowledge (such as syntax and scmantics and phenomena S11Ch as phoneme rcstoration) Cohen suggests a gestalt of high colllplexity which questions the adcquacy of the representation of the ",ord as a string of phoncmes.

One principlc of phonetic transcription has becn to use a set of symbols in the most economical and the most efficient way to represent the various uttcrances of a language. Linguists postuJate two levels, or a continuullI of Jevels, betwcen an abstract representation-a phonemic lcvel-and surface rcalizations-the phonetic levcl. A single phonemic sYllIbol, such I\S the phoneme [RI or the phoneme [1] could be associated with spectrographically very difTcreut looking sounds, dcpcndillg on the position of the phoneme in the syllable (e.g. iIIitinl, final, and in clustcr) and iII the word. It is Dot clear, for examplc, if aIl [R] and [1] allophoncs. shnre even pcrccptual equivalcnce or invnrinnt propcrtics at the acoustic level. Thcrc mny be Borne invariallt ClIC across su ch positiollS, but the lack of iIIvariallcc acros" clilT"r"nt Strikillg 8imilarities related to thc word 1\.'1 ail acoustic IImt havc boe/l obscrvcd for a 1I1II11l>c~r of related alld nnrclatcd lallgmlgcs su ch as l';lIglish, Swcdish, Frcllch, alld Japanesc: word-final syllablc lclIgthemng, word-init.ial collsollant Icngt.hcllillg, or worcl-initial allophones. Similarit.ics also cxtcnd to thc way wit.11 which the words arc rcgrouped into largcr units (Fo ri8CS associatcd with illiti.'1tion, Fo folll alld lellgt.hemllg with termination resctting of thc base-linc as boundary markerj see Vai8sière, 19S3a, for references). Such rcgularities mOlY come from similar ways of representing words in the mental lexicon, and from gcneral processcs in composing sentenccs that are independent of the languagc spokcn. Howcver, there is a lack of knowlcdge of a large number of languages for definitive conclusions.

The segmentaI fe.'1tures have to be interprctcd depending on their position in the sentence. The word as an unit imposes strong constrai/lts on the ûcoustic st.ructuring of the utterûnce. Analysis of the vclum behavior in selltcnces shows the velum to be typically higher in word illitial position and in prestressed position than in other positions for both the na8al consonants (which require an open vclopharyngeal port) and tHe oral consonants (requiring a close velopharyngeal port) (Vaissière, 19S3bj. A higher position of thc velnm corrcsponds to a greater tensing of the levator parntini. It may bc hypothesized that the difference in velum height between position al allophones may be due to the superposition of a common suprasegmcntal featurc, let say [+strong], corrcsponding to a greater tensing at word onsets, rather than to fluctnations associated with ench phoneme. With this hypothesis, the aspiration of the so-called tense stops consonants in English [p, t, k] should be considercd as being due to the samc factor as to the partial devoicing of word iuitiallax stops [b, d, g], the glottalization of word-initial vowels, thc higher position of the vemm in word-initial position, or the fact that vowels in scntence context (at least for French) secm to be, cetcris paribus, more precisely uttered if the strncturing of the sentence requires them to be uttered in the uppcr Fu register of the speaker.

A8piration, dcvoicing, glottalization, highcr position of the velum or more precise articulation have a common charactcristic: a greûter tensing of at least one of the articulators, the vocal folds, the velum (levator palitini), or the tongue. Such extra tensing mûY be contradictory or DOt to the articulatory rcqnircments of the underlying scgmental feûtnre(s). Less intra-and intcrspeakcr variability in velum height was obscrvcd, when thc velum w.'1S sllpposcd to bc low (nasal) nnd not suprascgmentally tensc (-strong), or high (oral) and snprasegmentally t.cnse (+strong), than in the combination (-f-nasal), (tstrong) or (oral), (-strong) (Vaissière. 19S3b). The

JACQUELINE VAISSIËRE

tensing is not ollly a function of position in the word, but n]so a nlllction of the relative importnnce of the word as ail iuformation carrying unit in the sentence colltext. If a word is sufficiellt]y strcsscd, it tends to be uttcred with nll higher fundamelltal frequcncy and Ip, t, k] m1\Y bc aspiratcd evcn in final or medial position. As a conscquence, segmentalalld sllprascgrncntal characteristics of speech shou]d be considercd as intimately conllcctcd.

Thercfore, their separate interpretation should be avoided.

Thus EMG data, articulatory positions, and the acoustic signal Illay be interpreted as the rcsult of a combination of bath scglllciltai fcaturcs and the sllprasegmcntal fralllcwork.

The rclntive ililport.'\nce of suprascgmcntal variables as cornparcd with the realizatiolls of the scgmcntal fcaturcs is spc1\kcr-depelldcnt, at least for the velum (Vnissière, 1983b). The Jack of clcctrornyographic, articu]atory, and acoustic invariance for the distinctive fcatures as discusscd in this volume may be partir explained first by context-dcpendent diITerenccs in suprasegmental variables, second by a speaker-dcpcndent way of combining suprasegmcntal influences and the segmentai fcatures, and third by ullpredicta.ble gestures unrelated to the colite nt of speech (sec my comments 011 speech ready gcstures, in Section lOb). Such a Jack of observcd invariance cloes not argue against the possible existcllce of invariance at a more centrallevel, beforc the illtegration of the segmentai and suprasegmental features into the en tire speech event.

CONCLUSION

The acoustic signal is intrinsically highly structured. An important determinant of this structure is the temporal variation of the rclative tenseness of the articulatory processes in the realization a{ a given segmentai featurc. The word plays a great foie in determining this structure and consequently it is a very important unit for interpreting the acoustic correlates of the distinctive features. However, the syllable, the phrase, and the sentence also play important foies; there are no unambiguous criteria to decide which unit is most important.

The word is an obvious building block in speech perception and production, but such a view is not incompatible with the use of the phoncmes as building blocks for constructing words. 80 instead of following Cohen's suggestion by "rcplacing" one building block (the phoneme) by another (the word), it may be more effective to search for an integration of the different acoustic units into a single framework.

Caution is advisable in cxtending the conclusions of studie. of nonsense words to real words. The distinctiveness of a feature at the lexicallevel may play an important role in its mental representation since therc scems to be a Datural tendency to reduce rcdundllncy in the communication process.

The decoding of the specch signal by the listcncr is known to bc 11 complex process involving various types of normnlization.

This proccss is far from being fully understood. But our prcscnt inability to deal with this problem should not lead to ignoring it in n10deling speech undcrstanding, to an Ilndercstirnation of what is actually contmned in the speech signal, and an ovcrcrnphasis of the foIe of syntactic, scmantic, and pragrnatic constrnints. No doubt, wllcn the words arc not in dicta-24a. VARIANCE AND INVARIANCE AT THE WaRD LEVEL 539 tion form but ernbedded in a context, the pronunciation of the word may become "sloppicr" probably because other sources of knowledge constrain the nurnber of possible words in cach position so that the full realization of the complete set of features may bccorne redundant. Nevertheless, literate listeners are perfectly able to discrimilliLte rnemlingful or nonsense words pronounced in isolation and to transcribe them phonetically. One of the problems with invariance is that rescarchers are still confused about what unit should be iilvariant or "lIIore" invariant than the others (featurcs, phonemes, allophones, or even words?) and at which level of the communication proccss (functionru invariance, perceptum cquivalence, acoustic invariance or articulatory control?). Thc question of invariance relllains the central problelll in speech research the importance of which is lIot reduccd by the acknowledgement of the ward as IJ.n important unit. Even ir no invlJ.riance will be round in the future, the search for invariance lIIay provide IJ. useful hypothesis and IJ.n adcquate frlJ.mework to studying speech.

  tcndcncy in contcmporary studics to scparatc thc scgmcntal aud suprascgmcntal information coutaincd in t.llc signal and to intcrprct th cm scparatcly. The rclativc duratiOlI of thc cvcnts ano thc funoamcntal frcquc/lcy co/lt()urs d(~livercd in parall(~1 with thc 8pcctr.ll illformation providc vcry u8cflll i/lformatioll for both decocling thc scgI11(~lIt.11 i/lformation, auo for dct.ccting woro uollllclaric8, word strcss, allo woro rcgroupillg. Ilowcvcr, thc word ha8 intrinsic sc~glllclltai .llld sllprasegmc/ltal Ch.ll..,ctcristics which app(~.'r collcomitantly.

24a. VARIANCE AND INVARIANCE AT THE WaRD LEVEL


Stéphane Rigobert

Noureddine Atalla

Franck C Sgard

Investigation of the convergence of the mixed displacement-pressure formulation for three-dimensional poroelastic materials using hierarchical elements

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

I. INTRODUCTION

The behavior of porous material is classically modeled using empirical formulas for locally reacting materials, [START_REF] Delany | Acoustic properties for fibrous absorbent materials[END_REF] or with accurate one-dimensional ͑1D͒ analytical methods based on Biot poroelasticity theory. [START_REF] Biot | The theory of propagation of elastic waves in a fluid-saturated porous[END_REF][START_REF] Allard | Propagation of Sound in Porous Media. Modeling Sound Absorbing Materials ͑Elsevier[END_REF][START_REF] Bolton | Normal incidence sound transmission through double panel[END_REF][START_REF] Ingard | Notes on Sound Absorption Technology ͑Noise Control Foundation[END_REF] At low frequencies, where the modal behavior of the system is important, several two-dimensional ͑2D͒ and three-dimensional ͑3D͒ finite elements models have been developed. The first models use, with different variations, the displacement fields in the solid and fluid phase of the porous material as variables. [START_REF] Kang | Finite element modeling of isotropic porous materials coupled with acoustical finite elements[END_REF][START_REF] Kang | A finite element model for sound transmission through foam-lined double-panel structures[END_REF][START_REF] Panneton | Numerical prediction of sound transmission through finite multilayer systems with poroelastic materials[END_REF][START_REF] Panneton | An efficient finite element scheme for solving the three-dimensional poroelasticity problem in acoustics[END_REF][START_REF] Coyette | A finite element model for predicting the acoustic transmission characteristics of layered structures[END_REF] They are referred to as ͕u, U͖ type formulations. While accurate, the ͕u, U͖ formulations lead to large frequency dependent systems and thus to cumbersome calculations for complex structures. To alleviate the problem, Atalla et al. [START_REF] Atalla | A mixed displacement-pressure formulation for poroelastic materials[END_REF] presented an exact mixed displacement-pressure ͕u, P͖ finite element formulation in three dimensions for a porous material, based on Biot's poroelastic equations. Debergue et al. [START_REF] Debergue | Boundary conditions for the weak formulation of the mixed ͑u,p͒ poroelasticity problem[END_REF] presented the boundary and coupling conditions for this new formulation. Atalla et al.' model presents the calculation of the response of the porous material in the form of a coupled fluidstructure problem. This formulation has been proven to give accurate results when compared to the ͕u, U͖ formulation. More importantly, by reducing the number of degrees of freedom involved ͑i.e., from six degrees of freedom per node to four͒ and by simplifying the coupling conditions between the porous medium and elastic or fluid media, the approach leads to important gains in memory and computation time requirements.

While the ͕u, P͖ formulation allows for a more efficient modeling to porous elastic media compared to the ͕u, U͖ formulation, yet both techniques suffer from important computational costs in practical applications. Indeed, because of the biphasic nature of poroelastic elements, mesh criterion used for elements describing monophasic media ͑solid, fluid͒ is a priori not valid for poroelastic elements. Using the ͕u, U͖ formulation, Dauchez et al. [START_REF] Dauchez | Convergence of poroelastic finite elements based on biot displacement formulation[END_REF] showed that linear poroelastic elements verify convergence rate of linear monophasic elements, according to each type of Biot wave. Classical mesh criterion, i.e., six linear elements per wavelength, provides a necessary condition to get reliable results. However, for real 3D deformations, the classical criterion gives indications for a minimal mesh but is insufficient because of locking of 3D linear elements and discrepancies in the wavelengths of the displacement fields of the two phases. An important refinement of the mesh is necessary to get satisfactory results. In particular, indicators related to fluid motion have been found very sensitive. As a consequence, the minimal number of elements required is difficult to predict. By nature, convergence of poroelastic elements is rather slower than convergence of either equivalent solid or fluid elements, because of the presence of two different scale phenomena. The discrepancy can be lowered when the behavior of the porous material is dominated by the motion of one phase. As a consequence, this slow convergence leads to large computational costs. Debergue [START_REF] Debergue | De ´veloppement d'une formulation mixte de ´placementpression pour les mate ´riaux poroe ´lastiques ͑Development of a mixed displacement-pressure formulation for poroelastic materials͒[END_REF] found similar results using the ͕u, P͖ formulation. Several recent attempts to alleviate the computing cost of the poroelastic formulations, using linear elements, have been unsuccessful. For example, Sgard et al. [START_REF] Sgard | A modal reduction technique for the finite element formulation of biot's poroelasticity equations in acoustics[END_REF][START_REF] Sgard | A modal reduction technique for the finite element formulation of biot's poroelasticity equations in acoustics applied to multilayered structures[END_REF] have investigated the use of a selective modal analysis to decrease the size of large finite elements models involving poroelastic materials. The method used a dual uncoupled basis of undamped modes associated with the skeleton in vacuo and the fluid phase occupying the volume of the porous material to approximate the ͕u, P͖ variables. However, this method converges poorly due to the importance of the damping brought in by the porous material. The same authors also proposed a mixed wave-finite element approach dedicated to the modeling of porous materials. [START_REF] Sgard | A mixed wave-finite element ap-proach for solving biot's poroelasticity equations in acoustics[END_REF] Starting from the weak ͕u, P͖ formulation for a porous material, the fields for each phase of the porous material were expanded at each node of the finite element mesh in terms of a finite number of wave functions. These functions are chosen to be plane waves solution of Biot's poroelasticity equations. This approach gives accurate results for a single porous material and leads to a considerable reduction of the size of the system to be solved. However, this system is no longer sparse and its construction is time consuming. Another way to solve the computing cost problem, would be to design specific models dedicated to particular applications. For example, Dauchez et al. [START_REF] Dauchez | Dissipation mechanisms in a porous layer bonded onto a plate[END_REF] pointed out that most of the extra-energy dissipation related to the foam coating of a plate was due to structural damping in the solid phase of the porous medium. Thus, the porous material could be reduced to its solid phase in such a configuration. However, this assumption can not be generalized, especially for a coating by a glasswool where viscous effects become quickly predominant as frequency increases. 19 In order to solve the convergence problem, it is proposed in this paper to investigate the performance of the ͕u, P͖ formulation using high-order hierarchical elements. It is known that these elements allow for higher order polynomial representation and thus eliminate problems such as locking. [START_REF] Babus | The p and hp version of the finite element method, basic principle and properties[END_REF][START_REF] Babus | On the rate of convergence of the finite element method[END_REF][START_REF] Szabo | Finite Element Analysis ͑Wiley[END_REF][START_REF] Szabo | Analytical and Computational Models of Shells presented at the Winter An[END_REF][START_REF] Co | Mode ´lisation vibroacoustique dans le domaine des moyennes fre ´quences par e ´le ´ments finis de type p ͑Modeling vibroacoustic problems in the medium frequency range using p elements͒[END_REF] Because of the higher order of approximation of hierarchical elements, a more accurate representation of the different type of deformations and coarser meshes, compared to classical finite elements, can be used. Therefore, a significant reduction of the number of degrees of freedom required for a correct modeling is expected. In the first part of the paper, a variant of the ͕u, P͖ formulation is recalled and the theory of hierarchical elements in this context is introduced. In particular, the treatment of the boundary conditions and excitations is detailed. Validation results are then presented to show the accuracy of the present approach. The performance of the hierarchical poroelastic formulation in terms of the needed number of degrees of freedom is underlined, and can be increased by the choice of different interpolation orders for the basis functions in the solid and fluid phase. The interest of this latter feature is shown and guidelines for the convergence of hierarchical poroelastic elements are presented from the study of a particular material.

II. THEORY

This section concentrates on the forced response of a single isotropic porous material. Biot-Allard equations of poroelasticy are considered for the description of the porous medium. The displacement of the solid phase and the interstitial pressure are chosen as variables.

A. The mixed displacement-pressure formulation

The starting point of the modeling of the porous medium is the weak ͕u, P͖ formulation ͑Atalla et al.͒. [START_REF] Atalla | A mixed displacement-pressure formulation for poroelastic materials[END_REF] The mixed displacement-pressure formulation for a poroelastic material occupying a volume ⍀ p reads

͵ ⍀ p ͓ < ˜S͑u͒:⑀ < S ͑ ␦u͒Ϫ ˜2 u•␦u͔d⍀ ϩ ͵ ⍀ p ͫ h 2 2 ˜22 ٌ P•ٌ␦ PϪ h 2 R ˜P␦ P ͬ d⍀ Ϫ␥ ˜͵⍀ p ␦ٌ͑P•u͒d⍀Ϫ ͵ ץ⍀ p h ͩ 1ϩ Q R ˜ͪ ϫ␦͑ Pu n ͒dSϪ ͵ ץ⍀ p ͓ < ˜t͑u͒•n͔•␦u dS Ϫ ͵ ץ⍀ p h͑U n Ϫu n ͒␦ Pϭ0, ͑1͒
where ץ⍀ p stands for the boundary of the poroelastic domain. u is the solid phase displacement, U and P are the displacement and the pressure of the fluid in the pores, respectively. < ˜S and ⑀ < S are, respectively, the stress and strain tensor of the solid phase in vacuo. The total stress tensor in the porous material < ˜t is related to < ˜S by the relation

< ˜S ϭ < ˜tϪh͓1ϩ(Q ˜/R ˜)͔ P1 < . ˜and ˜22
are the apparent density of the solid and the fluid phase, h is the porosity, ␥ ˜is a coupling factor between the two phases, Q ˜and R ˜are poroelastic coefficients. The expressions of the latter quantities can be found in Ref. 11.

Instead of using directly Eq. ͑1͒, a variant is presented here. The first surface integral in Eq. ͑1͒ can be turned into a volume integral using the second Green formula and given the following mathematical relation:

ٌ•͑ab͒ϭٌa•bϩaٌ•͑b͒. ͑2͒
Hence, assuming that the porous material is isotropic,

IϭϪ ͵ ץ⍀ p h ͩ 1ϩ Q R ˜ͪ ␦͑u n P ͒dS ϭϪ ͵ ⍀ p ٌ ͩ h ͩ 1ϩ Q R ˜ͪ ␦͑u n P ͒ ͪ d⍀ ϭϪh ͩ 1ϩ Q R ˜ͪ ͵ ⍀ p ␦͑u•ٌP͒d⍀ Ϫh ͩ 1ϩ Q R ˜ͪ ͵ ⍀ p ␦͑ Pٌu͒d⍀.

͑3͒

Using Eq. ͑3͒ in Eq. ͑1͒ the following expression of the weak ͕u, P͖ formulation is obtained

͵ ⍀ p ͓ < ˜S͑u͒:⑀ < S ͑ ␦u͒Ϫ ˜2 u•␦u͔d⍀ ϩ ͵ ⍀ p ͫ h 2 2 ˜22 ٌ P•ٌ␦ PϪ h 2 R ˜P␦ P ͬ d⍀ Ϫ ͫ ␥ ˜ϩh ͩ 1ϩ Q R ˜ͪ ͬ ͵ ⍀ p ␦ٌ͑P•u͒d⍀ Ϫh ͩ 1ϩ Q R ˜ͪ ͵ ⍀ p ␦͑ Pٌ•u͒d⍀ Ϫ ͵ ץ⍀ p ͓ < ˜t͑u͒•n͔•␦uϪ ͵ ץ⍀ p h͑U n Ϫu n ͒␦ PdSϭ0. ͑4͒
The weak form of the mixed pressure-displacement formulation ͓Eq. ͑4͔͒ presents two interesting features. First, the symmetric coupling terms between the two phases are expressed in terms of volume integrals over the domain ⍀ p . Second, the boundary integral terms involve the total stress tensor and the displacement flux at the boundary of the domain. Therefore, the application of coupling conditions with another domain, as well as boundary conditions and excitations, is achieved in a simpler manner compared to Eq. ͑1͒. 12

B. Hierarchical elements

In the following, the approximation of Eq. ͑4͒ using 8-nodes parallelipipedic volume hierarchical elements is considered. The so-called ''blending function method'' is generally used to map this parent element into the geometrical element. [START_REF] Szabo | Finite Element Analysis ͑Wiley[END_REF] Any point of the discretized subdomain with global coordinates (x,y,z) is located on a parent element by a set of local coordinates ͑,,͒. The relation between local and global coordinates is given by

ͭ x y z ͮ ϭ ͚ iϭ1 8 ͭ X i Y i Z i ͮ •N i ͑ ,, ͒ϩ ͚ jϭ1 12 ͭ l j x l j y l j z ͮ •l j 1 •l j 2 .

͑5͒

In Eq. ͑5͒, ͗X i ,Y i ,Z i ͘ stands for the global coordinates of node i. The functions N i are the classic linear basis functions used in finite elements and are related to node i on the parent element. In order to consider complex geometries for an element, additional functions related to the edge j are taken into account. These functions correspond to the term on the right-hand side of Eq. ͑5͒. l j 1 and l j 2 are linear functions depending on a single local coordinate ͑different for l j 1 and l j 2 )

, ,or. l j

x , l j y , and l j z are functions depending on the third local coordinate and used for the description of edge j. Further details are given in Ref. 22.

The next step is the interpolation of the fields on the mesh provided for the subdomain. On a given element, the variable q is approximated by

q͑,, ͒ϭ ͚ i N i ͑ , ͒q i ph ϩ ͚ j G j ͑ ,, ͒q j gen . ͑6͒
In Eq. ͑6͒, N i are called node modes and are identical to those used in Eq. ͑5͒. The associated amplitudes q i ph stand for the physical value of q at node i ͑the superscript ph means physical͒. Functions G j are additional shape functions classified in three categories: edge modes, face modes, and internal modes. The associated generalized amplitudes q j gen ͑the superscript gen means generalized͒ do not have simple physical meaning. Like node modes, edge and face modes are associated to a geometrical entity of the mesh. For example, an edge mode on an element is associated with a particular edge j of this element. Its value is 0 on all the other edges of the element. As for internal modes, their value is 0 on all the faces of the element. By considering each category of modes explicitly, Eq. ͑6͒ rewrites q͑,, ͒ϭ ͚

i N i ͑ , ͒q i ph ϩ ͚ j E j ͑ ,, ͒q j gen ϩ ͚ k F k ͑ ,, ͒q k gen ϩ ͚ t I t ͑ ,, ͒q t gen ,

͑7͒

where E j are edge modes, F k are face modes, and I t are internal modes. These functions are chosen to make complete polynomials of ascending order p. Namely, the basis functions of the hierarchical variables are constructed using Legendre polynomials. The number of basis functions depends on the interpolation order p. The selection process, detailed in Ref. 22, leads to the expression of the generalized shape functions; they are given in the Appendix for completeness. Note that the present approach considers a mapping of the parent element into the geometrical element, using only the functions N i in Eq. ͑5͒.

C. Numerical implementation

As previously stated, in the present approach, the weak integral formulation given by Eq. ͑4͒ is discretized using 8-node volume elements. The displacement of the solid phase u and the pressure in the pores P are expanded in terms of node modes and hierarchical shape functions. Hence, the field variables in a given element are written in the following form:

u e ϭ͓N s ͔͕u n ͖ e and p e ϭ͓N f ͔͕p n ͖ e ,

͑8͒

where ͓N s ͔ and ͓N f ͔ are the interpolation matrices on the considered element e. ͕u n ͖ e and ͕p n ͖ e stand for the physical and generalized amplitudes associated with the displacement of the solid phase and the pressure in the pores, respectively.

Discretization of the weak form

Substituting Eq. ͑8͒ into Eq. ͑4͒, one gets

͵ ⍀ p < S ͑ u͒:⑀ < S ͑ ␦u͒d⍀⇒͗␦u n ͓͘K͔͕u n ͖, ͑9͒ ͵ ⍀ p ˜u•␦u d⍀⇒͗␦u n ͓͘M ˜͔͕u n ͖, ͑10͒ ͫ ␥ ˜ϩh ͩ 1ϩ Q R ˜. ͪͬ ͵ ⍀ p ␦u•ٌPd⍀⇒͗␦u n ͓͘C ˜1͔͕P n ͖, ͑11͒ h ͩ 1ϩ Q R ˜ͪ ͵ ⍀ p Pٌ•␦u d⍀⇒͗␦u n ͓͘C ˜2͔͕P n ͖, ͑12͒ h 2 ˜22 ͵ ⍀ p ٌ P•ٌ␦ Pd⍀⇒͗␦ P n ͓͘H ˜͔͕P n ͖, ͑13͒ h 2 R ˜͵⍀ p P␦ Pd⍀⇒͗␦ P n ͓͘Q ˜͔͕P n ͖, ͑14͒
where ͕u n ͖ and ͕P n ͖ represent the solid and fluid phase physical and generalized degrees of freedom of the whole mesh. ͓K͔ and ͓ M ˜͔ are, respectively, the stiffness and mass matrix associated to the solid phase in vacuo. ͓H ˜͔ and ͓Q ˜͔ are, respectively, the kinetic and compressional energy matrices of the fluid phase. ͓C ˜1͔ and ͓C ˜2͔ stand for the coupling between the two phases of the porous medium. The notation ͓˜͔ indicates that the matrices are complex valued and frequency dependent. Note that compared to the implementation of the ͕u, P͖ formulation presented in Ref. 11, two new matrices have to be computed. In the theory of hierarchical elements, the basis functions on an element are built using Legendre polynomials. The computation of the matrices in Eqs. ͑9͒-͑14͒ is performed with the commonly used Gauss-Legendre integration scheme and enough integration points are chosen to ensure a correct approximation. This step is time consuming and is the main drawback of hierarchical elements. [START_REF] Co | Mode ´lisation vibroacoustique dans le domaine des moyennes fre ´quences par e ´le ´ments finis de type p ͑Modeling vibroacoustic problems in the medium frequency range using p elements͒[END_REF] To alleviate that problem, Hinnant [START_REF] Hinnant | Fast method of numerical quadrature for p-version finite element matrice[END_REF] proposed the ''vector quadrature'' integration scheme. However, this optimization procedure is not used in the present approach.

Substituting Eqs. ͑9͒-͑14͒ in the weak formulation leads to the following system to be solved:

ͩ ͓K͔Ϫ 2 ͓ M ˜͔ Ϫ͓C ˜1͔Ϫ͓C ˜2͔ Ϫ 2 ͓C ˜1͔ T Ϫ 2 ͓C ˜2͔ T ͓H ˜͔Ϫ 2 ͓Q ˜͔ ͪͭ u n P n ͮ ϭ ͭ F s F p ͮ ,

͑15͒

where the right-hand side of the equation denotes the loading vector for the porous material. The application of various loading conditions is described in the following paragraph.

The coupled system given by Eq. ͑15͒ is similar to the one obtained in classical finite elements.

Application of boundary and loading conditions

The boundary conditions and the loading terms involved in the ͕u, P͖ formulation were presented by Debergue et al. [START_REF] Debergue | Boundary conditions for the weak formulation of the mixed ͑u,p͒ poroelasticity problem[END_REF] Their application in the context of the theory of hierarchical elements is detailed hereafter. In the present paper, three boundary conditions are considered. The porous material can be ͑i͒ bonded onto a rigid wall, ͑ii͒ guided, or ͑iii͒ free ͑i.e., not subject to any boundary condition͒. In case ͑i͒, the displacement vector of the solid phase is set to 0. In case ͑ii͒, the normal component to the interface of the solid phase displacement is set to 0. In case ͑iii͒, the fluid pressure is set to 0 on the interface ͑note that this condition is approximate and its accuracy is discussed by Debergue 12 ͒. Therefore, for an element on the boundary, applying the boundary conditions mentioned above amounts to setting to 0 one or several components of the fields. This is achieved by constraining the amplitudes of the modes of the nodes, edges and faces lying on the boundary and related to the considered component of the fields. Practically, the coefficients of the matrices in Eqs. ͑9͒-͑14͒ relative to these degrees of freedom are not assembled. Note that, by construction, the value of internal modes are worth 0 on each face of the element and thus are not subject to the above mentioned conditions. In the configurations studied in this paper, two kinds of loads have been considered. First, the load is a rigid piston motion imposed on one side of the porous material. This leads to the following condition on the interface:

u"nϭu 0 , ͑ u n ϪU n ͒ϭ0, ͑16͒
where n is the normal to the interface, and u 0 the amplitude of the piston motion. Second, the porous medium is submitted to an acoustical excitation, which can be modeled as an imposed surface pressure of amplitude p 0 . The condition at the interface is then

pϭp 0 , < t •nϭϪp 0 n. ͑17͒
To be concise, the two considered kinds of excitations lead to an imposed value q 0 of a component q of the fields in the porous material. Considering the expression of q given by Eq. ͑6͒ the following conditions must be imposed on the interface:

͕q ph ͖ϭ͕q 0 ͖ and ͕q gen ͖ϭ0.

͑18͒

Equation ͑18͒ means that the values of q have to be set to q 0 at the nodes on the interface. Besides, the amplitudes of the hierarchical shape functions related to the edges or faces located on the interface are set to 0. By definition, the amplitudes of the internal modes are null on the sides of an element. The conditions of imposed degrees of freedom expressed by Eq. ͑18͒ are taken into account in the discretized weak formulation using discrete Lagrange multipliers.

In addition, when a surface pressure with amplitude p 0 is imposed, Eq. ͑17͒ leads to an elementary force vector given by

͕F s ͖ elem ϭϪ ͵ ץ⍀p elem p 0 n•͓N s ͔dS, ͑19͒
where ͓N s ͔ is the column vector containing the shape function related to the face if the considered element subject to the excitation. Note that only node modes, edge modes with interpolation order pϭ2, and face modes with interpolation order pϭ4 are taken into account in ͓N s ͔. In fact, the other modes of analytical expression G verify the following relation by construction:

͵ ץ⍀p elem G dSϭ0. ͑20͒
The proof for this assertion is given in the Appendix. Note that for all the excitations considered in this paper, ͕F p ͖ ϭ0 in Eq. ͑15͒.

Vibroacoustics indicators

In the result section, two vibroacoustic indicators are considered. For the solid phase of the porous material, the mean square velocity along the three directions is computed. Practically, the mean square velocity along direction i is given by

͗v i 2 ͘ϭ 2 2⍀ p ͵ ⍀ p ͉v i ͉ 2 d⍀ ϭ 2 2⍀ p ͗u n ph u n gen ͘*͓M i ͔ ͭ u n ph u n genͮ , ͑21͒
where ͗u n ph u n gen ͘ is the row vector containing the physical and generalized amplitudes relative to the solid phase displacement along direction i. ͓ M i ͔ contains the coefficients of a partition of the mass matrix ͓ M ˜͔ divided by ˜and related to the appropriate degrees of freedom. In Eq. ͑21͒, ͑*͒ means complex conjugate, and ⍀ p is the volume of the poroelastic domain.

For the fluid phase, the mean square pressure is computed using the following formula:

͗P 2 ͘ϭ 1 2⍀ p ͵ ⍀ p ͉P͉ 2 d⍀ ϭ 1 2⍀ p ͗P n ph P n gen ͘*͓Q 1 ͔ ͭ P n ph P n genͮ , ͑22͒
where ͗P n ph P n gen ͘ is the row vector of the physical and gen- eralized amplitudes related to the pressure in the pores. ͓Q 1 ͔ is the compression matrix of the fluid phase ͓Q ˜͔ divided by h 2 /R ˜.

III. NUMERICAL EXAMPLES

In the following, the validation of the ͕u,P͖ formulation using hierarchical elements is numerically assessed. The present approach is compared to a finite element code developed at the Universite ´de Sherbrooke and based on classical finite elements. This latter code has been validated both numerically and experimentally, elsewhere. [START_REF] Atalla | A mixed displacement-pressure formulation for poroelastic materials[END_REF][START_REF] Debergue | Boundary conditions for the weak formulation of the mixed ͑u,p͒ poroelasticity problem[END_REF] The accuracy of the results obtained by the use of hierarchical elements is underlined. The reduction of the number of degrees of freedom provided by the use of high order polynomials for the basis functions is highlighted, thus illustrating the performances of hierarchical poroelastic elements. An emphasis is made on this latter aspect by showing the benefit of using different interpolation orders for the basis functions of each phase. Finally, the convergence of the hierarchical poroelastic elements is investigated in the case of a particular material with several boundary conditions and excitations.

The studied problem consists of a single porous material with a rear face bonded onto a rigid wall ͑see Fig. 1͒. The boundary conditions on the lateral faces are either ͑i͒ free, or ͑ii͒ bonded. Two kinds of loads are used in the following, depending on the boundary conditions on the edges of the porous medium. In case ͑i͒, the porous medium is subjected to a rigid piston motion with a 10 Ϫ3 m amplitude normal displacement ͓hereby referred to as configuration ͑i͔͒. In case ͑ii͒ a surface pressure of amplitude 1 Pa is imposed on the front face ͓hereby referred to as configuration ͑ii͔͒.A5c m thick porous material with lateral dimensions 0.35 m ϫ0.22 m is considered. Among the vibroacoustic indicators presented in the preceding section, only the mean square velocity along z axis ͑denoted ͗V z 2 ͘) and the mean quadratic pressure ͑denoted ͗P 2 ͘) are studied. Actually the mean qua- dratic velocities along the lateral dimensions of the porous material are of limited interest in vibroacoustic applications.

A. Validation

In this subsection, the accuracy and the performance of the present approach is numerically assessed. In the following, the results obtained with hierarchical elements and finite elements are compared to a reference value. In the present paper, this reference is obtained by increasing considerably the interpolation orders for a fixed mesh when hierarchical elements are used, in order to be sure that convergence is reached. The same results can be obtained with classical finite elements but require a very refined mesh and thus huge computational resources. Hence, the classical finite element implementation has not been used for obtaining the reference values. In the following, the interpolation orders for the basis functions given a fixed mesh for the present approach or the classical finite element mesh are said to be acceptable when the response of the porous material is predicted within a tolerance from the reference response. This tolerance is fixed to 0.5 dB from the reference curve and 5% in the location of resonance peaks.

The first validation test consists in configuration ͑i͒. The material studied is a wool, UGW3, with characteristics given in Table I. ͗P 2 ͘ is computed as a function of frequency and is actually the indicator which proves to have the most difficulties to converge. The results obtained with the classical code and the present approach with different meshes are presented in Fig. 2. Note that the maximum difference between the different meshes, including classical finite elements, and the reference curve is less than 0.5 dB which is within the defined convergence criterion. For the modeling with hierarchical elements, parameters for the mesh, namely the number of elements and the order of the basis functions, are given in the following form n x* n y* n z elements (p s Ϫp f ) where n x , n y , n z denote the number of elements used in directions x, y, z, respectively, and p s , p f denote the interpolation order of the basis functions for the solid phase and the fluid phase, respectively. Figure 2 shows that the present approach leads to accurate results for all tested meshes. To test the ability of the hierarchical poroelastic element to model different materials in various configurations a foam, FM2, with characteristics given in Table I is considered in configuration ͑ii͒. ͗V z 2 ͘ is computed as a function of frequency; this indicator has the most difficulties to converge. The results are presented in Fig. 3. Once again, accurate results are obtained with the hierarchical poroelastic elements for different meshes.

The excellent results for the two previous validation tests show the validity of the present approach. In addition, hierarchical poroelastic elements prove to be very performant. Tables II and III summarize the number of degrees of freedom ͑dof͒ required to ensure convergence, either with classical finite elements or hierarchical elements. Each of the two configurations studied above are considered. For the present approach, the number of degrees of freedom is given for different meshes and the corresponding required interpolation orders. Actually, the present approach allows for an important decrease of the number of dof required to satisfy the convergence criterion defined previously ͑0.5 dB from the reference curve and 5% accuracy in the location of resonance peaks͒. Moreover, the best performance is achieved using a coarse mesh and a high interpolation order for the basis functions. This trend has already been pointed out for the modeling of elastic materials using hierarchical elements in static problems. [START_REF] Szabo | Finite Element Analysis ͑Wiley[END_REF] As seen in Tables II andIII, the interpolation order for the basis functions can be chosen to be different for the two phases of the porous material. Actually, the solid and the fluid phase exhibit different physical behaviors depending on the studied configuration. This feature of the hierarchical poroelastic elements is of great interest since convergence is ensured by choosing the adequate interpolation order for each phase separately. On the other hand, convergence in classical finite elements is ensured by refining the mesh, i.e., a refined description of both phases of the porous material. Hence, the former approach provides an efficient way for the reduction of the number of degrees of freedom required for a correct modeling of the material.

B. Convergence of the vibroacoustic indicators

In the following, foam FM2 is studied in configurations ͑i͒ and ͑ii͒. The geometry of the porous sample is identical to the one considered previously. The frequency range of interest is ͓0 Hz; 500 Hz͔. Biot's theory of poroelasticity 2 indicates that two compression waves ͑denoted P 1 and P 2 in the following͒ and a shear wave ͑denoted S͒ can propagate simultaneously in a porous medium. For each kind of wave, the ratio of the fluid phase displacement to the solid phase displacement can be computed. [START_REF] Allard | Propagation of Sound in Porous Media. Modeling Sound Absorbing Materials ͑Elsevier[END_REF] These ratio at 500 Hz for foam FM2 are given in Table IV. Table IV indicates that waves P 1 and S propagate in both phases. On the contrary, wave P 2 propagates mainly in the fluid phase, so that the two phases of the porous material can be considered as decoupled as far as wave P 2 is concerned. This point strongly influences the behavior of the porous materials according to the configuration studied. Foam FM2 has been previously studied in configuration ͑ii͒. Excellent values for ͗V z 2 ͘, repre- sented in Fig. 3, can be obtained with low interpolation order for the basis functions of the fluid phase. This means that a coarse description of the fluid phase suffices for the convergence of this indicator related to the solid phase. Indeed, this porous material in configuration ͑ii͒ is subjected to an acoustical excitation acting on both phases. The pressure in the pores and the stress tensor in the solid phase in vacuo are imposed at the interface where the material is excited. Because wave P 2 mainly propagates in the fluid phase, the fluid phase has a weak influence on the solid phase and thus does FIG. 3. ͗V z 2 ͘ for a foam FM2 in configuration ͑ii͒. not require a very accurate description for the convergence of ͗V z 2 ͘: the motion in the solid phase is due to wave P 1 and S. In comparison, indicator ͗P 2 ͘ for foam FM2 in the same configuration is represented in Fig. 4. Different interpolation orders have been used with a fixed 2ϫ2ϫ1 elements mesh for the present approach. Compared to ͗V z 2 ͘, ͗P 2 ͘ requires a noticeable lower interpolation order for the basis functions in the solid phase. This indicates that a coarse representation of the solid phase suffices for the convergence of this latter indicator. The influence of waves P 1 and S on the motion of the fluid phase is not significant. Actually the motion of the fluid phase is mainly due to wave P 2 .

Next, consider foam FM2 in configuration ͑i͒. As previously, ͗P 2 ͘ is plotted as a function of frequency in Fig. 5. It can be observed that a higher interpolation order for the basis functions of the solid phase in comparison with configuration ͑ii͒ is required for the convergence of ͗P 2 ͘. In this configu- ration, the porous is submitted to a rigid piston motion. For that kind of excitation, solid phase is directly excited and its motion is strongly transmitted to the fluid when the amplitude ratios are important for waves P 1 and S. Hence, in configuration ͑i͒ the motion of the fluid phase is mainly due to these latter waves. In conclusion, the two phases of a porous material can exhibit very different behaviors according to the excitation. This induces a convergence which highly depends on the configuration studied for the considered vibroacoustic indicators, namely ͗V z 2 ͘ or ͗P 2 ͘. Hence, it is difficult to conclude about the convergence of a set of indicators in a general manner.

C. Investigation of the convergence of hierarchical poroelastic elements

For 3D deformations, the classical convergence criterion gives indications for a minimal mesh but is insufficient because of locking of 3D linear elements and discrepancies of the fields in the two phases of the porous material. As stated in Sec. III A, the use of the theory of hierarchical elements together with the ͕u, P͖ formulation enables one to solve these latter problems and to get smaller linear systems to solve.

In this section, the convergence of hierarchical poroelastic elements is investigated for a particular material, wool UGW3. The aim of this study is to determine the interpolation order for the basis functions in the two phases of the porous material that allows for convergence given the element size. Initially, the motivation for this study was to de-rive a convergence criterion that could be reused for other porous materials. However, it appears that the results for wool UGW3 cannot be directly transposed and that only a convergence study performed on different materials would enable one to draw general trends for the convergence of hierarchical poroelastic elements. This task implies a huge amount of calculations. Therefore, only preliminary results are given here. This convergence study is however helpful to highlight the influence of the coupling between the two phases of wool UGW3 as it is shown thereafter and to confirm the observations made on foam FM2 in Sec. III B.

In the following, the mesh in the lateral dimensions of a 3D porous sample is studied. As a consequence, the number of elements along the thickness is chosen large enough so that only the lateral mesh influences the results. For a chosen indicator, the interpolation order in both phases can be represented as a function of the element size. In order to have nondimensional data in the abscissa, the element size is scaled to the wavelength of wave P 1 , P 2 or S at the highest frequency in the frequency range of interest. These wavelengths at 500 Hz are given in Table V. As these wavelengths are rather short, a wide range of element size to wavelength ratio is provided. The study that follows is carried out for each of the two configurations ͑i͒ and ͑ii͒ depicted in Sec. III A and for each vibroacoustic indicator (͗V z 2 ͘ and P 2 ) ac- cording to remarks in Sec. III B.

Configuration (i):

The interpolation orders for the solid phase ͑respectively, fluid phase͒ basis functions that ensure the convergence of the two indicators ͗P 2 ͘ and ͗V z 2 ͘ are plotted as a function of the ratio of the element size to P 1 , where P 1 denotes the wavelength of wave P 1 at 500 Hz.

The results are represented in Fig. 6 ͑respectively, Fig. 7͒.In these figures, fit curves are provided for clarity. Note that the wavelength of another kind of wave, e.g., S wave or P 2 , can be chosen to scale the element size simply by multiplying the FIG. 4. ͗P 2 ͘ for foam FM2 in configuration ͑ii͒. Considering a fixed mesh, it is seen in Figs. 6 and7 that the convergence of indicator ͗P 2 ͘ is not solely governed by the interpolation order for the basis functions in the fluid phase. Indeed, a sufficiently high interpolation order for the basis functions in the solid phase is required to ensure the convergence of ͗P 2 ͘. This interpolation is even higher than the one needed for the convergence of ͗V z 2 ͘. The opposite is observed regarding the interpolation order for the basis functions in fluid phase. This is because of important coupling phenomena between the two phases when the solid phase of the porous material is directly excited. Configuration (ii): Like in configuration ͑i͒, the interpolation orders required for the convergence of the two indicators ͗P 2 ͘ and ͗V z 2 ͘ are plotted as a function of the ratio of the element size to P 1 . These interpolation orders for the solid and the fluid phase are represented in Figs. 8 and9, respectively. As previously, fit curves are provided for clarity. When considering the basis functions in the solid phase, Fig. 8 shows that for a fixed element size to P 1 , the interpolation order required for the convergence of ͗P 2 ͘ is as high as the one for ͗V z 2 ͘, except for low size element to wavelength ratio. For the basis functions of the fluid phase, the interpolation order ensuring the convergence of ͗V z 2 ͘ is smaller than the one ensuring the convergence of ͗P 2 ͘.

Hence, from these observations, coupling phenomena between the two phases of the porous material appear to influence weakly the response of the porous material. This latter phenomenon can be explained by the fact that the porous material has bonded edges and is submitted to an acoustical excitation. To sum up the results found in this particular configuration, given a fixed mesh, the interpolation orders for the basis functions in the two phases of the porous material chosen to ensure simultaneously the convergence of ͗P 2 ͘ and ͗V z 2 ͘ exhibit the following feature: regarding the basis functions in the solid phase ͑respectively, fluid phase͒, the choice of the interpolation order ensuring the convergence of ͗V z 2 ͘ ͑respectively, ͗P 2 ͘) allows for the convergence of both indicators.

IV. CONCLUSION

The implementation of the ͕u, P͖ formulation for porous materials using the theory of hierarchical elements has been presented. The use or high order polynomials for the basis functions enables one to solve some of the problems met with linear poroelastic elements, among them numerical locking and the presence of two different scale phenomena. Hence, hierarchical poroelastic elements allows for the prediction of the forced response of a porous material using a reduced number of degrees of freedom. The reduction of the number of unknowns is all the more effective than a coarse mesh with high interpolation orders is provided. Besides, the performance of hierarchical elements is enhanced by the use of different interpolation orders for the basis functions of each phase separately. Considering a particular material, a preliminary study has been conducted in order to link the interpolation order of the basis functions in the two phases of the porous material to the ratio of the element size to the wavelength. It has been found that for a mechanical excitation, each phase happens to influence strongly the vibratory state of the other phase, showing the importance of the coupling between the two phases of the porous material. Regarding interpolation order for the basis functions in the solid phase, this implies, for example, to chose an interpolation order for the convergence of ͗P 2 ͘ larger than the one for ͗V z 2 ͘. This interaction phe- nomenon proves to be of least importance for an acoustical excitation. In this case, considering the basis functions in the solid phase, the interpolation order for the convergence of ͗V z 2 ͘ is found to be as large as for ͗P 2 ͘. To derive more general trends for the convergence of hierarchical poroelastic elements, further works involving the study of other kinds of materials are required. However, the presented examples show the difficulty of deriving a general meshing criterion for porous-elastic materials. More importantly, they highlight the importance of conducting a serious convergence study before any vibroacoustic study regarding the behavior of porous-elastic materials, based on using finite element codes.

APPENDIX

The shape functions involved in the approximation of the fields on one element ͓Eq. ͑6͔͒ are defined on a parent element, which is an 8-noded brick element represented in Fig. 10. On that element, the geometrical entities ͑nodes, edges, and faces͒ are assigned a number. From this follows the definition of the basis shape functions.

Node modes

A node mode of analytical expression N i , is related to node i on the parent element. N i is a linear function of the local coordinates and has value 1 at node i and 0 at the other nodes. Node modes are given by

N i ͑ ,, ͒ϭ 1 8 ͑ 1Ϫ ͒͑ 1Ϫ ͒͑ 1Ϫ ͒, N 2 ͑ ,, ͒ϭ 1 8 ͑ 1ϩ ͒͑ 1Ϫ ͒͑ 1Ϫ ͒, N 3 ͑ ,, ͒ϭ 1 8 ͑ 1ϩ ͒͑ 1ϩ ͒͑ 1Ϫ ͒, N 4 ͑ ,, ͒ϭ 1 8 ͑ 1Ϫ ͒͑ 1ϩ ͒͑ 1Ϫ ͒, N 5 ͑ ,, ͒ϭ 1 8 ͑ 1Ϫ ͒͑ 1Ϫ ͒͑ 1ϩ ͒, N 6 ͑ ,, ͒ϭ 1 8 ͑ 1ϩ ͒͑ 1Ϫ ͒͑ 1ϩ ͒, N 7 ͑ ,, ͒ϭ 1 8 ͑ 1ϩ ͒͑ 1ϩ ͒͑ 1ϩ ͒, N 8 ͑ ,, ͒ϭ 1 8 ͑ 1Ϫ ͒͑ 1ϩ ͒͑ 1ϩ ͒.

Edge modes

The edge i linking two nodes of the parent element is associated one or several edge modes. One particular edge mode of analytical expression E p has zero value on all the edges of the parent element except on edge j. On this particular edge, its expression is a polynomial with order p. The general expression of the edges modes for pу2 are given by

E p

͑A1͒

The relations verified by functions p are given further on.

Face modes

Each face i of the parent element is associated with face modes which are zero on all the faces of the parent element except for face i. Considering the interpolation order pу4, one or several couples of integers (k 1 ,k 2 ) verifying k 1 ϩk 2 ϭp are constituted. Note that the couple (k 1 ,k 2 ) is different from the couple (k 2 ,k 1 ). Each couple is associated a mode related to face i. The general expression for this mode F k 1 k 2 is given by

F k 1 k 2 1265 ϭ 1 2 k 1 ͑ ͒ k 2 ͑ ͒͑ 1Ϫ ͒, F k 1 k 2 2376 ϭ 1 2 k 7 1 ͑ ͒ k 2 ͑ ͒͑ 1Ϫ ͒, F k 1 k 2 1234 ϭ 1 2 k 1 ͑ ͒ k 2 ͑ ͒͑1Ϫ͒, F k 1 k 2 4158 ϭ 1 2 k 1 ͑ ͒ k 2 ͑ ͒͑ 1ϩ ͒, F k 1 k 2 3487 ϭ 1 2 k 1 ͑ ͒ k 2 ͑ ͒͑ 1ϩ ͒, F k 1 k 2 5678 ϭ 1 2 k 1 ͑ ͒ k 2 ͑ ͒͑1ϩ͒,
where the superscript klmn states for the face defined by the nodes k, l, m, and n.

Internal modes

For an interpolation order pу6, sets of three integers (k 1 ,k 2 ,k 3 ) verifying k 1 ϩk 2 ϩk 3 ϭp are considered. Note that the set (k 1 ,k 2 ,k 3 ) is different from the set (k 1 ,k 3 ,k 2 ) and (k 2 ,k 3 ,k 1 ). Each set is associated with an internal mode of analytical expression I k 1 k 2 k 3 given by

I k 1 k 2 k 3 ϭ k 1 ͑ ͒• k 2 ͑ ͒• k 3 ͑ ͒.

Features of functions k

The functions k introduced in the preceding sections satisfy some interesting properties. First, by construction, k ͑ Ϫ1 ͒ϭ k ͑ 1 ͒ϭ0. ͑A2͒

Also the first derivative of these functions are orthogonal, namely,

͵ Ϫ1 1 m Ј ͑ ͒ n Ј͑͒dϭ␦ mn .

͑A3͒

Finally, let's consider the analytical expression of a particular edge mode with an interpolation order k, for example, E k 12 . The integration of this mode on the parent element face pertaining to the nodes 1, 2, 3, 4 leads to

͵ 1234 E 12 dSϭ ͵ Ϫ1 1 ͑ 1Ϫ ͒d• ͵ Ϫ1 1 k ͑ ͒d. ͑A4͒
For kϾ2, the function k verify the relation

͵ Ϫ1 1 k ͑ ͒dϭͱ 2 2 jϩ1 ͫ 1 2kϩ1 ͓ kϩ1 ͑ 1 ͒Ϫ kϩ1 ͑ Ϫ1 ͔͒ Ϫ 1 2kϪ3 ͓ kϪ1 ͑ 1 ͒Ϫ kϪ1 ͑ Ϫ1 ͔͒ ͬ . ͑A5͒
Given Eq. ͑A2͒, the value of this integral is 0. Hence all the edge modes with an interpolation order kϾ2 satisfy Eq. ͑20͒. The extension of this relation for the face modes with interpolation order pϾ4 is straightforward.

FIG. 1 .

 1 FIG. 1. Configuration of the problem.

FIG. 5 .

 5 FIG. 5. ͗P 2 ͘ for foam FM2 in configuration ͑i͒.

FIG. 6 .

 6 FIG.6. Convergence criterion for the solid phase. Configuration ͑i͒. Given an element size for the mesh, an interpolation order in the solid phase that suffices for the convergence of ͗P 2 ͘ ͑respectively, ͗V z 2 ͘) is determined by numerical experiments. ᭝ ͑respectively, ᭺͒ stands for the corresponding couple of values interpolation order-element size to wavelength ratio.

FIG. 8 .

 8 FIG. 8. Convergence criterion for the solid phase. Configuration ͑ii͒. Given an element size for the mesh, an interpolation order in the solid phase that suffices for the convergence of ͗P 2 ͘ ͑respectively, ͗V z 2 ͘) is determined by numerical experiments. ᭺ ͑respectively, ᭝͒ stands for the corresponding couple of values interpolation order-element size to wavelength ratio.

TABLE II

 II

				2-2-1	3-3-1
			7-2	4-2	3-2
	Number of dof	2144	70	144	244

. Convergence of ͗V z 2 ͘ for FM2 in configuration ͑ii͒. Nb of elements interpolation order ͑solid-fluid͒ Classical code 1-1-1

TABLE I .

 I Characteristics of the materials.

					⌳	⌳Ј	S	N
	Material	h	kN s/m 4 ␣ ϱ	m	m	kg/m 3	kPa
	FM2	0.9	25	7.8 226	226	300 286 0.4 0.265
	UGW3 0.95	25	1.4	93.2 93.2 600	21 0	0.05

FIG.

2

. ͗P 2 ͘ for glasswool UGW3 in configuration ͑i͒.

TABLE III .

 III Convergence of ͗P 2 ͘ for UGW3 in configuration ͑i͒.

	Nb of elements		4-3-2	5-3-2	6-4-3	9-7-4
	interpolation order ͑solid-fluid͒	Classical code	3-2	3-2	3-2	2-2
	Number of dof	5214	895	1094	2397	4828

TABLE IV .

 IV Amplitude ratios for foam FM2 at 500 Hz. Amplitude ratio 0.998ϩ0.161ϫi Ϫ10.43ϩ1.70ϫi 0.97Ϫ1.14ϫi abscissa in Figs.6 and 7by the factor P 1 / S or P 1 / P 2 .

	Kind of wave	P 1	P 2	S

TABLE V .

 V Wave numbers for the different kinds of waves in wool UGW3 at 500 Hz.

	Kind of wave	P 1	P 2	S
	Wavelength ͑m͒	7.237ϫ10 Ϫ2	1.962ϫ10 Ϫ1	5.199ϫ10 Ϫ2

  where the superscript kl indicates that the considered edges is defined by node k and node l on the parent element. Functions p are built using orthogonal Legendre polynomials L p and are obtained by the formula

	12 ϭ 1 4 p ͑ ͒͑ 1Ϫ ͒͑ 1Ϫ ͒,
	E p 23 ϭ 1 4 p ͑ ͒͑1ϩ͒͑1Ϫ͒,
	E p 34 ϭ 1 4 p ͑ ͒͑ 1ϩ ͒͑ 1Ϫ ͒,
	E p 41 ϭ 1 4 p ͑ ͒͑1Ϫ͒͑1Ϫ͒,
	E p 15 ϭ 1 4 p ͑ ͒͑ 1Ϫ ͒͑ 1Ϫ ͒,
	E p 26 ϭ 1 4 p ͑ ͒͑ 1ϩ ͒͑ 1Ϫ ͒,
	E p 37 ϭ 1 4 p ͑ ͒͑ 1Ϫ ͒͑ 1ϩ ͒,
	E p 48 ϭ 1 4 p ͑ ͒͑ 1ϩ ͒͑ 1ϩ ͒,
	E p 56 ϭ 1 4 p ͑ ͒͑ 1Ϫ ͒͑ 1ϩ ͒,
	E p 67 ϭ 1 4 p ͑ ͒͑1ϩ͒͑1ϩ͒,
	E p 78 ϭ 1 4 p ͑ ͒͑ 1ϩ ͒͑ 1ϩ ͒,
	E p 85 ϭ 1 4 p ͑ ͒͑1Ϫ͒͑1ϩ͒,

p ͑ ͒ϭ 1 ͱ2 ͑ 2 jϪ1 ͒ ͑ L p ͑ ͒ϪL pϪ2 ͑ ͒͒ with pу2.


Noureddine Atalla

Raymond Panneton

Patricia Debergue

A mixed displacement-pressure formulation for poroelastic materials

Recently, finite element models based on Biot's displacement (u,U) formulation for poroelastic materials have been extensively used to predict the acoustical and structural behavior of multilayer structures. These models while accurate lead to large frequency dependent matrices for three-dimensional problems necessitating important setup time, computer storage and solution time. In this paper, a novel exact mixed displacement pressure (u,p) formulation is presented. The formulation derives directly from Biot's poroelasticity equations. It has the form of a classical coupled fluid-structure problem involving the dynamic equations of the skeleton in vacuo and the equivalent fluid in the rigid skeleton limit. The governing (u,p) equations and their weak integral form are given together with the coupling conditions with acoustic media. The numerical implementation of the presented approach in a finite element code is discussed. Examples are presented to show the accuracy and effectiveness of the presented formulation.

INTRODUCTION

The vibroacoustic performance of finite multilayer systems containing poroelastic materials is of utmost importance for noise control in automobiles, aircraft and several other engineering applications. In the absence of absorbing materials, the vibroacoustic response of complex multilayer structures are classically modeled using the finite element and the boundary element methods. To account for absorbing media, finite element formulations for sound absorbing materials have been developed. They range from simple approaches ͑impedance techniques: Beranek and Ver, 1992; equivalent fluid models [START_REF] Craggs | A finite element for rigid porous absorbing materials[END_REF][START_REF] Panneton | A finite element formulation for the vibro-acoustic behaviour of double plate structures with cavity absorption[END_REF] to sophisticated approaches based on the Biot theory ͑Kang and Bolton, 1995;[START_REF] Johansen | Finite element method for predicting the acoustical properties of porous samples[END_REF][START_REF] Coyette | A finite element model for predicting the acoustic transmission characteristics of layered structures[END_REF][START_REF] Panneton | Mode ´lisation nume ´rique par e ´le ´ments finis des structures complexes absorbantes[END_REF]Atalla, 1996, 1997a͒. For the general case of multilayer systems containing poroelastic media, the use of the sophisticated approaches becomes necessary. These approaches are mainly based on the classical displacement (u,U) formulation of Biot's poroelasticity equations ͑Biot, 1956͒. However, it has been shown that while accurate the (u,U) formulation has the disadvantage of requiring cumbersome calculations for large finite element models and spectral analyses ͑Panneton, 1996; [START_REF] Panneton | Mode ´lisation nume ´rique par e ´le ´ments finis des structures complexes absorbantes[END_REF]Atalla, 1996, 1997a͒. To alleviate these difficulties, a mixed formulation using the solid phase displacement and the interstitial fluid pressure as variables seems more suitable for large size finite element models since it requires only four degrees of freedom per node. However, few efforts have been devoted to this idea. Lately, a simplified mixed (u, p) formulation was proposed by Go ¨ransson ͑1995͒. However, in his formulation, Go ¨ransson neglected the strain couplings between the fluid and solid phases of the porous material. He solved the one-dimensional wave propagation in porous media using a weighted residual formulation based on Galerkin's method. In geomechanics, several authors have investigated mixed formulations ͑see Simon and Zienkiewics, 1986 for references͒. However, all these formulations are based on assumptions that neglect inertia coupling, which is not acceptable for vibroacoustic applications. Based on these investigations, it is widely believed in acoustics that a (u,p) formulation is not possible without limiting assumptions.

In this paper, an exact (u,p) formulation recently developed by the authors is presented ͑Atalla et al., 1996͒. In this new formulation, the classical Biot-Allard equations are rewritten in terms of the solid phase macroscopic displacement vector and the interstitial fluid phase macroscopic pressure. The new coupled equations have the advantage of recasting the poroelasticity equations in the classical form of a fluidstructure coupled problem wherein the coupling is between an equivalent elastodynamic equation for the skeleton in vacuo and an equivalent Helmholtz equation for the fluid in the limit of a rigid skeleton. In the first part of the paper, the classical formulation is recalled. Next, the new formulation is derived followed by its numerical implementation. The coupling conditions between the presented formulation and acoustic formulation are given. Finally, examples demonstrating the accuracy and efficiency of the approach are presented.

I. THEORY

The problem of interest deals with the prediction of the vibroacoustic response ͑dynamic and acoustic response͒ of multilayer structures made up of elastic, poroelastic and acoustic media. The poroelastic material may be bonded or unbonded to the structure. The classical assumptions concerning linear acoustic, elastic, and poroelastic wave propagation are assumed ͑Allard, 1993͒. Also, the air contained in the porous medium is initially at rest.

A. The Biot's displacement formulation

The modeling of the poroelastic material is based on the Biot's poroelasticity equations ͑Biot, 1956; Allard, 1993͒:

1 div = s ϭ 11 u ¨ϩ 12 U ¨ϩb ˜͑u ˙ϪU ˙͒ ͑1͒ div = f ϭ 22 U ¨ϩ 12 u ¨Ϫb ˜͑u ˙ϪU ˙͒,
where the tilde symbol indicates that the associated physical property is complex and frequency dependent. In Eq. ͑1͒, u and U denote the solid and fluid macroscopic displacement vectors, respectively. These displacements are average values in the sense of the Biot theory, i.e., the average volume displacements per unit area cross section. Densities 11 and 22 are mass coefficients which take into account the fact that the relative flow through the pores is not uniform. They are related to the mass density of the material from which the skeleton is made, s , and to the mass density of the intersti- tial fluid, 0 ,b y:

11 ϭ͑1Ϫh ͒ s Ϫ 12 , 22 ϭh 0 Ϫ 12 , ͑2͒
where h is the porosity of the poroelastic material. The coefficient 12 accounts for the interaction between the inertia forces of the solid and fluid phases; it is function of the tortuosity of the material ͑Allard, 1993͒. = s and = f denote the partial stress tensors associated with the skeleton particle and the macroscopic fluid particle, respectively. Tensor = f is related to the averaged pressure prevailing in the saturating fluid by = f ϭϪhp1 =, while tensor = s is given by (1Ϫh) times the averaged stress tensor prevailing in the matrix of the poroelastic material. Note that = s and = f are linearly related to the partial strain tensors = s and = f prevailing in the skeleton and the interstitial fluid, respectively. These relations define the stiffness coupling between the solid and the fluid phases and account for dissipation due to thermal exchanges between both phases in the material. Also, they include the structural damping associated with the frame. Finally, b ˜is a frequency dependent viscous damping coefficient accounting for viscous interaction forces. It can be related to the macroscopic flow resistivity of the porous material ͑Allard, 1993͒. A detailed discussion about the field quantities occurring in the Biot theory may be found in De Vries ͑1989͒ and Geerits ͑1996, 1997͒.

Using the solid and fluid displacement vectors (u,U)a s primary variables and assuming harmonic oscillations (e jt ), the weak integral form of the poroelasticity equations, Eq. ͑1͒, reads ͑Panneton and Atalla, 1996͒:

͵ ⍀ p ͑ = s ͑ u,U ͒: = s ͑ ␦u͒Ϫ 11 2 u.␦uϪ 12 2 U.␦u ϩ jb ˜͑uϪU ͒.␦u)d⍀ Ϫ ͵ ץ⍀ p ␦u.͓ = s ͑ u,U ͒.n͔dSϭ0, ͑3͒ ͵ ⍀ p ͑ = f ͑ u,U ͒: = f ͑ ␦U͒Ϫ 22 2 U.␦UϪ 12 2 u.␦U Ϫ jb ˜͑uϪU ͒.␦U)d⍀ Ϫ ͵ ץ⍀ p ␦U.͓ = f ͑ u,U ͒.n͔dSϭ0, ᭙͑␦u,␦U ͒,
where ␦u and ␦U denote admissible variations of u and U, respectively, and where ⍀ p and ץ⍀ p denote the medium domain and its boundary.

For the finite element implementation of Eq. ͑1͒,a n analogy with the three-dimensional elastic solid elements is used; however, this time, six degrees of freedom per node are used. They account for three displacement components of the solid phase and three displacement components of the fluid phase. However, because of the viscous and thermal dissipation mechanisms, the system's matrices are frequency dependent. In consequence, for large 3D multilayer structures, this formulation has the disadvantage of requiring cumbersome calculations for large finite element models and spectral analyses. To alleviate these difficulties, a mixed displacement pressure formulation is presented.

B. The mixed displacement pressure formulation

In the following, an exact mixed formulation valid for harmonic motion is given. It derives directly from Biot poroelasticity equations, Eq. ͑1͒. The developments assume that the porous material properties are homogeneous. Note that a more general time domain formulation valid for anisotropic materials is given by Gorog et al. ͑1997͒.

For harmonic motion (e jt ), introducing the effective densities

˜11 ϭ 11 ϩ b j ˜22 ϭ 22 ϩ b j ͑4͒ ˜12 ϭ 12 Ϫ b j , system ͑1͒ is rewritten: 2 ˜11 uϩ 2 ˜12 Uϩdiv = s ϭ0 ͑5͒ 2 ˜22 Uϩ 2 ˜12 uϪhٌ pϭ0.
Using the second equation in ͑5͒, the displacement vector of the fluid phase U is expressed in terms of the pressure p in the pores and in terms of the displacement vector of the solid phase particle u:

Uϭ h ˜22 2 ٌ pϪ ˜12 ˜22 u. ͑6͒
This last equation will be used to express system ͑5͒ in terms of the (u,p) variables.

The solid phase equation in terms of (u,p) variables

Using Eq. ͑6͒, the first equation in ͑5͒ transforms into:

2 ˜uϩh ˜12 ˜22 ٌ pϩdiv = s ϭ0, ͑7͒
where the following effective density has been introduced:

˜ϭ ˜11 Ϫ ͑ ˜12 ͒ 2 ˜22 . ͑8͒
Equation ͑7͒ is still dependent on the fluid phase displacement vector U because of the dependency = s ϭ = s (u,U). Indeed, the stress-strain relations in the Biot theory are given by ͑Allard, 1993͒:

= s ͑ u,U ͒ϭA ˜div u1 =ϩ2N = s ϩQ ˜div U1 = ͑9͒ Ϫhp1 =ϭR ˜div U1 =ϩQ ˜div u1 =,
in which 1 = denotes the identity tensor, the elastic coefficients A ˜and N correspond to the Lame ´coefficients for elastic solids, Q ˜is a coupling coefficient between the dilatation and stress of the two phases, and R ˜may be interpreted as the bulk modulus of the air occupying a fraction h of a unit volume of aggregate. R ˜is related to the bulk modulus K ˜f of the air in the pores by the scale factor h: R ˜ϭhK ˜f . The elastic coefficients A ˜, Q ˜and R ˜are related to K b , the bulk modulus of the skeleton in vacuo,t oK s , the bulk modulus of the elastic solid from which the skeleton is made and to K ˜f , the bulk modulus of the air in the pores by ͑Allard, 1993͒:

A ˜ϭ ͑ 1Ϫh ͓͒1ϪhϪK b /K s ͔K s ϩh͑K s /K ˜f ͒)K b 1ϪhϪK b /K s ϩhK s /K ˜f Ϫ 2 3 N, ͑10͒ Q ˜ϭ ͓1ϪhϪK b /K s ͔hK s 1ϪhϪK b /K s ϩhK s /K ˜f , ͑11͒ R ˜ϭ h 2 K s 1ϪhϪK b /K s ϩhK s /K ˜f , ͑12͒
where N is the complex in vacuo shear modulus of the skeleton taking into account structural damping. Since K ˜f takes into account the thermal effects in the pore, it has a frequency dependent complex amplitude, and so do the elasticity tensors. In consequence, in the numerical implementation of the (u,U) formulation, the stiffness matrix associated to the frame is frequency dependent, which limits the computational efficiency of this formulation ͑the corresponding matrix should be reconstructed at each computational fre-quency͒.

To eliminate the dependency = s ϭ = s (u,U), the first and second equation in Eq. ͑9͒ are combined to obtain:

= s ͑ u,U ͒ϭ ͩ A ˜Ϫ Q ˜2 R ˜ͪ div u1 =ϩ2N = s Ϫh Q R ˜p1 =. ͑13͒
Next, introducing the tensor = s defined by

= s ͑ u ͒ϭ ͩ A ˜ϪQ ˜2 R ˜ͪ div u1 =ϩ2N = s , ͑14͒
the solid phase stress tensor writes:

= s ͑ u,U ͒ϭ = s ͑ u ͒Ϫh Q R ˜p1 =. ͑15͒
Since tensor = s depends only on the solid phase displace- ment vector, Eq. ͑15͒ may be used to eliminate the dependency = s ϭ = s (u,U) in Eq. ͑7͒. Substituting ͑15͒ into ͑7͒, the solid phase equation is obtained in terms of the (u, p) variables:

div = s ͑ u ͒ϩ ˜2 uϩ␥ ˜ٌ pϭ0, ͑16͒
where the following notation is used:

␥ ˜ϭh ͩ ˜12 ˜22 Ϫ Q R ˜ͪ . ͑17͒
To grasp the physical meaning of Eq. ͑16͒, note from Eqs. ͑10͒ to ͑12͒ that the expression ͑14͒ for = s is equivalent to:

= s ͑ u ͒ϭ͑ K b Ϫ 2 3 N ͒div u1 =ϩ2N = s , ͑18͒
which is the classical stress-strain relation for an elastic solid. Moreover, since K b and N denote the bulk modulus and the shear modulus of the skeleton in vacuo, respectively, tensor = s is the stress tensor of the material in vacuo.I n other words, = s represents the stress tensor prevailing in the poroelastic material when the fluid is totally drained off. In consequence, the first two terms of Eq. ͑16͒ represent the elastodynamic equation of the material in vacuo. One valuable feature is that = s is independent of the bulk modulus of the air and hence does not depend on frequency thus increasing the efficiency of the numerical implementation.

The fluid phase equation in terms of (u,p) variables

Next, to derive the fluid phase equation in terms of (u, p) variables, the divergence of Eq. ͑6͒ is taken and gives:

div Uϭ h 2 ˜22 ⌬ pϪ ˜12 ˜22 div u. ͑19͒
Combining this equation with the second equation in ͑9͒, the fluid phase equation is obtained in terms of the (u,p) variables:

⌬ pϩ ˜22 R ˜2 pϩ ˜22 h 2 ␥ ˜2 div uϭ0, ͑20͒
where ␥ ˜is defined in ͑17͒. Note that Eq. ͑20͒ is exactly the classical equivalent fluid equation for absorbing media with a source term ͑Allard, 1993; [START_REF] Panneton | A finite element formulation for the vibro-acoustic behaviour of double plate structures with cavity absorption[END_REF] The first two terms of this equation may be obtained directly from Biot's equations in the limit of a rigid skeleton. Those two terms represent the dynamic behavior of the material when its frame is supposed motionless. It is well known that Eq. ͑20͒ models correctly the behavior of rigid porous materials and/or situations in which the material is not mechanically excited or directly bonded to a vibrating structure and that for frequencies higher than the decoupling frequency ͑Panneton and Atalla, 1996͒.

Poroelasticity equations in terms of (u,p) variables

Grouping Eqs. ͑16͒ and ͑20͒, the Biot poroelasticity equations in terms of (u, p) variables are given by:

div = s ͑ u ͒ϩw 2 ˜uϩ␥ ˜ٌ pϭ0 ͑21͒ ⌬ pϩ 2 ˜22 R ˜pϪ 2 ˜22 h 2 ␥ ˜div uϭ0.
This system exhibits the classical form of a fluid-structure coupled equation. However, the coupling is of a volume nature since the poroelastic material is a superposition in space and time of the elastic and fluid phases. The first two terms of the structure equation represent the dynamic behavior of the material in vacuo while the first two terms of the fluid equation represent the dynamic behavior of the fluid when the frame is supposed motionless. The third terms in both equations couple the dynamics of the two phases.

The next section discusses the implementation of this formulation in a finite element code. It will be shown that this formulation is symmetric. One obvious advantage of the presented formulation is that four degrees of freedom are needed rather than 6 for the classical (u,U) formulation.

C. The weak integral formulation

Let ␦u and ␦p be two admissible variations of the displacement field u and the pressure field p, respectively. Then, using Galerkin's procedure ͑Reddy, 1991͒, the following symmetric (u,p) coupled weak integral formulation is obtained after few algebraic manipulations:

͵ ⍀ p = s ͑ u ͒: = s ͑ ␦u͒d⍀Ϫ 2 ͵ ⍀ p ˜u.␦ud⍀ Ϫ ͵ ⍀ p ␥ ˜ٌp.␦ud⍀Ϫ ͵ ץ⍀ p ͓ =s .n͔.␦udSϭ0, ͑22͒ ͵ ⍀ p ͫ h 2 2 ˜22 ٌ p.ٌ␦ pϪ h 2 R ˜p␦ p ͬ d⍀Ϫ ͵ ⍀ p ␥ ˜ٌ␦ p.ud⍀ ϩ ͵ ץ⍀ p ͫ ␥ ˜un Ϫ h 2 ˜22 2 ץp ץn ͬ ␦pdSϭ0, ᭙͑␦u,␦ p ͒,
where ⍀ p and ץ⍀ p denote the poroelastic domain and its boundary. Note from Eq. ͑22͒ the volume nature of the symmetric coupling between the two phases.

Equation ͑22͒ exhibits several advantages of the presented formulation over the (u,U) formulation. First, the stiffness matrix associated with the solid phase is frequency independent. Second, it accounts naturally for poroelasticporoelastic interfaces; no essential boundary conditions are needed. Indeed, let indices 1 and 2 denote the two interfacing poroelastic media, respectively. The continuity equations are given by ͑Panneton and Atalla, 1996͒:

u 1 ϭu 2 h 1 ͑ u 1n ϪU 1n ͒ϭh 2 ͑ u 2n ϪU 2n ͒ ͑23͒ p 1 ϭp 2 = 1 t .nϭ = 2 t .n,
where the first condition ensures the continuity of the solid phase displacement vector. The second equation ensures the continuity of the relative mass flux across the boundary. The two last equations ensure the continuity of the total normal stresses and the interstitial fluid pressure, respectively. Note that the total stress tensor = t is the total stress tensor given by:

= t ͑ u,U ͒ϭ = s ͑ u,U ͒ϩ = f ͑ u,U ͒ϭ = s ͑ u,U ͒Ϫhp1 =, ͑24͒
which may be written in terms of the in-vacuo stress tensor = S using Eq. ͑15͒:

= t ͑ u,U ͒ϭ = s ͑ u,U ͒Ϫh ͫ 1ϩ Q R ˜ͬ p1 =. ͑25͒
Using the following two relations, derived easily from Eqs. ͑6͒, ͑16͒ and ͑17͒:

͑ 1Ϫh ͒u.nϩhU.nϭ ͫ 1ϪhϪh Q R ˜ͬ u.n Ϫ ͫ ␥ ˜u.nϪ h 2 ˜22 2 ץp ץn ͬ , ͑26͒ h͑uϪU ͒.nϭ ͫ ␥ ˜u.nϪ h 2 ˜22 2 ץp ץn ͬ ϩh ͫ 1ϩ Q R ˜ͬ u.n, ͑27͒
the coupling conditions given by Eq. ͑23͒ may be written in terms of the (u,p) variables as follows

u 1 ϭu 2 ͫ ␥ ˜1u 1n Ϫ h 1 2 ˜22 2 ץp 1 ץn ͬ ϩh ͫ 1ϩ Q ˜1 R ˜1 ͬ u 1n ϭ ͫ ␥ ˜2 u 2n Ϫ h 2 2 ˜22 2 ץp 2 ץn ͬ ϩh ͫ 1ϩ Q ˜2 R ˜2 ͬ u 2n ͑28͒ p 1 ϭp 2 = 1 s .nϪh ͫ 1ϩ Q ˜1 R ˜1 ͬ u 1n ϭ = 2 s .nϪh ͫ 1ϩ Q ˜2 R ˜2 ͬ u 2n .
Using the fourth relation, the first boundary integrals in Eq. ͑22͒ related to the interface ץ⍀ p between the two coupled domains reads:

I 1 ϭϪ ͵ ץ⍀ p h 1ͩ 1ϩ Q ˜1 R ˜1 ͪ p 1 ␦u 1n dS ϩ ͵ ץ⍀ p h 2ͩ 1ϩ Q ˜2 R ˜2 ͪ p 2 ␦u 2n dS.

͑29͒

Similarly, using the second and fourth relations of Eq. ͑28͒, the second boundary integral of Eq. ͑22͒ reduces to

I 2 ϭϪ ͵ ץ⍀ p h 1ͩ 1ϩ Q ˜1 R ˜1 ͪ u 1n ␦p 1 dS ϩ ͵ ץ⍀ p h 2ͩ 1ϩ Q ˜2 R ˜2 ͪ u 2n ␦p 2 dS.

͑30͒

The signs in both equations account for the direction of the normal vector which is assumed to be directed from medium 1 toward medium 2. Equations ͑29͒ and ͑30͒ show that the two poroelastic media will be coupled through symmetrical coupling terms. This coupling occurs naturally between the solid phase normal displacement of one medium and the interstitial fluid pressure of the other medium. In addition to these surface coupling terms, the kinematic relation u 1 ϭu 2 and p 1 ϭp 2 will have to be explicitly applied on ץ⍀ p . For the majority of poroelastic media used in acoustics K b /K S Ӷ1 ͑Allard, 1993͒. Thus since ͓see Eqs. ͑11͒ and ͑12͔͒:

h Q R ˜ϭ1ϪhϪ K b K s , ͑31͒
we have:

h ͩ 1ϩ Q R ˜ͪ Х1. ͑32͒
In this case, both integrals I 1 and I 2 fall to zero and only the kinematic relations have to be applied. These conditions are taken into account through assembling. This is in contrast with the (u,U) formulation which involves the imposition of the second boundary condition in ͑23͒ through kinematic relations.

Another advantage of the presented (u, p) formulation over the (u,U) formulation lies in its coupling with air media. Indeed, denoting by p a the pressure in the air layer, the continuity equations are given by ͑Panneton and Atalla, 1996͒:

1 0 2 ץp a ץn ϭ͑1Ϫh ͒u n ϩhU n = t .nϭϪp a 1 =.n ͑33͒
= s .nϭϪ͑1Ϫh͒p a 1 =.n; = f .nϭϪhp a 1 =.n.

The first equation ensures the continuity of the normal volume velocity ͑modified Euler's equation due to the porosity͒, while the second and third equations ensure the continuity of the normal stresses. Using Eqs. ͑25͒ and ͑26͒, the coupling conditions in Eq. ͑33͒ are proved to be equivalent to:

ͫ 1ϪhϪh Q R ˜Ϫͬ u n Ϫ ͫ ␥ ˜un Ϫ h 2 ˜22 2 ץp ץn ͬ ϭ 1 0 2 ץp a ץn pϭp a ͑34͒ = s •nϭ ͫ h Q R ˜Ϫ͑1Ϫh ͒ ͬ pn.
Note that the first and third conditions in ͑34͒ lead to symmetric coupling terms of the form ͓see Eq. ͑22͔͒:

I 1 ϭ ͵ ץ⍀ p ͩ 1ϪhϪ Q R ˜ͪ p␦ u n dS ͑35͒
and

I 2 ϭ ͵ ץ⍀ p ͩ 1ϪhϪ Q R ˜ͪ u n ␦pdS, ͑36͒
where ץ⍀ p denotes the poroelastic-air interface. Once again for the majority of porous materials used in acoustics, the approximation given by Eq. ͑32͒ applies and these coupling terms vanish. Thus, only the continuity of the pressure at the poroelastic-air interface should be accounted for. In this case, contrary to the (u,U) formulation, no coupling matrix is needed. This represents another advantage and time saving feature of the presented formulation. Finally, it is worth mentioning that a detailed description of the boundary, support and loading conditions of the presented (u, p) formulation together with corresponding validation examples will be presented in a sequel paper ͑Debergue et al., 1998͒.

II. NUMERICAL IMPLEMENTATION

In the presented work, the weak formulation ͓Eq. ͑22͔͒ has been discretized using linear volume elements with four degrees of freedom per node: the three solid macroscopic displacements u and the pressure p in the interstitial fluid. Two finite volume elements have been implemented: a six node wedge element and an eight node brick element. Accordingly, within a finite element, it is assumed that the solid phase displacement vector and the pressure can be represented in matrix form following u e ϭ͓N s ͔͕u n ͖ e and p e ϭ͓N f ͔͕p n ͖ e ,

͑37͒

where ͓N s ͔ and ͓N f ͔ are the element's shape functions used to approximate the solid phase displacement vector and the interstitial pressure within element ''e.'' ͕u n ͖ e and ͕p n ͖ e are respectively the nodal displacement and the nodal pressure variables.

Since the solid phase governing equation is similar to a 3D elastodynamic equation, an analogy with solid elements is used. Similarly, since the fluid phase equation is similar to the Helmholtz's equation, an analogy with acoustic finite elements is used.

Substituting Eq. ͑37͒ into Eq. ͑22͒, the following matrices are obtained:

͵ ⍀ p = s ͑ u ͒: = s ͑ ␦u͒d⍀⇒͗␦u n ͓͘K͔͕u n ͖, ͑38͒ ͵ ⍀ p ˜u.␦ud⍀⇒͗␦u n ͓͘M ˜͔͕u n ͖, ͑39͒ ͵ ⍀ p ␥ ˜ٌ p.␦udS⇒͗␦u n ͓͘C ˜͔͕p n ͖, ͑40͒ ͵ ⍀ p h 2 ˜22 ͑ ٌ p.ٌ␦p ͒d⍀⇒͗␦ p n ͓͘H ˜͔͕p n ͖, ͑41͒ ͵ ⍀ p h 2 R ˜p.␦ pd⍀⇒͗␦p n ͓͘Q ˜͔͕p n ͖, ͑42͒ ͵ ⍀ p ␥ ˜u.ٌ␦ pdS⇒͗␦p n ͓͘C ˜͔T ͕u n ͖, ͑43͒
where ͕u n ͖ and ͕p n ͖ represent the solid phase and the fluid phase global nodal variables, respectively. ͓ M ˜͔ and ͓K͔ represent equivalent mass and stiffness matrices for the solid phase, respectively. ͓H ˜͔ and ͓Q ˜͔ represent equivalent kinetic and compression energy matrices for the fluid phase, respectively. Finally, ͓C ˜͔ is a volume coupling matrix between the solid phase displacement variables and the fluid phase pressure variable. Note that contrary to a classical fluid-structure coupling problem, the coupling within a poroelastic system is of a volume nature. The numerical evaluation of the different matrices is classical and will not be detailed in this paper. Note finally that the discretization of the interface integrals depends on the boundary conditions of the system ͑ϭ loading conditions͒.

Substituting Eqs. ͑38͒ to ͑43͒ into Eq. ͑22͒, the following coupled system is formed:

ͩ ͓K͔Ϫ 2 ͓ M ˜͔ Ϫ͓C ˜͔ Ϫ 2 ͓C ˜͔T ͓H ˜͔Ϫ 2 ͓Q ˜͔ ͪͭ u n p n ͮ ϭ ͭ F s F p ͮ , ͑44͒
where ͕ F p F s ͖ denotes the loading vector for the poroelastic medium ͑these vectors depends on the nature of the excitation; they are not detailed here͒. Solution of system ͑44͒ leads simultaneously to the nodal displacements of the solid phase and to the nodal pressures of the fluid phase. Once more, note that system ͑44͒ exhibits the classical form of a fluid structure coupled system. This system is smaller compared to the corresponding system in the classical (u,U) formulation since four degrees of freedom are used rather than six. Moreover, for a homogeneous poroelastic medium, all the frequency dependency in matrices ͑38͒ to ͑43͒ show up as multiplication coefficients that may be taken out of the integrals. Thus, the matrices may be calculated and assembled, once for all, for a given mesh. The spectra of the frequency dependent coefficient are accounted for during the formation of the global system ͑44͒. This is in contrast with the (u,U) formulation where the frequency dependence of the stiffness matrix is implicit in the stress-strain relations, thus necessitating recalculation of the matrix at each frequency. Moreover, it is worth mentioning that the presented formulation while improving significantly the performances of the numerical solutions of the Biot's poroelasticity equations still suffers from the necessity of using a spectral approach necessitated by the nonlinearity of the frequency dependence of the final matrix system. At this stage of the research, two potential solutions may be mentioned. The first is the use of a set of acceptable approximations that lead to a linearisation of the matrix system ͑see Panneton and Atalla, 1997b͒. The second is to use the selective modal reduction technique developed by [START_REF] Sgard | A model reduction technique the finite-element formulation of Biot's poroelasticity equations in acoustics[END_REF] Finally, recall from the previous section that coupling conditions are accounted for naturally in the presented formulation compared to the (u,U) formulation. The above mentioned features of the presented approach lead to important savings in setup and solution time as it will be demonstrated in the following numerical examples.

III. NUMERICAL EXAMPLES

In the following, the presented (u,p) formulation is validated by comparisons with the (u,U) formulation for both two-dimensional and three-dimensional problems. Note that detailed validation examples for the (u,U) formulation and its interface with elastic and poroelastic media have been presented elsewhere ͑Panneton and Atalla, 1996.

A. Surface impedance for a single poroelastic layer

As a first example, the surface impedance of laterally infinite poroelastic materials predicted by the presented (u, p) finite element Biot model of Eq. ͑21͒ is compared with the prediction of the (u,U) formulation. The configuration under study is depicted in Fig. 1͑a͒. A laterally infinite glass wool layer, described in Table I, is bonded onto a rigid impervious wall. A normal incidence plane wave of unit amplitude excites the absorbing material. To simulate the laterally infinite extent with the finite element model, only the axial macroscopic displacements are considered, i.e., the lateral displacements are set to zero. For both formulations, the poroelastic domain is meshed using an eight node linear brick element. A mesh of 7ϫ7ϫ5 elements was used. The normal incidence surface impedance is calculated using the axial nodal solid and fluid components, u n and U n respectively, at the input surface for the unit acoustic pressure excitation by the following equation: For the formulation, the normal displacement of the fluid phase U n is calculated using Eq. ͑6͒. This equation necessitates the evaluation of the nodal values of the pressure gradient. The classical procedure relying on evaluation of the pressure gradients at Gauss integration points followed by a least-square procedure in order to get nodal values is used.

Z n ϭ 1 j͑hU n ϩ͑1Ϫh͒u n ͒ . ͑45͒
Note that the results of the (u,U) finite element Biot model have been validated by comparison with an exact analytical calculation ͑Allard, 1993͒. Figure 2 presents the real and imaginary parts of the surface impedance predicted by both models. An excellent agreement is observed. For this simple problem, the (u,p) formulation took approximately one-fifth of the time needed for the (u,U) formulation.

B. Surface impedance for a multilayer poroelastic material

To show the efficiency of the approach, the surface impedance of a laterally infinite multilayer poroelastic material is calculated. The configuration under study is depicted in Fig. 1͑b͒. A laterally infinite multilayer poroelastic material is bonded onto a rigid impervious wall. From the front face to the rear face, the multilayer consists of a blanket, a screen, foam A, and foam B. The material properties are listed in Table I. A normal incidence plane wave of unit amplitude excites the absorbing material. To simulate the laterally infi-nite extent with the finite element model, only the axial macroscopic displacements are considered, i.e., the lateral displacements are set to zero. Identical meshes using eight node linear brick elements are used for both formulations. The cross section of the material was meshed using 7ϫ7 elements. The blanket, the screen, foam A and foam B were meshed using five elements. Note that this problem has been validated elsewhere by comparison with an analytical calculation ͑Allard, 1993͒. Figure 3 compares the real and imaginary parts of the surface impedance predicted by both models. A perfect agreement is observed. For this problem, the (u, p) formulation took approximately one-fourth the time needed for the (u,U) formulation.

C. A rigid cavity with surface absorption

Next, a coupled acoustic-poroelastic problem is considered. The geometry of the problem is depicted in Fig. 4. A wall of a rigid cavity ͑dimensions: 0.35 m ϫ 0.22 m ϫ 0.01 m͒ is treated with 10 cm of the fiberglass of Table I. The fiberglass is totally bonded to the rigid wall. The system is excited with a point source positioned at a corner of the rigid cavity. The quadratic pressure in the cavity is calculated using the (u,U) and the (u, p) formulations. Identical meshes are used for both formulations. The cavity is discretized using six-node linear pentahedron acoustic elements with one FIG. 2. Comparison of the surface impedance calculated with the (u,U) and the (u, p) Biot finite element models for the laterally infinite glass wool layer of Table I. degree-of-freedom per node: the acoustic pressure. A 7ϫ7 ϫ7 mesh is used for the cavity. For the poroelastic domain, a7 ϫ 7 ϫ 3 mesh is used. Note that this problem has been validated elsewhere by comparison with an impedance model for the poroelastic material ͑Atalla and Panneton, 1996͒. The results of the calculation are given in Fig. 5. Once again, excellent agreement is found between the two approaches. The (u,p) formulation took approximately onetwelfth of the time needed with the (u,U) formulation. This is mostly due to the reduction of the number of degrees of freedom for the problem and the more natural way of accounting for the coupling between the poroelastic and the acoustic domain.

IV. CONCLUSION

The presented mixed (u, p) formulation for poroelastic materials has several features: ͑1͒ It is exact, in the sense that no new assumptions introduced beyond those governing Biot's poroelastic equations. ͑2͒ It leads to a classical coupled fluid-structure problem involving the dynamic equations of the skeleton in vacuo and the equivalent fluid in the rigid skeleton limit.

͑3͒ It handles naturally coupling conditions with acoustic and other pyroelastic media. ͑4͒ It involves four degrees of freedom per node element compared to six for the classical (u,U) formulation; important savings in computer storage and solution time are achieved. ͑5͒ The physical nature of its coupled equations and its variables allows for efficient numerical implementation; important savings in setup and solution time are achieved. ͑6͒ For the presented numerical simulations, the (u,p) was shown to be 5 to 12 times faster than the (u,U) formulation with identical accuracy.

FIG. 1 .

 1 FIG. 1. Geometry of the surface impedance problem. ͑a͒ Single poroelastic layer; and ͑b͒ multilayered poroelastic material.

FIG. 3 .

 3 FIG.3. Comparison of the surface impedance between the (u,U) and the (u, p) Biot finite element models for the laterally infinite multilayer.

FIG. 5 .

 5 FIG. 5. Comparison of the quadratic velocity at the poroelastic surface calculated with the (u,U) and the (u, p) Biot finite element models for the cavity problem of Fig. 4.

TABLE I .

 I Physical properties and dimensions of the poroelastic materials.

			l			N		⌳	⌳Ј	Thickness
		k s	(kg/m 3 )	(Ns/m 4 )	h	͑kPa͒		͑m͒	͑m͒	͑cm͒
	Glass wool 1.06	130	40 000	0.94 2200 ͑1ϩj0.1͒ 0	0.56ϫ10 Ϫ4 1.10ϫ10 Ϫ4	10
	Blanket	1.18	41	34 000	0.98 110͑1ϩj0.015͒ 0.3 0.60ϫ10 Ϫ4 0.87ϫ10 Ϫ4	0.4
	Screen	2.56	125	320ϫ10 4 0.80 1000 ͑1ϩj0.1͒ 0.3 0.06ϫ10 Ϫ4 0.24ϫ10 Ϫ4	0.08
	Foam A	2.52	31	87 000	0.97 55 ͑1ϩj0.055͒ 0.3 0.37ϫ10 Ϫ4 1.19ϫ10 Ϫ4	0.5
	Foam B	1.98	16	65 000	0.99 18 ͑1ϩj0.1͒	0.3 0.37ϫ10 Ϫ4 1.21ϫ10 Ϫ4	1.6
	Fiberglass 1.4	30	25 000	0.95 21 ͑1ϩj0.05͒	0.0 0.93ϫ10 Ϫ4 0.93ϫ10 Ϫ4	7.62

ACKNOWLEDGMENTS

The authors thank Professors Jean-Franc ¸ois Allard and M. A. Hamdi for various stimulating discussions during this work. Thanks are also expressed to Bombardier Inc. Canadair, N.S.E.R.C. and F.C.A.R. for their financial support.