T A Karatsori
email: theano.karatsori@grenoble-inp.fr

K Bennamane

G Ghibaudo
email: ghibaudo@minatec.inpg.fr

Novel On-Resistance based Methodology for MOSFET Electrical Characterization

Keywords: MOSFET, parameter extraction, On-resistance

A new methodology for MOSFET characterization making use of the on-resistance characteristics Ron(Vg,Vd)=Vd/Id(Vg,Vd) and associated derivatives dRon/dVg and dRon/dVd is proposed. This approach enables to eliminate the influence of source-drain series resistance Rsd not only in linear region but also in non-linear region of MOSFET operation. Therefore, it allows for intrinsic MOSFET parameter extraction free from source and drain series resistance.

Introduction

A precise extraction of MOSFET electrical parameters is needed for the proper engineering and optimization of CMOS technologies. Numerous methods have been proposed during the past decades for the MOSFET parameter extraction, mainly based on the Id(Vg,Vd) transfer or output drain current characteristics [START_REF] Krutsick | An improved model of MOSFET modeling and parameter extraction[END_REF][START_REF] Ghibaudo | New method for the extraction of MOSFET parameters[END_REF][START_REF] Taur | A new 'shift and ratio' method for MOSFET channel-length extraction[END_REF][START_REF] Katto | Device parameter extraction in the linear region of MOSFET's[END_REF][START_REF] Ortiz-Conde | A review of recent MOSFET threshold voltage extraction methods[END_REF][START_REF] Tanaka | Novel extraction method for size-dependent mobility based on BSIM3-like compact model[END_REF][START_REF] Fleury | A New Technique to Extract the Source/Drain Series Resistance of MOSFETs[END_REF][START_REF] Flandre | gm/Id method for threshold voltage extraction applicable in advanced MOSFETs with nonlinear behavior above threshold[END_REF][START_REF] Ferdousi | Improved MOSFET characterization technique for single channel length, scaled transistors[END_REF][START_REF] Karatsori | Full gate voltage range Lambert-function based methodology for FDSOI MOSFET parameter extraction[END_REF][START_REF] Henry | New Y-function based MOSFET parameter extraction method from weak to strong inversion range[END_REF][START_REF] Cretu | New ratio method for effective channel length and threshold voltage extraction in MOS transistors[END_REF]. In this context, the Y-function technique, where Yg=Id/√gm, gm being the transconductance, has proven very efficient for MOSFET parameter extraction in linear region, since it eliminates the impact of source-drain series resistance Rsd [START_REF] Ghibaudo | New method for the extraction of MOSFET parameters[END_REF][START_REF] Fleury | A New Technique to Extract the Source/Drain Series Resistance of MOSFETs[END_REF][START_REF] Henry | New Y-function based MOSFET parameter extraction method from weak to strong inversion range[END_REF][START_REF] Cretu | New ratio method for effective channel length and threshold voltage extraction in MOS transistors[END_REF], at least when Rsd is constant with gate voltage.

In this paper, we propose a new methodology for MOSFET characterization making use of the on-resistance characteristics Ron(Vg,Vd)=Vd/Id(Vg,Vd) and associated derivatives dRon/dVg and dRon/dVd. Indeed, this approach enables to eliminate the influence of Rsd not only in linear region but also in non-linear region of MOSFET operation.

Ron Methodology principles

Considering Ohm's law, the power dissipated in a MOSFET, P=IdVd, can be expressed as [START_REF] Xu | Power transfer-length method for full biasing contact resistance evaluation of organic fieldeffect transistors[END_REF],

2 0 2 2 d on d d d s d d I R I R I R V I P + + = = (1)
where Ron 0 is the intrinsic on-resistance, Rs and Rd are the source and drain series resistances, respectively. Therefore, it can be shown from (1) that the extrinsic on-resistance Ron is simply obtained by adding Rs and Rd to the intrinsic on-resistance, Ron 0 as [START_REF] Xu | Power transfer-length method for full biasing contact resistance evaluation of organic fieldeffect transistors[END_REF],

0 on d s on R R R R + + = (2)
As a result, Eq. ( 2) indicates that, as it was the case for MOSFET linear operation region [START_REF] Krutsick | An improved model of MOSFET modeling and parameter extraction[END_REF][START_REF] Ghibaudo | New method for the extraction of MOSFET parameters[END_REF][START_REF] Taur | A new 'shift and ratio' method for MOSFET channel-length extraction[END_REF][START_REF] Katto | Device parameter extraction in the linear region of MOSFET's[END_REF][START_REF] Ortiz-Conde | A review of recent MOSFET threshold voltage extraction methods[END_REF][START_REF] Tanaka | Novel extraction method for size-dependent mobility based on BSIM3-like compact model[END_REF][START_REF] Fleury | A New Technique to Extract the Source/Drain Series Resistance of MOSFETs[END_REF][START_REF] Flandre | gm/Id method for threshold voltage extraction applicable in advanced MOSFETs with nonlinear behavior above threshold[END_REF][START_REF] Ferdousi | Improved MOSFET characterization technique for single channel length, scaled transistors[END_REF][START_REF] Karatsori | Full gate voltage range Lambert-function based methodology for FDSOI MOSFET parameter extraction[END_REF][START_REF] Henry | New Y-function based MOSFET parameter extraction method from weak to strong inversion range[END_REF][START_REF] Cretu | New ratio method for effective channel length and threshold voltage extraction in MOS transistors[END_REF],

the use of on-resistance Ron allows generalizing the application of Ohm's law to the MOSFET non linear operation region. In turn, assuming as usual constant Rs and Rd, it implies that one could suppress the influence of the source and drain series resistance by dealing with its derivatives with respect to Vg, as it was previously done in ohmic linear regime [START_REF] Krutsick | An improved model of MOSFET modeling and parameter extraction[END_REF][START_REF] Ghibaudo | New method for the extraction of MOSFET parameters[END_REF][START_REF] Taur | A new 'shift and ratio' method for MOSFET channel-length extraction[END_REF][START_REF] Katto | Device parameter extraction in the linear region of MOSFET's[END_REF][START_REF] Ortiz-Conde | A review of recent MOSFET threshold voltage extraction methods[END_REF][START_REF] Tanaka | Novel extraction method for size-dependent mobility based on BSIM3-like compact model[END_REF][START_REF] Fleury | A New Technique to Extract the Source/Drain Series Resistance of MOSFETs[END_REF][START_REF] Flandre | gm/Id method for threshold voltage extraction applicable in advanced MOSFETs with nonlinear behavior above threshold[END_REF][START_REF] Ferdousi | Improved MOSFET characterization technique for single channel length, scaled transistors[END_REF][START_REF] Karatsori | Full gate voltage range Lambert-function based methodology for FDSOI MOSFET parameter extraction[END_REF][START_REF] Henry | New Y-function based MOSFET parameter extraction method from weak to strong inversion range[END_REF][START_REF] Cretu | New ratio method for effective channel length and threshold voltage extraction in MOS transistors[END_REF], but also with respect to Vd for any Vg and Vd bias condition.

It should be mentioned that this approach can be regarded as an attempt to generalize the concept of Y-function, originally derived for linear regime, to the non-linear operation region. Indeed, generalized Y-functions can be defined from the derivative of Ron with respect to Vg and Vd as:

g on d g dV dR V Y / - = (3)
d on d d dV dR V Y / = . ( 4
)
It is worth noting that Eq. ( 3) naturally reduces to the standard definition for the Y-function, Yg=Id/√gm=√β.√Vd.(Vg-Vt) for small Vd at strong inversion.

A. Long channel MOSFET at strong inversion

