A Urdaneta-Calzadilla

N Galopin

I Niyonzima

O Chadebec

B Bannwarth

G Meunier

A FEM-BEM coupling strategy for the modeling of magnetoelectric effects in composite structures

Keywords: FEM-BEM coupling, electro-magneto-mechanical, magnetoelectric composite, block Gauss-Seidel solver

This paper deals with the numerical modeling of devices based on magnetoelectric composite materials. These heterogeneous structures made of the mechanical association of piezoelectric and magnetostrictive materials display magneto-electric effects exceeding by several orders of magnitude the response of single-phase multiferroic materials. A coupling of the Finite Element Method (FEM) and the Boundary Element Method (BEM) is used to model the behavior of magnetic effects, while classical FEM formulations are used for the electrical and mechanical problems. This coupling of numerical methods allows avoiding considering a free space domain around the active domain, and thus to use a single mesh for the magnetic, mechanical and electrical problems. This results in a consequent reduction of the number of unknowns, which is accompanied by shorter computation times compared to a pure FEM approach. The final system of equations is solved by a block Gauss-Seidel type solver.

Introduction

Energy conversion in electrical transducers or actuators is based on electromagnetic interactions, which link the electromotive force to temporal variations of the magnetic flux density. However, these phenomena are sometimes difficult to exploit, especially for small devices submitted to very low frequency fields. The use of active composite structures, in particular, magnetoelectric (ME) composite structure can help address this issue [START_REF] Eerenstein | Multiferroic and magnetoelectric materials[END_REF].

The ME coupling consists in the existence of an electric polarization induced by a magnetization or, conversely, of a magnetization induced by an electric polarization. Materials with such properties have opened up possibilities of new applications in various fields, such as magnetic field sensors [START_REF] Liang | Roadmap on magnetoelectric materials and devices[END_REF], tunable radio-frequency magnetic filters [START_REF] Cheng | Recent development and status of magnetoelectric materials and devices[END_REF], antennas [START_REF] Dong | A portable very low frequency (vlf) communication system based on acoustically actuated magnetoelectric antennas[END_REF], gyrators [START_REF] Leung | Highly efficient solid state magnetoelectric gyrators[END_REF], energy harvesters [START_REF] Lafont | Magnetostrictive-piezoelectric composite structures for energy harvesting[END_REF], memory devices based on the principle of electric writingmagnetic reading [START_REF] Nan | Multiferroic magnetoelectric composites: Historical perspective, status, and future directions[END_REF][START_REF] Scott | Multiferroic memories[END_REF], biology and medicine [START_REF] Kopyl | Magnetoelectric effect: principles and applications in biology and medicine-a review[END_REF]. This interaction emerges as a material property and does not follow directly from Maxwell's equations. Past research has been conducted in order to obtain single-phase materials that simultaneously exhibit coupled magnetic and electrical ferroic orders, also known as multiferroics. Unfortunately, despite many efforts and with few exceptions [START_REF] Mundy | Atomically engineered ferroic layers yield a room-temperature magnetoelectric multiferroic[END_REF], the majority of single-phase multiferroic materials are made of complex oxides that exhibit very low coefficients and mainly at low temperatures [START_REF] Astrov | Magnetoelectric effect in chromium oxide[END_REF][START_REF] Vaz | Magnetoelectric coupling effects in multiferroic complex oxide composite structures[END_REF]. It has been shown that ordinary magnetic and electrical susceptibilities provide an upper bound on the coefficient for single-phase magnetoelectric materials [START_REF] Brown | Upper bound on the magnetoelectric susceptibility[END_REF]. This disadvantage has been successfully overcome by the fabrication of magnetoelectric composites which consist of coupled magneto-mechanical and electro-mechanical phases [START_REF] Eerenstein | Giant sharp and persistent converse magnetoelectric effects in multiferroic epitaxial heterostructures[END_REF][START_REF] Priya | Recent advancements in magnetoelectric particulate and laminate composites[END_REF], the resulting heterogeneous materials show large ME coefficients. The aim of such composites is to generate the intended magnetoelectric effect as a deformationinduced product property [START_REF] Van Suchtelen | Product properties: a new application of composite materials[END_REF], a property that their individual constituents do not have. In the direct ME effect, an applied magnetic field causes a deformation of the magneto-mechanical coupled phase which is transmitted to the electro-mechanical coupled phase. As a result, a strain-induced polarization modulation in the electric phase is obtained. While in the inverse effect, an applied electric field causes a deformation of the electro-mechanical phase, which is then transmitted to the magneto-mechanical phase (reverse effect). This results in strain-induced magnetization modulation.

The behavior, performance and effective properties of magnetoelectric composites depend on the material composition of each phase, the bonds between the different phases, their morphology and the electrical resistance of the composite. The electro-mechanical coupled phases can be composed of ferroelectric (FE) materials (BaTiO 3 , PbTiO 3 , Pb(Zr,Ti)TiO 3 ) or piezoelectric materials (PE) ((1-x)Pb(Mg1/3Nb2/3)O 3-x PbTiO 3 (PMN-PT), PZT), while the magneto-mechanical coupled phases can be made of ferrimagnetic (FI) (CoFe 2 O 4 , NiFe 2 O 4 ) or ferromagnetic (FM) materials (CoFeB, FeGa). These heterostructures can be either multiferroic (and also magnetoelectric), comprising FE/FM materials, or magnetoelectric but not multiferroic, i.e., having only one ferroic order, often with PE/FM materials. Particulate composites formed from a microscale mixture of FE and FI materials show magnetoelectric coefficient values between 1 and 500 mV cm -1 Oe -1 at low frequencies [START_REF] Islam | Effect of piezoelectric grain size on magnetoelectric coefficient of pb(zr0.52ti0.48)o3-ni0.8zn0.2fe2o4 particulate composites[END_REF], while larger values are obtained at mechanical resonance [START_REF] Zeng | Resonance magnetoelectric effect in bulk composites of lead zirconate titanate and nickel ferrite[END_REF]. In comparison, laminated composites composed of the same ferroelectric and ferrimagnetic materials show magnetoelectric coefficient values that are an order of magnitude higher. The highest magnetoelectric coupling coefficients (> 5 V Cm -1 Oe -1 ) are obtained by layered heterogeneous structures whose magneto-mechanical phase is composed of either a giant magnetostriction alloy (Terfenol-D: Tb 1-x Dy x Fe 2 ) [START_REF] Ryu | Effect of the magnetostrictive layer on magnetoelectric properties in lead zirconate titanate/terfenol-d laminate composites[END_REF] or an amorphous alloy with high magnetic permeability and piezomagnetic coefficients (Metglass: Fe-based alloys) [START_REF] Li | Enhancement in magnetic field sensitivity and reduction in equivalent magnetic noise by magnetoelectric laminate stacks[END_REF]. For the manufacture of laminated composites based on ceramics and alloys, co-sintering and bonding are the most popular methods.

The complex behavior of the magnetoelectric composite materials described above involves not only the definition of appropriate material models, but also the formulation and solution of fully coupled boundary-value problems, especially for the development of technological applications. Several approaches including analytical methods, semi-analytical methods and numerical methods have been applied to the prediction of the overall material properties and to the investigation of the coupling behavior of the magnetoelectric materials.

