A new methodology for incorporating the cutting deterioration of electrical sheets into electromagnetic finite-element simulation

This article proposes a novel methodology for incorporating electrical sheet cutting deterioration in electromagnetic finite-element simulations of energy conversion devices. While the existing methods account for the deterioration in the numerical integration either by increasing the mesh refinement or boosting the Gaussian quadrature order, the proposed method is based on the re-computation of quadrature weights and coordinates for a modeled deterioration, taking its explicit dependency into account. To validate the proposed method, numerical solutions are compared with electromagnetic analytical solutions in a beam geometry. A comprehensive analysis is then performed to evaluate the relative error, considering various model parameters. This analysis leads to a systematic procedure for selecting the optimal element size to achieve desired error levels. The procedure is successfully applied to a transformer geometry, and the computational performance of the proposed method is compared with the existing approaches through a time-stepping analysis. The results show that the proposed method is computationally more efficient than the existing approaches, and it eliminates the need to increase the mesh refinement or boost the order of the quadrature. It can be easily adapted for any type of deterioration profile.


Introduction
The use of electrical steel sheets in electromechanical energy conversion devices requires cutting the sheets into appropriate shapes by different techniques (e.g., punching, laser-cutting, waterjet cutting, and electrical discharge machining), which deteriorates the magnetization properties and increases the iron losses [1,2].To accurately simulate the energy conversion devices, magnetization, and iron loss models that include the cutting deterioration should be developed and identified based on experimental results and implemented into the finite-element (FE) simulation tools in a correct manner.
Continuous local material modeling is one of the popular approaches to model the deterioration of the magnetization of electrical steel sheets with a deterioration profile beginning from the cut edge towards the middle of the material.This profile was often modeled using exponential [3][4][5] or quadratic functions [6,7] parameterized as a function of the distance from the cut edge.A similar approach was used to include the effect of deterioration in the loss coefficients of empirical formulas (i.e., Bertotti's method [8] and Jordan's method [9]) [5,6,10,11].So far, successful results have been obtained in the material modeling stage with the use of these approaches compared to experimental findings.
While the continuous iron loss models accounting for the cutting deterioration have been commonly used at the post-processing stage of the FE simulations, there are two primary approaches for implementing the magnetization properties: (i) dividing the simulated geometry into smaller sections with distinct magnetic properties and (ii) modeling the magnetic properties as a continuous function of distance from the cut edge.While the first approach is easier to implement, the second approach offers a more accurate physical representation.Despite its advantages, implementing this approach poses several challenges, including numerical integration, appropriate selection of FE order, and increased computational time.
In the first approach, multiple regions or layers are created within the machine geometry and characterized by unique magnetization curves.This method was applied in the simulation of an induction machine, resulting in 5 regions in both the stator and rotor, as reported in [7].A similar approach was adopted in [2] for the simulation of a capacitor motor and a synchronous generator, though the number of layers was not specified.In [12], FE simulations were performed for a permanent magnet-assisted synchronous reluctance motor, including the effect of cutting on the magnetization of the stator material, using 3 distinct regions.Moreover, in [13], the effect of cutting on a permanent https://doi.org/10.1016/j.jmmm.2024.171843Received 2 May 2023; Received in revised form 14 November 2023; Accepted 9 February 2024 magnet synchronous motor was analyzed through FE simulations dividing the stator into 6 layers and applying specific magnetization curves to each layer.None of these studies provided information on the order of the elements or the selection of Gauss integration points.
The second approach to FE simulation involves the calculation of the distance to the cut edge, which is then used to obtain the local magnetization characteristics for each integration point in a continuous manner.To achieve enough accuracy, very fine meshes were used near the deteriorated edges with first-order elements in [11] and secondorder elements in [14], without further details on how the integration points and weights were selected.Due to the computational burden of using very fine meshes near the cut edge for certain geometries, the use of higher order elements with coarser meshes was proposed in [15] to reduce the computational time of the FE simulation of an induction motor with the inclusion of exponential deterioration profile.It was found that second-order polynomials were suitable for representing the flux density, and third-order elements were chosen accordingly.The use of an adequate number of integration points was mentioned, but specific details were not given.Better computational speed was achieved in [4] through the use of mixed-order elements.The use of the ''same'' number of Gauss integration points for all elements was mentioned, but no additional information was provided.
The literature review shows that there have been several studies on the implementation of cutting deterioration into FE machine simulation using either discrete modeling based on region subdivision or continuous modeling.Nevertheless, in these studies, the selection of the integration points for the quadrature has not been detailed, and the spatial dependency of the deterioration term on the numerical integration has been ignored.However, the deterioration term, whether exponential or quadratic, has an explicit dependency on the space and thus alters the form of the function to be integrated.Therefore, this dependency should be taken into account in the quadrature to ensure more precise numerical integration in a computationally efficient manner.
In this article, we propose a novel methodology for incorporating exponential deterioration in numerical integration based on the recomputation of the quadrature weights and coordinates for the modeled deterioration profile.The proposed method is first validated through comparisons with analytical solutions for a linear magnetostatic problem in a beam geometry.Comprehensive analyses are then conducted to examine the accuracy of the iron loss calculation, taking into account various parameters.In light of these analyses, a systematic procedure for selecting the optimal element size to achieve the desired error level is developed.This procedure is then applied to the 2-D FE simulation of a nonlinear problem in a transformer geometry.Circuit equations in the windings are coupled with the magnetostatic problem in the transformer core and then solved with a time-stepping analysis, where the computational efficiency of the proposed method compared to the existing ones is investigated.The proposed method is analyzed and validated for an exponential deterioration profile, but it is applicable to any form of deterioration profile and cutting method, with the need to identify the appropriate material coefficients.

