WO2005106720A2

WO2005106720A2 - Method of determining a hysteresis curve and apparatus

Info

Publication number: WO2005106720A2
Application number: PCT/GB2005/001587
Authority: WO
Inventors: Hasan Ekerol
Original assignee: Dunlop Aerospace Limited
Priority date: 2004-04-30
Filing date: 2005-04-26
Publication date: 2005-11-10
Also published as: GB2413656A; WO2005106720A3; GB0409767D0

Abstract

Known techniques for modelling hysteresis characteristics of magnetic materials involve generating approximate data corresponding to an initial magnetisation curve of the magnetic materials. A disadvantage of this approach is that no account is taken of hysteresis effects. The present invention overcomes these disadvantages by generating a model of a magnetic material, comprising a first order non-linear differential equation for a first order non-linear dynamic system, the equation comprising a non-linear gain function and a non-linear lag function. A plurality of measurements of at least one parameter of the magnetic material are then made, and the first order gain and lag functions are empirically determined using the plurality of measurements. The apparatus thus generates a model of the magnetic material of improved accuracy and hence utility.

Description

METHOD OF CALCULATING A HYSTERESIS CHARACTERISTIC AND HYSTERESIS CHARACTERISTIC MODELLING APPARATUS

The present invention relates to a method of calculating a hysteresis characteristic of the type, for example, used for finite element model analysis and/or the calculation of hysteresis losses. The present invention also relates to a hysteresis characteristic modelling apparatus for use with the method of calculating the hysteresis characteristic.

In the field of electromagnetic systems, when an initially demagnetised magnetic material is placed in a magnetic field of varying strength, H, and a resulting flux density, B is measured at intervals of the magnetic field strength, a plot of the resulting flux density, B, against the field strength, H, yields a familiar BH characteristic.

In this respect, the plot reveals that the flux density increases with increasing field intensity, first at a slow rate. As the field intensity is increased, the rate of increase of the flux density becomes more rapid. Eventually, no significant increase in flux density can be attained irrespective of the further increase in the field strength, the rate of increase decreasing. At this part of the curve, the flux density increases along a straight line only very gradually with a slope equal to the permeability of free space, μ₀ . At this part of the curve, the magnetic material is magnetically saturated and, in magnetic terms, behaves essentially like free space. The part of the characteristic curve leading up to the point of saturation is known as an initial magnetisation curve, and it is this part of the characteristic that is normally used in magnetic calculations.

When the field intensity is decreased from a maximum value, H_max , gradually to zero, the flux density does not return to zero following the initial magnetisation curve, but 'lags' behind the field intensity of the initial magnetisation curve and has a definite flux density value when the field intensity reaches zero. The material therefore retains some magnetism and essentially becomes a permanent magnet. This residual magnetism is termed "residual flux density" or "remanence". The value of the remanence depends upon the type of the material, the crystal structure of the material, and the maximum field intensity value, H_max , if the value is below the saturation value of the material. If a negative field intensity of varying strength, increasing in the negative direction, is now applied to the material after the magnetic field intensity has reached zero, the residual magnetism of the material reduces to zero when the negative field intensity attains a value, H₀ , known as the "coercive force" or simply the "coercivity". If the field intensity is reversed further to a maximum negative value, - H_max , and then back to zero, the material possesses a negative remanence. When the field intensity is increased in the positive direction again to the maximum value, H_max , the familiar BH characteristic, or hysteresis loop, is completed.

The form and size of the hysteresis loop changes in response to the frequency, ω, of the applied field intensity, for example, the loop would broaden if the frequency of the alternating field intensity were to be increased.

Since the introduction of finite element methods (FEM) to solve field problems, first in solid mechanics and later in the fields of electromagnetism and fluid mechanics, the modelling of rotary or linear electrical machines and devices, such as linear actuators and torque motors encountered in electromechanical hydraulic equipment, has become possible with increasing accuracy. In more recent years, commercially available FEM software packages have become available and are effective design tools to assist in the design of such devices. These packages are now able to calculate flux distributions accurately in such devices in two- or three dimensions, especially for calculations where the magnetic field intensity and the magnetic flux density are time-invariant.

However, where the magnetic field intensity and flux density vary with time, for example during transient or AC modelling of such devices, the hysteresis loop for the device being modelled is also required so that instantaneous permeability can be inferred from a

hysteresis loop given an instantaneous local field intensity vector, H . Some modelling software is able to perform transient and AC calculations, but using very approximate techniques.

Further, as can be seen from a hysteresis loop, energy imparted into the material during the application of the magnetic field is not totally recoverable, since the process known as "magnetic domain changes" undergone by the material is non-reversible. This irreversibility of the process results in so-termed "hysteresis loss" that presents itself in the form of heat. Hysteresis losses together with eddy current losses are an important issue in electrical machine design and therefore there is a strong desire to estimate these at the design stage as accurately as possible.

Also, known electrical machine design practice requires hysteresis losses to be calculated; usually, one of two methods is employed. In the first method, the inner loop area of the hysteresis loop for a given design, corresponding to the applied field intensity and the operating frequency of the applied field, is assumed to be proportional to the frequency and the full DC hysteresis loop area, i.e. Loss = k.f.(DC - loopArea), where k is a constant of proportionality and f is the frequency.

