COMPUTATIONAL AND MEASUREMENT METHODS FOR FERROMAGNETIC HYSTERESIS IN DIFFERING STEEL ALLOYS by KYLE ANDREW TOWNSEND B.S., University of Colorado Denver , 2016 A thesis submitted to the Faculty of the Graduate School of the University of Colorado in pa rtial fulfillment o f the requirements for the degree of Master of Science Electrical Engineering Program 2018
ii This thesis for the Master of Science degree by Kyle Andrew Townsend has been approved for the Electrical Engineering Program by Stephe n D. Gedney, Chair Mark Golkowski Vijay Harid Date: May 12, 2018
iii Townsend, Kyle Andrew ( M.S. , Electrical Engineering Program) Computational and Measurement Methods for Ferromagnetic Hysteresis in Differing Steel Alloys Thesis directed by Professor Stephe n D. Gedney ABSTRACT The objective of this thesis is to compute the magnetization induced in a non linear ferromagnetic material when exposed to slowly time varying magnetic fields. This thesis also explored the further development of an empirical non lin ear hysteretic model for magnetic differential susceptibility using the Schneider  cooperative hysteresis model. For the magneto stati c Volume Integral Equation (VIE) is used in tandem with a non linear Stepped Solver. Experimental methods have been developed to extract the fundamental physical parameters needed for the cooperative hysteresis model, as well as to validate the non linear LCN solution for the VIE. The form and content of this abstract are approved. I recommend its publication. Approved: Stephen D. Gedney
iv ACKNOWLEDGEMENTS I would like to express my deep and profound gratitude to Dr. Stephen Gedney, who has been my resea rch advisor and thesis advisor since my time under his tutelage began almost three years ago. His insight, patience, and thorough knowledge of electromagnetics has been invaluable in the continuation of my education. I would also like to thank Dr. Mark Gol kowski and Dr. Vijay Harid for serving as committee members on my thesis panel. I would also like to acknowledge the Office of Naval Research, who funded the CU Denver Magnetics Research Lab, and also made my degree possible. I would also like to thank my wife Ashley, who has shown continuous faith in my pursuit of higher education. Without her support, I would not have had the opportunity to pursue my degree. I am eternally indebted to her for her never ending support. Finally, I would like to dedicate my thesis to my children, Joule and Ari, who I hope can look to this project and realize that anything can be achieved with hard work and dedication.
v TABLE OF CONTENTS CHAPTER I. INTRODUCTION ................................ ................................ ................................ ...... 1 II. MAGSTR M ................................ ................................ ................................ ............. 3 Quasi magnetostatic Volume Integral Equation ................................ ........................ 3 Numerical Integration ................................ ................................ ................................ . 6 Mid Point Rule ................................ ................................ ................................ .... 6 Orthogonal Polynomials ................................ ................................ ...................... 8 Gauss Quadrature Rule ................................ ................................ ..................... 10 Curvilinear Coordinate Systems ................................ ................................ ............... 11 Differential Surface and Volume Elements Using Curvilinear Coordinates ............ 15 Differential Operators ................................ ................................ ............................... 17 Locally Corrected Nystr m (LCN) Method ................................ ............................. 20 Non linear Transient Solver ................................ ................................ ..................... 24 III. MATERIAL PARAMETERS AND DIFFERENTIAL SUSCEPTIBILITY ........... 28 Material Parameters ................................ ................................ ................................ .. 28 Demagnetization ................................ ................................ ................................ 28 Initial Susceptibility ................................ ................................ .......................... 32 Magnetic Saturation ................................ ................................ .......................... 35 Coercive Field ................................ ................................ ................................ ... 35
vi Cooperative Differential Susceptibility Model ................................ ........................ 36 Exponen tial Cooperative Differential Susceptibility Model ................................ .... 39 Parallel Differential Susceptibility Model ................................ ................................ 40 IV. MEASUREMENT METHODS ................................ ................................ ............... 42 Vibrating Sample Magnetometer ................................ ................................ ............. 42 Faraday Method ................................ ................................ ................................ ........ 44 V. VALIDATION ................................ ................................ ................................ ......... 47 Cooperative Model ................................ ................................ ............................ 47 Exponential Model ................................ ................................ ............................ 53 Comparison of Cooperative and Exponential Model ................................ ............... 55 Future Work ................................ ................................ ................................ .............. 57 LIST OF REFFERENCES ................................ ................................ ................................ ...... 60
1 CHAPTER I INTRODUCTION Full wave electromagnetic ( EM ) solvers are an indispensable tool for calculating electromagnetic field scattering effects for a variety of problems, ranging from simple to extremely complex. This thesis focuses on a custom E novel approach to calculate the non linear magnetization induced in ferromagnetic materials equation . With full wave EM solvers, ther e is almost always a need to control the approximation accuracy. The algorithm developed by Young and Gedney  allows for high fidelity models of complex structures and exhibits exponential convergence of the solution . With M for bulk magnetization, one can calculate the magnetic signature of an arbitrarily shaped ferromagnetic sample to high accuracy. This thesis also focuses on non linear magnetic hysteresis models developed by C. Schneider  . ferromagnetic material s, like initial susceptibility, saturation magnetization, and coercive field . The cooperative model accurately represents differential susceptibility over a broad range of magnetic field intensities . By implementing this non linear model in the EM software, one can calculate hysteretic internal magnetization, magnetic fields, and the radiation of fields beyond the sample itself. Finally, with the ability to calculate the magnetic moment of various steel types and differing geometries, verification of the models can be done. The CU Denver Magnetics Research Lab is equipped with custom hardware capable of measuring the magnetic moment
2 of th in steel disks and the induced magnetization of hollow steel pipes. This allows the comparison of numerical models to the physical ly measured data .
3 CHAPTER II MAGSTR M Quasi magnetostatic Volume Integral Equatio n Many numerical models exists for calculating electromagnetic (EM) problems for a variety of materials. For magnetostatic analysis , Young and Gedney  introduce d a L ocally C (LCN) formulation to solve the quasi magnetostatic volume integral equation ( VIE ) . Because of their inhomogeneous, anisotropic, and often nonlinear properties, VIE methods are ideal for magnetostatic modeli ng. I n order to solve the VIE, we must first create a volumetric mesh of our device. The VIE solution method is capable of handling a variety of polyhedra , such as tetrahedrons, prisms, and hexahedrons. One other benefit is that when only ferrous materials need to be meshed, and the surrounding non magnetic medium, including air, do not  . This thesis will focus solely on hexahedrons for the mesh elements, but the curvilinear coordinate system to be described can be extended t o handle other volume elements as well . First, consider a magnetic material with an inhomogeneous permeability Âµ( r ) occupying a volume region V , bound by surface S , with surface normal directed outward from V . Let the volume be placed inside of a magnetic field that is generated externally. With V inside of the magnetic field, induced magnetization will produce a demagnetizing field As such, the total magnetic field is represented as ( 2 . 1 ) where is the external applied field, and is the demagnetizing field. Assuming there are no externally applied currents on surface S, then the curl of is 0, or . As a
4 result , th e internal field can be represented as the gradient of some scalar potential We then show the magnetic flux density as ( 2 . 2 ) where Âµ o is permeability of free space, and is the magnetization vector in such that ( 2 . 3 ) (r) is the magnetic susceptibility which can be expressed as a function of the relative permeability as ( 2 . 4 ) G dictates that Thus, taking the divergence of equation ( 2 . 2 ) , t he left hand side of the equation becomes 0 . Also, since is solenoidal . Consequently, leads to ( 2 . 5 ) By using the vector identity ( 2 . 6 ) we can express the divergence of as ( 2 . 7 ) where is an equivalent magnetic charge density inside volume The expression is in the form e m, the solution for can be expressed as ( 2 . 8 )
5 where , the distance between the source point and observer point , tion , and is an equivalent magnetic surface charge density . By substituting ( 2 . 7 ) into ( 2 . 8 ) , and utilizing the fact that the discontinuity in the magnetization projected onto the outward normal vector over the entire surface is the equivalent magnetic surface current density ( 2 . 9 ) The volumetric integral term in ( 2 . 9 ) can be rewritten using the vector identity ( 2 . 10 ) C ombining ( 2 . 9 ) and ( 2 . 10 ) , the internal field can be rewritten as ( 2 . 11 ) Therefore, and can again be rewritten as ( 2 . 12 ) Lastly, by combining ( 2 . 1 ) ( 2 . 3 ) , ( 2 . 11 ) , and ( 2 . 12 ) , an equation for the H M VIE can be expressed as ( 2 . 13 ) w here ( 2 . 14 )
6 This derivation was followed largely from Young and Gedney  , and while other methods of derivation exist, this method is chosen for the calculations done in this thesis . Numerical Integration Mid Point Rule When calcul ating volume integrals for arbitrary geometries, the closed form solution is often unavailable or far too complex to be solved. Rather, one can approximate the solution to the integral using various numerical integration techniques to a high order of accur acy. Consider the integral T he integral I can be approximated via a weighted sum ( 2 . 15 ) where are the abscissa points in the range It is assumed that is finite and c ontinuous along the int erval T he lowest order approximation to the integral is mid point approximation . To this end, the interval is broken into M subsections . T he abscissae are then assigned at the mid point between x i and x i+1 , and sets the weight to be the width of the interval. Consequently , I can be rewritten as ( 2 . 16 )
7 For this case , the error converges as where The visualization of equation ( 2 . 16 ) can be seen in Figure 1 , where the arbitrary function shown in red is f(x) and the abscissae are chosen at the center points between x i and x i+1 . Figure 1 Visualizati on of mid point quadra ture rule . For example, let a smooth polynomial be ( 2 . 17 ) We can solve ( 2 . 17 ) analytically to find the exact solution, and examine the error convergence of the midpoint rule. First, the analytical solu tion to ( 2 . 17 ) is ( 2 . 18 ) Let the function be bounded by , which leads to an exact solution of . Next, by applying the midpoint quadrature rule in equation ( 2 . 16 ) , we can examine the error of the convergence by modifying . This error is shown in Table 1 .
8 As stated, using the unique case where the abscissa is selected at the midpoint of , the error in the integral solution converged as . This is shown as the number of abscissae is doubled, or equivalently i s halved, the error percentage is one fourth that of the previous width . While this method has its own applications, it will be shown that an Table 1 Evaluation of error when approximating polynomial function using midpoint quadrature rule. Number of Abscissae Approximation Value Error % Error 1 3 11.125 0.1250 1.1111 0.5 6 11.2188 0.0313 0.2778 0. 2 5 12 11.2422 0.0078 0. 0694 0. 1 25 24 11.248 0.0020 0.0174 Orthogonal Polynomials Orthogonal Poly nomials are specific classes of polynomials defined over the range such that ( 2 . 19 ) where is the weighting function, and is the Kronecker delt a function. This relationship is also known as the inner product between the two functions and I f the two polynomials are orthogonal, then we know that the inner product is zero when . A
9 special case of orthogonal polynomials are the Legendre Polynomials in which the weighting function is u nity and has the range . Legendre polynomial, , can be defined b y the recursive relation in ( 2 . 20 ) . Table 2 shows the actual polynomials evaluated up to N=4. ( 2 . 20 ) Table 2 Evaluation of Legendre polynomials up to 4th order j= 1 j=0 j=1 j =2 j =3 By definition, an n th order polynomial will have at most n distinct real roots. For the roots lie within the range If one chooses to use the n roots of the polynomial as the abscissae for an N point quadrature rule, then solving for the weights leads to a quadrature rule of degre e 2 N 1 and can be utilized for high order convergence when using the Gauss Quadrature r ule.
10 Gauss Quadrature Rule Gauss Quadrature rules are a set of quadrature rules that provide optimal convergence rates. The optimal convergence rate can be utilized by letting both the sampling points and respective weights be unknown s , hence leading to 2 N degrees of freedom. With 2 N degrees of freedom, one can expect that the error will converge exponentially as  . Starting with equation ( 2 . 15 ) , we can rewrite the weights as a weighting function, and limit the weights so that the quadrature rule is satisfied for all polynomials up to N 1. ( 2 . 21 ) Now, utilizing the fact that are orthogonal polynomials, the linear system of equations in ( 2 . 22 ) results. ( 2 . 22 ) where are the abscissae, or roots, of the N th polynomial. On the right side, we have that all terms e xcept the constant polynomial term are 0 because of the orthogonality of the polynomials developed in the previous section. It can be shown that each weight has the closed form solution  ( 2 . 23 )
11 where and is the fir st order derivative of the N th order Legendre polynomial. The use of the Gauss Quadrature r ule yields high order accuracy using numerical integration, which is what we are looking for when solving the right hand side of equation ( 2 . 13 ) . Curvilinear Coordinate Systems General meshing consists of fitted polyhedron that provide a discrete representation of an arbitrary geometry . General curvilinear polyhedron, suc h as hexahedron, prism, and tetrahedron, can be expressed via curvilinear coordinates. Further, the curved edges and convex surfaces can be represented by higher order, curvilinear polyhedron . Figure 2 Illustration of a genera l curvilinear coordinate system. The vectors are known as the unitary vectors. A curvilinear coordinate system is one in which the 3, not necessarily orthogonal, dimensions can represent each point P uniquely. Figure 2 shows curvilinear surfaces represented when one dimension of the curvilinear space is held constant. Also in this
12 coordinate system, one can describe a line by letting two of the coordinate values be held constant. As will be shown, use of the curvilinear coordinate allows for the rapid calculation of the volume element used in the volume integral of ( 2 . 13 ) For a given region, let be continuous, single valued functions  with each being a function of the Cartesian coordinates x, y, and z , ( 2 . 24 ) In general, are continuously differentiable, but occasionally it may contain singular points in which case care must be taken to avoid app lication of the general formulas. These functions are defined as general or curvilinear coordinates  . Now, let point be represented now as , whi ch is uniquely defined by ( 2 . 24 ) . We can now create the vector r which represents the vector of point P from an arbitrary origin in the given region , ( 2 . 25 ) Now a differential change in the vector r is represented by ( 2 . 26 ) If the differential change in r is along the continuous curve then that differential change is tangential to curve and the vectors known as unitary vectors can be defined as ( 2 . 27 )
13 The unitary vectors are not necessarily of unit length, but are rather a base system of reference for all vectors associated with point The l ength of the vectors are dependent on the nature of the general coordinates We can combine equations ( 2 . 26 ) and ( 2 . 27 ) to yield a compact form ( 2 . 28 ) The unitary vectors are used to define a differential volume V ( 2 . 29 ) Using this volume calculated as the triple scalar product , 3 new vectors can be defined as the reciprocal unitary vectors in ( 2 . 30 ) which are orthogonal to dual pair s of unitary vectors defined in ( 2 . 27 ) . The reciprocal vectors are defined as ( 2 . 30 ) The reciprocal unitary and unitary vectors satisfy the orthogonality relationship ( 2 . 31 ) where is the Kronecker delta function. The dif ferential length vector can then be expressed in terms of the reciprocal vectors and coordinates as ( 2 . 32 ) A common practice is to represent the scalar product of the reciprocal and unitary vectors by the metric tensor coefficients
14 ( 2 . 33 ) A general vector can be defined in the curvilinear coordinate system as the sum of either the unitary ve ctor or reciprocal unitary vector components, and can be represented as ( 2 . 34 ) where can be referred to as the contravariant component of F , and is known as the covariant component of F . Stratton  shows that using this derivation for curvilin ear coordinate systems, one can describe the geometrical properties of space completely with the metric coefficients or its reciprocal counterpart. We can then calculate a differential volume by ( 2 . 35 ) where ( 2 . 36 ) T he differential surface vector to a sur face where is constant can be defined as  ( 2 . 37 ) where is the surface Jacobian and the surface normal v ector . Consequently, to calculate the integral of vector through the surface , the integral becomes ( 2 . 38 ) Lastly, a volume integral can then be described as
15 ( 2 . 39 ) where is defined in equation ( 2 . 36 ) . Differential Surface and Volume Elements Using Curvilinear Coordinates Consider the p th order Lagrange type quadrilateral element in Figure 3 . Let the nodes p , equally spaced in the curvi linear coordinate frame to yield nodes in total. For example, a bilinear quadrilateral is defined by such that the element contains 4 total nodes. The discrete position vectors of the nodes are where and . Because the nodes are spaced equidistant along each axis, the coordinate in the local curvilinear space is located at . Next, a p th order one dimensional interpolation polynomial is introduced in the curvilinear coordinate space defined as  ( 2 . 40 ) I n the range , at the points , and and are Sylvester polynomials, defined as ( 2 . 41 ) Next, a two dimensional interpolation function space can be defined by the product of two one dimensional functions defined in ( 2 . 40 ) as ( 2 . 42 ) T he position vector of a point lying on the quadrilate ral patch is given by
16 ( 2 . 43 ) Figure 3 Cur vilinear quadrilateral surface defined with equally spaced nodes along and axes. With the position vector fully defined within the local curvilinear coordinate frame, the unitary vectors can be computed from ( 2 . 27 ) as  ( 2 . 44 ) By combining equations ( 2 . 40 ) ( 2 . 42 ) , one can solve the partial derivative in the right hand side of ( 2 . 44 ) as ( 2 . 45 )
17 where is the derivative of with res pect to its argument. A similar expression can be derived for T o complete the three dimensional space for the surface , a coordinate that is orthogonal to and should be used so the space can be completely defined. Per the definition in ( 2 . 37 ) , let the third unitary vector be defined as the unit normal ( 2 . 46 ) The same principles described here can be extended to the three dimensional hexahedral, triangul ar prism, pyramid, and tetrahedral elements. Differential Operators With the derivation of curvilinear coordinates that can be used for a polyhedral cell, derivation of the gradient and divergence will be paramount in the application of the Locally Correc The gradient of a scalar function will be a fixed vector defined in both magnitude and direction as the maximum rate of change with respect to the coordinates of To begin, we examine the change in incurred during displacement This change can be expressed as ( 2 . 47 ) where are the contravariant vector components of By combing ( 2 . 28 ) and ( 2 . 47 ) , the equality becomes
18 ( 2 . 48 ) The differential displacement is arbitrary, so ( 2 . 48 ) can now be written as  ( 2 . 49 ) To note, the expression for the gradient is expressed using the contravariant coordinate system, but can easily be expressed in the covariant coordinate system using the equality ( 2 . 50 ) which allows the transformation directly from covariant to contravariant coordinate systems. With the gradient fully defined, also necessary in the derivation of the LCN method is the divergence operator on a vector function The volume shown in Figure 4 should be used as reference. Consider the left and right su rfaces located at and The area of the face located at is expressed as The normal vector of the face is directed outwards from the volume, or to the left. The area of the face located at can be calculated in a similar fashion, except with the modification so that the normal is facing outward from the volume as well, or to the right. These two terms c an be summed in order to find the net contribution of the outward flux from the volume as ( 2 . 51 ) where the subscript after the bracket indicates the flux from each surface. For sufficiently small values of ( 2 . 51 ) can be approximated by using the linear term of a Taylor expansion
19 Figure 4 Volume element located in a curvilinear coordinate system  . ( 2 . 52 ) where is replaced with Now, the term can be rewritten as ( 2 . 53 ) by the relations in ( 2 . 31 ) , ( 2 . 34 ) , and ( 2 . 35 ) . Thus , the contribution from the surfaces can be expressed as ( 2 . 54 ) Similar derivations can be done for the remaining 4 surfaces of the volume. Lastly, we are measuring the flux contribution from the surfaces per unit volume and
20 take the limit such that ensuring that only the linear terms in the Taylor expansion remain. With this, we can now express the divergence of as ( 2 . 55 ) Locally Corrected Nystr m (LCN) Method T he magnetostatic VIE, is solved in this thesis via the Locally Corrected (LCN) method an efficient and simple technique for the discretization of the integral equations with regular kernels  . If the kernel does contain singularities, the must handle the singularities and near singularities that arise from arbitrary geometries. Let a volume be represented by a discretized mesh composed of curvilinear polyhedr a , such as hexahedron, prisms, and/or tetrahedron. A numerical quadrature rule, such as the Gauss Quadrature rule discussed previously is introduced over each individual cell allowing the volume integral in ( 2 . 13 ) to be approximated as  ( 2 . 56 ) where is the total number of volume elements in V , is the number of quadrature points on the i th cell for the k th vector component of are the quadrature abscissae, are the weights, and is the cell Jacobian defined in ( 2 . 36 ) and is evaluated at .
21 discretization , the unknown components are the magnetization vectors evaluated at every quadrature point in the cell. In order to solve for the unknowns, the discrete integral operator in ( 2 . 56 ) is sampled with a test vector at the quadrature points. The dis crete current is represented by three independent vectors, superimposed upon one another. These three independent vectors are chosen to be the set of unitary vectors of the local cell geometry to expand the magnetic current . Thus , the magnetic current along the k th unitary vector can be described as ( 2 . 57 ) where is the unknown coefficient, is the k th unitary vector, and is the cell Jacobian defined in ( 2 . 36 ) . Next, the set of test vectors must be selected in order to preserve the diagonal of the matrix. The reciprocal unitary vectors are a chosen because of the orthogonal relation of the reciprocal and unitary vectors. ( 2 . 58 ) The test vector sampled at the abscissa point is ( 2 . 59 ) With the definition of the source current along the k th un itary axis, ( 2 . 14 ) can now be expressed more accurately as  ( 2 . 60 ) w here and .
22 Given the choice of test vectors ( 2 . 59 ) and source vectors ( 2 . 57 ) , the diagonal term of the discrete integral operator in ( 2 . 56 ) can be represented as ( 2 . 61 ) As one may notice, ( 2 . 60 ) shows the possibility of a singularity occurring when , or the source node and field node are coincident , and , such that Also, in the very near region when R is very small, the quadrature rule converges slowly . This eliminates the high order con vergence of the quadrature rule. This can be mitigated by apply ing a local correction in the near region to account for the singular behavior when source and field points are coincident or very close in proximity. Locally corrected quadrature weights can t hen be computed specifically for each of the source cell and field point pairs within the small domain where is slowly converging. In light of this singularity, Young and Gedney  pose the followin g formulation to account for the singularities and near singularities. In short, a specialized quadrature rule is computed for the singularity in the integrand of ( 2 . 56 ) . Although the discretization of volume V can contain curvilinear prisms, tetrahedron, as well as hexahedron, the focus of  is on hexahedra. In the local coordinates of the volume cell (hexahedra), the magnetizati on vector components along the unitary vector are expanded using a mixed order basis function sets where ( 2 . 62 ) and
23 ( 2 . 63 ) as well as ( 2 . 64 ) are the local cur vilinear coordinates of the cell, and are the m th order Legendre polynomials. When the basis function sets described in ( 2 . 62 ) are used, is complete to order along coordinate and complete to order p along the other two curvilinear axes. Because of this property, we can represent the magnetic charge densit y as a linear combination of the terms  ( 2 . 65 ) which is complete to polynomial order p. Next, when the field is evaluated at test point , the moments are set equal for the locally corrected quadrature rule on the i th cell for the k th vector component, or ( 2 . 66 ) where are the corrected quadrature weights that ar e to be determined. As the test point and test vector change, both and the right hand side of ( 2 . 66 ) change. The application of ( 2 . 66 ) to the basis function set leads to a new linear system of equations  ( 2 . 67 )
24 is local to the k th vector component of on the i th cell , and is computed by solving the right hand side of ( 2 . 66 ) using adaptive quadrature , and the singularity extractio n methods discussed in  . Non linear Transient Solver Young and Gedney, et. al.  , developed an algorithm for a stepped nonlinear solution method capable for quasi magnetostatic analysis of mag netic materials, which are inherently non linear. The advantage of this algorithm is that it is not limited to the LCN method, and can be implemented using a number of models that support magnetic hysteresis. For materials with nonlinear properties that ar e nonhysteretic, other solvers are capable of converging quickly to the final solution. The problem with these methods is that they are not capable of handling hysteresis in materials. T hus the non linear transient solver developed in  does handle non linear hystere tic magnetization implicitly and is sufficient for the preferred analysis. T he LCN discretization for the quasi magnetostatic VIE presented in the previous section is used herein, however should be again noted that it is not required for the solver. To begin, let the LCN discretization be expressed in a matrix operator form as  ( 2 . 68 ) where is a n LCN source vector at the quadrature points, is a n vector of the unknown magnetization current densities at the quadrature points, and is a dense matrix representing the discretization of the integral operator in ( 2 . 13 ) . The magnetic susceptibility tensor is described as ( 2 . 69 )
25 The arises from the operation Methods already exist for rapid conversion of solutions when solving for nonlinear, nonhysteretic materials  ,  which work well when working with normal susceptibilities. Hysteretic materials however are more easily solved by using the differential susceptibility tensor which satisf ies ( 2 . 70 ) where the differentials are approximated by the small change s represented as Now, let k represent the step, then we can define a change in an LCN quantity ( 2 . 71 ) By applying ( 2 . 71 ) to the quantities in ( 2 . 68 ) , a discrete step equation is expressed as ( 2 . 72 ) With the description of the stepped solver, it should be noted that the term is now defined in terms of differential susceptibility, no longer normal susceptibility. As such, the differential susceptibil ity can be found from a combination of the simulation history as well as from With this, a hysteresis model must be provided a priori that will appropriately characterize the differential susceptibility of the material. As can be deduced from ( 2 . 72 ) , the discrete change in are calculated rather than the total. As such, the total ma gnetization and field are calculated as a secondary value, and are tallied through the entire simulation where The generalized algorithm  is shown in Table 3 .
26 The stepped solver algorithm is provides a time dependent solution for M , H , and for non linear hysteretic materials . It also inherently computes induced permanent magnetization. Secondly, the primary system matrix need only be computed once and stored, then the block diagonal portion of the system matrix that arises due to the differential susceptibility is computed for each iteration k . Also by using the LCN discretization of the quasi magnetostatic VIE, the choi ce of the test vectors Table 3 Algorithm for the Stepped Solver 1. 2. Compute the system matrix without material term 3. a. Fill material term usi ng b. Solve ( 2 . 72 ) for with applied excitation c. Compute and field where i. ii. d. as the reciprocal unitary vectors leads to the diagonalization of the material term. Therefore, prior to solving ( 2 . 72 ) , only the diagonal block must be updated at each time step. An iterative solution using the GMRES algorithm is used to solve ( 2 . 72 ) . Lastly, in order for the
27 calculation of in step 3c from Table 3 to remain valid, the step size in k must be chosen such that the loc ally linear approximation of the differential susceptibility will remain valid, as is the case in all locally linear approximations of nonlinear processes. Young  states that the optimum step size is dependent on the model a nd the excitation waveform, and it is not required that the step size be dependent on a time step. It should also be noted that in order to check for solution convergence, a smaller step can be used in a subsequent solve to ensure that the solution has con verged as expected.
28 Chapter III MATERIAL PARAMETERS AND DIFFER E NTIAL SU SCEPTIBILITY Material Parameters Demagnetization When measuring magnetization of ferromagnetic materials, it is important to start measurements from a known state of magnetization, whether that be at saturation or from a demagnetized state. There are multiple methods that can be used to demagnetize a ferromag netic sample. The most effective way to demagnetize a sample is to raise the temperature of the sample above the Curie temperature. This allows all of the magnetic moments in the sample to be overcome by the thermal oscillations at the atomic level, and as the sample returns to room temperature, the bulk magnetization orientation of the domains in the sample are no longer present. The Curie temperatures of some common magnetic elements can be found in Table 4 . When a macroscopic s ample of material is above the Curie temperature, the crystalline structure oscillates wildly enough that the magnetic moment of the individual atoms are no longer oriented above the microscopic scale, and thus the sample has a net zero magnetization on a macroscopic scale. This is not to say that the magnetic moments no longer exist, but rather that the net effect is that . When the sample begins to cool below the Curie temperature, the grains in the ferromagnetic sample begin to align and show magnetism, but the net magnetization again remains 0. This is known as thermal demagnetization, and is the best method for demagnetization.
29 Table 4 Curie Temperature for commonly studied magnetic materials ELEMENT NAME SYMBOL CURIE TEMPERATURE [K] IRON Fe 1043 K COBALT Co 1394 K NICKEL Ni 631 K The main problem with thermal demagnetization for large samples, or for steels that have strength characteristics based on their annealing process, is that applying a l arge amount heat uniformly and cooling it slowly and fairly uniformly is difficult to do properly . If the sample is not cooled uniformly, then a portion of the sample that cool s and has magnetic grains can create an internal magnetic field that has a net o rientation at the domain wall between the cooled and non cooled domains of the sample. This will , in effect , apply a field to the domain above the Curie temperature and cause further magnetization as the hot domain cools. This process can continue on and o n, thus nullifying the thermal demagnetization process and yielding magnetization in the sample. Another demagnetization process is known as the alternating field, or AC, demagnetization process. This process is one in which a decaying, alternating magnet ic field is applied to a sample until the sample has reached a demagnetized state. While this method is not theoretically as effective as the thermal demagnetization method, it is more than sufficient in almost all cases. Most VSM machines will give the de magnetization attenuation
30 as a percent of maximum magnetization, based on the sample volume and highest level of magnetization from the applied field. The AC demagnetization process is set up by applying a sinusoidal magnetic field with frequency and amplitude that exponentially decays as where is the decay rate. The AC demagnetization process works by aligning all of the domains inside the sample w ith the applied field to the point of saturation. While it is not required that the sample be driven into saturation, the process does work effectively when starting from a saturated state . Note , if the maximum applied field for the demagnetization process were to be doubled, the exponential decay model would still rapidly bring the starting value down to the original proposed value, without adding any significant operation time. Thus, any applied field that is large enough to drive the signal into a semi s aturated state would be adequate for the demag netization process. Assume the sample is illuminated to an external field starting at and This forces the orientation of the do mains in the sample to almost entirely flip orientation, except for a small portion that will not flip due to internal denucleation. This is done recursively until the final minimum applied field determined a priori is reached. After the process, the domai ns and grains are again oriented in such a way that the net magnetization of the sample is approximately zero, while on a local scale there will be some level of magnetization. For example, a thin mild steel alloy disk has a saturation magnetization of 1.6 7e6 [A/m]. If an applied field [A/m]. If an applied field is chosen so that the signal can be driven far
31 driven far into saturation, and a demagnetization threshold is chosen based on a minimum applied H field. If we let that minimum be , th en Table 5 shows the approximate number of field flips required to demagnetize the sample. Table 5 VSM demagnetization table showing estimated time for VSM sample demagnetization H init applied Decay rat e Time per flip Number of flips Demagnetization time (mm:ss) 10000[Oe] 98% 5 sec 456 38:00 10000[Oe] 95% 5 sec 180 15:00 10000[Oe] 90% 5 sec 88 7:20 In our lab, the only demagnetization process used to date is the AC demag netization method. Figure 5 shows the remnant magnetization of a thin steel disk made from MIL S 16216K , or mild steel. The disk itself is an 8mm diameter, 0.41mm thick sample that is placed transverse to the applied magnetic fiel d. The experiment was to examine the remnant magnetization of the disk after using two different decay rate values to see if there was any significant difference between the two. Each sample was demagnetized, the blue curve is the plot of the virgin curve with a decay constant and the red curve is a constant of The terminating criteria for the demagnetization process using the VSM is when the applied field is to be less than Th is is enforced by the software, and to date , has not been able to be modified.
32 Examining the magnetization of both samples after the 90% and 95% decay rate constants, it is shown that either process will yield a high quality demagnetization of the sample. Both samples are demagnetized to less than 0.05% saturation magnetization, which is deemed satisfactory for sample demagnetization. Figure 5 Plot of magnetization remaining after AC demag process using the VSM at 95% and 90% d ecay rates. The red curve is the 9 5 % decay rate, and the blue curve is the 9 0 % decay. A high fidelity demagnetization is extremely helpful when attempting to measure material parameters, particularly the initial susceptibility of the material. Initial S usceptibilit y An important physical parameter of a ferromagnetic material is the initial susceptibility This parameter can be described as the magnetic susceptibility of a material exposed to a small magnetic field when the mate rial is in a demagnetized, or de permed, state. For elemental magnetic materials such as nickel, cobalt, and iron, the initial
33 susceptibility can be calculated analytically, but the creation of alloys with real world application , physical measurement of th e initial susceptibility is the simplest method . In order to measure the initial susceptibility, the most efficient way is to apply a very small magnetic field to a ferromagnetic material, and measure the subsequent magnetization. The initial susceptibili ty can then be approximated as ( 3 . 1 ) where i s the resulting change in magnetization, and is the change in applied magnetic field. Note, should be as small as possible in order to remain in the reversible region of magnetization for the sample. Fo r example, using the data obtained in Figure 5 , the initial susceptibility can be calculated from the demagnetized state for both the 90% and the 95% AC demagnetization process. The results are shown in Table 6 . Table 6 Calculation of Initial Susceptibility from demagnetization in Figure 5 . 0.9 33.24 182.17 1125.50 11523.59 84.93 0.95 30.81 181.92 753.16 11496.76 81.07 By applying the approximation to the data in Figure 5 , the initial susceptibilit ies are larger than predicted and previously measured by Schneider  . To note, t here are currently software limitations in the VSM that have not been resolved yet that is preventing highly
34 accurate d emagnetization processes and measurements as of the time of this writing, so the Faraday method was used to measure the initial susceptibility of the MIL S 16216K steel. The design was to demagnetize a thin wall hollow pipe, the n apply a magnetic field starting from and stepping up to with The applied field was then decreased to using the same ste p size. This process was repeated 10 times, yielding the MH curve in Figure 6 . The shape of the curve shows that the relationship between MH is approximately linear, and by using the relation and averagi ng susceptibility, the initial susceptibility which is very close to that predicted and used by Schneider  . Figure 6 MH curve for applied fie ld on thin wall hollow pipe. Experiment was run 10 times to eliminate noise and increase the validity of the average .
35 Magnetic Saturation Magnetic saturation in ferromagnetic materials occurs when an increase in an applied external H field is not capabl e of increasing the magnetization of the material any more. This is seen on the MH curve in Figure 7 where the curve flattens out on the top and bottom of the hysteresis curve. Physically speaking, magnetization saturation is when all of the domains inside the sample are oriented in the same direction as that of the applied external field. Figure 7 shows a mild steel disk driven into saturation. At the top and bottom of the curve, around the magnetization measurement flattens out, which is the point which magnetic saturation occurs. Magnetic saturation can also be said to have been reached when the susceptibility of the material is 0, or This means tha t the material no longer has any crystalline elements in the structure capable of aligning with the applied magnetic field. Similar to initial susceptibility, saturation magnetization is dependent on multiple material properties  . Annealing, impurity location, impurity compounds, internal stresses, and temperature will all have an effect on the magnetic saturation value that is measured. As such, the simplest way to find this value is via physical measurement for specific m aterials. Coercive Field The coercive field is the field required to flip the magnetization inside of a hysteretic material. This occurs on an curve at the point in which magnetization flips from posit ive to negative, indicating that the magnetization has now switched direction, i.e. the poles have flipped positions.
36 Figure 7 Hysteresis (MH) plot of MIL S 16216K mild steel disk driven into magnetic saturation. Figure 8 shows the MH curve of the hollow pipe starting from a demagnetized state, and driven into positive saturation, then negative saturation. The coercive field can be seen where the blue curve crosses the curve. The field was applied axially to the hollow pipe, and was demagnetized prior to the test. Cooperative Differential Susceptibility Model The cooperative hysteresis model posed by Schneider  states that as the scale of observable effects increase from that at the atomic level to the macroscopic level, complexity decreases significantly and ferromagnetic hysteresis can be model using an empirical form. The model is based on material parameters and , which are defined as the initial magnetic susceptibility, saturation magnetization, and coercive field, respectively, and are described more in depth earlier in the chapter.
37 Figure 8 MH curve of MIL S 16216K st eel pipe driven into saturation . The coercive field remains the same as long as the reversal occurs in the saturated range. The normal susceptibility is defined as ( 3 . 2 ) where is the change in magnetization since the last reversal, a nd is the change in total field, and is the total internal field. The cooperative model expresses as the superposition of the reversible susceptibility, , and the irreversible susceptibility: ( 3 . 3 ) where or the change in magnetization since the last reversal. If beginning from a virgin state, then and . Otherwise, is the magnetization at the last
38 reversal, and H curve will form a hysteresis loop when H is cyclical. The slope of this loop is defined as the differential susceptibility  ( 3 . 4 ) The texture of the material, in the case of MIL S 16216K steel or mild steel, is a cubic p olycrystal . Bozorth  derived a function for in BCC polycrystals of iron where ( 3 . 5 ) and ( 3 . 6 ) Schneider states that is also found to be a function of magnetization normalized by the saturation magnetization as well  . The formulations are shown in ( 3 . 7 ) . ( 3 . 7 ) The strength of the cooperative model is for small to moderate applied magnetic fields, ferromagnetic hysteresis can be modeled accurately. The model is based on physical parameters of the material, and the model can also solve for the magnetiza tion throughout the discretization of the geometry . This allows for the researcher to examine the approximate occurs. Other characteristics can also be solved for, such as the permanent magnetization, magnetic signature of the sample, as well as the motion of the poles as applied field and magnetization changes. The downfall of the cooperative model is that when a sample is driven to fields where the internal field the model breaks down and begins to
39 calculate nonphysical values of magnetization for cells throughout the model. This will be discussed further in Chapter V where an example of a thin wall pipe will be driven to saturation and the curve will show the result of the failure . Exponential Cooperative Differential Susceptibility Model One of the limitations of the cooperative differential susceptibility model is that when the internal field is the cooperative model breaks down and becomes unstable. Dr. Carl Schneider developed a new form, known as the Exponential Cooperative Model, or more succinctly, the exponential model, to more accurately describe ferromagnetic hysteresis over a bro ad range of applied fields, including when the internal field is greater than the coercive field of the material  . T he exponential differential susceptibility model is defined as ( 3 . 8 ) where is the change in magnetization since the last reversal, is the sign of the difference in magnetization, on the virgin curve, otherwise, and is the coercive field. is the denucleation field, is the initial magnetic susceptibility, is the saturation magnetization, and is the polycrystalline texture, defined in ( 3 . 6 ) . As with the cooperative model, on the virgin curve T he exponential model was developed in order to be a physics based model capable of modeling the differential magnetic susceptibility s olely as a function of physical material parameters and current state of magnetization.
40 Also attempted in the exponential model was to include a viscous susceptibility term This term was superimposed with the exponential suscep tibility model in ( 3 . 8 ) such that ( 3 . 9 ) in an attempt to account for rate dependent effects of magnetization. The viscous term is approximated as ( 3 . 10 ) where is the thermal field which is linear in and is a time cons tant of the material, which is on the order of While the model was promising, the preliminary simulations the inverse field rate of change term when the applied field and coercive field inside the material were very close in magnitude. During the simulation, t he viscous term would dominate the exponential term because as Also, was held constant during the simulations in order to more accurately examine the impact of the term. Upon further discussion with the entire research team, the decision was made to remove the viscous susceptibility term and continue to study the other models more in depth. Parallel Differential Susceptibility Model Another model dev eloped by Dr. Schneider and Dr. Gedney is known as the Parallel model  . In the parallel model, the reluctance of a cooperative susceptibility and saturate susceptibility are added in parallel, along with a maximum susceptib ility term that prevents the model from creating nonphysical susceptibilities within the model. is the cooperative susceptibility, defined as
41 ( 3 . 11 ) the saturate susceptibility term , is defined as ( 3 . 12 ) L astly, is the maximum susceptibility for the specific material. As a note, the maximum susceptibility is related to D via the relation where is the demag factor for the sample. The demag factor is dependent on both the material, as well as the geometry  . is the saturate magnetization, is the denucleation magnetic field , is the initial magnetic susceptibility, is the magnetization at the last reversal, and The parallel differential susceptibility is then defined as a sum of parallel reluctances  ( 3 . 13 )
42 CHAPTER IV MEASUREMENT METHODS Vibrating Sample Magnetometer Figure 9 General VSM configuration showing sample vibration orientation inside the applied field and pick up coil locations. A vibrating sample magnetometer is a sophisticated apparatus that , in general , is used to measure the m agnetic properties of ferromagnetic materials. The general configuration can be seen in Figure 9 . A uniform magnetic field H is applied between two electromagnetic coils, and focused by poles, creating an approximately uniform mag netic field that the sample is placed inside. With the sample placed inside the magnetic field, the sample is then vibrated vertically using a mechanical oscillator at a desired frequency. This vibration causes
43 a temporal change in the magnetic flux B , whi ch in turn induces a voltage in the sensing coils proportional to the magnetic moment of the sample. In the case of ferromagnetic steel, it is known that the MH curve will show hysteresis , vious magnetization. This is one of the main areas of research in the University of Colorado Denver Magnetics Research Lab. With the VSM, we are capable of driving the applied magnetic field to very high values (approximately 15 kOe), which is well above t he required field to drive our thin disk samples into magnetic saturation. Figure 10 shows a hypothetical ferromagnetic hysteresis curve. is the maximum amount of magnetization that a given material can hold , or the saturation magnetization . It occurs when all domains inside the sample are aligned with an applied magnetic field. This value changes for all materials, and is based on the chemical composition of the material under test among other things . is the remnant magnetization of the sample when the applied magnetic field H is brought back to zero, and the coercive field is the field at which flips sign . Lastly, the virgin curve starts at a virgin state, where The outer curve in Figure 10 is referred to as a saturate loop. changing magnetic fi eld will induce a voltage in a closed circuit loop. Quantitatively, we can describe the effect in equation ( 4 . 1 ) . ( 4 . 1 ) Thus, the magnetic moment of the sample will induce a voltage in the pickup coils (both x and y oriented), and can be converted into a directionally dependent hysteresis curve. This
44 Figure 10 Hysteresis cur ve for ferromagnetic materials will allow for future study into geometrical dependencies of samples, but as of the time of this writing has not b een explored. One of the benefits of the VSM is that very accurate measurements of the magnetic moment can be taken and studied. This is very beneficial when measuring a sample to find the coer cive field, magnetic remanence , saturation m agnetization , and initial susceptibility . A disadvantage is that extracting local magnetization from the moment must make the assumption of uniform distribution of magnetization, which unfortunately, is generally not the case. Faraday Method A second wa y to measure magnetism is referred to here as the Faraday method, in which the flux in the steel of a hollow pipe is back calculated via a Faraday coil wrapped around a pipe
45 with turns. Figure 11 shows t he physical setup of CU Denver Magnetics Lab Faraday method measurement. The smaller coil located in the center of the apparatus is the Faraday coil in which the voltage is induced . T he larger, red solenoid on the outside is the source of the uniform magne tic field which illuminates the thin wall mild steel pipe. When in test, this will be centered over the sample. Figure 11 Solenoid and Faraday coil set up to measure magnetization of hollow steel pipe. The Faraday coil is wra pped tightly around the thin wall hollow steel pipe. From across a winding that is equal to the negative time varying magnetic flux cutting through the surface of the coil , ( 4 . 2 ) The Faraday coil contains N turns . T he flux cutting through each turn will induce the same voltage, then t he net voltage induced across the N turn winding is: ( 4 . 3 )
46 From equati on ( 2 . 2 ) , is decomposed into the and components . Then the induced voltage can b e expressed as ( 4 . 4 ) Assuming that the illuminating magnetic field in t he solenoid is constant through the cross section of the hollow pipe, we can again make an approximation of and as ( 4 . 5 ) where a is the inner radius of the pipe, and b is the outer radius. Now, the x directed magnetization can be solved for as ( 4 . 6 ) Thus, by taking the discrete integral of the voltage, we can approximate the magnetization through the cross section of the thin pipe. In order to measure the field internal to the pipe, a Hall probe is placed in the center of the pipe which will measure the intensity of the axial field, which is a superposition of the external ly applied field as well as the demagnetizing field created by the magnetization of the pipe . The axial magnetic field can be assumed to be constant through the pipe wall. Thus, the Hall probe measurement provides a good e stimate for and consequently ( 4 . 6 ) can be used to measure the axial magnetization induced in the cross section of the steel pipe.
47 Chapter V VALIDATION Cooperative Model One objective of this thesis, as well as the CU Denver Magnetics Research Lab is to implem ent non linear hysteresis models with in comparing measured result with those solved using the LCN method. The measured data was obtained by Sean Joyce, through collaboration with Dr. Carl Schneider and Dr. Stephen Gedney. With their assistance, comparison of measured and simulated data in house is now possible. To begin, consider a hollow pipe manufactured using MIL S 16216K steel, or mild steel. The pipe has an overall length of 203.2 mm, with an inner radius of 9.0 mm and an outer radius of 10.5 mm, yieldi ng a thickness of 1. 5 mm. The discretization of the pipe is displayed in Figure 12 . The mesh had 10 cells around the circumference of the pipe and only required one cell through the cross section. susceptibility saturation magnetization and a coercive field were used. The applied field is along the axial direction of the pipe, or direction. at the pipe center. To compute against measured data, the net magnetization cutting the plane divided by the cross se ctional area was computed . The total magnetic field, was also computed at the center of the pipe via an H to compare against measured data .
48 Figure 12 Quadratic, quad rilateral m esh discretization of thin hollow pipe to be used in Magstr m . The cylinder was assumed to be in a virgin state at time and was subjected to an applied field that increased linearly from to 1000 [A/m], then decreased linearly to 1000 [A/m], then returned linearly to the initial applied value of 1000 [A/m]. The discrete applied field step, was chosen to match the rate of the measured data. The measured data was ob tained using a rate slow enough to minimize the impact of viscous effect s and eddy currents inside the steel pipe , but fast enough to minimize thermal heating, and subsequent modifications of the physical parameters . The same field step size as the simulat ion were applied to the pipe, but the measured data contained a repeat of the positive and negative reversals for verification of magnetic moment measurement of the pipe. Figure 13 shows the overlay of measured magnetization data compared to that found perative model. As noted, when the field inside of the pipe
49 Figure 13 Comparison of Cooperative model to measured data using the Faraday method. Measured data is the red line, while simulated data is the blue line. remains less than the cooperative model can very accurately model hysteresis of ferromagnetic hysteresis when the material para meters are known, or can be approximated. Figure 14 shows the results of the same simulations results as that run for the comparison to measured data, except that the sample is driven well into saturation using the cooperative model. The internal field exceeds the experimentally determined instabilit y limit and the hysteresis curve shows this by calculating unrealistic susceptibility values, and the simulation goes unstable . Thus, for the measurement of initial susceptibility and low magnetic fields, the cooperative model pr ovides a high level of confide nce for simulation, but can go unstable for large magnetic fields. A new measurement setup is under development for the VSM, in which two thin steel disks , 26mm in diameter, are placed in parallel planes to each other, separat ed by a distance
50 Figure 14 Simulation results of MH curve of hollow pipe when limit of the cooperative model is excee ded. The model failure occurs at the blue line, where the susceptibility of the material becomes negative, and the model rapidly becomes unstable. of 4 .4 mm face to face in order to measure the total H field in between the two disks. The sample holders are printed using PLA on a 3D printer, which allows the two samples to be located in the center of the VSM, and a Hall probe is located between the two samples. This will allow the experimenter to more easily measure internal H fields, and calculate demagnetization factor s at a distance from the sample without having a physical setup . Figure 16 shows the magnetization of each cell in the discretized volume mesh after a total applied field of 2500 A/m with steps of 50 A/m using the cooperative model. The material parameters are , , and . It can be seen that at the center of the disk is where the largest magnitude magnetization resides . Also, as expected, the lower disk that is semi hidden from view in Figure 16 has symmetric magnetization to that of the top disk. To find the location of the effective magnetic poles, 101 H field probes were placed along the axis.
51 Figure 15 Magnetization of parallel double disk s ample calculated using cooperative model in Magstr m. Applied field of 2500 [A/m]. During the development of the model, it was proposed by Dr. Schneider that as the magnetization of a steel disk increases, the effective location of the magnetic poles shou ld be pushed outward slightly. Using the cooperative model for differential susceptibility, the same simulation done for Figure 15 was run with 101 H field probes located along the axis in order to find the sign flip of the demagnetizing field The point where sign flips indicates the effective pole between the two disks. Figure 16 shows the field through the center of the two disks. As predicted by Dr. Schneid er, the poles do move by . The amount of pole motion is not yet fully understood, but evidence of pole motion from the cooperative model gives further confidence in the validity of the model. Additional v erification of the cooperative model on oth er canonical models, including a 1m magnetic spherical shell, a hollow cylinder, and the TEAM workshop problem 13 was performed by Young et. al.  which also show v ery h i g h
52 Figure 16 H field located along x axis between two parallel disks . H field sign change indicates location of effective magnetic pole. The effective pole location moves outward from the remnant position while the disk is illuminat ed by an applied H field. Figure 17 Zoomed in zero crossing of H d crossing to show pole location movement predicted by Schneider and verified with the cooperative model. agreement with physical measurements using the coope rative model .
53 Exponential Model Verification and understanding of the exponential model is still under development as of the time of this composition. In order to verify the exponential model, a test was run on the same mild steel hollow cylinder that has been used throughout this thesis. An applied field of was applied axially to the cylinder from a demagnetized state, near saturation , then returned to This process was repeated for 3 full hysteresis lo ops. The red curve in Figure 18 shows the magnetization versus total field at the center of the pipe. Figure 18 Measured hysteresis curve for thin wall hollow steel pipe overlaid with Magstr m exponent ial model simulated data. acceptable, insofar as the model did not break down when the sample was driven into saturation and followed a hysteretic curve indicating that the stepp ed method and exponential model are accurately accounting for hysteresis. Visible tilting of the simulated hysteresis curve can still be seen, indicating that the differential susceptibility being calculated is lower than the measured data.
54 The exponentia l model is still in the developmental phase, which caused the difference in the hysteresis curves between measured and simulated data in Figure 18 . While the simulated and measured data are slightly ajar, the simulation was able t o remain stable, showing that the model has a high possibility of accurate hysteresis modeling when the sample is magnetically saturated. The magnetization process of the steel pipe is a complex process which is not fully understood yet, which makes the m odeling of ferromagnetic hysteresis a very complex problem that has yet to be solved. Figure 19 shows what are known as First Order Reversal Curves (FORCs). The plot is logarithmic for susceptibility, and linear for magnetization. Figure 19 First Order Reversal Curves (FORCs) showing differential susceptibility versus magnetization. The process of obtaining FORC data is to drive a sample into saturation, and reverse the applied field to a number of dif ferent points to examine differential susceptibility when reversing field close to the point of maximum susceptibility. From the curves, Dr. Schneider
55 has posed that an exponential model of differential susceptibility should be able to accurately describe the hysteresis process. This is the driving force behind the development of the exponential model. Comparison of Cooperative and Exponential Model With verification of the cooperative model at small magnetization values, one should expect that the exponen tial model would also approximate small magnetization values accurately. In order to do this, we can examine equations ( 3 . 4 ) and ( 3 . 8 ) . To do this, a Taylor series expansion of both will be done, and resulting terms compared. To begin, the exponential form is expanded, and eliminating all but the linear term as the y are assumed to be much less than 1. ( 5 . 1 ) Next, we expand the 1/(1 x) term using Taylor series again and expand to arrive at the approximation in ( 5 . 2 ) ( 5 . 2 ) where is the coercive field in the exponential model. This is differentiated because it is thought that this is a different physical parameter than initially thought. The cooperative model from ( 3 . 4 ) is expanded similarly , again using only the linear term. The result is shown in equation ( 5 . 3 ) ( 5 . 3 )
56 Direct comparison between the expansions of the cooperative and exponential model show that the two leading terms are the same, which shows that there are similarities between the two at low magnetizations, and both should be valid for low magnetizations. It was shown that for higher internal fields, the cooperative model breaks down, thus an in depth analysis of the remaining terms of the expansions is not considered . ( 5 . 4 ) With this comparison of the expanded differential susceptibilitie s, a comparison between the two models at low applied field is now performed wall hollow pipe . The parameters for each simulation are listed in Table 7 . Table 7 Material parameters for cooperative model versus exponential model simulation of thin wall hollow pipe . X i M s H c H sat Cooperative 68 1.66 e6 A/m 800 A/m N/A Exponential 72 1.6 6 e6 A/m 400 A/m 2500 A/m The simulation started from a de magnetized state, and the applied field was increased to 1500 [A/m] in 50 [A/ m/sec ] increments to 1500 [A/m] , then reversed to 1500 [A/m] using the same 50 [A/m] steps, then again returned to positive magnetization at 1500 [A/m]. Figure 20 show s interesting characteristics from each of the respective models.
57 Figure 20 Plot of magnetic hysteresis of thin wall mild steel pipe using cooperative and exponential model. Applied field of 2500 [A/m] applied in the axial direction of the pipe. Examine of the hysteresis curves in Figure 20 sh ows good agreement between the exponential and cooperative model. A second verification of the cooperative model, exponential model, and measured hollow pipe made of MIL S 16216K steel was performed by Gedney et. al.  which also shows good agreement between the two models and measured data. The model parameters used in the simulation are listed in Table 7 . In comparison to the simulation done in Figure 20 , it can be seen th at the cooperative field used in the simulation is smaller than that used in the cooperative model. This is an area that will have further investigation as to the cause, and possibly lead to modification of either the exponential model, or the development of another model. Future Work The development of the aforementioned differential hysteresis models continues to be a work in progress. The highly complex physics that occur in magnetization are still not fully understood, and as such, multiple models must be thoroughly test ed until a
58 Figure 21 Comparison of hysteresis curves from exponential, cooperative models overlaid with measured data. Measured data taken by Dr. Schneider. . solution is found, o r an analytical solution is derived. The development of the CU Denver Magnetics Research Lab will prove crucial in the advancement of the application of the models and allowing immediate comparison to physically measured data. The exponential model has be en modified multiple times in an attempt to more accurately fit the hysteretic data of both the pipe and the thin disks of different materials. To date, collaboration with Dr. Schneider has allowed our lab to increase the reliability and validity of our me asured data. The data now being acquired is being studied in depth by Dr. Schneider in an effort to further understand the underlying physical pr ocesses. The in house lab also allows the lab team to apply the new models to a variety of samples including pu re nickel, cobalt, as well as other steel alloys and compare the differential susceptibility calculated by the model to that of the measured data.
59 The development of the susceptibility models will continue to be refined, and in the next 1 2 years, magneti c stress effects are scheduled to be studied. The Faraday method described in this thesis was designed to be inserted into an MTS load frame, which is an apparatus capable of applying a wide range of compression and tension forces to the hollow pipe. This will allow the measurement of magnetization changes due to these stresses. those that contain a permanent magnetization. The VSM contains an oven that will allow th in disks to be subjected to both very high and extremely low temperatures that will allow measurement of magnetization property changes due to thermal effects. Understanding the thermal effects on magnetization for ferromagnetic materials is important beca use of the variety of conditions that sea bearing vessels can encounter.
60 LIST OF REFFERENCES  C. S. Schneider, "Maximum Susceptibility of Ferromagnetic Hysteresis," Maximum Susceptibility of Ferromagnetic Hysteresis, vol. 48, no. 1 1, pp. 3371 3374, 2012.  J. C. Young and S. D. Gedney, "A Locally Corrected Nystrom Formulation for the Magnetostatic Volume Integral Equation," IEEE Transactions on Magnetics, pp. 2163 2170, 2011.  S. Gedney, Numerical Quadrature Lecture, Univ ersity of Colorado Denver, 2016.  J. A. Stratton, Electromagnetic Theory, New York: McGraw Hill, 1941.  S. D. Gedney and J. C. Young, "The Locally Corrected Nystrom Method for Electromagnetics," in Computational Electromagnetics Recent Advances and Engineering and Applications , New York, Springer, 2013, pp. 149 198.  J. C. Young, S. D. Gedney, R. Adams, C. Schneider and C. Burgy, "A Stepped Nonlinear Solver for Nonlinear Magnetic Materials With Hysteresis," IEEE Transactions on Magnetics, v ol. 51, no. 6, 2015.  W. Hafla, A. Buchau and W. M. Rucker, "Consideration of Scalar Magnetic Hysteresis with the Integral Equation Method," IEEE Transactions on Magnetics, vol. 44, no. 6, pp. 910 913, 2008.  C. T. Kelley, Iterative Methods for Linear and Nonlinear Equations, Philadelphia: Society for Industrial and Applied Mathematics, 1995.  R. C. O'Handley, Modern Magnetic Materials Principles and Applications, John Wiley & Sons, Inc., 2000.  R. M. Bozorth, Ferromagnetism, Prince ton, N.J.: D. Van Nostrand, 1951.  J. C. Young, "Magneto Stress Analysis of Complex Structures," Technical Report EM 1, University of Kentucky, Department of Electrical and Computer Engineering, Lexington, 2016.  S. D. Gedney, C. Schneider and R. Geyer, "Measurement and Modeling Magnetic Material Properties of Steel," Technical Report CUDEN 1, University of colorado Denver, Dept. of Electrical Engineering, 2017.  D. X. Chen, J. A. Brug and R. B. Goldfarb, "Demagnetizing Factors for Cylind ers," IEEE Transactions on Magnetics, vol. 27, no. 4, pp. 3601 3619, 1991.
61  G. D. Mahan, "Crystal Magnetism | physics | britannica.com," Encyclopedia Brittanica, 2015. [Online]. Available: https://www.britannica.com/science/crystal/Magnetism#ref506523. [Accessed 15 March 2018].  R. D. Harrington and W. E. Case, "Calibration of Vibrating Sample Magnetometers," Journal of Research of the National Bureau of Standards, vol. 70c, no. 4, pp. 255 262, 1966.  R. M. Bozorth, Ferromagnetism, Piscataway: IEEE Press, 1993.  L. F. Canino, J. J. Ottusch, M. A. Stalzer, J. L. Visher and S. M. Wandzura, "Numerical Solution of the Helmholtz Equation in 2D and 3D Using a High Order Nystrom Discretization," Journal of Computational Physics, vol. 146, pp. 627 663, 1998.