d gt d L d V V V I - = β with 0 . . µ C L W ox = β (5)
where Vgt=Vg-Vt is the gate voltage overdrive, Cox the gate oxide capacitance per unit area, L the channel length, W the channel width and µ0 the low-field mobility. It is easy to show using the definition of the total resistance (Ron= Vd/Id) that the corresponding expressions for the on-resistance Ron and associated derivatives dRon/dVg and dRon/dVd are given by:

) 2 / .( 1 d gt L on V V R - = β (6) 2 2 . ) 2 / .( 1 L on d gt g L on R V V dV dR β β - = - - = (7) 4 2 2 . 2 ) 2 / .( 2 1 L on d gt d L on R V V dV dR β β = - = (8)
Note the proportionality of dRon/dVg and dRon/dVd with Ron 2 .

B. Short channel MOSFET at strong inversion

For short channel devices, using the classical MOSFET model accounting for saturation velocity effect, we have [START_REF] Taur | Fundamentals of modern VLSI[END_REF]:

L v V µ I I sat d L d d 0 1+ = ( 9
)
where vsat is the carrier saturation velocity. The associated expressions for the total on-resistance Ron (Ron= Vd/Id) and associated derivatives dRon/dVg and dRon/dVd can also be derived and now read,

) 1 .( 0 L v V µ R R sat d L on on + = ( 10
) ) 1 .( . ) 1 .( 0 2 0 L v V µ R L v V µ dV dR dV dR sat d L on sat d g L on g on + - = + = β (11) ) . . 2 1 .( . 2 0 0 2 L v R µ L v V µ R dV dR sat L on sat d L on d on β β + + = (12)
In this case, it should be noted that, unlike long channel devices, (dRon/dVd)/(-dRon/dVg) is no longer constant with Vd.

Results and discussion

In order to validate the above analysis, the on-resistance Ron(Vg,Vd) and the associated derivative characteristics have been reconstructed from Id(Vg,Vd) measurements performed on nMOSFETs from a 14nm FDSOI technology (as in [START_REF] Karatsori | All operation region characterization and modeling of drain and gate current mismatch in 14nm Fully Depleted SOI MOSFETs[END_REF] with gate width W=1µm and gate lengths from L=300nm down to L=60nm.

The impact of Rs and Rd was tested by adding external resistors (here Rext=400Ω) in series with source and drain electrodes. For example, Fig. 1 illustrates how the Id(Vg) characteristics measured in linear operation (Vd=20mV) are significantly modified by the insertion of an external series resistance, especially as the gate length L is reduced. Figure 2 confirms the fact that the onresistance Ron is well shifted by Rext from linear to non linear triode region as Vd varies from 0 to 0.5V. Moreover, Fig. 3 clearly indicates that dRon/dVg is independent of Rext for long and short channel devices. However, Fig. 4 reveals that, although dRon/dVd is independent of Rext for long channel devices, it is slightly varying with Rext for short devices. This fact could possibly be related to the assumptions we have considered. More precisely, in our theory low-field mobility, µ0, is considered independent of the gate voltage, Vg, and also the source and drain series resistances, Rs and Rd, do not depend on Vg.

As a result, Fig. 5 shows that the generalized Y-function characteristics Yg(Vg) are well immune to Rext as expected from Eqs (2)-( 3) whatever the channel length. This feature demonstrates that the Yg(Vg) function does provide information on the intrinsic MOSFET operation -free from series resistance effect-not only in ohmic region but also in non linear operation region i.e. for large drain voltage. Instead, for short channel, the generalized Y-function characteristics Yd(Vd) of Fig. 6 slightly varies with Rext, especially in triode region. Therefore, the Yd(Vd) characteristics are not perfectly immune to Rext, contrary to Yg(Vg) characteristics.

Figure 7 shows that the ratio (dRon/dVd)/(dRon/dVg) is lying around 0.5 from linear to non-linear region as predicted from Eqs ( 7)-( 8) for long channel devices, emphasizing the consistency of Ron model of Eqs ( 6)-( 8).

As suggested from Eq. ( 7) and shown in Fig. 8, the on-resistance Ron varies linearly with (dRon/dVg) 1/2 for long channel devices. This property could be used to extract the intrinsic Rsd of a single device after extrapolation of the straight line on y axis to zero abscissa as illustrated in Fig. 8.

It should be noted that the extracted value of Rsd (≈650Ω) by this method falls in agreement with the one extracted with standard techniques using several gate lengths as in [START_REF] Karatsori | All operation region characterization and modeling of drain and gate current mismatch in 14nm Fully Depleted SOI MOSFETs[END_REF]. For short channel devices, this property seems also to be applied as shown in Fig. 8, despite the difference in Eq. [START_REF] Henry | New Y-function based MOSFET parameter extraction method from weak to strong inversion range[END_REF].

However, in this case, the extracted Rsd value differs from the one extracted by standard techniques [START_REF] Karatsori | All operation region characterization and modeling of drain and gate current mismatch in 14nm Fully Depleted SOI MOSFETs[END_REF] and long channel case. Therefore, because of the extrapolation inaccuracy and limitation for short channel devices, we do not recommend this method for Rsd extraction on single short channel devices due to larger extrapolation error.

Figures 9 and 10 show the variations of (dRon/dVd)/(dRon/dVg) with Vd for various channel lengths and how they can be satisfactorily modeled by saturation velocity effect using Eqs ( 11) and [START_REF] Cretu | New ratio method for effective channel length and threshold voltage extraction in MOS transistors[END_REF]. However, the deviations at low drain voltage should be noted and are difficult to be interpreted by the simple model of Eqs ( 11) and [START_REF] Cretu | New ratio method for effective channel length and threshold voltage extraction in MOS transistors[END_REF].

In order to confirm the previous analysis and gain physical insight, we have performed calibrated modelling of the device characterstics using the Lambert function based compact model of [START_REF] Karatsori | Full gate voltage range Lambert-function based methodology for FDSOI MOSFET parameter extraction[END_REF] after inclusion of saturation velocity and series resistance effects.

Conclusions

A new methodology for MOSFET characterization using the on-resistance Ron(Vg,Vd) and associated derivatives dRon/dVg and dRon/dVd characteristics has been presented. This approach enables to eliminate the influence of Rsd not only in linear region but also in non-linear region of MOSFET operation. However, the dRon/dVd(Vd) characteristics are not perfectly independent of external series resistance -especially for short devices, as confirmed by simulation. Instead, the dRon/dVg(Vg) and, by turn, Yg(Vg) characteristics are found perfectly immune to source-drain series resistance, generalizing the use of Y-function to non linear operation region. This methodology should be very useful for MOSFET intrinsic parameter extraction over full gate and drain voltage operation range.

Figure captions

For

  long channel devices, based on the classical MOSFET model valid in strong inversion,

  Typical calibration curves obtained with this compact model for the transfer characteristics gm(Vg) and output chacracteristics Id(Vd) are shown in Figs 11 and 12 for linear and non linear regions for short channel device (L=60nm). As can be seen from Fig. 13, the modeled dRon/dVg(Vg) and dRon/dVd(Vd) characteristics allow to reproduce satisfactorily the experimental data of Figs 3 and 4 and, in particular, some dependence of dRon/dVd(Vd) curves with Rext. These simulations confirm that dRon/dVg(Vg) and, by turn, Yg(Vg) are perfectly immune to Rext, whereas dRon/dVd(Vd) and Yd(Vd) varies somehow with Rext, expecially in short channel devices.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Experimental Id(Vg) characteristics with Rext=400Ω (dashed line) and without Rext (solid

Fig. 3 .

 3 Fig. 3. dRon/dVg(Vd) with Rext=400Ω (dashed line) and without Rext (solid line) for Vg=1V

Fig. 4 .

 4 Fig. 4. dRon/dVd(Vd) with Rext=400Ω (dashed line) and without Rext (solid line) for Vg=1V

Fig. 5 .

 5 Fig. 5. Yg(Vg) with Rext=400Ω (dashed line) and without Rext (solid line) for various Vd and Vg=1V

Fig. 6 .

 6 Fig. 6. Yd vs √Vd with Rext=400Ω (dashed line) and without Rext (solid line) for Vg=1V (L=60nm

Fig. 7 .

 7 Fig. 7. (dRon/dVd)/(dRon/dVg) vs Vd with Rext=400Ω (dashed line) and without Rext (solid line) for

Fig. 8 .

 8 Fig.8. Ron vs (dRon/dVg) 1/2 as Vd is varying from 0 to 0.5V with Rext=400Ω (dashed line) and

Fig. 9 .

 9 Fig. 9. Experimental (dRon/dVd)/(dRon/dVg) vs Vd without Rext for various L (Vg=1V).

Fig. 10 .

 10 Fig. 10. Modeled (dRon/dVd)/(dRon/dVg) vs Vd for various L (Vg=1V, vsat=2.10 7 cm/s, µ0=140cm 2 /Vs, Vt=0.4V).

Fig. 11 .

 11 Fig. 11. Experimental (left) and modeled (right) gm(Vg) characteristics with Rext=400Ω and without

Fig. 12 .

 12 Fig. 12. Experimental (left) and modeled (right) Id(Vd) characteristics with Rext=400Ω and without

Fig. 13 .Fig. 1 .Fig. 3 .

 1313 Fig. 13. Modeled dRon/dVg and dRon/dVd vs Vd characteristics with Rext=400Ω and without Rext for

Fig. 4 .

 4 Fig. 4. dRon/dVd(Vd) with Rext=400Ω (dashed line) and without Rext (solid line) for Vg=1V

Fig. 5 .(

 5 Fig. 5. Yg(Vg) with Rext=400Ω (dashed line) and without Rext (solid line) for various Vd and Vg=1V

Acknowledgments

This work has been partially supported by 3D-MUSE European project under Grant agreement n° 780548.

Ron (Ω)


T Cazimajou

M Mouis

M Legallais

T T T Nguyen

C Ternon

Bassem Salem

Gérard Ghibaudo

Analysis of the Role of Inter-Nanowire Junctions on Current Percolation Effects in Silicon Nanonet Field-Effect Transistors

Keywords: Random network, silicon nanowire, electronic transport, compact model, percolation, energy barrier, Monte-

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.

substrate, which make them perfectly suitable to the 3D integration of sensors, with direct vertical wiring to the CMOS conditioning and readout electrics, or to an integration on flexible substrates. [START_REF] Gruner | Carbon nanonets spark new electronics[END_REF] The first NN structures have been made with other materials than silicon, for example carbon nanotubes (CNT) [START_REF] Snow | Random networks of carbon nanotubes as an electronic material[END_REF][START_REF] Kumar | Percolating Conduction in Finite Nanotube Networks[END_REF] and ZnO nanowires. [START_REF] Unalan | Zinc Oxide Nanowire Networks for Macroelectronic Devices[END_REF] These structures, especially CNT NN-FETs have been deeply studied with static, RF measurements [START_REF] Oueriemi | Wideband non-linear characteristics of random multi-walled carbon nanotubes network[END_REF] and low frequency noise measurements. [START_REF] Joo | Static and low frequency noise characterization of N-type random network of carbon nanotubes thin film transistors[END_REF] However, some characteristics of the CNT-based nanonets make them quite different from the silicon based nanonets. One important feature of the former comes from the respective roles of metallic and semiconducting CNTs on transport. [START_REF] Topinka | Charge Transport in Interpenetrating Networks of Semiconducting and Metallic Carbon Nanotubes[END_REF] Moreover, the contacts between CNTs have rather been modelled with Variable Range Hopping. [START_REF] Yanagi | Transport mechanisms in metallic and semiconducting single-wall carbon nanotube networks[END_REF] Conduction through a ZnO nanonet has been recently studied for chemiresistive detection applications with a focus on the analysis of the Id-Vd characteristic under the influence of barrier height modification at nanowire-nanowire junctions. [START_REF] Caicedo | Detection mechanism in highly sensitive ZnO nanowires network gas sensors[END_REF] Until recently, the main technological roadblock to the fabrication of Si-NN FETs has been the presence of nonconductive junctions between Si-NWs, due to the presence of native silicon dioxide layer at their interface. [START_REF] Heo | Large-scale assembly of silicon nanowire network-based devices using conventional microfabrication facilities[END_REF] For this reason, Si-NN was known to be insulating for channel lengths larger than NW length. [START_REF] Heo | Large-scale assembly of silicon nanowire network-based devices using conventional microfabrication facilities[END_REF] As a result, Si-NN FETs have been mainly used in devices featuring a maximal channel length close to the length of a single NW. For the same reason, there has not been yet deep studies of the transport through a random network of silicon nanowires, nor measurement of the electrical characteristics of the junctions. The detrimental effect of the native oxide at the junctions has been lately solved with a sintering at low temperature, compatible with CMOS integration [START_REF] Ternon | Low Temperature Processing to Form Oxidation Insensitive Electrical Contact at Silicon Nanowire/Nanowire Junctions[END_REF], leading to the first integration of Si NN. in a back gated transistor architecture. [START_REF] Legallais | Toward the integration of Si nanonets into FETs for biosensing applications[END_REF] This Si-NN FET has been extensively characterized through static measurements. It has been shown that a combination of percolation and dispersion effects was required to understand the dependence of static characteristics with channel length and with the density of nanowires in the nanonet. Especially, it has been explained how nanowires threshold voltage dispersion explains the variation of subthreshold slope and threshold voltage of the NN FET with geometry and NW density. [START_REF] Cazimajou | Electrical characteristics of silicon percolating nanonet-based Field Effect Transistors in the presence of dispersion[END_REF] Si-NN FET has been then modelled as a NW network with a random dispersion of threshold voltage. [START_REF] Cazimajou | Finite element simulation of 2D percolating silicon-nanonet field-effect transistor[END_REF] More recently, Si-NN FETs have been characterized through low frequency noise measurements. [START_REF] Cazimajou | Low-frequency Noise Characterization of Si Nanonet Field Effect Transistors[END_REF] In this paper, Si-NN FET were characterized through static electrical measurements at different temperatures. Experimentally, we found that ON-current and mobility were varying with temperature in significantly different ways, depending on channel length and nanowire density. These trends cannot be explained without assuming that there is a potential barrier at the interface of at least some NW/NW junctions. In this paper, the understanding of this behaviour is assessed with the help of Monte-Carlo simulations of Si-NN where internanowire junction resistance is accounted for. Our aim here is not to provide a fit to experimental results but to identify the origin of the observed trends.

FABRICATION

The process of Si-NN fabrication and their integration into back gate field effect transistor has been described elsewhere. [START_REF] Legallais | An innovative large scale integration of silicon nanowire-based field effect transistors[END_REF][START_REF] Legallais | Material engineering of percolating silicon nanowire networks for reliable and efficient electronic devices[END_REF][START_REF] Nguyen | Monolithic fabrication of nano-to-millimeter scale integrated transistors based on transparent and flexible silicon nanonets[END_REF] The NNs used in this study (Fig. 1-(a)) were self-assembled from Si-NWs grown by Vapor Liquid Solid chemical vapor deposition. NW length and diameter were respectively 7 µm ± 3 µm and 40 nm ± 7 nm. NNs can be fabricated with different NW densities. [START_REF] Serre | Fabrication of silicon nanowire networks for biological sensing[END_REF] For this paper, several NNs were fabricated with 4 different values of nanowire density dWires (from 0.28 NWs.µm -2 to 0.65 NWs.µm -2 ). They were then transferred on a Si / Si3N4 wafer. This material was chosen instead of the more commonly used SiO2 for technological reason. Indeed, junction sintering is prepared by the removal of the native oxide surrounding the nanowires by an HF-based Buffered Oxide Etching (BOE) step which would remove the back-gate oxide and short-circuit gate and channel. Thermal annealing at 400 °C was performed to improve electronic transport through the network, by means of NW-NW junctions sintering. [START_REF] Ternon | Low Temperature Processing to Form Oxidation Insensitive Electrical Contact at Silicon Nanowire/Nanowire Junctions[END_REF][START_REF] Legallais | An innovative large scale integration of silicon nanowire-based field effect transistors[END_REF] Then, Si-NNs were passivated with a 8 nm thick Al2O3 layer using Atomic Layer Deposition [START_REF] Legallais | Material engineering of percolating silicon nanowire networks for reliable and efficient electronic devices[END_REF] and finally contacted by Ni/Au metal pads using a lift-off technique associated with local removal of Si3N4 by HF treatment in contact pad openings. Field-effect control was obtained by using the Si / (200nm) Si3N4 substrate as a back gate. Finally, a last thermal annealing at 400 °C was realized to promote silicide formation at the interface between Si NWs and metal contact pads. [START_REF] Legallais | An innovative large scale integration of silicon nanowire-based field effect transistors[END_REF] Several device geometries were available with Source-Drain distance Lmask ranging from 5 µm to

EXPERIMENTAL PROTOCOL AND RESULTS

NN FETs were characterized using a Süsstec cryogenic probe station, where temperature could be varied between 80 K and 400 K. An HP 4155A semiconductor parameter analyzer was used to measure drain current (Id) as a function of the voltage applied on the backgate (Vg) with the source contact linked to the ground and a constant voltage (Vd = -4 V or Vd = -9 V) applied to drain. ). An increase of drain current with density is visible at constant temperature and constant source-drain distance, which can be explained by the increase of the number of paths bridging source and drain. This increase has been explained in [START_REF] Cazimajou | Electrical characteristics of silicon percolating nanonet-based Field Effect Transistors in the presence of dispersion[END_REF] for Si-NN FET and is a well-known result for metallic nanonets [START_REF] Mutiso | Integrating simulations and experiments to predict sheet resistance and optical transmittance in nanowire films for transparent conductors[END_REF] and CNT NN FETs. [START_REF] Kumar | Percolating Conduction in Finite Nanotube Networks[END_REF] This variation is nonlinear and can be modelled with a power law using percolation theory for a 2D stick network. [START_REF] Stauffer | Introduction to Percolation Theory[END_REF] However the detailed analysis of this dependence would require more statistical measurements and is outside the scope of this paper. Short devices (Lmask = 5 µm) have a S-D distance smaller than NW length (7 µm) and the channel consisted essentially in NWs linking directly S and D, with negligible role of NW-NW junctions. In contrast, with a S-D distance (Lmask = 50 µm or Lmask = 1000 µm) much larger than the average NW length (7 µm), the percolating conduction paths between source and drain involved several NWs, connected by sintered NW-NW junctions, in series. Comparison of Fig. 2-(a These curves were used to extract the main electrical parameters such as low field mobility µ0, subthreshold slope ideality factor n and threshold voltage Vth. Extraction was performed in most cases by fitting the experimental curves with a Lambert W function. [START_REF] Akkez | New parameter extraction method based on split C-V measurements in FDSOI MOSFETs[END_REF][START_REF] Karatsori | Full gate voltage range Lambert-function based methodology for FDSOI MOSFET parameter extraction[END_REF][START_REF] Cazimajou | Electrical characterization of percolating silicon nanonet FETs for sensing applications[END_REF] In the few cases when the threshold voltage was out of the measurement window, we used instead a Y-function based model for µ0 and Vth extraction. This was the case only for devices with 0.65 NWs.µm -2 and Lmask = 5 µm. The value of µ0 was obtained using Wmask and Lmask as channel width and length, respectively. This provided an apparent low field mobility value suitable for compact modelling. Because physical width is smaller than Wmask, the extracted µ0 is smaller than the physical mobility. This does not prevent variations to be analysed. Measurements were carried out on up to 5 devices of identical geometry for each value of Lmask and dWires.

In the following, we will in general show the results obtained for these identical devices, in order to keep the information about variability.

DISCUSSION

A. Extraction of activation energy

In short channel devices (Lmask = 5 µm), conduction paths are mainly made of nanowires which connect directly source and drain. The temperature variation of mobility can thus be considered as indicative of transport mechanisms in the nanowires. In contrast, with channel lengths much longer than mean nanowire length (which is the case with

Gate Voltage V g (V)

Increasing of T Lmask = 50 µm), current flows through a network of nanowires connected by inter-nanowire junctions. It is important to understand the role of these junctions in the percolating regime. This may allow getting some insight on their quality by proper analysis of experimental trends. Fig. 3 compares the temperature dependence of µ0 obtained for Lmask = 5 µm and Lmask = 50 µm, for a NW density dWires = 0.28 NWs.µm -2 . For Lmask = 5 µm, µ0 remained mostly constant or decreased slightly at low temperature, showing that the freezing of phonon scattering is screened and that scattering by defects is the prevailing scattering mechanism. [START_REF] Ando | Electronic properties of two-dimensional systems[END_REF] In contrast, it was observed that µ0 decreased as temperature decreased for long channel devices (Lmask = 50 µm) operating in the percolating regime, with inter-nanowire junctions present along the conduction paths. This trend could in principle be consistent with the presence of energy barrier heights at the junctions and it is the aim of this paper to understand to which extent this is the case. With a higher density of nanowires (dWires = 0.65 NWs.µm -2 ) and short channel devices (L = 5 µm), similar temperature dependence of µ0 was observed as µ0 remained almost independent of temperature (plots obtained for Lmask = 5 µm in Fig. 4). However this remained the case even with channels as long as 50 µm (plots for Lmask = 50 µm in Fig. 4). Thermal activation was recovered only for much longer devices (as shown in Fig. 5, with dWires = 0.65 NWs.µm -2 and Lmask = 1000 µm).

Altogether, the variation of µ0 with temperature was more important when Lmask increased and dWires decreased.

B. Extraction of the thermal activation energy of µ0

We assumed that the difference in the temperature dependence of long channel devices, with respect to short channel devices, is related to the presence of inter-nanowire junctions along the channel. As a first approximation, the junctions were considered as potential barriers with height Ea, as usually done for grain boundaries.

A typical model for taking account of mobility temperature dependence, where current is controlled by a single barrier would be of the form:

µ0 = K exp(-qEa / kbT) (1)
In this equation, K is a constant, kb the Boltzmann constant, q the electron charge in absolute value and T the temperature. The thermal activation energy Ea (in eV) can be extracted from the slope of an Arrhenius plot of µ0 (semilog plot as a function of 1/T). To first order it is equal to the energy barrier height and provides thus a very interesting insight into the governing mechanisms of current conduction in the device. Fig. 7 shows the Arrhenius plot of µ0 for one device with Lmask = 50 µm and dWires = 0.28 nanowires.µm -2 . From this plot, it is possible to extract a mean value of the activation energy Ea_F, describing the overall variation over the range of temperatures from 80 K to 200 K. For the device of Fig. 7, Ea_F = 24 meV. However, the experimental Arrhenius plots did not show a linear trend on the whole range of temperatures. In the example of Fig. 7, the activation energy extracted at low temperature is much lower than that extracted between 200 K and 250 K. To explain such a trend, it is necessary to assume that energy barrier height shows some variability from one junction to another. )

C. Activation energy dispersion

Activation energy dispersion can be explained by local variability of junction sintering. The influence of dispersion has been tentatively explained in [START_REF] Cazimajou | Low Temperature Electrical Characteristics of Si Nanonet Field Effect Transistors[END_REF] with an ideal model made of an infinity of junctions in parallel with a Gaussian distribution for the energy barrier at individual junctions. However this model is not general, nor would be a model with only junctions in series, as a real nanonet normally involves a mixture of series and parallel paths. For qualitative analysis we will discuss Ea_F variations while full simulation will be provided in next section.

On Fig. 8, Ea_F decreases with dWires and increases with Lmask.

In this case, activation energy of the path is the one of the junctions involved in the path with the higher value of activation energy. b) dWires increases. In this case, activation energy of the network can be only inferior to the one of the figure a). A decrease in Lmask has a similar effect, with the predominance of single conductive path for long Lmask, and of parallel paths for short Lmask.

On Fig. 6, Vth increases with temperature increase. For classical structure of FETs, this result is known and is linked to a variation of surface potential. [START_REF] Jeon | Low-temperature electrical characterization of junctionless transistors[END_REF] But this variation is not visible for devices with Lmask = 5 µm, so there is probably a supplementary contribution due to the junctions. Because of Ea dispersion for individual junctions, there is a variation of electrically active density with temperature (some insulating junctions becoming more conductive as temperature increases). Based on [START_REF] Cazimajou | Electrical characteristics of silicon percolating nanonet-based Field Effect Transistors in the presence of dispersion[END_REF] it is known than an increase of the density will lead to a more positive Vth, if there is a dispersion of Vth for individual nanowires. To summarize, an increase of the temperature will lead to an increase of the electrically active area, and thus to a more positive Vth.

In conclusion, taking account of an activation energy dispersion for individual junctions could explain the temperature dependence of µ0 and Vth and the variations of Ea_F with dWires and Lmask. However there is no simple model to confirm this qualitative description. A Monte-Carlo simulation of NN FET mobility as a function of temperature was thus performed in order to evaluate the influence of junction barrier heights and of their dispersion.

D. NN FET mobility Monte Carlo simulation for different temperatures

This simulation was based on the temperature dependent model of nanowires and junctions, based on experimental results.

Junctions are considered as temperature dependent resistance, which vary following Equation ( 2):

Gjunction_i = G0 exp(-Ea_i / kbT) (2)
In this equation Gjunction_i is the conductance of the junction i and Ea_i activation energy of the junction i. Ea_i will be randomly generated according to the chosen density of probability. The value of G0 was chosen so as to have a good order of magnitude for junction resistance, close to 10 MΩ at 300 K. For simplicity, threshold voltage dispersion for individual nanowires is not taken in account in the model. For a chosen Vg, all the nanowires of the nanonet were considered to have the same conductivity, the temperature dependence of which was neglected based on the temperature dependence found in short channel devices ("nanowires in parallel" regime).

Several nanonet structures were built, for each geometry (given values of Source-Drain distance Lmask and pad width Wmask), by randomly generating a number NNW of nanowires defined by the position of their centre (2D coordinates x0i and y0i) and their angle (θi) with respect to horizontal axis (taken along source-drain direction). The junctions between the nanowires were located. Then, a Kirchhoff matrix of conductance was filled, using the above mentioned models for nanowire segments and inter-nanowire junctions A similar matrix filling technique has been used for the study of junction dispersion and structural parameter influence on CNT-network conductance, however with a temperature independent purely resistive model for the junctions. [START_REF] Simoneau | Influence of statistical distributions on the electrical properties of disordered and aligned carbon nanotube networks[END_REF] Finally, this matrix was used to obtain the NN FET conductivity between source and drain. The current was calculated for a large negative gate voltage (ON-current Ion).

Here, low field mobility is proportional to Ion since threshold voltage dispersion was neglected so that we will discuss indifferently Ion and µ0 variations. [START_REF] Cazimajou | Electrical characteristics of silicon percolating nanonet-based Field Effect Transistors in the presence of dispersion[END_REF] In the next part, different densities of probability were tested for activation energy dispersion.

1) Gaussian dispersion

As a first assumption, activation energy dispersion was modelled as a Gaussian distribution with a mean value Ea0 and a standard deviation σEa. For a Gaussian dispersion of activation energy, Ion varied linearly with the inverse of temperature. The extracted Ea_F was closed to Ea0. Increase of σEa lead to a small decrease of Ea_F.

This figure does not fit the experimental data showed on Fig. 7. This can be due to the characteristics of the Gaussian distribution. Indeed, to keep some physical meaning, it had to be truncated to zero. If the standard deviation is not small enough compared to the mean value, this process generates a large number of null barrier height values, which would drive the distribution far from a strictly speaking Gaussian one. It was thus necessary to keep σEa small enough compared to Ea0. However, in this case, the amount of very small barrier heights is small and the chance to get a path that can connect source and drain through such low barrier height junctions runs towards zero. This explains why the plateau at low temperature is not observed in the Arrhenius plot of Id in the frame of this assumption.

2) Bimodal dispersion

