© The Author(s), Published by EDP Sciences, 2024

1 Introduction

Noise emission from modern industrial activities has turned into an important socio-environmental issue in our increasingly noise-polluted communities. In particular, the periodic mechanical operation of impactive production machines in industrial building constructions and machine halls [1–5] can generate repetitive transient radiated noise levels that often exceed the ever more restrictive legal regulations on noise exposure [6, 7]. Furthermore, the alternating stresses and residual vibrations caused by the collaborative action of the inertial forces of moving parts, the working forces of periodic impacts, and the dynamic interaction forces between the surrounding medium and machine base, can result in malfunctions or instantaneous failures of these machines and/or building structural damages. Therefore, effective energy-efficient control methods should be pursued for proper protection against the resultant noise and vibration [8–13].

Double-wall structures have frequently been employed in many noise control engineering applications due to their undisputable advantage over single-leaf structures in producing superior noise insulation in modern buildings, automotive vehicles, passenger trains, as well as aerospace and marine structures. Therefore, double-wall partitions sound proofing has attracted a great deal of attention recently [13–19]. Many practical predictive models have been developed and a wide variety of passive, semi-active, fully active, and hybrid control methods have been implemented. They include but are not limited to use of viscoelastically damped or stiffened floating floors, Tuned Mass Dampers (TMD), NVA, and multi-functional electromagnetic, electrostatic, piezoelectric, and triboelectric actuating systems. Numerous investigations have been reported. For example, Oyelade proposed a mathematical model based on the weighted residual Galerkin method for the vibroacoustic problem of sound transmission across a rectangular double-wall panel with an integrated magnetic stiffening connection that functions based on the interaction energy between two rectangular magnets attached to the clamped plates [20]. Shao et al. presented an overview of the nonlinear characteristics and prospects on the future research of the NVA, and application of the NVA for low-frequency noise control within three dimensional acoustic cavities such as vehicle interior cabin [21]. Mao used a shunted EMD connected between incident and radiating plate to enhance sound transmission loss through a double-plate structure in the low frequency region [22]. Wrona et al. proposed a semi-active control approach for double-panel low frequency noise barriers based on action of bistable links mounted between the incident and the radiating panels that can structurally couple or decouple [14]. Hasheminejad and Vesal presented an idealized two dimensional transient vibroacoustic model for hybrid active/semi-active impact sound transmission dampening through a smart inter-coupled ElectroRheological/Piezoelectric (ERF/PZT) double-beam floor structure [13]. Also, in a series of papers, Hasheminejad and Jamalpoor investigated the 3D steady-state vibroacoustic control of hybrid double-wall composite structures subjected to external acoustic excitations [17–19].

With rapid irreversible depletion of non-renewable resources, it is necessary to look for alternative sustainable and environment-friendly sources of energy that are compatible with the requirements of our modern society. Such systems can be used for powering a wide variety of modern applications, including Microelectro-Mechanical Systems (MEMS), smart buildings, autonomous Wireless Sensor Networks (WSN), Internet of Things (IoT), Global Positioning Systems (GPS), portable/wearable electronic devices/instruments, and remote environmental (structural) condition (health) monitoring sensors [23–26]. In this context, affordable Vibration Energy Harvesting (VEH) systems have widely been employed in intelligent living and industrial applications for converting ambient environmental kinetic energies into electrical energy using various energy conversion mechanisms [27–29]. In particular, simple and robust VEH floors or pavement slabs based on various transduction principles have been designed and tested for high power density kinetic-to-electrical energy conversion [30–34]. For example, Muñoz et al. built and tested a low-cost energy harvesting floor based on macro fiber composites and PZT disc transducers [35]. More recently, Wang et al. proposed a free-vibration-type piezoelectric beam energy harvesting floor structure using an impacting force amplification mechanism [36]. Also, Zhong et al. experimentally examined an energy harvesting floor system based on multiple layers (arrays) of clamped-clamped piezoelectric beam structures. Experiments were performed to find out the optimal external load resistance for maximum output power [37].

Dual-functional vibration damping and energy harvesting systems are capable of concurrent mitigation of structural motion and conversion of a portion of the absorbed energy into electrical power. Numerous dual-purpose controllable vibration suppression and energy regenerating systems have recently been developed [38–41]. TMDs coupled with electromagnetic energy harvesters have been employed for achieving simultaneous vibration control and energy regeneration, as demonstrated in previous studies [42–45]. Replacing linear TMDs with NVAs has proved to effectively enhance their ability to dampen vibrations and broaden the frequency range of vibration attenuation with lighter auxiliary mass [46]. NVAs with essential nonlinear stiffness, known as Nonlinear Energy Sinks (NES) in academic literature, have also garnered significant attention across various engineering applications [47]. Several studies have delved into the effectiveness of the NESs in vibration mitigation and kinetic energy absorption across various oscillatory systems, including single or multi-degree of freedom setups, as well as continuous configurations [48–53]. Additionally, a series of research endeavors have explored the energy harvesting capabilities of NESs. Remik et al., for example, coupled primary linear oscillators with lightweight NESs to improve energy harvesting efficiency under impulsive excitation [54, 55]. Integration of NVAs can serve as a dual purpose mechanism to simultaneously enhance the vibration suppression performance and improve energy harvesting capabilities [56–58]. In a numerical investigation, Blanchard et al. outfitted a 3D model of a submerged structure coupled with a NES. Their findings revealed that integrating the NES not only effectively dampens vortex-induced oscillations of the system but also enables the extraction of a substantial amount of kinetic energy from the structure’s motion [59].

Double-floor structures radiate undesirable sound waves that propagate into the adjacent rooms when subjected to impact loading. According to the above brief overview of recent literature, there is a large body of investigations that utilize a wide range of purely passive, semi-active, or entirely active techniques for effective mitigation of inter-floor impact sound transmission. Despite the acoustic discomfort, substantial vibrations observed in floor structures under impact loading offer a promising energy source for harvesting. Nevertheless, there exists a notable research gap in systematic analytical or numerical investigation concerning simultaneous inter-floor impact sound transmission control and energy harvesting. The major aim of current work is to fill this important breach in the literature. To do this, current study describes the conceptual development of an analytical model of a dual-functional double plate floor structure optimally designed for hybrid regenerative control of inter-floor impact sound transmission. Accordingly, the 3D transient acousto-elastic model of a cavity-coupled impact-driven piezo-sandwich double-floor structure that is mechanically interconnected through a NVA-based multi-resonant shunted electromagnetic vibration absorber is first formulated. Subsequently, the complete set of NVA and shunt circuit parameters for the composite system are simultaneously optimized through an efficient multi-objective optimization scheme based on a quadratic performance criterion with RMS value of the source room center-point pressure and harvested energy serving as the competing cost functions. The suggested methodology allows the noise control and architectural engineers to concurrently profit from advantageous features of Genetic Algorithm (GA)-optimized hybrid energy harvesting and acoustic suppression systems in a single platform. It can particularly assist in development of modern innovative VEH floor systems with increased energy productivity and superior broadband impact noise control characteristics [30–32, 37]. Finally, the disclosed wide range of time response simulation data can provide a reliable benchmark for assessment of stringently numerical and/or approximate procedures and also supplement the future experimental studies.

2 Problem description