Material modeling with cutting deterioration
In this paper, to facilitate the proposed incorporation methodology, we adopted continuous permeability and loss modeling approaches found in the existing literature.Initially, following the continuous local material modeling approach in [7], we express the single-valued reluctivity  as a function of magnetic flux density norm  and the shortest distance from the cut edges (, ) using two reluctivity curves, one for the case where the material would be fully undamaged  un () and one for the case where the material would be fully damaged  dam () such that  [17].(b) Measured data and modeled reluctivity curve for  un () and  dam () with the deterioration profile  −∕ r .The modeled reluctivity matches successfully with the measurements in [17] with an absolute relative error of less than 5%.
where  un (),  dam (), and  r are model parameters to be identified based on the experimental results.In line with the methodologies in [4,5], we used an exponential deterioration profile  −∕ r and assumed a continuous behavior throughout the sample for practical considerations.This approach removes the necessity of identifying a specific degradation depth, which has exhibited diverse treatment across various studies.The reluctivity curves  un () and  dam () in (1) are parameterized by Marrocco's equation proposed in [16], which expresses () as Here,  = [ 1  2  3  4 ] T is a vector of model parameters that needs to be identified.Consequently, the vectors of model parameters  un and  dam need to be identified for  un () and  dam (), respectively.These can be determined through a fitting procedure using experimental results of samples exposed to different levels of cutting deterioration.In this paper, we use the experimental results in [17], identify  un from the measurements of the EDM-cut sample, and use  un and the measurement results of the other samples to identify  dam and  r .The value of  r is identified as 1∕640 in m and the vector of model parameters obtained are given in Table 1.Geometries of the measured samples, the modeled reluctivity curves for these samples, and their comparison with the measured data are shown in Fig. 1.
Taking inspiration from the robustness and mathematical practicality of the approach introduced in [5], similar to the magnetization model used in Section 2, we express the total iron loss density  Fe under the sinusoidal excitation of frequency  and amplitude of flux density and  dy () in W kg −1 Hz −2 T −2 .(b) The measured and modeled losses at 10 Hz and 100 Hz with the identified coefficients.The modeled losses match successfully with the measurements in [17] with an absolute relative error of less than 8.5% for each sample at 10 Hz, 25 Hz, 50 Hz, and 100 Hz.
m using Jordan's two-component method with hysteresis and dynamic loss coefficients  hy () and  dy () [9]: where  hy,un ,  dy,un ,  hy,dam ,  dy,dam ,  hy , and  dy are model parameters need to be identified based on the experimental results.It should be noted here that although the dynamic loss components might not explicitly depend on the cutting, the interaction between the hysteresis and the classical eddy-current affects these components.
Similar to the identification procedure for the reluctivity model, we identify  hy,un and  dy,un from the measurements of the EDM-cut sample, and use those and the measurement results of the other samples to identify  hy,dam ,  dy,dam ,  hy , and  dy .Fig. 2 shows the identified coefficients and the modeled losses based on these coefficients in comparison with the measurement results.

