Three-Dimensional Multiferroic Structures under Time-Harmonic Loading

1  Introduction

In recent years, smart composite materials have found applications in ultrasonic imaging devices, space structures, and energy harvesting [1–4]. These smart materials are capable of altering their dimensions, shape, rigidity, and various mechanical properties when exposed to external temperature, electric field, mechanical field, and others.

Among the smart materials, multiferroic or magneto-electro-elastic (MEE) materials have energy conversion capacities [5]. While single-phase multiferroics have only weak magnetoelectric (ME) coupling, materials composed of piezoelectric (PE) and piezomagnetic (PM) phases can have strong ME coupling for converting energy from one form into another [6–8]. Their capacity to convert and manipulate energies makes them useful in diverse fields. From innovative sensors capable of detecting subtle changes in magnetic (or electric) fields to actuators designed to respond dynamically to external stimuli, these materials become important in diverse engineering fields. Their role in developing advanced transducers, efficient energy harvesting systems, and even cutting-edge medical devices underscores their significance in shaping the future of technology and engineering [9–11]. Malleron et al. [12] conducted an experimental study on magnetoelectric transducers for powering small biomedical devices. Li et al. [13] investigated the static behaviour of MEE materials in a hygrothermal environment using the multi-physical cell-based smoothed finite element method (MCS-FEM). Their study highlights the significance of employing MCS-FEM for solving multi-physical problems. Sasmal et al. [14] reviewed recent progress in flexible magnetoelectric composites and devices for next-generation wearable electronics, focusing on strategic fabrication techniques to improve the performance of flexible ME composites and devices.

To fully exploit the potential of these advanced materials, it is essential to accurately model their behaviour under various conditions. One of the most effective tools for this purpose is the fundamental solution, also called Green’s functions (GF). GFs are pivotal in various analytical and numerical techniques [15–17] as they provide a means to solve complex boundary-value problems involving different classes of materials such as MEE materials. Given the intricate coupling among magnetic, electric, and elastic fields in MEE solids, accurate determination and utilization of GFs become crucial. They help in understanding and predicting how these materials will respond to different stimuli, which is essential for designing applications [18–20]. For instance, in the development of advanced transducers and sensors, GFs enable precise calculations and optimizations necessary for high performance. This importance is particularly evident when smart materials are used in layered structures. The complexity of these structures requires precise solutions to boundary-value problems, which GFs can provide. Yang et al. [21] investigated the natural characteristics of anisotropic MEE multilayered plates using the analytical and finite element approach and drew their conclusion on three-layered plates with four different stacking sequences. Guided wave propagation in a multilayered MEE curved panel by the Chebyshev spectral elements method was discussed by Xiao et al. [22]. Vattree et al. [23] studied the semicoherent heterophase interfaces in 3D MEE multilayered composites under external loads. The appropriate structural interface conditions and the corresponding complicated boundary-value problem were solved using the mathematically elegant and computationally powerful Stroh formalism combined with the Fourier transform and dual-variable and position method (DVP). Singularity-free theory and adaptive finite element computations of arbitrarily shaped dislocation loop dynamics in 3D heterogeneous material structures are described by Vattré et al. [24]. Small-scale thermal analysis of piezoelectric–piezomagnetic Functionally Graded (FG) microplates using modified strain gradient theory was carried out by Hung et al. [25]. Kiran [26] studied thermal and hygrothermal buckling characteristics of porous magneto-electro-elastic skewed plates using third-order shear deformation theory. Additionally, Lyu et al. [27] conducted nonlinear dynamic modelling of geometrically imperfect MEE nanobeams made of functionally graded materials. Nonlinear isogeometric analysis of magneto-electro-elastic porous nanoplates was done by Phung et al. [28]. Circular loadings on the surface of an anisotropic and MEE half-space are investigated by Wang et al. [29]. Here, an analytical solution for a uniform and an indentation-type load is derived, distinguishing important features associated with the considered loadings. The different features of fracture analysis in PE and MEE materials are analyzed by Bui et al. [30] and Kiran et al. [31]. Various analytical and numerical methods have been proposed for solving the related problems in multilayered structures over the past decades [32–35]. In these studies, different real-life challenges such as imperfections, cracks, dynamic loads, and other factors are explored. These aspects are crucial when designing layered structures, as they significantly influence the performance and reliability of actual devices.

