# Structural dynamics of electric machine stators: Modelling guidelines and identification of three-dimensional equivalent material properties for multi-layered orthotropic laminates

Pierre Millithaler, Emeline Sadoulet-Reboul, Morvan Ouisse, Jean-Baptiste Dupont, Noureddine Bouhaddi

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Pierre Millithaler, Emeline Sadoulet-Reboul, Morvan Ouisse, Jean-Baptiste Dupont, Noureddine Bouhaddi. Structural dynamics of electric machine stators: Modelling guidelines and identification of three-dimensional equivalent material properties for multi-layered orthotropic laminates. Journal of Sound and Vibration, 2015, 348, pp.185 - 205. 10.1016/j.jsv.2015.03.010. hal-01512519

HAL Id: hal-01512519

https://hal.science/hal-01512519

Submitted on 23 Apr 2017

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Structural dynamics of electric machine stators: modelling guidelines and identification of three-dimensional equivalent material properties for multi-layered orthotropic laminates

P. Millithaler<sup>a,b,*</sup>, É. Sadoulet-Reboul<sup>b</sup>, M. Ouisse<sup>b</sup>, J.-B. Dupont<sup>a</sup>, N. Bouhaddi<sup>b</sup>

<sup>a</sup>Vibratec, 28 chemin du petit bois, BP36, F-69131 Écully Cedex, FRANCE

<sup>b</sup>Femto-S.T. Applied Mechanics, UMR CNRS 6174, 24 chemin de l'épitaphe, 25000 Besançon, FRANCE

## Abstract

Simulating the dynamic behaviour of heterogeneous finite-element structures such as electric motors often requires to homogenise the models in the first place. Current homogenisation methods do not always imply computing an equivalent homogeneous material's elasticity matrix and are often restrained to specific uses. In this document, a novel approach of equivalent material identification is developed for multi-layered orthotropic structures. A finite-element model of a 3D stratified structure is created, as well as its equivalent homogeneous medium. The dynamic behaviour of the homogeneous structure with the equivalent material identified by the new method is compared at low frequencies to the reference stack and to equivalent materials created using other existing homogenisation techniques. It is shown that this approach is more accurate than existing reference homogenisation methods. Applied to the magnetic core's finite-element model of a real laminated electric machine stator, the method enables simulating the experimental behaviour with good accuracy, without need of time-consuming model updating procedures.

**Keywords:** homogenisation, laminated structures, effective medium, electric machine stator, magnetic core, superelements

## 1. Introduction

In order to analyse a complex structure's dynamic behaviour, modelling its components may become difficult if they are numerous, small, or if some of the assembly properties are not known. This is one of the reasons why so-called "homogenisation" methods have been developed. They aim at recreating a given heterogeneous structure's behaviour by reducing the multiplicity of its components' properties, and also enable to mesh the structure independently from the sizes of the heterogeneities, that could have imposed numerous degrees of freedom in the models. Such methods are much desired for modelling composite materials, and especially for laminated structures. Kalamarov et al. [1] have thoroughly listed and compared over 200 studies about homogenisation and have assessed the pros and cons of various analytical methods and specific applications.

For industrial projects in structural dynamics as well as academic research involving finite-element simulations on heterogeneous structures, the need of efficient techniques able to model homogeneous

<footnote>*Corresponding author. Tel: +33 4 72 86 65 65</footnote>

<footnote>Email addresses: pierre.millithaler@vibratec.fr (P. Millithaler), emeline.sadoulet-reboul@univ-fcomte.fr (É. Sadoulet-Reboul), morvan.ouisse@univ-fcomte.fr (M. Ouisse), jean-baptiste.dupont@vibratec.fr (J.-B. Dupont), noureddine.bouhaddi@univ-fcomte.fr (N. Bouhaddi)</footnote>

<footnote>The nomenclature is detailed in Table 1, on page 2.</footnote>
<table><caption>Table 1: Nomenclature</caption><thead><tr><th>Symbols</th><th>Definitions</th></tr></thead><tbody><tr><td>F<sub>d,i</sub></td><td>reaction force on node i, direction d</td></tr><tr><td>T<sub>x</sub>, T<sub>y</sub>, T<sub>z</sub></td><td>translational DOFs</td></tr><tr><td>Δl<sub>d,i</sub></td><td>displacement of node i, direction d</td></tr><tr><td>Ψ</td><td>global/equivalent property Ψ</td></tr><tr><td>subscript n</td><td>referring to layer n</td></tr><tr><td>subscripts x, y, z</td><td>referring to directions x, y, z</td></tr><tr><td>DOF</td><td>degree of freedom</td></tr><tr><td>FE</td><td>finite-element</td></tr><tr><td>2D, 3D</td><td>two-dimension, three-dimension</td></tr><tr><td>B<sub>uN</sub></td><td>number of nodes on face u = N</td></tr><tr><td>U<sub>hN</sub></td><td>set of nodes on face h = N</td></tr></tbody></table>

equivalent structures both accurately and cost-effectively is great [2]. Therefore, there is a significant interest for simple tools yielding results directly usable for common finite-element software. Motivated by the current coming-up of hybrid or 100%-electric vehicles, the development of silent devices (as well as other noise, vibration and harshness (NVH) specifications) in the automotive industry involves finite-element modelling of heterogeneous structures such as electric motor stators [3]. In this perspective, the main objective of this article is to establish a method to determine representative equivalent material properties (elasticity matrices) for laminated structures, especially for the applicability to any kind of finite-element simulation (including dynamic responses, model updating, etc.) without being limited to specific cases.

Concerning stratified materials, many applications are made under the assumption of plane stresses and strains (for instance with laminated beams or shells), for which theory predicts static and dynamic behaviours with good accuracy (e.g. [4]). In addition, there also exist exact theories and solutions describing the vibration of stratified beams and plates, such as the works [5] and [6], as well as ready-to-use 2D laminated finite elements present in several software programmes (see e.g. [7]). A finite-element-based homogenisation technique taking into account viscoelastic properties in 2D-laminates has been proposed by Koishi et al. [8]. However, some other analyses can not be simplified by such assumptions – take the case of no dimension being negligible in the model – and have to be meshed in 3D.

A short review of some existing “3D-homogenisation” methods is made, as well as the fields of their applications. Such techniques are particularly relevant when e.g. parts of 3D finite-element models are multi-layered and need to be homogenised. First of all, the relations that may be the simplest for determining a homogeneous material equivalent to a heterogeneous structure, and that are used in many studies (including reference works in the field of composite materials, such as [4] and [6]), are weighted averages of the different components’ elastic constants, sometimes referred to as the “rule of mixtures”. As shown in [6], the expressions are built from the decompositions of the structure’s global stresses and strains according to each of the components.

Let us consider a laminated structure composed of N isotropic layers, stacked along the z-axis. For each layer n, the material’s corresponding density $ρ_n$, Young’s modulus $E_n$ and Poisson’s ratio $v_n$ are initially known, as well as its volume fraction:
in which 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. In its most general definition, matrix C stands is a fourth-order tensor [11], such as:

$$ \mathbf{C} = \begin{bmatrix} C_{1111} & C_{1122} & C_{1133} & C_{1123} & C_{1113} & C_{1112} \\ C_{2211} & C_{2222} & C_{2233} & C_{2223} & C_{2213} & C_{2212} \\ C_{3311} & C_{3322} & C_{3333} & C_{3323} & C_{3313} & C_{3312} \\ C_{2311} & C_{2322} & C_{2333} & C_{2323} & C_{2313} & C_{2312} \\ C_{1311} & C_{1322} & C_{1333} & C_{1323} & C_{1313} & C_{1312} \\ C_{1211} & C_{1222} & C_{1233} & C_{1223} & C_{1213} & C_{1212} \end{bmatrix}. \quad (10) $$

Begis et al. [10] propose a general definition of a Y-periodic structure's equivalent elasticity matrix from an asymptotic homogenisation approach. Considering the unit cell periodically repeated in the entire structure and represented in the set of directions $y = \{y_1, y_2, y_3\}$, Y-periodic vectors $W^{pq}(y)$ need to be calculated, as solutions of

$$ \frac{\partial}{\partial y_j} (C_{ijkh}(y) \cdot \varepsilon_{kh}(W^{pq})) = -\frac{\partial}{\partial y_j} C_{ijpq}(y) \quad (11) $$

