A FEM-Based Optimization Method for Driving Frequency of Contactless Magnetoelastic Torque Sensors in Steel Shafts

This paper presents a novel finite element method (FEM) of optimization for driving frequency in magneto-mechanical systems using contactless magnetoelastic torque sensors. The optimization technique is based on the generalization of the axial and shear stress dependence of the magnetic permeability tensor. This generalization creates a new possibility for the determination of the torque dependence of a permeability tensor based on measurements of the axial stress on the magnetization curve. Such a possibility of quantitative description of torque dependence of a magnetic permeability tensor has never before been presented. Results from the FEM-based modeling method were validated against a real magnetoelastic torque sensor. The sensitivity characteristics of the model and the real sensor show a maximum using a driving current of similar frequency. Consequently, the proposed method demonstrates the novel possibility of optimizing magnetoelastic sensors for automotive and industrial applications.

Soft magnetic materials the most commonly used as cores of contactless torque sensors utilizing magnetoelastic effects [26][27][28][29] are constructional steels [30][31][32]. However, research on the magnetoelastic effect caused by torque and possibilities of sensor optimization are significantly limited. Although the magnetoelastic effect caused by axial stresses is well described in the literature [33,34], information on the influence of shearing stresses generated by torque on magnetic characteristics of constructional steels is not readily available.
This paper addresses the gap in the state of the art by presenting the modeling results of the torque dependence of magnetic characteristics on drive shafts made of constructional steel. Following on from the previous work on the generalization of the axial and shear stress dependence of magnetic permeability tensors carried out in [35], this paper presents a novel approach to quantitative description of the torque dependences of magnetization process. Due to tensor description of magnetic permeability, the proposed approach enables considering the shear stress induced anisotropy of constructional steel magnetic permeability. Moreover, in the proposed model, the shear stress-induced anisotropy can be determined on the base of measurements of axial stress dependence of the magnetization curve [36]. This feature was not presented before in the literature. It should be highlighted that the possibility of determination of shear stress dependence of magnetic permeability tensor on the base of axial stress permeability dependence measurements is very useful from the practical point of view.
As a result, the proposed model may be efficiently applied for FEM-based modeling solutions. On the base of experimentally determined axial stress permeability dependences, the proposed model provides insight for understanding of magnetization processes under the influence of torque, as well as enables further contactless magnetoelastic torque sensors development. This development creates the possibility of optimizing magneto-mechanical systems using contactless magnetoelastic sensors in the area of modern industrial applications.

Finite Element Model for Development of Contactless Torque Sensors Utilizing Ferromagnetic Construction Steels
The description of the stress tensorσ dependence of relative magnetic permeabilityμ r tensor is the foundation of partial differential equation-based models of the magnetoelastic effect. Knowledge on the shear stress dependence of permeability is limited to specific systems and materials [37]. Moreover, axial and shear stress may occur in magnetomechanical systems simultaneously, which makes its description more complex.
It should be highlighted that, in the presence of stress induced anisotropy of magnetic permeability in materials, flux density B is not parallel to magnetizing field H. This phenomenon was neglected before [38]. In addition, due to sophisticated mechanical stress distribution as well as due to the simultaneous presence of axial and shear stress, only finite elements modelling (FEM) methods enable the comprehensive analysis of the stress dependence of the magnetization process in real systems.
Experimental results widely presented in the literature are mostly focused on the axial stress dependence of isotropic magnetic materials such as construction steels [38,39]. The axial stress dependence of relative permeability was previously estimated as a linear dependence given as follows [40]: where σ and µ r represent the axial stress (expressed in MPa) and relative permeability, respectively. However, it should be highlighted that a good agreement of Equation (1) with determination coefficient R 2 exceeding the value 0.94 was observed only for limited values of axial mechanical stresses σ. For the generalization of stress tensor dependence of relative magnetic permeability µ r tensor, the principal stress concept was applied [35]. Considering the Mohr's circle concept [41,42], each stress tensor can be reduced to the vector of axial stresses (principal stresses tensor σ P ) that is rotated to the new coordinate system by the rotation matrix R [36,43] Theoretical and experimental analyses presented in previous work indicate that effective axial stresses σ eff caused by the perpendicular stress σ ⊥ can be calculated according to the following dependence [42]: where ν is Poisson's ratio. The Poisson's ratio of constructional steels is typically equal to 0.3. Considering Equation (3), the efficient axial stresses influencing the relative magnetic permeabilityμ tensor can be estimated as: Consequently, the relative magnetic permeability µ r tensor can be calculated using the following equation: For the plain torque, only shear stress occurs in the stress tensor. In this case, Equation (2) can therefore be written as: Equation (7) is thus reduced to: Plain torque is rarely observed in real systems. For this reason, it is highly recommended to implement the Finite Element Modeling (FEM) systems where the stress tensor dependence of the relative magnetic permeability tensor is based on Equation (7).

