#### Abstract

Based on the symplectic structure of the Hamiltonian matrix, the precise integration method (PIM), and the Wittrick–Williams (W-W) algorithm, a generalized method for computing the dispersion curves of guided waves in multilayered anisotropic magneto-electro-elastic (MEE) structures for different types of mechanical, electrical, and magnetical boundaries is developed. A strictly theoretical analysis shows that the W-W algorithm cannot be applied directly to the MEE structure. This is because a block of the Hamiltonian matrix is not positive definite for MEE structures so that the eigenvalue count of the sublayer is not zero when the divided sublayer is sufficiently thin. To overcome this difficulty, based on the symplectic structure of the Hamiltonian matrix, a symplectic transformation is introduced to ensure that the W-W algorithm can be applied conveniently to solve wave propagation problems in multilayered anisotropic MEE structures. The application of the PIM based on the mixed energy matrix to solve the wave equation can ensure the stability and efficiency of the method, and all eigenfrequencies are found without the possibility of any being missed using the W-W algorithm. This research provides the necessary insight to apply the W-W algorithm in wave propagation and vibration problems of MEE structures.

#### 1. Introduction

Magneto-electro-elastic (MEE) materials, which can achieve the conversion between electrical and magnetic energy due to the coupled effect between electric and magnetic fields, have permeated every aspect of the modern technology. This feature promotes the wide application of the MEE in many fields of science and engineering [1–3], for example, sensors [4], smart devices [5], and nondestructive evaluation [6], which are closely related to the knowledge of wave propagation. Various techniques related to the wave propagation in MEE structures were developed by numerous scholars. Pan et al. [7–10] contributed many efforts to the wave propagation in MEE structures, and they derived a series of analytical expressions for various MEE structures. Ezzin et al. [11] investigated Love waves in a transversely isotropic piezoelectric layer bonded to an MEE semi-infinite space and derived explicit dispersion equations. Li et al. [12] investigated propagation behavior of the Bleusteine–Gulyaev waves in a functionally graded transversely isotropic MEE half-space and derived analytically the dispersion equations under electromagnetically open and shorted conditions. An analytical treatment based on the transfer matrix was presented by Chen et al. [13] for the propagation of harmonic waves in multilayered MEE plates. An exact solution for shear horizontal (SH) waves propagating in a transversely isotropic MEE was derived by Nie et al. [14].

Apart from analytical solutions, abundant numerical methods were also developed for wave propagation in layered MEE structures. Wu et al. [15] explored the dispersive behavior of the symmetric and antisymmetric Lamb waves in an infinite MEE plate. Using the ordinary differential equation and stiffness matrix method, Ezzin et al. [16] investigated the propagation behavior of SH waves in laminated MEE plates, and the effects of thickness ratio on phase velocity and group velocity were discussed. Based on Hamilton’s principle, Xiao et al. [17, 18] applied Chebyshev spectral element method to investigate the dispersion characteristics of guided waves in a functionally graded MEE plate and in a multilayered MEE curved panel, respectively. Chen et al. [19] employed the reverberation matrix method to compute the dispersion curves of functionally graded materials. Their method adopted a homogenous assumption of material parameters in each thin layer because the graded layer was divided into a sufficiently large number of thin layers. Yu et al. [20] employed Legendre orthogonal polynomial series method to study the dispersion behavior of guided waves in a layered MEE plate. Matar et al. [21] used the combined Legendre-Laguerre polynomial expansion method to compute dispersion curves and mode shapes in layered MEE composites. In addition, the study of periodic MEE structures has attracted wide attention in the past decades. Pang et al. [22,23] applied the transfer matrix method and the stiffness matrix method to investigate SH bulk/surface waves propagating in periodically layered infinite/semi-infinite MEE composites. Wang et al. [24] studied the wave propagation in periodic composites consisting of MEE plates using the plane-wave expansion method. Liu et al. [25] studied dispersion behavior and transmission coefficients of SH wave in a periodically layered piezoelectric structure and discussed the feature of bandgaps in periodic structures. Chen et al. [26] investigated dispersions and band structures of elastic waves in nanoscale periodic MEE structures based on the nonlocal theory. The localization factors and dispersion curves were computed using the transfer matrix method. Recently, the band structure and evanescent behavior of Bloch waves in periodic MEE structures were investigated by applying extended plane-wave expansion method and Bloch’s theory [1–3].

