# Numerical Exponential Decay to Dissipative Bresse System

Authors: M. L. Santos; Dilberto da S. Almeida Júnior
Journal: Journal of Applied Mathematics
 (2010)

---

## Abstract

We consider the Bresse system with frictional dissipative terms acting in all the equations. We show the exponential decay of the solution by using a method developed by Z. Liu and S. Zheng and their collaborators in past years. The numerical computations were made by using the finite difference method to prove the theoretical results. In particular, the finite difference method in our case is locking free.


---

## Body

## 1. Introduction and Main Results

The circular arch problem is also known as the Bresse system (see [1] for details). Elastic structures of the arches type are objects of study widely explored in engineering, architecture, marine engineering, aeronautics and others. In particular, the free vibrations about elastic structures is a function of their natural properties and are an important subject of investigation in engineering and also in mathematics. In the field of mathematical analysis is interesting to know properties which relate the behavior of the energy associated with the respective dynamic model. For feedback laws, for example, we can ask what conditions about dynamic model must be ensured to obtain the decay of the energy of solutions in the time t. In this sense, the property of stabilization has been studied for dynamic problems in elastic structures translated in terms of partial differential equations. A interesting property determines that the exponential decay with few feedback laws occurs only in an particular situation (see [2, 3]).Regarding the numerical aspects related to elastic structures, special attention was dedicated for a numerical pathology know as locking phenomenon about transverse shear force or simply shear locking (see [4–7]). The locking phenomenon is a numerical anomaly that affect some numerical methods when applied in elastic structures, such as finite element methods for linear shape functions. Basically, this obstacle numerical consists in an over-estimation about of the coefficient of rigidity and that obviously does not correspond to the real case. Therefore, the locking phenomenon is a numerical inconsistent. In this sense, any numerical method applied to elastic structures must be to avoid this numerical anomaly.Following the aim idea about deformation in elastic structures, we consider the Bresse system given by equations of motion(1.1)ρ1φtt=Qx+lN,ρ2ψtt=Mx-Q,ρ1wtt=Nx-lQ,
where we use N, Q and M to denote the axial force, the shear force and the bending moment, respectively. These forces are stress-strain relations for elastic behavior and given by(1.2)N=κ0(wx-lφ),Q=κ(φx+lw+ψ),M=bψx.
Here ρ1=ρA,ρ2=ρI,κ=k′GA,κ0=EA,b=EI,l=R-1 where ρ is the density of material, E is the modulus of elasticity, G is the shear modulus, k' is the shear factor, A is the cross-sectional area, I is the second moment of area of the cross-section and R is the radius of curvature. The functions w, φ, and ψ are the longitudinal, vertical and shear angle displacements, respectively (see Figure 1).Figure 1
The circular arch.From coupled equations (1.1) and (1.2) we obtain(1.3)ρ1φtt-κ(φx+ψ+lw)x-κ0l(wx-lφ)=0,in(0,L)×(0,T)ρ2ψtt-bψxx+κ(φx+ψ+lw)=0,in(0,L)×(0,T)ρ1wtt-κ0(wx-lφ)x+κl(φx+ψ+lw)=0,in(0,L)×(0,T).In the literature (1.3) are the equations for the theory of circular arch. For more details see [1].In this paper we will examine the issues concerning the asymptotic stabilization of Bresse system with frictional dissipative terms given by(1.4)ρ1φtt-κ(φx+ψ+lw)x-κ0l(wx-lφ)+γ1φt=0,in(0,L)×(0,T),(1.5)ρ2ψtt-bψxx+κ(φx+ψ+lw)+γ2ψt=0,in(0,L)×(0,T),(1.6)ρ1wtt-κ0(wx-lφ)x+κl(φx+ψ+lw)+γ3wt=0,in(0,L)×(0,T),
where γ1,γ2 and γ3 are positive constants. The initial conditions are(1.7)φ(·,0)=φ0,φt(·,0)=φ1,ψ(·,0)=ψ0,ψt(·,0)=ψ1,w(·,0)=w0,wt(·,0)=w1in]0,L[,
and we assume the Dirichlet boundary conditions(1.8)φ(0,t)=φ(L,t)=ψ(0,t)=ψ(L,t)=w(0,t)=w(L,t)=0in]0,T[.Remark 1.1.
IfR→∞, then l→0, and this model reduces to the well-known Timoshenko beam equations (see [1] for details).Related to the objectives of this paper there are few results in the literature. Recently, Liu and Rao [8] studied the asymptotic behavior of the circular arch in the context of linear thermoelasticity and they proved that the exponential decay occurs if and only if the speed of wave propagation that occur in the model are equal.Our main result is to prove that the Bresse system is exponentially stable in the presence of feedback laws and in the sequence we evidence our result by using the finite difference method of locking-free nature. This result was not considered for circular arch problem.The paper is organized as follows. In Section2 we establish the existence, regularity, and uniqueness of global solutions of (1.4)–(1.8). We use the semigroup techniques. In Section 3 we study the exponential decay of the strong solutions to system (1.4)–(1.8). We show the exponential decay of the solution by using a method developed by Liu and Zheng [9] and their collaborators in past years. Finally, in Section 4 some numerical aspects were considered. We use a particular discretization in finite difference method that is locking free (see [10, 11]) and our objective in this case is only to verify the numerical exponential decay for system (1.4)–(1.8). Questions about numerical analysis of numerical exponential decay as well the criterion stability for explicit time method in finite difference requires a more elaborate analysis.


## 2. The Semigroup Setting

In this section we will study the existence and uniqueness of strong and global solutions for the system (1.4)–(1.8) using the semigroup techniques. To give an accurate formulation of the evolution problem we are introducing the product Hilbert space(2.1)ℋ:=H01(0,L)×L2(0,L)×H01(0,L)×L2(0,L)×H01(0,L)×L2(0,L)
with norm given by(2.2)∥U∥ℋ2=∥(φ,φ̃,ψ,ψ̃,w,w̃)T∥ℋ2≡∫0Lρ1|φ̃|2+ρ2|ψ̃|2+ρ1|w̃|2+b|ψx|2+k|φx+ψ+lw|2+k0|wx-lφ|2dx.
Let U=(φ,φt,ψ,ψt,w,wt)T, and we define the operator 𝒜:D(𝒜)⊂ℋ→ℋ given by(2.3)𝒜=(010000kρ1∂x2-k0l2Iρ1-γ1ρ1Ikρ1∂x0(k+k0)lρ1∂x0000100-kρ2∂x0bρ2∂x2-kρ2I-γ2ρ2I-klρ2I0000001-(k0+k)lρ1∂x0lkIρ10k0ρ1∂x2-l2kρ1I-γ3ρ1I)
with  domain(2.4)D(𝒜):=(H01(0,L)∩H2(0,L)×H01(0,L))3.Therefore, system (1.4)–(1.8) is equivalent to(2.5)Ut=𝒜U,U(0)=U0,
where U0:=(φ0,φ1,ψ0,ψ1,w0,w1)T. It is not difficult to see that 𝒜 is a dissipative operator in the phase space ℋ. More precisely we have(2.6)ℛe(𝒜U,U)ℋ=-∫0L[γ1ρ1|φt|2+γ2ρ2|ψt|2+γ3ρ1|wt|2]dx≤0.
Under the above notation we can establish the following theorem.Theorem 2.1.
The operator𝒜 is the infinitesimal generator of C0-semigroup S(t) of contraction in ℋ.Proof.
Since𝒟(𝒜) is dense in ℋ and 𝒜 is a dissipative operator, to prove Theorem 2.1, it is sufficient to prove that 0 belongs to the resolvent set of 𝒜, that is, 0∈ρ(𝒜) (see [12, Theorem 4.3]). Let us take F=(f1,f2,f3,f4,f5,f6)∈ℋ, then there is an only U∈ℋ such that
(2.7)𝒜U=F.
It easily follows from (2.7) or the standard result on the linear elliptic equations that (2.7) has a unique solution (φ,ψ,w)∈(H01(Ω)∩H2(Ω))3. Therefore, 0∈ρ(𝒜) and (φ,φt,ψ,ψt,w,wt)∈𝒟(𝒜). The proof is now complete.


## 3. Exponential Decay

In this section we will prove that the semigroupS(t) on ℋ is exponentially stable. Here we will use necessary and sufficient conditions for C0-semigroups being exponentially stable in a Hilbert space. This result was obtained by Gearhart [13] and Huang [14], independently (see also [15]).Theorem 3.1.
LetS(t)=e𝒜t be a C0-semigroup of contractions on Hilbert space ℋ. Then S(t) is exponentially stable if and only if
(3.1)ρ(𝒜)⊇{iλ:λ∈ℝ}≡iℝ,(3.2)lim|λ|→∞¯∥(iλI-𝒜)-1∥ℒ(ℋ)<∞
hold, where ρ(𝒜) is the resolvent set of 𝒜.Theorem 3.2.
TheC0-semigroup of contractions (e𝒜t), t>0, generated by 𝒜, is exponentially stable.Proof.
To prove the exponential stability of(e𝒜t), it remains to verify the properties (3.1) and (3.2) of Theorem 3.1. First we prove that
(3.3)ρ(𝒜)⊇{iλ:λ∈ℝ}.
Suppose that the conclusion of (3.1) is not true. Then there is a β∈ℝ such that, β≠0, iβ∈σ(A) (spectrum of 𝒜) and since 𝒜-1 is compact, iβ must be an eigenvalue of 𝒜. Let U=(φ,φt,ψ,ψt,w,wt)T be with ∥U∥ℋ=1, such that 𝒜U=iβU. Using the definition of 𝒜 it follows that 𝒜U=iβU if and only if
(3.4)kρ1φxx-k0ρ1l2φ+kρ1ψx+(k+k0)ρ1lwx+β2φ=iγ1ρ1βφ,(3.5)bρ1ψxx-kρ2ψ-kρ2φx-kρ2lw+β2ψ=iγ2ρ2βψ,(3.6)k0ρ1wxx-kρ1l2w-(k+k0)ρ1lφx+kρ1lψ+β2w=iγ3ρ1w.
Multiplying equations (3.4) (3.5), and (3.6) by φ,ψ and w, respectively, and performing integration by parts on ]0,L[, we arrive at
(3.7)γ1β∫0Lφ2dx=0,γ2β∫0Lψ2dx=0,γ3β∫0Lw2dx=0,
which implies that φ=0,ψ=0,w=0 and consequently
(3.8)∫0Lb|ψx|2+k|φx+ψ+lw|2+k0|wx-lφ|2dx=0,
from where we can conclude that ∥U∥ℋ=0, which contradicts ∥U∥ℋ=1. This completes the proof of (3.1).
We now prove (3.2) by a contradiction argument again. Suppose that (3.2) is not true. Then there are a sequence Vn∈ℋ and a sequence βn∈ℝ such that(3.9)n∥Vn∥ℋ≤∥(IβnI-𝒜)-1Vn∥ℋ,∀n>0.
From iβn∈ρ(𝒜) we have that there exists a unique sequence Un∈𝒟(𝒜) such that
(3.10)iβnUn-𝒜Un=Vn,∥Un∥ℋ=1.
Denoting by Fn=iβnUn-𝒜Un it follows that
(3.11)∥Fn∥ℋ≤1n→0,asn→∞.
Let us denote by
(3.12)Un=(φn,φ̃n,ψn,ψ̃n,wn,w̃n)T,Fn=(f1n,f2n,f3n,f4n,f5n,f6n)T.
From Fn=iβnUn-𝒜Un we have the following equations in L2(0,L):
(3.13)iβnφn-φ̃n=f1n,(3.14)iρ1βnφ̃n-kφxxn+k0l2φn+γ1φ̃n-kψxn-(k+k0)lwxn=f2n,(3.15)iβnψn-ψ̃n=f3n,(3.16)iρ2βnψ̃n-bψxxn+kψn+kφxn+γ2ψ̃n+klwn=f4n,(3.17)iβnwn-w̃n=f5n,(3.18)iρ1βnw̃-k0wxxn+kl2wn+γ3w̃n+(k+k0)lφxn+lkψn=f6n.
Taking the inner product of (iβnI-𝒜)Un with Un in ℋ and then taking the real part, we obtain
(3.19)Re((iβnI-𝒜)Un,Un)ℋ=-∫0L[γ1ρ1|φ̃n|2+γ2ρ2|ψ̃n|2+γ3ρ1|w̃n|2]dx→0
from where we can conclude that
(3.20)φ̃n→0inL2(0,L),ψ̃n→0inL2(0,L),w̃n→0inL2(0,L).
Considering the equation
(3.21)βn∥Un∥2-i(𝒜Un,Un)ℋ=i(Fn,Un)ℋ,
then from (2.6)
(3.22)βn∥Un∥2+i∫0L[γ1ρ1|φ̃n|2+γ2ρ2|ψ̃n|2+γ3ρ1|w̃n|2]dx=i(Fn,Un)ℋ,
which implies that βn∥Un∥2→0, and then
(3.23)βn|φ̃n|2→0,βn|ψ̃n|2dx→0,βn|w̃n|2→0.
From
(3.24)φ̃n→0,ψ̃n→0,w̃n→0
we can conclude that
(3.25)βnφ̃n→0,βnψ̃n→0,βnw̃n→0.
From equations (3.14), (3.16) and (3.18) we obtain, after some calculations,
(3.26)∫0Lb|ψxn|2+k|φxn+ψn+lwn|2+k0|wxn-lφn|2dx→0.
In summary, we have proved that ∥Un∥ℋ→0 which contradicts ∥Un∥ℋ=1. Thus, the proof is complete.


## 4. Unlocking for Spatial Finite Difference Scheme

For numerical verification of the exponential decay of system (1.4)–(1.8) we use the total discretization in finite difference method. With respect to numerical schemes in this case, some aspects should be taking into account. The first one concerns to numerical phenomenon known as shear locking, which affects some numerical models applied to vibration problems in structures as shell, plates and beams. For the locking phenomenon (see [4–7]) for ample discussion for plane beams. For fast discussion about this numerical problem, we considerer the case of plane beams described by theory of Timoshenko and given by the following equations:(4.1)ρ1φtt-κ(φx+ψ)x=0,(4.2)ρ2ψtt-bψxx+κ(φx+ψ)=0.
are resulting of the equations (1.4)–(1.5) for l→0.By using the finite element method standard with linear shape functions, the rigidity coefficientb=EI is modified for(4.3)b*=b(1+κ12bh2)=EI[1+k′GE(hϵ)2],
where b=EI,κ=kGA,A=aϵ,I=aϵ3/12 considering a rectangular geometry with width a and thickness ϵ. In particular, for plane beams must be have ϵ<h where h is the spatial division (see ([10, 11] for details). Consequently, b*>b and clearly this value for rigidity coefficient does not correspond to the real case for h fixed. In reality, with b* the equation for angle rotation in (4.2) is rewrite as(4.4)ρ2ψtt-bψxx-κ12h2ψxx+κ(φx+ψ)=0
and the energy respective is(4.5)E(t)=12∫0L[ρ1|φt|2+ρ2|ψt|2]dx+12∫0L[κ|φx+ψ|2+b*|ψx|2]dx
and obviously this energy is different of the real case with b* = b. In general the numerical schemes in finite element method or finite difference method are locking free/over-estimation if not exist any additional term about the coefficients of the system, because the shear locking problem is basically a numerical anomaly characterized by a over-estimation about rigidity coefficients and dependent of parameter h.To avoid this numerical anomaly the finite element method can be used, however, special care must be taking in to account for choose of basis functions [6, 7]. Naturally, this shear locking/over-estimation can be to affect equations (1.4)–(1.6). See the studies by Loula et al. in [16–18] for numerical treatment in finite element for circular arch problem.In our case, we use the total discretization in finite difference method and to avoid shear locking/over-estimation we make a particular discretization for the functions of zero derivative such asκ0l2φ,κψ and κlw for a numerical operator of second order in relation to Δx.For our purposes, we use the space-time explicit method applied to equations (1.4)–(1.6) and we define Δx=L/(J+1), Δt=T/(N+1) for J,N∈ℕ and nets(4.6)x0=0<x1=Δx<⋯<xJ=JΔx<xJ+1=L,t0=0<t1=Δt<⋯<tN=NΔt<tN+1=T,
where xj=jΔx and tn=nΔt for j=0,1,2,…,J+1 and n=0,1,2,…,N+1. The numerical scheme consists in to find φjn+1,ψjn+1 and wjn+1 such that(4.7)ρ1∂¯t∂tφjn=κ∂¯x∂xφjn+κ∂x+∂¯x2ψjn+κ¯∂x+∂¯x2wjn-κ0l22(φj-1/2n+φj+1/2n)-γ1∂t+∂¯t2φjn,(4.8)ρ2∂¯t∂tψjn=b∂¯x∂xψjn-κ∂x+∂¯x2φjn-κ2(ψj-1/2n+ψj+1/2n)-κl2(wj-1/2n+wj+1/2n)-γ2∂t+∂¯t2ψjn,(4.9)ρ1∂¯t∂twjn=κ0∂¯x∂xwjn-κ¯∂x+∂¯x2φjn-κl2(ψj-1/2n+ψj+1/2n)-κl22(wj-1/2n+wj+1/2n)-γ3∂t+∂¯t2wjn
for j=1,2,…,J,n=1,2,…,N, with the following numerical operators of second order for a function u(x,t):(4.10)∂x+∂¯x2ujn=uj+1n-uj-1n2Δx,∂t+∂¯t2ujn=ujn+1-ujn-12Δt,∂¯x∂xujn=uj+1n-2ujn+uj-1nΔx2,∂¯t∂tujn=ujn+1-2ujn+ujn-1Δt2.Foruj-1/2n and uj+1/2n we denote the average of u(xj,tn) on the points (xj-1,tn),(xj,tn) and (xj+1,tn),(xj,tn), respectively. This approximation avoid any over-estimation about the coefficients of equations (1.4)–(1.6). Then, we have,(4.11)uj-1/2n+uj+1/2n2:=uj+1n+2ujn+uj-1n4.The boundary conditions in the numerical context are given by(4.12)φ0n=φJ+1n=0,ψ0n=ψJ+1n=0,w0n=wJ+1n=0,∀n=1,2,…,N,
and the discretizations to initial conditions are given by(4.13)φj0=φ(xj,0),φj1=φj-1+2Δtφt(xj,0),∀j=1,2,…,J,ψj0=ψ(xj,0),ψj1=ψj-1+2Δtψt(xj,0),∀j=1,2,…,J,wj0=w(xj,0),wj1=wj-1+2Δtwt(xj,0),∀j=1,2,…,J.The numerical energy associated for (4.7)–(4.13) is given by(4.14)En:=Δx2∑j=0J[ρ1(φjn+1-φjnΔt)2+ρ2(ψjn+1-ψjnΔt)2+ρ1(wjn+1-wjnΔt)2+b(ψj+1n+1-ψjn+1Δx)(ψj+1n-ψjnΔx)+κ0(wj+1n+1-wjn+1Δx-lφj+1n+1+φjn+12)(wj+1n-wjnΔx-lφj+1n+φjn2)+κ(φj+1n-φjnΔx+ψj+1n+ψjn2+lwj+1n+wjn2)(φjn+1-φjnΔt)2].Equations (4.7)–(4.9), in fact, are locking free because the numerical energy (4.14) is compatible with the continuous energy in (2.2), because the coefficients b,κ0 and κ are exactly those in (2.2), without any dependence with h. To verify this, we have the following proposition.Proposition 4.1.
For allΔt,Δx∈(0,L) the energy (4.14) of solutions of the discrete equations (4.7)–(4.9), with initial conditions (4.13) and any boundary conditions (4.12) is such that
(4.15)En≤E0,∀n=1,2,…,N.Proof.
The proof is more extensive and we have omitted it here. For a idea of the proof, we use the multiplicative techniques such as performed in [19, 20], that is we multiply the equations (4.7)–(4.9) by (1/2)(φjn+1-φjn-1),(1/2)(ψjn+1-ψjn-1) and (1/2)(wjn+1-wjn-1), respectively, and after we applied the discrete sum for j=1,…,J. Then after some properly simplifications and taking in to account the boundary conditions given in (4.12) we get (4.14).It is clear that discrete equations (4.7)–(4.9) are all consistent with 𝒪(Δx2,Δt2) for truncation error. Then, by Lax equivalence theorem the equations (4.7)–(4.9) are convergent if, and only if, they are stable. But another numerical limitation occurs for explicit time methods when applied to vibrations problem in elastic structures, with, say, respect to numerical stability. In particular, for explicit time method in finite difference applied to vibrations in Timoshenko beams, the prevailing numerical stability is given by Δt≤ϵ/3cs where cs=k′G/ρ. The limitation occurs when ϵ is small. The understanding and overcoming for this limitation was studied by Joseph P. Wright in [10, 11]. Naturally, the same numerical problem affect (4.7)–(4.9), in function of the dependence of thickness ϵ. Then, this problem should be the objective of study in another opportunity, and for our purpose of numerical verification to exponential decay, we consider ϵ fixed.For numerical example, we considerL=2πm, thickness ϵ=0.025m, width a=0.0040m, l=1/20,E=21×104N/m2,ρ=7850kg/m3,k′=5/6,r=0.29(Poissonratio), G=E/(2+2r) and the following initial conditions:(4.16)φ(xj,0)=ψ(xj,0)=w(xj,0)=0,wt(xj,0)=sin(μπxjL),φt(xj,0)=ψt(xj,0)=0,μ∈ℕ.First we reproduce the conservative case,γi=0,i=1,2,3 and for dissipative case, we consider γi=π,i=1,2,3. For the computational domain chosen 32 spatial points and 312 points in the time domain were given by T=1.5 seconds (see Figures 1, 2, 3, 4, 5, 6, and 7).Figure 2
Conservative case,μ=3.Figure 3
Conservative case,μ=3.Figure 4
Total dissipation,μ=3.Figure 5
Total dissipation,μ=3.Figure 6
Total dissipation,μ=5.Figure 7
Total dissipation,μ=5.With Figures8, 9, 10, 11, 12, and 13 we illustrate the simulation results for vibration in angle rotation ψ and the behavior of numerical energy in the cases conservative and dissipative. Of course the numerical energy (4.14) must be conservative in the absence of damping.Figure 8
Conservative case,μ=10.Figure 9
Conservative case,μ=10.Figure 10
Total dissipation case,μ=10.Figure 11
Total dissipation,μ=10.Figure 12
Conservative case.Figure 13
Dissipative case.



---

*Source: 848620-2010-08-09.xml*