Implementation of Proposed Model in Open-Source FEM-Oriented Software
The proposed model is implemented using an open-source software with GNU licenses (Version 3, Free Software Foundation, Boston, MA, USA) which can be adapted easily by small companies without access to expensive licenses. Other commercial software such as ANSYS or COMSOL can be employed to implement the proposed optimization method of stress sensors. In this investigation, the following software were used under the LINUX MINT [44] operating system-Octave 5.2.0 [45] (open-source MATLAB alternative), NETGEN [46], Elmer FEM [47] and ParaView [48]-to establish a modeling environment.
A schematic block diagram of the software interoperation is presented in Figure 1. ElmerGrid and ElmerSolver are modules of Elmer FEM software (Version 9.0, CSC-IT Center for Science, Espoo, Finland) which can be used independently. The geometry model is described in the NETGEN batch .geo file. The mesh is then computed using the Delaunay method and processed in ELMERGrid from a .gmsh file into Elmer mesh files. It should be stated that information exchange among modules is carried out using text files. These tools are beneficial for simplifying the troubleshooting and validation stages of this process.
The tetrahedral mesh of a magneto-mechanical system is presented in Figure 2. This mesh consists of a modelled shaft (1), a driving coil (2), and a sensing coil (3). The dimensions of the bounding box of the shaft with driving and sensing coils are x = y = 7.5 mm and z = 64 mm. The whole model was modelled in a sphere of radius r = 1 m. The mesh was built from 1,454,671 volume elements, 93,468 surface elements, and 140,202 edge elements in total. The sphere was assigned the properties of air at room temperature, and the copper and shaft coils were assigned the properties of structural steel. The relative magnetic permeability tensorμ r was calculated by using the procedure written in the FORTRAN programming language. The current density in the magnetizing coil and the torque load applied to the shaft were implemented within a .sif file in the MATC language as a source and boundary conditions for the model. The tetrahedral mesh of a magneto-mechanical system is presented in Figure 2. This mesh consists of a modelled shaft (1), a driving coil (2), and a sensing coil (3). The dimensions of the bounding box of the shaft with driving and sensing coils are x = y = 7.5 mm and z = 64 mm. The whole model was modelled in a sphere of radius r = 1 m. The mesh was built from 1,454,671 volume elements, 93,468 surface elements, and 140,202 edge elements in total. The sphere was assigned the properties of air at room temperature, and the copper and shaft coils were assigned the properties of structural steel. The relative magnetic permeability tensor was calculated by using the procedure written in the FORTRAN programming language. The current density in the magnetizing coil and the torque load applied to the shaft were implemented within a .sif file in the MATC language as a source and boundary conditions for the model.
The output signal on the real sensor is the voltage measured across a precise 1 kΩ resistor that is connected to the sensing coil. In the FEM modeling, the voltage was calculated from the base value of the current in the sensing coil while the resistance was modelled by measuring the resistance of a whole coil. This approach ensures the actual voltage drop across the resistor can be measured accurately.  The calculations were automated using scripts in Octave which generated batch .sif files and were run using the ElmerSolver. A total of 18 simulations were carried out for nine values of supply current with AC frequencies: 50 Hz, 250 Hz, 500 Hz, 1 kHz, 1.5 kHz, 2 kHz, 3 kHz, 4 kHz, and 5 kHz. The current amplitude used was 0.5 A for two values of torque load: 1.3 Nm and 0 Nm (no load), respectively. The results were then visualized using the ParaView software. Figure 3 shows the vectors of magnetic field strength before ( Figure 3a) and after ( Figure 3b) applying torque when the driving coil current frequency is equal to 5 kHz, the amplitude is 0.5 A, and the number of turns of the driving coil is 150. It is clearly seen that the direction of the vectors which were previously parallel to the axis of the shaft have now become helically skewed. The same phenomenon can be observed for magnetic flux density vectors as shown in Figure 4 where the driving current frequency, amplitude, and number of turns of the driving coil are set to 5 kHz, 0.5 A, and 200 respectively. The output signal on the real sensor is the voltage measured across a precise 1 kΩ resistor that is connected to the sensing coil. In the FEM modeling, the voltage was calculated from the base value of the current in the sensing coil while the resistance was modelled by measuring the resistance of a whole coil. This approach ensures the actual voltage drop across the resistor can be measured accurately.
The calculations were automated using scripts in Octave which generated batch .sif files and were run using the ElmerSolver. A total of 18 simulations were carried out for nine values of supply current with AC frequencies: 50 Hz, 250 Hz, 500 Hz, 1 kHz, 1.5 kHz, 2 kHz, 3 kHz, 4 kHz, and 5 kHz. The current amplitude used was 0.5 A for two values of torque load: 1.3 Nm and 0 Nm (no load), respectively. The results were then visualized using the ParaView software. Figure 3 shows the vectors of magnetic field strength H before ( Figure 3a) and after ( Figure 3b) applying torque when the driving coil current frequency is equal to 5 kHz, the amplitude is 0.5 A, and the number of turns of the driving coil is 150. It is clearly seen that the direction of the vectors which were previously parallel to the axis of the shaft have now become helically skewed. The same phenomenon can be observed for magnetic flux density vectors B as shown in Figure 4 where the driving current frequency, amplitude, and number of turns of the driving coil are set to 5 kHz, 0.5 A, and 200 respectively.
The calculations were automated using scripts in Octave which generated batch .sif files and were run using the ElmerSolver. A total of 18 simulations were carried out for nine values of supply current with AC frequencies: 50 Hz, 250 Hz, 500 Hz, 1 kHz, 1.5 kHz, 2 kHz, 3 kHz, 4 kHz, and 5 kHz. The current amplitude used was 0.5 A for two values of torque load: 1.3 Nm and 0 Nm (no load), respectively. The results were then visualized using the ParaView software. Figure 3 shows the vectors of magnetic field strength before ( Figure 3a) and after (Figure 3b) applying torque when the driving coil current frequency is equal to 5 kHz, the amplitude is 0.5 A, and the number of turns of the driving coil is 150. It is clearly seen that the direction of the vectors which were previously parallel to the axis of the shaft have now become helically skewed. The same phenomenon can be observed for magnetic flux density vectors as shown in Figure 4 where the driving current frequency, amplitude, and number of turns of the driving coil are set to 5 kHz, 0.5 A, and 200 respectively. The skewness of the magnetic field strength and flux density is connected with the direction of the easy axes of the stress induced anisotropy tensorμ r . These easy axes are observed by the rotation matrix R which was available during the Elmer Solver modeling. It should be highlighted that the observed results of the magnetic field strength and flux density modeling agree with previously presented experimental results in references [29,49]. Figure 5a,b show a cross-section view of the eddy current densities j before and after applying torque of 1.3 Nm. Similar to the modeling conditions presented in Figures 3 and 4, the driving current frequency, amplitude, and number of turns of the driving coil are 5 kHz, 0.5 A, and 150, respectively.
Before applying any torque, the current density vectors are tangent to the circle defining the circumference of the shaft and are contained in a plane perpendicular to the shaft axis. After applying a torque of 1.3 Nm, the eddy currents density vectors appearing on the shaft's surface have a greater value and are tangent to the thread line. Looking at the whole volume excluding the surface, there are smaller vectors with a direction that is parallel to the shaft axis running in the opposite way. This shows how the current circulates on the surface and returns to the center of the shaft.   It should be highlighted that eddy current flow in the parallel direction of the shaft subjected to torque has never been observed before. In spite of the fact that the amplitude of such parallel eddy currents is relatively small, this phenomenon may be practically utilized in the area of non-destructive testing of conductive shafts subjected to torque and shear stresses. Material discontinuities inside of a shaft may influence parallel eddy current distribution, which may be observed utilizing an eddy current tomography [50]. In such a case, the model presented in Section 2 creates the possibility of quantitative description of eddy current distribution and efficient estimation of a size and character of possible discontinuities in the core of a driving shaft subjected to torque.