Previous studies have explored the essential nature of MEE multilayered structures. Although these investigations can extract the most essential nature of MEE multilayered structures, they were mostly limited to the static response under vertical loading, employing relatively complicated methods. However, the dynamic responses and real-world applications necessitate more sophisticated approaches, which GFs can facilitate. Understanding the dynamic behavior of MEE multilayered structures is crucial for their integration into modern technological systems. The ability to accurately predict and control their performance under various operational conditions can significantly improve the design and functionality of smart devices. GFs, with their advanced capabilities, provide a robust framework for modelling and analyzing these complex interactions.

In this paper, we analyze the time-harmonic response of 3D multilayered transversely isotropic (TI) MEE layered solids induced by arbitrary time-harmonic loads in both horizontal and vertical directions (including the electric load). While previous studies have primarily focused on the static response, understanding the time-harmonic response is important for many applications that involve dynamic and oscillatory loading conditions [36]. Time-harmonic loads are particularly important because they simulate real-world conditions where the materials are subjected to time-dependent forces, such as in ultrasonic imaging, vibration control, and energy harvesting systems [37]. Computing the time-harmonic response allows for the accurate prediction of transient behaviors and the dynamic interaction between magnetic, electric, and mechanical fields within MEE materials, which is vital for optimizing the performance and efficiency of devices.

Recent advancements in the fundamental solutions for layered structures, including the Fourier-Bessel series (FBS) vector function [38] and the DVP method [16], enhance the analytical precision in the study of MEE materials. The FBS and DVP methods exhibit computational efficiency and robustness, making them particularly advantageous for addressing complex dynamic problems. These methods facilitate the accurate modelling and analysis of transient behaviors in multilayered structures. The FBS system of vector function is attractive, with the decoupled solutions (LM- and N-type) in each layer being simply expressed by the eigenvalues and eigenvectors. Additionally, the DVP method is effective for handling layering due to its stability and efficiency. Furthermore, in terms of the FBS system, the solved coefficients are discrete Love numbers [39], which can be saved and repeatedly used on the surface of the layered solid. Since the GFs are just the summation of the Love numbers (multiplied by the given base functions), they can be accurately and fast calculated. In this paper, the GFs corresponding to the surface loading over a circular area will be derived in detail, with numerical examples showing the influence of loading types, frequencies and layering on the solutions.

The results of this study have applications in several engineering fields. Specifically, understanding the time-harmonic responses of MEE materials can significantly enhance the design and performance of devices subjected to dynamic and oscillatory loading conditions. For example, accurate modelling of time-harmonic responses in ultrasonic imaging can lead to better image resolution and deeper tissue penetration. In vibration control, improved response predictions can enhance the stability and effectiveness of damping systems. In energy harvesting systems, the ability to predict and optimize the dynamic interaction between different fields can lead to more efficient energy conversion and storage solutions. By addressing these practical applications, this research contributes to the advancement of smart material technologies and their integration into various high-performance engineering systems.

2  Statement of the Problem

Fig. 1 depicts the 3D structure, which consists of an n-layered solid made of transversally isotropic (TI) -MEE solid materials. The Cartesian coordinate system (x, y, z) is attached to its (x, y)-plane on the surface and the layered solid in the positive z-direction. Each layer is ordered downward, with Layer 1 on the top (with z0 = 0). Layer j is bounded by two interfaces, z = zj–1 and z = zj, with thickness hj = zj –zj–1 (j = 1, ..., n). The last interface of Layer n, i.e., zn, rests on the homogeneous MEE half-space. We assume that the adjacent layers are perfectly bonded, and that the uniform time-harmonic loads in both vertical and horizontal directions are applied on the surface of the structure within a circular region of radius a (see Fig. 1). Since the loads are proportional to eiωt where ω = 2πf in Hz = 1/s is the angular frequency, the factor eiωt in the response will be omitted.

