Identification of representative anisotropic material properties accounting for friction and preloading effects: A contribution for the modeling of structural dynamics of electric motor stators

Simulating the dynamic behavior and determining equivalent material properties for anisotropic models, superelements or structures subjected to preloads or friction remains a challenging issue. Amongst other practical applications, modeling interactions between the steel sheets in industrial magnetic cores of electric motor stators is a complex task, as it requires anticipating behavioral heterogeneities in the structure, and possibly represents significantly costly operations for performing modal or dynamic response simulations. In this article, a method for identifying equivalent material properties to anisotropic structures is developed, which is able to take into account the influence of preloads and friction on the material properties, later used in structural dynamics simulations. The proposed approach can be used with superelements, converting stiffness matrices into elasticity matrices. The method is first applied to a triclinic model, and recreates its elasticity matrix with little derivation. Then, an equivalent linear material is computed for a continuous structure under preloading. Compared at low frequencies, the vibration behavior of the preloaded structure and its equivalent effective media are in good agreement. The operation is repeated with a laminated stack under preloading. Again, the dynamic behavior of the equivalent structure shows good accuracy compared to the initial preloaded stack. Finally, the magnetic core of an electric machine stator is modeled with equivalent anisotropic material properties, accounting for friction and preload in the yoke's and the teeth's steel sheets. The simulation of the structure's low-frequency radial vibration modes is satisfying, and shows improvement compared to orthotropic properties.


Introduction
For finite-element simulations, modelling heterogeneous structures strictly as they are in reality is sometimes both delicate and unrealistic if the heterogeneities are small or numerous: representing such models by equivalent homogeneous material properties may be a necessity in order to reduce the number of degrees of freedom in the models. In addition to this, constraints of cost-effectiveness and flexibility justify a great interest of the industry for effective methods able to recreate the dynamic behaviours of heterogeneous structures with equivalent material properties [12]. Some of the existing techniques able to determine homogeneous properties in order to model heterogeneous structures are called "homogenisation" methods.
A large number of homogenisation methods already exist in the cases of 2D, laminated, honeycomb or many other types of composite structures. A thorough review has been made by Kalamkarov et al. [14]. As for simulations on structures where simplifications are not possible, a 3-D approach needs to be preferred. This article focuses on this case.
There exist some cases in which equivalent isotropic or orthotropic materials are not accurate enough to recreate the behaviour of a given structure, with yet the same necessity of using representative elasticity matrices. In such a case, identifying elastic constants such as Young's moduli, shear moduli or Poisson's ratios is not possible in the case of anisotropic, or so-called "triclinic" structures, defined by 21 independent constants [3]. In this general case, the material's 21 independent coefficients form the elasticity matrix [C] in Hooke's law {σ } = [C] {ε}, where {σ } is the stress tensor and {ε} the strain tensor [5]. The entire linear system is: where the indices 1, 2 and 3 correspond to the respective directions x, y and z in a rectangular coordinate system, or r, θ and z in a cylindrical coordinate system.
Begis et al. [2] have developed a homogenisation method that can be applied to triclinic periodic structures, the principles of which were taken as a reference for other techniques. Although it yields equivalent elasticity matrices, the approach requires solving analytical equations that might make them delicate to implement for finite-element analyses (difficulty reported by Chung et al. [6]). A different analytical formulation of elastic coefficients has been proposed by Luciano and Barbero [16]. In spite of the ability to model equivalent anisotropic material properties, the expressions of only 6 coefficientsC i j amongst the elasticity matrix's 21 independent constants are detailed. Mathan and Siva Prasad [17] have developed a method of evaluation of equivalent material properties for a spiral-wound gasket and analysed its elasto-plastic behaviour. The principle is to average the stress-strain behaviour of a representative volume element over its volume by the means of independent load cases, in order to determine the equivalent compliance matrix's constants. However, the method is applied to the studied gasket only, and the structure itself is changed (angle of sealing ring) to ensure the independence of the load cases before averaging the results.
In some situations, the material properties of given structures can be affected by perturbations such as preload or friction. Therefore, modelling such a structure with equivalent properties requires taking the effects of the perturbations into account in the homogenisation or the identification process. Peillex et al. [20] successfully homogenised a model made of 10% of isotropic heterogeneities into an isotropic matrix, taking into account dynamic frictional contact conditions. By coupling their method to homogenised local friction coefficients, they were able to fairly approximate the stresses present in a heterogeneous model and their evolution through time. Smit et al. [22] as well as Yvonnet et al. [28] have tackled similar issues and taken into account possible non-linearities in their models. Yet both approaches analyse the behaviour of a non-linear, heterogeneous structure without determining any equivalent elasticity matrix, and no finite-element analyses have been performed to recreate the samples' mechanical behaviours.
In the field of laminates, Pirnat et al. [21] have developed a numerical model of a laminated stack with both tangential and normal contact conditions between the layers. The 14%-error in the prediction of the natural frequencies of a laminated stack was considered good by the authors, in addition to the fact that they also updated the interlayer contact parameters with experimental modal data. Focusing on electric machines, Kim and Kim [15] have shown experimentally that the first natural frequencies of a laminated rotor increased with the stack's pressure (i.e. preloading). Similarly, the works of Watanabe et al. [25] have outlined the fact that the natural frequencies of nonpurely-radial modes on a segmented-core stator tend to increase as the clamping force increases. Dias' experimental analyses [9] have confirmed this conclusion, with the observation that purely-radial (or "ovalisation") modes were not affected by any variations of the clamping force held by the weld beads (called tie rods). However, in spite of the new opportunities they represent to the knowledge in laminated structures' vibratory behaviours, these approaches are highly dependent on experimental data or model updating procedures and do not lead to any elasticity matrices that could represent the material properties of equivalent homogeneous structures.
The influence of preloads has also been analysed in the frame of railway dynamics. The simulations performed by Wu and Thomson [27] highlighted a strong influence of preloading on the studied railway's vibrations. The measurements of Kaewunruen and Remennikov [13] on railpads led to similar observations concerning the dynamic stiffness of all types of pads.
As for industrial projects in structural dynamics, a great interest is shown for accurate and cost-efficient modelling techniques applicable to common finite-element simulations such as modal bases. In particular, the automotive industry currently focuses on new technologies such as hybrid or 100%-electric powertrains, the stators of which are built on multi-layered magnetic cores with lateral weld beads [10,24,26]. Performing dynamic simulations on such heterogeneous structures is however highly dependent on costly and time-consuming experimental analyses. This is why understanding and being able to predict the vibratory behaviour of a stator's laminated magnetic core without experimental data is key to performing efficient noise and vibration simulations on entire electric motors. In addition to this, Van der Giet et al. [23] reported non-linearities in the mechanical behaviour of laminated stacks, as a result of their experimental analysis on a real electric machine.
Concerning finite-element simulations, Verma et al. [24] and Williams et al. [26] suggested to model the entire magnetic core with a single homogeneous isotropic material, what was applied by Dias [9]. An orthotropic material property identification method and modelling guidelines have been presented by Millithaler et al. [18] and allowed to simulate the modal behaviour of an electric motor stator with good accuracy, with an efficient zoning of the finiteelement model accounting for the influence of weld beads. The work presented by Millithaler et al. [18] also showed that the correlation accuracy between simulated and measured ovalisation modes on a real stator would have been decreased by a relative factor of 68.5% if the finite-element model had not been zoned. This article follows the work presented by Millithaler et al. [18], in which a method of 3D material property identification for multi-layered orthotropic laminates has been developed and where its applicability was compared to other existing 3D homogenisation techniques. In this paper, friction and preload will be taken into account phenomena in the computation of representative anisotropic materials for improving the accuracy of dynamic simulations on the magnetic core's finite-element model.
In order to compute representative elasticity matrices for modelling heterogeneous structures, several interesting methods have been reviewed. While some of them are able to model triclinic properties in finite-element models, others deal with perturbations and their non-linear influences on the behaviours of the studied structures. To the authors' knowledge, there currently exist no approaches able to take such external perturbations into account in the definition of equivalent, linear, representative elasticity matrices. Modelling such perturbations on a given structure generally requires performing non-linear simulations, which are in principle longer, more complicated and more restrictive (e.g. requiring special licences) than linear solutions. Approximating non-linear effects by representative linear material properties that could be used in any finite-element solution therefore offers interesting opportunities, that include analysing the influence of such effects on the material coefficients themselves, and integrating these material definitions into linear solutions such as model updating procedures. This is why a new method is proposed, based on finite elements, for modelling equivalent triclinic material properties for periodic structures and which is able to take into account the influence of perturbations on the representative elasticity matrix. Amongst the possible applications for this method, triclinic bodies as well as preloaded heterogeneous structures will be modelled with homogeneous models with equivalent material properties identified with the presented approach. Apart from the development itself, the main interest of this paper is to model multi-layered magnetic core bodies with representative homogeneous material properties that take into account the influence of preload and friction in the structure. As stated earlier, the idea is to approximate non-linear effects by representative linear material properties and allow performing simulations with linear models. Aiming at predicting the radial modes with best possible accuracy, the modal basis of the magnetic core's equivalent model will be compared to experimental data acquired on a real electric machine. This is why modal basis correlation criteria are preferred for all the proposed validation cases in this paper.

