An accurate approximate analysis of the large-amplitude asymmetric vibrations of a dielectric elastomer balloon

Abstract

The vibration model of the neo-Hookean DEB has a negative index nonlinearity which poses a solution problem for perturbative methods and most non-perturbative methods. In this paper, an accurate solution of the large-amplitude asymmetric vibrations of a dielectric elastomer balloon (DEB) made up of neo-Hookean material was derived using the continuous piecewise linearization method (CPLM), which was shown to handle the negative index nonlinearity smoothly. The accuracy of the CPLM was verified using published analytical results which showed that the CPLM was significantly more accurate than the published analytical results. The CPLM solution was used to conduct parametric analysis of the large-amplitude vibration of the DEB. The results revealed that increasing either of the applied voltage or static internal pressure results to lesser stiffness response, lower frequency of vibration, lower maximum kinetic energy and shorter range of initial stretch values for stable vibrations. Hence, the applied voltage and static internal pressure can be used as design parameters for neo-Hookean DEBs capable of stable large-amplitude vibrations.

Introduction

The exploration of electricity’s potential to deform solid materials has a rich history dating back to the late 18th century when scientists were experimenting with overcharged Leyden jars, which are devices that store electric charge. These jars were discovered by Georg von Kleist and Pieter van Musschenbroek of Leiden. The late 19th century saw Wilhelm Conrad Röntgen, who discovered X-rays, invent the first dielectric elastomer by coating a rubber sheet with metal electrodes. Röntgen used a natural rubber elastomer stripe, with 16cm width and 100 cm length, which he pre-stretched by a weight to approximately twice its initial length. He then charged the two surfaces of the elastomer with opposite charges by means of a corona discharge. Attractive forces between the opposite electric charges caused the rubber band to thin down in thickness and to expand in length resulting in a 10% increase of the initial length. ¹ It wasn’t until the late 20th century that dielectric elastomers began to garner significant attention as smart materials with potential applications as sensors, actuators, and generators. Pelrine et al. ² reported that films of dielectric elastomers, such as silicones, can be used to create electrical actuators. These films were coated with a flexible electrode material on both sides. When voltage was applied, the film was compressed and expanded due to electrostatic forces, resulting in strains of up to 30%–40%. It was found that pre-straining the film could enhance their performance. With silicone elastomers, actuated strains of up to 117% were achieved, and with acrylic elastomers, strains of up to 215% were achieved using both bi-axially and uni-axially pre-strained films.

In the field of electroactive polymer research, much study has been devoted to dielectric elastomers (DE). Dielectric elastomer actuators are polymer films sandwiched between two compliant electrodes. When a voltage difference is applied on the compliant electrodes, it causes a compression in thickness of the polymer film and a stretching in its area. Dielectric elastomer actuators possess several desirable qualities such as high energy density, low power requirement, large strain displacements and fast response time, making them suitable for a wide range of applications at the micro- and macro-scale. Applications of DE include soft robotics, synthetic muscle-like actuators, haptic devices, reciprocating or peristaltic pumps, loudspeakers, tunable diffraction grating for telecommunication and display devices, acoustic actuators, and adaptive optical devices.^3–5

In the last two decades research efforts to understand the quasi-static and dynamic response of dielectric elastomers have increased rapidly. Key to this endeavour is the accurate prediction of the constitutive behaviour of elastomers, which are inherently hyperelastic. The most common models used to predict the deformation of elastomers are the Ogden family of $n$ -parameter hyperelastic material models (Sharma et al. 2018), ⁴ which include the single parameter neo-Hookean model, the two parameter Mooney-Rivlin model and the two parameter Ogden model. Further advancements in the field include the three parameter models for the large-deformation of elastomers presented by Davidson and Goulbourne, ⁶ and the work of Mangan and Destrade ⁷ that examined some phenomenological model of hyperelasticity and concluded that the three parameter Gent-Gent model is accurate and versatile for elastomeric balloon inflation. Itskov et al. ⁸ proposed a new electromechanical constitutive model for dielectric elastomers based on the molecular chain statistics of elastomers. The constitutive model is an extended form of the classical neo-Hookean model and was shown to give good prediction of experimental data reported by Carpi et al. ⁹ and Kollosche et al. ¹⁰ Anssari-Benam and Bucchi ¹¹ developed a generalized neo-Hookean model with two model parameters that provided improved fits to experimental data compared to the Gent model. The generalized neo-Hookean model captures strain-stiffening effects just as well as the Gent model. However, when strain-stiffening effect is negligible, the classical one parameter neo-Hookean model is the simplest accurate hyperelastic model available and is commonly used for the dynamic analysis of DEB inflation. More details on the modelling progress of the constitutive behaviour of dielectric elastomer structures can be obtained from the recent review of Lu et al. ¹²

