Numerical Exponential Decay to Dissipative Bresse System

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.

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 𝑡. 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) where we use 𝑁, 𝑄 and 𝑀 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) Here 𝜌₁=𝜌𝐴,𝜌₂=𝜌𝐼,𝜅=𝑘^′𝐺𝐴,𝜅₀=𝐸𝐴,𝑏=𝐸𝐼,𝑙=𝑅⁻¹ where 𝜌 is the density of material, 𝐸 is the modulus of elasticity, 𝐺 is the shear modulus, 𝑘^′ is the shear factor, 𝐴 is the cross-sectional area, 𝐼 is the second moment of area of the cross-section and 𝑅 is the radius of curvature. The functions 𝑤, 𝜑, 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𝜌₁𝜑_𝑡𝑡−𝜅(𝜑_𝑥+𝜓+𝑙𝑤)_𝑥−𝜅₀𝑙(𝑤_𝑥𝜌−𝑙𝜑)=0,in(0,𝐿)×(0,𝑇)₂𝜓_𝑡𝑡−𝑏𝜓_𝑥𝑥+𝜅(𝜑_𝑥𝜌+𝜓+𝑙𝑤)=0,in(0,𝐿)×(0,𝑇)₁𝑤_𝑡𝑡−𝜅₀(𝑤_𝑥−𝑙𝜑)_𝑥+𝜅𝑙(𝜑_𝑥+𝜓+𝑙𝑤)=0,in(0,𝐿)×(0,𝑇). (1.3)

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𝜌₁𝜑_𝑡𝑡−𝜅(𝜑_𝑥+𝜓+𝑙𝑤)_𝑥−𝜅₀𝑙(𝑤_𝑥−𝑙𝜑)+𝛾₁𝜑_𝑡𝜌=0,in(0,𝐿)×(0,𝑇),₂𝜓_𝑡𝑡−𝑏𝜓_𝑥𝑥+𝜅(𝜑_𝑥+𝜓+𝑙𝑤)+𝛾₂𝜓_𝑡𝜌=0,in(0,𝐿)×(0,𝑇),₁𝑤_𝑡𝑡−𝜅₀(𝑤_𝑥−𝑙𝜑)_𝑥+𝜅𝑙(𝜑_𝑥+𝜓+𝑙𝑤)+𝛾₃𝑤_𝑡=0,in(0,𝐿)×(0,𝑇), (1.4)(1.5)(1.6) where 𝛾₁,𝛾₂ and 𝛾₃ are positive constants. The initial conditions are𝜑(⋅,0)=𝜑₀,𝜑_𝑡(⋅,0)=𝜑₁,𝜓(⋅,0)=𝜓₀,𝜓_𝑡(⋅,0)=𝜓₁,𝑤(⋅,0)=𝑤₀,𝑤_𝑡(⋅,0)=𝑤₁in]0,𝐿[, (1.7) and we assume the Dirichlet boundary conditions𝜑(0,𝑡)=𝜑(𝐿,𝑡)=𝜓(0,𝑡)=𝜓(𝐿,𝑡)=𝑤(0,𝑡)=𝑤(𝐿,𝑡)=0in]0,𝑇[. (1.8)

Remark 1.1. If 𝑅→∞, then 𝑙→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 Section 2 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ℋ∶=𝐻¹₀(0,𝐿)×𝐿²(0,𝐿)×𝐻¹₀(0,𝐿)×𝐿²(0,𝐿)×𝐻¹₀(0,𝐿)×𝐿²(0,𝐿) (2.1) with norm given by‖𝑈‖²_ℋ=‖(𝜑,𝜑,𝜓,𝜓,𝑤,𝑤)^𝑇‖²_ℋ≡^𝐿₀𝜌₁|𝜑|²+𝜌₂|𝜓|²+𝜌₁|𝑤|²+𝑏|𝜓_𝑥|²+𝑘|𝜑_𝑥+𝜓+𝑙𝑤|²+𝑘₀|𝑤_𝑥−𝑙𝜑|²𝑑𝑥. (2.2) Let 𝑈=(𝜑,𝜑_𝑡,𝜓,𝜓_𝑡,𝑤,𝑤_𝑡)^𝑇, and we define the operator 𝒜∶𝐷(𝒜)⊂ℋ→ℋ given by𝑘𝒜=010000𝜌₁𝜕²_𝑥−𝑘₀𝑙²𝐼𝜌₁−𝛾₁𝜌₁𝐼𝑘𝜌₁𝜕_𝑥0(𝑘+𝑘₀)𝑙𝜌₁𝜕_𝑥0−𝑘000100𝜌₂𝜕_𝑥0𝑏𝜌₂𝜕²_𝑥−𝑘𝜌₂𝐼−𝛾₂𝜌₂𝐼−𝑘𝑙𝜌₂𝐼0−000001(𝑘₀+𝑘)𝑙𝜌₁𝜕_𝑥0𝑙𝑘𝐼𝜌₁0𝑘₀𝜌₁𝜕²_𝑥−𝑙²𝑘𝜌₁𝐼−𝛾₃𝜌₁𝐼 (2.3) with domain𝐷(𝒜)∶=(𝐻¹₀(0,𝐿)∩𝐻²(0,𝐿)×𝐻¹₀(0,𝐿))³. (2.4)

Therefore, system (1.4)–(1.8) is equivalent to𝑈_𝑡=𝒜𝑈,𝑈(0)=𝑈₀, (2.5) where 𝑈₀∶=(𝜑₀,𝜑₁,𝜓₀,𝜓₁,𝑤₀,𝑤₁)^𝑇. It is not difficult to see that 𝒜 is a dissipative operator in the phase space ℋ. More precisely we haveℛ𝑒(𝒜𝑈,𝑈)_ℋ=−^𝐿₀[𝛾₁𝜌₁|𝜑_𝑡|²+𝛾₂𝜌₂|𝜓_𝑡|²+𝛾₃𝜌₁|𝑤_𝑡|²]𝑑𝑥≤0. (2.6) Under the above notation we can establish the following theorem.

Theorem 2.1. The operator 𝒜 is the infinitesimal generator of 𝐶₀-semigroup 𝑆(𝑡) 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 𝐹=(𝑓₁,𝑓₂,𝑓₃,𝑓₄,𝑓₅,𝑓₆)∈ℋ, then there is an only 𝑈∈ℋ such that 𝒜𝑈=𝐹. (2.7) It easily follows from (2.7) or the standard result on the linear elliptic equations that (2.7) has a unique solution (𝜑,𝜓,𝑤)∈(𝐻¹₀(Ω)∩𝐻²(Ω))³. Therefore, 0∈𝜌(𝒜) and (𝜑,𝜑_𝑡,𝜓,𝜓_𝑡,𝑤,𝑤_𝑡)∈𝒟(𝒜). The proof is now complete.

3. Exponential Decay

In this section we will prove that the semigroup 𝑆(𝑡) on ℋ is exponentially stable. Here we will use necessary and sufficient conditions for 𝐶₀-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. Let 𝑆(𝑡)=𝑒^𝒜𝑡 be a 𝐶₀-semigroup of contractions on Hilbert space ℋ. Then 𝑆(𝑡) is exponentially stable if and only if 𝜌(𝒜)⊇{𝑖𝜆∶𝜆∈ℝ}≡𝑖ℝ,lim_|𝜆|→∞‖(𝑖𝜆𝐼−𝒜)⁻¹‖_ℒ(ℋ)<∞ (3.1)(3.2) hold, where 𝜌(𝒜) is the resolvent set of 𝒜.

Theorem 3.2. The 𝐶₀-semigroup of contractions (𝑒^𝒜𝑡), 𝑡>0, generated by 𝒜, is exponentially stable.

Proof. To prove the exponential stability of (𝑒^𝒜𝑡), it remains to verify the properties (3.1) and (3.2) of Theorem 3.1. First we prove that 𝜌(𝒜)⊇{𝑖𝜆∶𝜆∈ℝ}. (3.3) Suppose that the conclusion of (3.1) is not true. Then there is a 𝛽∈ℝ such that, 𝛽≠0, 𝑖𝛽∈𝜎(𝐴) (spectrum of 𝒜) and since 𝒜⁻¹ is compact, 𝑖𝛽 must be an eigenvalue of 𝒜. Let 𝑈=(𝜑,𝜑_𝑡,𝜓,𝜓_𝑡,𝑤,𝑤_𝑡)^𝑇 be with ‖𝑈‖_ℋ=1, such that 𝒜𝑈=𝑖𝛽𝑈. Using the definition of 𝒜 it follows that 𝒜𝑈=𝑖𝛽𝑈 if and only if 𝑘𝜌₁𝜑_𝑥𝑥−𝑘₀𝜌₁𝑙²𝜑+𝑘𝜌₁𝜓_𝑥+(𝑘+𝑘₀)𝜌₁𝑙𝑤_𝑥+𝛽²𝜑=𝑖𝛾₁𝜌₁,𝑏𝛽𝜑𝜌₁𝜓_𝑥𝑥−𝑘𝜌₂𝜓−𝑘𝜌₂𝜑_𝑥−𝑘𝜌₂𝑙𝑤+𝛽²𝜓=𝑖𝛾₂𝜌₂𝑘𝛽𝜓₀,𝜌₁𝑤_𝑥𝑥−𝑘𝜌₁𝑙²𝑤−(𝑘+𝑘₀)𝜌₁𝑙𝜑_𝑥+𝑘𝜌₁𝑙𝜓+𝛽²𝑤=𝑖𝛾₃𝜌₁𝑤. (3.4)(3.5)(3.6) Multiplying equations (3.4) (3.5), and (3.6) by 𝜑,𝜓 and 𝑤, respectively, and performing integration by parts on ]0,𝐿[, we arrive at 𝛾₁𝛽^𝐿₀𝜑²𝛾𝑑𝑥=0,₂𝛽^𝐿₀𝜓²𝛾𝑑𝑥=0,₃𝛽^𝐿₀𝑤²𝑑𝑥=0, (3.7) which implies that 𝜑=0,𝜓=0,𝑤=0 and consequently ^𝐿₀𝑏|𝜓_𝑥|²+𝑘|𝜑_𝑥+𝜓+𝑙𝑤|²+𝑘₀|𝑤_𝑥−𝑙𝜑|²𝑑𝑥=0, (3.8) from where we can conclude that ‖𝑈‖_ℋ=0, which contradicts ‖𝑈‖_ℋ=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 𝑉_𝑛∈ℋ and a sequence 𝛽_𝑛∈ℝ such that𝑛‖𝑉_𝑛‖_ℋ≤‖(𝐼𝛽_𝑛𝐼−𝒜)⁻¹𝑉_𝑛‖_ℋ,∀𝑛>0. (3.9) From 𝑖𝛽_𝑛∈𝜌(𝒜) we have that there exists a unique sequence 𝑈_𝑛∈𝒟(𝒜) such that 𝑖𝛽_𝑛𝑈_𝑛−𝒜𝑈_𝑛=𝑉_𝑛,‖𝑈_𝑛‖_ℋ=1. (3.10) Denoting by 𝐹_𝑛=𝑖𝛽_𝑛𝑈_𝑛−𝒜𝑈_𝑛 it follows that ‖𝐹_𝑛‖_ℋ≤1𝑛→0,as𝑛→∞. (3.11) Let us denote by 𝑈_𝑛=(𝜑^𝑛,𝜑^𝑛,𝜓^𝑛,𝜓^𝑛,𝑤^𝑛,𝑤^𝑛)^𝑇,𝐹_𝑛=(𝑓^𝑛₁,𝑓^𝑛₂,𝑓^𝑛₃,𝑓^𝑛₄,𝑓^𝑛₅,𝑓^𝑛₆)^𝑇. (3.12) From 𝐹_𝑛=𝑖𝛽_𝑛𝑈_𝑛−𝒜𝑈_𝑛 we have the following equations in 𝐿²(0,𝐿): 𝑖𝛽_𝑛𝜑^𝑛−𝜑^𝑛=𝑓^𝑛₁,𝑖𝜌₁𝛽_𝑛𝜑^𝑛−𝑘𝜑^𝑛_𝑥𝑥+𝑘₀𝑙²𝜑^𝑛+𝛾₁𝜑^𝑛−𝑘𝜓^𝑛_𝑥−(𝑘+𝑘₀)𝑙𝑤^𝑛_𝑥=𝑓^𝑛₂,𝑖𝛽_𝑛𝜓^𝑛−𝜓^𝑛=𝑓^𝑛₃,𝑖𝜌₂𝛽_𝑛𝜓^𝑛−𝑏𝜓^𝑛_𝑥𝑥+𝑘𝜓^𝑛+𝑘𝜑^𝑛_𝑥+𝛾₂𝜓^𝑛+𝑘𝑙𝑤^𝑛=𝑓^𝑛₄,𝑖𝛽_𝑛𝑤^𝑛−𝑤^𝑛=𝑓^𝑛₅,𝑖𝜌₁𝛽_𝑛𝑤−𝑘₀𝑤^𝑛_𝑥𝑥+𝑘𝑙²𝑤^𝑛+𝛾₃𝑤^𝑛+(𝑘+𝑘₀)𝑙𝜑^𝑛_𝑥+𝑙𝑘𝜓^𝑛=𝑓^𝑛₆. (3.13)(3.14)(3.15)(3.16)(3.17)(3.18) Taking the inner product of (𝑖𝛽_𝑛𝐼−𝒜)𝑈^𝑛 with 𝑈^𝑛 in ℋ and then taking the real part, we obtain 𝑅𝑒((𝑖𝛽_𝑛𝐼−𝒜)𝑈^𝑛,𝑈^𝑛)_ℋ=−^𝐿₀[𝛾₁𝜌₁|𝜑^𝑛|²+𝛾₂𝜌₂|𝜓^𝑛|²+𝛾₃𝜌₁|𝑤^𝑛|²]𝑑𝑥→0 (3.19) from where we can conclude that 𝜑^𝑛→0in𝐿²𝜓(0,𝐿),_𝑛→0in𝐿²𝑤(0,𝐿),_𝑛→0in𝐿²(0,𝐿). (3.20) Considering the equation 𝛽_𝑛‖𝑈^𝑛‖²−𝑖(𝒜𝑈^𝑛,𝑈^𝑛)_ℋ=𝑖(𝐹^𝑛,𝑈^𝑛)_ℋ, (3.21) then from (2.6) 𝛽_𝑛‖𝑈^𝑛‖²+𝑖^𝐿₀[𝛾₁𝜌₁|𝜑^𝑛|²+𝛾₂𝜌₂|𝜓^𝑛|²+𝛾₃𝜌₁|𝑤^𝑛|²]𝑑𝑥=𝑖(𝐹^𝑛,𝑈^𝑛)_ℋ, (3.22) which implies that 𝛽_𝑛‖𝑈^𝑛‖²→0, and then 𝛽_𝑛|𝜑^𝑛|²→0,𝛽_𝑛|𝜓^𝑛|²𝑑𝑥→0,𝛽_𝑛|𝑤^𝑛|²→0. (3.23) From 𝜑^𝑛→0,𝜓^𝑛→0,𝑤^𝑛→0 (3.24) we can conclude that 𝛽_𝑛𝜑^𝑛→0,𝛽_𝑛𝜓^𝑛→0,𝛽_𝑛𝑤^𝑛→0. (3.25) From equations (3.14), (3.16) and (3.18) we obtain, after some calculations, ^𝐿₀𝑏|𝜓^𝑛_𝑥|²+𝑘|𝜑^𝑛_𝑥+𝜓^𝑛+𝑙𝑤^𝑛|²+𝑘₀|𝑤^𝑛_𝑥−𝑙𝜑^𝑛|²𝑑𝑥→0. (3.26) In summary, we have proved that ‖𝑈^𝑛‖_ℋ→0 which contradicts ‖𝑈^𝑛‖_ℋ=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:𝜌₁𝜑_𝑡𝑡−𝜅(𝜑_𝑥+𝜓)_𝑥𝜌=0,₂𝜓_𝑡𝑡−𝑏𝜓_𝑥𝑥+𝜅(𝜑_𝑥+𝜓)=0. (4.1)(4.2) are resulting of the equations (1.4)–(1.5) for 𝑙→0.

By using the finite element method standard with linear shape functions, the rigidity coefficient 𝑏=𝐸𝐼 is modified for𝐛^∗𝜅=𝑏(1+ℎ12𝑏²𝑘)=𝐸𝐼[1+^′𝐺𝐸(ℎ𝜖)²], (4.3) where 𝑏=𝐸𝐼,𝜅=𝑘𝐺𝐴,𝐴=𝑎𝜖,𝐼=𝑎𝜖³/12 considering a rectangular geometry with width a and thickness 𝜖. In particular, for plane beams must be have 𝜖<ℎ where ℎ is the spatial division (see ([10, 11] for details). Consequently, 𝐛^∗>𝑏 and clearly this value for rigidity coefficient does not correspond to the real case for ℎ fixed. In reality, with 𝐛^∗ the equation for angle rotation in (4.2) is rewrite as𝜌₂𝜓_𝑡𝑡−𝑏𝜓_𝑥𝑥−𝜅ℎ12²𝜓_𝑥𝑥+𝜅(𝜑_𝑥+𝜓)=0 (4.4) and the energy respective is1𝐸(𝑡)=2^𝐿₀[𝜌₁|𝜑_𝑡|²+𝜌₂|𝜓_𝑡|²1]𝑑𝑥+2^𝐿₀[𝜅|𝜑_𝑥+𝜓|²+𝐛^∗|𝜓_𝑥|²]𝑑𝑥 (4.5) and obviously this energy is different of the real case with 𝐛^∗ = 𝑏. 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 ℎ.

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 𝜅₀𝑙²𝜑,𝜅𝜓 and 𝜅𝑙𝑤 for a numerical operator of second order in relation to Δ𝑥.

For our purposes, we use the space-time explicit method applied to equations (1.4)–(1.6) and we define Δ𝑥=𝐿/(𝐽+1), Δ𝑡=𝑇/(𝑁+1) for 𝐽,𝑁∈ℕ and nets𝑥₀=0<𝑥₁=Δ𝑥<⋯<𝑥_𝐽=𝐽Δ𝑥<𝑥_𝐽+1𝑡=𝐿,₀=0<𝑡₁=Δ𝑡<⋯<𝑡_𝑁=𝑁Δ𝑡<𝑡_𝑁+1=𝑇, (4.6) where 𝑥_𝑗=𝑗Δ𝑥 and 𝑡_𝑛=𝑛Δ𝑡 for 𝑗=0,1,2,…,𝐽+1 and 𝑛=0,1,2,…,𝑁+1. The numerical scheme consists in to find 𝜑_𝑗^𝑛+1,𝜓_𝑗^𝑛+1 and 𝑤_𝑗^𝑛+1 such that𝜌₁𝜕_𝑡𝜕_𝑡𝜑^𝑛_𝑗=𝜅𝜕_𝑥𝜕_𝑥𝜑^𝑛_𝑗𝜕+𝜅^𝑥+𝜕_𝑥2𝜓^𝑛_𝑗+𝜅𝜕_𝑥+𝜕_𝑥2𝑤^𝑛_𝑗−𝜅₀𝑙²2(𝜑^𝑛_𝑗−1/2+𝜑^𝑛_𝑗+1/2)−𝛾₁𝜕_𝑡+𝜕_𝑡2𝜑^𝑛_𝑗𝜌₂,𝜕_𝑡𝜕_𝑡𝜓^𝑛_𝑗=𝑏𝜕_𝑥𝜕_𝑥𝜓^𝑛_𝑗𝜕−𝜅^𝑥+𝜕_𝑥2𝜑^𝑛_𝑗−𝜅2(𝜓^𝑛_𝑗−1/2+𝜓^𝑛_𝑗+1/2)−𝜅𝑙2(𝑤^𝑛_𝑗−1/2+𝑤^𝑛_𝑗+1/2)−𝛾₂𝜕_𝑡+𝜕_𝑡2𝜓^𝑛_𝑗𝜌₁,𝜕_𝑡𝜕_𝑡𝑤^𝑛_𝑗=𝜅₀𝜕_𝑥𝜕_𝑥𝑤^𝑛_𝑗−𝜅𝜕_𝑥+𝜕_𝑥2𝜑^𝑛_𝑗−𝜅𝑙2(𝜓^𝑛_𝑗−1/2+𝜓^𝑛_𝑗+1/2)−𝜅𝑙²2(𝑤^𝑛_𝑗−1/2+𝑤^𝑛_𝑗+1/2)−𝛾₃𝜕_𝑡+𝜕_𝑡2𝑤^𝑛_𝑗 (4.7)(4.8)(4.9) for 𝑗=1,2,…,𝐽,𝑛=1,2,…,𝑁, with the following numerical operators of second order for a function 𝑢(𝑥,𝑡):𝜕_𝑥+𝜕_𝑥2𝑢^𝑛_𝑗=𝑢^𝑛_𝑗+1−𝑢^𝑛_𝑗−1,𝜕2Δ𝑥^𝑡+𝜕_𝑡2𝑢^𝑛_𝑗=𝑢_𝑗^𝑛+1−𝑢_𝑗^𝑛−1,2Δ𝑡𝜕_𝑥𝜕_𝑥𝑢^𝑛_𝑗=𝑢^𝑛_𝑗+1−2𝑢^𝑛_𝑗+𝑢^𝑛_𝑗−1Δ𝑥²,𝜕_𝑡𝜕_𝑡𝑢^𝑛_𝑗=𝑢_𝑗^𝑛+1−2𝑢^𝑛_𝑗+𝑢_𝑗^𝑛−1Δ𝑡². (4.10)

For 𝑢^𝑛_𝑗−1/2 and 𝑢^𝑛_𝑗+1/2 we denote the average of 𝑢(𝑥_𝑗,𝑡_𝑛) on the points (𝑥_𝑗−1,𝑡_𝑛),(𝑥_𝑗,𝑡_𝑛) and (𝑥_𝑗+1,𝑡_𝑛),(𝑥_𝑗,𝑡_𝑛), respectively. This approximation avoid any over-estimation about the coefficients of equations (1.4)–(1.6). Then, we have,𝑢^𝑛_𝑗−1/2+𝑢^𝑛_𝑗+1/22𝑢∶=^𝑛𝑗+1+2𝑢^𝑛_𝑗+𝑢^𝑛_𝑗−14. (4.11)

The boundary conditions in the numerical context are given by𝜑^𝑛₀=𝜑^𝑛_𝐽+1=0,𝜓^𝑛₀=𝜓^𝑛_𝐽+1=0,𝑤^𝑛₀=𝑤^𝑛_𝐽+1=0,∀𝑛=1,2,…,𝑁, (4.12) and the discretizations to initial conditions are given by𝜑⁰_𝑗=𝜑(𝑥_𝑗,0),𝜑¹_𝑗=𝜑_𝑗⁻¹+2Δ𝑡𝜑_𝑡(𝑥_𝑗𝜓,0),∀𝑗=1,2,…,𝐽,_0𝑗=𝜓(𝑥_𝑗,0),𝜓¹_𝑗=𝜓_𝑗⁻¹+2Δ𝑡𝜓_𝑡(𝑥_𝑗𝑤,0),∀𝑗=1,2,…,𝐽,_0𝑗=𝑤(𝑥_𝑗,0),𝑤¹_𝑗=𝑤_𝑗⁻¹+2Δ𝑡𝑤_𝑡(𝑥_𝑗,0),∀𝑗=1,2,…,𝐽. (4.13)

The numerical energy associated for (4.7)–(4.13) is given by𝐸^𝑛∶=Δ𝑥2^𝐽_𝑗=0[𝜌₁(𝜑_𝑗^𝑛+1−𝜑^𝑛_𝑗)Δ𝑡²+𝜌₂(𝜓_𝑗^𝑛+1−𝜓^𝑛_𝑗)Δ𝑡²+𝜌₁(𝑤_𝑗^𝑛+1−𝑤^𝑛_𝑗)Δ𝑡²𝜓+𝑏(^𝑛+1𝑗+1−𝜓_𝑗^𝑛+1𝜓Δ𝑥)(^𝑛𝑗+1−𝜓^𝑛_𝑗)Δ𝑥+𝜅₀(𝑤^𝑛+1_𝑗+1−𝑤_𝑗^𝑛+1𝜑Δ𝑥−𝑙^𝑛+1𝑗+1+𝜑_𝑗^𝑛+12𝑤)(^𝑛𝑗+1−𝑤^𝑛_𝑗𝜑Δ𝑥−𝑙^𝑛𝑗+1+𝜑^𝑛_𝑗2)𝜑+𝜅(^𝑛𝑗+1−𝜑^𝑛_𝑗+𝜓Δ𝑥^𝑛𝑗+1+𝜓^𝑛_𝑗2𝑤+𝑙^𝑛𝑗+1+𝑤^𝑛_𝑗2)]. (4.14)

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 𝑏,𝜅₀ and 𝜅 are exactly those in (2.2), without any dependence with ℎ. To verify this, we have the following proposition.

Proposition 4.1. For all Δ𝑡,Δ𝑥∈(0,𝐿) 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 𝐸^𝑛≤𝐸⁰,∀𝑛=1,2,…,𝑁. (4.15)

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)(𝜑_𝑗^𝑛+1−𝜑_𝑗^𝑛−1),(1/2)(𝜓_𝑗^𝑛+1−𝜓_𝑗^𝑛−1) and (1/2)(𝑤_𝑗^𝑛+1−𝑤_𝑗^𝑛−1), respectively, and after we applied the discrete sum for 𝑗=1,…,𝐽. 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 𝒪(Δ𝑥²,Δ𝑡²) 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 Δ𝑡≤𝜖/√3𝑐_𝑠 where 𝑐_𝑠=𝑘^′𝐺/𝜌. 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 consider 𝐿=2𝜋m, thickness 𝜖=0.025m, width 𝑎=0.0040m, 𝑙=1/20,𝐸=21×10⁴N/m²,𝜌=7850kg/m³,𝑘^′=5/6,𝑟=0.29(Poissonratio), 𝐺=𝐸/(2+2𝑟) and the following initial conditions:𝜑(𝑥_𝑗,0)=𝜓(𝑥_𝑗,0)=𝑤(𝑥_𝑗𝑤,0)=0,_𝑡(𝑥_𝑗,0)=sin(𝜇𝜋𝑥_𝑗𝐿),𝜑_𝑡(𝑥_𝑗,0)=𝜓_𝑡(𝑥_𝑗,0)=0,𝜇∈ℕ. (4.16)

First we reproduce the conservative case, 𝛾_𝑖=0,𝑖=1,2,3 and for dissipative case, we consider 𝛾_𝑖=𝜋,𝑖=1,2,3. For the computational domain chosen 32 spatial points and 312 points in the time domain were given by 𝑇=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 Figures 8, 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.

Acknowledgments

The first author acknowledges the support of the CNPq (306338/2008-4) and FAPESPA-PA-Brazil. The authors are thankful to the referees of this paper for the valuable suggestions which improved the paper.