Development of the "Triclinic" method
As it has been explained by Millithaler et al. [18], orthotropic material properties imply that tension-compression and shear phenomena are not coupled. For a structure subjected to perturbations, these no-coupling assumptions may not be valid in the general case, and any linear elasticity matrix approximating its behaviour shall be therefore expressed as anisotropic, or "triclinic". The corresponding Hooke's law in the general case has been shown on Equation (1); in this case, the material can no longer be expressed with Young's moduli, shear moduli or Poisson's ratios [5]. Therefore, the method presented by Millithaler et al. [18], which aimed at identifying these elastic coefficients directly, is extended for application to triclinic materials and with effects of perturbations such as preload and friction.
First of all, the prerequisite steps for this method are: 1. The zones to be modelled with equivalent, homogeneous material properties are distinguished and analysed separately; 2. For a given zone, a sample exactly recreating the periodicity of the structure as well as the perturbations (if any) is created; 3. A stiffness matrix taking into account the influence of the perturbations is computed.

Computation of the stiffness matrix
In order to detail steps 2 and 3 of the sequence introduced above, a sample with two 8-node solid elements is used as an example as shown in Figure 1: the elements <1, 2, 3, 4, 101, 102, 103, 104> and <105, 106, 107, 108, 5, 6, 7, 8> are superimposed along the z-axis, and preloads are applied to the structure (red arrows). The two elements have to be separated if contact properties are to be taken into account (in such case, the interface nodes are doubled, although each node pair is perfectly coincident). If no contact or friction properties are modelled, the interface nodes are merged. The global cuboid's dimensions are L x , L y and L z , and its faces' respective areas A x (faces x = 0 and x = 1), A y (faces y = 0 and y = 1) and A z (faces z = 0 and z = 1). Also, the nodes located within an ellipse on the figure have the same coordinates. The description of external faces is simplified with '= 0' or '= 1' notations (referring to the limits of the structure's volume), in spite of the dimensions L x , L y and L z . For any types of finite-element simulations on a global structure (e.g. modal basis computation on a laminated core), a global stiffness matrix is computed by the solver for describing a given state (or increment) with the effects of perturbations such as preload or friction. Considering only a representative sample as illustrated in Figure 1, the corresponding stiffness matrix takes into account the same properties at the same state. The idea is therefore to approximate these properties with a linear homogeneous elasticity matrix describing the same state at which the stiffness matrix was output by the solver. Now, for the solution to be computed, the system has to be stabilised: the outer nodes 1, 2, 3, 4, 5, 6, 7 and 8 are respectively doubled with the nodes 11,12,13,14,15,16,17 and 18, at the exact same coordinates. The nodes 11 through 18 are then clamped, and each pair is linked with 3D stiffness elements (this is equivalent to linking a node-to-ground stiffness element to each of the structure's outer nodes 1 to 8). The numerical values of these stiffness elements have to be small enough to not perturb the entire structure's stiffness matrix values.
As an example, the stiffness matrix's diagonal components of a 1-mm-long cube made of isotropic steel with E = 210 GPa and ν = 0.25 have a value of 4.67 · 10 7 N · m −1 . In this case, using stiffness elements with a value of 1, 000 N · m −1 in every direction is considered negligible compared to the stiffness matrix's values, and does not perturb the structure's overall behaviour.
For the rest of this section, U x0 , U x1 , U y0 , U y1 , U z0 and U z1 are defined as the node sets corresponding to the respective faces x = 0, x = 1, y = 0, y = 1, z = 0 and z = 1. The case of a two-layer, 8-outer-node structure is detailed in Table 1.
The prestress effects σ 0 and σ 1 are then applied to the system along the z-axis, so that where the values F i stand for the loadings at the nodes 1 through 8. If modelled, the contact properties are then defined according to the physics of the interface (deformable bodies, friction, etc.). Creating the superelement with the twelve nodes 1, 2, 3, 4, 101, 102, 103, 104, 5, 6, 7 and 8 (DOFs T x , T y and T z ) is an efficient way to output a stiffness matrix 1 . It seems important to emphasise that the solver used has to take into account the influence of the perturbations over the values -if modelled -and export the corresponding stiffness matrix. In the case illustrated in Figure 1, the superelement created stands for a structure with two layers and four common interface nodes instead of two separate elements. The elasticity matrix computed with this new method is linear and approximates the non-linear perturbations to which the structure is subjected. In other terms, the stress-strain law it describes corresponds to the same situation in which the stiffness matrix was computed 2 .