In order to explain experimental results, we formulated another hypothesis, where activation energies would be distributed along a bimodal function, with 2 families of inter-nanowire junctions in the network. The presence of several families is not unconceivable as the quality of the sintering process may be non-uniform, depending on several parameters such as nanowire diameter, native oxide thickness and removal efficiency, among others. We retained the simplest picture where good sintering would be obtained in some places while a residual barrier would remain elsewhere.

In order to test this assumption against experimental trends, our Monte-Carlo simulations are now assuming that there is no energy barrier for a fraction of the inter-nanowire junctions while, for the remaining junctions, there is an energy barrier which follows a Gaussian law. The parameter P0 refers to the fraction of junctions with a zero value of Ea. P0 can vary between 0 % (all the junctions in the network have an activation energy which follows the Gaussian law) and 100 % (all energy barriers at inter-nanowire junctions are equal to 0 eV).

The operation mechanism is schematically described in Fig. 11. In this figure, the black dots represent the junctions which feature an energy barrier. At low temperature, the conduction paths which involve these junctions are blocked and the only conducting paths are those which involve purely resistive junctions (no energy barrier). Current flows through the sub-network of nanowires connected by such junctions and a plateau with low activation energy is seen. If P0 is too small or if the nanowire density is too low, there is a risk that this sub-network is below the percolation threshold and the low activation regime is not seen. At higher temperature, the other junctions start to be conductive and current can flow through the whole network, with a hierarchy of conducting paths which can be completely different from what it is at low temperature. We started the analysis by studying the influence of P0 on the temperature dependence of Ion (Figure 12). Five different random realizations were drawn by lot with the same channel length and nanowire density (Lmask = 30 µm and dWires = 0.13 nanowires.µm -2 ). They are differentiated by curve colours in the plots of Figure 12. The parameters of the junctions were then chosen randomly, with several values of P0 and the same Gaussian distribution for temperature activated junctions (Ea0 = 60 meV and σEa = 10 meV). It should be noted that for small value of P0, some Ion versus 1/T plots do not reach a plateau at low temperature. These behaviours can be explained by the fact than in such case, it may happen that there is no conduction path solely made of non-activated junctions. Then, even at low temperature, the temperature dependence of Ion is mainly driven by the junctions with high value of Ea.

