Jtam-A4.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

50, 3, pp. 743-753, Warsaw 2012
50th Anniversary of JTAM

NONLINEAR DYNAMIC ANALYSIS OF VISCOELASTIC MEMBRANES

DESCRIBED WITH FRACTIONAL DIFFERENTIAL MODELS

John T. Katsikadelis

School of Civil Engineering, National Technical University of Athens, Athens, Greece

e-mail: jkats@central.ntua.gr

The dynamic response of an initially flat viscoelastic membrane is investigated. The visco-
elastic model is described with fractional order derivatives. The membrane is subjected to
surface transverse and inplane dynamic loads. The governing equations are three coupled
second order nonlinear partial FDEs (fractional differential equations) of hyperbolic type in
terms of the displacement components. These equations are solved using the BEM for frac-
tional partial differential equations developed recently by Katsikadelis. Without excluding
other viscoelastic models, the herein employed material is the Kelvin-Voigt model with a
fractional order derivative. Numerical examples are presented which not only demonstrate
the efficiency of the solution procedure, but also give a better insight into this complicated
but very interesting response of structural viscoelastic membranes. It is worth noting that
in case of resonance, phenomena similar to those of the Duffing equation are observed.

Key words: viscoelastic membranes, nonlinear vibrations, nonlinear fractional differential
equations, analog equationmethod, boundary elements

1. Introduction

Membranes made of linear viscoelastic materials are extensively used in modern engineering
applications. Various models of the integer differential form have been proposed in order to
describe themechanical behavior of suchmaterials (e.g.Maxwell,Voigt,Kelvin,Zener).Recently,
many researchers have shown that viscoelastic models of the differential form with fractional
derivatives are in better agreement with the experimental results than the integer derivative
models (Stiassnie, 1979; Adolfsson et al., 2005; Meral et al., 2010).

The dynamic response of pure elastic and viscoelastic membranes, using differential consti-
tutive equations of an integer order derivative or hereditary integral type models, have been
examined by many investigators (Plaut, 1990; Koivurova and Pramila, 1997; Jevine and Mac-
kintosh, 2002;Wineman, 2007; Goncales et al., 2009). The use of fractional differentialmodels is
very limited and is restricted to viscoelastic beams (Galucio et al., 2004) and three dimensional
viscoelastic bodies (Schmidt andGaul, 2002). However, viscoelasticmembranes of the fractional
type derivative have not been analyzed to date. The reason is plausible as the response of such
membranes is described by a system of nonlinear partial FDEs for which an analytical solution
is out of question, while numerical methods were not available until lately. Recently, Katsika-
delis developed the AEM, which in conjunction with an integral equation approach provides an
efficient computational tool for solving linear and nonlinear FDEs, ordinary (Katsikadelis, 2009)
and partial (Katsikadelis, 2011). This method is general and paves the way to analyze com-
plicated systems whose response is described by such equations. It has been already employed
to solve several problems. Among them, the linear fractional diffusion wave equation in boun-
ded inhomogeneous anisotropic bodies (Katsikadelis, 2008), the response of the inhomogeneous
anisotropic viscoelastic bodies (Nerantzaki and Babouskos, 2011) and the nonlinear dynamic



744 J.T. Katsikadelis

response of viscoelastic plates (Babouskos and Katsikadelis, 2009; Katsikadelis and Babouskos,
2010).

Without excluding other fractional differential models (Nerantzaki and Babouskos, 2011),
the employed herein viscoelastic material is described by the Kelvin-Voigt type model with a
fractional order derivative











σx
σy
τxy











=
E

1−ν2









1 ν 0
ν 1 0

0 0
1−ν
2























εx+ηD
a
cεx

εy +ηD
a
cεy

γxy+ηD
a
cγxy















(1.1)

where E, ν are the elastic material constants, η the viscoelastic parameter and Dac the Caputo
fractional derivative of the order α defined as

Dαcu(t)=























1

Γ(m−α)

t
∫

0

u(m)(τ)

(t− τ)α+1−m
dτ for m−1<α<m

dm

dtm
u(t) for m=α