Alternatively, in the second method, the inner loop for the given design is assumed to be proportional to the applied frequency, f, and the peak flux density, B_P, used in the design. Peak flux density for the design is first estimated and the losses are inferred from loss data provided by magnetic material manufacturers in graphical form, the data normally being plotted in terms of loss vs. peak flux density at various experimental frequencies. Given the design flux density and the frequency, the loss can then be interpolated from these plots.

However, both of the above methods lack generality and involve extensive empiricism. Furthermore, in modern design practices involving FEM methods, a distribution of flux density needs to be calculated. Consequently, calculation of total loss, i.e. the combined losses of eddy current losses and hysteresis losses, could be a complicated and lengthy exercise, with erroneous results in any event, because correct account cannot be taken of the hysteresis effects in these models.

Clearly, without proper representation of the hysteresis loop for the material, these above losses cannot be correctly estimated and the transient response of a device correctly assessed. Hence, using current BH characteristic models, accurate assessment of the transient and AC performance of electromagnetic devices is not possible. Additionally, it is practically impossible to generate a good estimate of eddy current losses associated with electrical devices.

Although the magnetisation phenomena can be explained on the basis of what may be termed 'descriptive physics', which makes use of concepts such as magnetic domains mentioned above, domain wall motion, electron spin and magnetic moments, quantitative analysis is not possible by this route, because an analytical treatment of the physics of the process currently seems unattainable and the construction of a model for the BH characteristics of magnetic materials from first physical principles is also very difficult at present.

According to a first aspect of the present invention, there is provided a method of estimating a hysteresis characteristic of a physical system, the method comprising the steps of: modelling the hysteresis characteristic using an equation for a first order dynamic system; wherein the equation comprises: a non-linear gain function and a non-linear lag function; making a plurality of measurements of at least one parameter of the physical system, and empirically determining the first order gain and lag functions using the plurality of measurements.

Preferably, the physical system is an electromagnetic field system or a lumped parameter system.

Preferably, the physical system is a magnetic material.

The hysteresis characteristic may be a magnetic flux density-magnetic field intensity (BH) characteristic.

Preferably, the first order dynamic system is modelling using a differential equation. More preferably, the differential equation is: τ(H,B,ω) . A + B = G(H,B,ω) . H(t) at where τ(H,B,ω) is the lag function, and G(H,B,ω) is the gain function in general vectorial form.

Preferably, the gain function is determined by: measuring plurality of amplitude permeabilities at respective different frequencies; and calculating a plurality of respective ratios of the plurality of amplitude permeabilities to a DC permeability. More preferably, each of the respective ratios of the plurality of amplitude permeabilities to permeabilities of zero or substantially zero frequency is a gain factor. For example, DC permeabilities or permeabilities at frequencies close to 0 Hz, such as 1 to 5 Hertz, most preferably about 1 Hz. Very preferably, the gain factor (K_G) for a given frequency (ω) is calculated using:

where B is a peak magnetic flux density and H is a magnetic field intensity.

Preferably, the gain function is:

G(H,B,ω) = K_G ^DC = [l -cωⁿ exp[-(||H| - m|)^p /(SLO)] )^DC H where c, s, m, p and n are material dependent constants. However, the gain function may correspond to a frequency close to DC, such as 1 to 5 Hz.

The gain function may be provided by a look-up table. The look-up table may comprise a first dynamic gain corresponding to a first frequency and a second dynamic gain corresponding to a second frequency, wherein a dynamic gain corresponding to an intermediate frequency with respect to the first and second frequencies is determined by an interpolation technique.

Preferably, a value of the gain function for a predetermined, non-sinusoidal, field intensity is determined by the steps of: determining frequency and amplitude components of the predetermined field intensity; and referencing the frequency and amplitude components in the look-up table to construct a composite gain function. More preferably, the frequency and amplitude components are determined using Fourier analysis of the predetermined field intensity.

Preferably, the lag function is:

^■ B τ(H,B,ω) = ΔB _\ω_,t ΔH ^ΔH ΔΓ

ΛH ΛH

Preferably, a sinusoidal input field is provided having a magnitude of ( ) , wherein ( — ) =

ωH is provided. More preferably, a first set of peak flux density magnetisation data is used to B B calculate (— )^ω , and a second set of magnetisation data is used to calculate (TT)"'⁴ and

ΛB

( )^ω,t . The second set of magnetisation data may be initial magnetisation data. ΔH Preferably, the lag function is:

wherein a, b and q are empirical constants.

The equation for the physical system may have a vectorial form.

According to 'a second aspect of the present invention, there is provided a method of estimating a hysteresis loss, the method comprising the method of estimating the hysteresis characteristic as set forth in accordance with the first aspect of the present invention.

The method may further comprise the step of measuring an area within a hysteresis loop described by the hysteresis characteristic at a frequency of interest.

According to a third aspect of the present invention, there is provided a method of estimating an eddy current loss, the method comprising the method of estimating the hysteresis characteristic as set forth above in relation to the first aspect of the invention.