It is then interesting to know if this model can provide the correct trends for the variation of the activation energy as a function Lmask and dWires. For this study, the activation energy was extracted as experimentally as the mean value within the range of temperatures from 80 K and 200 K and the parameters for Ea distribution was kept the same for the random generation of energy barriers at inter-nanowire junctions (P0 = 50 %, Ea0 = 60 meV and σEa = 10 meV). P0 is chosen equal to 50 % in order to observe the influence of non activated junctions to the temperature dependence of Ion. Fig. 13 shows Ea_F variation as a function of dWires for different Lmask values. Ea_F obtained by simulation decreases with dWires and increases with Lmask at large NW density, which is the trend observed experimentally.

CONCLUSIONS

We studied the evolution with temperature of the electrical characteristics of nanonet-based FET devices with different source-drain distances and nanowire densities. The analysis was supported by a modelling of carrier transport through a random network of nanowires. The role of inter-nanowires junctions was evidenced. In order to assess their influence on transport, a simple model, which assumed the presence of energy barriers at the junctions, was proposed. This model was updated, with an account for the dispersion of energy barriers values, so as to capture the distinctive features

Figure 1 .

 1 Figure 1. (a) SEM image of a nanonet made of VLS grown Si-Nanowires after passivation step[START_REF] Legallais | Material engineering of percolating silicon nanowire networks for reliable and efficient electronic devices[END_REF]. (b) Top view SEM image of the fabricated devices[START_REF] Cazimajou | Electrical characteristics of silicon percolating nanonet-based Field Effect Transistors in the presence of dispersion[END_REF]

Figs. 2 -

 2 (a) and 2-(c) show typical transfer characteristics of two devices (with Lmask = 5 µm and Lmask = 50 µm, respectively) for different temperatures for the smallest studied density (0.28 NWs.µm -2 ). Figs.2-(b), 2-(d) and 2-(e) are plotting transfer characteristics of three devices (with Lmask = 5 µm, Lmask = 50 µm and Lmask = 1000 µm respectively) for different temperatures for the largest studied density (0.65 NWs.µm -2

Figure 2 .

 2 Fig.2-(c) reveals that the transfer characteristics dependence with temperature was strongly influenced by source drain

Figure 3 .

 3 Figure 3. Effective low field mobility as a function of temperature for a NW density of 0.28 NWs.µm -2 . Full symbols correspond to short devices (Lmask = 5 µm), while empty symbols correspond to long devices (Lmask = 50 µm).

Figure 4 .

 4 Figure 4. Effective low field mobility as a function of temperature for a NW density of 0.65 NWs.µm -2 . Full symbols correspond to short devices (Lmask = 5 µm), while empty symbols correspond to long devices (Lmask = 50 µm).

Figure 5 .

 5 Figure 5. Effective low field mobility as a function of temperature for a NW density of 0.65 NWs.µm -2 , and for a source-drain distance Lmask = 1000 µm.

Fig. 6

 6 Fig.6 plots the temperature dependence of Vth for Lmask = 50 µm and two different values of density, with Vth becoming more positive with temperature increase.

Figure 6 .

 6 Figure 6. Threshold Voltage Vth as a function of temperature for Lmask = 50 µm and two values of dWires. Filled symbols are for dWires = 0.28 nanowires.µm -2 and empty symbols are for dWires = 0.65 nanowires.µm -2 . The lines are only guides for the eyes.

Figure 7 .

 7 Figure 7. µ0 as a function of reciprocal temperature for one device with Lmask = 50 µm and dWires = 0.28 nanowires.µm -2 . Square symbols are for experimental data. Red line (resp. green line) are obtained by fit of experimental data with the equation (1) for the range of temperature from 80 K to 100 K (resp. 200 K to 250 K). Black line is obtained by a fit of experimental data with equation (1) for the range of temperatures from 80 K to 200 K.

Fig. 8 shows

 8 Fig.8 shows Ea_F as a function of dWires for different values of Lmask. Ea_F decreased with dWires and globally increased with Lmask.

Figure 8 .

 8 Figure 8. Activation energy Ea_F function of the nanowire density for different values of Lmask (extracted as from solid line in Fig.7)

Fig. 9 -

 9 (a) shows an ideal case where electrons can only pass through one conduction path made of nanowires in series. This ideal case can be only obtained for a small dWires (or a long Lmask). If there is a dispersion of Ea for individual junctions, Ea_path, energy of the path, will be equal to the maximum of Ea involved in the path. If dWires increases (or Lmask decreases), parallel conduction paths will be formed between source and drain (Fig.9-(b)). In this case, activation energy of the nanonet can only decrease, if one of the parallel paths has an activation energy smaller than Ea_path.

Figure 9 .

 9 Figure 9. a) Schematic of the ideal conduction path obtained for small dWires. In this case, activation energy of the path is the one of the junctions involved in the path with the higher value of activation energy. b) dWires increases. In this case, activation energy of the network can be only inferior to the one of the figure a). A decrease inLmask has a similar effect, with the predominance of single conductive path for long Lmask, and of parallel paths for short Lmask.

  Fig.10 shows the temperature dependence of Ion for Lmask = 30 µm and dWires = 0.13 NW.µm -2 and two Gaussian distribution functions, one with Ea0 = 120 meV and σEa = 40 meV and the other one with Ea0 = 120 meV and σEa = 20 meV.

Figure 10 .

 10 Figure 10. Ion simulated value for Vg = -20 V function of the invert of the temperature for 2 different distributions of Ea (Ea0 = 120 meV and σEa = 40 meV in black and Ea0 = 120 meV and σEa = 20 meV in red), for Lmask = 30 µm and dWires = 0.13 nanowires.µm -2 . 5 realizations are plotted for each distribution. Dotted lines correspond to the linear fit mean values.

Figure 11 .

 11 Figure 11. Schematic representation of the conduction paths across an NN FET as a function of temperature under the assumption that there is an energy barrier at the interface of some of the inter-nanowire junctions (marked with black dots) while the other ones have no energy barrier. The red lines are conductive paths, black lines are insulating, blue and green lines are highly resistive paths. At low temperature, current can flow only through the small number of conductive paths that involve nanowires separated from each other by purely resistive junctions (no energy barrier). As temperature increases, the junctions with small values of Ea become more conductive and current can flow through additional paths that involve junctions with small and large values of Ea.There is an increase of the electrically active density of nanowires with the temperature, without a change of the physical density of nanowires.