Determination of the elastic properties
Unlike finite elements for which materials would be clearly defined, the stiffness matrix of a superelement does not give any direct information about the materials it stands for. The matrix may however be imported by a solver to recreate the elastic properties of an entire structure, in our case composed of two 8-node solid elements and twelve nodes (the node pairs at the interface were merged during the superelement's creation process). It is therefore possible to compute the displacements ∆l x,i , ∆l y,i and ∆l z,i as well as the reaction forces F x,i , F y,i and F z,i at each node i.
As it has been shown, the 21 independent coefficients of an equivalent triclinic elasticity matrix C have to be identified. To do this, the proposed approach consists in computing the compliance matrix S such that S = C −1 . The general relation of Hooke's law, detailed in Equation (1), can be reversed to express the matrix S : To calculate all the constants of S , the method needs six series of simulations, ie. one per component of {ε}. The first series, called xx, refers to the deformed state ε xx . For better visibility, this development is made on the second series, called yy, referring to the deformed state ε yy , and repeated similarly for the deformed states ε xx and ε zz .
The second row of Equation (4) is then: Thus, this series needs six independent simulations in order to identify Equation (5)'s six constantsS 12 =S 21 ,S 22 , S 23 ,S 24 ,S 25 andS 26 . This series yy corresponds to a pure tension scheme, as shown in Figure 2: a static displacement δ y is enforced along +y to the nodes of the face y = 1, whereas plane contact constraints are applied to the face y = 0.
In other terms, the nodes of the face y = 0 are blocked in the direction y and left free in the directions x and z. To equilibrate the system, DOFs T x and T z are blocked at node 1, and DOF T z at node 2 (more details are given in the appendices, in Table A.10).
and σ (2a) For other uses such as local stress field computation, it can be noted that these expressions may be adapted to specific node sets (e.g. for determining stress fields at interfaces). The next five simulations of the series yy (corresponding to superscripts (2b) through (2 f )) obviously recreate the same pure tension scheme along the y-axis. The difference made between each simulation for the same series yy is in the boundary conditions of the scheme (unlike in the work of Mathan and Siva Prasad [17]), i.e. which DOFs are blocked or set free, while enforcing the same displacements to all simulations (2a) through (2 f ). The independence between each simulation's boundary conditions must then ensure that the 6-equation system (one equation similar to (5) for each simulation of the series) is of rank 6. Finally, the equation system obtained is where the 6 × 6 matrix is referred to as H (2) , and where: Now, solving the system yields the six coefficientsS 21 =S 12 throughS 26 . These steps are repeated with two new simulation series (xx and zz), in order to complete the first three rows of the compliance matrix S , as where the sub-matrices S A and S B are both of size 3 × 3. The next step consists in using the symmetry of S to transpose the submatrix S B : so that computing the last three rows of S only requires three simulations each instead of six. The series yz refers to the deformed state γ yz = 2 · ε yz . The corresponding row in Equation (4) is then: In a similar way as shown by Millithaler et al. [18], the analysis separates sliding shear from transverse shear, in the respective "Tricl1" and "Tricl2" scenarios, for which the identification of the compliance matrix's first three rows yet remains identical. Taking as example "Tricl2" scenario (i.e. transverse shear), the identification of the remaining rows is detailed in the appendices. Eventually, the compliance matrix S is fully constituted from all six simulation series, and thus stands for homogeneous material properties representative of the initial heterogeneous structure.

Validation
For the above-presented method to be applied, three validation cases are detailed throughout the following paragraphs, each presenting a specific use of the identification method: • a simple homogeneous structure composed of a triclinic material; • a preloaded homogeneous structure and a comparison with the method developed by Millithaler et al. [18]; • a preloaded laminated structure for which the method developed by Millithaler et al. [18] could not be applied.
In addition to these validation cases, the next section will detail the creation of a representative finite-element model for an industrial stator core made of laminated steel and where preloading and friction effects are accounted for.

Homogeneous triclinic sample
The method "Tricl1" has been first validated on a simple, homogeneous case. For this example, a cuboid 8-node element whose dimensions are L x = 20 mm, L y = 40 mm and L z = 60 mm was used, to which a triclinic material defined by the elasticity matrix was applied, but without any external perturbations. Given these input material properties, this application will attempt to show to what extent the elasticity matrix identified by the method corresponds to the input [C exp ]. In this case, each simulation has been made with enforced displacements of magnitude δ = 1 mm.
A linear static solution is initiated (each solution is completed within a few seconds), including output requests at all nodes in terms of displacements and reaction forces. The results of the xx-series lead to the system (12) (adapted to row 1 i.e. direction xx), whose numerical values are given in the following system: and where ε xx = 0.05. The 6 × 6 matrix H (1) has a rank of 6. As it can be seen in Equation (19), the independence of the boundary condition sets led to a number of null or negligible values in matrix H (1) . This is done so that each simulation in the series involves a minimum of terms in the compliance matrix. In particular, the diagonal termS 11 is identified alone in the first row of H (1) . Solving Equation (19) yields the first row of the equivalent compliance matrix, as shown in Equation (20): The values found here are very close to the initial compliance matrix, including the terms of coupling between tension-compression and shear, namelyS 14 ,S 15 andS 16 . The maximum relative discrepancy U rel max is such that which is very low and therefore shows that the computed results are close to the initial values. The other rows are computed in the same way. The symmetry of S is confirmed, and the values match; the validation is successful for this triclinic homogeneous element. Identical results are found with method "Tricl2".

Preloaded homogeneous structure 3.2.1. Global structure:
The second application of the "Triclinic" method is an analysis of a preloaded structure. The objective is to analyse the effects of the preload on the initial material's elasticity matrix, and to what extent a linear, homogeneous material could recreate the vibratory behaviour of the initial, preloaded structure. The FE model is a homogeneous cuboid of isotropic steel (Young's modulus E = 207 GPa, Poisson's ratio ν = 0.292 and density ρ = 7875 kg · m −3 ) of respective dimensions along x, y and z of 100 mm, 70 mm and 80 mm, and has 1,008 elements and 8,034 DOFs. Tension preloads are applied along direction y to the structure. The faces y = 0 mm and y = 70 mm are subjected to static forces of respective total magnitudes 9.81 · 10 6 N and −9.81 · 10 6 N, equally distributed on the faces' nodes, so that ±81.1 · 10 3 N is applied to each of these nodes along y. This value has been chosen to be voluntarily high to ensure observing notable effects on the responses, yet none of the yield or fracture limits are taken into account in the simulation: the material is assumed to never reach any of these limits while calculating the solutions. Therefore, the total pre-stress field σ p yy applied to the structure is: Also, a node-to-ground stiffness element is linked to each of the global cuboid's 8 outer nodes, with stiffness values of K ntg = 1 N · m −1 on each direction x, y and z. The initial structure is illustrated in Figure 3.

Equivalent material properties:
To apply the identification method and determine equivalent material properties, a sample is created from a few elements of the structure: 3 elements along y (48 DOFs), as shown in Figure 4. To recreate the stress field existing in the global structure, the sample is subjected to the same pre-stress field σ p yy . To stabilise the system, a node-to-ground stiffness element is linked to each of the sample's 8 outer nodes, with stiffness values of K ntg = 1 N · m −1 on each direction x, y and z.
It seems important to note that in spite of the initial structure's boundary conditions, applying the presented methods to identify equivalent materials must be made in "free" conditions, or in other words without any DOF constraints. Yet, a structure subjected to preloads needs to be stabilised, which explains the addition of node-to-ground elements to the sample (for this example, although other solutions may exist), and which is completely independent node-to-ground stiffness elements from the global structure's boundary conditions. For this material identification, the stiffness values of the node-toground elements have to be sufficiently high to enable computing the stiffness matrix, and as low as possible to be negligible compared to the matrix's values. In this case, it must be verified that their values after preloading are still negligible compared to the new stiffness matrix's values.
A 48 × 48 real and symmetric stiffness matrix is computed, accounting for the influence of the preload. By creating a new model with the sample's 16 nodes (and no elements), and importing the stiffness matrix as an external superelement, a linear static solution is initiated to apply the method according to the scenarios "Tricl1" and "Tricl2" presented in Section 2. Post-processing the results yields the elasticity matrices C ISO1 and C ISO2 (respectively corresponding to "Tricl1" and "Tricl2" methods), detailed in the Appendices, in the respective Equations (A.17) and (A.18). As a reference for comparisons, the elasticity matrix C stl , corresponding to steel (without preloading), is shown in Equation (A.19).
In order to compare it with these results, the method "Ortho1" developed by Millithaler et al. [18] (intended to be applied to orthotropic laminated structures) has been applied to the same superelement (computed from the isotropic steel sample under preloading), and yielded the following elastic coefficients:Ẽ x = 201 GPa,Ẽ y = 242 GPa, E z = 201 GPa,G zy = 116 GPa,G zx = 49.5 GPa,G xy = 89.1 GPa,ν yz = 0.313,ν xz = 0.298 andν xy = 0.259. Judging from the values of the matrices, the following observations can be made: • While C stl has a shape inherent to isotropic properties, this is not the case of C ISO1 and C ISO2 : the preloading effects have altered the initial material's isotropy; • Relatively low terms of coupling between tension-compression and shear have been determined by the two scenarios "Tricl1" and "Tricl2"; • In both matrices C ISO1 and C ISO2 , the tension preloading along y resulted in an increase of the coefficient C 22 from its value in C stl , which is the diagonal term of Hooke's law in direction yy. This is consistent with the expected stiffening effect from tension preloading [4].
Two equivalent homogeneous structures are then computed, with the same dimensions and the same density as the initial model.

Modal correlation analysis:
To evaluate the capacity of the equivalent material properties to recreate the behaviour of the preloaded structure, a state of modal correlation is calculated. To perform this, modal bases of the first 200 vibration modes are computed for the structures detailed in Table 2. The first 6 modes, describing the "suspension" related to the node-to-ground elements (all below 14 Hz), have been discarded. Compared to this, the lowest 7 th mode frequency amongst all the structures (and thus first to be correlated) is at 14,5 Hz.
The criteria defining the correlation consist in comparing each vibration mode of the first model to each of the second. When comparing the modes of two structures, several criteria can be referred to [19]. The approach we   propose in this article is to compare both natural frequencies and deformed shapes, to assert the extent to which one structure reproduces the modal behaviour of the second. These criteria are the same as in the work of Millithaler et al. [18]. Therefore, the relative frequency difference ∆ f (m e,i , m a, j ) between the natural frequency of the first structure's i-th mode m 1,i and the second structure's j-th mode m 2, j is computed with the following relation: where f 1,i et f 2, j are the natural frequencies corresponding to the modes m 1,i and m 2, j , respectively. In addition to this, the similarities between the deformed shapes of the modes m 1,i and m 2, j (respectively called φ 1,i and φ 2, j ) are determined according to the so-called MAC criterion (Modal Assurance Criterion), by the expression [1]: Then, the pairs of modes for which MAC values are highest are assembled, and are taken into account if the MAC values are above a fixed threshold. All the other mode pairs are discarded from the correlation process. For the correlation, the reference modal basis is "Prld", to which the other bases are compared. The MAC-threshold is fixed at 0% for pairing the modes (so that all the modes are paired and taken into account). The results of the correlation are gathered in Table 3. For N pm mode pairs in a given correlation, the entities |∆ f | and MAC are defined by the expressions: and where m q 1 and m q 2 are the modes composing the q-th pair. The results given in Table 3 can be summarised as below: • Both natural frequencies and deformed shapes are affected by the application of preloads, as the figures of column "Init" show; • Scenario "Tricl2" is capable of recreating the behaviour of the structure under preloads with good accuracy, while method "Tricl1" is much less efficient in this setting.
• In spite of the similarities of the shapes of C ISO1 and C ISO2 with that of an orthotropic material's elasticity matrix (for which no tension-compression/shear and shear/shear coupling terms exist), it can be clearly seen that using the "Ortho1" method is not adapted to such a setting.

Global structure:
For the last validation case, the same structure as in the work of Millithaler et al. [18] has been analysed: a laminated cuboid of 5,024 elements, 30,144 DOFs and respective dimensions along x, y and z of 210 mm, 110 mm and 60 mm. The stack's base cell is composed of 3 isotropic layers, the properties of which are detailed in Table 4, and is oriented along z. For each layer, E is the Young's modulus, ν the Poisson's ratio, ρ the density and l the thickness. Also, the volume fraction χ n of layer n is defined by the layer's volume V n and the base cell's total volume V cell , so that and  In this case, the structure is subjected to tension preloads along the stacking direction (z). The total loads on the top and bottom faces are respectively 13.2 · 10 7 N and −13.2 · 10 7 N, equally distributed on the faces' nodes, so that ±500 · 10 3 N is applied along z to each of these nodes. As before, this value has been chosen to be voluntarily high to ensure observing notable effects on the responses, yet none of the yield or fracture limits are taken into account in the simulation: the material is assumed to never reach any of these limits while calculating the solutions. Also, no contact conditions are taken into account between the different layers: the structure is assumed to experience no delamination. Similarly to the corresponding structure in the work of Millithaler et al. [18], a node-to-ground stiffness element is linked to each of the global cuboid's 8 outer nodes, with stiffness values of K ntg = 10 7 N · m −1 on each direction x, y and z. The global structure taken as reference is illustrated in Figure 5.

Equivalent material properties:
To apply the material homogenisation approach and determine an equivalent material, a sample is created from a few elements of the structure. In this application, the sample consists of the 3-layered base cell (48 DOFs) of which the entire model is composed, and is illustrated in Figure 6. To recreate the stress field existing in the global structure, the sample's 8 outer nodes are subjected to the same nodal loads (±500·10 3 N per node), as the thicknesses of the base cell layers are identical in the sample and in the global structure. To stabilise the system, a node-to-ground stiffness element is linked to each of the sample's 8 outer nodes, with stiffness values of K ntg = 1 N · m −1 on every direction x, y and z. node-to-ground stiffness elements  As before, the addition of node-to-ground elements to the sample is necessary to stabilise the system, and yet independent from the initial, global structure. Again, the stiffness values of the node-to-ground elements are negligible in comparison to the sample's stiffness matrix's.
A 48 × 48 stiffness matrix is computed (which is real and symmetric) and takes into account the influence of the preload. By creating a new model with the sample's 16 nodes (and no elements), and importing the stiffness matrix as an external superelement, a linear static solution is initiated to apply the methods "Tricl1" and "Tricl2" presented in Section 2. Post-processing the results yields the elasticity matrices C LMT1 and C LMT2 (respectively corresponding to "Tricl1" and "Tricl2" methods), detailed in the Appendices, in the respective Equations (A.20) and (A.21). Associated to them, the matrix C O1 composed of the elastic moduli determined in [18]

(in the case of an orthotropic material without perturbations) is recalled in Equation (A.22).
Judging from the values of the matrices, several observations can be made: • Non-negligible terms of coupling between tension-compression and shear have been determined by the two methods "Tricl1" and "Tricl2"; • The matrices C LMT1 and C LMT2 are both positive definite (their eigenvalues are all strictly positive), which is a necessary condition for a system to be stable [5]; • In both matrices C LMT1 and C LMT2 , the tension preloading along z resulted in an increase of the coefficients C 33 ,C 44 andC 55 from their values in C O1 , which are the diagonal terms of Hooke's law in each direction involving z (respectively zz, yz and xz). This is consistent with the expected stiffening effect from tension preloading [4];

Tricl1
Homogeneous (scenario "Tricl2"transverse shear) Tricl2 Table 5: Definition of the modal bases • In both matrices C LMT1 and C LMT2 , the directions x and y have similar coefficients, which shows that the laminated structure has in this case an equivalent behaviour in planes normal to the stacking direction; • In C LMT2 , the absolute values involving shear in directions yz and xz are significantly greater than in C LMT1 , which is similar to the comparison between the "Ortho2" and "Ortho1" cases in the work of Millithaler et al. [18].
The homogeneous structure to which the equivalent material of each method is applied has the same dimensions and the same total mass as the reference cuboid, and is made of 1-centimetre-long cubic, 8-node, solid elements. Therefore, the equivalent densityρ is calculated by the relation of weighted average: where ρ n is the density of each layer n.
A density value and an elasticity matrix fully define a material for the computation of a real modal basis: an equivalent homogeneous structure is created with each of the material properties respectively computed with methods "Tricl1" and "Tricl2". For each of them, the dimensions are identical to the reference matrix's, but one-centimetre-long cubic elements replace the three-layered base-cells that were homogenised.

Modal correlation analysis:
To evaluate the capacity of the equivalent material to recreate the behaviour of the preloaded structure, the modal correlations between the models are analysed. To perform this, modal bases of the first 50 modes are computed for the structures detailed in Table 5. The first 6 modes, describing the "suspension" related to the node-to-ground elements (all below 2,500 Hz), have been discarded.
For the correlation, the reference modal basis is "Prld", to which the other bases are compared. The paired modes for which MAC values are below 70% are discarded. Also, comparing the dynamic behaviours of the homogeneous global structures no longer requires involving suspension stiffness elements, even though their values in the base cell were small enough to have negligible impacts on the identified elastic properties. For including them in the laminated global model under preloads, equivalent values can be computed with the approach described in the Appendices. The results of the correlation are gathered in Table 6, where |∆ f | corresponds to the average of the frequency differences' absolute values, and MAC the mean MAC value of the paired modes.  The results given in Table 6 can be summarised as below: • Both natural frequencies and deformed shapes have been affected by the application of preloads, as the figures of column "Init" show; • In spite of the important behaviour difference induced by the application of the preloads, method "Tricl1" is capable of identifying 36 of the 44 modes; • However, the material from method "Tricl2" is not efficient to simulate the behaviour of the initial structure under preloading, as only 17 of the 44 modes are identified.
As shown by Millithaler et al. [18], it can be said that recreating the lower-frequency modes of laminated structures with homogeneous equivalent material properties requires identifying them with sliding shear simulations instead of transverse shear. On the contrary, the results of the analysis in Paragraph 3.2 show that identifying equivalent material properties for a continuous structure is much more accurate with transverse shear simulations.

Finite-element model
The ability of an electric machine stator's finite-element model with orthotropic material properties to simulate the modal behaviour of the corresponding real structure has been shown by Millithaler et al. [18], along with modelling guidelines that led to zoning the model. The same stator is considered in this section, and consists in a stack of several hundreds of steel sheets separated from each other by varnish. During its manufacturing process, weld beads are applied on the lateral side of the stack, while the magnetic core is placed under a press. When the pressure is released, the stack is held in one piece by the weld beads, while in the rest of the structure, the only bond between the sheets is the varnish. This is a source of heterogeneities in the behaviour of the entire structure. Therefore, the finite-element model of the magnetic core has been divided into several zones according to the distance to the weld beads (see Figure 7), with specific material properties associated to each of the zones. The same model as in the work of Millithaler et al. [18] is used in this section. It is made of 19,158 elements and 24,768 nodes (expressed in a cylindrical coordinate system of directions r, θ and z) for 12 teeth, its dimensions are 154 mm (length) and 245 mm (outer diameter). The sheets are stacked along the z-axis.

Equivalent material properties
The present analysis focuses on the interaction between the steel sheets and the possibility to take such effects into account in the equivalent material properties. Because of their proximity to the weld beads, the zone "prox" is assumed to experience no friction effects between the layers. Therefore, this zone is associated to the same orthotropic material as in the work of Millithaler et al. [18].
On the contrary, friction effects between the steel sheets are modelled in the epoxy layers for the zones "yoke" and "teeth". Compressive preload is taken into account in order to model the residual pressure remaining after the  Table 7: Stacking sequence -base cell manufacturing process, and is assumed to be homogeneous in each zone. The unit (or base) cell periodically repeated in this zone is thus a set of four 8-node solid elements superimposed along the z-axis. The base cell is illustrated in Figure 8, and the stacking sequence is detailed in Table 7. The same densityρ (acquired from measurements) is applied to the entire structure:ρ = 7, 750 kg · m −3 . The entire base cell is composed of 24 nodes and thus 72 DOFs. The idea is to identify equivalent material properties from a superelement in which the master-nodes describe the glued elements (nodes 1 to 18 -written in bold in Table 7) and the slave-nodes the interface (nodes 201 to 208 -in italic). In this article, a superelement's "master" or "external" DOFs correspond to the DOFs kept after condensing (that still exist in the superelement's reduced stiffness and mass matrices). On the contrary, "slave" or "internal" DOFs are not present after the reduction, as explained for instance in the so-called "Craig-Bampton" method [8], widely used in FE simulations. The pairs of interface nodes have the same coordinates. To stabilise the system, a node-to-ground stiffness element is linked to each of the outer nodes 1, 2, 3, 4, 5, 6, 7 and 8, with stiffness values of K ntg = 1 N · m −1 on every direction r, θ and z.
To the authors' knowledge, there are currently no data available concerning either prestress values on a magnetic core's tooth or friction between two steel sheets. In this application, the initial load applied to the stator during its production is equivalent to a mass of 2,500 kg (i.e. ≈ 25 kN), and therefore corresponds to a compressive prestress of 1.52 · 10 6 N · m −2 along the z-axis.
In order to evaluate the distribution of residual stresses in the different zones of the stator, the static stiffness values computed in the work of Millithaler et al. [18] are taken as references. The mean values were 3.0 · 10 8 N · m −1 in zone "prox", 2.5 · 10 8 N · m −1 in zone "yoke" (i.e. 83% of zone "prox") and 1.3 · 10 8 N · m −1 in zone "teeth" (i.e. 43% of zone "prox"). Using the same factors to represent the distribution of prestresses, a compression of 1.27 · 10 6 N · m −2 is applied on zone "yoke" and 6.58 · 10 6 N · m −2 on zone "teeth". Apart from this prestress, the identification of equivalent properties is made in the same way and with the same base cell for both zones. The contact properties are described by a Coulomb dry friction with a coefficient µ = 0.9, and occur at the interface described previously.
A 72 × 72 real symmetric stiffness matrix is computed and takes into account the influence of the preload and friction. Reducing it with the "Craig-Bampton" [8] method (widely used in FE simulations) with master-nodes 1 to 18 yields a new stiffness matrix (real, symmetric and of size 48 × 48). By creating a new model with the 16 masternodes (and no elements) and importing the reduced stiffness matrix as an external superelement, a linear static solution is initiated to apply the method "Tricl1" presented in Section 2 (sliding shear scenario). Post-processing the results yields the elasticity matrices C yoke and C teeth , detailed in the Appendices, in the respective Equations (A.23) and (A.24).
Judging from the values of matrices C yoke and C teeth , several observations can be made: • non-negligible inter-shear coupling terms have been determined for both zones, as well as non-negligible tension-shear coupling terms for zone "yoke"; • in rows 4 and 5 of both matrices, the diagonal terms correspond to the equivalent sliding shear stiffness values in directions θ − z and r − z; • in both matrices, the rows relative to series rr and θ θ (as well as θ z and rz) have similar values. Globally, this expresses the fact that directions r and θ are equivalent, which is consistent with the apparent symmetry of the base cell; • compared to the material properties of homogeneous isotropic steel (see Equation (A.19)), the diagonal values of both matrices C yoke and C teeth are lower than in matrix C stl , especially for rows 4 and 5. This seems consistent with the fact the base cells are subjected to friction instead of being glued elements (with common nodes); • compared to each other, the diagonal values of rows 4, 5 and 6 are higher in zone "yoke" than in zone "teeth". This indicates that the tightening effect due to prestress is more notable in the yoke's shear properties than in the teeth's.

