New Constitutive Matrix in the 3D Cell Method to Obtain a Lorentz Electric Field in a Magnetic Brake

In this work, we have obtained a new constitutive matrix to calculate the induced Lorentz electric current of in a conductive disk in movement within a magnetic field using the cell method in 3D. This disk and a permanent magnet act as a magnetic brake. The results obtained are compared with those obtained with the finite element method (FEM) using the computer applications Getdp and femm. The error observed is less than 0.1173%. Likewise, a second verification has been made in the laboratory using Hall sensors to measure the magnetic field in the proximity of the magnetic brake.


Introduction
The use of magnetic brakes has obvious advantages compared to brakes based on mechanical friction. The mechanical friction brakes have the risks of hydraulic fluid loss or the contamination of the fluid by cooling water, with a consequent loss of braking power, among others [1]. The use of permanent magnets in linear magnetic brakes is explained in detail in [2]. In [1], a mixed system with mechanical friction and a magnetic brake is analyzed, comparing the results obtained with an analytical equation and the finite element method (FEM) in 2D and 3D.
In most of the works consulted in the bibliography, approximate analytical equations or numerical methods such as the FEM are used. In the present work, we propose the finite formulation (FF) [3], and the cell method (CM) as its associated numerical method [4,5], to analyze this type of device.
FF works with global magnitudes associated with spatially oriented elements such as volumes, surfaces, lines and points of the discretized space, as well as temporal elements, instead of field magnitudes associated with the independent variables of spatial and temporal coordinates [6]. In FF, the equations of constitutive type-equations of the medium-are clearly differentiated from the topological type-equations of balance. The analysis of the magnetic brake with this methodology facilitates, to a great extent, the boundary conditions and continuity, when working with global magnitudes.
The contribution of this work consists on obtaining a new constitutive equation that relates the Lorentz current with the magnetic potential using CM. This result applies to the magnetic brake.
This consists of a copper disk with an angular velocity → W r , and a permanent magnet is placed in front it.

The Constitutive Matrix in the Cell Method
The electric field strength, → E, measured in a fixed coordinate system with respect to the laboratory, is related to the electric field strength → E v , measured in a coordinate system referring to each point of a moving conductor, with velocity → v , with respect to the laboratory (1) [7], where v 2 c 2 , and c is the speed of light, → B is the magnetic induction, → a is the magnetic vector-potential and ϕ is the electric scalar-potential, where If we work with a permanent magnet, then ∂ → a ∂t = 0, and (1) is simplified. The current density is shown in (2), where σ is the volumetric electric conductivity. The Equation (2) contains σ → grad ϕ .
These are the electric currents produced by the corresponding electric potentials. This term of the equation is developed with CM in [8]. The first term of Equation (2) in CM is Equation (3). The electric currents refer to the electric conductivity constitutive matrix, M σ . The differences of electric potentials are associated to the edges e i , i = 1:6, of the reference tetrahedron (see Figure 1a).
where the matrix M σ = σ v ε S i S j , i, j = 1 : 6, is function of the electric conductivity of each tetrahedron, its volume and the dot product of the surface vectors that correspond to the dual planes of primal edges e i and e j (see Figure 1a).
We start from the expression of the magnetic induction, (4), within each tetrahedron as a function of the magnetic fluxes associated with the faces of the reference tetrahedron [9].

