Violette Mulot

Hayat Benkhelifa

Didier PATHIER (c Fatou-Toutie Ndoye

Experimental and numerical characterization of food dehydration during freezing

Keywords: Model, food dehydration, weight loss, freezing

Using freezing in industry means understanding and controlling the process, particularly in terms of heat and water transfers. Water transfer is responsible for food dehydration during freezing which represents weight losses up to 6 % according to the product characteristics and the freezing process. A model of heat and mass transfer during the freezing process was developed. It takes into account not only heat transfer and crystallization but also water transfer, evaporation and sublimation. A methodology was established to identify some important and correlated model parameters: water activity, initial freezing temperature and ice fraction from the enthalpy-temperature curve of a product.

The model allows predicting the freezing time and weight loss for a given non-porous product (initial temperature, composition, geometry, dimensions, and physical properties) and a given freezing process (temperature, hygrometry, convective heat transfer coefficient). The numerical results were compared with experimental data obtained with Tylose in a pilot plant freezer from -30°C to -100°C.

INTRODUCTION

Freezing is an essential process in food preservation. It guarantees long preservation with relatively low impact on product quality. Freezing consists in reducing a product temperature below its freezing point in order to solidify a portion of its water. The water state change is the core of the freezing process because it decreases water activity which means water availability, thus, limiting the development of micro-organisms. Two types of transfers occur during the freezing of unwrapped products:

• Heat transfer: the product in contact with cold air releases energy to the surrounding environment. This leads to a temperature gradient in the product inducing a conductive heat transfer from the core to the surface and phase change (water crystallization) until reaching thermal equilibrium.

• Water transfer: the water vapor concentration in the surrounding air is lower than the water vapor concentration in the air in equilibrium with the product surface. Thus, a water transfer occurs between the product and its surroundings. As for heat transfer, a water concentration gradient appears in the product inducing a diffusive water transfer from the core to the surface and then to the surroundings where, during the pre-cooling stage, water evaporates from the surface. When ice is formed from the product surface, water diffusion is only significant in the non-frozen part of the product and ice sublimates from the surface causing the emergence of a dehydrated layer like in freeze-drying.

Water transfer induces food dehydration decreasing frozen food quality and weight. The weight loss depends on the product characteristics and the freezing process which represents both quality and financial losses.

Predictive models of freezing time and weight loss according to a product and a freezing process are useful tools. These models are developed from heat and mass transfer balances at unsteady state. These balances are coupled because of phase changes at the product surface: water evaporation and ice sublimation Analytical solving of transfer equations is very complex and supposes many simplifying assumptions. Indeed, heat transfer in the product is mainly conductive and the temperature gradient is the driving force, thereby Fourier's law is considered. Water transfer is essentially diffusive and the water concentration gradient in the product is the driving force. Fick's law is thus considered for water transfer.

Moreover, crystallization implies complexity in the solving of the heat transfer equation because of:

• A discontinuity zone at the interface between liquid and solid phase which leads to a sudden variation of food thermo-physical properties.

• A high amount of heat released due to the solidification of water (latent heat of solidification).

Numerical methods are usually chosen to solve the partial derivative equations of heat and water transfer during freezing. Many numerical models for prediction of freezing time and weight loss were developed, different from each other by the transfer hypotheses and the type of product. For discretization in time: either explicit [START_REF] Tocci | Numerical-models for the simulation of the simultaneous heat and masstransfer during food freezing and storage[END_REF], implicit [START_REF] Campañone | Finite-difference method to solve coupled heat and mass balances with simultaneous surface dehydration during freezing[END_REF] or Crank-Nicolson methods [START_REF] Campañone | Weight loss during freezing and storage of unpackaged foods[END_REF] were used. For discretization in space, various methods exist [START_REF] Pham | Modelling heat and mass transfer in frozen foods: a review[END_REF] : the finite difference method [START_REF] Campañone | Weight loss during freezing and storage of unpackaged foods[END_REF][START_REF] Campañone | Finite-difference method to solve coupled heat and mass balances with simultaneous surface dehydration during freezing[END_REF][START_REF] Tocci | Numerical-models for the simulation of the simultaneous heat and masstransfer during food freezing and storage[END_REF], the finite volume method [START_REF] Pham | Freezing time formulas for foods with low moisture content, low freezing point and for cryogenic freezing[END_REF] and the finite element method. [START_REF] Campañone | Weight loss during freezing and storage of unpackaged foods[END_REF] numerically solved transfer equations for two freezing stages considering 1D transfers, and food with regular shapes:

• The first freezing stage corresponds to the pre-cooling of the product surface from the initial temperature to the initial freezing temperature. During this stage, water diffuses from the product core to the surface where it evaporates.

• The second freezing stage occurs after the surface began freezing. A freezing front moves towards the thermal center. The product is divided into three zones: unfrozen, partially frozen and dehydrated. Water does not diffuse anymore through the product but ice starts sublimating from the surface to the surroundings creating a dehydrated layer with a sublimation front moving to the product thermal center. Some other numerical models were developed [START_REF] Andreasen | Modelling of Heat and Mass Transfer in Food Products[END_REF]Campañone et al., 2005a;[START_REF] Pham | Modelling heat and mass transfer in frozen foods: a review[END_REF][START_REF] Pham | Model for food desiccation in frozen storage[END_REF] with some differences in the assumptions. For example, [START_REF] Campañone | Weight loss during freezing and storage of unpackaged foods[END_REF] mentioned an ice fraction adsorbed by the product matrix that could not sublimate whereas [START_REF] Andreasen | Modelling of Heat and Mass Transfer in Food Products[END_REF] did not.

Once they are established, models are validated with experimental data of freezing time and weight loss during freezing. Meat and fish products were widely studied [START_REF] Boonsumrej | Effects of freezing and thawing on the quality changes of tiger shrimp (Penaeus monodon) frozen by air-blast and cryogenic freezing[END_REF][START_REF] Bustabad | Weight loss during freezing and the storage of frozen meat[END_REF][START_REF] Campanone | Monitoring of weight losses in meat products during freezing and frozen storage[END_REF][START_REF] Espinoza Rodezno | Cryogenic and air blast freezing techniques and their effect on the quality of catfish fillets[END_REF][START_REF] Lambrinos | The weight loss of frozen meat. Experimentation and prediction by using tylose[END_REF][START_REF] Pham | Model for food desiccation in frozen storage[END_REF]. Fewer studies focus on highly porous products like pre-baked bread (Hamdami et al., 2004a(Hamdami et al., , 2004b)).

From validated numerical models, simplified ones were also proposed [START_REF] Campañone | Weight loss during freezing and storage of unpackaged foods[END_REF]. All these models are decision tools for selecting one technology (characterized by the ranges of temperature and air velocity) relatively to another one in order to freeze a given product and limit the weight loss.

Numerical parametric analyzes permit to point out the influence of the freezing process parameters (temperature from -9°C to -45°C, flow velocity from 0.5 to 2.5 m s -1 and hygrometry from 40 to 52 %) on the weight loss and freezing time for a given product [START_REF] Campañone | Weight loss during freezing and storage of unpackaged foods[END_REF][START_REF] Campañone | Finite-difference method to solve coupled heat and mass balances with simultaneous surface dehydration during freezing[END_REF].

The aim of the present work is to develop a numerical model of heat and mass transfers in non-porous food during freezing. This model links three fundamental variables thanks to water phase diagram: water activity, initial freezing temperature and ice content. Only three parameters are necessary to predict all these variables and they are solely determined from the enthalpy curve Results of freezing time and weight loss obtained with the numerical model were compared with data from experiments with Tylose. Tylose has a constant composition and can be used as a beef meat model material for heat transfer although this is not really the case for mass transfer (meat membranes and fat layers are barriers for mass transfer which are not present in Tylose). These comparisons are presented for different freezing conditions not only in the mechanical cold domain (-30°C, -50°C) but also in the cryogenic freezing field (-100°C).

MATERIALS AND METHODS

1.1.

Experimental work

Freezing equipment

A freezing cabinet (Figure 1) was designed and developed for these experiments. It is able to recreate both mechanical and cryogenic freezing conditions in terms of temperature and gas flow velocity [START_REF] Mulot | Measurement of food dehydration during freezing: Development of an experimental tool and application to tylose[END_REF]. It consists of a closed cabinet where the cold flow is conducted through a 2.8 m long duct equipped with some perforated grids to homogenize it. Temperature is adjustable from room temperature to -100°C using a nitrogen injection and the gas flow velocity is adjustable from 0 to 9 m s - 1 thanks to a variable frequency fan. In the measuring zone, where samples are located, some preliminary measurements were made to validate the uniformity of temperature, gas flow velocity and convective heat transfer coefficient. A special device was used to estimate the convective heat transfer coefficients (h) at mid-height of the measuring zone. This device is composed of a copper plate (50×50 mm) which is heated by an electrical resistance and insulated underneath. Two type-T thermocouples measure the temperature of the copper surface and the gas temperature above it. A programmable controller performs temperature recording and flux calculation.

It permitted to establish the correlation between Nusselt (Nu), Reynolds (Re) and Prandtl (Pr) numbers characterizing the flow around plates in the measuring zone. Three flow velocities were studied: 4.2, 6.6 and 8.4 m s -1 for four temperature set points: 15°C, -30°C, -50°C and -100°C. Heat transfer coefficients were recorded for 120 s once the freezing cabinet was stabilized. From these results (Table 1), the convective heat transfer coefficient (h) as a function of temperature (T g ) and gas flow velocity (v s ) was deduced (Equation 1).

Nu = hL k ( ) = 0.066Re ( , ) . Pr ( ) (1)
where k g is the gas thermal conductivity and L the length of the sample (plate) placed in the freezing cabinet (L=80 mm).

Sample characteristics

Tylose plates were used (length=80 mm, width=50 mm, thickness=18 mm, mass=70.7-73.8 g). Their lateral faces were insulated with expanded polystyrene in order to consider only 1D heat and mass transfer (vertically, throughout the thickness of the product).

Sample initial temperature was uniform and was 20°C. Three samples were placed per freezing run in the freezing cabinet. One was dedicated to the sample core and surface temperature recording: two calibrated type-T thermocouples (precision 0.1°C) were inserted into a needle so that they could be accurately positioned at the sample core and just under the surface. The other two samples were dedicated to weight recording. Two scales were used: Mettler Toledo ® PG1003S and Sartorius ® QUINTIX3102-1S with a weighing precision of 0.01g. Scales were fitted with hangers for below-thebalance weighing. Thus, scales were kept outside the freezing cabinet while samples were in the cold zone.

Experimental procedure

Six freezing operating conditions were studied, three temperatures: -30°C, -50°C and -100°C and two gas flow velocities: 3.9 m s -1 and 7.7 m s -1 .

The freezing time was defined as the time required for thermal center of Tylose samples to reach -18°C.

Two types of experimental work were carried out:

• Firstly, Tylose dehydration due to freezing process was studied by measuring the total weight loss. Three replicates were made for each operating condition. Once the steady state was reached inside the freezing cabinet, process ventilation and nitrogen injection were switched off and samples were introduced in the freezing cabinet. The initial masses of samples suspended inside the cabinet were measured. Then, process ventilation and nitrogen injection were switched on until reaching the reference freezing time previously determined. At the end of freezing time, the freezing cabinet was switched off and the final masses of samples suspended inside the freezing cabinet were measured.

• Secondly, Tylose dehydration kinetics during freezing was studied by measuring weight loss at regular time steps for one freezing operating condition: -30°C/3.9 m s -1 (3 replicates). The protocol is similar to the previous one. However, four intermediate weighing were performed:

after 5, 10, 15, 30 and 45 minutes of freezing.

Modeling

A model was developed to predict heat and mass transfers in a non-porous food during freezing. This model assumes no change of density and size of food between the fresh and frozen state.

Water and heat transfer assumptions

Freezing can be considered like a succession of several periods.

The first step is the pre-cooling stage: the surface product temperature decreases from the initial temperature to the initial freezing temperature. Water is in liquid form and diffuses from the product core to the surface where it evaporates towards the outside environment. Heat is conducted throughout the product from the thermal center to the surface where it is partially released by convection in the environment and partially used for water evaporation. These heat and water transfer phenomena are represented in Figure 2a.

When product surface starts freezing, the second freezing step begins: water starts solidifying. Product is divided into three zones (Figure 2b):

• The first zone is the unfrozen zone (zone 1) where water has not begun freezing because temperature has not reached the initial freezing temperature. In this zone, water keeps diffusing from the core to the freezing front.

• The second zone from the freezing front to the sublimation front is the partially frozen zone (zone 2). The freezing front moves towards the thermal center of the product increasing the frozen part. It is assumed that water diffusion does not occur anymore in this zone because the ice crystal network hinders water migration.

• The third zone is a small dehydrated layer (zone 3) at the product surface due to ice sublimation. A sublimation front is slowly moving from the surface to the freezing front: ice crystals sublimation creates a porous dehydrated layer. In this dehydrated layer, water vapor diffuses through the air inside the pores until the product surface.

Water and heat transfer equations

Water transfer is essentially diffusive in the non-frozen part of the product, the driving force is the water concentration gradient in the product. Fick's law is considered for unidimensional water transfer (Equation 2).

∂C w ∂t = 1 x α ∂ ∂x x α D w ∂C w ∂x With C " = C #$ Y " (2)
where Cw is the water mass concentration, t is time, x is the location in the product, Dw is the water diffusivity and & is the geometric factor (& = 0 for an infinite plate, & = 1 for a cylinder and & = 2 in case of a sphere). Finally, Cdm is the dry matter concentration and Yw is the liquid water content in the product expressed on a dry basis. In the frozen part, it is considered that water does not diffuse anymore. In the dehydrated layer, water vapor diffuses from the sublimation front to the surface.

Heat transfer in the product is mainly conductive, the driving force is the temperature gradient.

Fourier's law for unidimensional heat transfer is used (Equation 3), considering that heat transfer by water diffusion is negligible.

ρCp *++ ∂T ∂t = 1 x α ∂ ∂x x α k ∂T ∂x With ρCp *++ = ρCp -C #$ L 1 dY 3 dT (3)
where ρ is the density, Cpapp is the apparent mass specific heat capacity, k is the product thermal conductivity, T is the product temperature, Cp is the mass specific heat capacity, Lf the latent heat of fusion and Yi the ice fraction in the product (dry basis).

Boundary conditions

Symmetry

As symmetrical shapes are considered for the products, the liquid water concentration gradient and the temperature gradient are zero at the symmetry axis or plan (Equation 4 and Equation 5).

∂C " ∂x =0 (4) ∂T ∂x =0 (5)

Surface and sublimation front

During the pre-cooling stage, product surface is not frozen yet, the difference between the water vapor concentration in equilibrium with product surface (Cw-veq) and in the surroundings (Cvs) causes a water evaporation flux (φevap) from the surface to the cold environment (Equation 6).

if T > 6if: -D " ∂C " ∂x = φ evap =h m :C w-veq -C ? @ (6)
where hm is the external mass transfer coefficient and Tif is the initial freezing temperature.

During the second period, when surface starts freezing, ice sublimates from the sublimation front creating a dehydrated layer of thickness e. The water vapor flux (φvap) resulting from the ice sublimation front diffuses throughout the gas contained in the pores of the dehydrated layer until the product surface and then, throughout the external boundary layer. The steady-state water vapor diffusion in the dehydrated layer is assumed: the variation of the water vapor mass in the dehydrated layer is neglected. φvap is calculated using Equation 7which involves the sum of the mass transfer resistances.

where Ybw is the fraction of bound water (dry basis), Ci-veq is the water vapor concentration in equilibrium with ice at the product surface, τ is the tortuosity of the dehydrated layer and ε the porosity. Dv is the water vapor diffusivity coefficient in the gas contained in the pores. It is assumed that partially frozen zone contains bound water, free water (Yw-Ybw) and ice whereas the dehydrated layer contains only bound water.

At the product surface, during the pre-cooling stage, a part of the heat conducted through the product is released by a convective transfer to the surroundings; the other part is consumed by water evaporation.

During the second freezing stage, the heat is conducted in the product from the thermal center to the sublimation front. At the sublimation front, a part of the heat is used by ice sublimation and unfrozen water evaporation; the other part is conducted through the dehydrated layer and then, it is released by convection to the surroundings. Equation 8 is used as boundary condition for the heat flux (φh-s), at the product surface for the pre-cooling stage (e=0 and de/dt=0) and at the sublimation front for the subcooling stage (φ evap = 0).

-k ∂T ∂x = φ h-s = T-T e k #A + 1 h +C dm :(Y w -Y bw )L v +Y i L s @ de dt +φ evap L v (8)
where h is the convective heat transfer coefficient, Tg is the temperature of the surroundings. Lv and Ls are respectively the latent heat of water evaporation and of ice sublimation and kdl is the thermal conductivity of the dehydrated layer.

Transfer parameters and Tylose properties

Parameters available in the literature

if T < T 31 : C #$ (Y " -Y C" + Y 3 ) de dt =φ *+ = C 3E FG -C ? eτ εD v + 1 h m (7)
To be able to solve equations introduced in part 1.2.2 and 1.2.3, several parameters need to be estimated. Table 2 summarizes these parameters and the selected methods for their estimation (all notations are described in the nomenclature).

The estimation of the thermal conductivity of the frozen product must be predicted accurately as the thermal conductivity of ice is significantly different to the conductivity of other components.

Therefore, three various models were tested in simulations: model 1 is and intermediate model presented by [START_REF] Gulati | Enabling computer-aided food process engineering: Property estimation equations for transport phenomena-based models[END_REF], model 2 was proposed by [START_REF] Pham | Effect of biot-number and freezing rate on accuracy of some food freezing time prediction methods[END_REF] to predict the thermal conductivity of Tylose gel according to the temperature and model 3 is a Maxwell model suggested by [START_REF] Cogné | Experimental data and modelling of thermal properties of ice creams[END_REF]. The last model takes into account the physical structure of food assuming that a component (ice) is dispersed in another one which is continuous (cryo-concentrated solution). A similar relation was proposed by [START_REF] Carson | Review of effective thermal conductivity models for foods[END_REF].

The tortuosity value of the dehydrated layer was also varied from 1.0 [START_REF] Campañone | Weight loss during freezing and storage of unpackaged foods[END_REF] to 2.0 in simulations. The tortuosity qualifies the complexity of the pore network and has an impact on the water vapor diffusivity throughout the dehydrated layer.

Water activity, initial freezing temperature and ice content correlation

Three other parameters used for the solving of heat and water transfer during freezing need to be determined: Cw-veq (or Ci-veq), Ybw, and Yi.

Cw-veq, which represents the water vapor concentration in equilibrium with liquid water at the product surface or sublimation location is calculated from Equation 9.

C "E FG = a " C vsat (9)
Cvsat is the saturation vapor concentration in equilibrium with pure liquid water, which depends on the product surface temperature. When the temperature is below 0°C, Cw-veq is still calculated using pure liquid water as reference state. Cvsat is then the vapor concentration in equilibrium with super-cooled pure liquid water as it is shown on Figure 3. The liquid phase of the product contains different solutes which makes that aw<1. When the product is partially frozen, the solute concentration in the liquid phase, which is directly related to the remaining liquid water content Yw, increases and aw decreases.

Therefore, the water activity is function of temperature and liquid water content and its variation is specific to each product. Nevertheless, a same expression can be used with parameters specific to each product. The expression of [START_REF] Iglesias | Prediction of the effect of temperature on water sorption isotherms of food material[END_REF] was chosen (Equation 10).

a " = exp L -exp(bT + c) (100Y " ) O P if T > 0°R (10) a " = exp -exp (c) (100Y " ) O if T < 0°R
b, c and r are constants, specific to a given product. To determine these parameters, a methodology was developed in this work, from equilibrium equation of pure water phase diagram (Figure 3) and from experimental data of enthalpy variation as a function of temperature. This methodology was applied to the case of Tylose. Parameter b was taken as -11.10 -3 °C-1 according to data found in literature [START_REF] Iglesias | Prediction of the effect of temperature on water sorption isotherms of food material[END_REF].