Comparison with experimental data
A modal basis is simulated in real domain between 0 and 10,000 Hz from the entire magnetic core's finite-element model. This simulated modal basis is compared with a set of 5 radial modes (same experimental modal basis as in the work of Millithaler et al. [18]), extracted from frequency response functions measured with an impact hammer on the magnetic core of a real stator (purely radial impact). The frequency response function measured on the radial direction at the impact point is shown in Fig. A.14 (placed in the appendices). These modes are sometimes referred to as "cylinder" or "ovalisation" modes, and can be clearly spotted by the peak magnitude values on the response function shown in Fig. A.14. Although their shapes are purely radial in theory, the ovalisation is not pure in practice due to boundary effects. This motivates accounting for friction and preload in the modelling process of the structure.
The resonance probability of ovalisation modes is yet strong due to the nature of the electromagnetic excitations occurring at the air gap (i.e. the inner tooth surface) while the machine is operating [10]. Therefore, they are critical for the acoustic behaviour of the entire stator: being able to predict them accurately is then of particular interest 3 .
The experimental settings of the measurements are presented on Fig. 9 and Fig. 10. The experimental mesh is composed of 108 degrees of freedom (36 points of 3 DOFs). The comparison of the simulated and experimental modal bases is presented in Fig. 11 and Table 8, where the columns "FEA" and "EMA" respectively refer to the mode frequencies in the FE model and in the experimental modal basis. The |∆ f | and MAC averages have been computed with Eq. (25) and (26) and shown in the bottom line of Table 8. In addition to this, the similarities between simulated and measured mode shapes are illustrated with the MAC-matrix shown in the appendices, in Fig. A.15.
In order to estimate the accuracy improvement of taking into account friction and preload effects on the structure, the results of Table 8 (called "Ortho1 + Tricl1") are compared with "Ortho1" averages (neglecting friction and preload effects). Furthermore, an additional model is built with equivalent properties computed with the so-called "rule of mixtures", based on weighted averages. This last model will be referred to as "WA", does not involve friction or preload effects, and consists in a single equivalent material property set for the entire structure, whose values are:Ẽ r =Ẽ θ = 205 GPa,Ẽ z = 132 GPa,G zθ =G zr = 51.2 GPa,G rθ = 82.1 GPa,ν θ z =ν rz = 0.16 andν rθ = 0.25. Definitions of both "Ortho1" and "WA" identification methods are detailed in the work of Millithaler et al. [18]. Finally, the comparison of the average values |∆ f | and MAC for the above-mentioned methods is shown in Table 9.
In Table 9, the 17% relative decrease of |∆ f | and only 4% relative decrease of MAC between "Ortho1 + Tricl1" and "Ortho1" scenarios show a significant improvement in the simulation accuracy of the above-presented ovalisation modes, of particular importance regarding electric motor acoustics [10]. This tends to compensate the complexity increase due to replacing explicit orthotropic elastic constants by a fully-defined elasticity matrix for the zone "teeth", although the identification method can be automatised with low-resource computer algorithms, and the computation costs for any FE simulations remain unchanged regardless of the type of materials. As for the "WA" scenario, it can be seen that modelling the stator core with a single equivalent material identified with the "rule of mixtures" leads to an overall frequency discrepancy which is more than twice as large as for "Ortho1 + Tricl1".

