Multiscale Finite Element Method for Coupling Analysis of Heterogeneous Magneto-Electro-Elastic Structures in Thermal Environment

Xinyue Li, Xiaolin Li^*, Hangran Yang

College of Construction Engineering, Jilin University, Changchun, China.

DOI: 10.4236/jamp.2024.129186                

Abstract

Magneto-electro-elastic (MEE) materials, a new type of composite intelligent materials, exhibit excellent multifield coupling effects. Due to the heterogeneity of the materials, it is challenging to use the traditional finite element method (FEM) for mechanical analysis. Additionally, the MEE materials are often in a complex service environment, especially under the influence of the thermal field with thermoelectric and thermomagnetic effects, which affect its mechanical properties. Therefore, this paper proposes the efficient multiscale computational method for the multifield coupling problem of heterogeneous MEE structures under the thermal environment. The method constructs a multi-physics field with numerical base functions (the displacement, electric potential, and magnetic potential multiscale base functions). It equates a single cell of heterogeneous MEE materials to a macroscopic unit and supplements the macroscopic model with a microscopic model. This allows the problem to be solved directly on a macroscopic scale. Finally, the numerical simulation results demonstrate that compared with the traditional FEM, the multiscale finite element method (MsFEM) can achieve the purpose of ensuring accuracy and reducing the degree of freedom, and significantly improving the calculation efficiency.

Keywords

Multiscale Finite Element Method, Magneto-Electro-Elastic, Multifield Coupling, Numerical Base Functions

1. Introduction

Magneto-electro-elastic (MEE) materials can convert mechanical energy, electrical energy, and magnetic energy into each other [1], and their electromagnetic coupling effect comes from the interaction between the piezoelectric phase and the piezomagnetic phase. The mechanism is as follows: after the piezoelectric phase under the action of the electric field produces mechanical deformation due to the inverse piezoelectric effect, the interface transmits this deformation to the piezomagnetic phase, and then the magnetic field is induced by the piezomagnetic effect in the piezomagnetic phase. The process can also occur in reverse [2]. The unique multifield coupling properties of MEE materials make them suitable for a wide range of applications in sensors, transducers, energy harvesters, etc. [3] [4]. As the application field of MEE materials expands, their service environment becomes increasingly complex and is significantly influenced by temperature changes in the environment. MEE materials have thermoelectric and thermomagnetic effects that produce electromagnetic polarization and deformation in the presence of a thermal field. Therefore, accurately solving the multi-physical coupling response of MEE structures under the influence of a thermal field is of great value for enhancing the performance of intelligent structural components.

To explore more deeply the intrinsic mechanism of multifield coupling in MEE materials under the action of a thermal field, researchers have worked on developing analytical solutions. Based on the theory of high-order shear deformation, Hamilton’s principle, and the constitutive equation of the thermo-mechanical-electro-magneto coupling effect, Zhang et al. [5] studied the buckling behavior of the functionally graded magneto-electro-elastic (FGMEE) cylindrical shell. Arefi et al. [6] proposed a simplified shear theory and non-local theory of normal deformation for flexural analysis of FGMEE nanobeams. However, the analytical solutions are generally only suitable for solving problems with simple structural geometry and boundary conditions, and they are difficult to apply to complex practical problems. To compensate for this limitation, a large number of numerical methods have been developed rapidly, and have gradually become the most effective means to analyze MEE structures.

Among them, the finite element method (FEM) is widely used because of its high efficiency and simplicity. Farai et al. [7] used the finite element method to analyze civil engineering structures, which makes the analysis of complex structural systems effective and accurate. Du et al. [8] performed numerical simulations using finite element analysis techniques to simulate and calculate the stress structure in the seismically excited near-source zone by establishing a numerical model of the explosion in the well. Yang et al. [9] investigated the response of the FGMEE structures under different loads by using the FEM and calculated the stresses, electric displacement, and magnetic induction intensity of the structures. Vinyas et al. [10]-[12] used FEM to solve the static response and natural frequency of MEE plates and analyzed the influence of the thermal field on structural response. Zenkour et al. [13] [14] investigated the mechanical characteristics of FGMEE hollow cylinders in terms of radial displacement, stress, and temperature. However, the FEM for solving composites with complex microstructures such as heterogeneous MEE structures under thermal field has an oversized matrix construction, which significantly increases the degrees of freedom of the finite element model, is computationally inefficient, and requires a large amount of computational resources. The complexity of the multi-physical field coupling relationship is difficult to describe, which further increases the difficulty of finite element analysis, making the problem more difficult to solve. It is necessary to develop the multiscale numerical calculation method to solve the multifield coupling problem in the thermal field more efficiently.

Related work on the multiscale finite element method (MsFEM) dates back to Babuska et al. and Hou et al. [15] [16]. Its primary concept revolves around constructing multiscale basis functions that enable the extraction of small-scale information within each coarse element through the localized resolution of boundary value problems. Yang et al. [17] developed an extended multiscale finite element (EMsFEM) based on the MsFEM to solve the thermo-elastic-plastic problem of heterogeneous materials. Xing et al. [18] used MsFEM to analyze periodic composite structures and compared the adaptability and similarity with the mathematical homogenization method and heterogeneous multiscale method. Fu et al. [19] used the MsFEM for geometrically nonlinear analysis of heterogeneous piezoelectric materials and demonstrated the effectiveness and efficiency of the method through numerical examples.

