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

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.

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

with initial conditions

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

$y(0,t)=\frac{\partial y(0,t)}{\partial x}=0,t>0.$

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

$\frac{{\partial }^{2}y(1,t)}{\partial x^2}=-\beta \frac{{\partial }^{2}y(1,t)}{\partial x\partial t},\frac{{\partial }^{3}y(1,t)}{\partial x^3}$ $=\alpha \frac{\partial y(1,t)}{\partial t},t>0.$

Notice that these last two equations indicate that the bending moment (i.e. $\frac{{\partial }^{2}y(1,t)}{\partial x^2}$) and lateral force (i.e. $\frac{{\partial }^{3}y(1,t)}{\partial x^3}$) of the beam at the free end ($x=1$), are respectively controlled by the feedback laws in rotating velocity $(-\beta \frac{{\partial }^{2}y(1,t)}{\partial x\partial t})$ and velocity $(\alpha \frac{\partial y(1,t)}{\partial t})$. Also, we assume that $\alpha$, $\beta$ are two positive constants involved in the boundary controls.

Assembling this together, we obtain the mathematical problem

(1) $\frac{{\partial }^{2}y(x,t)}{\partial t^2}+\frac{{\partial }^{4}y(x,t)}{\partial x^4}$ $=0,$

(2) $y(0,t)=\frac{\partial y(0,t)}{\partial x}$ $=0,$ $t>0,$

(3) $\frac{{\partial }^{2}y(1,t)}{\partial x^2}=-\beta \frac{{\partial }^{2}y(1,t)}{\partial x\partial t},\frac{{\partial }^{3}y(1,t)}{\partial x^3}$ $=\alpha \frac{\partial y(1,t)}{\partial t},$ $t>0,$

(4) $y(x,0)=y_0(x),\frac{\partial y(x,0)}{\partial t}$ $=v_0(x),$

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 C₀-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 $\alpha$ and $\beta$, 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 $\alpha$ and $\beta$ 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)=\frac{\partial y(.,t)}{\partial t}$ and consider the space

$$\mathscr{H}=\{w={(y,v)}^{\top }:y\in H_l^2(0,1),v\in L^2(0,1)\},$$

with

$$H_l^2(0,1)=\{y\in H^2(0,1):y(0)=\frac{\partial y(0)}{\partial x}=0\},$$

where $\top$ in superscript represents the transposition.

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

