COMSOL 6.4 - Vector Hysteresis Modeling

This model reproduces the TEAM (Testing Electromagnetic Analysis Method) problem 32, which aims to evaluate numerical methods for the simulation of anisotropic magnetic hysteresis. A hysteretic three-limbed laminated iron core is subject to a time-varying magnetic field generated by two coils. The Jiles–Atherton material model (available in the Magnetic Fields interface) is used to simulate the response of the material, reproducing published experimental and numerical data.

Model Definition

The geometry of the simulated experimental setup of TEAM problem 32 is represented in Figure 1.

Figure 1: The geometry of the device. The coils are colored in blue and the core in red.

The system is composed of a three-limbed magnetic core with two feeding coils on the two outer limbs. A low excitation frequency (10 Hz) and a finite lamination of the frame prevent skin effects in the core. The frame is composed of 5 layers with a thickness of 0.48 mm.

The applied magnetic field is mainly oriented in the xy-plane; the material is anisotropic and react differently to fields applied along the x or the y direction. In the experimental setup there are a series of pick-up coils used to accurately probe the magnetic field; these coils are not included in the model as point measurements can easily be made by means of direct numerical evaluations.

Ref. 1 details four analysis cases that differ in the applied excitations. The case represented numerically in this model is the third one, in which the two coils are excited with an AC source with a peak value of 14.5 V and in quadrature phase. The coils have a total DC resistance of 11.42 Ω which includes an externally applied resistance. The field generated by this setup is strong enough to drive the material to saturation, while the phase shift creates a rotating field at the junction between the central limb and the frame.

In the literature, experimental results have been compared favorably with vector hysteresis models (Ref. 1 and 4). This model follows Ref. 4 in using the empirical Jiles-Atherton magnetic hysteresis model to simulate the core material. The values of the parameters for the empirical model are presented in Table 1. For an anisotropic material, the parameters are all diagonal matrices; the table reports the values on the diagonal.

Table 1: Parameter for the Jiles-Atherton model.
Parameter	Symbol	Values on the diagonal
Saturation magnetization	Ms	1.31e6 A/m, 1.33e6 A/m, 1.31e6 A/m
Domain wall density	a	233.78 A/m, 172.856 A/m, 233.78 A/m
Pinning loss	k	374.975 A/m, 232.652 A/m, 374.975 A/m
Magnetization reversibility	c	736e-3, 652e-3, 736e-3
Inter-domain coupling	alpha	562e-6, 417e-6, 562e-6

The Jiles-Atherton model is particularly suitable for AC feeding and requires only a limited number of parameters: a and Ms control the slope of the hysteretic B-H curve respectively at zero field and at saturation; c and k control the strength of the hysteretic effects — with the limit of no hysteresis for c = 1 or for large k. The values presented in Table 1are taken from Ref. 1 and are obtained by fitting the model to experimental data.

Results and Discussion

Figure 2 shows the magnetic flux at two different time instants, t = 275 ms (top) and t = 300 ms (bottom), at which the current in respectively the left and the right coil is at the peak value. The images show how the magnetic field rotates in the xy-plane at the junction between the central limb and the outer frame.

Figure 2: Magnetic flux density at t = 275 ms (top) and t = 300 ms (bottom).

The hysteretic behavior can be displayed by plotting the magnetic flux density as a function of the magnetic field during one AC cycle (corresponding to one hysteresis loop). Figure 3 shows the hysteresis loop obtained by averaging the quantities on a cross section of the central limb.

Figure 3: Hysteresis (B-H) loop in the central limb.

Finally, a representation of instantaneous magnetization field is shown in Figure 4. In those figures the red and blue vectors represent the instantaneous fields respectively at the time t = 300 ms and t = 275 ms, when the fields are expected to be in quadrature. The plots highlight how the fields at the junction are rotating in the xy-plane.

Figure 4: Magnetization vector field at t = 275 ms (red) and t = 300 ms (blue).

The spatial distribution of hysteresis losses can also be calculated, as illustrated in Figure 5. Example hysteresis loops at different locations are presented in Figure 6.

For further details, see the modeling instructions below.

Figure 5: Time-averaged loss distribution in the core.

Figure 6: Two different hysteresis loops for two different locations.

Notes About the COMSOL Implementation

The application uses the Jiles-Atherton hysteresis model available in the physics interface. The anisotropic material is constructed starting from the default isotropic Jiles-Atherton material (available in the material library) and modifying the properties appropriately.

To obtain a good compromise between an accurate solution, robust convergence and efficient solving, the following settings are used:

•

A direct solver () is used instead of the default iterative solver. In order to solve a problem with a direct solver it is necessary to apply the feature.