The Constitutive Matrix in the Cell Method
The electric field strength, E ⃗ , measured in a fixed coordinate system with respect to the laboratory, is related to the electric field strength E ⃗ v , measured in a coordinate system referring to each point of a moving conductor, with velocity v ⃗, with respect to the laboratory (1) [7], where v 2 ≪c 2 , and is the speed of light, B ⃗ is the magnetic induction, ⃗ is the magnetic vector-potential and φ is the electric scalar-potential, where If we work with a permanent magnet, then ∂a ⃗ ∂t = 0, and (1) is simplified. The current density is shown in (2), where σ is the volumetric electric conductivity. The Equation (2) contains σ grad ⃗ φ .
These are the electric currents produced by the corresponding electric potentials. This term of the equation is developed with CM in [8]. The first term of Equation (2) in CM is Equation (3). The electric currents refer to the electric conductivity constitutive matrix, M σ . The differences of electric potentials are associated to the edges ei, i = 1:6, of the reference tetrahedron (see Figure 1a).
where the matrix M σ = σ ϵ v ε S i S j , i, j = 1:6, is function of the electric conductivity of each tetrahedron, its volume and the dot product of the surface vectors that correspond to the dual planes of primal edges ei and ej (see Figure 1a). We start from the expression of the magnetic induction, (4), within each tetrahedron as a function of the magnetic fluxes associated with the faces of the reference tetrahedron [9]. The second term of (2), which corresponds to the Lorentz term, is developed in this article with CM, where The integration of the differences of the electrical potential associated to both primal edges ei, i = 1:6, with the Lorentz term E ⃗ Lo = V ⃗ e × B ⃗ e is expressed in Equation (5): The second term of (2), which corresponds to the Lorentz term, is developed in this article with CM, where The integration of the differences of the electrical potential associated to both primal edges e i , i = 1:6, with the Lorentz term B e is expressed in Equation (5): Solving (5), the U i in the primal edges e i is expressed in (6), where With the six edges of the reference tetrahedron, the Lorentz voltages-vector is shown in (9) and V Lo is calculated using (10).
In the reference tetrahedron, the magnetic fluxes can be expressed as magnetic potentials associated to the edges of the tetrahedron [6]. This is developed as (11), Substituting (11) in (9), we get (12), which is the proposed equation in this work. Matrix C is the incidence matrix between the faces and edges of the reference tetrahedron (see Figure 1a), The current produced by the Lorentz voltage is shown in (13). The total current will be I = I φ + I Lo : where MV L is the new constitutive matrix proposed for the calculation of the Lorentz current, I Lo , depending on magnetic potentials. To validate the proposed equation, the problem of the magnetic brake (MB) is posed using FF. The MB is a disk, made of copper, which rotates around an orthogonal axis, y. On the upper side of the disk there is a neodymium permanent magnet (see Figure 1b). The velocity for any point of the disk is To solve the problem, the electric-potential ϕ and a magnetic-potential formulation are used. In this formulation, the continuity equation for electric current, (15), is used, where D is the face-volume incidence matrix in the dual mesh, Taking into account the fact that D = −G t , where G is the edge-node incidence matrix in the primal mesh (see Figure 1a), then Equations (3), (12) and (13), with the continuity equation for electric current, and all terms-including the term ∂ → a ∂t -remain as indicated in (16): The definition equation for the permanent magnet is Equation (17). M ν is the magnetic reluctivity matrix [4], F c is the vector of coercive magnetomotive forces supplied by the manufacturer of the magnet, and F is the vector of magnetomotive forces associated to the edges of the dual mesh: The constitutive equation in the air and the conductive material is (18). If Ampere's law (19) is applied to the permanent magnet, then it reads as (20). If the same law is applied to the air, it is shown as (21). Finally, applying Ampere's law to the conductive material, then we obtain (22), where C is the discrete curl operator in the dual mesh. Taking into account the fact that C = C t and φ = Ca, then Ampere's law and continuity equation for electric current in all domains both form a system of equations as in (23). This system of equations, in a sinusoidal steady state, has been programmed in C++, using the numerical software package PETSc [10].

Characterization of the Magnetic Brake
The magnetic brake consists of a disk of copper alloy and a permanent magnet located at a distance d (see Figure 2). The disk is coupled to a DC motor. The DC motor rotates at angular velocity → W r (see Figure 3). The permanent magnet, which is used in the numerical simulations and experiments at the laboratory, is modeled in the second quadrant, with H < 0 and B > 0. It is a linear model with two parameters, B r and µ. This model is adequate when rare earth permanents magnets are used. The characteristics of the permanent magnet are the following: the material us NdFeB, with a block form of 50.8 × 50.8 × 25.4 mm, a magnetization sense in axis y along the dimension 25.4 mm, the coating is nickel-plated (Ni-Cu-Ni), manufactured by sintering, the magnetization type is N40, the remanence B r is in the interval 1.26-1.29 T, the coercive H c is in the interval 860-955 kA/m, the intrinsic coercive H c ≥ 955 kA/m and the maximum energetic product is in the interval 303-318 kJ/m 3 .
The disk has a volumetric electric conductivity of σ = 4.1·10 7 S/m, its diameter is 315 mm and its thickness is 5 mm. When the equation system in (23) has been solved, the losses by Joule effect are obtained in the disk. These are calculated using a volume-integral as in (24), where → J is the volumetric density electric current indicated in (22). The braking torque, along the y axis, is calculated using a volume-integral as indicated in (25).