(b)
The cause of all the results presented in Figures 3-5 is the influence of the complex stress distribution introduced by the torsion on the magnetic permeabilityμ r tensor, which changes the anisotropy direction of the material [36]. For further validation, the source code for modeling the magnetoelastic torque sensors is available at Supplementary Materials: https://github.com/ostaszewska/ElmerFEM_Shaft_torque under the open-source MIT license. observed by the rotation matrix which was available during the Elmer Solver modeling. It should be highlighted that the observed results of the magnetic field strength and flux density modeling agree with previously presented experimental results in references [29,49]. Figure 5a,b show a cross-section view of the eddy current densities before and after applying torque of 1.3 Nm. Similar to the modeling conditions presented in Figures  3 and 4, the driving current frequency, amplitude, and number of turns of the driving coil are 5 kHz, 0.5 A, and 150, respectively. Before applying any torque, the current density vectors are tangent to the circle defining the circumference of the shaft and are contained in a plane perpendicular to the shaft axis. After applying a torque of 1.3 Nm, the eddy currents density vectors appearing on the shaft's surface have a greater value and are tangent to the thread line. Looking at the whole volume excluding the surface, there are smaller vectors with a direction that is parallel to the shaft axis running in the opposite way. This shows how the current circulates on the surface and returns to the center of the shaft.
It should be highlighted that eddy current flow in the parallel direction of the shaft subjected to torque has never been observed before. In spite of the fact that the amplitude of such parallel eddy currents is relatively small, this phenomenon may be practically