•

The discretization order for the magnetic vector potential A is set to use elements. The discretization order of the Jiles-Atherton auxiliary dependent variables is then automatically set to zero.

•

The scales of the dependent variables are set manually, in order to take advantage of the information on the maximum expected value of the magnetic field and the magnetization in the hysteretic material.

References

1. www.compumag.org/wp/

2. http://www.cadema.polito.it/team32

3. A.J. Bergqvist, “A Simple Vector Generalization of the Jiles-Atherton Model of Hysteresis,” IEEE Transactions on Magnetics, vol. 32, no. 5, p. 4213, 1996.

4. J.P.A. Bastos and N. Sadowski, Magnetic Materials and 3D Finite Element Modeling, CRC Press 2014.

5. S. Yan and J.-M. Jin, “Theoretical Formulation of a Time-Domain Finite Element Method for Nonlinear Magnetic Problems in Three Dimensions,” Progress In Electromagnetics Research, vol. 153, pp. 33–55, 2015.

ACDC_Module/Verifications/vector_hysteresis_modeling

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > .

Click .

Click

In the tree, select > .

Click

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
W	174.5[mm]	0.1745 m	Width of the core
H	180[mm]	0.18 m	Height of the core
w	30[mm]	0.03 m	Width of the central limb
h1	H-2*w	0.12 m	Height of the windows
w1	(W-3*w)/2	0.04225 m	Width of the windows
Th	5*0.48[mm]	0.0024 m	Thickness of the core
f	10[Hz]	10 Hz	Feeding voltage frequency
R_coil	11.42[ohm]	11.42 Ω	Coil resistance
p	1/f	0.1 s	Time period

Definitions

Step 1 (step1)

In the toolbar, click

and choose .

Variables 1

In the window, right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Unit	Description
Q_core	(mf.Hxd(mf.Bx,t))+(mf.Hyd(mf.By,t))+(mf.Hz*d(mf.Bz,t))	W/m³	Core hysteretic loss density
P_core	intop1((mf.Hxd(mf.Bx,t))+(mf.Hyd(mf.By,t))+(mf.Hz*d(mf.Bz,t)))		Core hysteretic loss
P_coils	mf.ICoil_1^2mf.RCoil_1 + mf.ICoil_2^2mf.RCoil_2		Coil resistive losses
P_rem	intop2((mf.Hxd(mf.Bx,t))+(mf.Hyd(mf.By,t))+(mf.Hz*d(mf.Bz,t)))		Remaining Power
P_input	mf.PCoil_1 + mf.PCoil_2		Input power

Integration 1 (intop1)

In the toolbar, click

and choose .

Step 1 (step1)

In the window, click .

In the window for , locate the section.

In the text field, type 0.5.

Click to expand the section. In the text field, type 1.

Geometry 1

Work Plane 1 (wp1)

In the toolbar, click

Work Plane 1 (wp1) > Plane Geometry

In the window, click .

Work Plane 1 (wp1) > Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type W.

In the text field, type H/2.

Locate the section. From the list, choose .

In the text field, type 3*H/4.

Work Plane 1 (wp1) > Rectangle 2 (r2)

In the toolbar, click

In the window for , locate the section.

In the text field, type w1.

In the text field, type h1/2.

Locate the section. From the list, choose .

In the text field, type -(w+w1)/2.

In the text field, type H-h1/4.

Work Plane 1 (wp1) > Rectangle 3 (r3)

In the toolbar, click

In the window for , locate the section.

In the text field, type w1.

In the text field, type h1/2.

Locate the section. From the list, choose .

In the text field, type (w+w1)/2.

In the text field, type H-h1/4.

Work Plane 1 (wp1) > Difference 1 (dif1)

In the toolbar, click

and choose .

Select the object only.

In the window for , locate the section.

Click to select the

toggle button for .

Select the objects and only.

Extrude 1 (ext1)

In the window, under > right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Distances (m)
Th/2

Cylinder 1 (cyl1)

In the toolbar, click

In the window for , locate the section.

In the text field, type w*0.7.

In the text field, type h1/2.

Locate the section. In the text field, type -w-w1.

In the text field, type H-h1/2.

Locate the section. From the list, choose .

Click to expand the section. In the table, enter the following settings:


Layer name	Thickness (m)
Layer 1	0.1*w

Delete Entities 1 (del1)

In the window, right-click and choose .

In the window for , locate the section.

From the list, choose .

Click the

button in the toolbar.

On the object , select Domains 3–5 only.

Copy 1 (copy1)

In the toolbar, click

and choose .

Select the object only.

In the window for , locate the section.