Numerical Simulations
Different numerical experiments have been developed using different angular velocities with an air-gap between the disk and the permanent magnet equivalent to 6 mm. The constitutive matrix of velocities proposed in (13), developed in CM and Equations (23), have been used in these simulations. The arising system, Equation (23), is not as symmetric as in FEM. We use numerical methods based in the subspace of Kryslov due to the fact that the matrices are sparse and of large dimensions. These algorithms are implemented in the software PETSc [10].
In particular, the linear solver employed is the generalized minimal residual algorithm (GMRES). Besides this, we precondition the matrix to improve the condition number of the matrix to reduce the number of iterations. This method is valid for non-symmetric systems. Relative and absolute tolerances with an order of magnitude of 10 −10 are enough to achieve convergence for the linear solver.

Numerical Simulations
Different numerical experiments have been developed using different angular velocities with an air-gap between the disk and the permanent magnet equivalent to 6 mm. The constitutive matrix of velocities proposed in (13), developed in CM and Equations (23), have been used in these simulations. The arising system, Equation (23), is not as symmetric as in FEM. We use numerical methods based in the subspace of Kryslov due to the fact that the matrices are sparse and of large dimensions. These algorithms are implemented in the software PETSc [10].
In particular, the linear solver employed is the generalized minimal residual algorithm (GMRES). Besides this, we precondition the matrix to improve the condition number of the matrix to reduce the number of iterations. This method is valid for non-symmetric systems. Relative and absolute tolerances with an order of magnitude of 10 −10 are enough to achieve convergence for the linear solver.

Numerical Simulations
Different numerical experiments have been developed using different angular velocities with an air-gap between the disk and the permanent magnet equivalent to 6 mm. The constitutive matrix of velocities proposed in (13), developed in CM and Equations (23), have been used in these simulations. The arising system, Equation (23), is not as symmetric as in FEM. We use numerical methods based in the subspace of Kryslov due to the fact that the matrices are sparse and of large dimensions. These algorithms are implemented in the software PETSc [10].
In particular, the linear solver employed is the generalized minimal residual algorithm (GMRES). Besides this, we precondition the matrix to improve the condition number of the matrix to reduce the number of iterations. This method is valid for non-symmetric systems. Relative and absolute tolerances with an order of magnitude of 10 −10 are enough to achieve convergence for the linear solver.
We have used tetrahedral elements because they are better for complex geometries. The number of elements depends on the model analyzed. For example, in a particular implementation, we used 261,868 elements, and the assembling and solution time was lower than 500 s, with an Intel core machine i7-3820, 3.6GHz and 32GB of RAM with four cores and eight threads of execution. The convergence has been studied, increasing the number of elements until a reasonable rate of convergence was obtained.
To avoid singularities in the problem solution, we control the dimensionless Péclet number to get a reasonable convergence, with manual mesh refinements, using a higher density of points where the linear velocity of the disk is higher.
The results obtained for the disk of copper alloy are shown in Figure 4a,b. In Figure 4a, the maximum current density in the disk is compared in CM and software package Getdp [11]. In Figure 4b, the losses by Joule's effect are compared. The proposed numerical simulations have been developed in CM. The reference, used to compare the results, is Getdp.
The software package Gmsh [12] is used in the mesh and the data visualization. In Figure 5, the normal force between the disk and permanent magnet is shown. The results obtained in CM are compared with the results obtained in FEM. The software package Gmsh [12] is used in the mesh and the data visualization. In Figure 5, the normal force between the disk and permanent magnet is shown. The results obtained in CM are compared with the results obtained in FEM.

Numerical Validation of the Simulations
Three experiments were developed and the results are compared with those obtained through CM and Getdp. These experiments are specified in Table 1.

C2:
Heat in the disk by eddy currents (see Figure 4b)