As mentioned above, the MsFEM can solve the mechanical problems of heterogeneous materials very well, but it has not, to our knowledge, been introduced into the analysis of MEE materials in thermal field environments. For this reason, in this paper, based on the basic principles of MsFEM, the heterogeneous MEE problem is solved on a macroscopic scale under the influence of the temperature field by capturing the heterogeneity in mechanics, electricity, and magnetic at the same time. Finally, the correctness and efficiency of the algorithm are verified by numerical examples.

The structure of this paper is as follows. Section 2 introduces the basic equation of MEE structure under the action of the thermal field. Section 3 introduces the calculation principle based on the MsFEM from microscopic calculation and macroscopic calculation. Section 4 is numerical verification, which proves the accuracy and efficiency of the multiscale method by numerical examples. The last part summarizes the work of this paper.

2. The Basic Formulations

The generalized equilibrium equations for MEE materials without considering body force, volume charge density, and volume current density:

$L_{\text{d}}^{\text{T}}\sigma \left(x\right)=0$ ,(1)

$L_{\text{b}}^{\text{T}}D\left(x\right)=0$ ,(2)

$L_{\text{b}}^{\text{T}}B\left(x\right)=0$ ,(3)

where $\sigma$ is the stress vector, $D$ is the electrical displacement vector, and $B$ is the magnetic field intensity vector. $L_{\text{d}}^{\text{T}}$ and $L_{\text{b}}^{\text{T}}$ are matrices of differential operators, whose expressions are as follows:

$L_{\text{d}}=\left[\begin{array}{cc}\frac{\partial }{\partial x}& 0\\ 0& \frac{\partial }{\partial y}\\ \frac{\partial }{\partial y}& \frac{\partial }{\partial x}\end{array}\right],L_{\text{b}}=\left[\begin{array}{c}\frac{\partial }{\partial x}\\ \frac{\partial }{\partial y}\end{array}\right]$ .(4)

The generalized geometric equations of the MEE materials are given as:

$S\left(x\right)=\frac{1}{2}{L}_{\text{d}}u\left(x\right)$ ,(5)

$E\left(x\right)=-L_{\text{b}}\varphi \left(x\right)$ ,(6)

$H\left(x\right)=-L_{\text{b}}\psi \left(x\right)$ ,(7)

where $S$ , $E$ , $H$ are the vectors of strain, electrical filed intensity, and magnetic flux density, respectively. $u$ , $\varphi$ , $\psi$ represent the vectors of displacement, electric potential, and magnetic potential, respectively.

The constitutive equations of MEE materials are given as:

$\sigma \left(x\right)=CS\left(x\right)-eE\left(x\right)-hH\left(x\right)-\beta \Delta T$ ,(8)

$D\left(x\right)=eS\left(x\right)+\kappa E\left(x\right)+\alpha H\left(x\right)+p\Delta T$ ,(9)

$B\left(x\right)=hS\left(x\right)+\alpha E\left(x\right)+\mu H\left(x\right)+\tau \Delta T$ ,(10)

where $C$ , $e$ and $h$ are the matrices of elastic coefficient, piezoelectric coefficient, and piezomagnetic coefficient, respectively. $\beta$ is the matrix of the thermal expansion coefficient. $\kappa$ , $\alpha$ , and $\mu$ denote dielectric coefficient matrix, magnetoelectric coefficient matrix, and magnetic permittivity matrix, respectively. $p$ and $\tau$ are the matrices of pyroelectric and pyromagnetic coefficients, respectively. The expressions are as follows:

$\begin{array}{l}C=\left[\begin{array}{ccc}{C}_{11}& C_{13}& 0\\ C_{13}& C_{33}& 0\\ 0& 0& C_{44}\end{array}\right],e=\left[\begin{array}{cc}0& e_{31}\\ 0& e_{33}\\ e_{15}& 0\end{array}\right],h=\left[\begin{array}{cc}0& q_{31}\\ 0& q_{33}\\ q_{15}& 0\end{array}\right],\\ \kappa =\left[\begin{array}{cc}{\kappa }_{11}& 0\\ 0& {\kappa }_{33}\end{array}\right],\alpha =\left[\begin{array}{cc}{\alpha }_{11}& 0\\ 0& {\alpha }_{33}\end{array}\right],\mu =\left[\begin{array}{cc}{\mu }_{11}& 0\\ 0& {\mu }_{33}\end{array}\right],\\ \beta =\left[\begin{array}{ccc}{C}_{11}& C_{13}& 0\\ C_{13}& C_{33}& 0\\ 0& 0& C_{44}\end{array}\right]\left[\begin{array}{c}{\beta }_1\\ {\beta }_2\\ 0\end{array}\right],p=\left[\begin{array}{c}{p}_1\\ p_3\end{array}\right],\tau =\left[\begin{array}{c}{\tau }_1\\ {\tau }_3\end{array}\right].\end{array}$ (11)