in Y. Then, the homogenised coefficients $\tilde{C}_{ijkh}$ are computed by:

$$ \tilde{C}_{ijkh} = \frac{1}{\text{vol}_Y} \int_Y [C_{ijkh}(y) - C_{ijpq}(y) \cdot \varepsilon_{kh}(W^{pq})] dy, \quad (12) $$

where vol$_Y$ is the actual volume of unit (base) cell. This approach was compared in [10] to the works [12] and [13] which addressed the same issues, and has been declared more accurate by the authors. Also, other studies were based on the same general expressions (see Eq. (11) and Eq. (12)) and have proposed similar analytical approaches – amongst others, the articles of Lukkassen et al. [14], Kalamkarov et al. [15], Hassan et al. [16] and Cheng et al. [17]. In spite of the appropriateness to describe periodic structures with equivalent elasticity matrices in the general case, this method may be complicated to apply due to the need of resolving condition (11) in the first place. This difficulty is outlined in [18], where the resolution techniques of the general expressions (11) and (12) are denoted to be “nebulous” by the authors.

However, this general homogenisation technique has been applied to laminated structures in [9], and led to “clear” expressions given for the same case as before: a structure of N isotropic layers stacked along the z-axis. The algorithm is referred to as “INRIA” (“Institut National de Recherche en Informatique et en Automatique”, i.e. french Institute for Research in Computer Science and Automation) in the entire article, and is detailed as following. Each layer's elasticity matrix $C_n$ is defined with Lamé's coefficients $\lambda_n$ and $\mu_n$ so that

$$ \mathbf{C}_n = \begin{bmatrix} \lambda_n + 2\mu_n & \lambda_n & \lambda_n & 0 & 0 & 0 \\ \lambda_n & \lambda_n + 2\mu_n & \lambda_n & 0 & 0 & 0 \\ \lambda_n & \lambda_n & \lambda_n + 2\mu_n & 0 & 0 & 0 \\ 0 & 0 & 0 & \mu_n & 0 & 0 \\ 0 & 0 & 0 & 0 & \mu_n & 0 \\ 0 & 0 & 0 & 0 & 0 & \mu_n \end{bmatrix}. \quad (13) $$

With the aid of the coefficients $I_0$, $I_1$, $I_2$ and $J_0$ such that
$$I_h = \sum_{n=1}^{N} \frac{(\lambda_n)^h}{\lambda_n + 2\mu_n} \cdot \chi_n \quad (14)$$

and

$$J_0 = \sum_{n=1}^{N} \frac{\chi_n}{\mu_n}, \quad (15)$$

the constants $\tilde{C}_{11}$ are computed with the global relations [9]:

$$\tilde{C}_{11} = \frac{I_1^2 - I_0 I_2}{I_0} + \sum_{n=1}^{N} (\lambda_n + 2\mu_n) \cdot \chi_n, \quad (16)$$

$$\tilde{C}_{22} = \tilde{C}_{11}, \quad (17)$$

$$\tilde{C}_{12} = \frac{I_1^2 - I_0 I_2}{I_0} + \sum_{n=1}^{N} \lambda_n \chi_n, \quad (18)$$

$$\tilde{C}_{13} = \tilde{C}_{23} = \frac{I_1}{I_0}, \quad (19)$$

$$\tilde{C}_{33} = \frac{1}{I_0}, \quad (20)$$

$$\tilde{C}_{44} = \tilde{C}_{55} = \frac{1}{J_0} \quad (21)$$

and

$$\tilde{C}_{66} = \sum_{n=1}^{N} \mu_n \chi_n. \quad (22)$$

The expressions, initially defined with an integration over the entire structure's thickness, have been here expressed under a discrete form. The laminated stack's equivalent elasticity matrix is therefore composed:

$$\tilde{\mathbf{C}} = \begin{bmatrix} \tilde{C}_{11} & \tilde{C}_{12} & \tilde{C}_{13} & 0 & 0 & 0 \\ \tilde{C}_{11} & \tilde{C}_{13} & 0 & 0 & 0 & 0 \\ \tilde{C}_{33} & 0 & 0 & 0 & 0 & 0 \\ & \tilde{C}_{44} & 0 & 0 & 0 & 0 \\ \text{sym.} & & \tilde{C}_{44} & 0 & 0 & 0 \\ & & & \tilde{C}_{66} & 0 & 0 \end{bmatrix}, \quad (23)$$

and by the way verifies the property of transverse isotropy

$$\tilde{C}_{11} = \tilde{C}_{12} + 2 \cdot \tilde{C}_{66}. \quad (24)$$

Identical expressions have been proposed in [19]. As for Cecchi and Sab [20], they have proposed a numerical homogenisation method for orthotropic structures based on a Reissner-Mindlin model, and
applied it to brick walls. Although efficient, the algorithm nonetheless requires computing correction factors amidst the numerous computation stages and therefore seems tricky to implement.

There exist other numerical models. The works of Watt et al. [21], Kamiński and Kleiber [22], and Magalhães Dourado [23] have aimed at identifying composite structures' elastic coefficients by the means of statistical methods, that seem complicated to use in simple, deterministic studies. As for Araújo et al. [24], they have proposed a numerical model for the identification of elastic properties of laminates, applied it to the dynamic behaviour of a stratified plate and compared the results with experimental measurements. Yet, this approach does not identify all 6 lines of Hooke's law elasticity matrix, and is therefore too restrictive to be applied to generic 3D structures.

Apart from homogenisation techniques, the recent development of several types of industrial materials and devices have motivated research in the comprehension and prediction of laminated structures' dynamic behaviours. A type of device currently under the spotlight is alternating electric machines, the stators of which are built from a multi-layered laminated structure, called the magnetic core. Amongst the studies made in this field, the works of Wang [25], Verma [26], Williams [27] or Le Besnerais [28] bring sensible notions of stator and magnetic core dynamics to the table, yet without exhibiting elastic properties or modelling techniques. As for modelling guidelines, the necessity of taking other heterogeneities such as weld beads into account in a homogenisation process of the magnetic core has, to the authors' knowledge, not yet been addressed in the literature.

As for the use of composite materials in commercial software, the implementation of 3D laminated structures is not documented or is restricted to specific non-linear or static studies. Only few works tackling homogenisation issues show the comparison of their results to those obtained with commercial software in such cases. Based on Barker's results [29] – stating that 3D-homogenisation of composite structures necessarily induces detrimental errors – Kuhlmann and Rolfes have developed their own 3D-laminated finite element [30], seemingly as accurate as MSC.MARC™'s. However, to the authors' knowledge this finite element has not yet been implemented into any software packages available on the market.

Finally, concerning experimental analyses, several approaches exist to measure the elastic behaviour of a structure. Hearmon [31] and Hayes [32] have detailed a few ways to identify a structure's entire elasticity matrix from measurements. But as this method requires several types of analysis and several types of samples, it seems very difficult to inspire an analogous application to finite-element models. This same difficulty also compromises the use of other methods, such as those developed in the articles of Pierron et al. [33] and Rikards et al. [34, 35, 36].

In this section, several types of approaches have been discussed for the application to laminated structures. Some of them are numerical or analytical/asymptotic homogenisation techniques, or deal with experimental applications, while others are based on finite-element models. The novel method proposed in this article belongs to this last category: it is based on finite elements. The emphasis is made on representative elasticity matrices, which constitute stress-strain laws for the considered materials and determine stiffness matrices for both static and dynamic analyses. The following sections will detail its applications to simulate the dynamic behaviours of multi-layered orthotropic laminates as well as the magnetic core of an electric motor stator. The last section will also investigate efficient modelling techniques for FE models of electric motor stators and compare the modal simulations to experimental data from an industrial structure.
## 2. Development of the “Orthotropic” method

In this section, a new approach aiming at determining the elastic properties of a heterogeneous structure is proposed. The equivalent material thus defined is assumed to be orthotropic (as explained in Section 1) and is characterised by nine elastic coefficients: $E_x, E_y, E_z, G_{xy}, G_{yz}, G_{zx}, G_{xy}, v_{yz}, v_{xz}$ and $v_{xy}$. By the means of a limited number of static simulations and by simple processing of the corresponding displacements in the structure, all nine elastic coefficients can be computed. The definition of shear in the case of laminates can be ambiguous due to the orientation of the structure, and is therefore addressed in this article.