Consider a simply-supported flexible piezo-sandwich double plate structure of dimension a × b that is coupled to three parallelepiped cavities of heights h_η (η = 1, 2, 3) with acoustically rigid flat boundaries without any prevailing scattering objects. The specifics of proposed configuration and the implemented coordinate systems (x, y, z_η; η = 1, 2, 3) are illustrated in Figure 1, with (ρ₀, c₀) respectively denoting the medium density (air) and speed of sound. The sandwich plates are assumed to be fabricated from an elastic supporting base panel symmetrically joined with two fully-electroded uniform piezoelectric layers that are themselves connected to multi-resonant RLC shunt circuits in parallel configuration. The multi-resonant shunts are well-known to be more effective in both energy harvesting and vibration damping compared to the customary single resonant shunts [40, 60, 61]. In particular, for a given stiffness ratio, the parallel shunt mechanism has been proven to outperform the series circuit in addition to increased robustness and lower sensitivity with regard to system parameter changes [40]. Similarly, the double-floor assembly, which essentially differentiates the source and receiving cavities, is assumed to be intercoupled at an arbitrary location (x₀, y₀) with a light-weight NVA of mass m, damping c, and linear & nonlinear stiffness constants (k_L, k_NL). A multi-resonant electromagnetic RLC shunt damper is installed in parallel to the NVA damper, while the upper piezoelectric-based plate (plate 1) is excited by a transient transverse impact force F(x, y, t).

Figure 1 Problem configuration.

2.1 Basic governing equations

Following the conventional methodology in modeling general structure/fluid interaction problems [19, 62], one begins with the standard three dimensional wave equation to model the transient pressure wave fields within the top/gap/lower enclosure fluids (see Fig. 1):

$$\frac{{\partial }^2{p}_{\eta }}{\partial x^2}+\frac{{\partial }^2{p}_{\eta }}{\partial y^2}+\frac{{\partial }^2{p}_{\eta }}{\partial z_{\eta }^2}=\frac{1}{c_0^2}\frac{{\partial }^2{p}_{\eta }}{\partial t^2},$$(1)

where 0 ≤ z_η ≤ h_η denotes the local vertical coordinate, and p_η (x, y, z_η, t) (η = 1, 2, 3) respectively signify the sound pressure within the source, gap, and receiving enclosures. Next, by using the classical Hamilton’s variational principle, Maxwell’s electrodynamics equations, along with Kirchhoff thin plate model (see Appendix A), the main governing equations for the piezo-sandwich plates may be formulated in the compact form:

$$\begin{array}{c}(Q_{11}^{\mathrm{h}}{I}_1^{\mathrm{h}}+2Q_{11}^{\mathrm{p}}{I}_1^{\mathrm{p}})\frac{{\mathrm{\partial }}^4{w}_i}{\mathrm{\partial }{x}^4}+(Q_{22}^{\mathrm{h}}{I}_1^{\mathrm{h}}+2Q_{22}^{\mathrm{p}}{I}_1^{\mathrm{p}})\frac{{\mathrm{\partial }}^4{w}_i}{\mathrm{\partial }{y}^4}+2(Q_{12}^{\mathrm{h}}{I}_1^{\mathrm{h}}+2Q_{12}^{\mathrm{p}}{I}_1^{\mathrm{p}}+2Q_{66}^{\mathrm{h}}{I}_1^{\mathrm{h}}+4Q_{66}^{\mathrm{p}}{I}_1^{\mathrm{p}})\frac{{\mathrm{\partial }}^4{w}_i}{\mathrm{\partial }{x}^2\mathrm{\partial }{y}^2}+({\rho }_{\mathrm{h}}{t}_{\mathrm{h}}+2{\rho }_{\mathrm{p}}{t}_{\mathrm{p}}){\ddot{w}}_i\\ -({\rho }_{\mathrm{h}}{I}_1^{\mathrm{h}}+2{\rho }_{\mathrm{p}}{I}_1^{\mathrm{p}})(\frac{{\mathrm{\partial }}^2{\ddot{w}}_i}{\mathrm{\partial }{x}^2}+\frac{{\mathrm{\partial }}^2{\ddot{w}}_i}{\mathrm{\partial }{y}^2})-2e_{31}{I}_2^{\mathrm{p}}\frac{{\mathrm{\partial }}^2{\psi }_i}{\mathrm{\partial }{x}^2}-2e_{32}{I}_2^{\mathrm{p}}\frac{{\mathrm{\partial }}^2{\psi }_i}{\mathrm{\partial }{y}^2}-(t_{\mathrm{h}}+t_{\mathrm{p}})\{e_{31}[{\delta }^{\mathrm{\prime }}(x)-{\delta }^{\mathrm{\prime }}(x-a)][H(y)-H(y-b)]\\ +e_{32}[H(x)-H(x-a)][{\delta }^{\mathrm{\prime }}(y)-{\delta }^{\mathrm{\prime }}(y-b)]\}{V}_{\mathrm{P}\mathrm{Z}\mathrm{T}}^{(i)}=f_i^{\mathrm{n}\mathrm{e}\mathrm{t}}(x,y,t),\\ \begin{array}{c}{\xi }_{11}{I}_3^{\mathrm{p}}\frac{{\mathrm{\partial }}^2{\psi }_i}{\mathrm{\partial }{x}^2}+{\xi }_{22}{I}_3^{\mathrm{p}}\frac{{\mathrm{\partial }}^2{\psi }_i}{\mathrm{\partial }{y}^2}-{\xi }_{33}\frac{{\pi }^2}{t_{\mathrm{p}}^2}{I}_4^{\mathrm{p}}{\psi }_i=e_{31}\frac{\pi I_2^{\mathrm{p}}}{t_{\mathrm{p}}}\frac{{\mathrm{\partial }}^2{w}_i}{\mathrm{\partial }{x}^2}+e_{32}\frac{\pi I_2^{\mathrm{p}}}{t_{\mathrm{p}}}\frac{{\mathrm{\partial }}^2{w}_i}{\mathrm{\partial }{y}^2},\end{array}\end{array}$$(2)

where δ'(·) is the first derivative of Dirac delta function, H(·) is the Heaviside function, w_i(x, y, t) (i = 1, 2) denote the transverse displacements of top and bottom PZT-incorporated plates, and ψ_i(x, y, t) (i = 1, 2), signify the internal electric potential induced in piezoelectric layers of each plate, respectively. Also, $V_{\mathrm{PZT}}^{\left(i\right)}(i=\mathrm{1,2})$ is the total electric voltage generated in the top and bottom piezo-sandwich panels, (ρ, t) refer to the layer density and thickness, and (Q, e, ξ) are the corresponding elastic, piezoelectric and dielectric parameters, while the subscripts (p, h) stand for the piezoelectric and elastic base layers, respectively. Furthermore, the coefficients $I_1^{\mathrm{h}}$ and $I_{\mathrm{1,2},\mathrm{3,4}}^{\mathrm{p}}$ are defined in Appendix A, and $f_{\mathrm{1,2}}^{\mathrm{net}}\left(x,y,t\right)$ are the net external force applied on the corresponding piezo-sandwich panel that will be stated later.

Now, utilizing the relevant Kirchhoff’s current and voltage laws, the dynamic equations governing the multi-mode shunt circuits in top/bottom PZT plates and the interconnected EM damper can respectively be derived in the form [40, 60, 61] (see Fig. 1):

Shunted circuits (piezo-plates):

$$\begin{array}{c}{I}_{i1}+I_{i2}+2C_{\mathrm{p}}{\dot{V}}_{\mathrm{PZT}}^{\left(i\right)}=I_{\mathrm{PZT}}^{\left(i\right)},\\ V_{\mathrm{ij}}+L_{\mathrm{ij}}{\dot{I}}_{\mathrm{ij}}+R_{\mathrm{ij}}{I}_{\mathrm{ij}}=V_{\mathrm{PZT}}^{\left(i\right)},\\ C_{\mathrm{ij}}{\dot{V}}_{\mathrm{ij}}=I_{\mathrm{ij}}, \left(i=\mathrm{1,2}, j=\mathrm{1,2}\right),\end{array}$$(3a)

Shunted circuits (EM damper):