Numerical Validation of the Simulations
Three experiments were developed and the results are compared with those obtained through CM and Getdp. These experiments are specified in Table 1.
In Table 2, we can see that all the metrics are different from the optimum. MAE is given in true measurements with its units. This metrics show the validity of FF-CM for our analysis.
To understand the differences of magnitude in MSE, RMSE and MAE, we have to observe that C1 is a large quantity-with a current density on the order of magnitude 10 7 (see Figure 4a)-and both C2 and C3 are small-with a heat in the disk and normal force on the order of magnitude 10 2 (see Figures 4b and 5). Therefore, C1 is apparently a large quantity and both C2 and C3 are small quantities. However, as these quantities are absolute errors they seem large. The optimum is 0 because this is the limit to which these magnitudes tend as we increase the number of elements in the mesh. Table 3 presents the evolution of the error and the execution time for the assembling and solution for Joule heating as the number of degrees of freedom increases in the system of Equations (23) for an angular velocity of 15 rad/s. With 46,723 degrees of freedom, we obtain an error less than 3%, but we have implemented, in this particular case, 233,471 degrees of freedom obtaining a negligible error in 324.30 s.

Experimental Validation of the Simulation
In this section, we obtain experimental measures to validate the cell method used in the simulation of the magnetic field in the proximities of a magnetic brake. Besides this, we experimentally analyze the electromagnetic torque generated in the brake. In order to develop the experimental system, it is necessary to calibrate the Hall sensors that measure the magnetic field. These sensors will be previously calibrated through reference magnets.

Calibration of the Hall Sensors
First of all, we will characterize the reference magnets, and then we will calculate the gain of the Hall sensors. For the measurements, we have used precision weights that have been verified in a precision balance showing a negligible deviation. Distances have been measured with a precision of hundredth of a millimeter, which is a very good precision for our experiments.
Reference magnet characterization: The experimental measurement of the magnetic field on the proximities of the copper disk is performed with three Hall sensors that allow the measurement of the magnetic induction in the three directions of the Cartesian axes x, y, z. These sensors are of the 49E linear type in a measurement range between −90 mT and +90 mT [29].
In order to calculate the gain of each sensor, we use a consistent magnetic field source with two magnets of NdFeB of cylindrical shape and a radius of 2 mm and a height of 5 and 40 mm, respectively, as shown in Figure 6a.
The first set of measurements is performed in order to determine, in an accurate way, the intensity of the coercitivity field, H c , of the magnets, because it is unknown.
In this way, the two magnets are located as in Figure 6a, and the equilibrium force between both is measured for different distances between 2 and 14 mm. The mass used to obtain the equilibrium force consists of different precision weights from 2 to 100 g and are situated on the 40 mm magnet as can be observed in Figure 6a. We have taken two sets of experimental data to assure the repeatability of the measures and they are called Exp-0 and Exp-1, as shown in Figure 6b.
The intensity of the coercitivity field, H c , which better fits the data, is H c = 1.0902 MA/m. These data have been obtained by successive MEF simulations in 2D and 3D, finding the repulsion force between the magnets. Figure 6b shows the best fits obtained through MEF using Getdp 3D software and femm 2D software which takes advantage of the axial symmetries of this problem. Besides this, working in 2D, we can use denser meshes for the same computational cost, obtaining more accurate results.
These fits confirm that MEF is an adequate tool to calculate the magnetic field in the proximities of magnetic fields generated by permanent magnets.
Calculation of the gain of the Hall sensors: Although we use three Hall sensors of the same type and manufacturer [30], the gain of each one is slightly different. In this section, we experimentally characterize the gain of each one. These fits confirm that MEF is an adequate tool to calculate the magnetic field in the proximities of magnetic fields generated by permanent magnets.
Calculation of the gain of the Hall sensors: Although we use three Hall sensors of the same type and manufacturer [30], the gain of each one is slightly different. In this section, we experimentally characterize the gain of each one.
As Figure 7a shows, we locate the magnet using a high precision CNC-Computer Numerical Control-at a particular distance from the Hall sensor. Measuring the electric potential difference in the Hall sensor generated by the magnetic field of the magnet, we can determinate the gain of the sensor. The magnetic field was calculated using femm MEF software as explained in the previous section.
The gains obtained for Hall sensors x, y, z are 20.0, 17.5, and 21 mV/mT, respectively, as can be seen in Figure 7b.  As Figure 7a shows, we locate the magnet using a high precision CNC-Computer Numerical Control-at a particular distance from the Hall sensor. Measuring the electric potential difference in the Hall sensor generated by the magnetic field of the magnet, we can determinate the gain of the sensor. The magnetic field was calculated using femm MEF software as explained in the previous section.
The gains obtained for Hall sensors x, y, z are 20.0, 17.5, and 21 mV/mT, respectively, as can be seen in Figure 7b.