(1.2)

where m is a positive integer. The advantage of this definition is that it permits the assignment
of initial conditions which have direct physical significance (Podlubny, 1999). Apparently, the
classical derivatives result fromEq. (1.2) for integer values of α. Thus, Eq. (1.1) for α=1 gives
the constitutive equation for the conventional Kelvin-Voigt model, while for α = 0 the pure
elastic model with the elastic modulus E∗ =E(1+η). The use of fractional differential models,
besides their simplicity to formulate the equations of structural viscoelastic systems, exhibit a
major advantage, namely they can describe more complicated integer order differential models
with fewparameters. Nevertheless, due to lack of efficient numericalmethods to solve the FDEs,
the use of such models, despite their advantages, was very limited to date. The development of
efficient numerical methods brings the fractional derivative models in the center of the arena.

The spatial discretization leads to a system of fractional Duffing type equations, which
include quadratic terms, beside the cubic ones. Therefore, similar resonance phenomena due
to a time harmonic excitation are observed. The nonlinear resonance of the membrane under
a harmonic force is studied by integrating the equations of motion over a long time, until
a steady-state solution is obtained. The stable solutions (attractors) can be obtained by an
appropriate selection of the initial conditions. However, this methodmay need a very long time
to reach a stable solution, while it is impossible to obtain the unstable ones. More efficient
methods, such as the harmonic balance method, the Poincare map and arc-length methods
for periodic solutions (Amabili, 2008) have been employed to obtain numerical results. The
dominant resonanceappearswhenthe excitation frequency is close to the lower natural frequency
of the membrane structure. Secondary resonances occur for specific values of the excitation
frequency, which may be lower (superharmonic resonance) or greater (subharmonic resonance)
than the first natural frequency. In all cases, the amplitude of vibration increases, while jump
and hysteretic phenomena appear due to the nonlinear character of the problem. Such behavior
hasbeen reportedpreviously inmany single andmulti-degree of freedomnonlinear systemsusing
analytical and approximate techniques (Nayfeh and Mook, 1979). A membrane of an arbitrary
shape is analyzed. Free vibrations as well as forced vibrations under step and harmonic loads
are studied. The obtained numerical results validate the efficiency of the solution method and
allow one to draw useful conclusions for the dynamic response of viscoelastic membranes.



Nonlinear dynamic analysis of viscoelastic membranes... 745

2. Problem statement and governing equations

We consider a thin flexible initially flat elastic membrane of thickness h and mass density ρ
consisting of a homogeneous linearly viscoelastic material occupying the two-dimensional, in
general multiply connected, domain Ω in the x,y plane. Themembrane is prestressed either by
the imposed displacement u∗n, v

∗

t or by external forces N
∗

n, N
∗

t acting along the boundary Γ .
Moderate large deflections result from nonlinear kinematic relations, which retain the square of
the slopes of the deflection surface, while the strain components remain still small compared
with the unity. This theory is good for considerably large deflections with the exception in the
vicinity of the boundary where the stress resultants of the finite deformation of the membrane
should be considered to satisfy the equilibrium and explain the folding near the edge. Thus, the
strain components are given as

εx =u,x+
1

2
w,2x εy = v,y+

1

2
w,2y γxy =u,y+v,x+w,xw,y (2.1)

where u = u(x,y,t), v = v(x,y,t) are the inplane (membrane) displacement components
and w = w(x,y,t) the transverse one. The membrane is subjected to the transverse load
pz = pz(x,y,t) and the inplane loads px = px(x,y,t) and py = py(x,y,t).
On the basis of Eqs. (1.1) and (2.1) the stress resultants are written as

N̄x =Nx+D
α
cNx N̄y =Ny +D

α
cNy N̄xy =Nxy +D

α
cNxy (2.2)

where

Nx =C
[

u,x+νv,y+
1

2
(w,2x+νw,

2
y )
]

Ny =C
[

νu,x+v,y+
1

2
(νw,2x+w,

2
y )
]

Nxy =C
1−ν
2
(u,y+v,x+w,xw,y )

(2.3)