$$\begin{array}{c}\left(L_{31}+L_{\mathrm{m}}\right){\dot{I}}_{31}+L_{\mathrm{m}}{\dot{I}}_{32}+{R_{31}{I}_{31}+R_{\mathrm{m}}\left(I_{31}+I_{32}\right)+V_{31}=V}_{\mathrm{EM}},\\ L_{\mathrm{m}}{\dot{I}}_{31}+\left(L_{32}+L_{\mathrm{m}}\right){\dot{I}}_{32}+{R_{32}{I}_{32}+R_{\mathrm{m}}\left(I_{31}+I_{32}\right)+V_{32}=V}_{\mathrm{EM}},\\ C_{3j}{\dot{V}}_{3j}=I_{3j},\hspace{1em}(j=\mathrm{1,2}),\end{array}$$(3b)

where V_i1, V_i2, I_i1, I_i2 (i = 1,2,3) are the voltages and electric currents induced in the associated capacitors and shunt circuit branches, and the general (R, C, L) parameters refer to the relevant electric circuits resistor, capacitor, and inductor, respectively. Also, (R_m, L_m) signify the internal resistance and inductance of the electromagnetic transducer, while V_EM is the induced electric voltage. Moreover, C_p = ξ₃₃ab/t_p denotes the equivalent electric capacitance of a single PZT layer, and $I_{\mathrm{PZT}}^{\left(\mathrm{1,2}\right)}$ are the electric current generated in the top and bottom PZT panels, which may be defined as the time derivative of accumulated electric charge as [63–66]:

$$I_{\mathrm{PZT}}^{\left(i\right)}=\frac{\mathrm{d}}{\mathrm{d}t}\left(\underset{A_{\mathrm{t}}}{\overset{}{\int }}{\mathbf{D}}^{\left(i\right)}.\mathbf{n} \mathrm{d}A+\underset{A_{\mathrm{b}}}{\overset{}{\int }}{\mathbf{D}}^{\left(i\right)}.\mathbf{n} \mathrm{d}A\right), (i=\mathrm{1,2}),$$(4)

where D⁽ⁱ⁾, (i = 1, 2) is the electric charge density vector of each plate, which is defined in equation (A6) in Appendix A. (A_t, A_b) refer to the surface area of the top and bottom electrode surfaces of PZT-based panels, and n is the corresponding surface normal vector. Moreover, the electric voltage and the mechanical force produced by the electromagnetic transducer can be defined as [40]:

$$\begin{array}{c}{V}_{\mathrm{EM}}=K_{\mathrm{v}}\left[{\dot{w}}_1\left(x_0,y_0,t\right)-{\dot{w}}_2\left(x_0,y_0,t\right)\right],\\ f_{\mathrm{EM}}=K_{\mathrm{f}}\left(I_{31}+I_{32}\right),\end{array}$$(5)

where (K_v, K_f) are the voltage and force constants of electromagnetic transducer, respectively.

Next, the net exterior loads on the top and bottom piezo-sandwich panels can respectively be expressed in the form (see Fig. 1):

$$\begin{array}{c}{f}_1^{\mathrm{net}}\left(x,y,t\right)=F\left(x,y,t\right)+\delta \left(x-x_0\right)\delta \left(y-y_0\right)\left[f_{\mathrm{a}}^{\left(1\right)}\left(t\right)-f_{\mathrm{EM}}\right]-\left({p_1|}_{z_1=0}-{p_2|}_{z_2=h_2}\right),\\ f_2^{\mathrm{net}}\left(x,y,t\right)=\delta \left(x-x_0\right)\delta \left(y-y_0\right)\left[f_{\mathrm{a}}^{\left(2\right)}\left(t\right)+f_{\mathrm{EM}}\right]-\left({p_2|}_{z_2=0}-{p_3|}_{z_3=h_3}\right),\end{array}$$(6)

where the actual point loads transmitted to the upper and lower piezo-sandwich plates through the NVA assembly, ($f_{\mathrm{a}}^{\left(1\right)}, f_{\mathrm{a}}^{\left(2\right)})$, are [67]:

$$\begin{array}{c}{f}_{\mathrm{a}}^{\left(1\right)}\left(t\right)=k_{\mathrm{L}}∆w_1+k_{\mathrm{NL}}∆w_1^3+c∆{\dot{w}}_1,\\ f_{\mathrm{a}}^{\left(2\right)}\left(t\right)=k_{\mathrm{L}}∆w_2+k_{\mathrm{NL}}∆w_2^3+c∆{\dot{w}}_2,\end{array}$$(7)

where ∆w_i = q(t) − w_i (x₀, y₀, t) (i = 1, 2), while the transverse motion of the NVA mass, q(t), follows:

$$m\ddot{q}\left(t\right)+f_{\mathrm{a}}^{\left(1\right)}\left(t\right)+f_{\mathrm{a}}^{\left(2\right)}\left(t\right)=0.$$(8)

2.2 Fluid-structure compatibility and modal equations of motion

Supposing completely rigid acoustical boundaries for the internal/external surfaces of the source, gap, and receiving cavities shown in Figure 1 (i.e. at x = 0, a; 0 ≤ y ≤ b, 0 ≤ z_η ≤ h_η (η = 1, 2, 3), y = 0, b; 0 ≤ x ≤ a, 0 ≤ z_η ≤ h_η (η = 1, 2, 3); and z₁ = h₁, z₃ = 0; 0 ≤ x ≤ a, 0 ≤ y ≤ b), while enforcing the compatibility of normal accelerations at the pertinent interfaces of each piezo-sandwich panel with the neighboring acoustic fluids (i.e. at z₁ = 0, z₂ = , h₂, z₃ = h₃; 0 ≤ x ≤ a, 0 ≤ y ≤ b), one has [68]:

$$\begin{array}{c}{\frac{\partial p_1}{\partial x}|}_{x=0}={\frac{\partial p_1}{\partial x}|}_{x=a}={\frac{\partial p_1}{\partial y}|}_{y=0}={\frac{\partial p_1}{\partial y}|}_{y=b}={\frac{\partial p_1}{\partial z_1}|}_{z_1=h_1}=0,\\ {\frac{\partial p_2}{\partial x}|}_{x=0}={\frac{\partial p_2}{\partial x}|}_{x=a} ={\frac{\partial p_2}{\partial y}|}_{y=0}={\frac{\partial p_2}{\partial y}|}_{y=b}=0,\\ {\frac{\partial p_3}{\partial x}|}_{x=0}={\frac{\partial p_3}{\partial x}|}_{x=a}={\frac{\partial p_3}{\partial y}|}_{y=0}={\frac{\partial p_3}{\partial y}|}_{y=b}={\frac{\partial p_3}{\partial z_3}|}_{z_3=0}=0,\\ {\frac{\partial p_1}{\partial z_1}|}_{z_1=0}=-{\frac{\partial p_2}{\partial z_2}|}_{z_2=h_2}=-{\rho }_0{\ddot{w}}_1\left(x,t\right), {\frac{\partial p_2}{\partial z_2}|}_{z_2=0}=-{\frac{\partial p_3}{\partial z_3}|}_{z_3=h_3}=-{\rho }_0{\ddot{w}}_2\left(x,t\right).\end{array}$$(9)

Furthermore, assuming simply support boundaries for the piezo-sandwich panels, the pertinent electric potentials and displacements may advantageously be written as functions of relevant modal components in the usual normal form:

$$\begin{array}{c}{w}_i\left(x,y,t\right)=\sum _{n=1}^{\infty }\sum _{m=1}^{\infty }{W}_{i,\mathrm{nm}}\left(t\right)\sin (\mathrm{n\pi }x/a)\sin (\mathrm{m\pi }y/b),\\ {\psi }_i\left(x,y,t\right)=\sum _{n=1}^{\infty }\sum _{m=1}^{\infty }{\Psi }_{i,\mathrm{nm}}\left(t\right)\sin (\mathrm{n\pi }x/a)\sin (\mathrm{m\pi }y/b), (i=1, 2),\end{array}$$(10)

where W_i,nm(t), Ψ_i,nm(t), (i = 1, 2) are unknown time-dependent modal coefficients.

At this point, noting the special form of the pressure boundary conditions in equation (9), the (triple) finite Fourier cosine transform can advantageously be applied to the three-dimensional acoustic pressure fields of the cavities, p_η(x, y, z_η, t), (η = 1, 2, 3), to obtain:

$$\begin{array}{c}{\mathcal{F}}_c\left\{p_{\eta }\left(x,y,z_{\eta },t\right)\right\}=P_{\eta }\left(i,j,k,t\right)=\frac{8}{\mathrm{ab}{h}_{\eta }}\underset{0}{\overset{h_{\eta }}{\int }}{\int }_0^b\underset{0}{\overset{a}{\int }}{p}_{\eta }\left(x,y,z_{\eta },t\right)\cos \left(\frac{\mathrm{i\pi }x}{a}\right)\cos \left(\frac{\mathrm{j\pi }y}{b}\right)\cos \left(\frac{\mathrm{k\pi }{z}_{\eta }}{h_{\eta }}\right)\mathrm{d}x\mathrm{d}y\mathrm{d}{z}_{\eta },\\ {\mathcal{F}}_c^{-1}\left\{P_{\eta }\left(i,j,k,t\right)\right\}=p_{\eta }\left(x,y,z_{\eta },t\right)=\frac{P_{\eta }\left(\mathrm{0,0},0,t\right)}{8}+\frac{1}{4}\sum _{i=1}^{\infty }{P}_{\eta }\left(i,\mathrm{0,0},t\right)\cos \left(\frac{\mathrm{i\pi }}{a}x\right)+\frac{1}{4}\sum _{j=1}^{\infty }{P}_{\eta }\left(0,j,0,t\right)\cos \left(\frac{\mathrm{j\pi }}{b}y\right)\\ +\frac{1}{4}\sum _{k=1}^{\infty }{P}_{\eta }\left(\mathrm{0,0},k,t\right)\cos \left(\frac{\mathrm{k\pi }}{h_{\eta }}{z}_{\eta }\right)+\frac{1}{2}\sum _{i=1}^{\infty }\sum _{j=1}^{\infty }{P}_{\eta }\left(i,j,0,t\right)\cos \left(\frac{\mathrm{i\pi }}{a}x\right)\cos \left(\frac{\mathrm{j\pi }}{b}y\right)+\frac{1}{2}\sum _{j=1}^{\infty }\sum _{k=1}^{\infty }{P}_{\eta }\left(0,j,k,t\right)\cos \left(\frac{\mathrm{j\pi }}{b}y\right)\cos \left(\frac{\mathrm{k\pi }}{h_{\eta }}{z}_{\eta }\right)\\ +\frac{1}{2}\sum _{i=1}^{\infty }\sum _{k=1}^{\infty }{P}_{\eta }\left(i,0,k,t\right)\cos \left(\frac{\mathrm{i\pi }}{a}x\right)\cos \left(\frac{\mathrm{k\pi }}{h_{\eta }}{z}_{\eta }\right)+\sum _{i=1}^{\infty }\sum _{j=1}^{\infty }\sum _{k=1}^{\infty }{P}_{\eta }\left(i,j,k,t\right)\cos \left(\frac{\mathrm{i\pi }}{a}x\right)\cos \left(\frac{\mathrm{j\pi }}{b}y\right)\cos \left(\frac{\mathrm{k\pi }}{h_{\eta }}{z}_{\eta }\right),\end{array}$$(11)

where ${\mathcal{F}}_c^{-1}$ signifies the classical Fourier cosine inverse transform. Furthermore, using the customary Galerkin approach, one may expediently utilize equations (6) and (10) in equation (2) in order to obtain the modal form of governing equations for the piezo-sandwich plates. Consequently, by using the relevant mode shape orthogonality relations, and after some basic manipulations, one arrives at the key system of ODEs for the top and bottom piezo-sandwich plates in the form:

$$\begin{array}{c}{\alpha }_{\mathrm{nm}}{W}_{1,\mathrm{nm}}+{\beta }_{\mathrm{nm}}{\ddot{W}}_{1,\mathrm{nm}}=\sin \left(\frac{\mathrm{n\pi }{x}_0}{a}\right)\sin \left(\frac{\mathrm{m\pi }{y}_0}{b}\right)\left[k_{\mathrm{L}}∆w_1+k_{\mathrm{NL}}∆w_1^3+c∆{\dot{w}}_1-k_f\left(I_{31}+I_{32}\right)\right]-{\mathrm{{\rm Y}}}_{\mathrm{mn}}^{\left(1\right)}\left(t\right)\\ +\left(t_{\mathrm{h}}+t_{\mathrm{p}}\right)\left(e_{31}\frac{\mathrm{nb}}{\mathrm{ma}}+e_{32}\frac{\mathrm{ma}}{\mathrm{nb}}\right)\left[1-{\left(-1\right)}^n\right]\left[1-{\left(-1\right)}^m\right]V_{\mathrm{PZT}}^{\left(1\right)}+{\int }_0^b{\int }_0^{a}F\left(x,y,t\right)\sin \left(\frac{\mathrm{n\pi }}{a}x\right)\sin \left(\frac{\mathrm{m\pi }}{b}y\right)\mathrm{d}x\mathrm{dy},\\ {\alpha }_{\mathrm{nm}}{W}_{2,\mathrm{nm}}+{\beta }_{\mathrm{nm}}{\ddot{W}}_{2,\mathrm{nm}}=\sin \left(\frac{\mathrm{n\pi }{x}_0}{a}\right)\sin \left(\frac{\mathrm{m\pi }{y}_0}{b}\right)\left[k_{\mathrm{L}}∆w_2+k_{\mathrm{NL}}∆w_2^3+c∆{\dot{w}}_2+k_f\left(I_{31}+I_{32}\right)\right]-{\mathrm{{\rm Y}}}_{\mathrm{mn}}^{\left(2\right)}\left(t\right)\\ +\left(t_{\mathrm{h}}+t_{\mathrm{p}}\right)\left(e_{31}\frac{\mathrm{nb}}{\mathrm{ma}}+e_{32}\frac{\mathrm{ma}}{\mathrm{nb}}\right)\left[1-{\left(-1\right)}^n\right]\left[1-{\left(-1\right)}^m\right]V_{\mathrm{PZT}}^{\left(2\right)},\end{array}$$(12)

where the expressions for coefficients (α_nm, β_nm) along with the time-dependent terms ${\mathrm{{\rm Y}}}_{\mathrm{mn}}^{\left(\mathrm{1,2}\right)}\left(t\right)$ are provided in Appendix B. Successive use of Fourier cosine transform for the sound pressure in equation (1), and applying the compatibility relations (9) as well as the modal expansions (10), gives:

$$\begin{array}{c}{\ddot{P}}_1\left(i,j,k,t\right)+{\Omega }_{1\mathrm{ijk}}^2{P}_1\left(i,j,k,t\right)=\sum _{j=1,j\ne m}^{\infty }\sum _{i=1,i\ne n}^{\infty }\frac{16{\rho }_0{c}_0^2}{\mathrm{ab}{h}_1}{\Lambda }_{\mathrm{nmij}}{\ddot{W}}_{1,\mathrm{nm}}\left(t\right),\\ {\ddot{P}}_2\left(i,j,k,t\right)+{\Omega }_{2\mathrm{ijk}}^2{P}_2\left(i,j,k,t\right)=\sum _{j=1,j\ne m}^{\infty }\sum _{i=1,i\ne n}^{\infty }\frac{16{\rho }_0{c}_0^2{\left(-1\right)}^{k+1}}{\mathrm{ab}{h}_2}{\Lambda }_{\mathrm{nmij}}{\ddot{W}}_{1,\mathrm{nm}}\left(t\right)+\sum _{j=1,j\ne m}^{\infty }\sum _{i=1,i\ne n}^{\infty }\frac{16{\rho }_0{c}_0^2}{\mathrm{ab}{h}_2}{\Lambda }_{\mathrm{nmij}}{\ddot{W}}_{2,\mathrm{nm}}\left(t\right),\\ {\ddot{P}}_3\left(i,j,k,t\right)+{\Omega }_{3\mathrm{ijk}}^2{P}_3\left(i,j,k,t\right)=\sum _{j=1,j\ne m}^{\infty }\sum _{i=1,i\ne n}^{\infty }\frac{16{\rho }_0{c}_0^2{\left(-1\right)}^{k+1}}{\mathrm{ab}{h}_3}{\Lambda }_{\mathrm{nmij}}{\ddot{W}}_{2,\mathrm{nm}}\left(t\right),\end{array}$$(13)