Analytical approaches based on Green's functions have been proposed, for example by Nan et al. [START_REF] Nan | Magnetoelectric effect in composites of piezoelectric and piezomagnetic phases[END_REF][START_REF] Nan | Magnetoelectricity of multiferroic composites[END_REF][START_REF] Nan | Magneticfield-induced electric polarization in multiferroic nanostructures[END_REF], Pan [START_REF] Pan | Three-dimensional green's functions in anisotropic magneto-electro-elastic bimaterials[END_REF], Wang and Shen [START_REF] Wang | The general solution of threedimensional problems in magnetoelectroelastic media[END_REF]. Elastodynamics methods combining the equation of motion of continuous media with mechanical and electrical boundary conditions have also been proposed, for example by Harshé et al. [START_REF] Harshe | Theoretical modelling of multilayer magnetoelectric composites[END_REF], Avellaneda and Harshé [START_REF] Avellaneda | Magnetoelectric effect in piezoelectric/magnetostrictive multilayer (2-2) composites[END_REF], Wu et al. [START_REF] Wu | Inequivalence of direct and converse magnetoelectric coupling at electromechanical resonance[END_REF], Muchenik and Barbero [START_REF] Muchenik | Charge, voltage, and workconversion formulas for magnetoelectric laminated composites[END_REF]. In the previous approach, the response of ME materials was studied assuming linear behaviors for the ferroelectric and ferromagnetic phases. These models have shown how the volume fractions of each phase, the connectivity [START_REF] Newnham | Connectivity and piezoelectric-pyroelectric composites[END_REF], as well as the piezoelectric, piezomagnetic and elastic properties participate in the ME coupling. The nonlinear response and stress dependence of the ME composite are addressed, for ex-ample, by Wang and Zhou [START_REF] Wang | A theoretical study of nonlinear magnetoelectric effect in magnetostrictive-piezoelectric trilayer[END_REF], Lin et al. [START_REF] Lin | An analytical nonlinear model for laminate multiferroic composites reproducing the dc magnetic bias dependent magnetoelectric properties[END_REF], Burdin et al. [START_REF] Burdin | Nonlinear magneto-electric effects in ferromagnetic-piezoelectric composites[END_REF], Shi et al. [START_REF] Shi | An analytical model for nonlinear magnetoelectric effect in laminated composites[END_REF] who considered the nonlinear anhysteretic behavior and the effect of stress in the constitutive relationships of the ferromagnetic phase. It results that pre-stress and a bias magnetic field improve the ME coefficients and the frequency-multiplying behavior of laminated composites. The interfaces between ferroelectric and ferromagnetic phases, which are not perfect and usually correspond to a layer of epoxy glue, are accounted for by means of an interface coupling factor [START_REF] Bichurin | Theory of low-frequency magnetoelectric coupling in magnetostrictivepiezoelectric bilayers[END_REF][START_REF] Wang | A theoretical study of nonlinear magnetoelectric effect in magnetostrictive-piezoelectric trilayer[END_REF]. Taking these interface effects into consideration allows for better prediction of the coupling behavior of ME composites.

Approaches based on the method of equivalent electrical circuits have also been proposed [START_REF] Dong | Longitudinal and transverse magnetoelectric voltage coefficients of magnetostrictive/piezoelectric laminate composite: theory[END_REF][START_REF] Dong | Equivalent circuit method for static and dynamic analysis of magnetoelectric laminated composites[END_REF]. They use an extension of Mason's model which allows establishing an equivalent electrical circuit of the ME composites, whose equivalent electrical parameters are established via a coupled equation of motion and the integration of the piezoelectric and magnetostrictive behavior laws. More specifically, these methods can be used to investigate the ME coefficient for different modes static and dynamic condititons and the electrical resistance load effect on the resonant ME coupling.

To describe and predict the behavior and effective properties of ME composites considering the microstructures and anisotropies of the materials, various micromechanical analyses were developed. Such analytical or semianalytical solutions are for example based on the homogenization and selfconsistent models [START_REF] Benveniste | Magnetoelectric effect in fibrous composites with piezoelectric and piezomagnetic phases[END_REF][START_REF] Kim | Micromechanical analysis of effective properties of magneto-electro-thermo-elastic multilayer composites[END_REF][START_REF] Jin | Macroscopic behavior prediction of multiferroic composites[END_REF][START_REF] Xu | A multi-scale and multi-field coupling nonlinear constitutive theory for the layered magnetoelectric composites[END_REF][START_REF] Zhang | Direct and converse nonlinear magnetoelectric coupling in multiferroic composites with ferromagnetic and ferroelectric phases[END_REF], Mori-Tanaka mean field theory [START_REF] Li | Micromechanics of magnetoelectroelastic composite materials: Average fields and effective behavior[END_REF][START_REF] Srinivas | The effective magnetoelectric coefficients of polycrystalline multiferroic composites[END_REF], variational asymptotic approach [START_REF] Tang | Variational asymptotic homogenization of heterogeneous electromagnetoelastic materials[END_REF][START_REF] Tang | Micromechanical modeling of the multiphysical behavior of smart materials using the variational asymptotic method[END_REF] or Eshelby's equivalent inclusion approach [START_REF] Huang | Analytical predictions for the magnetoelectric coupling in piezomagnetic materials reinforced by piezoelectric ellipsoidal inclusions[END_REF][START_REF] Li | Anisotropic coupled-field inclusion and inhomogeneity problems[END_REF].

All the previously described methods make it possible to study ME composites with trivial geometries such as composites with perfect ellipsoidal inclusions, perfect laminated structures or composites involving simple boundary value problems. Despite being computationally expensive, numerical tools are not restricted to specific topologies. FEM simulations were conducted for example, by Buchanan [START_REF] Buchanan | Layered versus multiphase magneto-electro-elastic composites[END_REF] and Galopin et al. [START_REF] Galopin | Finite element modeling of magnetoelectric sensors[END_REF] to study the multilayer and multiphase ME composites response, and by Lee et al. [START_REF] Lee | Effective properties of three-phase electro-magneto-elastic composites[END_REF], Avakian et al. [START_REF] Avakian | Nonlinear modeling and finite element simulation of magnetoelectric coupling and residual stress in multiferroic composites[END_REF] and Zhang et al. [START_REF] Zhang | Finite element analysis of the magnetoelectric effect on hybrid magnetoelectric composites[END_REF] in order to determine the effective properties of linear and non-linear multiphase ME composites. Nevertheless, a classical FE approach can become unsuitable for modeling ME devices with several disadvantages. One disadvantage is related to the stray magnetic field from the sources and ME composites. For simple geometries, at the boundary of the device, this stray field is linked to the geometry and proportional to the magnetization level. They can be accounted for by demagnetizing coefficients. This is no longer possible for more complex geometries or source configurations, and the treatment of an open boundary problem is required. Due to the decrease of the field away from the sources, an intuitive approach is to consider a sufficiently large but finite free space domain to respect the null field condition at infinity. Methods for treating open boundary problems have been developed including those involving infinite elements [START_REF] Chen | A review of finite element open boundary techniques for static and quasi-static electromagnetic field problems[END_REF][START_REF] Brunotte | Finite element modeling of unbounded problems using transformations: a rigorous, powerful and easy solution[END_REF] but they still require the presence of a free space domain which may lead to a huge mesh.

A second disadvantage of FEM is related to specific applications where the sources of the magnetic field, such as inductors, are located far from the active materials. An even larger free space domain containing both the field sources and the ME composite is needed. The free space domain can then be much larger than that associated with active materials. On the other hand, for problems with a huge free space domain compared to the active structure, the FEM leads to problems of accuracy and convergence [START_REF] Rakotoarison | Formal sensitivity computation of magnetic moment method[END_REF].

These problems can be avoided by coupling FEM with the BEM. This coupling of numerical methods is particularly well adapted to the numerical resolution of open boundary electromagnetic problems. For linear problems involving homogeneous materials embedded in the free space, only the boundary of the material domain needs to be discretized. The FEM-BEM approach also accounts for the nonlinear behavior of materials [START_REF] Meunier | Hybrid finite element boundary element solutions for three dimensional scalar potential problems[END_REF]. This paper proposes a FEM-BEM modeling of the electro-magneto-mechanically problem involving ME composites. This approach is particularly interesting since only the active domains associated with the materials need to be discretized and the free space domain is properly accounted for. In addition, a single mesh can be used for the coupled problem.

The outline of the paper is as follows. In section 2, after introducing the full set of continuous governing equations and constitutive laws, we develop weak and discrete forms of the coupled problem. In Section 3, we present the implemented iterative algorithm used for solving the multiphysics problem. In Section 4, we will apply the presented approach to the modeling of a three layer laminated PZT-5A/Terfenol-D energy harvester.

Theoretical framework of magnetoelectric coupling

