APPLICATION OF THE HOMOTOPY ANALYSIS METHOD FOR DETERMINING THE FREE VIBRATIONS OF BEAM

Publications

Share / Export Citation / Email / Print / Text size:

Architecture, Civil Engineering, Environment

Silesian University of Technology

Subject: Architecture, Civil Engineering, Engineering, Environmental

GET ALERTS

ISSN: 1899-0142

DESCRIPTION

7
Reader(s)
60
Visit(s)
0
Comment(s)
0
Share(s)

SEARCH WITHIN CONTENT

FIND ARTICLE

Volume / Issue / page

Related articles

VOLUME 11 , ISSUE 1 (March 2018) > List of articles

APPLICATION OF THE HOMOTOPY ANALYSIS METHOD FOR DETERMINING THE FREE VIBRATIONS OF BEAM

Krzysztof GROMYSZ * / Damian SŁOTA *

Citation Information : Architecture, Civil Engineering, Environment. Volume 11, Issue 1, Pages 61-71, DOI: https://doi.org/10.21307/ACEE-2018-006

License : (CC-BY-NC-ND 4.0)

Received Date : 24-October-2017 / Accepted: 19-February-2018 / Published Online: 01-April-2019

ARTICLE

ABSTRACT

Streszczenie

W artykule przedstawiono homotopijną metodę analizy równań opisujących drgania swobodne układu o jednym stopniu swobody opisane równaniami liniowymi oraz nieliniowymi.

Homotopijna metoda analizy daje możliwość poszukiwania rozwiązania szerokiego zakresu problemów opisanych za pomocą równań operatorowych. W artykule przedstawiono przykłady liczbowe w celu potwierdzenia dokładności i szybkiej zbieżności metody. W przypadku liniowym znane było dokładne rozwiązanie, dzięki czemu można było porównać je z rozwiązaniami uzyskanymi za pomocą analizy homotopijnej. Różnice między otrzymanymi roztworami były niewielkie. Zaletą badanej metody jest jednak to, że otrzymano przybliżone rozwiązania w postaci ciągłych funkcji, które można wykorzystać w dalszych analizach, w tym do interpretacji wyników badań doświadczalnych.

1. INTRODUCTION

Vibrations of the building structures can be divided into the forced ones and the free ones. The forced vibrations appear when the variable load acts on the structure. The free vibrations take place when the variable load does not act on the structure during the vibrations and the motion of the structure results from the initial conditions [14].

The vibrations registered while testing the real structures can be presented with the aid of continuous function. Such function is the polynomial of degree n inscribed into the quantized signal obtained during the research. By analyzing these polynomials one may determine the period of vibration and the damping parameters. In case of structures characterized by the constant stiffness and the damping force proportional to the vibration velocity, the period of vibration and the damping parameters are constant. If the previously mentioned conditions are not satisfied, we deal with the nonlinear vibrations [6].

Mathematical modeling of the building structures consists in deriving the differential equations including the linear and nonlinear parameters of the construction. With respect to the form of describing the vibrations of the real structures it is essential to get the solution of a mathematical model also in the form of continuous function. This function can be obtained by applying the homotopy analysis method. The homotopy analysis method has been developed in the 90s of the last century [15-19, 23]. This method gives the possibility to search for the solution of a wide range of problems described by means of the operator equations. In particular, it can be the ordinary and partial differential equations [3, 4, 9, 10, 18, 23, 28], fractional differential or integrodifferential equations [30, 31] (for some other applications of the fractional calculus see for example [11, 21, 24]) and also the integral equations [1, 2, 7, 8, 12, 13, 23]. Some theoretical results concerning the convergence of this method in case of differential equations can be found, among others, in papers [9, 16, 19, 20, 23, 25, 29], whereas in case of integral equations in papers [2, 7, 8, 23].

In the current paper we describe the application of the homotopy analysis method for determining the free vibrations of the beam in the linear and nonlinear case. The numerical examples are presented to confirm the exactness and fast convergence of the introduced method. In the linear case we know the exact solution, so we compare with it the solution obtained with the aid of homotopy analysis method. This case serves then to verify the exactness and convergence of the discussed method. In the nonlinear case the exact solution is unknown, therefore the approximate solution received by using the homotopy analysis method is compared with the approximate solution determined numerically by applying the Mathematica software [5].

2. HOMOTOPY ANALYSIS METHOD

We describe the concept of the homotopy analysis method by the example of its application for solving the ordinary differential equations of the second kind in the form

(1)
d2y(t)dt2+N1(t,y)+N2(t,y)=0,t(a,b),10.21307_ACEE-2018-006-eqn1.jpg

with the initial conditions

(2)
y(a)=ya,y(a)=va,10.21307_ACEE-2018-006-eqn2.jpg

where N1 is the linear operator, N2 is the non-linear operator, ya and va are the given numbers, while y is the unknown function. Let us denote

(3)
N(y)=d2y(t)dt2+N1(t,y)+N2(t,y).10.21307_ACEE-2018-006-eqn3.jpg

In the first step, we specify the homotopy operator ℋ as:

(4)
(Φ,p)(1p)L(Φ(t;p)y0(t))=p h N(Φ(t;p)),10.21307_ACEE-2018-006-eqn4.jpg

where p∈[0,1] is the embedding parameter, h≠ 0 denotes the convergence control parameter [17, 19, 22, 29], y means the initial approximation of the solution of problem 1 and L is the auxiliary linear operator with property L(0).

Considering the equation ℋ(Φp=0), we obtain the so-called zero-order deformation equation

(5)
(1p)L(Φ(t;p)y0(t))=p h N(Φ(t;p)).10.21307_ACEE-2018-006-eqn5.jpg