According to a fourth aspect of the present invention, there is provided a computer program element comprising computer program code means to make a computer execute the method as set forth above in relation to the first aspect of the present invention.

Preferably, the computer program element is embodied on a computer readable medium.

According to a fifth aspect of the present invention, there is provided a hysteresis characteristic modelling apparatus comprising: a processing unit arranged to generate a model of a hysteresis characteristic of a physical system, the model of the hysteresis characteristic comprising: an equation for a first order dynamic system, the equation including a first order, non-linear, gain function and a, non-linear, lag function; wherein the first order gain and lag functions are empirically determined functions based upon a plurality of measurements of at least one parameter of the physical system.

Preferably, the physical system is an electromagnetic field system. The physical system may be a ferromagnetic material. The hysteresis characteristic may be a magnetic flux density-magnetic field intensity (BH) characteristic.

The apparatus may further comprise a differential equation to represent the first order dynamic system.

Preferably, the differential equation is: HR τ(H,B,ω) • ( -) + B = G(R,B,ω) • H(t) dt where τ(H,B,ω) is the lag function, and G(H,B,ω) is the gain function.

The gain function may be:

G(H,B,ω)

where K_G is a gain factor, and c, s, m, p and n are material dependent constants. However, the gain function may correspond to a frequency close to DC, such as 1 to 5 Hz.

The apparatus may further comprise a look-up table constituting the gain function. The apparatus may employ a curve fitting algorithm in order to determine the gain function and/or the lag function. The curve fitting algorithm may be used to refine constants associated with the gain function and/or the lag function.

The look-up table may comprise a first dynamic gain corresponding to a first frequency and a second dynamic gain corresponding to a second frequency, wherein the processing means is arranged to calculate a dynamic gain corresponding to an intermediate frequency with respect to the first and second frequencies using an interpolation technique.

An impulse field intensity can be decomposed into frequency and amplitude components, for example, Fourier components. The look-up table may be used to obtain a composite gain function for each harmonic of the field intensity and frequency of the pulse. Alternatively or additionally, the composite gain function may be obtained by substituting components of the field intensity impulse into the algebraic gain function determined. Preferably, a value of the gain function for a predetermined, non-sinusoidal, field intensity are frequency and amplitude components of the predetermined field intensity referenced from the look-up table.

Preferably, the lag function is:

Preferably, the apparatus further comprises a sinusoidal input field having a magnitude of

(^r) . wherein (— ) = ωH.

More preferably, a first set of peak flux density magnetisation data is used to calculate B B ΔB

(— -)^ω , and a second set of magnetisation data is used to calculate (77) and (" 7τ)^ω' ■

The second set of magnetisation data may be initial magnetisation data.

Preferably, the lag function is:

wherein a, b and q are empirical constants.

Preferably, the equation for the system has a vectorial form.

Preferably, the processing unit is further arranged to estimate a hysteresis loss associated with the physical system.

Preferably, the processing unit is further arranged to estimate an eddy current loss associated with the physical system.

According to a sixth aspect of the invention, there is provided an apparatus for generating a finite element model, the apparatus comprising the hysteresis characteristic modelling apparatus as set forth above in relation to the fifth aspect of the invention. It is thus possible to provide a method of modelling a hysteresis characteristic with improved accuracy. It is also therefore possible to calculate hysteresis losses and/or eddy current losses with greater accuracy than with existing techniques. Further, it is possible to model electromechanical machines and/or devices with greater accuracy where time-invariant field intensities are employed.

At least one embodiment of the invention will now be described, by way of example only, with reference to the accompanying drawings, in which:

Figure 1 is a schematic diagram of an apparatus constituting an embodiment of the invention; Figure 2 is a plot of a hysteresis loop for a hypothetical linear ferromagnetic material; Figure 3 is a plot of a magnetisation curve at different frequencies; Figure 4 is a gain factor plot; Figure 5 is a plot of a predicted hysteresis loop for a sample material; Figure 6 is a flow diagram of a method for the apparatus of Figure 1 ; Figure 7 is a flow diagram of a set of steps for determining a gain function and not shown in Figure 6; Figure 8 is a flow diagram of a set of steps for determining a lag function and not shown in Figure 6; and Figure 9 is a flow diagram of a set of alternative steps to those of Figure 7 for determining the gain function and obtaining values for integration thereof.

Referring to Figure 1 , a hysteresis characteristic calculation apparatus comprises a computing device, such as a Personal Computer (PC) 100, having a processing unit 102 coupled to a data bus 104, the data bus 104 being coupled to a non-volatile memory 106, for example a Read Only Memory (ROM), and a volatile memory 108, for example a Random Access Memory (RAM). An input device 110 and an output device 112 are coupled to the_. processing unit 102, as well as a data storage device, such as a hard drive 114 via the data bus 104. In this example, the hard drive 114 stores computer program code to be executed by the PC 100 in order to model a hysteresis characteristic of a magnetic material, for example a ferromagnetic material. However, it should be appreciated that not only finite element systems can be modelled. Other systems, such as so-called lumped parameter systems can also be modelled. The use of the programmed PC 100 is optional, and the skilled person will appreciate that other in hardware implementations are possible. In order to better understand the operation of the program, the background theory upon which this practical implementation will also be included herein.

Consider a system to be modelled comprising a slender toroid having a large mean-diameter to cross-section ratio, produced from a hypothetical linear ferromagnetic material of a constant permeability, μ(H). If a single turn of a conductive sheet of cross-section dA is wrapped around the toroid, and a version of Ampere's law that includes a displacement current is applied, the following equation results:

d dL = J+ (^)]. dA (1 )