C =Eh/(1−ν2) is the membrane stiffness.
The governing equations result by taking the equilibrium of the membrane element in a

slightly deformed configuration. This yields

N̄x,x+N̄xy,y−ρhü=−px N̄y,y+N̄xy,x−ρhv̈=−py
N̄xw,xx+2N̄xyw,xy+N̄yw,yy−pxw,x−pyw,y−ρhẅ+ρhüw,x+ρhv̈w,y=−pz

(2.4)

Without restricting the generality, we consider displacement boundary conditions on Γ

un =u
∗

n ut = v
∗

t w=w
∗ (2.5)

Inserting Eqs. (2.2) combined with Eqs. (2.3) in Eqs. (2.4), we obtain equations governing the
response of the membrane in terms of the displacements in Ω

C
{[

u,x+νv,y+
1

2
(w,2x+νw,

2
y )
]

+ηDαc

[

u,x+νv,y+
1

2
(w,2x+νw,

2
y )
]}

,x

+C
1−ν
2
{(u,y+v,x+w,xw,y )+ηDαc (u,y+v,x+w,xw,y )},y−ρhü=−px

C
1−ν
2
{(u,y+v,x+w,xw,y )+ηDαc (u,y+v,x+w,xw,y )},x

+C
{[

νu,x+v,y+
1

2
(νw,2x+w,

2
y )
]

+ηDαc

[

νu,x+v,y+
1

2
(νw,2x+w,

2
y )
]}

,y−ρhv̈=−py
(2.6)

C
{[

u,x+νv,y+
1

2
(w,2x+νw,

2
y )
]

+ηDαc

[

u,x+νv,y+
1

2
(w,2x+νw,

2
y )
]}

w,xx

+C(1−ν){(u,y+v,x+w,xw,y )+ηDαc (u,y+v,x+w,xw,y )}w,xy
+C
{[

νu,x+v,y+
1

2
(νw,2x+w,

2
y )
]

+ηDαc

[

νu,x+v,y+
1

2
(νw,2x+w,

2
y )
]}

w,yy

−pxw,x−pyw,y−ρhẅ+ρhüw,x+ρhv̈w,y=−pz



746 J.T. Katsikadelis

Equations (2.6) are subjected besides boundary conditions (2.5) also to the initial conditions

u(x,y,0)= f1(x,y) u̇(x,y,0)= g1(x,y)

v(x,y,0)= f2(x,y) v̇(x,y,0)= g2(x,y)

w(x,y,0)= f3(x,y) ẇ(x,y,0)= g3(x,y)

(2.7)

Equations (2.6) constitute a system of three coupled nonlinear partial FDEs of hyperbolic
type that are solved using the BEM presented in Katsikadelis (2011), which for the problem at
hand is adjusted as in the following.

3. The BEM solution

Equations (2.6) are of the second order with regard to the spatial derivatives, thus according to
the principle of the analog equation they can be replaced by the equations

∇2u= b(1)(x, t) ∇2v= b(2)(x, t) ∇2w= b(3)(x, t) (3.1)

where x= {x,y}∈Ω, t> 0 and b(i)(x, t), i=1,2,3 represent timedependentfictitious sources,
unknown in the first instance. The solution to Eq. (3.1)1 is given in integral form (Katsikadelis,
2002)

εu(x, t) =

∫

Ω

u∗b dΩ−
∫

Γ

(u∗q−q∗u) ds x∈Ω∪Γ (3.2)

in which q = u,n; u
∗ = ℓnr/2π is the fundamental solution to Eq. (3.1) and q∗ = u,∗n its

derivative normal to the boundary with r= ‖ξ−x‖, x∈Ω∪Γ and ξ ∈Γ ; ε is the free term
coefficient (ε=1 if x∈Ω, ε=ω/2π if x∈Γ and ε=0 if x /∈Ω∪Γ) and ω is the interior angle
between the tangents of the boundary at the point x; ε=1/2 for points where the boundary is
smooth.Eq. (3.2) is solvednumericallyusing theBEM.Theboundary integrals areapproximated
using N constant boundary elements, whereas the domain integrals are approximated using M
linear triangular cells. The domain discretization is performedautomatically using theDelaunay
triangulation. Since the fictitious source is not defined on the boundary, the nodal points of the
triangles adjacent to the boundary are placed on their sides (Fig. 1). Thus, after discretization

Fig. 1. Boundary and domain discretization

and application Eq. (3.2) at the N boundary nodal points, we obtain

Hu=Gq+Ab(1) (3.3)

where H,G are N×N known coefficientmatrices originating from the integration of the kernel
functions on theboundaryand A is an N×M coefficientmatrix originating fromthe integration
of the kernel function on the domain cells; u, q are the vectors of the N nodal displacements



Nonlinear dynamic analysis of viscoelastic membranes... 747

and slopes and b(1) the nodal values of the fictitious source at the M domain nodal points,
all at time t. The derivatives of u(x, t) inside Ω are obtained by direct differentiation of Eq.
(3.2) with ε=1. Further, by applying Eq. (3.2) and its derivatives at the domain nodal points
andmaking use of Eq. (3.3) combined with the boundary conditions to eliminate the boundary
quantities, we can express u(t) and its derivatives at the domain nodal points in terms of the
fictitious source

u,pq (t)=S,pqb
(1)(t)+s(1),pq p,q=0,x,y (3.4)

Similarly, we obtain

v,pq (t)=S,pqb
(2)(t)+s(2),pq p,q=0,x,y (3.5)

and

w,pq (t)=S,pqb
(3)(t)+s(3),pq p,q=0,x,y (3.6)

where S,pq are known M×M matrices and s(i),pq known vectors. For homogeneous boundary
conditions, we have s(i),pq=0. The notation u,pq designates u,00=u,u,x0=u,x etc. Applying
now Eqs. (2.6) at the M domain nodal points and substituting the involved derivatives from
Eqs. (3.4)-(3.6), we obtain

M
(i)
b̈
(1)+F(i)(b(i),Dαcb

(i))= f(i) i=1,2,3 (3.7)

where M(i) are M ×M consistent mass matrices, F(i) are M × 1 vectors whose elements
depend nonlinearly on the arguments and f(i) are known load vectors. Using Eqs. (3.4)-(3.6),
initial conditions (2.7) for b(i) become

b
(i)(0)=S−1(fi−s(i)) ḃ(i)(0)=S−1gi (3.8)

Equations (3.7) constitute a system of 3M three-term nonlinear FDEs, which are solved using
the time step numerical procedure developed by Katsikadelis (2009). The use, however, of all
degrees of freedom, i.e. 3M, may be computationally costly and in some cases inefficient due

to the large number of the time dependent variables b
(i)
k
(t). To overcome this difficulty in this

investigation, thenumber of degrees of freedom is reducedusing theRitz transformation, namely

b(i) =Ψ(i)z(i) (3.9)

where z
(i)
k
(t), (k =1, . . . ,L<M) are new time dependent parameters and Ψ(i) is the M ×L

transformationmatrix. In this investigation the eigenmodes of the linear problem are employed
as Ritz vectors.

4. Numerical examples

4.1. Free and forced vibrations under suddenly applied uniform load

The dynamic response of the viscoelastic membrane shown in Fig. 2 has been studied. The
boundary of the domain is defined by the curve

r=
√
ab 4

√

√

√

√

√

√

√

(

cosθ
b

)2
+
(

sinθ
a

)2

(

cosθ
a

)2
+
(

sinθ
b

)2
0¬ θ¬ 2π



748 J.T. Katsikadelis

Fig. 2. Boundary and domain nodal points of the membrane

The membrane is prestressed by un = 0.2m in the direction normal to the boundary and
ut = 0m in the tangential direction. The employed data are: a = 3, b = 1.3; h = 0.002m,
ρ/h = 5000kg/m3, E = 1.1 · 105kN/m2, ν = 0.3. The results were obtained using N = 210
boundary elements and M =106 internal collocation points resulting from 164 linear triangular
domain cells.

Firstly, the influence of the inplane inertia forces is investigated. Themembrane is assumed
purely elastic and is subjected to a transverse load pz = 2H(t)kN/m

2. Figures 3a,b,c present
the time history of the transverse displacement w(0,0, t) at the center of the membrane, the
inplane displacement u(−1.741,0, t) and the membrane force Nx(−1.741,0, t), respectively, for
the cases where the membrane inertia forces are taken into account and neglected. Apparently,
the influence of the membrane inertia forces is negligible. Moreover, Figs. 4a,b show the deflec-
tion w(0,0, t) and the inplane displacement u(−1.741,0, t) for various numbers of linear modes
employed in the Ritz method for reduction of the degrees of freedom as compared with the
results of the unreduced system. It is noteworthy that the use of more than 20 modes changes
the results negligibly. This reduction is significant because it allows us to obtain accurate results
using a small number of degrees of freedom.

Fig. 3. Time history of (a) transverse displacements w(0,0, t), (b) inplane displacement u(−1.741,0, t)
and (c) membrane force N

x
(−1.741,0, t)



Nonlinear dynamic analysis of viscoelastic membranes... 749

Fig. 4. (a) Transverse displacement w(0,0, t) and (b) inplane displacement u(−1.741,0, t) for different
numbers of the linear modes for reduction of degrees of freedom

Next, the free vibrations of the elastic (η=0) and the viscoelastic membrane (η=0.2) are
studied. The initial conditions are w(x,y,0) the deflection of the membrane due to a uniform
static load pz =1.0kN/m

2 and ẇ(x,y,0)= 0. The results were obtained using 20 linear modes
for the Ritz reduction. Figure 5 shows the response of themembrane for the elastic (η=0) and
viscoelastic (α=1, η=0.2) material.

Fig. 5. Time history of (a) w(0,0, t), (b) u(−1.741,0, t), (c) N
x
(−1.741,0, t) for elastic (η=0) and

viscoelastic (α=1, η=0.2) material

Finally, the forced vibrations of the viscoelastic membrane are studied. The membrane is
subjected to a transverse load pz = 1.0H(t)kN/m

2. The results were obtained using 20 linear
modes for reduction. Figure 6 presents the timehistory of themembrane for various values of the
order of the fractional derivative of the fractional Kelvin-Voigt type model. This figure shows
how the fractional order influences the viscoelastic response. Figure 7a shows the deflection
w(0,0, t) for the elastic and viscoelastic material as compared with the static deflection, while
Fig. 7b depicts the phase plane of the deflection at the center for α=0.5, η=0.5.

4.2. Forced vibrations under harmonic load. Resonance

The dynamic response of the viscoelastic membrane of example 4.1 under the transverse
harmonic excitation pz = p0 sin(Ωt) has been studied.The resultswere obtained using N =210



750 J.T. Katsikadelis

Fig. 6. Time history of (a) w(0,0, t) and (b) N
x
(−1.741,0, t) for various values of the order α (η=0.5)

Fig. 7. (a) Time history of deflection at the center of a viscoelastic membrane (α=0.5, η=0.5) and
(b) phase

boundary elements and M =122 internal collocation points resulting from 193 linear triangular
cells and 10 linear mode shapes for the reduction of the degrees of freedom.

First, we study the response of the conventional Kelvin-Voigt viscoelastic model (α = 1)
when the external frequency Ω is close to the lower natural frequency of the linear membra-
ne, which was found ω = 0.916. Figure 8a presents the amplitude of vibration (steady state
response) at the center of the membrane for the nonlinear problem together with that of the

Fig. 8. Amplitude-frequency curves at the center of the membrane for the linear and nonlinear
problem (a) and for various values of the viscous parameter for the classical Kelvi-Voigtmodel (α=1)

linear undampedmembrane (the nonlinear terms have been neglected). It is observed that the
amplitude of vibrations increase as the external frequency Ω approaches the natural frequency.
However, the amplitude curves of the nonlinear problem move to right, indicating hardening
nonlinear behavior, while two stable solutions are observed for some values of the excitation
frequency. Moreover, due to the viscoelastic material, the vibrations are bounded contrary to
the linear membrane problem. Figure 8b presents the amplitude-frequency curves for various
values of the viscous parameter η. The amplitude of vibrations increases as the viscous para-
meter decreases. For η < 1 a jump phenomenon appears and two sable solutions are observed.



Nonlinear dynamic analysis of viscoelastic membranes... 751

Figure 9a presents the time history of the deflection at the center of the membrane for two
values of the external frequency, a little before and after the jump. The frequency is changed at
t=2400s. Figure 9b presents the influence of the initial conditions on the steady state response
of the membrane.

Fig. 9. Time history of the central deflection for (a) two values of the external frequency (α=1,
η=0.5, p0 =0.01) and (b) two values of the initial velocity (α=1, η=0.1,Ω=0.927, p0 =0.12)

Similar jump phenomena appear for certain values of the excitation frequency when the
amplitude p0 of the external harmonic load varies within a certain interval. This is shown in
Fig.10, where two stable solutions appear for two values of the external frequency (Ω1 =0.922,
Ω2 =0.927) with varying the amplitude p0.

Fig. 10. Amplitude of the response of the membrane as a function of the external amplitude p0

Fig. 11. Amplitude-frequency curves for: (a) various values of the fractional derivative α (η=5),
(b) two values of the fractional derivative α (η=1)

Next, we study the response of the fractional Kelvin-Voigt model. Figures 11a and 11b
present the amplitude of the steady state response at the center of the membrane for various
values of the order of the fractional derivative α and for two values of the viscous parameter
(η = 5 and η = 1). It is observed that the amplitude of the steady state response increases as



752 J.T. Katsikadelis

the fractional derivative decreases. For α = 0.2, η = 1 and for appropriate initial conditions,
two stable solutions appear.

5. Conclusions

The governing equations describing the nonlinear dynamic response of viscoelastic membranes
made of fractional type derivative viscoelastic materials are derived. The resulting nonlinear
partial fractional differential equations are solved using the BEM developed recently by Katsi-
kadelis. The membrane may have an arbitrary shape. Free and forced vibrations are studied.
The response under a harmonic excitation is also highlighted. Various important conclusions are
drawn pertaining to the nonlinear dynamic response of viscoelastic membranes. Among them:

• Themembrane inertia forces have a negligible effect.
• The semidiscretized equations of motion constitute a system of coupledDuffing type ordi-
nary FDEs. Therefore, similar resonance phenomena due to the time harmonic excitation
are observed.

• The jump phenomena are bounded and they are considerably influenced by the order
of the fractional derivative in the Kelvin-Voigt. Moreover, a unique stable position can
result for an appropriate value of the order of the fractional derivative. Thus, the use of
membranes made of viscoelastic materials enables them to avoid disastrous phenomena
due to resonance.

In closing, thepresentedmethodprovides an efficient computational tool to analyzenonlinear
viscoelasticmembranes described by realisticmodels and enables the investigator to understand
their complicated dynamic response.Moreover, the subject of nonlinear viscoelastic membranes
involves interesting and important applications, computational issues and appliedmathematics.

References

1. AdolfssonK., EnelundM., OlssonP., 2005,On the fractional ordermodel of viscoelasticity,
Mechanics of Time-Dependent Materials, 9, 15-34

2. Amabili M., 2008,Nonlinear Vibrations and Stability of Shells and Plates, Cambridge University
Press

3. BabouskosN.,Katsikadelis J.T., 2009,Nonlinear vibrations of viscoelastic plates of fractional
derivative type. An AEM solution,Open Mechanics Journal, 3, 25-44

4. GalucioA.C., Deu J.F., OhayonR., 2004, Finite element formulation of viscoelastic sandwich
beams using fractional derivative operators,Computational Mechanics, 33, 282-291

5. GoncalesP.B., SoaresR.M., PamplonaD., 2009,Nonlinear vibrations of a radially stretched
circular hyperelastic membrane, J. of Sound and vibration, 327, 231-248

6. Jevine A.J., Mackintosch F.C., 2002, Dynamics of viscoelastic membranes, Physical Review,
E 66, 061606

7. Katsikadelis J.T., 2002,Boundary Elements. Theory and Applications, Elsevier, Oxford, U.K.

8. Katsikadelis J.T., 2008, The fractional wave-diffusion equation in bounded inhomogeneous ani-
sotropic media. An AEM solution, [In:] Advances in Boundary Element Methods, G.D. Manolis,
D. Polyzos (Eds.), 255-276, Springer Science, Dordrecht, Netherlands

9. Katsikadelis J.T., 2009, Numerical solution of multi-term fractional differential equations,
ZAMM Zeitschrift für Angewandte Mathematik und Mechanik, 89, 7, 593-608

10. Katsikadelis J.T., 2011, The BEM for Numerical solution of partial fractional differential equ-
ations,Computers and Mathematics with Applications, 62, 891-901



Nonlinear dynamic analysis of viscoelastic membranes... 753

11. Katsikadelis J.T., Babouskos N., 2010, Post-buckling analysis of viscoelastic plates with frac-
tional derivativemodels,Engineering Analysis with Boundary Elements, 34, 1038-1048

12. Koivurova H., Pramila A., 1997, Nonlinear vibration of axially moving membrane by finite
element method,Computational Mechanics, 20, 573-581

13. Meral F.C., Royston T.J., Magin R., 2010, Fractional calculus in viscoelasticity: An experi-
mental study,Commun. Nonlinear Sci. Numer. Simulat., 15, 939-945

14. Nayfeh A.H., Mook D.T., 1979,Nonlinear Oscillations, New York,Wiley

15. Nerantzaki M.S., Babouskos N., 2011, Analysis of inhomogeneous anisotropic viscoelastic
bodies described by multi-parameter fractional differential constitutive models, Computers and
Mathematics with Applications, 62, 891-901

16. Plaut R., 1990, Parametric excitation of an inextensible, air-inflated, cylindrical membrane, Int.
J. Non-linear Mechanics, 25, 253-262

17. Podlubny I., 1999,Fractional Differential Equations, Academic Press, NewYork

18. Schmidt A., Gaul L., 2002, Finite element formulation of viscoelastic constitutive equations
using fractional time derivatives,Nonlinear Dynamics, 29, 1/4, 37-55

19. Stiassnie M., 1979, On the application of fractional calculus for the formulation of viscoelastic
models,Applied Mathematical Modeling, 3, 300-302

20. Wineman A., 2007, Nonlinear viscoelastic membranes, Computers and Mathematics with Appli-
cations, 53, 168-181

Analiza nieliniowej dynamiki lepko-sprężystych powłok opisanych modelem

o pochodnej ułamkowej

Streszczenie

W pracy przedstawiono analizę dynamiki płaskiej powłoki lepko-sprężystej po wstępnym ugięciu.
Reologiczny model materiału powłoki opisano równaniem o pochodnych ułamkowych. Powłokę pod-
dano płaskiemu, poprzecznemu obciążeniu zewnętrznemu. Zagadnienie dynamiki ujęto układem trzech
sprzężonych nieliniowych równań różniczkowych typu hiperbolicznego o pochodnych ułamkowych. Rów-
nania rozwiązano metodą elementów brzegowych (BEM) sformułowaną przez autora właśnie dla ukła-
dów opisanych pochodnymi ułamkowymi. W prezentowanej pracy założono, że powłoka wykonana jest
z lepko-sprężystego materiału Kelvina-Voigta o pochodnej ułamkowej, choć samo sformułowanie BEM
nie wykluczamożliwości analizy innychmodeli reologicznych. Przedstawiono kilka przykładów symulacji
numerycznych pokazujących efektywność zastosowanej metody rozwiązywania oraz dających lepszy ob-
raz interesującej, lecz skomplikowanej dynamiki lepko-sprężystych powłok.Warto również podkreślić, że
w warunkach rezonansowych odnotowano zjawiska podobne do obserwowanychw układzie Duffinga.

Manuscript received November 7, 2011, accepted for print February 6, 2012