Substituting p=0 we get, so L(Φ(t;0) – y0(t))=0. However, if we assume p=0, we obtain N(Φ(t; 1)) = 0. Therefore Φ(t; 1) = y(t), where y is the desired solution of equation 1. Thus, the change of parameter p from zero to one corresponds with the change from the trivial problem to the original problem (and thus the solution from y0 to y).

Taking the Maclaurin series of function Φ (with respect to parameter p) we get

(6)
(Φ(t;p)=m=0ym(t)pm,10.21307_ACEE-2018-006-eqn6.jpg

where y0(t) = Φ(t; 0) and

(7)
ym(t)=1m!mΦ(t;p)pm|p=0,m=1,2,3,10.21307_ACEE-2018-006-eqn7.jpg

If the above series is convergent for p = 1, then we obtain the required solution

(8)
y(t)=m=0ym(t).10.21307_ACEE-2018-006-eqn8.jpg

If we are not able to determine the sum of series in 8, then as an approximate solution of considered equation we can accept the partial sum of this series

(9)
y^n(t)=m=0nym(t).10.21307_ACEE-2018-006-eqn9.jpg

In order to determine the form of function ym we differentiate m times, with respect to parameter p, the left and the right hand side of relation 5. Next, we divide the result by m! and we substitute p = 0 which gives the so-called m th order deformation equation (m > 0):

(10)
L(ym(t)χmym1(t))=h R¯m(y¯m1,t),10.21307_ACEE-2018-006-eqn10.jpg

where y¯m1={y0(t),y1(t),,ym1(t)}10.21307_ACEE-2018-006-eqn60.jpg and

(11)
χm={1 m > 10 m 1,10.21307_ACEE-2018-006-eqn11.jpg

and

(12)
R¯m(y¯m1,t)=1(m1)!(m1pm1N(i=0yi(t)pi))|p=0.10.21307_ACEE-2018-006-eqn12.jpg

In view of the above definition of operator N (equation 3) we have for m= 1:

R¯1(y¯0,t)=N(y0(t))=d2y0(t)dt2+N1(t,y0)+N2(t,y0),10.21307_ACEE-2018-006-eqn13.jpg

whereas for m ≥ we get

R¯m(y¯m1,t)=d2ym1(t)dt2+N1(t,ym1)+1(m1)!(m1pm1N2(t,i=0yi(t)pi))|p=0.10.21307_ACEE-2018-006-eqn14.jpg

Thus, the equation 10 for m= 1 takes the form

(13)
L(y1(t))=h(d2y0(t)dt2+N1(t,y0)+N2(t,y0)),10.21307_ACEE-2018-006-eqn15.jpg

whereas for m ≥ 2 it takes the form

L(ym(t))=L(ym1(t))+h(d2ym1(t)dt2+N1(t,ym1)+1(m1)!(m1pm1N2(t,i=0yi(t)pi))|p=0).10.21307_ACEE-2018-006-eqn16.jpg

After selecting operator L the above two equations enable to determine recursively the successive functions ym and, in result, to obtain the solution 8.

As the operator L one can choose the operator of the highest order derivative occurring in the equation. It is the most natural and very often used procedure. This operator can also be selected in other way. In particular, Vajravelu and Van Gorder [27, 29] proposed three methods of selecting the auxiliary linear operator: method of linear partition matching, method of highest order differential matching, method of complete differential matching (see also [23]).

In the considered case we receive the simple form of equations by taking

(15)
L(y)=d2ydt2.10.21307_ACEE-2018-006-eqn17.jpg

Then we get for m=1:

(16)
d2y1(t)dt2=h(d2y0(t)dt2+N1(t,y0)+N2(t,y0)),10.21307_ACEE-2018-006-eqn18.jpg

and for m ≥ 2:

(17)
d2ym(t)dt2=d2ym1(t)dt2+h(d2ym1(t)dt2+N1(t,ym1)+1(m1)!(m1pm1N2(t,i=0yi(t)pi))|p=0).10.21307_ACEE-2018-006-eqn19.jpg

The last component in the above equation is nonlinear. We will be able to determine this nonlinear element by knowing the form of nonlinear operator N2.

As the initial approximation the best is to select a function satisfying the given initial conditions 2, for example in the form

(18)
y0(t)=ya+va(ta).10.21307_ACEE-2018-006-eqn20.jpg

But if one makes another choice of the linear operator (other than in the form of second order derivative), then one can connect the initial approximation with the selected linear operator and the given initial conditions. On can then determine y0 as the solution of linear equation L(y0) = 0 with the same initial conditions as in the original equation.

In order to ensure the uniqueness of solution of equations 16 and 17 we need to define the corresponding initial conditions. If we put to each equation the homogeneous initial conditions

(19)
ym(a)=0,ym(a)=010.21307_ACEE-2018-006-eqn21.jpg

for m ≥ 1, then the approximate solution ŷn,defined by relation 9, will satisfy the initial conditions 2.

In this way the solution of considered problem is reduced to the solution of a sequence of differential equations 16 and 17 with the homogeneous initial conditions 19. In these equations the sought functions appear only once within the symbol of second derivative. Thus the solution of these equations consist in the double integration of their right hand sides.

The proper selection of the convergence control parameter h affects the area of convergence of the created series and the rate of this convergence as well [19, 20, 26]. For determining this coefficient we can use the so-called “optimization method” [19], in which we define the following squared residual of governing equation

(20)
En(h)=ab(N[y^n(t)])2 dt.10.21307_ACEE-2018-006-eqn22.jpg

And next, the optimum value of the convergence control parameter is obtained by determining the minimum of this squared residual.

3. THE LINEAR HOMOGENEOUS CASE

Vibrations of the simply supported reinforced concrete beam of span l= 6 m and cross-section (b/h) 0.30 m/0.60 m can be described with the aid of a linear differential equation. Uniformly distributed mass m′ of the beam, under the assumption that the density of the beam material is ρ = 2500 kg/m3, is equal to

(21)
m=b h ρ=450 kg/m.10.21307_ACEE-2018-006-eqn23.jpg

Such beam can be considered as the set of segments of length dx and the kinetic energy of such elementary segment, treated as the element of mass mdx (Figure 1a), is equal to

(22)
mdx y˙2(x)2,10.21307_ACEE-2018-006-eqn24.jpg

Figure 1.

Reduction of the beam with the continuous distribution of mass into the system with one degree of freedom a) system with the continuous distribution of mass m′, b) substitute mass m lumped in the point in which the arrow of vibrations of the system with continuous distribution of mass occurs

