# A Numerical Method of the Euler-Bernoulli Beam with Optimal Local Kelvin-Voigt Damping

Authors: Xin Yu; Zhigang Ren; Qian Zhang; Chao Xu
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

This paper deals with the numerical approximation problem of the optimal control problem governed by the Euler-Bernoulli beam equation with local Kelvin-Voigt damping, which is a nonlinear coefficient control problem with control constraints. The goal of this problem is to design a control input numerically, which is the damping and distributes locally on a subinterval of the region occupied by the beam, such that the total energy of the beam and the control on a given time period is minimal. We firstly use the finite element method (FEM) to obtain a finite-dimensional model based on the original PDE system. Then, using the control parameterization method, we approximate the finite-dimensional problem by a standard optimal parameter selection problem, which is a suboptimal problem and can be solved numerically by nonlinear mathematical programming algorithm. At last, some simulation studies will be presented by the proposed numerical approximation method in this paper, where the damping controls act on different locations of the Euler-Bernoulli beam.


---

## Body

## 1. Introduction

LetT,L>0 be two positive constants. We denote by Q the product set (0,L)×(0,T). Consider a nonhomogeneous clamped elastic beam of length L, where one segment of the beam is made of a viscoelastic material with Kelvin-Voigt constitutive relation. By the Kirchhoff hypothesis, neglecting the rotatory inertia, the transversal vibration of the beam can be described by the following equation and boundary-initial conditions:
(1)ρw¨+(qw′′+Dw˙′′)′′=f,(x,t)∈Q,w(0,t)=w(L,t)=w′(0,t)=w′(L,t)=0,w(x,0)=w0(x),w˙(x,0)=w1(x),
where w represents transversal displacement of the beam, w0,w1∈L2(0,L) are given initial data, and the notations w˙ and w′ denote the derivatives with respect to the temporal variable and the spatial variable of w, respectively. Here f∈L2(Q) is an external applied distributed force, ρ is the linear mass density of the beam material, q is the flexural rigidity, and D is the Kelvin-Voigt damping coefficient. In this paper, we assume ρ,q∈L∞(0,L), D∈L∞(Q), and ρ,q≥c>0 for x∈(0,L) for some constant c. we say that the Kelvin-Voigt damping is globally distributed if the damping coefficient D≥c>0 on [0,L]; we say it is locally distributed if D≥c>0 only on some subinterval of [0,L] and D=0 elsewhere.Smart materials such as shape memory alloys and piezoceramics [1–3] have been applied for active vibration control of elastic structures. Accordingly, one can introduce control terms to the elastic systems such as the damping coefficients and Young’s moduli as well. Magnetorheological (MR) dampers [4–7] are one of the most promising new actuation mechanisms that use MR fluids to provide variable damping actuation for active control of structures. Because of their mechanical simplicity, high dynamic range, low power requirements, large force capacity, and robustness, these devices have been shown to mesh well with application demands and constraints to offer an attractive control method to structural vibration.In this paper, we will study the Euler-Bernoulli beam equation with optimal local Kelvin-Voigt damping. To be more specific, letω∈[0,L] be a subinterval and let χω be the characteristic function of ω, that is,
(2)χω(x)={1,x∈ω,0,otherwise;
We define
(3)U={u(·)∈L∞(0,T)∣R1≤u(t)≤R2fora.e.t∈[0,T]}
and assume the damping coefficient D has the following form:
(4)D(x,t)=χω(x)u(t),(x,t)∈Q,
where R1,R2∈R are two fixed constants and u∈U is a control function.LetV=H02(0,L) with the norm
(5)∥w∥V=(∫0L‍q(x)|w′′(x)|2dx)1/2,∀w∈V,
and H=Lρ2(0,L) with the norm
(6)∥v∥H=(∫0L‍ρ(x)|v(x)|2dx)1/2,∀v∈H.
Define H=V×H with the norm
(7)∥(wv)∥H=(∥w∥V2+∥v∥H2)1/2,∀(wv)∈H.
Then, H is Hilbert space and the energy of the beam at time t is
(8)∥(w(t)w˙(t))∥H2=∫0L‍q(x)|w′′(x,t)|2dx+∫0L‍ρ(x)|w˙(x,t)|2dx,
where w is the solution of (1). The optimal control problem that we will study is formulated as follows:
(9)(OCP)min⁡u∈U{J(u)=12∫0T‍∥(w(t)w˙(t))∥H2dt+12∫0T‍|u(t)|2dt∥(w(t)w˙(t))∥H2}
subject to the controlled equation (1), where w is the solution of (1). Throughout the paper, we will omit the notations t or x in the functions of t or x in the case that there is no risk to make any confusion.Due to the importance from both perspectives of mathematics and applied science, the control problem of various beam equations has been considered by many researchers [8–10]. The study of the Euler-Bernoulli beam is one of the most active research topics in control theory. In [11], the authors consider the vibration of the Euler-Bernoulli beam with Kelvin-Voigt damping distributed locally on any subinterval of the region occupied by the beam. By making use of the frequency domain method and the multiplier technique, they prove that the semigroup associated with the equation for the transversal motion of the beam is exponentially stable. In [12], the author studies the basis property and the stability of a distributed system described by a nonuniform Euler-Bernoulli beam equation under linear boundary feedback control. The Riesz basis property is presented and the exponential stability is concluded. In [13], stabilization of Euler-Bernoulli beam by means of the pointwise feedback force is considered. Both uniform and nonuniform energy decay may occur, which depend on the boundary conditions. There are some other related papers about the studies of Euler-Bernoulli beam equations [14–17].In this paper, we will study the numerical approximation of the optimal control problem(OCP), which is a nonlinear bilinear control problem with control constraint. Bilinear control problems are already studied by many researchers [18–20]. In our paper, we want to design a damping control numerically, which acts on local interval of the beam, such that the total energy of the beam and the control on a given time period is minimal. It appears that little work has been done on numerical methods for this problem. By the standard finite element method (FEM), problem (OCP) was firstly approximated by an optimal control problem (OCPh) governed by a system of ordinary differential equations. Then, using the control parameterization method [21], we will approximate the finite-dimensional problem by another suboptimal problem (OCPdh), which is a standard optimal parameter selection problem and can be solved numerically by nonlinear mathematical programming algorithm. At last, some simulation studies will be presented by the numerical method proposed in this paper.


## 2. The Semidiscrete Approximation by FEM

In this section, we will approximate the original optimal control problem(OCP) with FEM method. Noting that (1) involves the spatial derivative of four orders, the conforming FEM space should belong to H2(0,1). Consider the interval domain [0,L]. The triangulation Th of [0,L] divides [0,L] into a finite number of subintervals Ij=[xj-1,xj], j=1,2,…,N+1, using the grid points:
(10)0=x0<x1<⋯<xN<xN+1=L,
where we will call Ij=[xj-1,xj] the jth element and hj=xj-xj-1 the size of this element. The discretization parameter h is the maximum size of all hj, j=1,2,…,N+1. Associated with every triangulation Th, we define a finite-dimensional space as follows:
(11)Vh={(vh)′vh∈C1[0,L]∣vh∈P3(Ij)forj=1,2,…,N+1,vh(0)=vh(L)=(vh)′(0)=(vh)′(L)=0},
where P3(Ij) is the space of all polynomials of degree less than or equal to 3 over the subinterval Ij. Obviously, we have Vh∈H02(0,L). Thus, we can write
(12)Vh=span⁡{ϕi,∣i=1,2,…,2N},
where, for i=1,2,…,N,
(13)ϕi(xk)=δik,ϕi′(xk)=0,∀k=0,1,…,N+1,
and, for i=N+1,N+2,…,2N,
(14)ϕi(xk)=0,ϕi′(xk)=δi-N,k,∀k=0,1,…,N+1.Define a bilinear form(·,·)ρ over Lρ2(0,L)×Lρ2(0,L) by setting
(15)(f,g)ρ=∫0L‍ρ(x)f(x)g(x)dx,∀f,g∈Lρ2(0,L).
Define another two bilinear forms aq(·,·) and aω(·,·) over H02(Ω)×H02(Ω) by setting
(16)aq(f,g)=∫0L‍q(x)f′′(x)g′′(x)dx,∀f,g∈H02(Ω),(17)aω(f,g)=∫ω‍f′′(x)g′′(x)dx,∀f,g∈H02(Ω),
respectively. Obviously, the two bilinear forms (·,·)ρ and aq(·,·) are the inner products of Lρ2(0,L) and H02(Ω), respectively. Then, the finite element approximation of (1) consists in finding wh(t)=wh(·,t), which belongs to Vh for t∈[0,T], and satisfies
(18)(w¨h(t),vh)ρ+aq(wh(t),vh)+u(t)aω(w˙h,vh)=(f(t),vh)∀vh∈Vh,0<t≤T,wh(0)=w0h,w˙h(0)=w1h,
where (·,·) denotes the standard inner product of L2(0,L) and w0h,w1h∈Vh are the proper approximations of w0,w1 on Vh. In the following, we write
(19)wh(x,t)=∑j=12N‍Xj(t)ϕj(x),(20)wh(x,0)=w0h(x)=∑j=12N‍X0jϕj(x),w˙h(x,0)=w1h(x)=∑j=12N‍Y0jϕj(x).
Substituting (19) into (18) and taking vh=ϕi yield that
(21)∑j=12N‍(ϕi,ϕj)ρX¨j(t)+∑j=12N‍aq(ϕi,ϕj)Xj(t)+u(t)∑j=12N‍aω(ϕi,ϕj)X˙j(t)=(f(t),ϕi).
Moreover, in this paper, we take w0h,w1h as the Lρ2-projection approximations of w0,w1 on Vh; that is,
(22)(w0h,vh)ρ=(w0,vh)ρ,(w1h,vh)ρ=(w1,vh)ρ,∀vh∈Vh.
Then, substituting (20) into (22) and taking vh=ϕi yield that
(23)∑j=12N‍(ϕi,ϕj)ρX0j=(w0,ϕi)ρ,∑j=12N‍(ϕi,ϕj)ρY0j=(w1,ϕi)ρ.
Define
(24)X(t)=[Xi(t)]2N×1,M=[(ϕi,ϕj)ρ]2N×2N,K=[aq(ϕi,ϕj)]2N×2N,R=[aω(ϕi,ϕj)]2N×2N,X0=[X0i]2N×1,Y0=[Y0i]2N×1,F(t)=[(f(t),ϕi)]2N×1,η=[(w0,ϕi)ρ]2N×1,ξ=[(w1,ϕi)ρ]2N×1.
Thus, by (21) and (23), we can obtain the following system of controlled ordinary differential equations:
(25)MX¨(t)=-KX(t)-u(t)RX˙(t)+F(t),X(0)=X0,X˙(0)=Y0,
where X0=M-1η, Y0=M-1ξ. Let
(26)X˙=Y,Z=(XY).
We define
(27)G(t,Z(t),u(t))=(0I-M-1K-u(t)M-1R)Z(t)+(0M-1F(t)).
Then (25) can be rewritten as
(28)Z˙(t)=G(t,Z(t),u(t)),Z(0)=Z0,
where Z0=(X0Y0).By (19), a direct computation yields
(29)∥wh(t)∥V2=aq(wh(t),wh(t))=∑i=12N‍∑j=12N‍aq(ϕi,ϕj)Xi(t)Xj(t)=(X,KX)R2N,
where (·,·)R2N denotes the inner product of R2N and V=H02(0,L). Similarly, we have
(30)∥w˙h(t)∥H2=(w˙h(t),w˙h(t))ρ=∑i=12N‍∑j=12N‍(ϕi,ϕj)ρX˙i(t)X˙j(t)=(Y,MY)R2N,
where H=Lρ2(0,L). Define
(31)Π=(K00M).
Then, it follows from (29), (30), and (31) that
(32)J(u)≈Jh(u)=12∫0T‍∥(wh(t)w˙h(t))∥H2dt+12∫0T‍|u(t)|2dt=12∫0T‍(Z(t),ΠZ(t))R4Ndt+12∫0T‍|u(t)|2dt.
Thus, by combining (28) and (32), the semidiscrete approximation of problem (OCP) is formulated as follows:
(33)(OCPh)min⁡u∈UJh(u)
subject to (28).


## 3. Piecewise-Constant Control Approximation

In general, problem(OCPh) cannot be solved analytically. Using the control parameterization method, which has been successfully applied to provide numerical solutions for a wide variety of practical optimal control problems [21–24], we will approximate problem (OCPh) by a standard optimal parameter selection problem. This method involves approximating the control function by a piecewise-constant function with possible discontinuities at a set of preassigned switching points, which produces an approximation problem such that the solution of this approximation is a suboptimal solution to problem (OCPh).Letτk, k=0,1,…,d, be prefixed time knot points satisfying
(34)0=τ0<τ1<τ2<⋯<τd=T.
With piecewise-constant basis functions, the control input u for the problem (OCPh) is approximated over the kth control subinterval [τk-1,τk) as follows:
(35)u(t)≈ud(t)=μk,t∈[τk-1,τk),k=1,2,…,d,
where μk is the value of the control on the kth subinterval [τk-1,τk). Define
(36)Ud={μ=(μ1,μ2,…,μd)⊤∈Rd∣R1≤μk≤R2,k=1,2,…,d(μ1,μ2,…,μd)⊤}.
Then the approximate piecewise-constant control can be written as follows:
(37)ud(t∣μ)=∑k=1d‍μkχ[τk-1,τk)(t),t∈[0,T],μ∈Ud,
where χ[τk-1,τk) is the characteristic function of the interval [τk-1,τk), k=1,2,…,d. Substituting (37) into the dynamic system (28) yields that
(38)Z˙(t)=G(t,Z(t),ud(t∣μ))=∑k=1d‍G(t,Z(t),μk)χ[τk-1,τk)(t),t∈[0,T],Z(0)=Z0.
Let Zd(·∣μ) denote the solution of system (38) corresponding to μ∈Ud. Thus, from the problem (OCPh), we can obtain another parameter optimization problem, which is stated as follows:
(39)(OCPdh)min⁡u∈UdJdh(μ)
subject to (38), where
(40)Jdh(μ)=12∫0T‍(Zd(t∣μ),ΠZd(t∣μ))R4Ndt+12∫0T‍|ud(t∣μ)|2dt.After the parameterization of control, problem(OCPdh) involves a finite number of decision variables. Thus, it should be much easier to solve than problem (OCPh), which involves determining the value of a function at an infinite number of time points.


## 4. Variational Method for Solving Problem(OCPdh)

Problem(OCPdh) is an optimal parameter selection problem in the canonical form [24], which can be solved as nonlinear optimization problems using the SQP method. Standard SQP algorithm for nonlinear optimization exploits the gradient of the cost functional to generate search directions that lead to profitable areas of the search space [25, 26]. For the approximate problem (OCPdh), the cost functional is implicit function of the decision vector μ. Using the variational method [27, 28], we can compute this gradient and solve problem (OCPdh).For eachm=1,2,…,d, it follows from (38) that
(41)Zd(t∣μ)=Zd(τm-1∣μ)+∫τm-1t‍G(t,Zd(s∣μ),μm)ds,t∈[τm-1,τm).
Then for t∈[τm-1,τm), differentiating (41) with respect to μk yields that
(42)∂Zd(t∣μ)∂μk=∂Zd(τm-1∣μ)∂μk+∫τm-1t‍∂G(s,Zd(s∣μ),μm)∂Zd∂Zd(s∣μ)∂μkds,fork<m,∂Zd(t∣μ)∂μk=∂Zd(τm-1∣μ)∂μk+∫τm-1t‍∂G(s,Zd(s∣μ),μm)∂Zd∂Zd(s∣μ)∂μkds+∫τm-1t‍∂G(s,Zd(s∣μ),μm)∂uds,fork=m,∂Zd(t∣μ)∂μk=0,fork>m.
Here we have
(43)∂G(t,Zd(t∣μ),μm)∂Zd=(0I-M-1K-μmM-1R),∂G(t,Zd(t∣μ),μm)∂u=(000-M-1R)Zd(t∣μ),
for t∈[τm-1,τm). Define
(44)δkm={1,ifk=m,0,otherwise,δ^km={1,ifk≤m,0,otherwise.
Then, by (42) we obtain
(45)∂Zd(t∣μ)∂μk=δ^km∂Zd(τm-1∣μ)∂μk+δ^km∫τm-1t‍∂G(s,Zd(s∣μ),μm)∂Zd∂Zd(s∣μ)∂μkds+δkm∫τm-1t‍∂G(s,Zd(s∣μ),μm)∂uds,t∈[τm-1,τm).
Moreover, it is easy to see that
(46)∂Zd(0∣μ)∂μk=0.
For t∈[τm-1,τm), differentiating (45) with respect to t yields that
(47)ddt{∂Zd(t∣μ)∂μk}=δ^km∂G(t,Zd(t∣μ),μm)∂Zd∂Zd(t∣μ)∂μk+δkm∂G(t,Zd(t∣μ),μm)∂u.Now, we define(48)Γk(t∣μ)=∂Zd(t∣μ)∂μk,k=1,2,…,d.
Then, it follows from (46) and (47) that
(49)dΓk(t∣μ)dt=δ^km∂G(t,Zd(t∣μ),μm)∂ZdΓk(t∣μ)+δkm∂G(t,Zd(t∣μ),μm)∂u,t∈[τm-1,τm),m=1,2,…,d,
with the initial condition
(50)Γk(0∣μ)=0.
As a result, by using the chain rule, we can derive the gradient of Jdh(μ) with respect to μk, k=1,2,…,d, as follows:
(51)∂Jdh(μ)∂μk=∫0T‍(Zd(t∣μ),ΠΓk(t∣μ))R4Ndt+μk(τk-τk-1).
By incorporating these formulae into the SQP algorithm, we can solve the problem (OCPdh) numerically.


## 5. Numerical Simulations

In this section, we present some numerical simulation results by the approximation method presented in this paper. LetL=T=1, ρ(x)=q(x)=1, R1=0, and R2=1. Moreover, we take
(52)w0(x)=10x2(1-x)2,w1(x)=0,f(x,t)=(240-10π2x2(1-x)2)cos⁡(πt).
If the damping coefficient
(53)D(x,t)=0,(x,t)∈Q,
which implies that the system (1) is without any damping, the exact solution (see Figure 1) of (1) is
(54)w(x,t)=10cos⁡(πt)x2(1-x)2,(x,t)∈Q.
In this case, we have
(55)12∫0T‍∥(w(t)w˙(t))∥H2dt+12∫01‍|u(t)|2dt=20.39.Figure 1
The state without damping control; that is,u=0.In the following, we discuss the optimal control problem (OCP) with local Kelvin-Voigt damping acting on two different locations. We writeω1=[0.6,1] and ω2=[0.2,0.6]. By the approximation method presented in this paper, we can get the numerical results for the optimal control problem (OCP) with ω=ω1 and ω=ω2, which will be called problem (OCP)1 and problem (OCP)2, respectively. The numerical optimal controls for the two problems are presented by Figure 2. Moreover, the difference of w1 and w is shown in Figure 3, where w is defined by (54) and w1 is the numerical state function of (OCP)1. The difference of w2 and w has the same property, where w2 is the numerical state function of (OCP)2. For problem (OCP)1, the optimal value is 14.22 and for problem (OCP)2, the optimal value is 15.65, which means that the control effect of the location ω1 is better than the location ω2.Figure 2
The optimal controls for(OCP)1 and (OCP)2.Figure 3
The difference betweenw1 and w.



---

*Source: 982574-2014-06-23.xml*