Formulation
In a 2-D FE simulation, the Galerkin-discretized [18] electromagnetic field problem can be expressed as where  is the stiffness matrix,  contains the nodal vector potentials, and  is the load vector.Including the effect of cutting in the reluctivity using (1), the entries of the stiffness matrix   with the nodal shape functions   and   in domain  become which can be represented in two parts as The correct calculation of   requires the proper numerical integration of the  1, and  2, terms for all domains.As the reluctivity term in  1, does not depend explicitly on coordinates, the numerical integration can be accurately performed with 2-D Gaussian quadrature for triangular elements with the well-known weights and coordinates [19].However, the exponential weighting function in  2, has an explicit spatial dependency.Therefore, the same Gaussian quadrature used for the integration of  1, cannot yield an accurate result for  2, .
The weights and coordinates of the integration points should be recomputed taking the exponential deterioration function into account.
To address this issue, we propose a new approach, which will be detailed in the upcoming part.

Numerical integration with exponential deterioration
In a triangular domain  in  reference coordinate system with  =  ∶ 0 ≤ , ,  +  ≤ 1}, express the quadrature rule the integration of function  (, ) weighted with an exponential deterioration function  −(,)∕ r such that where  int is the number of integration points,   are the weights, and (  ,   ) are the coordinates of the integration points in the reference element.By definition, the quadrature should be exact for all polynomials in the complete polynomial space for degree  pol {    , 0 ≤ , ,  +  ≤  pol }, when the polynomials replace  (, ) [19].Following this definition, a nonlinear system of equations can be built by replacing  (, ) in ( 8) with each of polynomials in the complete polynomial space  1 (, ),  2 (, ), … ,   (, ) with  = ( pol+1 )( pol +2)∕2 as follows: If the terms on the left-hand side are computed accurately, the resulting nonlinear equation system can be solved to obtain  int ,   , and (  ,   ) for domain .Following this property, in the pre-computation stage of the FE simulations, we build the complete nonlinear system of equations in (9) for each domain of a meshed geometry by computing the left-hand side integrals numerically with an adaptive quadrature method.We then solve (9) by using a nonlinear solver based on Levenberg-Marquardt algorithm by keeping  int the same as in the Gaussian quadrature in [19].Next, we store the computed weights and coordinates, and then compute and store the shape functions and their derivatives at the computed coordinates for each domain.The entire process can be done in parallel in the pre-computation stage of the FE solution.
During the FE solution, we perform the numerical integration of  1, with the Gaussian quadrature weights and coordinates, and shape functions and their derivatives at those coordinates.We perform the numerical integration of  2, using the re-computed and stored weights and coordinates, and shape functions and their derivatives for the corresponding domain with the same  int as in  1, .Thus, for each FE domain, including the mapping between  reference coordinate system and  global coordinate system, the calculation of  2, is obtained as follows: where and  ′  ,  ′  are the shape functions at the re-computed coordinates.In the remainder of the paper, we adopt the following terminology: • Proposed method: The numerical integration of  1, with the Gaussian quadrature weights and coordinates, and shape functions and their derivatives at those coordinates; and the numerical integration of  2, with the re-computed and stored weights and coordinates, and shape functions and their derivatives at recomputed coordinates in the pre-computation stage.The same terminology is also used for the numerical integration of similar forms.
• Classical method: The numerical integration of  1, and  2, in a combined form as   in ( 5) with the Gaussian quadrature weights and coordinates, and shape functions and their derivatives at those coordinates.The same terminology is also used for the numerical integration of similar forms.

