Numerical Study of the Vibrations of Beams with Variable Stiffness under Impulsive or Harmonic Loading ()
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] :
∂2∂X2[EJ(X)∂2W∂X2]+μ(X)∂W2∂t2+c∂W∂t=q(X,t); (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(x,y,t) , 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≤x≤a ,
0≤y≤b ,
0≤t≤t1 . Considering the region
0≤x≤a ,
0≤y≤b , at each time level, we obtain a two-dimensional problem, with
0≤t≤∞ . 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(x,y,t) , for an approximate description of its change in time at
x=const ,
y=const we apply the algebraic polynomial
W(t)=a0+a1t+a2t2+⋯+antn, (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(t)=a0+a1t+a2t2, Wt=dWdt=a1+2a2t; Wtt=d2Wdt2=2a2 (3)
In which
Wt is the speed, and
Wtt the acceleration. The coefficients
a0,a1,a2 are determined from the initial conditions.
Let suppose that for
t=ti
W(ti)=Wi , and for
t=ti−1
W(ti−1)=Wi−1 ,
Wt(ti−1)=Wti−1 , then one can express the coefficients
a0,a1,a2 as functions of
Wi,Wi−1,Wti−1 .
For
t=ti−1 from Equation (2) we get:
W(ti−1)=Wi−1=a0+a1ti−1+a2t2i−1, Wt(ti−1)=Wti−1=a1+2a2ti−1. (4)
For
t=ti from Equation (1) we obtain:
W(ti)=Wi=a0+a1ti+a2t2i. (5)
Considering that
ti=ti−1+τi , where
τi -is the time step, solving Equations (4) and (5), simultaneously we find:
{a2=1τ2i(Wi−Wi−1)−Wti−1τi;a1=−2ti−1τ2i(Wi−Wi−1)+Wti−1(1+2ti−1τi);a0=Wi−1(1−t2i−1τ2i)+Wit2i−1τ2i−Wti−1(ti−1+t2i−1τi). (6)
Substituting the found values of the coefficients
a0,a1,a2 in Equation (2), taking into consideration Equation (3), we will get at
t=ti
Wti=−Wti−1−2τi(Wi−1−Wi), i=1,2,3,⋯ (7)
Wtti=−2τiWti−1−2τ2i(Wi−1−Wi); i=1,2,3,⋯ (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
¨W+ω2W=0 at
W1=0 admits as a solution
W=˙W1sin(ωt) , where
˙W1=Wt1 is a given initial velocity. Suppose
ω=W=1 , then
¨W+W=0 , otherwise considering Equation (3):
Wtt+W=0 (9)
and
W=sin(t) ,
Wtt=−sin(t) . Let write Equations (7) and (8) while
τ=cont in another form:
Vi=−Vi−1+2(Wi−Wi−1), (10)
ai=2(Wi−Wi−1−Vi−1); (11)
where
Vi=Wtiτ ,
ai=Wttiτ2 . By rewriting Equation (9) at
t=ti taking into consideration Equation (10) and Equation (11) leads to:
Wi=22+τ2(Vi−1+Wi−1). (12)
By solving Equations (10)-(12) simultaneously when dividing a quarter of the period into 32 equal parts along the length eighth equal parts (
τ=π2n , where n-the number of partitions), with the initial conditions:
W1=Wtt1=0 ;
Wt1=V1τ=1 (international units system), we obtain the following results at
t=π2 :
W=0.98037 ,
Wt=0.03692 ,
Wtt=−0.98038 , while the exact values are
W=1 ,
Wt=0 ,
Wtt=−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] .
{∂2m∂ξ2=−(p−ˉμvˉt ˉt−ˉcvˉt)∂2m∂ξ2=−γm (13)
In which
vˉt=∂v∂ˉt ,
vˉt ˉt=∂2v∂ˉt2 . Where the parameters are linked to original variables as:
ξ=Zl ;
m=Mq0l2 ;
v=WEJ0q0l4 ;
p=q(Z,t)q0 ;
ˉt=tl2√EJ0μ0 ;
ˉμ=μ(z)μ0 ;
γ=EJ0EJ(Z) ;
ˉc=cl2√EJ0μ0 . (14)
μ(z) being the mass per unit length of the beam;
μ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;
{vˉtij=−vˉti−1j−2ˉτi(vi−1j−vij),vˉt ˉtij=−2ˉτivˉti−1j−2ˉτ2i(vi−1j−vij), (15)
with
i=1,2,3,⋯ which is measured along the axis
ˉt , while
j=1,2,3,⋯,(n−1) is measured along the axis ξ.
ˉτi=τil2√EJ0μ0 is the dimensionless time step;
v,vˉt,vˉt ˉ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
ˉt=ˉti in the absence of discontinuities
m,v,p , and derivatives
v,p , is obtained by replacing p by
p−ˉμvˉt ˉt−ˉcvˉt , leading to:
mij−1−2mij+mij+1+hΔmξij=−h212[(p−ˉμvˉt ˉt−ˉcvˉt)ij−1+10(p−ˉμvˉt ˉt−ˉcvˉt)ij+(p−ˉμvˉt ˉt−ˉcvˉt)ij+1]. (16)
The finite difference equation approximating the second line of Equation (13) following [6] can be written as:
vij−1−2vij+vij+1=−h212[γj−1mij−1+10γjmij+γj+1mij+1]−γjh312Δmξij (17)
Accounting the system the set of Equation (15) into Equation (16) for a constant time step (
ˉτi=ˉτ ), leads to the recurrent equation for a regular point j in ith time layer as follows:
6ˉτ2h2(mij−1−2mij+mij+1)+6ˉτ2h2Δmξij−(dj−1vij−1+10djvij+dj+1vij+1)=−(dj−1vi−1j−1+10djvi−1j+dj+1vi−1j+1)−ˉτ(ej−1vˉti−1j−1+10ejvˉti−1j+ej+1vˉti−1j+1) −ˉτ22(pij−1+10pij+pij+1), (i=1,2,3,⋯),(j=1,2,3,⋯) (18)
with
d=ˉμ+ˉτ+ˉc ;
e=ˉμ+ˉτ2ˉ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
∂v∂ξ , the bending moment m, and the shear force
∂m∂ξ . 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 ith time layer as:
vij=vij(0) ,
mij=mij(0) , (19)
where
vij(0) and
mij(0) are the deflection and bending moment at the endpoint, respectively, which generally are functions of
ˉ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 ith temporal layer as follows:
mij=mij(0) ,
∂m∂ξ=mξij(0) . (20)
mξij(0) 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 ith time layer is given by [15] :
(∂m∂ξ)ij=mξij=1h(−mij+mij+1)+h12[5(p−ˉμvˉt ˉt−ˉcvˉt)ij+(p−ˉμvˉt ˉt−ˉcvˉt)ij+1] +h212(p−ˉμvˉt ˉt−ˉcvˉt)ξij, (21)
in which by the superscript ξ, we mean the partial derivative with respect to variable ξ. The above equation can be discretized to give:
mξij=1h(−mij+mij+1)+h12(5pij+pij+1)+h212pξij−h12(5ˉμjvˉt ˉtij+ˉμj+1vˉt ˉtij+1) −h212ˉμξjvˉt ˉtij−h212ˉμj(vˉt ˉt)ξij−h12(5ˉcjvˉtij+ˉcj+1vˉtij+1)−h212ˉcj(vˉt)ξij. (22)
Let us express
(vˉt)ξij in (22) in terms of
vˉt in a square parabola:
(vˉt)ξij=12h(−3vˉtij+4vˉtij+1−vˉtij+2). (23)
Inserting
(vˉt)ξij , above into we get:
mξij=1h(−mij+mij+1)+h12(5pij+pij+1)+h212pξij −h12[(3.5ˉμj+hˉμξj)vˉt ˉtij+(2ˉμj+ˉμj+1)vˉt ˉtij+1−0.5ˉμjvˉt ˉtij+2] −h12[3.5ˉcjvˉtij+(2ˉcj+ˉcj+1)vˉtij+1−0.5ˉcjvˉtij+2]. (24)
Which can be rewritten taking into consideration the set of equations (15) as:
mξij=1h(−mij+mij+1)−h6ˉτ2[(3.5dj+hdξj)vij+(2dj+dj+1)vij+1−0.5djvij+1] +h6ˉτ2[(3.5dj+hdξj)vi−1j+(2dj+dj+1)vi−1j+1−0.5djvi−1j+1] +h6ˉτ[(3.5ej+heξj)vˉti−1j+(2ej+ej+1)vˉti−1j+1−0.5ejvˉti−1j+1] +h12(5pij+pij+1)+h212p,ξij where i=1,2,3,⋯ (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(0) instead of
mξij in Equation (25).
· Beam with left end rigidly fixed
The boundary conditions for nod j of the restrained edge of the beam in the ith time layer read:
vij=vij(0) ,
(∂v∂ξ)ij=vξij=vξij(0) , (26)
where
vξij(0) is the specified rotation angle at the edge point. The second member of Equation (26) can be discretize for nod j in the ith time layer as follows [15] :
vξij=1h(−vij+vij+1)+h12(5γj+hγξj)mij+h12γj+1mij+1+h212γjmξij. (27)
Substituting (25) into (27), we get:
vξij=h3γj72ˉτ2[(3.5dj+hdξj+72ˉτ2h4γj)vij+(2dj+dj+1+72ˉτ2h4γj)vij+1−0.5djvij+2] +h3γj72ˉτ2[(3.5dj+hdξj)vi−1j+(2dj+dj+1)vi−1j+1−0.5djvi−1j+2] +h3γj72ˉτ2[(3.5ej+heξj)vˉti−1j+(2ej+ej+1)vˉti−1j+1−0.5ejvˉti−1j+2] +h12[(4γj+hγξj)mij+(γj+γj+1)mij+1]+h3γj144[5pij+pij+1+hpξij], i=2,3,⋯ (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(0) instead of
vξij 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ξ 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:
vij=0.9vij+120(vij−1+vij+1)+h2240(γj−1mij−1+10γjmij+γj+1mij+1)−h3γj240Δmξij; (29)
mij=0.9mij+120(mij−1+mij+1+hΔmξij)−h2120ˉτ2(dj−1vij−1+10djvij+dj+1vij+1) +h2120ˉτ2(dj−1vi−1j−1+10djvi−1j+dj+1vi−1j+1) +h2120ˉτ(ejvˉti−1j−1+10ejvˉti−1j+ej+1vˉti−1j+1)+h2240(pij−1+10pij+pij+1). (30)
The deflection at the free end of the beam is determined by remembering Equation (25) as:
vij=(2dj+dj+1)(3.5dj+hdξj)vij+1+0.5dj(3.5dj+hdξj)vij+2+6ˉτ2(3.5dj+hdξj)(mij+1−mij) +vi−1j+(2dj+dj+1)(3.5dj+hdξj)vi−1j+1−0.5dj(3.5dj+hdξj)vi−1j+2
+ˉτ(3.5ej+heξj)(3.5dj+hdξj)vˉti−1j+ˉτ(2ej+ej+1)(3.5dj+hdξj)vˉti−1j+1−0.5ˉτej(3.5dj+hdξj)vˉti−1j+2 +0.5ˉτ2(3.5dj+hdξj)(5pij+pij+1+hpξij)−6ˉτ2h(3.5dj+hdξj)mξij(0). (31)
The bending moment at the embedded end of the beam is determined using Equation (28), leading to:
mij=−γj+γj+14γj+hγξjmij+1−h2γj12(4γj+hγξj)(5pij+pij+1+hpξij) +h2γj6ˉτ2(4γj+hγξj)[(3.5dj+hdξj+72ˉτ2h4γj)vij+(2dj+dj+1−72ˉτ2h4γj)vij+1−0.5djvij+2] −h2γj6ˉτ2(4γj+hγξj)[(3.5dj+hdξj)vi−1j+(2dj+dj+1)vi−1j+1−0.5djvi−1j+2] −h2γj6ˉτ2(4γj+hγξj)[(3.5ej+heξj)vˉti−1j+(2ej+dej+1)vˉti−1j+1−0.5ejvˉti−1j+2]+12h(4γj+hγξj)vξij(0). (32)
The resulting formulas correspond to a regular grid with steps
ˉτ (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
ˉτ , as well as on the values of and , which are variables. The sequence of the algorithm implementation is as follows. If
ˉt=0(i=1)
vij=vij(0) ;
vˉtij=vˉtij(0) .
For
ˉt=ˉτ(i=2) the values
v2j ,
m2j are calculated from the solution of the above set of equations. The variables
vˉt2j are determined from Equation (17) at
ˉτj=ˉτ .
For the calculus,
v3j ,
m3j , the
v2j ,
m2j ,
vˉt2j 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] :
∂W(Z,t)∂t|t=0=S(Z)μ(Z). (33)
where μ is the linear mass of the beam.
Let us suppose that in Equation (14)
q0=S0l2√EJ0μ0 , where
S0 is the instantaneous impulse uniformly distributed over the entire length of the beam.
Then from the second line of Equation (13) we have
ω=W√EJ0μ0S0l2 ;
m=MS0√μ0EJ0 . We thus write down Equation (33) for dimensionless unknowns on one has
∂v∂t=ˉSˉμ. (34)
where
ˉS=S(Z)S0 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
ˉc=0 ,
p=0 ,
vˉt1=ˉS/ˉμ ,
v1=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
ˉc=0 ,
p≠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:
ˉt=0 ;
v=0 ;
∂(vt)=|vˉt=ˉSˉμ=1 . To this end Equations (29)-(31) are numerically implemented for varying time step
ˉτi=ˉτ and space step h, and for the following parameters:
γj−1=γj=γj+1=1 ;
ˉc=0 ;
dj−1=dj=dj+1=1 ;
ej−1=ej=ej+1=1 ;
pξij=pij−1=pij=pij+1=0 . (35)
In Figure 1, the changes of deflection and bending moment are plotted in the middle of the span (that is for
ξ=0.5 ) without taking into consideration the time dependency of the damping. Curves 1, 2, 3, 4, 5 are for (
h=110 ,
ˉτ=1100 ), (
h=112 ,
ˉτ=1144 ), (
h=116 ,
ˉτ=1256 ), (
h=120 ,
ˉτ=1400 ), (
h=130 ,
ˉτ=0.002π ), respectively. As can be seen, the deflection (curve) are fairly similar for all h and
ˉτ . 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 (
ξ=0.5 ).
bending moment | Grid size h ( ) | , |
| | | | |
mcp | 1.259 | 1.291 | 1.365 | 1.427 | 1.572 |
10v | 1.250 | 1.259 | 1.262 | 1.264 | 1.279 |
翻译:
Table 1. Values of the bending moment and deflection coefficients in the middle of the beam (
ξ=0.5 ).
ˉc≠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
ˉc=0.05 and forvarying values of h and
ˉτ . 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
γj=γj+1=1 ,
γξ=0 ,
vij=0 ,
pij−1=pij=pij+1=pξij=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=130 ,
ˉτ=0.002π . 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 (
ξ=0.5 ).
bending moment | Grid size h ( ) | , |
| | | | |
mcp | 1.255 | 1.287 | 1.361 | 1.422 | 1.567 |
10v | 1.246 | 1.254 | 1.257 | 1.259 | 1.269 |
翻译:
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=130 ,
ˉτ=0.002π , corresponds approximately to that at
t=T14 , (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=T12 , and the maximum value of the bending moment in the middle of the span is 1.342 and occurs later than
t=T12 . For an articulated beam, the maxima of the deflection and bending moment in the middle of the span appear approximately at
t=T5 , which is in agreement with results found in [33] , andaccordingly, are equal to
v(max)=0.1279 ,
mcp(max)=1.572 (for
h=130 ,
ˉτ=0.002π ).
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:
pij=sin(1.6πˉti), (36)
h | | | | |
| | | | | | | |
| 0.976 | 1.000 | 1.005 | 1.121 | 1.194 | 1.297 | 1.342 |
| 0.539 | 0.547 | 0.548 | 0.565 | 0.571 | 0.576 | 0.581 |
| −1.655 | −1.750 | −1.742 | −2.022 | −2.156 | −2.256 | −2.402 |
翻译:
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
ˉti=tiT ,
T=2πωmin , and
ωmin is the lowest natural frequency of the hinged beam.
Based on Equations (17) and (18) taking into consideration Equations (16), (22), (36)
γ=d=e=ˉμ=1 ,
Δmξ=0 , the maximum values of mcp 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
ˉt=0.425 for
h=14 ,
ˉτ=116.
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
r1+r2=1.0736 , the integral sum of the inertial forces
∑J=02369 , and the resultant harmonic load taken along the entire length of the beam
P=0.8197 at
h=14 ,
ˉτ=116 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
| Successive Approximation Method (SAM) | Newmark |
| h2 | |
h | | | | | |
| 0.1723 | 0.1948 | 0.2068 | 0.2117 | 0.2023 |
| 0.01768 | 0.02002 | 0.02124 | 0.02174 | 0.02085 |
翻译:
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 (
ξ=0.5 ).
h=14 ,
ˉτ=116 .
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 ith time layer:
{vˉtij=vˉti−1j+ˉτ(1−δ)vˉt ˉti−1j+δvˉt ˉtijvˉt ˉtij=(vij−vi−1j)αˉτ2−vˉti−1jαˉτ+(1−0.5α)vˉt ˉti−1j+δvˉt ˉtij (37)
where
δ≥0.50 ;
α≥0.25(0.5+δ)2 . If in the approximation of Equation (13) we use Equations (32) and (33), respectively, then for the case
ˉc=0 we get:
mij−1−2mij+mij+1=h2(ˉμj−1vij−1+10ˉμjvij+ˉμj+1vij+1)12αˉτ2−h2(ˉμj−1vi−1j−1+10ˉμjvi−1j+ˉμj+1vi−1j+1)12αˉτ2 −h2(ˉμj−1vˉti−1j−1+10ˉμjvˉti−1j+ˉμj+1vˉti−1j+1)12αˉτ +h2(ˉμj−1vˉt ˉti−1j+10ˉμjvˉt ˉti−1j+ˉμj−1vˉt ˉti−1j+1)12−h2(pij−1+10pij+pij+1)12 (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
δ=12 ,
α=14 ,
h=ˉτ=110 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 mcp and v in the middle of the span and m1,2 in the embedded ends of the beam are shown in Table 5. These values are obtained for
h=14 ,
h=16 ,
h=18 ,
h=110 and for
ˉτ=h2 . The last column of this table shows the results obtained by the analytical method at
h=16 ,
ˉτ=136 , 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 mcp, v, and m1,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.
| Successive Approximation Method (SAM) | Analytical method [29] |
| h2 | |
h | | | | | |
| 0.045190 | 0.045890 | 0.046890 | 0.047730 | 0.044100 |
| 0.002851 | 0.002853 | 0.002900 | 0.002946 | 0.002740 |
| −0.083290 | −0.087210 | −0.089550 | −0.09137 | −0.087300 |
翻译:
Table 5. Values of the bending moment and deflection coefficients in the middle and at the ends of the beam.
| | | | |
h | | | | |
mcp | 1.237 | 1.290 | 1.309 | 1.331 |
v | 0.1123 | 0.1162 | 0.1175 | 0.1176 |
翻译:
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(Z)=EJ0(1+0.1Zl)3, μ(Z)=μ0(1+0.1Zl), (39)
where
EJ0 and
μ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 mcp, v, in the middle of the span of an articulated beam. The maximum dimensionless bending moment equal to
mcp=1.331 at
h=112 ,
ˉτ=1144 occurs at around the time
ˉ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=112 ,
ˉτ=1144 , 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 mcp 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=112 ,
ˉτ=1144 .
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 (
ξ=0.5 ).
1 | 2 | 3 | 4 | 5 | 6 |
| h | mcp | v | m1 | m2 |
| | 0.9685 | 0.04759 | −1.355 | −1.598 |
| | 1.039 | 0.05122 | −1.542 | −1.733 |
| | 1.058 | 0.05225 | −1.685 | −1.882 |
| | 1.069 | 0.05290 | −1.759 | −1.953 |
翻译:
Table 7. Values of the bending moment and deflection coefficients in the middle and at the ends of the beam.
(
h=14 ,
ˉτ=116 ); (
h=16 ,
ˉτ=136 );
(h=18 ,
ˉτ=164 ); (
h=110 ,
ˉτ=1100 ); (
h=112 ,
ˉτ=1144 ).
mcp=McpS0√μ0EJ0 ,
v=WS0l2√EJμ0 ,
m1,2=M1,2S0√μ0EJ0 . (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
ξ=0 and
ξ=1 .
Comparison of the results shown in both Table 6 and Table 7 at
h=112 ,
ˉτ=1144 , it is obvious that the bending moment in the middle of the span (
mcp=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=112 ,
ˉτ=1144 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 mcp 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
mcp=0.1973 for
h=110 , and
ˉτ=1100 , which occurs approximately at
ˉ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 mcp 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=110 ,
ˉτ=1100 .
Let’s check the integral condition of the beam equilibrium. For this in a time
ˉt=0.425 , the dimensionless support reactions
rA+rB=1.055 , the integral sum of the inertial forces
∑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=110 ,
ˉτ=1100 . The results show that
∑J+P=1.0925 differs little from
| | | | |
h | | | | |
mcp | 0.17080 | 0.19100 | 0.19710 | 0.19730 |
v | 0.01560 | 0.01760 | 0.01830 | 0.01838 |
翻译:
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 (
ξ=0.5 ).
rA+rB=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=14 ,
ˉτ=116 ;
h=16 ,
ˉτ=136 ;
h=18 ,
ˉτ=164 ;
h=110 ,
ˉτ=1100 . 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
ξ=0 and
ξ=1 . Figure 9 shows the curves
1 | 2 | 3 | 4 | 5 | 6 |
| h | mcp | v | m1 | m2 |
| | 0.04658 | 0.002613 | −0.08083 | −0.09162 |
| | 0.04660 | 0.002629 | −0.08409 | −0.09486 |
| | 0.04685 | 0.002633 | −0.08567 | −0.09651 |
| | 0.04690 | 0.002641 | −0.08642 | −0.09728 |
翻译:
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 mcp, v in the middle of the beam span and
m1,m2 at its edges, respectively at
ξ=0 ,
ξ=1 . The oscillation process of the beam was studied over some time
t=0.5T . Curve 1 (mcp), curve 2 (v), curve 3 (m1) and curve 4 (m2) were plotted from the results obtained at
h=110 ,
ˉτ=1100 .
The integral condition of the beam equilibrium was checked. At the moment in time
ˉt=0.25 , (
mcp=0.04773 ,
v=0.002946 ,
m1=−0.08642 ,
m2=−0.09728 ) dimensionless support reactions
rA,rB have been calculated. The sum of those reactions gives
rA+rB=0.9648 and the integral sum of inertial forces
∑J=0.0110 , while the resultant harmonic force was determined, taken over the entire length of the beam
P=0.9150 (
h=110 ,
ˉτ=1100 ). The results show that
∑J+P=0.9260 is somewhat different from
rA+rB=
0.9648. The solution error is 4%.
Comparison of the results obtained here for
h=110 ,
ˉτ=1100 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.