C A "t

where: dL is the infinitesimal length along the central axis of the toroid placed at the mean diameter, and n is the unit vector normal to the surface.

The left hand side of equation (1 ) can be integrated along the toroid axis and the right hand side can be integrated over the conductor cross-section; the displacement current can be ignored in the conductor, but can be included in relation to the material. Also, the field intensity in the cross-section of the slender toroid can be assumed to be constant, and so equation (1) reduces to:

(2) L V dV ^l

In equation (2) above, dA_t is the cross-sectional area of the part of the material contained by the circular toroid axis.

Using the substitutions: D = εE and J = σE , and also assuming that at every instant in time, B = μH holds true, the following relationship linking the field intensity, H and the flux density, B can be obtained:

(!)(^) + B = μH (3) σ dt where: ε is the dielectric constant, σ is the conductivity of the material, ε (— ) has the dimensions of time; and σ μ is permeability.

Although equation (3) lacks generality, equation (3) suggests that the BH characteristics of a hypothetical ferromagnetic material with constant permeability, μ can be represented by a first order dynamic system governed by a differential equation of the form:

where: T is a time constant of the system, and μ can be regarded as a static gain of the system.

A sinusoidal field intensity of H = H₀ sin(ωt) is then applied to the system, and equation (4) is integrated with respect to time over a time period corresponding to the 1.25 of a cycle of the applied field. The results of the integration using the above expression for field intensity were plotted (Figure 2) as flux density, B vs. field strength, H. The ellipse described corresponds to a hysteresis loop of the hypothetical linear magnetic material of constant permeability.

Using Figure 2, a phase angle (lag), φ between the input, H and the output, B can be calculated from the expression:

φ = sirf¹(-) (5) a where: a is the ordinate of the maximum field intensity on the hysteresis loop, and b is the distance from the origin to the point where the ellipse intersects the flux density axis.

From Figure 2, it can be seen and/or concluded that the maximum flux density on the loop lies on the gain line of the plot. Also, the tip of the loop does not lie on the gain line, but rather deviates from the gain line; the higher the frequency of the field, the further the tip of the loop deviates from the gain line. Further, at only very low field intensity frequencies, when the phase angle is small, the tip of the loop lies in close proximity to the gain line, and as the frequency of the field intensity is increased, the maximum flux density remains on, but "slides" down, the gain line. It should be noted that if the frequency of the field intensity increases beyond a value ω₀ = 1/ τ , the dynamic gain, B_p / H (instantaneous peak permeability), of the system decreases rapidly. As the system is a linear dynamic system, the frequency value ω₀is considered the "bandwidth" of the system.

In relation to the loop described in Figure 1 , an inner loop can be obtained by reducing the magnitude of the field intensity, H₀, and the loop broadens when the frequency of the field intensity is increased.

In the light of the above observations, it can be seen that a first order lag equation produces all the characteristics that are observed in the hysteresis loop of real magnetic materials. However, the permeability of real ferromagnetic materials is not constant, but a function of the input field intensity, H, at a given frequency. The above initial model therefore has to be refined to model real ferromagnetic materials.

Referring to Figure 6, since magnetic material manufacturers provide data for magnetisation curves obtained (Step 600) at various frequencies, plotted in terms of a peak flux density, B_P and an applied field intensity, H, this data can be used to assist in refining the above initial model for a given material, such as a sample of such data for a material known as "HCR" (also known as Helcomax) available from Carpenter Technology (UK), produced in 0.1 mm strips. Alternatively, a characteristic curve (Figure 3) for the HCR material was plotted by applying a low amplitude fixed frequency AC field to the material, and then the amplitude of the field was increased in steps until saturation of the material was reached (Step 600). If data for this step is not available, a so-called permeameter can be interfaced with the PC 100 in order to make necessary measurements. For each step, the peak value of the flux density was measured and the process was then repeated at other frequencies for the field. The curves of Figure 3 can be regarded as constituting the gain curves of the material when regarded as a highly non-linear dynamic system, because the peak flux density, B_P, always lies on these curves. Furthermore, since different curves are obtained at different frequencies, the material is considered to approximate to different dynamic systems for different respective field intensity frequencies.

Consequently, when comparing this behaviour with the model for hypothetical constant permeability materials, it can be seen that real magnetic materials exhibit a non-linear gain characteristic that is not only a function of applied field (variable permeability), but also a function of the frequency of the applied field. Thus, the BH characteristic of ferromagnetic materials can be modelled as highly non-linear first order systems using first order, nonlinear, lag equations. Such an equation can be expressed in vector notation as follows:

- - - dR - - - τ(fi,B,ω) . -) + B = G(H,B,ω) . fl(t) (6) dt

