Numerical Study of the Vibrations of Beams with Variable Stiffness under Impulsive or Harmonic Loading

Moussa Sali^1,2*, Fabien Kenmogne², Jean Bertin Nkibeu³, Abdou Njifenjou^4,5

¹Laboratory of Materials, Mechanics and Civil Engineering, National Higher Polytechnic School of Maroua, University of Maroua, Maroua, Cameroon.
²Department of Civil Engineering, Advanced Teachers Training College of the Technical Education, University of Douala, Douala, Cameroon.
³Laboratory Engineering Civil and Mechanics, National Advanced School of Engineering, University of Yaoundé 1, Yaoundé, Cameroon.
⁴Laboratory of Energy, Materials, Modeling and Method (E3M), National Advanced School of Engineering, University of Douala, Douala, Cameroon.
⁵Laboratory of Mathematical Engineering and Information Systems, National Advanced School of Engineering, University of Yaoundé 1, Yaoundé, Cameroon.

DOI: 10.4236/wjet.2024.122026                

Abstract

The behavior of beams with variable stiffness subjected to the action of variable loadings (impulse or harmonic) is analyzed in this paper using the successive approximation method. This successive approximation method is a technique for numerical integration of partial differential equations involving both the space and time, with well-known initial conditions on time and boundary conditions on the space. This technique, although having been applied to beams with constant stiffness, is new for the case of beams with variable stiffness, and it aims to use a quadratic parabola (in time) to approximate the solutions of the differential equations of dynamics. The spatial part is studied using the successive approximation method of the partial differential equations obtained, in order to transform them into a system of time-dependent ordinary differential equations. Thus, the integration algorithm using this technique is established and applied to examples of beams with variable stiffness, under variable loading, and with the different cases of supports chosen in the literature. We have thus calculated the cases of beams with constant or variable rigidity with articulated or embedded supports, subjected to the action of an instantaneous impulse and harmonic loads distributed over its entire length. In order to justify the robustness of the successive approximation method considered in this work, an example of an articulated beam with constant stiffness subjected to a distributed harmonic load was calculated analytically, and the results obtained compared to those found numerically for various steps (spatial h and temporal $\overline{\tau }$ ) of calculus, and the difference between the values obtained by the two methods was small. For example for ($h=1/8$ , $\overline{\tau }=1/64$ ), the difference between these values is 17%.

Keywords

Successive Approximations Method, Direct Integration, Differential Equations, Beams of Variable Stiffness, Quadratic Parabola, Impulse and Harmonic Loads

1. Introduction

Beam structures are commonly used in many technological and industrial sectors, particularly in the construction industry [1] . Variable stiffness beams are frequently used to save beam materials or to make beams lighter [2] [3] , particularly when they are used in aeronautical or automobile manufacturing as well as in Civil engineering [4] . In civil engineering, they are used as foundation beams, columns of certain buildings, electrical transport pillars and chimneys. By carefully designing the stiffness distribution along the beam, a substantial increasing in the swelling of the beam can be achieved over its constant stiffness [5] .

During the operating process, these structures are subjected to both the action of static and dynamic loads, and it is important to develop appropriate models capable of capturing the complicated behavior of variable stiffness of non-uniform beams under various loading conditions [6] , these models must be accurate, efficient and simple to implement. It is obvious that analytical solutions offer reliability for benchmarking purposes and as preliminary design tools [6] . Numerical methods like the finite difference method [7] [8] , the differential transformed method… [9] , have been developed for the solution of complex problems. Despite the practical importance of certain models found in the literature (notably in [10] - [18] ), works to find improved formulations remain open and deserve particular attention.

The calculation of the structures in civil engineering requires the implementation of increasingly sophisticated numerical tools and methods to model the mechanical behavior and take into account the specificities of these structures. Their calculation by analytical methods remains very tedious and voluminous (see for example [10] [12] [16] [17] [18] [19] ), this is why numerical methods (see [8] [11] [14] [20] [21] [22] [23] ) are necessary and more efficient. Among them, one has the finite element method which is the most widely used, although it presents several difficulties such as forming the stiffness matrix and tightening the mesh around specific areas. In addition to this method, one has the finite difference method (MDF) which uses fictitious points in the boundary zones, and the results are imprecise for complex problems [8] . Some authors have made dazzling developments in the numerical resolutions of partial differential equations [24] - [28] . For this, some of them use the Caputo fractional operator [24] [27] [28] , while others call on the Riesz-Feller operator and use the Jacobi spectral collocation method combined with the trapezoidal rule to reduce Lévy-Feller fractional advection dispersion equation to a system of algebraic equations [25] [26] . Other interesting approaches, similar to the finite difference method, have been developed to solve similar problems [11] [15] [20] [21] [29] [30] , among which one is the successive approximation method (SAM). As outlined in [31] , the dynamic problems are solved by decomposing the motion into eigenmodes or by direct integration methods applied on differential equations of motion along the time axis. In the present work, we use SAM because, as we will see further, the implementation of its algorithm is simple. Furthermore, it is not necessary to first construct the stiffness matrix (elementary or global). It is reliable, stable and accurate compared to reference results. In addition, the calculation time is reduced. The application of the SAM requires knowledge of ordinary differential equations or higher order partial differential equations of the models, which must be transformed into a system of 2 or lower order differential equations or partial differential equations. Let us outline that the beams which will be examined in this work are exclusively beams whose section is a continuous function of the position. The differential equation of the oscillations of a bent beam of variable section with a distributed mass is generally in the form [33] :

$\frac{{\partial }^2}{\partial X^2}\left[EJ\left(X\right)\frac{{\partial }^{2}W}{\partial X^2}\right]+\mu \left(X\right)\frac{\partial W^2}{\partial t^2}+c\frac{\partial W}{\partial t}=q\left(X,t\right);$ (1)

where E is the Young’s modulus, J(X) is the quadratic moment of the cross section of the beam in a given section, μ(X) is the linear mass of the beam in a given section, c is the damping coefficient of the material of the beam, and q(X,t) the dynamic load in a given section. In order to solve Equation (1) by the SAM, after discretizing this equation in the integration domain, the parabolic spline function is used to approximate it of order two. The unknown function, and its first and second derivatives are continuous in the elements of the mesh but assumed can be discontinuous at the boundaries between the elements. Unlike the traditional finite difference method, the successive approximation method does not use fictitious nodes to deal with boundary conditions.

This paper aims to show how SAM could be used to develop algorithms for calculating isotropic beams of constant or variable thickness, subjected to the action of impulsive or harmonic forces. Thus, the paper is organized as follows. The first section performs the methodology and tools which is divided into three subsections. In the first one, we present the direct integration technics used to solve the differential equations governing the models. The second is devoted to modeling the equations of the models, and to the description of the boundary conditions, a set of algebraic equations are thus created. In the last subsection, we implement the calculation algorithm. Then we devote section 3 to the validation of the approach considered based on the method of successive approximations. For this purpose, tests are carried out on some well-known problems. Finally, in section 4, the conclusion is presented.

2. Methodology and Tools

2.1. Method of Successive Approximations (MSA) for Solving of Partial Differential Equations of Beams

When analyzing dynamic equilibrium equations, direct integration methods are usually used [32] . Square and cubic spline methods were used to substitute the desired function into the dynamic equation [33] , while the SAM [34] was also used. The main idea of the SAM consists of substituting the desired function and its derivatives with a polynomial (spline) of the same type, for example the cubic spline. These equations are obtained when we link the domain of integration of partial differential equations and in the limits of each element we substitute the latter by splines. The SAM allows solving problems relating to a system of second-order differential equations with partial or simple derivatives. Unlike the traditional finite difference method, the SAM does neither need the creation of fictitious points to describe the support conditions, nor the refinement of the mesh in specific areas. The SAM makes it possible to successfully solve problems of elastic calculation of systems of bars, plates, shells, as well as nonlinear problems. In this paper, we use a parabolic spline to approximate the dynamic equations over time. In the follow we will consider this issue in more detail. If the desired function is a function of three variables $F\left(x,y,t\right)$ , then numerical integration can be performed, representing the problem as two-dimensional and three-dimensional. Considering the time axis t as one of the coordinate axes, we obtain a three-dimensional problem, where $0\le x\le a$ , $0\le y\le b$ , $0\le t\le t_1$ . Considering the region $0\le x\le a$ , $0\le y\le b$ , at each time level, we obtain a two-dimensional problem, with $0\le t\le \infty$ . In the follow, we will focus on the two-dimensional formulation of the problem (for a beam—one-dimensional element), in which on the contrary to three-dimensional formulation, the time for studying the dynamic process is not limited, and the two-dimensional formulation of matrices are also used.

