Forced Vibration Behaviour of Elastically Constrained Graphene Origami-Enabled Auxetic Metamaterial Beams

Abstract

This paper explores the vibration behaviour of an elastically constrained graphene origami-enabled auxetic metamaterial beam subject to a harmonic external force. The effective mechanical properties of the metamaterial are approximated using a micromechanical model trained via a genetic algorithm provided in the literature. The three coupled equations of motion are solved numerically; a set of trigonometric functions is used to approximate the displacement components. The accuracy of the proposed model is confirmed by comparing it with the natural frequencies of a simplified non-metamaterial structure available in the literature. Following this validation, the investigation extends to investigate the forced vibration response of the metamaterial beam, examining the influence of the graphene origami distribution pattern and content, graphene folding degree, linear and shear layer stiffness, and geometrical parameters on the dynamic behaviour of the structure. The results generally highlight the considerable effect of the shear layer, modelled as a Pasternak foundation, on the vibration behaviour of the elastically constrained metamaterial beams.

1. Introduction

Beams, which are fundamental components in many engineering designs, often experience various vibrational loads due to environmental or operational conditions. To ensure safety and performance, it is essential to analyse their vibration behaviour [1,2,3,4].

The recent advent of metamaterials has provided an opportunity for better vibration design of systems. These structures, designed with properties not usually found in conventional structures, are becoming highly valuable in engineering because they can be tailored to display unique behaviours, for instance in vibration isolation, through their unique characteristics such as negative thermal expansion, negative Poisson’s ratio, etc. [5,6,7,8]. Such characteristics make metamaterials particularly useful for improving how structures respond to vibrations. By using metamaterial beams, we can achieve better control over vibration frequencies and improve the overall performance of engineering systems [9,10,11,12,13].

External constraints on beams and plates are often modelled as elastic foundations. The Winkler and Pasternak models are two of the most well-known methods for representing elastic foundations. The Winkler model treats the foundation as a series of independent springs, where each point reacts only to the local deflection of the beam. The Pasternak model improves this by introducing shear interaction between points, offering a more realistic representation of how elastic media behave in real-world conditions. Together, these models allow for more accurate analysis of structures on elastic supports, making them particularly useful for studying vibration. Kacar et al. [14], for example, studied the free vibration behaviour of beams resting on a variable Winkler elastic foundation via the differential transform method. Kumar [15] examined the free vibration analysis of composite beams resting on a variable Pasternak foundation using the Rayleigh–Ritz method. Avcar and Mohammed [16] presented the free vibration analysis of functionally graded beams resting on a Winkler–Pasternak foundation using the separation of variables method. The forced vibrations of an elastically constrained cracked double-beam system interconnected by a viscoelastic layer were studied by Cehn et al. [17] via the development of a closed-form solution. A comprehensive study for the bending, buckling, and free vibration behaviours of materially imperfect beams resting on elastic foundation was performed in [18].

The integration of metamaterial beams with Winkler–Pasternak elastic constraints subject to external dynamic loads represents a promising fundamental research. While beams on elastic foundations have been studied extensively, metamaterial beams on elastic foundations remain largely unexplored. This presents a gap in the current understanding, particularly when it comes to forced vibrations, which are highly relevant in many practical applications; the forced vibration of beams supported by elastic foundations is of great importance in engineering fields such as civil and mechanical engineering. Metamaterials are also increasingly being used in advanced engineering sectors such as automotive and civil engineering due to their ability to control vibrations and improve the resilience of structures [19,20,21].

The aim of this paper is to address this research gap by conducting an in-depth analysis of the forced vibration behaviour of metamaterial beams constrained by a Winkler–Pasternak foundation. The metamaterial beam model has been adopted from our previous study [10] and is modified to take into account the effect of the Winkler–Pasternak foundation and a harmonic external force. Then, it is used to study how the properties of the metamaterial beam and the stiffness of the foundation influence key vibration characteristics and the time-dependent central deflection under the excitation force. The Winkler–Pasternak foundation is modelled as an elastic support, and a range of parametric studies are conducted to investigate the effects of varying foundation stiffness, metamaterial properties, and external forcing on the vibration behaviour of the beam.

2. Forced Vibration Metamaterial Beam Model

Metamaterial Beam

