A coupled electromagnetic-thermomechanical approach for the modeling of electric motors
AA coupled electromagnetic-thermomechanical approach for the modeling ofelectric motors
N. Hanappier a , E. Charkaluk a,b , N. Triantafyllidis a,b,c, ∗ a Laboratoire de M´ecanique des Solides (CNRS UMR 7649), Ecole Polytechnique, Institut Polytechnique de Paris b D´epartement de M´ecanique, ´Ecole Polytechnique, Route de Saclay, Palaiseau 91128, FRANCE c Aerospace Engineering Department & Mechanical Engineering Department (emeritus),The University of Michigan, Ann Arbor, MI 48109-2140, USA
Abstract
Future developments of lighter, more compact and powerful motors – driven by environmental andsustainability considerations in the transportation industry – involve higher stresses, currents and electro-magnetic fields. Strong couplings between mechanical, thermal and electromagnetic effects will consequentlyarise and a consistent multiphysics modeling approach is required for the motors’ design. Typical simulations– the bulk of which are presented in the electrical engineering literature – involve a stepwise process, wherethe resolution of Maxwell’s equations provides the Lorentz and magnetic forces which are subsequently usedas the external body forces for the resolution of Newton’s equations of motion.The work presented here proposes a multiphysics setting for the boundary value problem of electricmotors. Using the direct approach of continuum mechanics, a general framework that couples the electro-magnetic, thermal and mechanical fields is derived using the basic principles of thermodynamics. Particularattention is paid to the derivation of the coupled constitutive equations for isotropic materials under smallstrain but arbitrary magnetization. As a first application, the theory is employed for the analytical mod-eling of an idealized asynchronous motor for which we calculate the electric current, magnetic, stress andtemperature fields as a function of the applied current and slip parameter. The different components of thestress tensor and body force vector are compared to their purely mechanical counterparts due to inertia,quantifying the significant influence of electromagnetic phenomena.
Keywords:
Coupled thermo-mechanical and electromagnetic processes, Electric motors, Continuummechanics, Analytical solutions ∗ Corresponding author: [email protected]
Preprint submitted to JMPS Septembre 9, 2020 a r X i v : . [ phy s i c s . c l a ss - ph ] O c t ontents1 Introduction 22 Boundary value problem for electric motors 4 The increasing importance and market share of hybrid and purely electric vehicles, in the quest to reducetheir carbon footprint, urges the electric motor industry to develop higher performance products with reducedmanufacturing costs. New goals are set by various government agencies and industrial associations (L´opezet al., 2019) in terms of efficiency, reliability, power losses, power density, higher rotation velocity andreduced weight. Novel electric motor designs are needed to overcome these technological challenges in orderto comply with the aforementioned technical objectives and appropriate modeling tools must be developed.Modeling of electric motors has in the past been a topic studied predominantly by the electrical engi-neering community. The focus has been on the calculation of the magnetic field and resulting torque andiron losses for different motor designs using both analytical, (e.g. see: Boules (1984); Zhu et al. (1993);Lubin et al. (2011)) and numerical (e.g. see: Chari and Silvester (1971); Silvester et al. (1973); Abdel-Razek2t al. (1982); Arkkio (1987); Huppunen et al. (2004)) methods. A particular class of analytical methods,termed “subdomain methods” (e.g. see: Devillers et al. (2016)) constitute an approximate but efficient toolfor evaluating the magnetic characteristics of motor concepts at the preliminary design stage.In the late 90s, stress calculations in electric motors have appeared as a result of noise and vibrationsconcerns. As pointed out by Reyne et al. (1987), the first difficulty encountered was the evaluation of theelectromagnetic body forces, for which various authors gave different expressions, due to the absence of aconsistent continuum electrodynamics framework. The multiplicity of the different formulations, direct aswell as variational, for the thermomechanical-electromagnetic continuum, is a source of confusion. Different(albeit equivalent) expressions for the Maxwell stress and electromagnetic body forces can be obtained andare thus responsible for the difficulty in the correct modeling of stresses in electric motors. For furtherdiscussion on this issue, the interested reader is referred to the article by Kankanala and Triantafyllidis(2004) and book by Hutter et al. (2007). The first FEM computations for stresses in electric motors useda stepwise, uncoupled, approach: electric currents and magnetic fields where calculated using a purelyelectromagnetic model; the electromagnetic body force vector was then introduced as the external bodyforce in a purely mechanical model to calculate the resulting stress state (e.g. see: Reyne et al. (1988);Javadi et al. (1995)).The above-described approximate methods are inadequate to deal with the true multiphysics nature ofthe electric motor problem. The magnetic fields and currents generate the body forces driving the motor.These become even more important in the ferrous materials with high magnetic susceptibility that are usedto enhance and channel the magnetic flux for improved motor performance. Moreover, these materials haveintrinsic strongly coupled magnetic and mechanical behavior, with the material magnetization influencingthe stress state via the “ magnetostriction ” phenomenon and the stress state of the material also impactingits magnetization via “ inverse magnetostriction (Daniel et al., 2020). Moreover, strong currents influencetemperature due to ohmic effects and so on. Recognizing these issues, recent work by Fonteyn et al. (2010);Fonteyn et al. (2010a,b) takes into account magnetoelastic coupling effects for the numerical stress calculationin electric motors. However several approximations are used (e.g. a small strain approximation involvingnon frame-indifferent invariants and the angular momentum balance principle is not imposed), thermal fieldsare not considered and resulting stresses are not compared to inertial terms, motivating the present study.The goal of this work is a thermodynamically consistent formulation that couples the electromagnetic,thermal and mechanical effects for the boundary value problem of electric motors. On the theoreticalside, general continuum mechanics theories coupling thermomechanical and electromagnetic effects in solidsstarted back in the 1950s and 1960s. Although a literature review is beyond the scope of this study, a fewcomments are helpful to put in perspective the present work. As in Fonteyn et al. (2010a,b), the modelingapproach followed here is the “ direct ” method which uses conservation laws of continuum mechanics andthe thermodynamics procedure introduced by Coleman and Noll (1963) to obtain the problem’s governingequations and constitutive laws; a very readable account is presented in the book by Kovetz (2000). Forthe electric motor applications of interest the “ eddy current ” simplification of the problem is adopted (seeHiptmair and Ostrowski (2005) for a justification in linear materials) that neglects electric polarizationand displacement currents for low frequency electric fields. This theory is subsequently used to obtainthe analytical solution of an idealized asynchronous motor for which we calculate the electric current,magnetic, stress and temperature fields. The stress tensor and body force vector are compared to their purelymechanical counterparts due to inertia, quantifying the significant influence of electromagnetic phenomena,a novelty in this area to the best of the author’s knowledge.The presentation is organized as follows: following this introduction in Section 1, the general formulationfor the boundary value problem for electric motors is given in Section 2, where particular attention is paid tothe derivation of the coupled constitutive equations for isotropic materials under small strain but arbitrarymagnetization. The analytical model of an idealized asynchronous motor is presented in Section 3, where Also applied to the modeling of other electromagnetic problems such as Magneto-Rheological-Elastomers (e.g. seeKankanala and Triantafyllidis (2004); Dorfmann and Ogden (2003)) or Electro-Magnetic Forming processes (e.g. see Thomasand Triantafyllidis (2009)). Other applications use this approximation, such as “
Electromagnetic Forming ”; e.g. see Thomas and Triantafyllidis (2009).
3e calculate the magnetic field for the rotor and airgap in addition to the temperature field, the magnetic,total and elastic stresses in the rotor and the torque as a function of the applied current and the slipparameter (equivalent to the mechanical torque). The results for three different rotor materials (electricsteel, copper and aluminum) using realistic geometric and operational regime values and material parametersare presented in Section 4 and the work is concluded with a critical review and suggestions for future workin Section 5. The detailed derivations of the constitutive laws for isotropic materials under small strainbut arbitrary magnetization are given in Appendix A, detailed expressions for some elastic stress solutioncomponents are given in Appendix B while the determination of the magnetostrictive coefficient is presentedin Appendix C.
2. Boundary value problem for electric motors
The general formulation of the coupled electromagnetic-thermomechanical boundary value problem forelectric motors is presented in this section. Coordinate-free (dyadic) continuum mechanics notation is usedwith bold scripts referring to tensors, regular scripts to scalars. Eulerian fields are written using lowercaseletters, while capital letters are used for their Lagrangian counterparts. A superposed dot . f denotes the totaltime derivative of field f . The method adopted is the current configuration, direct approach of continuummechanics and tacitly assumes adequate smoothness of the fields involved. Unless stated otherwise, allfield quantities are functions of the current position x and time t . Although the governing equations forelectromagnetic continua are known (see Kovetz (2000)), for self-sufficiency and clarity of the work a briefpresentation is given in this section. The general equations of the problem can be distinguished in three groups, as presented in the subsec-tions below: electromagnetics (Gauss and Amp`ere), mechanics (conservation of mass, balance of linear andangular momenta) and thermodynamics (conservation of energy and entropy inequality). In the applica-tions considered the interfaces encountered are not moving with respect to matter, since they are either afree surface boundary or an interface between two different materials and hence in the sequel the interfacevelocity is the material velocity at the interface: v s = . x . relates the electric displacement d to the free electric charge density q . The differ-ential equation and the associated interface condition (in the absence of surface charges) are ∇ · d = q ; n · (cid:74) d (cid:75) = 0 , (2.1)where (cid:74) f (cid:75) denotes the jump of field f across a boundary/interface surface with an outward normal n . Maxwell-Amp`ere law links the h-field (sometimes also called the magnetic field ) h to the time-rate ofthe electric displacement d and the free total current density j . The differential equation and the associatedinterface condition are ∇ × h = ∂ d ∂t + j ; n × (cid:74) h (cid:75) + ( v s · n ) (cid:74) d (cid:75) = κ , (2.2)where v s is the velocity of the interface and κ the corresponding surface current density . The current density j consists of the conduction current density (cid:106) augmented by the convection of free electric charges q . x , i.e. j = (cid:106) + q . x . Maxwell-Faraday law relates the electric field e to the time-rate of the magnetic field (sometimes alsocalled the magnetic flux ) b . The differential equation and the associated interface condition are ∇ × e = − ∂ b ∂t ; n × (cid:74) e (cid:75) − ( v s · n ) (cid:74) b (cid:75) = 0 . (2.3)4 o magnetic monopole law confirms the absence of signed magnetic charges (monopoles) – hence thezero in its right-hand side, as compared to the Maxwell-Gauss law in (2.1). The corresponding differentialequation and the associated interface condition are ∇ · b = 0 ; n · (cid:74) b (cid:75) = 0 . (2.4)It should be mentioned here that the first set of two equations – Maxwell-Gauss and Maxwell-Amp`ere –result in the charge conservation principle ( ∇ · j + ∂q/∂t = 0) which thus need not be additionally enforced. Aether frame principle connects the fields ( d , h ) to ( e , b ). For the electric motor applications of interest,the polarization of the material is assumed negligible, in contrast to its magnetization (electric motorsinclude magnets and high permeability materials). The corresponding relations are d = (cid:15) e , h = 1 µ b − m , (2.5)where m is the magnetization (per unit volume) of the material and (cid:15) and µ are respectively the electricpermittivity and the magnetic permeability of free space. is described by the following differential equation ρ = ρJ = ⇒ . ρ + ρ ( ∇ · . x ) = 0 , (2.6)where ρ and ρ are respectively the reference and current mass densities and J ≡ det( ∂ x /∂ X ) the vol-ume change. In the absence of a discontinuity propagating in the continuum the corresponding inter-face/boundary condition gives no additional information. Linear momentum balance requires the introduction of the generalized electromagnetic-mechanical mo-mentum density g (e.g. Kovetz (2000)) – instead of . x for the purely mechanical problems – to be determinedsubsequently and gives the following differential equation and boundary/interface condition in the absenceof mechanical surface tractions ρ . g = ∇ · σ + ρ f , n · (cid:74) σ (cid:75) = 0 . (2.7)The body force per unit mass f contains only external, purely mechanical body forces , typically gravity.Electromagnetic forces are embedded in the total Cauchy stress σ and in g . Angular momentum balance
The non-reciprocity of actions/reactions in an electromagnetic - thermome-chanical continuum implies an asymmetric stress tensor, thus requiring the introduction of the generalizedmomentum g , resulting in the following relation for the asymmetric total stress σ ρ . x ∧ g = σ − σ T . (2.8)As a check we note that for a purely mechanical theory where g = . x , the Cauchy stress tensor is symmetric. Thus far the form of Maxwell laws and mechanics laws have the same expressions as in their correspondingpurely electromagnetic and purely thermomechanical counterparts; no electromagnetic body forces or bodytorques have been postulated. The coupling comes through the energy balance by adding an electromagneticenergy flux to the mechanical and thermal contributions, which allows us to find the missing constitutiveinformation involving the electromagnetic - thermomechanical coupling terms. We denote by ε the totalspecific energy of the continuum (i.e. mechanical, electromagnetic and thermal) and by η the specific entropyof the continuum , each defined at a point x and time t . The wedge product of two vectors a and b is an antisymmetric rank two tensor, defined by a ∧ b ≡ ab − ba . nergy conservation for the generalized electromagnetic-thermomechanical continuum in local form andits associated boundary condition give ρ . ε = ∇ · (cid:0) σ · . x − q − (cid:101) × (cid:104) (cid:1) + ρ ( f · . x + r ) ; n · (cid:74) σ · . x − q − (cid:101) × (cid:104) (cid:75) = 0 , (2.9)where r is the internal heat source per unit mass, q is the heat flux and (cid:101) × (cid:104) – also termed the Poyntingvector – is the electromagnetic energy flux , both fluxes leaving the continuum (hence their minus signs).The Poynting vector is the cross product of the electromotive force (cid:101) by the magnetotomotive force (cid:104)(cid:101) ≡ e + . x × b , (cid:104) ≡ h − . x × d . (2.10)Let η denote the specific entropy of the continuum at a point x and time t . Entropy production inequality , written here in terms of the continuum’s dissipation (cid:68) in local form andthe associated boundary condition are (cid:68) ≡ ρT . η − ρr + T ∇ · (cid:16) q T (cid:17) ≥ n · (cid:114) q T (cid:122) ≥ , (2.11)where T denotes the continuum’s absolute temperature field. Note that the adiabatic entropy source andthe adiabatic entropy flux have the same expressions as for the classical thermomechanics model: ρr/T and − ( n · q ) /T but η and q may now also depend on the electric and magnetic fields ( e , b ).The stage is now set to exploit the requirement of a positive dissipation by applying the method ofColeman and Noll (Coleman and Noll, 1963) in order to obtain the problem’s constitutive relations. Instead of working with the total specific energy of the continuum ε , following Kovetz (2000) we introducethe specific free energy of the solid ψ , a function of the thermodynamic state variables: . x , F ≡ ∂ x /∂ X the solid’s deformation gradient , b , (cid:101) , T , ∇ T and following Thomas and Triantafyllidis (2009) ξ , a set of internal variables associated with the mechanical and magnetic dissipative processes in the solid ψ ( . x , F , b , (cid:101) , T, ∇ T, ξ ) ≡ ε − T η − g · . x + 12 . x · . x − ρ (cid:20) (cid:15) e · e + 12 µ b · b − (cid:15) ( e × b ) · . x (cid:21) . (2.12)Using the Coleman and Noll procedure, constitutive relations are deduced for η, m , q , j , g , σ and ˙ ξ , interms of the thermodynamic state variables. These relations are distinguished in two categories: necessaryconstitutive relations (equalities) obtained from reversible restrictions – involving terms multiplying the ratesor gradients of the state variables that can assume arbitrary values – and sufficient constitutive relations(inequalities) deduced from non-reversible restrictions , and more specifically from the dissipation inequality,once its reversible terms are removed. Constitutive equalities give the following results (in addition to ∂ψ/∂ (cid:101) = ∂ψ/∂ ( ∇ T ) = ∂ψ/∂ . x = ) σ = ρ F · (cid:18) ∂ψ∂ F (cid:19) T + (cid:15) (cid:16) ee −
12 ( e · e ) I (cid:17) + 1 µ (cid:16) bb −
12 ( b · b ) I (cid:17) − (cid:16) bm − ( b · m ) I (cid:17) + . x (cid:15) ( e × b ) , m = − ρ ∂ψ∂ b , g = . x + 1 ρ (cid:15) ( e × b ) , η = − ∂ψ∂T . (2.13)Using the above results, in combination with (2.9) and (2.12), the dissipation inequality (2.11) yields (cid:68) = − ρ ∂ψ∂ ξ · . ξ + (cid:106) · (cid:101) − q T · ( ∇ T ) ≥ . (2.14) Note that in the absence of electromagnetic fields, g reduces to . x and ψ to the Helmholtz specific free energy ψ = (cid:117) − T η ,with (cid:117) = ε − / . x · . x ) the internal energy of the system, as expected from classical thermo-mechanics. onstitutive inequalities At this point no further details can be given about a generalized
Ohm’s law forthe conduction current density (cid:106) and a generalized
Fourier’s law for the heat flux q , on how they dependon the thermodynamic state variables, other than (2.14) has to be satisfied by (cid:106) = ˆ (cid:106) ( F , b , T, ∇ T, ξ , (cid:101) ) , q = ˆ q ( F , b , T, ∇ T, ξ , (cid:101) ) , (2.15)where it is assumed for simplicity that these vector fields are independent on . x . The well known formsof these relations require further assumptions about linearity and decoupling between different physicalmechanisms and will be discussed in Subsection 2.5.Using the above-obtained constitutive results from (2.13), we are now in a position to give a more concisethan in (2.12) expression for the solid’s free energy ρψ ( F , b , T, ξ ) = ρε − ρT η − ρ . x · . x − (cid:20) (cid:15) e · e + 12 µ b · b (cid:21) . (2.16)The above expression has a clear physical interpretation: the solid’s free energy density (per unit currentvolume) ρψ is obtained from the corresponding total energy density ρε of the continuum by subtracting thethermal contribution, the kinetic energy of the solid and the energy of the electromagnetic field.One final restriction must be recalled, that of material frame indifference which dictates the objectivity of ψ , i.e. its invariance under all translations and rigid body rotations of its arguments, dictating that ψ = ˆ ψ ( C , B , T, ξ ) ; B ≡ b · F , C ≡ F T · F , (2.17)where C the right Cauchy-Green deformation tensor. As it turns out, the use of ˆ ψ is the most convenientfor expressing the stress tensor and its subsequent simplification for small strains. An alternative formulation of the two last Maxwell laws, (2.3) and (2.4), involves the introduction of anelectric scalar potential φ and a magnetic vector potential ae = − ∇ φ − ∂ a ∂t ; b = ∇ × a . (2.18)As defined, the two potentials φ and a are not unique and a gauge condition needs to be additionallyenforced, such as the Coulomb gauge : ∇ · a = 0. For the problem at hand, the potential formulation leadsto a lower number of unknowns, thus justifying its introduction in (2.18). A convenient approximation for certain applications of electromagnetism (electric motors, electromag-netic forming etc.) is the eddy current approximation , which consists of ignoring the electric energy of theproblem as compared to its magnetic counterpart (e.g. Thomas and Triantafyllidis (2009)). This assumptionneglects the free electric charges (and hence Gauss’ equation (2.1)) and results in ignoring the displacementcurrent ∂ d /∂t and the convection of electric charges q . x , and hence j = (cid:106) , in Maxwell-Amp`ere’s law (2.2).The simplified Maxwell-Amp`ere’s law and boundary condition, recalling also (2.15) , reduce to ∇ × h = (cid:106) ; n × (cid:74) h (cid:75) = κ . (2.19)Note that the approximate charge conservation is ∇ · (cid:106) = 0, which is automatically satisfied given (2.19) .For the mechanical governing equations, the eddy current approximation implies that electric field termscan be ignored compared to their magnetic counterparts in the expression for the stress tensor in (2.13) andin the linear momentum density in (2.13) , which now reduces to the classical mechanics condition g = . x . No further assumption is made here about the constitutive equation for the internal variables.
7s a consequence, the angular momentum balance (2.1.2) now requires a symmetric total stress σ , as foundin (2.21). Note that other related works (see Fonteyn et al. (2010); Fonteyn et al. (2010a,b)) do not have asymmetric total stress.Taking into account (2.19), the simplified version of the linear momentum law (2.7) is rewritten as ρ .. x = ∇ · (cid:32) ρ F · ∂ ˆ ψ∂ C · F T (cid:33) + (cid:106) × b + m × ( ∇ × b ) + ( ∇ · m ) b + ρ f ; n · (cid:74) σ (cid:75) = 0 , (2.20)where (cid:106) × b are the Lorentz body forces, followed by the magnetic and the mechanical body forces.The constitutive equalities under the eddy current approximation, written in terms of the specific freeenergy density ˆ ψ in (2.17) take the form σ = 2 ρ F · ∂ ˆ ψ∂ C · F T + 1 µ (cid:16) bb −
12 ( b · b ) I (cid:17) − (cid:16) mb + bm − ( b · m ) I (cid:17) , m = − ρ F · ∂ ˆ ψ∂ B . (2.21)The remaining constitutive relations, i.e. Ohm’s and Fourier’s laws (2.15), the entropy constitutive equalityin (2.13) and the dissipation inequality in (2.14) remain unaltered.One more simplification is made possible by the eddy current approximation, consistent with ignoringthe electric energy of the system (and hence Gauss’s law), which allows the potential formulation for theelectric field e to be expressed only in terms of the magnetic potential vector ae = e app − ∂ a ∂t , (2.22)where e app is an externally applied electric field (typically to the coil that drives the system, e.g. see Thomasand Triantafyllidis (2009)). The magnetic field is still given by b = ∇ × a as in (2.18) and the electromotiveforce remains (cid:101) = e + . x × b , as defined in (2.10). The eddy current boundary value problem formulated thus far is general, for it accounts for nonlinearmagnetic and mechanical material response, both constitutive and kinematic (finite strains), as well asdissipative phenomena, i.e. plasticity, magnetic hysteresis etc., to be described by the evolution laws forthe internal variables ξ . We assume that a typical electric motor in its steady-state regime experiences onlysmall strains while it can also sustain large magnetizations, often up to saturation level. The implications ofthese restrictions on the selected specific free energy and the resulting expressions for the constitutive lawsare given progressively below, as more assumptions are introduced from one step to the next. Absence of mechanical and magnetic dissipation
We consider material behavior that does not includeplasticity or magnetic hysteresis, so internal variables ξ are not required for material description. As aresult, the specific free energy is a function of strain, magnetic field and temperature: ˆ ψ ( C , B , T ). Material isotropy
Isotropy of the material response implies that its specific free energy is a function ofsix invariants (and temperature), i.e. ˆ ψ ( C , B , T ) = ˆ ψ ( I , I , I , J , J , J , T ), where I i are the invariants ofthe right Cauchy-Green tensor C and J i are the coupled magneto-mechanical invariants of C and B . Decoupling of physical phenomena
It is assumed that thermo-mechanical, thermo-magnetic couplingscan be neglected, resulting in a separate thermal contribution ˆ ψ th constructed under the assumption of aconstant specific heat coefficient c (cid:15) . It is further assumed that, in the absence of magnetic fields, the freeenergy of the solid is ˆ ψ e ( I , I , I ) and that the magneto-mechanical coupling is described by the magneticinteraction energy ˆ ψ m ( J , J , J ).ˆ ψ ( C , B , T ) = ˆ ψ e ( I , I , I ) + ˆ ψ m ( J , J , J ) + ˆ ψ th ; ˆ ψ th = − c (cid:15) T [ln ( T /T ) − ,I = tr( C ) , I = (tr( C ) − tr( C · C )) , I = det( C ) ,J = B · C − · B , J = B · B , J = B · C · B , (2.23)8here T is a reference temperature.The implication of isotropy and decoupling on the generalized Ohm and Fourier laws in (2.15) is discussednext. We assume that the conduction current density (cid:106) depends solely on the electromotive force (cid:101) andthat the heat flux q is only a function of the temperature gradient ∇ T (cid:106) = γ ( (cid:107) (cid:101) (cid:107) ) (cid:101) ; q = − k ( (cid:107) ∇ T (cid:107) ) ∇ T , (2.24)where the scalar electrical conductivity γ ( (cid:107) (cid:101) (cid:107) ) > thermal conductivity k ( (cid:107) ∇ T (cid:107) ) >
0, asdictated by the dissipation inequality (2.14). The norm-dependence of these two scalar quantities is due tomaterial isotropy.
Small strain approximation
For the electric motor applications of interest here, we adopt the small strainapproximation, i.e. (cid:107) (cid:15) (cid:107) (cid:28)
1, where (cid:15) ≡ (1 / ∇ u + u ∇ ). Using Taylor series expansions in (cid:15) about thereference configuration of the quantities involved up to first order in (cid:15) and neglecting terms of order (cid:15) b ,we obtain a total stress σ as the sum of a purely elastic part e σ ( (cid:15) ) and a purely magnetic part m σ ( b ) σ = e σ + m σ ; e σ ≡ λ tr( (cid:15) ) I + 2 G (cid:15) , m σ ≡ µ (cid:20) bb −
12 ( b · b ) I (cid:21) − χ ( (cid:107) b (cid:107) ) µ ( (cid:107) b (cid:107) ) [ bb − ( b · b ) I ] + Λ( (cid:107) b (cid:107) ) µ ( (cid:107) b (cid:107) ) bb , m = χ ( (cid:107) b (cid:107) ) µ ( (cid:107) b (cid:107) ) b ; χ ( (cid:107) b (cid:107) ) µ ( (cid:107) b (cid:107) ) = − ρ (cid:34) ∂ ˆ ψ m ∂J + ∂ ˆ ψ m ∂J + ∂ ˆ ψ m ∂J (cid:35) C = I , Λ( (cid:107) b (cid:107) ) µ ( (cid:107) b (cid:107) ) = 2 ρ (cid:34) ∂ ˆ ψ m ∂J + 2 ∂ ˆ ψ m ∂J (cid:35) C = I , (2.25)where χ ( (cid:107) b (cid:107) ) is the material’s magnetic susceptibility , µ ( (cid:107) b (cid:107) ) = µ [1 + χ ( (cid:107) b (cid:107) )] its magnetic permeability andΛ( (cid:107) b (cid:107) ) a magneto mechanical coupling coefficient . It is important to note that at this stage our isotropicmaterial model is valid for small strains but arbitrary magnetization – the typical case of interest in magneticmotors – and that the corresponding magnetic susceptibility, magnetic permeability and magnetomechanicalcoupling coefficient are functions of the norm of the magnetic field b (due to isotropy). We should alsomention another consequence of small strain: the density equals its reference counterpart, i.e. ρ = ρ ,thus justifying its appearance (2.25). A remark is in order at this point about the expressions presented in(2.25); they differ from similar expressions presented by other authors (e.g. Aydin et al. (2017); Fonteynet al. (2010)) in view of our use of the objective invariants J k in our linearization procedure instead of theirsimplified, non-objective counterparts. The interested reader can find the details of these lengthy derivationsin Appendix A.
3. Application to an idealized asynchronous electric motor
This section pertains to the steady-state regime solution of an idealized, asynchronous electric motor,consisting of a cylindrical rotor and stator, as an application of the theory developed in Subsection 2.4.The solid cylindrical rotor geometry adopted here for the sake of the analytical treatment of the boundaryvalue problem, although uncommon in typical induction motors that have slots for conducting wires, isused for high frequency applications (see Gieras and Saari (2012)). The novelty here lies in the analyticalcomputation of the different body forces, stresses and temperature fields, performed using classical methodsof elasticity. The results obtained show how the analytical magnetic field computations presented by theelectrical engineering community (e.g. Lubin et al. (2011); Gieras and Saari (2012), can be complemented bymechanics. An added advantage of this simplified analytical model is its use as a benchmark for verificationin numerical codes. The small strain constitutive expressions that include terms order (cid:15) b and the justification for the omission of these terms in(2.25) are given in Appendix A. In the completely analogous – e → b , p → m , ε → µ − – electroelastic problems neglectingthe coupling terms is justified by assuming the small strain is of the same order as the square of the moderate electric fields,e.g. see Tian et al. (2012); Lefevre and Lopez-Pamies (2017). The elastic part of the free energy ˆ ψ e is independent of the magnetic field; upon linearization at C = I one obtains theclassical Lame constants λ and G appearing in (2.25). This coefficient gives the curvature of the strain vs magnetic field in a stress-free uniaxial magnetostriction experiment.
9o allow for an analytical solution, the motor geometry and the material behavior are considerably sim-plified using a 2D, plane strain framework and a homogeneous, linearized material response. The magneticsusceptibility χ and permeability µ , the magneto mechanical coupling coefficient Λ, the electrical conductiv-ity γ , the thermal conductivity k , the Lam´e constants λ , G and the mass density ρ are all given constants.Details for the setting of the corresponding boundary value problem are given below, where the unknownfields to be determined are the scalar magnetic potential a ( a = a e z ) in the rotor and the airgap, the rotor’stemperature field T and the Airy stress potential φ of the elastic stress field e σ . Figure 1: Cross-section of the simplified electric motor, indicating rotor, airgap and stator domains and corresponding frames.
The cross-section of the simplified induction motor is shown in Figure 1; the motor is considered in-finitely long in the normal to the plane and under plane strain conditions. It is composed of a cylindri-cal ferromagnetic rotor (domain (cid:68) : 0 ≤ r ≤ R ), surrounded by a cylindrical tubular stator (domain (cid:68) : R ≤ r ≤ R ), separated by an airgap (domain (cid:68) : R ≤ r ≤ R ). Two different polar coordinatesystems are used: the stator’s fixed reference frame (cid:83) ( r, θ s , z ) and the rotor’s moving frame (cid:82) ( r, θ, z ), where θ ≡ θ s − Ω t , with Ω the clockwise angular velocity of the rotor, as shown in Figure 1.Following Lubin et al. (2011), the motor is loaded by a current sheet of surface density κ perpendicularto the plane located on the internal radius of the stator. This current sheet models typical stator coils orwindings supplied by a poly-phased (usually three-phased) alternating electric current of angular frequency ω . The coils or windings are organized in p pairs per phase and the applied surface current density is κ = κ cos( pθ s − ωt ) e z , (3.1)with κ the oscillation’s amplitude in A/m . This current sheet rotates around the z -axis at the angularfrequency ω/p . It creates a rotating magnetic field of the same angular frequency, which triggers inducedcurrents at the rotor. The interaction of the induced currents in the rotor with the magnetic field createLorentz forces that result in the rotor spinning at an angular frequency Ω. Given that the phenomenonrelies on induction, an angular frequency differential exists between the stator field and the rotor: Ω < ω/p .We thus define the relative angular frequency ω r together with the slip parameter sω r = ω − p Ω ; s ≡ ω r ω , (3.2) Not to be confused with the electric potential, which is no longer needed. For simplicity only the fundamental time harmonic of the current supply is considered here. ω and Ω are constants, since the steady-state response of the motor is modeled.Some additional assumptions are necessary to solve the problem. i) Infinite permeability, rigid stator It is assumed that the stator’s strains are negligible – thus guaran-teeing a constant radius current sheet – and that it has an infinite permeability, i.e. µ −→ ∞ , resulting ina negligible stator h -field r > R : h = ( ∇ × a ) /µ ≈ . (3.3) ii) Constant temperature airgap The air in the airgap is assumed to be maintained at a constant tem-perature T a by forced ventilation. Due to ohmic losses the rotor temperature rises, but a convective heatexchange discharges its excess heat in the airgap. The corresponding radiation condition is r = R : q · e r = − k ( ∇ T ) · e r = h c ( T ( R ) − T a ) , (3.4)where h c is the convection coefficient and T is the rotor temperature field. iii) No external mechanical body forces No purely mechanical body forces, introduced in (2.7) are consid-ered, i.e. f = , since gravity effects are assumed negligible compared to inertia and magnetic contributions. iv) Constant velocity and acceleration Assuming a small slip s ( ω r (cid:28) Ω) and a small vibration amplitude,we can ignore the rates of the displacement . u and .. u in the velocity and acceleration terms, by keeping onlytheir Ω-dependent contributions, thus considerably simplifying the resulting algebra . x ≈ r Ω e θ , .. x ≈ − r Ω e r . (3.5)One consequence is a constant inertia term − ρ r Ω e r in the linear moment balance (2.20). The otherconsequence of (3.5) are the simpler expressions of the electromotive intensity (cid:101) defined in (2.10) and thematerial time derivative ˙ T , when expressed in the moving rotor frame (recall θ ≡ θ s − Ω t ) (cid:101) = (cid:20) − ∂ a ∂t + . x × b (cid:21) (cid:83) = − (cid:20) ∂ a ∂t + Ω ∂ a ∂θ s (cid:21) (cid:83) = − ∂ a ∂t (cid:12)(cid:12)(cid:12)(cid:12) (cid:82) , ˙ T = (cid:20) ∂T∂t + . x · ( ∇ T ) (cid:21) (cid:83) = (cid:20) ∂T∂t + Ω ∂T∂θ s (cid:21) (cid:83) = ∂T∂t (cid:12)(cid:12)(cid:12)(cid:12) (cid:82) . (3.6)Henceforth all equations are written in the rotor frame (cid:82) and all field quantities are functions of ( r, θ, t ).These governing equations and boundary conditions for the idealized, 2D motor are summarized below.( r, θ ) ∈ (cid:68) : ∇ × b = µγ (cid:101) ; ( R , θ ) : e r × (cid:74) h (cid:75) = , e r · (cid:74) b (cid:75) = , ( r, θ ) ∈ (cid:68) : ∇ × b = ; ( R , θ ) : e r × h = κ , κ = κ cos( pθ − ω r t ) e z , ( r, θ ) ∈ (cid:68) : ρ c (cid:15) ∂T∂t − k ∇ T = γ (cid:101) · (cid:101) ; ( R , θ ) : e r · [ k ( ∇ T )] = − h c ( T − T a ) , ( r, θ ) ∈ (cid:68) : ∇ · ( e σ + m σ ) = − ρ r Ω e r ; ( R , θ ) : e r · ( e σ + (cid:114) m σ (cid:122) ) = , ( r, θ ) ∈ (cid:68) ∪ (cid:68) : a = a ( r, θ, t ) e z , (cid:101) = − ∂ a ∂t , b = ∇ × a = 1 r ∂a∂θ e r − ∂a∂r e θ . (3.7)In the boundary condition for the equilibrium equation (3.7) there is no elastic stress field in the airgap e σ = , in contrast to the magnetic stress field m σ that exists in both the rotor and the airgap. The temperature field is only defined for the rotor, the T notation is not used and the subscript 1 is left out as superfluous. The elastic stress field is only defined in the rotor, the e σ notation is not used as unnecessary. ) External torque applied at rotor’s center To balance the moment produced by the shear stresses, itis assumed that an external mechanical torque is applied at the center line of the rotor ( r = 0) along the z -axis. The resulting torque per unit rotor length (cid:84) e z is (cid:84) = r (cid:90) π σ rθ ( r, θ )d θ , (3.8)and will be shown to be a constant, function of the relative angular frequency (cid:84) ( ω r ) with (cid:84) (0) = 0. To guide the physical interpretation of the results, the following dimensionless variables and parametersof the problem are introduced rR → r, ω r t → t, aµ κ R → a, k ( T − T a ) γω r ( µ κ ) R → T, σ ρ R Ω → σ , ζ ≡ ( R − R ) /R . (3.9)Henceforth, for simplicity the dimensionless variables and field quantities of the problem, r, t, a, T, σ aredenoted by the same symbol as their dimensioned counterparts.The governing equations and the associated interface and boundary conditions (in the rotor frame) aregiven below, starting with the magnetic potential a ∇ a = α ∂a ∂t , α ≡ µγω r R ; 0 ≤ r ≤ ,∂a ∂r = (1 + χ ) ∂a ∂r , ∂a ∂θ = ∂a ∂θ ; r = 1 , ∇ a = 0 ; 1 ≤ r ≤ ζ ,∂a ∂r = cos( pθ − t ) ; r = 1 + ζ . (3.10)The governing equation and boundary condition for the rotor’s temperature field T are (cid:70) − ∂T∂t − ∇ T = (cid:18) ∂a ∂t (cid:19) , (cid:70) ≡ kρ c (cid:15) ω r R ; 0 ≤ r ≤ , (cid:66) ∂T∂r + T = 0 , (cid:66) ≡ kR h c ; r = 1 , (3.11)with (cid:70) and (cid:66) the “ Fourier ” and “
Biot ” dimensionless coefficients respectively.Finally, the governing equations and boundary conditions for the rotor’s elastic stress field e σ are ∇ · e σ = f , f ≡ s j α ∂ a ∂t × ( ∇ × a ) − s m ∇ (cid:0) (cid:107) ∇ × a (cid:107) (cid:1) − r e r ; 0 ≤ r ≤ ,s j ≡ s χ , s m ≡ s χ + Λ1 + χ , s ≡ µ κ ρ R Ω . e σ rr = s (cid:34)(cid:18) ∂a ∂θ (cid:19) − (cid:18) ∂a ∂r (cid:19) (cid:35) − ( s j s m ) (cid:18) ∂a ∂θ (cid:19) + ( s j − s m ) (cid:18) ∂a ∂r (cid:19) ; r = 1 , e σ rθ = − s ∂a ∂θ ∂a ∂r + s j ∂a ∂θ ∂a ∂r ; r = 1 . (3.12) Only the radius of each domain of validity is recorded, since in all domains the angle θ ∈ [0 , π ) and the time t ∈ R + . Henceforth the rotor’s body force is denoted by f , taking the symbol used in (2.7) for the purely mechanical body force. is an equivalent of the “ Stuart ” number for magnetic fluids and gives the ratio of Maxwell over inertiastress magnitudes. The dimensionless coefficients s j and s m appearing in the expressions for the total stressin the rotor σ (sum of the elastic e σ and the magnetic m σ components respectively) depend on its magneticproperties while the total stress tensor in the airgap σ (Maxwell stress in vacuum) depends only on s .The corresponding expressions for the magnetic field and the total stress in each domain are given by σ = e σ + m σ , m σ = s j b b + ( s m − s j b · b ) I , b = ∇ × a ; 0 ≤ r ≤ , σ = m σ = s [ b b −
12 ( b · b ) I ] , b = ∇ × a ; 1 ≤ r ≤ ζ . (3.13)We first solve (3.10) to find the magnetic potential a , thus obtaining the ohmic dissipation for the heatequation (3.11), which is then used to determine the rotor’s temperature field T . The magnetic potentialgives the body forces for the linear momentum balance in (3.12), thus providing the rotor’s elastic field e σ . Solving the linear problem in (3.10) subject to the harmonic loading in (3.1), is more efficiently done inthe complex domain, where the magnetic potential a k ( r, θ, t ) takes the form a k ( r, Θ) = (cid:60) { ¯ a k ( r ) exp( − i Θ) } = A k ( r ) cos Θ + B k ( r ) sin Θ , k = 1 , ≡ pθ − t , (3.14)where ¯ a k ( r ) = A k ( r ) + iB k ( r ) is the complex magnetic potential amplitude that depends only on r .In the rotor domain, (3.10) results in a Bessel differential equation for the complex amplitude ¯ a ( r ) r d ¯ a d r + r d¯ a d r + (¯ α r − p )¯ a = 0 = ⇒ ¯ a ( r ) = ¯ AJ p (¯ αr ) ; ¯ α ≡ − iα , (3.15)where the constant α is defined in (3.10) and J p denotes a Bessel function of the first kind. The aboveexpression for ¯ a accounts for the fact that there is no singularity in r = 0, and hence explains the absenceof a Bessel function of the second kind in the general solution.In the airgap domain, (3.10) gives a Laplace equation for the complex amplitude ¯ a ( r ) r d ¯ a d r + r d¯ a d r − p ¯ a = 0 = ⇒ ¯ a ( r ) = ¯ Br p + ¯ Cr − p . (3.16)The complex-valued constants ¯ A , ¯ B and ¯ C appearing in (3.15) and (3.16) are determined using the interfaceand boundary conditions in (3.10), and are found to be¯ A = 2 (cid:104) J p (¯ α ) + ¯ (cid:103) , ¯ B = (cid:104) (cid:20) − (cid:103) [ J p (¯ α ) + ¯ (cid:103) ][1 + (1 + ζ ) − p ] (cid:21) , ¯ C = (cid:104) (cid:20) − (cid:103) [ J p (¯ α ) + ¯ g ][1 + (1 + ζ ) p ] (cid:21) , ¯ (cid:103) ≡ (cid:20) J p (¯ α ) − ¯ αp J p +1 (¯ α ) (cid:21) (cid:20) (1 + ζ ) p + (1 + ζ ) − p (1 + ζ ) p − (1 + ζ ) − p (cid:21)
11 + χ , (cid:104) ≡ (1 + ζ ) p [(1 + ζ ) p − (1 + ζ ) − p ] . (3.17)Using (3.17), the sought real amplitudes A k ( r ) and B k ( r ) in (3.14) are given in terms of their complexcounterparts found in (3.15) and (3.16), i.e. A k ( r ) = (cid:60) { ¯ a k ( r ) } , B k ( r ) = (cid:61) { ¯ a k ( r ) } ; k = 1 , Complex quantities are henceforth denoted by an overbar (¯). .4. Rotor temperature From the linearity of the governing equations for the temperature field in (3.11) and the magneticpotential solution in the rotor in (3.15), the forcing term in the conduction equation is found to be:( ∂a /∂t ) = 0 . B ( r )) + ( A ( r )) ] + 0 . B ( r )) − ( A ( r )) ] cos(2Θ) − A ( r ) B ( r ) sin(2Θ). The useof superposition and complex formulation lead to the following rotor temperature field T ( r, θ, t ) T ( r, Θ) = T ( r ) + (cid:60) (cid:8) ¯ T ( r ) exp( − i (cid:9) ; Θ ≡ pθ − t , (3.18)where the function T ( r ) is real and ¯ T ( r ) is complex. The real function T ( r ) is found from (3.11) to bed T d r + 1 r d T d r = − B ( r ) + A ( r )2 = ⇒ T ( r ) = c − (cid:90) r (cid:18) r (cid:90) r [ B + A ] r d r (cid:19) d r , (3.19)with the unknown constant c to be determined from the boundary condition.Solving for the complex function ¯ T ( r ) is reduced to solving a Bessel differential equation with a forcingterm through the superposition of a homogeneous and a particular solution ¯ T p ( r ), as follows r d ¯ T d r + r d ¯ T d r + (cid:0) ¯ β r − (2 p ) (cid:1) ¯ T = r ¯ a ⇒ ¯ T ( r ) = ¯ cJ p ( ¯ βr ) + ¯ T p ( r ) ; ¯ β ≡ − i (cid:70) , ¯ T p ( r ) = π (cid:20) − J p ( ¯ βr ) (cid:90) r Y p ( ¯ βr )¯ a r d r + Y p ( ¯ βr ) (cid:90) r J p ( ¯ βr )¯ a r d r (cid:21) , (3.20)where the unknown constant ¯ c in the homogeneous part of the solution will be specified from the boundarycondition. In solving (3.20) we made use of the fact that the solution is bounded at r = 0, and hence thereis no contribution from the Bessel function of the second kind Y p to the homogeneous part of the solution.However, Y p does enter under the integrals in the expressions for the particular solution ¯ T p ( r ) as seen above.Finally, the boundary condition at r = 1 in (3.11) splits into two boundary conditions: one for T ( r )that gives c and the other for ¯ T ( r ) that provides ¯ cc = 12 (cid:20)(cid:90) (cid:18) r (cid:90) r [( B ( r )) + ( A ( r )) ] r d r (cid:19) d r + (cid:66) (cid:90) [( B ( r )) + ( A ( r )) ] r d r (cid:21) , ¯ c = π (cid:34)(cid:90) Y p ( ¯ βr )¯ a ( r ) r d r − Y p ( ¯ β ) + (cid:66) ¯ βY (cid:48) p ( ¯ β ) J p ( ¯ β ) + (cid:66) ¯ βJ (cid:48) p ( ¯ β ) (cid:90) J p ( ¯ βr )¯ a ( r ) r d r (cid:35) , (3.21)where J (cid:48) p and Y (cid:48) p denote the derivatives of the first and second kind Bessel functions of order 2 p withrespect to their argument.Having determined T ( r ) and ¯ T ( r ), one can find from (3.18) the rotor temperature field T ( r, Θ).
The principle of superposition is used again for determining the rotor’s elastic stress field e σ . Recallingthe definitions for f in (3.12) and the solution for the magnetic potential a in (3.14) and (3.15), the bodyforces can be expressed as f ( r, Θ) = N ( r ) + ∇ V ( r, Θ) , where N ( r ) is not derivable from a potential Given the electromagnetic part of the forcing m f = − ∇ · m σ in (3.12), it is tempting to choose e σ = − m σ as a particular solutionto the electromagnetic forcing m f . However, this particular solution is ineligible as it does not satisfy the compatibility condition(see Barber (2009)), thus leading to the proposed approach. V ( r, Θ). f = N + ∇ V ; N = − s j α pr ( A + B ) e θ ; V ( r, Θ) = V ( r ) + V cs ( r, Θ) , V cs = V s ( r ) sin(2Θ) + V c ( r ) cos(2Θ) ,V ( r ) = − r s j α (cid:90) r ( A B (cid:48) − A (cid:48) B )d r − s m (cid:18) p r ( A + B ) + ( A (cid:48) + B (cid:48) ) (cid:19) ,V c ( r ) = − s j α A B − s m (cid:18) p r ( B − A ) + ( A (cid:48) − B (cid:48) ) (cid:19) ,V s ( r ) = s j α A − B )2 − s m (cid:18) p r A B + A (cid:48) B (cid:48) (cid:19) . (3.22)Consequently, the rotor’s elastic stress field e σ is decomposed as follows e σ ( r, Θ) = e σ N ( r ) + e σ V ( r, Θ) + e σ h ( r, Θ) , ∇ · e σ N = N , ∇ · e σ V = ∇ V , ∇ · e σ h = , (3.23)where each one of the constituent fields e σ N , e σ V , e σ h corresponds, in view of (2.25), to a compatible elasticstrain field, i.e. derivable from a displacement field. By abuse of terminology we call these elastic stressfields elastically compatible .Using the expression for N ( r ) from (3.22), an elastically compatible particular solution for e σ N ( r ) isfound by solving the tangential equilibrium ODE,d e σ N rθ d r + 2 r e σ N rθ = − s j α pr ( A + B ) = ⇒ e σ N rθ = − s j α pr (cid:90) r r ( A + B )d r . (3.24)An elastically compatible particular solution for the elastic stress field e σ V is found using the Airy stressfunction method in polar coordinates (see Barber (2009)). The components of e σ V can be expressed in termsof a stress potential φ V as follows e σ V rr = 1 r ∂φ V ∂r + 1 r ∂ φ V ∂ Θ + V , e σ V θθ = ∂ φ V ∂r + V , e σ V rθ = − ∂∂r (cid:18) r ∂φ V ∂ Θ (cid:19) ; ∇ φ V = − − ν − ν V . (3.25)The stress potential φ V is found (see footnote 17), by solving the Laplacian in (3.25) with the help of (3.22) φ V ( r, Θ) = − − ν − ν (cid:18)(cid:90) r r (cid:90) r V r d r d r + r p p (cid:90) r V cs r − p +1 d r − r − p p (cid:90) r V cs r p +1 d r (cid:19) . (3.26)The components of the elastically compatible homogeneous solution stress field e σ h are expressed in termsof the potential φ he σ h rr = 1 r ∂φ h ∂r + 1 r ∂ φ h ∂ Θ , e σ h θθ = ∂ φ h ∂r , e σ V rθ = − ∂∂r (cid:18) r ∂φ V ∂ Θ (cid:19) ; ∇ φ h = 0 . (3.27) Because we look for a particular solution only, integration constants are discarded. φ h in (3.27) we obtain φ h ( r, Θ) = Φ r θ + (cid:0) Φ c r p + Φ c r p +2 (cid:1) cos(2Θ) + (cid:0) Φ s r p + Φ s r p +2 (cid:1) sin(2Θ) . (3.28)The final expressions for e σ h , e σ V are obtained from (3.27) and (3.27). The six constants Φ , Φ , Φ c , Φ c , Φ s , Φ s are determined by the r = 1 boundary conditions in (3.12) e σ V rr (1) + e σ h rr (1) = s (cid:34)(cid:18) ∂a ∂θ (cid:19) − (cid:18) ∂a ∂r (cid:19) (cid:35) − ( s j s m ) (cid:18) ∂a ∂θ (cid:19) + ( s j − s m ) (cid:18) ∂a ∂r (cid:19) , e σ N rθ (1) + e σ V rθ (1) + e σ h rθ (1) = − s ∂a ∂θ ∂a ∂r + s j ∂a ∂θ ∂a ∂r . (3.29)From the decomposition in radial, cosine and sine terms (see footnote 18), result three equations for thenormal and three equations for the tangential boundary conditions, thus uniquely determining the soughtconstants. The full expressions for the stress at the rotor (elastic and magnetic components) can be thendetermined from (3.13) and (3.23) but are too cumbersome to be recorded here; the components of e σ V and e σ h are given in Appendix B. We are now in a position to give the expression for the torque/unit length (cid:84) . Recalling (3.8) and usingthe results for the stress field obtained above, one has (cid:84) = 4 πρ Ω R s p (cid:18) (cid:104) (cid:107) J p (¯ α ) + ¯ (cid:103) (cid:107) (cid:19) (cid:61) { J p (¯ α )¯ α ∗ J p +1 (¯ α ) ∗ } , (3.30)where ( ) ∗ denotes complex conjugation. This result gives the torque in terms of geometry, applied current(poles, amplitude and frequency), magnetic and electric properties and density of the rotor. Remarkably, (cid:84) is independent of the mechanical properties of the rotor, i.e. its shear modulus G and Poisson ratio ν .As the torque is a function of slip velocity ω r , it is instructive to find from (3.30) the initial slope of the (cid:84) ( ω r ) curve. Using asymptotics of the Bessel functions with respect to ¯ α for (cid:107) ¯ α (cid:107) = α = ω r µγR << (cid:84) ≈ ω r πγp (1 + p ) (cid:20) µ κ (1 + χ )(1 + ζ ) R (1 + ζ ) p + (1 + ζ ) − p + (1 + χ )[(1 + ζ ) p − (1 + ζ ) − p ] (cid:21) + O ( ω r ) . (3.31)One should keep in mind that the above expression gives only the initial slope of the (cid:84) ( ω r ) curve, butdepending on the problem, the range of validity of this linear approximation can be very small.
4. Results and discussion
Although we solve an idealized motor, the results presented here correspond to materials, geometries andoperating parameters found in the electrical engineering literature. The dimensionless quantities introducedin (3.9) allow a direct comparison of the results to related physically meaningful quantities. The solution is extracted from the general Michell (1899) solution. Given the form of the body forces and boundaryconditions, only those terms consistent with a solution of the form e σ = e σ ( r ) + e σ c ( r ) cos(2Θ) + e σ s ( r ) sin(2Θ) are kept. Also,all terms leading to stress singularities in r = 0 are excluded except for the Φ θ term required for the torque at r = 0. .1. Material, geometry and operating parameters The motor geometry and operating parameters used in the calculations are shown in Table 1. The studycovers three materials typically found in electric motors: electrical steel, copper and aluminum. Despite thedifferent motor architecture, the same values as in Lubin et al. (2011) are used whenever possible. The peakvalue of the current sheet is presently reduced to 1 . × A/m – from 8 × A/m in Lubin et al. (2011) – inorder to keep the maximum value of the magnetic field in the steel rotor below saturation, phenomenonnot accounted for here.Unfortunately, not all needed parameters can be found for a particular electric steel, thus requiring theuse of experimental data from the open literature for comparable materials. The value for the magneto-mechanical coupling coefficient Λ is fitted from Aydin et al. (2017), for the no-prestressed case, as detailedin Appendix C. A typical value for the magnetic susceptibility χ = 4000 for electric steel is adopted, whilethe elastic constants ν and E are taken from Belahcen et al. (2006). The rest of the material parameters –not given in Belahcen et al. (2006) and Aydin et al. (2017) – are taken from the open literature, as it is alsodone for the case of copper and aluminum, where we assume negligible magnetic effects ( χ = Λ = 0).The base case motor, which serves as a benchmark, is made of electric steel, has an airgap parameter ζ = 0 .
05 and a slip parameter s = 0 .
02. The rest of the geometric and operating parameters are kept fixed,independently of the rotor material, as shown below in Table 1.
Geometry
Rotor radius R ζ = ( R − R ) /R .
05 (base case)Number of pole pairs p Operating parameters
Peak value of current sheet κ . × A/mAngular velocity of current supply ω π rad/sSlip s = ω r /ω
2% (base case)External temperature T a ◦ CConvection coefficient h c
40 W/m /K Material properties Electrical steel Copper Aluminum
Electric conductivity γ . × S/m 5 . × S/m 3 . × S/mMagnetic susceptibility χ , ≈ ≈ − , ≈ ≈ ρ ,
650 kg/m ,
940 kg/m ,
700 kg/m Young’s modulus E × Pa 117 × Pa 69 × PaPoisson ratio ν .
34 0 .
33 0 . c (cid:15)
480 J/kg/K 385 J/kg/K 921 J/kg/KThermal conductivity k
45 W/m/K 397 W/m/K 225 W/m/K
Table 1: Motor geometry, operating parameters and rotor material properties
As discussed in Subsection 3.1, the equations are solved in the rotor frame (cid:82) and all field quantities arefunctions of ( r, Θ), where Θ = pθ − t and p the motor pole number (here taken p = 2). The results here are The chosen current sheet amplitude results in a maximum magnetic field of about 1 . T for the base case motor (steelrotor), roughly corresponding to the onset of magnetic field saturation for typical electrical steels (e.g. M400-50A), see Rekiket al. (2014).
17 snapshot of these rotating fields at t = 0 and are presented by plotting the corresponding field quantityat ( r, θ ). Magnetic field calculations for realistic geometries are routine for the electrical engineering community.The results for the current simple motor geometry are presented here solely for the purpose of explainingthe resulting force and strain fields.The magnetic field plots in Figures 2 and 3 show the contours of the dimensionless (normalized by µ κ )magnetic field (cid:107) b (cid:107) = ( b r + b θ ) / for three different values of the slip parameter s = 0 . , . , .
10 in thecase of a steel rotor with an airgap parameter ζ = 0 .
05. Notice that the magnetic field increases away fromthe center and peaks in a localized zone near the rotor periphery. As the slip s (equivalently the relativevelocity ω r ) increases, the localized high magnetization zone narrows, (e.g. see Jackson (1999) that the skindepth δ = (2 /γω r µ ) / ). The four localized magnetic field zones are a result of the number of poles ( p = 2).The high permeability of the rotor material ( χ = 4000 for electric steel) drastically increases its magneticfield, thus masking the variations of the considerably smaller – by one order of magnitude – strength of themagnetic field in the airgap in Figures 2. To remedy this, Figure 3 shows only the airgap magnetic field(hiding the rotor magnetic field) for the s = 0 .
02 slip motor of Figure 2(a). (a) (cid:107) b (cid:107) for s = 0 .
02 (b) (cid:107) b (cid:107) for s = 0 .
05 (c) (cid:107) b (cid:107) for s = 0 . Figure 2: Magnetic field norm (cid:107) b (cid:107) for a steel rotor (normalized by µ κ ), for different values of the slip parameter s .Figure 3: Magnetic field norm (cid:107) b (cid:107) in the airgap region (normalized by µ κ ) for the base case motor in Figure 2(a). The influence of changing motor geometry is presented in Figure 4 for three different airgap parameters ζ = 0 , , . , .
10 in a steel rotor and a slip value s = 0 .
02. As expected, Reducing the airgap size doesnot affect the distribution of the magnetic field, but increases drastically the maximum strength of the field.18 a) (cid:107) b (cid:107) for ζ = 0 .
02 (b) (cid:107) b (cid:107) for ζ = 0 .
05 (c) (cid:107) b (cid:107) for ζ = 0 . Figure 4: Magnetic field norm (cid:107) b (cid:107) for a steel rotor (normalized by µ κ ), for different values of the airgap parameter ζ . Comparison of the magnetic fields for different rotor materials is presented in Figure 5, where the resultsfor the high magnetic susceptibility steel are contrasted to the non-magnetic copper and aluminum rotors.The slip and airgap parameters are kept at their default value s = 0 . ζ = 0 . (a) (cid:107) b (cid:107) – steel (b) (cid:107) b (cid:107) – copper (c) (cid:107) b (cid:107) – aluminum Figure 5: Magnetic field norm (cid:107) b (cid:107) (normalized by µ κ ), for different rotor materials in motors with s = 0 . , ζ = 0 . Notice that for both the copper and aluminum rotors the maximum value of the magnetic field is twoorders of magnitude less than in steel. One can also observe that the normalized magnetic field for aluminumand copper reaches its maximum value at the rotor boundary, given the absence of magnetization in thesematerials. The slightly larger extent for the maximum magnetic field zone for the copper rotor, is attributedto its higher electrical conductivity which results in higher induced currents than in aluminum.
The full-field dimensionless temperature ( T − T a ) /T a → T ( r, Θ) for the base case steel rotor is presentedin Figure 6; the normalization with respect to the reference temperature T a adopted here as a more physicallymeaningful choice. Since the mean field dominates, the Θ-dependent variations are completely masked bythe scale used to plot Figure 6(a). The Θ-dependent variation (cid:60) (cid:8) ¯ T ( r ) exp( − i (cid:9) , whose amplitude is fourorders of magnitude lower than the mean, is plotted by itself in Figure 6(b). According to the values givenin Table 1 for the thermal characteristics of the idealized motor, the almost uniform temperature increase ofthe rotor is a mere 0 . ◦ C from an ambient airgap temperature of 20 ◦ C , with the maximum temperatureoccurring at the center. Here T denotes absolute temperature in ◦ K and not its normalized counterpart defined in (3.9). a) Full temperature field T ( r, Θ) (b) Θ-variation (cid:60) (cid:8) ¯ T ( r ) exp( − i (cid:9) Figure 6: Normalized temperature increase for the steel rotor (base case); (a) full field and (b) angular variation.
The influence of the rotor material on the dimensionless temperature increase T ( r, Θ) in the base casemotor is presented next in Figure 7. In comparing the results for steel in (a), copper in (b) and aluminumin (c), we notice that the temperature increase is almost uniform over the rotor, with the highest increase0 . ◦ C occurring in steel, 0 . ◦ C for copper and 0 . ◦ C for aluminum. (a) T ( r, Θ) – steel (b) T ( r, Θ) – copper (c) T ( r, Θ) – aluminum
Figure 7: Normalized temperature increase T for the base case motor: (a) steel, (b) copper and (c) aluminum rotors. Ohmic dissipation is the sole dissipation mechanism considered, as discussed in the first remark ofSubsection 2.5 and depends on the relative frequency ω r . The relatively low frequency used (about 1 Hz ,we consider ω r at 2% slip) explains the very low temperature increase found here. The dimensionless current density field j = j z = − γ ( ∂a/∂t ), (normalized by κ /R )for the base case motor is presented in Figure 8 for steel (a), copper (b) and aluminum (c) rotors, respec-tively. The currents for steel are forming thin plumes near the rotor surface because of the high magneticpermeability that concentrates the magnetic field at the rotor-airgap interface – see Figure 5(a) – limitingits penetration into the rotor. The current distribution for copper and aluminum rotors is very similar,given the absence of magnetization. Notice in Figure 8 that the maximum current values for steel are thelowest while the corresponding ones for copper are the highest, as expected by the different rotor materialconductivities according to Table 1. 20 a) j ( r, Θ) - steel (b) j ( r, Θ) - copper (c) j ( r, Θ) - aluminum
Figure 8: Current density j (normalized by κ /R ), for the base case motor: (a) steel, (b) copper and (c) aluminum rotors. Lorentz, magnetization and magnetostricive body forces
The different components of the magnetic bodyforce m f , defined as the divergence of the magnetic stress m σ in (2.25), are m f ≡ ∇ · m σ = j × b + m · ( b ∇ ) + Λ µ b · ( ∇ b ) ; m = χµ b , (4.1)where µ = µ (1 + χ ). The three different magnetic body force components in (4.1) are: the Lorentz bodyforce: j × b , a magnetization body force: m · ( b ∇ ) and a magnetostriction force: (Λ /µ ) b · ( ∇ b ). The last twocomponents are absent in non-magnetic copper and aluminum ( χ ≈ Λ ≈ ρ R Ω ), for the base case motor with a steel rotor case.The first important observation is that the Lorentz forces are negligible, with their maximum value of theorder of 1% of the inertial forces. A straightforward dimensional analysis indicates (cid:107) j (cid:107) ≈ (cid:107) b (cid:107) / ( µR ), giving (cid:107) j × b (cid:107) ≈ (cid:107) b (cid:107) / ( µR ) for the Lorentz component of the body force, compared to the magnetic χ (cid:107) b (cid:107) / ( µR )and magnetostrictive Λ (cid:107) b (cid:107) / ( µR ) components. (a) Lorentz (cid:107) j × b (cid:107) (b) Magnetization (cid:107) m · ( b ∇ ) (cid:107) (c) Magnetostriction (cid:107) (Λ /µ ) b · ( ∇ b ) (cid:107) Figure 9: Comparison of the different magnetic body forces (normalized by ρ R Ω ) for the base case motor with a steel rotor. Observe that the magnetization force is larger than its inertial counterpart – up to approximately fortytimes at the rotor’s edge due to the highest magnetic field gradients there, according to Figure 5(b) –pointing to the importance of accounting for magnetization body forces in electric motor models. Themagnetostrictive forces are not negligible and peak at about 160% of their inertial counterpart (or about5% of the maximum magnetization forces), a somewhat surprising result in view of the same order χ and Λcoefficients from Table 1 but explained by the different expressions for the corresponding forces in (4.1).21 a) (cid:107) m f (cid:107) – steel (b) (cid:107) m f (cid:107) = (cid:107) j × b (cid:107) – copper (c) (cid:107) m f (cid:107) = (cid:107) j × b (cid:107) – aluminum Figure 10: Comparison of the total magnetic body force (cid:107) m f (cid:107) (normalized by ρ R Ω ) for the base case motor with steel,copper and aluminum rotors. Notice that the magnetic body force is the Lorentz force j × b for the two non-magnetic materials. The results in Figure 10 compare the magnetic body force (normalized by ρ R Ω ) of the base motorfor the different rotor materials. Recall that the magnetic body force is just the Lorentz force for thecopper and aluminum rotors, in view of their negligible magnetic properties. We emphasize again the ordersof magnitude difference in the magnetic body force between the magnetic (steel) and the non-magnetic(copper, aluminum) materials. The Lorentz forces for the copper and aluminum rotor cases are comparable,given their close electric conductivity (see Table Table 1). Notice however that although the maximumcurrent density is higher in the better conducting copper, the corresponding maximum Lorentz force ishigher for the aluminum rotor. In order to better assess the influence of the electromagnetic effects on the total σ and elastic e σ stresses,we propose to compare them to the purely mechanical (only inertial body forces applied), plane strain elasticstress solution i σ for the spinning rotor of the base case motor under angular velocity Ω, a straightforwardlinear elasticity calculation resulting in the following stress field i σ rr = ρ R Ω (cid:18) − ν − ν − − ν − ν r (cid:19) , i σ rθ = 0 , i σ θθ = ρ R Ω (cid:18) − ν − ν − ν − ν r (cid:19) . (4.2)The maximum value for i σ rr and i σ θθ is [ ρ (3 − ν ) / − ν )]( R Ω) and occurs at the rotor’s center r = 0.For a more meaningful comparison to the purely mechanical stresses due to inertial effects, all future stressresults are normalized by this maximum value, instead of ρ ( R Ω) used thus far.22 a) σ rr ( r, Θ) (b) σ rθ ( r, Θ) (c) σ θθ ( r, Θ) Figure 11: Dimensionless total stresses in rotor and airgap (normalized by the maximum inertial stress): (a) normal, (b) shearand (c) hoop, for the base case steel motor.
The normalized total stress components for the base case motor with the steel rotor are presented inFigure 11, – with the stress fields shown both in the rotor and the airgap – where one can see the continuityof the normal σ rr and shear σ rθ components at the rotor-airgap interface.The total normal stress σ rr is always positive, never exceeding the maximum, purely inertial value, asseen in Figure 11(a). It monotonically increases away from the rotor’s edge and reaches its maximum atthe center, region where the electromagnetic effects are negligible, in contrast to the rotor’s edge. The totalshear stress σ rθ varies symmetrically between approximately ±
5% of the maximum (normal) inertial stress ,following the angular pattern imposed by the cos(2Θ) and sin(2Θ) terms. Also notice in Figure 11(b) thesingularity in r = 0 – truncated in the figure – due to the external torque applied there. The total hoopstress σ θθ is positive in most of the central domain, where inertial effects dominate, with the same maximumvalue as for the purely inertial case. The influence of the magnetic field is however evident on the rotor’sedge, where a compressive stress of the same absolute value as the maximum inertial stress does appear. (a) e σ rr ( r, Θ) (b) e σ rθ ( r, Θ) (c) e σ θθ ( r, Θ) Figure 12: Dimensionless elastic stresses in rotor (normalized by the maximum inertial stress): (a) normal, (b) shear and (c)hoop, for the base case steel motor.
The normalized elastic stress e σ components in the rotor are given in Figure 12 and differ significantlyfrom their total stress counterparts σ , as a simple comparison between Figure 11 and Figure 12 shows. Theelastic stress components are approximately their inertial counterparts i σ , given by (4.2), due to the weak The rotor has no shear stresses for the purely inertial loading; plotting the shear stress over the maximum value of theinertial stress (which corresponds to the radial and hoop stresses) allows the comparison of its magnitude with respect to thenormal stresses.
Figure 13: Dimensionless torque (cid:84) (normalized by πρ R Ω s ) vs slip coefficient s = ω r /ω for the base case motor with threedifferent rotor materials. The torque (cid:84) , normalized by πρ R Ω s , is plotted in Figure 13 as a function of the slip coefficient s . For low κ values (where the magnetic field remains below the saturation level for steel for all slip valuesconsidered ) the steel rotor shows higher torque than its copper and aluminum counterparts across almostall the slip range, only slightly dominated by the copper rotor in a region around 5 −
10% slip.For high κ values, the monotonic increase of the torque as a function of slip for steel – due to its linearmagnetic response – is misleading, as saturation may occur, which is not accounted for in the model. In thebase case motor, the magnetic field for the steel rotor is already close to saturation for κ = 1 . × , s = 2%with a value of 1 .
3T (see Figure 2(a)). In this case, it is expected that due to magnetic saturation, thesteel torque-slip curve above s = 2% should be reaching a maximum torque, as is the case for the copperand aluminum rotors. For s = 5% or higher, the copper motor would produce a larger torque than its steelcounterpart.
5. Conclusion
Using the direct approach of continuum mechanics, based on Kovetz (2000), a general framework thatcouples the electromagnetic, thermal and mechanical effects is derived and subsequently applied to formulatethe boundary value problem for electric motors. Particular attention is paid to the derivation of the coupledconstitutive equations for isotropic materials under small strain but arbitrary magnetization. As a firstapplication, the theory is employed for the analytical modeling of an idealized asynchronous motor forwhich we calculate the magnetic, thermal, stress fields and its torque. To better assess the influence ofmagnetization on stresses, three different rotor materials are examined: electric steel, copper and aluminum The normalization quantity is the product of the rotor’s area πR by the electromagnetic stress term ρ R Ω s = µ κ . As shown in (4.2), the peak value of the magnetic field increases with slip.
ACKNOWLEDGMENTS
The work of N. H. is supported by a Fellowship from the
Andr´e Citro¨en Chair of the Ecole Polytechnique.
6. ReferencesReferences
Abdel-Razek, A., Coulomb, J., Feliachi, M., Sabonnadiere, J., 1982. Conception of an air-gap element for the dynamicanalysis of the electromagnetic field in electric machines. IEEE Transactions on Magnetics 18, 655–659. URL: https://ieeexplore.ieee.org/abstract/document/1061898 .Arkkio, A., 1987. Analysis of induction motors based on the numerical solution of the magnetic field and circuit equations ,97URL: http://urn.fi/urn:nbn:fi:tkk-001267 .Aydin, U., Rasilo, P., Martin, F., Singh, D., Daniel, L., Belahcen, A., Rekik, M., Hubert, O., Kouhia, R., Arkkio, A.,2017. Magneto-mechanical modeling of electrical steel sheets. Journal of Magnetism and Magnetic Materials 439, 82– 90. URL: , doi: https://doi.org/10.1016/j.jmmm.2017.05.008 . arber, J., 2009. Elasticity. Solid Mechanics and Its Applications, Springer Netherlands. URL: https://books.google.fr/books?id=5M9j319PbKMC .Belahcen, A., Fonteyn, K., Fortino, S., Kouhia, R., 2006. A coupled magnetoelastic model for ferromagnetic materials. Proc.of the IX Finnish Mechanics Days. von Hertzen R., Halme T.(eds.) , 673–682.Boules, N., 1984. Two-dimensional field analysis of cylindrical machines with permanent magnet excitation. IEEE Transactionson Industry Applications IA-20, 1267–1277.Brown, W.F., 1966. Magnetoelastic Interactions. Springer-Verlag, New York.Chari, M.V.K., Silvester, P., 1971. Analysis of turboalternator magnetic fields by finite elements. IEEE Transactions on PowerApparatus and Systems PAS-90, 454–464. URL: https://ieeexplore.ieee.org/abstract/document/4074358 .Coleman, B.D., Noll, W., 1963. The thermodynamics of elastic materials with heat conduction and viscosity. Archive forRational Mechanics and Analysis 13, 167–178. URL: https://doi.org/10.1007/BF01262690 , doi: .Daniel, L., Bernard, L., Hubert, O., 2020. Multiscale Modeling of Magnetic Materials. Elsevier.Daniel, L., Hubert, O., 2009. An analytical model for the δ e effect in magnetic materials. The European Physical JournalApplied Physics 45, 31101.Daniel, L., Hubert, O., Ossart, F., Billardon, R., 2003. Experimental analysis and multiscale modelling of the anisotropicmechanical and magnetostrictive behaviours of electrical steels, in: Journal de Physique IV (Proceedings), EDP sciences.pp. 247–254.Devillers, E., Le Besnerais, J., Lubin, T., Hecquet, M., Lecointe, J., 2016. A review of subdomain modeling techniques inelectrical machines: Performances and applications, in: 2016 XXII International Conference on Electrical Machines (ICEM),pp. 86–92. URL: https://ieeexplore.ieee.org/abstract/document/7732510 .Dorfmann, A., Ogden, R., 2003. Magnetoelastic modelling of elastomers. European Journal of Mechanics - A/Solids 22, 497– 507. URL: , doi: https://doi.org/10.1016/S0997-7538(03)00067-6 .Fonteyn, K., Belahcen, A., Kouhia, R., Rasilo, P., Arkkio, A., 2010a. Fem for directly coupled magneto-mechanical phenomenain electrical machines. IEEE Transactions on Magnetics 46, 2923–2926.Fonteyn, K.A., Belahcen, A., Rasilo, P., Kouhia, R., Arkkio, A., 2010b. Contribution of maxwell stress in air on the deformationsof induction machines, in: 2010 International Conference on Electrical Machines and Systems, pp. 1749–1753.Fonteyn, K.A., et al., 2010. Energy-based magneto-mechanical model for electrical steel sheets .Gieras, J.F., Saari, J., 2012. Performance calculation for a high-speed solid-rotor induction motor. IEEE Transactions onIndustrial Electronics 59, 2689 – 2700.Hiptmair, R., Ostrowski, J., 2005. Coupled boundary-element scheme for eddy-current computation. Journal of EngineeringMathematics 51, 231–250. URL: https://doi.org/10.1007/s10665-004-2116-3 , doi: .Huppunen, J., et al., 2004. High-speed solid-rotor induction machine–electromagnetic calculation and design URL: https://lutpub.lut.fi/handle/10024/36551 .Hutter, K., Ven, A.A., Ursescu, A., 2007. Electromagnetic Field Matter Interactions in Thermoelasic Solids and Viscous Fluids.volume 710. Springer.Jackson, J.D., 1999. Classical Electrodynamics, 3rd ed. John Wiley & Sons, Inc.Javadi, H., Lef`evre, Y., Cl´enet, S., Mazenc, M., 1995. Electro-magneto-mechanical characterizations of the vibration ofmagnetic origin of electrical machines. IEEE transactions on magnetics 31, 1892–1895. URL: https://ieeexplore.ieee.org/abstract/document/376408 .Kankanala, S., Triantafyllidis, N., 2004. On finitely strained magnetorheological elastomers. Journal of the Mechanicsand Physics of Solids 52, 2869 – 2908. URL: ,doi: https://doi.org/10.1016/j.jmps.2004.04.007 .Kovetz, A., 2000. Electromagnetic theory. volume 975. Oxford University Press Oxford.Lefevre, V., Lopez-Pamies, O., 2017. Homogenization of elastic dielectric composites with rapidly oscillating passive and activesource terms. SIAM Journal on Applied Mathematics 77, 1962–1988.L´opez, I., Ibarra, E., Matallana, A., Andreu, J., Kortabarria, I., 2019. Next generation electric drives for hev/ev propulsionsystems: Technology, trends and challenges. Renewable and Sustainable Energy Reviews 114, 109336. URL: , doi: https://doi.org/10.1016/j.rser.2019.109336 .Lubin, T., Mezani, S., Rezzoug, A., 2011. Analytic calculation of eddy currents in the slots of electrical machines: Applicationto cage rotor induction motors. IEEE Transactions on Magnetics 47, 4650–4659. doi: .Michell, J.H., 1899. On the Direct Determination of Stress in an Elastic Solid, with application to the Theory of Plates. Pro-ceedings of the London Mathematical Society s1-31, 100–124. URL: https://doi.org/10.1112/plms/s1-31.1.100 , doi: , arXiv:http://oup.prod.sis.lan/plms/article-pdf/s1-31/1/100/4637069/s1-31-1-100.pdf .Pao, Y.H., Yeh, C.S., 1973. A linear theory for soft ferromagnetic elastic solids. International Journal of Engineering Science11, 415 – 436. URL: , doi: https://doi.org/10.1016/0020-7225(73)90059-1 .Rekik, M., Hubert, O., Daniel, L., 2014. Influence of a multiaxial stress on the reversible and irreversible magnetic behaviourof a 3% si-fe alloy. International Journal of Applied Electromagnetics and Mechanics 44, 301–315.Reyne, G., Sabonnadiere, J., Coulomb, J., Brissonneau, P., 1987. A survey of the main aspects of magnetic forces andmechanical behaviour of ferromagnetic materials under magnetisation. IEEE Transactions on Magnetics 23, 3765–3767.URL: https://ieeexplore.ieee.org/abstract/document/1065518 .Reyne, G., Sabonnadiere, J.C., Imhoff, J.F., 1988. Finite element modelling of electromagnetic force densities in dc machines.IEEE Transactions on Magnetics 24, 3171–3173. URL: https://ieeexplore.ieee.org/abstract/document/92371 .Silvester, P., Cabayan, H.S., Browne, B.T., 1973. Efficient techniques for finite element analysis of electric machines. EEE Transactions on Power Apparatus and Systems PAS-92, 1274–1281. URL: https://ieeexplore.ieee.org/abstract/document/4075205 .Thomas, J., Triantafyllidis, N., 2009. On electromagnetic forming processes in finitely strained solids: Theory and examples.Journal of the Mechanics and Physics of Solids 57, 1391 – 1416. URL: , doi: https://doi.org/10.1016/j.jmps.2009.04.004 .Tian, L., Tevet-Deree, L., DeBotton, G., Bhattacharya, K., 2012. Dielectric elastomer composites. Journal of the Mechanicsand Physics of Solids 60, 181–198.Zhu, Z.Q., Howe, D., Bolte, E., Ackermann, B., 1993. Instantaneous magnetic field distribution in brushless permanentmagnet dc motors. i. open-circuit field. IEEE Transactions on Magnetics 29, 124–135. URL: https://ieeexplore.ieee.org/abstract/document/195557/references . Appendix A. Isotropic, small strain, arbitrary magnetization constitutive laws
The derivation of the constitutive laws for an isotropic magnetoelastic material for small strain (cid:15) , but ar-bitrary magnetic field b , although straightforward requires lengthy calculations. Although such calculationshave been presented in the literature a long time ago by Pao and Yeh (1973), following the early works onmagnetoelasticity by Brown (1966), a direct comparison with our results is not possible due to the differentformulations adopted (e.g. different independent variables of the free energy densities, different definitionsof total stress etc.). Moreover, such derivations are not always done consistently in the available literature; alinearized version of the invariants is often considered, thus violating the objective nature of the free energysince the small strain tensor (cid:15) is not objective.Derivations are presented here for two different scenarios: the first assumes the most general form ofHelmholtz free energy ˆ ψ ( I k , J k , T ) and the second is based on the decoupled form ˆ ψ = ˆ ψ e ( I k ) + ˆ ψ m ( J k ) +ˆ ψ th ( T ) proposed in (2.23). In both cases terms in (cid:15) b are kept, providing a more general result than the onepresented in (2.25).i) General form of free energy ψ = ˆ ψ ( I , I , I , J , J , J , T ) Recall that the current configuration expres-sions for the magnetization and total stress in (2.21) are found by differentiating the Helmoltz free energyˆ ψ ( C , B , T ). In the case of an isotropic material ˆ ψ ( C , B , T ) = ˆ ψ ( I , I , I , J , J , J , T ) whose invariants areexpressed in terms of the right Cauchy-Green tensor C ≡ F T · F and B ≡ b · F according to (2.23).Applying the chain rule of differentiation to the expressions in (2.21), one obtains m = − ρ √ I (cid:32) ∂ ˆ ψ∂J I + ∂ ˆ ψ∂J c + ∂ ˆ ψ∂J c (cid:33) · b , σ = 2 ρ √ I (cid:20) ∂ ˆ ψ∂I c + ∂ ˆ ψ∂I (tr( c ) c − c ) + ∂ ˆ ψ∂I det( c ) I − ∂ ˆ ψ∂J bb + ∂ ˆ ψ∂J ( c · b )( c · b ) (cid:21) ++ 1 µ (cid:16) bb −
12 ( b · b ) I (cid:17) − (cid:16) mb + bm − ( b · m ) I (cid:17) , (A.1)where the left Cauchy-Green tensor c ≡ F · F T appears naturally in the constitutive relations (A.1). Thesubsequent algebra of small strain linearization is considerably simplified by noting that the invariantsinvolved can be alternatively expressed in terms of c and b as follows I = tr( c ) , I = 12 (tr( c ) − tr( c · c )) , I = det( c ) ; c ≡ F · F T ,J = b · b = (cid:107) b (cid:107) , J = b · c · b , J = b · c · b . (A.2)Expanding the expressions in (A.1) about c = I up to the first order in the small strain tensor (cid:15) ≡ (1 / ∇ u + u ∇ ), for (cid:107) (cid:15) (cid:107) (cid:28)
1, we obtain up to O ( (cid:107) (cid:15) (cid:107) ) m ≈ m ( c = I , b , T ) + ∂ m ∂ c (cid:12)(cid:12)(cid:12) c = I : (cid:15) , σ ≈ σ ( c = I , b , T ) + ∂ σ ∂ c (cid:12)(cid:12)(cid:12) c = I : (cid:15) ; c − I ≈ (cid:15) . (A.3)27fter lengthy algebraic manipulations of (A.1) and (A.3), the following expression for the magnetization m is found involving the scalar quantities ζ i ( (cid:107) b (cid:107) ) , i = 1 , · · · , m = ζ b + ζ tr( (cid:15) ) b + ζ ( b · (cid:15) · b ) b + ζ (cid:15) · b ; ζ ( (cid:107) b (cid:107) ) ≡ − ρ (cid:34) ∂ ˆ ψ∂J + ∂ ˆ ψ∂J + ∂ ˆ ψ∂J (cid:35) c = I ,ζ ( (cid:107) b (cid:107) ) ≡ − ζ ( (cid:107) b (cid:107) ) − ρ (cid:20) ∂∂I + 2 ∂∂I + ∂∂I (cid:21) (cid:34) ∂ ˆ ψ∂J + ∂ ˆ ψ∂J + ∂ ˆ ψ∂J (cid:35) c = I ,ζ ( (cid:107) b (cid:107) ) ≡ − ρ (cid:20) ∂∂J + 2 ∂∂J (cid:21) (cid:34) ∂ ˆ ψ∂J + ∂ ˆ ψ∂J + ∂ ˆ ψ∂J (cid:35) c = I ,ζ ( (cid:107) b (cid:107) ) ≡ − ρ (cid:34) ∂ ˆ ψ∂J + 2 ∂ ˆ ψ∂J (cid:35) c = I . (A.4)The corresponding small strain linearization expressions yield a total stress σ as the sum of an elastic e σ , a magnetic m σ and a magnetostrictive ms σ (involving terms of the order (cid:15) b ) component σ = e σ + m σ + ms σ ; e σ ≡ λ tr( (cid:15) ) I + 2 G (cid:15) , m σ ≡ µ (cid:20) bb −
12 ( b · b ) I (cid:21) − ζ [ bb − ( b · b ) I ] − ζ bb , ms σ ≡ Σ I + [Σ bb + ζ ( b · b ) I ]tr( (cid:15) ) + [Σ I + Σ bb + ζ ( b · b ) I ]( b · (cid:15) · b ) + Σ [( b · (cid:15) ) b + b ( (cid:15) · b )] ,λ ( (cid:107) b (cid:107) ) ≡ ρ (cid:34) ∂ ˆ ψ∂I − ∂ ˆ ψ∂I (cid:35) c = I + 4 ρ (cid:34)(cid:18) ∂∂I + 2 ∂∂I + ∂∂I (cid:19) (cid:32) ∂ ˆ ψ∂I + 2 ∂ ˆ ψ∂I + ∂ ˆ ψ∂I (cid:33)(cid:35) c = I ,G ( (cid:107) b (cid:107) ) = 2 ρ (cid:34) ∂ ˆ ψ∂I + ∂ ˆ ψ∂I (cid:35) c = I , Σ ( (cid:107) b (cid:107) ) ≡ ρ (cid:34) ∂ ˆ ψ∂I + 2 ∂ ˆ ψ∂I + ∂ ˆ ψ∂I (cid:35) c = I , Σ ( (cid:107) b (cid:107) ) ≡ − ζ ( (cid:107) b (cid:107) ) − ζ ( (cid:107) b (cid:107) ) + Σ ( (cid:107) b (cid:107) ) , Σ ( (cid:107) b (cid:107) ) = ζ ( (cid:107) b (cid:107) ) + 4 ρ (cid:34)(cid:18) ∂∂J + 2 ∂∂J (cid:19) (cid:32) ∂ ˆ ψ∂I + 2 ∂ ˆ ψ∂I + ∂ ˆ ψ∂I (cid:33)(cid:35) c = I , Σ ( (cid:107) b (cid:107) ) ≡ − ζ ( (cid:107) b (cid:107) ) + 4 ρ (cid:34) ∂ ˆ ψ∂J (cid:35) c = I , Σ ( (cid:107) b (cid:107) ) ≡ − ζ ( (cid:107) b (cid:107) ) + 4 ρ (cid:34)(cid:18) ∂∂J + 2 ∂∂J (cid:19) (cid:32) ∂ ˆ ψ∂J + 2 ∂ ˆ ψ∂J (cid:33)(cid:35) c = I , (A.5) A further simplification can be made for small strains in the expression of ζ : since | − ζ tr( (cid:15) ) b | << | − ζ b | , one has ζ ≈ − ρ (cid:104) ∂∂I + 2 ∂∂I + ∂∂I (cid:105) (cid:104) ∂ ˆ ψ∂J + ∂ ˆ ψ∂J + ∂ ˆ ψ∂J (cid:105) c = I . λ ( (cid:107) b (cid:107) )and G ( (cid:107) b (cid:107) ) plus five more scalars Σ i ( (cid:107) b (cid:107) ) , i = 0 , · · · , . This expansion proves that in a first orderapproximation in (cid:15) , the coefficients in the expressions for m and σ depend solely on || b || . The fact that λ and G – and hence the Young’s modulus E – may depend on || b || is referred to as the ∆ E effect (see e.g.Daniel and Hubert (2009)). ii) Decoupled form of the free energy ψ = ˆ ψ e ( I , I , I ) + ˆ ψ m ( , J , J , J ) + ˆ ψ th ( T ) Under the additionalhypothesis of additive decomposition for the specific free energy in (2.23), one obtains the simplification ζ = − ζ yielding from (A.4) the following expression for the magnetization mm = ζ [1 − tr( (cid:15) )] b + ζ ( b · (cid:15) · b ) b + ζ (cid:15) · b ,ζ ( (cid:107) b (cid:107) ) = − ρ (cid:34) ∂ ˆ ψ m ∂J + ∂ ˆ ψ m ∂J + ∂ ˆ ψ m ∂J (cid:35) c = I ,ζ ( (cid:107) b (cid:107) ) = − ρ (cid:20) ∂∂J + 2 ∂∂J (cid:21) (cid:34) ∂ ˆ ψ m ∂J + ∂ ˆ ψ m ∂J + ∂ ˆ ψ m ∂J (cid:35) c = I ,ζ ( (cid:107) b (cid:107) ) = − ρ (cid:34) ∂ ˆ ψ m ∂J + 2 ∂ ˆ ψ m ∂J (cid:35) c = I . (A.6)The corresponding expressions for the elastic e σ , magnetic m σ and magnetostrictive ms σ components of thetotal stress σ simplify from their corresponding counterparts in (A.5) into e σ = λ tr( (cid:15) ) I + 2 G (cid:15) , m σ = 1 µ (cid:20) bb −
12 ( b · b ) I (cid:21) − ζ [ bb − ( b · b ) I ] − ζ bb , ms σ = [Σ bb + ζ ( b · b ) I ]( b · (cid:15) · b ) + Σ [( b · (cid:15) ) b + b ( (cid:15) · b ))] , (A.7)where the scalars ζ , ζ , ζ are given in (A.6) and Σ and Σ given in (A.5) but with ˆ ψ replaced by ψ m . Inderiving (A.7) from (A.5) under the decoupling hypothesis, the pre-stress Σ and the corresponding Lam´ecoefficients λ, G are now constants independent of the magnetic field b . It is further assumed that the elasticprestress Σ = 0. Five functions of (cid:107) b (cid:107) are thus need to characterize the response of an isotropic, smallstrain, decoupled-energy, magnetoelastic material: ζ , ζ , ζ , Σ , Σ .A final remark is in order here to connect the above results to the constitutive equation in (2.25) thatneglects the magnetostrictive stress component ms σ . The reason for this simplification is that for small strains( (cid:107) (cid:15) (cid:107) (cid:28)
1) and assuming that the constants appearing in m σ and ms σ are of the same order of magnitude,one deduces that (cid:107) ms σ (cid:107) (cid:28) (cid:107) m σ (cid:107) . In the field of dielectric elastomers – a completely analogous problem where e → b , p → m , ε → µ − – similar results that neglect the coupled terms are justified under the typicalhypothesis of small strain and moderate electric field: (cid:15) = O ( ζ ), e = O ( √ ζ ), where ζ a vanishingly smallparameter (e.g. see Tian et al. (2012); Lefevre and Lopez-Pamies (2017)). The two coefficients ζ and ζ needed for the determination of m σ are related to the magnetic susceptibility χ ( (cid:107) b (cid:107) ) and magnetostrictivecoefficient Λ( (cid:107) b (cid:107) ) by: ζ ( (cid:107) b (cid:107) ) = χ ( (cid:107) b (cid:107) ) / [ µ (1 + χ ( (cid:107) b (cid:107) ))] and ζ ( (cid:107) b (cid:107) ) = − (cid:107) b (cid:107) ) / [ µ (1 + χ ( (cid:107) b (cid:107) ))]. A further simplification is possible for small strains: since terms in ζ (cid:15)bb (respectively ζ (cid:15)bb ) are negligible in frontof terms in ζ bb (respectively ζ bb ), one obtains Σ ≈ − ζ + Σ , Σ ≈ ρ (cid:104)(cid:16) ∂∂J + 2 ∂∂J (cid:17) (cid:16) ∂ ˆ ψ∂I + 2 ∂ ˆ ψ∂I + ∂ ˆ ψ∂I (cid:17)(cid:105) c = I andΣ ≈ ρ (cid:104) ∂ ˆ ψ∂J (cid:105) c = I ppendix B. Particular and homogeneous solution elastic stress fields From (3.25) and (3.26), the particular solution stress field e σ V components are e σ V rr = V −
12 1 − ν − ν (cid:18) r (cid:90) r rV d r + (2 p − r p − (cid:90) r V cs r p − d r − p + 1 r p +2 (cid:90) r r p +1 V cs d r (cid:19) e σ V θθ = νV − ν + 12 1 − ν − ν (cid:18) r (cid:90) r rV d r − (2 p − r p − (cid:90) r V cs r p − d r + 2 p + 1 r p +2 (cid:90) r r p +1 V cs d r (cid:19) e σ V rθ = 12 1 − ν − ν (cid:18) (2 p − r p − (cid:90) r V ∗ cs r p − d r + 2 p + 1 r p +2 (cid:90) r r p +1 V ∗ cs d r (cid:19) (B.1)where V ∗ cs ≡ V s cos(2Θ) − V c sin(2Θ) and the V potential components are given by (3.22).From (3.27) and (3.28), the homogeneous solution stress field e σ h components are e σ h rr = Φ + (cid:0) (2 p − p )Φ c r p − + (2 p + 2 − p )Φ c r p (cid:1) cos(2Θ)+ (cid:0) (2 p − p )Φ s r p − + (2 p + 2 − p )Φ s r p (cid:1) sin(2Θ) e σ h θθ = Φ + (cid:0) p (2 p − c r p − + (2 p + 2)(2 p + 1)Φ c r p (cid:1) cos(2Θ)+ (cid:0) p (2 p − s r p − + (2 p + 2)(2 p + 1)Φ s r p (cid:1) sin(2Θ) e σ h rθ = Φ r − (cid:0) p (2 p − s r p − + 2 p (2 p + 1)Φ s r p (cid:1) cos(2Θ)+ (cid:0) p (2 p − c r p − + 2 p (2 p + 1)Φ c r p (cid:1) sin(2Θ) (B.2)Application of the stress boundary condition in (3.29) provides the six Φ constants of integration in (B.2). Appendix C. Experimental determination of the magneto-mechanical coupling coefficient
Of all the material constants required for the constitutive model in (2.25) only the magneto-mechanicalcoupling coefficient Λ in (2.25) is not readily available and needs to be found from experiments. Its determi-nation is based here on results presented by Aydin et al. (2017) who provide analytical calculations as wellas experimental data from Rekik et al. (2014), for the uniaxial magnetostriction vs. the magnetic field forelectrical steel samples under different levels of mechanical prestress; a schematic of the setup is depicted inFigure C.1 based on the description of the typical experimental setup from Belahcen et al. (2006).
Figure C.1: Schematics of the magnetostriction setup.
30 thin plate of electrical steel is subjected to an external magnetic field b e along its axial direction,resulting in an axial magnetic field b = (1 + χ ) b (assumed uniform) inside the specimen, where χ is thematerial’s magnetic susceptibility . The plate is also subjected to an externally applied uniaxial stress σ ext e e and hence the total stress σ is the sum of the applied stress and the Maxwell stress in vacuum dueto the magnetic field b σ = e σ + m σ = σ ext e e + 1 µ (cid:20) b b −
12 ( b · b ) I (cid:21) (C.1)where the expressions for the elastic and magnetic part of the total stress are given by (2.25). The corre-sponding strain and the stress fields in the plate are assumed uniform with edge effects near the corners andedges of the plate neglected.Consequently the resulting axial strain (cid:15) is made of an elastic component σ ext /E plus a componentproportional to the square of the magnetic field strength ζ m ( b ) , where the curvature coefficient ζ m de-pends on the magnetic constants (susceptibility χ and magneto-mechanical coupling Λ). A straightforwardcalculation from (C.1) and (2.25), considering that the specimen’s lateral strain is (cid:15) = (cid:15) , gives twoindependent equations( λ + 2 G ) (cid:15) + 2 λ(cid:15) = σ ext + ( b ) µ (cid:2) − (1 + χ ) − χ ) (cid:3) , λ(cid:15) + 2( λ + G ) (cid:15) = − ( b χ ) µ , (C.2)where the Lam´e constants are given in terms of Young’s modulus E and Poisson ratio ν by G = E/ ν )and λ = νE/ (1 + ν )(1 − ν ). From (C.2) one obtains the sought relation between the axial strain, theexternal stress and the magnetic field as well as the expression for the curvature coefficient ζ m (cid:15) = σ ext E + ζ m ( b ) ; ζ m = ζ mχ + ζ m Λ , ζ mχ ≡ − ( − ν ) χ + χEµ (1 + χ ) , ζ m Λ ≡ − Λ Eµ (1 + χ ) . (C.3)In decomposing the curvature ζ m into a magnetic susceptibility ζ mχ and a magneto-mechanical ζ m Λ com-ponent we follow the approach of Daniel et al. (2003) , where the coefficients ζ mχ and ζ m Λ correspondrespectively to the magnetic susceptibility χ and the magneto-mechanical coupling Λ parts of the magneticstress m σ defined in (2.25).For the no external stress case ( σ ext = 0) the data from Aydin et al. (2017), which are based onthe approach adopted in Daniel et al. (2003), provide the same magneto-mechanical coupling curvature ζ m Λ = 2 10 − T − for the two materials analyzed. Unfortunately, the values for ν associated to thesematerials are not reported there. We assume typical values for steel: ν = 0 . E = 183GPa and a magneticsusceptibility χ = 4 10 , resulting in Λ ≈ − . × which is used in our calculations, as seen in Table 1. The materials used for holding the plate have no magnetic properties. In Daniel et al. (2003) and subsequent work by this research group by “ pure magnetostrictive ” strains they refer to thestrains due to the magneto-mechanical coupling Λ.” strains they refer to thestrains due to the magneto-mechanical coupling Λ.