Considering the deflection as the required function $W\left(x,y,t\right)$ , for an approximate description of its change in time at $x=const$ , $y=const$ we apply the algebraic polynomial

$W\left(t\right)=a_0+a_{1}t+a_2{t}^2+\cdots +a_n{t}^n,$ (2)

which has $n+1$ coefficients and, by appropriately selecting them, $n+1$ conditions will be satisfied.

· Case $n=2$ :

For the case $n=2$ , one has:

$W\left(t\right)=a_0+a_{1}t+a_2{t}^2,W^t=\frac{\text{d}W}{\text{d}t}=a_1+2a_{2}t;W^{tt}=\frac{{\text{d}}^{2}W}{\text{d}{t}^2}=2a_2$ (3)

In which $W^t$ is the speed, and $W^{tt}$ the acceleration. The coefficients $a_0,a_1,a_2$ are determined from the initial conditions.

Let suppose that for $t=t_i$ $W\left(t_i\right)=W_i$ , and for $t=t_{i-1}$ $W\left(t_{i-1}\right)=W_{i-1}$ , $W^t\left(t_{i-1}\right)=W_{i-1}^t$ , then one can express the coefficients $a_0,a_1,a_2$ as functions of $W_i,W_{i-1},W_{i-1}^t$ .

For $t=t_{i-1}$ from Equation (2) we get:

$W\left(t_{i-1}\right)=W_{i-1}=a_0+a_1{t}_{i-1}+a_2{t}_{i-1}^2,W^t\left(t_{i-1}\right)=W_{i-1}^t=a_1+2a_2{t}_{i-1}.$ (4)

For $t=t_i$ from Equation (1) we obtain:

$W\left(t_i\right)=W_i=a_0+a_1{t}_i+a_2{t}_i^2.$ (5)

Considering that $t_i=t_{i-1}+{\tau }_i$ , where ${\tau }_i$ -is the time step, solving Equations (4) and (5), simultaneously we find:

$\{\begin{array}{l}{a}_2=\frac{1}{{\tau }_i^2}\left(W_i-W_{i-1}\right)-\frac{W_{i-1}^t}{{\tau }_i};\\ a_1=-\frac{2t_{i-1}}{{\tau }_i^2}\left(W_i-W_{i-1}\right)+W_{i-1}^t\left(1+\frac{2t_{i-1}}{{\tau }_i}\right);\\ a_0=W_{i-1}\left(1-\frac{t_{i-1}^2}{{\tau }_i^2}\right)+W_i\frac{t_{i-1}^2}{{\tau }_i^2}-W_{i-1}^t\left(t_{i-1}+\frac{t_{i-1}^2}{{\tau }_i}\right).\end{array}$ (6)

Substituting the found values of the coefficients $a_0,a_1,a_2$ in Equation (2), taking into consideration Equation (3), we will get at $t=t_i$

$W_i^t=-W_{i-1}^t-\frac{2}{{\tau }_i}\left(W_{i-1}-W_i\right),i=1,2,3,\cdots$ (7)

$W_i^{tt}=-\frac{2}{{\tau }_i}{W}_{i-1}^t-\frac{2}{{\tau }_i^2}\left(W_{i-1}-W_i\right);i=1,2,3,\cdots$ (8)

if nod 1 is taken as the beginning. Thus, based on the square parabola approximation, we obtain recurring formulas to determine the speed and acceleration at each time layer. Let us consider the application of the obtained Equations (7) and (8) for a specific example of a system with one degree of freedom. The differential equation of free oscillations in the form $\ddot{W}+{\omega }^{2}W=0$ at $W_1=0$ admits as a solution $W={\stackrel{\dot{}}{W}}_1\sin \left(\omega t\right)$ , where ${\stackrel{\dot{}}{W}}_1=W_1^t$ is a given initial velocity. Suppose $\omega =W=1$ , then $\ddot{W}+W=0$ , otherwise considering Equation (3):

$W^{tt}+W=0$ (9)

and $W=\sin \left(t\right)$ , $W^{tt}=-\sin \left(t\right)$ . Let write Equations (7) and (8) while $\tau =cont$ in another form:

$V_i=-V_{i-1}+2\left(W_i-W_{i-1}\right),$ (10)

$a_i=2\left(W_i-W_{i-1}-V_{i-1}\right);$ (11)

where $V_i=W_i^t\tau$ , $a_i=W_i^{tt}{\tau }^2$ . By rewriting Equation (9) at $t=t_i$ taking into consideration Equation (10) and Equation (11) leads to:

$W_i=\frac{2}{2+{\tau }^2}\left(V_{i-1}+W_{i-1}\right).$ (12)

By solving Equations (10)-(12) simultaneously when dividing a quarter of the period into 32 equal parts along the length eighth equal parts ( $\tau =\frac{\pi }{2n}$ , where n-the number of partitions), with the initial conditions: $W_1=W_1^{tt}=0$ ; $W_1^t=\frac{V_1}{\tau }=1$ (international units system), we obtain the following results at $t=\frac{\pi }{2}$ : $W=0.98037$ , $W^t=0.03692$ , $W^{tt}=-0.98038$ , while the exact values are $W=1$ , $W^t=0$ , $W^{tt}=-1$ . Thus we have satisfactory accuracy.

2.2. Problem Model Equation and Boundary Conditions

2.2.1. Problem Formulation and Model Equation

Let us consider in this subsection the dynamical behavior of a damping elastic beam with variable stiffness under the influence of an arbitrary dynamic load, which is governed by the following partial differential equations in dimensionless quantities describing transverse vibrations [33] .

$\{\begin{array}{l}\frac{{\partial }^{2}m}{\partial {\xi }^2}=-\left(p-\bar{\mu }{v}^{\bar{t}\bar{t}}-\bar{c}{v}^{\bar{t}}\right)\\ \frac{{\partial }^{2}m}{\partial {\xi }^2}=-\gamma m\end{array}$ (13)

In which $v^{\bar{t}}=\frac{\partial v}{\partial \bar{t}}$ , $v^{\bar{t}\bar{t}}=\frac{{\partial }^{2}v}{\partial {\bar{t}}^2}$ . Where the parameters are linked to original variables as:

$\xi =\frac{Z}{l}$ ; $m=\frac{M}{q_0{l}^2}$ ; $v=\frac{WE{J}_0}{q_0{l}^4}$ ; $p=\frac{q\left(Z,t\right)}{q_0}$ ; $\bar{t}=\frac{t}{l^2}\sqrt{\frac{E{J}_0}{{\mu }_0}}$ ; $\bar{\mu }=\frac{\mu \left(z\right)}{{\mu }_0}$ ; $\gamma =\frac{E{J}_0}{EJ\left(Z\right)}$ ; $\bar{c}=\frac{c{l}^2}{\sqrt{E{J}_0{\mu }_0}}$ . (14)

$\mu \left(z\right)$ being the mass per unit length of the beam; ${\mu }_0$ the linear mass of the beam at its fixed point, and c the dissipation coefficient.

Let us solve Equation (13) using the SAM, which was already used in references [11] and [33] to solve problem with constant stiffness. Remembering this method, and accounting to both Equations (7) and (8), one has;

$\{\begin{array}{l}{v}_{ij}^{\bar{t}}=-v_{i-1j}^{\bar{t}}-\frac{2}{{\bar{\tau }}_i}\left(v_{i-1j}-v_{ij}\right),\\ v_{ij}^{\bar{t}\bar{t}}=-\frac{2}{{\bar{\tau }}_i}{v}_{i-1j}^{\bar{t}}-\frac{2}{{\bar{\tau }}_i^2}\left(v_{i-1j}-v_{ij}\right),\end{array}$ (15)

with $i=1,2,3,\cdots$ which is measured along the axis $\bar{t}$ , while $j=1,2,3,\cdots ,\left(n-1\right)$ is measured along the axis ξ. ${\bar{\tau }}_i=\frac{{\tau }_i}{l^2}\sqrt{\frac{E{J}_0}{{\mu }_0}}$ is the dimensionless time step; $v,v^{\bar{t}},v^{\bar{t}\bar{t}}$ are the dimensionless deflection, speed and acceleration, respectively.

The difference equation of the SAM, which approximates the differential Equation (13) on a uniform grid $\bar{t}={\bar{t}}_i$ in the absence of discontinuities $m,v,p$ , and derivatives $v,p$ , is obtained by replacing p by $p-\bar{\mu }{v}^{\bar{t}\bar{t}}-\bar{c}{v}^{\bar{t}}$ , leading to:

$\begin{array}{l}{m}_{ij-1}-2m_{ij}+m_{ij+1}+h\Delta m_{ij}^{\xi }\\ =-\frac{h^2}{12}\left[{\left(p-\bar{\mu }{v}^{\bar{t}\bar{t}}-\bar{c}{v}^{\bar{t}}\right)}_{ij-1}+10{\left(p-\bar{\mu }{v}^{\bar{t}\bar{t}}-\bar{c}{v}^{\bar{t}}\right)}_{ij}+{\left(p-\bar{\mu }{v}^{\bar{t}\bar{t}}-\bar{c}{v}^{\bar{t}}\right)}_{ij+1}\right].\end{array}$ (16)

The finite difference equation approximating the second line of Equation (13) following [6] can be written as:

$v_{ij-1}-2v_{ij}+v_{ij+1}=-\frac{h^2}{12}\left[{\gamma }_{j-1}{m}_{ij-1}+10{\gamma }_j{m}_{ij}+{\gamma }_{j+1}{m}_{ij+1}\right]-\frac{{\gamma }_j{h}^3}{12}\Delta m_{ij}^{\xi }$ (17)

Accounting the system the set of Equation (15) into Equation (16) for a constant time step ( ${\bar{\tau }}_i=\bar{\tau }$ ), leads to the recurrent equation for a regular point j in i^th time layer as follows:

$\begin{array}{l}\frac{6{\bar{\tau }}^2}{h^2}\left(m_{ij-1}-2m_{ij}+m_{ij+1}\right)+\frac{6{\bar{\tau }}^2}{h^2}\Delta m_{ij}^{\xi }-\left(d_{j-1}{v}_{ij-1}+10d_j{v}_{ij}+d_{j+1}{v}_{ij+1}\right)\\ =-\left(d_{j-1}{v}_{i-1j-1}+10d_j{v}_{i-1j}+d_{j+1}{v}_{i-1j+1}\right)-\bar{\tau }\left(e_{j-1}{v}_{i-1j-1}^{\bar{t}}+10e_j{v}_{i-1j}^{\bar{t}}+e_{j+1}{v}_{i-1j+1}^{\bar{t}}\right)\\ -\frac{{\bar{\tau }}^2}{2}\left(p_{ij-1}+10p_{ij}+p_{ij+1}\right),\left(i=1,2,3,\cdots \right),\left(j=1,2,3,\cdots \right)\end{array}$ (18)

with $d=\bar{\mu }+\bar{\tau }+\bar{c}$ ; $e=\bar{\mu }+\frac{\bar{\tau }}{2}\bar{c}$ . Both Equations (16) and (18) are solved together with the boundary and initial conditions of the problem applied on the deflection v, the rotation angle $\frac{\partial v}{\partial \xi }$ , the bending moment m, and the shear force $\frac{\partial m}{\partial \xi }$ . The problem is solved here for certain variants of boundary conditions which are described as follows.

2.2.2. Boundary Conditions

· Beam simply supported at left end

Let us write Equation (13) for node j of the hinged support of the beam at the i^th time layer as:

$v_{ij}=v_{ij\left(0\right)}$ , $m_{ij}=m_{ij\left(0\right)}$ , (19)

where $v_{ij\left(0\right)}$ and $m_{ij\left(0\right)}$ are the deflection and bending moment at the endpoint, respectively, which generally are functions of $\bar{t}$ .

· Beam with free left-hand end

Let’s rewrite the boundary conditions of the free end of the beam for node j of the beam at the i^th temporal layer as follows:

$m_{ij}=m_{ij\left(0\right)}$ , $\frac{\partial m}{\partial \xi }=m_{ij\left(0\right)}^{\xi }$ . (20)

$m_{ij\left(0\right)}^{\xi }$ being the shear force at the edge point. In order to approximate the second member of Equation (20) by the difference equations of SAM, the equation at nod j at the i^th time layer is given by [15] :

$\begin{array}{l}{\left(\frac{\partial m}{\partial \xi }\right)}_{ij}=m_{ij}^{\xi }\\ =\frac{1}{h}\left(-m_{ij}+m_{ij+1}\right)+\frac{h}{12}\left[5{\left(p-\bar{\mu }{v}^{\bar{t}\bar{t}}-\bar{c}{v}^{\bar{t}}\right)}_{ij}+{\left(p-\bar{\mu }{v}^{\bar{t}\bar{t}}-\bar{c}{v}^{\bar{t}}\right)}_{ij+1}\right]\\ +\frac{h^2}{12}{\left(p-\bar{\mu }{v}^{\bar{t}\bar{t}}-\bar{c}{v}^{\bar{t}}\right)}_{ij}^{\xi },\end{array}$ (21)

in which by the superscript ξ, we mean the partial derivative with respect to variable ξ. The above equation can be discretized to give:

$\begin{array}{c}{m}_{ij}^{\xi }=\frac{1}{h}\left(-m_{ij}+m_{ij+1}\right)+\frac{h}{12}\left(5p_{ij}+p_{ij+1}\right)+\frac{h^2}{12}{p}_{ij}^{\xi }-\frac{h}{12}\left(5{\bar{\mu }}_j{v}_{ij}^{\bar{t}\bar{t}}+{\bar{\mu }}_{j+1}{v}_{ij+1}^{\bar{t}\bar{t}}\right)\\ -\frac{h^2}{12}{\bar{\mu }}_j^{\xi }{v}_{ij}^{\bar{t}\bar{t}}-\frac{h^2}{12}{\bar{\mu }}_j{\left(v^{\bar{t}\bar{t}}\right)}_{ij}^{\xi }-\frac{h}{12}\left(5{\bar{c}}_j{v}_{ij}^{\bar{t}}+{\bar{c}}_{j+1}{v}_{ij+1}^{\bar{t}}\right)-\frac{h^2}{12}{\bar{c}}_j{\left(v^{\bar{t}}\right)}_{ij}^{\xi }.\end{array}$ (22)

Let us express ${\left(v^{\bar{t}}\right)}_{ij}^{\xi }$ in (22) in terms of $v^{\bar{t}}$ in a square parabola:

${\left(v^{\bar{t}}\right)}_{ij}^{\xi }=\frac{1}{2h}\left(-3v_{ij}^{\bar{t}}+4v_{ij+1}^{\bar{t}}-v_{ij+2}^{\bar{t}}\right).$ (23)

Inserting ${\left(v^{\bar{t}}\right)}_{ij}^{\xi }$ , above into we get:

$\begin{array}{c}{m}_{ij}^{\xi }=\frac{1}{h}\left(-m_{ij}+m_{ij+1}\right)+\frac{h}{12}\left(5p_{ij}+p_{ij+1}\right)+\frac{h^2}{12}{p}_{ij}^{\xi }\\ -\frac{h}{12}\left[\left(3.5{\bar{\mu }}_j+h{\bar{\mu }}_j^{\xi }\right)v_{ij}^{\bar{t}\bar{t}}+\left(2{\bar{\mu }}_j+{\bar{\mu }}_{j+1}\right)v_{ij+1}^{\bar{t}\bar{t}}-0.5{\bar{\mu }}_j{v}_{ij+2}^{\bar{t}\bar{t}}\right]\\ -\frac{h}{12}\left[3.5{\bar{c}}_j{v}_{ij}^{\bar{t}}+\left(2{\bar{c}}_j+{\bar{c}}_{j+1}\right)v_{ij+1}^{\bar{t}}-0.5{\bar{c}}_j{v}_{ij+2}^{\bar{t}}\right].\end{array}$ (24)

Which can be rewritten taking into consideration the set of equations (15) as:

$\begin{array}{l}{m}_{ij}^{\xi }=\frac{1}{h}\left(-m_{ij}+m_{ij+1}\right)-\frac{h}{6{\bar{\tau }}^2}\left[\left(3.5d_j+h{d}_j^{\xi }\right)v_{ij}+\left(2d_j+d_{j+1}\right)v_{ij+1}-0.5d_j{v}_{ij+1}\right]\\ +\frac{h}{6{\bar{\tau }}^2}\left[\left(3.5d_j+h{d}_j^{\xi }\right)v_{i-1j}+\left(2d_j+d_{j+1}\right)v_{i-1j+1}-0.5d_j{v}_{i-1j+1}\right]\\ +\frac{h}{6\bar{\tau }}\left[\left(3.5e_j+h{e}_j^{\xi }\right)v_{i-1j}^{\bar{t}}+\left(2e_j+e_{j+1}\right)v_{i-1j+1}^{\bar{t}}-0.5e_j{v}_{i-1j+1}^{\bar{t}}\right]\\ +\frac{h}{12}\left(5p_{ij}+p_{ij+1}\right)+\frac{h^2}{12}{p}_{ij}^{,\xi }\text{where}i=1,2,3,\cdots \end{array}$ (25)

In which j is constant integer. To take into consideration the condition of the free edge of the beam in the above equation, it is necessary to substitute the specified the value of $m_{ij\left(0\right)}^{\xi }$ instead of $m_{ij}^{\xi }$ in Equation (25).

· Beam with left end rigidly fixed

The boundary conditions for nod j of the restrained edge of the beam in the i^th time layer read:

$v_{ij}=v_{ij\left(0\right)}$ , ${\left(\frac{\partial v}{\partial \xi }\right)}_{ij}=v_{ij}^{\xi }=v_{ij\left(0\right)}^{\xi }$ , (26)

where $v_{ij\left(0\right)}^{\xi }$ is the specified rotation angle at the edge point. The second member of Equation (26) can be discretize for nod j in the i^th time layer as follows [15] :

$v_{ij}^{\xi }=\frac{1}{h}\left(-v_{ij}+v_{ij+1}\right)+\frac{h}{12}\left(5{\gamma }_j+h{\gamma }_j^{\xi }\right)m_{ij}+\frac{h}{12}{\gamma }_{j+1}{m}_{ij+1}+\frac{h^2}{12}{\gamma }_j{m}_{ij}^{\xi }.$ (27)

Substituting (25) into (27), we get:

$\begin{array}{l}{v}_{ij}^{\xi }=\frac{h^3{\gamma }_j}{72{\bar{\tau }}^2}\left[\left(3.5d_j+h{d}_j^{\xi }+\frac{72{\bar{\tau }}^2}{h^4{\gamma }_j}\right)v_{ij}+\left(2d_j+d_{j+1}+\frac{72{\bar{\tau }}^2}{h^4{\gamma }_j}\right)v_{ij+1}-0.5d_j{v}_{ij+2}\right]\\ +\frac{h^3{\gamma }_j}{72{\bar{\tau }}^2}\left[\left(3.5d_j+h{d}_j^{\xi }\right)v_{i-1j}+\left(2d_j+d_{j+1}\right)v_{i-1j+1}-0.5d_j{v}_{i-1j+2}\right]\\ +\frac{h^3{\gamma }_j}{72{\bar{\tau }}^2}\left[\left(3.5e_j+h{e}_j^{\xi }\right)v_{i-1j}^{\bar{t}}+\left(2e_j+e_{j+1}\right)v_{i-1j+1}^{\bar{t}}-0.5e_j{v}_{i-1j+2}^{\bar{t}}\right]\\ +\frac{h}{12}\left[\left(4{\gamma }_j+h{\gamma }_j^{\xi }\right)m_{ij}+\left({\gamma }_j+{\gamma }_{j+1}\right)m_{ij+1}\right]+\frac{h^3{\gamma }_j}{144}\left[5p_{ij}+p_{ij+1}+h{p}_{ij}^{\xi }\right],i=2,3,\cdots \end{array}$ (28)

In order to take into account the pinching conditions of the edge of the beam, it is necessary to substitute the specified value $v_{ij}^{\xi }\left(0\right)$ instead of $v_{ij}^{\xi }$ into Equation (28), which is written for the left edge of the beam. Otherwise for the right edge, the above equation is written in “mirror image” with the replacement of i+1, and i+2 by i-1, and i-2, respectively while $v^{\xi }$ is inverted.

2.3. Implementation of an Algorithm of the Calculation

If the two ends of the beam are rigidly fixed, the calculation reduces to the joint solution of systems of equations such as Equations (17), (18), (28), taking into account the initial conditions.

The finite difference Equation (17) and Equation (18) obtained above, written for all grid nodes, associated to the boundary conditions and accounting to the initial conditions allow to determine v and mat any design point of the beam. The set of Equations (17), (18), (27), (28), associated with the initial conditions form a closed system of algebraic equations, which will be solved using the iterative Gauss-Seidel algorithm.

For this purpose, the set of Equation (17) and Equation (18) is rewritten in the following form:

$v_{ij}=0.9v_{ij}+\frac{1}{20}\left(v_{ij-1}+v_{ij+1}\right)+\frac{h^2}{240}\left({\gamma }_{j-1}{m}_{ij-1}+10{\gamma }_j{m}_{ij}+{\gamma }_{j+1}{m}_{ij+1}\right)-\frac{h^3{\gamma }_j}{240}\Delta m_{ij}^{\xi };$ (29)

$\begin{array}{c}{m}_{ij}=0.9m_{ij}+\frac{1}{20}\left(m_{ij-1}+m_{ij+1}+h\Delta m_{ij}^{\xi }\right)-\frac{h^2}{120{\bar{\tau }}^2}\left(d_{j-1}{v}_{ij-1}+10d_j{v}_{ij}+d_{j+1}{v}_{ij+1}\right)\\ +\frac{h^2}{120{\bar{\tau }}^2}\left(d_{j-1}{v}_{i-1j-1}+10d_j{v}_{i-1j}+d_{j+1}{v}_{i-1j+1}\right)\\ +\frac{h^2}{120\bar{\tau }}\left(e_j{v}_{i-1j-1}^{\bar{t}}+10e_j{v}_{i-1j}^{\bar{t}}+e_{j+1}{v}_{i-1j+1}^{\bar{t}}\right)+\frac{h^2}{240}\left(p_{ij-1}+10p_{ij}+p_{ij+1}\right).\end{array}$ (30)

The deflection at the free end of the beam is determined by remembering Equation (25) as:

$\begin{array}{c}{v}_{ij}=\frac{\left(2d_j+d_{j+1}\right)}{\left(3.5d_j+h{d}_j^{\xi }\right)}{v}_{ij+1}+\frac{0.5d_j}{\left(3.5d_j+h{d}_j^{\xi }\right)}{v}_{ij+2}+\frac{6{\bar{\tau }}^2}{\left(3.5d_j+h{d}_j^{\xi }\right)}\left(m_{ij+1}-m_{ij}\right)\\ +v_{i-1j}+\frac{\left(2d_j+d_{j+1}\right)}{\left(3.5d_j+h{d}_j^{\xi }\right)}{v}_{i-1j+1}-\frac{0.5d_j}{\left(3.5d_j+h{d}_j^{\xi }\right)}{v}_{i-1j+2}\end{array}$

$\begin{array}{c}+\frac{\bar{\tau }\left(3.5e_j+h{e}_j^{\xi }\right)}{\left(3.5d_j+h{d}_j^{\xi }\right)}{v}_{i-1j}^{\bar{t}}+\frac{\bar{\tau }\left(2e_j+e_{j+1}\right)}{\left(3.5d_j+h{d}_j^{\xi }\right)}{v}_{i-1j+1}^{\bar{t}}-\frac{0.5\bar{\tau }{e}_j}{\left(3.5d_j+h{d}_j^{\xi }\right)}{v}_{i-1j+2}^{\bar{t}}\\ +\frac{0.5{\bar{\tau }}^2}{\left(3.5d_j+h{d}_j^{\xi }\right)}\left(5p_{ij}+p_{ij+1}+h{p}_{ij}^{\xi }\right)-\frac{6{\bar{\tau }}^2}{h\left(3.5d_j+h{d}_j^{\xi }\right)}{m}_{ij\left(0\right)}^{\xi }.\end{array}$ (31)

The bending moment at the embedded end of the beam is determined using Equation (28), leading to:

$\begin{array}{c}{m}_{ij}=-\frac{{\gamma }_j+{\gamma }_{j+1}}{4{\gamma }_j+h{\gamma }_j^{\xi }}{m}_{ij+1}-\frac{h^2{\gamma }_j}{12\left(4{\gamma }_j+h{\gamma }_j^{\xi }\right)}\left(5p_{ij}+p_{ij+1}+h{p}_{ij}^{\xi }\right)\\ +\frac{h^2{\gamma }_j}{6{\bar{\tau }}^2\left(4{\gamma }_j+h{\gamma }_j^{\xi }\right)}\left[\left(3.5d_j+h{d}_j^{\xi }+\frac{72{\bar{\tau }}^2}{h^4{\gamma }_j}\right)v_{ij}+\left(2d_j+d_{j+1}-\frac{72{\bar{\tau }}^2}{h^4{\gamma }_j}\right)v_{ij+1}-0.5d_j{v}_{ij+2}\right]\\ -\frac{h^2{\gamma }_j}{6{\bar{\tau }}^2\left(4{\gamma }_j+h{\gamma }_j^{\xi }\right)}\left[\left(3.5d_j+h{d}_j^{\xi }\right)v_{i-1j}+\left(2d_j+d_{j+1}\right)v_{i-1j+1}-0.5d_j{v}_{i-1j+2}\right]\\ -\frac{h^2{\gamma }_j}{6{\bar{\tau }}^2\left(4{\gamma }_j+h{\gamma }_j^{\xi }\right)}\left[\left(3.5e_j+h{e}_j^{\xi }\right)v_{i-1j}^{\bar{t}}+\left(2e_j+d{e}_{j+1}\right)v_{i-1j+1}^{\bar{t}}-0.5e_j{v}_{i-1j+2}^{\bar{t}}\right]+\frac{12}{h\left(4{\gamma }_j+h{\gamma }_j^{\xi }\right)}{v}_{ij\left(0\right)}^{\xi }.\end{array}$ (32)

The resulting formulas correspond to a regular grid with steps $\bar{\tau }$ (in time) and h (in space), d, e, and γ are variables. As outlined in [33] [34] , the iterative process converge weather the coefficient of unknowns in the above Equations (29)-(32) are all less than one. As shown in the above equations, this condition depends both on the choice of h and $\bar{\tau }$ , as well as on the values of and , which are variables. The sequence of the algorithm implementation is as follows. If $\bar{t}=0\left(i=1\right)$ $v_{ij}=v_{ij\left(0\right)}$ ; $v_{ij}^{\bar{t}}=v_{ij\left(0\right)}^{\bar{t}}$ .

For $\bar{t}=\bar{\tau }\left(i=2\right)$ the values $v_{2j}$ , $m_{2j}$ are calculated from the solution of the above set of equations. The variables $v_{2j}^{\bar{t}}$ are determined from Equation (17) at ${\bar{\tau }}_j=\bar{\tau }$ .

For the calculus, $v_{3j}$ , $m_{3j}$ , the $v_{2j}$ , $m_{2j}$ , $v_{2j}^{\bar{t}}$ found above are taken into account and the process is repeated sequentially.

Considering the dynamic behavior of a beam under the influence of a distributed instantaneous impulse, meaning that at the initial instant the system acquires the greatest speed [35] :

${\frac{\partial W\left(Z,t\right)}{\partial t}|}_{t=0}=\frac{S\left(Z\right)}{\mu \left(Z\right)}.$ (33)

where μ is the linear mass of the beam.

Let us suppose that in Equation (14) $q_0=\frac{S_0}{l^2}\sqrt{\frac{E{J}_0}{{\mu }_0}}$ , where $S_0$ is the instantaneous impulse uniformly distributed over the entire length of the beam.

Then from the second line of Equation (13) we have $\omega =\frac{W\sqrt{E{J}_0{\mu }_0}}{S_0{l}^2}$ ; $m=\frac{M}{S_0}\sqrt{\frac{{\mu }_0}{E{J}_0}}$ . We thus write down Equation (33) for dimensionless unknowns on one has

$\frac{\partial v}{\partial t}=\frac{\bar{S}}{\bar{\mu }}.$ (34)

where $\bar{S}=\frac{S\left(Z\right)}{S_0}$ is a dimensionless instantaneous impulse distributed according to an arbitrary law. So to calculate a beam with variable rigidity subjected to the action of an instantaneous impulse without taking damping into account, it suffices to set $\bar{c}=0$ , $p=0$ , $v_1^{\bar{t}}=\bar{S}/\bar{\mu }$ , $v_1=0$ in Equations (28)-(30).

The calculation of beams with variable rigidity subjected to the action of a short-term load without taking into account damping is carried out using Equations (26), (28), (29), (33) if we put $\bar{c}=0$ , $p\ne 0$ , and take into account the given initial conditions.

Based on the developed algorithms, personal computer programs in FORTRAN [35] were compiled for static and dynamic calculations of beams.

3. Results and Discussion

3.1. Beam Embedded on Its Two Supports under the Action of a Uniformly Distributed Instantaneous Impulse

As the first test problem, let us consider a beam embedded on its two supports, with constant stiffness under the action of a uniformly distributed instantaneous impulse, for which a solution is well known [33] , with as initial conditions kept as: $\bar{t}=0$ ; $v=0$ ; $\partial \left(\frac{v}{t}\right)=|v^{\bar{t}}=\frac{\bar{S}}{\bar{\mu }}=1$ . To this end Equations (29)-(31) are numerically implemented for varying time step ${\bar{\tau }}_i=\bar{\tau }$ and space step h, and for the following parameters:

${\gamma }_{j-1}={\gamma }_j={\gamma }_{j+1}=1$ ; $\bar{c}=0$ ; $d_{j-1}=d_j=d_{j+1}=1$ ;

$e_{j-1}=e_j=e_{j+1}=1$ ; $p_{ij}^{\xi }=p_{ij-1}=p_{ij}=p_{ij+1}=0$ . (35)

In Figure 1, the changes of deflection and bending moment are plotted in the middle of the span (that is for $\xi =0.5$ ) without taking into consideration the time dependency of the damping. Curves 1, 2, 3, 4, 5 are for ( $h=\frac{1}{10}$ , $\bar{\tau }=\frac{1}{100}$ ), ( $h=\frac{1}{12}$ , $\bar{\tau }=\frac{1}{144}$ ), ( $h=\frac{1}{16}$ , $\bar{\tau }=\frac{1}{256}$ ), ( $h=\frac{1}{20}$ , $\bar{\tau }=\frac{1}{400}$ ), ( $h=\frac{1}{30}$ , $\bar{\tau }=\frac{0.002}{\pi }$ ), respectively. As can be seen, the deflection (curve) are fairly similar for all h and $\bar{\tau }$ . Table 1 shows the maximum values of bending moment and deflection in the middle of the beam span for different partitions, while the best result is that of the last column, which coincides well with the results found in [25] .

3.2. Beam Articulated on Its One End

As the second test problem, one considers the articulated beam under the action of a uniformly distributed instantaneous impulse, with nonzero damping, that is

Figure 1. Variation of deflection and bending moment coefficients in the middle of the beam span ( $\xi =0.5$ ).

Table 1. Values of the bending moment and deflection coefficients in the middle of the beam ( $\xi =0.5$ ).

$\bar{c}\ne 0$ . The other parameters are given in Equation (35). The results found here are given in Table 2, giving as in the above Subsection the maximum values of the bending moment and deflection in the middle of the span of the beam for $\bar{c}=0.05$ and forvarying values of h and $\bar{\tau }$ . As one can see the results found differ little from those shown in Table 2, which is due to the effect of damping.

3.3. Beam with Constant Stiffness Embedded on Its Two Ends Subjected to the Action of a Uniformly Distributed Impulse

Let us now consider a beam of constant stiffness embedded on its two ends, subjected to the action of a uniformly distributed impulse. In this case, in addition to Equations (29)-(31) is also necessary to write Equation (32) with ${\gamma }_j={\gamma }_{j+1}=1$ , ${\gamma }^{\xi }=0$ , $v_{ij}=0$ , $p_{ij-1}=p_{ij}=p_{ij+1}=p_{ij}^{\xi }=0$ . As results, Figure 2 shows the curves of the deflection and bending moments as in Subsection 3.1, in which curves 1 and 3 describe the variation of bending moment and deflection in the middle of the span, respectively, while curve 2 shows the changes in the bending moment in the embedment at time and space steps $h=\frac{1}{30}$ , $\bar{\tau }=\frac{0.002}{\pi }$ . The maximum value of the dimensionless bending moment shown Table 2. Values of the bending moment and deflection coefficients in the middle of the beam ( $\xi =0.5$ ).

Figure 2. Variation of the bending moments (curve 1) and deflection coefficients (curve 3) in the middle and at the fixed end (curve 2) of the beam.

in Table 3 arising on the support is 2.402 for $h=\frac{1}{30}$ , $\bar{\tau }=\frac{0.002}{\pi }$ , corresponds approximately to that at $t=\frac{T}{14}$ , (where T is fundamental period of the oscillations of an articulated beam). It is obvious that the maximum value of the deflection (see Table 3) in the middle of the span is 0.0581 which occurs at the time $t=\frac{T}{12}$ , and the maximum value of the bending moment in the middle of the span is 1.342 and occurs later than $t=\frac{T}{12}$ . For an articulated beam, the maxima of the deflection and bending moment in the middle of the span appear approximately at $t=\frac{T}{5}$ , which is in agreement with results found in [33] , andaccordingly, are equal to $v_{\left(\max \right)}=0.1279$ , $m_{cp\left(\max \right)}=1.572$ (for $h=\frac{1}{30}$ , $\bar{\tau }=\frac{0.002}{\pi }$ ).

3.4. Beam of Constant Stiffness, Which Is Subjected to the Action of Harmonic Load Distributed along Its Entire Length

Here we consider an articulated beam of constant stiffness, subjected to the action of harmonic load distributed along its entire length as shown in Figure 3:

$p_{ij}=\sin \left(1.6\pi {\bar{t}}_i\right),$ (36)

Table 3. Values of the bending moment and deflection coefficients in the middle and at the ends of the beam.

Figure 3. Beam resting on two articulated supports and subjected to a harmonic load.

where ${\bar{t}}_i=\frac{t_i}{T}$ , $T=\frac{2\pi }{{\omega }_{\min }}$ , and ${\omega }_{\min }$ is the lowest natural frequency of the hinged beam.

Based on Equations (17) and (18) taking into consideration Equations (16), (22), (36) $\gamma =d=e=\bar{\mu }=1$ , $\Delta m^{\xi }=0$ , the maximum values of m_cp and v in time at the center of the considered beam are obtained. The obtained values as shown in Table 5 appear approximately at time $t=0.425T$ . The third and fourth lines of Table 4 give the maximum dimensionless bending moment and deflection in the middle of the beam span, respectively. The last column this same table shows the results obtained using the analytical method. Results close to them were obtained by the equations of SAM at $\bar{t}=0.425$ for $h=\frac{1}{4}$ , $\bar{\tau }=\frac{1}{16}$.

Similarly we estimated the error of the numerical solution of the problem from the integral condition of the beam equilibrium. For this purpose, the sum of dimensionless support reactions $r_1+r_2=1.0736$ , the integral sum of the inertial forces ${\displaystyle \sum J}=02369$ , and the resultant harmonic load taken along the entire length of the beam $P=0.8197$ at $h=\frac{1}{4}$ , $\bar{\tau }=\frac{1}{16}$ was determined. According to the results obtained, it follows that the solution error is 1.5%.

Figure 4 shows the curves of changes in dimensionless parameters m and v in the middle of the beam span at $t=0.50T$ . Curve 1 (dimensionless bending moment), curve 2 (dimensionless deflection) according to the results obtained at

Table 4. Values of the bending moment and deflection coefficients in the middle of the beam.

Figure 4. Curves of variation of deflection and bending moment coefficients in the middle of the beam span ( $\xi =0.5$ ).

$h=\frac{1}{4}$ , $\bar{\tau }=\frac{1}{16}$ .

To clarify the results obtained above, we will solve this problem using the Newmark method when approximating the equations of dynamics [31] . Let us write these equations in dimensionless form for nod j at the i^th time layer:

$\{\begin{array}{l}{v}_{ij}^{\bar{t}}=v_{i-1j}^{\bar{t}}+\bar{\tau }\left(1-\delta \right)v_{i-1j}^{\bar{t}\bar{t}}+\delta v_{ij}^{\bar{t}\bar{t}}\\ v_{ij}^{\bar{t}\bar{t}}=\frac{\left(v_{ij}-v_{i-1j}\right)}{\alpha {\bar{\tau }}^2}-\frac{v_{i-1j}^{\bar{t}}}{\alpha \bar{\tau }}+\left(1-\frac{0.5}{\alpha }\right)v_{i-1j}^{\bar{t}\bar{t}}+\delta v_{ij}^{\bar{t}\bar{t}}\end{array}$ (37)

where $\delta \ge 0.50$ ; $\alpha \ge 0.25{\left(0.5+\delta \right)}^2$ . If in the approximation of Equation (13) we use Equations (32) and (33), respectively, then for the case $\bar{c}=0$ we get:

$\begin{array}{l}{m}_{ij-1}-2m_{ij}+m_{ij+1}\\ =\frac{h^2\left({\bar{\mu }}_{j-1}{v}_{ij-1}+10{\bar{\mu }}_j{v}_{ij}+{\bar{\mu }}_{j+1}{v}_{ij+1}\right)}{12\alpha {\bar{\tau }}^2}-\frac{h^2\left({\bar{\mu }}_{j-1}{v}_{i-1j-1}+10{\bar{\mu }}_j{v}_{i-1j}+{\bar{\mu }}_{j+1}{v}_{i-1j+1}\right)}{12\alpha {\bar{\tau }}^2}\\ -\frac{h^2\left({\bar{\mu }}_{j-1}{v}_{i-1j-1}^{\bar{t}}+10{\bar{\mu }}_j{v}_{i-1j}^{\bar{t}}+{\bar{\mu }}_{j+1}{v}_{i-1j+1}^{\bar{t}}\right)}{12\alpha \bar{\tau }}\\ +\frac{h^2\left({\bar{\mu }}_{j-1}{v}_{i-1j}^{\bar{t}\bar{t}}+10{\bar{\mu }}_j{v}_{i-1j}^{\bar{t}\bar{t}}+{\bar{\mu }}_{j-1}{v}_{i-1j+1}^{\bar{t}\bar{t}}\right)}{12}-\frac{h^2\left(p_{ij-1}+10p_{ij}+p_{ij+1}\right)}{12}\end{array}$ (38)

which written in conjunction with Equation (17) taking into consideration Equations (18) and (37) and the initial conditions for the design nods beams

allow you to determine m, v. For the above problem with $\delta =\frac{1}{2}$ , $\alpha =\frac{1}{4}$ , $h=\bar{\tau }=\frac{1}{10}$ we get $m=0.2023$ , $v=0.02085$ which confirms the results obtained.

3.5. Beam with Two Embedded Ends

This is the beam from the previous example with two embedded ends. The maximum values of m_cp and v in the middle of the span and m_1,2 in the embedded ends of the beam are shown in Table 5. These values are obtained for $h=\frac{1}{4}$ , $h=\frac{1}{6}$ , $h=\frac{1}{8}$ , $h=\frac{1}{10}$ and for $\bar{\tau }=h^2$ . The last column of this table shows the results obtained by the analytical method at $h=\frac{1}{6}$ , $\bar{\tau }=\frac{1}{36}$ , leading that the difference between the results by SAM and the analytical ones for bending moments in the middle of the span and at the embedded edge of the beam is, 4%, 0.1%, respectively for deflection—4%.

For the case under consideration, the maximum values of the moments and deflections occur approximately at the time $t=0.25T$ . Curves 1, 2, 3 in Figure 5 correspond to the results given in the second column of Table 6, which reflect the change in m_cp, v, and m_1,2 at the period $t=0.50T$ .

3.6. Articulated Beam of Variable Stiffness, Which Is Subjected to the Action of an Instantaneous Impulse Uniformly Distributed over Its Entire Length

Consider an articulated beam of variable stiffness, which is subjected to the action of an instantaneous impulse uniformly distributed over its entire length.

Table 5. Values of the bending moment and deflection coefficients in the middle and at the ends of the beam.

Table 6. Values of the bending moment and deflection coefficients in the middle of the beam.

Figure 5. Variation of deflection and bending moment coefficients in the middle and at the ends of the beam.

The beam has rigidity:

$EJ\left(Z\right)=E{J}_0{\left(1+\frac{0.1Z}{l}\right)}^3,\mu \left(Z\right)={\mu }_0\left(1+\frac{0.1Z}{l}\right),$ (39)

where $E{J}_0$ and ${\mu }_0$ , are the bending stiffness and the mass of a unit length of the beam at $Z=0$ , respectively.

Table 6 gives the maximum values in time m_cp, v, in the middle of the span of an articulated beam. The maximum dimensionless bending moment equal to $m_{cp}=1.331$ at $h=\frac{1}{12}$ , $\bar{\tau }=\frac{1}{144}$ occurs at around the time $\bar{t}=0.10$ , and deflection equal $v=0.1176$ is obtained at around the same period.

The last and the penultimate lines of Table 6 give dimensionless deflection and bending moment in the middle of the beam span, respectively. Comparison of these results with the results of Table 1 at $h=\frac{1}{12}$ , $\bar{\tau }=\frac{1}{144}$ , shows that the bending moment in the middle of the span of a beam of variable stiffness is 3% greater than that of a beam of constant stiffness, and the deflection is 7% less than that of a beam of constant stiffness. These differences may come from inertia forces.

Figure 6 shows the curves of variation in dimensionless m_cp and v in the middle of the span of the beam at the period $t=0.25T$ . Curve (dimensionless bending moment), curve (dimensionless deflection) are plotted according to the results obtained at $h=\frac{1}{12}$ , $\bar{\tau }=\frac{1}{144}$ .

3.7. Beam from the Previous Example with Rigidly Fixed Ends

Here is the beam from the previous example with rigidly fixed ends.

Table 7 shows the maximum values in a time of dimensionless moments and deflection in the middle of the span and at the ends of the beam. The results were obtained for

Figure 6. Variation of deflection and bending moment coefficients in the middle of the beam span ( $\xi =0.5$ ).

Table 7. Values of the bending moment and deflection coefficients in the middle and at the ends of the beam.

( $h=\frac{1}{4}$ , $\bar{\tau }=\frac{1}{16}$ ); ( $h=\frac{1}{6}$ , $\bar{\tau }=\frac{1}{36}$ ); $(h=\frac{1}{8}$ , $\bar{\tau }=\frac{1}{64}$ ); ( $h=\frac{1}{10}$ , $\bar{\tau }=\frac{1}{100}$ ); ( $h=\frac{1}{12}$ , $\bar{\tau }=\frac{1}{144}$ ).

$m_{cp}=\frac{M_{cp}}{S_0}\sqrt{\frac{{\mu }_0}{E{J}_0}}$ , $v=\frac{W}{S_0{l}^2}\sqrt{EJ{\mu }_0}$ , $m_{1,2}=\frac{M_{1,2}}{S_0}\sqrt{\frac{{\mu }_0}{E{J}_0}}$ . (40)

The third and fourth columns of this table give dimensionless bending moment and deflection, respectively, in the middle of the beam span, and the fifth and sixth columns, dimensionless bending moments at the edges, at $\xi =0$ and $\xi =1$ .

Comparison of the results shown in both Table 6 and Table 7 at $h=\frac{1}{12}$ , $\bar{\tau }=\frac{1}{144}$ , it is obvious that the bending moment in the middle of the span ( $m_{cp}=1.069$ ) for a restrained beam is 19.7% less than that of a hinged beam, and the deflection $v=0.05290$ is 55% less.

Comparison of the results obtained here with the results of Table 3 at $h=\frac{1}{12}$ , $\bar{\tau }=\frac{1}{144}$ it appears that the bending moment in the middle of the span of an embedded beam of variable stiffness is 6% higher than that of a beam of constant stiffness, and the deflection is 3.6% less than that of a beam of constant stiffness. Figure 7 shows the time variation of the deflection and bending moments in the middle of the span and at the restrained edges of the considered beam under the action of the uniformly distributed momentum.

3.8. Pivotally Supported Beam

Here a pivotally supported beam is considered with stiffness μ defined by Equation (37) and loaded as (35). Table 8 shows the maximum values m_cp and v at

Figure 7. Variation of deflection and bending moment coefficients in the middle and at the ends of the beam.

the middle of the beam span. The maximum dimensionless bending moment equal to $m_{cp}=0.1973$ for $h=\frac{1}{10}$ , and $\bar{\tau }=\frac{1}{100}$ , which occurs approximately at $\bar{t}=0.425$ , which is 7.3% less than the corresponding moment in a beam of constant stiffness.

The deflection equal to $v=0.01838$ occurs at approximately the same time and is less than the corresponding deflection in a beam of constant stiffness by 18.3%. The last and penultimate lines of Table 8 give dimensionless deflection and bending moment in the middle of the beam span, respectively. Figure 8 shows the curves of variation of dimensionless variables m_cp and v in the middle of the beam at approximately the time $t=0.5T$ . Curve 1 (dimensionless bending moment), curve 2 (dimensionless deflection) which are plotted according to the results obtained for $h=\frac{1}{10}$ , $\bar{\tau }=\frac{1}{100}$ .

Let’s check the integral condition of the beam equilibrium. For this in a time $\bar{t}=0.425$ , the dimensionless support reactions $r_A+r_B=1.055$ , the integral sum of the inertial forces ${\displaystyle \sum J}=0.1411$ , and the resultant of the harmonic force taken along the entire length of the beam $P=0.9514$ were determined at $h=\frac{1}{10}$ , $\bar{\tau }=\frac{1}{100}$ . The results show that ${\displaystyle \sum J}+P=1.0925$ differs little from

Table 8. Values of the bending moment and deflection coefficients in the middle of the beam.

Figure 8. Curves of variation of deflection and bending moment coefficients in the middle of the beam span ( $\xi =0.\text{5}$ ).

$r_A+r_B=1.055$ . The numerical solution error is 1.6%.

Comparison of Table 9 with the results of Table 5 at h = 1, t = 2 shows that the bending moment in the middle of the span of a beam of variable stiffness is 1% less than that of a beam of constant stiffness, and the deflection is 13% less than a beam of constant stiffness.

3.9. Beam with Two Restrained Ends

Here we consider a beam with two restrained ends, stiffness and load which are defined according to Equation (35). Table 9 shows the maximum values in time of dimensionless moments and deflection in the middle of the span and at the ends of the beam.

The results in Table 8 were obtained for $h=\frac{1}{4}$ , $\bar{\tau }=\frac{1}{16}$ ; $h=\frac{1}{6}$ , $\bar{\tau }=\frac{1}{36}$ ; $h=\frac{1}{8}$ , $\bar{\tau }=\frac{1}{64}$ ; $h=\frac{1}{10}$ , $\bar{\tau }=\frac{1}{100}$ . In which the third and fourth columns give dimensionless bending moment and deflection in the middle of the beam span, respectively, and the fifth and sixth columns give dimensionless bending moments at the edges, respectively, at $\xi =0$ and $\xi =1$ . Figure 9 shows the curves

Table 9. Values of the bending moment and deflection coefficients in the middle and at the ends of the beam.

Figure 9. Curves of variation of deflection and bending moment coefficients in the middle and at the ends of the beam.

of time variation of dimensionless variables m_cp, v in the middle of the beam span and $m_1,m_2$ at its edges, respectively at $\xi =0$ , $\xi =1$ . The oscillation process of the beam was studied over some time $t=0.5T$ . Curve 1 (m_cp), curve 2 (v), curve 3 (m₁) and curve 4 (m₂) were plotted from the results obtained at $h=\frac{1}{10}$ , $\bar{\tau }=\frac{1}{100}$ .

The integral condition of the beam equilibrium was checked. At the moment in time $\bar{t}=0.25$ , ( $m_{cp}=0.04773$ , $v=0.002946$ , $m_1=-0.08642$ ,

$m_2=-0.09728$ ) dimensionless support reactions $r_A,r_B$ have been calculated. The sum of those reactions gives $r_A+r_B=0.9648$ and the integral sum of inertial forces ${\displaystyle \sum J}=0.0110$ , while the resultant harmonic force was determined, taken over the entire length of the beam $P=0.9150$ ( $h=\frac{1}{10}$ , $\bar{\tau }=\frac{1}{100}$ ). The results show that ${\displaystyle \sum J}+P=0.9260$ is somewhat different from $r_A+r_B=$ $0.9648$. The solution error is 4%.

Comparison of the results obtained here for $h=\frac{1}{10}$ , $\bar{\tau }=\frac{1}{100}$ with the results in Table 6 show that the bending moment in the middle of the span of a beam of variable stiffness is less than the corresponding bending moment in a beam of constant stiffness by 1.8%, and the deflection is 11.5% less than a beam of constant stiffness.

4. Conclusion

The SAM was applied to calculate isotropic beam structures with variable stiffness modeled by a system of partial differential equations with various boundary conditions. This SAM used the cubic parabola and the results were compared to those obtained using the quadratic parabola. The hypothesis of constant and variable thickness was made as well as that of impulsive and harmonic loads. The reliability and stability of the considered approach are confirmed by the resolution of nine dynamic calculation test problems for beams with constant and variable stiffness. The calculation was carried out for each beam without damping and with damping. As expected, the damping effect is noticeable. The results obtained by taking into account the damping are lower than those obtained without damping. We carried out the same calculations for beams with variable stiffness and we arrived at the same conclusions. The results obtained in the different examples show good convergence. The convergence of the solutions was studied numerically both on test examples and on new problems. The integral condition of the balance of the beam was also verified for each beam. The curves illustrate the variation of the bending moment and deflection coefficients at a fixed node during the time grid is constant. Different nodes were inspected and one can notice that on the curves, the value of the coefficients reaches the maximum in the time interval [0.15T - 0.25T], which clearly shows the numerical convergence of the method. Despite the fact that the proposed method is proven to be effective for beams, its application to the calculation of massive bodies (3D) is necessary, and constitutes a perspective for future investigations.

Data Availability

The authors confirm that the data supporting the finding of this study are available within the article and/or its supplementary materials.

Acknowledgements

Thanks to Professor Gabbassov Radek Fatihovich for his precious help.

Funding

Cameroonian Government.

Conflicts of Interest

We declare that there is no conflict of interest regarding the publication of this paper.

References

[1]	Jewett, J.L. and Carstensen, J.V. (2019) Topology-Optimized Design, Construction and Experimental Evaluation of Concrete Beams. Automation in Construction, 102, 59-67. https://doi.org/10.1016/j.autcon.2019.02.001
[2]	Manickam, G., Polit, O., Balaji, L., Kumar, M.A. and Dineshkumar, S. (2023) Variable-Stiffness Curved Laminated-Beams by Curvilinear Fibers with Arbitrarily Layup—Vibrational Features by Sine-Based Higher-Order Beam Model with Renewed-Constitutive Relations and Improved-Kinematics. Composite Structures, 15, Article ID: 117514. https://doi.org/10.1016/j.compstruct.2023.117514
[3]	Sayyad, A.S. and Ghugal, Y.M. (2017) Bending, Buckling and Free Vibration of Laminated Composite and Sandwich Beams: A Critical Review of Literature. Composite Structures, 171, 486-504. https://doi.org/10.1016/j.compstruct.2017.03.053
[4]	Yadav, S. (2019) Design and Static Structural Analysis of Aircraft Floor Beam. GRD Journals—Global Research and Development Journal for Engineering, 4, 6-10.
[5]	Bang, Y.D.Y., Sali, M., Noah, P.M., Kenmogne, F., Adoum, D.A. and Ateba, A. (2021) Prismatic Bending Shell with Variable Rigidity. International Journal of Advanced and Multidisciplinary Engineering Science, 4, 1-6.
[6]	Masjedi, P.K. and Weaver, P.M. (2021) Variable Stiffness Composite Beams Subject to Non-Uniformly Distributed Loads: An Analytical Solution. Composite Structures, 256, Article ID: 112975. https://doi.org/10.1016/j.compstruct.2020.112975
[7]	Katsikadelis, J.T. and Tsiatas, G.C. (2004) Non-Linear Dynamic Analysis of Beams with Variable Stiffness. Journal of Sound and Vibration, 270, 847-863. https://doi.org/10.1016/S0022-460X(03)00635-7
[8]	Sali, M., Frédéric, L., Hamandjoda, O. and Raidandi, D. (2018) Calculation of Plates on Elastic Foundation by the Generalized Equations of Finite Difference Method. The International Journal of Engineering and Science, 7, 32-38.
[9]	Kenmogne, F. (2015) Generalizing of Differential Transforms Method for Solving Nonlinear Differential Equations. Journal of Applied & Computational Mathematics, 4, Article ID: 1000196. https://doi.org/10.4172/2168-9679.1000196
[10]	Timoshenko, S.P. (1972) Theory of Elasticity. NAHUKOVA DUMKA, Kiev.
[11]	Gabassov, R.F., Gabassov, A.R. and Filatov, V.V. (2008) Numerical Construction of Discontinuous Solutions to the Problems of Constructive Mechanics. ABC Book Publishers, Moscow.
[12]	Vardanian, G.S., Andreev, V.I., Atapov, N.M. and Gorshkov, A.A. (1995) Material Resistance with Theory of Elasticity and Plasticity. ASB Publishing, Moscow.
[13]	Sali, M., Filatov, V.V. and Njomoue, P.A. (2019) Successive Approximation Method for the Calculation of Bent-Compressed Beams of Constant Stiffness. International Journal of Advances in Scientific Research and Engineering (IJASRE), 5, 7-12. https://doi.org/10.31695/IJASRE.2019.33242
[14]	Gabbassov, R.F. (1989) Numerical Solution of Problems of Construction Mechanics with Discontinuous Parameters. Master’s Thesis, Moscow Institute of Physics and Technology, Dolgoprudny.
[15]	Sali, M. and Chills, A.J. (2019) Static Calculation of Bending Beams of Varying Thickness Using the Finite Difference Equations of Successive Approximation Method. International Journal of Advances in Scientific Research and Engineering, 5, 240-248. https://doi.org/10.31695/IJASRE.2019.33176
[16]	Itskovich, M., Minin, L.S. and Vinokurov, A.I. (2001) Guide to Solving Problems on the Resistance of Materials. 3rd Edition, Higher School, Moscow.
[17]	Wilson, E.L. (2002) Three-Dimensional Static and Dynamic Analysis of Structures, a Physical Approach. Third Edition, Computers and Structures, Inc., Berkeley.
[18]	Ananin, A.I. (1969) On the Theory of Transverse Vibrations of Variable Section Bars. Izvestiya VUZov: Building and Architecture, No. 11, 45-50.
[19]	Korenev, B.G. and Bukeikhanov, S.R. (1969) On the Vibration of Bars of Variable Section. Structural Mechanics and Calculation of Structures, No. 6, 43-48.
[20]	Korenev, B.G., Laschenikov, B.Y., Leontiev, N.N., Gabbasov, R.F. and Szilard, R. (1984) Numerical Methods of Structural Mechanics. Volume I. Rod Systems. New Books Abroad, Series B, No. 3, 97-100.
[21]	Alexandrov, A.V. (1961) Numerical Solution of Linear Differential Equations Using the Differentiation Matrix. MIIT, Moscow, 253-266.
[22]	Jafari, H., Firoozjaee, M.A. and Johnston, S.J. (2020) An Effective Approach to Solving a System of Fractional Differential Equations. Alexandria Engineering Journal, 59, 3213-3219. https://doi.org/10.1016/j.aej.2020.08.015
[23]	Vijayakumar, V., Panda, S.K., Nisar, K.S. and Baskonus, H.M. (2021) Results on Approximate Controllability Results for Second-Order Sobolev-Type Impulse Neutral Differential Evolution Inclusions with Infinite Delay. Numerical Methods for Partial Differential Equations, 37, 1200-1221. https://doi.org/10.1002/num.22573
[24]	Odibat, Z. and Baleanu, D. (2020) Numerical Simulation of Initial Value Problems with Generalized Caputo-Type Fractional Derivatives. Applied Numerical Mathematics, 156, 94-105. https://doi.org/10.1016/j.apnum.2020.04.015
[25]	Shiri, B. and Baleanu, D. (2020) Collocation Method for Terminal Value Problems of Tempered Fractional Differential Equations. Applied Numerical Mathematics, 156, 385-395. https://doi.org/10.1016/j.apnum.2020.05.007
[26]	Sweilam, N.H. and Abdou Hassan, M.M. (2020) Efficient Method for Factional Lévy-Feller Advection-Dispersion Equation Using Jacobi Polynomials. Progress in Fractional Differentiation and Applications, 6, 115-128. https://doi.org/10.18576/pfda/060204
[27]	Abdullaev, O.K. (2020) Some Problems for the Degenerate Mixed Type Equation Involving Caputo and Atangana-Baleanu Operators Fractional Order. Progress in Fractional Differentiation and Applications, 6, 101-114. https://doi.org/10.18576/pfda/060203
[28]	Pachpatte, D.B. (2020) A Generalized Gronwall Inequality for Caputo Fractional Dynamic Delta operator. Progress in Fractional Differentiation and Applications, 6, 129-136. https://doi.org/10.18576/pfda/060205
[29]	Njifenjou, A., Donfack, H. and Moukouop-Nguena, I. (2013) Analysis on General Meshes of Discrete Duality Finite Volume Method for Subsurface Flow Problems. Computational Geosciences, 17, 391-415. https://doi.org/10.1007/s10596-012-9339-6
[30]	Filatov, V.V. (1999) The Calculation of Compressed-Curved Beams and Slabs on a Non-Continuous Elastic Base. Ph.D. Thesis, Moscow State University of Civil Engineering, Moscow.
[31]	Bate, K. and Wilson, E. (1982) Numerical Methods of Analysis and the Method of Finite Elements. Stroyizdat, Moscow.
[32]	Clough, R. and Penzien, D. (1979) Dynamics of Structures. Third Edition, Computers & Structures, Inc., University Ave., Berkeley.
[33]	Nikkhoo, A. and Kananipour, H. (2014) Numerical Solution for Dynamic Analysis of Semicircular Curved Beams Acted upon by Moving Loads. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, 228, 2314-2322. https://doi.org/10.1177/0954406213518908
[34]	Gabbasov, R.F. and Nizomov, D.N. (1985) Numerical Solution of Some Dynamic Problems of Structural Mechanics. Structural Mechanics and Calculation of Structures, No. 6, 51-54.
[35]	Daniel, D. and Dorn McCracken, W.S. (1977) Numerical Methods and Programming in FORTRAN. John Wiley and Sons, New York.