In this section, the FEM-BEM coupling used to model the open boundary magnetoelectric problem is described. First, the multiphysics problem is introduced and the electro-mechanical and magneto-mechanical constitutive laws are defined. The weak and discrete formulations as well as the coupling between the FEM and the BEM are finally detailed. For the ease of reading, the basic magnetic, electric and mechanical quantities are summarized in Table 1.

V Q s Electric surface charge density C m -2 E Electric field V m -1 D Electric displacement C m -2 ϕ red Magnetic reduced scalar potential A J s Electric current density A m -2 H Magnetic field A m -1 B
Magnetic flux density T

Electro-magneto-mechanical problem description

We consider a domain Ω m = Ω pe ∪Ω pm and Ω 0 the domain exterior to Ω m , with Ω 0 having the properties of vacuum. The bond at the interface between the magnetostrictive and piezoelectric bodies is assumed to be perfect.

Electric and magnetic fields are computed from Maxwell's equation in static conditions,

∇ • B = 0 ∀x ∈ Ω, (1)
∇ × H = J s ∀x ∈ Ω, (2)
∇ • D = 0 ∀x ∈ Ω, (3)
∇ × E = 0 ∀x ∈ Ω, (4)
where the electric current density J s is taken as the input of the problem and considered null inside the active material and ∇• and ∇× indicate the divergence and curl operators respectively. Electric and magnetic equations are complemented by the conservation equation of linear momentum also in static conditions:

∇ • T + f = 0 ∀x ∈ Ω m , ( 5
)
where the body force f is also taken equal to zero inside the active material.

The boundary ∂Ω m of the mechanical domain Ω m = Ω pe ∪ Ω pm is partitioned as ∂Ω m = ∂Ω u m ∪ ∂Ω t m with ∂Ω u m ∩ ∂Ω t m = ∅, and the conditions for the displacements u s and the surface traction f T are given by:

u = u s ∀x ∈ ∂Ω u m , (6)
T • n = f T ∀x ∈ ∂Ω t m , (7)
where n denotes the outward unit normal vector to the surface ∂Ω. Similarly, the boundary ∂Ω e of the electrical domain Ω e is partitioned such that ∂Ω e = ∂Ω φ e ∪ ∂Ω d e and ∂Ω φ e ∩ ∂Ω d e = ∅, with the conditions for the electric potential φ 0 and the surface charge density Q s :

φ = φ 0 ∀x ∈ ∂Ω φ e , ( 8
)
D • n = -Q s ∀x ∈ ∂Ω d e . ( 9
)
In a pure FEM approach, the magnetic domain Ω contains the active materials and the free space domain Ω 0 , now of finite size, i.e., Ω = Ω m ∪ Ω 0 and Ω m = Ω pe ∪ Ω pm the magnetic subdomain associated with the active materials. Its boundary ∂Ω 0 would be the boundary of the air region partitioned such that ∂Ω 0 = ∂Ω h 0 ∪ ∂Ω b 0 and ∂Ω h 0 ∪ ∂Ω b 0 = ∅, with the conditions for the magnetic flux density Φ 0 and the magnetic field:

B • n = -Φ 0 ∀x ∈ ∂Ω b 0 (10) H × n = 0 ∀x ∈ ∂Ω h 0 (11)
In all generality, mixed boundary conditions can be applied on ∂Ω 0 . How- ever, in practice, B • n = 0 is often imposed in ∂Ω 0 . In FEM, Ω 0 is supposed to approach the open infinite free space domain surrounding the physical device so for the FEM approach to be accurate for the magnetic problem, a large enough domain Ω 0 must be considered. Pure FEM formulations of the magnetoelectric effect can be found in [58] [59]. The FEM-BEM coupling does not need such approximations, and the decay of magnetic field is properly taken into account without the need of explicitly considering a finite free space domain.

∂Ω d e -Q s n ∂Ω v e φ 0 Ω pe Ω pm (a) Electrical domain Ω pm ∂Ω t m ∂Ω u m f T u s Ω pe (b) Mechanical domain ∂Ω h 0 ∂Ω b 0 Ω 0 ϕ 0 J s ∂Ω ϕ mg ϕ 0 Ω pe Ω pm (c) Magnetic domain

Constitutive laws

Linear behavior laws are considered for both the electro-mechanical and the magneto-mechanical phases. Although the magneto-mechanical behavior is strongly non-linear, it is possible to describe it using a linearized piezomagnetic behavior law [START_REF]Handbook of Giant Magnetostrictive Materials[END_REF] obtained by considering a magneto-mechanical polarization state composed by a polarizing magnetic field and a mechanical prestress and by working under conditions which respects small variations around this state. An analogous approach was applied to the linear piezoelectric behavior in [START_REF][END_REF].

Assuming the existence of a thermodynamic potential, it follows that the constitutive relations depend on the choice of state variables. In this context, despite the heterogeneous nature of the considered structures, both the piezoelectric and piezomagnetic relations can be derived from the potential [START_REF] Vaz | Magnetoelectric coupling effects in multiferroic complex oxide composite structures[END_REF],

H me (E, H, S) = 1 2 S : C E,H : S H mech - 1 2 E • ε S • E H diel - 1 2 H • µ S • H H magn -E • e : S H piel -H • q : S H pimg , ( 12
)
where ":" denotes the double dot product. The different terms of ( 12) are the purely mechanical (H mech ), the purely dielectric (H diel ), the purely magnetic (H magn ), the piezoelectric (H piel ) and the piezomagnetic (H pimg ) energies.

In [START_REF] Vaz | Magnetoelectric coupling effects in multiferroic complex oxide composite structures[END_REF], C, ε, µ, e and q represent the fourth-order elasticity tensor, the second-order permittivity tensor, the second-order permeability tensor, the third-order piezoelectric tensor and the third-order piezomagnetic tensor, respectively. They are derived from the potential function H me as:

C = ∂ 2 SS H me , ε = ∂ 2 EE H me , µ = ∂ 2 HH H me , e = ∂ 2 SE H me , q = ∂ 2 SH H me , (13)
The behavior relations of the linear electro-magneto-mechanical coupling are then defined by:

T = -t e • E -t q • H + C E,H : S, ( 14
)
D = ε S • E + e : S, (15)
B = µ S • H + q : S, ( 16
)
where t (•) denotes the transpose operator. In ( 14)-( 16), the piezoelectric tensor e is taken to be null in the piezomagnetic phase and the piezomagnetic tensor q is taken to be null in the piezoelectric phase. The extrinsic nature of the strain-induced magnetoelectric effect considered here is highlighted by the absence of an eventual explicit magnetoelectric coefficient α linking the electric and magnetic fields in equations ( 15) and ( 16).

Choice of resolution variables 2.3.1. Variable of the mechanical problem

The deformations brought into play by the active materials considered are relatively low. Under the small deformation hypothesis, the linear strain tensor is defined by:

S = sym (∇u) = 1 2 ∇u + t ∇u , ( 17
)
where ∇ denotes the gradient operator. We chose u as the solving variable for the mechanical problem, as standard for FEM applied to mechanics.

Variable of the electrical problem

From (4) in a simply connected domain, the electric field can be expressed as,

E = -∇φ. (18)
The computation of the electric field via the scalar electric potential φ is a standard procedure in electrostatics. It involves a state variable, the electric voltage, with a physical meaning and that can be directly imposed as a boundary condition.

Variable of the magnetic problem

According to Helmholtz decomposition, the magnetic field H can be decomposed into two fields [START_REF] Ryu | Effect of the magnetostrictive layer on magnetoelectric properties in lead zirconate titanate/terfenol-d laminate composites[END_REF]: H 0 , the field created by the electric current density J s and H red the magnetic field created by the magnetized matter:

H = H 0 + H red . ( 19
)
Assuming that all electric currents are external to the active materials, it follows that:

∇ × H 0 = J s . ( 20
)
From ( 19)- [START_REF] Li | Enhancement in magnetic field sensitivity and reduction in equivalent magnetic noise by magnetoelectric laminate stacks[END_REF] together with (2), we get:

∇ × H red = 0. (21)
Therefore, for a simply connected domain there exists a reduced scalar field ϕ red such that:

H red = -∇ϕ red , ( 22
)
with ϕ red the magnetic reduced scalar potential. This potential is the independent variable of the magnetic problem. There are some clear advantages to the use of ϕ red as a resolution variable than a vector potential: it leads to shape functions which are nodal/scalar quantities and their gradients are easier to compute and integrate than vector shape functions and their curl. This choice of resolution variable can however lead to cancellation errors for problems with high magnetic permeability [START_REF] Phung | Automatic cuts for magnetic scalar potential formulations[END_REF]. This problem can be avoided by projecting H 0 onto edge elements [START_REF] Biro | Computation of 3-d magnetostatic fields using a reduced scalar potential[END_REF]. Such methods won't be necessary as the materials in question have low magnetic permeability.

Weak formulations of the electro-magneto-mechanical problem

From the mechanical and electromagnetic governing equations ( 1)-( 5) and the constitutive laws ( 14)-( 16) three sub-problems arise: an electrical, a mechanical and a magnetic problem which will be discretized separately. Boundary value problems associated with the electrical and mechanical subproblems are solved using the FE method, while the open-boundary value problem associated with the magnetic sub-problem is solved using FEM-BEM. In the following subsections, weak formulations of the electrical and mechanical coupled sub-problems are briefly recalled and the weak formulation of the magnetic problem is given in detail.

Weak form of the mechanical problem

Considering the coupled behavior law [START_REF] Eerenstein | Giant sharp and persistent converse magnetoelectric effects in multiferroic epitaxial heterostructures[END_REF], the decomposition of the magnetic field [START_REF] Ryu | Effect of the magnetostrictive layer on magnetoelectric properties in lead zirconate titanate/terfenol-d laminate composites[END_REF] and an appropriate virtual mechanical displacement vector field δu, which fulfills the homogeneous condition δu = 0 on ∂Ω u m , the weak formulation of the balance of linear momentum ( 5) is given by: find u such that:

Ωm δS : C : S dΩ - Ωpe δS : t e • E dΩ - Ωpm δS : t q • H red dΩ = Ωpm δS : t q • H 0 dΩ + Ωm δu • f dΩ + ∂Ω t m δu • f T d∂Ω m ∀δu, (23)
where, by definition, δS(x) = sym (∇δu(x)).

Weak form of the electrical problem

Though a FEM-BEM approach could be used for the electric problem with the surface term that accounts for the leaks of the displacement field, for our particular application, this coupling can be avoided for two reasons. Firstly, piezoelectric materials have high permittivity which allow them to canalize the displacement field inside the material domain, leading to low values of the stray field. Secondly, the piezoelectric material is in contact with two electrodes which can be regarded as equipotential surfaces, each one fixing the potential in a surface region of the piezoelectric material. The boundary of the active material can then be partitioned, without much loss of accuracy, into fixed potential regions and regions with no leaks, thus making the problem well posed for a FEM approach.

Considering the coupled behavior law [START_REF] Priya | Recent advancements in magnetoelectric particulate and laminate composites[END_REF] and an appropriate virtual electric scalar potential field δφ, which fulfills the homogeneous condition δφ = 0 on ∂Ω φ m , the FEM weak form of (3) reads: Find φ such that:

Ωm δE • ε S • E dΩ + Ωpe δE • e : S dΩ = 0 ∀δφ, (24)
with δE = -∇δφ.

Weak form of the magnetic problem

To establish the FEM-BEM weak formulation of the magnetic problem, the magnetic domain Ω is subdivided into an open exterior Ω 0 and an interior Ω m domains (Figure 2). The interface between the two magnetic subdomains, the boundary of the active material, is denoted by ∂Ω m .

Magnetic formulation inside the active material

To obtain a magnetic formulation inside the active domain Ω m , we consider the Maxwell-Thomson equation (1) weighted by an appropriate virtual magnetic scalar potential field δϕ and integrated by part on the domain Ω m :

Ωm ∇δϕ • B dΩ m - ∂Ωm δϕ B n d∂Ω m = 0 ∀δϕ, (25)
where B n = B • n with n the unitary normal vector perpendicular to the elementary surface d∂Ω m . Let's now replace B by the constitutive law ( 16) and H by the decomposition [START_REF] Ryu | Effect of the magnetostrictive layer on magnetoelectric properties in lead zirconate titanate/terfenol-d laminate composites[END_REF]. This results in the weak formulation of the Gauss-Thomson law (1) in the domain Ω m :

∂Ω m = ∂Ω 0 J s Ω 0 Ω pe Ω pm
Find (ϕ red , B n ) such that:

Ωm ∇δϕ • µ S • ∇ϕ red dΩ m - Ωpm ∇δϕ • q : S dΩ m + ∂Ωm δϕ B n d∂Ω m = Ωm ∇δϕ • µ S • H 0 dΩ m ∀δϕ. (26)
In a classical FEM approach, a sufficiently large free space domain is considered in order to neglect the surface term B n . A FEM-BEM coupling can take into account leakage fields and thus not have to explicitly consider a free space region.

Magnetic formulation in the free space domain

Let's consider magnetic Gauss law equation [START_REF] Eerenstein | Multiferroic and magnetoelectric materials[END_REF]. In Ω 0 the permeability is uniform, linear and taken equal to be the vacuum permeability. In this case, the outer magnetic reduced scalar potential fulfills the Laplace equation [START_REF] Bastos | Electromagnetic Modeling by Finite Element Methods[END_REF]:

∆ϕ red = 0. ( 27
)
Green's third identity can then be applied to ϕ red on a closed surface ∂Ω contained within the free space domain Ω 0 :

c ϕ red = ∂Ω ϕ red ∂G ∂n d∂Ω - ∂Ω G ∂ϕ red ∂n d∂Ω, ( 28
)
with G the Green's kernel fundamental solution of the Laplace equation defined as:

G(r) = 1 4πr in 3-D , ( 29
)
∂ϕ red ∂n = ∇ϕ red • n, ( 30
)
and c equal to

c = Θ 4π in 3-D (31)
with Θ the solid angle of the given point where this equation is formulated onto ∂Ω, in particular c = 1 2 if the equation is formulated on a smooth surface. In order to couple fields at the interface ∂Ω m between the active material and the free space domain, it is preferable to introduce the quantity B n in the strong formulation. Indeed, in the free space domain, this quantity is given by:

B n = µ 0 (H 0 • n) - ∂ϕ red ∂n , ( 32
)
which is equivalent to:

∂ϕ red ∂n = H 0n - B n µ 0 , ( 33
)
with H 0n = H 0 • n. Introducing this last relation into the previous strong formulation [START_REF] Wu | Inequivalence of direct and converse magnetoelectric coupling at electromechanical resonance[END_REF] we get:

c ϕ red = ∂Ω ∂G ∂n ϕ red d∂Ω - ∂Ω G H 0n - B n µ 0 d∂Ω. ( 34
)
Rearranging the terms, the strong formulation of the magnetic problem in the free space domain becomes:

-c ϕ red + ∂Ω G B n µ 0 d∂Ω + ∂Ω ϕ red ∂G ∂n d∂Ω = ∂Ω G H 0n d∂Ω. ( 35
)
As B n and ϕ red are considered continuous across ∂Ω m , the weak formulation in the free space domain can then be obtained by projecting [START_REF] Bichurin | Theory of low-frequency magnetoelectric coupling in magnetostrictivepiezoelectric bilayers[END_REF] onto an appropriate virtual magnetic scalar potential field δϕ 0 associated with a closed surface ∂Ω corresponding here to the external surface ∂Ω m of the active material domain:

Find (ϕ red , B n ) such that: ∂Ωm δϕ 0 ∂Ωm ϕ red ∂G ∂n d∂Ω m d∂Ω m + ∂Ωm δϕ 0 ∂Ωm G B n µ 0 d∂Ω m d∂Ω m - ∂Ωm δϕ 0 c ϕ red d∂Ω m = ∂Ωm δϕ 0 ∂Ωm G H 0n d∂Ω m d∂Ω m ∀δϕ 0 . (36)
In the absence of surface current density on ∂Ω m , the interface conditions for the magnetic fields along ∂Ω m read:

[ H × n ] ∂Ωm = 0, (37)
[

B • n ] ∂Ωm = 0, (38)
where [•] ∂Ω denotes the jump across the surface ∂Ω m . Taking ϕ red continuous across ∂Ω m ensures [START_REF] Dong | Equivalent circuit method for static and dynamic analysis of magnetoelectric laminated composites[END_REF]. The continuity of the normal magnetic induction field across ∂Ω m , translated by [START_REF] Benveniste | Magnetoelectric effect in fibrous composites with piezoelectric and piezomagnetic phases[END_REF], allows us to couple both magnetic formulations [START_REF] Harshe | Theoretical modelling of multilayer magnetoelectric composites[END_REF] and [START_REF] Dong | Longitudinal and transverse magnetoelectric voltage coefficients of magnetostrictive/piezoelectric laminate composite: theory[END_REF] which concludes the FEM-BEM coupling.

Computation of the source field

In a pure FEM approach, the source field is obtained in a pre-resolution on a mesh containing the entire computational domain. For problems with a large free space, this resolution increases the overall computational cost of the multiphysics problem. With the FEM-BEM coupling, the computation of the source field is performed by the use of the Biot-Savart law (39) which involves the numerical and/or analytical [START_REF] Rubeck | Analytical calculation of magnet systems: Magnetic field created by charged triangles and polyhedra[END_REF] integration of ∇G.

H 0 = ∇G × J s dΩ. ( 39
)
Hence, there will be no need to consider a free space domain and only the conductors carrying current and the studied device need to be meshed in order to compute the source field. In addition, from a numerical point of view, the computational cost of calculating the source field is independent of the distance between the field source and the active material domain. This is not the case for the pure FEM approach, for which a free space domain big enough to contain both domains has to be explicitly considered and meshed. This could mean a huge difference in computation times between the FEM and the proposed method.

Discrete formulations of the electro-magneto-mechanical problem

The following FEM-BEM formulations will be derived in vector-matrix notation, which results in a closed, efficient description of the implementation. Henceforth, all vectors and matrices will be signified with an underline as •. In addition, we will use Voigt's notation for the Cauchy stress tensor and the linear strain tensor: indicates the approximated domain and Ω e a finite element. The fields as well as their virtual counterparts are approximated element-wise by means of,