Specificities of the methods
Several details about the properties of the method presented above have to be added. For both the "orthotropic" and "triclinic" methods, creating a sample whose global geometry is not a cuboid highly complicates the establishment of the solution. Indeed, if one of the element's faces is not a perfect rectangle (for instance if one of its angles differs from 90°), the reaction forces computed on neighbouring faces will be misused. This tendency is particularly important for the identification of orthotropic and isotropic materials in unknown structures. For similar reasons, using the "triclinic" method on homogeneous, orthotropic structures may lead to fewer variables than equations while computing the matrix H (see Equation (12)), and therefore to a rank lower than 6. In such cases, pseudo-inverting H may retrieve the expected values, including null ones.   The elements' geometric properties have another influence on the method. In the case of heterogeneous structures, for instance a laminated composite made of isotropic layers, the elastic constants E, G and ν vary discontinuously in the stacking direction, thus implying heterogeneities and discontinuities in the computed reaction forces. Indeed, while applying the "triclinic" method to a stack made of isotropic layers without perturbations, the simulations of the series xx generate important and uneven reaction forces at the nodes constrained with plane contacts on faces y = 0 and z = 0, because of the layers' different materials. This comes out as non-negligible coupling terms between tension-compression and shear, and therefore contradicts the assumption of orthotropy due to the system's symmetries [5,2]. This is why a laminated structure presenting a priori no coupling between tension-compression and shear needs to be analysed with the "orthotropic" method presented in the work of Millithaler et al. [18]. In practice, applying a single pure tension simulation to a cuboid structure would confirm the existence of non-null coupling terms between tension-compression and shear in the equivalent material, and would therefore suggest using either the "orthotropic" or the "triclinic" methods.
Eventually, the last point concerns the ability of both the "orthotropic" and "triclinic" methods to be applied to superelements. In a structure to be reduced into a superelement and taken as sample for material identification, it must be clear that none of the slave DOFs may be located on external surfaces, or correspond to any of the simulations' constrained DOFs (as listed in Table A.10). Otherwise, the boundary conditions have to be taken into account before reducing the structure, which therefore requires as many independent superelements as independent simulation  Table 9: Accuracy improvement schemes. All other nodes, not located on external faces, may be taken as slave DOFs and reduced, without hindering the identification process. The "Triclinic" method can be summarised in the following sequence: for a continuous anisotropic model, a continuous structure subjected to preloading and a preloaded laminated stack have been exhibited. As a result, it has been shown that anisotropic materials can be identified with little derivation on continuous structures if no perturbations are applied. The other analyses have shown that preloading effects alter the initial symmetries in the material properties on a 3D isotropic model, and induce couplings between tensioncompression and shear in the equivalent material properties of a multi-layered laminated stack. It has been observed that under preloads, anisotropic continuous structures require transverse shear simulations, whereas sliding shear identification scenarios are more accurate for recreating the behaviour of laminated models. Also, the modal behaviour of an electric motor stator's laminated magnetic core has been simulated with equivalent anisotropic material properties that accounted for the friction behaviour under compression preloads between steel sheets in the teeth and the yoke. Low-frequency ovalisation modes have been computed and showed good accuracy in comparison to experimental data from a real stator. In comparison to simpler modelling approaches, it has also been shown that dividing the model into several zones and taking into account the effects of preloading and friction in the case of the stator core led to accuracy improvements for the simulation of ovalisation modes. An accuracy improvement has been also observed in comparison to orthotropic properties. This new identification method raises hopes to improve the current prediction capacities to perform noise and vibration simulations on multi-layered magnetic cores, without needing to rely on experimental data from costly prototypes and time-consuming model updating procedures. The ability of this new method to be applied to superelements, and therefore estimate the influence of perturbations on the material properties, presents a "conversion" opportunity from stiffness matrices to elasticity matrices.