(5) $$\{\begin{array}{cc}\frac{d}{dt}w(t)=𝒜w(t),\hfill & \\ w(0)=w_0\in \mathscr{H}.\hfill & \end{array}$$

With $w(t)={(y(.,t),\frac{\partial y(.,t)}{\partial t})}^{\top }$, $w(0)={(y_0,v_0)}^{\top }$ for all $t>0$. Note that $𝒜$ is a linear operator, with domain

(6) $D(𝒜)=$

$=\{(f,g)\in \left(H_l^2(0,1)\cap H^4(0,1)\right)\times H_l^2(0,1):\frac{{\partial }^{2}f(1)}{\partial x^2}=-\beta \frac{\partial g(1)}{\partial x},\frac{{\partial }^{3}f(1)}{\partial x^3}=\alpha g(1)\}$

and defined by

(7) $$𝒜\left(\begin{array}{c}f\\ g\end{array}\right)=\left(\begin{array}{c}g\\ -\frac{{\partial }^{4}f}{\partial x^4}\end{array}\right).$$

2.1. Discretization

2.1.1. Time discretization.

Let $T$ be a strictly positive real. The interval $[0,T]$ is discretized into $N_t$ intervals of the same length $\mathrm{\Delta }t$. Then, at each moment $t_n=n\mathrm{\Delta }t$, we note respectively $y^n$ and $v^n$ the approximate values of $y(t_n,x)$ and $v(t_n,x).$ Thus, we have $w^n=w(t_n)=(y^n,v^n)$ with $0\le n\le N_t$.
Using the Crank-Nicholson scheme for the time discretization of the problem (5), we obtain the following system

(8) $\{\begin{array}{cc}\frac{y^{n+1}-y^n}{\mathrm{\Delta }t}=\frac{v^{n+1}+v^n}{2},\hfill & \\ \frac{v^{n+1}-v^n}{\mathrm{\Delta }t}+\frac{1}{2}\left(\frac{{\partial }^4}{\partial x^4}\left(y^{n+1}+y^n\right)\right)=0,\hfill & \\ \text{with}{y}^0=y_0\text{and}{v}^0=v_0.\hfill & \end{array}$

Moreover, for all $\mathrm{\Psi }\in L^2(0,1)$ and $\varphi \in H_l^2(0,1),$ the variational formulation of the problem (8), is written

(9) $\{\begin{array}{cc}{\displaystyle {\int }_0^1}\frac{y^{n+1}-y^n}{\mathrm{\Delta }t}\mathrm{\Psi }(x)𝑑x={\displaystyle {\int }_0^1}\frac{v^{n+1}+v^n}{2}\mathrm{\Psi }(x)𝑑x,\hfill & \\ {\displaystyle {\int }_0^1}\frac{v^{n+1}-v^n}{\mathrm{\Delta }t}\varphi (x)𝑑x+\frac{1}{2}{\displaystyle {\int }_0^1}\left(\left(\frac{{\partial }^2{y}^{n+1}}{\partial x^2}+\frac{{\partial }^2{y}^n}{\partial x^2}\right)\frac{{\partial }^2\varphi (x)}{\partial x^2}\right)𝑑x\hfill & \\ +\frac{\alpha }{2}\left(v^{n+1}(1)+v^n(1)\right)\varphi (1)+\frac{\beta }{2}\left(\frac{\partial v^{n+1}(1)}{\partial x}+\frac{\partial v^n(1)}{\partial x}\right)\frac{\partial \varphi (1)}{\partial x}=0.\hfill & \end{array}$

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 $N_e$ intervals $[x_i,x_{i+1}]$ of the same length $h=x_{i+1}-x_i$ with $x_i=ih$ and $N_{e}h=1$. On each node $x_i$, we associate the approximate values of $y$ and its derivative at the point $x_i$.

Thus, in the same way as in [5], consider the functions ${\left(\varphi \right)}_{1\le k\le 4}$ as reference functions, from which the polynomial functions, ${\phi }_1^i(x)$ and ${\phi }_2^i(x)$, are constructed such that:

For $i=1,\mathrm{\dots },N_e-1,$

For $i=N_e,$

with ${\phi }_1^i(0)={\phi }_2^i(0)=\frac{\partial {\phi }_1^i}{\partial x}(0)=\frac{\partial {\phi }_2^i}{\partial x}(0)=0.$

We thus obtain a $2N_e$-tuple $({\phi }_1^1,{\phi }_2^1,{\phi }_1^2,{\phi }_2^2,\mathrm{\dots },{\phi }_1^{N_e},{\phi }_2^{N_e})$ which we simply denote in the following ${\left({\mathrm{\Phi }}_i\right)}_{1\le i\le 2N_e}$ and which constitutes a basis of dimension $2N_e$ of the interpolation space $V_h$. With

$$V_h=\{u\in 𝒞\left([0,1]\right):u_{\mid [x_i,x_{i+1}]}\in P_3,1\le i\le N_e,u(0)=\frac{\partial u}{\partial x}(0)=0\}.$$

Thus, we have

(10) $y^n={\displaystyle \sum _{i=1}^{2N_e}}{y}_i^n{\mathrm{\Phi }}_i(x),v^n={\displaystyle \sum _{i=1}^{2N_e}}{v}_i^n{\mathrm{\Phi }}_i(x)$

and

(11) $Y^n$ $=[y_1^n,y_2^n,\mathrm{\dots },y_{2N_e-1}^n,y_{2N_e}^n],$

(12) $V^n$ $=[v_1^n,v_2^n,\mathrm{\dots },v_{2N_e-1}^n,v_{2N_e}^n],$

(13) $W^n$ $={\left[Y^n{V}^n\right]}^{\top },$

with respectively $y_i^n$ and $v_i^n$, the approximate values of $y^n(x_i)$ and $v^n(x_i).$

By replacing the expressions (10) in (9), we obtain, for all $i,j=1,\mathrm{\dots },2N_e$, the following numerical scheme

(14) $\{\begin{array}{cc}{\displaystyle \sum _{i=1}^{2N_e}}\left(2y_i^{n+1}-\mathrm{\Delta }t{v}_i^{n+1}\right){\displaystyle {\int }_0^1}{\mathrm{\Phi }}_i(x){\mathrm{\Phi }}_j(x)𝑑x={\displaystyle \sum _{i=1}^{2N_e}}\left(2y_i^n+\mathrm{\Delta }t{v}_i^n\right){\displaystyle {\int }_0^1}{\mathrm{\Phi }}_i(x){\mathrm{\Phi }}_j(x)𝑑x\hfill & \\ {\displaystyle \sum _{i=1}^{2N_e}}(2v_i^{n+1}{\displaystyle {\int }_0^1}{\mathrm{\Phi }}_i(x){\mathrm{\Phi }}_j(x)dx+\mathrm{\Delta }t{y}_i^{n+1}{\displaystyle {\int }_0^1}\frac{{\partial }^2{\varphi }_i(x)}{\partial x^2}\frac{{\partial }^2{\varphi }_j(x)}{\partial x^2}dx\hfill & \\ +\mathrm{\Delta }t\alpha v_i^{n+1}{\mathrm{\Phi }}_i(1){\mathrm{\Phi }}_j(1)+\mathrm{\Delta }t\beta v_i^{n+1}\frac{\partial {\varphi }_i(1)}{\partial x}\frac{\partial {\varphi }_j(1)}{\partial x})=\hfill & \\ ={\displaystyle \sum _{i=1}^{2N_e}}(2v_i^n{\displaystyle {\int }_0^1}{\mathrm{\Phi }}_i(x){\mathrm{\Phi }}_j(x)dx-\mathrm{\Delta }t{y}_i^n{\displaystyle {\int }_0^1}\frac{{\partial }^2{\varphi }_i(x)}{\partial x^2}\frac{{\partial }^2{\varphi }_j(x)}{\partial x^2}dx\hfill & \\ -\mathrm{\Delta }t\alpha v_i^n{\mathrm{\Phi }}_i(1){\mathrm{\Phi }}_j(1)-\mathrm{\Delta }t\beta v_i^n\frac{\partial {\varphi }_i(1)}{\partial x}\frac{\partial {\varphi }_j(1)}{\partial x}).\hfill & \end{array}$

2.2. Stability of the numerical scheme

For all, fixed point $x_j$ with $1\le j\le 2N_e$, let us set

(15) $a={\displaystyle {\int }_0^1}{\left({\mathrm{\Phi }}_j(x)\right)}^2𝑑x,b={\displaystyle {\int }_0^1}{\left(\frac{{\partial }^2{\varphi }_j(x)}{\partial x^2}\right)}^2𝑑x,c={\left({\mathrm{\Phi }}_j(1)\right)}^2,d={\left(\frac{\partial {\varphi }_j(1)}{\partial x}\right)}^2.$

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

(16)

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

(17)

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

(18) $\left[\begin{array}{c}{Y}_j^n\\ \\ V_j^n\end{array}\right]=\left[\begin{array}{c}{e}^{ik{x}_j}{\widehat{Y}}_j^n\\ \\ e^{ik{x}_j}{\widehat{V}}_j^n\end{array}\right].$

The system (17) becomes

(19)

The eigenvalues associated with the amplification matrix (19) are ${\lambda }_1=1$ and ${\lambda }_2=\frac{-\mathrm{\Delta }tb+\frac{-4a}{\mathrm{\Delta }t}+2(\alpha c+\beta d)}{\mathrm{\Delta }tb+\frac{4a}{\mathrm{\Delta }t}+2(\alpha c+\beta d)}$. It is clear that these eigenvalues satisfy the Von-Neumann conditions. Indeed, we have $\mid {\lambda }_i\mid \le 1,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,\mathrm{\dots },2N_e$, consider the matrices

${𝕄}_{i,j}$ $={\displaystyle {\int }_0^1}{\mathrm{\Phi }}_i{\mathrm{\Phi }}_j𝑑x,$

${𝕊}_{i,j}$ $=\alpha {\mathrm{\Phi }}_i(1){\mathrm{\Phi }}_j(1)+\beta {({\varphi }_i)}_x(1){({\mathrm{\Phi }}_j)}_x(1),$

${𝕂}_{i,j}$ $={\displaystyle {\int }_0^1}{({\mathrm{\Phi }}_i)}_{xx}{({\mathrm{\Phi }}_j)}_{xx}𝑑x.$

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

(20) $P{W}^{n+1}$ $=Q{W}^n,$

$W^0$ $=W_0$

with

$I_{2N_e}$ denotes the unit matrix of order $2N_e.$
For numerical simulations, we consider as initial conditions:

$y_0(x)$ $=-3x^2+x^3,$

$v_0$ $=0.$

with

$T$ $=10,$

$N_e$ $=70,$

$\mathrm{\Delta }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 $\alpha$ and $\beta$, in order to highlight their influences on the stabilization of the beam. The numerical simulations are carried out with MATLAB.

case 1: We set $\alpha =0.3$ and we vary $\beta .$ case 2: We set $\beta =0.3$ and we vary $\alpha .$

Figure 1. Position for different values of $\beta .$

Figure 2. Rotation angle for different values of $\beta .$

Figure 3. Position for different values of $\alpha .$

Figure 4. Rotation angle for different values of $\alpha .$

3.2. Energy stabilization

Multiply the equation (1) by $\frac{\partial y(x,t)}{\partial t}\in H_l^2(0,1).$ After two integrations by parts, we obtain

${\displaystyle {\int }_0^1}\frac{{\partial }^{2}y(x,t)}{\partial t^2}\frac{\partial y(x,t)}{\partial t}𝑑x+{\displaystyle {\int }_0^1}\frac{{\partial }^{2}y(x,t)}{\partial x^2}\frac{{\partial }^2}{\partial x^2}\left(\frac{\partial y(x,t)}{\partial t}\right)𝑑x+\alpha {\left(\frac{\partial y(1,t)}{\partial t}\right)}^2+\beta {\left(\frac{{\partial }^{2}y(1,t)}{\partial x\partial t}\right)}^2=0.$

So,

(21) $\frac{dE(t)}{dt}+\alpha {\left(\frac{\partial y(1,t)}{\partial t}\right)}^2+\beta {\left(\frac{{\partial }^{2}y(1,t)}{\partial x\partial t}\right)}^2=0.$

Or,

(22) $E(t)=\frac{1}{2}{\displaystyle {\int }_0^1}{\left(\frac{\partial y(x,t)}{\partial t}\right)}^2𝑑x+\frac{1}{2}{\displaystyle {\int }_0^1}{\left(\frac{{\partial }^{2}y(x,t)}{\partial x^2}\right)}^2𝑑x.$

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

(23) $\frac{dE(t)}{dt}\le 0.$

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

Considering the Crank-Nicholson scheme, we have for $0\le n\le N_t,$ the following discrete energy:

(24) $E^n=\frac{1}{2}{\displaystyle {\int }_0^1}{\left(v^n\right)}^2𝑑x+\frac{1}{2}{\displaystyle {\int }_0^1}{\left(\frac{{\partial }^2{y}^n}{\partial x^2}\right)}^2𝑑x.$

We then obtain the following result:

Proposition 1.

For any integer n, there is a positive real ${\eta }_n$ such that

(25) $\frac{E^{n+1}-E^n}{\mathrm{\Delta }t}+{\eta }_n=0.$

Proof.

According to (24), we have

(26) $E^{n+1}-E^n$ $=\frac{1}{2}{\displaystyle {\int }_0^1}\left({\left(v^{n+1}\right)}^2-{\left(v^n\right)}^2\right)𝑑x+\frac{1}{2}{\displaystyle {\int }_0^1}\left({\left(\frac{{\partial }^2{y}^{n+1}}{\partial x^2}\right)}^2-{\left(\frac{{\partial }^2{y}^n}{\partial x^2}\right)}^2\right)𝑑x.$

For $x\in [0;1]$, we replace in the first equation of the variational formulation (9), $\mathrm{\Psi }(x)$ by $\left(v^{n+1}-v^n\right)$. We obtain the following equation

(27) ${\displaystyle {\int }_0^1}\frac{y^{n+1}-y^n}{\mathrm{\Delta }t}\left(v^{n+1}-v^n\right)𝑑x={\displaystyle {\int }_0^1}\frac{{\left(v^{n+1}\right)}^2-{\left(v^n\right)}^2}{2}𝑑x.$

Furthermore, considering $\varphi =y^{n+1}$ as a test function in the second equation of the problem (9), we obtain the following equation:

(28) ${\displaystyle {\int }_0^1}$ $\frac{v^{n+1}-v^n}{\mathrm{\Delta }t}{y}^{n+1}dx+\frac{1}{2}{\displaystyle {\int }_0^1}\left(\left(\frac{{\partial }^2{y}^{n+1}}{\partial x^2}+\frac{{\partial }^2{y}^n}{\partial x^2}\right)\frac{{\partial }^2{y}^{n+1}}{\partial x^2}\right)𝑑x$

$+\frac{\alpha }{2}\left(v^{n+1}(1)+v^n(1)\right)y^{n+1}(1)+\frac{\beta }{2}\left(\frac{\partial v^{n+1}(1)}{\partial x}+\frac{\partial v^n(1)}{\partial x}\right)\frac{\partial y^{n+1}(1)}{\partial x}=0.$

Analogously, if we choose as test function $\varphi =y^n$ in the second equation of the system (9), we obtain:

(29) ${\displaystyle {\int }_0^1}$ $\frac{v^{n+1}-v^n}{\mathrm{\Delta }t}{y}^{n}dx+\frac{1}{2}{\displaystyle {\int }_0^1}\left(\left(\frac{{\partial }^2{y}^{n+1}}{\partial x^2}+\frac{{\partial }^2{y}^n}{\partial x^2}\right)\frac{{\partial }^2{y}^n}{\partial x^2}\right)𝑑x$

$+\frac{\alpha }{2}\left(v^{n+1}(1)+v^n(1)\right)y^n(1)+\frac{\beta }{2}\left(\frac{\partial v^{n+1}(1)}{\partial x}+\frac{\partial v^n(1)}{\partial x}\right)\frac{\partial y^n(1)}{\partial x}=0.$

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

(30) $y^{n+1}-y^n=\frac{\mathrm{\Delta }t}{2}\left(v^{n+1}+v^n\right)$

We obtain

(31) $\frac{E^{n+1}-E^n}{\mathrm{\Delta }t}+{\eta }_n=0,\text{for all}n\in \mathbb{N}.$

With

(32) ${\eta }_n=\alpha {\left(\frac{v^{n+1}(1)+v^n(1)}{2}\right)}^2+\beta {\left(\frac{\partial }{\partial x}\left(\frac{v^{n+1}(1)+v^n(1)}{2}\right)\right)}^2.$

∎

Therefore we have the following result:

Corollary 2.

For all $n\ge 1,$ we have

$$E^{n+1}\le E^n\le E^0$$

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) $E^{n+1}=E^n-\frac{\mathrm{\Delta }t}{4}\left[\alpha {\left(v^{n+1}(1)+v^n(1)\right)}^2+\beta {\left(\frac{\partial }{\partial x}\left(v^{n+1}(1)+v^n(1)\right)\right)}^2\right],$

$0\le n\le N_t$, with

(34) $E^0=\frac{1}{2}{\displaystyle {\int }_0^1}{\left(v^0\right)}^2𝑑x+\frac{1}{2}{\displaystyle {\int }_0^1}{\left(\frac{{\partial }^2{y}^0}{\partial x^2}\right)}^2𝑑x.$

Now let us represent the energy of the system (1)–(4) for different values of $\alpha$ and $\beta .$

case 1: We set $\alpha =0.08$ and we vary $\beta .$

Figure 5. Energy dissipation for different values of $\beta$.

case 2: We set $\beta =0.08$ and we vary $\alpha .$

Figure 6. Energy dissipation for different values of $\alpha$.

3.3. Discussion of results

For a fixed value of $\alpha$ and different values of $\beta$, 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 $\beta$ parameter. However, for a fixed value of $\beta$ and for different values of $\alpha$, 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 $\beta$, the energy decreases and cancels out faster than when the value of parameter $\alpha$ 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) ${\lambda }_n=\left(-\left(\frac{1}{\beta }+\alpha \right)+𝒪\left(\frac{1}{n^2}\right)\right)+i\left({\pi }^2{\left(n-\frac{1}{4}\right)}^2\right)-\frac{2\alpha \beta -1-{(\alpha \beta )}^2}{2\pi {\beta }^{2}n}+𝒪\left(\frac{1}{n^2}\right),n\in {\mathbb{N}}^{\ast }$

and

$$\underset{n\to +\mathrm{\infty }}{lim}\Re 𝔢({\lambda }_n)=-\left(\frac{1}{\beta }+\alpha \right).$$

Notice that $-\left(\frac{1}{\beta }+\alpha \right)$ 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 $\varphi \in H_l^2(0,1)$, we obtain after two integrations by parts the following weak formulation:

(36) ${\displaystyle {\int }_0^1}\frac{{\partial }^{2}y(x,t)}{\partial t^2}\varphi (x)𝑑x+{\displaystyle {\int }_0^1}\frac{{\partial }^{2}y(x,t)}{\partial x^2}\frac{{\partial }^2\varphi (x)}{\partial x^2}+\alpha \frac{\partial y(1,t)}{\partial t}\varphi (1)+\beta \frac{{\partial }^{2}y(1,t)}{\partial x\partial t}\frac{\partial \varphi (1)}{\partial x}=0.$

Also, by separation of variables, the approximate solution in the base ${\left({\mathrm{\Phi }}_i\right)}_{1\le i\le 2N_e},$ can be written as follows:

(37) $y(x,t)={\displaystyle \sum _{i=1}^{2N_e}}{y}_i(t){\mathrm{\Phi }}_i(x).$

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

$$Y={\left[y_1{y}_2\mathrm{\dots }{y}_{2N_e}\right]}^{\top }.$$

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

(38) $$𝕄Y_{tt}+𝕊Y_t+𝕂Y=0.$$

Let $V=\frac{dY}{dt}$ and $W={\left[YV\right]}^{\top }$. Then the equation (38) becomes $\frac{dW(t)}{dt}=LW(t)$ with

the spectrum matrix.

Now, by implementing the spectrum matrix $L$, we represent the spectrum of the system (1)–(4) for different values of $\alpha$ and $\beta .$ 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 $\alpha =0$ and we vary $\beta .$

Figure 7. Effect of the $\beta$ parameter on the spectrum.

Figure 8. Effect of the parameter $\beta$ on the asymptote of the spectrum.

case 2: We set $\alpha =0.08$ and we vary $\beta .$

Figure 9. Effect of the $\beta$ parameter on the spectrum.

Figure 10. Effect of the parameter $\beta$ on the asymptote of the spectrum.

case 3: We set $\beta =0.08$ and we vary $\alpha .$

Figure 11. Effect of the $\alpha$ parameter on the spectrum.

Figure 12. Effect of the parameter $\alpha$ 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 $\alpha$, when we increase the $\beta$ 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 $\beta$ significantly improves the energy decay rate and stability of the system.

However, Figs. 12 and 12 show that for a fixed value of $\beta$, increasing the values of $\alpha$ 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 $\alpha$ 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 $\beta$ has a greater influence on the optimal rate of energy decay and the stability of the system

4.2. Sensitivity of Parameters $\alpha$ and $\beta$ on the Stability

In this part, we objectively study the influence of the parameters $\alpha$ and $\beta$ on the stabilization of the beam. So, for different values of $\alpha$ and $\beta$, we note in the table below, the damping time $\delta \tau$ of the beam vibrations (i.e. time required to return to steady state $y\equiv 0$).

The first-order Sobol sensitivity [6] is written

(39) $$S_i=\frac{𝕍ar\left[𝔼(\delta \tau |Y_i)\right]}{𝕍ar(\delta \tau )},$$

where $𝕍ar(\delta \tau )$ and $𝔼(\delta \tau |Y_i)$) represent respectively the variance of $\delta \tau$ and the conditional expectation of $\delta \tau$ obtained by setting $Y_i=\alpha$ and $Y_i=\beta$ 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 $\alpha$ and $\beta$.

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

$S_{\alpha }$ $\simeq 0.36$

$S_{\beta }$ $\simeq 0.64.$

According to the previous result, we have . We can deduce that the control parameter $\beta$, has more impact on the stabilization of the beam compared to the parameter $\alpha .$

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 $\beta$ 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 $\alpha$. This result was then confirmed by a statistical study of the sensitivity of the control parameters.