Experimental Validation of Proposed Model
The experimental validation of the proposed model was performed on a 7250R drive shaft (Traxxas, McKinney, TX, USA) which is made of constructional steel. This type of shaft is generally used for electric vehicle models on a 1:16 scale. The geometry and dimensions of the shaft are shown in Figure 6a. For the purpose of this experiment, two coils were wound on the shaft: the magnetizing coil (where the length of the entire shaft is 150 turns) and a sensing coil of 50 turns. The coil arrangement diagram is presented in Figure 6b where numbers '1', '2', and '3' denote the shaft, magnetizing coil, and sensing coil, respectively, and the given dimensions are in millimeters. One end of the shaft was anchored, and torque was applied to the other end by hanging weights on a lever arm.
All measurements were done with the measurement system setup which is presented in the schematic block diagram shown in Figure 7. The experimental setup consists of the following: shaft is generally used for electric vehicle models on a 1:16 scale. The geometry and dimensions of the shaft are shown in Figure 6a. For the purpose of this experiment, two coils were wound on the shaft: the magnetizing coil (where the length of the entire shaft is 150 turns) and a sensing coil of 50 turns. The coil arrangement diagram is presented in Figure 6b where numbers '1', '2', and '3' denote the shaft, magnetizing coil, and sensing coil, respectively, and the given dimensions are in millimeters. One end of the shaft was anchored, and torque was applied to the other end by hanging weights on a lever arm. All measurements were done with the measurement system setup which is presented in the schematic block diagram shown in Figure 7. The experimental setup consists of the following: A dedicated software was installed on the PC computer to control the signal generator and collect data simultaneously in the LabView environment (National Instruments, Austin, TX, USA).
The experiment was carried out using 33 frequency values ranging from 50 Hz up to 5 kHz that were increased gradually in an adaptive pattern. The amplitude of driving current was constant and set to 0.5 A Peak . For each frequency change, the voltage induced on the sensing coil was measured for the shaft without torque (U detM0 (f )) and when the torque value was equal to 1.3 Nm (U detM1.3 (f )). Two values of torque M were selected to linearize the torque dependence of the output signal and to simplify the model for phenomenon explanation. The sensitivity value S(f ) for each frequency was calculated according to the following equation: (10) where: and In Equation (12), ∆U det_max is the maximum computed value of ∆U det ( f ) and ∆U det_min is the minimum value. Figure 8 presents the experimental results obtained for the sensitivity value S(f ) in comparison to the FEM modeling, while Table 1 provides detailed experimental results.
The output voltage and flux density characteristics were similar, which is due to low magnetizing fields used, and thus the steel shaft was working in the initial part of the magnetization curve. For practical engineering applications, the output voltage sensitivity characteristic is further discussed.
It should be highlighted that stress dependence of magnetic properties of materials (such as magnetic permeability tensor) is nonlinear for a wide range of axial or shear stresses. This nonlinearity is the main reason for differences between the results of experiments and modeling presented in Figure 8. Moreover, the decrease in magnetoelastic sensitivity for higher frequencies is connected with the strong influence of eddy current skin effect in conductive bulk elements made of constructional steel. This phenomenon is especially visible in the case of real shafts. For driving frequencies exceeding 5 kHz, the sensitivity of the sensor radically degrades due to the roughness of a surface of material and non-uniform distribution of surface stresses. The experiment was carried out using 33 frequency values ranging from 50 Hz up to 5 kHz that were increased gradually in an adaptive pattern. The amplitude of driving current was constant and set to 0.5 APeak. For each frequency change, the voltage induced on the sensing coil was measured for the shaft without torque (UdetM0(f)) and when the torque value was equal to 1.3 Nm (UdetM1.3(f)). Two values of torque M were selected to linearize the torque dependence of the output signal and to simplify the model for phe- On the other hand, in both scenarios, the maximum sensitivity is observed to be approximately 1.2 kHz. As a result, in determining the optimal sensor driving frequency, important similarities can be observed between the obtained experimental results and the proposed FEM-based modeling results.  It should be highlighted that stress dependence of magnetic properties of materials (such as magnetic permeability tensor) is nonlinear for a wide range of axial or shear stresses. This nonlinearity is the main reason for differences between the results of experiments and modeling presented in Figure 8. Moreover, the decrease in magnetoelastic sensitivity for higher frequencies is connected with the strong influence of eddy current skin effect in conductive bulk elements made of constructional steel. This phenomenon is especially visible in the case of real shafts. For driving frequencies exceeding 5 kHz, the sensitivity of the sensor radically degrades due to the roughness of a surface of material and non-uniform distribution of surface stresses.
On the other hand, in both scenarios, the maximum sensitivity is observed to be approximately 1.2 kHz. As a result, in determining the optimal sensor driving frequency, important similarities can be observed between the obtained experimental results and the proposed FEM-based modeling results.

Conclusions
The presented results in this paper have confirmed the possibility of modeling the torque dependence of magnetization processes in steel shafts. Considering previously presented work on the generalization of stress dependence of a magnetic permeability tensor, the proposed generalization represents the dynamic magnetization process that is connected to the eddy currents in a magnetized shaft. From modeling, it is achievable to observe the helical anisotropy in the shaft that is subjected to torque and the axial eddy currents in this shaft. It should be highlighted that research on the axial component in eddy current distribution in magnetoelastic shafts that are subjected to torque has not been presented previously in the literature.
The FEM-based modeling results have confirmed the prospect of determining the optimal driving frequency for magnetoelastic torque sensors with a core made of construction steel. The modeling and experimental test results both showed an optimal frequency of approximately 1.2 kHz. It is noted, however, that the optimal driving frequency of magnetoelastic sensors may be strongly dependent on its geometry and the type of steel used for the shaft. Therefore, modeling should be carried out specifically for any given material using the specific geometry of the magnetoelastic sensing shaft.
In future works, better modelling accuracy could be achieved with the use of an experimentally obtained model of the Anhysteretic Magnetization curve for given material, due to the recent achievements in this area [51].