AppendixA. Appendices
The following subsections gather some data and applications that were not shown in this article's main matter.

AppendixA.1. Determination of equivalent material properties for transverse shear scenarios
Following the identification steps presented in the Subsection "Determination of the elastic properties", the following paragraphs detail the determination of the remaining compliance coefficients by transverse shear schemes (therefore corresponding to the "Tricl2" method). The fourth series of simulations corresponds to a transverse shear scheme y − z, as shown in Figure A.12: it combines an enforced displacement δ z along +z at the nodes of the face y = 1, and plane contact constraints on the same face y = 1 in order to generate pure shear. In this example also, DOFs T x , T y and T z are blocked at node 1, and DOFs T y and T z at the other nodes of face y = 0. As is the case for the previous series, Equations (6) through (11) yield the stress values needed for this simulation. Also, new boundary conditions enable creating two new independent simulations to complete this series (more details are given in Table A.10). The system obtained is thus: The 3 × 6 matrix H (4) can be divided into the two submatrices σ (4TC) and σ (4SR) , respectively consisting of the terms of tension-compression and shear: In detail, they stand for with the help of which the remaining unknown coefficientsS 44 ,S 45 andS 46 may be computed: At last, the matrix S is completed with the last two series, namely xz and xy, from which the symmetryS i j =S ji has to be verified again. The detail of all boundary condition sets is given in Table A. 10.