where the coefficients Λ_nmij and ${\Omega }_{\eta \mathrm{ijk}}^2, (\eta =\mathrm{1,2},3)$ are defined in Appendix B. Finally, the key unknown modal coefficients {W_i,nm(t), P_η(i, j, k, t} along with NVA displacement, q(t), can be determined through simultaneous solution of system of coupled nonlinear ordinary differential equations (8), (12), and (13). Subsequently, the transverse displacement of each plate and acoustic pressure at any desired location can readily be obtained by using the series expansions (8)and 2nd of equations (11), respectively.

3 Optimization algorithm

In the present work, the Non-dominated Sorting Genetic Algorithm (NSGA-II) that iteratively pursues a set of acceptable trade-offs of two conflicting objectives [69, 70], is numerically implemented for simultaneous optimization of all system parameters. Total harvested energy, as one of the competing objectives, is determined by integrating the harvesting power over time, while the instantaneous power extracted by the transduction mechanism is often represented by the energy dissipated across an equivalent resistive element in the coupled circuit [42, 71–73]. The purpose of the optimization procedure is to find optimized values for the target variables $\overline{\mathbf{X}}$ = [m, k_L, k_NL, c, R_ij, C_ij, L_ij] (i = 1,2,3, j = 1,2) in order to maximize the total harvested energy, $J_{\mathrm{E}}\left(\overline{\mathbf{X}}\right),$ and minimize the root mean square of receiving room center-point sound pressure level (SPL), $J_{\mathrm{P}}\left(\overline{\mathbf{X}}\right),$ as mathematically described by the following cost functions:

$$\begin{array}{c}{J}_{\mathrm{E}}\left(\overline{\mathbf{X}}\right)={\int }_0^T\sum _{i=1}^3\sum _{j=1}^2{R}_{\mathrm{ij}}{I}_{\mathrm{ij}}^2\mathrm{d}t,\\ J_{\mathrm{P}}\left(\overline{\mathbf{X}}\right)=\sqrt{\frac{1}{T}{\int }_0^T{p}_3^2(\frac{a}{2},\frac{b}{2},\frac{h_3}{2},t)\mathrm{d}t},\end{array}$$(14)

Here, the choice of minimizing the receiving room center-point SPL as an optimization target is well-justified due to its critical role as a representative location that captures the symmetrical sound distribution, offers valuable initial insights into overall sound trends, and, when effectively optimized, potentially yields noticeable enhancements in the whole receiving room’s average SPL. Following the standard procedure [69, 70], the main concept of NSGA-II is to find a nondominated set of potential candidates for the target variables (Pareto front), and make improvements at each iteration step until the termination condition is met. Once the optimized values are obtained, they will be used as the input parameters for the main simulation results. The main steps of the optimization procedure are summarized in the flowchart of Figure 2. Also, the input parameters of optimization algorithm are given in Table 1. Lastly, the implemented optimization routes include: Maximum Pressure Mitigation-Open Circuit (MPMO), Maximum Pressure Mitigation-Shunted (MPMS), Maximum Energy Harvesting-Shunted (MEHS), and Dual-Objective Optimization System (DOOS).

Figure 2 Flowchart of the main steps in the NSGA-II optimization procedure.

4 Numerical results

Having noted the relatively wide variety of physical parameters employed in the above presented relatively complex acousto-elastic formulation, while realizing our computational power restrictions, we should confine our attention to some specific numerical examples. The material and geometrical parameters of the proposed double piezo-sandwich plate-cavity system used for the main numerical simulations are listed in Table 2. A dedicated code is developed in order to numerically solve the coupled non-linear ODE equations (8), (12), and (13). Convergence of numerical results was assured in a simple trial and error fashion, i.e., by gradually increasing the truncation constants and seeking for numerical consistency of the calculated solutions. Taking up to 10 acoustics and structural modes (n, m = 1, 2, 3, …, N_max = 10; i, j, k = 0, 1, 2, …, M_max = 9) (i.e., a maximum of 1000 acoustic modes along with 100 structural modes) were seen to produce accurate results (the convergence plots are omitted for briefness). The NVA/EMD-attachment point is a critical factor for its vibration mitigation performance. Evidently, as one gets closer to the mid-plate location (i.e., x₀ = 0.5a, y₀ = 0.5b), the more pronounced will be the vibration reduction effects.

Before presenting the key simulation results, the overall accuracy of the formulation is rigorously verified through finite element simulation. Accordingly, the time response evolution of center-point acoustic pressures in the source/receiving rooms and the center-point transverse plate displacements as well as that of the NVA are calculated for the open circuit configuration of the system with the following input parameters: The elastic host layer is composed of aluminum, while PZT-4 is employed for the PZT sandwich layers, the plane stress properties of which are detailed in [74]. In addition, a = 0.8 m, b = 0.5 m, h₁ = h₃ = 1 m, h₂ = 0.4 m; x₀ = 0.5 m, y₀ = 0.3 m, m = 1 kg, c = 0.2 Ns/m, K_v = K_f = 0, k_L = 2 kN/m and k_NL = 7 kN/m³, with the remaining properties as listed in Table 2. Here, a concentrated impulsive force on the top plate is applied at (x₁ = 0.6 m, y₁ = 0.2 m) with the mathematical form:

$$F\left(x,y,t\right)=-\delta \left(x-x_1\right)\delta \left(y-y_1\right)f\left(t\right)=-\delta \left(x-x_1\right)\delta \left(y-y_1\right)\{\begin{array}{c}\left(3.3\times {10}^4\right)\mathrm{N}\hspace{0.5em}t\le t_0\\ 0 { t}_0<t\end{array}$$(15)

where t₀ = 0.15 ms. In the FEM model, the physics modules of “Solid Mechanics,” “Pressure Acoustics, Transient” and “Electrostatics” are adopted as the main components of the numerical simulation. The real-time multi-physics interactions are realized by coupling the “Piezoelectric Effect” and “Acoustic-Structure Boundary” modules. The general geometry is composed of nine domains (i.e., three acoustic cavities, two elastic host layers of the plates, and four piezoelectric layer domains in the top/bottom plates). The attachment of the NVA to the double plate structure is arranged with a concentrated loading of an equivalent magnitude on each panel at the connection point based on equation (7), where the nonlinear dynamics of the NVA mass is coupled to the system governing equations through an “ODE Interface”. To perform the meshing, a web of rectangles is mapped on one face of the model, which is then swept through the whole body to yield 69120 cuboid elements with 2189821 degrees of freedom, as depicted in Figure 3. Lastly, the “Time Dependent Study” is conducted to obtain the numerical solution of the problem. Here, it should be noted that the FEM model employs a 3D elasticity methodology, which is different from the Kirchhoff thin plate assumption of the developed mathematical model. However, similar results are anticipated from the two approaches due to the thinness of the plates. The results, as shown in Figure 4, demonstrate good harmony with those calculated via the 3D-FEM methodology constructed in COMSOL Multiphysics (version 6.1) package [75].

Figure 3 3D finite element model developed in COMSOL Multiphysics software.

Figure 4 Time histories of center-point acoustic pressures in the source and receiving rooms and the center-point transverse plate & NVA displacements for a concentrated impulsive force applied on the top floor panel.

Figure 5 depicts the collection of Pareto optimum points in the objective-space [$J_{\mathrm{P}}\left(\overline{\mathbf{X}}\right),J_{\mathrm{E}}\left(\overline{\mathbf{X}}\right)]$ based on the NSGA-II algorithm for the shunted bare system (i.e., the shunted double PZT-incorporated plates in absence of NVA and/or EM damper), shunted NVA-equipped system, and shunted hybrid NVA/EMD-equipped system under a single central point impulsive force acting on the top plate at t = 0 for 0 ≤ t ≤ 0.2 s. Here, it is clear that all acquired nondominated solutions are consistently distributed and converge to the Pareto-optimal frontiers of the problem. Also, in each configuration, the “Green” circular marker completely dominates the other Pareto solutions for MPMS, while the “Red” marker dominates for MEHS. The “yellow” marker is roughly selected in order to represent a DOOS system that can achieve a reasonable trade-off between maximum pressure mitigation and maximum energy harvesting. Furthermore, the last three main columns in Table 3 presents the numerical values of the design parameters associated with the selected Pareto optimal solutions that were marked by the green (MPMS), red (MEHS), and yellow (DOOS) circular markers in Figure 5, respectively. The numerical values listed in the first main column of Table 3 are associated with the optimal MPMO systems. Lastly, the superior performance of the shunted hybrid NVA/EMD-equipped system in comparison to the other two configurations (in view of both maximum pressure mitigation and maximum energy harvesting) is apparent in the last row subplot of Figure 5. Here, a considerably more widely distributed Pareto solutions is observed in comparison to the sharper first row Pareto subplot associated with the bare system which is less sensitive to energy harvesting and/or pressure mitigation.

Figure 5 Pareto optimum points in the [$J_{\mathrm{P}}\left(\overline{\mathbf{X}}\right),J_{\mathrm{E}}\left(\overline{\mathbf{X}}\right)]$ objective-space for the shunted bare, shunted NVA-equipped, and shunted hybrid NVA/EMD-equipped configurations.

Figure 6 presents the time evolutions of different energy components in the optimized MPMO, MPMS, MEHS, and DOOS systems for a train of five central point impulsive forces acting on the top plate at t = 0, 0.2, 0.4, 0.6, 0.8 s (see equation (15)), based on the optimal input parameters listed in Table 3. Here, $J_{\mathrm{E}}\left(\overline{\mathbf{X}}\right)$ is the harvested energy (green lines), the input energy (black lines) refer to the mechanical energy injected into the system by the applied train of center-point impulsive loads, the absorbed energy (purple lines) denotes the total amount of energy absorbed/stored eighter in the mass-spring-damper or in the internal resistance and inductance of the EMD, and the residual energy (red lines) is the energy remaining in the coupled vibroacoustic system (i.e., the input energy minus the absorbed and harvested energies). According to the open circuit case results illustrated in Figure 6, the residual energy matches the input energy for the bare system configuration in the absence of any harvesting or absorbing mechanism. Introducing an optimized NVA results in energy absorption in the system, which manifests in the form of kinetic energy of the mass, potential energy stored in the spring, or dissipated energy in the damper. Further incorporation of an EMD into the system facilitates additional energy absorption, either stored in the inductance or dissipated as heat in the resistance, leading to a slight increase in absorbed energy.

Figure 6 Time evolutions of various energy components in the NSGA-optimized MPMO, MPMS, MEHS, and DOOS systems based on the calculated optimal target parameter values.

Table 4 presents the numerical values of total harvested energy, absorbed energy, and residual energy associated with the optimized MPMO, MPMS, MEHS, and DOOS systems for the bare, NVA-equipped, and hybrid NVA/EMD-equipped configurations as calculated at t = 0.2, 0.6 and 1.0 s. It is worthwhile to further discuss some key points here. In MPMS system for example, which prioritizes acoustic pressure mitigation rather than emphasizing energy absorption or harvesting, introduction of an NVA naturally increases the absorbed energy, as expected. However, its direct impact on harvested energy is not as straightforward to assess. As indicated in Table 4, the presence of the NVA slightly decreases harvested energy. This can be attributed to the energy absorbed by the NVA, leaving less kinetic energy in the system available for harvesting. On the other hand, while the presence of the EMD does not substantially impact the absorbed energy, it increases the harvested energy by introducing an additional harvesting shunt circuit. This effect can be better observed in Table 4, where, for example at t = 0.2, integrating the EMD into the NVA-equipped system boosts the harvested energy from 8.8 to 12.0 mJ (~36% increase). In contrast, the primary objective in the MEHS scenario is to maximize the harvested energy without addressing sound pressure mitigation or vibration suppression. Unlike its role in MPMO and MPMS cases, the NVA no longer aims to increase the absorbed energy in this scenario. Instead, it functions as an additional degree of freedom attached to the system, the parameters of which are specifically tailored to enhance the harvested energy. The impact of NVA on energy harvesting becomes more apparent when it is integrated into the bare system in order to maximize the harvested energy. In this configuration, the influence of the optimized NVA induces a shift in system dynamics and provides a notable boost in harvested energy, increasing from 10.2 to 13.3 mJ at 0.2 s according to Table 4 (~30% increase).

Furthermore, the bar-chart presented in Figure 7, which is constructed based on the data presented in Table 4, clearly illustrates the notable energy harvesting performance advantages of the MEHS-and DOOS-based systems, particularly for the hybrid NVA/EMD-equipped configuration in comparison to the MPMS system (i.e., up to 69% improvement in the total harvested energy for the hybrid MEHS system compared to the bare MPMS system at t = 1.0 s). Here, it should be mentioned that the objective in MPMS case is to minimize the acoustic pressure, and the parameters of the system are optimized to achieve this goal through absorbing or harvesting the energy. When the NVA is added to the structure of the bare system, although the amount of harvested energy slightly decreases (e.g. from 38.9 to 36.7 mJ at 0.6 s) but it also results in a notable amount of energy that is absorbed in the NVA. Therefore, addition of the NVA effectively assists in reducing the residual energy of the system (e.g. from 74.8 to 56.8 mJ at 0.6 s) and approach the objective of the optimization.

Figure 7 Bar-chart plot associated with the total harvested energy with different configurations.

Finally, Figure 8 compares the time evolutions of center-point SPL in the source/receiving rooms for the optimized MPMO, MPMS, MEHS, and DOOS systems associated with the bare, NVA-equipped, and hybrid NVA/EMD-equipped configurations. The proposed dual-functional strategy is adapted to efficiently mitigate acoustic pressure while concurrently harvesting energy, as evidenced by Figures 6 and 8. It is worthwhile to note that conventional active damping methods rely on external energy resources to mitigate acoustic pressure through reducing vibrations. Hence, the energy harvesting capability of the suggested approach makes it particularly appealing for applications, where access to energy resources is limited. Table 5 presents the RMS values of center-point acoustic pressure in the receiving room for the optimized MPMO, MPMS, MEHS, and DOOS systems associated with the bare, NVA-equipped, and hybrid NVA/EMD-equipped configurations as calculated at t = 0.2, 0.6 and 1.0 s. Moreover, the bar-chart presented in Figure 9, which is constructed using the data listed in Table 5, visibly illustrates the superior pressure mitigation performance advantages of the MPMS and DOOS systems, particularly for the hybrid NVA/EMD-equipped configuration, in comparison to the other two (MEHS, and MPMO) systems (i.e., up to 96% reduction in the SPL for the hybrid MPMS system compared to the bare MPMO system).

Figure 8 Time evolutions of center-point acoustic pressures in the source and receiving rooms for the NSGA-optimized MPMO, MPMS, MEHS, and DOOS systems associated with the bare, NVA-equipped, and hybrid NVA/EMD-equipped configurations.

Figure 9 Bar-chart plot associated with RMS value of acoustic pressure for different configurations.

5 Conclusions

This study describes the conceptual development and analytical modeling of a novel high-performance vibration-based hybrid synergetic electromagnetic/piezoelectric double plate floor structure for concurrent energy harvesting and inter-floor impact sound transmission control. An efficient multi-objective evolutionary algorithm (NSGA-II), which iteratively pursues a set of acceptable Pareto-optimal solutions, is implemented to balance an optimal trade-off between maximum total harvested energy and minimum receiving room center-point RMS acoustic pressure as two key competing objectives. Intensive numerical simulations reveal that the dual-functional elasto-acoustic response suppression and energy extraction performance can be substantially improved with the NSGA-optimized multi-mode shunted hybrid PZT/NVA/EMD-equipped energy transduction mechanism as compared to the stand-alone (single element) piezoelectric conversion devices. In particular, the so-named MPMS system is observed to decisively outperform the other Pareto solutions when the impact noise transmission control is of main concern. Also, the MEHS system is seen to dominate when the highest total energy extraction is desired. On the other hand, the DOOS system, which benefits from advantageous features of both systems, may intuitively be selected to achieve a reasonable trade-off between maximum pressure mitigation and maximum energy harvesting. Furthermore, the notable energy harvesting performance advantages of the MEHS-and DOOS-based systems, particularly for the hybrid PZT/NVA/EMD-equipped configuration in comparison to the MPMS system are observed (i.e., up to (69%, 58%) improvements in the total harvested energy for the hybrid (MEHS, DOOS) systems compared to the bare MPMS system). Moreover, the superior pressure mitigation performance benefits of the MPMS and DOOS systems, especially for the hybrid PZT/NVA/EMD-equipped configuration, in comparison to the other two (MEHS, and MPMO) systems are noted (i.e., about (96%, 93%) reductions in the SPL for the hybrid (MPMS, DOOS) systems compared to the bare MPMO system). In conclusion, the obtained results are very promising as the proposed hybrid NVA-coupled piezoelectric-electromagnetic transduction mechanisms are established to provide a novel approach for effective dual-functional impact sound transmission control and vibration energy harvesting with particular applications in marine and heavy industrial building installations. Future research may include incorporation of the proposed procedure in design of smart structures with sophisticated controllers as well as further optimization of the structural and controller design parameters in order to achieve higher power extraction and lower SPLs.

Conflicts of interest

Authors declare no conflict of interests.

Data availability statement

Data are available on request from the authors.

Appendix A

In this appendix, the governing equations of motion for the three-layer sandwich PZT-plate are derived. Accordingly, the local vertical coordinates for the piezo-sandwich plate layers are fixed with respect to their respective midplanes as follows (see Fig. 1):

$${\overline{z}}_1=z-\frac{1}{2}\left(t_{\mathrm{h}}+t_{\mathrm{p}}\right), {\overline{z}}_2=z, {\overline{z}}_3=z+\frac{1}{2}\left(t_{\mathrm{h}}+t_{\mathrm{p}}\right),$$(A1)

where the subscripts (1, 2, 3) signify the upper piezoelectric, elastic base and lower piezoelectric layers, respectively. Next, using the Kirchhoff-Love assumptions for thin plates, the respective layer displacement components can be stated in the form [76]:

$$\begin{array}{c}{u}_{\mathrm{i}}\left(x,y,{\overline{z}}_i,t\right)=-{\overline{z}}_i\frac{\partial w_1\left(x,y,t\right)}{\partial x},\\ v_{\mathrm{i}}\left(x,y,{\overline{z}}_i,t\right)=-{\overline{z}}_i\frac{\partial w_1\left(x,y,t\right)}{\partial y},\\ w_{\mathrm{i}}\left(x,y,{\overline{z}}_i,t\right)=w_1\left(x,y,t\right),\end{array}$$(A2)

where i = 1, 2, 3; and w(x, y, t) is the vertical plate displacement, and with presumption of plane stress, while utilizing equations (A1), (A2), the relevant stress-strain relationships for each layer can be expressed in the form [77, 78]:

$$\begin{array}{c}{\epsilon }_{\mathrm{xx}}^{\left(i\right)}=\frac{\partial u_i}{\partial x}=-{\overline{z}}_i\frac{{\partial }^2{w}_1}{\partial x^2}, {\sigma }_{\mathrm{xx}}^{\left(i\right)}=Q_{11}^{\left(i\right)}{\epsilon }_{\mathrm{xx}}^{\left(i\right)}-e_{31}^{\left(i\right)}{E}_z^{\left(i\right)},\\ {\epsilon }_{\mathrm{yy}}^{\left(i\right)}=\frac{\partial v_i}{\partial y}=-{\overline{z}}_i\frac{{\partial }^2{w}_1}{\partial y^2}, {\sigma }_{\mathrm{yy}}^{\left(i\right)}=Q_{22}^{\left(i\right)}{\epsilon }_{\mathrm{yy}}^{\left(i\right)}-e_{32}^{\left(i\right)}{E}_z^{\left(i\right)},\\ {\gamma }_{\mathrm{xy}}^{\left(i\right)}=\frac{\partial u_i}{\partial y}+\frac{\partial v_i}{\partial x}=-2{\overline{z}}_i\frac{{\partial }^2{w}_1}{\partial \mathrm{x\partial }y}, {\tau }_{\mathrm{xy}}^{\left(i\right)}=Q_{66}^{\left(i\right)}{\gamma }_{\mathrm{xy}}^{\left(i\right)},\end{array}$$(A3)

where i = 1, 2, 3; and (Q, e) refer to the pertinent elastic/piezoelectric parameters. Moreover, the electric field intensity E components in the piezo-skin layers may be stated in the usual gradient-potential form [79, 80]:

$$\begin{array}{c}{E}_x^{\left(i\right)}=-\frac{\partial {\phi }_{\mathrm{i}}}{\partial x},\\ E_y^{\left(i\right)}=-\frac{\partial {\phi }_{\mathrm{i}}}{\partial y},\\ E_z^{\left(i\right)}=-\frac{\partial {\phi }_{\mathrm{i}}}{\partial {\overline{z}}_i},\end{array}$$(A4)

where i = 1, 3. With subsequent use of Maxwell’s electrodynamics equations, the electric potential in the piezo-sandwich panel may be defined with a cosine variation in the thickness direction as [81,82]:

$${\phi }_{\mathrm{i}}\left(x,y,{\overline{z}}_i,t\right)=-\cos \left(\pi {\overline{z}}_i/t_{\mathrm{p}}\right){\psi }_i\left(x,y,t\right)+\frac{1}{t_{\mathrm{p}}}{\overline{z}}_i{V}_{\mathrm{PZT}}^{\left(1\right)},$$(A5)

Furthermore, the electric displacement components of the piezoelectric skins may correspondingly be stated as:

$$\begin{array}{c}{D}_x^{\left(i\right)}={\xi }_{11}{E}_x^{\left(i\right)},\\ D_y^{\left(i\right)}={\xi }_{22}{E}_y^{\left(i\right)},\\ D_z^{\left(i\right)}=e_{31}^{\left(i\right)}{\epsilon }_{\mathrm{xx}}^{\left(i\right)}+e_{32}^{\left(i\right)}{\epsilon }_{\mathrm{yy}}^{\left(i\right)}+{\xi }_{33}{E}_z^{\left(i\right)},\end{array}$$(A6)

where ξ refer to the dielectric parameter.

Now, one can enforce the well-known Hamilton’s variational principle to derive the dynamic model for the piezo-sandwich plates [76]:

$$\underset{t_1}{\overset{t_2}{\int }}\delta (T+W-U)\mathrm{d}t=0,$$(A7)

where the variations of kinetic & strain energies, and the external work are respectively written as:

$$\begin{array}{c}\delta T=\underset{A}{\overset{}{\iint }}{\int }_{\frac{t_{\mathrm{h}}}{2}}^{\frac{t_{\mathrm{h}}}{2}+t_{\mathrm{p}}}{\rho }_{\mathrm{p}}\left({\dot{u}}_1\delta {\dot{u}}_1+{\dot{v}}_1\delta {\dot{v}}_1+{\dot{w}}_1\delta {\dot{w}}_1\right)\mathrm{d}z\mathrm{d}A+\underset{A}{\overset{}{\iint }}{\int }_{-\frac{t_{\mathrm{h}}}{2}}^{\frac{t_{\mathrm{h}}}{2}}{\rho }_{\mathrm{h}}\left({\dot{u}}_2\delta {\dot{u}}_2+{\dot{v}}_2\delta {\dot{v}}_2+{\dot{w}}_2\delta {\dot{w}}_2\right)\mathrm{d}z\mathrm{d}A\\ +\underset{A}{\overset{}{\iint }}{\int }_{-\frac{t_{\mathrm{h}}}{2}-t_{\mathrm{p}}}^{-\frac{t_{\mathrm{h}}}{2}}{\rho }_{\mathrm{p}}\left({\dot{u}}_3\delta {\dot{u}}_3+{\dot{v}}_3\delta {\dot{v}}_3+{\dot{w}}_3\delta {\dot{w}}_3\right)\mathrm{d}z\mathrm{d}A,\\ \delta U=\underset{A}{\overset{a}{\iint }}\left\{{\int }_{\frac{t_{\mathrm{h}}}{2}}^{\frac{t_{\mathrm{h}}}{2}+t_{\mathrm{p}}}\left[{\sigma }_{\mathrm{xx}}^{\left(1\right)}\delta {\epsilon }_{\mathrm{xx}}^{\left(1\right)}+{\sigma }_{\mathrm{yy}}^{\left(1\right)}\delta {\epsilon }_{\mathrm{yy}}^{\left(1\right)}+{\tau }_{\mathrm{xy}}^{\left(1\right)}\delta {\gamma }_{\mathrm{xy}}^{\left(1\right)}-D_x^{\left(1\right)}\delta E_x^{\left(1\right)}-D_y^{\left(1\right)}\delta E_y^{\left(1\right)}-D_z^{\left(1\right)}\delta E_z^{\left(1\right)}\right]\mathrm{d}z+{\int }_{-\frac{t_{\mathrm{h}}}{2}}^{\frac{t_{\mathrm{h}}}{2}}\left[{\sigma }_{\mathrm{xx}}^{\left(2\right)}\delta {\epsilon }_{\mathrm{xx}}^{\left(2\right)}+{\sigma }_{\mathrm{yy}}^{\left(2\right)}\delta {\epsilon }_{\mathrm{yy}}^{\left(2\right)}+{\tau }_{\mathrm{xy}}^{\left(2\right)}\delta {\gamma }_{\mathrm{xy}}^{\left(2\right)}\right]\mathrm{d}z+{\int }_{-\frac{t_{\mathrm{h}}}{2}-t_{\mathrm{p}}}^{-\frac{t_{\mathrm{h}}}{2}}\left[{\sigma }_{\mathrm{xx}}^{\left(3\right)}\delta {\epsilon }_{\mathrm{xx}}^{\left(3\right)}+{\sigma }_{\mathrm{yy}}^{\left(3\right)}\delta {\epsilon }_{\mathrm{yy}}^{\left(3\right)}+{\tau }_{\mathrm{xy}}^{\left(3\right)}\delta {\gamma }_{\mathrm{xy}}^{\left(3\right)}-D_x^{\left(3\right)}\delta E_x^{\left(3\right)}-D_y^{\left(3\right)}\delta E_y^{\left(3\right)}-D_z^{\left(3\right)}\delta E_z^{\left(3\right)}\right]\mathrm{d}z\right\}\mathrm{d}A,\\ \delta W=\underset{A}{\overset{a}{\iint }}{f}_1^{\mathrm{net}}\delta w_1 \mathrm{d}A,\end{array}$$(A8)

in which $f_1^{\mathrm{net}}$ signifies the net external load acting on the top piezo-sandwich plate. After direct substitution of equations (A3) to (A6) into equation (A8) and the subsequent outcome into equation (A7), and going through the necessary integration by parts, the dynamic equations of motion for the top piezo-sandwich plate (plate 1) is determined as [63, 64]:

$$\begin{array}{c}\mathit{\delta }{\mathit{w}}_{\mathbf{1}}:\left(Q_{11}^{\mathrm{h}}{I}_1^{\mathrm{h}}+2Q_{11}^{\mathrm{p}}{I}_1^{\mathrm{p}}\right)\frac{{\partial }^4{w}_1}{\partial x^4}+\left(Q_{22}^{\mathrm{h}}{I}_1^{\mathrm{h}}+2Q_{22}^{\mathrm{p}}{I}_1^{\mathrm{p}}\right)\frac{{\partial }^4{w}_1}{\partial y^4}+2\left(Q_{12}^{\mathrm{h}}{I}_1^{\mathrm{h}}+2Q_{12}^{\mathrm{p}}{I}_1^{\mathrm{p}}+2Q_{66}^{\mathrm{h}}{I}_1^{\mathrm{h}}+4Q_{66}^{\mathrm{p}}{I}_1^{\mathrm{p}}\right)\frac{{\partial }^4{w}_1}{\partial x^2\partial y^2}+\left({\rho }_{\mathrm{h}}{t}_{\mathrm{h}}+2{\rho }_{\mathrm{p}}{t}_{\mathrm{p}}\right){\ddot{w}}_1-\left({\rho }_{\mathrm{h}}{I}_1^{\mathrm{h}}+2{\rho }_{\mathrm{p}}{I}_1^{\mathrm{p}}\right)\left(\frac{{\partial }^2{\ddot{w}}_1}{\partial x^2}+\frac{{\partial }^2{\ddot{w}}_1}{\partial y^2}\right)-2e_{31}{I}_2^{\mathrm{p}}\frac{{\partial }^2{\psi }_1}{\partial x^2}-2e_{32}{I}_2^{\mathrm{p}}\frac{{\partial }^2{\psi }_1}{\partial y^2}-\left(t_{\mathrm{h}}+t_{\mathrm{p}}\right)\left\{e_{31}\left[{\delta }^{\text{'}}\left(x\right)-{\delta }^{\text{'}}\left(x-a\right)\right]\left[H\left(y\right)-H\left(y-b\right)\right]+e_{32}\left[H\left(x\right)-H\left(x-a\right)\right]\left[{\delta }^{\text{'}}\left(y\right)-{\delta }^{\text{'}}\left(y-b\right)\right]\right\}{V}_{\mathrm{PZT}}^{\left(1\right)}=f_1^{\mathrm{net}}\left(x,y,t\right),\\ \mathit{\delta }{\mathit{\psi }}_{\mathbf{1}}:{\xi }_{11}{I}_3^{\mathrm{p}}\frac{{\partial }^2{\psi }_1}{\partial x^2}+{\xi }_{22}{I}_3^{\mathrm{p}}\frac{{\partial }^2{\psi }_1}{\partial y^2}-{\xi }_{33}\frac{{\pi }^2}{t_{\mathrm{p}}^2}{I}_4^{\mathrm{p}}{\psi }_1=e_{31}\frac{\pi I_2^{\mathrm{p}}}{t_{\mathrm{p}}}\frac{{\partial }^2{w}_1}{\partial x^2}+e_{32}\frac{\pi I_2^{\mathrm{p}}}{t_{\mathrm{p}}}\frac{{\partial }^2{w}_1}{\partial y^2},\end{array}$$(A9)

where the coefficients, $I_{\mathrm{1,2},\mathrm{3,4}}^{\mathrm{p}}$ and $I_1^{\mathrm{h}}$ are given as: