Influence of control parameters on the stabilization of an Euler-Bernoulli flexible beam

André-Pascal Abro Goh1, Kidjégbo Augustin Touré2
and Gossrin Jean-Marc Bomisso3
(Date: December 27, 2023; accepted: June 07, 2024; published online: July 11, 2024.)
Abstract.

In this work, we numerically study the influence of control parameters on the stabilization of a flexible Euler-Bernoulli beam fixed at one end and subjected at the other end to a force control and a moment control proportional respectively to velocity and rotating velocity. First, we analyze the displacement stabilization and the asymptotic behavior of the beam energy using a stable numerical scheme, resulting from the Crank-Nicholson algorithm for time discretization and the finite element method based on the approximation by Hermite’s cubic polynomial functions, for discretization in space. Then, by means of the finite element method, we represent the spectrum of the operator associated with this beam problem and we carry out a qualitative study of the locus of the eigenvalues according to the positive control parameters. From these studies we conclude that rotating velocity control has more effect on the stabilization of the beam compared to velocity control. Finally, this result is confirmed by a sensitivity study on the control parameters involved in the stabilization of the beam.

Key words and phrases:
beam equation, finite element method, asymptotic stability.
2005 Mathematics Subject Classification:
35B35, 65N30, 35B40
1Département de Science et Technique, Université Alassane Ouattara, Bouaké, Côte d’Ivoire, e-mail: esaie54v1@gmail.com.
2Départment de Mathematiques et Informatique, Institut National Polytechnique Houphouet Boigny, Yamoussoukro, Côte d’Ivoire, e-mail: latoureci@gmail.com.
3UFR Sciences Fondamentales et Appliquées, Université Nangui Abrogoua D’Abobo-Adjamé, Abidjan, Côte d’Ivoire, e-mail: bogojm@yahoo.fr.

1. Introduction

Boundary feedback stabilization of systems modeled by an Euler-Bernoulli beam equation is an important area of research for engineers and mathematicians alike. Many mechanical systems, such as telecommunications antennas and flexible robot arms, can be modeled using Euler-Bernoulli beam equations. In order to control the asymptotic dynamic behavior of these systems, a control law is associated with them. The aim is to promote energy dissipation in order to stabilize the system within a reasonable time.

In our study, we consider a flexible Euler-Bernoulli beam fixed at one end and subjected, at the other end, to a force control in velocity and a moment control in rotational velocity. The beam satisfies the Euler-Bernoulli equation

2y(x,t)t2+4y(x,t)x4=0,  0<x<1,t>0,

with initial conditions

y(x,0)=y0(x),y(x,0)t=v0(x),  0<x<1.

At the embedded end (x=0), we have the following boundary conditions

y(0,t)=y(0,t)x=0,t>0.

At the end x=1, we apply a force control in velocity and a moment control in rotating velocity

2y(1,t)x2=β2y(1,t)xt,3y(1,t)x3 =αy(1,t)t,t>0.

Notice that these last two equations indicate that the bending moment (i.e. 2y(1,t)x2) and lateral force (i.e. 3y(1,t)x3) of the beam at the free end (x=1), are respectively controlled by the feedback laws in rotating velocity (β2y(1,t)xt) and velocity (αy(1,t)t). Also, we assume that α, β are two positive constants involved in the boundary controls.

Assembling this together, we obtain the mathematical problem

(1) 2y(x,t)t2+4y(x,t)x4 =0, 0<x<1,t>0,
(2) y(0,t)=y(0,t)x =0, t>0,
(3) 2y(1,t)x2=β2y(1,t)xt,3y(1,t)x3 =αy(1,t)t, t>0,
(4) y(x,0)=y0(x),y(x,0)t =v0(x), 0<x<1,

where y(x,t) is the transverse deflection of the beam at position x and time t. Also, the bending stiffness function, mass density function and beam length are taken to be equal to unity.

In [2], the system (1)–(4) is formulated as an evolution problem and it is proven that this problem is well posed in the sense of C0-semigroups of contractions. Also, Mensah E. Patrice in [7], studies the spectrum of the differential operator appearing in the exponential stabilization of this system. By means of a numerical scheme obtained by the finite difference method, the author studies the locus of eigenvalues as a function of the positive feedback parameters α and β, carrying out a qualitative study of the position of the spectrum with respect to a vertical asymptote.

In the case of our study, the aim is to determine the influence of each control parameter on the stabilization of this beam system and, more importantly, to show which one has the greatest influence. This will enable us to calibrate these parameters appropriately, so as to stabilize the system more rapidly.

We begin by developing a numerical scheme equivalent to the system (1)–(4), using the Crank-Nicholson algorithm for discretization in time and the finite element method based on approximation by cubic Hermite functions for discretization in space (see e.g. [4], [8]). Using this numerical scheme, we graphically analyze the influence of the positive feedback parameters α and β on displacement stabilization and the asymptotic behavior of beam energy. Then, using the finite element method, we represent the spectrum of the differential operator associated with the system (1)–(4) in order to qualitatively describe the impact of the control parameters on the positioning of the spectrum and its asymptote. In this way, we can observe the effect of each parameter on the energy decay rate and stabilization of the system under study. Finally, the conclusions of these graphical analyses are confirmed by a sensitivity study on these control parameters.

2. Numerical approximation of the system (1)–(4).

Let v(.,t)=y(.,t)t and consider the space

={w=(y,v):yHl2(0,1),vL2(0,1)},

with

Hl2(0,1)={yH2(0,1):y(0)=y(0)x=0},

where in superscript represents the transposition.

The system (1)–(4) can be written as the evolution problem

(5) {ddtw(t)=𝒜w(t),w(0)=w0.

With w(t)=(y(.,t),y(.,t)t), w(0)=(y0,v0) for all t>0. Note that 𝒜 is a linear operator, with domain

(6) D(𝒜)=
={(f,g)(Hl2(0,1)H4(0,1))×Hl2(0,1):2f(1)x2=βg(1)x,3f(1)x3=αg(1)}

and defined by

(7) 𝒜(fg)=(g4fx4).

2.1. Discretization

2.1.1. Time discretization.

Let T be a strictly positive real. The interval [0,T] is discretized into Nt intervals of the same length Δt. Then, at each moment tn=nΔt, we note respectively yn and vn the approximate values of y(tn,x) and v(tn,x). Thus, we have wn=w(tn)=(yn,vn) with 0nNt.
Using the Crank-Nicholson scheme for the time discretization of the problem (5), we obtain the following system

(8) {yn+1ynΔt=vn+1+vn2,vn+1vnΔt+12(4x4(yn+1+yn))=0,withy0=y0andv0=v0.

Moreover, for all ΨL2(0,1) and ϕHl2(0,1), the variational formulation of the problem (8), is written

(9) {01yn+1ynΔtΨ(x)𝑑x=01vn+1+vn2Ψ(x)𝑑x,01vn+1vnΔtϕ(x)𝑑x+1201((2yn+1x2+2ynx2)2ϕ(x)x2)𝑑x+α2(vn+1(1)+vn(1))ϕ(1)+β2(vn+1(1)x+vn(1)x)ϕ(1)x=0.

2.1.2. Discretization in space.

For discretization in space, we use the finite element method with Hermite’s cubic polynomial functions.
First, let us subdivide the interval [0,1] into Ne intervals [xi,xi+1] of the same length h=xi+1xi with xi=ih and Neh=1. On each node xi, we associate the approximate values of y and its derivative at the point xi.

Thus, in the same way as in [5], consider the functions (ϕ)1k4 as reference functions, from which the polynomial functions, φ1i(x) and φ2i(x), are constructed such that:

For i=1,,Ne1,

φ1i(x)={ϕ3i1(x),ifx[xi1,xi],ϕ1i(x),ifx[xi,xi+1],0,otherwise.
φ2i(x)={ϕ4i1(x),ifx[xi1,xi],ϕ2i(x),ifx[xi,xi+1],0,otherwise.

For i=Ne,

φ1Ne(x)={ϕ3Ne1(x)ifx[xNe1,1],0,otherwise.
φ2Ne(x)={ϕ4Ne1(x),ifx[xNe1,1],0,otherwise.

with φ1i(0)=φ2i(0)=φ1ix(0)=φ2ix(0)=0.

We thus obtain a 2Ne-tuple (φ11,φ21,φ12,φ22,,φ1Ne,φ2Ne) which we simply denote in the following (Φi)1i2Ne and which constitutes a basis of dimension 2Ne of the interpolation space Vh. With

Vh={u𝒞([0,1]):u[xi,xi+1]P3,1iNe,u(0)=ux(0)=0}.

Thus, we have

(10) yn=i=12NeyinΦi(x),vn=i=12NevinΦi(x)

and

(11) Yn =[y1n,y2n,,y2Ne1n,y2Nen],
(12) Vn =[v1n,v2n,,v2Ne1n,v2Nen],
(13) Wn =[YnVn],

with respectively yin and vin, the approximate values of yn(xi) and vn(xi).

By replacing the expressions (10) in (9), we obtain, for all i,j=1,,2Ne, the following numerical scheme

(14) {i=12Ne(2yin+1Δtvin+1)01Φi(x)Φj(x)𝑑x=i=12Ne(2yin+Δtvin)01Φi(x)Φj(x)𝑑xi=12Ne(2vin+101Φi(x)Φj(x)dx+Δtyin+1012ϕi(x)x22ϕj(x)x2dx+Δtαvin+1Φi(1)Φj(1)+Δtβvin+1ϕi(1)xϕj(1)x)==i=12Ne(2vin01Φi(x)Φj(x)dxΔtyin012ϕi(x)x22ϕj(x)x2dxΔtαvinΦi(1)Φj(1)Δtβvinϕi(1)xϕj(1)x).

2.2. Stability of the numerical scheme

For all, fixed point xj with 1j2Ne, let us set

(15) a=01(Φj(x))2𝑑x,b=01(2ϕj(x)x2)2𝑑x,c=(Φj(1))2,d=(ϕj(1)x)2.

Then, the system (14) can be written as the following vector equation:

(16) [1Δt212aΔt+αc+βdb][Yjn+1Vjn+1]=[1Δt212aΔt+αc+βdb][YjnVjn].

By multiplying (16) by the inverse of the left matrix, we obtain:

(17) [Yjn+1Vjn+1]=12aΔt+αc+βdb+Δt2[2aΔt+αc+βdb+Δt2Δt(αc+βd)b02aΔt+αc+βdbΔt2][YjnVjn].

Now, to determine the associated Von-Neumann stability condition, let us consider

(18) [YjnVjn]=[eikxjY^jneikxjV^jn].

The system (17) becomes

(19) [Y^jn+1V^jn+1]=[12Δt(αc+βd)+Δtb4aΔt+2(αc+βd)0Δtb+4aΔt+2(αc+βd)Δtb+4aΔt+2(αc+βd)][Y^jnV^jn].

The eigenvalues associated with the amplification matrix (19) are λ1=1 and λ2=Δtb+4aΔt+2(αc+βd)Δtb+4aΔt+2(αc+βd). It is clear that these eigenvalues satisfy the Von-Neumann conditions. Indeed, we have λi1,i=1,2. We deduce from this that the numerical scheme (14) is stable.

3. Effects of control parameters on beam stabilization

3.1. Stabilization in motion

For all i,j=1,,2Ne, consider the matrices

𝕄i,j =01ΦiΦj𝑑x,
𝕊i,j =αΦi(1)Φj(1)+β(ϕi)x(1)(Φj)x(1),
𝕂i,j =01(Φi)xx(Φj)xx𝑑x.

The system (14) can then be written as follows:

(20) PWn+1 =QWn,
W0 =W0

with

P=[2I2NeΔtI2NeΔtK2M+ΔtS],
Q=[2I2NeΔtI2NeΔtK2MΔtS].

I2Ne denotes the unit matrix of order 2Ne.
For numerical simulations, we consider as initial conditions:

y0(x) =3x2+x3,
v0 =0.

with

T =10,
Ne =70,
Δt =0,01.

We represent the position and the angle of rotation of the beam at the free end (i.e. in x=1), for different values (taken arbitrarily) of the parameters α and β, in order to highlight their influences on the stabilization of the beam. The numerical simulations are carried out with MATLAB.

case 1: We set α=0.3 and we vary β. case 2: We set β=0.3 and we vary α.

Refer to caption
Figure 1. Position for different values of β.
Refer to caption
Figure 2. Rotation angle for different values of β.
Refer to caption
Figure 3. Position for different values of α.
Refer to caption
Figure 4. Rotation angle for different values of α.

3.2. Energy stabilization

Multiply the equation (1) by y(x,t)tHl2(0,1). After two integrations by parts, we obtain

012y(x,t)t2y(x,t)t𝑑x+012y(x,t)x22x2(y(x,t)t)𝑑x+α(y(1,t)t)2+β(2y(1,t)xt)2=0.

So,

(21) dE(t)dt+α(y(1,t)t)2+β(2y(1,t)xt)2=0.

Or,

(22) E(t)=1201(y(x,t)t)2𝑑x+1201(2y(x,t)x2)2𝑑x.

We deduce from the relations (21) and (22) that

(23) dE(t)dt0.

Thus, the energy E(t) defines a Lyapunov function and satisfies 0E(t)E(0)<.

Considering the Crank-Nicholson scheme, we have for 0nNt, the following discrete energy:

(24) En=1201(vn)2𝑑x+1201(2ynx2)2𝑑x.

We then obtain the following result:

Proposition 1.

For any integer n, there is a positive real ηn such that

(25) En+1EnΔt+ηn=0.
Proof.

According to (24), we have

(26) En+1En =1201((vn+1)2(vn)2)𝑑x+1201((2yn+1x2)2(2ynx2)2)𝑑x.

For x[0;1], we replace in the first equation of the variational formulation (9), Ψ(x) by (vn+1vn). We obtain the following equation

(27) 01yn+1ynΔt(vn+1vn)𝑑x=01(vn+1)2(vn)22𝑑x.

Furthermore, considering ϕ=yn+1 as a test function in the second equation of the problem (9), we obtain the following equation:

(28) 01 vn+1vnΔtyn+1dx+1201((2yn+1x2+2ynx2)2yn+1x2)𝑑x
+α2(vn+1(1)+vn(1))yn+1(1)+β2(vn+1(1)x+vn(1)x)yn+1(1)x=0.

Analogously, if we choose as test function ϕ=yn in the second equation of the system (9), we obtain:

(29) 01 vn+1vnΔtyndx+1201((2yn+1x2+2ynx2)2ynx2)𝑑x
+α2(vn+1(1)+vn(1))yn(1)+β2(vn+1(1)x+vn(1)x)yn(1)x=0.

Now making the difference between (28) and (29) and considering the relation

(30) yn+1yn=Δt2(vn+1+vn)

We obtain

(31) En+1EnΔt+ηn=0,for alln.

With

(32) ηn=α(vn+1(1)+vn(1)2)2+β(x(vn+1(1)+vn(1)2))2.

Therefore we have the following result:

Corollary 2.

For all n1, we have

En+1EnE0

We notice that the energy of the discrete system decreases and is controlled by that of the initial data. From Proposition 1, we deduce the following equality:

(33) En+1=EnΔt4[α(vn+1(1)+vn(1))2+β(x(vn+1(1)+vn(1)))2],

0nNt, with

(34) E0=1201(v0)2𝑑x+1201(2y0x2)2𝑑x.

Now let us represent the energy of the system (1)–(4) for different values of α and β.

case 1: We set α=0.08 and we vary β.

Refer to caption
Figure 5. Energy dissipation for different values of β.

case 2: We set β=0.08 and we vary α.

Refer to caption
Figure 6. Energy dissipation for different values of α.

3.3. Discussion of results

For a fixed value of α and different values of β, we see in Figs. 2 and 2 that vibrations at the free end of the beam cancel out at different times. In particular, they cancel faster as we increase the value of the β parameter. However, for a fixed value of β and for different values of α, we see in Figs. 4 and 4 that vibrations at the free end of the beam cancel out at almost the same time. Furthermore, in Figs. 5 and 6, we see that when we increase the value of parameter β, the energy decreases and cancels out faster than when the value of parameter α is increased.

The moment control in rotation velocity therefore seems to have more impact on the stabilization of the beam and energy dissipation.

4. Influence of control parameters on the system spectrum.

To begin, let us recall the following theoretical result.

Theorem 3 ([7]).

The eigenvalues of the unbounded operator 𝒜 which governs the system (1)–(4) take asymptotically the form

(35) λn=((1β+α)+𝒪(1n2))+i(π2(n14)2)2αβ1(αβ)22πβ2n+𝒪(1n2),n

and

limn+𝔢(λn)=(1β+α).

Notice that (1β+α) is the vertical asymptote of the spectrum in the complex plane.

Recall that for the numerical study of the spectrum, we use the finite element method with the cubic polynomial functions of Hermite defined previously.

When we multiply the equation (1) by the test function ϕHl2(0,1), we obtain after two integrations by parts the following weak formulation:

(36) 012y(x,t)t2ϕ(x)𝑑x+012y(x,t)x22ϕ(x)x2+αy(1,t)tϕ(1)+β2y(1,t)xtϕ(1)x=0.

Also, by separation of variables, the approximate solution in the base (Φi)1i2Ne, can be written as follows:

(37) y(x,t)=i=12Neyi(t)Φi(x).

We denote Y, a vector representation of the function y(x,t) defined as follows:

Y=[y1y2y2Ne].

Thus, the equation (36) becomes the following second order ordinary differential equation:

(38) 𝕄Ytt+𝕊Yt+𝕂Y=0.

Let V=dYdt and W=[YV]. Then the equation (38) becomes dW(t)dt=LW(t) with

L=[O𝕀𝕄1𝕂𝕄1𝕊]

the spectrum matrix.

Now, by implementing the spectrum matrix L, we represent the spectrum of the system (1)–(4) for different values of α and β. Note that for greater numerical precision, eigenvalues are calculated using the Advanpix multiprecision toolbox [3], quadruple precision (34 decimal digits) computation.

case 1: We set α=0 and we vary β.

Refer to caption
Figure 7. Effect of the β parameter on the spectrum.
Refer to caption
Figure 8. Effect of the parameter β on the asymptote of the spectrum.

case 2: We set α=0.08 and we vary β.

Refer to caption
Figure 9. Effect of the β parameter on the spectrum.
Refer to caption
Figure 10. Effect of the parameter β on the asymptote of the spectrum.

case 3: We set β=0.08 and we vary α.

Refer to caption
Figure 11. Effect of the α parameter on the spectrum.
Refer to caption
Figure 12. Effect of the parameter α on the asymptote of the spectrum.

4.1. Discussion of results

Numerical simulations of the spectrum have led to the following observations.

According to Figs. 10 and 10, for fixed values of α, when we increase the β value, the spectrum and asymptote shift significantly from left to right in the left complex half-plane. In addition, the number of low-frequency eigenvalues (those to the right of the asymptote ) decreases. It can therefore be said that β significantly improves the energy decay rate and stability of the system.

However, Figs. 12 and 12 show that for a fixed value of β, increasing the values of α leads to a slight shift in the spectrum and asymptote of the left complex half-plane. In addition, the number of eigenvalues at low frequencies remains the same. This means that the α parameter has little impact on the optimal energy decay rate and system stability.

We can therefore deduce that the moment control of the rotational velocity β has a greater influence on the optimal rate of energy decay and the stability of the system

4.2. Sensitivity of Parameters α and β on the Stability

In this part, we objectively study the influence of the parameters α and β on the stabilization of the beam. So, for different values of α and β, we note in the table below, the damping time δτ of the beam vibrations (i.e. time required to return to steady state y0).

The first-order Sobol sensitivity [6] is written

(39) Si=𝕍ar[𝔼(δτ|Yi)]𝕍ar(δτ),

where 𝕍ar(δτ) and 𝔼(δτ|Yi)) represent respectively the variance of δτ and the conditional expectation of δτ obtained by setting Yi=α and Yi=β in (39).

α β     0.8     1     5     10     15     20     30
0.1 5.53 6.69 17.07 24.76 26.6 27.62 31.41
0.5 3.67 4.61 8.39 9.17 9.51 9.93 9.86
0.85 4.28 4.18 10.81 15.17 15.53 15.19 12.08
1 4.61 4.91 11.79 16.85 18.22 18.46 18.08
1.5 5.70 5.6 14.39 21.85 24.51 26.77 27.48
2 6.69 7.49 19.44 25.62 28.86 30.88 37.79
Table 1. Beam vibration damping time (in seconds) for different values of α and β.

The calculation of the first order Sobol sensitivity indices gives us the following result:

Sα 0.36
Sβ 0.64.

According to the previous result, we have Sα<Sβ. We can deduce that the control parameter β, has more impact on the stabilization of the beam compared to the parameter α.

5. Conclusion

In this work, we were interested in the influence of control parameters on the stabilization of a flexible Euler-Bernoulli beam, fixed at one end and subjected at the other end to a punctual force control and moment control proportional respectively to the velocity and the rotation velocity. We have developed a stable numerical scheme by means of which, we have shown by graphical analysis, that the moment control in rotation velocity β has more influence on the stabilization of our model in displacement, in energy and on the positioning of the spectrum, in relation to the force control in velocity α. This result was then confirmed by a statistical study of the sensitivity of the control parameters.

References