Figure 12 .

 12 Figure 12. Id simulated value for Vg = -20 V function of the invert of the temperature for different values of P0, a configuration Lmask = 30 µm and dWires = 0.13 nanowires.µm -2 and a chosen Gaussian Law with Ea0 = 60 meV and σEa = 10 meV. P0 can vary between 0 % (all the junctions in the network have an activation energy which follows the Gaussian law) and 100 % (all energy barriers at inter-nanowire junctions are equal to 0 eV)

Figure 13 .

 13 Figure 13. Ea_F extracted from the simulation as a function of dWires and Lmask for P0 = 50 %, Ea0 = 60 meV and σEa = 10 meV.

  observed experimentally. It was then implemented within a Monte-Carlo description of the Si-NN. Comparison of simulation results with experimental trends steered the analysis towards a picture where a significant proportion of the junctions would be purely resistive, while the others would have an energy barrier. It is concluded that, with proper analysis of electrical characteristics dependence with temperature, source-drain distance and nanowire density, it is possible to get some insight, and thus to provide some technological feedback, about the role of inter-nanowire junction quality in nanonet-based devices.

ACKNOWLEDGMENT

This work has received funding from the EU H2020 RIA project Nanonets2Sense under grant agreement n°688329. This work was partly supported by the French RENATECH network. This work benefited from the facilities and expertise of the OPE)N(RA characterization platform of FMNT (FR 2542, fmnt.fr) supported by CNRS, Grenoble INP and UGA.


W Larbi

Numerical modeling of sound and vibration reduction using viscoelastic materials and shunted piezoelectric patches

Keywords: Vibroacoustic, Sandwich structure, Piezoelectric shunt, Finite Element, Boundary Element, Modal reduction

In this paper, a reduced order coupled finite element/boundary element method (FEM/BEM) for control of noise radiation and sound transmission of vibrating structure by viscoelastic and passive piezoelectric treatments is presented. The system consists of a sandwich structure composed of elastic faces and a viscoelastic core (with surface mounted piezoelectric patches) and coupled to external/internal acoustic domains. The passive shunt damping strategy is employed for vibration attenuation in the low-frequency range while dissipative viscoelastic materials are used for the higher-frequency domain. The originality of the present paper lies (i) in evaluating the classically used FEM/BEM methods for structural-acoustics problems when taking account smart systems at the fluid-structure interfaces and viscoelastic interlayers and (ii) in development of a reduced-order model, based on a normal mode expansion, able to capture the main characteristics of the system dynamical behavior.

Introduction

Sandwich structures with a viscoelastic layers are commonly used in many applications for vibration damping and noise control. In such structures, the main energy loss mechanism is due to the transverse shear of the viscoelastic core. However, accurate modeling of structures with viscoelastic materials is difficult because the measured dynamic properties of viscoelastic material are frequency and temperature dependent. This motivated several authors to develop accurate numerical methods of modeling the effects of viscoelastic damping mechanisms which introduce frequency dependence [START_REF] Ferreira | A finite element model using a unified formulation for the analysis of viscoelastic sandwich laminates[END_REF][START_REF] Madeira | Multiobjective design of viscoelastic laminated composite sandwich panels[END_REF][START_REF] Araújo | Finite element model for hybrid active-passive damping analysis of anisotropic laminated sandwich structures[END_REF][START_REF] Araújo | Optimal design for active damping in sandwich structures using the Direct MultiSearch method[END_REF]. A review of these methods can be found in [START_REF] Vasques | Viscoelastic damping technologiesPart I: Modeling and finite element implementation[END_REF]. Concerning the application of these structures in noise and vibration attenuation, we can cite for example [START_REF] Dijckmans | Sound transmission through finite lightweight multilayered structures with thin air layers[END_REF][START_REF] Akrout | Vibro-acoustic behaviour of laminated double glazing enclosing a viscothermal fluid cavity[END_REF]. In [START_REF] Dijckmans | Sound transmission through finite lightweight multilayered structures with thin air layers[END_REF], the measured sound transmission loss of multilayered structures is compared with transfer matrix method results (assuming infinite layers) and a wave based model (taking into account finite dimensions) to show the importance of the finite dimensions in a broad frequency range. The effects of viscothermal fluid in laminated double glazing are investigated in [START_REF] Akrout | Vibro-acoustic behaviour of laminated double glazing enclosing a viscothermal fluid cavity[END_REF] using a finite element approach. However, at low frequency, the acoustic performance of this type of systems is greatly deteriorated and the viscoelastic layer is not effective in treating the fall of sound transmission loss and the vibration reduction. The aim of this reported research work is to reduce the vibration and the sound transmission at low resonance frequencies by a passive piezoelectric shunt technique through a fully coupled finite element/boundary element modeling of the problem. In this technology, the structure is equipped with piezoelectric patches that are connected to a passive electrical circuit, called a shunt. The piezoelectric patches transform mechanical energy of the vibrating structure into electrical energy, which is then dissipated by Joule heat in the shunt circuits. Several shunt circuits can be considered: the classical R-and RL-shunts, proposed by Hagood and Von Flotow [START_REF] Hagood | Damping of structural vibrations with piezoelectric materials and passive electrical network[END_REF] and improvements of those techniques, using several piezoelectric elements [START_REF] Dellisola | Passive damping of beam vibrations through distributed electric networks and piezoelectric transducers: prototype design and experimental validation[END_REF][START_REF] Collet | Wave motion optimization in periodically distributed shunted piezocomposite beam structures[END_REF][START_REF] Casadei | Broad band vibration control through periodic arrays of resonant shunts: experimental investigation on plates[END_REF], active fiber composites [START_REF] Belloli | Structural vibration control via RL-shunted active fiber composites[END_REF] or adaptive shunts [START_REF] Niederberger | Adaptive multi-mode resonant piezoelectric shunt damping[END_REF], and recently semi-passive techniques, commonly known as switch techniques [START_REF] Guyomar | Wave reflection and transmission reduction using piezoelectric semipassive nonlinear technique[END_REF][START_REF] Ducarne | Structural vibration reduction by switch shunting of piezoelectric elements: modeling and optimization[END_REF][START_REF] Larbi | Structural-acoustic vibration reduction using switched shunt piezoelectric patches: a finite element analysis[END_REF]. As those techniques are passive (or semi-passive if some electronic components have to be powered), a critical issue is that their performances, in terms of damping efficiency, directly depend on the electromechanical coupling between the host structure and the piezoelectric elements, which has to be maximized and necessitates the development of predictive models.

The present work concerns the numerical modeling of noise and vibration reduction of laminated structures with viscoelastic interlayers by using shunted piezoelectric elements. The frequency domains of interest are the low and medium frequency ranges. The aim is to propose an efficient reduced order coupled finite element/boundary element model able to predict the shunt damping around the resonance frequencies of the system. In the first part of this paper, a finite element formulation (FEM) of sandwich structures with viscoelastic core and equipped with shunted piezoelectric patches is presented. This formulation involves structural displacement in the structure (sandwich structure with piezoelectric elements), acoustic pressure in the fluid cavity and the electric charge and voltage between the electrodes in the piezoelectric patches. The charge/voltage variables are intrinsically adapted to include any external electrical circuit into the electromechanical problem and to simulate the effect of shunt damping techniques. Moreover, since the elasticity modulus of the viscoelastic core is complex and frequency dependent, this formulation is complex and nonlinear in terms of frequency. The direct solution of this problem can only be considered with models which do not imply a prohibitive number of degrees of freedom. This has severe limitations in attaining adequate accuracy and wider frequency ranges of interest. A reduced order-model is then proposed to solve the problem at a lower cost in the second part of this paper. The proposed methodology, based on a normal mode expansion, requires the computation of the uncoupled structural and acoustic modes. The uncoupled structural modes are the real and undamped modes of the sandwich structure without fluid pressure loading at fluid-structure interface and with short-circuited and open-circuited patches, whereas the uncoupled acoustic modes are the cavity modes with rigid wall boundary conditions at the fluid-structure interface. It is shown that the projection of the full-order coupled finite element model on the uncoupled bases, leads to a reduced order model in which the main parameters are the classical fluid-structure coefficient, the residual stiffness complex coupling factors and the electromechanical coupling factors. Because of its reduced size, this model is proved to be very efficient for simulations of steady-state and frequency analyses of the fully coupled visco-electro-mechanical-acoustic system and the computational effort is significantly reduced. Note that the computing of eigenmodes in the medium frequency range presents no difficulty in this work. Indeed, the new numerical computing tools allow easy access to these modes which has not been the case before. As a next step, the direct boundary element method (BEM) is used for modeling the scattering/radiation of sound by the structure coupled to an external acoustic domain. The BEM is derived from a boundary integral equation involving the surface pressure and normal acoustic velocity at the boundary of the acoustic domain. The coupled FEM-BEM model is obtained by using a compatible mesh at the fluid-structure interface. In the last part of the paper, a numerical example is presented in order to validate and analyze results computed from the proposed formulations.

Numercial modeling of the structure

Consider a thin sandwich structure, denoted Ω E , made of elastic faces and viscoelastic core (Figure 1). A prescribed force density F d is applied to the external boundary Γ σ of Ω E and a prescribed displacement u d is applied on a part Γ u of Ω E . The structure is coupled to internal and external fluids denoted Ω F i and Ω F e , respectively. These fluids are considered homogeneous, inviscid and compressible and the gravity effects are neglected. We denote by Σ i and Σ e the internal and external fluid-structure interfaces, respectively, such that Σ i ∪ Σ e is a partition of the all fluid-structure interface Σ. Moreover, n S and n F are the normal units external to the solid and fluid domains, respectively. In order to achieve maximum vibration dissipation and acoustic radiation attenuation of selected modes, the passive piezoelectric shunt damping technique is used. Thus, a set of P piezoelectric patches are bonded to the structure surface and connected to resistive or resonant shunt circuits. In this technology, the piezoelectric patches convert a fraction of mechanical energy associated with the structure vibration into electrical energy, which is dissipated by heat through the resistor in the shunt circuits. Each piezoelectric patch has the shape of a plate with its upper and lower surfaces covered with very thin layer electrodes. The pth patch, p ∈ {1, • • • , P }, occupies a domain Ω (p) such that (Ω E , Ω (1) , • • • , Ω (P ) ) is a partition of the all structure domain Ω S . Moreover, we denote by R (p) and L (p) the resistance and the inductance of the resonant shunt circuit connected to the pth patch.

Γ u F d External Fluid Internal Fluid R (1) L (1) R (p) L (p) Sandwich Structure Viscoelastic core Elastic faces Shunted Piezo. Patch Ω E Ω (p) Ω Fi Ω Fe Γ σ Σe Σi

Modeling of the viscoelastic core

According to the theory of linear viscoelasticity [START_REF] Christensen | Theory of Viscoelasticity: An Introduction[END_REF] applied to anisotropic viscoelastic material, the one-dimensional stress σ is related to the strain history ε by a characteristic relaxation function G(t):

σ(t) = t -∞ G(t -τ ) ∂ε ∂τ dτ (1)
Under steady-state harmonic oscillation conditions and using the Fourier transform, the constitutive equation becomes

σ * (ω) = G * (ω)ε * (ω) ( 2
)
where ω is the angular frequency, G * (ω) is the so-called complex modulus characterizing the viscoelastic material, which can be expressed as

G * (ω) = G (ω)(1 + iη(ω)) (3)
with G (ω) = Re(G * (ω)) the storage modulus, η(ω) = Im(G * (ω))/Re(G * (ω)) the loss factor and i 2 = -1.

Many viscoelastic models aim at representing the frequency-dependence of viscoelastic materials' properties by a complex shear modulus G * , which depends on frequency and a series of parameters to be identified on master curves. The Poisson coefficient is commonly supposed constant in the studied frequency range, so that the Young modulus and the shear modulus are proportional. Some of the most well-known models are the Golla-Hughes-McTavish (GHM) [START_REF] Golla | Dynamics of viscoelastic structures-a time domain finite element formulation[END_REF], the anelastic displacement field (ADF) [START_REF] Lesieutre | Time domain model ingof linear viscoelastic ityusinganelastic displacement fields[END_REF] models, and the fractional derivative model [START_REF] Bagley | Fractional calculus -a different approach to analysis of viscoelastically damped structures[END_REF][START_REF] Galucio | Finite element formulation of viscoelastic sandwich beams using fractional derivative operators[END_REF]. Their corresponding expressions of the complex shear modulus are presented in Table 1 in terms of various parameters.

The frequency dependent material properties of the polyvinyl butyral (PVB) were determined by doing a Havriliak-Negami fit of data found using a Dynamic Mechanical Thermal Analyzer. The fractional derivative model enables a good representation of the viscoelastic properties of the PVB (Figure 2) with a minimum number of parameters (G0 = 0.479 MPa

Viscoelastic model

Complex shear modulus

GHM G * (ω) = G 0 1 + J j=1 α j -ω 2 +2iξ j ω j ω -ω 2 +2iξ j ω j ω+ω 2 j ADF G * (ω) = G 0 1 + J j=1 ∆ j ω 2 +iωΩ j ω 2 +Ω 2 j Fractional derivative G * (ω) = G ∞ + (G 0 -G ∞ ) [1 + (iωτ ) 1-α ] -β
Table 1: Complex shear modulus expressions of GHM, ADF and fractional models. and G ∞ = 0.235 GPa are the asymptotic values of the shear modulus when ω -→ 0 and ω -→ ∞, respectively, τ = 0.3979 is the relaxation time and α = 0.46 and β = 0.1946 are the parameters of the fractional derivative). The choice of the fractional derivative model minimizes the number of variables which is of practical importance for future optimization of the viscoelastic properties. Note that this material will be used for the numerical example of this paper.

Modeling of the piezoelectric patches

For each piezoelectric patch, a set of hypotheses, which can be applied to a wide spectrum of practical applications, are formulated:

• The piezoelectric patches are thin, with a constant thickness, denoted h (p) for the pth patch;

• The thickness of the electrodes is much smaller than h (p) and is thus neglected;

• The piezoelectric patches are polarized in their transverse direction (i.e. the normal direction to the electrodes).

Under those assumptions, the electric field vector, E, can be considered normal to the electrodes and uniform in the piezoelectric patch [START_REF] Thomas | Vibrations of an elastic structure with shunted piezoelectric patches: efficient finite elements formulation and electromechanical couplings coefficients[END_REF], so that for all p ∈ {1, • • • , P }:

E = - V (p) h (p) n in Ω (p) (4)
where V (p) is the potential difference between the upper and the lower electrode surfaces of the pth patch which is constant over Ω (p) and n is the normal unit vector to the surface of the electrodes. Moreover, the linear piezoelectric constitutive equations are:

σ = c ε -e T E (5) D = e ε + E (6)
where ε = 1 2 ∇u + ∇ T u is the linearized deformation tensor, u is the mechanic displacement field, σ is the stress tensor, D denotes the electric displacement, c denotes the elastic moduli at a constant electric field, e denotes the piezoelectric constants, and denotes the dielectric permittivities at constant strain [START_REF] Gaudenzi | Smart structures: physical behaviour, mathematical modelling and applications[END_REF].

Finite element formulation of sandwich structure with piezoelectric shunt systems 2.3.1. System equations in the frequency domain

The local equations of sandwich structure with piezoelectric patches and submitted to an acoustic pressure p are

div σ + ω 2 ρ S u = 0 in Ω S (7a) σ n S = F d on Γ σ (7b) u = 0 on Γ u (7c) σ n S = p n F on Σ i (7d) div D = 0 in Ω (p) (8a) D • n S = Q d on Γ (p) D (8b)
where ρ S is the structure mass density. Equation (7a) corresponds to the elastodynamic equation in the absence of body force; Equations (7b) and (7c) are the classical mechanical boundary conditions; Equation (7d) results from the action of pressure forces exerted by the internal fluid on the structure; Equation (8a) corresponds to the electric charge equation for a dielectric medium; Equation (8b) corresponds to the electric boundary condition for the pth patch defined by a prescribed surface density of the electric charge

Q d on Γ (p) D .
Note that due to the properties of the viscoelastic core, the stress tensor of the structure is complex and frequency dependent. It will be noted σ * (u, ω).

Variational formulation

By considering successively each of the P +1 subdomains (Ω E , Ω (1) , • • • , Ω (P ) ), the variational formulation of the sandwich structure/piezoelectric-patches coupled system can be written in terms of the structural mechanical displacement u, the electric potential difference V (p) constant in each piezoelectric patch, and the internal fluid pressure p at the fluid-structure interface Σ i [START_REF] Larbi | Vibration of axisymmetric composite piezoelectric shells coupled with internal fluid[END_REF][START_REF] Larbi | Finite element formulation of smart piezoelectric composite plates coupled with acoustic fluid[END_REF][START_REF] Deü | Piezoelectric structural acoustic problems: Symmetric variational formulations and finite element results[END_REF]]