10.21307_ACEE-2018-006-f001.jpg

where (x) is the velocity of point of coordinate x (in the beam). In case of the simply supported beam wrapped by the mass m′ the deflection line y(x) is symmetric with respect to the center of gravity and for 0x12l10.21307_ACEE-2018-006-eqn61.jpg is expressed by relation

(23)
y(x)=165ym(xl2x3l3+x4l4),10.21307_ACEE-2018-006-eqn25.jpg

where ym is the largest deflection of the beam (Figure 1a).

Velocity of every point of coordinate x is equal to

(24)
y˙(x)=165y˙m(xl2x3l3+x4l4),10.21307_ACEE-2018-006-eqn26.jpg

where y˙m=dymdt10.21307_ACEE-2018-006-eqn62.jpg and the kinetic energy of the whole beam is equal to

(25)
Ek=20l/2my˙2(x)2 dx.10.21307_ACEE-2018-006-eqn27.jpg

By substituting 24 into 25 and by calculating the integrals we get

(26)
Ek=39687875mly˙m22=0.504 mly˙m22.10.21307_ACEE-2018-006-eqn28.jpg

The kinetic energy described by equation 26 corresponds to the kinetic energy of the system with one degree of freedom and the substitute mass equal to

(27)
m=0.504 ml=1360.8 kg10.21307_ACEE-2018-006-eqn29.jpg

focused in the deflection, that is in the point in which the deflection of a beam is the largest (Figure 1b). Stiffness of the substitute system is equal to

(28)
k=48E Il3,10.21307_ACEE-2018-006-eqn30.jpg

where E= 27 GPa is the modulus of elasticity of the beam material and I is the moment of inertia of the cross-section which is equal to b h312=5.4·103m410.21307_ACEE-2018-006-eqn63.jpg Thus we have k= 32.4 MN/m.

Equation describing the first free vibration of the reinforced concrete beam of span l= 6 m and cross-section (b/h) 0.30 m/0.60 m is then of the form

(29)
m y¨(t)+k y(t)=0,10.21307_ACEE-2018-006-eqn31.jpg

where m= 1360.8 kg denotes the substitute mass of the beam and k= 32.4 MN/m means the stiffness of the beam in the middle of its span.

Equation 29, after dividing by m, takes the form

(30)
y¨(t)+ω02 y(t)=0,10.21307_ACEE-2018-006-eqn32.jpg

where ω0=km=140.8573 rad/s10.21307_ACEE-2018-006-eqn64.jpg is the natural frequency of non-damping vibrations. The above equation must be completed with the initial conditions

(31)
y(0)=ya,y˙(0)=va.10.21307_ACEE-2018-006-eqn33.jpg

Example 1.

We show now how the described above homotopy analysis method may be used for solving the equation 30 with conditions 31. For the value ω0 = 140.9 rad/s and the initial conditions ya=0.01 m and va=0 m/s, the exact solution is given by function

ye(t)=0.01cos(140.9 t).10.21307_ACEE-2018-006-eqn34.jpg

If we take as the initial approximation the function satisfying the initial conditions

y0(t)=ya+va t=0.01,10.21307_ACEE-2018-006-eqn35.jpg

then in the first steps of the method we get

y1(t)=99.2641 h t2,y2(t)=99.2641 (h t2+h2 t2)+164223h2 t4,10.21307_ACEE-2018-006-eqn36.jpg

Figure 3 presents the squared residual for n= 9. Numerically determined optimum value of the convergence control parameter is equal to −0.950692782. In Table 1 there are collected the errors in reconstructing the exact solution for the successive approximate solutions ŷn, n∈ {1,2, …, 15}. The presented errors are calculated in the supremum norm

Λn=yey^n= sup |ye(t)y^n(t)|,t(a,b)10.21307_ACEE-2018-006-eqn37.jpg

Figure 2.

Oscillator of parameters k and m as the substitute system for the reinforced concrete beam

10.21307_ACEE-2018-006-f002.jpg
Figure 3.

The squared residual E9

10.21307_ACEE-2018-006-f003.jpg
Table 1.

Values of errors in the reconstruction of the exact solution (Δn = ||yeŷn||)

10.21307_ACEE-2018-006-tbl1.jpg

where ye means the exact solution and ŷn denotes the approximate solution. At first, for n= 2 the errors increase but afterwards they quickly decrease. So, for n= 5 the errors are about 10–2, for n= 10 they are about 10–8 and for n= 15 the errors are at the level of 10–14. Next the error decreases a little bit slower, so for n= 20 the error is equal to 3.747 · 10–16. In Figure 4 the absolute error |ye(t) – ŷn(t)| for n= 6 and n= 10 is plotted. As revealed by the above results, starting from some moment the errors quickly decrease when the number of components in the sum 9 increases.

Figure 4.

Distribution of error (|ye (x) – ŷn (x)|) of the exact solution approximation for n = 6 (a) and n = 10 (b)

10.21307_ACEE-2018-006-f004.jpg

Example 2.

In the second example we execute the calculations for the changed initial conditions by taking ya = 0.01 m and va=0.2 m/s. Value of parameter ω0 is not changed (ω0=140.9 rad/s). In this case we also know the exact solution given by the function

ye(t)=0.01cos(140.9 t)+0.00141945sin(140.9 t).10.21307_ACEE-2018-006-eqn38.jpg

As the initial approximation we take the function fulfilling the initial conditions

y0(t)=ya+va t=0.01+0.2 t.10.21307_ACEE-2018-006-eqn39.jpg

Then we obtain successively

y1(t)=99.2641 h t2+661.76 h t3,y2(t)=99.2641 (h t2+h2 t2)+661.76 (h t3+h2 t2)+164223 h2 t4+656890 h2 t5,10.21307_ACEE-2018-006-eqn40.jpg

Figure 5 displays the squared residual for n = 9.

Figure 5.

The squared residual E9

10.21307_ACEE-2018-006-f005.jpg

Numerically computed optimum value of the convergence control parameter is equal to –0.9580711. In Table 2 there are compiled the errors in reconstructing the exact solution for the successive approximate solutions ŷn, n∈ {1,2, …, 15}. Similarly as in the previous example, at first the errors increase but afterwards they quickly decrease. In Figure 6 there is plotted the absolute error of reconstructing the exact solution for n= 10 and n= 15.

Table 2.

Values of errors in the reconstruction of the exact solution (Δn = ||yeŷn||)

10.21307_ACEE-2018-006-tbl2.jpg
Figure 6.

Distribution of error (|ye (x) – ŷn (x)|) of the exact solution approximation for n=10 (a) and n=15 (b)

10.21307_ACEE-2018-006-f006.jpg

4. THE LINEAR NONHOMOGENEOUS CASE

The initial condition

(32)
y(0)=0,y˙(0)=va10.21307_ACEE-2018-006-eqn41.jpg

of the free vibrations can be achieved by dropping on the beam from height h the body of mass m1 (the beater – Figure 7a). When the impact is of plastic character then the structure and the beater would perform the free vibrations (Figure 7b). Assuming m1 = 100 kg and h= 0.5 m velocity of the beater is equal to

(33)
v1=2 g h=3.13 m/s,10.21307_ACEE-2018-006-eqn42.jpg

Figure 7.

Impact of plastic character implying the initial conditions y(0)=0 and ẏ (0)=va (a, b – described in text)

10.21307_ACEE-2018-006-f007.jpg

where g is the acceleration of gravity. Basing on the momentum conservation principle the velocity va of the system after the strike is equal to

(34)
va=m1v1m1+m=0.181 m/s.10.21307_ACEE-2018-006-eqn43.jpg

After the impact of the plastic character the motion of the oscillator is executed according to the equation

(35)
(m1+m)y¨(t)+k y(t)=m1 g,10.21307_ACEE-2018-006-eqn44.jpg

the solution of which, after taking into account the initial conditions

(36)
y(0)=0,y˙(0)=va,10.21307_ACEE-2018-006-eqn45.jpg

is given by the function

(37)
y(t)=vaω01sin(ω01 t)+m1 gk(1cos(ω01 t)),10.21307_ACEE-2018-006-eqn46.jpg

where

ω01=km1+m=136.73 rad/s.10.21307_ACEE-2018-006-eqn47.jpg

Example 3.

Now we deal with the solution of equation 35 with conditions 36, for the following parameters m1 = 100 kg, m = 1633 kg, k= 32.4 MN/m, g=9.81 m/s 2, and the initial conditions ya=0 m and va=0.181 m/s. This time by using our method we get in turn

y1(t)=0.283035 h t2+563.159 h t3,y2(t)=0.283035 (h t2h2 t2)+563.159 (h t3+h2 t3)440.967 h2 t4+526438 h2 t5.10.21307_ACEE-2018-006-eqn48.jpg

Numerically determined optimum value of the convergence control parameter is equal to –0.77670315. Table 3 compiles the errors in reconstructing the exact solution for the successive approximate solutions ŷn, n ∈ {1,2, …, 15}. In Figure 8 the absolute error of reconstructing the exact solution for n= 10 and n= 15 is presented. Similarly as before the errors quickly decrease.

Table 3.

Values of errors in the reconstruction of the exact solution (Δn = ||yeŷn||)

10.21307_ACEE-2018-006-tbl3.jpg
Figure 8.

Distribution of error (|ye (x) – ŷn (x)|) of the exact solution approximation for n=10 (a) and n=15 (b)

10.21307_ACEE-2018-006-f008.jpg

5. NONLINEAR CASE

Nonlinearity in the mechanical dynamic system can occur in result of nonlinearity of the elastic force S, the value of which does not depend linearly on the displacement y. Characteristic of this force can be presented in the form of the well known in literature Duffing’s proposal

(38)
s(y)=k y±k^ y3.10.21307_ACEE-2018-006-eqn49.jpg

The plus sign preceding k^10.21307_ACEE-2018-006-eqn65.jpg means that the stiffness increases together with the displacement y, whereas the minus sign means that it decreases. Nonlinearity is the greater, the greater is the displacement. Equation describing such nonlinear free vibrations with no damping takes the form

(39)
m y¨(t)+k y(t)±k^ y3(t)=0,10.21307_ACEE-2018-006-eqn50.jpg

with the initial conditions