The study taken as an example in this section is a stack along the z-axis of three isotropic thick layers, for which the theory predicts a global transversely isotropic behaviour [11]. The structure is a cuboid $<1, 2, 3, 4, 5, 6, 7, 8>$ composed of three 8-node solid elements, as shown in Fig. 1. The nodes 101, 102, 103, 104, 105, 106, 107 and 108 represent the respective interfaces with their two neighbouring elements. Also, the coordinate system for the whole structure is global (i.e. rectangular unitary system) and its origin taken as node 1. The 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$). The descriptions of external faces is made with '$= 0$' or '$= 1$' notations (referring to the limits of the structure's volume), despite $L_x, L_y$ and $L_z$.

<figure><img src="image_3.png" /><figcaption>Figure 1: Structure geometry</figcaption></figure>

The first simulation is a pure tension along the x-axis shown on Fig. 2: a static displacement $\delta_x$ is enforced along $-x$ to the nodes 1, 4, 5, 8, 101, 104, 105 and 108 (face $x = 0$), while the same displacement $\delta_x$ is enforced along $+x$ to the nodes 2, 3, 6, 7, 102, 103, 106 and 107 (face $x = 1$). In order to stabilise the system, plane contact constraints are applied to the faces $y = 0$ and $z = 0$.

By computing the nodal displacements $\Delta l_{x,i}, \Delta l_{y,i}$ and $\Delta l_{z,i}$ as well as the reaction forces $F_{x,i}, F_{y,i}$ and $F_{z,i}$ at every node $i$, the stress $\sigma_{xx}$ can be found using the relation

$$ \sigma_{xx} = \frac{1}{A_x} \sum_{i \in \mathcal{U}_{x1}} \frac{F_{x,i}}{2}, \qquad (25) $$
<figure><img src="image_4.png" /><figcaption>Figure 2: Pure tension along x</figcaption></figure>

where $\mathcal{U}_{x1}$ is the set of nodes located on face $x = 1$. Therefore, the Young's modulus $\tilde{E}_x$ in the x-direction may be found with [37]

$$ \tilde{E}_x = \sigma_{xx} \cdot \frac{L_x}{\delta_x}, \qquad (26) $$

where the ratio $L_x/\delta_x$ is related to the normal strain $\epsilon_{xx}$ by

$$ \frac{L_x}{\delta_x} = \frac{1}{\epsilon_{xx}}. \qquad (27) $$

The Poisson's ratios $\tilde{v}_{xy}$ and $\tilde{v}_{xz}$ can then be computed with the mean displacements $\overline{\Delta l_y}(\mathcal{U}_{y1})$ and $\overline{\Delta l_z}(\mathcal{U}_{z1})$, such as

$$ \overline{\Delta l_y}(\mathcal{U}_{y1}) = \frac{1}{\mathcal{B}_{y1}} \sum_{i \in \mathcal{U}_{y1}} \Delta l_{y,i} \qquad (28) $$

and

$$ \overline{\Delta l_z}(\mathcal{U}_{z1}) = \frac{1}{\mathcal{B}_{z1}} \sum_{i \in \mathcal{U}_{z1}} \Delta l_{z,i}, \qquad (29) $$

where $\mathcal{B}_{y1}$ and $\mathcal{B}_{z1}$ respectively refer to the numbers of nodes on the faces $y = 1$ and $z = 1$ (e.g. for a stack with three layers, $\mathcal{B}_{x1} = \mathcal{B}_{y1} = 8$ and $\mathcal{B}_{z1} = 4$). This finally yields:

$$ \tilde{v}_{xy} = \frac{-\overline{\Delta l_y}(\mathcal{U}_{y1})}{2 \cdot \delta_x} \cdot \frac{L_x}{L_y} \qquad (30) $$
and

$$ \tilde{\mathbf{V}}_{xz} = \frac{-\overline{\Delta l_z^{(\mathcal{U}_{z1})}}}{2 \cdot \delta_x} \cdot \frac{L_x}{L_z}. \qquad (31) $$

Two similar simulations along y and z yield the global structure's remaining Young's moduli and Poisson's ratio:

$$ \tilde{E}_y = \frac{L_y}{A_y \cdot \delta_y} \sum_{i \in \mathcal{U}_{y1}} \frac{F_{y,i}}{2}, \qquad (32) $$

$$ \tilde{E}_z = \frac{L_z}{A_z \cdot \delta_z} \sum_{i \in \mathcal{U}_{z1}} \frac{F_{z,i}}{2} \qquad (33) $$

and

$$ \tilde{\mathbf{V}}_{yz} = \frac{-\overline{\Delta l_z^{(\mathcal{U}_{z1})}}}{2 \cdot \delta_y} \cdot \frac{L_y}{L_z}. \qquad (34) $$

Finally, the reciprocal Poisson's ratios $\tilde{\mathbf{V}}_{yx}$, $\tilde{\mathbf{V}}_{zx}$ and $\tilde{\mathbf{V}}_{zy}$ must verify the symmetry of the global elasticity matrix [11], such that

$$ \frac{\tilde{\mathbf{V}}_{ij}}{\tilde{\mathbf{E}}_i} = \frac{\tilde{\mathbf{V}}_{ji}}{\tilde{\mathbf{E}}_j}. \qquad (35) $$

The shear moduli $\tilde{G}_{zy}$, $\tilde{G}_{zx}$ and $\tilde{G}_{xy}$ may be found with shear simulations. However, some attention must be paid to defining shear in the case of non-isotropic materials and particularly for laminated composites. Although an orthotropic structure's shear moduli are often estimated without explicit consideration of either sliding or transverse shear configurations ([10, 4, 11]), it may be observed in practice that the two behaviours are not equivalent in general. This is why, in this document, the analysis separates sliding shear (illustrated on Fig. 3) from transverse shear (illustrated on Fig. 4) in the respective “Ortho1” and “Ortho2” scenarios. It may be noted that the identification of Young's moduli and Poisson's ratios remains identical in both scenarios.

The next paragraphs detail the determination of shear moduli by sliding shear schemes, which therefore correspond to the “Ortho1” scenario (see Fig. 3). The simulation combines an enforced displacement $\delta_y$ applied along $-y$ at nodes 1, 2, 3, 4 (face $z=0$), the same displacement $\delta_y$ along $+y$ at nodes 5, 6, 7, 8 (face $z=1$) and plane contact constraints on the same faces $z=0$ and $z=1$ in order to generate pure shear. To stabilise the system, the $T_x$ degree of freedom at node 1 is fixed.

The stress $\sigma_{zy}$ is defined by the relation

$$ \sigma_{zy} = \frac{1}{A_z} \sum_{i \in \mathcal{U}_{z1}} \frac{F_{y,i}}{2}, \qquad (36) $$

where $\mathcal{U}_{z1}$ is the set of the nodes located on the face $z=1$. Finally, this yields the shear modulus $\tilde{G}_{zy}$ [37]:

$$ \tilde{G}_{zy} = \sigma_{zy} \cdot \frac{L_z}{\delta_y}, \qquad (37) $$

where the ratio $L_z/\delta_y$ is related to the shear strain $\epsilon_{zy}$ by
<figure><img src="image_5.png" /><figcaption>Figure 3: Pure sliding shear z-y (Ortho1)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4: Pure transverse shear y-z (Ortho2)</figcaption></figure>

$$ \frac{L_z}{\delta_y} = \frac{1}{\gamma_{zy}} = \frac{1}{2 \cdot \varepsilon_{zy}}. \qquad (38) $$

The two other shear moduli are then computed in a similar way and are defined by the following expressions:

$$ \tilde{G}_{zx} = \frac{L_z}{A_z \cdot \delta_x} \sum_{i \in \mathcal{U}_{z1}} \frac{F_{x,i}}{2} \qquad (39) $$

and

$$ \tilde{G}_{xy} = \frac{L_x}{A_x \cdot \delta_y} \sum_{i \in \mathcal{U}_{x1}} \frac{F_{y,i}}{2}. \qquad (40) $$
### 3. Comparison with other existing methods

The results of the “orthotropic” method presented in Section 2 have been compared to those obtained from some other methods that could be applied to such structures: three isotropic layers stacked along the z-axis, whose properties are gathered into Table 2. The structure’s base cell is a 1-centimetre-long cube for which the layers are organised as following (for increasing z): titanium ($l_1$ = 0.4cm), polypropylene ($l_2$ = 0.2cm) and steel ($l_3$ = 0.4cm), where $l_i$ represents the layer’s thickness.

<table><caption>Table 2: Initial stack</caption><thead><tr><th></th><th>Steel</th><th>Polypropylene</th><th>Titanium</th></tr></thead><tbody><tr><td><i>E</i> [GPa]</td><td>207</td><td>2.0</td><td>121</td></tr><tr><td><i>v</i> [–]</td><td>0.25</td><td>0.40</td><td>0.34</td></tr><tr><td><i>ρ</i> [kg m<sup>-3</sup>]</td><td>7875</td><td>1200</td><td>4430</td></tr><tr><td><i>l</i> [cm]</td><td>0.40</td><td>0.20</td><td>0.40</td></tr></tbody></table>

The other methods being compared are the homogenisation algorithm developed in [9] and the so-called “weighted averages”, that is to say the two approaches detailed in Section 1. A larger structure was built from this base cell, standing for a reference for comparisons. The global structure is shaped as a cuboid, has 5,024 elements, 30,144 DOFs; its dimensions along *x*, *y* and *z* are respectively 210 mm, 110 mm and 60 mm, so that the unit cell is reproduced periodically in every direction.

Also, a node-to-ground 3-D stiffness element is linked to each of the global cuboid’s 8 outer nodes, with stiffness values of $K^{\text{ntg}} = 10^7 \text{ N} \cdot \text{m}^{-1}$ on every direction (*x*, *y* and *z*). The global laminated structure taken as a reference is illustrated on Fig. 5.

<figure><img src="image_7.png" /><figcaption>Figure 5: Global laminated reference structure</figcaption></figure>

The homogeneous finite-element models based on each of the equivalent materials have the same dimensions and the same total mass as the reference structure and are made of 1-centimetre-long cubic 8-node solid elements.
The elastic coefficients of the equivalent materials respectively computed with each of the three methods are compared in Table 3. “Ortho1” and “Ortho2” refer to the scenarios developed in Section 2, “INRIA” to the homogenisation algorithm detailed in [9], and “WA” to the weighted averages. Also, the elasticity matrix $\tilde{\mathbf{C}}^{\text{O1}}$ equivalent to the material computed with method “Ortho1” corresponds to:

$$ \tilde{\mathbf{C}}^{\text{O1}} = \begin{bmatrix} 164 & 59.3 & 9.63 & 0 & 0 & 0 \\ & 164 & 9.63 & 0 & 0 & 0 \\ & & 19.4 & 0 & 0 & 0 \\ & & & 3.34 & 0 & 0 \\ \text{sym.} & & & & 3.34 & 0 \\ & & & & & 51.3 \end{bmatrix} \cdot 10^9 . \quad (41) $$

For each case, the mass matrix of the homogeneous equivalent structure is computed from the relation (2) (weighted averages). Apart from the reference structure, every mass matrix is therefore identical. Finally, the row “CC” compares the computation costs (in seconds) for each of the methods used. The finite-element solver NASTRAN™ v. 2013 is used for scenarios “Ortho1” and “Ortho2” (with a dual-core 2.5 GHz processor and 16 GB of RAM), whereas results for “INRIA” and “WA” are computed with Scilab v. 5.4.1.

<table><caption>Table 3: Comparison of elastic coefficients</caption><thead><tr><th></th><th>Ortho1 (sliding shear)</th><th>Ortho2 (transverse shear)</th><th>INRIA</th><th>WA</th></tr></thead><tbody><tr><td>$\tilde{E}_x$ [GPa]</td><td>132</td><td>132</td><td>137</td><td>132</td></tr><tr><td>$\tilde{E}_y$ [GPa]</td><td>132</td><td>132</td><td>137</td><td>132</td></tr><tr><td>$\tilde{E}_z$ [GPa]</td><td>18.7</td><td>18.7</td><td>19.0</td><td>9.50</td></tr><tr><td>$\tilde{G}_{zy}$ [GPa]</td><td>3.34</td><td>51.3</td><td>3.40</td><td>3.40</td></tr><tr><td>$\tilde{G}_{zx}$ [GPa]</td><td>3.34</td><td>51.3</td><td>3.40</td><td>3.40</td></tr><tr><td>$\tilde{G}_{xy}$ [GPa]</td><td>51.3</td><td>51.3</td><td>51.3</td><td>50.0</td></tr><tr><td>$\tilde{v}_{yz}$ [-]</td><td>0.326</td><td>0.326</td><td>0.313</td><td>0.023</td></tr><tr><td>$\tilde{v}_{xz}$ [-]</td><td>0.326</td><td>0.326</td><td>0.313</td><td>0.023</td></tr><tr><td>$\tilde{v}_{xy}$ [-]</td><td>0.343</td><td>0.343</td><td>0.338</td><td>0.316</td></tr><tr><td>$\tilde{\rho}$ [kg m<sup>-3</sup>]</td><td>5,160</td><td>5,160</td><td>5,160</td><td>5,160</td></tr><tr><td>CC [s]</td><td>8</td><td>8</td><td>&lt;1</td><td>&lt;1</td></tr></tbody></table>

It can be seen that the computation costs are all below 10 seconds, although they do not account for the time necessary to prepare the solutions. Independently from this preparation time, the presented method can be automatised and be adapted to any types of structures, at the cost of quick changes to perform (e.g. element thicknesses).
### 3.1. Modal correlation between the two structures

A convenient and accurate way to evaluate the extent to which a given structure is able to recreate the dynamic behaviour of a reference is to compare vibration modes. These modes constitute a so-called modal basis and describe how and at which frequencies the structure generates a resonance in response of applied loads. Being able to compute the modal basis with accuracy is thus necessary to ensure a good representativity of the finite-element model for dynamic simulations.

The comparison of two modal bases (usually truncated to the first, low-frequency modes) can be made according to several criteria [38] and is part of model updating processes. One convenient way to estimate such modal correlation is to compare each vibration mode of the first model to each of the second. Therefore, two matrices are computed: $\Delta f$ and $\mathbf{MAC}$, of size $M_1 \times M_2$, and where $M_1$ and $M_2$ respectively stand for the numbers of vibration modes considered for the first and the second structures. Each component $\Delta f(m_{e,i}, m_{a,j})$ of the matrix $\Delta f$ expresses the relative difference between the natural frequency of the first structure's i-th vibration mode $m_{1,i}$ and the second structure's j-th mode $m_{2,j}$, and is defined by the relation:

$$ \Delta f(m_{1,i}, m_{2,j}) = \frac{f_{1,i} - f_{2,j}}{f_{1,i}}, \qquad (42) $$

where $f_{1,i}$ et $f_{2,j}$ are the natural frequencies respectively corresponding to the modes $m_{1,i}$ and $m_{2,j}$. The second matrix, $\mathbf{MAC}$, expresses the similarities between the deformed shapes of the modes $m_{1,i}$ and $m_{2,j}$ (respectively called $\phi_{1,i}$ and $\phi_{2,j}$), according to the so-called MAC criterion (Modal Assurance Criterion). Its components $\mathbf{MAC}(m_{1,i}, m_{2,j})$ are defined by the expression [39]:

$$ \mathbf{MAC}(m_{1,i}, m_{2,j}) = \frac{|\{\phi_{1,i}\}^{\top} \{\phi_{2,j}\}|^2}{\{\phi_{1,i}\}^{\top} \{\phi_{1,i}\} \{\phi_{2,j}\}^{\top} \{\phi_{2,j}\}}. \qquad (43) $$

With these notions, two models perfectly correlated are defined by a matrix $\Delta f$ in which every diagonal component is at 0%, and a matrix $\mathbf{MAC}$ in which every diagonal term is at 100%, and the others at 0%. Finally, the pairs of vibration modes for which MAC values are highest are assembled, and are taken into account for the correlation if the MAC values are above a fixed threshold.

### 3.2. Results

The correlation between each of the three homogeneous structures (to which the materials computed by each method have been applied) and the reference model have been compared. To do this, the first 50 vibration modes have been computed for each case, amongst which the first 6 modes (between 413 Hz and 853 Hz) describing the “suspension” related to the node-to-ground elements have been discarded. For the reference structure, the 7th mode (and thus the first to be correlated) is at 3,073 Hz, while amongst all homogeneous structures, the 7th mode of lowest frequency is at 2,947 Hz. The paired modes for which MAC values were below 70% were discarded too. An example of paired deformed shapes between the reference structure and the homogeneous model (computed with scenario “Ortho1”) is given on Fig. 6.

Structure (a) shows the deformation of the individual layers while structure (b) shows “flatter” boundaries due to bigger elements made of identical material properties. The two structures are paired in (c) even though the elements are of different sizes and numbers in (a) and (b).

Then, the results of the correlation are plotted on Fig. 7 and summed up in Table 4 for the first 44 modes of the reference structure. For $N_{pm}$ mode pairs in a given correlation, the entities $|\Delta f|$ and $\overline{\mathbf{MAC}}$ are defined by the expressions:
<figure><img src="image_8.png" /><figcaption>Figure 6: Compared models’ deformed shapes ((a): reference, (b): homogeneous, (c): mode pair)</figcaption></figure>

$$ \overline{|\Delta f|} = \frac{1}{N_{pm}} \cdot \sum_{q=1}^{N_{pm}} |\Delta f(m_1^q, m_2^q)| \qquad (44) $$

and

$$ \overline{MAC} = \frac{1}{N_{pm}} \cdot \sum_{q=1}^{N_{pm}} MAC(m_1^q, m_2^q), \qquad (45) $$

where $m_1^q$ and $m_2^q$ are the modes composing the q-th pair. The deformed shapes of the first 9 paired modes used to illustrate the results on Fig. 7 are detailed in the Appendix, on Fig. A.1.

It can be clearly seen that the method “Ortho1” proposed in Section 2 is more accurate than the weighted averages, and of equivalent quality to “INRIA” in such a setting, whereas using “Ortho2” leads to more significant errors. Although for this stacking setting, the values of $|\Delta f|$ are all relatively high, still 39 modes of the 44 are identified by the method “Ortho1”.

Finally, “Ortho1” presents a new opportunity, that none of the approaches analysed in Section 1 did. Indeed, this type of static simulation is not restrained to applications on finite-element models including distinctive elements and materials, but may also be used with structures on which no information about constitutive materials is available. For example, a superelement whose master-nodes describe the geo-
<figure><img src="image_9.png" /><figcaption>Figure 7: Comparison of frequency differences with reference model (MAC > 70%) - captions “INRIA”, “WA”, “Ortho1” and “Ortho2” respectively refer to the corresponding identification methods</figcaption></figure>

metry of a cuboid could be analysed with the same method, leading to the construction of an equivalent homogeneous material that would recreate the superelement's stiffness behaviour.

## 4. Influence of stacking settings in a beam

This second application aims at analysing the influence of the number of layers in a cantilever beam. For this, four FE models of laminated cantilever beams are created from the same base cell: three layers of isotropic materials, either stacked along the beam's length (models 1 and 2) or in a transverse direction (models 3 and 4). Beam 1 is composed of more layers than beam 2, but of similar thicknesses: beam 2 is shorter. To the contrary, the global dimensions of beams 3 and 4 are identical, while beam 4 is made of more layers: each of them is therefore thinner. Finally, each beam is spanned along z, and its dimensions along x and y are respectively 50 and 30 mm. The four configurations are shown on Fig. 8, and detailed in Table 5.

The properties of the base cell are gathered in Table 6, in which E stands for the material's Young's modulus, v for its Poisson's ratio and χ for the layer's volume fraction (see Eq. (1)).

Applying the method “Ortho1” to this base cell yields equivalent material properties to models 1 and 2 (layers stacked along z):

$$ \tilde{E}_x = \tilde{E}_y = 139 \text{ GPa}, $$

$$ \tilde{E}_z = 14.6 \text{ GPa}, $$

$$ \tilde{G}_{zy} = \tilde{G}_{zx} = 3.20 \text{ GPa}, $$

$$ \tilde{G}_{xy} = 53.7 \text{ GPa}, $$

$$ \tilde{V}_{yz} = \tilde{V}_{xz} = 0.333, $$
<table><caption>Table 4: Correlation of the reference structure's first 44 modes above 853 Hz</caption><thead><tr><th></th><th>Ortho1 (sliding shear)</th><th>Ortho2 (transverse shear)</th><th>INRIA</th><th>WA</th></tr></thead><tbody><tr><td>Nb paired modes</td><td>39</td><td>21</td><td>37</td><td>29</td></tr><tr><td>Δf [%]</td><td>9.29</td><td>58.6</td><td>9.86</td><td>15.3</td></tr><tr><td>MAC [%]</td><td>93.1</td><td>83.7</td><td>93.8</td><td>89.9</td></tr></tbody></table>

<table><caption>Table 5: Details of laminated beams</caption><thead><tr><th>Configuration ID</th><th>Nodes</th><th>Elements</th><th>DOFs</th><th>Stratification along</th><th>Dimension along z [mm]</th><th>Total number of layers</th></tr></thead><tbody><tr><td>1</td><td>4,776</td><td>2,970</td><td>28,656</td><td>z</td><td>1,200</td><td>198</td></tr><tr><td>2</td><td>1,032</td><td>630</td><td>6,192</td><td>z</td><td>252</td><td>42</td></tr><tr><td>3</td><td>4,824</td><td>3,000</td><td>28,944</td><td>y</td><td>1,200</td><td>3</td></tr><tr><td>4</td><td>37,386</td><td>30,000</td><td>224,316</td><td>y</td><td>1,200</td><td>30</td></tr></tbody></table>

and

$$ \tilde{V}_{xy} = 0.293. $$

As the stacking sequences are identical in all cases, the equivalent material properties for models 3 and 4 (layers stacked along y) are defined by the same elastic constants, where indices y and z are switched.

Each of the four beams is associated with a homogeneous one of identical dimensions, but made of only one, equivalent material. A base of the first 72 vibration modes above 0 Hz is then computed for these 8 finite-element models. The modal bases of the initial models are compared to those of the corresponding equivalent structures using the same criteria as described in Subsection 3.1, with a MAC threshold at 0% for pairing the modes (so that all the modes are paired and taken into account). The results are gathered in Table 7.

Table 7 shows that the correlation is globally much better than for Section 3 (in which the objective was to compare the method with other existing ones). Also, the values of |Δf| are significantly lower for the beams with many layers (models 1 and 4) than for the others, showing that in these cases, the equivalent materials are able to recreate the initial structure's behaviour with good precision. Observing good overall
<figure><img src="image_10.png" /><figcaption>Figure 8: Configurations of laminated beams (a): configuration 1, (b): configuration 2, (c): configuration 3, (d): configuration 4)</figcaption></figure>