The saturation vapor pressure in equilibrium with pure liquid water (Psat.w(T)) and ice (Psat.i(T)) are expressed as a function of temperature (Figure 3) thanks to Clausius-Clapeyron expressions using the triple point T0= 0°C as reference (Equation 11 and Equation 12). In these equations, R is the universal gas constant and Mw is the water molar mass.

From Equation 12and the ideal gas equation, Cvsat is calculated with Equation 13.

C vsat = M " P ?*S." R(T + 273.15) (13)
As pure liquid water is used as the reference for the calculation of water activity, the initial freezing point is defined as when water activity for the unfrozen food equals the ratio of ice and water saturation vapor pressure (Figure 3) which is expressed by Equation 14.

a "(X YZ ) = P ?*S.3( [\ ) P ?*S."( [\ ) (14) Yw0 is the initial water content of the product on a dry basis.

Thanks to Equation 10, 11, 12 and 14 parameter r can be expressed as a function of c and Tif (Tif=-1.0°C for Tylose) with Equation 15.

r = ln -273.15 exp(c) R(T 31 + 273.15) L 1 M " T 31 ln(100Y " ) (15)
Below the initial freezing point (T<Tif), the rate of ice crystallization being fast relative to the rate of heat transfer throughout the product thickness, the food is locally considered in equilibrium. Thus, in the partially frozen zone, the remaining liquid food solution is in equilibrium with ice. For a given temperature, activity of the water in the liquid phase is equal of the activity of the water in the ice phase [START_REF] Fennema | Water activity at subfreezing temperatures[END_REF][START_REF] Storey | The equilibrium water vapour pressure of frozen cod[END_REF]. Thus, water activity is defined as the ratio of ice and water saturation vapor pressures and is independent of the composition of the food (Figure 3) which is expressed by Equation 16.

a " (X Y ) = P ?*S.3( ) P ?*S."( ) where Y " = Y " -Y 3 (16) Equation 16 combined with Equation 10, 11 and 12 leads to the following equation (Equation 17) which allows estimating the ice fraction on a dry basis Yi in the product as a function of temperature.