Zhu et al. ¹³ investigated the dynamic behaviour of a dielectric elastomer balloon (DEB) based on the neo-Hookean hyperelastic model. They found that the DEB resonates at multiple values of the frequency of excitation giving rise to super-harmonic, harmonic, and sub-harmonic responses. They also predicted the existence of jump and hysteresis in the amplitude due to continuous variation of the excitation frequency. Yong et al. ¹⁴ investigated the dynamics of a thick-walled dielectric elastomer spherical shell and derived an explicit analytical equation of motion for the spherical shell based on simple geometrical and spherical capacitor assumptions. Liang and Cai ¹⁵ investigated the shape bifurcation of a spherical DEB under the actions of internal pressure and electric voltage. Using linear perturbation analysis, they demonstrated that a spherical DEB can bifurcate to a non-spherical DEB under certain electromechanical loading conditions. From the analysis of the variation of the free energy of the DEB, they showed that the DEB can be energetically stable or unstable depending on the electromechanical loading conditions. Chen et al. ¹⁶ studied the dynamic electromechanical instability of a DEB subject to parametric voltage excitations. They observed that the DEB experienced dynamic electromechanical instabilities over a wide range of excitation frequencies, in the form of snap-through or snap-back, provided a critical voltage was exceeded. Jin and Huang ¹⁷ proposed a stochastic differential equation (SDE) to examine the random response of a DEB. They derived an analytical solution based on stochastic averaging technique and the results of the analytical solution agree with results of Monte Carlo simulation. Wang et al. ¹⁸ presented a theoretical model to investigate the effect of strain-stiffening on the nonlinear vibration of a circular DE membrane subjected to electromechanical loading. They applied the Gent model for elastomers and investigated the free vibration of the circular DE membrane under parametric perturbations due to applied sinusoidally alternating voltage. They solved the vibration model numerically and predicted the existence of bifurcations, bi-stable vibrations and chaotic vibrations. Liang and Cai ¹⁹ developed a finite element model to study the occurrence a new instability mode that can cause localized failure. The new instability results to a localized bulging-out of the DEB, which is consistent with the experimental results in Li et al. ²⁰ Lv et al. ²¹ proposed a theoretical model for a balloon actuator taking into account the effects of strain stiffening and damping force. They investigated the occurrence of snap-through or snap-back, bi-stable vibrations and chaotic vibrations in the presence of viscous damping. Liu and Zhou ²² investigated the nonlinear oscillation of viscoelastic dielectric elastomers (DEs) subject to alternating voltage. They derived the vibration model based on Lagrange method, neo-Hookean constitutive behaviour and non-equilibrium thermodynamics. The vibration model was solved using a combined shooting and arc-length continuation numerical method. Their results showed the existence of bifurcation, jump phenomenon, super-harmonic resonance, softening behaviour and hysteresis. Mao et al. ²³ exploited the biasing fields to tune the three-dimensional small-amplitude free vibrations of a soft thick-walled electroactive spherical balloon. They applied both the neo-Hookean and Gent electroelastic models and showed that their analytical solutions account for radial breathing, torsional and any other spheroidal modes of vibration. Alibakhshi and Heidari ²⁴ applied the Gent-Gent hyperelastic model to examine the nonlinear vibration of a DEB with strain-stiffening and second invariant of the Cauchy-Green deformation tensor. They identified the chaotic interval for critical parameters of the DEB and determined that the second invariant parameter could suppress the chaotic motion of the DEB. Alibakhshi et al. ²⁵ used the generalized neo-Hookean model of Anssari-Benam and Bucchi ¹¹ to investigate the nonlinear vibration and stability of a spherical DEB. They determined that the molecular structure of the DEB material can influence the critical voltage for pull-in or snap through instability, chaos and quasi-periodicity. They concluded that the DEB tends to be more stable and less prone to chaos when strain stiffening occurs early in the deformation of the DEB.

The studies on the nonlinear vibrations of DEBs discussed so far were conducted using numerical methods. In contrast, only few studies have been conducted using approximate analytical techniques. The main reason is that the hyperelastic model of the DEB has a negative power nonlinearity (or non-polynomial nonlinearity) that is not amenable to most approximate analytical methods, especially those methods that rely on asymptotic series expansion. However, some successful attempts have been made to use approximate analytical techniques because of the benefit of providing additional insight of the vibration response that may not be obvious from the numerical solution. Tang et al. ²⁶ used the Newton Harmonic Balance (NHB) method to estimate the dynamic response of a freely vibrating DEB with neo-Hookean hyperelastic model. They compared the NHB results with numerical results and found that the second-order NHB gave accurate results for a stretch ratio of up to 2.5. Beyond this stretch ratio, they recommended the use of higher order approximations that are algebraically tedious to derive. Follow-up studies on the free vibration of the DEB in the presence of strain-stiffening effect were conducted in Tang et al. ²⁷ The NHB was used to examine the frequency-amplitude response, which showed that the frequency reduces as the pressure or voltage increases. They also investigated the stability of the DEB and the analysis revealed the presence of a stable centre and a saddle node. Alibakhshi and Heidari ²⁸ applied the multiple scales method (MMS) to estimate the vibration response of a DEB modelled as a neo-Hookean material. They used the Taylors’ series expansion to transform the vibration model to a nonlinear ODE with quadratic and cubic nonlinearities, which enabled solution using MMS. They predicted the occurrence of jump phenomenon and softening behaviour in the DEB. Wang et al. ²⁹ investigated a viscously damped parametrically excited DEB vibration model using the incremental harmonic balance (IHB) method. They applied the neo-Hookean model to describe the deformation of the DEB and combined the IHB with the arc-length continuation technique to obtain periodic solutions of the DEB. Their results showed an improved accuracy over the classical harmonic balance method.

The existing approximate analytical solutions are limited to weakly or moderately nonlinear responses and are not able to capture accurately the strong nonlinear asymmetric oscillations inherent in the large-amplitude vibrations of the DEB. ²⁶ Estimating the large-amplitude vibrations of the DEB is important for predicting the vibration responses occurring close to the snap-through failure conditions. Under these conditions, strong nonlinear effects are at play due to hyperelastic effects. This underscores the need for the derivation of more accurate approximate analytical solutions capable of capturing the strong nonlinear hyperelastic effects during large-amplitude vibrations. In this article, the continuous piecewise linearization method (CPLM) was applied to obtain accurate periodic solutions for a large-amplitude DEB with neo-Hookean deformation behaviour. The CPLM was chosen because of its simplicity and ability of produce very accurate results irrespective of the complex nature of the restoring force of the vibrating system. The derived CPLM solution was validated by comparing with results of standard numerical solutions and the NHB. Furthermore, a stability analysis was conducted to examine the stable and unstable equilibrium points, and then, a parametric analysis was conducted to analyse the effects of the static internal pressure and DC voltage on the large-amplitude vibration response of the DEB. The remainder of the article is arranged as follows. In Section ‘System description and model derivation’, a brief description of the DEB and derivation of the vibration model was presented. In section ‘Solution methodology’, the CPLM solution algorithm was presented and applied to solve the DEB model whereas the results were discussed in Section ‘Discussion of results’. The conclusions of the study and areas for further improvement were discussed in Section ‘Conclusions’.

System description and model derivation

The dynamic inflation of dielectric elastomers with neo-Hookean deformation behaviour and having different shapes has been extensively studied. Figure 1 shows a spherical balloon of radius $R$ and thickness $H$ in the reference state. The balloon membrane is made of a dielectric elastomer of density $\rho$ that is incompressible and undergoing dynamic inflation in an isothermal environment. Both sides of the membrane are coated with compliant electrode materials and the DEB is assumed to be subject to static pressure and voltage. In the presence of an applied voltage, the electrostatic force causes the membrane to reduce in thickness while the balloon expands in volume and surface area. ¹³

Figure 1.

The state of a spherical DEB before and after applying static pressure and DC voltage.

In the presence of an applied pressure ( $P$ ) and an electric potential difference with electrode charges $+Q$ and $-Q$ , the balloon deforms homogenously from an initial radius $R$ to a final radius $r$ . Given that the deformation of the DEB is homogeneous, the stretch is given by:

$$\lambda =\frac{r}{R}$$

(1)

On the other hand, the electric displacement is defined as the charge per unit area of the spherical DEB and is given by the equation:

$$D=\frac{Q}{4\pi r^2}$$

(2)

The work done on the system consists of the electric work due to electric charge variation, the isobaric work due to variation in the radius at constant pressure, and the inertia work due to variation in the radius of the DEB mass. When the charge on the electrode has a small variation, $\delta Q$ , the voltage ( $\mathrm{\Phi }$ ) does a work that is given as $\mathrm{\Phi }\delta Q$ , Also, the pressure ( $P$ ) does a work due to a small variation in the radius ( $\delta r$ ) that is given as $4\pi r^{2}P\delta r$ , and simultaneously the inertia force does a work that is given as $-4\pi R^{2}H\rho (d^{2}r/d{t}^2)\delta r$ . Due to energy conservation of the system, the change in the energy of the system is equal to the work done on the system. Hence,

$$4\pi R^{2}H\delta W=\mathrm{\Phi }\delta Q+4\pi r^{2}P\delta r-4\pi R^{2}H\rho \frac{d^{2}r}{d{t}^2}\delta r$$

(3)

where $\delta W$ is a small variation in the energy density (or free energy function). Using equations (1) to (3), we have:

$$\delta W=\frac{\mathrm{\Phi }}{H}{\lambda }^2\delta D+(\frac{2\mathrm{\Phi }D\lambda }{H}+\frac{\mathrm{PR}}{H}{\lambda }^2-\rho R^2\frac{d^2\lambda }{d{t}^2})\delta \lambda$$

(4)

From thermodynamics considerations, ¹³ the free energy of system is determined by two independent variables, the stretch ( $\lambda$ ) and the electric displacement ( $D$ ). Therefore, the small variation of free energy can be expressed as:

$$\delta W(\lambda ,D)=\frac{\partial W(\lambda ,D)}{\partial D}\delta D+\frac{\partial W(\lambda ,D)}{\partial \lambda }\delta \lambda$$

(5)

From equation (4) and (5) the following partial derivatives are applicable:

$$\frac{\partial W\left(\lambda ,D\right)}{\partial \lambda }=\frac{2\mathrm{\Phi }D\lambda }{H}+\frac{PR}{H}{\lambda }^2-\rho R^2\frac{d^2\lambda }{d{t}^2}$$

(6a)

$$\frac{\partial W(\lambda ,D)}{\partial D}=\frac{\mathrm{\Phi }}{H}{\lambda }^2$$

(6b)

Equation (6a) represents the force balance while equation (6b) represents the electrostatic balance.

The free energy function of the dielectric elastomer based on the neo-Hookean constitutive behaviour is given as ²⁷ :

$$W(\lambda ,D)=\frac{\mu }{2}(2{\lambda }^2+{\lambda }^{-4}-3)+\frac{D^2}{2\epsilon }$$

(7)

where $\mu$ is the shear modulus and ε is the permittivity, which is a constant that is independent of the deformation. The first part of equation (7) in bracket represents the elastic energy and the second part is the dielectric energy. Substituting equation (7) into equations (6a) and (6b), and eliminating $D$ gives:

$$\frac{d^2\lambda }{d{\tau }^2}+f\left(\lambda \right)=0$$

(8)

where $f\left(\lambda \right)$ is the conservative restoring force given as:

$$f\left(\lambda \right)=2\lambda -2{\lambda }^{-5}-\frac{\mathrm{PR}}{\mu H}{\lambda }^2-2\frac{\epsilon {\mathrm{\Phi }}^2}{\mu H^2}{\lambda }^3$$

(9)

$\tau =t/R\sqrt{\rho /\mu }$ is the dimensionless time, $\mathrm{PR}/\mu H$ is the dimensionless pressure and $\epsilon {\mathrm{\Phi }}^2/\mu H^2$ is the dimensionless electric potential difference. The initial conditions to equation (8) are $\lambda \left(0\right)=A>0$ and $\lambda \prime \left(0\right)=0$ where the prime denotes differentiation with respect to $\tau$ and $A$ is the initial stretch.

The restoring force of the DEB is a function of three independent variables, that is $f(\lambda ,P,\mathrm{\Phi })$ , where the pressure and voltage are generally time-dependent functions. Nonetheless, the pressure and voltage are constant when the restoring force is purely conservative and the restoring force becomes a function of stretch only, that is $f\left(\lambda \right)$ . It is obvious that equation (9) is a mixed parity function which implies that the DEB undergoes asymmetric vibrations in general. Equation (9) is a complex restoring force with negative power nonlinearity. The latter makes equation (8) intractable to conventional perturbation methods, except if an approximate restoring force is derived in terms of a polynomial series with positive indices. ²⁸ Such series approximation of the restoring force introduces a first level error even before applying the approximate analytical method. In contrast, the CPLM algorithm requires no approximation of the restoring force and can work with non-polynomial functions of any complexity. Hence, we demonstrate how the CPLM can be applied to solved equations (8) and (9) in Section ‘Solution Methodology’.

Solution methodology

CPLM solution for neo-Hookean DEB vibration model

The CPLM is a non-perturbative iterative approximate analytical technique for the periodic solution of nonlinear conservative systems. The methodology of the CPLM was initially proposed by Big-Alabo et al. ³⁰ in a similar algorithm called the force-indentation linearisation method (FILM). Nevertheless, the application of the methodology to derive the periodic response of nonlinear conservative oscillators, which is now called the CPLM, was first carried out by Big-Alabo. ³¹ The CPLM is based on the concept of piecewise discretization and linearization of the nonlinear restoring force. The linearized restoring force results to a linearized equation of motion whose analytical solution gives an approximate solution of the nonlinear model over the range of discretization. The iteration of the CPLM procedure automatically updates the solution of the linearized equation from one discretization to the next, thereby producing an unconditionally-stable solution for any nonlinear conservative oscillator. ³² The CPLM has been used to study the asymmetric vibrations of a ship roll motion ³³ and a functionally graded microbeam ³⁴ , but it has not been applied for an asymmetric oscillator with negative power nonlinearity as in the present case.

For conservative systems undergoing asymmetric vibrations, the discretization of the nonlinear restoring force is done over a half cycle, that is, ( $B<\lambda <A$ ), where $A$ and $B$ are the amplitudes on either side of the centre of vibration. In the case of the DEB vibrations, $|A-{\lambda }_c|\ne |{\lambda }_c-B|$ where ${\lambda }_c$ is the centre of vibration, which can be determined from the positive root of the nonlinear equation:

$${\lambda }^8+\frac{\mathrm{PRH}}{2\epsilon {\mathrm{\Phi }}^2}{\lambda }^7-\frac{\mu H^2}{\epsilon {\mathrm{\Phi }}^2}{\lambda }^6+\frac{\mu H^2}{\epsilon {\mathrm{\Phi }}^2}=0$$

(10)

At the amplitudes, the stretch rate $(\lambda \prime )$ is zero. The amplitude $A$ is the initial stretch while the other amplitude $B$ is the positive root of the nonlinear equation obtained from the condition, $\lambda \prime =0$ . Therefore,

$$B^8+\frac{2\mathrm{PRH}}{3\epsilon {\mathrm{\Phi }}^2}{B}^7-\frac{2\mu H^2}{\epsilon {\mathrm{\Phi }}^2}{B}^6-\frac{\mu H^2}{\epsilon {\mathrm{\Phi }}^2}=0$$

(11)

Typical variations of the asymmetric amplitude ( $B$ ) and the asymmetric ratio ( $\sigma$ ) with the initial stretch ( $A$ ) are shown in Figure 2 for $1.10\le A\le 2.91$ . The asymmetric ratio defines the asymmetry in the vibrations of the DEB and it is given as: $\sigma =(A-{\lambda }_c)/(A-B)$ . Its value ranges from 0.50 to 1.0 with $\sigma =0.50$ representing a symmetric vibration and $\sigma >0.50$ representing an asymmetric vibration. As $\sigma \to 1.0$ , the vibrations become more asymmetric. Figure 2 reveals that the centre of vibration is independent of the initial stretch and the asymmetric amplitude approaches an asymptotic limit as the initial stretch approaches the limiting stretch from which point the DEB experiences electromechanical instability.

Figure 2.

Variation of asymmetric amplitude and asymmetric ratio with initial stretch: $a=\mathrm{PR}/\mu H=0.1$ and $b=\epsilon {\mathrm{\Phi }}^2/\mu H^2=0.1$ .

The roots of equation (10) give the equilibrium points ( ${\lambda }_e$ ) of the DEB, which may be stable or unstable. The stability condition of the equilibrium points can be determined from the gradient of the conservative restoring force given as:

$$\frac{\mathrm{df}}{d\mathit{\lambda }}=g\left(\lambda \right)=2+10{\lambda }^{-6}-2\frac{\mathrm{PR}}{\mu H}\lambda -6\frac{\epsilon {\mathrm{\Phi }}^2}{\mu H^2}{\lambda }^2$$

(12)

The stability condition is such that if $g\left({\lambda }_e\right)>0$ , then ${\lambda }_e$ is stable and ${\lambda }_e={\lambda }_c$ is a centre of vibration. Conversely, if $g\left({\lambda }_e\right)\le 0$ then ${\lambda }_e$ is unstable. The unstable points represent the limiting stretch which is the maximum value of the initial stretch below which stable vibrations occurs and from which no vibrations occur. Figure 3 shows 3-D stability plots of the DEB for various values of the non-dimensional static pressure ( $a=\mathrm{PR}/\mu H$ ) and the non-dimensional electric potential ( $b=\epsilon {\mathrm{\Phi }}^2/\mu H^2$ ). The plots reveal that the neo-Hookean DEB has only one stable equilibrium point for a given ( $a$ , $b$ ). Figures 3(a) and (b) illustrate the regions of stable and unstable points respectively, whereas Figure 3(c) shows that the neo-Hookean DEB undergoes a saddle-node bifurcation resulting in a transition from stable equilibrium to unstable equilibrium and vice versa.

Figure 3.

Stability analysis of neo-Hookean DEB vibration model: (a) stable equilibrium, (b) unstable equilibrium, (c) saddle-node bifurcation.

Conventionally, the discretization procedure of the CPLM is based on $n$ equal divisions of the solution domain. ³² For the DEB vibrations, the discretization range is a half-cycle and the step size for each discretization is given as: $\mathrm{\Delta }\lambda =|A-B|/n$ . Accordingly, the linearized force within each discretization could be expressed as ³⁴ :

$$f_{\mathrm{rs}}\left(\lambda \right)=\pm \left|K_{\mathrm{rs}}\right|(\lambda -{\lambda }_r)+F_r$$

(13)

where r and $s=r+1$ represent the initial and the end points of each discretization respectively, $F_r=f\left({\lambda }_r\right)$ and $K_{\mathrm{rs}}$ is the linearized stiffness given as:

$$\begin{array}{c}{K}_{\mathrm{rs}}=\frac{f\left({\lambda }_s\right)-f({\lambda }_r)}{{\lambda }_s-{\lambda }_r}\hfill \\ =\frac{2({\lambda }_s-{\lambda }_r)-({\lambda }_s^{-5}-{\lambda }_r^{-5})-\frac{\mathrm{PR}}{\mu H}({\lambda }_s^2-{\lambda }_r^2)-2\frac{\epsilon {\mathrm{\Phi }}^2}{\mu H^2}({\lambda }_s^3-{\lambda }_r^3)}{{\lambda }_s-{\lambda }_r}\hfill \end{array}$$

(14)

Substituting eq. (13) in (9) gives the linearized equation of motion for each discretization as:

$$\frac{d^2\lambda }{d{\tau }^2}\pm \left|K_{\mathrm{rs}}\right|\lambda =\pm \left|K_{\mathrm{rs}}\right|{\lambda }_r-F_r$$

(15)

The solution to equation (15) depends on whether the sign before $\left|K_{\mathrm{rs}}\right|$ is positive or negative. If $K_{\mathrm{rs}}>0$ then the sign is positive while $K_{\mathrm{rs}}<0$ implies that the sign is negative.

Positive linearized stiffness, $K_{\mathrm{rs}}>0$

The solution to equation (15) when $K_{\mathrm{rs}}>0$ is given as:

$$\lambda \left(\tau \right)=R_{\mathrm{rs}}\sin ({\omega }_{\mathrm{rs}}\tau +{\phi }_{\mathrm{rs}})+C_{\mathrm{rs}}$$

(16)

where ${\omega }_{\mathrm{rs}}=\sqrt{\left|K_{\mathrm{rs}}\right|}$ ; $C_{\mathrm{rs}}={\lambda }_r-F_r/\left|K_{\mathrm{rs}}\right|$ ; and $R_{\mathrm{rs}}=\sqrt{{({\lambda }_r-C_{\mathrm{rs}})}^2+{({\lambda }_r\prime /{\omega }_{\mathrm{rs}})}^2}$ . The initial conditions and other parameters are determined based on the oscillation stage. For the first half-cycle that moves from $A$ to $B$ the initial conditions are ${\lambda }_r=A-r\mathrm{\Delta }\lambda$ and

$$\begin{array}{c}\multicolumn{1}{c}{{\lambda }_r\prime \left(0\right)=-\sqrt{|2\underset{A}{\overset{{\lambda }_r}{\int }}-f\left(\lambda \right)d\lambda |}}\\ =-\sqrt{|2(A^2-{\lambda }_r^2)+(A^{-4}-{\lambda }_r^{-4})-\frac{2\mathrm{PR}}{3\mu H}(A^3-{\lambda }_r^3)-\frac{\epsilon {\mathrm{\Phi }}^2}{\mu H^2}(A^4-{\lambda }_r^4)|}\hfill \end{array}$$

The other solution parameters for this half-cycle are:

$${\phi }_{\mathrm{rs}}=\{\begin{array}{c}\pi +\mathrm{ta}{n}^{-1}\left[\frac{{\omega }_{\mathrm{rs}}({\lambda }_r-C_{\mathrm{rs}})}{{\lambda }_r\prime }\right],\mathrm{if}{\lambda }_r\prime <0\hfill \\ \pi /2,\mathrm{if}{\lambda }_r\prime =0\hfill \end{array}$$

(17)

$$\mathrm{\Delta }\tau =\{\begin{array}{c}\frac{1}{{\omega }_{\mathrm{rs}}}[0.5\pi +\mathrm{co}{s}^{-1}\left(\frac{{\lambda }_s-C_{\mathrm{rs}}}{R_{\mathrm{rs}}}\right)-{\phi }_{\mathrm{rs}}],\mathrm{if}({\lambda }_s-C_{\mathrm{rs}})<R_{\mathrm{rs}}\hfill \\ \frac{1}{{\omega }_{\mathrm{rs}}}(0.5\pi -{\phi }_{\mathrm{rs}}),\mathrm{if}({\lambda }_s-C_{\mathrm{rs}})\ge R_{\mathrm{rs}}\hfill \end{array}$$

(18)

where $\mathrm{\Delta }\tau$ is the time range for each discretization. The time at the end of each discretization is ${\tau }_s={\tau }_r+\mathrm{\Delta }\tau$ and the end conditions ${\lambda }_s$ and ${\lambda }_s\prime$ are determined by replacing $r$ with $s$ in the formulas for the initial conditions.

Basically, the CPLM derives the corresponding time solutions for points on the closed loop of the exact state-space solution. In other words, the CPLM generates the solution array ( ${\tau }_s$ , ${\lambda }_s$ , ${\lambda }_s\prime$ ) at the endpoint of each discretization and this can be used to produce the displacement and velocity histories. In this sense, the CPLM can be implemented as a point-wise solution. However, it is noteworthy that the CPLM solution for each discretization is applicable over a time range and not just a point. The point-wise and range-wise attributes of the CPLM and its unique qualities that distinguish it from conventional numerical techniques have been discussed by Big-Alabo and Ossia. ³²

On the other hand, the initial conditions for the second half-cycle that moves from $B$ to $A$ are ${\lambda }_r=B+r\mathrm{\Delta }\lambda$ and

$$\begin{array}{c}\multicolumn{1}{c}{{\lambda }_r\prime =+\sqrt{|2\underset{A}{\overset{{\lambda }_r}{\int }}-f\left(\lambda \right)d\lambda |}}\\ =\sqrt{|2(A^2-{\lambda }_r^2)+(A^{-4}-{\lambda }_r^{-4})-\frac{2\mathrm{PR}}{3\mu H}(A^3-{\lambda }_r^3)-\frac{\epsilon {\mathrm{\Phi }}^2}{\mu H^2}(A^4-{\lambda }_r^4)|}\hfill \end{array}$$

(19)

The other solution parameters for this half-cycle are:

$${\phi }_{\mathrm{rs}}=\{\begin{array}{c}\mathrm{ta}{n}^{-1}\left[\frac{{\omega }_{\mathrm{rs}}({\lambda }_r-C_{\mathrm{rs}})}{{\lambda }_r\prime }\right],\mathrm{if}{\lambda }_r\prime >0\hfill \\ -\pi /2,\mathrm{if}{\lambda }_r\prime =0\hfill \end{array}$$

(20)

$$\mathrm{\Delta }\tau =\{\begin{array}{c}\frac{1}{{\omega }_{\mathrm{rs}}}[0.5\pi -\mathrm{co}{s}^{-1}\left(\frac{{\lambda }_s-C_{\mathrm{rs}}}{R_{\mathrm{rs}}}\right)-{\phi }_{\mathrm{rs}}],\mathrm{if}({\lambda }_s-C_{\mathrm{rs}})<R_{\mathrm{rs}}\hfill \\ \frac{1}{{\omega }_{\mathrm{rs}}}(0.5\pi -{\phi }_{\mathrm{rs}}),\mathrm{if}({\lambda }_s-C_{\mathrm{rs}})\ge R_{\mathrm{rs}}\hfill \end{array}$$

(21)

Negative linearized stiffness, $K_{\mathrm{rs}}<0$

The solution to equation (15) when $K_{\mathrm{rs}}<0$ is given as:

$$\lambda \left(\tau \right)=A_{\mathrm{rs}}{e}^{{\omega }_{\mathrm{rs}}\tau }+B_{\mathrm{rs}}{e}^{-{\omega }_{\mathrm{rs}}\tau }+C_{\mathrm{rs}}$$

(22)

where ${\omega }_{\mathrm{rs}}=\sqrt{\left|K_{\mathrm{rs}}\right|}$ ; $C_{\mathrm{rs}}={\lambda }_r+F_r/\left|K_{\mathrm{rs}}\right|$ ; $A_{\mathrm{rs}}$ and $B_{\mathrm{rs}}$ are integration constants that were derived based on the initial conditions (same as in $K_{\mathrm{rs}}>0$ ) as:

$$A_{\mathrm{rs}}=\frac{1}{2}({\lambda }_r+\frac{{\lambda }_r\prime }{{\omega }_{\mathrm{rs}}}-C_{\mathrm{rs}})$$

(23a)

$$B_{\mathrm{rs}}=\frac{1}{2}({\lambda }_r-\frac{{\lambda }_r\prime }{{\omega }_{\mathrm{rs}}}-C_{\mathrm{rs}})$$

(23b)

The time range for each discretization was obtained by applying the end conditions to the solution in equation (22). Thus,

$$\begin{array}{c}\mathrm{\Delta }\tau =\hfill \\ \{\begin{array}{c}\frac{1}{{\omega }_{\mathrm{rs}}}\mathrm{In}\left[\frac{({\lambda }_s-C_{\mathrm{rs}})+\mathrm{sgn}({\lambda }_r\prime )\sqrt{{({\lambda }_s-C_{\mathrm{rs}})}^2-4A_{\mathrm{rs}}{B}_{\mathrm{rs}}}}{2A_{\mathrm{rs}}}\right],{({\lambda }_s-C_{\mathrm{rs}})}^2>4A_{\mathrm{rs}}{B}_{\mathrm{rs}}\hfill \\ \frac{1}{{\omega }_{\mathrm{rs}}}\mathrm{In}\left[\frac{({\lambda }_s-C_{\mathrm{rs}})}{2A_{\mathrm{rs}}}\right],{({\lambda }_s-C_{\mathrm{rs}})}^2\le 4A_{\mathrm{rs}}{B}_{\mathrm{rs}}\hfill \end{array}\hfill \end{array}$$

(24)

where $\mathrm{sgn}\left(x\right)$ is the signum function.

Finally, the nonlinear frequency of the asymmetric vibrations of the DEB can be determined as:

$${\omega }_{\mathrm{CPLM}}=\frac{2\pi }{T_{\mathrm{CPLM}}}=\frac{\pi }{{\sum }_{i=1}^n\mathrm{\Delta }{\tau }_i}$$

(25)

The CPLM procedure discussed above was developed into a pseudocode algorithm (see Appendix) which was implemented as a bespoke Mathematica code to simulate results.

Evaluation of the strength of the nonlinearity of the DEB vibration model

Nonlinear oscillators exhibit various degrees of nonlinearity depending on the values of the stiffness parameters and the amplitude of vibration. The degree of nonlinearity is important because it can inform the choice of the approximate analytical solution that should be applied. The reason is that many approximate analytical solutions do not give accurate prediction of the vibration response for stiff or strong nonlinear systems. The degree of nonlinearity is usually categorized as weak, moderate or strong nonlinearity. For most of the literature, such distinctions have been made arbitrarily and usually on the basis of what is considered a small- or large-amplitude vibration. Big-Alabo ³⁵ proposed a simple energy-based criterion to objectively evaluate the degree of nonlinearity of a nonlinear oscillator based on its amplitude of vibration and stiffness parameters. The criterion is based on comparing the potential energy of the nonlinear oscillator with that of an equivalent linear oscillator. The relative error between these two energy values was applied to determine the degree of nonlinearity.

The asymmetric vibration of the DEB has displacement limits at $A$ and $B$ where $A>{\lambda }_c>B>0$ . Since $f\left({\lambda }_c\right)=0$ , then $f\left(A\right)>0$ and $f\left(B\right)<0$ . The implication is that the area under the compliance curve in the stretch range, ${\lambda }_c\le \lambda \le A$ , is positive and the area in the stretch range, $B\le \lambda \le {\lambda }_c$ , is negative. Therefore, the total potential energy during a half-cycle vibration (i.e., $B\le \lambda \le A$ ) is:

$$V_{\mathrm{AB}}=V\left(A\right)+V\left(B\right)-2V\left({\lambda }_c\right)=2[V\left(A\right)-V\left({\lambda }_c\right)]$$

(26)

where

$$V\left(\lambda \right)={\lambda }^2+\frac{1}{2}{\lambda }^{-4}-\frac{\mathrm{PR}}{3\mu H}{\lambda }^3-\frac{\epsilon {\mathrm{\Phi }}^2}{2\mu H^2}{\lambda }^4$$

(27)

The equivalent linear oscillator can be written as:

$$\frac{d^2\lambda }{d{\tau }^2}+K_{\mathrm{eq}}\lambda =K_{\mathrm{eq}}{\lambda }_c$$

(28)

where $K_{\mathrm{eq}}$ is the equivalent linear stiffness that approximates the asymmetric nonlinear stiffness. Considering that the DEB vibrations are characterized by a non-zero centre of vibration ( ${\lambda }_c$ ) and asymmetric vibrations, $K_{\mathrm{eq}}$ can be expressed as a bi-linear stiffness as follows ³⁵ :

$$\begin{array}{c}\multicolumn{1}{c}{K_{\mathrm{eq}}=\{\begin{array}{c}\frac{f\left(A\right)-f\left({\lambda }_c\right)}{A-{\lambda }_c}\mathrm{for}{\lambda }_c\le \lambda \le A\hfill \\ \multicolumn{1}{c}{\frac{f\left({\lambda }_c\right)-f\left(B\right)}{{\lambda }_c-B}\mathrm{for}B\le \lambda \le {\lambda }_c}\end{array}}\end{array}$$

(29)

Since $f\left({\lambda }_c\right)=0$ , we have:

$$\begin{array}{c}\multicolumn{1}{c}{K_{\mathrm{eq}}=\{\begin{array}{c}\frac{f\left(A\right)}{(A-{\lambda }_c)}\mathrm{for}{\lambda }_c\le \lambda \le A\hfill \\ \multicolumn{1}{c}{-\frac{f\left(B\right)}{({\lambda }_c-B)}\mathrm{for}B\le \lambda \le {\lambda }_c}\end{array}}\end{array}$$

(30)

Therefore, the potential energy of the equivalent bi-linear oscillator is:

$$\begin{array}{c}\multicolumn{1}{c}{V_{\mathrm{eqi}}=\{\begin{array}{c}\frac{f\left(A\right)(A-{\lambda }_c)}{2}\mathrm{for}{\lambda }_c\le \lambda \le A\hfill \\ \multicolumn{1}{c}{\frac{f\left(B\right)(B-{\lambda }_c)}{2}\mathrm{for}B\le \lambda \le {\lambda }_c}\end{array}}\end{array}$$

(31)

where $i=1$ for ${\lambda }_c\le \lambda \le A$ and $i=2$ for $B\le \lambda \le {\lambda }_c$ . Figure 5 shows the potential energies of the linear approximations ( $V_{\mathrm{eq}1}$ and $V_{\mathrm{eq}2}$ ) in comparison with the actual potential energy, $V\left(A\right)$ , in each quarter. The difference, $V\left(A\right)-V_{\mathrm{eqi}}$ , gives the error in using a linear approximation for each quarter cycle. Therefore, the relative error of the bi-linear approximation was determined as the arithmetic mean of the relative errors for each quarter cycle as follows:

$${ϵ}_r=\frac{1}{2}(|1-\frac{f\left(A\right)(A-{\lambda }_c)}{2[V\left(A\right)-V\left({\lambda }_c\right)]}|+|1-\frac{f\left(B\right)(B-{\lambda }_c)}{2[V\left(A\right)-V\left({\lambda }_c\right)]}|)\times 100\%$$

(32)

According to the energy-based criterion ³⁵ : ${ϵ}_r\le 0.25\%$ implies that the system is linear, $0.25\%<{ϵ}_r\le 2.5\%$ implies that the system is weakly nonlinear, $2.5\%<{ϵ}_r\le 25\%$ implies that the system is moderately nonlinear and ${ϵ}_r>25\%$ implies that the system is strongly nonlinear.

Typical plots of the compliance response and the potential function of the neo-Hookean DEB vibration model are shown in Figure 4. These plots show the strong nonlinear and asymmetric nature of the DEB vibrations. Figure 5 plots the relative error to evaluate the nonlinear strength of the DEB vibration model based on equation (32). If $a=\mathrm{PR}/\mu H=0.1$ and $b=\epsilon {\mathrm{\Phi }}^2/\mu H^2=0.1$ , then the DEB experiences a moderate nonlinearity for initial stretch values in the range $1.10\le A<1.40$ and a strong nonlinearity when the initial stretch is $A>1.40$ . As $A\to 2.91$ , the nonlinearity becomes stronger and the periodic solution is more difficult to obtain accurately using approximate analytical methods. ²⁶ However, this article demonstrates that the CPLM solution presented in Section ‘CPLM solution for neo-Hookean DEB vibration model’ is one exception that can produce accurate periodic solutions for the large-amplitude vibrations of the DEB.

Figure 4.

(a) Compliance response and (b) potential function of the neo-Hookean DEB vibration model: $a=0.1$ and $b=0.1$ .

Figure 5.

Estimating the nonlinear strength of the neo-Hookean DEB vibration model: $a=0.1$ and $b=0.1$ . $V_{\mathrm{ab}}=$ exact PE during a quarter cycle (see equation (26)), $V_{\mathrm{eq}1}=$ linear approximation of PE for ${\lambda }_c\le \lambda \le A$ , $V_{\mathrm{eq}2}=$ linear approximation of PE for $B\le \lambda \le {\lambda }_c$ .

Discussion of results

In this section, the results of the present study were discussed. In the first instance, the accuracy of the derived CPLM solution was tested using conventional numerical solutions (Runge-Kutta method) as a reference solution. Also, comparison with the NHB was conducted. Then, the effects of the static pressure and voltage on the large-amplitude vibration were discussed. In all cases, the input parameters used were $a=\mathrm{PR}/\mu H=0.1$ and $b=\epsilon {\mathrm{\Phi }}^2/\mu H^2=0.1$ except where it was stated otherwise.

Verification of accuracy of CPLM solution

The convergence test was carried out to determine the optimum number of discretization that produces less than 0.1% relative error in the CPLM time period estimate during large-amplitude vibrations. The results of the test for $A=2.8$ are plotted in Figure 6. The CPLM converges to an error of 0.1% at $n=125$ and 0.073% at $n=150$ . Hence, we used $n=150$ for all subsequent investigations. To test the accuracy of using $n=150$ to predict the vibration response of the DEB, a comparison was made with the second-order Newton Harmonic Balance (NHB) method ²⁶ and numerical solution as shown in Figures 7 to 10. In all cases, the CPLM was several orders more accurate than the NHB. While the NHB performs well for small- and moderate-amplitude vibrations (see Figures 7 and 8), it produces significant error for large- and very large-amplitude vibrations (see Figures 9 and 10).

Figure 6.

Convergence of the CPLM time period solution.

Figure 7.

Comparison of CPLM and NHB for small-amplitude ( $A=1.10$ ) vibration of the DEB: (a) Displacement history, (b) Phase plot, (c) Error analysis of the first cycle displacement history.

Figure 8.

Comparison of CPLM and NHB for moderate-amplitude ( $A=1.50$ ) vibration of the DEB: (a) Displacement history, (b) Phase plot, (c) Error analysis of the first cycle displacement history.

Figure 9.

Comparison of CPLM and NHB for large-amplitude ( $A=2.50$ ) vibration of the DEB: (a) Displacement history, (b) Phase plot, (c) Error analysis of the first cycle displacement history.

Figure 10.

Comparison of CPLM and NHB for very large-amplitude ( $A=2.80$ ) vibration of the DEB: (a) Displacement history, (b) Phase plot, (c) Error analysis of the first cycle displacement history.

It can be observed that the error in the CPLM displacement solution in Figures 8 to 10 peaks around the asymmetric amplitude. When $\lambda \left(t\right)<1.0$ , the relative error in the CPLM solution becomes sensitive because small deviations in the absolute error show up as significant relative errors due to division by a number less than 1.0. The closer $\lambda \left(t\right)$ gets to zero the more sensitive is the error in the CPLM solution. This explanation accounts for the peak relative error of the CPLM in Figures 8 to 10 occurring around the asymmetric amplitude.

Parametric analysis of the vibration response of the neo-Hookean DEB

The effect of the critical parameters on the large-amplitude vibration response of the neo-Hookean DEB was examined. Specifically, the effects of the static internal pressure and DC voltage on the large-amplitude vibration response of the DEB were investigated. For this purpose, an initial stretch of 2.5 was used and the results of the parametric analysis are presented in Figures 11 to 13. In the absence of applied electric voltage, the DEB material experiences softening and becomes less stiff as the static internal pressure increases. This results to a lower frequency of vibration (Figure 11(a)) and a lower maximum kinetic energy (Figure 11(a)). In the case of a constant applied voltage (Figure 12), the same behaviour was observed regarding the influence of internal pressure on the large-amplitude vibration response, except that the observed behaviour is more pronounced. Similarly, increasing the applied voltage in the absence of internal pressure leads to more softening effects in the large-amplitude vibration response, whereas the introduction of a constant pressure makes the softening effects more obvious (Figure 13).

Figure 11.

Influence of static pressure on the large-amplitude ( $A=2.50$ ) vibration response of the DEB with no applied electric voltage: (a) Displacement history, (b) Velocity history and (c) Phase plot.

Figure 12.

Influence of static pressure on the large-amplitude ( $A=2.50$ ) vibration response of the DEB with constant applied electric voltage: (a) Displacement history, (b) Velocity history and (c) Phase plot.

Figure 13.

Influence of applied electric voltage on the large-amplitude ( $A=2.50$ ) vibration response of the DEB with constant static pressure: (a) Displacement history, (b) Velocity history and (c) Phase plot.

Furthermore, it was observed that increasing the pressure and voltage reduces the range of values of the initial stretch that permits stable vibrations. This is because the limiting stretch shifts towards the centre of vibration as the pressure and voltage increases. For instance, if $a=0.1$ and $b=0.1$ the centre of vibration is 1.029 and the initial stretch range for stable vibrations is $1.029<A<2.920$ , whereas $a=0.3$ and $b=0.1$ gives 1.054 and $1.054<A<2.494$ respectively, and $a=0.3$ and $b=0.3$ gives 1.151 and $1.151<A<1.517$ respectively. The direction of shift of the limiting stretch causes the DEB to experience stable vibrations within a more restricted range of initial stretch. At the point (i.e., values of $a$ and $b$ ) where the centre of vibration collides with the limiting stretch a saddle node bifurcation occurs and the DEB transits from stable vibrations to instability. The implication is that the applied voltage and static internal pressure can be used as control parameters for the design of stable DEB systems capable of large-amplitude vibrations.

Conclusions

The study applied the continuous piecewise linearization method to examine the large-amplitude vibration of a dielectric elastomer balloon made up of neo-Hookean material. The vibration model of the DEB is highly nonlinear and is characterized by a negative index nonlinearity. Consequently, most approximate analytical methods are either unable to solve the neo-Hookean DEB vibration model directly or they produce approximate solutions that are only accurate for small- to moderate-amplitude vibrations. However, the CPLM can solve vibration models with complex nonlinearity and it was applied to handle the negative index nonlinearity. The CPLM solution that was derived for the neo-Hookean DEB vibration model produced very accurate results that are several orders more accurate than the corresponding NHB results for small- to large-amplitude vibrations. Therefore, the CPLM is recommended for the vibration analysis of other DEB models, for example, models accounting for strain-stiffening effect during large-amplitude vibrations.

A qualitative treatment of the neo-Hookean DEB model was conducted and the stability analysis showed that the stable vibrations of the neo-Hookean DEB are all uni-centric. Also, the DEB experiences a saddle node bifurcation that results in transition from stable vibrations to instability and vice versa. An analysis of the strength of nonlinearity of the DEB showed that the DEB exhibits strong nonlinear vibration for most of the initial stretch that produces stable vibrations. The nonlinearity becomes stronger as the initial stretch approaches the limiting stretch thereby producing a more asymmetric vibration response. The asymmetry in the vibration response can be defined in terms of the ratio of the deviation of the initial stretch from the centre of vibration to the difference between both amplitudes. This asymmetric ratio ranges from 0.5 to 1.0. A value of 0.5 represents a symmetric response while a value greater than 0.5 represents an asymmetric response. The closer the ratio is to 1.0, the more asymmetric the response and the stronger the nonlinearity.

The neo-Hookean DEB vibration model has two critical parameters which are the applied electric voltage and the static internal pressure. Investigations to examine the effects of these two parameters on the vibration response of the DEB revealed that an increase in either or both of these parameters reduces the effective stiffness of the DEB thereby producing a lower frequency response with lower maximum kinetic energy. In addition, the limiting initial stretch moves closer to the centre of vibration as either of the two parameters are increased. This results to a more restricted range of initial stretch values that can yield stable vibrations. Hence, the two parameters can be used as control parameters for the design of stable large-amplitude vibrations of the neo-Hookean DEB.

This study provides a methodological framework that can be used to conduct future studies to accurately investigate the stable dynamic response of DEB vibrations. Additionally, more accurate non-polynomial hyperelastic models incorporating strain-stiffening effects and non-local constitutive behaviour can be applied to conduct vibration analysis using the CPLM.