These numerical methods can be broadly classified into two groups. One group is based on analytical models, such as the transfer matrix method. In principle, these methods are exact because mathematical expressions are free from approximations, but in numerical calculation here exist many difficulties that are associated with the numerical overflow and the loss of precision at high frequencies. Meanwhile, a disadvantage is that these methods involve a challenging problem of seeking the roots from transcendental eigenequations, which means that its solution is usually obtained by the intensive search technique. The other group is based on discrete models, such as the finite-element method. For these methods, the substructuring technique is usually used to divide each layer into a large number of thin layers, then the displacements are approximate via interpolation polynomials. These discrete methods avoid seeking the roots from the transcendental eigenequation by solving a quadratic eigenvalue problem. In spite of this advantage, it should be highlighted that, due to the need for the discretization, these methods give rise to more degrees of freedom of system than previous methods.

The precise integration method (PIM) was proposed by Zhong [27] to solve accurately the sets of first-order ordinary differential equations with specified two-point boundary value conditions for space domain problems, or with specified initial value conditions for time-domain problems. The combination of the PIM and the Wittrick–Williams (W-W) algorithm [28–31] has used extensively to study wave propagation in layered media [32–36]. Numerous results shown that the method for solving wave propagation problem not only can avoid seeking the roots from the nonlinear transcendental eigenequation directly but also can guarantee that the calculations have good stability and high accuracy. Later, numerous researches show that the method based on the mixed energy matrix was more stable compared to those based on the transfer and stiffness matrices [36]. A key process for the application of the W-W algorithm is to determine the eigenvalue count *J* of the structure, and it is usually difficult to directly solve *J*. To this end, a substructure technique is usually used to divide each layer into sufficient number of sublayers, so that the eigenvalue count of the sublayer with both ends clamped is zero. Then, by gradually condensing sublayers, the eigenvalue count of the whole structure is obtained. However, for a piezoelectric structure, the W-W algorithm is no longer suitable because the eigenvalue count of sublayer with both ends mechanically clamped, electrically short and magnetically short is zero even if each layer is divided into sufficient number of sublayers. For this reason, a symplectic transformation was adopted to deal with guided wave propagation in multilayered anisotropic piezoelectric structures [37]. After performing symplectic transformation, the W-W algorithm combined to the PIM could be applied to compute the eigenfrequencies of waves in the piezoelectric structure. Obviously, similar problems exist in MEE structures as well when the W-W algorithm is applied to compute the wave dispersion. Compared to piezoelectric cases, solutions to the wave propagation problem in MEE structure are much more complicated because of its complex boundary conditions and high-dimensional governing equations. In this paper, based on symplectic transformation, a generalized method for dispersion analysis of waves in multilayered anisotropic MEE structures is developed. The performance of the method is verified by several numerical examples.

#### 2. Statement of the Problem and Basic Equations

Consider a multilayered MEE structure consisting of *l* layers with different material properties and layer thicknesses, as illustrated in Figure 1. The origin of the Cartesian coordinate system is set on the surface, where the *z*-axis is along the depth direction. The structure extends infinitely in the *x* and *y* directions. The *k*th layer is bounded by the interfaces and . The thickness for the *k*th layer is (*k* = 1, 2,…, *l*), and the total thickness is . Without the loss of generality, it is assumed to the wave propagating along the *x-y* plane and the incident angle measured from the positive *x*-axis in the clockwise direction.

