Silesian University of Technology
Subject: Architecture, Civil Engineering, Engineering, Environmental
ISSN: 1899-0142
SEARCH WITHIN CONTENT
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
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.
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].
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
with the initial conditions
where N_{1} is the linear operator, N_{2} is the non-linear operator, y_{a} and v_{a} are the given numbers, while y is the unknown function. Let us denote
In the first step, we specify the homotopy operator ℋ as:
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
Substituting p=0 we get, so L(Φ(t;0) – y_{0}(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 y_{0} to y).
Taking the Maclaurin series of function Φ (with respect to parameter p) we get
where y_{0}(t) = Φ(t; 0) and
If the above series is convergent for p = 1, then we obtain the required solution
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
In order to determine the form of function y_{m} 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):
where $${\overline{y}}_{m-1}=\{{y}_{0}(t),{y}_{1}(t),\dots ,{y}_{m-1}(t)\}$$ and
and
In view of the above definition of operator N (equation 3) we have for m= 1:
whereas for m ≥ we get
Thus, the equation 10 for m= 1 takes the form
whereas for m ≥ 2 it takes the form
After selecting operator L the above two equations enable to determine recursively the successive functions y_{m} 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
Then we get for m=1:
and for m ≥ 2:
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 N_{2}.
As the initial approximation the best is to select a function satisfying the given initial conditions 2, for example in the form
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 y_{0} as the solution of linear equation L(y_{0}) = 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
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
And next, the optimum value of the convergence control parameter is obtained by determining the minimum of this squared residual.
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/m^{3}, is equal to
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 m′ dx (Figure 1a), is equal to
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 $$0\le x\le {\displaystyle \frac{1}{2}}l$$ is expressed by relation
where y_{m} is the largest deflection of the beam (Figure 1a).
Velocity of every point of coordinate x is equal to
where $${\dot{y}}_{m}={\displaystyle \frac{d{y}_{m}}{dt}}$$ and the kinetic energy of the whole beam is equal to
By substituting 24 into 25 and by calculating the integrals we get
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
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
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 $$\frac{b\text{}{h}^{3}}{12}}=5.4\xb7{10}^{-3}{\text{m}}^{4$$ 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
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
where $${\omega}_{0}=\sqrt{{\displaystyle \frac{k}{m}}=140.8573}\text{rad/s}$$ is the natural frequency of non-damping vibrations. The above equation must be completed with the initial conditions
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 y_{a}=0.01 m and v_{a}=0 m/s, the exact solution is given by function
If we take as the initial approximation the function satisfying the initial conditions
then in the first steps of the method we get
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
where y_{e} 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 |y_{e}(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.
In the second example we execute the calculations for the changed initial conditions by taking y_{a} = 0.01 m and v_{a}=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
As the initial approximation we take the function fulfilling the initial conditions
Then we obtain successively
Figure 5 displays the squared residual for n = 9.
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.
The initial condition
of the free vibrations can be achieved by dropping on the beam from height h the body of mass m_{1} (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 m_{1} = 100 kg and h= 0.5 m velocity of the beater is equal to
where g is the acceleration of gravity. Basing on the momentum conservation principle the velocity v_{a} of the system after the strike is equal to
After the impact of the plastic character the motion of the oscillator is executed according to the equation
the solution of which, after taking into account the initial conditions
is given by the function
where
Now we deal with the solution of equation 35 with conditions 36, for the following parameters m_{1} = 100 kg, m = 1633 kg, k= 32.4 MN/m, g=9.81 m/s ^{2}, and the initial conditions y_{a}=0 m and v_{a}=0.181 m/s. This time by using our method we get in turn
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.
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
The plus sign preceding $$\widehat{k}$$ 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
with the initial conditions
Referring to the previous example when we assumed the mass m= 1633 kg and stiffness k= 32.4 MN/m, the parameter $$\widehat{k}$$ can take the values about 10^{9} with the amplitude of variable y up to 0.002 m. Certainly the value of $$\widehat{k}$$ 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:
In our case the physically real value of c (determined from the engineering relations, ζ= 0.05 – damping ratio) is equal to
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, $$\widehat{k}={10}^{9}$$ and v_{a}=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 y_{0} the function satisfying the initial conditions
we get, by using the homotopy analysis method, the function
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
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).
In the next example we solve the equation 41 with initial conditions 40 for m = 1633 kg, k = 32.4 MN/m, $$\widehat{k}={10}^{9}$$, c= 23·10^{3} kg/s and v_{a} = 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 y_{0}, that is
we get successively the functions
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
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}.
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.
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