Applications and results
In this section, we apply the proposed method to two different geometries and compare it against the application of the classical method.A short description of the applications in the studied geometries is as follows: (i) Beam geometry: We model the deterioration in the magnetization using a linear material law in a magnetostatic case.We derive an analytical solution and use it as a reference to evaluate the accuracy of the FE solutions and numerical integration methods.(ii) Transformer: We model deterioration in both the magnetization and iron losses using a nonlinear material law in 2-D FE simulations.We select the element size based on the results from the beam geometry and couple the magnetostatic problem in the transformer core with circuit equations in the windings.Then, we perform a time-stepping analysis to evaluate the accuracy of numerical integration methods in calculating iron losses and investigate the computational performance of these methods.

Beam geometry
A beam geometry with two deteriorated edges in the -plane is studied (Fig. 3(a)).The vectors of magnetic flux density  = ()  and magnetic field strength  = ()  are aligned in the -direction.Both edges of the beam are considered to be deteriorated by cutting and the minimum distance to the closest cut edge is used for (, ) in the material modeling (see (1)).

Analytical and numerical solutions
The magnetic flux  per thickness (Wb/m) is enforced with a predefined average flux density  p such that In the absence of current density  , Ampere's law is Using the magnetization model in (1) for the constitutive relationship between  and  for the linear material case (with a generic representation of decay constant  to distinguish it from the identified decay constant  r ), the differential equation in (13) becomes Due to the symmetry reasons, it is sufficient to solve (14) only for the positive region of the geometry (0 <  < ), for which the analytical solution can be expressed as where  c is a constant and can be obtained by inserting (15) into (12) and solving for  c .The analytical solution is exact and considered to be a reference for the comparisons.
For the numerical solution, a homogeneous Dirichlet condition is applied on the deteriorated edges by setting the nodal vector potentials to ( = ±) = ±∕2.The problem in Fig. 3(a) with the dimensions ℎ =  = 10 mm is simulated for meshes with both second-order and third-order triangular elements, labeled as  el = 2 and  el = 3, for a wide range of element sizes  size = [ ∕8,  ] (see Fig. 3(b) for  size = ∕8) and a wide range of decay constants  = [ 1∕10 4 , 1∕10 ] m.With the proposed method, we keep the order of quadrature  quad at the minimum required order  quad = 2( el − 1) (i.e.,  quad = 2 for  el = 2 and  quad = 4 for  el = 3) and choose the corresponding number of integration points  int (i.e.,  int = 3 for  quad = 2 and  int = 6 for  quad = 3).With the classical method, we boost  quad gradually to reach to a similar accuracy as the proposed method, which is achieved for  quad = 8 and correspondingly  int = 16.
Fig. 4 shows an example comparison of the FE simulation results obtained using classical and proposed methods, with the analytical As Fig. 4 shows, the FE simulation results have a high level of accuracy compared to the analytical solution, which validates the correct implementation of the FE algorithm and the applicability of the proposed method.