In the text field, type (w+w1)*2.

Block 1 (blk1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 2*W.

In the text field, type H.

In the text field, type 3*w.

Locate the section. In the text field, type -W.

In the toolbar, click

Click the

button in the toolbar.

Click the

button in the toolbar to get a better view.

The volume around the components is modeled using the feature. This applies a small but nonzero conductivity in 3D simulations to obtain consistent equations.

Magnetic Fields (mf)

Apply a on the antisymmetry cut to set the appropriate boundary condition (zero tangential magnetic field). The default is the correct boundary condition for the symmetry cut boundaries (zero normal magnetic field).

Symmetry Plane 1

In the toolbar, click

and choose .

In the window, click .

Select Boundaries 5, 9, 11, 16, 21, 30, 37, 42, 44, and 49 only.

In the window for , locate the section.

From the list, choose (all the boundaries at x = 0).

Ampère’s Law in Solids 1

In the toolbar, click

and choose .

Select Domain 3 only (the core).

It might be easier to select the correct domain by using the window. To open this window, in the toolbar click and choose . (If you are running the cross-platform desktop, you find in the main menu.)

In the window for , locate the section.

From the list, choose .

Domain Coil 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

In the Vcoil text field, type 14.5[V]*sin(2*pi*f*t)*step1(f*t).

Locate the section. In the N text field, type 90.

From the list, choose .

In the Rcoil text field, type R_coil.

Select Domains 5 and 6 only.

Geometry Analysis 1

In the window, click .

In the window for , click to expand the section.

In the FL text field, type 2.

In the FA text field, type 2.

Input 1

In the window, expand the node, then click .

Select Boundary 51 only.

Geometry Analysis 1

In the window, click .

Output 1

In the toolbar, click

and choose .

Select Boundary 35 only.

Domain Coil 2

In the toolbar, click

and choose .

Select Domains 2 and 4 only.

In the window for , locate the section.

From the list, choose .

In the Vcoil text field, type 14.5[V]*cos(2*pi*f*t)*step1(f*t).

Locate the section. In the N text field, type 90.

From the list, choose .

In the Rcoil text field, type R_coil.

Geometry Analysis 1

In the window, click .

In the window for , locate the section.

In the FL text field, type 2.

In the FA text field, type 2.

Input 1

In the window, expand the node, then click .

Select Boundary 26 only.

Geometry Analysis 1

In the window, click .

Output 1

In the toolbar, click

and choose .

Select Boundary 6 only.

Apply a feature to improve the stability of the computation and in order to use a direct solver.

Gauge Fixing for A-Field 1

In the toolbar, click

and choose .

Using lower order shape functions improves the robustness of the solution process for a nonlinear problem such as the Jiles–Atherton hysteresis model. Lowering the order also reduces the size of the problem, making it easier to solve with and a direct solver.

In the window, click .

In the window for , click to expand the section.

From the list, choose .

Materials

Coil

In the window, under right-click and choose .

Select Domains 2 and 4–6 only.

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Relative permeability	mur_iso ; murii = mur_iso, murij = 0	1	1	Basic
Relative permittivity	epsilonr_iso ; epsilonrii = epsilonr_iso, epsilonrij = 0	1	1	Basic

In the text field, type Coil.

The following steps create the material for the Jiles–Atherton hysteresis model. First add the isotropic default material, then modify it to make it anisotropic.

Add Material

In the toolbar, click

to open the window.

Go to the window.

In the tree, select > .

Right-click and choose .

In the toolbar, click

to close the window.

Materials

Jiles–Atherton Hysteretic Material (mat2)

Select Domain 3 only.

Start by specifying the basic material properties. Since the – will be used for the magnetic behavior of the material it is not necessary to specify a magnetic permeability.

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Electric conductivity	sigma_iso ; sigmaii = sigma_iso, sigmaij = 0	1	S/m	Basic
Relative permittivity	epsilonr_iso ; epsilonrii = epsilonr_iso, epsilonrij = 0	1	1	Basic

Now update the parameters specific to the Jiles–Atherton model. For each parameter in the table, perform the following steps:

In the window, expand the node, then click .

Click the corresponding row.

Click the button below the table.

Choose and enter the diagonal elements according to the following table:


Parameter	Values on the diagonal
Saturation magnetization	1.31e6[A/m], 1.33e6[A/m], 1.31e6[A/m]
Domain wall density	233.78[A/m], 172.856[A/m], 233.78[A/m]
Pinning loss	374.975[A/m], 232.652[A/m], 374.975[A/m]
Magnetization reversibility	736e-3, 652e-3, 736e-3
Interdomain coupling	562e-6, 417e-6, 562e-6