<table><caption>Table 6: Properties of the beam's base cell</caption><thead><tr><th></th><th>Steel</th><th>Epoxy</th><th>Steel</th></tr></thead><tbody><tr><th>E [GPa]</th><td>207</td><td>3.0</td><td>207</td></tr><tr><th>ν [-]</th><td>0.25</td><td>0.37</td><td>0.25</td></tr><tr><th>χ [-]</th><td>1/3</td><td>1/3</td><td>1/3</td></tr></tbody></table>

correlation for models 1 and 4 and poor coefficients for models 2 and 3 is in agreement to the remark of Lukkassen et al. [14], stating that their homogenisation expression is closer to the initial, heterogeneous behaviour when the dimensions of the unit cell are small in comparison to the entire structure's.

Eventually, the homogeneous material able to recreate the initial structure's modal behaviour with the best accuracy corresponds to the case of stratification in a direction transverse to the beam's length.

## 5. Electric machine stators: experimental-numerical application

Modelling the stator of an electric machine, which is generally a laminated steel stack, is a very interesting application for the studies of heterogeneous structures. Understanding the dynamic behaviour of an electric machine stator is a key issue in the prediction of the electric machine's noise and vibration simulation and prediction [2]. As outlined in Section 1, a stator is built from a multi-laminated magnetic core, on which windings are placed. This section will detail the simulation of a magnetic core's modal behaviour with the aid of the above-presented material property identification method, and will propose several modelling guidelines for this type of structures.