AppendixA.2. Additional details for validation
Following the validation case described in Subsection "Preloaded laminated structure", the equivalent values of the node-to-ground elements after preloading can be computed from the stiffness matrix of the global structure under preloads. The idea is to select a corner node, to which a node-to-ground stiffness element is linked, and to extract the values involving both its neighbouring nodes (except for the node at the interface of two different materials) and itself. As an example, the corner node 1 and its direct neighbours 2 and 4 are selected (see Figure  When assembling the stiffness matrices of different elements in a system, the values corresponding to identical DOFs present on several elements are added [7]. As nodes 2 and 4 are each linked to two elements having identical stiffness matrices (same dimensions and material), the diagonal values in the global stiffness matrix K f of the free system (without either node-to-ground stiffness elements or boundary conditions) verify the equations: In the case of a preloaded structure, a way of computing equivalent stiffness values K LMT x , K LMT y and K LMT z for the node-to-ground elements is by averages. Therefore, from the preloaded system's global stiffness matrix K p , we have: where K p is the stiffness matrix of the global system under preloads. For the system studied, the final values are: K LMT x = 20.4 · 10 6 N · m −1 , K LMT y = 19.4 · 10 6 N · m −1 and K LMT z = 120 · 10 6 N · m −1 . This shows that the influence of the preloads affects the stiffness elements in a similar way along the directions x and y, and that the equivalent suspension stiffness values are greater than the original ones (K ntg = 10.0 · 10 6 N · m −1 ).

AppendixA.4. Frequency response function measured on the magnetic core
A frequency response function, obtained in the measurements described in Section "Electric machine stators: finite-element modelling accounting for frictional effects", is shown in Fig. A.14. The magnitude of the signal is an acceleration per unit force. This function has been obtained with a radial excitation by impact hammer and shows the response on the radial direction and the impact point.

AppendixA.5. MAC coefficients between simulated and measured mode shapes
The correlation presented in Section "Electric machine stators: finite-element modelling accounting for frictional effects" involved both simulated and measured modal bases. The similarities between mode shapes have been expressed with the aid of so-called MAC-coefficients, defined with Eq. (24). The corresponding matrix comparing simulated and measured modes up to 6,550 Hz is illustrated in Fig. A.15, where only coefficients above 40% are considered for better readability.