- Email: [email protected]
A modified variational iteration method for the analysis of viscoelastic beams
Accepted Manuscript
A modified variational iteration method for the analysis of viscoelastic beams Olga Martin PII: DOI: Reference:
S0307-904X(16)30212-8 10.1016/j.apm.2016.04.011 APM 11132
To appear in:
Applied Mathematical Modelling
Received date: Revised date: Accepted date:
25 April 2015 4 March 2016 13 April 2016
Please cite this article as: Olga Martin , A modified variational iteration method for the analysis of viscoelastic beams, Applied Mathematical Modelling (2016), doi: 10.1016/j.apm.2016.04.011
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Highlights
CR IP T
An algorithm based on VIM and Laplace transformation is created for the dynamic study of viscoelastic beams. It is avoided the definition of initial approximation in the using of VIM, which is a difficult step for nonlinear problems. There are no errors, which usually accompany the methods based on the discretization of time interval.
AC
CE
PT
ED
M
AN US
ACCEPTED MANUSCRIPT
A modified variational iteration method for the analysis of viscoelastic beams
CR IP T
Professor Olga Martin Applied Sciences Faculty Polytechnic University of Bucharest [email protected]
Abstract
AN US
Based on a constitutive law in a hereditary integral form, a mathematical model for dynamic analysis of the isotropic linear viscoelastic beams is presented. To solve the governing equation for these structures subjected to a distributed load is created an accuracy and computational efficiency algorithm that uses the Galerkin’s method and a modified form of the variational iteration method (VIM) for time-domain equations. Numerical results for both quasi-static and dynamic analysis are presented and these will be accompanied by the graphical representations.
Key Words: viscoelastic beam, Euler-Bernoulli beam theory, Galerkin method, variational iteration
M
method, Laplace transform, correspondence principle.
1 Introduction
ED
Viscoelastic materials, such as polymers, rubbers and metals at elevated temperature, are frequently present today in the engineering design due to their ability to dampen out the structural vibrations. There
PT
are remarkable theoretical studies on viscoelastic materials such as those of Flugge [8], Chritensen [4], Reddy [19], Kennedy [13], Krevelen [12], Cristescu [7] and Crawford [6]. Mathematical models, accompanied by the numerical methods were proposed to predict, as accurately, the response of
CE
viscoelastic components which depend on their hereditary loading and its service conditions [1], [3], [10], [14], [16], [21] - [24]. The methods of analysis of viscoelastic structures that occur often in the literature
AC
are the difference method and finite element method [5], [11], [12], [17], [20], [23], [24]. In this paper, the governing equation is determined for a simply supported viscoelastic beam under a
uniform distributed load using Euler-Bernoulli beam theory. This equation is accompanied by a constitutive law presented in a hereditary integral form. In order to obtain the quasi-static exact solution (i.e. the solution ignoring inertia effects), the correspondence principle was used [15], [24]. This principle relates mathematically the solution of a linear, viscoelastic boundary value problem to an analogous problem of an elastic body of the same geometry and under the same initial boundary conditions. It has to
ACCEPTED MANUSCRIPT
be mentioned that not all problems can be solved by this principle, but only those for which the boundary conditions do not vary with the time. Then, applying d’Alembert principle, this structure has been analyzed in the dynamic case by an algorithm that approaches the Galerkin’s method for the spatial domain and a modified variational iteration method, [9], for the temporal domain. The theoretical results
2 Mathematical model
CR IP T
have been applied both for quasi-static and dynamic responses of an isotropic linear viscoelastic beam.
Let us consider a simply supported beam with the span L and transverse section S, subjected to uniformly distributed forces 𝑝̅ . In Euler’s theory, the geometrical equation is given by
2 w( x, t ) z x2
(1)
AN US
z
where x and z represent the axial and transverse coordinates, ρ is the curvature radius, w(x,t) – the deflection and ε- the strain.
M
A constitutive law in the hereditary integral form, [22], [25] is proposed: t
( x, t ) E (0) ( x, t )
(2)
ED
0
dE (t ) ( x, ) d d
where E(t) is the relaxation modulus for the beam material and σ the stress corresponding to the strain ε.
CE
PT
Accordance with the d’Alembert’s principle, the governing equation is of the form
2 M 2w p A x2 t2
(3)
AC
where the bending moment M(x,t) is written as
M ( x, t ) ( x, t ) z dydz
(4)
S
Using (2) the bending moment in a beam of constant cross section is equal to
2 w( x, t ) dE (t ) 2 w( x, ) I d 0 d x2 x2 t
M ( x, t ) E (0) I
(5)
ACCEPTED MANUSCRIPT
where I is the area moment of inertia given by ∬𝑆 𝑧 2 𝑑𝑦𝑑𝑧. Thus, after the substitution of (5) in (3), the following motion equation in terms of the transverse deflection is obtained:
2 w( x, t ) 4 w( x, t ) dE (t ) 4 w( x, ) E ( 0 ) I I d p 0 d t 2 x4 x4 t
3 Mixed method for solving the motion equation
CR IP T
A
(6)
To solve the equation (6) for a simply supported beam subjected to a uniformly distributed loading 𝑝̅ , we express the transverse deflection w by an expansion of the form n
wn ( x, t ) ai (t ) i ( x)
(7)
AN US
i 1
where φi(x) is the i th shape function and ai (t) is the corresponding time-dependent amplitude. For the spatial domain will be used the Galerkin’s method and then, for the time domain will be proposed a new variational iteration method accompanied by the techniques of the Laplace transform.
M
The shape functions are chosen to be linearly independent, orthonormal L
1 i j 0 i j
ED
i ( x) j ( x) dx 0
(8)
PT
and must satisfy all boundary condition for the convergence of Galerkin’s method. Although wn satisfies the boundary conditions, in generally, does not satisfy equation (6). If the expansion
CE
(7) is substituted into (6) will result the residual function
AC
n n Rn ( x, t ) A ai(t ) i ( x) E (0) I ai (t ) i( 4) ( x) i 1 i 1 t
I 0
dE (t ) n ai ( ) i( 4) ( x) d p d i 1
(9)
Now consider the shape functions of the form
i ( x)
2 i x sin , i 1, 2,...., n L L
(10)
ACCEPTED MANUSCRIPT
for which
i i ( x ) i i ( x ) L 4
i( 4) ( x)
(11)
The Galerkin’s method requires that the residual to be orthogonal to each of the shape functions, so
R ( x, t ) n
j
(x) d x d t 0, j 1, 2,..., n
(12)
CR IP T
where [0, L] [0, t ] . This conditions leads to n equations verified by the functions aj(t):
dE (t ) a j ( ) d j I E (0)a j (t ) p j ( x)dx d 0 0 t
L
Aaj (t ) j I or
dE (t ) a j ( ) d j a j (t ) Pj d 0 t
AN US
aj (t ) j where
(13)
(14)
4 L jI j I E (0) 1 j , , , P p j ( x)dx j j j A A A 0 L
j
d kaj
L
t 0
0
ED
dt k
M
The initial conditions are the following
k wn ( x,0) j ( x)dx, k 0, 1, 2,.... t k
(15)
The functions aj(t) are determined independent of one another by (14) – (15). Finally, the value of the
PT
transverse deflection w(x,t) will be found by (7).
CE
3. 1. The variational iteration method (VIM) To obtain the solution of the integral-differential equation (14) will rewrite it using the linear operator L
AC
and the nonlinear operator N:
where
La j (t ) Na j (t ) F j 0
(16)
La j (t ) a j (t ) j a j (t ) dE (t ) a j ( )d , F j Pj d 0 t
N a j (t ) j
According to the variational iteration method (VIM), [9] the correction functional in t direction is defined as:
ACCEPTED MANUSCRIPT
t
a j ,m1 (t ) a j ,m (t ) ( )( La j ,m ( ) Na~ j ,m ( ) F j )d
(17)
0
To obtain the successive approximations a j ,m1 will determine the Lagrange multiplier λ, using the variational techniques. The function 𝑎̃𝑗,𝑚 has a restricted variation in t direction, which means
CR IP T
𝛿𝑎̃𝑗,𝑚 = 0, [26]. Remind that a j ,m (0) 0, u j ,m (0) 0 and F 0. Therefore, it will get, via the integration by parts, that
dE ( ) ~ a j ,m1 (t ) a j ,m (t ) ( ) a j,m ( ) j a j ,m ( )d j a j ,m ( ) F j d d 0 0 t
t 1 - ( ) a j ,m (t ) ( ) t a j ,m (t ) ( ) j ( ) a j ,m ( ) d 0 t 0
So, λ(ξ) verifies the following conditions
( )
t
0, 1 ( )
t
M
and the Lagrange multiplier will be identified as
ED
( )
AN US
1
j
(18)
0, ( ) j ( ) 0
(19)
(20)
sin j t
The iteration formula (17) becomes
dE ( ) a j ,m1 (t ) a j ,m (t ) sin j t a j,m ( ) j a j ,m ( ) d j a j ,m ( ) F j d d j 0 0 t
CE
PT
1
m = 0, 1,…
AC
with the initial conditions:
d k a j ,m (t ) dt k
t 0
0, m, k 0, 1, 2,....
Applying the integration by parts twice for the first term of the integral and (19), will get
(21)
ACCEPTED MANUSCRIPT
a j ,m 1 (t ) a j ,m (t )
1
j
sin j ( t ) a j ,m ( ) t0 cos j ( t ) a j ,m ( ) t0
t
j
t
1
dE ( t ) d j a j ,m ( ) sin j t d sin ( t ) a ( ) d j j , m d 0 j 0 0 j a j ,m ( ) sin j t d
j
t
F
j
Finally, (21) becomes
j
dE ( t ) 1 sin ( t ) a j ,m ( ) d d j d j 0 j 0 t
t
F
j
sin j ( t ) d
0
AN US
a j ,m1 (t )
sin j ( t ) d
0
CR IP T
0
t
(22)
In the case of an exponential form for the relaxation module: t
E (t ) m1 m2 e , m1 , m2 , 0 the iteration formula (22) becomes:
(23)
P e d j 1 cos( j t ) a j ,m1 (t ) sin ( t ) a ( ) d j j , m 0 j j 0
j m2 .
M
(24)
PT
where h j
t
ED
hj
CE
3. 2. Solution of the governing equation for the amplitude 𝒂𝒋 (𝒕) In the case of nonlinear problems, the iterative formula (22) defines a slow process, requiring a longer period of time, even if we use a modern computing technique. On the other hand, the sequence
AC
{𝑎𝑗,𝑚 (𝑡)}𝑚 defined by (21) or (24), converges to the exact solution of the equation (14), if the initial approximation 𝑎𝑗,0 (𝑡) verifies the initial conditions. In this paper, we propose to remove these inconveniences by associating Laplace transform techniques to MIV. Let us now consider that the sequence {𝑎𝑗,𝑚 (𝑡)}𝑚 converges to the solution 𝑎𝑗 (𝑡) of the equation (14). Therefore, the iterative process (22) ceases when 𝑎𝑗,𝑚 (𝑡) = 𝑎𝑗,𝑚+1 (𝑡) with an error less than a value 𝜀̃ very
ACCEPTED MANUSCRIPT
small and the solution 𝑎𝑗 (𝑡) will then be equal to
a j (t ) a j ,m (t ) a j ,m1 (t )
(25)
The formula (22) will become of the form
j
dE ( t ) Pj sin ( t ) a j ( )d d 1 cos( j t ) j d j j 0 0 t
(26)
CR IP T
a j (t )
Theorem. If the function 𝑎𝑗 : [0, 𝑇] → 𝐑, 𝐷 = [0, 𝑇], 𝑎𝑗 ∈ 𝐶 2 (𝐷), verifies the integral equation (26), then it is the solution of the integral-differential equation (14).
0
dE ( t ) a j ( )d . Then, the derivatives of the function a j (t ) will be d
AN US
□ Let be f ( )
t
a j (t ) j cos ( j ( t )) f ( )d Pj 0 t
sin ( j t )
j
a j (t ) j j sin ( j ( t )) f ( )d j f (t ) Pj cos( j t )
M
0
Introducing this last derivative in the left side of the equation (14) is found that
t
ED
j j sin ( j ( t )) f ( )d j 0
0
j t dE ( t ) a j ( )d j sin d j 0
PT
j
0
dE ( t ) a j ( )d Pj cos( j t ) d
j
( t ) f ( )d
Pj
j
1 cos (
j t ) Pj
CE
Thus, the theorem is proved. □
AC
Next, using the Laplace transform of (26), where the relaxation module is defined by the exponential form (23), will find that
a j (t )
hj
j
sin t
0
Pj j ( t ) e a j ( )d d 1 cos( j t ) 0 j
(27)
L If Aj(p) and Bj(p) are the Laplace transforms ( ) of 𝑎𝑗 (𝑡) and f(ξ), respectively, then, applying the
Convolution Theorem the following result is given by
ACCEPTED MANUSCRIPT
sin
t
L j ( t ) f ( )d
0
j p2 j
B j ( p)
where
e
L a j ( )d B j ( p) A j ( p)
0
So, the equation (27) becomes
A j ( p)
j
j
j p j 2
Pj j 1 m2 / A( p) 2 p 1/ j p j p
AN US
Finally,
A j ( p)
m2 A j ( p ) m2 / p 1/ p 1
CR IP T
f ( )
m2
Pj ( p 1)
p( p p j p j m2 j ) 3
2
, Pj
p j
Using the Newton iterative formula, a real root of the polynomial of degree three is found. Then, the partial fraction decomposition and the inverse Laplace transform for Aj(p) will lead to the solution a j (t )
ED
M
of the problem (14) – (15).
4 Computational results
PT
A simply supported beam under the uniform distributed load 𝑝̅ = 4 N/m, which is applied as a creep load at t = 0 (the load is applied suddenly at t = 0 and then held constant), is considered. The length of the beam is L = 4 m, width b = 0.08 m and height h = 0.23 m. These input data lead to the moment of inertia
CE
of the rectangular section: I bℎ3 /12 8 10−5 𝑚4. The polymeric material is taken to have the density of 1565 kg/𝑚3 . For this numerical example, we employ the standard three parameter solid model that
AC
consists in a serial connection of a spring and a Kelvin - Voigt model, [16]. The relaxation modulus is
E (t ) 1.96 10 7 7.84 10 7 e t / 2.24 N / m 2
(29)
with t in seconds and the corresponding creep compliance has the form:
D(t ) 0.51 10 7 0.408 10 7 e t / 11.2
(30)
Quasi - static case In the quasi-static case the inertial forces are ignored. An exact solution for the creep loading applied at
ACCEPTED MANUSCRIPT
t = 0 can be computed using the correspondence principle. The transverse displacements will be of the following form
we ( x, t ) w ( x) D(t )
(31)
elasticity E = 1, [25], [19]:
CR IP T
where D(t) is given in (30) and 𝑤 ̅ is the solution for a similar elastic structure that has the modulus of
px( x 3 2 Lx 2 L3 ) w ( x) 24 I Dynamic analysis
a j (t )
m2 j
t
j
sin
0
AN US
For above input data, the integral equation (27) becomes:
Pj e ( t ) a j ( )d d 1 cos( j t ) j 0 j
(33) where
M
m2 j 2.24 37 j 4 82.88 j 4 , j 104 j 4 , Pj
(32)
0.2 j
ED
Using (28), the Laplace transform of (33) is equal to
PT
A j ( p)
CE
or
A j ( p)
0.2 (2.24 p 1) jp (2.24 p p 2 233 p j 4 21.12 j 4 ) 3
0.2 p 0.45 3 j p ( p 0.45 p 2 104 p j 4 9.43 j 4 )
(34)
AC
This function A j ( p) has the same form as that obtained in [16] by direct application of the Laplace transformation to the integral-differential equation (14): t
a j (t ) 37 j 4 e (t ) / 2.24 a j ( ) d 104 j 4 a j (t ) 0
0.1 [(1) j 1 1] j
(35)
We remark from the form of (35) that a j (t ) are zero for even values of j. Hence, only odd values of j need to be considered in the finding of the Galerkin’s solution.
ACCEPTED MANUSCRIPT
Finally, calculating the inverse Laplace transform for A j ( p) , the solution a j (t ) of the problem (14) – (15) is determined. The transverse deflections w are approximated by (7) for n = 9: 9
w ( x, t ) ai (t ) i ( x)
CR IP T
i 1
Form of the equation (35) leads to the conclusion that as j increases, the importance of the j-th shape will decrease and a total of five terms with odd indices in the above sum will be enough to find an exact solution. The quasi-static and dynamic results for the midpoint transverse deflection are shown in Fig. 1. The
(denoted by t1) in step i being equal to 15i/300.
AN US
numerical solutions are computed for the time interval between zero and 15 seconds, the time value
The dynamic results produced using VIM and Laplace transform with Galerkin's analysis are graphically represented by curve WL1 for uniform distributed loading 𝑝̅ = 4 N/m and WL2 for 𝑝̅ = 8 N/m; we1 is the deflection (31) for 𝑝̅ = 4 N/m and we2 for 𝑝̅ = 8 N/m.
M
It notes that the dynamic results oscillate around of the quasi-static results and the amplitude of vibrations decrease with the increasing of time, due to the presence of the viscoelastic damping. The value of t, denoted
PT
ED
by T, at which the amortization process ends in the section x, will be determined from the condition: WL1(x,T) = we1(x,T)
It is interesting to remark that the numerical results presented in Fig. 1 are in very good agreement with
AC
CE
those obtained by Payette and Reddy, [18].
CR IP T
ACCEPTED MANUSCRIPT
0.01
WL2 2 t1i we2 2 t1i
we1 2 t1i
AN US
WL1 2 t1i 3
510
M
0 0
100
200
300
ED
i
AC
CE
PT
Fig. 1. Quasi-static and dynamic results for the midpoint transverse deflection of the beam
ACCEPTED MANUSCRIPT
WL2 2 t1i we2 2 t1i
0.01
CR IP T
WL2 3.t1i we2 3.t1i
we2 3.5 t1i
3
510
0 0
100
AN US
WL2 3.5 t1i
200
300
i
ED
M
Fig. 2. Quasi-static and dynamic results for the beam sections: x = 2 m, x = 3 m and x = 3.5 m
PT
The Fig. 2 shows the dynamic and quasi-static results for x = 2 m, x = 3 m and x = 3.5 m for 𝑝̅ = 8 N/m and time less than or equal to 15 seconds. It can be seen the increasing of deflection with the increasing of
AC
CE
x, where x is between zero and 2. Then, because symmetry is found that w (2 + x, t) = w (2 - x, t), x [0, 2]
ACCEPTED MANUSCRIPT
6 Conclusions The results obtained in the present study can be summarized as follows: (1) A mixed algorithm based on Galerkin's method, VIM and the techniques of Laplace transformation was created for the dynamic analysis of viscoelastic beams.
CR IP T
(2) Using VIM, the Laplace transform is applied to a simple integral equation (26) and not for the integral - differential equation (14). Numerical example shows that new method leads to the same Laplace transform of the solution a j (t ) , as that obtained by applying the Laplace transform to the differential equation (14), [16].
(3) In this new algorithm is avoided the definition of initial approximation, a j , 0 (t ) , which is a
AN US
difficult
step in using VIM for nonlinear problems.
(4) In the last decade, the viscoelastic structures are analyzed using the finite element models for the spatial domain. Unfortunately, the methods based on the discretization of time interval insert errors
M
that depend on the choice of the time step. Such errors do not accompany the results presented in this paper.
ED
(5) The number of calculations and the errors generated by the iterative process defined by VIM are significantly reduced if (22) is replaced with the integral equation (26), which is solved using the
PT
Laplace transform techniques.
(6) The method proposed in this paper is a powerful tool for the dynamic analysis of viscoelastic beams
CE
that appear frequently in engineering design today. Theoretical results of this study can be extended
AC
easily to the study of more complex structures.
ACCEPTED MANUSCRIPT
References
AC
CE
PT
ED
M
AN US
CR IP T
[1] Argyris,J., Doltsini, I., Silva, V.D., Constitutive modelling and computation of non-linear viscoelastic solids. Part I: Rheological models and numerical integration techniques. Comput. Meth. Appl. Mech. Eng, 88 (1991) 135 -163. [2] Barari, A., Kalijib, A., Ghadimic, M. and Domairry,G., Non-linear vibration of Euler-Bernoulli beams, Latin American Journal of Solids and Structures, 8 (2011)139 – 148. [3] Chen, Q., Chan, Y. W., Integral finite element method for dynamical analysis of elastic-viscoelastic composite structures,” Computers and Structures, 74 (2000) 51–64. [4] Christensen, R.M., Theory of Viscoelasticity, Academic Press, New York, 1982. [5] Chazal,C., Pitti,R., Integral approach for time dependent material using finite element method, Journal of Theoretical and Appl. Mech, 49 (2011) 1029-1048. [6] R. J. Crawford, Plastics Engineering, Elsevier Butterworth – Heinemann, Oxford, 1998. [7] Cristescu, N.D. , Dynamic plasticity, North Holland Pub., Amsterdam, 1967 [8] Flugge, W., Viscoelasticity, Springer, Berlin, Heidelberg, 2nd edition, 1975. [9] He, J. H., Variational iteration method for autonomous ordinary differential systems, Appl. Math. Comput., 114 (2000) 115-123. [10] He, J. H., Comparison of homotopy perturbation method and homotopy analysis method, Appl Math Comp., 156 (2004), 527-539. [11] Janovsky, V., Shaw,V.S., Warby, M.K. , Whiteman, J.R., Numerical methods for treating problems of viscoelastic isotropic solid deformation, J. Comput. Appl. Math., 63 (1995) 91–107. [12] Van Krevelen, D.W., Properties of Polymers, Elsevier, Amsterdam, 1990. [13] Kennedy, T. C., Nonlinear viscoelastic analysis of composite plates and shells. Composite Structures, 41 (1998), 265–272. [14] Lee, U., Oh, H., Dynamic of an axially moving viscoelastic beam subject to axial tension, Int. J. Solids Struct. 42, (2005) 2381-2398. [15] Martin,O., Propagation of elastic-plastic waves in bars, Proceedings of the European Computing Conference, Springer Science, Lecture Notes in Electrical Engineering, 28 (2009), Part II, cap.10, 1011-1024. [16] Martin, O., Quasi-static and dynamic analysis for viscoelastic beams with the constitutive equation in a hereditary integral form, Annals of the University of Bucharest (mathematical series) 5 (LXIII) (2014), 329–341. [17] Park, S.W., Kim, Y.R., Interconversion between relaxation modulus and creep compliance for viscoelastic solids, J. Mater. Civil Eng., 11 (1999), 76 – 82. [18] Payette, G. S., Reddy, J. N., Nonlinear quasi-static finite element formulations for viscoelastic Euler– Bernoulli and Timoshenko beams, International Journal for Numerical Methods in Biomedical Engineering, 26 (2010),1736 – 1755. [19] Pissarenko, G., Yakovlev, A., Matveev, V., Aide-mémoire de résistance des matériaux, Editions de Moscou, 1979. [20] Reddy, J.N. , An Introduction to Continuum Mechanics with Applications Cambridge University Press, 2008. [21] Sarma B.S., Varadan, T.K., Lagrange-type formulation for finite element analysis of nonlinear beam vibrations, J. Sound Vibrations., 86 (1983), 61–70. [22] Shaw, S., Warby, M. K., Whiteman, J. R., A comparison of hereditary integral and internal variable approaches to numerical linear solid viscoelasticity, In Proc. XIII Polish Conf. on Computer Methods in Mechanics, Poznan, (1997). [23] Theodore, R.J., Arakeri, J.H., Ghosal, A., The modeling of axially translating flexible beams, J. Sound Vib., 191 (1996) 363-376.
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
AN US
CR IP T
[24] Zhang Neng-hui, Cheng Chang-jun, Nonlinear mathematical model of viscoelastic thin plates with its Applications, Comput. Methods Appl Mech. Engng., 165 (1998) 307–319. [25] Yang, T. Y., Lianis, G., Large Displacement Analysis of Viscoelastic Beams and Frames by the Finite-Element Method, Journal of Applied Mechanics, 41 (1974), 635-657. [26] Wang, S. Q., He, J. H., Variational iteration method for solving integro-differential equations, Phys. Lett. A, 367 (2007), 188-191.
A modified variational iteration method for the analysis of viscoelastic beams Olga Martin PII: DOI: Reference:
S0307-904X(16)30212-8 10.1016/j.apm.2016.04.011 APM 11132
To appear in:
Applied Mathematical Modelling
Received date: Revised date: Accepted date:
25 April 2015 4 March 2016 13 April 2016
Please cite this article as: Olga Martin , A modified variational iteration method for the analysis of viscoelastic beams, Applied Mathematical Modelling (2016), doi: 10.1016/j.apm.2016.04.011
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Highlights
CR IP T
An algorithm based on VIM and Laplace transformation is created for the dynamic study of viscoelastic beams. It is avoided the definition of initial approximation in the using of VIM, which is a difficult step for nonlinear problems. There are no errors, which usually accompany the methods based on the discretization of time interval.
AC
CE
PT
ED
M
AN US
ACCEPTED MANUSCRIPT
A modified variational iteration method for the analysis of viscoelastic beams
CR IP T
Professor Olga Martin Applied Sciences Faculty Polytechnic University of Bucharest [email protected]
Abstract
AN US
Based on a constitutive law in a hereditary integral form, a mathematical model for dynamic analysis of the isotropic linear viscoelastic beams is presented. To solve the governing equation for these structures subjected to a distributed load is created an accuracy and computational efficiency algorithm that uses the Galerkin’s method and a modified form of the variational iteration method (VIM) for time-domain equations. Numerical results for both quasi-static and dynamic analysis are presented and these will be accompanied by the graphical representations.
Key Words: viscoelastic beam, Euler-Bernoulli beam theory, Galerkin method, variational iteration
M
method, Laplace transform, correspondence principle.
1 Introduction
ED
Viscoelastic materials, such as polymers, rubbers and metals at elevated temperature, are frequently present today in the engineering design due to their ability to dampen out the structural vibrations. There
PT
are remarkable theoretical studies on viscoelastic materials such as those of Flugge [8], Chritensen [4], Reddy [19], Kennedy [13], Krevelen [12], Cristescu [7] and Crawford [6]. Mathematical models, accompanied by the numerical methods were proposed to predict, as accurately, the response of
CE
viscoelastic components which depend on their hereditary loading and its service conditions [1], [3], [10], [14], [16], [21] - [24]. The methods of analysis of viscoelastic structures that occur often in the literature
AC
are the difference method and finite element method [5], [11], [12], [17], [20], [23], [24]. In this paper, the governing equation is determined for a simply supported viscoelastic beam under a
uniform distributed load using Euler-Bernoulli beam theory. This equation is accompanied by a constitutive law presented in a hereditary integral form. In order to obtain the quasi-static exact solution (i.e. the solution ignoring inertia effects), the correspondence principle was used [15], [24]. This principle relates mathematically the solution of a linear, viscoelastic boundary value problem to an analogous problem of an elastic body of the same geometry and under the same initial boundary conditions. It has to
ACCEPTED MANUSCRIPT
be mentioned that not all problems can be solved by this principle, but only those for which the boundary conditions do not vary with the time. Then, applying d’Alembert principle, this structure has been analyzed in the dynamic case by an algorithm that approaches the Galerkin’s method for the spatial domain and a modified variational iteration method, [9], for the temporal domain. The theoretical results
2 Mathematical model
CR IP T
have been applied both for quasi-static and dynamic responses of an isotropic linear viscoelastic beam.
Let us consider a simply supported beam with the span L and transverse section S, subjected to uniformly distributed forces 𝑝̅ . In Euler’s theory, the geometrical equation is given by
2 w( x, t ) z x2
(1)
AN US
z
where x and z represent the axial and transverse coordinates, ρ is the curvature radius, w(x,t) – the deflection and ε- the strain.
M
A constitutive law in the hereditary integral form, [22], [25] is proposed: t
( x, t ) E (0) ( x, t )
(2)
ED
0
dE (t ) ( x, ) d d
where E(t) is the relaxation modulus for the beam material and σ the stress corresponding to the strain ε.
CE
PT
Accordance with the d’Alembert’s principle, the governing equation is of the form
2 M 2w p A x2 t2
(3)
AC
where the bending moment M(x,t) is written as
M ( x, t ) ( x, t ) z dydz
(4)
S
Using (2) the bending moment in a beam of constant cross section is equal to
2 w( x, t ) dE (t ) 2 w( x, ) I d 0 d x2 x2 t
M ( x, t ) E (0) I
(5)
ACCEPTED MANUSCRIPT
where I is the area moment of inertia given by ∬𝑆 𝑧 2 𝑑𝑦𝑑𝑧. Thus, after the substitution of (5) in (3), the following motion equation in terms of the transverse deflection is obtained:
2 w( x, t ) 4 w( x, t ) dE (t ) 4 w( x, ) E ( 0 ) I I d p 0 d t 2 x4 x4 t
3 Mixed method for solving the motion equation
CR IP T
A
(6)
To solve the equation (6) for a simply supported beam subjected to a uniformly distributed loading 𝑝̅ , we express the transverse deflection w by an expansion of the form n
wn ( x, t ) ai (t ) i ( x)
(7)
AN US
i 1
where φi(x) is the i th shape function and ai (t) is the corresponding time-dependent amplitude. For the spatial domain will be used the Galerkin’s method and then, for the time domain will be proposed a new variational iteration method accompanied by the techniques of the Laplace transform.
M
The shape functions are chosen to be linearly independent, orthonormal L
1 i j 0 i j
ED
i ( x) j ( x) dx 0
(8)
PT
and must satisfy all boundary condition for the convergence of Galerkin’s method. Although wn satisfies the boundary conditions, in generally, does not satisfy equation (6). If the expansion
CE
(7) is substituted into (6) will result the residual function
AC
n n Rn ( x, t ) A ai(t ) i ( x) E (0) I ai (t ) i( 4) ( x) i 1 i 1 t
I 0
dE (t ) n ai ( ) i( 4) ( x) d p d i 1
(9)
Now consider the shape functions of the form
i ( x)
2 i x sin , i 1, 2,...., n L L
(10)
ACCEPTED MANUSCRIPT
for which
i i ( x ) i i ( x ) L 4
i( 4) ( x)
(11)
The Galerkin’s method requires that the residual to be orthogonal to each of the shape functions, so
R ( x, t ) n
j
(x) d x d t 0, j 1, 2,..., n
(12)
CR IP T
where [0, L] [0, t ] . This conditions leads to n equations verified by the functions aj(t):
dE (t ) a j ( ) d j I E (0)a j (t ) p j ( x)dx d 0 0 t
L
Aaj (t ) j I or
dE (t ) a j ( ) d j a j (t ) Pj d 0 t
AN US
aj (t ) j where
(13)
(14)
4 L jI j I E (0) 1 j , , , P p j ( x)dx j j j A A A 0 L
j
d kaj
L
t 0
0
ED
dt k
M
The initial conditions are the following
k wn ( x,0) j ( x)dx, k 0, 1, 2,.... t k
(15)
The functions aj(t) are determined independent of one another by (14) – (15). Finally, the value of the
PT
transverse deflection w(x,t) will be found by (7).
CE
3. 1. The variational iteration method (VIM) To obtain the solution of the integral-differential equation (14) will rewrite it using the linear operator L
AC
and the nonlinear operator N:
where
La j (t ) Na j (t ) F j 0
(16)
La j (t ) a j (t ) j a j (t ) dE (t ) a j ( )d , F j Pj d 0 t
N a j (t ) j
According to the variational iteration method (VIM), [9] the correction functional in t direction is defined as:
ACCEPTED MANUSCRIPT
t
a j ,m1 (t ) a j ,m (t ) ( )( La j ,m ( ) Na~ j ,m ( ) F j )d
(17)
0
To obtain the successive approximations a j ,m1 will determine the Lagrange multiplier λ, using the variational techniques. The function 𝑎̃𝑗,𝑚 has a restricted variation in t direction, which means
CR IP T
𝛿𝑎̃𝑗,𝑚 = 0, [26]. Remind that a j ,m (0) 0, u j ,m (0) 0 and F 0. Therefore, it will get, via the integration by parts, that
dE ( ) ~ a j ,m1 (t ) a j ,m (t ) ( ) a j,m ( ) j a j ,m ( )d j a j ,m ( ) F j d d 0 0 t
t 1 - ( ) a j ,m (t ) ( ) t a j ,m (t ) ( ) j ( ) a j ,m ( ) d 0 t 0
So, λ(ξ) verifies the following conditions
( )
t
0, 1 ( )
t
M
and the Lagrange multiplier will be identified as
ED
( )
AN US
1
j
(18)
0, ( ) j ( ) 0
(19)
(20)
sin j t
The iteration formula (17) becomes
dE ( ) a j ,m1 (t ) a j ,m (t ) sin j t a j,m ( ) j a j ,m ( ) d j a j ,m ( ) F j d d j 0 0 t
CE
PT
1
m = 0, 1,…
AC
with the initial conditions:
d k a j ,m (t ) dt k
t 0
0, m, k 0, 1, 2,....
Applying the integration by parts twice for the first term of the integral and (19), will get
(21)
ACCEPTED MANUSCRIPT
a j ,m 1 (t ) a j ,m (t )
1
j
sin j ( t ) a j ,m ( ) t0 cos j ( t ) a j ,m ( ) t0
t
j
t
1
dE ( t ) d j a j ,m ( ) sin j t d sin ( t ) a ( ) d j j , m d 0 j 0 0 j a j ,m ( ) sin j t d
j
t
F
j
Finally, (21) becomes
j
dE ( t ) 1 sin ( t ) a j ,m ( ) d d j d j 0 j 0 t
t
F
j
sin j ( t ) d
0
AN US
a j ,m1 (t )
sin j ( t ) d
0
CR IP T
0
t
(22)
In the case of an exponential form for the relaxation module: t
E (t ) m1 m2 e , m1 , m2 , 0 the iteration formula (22) becomes:
(23)
P e d j 1 cos( j t ) a j ,m1 (t ) sin ( t ) a ( ) d j j , m 0 j j 0
j m2 .
M
(24)
PT
where h j
t
ED
hj
CE
3. 2. Solution of the governing equation for the amplitude 𝒂𝒋 (𝒕) In the case of nonlinear problems, the iterative formula (22) defines a slow process, requiring a longer period of time, even if we use a modern computing technique. On the other hand, the sequence
AC
{𝑎𝑗,𝑚 (𝑡)}𝑚 defined by (21) or (24), converges to the exact solution of the equation (14), if the initial approximation 𝑎𝑗,0 (𝑡) verifies the initial conditions. In this paper, we propose to remove these inconveniences by associating Laplace transform techniques to MIV. Let us now consider that the sequence {𝑎𝑗,𝑚 (𝑡)}𝑚 converges to the solution 𝑎𝑗 (𝑡) of the equation (14). Therefore, the iterative process (22) ceases when 𝑎𝑗,𝑚 (𝑡) = 𝑎𝑗,𝑚+1 (𝑡) with an error less than a value 𝜀̃ very
ACCEPTED MANUSCRIPT
small and the solution 𝑎𝑗 (𝑡) will then be equal to
a j (t ) a j ,m (t ) a j ,m1 (t )
(25)
The formula (22) will become of the form
j
dE ( t ) Pj sin ( t ) a j ( )d d 1 cos( j t ) j d j j 0 0 t
(26)
CR IP T
a j (t )
Theorem. If the function 𝑎𝑗 : [0, 𝑇] → 𝐑, 𝐷 = [0, 𝑇], 𝑎𝑗 ∈ 𝐶 2 (𝐷), verifies the integral equation (26), then it is the solution of the integral-differential equation (14).
0
dE ( t ) a j ( )d . Then, the derivatives of the function a j (t ) will be d
AN US
□ Let be f ( )
t
a j (t ) j cos ( j ( t )) f ( )d Pj 0 t
sin ( j t )
j
a j (t ) j j sin ( j ( t )) f ( )d j f (t ) Pj cos( j t )
M
0
Introducing this last derivative in the left side of the equation (14) is found that
t
ED
j j sin ( j ( t )) f ( )d j 0
0
j t dE ( t ) a j ( )d j sin d j 0
PT
j
0
dE ( t ) a j ( )d Pj cos( j t ) d
j
( t ) f ( )d
Pj
j
1 cos (
j t ) Pj
CE
Thus, the theorem is proved. □
AC
Next, using the Laplace transform of (26), where the relaxation module is defined by the exponential form (23), will find that
a j (t )
hj
j
sin t
0
Pj j ( t ) e a j ( )d d 1 cos( j t ) 0 j
(27)
L If Aj(p) and Bj(p) are the Laplace transforms ( ) of 𝑎𝑗 (𝑡) and f(ξ), respectively, then, applying the
Convolution Theorem the following result is given by
ACCEPTED MANUSCRIPT
sin
t
L j ( t ) f ( )d
0
j p2 j
B j ( p)
where
e
L a j ( )d B j ( p) A j ( p)
0
So, the equation (27) becomes
A j ( p)
j
j
j p j 2
Pj j 1 m2 / A( p) 2 p 1/ j p j p
AN US
Finally,
A j ( p)
m2 A j ( p ) m2 / p 1/ p 1
CR IP T
f ( )
m2
Pj ( p 1)
p( p p j p j m2 j ) 3
2
, Pj
p j
Using the Newton iterative formula, a real root of the polynomial of degree three is found. Then, the partial fraction decomposition and the inverse Laplace transform for Aj(p) will lead to the solution a j (t )
ED
M
of the problem (14) – (15).
4 Computational results
PT
A simply supported beam under the uniform distributed load 𝑝̅ = 4 N/m, which is applied as a creep load at t = 0 (the load is applied suddenly at t = 0 and then held constant), is considered. The length of the beam is L = 4 m, width b = 0.08 m and height h = 0.23 m. These input data lead to the moment of inertia
CE
of the rectangular section: I bℎ3 /12 8 10−5 𝑚4. The polymeric material is taken to have the density of 1565 kg/𝑚3 . For this numerical example, we employ the standard three parameter solid model that
AC
consists in a serial connection of a spring and a Kelvin - Voigt model, [16]. The relaxation modulus is
E (t ) 1.96 10 7 7.84 10 7 e t / 2.24 N / m 2
(29)
with t in seconds and the corresponding creep compliance has the form:
D(t ) 0.51 10 7 0.408 10 7 e t / 11.2
(30)
Quasi - static case In the quasi-static case the inertial forces are ignored. An exact solution for the creep loading applied at
ACCEPTED MANUSCRIPT
t = 0 can be computed using the correspondence principle. The transverse displacements will be of the following form
we ( x, t ) w ( x) D(t )
(31)
elasticity E = 1, [25], [19]:
CR IP T
where D(t) is given in (30) and 𝑤 ̅ is the solution for a similar elastic structure that has the modulus of
px( x 3 2 Lx 2 L3 ) w ( x) 24 I Dynamic analysis
a j (t )
m2 j
t
j
sin
0
AN US
For above input data, the integral equation (27) becomes:
Pj e ( t ) a j ( )d d 1 cos( j t ) j 0 j
(33) where
M
m2 j 2.24 37 j 4 82.88 j 4 , j 104 j 4 , Pj
(32)
0.2 j
ED
Using (28), the Laplace transform of (33) is equal to
PT
A j ( p)
CE
or
A j ( p)
0.2 (2.24 p 1) jp (2.24 p p 2 233 p j 4 21.12 j 4 ) 3
0.2 p 0.45 3 j p ( p 0.45 p 2 104 p j 4 9.43 j 4 )
(34)
AC
This function A j ( p) has the same form as that obtained in [16] by direct application of the Laplace transformation to the integral-differential equation (14): t
a j (t ) 37 j 4 e (t ) / 2.24 a j ( ) d 104 j 4 a j (t ) 0
0.1 [(1) j 1 1] j
(35)
We remark from the form of (35) that a j (t ) are zero for even values of j. Hence, only odd values of j need to be considered in the finding of the Galerkin’s solution.
ACCEPTED MANUSCRIPT
Finally, calculating the inverse Laplace transform for A j ( p) , the solution a j (t ) of the problem (14) – (15) is determined. The transverse deflections w are approximated by (7) for n = 9: 9
w ( x, t ) ai (t ) i ( x)
CR IP T
i 1
Form of the equation (35) leads to the conclusion that as j increases, the importance of the j-th shape will decrease and a total of five terms with odd indices in the above sum will be enough to find an exact solution. The quasi-static and dynamic results for the midpoint transverse deflection are shown in Fig. 1. The
(denoted by t1) in step i being equal to 15i/300.
AN US
numerical solutions are computed for the time interval between zero and 15 seconds, the time value
The dynamic results produced using VIM and Laplace transform with Galerkin's analysis are graphically represented by curve WL1 for uniform distributed loading 𝑝̅ = 4 N/m and WL2 for 𝑝̅ = 8 N/m; we1 is the deflection (31) for 𝑝̅ = 4 N/m and we2 for 𝑝̅ = 8 N/m.
M
It notes that the dynamic results oscillate around of the quasi-static results and the amplitude of vibrations decrease with the increasing of time, due to the presence of the viscoelastic damping. The value of t, denoted
PT
ED
by T, at which the amortization process ends in the section x, will be determined from the condition: WL1(x,T) = we1(x,T)
It is interesting to remark that the numerical results presented in Fig. 1 are in very good agreement with
AC
CE
those obtained by Payette and Reddy, [18].
CR IP T
ACCEPTED MANUSCRIPT
0.01
WL2 2 t1i we2 2 t1i
we1 2 t1i
AN US
WL1 2 t1i 3
510
M
0 0
100
200
300
ED
i
AC
CE
PT
Fig. 1. Quasi-static and dynamic results for the midpoint transverse deflection of the beam
ACCEPTED MANUSCRIPT
WL2 2 t1i we2 2 t1i
0.01
CR IP T
WL2 3.t1i we2 3.t1i
we2 3.5 t1i
3
510
0 0
100
AN US
WL2 3.5 t1i
200
300
i
ED
M
Fig. 2. Quasi-static and dynamic results for the beam sections: x = 2 m, x = 3 m and x = 3.5 m
PT
The Fig. 2 shows the dynamic and quasi-static results for x = 2 m, x = 3 m and x = 3.5 m for 𝑝̅ = 8 N/m and time less than or equal to 15 seconds. It can be seen the increasing of deflection with the increasing of
AC
CE
x, where x is between zero and 2. Then, because symmetry is found that w (2 + x, t) = w (2 - x, t), x [0, 2]
ACCEPTED MANUSCRIPT
6 Conclusions The results obtained in the present study can be summarized as follows: (1) A mixed algorithm based on Galerkin's method, VIM and the techniques of Laplace transformation was created for the dynamic analysis of viscoelastic beams.
CR IP T
(2) Using VIM, the Laplace transform is applied to a simple integral equation (26) and not for the integral - differential equation (14). Numerical example shows that new method leads to the same Laplace transform of the solution a j (t ) , as that obtained by applying the Laplace transform to the differential equation (14), [16].
(3) In this new algorithm is avoided the definition of initial approximation, a j , 0 (t ) , which is a
AN US
difficult
step in using VIM for nonlinear problems.
(4) In the last decade, the viscoelastic structures are analyzed using the finite element models for the spatial domain. Unfortunately, the methods based on the discretization of time interval insert errors
M
that depend on the choice of the time step. Such errors do not accompany the results presented in this paper.
ED
(5) The number of calculations and the errors generated by the iterative process defined by VIM are significantly reduced if (22) is replaced with the integral equation (26), which is solved using the
PT
Laplace transform techniques.
(6) The method proposed in this paper is a powerful tool for the dynamic analysis of viscoelastic beams
CE
that appear frequently in engineering design today. Theoretical results of this study can be extended
AC
easily to the study of more complex structures.
ACCEPTED MANUSCRIPT
References
AC
CE
PT
ED
M
AN US
CR IP T
[1] Argyris,J., Doltsini, I., Silva, V.D., Constitutive modelling and computation of non-linear viscoelastic solids. Part I: Rheological models and numerical integration techniques. Comput. Meth. Appl. Mech. Eng, 88 (1991) 135 -163. [2] Barari, A., Kalijib, A., Ghadimic, M. and Domairry,G., Non-linear vibration of Euler-Bernoulli beams, Latin American Journal of Solids and Structures, 8 (2011)139 – 148. [3] Chen, Q., Chan, Y. W., Integral finite element method for dynamical analysis of elastic-viscoelastic composite structures,” Computers and Structures, 74 (2000) 51–64. [4] Christensen, R.M., Theory of Viscoelasticity, Academic Press, New York, 1982. [5] Chazal,C., Pitti,R., Integral approach for time dependent material using finite element method, Journal of Theoretical and Appl. Mech, 49 (2011) 1029-1048. [6] R. J. Crawford, Plastics Engineering, Elsevier Butterworth – Heinemann, Oxford, 1998. [7] Cristescu, N.D. , Dynamic plasticity, North Holland Pub., Amsterdam, 1967 [8] Flugge, W., Viscoelasticity, Springer, Berlin, Heidelberg, 2nd edition, 1975. [9] He, J. H., Variational iteration method for autonomous ordinary differential systems, Appl. Math. Comput., 114 (2000) 115-123. [10] He, J. H., Comparison of homotopy perturbation method and homotopy analysis method, Appl Math Comp., 156 (2004), 527-539. [11] Janovsky, V., Shaw,V.S., Warby, M.K. , Whiteman, J.R., Numerical methods for treating problems of viscoelastic isotropic solid deformation, J. Comput. Appl. Math., 63 (1995) 91–107. [12] Van Krevelen, D.W., Properties of Polymers, Elsevier, Amsterdam, 1990. [13] Kennedy, T. C., Nonlinear viscoelastic analysis of composite plates and shells. Composite Structures, 41 (1998), 265–272. [14] Lee, U., Oh, H., Dynamic of an axially moving viscoelastic beam subject to axial tension, Int. J. Solids Struct. 42, (2005) 2381-2398. [15] Martin,O., Propagation of elastic-plastic waves in bars, Proceedings of the European Computing Conference, Springer Science, Lecture Notes in Electrical Engineering, 28 (2009), Part II, cap.10, 1011-1024. [16] Martin, O., Quasi-static and dynamic analysis for viscoelastic beams with the constitutive equation in a hereditary integral form, Annals of the University of Bucharest (mathematical series) 5 (LXIII) (2014), 329–341. [17] Park, S.W., Kim, Y.R., Interconversion between relaxation modulus and creep compliance for viscoelastic solids, J. Mater. Civil Eng., 11 (1999), 76 – 82. [18] Payette, G. S., Reddy, J. N., Nonlinear quasi-static finite element formulations for viscoelastic Euler– Bernoulli and Timoshenko beams, International Journal for Numerical Methods in Biomedical Engineering, 26 (2010),1736 – 1755. [19] Pissarenko, G., Yakovlev, A., Matveev, V., Aide-mémoire de résistance des matériaux, Editions de Moscou, 1979. [20] Reddy, J.N. , An Introduction to Continuum Mechanics with Applications Cambridge University Press, 2008. [21] Sarma B.S., Varadan, T.K., Lagrange-type formulation for finite element analysis of nonlinear beam vibrations, J. Sound Vibrations., 86 (1983), 61–70. [22] Shaw, S., Warby, M. K., Whiteman, J. R., A comparison of hereditary integral and internal variable approaches to numerical linear solid viscoelasticity, In Proc. XIII Polish Conf. on Computer Methods in Mechanics, Poznan, (1997). [23] Theodore, R.J., Arakeri, J.H., Ghosal, A., The modeling of axially translating flexible beams, J. Sound Vib., 191 (1996) 363-376.
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
AN US
CR IP T
[24] Zhang Neng-hui, Cheng Chang-jun, Nonlinear mathematical model of viscoelastic thin plates with its Applications, Comput. Methods Appl Mech. Engng., 165 (1998) 307–319. [25] Yang, T. Y., Lianis, G., Large Displacement Analysis of Viscoelastic Beams and Frames by the Finite-Element Method, Journal of Applied Mechanics, 41 (1974), 635-657. [26] Wang, S. Q., He, J. H., Variational iteration method for solving integro-differential equations, Phys. Lett. A, 367 (2007), 188-191.