Figure 1: Schematic of the 3D structure composed of TI-MEE layered solid subjected to time-harmonic loads on its surface in (a), three types of time-harmonic loads, i.e., vertical mechanical load pz, vertical electric displacement dz and horizontal mechanical load (px, py) applied on the surface plane (x, y) of the structure within the circle of radius a in (b)

2.1 Governing Equations

The problem presented above can be solved in terms of either Cartesian or cylindrical coordinate systems. Here, we employ the FBS vector system to solve the problem. Thus, the equations of motion and the Gauss law in the MEE solid without body forces, electric charge, and magnetic monopoles (proportional to eiωt) are given by [20]

∂σrr∂r+1r∂σrθ∂θ+∂σrz∂z+σr−σθ2+ρω2ur=0∂σrθ∂r+1r∂σθθ∂θ+∂σθz∂z+2σrθr+ρω2uθ=0∂σrz∂r+1r∂σθz∂θ+∂σzz∂z+σrzr+ρω2uz=0∂Dr∂r+1r∂Dθ∂θ+∂Dz∂z+Drr=0∂Br∂r+1r∂Bθ∂θ+∂Bz∂z+Brr=0(1)

where σij (N/m2) (i, j = r, θ, z) is mechanical stresses, Dj (C/m2) is electric displacement, and Bj (T = Wb/m2) is magnetic inductions; ρ (kg/m3) and ui (m) are mass density and elastic displacements.

The constitutive relation for TI-MEE materials (with its symmetry axis along z-axis) is given by [40]

σrr=c11γrr+c12γθθ+c13γzz−e31Ez−q31Hzσθθ=c12γrr+c11γθθ+c13γzz−e31Ez−q31Hzσzz=c13γrr+c13γθθ+c33γzz−e33Ez−q33Hzσrz=2c44γrz−e15Er−q15Hθσθz=2c44γθz−e15Eθ−q15Hrσrθ=2c66γrθDr=2e15γrz+η11Er+d11HrDθ=2e15γθz+η11Eθ+d11HθDz=e31γrr+e31γθθ+e33γzz+η33Ez+d33HzBr=e15γrz+η11Er+μ11HrBθ=e15γθz+η11Eθ+μ11HθBz=e31γrr+e31γθθ+e33γzz+η33Ez+μ33Hz(2)

where cij (N/m2), eij (C/m2), and εij (C2/N.m2) are the elastic moduli, PE tensor, and dielectric permittivity tensor, respectively; qij (N/A.m), dij (C/A.m), and μij (C2/N.m2) are the PM tensor, ME tensor and magnetic permeability tensor, respectively; γij (dimensionless), Ei (V/m) and Hi (A/m) are the strain tensor, electric field and magnetic field.

The strain components γij, electric field Ej and magnetic field Hj are related to the mechanical displacements ui, electric potential ϕ (V) and magnetic potential ψ (A) by [40]

γrr=∂ur∂r,γθθ=∂uθr∂r+urr,γzz=∂uz∂z,γθz=12(∂uθ∂z+∂uzr∂z)γrz=12(∂uz∂r+∂ur∂z),γrθ=12(∂urr∂θ+∂uθ∂r−uθr)Er=−∂ϕ∂r,Eθ=−1r∂ϕ∂θ,Ez=−∂ϕ∂zHr=−∂ψ∂r,Hθ=−1r∂ψ∂θ,Hz=−∂ψ∂z(3)

2.2 Continuity Condition on any Interface z=zk (Except for the Last Half-Space)