### 5.1. *Experimental data*

The structure studied in this section is the magnetic core of a 12-tooth switched-reluctance machine (without windings or resin). It consists in 400 steel sheets of thickness 360 µm, separated from each other by 3-µm-thick varnish (epoxy) layers. Its dimensions are 154 mm (length) and 245 mm (outer diameter).
<table><caption>Table 7: Correlation of each beam's first 72 vibration modes</caption><thead><tr><th></th><th>1</th><th>2</th><th>3</th><th>4</th></tr></thead><tbody><tr><td>Nb paired modes</td><td>72</td><td>72</td><td>72</td><td>72</td></tr><tr><td>|Δf| [%]</td><td>3.19</td><td>11.1</td><td>9.51</td><td>1.38</td></tr><tr><td>MAC [%]</td><td>98.9</td><td>59.6</td><td>77.3</td><td>99.5</td></tr></tbody></table>

By the means of a shock-hammer analysis, frequency response functions are measured for each of the experimental mesh's 108 degrees of freedom (36 points of 3 DOFs). From these response functions, 5 purely radial vibration modes could be extracted (amongst others). These modes, sometimes referred to as “ovalisation” modes are particularly critical for the acoustic behaviour of the entire electric machine [3]; predicting them with good accuracy is therefore of particular interest. Pictures of the experimental setting and mesh are shown on Figure 9 and Figure 10, respectively. This experimental modal basis stands for a reference in this section.

<figure><img src="image_11.png" /><figcaption>Figure 9: Experimental setting</figcaption></figure>

## 5.2. Finite-element model

Although it may differ from one type of machine to another, the production process for this magnetic core consisted in coating steel layers with insulating varnish (in order to prevent eddy currents from taking place in the structure and therefore dissipating energy) and then simply piling them one onto the others. The stack was then placed under a press, and weld beads were applied on the lateral face while the pressure was maintained. Two examples of finite-element models representing this structure are shown on Figure 11.

The “initial” model details the structure as it really is: a stack of 400 isotropic steel layers separated from each other by isotropic epoxy layers, and is made of 493,164 elements and 618,426 nodes. The details of the layers is given in Table 8. Even with powerful computational resources, performing simulations
<figure><img src="image_12.png" /><figcaption>Figure 10: Experimental mesh (undeformed)</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 11: Finite-element models of the magnetic core (axis along z)</figcaption></figure>

on such a structure is time-consuming and therefore cost-prohibitive. This is why another finite-element model is generated, from the same mesh in the top face, but whose elements are extruded with the same dimensions throughout the structure's length (axis z). The new FE model is called “homogeneous” as it is made of only one homogeneous material (to be determined in the following paragraphs) throughout axis z instead of the details of the layers. Its 19,158 elements and 24,768 nodes make it more appropriate than “initial” for computing modal bases. Also, due to the revolution symmetry of the structure, the finite-element models as well as the material properties are expressed in a cylindrical coordinate system of directions r, θ and z, where z is the stacking direction of the layers.

### 5.3. Modelling guidelines

The weld beads applied on the stack's lateral face are necessary to hold the entire structure in one piece, by imposing a mechanical bond onto the whole length. In the rest of the structure, the sheets are bound together only by the varnish. The pressure applied to the stack is maintained after manufacturing near the weld beads and decreases with the distance in the rest of the structure. Although they are only a few micrometres thick and mechanical properties may vary from one varnish type to another, the local
<table><caption>Table 8: Details of the “initial” model’s layers</caption><thead><tr><th></th><th>Steel</th><th>Epoxy</th></tr></thead><tbody><tr><td><i>E</i> [GPa]</td><td>207</td><td>3.0</td></tr><tr><td><i>v</i> [−]</td><td>0.25</td><td>0.37</td></tr><tr><td><i>ρ</i> [kg·m<sup>−3</sup>]</td><td>7,875</td><td>1,300</td></tr><tr><td><i>l</i> [µm]</td><td>360</td><td>3</td></tr></tbody></table>

behaviour of the varnish layers and thus the interaction between the steel sheets are very likely to be dependent on the residual pressure in the stack, and therefore on the proximity to the weld beads (where the pressure is maintained) or to the free edges. This explains the necessity of dividing the entire FE model into several zones.

The analysis proposed for dividing the model efficiently is the computation of static stiffness at each of the “initial” model’s top-face nodes, in order to estimate the tightness of the structure in relation to the position on the face and the distance to the weld beads. The entire model is clamped (DOFs $Tr$, $Tθ$ and $Tz$ fixed as no rotations are considered on 3D elements) at the nodes on the weld beads. Displacements $Δl_z$ (along $z$) are computed in response of a 1-N static force applied along $z$ at each node (except on the weld beads). The static stiffness values $\mathcal{K}_z$ are then found by:

$$ \mathcal{K}_z = \frac{1}{|\Delta l_z|} \qquad (46) $$

The results are plotted on Figure 12.

<figure><img src="image_14.png" /><figcaption>Figure 12: Distribution of $\mathcal{K}_z$ (in N·m⁻¹)</figcaption></figure>

Judging from the results shown on Figure 12, four zones have been drawn according to the stiffness values:

1. zone “prox” (elements near the weld beads, with $1.5 \times 10^8 < \mathcal{K}_z < 4.6 \times 10^8$ N·m⁻¹);
2. zone “yoke” (with $1.3 \times 10^8 < \mathcal{R}_z < 3.7 \times 10^8 \text{ N} \cdot \text{m}^{-1}$);

3. zone “teeth-W” (elements on the teeth next to the weld beads, with $7.0 \times 10^7 < \mathcal{R}_z < 3.8 \times 10^8 \text{ N} \cdot \text{m}^{-1}$);

4. zone “teeth-Y” (elements on the teeth next to the yoke, with $6.7 \times 10^7 < \mathcal{R}_z < 3.6 \times 10^8 \text{ N} \cdot \text{m}^{-1}$).

Therefore, the “homogeneous” model should be divided into corresponding zones. It seems also relevant to associate zones 3 and 4 into one same “teeth” material definition. Finally, the finite-element model of the magnetic core which will be used for the dynamic simulations is shown on Figure 13, where each zone corresponds to specific material properties.

<figure><img src="image_15.png" /><figcaption>Figure 13: Zoning the magnetic core’s finite-element model (axis along z)</figcaption></figure>

It seems important to note that the entire model is composed of 3D, solid elements, which need not necessarily be cuboidal. Although the identification method is applied to base cells made of 8-node cuboid elements (particularly convenient for modelling and homogenising stacked structures), meshing the equivalent structure remains under the control of the user: the equivalent material may be applied to any type of mesh, as long as the global geometry is kept unchanged. Therefore, the stator could have been meshed with any type of solid element (e.g. four-, five- or six-sided), even mixed, as this choice has negligible influence over the global structure’s vibrational behaviour. In this case, creating the mesh with a base of six-sided elements is a good compromise between representativity and total number of DOFs.

### 5.4. Equivalent material properties

The first zone under the spotlight is “prox”, gathering the elements which are closest to the weld beads. In this zone, the proximity to the weld beads ensures the stator’s tightest cohesion between the steel sheets, and implies therefore the best regularity in the successive varnish layers’ thicknesses. This means that the cell <half-thick steel sheet; varnish sheet; half-thick steel sheet> is repeated regularly through the stator’s length (dimension along z), with the same thicknesses everywhere. The corresponding finite-element model which will be used for material identification is therefore a 3-layered cuboid with a square base of length 400 µm, as illustrated on Figure 14. Its layers’ Young’s moduli E, Poisson’s ratios v and thicknesses l are gathered in Table 9. In the whole FE model, the same density $\tilde{\rho}$ is applied to all the zones and directly taken from the measurements on the stator:

$$ \tilde{\rho} = \frac{m^{\text{tot}}}{V^{\text{tot}}} = 7,750 \text{ kg} \cdot \text{m}^{-3}, $$
<figure><img src="image_16.png" /><figcaption>Figure 14: Base cell</figcaption></figure>

<table><caption>Table 9: Properties of the base cell of zone “prox”</caption><thead><tr><th></th><th>Steel</th><th>Epoxy</th><th>Steel</th></tr></thead><tbody><tr><th>E [GPa]</th><td>207</td><td>3.0</td><td>207</td></tr><tr><th>v [-]</th><td>0.25</td><td>0.37</td><td>0.25</td></tr><tr><th>l [µm]</th><td>180</td><td>3</td><td>180</td></tr></tbody></table>

where $m^{tot}$ and $V^{tot}$ respectively refer to the stator’s total mass and volume.

Applying method “Ortho1” to the base cell of zone “prox” yields an equivalent material, whose elastic constants are:

$$ \tilde{E}_r^{prox} = \tilde{E}_\theta^{prox} = 205 \text{ GPa}, $$

$$ \tilde{E}_z^{prox} = 157 \text{ GPa}, $$

$$ \tilde{G}_{z\theta}^{prox} = \tilde{G}_{zr}^{prox} = 51.2 \text{ GPa}, $$

$$ \tilde{G}_{r\theta}^{prox} = 82.1 \text{ GPa} $$

and

$$ \tilde{V}_{\theta z}^{prox} = \tilde{V}_{rz}^{prox} = \tilde{V}_{r\theta}^{prox} = 0.25. $$

Concerning the other zones, to the authors’ knowledge there does not exist any analytical or numerical method able to describe with precision the variation of the elastic properties with the distance to the weld
beads. This is why a good experience in model updating of electric machine stators is necessary to describe the behaviours of the different zones. In this case, the factors 3/4 and 1/2 have been applied to the shear moduli of the respective zones “yoke” and “teeth”, according to their distance to the weld beads. The other coefficients remain unchanged. Eventually, the material properties are detailed in Table 10.

<table><caption>Table 10: Comparison of elastic coefficients</caption><thead><tr><th></th><th>prox</th><th>yoke</th><th>teeth</th><th>weld beads</th></tr></thead><tbody><tr><td>E&#773;<sub>r</sub> [GPa]</td><td>205</td><td>205</td><td>205</td><td rowspan="3">207</td></tr><tr><td>E&#773;<sub>&theta;</sub> [GPa]</td><td>205</td><td>205</td><td>205</td></tr><tr><td>E&#773;<sub>z</sub> [GPa]</td><td>157</td><td>157</td><td>157</td></tr><tr><td>G&#773;<sub>z&theta;</sub> [GPa]</td><td>51.2</td><td>38.4</td><td>25.6</td><td rowspan="3">-</td></tr><tr><td>G&#773;<sub>zr</sub> [GPa]</td><td>51.2</td><td>38.4</td><td>25.6</td></tr><tr><td>G&#773;<sub>r&theta;</sub> [GPa]</td><td>82.1</td><td>61.6</td><td>41.1</td></tr><tr><td>v&#773;<sub>&theta;z</sub> [-]</td><td>0.25</td><td>0.25</td><td>0.25</td><td rowspan="3">0.288</td></tr><tr><td>v&#773;<sub>rz</sub> [-]</td><td>0.25</td><td>0.25</td><td>0.25</td></tr><tr><td>v&#773;<sub>r&theta;</sub> [-]</td><td>0.25</td><td>0.25</td><td>0.25</td></tr><tr><td>p&#773; [kg m<sup>-3</sup>]</td><td>7,750</td><td>7,750</td><td>7,750</td><td>7,750</td></tr></tbody></table>

## 5.5. Comparison with experimental data

Computing the FE model’s modal basis and comparing it to the experimental data with the criteria described in Subsection 3.1 leads to the correlation state presented in Figure 15 and Table 11, where the columns “FEA” and “EMA” respectively refer to the mode frequencies in the FE model and in the experimental modal basis. The $|\Delta f|$ and $MAC$ averages have been computed in the bottom line (cf. definitions in Subsection 3.2) of Table 11. Only mode pairs for which MAC values were over 60% were taken into account.

Judging from these results, it can be observed that the behaviour of the finite-element model computed with the method “Ortho1” is in good agreement with the measured natural frequencies, and that the model’s deformed shapes correspond to the experimental ones fairly well.

The results of Table 11 are compared with the case “NZ” (no zoning: the same “prox” material properties applied to the whole structure), as shown in Table 12.

As in the previous sections, it can be seen that method “Ortho1” generates equivalent material properties that simulate the dynamic behaviour of the entire structure accurately. The necessity of zoning the model and adapting the shear coefficients according to the elements’ distances to the weld beads can be clearly seen as well, as the average frequency discrepancy is increased by 68.5%. The modelling guidelines that have been proposed should be therefore followed for other geometries of magnetic core structures.

In this application especially, the initial correlations of simulated and experimental modal bases on two different electric machine stators were particularly promising for the prediction of such structures’ modal behaviours. Judging from the time needed for the thorough experimental-based updating of an electric
<figure><img src="image_17.png" /><figcaption>Figure 15: Correlation of FE and experimental modal bases - mode shape pairs (experimental frequencies) - blue lines: simulation, red dots: measurement</figcaption></figure>

machine stator's FE model, the stakes of an accurate, initial modelling of such structures in a relatively simple and fast way are both high and numerous [3].

## 6. Conclusion

In this document, a new method for identifying equivalent material properties to laminated orthotropic structures was developed. As a necessity in order to reduce the number of degrees of freedom in the finite-element models of heterogeneous structures such as electric motor stators, the identification method enabled modelling them with equivalent representative homogeneous material properties. This new method has been applied to various domains involving multi-layered material simulations. It has been shown that analysing laminated structures requires a “sliding shear” approach rather than a “transverse shear” one. The corresponding “Ortho1” approach is globally more accurate than existing reference homogenisation methods for such structures, yields good simulation results when used in an experimental-numerical application, and unlike the others, can also be applied to superelements. This finite-element-based identification method also allows quick parametric sensitivity computations on heterogeneous structures.

An analysis of stacking settings in laminated cantilever beams has shown that the equivalent orthotropic material created with this new technique was able to recreate the initial vibrational behaviour with good precision for multi-layered structures. Eventually, the vibrational behaviour of the magnetic core of an electric motor stator has been simulated by the means of a finite-element model with equivalent homogeneous material properties. The simulated ovalisation modes have been compared to experimental data measured on a real stator, and showed good agreement. In addition to that, the zoning guidelines proposed for finite-element models of magnetic cores have been proven efficient. It is hoped that this identifica-
<table><caption>Table 11: Correlations of FE and experimental modal bases</caption><thead><tr><th>Mode pair</th><th>FEA [Hz]</th><th>EMA [Hz]</th><th>Δf [%]</th><th>MAC [%]</th><th>Mode description</th></tr></thead><tbody><tr><td>1</td><td>823.5</td><td>798.6</td><td>3.12</td><td>63.8</td><td>cylinder mode, order 2</td></tr><tr><td>2</td><td>2179.6</td><td>2125.5</td><td>2.55</td><td>62.3</td><td>cylinder mode, order 3</td></tr><tr><td>3</td><td>3855.0</td><td>3747.7</td><td>2.86</td><td>90.6</td><td>cylinder mode, order 4</td></tr><tr><td>4</td><td>5529.3</td><td>5330.1</td><td>3.74</td><td>82.5</td><td>cylinder mode, order 5</td></tr><tr><td>5</td><td>6593.4</td><td>6286.5</td><td>4.88</td><td>74.3</td><td>cylinder mode, order 0</td></tr><tr><td>Averages |Δf| and MAC</td><td></td><td></td><td>3.43</td><td>74.7</td><td></td></tr></tbody></table>

<table><caption>Table 12: Influence of zoning</caption><thead><tr><th></th><th>Ortho1</th><th>NZ</th></tr></thead><tbody><tr><td>|Δf| [%]</td><td>3.43</td><td>5.78</td></tr><tr><td>MAC [%]</td><td>74.7</td><td>74.5</td></tr></tbody></table>

tion method will eventually help improve further the current knowledge about the behaviour of laminated structures and in particular of electric machine stators. Also, substituting time-consuming model updating procedures by a simple, direct method could even replace the manufacture of costly prototypes and associated measurements. Some other modelling aspects which are still based on practical experience will be discussed in future work.

