General introduction

Currently, the digital "nomads" wireless technologies have attracted significant researchers in the international scientific community, to the point that we now speak of "Internet of Everything" (IoE). The IoE is based on the idea that identifiable objects are located and controlled via the Internet. To achieve this goal, it is necessary to design embedded systems in millimeter/micrometer scales composed of wireless sensor nodes while overcoming a major drawback that is the excessive use of batteries, produced by the large number of power supply sensors. The problem is that the batteries must be changed or requires the use of chargers because their lifetimes are limited and that they are made with polluting component elements. To reduce this excessive use of pollutants and to obtain autonomy, each wireless sensor should be supplied by green energy harvesting techniques. Among the most proposed and studied solutions for micro-systems, we find essentially transducers based on mechanical vibrations using piezoelectric materials or electromagnetic energy from small coils or rectenna antennas. Use of mechanical vibrations as an exclusive excitation source is limited of ambient vibration areas and the recoverable electromagnetic energy for micro systems is often low to allow of a useful supply. One solution would be to get simultaneously both energies using materials sensitive to the electromagnetic field and the mechanical vibration such as magnetoelectric materials (ME) that combine the magnetostrictive (change of mechanical stress under an applied magnetic field and reciprocally) and piezoelectric (change of mechanical stress under an applied electric field and reciprocally) effects. Although early studies have started in the 1970s with notably the discovery of Terfenol-D ( Tb 1-x Dy x Fe 2 ), it was not until the early 2000s to see their interest arise in the international scientific community with the emergence of new magnetostrictive materials such as Metglas and piezoelectric materials including PMN-PT (Pb(Mg,Nb) O 3 -PbTi O 3 ), PZT-4/5/8 (Pb(Zr,Ti) O 3 ), and BTO (BaTiO 3 ). Homogenous analytical methods were developed to estimate ME bulk materials according to different polarizations (transverse-transverse and transverse-longitudinal) depending on the combination parameters (type, number, and thickness of layers). To validate the simulation results, bi-and tri-ME layer bulk composites were fabricated by sticking the magnetostrictive and piezoelectric materials with adhesive layers of Epoxy-type. Experimental results of ME coefficients have confirmed the possibility to obtain a few of V/(cm•Oe) in no-resonant regime and few tens of V/(cm•Oe) in resonant regime. In case of classical laminate bulk material (Terfenol-D/PZT/Terfenol-D), the delivered powers into optimal impedance are in the order of mW/ cm 3 . Nevertheless, it has been noticed that due to fatigue and inhomogeneous thicknesses of adhesive layers, the mechanical coupling can be degraded over time. This affects the magnetostrictive response and consequently the effective delivered power.

Since 2010, many studies have been devoted to the design and the miniaturization of new ME composite materials coupling the giant magnetostrictive and piezoelectric effects. Thanks to thin film deposit processes, new robust combinations such as the new alloys CFO (CoFe 2 O 4 ), FeONi, FeGa and FeCo can be created laying aside rare earth materials. At this scale level, the previous homogenous analytical methods do not allow of accurate modelling of the coupling phenomena and do not take into account the adhesive layers and the mechanical effects of electrodes.

It is in this context that our research work at the Electronic and Electromagnetism Laboratory (L2E) propose a contribution to modelling of magnetoelectric composites for energy harvesting. This thesis is composed of four chapters. The first chapter reproduces in a general context properties on magnetoelectric composites. After a brief history, the magnetoelectric effect applications and the modelling and characterizations of magnetoelectric materials and devices are presented.

The second chapter presents the different numerical simulation tools to be used and developed in this thesis. Firstly, an analytical numerical method based on 0D-assumption modelling in static regime is presented. Secondly an analytical numerical method based on 1D-modelling using an equivalent electrical circuit is studied. Thirdly, a 2D multiphysics code based on the finite element method is presented.

The third chapter investigates the deliverable output of a magnetoelectric laminate composites composed of Terfenol-D/PZT-5A/Terfenol-D materials in considering electrical load. The model has been developed for both harmonic and transient cases in considering a SSD technique.

The fourth chapter presents a potential application in the biomedical domain and show measurement realization on a bilayer magnetoelectric laminate composites composed of Terfenol-D/PZT-5H.

Chapter 1. General context -State of Art 1.1 Magnetoelectric effect and materials

Magnetoelectric (ME) effects exist in multifunctional active materials, they refer to the electric polarization induced by applied magnetic field, and the magnetization induced by applied electric field. The effects can be linear or non-linear with respect to the external fields. In ambient temperature, they can be measured and defined as 𝛼 𝑉 = 𝜕𝑉 𝜕𝐻 ⁄ (called the ME voltage coefficient), or 𝛼 𝐸 = 𝜕𝐸 𝜕𝐻 ⁄ (called the ME electric field coefficient), with V, E and H denoting the induced voltage, induced electric field and applied magnetic field, respectively. Researches on ME materials date back to the 1890s. Such materials now widely exist in various applications, in the form of either single-or multi-phase, with a tendency of gradually increasing performances.

Single-phase ME materials

In single-phase ME materials, the crystal symmetry plays a key role that determines ME effects in the material, which only shows up below a certain temperature, called the Curie temperature (also called the Neel temperature in some cases). A non-exhaustive history is briefly listed as follows.

In 1894, Curie stated "Les conditions de symétrie nous montrent qu'il pourrait se faire qu'un corps à molecules dissymétriques polarise dialectiquement lorsqu'on le place dans un champ magnétique ... Et peut-être polarisé magnétiquement lorsqu'on le place dans un champ électrique" [START_REF] Curie | Sur la symétrie dans les phénomènes physiques[END_REF], which firstly revealed the existence of intrinsic magnetoelectric effects in single-phase materials, due to asymmetric molecule structures.

In 1958, Landau and Lifshitz proved the possibility to obtain ME effects for certain crystals, according to their crystal symmetry [START_REF] Landue | Electrodynamics of Continuous Media[END_REF].

In 1960, Dzyaloshinskii [START_REF] Dzyaloshinskii | On the magneto-electrical effects in antiferromagnetics[END_REF] discovered ME effects in antiferromagnetic Cr 2 O 3 . Later, Astrov [START_REF] Astrov | The magnetoelectric effects[END_REF] experimentally demonstrated the electrically induced ME effect in the same material, and reported a ME coefficient of 20mV/(cm• Oe) (under the temperature range of 80 -330K, and an applied alternating electric field at 10kHz frequency in the orientation of the symmetrical axis) .

In 1963, Foner and Hanabusa [START_REF] Foner | Magnetoelectric effects in Cr 2 O 3 and (Cr 2 O 3 ) 0.8 (Al 2 O 3 ) 0.2[END_REF] investigated the effect of temperature on the critical magnetic field (for ME effects) of both Cr 2 O 3 and (Cr 2 O 3 ) 0.8 (Al 2 O 3 ) 0.2 . They also found that a greater ME coefficient could be obtained when the material had been annealed under magnetic fields.

In 1963, Shtrikman and Treves [START_REF] Shtrikman | Observation of the magnetoelectric effect in Cr 2 O 3 powders[END_REF] also discovered ME effects in polycrystalline powder compacts that were annealed in electric and magnetic environments.

In 1969, Hornreich [START_REF] Hornreich | The magnetoelectric effect: some likely candidates[END_REF] predicted a family of potential ME materials, including DyFeO 3 , Fe 2 TeO 6 , Cr 2 TeO 6 ,FeCrWO 6 , Cr 2 WO 6 , V 2 WO 6 , FeSb 2 O 4 , Ca 2 FeAlO 5 , Eu 3 O 4 and β -FeNaO 2 -FeNaO 2 .

In 2003, Ramesh's group successfully manufactured an epitaxial multiferroic material -BiFeO 3 (bismuth ferrite), which possesses both "ferroic" characteristics (ferromagnetism and ferroelectricity), and thus shows significant ME coefficients at room temperature, being orders of magnitudes larger than other single-phase ME materials [START_REF] Wang | Epitaxial BiFeO 3 multiferroic thin film heterostructures[END_REF].

Between 2003 and 2006, ME effects were also found in antiferromagnetic, weak ferromagnetic and ferromagnetic materials like TbFeO 3 and Fe x Ga 2-x O 3 . See [START_REF] Schmid | Magnetoelectric effects in insulating magnetic materials[END_REF], [START_REF] Spaldin | The renaissance of magnetoelectric multiferroics[END_REF], [START_REF] Fiebig | Revival of the magnetoelectric effect[END_REF] and [START_REF] Eerenstein | Multiferroic and magnetoelectric materials[END_REF].

It should be noted that almost all aforementioned single-phase ME materials are limited in realistic applications, as a result of their very low Curie temperatures, and weak ME coefficients. Consequently, exploiting ME applications has been largely based on the development of new composite materials that elaborate and enhance the properties of single-phase ones.

Two-phase ME composite materials

Recently, two-phase ME materials, also called composite materials, have sparked increasingly interest in various applications, thanks to their promising properties. On one hand, composite materials exhibit colligative properties (in terms of e.g. density, compliance, etc.) of constituent phases. In such cases, properties can be determined quantitatively by averaging properties of the constituent components using volume fractions. On the other hand, composite materials exhibit enhanced product properties of their constituent single phases. Take the two-phase ME composite consisting of ferromagnetic (M) and ferroelectric (E) materials, for example. The M phase converts energy between magnetic and elastic forms (i.e. magnetostrictive effects), whereas the E phase converts energy between electric and elastic forms (i.e. piezoelectric effects). When combined together, energy can be converted between magnetic and electric forms, via the media of elastic bonding (as depicted in Fig. 1.1), thus, exhibiting new properties that do not belong to any of its individual constituent phases. In 1972, van Suchtelen initialized the concept of product property of ferromagnetic and ferroelectric material composites [START_REF] Van Suchtelen | Product properties: a new application of composite materials[END_REF].

In 1978, van den Boomgaard and Born summarized the key aspects to prepare performing ME composites [START_REF] Van Den Boomgaard | A sintered magnetoelectric composite material BaTiO 3 -Ni(Co, Mn)Fe 3 O 4[END_REF]: (i) constituent phases should be in equilibrium, (ii) there should be no mismatching between grains, (iii) for those consisting of magnetostrictive and piezoelectric materials, the involved magnetostrictive and piezoelectric coefficients should be large, (iv) there should be no leak of accumulated charge through magnetostrictive or piezoelectric phases, and (v) appropriate poling procedures. [START_REF] Li | Study of magnetoelectric composites for sensor and transducer applications[END_REF] At the same time, magnetoelectric composite materials consisting of magnetostrictive and piezoelectric phases can be categorized with respect to the geometry structure formed by their constituent phases. According to Newnham et al, most fabrication processes of the ceramic composites are based on an initial mixing of starting powders, followed by e.g. pressing, sintering, densification and solidification, to a net-shape. As shown in Fig. 1.3, different categories refer to the volume fraction of one constituent phase in the entire composite. If the concentration of one phase is so small that the phase turns out to consist of isolated particles in a matrix, and this composite will be referred to as a 0-3 type (one phase interconnected in zero dimensions and the other in three). If the volume fraction of one phase increases to reach an initial percolation limit, the composite will be classified as a 1-3 type (for example, the cylinder matrix form). If the volume fraction of one phase crosses another percolation limit and starts to be interconnected in two dimensions throughout the entire composite, it results in a 2-2 type. For the category of 0-3 bulk composites, several limitations exist, which leads to issues of reproducibility and reliability, as listed as follows.

1) chemical reactions which is difficult to prevent from taking place due to the sintering process, in which high temperatures are involved.

2) dielectric breakdown of the magnetostrictive phase, during the poling process for the piezoelectric phase.

3) process caused mechanical defects like cracks, pores, etc. resulting in poor elastic bonding between constituent phases.

Besides the deficiencies mentioned above, another significant limitation for development and utilisation of the bulk 0-3 composite is the lack of reliable theoretical models, with which precise prediction of the effective properties of this type of inhomogeneous structure can be achieved. Until now, only some approximate approaches have been developed to model this two-phase composite, such as the matrix systems and two-component mixtures [START_REF] Bichurin | Modeling of Magnetoelectric Effects in Composites[END_REF]. However, they also possess their specific applicability limits. Therefore, more details about further studies of bulk 0-3 composite will not be presented in the rest of the dissertation.

Laminated 2-2 composites

Composites belonging to the 2-2 category are also referred to as laminated composites. They have better properties than their bulk counterparts due to their configuration, such as easier manufacturing control and accurate volume fraction tuning. Nevertheless, the co-sintering or adhesive bonding between each laminate composite inevitably leads to defaults in the interfaces, which involves mechanical losses.

For oxide magnetic material based laminate ME composites, there are two main deficiencies: (i) interphase diffusions resulting from the co-firing process, which deteriorates the coupling effect, and (ii) micro cracks or pores generated due to thermal expansion in the cooling process.

Rod matrix 1-3 type composites

The 1-3 category is also widely employed in various applications. This configuration contains a rod (in form of square rod, cylinder, etc.) array of piezoelectric materials embedded in continuous magnetic phases or magnetic/epoxy mixture matrix structure [START_REF] Lam | Frequency response of magnetoelectric 1-3-type composites[END_REF]. This kind of composite can be fabricated using a dice-and-fill technique consisting of the following steps: (i) poled piezoelectric ceramic sample (generally in disk form) cut by two mutually perpendicular sets of precise saw to form the array of piezoelectric rods; (ii) solution with magnetic phase particles and epoxy slurry is filled into the diced piezoelectric sample; (iii) a hardening process is imposed on the slurry-filled sample, and after polishing and silver painting electrodes on the top and bottom surfaces, the 1-3 type composite is fabricated. The volume ratio of the piezoelectric phase can be controlled by modifying the periodicity of the rod array. The 1-3 type composite may also exhibit giant ME effects while eliminating the brittleness defect of 2-2 composite structures. See e.g. [START_REF] Lam | Frequency response of magnetoelectric 1-3-type composites[END_REF], [START_REF] Shi | Magnetoelectric effect of Pb(Zr, Ti)O 3 rod arrays in a (Tb,Dy)Fe 2 /epoxy medium[END_REF].

Magneto-mechanical effect and material 1.2.1 Magneto-mechanical coupling

Magneto-mechanical couplings in ME composites involve stress, strain, magnetic field and magnetization. Such couplings exist in a large majority of ferromagnetic materials and rare earth element materials (e.g. Terfenol-D). They can be termed as either the Joule effect or the Villari effect. The former implies mechanical elongation or contraction of the material, in response to external magnetic fields. This effect was discovered by James Prescott Joule in iron bars [START_REF] Joule | On the effects of magnetism upon the dimensions of iron and steel bars[END_REF]. The Villari effect, on the other hand, implies changes in magnetization due to mechanical stress, which was found by Villari [START_REF] Villari | Change of magnetization by tension and electric current[END_REF].

Magneto-mechanical couplings stem from interactions of electric spins (i.e. orbital angular momenta) with the crystal lattice. In fact, lattice dimensions and spacing are strongly influenced by the magnetization, which is caused by alignment of large amount of magnetic moments (due to a coupling between the orbital angular momenta and electric spins exhibited in most of the lanthanide series rareearth elements), forming magnetic domains along the easy axis of magnetization (called the Weiss mean field). Two steps are involved during the process. Firstly, displacements of walls (that separate magnetic domains), which are dominant under weak field, and rotations of domains, which are dominant under strong fields. Such phenomena for single crystals are illustrated in Fig. 1.4. The material attains saturation when all magnetic moments are aligned (and thus, forming a single domain).

In this case, the generated strain also attains its maximum, called the saturation magnetostriction that is denoted as λ s . It should be noted that several factors play important role in such phenomena, including temperature, magnetic fields, and stresses.

Figure 1.4 Magnetization process of a single crystal cubic [START_REF] Nguyen | Modélisation par éléments finis de matériaux composites magnéto-électriques[END_REF] As illustrated in Fig. 1.5, spontaneous magnetostrictions (denoted by λ 0 ) occur, even under null external magnetic field, when the material is under Curie temperature Tc. Below Tc, with the increasing of magnetic fields, magnetic moments experience the isotropic, disordered paramagnetic, and ordered ferromagnetic states in order. At the final stage, saturation magnetostrictions are achieved. Above Tc, no magnetostrictive effect exhibits. This is because in this case the thermal energy overcomes the Weiss mean field, thereby alignment of moments is destroyed, and thus formation of domains is impossible. As shown in Fig. 1.6, when applying (either compressive or extensive) stresses to magnetostrictive materials along a certain direction, the distribution of magnetic moments can be changed from random or isotropic states into anisotropic ones. It is worth mentioning that the relations between excitation magnetic fields and magnetostrictions are intrinsically nonlinear. Indeed, the material behaviour has nonlinear dependencies with respect to both mechanical and magnetic quantities, as shown in Fig. 1.7 where measurements for Terfenol-D are depicted [START_REF] Nguyen | Modélisation par éléments finis de matériaux composites magnéto-électriques[END_REF]. Nonetheless, when the material operates under static bias (mechanical and magnetic) conditions super-positioned with dynamic magnetic fields of small magnitudes, the behaviour can be considered as linear. In this case, it is referred to as piezomagnetic.

Magnetostrictive material

Magnetostrictive materials refer to a branch of ferromagnetic materials that exhibit magnetostrictive effects. In practice, they are utilized for energy conversion between magnetic and elastic forms, which is depicted in the diagram of Fig. 1.8. Since the discovery of magnetostrictive effects in "iron bar" by Joule, magnetostrictive materials that had been researched were limited within transition metals e.g. Nickel (Ni), Cobalt (Co), and Permalloy (Ni-Fe-Mo). Such materials had also been widely used in telephones and sonar systems during the 20 th century [START_REF] Hunt | Electroacoustics: the analysis of transduction and its historical background[END_REF]. Further applications were limited due to the low saturation magnetostriction of these materials (normally 10 to 100 ppm [START_REF] Restorff | Magnetostrictive materials and devices[END_REF]). During the 1960s, it was discovered, with the help of the technology of neutron diffraction measurements, that certain heavy rare-earth metal elements possess promising magnetic properties, including: (i) more sophisticated spin structures (compared with classical ferromagnetic and antiferromagnetic materials), (ii) significant magnetic anisotropies, and (iii) huge magnetic moments (resulting from parallel couplings of large orbital and spin angular momenta) that were an order of magnitude larger than those of Ni and Fe [START_REF] Li | Study of magnetoelectric composites for sensor and transducer applications[END_REF]. Then, in 1971 witnessed a historic breakthrough in the development of magnetostrictive materials: the US Naval Ordnance Laboratory fabricated an alloy consisting of metal elements Tb, Dy and Fe, which was later named Terfenol-D. The latter exhibits a giant magnetostriction as large as (up to) 2000 ppm in ambient temperature; that is nearly 100 times larger than those of previously developed transition metals and alloy materials [START_REF] Cedell | New magnetostrictive alloy for rapid conversion electric energy to mechanical motion[END_REF].

The stoichiometry of Terfenol-D can be represented as Tb x Dy 1-x Fe y , in which the sample for x = 0.3 and y = 1.92 has been commercialized by some US companies (among which is the ETREMA Products). The giant magnetostriction of Terfenol-D stems from the anisotropic property of its cubic crystallographic structures. To be more specific, the structure features in certain special crystal orientations, which favours the alignment of magnetic domains when magnetic fields are applied along (or in the adjacency of) these crystal orientations. As a result, the saturation magnetization becomes much more significant. More recently, a Terfenol-D based composite has been created. This composite consists of Terfenol-D particles and insulating polymer matrix. It has high resistivity and operates over a broader frequency regime, without suffering from mechanical brittleness and large eddy current losses (at high frequencies) of Terfenol-D [START_REF] Sandlund | Dynamic characterization of 1-3 type Terfenol-D/epoxy magnetostrictive composite[END_REF][START_REF] Ruiz De Angulo | Magnetostrictive properties of polymer bonded Terfenol-D[END_REF]. In the following table, we depict characteristic constants of typical magnetostrictive materials, including Terfenol-D.

Table1.1 Saturation magnetostriction and magnetic induction of some magnetostrictive materials [START_REF] Nguyen | Modélisation par éléments finis de matériaux composites magnéto-électriques[END_REF] 1.3 Electro-mechanical effect and materials

Electro-mechanical coupling

Electro-mechanical couplings refer to piezoelectric effects, in which mechanical quantities like stress and strain, as well as electrical quantities like polarization and electric field, are involved. Piezoelectric effects exist in many dielectric materials; they can be divided into the direct piezoelectric effect, and the converse piezoelectric effect.

Figure 1.9 Microscopic origin of the piezoelectricity [START_REF] Nguyen | Modélisation par éléments finis de matériaux composites magnéto-électriques[END_REF] The direct effect was discovered by Pierre Curie in the 1880s. It implies that when piezoelectric materials are applied with mechanical stresses, electric polarization will be developed inside the material. The relationship between resulting electrical outputs and mechanical inputs are usually referred to as the laws of piezoelectricity. The underlying mechanism of the direct piezoelectric effect is that when piezoelectric materials are under the Curie temperature, there is no symmetry centre (see Fig. 1.9) within its crystal structures. Hence, applying stresses leads to unbalance between positive and negative charges, which, in turn, generates electric polarization, or electric charges on surface of the material. On the other hand, when the material is beyond the Curie temperature, the inner structures are centrosymmetric. In other words, materials are isotropic, for which applying stresses has no effect on polarization. Meanwhile, it is noted that the orientation of applied stresses (with respect to the poling direction) has significant impact on the generated charges. This phenomenon is referred to as working modes, which are illustrated in Fig. 1.10 (where we note the effects of a tensile or compressive stress, as well as orientations of stresses). The converse piezoelectric effect is also termed as the Villard effect. It implies that applying external electric fields on piezoelectric materials causes changes in volume of the latter. As illustrated in Fig. 1.11, for piezoelectric materials with a given poling direction, electric fields along the same direction lead to extension of the material along the poling direction, whereas electric fields along the opposite direction lead to contraction along the poling direction. Meanwhile, it is also worth mentioning that the pre-polarization (as depicted in Fig. 1.12) plays an important role in the functioning of piezoelectric materials. In fact, for piezoelectric materials at their original states (i.e. no electric fields applied), electric domains that are separated by walls along their individual spontaneous polarization directions, are randomly distributed, therefore, exhibiting no macroscopic polarization. When weak electric fields are applied, piezoelectricity takes place, which is mainly due to crystalline asymmetry of inner structures (i.e. the so-called intrinsic piezoelectricity).

When the applied field gets stronger (than the coercive field Ec of the material), dipole domains align themselves proximately parallel to the field, thereby domains of orientations close to the field experience growing in volume while volumes of others are reduced. It is at the same time that macroscopic mechanical deformation occurs. The deformation (and thus, the polarization, too) reaches its maximum when the applied field overpasses certain level of magnitudes, which is referred to as saturation. On the other hand, removing the applied field does not lead to vanishing of polarization, and the remaining part is called the retentive polarization (denoted by Pr). In fact, the latter is of critical importance for enhancing performance of piezoelectric materials, as it yields significant coupling effects in the material (see blue curves in the following figure that are generated due to Pr).

Piezoelectric material

In analogue to magnetostrictive materials, piezoelectric materials refer to a branch of ferroelectric materials that exhibit piezoelectric effects. As shown in Fig. 1.13, conversion between electric and mechanical energies is involved inside piezoelectric materials. Nowadays, piezoelectric materials are largely seen in a wide range of applications. For instance, the early found single piezoelectric crystals, including quartz, Rochelle salt, and tourmaline are extensively employed in resonators [START_REF] Parzen | Design of Crystal and Other Harmonic Oscillators[END_REF]. Polycrystalline piezoelectric materials such as the lead zirconate titanate (PZT), barium titanate (BiFeO 3 ), and lead metaniobate (PMN), on the other hand, possessing relatively higher coupling coefficients, are broadly utilized in acoustic transducers [START_REF] Medill | A miniature piezoelectric microphone[END_REF], motors and piezoelectric transformers [START_REF] Crawford | Piezoelectric ceramic transformers and filters[END_REF]. In Table 1.2 we depict piezoelectric coefficients of typical piezoelectric materials. 1.4 Magnetoelectric effect applications

Magnetic field sensors

ME composite materials that inherit virtues from both piezoelectric and magnetostrictive materials are widely used in various applications. A typical one is for measurements of magnetic fields. Thanks to the significant piezoelectric and magnetostrictive coefficients of its constituent phases, magnetic field sensors made of ME composites normally have outstanding sensitivity. Additionally, outputs can be conveniently read out in terms of electric quantities. In the following, we introduce magnetic field sensors in two categories: those used for static fields measuring and others for dynamic fields measuring.

Static magnetic field sensor

The first category of sensors can be used to measure static magnetic fields, such as the earth magnetic field. They are based on the nonlinearity of the magnetostrictive phase of ME composites. Under static magnetic fields, strains are developed in the magnetostrictive phase, which are then transferred to the piezoelectric phase and lead to generation of electric charges on surface of the latter. The magnetic to electric signal relationship can be quantified based on material properties of constituent phases. Hence, values of magnetic fields can be determined by the output electric signals. In practice, several types of configurations and mechanisms can be employed, in order to ensure satisfying sensibility and range of measurements. An example of static magnetic field sensor is depicted in Fig. 1.14. It consists of three bilayer ME structures of different optimal bias static magnetic fields (with respect to the maximum output voltage). We see the combining effects of the three layers on the most right of the figure.

Dynamic magnetic field sensor

The second category of sensors is for measurements of dynamic magnetic fields [START_REF] Dong | Ultrahigh magnetic field sensitivity in laminates of Terfenol-D and Pb(Mg 1/3 Nb 2/3 )O 3 crystals[END_REF]. An example of a trilayer ME structure is displayed in the figure below. This sensor is able to detect weak signals (between the range of 10 -3 to 10 -11 Tesla). Furthermore, it is able to measure magnetic fields over a wide range of frequencies (to be more precise, from one kHz to several hundreds of kHz).

Energy harvesting applications

With the rapid development of the electronic and telecommunication industries, the problem of regular energy dependence is increasing: (i) the energy source is still relative in shortage and the waste amount is great. (ii) most of the embedded batteries in use have a short lifetime; (iii) to fulfil the functions and the charging of the energy source depend highly on electricity and are quite wire-rely in most fields. In order to face the energy challenge, it is now necessary to consider new technologies such as the green energy harvesting to provide efficient solutions. The main principle is to get maximum diversity of usage in ambient energy sources including light, body or machine heat, electromagnetic field, motion and vibration of mechanics and so on. Various energy harvesting devices and systems have been proposed and realized [START_REF]Energy Harvesting for No-Power Embedded Systems[END_REF]. Among the most proposed and studied solutions for micro-systems, we find essentially transducers based on mechanical vibrations using piezoelectric materials or electromagnetic energy from small coils or antennas. Use of mechanical vibrations as an exclusive excitation source is limited in ambient vibration areas, and the recoverable electromagnetic energy for micro systems is often low to allow of a useful supply. One solution would be to get simultaneously both energies using materials sensitive to the electromagnetic field and the mechanical vibration such as magnetoelectric materials (ME) that combine the magnetostrictive and piezoelectric effects.

Modelling and Characterizations of magnetoelectric materials and devices 1.5.1 Theoretical modelling methods

During these last 40 years, the researches and publications about the magnetoelectric effect and materials have been made with large and rapid growing amount. However, the investigations concerning the numerical modelling have played a relatively minor role throughout all the researches. It is evident that a precise and efficient modelling is able to provide rapid and reliable guides for the design and optimization of magneto-electric applications and could thus largely reduce the trial costs for the development of the devices. Attribute also to the fact that the rapid development of the computer technologies (large memory, high speed CPU, etc.) has made the implementation of the numerical calculates much more faster and easier than before, it is therefore significant and even urgent to establish and improve more rigorous numerical modelling approaches, which will be quite useful for the practical fabrications.

The modelling of the ME effect calls for macroscopic approaches based on the constitutive laws and equilibrium equations of the materials studied. However, the complex multi-physics phenomena that occur inside this type of materials (including, for instance, the nonlinear magneto-mechanical or electro-mechanical coupling characteristics, the electromechanical or magnetic-acoustic frequency effects and the dynamic electromagnetic excitation based on the Maxwell equations) result in a big difficulty in creation of the macroscopic modelling methods currently. The major object of this dissertation is the establishment and improvement of such efficient approaches, as well as the practice of them to guide the design and development of the devices with the ME materials as their core components.

As most of the physical problems can be described mathematically by the partial differential equations (PDE) and resolved by numerical methods through discretization of the equations system, the method of discretization is significantly necessary. The finite element method (FEM) is one of the macroscopic approaches of discretization which operates as replacing the PDE defined on a continuous problem area with a system of discrete equations to be solved by numerical methods. It is especially adapted for the physical problems with complex geometries and multi-physical phenomena couplings. There exist different sorts of FEM such as the variation method and the weighted residual method. The former operates generally as resolving the PDE by minimizing the energy functions of the problem system, but it is limited in usage for complex physical couplings. The latter operates based directly on the equilibrium equations and boundary conditions of the problems, thus it is more convenient and will be therefore the first choice for the modelling methods utilised in this dissertation.

There also exist already plenty of commercial software, such as ANSYS and COMSOL, which are capable of complex finite element calculation, but their programme code sources do not usually corporate adapted model, especially for resolution of nonlinear cases. Thus, another important work in the dissertation is to develop new programs for the studied complicated problems by using certain coding platforms such as MATLAB or C++.

There have been many researches upon FEM modelling for magneto-mechanical and electromechanical coupling problems. For the first type of problems, a model for static nonlinear magnetomechanical problems has been developed by Azoum by means of decomposing the total deformation of the material into an original mechanical deformation and a magnetostriction with a quadratic relation representing the dependence of the magnetostriction upon the magnetic induction [START_REF] Azoum | Contribution à la modélisation numérique de phénomènes magnéto-élastiques. Étude de dispositifs à base de matériaux magnétostrictifs[END_REF]. And a plate-form of the magnetostrictive materials has been developed by Galopin in order to integrate the relevant magneto-mechanical constitutive law into his static ME model [START_REF] Galopin | Modélisation et caractérisation de matériaux actifs pour la conception de dispositifs magnéto-électriques[END_REF]. The simulation results are compared with those of experiments from Ueno and Higuchi [START_REF] Ueno | Magnetic sensor for high temperature using a laminate composite of magnetostrictive material and piezoelectric material[END_REF]. And for the second type of problems, a model for both static and dynamic problems of multilayer piezoelectric material has been developed by Piefort by setting the mechanical deformation and electric field as the state variables and integrating the linear electro-mechanical constitutive laws [START_REF] Piefort | Finite Element Modelling of Piezoelectric Active Structures[END_REF]. The constitutive laws based on the Helmholtz free energy function was proposed by Belahcen and the relevant model was applied on calculation for the electric machine vibration problems [START_REF] Belahcen | Magnetoelasticity, magnetic forces and magnetostriction in electrical machines[END_REF].

In addition, the magnetoelectric effect in harmonic behavior has been studied by Liu and Wan [START_REF] Liu | Numerical modeling of magnetoelectric effect in a composite structure[END_REF] by simulating a piezoelectricity module in ANSYS under the assumption that the magneto-mechanical and electro-mechanical couplings are linear in the original model. Then, the deformation, the magnetic field and the electric field are all set as the state variables and the Foucault current is ignored. Next, both a static and a dynamic magnetic field of much smaller magnitude are applied on the material studied to obtain a larger output voltage, and a curve of the mechanical strain in function of the applied static magnetic field is obtained by experiment for the purpose of integrating the nonlinear magnetoelectric characteristics into the model. Thus at each of the static magnetic field values, a coupling coefficient can be obtained around this magnetization point and utilized for the dynamic analysis. The results demonstrated that the resonance frequency obtained by the model is very close to that of the mechanical resonance of the material, and the output voltage is dependent of the applied static magnetic field for the reason of the considered nonlinear magnetic-mechanical coupling explained by the model [START_REF] Liu | Effect of magnetic bias field on magnetoelectric coupling in magnetoelectric composites[END_REF].

The equivalent circuit method (ECM) is another important analytical modeling approach for the resolution of harmonic problems with ME effect. The working principle of the method is to replace the studied ME system with a pseudo-circuit module that can describe the magneto-electro-mechanical couplings on the basis of the equilibrium equations and behavior laws. Thus the ME coefficient of the system can be obtained simply through calculation of the formulations from ECM, which is of much less complex and of less time expense than the implementation of FEM. Many research studies have been made to develop the ECM formulations and the results are compared with from the experiments [START_REF] Dong | Longitudinal and Transverse Magnetoelectric Voltage Coefficients of Magnetostrictive/Piezoelectric Laminate Composite: Theory[END_REF][START_REF] Yang | Resonant magnetoelectric response of magnetostrictive/piezoelectric laminate composite in consideration of losses[END_REF]. However, it should be importantly noted that most of the ECM presented in the literatures are based on certain simplified 1-dimension assumptions such as the neglect of the mechanical variable components in transverse direction of the composite, and the electric and magnetic variable components are uniformly distributed along the transverse direction. All the assumptions inevitably lead to strong restrictions for the ECM and obvious discrepancies in values from the results of ECM, FEM and experiments for problems with some complex composite configurations.

Experimental Characterization

For the purpose of acquiring actual property characterizations of the materials to be studied, and subsequently of verifying the availability of the current and modified theoretical models, measurement experiments need to be designed and fulfilled. The results will then be compared to those from the theoretical models for further improvement to better fit the measurement.

Under the limited condition of the laboratory and research circumstance, some fundamental characterization measurements of the studied magnetostrictive materials, piezoelectric materials and their composite in some simple forms will be carried out, such as the DC bias and frequency dependent ME coefficients of the bilayer MP composite structure. In the measurement setup, the composite samples will be mechanically supported and fixed by a plastic framework; the external applied magnetic field will be provided by a pair of permanent magnets and copper search coil around the sample via a dynamic signal generator; the output electrical parameter results such as the voltage will be displayed by a dynamic signal analyzer and recorded for further calculation and analysis. Finally with the experimental results the property characterization curves will thus be calculated and presented to make comparison with the theoretical results. The agreements and discrepancies will be discussed and interpreted.

Conclusion

In this chapter, the working mechanisms of magnetostrictive and piezoelectric materials are discussed in detail and the definition as well as the development history of the magnetoelectric effect and materials have been introduced. The effect and composite materials also possess a great potential and wide application in the study of micro-scale magnetoelectric devices. However, it is also important to note the shortcomings of the current theoretical models and the necessity to achieve more suitable and more rigorous modeling. The principle objective of this dissertation is to develop a new model with FEM to solve the static nonlinear magnetostrictive problems. And subsequently, the problems in harmonic area with the frequency effect can be simply treated as linear around one point of magnetization and then resolved by utilizing the ECM. The improvement of ECM will be furtherly studied by taking into account the 2D/3D effects to solve some expansion of ME problems based on various configurations such as the multilayer composites. Most of them can be predicted to be quite applicable for the development and optimization of magnetoelectric materials and devices.

Chapter 2. Analytical and Numerical Modelling of Magnetoelectric Composites

Introduction

Investigation of magnetoelectric problems requires the development of numerical and analytical methods that combine the mechanical and electromagnetism governing physics equations with the electro-magneto-mechanical laws. This chapter proposes to expose the methods investigated in this thesis that address the magnetoelectric composites in various configurations. For that, the mechanical and electromagnetism governing physics equations and the electro-magneto-mechanical laws are presented firstly. Secondly, we present methods that must be distinguished, according to distinct groups. The first group uses analytical methods (simplified analytical resolution using the 'matrix' and equivalent circuit method) with the homogenous quasi-static approximation and the second group uses a numerical method based on the finite element method (FEM).

The mechanical-electric-magnetic properties and constants parameters of magnetostrictive and piezoelectric materials used here are given in the Appendix A.

The development of the modelling methods aims to obtain accurate values of the output coefficients (ME coefficient, output power, etc.) of laminated magnetostrictive (simply noted as 'M', the same as following)-piezoelectric (simply noted as 'P', the same as following) composites under both static and dynamic magnetic excitations in various magnetization-polarization modes. With these theoretical results the effect of the important parameters of the ME composite on the output coefficients would be demonstrated, and thus the optimization for design of the devices utilising the composite would be studied in a more sensible way.

Electromagnetic and mechanical governing equations and constitutive laws

The electromagnetism equilibrium equations are given by the Maxwell equations:

div𝑩 = 0 (2-2-1) rot𝑯 = 𝑱 + 𝜕𝑫 𝜕𝑡 (2-2-2) div𝑫 = 𝜌 (2-2-3) rot𝑬 = - 𝜕𝑩 𝜕𝑡 (2-2-4)
where B is the magnetic induction, H the magnetic field, J the current density, D the electric displacement, E the electric field and ρ the charge density.

These equations are completed with the following constitutive laws:

𝑩 = 𝜇𝑯 (2-2-5) 𝑫 = 𝜀𝑬 (2-2-6) 𝑱 = 𝜎𝑬 (2-2-7)
where 𝜇 is the magnetic permeability, 𝜀 the dielectric permittivity, 𝜎 the electric conductivity. If the material is linear, all these coefficients are constant. The non-linear case of the permeability will be introduced in section 2.3. In addition, the displacement current in the composite is negligible and ignored in our study.

From the equations (2-2-1) and (2-2-4), a magnetic vector potential 𝒂 and an electric scalar potential 𝑉 can be introduced by

𝑩 = rot(𝒂) (2-2-8) 𝑬 = -grad(𝑉) - 𝜕𝐚 𝜕𝑡 (2-2-9)
In our case, the eddy current and the displacement current will not be take into account, thus 𝑬 = -grad(𝑉) and rot𝑯 = 𝑱.

For the study of mechanical properties, the elastodynamic equilibrium equation for a continuous medium is given by:

𝑑𝑖𝑣𝑻 + 𝒇 = 𝜌 𝑚 𝜕 2 𝒖 𝜕𝑡 2 (2-2-10)
where 𝑻 is the mechanical stress, a 2-order tensor, 𝒇 the external volume force, 𝜌 𝑚 the density of material and 𝒖 the displacement vector defined in a 3D orthogonal coordinates a

𝒖 = [ 𝑢 𝑥 𝑢 𝑦 𝑢 𝑧] 𝑡 (2-2-11)
This equation is completed with the following compatibility relation between the strain 𝑺 and the displacement 𝒖:

𝑺 = 𝒟𝒖 (2-2-12)
where 𝒟 is a gradient transformation operator described in 3D orthogonal coordinates as

𝒟 = [ 𝜕 𝜕𝑥 0 0 0 𝜕 𝜕𝑦 0 0 0 𝜕 𝜕𝑧 0 𝜕 𝜕𝑧 𝜕 𝜕𝑦 𝜕 𝜕𝑧 0 𝜕 𝜕𝑥 𝜕 𝜕𝑦 𝜕 𝜕𝑥 0 ] (2-2-13)
The mechanical stress 𝑻 and strain 𝑺 can be written with the Voigt notation as: The magnitudes of the magnetic, electrical and mechanical variables in connected mediums are generally discontinuous due to their different properties. By taking an infinitesimal on the boundary and utilising the equilibrium equations and the Gauss's law, the interface continuous conditions of the physical variables can be obtained as:

𝑻 = [
(𝑩 𝟐 -𝑩 𝟏 ) • 𝐧 = 𝟎 (𝑯 𝟐 -𝑯 𝟏 ) × 𝐧 = 𝐉 𝐬 (𝑫 𝟐 -𝑫 𝟏 ) • 𝐧 = 𝜌 𝐬 (𝑬 𝟐 -𝑬 𝟏 ) × 𝐧 = 𝟎 (2-2-18) (𝑻 𝟐 -𝑻 𝟏 ) • 𝐧 = 𝐟 𝐬 (𝑺 𝟐 -𝑺 𝟏 ) × 𝐧 = 𝟎
where 𝐉 𝐬 , 𝜌 𝐬 and 𝐟 𝐬 indicate the surface current density, the surface free charge density and the surface force on the boundary surface between the two mediums, 𝐧 indicates the normal vector directing from the medium '1' to '2'. Under the circumstance of the problems to be studied in this dissertation, all the three surface source terms are zero. In this way the interface continuous conditions demonstrate physically that the magnetic induction (B), the electric displacement (D) and the stress (T) are continuous in the normal direction, and the magnetic field (H), the electric field (E) and the strain (S) are continuous in the tangential direction between the two mediums.

The constitutive laws for the magnetoelectric composite material are the combinations of the electromechanical law in piezoelectric material and the magneto-mechanical law in magnetostrictive material.

The combinations between the mechanical-electro-magnetic parameters of the composite material are illustrated in the Fig. 2.1. The read/purple and blue/yellow arrows are respectively the different possible electro-mechanical and the magneto-mechanical combinations in which the constitutive can be established. The constitutive laws can be written, respectively, in the piezoelectric material [44] 𝑫 = 𝜀 𝑆 𝑬 + 𝑒𝑺, 𝑻 = 𝑐 𝐸 𝑺 -𝑒 𝑡 𝑬 , (2-2-19-a)

and in the magnetostrictive material [START_REF] Engdahl | Handbook of Giant Magnetostrictive Materials[END_REF]:

𝑯 = 𝑣 𝑆 𝑩 -ℎ𝑺, 𝑻 = 𝑐 𝐵 𝑺 -ℎ 𝑡 𝑩, (2-2-19-b)
where 𝜀 𝑆 and 𝑣 𝑆 = (𝜇 𝑆 ) -1 are respectively the permittivity and the reluctivity under constant strain (noted by the superscription s), 𝑒 the piezoelectric coupling coefficient in C•m -2 , ℎ = 𝑞𝑣 𝑆 where 𝑞 is the piezomagnetic coupling coefficient in N • A -1 m -1 , 𝑐 𝐸 and 𝑐 𝐵 are the piezoelectric and magnetostrictive stiffness tensors under constant electric field (noted by the superscription E) and magnetic induction(noted by the superscription B), respectively.

The constitutive laws can also be rewritten in another form [44,[START_REF] Engdahl | Handbook of Giant Magnetostrictive Materials[END_REF]

𝑫 = 𝜀 𝑇 𝑬 + 𝑑 𝑝 𝑻, 𝑺 = 𝑠 𝐸 𝑻 + 𝑑 𝑝 𝑡 𝑬 (2-2-20-a) 𝑩 = 𝜇 𝑇 𝑯 + 𝑑 𝑚 𝑻, 𝑺 = 𝑠 𝐻 𝑻 + 𝑑 𝑚 𝑡 𝑯 (2-2-20-b)
where 𝜀 𝑇 and 𝑣 𝑇 = 𝜇 𝑇 -1 are respectively the permittivity and the reluctivity under constant stress, 𝑑 𝑝 = 𝑒(𝑐 𝐸 ) -1 and 𝑑 𝑚 = 𝑞(𝑐 𝐻 ) -1 are the piezoelectric andpiezomagnetic coupling in C/N, 𝑠 = 𝑐 -1 is the compliance tensor, and 𝑠 𝐸 = (𝑐 𝐸 ) -1 and 𝑠 𝐵 = (𝑐 𝐵 ) -1 represent the tensors under constant electric field and magnetic induction, respectively.

The relations between the expressions under constant strain and constant stress are given with the coupling factor 𝑘 𝑝 and 𝑘 𝑚 [44,[START_REF] Engdahl | Handbook of Giant Magnetostrictive Materials[END_REF]:

𝜀 𝑆 = 𝜀 𝑇 -𝑑 𝑝 𝑐 𝐸 𝑑 𝑝 𝑡 = 𝜀 𝑇 (1 - 𝑑 𝑝 𝑐 𝐸 𝑑 𝑝 𝑡 𝜀 𝑇 ) = 𝜀 𝑇 (1 -𝑘 𝑝 2 ) (2-2-21-a)
with 𝑘 𝑝 = √ 𝑑 𝑝 𝑐 𝐸 𝑑 𝑝 𝑡 𝜀 𝑇

and

𝜇 𝑆 = 𝜇 𝑇 -𝑑 𝑚 𝑐 𝐻 𝑑 𝑚 𝑡 = 𝜇 𝑇 (1 - 𝑑 𝑚 𝑐 𝐻 𝑑 𝑚 𝑡 𝜇 𝑇 ) = 𝜇 𝑇 (1 -𝑘 𝑚 2 ) (2-2-21-b) with 𝑘 𝑚 = √ 𝑑 𝑚 𝑐 𝐻 𝑑 𝑚 𝑡 𝜇 𝑇 and c 𝐻 = c 𝐵 -𝑞 𝑣 𝑆 𝑞 𝑡
In orthotropic material case and in accordance with the local coordinates of the magnetostrictive and piezoelectric materials, the following matrix forms are obtained [44,[START_REF] Engdahl | Handbook of Giant Magnetostrictive Materials[END_REF]

𝑑 33 𝑝 0 0 0 ] [ 𝐸 1 𝐸 2 𝐸 3 ] (2-2-22) [ 𝐷 1 𝐷 2 𝐷 3 ] = [ 0 0 0 0 𝑑 15 𝑝 0 0 0 0 𝑑 24 𝑝 0 0 𝑑 31 𝑝 𝑑 31 𝑝 𝑑 33 𝑝 0 0 0 ] [ 𝑇 1 𝑇 2 𝑇 3 𝑇 4 𝑇 5 𝑇 6 ] 𝑝 + [ 𝜀 11 𝑇 0 0 0 𝜀 11 𝑇 0 0 0 𝜀 33 𝑇 ] [ 𝐸 1 𝐸 2 𝐸 3 ] (2-2-23)
and for the magnetostrictive materials:

[ 𝑆 1 𝑆 2 𝑆 3 𝑆 4 𝑆 5 𝑆 6 ] 𝑚 = [ 𝑠 11 𝑚 𝑠 12 𝑚 𝑠 12 𝑚 𝑠 22 𝑚 𝑠 13 𝑚 0 𝑠 23 𝑚 0 0 0 0 0 𝑠 13 𝑚 𝑠 23 𝑚 0 0 𝑠 33 𝑚 0 0 𝑠 44 𝑚 0 0 0 0 0 0 0 0 0 0 0 0 𝑠 55 𝑚 0 0 𝑠 66 𝑚 ] [ 𝑇 1 𝑇 2 𝑇 3 𝑇 4 𝑇 5 𝑇 6 ] 𝑚 + [ 0 0 𝑑 31 𝑚 0 0 𝑑 31 𝑚 0 0 𝑑 15 𝑚 0 0 𝑑 24 𝑚 0 0 𝑑 33 𝑚 0 0 0 ] [ 𝐻 1 𝐻 2 𝐻 3 ] (2-2-24) [ 𝐵 1 𝐵 2 𝐵 3 ] = [ 0 0 0 0 𝑑 15 𝑚 0 0 0 0 𝑑 24 𝑚 0 0 𝑑 31 𝑚 𝑑 31 𝑚 𝑑 33 𝑚 0 0 0 ] [ 𝑇 1 𝑇 2 𝑇 3 𝑇 4 𝑇 5 𝑇 6 ] 𝑚 + [ 𝜇 11 𝑇 0 0 0 𝜇 11 𝑇 0 0 0 𝜇 33 𝑇 ] [ 𝐻 1 𝐻 2 𝐻 3 ] (2-2-25)
where the superscript p and m stand for the piezoelectric material and the magnetostrictive material, re spectively.

The constitutive relations will be employed in the following modelling as fundamental formulations. T he coupling piezomagnetic and piezoelectric coefficients will be treated as constant in the linear cases and some as parameters to be determined in the nonlinear cases. The purpose of the modelling is to fin ally determine the ME coefficient.

Analytical methods

In small signal analysis, the magneteoelectric composites can be analytically modelled by a homogenous approximation either by a resolution with a simplified analytical method or by a magneto-elastic-electric equivalent circuit model.

Simplified analytical method in static regime

The simplified analytical method presented here is a 0D-modelling that has been used in [START_REF] Bayrashev | Low frequency wireless powering of microsystems using piezoelectricmagnetostrictive laminate composites[END_REF] (denoted as the matrix resolution method) to investigate the output voltage of the magnetoelectric composite in LT and TT modes shown in Fig. 2.3. In L-T mode, the magnetostrictive material is magnetized along the longitudinal direction whereas the piezoelectric material is polarized along the transversal direction.

In TT-mode, the magnetostrictive material and piezoelectric material are magnetized and polarized along the transversal direction.

LT-mode TT-mode By convention, the direction 3 is aligned with the excitation magnetic field by adjusting the local coordinates of the materials from Fig. 2.2. Thus, for the mode TT, the following conditions produced by the 0D-modelling magneto-electric assumption (all the magnetic and electric variables in the nonmagnetization or non-polarization directions are neglected) and the open circuit condition (the coupling current through the piezoelectric material along the polarization direction is zero) is imposed:

𝐻 1 = 𝐻 2 = 0, and 𝐷 3 = 0 (2-3-3)
Moreover, the Newton's third law involves:

𝑣𝑇 1 𝑝 = (𝑣 -1)𝑇 1 𝑚 , 𝑣𝑇 2 𝑝 = (𝑣 -1)𝑇 2 𝑚 and 𝑇 3 𝑚 = 𝑇 3 𝑝 = 0 (2-3-4)
where 𝑣 = In assuming both materials expand and contract together, the following conditions are respected:

𝑆 1 𝑚 = 𝑆 1 𝑝 , 𝑆 2 𝑚 = 𝑆 2 𝑝 (2-3-5)
In this way, in combining all conditions the following matrix system is obtained:

𝐾𝑋 = 𝑏 (2-3-6)
where

, 𝑋 = [ 𝑇 1 𝑚 𝑇 2 𝑚 𝐻 3 𝐸 3 ]
and 𝑏 = [

0 0 𝐵 3 0 ]
where the coefficient 𝛽 is defined as 𝛽 = 𝑣 1-𝑣 ..

The obtained solution 𝑋 = 𝐴 -1 𝑌 can be exploited to find the magnetoelectric coefficients as:

𝛼 33 𝐸 = 𝐸 3 𝐻 3 = 𝑋(3) 𝑋(4) = 2𝑑 31 𝑚 𝑑 31 𝑝 (1-𝑣) 2𝑑 31 𝑝 2 (1-𝑣)-𝜀 33 𝑇 (𝑠 11 𝑚 +𝑠 12 𝑚 )𝑣-𝜀 33 𝑇 (𝑠 11 𝑝 +𝑠 12 𝑝 )(1-𝑣) (2-3-7)
which can also be written in terms of the voltage magnetoelectric coefficient:

𝛼 33 𝑉 = 𝑡 𝑝 𝛼 33 𝐸 = 2𝑑 31 𝑚 𝑑 31 𝑝 (1-𝑣)𝑣𝑡 𝑙𝑎𝑚 2(𝑑 31 𝑝 ) 2 (1-𝑣)-𝜀 33 𝑇 (𝑠 11 𝑚 +𝑠 12 𝑚 )𝑣-𝜀 33 𝑇 (𝑠 11 𝑝 +𝑠 12 𝑝 )(1-𝑣) (2-3-8)
These expressions can also be found in [START_REF] Bichurin | Theory of low-frequency magnetoelectric coupling in magnetostrictive-piezoelectric bilayers[END_REF].

The voltage magnetoelectric coefficient 𝛼 33 𝑉 can also be written as (2-3-9) in considering the ratio

𝑛 = 2𝑡 𝑚 𝑡 𝑙𝑎𝑚
, in the trilayer case 𝑣 = 1 -𝑛.

𝛼 33 𝑉 = 2𝑑 31 𝑚 𝑑 31 𝑝 (1-𝑛)𝑛𝑡 𝑙𝑎𝑚 2(𝑑 31 𝑝 ) 2 𝑛-𝜀 33 𝑇 (𝑠 11 𝑚 +𝑠 12 𝑚 )(1-𝑛)-𝜀 33 𝑇 (𝑠 11 𝑝 +𝑠 12 𝑝 )𝑛 (2-3-9)
In employing the same procedure, the voltage magnetoelectric coefficients concerning the others modes (LT, TL and LL) can be obtained in employing the Bond's transformation on the formulations from (2-2-22) to (2-2-25) [START_REF] Bond | The mathematics of the physical properties of crystals[END_REF]. For instance, the matrix forms for the magnetostrictive material part of the LT-mode becomes: (2-3-12)

[ 𝑆 1 𝑆 2 𝑆 3 𝑆 4 𝑆
Figure 2.4 shows the simulation results for a MPM (magnetostrictive-piezoelectric-magnetostrictive) composite composed of Terfenol-D/PZT-5A/Terfenol-D. We can notice that ME voltage coefficients in LT mode is significantly higher relative to the TT mode with a maximum ratio n around 2/3. This difference is explained by the fact that the coefficient is 𝑑 33 𝑚 that is more greater that 𝑑 31 𝑚 . This part has shown that the Simplified Analytical Method in static regime is a quite useful method that gives the analytical expressions of the magnetoelectric coefficients. The ME coupling in ME laminates is strongly related to their working modes in which the value of ME voltage coefficient in LT mode is significantly higher than the TT mode. The maximum coupling coefficient depends on the volume ratio and the composite materials used.

Equivalent circuit method

The equivalent circuit method for the magnetoelectric composite uses an extension of the well-known 1D electrical Mason's model of piezoelectric material in combining the magnetic-mechanical coupling with the mechanical-electrical coupling of the piezoelectric. It was recently employed in literatures [START_REF] Dong | Longitudinal and Transverse Magnetoelectric Voltage Coefficients of Magnetostrictive/Piezoelectric Laminate Composite: Theory[END_REF][START_REF] Yang | Resonant magnetoelectric response of magnetostrictive/piezoelectric laminate composite in consideration of losses[END_REF][START_REF] Dong | Equivalent circuit method for static and dynamic analysis of magnetoelectric laminated composites[END_REF] to investigate the voltage magnetoelectric coefficient for different modes in static and dynamic regimes. The procedures based for the LT and TT modes are presented in Figures 2.6 and 2.7. Considering the 3D composite structure, if the length of the composite bulk is much longer than the width and thickness, we can approximately only take into account the components in z-direction of the mechanical variables(stress and strain) and that of the electric displacement in polarization direction(direction '3' in local coordinates).

If an infinitesimal longitudinal length 'Δz' is taken into consideration that is as much smaller than the total length of the bulk, then the continuous mechanical equilibrium equation in dynamic regime (2-2-10) for the LT mode can be rewritten as:

𝜕 2 𝑢 𝜕𝑡 2 (2∆𝑚 1 + ∆𝑚 2 ) = ∆𝑇 33 𝑚 • 2𝐴 1 + ∆𝑇 11 𝑝 • 𝐴 2 = 𝜕 2 𝑢 𝜕𝑡 2 (2𝜌 𝑚 𝐴 1 ∆𝑧 + 𝜌 𝑝 𝐴 2 ∆𝑧) (2-3-13)
where 𝐴 1 and 𝐴 2 are the cross-sectional areas of the magnetostrictive layer noted as M-layer, same after) and the piezoelectric layer noted as P-layer, same after) respectively, ∆𝑚 1 = 𝜌 𝑚 𝐴 1 ∆𝑧 and ∆𝑚 2 = 𝜌 𝑝 𝐴 2 ∆𝑧 are the volume masses of the M-layer and P-layer respectively, ∆𝑇 33 𝑚 and ∆𝑇 11 𝑝 are the normal stresses by unit along the z-direction of the M-layer and P-layer respectively.

After simplification we obtain:

𝜕 2 𝑢 𝜕𝑡 2 • 2𝜌 𝑚 𝐴 1 +𝜌 𝑝 𝐴 2 𝐴 = 2𝐴 1 𝐴 ∆𝑇 33 𝑚 ∆𝑧 + 𝐴 2 𝐴 ∆𝑇 11 𝑝 ∆𝑧 (2-3-14)
where 𝐴 = 2𝐴 1 + 𝐴 2 is the total cross-sectional area of the composite.

By defining the thickness ratio 𝑛 = 2𝐴 1 /𝐴 and the average density 𝜌̅ = , therefore we obtain the equation by defining k the wave number and 𝜔 the angular frequency of vibration:

𝜕 2 𝑢 𝜕𝑡 2 = 𝑣̅ 2 𝜕 2 𝑢 𝜕𝑧 2 = ( 𝜔 𝑘 ) 2 𝜕 2 𝑢 𝜕𝑧 2 (2-3-16)
The general solution of the equation above is 𝑢(𝑧) = 𝐴𝑐𝑜𝑠(𝑘𝑧) + 𝐵𝑠𝑖𝑛(𝑘𝑧) , and with the limit conditions 𝑢 1 ̇= 𝑢̇(0) = 𝑗𝜔𝑢(0), 𝑢 2 ̇= 𝑢̇(𝐿) = 𝑗𝜔𝑢(𝐿), we obtain [START_REF] Bayrashev | Low frequency wireless powering of microsystems using piezoelectricmagnetostrictive laminate composites[END_REF]:

𝑢(𝑧) = 𝑢 1 j𝜔 𝑐𝑜𝑠(𝑘𝑧) + 𝑢 2 ̇-𝑢 1 ̇𝑐𝑜𝑠 (𝑘𝐿) 𝑗𝜔𝑠𝑖𝑛 (𝑘𝐿) 𝑠𝑖𝑛 (𝑘𝑧) (2-3-17)
where L is the length of the composite.

According to 𝑆(𝑧) = 𝜕𝑢(𝑧) 𝜕𝑧

we can obtain the strains on the boundary of the bulk:

𝑆(0) = 𝑢 2 ̇-𝑢 1 ̇𝑐𝑜𝑠 (𝑘𝐿) 𝑗𝑣 ̅𝑠𝑖𝑛 (𝑘𝐿) , 𝑆(𝐿) = 𝑢 2 ̇𝑐𝑜𝑠(𝑘𝐿)-𝑢 1 j𝑣 ̅𝑠𝑖𝑛 (𝑘𝐿) (2-3-18)
Considering the external applied forces shown in Figure 2.8 along the longitudinal direction, the equilibrium law gives:

𝐹 1 = -2𝐴 1 𝑇 33 𝑚 (0) -𝐴 2 𝑇 11 𝑝 (0), 𝐹 2 = -2𝐴 1 𝑇 33 𝑚 (𝐿) -𝐴 2 𝑇 11 𝑝 (𝐿) (2-3-19)
Figure 2.8 Composite laminates with external applied surface forces on z=0 and z=L [START_REF] Su | Modélisation par éléments finis des matériaux composites magnéto-électriques[END_REF] By introducing the constitutive laws in the strain expressions we conclude that:

𝐹 1 = 𝑗𝐴𝜌̅ 𝑣̅ 𝑡𝑎𝑛 ( 𝑘𝐿 2 ) 𝑢 1 ̇+ 𝐴𝜌 ̅𝑣 ̅ 𝑗𝑠𝑖𝑛(𝑘𝐿) (𝑢 1 ̇-𝑢 2 ̇) + 2𝐴 1 𝑠 33 𝐻 𝑑 33 𝑚 𝐻 3 + 𝑤 𝑠 11 𝐸 𝑑 31 𝑝 𝑉 (2-3-20) 𝐹 2 = -𝑗𝐴𝜌̅ 𝑣̅ 𝑡𝑎𝑛 ( 𝑘𝐿 2 ) 𝑢 2 ̇+ 𝐴𝜌 ̅𝑣 ̅ 𝑗𝑠𝑖𝑛(𝑘𝐿) (𝑢 1 ̇-𝑢 2 ̇) + 2𝐴 1 𝑠 33 𝐻 𝑑 33 𝑚 𝐻 3 + 𝑤 𝑠 11 𝐸 𝑑 31 𝑝 𝑉 (2-3-21)
where w is the width of the bulk and V the voltage between the electrodes with 𝑉 = 𝑡 𝑝 𝐸 3 = 𝐴 2 𝑤 𝐸 3 .

By defining the corresponding coefficients in the expressions of the forces as 𝑍 1 , 𝑍 2 , 𝜑 𝑚 , 𝜑 𝑝 the mechanical impedances, the magnetostrictive coupling coefficient and the piezoelectric coupling coefficient, and we rewrite the expressions as:

𝐹 1 = 𝑍 1 𝑢 1 ̇+ 𝑍 2 (𝑢 1 ̇-𝑢 2 ̇) + 𝜑 𝑚 𝐻 3 + 𝜑 𝑝 𝑉 (2-3-22) 𝐹 2 = -𝑍 1 𝑢 2 ̇+ 𝑍 2 (𝑢 1 ̇-𝑢 2 ̇) + 𝜑 𝑚 𝐻 3 + 𝜑 𝑝 𝑉 (2-3-23)
The coupling current 𝐼 𝑝 produced by the piezoelectric layer can be calculated as [START_REF] Wang | Electrical resistance load effect on magnetoelectric coupling of magnetostrictive piezoelectric laminated composite[END_REF]:

𝐼 𝑝 = 𝑗𝜔𝐶 0 𝑉 -𝜑 𝑝 (𝑢 1 ̇-𝑢 2 ̇) (2-3-24)
where 𝐶 0 = Under free-free boundary condition where 𝐹 1 = 𝐹 2 = 0 the two ends can be shorted to the 'ground' and the equivalent circuit can be simplified as in Figure 2.10: where

𝑍 = 𝑍 1 //𝑍 1 + 𝑍 2 = - 1 2 𝑗𝐴𝜌̅ 𝑣̅ 𝑡𝑎𝑛 ( 𝑘𝐿 2 ) + 𝐴𝜌 ̅𝑣 ̅ 𝑗𝑠𝑖𝑛(𝑘𝐿) = - 1 2 𝑗𝐴𝜌̅ 𝑣̅ 𝑐𝑜𝑡 ( 𝑘𝐿 2 )
In open circuit condition, the induced current 𝐼 𝑝 = 0, thus 𝑗𝜔𝐶 0 𝑉 = 𝜑 𝑝 (𝑢 1 ̇-𝑢 2 ̇), and by corporating the expressions of the mechanical voltage we can finally conclude that

𝛼 31 𝑉 = | 𝑑𝑉 𝑑𝐻 3 | = | 𝜑 𝑚 𝜑 𝑝 𝜑 𝑝 2 +𝑗𝜔𝐶 0 𝑍 | (2-3-25)
As mentioned previously, in practice, due to the mechanical damping and the electric and magnetic dissipation losses, the impedance Z is a RLC circuit such as [START_REF] Dong | Equivalent circuit method for static and dynamic analysis of magnetoelectric laminated composites[END_REF] and 𝑄 𝑚𝑒𝑐ℎ is the mechanical quality factor of the composite, which is an important parameter characterizing the damping effect and bandwidth of the material.

𝑍 = 𝑅 𝑚 + 𝑗𝜔𝐿 𝑚 + 1 𝑗𝜔𝐶 𝑚 , where 𝑅 𝑚 = 𝜋𝑍 0 8𝑄 𝑚𝑒𝑐ℎ , 𝐿 𝑚 = 𝜋𝑍 0 8𝜔 𝑟 , 𝐶 𝑚 = 8 𝜋𝑍 0 𝜔 𝑟 , 𝑍 0 = 𝐴𝜌̅ 𝑣̅ all presented in
The total quality factor 𝑄 𝑚𝑒𝑐ℎ of the composite can be expressed with quality factors 𝑄 𝑚-𝑚𝑒𝑐ℎ and 𝑄 𝑝-𝑚𝑒𝑐ℎ from the magnetostrictive and piezoelectric materials:

𝟏 𝑸 𝒎𝒆𝒄𝒉 = 𝒏 𝑸 𝒎-𝒎𝒆𝒄𝒉 + 𝟏-𝒏 𝑸 𝒑-𝒎𝒆𝒄𝒉 (2-3-26) 𝛼 31 𝑉 = 8𝑄 𝑚𝑒𝑐ℎ 𝜋 2 (𝛼 31 𝑉 ) 𝑙𝑜𝑤 (2-3-27)
where (𝛼 31 𝑉 ) 𝑙𝑜𝑤 signifies the voltage coefficient in low frequency. It can be seen that the two methods display great frequency-dependent result accordance, with a small resonant frequency shift and little higher resonant ME coefficient value for the equivalent circuit method.

In the mechanical resonant frequency regime ( 𝜔 = 𝜔 𝑟 = 𝜋𝑣 ̅ 𝐿 ), the impedance vanishes since 𝑍 = -

1 2 𝑗𝐴𝜌̅ 𝑣̅ 𝑐𝑜𝑡 ( 𝜔 𝑟 𝐿 2𝑣 ̅ ) = - 1 2 𝑗𝐴𝜌̅ 𝑣̅ 𝑐𝑜𝑡 ( 𝜋 2 ) = 0
Thus, in LT-mode the ME coefficient achieves the peak value as This part has shown that the equivalent circuit method is quite useful for both static and dynamic analysis of magnetoelectric (ME) laminates, especially for electromechanical resonance analysis.

(𝛼 𝑣 ) 𝑚𝑎𝑥 = | 𝑑𝑉 𝑑𝐻 3 | = | 𝜑 𝑚 𝜑 𝑝 | = | 2𝑡 𝑚 𝑑
However, the simplified analytical method and equivalent circuit method are not adapted to include the nonlinearity of piezomagnetic coefficients, or to take into account the physical mechanical impact of the resin (like Epoxy) that sticks the magnetrostrictive and piezoelectric layers. In these conditions, the achievement of a rigorous numerical modelling is essential for the design and optimization of ME devices composed of laminated layers.

The constitutive relations of magnetostrictive materials are fundamentally nonlinear and due to these nonlinear material properties, modelling of the system will also become nonlinear. Consequently it is essential to consider a nonlinear model that takes into account the properties of magnetostrictive materials. The first part of the remainder of this chapter start in providing the essential conditions to perform a 2D modelling under the stress plan conditions while taking into account the electrical and magnetic assumptions of a ME composite. The second part introduces the 2D finite element formulation of the field problem that combines the nonlinear material properties of the magnetostrictive layers as well as the electrical circuit load effect.

FEM modelling of the field problem in 2D

Due to the complexity of the partial differential constitutive equations of the studied material and of the geometry, it is necessity to introduce appropriate numerical resolution methods. The finite element method (FEM) is employed to perform the discretization of the solution domain for the studied ME composite problem. In this section we expose the 2D finite element approach to investigate the ME energy transducer presented in Figure 2.12. It is a trilayer magnetostrictive/piezoelectric laminated composite including an electrical load representing the conditioning circuit connected to the electrodes of the piezoelectric layer. The magnetostrictive material is magnetized along the longitudinal direction whereas the piezoelectric material is polarized along the transversal direction (i.e. L-T mode). It is to be noticed that the configuration of figure 2.12 is only given in title of example and the FEM formulation presented in this section is general and valid for any configuration modes and more complicated structures.

Establishment of the 2D tensors

The symmetry of the laminate structures problem satisfies the plane theory of elasticity that reduces the problem in a 2D problem. For that, two conditions are suitable: Either 2D plane stress conditions or 2D plane strain conditions. In the 2D plane stress, the geometry of the structure is essentially that of a plate with one dimension much smaller than the other, whereas for the 2D plane strain conditions the direction (z-coordinate) of the structure in one direction is very large in comparison with in other directions (x and y-coordinate axes). In both cases, the stiffness tensors can be written as : If 𝑚 = 1 the stress plane is respected whereas when 𝑚 = 0 the strain plane is respected.

𝑐 2𝐷 = 𝜙 [ 1 𝜗 1-𝑚𝜗 0 𝜗 1-𝑚𝜗 1 0 0 0 1-𝜗-𝑚𝜗 2(1-𝑚𝜗)]
Here, the modelling has been employed under the 2D plane stress conditions for the composite to be studied is a thin structure as illustrated in Figure 2.13. It can be noticed that the reduced stiffness tensor (2-4-1) is only valid for an isotropic case and cannot to be employed with an orthotropic case. To establish the reduced stiffness tensor of an orthotropic case, consider the following constitutive laws of the piezoelectric with linear constants.

𝑻 = 𝑐 𝐸 𝑺 -𝑒 𝑡 𝑬 (2-4-2) 𝑫 = 𝜀 𝑆 𝑬 + 𝑒𝑺 (2-4-3)
In matrix form, we have: (2-4-4)

[ 𝐷 1 𝐷 2 𝐷 3 ] = [ 0 0 0 0 𝑒 15 0 0 0 0 𝑒 24 0 0 𝑒 31 𝑒 32 𝑒 33 0 0 0 ] [ 𝑆 1 𝑆 2 𝑆 3 𝑆 4 𝑆 5 𝑆 6 ] + [ 𝜀 11 𝑆 0 0 0 𝜀 11 𝑆 0 0 0 𝜀 33 𝑆 ] [ 𝐸 1 𝐸 2 𝐸 3 ] (2-4-5)
For the elasticity part, the 2D plane stress conditions involve according to the global coordinates the following conditions:

𝑇 3 = 𝑇 4 = 𝑇 5 = 0 (2-4-6) 𝑆 4 = 𝑆 5 = 0 (2-4-7) 𝑆 1 = -𝜗(𝑇 1 + 𝑇 2 )/𝔈 (2-4-8)
The problem can be solved in exploiting only the tangential components of the magnetic induction (𝐵 2 ,𝐵 3 ) and the electrical field (𝐸 2 ,𝐸 3 ) in the working plan yOz. Their normal components 𝐵 1 and 𝐸 1

(perpendicular to the working plan) are assumed invariant where here we considered 𝐵 1 = 0 and 𝐸 1 = 0. { 𝑇 1 = 𝑐 where 𝑓 = 𝑒 or 𝑞 for respectively the piezoelectric and magnetostrictive materials.

The permittivity or the permeability is given by:

𝑔 = [ 𝑔 11 𝑈 + 𝑓 13 2 𝑐 33 𝐸 0 0 𝑔 33 𝑈 + 𝑓 33 2 𝑐 33 𝐸 ] = [ 𝑔 11 𝑈 ̅̅̅̅̅ 0 0 𝑔 33 𝑈 ̅̅̅̅̅ ] (2-4-15)
with 𝑔 = 𝜀 or 𝜇.

In general manner, the coefficient elements can be rewritten as:

𝑐 𝑎𝑏 ̅̅̅̅ = 𝑐 𝑎𝑏 -𝑐 𝑎3 𝑐 3𝑏 /𝑐

FEM formulation 2.4.2.1 Fundamental formulations

The finite element formulation (FEM) of a ME composite is derived in combining the elastodynamic and electromagnetic governing physics equations [START_REF] Nguyen | Modélisation par éléments finis de matériaux composites magnéto-électriques[END_REF][START_REF] Nguyen | Finite element modeling of magnetic field sensors based on nonlinear magnetoelectric effect[END_REF][START_REF] Nguyen | Finite Element Harmonic Modeling of Magnetoelectric Effect[END_REF][START_REF] Besbes | A generalized finite element model of magnetostriction phenomena[END_REF][START_REF] Talleb | Finite Element Modeling of a Magnetoelectric Energy Transducer Including the Load Effect[END_REF][START_REF] Talleb | Finite element modeling of magnetoelectric laminate composites in considering nonlinear and load effects for energy harvesting[END_REF] presented in the first part of this chapter and together with the systems (2-4-20) and (2-4-21).

{ div𝑻 + 𝒇 = 𝜌 𝑣 𝜕 𝑡 2 𝒖 curl𝑯 = 𝒋 𝒔 + 𝜕 𝑡 𝑫, div𝑩 = 0 curl𝑬 = -𝜕 𝑡 𝑩, div𝑫 = 𝜌 (2-4-20) { 𝑻 = c 𝑋 𝑺 -e t 𝑬 -h t 𝑩 𝑫 = -e𝑺 + ε 𝑆 𝑬 𝑯 = -h𝑺 + 𝑣 𝑆 𝑩 (2-4-21)
with c 𝑋 = c 𝐸 or c 𝐵 in according to the considered material.

The links between the mechanical strain 𝑺, the electric field 𝑬 and the magnetic induction 𝑩 with the working variables, i.e. the mechanical displacement 𝒖, the electric potential scalar 𝑉 and the magnetic vector potential 𝒂, are respectively gathered by (2-4-22).

{ 𝑬 = -grad𝒗 𝑩 = curl𝒂 𝑺 = 1 2 (grad + grad 𝑇 )𝒖 (2-4-22)
As 𝐵 1 = 0, 𝐵 2 ≠ 0 and 𝐵 3 ≠ 0 just the normal component 𝒂 𝟏 of the magnetic potential vector 𝒂 in the x direction is considered, and the electrostatic and magnetostatic formulations (no free charges 𝜌 = 0, curl𝒆 = 0 and curl𝒉 = 𝒋 𝒔 )are applied. For that, the operator curl is degenerated to the operator gradient by the relation curl=r * . grad, where r * = [ 0 -1 1 0 ] is a rotation matrix.

In considering just the tangential components of the magnetic induction (𝐵 2 ,𝐵 3 ) and the electrical field (𝐸 2 ,𝐸 3 ) in the working plan yOz, we obtain :

{ 𝑬 = 𝐺 𝑣 {𝑉} 𝑩 = 𝐺 𝑎 {𝑎 𝑧 } 𝑺 = 𝐺 𝑢 {𝒖} (2-4-23) { {𝑉} = 𝑁 𝑣 𝑉 {𝑎 𝑧 } = 𝑁 𝑎 𝑎 𝑧 {𝒖} = 𝑁 𝑢 𝒖 (2-4-24)
with 𝐺 𝑣 = grad𝑁 𝑣 , 𝐺 𝑎 = 𝑟 * . grad𝑁 𝑎 and 𝐺 𝑢 = 1 2 (grad𝑁 𝑢 + grad 𝑇 𝑁 𝑢 ), where 𝑁 𝑣 , 𝑁 𝑎 and 𝑁 𝑢 are the shape functions associated with each node which satisfy the properties:

∑ 𝑁 𝑖 𝑚 𝑖=1
= 1, 𝑁 𝑗 (𝜁 𝑖 ) = 𝛿 𝑖𝑗 , 𝑚 is the number of nodes by mesh element, 𝜁 𝑖 the natural coordinates and the 𝛿 𝑖𝑗 . Kronecker delta.

In our modelling, linear triangular element (𝑚 =3 nodes) is used for the discretization of all fields. In this case :

𝒖 = [ 𝑢 𝑥 𝑢 𝑦 ] = [ 𝑁 1 0 𝑁 2 0 𝑁 3 0 0 𝑁 1 0 𝑁 2 0 𝑁 3 ] [ 𝑢 𝑥𝑚1 𝑢 𝑦𝑚1 𝑢 𝑥𝑚2 𝑢 𝑦𝑚2 𝑢 𝑥𝑚3 𝑢 𝑦𝑚3] = [𝑁 𝑢 ]{𝒖} (2-4-25) 𝑉 = [𝑁 1 𝑁 2 𝑁 3 ] [ 𝑉 𝑚1 𝑉 𝑚2 𝑉 𝑚3 ] = [𝑁 𝑣 ]{𝑣} (2-4-26) 𝑎 𝑧 = [𝑁 1 𝑁 2 𝑁 3 ] [ 𝑎 𝑚1 𝑎 𝑚2 𝑎 𝑚3 ] = [𝑁 𝑎 ]{𝑎 𝑧 } (2-4-27)
with 𝑁 𝑣 = 𝑁 𝑎 = 𝑁.

The finite element formulation after the FEM discretization of the coupling system equation is:

𝜕 𝑡 2 [ℳ]{𝒳} + 𝜕 𝑡 [𝒞]{𝒳} + [𝒦]{𝒳} = {ℱ} (2-4-28)
where The parameters β (in s) and α (in s -1 ) are Rayleigh's damping coefficients. Generally, in a mechanical oscillation system, the resonant angular frequency with damping effect 𝜔 𝑟 and the frequency without damping effect 𝜔 𝑛 have relation as:

𝜔 𝑟 = 𝜔 𝑛 √1 -𝜍 2 (2-4-29)
where 𝜍 is called the damping coefficient and 𝑄 𝑚𝑒𝑐ℎ = 1/(2 𝜍),

The parameters β and α depend on the energy dissipation characteristic of the structure. Currently, these values cannot be obtained through direct calculation and they must be measured. Then, with the The unknown vector {𝒳} includes the nodal displacement {𝒖} , electrical potential {𝑉} and the magnetic vector potential {𝑎}. Here, no external body force (𝐹 = 0) and, under the open-circuit condition, no electric charge 𝑄 𝑛 brought to the electrodes are considered. Instead of imposing a current density source 𝐽 𝑠 in the excitation vector, the magnetic excitation is implemented in considering non-homogenous Dirichlet conditions on the magnetic vector potential 𝑎 1 in the boundaries of the problem domain Ω.

The elements of each submatrice in (2-4-28) are :

[ 𝐾 𝑢𝑢 𝐾 𝑝𝑝 𝐾 𝑎𝑎 ] = ∑ ∫ [ [𝐺 𝑢 ] 𝑡 𝑐 𝑎 [𝐺 𝑢 ] [𝐺 𝑝 ] 𝑡 𝜀 𝑆 [𝐺 𝑝 ] [𝐺 𝑎 ] 𝑡 𝑣 𝑆 [𝐺 𝑎 ] ] Ω 𝑒 𝑒 𝑑Ω (2-4-32) [ 𝐾 𝑢𝑝 𝐾 𝑢𝑎 ] = ∑ ∫ [ [𝐺 𝑢 ] 𝑡 𝑒 𝑡 [𝐺 𝑝 ] [𝐺 𝑢 ] 𝑡 ℎ 𝑡 [𝐺 𝑎 ] ] Ω 𝑒 𝑒

𝑑Ω

(2-4-33)

Boundary conditions

The magnetic, electric and elastic boundary conditions of the solution domain associated in the resolution of the system in 2D are, respectively, illustrated in Figure 2.16.

Mechanical conditions:

To immobilize the composite ME two fixed mechanical displacements are imposed by the Dirichlet condition u y =u z =0 on the middle (bottom and top) in each magnetostrictive layers.

Magnetic conditions:

The external magnetic excitation source 𝐻 𝑒𝑥𝑡 is composed of a small signal harmonic ac field 𝐻 𝑎𝑐 around a magnetization bias 𝐻 𝑑𝑐 . They are separately modeled in the FEM in static and dynamic regimes. The static regime with the magnetic bias 𝐻 𝑑𝑐 allows determining the incremental characteristics of the nonlinear materials, which will be employed in the dynamic regime. The magnetic excitation is implemented in the excitation vector {ℱ} in considering the non-homogenous Dirichlet conditions, 𝒂(Γ 𝑦𝑧 𝑚𝑎𝑥 ) = a 0 /2 and 𝒂(Γ 𝑦𝑧 𝑚𝑖𝑛 ) = -a 0 /2 as shown in Figure 2 where Ω represents the domain bounded externally by Γ.

In this way, the small signal ac field 𝐻 𝑎𝑐 is:

𝐻 𝑎𝑐 = 𝐵 0 /𝜇 0 = 𝑎 0 /(∆ 𝑧 𝜇 0 ) (2-4-35)
For example, in considering 𝐻 𝑎𝑐 = 79.57 A/m, i.e. 1 Oe and ∆ 𝑧 = 16 𝑚𝑚, we obtain 𝑎 0 = 1.6µT.m.

Electrical conditions:

To guarantee a solution of the electric potential {V} the Dirichlet condition 𝑉(Γ 𝑦𝑧 ) = 0 is applied on the out boundary of the solution domain.

FEM Simulation results

The 2D FEM modelling, presented previously, has been performed in employing the automatic and adaptive meshing of the PDE tool product by MATLAB @ . Figure 2.17

Electric potential scalar Electric field

From the distribution figures above, it can be seen that for the magnetic phenomena, the magnitude of the potential field becomes smaller when approaching near the composite in both modes, and the magnetic field lines are concentrated by the magnetostrictive layers. And for the mechanical phenomena, the general deformation of the composite in the LT mode is about 10 times larger than that in the TT mode in the longitudinal direction and about 3 times larger in the thickness direction. As for the electric potential, the general magnitude is approximately more than 5 times larger in the LT mode than the TT mode, also a reversion of the polarization and electric field directions is observed in the TT mode. In general, we can find that with magnetic excitation of the same value for the composites of the same materials and sizes, the LT mode shows a greater effect on the deformation and output of the composite than the TT mode, as well as a larger ME coefficient.

Nonlinear static case

As previously stated, the magnetoelectric composites operate under a composite external magnetic field excitation: a static biasing field H dc and a small amplitude alternative field H ac . Due to the nonlinear property of magnetostrictive materials, the change of static biasing field allows determining the optimal operation point for which the magnetoelectric coupling coefficient for the small signal field maximizes. In order to determine the physical properties of the materials around the operation point, the so called incremental characteristics, the system equation (2-4-28) needs to be solved in nonlinear regime for static biasing magnetic excitation.

Modelling of nonlinear piezomagnetic coupling

When considering the case of linear elastic property (constant elastic coefficients) and nonlinear magnetic property, the expression of the stress tensor is:

𝑇 = 𝑐 𝐵 (𝑆 -𝑆 𝜇 (𝑩)) = 𝑐 𝐵 𝑆 -𝑇 𝜇 (𝑩) (2-4-36)
where 𝑆 𝜇 (𝑩) and 𝑇 𝜇 (𝑩) are respectively the "coercitive" strain and stress tensor induced by 𝑩 The expression of the "coercitive" strain tensor 𝑆 𝜇 (𝑩) of a single cubic crystal is the Maxwell magnetic stress tensor with six independent components and can be described by a quadratic model of the magnetic induction 𝑩 as [START_REF] Hirsinger | Etude des déformations magnéto-élastiques dans les matériaux ferromagnétiques doux. Application à l'étude des déformations d'une structure de machine électrique[END_REF]:

𝑆 𝜇 (𝑩) = 𝑆 𝑠𝑎𝑡 𝜇 ( 𝑩 𝑩 𝒔𝒂𝒕 ) 2 (2-4-37)
Assuming homogeneous material properties and quasi-stationary conditions, the total magnetic force acting on a body with volume V enclosed by a surface S can be written as:

𝐹 𝐵 = ∭ 𝑓 𝐵 𝑑𝑉 = ∭ 𝐽 × 𝐵𝑑𝑉 = ∭ rot𝐻 × 𝐵𝑑𝑉 = ∬ [𝑯(𝑯. n) -𝜇 1 2 𝐻 2 n] 𝑑𝑉 (2-4-38)
It can also be rewritten as the divergence of the coercitive stress tensor:

𝐹 𝐵 = ∭ 𝑑𝑖𝑣𝑇 𝜇 𝑑𝑉 = ∬ 𝑇 𝜇 . 𝑛𝑑𝑆 (2-4-39)
Comparing both expressions, we obtain:

𝑇 𝜇 (𝑯) = [ 1 2 (𝐻 1 2 -𝐻 2 2 -𝐻 3 2 ) 𝐻 1 𝐻 2 𝐻 1 𝐻 3 𝐻 1 𝐻 2 1 2 (𝐻 2 2 -𝐻 3 2 -𝐻 1 2 ) 𝐻 2 𝐻 3 𝐻 1 𝐻 3 𝐻 2 𝐻 3 1 2 (𝐻 3 2 -𝐻 1 2 -𝐻 2 2 ) ] (2-4-40)
As shown by Hirsinger [START_REF] Buiron | A multiscale model for magneto-elastic couplings[END_REF][START_REF] Hirsinger | Etude des déformations magnéto-élastiques dans les matériaux ferromagnétiques doux. Application à l'étude des déformations d'une structure de machine électrique[END_REF], this tensor can be reduced with 3 independent components when the direction of the magnetic induction 𝑩 is collinear to the easy magnetisation axe. In addition, the isochoric principle ((tr(𝑇 𝜇 (𝑩)) = 0) is respected.

Hirsinger model

The reduced quadratic model proposed by Hirsinger [START_REF] Buiron | A multiscale model for magneto-elastic couplings[END_REF][START_REF] Hirsinger | Etude des déformations magnéto-élastiques dans les matériaux ferromagnétiques doux. Application à l'étude des déformations d'une structure de machine électrique[END_REF] is given by:

𝑆 𝜇 (𝑩) = [ 𝑠 // 𝜇 (𝑩) 0 0 0 𝑠 ⊥ 1 𝜇 (𝑩) 0 0 0 𝑠 ⊥ 2 𝜇 (𝑩) ] = 𝛽 0 𝑩 2 [ 1 0 0 0 1 2 0 0 0 1 2 ] (2-4-41)
where 𝑠 // 𝜇 (𝑩) and 𝑠 ⊥ 𝜇 (𝑩) are the magnetostrictive strain saturation in the parallel and orthogonal directions of the magnetic induction 𝑩. 𝛽 0 is a coefficient determined by experimentations.

In employing a coordinate transformation from the special coordinate into a global orthogonal one, the constitutive relation in Einstein notation can be rewritten as in dependence of 𝑩:

𝑠 𝑘𝑙 𝜇 = 𝛽 0 2 (3𝐵 𝑘 𝐵 𝑙 -𝛿 𝑘𝑙 ‖𝑩‖ 2 ) (2-4-42)
Recent works in the GeePs Lab proposed to modify this constitutive relation either in using higher polynomial functions

1 2 ∑ 𝛽 𝑛 𝑁 𝑛=0
instead of 𝛽 0 [START_REF] Azoum | Contribution à la modélisation numérique de phénomènes magnéto-élastiques. Étude de dispositifs à base de matériaux magnétostrictifs[END_REF][START_REF] Galopin | Modélisation et caractérisation de matériaux actifs pour la conception de dispositifs magnéto-électriques[END_REF][START_REF] Galopin | Finite element modelling of magneto electric sensors[END_REF] or in considering that 𝑩 and the magnetisation 𝑴 are not collinear [START_REF] Nguyen | Modélisation par éléments finis de matériaux composites magnéto-électriques[END_REF][START_REF] Nguyen | Finite element modeling of magnetic field sensors based on nonlinear magnetoelectric effect[END_REF].

𝑠 𝑘𝑙 𝜇 = 1 2 ∑ 𝛽 𝑛 𝑁 𝑛=0 (3𝐵 𝑘 𝐵 𝑙 -𝛿 𝑘𝑙 ‖𝑩‖ 2 ) (2-4-43) 𝑠 𝑘𝑙 𝜇 = 𝛽 0 2𝜇 𝑜 2 (3𝐵 𝑘 𝐵 𝑙 -𝛿 𝑘𝑙 𝑏 𝑖 𝑏 𝑖 ) ‖𝑴‖ 2 ‖𝑩‖ 2
(2-4-44)

In our study, only the constitutive relation (2-4-43) is used.

Effective piezomagnetic coefficients in isotropic case

For convenience and as proposed in [START_REF] Nguyen | Modélisation par éléments finis de matériaux composites magnéto-électriques[END_REF][START_REF] Azoum | Contribution à la modélisation numérique de phénomènes magnéto-élastiques. Étude de dispositifs à base de matériaux magnétostrictifs[END_REF][START_REF] Galopin | Modélisation et caractérisation de matériaux actifs pour la conception de dispositifs magnéto-électriques[END_REF], first we study to obtain the expression of the "coercitive" stress in isotropic case. Consider the Hook's law expressed as:

𝑇 𝑖𝑗 = 𝜆(𝑆 𝑘𝑘 -𝑆 𝑘𝑘 𝜇 (𝐵))𝛿 𝑖𝑗 + 2𝜇 * (𝑆 𝑗𝑖 -𝑆 𝑖𝑗 𝜇 (𝐵)) (2-4-45)
Where

𝜇 * = 𝔈 2(1+𝜗) or 𝜇 * = 𝑐 66 𝐵 = 𝑐 11 𝐵 -𝑐 12 𝐵 2
is the Lamé coefficient .

As mentioned previously, the magnetostriction phenomenon is assumed to be isochoric, i.e. 𝑆 𝑘𝑘 𝜇 = 0 We obtain:

𝑇 𝑖𝑗 = 𝜆𝑆 𝑘𝑘 𝛿 𝑖𝑗 + 2𝜇 * 𝑆 𝑗𝑖 -2𝜇 * 𝑆 𝑖𝑗 𝜇 (𝐵) = 𝑐 𝑗𝑗𝑘𝑙 𝑆 𝑗𝑖 -2𝜇 * 𝑆 𝑖𝑗 𝜇 (𝐵) = 𝑐 𝐵 𝑆 -2𝜇 * 𝑆 𝜇 (𝐵) (2-4-46)
By comparison with (2-4-37) the "coercitive" stress 𝑇 𝜇 (𝐵) is given by :

𝑇 𝜇 (𝐵) = 2𝜇 * 𝑆 𝜇 (𝐵) (2-4-47)
The expression 𝑆 𝜇 (𝐵) can be given with the Hirsinger model in which the component 𝐵 1 = 0

𝑇 𝜇 = [ 𝑇 1 𝜇 𝑇 2 𝜇 𝑇 3 𝜇 𝑇 4 𝜇 𝑇 5 𝜇 𝑇 6 𝜇 ] = 2𝜇 * 𝛽 𝑜 [ 𝐵 2 2 - 𝐵 3 2 2 𝐵 3 2 - 𝐵 2 2 2 - 𝐵 2 2 +𝐵 3 2 2 3 2 𝐵 2 𝐵 3 0 0 ] (2-4-48)
where 𝛽 0 is deduced from experimental results .

In the linear case we have 𝑻 = 𝑐 𝐵 𝑺 -ℎ 𝑡 𝑩. Substituting ℎ 𝑡 𝑩 by 𝑇 𝜇 (𝐵) of the expression we have (2-4-49)

where the nonlinear effective piezomagnetic coefficients are function of 𝑩.

𝑞 33 = 2𝜇 * 𝛽 0 (𝐵 2 2 - 𝐵 3 2
2 )

𝜈 33 𝑆 𝐵 2

(2-4-50)

𝑞 31 = 2𝜇 * 𝛽 0 (𝐵 3 2 - 𝐵 2 2
2 )

𝜈 33 𝑆 𝐵 2
(2-4-51)

𝑞 24 = 𝜇 * 𝛽 0 3𝐵 3 𝜈 11 𝑆
(2-4-52)

Effective piezomagnetic coefficients in orthotropic case

To obtain the expression of the "coercive" stress in orthotropic case, we consider the Hook's law expressed in matrix form with the Hirsinger model in which the component 𝐵 1 = 0 in the 2D case. In this case, we have: And we can obtain:

𝑇 𝜇 = [ 𝑇 1
𝑇 1 𝜇 = 𝛽 0 [𝑐 33 𝐵 (𝐵 3 2 - 𝐵 2 2
2 ) + 𝑐 In this part it was shown that the nonlinear magnetostritive behaviour established by the piezomagnetic coupling could be established in considering the nonlinear magnetic property. To take into account the nonlinear behaviours in the FEM code, it essential to use a magnetic B-H, or M-H curve that is usually obtained by an experimental result or by analytical model such as the Jile-Atherton model . This relation between the piezomagnetic coupling and the magnetic curve could be achieved with the expression of the reluctivity. The latter can be expressed by different analytical models such as the Brauer Model [START_REF] Talleb | Finite Element Modeling of a Magnetoelectric Energy Transducer Including the Load Effect[END_REF] or the exponential series model [START_REF] Hwang | Exponential series for B/H curve modeling[END_REF]. The magneto-elastic hysteresis loops behaviour can be modelled with the Jile-Atherton model [START_REF] Jiles | Theory of ferromagnetic hysteresis[END_REF].

Modelling of magnetic nonlinearity

The Jile-Atherton model for a magnetostrictive material

Introduced in 1984, the Jile-Atherton model describes the hysteretic loops of magnetization of paramagnetic materials as an anhysteretic magnetization in the ferromagnetic materials according to the physical parameters of the magnetic materials [START_REF] Jiles | Theory of ferromagnetic hysteresis[END_REF].

The effective magnetic field influencing on magnetic moments within the material can be described by the following equation:

𝑯 𝒆 = 𝑯 + 𝛼𝑴 (2-4-60)
where 𝑯 represents the applied magnetic field to the material, 𝑴 the total magnetization and 𝛼 is the m ean field parameter representing inter domain coupling according to the Bloch model.

As shown in [START_REF] Zheng | A nonlinear constitutive model of magneto-thermo-mechanical coupling for giant magnetostrictive materials[END_REF][START_REF] Jin | A nonlinear magneto-thermo-elastic coupled hysteretic constitutive model for magnetostrictive alloys[END_REF][START_REF] Marcelo | Coupled Structural-Magnetic Strain and Stress Model for Magnetostrictive Transducers[END_REF][START_REF] Calkins | Energy-based hysteresis model for magnetostrictive transducers[END_REF], for magnetostrictive materials such as the Terfenol-D, the effective magnetic field must be completed by a magnetostrictive parameter 𝛼 𝑀 . from the magnetomechanical stress anisotropies.

𝑯 𝒆 = 𝑯 + (𝛼 + 𝛼 𝑀 )𝑴 = 𝑯 + 𝛼 ̃𝑴 (2-4-61) with 𝛼 ̃= 𝛼 + 𝛼 𝑀 .

In considering that the magnetostrictive material is subjected to an axially applied magnetic field, the parameter 𝛼 𝑀 in 1-D assumption can be expressed [START_REF] Zheng | A nonlinear constitutive model of magneto-thermo-mechanical coupling for giant magnetostrictive materials[END_REF][START_REF] Jin | A nonlinear magneto-thermo-elastic coupled hysteretic constitutive model for magnetostrictive alloys[END_REF] in ambient temperature as: where 𝜆 𝑠 and 𝑀 𝑠 are, respectively, the maximum magnetostriction and the saturation magnetisation. 𝑻 𝒔 is a reference stress and 𝑻 𝒐 is the applied pre-stress.

𝛼 𝑀 = { 2𝑻 𝒐 -
The different derivation processes to get the analytical formulations of the nonlinear constitutive models are given in Appendix B.

The nonlinear magnetization 𝑴 is decomposed of its reversible component 𝑴 𝒓𝒆𝒗 and its irreversible c omponent 𝑴 𝒊𝒓𝒓 :

𝑴 = 𝑴 𝒓𝒆𝒗 + 𝑴 𝒊𝒓𝒓 (2-4-63)
The relationship between these two components and the anhysteretic magnetization 𝑴 𝒂𝒏 is given by :

𝑴 𝒓𝒆𝒗 = 𝑐(𝑴 𝒂𝒏 -𝑴 𝒊𝒓𝒓 ) (2-4-64-a) In other words, 𝑴 = 𝑐𝑴 𝒂𝒏 + (1 -𝑐)𝑴 𝒊𝒓𝒓 (2-4-64-b)
where 𝑐 is a reversibility coefficient to be determined by measurements.

The anhysteretic magnetization 𝑴 𝒂𝒏 can be described with the Langevin function:

𝑴 𝒂𝒏 = 𝑴 𝒔 • (coth ( 𝑯 𝒆 𝑎 ) - 𝑎 𝑯 𝒆 ) (2-4-65)
where 𝑎 = 𝑁𝑘 𝐵 𝑇 𝜇 𝑜 𝑀 𝑠 , 𝑘 𝐵 is the Boltzmann's constant, 𝑇 the ambient temperature, 𝑀 𝑠 the saturation magnetisation and 𝑁 a constant. (2-4-67)

The expression

𝜕𝑴 𝒂𝒏 𝜕𝑯 𝒆 = 𝑴 𝒔 • (1 -coth 2 ( 𝑯 𝒆 𝑎 ) -( 𝑎 𝑯 𝒆 ) 2 ) (2-4-68)
The coefficients 𝛼, 𝑘, 𝛿 are constant parameters of the material to be determined by experiments.

Simulation results for the model

To solve the above ordinary differential equations to obtain the 𝑴(𝑯) dependence, the Runge-Kutta method has been developed and employed. [START_REF] Sandlund | Dynamic characterization of 1-3 type Terfenol-D/epoxy magnetostrictive composite[END_REF] shows the characteristic curves of the magnetisation and the magnetic induction of the mean curve with the applied magnetic field under different pre-stress 𝑻 𝒐 . These curves can be implemented in the FEM model using the incremental procedure presented in the next section. The characteristic behaviour of the piezomagnetic coefficient 𝑑 33 in a 1-D rod structure can be obtained in using the following formulation [START_REF] Zheng | A nonlinear constitutive model of magneto-thermo-mechanical coupling for giant magnetostrictive materials[END_REF][START_REF] Jin | A nonlinear magneto-thermo-elastic coupled hysteretic constitutive model for magnetostrictive alloys[END_REF] in 1-D assumption as (where the formulation of ∂𝑴 ∂𝑯

(2-4-67) is used.): The next section explains the incremental method to solve the nonlinear FEM problem.

𝑑 33 = { [1-tanh(𝐓 𝐨 𝐓 𝐬 ⁄ )]2𝑀𝜆 𝑠 𝑀 𝑠 2 𝜕𝑴 𝜕𝑯 𝑻 𝒐 𝑻 𝒔 ⁄ ≥ 0 [1-

Incremental piezomagnetic coefficients

To solve the finite element equation system (2-4-28) with the nonlinear couplings, it is necessary to introduce and employ some effective numerical method. With the incremental method we need to obtain, the global incremental piezomagnetic coefficients for the nonlinear problem in static case with a piecewise-linear implementation of the constitutive law that can be solved with an iterative method such as the Newton-Raphson process method exposed with details in Appendix C. It is a robust implicit method with a high convergence speed.

With a piecewise-linear implementation, the stress and field increments are related to the nodal degrees of freedoms by: The static problem is solved in considering the incremental system:

∆𝑯 = -
𝐾∆𝒳 = ∆ℱ (2-4-73) [ 𝐾 𝑢𝑢 𝐾 𝑢𝑎 𝐾 𝑎𝑢 𝐾 𝑎𝑎 ] { ∆𝑼 ∆𝑨 } = { 0 ∆𝑱 𝒔 } (2-4-74)
Figure 2.26 shows a the flowchart of the piecewise-linear solution procedure: It should be noted that the material model must be based on the incremental reluctivity and piezomagnetic coefficients.

For the calculation convenience, the following polynomial model of 𝑯 = 𝜈(𝑩)𝑩 is used to fit the curve B-H obtained from the Jile-Atherton model [START_REF] Hwang | Exponential series for B/H curve modeling[END_REF]:

𝑯 = 𝜈(𝑩)𝑩 = 𝑘 𝑜 + ∑ 𝑘 𝑖 exp(𝑡 𝑖 𝑩) 𝑛 𝑖=2
(2-4-75)

With 𝑛 the polynomial degree and the coefficient, 𝑡 𝑖 is calculated according to the points selected in the B-H curve:

𝑡 𝑖 = ln(𝐻 𝑖-1 )-ln(𝐻 𝑖 ) 𝐻 𝑖-1 -𝐻 𝑖 (2-4-76)
The coefficient 𝑘 𝑖 can be determined in considering the matrix systems:

[ 1 exp(𝑡 2 𝐵) ⋯ exp(𝑡 2 𝐵) ⋮ 1 ⋮ exp(𝑡 𝑛 𝐵) ⋮ exp(𝑡 𝑛 𝐵) ] [ 𝑘 𝑜 ⋮ 𝑘 𝑛 ] = [ 𝜈(𝐵 1 ) ⋮ 𝜈(𝐵 𝑛 ) ] (2-4-77) [ 𝑘 𝑜 ⋮ 𝑘 𝑛 ] = [ 1 exp(𝑡 2 𝐵) ⋯ exp(𝑡 2 𝐵) ⋮ 1 ⋮ exp(𝑡 𝑛 𝐵) ⋮ exp(𝑡 𝑛 𝐵) ] -1 [ 𝜈(𝐵 1 ) ⋮ 𝜈(𝐵 𝑛 ) ]
(2-4-78)

The finite element analysis with the Newton-Raphson process method requires an equation of 𝜈(𝑩)

and an equation for the derivative 𝑑𝑯 𝑑𝑩 as function of B to calculate the tangent stiffness matrix.

For that, we have: The incremental piezomagnetic coefficients are obtained in considering the isotropic model [START_REF] Nguyen | Modélisation par éléments finis de matériaux composites magnéto-électriques[END_REF][START_REF] Galopin | Modélisation et caractérisation de matériaux actifs pour la conception de dispositifs magnéto-électriques[END_REF], thus:

𝜈(𝑩) =
ℎ 𝑒 𝑡 = 𝜕𝑇 𝑘𝑙 𝜇 (𝐵) 𝜕𝐵 𝑖 = 3𝛽 𝑜 (𝛿 𝑘𝑖 𝐵 𝑙 + 𝛿 𝑖𝑙 𝐵 𝑘 ) -2𝛿 𝑘𝑖 𝐵 𝑖 (2-4-82)
where 𝑇 𝑘𝑙 𝜇 (𝐵) = 2𝜇 * 𝑆 𝑘𝑙 𝜇 (𝐵) = 𝜇 * 𝛽 0 (3𝐵 𝑘 𝐵 𝑙 -𝛿 𝑘𝑙 ‖𝑩‖ 2 )

In the matrix form, we obtain with the considered coordinate the incremental piezomagnetic coefficients: The elementary stiffness matrix c 𝑒 𝐵 is then calculated with:

ℎ 𝑒 𝑡 = ∆𝑇 𝑘𝑙 𝜇 (𝐵) ∆𝐵 = 𝜕𝑇 𝑘𝑙 𝜇 (𝐵) 𝜕𝐵 𝑖 = 𝜇 * 𝛽 0 [ 2𝐵 3 -𝐵 2 -𝐵 3 2𝐵 2 - 3 2 𝐵 2 - 3 2 𝐵 3 ] ( 2
c 𝑒 𝐵 = ∂𝑻 ∂𝑺 = 𝜕𝑻 𝜕𝑩 𝜕𝑩 𝜕𝑺
(2-4-85)

Figure 2.28 shows the piezomagnetic coefficients obtained with the FEM method through the nonlinear model for a trilayer Terfenol-D/PZT-5A/ Terfenol-D in LT-mode. With the employed model, the optimal dc bias is around 100 Oe.

q coefficients (Cm -2 ) d=c -1 q coefficients (nm/A) Figure 2.28 Piezomagnetic coefficients under the applied magnetic field in Oe unity

Dynamic small signal regime

As illustrated in Figure 2.29, for dynamic signal ac field 𝐻 𝑎𝑐 around a magnetization bias 𝐻 𝑑𝑐 , the system can be solved linearly around an external dc bias point. The magnetostrictive coefficients as well as the reluctivity are incremental values and determined according to the operation point obtained in the nonlinear static FEM analysis. In this way, the system in harmonic regime (𝜕 𝑡 → 𝑗𝜔) can be solved in terms of complex variables as:

[𝐾 ̃]{𝒳 ̅ } = {ℱ ̅ } (2-4-86)
where [𝐾 ̃] denotes the small signal linear matrix around a magnetization point.

[

𝐾 𝑢𝑢 -𝜔 2 𝑀 + 𝑗𝜔𝐶 𝑢𝑢 𝐾 𝑢𝑝 𝐾 𝑢𝑎 𝐾 𝑝𝑢 𝐾 𝑝𝑝 0 𝐾 𝑎𝑢 0 𝐾 𝑎𝑎 ] [ {𝑢} {𝑉} {𝑎} ] = [ 0 0 𝑎 𝑜 ] (2-4-87)
It is well known that the performances of a piezoelectric material depend of the impedance electrical load connected to the electrodes. The next section explains how to take into account the effect of the impedance electrical load in the FEM model.

Coupling with the electric circuit load equation

Since the piezoelectric layer is dielectric, there is no free charge inside it. All the free charges are confined to the electrodes. The current I produced and flowing toward the electrical impedance load is the time derivative of the total charge {Q}, 𝐼 = 𝜕 𝑡 {𝑄}. Thus, this last can be directly obtained with the discretization of the circuit equation expressed by:

{𝑉}𝐾 𝑝𝑞 -𝑍𝜕 𝑡 {𝑄} = 0 (2-4-88)
where 𝐾 𝑝𝑞 is an incident vector in which the elements are respectively equal to 1 or -1 depending on whether the related node is associated to the top electrode or the bottom one, otherwise with the value 0. In this process, the two electrodes (top and bottom) have been considered, respectively, as equipotential and the incident vector 𝐾 𝑝𝑞 means that the nodes associated to each electrode are assigned with the same degrees of freedom. The electric charge 𝑄 on the electrodes becomes an unknown of the problem.

In 2D case considered here, the width 𝐿 𝑥 in the x-direction is not taking into account. Actually, this width is considered as normalized in according to [1m]. The current flowing toward the electrical impedance 𝐼 is thus proportional of the surface 𝒮 = 𝐿 𝑦 × [1𝑚] of the electrodes since 𝑄 = 𝐶𝑉, where 𝐶 is the capacitance between the electrodes of the piezoelectric is given by :

𝐶 = 𝜀 𝑆 𝒮 𝑡 𝑝 (2-4-89)
For a considered distance 𝐿 𝑥 , the capacitance 𝐶 is then normalized as : (2-4-93)

𝐶 ̅ = 𝐶 𝐿 𝑥 = 𝜀 𝑆 𝒮 𝑡 𝑝 𝐿 𝑥
The elements of 𝐾 𝑞𝑝 and 𝐾 𝑝𝑞 are respectively equal to 1 or -1 if a node is associated with the top electrode or the bottom electrode, otherwise its elements equal to 0.

Effect of the complex impedance on the damping losses

Consider the first constitutive law of (2-4-21) in which the electric field 𝑫 is substituted by displacement field 𝑫:

𝑻 = 𝑐 𝑋 𝑺 -𝑒 𝑡 (𝜀 𝑆 )𝑫 -ℎ 𝑡 𝑩 (2-4-94)
The displacement field 𝑫 is given by the Gauss law:

∬ 𝑫𝑑𝑆 = 𝑄 Ω (2-4-95)
In harmonic regime the current 𝐼 produced and flowing toward the electrical impedance load is given by: 𝐼 = 𝑗𝜔𝑄 (2-4-96) In this way, we have:

𝑫 = 𝐼 𝑗𝜔𝒮 = 𝑉 𝑗𝜔𝑍𝒮 (2-4-97)
where 𝒮 represents the electrode surface defined. Thus, the constitutive law becomes

𝑻 ̃= 𝑐 𝑎 𝑺 ̃-𝑒 𝑡 𝑉 𝜀 𝑆 𝑗𝜔𝑍𝒮 -ℎ 𝑡 𝑩 ̃= 𝑐 𝑎 𝑺 ̃-ℎ 𝑡 𝑩 ̃ (2-4-98)
where the strain can decomposed in real and imaginary parts:

𝑺 ̃= 𝑺 ̃-𝑒 𝑡 𝑉 𝜀 𝑆 𝑗𝜔𝑍𝒮 = 𝑺 ̃𝒓 + 𝑗𝑺 ̃𝒊 (2-4-99)

Simulation results

Dynamic response under open circuit condition

Conclusion

This chapter has presented the modelling using analytical methods (simplified analytical resolution using the 'matrix' in 0D assumptions and equivalent circuit method in 1D) with the homogenous quasi-static approximation and the second group and using a numerical method based on the finite element method (FEM) in 2D. Unlike analytical methods, the FEM code takes into account the magnetic nonlinearity of the magnetostrictive material as well as the electrical impact when the structure is loaded by an electrical impedance.

Chapter 3. Assessment of ME composite performances 3.1 Introduction

In Chapter 2, it has been shown that a ME composite bulk composed of Terfenol-D/PZT-5A/Terfenol-D materials can deliver around 10 volts in LT-mode to 50 volts in LL-mode for a volume of 0.252 cm3 (with the length of 14mm, the thickness of 3mm and the width of 6mm) under an open circuit condition. This puts in evidence the potential of ME composites to work as energy transducer for energy harvesting or wireless powering for autonomous micro embedded systems. Thus, in this chapter, we examines the ME composites as energy transducer in the small signal harmonic regime.

Performances of a ME composite as energy transducer

The energy harvesting process for a ME energy transducer is illustrated in For most practical WSN applications, the use of a small duty cycle of the sensor node's operation is a common method to reduce its power consumption and therefore extend the lifetime of the WSN. During each sleeping time 𝑇 𝑠𝑙𝑒𝑒𝑝 , the energy load process in the storage stage is applied to obtain the necessary required voltage to active the WSN.

In fact, the optimal impedance depends on the magneto-electro-mechanical coupling factor of the ME composite transducer. Usually, an Impedance Matching Circuit (IMC) such as Boost Converter (composed of MOSFET, diode and inductor) must be placed between the transducer and the rectified diode bridge for impedance matching in order to maximize energy harvested under different operating conditions. Obviously, the addition of an IMC increases the current consumption, i.e. increase the average power consumption. In considering a consumption of 15 mA when an IMC is activated, the Table 3.2 gives the corrected average power for the three different duty cycles in considering in the equation (3-2-1) a total current of 40 mA (25 mA for the WSN and 15 mA for the IMC) during the activation mode.

Duty

Output deliverable power under different modes

To demonstrate that a ME composite can work as an energy transducer by wireless powering; the power requirements in according to the different modes of the ME composite (LL, TL, LT, and TT) must be determined. Thus, this part shows the mode configuration that supply the maximal active output deliverable power in considering the composite Terfenol-D/PZT-5A/Terfenol-D investigated in Chapter 2. Thus, it appears that the adequate solution for a duty cycle of 1% is the LT-mode since as shown in Figure 3.4 (a) and (b) the deliverable power reaches 800 µW (>500 µW) for optimal impedance close to 18 kOhm.

Here, we consider the ideal case where the impedance matching would be realized under a duty cycle fixed at 1 %. In this case and in according to the Table 3.2, with a potential deliverable power of 800 µW, the ME composite in LT-mode remains a possible solution to active the WSM + IMC.

Electrical equivalent circuit model

In this part we propose an extracted electric elements method to establish an electrical equivalent circuit model thank to the FEM results. This model can be used to estimate the effect of the electric on the output deliverable power or for comparisons with experimental results from an impedance analyzer or an oscilloscope.

As specified in Chapter 2, a magneto-elastic-electric equivalent circuit such as that in Figure 3.5-a or a simplified circuit as in Figure 3.

𝑍 𝑇 = 𝑍 𝐶 𝑜 //𝑍 𝑚 ′ = 𝑍 𝐶 𝑜 𝑍 𝑚 ′ 𝑍 𝐶 𝑜 +𝑍 𝑚 ′ (3-2-2)
with 𝑍 𝑚 ′ = 𝑅 𝑚 ′ + 𝑗𝐿 𝑚 ′𝜔 + we can write:

(1 -

𝑉 𝑜𝑢𝑡 𝑉 𝑖𝑛 ) = 𝑍 𝑚 ′ 𝑍 𝐶 𝑜 +𝑍 𝑚 ′ (3-2-3)
Furthermore 𝑉 𝑖𝑛 = 𝛾𝐻 𝑎𝑐 and 𝛼 ̃𝑉 = 𝑉 𝑜𝑢𝑡 /𝐻 𝑎𝑐 , that involves:

𝑉 𝑜𝑢𝑡 𝑉 𝑖𝑛 = 𝛼 ̃𝑉 𝛾 (3-2-4)
Thus, we can rewrite the internal impedance 𝑍 𝑇 as:

𝑍 𝑇 = (1 - 𝛼 ̃𝑉 𝛾 ) 𝑍 𝐶 𝑜 = 𝑅 𝑇 + 𝑗𝑋 𝑇 (3-2-5)
where the real part 𝑅 𝑇 and the imaginary part 𝑋 𝑇 are both function of the frequency. The change of the amplitude of the internal impedance with the frequency is illustrated in Figure 3.6 obtained in combining the FEM result with the relation (3-2-5). The frequencies 𝑓 1 and 𝑓 2 are the wellknown resonance and the antiresonance frequencies existing into a piezoelectric layer. They can respectively be expressed as:

𝑓 2 = 1 2𝜋√𝐿 𝑚 ′𝐶 𝑚 ′ (3-2-6) 𝑓 1 2 = 𝑓 2 2 (1 + 𝐶 𝑚 ′ 𝐶 𝑜 ) (3-2-7)
In this way, the electrical parameter 𝐿 𝑚 ′ and 𝐶 𝑚 ′ can be extracted:

𝐶 𝑚 ′ = 𝐶 𝑜 ( 𝑓 1 2 𝑓 2 2 -1) (3-2-8) 𝐿 𝑚 ′ = 1 4𝜋 2 𝑓 2 2 𝐶 𝑚 ′ (3-2-9)
The electrical parameter 𝑅 𝑚 ′ can be extracted in considering that

𝑍 𝑚 ′ (𝜔 2 ) = 𝑅 𝑚 ′ 𝑍 𝑇 (𝜔 2 ) = 𝑍 𝐶 𝑜 (𝜔 2 )//𝑅 𝑚 ′ (3-2-10) 𝑌 𝑇 (𝜔 2 ) = 1 𝑍 𝑇 (𝜔 2 ) = 1 𝑍 𝐶 𝑜 (𝜔 2 ) + 1 𝑅 𝑚 ′ (3-2-11) 𝑅 𝑚 ′ = 1 𝑟𝑒𝑎𝑙(𝑌 𝑇 (𝜔 2 )) (3-2-12)
Here, we find for the trilayer Terfenol-D/PZT-5A/Terfenol-D investigated in chapter 2 in LT-mode: 𝑅 𝑚 ′ = 278,28 Ω, 𝐶 𝑚 ′ = 35.22 pF and 𝐿 𝑚 ′ = 0.1345 H

In order to validate the procedure, a comparison in the dynamic regime between this extracted electric elements method and the FEM simulation is shown in Figure 3.7. We can confirm the good concordance between both methods. Thus, the extracted electric elements method can be explored to establish electrical equivalent circuit models that can be integrated in circuit simulators (e.g. VHDL-AMS, Verilog A, Spice, etc.).

Establishment of the optimal electrical load

The deliverable output power can be calculated with the output current 𝐼 𝑜𝑢𝑡 in using the equivalent Norton equivalent circuit represented in Figure 3.8: If the load is a complex impedance 𝑍 = 𝑅 𝐿 + 𝑗𝑋 𝐿 , the deliverable output power is given by :

𝐼 𝑜𝑢𝑡 = 𝐼 𝑖𝑛 𝑍 𝑇 𝑍 𝑇 +𝑍 (3-2-13)
𝑃 = 𝑅 𝐿 𝐼 𝑜𝑢𝑡 2 2 = | 𝑅 𝑇 +𝑗𝑋 𝑇 𝑅 𝑇 +𝑗𝑋 𝑇 +𝑅 𝐿 +𝑗𝑋 𝐿 | 2 𝑅 𝐿 2 𝐼 𝑖𝑛 2 = 𝑅 𝑇 2 +𝑋 𝑇 2 (𝑅 𝑇 +𝑅 𝐿 ) 2 +(𝑋 𝑇 +𝑋 𝐿 ) 2 𝑅 𝐿 2 𝐼 𝑖𝑛 2 (3-2-14)
It is well known that the deliverable power is maximal when the load is the complex conjugate of the source impedance, namely

𝑍 𝐿𝑜𝑎𝑑 = 𝑅 𝑇 -𝑗𝑋 𝑇 𝑃 = 𝑅 𝑇 2 +𝑋 𝑇 2 8𝑅 𝑇 𝐼 𝑖𝑛 2 (3-2-15)
If the load is resistive 𝑍 = 𝑅 𝐿 , the deliverable output power is given by : The optimal resistive load 𝑅 𝐿𝒐𝒂𝒅 is found in solving 𝜕𝑃 𝜕𝑅 𝐿 = 0 :

𝑃 = 𝑅
𝑹 𝑳𝒐𝒂𝒅 = √ 𝑹 𝑻 𝟐 + 𝑿 𝑻 𝟐 = |𝒁 𝑻 | (3-2-18)
Thus, in comparing the values of 𝒁 𝑻 at 75 kHz at 30 kOhm in Figure 3.6 with the optimal value of 18 kOhm found in Figure 3.4-(a), we can conclude that the magnetoelectric parameter 𝜸 is close to 1.5

[V•m/A].

Transient dynamic response with a non-linear electrical load

As mentioned previously, the harmonic regime works with the assumption of small signal and linear material. However, in the case of a nonlinear signal in harmonic regime due to, for instance, a nonlinear electrical load the transient analysis becomes necessary to predict the energy transfer. This is, for instance, the case that a synchronized switch damping technique (SSD) is employed [START_REF] Guyomar | Toward energy harvesting using active materials and conversion improvement by nonlinear processing[END_REF][START_REF] Li | Structural vibration damping with synchronized energy transfer between piezoelectric patches[END_REF][START_REF] Corcolle | Modeling of a beam structure with piezoelectric materials: introduction to SSD techniques[END_REF].

To correctly model the ME energy transducer with such an electronic circuit, the previously developed ME model equations needs to be solved in time domain. The use of a magneto-elastic-electric equivalent circuit model is possible but in contrast to the FEM, it does not take into account rigorously the mechanical and electrical impacts when the structure is loaded by electrical impedance.

In this work, the Newmark method is employed to calculate the transient dynamic response. For each time step, the system equation to be solved becomes:

[𝐾 ̅ ]{𝒳 𝑡+∆𝑡 } = {𝑅 𝑡+∆𝑡 } with [𝐾 ̅ ] = [ℳ] + ∆𝑡𝑎[𝒞] + ∆𝑡 2 𝑏 2 [𝒦] (3-2-19
) where

{𝑅 𝑡+∆𝑡 } = ∆𝑡 2 2 𝑏{ℱ 𝑡+∆𝑡 } + [ℳ] ({𝒳 𝑡 } + ∆𝑡{𝒳 𝑡 ̇} + ∆𝑡 2 2 (1 -𝑏){𝒳 𝑡 ̈}) + [𝒞] (∆𝑡𝑎{𝒳 𝑡 } + ∆𝑡 2 2 (2𝑎 - 𝑏){𝒳 𝑡 ̇} + ∆𝑡 3 2 (𝑎 -𝑏){𝒳 𝑖 ̈}) with {𝒳 ̈𝑡+∆𝑡 } = 2 𝑏∆𝑡 2 ({𝒳 𝑡+∆𝑡 } -{𝒳 ∆𝑡 }) - 2 𝑏∆𝑡 {𝒳 ∆𝑡 ̇} -( 1 𝑏 -1) {𝒳 ̈∆𝑡 } {𝒳 ̇𝑡+∆𝑡 } = {𝒳 ∆𝑡 ̇} + ∆𝑡 ((1 -𝑎){𝒳 ̈∆𝑡 } + 𝑎{𝒳 ̈𝑡+∆𝑡 }) , 𝑎 = 1, 𝑏 = 1 2 (𝑎 + 1 2 ) 2
The objective of all the SSD techniques is to maximize the energy of the system by modification of energy equilibrium of the system. For that, consider the following transient system of the energy transducer solved by the Newmark method:

[ 𝑀 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ] [ 𝑢V Q ä3 ̈] + [ 𝐶 𝑢𝑢 0 0 0 0 0 0 0 0 0 0 0 𝑅 0 0 0 ] [ 𝑢V Q ȧ3 ̇] + [ 𝐾 𝑢𝑢 𝐾 𝑢𝑝 0 𝐾 𝑢𝑎 𝐾 𝑝𝑢 𝐾 𝑝𝑝 𝐾 𝑝𝑞 0 0 𝐾 𝑎𝑢 𝐾 𝑞𝑝 0 0 0 0 𝐾 𝑎𝑎 ] [ 𝑢 𝑉 𝑄 𝑎 1 ] = [ 0 0 0 𝐽 𝑠 ] (3-2-20)
where 𝐽 𝑠 is the density of the source current. But as mentioned in chapter 2 (see 2.4.2.1.) the magnetic excitation can be implemented in considering the Dirichlet conditions on the magnetic vector potential 𝑎 1 in the boundaries of the problem domain Ω.

In considering the open circuit condition (R≈∞), the coupled system is solved with the following relations: The provided energy is:

∫ 𝐹 𝑒𝑓𝑓 𝑢ṫ 𝑜 𝑑𝑡 = 𝐹 𝑜 𝑢(𝑡) (3-2-24)
where 𝐹 𝑜 is a primitive in time of 𝐹 𝑒𝑓𝑓 .

The provided energy 𝐹 𝑜 𝑢(𝑡) is distributed into the kinetic energy This latter corresponds to the mechanical energy, which is converted into electrical energy composed by the electrostatic energy During a periodic oscillation, the kinetic energy, potential elastic energy and the electrostatic energy are null. Thus, during a periodic oscillation, the energy equilibrium can be simplified with the mechanical damping and the electrical load:

∫ 𝐶 𝑢𝑢 𝑢̇2 𝑡 𝑜 𝑑𝑡 + ∫ 𝑉𝐼 𝑡 𝑜 𝑑𝑡 = 𝐹 𝑜 𝑢(𝑡) (3-2-26)
The maximization of this energy leads to minimization of the mechanical energy in the structure. Thus, the usefulness of a SSD control approaches is to maximize this energy in putting the voltage and displacement velocity in phase.

Several techniques exist among them one imposes an inductance, which leads to the so-called Synchronized Switch Damping on Inductor (SSDI) technique [START_REF] Guyomar | Toward energy harvesting using active materials and conversion improvement by nonlinear processing[END_REF][START_REF] Li | Structural vibration damping with synchronized energy transfer between piezoelectric patches[END_REF][START_REF] Corcolle | Modeling of a beam structure with piezoelectric materials: introduction to SSD techniques[END_REF] as presented by the both possible cases in Figure 3.9. The case 1 uses a SSDI technique with L and R in parallel whereas the case 2 uses a SSDI technique with L and R in series. For both cases, the switch is open, except when zero displacement velocities 𝑢̇𝑡 +∆𝑡 occur. In this case, the switch is locked in a very short period that puts in phase the output voltage and displacement velocity. In this short period 𝑡 𝑖 = 𝜋√𝐿𝐶 𝑜 the clamped capacitor 𝐶 𝑜 and the added inductance L represents an oscillator. , we obtain :

-𝐿𝑄 ̈+ 𝐿 𝑅 𝑉 ̇𝐾𝑝𝑞 + 𝑉𝐾 𝑞𝑝 = 0 (3-2-27)
The system equation (3-1-20) to be solved becomes:

[ 𝑀 0 0 0 0 0 0 0 0 0 0 0 -𝐿 0 0 0 ] [ 𝑢V Q ä3 ̈] + [ 𝐶 𝑢𝑢 0 0 0 0 0 0 0 0 0 𝐿/𝑅 0 0 0 0 0 ] [ 𝑢V Q ȧ3 ̇] + [ 𝐾 𝑢𝑢 𝐾 𝑢𝑝 0 𝐾 𝑢𝑎 𝐾 𝑝𝑢 𝐾 𝑝𝑝 𝐾 𝑝𝑞 0 0 𝐾 𝑎𝑢 𝐾 𝑞𝑝 0 0 0 0 𝐾 𝑎𝑎 ] [ 𝑢 𝑉 𝑄 𝑎 3 ] = [ 0 0 0 𝐼 𝑛 ]
In case 2, the differential circuit equation is introduced in the system through considering the voltage equation 𝑢(𝑡) = 𝑢 𝑅 (𝑡) + 𝑢 𝐿 (𝑡), with 𝑢 𝑅 = 𝑅𝑖 and 𝑢 𝐿 = 𝐿𝜕 𝑡 𝑖 𝐿 . Thus, in considering 𝜕 𝑡 𝑢 = 𝑅𝜕 𝑡 𝑖 + 𝐿𝜕 𝑡 (𝜕 𝑡 𝑖), we obtain :

-𝐿𝑄 ̈-𝑅𝑄 ̇+ 𝑉 ̇𝐾𝑞𝑝 = 0 (3-2-28)
The system equation (3-2-20) to be solved becomes:

[ 𝑀 0 0 0 0 0 0 0 0 0 0 0 -𝐿 0 0 0 ] [ uV Q ä3 ̈] + [ 𝐶 𝑢𝑢 0 0 0 0 0 0 0 0 0 𝐾 𝑞𝑝 0 -𝑅 0 0 0 ] [ 𝑢V Q ȧ3 ̇] + [ 𝐾 𝑢𝑢 𝐾 𝑢𝑝 0 𝐾 𝑢𝑎 𝐾 𝑝𝑢 𝐾 𝑝𝑝 𝐾 𝑝𝑞 0 0 𝐾 𝑎𝑢 0 0 0 0 0 𝐾 𝑎𝑎 ] [ 𝑢 𝑉 𝑄 𝑎 3 ] = [ 0 0 0 𝐼 𝑛 ] (3-2-29)
In both cases, as shown in Figure 3. 10 [52], when the SSDI harvesting technique is activated (i.e. time > 0.13ms) the output voltage reaches a maximal value and decreases until a steady state and the deliverable energy stored in the load ∫ 𝑉𝐼𝑑𝑡 is strongly gradually increased.

Multilayer ME composite materials

Although the above results are encouraging, the bonding layer between adjacent layers inherently limits the performances of laminate materials. This reduces the quality of the mechanical contact between the layers and thus the ME coupling, and the lifetime of the device due to fatigue, as well. This is particularly true for operation of bulk sensors at their resonance frequencies of tens of kHz. Such high resonant frequencies can also result in significant eddy current losses due to the conductive nature of the magnetostrictive alloys and further lower the efficiency of energy conversion. Furthermore, currently used magnetostrictive material, such as Terfenol-D, requires large magnetic bias fields 100-300 Oe to achieve highest ME coefficients.

To address this issue, recent studies using deposition process technologies (such as magnetron sputtering or vaporizer) have focused on emergence of a new generation of ME composites constituted of flexible thin films such as NFZO/PZT or CFO/PZT [START_REF] Kambale | Current Status of Magnetoelectric Composite Thin/Thick Films[END_REF][START_REF] Ortega1 | Multifunctional magnetoelectric materials for device applications[END_REF]. However, the epitaxial deposition processes are limited and the thickness of each thin film is usually around a few µm [START_REF] Kambale | Current Status of Magnetoelectric Composite Thin/Thick Films[END_REF][START_REF] Ortega1 | Multifunctional magnetoelectric materials for device applications[END_REF], that involves a low energy density by piezoelectric layer. It is also well known that performance (greater current or voltage output) of a structure composed of multilayer piezoelectric can be adjusted with either a series or parallel electrical connection. In all cases, the advantage of a multilayer piezoelectric structure is that it produces a higher power than a single-layer generator having the same total thickness [START_REF] Nath | Enhanced piezoelectric response from barium strontium titanate multilayer films[END_REF][START_REF] Sun | Modeling for Piezoelectric Stacks in Series and Parallel[END_REF][START_REF] Zhu | A Bimorph Multilayer Piezoelectric Vibration Energy Harvester[END_REF]. The interest of multilayer MEC is apparent since in according to the electrical connection type (series or parallel) it is possible to obtain a high output voltage to operate such as an efficient energy transduced or to provide a high output current to operate such as efficient magnetic sensor. For instance, in [START_REF] Yuri | Ferrite-Piezoelectric Multilayers for Magnetic Field Sensors[END_REF] the MEC multilayer sample (operating in LT-mode) is 3x6 mm 2 with respectively 11 NZFO layers and 10 PZT layers (with each layer being of 18 µm), which provides a magnetoelectric coefficient close to 1 mV/Oe under resonance. In this context, establishment of models are required to compare the performances between a trilayer and a multilayer MEC with the same material composition and an unchanged dimension. Specific attention should be focused on the effects when an alternative opposing polarization is applied between the superposed piezoelectric layers in according to the operating magnetization-polarization mode.

Multilayer ME composite problems description and FEM modelling

The modeling of the multilayer ME composites can be described in terms of the relationship between the ME coefficient and the physical parameters of a multilayer laminate composite structure, consisting of alternatively bonded magnetostrictive (M) layers and piezoelectric (P) ones under an external magnetic field 𝐻 𝑒𝑥𝑡 . The Dirichlet boundary conditions of all mechanical, electrical and magnetic will be the same as those for the trilayer MPM composite problem, which have been presented in Chapter 2.

As mentioned previously, there are four possible ME modes, which are combinations of two types of magnetostriction modes in M-layer and two polarization modes, denoted as L-L, L-T, T-L and L-L mode, respectively. And for each ME modes, there are two electrical connectivity types either in parallel or in series. The typical modes are displayed in Figure 3.11 (the electrodes and wires are marked with the light yellow color): L-T mode: L-L mode:

Equivalent circuit method for multilayer ME composite problems

The equivalent circuit models (ECM) for MEC in LL and LT mode are recapitulated in Figure 3.13 and the expressions of their magnetoelectric coefficients are given by (3-3-1) and (3-3-2). Details of the derivation are given in the Appendix D. It can be noticed that equivalent circuit models are the same for both the trilayer and multiplayer situations but the circuit parameters depends on the number of layers.

Consider the subdivision of a trilayer MEC in 𝑚 magnetostrictive layers and (𝑚-1) piezoelectric layers. The total thickness is maintained with this condition:

𝑡 = 2𝑡 𝑚 + 𝑡 𝑝 = 𝑚𝑡 𝑚 ′ + (𝑚 -1)𝑡 𝑝 ′ (3-3-3)
where 𝑡 𝑝 , 𝑡 𝑚 , 𝑡 𝑝 ′, 𝑡 𝑚 ′ are respectively the thicknesses of the piezoelectric and magnetostrictive layers for the trilayer and multilayer, both illustrated in Figure 3.14, where 𝑤 is the normalized width.

In considering the ratio 𝛿 = 𝑡 𝑝 𝑡 𝑝 ′ between the thicknesses of each single piezoelectric layer in the trilayer and the multilayer cases, the total equivalent clamped capacitor 𝐶 𝑇 for the series and parallel electrical connections are given by the following expressions in Table 3 In considering the ratio 𝛾 = 𝑡 𝑚 𝑡 𝑚 ′ , (that is different to 𝛿 = 𝑡 𝑝 𝑡 𝑝 ′ ) the magnetomechanical and electromechanical coupling coefficients 𝜑 𝑚 , 𝜑 𝑝 and the mechanical impedance 𝑍 𝑚 are also function of the change of layers in accordance with the configuration applied. Table 3 Thus, the magnetoelectric coefficients (3-3-1) and (3-3-2) are rewritten for a multilayer case as:

𝛼 𝑣 ̃= | 𝑈 𝐻 𝑎𝑐 | mode 𝐿𝐿 = | 𝜑 𝑝′ 𝜑 𝑚′ 𝑗𝑍 𝑚 𝐶 𝑇 𝜔 | (3-3-4) 𝛼 𝑣 ̃= | 𝑈 𝐻 𝑎𝑐 | mode 𝐿𝑇 = | 𝜑 𝑝 ′𝜑 𝑚 ′ 𝜑 𝑝 ′ 2 +𝑗𝑍 𝑚 𝐶 𝑇 𝜔 | (3-3-5)

Results and comparisons

The studied example is a MEC composed of a Terfenol-D and PZT-5A laminated composites. The device has the following dimensions: l=14mm, t=3mm (t m =t p =1mm for the trilayer case). The material property parameters are given in the Appendix A.

The LT mode as shown in Figure 3.11 (a) is firstly considered. The FEM and the ECM simulation results of the magnetoelectric coupling coefficient for parallel and series connections at the resonance frequency are compared in Figure 3.15, when the number of layers increases. A good concordance between the two methods is observed. In addition, it can be noticed that with the series connection, the coupling coefficient increases slightly with the increase of the number of layers. This can be explained by the diminution of the clamped capacitor 𝐶 𝑇 according to the ECM formula in Table 3.3 which results in an increase of 𝛼 𝑣 ̃ according to equation (3-3-5). On the other hand, in parallel connection, 𝐶 𝑇 increases with the number of layers, the effect is inversed and the coupling coefficient 𝛼 𝑣 ̃ diminishes quickly, according to equation (3-3-5). The voltage provided by the parallel connection diminishes but instead, the transducer allows furnishing the greater current.

Series electrical connection

Parallel electrical connection For the LL mode, two configurations are considered. In the first configuration (configuration 1), all the piezoelectric layers are polarized in the same direction, as shown in Fig. 3.11. In the second configuration (configuration 2), the polarization of the piezoelectric layers is altered for the neighboring layers, while the magnetization of magnetostrictive layers remain the same, as illustrated lately in Figure 3.19 and Figure 3.20. Figure 3.17 shows the comparison between the FEM and the ECM simulation results of the series and parallel connections. Contrasting in the LT mode, we can notice a good concordance just with the configuration 1 in the parallel electrical connection. In the case of series connection for configuration 1, as previously stated, there is a dissymmetry of the field distribution due to the connection of the positive and negative electrodes situated on the opposite sides. This dissymmetry cannot be reflected by the ECM because of the one dimension assumption in the analytical formula. The electric potential distributions for a trilayer case and for the configuration 1 are given in Figure 3.18. The dissymmetry of the series connection is clearly illustrated.

Series electrical connection Parallel electrical connection In the case of the configuration 1, the series connection presents also a technical realization difficulty.

For that reason, the configuration 2 by alternating piezoelectric polarization makes more sense. However, as shown by the electric potential distribution given in Figure 3.19 and Figure 3.20, the field results depends on the number of layers. The symmetry of the field distribution is exhibited only for odd number of piezoelectric layers while for even number of piezoelectric layers the field distribution is dissymmetrical. This explains why the computed magnetoelectric coupling coefficient by the FEM zigzags with the number of layers for both series and parallel connections. It can be noticed that, as the ECM does not take into account the alternating polarization of piezoelectric layers, it is not valid for the configuration 2.

Conclusion

In conclusion, the finite element method has been applied for both harmonic and transient cases and successfully applied to an energy transducer composed of Terfenol-D/PZT/Terfenol-D materials. The harmonic model allows to determine the resonance frequencies and the optimal loads to maximize the deliverable power. The transient model using the Newmark method provides a useful tool to study the energy transfer when a conditioning circuit such as SSDI harvesting technique is employed. An equivalent circuit method has been developed to analyze multilayer MEC of two typical magnetization-polarization modes (LL and LT) with different electrical connections (series and parallel). The results demonstrate two important conclusions concerning the change of the ME coefficient in function of the number of the deposited layers in the composites: first, the change tendencies between the results utilizing FEM and the ECM remain quite similar only for some simple configurations with symmetrical field distributions, which denotes the validation of ECM for these cases; second, the change tendency (increase or decrease with increasing layers number) of the ME coefficient depends highly on the electrodes connection type, which can be analytically explained that the variation of the clamped capacitance and the internal impedance of the piezoelectric layers play an important role in these results.

Chapter 4. Prospective application of ME composites as energy transducer 4.1 Introduction

This section proposes a prospective application of the ME composites as energy transducers. For that, a potential application in biomedical domain is proposed in first part. To emphasize this proposition some measurements on a bilayer ME composite are presented in the second part.

Potential application in biomedical domain

Remote monitoring of medical physical data, such as the pressure or temperature of a patient suffering from a disease, would constitute a major advance in the medical domain, since it would improve the quality of the life and the medical service. These considerations motivate researchers to develop novel autonomous wireless sensors while overcoming two major drawbacks as for the powering of the sensors and the IC circuits. The first is the use of batteries, which are pollutant, and with limit lifetime. The wireless powering or energy harvesting solutions must replace them. The second concerns the insertion in-body that must be done by a minimally invasive technique, for example the incorporation of a catheter through the femoral artery for monitoring the temperature or blood pressure. The common solutions for wireless energy harvesting and data communication for human implants are the piezoelectric transducers [START_REF] Weiss | Active Medical Implant[END_REF], the electromagnetic RF transducers with antennas operating at 2.45 GHz [START_REF] Mizuno | Development of an implanted helical folded dipole antenna for 2.45 GHz applications[END_REF] or 900 MHz [START_REF] Basari | A simple folded dipole antenna for medical implant communications at 900 MHz band[END_REF], and the near field communication (NFC) in using inductive coupling circuits on small coils operating at resonance in the radio-frequency (RF) of 1-300 MHz range [START_REF] Guo | Inductive wireless power transmission for implantable devices[END_REF]. The piezoelectric transducers are not suitable in-body insertion case because the low amplitude and the low frequency of the vibrations cannot produce sufficient power. The electromagnetic RF transducers are not suitable since the absorption of the electromagnetic radiation by the human body is significant above 300 MHz. The NFC technologies such as the Radio Frequency Identification (RFID) passive devices (batteryless) are suitable if the working frequency is included in the Industrial, Scientific and Medical (ISM) band. Although the use of passive RFID armbands at 13.56 MHz (within the ISM band) to collect data from an identified patient is common, their extension as a miniaturized energy transducer to supply an implanted biomedical sensor (such as telemetric pressure sensors) in human body seems limited by the coupling factor because its communication distance is usually limited to 4 cm [START_REF] Hierold | Implantable low power integrated pressure sensor system for minimal invasive telemetric patient monitoring[END_REF][START_REF] Ginggen | A Telemetric Pressure Sensor System for Biomedical Applications[END_REF][START_REF] Schnakenberg | Initial investigations on systems for measuring intraocular pressure[END_REF][START_REF] Mokwa | Medical implants based on microsystems[END_REF]. Thus, in order to power an implanted sensor, the ME trilayer composites seem an alternative energy collector other than the usual RFID to reach a vicinity contactless distance communication greater that 4 cm.

In considering that the IC chip works with the power requirements previously exposed, namely at least 0.8 mW, the advantages of use the trilayer Terfenol-D/PZT-5A/Terfenol-D having the dimensions L x =14mm, L y =3mm (t m =t p =1mm), L z =3mm are :

-Incorporation into a medical catheter for the minimally invasive purpose could be guaranteed as illustrated in Figure 4.1.

-Low absorption by the human body by the low frequencies at frequency f<100 kHz delivered by an external source (reader excitation) distance.

Ferrite solenoid as external reader

Considering that Helmholtz coils are placed in the human thorax to produce the biasing static field.

The external reader will be used just to produce the AC field H ac . The ideal solution would be a compact antenna. A recent study [START_REF] O'handley | Improved Wireless, Transcutaneous Power Transmission for In Vivo Applications[END_REF] has proposed to use the solution shown in Figure 4.2. This solution is based on a long ferrite solenoid surrounded by a simple coil. According to authors, it is possible to obtain at least 2 Oe at 10 cm. This value would ensure a quasi-continuous voltage required to power the IC for a distance greater than 4 cm (limitation of the RFID technique).

As shown in Figure 2.33 in chapter 2 for a trilayer Terfenol-D/PZT-5A/Terfenol-D at resonance the ME coefficient is close 10 V/Oe. Thus, in considering a field level in the order of δ𝐻 𝑎𝑐 =2 Oe the induced voltage between two electrodes of the piezoelectric reaches 20V. In this case, the output deliverable power at the resonance frequency is respectively close to 3.2 mW. This result put in evidence that the ME transducer in L-T mode could deliver the power necessary for a biomedical IC chip with a distance communication greater than 10 cm.

Exposition limits of the magnetic field

The exposition limits of the magnetic field (static and dynamic regimes) must respect the guidelines established by the International Commission on Non-Ionizing Radiation Protection (ICNIRP) [START_REF]International Commission on Non-Ionizing Radiation Protection (ICNIRP)[END_REF].

Concerning the static field, the ICNIRP recommends in the Table 4.1 that static fields on humans, acute exposure of the general public should not exceed 400 mT but for specific work applications such as a biomedical context that "the limit on exposure is set at 2 T, to prevent vertigo, nausea and other sensations, but for specific work applications, when the environment is controlled and appropriate work practices are implemented, then exposure up to 8 T is acceptable".

Table 4.1 Limits of exposure to static magnetic fields [START_REF]International Commission on Non-Ionizing Radiation Protection (ICNIRP)[END_REF] Concerning the dynamic field working at frequency 75 kHz, the reference level shown in Table 4.2 is 80 A/m i.e. 1 Oe for occupational exposure such as a biomedical context [START_REF]International Commission on Non-Ionizing Radiation Protection (ICNIRP)[END_REF].

Table 4.2 Reference levels for occupational exposure to time-varying electric and magnetic fields [START_REF]International Commission on Non-Ionizing Radiation Protection (ICNIRP)[END_REF] Exceeding the reference value of 1 Oe is not proscribed, but in this case the limits of exposure (base restrictions) must be positioned in according to value of the current density for all tissue of head and body (must not exceed 150 mA.m -2 at 75 kHz for instance).

Measurement of a bilayer ME composite

Measurement bench set-up

In order to estimate the real performance of a material composite ME as a potential energy transducer, this section shows the experimental data from the bilayer ME composite composed of Terfenol-D/PZT-5H represented and described in As shown in Figure 4.7, a magnetic solenoid coil composed of 100 turns is used to generate the externally small signal harmonic ac field 𝐻 𝑎𝑐 . The magnetic field in the middle of the solenoid can be approximated with the following expression:

𝐻 = 𝑁𝐼 2𝑎 𝑎 √𝑎 2 +𝑅 2 (4-1)
where 2𝑎 and 𝑅 represent respectively the length and the radius of the solenoid coil. 𝐼 is the crossing current along the coil that can be measured in placing a multimeter in series.

Static and dynamic measure responses

The measurements are achieved with the experimental set-up shown in Figure 4.8 that was realized by Kevin Malleron, a PhD student at the L2E Lab. The ME composite is placed inside the solenoid coil that is supplied by a function generator (HAMEG HM8 150 from Rhode&Schwarz) with 20 V. In low frequency (10 Hz 1 kHz), the magnetic field inside the solenoid coil can be estimated is connecting in series a multimeter (FLUKE 179) to measure the crossing current 𝐼 at each work frequency. Nevertheless for frequencies greater than 1 kHz must be variable resistor is placed in series with the function generator to stabilize the value of the crossing current 𝐼 around 10 mA to impose an AC magnetic field of 1 Oe inside the solenoid. A numerical oscilloscope (Agilent Tech MSO7054A) is used to display and to measure the output voltage from the ME composite (between both electrodes of the PZT-5H).

Deliverable output power

To estimate the deliverable output power, a variable resistor R has been placed between both electrodes to establish the value of the internal impedance 𝑍 𝑇 of the ME composite. For that, as illustrated in Figure 4.11 the principle of the tension divider bridge is used namely that 𝑅 = 𝑍 𝑇 when the output voltage is divided by two. where or 𝑉 𝑖𝑛 represents the stress force due to the Terfenol-D as a result of the externally applied magnetic field 𝐻 𝑎𝑐 Due to the mechanical losses, the measured value is close to 300 Ohm. Thus, the estimation of the deliverable output power for this ME composite with an excitation AC dynamic magnetic field of 1 Oe is:

𝑃 𝑜𝑢𝑡 = 𝑉 𝑜𝑢𝑡 2 2𝑍 𝑇 ≈ 0.36 600 = 600 µ𝑊
This power estimation is close to 800 µW founded with the ME composite Terfenol-D/PZT-5A/ Terfenol-D in the first section of the chapter 3. For an excitation AC dynamic magnetic field of 2 Oe the power will be close to 1 mW that is sufficient to supply a WSN sensor.

Conclusion

This section has proposed a prospective application of ME composites Terfenol-D/PZT-5A/ Terfenol-D as energy transducer in a biomedical domain that would respect the exposition limits established by the International Commission on Non-Ionizing Radiation Protection. Measurements on a bilayer ME composite Terfenol-D/PZT-5H results have highlighted this possible prospective application in exhibiting an deliverable output power close to 600 µW under an excitation AC dynamic magnetic field of 1 Oe.

General conclusion

The work presented in this thesis is a contribution to modelling of ME composites for energy harvesting. Two parts have been mainly studied in this work: the analytical and numerical modelling for which we have carefully analyzed the static and dynamic behaviors of different ME structures, and a practical application oriented to highlight the possible assessment of ME structures as energy transducers.

Thus, the first chapter has familiarized us to ME structures by introducing the magnetoelectric effect origin according to the phases of a composite composed of laminated piezoelectric and magnetostrictive layers. The presentation of several practical applications showed how the ME structures have become today major assets in the improvement and design of new electronic devices such as the energy transducers.

The second chapter focused on a description of the various simulation tools and used in this dissertation to describe the behavior of ME structures composed of laminated piezoelectric and magnetostrictive layers. Firstly, two analytical methods in 0D and 1D assumptions have been addressed and detailed. Secondly, a description of a multiphysics code based on the finite element method in 2D using a discretization of quantities of electromagnetic and elastic has been performed. It has been shown that the involved fields, as well as coupling effects between these fields, can be conveniently modelled in this FEM framework. Particularizing electrical circuit taking into account the electrical impact when the structure is loaded by an electrical impedance has been added into the FEM formulations. For constitutive modelling, the linear piezoelectricity for piezoelectric materials and the nonlinear magnetization and magnetostriction for magnetostrictive materials have been considered. To deal with nonlinearity a piecewise linear strategy was adopted in the FEM code.

In the third chapter we highlighted in the first part that the LT-mode is the ideal configuration to obtain the sufficient deliverable output power to supply a wireless sensor node in working with 1% duty-cycle. In the second part, a transient model using Newmark method has been employed to study the energy transfer when a conditioning circuit such as SSDI harvesting technique is used. In the last part, an equivalent circuit method to analyze multilayer MEC in LL and LT modes with series and parallel connections has been detailed and investigated. The results showed that the LT mode in series connection could be used to conserve the performances of ME composites constituted of thin films.

The fourth chapter has proposed a potential application of a ME energy transducer to supply an implanted biomedical sensor with a specific external coil while respecting the magnetic field limits exposition imposed by the International Commission on Non-Ionizing Radiation Protection. A set-up experimentation to measure the ME coefficients in static and dynamic regimes has been presented and measurements on a bilayer composite have been performed. The estimated output deliverable power close to 600 µW has confirmed the potential possibility to supply a wireless sensor node in working with 1% duty-cycle.

Finally, the work of this thesis is a first step to open the way to contribute to bridge the gap between the physical concepts and the design optimization of devices by providing tools for modelling and designing of devices such as ME energy transducers according to imposed specifications (desired application, choice of materials, dimensions imposed, etc.). The health domain of remote power supply is a perfect example, since the proposed energy transducer can meet the constraints on dimensions (the micrometer and millimeter scales presage no-invasive implantation) by adjusting the magnetic field amplitude at a fixed distance of the reader or by adjusting its distance for fixed magnetic field amplitude. The environment effects (temperature, biocompatibility, pressure...etc.) must be addressed in future works in the FEM numerical code. where 𝑠 𝑖𝑗𝑘𝑙 is the elastic compliance tensor, 𝑚 𝑖𝑗𝑘𝑙 the field magnetostrictive modulus tensor, 𝑟 𝑖𝑗𝑘𝑙𝑚𝑛 the field magnetoelastic modulus tensor, 𝜇 𝑘𝑙 the permeability tensor.

𝐺 = 𝐺

In one-dimensional problem, the constitutive relations can be simplified as

𝜀 = 𝑠𝜎 + 𝑚𝐻 2 + 𝑟𝜎𝐻 2 (B-7) 𝐵 = 𝜇𝐻 + 𝑚𝜎𝐻 + 𝑟𝜎 2 𝐻 (B-8)
and the magnetostrictive modulus 𝑚 and magnetoelastic modulus 𝑟 can be calculated as

𝑚 = 𝑑 ̃0 2𝐻 ̃0 (B-9) 𝑟 = 1 𝜎 ( 𝑑 2 𝐻 ̃-𝑑 ̃0 2𝐻 ̃0) (B- 10
)
where 𝑑 = (𝜕𝜀 𝜕𝐻 ⁄ )| 𝜎 indicates the piezomagnetic coefficient, and the waves above d and H indicate the peak values of these two parameters and the subscript "0" indicate the values under no external pre-stress, which can all be obtained based on relevant experiment results.

In three-dimensional problem, for an isotropic material as example, the modulus must be calculated respectively as (where 𝛿 𝑖𝑗 is the Kronecker delta) by measuring 𝑚 1111 and 𝑚 1122 denoting respectively the strains caused by magnetic field of per unit in the directions along and perpendicular to the external magnetic field, and also by characterizing the dimensionless coefficient C in terms of the magnetoelastic modulus in the one-dimensional cases.

By comparing the calculation results with those from experiments, it can be found from Fig. B.1 [START_REF] Carman | Nonlinear constitutive relations for magnetostrictive materials with applications to 1-D problems[END_REF] that the SS model presents good agreements only under the region of low and moderate magnetic fields and high pre-stress. The obvious deficiencies are that it cannot describe the saturation magnetostriction under intensive magnetic fields and it becomes unavailable when applied small prestresses. where 𝑘 = 1/𝐻 ̃ is a relaxion parameter to make independent variable of the hyperbolic function dimensionless, and 𝐻 ̃ the external magnetic field corresponding to the peak piezomagnetic coefficient defined in the previous section.

Then, by utilising the Legendre transform relations showed in the previous section, the constitutive laws can be obtained for one-dimensional problems as [92], one can see that the magnetostriction saturations in the magnetostrictive material in the region of intense magnetic fields can be simulated by this model, but there are large discrepancies between the experimental and predicted results (reach to overestimation of about 40% when under large fields and compressive pre-stress of 65.4 MPa). Then, based on the decomposition of the total strain as the sum of the elastic strain 𝜀 𝑒 and the mechanical-magnetic coupling strain 𝜀 𝐻 (which can be given by integration of the expression of d) and the magnetic flux density as the sum of the magnetization part and the mechanical-magnetic coupling part, the expression of the Gibbs free energy can be derived by the similar method of Hom et al [START_REF] Hom | A fully coupled constitutive model for electrostrictive ceramic materials[END_REF] for electrostrictive materials as (where 𝐺 0 (𝐻) = ∫ 𝜇𝐻𝑑𝐻 𝐻 0

𝜀 = 𝑠𝜎 + 1 𝑘 2 𝑚𝑡𝑎𝑛ℎ 2 (𝑘𝐻) + 1 𝑘 2 𝑟𝜎𝑡𝑎𝑛ℎ 2 (𝑘𝐻) (B-14) 𝐵 = 𝜇𝐻 +
): 𝐺 = 𝐺 0 (𝐻) + [92], one can see that this model can simulate the experimental curves more effectively than the two above models and capture the characteristics of the mechanism of the magnetiomechanical coupling in a one-dimension case problem for the magnetostrictive material, such as the saturation magnetostriction under high magnetic fields and the notable magnetoelastic dependence on the external stress. But the predicted results are about 30% underestimated in the region of high magnetic fields when the compressive pre-stress is 65.4 MPa.

FigureB.3 Magnetostriction curves by experiments and the DDS model

The Zheng-Liu model

This model, proposed by Xiaojing Zheng and Xin'en Liu based on a more adequate capture of the nonlinear main magnetoelastic coupling characteristics in the magnetostrictive materials [START_REF] Liu | A nonlinear constitutive model for magnetostrictive materials[END_REF][START_REF] Liu | A nonlinear constitutive model for Terfenol-D rods[END_REF], can overcome the deficiencies presented by the models above. The model is built by taking into account some important modifications and employments, such as reserving more high-order terms in the Taylor series expansion of the Gibbs free energy function, relating the pre-stress induced nonlinear strain to the magnetic domain rotations and describing the magnetization curve of no external stress by using the Langevin function. To be detailed based on the total differentiation of the Gibbs free energy function with stress 𝜎 𝑖𝑗 and magnetization 𝑀 𝑘 as the independent variables without considering the heat effect term:

𝑑𝐺 = -𝜀 𝑖𝑗 𝑑𝜎 𝑖𝑗 + 𝜇 0 𝐻 𝑘 𝑑𝑀 𝑘 (B-22)
and based on the decomposition of the total strain 𝜀 𝑖𝑗 as the sum of the elastic strain produced by a pre-stress 𝜀 𝑖𝑗 0 (𝜎 𝑘𝑙 ), which can be obtained by mechanical constitutive equations, and the magnetoelastic coupling strain 𝜀 𝑖𝑗 1 (𝑀 𝑘 , 𝜎 𝑚𝑛 ), which can be obtained from the coupling terms of the function of G, and also the analogous decomposition of the total magnetic field 𝐻 𝑘 = 𝐻 𝑘 0 (𝑀 𝑙 ) + 𝐻 𝑘 1 (𝑀 𝑙 , 𝜎 𝑚𝑛 ), the constitutive relations can be derived as: In which, the fourth-order isotropic tensors and the inverse function for 3-D case can be expressed as: From the comparison between the results from the numerical model and experiments for a 1-D Terfenol-D in Fig. B.4 [START_REF] Liu | A nonlinear constitutive model for magnetostrictive materials[END_REF], one can see that the results have a great coincidence between each other under various applied compressive pre-stress in both low and high magnetic field regions, which demonstrates a wider applicability and higher accuracy than the models above, and also a more effective characterization for the effect of the axial pre-stress or in-plane residual stress on the material characteristic curves in the further studies of the literature.

𝜀 𝑖𝑗 =
𝑆 𝑖𝑗𝑘𝑙 =

B.2 Hysteresis models

The Preisach model is one of the earliest tools used to model the hysteresis of ferromagnetic materials by employing a hysteresis operator to take the value ±1 depending on the current state and history based on the previous state, and improved by taking into account the stress effect recently to model the magnetoelastic behavior of magnetostrictive materials. However, they are purely mathematical and do not address the underlying physics of the problem [START_REF] Suzuki | Magnetoelastic behaviour of ferromagnetic materials using stress dependent Preisach model based on continuum theory[END_REF].

The original Jile-Atherton model has also been generalized to a vector model [START_REF] Bergqvist | A simple vector generalization of the Jiles-Atherton model of hysteresis[END_REF], and extended to incorporate the effect of stress by Li and Jiles [98], the magnetocrystalline anisotropy [START_REF] Raghunathan | Generalized form of anhysteretic magnetization function for Jiles-Atherton theory of hysteresis[END_REF] and the temperature effects [START_REF] Raghunathan | Theoretical model of temperature dependence of hysteresis based on mean field theory[END_REF].

There are also many hysteresis models such as those proposed by Smith et al [START_REF] Smith | A unified framework for modelling hysteresis in ferroic materials[END_REF], Falk [START_REF] Falk | Landau theory and martensitic phase transitions[END_REF] and Armstrong [START_REF] Armstrong | An incremental theory of magnetoelastic hysteresis in pseudo-cubic ferromagnetostrictive alloys[END_REF]. They all have their particular advantages and limitations respectively. (C-4)

𝑎
𝜕𝐹(𝑋) 𝜕𝑋 = [ 𝜕𝐹 1 𝜕𝑥 1 𝜕𝐹 1 𝜕𝑥 2 𝜕𝐹 2 𝜕𝑥 1 𝜕𝐹 2 𝜕𝑥 2 ⋯ 𝜕𝐹 1 𝜕𝑥 𝑛 ⋯ 𝜕𝐹 2 𝜕𝑥 𝑛 ⋯ ⋯ 𝜕𝐹 𝑛 𝜕𝑥 1 𝜕𝐹 𝑛 𝜕𝑥 2 ⋯ ⋯ ⋯ 𝜕𝐹 𝑛 𝜕𝑥 𝑛 ] (C-5)
We assume that 𝑋 𝑘 indicates the value of the unknown vector X in the step k of the iteration process. The detailed calculation procedure is as follows:

Step 1: Set the initial value of the unknown variable vector to be zero as 𝑋 0 = 0

Step 2: According to the Newton's equation [START_REF] Sun | A Practical Numerical Analysis[END_REF] Step 4: Take 𝑋 = 𝑋 𝑘+1 as the solution of the problem 𝐴𝑋 -𝑏 = 0 and terminate the iteration.

The process of iteration can be understood in a visual way. Since it is displayed that in every step of the iteration we calculate a new value of 𝑩 by utilizing the value of the last step, which can be reasonably considered as a new value is obtained through a function F(𝑩), therefore the entire iteration procedure can be equivalently treated as solving an equation in form of F(𝑩)= 𝑩 (in which it's not necessary to know the exact expression of F(𝑩)). When we start from that 𝑩 is zero, the iteration process can be visually presented as: directly. Hence, we develop an approach to obtain an approximate calculation result of the matrix.

For each element in the tangent stiffness matrix expression (C-5) with 1 ≤ 𝑖, 𝑗 ≤ 𝑛 , we can derive that Thus, we obtain the values of the matrix A, by using which we can subsequently obtain the values of the partial derivative terms in (C-7) and finally calculate the tangent stiffness matrix.

Appendix D. Multilayer analytical modelling

For an ME multilayer composite in 2D case, by applying the plane stress condition and the 1-D assumptions as introduced in the Chapter 2.

In LT mode:

Assuming to apply two surface forces along the longitudinal direction on both ends of the composite, we can obtain their expressions by applying the constitutive laws (2-3-10) for LT mode in 2D case as:

𝐹(0) = -𝑚𝑡 𝑚 ′𝑤𝑇
In the series connectivity type, the output voltage 𝑈 = (𝑚 -1)𝑡 𝑝 ′𝐸 2 , thus 𝜑 𝑝 ′ = 𝑤𝑑 31 𝑝 𝑠 11 𝐸 , and the coupling current through the piezoelectric layers can be calculated as [START_REF] Dong | Longitudinal and Transverse Magnetoelectric Voltage Coefficients of Magnetostrictive/Piezoelectric Laminate Composite: Theory[END_REF]: In LL mode:

𝐼 𝑝 =
In the series electrical connectivity type, the output voltage can be expressed as: Thus assuming to apply two surface forces along the longitudinal direction on both ends of the composite, we can obtain their expressions by applying the constitutive laws for LL mode in 2D case as: And the coupling current induced through the piezoelectric layers can be calculated as [START_REF] Dong | Longitudinal and Transverse Magnetoelectric Voltage Coefficients of Magnetostrictive/Piezoelectric Laminate Composite: Theory[END_REF]: In the parallel electrical connectivity type, the output voltage can be expressed as: Thus assuming to apply two surface forces along the longitudinal direction on both ends of the composite, we can obtain their expressions by applying the constitutive laws for LL mode in 2D case as: And the coupling current induced through the piezoelectric layers can be calculated as [START_REF] Dong | Longitudinal and Transverse Magnetoelectric Voltage Coefficients of Magnetostrictive/Piezoelectric Laminate Composite: Theory[END_REF]:

𝑈 = ∫ 𝐸 1 𝑑𝑥
𝐹(0) = -
𝐼
𝑈 = ∫
𝐹(0) = -
𝐼

Figure 1 . 1

 11 Figure 1.1 Composite material conversion mechanisms

Figure 1 . 2

 12 Figure 1.2 Evolution of magnetoelectric materials[START_REF] Li | Study of magnetoelectric composites for sensor and transducer applications[END_REF]

Figure 1 . 3

 13 Figure 1.3 Schematic configurations for different connectivity types of two-phase composite materials: (a) 3-0 type particulate structure, (b) 2-2 type laminate structure, (c) 1-3 type rod matrix structure

Figure 1 . 5

 15 Figure 1.5 Magnetostriction creation with applied magnetic field under different temperature regions[START_REF] Li | Study of magnetoelectric composites for sensor and transducer applications[END_REF]

Figure 1 . 6

 16 Figure 1.6 Magnetization with positive magnetostriction under an uniaxial mechanical stress[START_REF] Nguyen | Modélisation par éléments finis de matériaux composites magnéto-électriques[END_REF]

Figure 1 . 7

 17 Figure 1.7 Example of curves of magnetic induction and magnetostriction of Terfenol-D [24]

Figure 1 . 8

 18 Figure 1.8 Magnetostrictive material conversion mechanisms

Figure 1 .

 1 Figure 1.10 Direct piezoelectric effect[START_REF] Li | Study of magnetoelectric composites for sensor and transducer applications[END_REF]

Figure 1 .

 1 Figure 1.11 Converse piezoelectric effect[START_REF] Li | Study of magnetoelectric composites for sensor and transducer applications[END_REF]

Figure 1 .

 1 Figure1.12 Polarization process[START_REF] Nguyen | Modélisation par éléments finis de matériaux composites magnéto-électriques[END_REF]

Figure 1 .

 1 Figure 1.13 Piezoelectric material conversion mechanisms

Figure 1 .

 1 Figure 1.14 Static magnetic field sensor[START_REF] Nguyen | Modélisation par éléments finis de matériaux composites magnéto-électriques[END_REF]

Figure 1 .

 1 Figure 1.15 Dynamic magnetic field sensor[START_REF] Nguyen | Modélisation par éléments finis de matériaux composites magnéto-électriques[END_REF]

Figure 1 .

 1 Figure 1.16 Energy source problems among current world [from Internet]

Figure 1 .

 1 Figure 1.17 Energy harvesting devices and systems [from Internet]

Figure 2 . 1

 21 Figure 2.1 Multiphysic coupling of a ME composite

Figure 2 . 2

 22 Figure 2.2 Local coordinates. The index 3 indicates the magnetization and polarization direction of the magnetostrictive and piezoelectric materials

Figure 2 . 3

 23 Figure 2.3 Illustration of the LT and TT modes for a trilayer ME composite

1 )

 1 𝑡 𝑝 +𝑚𝑡 𝑚 in which 𝑡 𝑚 , 𝑡 𝑝 represent the thicknesses of the single magnetostrictive and piezoelectric material layers in the laminated ME multilayer composite, and 𝑚 denotes the number of magnetostrictive layers, namely 𝑚 = 2 for the trilayer case study in this chapter.

Figure 2 . 4 Figure 2 . 5

 2425 Figure 2.4 Simulation results LT and TT modes for a MPM composite of Terfenol-D/PZT-5A/Terfenol-D

Figure 2 . 6 Figure 2 . 7

 2627 Figure 2.6 Illustrations of 3D composite structure in LT mode[START_REF] Dong | Longitudinal and Transverse Magnetoelectric Voltage Coefficients of Magnetostrictive/Piezoelectric Laminate Composite: Theory[END_REF]

  the static capacitance of the piezoelectric layer, 𝜑 𝑝 = 𝑤 𝑑 31𝑝 𝑠 11 𝐸 and for respectively the LT-mode and the TT-mode 𝜑 𝑚_𝐿𝑇 = 𝑤𝑡 𝑚 𝑑 33𝑚 𝑠 33 𝐻 and 𝜑 𝑚_𝑇𝑇 = 𝑤𝑡 𝑚 𝑑 31𝑚 𝑠 11 𝐻 . The MPM composite can been represented with the equivalent electric circuit shown in Figure 2.9 (Mason's model), in which 𝐹 1 , 𝐹 2 are 'mechanical voltages', 𝑍 1 , 𝑍 2 are 'electrical impedances' and 𝑢 1 ̇, 𝑢 2 ̇ are electrical currents.

Figure 2 . 9

 29 Figure2.9 Equivalent electric circuit model[START_REF] Dong | Longitudinal and Transverse Magnetoelectric Voltage Coefficients of Magnetostrictive/Piezoelectric Laminate Composite: Theory[END_REF]

Figure 2 .

 2 Figure 2.10 Simplified equivalent circuit model under free-free boundary condition[START_REF] Dong | Longitudinal and Transverse Magnetoelectric Voltage Coefficients of Magnetostrictive/Piezoelectric Laminate Composite: Theory[END_REF]

Figure 2 .

 2 Figure 2.11 shows for the LT-mode the dynamic response of the voltage ME coefficient for the same previous MPM composite (Terfenol-D/PZT-5A/Terfenol-D, 𝑄 𝑚𝑒𝑐ℎ = 200, 𝐿 = 14 𝑚𝑚, 𝑡 𝑝 = 1 𝑚𝑚, 𝑡 𝑚 = 1 𝑚𝑚).

Figure 2 .

 2 Figure 2.11 Dynamic response of a trilayer composite Terfenol-D/PZT-5A/Terfenol-D in LT mode

Figure 2 .

 2 Figure 2.12 Illustration of the 2D trilayer MPM composites in L-T mode

  (1-𝜗-𝑚𝜗) 𝔈 and 𝜗 are respectively the Young's module and the Poisson coefficient that are related to the compliance constants by 𝔈 =

Figure 2 .

 2 Figure 2.13 Illustration of the stress plan conditions

Figure 2 .

 2 Figure 2.14 shows the mentioned mechanical, electric and magnetic conditions in according to the global coordinates.

Figure 2 .

 2 Figure 2.14 Illustration of the mechanical, electric and magnetic conditions

Where

  [ℳ], [𝒞] , [𝒦] are respectively, the electro-magneto-mechanical mass, damping and stiffness matrix {𝒳} = [𝒖 𝑉 𝑎] 𝑇 , {ℱ} = [𝐹 𝑄 𝑛 𝐽 𝑠 ] 𝑇 ,𝐾 𝑝𝑢 = 𝐾 𝑢𝑝 𝑡 , 𝐾 𝑎𝑢 = 𝐾 𝑢𝑎 𝑡 , 𝐶 𝑢𝑢 = 𝛽𝐾 𝑢𝑢 + 𝛼𝑀.

2 ( 2 - 4 - 30 ) 1 𝜔

 224301 relations of 𝜔 𝑛 = √𝑘/𝑀 and 𝜔 𝑟 = √𝑘 𝑀 In considering a low modal damping coefficient 𝜍, we approximate 𝜔 𝑟 ≈ 𝜔 𝑛 , where 𝜔 𝑟 = 2𝜋𝑓 𝑟 can be analytically estimated in LT mode by :In the case of a ME composite Terfenol-D/PZT-5A/Terfenol-D with the dimensions: 𝐿 = 14 𝑚𝑚, 𝑡 𝑝 = 1 𝑚𝑚, 𝑡 𝑚 = 1 𝑚𝑚, the theoretical frequency resonance is close to 75 kHz. Thus, to impose a mechanical quality factor 𝑄 𝑚𝑒𝑐ℎ close to 200, the 𝛽 =

Figure 2 .Figure 2 .

 22 Figure 2.15 Illustration in 3D of the domain problem with the boundary conditions

. 16 .

 16 In fact, the magnetic flux Φ defined by(2-4-34) can be degenerated in 2D as 𝑎 0 = 𝐵 0 ∆ 𝑧 with ∆ 𝑧 = (𝑧 𝑚𝑎𝑥 -𝑧 𝑚𝑖𝑛 ) and 𝐵 0 = 𝐻 𝑎𝑐 𝜇 0 where 𝜇 0 represents the vacuum permeability.

  shows the structure for which the results will be presented. As previously, we consider the case Terfenol-D/PZT-5A/Terfenol-D where their magneto-electro-elastic proprieties are given in the Appendix A and the dimensions: 𝐿=14 mm, 𝑡 𝑝 =1 mm, 𝑡 𝑚 =1 mm, the theoretical frequency resonance is close to 75 kHz, and to impose a mechanical quality factor 𝑄 𝑚𝑒𝑐ℎ close to 200, the 𝛽= 1 𝜔 𝑟 𝑄 𝑚𝑒𝑐ℎ ≈1.10× 10 -8 (𝑠). The magnetic, mechanical and electric distribution fields for the LT-mode and TT-mode in the linear and static case are presented inFigures 2.18

  to 2.20.

Figure 2 .Figure 2 .Figure 2 .

 222 Figure 2.17 The studied ME composite

Figure 2 .

 2 Figure 2.20 Electrical potential scalar and electrical field distributions

  usually simplified with the isotropic case in assuming that 𝑐 = 0 and 𝑴 𝒓𝒆𝒗 = 0[START_REF] Marcelo | Coupled Structural-Magnetic Strain and Stress Model for Magnetostrictive Transducers[END_REF][START_REF] Calkins | Energy-based hysteresis model for magnetostrictive transducers[END_REF][START_REF] Du | Dynamic Jiles-Atherton Model for Determining the Magnetic Power Loss at High Frequency in Permanent Magnet Machines[END_REF]

Figure 2 .

 2 [START_REF] Cedell | New magnetostrictive alloy for rapid conversion electric energy to mechanical motion[END_REF] displays an example of the simulation results of the magnetisation and the magnetic induction (𝑩 = 𝜇 𝑜 (𝑯 + 𝑴)) obtained for Terfenol-D[START_REF] Calkins | Energy-based hysteresis model for magnetostrictive transducers[END_REF] in which its material parameters are given in

Figures 2 .Figure 2 .Figure 2 .

 222 Figures 2.22 Magnetisation and the magnetic induction in function of the applied magnetic field and under different pre-stress 𝑻 𝒐

Figure 2 .

 2 Figure 2.25 Piezomagnetic coefficient d 33 with the applied magnetic field in Oe unity and under different prestress 𝑻 𝒐

Figure 2 .

 2 Figure 2.26 Flowchart of the piecewise-linear solution procedure

Figure 2 .

 2 Figure 2.27 shows the efficient of the proposed modelling considering 𝑛 = 6 points in the curve.

Figure 2 .

 2 Figure 2.27 Mean curve H(B)=v(B).B reconstructed with polynomial model

Figure 2 .

 2 Figure 2.29 Illustration of the dynamic small signal regime

( 2 - 4 - 90 ).

 2490 In this case, the current 𝐼 will be also normalized by 𝐼 ̅ = 𝜕 𝑡 {𝐶 ̅ 𝑉} = 𝜕 𝑡 { To conserve the physical result 𝑉 = 𝑍𝐼 = 𝑍𝜕 𝑡 {𝐶𝑉}, it is necessary to normalize also the impedance 𝑍 by 𝑍 ̅ = 𝑍𝐿 𝑥 . is thus solved with the following matrices and vectors:

Figure 2 .

 2 Figure 2.30 compares under open circuit condition in considering (Z= 1 MΩ) the dynamic response between the FEM and the equivalent circuit method presented in the first part in this chapter 2 for the trilayer case Terfenol-D/PZT-5A/Terfenol-D in the L-T mode with the configuration showed as the first one of the four in Fig.2.31. The simulations have been performed with an externally applied small signal ac magnetic field δH ac = 1Oe around an optimal magnetization bias close to H dc =200 Oe with the sizes 𝐿=14 mm, 𝑡 𝑝 =1 mm, 𝑡 𝑚 =1 mm and the parameter values in Appendix A. We can notice that the frequency resonance of 75 kHz of the FEM result corresponds to theoretical value 74.3 kHz given by (2-4-31).

Figure 2 . 2 . 4 . 4 . 3 . 2 Figure 2 .Figure 2 . 2 . 4 . 4 . 3 . 3 Figure 2 .

 22443222244332 Figure 2.30 Dynamic responses for the trilayer in L-T mode

Figure 2 .

 2 Figure 2.33 ME coefficient in function of frequency under different load resistances

Figure 3 . 1 .

 31 The generated AC voltage is rectified into DC voltage by the full-wave rectifying diode bridge, and the DC voltage is stored temporarily (by a capacitor) in a storage stage. The Figure3.2 shows an electrical simulation performed with the Simulink tool "Simscape Power Systems" of MATLAB ® in considering an AC input voltage working at 75 kHz (in the LT mode) with 10 volt peak-to-peak, a rectified diode bridge composed of diodes (Forward voltage at 0.6 V) and a storage stage composed by a capacitor of 22 µF. We notice that the DC output voltage reaches promptly (in 0.6 ms) 4 V that corresponds to the typical value required to power supply (active mode) the micro Integrated Circuit (IC) chip or Wireless Sensor Nodes (WSN)[START_REF] Beeby | Energy Harvesting for Autonomous Systems[END_REF][START_REF] Tan | Energy Harvesting Autonomous Sensor Systems: Design, Analysis, and Practical Implementation[END_REF].

Figure 3 . 1

 31 Figure 3.1 Illustration of the energy harvesting process

Figure 3 . 2 Figure 3 . 3

 3233 Figure 3.2 Electrical simulation performed with Simulink tool Simscape Power Systems of MATLAB

Figure 3 .

 3 Figure 3.4 shows the FEM simulation results of the active output deliverable power (defined as 𝑃 = 𝑉 𝑜𝑢𝑡 2 2𝑅 ⁄ ) from the Terfenol-D/PZT-5A/Terfenol-D investigated in chapter 2 with an input δ𝐻 𝑎𝑐 of 1Oe for different modes (respectively at their resonance frequencies) as function of the resistive load 𝑅 connected between both electrodes. This resistive load 𝑅 represents symbolically the real part of the load impedance of the rectified and stored circuit. Although the TL-mode and LL-mode have significant ME voltage coefficients in respect to the LT-mode, their deliverable powers are, as shown in Figure 3.4 (c) and (d), very low (few µW) due to the high internal impedance close to 40 MOhm.Thus, it appears that the adequate solution for a duty cycle of 1% is the LT-mode since as shown in Figure3.4 (a) and (b) the deliverable power reaches 800 µW (>500 µW) for optimal impedance close to 18 kOhm.

Figure 3 . 4

 34 Figure 3.4 Output voltage and deliverable power as function of the resistive load in according to different modes

Figure 3 . 5

 35 Figure 3.5 Equivalent circuit magnetoelectric generator and Thevenin model. 𝑍 𝑇 represents the complex internal impedance 𝑍 𝑚 ′//𝑍 𝐶 𝑜

Figure 3 . 6

 36 Figure 3.6 Internal impedance Z T versus frequency under different 𝛾 values

Figure 3 . 7

 37 Figure 3.7 Comparison between the FEM simulation and the extraction method

Figure 3 . 8

 38 Figure 3.8 Equivalent circuit magnetoelectric generator in Norton model.

{

  𝑀𝑢̈+ 𝐶 𝑢𝑢 𝑢̇+ 𝐾 𝑢𝑢 𝑢 + 𝐾 𝑢𝑝 𝑉 + 𝐾 𝑢𝑎 𝑎 1 = 0 𝐾 𝑝𝑢 𝑢 + 𝐾 𝑝𝑝 𝑉 = 0 𝐾 𝑎𝑢 𝑢 + 𝐾 𝑎𝑎 𝑎 1 = 𝐽 𝑠 (3-2-21) Thus, in combining the expression of 𝑎 1 = 𝐾 𝑎𝑎 -1 (𝐽 𝑠 -𝐾 𝑎𝑢 𝑢) with 𝑢 = -𝐾 𝑝𝑢 -1 𝐾 𝑝𝑝 𝑉 in the first equation we obtain the equivalent system: 𝑀𝑢̈+ 𝐶 𝑢𝑢 𝑢̇+ 𝐾 𝑢𝑢 𝑢 + (𝐾 𝑢𝑝 + 𝐾 𝑢𝑝 )𝑉 = 𝐹 𝑒𝑓𝑓 (3-2-22)with 𝐾 𝑢𝑝 = 𝐾 𝑢𝑎 𝐾 𝑎𝑎 -1 𝐾 𝑎𝑢 𝐾 𝑝𝑢 -1 𝐾 𝑝𝑝 and where 𝐹 𝑒𝑓𝑓 = -𝐾 𝑢𝑎 𝐾 𝑎𝑎 -1 𝐽 𝑠 represents the effective mechanical external force applied due to the external magnetic excitation 𝐽 𝑠 .The establishment of energy equilibrium can be obtained in multiplying both sides of (3-2-22) by the displacement velocity 𝑢̇ and integrating over the time variable[START_REF] Guyomar | Toward energy harvesting using active materials and conversion improvement by nonlinear processing[END_REF][START_REF] Li | Structural vibration damping with synchronized energy transfer between piezoelectric patches[END_REF][START_REF] Corcolle | Modeling of a beam structure with piezoelectric materials: introduction to SSD techniques[END_REF]:

  𝑢 2 (𝑡) and the transmitted energy ∫ (𝐾 𝑢𝑝 + 𝐾 𝑢𝑝 )𝑉𝑢ṫ 𝑜 𝑑𝑡.

1 2 1 2

 11 𝐶 𝑜 𝑉 2 corresponding to the energy stored on the piezoelectric capacitance 𝐶 𝑜 and the energy absorbed by the electrical load ∫ 𝑉𝐼 𝑡 𝑜 𝑑𝑡. ∫ (𝐾 𝑢𝑝 + 𝐾 𝑢𝑝 )𝑉𝑢ṫ 𝑜 𝑑𝑡 = 𝐶 𝑜 𝑉 2 + ∫ 𝑉𝐼 𝑡 𝑜 𝑑𝑡 (3-2-25)

Case 1 :

 1 SSDI technique with L and R in parallel Case 2: SSDI technique with L and R in series Figure 3.9 SSDI techniques To implement the SSDI in the FEM multiphysics model, it is necessary to write the differential circuit equation in the transient system. For that, the relation 𝑖(𝑡) = 𝜕 𝑡 𝑄 between the current and the charge is considered. In case 1, the differential circuit equation is introduced in the system through considering the current equation 𝑖(𝑡) = 𝑖 𝑅 (𝑡) + 𝑖 𝐿 (𝑡) , with 𝑖 𝑅 = 𝑈 𝑅 and 𝜕 𝑡 𝑖 𝐿 = 𝑈 𝐿 . Thus, in considering 𝜕 𝑡 𝑖 =

Figure 3 .

 3 Figure 3.10 Transient response of output voltage and deliverable energy[START_REF] Nguyen | Finite Element Harmonic Modeling of Magnetoelectric Effect[END_REF]

Figure 3 .

 3 Figure 3.11 Standard configurations in series and parallel electrical connections

Figure 3 .

 3 Figure 3.12 shows an illustration of the mesh structure employed in the FEM simulation.

Figure 3 .

 3 Figure 3.12 FEM meshing for a multilayer composite

Figure 3 .Figure 3 .

 33 Figure 3.13 Equivalent circuit model in LT and LL modes

Figure 3 . 15 :

 315 Figure 3.15: Results comparison in LT mode.

Figure 3 .Figure 3 .

 33 Figure 3.16 shows a comparison between the distributions of the electric scalar potential for the trilayer case and for the multilayer case with series and parallel electrical connections.

Figure 3 .

 3 Figure 3.17 Results comparison in LL mode. Configuration 1 is the standard one and configuration 2 is the one with an alternative opposing polarization between each neighboring piezoelectric layer

Figure 3 .

 3 Figure 3.19 Series electrical connection in according to the configuration 2 in LL mode

Figure 4 . 1

 41 Figure 4.1 Illustration of the ME composite energy transducer in a catheter and inside the human body

Figure 4 . 2

 42 Figure 4.2 Ferrite solenoid proposed in [87] to be used as external reader

Figure 4 . 3 .

 43 As shown in Appendix A, the material characteristics of the PZT-5H are close to the PZT-5A. The layers are bounded with a strong epoxy resin and the PZT-5H is polarized transversally while the Terfenol-D is magnetized longitudinally. In this way, the composite ME is employed in L-T mode. The PZT-5H has been chose longer than the Terfenol-D in order to perform the connection between both electrodes of the PZT-5H. Unfortunately, the commercial providers of the PZT-5H (of The Roditi International Corporation Ltd) and the Terfenol-D (of ETREMA PRODUCTS INC, and Subsidiary of Edge Technologies Inc.) do not provide values of constants for their electrical, magnetic and elastic material characteristics.

Figure 4 . 3

 43 Figure 4.3 The considered bilayer ME composite

Figure 4 .

 4 Figure 4.4 shows the FEM simulation results for the considered ME composite in according to the model of the Terfenol-D shown in Chapter 2 and with the material characteristics of the PZT-5H presented in Annexe A. The resonance frequency is close to 70 kHz with a maximum ME coefficient reaching 20 V/Oe. The first local peak around 60 kHz is a parasite resonance mode caused by the difference of length between the PZT-5H (20 mm) and the Terfenol-D (15 mm). It must be noticed that the effect of the epoxy resin was not considered in the simulation.

Figure 4 . 4 FEM

 44 Figure 4.4 FEM simulation results of the bilayer ME composite

  (a) Experimental set-up (b) FEM simulation of the configuration Figure 4.5 Realisation of the dc magnetic field source

Figure 4 .

 4 Figure 4.6 shows the magnetic field 𝐻 𝑑𝑐 in Oe versus the distance d between both sources. We can note that the magnitudes of the magnetic field are sufficient to magnetize the Terfenol-D (usually ≥ 200 Oe as shown in chapter 2) for distances between 40 mm and 80 mm.

Figure 4 . 6

 46 Figure 4.6 Evolution of the dc magnetic field 𝐻 𝑑𝑐 in function of the distance between both sources

Figure 4 . 7

 47 Figure 4.7 Support of the solenoid coil for the ME composite

Figure 4 . 8

 48 Figure 4.8 Experimental set-up

Figure 4 .

 4 Figure 4.9 shows the ME coefficient measured at low frequency (1 kHz) under the variation of the distance d between both magnetic sources. To respect the definition of the ME coefficient namely V/Oe, for each measure, the output voltage V out must be rectified in normalizing the magnetic field inside the solenoid coil (estimated by the crossing current), in respect to 1 Oe. The maximal ME coefficient is obtained for a distance close to 55 mm. In according to Figure 4.6, this distance corresponds to a magnetization bias 𝐻 𝑑𝑐 close to 350 Oe.

Figure 4 . 9

 49 Figure 4.9 ME coefficient at low frequency in function of the distance between both sources

Figure 4 .

 4 Figure 4.10 ME coefficient in function of frequency under an ac magnetic field of 1 Oe

Figure 4 .

 4 Figure 4.11 Principle to determinate the value of the internal impedance

Figure B. 1

 1 Figure B.1 Magnetostriction curves by experiments and the SS model

  modulus 𝑚 and magnetoelastic modulus 𝑟 can be obtained as of comparisons between the results from experiments and the HT model (Fig.B.2)

Figure B. 2

 2 Figure B.2 Magnetostriction curves by experiments and the HT model

  equations can be obtained by using the Legendre transform as: of comparisons between the results from experiments and the DDS model (Fig.B.3)

Figure B. 4

 4 Figure B.4 Comparison of the predicted magnetostrictive strain curves with the experimental data for a Terfenol-D rod under various compressive pre-stress levels (Hysteresis loops: experimental; solid lines: predicted)

5 .

 5 we have 𝐹(𝑋 𝑘+1 ) -𝐹(𝑋 𝑘 ) = 𝜕𝐹(𝑋) 𝜕𝑋 | 𝑋=𝑋 𝑘 • (𝑋 𝑘+1 -𝑋 𝑘 ), and by assuming 𝐹(𝑋 𝑘+1 ) = 0, we can obtain the value of the vector X in the next iteration step by calculating as 𝑋 𝑘+1 = 𝑋 𝑘 -( 𝜕𝐹(𝑋) 𝜕𝑋 | 𝑋=𝑋 𝑘 ) -1 • 𝐹(𝑋 𝑘 )Step 3: Check if the relative error expression in this iteration step satisfies the condition 𝜀 = If not, return to Step 2 to calculate the next vector X value otherwise go toStep 4.

Figure C. 1

 1 Figure C.1 An image indicating the procedure of iteration in the problem

3 . 2 ,

 32 we can obtain the wave equation of the structure based on the mechanical equilibrium equation (2-2-10) as derived in the Chapter 2.3.2: 𝑢 𝑥 indicates the mechanical displacement of the composite along the longitudinal (L) direction (x), ω the angular frequency of the excitation signal, k the wave number of the composite longitudinal vibration, mechanical wave velocity along L direction and 𝑛 = 𝑚𝑡 𝑚 ′ 𝑚𝑡 𝑚 ′ +(𝑚-1)𝑡 𝑝 ′ , 𝜌̅ = 𝑛𝜌 𝑚 + (1 -𝑛)𝜌 𝑝 , i ,j=1 or 3 indicate the magnetization and polarization direction is transverse or longitudinal, respectively. The wave equation above has the general solution as: 𝑢 1 ̇= 𝑗𝜔𝑢(0), 𝑢 2 ̇= 𝑗𝜔𝑢(𝐿), and 𝑗 = √-1And we can obtain the strains on the boundary of the composite:

𝑠 11 𝐸

 11 By applying the free-free boundary condition (F(0)=F(L)=0) and the open circuit condition (𝐼 𝑝 = 0) we can derive the expression of ME coefficient in LT mode as:

  1 ̇-𝑢 2 ̇)) = -𝜑 𝑝 ′(𝑢 1 ̇-𝑢 2 ̇) + 𝑗𝜔𝐶 𝑇 𝑈 (D-12)By applying the free-free boundary condition (F(0)=F(L)=0) and the open circuit condition (𝐼 𝑝 = 0) we can derive the expression of ME coefficient in LL mode in series connectivity as:

  ̇-𝑢 2 ̇)) = -𝜑 𝑝 ′(𝑢 1 ̇-𝑢 2 ̇) + 𝑗𝜔𝐶 𝑇 𝑈 (D-17)By applying the free-free boundary condition (F(0)=F(L)=0) and the open circuit condition (𝐼 𝑝 = 0) we can derive the expression of ME coefficient in LL mode in series connectivity as:

Table 1 .

 1

2 Microscopic origin of the piezoelectricity

[START_REF] Nguyen | Modélisation par éléments finis de matériaux composites magnéto-électriques[END_REF]

  𝑇 𝑥𝑥 𝑇 𝑦𝑦 𝑇 𝑧𝑧 𝑇 𝑦𝑧 𝑇 𝑥𝑧 𝑇 𝑥𝑦 ] 𝑡

							(2-2-14)
	𝑺 = [𝑆 𝑥𝑥 𝑆 𝑦𝑦 𝑆 𝑧𝑧 2𝑆 𝑦𝑧 2𝑆 𝑥𝑧 2𝑆 𝑥𝑦 ] 𝑡	(2-2-15)
	They are linked by
			𝑻 = 𝑐 𝑺				(2-2-16)
	where 𝑐 is the stiffness tensor in N•m -2 . For example, in the orthotropic material studied in this
	dissertation, the tensor can be expressed in matrix form as
		𝑐 11 𝑐 12 𝑐 13 0	0	0
		𝑐 12 𝑐 22 𝑐 23 0	0	0
	𝑐 =	𝑐 13 𝑐 23 𝑐 33 0 0 0 0 𝑐 44 0 0	0 0	(2-2-17)
		0	0	0	0 𝑐 55 0
	[	0	0	0	0	0 𝑐 66 ]

  2𝜌 𝑚 𝐴 1 +𝜌 𝑝 𝐴 2

	𝜌̅	𝜕 2 𝑢 𝜕𝑡 2 = 𝑛	𝜕𝑇 33 𝑚 𝜕𝑧 + (1 -𝑛)	𝜕𝑇 11 𝑝 𝜕𝑧 =	𝑛 𝑠 33 𝐻	𝜕𝑆 33 𝑚 𝜕𝑧 +	1-𝑛 𝑠 11 𝐸	𝜕𝑆 11 𝑝 𝜕𝑧 =	𝑛 𝑠 33 𝐻	𝜕 2 𝑢 𝜕𝑧 2 +	1-𝑛 𝑠 11 𝐸	𝜕 2 𝑢 𝜕𝑧 2	(2-3-15)
	And according to the propagation principle the average sound velocity of the material can be
	expressed respectively for the LT-mode and the TT-mode as 𝑣 𝐿𝑇 ̅̅̅̅̅ = √ 𝑠 33 𝑛 𝐻 + 𝜌 ̅ 1-𝑛 𝑠 11 𝐸	and 𝑣 𝐿𝑇 ̅̅̅̅̅ = √ 𝑠 11 𝑛 𝐻 + 𝜌 ̅ 1-𝑛 𝑠 11 𝐸

𝐴

= 𝑛𝜌 𝑚 + (1 -𝑛)𝜌 𝑝 and by substituting the constitutive laws we can obtain:

  For the LT (𝛼 31 𝑉 ) 𝑙𝑜𝑤 and TT (𝛼 33 𝑉 ) 𝑙𝑜𝑤 modes we obtain with 𝑡 𝑝 = 𝑡 𝑙𝑎𝑚 (1 -𝑛) where 𝑡 𝑙𝑎𝑚 = 𝑡 𝑝 + 2𝑡 𝑚 is the total thickness of the ME composite.

						(𝛼 31 𝑉 ) 𝑙𝑜𝑤 = | 𝛾 1 (𝑑 31 2𝑑 33 𝑚 𝑑 31 𝑝 (𝑛-1)𝑛𝑡 𝑙𝑎𝑚 𝑝 ) 2 -(𝑛(𝑠 11 𝑝 -𝑠 33 𝑚 )+𝑠 33 𝑚 )𝜀 33 𝑇 |	(2-3-30)
	with γ 1 = (𝑛 -1)	𝑚 𝑠 33 𝑠 11 𝑝
						(𝛼 33 𝑉 ) 𝑙𝑜𝑤 = |	2𝑑 31 𝑚 𝑑 31 𝑝 (𝑛-1)𝑛𝑡 𝑙𝑎𝑚 2 γ 2 -(𝑛(𝑠 11 𝑝 ) (𝑑 31 𝑝 -𝑠 11 𝑚 )+𝑠 11 𝑚 )𝜀 33 𝑇 |	(2-3-31)
	with γ 2 = (𝑛 -1)	𝑚 𝑠 11 𝑠 11 𝑝
												33 𝑚 𝑠 11 𝐸 𝑑 31 𝑝 𝑠 33 𝐻	|	(2-3-28)
	In the low frequency case we have cot ( 𝑘𝐿 2 ) ≈	𝑘𝐿 2 , in other words the impedance 𝑍 is approximated as
	𝑍 ≈	-𝑗	1 2 𝜌 ̅𝑣 ̅𝐴 𝑘𝐿 2	=	-𝑗𝜌 ̅𝑣 ̅𝐴 𝑘𝐿 .
	Thus in replacing 𝑗𝑍𝐶 𝑜 𝜔 by	𝜌 ̅𝑣 ̅ 2 𝐴 𝐿 𝐶 𝑜 . in the expression (2-3-25) we find the expression of the voltage
	magnetoelectric coefficients in low frequency as:
						(𝛼 31 𝑉 ) 𝑙𝑜𝑤 = |	𝜑 𝑝	𝜑 𝑝 𝜑 𝑚 𝐿 2 𝐿+𝜌 ̅𝑣 ̅ 2 𝐴𝐶 𝑜	| = | 𝜑 𝑝	𝜑 𝑝 𝜑 𝑚 2 + 𝜌 ̅𝑣 ̅ 2 𝐴 𝐿	𝐶 𝑜	| = |	𝜑 𝑝	𝜑 𝑝 𝜑 𝑚 𝐿 2 𝐿+𝜌 ̅𝑣 ̅ 2 𝐴𝐶 𝑜	|	(2-3-29)

Table 2

 2

	Parameter	Value
	𝝀 𝒔	995 ppm
	𝑻 𝒐	1 ksi≈6.9 MPa
	𝑴 𝒔	7,5x 10 5 A/m
	𝒂	7012 A/m
	𝒄	0.18
	𝒌	3942 A/m
	𝜶 ̃	-1,17 10 -2
	Table 2.1 Material parameters for a rod Terfenol-D [64]
	Figures 2.

.1.

Figure 2.21 Magnetization and magnetic induction under an applied magnetic field

Table 3 .

 3

	Duty cycle	10%	1%	0.1%
	〈𝑃〉	5 mW	500 µW	62 µW

1 gives the average power for three different duty cycles.

Table 3 .

 3 1. Average power consumed by a Wireless Sensor Node

Table 3 .

 3 2. Average power in considering an Impedance Matching Circuit

	cycle	10%	1%	0.1%
	〈𝑃〉	8 mW	800 µW	92 µW

Table 3 .

 3 .3. 𝐶 𝑜 represents the clamped capacitor of the piezoelectric layer in the trilayer case. 3. Total equivalent clamped capacitor 𝐶 𝑇 in LT and LL modes

	Electrical connection type	LT mode, 𝐶 𝑜 =	𝜀 33 𝑆 𝑙𝑤 𝑡 𝑝	LL mode, 𝐶 𝑜 =	𝑡 𝑝 𝑤 𝑙𝛽 33 ̅̅̅̅̅
	Series			𝐶 𝑇 =	𝛿𝐶 𝑜 (𝑚 -1)	𝐶 𝑇 =	𝐶 𝑜 𝛿(𝑚 -1)
	Parallel			𝐶 𝑇 = (𝑚 -1)𝛿𝐶 𝑜	𝐶 𝑇 =	(𝑚 -1) 𝛿	𝐶 𝑜
	where 𝜀 33 𝑆 = 𝜀 33 𝑇 -	𝑑 31 𝑝 2 𝑠 11 𝐸 , 𝛽 33 ̅̅̅̅ =	𝐸 𝑠 33 𝑇 𝑠 33 𝜀 33 𝐸 -𝑑 33 𝑝 2

Table 3 .

 3 .[START_REF] Astrov | The magnetoelectric effects[END_REF] gives the adjustment of expressions of 𝜑 𝑚 , 𝜑 𝑝 . The expressions of 𝑍 𝑚 are given in the Appendix D. 4. Adjustment on the magnetomechanical and electromechanical coupling coefficients in multilayer case

	Mode	Series electrical connection	Parallel electrical connection
	LT	𝜑 𝑝 ′ = 𝜑 𝑝 = 𝑤	𝑑 31𝑝 𝑠 11 𝐸	𝜑 𝑝 ′ = 𝜑 𝑝 = 𝑤	𝑑 31𝑝 𝑠 11 𝐸

  By utilising the Legendre transform relations 𝜀 𝑖𝑗 = -(𝜕𝐺/ 𝜕𝜎 𝑖𝑗 )| 𝐻 and 𝐵 𝑘 = -(𝜕𝐺/ 𝜕𝐻 𝑘 )| 𝜎 from equation (B-4) and truncating the polynomial expansion according to the conclusions from the experimental results of a giant magnetostrictive compound material Tb 0.27 Dy 0.73 Fe 1.95 that the strain is the function of the magnetic field with only even powers and thus the material constant tensors of odd orders, such as the three order piezomagnetic coefficient tensor, must vanish from the series expansion, the constitutive relations can be derived as:𝜀 𝑖𝑗 =𝑠 𝑖𝑗𝑘𝑙 𝜎 𝑘𝑙 + 𝑚 𝑖𝑗𝑘𝑙 𝐻 𝑘 𝐻 𝑙 + 𝑟 𝑖𝑗𝑘𝑙𝑚𝑛 𝜎 𝑘𝑙 𝐻 𝑚 𝐻 𝑛 (B-5) 𝐵 𝑘 = 𝜇 𝑘𝑙 𝐻 𝑙 + 𝑚 𝑘𝑙𝑚𝑛 𝜎 𝑚𝑛 𝐻 𝑙 + 𝑟 𝑘𝑙𝑚𝑛𝑝𝑞 𝜎 𝑚𝑛 𝜎 𝑝𝑞 𝐻 𝑙 (B-6)

		0 +	𝜕𝐺 𝜕𝜎 𝑖𝑗	∆𝜎 𝑖𝑗 +	𝜕𝐺 𝜕𝐻 𝑘	Δ𝐻 𝑘 +	1 2	𝜕 2 𝐺 𝜕𝜎 𝑖𝑗 𝜕𝜎 𝑘𝑙	∆𝜎 𝑖𝑗 ∆𝜎 𝑘𝑙 +	1 2	𝜕 2 𝐺 𝜕𝜎 𝑖𝑗 𝜕𝐻 𝑘	∆𝜎 𝑖𝑗 ∆𝐻 𝑘 +	1 2	𝜕 2 𝐺 𝜕𝐻 𝑙 𝜕𝐻 𝑘	∆𝐻 𝑙 ∆𝐻 𝑘 +
	1 3!	𝜕 3 𝐺 𝜕𝜎 𝑖𝑗 𝜕𝜎 𝑘𝑙 𝜕𝜎 𝑚𝑛	∆𝜎 𝑖𝑗 ∆𝜎 𝑘𝑙 ∆𝜎 𝑚𝑛 +	1 3!	𝜕 3 𝐺 𝜕𝜎 𝑖𝑗 𝜕𝜎 𝑘𝑙 𝜕𝐻 𝑚	∆𝜎 𝑖𝑗 ∆𝜎 𝑘𝑙 ∆𝐻 𝑚 +	1 3!	𝜕 3 𝐺 𝜕𝜎 𝑖𝑗 𝜕𝐻 𝑘 𝜕𝐻 𝑙	∆𝜎 𝑖𝑗 ∆𝐻 𝑘 ∆𝐻 𝑙 +
													1 3!	𝜕 3 𝐺 𝜕𝐻 𝑚 𝜕𝐻 𝑘 𝜕𝐻 𝑙	∆𝐻 𝑚 ∆𝐻 𝑘 ∆𝐻 𝑙 + ⋯	(B-4)

  𝑆 𝑖𝑗𝑘𝑙 𝜎 𝑘𝑙 + (𝑚 𝑖𝑗𝑘𝑙 -

	𝑆 𝑖𝑗𝑚𝑛 (𝑚) 𝜎 𝑚𝑛 𝑀 𝑠 2	𝛿 𝑘𝑙 ) 𝑀 𝑘 𝑀 𝑙		(B-23)
	𝐻 𝑘 = 𝑓 𝑘 -1 (𝑀 𝑙 ) -𝜇 0 -1 (2𝑚 𝑖𝑗𝑘𝑙 𝜎 𝑖𝑗 -	𝑆 𝑖𝑗𝑚𝑛 (𝑚) 𝜎 𝑖𝑗 𝜎 𝑚𝑛 𝑀 𝑠 2	𝛿 𝑘𝑙 ) 𝑀 𝑙	(B-24)

  𝜎 𝑠 = 𝜆 𝑠 𝐸 𝑠 𝐸 (𝐸 𝑠 -𝐸) ⁄ where 𝐸 𝑠 denotes the intrinsic Young's Modulus, 𝑘 = 𝜒 𝑚 𝑀 𝑠 ⁄ where 𝜒 𝑚 is the initial susceptibility and 𝑀 𝑠 is the saturation magnetization, 𝑓(𝑥) from the inverse 𝑓 𝑘 -1 (𝑀 𝑙 ) has various choices, such as the hyperbolic tangent function 𝑓(𝑥) = 𝑡𝑎𝑛ℎ (𝑥), which is used in the D-H model, and the Langevin function 𝑓(𝑥) = 𝑐𝑜𝑡ℎ(𝑥) -1/𝑥, which is adopted in the Zheng-Liu model due to its clear physical background based on the Boltzmann statistics[START_REF] Calkins | Energy-based hysteresis model for magnetostrictive transducers[END_REF] and possibility of better description of the magnetization curve.

	1 𝐸 [	1+𝜈 2 (𝛿 𝑖𝑘 𝛿 𝑗𝑙 + 𝛿 𝑖𝑙 𝛿 𝑗𝑘 ) -𝜈𝛿 𝑖𝑗 𝛿 𝑘𝑙 ]	(B-25)
	𝑆 𝑖𝑗𝑚𝑛 (𝑚) =	𝜆 𝑠 𝜎 𝑠	[	3 4 (𝛿 𝑖𝑚 𝛿 𝑗𝑛 + 𝛿 𝑖𝑛 𝛿 𝑗𝑚 ) -	1 2 𝛿 𝑖𝑗 𝛿 𝑚𝑛 ]	(B-26)
	𝑚 𝑖𝑗𝑘𝑙 =	𝜆 𝑠 𝑀 𝑠 2 [	3 4 (𝛿 𝑖𝑘 𝛿 𝑗𝑙 + 𝛿 𝑖𝑙 𝛿 𝑗𝑘 ) -	1 2	𝛿 𝑖𝑗 𝛿 𝑘𝑙 ]	(B-27)
	𝑓 𝑘 -1 (𝑀 𝑙 ) =	1 𝑘𝑀 𝑓 -1 ( 𝑀 𝑠 𝑀	) 𝛿 𝑘𝑙 𝑀 𝑙	(B-28)

where E is the Young's Modulus, 𝜈 the Poisson's Ratio, 𝜆 𝑠 the saturation magnetostrictive strain,

  11 𝑥 1 + 𝑎 12 𝑥 2 + ⋯ + 𝑎 1𝑛 𝑥 𝑛 -𝑏 1 𝑎 21 𝑥 1 + 𝑎 22 𝑥 2 + ⋯ + 𝑎 2𝑛 𝑥 𝑛 -𝑏 2 ⋯ 𝑎 𝑛1 𝑥 1 + 𝑎 𝑛2 𝑥 2 + ⋯ + 𝑎 𝑛𝑛 𝑥 𝑛 -𝑏 𝑛

		𝐹 1 (𝑥 1 , 𝑥 2 , ⋯ , 𝑥 𝑛 )
	] = [	𝐹 2 (𝑥 1 , 𝑥 2 , ⋯ , 𝑥 𝑛 ) ⋯ ]
		𝐹 𝑛 (𝑥 1 , 𝑥 2 , ⋯ , 𝑥 𝑛 )

  Due to the difficulty to obtain the exact expressions of the terms Thus for each 1 ≤ 𝑖, 𝑗, 𝑚 ≤ 𝑛 , we can use 𝑋 𝑘,𝑗1 = [ 𝑥 1 𝑥 2 … 𝑥 𝑗 + ∆ … 𝑥 𝑛 ] 𝑇 and 𝑋 𝑘,𝑗2 = [ 𝑥 1 𝑥 2 … 𝑥 𝑗 -∆ … 𝑥 𝑛 ] 𝑇 to calculate the subsequent variables along the chain as: 𝑋 𝑘,𝑗1 , 𝑋 𝑘,𝑗2 → 𝐵 ̅ (𝑋 𝑘,𝑗1 ), 𝐵 ̅ (𝑋 𝑘,𝑗2 ) → 𝑞 (𝐵 ̅ (𝑋 𝑘,𝑗1 )) , 𝑞 (𝐵 ̅ (𝑋 𝑘,𝑗2 )) , 𝑣 𝑆 (𝐵 ̅ (𝑋 𝑘,𝑗1 )) , 𝑣𝑆 (𝐵 ̅ (𝑋 𝑘,𝑗2 )) → [ 𝐾 𝑢𝑢 (𝑞 (𝐵 ̅ (𝑋 𝑘,𝑗1 )) , 𝑣 𝑆 (𝐵 ̅ (𝑋 𝑘,𝑗1 ))) 𝐾 𝑢𝑎 (𝑞 (𝐵 ̅ (𝑋 𝑘,𝑗1 ))) 𝐾 𝑎𝑢 (𝑞 (𝐵 ̅ (𝑋 𝑘,𝑗1 ))) 𝐾 𝑎𝑎 (𝑣 𝑆 (𝐵 ̅ (𝑋 𝑘,𝑗1 ))) 𝑞 (𝐵 ̅ (𝑋 𝑘,𝑗2 )) , 𝑣 𝑆 (𝐵 ̅ (𝑋 𝑘,𝑗2 ))) 𝐾 𝑢𝑎 (𝑞 (𝐵 ̅ (𝑋 𝑘,𝑗2 ))) 𝐾 𝑎𝑢 (𝑞 (𝐵 ̅ (𝑋 𝑘,𝑗2 ))) 𝐾 𝑎𝑎 (𝑣 𝑆 (𝐵 ̅ (𝑋 𝑘,𝑗2 ))) ] → 𝑎 𝑖𝑚 (𝑋 𝑘,𝑗1 ), 𝑎 𝑖𝑚 (𝑋 𝑘,𝑗2 )

									𝜕𝑎 𝑖𝑚 𝜕𝑥 𝑗	≈	𝑎 𝑖𝑚 (𝑥 𝑗 +∆)-𝑎 𝑖𝑚 (𝑥 𝑗 -∆) 2∆	(C-7)
															],
				[	𝐾 𝑢𝑢 ((C-8)
	𝜕𝐹 𝑖 𝜕𝑥 𝑗	= 𝑥 1	𝜕𝑎 𝑖1 𝜕𝑥 𝑗	+ 𝑥 2	𝜕𝑎 𝑖2 𝜕𝑥 𝑗	+ ⋯ + 𝑥 𝑗	𝜕𝑎 𝑖𝑗 𝜕𝑥 𝑗	+ 𝑎 𝑖𝑗 + 𝑥 𝑗+1	𝜕𝑎 𝑖,𝑗+1 𝜕𝑥 𝑗	+ ⋯ + 𝑥 𝑛	𝜕𝑎 𝑖𝑛 𝜕𝑥 𝑗	-	𝜕𝑏 𝑖 𝜕𝑥 𝑗	= 𝑎 𝑖𝑗 + ∑ 𝑛 𝑚=1 (𝑥 𝑚 •
	𝜕𝑎 𝑖𝑚 𝜕𝑥 𝑗	)												(C-6)
														𝜕𝑎 𝑖𝑚
														𝜕𝑥 𝑖

, we replace the partial derivative terms with the central difference quotients by using an infinitesimal ∆ when calculating:

  𝑢 2 ̇+ 𝑍 2 (𝑢 1 ̇-𝑢 2 ̇) + 𝜑 𝑚 ′𝐻 𝑎𝑐 + 𝜑 𝑝 ′𝑈 (D-5)

	= 𝑤 (-𝑚𝑡 𝑚 ′ (	𝑢 2 ̇-𝑢 1 ̇𝑐𝑜𝑠(𝑘𝐿) 𝑗𝑣̅ 𝑠𝑖𝑛(𝑘𝐿) • 𝑠 33 𝐻 -	𝑑 33 𝑚 𝑠 33 𝐻 𝐻 𝑎𝑐 ) -(𝑚 -1)𝑡 𝑝 ′ (	𝑢 2 ̇-𝑢 1 ̇𝑐𝑜𝑠(𝑘𝐿) 𝑗𝑣̅ 𝑠𝑖𝑛(𝑘𝐿) • 𝑠 11 𝐸 -	𝑝 𝑑 31 𝑠 11 𝐸 𝐸 2 ))
	= (	𝑚𝑡 𝑚 ′𝑤 𝑗𝑣̅ 𝑠𝑖𝑛(𝑘𝐿) • 𝑠 33 𝐻 + 𝑚𝑡 𝑚 ′𝑤𝑑 33 (𝑚 -1)𝑡 𝑝 ′𝑤 𝑗𝑣̅ 𝑠𝑖𝑛(𝑘𝐿) • 𝑠 11 𝐸 ) (𝑢 1 ̇-𝑢 2 ̇) + 𝑡𝑎𝑛 ( 𝑘𝐿 2 𝑚 𝑠 33 𝐻 𝐻 𝑎𝑐 + 𝑝 (𝑚 -1)𝑡 𝑝 ′𝑤𝑑 31 𝑠 11 𝐸 𝐸 2	) (	𝑚𝑡 𝑚 ′𝑤 𝑗𝑣̅ 𝑠 33 𝐻 +	(𝑚 -1)𝑡 𝑝 ′𝑤 𝐸 𝑗𝑣̅ 𝑠 11	) 𝑢 1 +
								= 𝑍 1 𝑢 1 ̇+ 𝑍 2 (𝑢 1 ̇-𝑢 2 ̇) + 𝜑 𝑚 ′𝐻 𝑎𝑐 + 𝜑 𝑝 ′𝑈	(D-4)
			𝐹(𝐿) = -𝑚𝑡 𝑚 ′𝑤𝑇 11 𝑚 (𝐿) -(𝑚 -1)𝑡 𝑝 ′𝑤𝑇 11 𝑝 (𝐿)
				= -𝑚𝑡 𝑚 ′𝑤 (	𝑆 11 𝑚 (𝐿) 𝑠 33 𝐻	-	𝑑 33 𝑚 𝑠 33 𝐻 𝐻 𝑎𝑐 ) -(𝑚 -1)𝑡 𝑝 ′𝑤 (	𝑆 11 𝑝 (𝐿) 𝑠 11 𝐸	-	𝑝 𝑑 31 𝑠 11 𝐸 𝐸 2 )
	= 𝑤 (-𝑚𝑡 𝑚 ′ (	𝑢 2 ̇𝑐𝑜𝑠(𝑘𝐿) -𝑢 1 j𝑣̅ 𝑠𝑖𝑛(𝑘𝐿) • 𝑠 33 𝐻 -𝑑 33 𝑚 𝑠 33 𝐻 𝐻 𝑎𝑐 ) -(𝑚 -1)𝑡 𝑝 ′ (	𝑢 2 ̇𝑐𝑜𝑠(𝑘𝐿) -𝑢 1 j𝑣̅ 𝑠𝑖𝑛(𝑘𝐿) • 𝑠 11 𝐸 -𝑝 𝑑 31 𝑠 11 𝐸 𝐸 2 ))
	= (	𝑚𝑡 𝑚 ′𝑤 𝑗𝑣̅ 𝑠𝑖𝑛(𝑘𝐿) • 𝑠 33 𝐻 + 𝑚𝑡 𝑚 ′𝑤𝑑 33 (𝑚 -1)𝑡 𝑝 ′𝑤 𝑗𝑣̅ 𝑠𝑖𝑛(𝑘𝐿) • 𝑠 11 𝐸 ) (𝑢 1 ̇-𝑢 2 ̇) -𝑡𝑎𝑛 ( 𝑘𝐿 2 𝑚 𝑠 33 𝐻 𝐻 𝑎𝑐 + 𝑝 (𝑚 -1)𝑡 𝑝 ′𝑤𝑑 31 𝑠 11 𝐸 𝐸 2	) (	𝑚𝑡 𝑚 ′𝑤 𝑗𝑣̅ 𝑠 33 𝐻 +	(𝑚 -1)𝑡 𝑝 ′𝑤 𝐸 𝑗𝑣̅ 𝑠 11	) 𝑢 2 +
	= -𝑍 1 where 𝑣̅ = √ 𝑛 𝑠 33 1-𝑛 𝐻 + 𝑠 11 𝐸 𝜌 ̅ , 𝑍 1 = 𝑗𝑤(𝑚𝑡 𝑚 ′ + (𝑚 -1)𝑡 𝑝 ′)𝜌̅ 𝑣̅ 𝑡𝑎𝑛 ( 𝑘𝐿 2	), 𝑍 2 =	𝑤(𝑚𝑡 𝑚 ′+(𝑚-1)𝑡 𝑝 ′)𝜌 ̅𝑣 ̅ 𝑗𝑠𝑖𝑛(𝑘𝐿)
				11 𝑚 (0) -(𝑚 -1)𝑡 𝑝 ′𝑤𝑇 11 𝑝 (0)
				= -𝑚𝑡 𝑚 ′𝑤 (	𝑆 11 𝑚 (0) 𝑠 33 𝐻	-	𝑑 33 𝑚 𝑠 33 𝐻 𝐻 𝑎𝑐 ) -(𝑚 -1)𝑡 𝑝 ′𝑤 (	𝑆 11 𝑝 (0) 𝑠 11 𝐸	-	𝑝 𝑑 31 𝑠 11 𝐸 𝐸 2 )

  In the parallel connectivity type, the output voltage 𝑈 = 𝑡 𝑝 ′𝐸 2 , thus 𝜑 𝑝 ′ =

				𝐼 𝑝 =	(𝑚-1)𝑑 𝑑𝑡	∫ 𝑤𝐷 2 𝑑𝑥 𝐿 0	= 𝑤(𝑚 -1)	𝑑 𝑑𝑡 ∫ (𝑑 31 𝑝 ( 𝑆 11 𝑝 (𝑥) 𝑠 11 𝐸 𝐿 0	-	𝑑 31 𝑝 𝑠 11 𝐸 𝐸 2 ) + 𝜀 33 𝑇 𝐸 2 ) 𝑑𝑥	= 𝑤(𝑚 -
	1)	𝑑 𝑑𝑡 ∫ (𝑑 31 𝑝 𝑆 11 𝑝 (𝑥) 𝑠 11 𝐸 𝐿 0	+ 𝜀 33 𝑆 𝐸 2 )𝑑𝑥	=	𝑤(𝑚-1)𝑑 31 𝑝 𝑠 11 𝐸	𝑑(𝑢(𝐿)-𝑢(0)) 𝑑𝑡	+ 𝑗𝜔 • 𝑤(𝑚 -1)𝐿𝜀 33 𝑆 𝐸 2 = -𝜑 𝑝 ′(𝑢 1 ̇-𝑢 2 ̇) +
										𝑗𝜔𝐶 𝑇 𝑈	(D-7)
	where 𝐶 𝑇 =	𝑤(𝑚-1)𝐿𝜀 33 𝑆 𝑡 𝑝 ′	, 𝜀 33 𝑆 = 𝜀 33 𝑇 -	𝑑 31 𝑝 2
			𝑑 𝑑𝑡 ∫ 𝑤𝐷 2 𝑑𝑥 𝐿 0	= 𝑤	𝑑 𝑑𝑡 ∫ (𝑑 31 𝑝 ( 𝑆 11 𝑝 (𝑥) 𝑠 11 𝐸 𝐿 0	-	𝑑 31 𝑝 𝑠 11 𝐸 𝐸 2 ) + 𝜀 33 𝑇 𝐸 2 )𝑑𝑥	= 𝑤	𝑑 𝑑𝑡 ∫ (𝑑 31 𝑝 𝑆 11 𝑝 (𝑥) 𝑠 11 𝐸 𝐿 0	+ 𝜀 33 𝑆 𝐸 2 )𝑑𝑥	=
	𝑤𝑑 31 𝑝 𝑠 11 𝐸	𝑑(𝑢(𝐿)-𝑢(0)) 𝑑𝑡	+ 𝑗𝜔 • 𝑤𝐿𝜀 33 𝑆 𝐸 2 = -𝜑 𝑝 ′(𝑢 1 ̇-𝑢 2 ̇) + 𝑗𝜔𝐶 𝑇 𝑈	(D-6)
	where 𝐶 𝑇 =	𝑤𝐿𝜀 33 𝑆 (𝑚-1)𝑡 𝑝 ′ , 𝜀 33 𝑆 = 𝜀 33 𝑇 -	𝑑 31 𝑝 2 𝑠 11 𝐸
										𝑝 𝑤(𝑚-1)𝑑 31 𝑠 11 𝐸	, and the coupling
	current through the piezoelectric layers can be calculated as [40]:

  𝐶 𝑇 ) )(𝑢 1 ̇-𝑢 2 ̇) + 𝜑 𝑚 ′𝐻 𝑎𝑐 + 𝜑 𝑝 ′𝑈 (D-10) 𝑗𝜔(-𝐶 𝑇 ) )(𝑢 1 ̇-𝑢 2 ̇) + 𝜑 𝑚 ′𝐻 𝑎𝑐 + 𝜑 𝑝 ′𝑈(D-11)

	= -	𝑤𝑚𝑡 𝑚 ′ 𝑠 33 𝐻	•	𝑢 2 ̇-𝑢 1 ̇𝑐𝑜𝑠(𝑘𝐿) 𝑗𝑣̅ 𝑠𝑖𝑛(𝑘𝐿)	+	𝑚𝑡 𝑚 ′ 𝑤𝑑 33 𝑚 𝑠 33 𝐻	𝐻 𝑎𝑐
					+	(𝑚 -1)𝑡 𝑝 ′ 𝑤𝛽 33 ̅̅̅̅ 𝑑 33 𝑝 𝑠 33 𝐸	(	1 ̅̅̅̅ (𝑚 -1)𝐿 𝛽 33	𝑈 -	𝑑 33 𝑝 𝑠 33 𝐸 𝐿	•	1 𝑗𝜔	(𝑢 1 ̇-𝑢 2 ̇)) -	(𝑚 -1)𝑡 𝑝 ′ 𝑤 𝑠 33 𝑑 33 𝑝 2 𝐸 -𝜀 33 𝑇
					•	𝑢 2 ̇-𝑢 1 ̇𝑐𝑜𝑠(𝑘𝐿) 𝑗𝑣̅ 𝑠𝑖𝑛(𝑘𝐿)
										𝜑 𝑝 ′ 2 𝑗𝜔(-𝐹(𝐿) = -𝑚𝑡 𝑚 = 𝑍 1 𝑢 1 ̇+ (𝑍 2 + ′ 𝑤𝑇 11 𝑚 (𝐿) -(𝑚 -1)𝑡 𝑝 ′ 𝑤𝑇 11 𝑝 (𝐿)
									= -𝑚𝑡 𝑚 ′ 𝑤 (	𝑆 11 𝑚 (𝐿) 𝑠 33 𝐻	-	𝑑 33 𝑚 𝑠 33 𝐻 𝐻 𝑎𝑐 ) -(𝑚 -1)𝑡 𝑝 ′ 𝑤 (	𝑆 11 𝑝 (𝐿) 𝑠 33 𝐸	-	𝑝 𝑑 33 𝑠 33 𝐸 𝐸 1 )
	= 𝑤 (-𝑚𝑡 𝑚 ′ (	𝑆 11 𝑚 (𝐿) 𝑠 33 𝐻	-	𝑑 33 𝑚 𝑠 33 𝐻 𝐻 𝑎𝑐 ) -(𝑚 -1)𝑡 𝑝 (-𝛽 33 ̅̅̅̅	𝑑 33 𝑝 𝑠 33 𝐸 𝐷 1 + ( 𝑠 33 1 𝐸 + 𝛽 33 ̅̅̅̅ (	𝑝 𝑑 33 𝑠 33 𝐸 )	2	) 𝑆 11 𝑝 (𝐿)))
	= -	𝑤𝑚𝑡 𝑚 ′ 𝑠 33 𝐻	•	𝑢 2 ̇𝑐𝑜𝑠(𝑘𝐿) -𝑢 1 j𝑣̅ sin(𝑘𝐿) • 𝑠 33 𝐻 + 𝑚𝑡 𝑚 𝑚 ′ 𝑤𝑑 33 𝑠 33 𝐻	𝐻 𝑎𝑐
					+	(𝑚 -1)𝑡 𝑝 ′ 𝑤𝛽 33 ̅̅̅̅ 𝑑 33 𝑝 𝑠 33 𝐸	(	1 ̅̅̅̅ (𝑚 -1)𝐿 𝛽 33	𝑈 -	𝑑 33 𝑝 𝑠 33 𝐸 𝐿	•	1 𝑗𝜔	(𝑢 1 ̇-𝑢 2 ̇)) -	(𝑚 -1)𝑡 𝑝 ′ 𝑤 𝑠 33 𝑑 33 𝑝 2 𝐸 -𝜀 33 𝑇
					•	𝑢 2 ̇𝑐𝑜𝑠(𝑘𝐿) -𝑢 1 j𝑣̅ sin(𝑘𝐿) • 𝑠 33 𝐻
											= -𝑍 1 𝑢 2 ̇+ (𝑍 2 +	𝜑 𝑝 ′ 2
	where 𝑣̅ = √ 𝑠 33 𝑛 𝐻 + 𝜌 ̅ 1-𝑛 𝑠 33 𝐸 𝜑 𝑚 ′ = 𝑚𝑡 𝑚 ′ 𝑤𝑑 33 𝑚 𝑠 33 𝐻 , 𝜑 𝑝 ′ =	𝑚𝑡 𝑚 ′ 𝑤𝑇 11 𝑚 (0) -(𝑚 -1)𝑡 𝑝 ′ 𝑤𝑇 11 𝑝 (0) ′ 𝑤 ( = -𝑚𝑡 𝑚 𝑆 11 𝑚 (0) 𝑠 33 𝐻 -𝑑 33 𝑚 𝑠 33 𝐻 𝐻 𝑎𝑐 ) -(𝑚 -1)𝑡 𝑝 ′ 𝑤 ( , 𝑍 1 = 𝑗𝑤(𝑚𝑡 𝑚 ′ + (𝑚 -1)𝑡 𝑝 ′ )𝜌̅ 𝑣̅ 𝑡𝑎𝑛 ( 𝑘𝐿 2 ) , 𝑍 2 = 𝑆 11 𝑝 (0) 𝑠 33 𝐸 𝑤(𝑚𝑡 𝑚 -𝑝 ′ +(𝑚-1)𝑡 𝑝 ′ )𝜌 ̅𝑣 ̅ 𝑗𝑠𝑖𝑛(𝑘𝐿) 𝑑 33 𝑠 33 𝐸 𝐸 1 ) 𝑡 𝑝 ′ 𝑤 𝐿 𝑑 33 𝑝 𝑠 33 ′ 𝑤𝑡 𝑝 𝐸 , 𝐶 𝑇 = (𝑚-1)𝐿𝛽 33 ̅̅̅̅̅	,
	= 𝑤 (-𝑚𝑡 𝑚 ′ (	𝑆 11 𝑚 (0) 𝑠 33 𝐻	-	𝑑 33 𝑚 𝑠 33 𝐻 𝐻 𝑎𝑐 ) -(𝑚 -1)𝑡 𝑝 ′ (-𝛽 33 ̅̅̅̅	𝑑 33 𝑝 𝑠 33 𝐸 𝐷 1 + ( 𝑠 33 1 𝐸 + 𝛽 33 ̅̅̅̅ (	𝑝 𝑑 33 𝑠 33 𝐸 )	2	) 𝑆 11 𝑝 (0)))

  𝑗𝜔(-𝐶 𝑇 ) )(𝑢 1 ̇-𝑢 2 ̇) + 𝜑 𝑚 ′𝐻 𝑎𝑐 + 𝜑 𝑝 ′𝑈 (D-15)𝑗𝜔(-𝐶 𝑇 ) )(𝑢 1 ̇-𝑢 2 ̇) + 𝜑 𝑚 ′𝐻 𝑎𝑐 + 𝜑 𝑝 ′𝑈 (D-16)

											𝐹(𝐿) = -𝑚𝑡 𝑚 ′ 𝑤𝑇 11 𝑚 (𝐿) -(𝑚 -1)𝑡 𝑝 ′ 𝑤𝑇 11 𝑝 (𝐿)
				= -𝑚𝑡 𝑚 ′ 𝑤 (	𝑆 11 𝑚 (𝐿) 𝑠 33 𝐻	-	𝑑 33 𝑚 𝑠 33 𝐻 𝐻 𝑎𝑐 ) -(𝑚 -1)𝑡 𝑝 ′ 𝑤 (	𝑆 11 𝑝 (𝐿) 𝑠 33 𝐸	-	𝑝 𝑑 33 𝑠 33 𝐸 𝐸 1 )
	= 𝑤 (-𝑚𝑡 𝑚 ′ (	𝑆 11 𝑚 (𝐿) 𝑠 33 𝐻	-	𝑑 33 𝑚 𝑠 33 𝐻 𝐻 𝑎𝑐 ) -(𝑚 -1)𝑡 𝑝 ′ (-𝛽 33 ̅̅̅̅	𝑑 33 𝑝 𝑠 33 𝐸 𝐷 1 + ( 𝑠 33 1 𝐸 + 𝛽 33 ̅̅̅̅ (	𝑝 𝑑 33 𝑠 33 𝐸 )	2	) 𝑆 11 𝑝 (𝐿)))
	= -	𝑤𝑚𝑡 𝑚 ′ 𝑠 33 𝐻	•	𝑢 2 ̇𝑐𝑜𝑠(𝑘𝐿) -𝑢 1 j𝑣̅ 𝑠𝑖𝑛(𝑘𝐿) • 𝑠 33 𝐻 + 𝑚𝑡 𝑚 𝑚 ′ 𝑤𝑑 33 𝑠 33 𝐻	𝐻 𝑎𝑐
								+	(𝑚 -1)𝑡 𝑝 ′ 𝑤𝛽 33 ̅̅̅̅ 𝑑 33 𝑝 𝑠 33 𝐸	(	1 𝛽 33 ̅̅̅̅ 𝐿	𝑈 -	𝑑 33 𝑝 𝑠 33 𝐸 𝐿	•	1 𝑗𝜔	(𝑢 1 ̇-𝑢 2 ̇)) -	(𝑚 -1)𝑡 𝑝 ′ 𝑤 𝑠 33 𝑑 33 𝑝 2 𝐸 -𝜀 33 𝑇
								•	𝑢 2 ̇𝑐𝑜𝑠(𝑘𝐿) -𝑢 1 j𝑣̅ 𝑠𝑖𝑛(𝑘𝐿) • 𝑠 33 𝐻
	= -𝑍 1 𝑢 2 ̇+ (𝑍 2 +	𝜑 𝑝 ′ 2
	where 𝑣̅ = √ 𝑠 33 𝑛 𝐻 + 𝜌 ̅ 1-𝑛 𝑠 33 𝐸	, 𝑍 1 = 𝑗𝑤(𝑚𝑡 𝑚 ′ + (𝑚 -1)𝑡 𝑝 ′ )𝜌̅ 𝑣̅ 𝑡𝑎𝑛 ( 𝑘𝐿 2 ) , 𝑍 2 =	𝑤(𝑚𝑡 𝑚 ′ +(𝑚-1)𝑡 𝑝 ′ )𝜌 ̅𝑣 ̅ 𝑗𝑠𝑖𝑛(𝑘𝐿)	, 𝜑 𝑚 ′ =
	𝑚𝑡 𝑚 ′ 𝑤𝑑 33 𝑚 𝑠 33 𝐻	, 𝜑 𝑝 ′ =	(𝑚-1)𝑡 𝑝 ′ 𝑤 𝐿	𝑑 33 𝑝 𝑠 33 𝐸 , 𝐶 𝑇 =	𝑤𝑡 𝑝 ′ (𝑚-1) 𝐿𝛽 33 ̅̅̅̅̅
											𝑚𝑡 𝑚 ′ 𝑤𝑇 11 𝑚 (0) -(𝑚 -1)𝑡 𝑝 ′ 𝑤𝑇 11 𝑝 (0)
				= -𝑚𝑡 𝑚 ′ 𝑤 (	𝑆 11 𝑚 (0) 𝑠 33 𝐻	-	𝑑 33 𝑚 𝑠 33 𝐻 𝐻 𝑎𝑐 ) -(𝑚 -1)𝑡 𝑝 ′ 𝑤 (	𝑆 11 𝑝 (0) 𝑠 33 𝐸	-	𝑝 𝑑 33 𝑠 33 𝐸 𝐸 1 )
	= 𝑤 (-𝑚𝑡 𝑚 ′ (	𝑆 11 𝑚 (0) 𝑠 33 𝐻	-	𝑑 33 𝑚 𝑠 33 𝐻 𝐻 𝑎𝑐 ) -(𝑚 -1)𝑡 𝑝 ′ (-𝛽 33 ̅̅̅̅	𝑑 33 𝑝 𝑠 33 𝐸 𝐷 1 + ( 𝑠 33 1 𝐸 + 𝛽 33 ̅̅̅̅ (	𝑝 𝑑 33 𝑠 33 𝐸 )	2	) 𝑆 11 𝑝 (0)))
	= -	𝑤𝑚𝑡 𝑚 ′ 𝑠 33 𝐻	•	𝑢 2 ̇-𝑢 1 ̇𝑐𝑜𝑠(𝑘𝐿) 𝑗𝑣̅ 𝑠𝑖𝑛(𝑘𝐿)	+	𝑚𝑡 𝑚 ′ 𝑤𝑑 33 𝑚 𝑠 33 𝐻	𝐻 𝑎𝑐
								+	(𝑚 -1)𝑡 𝑝 ′ 𝑤𝛽 33 ̅̅̅̅ 𝑑 33 𝑝 𝑠 33 𝐸	(	1 𝛽 33 ̅̅̅̅ 𝐿	𝑈 -	𝑑 33 𝑝 𝑠 33 𝐸 𝐿	•	1 𝑗𝜔	(𝑢 1 ̇-𝑢 2 ̇)) -	(𝑚 -1)𝑡 𝑝 ′ 𝑤 𝑠 33 𝑑 33 𝑝 2 𝐸 -𝜀 33 𝑇
								•	𝑢 2 ̇-𝑢 1 ̇𝑐𝑜𝑠(𝑘𝐿) 𝑗𝑣̅ 𝑠𝑖𝑛(𝑘𝐿)
	= 𝑍 1 𝑢 1 ̇+ (𝑍 2 +		𝜑 𝑝 ′ 2

Relative permittivity: 𝜀 11 = 𝜀 22 = 3130, 𝜀 33 = 3400

Terfenol-D [91]:

Relative permeability: 𝜇 11 = 𝜇 22 = -469.5, 𝜇 33 = 124.9

𝐁𝐚𝐓𝐢𝐎 𝟑 [START_REF] Ramirez | [END_REF]:

Elastic compliances (pm [START_REF] Carman | Nonlinear constitutive relations for magnetostrictive materials with applications to 1-D problems[END_REF] by means of truncating the Taylor series expansion of the Gibbs free energy function after some high order terms. Briefly presented, the formulations are derived based on the energy balance equation ignoring the heat effect terms among an electric-magnetic body

where 𝑓 𝑖 , 𝑡 𝑖 , 𝑢 𝑖 are the components of body force, surface force and the displacement vector, respectively, the dots above represents differentiation with respect to time, 𝜌 the mass density, 𝑈 the internal energy density per unit mass, Φ = -∇ • (𝑬 × 𝑯) the electric-magnetic energy density per unit time where 𝑬 is the electric field intensity vector and 𝑯 the magnetic field intensity vector.

By simplify the balance equation using the mass conservation and momentum conservation equations under quasi-magnetostatic and small distortion assumptions, one can obtain

where 𝜀 𝑖𝑗 and 𝜎 𝑖𝑗 are the components of the strain tensor and stress tensor, respectively.

Then, based on the expression of the Gibbs free energy without the heat effect term 𝐺 = 𝑈 -𝐵 𝑘 𝐻 𝑘 and by taking the magnetic field and stress as the independent variables, one can derive the total derivative of Gibbs free energy as 𝑑𝐺 = -𝜀 𝑖𝑗 𝑑𝜎 𝑖𝑗 -𝐵 𝑘 𝑑𝐻 𝑘 (B-3) and thus the Taylor series expansion of Gibbs free energy function as

Appendix C. Modified Newton-Raphson method

In the nonlinear case, since the coefficients 𝑣 𝑆 , 𝑞 and c 𝐵 = c 𝐻 + 𝑞 𝑡 𝑣 𝑆 𝑞 are all functions of the 𝑩 and as well as those of the unknown variable 𝒂 according to equations (2-4-57) to (2-4-59) and (2-4-81), the stiffness matrix components 𝐾 𝑢𝑢 , 𝐾 𝑢𝑎 , 𝐾 𝑎𝑢 and 𝐾 𝑎𝑎 in the equation (2-4-28) are thus not constant. Therefore, in order to solve the problem equations system an iteration method must be implemented.

In this paper's study, we choose the Newton-Raphson method with some modifications for the convenience of calculating the tangent stiffness matrix in the problem resolution. Some definitions of the magnetostrictive phase part of the equations system are made as follows before the method is employed: