Dissipative Numerical Method for a Flexible Euler-Bernoulli Beam with a Force Control in Rotation and Velocity Rotation

In this paper, we study a flexible Euler-Bernoulli beam clamped at one end and subjected to a force control in rotation and velocity rotation. We develop a finite element method, stable and convergent which preserves the property of time decay of energy in the continuous case. We prove firstly the existence and uniqueness of the weak solution. Then, we discretize the system in two steps: in the first step, a semi-discrete scheme is obtained for discretization in space and, in the second step, a fully-discrete scheme is obtained for discretization in time by the Crank-Nicolson scheme. At each step of the discretization, the a-priori error estimates are obtained.


Introduction
In this work, we study a dissipative numerical property by the finite element method for a flexible Euler-Bernoulli beams with a force control in rotation and velocity rotation.The dynamic system that models the mechanical phenomenon that changes over time is described as follows: w (0, t) = w x (0, t) = 0, t > 0, (2) w(x, t) stands for the transverse displacement of the beam at the position x and time t.The subscripts t and x denote derivatives with respect to time t and position x respectively.Moreover, −w xxxx (x, t) dx is the total lateral force acting on a slice of the beam of length dx, located at position x and time t and w xx (1, t) is the force in rotation acting on the rigid body from the beam at the time t.The nonnegative constants α and β are the feedback gains that can be tuned in practice.
For simplicity sake, the flexural rigidity function, the mass density function of the beam and the length of the beam are assumed to be unity.Moreover, the following notation v = w t (velocity of the beam) will be used in the sequel.
In (Touré, Koua & Diop, 2016), it has been proved that the system (1)-( 4) is well posed in the sense of C 0 -semigroup of contractions.In order to perform the stability analysis of this system, the authors formulate the problem as an evolution problem first.Also, from Shkalikov's method (Shkalikov, 1986), a spectral analysis of the operator and the property of the Riesz basis were studied to derive the exponential stability of the system.
Furthermore, it should be noted that for all t, in our case, the total mechanical energy ε : R + → R + of the system (1)-( 4) is given by which decreases over time.Indeed, the time derivative of the energy functional ε(t) along the classical solutions of (1)-( 4) read as follows: because α ≥ 0. The right hand side of (6) serves as a motivation in the design of the control αw xt (1, t) + βw x (1, t), ensures the energy decay of the system in time.Moreover, according to Theorem 3.5 of (Touré, Koua & Diop, 2016) whose proof is based on an idea of (Guo, 2002), the system (1)-( 4) is exponentially stable for any β > 0 and α ≥ 0. Our main contribution is to develop a convergent numerical method which faithfully reproduces some properties of this problem such as stability and energy decay.
The rest of the paper is organized as follows.In section 2, from the weak formulation, we show the existence, uniqueness and higher regularity of the weak solution.In section 3 and in section 4, we develop by finite element method, a numerical method for the system (1)-( 4) which conserves the dissipativity property.
2. Existence, Uniqueness and Higher Regularity of the Weak Solution

}
where Then, we also introduce the following functional space: and for energy space, the following hilbert space : where the superscript T stands for the transpose.In the space χ, we define the inner-product: < y, y >= 1 2 where y = (w, v) T ∈ χ and y = ( w, v) T ∈ χ.We denote by ∥.∥ χ the associated norm.Next, we define an unbounded linear operator A : D (A) ⊂ χ→χ as follows: where D (A) , the domain of operator A is as follows Now we can write our problem as a first order evolution problem where y (t) = (w (., t) , v (., t)) T , y (0) = (w 0 , v 0 ) T for all t > 0.
Notice that the contractivity of the semigroup also implies that ∥.∥ χ is a good candidate for the Lyapunov functional for (12).Let the functional ℓ : χ → R defined as follows Analogously as in ( 6), for all classical solutions y it follows that: hence time evolution of the Lyapunov functional along the classical solutions is non-increasing.Furthermore, from Theorem 1, the decay of energy along the classical solutions can be extended to mild solutions : Theorem 3 Assume that y(t) is the mild solution of (12) for all y 0 ∈ χ.Then y(t) → 0 in χ when t → ∞.
Remark 1 The right hand side of ( 14) is formed by a control variable.Thus, d dt ∥y∥ 2 χ = 0 does not imply y = 0 through (12).Now write the system of equations ( 1)-(4) in the weak form.

Weak Formulation
Let ϕ ∈ V. Multiplying (1) by ϕ, integrating over [0, 1] and taking into account the given boundary conditions (2)-(4), we have It seeks to define a weak solution of (1)-(4).But first, make an appropriate choice of spaces.We follow an idea used in (Banks & Rosen, 1987).Let the Hilbert space Y = R 2 × L 2 (0, 1) with the following inner product: We also define the following Hilbert space X = R 2 × V = {y = (w(1), w x (1), w); w ∈ V} with the inner product It can easily be shown that X is densely embedded in Y and suppose that the canonical injection of X into Y is continuous.Therefore, taking Y as a pivot space we obtain a Gelfand triple : where X ′ is the dual of X.Consider the following bilinear forms: In the following definition, the bilinear form < ., .> X,X ′ is the duality pairing between X and X ′ , which is a natural extension of the inner product in Y.
Definition 4 Let T > 0 be fixed.We say that w = (w(1), w x (1), w) is a weak solution of problem (1)-( 4) for almost everywhere on t ∈ (0, T ) and for all ϕ ∈ X, with the following initial conditions w(0) = w 0 = (w 0 (1), (w 0 ) x (1), w 0 ) ∈ X (19) Remark 2 The formulation ( 18) is equivalent to equality (15) if w ∈ H 2 (0, T, X).Furthermore, in the expression (19), the first two components of the right hand side are the boundary traces of w 0 ∈ V.But for (v 0 )(1) and (v 0 ) x (1) in ( 20), they are given in addition to the function v 0 and not as its trace.Here, the term u tx (1) also need to be considered.Then, the bilinear form a 2 (., .)with the first order boundary term in t requires a slight generalization of the standard theory (as presented for example in chapter 3 section 8 of (Lions & Magenes, 1968) or again in section 7.2 of (Evans, 1998)).
Recall the following Lemmas provided in Theorems 3.1 pp.23 and 6.2 pp.34 of (Lions & Magenes, 1968) where the definition of the intermediate spaces is given and which will be useful in proving theorem of existence of weak solution.
Notice that, in the following, [X, Y] θ , with 0 ≤ θ ≤ 1 is the intermediate space defined as in chapter 1 pp.11 − 13 of (Lions & Magenes, 1968), for X and Y Hilbert spaces, X ⊂ Y, X dense in Y with continuous injection by means of domains of positive self-adjoint operators.Remark also that for θ Lemma 5 Let X and Y be two Hilbert spaces, such that X is dense and continuously embedded in Y. Assume that w ∈ L 2 (0, T ; X) and v ∈ L 2 (0, T ; Y).Then w ∈ C([0, T ] ; [X, Y] 1 2 ), after, possibly, a modification on a set of measure zero.Since X ⊂ [X, Y] θ ⊂ Y, each space being dense in the following, we have by duality (without any identification between space and its dual) for θ ∈ ]0, 1[ : , each space being dense in the following.We have the following duality theorem: Lemma 6 Let X and Y be two Hilbert spaces, such that X is dense and continuously embedded in Y.For all θ ∈ ]0, 1[, [X, Y] ′ θ = [Y ′ , X ′ ] 1−θ holds.We will also use the following result: Theorem 7 Let V be a subspace of H 2 (0, 1).Then there exists a infinite sequence of functions {ϕ i } ∞ i=1 such that Consider the following boundary value problem : Assume that f ∈ L 2 (0, 1).Let ϕ ∈ V. Multiplying by ϕ, integrating twice by parts and taking into account the given boundary conditions yields: The fact that f ∈ L 2 (0, 1) ensures the continuity of the linear form.Set Then a 3 is symmetric bilinear form, bounded and coercive on V.The Lax-Milgram theorem allows us to conclude the existence and uniqueness of the solution w of (21).Then, there exists a unique weak solution w in V such that w = B −1 ( f ) with B −1 : L 2 (0, 1) → L 2 (0, 1).We remark that B −1 is obviously linear and bounded.Furthermore, there exists a constant and since V is compactly embedded in L 2 (0, 1) then B −1 is compact.It remains to show that it's symmetric.Let f, g ∈ L 2 (0, 1) and denote w = B −1 ( f ) and v = B −1 (g).By straightforward calculation and integration by parts, we finally obtain Analogously, we get Thus, a 3 is symmetric and we have Hence B −1 is symmetric.In view of compactness and symmetry properties of the bounded linear operator B −1 , there exists a countable orthonormal basis {ϕ i } ∞ i=1 of L 2 (0, 1) constitued of eigenvectors B −1 .Furthermore, these eigenvectors are functions of V according to definition of B −1 .Moreover, from the weak formulation, one can see that the basis {ϕ i } ∞ i=1 is an orthogonal basis of V with respect to the inner product a 3 (., .). 2

Existence of the Weak Solution
Theorem 8 There exists a weak solution w of the equivalent weak formulation (18) such that : ). ( 26) Proof.This proof is based on the Faedo-Galerkin's method and is an adaption of the proof of Theorem 8.1 pp.287-290 in (Lions & Magenes, 1968).According to Theorem 7, there exists by extension an infinite sequence of functions that is an orthogonal basis for X and an orthonormal basis for Y. Consider such a sequence.Introduce the following finite dimensional spaces spanned by Step 1 (Construction of approximate solutions): We seek w = w m (t) ∈ V m the approximate solution of the problem.Then w is in the form: where m) are solutions of the formulation (18) on V m .For a fixed m ∈ N, we have: The approximate differentials equations system (28) is completed with the initial conditions: with α im = g im (0) and β im = (g im ) t (0).According to standard existence theory for ordinary differential equations, we are insured of the existence of solution w m ∈ C 2 ([0, T ] ; X) of ( 28)-( 30) for 0 ≤ t ≤ T.
Step 2 (A-priori estimates on approximate solutions): Let E : R × X → R be the analogue of the Lyapunov functional as defined by ( 13): Assuming that there exists a solution w m ∈ C 2 ([0, τ] ; V m ) to ( 28) on some interval [0, τ] and taking ϕ = ( w m ) t in (28), a straightforward calculation yields Dissipation of the functional E corresponds to the decay in ( 14) for the classical solution.This implies uniform boundedness of the solution on [0, τ]: It remains to find set in which ( w m ) tt is bounded.Considering the results ( 34)-( 35), it is shown that for all ϕ ∈ X : where M is a positive constant which does not depend on m.Now, let m ∈ N be fixed.Furthermore, let ϕ ∈ X and 28) and ( 36), we have: Step 3 (Passage to the limit): According to the Eberlein-Šmulian Theorem (see Brezis, 2011, e.g.), we can extract weakly convergent subsequences Furthermore, (40) yields where µ j ∈ L 2 (0, T ; R) and for all m l ≥ m 0 , the formulation (28) becomes: Therefore, passing on to the limit in ( 44) for m = m l , when l → ∞ and using the convergence results ( 39)-( 41), one obtains: Then, one obtains < w tt , φ > X,X ′ +a 1 ( w, φ)+a 2 ( w t , φ) = 0 a.e on [0, T ] for all φ ∈ L 2 (0, T ; X).Since the functions φ of the form (43) are dense in L 2 (0, T ; X), w is therefore the solution of the weak formulation.For the additional regularity, from the construction of the weak solution and due to (34)-( 35), w satisfies (24).Furthermore w satisfies (25) using the Lemma 5, after, possibly a modification on a set of measure zero, and the regularity (26) follows from Lemma 5 and Lemma 6. 2

Uniqueness of the Weak Solution
Theorem 9 The solution w of the weak formulation (18) with the initial conditions (19)-( 20) is unique.
Proof.This proof of uniqueness is an adaption of the proof of Theorem 8.1 pp.290-291 in (Lions & Magenes, 1968).

Higher Regularity Results
Before stating the theorem of the stronger continuity of the weak solution, recall the lemma 8.1 of chapter 3 pp.297 of (Lions & Magenes, 1968) which will be used in the proof of this theorem.
Lemma 10 Let X and Y two Banach spaces, X ⊂ Y with continuous injection, the space X being reflexive.We set: which denotes the space of weakly continuous functions with values in Y. Thus we get Proof.For the proof, the reader is referred to chapter 3 pp.297−298 in (Lions & Magenes, 1968). 2 Theorem 11 The weak solution w of (18)-( 20) satisfies w ∈ C([0, T ] ; X), (52) after possibly a modification on a set of measure zero.
Proof.This proof is an adaption of standard strategies in section 8.4 of (Lions & Magenes, 1968) pp.297-301 and in section 2.4 of (Temam, 1988).Using Lemma 10, it results from ( 26)-( 27) that w ∈ C w ([0, T ] ; X).Furthermore, ( 24) and ( 26) imply w t ∈ C w ([0, T ] ; Y).Now, we use a common technique in functional analysis, specifically in distribution theory, to move from a problem of generalized functions to a restriction of regular functions easier to handle.Let a scalar cutoff function ξ ∈ C ∞ (R) be fixed such that ξ(x) = 1 if x ∈ J ⊂⊂ [0, T ] and ξ(x) = 0 else.The function ξ w is then compactly supported.Let η ε be a standard mollifier in time.For example, the function η ε may be given by η ε (t ) where In addition, w ε converges to w in X and ( w ε ) t converges to w t a.e in H for all element on J. Hence, E(t, w ε ) converges to E(t, u) a.e on J. Since w ε is smooth, a straightforward calculation on J gives: Passing to the limit, one obtains when ε → 0 : in the sense of distributions on J. Since J was arbitrary, (54) holds on all compact subintervals of [0, T ].Then, let t ∈ [0, ∞[ be fixed, and let lim n−→∞ t n = t.Let the sequence ν n be defined by Thus, we have: Since w, w t are weakly continuous and E is continuous in t, we have, passing to the limit in (56) : Thus we get w ∈ C([0, T ] ; X) and w t ∈ C([0, T ] ; Y). 2 In the next sections, the goal is to develop a numerical method for (1)-( 4) in such a way that the decay of the Lyapunov function is preserved.The first step of this method is the discretization of the system in space to obtain the semi-discrete scheme, and then in time, in order to get the fully-discrete scheme.

Piecewise Cubic Hermite Polynomials
Our numerical work consist to construct an appropriate piecewise space of C 1 -functions on In particular, if the subdivision is uniform, let us denote the step length by h = 1/P and Let us find cubic polynomial functions N i j , i = 0, 1, ..., P − 1 and j = 1, 2, 3, 4 that satisfy the following conditions : Let us use the following affine transformation: allowing us to manipulate all operations on [0, 1].The intermediate variables θ i are called local coordinates.Then under these coordinates, the unknown functions N j (θ) must satisfy the following boundary conditions: Let us find explicit expressions of the functions N j (θ).For all θ ∈ [0, 1] , for all j = 1, 2, 3, 4, the functions N j (θ) are polynomials of degree 3, therefore of the form N j (θ) = aθ 3 + bθ 2 + cθ + d.We remark that 1 is a double root of N 1 .Then N 1 (θ) = (θ − 1) 2 (aθ + b).Using the two remaining conditions, we find a = 2 and b = 1.Hence we get Analogously, we find

Figure 2. Hermitian Polynomial functions
By extension by 0 on Z − Z i , for all i = 0, ..., P − 1, we define the Hermitian functions on Z : By concatenating the Hermitian functions obtained above, we can obtain the basic functions on the different supports as follows: with the support Z 0 , ϕ 0 x) and with the support Z P , ϕ P 1 (x) = N P−1 3 (x), ϕ P 2 (x) = N P−1 4 (x).The set B = {ϕ k l , k = 0, ..., P, l = 1, 2} forms a basis that generates a subspace of V of dimension 2P + 2 denoted by V h .With the separation of variables, the approximate solution w h ∈ V h which we seek can be written as follows: Thus, since w 0 h = w 0 h = 0, the N-dimensional space is as follows:

Semi-discrete Scheme: Space Discretization
Since rounding errors are cumulative because of the high partial derivative terms (the spatial derivatives being of order 4) in the equation of the beam (1), we propose a finite element scheme semi-discretized in space.In fact, a variational approach allows us to reduce the degree of the derivatives by the integration by parts.Let ϕ j , j = 1, ..., N be fixed basis for V h .The semi-discrete solution ) is defined as the solution of the finite element method: for all j = 1, ..., N and t > 0, which solves the initial conditions Equation ( 62) is a second order ODE-system in time.By separation of variables, its solution can be written in the following form: where W is a vector representation of the function w h defined as follows: Equation ( 62) is equivalent to the following equation: M is the mass matrix and K is the rigidity matrix.The corresponding matrices M, S and K are given by: The matrix K is symmetric, defined and positive because β > 0 and therefore K is invertible.Since the matrix M is also symmetric, defined and positive , this implies the existence and the uniqueness of the solution of problem ( 62)-( 64).Note that M and K are tridiagonal matrices by blocks while S is diagonal.The calculation of elements of S is trivial because all the elements of S are zero except one nonzero element S N,N = α with N = 2P.
Values of elements of matrices M and K

Dissipativity of the Semi-discrete Scheme
In order to show that the scheme given by ( 62)-( 64) is dissipative, first a time dependent energy functional E for a trajectory w ∈ C 2 ([0, ∞[ , V) is defined as analogous of the Lyapunov functional (13).
Theorem 12 Let w h ∈ C([0, ∞[ ; V h ) solution of (62)-( 64).Then we get: Proof.Derive for all t > 0, E(t, w h ).We have: Using ( 62) with the test function ϕ h = (w h ) t , we obtain: Then, we get the result (67). 2 Remark 3 Note that the property of dissipativity theorem of the norm was written independently of the basis ϕ j , j = 1, ..., N and can be applied to any choice of the subspace V h ⊂ V.

A-priori Error Estimates
In this subsection, the a-priori error estimates for the semi-discrete solution approximation (62) are obtained.We will use a common method used in (Choo, Chung & Kannan, 2002) to obtain error estimates.Of course, we will adapt the method used in this article to our problem.The projection of weak solution w to V h on H 2 (0, 1) denoted by w is defined as follows: We set G = {w ∈ H 4 (0, 1), w(0) = w x (0) = 0}.Assume for later that: Then, from the lemma 2.1 of (Choo, Chung & Kannan, 2002), we have the following estimations almost every in t: Now, we give an important result for the convergence of the semi-discrete scheme: Theorem 13 Let V h the space of cubic Hermite polynomials.Assume the expressions (70)-( 72).The following error estimate holds for w h ∈ C 2 ([0, T ] , V h ) solving (62) is given by: ) ) .
Furthermore, if w h0 and v h0 are respectively Hermite interpolations of w 0 and of v 0 , then there exists a positive constant C such that: ) . (77) Proof.The error of semi-discrete solution w h is defined as e h = w h − w.We remark that e h is an element of V h .Then, substituting w h = e h + w in (62), we get: for all ϕ ∈ V h .Furthermore, w is the projection of weak solution w on discret espace V h .Therefore using again (62), we have the following equation: for all ϕ ∈ V h .Taking now ϕ = (e h ) t ∈ V h .Thereby, (78) becomes : Thus we get: This implies that: 1 2 Or again : By integrating (81) in the time direction ie on t ∈ [0, T ] , one obtains: By one integration by parts ∫ t 0 (w − w) xx (e h ) txx dτ, (82) becomes finally: Using Cauchy Schwarz's inequality to (83) yields where C 1 , C 2 , C 3 and C 4 are positive constants.

Fully-discrete Scheme
In this section, in order to obtain a fully discretized scheme, we discretize in time of semi-discretized system (62) in such a way that the dissipation of energy is preserved.To achieve this goal, the system ( 62) is written as a system of ordinary differential equations of first order.Then, the Crank-Nicolson scheme obtained is used to demonstrate dissipativity of the numerical scheme.Finally, the a-priori estimates are obtained.

Crank-Nicolson Scheme
Let L be positive integer.Here, the interval [0, T ] is discretized into L equidistant subintervals.Let k = T/L denotes the size of time discretization and t n = nk where n = 0, 1, ..., L represent the nodes of the discretization.In order to rewrite the semi-discretized scheme (62) as a differential equation of the first order, introduce the element T be its vector representation in the basis of Hermitian cubic polynomial space.Note that the solution w h of semi-discretized scheme (62) becomes for the full discrete scheme a vector of the form y h = [w h v h ] T .Furthermore, similar to (13), the natural norm of y h = y h (t) is defined as follows: Let now y n = [w n v n ] be the approximate solution of y h in time t = t n .Let again W n = W(x, t n ) and V n = V(x, t n ) the vector representations in basis respectively of w n and v n .For the time discretization of ( 62), the Crank-Nicolson scheme is used.Then we get: for all ϕ h ∈ V h .Furthermore, the vector equation (66) becomes: which is equivalent to: (84) and ( 87) give us the following system of equation: PY n+1 = QY n where P and Q are block matrices defined as follows:

Dissipativity of Numerical Scheme
Now, we show that the fully discrete scheme ( 84) and ( 86) dissipates the norm (energy) in time.
Theorem 14 For all n = 0, 1, ..., L, L and k the positive integers, we get Proof.We have Multiplying (84) by v n+1 − v n and integrating over [0, 1], we obtain: Taking ϕ h = w n+1 ∈ V h in (85) yields: 1 2 x (1).( 90) Taking now ϕ h = w n ∈ V h in (85).Hence we get: Thus, using ( 89)-( 91) and again (84), we obtain the result (88). 2 Remark 5 The norm dissipates in time : ∥y n+1 ∥ 2 ≤ ∥y n ∥ 2 .This decay of the norm when k −→ 0 corresponds to the decay of the norm (14) in the continuous case and with the norm (67) in the semi-discrete case.However, if the beam is not controlled ie when α = β = 0 then ∥y n+1 ∥ = ∥y n ∥, therefore the norm ∥y n ∥ is constant for all n = 0, 1, ..., L where L is a positive integer.Also, notice that the Crank Nicolson scheme (84), the expression(86) and the norm dissipation property from Theorem 14 of the norm were written independently of the basis { ϕ j } .Therefore, this property of dissipativity can be applied to any choice of the subspace V h ⊂ V.

A-priori Error Estimates
Assume that w ∈ H 4 (0, T ; V).Let w ∈ V h be defined as the projection of the weak solution w on V h such that a 1 (w Moreover, let w e := w − w denote the error of the projection.Assume also that: Then due to ( Strang & Fix, 1973), we have the following estimations: T denotes the weak solution of (18) at time t = t n .This approximation is defined by y n = [w n v n ] T , the n-th iteration of the fully-discrete scheme of Crank Nicolson.Thus, the approximation error is defined by Therefore, the second order error estimate both in space and time of the fully discrete scheme is obtained in the following Theorem.
We have the following estimate: In order to obtain the estimate of Q n 2 , we need to rewrite the second term of Q n 2 (e n 2 ).Integrating twice by parts over [0, 1] and assuming that e n 2 (0) = e n 2x (0) = 0, we have: ) .