For an anisotropic MEE solid, the constitutive equation in linear approximation can be expressed by equation [38].where and are the stress and strain tensors; and are electric displacement and magnetic induction tensors; and are the tensors of the electric and magnetic fields; , , and are the elastic, piezoelectric, and piezomagnetic coefficient tensors; , , and dielectric permittivity, magnetic permeability, and magneto-electric coefficient tensors; and the superscript *T* represents the transpose. In the absence of the body force, electric charge, and magnetic charge, the mechanical, electric, and magnetic governing equations can be expressed bywhere is the displacement vector, *t* denotes time, is the density, and is the nabla operator. The generalized geometric equations are described bywhere and are the electric and magnetic potentials. Note that the detailed forms of the physical quantities and material constants for elastic, electric, and magnetic fields are presented in Appendix A. To ensure stability of the system, the total internal energy is more than zero [39], such that the magneto-electric matrixand elastic coefficient matrix are both positive definite, as shown in Appendix A.

The plane-wave solution is considered aswhere ; and are the generalized displacement vectors in the time-space and frequency-wavenumber fields, respectively; represents the wave circular frequency; and are the two components of the wave vector in the *x* and *y* directions, described separately as and ; and is the magnitude of the wave vector along the propagation direction. Substituting equations (1) and (3) into equation (2) and then using equation (5), we obtainwhere and denote the first and second derivatives of with respect to *z*; and

Here, the superscript H denotes a Hermitian transpose. It can be seen that and are both real symmetric matrices, and and are both pure imaginary matrices.

#### 3. State-Space Formulation and Symplectic Transformation

##### 3.1. State-Space Formulation

Introducing a dual vector satisfied by and defining the state vector , from equation (6) the following equation of state can be given aswhere

It is easy to verify that is a Hamiltonian matrix, satisfied by , where **J** is a unit symplectic matrix defined byand denotes a identify matrix.

In the aforementioned paragraph, the governing equation of the wave motion was converted into a first-order ordinary differential equations, which can be accurately solved using the PIM [27]. Then, combined to the W-W algorithm, all eigenfrequencies of waves in layered elastic media can be calculated [33, 34]. The PIM is used to ensure the precision and stability of computation, and the W-W algorithm is used to ensure that all eigenfrequencies are found without the possibility of any being missed. However, the method could not be applied directly to wave propagation problems in MEE structures. This is due primarily to the fact that a block matrix of the Hamiltonian matrix in equation (10) is not positive definite.

For purely elastic structures, according to the positive definiteness of the elastic coefficient matrix **c**, from the form of in equation (7) and the relation of in equation (11), we conclude that is positive definite. However, for MEE structures, matrices **c** and are both positive definite (see Appendix A), so defined in equation (7) is not positive definite. Then, according to the relation of , it is evident that is not positive definite. The computational formulas of the W-W algorithm require that must be a positive definite matrix. In the following section, a symplectic transformation is introduced to overcome this difficulty.

##### 3.2. Symplectic Transformation

Introducing the following mathematical transformation,in which **T** and **v** can be expressed as

It can be easily confirmed that **T** is satisfied by , so it is a symplectic matrix [40]. Substituting equation (13) into equation (10) gives

According to the properties of the Hamiltonian matrix and symplectic theory [40], if **T** is a symplectic matrix, **H** is also a Hamiltonian matrix, which can be written in the following form:

After performing the aforementioned symplectic transformation, the new matrix **D** will be positive definite. The proof of positive definitiveness of matrix **D** is given in Appendix B.

#### 4. PIM and W-W Algorithm for Multilayered MEE Structures

After performing the symplectic transformation, this section will be devoted to the derivation of formulas for calculating wave propagation problems in multilayered MEE structures based on the PIM and the W-W algorithm. Those formulas may be considered as an extension of the purely elastic structure, and they can be found in previously published literature [33, 34, 36].

The PIM [27] is first used to divide the *k*th layer into sublayers with equal thickness , from which two arbitrary adjacent sublayers and are selected. Then, the solution of equation (10) based on the mixed energy matrix form gives the following equations:in which , , and , , are mixed energy matrices for the sublayer and ; , , and , , are the displacement and dual vectors at the interfaces , and . Using equations (17) and (18) and eliminating and , the following equation for the sublayer is obtained as follows:where

For the sake of brevity, assuming that F and C represent the boundary conditions of mechanically free and clamped, and that O and S represent the boundary conditions of electrically or magnetically open and short. Then, a symbol composed of three letters is used to describe the boundary condition, for example, COO denotes the boundary condition of mechanically clamped, electrically open, and magnetically open. Assume that and are the eigenvalue counts for the sublayers and with both ends mechanically clamped, electrically open, and magnetically open (COO-COO). Then, according to the W-W algorithm [28], the eigenvalue count for sublayer with both ends mechanically clamped, electrically open and magnetically open can be calculated by [29]where is the sign count of the matrix #, which is defined as the number of the negative eigenvalues of the matrix # and is equal to the number of changes in sign of the Sturm sequence [41].

Repeated application of equations (20) and (21) times, the eigenvalue count and mixed energy matrices , and for the *k*th layer can be obtained. Then, again repeated application of equations (20) and (21) *l* − 1 times, the eigenvalue count and global mixed energy matrices , and for all *l* layers can be obtained.

Note that the mixed energy matrices for a sublayer with thickness can be calculated using the Taylor series, and the computational formula can be found from literature [33,36]. Moreover, since the matrix **D** is positive definite, sufficiently large can guarantee for a sublayer with thickness .

#### 5. Computing Eigenfrequencies for Different Boundary Conditions

In the aforementioned section, the procedure for computing the eigenvalue count for a multilayered structure with both ends mechanically clamped, electrically open, and magnetically open has been given. In this section, five types of mechanical, electrical, and magnetical boundary conditions are considered (see Figure 2) and under these boundary conditions, the computational formulas of the eigenvalue count of the multilayered structure are also given. Finally, all eigenfrequencies can be calculated based on the concept of the eigenvalue count.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

On the boundary surfaces and , commonly encountered mechanical, electrical, and magnetical boundary conditions are considered as follows: Case I: The boundaries at both surfaces are mechanically clamped, electrically open, and magnetically open (COO-COO), that is, The eigenvalue count for this case has been given in Section 4, that is, . Case II: The boundaries at both surfaces are mechanically free, electrically open, and magnetically open (FOO-FOO), that is, According to the W-W algorithm, the eigenvalue count for this case is Case III: The boundary at the surface is mechanically clamped, electrically short and magnetically short, and the boundary at the surface is mechanically free, electrically open, and magnetically open (CSS-FOO), that is, According to the W-W algorithm, the eigenvalue count for this case is Case IV: The boundary at the surface is mechanically clamped, electrically short and magnetically short, and the boundary at the surface is mechanically free, electrically short, and magnetically short (CSS-FSS), that is, According to the W-W algorithm, the eigenvalue count for this case is Case V: For a periodically infinite multilayered structure, the state vectors and at both ends of a unit cell satisfy the Bloch’s theorem [42], that is, where is the Bloch wavenumber in periodically infinite multilayered structures. According to the W-W algorithm, the eigenvalue count for this case is

So far, the eigenvalue count for the whole structure with a certain boundary condition has been obtained for a given wavenumber and trial frequency . Then, based on the concept of the eigenvalue count, all eigenfrequencies can be computed using the bisection method [28,36]. Note that the general anisotropic case are derived in this paper, but the other special cases (e.g., antiplane and in-plane wave problems) can be easily recovered by removing some items of wave equation and reducing the size of matrices.

#### 6. Numerical Results and Discussion

Three numerical examples are considered in this section. The focus is on examining the numerical performance of the method. In Section 6.1, an example for a three-layered MEE plate is first presented to verify the correctness of this method, which involves the comparison with previously published results. In Section 6.2, this method is used to compute the dispersion curves for a strongly anisotropic four-layered MEE plate. In Section 6.3, an example is presented to compute the dispersion curves of a periodically multilayered finite MEE structure and bandgaps of a periodically infinite multilayered MEE structure. The parameters of materials BaTiO_{3}, CoFe_{2}O_{4}, Bimaterial A, and Bimaterial B used in examples are listed in Table 1 [8,13,43].

##### 6.1. A Multilayered Transversely Isotropic MEE Structure

A three-layered Sandwich plate with BaTiO_{3}/CoFe_{2}O_{4}/BaTiO_{3} (B/F/B) configuration is considered in this example, and all three layers have equal thickness. The objective of this example is to verify the correctness of this method compared to the published results, and then to examine the effects of varying boundary condition and stacking sequence to the plate on dispersion curves. The wave propagation direction is assumed to be along the *x*-axis, that is, . To facilitate the comparison, the dimensionless wavenumber and dimensionless phase velocity are defined as and ; and and are the maximum elastic constant and maximum density in the system.

The dispersion curves for the B/F/B plate for four different types of boundary conditions are shown in Figure 3, where subfigures (a), (b), (c), and (d) represent the results for the boundary conditions of FOO-FOO, COO-COO, CSS-FOO, and CSS-FSS. Generally, in the case of anisotropic MEE material [44,45], all waves with respect to the mechanical, electrical and magnetical components are coupled and in order to identify them one needs to solve a tenth order polynomial characteristic equation. For special material symmetry directions, such as isotropic or transversely isotropic materials, the wave motion can be decoupled into two sets of independent equations. One of the equations is a second-order polynomial equation corresponding to the antiplane wave with respect to the component . Another equation is an eighth-order polynomial equation corresponding to the in-plane wave with respect to the component *u*, , , and . Since BaTiO_{3} and CoFe_{2}O_{4} are both transversely isotropic materials, the wave motion can be decoupled into two sets of independent equations, corresponding to antiplane wave with respect to the component and in-plane wave with respect to the components *u*, , , and . Figure 3(a) is in good agreement with the results given in Figures 4(c) of Ref. [13], which confirms the correctness of this method. From Figures 3(a) and 3(b), it can be seen that the dispersion curves have significant changes when the boundary condition changes from FOO-FOO to COO-COO. However, compared Figures 3(c) and 3(d), it can be observed that the dispersion curves are only slight variations when the boundary condition changes from CSS-FOO to CSS-FSS. This indicates that the effects of the mechanical boundary on dispersion curves are more pronounced than that of electrical and magnetical boundaries.

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(d)**

Next, the effects of varying stacking sequence to plate on dispersion curves are examined. The boundary conditions of the top and bottom surfaces are specified to FOO-FOO. The plates with stacking sequences B/B/B, B/F/B, F/B/F, and F/F/F are considered, respectively. The dispersion curves for the first four modes are shown in Figures 4(a)–4(d), respectively, where the solid, dotted, dash-dotted, and dashed curves denotes the dispersion curves for stacking sequences B/B/B, B/F/B, F/B/F, and F/F/F, respectively. Figure 4(a) shows that for the first mode, the dispersion curves for the four different stacking sequences are close to each other. However, from Figures 4(b)–4(d), it is clear that for the second, third, and fourth modes, the dispersion curves are quite different for the four different stacking sequences.

To further demonstrate the performance of the proposed methods, the results obtained from the proposed methods are compared with those obtained from the semianalytical finite-element (SAFE) method [46]. Using the linear element along the *z*-axis direction, the semidiscretized equation can be obtained byin which stiffness matrix **K** and mass matrix **M** can be obtained by assembling the element counterparts. Then, equation (29) can be solved by using the *eigs* function in MATLAB for a given wavenumber . At , the dimensionless frequencies obtained from this method and the SAFE method discretized each layer with 0.001 m thickness sublayer are listed in Table 2. The results show very good agreement between the results obtained from the proposed method and the SAFE method.

##### 6.2. A Multilayered Strongly Anisotropic MEE Structure

The aforementioned example focuses on the transversely isotropic material. To confirm the feasibility of this method for the application of more generally anisotropic problem, a strongly anisotropic multilayered MEE structure is further considered in this section. The structure consists of four layers with stacking sequences of materials BaTiO_{3}, CoFe_{2}O_{4}, Bimaterial A, and Bimaterial B. Dispersion curves of the first 10 modes for are shown in Figure 5, where (a), (b), (c), and (d) denote the results for the boundary conditions of FOO-FOO, COO-COO, CSS-FOO, and CSS-FSS, respectively. Comparing Figure 5(a) with Figure 5(b), we observe that the dispersion curves are significantly different when the boundary condition changes from FOO-FOO to COO-COO. However, comparing Figure 5(c) with Figure 5(d), we observe that the dispersion curves have only slight variations when the boundary condition changes from CSS-FOO to CSS-FSS. Therefore, this example obtains similar conclusion to Section 6.1, that is, the effects of mechanical boundary are more significant than that of electrical and magnetical on dispersion curves.

**(a)**

**(b)**

**(c)**

**(d)**

For the dispersion diagram shown in Figure 5(a), at , the first three eigenfrequencies are rad·*μ*s^{−1}, rad·*μ*s^{−1} and rad·*μ*s^{−1}, and their corresponding eigenfunctions are shown in Figures 6 and 7. It can be seen from Figures 6 and 7 that all of mechanical, electrical, and magnetical components are complex functions, and their real and imaginary parts are separately plotted in (a)–(e) and in (f)–(g). From Figures 6 and 7, we observe that all these components *u*, , , , and are not identically zero, which indicates that the wave motion is coupled with regard to all components. Furthermore, at 0.1, 1, 5, and 10, the eigenfrequencies in the first five modes for four types of boundary conditions are listed in Table 3, which are provided for examination and reference by other researchers. In addition, to confirm the accuracy of this method, at 0.1, 1, 5, and 10, the first five frequencies of this method and the SAFE method discretized each layer with 0.001 m thickness sublayer for four types of boundary conditions are listed in Table 3.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

**(i)**

**(j)**

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

**(i)**

**(j)**

##### 6.3. A Periodically Multilayered Anisotropic MEE Structure

In this example, the results for both periodically multilayered finite and infinite MEE structures obtained from the presented method are presented. The unit cell consists of four layers with materials BaTiO_{3}, CoFe_{2}O_{4}, Bimaterial A and Bimaterial B, and their thicknesses are the same. The incident angle is assumed to be .

The dispersion curves of the first 10 modes for the periodically finite MEE structure for different types of boundary conditions and different number of unit cells are shown in Figure 8, where the dashed, solid, and dotted curves denote the phase velocities for the periodic multilayered plate composed of 1, 10, and 100 unit cells, respectively. The subfigures (a)–(d) represent the dispersion curves for the boundary conditions of FOO-FOO, COO-COO, CSS-FOO, and CSS-FSS. From Figure 8, we observe that the phase velocities for the structure consisting of 1 unit cell are significantly different from that consisting of 10 unit cells. However, the phase velocities for the structure consisting of 10 and 100 unit cells have small changes.

**(a)**

**(b)**

**(c)**

**(d)**

Finally, we investigate a periodic multilayered structure arranged by an infinite sequence of unit cell, and the unit cell is identified to that given in the finite periodic structure. Wave propagation in a periodical structure can exhibit a characteristic feature of band structure, known as passbands or Bragg bandgaps, within which the waves can propagate to infinity. In complementary frequency bands, known as stopbands or locally resonant bandgaps, the waves are effectively attenuated. Figure 9 shows the dispersion curves and band structures of the periodically infinite multilayered structure in the first 10 modes, where (a), (b), (c), and (d) denote the results for 0.1, 1, 5 and 10, respectively. The bandgaps are closely related to . When , the first bandgap starts from zero frequency. Only when , the first bandgap does not start from zero frequency. From Figure 9, it can be observed that as increases, the width of the first stopband becomes large and the number of stopbands exhibits an increasing trend.

**(a)**

**(b)**