Analysis
To determine how accurately the numerical integration method reflects on the increase of the iron losses in a deteriorated material due to the variation of the square of the RMS flux density  2 rms (as  Fe ∝  2 rms ), the variation from the square of the RMS flux density of the undamaged case  2  un for both the FE-simulated cases and the corresponding analytical solution is calculated.This is represented by the following equations: where  2 FE and  2 A are the squares of the RMS flux densities obtained from the FE and analytical solutions, respectively.Then, the relative error  between the FE-simulated and analytical solutions is calculated by Fig. 5 illustrates  calculated for each FE-simulated case, which is performed using both proposed and classical methods for the numerical integration.As the deterioration profile is a function of both decay constant and distance from the cut edge, the data in the plots is given per length , which is 10 mm for the simulated geometry in Fig. 3(a), and this allows easy adaptation to similar geometries with different .To show a range where the absolute relative error is acceptable, e.g.|| < 5%, red dashed lines are used in the plots correspondingly.Fig. 5 shows that for a fixed decay constant ∕,  decreases as the element size  size ∕ decreases.This is expected as denser meshes provide greater accuracy in the FE simulations.Additionally, as the decay constant ∕ decreases, the value of the required element size  size ∕ to reach the acceptable accuracy decreases.The reason is that as the decay constant ∕ decreases, the decay becomes steeper, which can be captured with the use of smaller elements.In the same manner, for a fixed element size  size ∕,  decreases as decay constant ∕ increases.
For a fixed decay constant ∕ and utilizing the same method for numerical integration, whether proposed or classical, FE simulations performed with the third order triangular elements ( el = 3) always yield better accuracy than those with the second order triangular elements ( el = 2) due to a better representation of the approximated quantities with the shape functions.It is observed that the same accuracy level can be reached with the use of much coarser meshes.For instance, with the proposed method, while the absolute relative error of less than 5% for the case ∕ = 1∕(2.5 × 10 1 ) is reached when  size ∕ = 1∕8 with  el = 2, the use of  size ∕ = 1∕2 is sufficient for the case with  el = 3.
For a fixed order of triangular elements  el and using the same order of quadrature  quad , the proposed method typically yields better accuracy than the classical method for the same simulated case.In most cases, the same level of accuracy can be achieved with the coarser meshes.For instance, when  el = 2 and  quad = 2, the proposed method reaches an absolute relative error of less than 5% for the case ∕ = 1∕(5 × 10 1 ) when  size ∕ = 1∕8, while the classical method requires the use of  size ∕ = 1∕16 to achieve the same accuracy.
For a fixed order of triangular elements  el and using the classical method, increasing the order of quadrature  quad and correspondingly the number of integration points  int increases the accuracy for the same simulated case.As the results show, the use of proposed method with for instance  el = 2 and  quad = 2 and the use of classical method with  el = 2 and  quad = 8 yield a similar level of accuracy, and require a similar  size ∕ to reach a specific accuracy.However, for  quad < 8, the proposed method predominantly yields better accuracy than the classical method for the same simulated case.

Optimal element size selection
Next, we investigate the threshold element size  size ∕ required to achieve || < 5% as a function of the decay constant ∕ for each simulated case in Fig. 5 as well as the identified material model decay constants in Section 2. To achieve this, in addition to the simulated cases in Fig. 5, simulations for the identified material model decay constants ( r = (1∕640) m,  hy = (1∕3600) m, and  dy = (1∕240) m) are performed.For a fixed decay constant ∕, the values of the largest element size  size ∕ where || < 5% and the smallest element size  size ∕ where || > 5% are then extracted.Afterward, these values are used in the interpolation to obtain the threshold  size ∕ where || = 5%.Exceptionally, if || = 5% is already reached when  size ∕ = 1, the threshold is also set to  size ∕ = 1.Fig. 6 illustrates the required element size  size ∕ as a function of the decay constant ∕.The identified model decay constants  r ,  hy , and  dy per  are shown with the red dashed lines.
Fig. 6 shows that for the second order triangular elements ( el = 2), both the classical method with  quad = 8 and the proposed method have similar requirements for the element size  size ∕ for any value of decay constant ∕.The classical method with  quad = 2 has similar requirements for  size ∕ as the others when ∕ > 3×10 −1 , but requires a smaller  size ∕ when ∕ ≤ 3 × 10 −1 .For the third order triangular elements ( el = 3), both the proposed method and the classical method with  quad = 8 have similar requirements for the element size  size ∕ when ∕ ≤ 2 × 10 −2 and ∕ ≥ 10 −1 .When 2 × 10 −2 < ∕ < 10 −1 , the proposed method has a higher threshold for the element size  size ∕, and thus requires a coarser mesh.The classical method with  quad = 2 has similar requirements for the element size  size ∕ as the others when ∕ > 2 × 10 −1 , but requires a denser mesh and a smaller element size  size ∕ when ∕ ≤ 2 × 10 −1 .
Although this whole analysis is performed with a linear material, the presented scenarios in Fig. 5 can be employed to select the element size for a given accuracy in various applications.In the next section, we will adopt this procedure for selecting the required element size for transformer applications.