(40)
y(0)=0,y˙(0)=va.10.21307_ACEE-2018-006-eqn51.jpg

Referring to the previous example when we assumed the mass m= 1633 kg and stiffness k= 32.4 MN/m, the parameter k^10.21307_ACEE-2018-006-eqn66.jpg can take the values about 109 with the amplitude of variable y up to 0.002 m. Certainly the value of k^10.21307_ACEE-2018-006-eqn67.jpg can be increased and then the nonlinearity effect would be greater, however, for the bigger amplitudes of y, the problem loses its physical sense.

The free vibrations with the damping include additionally one more term by the first derivative denoted usually as c:

(41)
m y¨(t)+c y˙+k y(t)±k^ y3(t)=0.10.21307_ACEE-2018-006-eqn52.jpg

In our case the physically real value of c (determined from the engineering relations, ζ= 0.05 – damping ratio) is equal to

(42)
c=2ζk m=23103 kg/s.10.21307_ACEE-2018-006-eqn53.jpg

Example 4.

Proceeding with the computational example we discuss the solution of equation 39 with initial conditions 40 for m = 1633 kg, k= 32.4 MN/m, k^=10910.21307_ACEE-2018-006-eqn68.jpg and va=0.2 m/s, under the assumption that together with the increase of variable y the stiffness decreases (therefore we use the minus sign in the proper equation). This time we do not know the exact solution of considered equation. That is way the approximate solution obtained by applying the homotopy analysis method will be compared with the approximate solution received numerically with the aid of Mathematica software. Assuming as the initial approximation y0 the function satisfying the initial conditions

y0(t)=va t=0.2 t,10.21307_ACEE-2018-006-eqn54.jpg

we get, by using the homotopy analysis method, the function

y1(t)=661.359 h t3244.948 h t5,y2(t)=661.359 (h t3+h2 t3)244.948 h t5+655850 h2 t51.27285106 h2 t7+249998 h2 t9,10.21307_ACEE-2018-006-eqn55.jpg

The optimum value of the convergence control parameter is equal to –0.7227439. The plot of the squared residual for n= 20 is presented in Figure 9. Figure 10 shows the plot of residual

Rsn(t)=|1633y^n(t)+32.4106y^n(t)109(y^n(t))3|10.21307_ACEE-2018-006-eqn56.jpg

Figure 9.

The squared residual E20

10.21307_ACEE-2018-006-f009.jpg
Figure 10.

Plot of the residual Rsn for n=10 (a) and n=20 (b)

10.21307_ACEE-2018-006-f010.jpg

for n= 10 and n= 20, that is the plot of error of satisfying the equation 39 by the approximate solution ŷn. The above error decreases very quickly when the number of elements in the approximate solution increases (with the increase of number n). So, for n= 5 the maximum of this error takes the value 67.72, for n= 10 this maximum does not exceed the value 0.45, for n= 15 the value 0.081, and finally for n= 20 it does not exceed the value 0.0003. The initial conditions are fulfilled exactly by every approximate solution.

Figure 11 presents the comparison of approximate solutions obtained by applying the homotopy analysis method (solid line) and the numerical methods available in the Mathematica software (dots). The differences are slight even for the small number n of elements (see Table 4).

Figure 11.

Comparison of the approximate solutions (solid line – solution ŷ20 obtained by the homotopy analysis method, dots – approximate solution obtained numerically in Mathematica software)

10.21307_ACEE-2018-006-f011.jpg
Table 4.

Maximal absolute differences (Δn) between the approximate solutions obtained by the homotopy analysis method and the numerical methods available in Mathematica software

10.21307_ACEE-2018-006-tbl4.jpg

Example 5.

In the next example we solve the equation 41 with initial conditions 40 for m = 1633 kg, k = 32.4 MN/m, k^=10910.21307_ACEE-2018-006-eqn69.jpg, c= 23·103 kg/s and va = 0.2 m/s, under the assumption that when the coordinate y increases then the stiffness decreases (so we include the minus sign in the proper equation). This time again we do not know the exact solution of investigated equation. Taking, similarly as before, the function satisfying the initial conditions as the initial approximation y0, that is

y0(t)=va t=0.2 t,10.21307_ACEE-2018-006-eqn57.jpg

we get successively the functions

y1(t)=1.40857 h t2+661.359 h t3244.948 h t5,y2(t)=1.40857 (h t2+h2 t2)+661.359 h t3+667.973 h2 t3+4657.87 h2 t4244.948 h t5+655850 h2 t54025.32 h2 t61.27285106 h2 t7+249998 h2 t9,10.21307_ACEE-2018-006-eqn58.jpg

Optimum value of the convergence control parameter is equal to –0.686325. The plot of the squared residual for n= 20 is shown in Figure 12. In Figure 13 there is displayed the plot of residual

R^sn(t)=|1633y^n(t)+23103y^n(t)+32.4106y^n(t)109(y^n(t))3|10.21307_ACEE-2018-006-eqn59.jpg

Figure 12.

The squared residual E20

10.21307_ACEE-2018-006-f012.jpg
Figure 13.

Plot of the residual R^sn10.21307_ACEE-2018-006-eqn70.jpg for n=10 (a) and n=20 (b)

10.21307_ACEE-2018-006-f013.jpg

for n= 10 and n= 20, it means the plot of error of satisfying the equation 41 by the approximate solution ŷn. Similarly as in the previous example, the above error decreases very quickly when the number of elements in the approximate solution increases. For n= 5 the maximum of this error takes the value 66.37 for n= 10 this maximum does not exceed the value 0.78, for n= 15 the value 0.12 and, in turn, for n= 20 it does not exceed the value 0.00052. The initial conditions are satisfied exactly by every approximate solution.