Measures of the Magnetic Field and Torque
As previously mentioned, we use a disk made of a copper alloy, of diameter 315 mm and thickness 5 mm that can rotate in both directions. This disk is moved by a DC motor regulated in speed by the electric potential of the armature.
The effective conductivity of the copper alloy has been obtained by minimizing the error between the measured and simulated magnetic field, varying the electric conductivity between 1 × 10 7 and 8 × 10 7 S/m, with an angular velocity of the disk of 20.49 rad/s, as Figure 8a shows. The obtained value for

Measures of the Magnetic Field and Torque
As previously mentioned, we use a disk made of a copper alloy, of diameter 315 mm and thickness 5 mm that can rotate in both directions. This disk is moved by a DC motor regulated in speed by the electric potential of the armature.
The effective conductivity of the copper alloy has been obtained by minimizing the error between the measured and simulated magnetic field, varying the electric conductivity between 1 × 10 7 and 8 × 10 7 S/m, with an angular velocity of the disk of 20.49 rad/s, as Figure 8a shows. The obtained value for a minimum error corresponds to 4.1 × 10 7 S/m. With this value, we have implemented all the numeric simulations from now on.
At a variable distance from the disk and with magnetization in the y direction, we have placed a permanent magnet. The magnet is a brick with the dimensions 50.8 × 50.8 × 25.4 mm. As shown in Figure 8b, we have also placed three Hall sensors in a trihedral trirectangle which determines a unique vertex. The trihedral can be displaced in the x, y, z directions.
Then, we analyze the magnetic field vector as a function of the angular velocity of the disk in the points 1, 2, 3 and 4, as can be seen in Figure 8b.
Experimental results and magnetic field discussion: We have taken eight measurements for different values of the angular velocity of the magnetic field in points 1, 2, 3 and 4 of Figure 8b. Four have been taken with the disk rotating clockwise and four anticlockwise. Figure 9 represents the module of magnetic induction measured by the Hall transducers and the numeric simulations of the angular anticlockwise velocity, Wr-. We observe that there is a great coincidence between experimental values and those obtained by numeric simulation because the discrete values are almost coincident with the continuous curves that represent the simulations. Besides this, the curves are increasing: that is to say, the magnetic induction increases as the angular velocity increases. Figure 8b, we have also placed three Hall sensors in a trihedral trirectangle which determines a unique vertex. The trihedral can be displaced in the x, y, z directions.
Then, we analyze the magnetic field vector as a function of the angular velocity of the disk in the points 1, 2, 3 and 4, as can be seen in Figure 8b.
Experimental results and magnetic field discussion: We have taken eight measurements for different values of the angular velocity of the magnetic field in points 1, 2, 3 and 4 of Figure 8b. Four have been taken with the disk rotating clockwise and four anticlockwise.  Figure 9 represents the module of magnetic induction measured by the Hall transducers and the numeric simulations of the angular anticlockwise velocity, Wr-. We observe that there is a great coincidence between experimental values and those obtained by numeric simulation because the discrete values are almost coincident with the continuous curves that represent the simulations. Besides this, the curves are increasing: that is to say, the magnetic induction increases as the angular velocity increases. On the other hand, Figure 10 represents the module of the magnetic induction, measured with the Hall sensors at the same points, and the numeric simulations, as functions of the angular   Figure 9 represents the module of magnetic induction measured by the Hall transducers and the numeric simulations of the angular anticlockwise velocity, Wr-. We observe that there is a great coincidence between experimental values and those obtained by numeric simulation because the discrete values are almost coincident with the continuous curves that represent the simulations. Besides this, the curves are increasing: that is to say, the magnetic induction increases as the angular velocity increases. On the other hand, Figure 10 represents the module of the magnetic induction, measured with the Hall sensors at the same points, and the numeric simulations, as functions of the angular On the other hand, Figure 10 represents the module of the magnetic induction, measured with the Hall sensors at the same points, and the numeric simulations, as functions of the angular clockwise velocity, Wr+. Likewise, we observe that there is a great coincidence between the experimental values and those obtained by numeric simulation. However, these curves are decreasing: that is to say, the magnetic induction decreases as the angular velocity increases. This behavior is due to the asymmetry of the currents induced in the disk when the rotating direction of the angular velocity changes (see Figure 11). clockwise velocity, Wr+. Likewise, we observe that there is a great coincidence between the experimental values and those obtained by numeric simulation. However, these curves are decreasing: that is to say, the magnetic induction decreases as the angular velocity increases. This behavior is due to the asymmetry of the currents induced in the disk when the rotating direction of the angular velocity changes (see Figure 11).    Experimental results and electromagnetic torque: This set of experiments consists of the measurements of the losses by the Joule effect made in the laboratory. These measurements have been obtained by subtracting the total electrical power at the input of the DC motor from the Joule losses in the armature windings and the mechanical losses in the non-load regimen, without the magnetic brake. With these calculations, the electromagnetic torque for braking is obtained for a determinate angular velocity (see Figure 12). The characteristics of the DC motor are the following: AEG, direct current, armature 220 V-2.2 A, field 220 V, 0.5 kW, 1400 rpm, IP E22. Figure 12 shows the electromagnetic torque obtained by experimentation and simulation in the copper disk for an angular velocity from 0 to 35 rad/s. The magnets are situated at distances of 9, 10 and 11 mm. Experimental results are almost coincident with numeric simulations. Besides this, as expected, we can observe that as the angular velocity increases, the torque increases. It can also be observed that as the proximity to the disk increases, the torque also increases.
As there is no induced current when the angular velocity is zero, the magnetic torque converges to zero in both experimental and simulation results. Figures 9, 10 and 12 show the experimental and simulation results. As can be seen, the differences are small. Figures 9 and 10 show the calculations and measurements of magnetic fields and Figure 12 shows the calculations and measurements of the electromagnetic torque. Experimental results and electromagnetic torque: This set of experiments consists of the measurements of the losses by the Joule effect made in the laboratory. These measurements have been obtained by subtracting the total electrical power at the input of the DC motor from the Joule losses in the armature windings and the mechanical losses in the non-load regimen, without the magnetic brake. With these calculations, the electromagnetic torque for braking is obtained for a determinate angular velocity (see Figure 12). The characteristics of the DC motor are the following: AEG, direct current, armature 220 V-2.2 A, field 220 V, 0.5 kW, 1400 rpm, IP E22. Figure 12 shows the electromagnetic torque obtained by experimentation and simulation in the copper disk for an angular velocity from 0 to 35 rad/s. The magnets are situated at distances of 9, 10 and 11 mm. Experimental results are almost coincident with numeric simulations. Besides this, as expected, we can observe that as the angular velocity increases, the torque increases. It can also be observed that as the proximity to the disk increases, the torque also increases.
As there is no induced current when the angular velocity is zero, the magnetic torque converges to zero in both experimental and simulation results. Figures 9, 10 and 12 show the experimental and simulation results. As can be seen, the differences are small. Figures 9 and 10 show the calculations and measurements of magnetic fields and Figure 12 shows the calculations and measurements of the electromagnetic torque.

Conclusions
A constitutive matrix has been developed into the cell method in 3D. The Lorenz electric current has been obtained using the magnetic potential.
This methodology has been applied to a disk that moves inside a magnetic field. The results obtained have been compared with the finite element method, using the free software packages Getdp and femm as a reference. The error is less than 0.1173%.
A second verification has been made in the laboratory using Hall sensors to measure the

Conclusions
A constitutive matrix has been developed into the cell method in 3D. The Lorenz electric current has been obtained using the magnetic potential.
This methodology has been applied to a disk that moves inside a magnetic field. The results obtained have been compared with the finite element method, using the free software packages Getdp and femm as a reference. The error is less than 0.1173%.
A second verification has been made in the laboratory using Hall sensors to measure the magnetic field in the proximity of the magnetic brake.