Since we have considered that the internal interfaces between layers are perfectly bonded, the mechanical displacements, electric potential, magnetic potential, tractions, electric displacements, and magnetic induction are continuous across it (e.g., at z = zk). Namely, we have

ui(r,θ,zk+)=ui(r,θ,zk−),i=r,θ,zϕ(r,θ,zk+)=ϕ(r,θ,zk−)ψ(r,θ,zk+)=ψ(r,θ,zk−)σiz(r,θ,zk+)=σiz(r,θ,zk−)Dz(r,θ,zk+)=Dz(r,θ,zk−)Bz(r,θ,zk+)=Bz(r,θ,zk−)(4)

In the last homogeneous TI-MEE half-space, it is required that the solution remains finite.

2.3 Fundamental Solution in Terms of FBS System and DVP Method

2.3.1 FBS System of Vector Functions

To find the solution to the aforementioned boundary-value problem, the Fourier-Bessel series (FBS) vector functions are introduced [38]

L(r,θ;λmk)=ezS(r,θ;λmk)M(r,θ;λmk)=∇S=(∂rer+r−1∂θeθ)SN(r,θ;λmk)=∇×(Sez)=(r−1∂θer−∂reθ)S(5)

where er, eθ and ez are the unit vectors along r-, θ-, and z-directions of the cylindrical system; λmk is the k-th zero of the Bessel function of order m, scaled by a large value R, i.e., Jm(λmkR)=0 with deformation vanishes when r > R. Notice that the vector systems L, M and N are orthogonal.

The scalar function in Eq. (5) is defined by

S(r,θ;λmk)=Jm(λmkr)eimθ/2π;m=0, ±1, ±2,…(6)

which satisfies the following Helmholtz equation:

∂2S∂r2+∂Sr∂r+∂2Sr2∂θ2+λmk2S=0(7)

Due to the orthogonal property of the FBS system (i.e., Eq. (5)), u, ϕ and ψ can be expanded as [38]

u(r,θ,z)=urer+uθeθ+uzez=∑m∑k[UL(m,k;z)L(r,θ;λmk)+UM(m,k;z)M(r,θ;λmk)+UN(m,k;z)N(r,θ;λmk)](8)

ϕ(r,θ,z)=∑m∑kΦ(z)S(r,θ;λmk)(9)

ψ(r,θ,z)=∑m∑kΨ(z)S(r,θ;λmk)(10)

Also, t (at z = constant plane), D and B can be expanded as [38]

t(r,θ,z)=σrzer+σθzeθ+σzzez=∑m∑k[TL(m,k;z)L(r,θ;λmk)+TM(m,k;z)M(r,θ;λmk)+TN(m,k;z)N(r,θ;λmk)](11)

D(r,θ,z)=Drer+Dθeθ+Dzez=∑m∑k[DL(m,k;z)L(r,θ;λmk)+DM(m,k;z)M(r,θ;λmk)+DN(m,k;z)N(r,θ;λmk)](12)

B(r,θ,z)=Brer+Bθeθ+Bzez=∑m∑k[BL(m,k;z)L(r,θ;λmk)+BM(m,k;z)M(r,θ;λmk)+BN(m,k;z)N(r,θ;λmk)](13)

where Um, Φ, Ψ, Tm, Dm, Bm (m = L, M and N) are the expansion coefficients. Notice that, upon determining these coefficients, the corresponding physical quantities can be obtained by taking the simple summation in Eqs. (8)–(13).

We now derive the ordinary differential equations for these expansion coefficients and their solutions. First, from the constitutive relations, we find

TL=−λ2c13UM+c33UL′+e33Φ′+q33Ψ′TM=c44(UL+UM′)+e15Φ+q15ΨTN=c44UN′(14)

DL=−λ2e31UM+e33UL′−η33Φ′−d33Ψ′(15a)

DM=e15U+e15UM′−η11Φ−d11Ψ(15b)

DN=e15UN′(15c)

BL=−λ2q31UM+q33UL′−d33Φ′−μ33Ψ′(16a)

BM=q15UL+q15UM′−d11Φ−μ11Ψ(16b)

BN=q15UN′(16c)

where a prime represents the derivative of the variable with respect to z.

Eqs. (14), (15a) and (16a) contain 10 expansion coefficients (i.e., UL, UM, UN, Φ, Ψ, TL, TM, TN, DL and BL) and need to be combined with the following 3 relations (from Eqs. (1)–(3)) to solve them

(c44+c13)UL′−λ2c11UM+c44UM′′+(e15+e31)Φ′+(q15+q31)Ψ′+ρω2UM=0(λ2UNc66−c44UN′′−ρω2UN)=0−λ2c44(UL+UM′)−λ2c13UM′+c33UL′′−λ2e15Φ+e33Φ′′−λ2q15Ψ+q33Ψ′′+ρω2UL=0e33UL′′−η33Φ′′−d33Ψ′′−λ2(e15+e31)UM′−λ2e15UL+λ2η11Φ+λ2d11ΨL=0q33UL′′−d33Φ′′−μ33Ψ′′−λ2(q15+q31)UM′−λ2q15UL+λ2d11Φ+λ2μ11Ψ=0(17)

Now Eqs. (14), (15a), (16a) and (17) contain a total of ten equations for the ten coefficients. They can be converted into a linear system of first-order differential equations. Furthermore, both the LM-type and N-type systems are decoupled, with the N-type being purely elastic, as listed below:

ddz[UNTN/λ]=λmk[01/c44c66−ρω2/λmk20][UNTN/λ](18)

ddz[ULM(z)TLM(z)]=λmk[R1−1R2TR1−1−R2R1−1R2T+R3+R4−R2R1−1][ULM(z)TLM(z)](19)

where

U(z)=[ULλmkUMΦΨ]T;T(z)=[TL/λmkTMDL/λmkBL/λmk]T

and

R1=[c330e33q330c4400e330−η33−d33q330−d33−μ33],R2=[0−c4400c130e31q310−e15000−q1500],R3=[c440e15q150c1100e150−η11−d11q330−d11−μ11],R4=−ρω2λmk2[1000010000000000]

Furthermore, the N-type solution is zero if the applied load is axis-symmetric.

2.3.2 Layer Solution

Solutions of Eqs. (18) and (19) can be simply expressed in terms of their eigenvalues and eigenvectors. The general solution in any layer, say in Layer j bounded by zj-1 and zj, i.e., zj−1≤z≤zj(1≤j≤n), can be found as below:

[λmkUN(z)TN(z)]=[11c44sN−c44sN][eλmksN(z−zj)00e−λmksN(z−zj−1)][cN+cN−](20)

where

sN=(c66−ρω2/λmk2)/c44

has a positive real part.

[ULM(z)TLM(z)]=[E11E12E21E22][⟨eλmks11(z−zj)⟩00⟨eλmks22(z−zj−1)⟩][c+c-](21)

where [E] is the eigenmatrix and

⟨eλmks11z⟩=diag[eλmks1zeλmks2zeλmks3zeλmks4z];⟨eλmks22z⟩=diag[eλmks5zeλmks6zeλmks7zeλmks8z]

with the first four eigenvalues (s1–s4) having positive real parts and the other four having negative real parts.

2.3.3 DVP Matrix Method

To find the solution in the multilayered TI-MEE solid, the DVP method is adopted. Thus, by following Pan et al. [38], we have the following layer relation:

[λmkUN(zj−1)TN(zj)]=[N11jN12jN21jN22j][λmkUN(zj)TN(zj−1)][ULM(zj−1)TLM(zj)]=[S11jS12jS21jS22j][ULM(zj)TLM(zj−1)](22)

where

[N11jN12jN21jN22j]=[e−λmksNhj1sNc44−sNc44e−λmksNhj][1e−λmksNhjsNc44e−λmksNhj−sNc44]−1[S11jS12jS13jS14j]=[E11⟨e−λmks1hj⟩E12E21E22⟨eλmks2hj⟩][E11E12⟨eλmks2hj⟩E21⟨e−λmks1hj⟩E22]−1

In a similar manner, the expansion coefficients for Layer j+1 can be obtained. Thus, by assuming the continuity condition between the two adjacent layers (i.e., Eq. (4)), we can determine the recursive relation between Layer j and Layer j+1, which is given by

[λmkUN(zj−1)TN(zj+1)]=[Nj:j+1][λmkUN(zj+1)TN(zj−1)][ULM(zj−1)TLM(zj+1)]=[Sj:j+1][ULM(zj+1)TLM(zj−1)](23)

where

[Nj:j+1]=[N11jN11j+1+N11jN12j+1(1−N21jN12j+1)−1N21jN11j+1N12j+N11jN12j+1(1−N21jN12j+1)−1N22jN21j+1+N22j+1(1−N21jN12j+1)−1N21jN11j+1N22j+1(1−N21jN12j+1)−1N22j][Sj:j+1]=[S11jS11j+1+S11jS12j+1(I−S21jS12j+1)−1S21jS11j+1S12j+S11jS12j+1(I−S21jS12j+1)−1S22jS21j+1+S22j+1(I−S21jS12j+1)−1S21jS11j+1S22j+1(I−S21jS12j+1)−1S22j]

The recursive relation (i.e., Eq. (23)) can be repeatedly applied in structure as long as the interfaces between the layers are perfect.

2.3.4 Boundary Conditions on the Surface

The boundary condition at z = z0 (Fig. 1) (in terms of the Cartesian coordinate system and proportional to eiωt) is given by [38]

σxz={pxr≤a0r>0},σyz={pyr≤a0r>0},σzz={pzr≤a0r>0},Dz={dzr≤a0r>0}(24)

In terms of the FBS system, we have obtained the following coefficients due to the applied loads on the surface of the structure [38]

TL(z0)=2πpzJ1(λ0ka)aλ0kN1,m=0DL(z0)=2πdzJ1(λ0ka)aλ0kN1,m=0TM(z0)=π2aJ1(λ1ka)N2(±px−ipy),m=±1TN(z0)=π2aJ1(λ1ka)N2(−ipx∓py),m=±1(25)

We consider the following three cases separately one by one here.

Case 1: Uniform vertical load pz. For this case, we have

TLM(z0)=[2πpzJ1(λ0ka)aλ0kN1000]t

Case 2: Uniform vertical electric displacement dz. For this case, we have

TLM(z0)=[002πdzJ1(λ0ka)aλ0kN10]t

Case 3: Uniform horizontal load (px, py) in both x- and y-directions. For this case, we have TLM(z0)=[0π2aJ1(λ1ka)N2(±px−ipy)00]t, for m = ±1

2.3.5 Expansion Coefficients on the Surface

With the surface traction coefficients, we apply the recursive relation for propagating the expansion coefficients from the surface z0 to the last layer interface zn to obtain the following relations:

[λmkUN(z0)TN(zn)]=[N11j:nN12j:nN21j:nN22j:n][λmkUN(zn)TN(z0)][ULM(z0)TLM(zn)]=[S11j:nS12j:nS21j:nS22j:n][ULM(zn)TLM(z0)](26)

On the last interface z = zn, the unknowns in Eq. (26) are related as

TN(zn)=[E22(E12)−1]UN(zn)TLM(zn)=[E22(E12)−1]ULM(zn)(27)

where Eij (Eij) are the eigenvector elements (submatrices) in the last MEE half-space.

Combined with Eq. (27), we can solve, from Eq. (26), the coefficients of displacements, electric potential, and magnetic potential on the surface of the layered structure, as given below.

Due to the applied vertical load (i.e., m = 0), we have only LM-type system, for which the expansion coefficients on the surface is given by

ULM(z0)=[S11j:n(E22(E12)−1−S21j:n)−1S22j:n+S12j:n]TLM(0)(z0)(28)

where

TLM(z0)=[2πpzJ1(λ0ka)aλ0kN1000]t;TLM(z0)=[0002πdzJ1(λ0ka)aλ0kN1]t

for the vertical mechanical load and vertical electric displacement, respectively.

Due to the applied horizontal load (i.e., for m = ±1) we have both N-type and LM-type systems, and for those, the coefficients on the surface are given by

{λ1kUN(z0)}=[N11j:n(E22(E12)−1−N21j:n)−1N22j:n+N12j:n]TN(1)(z0)ULM(z0)=[S11j:n(E22(E12)−1−S21j:n)−1S22j:n+S12j:n]TLM(1)(z0)(29)

and

{λ1kUN(z0)}=[N11j:n(E22(E12)−1−N21j:n)−1N22j:n+N12j:n]TN(−1)(z0)ULM(z0)=[S11j:n(E22(E12)−1−S21j:n)−1S22j:n+S12j:n]TLM(−1)(z0)(30)

where

TLM(1)(z0)=[0π2aJ1(λ1ka)N2(px−ipy)00]tTLM(−1)(z0)=[0π2aJ1(λ1ka)N2(−px−ipy)00]tTN(1)(z0)=π2aJ1(λ1ka)N2(px−ipy)TN(−1)(z0)=π2aJ1(λ1ka)N2(−px−ipy)

Once the coefficients on the surface are obtained, we can find the solution at any point on the surface of the structure by carrying out the summation (i.e., Eqs. (8)–(13)). To find the solutions at any specific z level, we just need to replace z0 with z in Eq. (26). Since coefficients are discreet numbers, they represent the exact Love numbers [39]. These Love numbers are advantageous to simplify the computation process, as they can be pre-calculated and stored for repeated use in calculating the GFs [41].

3  Numerical Studies

3.1 Validation and Convergence of the Problem

Before presenting numerical results, the present field solutions are validated against existing ones and found to be consistent with the solution reported by Pan et al. [38], Chu et al. [40] and Qi et al. [42]. Fig. 2 illustrates the accuracy of the current method as compared with those of Chu et al. [40] and Qi et al. [42].

Figure 2: Variation of electric displacement D¯z with respect to depth due to the electrical loading d¯z (=dz/pz) on the surface of the half-space structure [42] (a), Variation of stress σ¯zz with respect to the thickness of the surface layer due to applied mechanical load px on the surface of a two-layered structure for different field quantitates at the interface (r/a = 0 and r/a = 1) [40] (b)

A loading radius of a = 1 nm is located on the surface with its center at the origin of the coordinate system. In the calculation, a small material damping factor, i.e., β = 0.0001, is applied to the material. Namely, ckj, ekj, qkj, ηkj, μkj and dkj are replaced by ckj (1 + 2iβ), ekj (1 + 2iβ), qkj (1 + 2iβ), ηkj (1 + 2iβ), μkj (1 + 2iβ) and dkj (1 + 2iβ) during the numerical computation. Furthermore, R in the FBS is fixed at 1500 nm (i.e., R = 1500 a) and M = 9000. In presenting the numerical findings in the paper, dimensionless parameters are used, as defined below.

x¯=x/a; z¯=z/a; a¯=a/a; p¯x=px/cmax; dz¯=dz/emax; pz¯=pz/cmax; ux¯=ux/a; uz¯=uz/a; ϕ¯=ϕ/ϕmax; Ψ¯=Ψ/Ψmax; ω¯=ωa(ρmax/cmax); σ¯zz=σzz/cmax; Dz¯=Dz/emax; Bz¯=Bz/qmax.