Figure 14 shows the comparison of approximate solutions obtained by applying the homotopy analysis method (solid line) and the numerical methods available in the Mathematica software (dots). The differences are slight even for small number n of components. For n= 5 the absolute value of the difference does not exceed the value 2.56 · 10–4, for n= 10 it does not exceed 3.12 · 10–6, for n= 15 the value 3.61 · 10–8, whereas for n= 20 it is not higher than 4.2 · 10–9.

Figure 14.

Comparison of the approximate solutions (solid line – solution ŷ20 obtained by the homotopy analysis method, dots – approximate solution obtained by the numerical method available in Mathematica software)

10.21307_ACEE-2018-006-f014.jpg

6. CONCLUSIONS

In the paper we have presented the application of the homotopy analysis method for determining the free vibrations of the reinforced concrete beam in the linear and nonlinear case. In the investigated method the series is created, the successive terms of which are derived by calculating the proper integrals obtained from the previous terms. The formed series is in general quickly convergent. For the optimal selection of the convergence control parameter computation of just several first terms of the series ensures a very good approximation of the sought solution. Performed calculations show that the proposed method is effective in solving the problems under consideration.

The presented computational examples confirm the precision and the fast convergence of investigated method. In the linear case we knew the exact solutions, so we could compare with them the solutions obtained with the aid of homotopy analysis method. This case served then to verify the exactness and convergence of the discussed method. In the nonlinear case, when the exact solution was unknown, we compared the approximate solutions received by using the homotopy analysis method with the approximate solutions determined numerically by applying the Mathematica software. Differences between the obtained solutions were slight. However, the advantage of the examined method is that we receive here the approximate solution in the form of continuous function which can be used then in a further analysis or to perform various simulations. Whereas when we apply the numerical methods, we get the discrete set of values which may be potentially used later for approximating the appropriate functions.

The executed calculations indicate that it is possible to obtain the mathematical solution of the model describing the free vibrations of the reinforced concrete beam in the form of polynomial. The obtained results can be compared with the results of experimental research also approximated by polynomials. Such comparison will give the possibility to verify the taken model of structure.

References


  1. S. Abbasbandy and E. Shivanian, (2011). A new analytical technique to solve Fredholm’s integral equations, Numer. Algor. 56, 27–43.
    [CROSSREF]
  2. R. Brociek, E. Hetmaniok, J. Matlak, and D. Słota, (2016). Application of the homotopy analysis method for solving the systems of linear and nonlinear integral equations, Math. Model. Anal. 21, 350–370.
    [CROSSREF]
  3. D.S. Chauhan, R. Agrawal, and P. Rastogi, (2012). Magnetohydrodynamic slip flow and heat transfer in a porous medium over a stretching cylinder: homotopy analysis method, Numer. Heat Transfer A 62, 136–157.
  4. T. Fan and X. You, (2013). Optimal homotopy analysis method for nonlinear differential equations in the boundary leyer, Numer. Algor. 62, 337–354.
    [CROSSREF]
  5. H. Gliński, R. Grzymkowski, A. Kapusta, and D. Słota, (2012) Mathematica 8, WPKJS, Gliwice (in Polish).
  6. K. Gromysz, (2013). Examination of the Reinforced Concrete Plates Loaded Temporarily, Periodicaly and Kinematically, Monograph no 452, Silesian University of Technology Press, Gliwice (in Polish).
  7. E. Hetmaniok, I. Nowak, D. Słota, and R. Wituła, (2015). Convergence and error estimation of homotopy analysis method for some type of nonlinear and linear integral equations, J. Numer. Math. 23, 331–344.
    [CROSSREF]
  8. E. Hetmaniok, D. Słota, T. Trawiński, and R. Wituła, (2014). Usage of the homotopy analysis method for solving the nonlinear and linear integral equations of the second kind, Numer. Algor. 67, 163–185.
    [CROSSREF]
  9. E. Hetmaniok, D. Słota, R. Wituła, and A. Zielonka, (2015). An analytical method for solving the two-phase inverse Stefan problem, Bull. Pol. Ac.: Tech 63(3), 583–590.
  10. E. Hetmaniok, D. Słota, R. Wituła, and A. Zielonka, (2015). Solution of the one-phase inverse Stefan problem by using the homotopy analysis method, Appl. Math. Modelling 39, 6793–6805.
    [CROSSREF]
  11. T. Kaczorek, (2016). A new approach to the realization problem for fractional discrete-time linear systems, Bull. Pol. Ac.: Tech 64(1), 9–14.
  12. Y. Khan, M. Fardi, (2015). A new efficient multi-parametric homotopy approach for two-dimensional Fredholm integral equations of the second kind, Hacettepe J. Math. Stat. 44, 93–99.
  13. Y. Khan, K. Sayevand, M. Fardi, M. Ghasemi, (2014). A novel computing multi-parametric homotopy approach for system of linear and nonlinear Fredholm integral equations, Appl. Math. Comput. 249, 229–236.
  14. R. Lewandowski, (2006). Dynamics of the Building Structures, Poznań University of Technology Press, Poznań 2006 (in Polish).
  15. S. Liao, (1998). Homotopy analysis method: a new analytic method for nonlinear problems, Appl. Math. Mech. – Engl. Ed. 19, 957–962.
    [CROSSREF]
  16. S. Liao, (2003). Beyond Perturbation: Introduction to the Homotopy Analysis Method, Chapman and Hall–CRC Press, Boca Raton.
    [CROSSREF]
  17. S. Liao, (2009). Notes on the homotopy analysis method: some definitions and theorems, Commun. Nonlinear Sci. Numer. Simulat. 14, 983–997.
    [CROSSREF]
  18. S. Liao, (2010). An optimal homotopy-analysis approach for strongly nonlinear differential equations, Commun. Nonlinear Sci. Numer. Simulat. 15, 2003–2015.
    [CROSSREF]
  19. S. Liao, (2012). Homotopy Analysis Method in Nonlinear Differential Equations, Springer/Higher Education Press, Berlin/Beijing.
    [CROSSREF]
  20. Z.M. Odibat, (2010). A study on the convergence of homotopy analysis method, Appl. Math. Comput. 217, 782–789.
  21. P. Ostalczyk, (2015). On simplified forms of the fractional-order backward difference and related fractional-order linear discrete-time system description, Bull. Pol. Ac.: Tech 63(2), 423–433.
  22. M. Russo and R.A. Van Gorder, (2013). Control of error in the homotopy analysis of nonlinear Klein-Gordon initial value problems, Appl. Math. Comput. 219, 6494–6509.
  23. D. Słota, (2011). Homotopy Analysis Method and the Examples of Its Applications, Wyd. Pol. Śl., Gliwice (in Polish).
  24. W. Sumelka, (2016). Fractional calculus for continuum mechanics – anisotropic non-locality, Bull. Pol. Ac.: Tech 64 (2), 361–372.
  25. M. Turkyilmazoglu, (2010). A note on the homotopy analysis method, Appl. Math. Lett. 23, 1226–1230.
    [CROSSREF]
  26. M. Turkyilmazoglu, (2010). Convergence of the homotopy analysis method, arXiv, 1006.4460v1.
  27. K. Vajravelu, R.A. Van Gorder, (2012). Nonlinear Flow Phenomena and Homotopy Analysis. Fluid Flow and Heat Transfer, Springer/Higher Education Press, Berlin/Beijing.
    [CROSSREF]
  28. R.A. Van Gorder, (2012). Control of error in the homotopy analysis of semi-linear elliptic boundary value problems, Numer. Algor. 61, 613–629.
    [CROSSREF]
  29. R.A. Van Gorder and K. Vajravelu, (2009). On the selection of auxiliary functions, operators, and convergence control parameters in the application of the homotopy analysis method to nonlinear differential equations: a general approach, Commun. Nonlinear Sci. Numer. Simulat. 14, 4078–4089.
    [CROSSREF]
  30. X. Zhang, B. Tang, and Y. He, (2011). Homotopy analysis method for higher-order fractional integrodifferential equations, Comput. Math. Appl. 62, 3194–3203.
    [CROSSREF]
  31. M. Zurigat, S. Momani, Z. Odibat, and A. Alawneh, (2010). The homotopy analysis method for handling systems of fractional differential equations, Appl. Math. Modelling 34, 24–35.
    [CROSSREF]
XML PDF Share

FIGURES & TABLES

Figure 1.

Reduction of the beam with the continuous distribution of mass into the system with one degree of freedom a) system with the continuous distribution of mass m′, b) substitute mass m lumped in the point in which the arrow of vibrations of the system with continuous distribution of mass occurs

Full Size   |   Slide (.pptx)

Figure 2.

Oscillator of parameters k and m as the substitute system for the reinforced concrete beam

Full Size   |   Slide (.pptx)

Figure 3.

The squared residual E9

Full Size   |   Slide (.pptx)

Figure 4.

Distribution of error (|ye (x) – ŷn (x)|) of the exact solution approximation for n = 6 (a) and n = 10 (b)

Full Size   |   Slide (.pptx)

Figure 5.

The squared residual E9

Full Size   |   Slide (.pptx)

Figure 6.

Distribution of error (|ye (x) – ŷn (x)|) of the exact solution approximation for n=10 (a) and n=15 (b)

Full Size   |   Slide (.pptx)

Figure 7.

Impact of plastic character implying the initial conditions y(0)=0 and ẏ (0)=va (a, b – described in text)

Full Size   |   Slide (.pptx)

Figure 8.

Distribution of error (|ye (x) – ŷn (x)|) of the exact solution approximation for n=10 (a) and n=15 (b)

Full Size   |   Slide (.pptx)

Figure 9.

The squared residual E20

Full Size   |   Slide (.pptx)

Figure 10.

Plot of the residual Rsn for n=10 (a) and n=20 (b)

Full Size   |   Slide (.pptx)

Figure 11.

Comparison of the approximate solutions (solid line – solution ŷ20 obtained by the homotopy analysis method, dots – approximate solution obtained numerically in Mathematica software)

Full Size   |   Slide (.pptx)

Figure 12.

The squared residual E20

Full Size   |   Slide (.pptx)

Figure 13.

Plot of the residual R^sn for n=10 (a) and n=20 (b)

Full Size   |   Slide (.pptx)

Figure 14.

Comparison of the approximate solutions (solid line – solution ŷ20 obtained by the homotopy analysis method, dots – approximate solution obtained by the numerical method available in Mathematica software)

Full Size   |   Slide (.pptx)

REFERENCES

  1. S. Abbasbandy and E. Shivanian, (2011). A new analytical technique to solve Fredholm’s integral equations, Numer. Algor. 56, 27–43.
    [CROSSREF]
  2. R. Brociek, E. Hetmaniok, J. Matlak, and D. Słota, (2016). Application of the homotopy analysis method for solving the systems of linear and nonlinear integral equations, Math. Model. Anal. 21, 350–370.
    [CROSSREF]
  3. D.S. Chauhan, R. Agrawal, and P. Rastogi, (2012). Magnetohydrodynamic slip flow and heat transfer in a porous medium over a stretching cylinder: homotopy analysis method, Numer. Heat Transfer A 62, 136–157.
  4. T. Fan and X. You, (2013). Optimal homotopy analysis method for nonlinear differential equations in the boundary leyer, Numer. Algor. 62, 337–354.
    [CROSSREF]
  5. H. Gliński, R. Grzymkowski, A. Kapusta, and D. Słota, (2012) Mathematica 8, WPKJS, Gliwice (in Polish).
  6. K. Gromysz, (2013). Examination of the Reinforced Concrete Plates Loaded Temporarily, Periodicaly and Kinematically, Monograph no 452, Silesian University of Technology Press, Gliwice (in Polish).
  7. E. Hetmaniok, I. Nowak, D. Słota, and R. Wituła, (2015). Convergence and error estimation of homotopy analysis method for some type of nonlinear and linear integral equations, J. Numer. Math. 23, 331–344.
    [CROSSREF]
  8. E. Hetmaniok, D. Słota, T. Trawiński, and R. Wituła, (2014). Usage of the homotopy analysis method for solving the nonlinear and linear integral equations of the second kind, Numer. Algor. 67, 163–185.
    [CROSSREF]
  9. E. Hetmaniok, D. Słota, R. Wituła, and A. Zielonka, (2015). An analytical method for solving the two-phase inverse Stefan problem, Bull. Pol. Ac.: Tech 63(3), 583–590.
  10. E. Hetmaniok, D. Słota, R. Wituła, and A. Zielonka, (2015). Solution of the one-phase inverse Stefan problem by using the homotopy analysis method, Appl. Math. Modelling 39, 6793–6805.
    [CROSSREF]
  11. T. Kaczorek, (2016). A new approach to the realization problem for fractional discrete-time linear systems, Bull. Pol. Ac.: Tech 64(1), 9–14.
  12. Y. Khan, M. Fardi, (2015). A new efficient multi-parametric homotopy approach for two-dimensional Fredholm integral equations of the second kind, Hacettepe J. Math. Stat. 44, 93–99.
  13. Y. Khan, K. Sayevand, M. Fardi, M. Ghasemi, (2014). A novel computing multi-parametric homotopy approach for system of linear and nonlinear Fredholm integral equations, Appl. Math. Comput. 249, 229–236.
  14. R. Lewandowski, (2006). Dynamics of the Building Structures, Poznań University of Technology Press, Poznań 2006 (in Polish).
  15. S. Liao, (1998). Homotopy analysis method: a new analytic method for nonlinear problems, Appl. Math. Mech. – Engl. Ed. 19, 957–962.
    [CROSSREF]
  16. S. Liao, (2003). Beyond Perturbation: Introduction to the Homotopy Analysis Method, Chapman and Hall–CRC Press, Boca Raton.
    [CROSSREF]
  17. S. Liao, (2009). Notes on the homotopy analysis method: some definitions and theorems, Commun. Nonlinear Sci. Numer. Simulat. 14, 983–997.
    [CROSSREF]
  18. S. Liao, (2010). An optimal homotopy-analysis approach for strongly nonlinear differential equations, Commun. Nonlinear Sci. Numer. Simulat. 15, 2003–2015.
    [CROSSREF]
  19. S. Liao, (2012). Homotopy Analysis Method in Nonlinear Differential Equations, Springer/Higher Education Press, Berlin/Beijing.
    [CROSSREF]
  20. Z.M. Odibat, (2010). A study on the convergence of homotopy analysis method, Appl. Math. Comput. 217, 782–789.
  21. P. Ostalczyk, (2015). On simplified forms of the fractional-order backward difference and related fractional-order linear discrete-time system description, Bull. Pol. Ac.: Tech 63(2), 423–433.
  22. M. Russo and R.A. Van Gorder, (2013). Control of error in the homotopy analysis of nonlinear Klein-Gordon initial value problems, Appl. Math. Comput. 219, 6494–6509.
  23. D. Słota, (2011). Homotopy Analysis Method and the Examples of Its Applications, Wyd. Pol. Śl., Gliwice (in Polish).
  24. W. Sumelka, (2016). Fractional calculus for continuum mechanics – anisotropic non-locality, Bull. Pol. Ac.: Tech 64 (2), 361–372.
  25. M. Turkyilmazoglu, (2010). A note on the homotopy analysis method, Appl. Math. Lett. 23, 1226–1230.
    [CROSSREF]
  26. M. Turkyilmazoglu, (2010). Convergence of the homotopy analysis method, arXiv, 1006.4460v1.
  27. K. Vajravelu, R.A. Van Gorder, (2012). Nonlinear Flow Phenomena and Homotopy Analysis. Fluid Flow and Heat Transfer, Springer/Higher Education Press, Berlin/Beijing.
    [CROSSREF]
  28. R.A. Van Gorder, (2012). Control of error in the homotopy analysis of semi-linear elliptic boundary value problems, Numer. Algor. 61, 613–629.
    [CROSSREF]
  29. R.A. Van Gorder and K. Vajravelu, (2009). On the selection of auxiliary functions, operators, and convergence control parameters in the application of the homotopy analysis method to nonlinear differential equations: a general approach, Commun. Nonlinear Sci. Numer. Simulat. 14, 4078–4089.
    [CROSSREF]
  30. X. Zhang, B. Tang, and Y. He, (2011). Homotopy analysis method for higher-order fractional integrodifferential equations, Comput. Math. Appl. 62, 3194–3203.
    [CROSSREF]
  31. M. Zurigat, S. Momani, Z. Odibat, and A. Alawneh, (2010). The homotopy analysis method for handling systems of fractional differential equations, Appl. Math. Modelling 34, 24–35.
    [CROSSREF]

EXTRA FILES

COMMENTS