Transformer
A transformer geometry in the -plane is studied through 2-D FE simulations (Fig. 7(a)).In these simulations, the deterioration in the  magnetization and iron losses are modeled using the presented models in Section 2. The dimensions of the simulated transformers are chosen such that ℎ in = 20 mm, ℎ out = 40 mm, and  = 10 mm.All edges of the transformer are considered to be deteriorated by cutting and the minimum distance to the closest cut edge is used for (, ) in the material modeling (see ( 1) and ( 3)).

Simulation results and analysis
In order to evaluate the accuracy of the numerical integration methods in the computation of iron losses and compare their computational performance, a time-stepping analysis is carried out whereby the circuit equations in the transformer windings are solved in a coupled manner with the magnetostatic problem in the magnetic core.Sinusoidal voltages with an amplitude of 240 √ 2 V and a frequency of 50 Hz are applied to the primary and secondary windings of the transformer, and the number of turns in each winding is adjusted to yield a current density  = 5 × 10 6 A∕m 2 .The input parameters used in the circuit equations of both primary and secondary windings are listed in Table 2.
The deterioration in the magnetization is included in the FE solution of the electromagnetic field problem using the formulation in (1) with the identified parameters in Fig. 1(b).The total iron losses  Fe are computed in a post-processing manner after obtaining the FE solution utilizing the iron loss model in (3) with the identified parameters in Fig. 2(a).In order to account for the contribution of each harmonic in the losses, Fourier harmonic analysis is performed for each element.Following a similar approach to [20,21], total hysteresis and dynamic losses  hy and  dy over the volume  are calculated such that where  s is the supply frequency,  is the harmonic number, and  m, is the peak value of the flux density norm for the th harmonic for the corresponding element.Each of the integrals is in a similar form as (5) and thus can be numerically integrated with both the proposed and classical methods.
To accurately estimate the iron losses within a specified error range, proper selection of  size is crucial in FE simulations.Proper selection of  size enables accurate numerical integration of each term with the deterioration decay constant parameters used in the modeling ( r = (1∕640) m,  hy = (1∕3600) m, and  dy = (1∕240) m).The most accurate estimation can be reached by first finding the required element size  size for each deterioration decay constant parameter and then choosing the minimum among those to account for the worst case.Following this property, for each simulated case, element size  size can be selected based on the corresponding threshold value for  hy = (1∕3600) m.The simulations are performed for both second and third-order triangular elements ( el = 2 and  el = 3) using the proposed and classical methods for the numerical integration within the FE simulation and also in the post-processing for the computation of the total iron losses over the volume  Fe .Table 3 shows the details of the simulated cases, referred to as Cases 1-6 in the upcoming parts, and the selected values of element sizes  size .Fig. 7(b) shows an example of the flux density distribution obtained from the simulations using the proposed method with  size = ∕6,  el = 2, and  quad = 2.
As there is no analytical solution for the nonlinear case to use as a reference, a very dense mesh with  size = ∕12 is simulated using the proposed method with the third order triangular elements ( el = 3). hy and  dy are computed at the post-processing stage and then considered as a reference, i.e,  hy,ref and  dy,ref .Similar results, with an absolute error of less than 0.1%, are obtained from the simulations of the same mesh discretization using the classical method with  el = 3 and  quad = 8, which verifies the feasibility of the reference selection.Then, absolute relative errors for hysteresis loss  hy and dynamic loss  dy for each simulated case in Table 3 are calculated as follows: where  hy,FE and  dy,FE denote the hysteresis and dynamic losses for the FE-simulated cases.Fig. 8 shows the calculated  hy and  dy .
As shown in Fig. 8,  hy and  dy are below 2% for all simulated cases, confirming the validity of the element size  size procedure presented in Section 4.1.It is observed that  hy is slightly higher than  dy , which is expected as  hy <  dy .The best accuracy is achieved when  el = 3.

