Nonlinear Resonances in 3D Printed Structures

Abstract

Nonlinear resonators can have advantages over linear designs including increased sensitivity towards changes in their physical properties and environment, and high quality factors which make them attractive in applications such as mass/chemical sensors or signal filters. Designing nonlinear structures, however, requires much understanding of nonlinear behavior characteristics of structures. Similarly, the proliferation of 3D or additive manufacturing/printing capabilities has opened the doors to deploying nonlinear resonators on scales not possible earlier. However, to obtain consistent nonlinear dynamic performance the designer must perform a careful analysis to explore the existence and repeatability of desired nonlinear behavior. Also, the use of 3D printing with the associated substrate material properties poses its own challenges in regards to device simulation in view of the fact that most of the traditional literature on nonlinear resonators assumes linear material stiffness. In this chapter, the authors discuss computational design methods for structural design, and specifically study the case of 1:2 internal resonances in resonators made of nonlinear (hyperelastic) materials. The design methods allow for development of large number of candidate resonator designs without a required significant nonlinear structural design experience, and the study of the dynamic response of the resonators provides a glimpse in to the 1:2 nonlinear internal resonance exhibited by the candidate resonators.

1. Introduction

The development of the micro- and nano-electronics industry coupled with the advances in semiconductor manufacturing techniques led to an interest in developing and applying micro- and nano-electromechanical systems (MEMS and NEMS) [1, 2]. Due to the small length scales of such devices and the popular modes of actuation employed by designers for MEMS or NEMS, such as electrostatic actuation, it led to the observations that nonlinear effects in many of these devices were the norm rather than the exception. This realization led to the study of effects of nonlinearities on the dynamic response of MEMS and NEMS in their various modes of operation, the effects in some cases being detrimental to linear design performance, and in some cases being beneficial to performance [3]. The study of nonlinear dynamics is an area of research with long history in structural and mechanical systems [4]. Several attempts have been made to incorporate the nonlinear dynamic effects and use the plethora of associated phenomena into operating mechanisms of MEMS devices particularly as mass and/or chemical sensors and filters [3]. While such nonlinear devices have several advantages over linear ones with the same functionality in terms of measurement resolution, there have been challenges in making sure that the designed topologies are “in-tune” with the semiconductor manufacturing processes. At a conceptual level, researchers in nonlinear dynamics often work with lumped-parameter models [3] and many interesting applications have been considered [5, 6] for systems characterized by a single degree of freedom. For systems with more than one degree of freedom, one of the most compact representations is a system exhibiting a nonlinear 1:2 internal resonance, the spring-pendulum system [4, 7]. While this system lends itself very conveniently to a systematic analytical study [8], it is a relatively hard task to replicate it physically at micro- or nano-scales. This is even more so when the structural components fabricated involve two- or three-dimensional elastic structures.

3D printing or additive manufacturing offers several appealing advantages in terms of building devices based on nonlinear dynamic principles. The fabrication processes and dimensions are such that there is better repeatability with complex mechanical designs, as well as initial prototyping and low volume production costs come without the need for significant capital expenditure [9, 10]. Recently studies have been reported for prototyping nonlinear vibratory components made with 3D printers having a feature size of less than 1 mm [11]. While the dimensions of polymeric resonators created by 3D printers may limit their measurement resolution and frequency range of operations, the manufacturing process itself is more repeatable for certain kinds of resonators. However, using 3D printing for producing nonlinear resonators also comes with its own challenges. While micro- and nano-resonators operate in an environment with many kinds of nonlinearities, 3D printed structures have to largely rely on only two sources of nonlinearities, namely, geometric and material nonlinearities unless controlled material variations or composite structures are explicitly introduced in fabrication. Fortunately, as was demonstrated by Tripathi and Bajaj [12, 13], both geometric nonlinearities due to finite deformations and material nonlinearities due to nonlinear hyperelastic properties of the 3D printed material can produce nonlinear dynamic effects such as 1:2 internal resonance.

A 1:2 internal resonance is a popular mechanism exhibited and employed by many nonlinear dynamics based resonators [14]. Internal resonance in a structure refers to the energy transfer that occurs between two modes of the structure when their natural frequencies are almost commensurable and the structure has some appropriate nonlinearity. For example, for 1:2 internal resonance, if the two modes of a structure have their natural frequencies close to the ratio 1:2 and harmonic excitation of the higher mode is above a certain threshold, energy can be transferred from the resonant response of the higher mode to the lower frequency mode in the presence of quadratic nonlinearities. The mathematical description of a 1:2 internal resonance and the dynamics is well established [4, 7, 8]. For the purposes of evaluating the suitability of using 3D printing to produce resonators exhibiting 1:2 internal resonances, it is important to demonstrate that the dynamic response equations for the resonators exhibit the same mathematical characteristics as the cardinal examples.

In the context of elastic structures exhibiting various internal resonances, the present work focuses on elastic plate-type structures [4, 15]. A few representative works on different aspects of nonlinear vibrations of rectangular plates with internal resonances are [16, 17, 18]. In general, isotropic plates with different simple boundary conditions do not exhibit any commensurate frequencies unless there exists some type of symmetry of the structure. Some works have considered optimization of geometry and material distribution to affect frequency distributions as well as internal resonances [19, 20]. A systematic approach is based on the concepts in “topology optimization” [21]. Applications of topology or shape optimization have now appeared in the literation on nonlinear dynamics as well, with the works in [22, 23, 24] focusing on general one-dimensional elastic systems whereas the works in [12, 13] focusing on plate structures with internal resonances. The overall goal is to tailor the system’s dynamic response to some desired form for appropriate external excitations.

In the present study, particular classes of resonator designs consisting of rectangular plates with cutouts which can be easily fabricated using 3D printing are analyzed for their nonlinear dynamic response. To obtain a suitable resonator design with commensurable (1:2) natural frequencies, a parametric optimization process which varies the sizes of the cutouts is employed. The natural frequencies themselves are computed using linear finite element analysis (FEA). The resonators are assumed to be made of a Mooney-Rivlin hyperelastic material [25] which is anticipated to provide the material nonlinearity necessary to produce 1:2 internal resonances. Once the optimization process is able to provide a candidate structure, the mode shapes obtained by the finite element analysis are used to build a reduced order model of the resonator displacements. This displacement field can then be used to derive the kinetic and strain energies of the structure which can provide the system Lagrangian. This Lagrangian is then averaged and subjected to the Euler-Lagrange conditions to derive the slow-amplitude equations of motion of the structure that provide the dynamic steady state response.

This work has two following sections: Section 2 describes the design and optimization process which leads to a desired candidate structure. It discusses the aspects of the hyperelastic material model as well as the use of mode shapes to construct the reduced order model. Section 3 describes the development of the structure’s Lagrangian, the extraction of the nonlinear equations of motion for the modal amplitudes, and the steady state dynamic response of the system under harmonic excitation of the higher frequency mode. Section 4 contains some concluding remarks for this work.

2. Candidate structure synthesis

The principal objective of the structural synthesis proves is to obtain a resonator design with commensurable natural frequencies. As it is difficult to come up with such a structure by just relying on the researcher’s experience, a computational optimization method is proposed to design the candidate resonators. For 1:2 internal resonances, the desired frequency relation between the two modes taking part in the energy transfer can be expressed as:

$$\frac{{\omega }_1}{{\omega }_2}=\frac{1}{2}$$E1

where ω₁ and ω₂ are the natural frequencies of the lower and the higher mode, respectively. To obtain such a candidate resonator, the natural frequency requirement represented by Eq. (1) can be formulated as an optimization problem. An example optimization problem for such a task can be represented by

$$\mathit{minimize},c\left(\omega \right)=\left(\frac{1}{2}-\frac{{\omega }_1}{{\omega }_2}\right)$$E2

Thus, the optimization process attempts to minimize the deviation of the two natural frequencies from the perfect 1:2 natural frequency ratio. Solving the optimization problem posed by Eq. (2) would lead to a structure with two of its natural frequencies close to the ratio of 1:2 which is a major requirement for resonators exhibiting 1:2 internal resonance. In this study, two methods for solving the optimization problem are discussed. The first method is a topology optimization method based on simple isotropic material with penalization (SIMP) model and the method of moving asymptotes (MMA) [21]. The second method is a parametric optimization method in which a starting parameterized base structure is chosen whose topology is similar to the final desired candidate structure. Then this base structure is optimized by a nonlinear quadratic programming process to produce a viable candidate structure.

2.1 Topology optimization with SIMP

Topology optimization techniques have been widely used to solve a broad range of structural optimization problems. While quite versatile, an occasional drawback against topology optimization generated optimal topologies has been the difficulty of their reproduction using conventional manufacturing processes. In this regard 3D printing is eminently suited to produce topologically optimized design as both techniques are adept at producing extruded structures with complex planar patterns. Topology optimization methods are based on finite element discretization of the design spaces. In the context of designing candidate hyperelastic resonators for 1:2 internal resonance, the design space can be discretized with finite elements and the density and material stiffness of the i^th element in the design space can be written using the SIMP formulation as

$${\mathit{\rho }}_{\mathit{i}}={\mathit{\rho }}_{\mathit{min}}+{\mathit{x}}_{\mathit{i}}^{\mathit{n}\mathbf{1}}{\mathit{\rho }}_{\mathbf{0}}$$E3

$${\mathit{E}}_{\mathit{i}}={\mathit{E}}_{\mathit{min}}+{\mathit{x}}_{\mathit{i}}^{\mathit{n}\mathbf{1}}{\mathit{E}}_{\mathbf{0}}$$E4

where, ${\mathit{\rho }}_{\mathbf{0}}$ is the material density, ${\mathit{E}}_{\mathbf{0}}$ is the material stiffness and ${\mathit{x}}_{\mathit{i}}$ is the design variable which varies between 0 and 1. ${\mathit{\rho }}_{\mathit{min}}$ and ${\mathit{E}}_{\mathit{min}}$ are infinitesimal constants to prevent numerical singularities in case ${\mathit{x}}_{\mathit{i}}$ becomes equal to zero. The exponents, n1 and n2 are usually chosen larger than three so as to penalize any intermediate values of the design variable. As can be observed from Eqs. (3) and (4), a value of ${\mathit{x}}_{\mathit{i}}=\mathbf{1}$, implies that material is present and ${\mathit{x}}_{\mathit{i}}=\mathbf{0}$ implies presence of a void. Any intermediate value of ${\mathit{x}}_{\mathit{i}}$ would produce non-physical results which is why a high value of exponents n1 and n2 are chosen to initially penalize intermediate values followed by filtering process at the end of optimization in which intermediate properties are taken to one extreme or the other depending on the filter design.

As an example of the topology optimization based resonator generation process, consider the structure shown in Figure 1. This structure is a rectangular plate which is constrained at its bottom edge. This plate is assigned Mooney-Rivlin material properties and is meshed with four node planar elements as it is assumed that the plate is undergoing vibrations in its plane.

Figure 1.

The base structure used as a starting point for the topology optimization process. The red line indicates the edge that is fixed.

The ratio of the first two planar natural frequencies of the base structure shown in Figure 1 was 3.3. The aim of the topology optimization process is to fill the central cavity of the starting structure so that its first two planar natural frequencies are in the ratio close to 1:2. Thus, the design space is discretized with finite elements and the optimization problem posed by Eq. (2) is solved using the method moving asymptotes (MMA) which yields the optimized structure shown in Figure 2.

Figure 2.

The optimized structure obtained after applying the topology optimization process to the base structure shown in Figure 1.

The ratio of the first two planar natural frequencies of the optimized structure shown in Figure 2 was 1.99. Thus, the topology optimization process was successful in bringing the natural frequencies of interest close to the ratio of 1:2. As the optimization process uses finite elements to compute the natural frequencies of the structure, the mode shapes of the optimized structure also become available and are shown in Figure 3.

Figure 3.

Mode shapes of the optimized structure shown in Figure 2. (a) Lower Mode (Mode 1 of the structure) (b) Upper Mode (Mode 2 of the structure).

2.2 Parametric optimization

Parametric optimization process is a simple but powerful tool which can also be used to generate various candidate structures for 1:2 internal resonance. As an example to illustrate the essential aspects of this procedure, consider the base structure shown in Figure 4. This base structure consists of a rectangular cantilever plate with two cutouts.

Figure 4.

The base structure used as a starting point for the optimization process in parametric optimization. The red line indicates the edge that is fixed.

This plate can be assigned Mooney-Rivlin material properties and meshed with four node shell elements. In this study, Abaqus software is used to compute the natural frequencies of the base structure with the frequencies of interest being the second and third natural frequencies respectively. For the base structure shown in Figure 4, the ratio between the natural frequencies of the higher and lower mode of interest (third and second natural frequencies, respectively) was computed as 2.4. This base structure was then subjected to an optimization process with the objective function being described by Eq. (2). The design parameters for this optimization were the cut-out size and positions on the cantilever plate and the optimization was performed by a sequential quadratic programming method. The optimized structure obtained by applying the optimization process on the base structure is shown in Figure 5.

Figure 5.

The optimized structure obtained after applying the parametric optimization process to the base structure shown in Figure 4.

The ratio between the third and second natural frequencies of the optimized structure shown in Figure 5 was 2.0. Thus, the optimization method was able to successfully bring the natural frequencies of the structure close to the desired ratio of 1:2. The mode shapes of the optimized structure also become available from the finite element model and are shown in Figure 6. These mode shapes can then be used to construct a reduced order model for the system which will be used to develop the nonlinear dynamic response for the candidate structure.

Figure 6.

Mode shapes of the optimized structure shown in Figure 5. (a) Lower Mode (Mode 2 of the structure) (b) Upper Mode (Mode 3 of the structure).

The method of parametric optimization allows for development a wide range of topologies which can each potentially exhibit 1:2 internal response. The optimal topology obtained depends on the starting structure, reflecting the local optimal nature of the solution. For example, consider the starting structure shown in Figure 7, the ratio between the natural frequencies of the higher and lower mode of interest (third and second natural frequencies, respectively) was computed as 1.6. Also note that the boundary conditions in this case involve fixing the resonator along both of its vertical sides. In this particular case, the optimization parameters were the size and location of the three circular cutouts. After performing the shape optimization process again using the sequential quadratic programming method, the optimized structure obtained is shown in Figure 8. The ratio between the third and second natural frequencies of the optimized structure shown in Figure 8 was 2.0. Thus, the examples of Figures 5 and 8 illustrate the possibilities of generating a large number of examples with different topologies as candidate resonators for 1:2 internal resonance. The mode shapes of the optimized structure shown in Figure 8 are shown in Figure 9.

Figure 7.

The base structure used as a starting point for the parameter optimization process. The red lines indicate the edges that are fixed.

Figure 8.

The optimal structure obtained after implementing the parametric optimization process on the base structure shown in Figure 7.

Figure 9.

Mode shapes of the optimized structure shown in Figure 8. (a) Lower Mode (Mode 2 of the structure) (b) Upper Mode (Mode 3 of the structure).

3. Nonlinear dynamic response

For the development of the nonlinear dynamic response of a 3D printed structure, consider the structure shown in Figure 10. This structure was designed using the simple iterative optimization procedure detailed in Section 2, and the resulting structure’s modes 2 and 3 are in near internal resonance of 1:2. Thus, the frequency ratio achieved was 2:0005. The resonator was then fabricated using 3D printing machine Stratsys Dimension 1200es. This structure has the ratio of its second and third natural frequencies as ∼2.0. The two transverse modes of interest for this candidate structure are shown in Figure 11. Using the mode shapes shown in Figure 11, assuming that the structure is subjected to a base excitation, and using the Kirchhoff plate theory [15], the displacement at any point on this rectangular plate can we written as:

$$u\left(X,Y,Z,t\right)=\epsilon \left(A_1\left(t\right)\left(u_{01}\left(X,Y\right)-Z\frac{\partial w_1\left(X,Y\right)}{\partial X}\right)+A_2\left(t\right)\left(u_{02}\left(X,Y\right)-Z\frac{\partial w_2\left(X,Y\right)}{\partial X}\right)\right)$$E5

$$v\left(X,Y,Z,t\right)=\epsilon \left(A_1\left(t\right)\left(v_{01}\left(X,Y\right)-Z\frac{\partial w_1\left(X,Y\right)}{\partial Y}\right)+A_2\left(t\right)\left(v_{02}\left(X,Y\right)-Z\frac{\partial w_2\left(X,Y\right)}{\partial Y}\right)\right)$$E6

$$w\left(X,Y,Z,t\right)=\epsilon \left(A_1\left(t\right)w_1\left(X,Y\right)+A_2\left(t\right)w_2\left(X,Y\right)\right)+w_f\left(t\right)$$E7

where, u, v and w are the displacements in the X-, Y- and Z-directions, respectively, A₁ and A₂ are the modal amplitudes, u₀₁; u₀₂ and v₀₁; v₀₂ are the independent in-plane modal displacements in the X and Y directions, and w₁ and w₂ are the corresponding modal displacements (the mode shapes) in the Z- or the transverse direction. The base excitation is applied in the transverse direction (Z-direction) and is denoted by w_f which only depends on time, and ε is a small dimensionless parameter to keep track of the significant terms in the system response. It is assumed that the displacement field can be written with linear superposition of the two modes because in the canonical 1:2 internal resonance form, the higher mode is directly excited by external means and the system nonlinearities can cause an energy transfer between the higher mode and the lower mode. All other modes, if present, will see their contribution to the displacement field decay over time in the presence of damping as they are neither directly excited, nor excited by internal energy transfer.

Figure 10.

The candidate structure for which the nonlinear dynamic transverse response is to be developed subject to a base excitation.

Figure 11.

The mode shapes of the candidate structure shown in Figure 10. The modes are (a) mode 2 and (b) mode 3 respectively of the candidate resonator.

The displacement field given in Eqs. (5)–(7) can now be used to write the kinetic and strain energy of the candidate structure. The kinetic energy, T is given by

$$T=\underset{0}{\overset{V}{\int }}\frac{1}{2}\rho \left({\dot{u}}^2+{\dot{v}}^2+{\dot{w}}^2\right)\mathit{\text{dXdYdZ}}$$E8

where the dot (.) represents the derivative of the displacements with respect to time and ρ is the material density. For the time derivatives of the displacement, it must be noted that the modal amplitudes will be differentiated and the mode shapes can be treated as constants as they do not depend on time. In a similar manner, for a Mooney-Rivlin material, the strain energy, U can be written as,

$$U=\underset{0}{\overset{V}{\int }}\left\{C_{10}\left({\overline{I}}_1-3\right)+C_{01}\left({\overline{I}}_2-3\right)+\frac{1}{d}{\left(J-1\right)}^2\right\}\mathit{\text{dXdYdZ}}$$E9

where, ${\overline{I}}_1$ and ${\overline{I}}_2$ are the first and second deviatoric invariants, respectively, of the Left Cauchy Green deformation tensor B, J is the determinant of the deformation gradient, F, and C₁₀, C₀₁ and d are the material constitutive parameters. The deformation gradient F is derived from the displacement field of the structure. The relationship between the original coordinates of a point on the plate and the deformed coordinates can be written as,

$$x=X+u$$E10

$$y=Y+v$$E11

$$z=Z+w$$E12

Then the deformation gradient F, can we written as,

$$F=\left[\begin{array}{ccc}1+\frac{\partial u}{\partial X}& \frac{\partial u}{\partial Y}& \frac{\partial u}{\partial Z}\\ \frac{\partial v}{\partial X}& 1+\frac{\partial v}{\partial Y}& \frac{\partial v}{\partial Z}\\ \frac{\partial w}{\partial X}& \frac{\partial w}{\partial Y}& 1\end{array}\right]$$E13

The left Cauchy Green deformation tensor is given by

$$B=F{F}^T$$E14

The deviatoric strain invariants of the left Cauchy Green deformation tensor B are written as

$${\overline{I}}_1=J^{-\frac{2}{3}}{I}_1$$E15

$${\overline{I}}_2=J^{-\frac{4}{3}}{I}_2$$E16

where J is the determinant of the deformation tensor given by Eq. (13), and the strain invariants I₁ and I₂ are given by

$$I_1=\mathit{tr}\left(B\right)$$E17

$$I_2=\frac{1}{2}\left(\mathit{tr}{\left(B\right)}^2-\mathit{tr}\left(B^2\right)\right)$$E18

where tr(B) refers to the trace of the matrix B. Note that using Eqs. (5)–(7) and (9)–(14), the strain energy of the structure can also be computed. Once the expressions of both the kinetic energy and the strain energy are available, the Lagrangian of the resonator can be written as

$$L=T-U$$E19

This Lagrangian from Eq. (19) will be a nonlinear function of the modal amplitudes owing to the nonlinear nature of the strain energy potential given in Eq. (9). The base excitation of the structure is now assumed to be of the form

$${\dot{w}}_f={\epsilon }^2{V}_{B}sin\left(\mathit{\Omega t}\right)$$E20

where V_B is the amplitude of the base excitation velocity and Ω is the excitation frequency. For 1:2 internal resonance, the excitation frequency can be near the lower or the upper natural frequency. In case of subharmonic external resonance, the external frequency is assumed to be close to the upper natural frequency so that the higher mode is resonantly excited [4]. The difference between the excitation frequency and the second natural frequency is known as the external mistuning σ₂ defined as

$$\Omega ={\omega }_2+\epsilon {\sigma }_2$$E21

Similarly, another mistuning parameter, the internal mistuning σ₁, is introduced to take into account the deviation of the two interacting natural frequencies from the perfect 1:2 ratio, that is,

$${\omega }_2=2{\omega }_1+\epsilon {\sigma }_1$$E22

To further study the nonlinear dynamic response of the structure for small nonlinear motions, and to formulate the application of the method of averaging [4], the modal amplitudes (for Eqs. (5)–(7)) can be written as

$${\mathit{A}}_{\mathbf{1}}\left(\mathit{t}\right)={\mathit{p}}_{\mathbf{1}}\left(\mathit{\epsilon t}\right)\mathbf{cos}\left(\frac{\mathbf{\Omega }}{\mathbf{2}}\mathbf{t}\right)+{\mathit{q}}_{\mathbf{1}}\left(\mathit{\epsilon t}\right)\mathbf{sin}\left(\frac{\mathbf{\Omega }}{\mathbf{2}}\mathbf{t}\right)$$E23

$${\mathit{A}}_{\mathbf{2}}\left(\mathit{t}\right)={\mathit{p}}_{\mathbf{2}}\left(\mathit{\epsilon t}\right)\mathbf{cos}\left(\mathbf{\Omega t}\right)+{\mathit{q}}_{\mathbf{2}}\left(\mathit{\epsilon t}\right)\mathbf{sin}\left(\mathbf{\Omega t}\right)$$E24

where p_i and q_i are amplitude components which vary on a slow time scale τ = εt, as defined in [4, 8]. Using the expressions for amplitudes in Eqs. (23) and (24), the time derivatives of the amplitudes can be written as

$${\dot{A}}_1=\epsilon \left(p_1^{\prime }\left(\mathit{\epsilon t}\right)cos\left(\frac{\mathrm{\Omega }}{2}\mathrm{t}\right)+q_1^{\prime }\left(\mathit{\epsilon t}\right)sin\left(\frac{\mathrm{\Omega }}{2}\mathrm{t}\right)\right)+\frac{\mathrm{\Omega }}{2}\left(-p_1\left(\mathit{\epsilon t}\right)sin\left(\frac{\mathrm{\Omega }}{2}\mathrm{t}\right)+q_1\left(\mathit{\epsilon t}\right)cos\left(\frac{\mathrm{\Omega }}{2}\mathrm{t}\right)\right)$$E25

$${\dot{A}}_2=\epsilon \left(p_2^{\prime }\left(\mathit{\epsilon t}\right)cos\left(\mathrm{\Omega t}\right)+q_2^{\prime }\left(\mathit{\epsilon t}\right)sin\left(\mathrm{\Omega t}\right)\right)+\frac{\mathrm{\Omega }}{2}\left(-p_2\left(\mathit{\epsilon t}\right)sin\left(\mathrm{\Omega t}\right)+q_2\left(\mathit{\epsilon t}\right)cos\left(\mathrm{\Omega t}\right)\right)$$E26

where a prime ($\prime$) denotes a derivative with respect to the slow time τ. Now the Lagrangian given in Eq. (19) is averaged over the period of oscillation$=\frac{4\pi }{\Omega }$. The slow time amplitudes are treated as constants for this averaging operation and also only terms till O(${\epsilon }^3)$ are retained in the Lagrangian as the cubic nonlinear terms are sufficient to capture the 1:2 internal resonance of the structure. The effects of internal resonance are essentially captured by quadratic nonlinear terms in the equations of motion. The averaged Lagrangian is defined by

$$⟨L⟩={\int }_0^{\frac{4\pi }{\Omega }}\left(T-U\right)\mathit{dt}$$E27

Subjecting the averaged Lagrangian shown in Eq. (27) to the Euler-Lagrange conditions ($\frac{d}{\mathit{d\tau }}\frac{\partial ⟨L⟩}{\partial p_i^{\prime }}-\frac{\partial ⟨L⟩}{\partial p_i}=0$, $\frac{d}{\mathit{d\tau }}\frac{\partial ⟨L⟩}{\partial q_i^{\prime }}-\frac{\partial ⟨L⟩}{\partial q_i}=0$, i = 1, 2) provides the following equations of motion for the slow time amplitudes

$$p_1^{\prime }+{\zeta }_1{p}_1+\left(\frac{{\sigma }_1+{\sigma }_2}{2}\right)q_1+{\Lambda }_1{\omega }_1\left(p_2{q}_1-p_1{q}_2\right)=0$$E28

$$q_1^{\prime }+{\zeta }_1{q}_1-\left(\frac{{\sigma }_1+{\sigma }_2}{2}\right)p_1+{\Lambda }_1{\omega }_1\left(p_1{p}_2+q_1{q}_2\right)=0$$E29

$$p_2^{\prime }+{\zeta }_1{p}_2+\left({\sigma }_2\right)q_2-2{\Lambda }_2{\omega }_2{p}_1{q}_1=0$$E30

$$q_2^{\prime }+{\zeta }_1{q}_2-\left({\sigma }_2\right)p_2+{\Lambda }_2{\omega }_2\left(p_1^2-q_1^2\right)={\Lambda }_3{V}_B$$E31

where a prime ($\prime$) denotes a derivative with respect to the slow time τ, and.

Λ_i, i = 1, 2, 3, are constants which come from the averaged Lagrangian of the structure and depend on the material constitutive parameters (Eq. (9)) and mode shapes. The modal damping terms ζ₁ and ζ₂ were introduced in Eqs. (28)–(31) to represent damping in the system to prevent the solutions from becoming unbounded. The expressions in Eqs. (28)–(31) are the same as the expressions in a standard 1:2 internal resonance system [4, 8], thus demonstrating that 3D printed rectangular plates with cutouts are able to exhibit nonlinear dynamic phenomena such as 1:2 internal resonance provided the constants Λ_i, i = 1, 2, 3, are non-zero. To make the analysis of Eqs. (28)–(31) a little more tractable, the following variable transformations are defined

$$p_1=a_{1}cos\left({\beta }_1\right),q_1=a_1\mathrm{s}\text{in}\left({\beta }_1\right)$$E32

$$p_2=a_{2}cos\left({\beta }_2\right),q_2=a_2\mathrm{s}\text{in}\left({\beta }_2\right)$$E33

where $a_i$’s are the amplitudes and ${\beta }_i$ are the phase angles. With the transformations introduced in Eqs. (32) and (33), the equations of motion for modal amplitudes become

$$a_1^{\prime }=-{\zeta }_1{a}_1-{\Lambda }_1{\omega }_1{a}_1{a}_2\mathit{sin}\left(2{\beta }_1-{\beta }_2\right)$$E34

$$a_1{\beta }_1^{\prime }=\left(\frac{{\sigma }_1+{\sigma }_2}{2}\right)a_1-{\Lambda }_1{\omega }_1{a}_1{a}_2\mathit{cos}\left(2{\beta }_1-{\beta }_2\right)$$E35

$$a_2^{\prime }=-{\zeta }_2{a}_2-{\Lambda }_2{\omega }_2{a}_1^2\mathit{sin}\left({\beta }_2-2{\beta }_1\right)+{\Lambda }_3{V}_{B}sin\left({\beta }_2\right)$$E36

$$a_2{\beta }_2^{\prime }={\sigma }_2{a}_2-{\Lambda }_2{\omega }_2{a}_1^2\mathit{cos}\left(2{\beta }_1-{\beta }_2\right)+{\Lambda }_3{V}_{B}cos\left({\beta }_2\right)$$E37

The steady-state solutions for the system of equations described by Eqs. (34)–(37) can be obtained by setting $a_i^{\prime }=0$and ${\beta }_i^{\prime }=0$. These equations can be solved for steady-state solutions to give single-mode (only second modal amplitude a₂ is non-zero) and coupled-mode solutions (both first and second mode amplitudes are non-zero, that is, a₁ ≠ 0 and a₂ ≠ 0). The coupled-mode solution is the main nonlinear response as it implies energy transfer from the higher to lower mode when the higher mode is directly excited in the presence of quadratic nonlinearities. Plots in Figure 12 give a representative set of steady-state solutions for the single and coupled-mode response for the structure shown in Figure 10 with zero damping. Note that the coupled-mode response is slightly asymmetric about ${\sigma }_2=0$ axis for the structure as some minimal internal mistuning exists (${\sigma }_1\ne 0)$ due to the fact that $\frac{{\omega }_2}{{\omega }_1}=$ 2:0005. For perfect internal resonance, the coupled-mode solutions will be completely symmetric around ${\sigma }_2=0$. The non-zero first mode arises as a result of the subcritical pitchfork bifurcations (at P₁ and P₂) from the single mode solution consisting of only the second mode (The lower mode amplitude is zero). A more detailed study of the stability of the solutions for the single and coupled mode results is provided in [13].

Figure 12.

Non-linear response of the 3D printed structure shown in Figure 10 to a transverse harmonic base excitation. The plots are for the amplitudes of the two interacting modes for both the single-mode and coupled-mode response. Note that ${\sigma }_1\ne \mathit{0}$ though very small and the modal damping is low (ζ₁ and ζ₂ equal to 0.05).

As is clear from Eqs. (34)–(37), the modal amplitudes depend upon many parameters. Some of the more interesting of these are the internal mistuning σ₁ and the modal damping terms ${\zeta }_1$ and ${\zeta }_2$. The effect of change in internal mistuning is shown in Figure 13. As can be observed from Figure 13, changes in internal mistuning can result in the coupled-mode motion to lose existence and disappear, that is, the modal interaction is lost for large internal mistuning. Thus, in actual physical systems deviation of natural frequencies of participating modes from the perfect 1:2 ratio can cause the non-trivial coupled mode response to not manifest itself.

Figure 13.

Non-linear response curves of the 3D printed (hyperelastic) structure shown in Figure 10 for (a) large negative internal mistuning, and (b) large positive internal mistuning.

Figure 14 shows the effect of damping coefficients on the nonlinear response curves obtained using Eqs. (34)–(37). Increasing damping coefficients ${\zeta }_1$ and ${\zeta }_2$ has interesting effects on the overall nonlinear response. As can be observed from Figure 14a, increasing damping leads to a reduction in the frequency range in which the nonlinear coupled-mode response is observed.

Figure 14.

Non-linear response curves of the hyperelastic structure in Figure 10 for representative low (red) and higher (blue) damping coefficients. The two figures are for (a) mode 1 amplitudes (b) mode 2 amplitudes. The points T1 and T2 are turning point bifurcations.

4. Summary and conclusions

This work explored the possibility of synthesizing 3D printed hyperelastic plate structure exhibiting 1:2 internal resonances. 3D printing occupies a potential sweet spot in terms of dimensional capabilities and repeatability to produce nonlinear resonators which can be used as vibration absorbers, sensors, or for signal processing applications. The synthesis methodology allows for designing a large set of designs meeting the desired internal resonance conditions resulting in complex modal coupling and energy transfer between modes of the structure.

While the nonlinear dynamical response studied here was focused on 3D printed cantilever plates that exhibited 1:2 internal resonances on account of material nonlinearities, the methodology can be easily applied to other boundary conditions and internal resonances, as well as for structures with geometric nonlinearities caused by finite deformations of plates.

Conflict of interest

The authors declare no conflict of interest.