- Email: [email protected]
Post-buckling analysis of viscoelastic plates with fractional derivative models
Engineering Analysis with Boundary Elements 34 (2010) 1038–1048
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
Post-buckling analysis of viscoelastic plates with fractional derivative models J.T. Katsikadelis n, N.G. Babouskos School of Civil Engineering, National Technical University of Athens, Athens GR-15773, Greece
a r t i c l e in f o
a b s t r a c t
Article history: Received 17 February 2010 Accepted 26 June 2010 Available online 23 July 2010
The post-buckling response of thin plates made of linear viscoelastic materials is investigated. The employed viscoelastic material is described with fractional order time derivatives. The governing equations, which are derived by considering the equilibrium of the plate element, are three coupled nonlinear fractional partial evolution type differential equations in terms of three displacements. The nonlinearity is due to nonlinear kinematic relations based on the von Ka´rma´n assumption. The solution is achieved using the analog equation method (AEM), which transforms the original equations into three uncoupled linear equations, namely a linear plate (biharmonic) equation for the transverse deflection and two linear membrane (Poisson’s) equations for the inplane deformation under fictitious loads. The resulting initial value problem for the fictitious sources is a system of nonlinear fractional ordinary differential equations, which is solved using the numerical method developed recently by Katsikadelis for multi-term nonlinear fractional differential equations. The numerical examples not only demonstrate the efficiency and validate the accuracy of the solution procedure, but also give a better insight into this complicated but very interesting engineering plate problem & 2010 Elsevier Ltd. All rights reserved.
Keywords: Plates Large deflections Buckling Analog equation method Boundary element method Viscoelasticity Fractional derivatives
1. Introduction The stability and the post-buckling response of structures made of viscoelastic materials have attracted the interest in recent years. Viscoelastic materials, such as polymers, exhibit both viscous and elastic behavior and are extensively used in many modern engineering applications like aircrafts, ships, buildings and space structures due to the advantage of light weight and high strength. The buckling of structures subjected to external compressive loads is initiated by small imperfections which are unavoidable in practice. In the case of elastic material, when the loads exceed a critical value, the structure becomes unstable and reaches a new stable position instantaneously. However, for viscoelastic materials, buckling is not a classical stability problem, since it is an evolution phenomenon. As the loads exceed a critical value, the structure becomes unstable and needs some time to reach its final stable position, which depends on the viscoelastic characteristics of the material and the value of the initial imperfection. Moreover, for some type of viscoelastic materials with initial elastic response and when the external loads remain within certain values creep buckling occurs. In this case the structure is unstable but the displacements increase slowly with time. The time needed for high deformation rates is large enough extending the useful life of the structure [1,2].
n
Corresponding author. E-mail addresses: [email protected] (J.T. Katsikadelis), [email protected] (N.G. Babouskos). 0955-7997/$ - see front matter & 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2010.07.003
There have been proposed various linear and nonlinear models for the mechanical behavior of these materials, such as Maxwell, Kelvin–Voigt, Zener and multi-parameter models. The constitutive equations can be presented in differential or in hereditary integral form. Recently, many researchers have shown that differential viscoelastic models with fractional derivatives are in better agreement with the experimental results than the integer derivative models [3–6]. Adolfsson et al. [6] proved that complicated multi-parameter integer viscoelastic models converge to fractional models with one single parameter. While the post-buckling of viscoelastic beams using integral or integer derivative models has been investigated by many authors [2,7–12] using analytic or approximate solution techniques, the literature on buckling of viscoelastic plates is rather limited. The reason is that the governing equations are linear or nonlinear evolution partial differential equations whose solution is difficult to obtain, especially for plates of arbitrary shape. Hewitt and Mazumdart [13] studied the stability of linear viscoelastic isotropic plates of arbitrary geometry using the method of constant deflection contours. Shalev and Aboudi [14] investigated the post-buckling response of laminated viscoelastic plates using Laplace transform and taking into account shear deformation. Touati and Cederbaum [15] studied the post-buckling of nonlinear viscoelastic laminated plates using the integral representation formulation. Chen and Cheng [16] studied the instability of nonlinear viscoelastic plates using constitutive equations in integral formulation. In most cases the plate problem is solved using Galerkin method in conjunction with a numerical integration scheme.
J.T. Katsikadelis, N.G. Babouskos / Engineering Analysis with Boundary Elements 34 (2010) 1038–1048
With regard to the fractional derivative models the viscoelastic problems treated are quite few. The reason is obvious. The equations describing the response of the systems are complicated linear or nonlinear fractional partial differential equations whose solution is very difficult to obtain. There have been recently some rather ad hoc solutions for the analysis of viscoelastic beams and ¨ plates. Grunwald formalism for the fractional order derivatives have been used by Galucio et al. [17], who analyzed viscoelastic sandwich beams with FEM and Schmidt and Gaul [18], who analyzed three dimensional viscoelastic bodies using FEM. Stankovic and Atanackovic [19] studied the vibrations of an axially compressed viscoelastic rod using the fractional Kelvin– Voigt model. Rossikhin and Shitikova [20,21] treated linear and nonlinear vibrations of rectangular viscoelastic plates with damping modeled with Riemann–Liouville fractional derivative using an approximate analytical solution. Recently Katsikadelis developed the AEM for solving linear and nonlinear fractional differential equations [22]. This method is general and has been already used to solve the linear fractional diffusion–wave equation in bounded inhomogeneous anisotropic bodies [23], the nonlinear dynamic response of viscoelastic plates and membranes [24,25]. In this paper this method is employed to study the post-buckling response of thin viscoelastic plates. Without excluding other models, the employed herein model is the generalized Kelvin–Voigt model with fractional order derivatives. The plate may have an arbitrary shape and be subjected to all types of conventional boundary conditions including also elastic support. For the two dimensional elasticity, the stress–strain relations can be reduced to the form [24] 38 9 ex > > = < 6 7 n 1 0 7 ey sy ¼ E 2 6 4 5> 1 n > > 1 n ; :t ; :g > 0 0 xy xy 2 2 38 a 9 1 n 0 D e > < c x > = ZE 6 0 7 6n 1 7 Dac ey þ 4 5 2 1n > 1n : Da g > ; 0 0 c xy 2 9 8 s = > < x>
2
1
n
1039
deformation under fictitious loads that are unknown in the first instance. Application of the D/BEM (Domain BEM) yields the initial value problem for the fictitious sources which is a system of nonlinear fractional ordinary differential equations. This system is solved utilizing the numerical method developed recently by Katsikadelis for multi-term linear and nonlinear fractional differential equations [22]. The post-buckling response of several viscoelastic plates is studied. The obtained numerical results demonstrate the efficiency of the developed solution method and validate its accuracy. Moreover the conclusions drawn from this investigation permit a better insight into this complicated but very interesting engineering plate problem.
2. Governing equations 2.1. The nonlinear plate problem We consider a thin elastic plate of uniform thickness h occupying the two dimensional multiply connected domain O of the xy plane with boundary G [Ki¼ 0 Gi (Fig. 1). The nonintersecting curves Gi (i¼0, 1, 2, y, K) may be piece-wise smooth. The boundary may be simply supported, clamped, free or elastically supported with transverse stiffness kT ðxÞ and rotational stiffness kR ðxÞ x : ðx,yÞ A G. The plate is subjected to inplane loads px, py (body forces). Moreover, along the movable edges inplane forces Nn , Nt may act. The von Ka´rma´n assumption for the kinematic relation is adopted. That is 1 2
1 2
ex ¼ u, x þ w, 2x ey ¼ v, y þ w, 2y gxy ¼ u, y þ v, x þw, x w, y
0
ð1Þ
where E, n are the elastic material constants, Z is the viscous parameter and Dac the Caputo fractional derivative of order a defined as 8" # Z t > 1 uðmÞ ðtÞ > > > , m1 o a o m d t < GðmaÞ a þ 1m 0 ðttÞ ð2Þ Dac uðtÞ ¼ > dm > > > m¼a : m uðtÞ dt where m is a positive integer. The advantage of this definition for the fractional derivative is that it permits the assignment of initial conditions which have direct physical significance [26]. Apparently, the classical derivatives result for integer values of a. For a ¼1 the conventional Kelvin–Voigt model is obtained while for a ¼0 the material is purely elastic with Young modulus E* ¼E(1+ Z). The governing equations result by taking the equilibrium of the plate element in a slightly deformed configuration and using the constitutive equations (1). They are three coupled nonlinear fractional partial differential equations in terms of displacements. The nonlinearity is due to nonlinear kinematic relations based on the von Ka´rma´n assumption. The concept of the analog equation is applied [27,28] to convert the original nonlinear fractional partial equations into three uncoupled linear equations, namely a linear plate (biharmonic) equation for the transverse deflection and two linear membrane (Poisson’s) equations for the inplane
ð3a; b; cÞ
where u ¼ uðx,y,tÞ, v ¼ vðx,y,tÞ are the membrane and w ¼ wðx,y,tÞ the transverse displacements. On the basis of the Eqs. (1) and (3) the stress resultants are written as 1 2 w, x þ nw, 2y Nx ¼ C u, x þ nv, y þ 2 1 ð4aÞ þ ZC Dac u, x þ nDac v, y þ Dac w, 2x þ nw, 2y 2 1 Ny ¼ C nu, x þ v, y þ nw, 2x þw, 2y 2 1 a þ ZC nDc u, x þDac v, y þ Dac nw, 2x þ w, 2y 2
x
Γ0
ð4bÞ
ss *
Nn px (x, y)
y
c
Γk
f hole py (x, y)
*
Nt
(Ω) corner ss Fig. 1. Plate geometry and supports (c¼clamped, ss ¼ simply supported, f¼ free).
1040
J.T. Katsikadelis, N.G. Babouskos / Engineering Analysis with Boundary Elements 34 (2010) 1038–1048
1n u, y þ v, x þ w, x w, y 2 1n a 1 Dc u, y þ Dac v, x þ Dac w, x w, y þ ZC 2 2
Nxy ¼ C
ð4cÞ
Mx ¼ D w, xx þ nw, yy ZD Dac w, xx þ nDac w, yy
My ¼ D nw, xx þw, yy ZD nDac w, xx þDac w, yy
ð5aÞ
The star designates a prescribed quantity. In Eq. (8), Vw is the equivalent elastic shear force, Mw is the elastic normal bending moment; Tw is the elastic twisting moment along the boundary and ½Twk is its discontinuity jump at corner k. The operators producing these quantities are given as [29] " !# @ 2 @ @2 @ k V ¼ D r þ ð1nÞ ð10aÞ @n @s @s@n @s
ð5bÞ "
Mxy ¼ Dð1nÞw, xy ZDð1nÞDac w, xy
ð5cÞ
where C ¼Eh/(1 n2) is the membrane stiffness and D ¼ Eh3 =12ð1n2 Þ the flexural stiffness of the plate. The governing equations result by taking the equilibrium of the plate element in a slightly deformed configuration. This yields @2 Mxy @2 Myy @2 Mxx @ Nx w, x þ Nxy w, y þ þ2 þ 2 @x @x@y @x @x2 @ Nxy w, x þNy w, y ¼ 0 þ @y
ð6aÞ
2
M ¼ D r ð1nÞ
T ¼ Dð1nÞ
!# @2 @ þ k @n @s2
ð10bÞ
! @2 @ k @s @s@n
ð10cÞ
where k ¼ k(s) is the curvature of the boundary and n, s the intrinsic boundary coordinates, i.e. the coordinates along the normal and tangential directions at a boundary point. 2.2. The linear plate problem
@Nx @Nxy þ þpx ¼ 0 @x @y
ð6bÞ
@Nxy @Ny þ þpy ¼ 0 @x @y
ð6cÞ
Using Eqs (4)–(6) we obtain the plate equations in terms of displacements in O as i. For the transverse deflection 4 4 Dr wþ ZDDac r wNx w, xx 2Nxy w, xy Ny w, yy þ px w, x þpy w, y ¼ 0 ð7aÞ ii. For the inplane deformation
1þn 2 r2 u þ u, x þ v, y , x þw, x w, xx þw, yy 1n 1n 1þ n 1þn 2 w, xy w, y þ ZDac r u þ u, x þv, y , x þ 1n 1n
2 1þn px ¼0 w, xx þ w, yy þ w, xy w, y þ þw, x 1n 1n Gh
4
Dr wlðNx w, xx þ2Nxy w, xy þ Ny w, yy Þ ¼ 0 Vw þ lðNn* w, n þ Nt* w, t Þ þkT w ¼ 0 Mw þkR w, n ¼ 0 ð7bÞ
1þn 2 ðu, x þv, y Þ, y þ w, y w, yy þ w, xx 1n 1n 1þ n 1þn w, xy w, x þ ZDac r2 vþ ðu, x þ v, y Þ, y þ 1n 1n
py 2 1þn ¼0 w, yy þ w, xx þ w, xy w, x þ þw, y 1n 1n Gh
r2 v þ
The plate problem is linearized if the stretching of the middle surface due to bending is neglected, that is if it is set w, 2x ¼ w, 2y ¼ w, x w, y ¼ 0 in Eq. (4). This implies that, although the inplane forces contribute to bending, they are not influenced by it. Thus the nonlinear terms in Eq. (7) vanish. The boundary conditions remain the same for both problems. In this case the two problems are uncoupled, hence the inplane forces Nx, Ny, Nxy are obtained from the solution of the Eq. (7b), (7c), which are solved independently. For linear elastic buckling Z ¼0 and the eigenvalue problem is derived as
ð7cÞ
ðkÞ kðkÞ T w ½Twk ¼ 0
or
w, n ¼ 0 or
or
in O
w¼0
on G
on G
wðkÞ ¼ 0
at corner point k
ð11Þ ð12aÞ ð12bÞ ð12cÞ
where l is the load parameter. The linear eigenvalue problem is also solved in this investigation, because the mode shapes (eigenvectors) of this problem are used as the initial imperfection for the post-buckling analysis. Besides, the linear elastic buckling loads are needed to compare the elastic with the viscoelastic plate response. The AEM solution for the linear and the nonlinear plate problems is presented below.
The associated boundary conditions result as Vw þ ZDac ðVwÞ þ Nn* w, n þ Nt* w, t þ kT w ¼ 0 Mw þ ZDa ðMwÞ þ k
R w, n
c
¼0
ZDa ½Tw
ðkÞ kðkÞ T w ½Twk
c
k
or ¼0
w,n ¼ 0 or
ðkÞ
w
or
w¼0
on G ¼0
on G
ð8aÞ 3. The AEM solution ð8bÞ
at corner point k ð8cÞ
Nn ¼ Nn*
or
un ¼ u*n
Nt ¼ Nt*
or
ut ¼ u*t
on G
ð8dÞ
The initial boundary value problem (7a), (8a)–(8c), (9a) for the response of the plate is solved using the AEM. The analog equation for the problem at hand is
r4 w ¼ bðx,tÞ, x ¼ fx,yg A O on G
ð8eÞ
Moreover, the viscoelastic plate is subjected to the initial conditions wðx,0Þ ¼ gðxÞ,
3.1. The plate problem
uðx,0Þ ¼ g1 ðxÞ,
vðx,0Þ ¼ g2 ðxÞ, in O
ð9a; b; cÞ
ð13Þ
where bðx,tÞ represents the time dependent fictitious load. Eq. (13) is a quasi-static equation, that is the time appears as a parameter, and it can be solved with the boundary conditions (8a)–(8c) at any instant t using the BEM. Thus, the solution at a
J.T. Katsikadelis, N.G. Babouskos / Engineering Analysis with Boundary Elements 34 (2010) 1038–1048
point xA O is obtained in integral form as Z Z wðx,tÞ ¼ w* bdO þ ðw* Vw þ w, n Mw* w* , n MwwVw* Þds O
X
1
ð14Þ
k
which for x A G yields the following two boundary integral equations for points where the boundary is smooth Z Z 1 wðx,tÞ ¼ w* bdO þ ðw* Vw þ w, n Mw* w* , n MwwVw* Þds 2 O G X w* ½Tww½Tw* ð15Þ
y
w*1 bdO þ
Z
M domain points
6 5
4
Fig. 3. Boundary and domain discretization.
ðw*1 Vw þ w, n Mw*1 w*1 , n MwwVw*1 Þds
O G X w*1 ½Tww½Tw*1 k
ð16Þ
k
in which w* ¼ w* ðx,yÞ, x,y A G, is the fundamental solution and w*1 its normal derivative at point x A G, i.e. 1 2 r ln r 8p
1 2 1 r ln r , n ¼ rr, n ð2 ln r þ 1Þ w*1 ¼ 8p 8p
w* ¼
ð17a; bÞ
n is the unit normal vector to the boundary at point x, whereas n is the unit normal vector to the boundary at the integration point y and r ¼ :xy: (see Fig. 2). Eqs. (15) and (16) can be used to establish the unspecified boundary quantities. They are solved numerically using the BEM. The boundary integrals are approximated using N constant boundary elements. Although the domain integrals can be approximated by radial basis functions and using the DRM to convert them to boundary integrals, maintaining thus the pure boundary character of the method [28], we approximate them with M linear triangular elements. This is more practical because the domain discretization is performed automatically using the Delaunay triangulation utilizing ready-touse subroutines. In this case, attention should be paid to place the nodal points of the triangles adjacent to the boundary on their sides, because the fictitious source is not defined on the boundary (see Fig. 3). Thus, after discretization and application of the boundary integral Eqs. (15) and (16) at the N boundary nodal points and Eq. (15) at the Nc corner points we obtain 8 9 8 9 > > = = H wc ¼ G R þAb ð18Þ > > : w, > ; : > ; M n
where H, G are (2N+ Nc) (2N + Nc) known coefficient matrices originating from the integration of the kernel functions on the boundary elements; A is an (2N + Nc) M coefficient matrix originating from the integration of the kernel function on the domain triangular cells; w,wc ,w, n are the vectors of the N boundary nodal displacements, Nc corner displacements and N boundary nodal normal slopes, respectively, and V, R, M are the vectors of the N nodal values of effective shear force, Nc concentrated corner forces and N nodal values of the normal bending moment. Finally, b is the vector of the M nodal values of the fictitious source inside O. Eq. (18) constitutes a system of 2N + Nc equations for 4N +2Nc +M unknowns. Additional 2N + Nc equations are obtained from the boundary conditions. Thus, the boundary conditions (8a)–(8c), when applied at the N boundary nodal points and the Nc corner points yield the set of equations
a1 w þ a2 w, n þ a3 V ¼ 0 b1 w, n þ b2 M ¼ 0 c1 wc þc2 R ¼ 0 ð19a; b; cÞ where a1, a2, a3, b1, b2, c1, c2 are known coefficient matrices. Note that Eq. (19a) has resulted after approximating the derivative w,t in Eq. (8a) with a finite difference scheme. Eqs. (18) and (19) can be combined and solved for the boundary quantities w, wc , w, n , V, R, M in terms of the fictitious load b. Subsequently, these expressions are used to eliminate the boundary quantities from the discretised counterpart of Eq. (14). Thus we obtain the following representation for the deflection wðx,tÞ ¼
M X
bk ðtÞWk ðxÞ
xAO
ð20Þ
k¼1
The derivatives of wðxÞ at points x inside O are obtained by direct differentiation of Eq. (14). Thus, we obtain after elimination of the boundary quantities w, pqr ðx,tÞ ¼
M X
bk ðtÞWk , pqr ðxÞ,
p,q,r ¼ 0,x,y
xAO
ð21Þ
k¼1
t
Γ
N boundary points
x
k
k
Z
Nc corner points
G
w* ½Tww½Tw*
k
1 w, n ðx,tÞ ¼ 2
3
2
1041
n y
r = ||x - y||
r = ||x - y||
where Wk , pqr are known functions. Note that the above notation implies w,000 ¼w, w,0y0 ¼w,y, etc. 3.2. The plane stress problem
x x
y
a x
corner k Fig. 2. BEM notation.
Noting that Eqs. (7b) and (7c) are of the second order their analog equations are obtained using the Laplace operator. This yields
r2 u ¼ b1 ðx,tÞ r2 v ¼ b2 ðx,tÞ
ð22a; bÞ
The integral representation of the solution of Eq. (22a) is Z Z euðxÞ ¼ v* b1 dO ðv* qq* uÞds x A O [ G ð23Þ O
G
1042
J.T. Katsikadelis, N.G. Babouskos / Engineering Analysis with Boundary Elements 34 (2010) 1038–1048
in which q ¼u,n; v* ¼ ‘nr=2p is the fundamental solution to Eq. (22a) and q* ¼ v, *n its derivative normal to the boundary; r ¼ :yx: x A O [ G and y A G; e is the free term coefficient (e ¼1 if x A O, e ¼1/2 if x A G and e ¼0 if x= 2O [ G). Using the BEM with constant boundary elements and linear triangular domain elements and following the same procedure applied for the plate equation, we obtain the following representation for the inplane displacement u and its derivatives: u, pq ðx,tÞ ¼
M X
bð1Þ ðtÞUkð1Þ , pq ðxÞ þ k
k¼1
M X
bð2Þ ðtÞUkð2Þ , pq ðxÞ k
k¼1
þ U0 , pq ðxÞ, xA O, p,q ¼ 0,x,y
ð24Þ
Eq. (28) yields the buckling load parameter l and the buckling mode shapes. Buckling occurs due to an imperfection of the plate, which is given as a small initial deflection. In our examples the form of the first buckling mode of the linear elastic plate is taken as initial deflection. Subsequently, we increase the inplane loads and we study the response of the plate. When the inplane loads become greater than a critical value instability occurs, the displacements increase with time and the plate reaches a new stable position. The initial value problem Eqs. (26) and (27) is solved using a numerical method which has been developed recently by Katsikadelis for multi-term linear and nonlinear fractional differential equations [22].
Similarly, we obtain for the displacement v v, pq ðx,tÞ ¼
M X
bð1Þ ðtÞVkð1Þ , pq ðxÞ þ k
k¼1
M X
5. Examples
bð2Þ ðtÞVkð2Þ , pq ðxÞ k
k¼1
þ V0 , pq ðxÞ, x A O, p,q ¼ 0,x,y
ð25Þ
Ukð1Þ , Ukð2Þ , Vkð1Þ , Vkð2Þ , U0 , V0 are known functions; U0 , V0 result from the nonhomogeneous boundary conditions. Note that u,00 ¼u and v,00 ¼v.
4. The final step of the AEM Eqs. (21), (24) and (25) give the displacements w(x,t), u(x,t), v(x,t) and their derivatives provided that the three fictitious ð1Þ ð2Þ sources bðtÞ, b ðtÞ, b ðtÞ are first established. This is achieved by working as following. Collocating the fractional PDEs (7) at the M internal nodal points and substituting the expressions for the transverse deflection, Eq. (21), and the membrane displacements, Eqs. (24) and (25), we obtain the following system of 3M nonlinear fractional equations for bk ðtÞ, bð1Þ ðtÞ, bð2Þ ðtÞ, (k¼1, y, M) k k Fðb,b ,b ,Dac b,Dac b ,Dac b Þ ¼ P ð1Þ
ð2Þ
ð1Þ
ð2Þ
ð26Þ
where P is the vector of the external loads. Collocating the initial conditions Eq. (9) at the M internal nodal points and using Eqs. (21), (24) and (25), the initial conditions for the fictitious loads bk ðtÞ, bð1Þ ðtÞ, bð2Þ ðtÞ are obtained as k k n oT ð1Þ bð0Þ ¼ W1 g, b ð0Þ ¼ S1 g1 g2 þH1 , n oT ð2Þ þ H2 ð27Þ b ð0Þ ¼ S2 g1 g2 where W is N N known matrix, S1, S2 are N 2N known matrices and H1, H2 known vectors. The solution of the linear eigenvalue problem, Eqs. (11) and (12), is obtained following a similar procedure. Thus, collocating Eq. (11) at the M internal nodal points and substituting the expressions for the transverse deflection and its derivatives, Eq. (21), we obtain the following eigenvalue problem: ðKlGÞb ¼ 0
ð28Þ
Example 1. The post-buckling response of the square plate of thickness h¼0.1 m is investigated (Fig. 4). All edges are simply supported and movable in the plane of the plate. The plate is subjected to an inplane line load along two opposite sides. The results were obtained with N ¼196 boundary points and M ¼97 internal points resulting from 145 triangular cells (Fig. 5). The parameters of the viscoelastic material are E¼ 21 106 kN/m2, n ¼0.25, Z ¼0.5. The computed linear buckling load of the elastic plate is Pel ¼4681 kN/m (exact Pel ¼4606 kN/m). The initial deflection has the form of the first buckling mode of the elastic plate. First we analyze a viscoelastic plate with the conventional Kelvin–Voigt model (a ¼ 1). Fig. 6 presents the time history of the
w = Mw = 0 Nn = Nt = 0 4
w = Mw = 0
w = Mw = 0
Nt = 0
4 Nt = 0 Nn = -P
Nn = -P w = Mw = 0 Nn = Nt = 0
Fig. 4. Geometry and boundary conditions of plate in Example 1.
4
3
2
where K, G are known N N matrices given as Kik ¼ Ddik
i Gik ¼ Nxi Wki , xx þ 2Nxy Wki , xy þ Nyi Wki , yy
ð29Þ
with i, k¼1, y, M, dik the Kronecker delta and Wki , xx ¼ Wk , xx ðxi Þ, Wki , xy ¼ Wk , xy ðxi Þ, Wki , xx ¼ Wk , yy ðxi Þ; b is a vector with the nodal values of the fictitious load, which now does not depend on time. The inplane forces Nx, Ny, Nxy are established by solving independently the linear plane stress problem (see Section 2.2) using the direct BEM [30]. The solution of the eigenvalue problem
1
0
0
1
2
3
4
Fig. 5. Boundary and domain discretization of plate in Example 1.
J.T. Katsikadelis, N.G. Babouskos / Engineering Analysis with Boundary Elements 34 (2010) 1038–1048
0.014
-4200 P = 4800
P = 4700 P = 4650
0.012 0.01
P = 5000
P = 5200
-4600
0.008
Nx
w/h
1043
0.006
-5000
0.004 0.002 0 0
20
40 60 time (sec)
80
-5400
100
0
50
100
150
200
time (sec)
Fig. 6. Time history of the deflection at the center of the plate for a ¼ 1 just before and after the elastic buckling load.
Fig. 9. Time history of the membrane force Nx at point (1.4, 2) for a ¼ 1.
5400
1 P = 4800
P = 5000
P = 5200
5300
0.8
5200 5100
w/h
load P
0.6
0.4
5000 t = 20 sec t = 40 sec t = 60 sec t = 100 sec t = 160 sec Elastic with E
4900 4800
0.2
4700 4600
0 0
50
100 time (sec)
150
200
0
0.2
0.4
0.6
0.8
1
w/h
Fig. 7. Time history of the deflection at the center of the plate for a ¼1.
Fig. 10. Load–deflection curves at the center of the plate at different values of time for a ¼ 1.
150
observed that viscoelastic plates become unstable when the inplane load reaches the linear elastic buckling load and displacements increase with time until they reach a final value, in contrast to the elastic plates, which reach the final stable position immediately. The viscoelastic plate has a softer performance in buckling than the elastic plate. Fig. 10 presents the load–deflection curves for different values of time and the results are compared with the post-buckling path of the elastic plate [31]. It is observed that after certain time the viscoelastic post-buckling path approaches the elastic post-buckling path. Next, we study the response of the generalized Kelvin–Voigt model. Fig. 11 presents the time history of the deflection at the center of the plate for various values of a and for two cases of the inplane line load, just before and after the linear buckling load of the elastic plate. It is observed that in all cases the displacements increase as the inplane load becomes greater than the elastic buckling load, which means that the plate becomes unstable. However, for a o1 the displacements increase slowly and only when the load becomes greater than another critical value P 0 (P 0 oPel) the displacements increase rapidly. Such a gradual buckling, which occurs for loads Pel oPoP 0 , is known in literature as creep buckling and it is also observed in integer derivative viscoelastic models with initial elastic response [1,2]. Fig. 12 presents the time history of the
P = 4800
P = 5000
P = 5200
Mx
100
50
0 0
50
100
150
200
time (sec) Fig. 8. Time history of the bending moment at the center of the plate for a ¼ 1.
deflection at the center of the plate for two values of the inplane load P just before and after the buckling load Pel of the elastic plate. Figs. 7–9 present the time history of the deflection and the stress resultants for various values of the inplane load. It is
1044
J.T. Katsikadelis, N.G. Babouskos / Engineering Analysis with Boundary Elements 34 (2010) 1038–1048
6.2
7500
x 10-3
7400
α = 0.8 α = 0.5
6
7300 load P
α = 0.2 w/h
5.8
7200 7100
5.6
7000
α = 0.2
5.4
α = 0.5
6900
α = 0.8 0
t = 0.2 sec t = 0.5 sec t = 2 sec t = 50 sec elastic with E* = (1+η)E
2
4
0
6 time (sec)
8
0.1
0.2
0.3
0.4
10
0.7
7500
α=1 α = 0.2 α = 0.01 Elastic with E Elastic with E* =(1+η)E
7000 P = 7000
6500 load P
P = 7100 P = 7300 P = 7400
0.6
0.6
Fig. 14. Load–deflection curves at the center of the plate at different values of time for a ¼ 0.01.
Fig. 11. Time history of the deflection at the center of the plate for various values of a just before and after the elastic buckling load (solid line: P ¼ 4700, dashed line: P¼ 4650).
0.8
0.5
w/h
6000
w/h
5500 0.4
5000 4500
0.2
0
0 0
0.5
1 time (sec)
1.5
2
0.2
0.4 w/h
0.6
0.8
Fig. 15. Load–deflection curves for various values of the order a of the fractional derivative and comparison with the elastic post-buckling paths.
y
Fig. 12. Time history of the deflection at the center of the plate for a ¼0.01.
w = Mw = 0 u n = ut = 0 Nn = -P Nt = 0 w = Mw = 0
7000
x
load P
6500
6000 t = 5 sec t = 10 sec t = 20 sec t = 30 sec t = 100 sec t > 1000 sec
5500
5000 0
0.2
0.4
0.6
0.8
1
1.2
P
1.4
w/h Fig. 13. Load–deflection curves at the center of the plate at different values of time for a ¼ 0.2.
2 5 Fig. 16. Geometry and boundary conditions of plate in Example 2.
J.T. Katsikadelis, N.G. Babouskos / Engineering Analysis with Boundary Elements 34 (2010) 1038–1048
1045
4500
deflection at the center of the plate of various values of the inplane load for a ¼0.01. Figs. 13 and 14 present the load–deflection curves for different values of time for fractional order derivative a ¼0.2
4000
load P
5 3500 t = 2 sec
3000
2.5
t = 10 sec t = 100 sec
2500
t = 200 sec
0
t > 500 sec 2000 0
-2.5
0.1
0.2
0.3
0.4
0.5
0.6
w/h Fig. 20. Load–deflection curves at point r¼ 3.5 for various values of time (a ¼0.2).
-5 -5
-2.5
0
2.5
5
5000
Fig. 17. Boundary and domain discretization annular plate in Example 2.
4800
x
load P
1.3
10-3
1.2
t = 0.4 sec 4400
α = 0.5 α=1
t = 1 sec t = 2 sec
α = 0.2
t = 10 sec
4200
w/h
1.1
4600
t > 50 sec
1
4000 0
α = 0.5 5
10 time (sec)
15
20
Fig. 18. Time history of deflection at point r ¼3.5 for various values of a and for two cases of the inplane load, just before and after the elastic buckling load (solid line: P¼ 2250, dashed line: P ¼2200).
load P
3500
3000 t = 1sec t = 10sec
2500
t = 20 sec t > 100 sec 2000 0
0.1
0.2
0.3
0.4
Fig. 21. Load–deflection curves at point r ¼3.5 for various values of time (a ¼0.01).
0.8 0
0.2 w/h
α = 0.2
α=1
0.9
0.1
0.3
0.4
0.5
0.6
w/h Fig. 19. Load–deflection curves at point r ¼3.5 for various values of time (a ¼ 1).
and 0.01, respectively. It is observed that the critical value P 0 , above which high deformation rates appear, increases as the order of the fractional derivative decreases, and for a-0 the viscoelastic postbuckling response approaches the elastic post-buckling response with a modulus E*¼ (1+ Z)E. Finally, Fig. 15 presents the load– deflection curves after considerable time has elapsed. It is observed that for a ¼0.2 the critical load is P 0 ¼5400 and for a ¼0.01 is * ¼ 7021 of the P 0 ¼6950 very close to the elastic buckling load Pel plate with modulus E*. Example 2. The post-buckling response of the annular plate of Fig. 16 is investigated. The boundary conditions and the loading are shown in Fig. 16. The plate is simply supported at both boundaries. The inplane displacements are restrained along the external boundary. The results were obtained with N¼300 boundary points and M¼130 internal points resulting from 190 triangular cells (Fig. 17). The parameters of the viscoelastic material are E¼ 11 106 kN/m2, n ¼0.25, Z ¼1. The thickness of the plate is h¼0.1 m. The buckling load of the elastic plate is Pel ¼2224 kN/m (a FEM solution with 2000 elements yields Pel ¼ 2195 kN/m). Fig. 18 presents the time history of the deflection at point r¼3.5 for various values of a and for two cases of the inplane line load, just before and after the buckling load of the elastic plate. It is observed that in all cases the plate becomes unstable as the inplane load becomes greater than the elastic buckling load. For a o1 creep buckling occurs. Figs. 19–21
1046
J.T. Katsikadelis, N.G. Babouskos / Engineering Analysis with Boundary Elements 34 (2010) 1038–1048
present the load–deflection curves at point r¼3.5 for different values of time and for a ¼1, 0.2, 0.01. Finally, Fig. 22 presents the load–deflection curves for various values of the order of fractional derivative after considerable time has elapsed. The results are compared with the post-buckling path of the elastic plate with elastic constants E and E*¼E(1+ Z).
2
1.5 Example 3. The post-buckling response of the rectangular plate of Fig. 23 is investigated. The plate is clamped along side 1–4, simply supported movable along sides 1–2 and 3–4 and free to move along side 2–3. The plate is compressed along sides 1–2 and 3–4 with an inplane line load P. The results were obtained with N ¼200 boundary points and M ¼93 internal points resulting from 136 triangular cells (Fig. 24). The thickness of the plate is h¼0.01 m. The parameters of the viscoelastic material are E¼29 106 kN/m2, n ¼0.3, Z ¼1. The buckling load of the elastic plate is Pel ¼36.2 kN/m (a FEM solution with 200 elements yields
1
0.5
5000 4500 α=1 α = 0.5 α = 0.2 α = 0.01 Elastic with E Elastic with E* = E(1+η)
load P
4000 3500
0 0
0.5
1
Fig. 24. Boundary and domain discretization of plate in Example 3.
3000
50 2500
α=1 2000 0
0.1
0.2
0.3
0.4
α = 0.3
0.5
45
Fig. 22. Load–deflection curves at point r¼ 3.5 for various values of the order a of the fractional derivative and comparison with the elastic post-buckling paths.
Elastic
load P
w/h
40
y Nn = -P Nt = 0 w = Mw = 0
35 0
1
4
1
1.5
2
w/h Fig. 25. Load–deflection curves for two cases of the order a of the fractional derivative and comparison with the elastic post-buckling path.
3
un = ut = 0 2 w = w,n = 0
Nn = Nt = 0 Vw = Mw = 0
1
0.5
2
x
35.44 kN/m). Fig. 25 presents the load–deflection curves at point (1, 1) for two cases of the order a ¼1, 0.3, after considerable time has elapsed and the results are compared to the post-buckling path of the elastic plate. For a ¼0.3 the critical load is P 0 ¼ 39 kN/m, above which the displacements of the plate increase rapidly. Fig. 26a and b presents the time history of the deflection at point (1, 1) and the inplane displacement at point (0.16, 1) for various values of the initial imperfection (P ¼40 kN/m).
6. Conclusions
Nn = -P Nt = 0 w = Mw = 0 Fig. 23. Geometry and boundary conditions of plate in Example 3.
The post-buckling response of thin plates made of linear viscoelastic material has been investigated in this work. The behavior of the viscoelastic material is described by the
J.T. Katsikadelis, N.G. Babouskos / Engineering Analysis with Boundary Elements 34 (2010) 1038–1048
0
2
x 10-5
w0/h = 0.07
-0.2
1.5
w0/h = 0.007 -0.4
1
w0/h = 0.0007
u
w/h
1047
w0/h = 0.07
-0.6
w0/h = 0.007
0.5
-0.8
w0/h = 0.0007
-1
0 0
20
40
60 80 time (sec)
100
120
140
0
20
40
60 80 time (sec)
100
120
140
Fig. 26. Time history of the deflection w at point (1, 1) and the inplane displacement u at point (0.16, 1) for various values of the initial imperfection (P ¼ 40 kN/m).
generalized Kelvin–Voigt model with fractional order derivative. The solution of the derived governing equations is achieved using the AEM. This method converts the coupled nonlinear fractional PDEs, describing the response of the viscoelastic plate, into uncoupled linear fractional PDEs, which are subsequently solved by the D/BEM. The discretization procedure yields an initial value problem for simultaneous ordinary FDEs, which is solved by a new method for nonlinear multi-term fractional differential equations developed recently by Katsikadelis. The influence of the viscoelastic material on the post-buckling response of the plate has been studied. The conclusions that can be drawn from this study are summarized as:
(a) The viscoelastic plate becomes unstable when the inplane load reaches the elastic buckling load, which is obtained by the linear eigenvalue problem. Then the post-critical response is studied by gradually increasing the inplane load. (b) While the post-buckling response for an elastic material is initiated when the inplane load reaches the critical value, for a viscoelastic material the post-buckling response starts after a certain time, which depends on the viscoelastic parameters and the value of the initial imperfection. (c) For the conventional Kelvin–Voigt model (a ¼1) the viscoelastic post-buckling path converges to the elastic postbuckling path (Z ¼0). Practically, after certain time, which depends on Z and the value of the initial imperfection, the elastic and viscoelastic post-buckling paths become identical. (d) For the generalized Kelvin–Voigt model (0o a o1) creep buckling occurs. For certain values of the external load Pel oPoP 0 the plate becomes unstable but the displacements increase slowly. When the external load exceeds the critical value P 0 high deformation rates appear and the plate reaches a new stable position. For 0 o a o1, after certain time, the viscoelastic post-buckling path lies between the elastic paths with moduli E and E* ¼E(1+ Z). (e) As a-0 the viscoelastic post-buckling response of the plate approaches an elastic post-buckling response with a modulus E* ¼E(1+ Z). (f) The fractional type viscoelastic model with 0 o a o1 can be visualized as a model of combined springs and dashpots,
whose contribution depends on the value of a. This is reflected in the creep buckling response of the analyzed structures, as it happens with more complicated integer derivative models with initial elastic response.
References [1] Flugge W. Viscoelasticity. Blaisdell Publishing Company; 1967. [2] Vinogradov AM. Nonlinear effects in creep buckling analysis of columns. Journal of Engineering Mechanics 1985;111:757–67. [3] Stiassnie M. On the application of fractional calculus for the formulation of viscoelastic models. Applied Mathematics Modeling 1979;3:300–2. [4] Haneczok G, Weller M. A fractional model of viscoelastic relaxation. Materials Science and Engineering A 2004;370:209–12. [5] Ramirez LES, Coimbra CFM. A variable order constitutive relation for viscoelasticity. Annals of Physics (Leipzig) 2007;16:543–52. [6] Adolfsson K, Enelund M, Olsson P. On the fractional order model of viscoelasticity. Mechanics of Time-Dependent Materials 2005;9:15–34. [7] Distefano JN, Sackman JL. On the stability of initially imperfect, nonlinearly viscoelastic column. International Journal of Solids and Structures 1968;4:341–54. [8] Shirahatti US, Sinha SC. On the stability of perfect viscoelastic columns. Journal of Sound and Vibration 1994;174:57–68. [9] Touati D, Cederbaum G. Post buckling analysis of imperfect nonlinear viscoelastic columns. International Journal of Solids and Structures 1997;34:1751–60. [10] Chang WP. Creep buckling of nonlinear viscoelastic columns. Acta Mechanica 1986;60:199–215. [11] Vinogradov AM. Buckling of viscoelastic beam columns. AIAA Journal 1987;25:479–83. [12] Drawshi M, Cederbaum G. Stability of multiloaded viscoelastic nonlinear beams. Computers & Structures 1993;46(2):215–8. [13] Hewitt JS, Mazumdart J. Buckling of viscoelastic plates. AIAA Journal 1977;15:451–2. [14] Shalev D, Aboudi J. Postbuckling analysis of viscoelastic laminated plates using higher-order theory. International Journal of Solids and Structures 1991;27:1747–55. [15] Touati D, Cederbaum G. Post buckling of nonlinear viscoelastic imperfect laminated plates part II: structural analysis. Computers & Structures 1998;42:43–51. [16] Chen L-Q, Cheng C-J. Instability of nonlinear viscoelastic plates. Applied Mathematics and Computation 2005;162:1453–63. [17] Galucio AC, Deu J –F, Ohayon R. Finite element formulation of viscoelastic sandwich beams using fractional derivative operators. Computational Mechanics 2004;33:282–91. [18] Schmidt A, Gaul L. Finite element formulation of viscoelastic constitutive equations using fractional time derivatives. Nonlinear Dynamics 2002; 29(1–4):37–55. [19] Stankovic B, Atanackovic TM. Dynamics of a rod made of generalized Kelvin– Voigt model visco-elastic material. Journal of Mathematical Analysis and Applications 2002;268:550–63.
1048
J.T. Katsikadelis, N.G. Babouskos / Engineering Analysis with Boundary Elements 34 (2010) 1038–1048
[20] Rossikhin YuA, Shitikova MV. Analysis of damped vibrations of linear viscoelastic plates with damping modelled with fractional derivatives. Signal Processing 2006;86:2703–11. [21] Rossikhin YuA, Shitikova MV. Analysis of free non-linear vibrations of a viscoelastic plate under the conditions of different internal resonances. International Journal of Non-Linear Mechanics 2006;41:313–25. [22] Katsikadelis JT. Numerical solution of multi-term fractional differential ¨ Angewandte Mathematik und Mechanik equations. ZAMM Zeitschrift fur 2009;89(7):593–608. [23] Katsikadelis JT. The fractional wave-diffusion equation in bounded inhomogeneous anisotropic media. An AEM solution. In: Manolis GD, Polyzos D, editors. Advances in boundary element methods: a volume to honor professor dimitri beskos. Dordrecht , Netherlands: Springer Science; 2009. p. 255–276. [24] Babouskos NG, Katsikadelis JT. Nonlinear vibrations of viscoelastic plates of fractional derivative type : an aem solution. The Open Mechanics Journal 2010;4:8–20. [25] Katsikadelis JT. Nonlinear vibrations of viscoelastic membranes of fractional derivative type. In: Sapountzakis EJ, Aliabadi MH, editors. Advances in
[26] [27]
[28]
[29]
[30] [31]
boundary element techniques X, proceedings of the international conference on boundary element techniques BeTeq’09. UK: ECltd; 2009. p. 7–18. Podlubny I. Fractional differential equations. New York: Academic Press; 1999. Katsikadelis JT. The analog equation method—a powerful BEM-based solution technique for solving linear and nonlinear engineering problems. In: Brebbia CA, editor. Boundary element method XVI. Southampton: Computational Mechanics Publications; 1994. p. 167–82. Katsikadelis JT. A generalized Ritz method for partial differential equations in domains of arbitrary shape using global shape functions. Engineering analysis with boundary elements 2007;32:353–67. Katsikadelis JT. Special methods for plate analysis. In: Beskos DE, editor. Boundary element analysis of plates and shells. Springer-Verlag; 1991. p. 221–311. Katsikadelis JT. Boundary elements: theory and applications. Amsterdam– London: Elsevier; 2002. Katsikadelis JT, Babouskos N. The post-buckling analysis of plates. A BEM based meshless variational solution. Facta Universitatis, Series Mechanics, Automatic Control and Robotics 2007;6:113–8.
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
Post-buckling analysis of viscoelastic plates with fractional derivative models J.T. Katsikadelis n, N.G. Babouskos School of Civil Engineering, National Technical University of Athens, Athens GR-15773, Greece
a r t i c l e in f o
a b s t r a c t
Article history: Received 17 February 2010 Accepted 26 June 2010 Available online 23 July 2010
The post-buckling response of thin plates made of linear viscoelastic materials is investigated. The employed viscoelastic material is described with fractional order time derivatives. The governing equations, which are derived by considering the equilibrium of the plate element, are three coupled nonlinear fractional partial evolution type differential equations in terms of three displacements. The nonlinearity is due to nonlinear kinematic relations based on the von Ka´rma´n assumption. The solution is achieved using the analog equation method (AEM), which transforms the original equations into three uncoupled linear equations, namely a linear plate (biharmonic) equation for the transverse deflection and two linear membrane (Poisson’s) equations for the inplane deformation under fictitious loads. The resulting initial value problem for the fictitious sources is a system of nonlinear fractional ordinary differential equations, which is solved using the numerical method developed recently by Katsikadelis for multi-term nonlinear fractional differential equations. The numerical examples not only demonstrate the efficiency and validate the accuracy of the solution procedure, but also give a better insight into this complicated but very interesting engineering plate problem & 2010 Elsevier Ltd. All rights reserved.
Keywords: Plates Large deflections Buckling Analog equation method Boundary element method Viscoelasticity Fractional derivatives
1. Introduction The stability and the post-buckling response of structures made of viscoelastic materials have attracted the interest in recent years. Viscoelastic materials, such as polymers, exhibit both viscous and elastic behavior and are extensively used in many modern engineering applications like aircrafts, ships, buildings and space structures due to the advantage of light weight and high strength. The buckling of structures subjected to external compressive loads is initiated by small imperfections which are unavoidable in practice. In the case of elastic material, when the loads exceed a critical value, the structure becomes unstable and reaches a new stable position instantaneously. However, for viscoelastic materials, buckling is not a classical stability problem, since it is an evolution phenomenon. As the loads exceed a critical value, the structure becomes unstable and needs some time to reach its final stable position, which depends on the viscoelastic characteristics of the material and the value of the initial imperfection. Moreover, for some type of viscoelastic materials with initial elastic response and when the external loads remain within certain values creep buckling occurs. In this case the structure is unstable but the displacements increase slowly with time. The time needed for high deformation rates is large enough extending the useful life of the structure [1,2].
n
Corresponding author. E-mail addresses: [email protected] (J.T. Katsikadelis), [email protected] (N.G. Babouskos). 0955-7997/$ - see front matter & 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2010.07.003
There have been proposed various linear and nonlinear models for the mechanical behavior of these materials, such as Maxwell, Kelvin–Voigt, Zener and multi-parameter models. The constitutive equations can be presented in differential or in hereditary integral form. Recently, many researchers have shown that differential viscoelastic models with fractional derivatives are in better agreement with the experimental results than the integer derivative models [3–6]. Adolfsson et al. [6] proved that complicated multi-parameter integer viscoelastic models converge to fractional models with one single parameter. While the post-buckling of viscoelastic beams using integral or integer derivative models has been investigated by many authors [2,7–12] using analytic or approximate solution techniques, the literature on buckling of viscoelastic plates is rather limited. The reason is that the governing equations are linear or nonlinear evolution partial differential equations whose solution is difficult to obtain, especially for plates of arbitrary shape. Hewitt and Mazumdart [13] studied the stability of linear viscoelastic isotropic plates of arbitrary geometry using the method of constant deflection contours. Shalev and Aboudi [14] investigated the post-buckling response of laminated viscoelastic plates using Laplace transform and taking into account shear deformation. Touati and Cederbaum [15] studied the post-buckling of nonlinear viscoelastic laminated plates using the integral representation formulation. Chen and Cheng [16] studied the instability of nonlinear viscoelastic plates using constitutive equations in integral formulation. In most cases the plate problem is solved using Galerkin method in conjunction with a numerical integration scheme.
J.T. Katsikadelis, N.G. Babouskos / Engineering Analysis with Boundary Elements 34 (2010) 1038–1048
With regard to the fractional derivative models the viscoelastic problems treated are quite few. The reason is obvious. The equations describing the response of the systems are complicated linear or nonlinear fractional partial differential equations whose solution is very difficult to obtain. There have been recently some rather ad hoc solutions for the analysis of viscoelastic beams and ¨ plates. Grunwald formalism for the fractional order derivatives have been used by Galucio et al. [17], who analyzed viscoelastic sandwich beams with FEM and Schmidt and Gaul [18], who analyzed three dimensional viscoelastic bodies using FEM. Stankovic and Atanackovic [19] studied the vibrations of an axially compressed viscoelastic rod using the fractional Kelvin– Voigt model. Rossikhin and Shitikova [20,21] treated linear and nonlinear vibrations of rectangular viscoelastic plates with damping modeled with Riemann–Liouville fractional derivative using an approximate analytical solution. Recently Katsikadelis developed the AEM for solving linear and nonlinear fractional differential equations [22]. This method is general and has been already used to solve the linear fractional diffusion–wave equation in bounded inhomogeneous anisotropic bodies [23], the nonlinear dynamic response of viscoelastic plates and membranes [24,25]. In this paper this method is employed to study the post-buckling response of thin viscoelastic plates. Without excluding other models, the employed herein model is the generalized Kelvin–Voigt model with fractional order derivatives. The plate may have an arbitrary shape and be subjected to all types of conventional boundary conditions including also elastic support. For the two dimensional elasticity, the stress–strain relations can be reduced to the form [24] 38 9 ex > > = < 6 7 n 1 0 7 ey sy ¼ E 2 6 4 5> 1 n > > 1 n ; :t ; :g > 0 0 xy xy 2 2 38 a 9 1 n 0 D e > < c x > = ZE 6 0 7 6n 1 7 Dac ey þ 4 5 2 1n > 1n : Da g > ; 0 0 c xy 2 9 8 s = > < x>
2
1
n
1039
deformation under fictitious loads that are unknown in the first instance. Application of the D/BEM (Domain BEM) yields the initial value problem for the fictitious sources which is a system of nonlinear fractional ordinary differential equations. This system is solved utilizing the numerical method developed recently by Katsikadelis for multi-term linear and nonlinear fractional differential equations [22]. The post-buckling response of several viscoelastic plates is studied. The obtained numerical results demonstrate the efficiency of the developed solution method and validate its accuracy. Moreover the conclusions drawn from this investigation permit a better insight into this complicated but very interesting engineering plate problem.
2. Governing equations 2.1. The nonlinear plate problem We consider a thin elastic plate of uniform thickness h occupying the two dimensional multiply connected domain O of the xy plane with boundary G [Ki¼ 0 Gi (Fig. 1). The nonintersecting curves Gi (i¼0, 1, 2, y, K) may be piece-wise smooth. The boundary may be simply supported, clamped, free or elastically supported with transverse stiffness kT ðxÞ and rotational stiffness kR ðxÞ x : ðx,yÞ A G. The plate is subjected to inplane loads px, py (body forces). Moreover, along the movable edges inplane forces Nn , Nt may act. The von Ka´rma´n assumption for the kinematic relation is adopted. That is 1 2
1 2
ex ¼ u, x þ w, 2x ey ¼ v, y þ w, 2y gxy ¼ u, y þ v, x þw, x w, y
0
ð1Þ
where E, n are the elastic material constants, Z is the viscous parameter and Dac the Caputo fractional derivative of order a defined as 8" # Z t > 1 uðmÞ ðtÞ > > > , m1 o a o m d t < GðmaÞ a þ 1m 0 ðttÞ ð2Þ Dac uðtÞ ¼ > dm > > > m¼a : m uðtÞ dt where m is a positive integer. The advantage of this definition for the fractional derivative is that it permits the assignment of initial conditions which have direct physical significance [26]. Apparently, the classical derivatives result for integer values of a. For a ¼1 the conventional Kelvin–Voigt model is obtained while for a ¼0 the material is purely elastic with Young modulus E* ¼E(1+ Z). The governing equations result by taking the equilibrium of the plate element in a slightly deformed configuration and using the constitutive equations (1). They are three coupled nonlinear fractional partial differential equations in terms of displacements. The nonlinearity is due to nonlinear kinematic relations based on the von Ka´rma´n assumption. The concept of the analog equation is applied [27,28] to convert the original nonlinear fractional partial equations into three uncoupled linear equations, namely a linear plate (biharmonic) equation for the transverse deflection and two linear membrane (Poisson’s) equations for the inplane
ð3a; b; cÞ
where u ¼ uðx,y,tÞ, v ¼ vðx,y,tÞ are the membrane and w ¼ wðx,y,tÞ the transverse displacements. On the basis of the Eqs. (1) and (3) the stress resultants are written as 1 2 w, x þ nw, 2y Nx ¼ C u, x þ nv, y þ 2 1 ð4aÞ þ ZC Dac u, x þ nDac v, y þ Dac w, 2x þ nw, 2y 2 1 Ny ¼ C nu, x þ v, y þ nw, 2x þw, 2y 2 1 a þ ZC nDc u, x þDac v, y þ Dac nw, 2x þ w, 2y 2
x
Γ0
ð4bÞ
ss *
Nn px (x, y)
y
c
Γk
f hole py (x, y)
*
Nt
(Ω) corner ss Fig. 1. Plate geometry and supports (c¼clamped, ss ¼ simply supported, f¼ free).
1040
J.T. Katsikadelis, N.G. Babouskos / Engineering Analysis with Boundary Elements 34 (2010) 1038–1048
1n u, y þ v, x þ w, x w, y 2 1n a 1 Dc u, y þ Dac v, x þ Dac w, x w, y þ ZC 2 2
Nxy ¼ C
ð4cÞ
Mx ¼ D w, xx þ nw, yy ZD Dac w, xx þ nDac w, yy
My ¼ D nw, xx þw, yy ZD nDac w, xx þDac w, yy
ð5aÞ
The star designates a prescribed quantity. In Eq. (8), Vw is the equivalent elastic shear force, Mw is the elastic normal bending moment; Tw is the elastic twisting moment along the boundary and ½Twk is its discontinuity jump at corner k. The operators producing these quantities are given as [29] " !# @ 2 @ @2 @ k V ¼ D r þ ð1nÞ ð10aÞ @n @s @s@n @s
ð5bÞ "
Mxy ¼ Dð1nÞw, xy ZDð1nÞDac w, xy
ð5cÞ
where C ¼Eh/(1 n2) is the membrane stiffness and D ¼ Eh3 =12ð1n2 Þ the flexural stiffness of the plate. The governing equations result by taking the equilibrium of the plate element in a slightly deformed configuration. This yields @2 Mxy @2 Myy @2 Mxx @ Nx w, x þ Nxy w, y þ þ2 þ 2 @x @x@y @x @x2 @ Nxy w, x þNy w, y ¼ 0 þ @y
ð6aÞ
2
M ¼ D r ð1nÞ
T ¼ Dð1nÞ
!# @2 @ þ k @n @s2
ð10bÞ
! @2 @ k @s @s@n
ð10cÞ
where k ¼ k(s) is the curvature of the boundary and n, s the intrinsic boundary coordinates, i.e. the coordinates along the normal and tangential directions at a boundary point. 2.2. The linear plate problem
@Nx @Nxy þ þpx ¼ 0 @x @y
ð6bÞ
@Nxy @Ny þ þpy ¼ 0 @x @y
ð6cÞ
Using Eqs (4)–(6) we obtain the plate equations in terms of displacements in O as i. For the transverse deflection 4 4 Dr wþ ZDDac r wNx w, xx 2Nxy w, xy Ny w, yy þ px w, x þpy w, y ¼ 0 ð7aÞ ii. For the inplane deformation
1þn 2 r2 u þ u, x þ v, y , x þw, x w, xx þw, yy 1n 1n 1þ n 1þn 2 w, xy w, y þ ZDac r u þ u, x þv, y , x þ 1n 1n
2 1þn px ¼0 w, xx þ w, yy þ w, xy w, y þ þw, x 1n 1n Gh
4
Dr wlðNx w, xx þ2Nxy w, xy þ Ny w, yy Þ ¼ 0 Vw þ lðNn* w, n þ Nt* w, t Þ þkT w ¼ 0 Mw þkR w, n ¼ 0 ð7bÞ
1þn 2 ðu, x þv, y Þ, y þ w, y w, yy þ w, xx 1n 1n 1þ n 1þn w, xy w, x þ ZDac r2 vþ ðu, x þ v, y Þ, y þ 1n 1n
py 2 1þn ¼0 w, yy þ w, xx þ w, xy w, x þ þw, y 1n 1n Gh
r2 v þ
The plate problem is linearized if the stretching of the middle surface due to bending is neglected, that is if it is set w, 2x ¼ w, 2y ¼ w, x w, y ¼ 0 in Eq. (4). This implies that, although the inplane forces contribute to bending, they are not influenced by it. Thus the nonlinear terms in Eq. (7) vanish. The boundary conditions remain the same for both problems. In this case the two problems are uncoupled, hence the inplane forces Nx, Ny, Nxy are obtained from the solution of the Eq. (7b), (7c), which are solved independently. For linear elastic buckling Z ¼0 and the eigenvalue problem is derived as
ð7cÞ
ðkÞ kðkÞ T w ½Twk ¼ 0
or
w, n ¼ 0 or
or
in O
w¼0
on G
on G
wðkÞ ¼ 0
at corner point k
ð11Þ ð12aÞ ð12bÞ ð12cÞ
where l is the load parameter. The linear eigenvalue problem is also solved in this investigation, because the mode shapes (eigenvectors) of this problem are used as the initial imperfection for the post-buckling analysis. Besides, the linear elastic buckling loads are needed to compare the elastic with the viscoelastic plate response. The AEM solution for the linear and the nonlinear plate problems is presented below.
The associated boundary conditions result as Vw þ ZDac ðVwÞ þ Nn* w, n þ Nt* w, t þ kT w ¼ 0 Mw þ ZDa ðMwÞ þ k
R w, n
c
¼0
ZDa ½Tw
ðkÞ kðkÞ T w ½Twk
c
k
or ¼0
w,n ¼ 0 or
ðkÞ
w
or
w¼0
on G ¼0
on G
ð8aÞ 3. The AEM solution ð8bÞ
at corner point k ð8cÞ
Nn ¼ Nn*
or
un ¼ u*n
Nt ¼ Nt*
or
ut ¼ u*t
on G
ð8dÞ
The initial boundary value problem (7a), (8a)–(8c), (9a) for the response of the plate is solved using the AEM. The analog equation for the problem at hand is
r4 w ¼ bðx,tÞ, x ¼ fx,yg A O on G
ð8eÞ
Moreover, the viscoelastic plate is subjected to the initial conditions wðx,0Þ ¼ gðxÞ,
3.1. The plate problem
uðx,0Þ ¼ g1 ðxÞ,
vðx,0Þ ¼ g2 ðxÞ, in O
ð9a; b; cÞ
ð13Þ
where bðx,tÞ represents the time dependent fictitious load. Eq. (13) is a quasi-static equation, that is the time appears as a parameter, and it can be solved with the boundary conditions (8a)–(8c) at any instant t using the BEM. Thus, the solution at a
J.T. Katsikadelis, N.G. Babouskos / Engineering Analysis with Boundary Elements 34 (2010) 1038–1048
point xA O is obtained in integral form as Z Z wðx,tÞ ¼ w* bdO þ ðw* Vw þ w, n Mw* w* , n MwwVw* Þds O
X
1
ð14Þ
k
which for x A G yields the following two boundary integral equations for points where the boundary is smooth Z Z 1 wðx,tÞ ¼ w* bdO þ ðw* Vw þ w, n Mw* w* , n MwwVw* Þds 2 O G X w* ½Tww½Tw* ð15Þ
y
w*1 bdO þ
Z
M domain points
6 5
4
Fig. 3. Boundary and domain discretization.
ðw*1 Vw þ w, n Mw*1 w*1 , n MwwVw*1 Þds
O G X w*1 ½Tww½Tw*1 k
ð16Þ
k
in which w* ¼ w* ðx,yÞ, x,y A G, is the fundamental solution and w*1 its normal derivative at point x A G, i.e. 1 2 r ln r 8p
1 2 1 r ln r , n ¼ rr, n ð2 ln r þ 1Þ w*1 ¼ 8p 8p
w* ¼
ð17a; bÞ
n is the unit normal vector to the boundary at point x, whereas n is the unit normal vector to the boundary at the integration point y and r ¼ :xy: (see Fig. 2). Eqs. (15) and (16) can be used to establish the unspecified boundary quantities. They are solved numerically using the BEM. The boundary integrals are approximated using N constant boundary elements. Although the domain integrals can be approximated by radial basis functions and using the DRM to convert them to boundary integrals, maintaining thus the pure boundary character of the method [28], we approximate them with M linear triangular elements. This is more practical because the domain discretization is performed automatically using the Delaunay triangulation utilizing ready-touse subroutines. In this case, attention should be paid to place the nodal points of the triangles adjacent to the boundary on their sides, because the fictitious source is not defined on the boundary (see Fig. 3). Thus, after discretization and application of the boundary integral Eqs. (15) and (16) at the N boundary nodal points and Eq. (15) at the Nc corner points we obtain 8 9 8 9 > > =
where H, G are (2N+ Nc) (2N + Nc) known coefficient matrices originating from the integration of the kernel functions on the boundary elements; A is an (2N + Nc) M coefficient matrix originating from the integration of the kernel function on the domain triangular cells; w,wc ,w, n are the vectors of the N boundary nodal displacements, Nc corner displacements and N boundary nodal normal slopes, respectively, and V, R, M are the vectors of the N nodal values of effective shear force, Nc concentrated corner forces and N nodal values of the normal bending moment. Finally, b is the vector of the M nodal values of the fictitious source inside O. Eq. (18) constitutes a system of 2N + Nc equations for 4N +2Nc +M unknowns. Additional 2N + Nc equations are obtained from the boundary conditions. Thus, the boundary conditions (8a)–(8c), when applied at the N boundary nodal points and the Nc corner points yield the set of equations
a1 w þ a2 w, n þ a3 V ¼ 0 b1 w, n þ b2 M ¼ 0 c1 wc þc2 R ¼ 0 ð19a; b; cÞ where a1, a2, a3, b1, b2, c1, c2 are known coefficient matrices. Note that Eq. (19a) has resulted after approximating the derivative w,t in Eq. (8a) with a finite difference scheme. Eqs. (18) and (19) can be combined and solved for the boundary quantities w, wc , w, n , V, R, M in terms of the fictitious load b. Subsequently, these expressions are used to eliminate the boundary quantities from the discretised counterpart of Eq. (14). Thus we obtain the following representation for the deflection wðx,tÞ ¼
M X
bk ðtÞWk ðxÞ
xAO
ð20Þ
k¼1
The derivatives of wðxÞ at points x inside O are obtained by direct differentiation of Eq. (14). Thus, we obtain after elimination of the boundary quantities w, pqr ðx,tÞ ¼
M X
bk ðtÞWk , pqr ðxÞ,
p,q,r ¼ 0,x,y
xAO
ð21Þ
k¼1
t
Γ
N boundary points
x
k
k
Z
Nc corner points
G
w* ½Tww½Tw*
k
1 w, n ðx,tÞ ¼ 2
3
2
1041
n y
r = ||x - y||
r = ||x - y||
where Wk , pqr are known functions. Note that the above notation implies w,000 ¼w, w,0y0 ¼w,y, etc. 3.2. The plane stress problem
x x
y
a x
corner k Fig. 2. BEM notation.
Noting that Eqs. (7b) and (7c) are of the second order their analog equations are obtained using the Laplace operator. This yields
r2 u ¼ b1 ðx,tÞ r2 v ¼ b2 ðx,tÞ
ð22a; bÞ
The integral representation of the solution of Eq. (22a) is Z Z euðxÞ ¼ v* b1 dO ðv* qq* uÞds x A O [ G ð23Þ O
G
1042
J.T. Katsikadelis, N.G. Babouskos / Engineering Analysis with Boundary Elements 34 (2010) 1038–1048
in which q ¼u,n; v* ¼ ‘nr=2p is the fundamental solution to Eq. (22a) and q* ¼ v, *n its derivative normal to the boundary; r ¼ :yx: x A O [ G and y A G; e is the free term coefficient (e ¼1 if x A O, e ¼1/2 if x A G and e ¼0 if x= 2O [ G). Using the BEM with constant boundary elements and linear triangular domain elements and following the same procedure applied for the plate equation, we obtain the following representation for the inplane displacement u and its derivatives: u, pq ðx,tÞ ¼
M X
bð1Þ ðtÞUkð1Þ , pq ðxÞ þ k
k¼1
M X
bð2Þ ðtÞUkð2Þ , pq ðxÞ k
k¼1
þ U0 , pq ðxÞ, xA O, p,q ¼ 0,x,y
ð24Þ
Eq. (28) yields the buckling load parameter l and the buckling mode shapes. Buckling occurs due to an imperfection of the plate, which is given as a small initial deflection. In our examples the form of the first buckling mode of the linear elastic plate is taken as initial deflection. Subsequently, we increase the inplane loads and we study the response of the plate. When the inplane loads become greater than a critical value instability occurs, the displacements increase with time and the plate reaches a new stable position. The initial value problem Eqs. (26) and (27) is solved using a numerical method which has been developed recently by Katsikadelis for multi-term linear and nonlinear fractional differential equations [22].
Similarly, we obtain for the displacement v v, pq ðx,tÞ ¼
M X
bð1Þ ðtÞVkð1Þ , pq ðxÞ þ k
k¼1
M X
5. Examples
bð2Þ ðtÞVkð2Þ , pq ðxÞ k
k¼1
þ V0 , pq ðxÞ, x A O, p,q ¼ 0,x,y
ð25Þ
Ukð1Þ , Ukð2Þ , Vkð1Þ , Vkð2Þ , U0 , V0 are known functions; U0 , V0 result from the nonhomogeneous boundary conditions. Note that u,00 ¼u and v,00 ¼v.
4. The final step of the AEM Eqs. (21), (24) and (25) give the displacements w(x,t), u(x,t), v(x,t) and their derivatives provided that the three fictitious ð1Þ ð2Þ sources bðtÞ, b ðtÞ, b ðtÞ are first established. This is achieved by working as following. Collocating the fractional PDEs (7) at the M internal nodal points and substituting the expressions for the transverse deflection, Eq. (21), and the membrane displacements, Eqs. (24) and (25), we obtain the following system of 3M nonlinear fractional equations for bk ðtÞ, bð1Þ ðtÞ, bð2Þ ðtÞ, (k¼1, y, M) k k Fðb,b ,b ,Dac b,Dac b ,Dac b Þ ¼ P ð1Þ
ð2Þ
ð1Þ
ð2Þ
ð26Þ
where P is the vector of the external loads. Collocating the initial conditions Eq. (9) at the M internal nodal points and using Eqs. (21), (24) and (25), the initial conditions for the fictitious loads bk ðtÞ, bð1Þ ðtÞ, bð2Þ ðtÞ are obtained as k k n oT ð1Þ bð0Þ ¼ W1 g, b ð0Þ ¼ S1 g1 g2 þH1 , n oT ð2Þ þ H2 ð27Þ b ð0Þ ¼ S2 g1 g2 where W is N N known matrix, S1, S2 are N 2N known matrices and H1, H2 known vectors. The solution of the linear eigenvalue problem, Eqs. (11) and (12), is obtained following a similar procedure. Thus, collocating Eq. (11) at the M internal nodal points and substituting the expressions for the transverse deflection and its derivatives, Eq. (21), we obtain the following eigenvalue problem: ðKlGÞb ¼ 0
ð28Þ
Example 1. The post-buckling response of the square plate of thickness h¼0.1 m is investigated (Fig. 4). All edges are simply supported and movable in the plane of the plate. The plate is subjected to an inplane line load along two opposite sides. The results were obtained with N ¼196 boundary points and M ¼97 internal points resulting from 145 triangular cells (Fig. 5). The parameters of the viscoelastic material are E¼ 21 106 kN/m2, n ¼0.25, Z ¼0.5. The computed linear buckling load of the elastic plate is Pel ¼4681 kN/m (exact Pel ¼4606 kN/m). The initial deflection has the form of the first buckling mode of the elastic plate. First we analyze a viscoelastic plate with the conventional Kelvin–Voigt model (a ¼ 1). Fig. 6 presents the time history of the
w = Mw = 0 Nn = Nt = 0 4
w = Mw = 0
w = Mw = 0
Nt = 0
4 Nt = 0 Nn = -P
Nn = -P w = Mw = 0 Nn = Nt = 0
Fig. 4. Geometry and boundary conditions of plate in Example 1.
4
3
2
where K, G are known N N matrices given as Kik ¼ Ddik
i Gik ¼ Nxi Wki , xx þ 2Nxy Wki , xy þ Nyi Wki , yy
ð29Þ
with i, k¼1, y, M, dik the Kronecker delta and Wki , xx ¼ Wk , xx ðxi Þ, Wki , xy ¼ Wk , xy ðxi Þ, Wki , xx ¼ Wk , yy ðxi Þ; b is a vector with the nodal values of the fictitious load, which now does not depend on time. The inplane forces Nx, Ny, Nxy are established by solving independently the linear plane stress problem (see Section 2.2) using the direct BEM [30]. The solution of the eigenvalue problem
1
0
0
1
2
3
4
Fig. 5. Boundary and domain discretization of plate in Example 1.
J.T. Katsikadelis, N.G. Babouskos / Engineering Analysis with Boundary Elements 34 (2010) 1038–1048
0.014
-4200 P = 4800
P = 4700 P = 4650
0.012 0.01
P = 5000
P = 5200
-4600
0.008
Nx
w/h
1043
0.006
-5000
0.004 0.002 0 0
20
40 60 time (sec)
80
-5400
100
0
50
100
150
200
time (sec)
Fig. 6. Time history of the deflection at the center of the plate for a ¼ 1 just before and after the elastic buckling load.
Fig. 9. Time history of the membrane force Nx at point (1.4, 2) for a ¼ 1.
5400
1 P = 4800
P = 5000
P = 5200
5300
0.8
5200 5100
w/h
load P
0.6
0.4
5000 t = 20 sec t = 40 sec t = 60 sec t = 100 sec t = 160 sec Elastic with E
4900 4800
0.2
4700 4600
0 0
50
100 time (sec)
150
200
0
0.2
0.4
0.6
0.8
1
w/h
Fig. 7. Time history of the deflection at the center of the plate for a ¼1.
Fig. 10. Load–deflection curves at the center of the plate at different values of time for a ¼ 1.
150
observed that viscoelastic plates become unstable when the inplane load reaches the linear elastic buckling load and displacements increase with time until they reach a final value, in contrast to the elastic plates, which reach the final stable position immediately. The viscoelastic plate has a softer performance in buckling than the elastic plate. Fig. 10 presents the load–deflection curves for different values of time and the results are compared with the post-buckling path of the elastic plate [31]. It is observed that after certain time the viscoelastic post-buckling path approaches the elastic post-buckling path. Next, we study the response of the generalized Kelvin–Voigt model. Fig. 11 presents the time history of the deflection at the center of the plate for various values of a and for two cases of the inplane line load, just before and after the linear buckling load of the elastic plate. It is observed that in all cases the displacements increase as the inplane load becomes greater than the elastic buckling load, which means that the plate becomes unstable. However, for a o1 the displacements increase slowly and only when the load becomes greater than another critical value P 0 (P 0 oPel) the displacements increase rapidly. Such a gradual buckling, which occurs for loads Pel oPoP 0 , is known in literature as creep buckling and it is also observed in integer derivative viscoelastic models with initial elastic response [1,2]. Fig. 12 presents the time history of the
P = 4800
P = 5000
P = 5200
Mx
100
50
0 0
50
100
150
200
time (sec) Fig. 8. Time history of the bending moment at the center of the plate for a ¼ 1.
deflection at the center of the plate for two values of the inplane load P just before and after the buckling load Pel of the elastic plate. Figs. 7–9 present the time history of the deflection and the stress resultants for various values of the inplane load. It is
1044
J.T. Katsikadelis, N.G. Babouskos / Engineering Analysis with Boundary Elements 34 (2010) 1038–1048
6.2
7500
x 10-3
7400
α = 0.8 α = 0.5
6
7300 load P
α = 0.2 w/h
5.8
7200 7100
5.6
7000
α = 0.2
5.4
α = 0.5
6900
α = 0.8 0
t = 0.2 sec t = 0.5 sec t = 2 sec t = 50 sec elastic with E* = (1+η)E
2
4
0
6 time (sec)
8
0.1
0.2
0.3
0.4
10
0.7
7500
α=1 α = 0.2 α = 0.01 Elastic with E Elastic with E* =(1+η)E
7000 P = 7000
6500 load P
P = 7100 P = 7300 P = 7400
0.6
0.6
Fig. 14. Load–deflection curves at the center of the plate at different values of time for a ¼ 0.01.
Fig. 11. Time history of the deflection at the center of the plate for various values of a just before and after the elastic buckling load (solid line: P ¼ 4700, dashed line: P¼ 4650).
0.8
0.5
w/h
6000
w/h
5500 0.4
5000 4500
0.2
0
0 0
0.5
1 time (sec)
1.5
2
0.2
0.4 w/h
0.6
0.8
Fig. 15. Load–deflection curves for various values of the order a of the fractional derivative and comparison with the elastic post-buckling paths.
y
Fig. 12. Time history of the deflection at the center of the plate for a ¼0.01.
w = Mw = 0 u n = ut = 0 Nn = -P Nt = 0 w = Mw = 0
7000
x
load P
6500
6000 t = 5 sec t = 10 sec t = 20 sec t = 30 sec t = 100 sec t > 1000 sec
5500
5000 0
0.2
0.4
0.6
0.8
1
1.2
P
1.4
w/h Fig. 13. Load–deflection curves at the center of the plate at different values of time for a ¼ 0.2.
2 5 Fig. 16. Geometry and boundary conditions of plate in Example 2.
J.T. Katsikadelis, N.G. Babouskos / Engineering Analysis with Boundary Elements 34 (2010) 1038–1048
1045
4500
deflection at the center of the plate of various values of the inplane load for a ¼0.01. Figs. 13 and 14 present the load–deflection curves for different values of time for fractional order derivative a ¼0.2
4000
load P
5 3500 t = 2 sec
3000
2.5
t = 10 sec t = 100 sec
2500
t = 200 sec
0
t > 500 sec 2000 0
-2.5
0.1
0.2
0.3
0.4
0.5
0.6
w/h Fig. 20. Load–deflection curves at point r¼ 3.5 for various values of time (a ¼0.2).
-5 -5
-2.5
0
2.5
5
5000
Fig. 17. Boundary and domain discretization annular plate in Example 2.
4800
x
load P
1.3
10-3
1.2
t = 0.4 sec 4400
α = 0.5 α=1
t = 1 sec t = 2 sec
α = 0.2
t = 10 sec
4200
w/h
1.1
4600
t > 50 sec
1
4000 0
α = 0.5 5
10 time (sec)
15
20
Fig. 18. Time history of deflection at point r ¼3.5 for various values of a and for two cases of the inplane load, just before and after the elastic buckling load (solid line: P¼ 2250, dashed line: P ¼2200).
load P
3500
3000 t = 1sec t = 10sec
2500
t = 20 sec t > 100 sec 2000 0
0.1
0.2
0.3
0.4
Fig. 21. Load–deflection curves at point r ¼3.5 for various values of time (a ¼0.01).
0.8 0
0.2 w/h
α = 0.2
α=1
0.9
0.1
0.3
0.4
0.5
0.6
w/h Fig. 19. Load–deflection curves at point r ¼3.5 for various values of time (a ¼ 1).
and 0.01, respectively. It is observed that the critical value P 0 , above which high deformation rates appear, increases as the order of the fractional derivative decreases, and for a-0 the viscoelastic postbuckling response approaches the elastic post-buckling response with a modulus E*¼ (1+ Z)E. Finally, Fig. 15 presents the load– deflection curves after considerable time has elapsed. It is observed that for a ¼0.2 the critical load is P 0 ¼5400 and for a ¼0.01 is * ¼ 7021 of the P 0 ¼6950 very close to the elastic buckling load Pel plate with modulus E*. Example 2. The post-buckling response of the annular plate of Fig. 16 is investigated. The boundary conditions and the loading are shown in Fig. 16. The plate is simply supported at both boundaries. The inplane displacements are restrained along the external boundary. The results were obtained with N¼300 boundary points and M¼130 internal points resulting from 190 triangular cells (Fig. 17). The parameters of the viscoelastic material are E¼ 11 106 kN/m2, n ¼0.25, Z ¼1. The thickness of the plate is h¼0.1 m. The buckling load of the elastic plate is Pel ¼2224 kN/m (a FEM solution with 2000 elements yields Pel ¼ 2195 kN/m). Fig. 18 presents the time history of the deflection at point r¼3.5 for various values of a and for two cases of the inplane line load, just before and after the buckling load of the elastic plate. It is observed that in all cases the plate becomes unstable as the inplane load becomes greater than the elastic buckling load. For a o1 creep buckling occurs. Figs. 19–21
1046
J.T. Katsikadelis, N.G. Babouskos / Engineering Analysis with Boundary Elements 34 (2010) 1038–1048
present the load–deflection curves at point r¼3.5 for different values of time and for a ¼1, 0.2, 0.01. Finally, Fig. 22 presents the load–deflection curves for various values of the order of fractional derivative after considerable time has elapsed. The results are compared with the post-buckling path of the elastic plate with elastic constants E and E*¼E(1+ Z).
2
1.5 Example 3. The post-buckling response of the rectangular plate of Fig. 23 is investigated. The plate is clamped along side 1–4, simply supported movable along sides 1–2 and 3–4 and free to move along side 2–3. The plate is compressed along sides 1–2 and 3–4 with an inplane line load P. The results were obtained with N ¼200 boundary points and M ¼93 internal points resulting from 136 triangular cells (Fig. 24). The thickness of the plate is h¼0.01 m. The parameters of the viscoelastic material are E¼29 106 kN/m2, n ¼0.3, Z ¼1. The buckling load of the elastic plate is Pel ¼36.2 kN/m (a FEM solution with 200 elements yields
1
0.5
5000 4500 α=1 α = 0.5 α = 0.2 α = 0.01 Elastic with E Elastic with E* = E(1+η)
load P
4000 3500
0 0
0.5
1
Fig. 24. Boundary and domain discretization of plate in Example 3.
3000
50 2500
α=1 2000 0
0.1
0.2
0.3
0.4
α = 0.3
0.5
45
Fig. 22. Load–deflection curves at point r¼ 3.5 for various values of the order a of the fractional derivative and comparison with the elastic post-buckling paths.
Elastic
load P
w/h
40
y Nn = -P Nt = 0 w = Mw = 0
35 0
1
4
1
1.5
2
w/h Fig. 25. Load–deflection curves for two cases of the order a of the fractional derivative and comparison with the elastic post-buckling path.
3
un = ut = 0 2 w = w,n = 0
Nn = Nt = 0 Vw = Mw = 0
1
0.5
2
x
35.44 kN/m). Fig. 25 presents the load–deflection curves at point (1, 1) for two cases of the order a ¼1, 0.3, after considerable time has elapsed and the results are compared to the post-buckling path of the elastic plate. For a ¼0.3 the critical load is P 0 ¼ 39 kN/m, above which the displacements of the plate increase rapidly. Fig. 26a and b presents the time history of the deflection at point (1, 1) and the inplane displacement at point (0.16, 1) for various values of the initial imperfection (P ¼40 kN/m).
6. Conclusions
Nn = -P Nt = 0 w = Mw = 0 Fig. 23. Geometry and boundary conditions of plate in Example 3.
The post-buckling response of thin plates made of linear viscoelastic material has been investigated in this work. The behavior of the viscoelastic material is described by the
J.T. Katsikadelis, N.G. Babouskos / Engineering Analysis with Boundary Elements 34 (2010) 1038–1048
0
2
x 10-5
w0/h = 0.07
-0.2
1.5
w0/h = 0.007 -0.4
1
w0/h = 0.0007
u
w/h
1047
w0/h = 0.07
-0.6
w0/h = 0.007
0.5
-0.8
w0/h = 0.0007
-1
0 0
20
40
60 80 time (sec)
100
120
140
0
20
40
60 80 time (sec)
100
120
140
Fig. 26. Time history of the deflection w at point (1, 1) and the inplane displacement u at point (0.16, 1) for various values of the initial imperfection (P ¼ 40 kN/m).
generalized Kelvin–Voigt model with fractional order derivative. The solution of the derived governing equations is achieved using the AEM. This method converts the coupled nonlinear fractional PDEs, describing the response of the viscoelastic plate, into uncoupled linear fractional PDEs, which are subsequently solved by the D/BEM. The discretization procedure yields an initial value problem for simultaneous ordinary FDEs, which is solved by a new method for nonlinear multi-term fractional differential equations developed recently by Katsikadelis. The influence of the viscoelastic material on the post-buckling response of the plate has been studied. The conclusions that can be drawn from this study are summarized as:
(a) The viscoelastic plate becomes unstable when the inplane load reaches the elastic buckling load, which is obtained by the linear eigenvalue problem. Then the post-critical response is studied by gradually increasing the inplane load. (b) While the post-buckling response for an elastic material is initiated when the inplane load reaches the critical value, for a viscoelastic material the post-buckling response starts after a certain time, which depends on the viscoelastic parameters and the value of the initial imperfection. (c) For the conventional Kelvin–Voigt model (a ¼1) the viscoelastic post-buckling path converges to the elastic postbuckling path (Z ¼0). Practically, after certain time, which depends on Z and the value of the initial imperfection, the elastic and viscoelastic post-buckling paths become identical. (d) For the generalized Kelvin–Voigt model (0o a o1) creep buckling occurs. For certain values of the external load Pel oPoP 0 the plate becomes unstable but the displacements increase slowly. When the external load exceeds the critical value P 0 high deformation rates appear and the plate reaches a new stable position. For 0 o a o1, after certain time, the viscoelastic post-buckling path lies between the elastic paths with moduli E and E* ¼E(1+ Z). (e) As a-0 the viscoelastic post-buckling response of the plate approaches an elastic post-buckling response with a modulus E* ¼E(1+ Z). (f) The fractional type viscoelastic model with 0 o a o1 can be visualized as a model of combined springs and dashpots,
whose contribution depends on the value of a. This is reflected in the creep buckling response of the analyzed structures, as it happens with more complicated integer derivative models with initial elastic response.
References [1] Flugge W. Viscoelasticity. Blaisdell Publishing Company; 1967. [2] Vinogradov AM. Nonlinear effects in creep buckling analysis of columns. Journal of Engineering Mechanics 1985;111:757–67. [3] Stiassnie M. On the application of fractional calculus for the formulation of viscoelastic models. Applied Mathematics Modeling 1979;3:300–2. [4] Haneczok G, Weller M. A fractional model of viscoelastic relaxation. Materials Science and Engineering A 2004;370:209–12. [5] Ramirez LES, Coimbra CFM. A variable order constitutive relation for viscoelasticity. Annals of Physics (Leipzig) 2007;16:543–52. [6] Adolfsson K, Enelund M, Olsson P. On the fractional order model of viscoelasticity. Mechanics of Time-Dependent Materials 2005;9:15–34. [7] Distefano JN, Sackman JL. On the stability of initially imperfect, nonlinearly viscoelastic column. International Journal of Solids and Structures 1968;4:341–54. [8] Shirahatti US, Sinha SC. On the stability of perfect viscoelastic columns. Journal of Sound and Vibration 1994;174:57–68. [9] Touati D, Cederbaum G. Post buckling analysis of imperfect nonlinear viscoelastic columns. International Journal of Solids and Structures 1997;34:1751–60. [10] Chang WP. Creep buckling of nonlinear viscoelastic columns. Acta Mechanica 1986;60:199–215. [11] Vinogradov AM. Buckling of viscoelastic beam columns. AIAA Journal 1987;25:479–83. [12] Drawshi M, Cederbaum G. Stability of multiloaded viscoelastic nonlinear beams. Computers & Structures 1993;46(2):215–8. [13] Hewitt JS, Mazumdart J. Buckling of viscoelastic plates. AIAA Journal 1977;15:451–2. [14] Shalev D, Aboudi J. Postbuckling analysis of viscoelastic laminated plates using higher-order theory. International Journal of Solids and Structures 1991;27:1747–55. [15] Touati D, Cederbaum G. Post buckling of nonlinear viscoelastic imperfect laminated plates part II: structural analysis. Computers & Structures 1998;42:43–51. [16] Chen L-Q, Cheng C-J. Instability of nonlinear viscoelastic plates. Applied Mathematics and Computation 2005;162:1453–63. [17] Galucio AC, Deu J –F, Ohayon R. Finite element formulation of viscoelastic sandwich beams using fractional derivative operators. Computational Mechanics 2004;33:282–91. [18] Schmidt A, Gaul L. Finite element formulation of viscoelastic constitutive equations using fractional time derivatives. Nonlinear Dynamics 2002; 29(1–4):37–55. [19] Stankovic B, Atanackovic TM. Dynamics of a rod made of generalized Kelvin– Voigt model visco-elastic material. Journal of Mathematical Analysis and Applications 2002;268:550–63.
1048
J.T. Katsikadelis, N.G. Babouskos / Engineering Analysis with Boundary Elements 34 (2010) 1038–1048
[20] Rossikhin YuA, Shitikova MV. Analysis of damped vibrations of linear viscoelastic plates with damping modelled with fractional derivatives. Signal Processing 2006;86:2703–11. [21] Rossikhin YuA, Shitikova MV. Analysis of free non-linear vibrations of a viscoelastic plate under the conditions of different internal resonances. International Journal of Non-Linear Mechanics 2006;41:313–25. [22] Katsikadelis JT. Numerical solution of multi-term fractional differential ¨ Angewandte Mathematik und Mechanik equations. ZAMM Zeitschrift fur 2009;89(7):593–608. [23] Katsikadelis JT. The fractional wave-diffusion equation in bounded inhomogeneous anisotropic media. An AEM solution. In: Manolis GD, Polyzos D, editors. Advances in boundary element methods: a volume to honor professor dimitri beskos. Dordrecht , Netherlands: Springer Science; 2009. p. 255–276. [24] Babouskos NG, Katsikadelis JT. Nonlinear vibrations of viscoelastic plates of fractional derivative type : an aem solution. The Open Mechanics Journal 2010;4:8–20. [25] Katsikadelis JT. Nonlinear vibrations of viscoelastic membranes of fractional derivative type. In: Sapountzakis EJ, Aliabadi MH, editors. Advances in
[26] [27]
[28]
[29]
[30] [31]
boundary element techniques X, proceedings of the international conference on boundary element techniques BeTeq’09. UK: ECltd; 2009. p. 7–18. Podlubny I. Fractional differential equations. New York: Academic Press; 1999. Katsikadelis JT. The analog equation method—a powerful BEM-based solution technique for solving linear and nonlinear engineering problems. In: Brebbia CA, editor. Boundary element method XVI. Southampton: Computational Mechanics Publications; 1994. p. 167–82. Katsikadelis JT. A generalized Ritz method for partial differential equations in domains of arbitrary shape using global shape functions. Engineering analysis with boundary elements 2007;32:353–67. Katsikadelis JT. Special methods for plate analysis. In: Beskos DE, editor. Boundary element analysis of plates and shells. Springer-Verlag; 1991. p. 221–311. Katsikadelis JT. Boundary elements: theory and applications. Amsterdam– London: Elsevier; 2002. Katsikadelis JT, Babouskos N. The post-buckling analysis of plates. A BEM based meshless variational solution. Facta Universitatis, Series Mechanics, Automatic Control and Robotics 2007;6:113–8.