Computational time analysis
The FE simulations are performed for all the cases in Table 3 using a time-stepping analysis for two supply periods, discretizing each period into 400 steps.The simulations are performed in MATLAB software using Intel ® Core TM i5-11600 @ 2.80 GHz processor with 32 GB RAM.The iterations are set to stop when the norm of the nodal vector potentials reaches below 1 × 10 −5 .During the simulations, CPU times for pre-computation, time-stepping, and post-processing stages are recorded.Table 4 shows the recorded CPU time data for each simulated case, denoted as Cases 1-6 (see Table 3 for the details of the cases).Each plot shows the calculations of the simulations conducted with both second-order and third-order triangular elements, labeled as  el = 2 and  el = 3.For the same triangular element order  el value, 3 different combinations of the quadrature order  quad and numerical integration method are simulated, which makes 6 cases in total (see Table 3).Table 4 shows that the majority of computational time is spent in the time-stepping stage.The pre-computation stage takes less than 1 s and can be considered to be negligible in the analysis of total computational time.For the simulations with  el = 2 (Cases 1-3), the CPU time recorded in the time-stepping stage of the cases using the classical method with  quad = 8 (Case 2) and the proposed method (Case 3) is significantly less than that of the case using the classical method with  quad = 2 (Case 1).Similarly, for the simulations with  el = 3 (Cases 4-6), the CPU time recorded in the time-stepping stage of the cases using the classical method with  quad = 8 (Case 5) and the proposed method (Case 6) is significantly less than that of the case using the classical method with  quad = 4 (Case 4).The reason behind the significant differences can be directly attributed to the increased degree of freedom due to denser meshes.
The use of the proposed method leads to a 9.9% and 12% less computational time in the time-stepping stage for the simulations with  el = 2 and  el = 3 respectively, compared to the use of the classical method with  quad = 8.Although the post-processing stage is longer when the proposed method is used, it is significantly shorter than the time-stepping stage and does not significantly impact the total computational time.As a result, the overall difference in the computational time becomes 9.4% and 10.4% for the simulations with  el = 2 and  el = 3 respectively.
Comparison of the simulated cases with  el = 2 and  el = 3 (i.e., Cases 1, 2, 3 compared to Cases 4, 5, 6, respectively) shows that using  el = 3 results in shorter time, despite increased degree of freedom that leads to longer computational for fixed  size .However, as  el = 2 requires a smaller  size shown in Table to achieve the specified accuracy, it results in higher computational time in the studied cases.