Consider an elastically constrained graphene origami-enabled auxetic metamaterial beam as shown in Figure 1. The multilayer system is characterised by a length L, a width b, and a thickness h. In the current study, the concentration of graphene origami is distributed with respect to two different patterns, Pattern-U and Pattern-A, with an even number of layers, mathematically, where [22]

$$\mathrm{U}-\mathrm{Pattern}:V_{\mathrm{Gr}}^{(k)}=V_{\mathrm{Gr}},$$

(1)

$$\mathrm{A}-\mathrm{Pattern}:V_{\mathrm{Gr}}^{(k)}=V_{\mathrm{Gr}}(2k-1)/N_{\mathrm{L}},$$

(2)

where the volume fraction of the graphene origami is $V_{\mathrm{Gr}}^{(k)}$, and N_L is the total layer number. The $V_{\mathrm{Gr}}$ (total volume fraction of the graphene origami) is defined by [22]

$$V_{\mathrm{Gr}}={\displaystyle \frac{W_{\mathrm{Gr}}}{W_{\mathrm{Gr}}+({\rho }_{\mathrm{Gr}}/{\rho }_{\mathrm{Cu}})(1-W_{\mathrm{Gr}})}},$$

(3)

where [W_Gr, ρ_Gr] is the [weight fraction, density] of nanofiller, and ρ_Cu is the matrix density. The effective mechanical properties of metamaterial are determined by [23]

$${\rho }_i=({\rho }_{\mathrm{Gr}}{V}_{\mathrm{Gr}}+{\rho }_{\mathrm{Cu}}{V}_{\mathrm{Cu}})\times f_{\rho }(V_{\mathrm{Gr}},T),$$

(4)

$${\nu }_i=({\nu }_{\mathrm{Gr}}{V}_{\mathrm{Gr}}+{\nu }_{\mathrm{Cu}}{V}_{\mathrm{Cu}})\times f_{\nu }(H_{\mathrm{Gr}},V_{\mathrm{Gr}},T),$$

(5)

$$E_i={\displaystyle \frac{1+\xi \eta V_{\mathrm{Gr}}}{1-\eta V_{\mathrm{Gr}}}}{E}_{\mathrm{Cu}}\times f_E(H_{\mathrm{Gr}},V_{\mathrm{Gr}},T),$$

(6)

where [ρ_i, ν_i, E_i] is [mass density, Poisson’s ratio, Young’s modulus] of the graphene origami-enabled auxetic metamaterial, and the modification functions are denoted by f_ν, f_ρ, and f_E. Geometrical parameters (η; ξ) are formulated as [23]

$$\eta ={\displaystyle \frac{(E_{\mathrm{Gr}}/E_{\mathrm{Cu}})-1}{(E_{\mathrm{Gr}}/E_{\mathrm{Cu}})+\xi }},$$

(7)

$$\xi =2({\displaystyle \frac{{\mathcal{l}}_{\mathrm{Gr}}}{t_{\mathrm{Gr}}}}),$$

(8)

in which Young’s moduli of the matrix and the graphene are shown by E_Cu and E_Gr, respectively; ℓ_Gr and t_Gr are the length and the thickness of the graphene, and [23]

$$\begin{array}{l}{f}_{\nu }(H_{\mathrm{Gr}},V_{\mathrm{Gr}},T)=1.01-1.43V_{\mathrm{Gr}}+0.165({\displaystyle \raisebox{1ex}{$T$}\!\left/ \!\raisebox{-1ex}{$T_0$}\right.})-16.8H_{\mathrm{Gr}}{V}_{\mathrm{Gr}}\\ -1.1H_{\mathrm{Gr}}{V}_{\mathrm{Gr}}({\displaystyle \raisebox{1ex}{$T$}\!\left/ \!\raisebox{-1ex}{$T_0$}\right.})+16H_{\mathrm{Gr}}^2{V}_{\mathrm{Gr}}^2,\end{array}$$

(9)

$$f_{\rho }(V_{\mathrm{Gr}},T)=1.01-2.01V_{\mathrm{Gr}}^2-0.0131({\displaystyle \raisebox{1ex}{$T$}\!\left/ \!\raisebox{-1ex}{$T_0$}\right.}),$$

(10)

$$\begin{array}{l}{f}_E(H_{\mathrm{Gr}},V_{\mathrm{Gr}},T)=1.11-1.22V_{\mathrm{Gr}}-0.134({\displaystyle \raisebox{1ex}{$T$}\!\left/ \!\raisebox{-1ex}{$T_0$}\right.})+0.559V_{\mathrm{Gr}}({\displaystyle \raisebox{1ex}{$T$}\!\left/ \!\raisebox{-1ex}{$T_0$}\right.})\\ -5.5H_{\mathrm{Gr}}{V}_{\mathrm{Gr}}+38H_{\mathrm{Gr}}{V}_{\mathrm{Gr}}^2-20.6H_{\mathrm{Gr}}^2{V}_{\mathrm{Gr}}^2,\end{array}$$

(11)

where T/T₀ = 1, which is assumed to represent room temperature (300 K), and H_Gr denotes the folding degree of graphene origami.

The equations of motion used in this paper have been adopted from our previous study [10] and have been modified to account for the effect of the Winkler–Pasternak foundation and a harmonic external force; this gives

$$\begin{array}{l}\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\frac{E_i(k)}{1-{\nu }_i^2(k)}dz}}\right\}{\displaystyle \frac{{\partial }^{2}u}{\partial x^2}}+\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\frac{E_i(k)}{1-{\nu }_i^2(k)}zdz}}\right\}{\displaystyle \frac{{\partial }^2\phi }{\partial x^2}}\\ -c_1\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\frac{E_i(k)}{1-{\nu }_i^2(k)}{z}^{3}dz}}\right\}{\displaystyle \frac{{\partial }^2\phi }{\partial x^2}}-c_1\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\frac{E_i(k)}{1-{\nu }_i^2(k)}{z}^{3}dz}}\right\}{\displaystyle \frac{{\partial }^{3}w}{\partial x^3}}\\ =\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\rho (k)dz}}\right\}{\displaystyle \frac{{\partial }^{2}u}{\partial t^2}}+\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\rho (k)zdz}}\right\}{\displaystyle \frac{\partial \phi }{\partial t^2}}\\ -c_1\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\rho (k)z^{3}dz}}\right\}{\displaystyle \frac{\partial \phi }{\partial t^2}}-c_1\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\rho (k)z^{3}dz}}\right\}{\displaystyle \frac{{\partial }^{3}w}{\partial x\partial t^2}},\end{array}$$

(12)

$$\begin{array}{l}\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\frac{E_i(k)}{2+2{\nu }_i(k)}dz}}\right\}{\displaystyle \frac{\partial \phi }{\partial x}}-c_2\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\frac{E_i(k)}{2+2{\nu }_i(k)}{z}^{2}dz}}\right\}{\displaystyle \frac{\partial \phi }{\partial x}}\\ -c_2\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\frac{E_i(k)}{2+2{\nu }_i(k)}{z}^{2}dz}}\right\}{\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}+\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\frac{E_i(k)}{2+2{\nu }_i(k)}dz}}\right\}{\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}\\ -c_2\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\frac{E_i(k)}{2+2{\nu }_i(k)}{z}^{2}dz}}\right\}{\displaystyle \frac{\partial \phi }{\partial x}}+c_2^2\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\frac{E_i(k)}{2+2{\nu }_i(k)}{z}^{4}dz}}\right\}{\displaystyle \frac{\partial \phi }{\partial x}}\\ +c_2^2\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\frac{E_i(k)}{2+2{\nu }_i(k)}{z}^{4}dz}}\right\}{\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}-c_2\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\frac{E_i(k)}{2+2{\nu }_i(k)}{z}^{2}dz}}\right\}{\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}\\ +c_1\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\frac{E_i(k)}{1-{\nu }_i^2(k)}{z}^{3}dz}}\right\}{\displaystyle \frac{{\partial }^{3}u}{\partial x^3}}+c_1\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\frac{E_i(k)}{1-{\nu }_i^2(k)}{z}^{4}dz}}\right\}{\displaystyle \frac{{\partial }^3\phi }{\partial x^3}}\\ -c_1^2\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\frac{E_i(k)}{1-{\nu }_i^2(k)}{z}^{6}dz}}\right\}{\displaystyle \frac{{\partial }^3\phi }{\partial x^3}}-c_1^2\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\frac{E_i(k)}{1-{\nu }_i^2(k)}{z}^{6}dz}}\right\}{\displaystyle \frac{{\partial }^{4}w}{\partial x^4}}\\ -k_{\mathrm{W}}w+k_{\mathrm{P}}{\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}+q=\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\rho (k)dz}}\right\}{\displaystyle \frac{{\partial }^{2}w}{\partial t^2}}+c_1\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\rho (k)z^{3}dz}}\right\}{\displaystyle \frac{{\partial }^{3}u}{\partial x\partial t^2}}\\ +c_1\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\rho (k)z^{4}dz}}\right\}{\displaystyle \frac{{\partial }^3\phi }{\partial x\partial t^2}}-c_1^2\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\rho (k)z^{6}dz}}\right\}{\displaystyle \frac{{\partial }^3\phi }{\partial x\partial t^2}}\\ -c_1^2\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\rho (k)z^{6}dz}}\right\}{\displaystyle \frac{{\partial }^{4}w}{\partial x^2\partial t^2}},\end{array}$$

(13)

$$\begin{array}{l}\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\frac{E_i(k)}{1-{\nu }_i^2(k)}zdz}}\right\}{\displaystyle \frac{{\partial }^{2}u}{\partial x^2}}+\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\frac{E_i(k)}{1-{\nu }_i^2(k)}{z}^{2}dz}}\right\}{\displaystyle \frac{{\partial }^2\phi }{\partial x^2}}\\ -c_1\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\frac{E_i(k)}{1-{\nu }_i^2(k)}{z}^{4}dz}}\right\}{\displaystyle \frac{{\partial }^2\phi }{\partial x^2}}-c_1\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\frac{E_i(k)}{1-{\nu }_i^2(k)}{z}^{4}dz}}\right\}{\displaystyle \frac{{\partial }^{3}w}{\partial x^3}}\\ -\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\frac{E_i(k)}{2+2{\nu }_i(k)}dz}}\right\}\phi +c_2\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\frac{E_i(k)}{2+2{\nu }_i(k)}{z}^{2}dz}}\right\}\phi \\ +c_2\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\frac{E_i(k)}{2+2{\nu }_i(k)}{z}^{2}dz}}\right\}{\displaystyle \frac{\partial w}{\partial x}}-\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\frac{E_i(k)}{2+2{\nu }_i(k)}dz}}\right\}{\displaystyle \frac{\partial w}{\partial x}}\\ +c_2\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\frac{E_i(k)}{2+2{\nu }_i(k)}{z}^{2}dz}}\right\}\phi -c_2^2\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\frac{E_i(k)}{2+2{\nu }_i(k)}{z}^{4}dz}}\right\}\phi \\ -c_2^2\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\frac{E_i(k)}{2+2{\nu }_i(k)}{z}^{4}dz}}\right\}{\displaystyle \frac{\partial w}{\partial x}}+c_2\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\frac{E_i(k)}{2+2{\nu }_i(k)}{z}^{2}dz}}\right\}{\displaystyle \frac{\partial w}{\partial x}}\\ -c_1\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\frac{E_i(k)}{1-{\nu }_i^2(k)}{z}^{3}dz}}\right\}{\displaystyle \frac{{\partial }^{2}u}{\partial x^2}}-c_1\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\frac{E_i(k)}{1-{\nu }_i^2(k)}{z}^{4}dz}}\right\}{\displaystyle \frac{{\partial }^2\phi }{\partial x^2}}\\ +c_1^2\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\frac{E_i(k)}{1-{\nu }_i^2(k)}{z}^{6}dz}}\right\}{\displaystyle \frac{{\partial }^2\phi }{\partial x^2}}+c_1^2\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\frac{E_i(k)}{1-{\nu }_i^2(k)}{z}^{6}dz}}\right\}{\displaystyle \frac{{\partial }^{3}w}{\partial x^3}}\\ =\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\rho (k)zdz}}\right\}{\displaystyle \frac{{\partial }^{2}u}{\partial t^2}}+\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\rho (k)z^{2}dz}}\right\}{\displaystyle \frac{{\partial }^2\phi }{\partial t^2}}-c_1\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\rho (k)z^{3}dz}}\right\}{\displaystyle \frac{{\partial }^{2}u}{\partial t^2}}\\ -2c_1\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\rho (k)z^{4}dz}}\right\}{\displaystyle \frac{{\partial }^2\phi }{\partial t^2}}-c_1\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\rho (k)z^{4}dz}}\right\}{\displaystyle \frac{{\partial }^{3}w}{\partial x\partial t^2}}\\ +c_1^2\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\rho (k)z^{6}dz}}\right\}{\displaystyle \frac{{\partial }^2\phi }{\partial t^2}}+c_1^2\left\{{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\rho (k)z^{6}dz}}\right\}{\displaystyle \frac{{\partial }^{3}w}{\partial x\partial t^2}},\end{array}$$

(14)

where k_w and k_P are the linear layer (spring) stiffness and shear layer stiffness, respectively, and q is the dynamic force; the rest of the coefficients are given in [10]. The foundation parameters are presented in the following dimensionless form:

$$K_{\mathrm{W}}={\displaystyle \frac{k_{\mathrm{W}}D}{L^4}}\mathrm{and}{K}_{\mathrm{P}}={\displaystyle \frac{k_{\mathrm{P}}D}{L^2}}\mathrm{for}D={\displaystyle \frac{E_{\mathrm{Cu}}{h}^3}{12(1-{\nu }_{\mathrm{Cu}}^2)}}.$$

(15)

3. Solution Procedure

To date, numerous numerical and analytical methods have been proposed for the static and dynamic analyses of different continuous systems [24,25,26,27,28]. This section details the solution procedure for the forced vibration problem introduced in Section 2.

The solution to the forced vibration problem is formulated via a set of trigonometric functions; these series are utilised to approximate the displacement components while ensuring that the boundary conditions for all simply supported edges are satisfied.

$$u(x,t)={\displaystyle \sum _{m=1}^{\infty }{U}_m\cos \frac{m\pi x}{L}}\sin (\Omega t),$$

(16)

$$w(x,t)={\displaystyle \sum _{m=1}^{\infty }{W}_m\sin \frac{m\pi x}{L}}\sin (\Omega t),$$

(17)

$$\phi (x,t)={\displaystyle \sum _{m=1}^{\infty }{\Phi }_m\cos \frac{m\pi x}{L}}\sin (\Omega t),$$

(18)

where Ω denotes the excitation frequency. The dynamic load is given by the following expression:

$$q={\displaystyle \sum _{m=1}^{\infty }{Q}_m}\sin {\displaystyle \frac{m\pi x}{L}}\sin (\Omega t),$$

(19)

where

$$Q_m=q_0,$$

(20)

with q₀ being distributed-load intensity. The utilisation of trigonometric series provides an efficient and systematic approach to representing displacement components, offering reduced computational complexity compared to methods such as the finite element or finite difference techniques. However, the methodology also has limitations. The use of trigonometric series inherently restricts their application to problems with boundary conditions that align with the assumed series. For more complex geometries or boundary conditions, alternative numerical techniques may be required.

4. Numerical Results

For validation, the natural frequencies obtained from our simulations for a simplified version of the current model (excluding the metamaterial and multilayered effects, the external force, and the elastic foundation) are compared with the results reported in Ref. [29]. As illustrated in Figure 2, the comparison demonstrates excellent agreement.

In the present study, the auxetic metamaterial beam is fabricated using a copper (Cu) matrix, characterised by a Young’s modulus E_Cu = 65.79 GPa, Poisson’s ratio ν_Cu = 0.387, and mass density ρ_Cu = 8800 kg/m³ at room temperature; graphene has material properties with a Young’s modulus E_Gr = 929.57 GPa, Poisson’s ratio ν_Gr = 0.220, and mass density ρ_Gr = 1800 kg/m³ at room temperature. The graphene layers used in this study have a length of ℓ_Gr = 8.376 nm and a thickness of t_Gr = 0.34 nm, as detailed in Refs. [23,30].

Depicted in Figure 3 are the phase plots of the beam (the centre point) subject to a harmonic force for (a) a metallic matrix (non-metamaterial), (b) a metamaterial with Pattern-U, and (c) a metamaterial with Pattern-A; each panel displays three phase plots: namely, with no foundation (K_W = K_P = 0.0), with a Winkler foundation (K_W = 10.0, K_P = 0.0), and with a Pasternak foundation (K_W = 10.0, K_P = 1.0). The results show that the addition of the elastic foundations, specifically the Winkler and Pasternak types, reduces the transverse deflection compared to the case without any foundation. This reduction occurs because the elastic foundations provide additional stiffness and resistance to deformation, effectively suppressing larger deflections. This phenomenon aligns with the theoretical predictions, where the elastic foundation parameters enhance the beam’s rigidity, thereby suppressing deflections under harmonic loading. The maximum deflection difference between the no-foundation and Winkler-foundation (panel (a)) cases is larger than those of the metamaterial ones (panels (b) and (c)). Comparing panels (b) and (c), i.e., the metamaterial cases, highlights that the effect of the Winkler foundation on the maximum deflection for Pattern-U is larger (when compared to Pattern-A).

Figure 4 depicts the time history response of a graphene origami-enabled auxetic metamaterial beam for different degrees of graphene folding (H_Gr) for a fixed graphene weight fraction (W_Gr = 2.5 wt%). For H_Gr = 20% (a), both Pattern-U and Pattern-A exhibit similar amplitudes and harmonic behaviour, indicating slight differences in the responses. At H_Gr = 50% (b), Pattern-A shows a lower amplitude compared to Pattern-U; this trend becomes more visible for H_Gr = 100%, as seen in panel (c). Furthermore, normalised transverse deflection, in general, increases with the graphene folding degree.

Figure 5 shows the phase diagrams plotted for various graphene weight fractions (W_Gr) for different values of the graphene folding degree (H_Gr): (a) H_Gr = 10%, (b) H_Gr = 30%, (c) H_Gr = 60%, and (d) H_Gr = 100%. For smaller values of the graphene folding degree (H_Gr ≤ 60%), in general, the transverse deflection decreases as the graphene weight fraction increases. However, in panel (d), where the folding degree is H_Gr = 100%, the normalised transverse deflection initially increases with W_Gr, reaching a maximum at W_Gr = 0.25 wt%, before decreasing at higher weight fractions, as seen by the closed-loop trajectories. This behaviour suggests an optimal graphene weight fraction at which the transverse deflection is maximised, likely due to the interplay between the material stiffness and flexibility, where the concentration of graphene initially enhances flexibility, leading to larger deflections, while higher concentrations lead to an increased stiffness, reducing the deflection amplitude.

Depicted in Figure 6 are the time history responses of the graphene origami-enabled auxetic metamaterial beam, illustrating the influences of the two-parameter elastic foundation coefficients on its forced dynamic behaviour. Panel (a) presents the normalised transverse deflection of the beam versus the linear layer stiffness (K_W). The results indicate that increasing the linear layer stiffness generally reduces the normalised transverse vibration amplitude. Panel (b) shows the normalised transverse deflection against the shear layer stiffness (K_P), with K_W fixed at 10. Similarly, the normalised transverse deflection decreases as the shear layer stiffness increases.

The dynamic response of an elastically constrained graphene origami-enabled auxetic metamaterial beam under different excitation frequencies is depicted in Figure 7 for both Pattern-U and Pattern-A. An increase in the excitation frequency, in general, leads to a rise in the amplitude of the oscillations. As the forcing frequency approaches the linear resonant fundamental frequency (ω₁), the transverse vibration amplitude becomes larger.

5. Conclusions

The forced vibrational behaviour of elastically constrained graphene origami-enabled auxetic metamaterial beams subjected to harmonic external forces has been investigated. The three coupled equations of motion were adopted from the literature and were modified to take into account external forcing and constraint effects. These equations were numerically solved using a set of trigonometric functions and then validated by comparing the natural frequencies with those of non-metamaterial beams available in the literature. Extensive numerical simulations were subsequently performed to analyse the forced vibration response of the metamaterial beam under harmonic excitations, focusing on the influence of the graphene origami content, distribution patterns, folding degree, and the elastic foundation linear and shear-layer stiffnesses. For the elastically constrained metamaterial beam studied in this study, the following conclusions have been reached:

• Similar to non-metamaterial beams, in general, having elastic constraints reduces the forced vibration amplitude of the metamaterial beam.
• At a fixed graphene weight fraction (W_Gr = 2.5 wt%), both Pattern-U and Pattern-A exhibit similar amplitudes with a minor difference and harmonic behaviour for H_Gr = 20%. However, for H_Gr = 50%, Pattern-A shows a lower amplitude compared to Pattern-U, and this difference becomes more dominant as H_Gr increases to 100%.
• For graphene folding degrees ≤ 60% (for the H_Gr cases studied in this paper), the transverse deflection generally decreases as the graphene weight fraction (W_Gr) increases. However, for H_Gr = 100%, the transverse vibration amplitude initially increases with W_Gr, and then starts decreasing thereafter.
• For the excitation frequencies lower than the fundamental linear frequency, the amplitude increases with the forcing frequency; however, no phase-shift is observed.