3. The Basic Principle of MsFEM

The base function of the FEM is generally in analytical form, which is only related to the coordinates of the unit nodes and has nothing to do with the material parameters, so it cannot take into account the variability of the material properties within the unit. The MsFEM adopts a double-mesh system. Following the discretization of the structure into coarse-grid elements, the coarse grid domain undergoes further refinement into a sub-grid system that is capable of capturing material heterogeneity in detail. Then the multiscale base function is constructed by solving the equilibrium equation under the specific boundary conditions of the coarse mesh so that the solution on the coarse-grid complies with the required precision specifications. Thereby, it solves the problem of large computational volume and long solution time when solving heterogeneous materials by the FEM.

The multi-physical numerical base functions (displacement, electric potential, and magnetic potential multiscale base functions) constructed by MsFEM embed microscale heterogeneity into it, taking into account the complex structure and property changes inside the material at the microscale. It is reflected in the simulation results of macroscopic scale through the assembly of macroscopic matrix, so that the original problem can be solved directly on the macroscopic scale. To achieve this goal, the method constructs a multiscale model by establishing an association between a macroscopic model, which focuses on the overall mechanical properties of the material, and a microscopic model, which details the internal microstructure and properties of the material. This multiscale model construction method can simultaneously consider the behavior of materials at different scales.

3.1. Construction of Numerical Base Functions

When using the MsFEM to analyze the mechanical-thermal-electro-magnetic coupling problem of MEE materials, the whole solution domain is divided into coarse meshes first, and the fine meshes are further divided within each coarse meshes so that the fine meshes can fully analyze the heterogeneity of materials. Appropriate boundary conditions are then defined for the local problem on the fine grid, and the solution of the local problem is used to construct the numerical base function.

Under the elastic small deformation hypothesis, three kinds of multiscale base functions which can reflect the internal physical characteristics of the element are constructed: displacement base function ${\overline{N}}^u$ and ${\overline{N}}^v$ , electric potential base function ${\overline{N}}^{\varphi }$ , and magnetic potential base function ${\overline{N}}^{\psi }$ . The y-direction displacement base function is constructed similarly to the x-direction, and the magnetic potential base function is constructed similarly to the electric potential. Thus, this paper only describes the construction process of ${\overline{N}}^u$ and ${\overline{N}}^{\varphi }$ .

1) Displacement multiscale base functions

The displacement field is a vector field and the deformations in all directions interact with each other. According to the basic principles of MsFEM, considering the coupling effect between the displacement field in each coordinate direction when the solid deformation introduces the coupling additional term, then the displacement of any nodes of the fine mesh within the coarse mesh can be expressed as:

${\left[u_{11}{u}_{21}\cdots u_{2\tilde{n}}\right]}^{\text{T}}={\left[R_1^{u\mathrm{T}}{R}_2^{u\mathrm{T}}\cdots R_{\tilde{n}}^{u\mathrm{T}}\right]}^{\text{T}}{\left[{\overline{u}}_{11}{\overline{u}}_{21}\cdots {\overline{u}}_{2M}\right]}^{\text{T}}$,(12)

$R_i^u=\left[\begin{array}{ccccc}{\overline{N}}_{111}^u\left(i^{\prime }\right)& {\overline{N}}_{121}^u\left(i^{\prime }\right)& \cdot \cdot \cdot & {\overline{N}}_{11M}^u\left(i^{\prime }\right)& {\overline{N}}_{12M}^u\left(i^{\prime }\right)\\ {\overline{N}}_{211}^u\left(i^{\prime }\right)& {\overline{N}}_{221}^u\left(i^{\prime }\right)& \cdot \cdot \cdot & {\overline{N}}_{21M}^u\left(i^{\prime }\right)& {\overline{N}}_{22M}^u\left(i^{\prime }\right)\end{array}\right]$ ,(13)

where $u_{ij}$ denotes the fine-grid nodes displacement. ${\overline{u}}_{ij}$ represents the coarse-grids node displacement. ${\overline{N}}_{ijk}^u$ denotes the displacement of coarse grid nodes in the j direction when unit displacement is applied to the coarse-grid node k along the i direction. M represents the number of coarse-grid nodes. $\tilde{n}$ is a coarse grid contains the number of fine-grid nodes; $i^{\prime }=1,2,\cdots ,\tilde{n}$ .

The displacement multiscale base function must be satisfied:

$\{\begin{array}{l}L{\overline{N}}_{\dot{k}}^u=0,\text{on}{\Omega }_c\\ {\overline{N}}_{\dot{k}}^u={\tilde{N}}_{\dot{k}}^u,\text{on}\partial {\Omega }_c\end{array}$,(14)

where the differential operator L satisfies:

$Lu=\text{div}\left(C:\frac{1}{2}\left(\nabla u+{\left(\nabla u\right)}^{\text{T}}\right)\right)$ ,(15)

and ${\Omega }_c$ represents the region of coarse-grids. $\partial {\Omega }_c$ represents the boundary of coarse-grids. ${\overline{N}}_{\dot{k}}^u$ represents the displacement numerical base function on a coarse grid $\dot{k}$ . ${\tilde{N}}_{\dot{k}}^u$ represents the boundary condition of base functions. $\dot{k}=1,2,\cdots ,M$ .

The base function is numerically generated by computing the MEE material equations within a single-cell, where specific boundary conditions are applied. In this paper, linear boundary conditions and periodic boundary conditions are used to construct the base function.

Figure 1 demonstrates the use of linear boundary conditions in constructing the displacement base function. A forward unit displacement of x is applied at node 1. The displacement in the x-direction of the boundary ${\Gamma }_{23}$ and ${\Gamma }_{34}$ of the coarse-grids is fixed, and the displacement in the y-direction of all the boundaries of the coarse-grids is fixed. The x-direction displacements on the boundaries ${\Gamma }_{12}$ and ${\Gamma }_{14}$ change linearly from one to zero.

Figure 1. Using linear boundary conditions for constructing the displacement base function.

For materials characterized by periodic microscopic structures, the displacement multiscale base functions can be constructed using periodic boundary conditions, which are different from the forced constraints imposed by linear boundary conditions, as depicted in Figure 2. Periodic boundary conditions impose dynamic constraints to ensure more natural deformation of the sub-grids. The displacement increment $\Delta u$ is the specified constant. The value is assigned unity at node 1 and zeroed out at nodes 2 and 4. Then it varies linearly along the boundaries ${\Gamma }_{12}$ and ${\Gamma }_{14}$ . Meanwhile, the degrees of freedom at the macroscopic node 3 are fixed in both x- and y-directions to prevent rigid body motion.

Figure 2. Using periodic boundary conditions for constructing the displacement base function.

2) Electric potential multiscale base functions

The electric potential is a scalar quantity, characterized solely by its magnitude and devoid of direction. The electric potential of any nodes of the fine-grid is expressed as:

${\left[{\psi }_1{\psi }_2\cdots {\psi }_{\tilde{n}}\right]}^{\text{T}}={\left[R_1^{\psi \text{T}}{R}_2^{\psi \text{T}}\cdots R_{\tilde{n}}^{\psi \text{T}}\right]}^{\text{T}}{\left[{\overline{\psi }}_1{\overline{\psi }}_2\cdots {\overline{\psi }}_M\right]}^{\text{T}}$,(16)

$R_{\tilde{i}}^{\psi }=\left[{\overline{N}}_1^{\psi }\left(\tilde{i}\right){\overline{N}}_2^{\psi }\left(\tilde{i}\right)\cdots {\overline{N}}_M^{\psi }\left(\tilde{i}\right)\right],\tilde{i}=1,2,\cdots ,\tilde{n}$,(17)

where ${\psi }_{ij}$ represents the electric potential of fine-grid nodes. ${\overline{\psi }}_{ij}$ represents the electric potential of coarse-grid nodes.

The electric potential multiscale base function must be satisfied:

$\{\begin{array}{l}\overline{L}{\overline{N}}_{\tilde{k}}^{\psi }=0,\text{on}{\Omega }_c\\ {\overline{N}}_{\tilde{k}}^{\psi }={\tilde{N}}_{\tilde{k}}^{\psi },\text{on}\partial {\Omega }_c\end{array}$,(18)

where the differential operator $\overline{L}$ satisfies:

$\overline{L}\varphi =-{\nabla }^{\text{T}}\left(\kappa \nabla \varphi \right)$,(19)

and ${\overline{N}}_{\tilde{k}}^{\varphi }$ represents the numerical base function of electric potential on coarse-grid $\tilde{k}$ ; ${\tilde{N}}_{\tilde{k}}^{\varphi }$ is the boundary conditions of the numerical base functions; $\tilde{k}=1,2,\cdots ,M$ .

Figure 3 demonstrates the use of linear boundary conditions in constructing the electric potential base function. The unit electric potential is applied at macroscopic node 1 of the coarse grid. The electric potential on the boundary is grounded, and the electric potential value on the boundary ${\Gamma }_{23}$ , ${\Gamma }_{34}$ changes linearly from one to zero along the direction 14 and 12, respectively.

Figure 3. Using linear boundary conditions for constructing the electric potential base function.

Figure 4 demonstrates the use of periodic boundary conditions in constructing the electric potential base function. The degrees of freedom in the x and y directions of node 3 are fixed and the electric potential is grounded. To ensure the periodicity of electrical characteristics of adjacent coarse-grid elements, dynamic constraints are imposed on the nodes of two sets of opposite boundaries of coarse-grid elements. The electric potential increment $\Delta \varphi$ is the specified constant. The value is assigned unity at node 1 and zeroed out at nodes 2 and 4. Then it varies linearly along the boundaries ${\Gamma }_{12}$ and ${\Gamma }_{14}$ .

Figure 4. Using periodic boundary conditions for constructing the electric potential base function.

3.2. Macroscopic Computation

Based on Green’s formulation and Galerkin method using Equations (12) and (16), the system equation of MEE structure in the thermal field can be obtained as follows:

$K_{uu}u+K_{u\varphi }\varphi +K_{u\psi }\psi =F_u+F_T$ ,(20)

$K_{\varphi u}u-K_{\varphi \varphi }\varphi -K_{\psi \psi }\psi =-F_{pe}$ ,(21)

$K_{\psi u}u-K_{\psi \varphi }\varphi -K_{\psi \psi }\psi =-F_{pm}$ ,(22)

where $F_u$ is the static force applied to the structure. $F_T$ , $F_{pe}$ , $F_{pm}$ are temperature load, thermoelectric load, and thermomagnetic load, respectively. The expressions are as follows:

$\begin{array}{l}{F}_u=Q_e^{u\text{T}}{f}^e,F_T=Q_e^{u\text{T}}{f}_T^e\\ F_{pe}=Q_e^{\varphi \text{T}}{f}_{pe}^e,F_{pm}=Q_e^{\psi \text{T}}{f}_{pm}^e,\end{array}$ (23)

and $K_{uu}$ , $K_{u\varphi }$ , $K_{\varphi u}$ , $K_{\varphi \varphi }$ , $K_{u\psi }$ , $K_{\psi u}$ , $K_{\psi \varphi }$ , and $K_{\psi \psi }$ are the generalized stiffness matrices of coarse-grid elements, and the expressions are:

$\begin{array}{l}{K}_{uu}=Q_e^{u\text{T}}{K}_{uu}^e{Q}_e^u,K_{u\varphi }=Q_e^{u\text{T}}{K}_{u\varphi }^e{Q}_e^{\varphi },K_{\varphi u}=Q_e^{\varphi \text{T}}{K}_{\varphi u}^e{Q}_e^u,\\ K_{\varphi \varphi }=Q_e^{\varphi \text{T}}{K}_{\varphi \varphi }^e{Q}_e^{\varphi },K_{u\psi }=Q_e^{u\text{T}}{K}_{u\psi }^e{Q}_e^{\psi },K_{\psi u}=Q_e^{\psi \text{T}}{K}_{\psi u}^e{Q}_e^u,\\ K_{\varphi \psi }=Q_e^{\varphi \text{T}}{K}_{\varphi \psi }^e{Q}_e^{\psi },K_{\psi \varphi }=Q_e^{\psi \text{T}}{K}_{\psi \varphi }^e{Q}_e^{\varphi },K_{\psi \psi }=Q_e^{\psi \text{T}}{K}_{\psi \psi }^e{Q}_e^{\psi }.\end{array}$ (24)

In the Equation (23), $f^e$ , $f_T^e$ , $f_{pe}^e$ , $f_{pm}^e$ are the fine-grid unit load arrays, and the expressions are:

$\begin{array}{l}{f}^e={\displaystyle {\int }_{{\Omega }^e}{N}^{u\text{T}}f\text{d}\Omega }+{\displaystyle {\int }_{{\Gamma }_{\sigma }}{N}^{u\text{T}}\overline{t}\text{d}\Gamma },f_T^e={\displaystyle \sum _{i=1}^{N_n}{B}_u^{\text{T}}C\beta \Delta T{A}_i},\\ f_{pe}^e={\displaystyle \sum _{i=1}^{N_n}{B}_{\varphi }^{\text{T}}p\Delta T{A}_i},f_{pm}^e={\displaystyle \sum _{i=1}^{N_n}{B}_{\psi }^{\text{T}}\tau \Delta T{A}_i},\end{array}$(25)

where i represents the number of fine-grid elements. A_i represents the area of fine-grid elements. $\Delta T$ is the amount of temperature change. $B_u$ , $B_{\varphi }$ , $B_{\psi }$ represents the strain transfer matrices. And $Q_e^u$ , $Q_e^{\varphi }$ , $Q_e^{\psi }$ serve as transformation matrices, facilitating the interconversion between the displacement, electric potential, and magnetic potential values at the nodes of coarse and fine grids. The expressions are as follows:

$Q_e^u={\left[R_{e1}^{u\text{T}}{R}_{e2}^{u\text{T}}\cdots R_{e\overline{n}}^{u\text{T}}\right]}^{\text{T}}$,(26)

$Q_e^{\varphi }={\left[R_{e1}^{\varphi \text{T}}{R}_{e2}^{\phi \text{T}}\cdots R_{e\overline{n}}^{\varphi \text{T}}\right]}^{\text{T}}$,(27)

$Q_e^{\psi }={\left[R_{e1}^{\psi \text{T}}{R}_{e2}^{\psi \text{T}}\cdots R_{e\overline{n}}^{\psi \text{T}}\right]}^{\text{T}}$.(28)

In the Equation (24), $K_{uu}^e$ , $K_{u\varphi }^e$ , $K_{\varphi u}^e$ , $K_{\varphi \varphi }^e$ , $K_{u\psi }^e$ , $K_{\psi u}^e$ , $K_{\psi \varphi }^e$ , and $K_{\psi \psi }^e$ are the stiffness matrices of fine-grid elements, and the expression are:

$\begin{array}{l}{K}_{uu}^e={\displaystyle {\int }_{{\Omega }^e}{B}_u^{\text{T}}C{B}_u\text{d}\Omega },{{\rm K}}_{u\varphi }^e={\displaystyle {\int }_{{\Omega }^e}{B}_u^{\text{T}}e{B}_{\varphi }\text{d}\Omega },\\ K_{\varphi u}^e={\displaystyle {\int }_{{\Omega }^e}{B}_{\varphi }^{\text{T}}e{B}_u\text{d}\Omega },K_{\varphi \varphi }^e=-{\displaystyle {\int }_{{\Omega }^e}{B}_{\varphi }^{\text{T}}\kappa B_{\varphi }\text{d}\Omega },\\ K_{u\psi }^e={\displaystyle {\int }_{{\Omega }^e}{B}_u^{\text{T}}h{B}_{\psi }\text{d}\Omega },K_{\psi u}^e={\displaystyle {\int }_{{\Omega }^e}{B}_{\psi }^{\text{T}}h{B}_u\text{d}\Omega },\\ K_{\varphi \psi }^e=-{\displaystyle {\int }_{{\Omega }^e}{B}_{\varphi }^{\text{T}}\alpha B_{\psi }\text{d}\Omega },K_{\psi \varphi }^e=-{\displaystyle {\int }_{{\Omega }^e}{B}_{\psi }^{\text{T}}\alpha B_{\varphi }\text{d}\Omega },\\ K_{\psi \psi }^e=-{\displaystyle {\int }_{{\Omega }^e}{B}_{\psi }^{\text{T}}\mu B_{\psi }\text{d}\Omega },\end{array}$ (29)

where ${\Omega }^e$ is the corresponding unit area.

4. Numerical Examples

In this section, two numerical examples are used to compare the solution obtained by the MsFEM with the solution obtained by the standard FEM at the fine scale to verify the effectiveness and high efficiency of the MsFEM. The materials related parameters adopted by all structures are as follows:

(Elastic~C_ij~N∙m⁻², Dielectric~t_ij~C∙V∙m⁻¹, Magnetic~d_ij~N∙s²∙C⁻², Piezoelectric~e_ij~C∙m⁻², Magnetoelectric~m_ij~N∙s∙V∙C⁻¹, Piezomagnetic~h_ij~N∙A∙m⁻¹, Thermal expansion~β_i~K⁻¹, Pyroelectric~p_i~C∙m⁻²∙K⁻¹, Pyromagnetic~τ_i~N∙A⁻¹∙m⁻¹∙K⁻¹, Density~ρ~kg∙m⁻³)

$\begin{array}{l}{C}_{11}=166.0\times {10}^9,C_{12}=77.0\times {10}^9,C_{13}=78.0\times {10}^9,C_{33}=162.0\times {10}^9,\\ C_{44}=43.0\times {10}^9,t_{11}=11.2\times {10}^{-9},t_{33}=12.6\times {10}^{-9},d_{11}=0.05\times {10}^{-4},\\ d_{33}=0.1\times {10}^{-4},e_{31}=-4.4,e_{33}=18.6,e_{15}=11.6,m_{11}=5.0\times {10}^{-12},\\ m_{33}=3.0\times {10}^{-12},h_{31}=580.30,h_{33}=699.70,h_{15}=550.0,\\ {\beta }_1={\beta }_2=10.8\times {10}^{-6},p_3=-3.5\times {10}^{-7},{\tau }_3=-36\times {10}^{-5},\rho =\mathrm{5730.0.}\end{array}$

In the calculation results, “FEM” represents the solution of standard FEM on a fine-grid, “MsFEM-L” represents the MsFEM solution based on linear boundary conditions, and “MsFEM-P” represents the MsFEM solution based on periodic boundary conditions.

Example 1: A FGMEE clamped-clamped (C-C) beam under the temperature variation $\Delta T=20\text{K}$ is shown in Figure 5. The geometric size of the beam is L = 30 mm and the width H = 2 mm, which is a plane stress problem. The constraints at both clamped sides are $u_x=u_y=0$ and the constraints on the bottom boundary AD of the structure are $\varphi =\psi =0$ . The structural material is BaTio₃-CoFe₂O₄, and the properties of the material vary gradient along the y-direction and follow the exponential distribution $f\left(y\right)={\text{e}}^{\xi \left(y/h\right)}$ Exponential factor is $\xi =0.2$ .

Figure 5. Geometry and grid division of the FGMEE C-C beam.

In multiscale calculations, the entire model is divided into 90 × 6 coarse-grids, and each coarse-grid is divided into 5 × 5 fine-grids. The standard FEM (180 × 12) and the MsFEM (90 × 6) are used to solve the multi-physical coupling problem of the FGMEE C-C beam of the heterogeneous mechanical-thermo-electro-magneto. The generalized displacements (displacements u_x and u_y, electric potential $\varphi$ , and magnetic potential $\psi$ ) of upper boundary AD are obtained. Then, the solution results of standard FEM and MsFEM were compared with the solution in the literature [20], as shown in Figure 6. It can be seen that the results of the three algorithms are consistent with the analytical solutions in the literature, which verifies the correctness of the three algorithms and programs and the rationality of using the FEM results as the reference solutions. At the same time, the maximum errors of displacements u_x and u_y, electric potential $\varphi$ , and magnetic potential are 3.71%, 2.54%, 0.92%, and 3.44%, respectively, which shows that this method has high calculation accuracy.

(a) (b)

(c) (d)

Figure 6. Generalized displacements at edge AD of the FGMEE C-C beam. (a) displacement u_x; (b) displacement u_y; (c) electrical potential ϕ; (d) magnetic potential ψ.

Under identical hardware conditions, the calculation time of the three algorithms is compared, as shown in Table 1. MsFEM-L and MsFEM-P take only 42.54% and 42% of the time of FEM. Therefore, MsFEM significantly improves the efficiency of calculation, and this advantage will be more obvious when dealing with large degrees of freedom structures.

Table 1. Comparing the duration of calculations for the FGMEE C-C beam.

Method	FEM (180 × 12)	MsFEM-L (90 × 6)	MsFEM-P (90 × 6)
Elapsed (second)	38.12	16.22	16.01

Example 2: A MEE cantilever beam with holes under the temperature variation $\Delta T=10\text{K}$ is shown in Figure 7. The geometric size of the beam is L = 15 mm and the width H = 1.2 mm, which is a plane stress problem. In multiscale calculation, the entire model is divided into 50 × 4 coarse-grids, and each coarse-grid is divided into 52 fine-grids. The structural constraints of the left fixed end are $u_x=u_y=\varphi =\psi =0$ .

Figure 7. Geometry and grid division of the MEE cantilever beam with holes.

FEM (200 × 52), MsFEM-L (50 × 4), and MSFEMP (50 × 4) were used to solve the MEE cantilever beam, respectively. The generalized displacements (displacements u_x and u_y, electric potential $\varphi$ , and magnetic potential $\psi$ ) of AD edge obtained by three algorithms are shown in Figure 8. The finite element solution is used as the reference solution. The maximum errors of displacements u_x and u_y in multiscale results are 0.71% and 3.17%. The maximum errors of the electric potential and magnetic potential are 4.77% and 0.78%. Hence, the high precision of the multiscale analysis algorithm can be verified. From the comparison of these results, it can be seen that the MsFEM can capture the static response of each physical field more accurately.

Under identical hardware conditions, the calculation time of the three algorithms is compared, as shown in Table 2. The time of MsFEM-L and MsFEM-P is only 45.61% and 43.62% of that of FEM. From the above results, it can be seen that the MsFEM significantly reduces the computation amount and computation time, which verifies the efficiency of the multiscale algorithm.

(a) (b)

(c) (d)

Figure 8. Generalized displacements at edge AD of the MEE cantilever beam with holes. (a) displacement u_x; (b) displacement u_y; (c) electrical potential ϕ; (d) magnetic potential ψ.

Table 2. Comparing the duration of calculations for the MEE cantilever beam with holes.

Method	FEM (200 × 52)	MsFEM-L (50 × 4)	MsFEM-P (50 × 4)
Elapsed (second)	42.03	19.17	18.33

5. Conclusion

In this paper, the MsFEM is developed to solve the problem of multi-field coupling statics of heterogeneous MEE structures in thermal environment. After constructing three kinds of numerical base functions (displacement, electric potential, and magnetic potential base functions), the macroscopic equivalent stiffness matrix and load matrix of the mechanical-thermal-electro-magnetic coupling problem are derived. The heterogeneous MEE structure is discretized on the coarse scale, and the original problem is solved on the macroscopic scale, which saves a lot of computing resources. Finally, the correctness of the algorithms and programs as well as the high precision and efficiency of the algorithms in solving the heterogeneous MEE thermal-force-electric-magnetic coupling problem are verified by solving the problems of FGMEE C-C beams and MEE cantilever beams with holes.

Funding

This work was supported by the National Natural Science Foundation of China (Grant Numbers. 42102346, 42172301).

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Boomgaard, J.V.D., Van Run, A.M.J.G. and Van Suchtelen, J. (1976) Magnetoelectricity in Piezoelectric-Magnetostrictive Composites. Ferroelectrics, 10, 295-298. https://doi.org/10.1080/00150197608241997
[2]	Liu, J.X. (2010) Propagation of Elastic Waves in Magneto-Electro-Elastic Solids. Journal of Shijiazhuang Tiedao University (Natural Science Edition), 23, 5. https://doi.org/10.13319/j.cnki.sjztddxxbzrb.2010.02.001
[3]	Fetisov, L.Y., Serov, V.N., Chashin, D.V., Makovkin, S.A., Srinivasan, G., Viehland, D., et al. (2017) A Magnetoelectric Sensor of Threshold DC Magnetic Fields. Journal of Applied Physics, 121, Article ID: 154503. https://doi.org/10.1063/1.4981533
[4]	Liu, H. and Lv, Z. (2018) Uncertain Material Properties on Wave Dispersion Behaviors of Smart Magneto-Electro-Elastic Nanobeams. Composite Structures, 202, 615-624. https://doi.org/10.1016/j.compstruct.2018.03.024
[5]	Lang, Z. and Xuewu, L. (2013) Buckling and Vibration Analysis of Functionally Graded Magneto-Electro-Thermo-Elastic Circular Cylindrical Shells. Applied Mathematical Modelling, 37, 2279-2292. https://doi.org/10.1016/j.apm.2012.05.023
[6]	Arefi, M. and Zenkour, A.M. (2016) A Simplified Shear and Normal Deformations Nonlocal Theory for Bending of Functionally Graded Piezomagnetic Sandwich Nanobeams in Magneto-Thermo-Electric Environment. Journal of Sandwich Structures & Materials, 18, 624-651. https://doi.org/10.1177/1099636216652581
[7]	Farai, K.K. and Zhang, H. (2024) Global Sensitivity Analysis on Civil Engineering Structures (Bridges) Using the Homotopy Finite Element (HFE) Method. Open Access Library Journal, 11, 1-13. https://doi.org/10.4236/oalib.1111678
[8]	Du, J., Wu, J., Jing, L., Li, S. and Zhang, Q. (2023) Seismic Wavelet Analysis Based on Finite Element Numerical Simulation. Journal of Geoscience and Environment Protection, 11, 220-228. https://doi.org/10.4236/gep.2023.116015
[9]	Yang, F., Zhang, D., Li, L., Han, X., Zhong, Z. and Jiang, S. (2008) The Effective Properties of Smart Composites with Linear Coupling Behaviors. International Journal of Mechanics and Materials in Design, 4, 255-263. https://doi.org/10.1007/s10999-008-9059-1
[10]	Vinyas, M. and Harursampath, D. (2020) Computational Evaluation of Electro-Magnetic Circuits’ Effect on the Coupled Response of Multifunctional Magneto-Electro-Elastic Composites Plates Exposed to Hygrothermal Fields. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, 235, 2832-2850. https://doi.org/10.1177/0954406220954485
[11]	Vinyas, M. and Kattimani, S.C. (2017) Hygrothermal Analysis of Magneto-Electro-Elastic Plate Using 3D Finite Element Analysis. Composite Structures, 180, 617-637. https://doi.org/10.1016/j.compstruct.2017.08.015
[12]	Vinyas, M., Kattimani, S.C. and Joladarashi, S. (2018) Hygrothermal Coupling Analysis of Magneto-Electroelastic Beams Using Finite Element Methods. Journal of Thermal Stresses, 41, 1063-1079. https://doi.org/10.1080/01495739.2018.1447856
[13]	Zenkour, A.M. and Abbas, I.A. (2013) Magneto-Thermoelastic Response of an Infinite Functionally Graded Cylinder Using the Finite Element Method. Journal of Vibration and Control, 20, 1907-1919. https://doi.org/10.1177/1077546313480541
[14]	Abbas, I.A. and Zenkour, A.M. (2013) LS Model on Electro-Magneto-Thermoelastic Response of an Infinite Functionally Graded Cylinder. Composite Structures, 96, 89-96. https://doi.org/10.1016/j.compstruct.2012.08.046
[15]	Babuska, I., Caloz, G. and Osborn, J.E. (1992) Special Finite Element Methods for a Class of Second Order Elliptic Problems with Rough Coefficients. SIAM Journal on Numerical Analysis, 31, 945-981.
[16]	Hou, T.Y. and Wu, X. (1997) A Multiscale Finite Element Method for Elliptic Problems in Composite Materials and Porous Media. Journal of Computational Physics, 134, 169-189. https://doi.org/10.1006/jcph.1997.5682
[17]	Yang, D.S., Zhang, H.W., Zhang, S. and Lu, M.K. (2014) A Multiscale Strategy for Thermo-Elastic Plastic Stress Analysis of Heterogeneous Multiphase Materials. Acta Mechanica, 226, 1549-1569. https://doi.org/10.1007/s00707-014-1269-z
[18]	Xing, Y.F., Yang, Y. and Wang, X.M. (2010) A Multiscale Eigenelement Method and Its Application to Periodical Composite Structures. Composite Structures, 92, 2265-2275. https://doi.org/10.1016/j.compstruct.2009.08.006
[19]	Fu, P., Liu, H. and Chu, X. (2017) An Efficient Multiscale Computational Formulation for Geometric Nonlinear Analysis of Heterogeneous Piezoelectric Composite. Composite Structures, 167, 191-206. https://doi.org/10.1016/j.compstruct.2017.02.005
[20]	Ren, S., Meng, G., Cheng, F. and Zhou, L. (2020) Transient Responses of Functionally Graded Magneto-Electro-Elastic Structures with Holes in Thermal Environment Using Stabilized Node-Based Smoothed Radial Point Interpolation Method. International Journal of Mechanical Sciences, 185, Article ID: 105870. https://doi.org/10.1016/j.ijmecsci.2020.105870