**(c)**

**(d)**

#### 7. Conclusions

By combining the PIM with the W-W algorithm, based on the symplectic structure of the Hamiltonian matrix, a generalized method for calculating dispersion curves of waves in multilayered anisotropic MEE structures was proposed. To facilitate the application of the W-W algorithm in MEE structures, based on the symplectic structure of the Hamiltonian matrix, a symplectic transformation is introduced. The feasibility of the W-W algorithm after performing the symplectic transformation is demonstrated by a strictly theoretical analysis. The method combines the advantages of the PIM with those of the W-W algorithm so that the computation is accurate and stable, and no eigenfrequency is missed. Several typical examples are presented to verify the numerical performance of this method. In examples, the dispersion curves of waves in multilayered anisotropic MEE structures are calculated by using this method. It was shown that this method is highly accurate by comparing with the SAFE method.

The proposed method can provide the basis of theory for the application of the PIM and the W-W algorithm to solve the wave propagation and dynamic problems of MEE structure. The method is not restricted to the specific material properties, the layer number of unit cells, and the complexity of the unit cell. Based on the idea of the proposed method, it can be easily applied to more complex problems. For example, by combining with the finite-element method, the proposed method can be extended to solve dispersions or bandgaps of finite, infinite, and semi-infinite anisotropic MEE phononic crystals.

#### Appendix

#### A. The Physical Quantities and Material Constants For Elastic, Electric, and Magnetic Fields

The physical quantities for elastic, electric, and magnetic fields are

The coefficient matrices for elastic, electric and magnetic fields are

For an MEE solid, the coupled potential energy is

The total internal energy function *U* can be given by

in which

To ensure the stability of the system, it requires . Therefore, the elastic coefficient matrix **c** and magneto-electric coupled matrix are both positive definite.

#### B. Proof for the Positive Definitiveness of D

For the convenience of the proof, a 5×5 matrix **X** representing , , , , , , or , etc. is partitioned with regard to mechanics and coupled electricity magnetism parts, which are represented by the superscripts *M* and *E*, that is,where the sizes of the block matrices , , , and are 3×3, 3×2, 2×3, and 2×2, respectively.

Using equations (10), (14), (15), (16), and (B.1), the matrix **D** can be given bywhere

The matrix **D** defined in equation (B.2) can also be written in the following form:

Obviously, the positive definitiveness of **D** is equivalent to that and are both positive definite. According to equation (B.1) and the expression of given in equation (7), and are positive and negative definite, then from equation (B.3), we can obtain that is positive definite. Therefore, to prove the positive definitiveness of **D**, we only need to prove that is positive definite.

Let

Substituting and given in equation (B.3) into equation (B.5) and simplifying the obtained results gives

Then, by using the relation of defined in equation (11) and the inverse of block matrix [47], we have

So, equation (B.6) can be further expressed as

According to equation (7)–(9) and (B.1), we have

with

Substituting equation (B.9) into (B.8) gives

Defining a 6 × 6 matrix where the sizes of submatrices , , and and are 4 × 4, 4 × 2, 2 × 4, and 2 × 2, respectively, and using the inverse of block matrix [47], from equation (B.12) we have

And thus, equation (B.11) can be written as

By exchanging columns and rows of the matrix defined in equation (B.12), can be changed as the magneto-electric coupled matrix , that is,

In Appendix A, it has been demonstrated that is positive definite, so is also positive definite, such that defined in equation (B.13) is positive definite. Furthermore, from equation (B.11) we conclude that the matrix defined in equation (B.11) is positive definite. Finally, by combining with the positive definiteness of , it can conclude that **D** is positive definite.

#### Data Availability

The data used to support the findings of this study are included within the article.

#### Conflicts of Interest

The authors declare no conflicts of interest.

#### Acknowledgments

The authors are grateful for the support of the National Natural Science Foundation of China (Grant nos. 11972107 and 91748203) and the Fundamental Research Funds for the Central Universities (Grant no. DUT2019TD37).