For more complicated situations however, notably in case of non-negligible preloading, or when couplings exist between tension-compression and shear in the global structure, orthotropic equivalent materials might not be accurate enough to simulate the expected behaviour of a given structure. Future work will present another method for the identification of equivalent material properties for anisotropic or triclinic structures and superelements, accounting for the influence of external perturbations such as friction and preload.

## Appendix A. First 9 paired modes

The deformed shapes of the first 9 modes paired in Subsection 3.2 are shown on Fig. A.1. For each mode pair, the reference structure is drawn in red dots, and the homogeneous model in blue solid elements.

## References

[1] A. Kalamkarov, I. V. Andrianov, V. V. Danishevs'kyy, Asymptotic homogenization of composite materials and structures, Applied Mechanics Reviews 62.
<figure><img src="image_18.png" /><figcaption>Figure A.1: First 9 paired modes (reference structure's natural frequencies)</figcaption></figure>
[2] L. Humbert, P. Pellerey, S. Cristaudo, Electromagnetic and structural coupled simulation to investigate NVH behavior of an electrical automotive powertrain, SAE International (12SNVH-0002) (2012) pp. 1–12.

[3] J.-B. Dupont, P. Bouvet, L. Humbert, Vibroacoustic simulation of an electric motor, methodology and focus on the structural FEM representativity, in: International Conference on Electrical Machines, Marseille, pp. 1–7, 2012.

[4] J.-M. Berthelot, Composite materials, mechanical behavior and structural analysis, Springer, 1999.

[5] S. Vel, R. Batra, Three-dimensional exact solution for the vibration of functionally graded rectangular plates, Journal of Sound and Vibration 272 (3) (2004) pp. 703–730.

[6] J. N. Reddy, Mechanics of Laminated Composite Plates and Shells, Theory and Analysis, 2nd Edition, CRC Press, 2004.

[7] E. Fasanella, K. Jackson, K. Lyle, Finite Element Simulation of a Full-Scale Crash Test of a Composite Helicopter, Journal of the American Helicopter Society 47 (3) (2002) pp. 156–168.

[8] M. Koishi, M. Shiratori, T. Miyoshi, K. Kabe, Homogenization method for dynamic viscoelastic analysis of composite materials, JSME International Journal-Series A-Mechanics and Material Engineering 40 (3) (1997) pp. 306–312.

[9] D. Begis, G. Duvaut, A. Hassim, Homogénéisation par éléments finis des modules de comportements élastiques de matériaux composites [Finite-element homogenisation of composite materials’ elastic constants], Tech. Rep. 101, Institut National de Recherche en Informatique et Automatique, 1981.

[10] D. Begis, A. Bestagno, G. Duvaut, A. Hassim, M. Nuc, A new method of computing global elastic moduli for composite materials, Tech. Rep. 195, Institut National de Recherche en Informatique et Automatique, 1983.

[11] Y. Chevalier, J. V. Tuong, Mechanics of viscoelastic materials and wave dispersion, Iste/Wiley, 2010.

[12] A. Puck, Grundlagen der Spannungs- und Verformungs Analyse [Basic notions of tension and distortion analysis], Dipl. Ing. Kunststoffe, Bd 57 (4).

[13] S. W. Tsai, J. C. Halpin, N. J. Pagano, R. T. Schwartz, Composite materials workshop, Technomic, 1968.

[14] D. Lukkassen, L.-E. Persson, P. Wall, Some engineering and mathematical aspects on the homogenization method, Composites Engineering 5 (5) (1995) pp. 519–531.

[15] A. L. Kalamkarov, E. M. Hassan, A. V. Georgiades, M. A. Savi, Asymptotic homogenization model for 3D grid-reinforced composite structures with generally orthotropic reinforcements, Composite structures 89 (2) (2009) pp. 186–196.

[16] E. M. Hassan, A. V. Georgiades, M. A. Savi, A. L. Kalamkarov, Analytical and numerical analysis of 3D grid-reinforced orthotropic composite structures, International Journal of Engineering Science 49 (7) (2011) pp. 589–605.

[17] G.-D. Cheng, Y.-W. Cai, L. Xu, Novel implementation of homogenization method to predict effective properties of periodic materials, Acta Mechanica Sinica 29 (4) (2013) pp. 550–556.
[18] P. W. Chung, K. K. Tamma, R. R. Namburu, Asymptotic expansion homogenization for heterogeneous media - computational issues and applications, Composites Part A - Applied Science and Manufacturing 32 (9) (2001) pp. 1291-1301.

[19] G. Milton, The theory of composites, vol. 6, Cambridge University Press, 2002.

[20] A. Cecchi, K. Sab, A homogenized Reissner-Mindlin model for orthotropic periodic plates - Application to brickwork panels, International journal of solids and structures 44 (18) (2007) pp. 6055-6079.

[21] J. P. Watt, G. F. Davies, R. J. O'Connell, The elastic properties of composite materials, Reviews of Geophysics 14 (4) (1976) pp. 541-563.

[22] M. Kamiński, M. Kleiber, Perturbation based stochastic finite element method for homogenization of two-phase elastic composites, Computers & Structures 78 (6) (2000) pp. 811-826.

[23] N. Magalhães Dourado, J. Cardoso Xavier, J. Lopes Morais, On the identification of linear elastic mechanical behaviour of orthotropic materials using evolutionary algorithms, in: IUTAM Symposium on Evolutionary Methods in Mechanics, Springer, pp. 251-264, 2004.

[24] A. L. Araújo, M. Soares, M. J. Moreira de Freitas, P. Pedersen, J. Herskovits, Combined numerical-experimental model for the identification of mechanical properties of laminated structures, Composite Structures 50 (4) (2000) pp. 363-372.

[25] H. Wang, K. Williams, Effects of laminations on the vibrational behaviour of electrical machine stators, Journal of sound and vibration 202 (5) (1997) pp. 703-715.

[26] S. P. Verma, R. K. Singal, K. Williams, Vibration behaviour of stators of electrical machines, part 1, Theoretical study, Journal of sound and vibration 115 (1) (1987) pp. 1-12.

[27] K. Williams, R. K. Singal, S. P. Verma, Vibrations of long and short laminated stators of electrical machines Part 2, Results for long stators, Journal of sound and vibration 129 (1) (1989) pp. 15-29.

[28] J. Le Besnerais, Reduction of magnetic noise in PWM-supplied induction machines, low-noise design rules and multi-objective optimisation, Ph.D. thesis, Laboratoire d'Électricité et d'Électronique de Puissance de Lille, École Centrale de Lille, 2008.

[29] R. M. Barker, L. Fu-Tien, J. R. Dana, Three-dimensional finite-element analysis of laminated composites, Computers & Structures 2 (5) (1972) pp. 1013-1029.

[30] G. Kuhlmann, R. Rolfes, A hierarchic 3D finite element for laminated composites, International journal for numerical methods in engineering 61 (1) (2004) pp. 96-116.

[31] R. F. S. Hearmon, The elastic constants of anisotropic materials, Reviews of Modern Physics 18 (3) (1946) p. 409.

[32] M. Hayes, A simple statical approach to the measurement of the elastic constants in anisotropic media, Journal of Materials Science 4 (1) (1969) pp. 10-14.

[33] F. Pierron, S. Zhavoronok, M. Grédiac, Identification of the through-thickness properties of thick laminated tubes using the virtual fields method, International Journal of Solids and Structures 37 (32) (2000) pp. 4437-4453.
[34] R. Rikards, A. Chate, W. Steinchen, A. Kessler, A. K. Bledzki, Method for identification of elastic properties of laminates based on experiment design, Composites Part B: Engineering 30 (3) (1999) pp. 279–289.

[35] R. Rikards, A. Chate, G. Gailis, Identification of elastic properties of laminates based on experiment design, International Journal of Solids and Structures 38 (30) (2001) pp. 5097–5115.

[36] R. Rikards, Identification of mechanical properties of laminates, Modern Trends in Composite Lam- inates Mechanics (448) (2004) p. 181.

[37] E. J. Hearn, Mechanics of Materials, vol. 1 - 2, Pergamon Press, 1985.

[38] J. Mottershead, M. Friswell, Model updating in structural dynamics, a survey, Journal of sound and vibration 167 (2) (1993) pp. 347–375.

[39] R. Allemang, The modal assurance criterion - twenty years of use and abuse, Journal of Sound and Vibration 37 (8) (2003) pp. 14–23.