Ω S σ * : ε(δu) dv + P p=1 V (p) h (p) Ω (p) e (p)T n (p) : ε(δu) dv -ω 2 Ω S ρ S u • δu dv - Σ i p n F • δu ds = Γσ F d • δu ds ∀δu ∈ C * u (9)
where the admissible space C * u is defined by (p) defines the capacitance of the pth piezoelectric patch (S (p) being the area of the patch and (p) [START_REF] Trindade | On the finite element characterization of d31 macro-fibre composites effective properties[END_REF] the piezoelectric material permittivity in the normal direction to the electrodes) and Q (p) is the global charge in one of the electrodes (see [START_REF] Thomas | Vibrations of an elastic structure with shunted piezoelectric patches: efficient finite elements formulation and electromechanical couplings coefficients[END_REF]).

C * u = {u ∈ C u | u = 0 on Γ u }. C u being the admissible space of regular functions u in Ω S . P p=1 δV (p) C (p) V (p) - P p=1 δV (p) h (p) Ω (p) e (p) : ε(u) •n (p) dv = P p=1 δV (p) Q (p) ∀δV (p) ∈ R (10) where C (p) = (p) 33 S (p) /h

Finite element discretization

After discretization of the previous variational formulation by the finite element method and using the following additional relation between electrical potential differences and electric charges due to the shunt circuits:

-ω 2 LQ -iωRQ + V = 0 (11)
we find the following matrix equation:

K * u (ω) + C uV K -1 V C T uV C uV K -1 V -C up K -1 V C T uV K -1 V 0   U Q P   -iω 0 0 0 0 R 0   U Q P   -ω 2 M u 0 0 0 L 0   U Q P   = F 0 ( 12
)
where

Q = (Q (1) Q (2) • • • Q (P ) ) T and V = (V (1) V (2) • • • V (P )
) T are the column vectors of electric charges and potential differences;

R = diag R (1) R (2) • • • R (P ) and L = diag L (1) L (2) • • • L (P )
are the diagonal matrices of the resistances and inductances of the patches; U and P are the vectors of nodal values of u and p; M u and K * u (ω) are the mass and stiffness matrices of the structure (sandwich structure and piezoelectric patches); C uV is the electric mechanical coupled stiffness matrix; ) is a diagonal matrix filled with the P capacitances of the piezoelectric patches; C up is the fluidstructure coupled matrix; F is the applied mechanical force vector. Note that the stiffness matrix K * u (ω) is complex and also frequency dependent.