where: τ(H,B,ω) is a non-linear lag function that is a function of instantaneous field, fi , and flux density, B , vectors, and the frequency of the applied field, ω, and

G(H,B,ω) is the non-linear gain function, also a function of the field intensity, fi , and the flux density, B , vectors and the frequency of the applied field, ω.

Equation (6) constitutes a dynamic vector field describing the dynamic relationship between

the magnetic field vector, fi , and the flux density vector, B , and is also applicable to non- isotropic media if the components of the lag and the gain functions can be identified.

In order to determine the non-linear gain and lag functions of equation (6), semi-empirical techniques are, in this example, employed.

Typically, an initial magnetisation curve of a ferromagnetic material has an S-shaped characteristic curve, although the individual form of the characteristic curve may differ significantly depending on the type of the material.

Turning to Figure 7, the peak flux density B_P vs. field intensity data collected is provided to the PC 100 as a data file and subsequently loaded. In the present example, data for frequencies of 50, 400, 1600, 3200Hz were provided by the material supplier, but the DC magnetisation curve can not be provided for the obvious reason that they are unattainable, since the applied field frequency would have to be reduced to zero. Instead, a hypothetical DC characteristic was generated (Step 700) by the PC 100 by extrapolating data to zero frequency, so as to be able to determine the DC gain characteristics at zero frequency.

In order to determine the gain function to be used as a dynamic gain curve at any desired frequency, the amplitude permeabilities of the material obtained from the peak flux density distributions, and measured at different frequencies, and their ratios to the DC permeability were calculated by the PC 100 (Step 702). The algebraic expression for the ratio of the amplitude permeabilities to the DC permeability is:

where: K_G is a gain factor.

The gain factor K_G was then plotted (Figure 4) by the PC 100 against the field intensity, H for different frequencies, ω, where it could be seen that the curves had inverted Gaussian- shaped characteristics. The PC 100 then implemented a curve fitting algorithm (Step 704) , trying a number of known equations describing respectively differently curves until, in this example, the following approximate empirical function constituting equation (6) was found:

,B_C

G(H,B,ω) = K_Gφ^DC = [l - cωⁿ exp[-(|H| - m|)^p /(sω)]]( )^DC (8) H H where: c, s, m, p and n are material-dependent constants that can be determined empirically from experimental data, and K_G is the gain factor represented by the term within the square brackets.

Turning to the lag function, τ(H,B,ω) (Figure 8), a BH loop associated with equation (6), dB having the empirical gain function, [G(H,B,ω)]^ω , has a slope ( — )^ω at frequency, ω, at a dt particular instant in time. Consequently, equation (6) can be solved for the lag function using experimental magnetisation data:

ΔH For a sinusoidal AC input field, the magnitude of ( ) is equal to ωH and so equation (9)