Discussions
The results of this study indicate that the utilization of the proposed method enabled the use of larger  size and thus coarser meshes in comparison to the classical method for the same values of  el and  quad (as shown in Fig. 6 and Table 3), while still achieving a predefined level of accuracy (as shown in Fig. 8).This led to a significant decrease in the computational time required for the time-stepping analysis conducted for the transformer application.In the simulated cases, 10-13 times difference was observed in the total computational time of a timestepping analysis, including the time spent in pre-computation and post-processing stages (as shown in Table 4).
For the same triangular element order  el , when the use of the proposed method is compared with to use of the classical method with boosted Gaussian quadrature order  quad = 8, it was seen that they required similar element size  size and thus similar mesh discretization (as shown in Fig. 6 and Table 3) to achieve a predefined level of accuracy (as shown in Fig. 8).Despite the use of the same  size , the time-stepping analyses for the transformer showed that the use of the proposed method resulted in up to 10.4% shorter total computational time for the simulated cases (as shown in Table 4) with the single-valued material law.
In addition to this difference in the computational time, there are two main drawbacks of boosting the Gaussian quadrature order  quad and corresponding number of integration points  int .The first drawback is that in the case of a hysteretic material law based on the hysteron decomposition, it results in a much poorer performance compared to the use of the proposed method in terms of both computational time and memory allocation.The second drawback is the uncertainty in the selection of the required quadrature order  quad .In the simulated cases of this study, a similar level of accuracy with the proposed method was reached by boosting the quadrature order  quad up to 8, but for any other study, this selection might be ambiguous.Nonetheless, with the proposed method, the minimum required quadrature order  quad based on the triangular element order  el can always be employed.

Conclusion
In this paper, a novel methodology to incorporate the cutting deterioration in FE simulation is proposed.The findings of the study have proved the accuracy of the proposed method as well as its computational efficiency compared to the existing methods.While this study focused on incorporating exponential deterioration, the method can be adapted to any type of deterioration profile and cutting method.Additionally, the presented systematic approach for selecting the required element size to reach a specific predefined level of accuracy is expected to provide a perspective for future studies.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 2 .
Fig. 2. (a) Identified hysteresis and dynamic loss coefficients  hy () in W kg −1 Hz −1 T −2 and  dy () in W kg −1 Hz −2 T −2 .(b)The measured and modeled losses at 10 Hz and 100 Hz with the identified coefficients.The modeled losses match successfully with the measurements in[17] with an absolute relative error of less than 8.5% for each sample at 10 Hz, 25 Hz, 50 Hz, and 100 Hz.

Fig. 3 .Fig. 4 .
Fig. 3. (a) Studied beam geometry in the -plane with the deteriorated edges along the  axis.(b) Mesh of the simulated beam geometry for ℎ =  = 10 mm and  size = ∕8.

Fig. 5 .
Fig. 5.The figure illustrates the relative error, denoted by , between the analytical reference solution and numerical solution in the calculation of the variation of  2 m from the undamaged case for the beam geometry in Fig. 3.The simulations are obtained using both (a) second-order and (b) third-order triangular elements, labeled as  el = 2 and  el = 3.For both  el = 2 and  el = 3, the results are obtained for a wide range of ∕ and  size ∕ values using proposed and classical methods, with the latter utilizing different orders of quadrature  quad and corresponding number of integration points  int .The used values of  el ,  quad , and  int are indicated inside the figures.The red dashed lines inside the figures indicate where  = ±5%.

Fig. 7 .
Fig. 7. (a) Simulated 2-D transformer geometry.ℎ in = 20 mm, ℎ out = 40 mm, and  = 10 mm.It should be noted that all edges of the transformer are assumed to be deteriorated by cutting.(b) Flux density distribution obtained from the FE simulation using proposed method with  el = 2,  quad = 2, and  size = ∕6.

Fig. 8 .
Fig.8.Calculated  hy and  dy .Each plot shows the calculations of the simulations conducted with both second-order and third-order triangular elements, labeled as  el = 2 and  el = 3.For the same triangular element order  el value, 3 different combinations of the quadrature order  quad and numerical integration method are simulated, which makes 6 cases in total (see Table3).

Table 1
Identified model parameters of the reluctivity curves.

Table 2
Parameters used in the circuit equations.

Table 3
Details of the simulated cases.

Table 4
Details of the CPU time data of the simulated cases.