Y i =Y w0 - 1 100 _ L f M w exp(c)R 1 T+273.15 - 1 273.15 `-1 r (17)
Nevertheless, the bound water can not crystallize so that the maximal value of ice content corresponds to the initial free water content in the product (Equation 18).

Y i a Y w0 -Y bw (18)
At a given temperature, Tylose contains carbohydrates and ash (mass fractions Xch=23 % and Xa=0.6%

, specific heat capacity Cpch and Cpa), free and bound liquid water (mass fraction Xw=76.4 %, specific heat capacity Cpw) and ice (mass fraction Xi, specific heat capacity Cpi). The ice mass fraction is obtained from Equation 17. The specific heat capacity of Tylose at this temperature is then obtained as the weighted average of the heat capacity of the components (Equations 19).

Cp = b :X d Cp d @ d (19)
From Tylose mass specific heat capacity, Tylose enthalpy is then calculated as a function of temperature (Equation 20). The reference temperature was chosen at Tref=-40°C (Hcal(-40)=0).

H f*A( ) = g Cp ( ) dT -X 3( ) L 1 hi\ (20)
where Xi is the ice fraction in the product on a humid basis.

The calculated enthalpies (Hcal) were then compared with experimental ones (Hexp) [START_REF] Cleland | Assessment of Freezing Time Prediction Methods[END_REF] calculating the residues with Equation 21.

residues = g :H f*A( ) -H Fk+( ) @ l dt ml °n mEo °n (21)
By minimizing the residues, following parameter values were obtained for Tylose: c=4.0, Ybw=0.61 (Xbw=0.14), r=1.5 and b=-11.10 -3 °C-1 . Figure 4 shows the calculated enthalpies considering these parameter values in comparison with the experimental enthalpy curve.

Finally, for the numerical solving, the finite volume method was used for the spatial discretization. The mesh was refined near the product surface and a predictor corrector method was used for time integration. The problem was solved using the Matlab ® software.

Model without considering internal water transfer

The complete model presented above was compared to a simplified one which does not consider internal water transfer. The same equation (Equation 3) was used for heat transfer by conduction inside the product. At the product surface, the heat transfer by convection and the refrigeration effect of water evaporation or ice sublimation (assuming pure water) was expressed with Equation 22.

-k ∂T ∂x = h:T -T @ + φ vap L E?

with φ vap = h $ (C ?*S -C ? )

(22

)
Where Cvsat is the vapor concentration in equilibrium with pure water or ice and Lv-s is the latent heat of evaporation (T>Tif) or sublimation (T<Tif).

RESULTS AND DISCUSSION

Validation against experimental data

Experiments done with Tylose samples provided data on total relative weight loss and core temperature kinetics during the freezing process for 6 freezing operating conditions.

Simulations were also done with the same Tylose sample characteristics and the same 6 freezing process conditions. For experiments as well as for simulations, process was stopped when product core temperature reaches -18°C. Thanks to Equation 1, the experimental convective heat transfer coefficient h corresponding to each operating condition was calculated (L=80 mm). These values were used as input freezing process parameters for simulations.

Sensitivity analysis of the numerical model

Effect of tortuosity, thermal conductivity model and water diffusion coefficient

Table 3 presents, for the freezing condition -30°C/3.9 m s -1 , the numerical results for the freezing time and the relative weight loss according to the tortuosity value of the dehydrated layer (from 1.0 to 2.0)

and to the model used to estimate the thermal conductivity of the frozen product. As it is shown on Table 3, the thermal conductivity model for the frozen product has a relatively small impact on the total freezing time. Indeed, the Biot number for the frozen product is below 1 (around 0. 4 in this case) which means that the internal resistance of the Tylose plate is lower than the external resistance due to the small thickness of the plate. Moreover, model 1 estimates thermal conductivities only 8 % higher on average between -1.0°C and -20°C than model 2 and 3; model 2 and 3 are very close (1.3 % variation).

For tortuosity values, results indicate a large effect on the weight loss: increasing the tortuosity value significantly decreases the relative weight loss because the pore network is more complex and slows down the water vapor diffusivity throughout the dehydrated layer and so ice sublimation at the sublimation front. Thus, when tortuosity rises, the vaporization flux at the sublimation front lowers.

Since water vaporization is endothermic, it speeds up the cooling and so, if the sublimation flux decreases, the freezing time increases.

These tables allow selecting the reference value for the tortuosity and the reference model to estimate the thermal conductivity. The grey cases on Table 3 correspond to the selected values: the reference tortuosity is 1.75 and the conductivity reference model is the one propose by [START_REF] Pham | Effect of biot-number and freezing rate on accuracy of some food freezing time prediction methods[END_REF] Figure 5 presents the numerical and experimental temperature variation at product surface and core during freezing. Numerical curves were obtained using the reference tortuosity and the reference model of thermal conductivity. Figure 5 reveals that for the core temperature, the model slightly underestimates the pre-cooling stage whereas the crystallization and sub-cooling period are well estimated. When introducing samples in the freezing cabinet, ventilation and nitrogen injection were stopped during the first minute of freezing. This can explained the small offset between experiments and simulations. For the surface temperature, temperature could not be measured exactly at the surface.

The thermocouple was located approximately between 1 and 2 mm from the surface. The experimental temperature is close to the predicted one for a depth of 2 mm.

For the freezing condition -30°C/3.9 m s -1 , Figure 6 shows the experimental relative weight loss after 5, 10, 15, 30 and 45 minutes of freezing and the numerical prediction (τ.=1.75 and reference model for the thermal conductivity). The simulation is run until 45 min for this comparison. The numerical relative weight loss kinetics during freezing slightly overestimates the experimental one: +0.18 % after 5 min,+0.20 % after 10 min, +0.19 % after 15 min, +0.12 % after 30 min and +0.04 % at the end of freezing. However, the same trend is obtained: the weight loss rate is high at the beginning (precooling stage) and gradually slows down during the freezing and the sub-cooling stage.

A sensitivity analysis according to the water diffusivity coefficient, multiplying it by 2 or dividing it by 2 with respect to the reference model given in Table 2, showed that it has a small impact on the total weight loss (when core temperature reaches -18°C). The largest impact is for the slowest freezing condition: -30°C/3.9 m s -1 , where dividing it by 2 decreases the total weight loss by 1.9 % and multiplying it by 2 increases it by 3.7 %. Indeed, for this freezing condition, the precooling stage is the longest, nevertheless, the predicted water activity at the product surface remains close to one during this period. Moreover, the predicted water content at the product surface decreases almost in the same way and is 67.8 % at the end of the precooling using the reference diffusivity, 65.4 % dividing the diffusivity by 2 and 68.3 % multiplying it by 2. The impact of the water activity coefficient is becoming even less noticeable for faster freezing rate as the product starts freezing almost instantaneously (no more water diffusion in the frozen product).

Comparison with the simplified model (without considering internal water transfer)

Table 4 compares results of total relative weight losses and freezing times obtained experimentally with those predicted using the developed complete model (reference values for tortuosity and thermal conductivity model) and those predicted with the simplified model which does not consider internal water transfer (presented in Section 1.2.5).

For weight losses, results are very different. The weight loss is overestimated by the simplified model for slow freezing (-30°C and -50°C) in comparison with the experimental data and results predicted by the complete model. In this case, the resistance of water vapor transfer through the dehydrated layer cannot be neglected. For freezing at -100°C, the weight loss is still overestimated by the simplified model in comparison with the complete one, but, it is closer to experimental data. In this case, the resistance of the dehydrated layer (thickness or tortuosity) may be overestimated by the complete model.

For freezing times, the simplified model gives slightly lower values than the complete model at -30°C because it overestimates the cooling effect of vaporization. For low freezing temperatures (-50°C and -100°C), the values obtained with the simplified model are close to the prediction of the complete model because the dehydrated layer is finer.

In conclusion, a model including internal water transfer, notably the resistance of water vapor transfer through the dehydrated surface layer, is necessary to accurately predict the weight loss but does not improve significantly the prediction of freezing time.

Total relative weight loss and freezing time: experimental and numerical results

Figure 7 compares the experimental data of total relative weight losses (when core temperature reaches -18°C) and the numerical ones (τ =1.75 and reference model of thermal conductivity for the frozen product) according to the three freezing temperatures (-30°C, -50°C and -100°C) and the two gas flow velocities (3.9 m s -1 and 7.7 m s -1 ). The average absolute deviation between the numerical predictions and experimental values is 0.07. The difference is very low for the freezing conditions of -30°C and -50°C. Nevertheless, the gap is a little larger at -100°C. The models tends to underestimate the total weight loss. This can be explained both by the lack of accuracy of the experimental data for total relative weight losses and of model parameters at this low freezing temperature. Regarding the experimental measurements, the time during which ventilation and nitrogen injection are switched off (about 1 min), when samples are introduced and initial weighing is done, is relatively long compared to the freezing time at this temperature (14 min for 3.9 m s -1 and 10 min for 7.7 m s -1 ). This may introduce some measurement artifacts and explain the larger difference between experimental and numerical results at -100°C.

Figure 8 presents the experimental freezing times for the product core temperature to reach -18°C in comparison with numerical ones (τ =1.75 and reference model of thermal conductivity for the frozen product) according to the three freezing temperatures (-30°C, -50°C and -100°C) and the two gas flow velocities (3.9 m s -1 and 7.7 m s -1 ). It appears that the experimental freezing time is always longer than the numerical one. This can be explained by the fact that, experimentally, the first minute of freezing is done without ventilation and nitrogen injection. Moreover, the super-cooling is not taken into account in the model, it is assumed that the system is always at equilibrium. In addition, product expansion due to freezing which is about 5 % is not included to the model, this can also explain an underestimation of the freezing time. Indeed, the thermal resistance of the frozen product is therefore underestimated.

For the validation of the dehydrated layer thickness, some measurements were done on frozen sample images obtained by X-ray micro-tomography. Nevertheless, resolution of X-ray micro-tomography was not sufficient to analyze and accurately distinguish pores forming the dehydrated layer.

Influence of operating freezing conditions on predicted temperature and dehydration evolutions

Simulations were run (τ =1.75 and reference model of thermal conductivity for the frozen product)

until the product core temperature reach -18°C for the 6 studied operating conditions. Figure 9 shows the influence of freezing temperature and flow velocity on product core temperature, surface temperature, relative weight loss and dehydrated layer thickness variation over freezing time. On Figure 9a, a plateau can be observed for the core temperature when it reaches Tif (-1.0°C).

Both freezing temperature and flow velocity have a great influence on the product freezing kinetics and on the progress of the dehydration front. Decreasing the gas temperature or increasing the flow velocity (higher h) accelerates the heat transfer kinetics and thus the product freezing kinetics. It takes less time for the thermal center to reach -18°C and therefore time for water evaporation and ice sublimation is shorter. Moreover, for a lower gas temperature and a higher flow velocity, the product surface temperature decreases faster (Figure 9b): therefore, the water vapor concentration at the surface also decreases reducing the vapor flux from the product surface. Both phenomena mentioned previously lead to lower weight losses and thinner dehydrated layers (Figure 9c and9d). The less efficient freezing condition: -30°C/3.9 m s -1 corresponds to a standard mechanical freezing process, it leads to a relative weight loss of 1.62 % and a dehydrated layer thickness of 0.318 mm. The freezing condition -50°C/7.7 m s -1 corresponds to an efficient mechanical freezing process, it halves the relative weight loss (0.87 %)

and the dehydrated layer thickness is 0.130 mm. The most efficient freezing condition is -100°C/7.7 m s -1 which is a cryogenic freezing process, the relative weight loss is 0.28 %, it represents a third of the weight loss get for the efficient mechanical process and 5.8 times less than the less efficient mechanical process. The dehydrated layer thickness is only 0.039 mm.

Figure 10 presents the total relative weight loss (when the product core temperature is at -18°C) due to water evaporation and due to ice sublimation (and unfrozen water evaporation) according to the 6 freezing operating conditions. Figure 10 shows that the relative weight loss due to water evaporation (WLevap) during the first freezing stage is very low and negligible for freezing conditions -50°C/7.7 m s - 1 (WLevap=0.07 %); -100°C/3.9 (WLevap=0.05 %) and -100°C/7.7 m s -1 (WLevap=0.03 %) because as it is shown on Figure 9b, product surface freezes almost instantly being in contact with the very cold and highly ventilated environment.

CONCLUSIONS

A numerical model for predicting freezing time, weight loss during freezing and thickness of a dehydrated layer at the product surface for unwrapped non-porous food products was developed. A methodology was established to accurately determine and link expressions to calculate water activity, initial freezing temperature and ice fraction formed in the product. These correlations were set thanks to Clausius-Clapeyron phase equilibrium equations and thanks to experimental data related to the variation of the product enthalpy as a function of temperature.

Numerical results of relative weight loss were compared to experimental ones obtained with Tylose samples for 6 operating conditions in the mechanical and cryogenic freezing field. Results showed a good agreement between numerical and experimental data.

This numerical model allows studying the influence of the freezing operating conditions: freezing temperature and gas flow velocity on the product freezing kinetics, the relative weight loss and the thickness of the dehydrated layer due to ice sublimation at its surface.

This predictive tool is reliable and can be applied for various food materials both for mechanical and cryogenic freezing conditions.

NOMENCLATURE

  , r Parameters for water activity expression (b in °C-1 ) Cdm Dry matter concentration in the product (kg of dry matter m -3 of product) Cp Product mass specific heat capacity (J kg -1 °C-1

Figure 4 :Figure 5 :Figure 6 :Figure 8 :Figure 10 :

 456810 Figure 1: Operating diagram of the freezing cabinet-Top view

  Figure 1: Operating diagram of the freezing cabinet-Top view .............................................................. Figure 2: Heat and water transfers during pre-cooling stage (a) and freezing stage (b) (1: unfrozen zone; 2: partially frozen zone; 3: Dehydrated layer) ............................................................................... Figure 3: Pure water phase diagram ....................................................................................................... Figure 4: Calculated enthalpies and experimental enthalpy curve according to temperature .............. Figure 5: Numerical and experimental surface and core temperatures variations during freezing (operating condition -30°C/3.9 m s -1 ) ......................................................................................................Figure 6: Experimental and numerical relative weight loss variations over freezing time for the freezing condition -30°C/3.9 m s -1 ........................................................................................................... Figure 7: Experimental total relative weight losses against numerical ones according to freezing temperatures and for a flow velocity of 3.9 m s -1 (a) and 7.7 m s -1 (b) ................................................... Figure 8: Experimental freezing times for the core temperature to reach -18°C against numerical ones according to freezing temperatures and for a flow velocity of 3.9 m s -1 (a) and 7.7 m s -1 (b) ................ Figure 9: Core temperature (a), surface temperature (b), relative weight loss (c) and dehydrated layer thickness (d) according to freezing time for the 6 freezing operating conditions ................................ Figure 10: Predicted total relative weight loss (water evaporation + ice sublimation) according to freezing operating conditions for Tylose ...............................................................................................

						(kg m -2 s -1 )
	L	ρ ρ d	Length of the sample (plate) (m) Product density (kg m -3 ) Density of the component j (kg m -3 )	Lf	ρ g τ	-1 °C-1 ) Latent heat of fusion (J kg -1 ) Gas Density (kg m -3 ) Tortuosity of the dehydrated layer
	Lv-s	Latent heat of evaporation or sublimation (J kg -1 )	Ls	Latent heat of ice sublimation (J kg -1 )
	Lv	Latent heat of water evaporation (J kg -1 )	Mw	Water molar mass (kg mol -1 )
	Nu	Nusselt number	Pr	Prandtl number
	Psat.i Saturation pressure of water vapor in	Psat.w	Saturation pressure of water vapor in
			equilibrium with ice (Pa)			equilibrium with liquid water (Pa)
	R		Universal gas constant (8.314 J K -1 mol -1 )	Re	Reynolds number
	Sc	Schmidt number	t		Time (s)
	T		Product temperature (°C)	Tg	Temperature of the surrounding gas
						(°C)
	Tif	Initial freezing temperature (°C)	Tm	Mean temperature between the product
						surface and the
	Tref	Reference temperature (°C)	Ts		Temperature of the product surface (°C)
		v d	Volumetric fraction of the component j	Vpore	Volume of pores (m 3 )
	vs	Flow velocity in the surrounding gas (m s -1 )	Vtot	Total volume (m 3 )
	x		Location in the product (m)	Xi	Mass fraction of ice
	Xj	Mass fraction of the component j	Xbw	Mass fraction of bound water
	Yi	Ice content (kg of ice kg -1 of dry matter)	Yw	Liquid water content (kg of water kg -1 of dry matter)
	Ybw Bound water content in dry basis (kg of water kg -1 of dry matter)	Yw0	Initial water content of the product in dry basis (kg of water kg -1 of dry matter)
	α		Geometric factor (& = 0 for an infinite plate, sphere) & = 1 for a cylinder and & = 2 in case of a	ε #A	Porosity of the dehydrated layer
	μ	(Pa s) Gas dynamic viscosity	φevap	Water evaporation flux (kg m -2 s -1 )
	φh-s	Heat flow at product surface (W m -2 )	φvap	Vaporization (ice sublimation and
						unfrozen water evaporation) flux

Table 1 :

 1 Heat transfer coefficient (HTC) measurements ......................................................................... Table 2: Estimation of heat and mass transfer parameters .................................................................... Table 3 : Numerical freezing time (min) and relative weight loss (%) according to the tortuosity value and the model used to estimate the thermal conductivity of the frozen product for the freezing condition -30°C, 3.9 m s -1 ........................................................................................................................ Table 4: Total relative weight losses and freezing times obtained experimentally (Exp.), predicted with the developed model (Model) and predicted with the simplified model (Simp. Model) ...............

Table 1 : Heat transfer coefficient (HTC) measurements HTC (W m -² K -1 ) Temperature set point (°C)

 1

	Velocity set point (m s -1 )	15	-30	-50	-100
	4.2	61 ± 2	67 ± 6	69 ± 6	83 ± 6
	6.6	89 ± 4 93 ± 10 97 ± 10 123 ± 9
	8.4	108 ± 7 112 ± 14 120 ± 13 159 ± 13

Table 3 : Numerical freezing time (min) and relative weight loss (%) according to the tortuosity value and the model used to estimate the thermal conductivity of the frozen product for the freezing condition -30°C, 3.9 m s -1 Experimental freezing time = 43 ± 3 min

 3

		Tortuosity τ
	1.0	1.25	1.5	1.75	2.0

Table 4 : Total relative weight losses and freezing times obtained experimentally (Exp.), predicted with the developed model (Model) and predicted with the simplified model (Simp. Model)

 4

			Relative weight loss (%)	Freezing time (min)
	Temperature (Tg, °C)	Heat transfer coefficient (h, W m -2 K -1 )	Exp.	Model	Simp. model	Exp. Model	Simp. model
	-30	55 94	1.61 ± 0.05 1.56 ± 0.07	1.62 1.45	2.65 43 ± 3 2.81 29 ± 1	41.3 27.7	37.1 25.2
	-50	58 99	0.98 ± 0.11 0.85 ± 0.15	1.10 0.87	1.57 30 ± 2 1.32 21 ± 1	22.9 15.5	22.2 15.1
	-100	67 116	0.71 ± 0.10 0.46 ± 0.03	0.43 0.28	0.50 14 ± 1 0.34 11 ± 1	10.6 7.2	10.6 7.2

ACKNOWLEDGEMENTS

The authors gratefully acknowledge the financial support granted by the CIFRE department (Conventions Industrielles de Formation par la Recherche) from the French ANRT (Association Nationale Recherche et Technologie).

For T> T LM k= N k j


Diala Bitar
email: diala.bitar@entpe.fr

Emmanuel Gourdon
email: emmanuel.gourdon@entpe.fr

Claude-Henri Lamarque
email: claude.lamarque@entpe.fr

Manuel Collet
email: manuel.collet@ec-lyon.fr

C.-H Lamarque

Shunt loudspeaker using nonlinear energy sink

Keywords: Electroacoustic absorber, multiple time scales, nonlinear energy sink, nonlinear passive shunt

diffusion de documents scientifiques 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.

Introduction

Excessive use of transportation, machinery, sound systems and so on, in parallel with the high demands of sound comfort and health care, encouraged researchers to achieve the noise reduction challenge. Several procedures are used in order to reduce noise at the propagation and reception paths. For instance, classical passive systems as sound absorbing materials such as porous and adaptive liners have been found to be very efficient at high frequencies in controlling room noise in building technologies [START_REF] Arenas | Recent trends in porous sound-absorbing materials[END_REF]. However, these solutions are less efficient at mid and low frequencies in addition to their size. One of the best solutions that can be employed for these frequencies is to couple the treated system with linear resonating one, playing the role of an absorber.

As an example, we can mention the Dynamic Vibration Absorbers (DVA) such as Helmholtz resonators [START_REF] Helmholtz | Die Lehre von den Tonempfindungen als physiologische Grundlage für die Theorie der Musik[END_REF] and Frahm dampers [START_REF] Frahm | Device for damping vibrations of bodies[END_REF] which consist of a tuned spring-mass systems. Although these techniques can give perfect absorption, they can only treat very narrow frequency band [START_REF] Zilletti | Optimisation of dynamic vibration absorbers to minimise kinetic energy and maximise internal power dissipation[END_REF]. Later on, the concept of active noise absorption has been introduced based on controlling a secondary source in order to cancel the primary sound wave. This concept opened the way to a large variety of active sound absorption techniques.

We can mention feedback control based on sound pressure control of a loudspeaker [START_REF] Olson | Electronic sound absorber[END_REF][START_REF] Fleming | Control of resonant acoustic sound fields by electrical shunting of a loudspeaker[END_REF], the direct impedance control combining both sound pressure and diaphragm velocity sensing [START_REF] Bustamante | An adaptive controller for the active absorption of sound[END_REF] allowing a broad frequency control. In addition, it has been proposed to employ a pure electrical network as a second source to maximize the absorption efficiency under the concept of shunt loudspeaker [START_REF] Lissek | Electroacoustic absorbers: Bridging the gap between shunt loudspeakers and active sound absorption[END_REF]. For low-frequency noise absorption, a passive shunt electric circuit connected to the coil and powered by a DC source was employed by Zhang et al. [START_REF] Zhang | Thin broadband noise absorption through acoustic reactance control by electro-mechanical coupling without sensor[END_REF] allowing to modify the mechanical impedance of the coil, where Tao et al. [START_REF] Tao | Sound absorption of a finite micro-perforated panel backed by a shunted loudspeaker[END_REF] implemented a Micro-perforated panel backed by a shunted loudspeaker. Moreover, in a Rijke-tube an electroacoustic control device was experimentally and numerically tested as an alternative damping approach for suppressing the thermoacoustic instability [START_REF] Zhang | Electroacoustic control of rijke tube instability[END_REF].

Recent studies have shown that the performances of DVAs can be improved by employing nonlinear strategies. In fact, it has been shown in [START_REF] Lee | Suppression aeroelastic instability using broadband passive targeted energy transfers, part 1: Theory[END_REF][START_REF] Gendelman | Enhanced passive targeted energy transfer in strongly nonlinear mechanical oscillators[END_REF][START_REF] Sapsis | Effective stiffening and damping enhancement of structures with strongly nonlinear local attachments[END_REF][START_REF] Benarous | Nonlinear energy sink with combined nonlinearities: Enhanced mitigation of vibrations and amplitude locking phenomenon[END_REF] that employing nonlinear absorbers enables to absorb more efficiently the sound or vibrations at higher frequency band. This conceptual approach is based on irreversible energy transfer from a primary linear structure and nonlinear oscillator with small mass. This phenomenon is called targeted energy transfer or energy pumping, where the Nonlinear Energy Sink (NES) may serve as a nonlinear absorber [START_REF] Gourdon | Contribution to efficiency of irreversible passive energy pumping with a strong nonlinear attachment[END_REF][START_REF] Vakakis | Nonlinear Targeted Energy Transfer in Mechanical and Structural Systems[END_REF]. This procedure has achieved success in the mechanical field, where nonlinear phenomena were exploited to reduce vibrations [START_REF] Gourdon | Nonlinear energy pumping under transient forcing with strongly nonlinear coupling: Theoretical and experimental results[END_REF]. This concept was applied in the field of acoustics, providing a new technique of passive absorption for low frequencies [START_REF] Bellet | Experimental study of targeted energy transfer from an acoustic system to a nonlinear membrane absorber[END_REF].

However, this previous work [START_REF] Bellet | Experimental study of targeted energy transfer from an acoustic system to a nonlinear membrane absorber[END_REF] deals with very low frequencies with an absorber that is not essentially nonlinear. For this purpose, we suggest to switch a loudspeaker from its primary quality of diffusing sound and convert it into an electroacoustic absorber as a technical solution to achieve the desired nonlinear acoustical behavior. The general concept is to control the acoustical impedance a membrane, which can be achieved by tuning the dynamic of the loudspeaker membrane by connecting an appropriate electrical load. This approach allows reaching almost prefect absorption within a narrow frequency bandwidth around the resonance by employing a simple optimal passive electric resistance [START_REF] Lissek | Electroacoustic absorbers: Bridging the gap between shunt loudspeakers and active sound absorption[END_REF]. In addition, it has been shown that the bandwidth of control can be increased by choosing appropriate feedback gains in a combined pressure-velocity-feedback or by coupling the primary loudspeaker to a combination of several resonators [START_REF] Rivet | Multi-degree-of-freedom lowfrequency electroacoustic absorbers through coupled resonators[END_REF]. Hence, the goal of the present work consists in performing a purely acoustic nonlinear absorber acting on low frequency domain, by connecting to the loudspeaker terminals an electric nonlinear circuit that can be controlled. To our knowledge, it is the first attempt to apply the concept of targeted energy transfer to an electroacoustic absorber. Hence, this study aims at first to test the feasibility of the concept, where the non-optimal nonlinear results were compared to the ones obtained via an optimal linear resistor.

The present work is organized as follows: In section 2, we describe the electroacoustic absorber consisting of loudspeaker shunted with an electrical circuit. In section 3, an electroacoustic absorber consisting of a loudspeaker coupled with a nonlinear electrical shunt circuit at its terminals has been considered. The equivalent model is described by a 2 dofs system, including the main linear oscillator characterizing the dynamics of the loudspeaker membrane, which is linearly coupled to a cubic NES. In section 4, the Slow Invariant Manifold (SIM) approach has been employed. The differential system has been treated analytically for the case of 1:1 resonance at various time scales [START_REF] Gendelman | Targeted energy transfer in systems with external and self-excitation[END_REF]. The present methodology allows to identify the equilibrium points and fold singularities at the first slow time scale. In addition, it provides a predictive tool to design an optimize nonlinear electric circuit according to the energy exchange process between both systems. Section 5 provides an analytical expression allowing the identification of the different mechanisms of energy exchange in the case of free oscillations. In section 6, a physical application is considered presenting the advantage of combining a nonlinear electric shunt circuit to the loudspeaker terminals. Several analytical and numerical results are given in the case of free and forced vibrations, allowing a global understanding of the presented approach and the possible behaviors that can be obtained. Finally, conclusions and perspectives are given in section 7.

Electroacoustic absorber description

At the beginning, we consider an electroacoustic absorber (see Figure 1) following the same concept of Lissek et al. in [START_REF] Lissek | Electroacoustic absorbers: Bridging the gap between shunt loudspeakers and active sound absorption[END_REF][START_REF] Lissek | Optimization of electric shunt resonant circuits for electroacoustic absorbers[END_REF]. The physical model provides insight into the loudspeaker parameters that can be affected by electrical means. The mechanical dynamics of a loudspeaker diaphragm for small amplitudes and below the first modal frequency can be modeled by the following differential equation derived from Newton's second law

SP (t) = M ms vs (t) + R ms v s (t) + 1 C ms 1 C ms + ρc 2 S 2 V b v s (t)dt -Bl i(t). (1)
M ms , R ms and C ms are respectively the mass, the mechanical resistance92 and the compliance of the moving bodies of the loudspeaker. C ms is dened93 is the outgoing diaphragm velocity, Bl is the force factor of the moving-coil transducer with B representing the magnetic field magnitude and l the length of the wire in the voice coil. i(t) is the driving current, Bl i(t) the Laplace force induced by the current circulating through the coil. P (t) is the overall external sound pressure acting at the outer surface of the loudspeaker in P a.

so that 1/C mc = 1/C ms + ρc 2 S 2 /V b is
In the following, we will denote

C mc = (1/C ms + ρc 2 S 2 /V b ) -

Electrical nonlinear shunt resistor

The use of nonlinearity in the form of an electrical circuit permits the existence of quasiperiodic regimes over a range of external forcing terms, which provide efficient vibration suppression [START_REF] Gendelman | Quasiperiodic energy pumping in coupled oscillators under periodic forcing[END_REF]. Based on this assumption, and in order to design a nonlinear shunt resistor, we propose to connect an electric circuit enabling to obtain a pure cubic nonlinearity. As mentioned in [START_REF] Gendelman | Quasiperiodic energy pumping in coupled oscillators under periodic forcing[END_REF] it is possible to obtain electronically a cubic nonlinearity using multipliers. Based on this concept, we propose to connect the following example of a nonlinear shunt circuit (see Figure 3, where the crossed boxes represent quadratic multipliers ) to the transducer terminals.

One must mention that since we are connecting an unstable circuit, this kind of system becomes complicated to stabilize electronically. This is just an example of implementation; however, the practical realization may depend on weighting terms that will be considered as control parameters. The electrical nonlinear shunt circuit includes an inductor L c , a resistor R c playing the role of the NES linear damping coefficient and a capacitance C providing a pure cubic stiffness due the use of the multipliers. Then, the electrical voltage e(t) can be expressed as follows

e(t) = -L c di dt -R c i + V c -V c -kV 3 c ( 3
)
where

k = k 1 k 2 (r 1 + r 2 )(r 3 + r 4 ) r 1 r 2
is the nonlinear coefficient related to the multipliers connections with 1kΩ ≤ r 1 , r 3 and r 2 , r 4 ≤ 100kΩ. The ratios k 1 and k 2 are homogeneous to the inverse of voltage (V -1 ). Referring to the principle of operational amplifiers, the current I C , through the capacitor C is given by

I C = C dV c dt (4)
Then, the driving current i can be expressed as follows

i(t) = -(-I C ) = C dV c dt (5)
After replacing Equations [START_REF] Frahm | Device for damping vibrations of bodies[END_REF][START_REF] Olson | Electronic sound absorber[END_REF] into Equations (1,2), the dynamics of the loudspeaker can be described by the following system of differential equations coupling the electrical and the mechanical parts

       M ms ẍ(t) + R ms ẋ(t) + 1 C mc x(t) -CBl dV c dt = SP (t) C(L e + L c ) d 2 V c dt 2 + C(R e + R c ) dV c dt + kV 3 c + Bl ẋ(t) = 0 (6)
Considering the fact that the front of the diaphragm is subjected to a From a global point of view, the circuit is not passive in the way that very few tenth of watt are used to power the amplifiers in the electronics.

However, if we consider the electromechanical part on its own, the system is considered passive (negative energy dissipation) or semi-passive (the average of the exchanged power is negative) enabling to determine the stability of the system. In the mechanical part, the pumped energy is dissipated by joule heating effect: we implement to dissipate. The aim of this paper is to show the effects of a new (passive) nonlinear control that is not necessarily optimized.

Analytical developments

An analytical treatment suitable for this kind of systems is detailed below.

It allows to detect time multi-scale energy pumping between the primary system that describes the displacement of the loudspeaker and the shunt nonlinear circuit. This procedure can be divided into several steps: re-scaling, complexification of the system and keeping first harmonics and embedding time into fast and slow time scales. This schema permits the detection of the invariant manifold of the system at fast time scale, in addition to the equilibrium and fold singularities identification of the obtained reduced order system at slow time scales.

Rescaled system

For convenience and equation simplicity, we introduce the following nondimensional variable T = ω 0 t with ω 0 = 1/M ms C mc representing the natural frequency of the main system (Loudspeaker). Then, we write the linear and nonlinear coefficients with the following corresponding scaling with respect to realistic physical parameters where

     d dt = d dT dT dt = ω 0 d dT d 2 dt 2 = ω 2 0 d 2 dT 2 (7)
Then the system takes the following form

       d 2 x dT 2 + R ms C mc ω 0 dx dT + x(T ) -CBlC mc ω 0 dV c dT = SA m C mc cos ( ω ω 0 T ) Cω 2 0 (L e + L c ) d 2 V c dT 2 + Cω 0 (R e + R c ) dV c dT + kV 3 c + Blω 0 dx dT = 0 (8)
In what follows, we will use the complexification method of Manevitch and the multiple scale approach to detect the system invariant manifold at the fast time scale, and the system behavior at slower time scales. To do so, one must write the second equation of System [START_REF] Lissek | Electroacoustic absorbers: Bridging the gap between shunt loudspeakers and active sound absorption[END_REF] in a manner to have very small mass of the NES compared to the loudspeaker. According to the system physical parameters, the sum of inductance (L e + L c ) is smaller than 1 H. Then, ε = (L e + L c )/L 0 with L 0 = 1 H can be considered as our small dimensionless parameter; describing the ratio of mass between the loudspeaker and the NES. Then, dividing the second equation by Cω 2 0 L 0 , System (8) can be written as

       d 2 x dT 2 + R ms C mc ω 0 dx dT + x(T ) -CBlC mc ω 0 dV c dT = SA m C mc cos ( ω ω 0 T ) ε d 2 V c dT 2 + (R e + R c ) ω 0 L 0 dV c dT + k Cω 2 0 L 0 V 3 c + Bl Cω 0 L 0 dx dT = 0 (9)
The loudspeaker damping coefficient R ms C mc ω 0 is smaller than one, so it can be expressed as ελ. In addition, the NES damping coefficient (R e + R c )/(ω 0 L 0 ) can be expressed as εγ by choosing an appropriate control resistor R c . Moreover, the coupling terms CBlC mc ω 0 and Bl/(Cω 0 L 0 ) can be written as εα and εη respectively due primarily to the value of Bl and the chosen capacitance C. The value of the nonlinear cubic coefficient k/(Cω 2 0 L 0 ) can be adjusted by both k and C and then be expressed as εξ. On the other hand, the external forcing term SA m C mc is of order 1, so it can be represented by εf Then, the scaled differential system (9) becomes

       d 2 x dT 2 + ελ dx dT + x(T ) -εα dV c dT = εf cos ( ω ω 0 T ) ε d 2 V c dT 2 + γ dV c dT + ξV 3 c + η dx dT = 0 (10)
In order to examine the damped nonlinear normal modes of system [START_REF] Tao | Sound absorption of a finite micro-perforated panel backed by a shunted loudspeaker[END_REF] in the vicinity of 1:1 resonance, both variables x and V c are supposed to have frequency close to Ω = ω/ω 0 . The complex variables of Manevitch are introduced [START_REF] Manevitch | The description of localized normal modes in a chain of nonlinear coupled oscillators using complex variables[END_REF], which enable us to apply the truncated Fourier series later on (i = √ -1):

     dx dT + iΩx = ϕ 1 e iΩT dV c dT + iΩV c = ϕ 2 e iΩT (11)

Complexification-averaging procedure

If the damped nonlinear normal modes for the case of 1:1 resonance is considered, then ϕ 1 and ϕ 2 which defines the coefficients of the Fourier series expansion of x and V c may be considered as functions of slow time scales ε j T (j ≥ 1). In order to preserve only first harmonic e iΩT for each equation of System [START_REF] Tao | Sound absorption of a finite micro-perforated panel backed by a shunted loudspeaker[END_REF], we apply the Galerkin method and a truncated Fourier series.

Then we obtain the following system

     φ1 -i (1 -Ω 2 ) 2Ω ϕ 1 + ε λ 2 ϕ 1 -ε α 2 ϕ 2 = ε 1 2 f ε φ2 + i Ω 2 ϕ 2 + γ 2 ϕ 2 -i ξ 2Ω ϕ 2 G(|ϕ 2 | 2 ) + η 2 ϕ 1 = 0 (12) with G(|ϕ 2 | 2 ) = 3 4Ω 2 |ϕ 2 | 2
representing the resonant term.

Multiple-scale analysis

System (12) may be analyzed by using a multiple time scales approach, where we introduce the slow times τ 1 , τ 2 , . . . and the fast time τ 0 as

τ 0 = T τ j = ε j T j = 1, 2, . . . (13)
so that

d dT = ∂ ∂τ 0 + ε ∂ ∂τ 1 + ε 2 ∂ ∂τ 2 + . . . (14)
We will analyze the system behavior around 1:1 resonance by using a multiple scales approach assuming that

Ω = 1 + εσ (15)
with σ a detuning parameter. At the order ε 0 , resonant terms at τ 0 scale in system (12) yields

∂ϕ 1 ∂τ 0 = 0 ⇒ ϕ 1 = ϕ 1 (τ 1 , τ 2 , . . . ) (16)
The previous equation shows that the function ϕ 1 is constant during the time τ 0 = T , which means that in the principal approximation ϕ 1 depends on the slow times τ 1 , τ 2 , . . . However, at the order ε 1 the first equation of system [START_REF] Lee | Suppression aeroelastic instability using broadband passive targeted energy transfers, part 1: Theory[END_REF] gives

       ∂ϕ 1 ∂τ 1 + iσϕ 1 + λ 2 ϕ 1 - α 2 ϕ 2 = 1 2 f ∂ϕ 2 ∂τ 0 + i 1 2 ϕ 2 + γ 2 ϕ 2 + η 2 ϕ 1 -i ξ 2 ϕ 2 G(|ϕ 2 | 2 ) = 0 (17) (18)
Fixed points corresponds to the behavior of ϕ 2 for τ 0 → +∞ which means ∂ϕ 2 ∂τ 0 = 0 can be obtained from Equation ( 18):

γ + i(1 -ξG(|ϕ 2 | 2 )) ϕ 2 = -ηϕ 1 . (19)
Equation ( 19) presents an asymptotic equilibrium governed by a manifold called Slow Invariant Manifold (SIM). By writing ϕ 1 and ϕ 2 into their polar form as ϕ 1 = N 1 e iδ 1 and ϕ 2 = N 2 e iδ 2 , Equation ( 19) can be written and reduced into the following form:

ηN 1 e i(δ 1 -δ 2 ) = -iN 2 -γN 2 + iξN 2 G(N 2 2 ). ( 20
)
Taking into account the definitions of the complex variables ϕ 1 and ϕ 2 , their modules N 1 and N 2 can be expressed in terms of the initial variables x and V c of the physical system as

       N 1 = |ϕ 1 | = dx dT 2 + Ω 2 x 2 N 2 = |ϕ 2 | = dV c dT 2 + Ω 2 V 2 c ( 21
)
After separating the real and imaginary parts of Equation [START_REF] Rivet | Multi-degree-of-freedom lowfrequency electroacoustic absorbers through coupled resonators[END_REF], one can reach the following system

γN 2 = -ηN 1 cos(δ 1 -δ 2 ), (1 -ξG(N 2 2 ))N 2 = -ηN 1 sin(δ 1 -δ 2 ). (22) (23)
Taking the square of both sides of equations ( 22) and ( 23) and summing them, we obtain

N 1 = N 2 η γ 2 + 1 -ξG(N 2 2 ) 2 . ( 24
)
Local extrema of Equation ( 24) can be revealed using the following criteria

∂N 2 1 ∂N 2 2 = 0, ( 25
)
which finally yields the equation

γ 2 + (1 -ξG(N 2 2 ))(1 -ξG(N 2 2 ) -2ξN 2 2 G (N 2 2 )) = 0. ( 26
)
Equation ( 26) can be verified when the damping coefficient γ associated to the nonlinear equation is smaller than a critical value γ c which is equal to 1/ √ 3 [START_REF] Gendelman | Quasiperiodic energy pumping in coupled oscillators under periodic forcing[END_REF]. In addition, depending on the value of γ the system can have one, two or three fixed points. Moreover, for γ < γ c there exists two local extrema called N 21 and N 22 that can be written as follows

N 21 = 2 3 2 ξ - 1 -3γ 2 ξ , N 22 = 2 3 2 ξ + 1 -3γ 2 ξ . ( 27
)
It has been proved that the stable area of the fixed points in the reduced linearized system at the τ 0 time scale can be defined by the following criterion

[25]

γ 2 + (1 -ξG(N 2 2 ))(1 -ξG(N 2 2 ) -2ξN 2 2 G (N 2 2 )) ≥ 0 (28)
According to the previous inequality, stable fixed points corresponds to

N 2 < N 21 and N 2 > N 22
. By substituting the stable solution ϕ 2 with respect to the fast time scale τ 0 generated in Equation [START_REF] Bellet | Experimental study of targeted energy transfer from an acoustic system to a nonlinear membrane absorber[END_REF] into Equation ( 17), writing ϕ 2 into its polar form as ϕ 2 = N 2 (τ 1 )e iδ 2 (τ 1 ) and separating the real and imaginary parts of the resulting equation one finally obtains the following system

             ∂N 2 ∂τ 1 = BF -DE AD -CB = f 1 (N 2 , δ 2 ) g(N 2 ) , ∂δ 2 ∂τ 1 = CE -AF AD -CB = f 2 (N 2 , δ 2 ) g(N 2 ) . (29)
Explicit values of A, B, C, D, E and F are given in Appendix A. The approximation of the fast time scale τ 0 , enables Equations [START_REF] Vakakis | Nonlinear Targeted Energy Transfer in Mechanical and Structural Systems[END_REF][START_REF] Gourdon | Nonlinear energy pumping under transient forcing with strongly nonlinear coupling: Theoretical and experimental results[END_REF] to illustrate the response of the system regarding the initial conditions, when the system approaches the resonance manifold. In addition, the approximation of N 2 and δ 2 represented by System (29), describes the slow-time τ 1 evolution of the system at the manifold.

Equilibrium points of the system are defined by

f 1 (N 2 , δ 2 ) = f 2 (N 2 , δ 2 ) =
0 and g(N 2 ) = 0 while fold singularities are given by

f 1 (N 2 , δ 2 ) = f 2 (N 2 , δ 2 ) = g(N 2 ) = 0.
We should note that the stability of g(N 2 ) can be treated as same as the invariant manifold (see Equation ( 24)). In fact, g(N 2 ) = 0 provides both extrema N 21 and N 22 which are called the fold lines of the system.

Observation of free oscillations

In order to illustrate the mechanism of Targeted Energy Transfer (TET) between the primary system and the NES according to the slow time τ 1 , we start by setting Z 1 = N 2 1 and Z 2 = N 2 2 . Hence, in the case of free oscillations and by using the first equation of System (29), ∂Z 2 /∂τ 1 can be written in the following form

∂Z 2 ∂τ 1 = - αγη + λ γ 2 + 1 -ξG(Z 2 ) 2 γ 2 + 1 -ξG(Z 2 ) 1 -2Z 2 ξG (Z 2 ) -ξG(Z 2 ) Z 2 . ( 30
)
The integral of the previous equation can be computed explicitly as below

R(τ 1 , Z 2 ) = Cte (31) with R(τ 1 , Z 2 ) = 2λ αη + γλ αγη + λ + γ 2 λ τ 1 + 2αη γλ arctan λ γ(αη + γλ) (1 - 3 4 ξZ 2 ) + αη + γλ 2(1 + γ 2 )λ log (Z 2 ) + 3αγη + 2(1 + γ 2 )λ log 16αγη + λ 16γ 2 + (4 -3ξZ 2 ) 2 (32)
and Cte is an integration constant that can be related to the initial conditions. In fact, Z 1 and Z 2 reflect the energy stored in the relative displacement of the loudspeaker membrane and the voltage across the capacitor in the nonlinear shunt circuit respectively.

Results and discussions

Validation of the approach

In order to validate the obtained analytical approaches, we use a direct numerical integration method. The function ode45 implemented in MATLAB c which is based on an explicit Runge-Kutta (4,5) formula is employed in order to solve the dimensional differential system [START_REF] Tao | Sound absorption of a finite micro-perforated panel backed by a shunted loudspeaker[END_REF] and the complex averaged one [START_REF] Lee | Suppression aeroelastic instability using broadband passive targeted energy transfers, part 1: Theory[END_REF].

We consider the scaled parameters: ε = 0.01, λ = 0.1, α = 0.1, γ = 0.35, η = 0.1 and ξ = 0.5 and the following initial conditions: x(0) = V c (0) = Vc (0) = 0 and ẋ(0) = 7. For free oscillations (f = 0), Figure 4 shows a comparison of N 1 and N 2 between numerical solutions obtained by direct integration of the initial differential system [START_REF] Tao | Sound absorption of a finite micro-perforated panel backed by a shunted loudspeaker[END_REF] and the averaged one [START_REF] Lee | Suppression aeroelastic instability using broadband passive targeted energy transfers, part 1: Theory[END_REF].

Curves shown in Figure 4 includes three various phenomena stages [START_REF] Weiss | Dynamical behavior of a mechanical system including saint-venant component coupled to a non-linear energy sink[END_REF]. The first one begins from the initial point of the simulation until almost τ 0 = 550, where the NES maximal amplitude decreases slowly and the main system decreases rapidly. In the second stage and around τ 0 = 550, the system reaches a local extrema during which N 2 decreases abruptly and N 1 is almost constant. Finally, during the third stage (τ 0 > 550) both amplitudes decreases almost linearly. Since the numerical results of the very first differential system (12) includes all harmonics, when we kept the first harmonic and truncated the highest ones for the averaged system (10), we can remark the good agreement between both results. At τ 0 time scales, fixed points of the system are represented by the Slow Invariant Manifold (SIM) shown in Figure 5. Full and dashed lines correspond to stable and unstable parts of the τ 0 invariant respectively with N 21

and N 22 defining the extrema according to the N 2 variable. In fact, the SIM corresponds to a geometrical representation of all types of system behaviors in the slow time scales regardless the forcing term. In the following sections, we will provide a detailed study about the different dynamics of the system during free and forced vibrations. For free oscillations, the analytical energy trends can be obtained by solving Equation (31) with a corresponding constant C, providing an illustration of the possible scenarios of TET according to the initial conditions [START_REF] Nguyen | Design criteria for optimally tuned nonlinear energy sinks-part 1: transient regime[END_REF]. Figure 6 shows that two stable possibilities of energy pumping for the cases of free oscillations can be reached. More precisely, for initial conditions corresponding to a fixed point on the lower branch of the SIM, the mechanism of energy shows that this case is not efficient for TET. However, for initial conditions corresponding to a fixed point on the upper branch of the SIM, Figure 7 demonstrate that the process of energy is favorable for optimal NES design. It can be concluded that in order to activate the NES for efficient TET, initial conditions must be chosen in a manner that N 1 (0) is greater that the maximum N 21 ordinate.

Application

In order to study the system behavior when coupling a nonlinear shunt resonator, we consider the physical parameters as Lissek et al. [START_REF] Lissek | Electroacoustic absorbers: Bridging the gap between shunt loudspeakers and active sound absorption[END_REF] from the electrodynamic moving coil to the enclosure volume which are listed in Table 1. Among the large set of physical parameters that can be considered in the nonlinear shunt circuit, additional criteria are required. For instance, according to our choice of scaling and since the NES damping constant γ must be smaller than 1/ √ 3 as detailed previously, the choice of the control resistor R c must verify the following inequality

0 < R c + R e < ε ω 0 √ 3 . ( 33
)
Initial conditions were considered in a manner that the loudspeaker membrane is subjected to an initial velocity ẋ(0) = 0. The rest of the initial conditions are assumed to be null: x(0) = V c (0) = Vc (0) = 0, which corresponds to N 2 (0) = 0 and N 2 (0) = N 0 = ẋ(0). N 0 is chosen in a manner to be above the ordinate of the SIM maximum (N 0 > N 11 ) to examine a situation where the NES is active to control oscillations. Then, we set the control parameters of the nonlinear electric shunt in Table 2.

Results for parameters of Test 1

In what follows, we will study the system behaviors for the physical parameters listed in Table 2 untitled Test 1 for and forced oscillations.

Free oscillations

Figure 8 shows N 1 vs N 2 , with a comparison between the analytical invariant manifold given in Equation (24) (dashed-dotted curve) and the numerical integration of the differential system (10) (full line). Remarkably, the behavior of the system amplitudes N 1 and N 2 tracks the SIM. More precisely, full curve undergoes an immediate jump and oscillates around the lower stable branch until its amplitude decreases quasi-linearly.

Figure 8: Invariant manifold of the system at the fast time scale τ 0 (black curve) and its corresponding numerical results (green curve) for the system under free oscillations P (t) = 0 and σ = 0. Numerical results correspond to the integration of System (10) versus t for t ∈ [0, 0.15] with the initial conditions x(0) = V c = Vc = 0 and ẋ(0) = 1.9 × 10 -3 m s -1 .

Figure 9: Histories of displacement in terms of time t(s) obtained by numerical integration of the system under free vibrations, σ = 0 with the initial conditions corresponding to x(0) = V c = Vc = 0 and ẋ(0) = 1.9 × 10 -3 m s -1 . Dashed and full lines denote the cases of open circuit and shunt nonlinear circuit respectively.

In fact, when choosing initial conditions that corresponds to a point above N 21 , energy flows immediately and encounters a sudden breakdown toward the stable lower branch. In other words, the energy in the loudspeaker membrane is irreversibly streamed toward the nonlinear shunt circuit until it converges to zero.

The corresponding velocity histories versus time are represented in full line in Figure 9. Dashed curve represents the velocity histories in terms of time for the case of open circuit, which means when the electroacoustic transducer is not coupled with any electric charge. Clearly, the NES permits reduction of oscillations amplitudes.

Primary system under external forcing

Unlike the previous investigations, the loudspeaker diaphragm is now subjected to a periodically varying sound pressure P (t). As detailed previously, at the slow time scale τ 1 , all achievable system dynamics can be obtained by solving [START_REF] Lee | The jump phenomenon effect on the sound absorption of a nonlinear panel absorber 36 and sound transmission loss of a nonlinear panel backed by a cavity[END_REF]. More precisely, equilibrium points which corresponds to periodic regimes are obtained when

f 1 (δ 2 , N 2 ) = f 2 (δ 2 , N 2 ) = 0 and g(N 2 ) = 0,
where fold singularities corresponds to quasi-periodic regimes are achievable

when f 1 (δ 2 , N 2 ) = f 2 (δ 2 , N 2 ) = g(N 2 ) = 0.
As the scaled forcing term f is proportional to the pressure amplitude A m which is defined as a function of the sound pressure level L p , this later is considered as a main parameter of the nonlinear study. For instance, considering σ = 0 which corresponds to the resonance frequency ω 0 , Figure 10 For L p = 120 dB and σ = 0, Figure 10 shows the existence of one equilibrium point represented in triangle, corresponding to the intersection of full line and dash-dotted curves labeled 2 far the fold singularities. |dV C /dT + iΩV c | (Green curve) obtained by numerical integration of System (10) compared to the invariant manifold. We can see that the system oscillates around the lower stable branch of the invariant until being attracted by the fixed point. The corresponding velocity time series is represented in Figure 12, where we compare the histories of ẋ(t) with and without coupling in addition to the time histories of Vc (t). The periodic regimes obtained numerically validate our analytical predictions illustrated by the equilibrium point in Figure 10. In fact, the energy was transmitted from the loudspeaker membrane under periodic sound pressure to the non-excited nonlinear shunt circuit, with the same frequency defining the 1:1 resonance.

Figure 11: System invariant manifold and its corresponding numerical result for L p = 120 dB and σ = 0. The numerical result corresponds to the integration of System [START_REF] Tao | Sound absorption of a finite micro-perforated panel backed by a shunted loudspeaker[END_REF] with the initial conditions corresponding to x(0) = V c = Vc = 0 and ẋ(0) = 1.9 × 10 -3 m s -1 . In this case, the system is attracted by the periodic regime that corresponds to the fixed point marked with a triangle in Figure 10. Figure 13: System invariant manifold and its corresponding numerical result for L p = 132 dB and σ = 0. The numerical result corresponds to the integration of System [START_REF] Tao | Sound absorption of a finite micro-perforated panel backed by a shunted loudspeaker[END_REF] with the initial conditions corresponding to x(0) = V c = Vc = 0 and ẋ(0) = 1.9 × 10 -3 m s -1 . Here, the system present a SMR illustrated in Figure 10 with two dots.

For L p = 132 dB and σ = 0, we can deduce from the dash-dotted and solid curves labeled 5 in Figure 10 that the system possesses one unstable equilibrium point marked with a stars between the fold lines and two fold singularities denoted by dots. Hence, the corresponding numerical result and the analytical invariant are shown in Figure 13. This curve refers to a Strongly Modulated Response (SMR) of the system, including jump phenomena between both fold lines. In fact, for similar behavior both the loudspeaker membrane and the nonlinear shunt circuit are beating, which is very favorable for passive energy harvesting control applications. Energy exchange process of the SMR is qualitatively well described by the jump phenomenon in Figure 13. Besides, the values of N 2 seem to be perfectly appropriate when the jumps occur. However, some quantitative errors appear for the values of N 1 due to several factors. Among them we can mention the choice of ε, the numerical inaccuracies besides the choice of the Manevitch complex variables which takes into account the first harmonic only.

Figure 14 shows the corresponding numerical results using a Runge-Kutta schema according to the time t ∈ [0, 0.6]. In addition, Figure 15 shows the Poincaré's section after eliminating the data corresponding to the transient regime of both primary structure and attached nonlinear circuit. These figures confirm the obtained quasi-periodic regime. In fact, an important amount of vibrational energy is rapidly transferred to the NES and damped out in a quasi-periodic regime. Thus, the attenuation of vibrational energy of the loudspeaker arise with good efficiency in the vicinity of the natural frequency ω 0 .

(a) (b) show the results for test 2 for for L p = 115 dB and for L p = 124 dB respectively. Upper figures represent the system manifold invariant and their corresponding numerical results with σ = 0. The numerical results correspond to the integration of System [START_REF] Tao | Sound absorption of a finite micro-perforated panel backed by a shunted loudspeaker[END_REF] with the initial conditions corresponding to x(0) = V c = Vc = 0 and ẋ(0) = 4 × 10 -5 m s -1 . The middle and lower figures show the comparison of velocity histories as a function of time t(s) obtained by numerical integration of the studied system without coupling (black dotted line) and with nonlinear coupling (red solid line). Here, the system present a periodic regime for (a) and a SMR illustrated for (b).

Acoustic admittance

In order to examine the influence of coupling a passive nonlinear shunt circuit to the transducers terminals, it is necessary to have an understanding of the sound wave propagation. The normalized acoustic admittance presented at the loudspeaker diaphragm provides information about the amplitude of the reflection for a given frequency, which can be written as follows

Y m (ω) = -ρc V (ω) P t (ω) , (34)
with V (w) and P t (w) are the Fourier transforms of v(t) and p t (t) respectively. The minus sign appears as V (w) is defined to describe the outgoing diaphragm velocity.

According to the previous analytical investigations detailed in Section 6.3.2 at the τ 1 time scale and for the case of 1:1 resonance, two different regimes can be distinguished. The periodic or the quasi-periodic regimes are reached according to the given values of the forcing term and frequency.

Then, in order to represent the obtained normalized admittance for a given sound pressure level we followed the process below.

For each ω in the frequency domain we integrate the differential system using an explicit Runge-Kutta (4,5) algorithm implemented in Matlab c from t 0 = 0 to t f = N T (N denoting the number of periods T = 2π/ω) with y 0 = [0, 0, 0, 0] as initial conditions.

Then, in the permanent regime, a study of Poincaré section has been applied. More precisely, we calculate the value of ∆(P m ) = |P m+1 ( ẋ) -P m ( ẋ)| for a certain number M of periods, where P m ( ẋ) corresponds to the intersection of the trajectory ẋ( ẋ0 , t 0 ; t) with the plane t = t 0 + mT with ẋ0 is an initial velocity at time t = t 0 and m ∈ 1, 2, . . . , M . In the case when ∆(P m ) converges to zero in the stationary regime, a periodic regime can be identified. Therefore, the membrane velocity ẋ(t) can be expressed as a Fourier series in the following form

ẋ(t) = a 0 + ∞ n=1 (a n cos (nωt) + b n sin (nωt)), (35)
and the normalized admittance is then represented by its L 2 -norm which

corresponds to (ρc/A m ) n (a 2 n + b 2 n ).
Otherwise, a quasi-periodic regime occurs and so the velocity is chosen to be represented by the absolute value of the peaks maximum amplitude. Figure 17 shows the normalized admittance Y m according to frequency f (Hz). Black curve line represents the case when the electroacoustic transducer is not connected to any electric load, where its corresponding absorption coefficient present a maximum value α max which is smaller than 1 at the resonance. Stars dashed curve illustrates the case when plugging a positive optimal resistor R opt that can be obtained by the following formula [8]

R opt = (Bl) 2 Z mc -R ms -R e with Z mc = ρcS (36)
Remarkably, the shunt resistor permits a significant decrease in the normalized admittance with a perfect absorption at the resonance. However, this approach is limited to a narrow range of frequency with no possible broadening control of the bandwidth.

Regarding the behavior of the system when connecting a nonlinear electric shunt circuit; the curve illustrated in crosses represents the L 2 -norm of the normalized acoustic admittance when reaching a periodic regime. Quasiperiodic regime illustrating the maximum of peaks absolute values is represented with dark circles. Remarkably, it exists in the vicinity of the resonance frequency ω 0 and can be identified by the vertical black dashed lines.

In the vicinity of 1:1 resonance, an optimal response frequency of the system can be identified through a selected threshold. It corresponds to the maximum of energy that the primary system can reach during an energy exchange process with the NES. Thus, the optimal design defined in terms of normalized admittance is represented by the horizontal dash-dotted light line.

Remarkably, the added passive nonlinear shunt circuit allowed a significant decrease of the admittance, principally at the vicinity of the resonance frequency where the TET prevents the velocity to exceed a certain amplitude. Moreover, we can identify that the frequency bandwidth undergoes a 45% of relative increase.

The sound absorption coefficient is usually used as a performance indicator for linear studies to optimize its bandwidth. However, underlining the fact that with the nonlinear coupling, this coefficient may not be necessarily appropriate since the interest of the nonlinear targeted energy transfer is mainly to highly-energetic non-stationary regimes. Moreover, it does not reflect the different responses of such systems (strong modulated responses, modal exchanges. . . ). Then, we choose to apply the following expression of the sound absorption coefficient α, defined for periodic regimes only as

α(ω) = 1 -|r(ω)| 2 (37)
with |r(ω)| defining the magnitude of the following reflection coefficient

r(ω) = 1 -Y m (ω) 1 + Y m (ω) . ( 38
)
As done by some authors in the nonlinear field ( [START_REF] Lee | Analytic formulation for the sound absorption of a panel absorber under the effects of microperforation, air pumping, linear vibration and nonlinear vibration[END_REF][START_REF] Lee | The jump phenomenon effect on the sound absorption of a nonlinear panel absorber 36 and sound transmission loss of a nonlinear panel backed by a cavity[END_REF]) we choose to illustrate certain properties of the nonlinear shunt circuits compared to classical linear ones, highlighting the performance of the sound absorption just for periodic regime (the sound absorption coefficient is not defined for quasiperiodic ones). Figure 18 shows a comparison of the sound absorption coefficient between the different treated cases (without coupling, with linear optimal resistor and nonlinear shunt circuits). We should note that for the nonlinear shunt circuit using the physical parameters of Test 2, the sound absorption was not defined during quasi-periodic regimes and thus not illustrated between the vertical two dashed lines. Figure 18 highlights the fact that we managed to increase the controlled absorption bandwidth and obtain better absorption for low frequencies without an optimization study for now, even if we are less efficient than an optimal linear resistor around the resonance frequency.

We should highlight he fact that the absorption bandwidth can be controlled by choosing suitable passive nonlinear control parameters R c , L c , C and k.

Remarks and interpretations

The present approach uses the complexification of Manevitch for the case of 1:1 resonance. Then, in the vicinity of ω 0 this approach presents a predictive tool allowing an analytical identification of the desirable regime for a give sound pressure level and a selected frequency. An optimal NES design can be identified when the system response does not exceeded the chosen threshold. However, for the present physical parameters, Figure 17 displays a slight difference between the admittance maximum amplitude in the quasiperiodic regime and the horizontal threshold. This gap can be due to several factors, among which we can mention the choice of the physical parameters.

In addition, we used the classical Manevitch variables which only considers the first harmonics. An extended version of the complex variables is under investigations, for a finer approximation of the result.

Conclusion

A passive nonlinear passive shunt was coupled to the loudspeaker's terminals, playing the role of an electroacoustic absorber. The resulting structure was described as a two dofs system, which was analytically treated using an invariant manifold approach for the 1:1 resonance. The equilibrium points and fold singularities were detected, allowing a predictive tool to design the nonlinear passive shunt during the energy exchange process. The preliminary results show that this approach is able to vary the acoustic properties of the loudspeaker, with a relative reduction of the normalized admittance and a broadening of the absorption frequency range.

Further work is dedicated to the design and optimization process of the passive nonlinear shunt, which allows larger frequency bandwidth and higher acoustic absorption performances. Moreover, we aim to extend the present analytical approach for a finer investigations with larger frequency range(1:2, 2:1,. . . resonances). In practice, we aim to implement an experimental setup in order to validate our analytico-numerical results.

Figure 1 :

 1 Figure 1: Shunt loudspeaker [22]

  the equivalent compliance of the en-94 closed loudspeaker (due to the closed-box at the rear side of the loudspeaker). V b represents the enclosure volume, S the diaphragm area, ρ the density of the medium and c the celerity of sound in the medium. The mechanical moving part of the loudspeaker made of the suspended diaphragm and coil is assimilated to a simple mass-spring-damper system with M ms describing the mass, C ms accounting for the surround suspension and the spider and R ms the mechanical resistance. Equation (1) describes the motion of the closed-box loudspeaker with S denoting the effective piston area and V b the volume of the cabinet. Then, the reaction of the fluid acting on the rear face is modelled by ρc 2 S 2 /V b as a mechanical compliance of the cabinet, where ρ representing the density and c the celerity of sound in the medium. v s (t)

Figure 2 :

 2 Figure 2: Circuit representation of an electrodynamic loudspeaker [8].

Figure 3 :

 3 Figure 3: Nonlinear resonant circuit shunting.

1 T T 0 P

 10 periodically varying sound pressure P (t) = A m cos (ωt) with ω the angular frequency. L p = 20 log 10 (p/p ref ) dB represents the sound pressure level with p = (t) 2 dt the root mean square sound pressure and P p ref = 20 µPa is the reference sound pressure in the air.

Figure 4 :

 4 Figure 4: Numerical and analytical results for N 1 and N 2 during free vibrations.

Figure 5 :

 5 Figure 5: τ 0 invariant where stable and unstable parts are represented in full and dashed lines respectively. N 21 and N 22 define the extrema according to the N 2 variable.

Figure 6 :

 6 Figure 6: Analytical energy trends obtained by solving Equation (31) for a corresponding constant C.

Figure 7 :

 7 Figure 7: (a) Efficient passive TET to the NES for the case of free oscillations with initial conditions corresponding to a fixed point on the upper branch of the SIM. (b) and (c) show a comparison of energy flow between the numerical integration of the initial differential System (10) in full lines and the analytical results obtained by solving Equation (31) in dashed lines.

  represent the cartography of the achievable dynamics of the system for several values of sound pressure levels (L p ). Horizontal dashed black lines define the stability boundaries N 21 and N 22 . dash-dotted and full line curves correspond respectively to f 1 (δ 2 , N 2 ) = 0 and f 2 (δ 2 , N 2 ) = 0. According to the value of L p , the positions and the number of fold singularities and equilibrium points change so as their stability. More precisely, when full line and dash-dotted curves intersect each other far from the fold lines in black (N 21 , N 22 ) we can deduce the existence of an equilibrium points which corresponds to a periodic regime. Otherwise, when they intersect on a fold line, fold singularities corresponding to quasi-periodic regime can be identified.The color shade together with the numbers from 1 to 6 correspond to the pressure levels L p = {116, 120, 124, 128, 132, 136} (dB) respectively. The stars and the triangle illustrate the positions of equilibrium points for each L p and the dots denote fold singularities corresponding to L p = 132 dB.

Figure 10 :

 10 Figure 10: Positions of equilibrium points marked with stars for the system under several values of sound pressure levels L p = {116, 120, 124, 128, 132, 136} (dB) for the detuning parameter σ = 0. g(N 2 ) = 0 corresponds to the dotted horizontal black lines N 21 and N 22 . f 1 (δ 2 , N 2 ) = 0 and f 2 (δ 2 , N 2 ) = 0 are represented in dash-dotted and full line curves respectively.

  Figure 11 shows the evolution of numerical values N 1 = |dx/dT + iΩx|, N 2 =

Figure 12 :

 12 Figure 12: (a) Comparison of velocity histories as a function of time t(s) obtained by numerical integration of the system under L p = 120 dB between the master system (The loudspeaker) coupled with the NES (the shunt nonlinear circuit) represented in solid lines and without the NES in dotted line. (b) Histories of the voltage derivative in terms of time t(s).

Figure 14 :

 14 Figure 14: (a) Comparison of velocity histories as a function of time t(s) obtained by numerical integration of the system under L p = 132 dB between the master system (The loudspeaker) coupled with the NES (the shunt nonlinear circuit) represented in solid line and without NES in dotted line. (b) Histories of the voltage derivative in terms of time t(s).

Figure 15 : 6 . 4 .Figure 16 :

 156416 Figure 15: Poincaré's section of the primary structure (Shunted loudspeaker) and the NES (nonlinear eclectic shunt circuit) showing a quasi-periodic vibrational regime.

Figure 17 :

 17 Figure 17: Comparison of the normalized admittance as function of frequency between the cases of open circuit, optimal linear resistor and the shunt nonlinear circuits. Test 1 for L p = 132 dB and Test 2 for L p = 115 dB both for k = 3468 V -2 .

Figure 18 :

 18 Figure 18: Comparison of the sound absorption coefficient as function of frequency between the cases of open circuit, optimal linear resistor, the shunt nonlinear circuits Test 1 for L p = 132 dB and Test 2 for L p = 115 dB both for k = 3468 V -2 .

Table 1 :

 1 Visaton al 170 electrodynamic moving coil loudspeaker parameters

	Parameter	Symbol Value	Unit
	DC resistance	R e	5.6	Ω
	Voice coil inductance	L e	0.9	mH
	Force factor	Bl	6.9	NA -1
	Moving mass	M ms	13	g
	Mechanical resistance	R ms	0.8	Nm -1 s
	Mechanical compliance	C ms	1.2	mmN -1
	Effective area	S	133	cm 2
	Sound celerity in the air	c	343	ms -1
	Enclosure volume	V b	10	dm 3
	Air density	ρ	1.18 Kg m -3

Table 2 :

 2 Configuration presenting the operational amplifiers parameters

	Parameter Resistance Inductance Capacitance Nonlinear coefficient
		R c (Ω)	L c (mH)	C (mF)	k (V -2 )
	Test 1	1	50	5	8160
	Test 2	-4.7	5	5	8160

Acknowledgment This work was performed within the framework of the LABEX CeLyA (ANR-10-LABX-0060) of Université de Lyon, within the program " Investissements d'Avenir (ANR-16IDEX-0005) operated by the French National Research Agency (ANR).

Appendix A.

Explicit values of variables in System (29) can be written as follows

Nonlinear Dynamics, 69(1):99-116, Jul 2012.

651