above becomes: ( P \ω _ ( \ω,t τ(H,B,ω) = -ϋ- ^H (10) ,ΔB _V ω( — ) ,ω,t

If the data for the full hysteresis loop and the initial magnetisation curve (curve obtained during the application of the first quarter of the cycle with the material at its demagnetised state) for the material is available in addition to the magnetisation curves of Figure 3(the peak flux density distribution), equation (10) above could be plotted, the denominator being obtained through differentiation of the initial magnetisation curve. Consequently, the nature of the lag function can be determined for the particular material at the particular frequency ω at which the experimental data was obtained. In order to properly determine the nature of the lag function, calculations should be made in respect of sets of data corresponding to at least 3 different frequencies, respectively.

Alternatively (Figure 8), in order to determine the lag function of equation (10), two consecutive magnetisation curves were used, by the PC 100, from the magnetisation data obtained or measured initially (in Step 600). The first curve of the two curves was used to D calculate (— — )^ω (Step 800), whilst the second magnetisation curve was used to H B ΔB calculate (— )^ω,t (Step 802) and ( )^ω,t (Step 804), because the initial magnetisation curves H ΔH of the magnetisation curves used were expected to be similar in form to those obtained from peak flux density measurements.

R R ΛR

The values of (^-)^ω , (-)^ω,t , and ( )^ω,t were then substituted by the PC 100 into equation H H ΔH

(10) and the lag function of equation (10) was plotted (Step 806) for different values of field intensity, H. Thereafter, a function was then selected, by the curve fitting algorithm, to best reflect the curves:

where: (— ) is the instantaneous dynamic amplitude permeability, and H a, b and q are empirical constants. The curve function selected by the curve fitting algorithm was then modified by the algorithm choosing (Step 808) an optimum value for the factor b/ω^q. The curve fitting algorithm then selects (Step 810), an, initially, small tentative value for the empirical constant 'a', and then refines estimates of 'a' until an optimum value is obtained. Using the empirical constants determined, the lag function can then be formed (Step 810).

It should be noted at this point that the lag function is considered to be not only a function of frequency of applied field intensity, but also of the magnetic state of the material, i.e. the p instantaneous dynamic permeability, μ' = (— ). H

By substituting (step 602) the expressions derived above for the gain and lag functions into equation (6), we have an expression for BH characteristics for an isotropic material, or in a suitable form for modelling by a lumped parameter technique: )^DCH(t) (12)

Equation (12) has been expressed in terms of magnitudes of vectors, since the vectors have the same direction in the case of the isotropic material.

For anisotropic media, equation (12) above can be expressed in the following vectorial form, provided components of the magnetisation curves were available for all directions required as well as data to calculate and the empirical constants a, b, c, m, n, p, q and s associated with the components:

where 1 is a three component unit vector.

Whilst the above example derivations of the gain and lag functions have been described in the context of semi-empirical techniques, the skilled person will appreciate that other forms of algebraic functions can also be derived. However, it should be noted that the lag function should be considered to be a function of not only the instantaneous magnetic state of the material, i.e. the instantaneous permeability, μ' = (B /H) , but also of the frequency of the applied field, ω; the lagging of the flux density behind the field intensity is affected by the magnetic state of the material and hence the dependency on the instantaneous permeability, μ' = (B / H) . This is also true of the gain function.

In order to obtain BH characteristics of the material used in this example, the PC 100 integrates equation (12) to obtain an expression for flux density, B, in terms of field intensity, H. In this example, a numerical integration technique (Step 604) was employed by the PC 100, such as, a 4^th Order Runga-Kutta algorithm. Thereafter, the integrated equation (12) was plotted (Figure 5) for frequencies of 1, 5, 50, 400 and 1600 Hz (Step 606).

The curve fitting algorithm was then executed again in order to refine (Step 608) the calculated values for the constants a and b in the lag function, by using the magnetisation data initially obtained in Step 600.

It should be appreciated that the above model for the gain function is not unique and that other forms of algebraic expressions, such as a series, exist. Further, whilst the above technique for deriving the above gain function has been described, it should be appreciated that other suitable known techniques can be employed.

For example (Figure 9), if extensive magnetisation data is available, steps 602, 604, 700, 702, 704 can be replaced by use of a data file generated (Step 900) so as to represent the gain function, but as a look-up table rather than as a function. The data file would therefore contain a numerical gain function, G(H,B,ω) = (B_p /H)^ω , tabulated against the field intensity,

H and the frequency, ω. During the integration of equation (11), the look-up table would be used in conjunction with a suitable two-dimensional interpolation technique (Step 902) in order to determine the dynamic gain for intermediate frequencies not present in the look-up table, prior to integration (Step 904). This technique is particularly useful where algebraic functions cannot be found to represent the gain function.

Where arbitrary input field intensities in non-sinusoidal form are used, for example impulse input field intensities, the look-up table is also useful. In such situations, the input field intensity can be decomposed into Fourier components for frequency and amplitude, and then referenced in the look-up table, thereby constructing a composite gain function.

Alternatively, if insufficient data is available, the impulse field intensities can nevertheless be decomposed into Fourier components, and the harmonic amplitudes and frequencies used in conjunction with equation (8), the algebraic gain function, in order to represent the composite gain function.

After determining the composite gain function by either of the above techniques, the composite gain function can be used, in respect of the frequency to which the composite gain function corresponds, for any field intensity amplitude when solving Maxwell's equations in the manner to be described later herein. However, for a different fundamental harmonic frequency, i.e. a main frequency of the field, the composite gain function needs to be recalculated and the coefficient 'b' of the lag function needs to be scaled down by a ratio of the first harmonic amplitude to the amplitude of the harmonic being used to recalculate the composite gain function.

Having determined a suitable mathematical model for the BH characteristics of a magnetic material, it is desirable to use the calculation apparatus to apply the above model in the calculation of hysteresis losses.

In AC magnetisation of magnetic materials, the area bounded by the hysteresis loop traced during a given cycle gives the hysteresis loss per unit volume of the core per cycle.

In order to calculate the hysteresis loss, the PC 100 is arranged to calculate instantaneous

local field intensity and the flux density vectors, fi and B , using Maxwell's equations and use the calculated vectors in conjunction with equation (13) described above.

Assuming each component of the field intensity vector, fi , imparts a certain amount of

energy into the material proportional to the magnitude of the field intensity vector, fi , each pair of field intensity/flux density vectors would trace out a separate hysteresis loop. The sum of areas of each loop described by each pair of BH vectors yields a total loss per cycle per unit volume for every cycle. For example, an expression for the area enclosed by the

loop traced out by one pair of x -components of H and B vectors is:

A_x = §dH_xdB_x (14) BH

Hence, for 3-dimensional space, the total area is the sum total of the loss during one cycle per unit volume of the material:

A = §dH_xdB_x +§ dH_ydB_y + dH₂dB_z (15) A_χ A_y A_z

Or, in vectorial form:

Direct application of this integral is, in practice, difficult owing to the fact that direct application requires the existence of a priori knowledge of the BH characteristic contour, i.e. data for a complete cycle after a first quarter cycle is required. Instead, in this example, the integral of equation (16) is reduced to a line integral that can be more easily handled by simultaneous integration with one of the field equations, for example equation (12) or equation (13), over one cycle following the first quarter cycle.

Transformation of the integral of equation (16) can be achieved by first considering each individual integral separately, i.e. as in equation (15), and applying Green's theorem, more particularly a special case of Green's theorem, to these integrals, thereby transforming each integral into a respective line integral: fjdxdy = c[(ydx + 2xdy) (17) A C

Using the special case of Greens theorem above, equation (15) can be re-written as: A_BH = §dfi.dB = r|(B_xdH_x +2H_xdB_x)+ cf(B_ydH_y +2H_ydB_y)+ ^(B_zdH_z + 2H_zdB₂) A BH_X BH_y BH₂ (18) Or, in vectorial form:

A_BH = φdH »dB = < (B « dfi+ 2H« dB) (19) BH BH

For simplicity and speed of processing, the PC 100 integrates the component integrals of equation (18) separately using any suitable numerical integration technique instead of performing the integral of equation (19). However, to make integration possible, Maxwell's equations are solved concurrently to determine the spatial distribution of the field intensity with time and the previously determined BH model of equation (11) is used to obtain corresponding values of flux density. Field intensity and flux density values are calculated over a time period corresponding to one and a quarter cycles.

For the purposes of clarity and ease of understanding, a simple calculation of hysteresis loss using an isotropic material will now be described. However, it should be appreciated that hysteresis losses can be calculated for more complex materials.

Using a Graphical User Interface (GUI) provided by the PC 100, a user of the PC 100 can select an option to indicate to the PC 100 that a material being analysed is an isotropic

material. For isotropic materials, flux density and field intensity vectors fi , B have a same direction and so hysteresis loss is given by the following expression:

Loss = c[(BdH + 2HdB) (20) BH where B and H are magnitude of associated flux density and field intensity vectors, respectively. The PC 100, as explained above, simultaneously integrates equation (20) and equation (12), whilst assuming a sinusoidal field intensity, to yield: a + <^~ -cωⁿ exp[-(||H| - m|)^p /(sω)]]fe^DCH ω^M )( Α + = [l ₀ sin(ωt) (21) dt n ' i H

In relation to equation (20), integration is carried out after the first quarter cycle, i.e. the first quarter cycle is skipped, and the integration process is started from a point where the field intensity is at a maximum value, the integration being carried out over one cycle after the first quarter cycle. The PC 100 then calculates the hysteresis loss as a function of field intensity for the full loop, and then, if desired, as a function of the peak flux density, B_P, at different applied frequencies.

Some or all of the above functionality can be integrated or interfaced with a suitable FEM software package or a lumped parameter system modelling package to obtain improved modelling of devices, for example, electrical machines, such as rotary or linear electrical machines.

Alternative embodiments of the invention can be implemented as a computer program product for use with a computer system, the computer program product being, for example, a series of computer instructions stored on a tangible data recording medium, such as a diskette, CD-ROM, ROM, or fixed disk, or embodied in a computer data signal, the signal being transmitted over a tangible medium or a wireless medium, for example, microwave or infrared. The series of computer instructions can constitute all or part of the functionality described above, and can also be stored in any memory device, volatile or non-volatile, such as semiconductor, magnetic, optical or other memory device.

Claims

1. A method of estimating a hysteresis characteristic of a physical system, the method comprising the steps of: modelling the hysteresis characteristic using an equation for a first order dynamic system; wherein the equation comprises: a first order, non-linear, gain function and a first order, non-linear, lag function; making a plurality of measurements of at least one parameter of the physical system, and empirically determining the first order gain and lag functions using the plurality of measurements.

2. A method as claimed in Claim 1 , wherein the physical system is an electromagnetic field system or a lumped parameter system.

3. A method as claimed in Claim 1 or Claim 2, wherein the physical system is a magnetic material.

4. A method as claimed in any one of Claims 1 , 2 or 3, wherein the hysteresis characteristic is a magnetic flux density-magnetic field intensity (BH) characteristic.

5. A method as claimed in any one of the preceding claims, wherein the first order dynamic system is modelling using a differential equation.

6. A method as claimed in Claim 5, wherein the differential equation is: τ(fi,B,ω) • (— ) + B = G(H,B,ω) • H(t) dt where τ(H,B,ω) is the lag function, and G(H,B,ω) is the gain function in general vectorial form.

7. A method as claimed in Claim 1 , wherein the gain function is determined by: measuring plurality of amplitude permeabilities at respective different frequencies; and calculating a plurality of respective ratios of the plurality of amplitude permeabilities to a DC permeability.

8. A method as claimed in Claim 7, wherein each of the respective ratios of the plurality of amplitude permeabilities to permeabilities of zero or substantially zero frequency is a gain factor.

9. A method as claimed in Claim 8, wherein the gain factor (K_G) for a given frequency (ω) is calculated using:

where B is a magnetic flux density and H is a magnetic field intensity.

10. A method as claimed in any one of the preceding claims, wherein the gain factor is:

G(H,B, ω) = K_G ^DC = [1 - cωⁿ exp[-(||H| - m|)^p /(sω)]] DC where c, s, m, p and n are material dependent constants.

11. A method as claimed in any one of the preceding claims, wherein the gain function is provided by a look-up table.

12. A method as claimed in Claim 11 , wherein the look-up table comprises a first dynamic gain corresponding to a first frequency and a second dynamic gain corresponding to a second frequency, wherein a dynamic gain corresponding to an intermediate frequency with respect to the first and second frequencies is determined by an interpolation technique.

13. A method as claimed in any one of the preceding claims, wherein a value of the gain function for a predetermined, non-sinusoidal, field intensity is determined by the steps of: determining frequency and amplitude components of the predetermined field intensity; and referencing the frequency and amplitude components in the look-up table to construct a composite gain function.

14. A method as claimed in Claim 13, wherein the frequency and amplitude components are determined using Fourier analysis of the predetermined field intensity.

15. A method as claimed in any one of the preceding claims, wherein the lag function is:

16. A method as claimed in Claim 15, further comprising a sinusoidal input field having a magnitude of ( ) , wherein ( ) = ωH. ΔΓ ΔΓ

17. A method as claimed in Claim 16, wherein a first set of peak flux density p magnetisation data is used to calculate (— -)^ω , and a second set of magnetisation

R ΛR data is used to calculate (— )^ω'⁴ and ( )^ω,t .

18. A method as claimed in any one of the preceding claims, wherein the lag function is:

wherein a, b and q are empirical constants.

19. A method as claimed in any one of the preceding claims, wherein the equation for the physical system has a vectorial form.

20. A method of estimating a hysteresis loss, the method comprising the method of estimating the hysteresis characteristic as claimed in any one of the preceding claims.

21. A method as claimed in Claim 20, further comprising the step of measuring an area within a hysteresis loop described by the hysteresis characteristic at a frequency of interest.

22. A method of estimating an eddy current loss, the method comprising the method of estimating the hysteresis characteristic as claimed in any one of Claims 1 to 19.

23. A computer program element comprising computer program code means to make a computer execute the method as claimed in any one of the preceding claims.

24. A computer program element as claimed in Claim 23, embodied on a computer readable medium.

25. A hysteresis characteristic modelling apparatus comprising: a processing unit arranged to generate a model of a hysteresis characteristic of a physical system, the model of the hysteresis characteristic comprising: an equation for a first order dynamic system, the equation including a first order, nonlinear, gain function and a, non-linear, lag function; wherein the first order gain and lag functions are empirically determined functions based upon a plurality of measurements of at least one parameter of the physical system.

26. An apparatus as claimed in Claim 25, wherein the physical system is an electromagnetic field system.

27. An apparatus as claimed in Claim 25 or Claim 26, wherein the physical system is a ferromagnetic material.

28. An apparatus as claimed in any one of Claims 25, 26 or 27, wherein the hysteresis characteristic is a magnetic flux density-magnetic field intensity (BH) characteristic.

29. An apparatus as claimed in any one of Claims 25 to 28, further comprising a differential equation to represent the first order dynamic system.

30. An apparatus as claimed in Claim 29, wherein the differential equation is: - dR - τ(H.B,ω) . Α + B = G(H,B,ω) • H(t) dt where τ(H,B,ω) is the lag function, and G(H,B,ω) is the gain function.

31. An apparatus as claimed in any one of Claims 25 to 30, wherein the gain function is:

G(H,B,ω) = K_G ^DC = [1 - cωⁿ exp[-(||H| - m|)^p /(sω)]] ^DC where K_G is a gain factor, and c, s, m, p and n are material dependent constants.

32. An apparatus as claimed in any one of the Claims 25 to 31 , further comprising a lookup table constituting the gain function.

33. An apparatus as claimed in Claim 32, wherein the look-up table comprises a first dynamic gain corresponding to a first frequency and a second dynamic gain corresponding to a second frequency, wherein the processing means is arranged to calculate a dynamic gain corresponding to an intermediate frequency with respect to the first and second frequencies using an interpolation technique.

34. An apparatus as claimed in any one of Claims 25 to 33, wherein a value of the gain function for a predetermined, non-sinusoidal, field intensity are frequency and amplitude components of the predetermined field intensity referenced from the lookup table.

35. An apparatus as claimed in any one of Claims 25 to 34, wherein the lag function is: (-__JL)^ωH - B τ(H.B,ω) = -Ji ΔB _t ΔH ΔH ^ Δt

36. An apparatus as claimed in Claim 35, further comprising a sinusoidal input field ΛH ΛH having a magnitude of ( ) , wherein ( ) = ωH.

37. An apparatus as claimed in Claim 36, wherein a first set of peak flux density B_p magnetisation data is used to calculate {^~τr , and a second set of magnetisation

R ΔB data is used to calculate (-)"^•' and ( )^ω'' . ΔH^;

38. An apparatus as claimed in any one of Claims 25 to 37, wherein the lag function is:

wherein a, b and q are empirical constants

39. An apparatus as claimed in any one of Claims 25 to 38, wherein the equation for the system has a vectorial form.

40. An apparatus as claimed in any one of Claims 25 to 39, wherein the processing unit is further arranged to estimate a hysteresis loss associated with the physical system.

41. An apparatus as claimed in any one of Claims 25 to 40, wherein the processing unit is further arranged to estimate an eddy current loss associated with the physical system.

42. An apparatus for generating a finite element model, the apparatus comprising the hysteresis characteristic modelling apparatus as claimed in any one of Claims 25 to 41.