K V = diag C (1) C (2) • • • C (P

Reduced order model

We present in this section a reduced-order formulation of the discretized problem obtained in the previous section. This formulation consists in expanding the structural displacement over the in vacuo real structure modes in short-circuited and open-circuited configurations.

In a first phase, the first N s eigenmodes of the structure in vacuo with all patches short circuited are obtained from [3]

K * u (ω) -ω 2 M u U = 0 (13)
Due to the frequency dependence of the stiffness matrix, this eigenvalue problem is complex and nonlinear. It is assumed that vibrations of the damped structure can be represented in terms of the real modes of the associated undamped system if appropriate damping terms are inserted into the uncoupled modal equations of motion [START_REF] Fan | Transient vibration and sound radiation of a rectangular plate with viscoelastic boundary supports[END_REF][START_REF] Bouayeda | Finite element analysis of the dynamic behavior of a laminated windscreen with frequency dependent viscoelastic core[END_REF]. Thus, the complex stiffness matrix is decomposed in the sum of two matrices:

K * u (ω) = K u0 + δK * u (ω) (14)
where K u0 is the real and frequency-independent stiffness matrix calculated with a constant Young's modulus of the viscoelastic core and δK * u (ω) is the residual stiffness matrix.

The ith real eigenmode is obtained from the following equation

K u0 -ω 2 si M u Φ si = 0 (15)
where (ω si , Φ si ) are the natural frequency and eigenvector for the ith structural mode in short circuited configuration. These modes verify the following orthogonality properties

Φ T si M u Φ sj = δ ij and Φ T si K u0 Φ sj = ω 2 si δ ij ( 16
)
where δ ij is the Kronecker symbol and Φ si have been normalized with respect to the structure mass matrix.

In a second phase, the first N s real eigenmodes of the structure with all patches open circuited are obtained from the following equation

(K u0 + C uV K -1 V C T uV ) -ω 2 i M u Φ i = 0 (17)
where ( ω i , Φ si ) are the natural frequency and eigenvector for the ith structural mode in open circuited configuration. These modes verify the following orthogonality properties

Φ T i M u Φ j = δ ij and Φ T i (K u0 + C uV K -1 V C T uV ) Φ j = ω 2 i δ ij (18)
By introducing the modal matrix Φ s formed by both the short-circuited and open-circuited modes:

Φ s = Φ Φ = Φ 1 • • • Φ Ns Φ 1 • • • Φ Ns (19)
of size (M s × 2N s ), M s is the total number of degrees of freedom in the finite elements model associated to the structure, the displacement U is sought as

U = Φ s q s (ω) (20)
where the vector q s is the unknown modal amplitudes.

By applying the Ritz-Galerkin projection method, which consists of substituting Equation [START_REF] Guyomar | Wave reflection and transmission reduction using piezoelectric semipassive nonlinear technique[END_REF] into Equation [START_REF] Deü | Piezoelectric structural acoustic problems: Symmetric variational formulations and finite element results[END_REF] and premultiplying the first row by Φ T s , we obtain the reduced matrix system,

-ω 2 Φ T s M u Φ s 0 0 0 L 0   q s Q P   + iω 0 0 0 0 R 0   q s Q P   + Φ T s K * u (ω) + C uV K -1 V C T uV Φ s Φ T s C uV K -1 V -Φ T s C up K -1 V C T uV Φ s K -1 V 0   q s Q P   = Φ T s F 0 (21)
Remarks:

• The projection bases widely used in the literature are obtained using only the short-circuited or open-circuited eigenvalue problems. We have shown that the proposed new modal projection technique is more efficient than the classical bases in terms of prediction of vibration reduction with shunted piezoelectric patches [START_REF] Larbi | Vibration reduction of an elastic structure using shunted piezoelectric systems: New reduced order finite element models[END_REF]. This efficiency is related to the introduction of additional electromechanical coupling terms in the reduced system, which allows better modeling of the effects of the piezo-patches on the system response.

• The proposed combined basis requires two modal calculations to find short-circuit and open-circuit modes. The short-circuited and opencircuited eigenvectors are not mutually orthogonal with respect to the mass and stiffness matrices of the problem. Thus, the resulting basis is not orthogonal. The Gram-Schmidt algorithm is applied in this work to orthogonalize these eigenvectors and build well-conditioned matrices.

• The change of variables from physical coordinates U to modal coordinates q s and the modal truncation concern only the mechanical displacement and not the electrical unknown field Q. This is justified by the fact that it is unnecessary to reduce Q since its size is small and represents the number of piezoelectric patches.

Numerical modeling of fluid

This section deals with the numerical modeling of the fluid. The internal fluid is modeled using the finite element method. The external fluid is modeled using the boundary element method.

Equations in the frequency domain

The governing equations of the exterior/interior acoustic fluid are

∆p + ω 2 c 2 F p = 0 in Ω F (22a) ∂p ∂n = ρ F ω 2 u • n F on Σ (22b) lim r→∞ r( ∂p ∂r + ikp) = θ( 1 r ) on Γ ∞ (22c)
where ρ F is the mass density of the fluid, c F is the sound velocity, k = ω/c F is the wave number and the function θ( 1 r ) tends to zero when r tends to infinity.

Equation (22a) represents the Helmholtz equation; Equation (22b) is the kinematic interface fluid-structure condition on Σ; Equation (22c) represents the Sommerfeld condition at infinity (used for the external fluid).

Finite element formulation for internal fluid

Let δp be the test function, associated to p, belonging to the admissible space C p . Multiplying Equation (22a) by δp ∈ C p , applying Green's formula, and finally taking Equation (22b) into account, we obtain:

1 ρ F Ω F i ∇p • ∇δp dv -ω 2 Σ i u • nδp ds - ω 2 ρ F c 2 F Ω F i p δp dv = 0 ∀δp ∈ C p (23)
After discretizing by the finite element method the bilinear forms in Equation [START_REF] Larbi | Structural-acoustic vibration reduction using switched shunt piezoelectric patches: a finite element analysis[END_REF], we obtain the following matrix system:

K p P -ω 2 C T up U -ω 2 M p P = 0 (24)
where M p and K p are the mass and stiffness matrices of the internal fluid.

This last equation is projected on a truncated basis constructed from the first N f acoustic modes. These modes are calculated with rigid boundary conditions from the following equation

K p -ω 2 f i M p Φ f i = 0 (25)
where (ω f i , Φ f i ) are the natural frequency and eigenvector for the ith acoustic mode. These modes verify the following orthogonality properties

Φ T f i M p Φ f j = δ ij and Φ T f i K p Φ f j = ω 2 f i δ ij ( 26
)
where Φ f i have been normalized with respect to the fluid mass matrix.

By introducing the matrix

Φ f = Φ f 1 • • • Φ f N f of size (M f × N f
) corresponding to the acoustic basis (M f is the total number of degrees of freedom in the finite elements model associated to the internal acoustic domain), the pressure is sought as

P = Φ f q f (ω) (27)
where the vector

q f = q f 1 • • • q f N f
T is the vector of modal amplitudes of the fluid pressure.

Substituting this relation into Equation [START_REF] Larbi | Finite element formulation of smart piezoelectric composite plates coupled with acoustic fluid[END_REF] after pre-multiplying it by Φ T f and taking into account Equation ( 20), we obtain

Φ T f K p Φ f q f -ω 2 Φ T f C T up Φ s q s -ω 2 Φ T f M p Φ f q f = 0 (28)
This matrix equation represents the reduced order model of the internal acoustic problem. If only few modes are kept for the projection, the size of this reduced order model (N f × N f ) is much smaller than the initial one (M f × M f ).

Boundary element formulation for external fluid

The boundary element formulation is used for the exterior acoustic domain. The Helmholtz equation is valid for the pressure p at the arbitrary collocation point x within the acoustic domain Ω F e . A weak form of this equation is obtained by weighting with the fundamental solution:

G(x, y) = e ik|x-y| 4π|x -y| (29)
where |x -y| denotes the distance between an arbitrary point x and the load source point y.

Applying Green's second theorem, the Helmholtz equation can be transformed into a boundary integral equation, which can be expressed as follows

c(x)p(x) = Σe p(y) ∂G(x, y) ∂n y dS - Σe ∂p(y) ∂n y G(x, y)dS (30)
where

c(x) =        1 x in fluid domain 1 2
x on smooth boundary of fluid domain

Ω(x) 4π
x on nonsmooth boundary of fluid domain 0 x outside fluid domain [START_REF]Smart Material Corp. Macro-Fibre Composites Technical Data, www.smart-material[END_REF] and

Ω(x) = 4π + Σe ∂(|x -y| -1 )
∂n y dS is the solid angle seen from x. Note that the value c(x) = 1 2 is valid if the surface Σ e is assumed to be closed and sufficiently smooth, i.e. there is a unique tangent to Σ e at every x ∈ Σ e . For the general case, where a nonunique tangent plane exists at x ∈ Σ e , we use c(x) = Ω(x) 4π (for example, when x is lying on a corner or an edge).

The fluid boundary is divided into N quadrilateral elements (∂Ω F = N j=1 S j ) and Equation ( 30) is discretized. After using the relation between the acoustic pressure and the fluid normal velocity ∂p ∂n = -iρ F ωv (where v = v F • n F ), the discrete Helmholtz equation can be written for any point x defined by the node i as

c i p i = N j=1 S j p j (y) ∂G(x i , y) ∂n y dS + iρ F ω N j=1 S j v j (y)G(x i , y)dS (32)
For each quadrilateral element j, the pressure p j (y) and the normal velocity v j (y) can be expressed as a function of their nodal values

p j (y) = 4 k=1 N k p k j = Np j , v j (y) = 4 k=1 N k v k j = Nv j ( 33
)
and Equation (32) becomes

c i p i = N j=1 4 k=1 S j N k ∂G(x i , y j ) ∂n y dSp k j + iρ F ω N j=1 4 k=1 S j N k G(x i , y j )dSv k j ( 34
) or in the following form

c i p i = N j=1 4 k=1 H k ij p k j + iρ F ω N j=1 4 k=1 G k ij v k j ( 35
)
where

H k ij = S j N k ∂G(x i , y j ) ∂n y dS , G k ij = S j N k G(x i , y j )dS ( 36
)
We now place the point x i at each of nodal points on the boundary y j successively, which is known as a "collocation point". We obtain

c i δ ij p j = N j=1 4 k=1 H k ij p k j + iρ F ω N j=1 4 k=1 G k ij v k j (37)
When the collocation scheme is repeated for all nodal points N n of the boundary element mesh, a set of N n expressions in the nodal field variables is obtained which can be assembled into the following matrix equation

HP Σe = iρ F ωGv Σe ( 38
)
where P Σe and v Σe are the vectors with sound pressure and velocity in the normal direction to the boundary surface at the nodal position of the boundary element mesh.

The boundary conditions given in Equation ( 22b) can be expressed in discretized form

v Σe = iωTU ( 39
)
where T is the global coupling matrix that transforms the nodal normal displacement of the structure to the normal velocity of the acoustic fluid at the interface. Substituting this relation into Equation (38), we obtain

HP Σe = -ρ F ω 2 GTU (40)
4. Reduced FE/BE formulation for the fluid-structure coupled problem with piezoelectric shunt systems and viscoelastic materials By combining Equation [START_REF] Hagood | Damping of structural vibrations with piezoelectric materials and passive electrical network[END_REF] with Equations ( 28) and (40), we find the following coupled FE/BE matrix equation

    Φ T s K * u (ω) + C uV K -1 V C T uV Φ s Φ T s C uV K -1 V -Φ T s C up Φ f -Φ T s C up K -1 V C T uV Φ s K -1 V 0 0 0 0 Φ T f K p Φ f 0 0 0 0 H         q s Q q f P Σe     +iω     0 0 0 0 0 R 0 0 0 0 0 0 0 0 0 0         U Q q f P Σe     -ω 2     Φ T s M u Φ s 0 0 0 0 L 0 0 Φ T f C T up Φ s 0 Φ T f M p Φ f 0 ρ F GT 0 0 0         q s Q q f P Σe     =     Φ T s F 0 0 0     (41)
This matrix equation represents the reduced order model of the fully coupled visco-electro-mechanical-acoustic system.

The matrices H and G, associated to the BE formulation for the external fluid, are complex-valued, nonsymmetric, full and frequency-dependent. This represents a main drawback compared to the FEM where matrices are sparse, symmetric, and the frequency dependency can be taken out of the elementary matrices (case of the acoustic interior problem). Moreover, the stiffness matrix of the structure K * u is complex and frequency-dependent. Consequently, the matrices of Equation (41) must be evaluated at each frequency and the dynamic response is computed by a direct method. Despite the reduction of the degrees-of-freedom of the structure and the internal fluid in the system, the use of the BEM for the external fluid can prove to be costly. A possible acceleration of the direct solution can be performed using frequency interpolation methods [START_REF] Atalla | Finite element and boundary methods in structural acoustics and vibration[END_REF] but is beyond the scope of this paper. Furthermore, because the Greens function is a function of 1/r, a special numerical treatment is necessary when the collocation point x coincides with the interpolation point y (see for example [START_REF] Atalla | Finite element and boundary methods in structural acoustics and vibration[END_REF] for more details about the BEM implementation techniques).

Remarks: acoustic indicators

In order to evaluate the vibroacoustic performance and the vibration reduction properties of the viscoelastic and piezoelectric materials, the radiated sound power (Π) in the external acoustic domain and the mean quadratic normal velocity V 2 are used as vibroacoustic indicators in this work.

• The sound intensity I at every point Q on the vibrating surface may be calculated from

I = 1 2 Re {p(Q)v * (Q)} = 1 2 Re p(Q)(iωu(Q) • n S ) * (42)
where Re denotes the real part of the expression in parentheses and * denotes the complex conjugate.

• The intensity is the sound power radiated per unit area of the transmission; the total radiated sound power Π is found by integrating the intensity over the surface Σ e :

Π = Σe I ds = ω 2 Im U * T C up P Σe = ω 2 Im q * T s Φ T s C up P Σe (43)
• The mean quadratic normal velocity V 2 is given by

V 2 = 1 2Σ i Σ i v 2 ds = ω 2 2Σ i w * Mu w ( 44
)
where w is the vector of nodal normal displacements of the structure and Mu is the structure mass matrix reduced to the dofs of interest.

Numerical example

In this last section, numerical results obtained with a Matlab program developed by the author are proposed in order to validate and analyze the results calculated from the proposed formulation.

We consider a flexible sandwich plate (in-plane (x, y) dimensions are A= 1.5 m, B=1.25 m) composed of two identical glass plates bonded together by a PVB interlayer. The thickness of the outer and inner glass ply is 3 mm and those of the PVB interlayer is 1.14 mm (figure 3). The inner face of the plate is coupled to a 3D hexahedric acoustic cavity of size A, B and C=1 m along the directions x, y, and z, respectively. The cavity is completely filled with air (density ρ=1.2 kg/m 3 and speed of sound c=340 m/s). The cavity walls are rigid except the top one, which is in contact with the sandwich plate.

The outer face of the plate is in an infinite baffle, radiating into an infinite half-space acoustic domain filled with air. The glass ply is modeled as linear elastic material (density 2500 kg/m 3 , Young's modulus 72 GPa, and Poisson ratio 0.22). The material properties of the PVB are both thermal and frequency dependent. The complex PVB shear modulus is represented by the fractional derivative model given in section 2.1 (the Poisson ratio is 0.4 and density is 999 kg/m 3 ). A normal incidence plane wave with an amplitude of 1 N/m 2 is applied to the outer face of the sandwich plate as the only external force to the system.

In this example, the structure and the acoustic cavity are represented by FE reduced order model obtained by projection the mechanical displacement unknown on a truncated basis composed by the first real "in vacuo" normal modes and the internal pressure unknown on a truncated basis composed by the first acoustic modes with rigid boundaries conditions.

A modal damping coefficient η= 0.001 is applied for the computation of the dynamic response of the plate with and without viscoelastic interlayer. This coefficient represents the damping in the glass layers independently of the viscoelastic damping. It is applied for all eigenmodes in the selected structural modal basis. The sound radiated by the structure in the external acoustic domain is represented by the BEM.

Concerning the finite element discretization, we have used, for the structural part, 10×10 rectangular elements for the sandwich plate developed by the author in [START_REF] Larbi | Vibroacoustic behavior of a laminated plate with frequency dependent viscoelastic core: finite element and experimental analysis[END_REF]. This element is based on the layerwise theory with a first order shear deformation in each layer. The skins are described according to the Kirchhoff-Love theory with a correction which takes into account the rotational influence of the transversal shearing in the core. A Mindlin model is used to describe the displacement field of the core. The element has four nodes with seven degrees-of-freedom (dof) per node. The in-plane displacements and the rotations of the core are discretized by conforming bi-linear Lagrange shape functions, while the transverse displacement and rotation of the skins are discretized by non-conforming cubic Hermite shape functions. This choice was proved to be efficient compared to other FE models and well adapted to structural-acoustic applications. The acoustic cavity is discretized using 10×10×10 hexahedric elements with one degree-of-freedom per node, corresponding to the acoustic pressure. The structural and acoustic meshes are compatible at the interface, and the fluid-structure coupling is realized through the C up matrix.

For the BEM formulation of the external acoustic fluid domain and the coupling with the sandwich plate, besides the FE discretization of the plate, only the boundary Σ e of the external fluid domain is discretized with boundary elements. Notice that the BE nodes on this part must coincide with FE nodes.

Influence of internal fluid

In order to study the influence of the acoustic cavity on the vibrations of the plate, we consider an equivalent glass plate with the same mass (without viscoelastic interlayer and obtained by adjusting its thickness) coupled to the acoustic cavity. Table 2 presents the first ten eigenfrequencies in three cases: (i) the 3D rigid acoustic cavity, (ii) the plate without fluid and (iii) the plate/acoustic-cavity coupled system in the short circuited case. All results are computed from the finite element formulation presented in [START_REF] Larbi | Vibroacoustic behavior of a laminated plate with frequency dependent viscoelastic core: finite element and experimental analysis[END_REF].

The sixth and ninth frequencies of the coupled system are associated with the first two acoustics modes in the rigid cavity lower than 280 Hz, while the other frequencies correspond to the first eight vibration modes of the structure. This can be confirmed by comparing the mode shapes in case (iii) with those obtained in case (i) or case (ii), which are not shown here for the sake of brevity. Moreover, as expected, the natural frequencies of the coupled modes (structure dominated) are lower than those for the structure in vacuum (except for the first mode) due to the added-mass effect of the fluid.

Fluid

Influence of the viscoelastic layer

A comparison between a simple glass and a laminated glass with PVB interlayer with an equivalent surface mass is shown in Figures 4 and5. Calculation was limited to 1000 Hz maximum. This comparison shows that laminated glass has a much lower acoustic radiation compared to conventional glass at resonance frequencies due to the effect of the viscoelastic layer. The reduction of sound radiation power is around 5 dB in lower frequencies and around 10 dB in higher frequencies. For each mode, a frequency shift is observed (for example, the frequency of the first coupled mode is about 33.22 Hz without PVB and about 38 Hz with the PVB interlayer) .

In fact, at low frequencies, the viscoelastic material is soft and the damping is small. At higher frequencies, the stiffness decreases rapidly and the damping is highest. Moreover, flexural vibrations cause shear strain in the viscoelastic core which dissipates energy and reduces vibration and noise radiation. Note that the thickness of the viscoelastic layer has a significant influence in terms of attenuation.

Influence of the shunted piezoelectric patches

In order to reduce the vibrations and the sound radiations around the first resonance of the sandwich plate, four very thin identical piezoelectric patches (length 75 cm, width 5 cm, thick 0.5 mm ) are mounted on the plate (free of any mechanical charge) and tuned to an RL shunt circuit (see figure 3). The patch is composed of piezoelectric Macro Fiber Composite (MFC). The material properties of the MFC are: E 1 =30 GPa, E 2 =15.5 GPa, E 3 =15.5 GPa, G 13 =10.7 GPa, G 23 =10.7 GPa, G 12 =5.7 GPa, ν 13 =0.4, ν 12 =0.35, d 33 = 4.18 × 10 -10 m/V, d 32 = d 31 = -1.98 × 10 -10 m/V, 33 / 0 =1558,6 ( 0 =8.854 10 -12 F/m is the permittivity of free space), ρ =4700 kg/m 3 [START_REF]Smart Material Corp. Macro-Fibre Composites Technical Data, www.smart-material[END_REF][START_REF] Trindade | On the finite element characterization of d31 macro-fibre composites effective properties[END_REF]. Note that the e ij piezoelectric constants (in Coulomb/m2) can be deduced from the d ij piezoelectric constants (expressed in m/V or Coulomb/Newton) using the elastic compliances at constant electric field [START_REF] Rogacheva | The theory of piezoelectric shells and plates[END_REF].

The resistance R and the inductance L can be adjusted and properly chosen to maximize the damping effect of this particular mode. The optimal resistance and inductance of the ith mode for a series resonant shunt are given by [START_REF] Hagood | Damping of structural vibrations with piezoelectric materials and passive electrical network[END_REF]:

R opt = 2k 2 eff,i Cω i (1 + k 2 eff,i ) (45a)
L opt = 1 Cω 2 i (1 + k 2 eff,i ) (45b)
where ω i is the short circuit natural frequency of the ith mode, C is the capacitance of the piezoelectric patch and k eff,i is the effective electromechanical modal coupling factor (EEMCF), characterizing the energy exchanges between the mechanical structure and the piezoelectric patches and defined by:

k 2 ef f,i = ω2 i -ω 2 i ω 2 i ( 46
)
and where ω i is the open-circuited natural frequency. For the chosen mode, the optimal values of the shunt electrical parameters are then, R=429 Ω and L=83 H. Note that the portions of the plate covered by the piezoelectric patches have been modeled according to the first-order shear deformation laminated theory and only one electrical degree of freedom is used to represent the electrical charge Q in each patch. Figures 6 present the mean quadratic normal velocity and the sound power radiated by the plate in the external acoustic domain with and without piezoelectric shunt. The system response is obtained with the proposed FEM/BEM reduced order approach. The unknown mechanical displacements are projected on the basis formed by the first 20 real short-circuited modes and the first 20 real open-circuited modes. These modes are associated to the in-vacuo sandwich plate with viscoelastic interlayer and equipped with the piezoelectric patches. They are calculated from Equations ( 15) and (17) using the real and constant shear storage modulus G ∞ . The unknown pressures are projected on the basis formed by the first 20 acoustic modes of the internal fluid in rigid cavity calculated from Equation [START_REF] Larbi | Vibration reduction of an elastic structure using shunted piezoelectric systems: New reduced order finite element models[END_REF]. It can be seen that the resonant magnitude of the first mode has been significantly reduced (around 18 dB). In fact, the strain energy contained in the piezoelectric material is converted into electrical energy and hence dissipated into heat using the RL shunt device.

Conclusions

In this work, a coupled finite element/boundary element formulation (FEM/BEM) used to model the reduction of noise and vibrations by viscoelastic materials and piezoelectric techniques is presented. The multiphysics studied system is constituted of a sandwich structure with a viscoelastic core, equipped with shunted piezoelectric patches and coupled to interior and exterior acoustic domains. A reduced-order FE model, based on normal mode expansion, is developed for modeling the internal fluid-structure coupled problem. The proposed methodology requires the computation of the uncoupled eigenmodes of the undamped structure with short-circuited and open-circuited conditions and the acoustic eigenmodes of the fluid in a rigid cavity. Despite its reduced size, this model is proved to be very efficient for simulations of dynamic analysis of structural-acoustic coupled systems with viscoelastic interlayers and with shunt damping when appropriate damping terms are inserted into the modal equations of motion. Further investigations will concern the introduction of passive dissipation into the fluid and experimental aspects.

Figure 1 :

 1 Figure 1: Studied coupled problem.

Figure 2 :

 2 Figure 2: Complex shear modulus of PVB determined with a fractional derivative model.

18

 18

Figure 3 :

 3 Figure 3: Sandwich structure equipped with piezoelectric shunted patches and coupled to internal/external acoustic domains.

Figure 4 :

 4 Figure 4: The mean quadratic normal velocity and the sound pressure level at the cavity center with and without viscoelastic layer.

Figure 5 :

 5 Figure 5: Radiated sound power of the plate in the exterior acoustic domain with and without viscoelastic layer.

Figure 6 :

 6 Figure 6: The mean quadratic normal velocity and the radiated sound power of the plate with and without shunted piezoelectric patches.

Table 2 :

 2 Computed frequencies (Hz) of the structural-acoustic coupled system.