T = t (
u h , δu h = n node I=1 N I u u I , δu I = N e u u e , δu e , ( 41
)
φ h , δφ h = n node I=1 N I φ φ I , δφ I = N e φ φ e , δφ e , ( 42
)
ϕ h , δϕ h = n node I=1 N I ϕ ϕ I , δϕ I = N e ϕ ϕ e , δϕ e , ( 43
)
where u I , φ I and ϕ I denote, respectively, the discrete nodal displacement, the discrete nodal electric potential and the discrete nodal magnetic potential at node I, and u e , φ e , ϕ e are the associated element vectors of unknowns. In the same way, δu I , δφ I , δϕ I , δu e , δφ e and δϕ e denote the respective discrete virtual counterparts. N I u , N I φ and N I ϕ are shape functions associated with the node I and N e

• the corresponding matrix of shape functions. n node defines the number of nodes per element. In this framework, the mechanical strain, electrical and magnetic fields as well as their respective virtual counterparts as:

S h , δS h = n node I=1 B I u u I , δu I = B e u u e , δu e , ( 44
)
E h , δE h = n node I=1 B I φ φ I , δφ I = B e φ φ e , δφ e , ( 45
)
H h , δH h = n node I=1 B I ϕ ϕ I , δϕ I = B e ϕ ϕ e , δϕ e , (46)
with B I u , B I φ , B I ϕ , B e u , B e φ and B e ϕ the node-wise and element-wise B-matrices containing the Cartesian derivatives of the shape functions involved in the discrete form of ( 17), [START_REF] Zeng | Resonance magnetoelectric effect in bulk composites of lead zirconate titanate and nickel ferrite[END_REF] and [START_REF] Nan | Magnetoelectricity of multiferroic composites[END_REF].

Discrete form of the mechanical problem

The substitution of the approximations ( 41)- [START_REF] Tang | Micromechanical modeling of the multiphysical behavior of smart materials using the variational asymptotic method[END_REF] in the weak form (23) of the mechanical problem yields the following discrete form:

n elem e=1 t δu e            Ωe t B e u C
           = 0 ∀ δu e . ( 47
)

Discrete form of the electrical problem

In the same way, the substitution of approximations in the weak form [START_REF] Pan | Three-dimensional green's functions in anisotropic magneto-electro-elastic bimaterials[END_REF] of the electrical problem yields its discrete representation: (48)

Discrete form of the magnetic problem

The reduced magnetic scalar potential and the linear strain vector in the weak formulation [START_REF] Harshe | Theoretical modelling of multilayer magnetoelectric composites[END_REF] are discretized using ( 43) and [START_REF] Srinivas | The effective magnetoelectric coefficients of polycrystalline multiferroic composites[END_REF]. For the boundary term, from the discrete point of view, it is the set composed by the restriction to the facets belonging to ∂Ω h of the elements used for the discretization of Ω h . The interpolation of the reduced magnetic scalar potential is thus realized by the Nϕ shape functions associated to the restriction to the boundary of the N ϕ shape functions of the elements of Ω h . Similarly, for the normal component B n of the magnetic induction field, which, given the continuity condition, is interpolated by 0-order shape functions N0 . elem the number of facet elements belonging to ∂Ω h . The discretization of the BEM weak form (36) from the BEM formulation is performed with the same shape functions restricted to the boundary ∂Ω h . Here, as [START_REF] Dong | Longitudinal and transverse magnetoelectric voltage coefficients of magnetostrictive/piezoelectric laminate composite: theory[END_REF] is projected onto the Gauss points of the discretized surface, which is smooth at these positions, the coefficient c is taken equal to 1 2 . The resulting discretized form is:

n f elem e f =1 t δϕ e f 0              ∂Ωm t N e f 0 ∂Ωm ∂G ∂n N e f ϕ d∂Ω d∂Ω - 1 2 ∂Ωm t N e f 0 N e f ϕ d∂Ω k e f ∇nG ϕ e f red ∂Ωm t N e f 0 ∂Ωm G µ 0 N e f 0 d∂Ω d∂Ω k e f G B n e f - ∂Ω t N e f 0 ∂Ωm G H e f 0n d∂Ω d∂Ω sh ∂Ωm mag e f              = 0 ∀ δϕ e f .
(50) Contrary to the discrete formulations ( 47)- [START_REF] Buchanan | Layered versus multiphase magneto-electro-elastic composites[END_REF] where local interaction between elements translates into sparse matrices, the discretization (50) of the BEM formulation results into full matrices due to the double integration on all surface elements. From a numerical point of view, sparse matrices are stored in compressed row format (CRS) whereas BEM matrices are stored in full matrix format. Matrix compression techniques such as Fast Multipole Method (FMM) [START_REF] Greengard | A fast algorithm for particle simulations[END_REF] or Hybrid Cross Approximation (HCA) [START_REF] Börm | Hybrid cross approximation of integral operators[END_REF] can be applied to BEM matrices in order to reduce storage space.

Iterative multiphysics resolution of the block matrix assembly

A standard assembly procedure of the element matrices k e

• and the element right-hand sides sh e

• is applied to ( 47)-( 50) to generate the global matrix and right-hand side:

K • = n elem e=1 k e • and Sh • = n elem e=1 sh e • , (51)
where

n elem e=1
is an appropriate assembly operators. The two magnetic equations are treated as block and the fully coupled problem is made of a 3 × 3 block of matrices representing the magnetic, the mechanical and the electrical problems. The resulting assembled system of equations is presented in [START_REF] Avakian | Nonlinear modeling and finite element simulation of magnetoelectric coupling and residual stress in multiferroic composites[END_REF] where single-physics matrices are represented in color, sparse matrices are represented in light colors whereas full matrices are represented in dark colors.

                 0 K ϕbn 0 K ∇ n G K G 0 0 K uu K uφ K φu K φφ K ϕϕ K ϕu 0 K uϕ 0                              ϕ red B n u φ             =             Sh Ωm mag Sh ∂Ωm mag Sh meca Sh elec             (52)
This overall matrix system is not easy to solve using a single solver. Indeed, the full system is not symmetric and is made of both sparse (FEM) and full (BEM) matrices. Solvers are in general adapted to one type of matrix and not both. Additionally, there is a large scaling difference between coupling coefficients of the matrix K • . Indeed, elements of the stiffness matrix are computed using coefficients of the stiffness tensor of the order of 10-100 GPa, whereas the electric permittivity and magnetic permeability used to generate the electric and magnetic blocks are of the order of 10 -8 F/m and 10 -6 H/m, respectively. This difference in coupling coefficients translates into a big difference of eigenvalues of the global matrix K • , and therefore in a poor conditioning number of this matrix.

To overcome these problems, we implemented a block Gauss-Seidel scheme where the resolution of the global coupled system is treated as the resolution of a set of sub-problems [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF] [START_REF] Poblet-Puig | Modal-based prediction of sound transmission through slits and openings between rooms[END_REF]. As the global matrix of each sub-system is a diagonal block of the global matrix, they are equivalent to single-physic problems. The coupling is then introduced via the second hand term. The sub-system to be solved is given by ( 53):

[K ii ] { x n+1 i } = {Sh i } - i-1 j=1 [K ij ] { x n+1 j } - nb prob j=i+1 [K ij ] { x n j } (53)
where [K ii ] relates to single-physic matrices, [K ij ] to coupling matrices and { x n+1 i } to the solution at step n + 1 of the block Gauss-Seidel algorithm of the sub-system i. To solve each subsystem, a dedicated solver is used depending on the nature of the sub-problem to be solved. MUMPS direct solvers [START_REF] Amestoy | Multifrontal parallel distributed symmetric and unsymmetric solvers[END_REF] are used for the FE sub-problems associated with sparse matrices, with the LU decomposition performed once and for all at the beginning of the iterative resolution. The single-physics iterations for these problems are therefore fast.

The magnetic problem contains both sparse and full matrices and is therefore more difficult to solve. Direct solvers are prohibitive in terms of computing time, and iterative solvers such as GMRES are therefore preferred. A block preconditioner is used. For the FEM block, an incomplete LU (ILU) is used with a shift of the diagonal in order to avoid the singularity of the matrix as no reference potential is imposed. For the BEM blocs, if no compression technique is considered, a Jacobi preconditioner is used. HCA compression can also be used to obtain an approximation of both matrices K G and K ∇ n G . In this case, HCA-LU decomposition of K G can be obtained and used as a preconditioner for the BEM block, leading to faster convergence. If a nonlinear law µ(H) is considered, then a Newton-Raphson type solver for the magnetic problem can be easily integrated in the global resolution loop.

Validations and numerical results

Validation of the purely magnetic formulation

Firstly, we validated the magnetostatic formulation considering a magnetic sphere of radius of 1 mm discretized into increasingly finer meshes using the software GMSH [START_REF] Geuzaine | Gmsh : a three-dimensional finite element mesh generator with built-in pre-and post-processing facilities[END_REF]. The analytical solution for the magnetization of a sphere of uniform relative permeability µ r under an uniform exciting field H 0 is uniform and given by [START_REF] Durand | Électrostatique et magnétostatique[END_REF]:

H exact = 3 H 0 µ r + 2 . ( 54
)
Taking µ r = 10 and H 0 = 50 kA m -1 along the z-axis, the resulting magnetic field is also oriented along the z-direction, uniform and equal to H = 12.5 kA m -1 . Figure 3 shows the L 2 error in the computed magnetic field, ε L 2 (H) between the FEM-BEM solution H and the analytical solution H exact , with,

ε L 2 (X) = Ω ∥X -X exact ∥ 2 dΩ Ω ∥X exact ∥ 2 dΩ , ( 55
)
as a function of the number of degrees of freedom (DoFs). With the finest mesh tested of 69,386 we obtained an average field of 12, 499 A m -1 clearly showing that the FEM-BEM magnetostatic formulation is very accurate. The characteristics of this mesh are given in Table 2.

Validation of the magneto-mechanical coupling

To validate the magneto-mechanical formulation, we considered an unconstrained magnetostrictive sphere with a constant permeability µ r under a uniform magnetic field along the (z) direction, neglecting magnetic forces. Once again, the analytical magnetic field is given by [START_REF] Chen | A review of finite element open boundary techniques for static and quasi-static electromagnetic field problems[END_REF]. Defining the relative magnetic permeability µ r as the ratio between H and B (i.e., µ r = B/(µ 0 H)), we obtain the following system of equations used to obtain the analytical solution:

T = C : S -q t • H = 0, (56)
(2µ 0 + µ S ) H + q : S = 3µ 0 H 0 . ( 57
)
Using material parameters of Table 4 and coupling tensor structure from [START_REF] Do | Homogenization of Magnetoelectric 0-3 Type Composites by 3-D Multiphysics Finite-Element Modeling[END_REF], we compared analytical results to results of FEM-BEM with increasingly finely meshed spheres and obtained the curves of Figure 4 showing the L 2 error of the FEM-BEM solution vs the analytical solution, presented in Table 5. Both the magnetic and mechanical solutions converge, the characteristics of the finest mesh are given in Table 3. The magnetic and mechanical problems were solved using GMRES which can handle singular systems resulting from the imposition of the Neumann conditions on the entire boundary of the

Parameter

Value Young modulus (GPa) 100 Poisson coefficient 0.3 q 13 (N A -1 m -1 ) -30

q 33 (N A -1 m -1 ) 200 q 24 (N A -1 m -1 ) 60 q 15 (N A -1 m -1 ) 150 H 0 (A m -1 ) 50 • 10 3 µ r 10
mechanical problem. The tolerance of the mechanical solver was set to 10 -11 . All the multiphysics resolutions took 6 iterations to get a convergence of the block Gauss-Seidel algorithm inferior to 10 -6 for each single-physics solution and 8 iterations to get a tolerance of 10 -10 , and this despite the tolerance of the discretized mechanical problem being 10 -7 . Figure 5 shows the relative difference between partial solutions of the magnetic and mechanical problems as the function of the iteration number of the block Gauss-Seidel algorithm for the case of the most finely meshed sphere. At each block Gauss-Seidel iteration, the magnetic problem was solved first, followed by the mechanical problem.

Simplified analytical solution

To validate the fully coupled problem, we consider a composite structure made of a piezoelectric layer poled along the z-direction sandwiched between two electrodes, a reference electrode at potential 0 V and a floating one at unknown voltage resulting from the deformation of the piezoelectric layer. On top and on the bottom of the piezoelectric phase, two layers of magne- tostrictive materials poled along the x-direction drive the deformation of the piezoelectric layer. Figure 6 shows the model geometry used in numerical simulations.

The system works as follows : the source field H 0 produces the elongation of the magnetostrictive phase, which transmits the mechanical deformation to the piezoelectric phase thanks to the bonding between these two phases and the poling direction of the piezoelectric phase. A potential difference appears between the electrodes.

We compared the FEM-BEM results to an analytical solution obtained with the same assumptions used in [START_REF] Loyau | Magnetoelectric effect in layered ferrite/PZT composites. Study of the demagnetizing effect on the magnetoelectric behavior[END_REF] but with a different set of state variables, i.e., the same strain in the piezoelectric and piezomagnetic phases, zero strain along the z-axis, zero electric field along the z-direction inside the piezoelectric phase and zero current between electrodes, which corresponds to relations [START_REF] Nguyen | Finite element modeling of magnetic field sensors based on nonlinear magnetoelectric effect[END_REF].

where d = 1 mm is the thickness of the piezoelectric layer, and γ = 1/3 the volume fraction of the piezoelectric layer. We will solve this set of equations together with the previously discussed but rewritten constitutive laws ( 14)-( 16):

       T e 11 =
where the superscript • m refers to the value of fields inside the magnetostrictive phase and • e inside the piezoelectric phase. The solution to these equations is:

V out = e 31 q 12 + q 11 d (ε 33 (C H 11 + γ C E 11 + C H 12 + γ C E 12 ) + 2e 2 31 ) H 1 (60)
In order to compare this analytical solution to our simulations, we considered isotropic materials and the same mechanical properties for both the piezoelectric and magnetostrictive materials. We also didn't take into account demagnetizing fields, i.e., H 1 = H 0 which we took from 100 A m -1 to 100 A m -1 , while not exact, this is an acceptable approximation for materials with low relative permeability.

For simulations to be closer to the conditions of the analytical solution, we assume all the coefficients in the piezoelectric and piezomagnetic tensors to be null except of those appearing in [START_REF]Handbook of Giant Magnetostrictive Materials[END_REF], the coefficients used are presented in Table 6. Although close to the coefficients of PZT-5A and Terfenol-D, these coefficients are not expected to be accurate but to allow us to roughly validate our formulation. The slope of the curve in Figure 7 represents the DC magnetoelectric coefficient, or the output voltage divided by the exciting magnetic field, its value is presented in Table 7. For the analytic case, it is given by:

α = V out H 0 = e 31 q 12 + q 11 d ε 33 (C H 11 + γC E 11 + C H 12 + γC E 12 ) + 2e 2 31 (61)
Because of the demagnetizing field, the total field inside the magnetoelectric is lower than the source field (µ r > 1). So, in theory, supposing H 1 = H 0 overestimates the value of the output voltage and supposing equal strain in both piezoelectric and piezomagnetic phases also overestimates the output voltage. Indeed, in the simulations, it is the magnetic field which drives the deformation of the piezoelectric phase. Also, the 3D nature of the simulation (and therefore of the strain) vs the 2D strain hypothesis of the analytical formula also means that the output voltage of the fully modelled magnetoelectric should be lower than the analytical formula. A factor contributing to the high value of the output voltage of the magnetoelectric device with respect to the analytical solution is edge effects due to perfect corners in our geometry and linear material laws in our simulation. Otherwise, solutions are in very good agreement. The magnetic, electric and mechanical solutions resulting from FEM-BEM simulations for a source field of 16 kA m -1 are presented in Figures 8a-8d. In Figure 8a we see the amplified deformation of the device, as expected, more important in the magnetostrictive layer, in 8b we see µH, in Figure 8c the electric field appearing in the piezoelectric phase and in Figure 8d the electric potential, which translates into a voltage between the electrodes.

Comparison with the FEM

In order to position and compare the proposed approach to the FEM, a pure FEM approach was implemented in the software GetDP 3.5.0 [START_REF] Dular | A general environment for the treatment of discrete problems and its application to the finite element method[END_REF] -5 q 11 (N A -1 m -1 ) 200 q 12 (N A -1 m -1 ) -30 relative ε 33 (piezoelectric layer) 1800 relative ε 33 (piezomagnetic layer) 1 µ r (piezoelectric layer) 5 µ r (piezomagnetic layer) 9.5 A computer equipped with an Intel(R) Core(TM) i5-9500 CPU processor at 3.00GHz, with 48 GB RAM, was used for both FEM and FEM-BEM simulations. The study domains are those of Figure 1, therefore, a free space region and an infinite box were considered and the leaks of the magnetic field taken equal to zero at the boundary of the infinite box. The arrangement of both the coil and magnetoelectric device is shown in Figures 9a-9c. It consists of a magnetoelectric laminate structure identical to the one presented in Figure 6 and an un-centered coil. The source field was computed by the FEM by solving the preconditioned (LU decomposition) discretized form of the following weak formulation: Find H 0 such that:

Ω 0 ∇ × δH 0 • ∇ × H 0 dΩ 0 = Ω 0 ∇ × δH 0 • J s dΩ 0 , ∀δH 0 , ( 62
)
which comes from projecting (2) onto a virtual vector field ∇ × δH 0 and integrating by parts. Adding a Coulomb gauge [START_REF] Ledger | hp-finite element simulation of threedimensional eddy current problems on multiply connected domains[END_REF] the weak form of the computing of the source field reads: Find (H 0 , ϕ ′ ):

Ω 0 ∇ × δH 0 • ∇ × H 0 dΩ 0 + Ω 0 δH 0 • ∇ϕ ′ dΩ 0 = Ω 0 ∇ × δH 0 • J s dΩ 0 , ∀δH 0 , ( 63
) Ω 0 ∇δϕ ′ • H 0 dΩ 0 = 0, ∀δϕ ′ (64)
with the boundary condition ϕ ′ = 0 at ∂Ω 0 . Equations ( 63) and (64) were discretized using edge elements for H 0 and nodal elements for ϕ ′ . Having computed the source field, we obtained the discretized system describing the behavior of the magnetoelectric device presented in (65), 1-st order FEM were used. The same notations that in [START_REF] Avakian | Nonlinear modeling and finite element simulation of magnetoelectric coupling and residual stress in multiferroic composites[END_REF] were used for the FEM matrices, besides K ′ ϕϕ and Sh ′ Ω 0 mag which are now integrated over the domain Ω. All material parameters were taken from [START_REF] Do | Homogenization of Magnetoelectric 0-3 Type Composites by 3-D Multiphysics Finite-Element Modeling[END_REF]. The preconditioned (LU decomposition) linear system of equations (65) was solved using MUMPS direct solver [START_REF] Amestoy | Multifrontal parallel distributed symmetric and unsymmetric solvers[END_REF].

         K ϕu 0 K uu K uφ K φu K φφ K ′ ϕϕ K uϕ 0                  ϕ red u φ         =         Sh ′ Ω 0 mag Sh meca Sh elec         (65)
We considered two geometries and meshes differing in the size of the free space box and the infinite box (they will be referred as FEM 1 and FEM 2 respectively). The mesh of the magnetoelectric device and the coil is presented in Figures . and remained the same for the two simulations, they were produced with the software GMSH [START_REF] Geuzaine | Gmsh : a three-dimensional finite element mesh generator with built-in pre-and post-processing facilities[END_REF]. A coil centered at (0,0,0), of length of 3 mm with an inner radius of 15 mm and an outer radius of 18 mm centered around the point (30, 10, 0) mm. In the first case, a free space region with a radius of 316 mm and an infinite box with a radius of 474 mm were considered. The full mesh is presented in Figure 10.

For the first case (FEM 1), the magnetoelectric problem was applied to a mesh consisting of 83,921 nodes and 514,887 elements. For the computation of the source field, the matrix system to be solved involved 664,821 DoFs while the magnetoelectric problem involved 103,190 DoFs. The comparison between the resolution times, including the assembly, is presented in Table 8. We see that there is a major difference in the number of DoFs and in the total computing time of the multiphysic solution. The difference in the computing method of the source field contributes to the large difference in the total computing time due to the large number of DoFs of the magnetostatic resolution leading to the source field. In the FEM approach it is computed in the whole domain Ω, contrary to the FEM-BEM approach, in which it is only computed at the Gauss points of the magnetoelectric domain Ω m . Also, the FEM-BEM approach allocated a maximum of 3.1 GB of RAM vs 35.8 GB for FEM 1. Indeed, the storage space of the source field is considerably lower for the FEM-BEM approach as it is only stored at the Gauss points of Ω m , whereas for the FEM approach, it has to be computed and stored in the whole domain Ω. Also, HCA compression applied to the BEM matrices allowed for a compression ratio of 30 % on both matrices.

In Figure 11 we see the plots of the source field along the length of the laminate structure at its center, in Figure 12 the plot of the magnetic field and 9: Comparison between the maximum relative difference of the source field along (z), magnetic field along (z) and output voltage between the FEM approach and the proposed FEM-BEM approach.

Approach H 0 • e z (A m -1 ) H • e z (A m -1 ) output voltage (V) FEM 1 7.3 % 6.7 % 5.8 % FEM 2 6.9 % 6.2 % 5.2 % in Figure 13 the plot of the displacement. The output voltage of the device was of 43.78 mV for the FEM-BEM approach, vs 46.26 mV for the FEM approach, amounting to a difference of 5.78 % between the two approaches. While being very similar, we observe minor differences in the computed source field, for example, a maximum difference of its (z) component of 7.3 % at the center of the magnetoelectric device. Similar maximum relative differences are found in the other presented quantities, some of these differences are presented in Table 9. We see that the FEM approach tends to overestimate the computed source field with respect to the FEM-BEM approach, and therefore, we observe this trend over all the computed fields. For the FEM-BEM approach, the block Gauss-Seidel resolution algorithm converged in 14 iteration to a tolerance of 10 -10 amounting to a total resolution time of 5.7 s.

In order to obtain more precise results, a second mesh with a free space region of radius 443 mm and an infinite box of radius 664 mm was consid- ered. This mesh contains 108,927 nodes and 672,713 elements. We observe better agreement with the FEM-BEM solution, as seen in Table 9, but the computational cost of the complete problem, including the computing of the source field, is considerably higher: 2.5 times the computing time of FEM 1.

A notable difference was nevertheless seen in the amplitude of the edge effects, which were a lot more important in the FEM simulations than in the FEM-BEM approach: the FEM 1 solution gave a field with a norm of 317.9 A m -1 at the center of the element containing the edge closer to the coil vs 135.1 at the same point for the FEM-BEM approach, while being very similar elsewhere.

ME composite in an inhomogeneous field

Finally, we study the magnetoelectric composite structure shown in Figure 14 and fed by a volume inductor with coefficients from [START_REF] Do | Homogenization of Magnetoelectric 0-3 Type Composites by 3-D Multiphysics Finite-Element Modeling[END_REF]. The coil was purposely placed in such a way that the source field seen by the composite has no symmetries and cannot be approximated by a uniform field. the single-physics problems. After 15 iterations, the relative convergence of 10 -10 is achieved. In this particular test, the output voltage of the magnetoelectric is of -0.277 V.

Conclusion

In this paper, we proposed a FEM-BEM coupling strategy for the description of 3D magnetoelectric effects in composite structures. This coupling of numerical method allowed us only not to explicitly consider a free space domain, and thus to use a single mesh for the three sub-problems thus significantly reducing the number of DoFs. We validated our magnetostatic and magneto-mechanical formulation against analytical solutions. We succesfully validated our full coupled problem against a rough analytical solution and two FEM simulations with good agreement between the computed fields, lesser edge effects and a considerably shorter computation time. We also tested the presented full formulation of the magnetoelectric composite structure on a three layer device driven by a coil, with a non-linear material law for the magnetostrictive material and explored the performance of the block Gauss-Seidel algorithm for solving multiphysics problems in the previously described situation.

Figure 1 :

 1 Figure 1: Representation of the study domains of the electro-magneto-mechanical problem by the FEM.

Figure 2 :

 2 Figure 2: Subdomains of the magnetic problem based on a FEM-BEM formulation. The FEM is associated with the bounded domain Ω m = Ω pe ∪ Ω pm and the BEM with the exterior infinite domain Ω 0 corresponding to an infinite free space. The FEM and BEM are coupled on the common boundary ∂Ω m .

Figure 3 :

 3 Figure 3: L 2 error of the magnetic field H, ε L 2 (H) vs the analytical solution as a function of the number of DoFs of the discretized problem.

Figure 4 :

 4 Figure 4: ε L 2 (H), ε L 2 (S 11 ) and ε L 2 (S 33 ) vs the analytical solutions as function of the number of degrees of freedom of the discretized problem.

Table 5 :Figure 5 :

 55 Figure 5: Convergence of magnetic (o) and mechanical (x) solutions vs iteration number of the block Gauss-Seidel algorithm.

Figure 6 :

 6 Figure 6: Diagram of the composite structure (3×6×14 mm). For all simulations, it was discretized in 12×14×28 elements

Figure 7 :Figure 8 :

 78 Figure 7: Output voltage of the magnetoelectric composite vs source field. (x) corresponds to the analytical solution and (o) to the FEM-BEM results

Figure 9 :

 9 Figure 9: Mesh of the magnetoelectric device and the coil, seen from the three axes

Figure 10 :

 10 Figure 10: Mesh of the magnetoelectric device and the coil for the FEM 1 case.

Figure 11 :

 11 Figure 11: Comparison between the source fields H 0 between the FEM 1 and FEM-BEM approaches along the (z) component of the position inside the magnetoelectric device at its center.

Figure 13 :

 13 Figure 13: Comparison between the resulting displacements between the FEM 1 and FEM-BEM approaches along the (z) component of the position inside the magnetoelectric device at its center.

Figure 14 :

 14 Figure 14: The studied device is made of a coil conductor around the previously described composite structure.

Figure 15 :

 15 Figure 15: Convergence of magnetic (o), mechanical (x) and electric (+) solutions vs iteration number of the block Gauss-Seidel algorithm.

Table 1 :

 1 Magnetic, electric and mechanical fields and their SI-Units.

	Symbol Description	SI-unit
	u	Mechanical displacement	m
	f	Body forces	N
	S	Linear strain	-
	T	Cauchy stress	Pa
	φ	Electric scalar potential

  T 11 , T 22 , T 33 , T 12 , T 23 , T 13 ) and S = t (S 11 , S 22 , S 33 , 2S 12 , 2S 23 , 2S 13 ) . (40) The domain Ω m is discretized into a number n elem of finite elements resulting in a discrete counterpart such that: Ω ≈ Ω h =

	n elem	Ω e , where (•) h
	e=1

Table 2 :

 2 Size of the matrices of the finest mesh of the magnetostatic and magnetomechanical validations

	ϕ red	B n
	# of DoFs 49,876 9,757

Table 3 :

 3 Size of the matrices of the finest mesh of the magnetostatic and magnetomechanical validations

	ϕ red	B n	u
	# of DoFs 49,876 9,757 149,628

Table 4 :

 4 Material parameters used for the modeling of the piezomagnetic sphere.

Table 6 :

 6 Coefficients used for the comparison between simulation and the analytical solution

	Parameter	Value
	Young modulus (GPa)	70.3
	Poisson coefficient	0.345
	e 31 (C m -2 )

Table 7 :

 7 Comparison of magnetoelectric coefficients between simulation and analytical solution

	Analytical solution 2.3397e-04 mV A -1
	Simulation	2.2530e-04 mV A -1

Table 8 :

 8 Comparison between DoFs used for the computation of the source field, DoFs of the magnetoelectric problem (having computed the source field) and total resolution time, including the computation of the source field, between the two FEM simulations and the proposed FEM-BEM approach.

	Approach	FEM-BEM	FEM 1	FEM 2
	# nodes	5278	83,921	108,927
	# DoFs for computing H 0	-	664,821	869,672
	# DoFs of magnetoelectric problem	28,613	103,190	128,683
	total computing time	3 min 32 s 21 h 17 min 47 s 52 h 30 min 45 s
	maximum allocated memory	3.1 GB	35.8 GB	43.8 GB

Table

Acknowledgment

This work is supported by the Centre National de la Recherche Scientifique (CNRS) and the Université Grenoble Alpes, France.

The considered inductor has the following dimensions: a height of 2 mm, an inner radius of 7 mm and an outer radius of 9 mm. The total current flowing though the inductor was set to 100 A. The magnetostrictive phase is considered a nonlinear magnetic material with a relative permeability at the origin equal 100 and a saturation value of 2 T. We use MUMPS to solve the mechanical and the electrical problems. We also use a Newton-Raphson scheme with a tolerance of 10 -10 for the nonlinear magnetic problem, with GMRES for solving the resulting linearized system at each nonlinear iteration. The convergence of the magnetic problem was always achieved after two iterations of the Newton-Raphson. The stopping criteria of the Gauss-Seidel loop was set to relative 10 -10 difference between successive solutions. Within each Gauss-Seidel iteration, the magnetic problem was solved first, followed by the mechanical and the electrical problems. The convergence rate of individual single-physics problems is presented in Figure 15.

In Figure 15, from the 1 st to the 4 th iterations, we observe different convergence rates and beyond the 4 th , the rate of convergence is the same for all