# Conservative Semidiscrete Difference Schemes for Timoshenko Systems

Authors: D. S. Almeida Júnior
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

We present a parameterized family of finite-difference schemes to analyze the energy
properties for linearly elastic constant-coefficient Timoshenko systems considering shear deformation
and rotatory inertia. We derive numerical energies showing the positivity, and the energy conservation
property and we show how to avoid a numerical anomaly known as locking phenomenon on shear force.
Our method of proof relies on discrete multiplier techniques.


---

## Body

## 1. Introduction

The importance of beam theories is well known in the world of engineering and mathematics. The most commonly used beam theories of linear elastic structures are the Euler-Bernoulli beam, the Rayleigh beam, and the Timoshenko beam, with the last being a hyperbolic partial differential equation that has many applications in cutting-edge technologies for flexible structures, such as robotics. The governing equations of Timoshenko beams are considered an improvement over the Euler-Bernoulli and Rayleigh beams since shear deformation is taken into account. The work of Han et al. [1] offers an excellent description of these models for transversely vibrating uniform beams.There is a large number of publications concerning the study of Timoshenko systems [2]. In general, this has been always a description through a system of second-order differential equations, in which the vibration amplitude and the angle due to pure bending were the searched functions. One instructive introduction in [3] contains important references to the early progress in investigating the solutions that represent the mechanical deformations of the Timoshenko model.The 1D equations to the theory of Timoshenko for plane beams are given by(1)
ρ
A
φ
t
t
(
x
,
t
)
=
S
x
(
x
,
t
)
,
ρ
I
ψ
t
t
(
x
,
t
)
=
M
x
(
x
,
t
)
-
S
(
x
,
t
)
,
where t is the time, x is the distance along the center line of beam, φ is the transverse displacement, ψ is the rotation of the neutral axis due to bending, ρ is the mass density of material, M is the bending moment, S is the transverse shear force, A is the cross-section area, and I is the second moment of cross-section area. The relationship between stress and stretch for elastic behavior of the elastic beam is given by
(2)
M
(
x
,
t
)
=
E
I
ψ
x
(
x
,
t
)
,
S
(
x
,
t
)
=
k
A
G
(
φ
x
(
x
,
t
)
+
ψ
(
x
,
t
)
)
,
where E is Young's modulus, G is the modulus of rigidity, and k is the transverse shear factor.Therefore, taking into account (2) and (1), Timoshenko [2] established the following differential equations that incorporate the effects of transverse shear in the cross section of the beam of length L:
(3)
ρ
1
φ
t
t
-
κ
(
φ
x
+
ψ
)
x
=
0
,
in
(
0
,
L
)
×
(
0
,
T
)
,
(4)
ρ
2
ψ
t
t
-
b
ψ
x
x
+
κ
(
φ
x
+
ψ
)
=
0
,
in
(
0
,
L
)
×
(
0
,
T
)
,
where ρ
1
=
ρ
A
,
κ
=
k
G
A
,
ρ
2
=
ρ
I, and b
=
E
I.Now let us mention some numerical aspects concerning 1D Timoshenko systems. Space semidiscretizations are common in many works dealing with Timoshenko beams [4–8]. In general, they arise in the design of numerical schemes to compute approximated solutions of continuous models. For example, a common approach consists in choosing a spatial discretization taking into account the application of finite element methods on partial differential equations. The results are systems of ordinary differential equations in time t.It is well known that standard Galerkin finite element method using equal-order piecewise linear approximations for the rotationsψ and displacement φ yields locking phenomenon. This means that this method produces unsatisfactory numerical results when the thickness is very small. From the point of view of the numerical analysis, the locking phenomenon usually shows itself in that the a priori error estimates for these methods depend on the thickness of the structure in such a way that it degenerates when this parameter becomes small. We recommend the notable works by Arnold in [9], Hughes et al. [10], Prathap and Bhashyam [11], Li [12], Loula et al. [13, 14], and Reddy [15] to a deep discussion of the subject. In particular, the locking phenomenon is characterized by an overestimation on rigidity coefficient given by
(5)
b
h
=
b
(
1
+
κ
12
b
h
2
)
,
in the case of standard finite element method where h is the mesh-size (see [11]).For a finiteh the overestimation b
h is grossly exaggerated by the shear stiffness term and it modifies the bending rigidity. In particular, for plane beams of rectangular geometry with width a and thickness ϵ we get A
=
a
ϵ
,
I
=
a
ϵ
3
/
12 and b
h is rewritten as
(6)
b
h
=
E
I
[
1
+
k
′
G
E
(
h
ϵ
)
2
]
.We can see the locking on shear force in the limit caseϵ
→
0. To avoid it we can take the limit h
→
0, but naturally there exist limitations to computational cost.This paper is mainly concerned with theoretical analysis of finite-difference schemes applied to Timoshenko equations (3)-(4) aiming to identify the locking numbers. To identify them, we use the energy method adapted for our numerical approach. We do not attempt to converge questions for the semidiscrete schemes treated here. We refer the readers to the excellent survey of Zuazua [16] where semidiscrete hyperbolic problems have been considered.The remainder of the paper is organized as follows. In Section2 we established the energy conservation property to continuous model of Timoshenko beams. In Section 3 we introduce numerical schemes in spatial finite difference applied to Timoshenko equations and we show how the locking numbers appear. To do this, we use the θ-scheme. Finally, in Section 4 we conclude our work with future perspectives.


## 2. Energy Conservation Property of the Timoshenko Equations

In this section we describe the energy conservation property with respect to energy functional of the Timoshenko equations.For our purposes we consider here the following boundary conditions of more practical interest that are given byfree-free:(7)
κ
(
φ
x
+
ψ
)
(
0
,
t
)
=
b
ψ
x
(
0
,
t
)
=
0
,
κ
(
φ
x
+
ψ
)
(
L
,
t
)
=
b
ψ
x
(
L
,
t
)
=
0
,
built in-free:(8)
φ
(
0
,
t
)
=
ψ
(
0
,
t
)
=
0
,
κ
(
φ
x
+
ψ
)
(
L
,
t
)
=
b
ψ
x
(
L
,
t
)
=
0
,
hinged-hinged:
(9)
φ
(
0
,
t
)
=
ψ
x
(
0
,
t
)
=
0
,
φ
(
L
,
t
)
=
ψ
x
(
L
,
t
)
=
0
,
built in-hinged:
(10)
φ
(
0
,
t
)
=
ψ
(
0
,
t
)
=
0
,
φ
(
L
,
t
)
=
ψ
x
(
L
,
t
)
=
0
,
built in-built in:
(11)
φ
(
0
,
t
)
=
ψ
(
0
,
t
)
=
φ
(
L
,
t
)
=
ψ
(
L
,
t
)
=
0
,for all t
≥
0. Additionally we consider initial conditions given by
(12)
φ
(
·
,
0
)
=
φ
0
(
·
)
,
φ
t
(
·
,
0
)
=
φ
1
(
·
)
,
ψ
(
·
,
0
)
=
ψ
0
(
·
)
,
ψ
t
(
·
,
0
)
=
ψ
1
(
·
)
,
∀
x
∈
(
0
,
L
)
.There exists an important nonlinear functional related to Timoshenko equations. It is an energy functional composed of potential and kinetic energy. Indeed, taking into consideration any boundary conditions (7)–(11) the energy of solutions for Timoshenko equations (3)-(4) is given by
(13)
E
(
t
)
:
=
ρ
1
2
∫
0
L
‍
|
φ
t
|
2
d
x
+
ρ
2
2
∫
0
L
‍
|
ψ
t
|
2
d
x
+
b
2
∫
0
L
‍
|
ψ
x
|
2
d
x
+
κ
2
∫
0
L
‍
|
φ
x
+
ψ
|
2
d
x
.As usual, we can find that if we multiply formally the first equation in (3) by φ
t we get
(14)
ρ
1
∫
0
L
‍
φ
t
t
φ
t
d
x
-
κ
∫
0
L
‍
(
φ
x
+
ψ
)
x
φ
t
d
x
=
0
,
and then
(15)
ρ
1
2
d
d
t
∫
0
L
‍
φ
t
2
d
x
-
κ
(
φ
x
+
ψ
)
φ
t
|
0
L
+
κ
∫
0
L
‍
(
φ
x
+
ψ
)
φ
t
x
d
x
=
0
.On the other hand, if we multiply formally the second equation in (4) by ψ
t we get
(16)
ρ
2
2
d
d
t
∫
0
L
‍
ψ
t
2
d
x
-
b
ψ
x
ψ
t
|
0
L
+
b
2
d
d
t
∫
0
L
‍
ψ
x
2
d
x
+
κ
∫
0
L
‍
(
φ
x
+
ψ
)
ψ
t
d
x
=
0
.If we add up (15) and (16) we get
(17)
d
d
t
[
ρ
1
2
∫
0
L
‍
φ
t
2
d
x
+
ρ
2
2
∫
0
L
‍
ψ
t
2
d
x
n
n
l
+
b
2
∫
0
L
‍
ψ
x
2
d
x
+
κ
2
∫
0
L
‍
|
φ
x
+
ψ
|
2
d
x
]
-
b
ψ
x
ψ
t
|
0
L
-
κ
(
φ
x
+
ψ
)
φ
t
|
0
L
=
0
,
and then
(18)
d
d
t
E
(
t
)
-
b
ψ
x
ψ
t
|
0
L
-
κ
(
φ
x
+
ψ
)
φ
t
|
0
L
=
0
.Taking any boundary conditions (7)–(11) we obtain that
(19)
E
(
t
)
=
ρ
1
2
∫
0
L
‍
|
φ
1
|
2
d
x
+
ρ
2
2
∫
0
L
‍
|
ψ
1
|
2
d
x
+
b
2
∫
0
L
‍
|
ψ
0
x
|
2
d
x
+
κ
2
∫
0
L
‍
|
φ
0
x
+
ψ
0
|
2
d
x
=
E
(
0
)
,
∀
t
≥
0
.We call this theenergy conservation property and the Timoshenko equations are purely conservative. Therefore, it is interesting and convenient that numerical schemes applied to Timoshenko equations preserve consistently the property (19). That is to say, the numerical energy must be nonlinear and positive and it will obey theenergy conservation property.


## 3. Finite-Difference Semidiscretizations and Properties

In this section, we introduce numerical schemes in finite-difference semi-discretization applied to 1D Timoshenko equations.Firstly, we present the equations used by Wright [8] and we obtain the associated energy (locking free) by using the energy method. This method consists in the use of the discrete multiplier techniques as described by Infante and Zuazua [17].Secondly, we consider theθ-scheme applied to angle rotation ψ and we determine the values of θ to achieve the positivity of the associated energy and to avoid the number b
h in (6). Again, we use the discrete multiplier techniques.For our purposes we considerJ an integer nonnegative and h
=
L
/
(
J
+
1
) a spatial subdivision of the interval (
0
,
L
) given by
(20)
0
=
x
0
<
x
1
<
⋯
<
x
J
<
x
J
+
1
=
L
,
with x
j
=
j
h each node of the mesh. We write the following finite-difference semidiscretization to (3)-(4) as
(21)
ρ
1
φ
j
′′
-
κ
Δ
h
φ
j
-
κ
ψ
j
+
1
-
ψ
j
-
1
2
h
=
0
,
(22)
ρ
2
ψ
j
′′
-
b
Δ
h
ψ
j
+
κ
φ
j
+
1
-
φ
j
-
1
2
h
+
κ
ψ
j
+
1
+
2
ψ
j
+
ψ
j
-
1
4
=
0
,
where we use
(23)
Δ
h
φ
j
:
=
φ
j
+
1
-
2
φ
j
+
φ
j
-
1
h
2
,
Δ
h
ψ
j
:
=
ψ
j
+
1
-
2
ψ
j
+
ψ
j
-
1
h
2
,
to discrete Laplacian in one spatial dimension. Here “′” denotes derivation with respect to time. Moreover, we use φ
j
=
φ
j
(
t
)
,
ψ
j
=
ψ
j
(
t
) for all j
=
1,2
,
…
,
J and t
>
0 to denote the approximate values of the continuous solutions φ
(
j
h
,
t
) and ψ
(
j
h
,
t
) on the mesh, respectively. The homogeneous boundary conditions (7)–(11) in semidiscrete setting are given by(i) free-free:(24)
κ
[
φ
1
-
φ
0
h
+
ψ
1
+
ψ
0
2
]
=
0
,
κ
[
φ
J
+
1
-
φ
J
h
+
ψ
J
+
1
+
ψ
J
2
]
=
0
,
b
ψ
1
-
ψ
0
h
-
κ
h
2
[
φ
1
-
φ
0
h
+
ψ
1
+
ψ
0
2
]
=
0
,
b
ψ
J
+
1
-
ψ
J
h
+
κ
h
2
[
φ
J
+
1
-
φ
J
h
+
ψ
J
+
1
+
ψ
J
2
]
=
0
,(ii) built in-free:(25)
φ
0
=
0
,
κ
[
φ
J
+
1
-
φ
J
h
+
ψ
J
+
1
+
ψ
J
2
]
=
0
,
ψ
0
=
0
,
b
ψ
J
+
1
-
ψ
J
h
+
κ
h
2
[
φ
J
+
1
-
φ
J
h
+
ψ
J
+
1
+
ψ
J
2
]

=

0
,(iii) hinged-hinged:(26)
φ
0
=
φ
J
+
1
=
0
,
b
ψ
1
-
ψ
0
h
-
κ
h
2
[
φ
1
-
φ
0
h
+
ψ
1
+
ψ
0
2
]
=
0
,
b
ψ
J
+
1
-
ψ
J
h
+
κ
h
2
[
φ
J
+
1
-
φ
J
h
+
ψ
J
+
1
+
ψ
J
2
]
=
0
,(iv) built in-hinged:(27)
φ
0
=
φ
J
+
1
=
0
,
ψ
0
=
0
,
b
ψ
J
+
1
-
ψ
J
h
+
κ
h
2
[
φ
J
+
1
-
φ
J
h
+
ψ
J
+
1
+
ψ
J
2
]
=
0
,(v) built in-built in:(28)
φ
0
=
φ
J
+
1
=
ψ
0
=
ψ
J
+
1
=
0
,∀
t
∈
(
0
,
T
). The initial conditions are discretized as follows:
(29)
φ
j
(
0
)
=
φ
j
0
,
φ
j
′
(
0
)
=
φ
j
1
,
ψ
j
(
0
)
=
ψ
j
0
,
ψ
j
′
(
0
)
=
ψ
j
1
,
j
=
0,1
,
2
,
…
,
J
+
1
.Equations (21)-(22) were the object of the study in the works of Wright [7, 8] and the references contained therein. In particular, these equations can be obtained by assembling the stiffness and (lumped) mass matrices of a two-noded beam element which has linear shape functions for rotation ψ and displacement φ. Moreover, these equations are locking free; that is, the associated numerical energy is free of the overestimation (6). See Proposition 1 below.The energy of system (21)–(29) is given by
(30)
E
h
(
t
)
:
=
h
2
∑
j
=
1
J
‍
ρ
1
|
φ
j
′
|
2
+
h
2
∑
j
=
1
J
‍
ρ
2
|
ψ
j
′
|
2
+
h
2
∑
j
=
0
J
‍
[
b
|
ψ
j
+
1
-
ψ
j
h
|
2
n
n
n
l
l
l
l
+
κ
|
φ
j
+
1
-
φ
j
h
+
ψ
j
+
1
+
ψ
j
2
|
2
]
,
which is a discretization of the continuous energy E
(
t
) given in (13). Note that E
h
(
t
) is free of the overestimation (5).Proposition 1 (conservation of energy).
For anyh
>
0 and (
φ
j
(
t
)
,
φ
j
′
(
t
)
,
ψ
j
(
t
)
,
ψ
j
′
(
t
)
) solutions of (21)-(22) and for any boundary conditions (24)–(28) and initial conditions (29), the energy is given by
(31)
E
h
(
t
)
:
=
h
2
∑
j
=
1
J
‍
ρ
1
|
φ
j
′
|
2
+
h
2
∑
j
=
1
J
‍
ρ
2
|
ψ
j
′
|
2
+
h
2
∑
j
=
0
J
‍
[
b
|
ψ
j
+
1
-
ψ
j
h
|
2
n
n
n
l
l
l
l
l
+
κ
|
φ
j
+
1
-
φ
j
h
+
ψ
j
+
1
+
ψ
j
2
|
2
]
,
and it satisfies
(32)
E
h
(
t
)
=
E
h
(
0
)
,
∀
t
∈
[
0
,
T
]
.Proof.
To prove our affirmation we use the energy method in the following manner: we multiply (21)-(22) by using the discrete multipliers h
φ
j
′
(
t
) and h
ψ
j
′
(
t
), respectively, and we take the sum over j
=
1,2
,
…
,
J. Then we obtain
(33)
ρ
1
h
2
d
d
t
∑
j
=
1
J
‍
|
φ
j
′
|
2
=
κ
h
∑
j
=
1
J
‍
[
φ
j
+
1
-
2
φ
j
+
φ
j
-
1
h
2

+

ψ
j
+
1
-
ψ
j
-
1
2
h
]
φ
j
′
,
ρ
2
h
2
d
d
t
∑
j
=
1
J
‍
|
ψ
j
′
|
2
=
b
h
∑
j
=
1
J
[
ψ
j
+
1
-
2
ψ
j
+
ψ
j
-
1
h
2
]
ψ
j
′
‍
-
κ
h
∑
j
=
1
J
‍
φ
j
+
1
-
φ
j
-
1
2
h
ψ
j
′
-
κ
h
∑
j
=
1
J
‍
ψ
j
+
1
+
2
ψ
j
+
ψ
j
-
1
4
ψ
j
′
.
Some elementary calculations show that(34)
κ
h
∑
j
=
1
J
‍
[
φ
j
+
1
-
2
φ
j
+
φ
j
-
1
h
2
+
ψ
j
+
1
-
ψ
j
-
1
2
h
]
φ
j
′
=
κ
h
h
∑
j
=
1
J
‍
[
φ
j
+
1
-
φ
j
h
+
ψ
j
+
1
+
ψ
j
2
]
φ
j
′
-
κ
h
h
∑
j
=
1
J
‍
[
φ
j
-
φ
j
-
1
h
+
ψ
j
+
ψ
j
-
1
2
]
φ
j
′
=
κ
h
h
∑
j
=
1
J
‍
φ
j
+
1
-
φ
j
h
φ
j
′
+
κ
h
h
∑
j
=
1
J
‍
ψ
j
+
1
+
ψ
j
2
φ
j
′
-
κ
h
h
∑
j
=
0
J
-
1
‍
φ
j
+
1
-
φ
j
h
φ
j
+
1
′
-
κ
h
h
∑
j
=
0
J
-
1
‍
ψ
j
+
1
+
ψ
j
2
φ
j
+
1
′
=
κ
h
h
∑
j
=
0
J
‍
φ
j
+
1
-
φ
j
h
φ
j
′
+
κ
h
h
∑
j
=
0
J
‍
ψ
j
+
1
+
ψ
j
2
φ
j
′
-
κ
h
h
∑
j
=
0
J
‍
φ
j
+
1
-
φ
j
h
φ
j
+
1
′
-
κ
h
h
∑
j
=
0
J
‍
ψ
j
+
1
+
ψ
j
2
φ
j
+
1
′
-
κ
φ
1
-
φ
0
h
φ
0
′
-
κ
ψ
1
+
ψ
0
2
φ
0
′
+
κ
φ
J
+
1
-
φ
J
h
φ
J
+
1
′
+
κ
ψ
J
+
1
+
ψ
J
2
φ
J
+
1
′
,
and then we can write(35)
ρ
1
h
2
d
d
t
∑
j
=
1
J
‍
|
φ
j
′
|
2
=
-
κ
h
2
d
d
t
∑
j
=
0
J
‍
|
φ
j
+
1
-
φ
j
h
|
2
︸
:
=
𝒮
0
h
(
t
)
-
κ
h
∑
j
=
0
J
‍
ψ
j
+
1
+
ψ
j
2
φ
j
+
1
′
-
φ
j
′
h
︸
:
=
𝒮
1
h
(
t
)
-
κ
[
φ
1
-
φ
0
h
+
ψ
1
+
ψ
0
2
]
φ
0
′
+
κ
[
φ
J
+
1
-
φ
J
h
+
ψ
J
+
1
+
ψ
J
2
]
φ
J
+
1
′
.
Analogously,(36)
b
h
∑
j
=
1
J
‍
[
ψ
j
+
1
-
2
ψ
j
+
ψ
j
-
1
h
2
]
ψ
j
′
-
κ
h
∑
j
=
1
J
‍
φ
j
+
1
-
φ
j
-
1
2
h
ψ
j
′
-
κ
h
∑
j
=
1
J
‍
ψ
j
+
1
+
2
ψ
j
+
ψ
j
-
1
4
ψ
j
′
=
-
b
h
2
d
d
t
∑
j
=
0
J
‍
|
ψ
j
+
1
-
ψ
j
h
|
2
-
κ
h
2
d
d
t
∑
j
=
0
J
‍
|
ψ
j
+
1
+
ψ
j
2
|
2
-
b
ψ
1
-
ψ
0
h
ψ
0
′
+
b
ψ
J
+
1
-
ψ
J
h
ψ
J
+
1
′
-
κ
h
∑
j
=
0
J
‍
φ
j
+
1
-
φ
j
h
ψ
j
′
+
ψ
j
+
1
′
2
︸
≔
𝒮
2
h
(
t
)

+

κ
h
2
[
φ
1
-
φ
0
h

+

ψ
1
+
ψ
0
2
]
ψ
0
′
+
κ
h
2
[
φ
J
+
1
-
φ
J
h
+
ψ
J
+
1
+
ψ
J
2
]
ψ
J
+
1
′
.
Adding the functionals𝒮
i
h
(
t
)
,
i
=
0,1
,
2, and collecting the terms above we can write
(37)
ρ
1
h
2
d
d
t
∑
j
=
1
J
‍
|
φ
j
′
|
2
+
ρ
2
h
2
d
d
t
∑
j
=
1
J
‍
|
ψ
j
′
|
2
=
-
h
2
κ
d
d
t
∑
j
=
0
J
‍
[
(
φ
j
+
1
-
φ
j
h
)
+
(
ψ
j
+
1
+
ψ
j
2
)
]
2
-
b
h
2
d
d
t
∑
j
=
0
J
‍
|
ψ
j
+
1
-
ψ
j
h
|
2
-
κ
[
φ
1
-
φ
0
h
+
ψ
1
+
ψ
0
2
]
φ
0
′
+
κ
[
φ
J
+
1
-
φ
J
h
+
ψ
J
+
1
+
ψ
J
2
]
φ
J
+
1
′
-
b
ψ
1
-
ψ
0
h
ψ
0
′
+
b
ψ
J
+
1
-
ψ
J
h
ψ
J
+
1
′
+
κ
h
2
[
φ
1
-
φ
0
h
+
ψ
1
+
ψ
0
2
]
ψ
0
′
+
κ
h
2
[
φ
J
+
1
-
φ
J
h
+
ψ
J
+
1
+
ψ
J
2
]
ψ
J
+
1
′
.
We note that the boundary conditions (24)–(28) result naturally from this procedure. Therefore, for any boundary conditions (24)–(28) we get
(38)
d
d
t
E
h
(
t
)
=
0
⟹
E
h
(
t
)
=
E
h
(
0
)
,
∀
t
∈
[
0
,
T
]
,
and therefore the system (21)–(29) is conservative in agreement with continuous counterpart.


### 3.1. The Locking Number fromθ-Scheme

In general, discrete models may not share the same qualitative behavior as their continuous counterparts. This is the case of locking phenomenon on shear force in Timoshenko beams characterized by overestimation (6).Here we show forθ-scheme how locking number appears. We assume a convex linear combination to angle rotation ψ in the following manner: we replace ψ
j
(
t
) by
(39)
θ
ψ
j
+
1
(
t
)
+
(
1
-
2
θ
)
ψ
j
(
t
)
+
θ
ψ
j
-
1
(
t
)
,
∀
t
∈
[
0
,
T
]
.Then we establish a difference scheme for (3)-(4) as follows:
(40)
ρ
1
φ
j
′′
-
κ
Δ
h
φ
j
-
κ
ψ
j
+
1
-
ψ
j
-
1
2
h
=
0
,
ρ
2
ψ
j
′′
-
b
Δ
h
ψ
j
+
κ
φ
j
+
1
-
φ
j
-
1
2
h
+
κ
(
θ
ψ
j
+
1
+
(
1
-
2
θ
)
ψ
j
+
θ
ψ
j
-
1
)
=
0
,
where θ
∈ is a positive parameter. For simplicity we consider homogeneous Dirichlet boundary conditions
(41)
φ
0
(
t
)
=
φ
J
+
1
(
t
)
=
ψ
0
(
t
)
=
ψ
J
+
1
(
t
)
=
0
,
∀
t
≥
0
.The energy of (40)-(41) is given by
(42)
E
θ
h
(
t
)
:
=
h
2
∑
j
=
0
J
‍
[
ρ
1
|
φ
j
′
|
2
+
ρ
2
|
ψ
j
′
|
2
+
b
θ
h
|
ψ
j
+
1
-
ψ
j
h
|
2
n
n
n
l
l
l
+
κ
|
φ
j
+
1
-
φ
j
h
+
ψ
j
+
1
+
ψ
j
2
|
2
]
,
where we have that
(43)
b
θ
h
=
b
[
1
+
κ
b
h
2
(
1
4
-
θ
)
]
.It is clear that the positivity ofE
θ
h
(
t
) holds for θ
∈
[
0,1
/
4
] and only for θ
=
1
/
4 we have that b
(
1
/
4
)
h
=
b. To see how the locking number b
θ
h for any θ
∈
[
0,1
/
4
[ affects the numerical solution, again we consider a rectangular beam, where A
=
a
ϵ
,
I
=
a
ϵ
3
/
12 with width a and thickness ϵ and then
(44)
b
θ
h
=
b
[
1
+
κ
b
h
2
(
1
4
-
θ
)
]
=
E
I
[
1
+
12
k
′
G
E
(
h
ϵ
)
2
(
1
4
-
θ
)
]
.Clearly for0
≤
θ
<
1
/
4 we do not have accurate numerical approximations when the beam thickness is small. This is the case of θ
=
1
/
6 (equivalent to standard finite element method) such that ψ
(
x
j
,
t
)
≈
(
ψ
j
+
1
(
t
)
+
4
ψ
j
(
t
)
+
ψ
j
-
1
(
t
)
)
/
6 and the same problem occurs with standard finite difference given by ψ
(
x
j
,
t
)
≈
ψ
j
(
case
θ
=
0
).Proposition 2.
For anyh
>
0 and (
φ
j
(
t
)
,
φ
j
′
(
t
)
,
ψ
j
(
t
)
,
ψ
j
′
(
t
)
) solutions of (40)-(41) we have
(45)
E
θ
h
(
t
)
=
E
θ
h
(
0
)
,
∀
t
∈
[
0
,
T
]
.Proof.
In what follows we use the results from Proposition1. Then we have
(46)
∑
m
=
0
2
‍
𝒮
m
h
(
t
)
=
-
κ
h
2
d
d
t
∑
j
=
0
J
‍
[
(
φ
j
+
1
-
φ
j
h
)
+
(
ψ
j
+
1
+
ψ
j
2
)
]
2
+
κ
h
2
d
d
t
∑
j
=
0
J
‍
(
ψ
j
+
1
+
ψ
j
2
)
2
.
On the other hand,(47)
κ
h
∑
j
=
1
J
‍
[
θ
ψ
j
+
1
+
(
1
-
2
θ
)
ψ
j
+
θ
ψ
j
-
1
]
ψ
j
′
=
κ
θ
h
∑
j
=
1
J
‍
[
ψ
j
+
1
-
2
ψ
j
+
ψ
j
-
1
]
ψ
j
′
+
κ
h
∑
j
=
1
J
‍
ψ
j
ψ
j
′
,
=
-
κ
θ
h
3
2
d
d
t
∑
j
=
0
J
‍
|
ψ
j
+
1
-
ψ
j
h
|
2
+
κ
h
∑
j
=
1
J
‍
ψ
j
ψ
j
′
.
Therefore,(48)
h
2
d
d
t
∑
j
=
0
J
‍
[
ρ
1
|
φ
j
′
|
2
+
ρ
2
|
ψ
j
′
|
2
]
=
-
h
2
κ
d
d
t
∑
j
=
0
J
‍
[
(
φ
j
+
1
-
φ
j
h
)
+
(
ψ
j
+
1
+
ψ
j
2
)
]
2
-
(
b
-
κ
θ
h
2
)
h
2
d
d
t
∑
j
=
0
J
‍
|
ψ
j
+
1
-
ψ
j
h
|
2
+
κ
h
2
d
d
t
∑
j
=
0
J
‍
|
ψ
j
+
1
+
ψ
j
2
|
2
-
κ
h
2
d
d
t
∑
j
=
0
J
‍
ψ
j
2
︸
:
=
𝒥
h
(
t
)
Now𝒥
h
(
t
) is rewritten as
(49)
𝒥
h
(
t
)
=
κ
h
2
d
d
t
∑
j
=
0
J
‍
|
ψ
j
+
1
+
ψ
j
2
|
2
-
κ
h
2
d
d
t
∑
j
=
0
J
‍
(
ψ
j
+
1
-
ψ
j
)
2
+
κ
h
2
d
d
t
∑
j
=
0
J
‍
[
ψ
j
+
1
2
-
2
ψ
j
+
1
ψ
j
]
=
-
κ
h
3
2
d
d
t
∑
j
=
0
J
‍
(
ψ
j
+
1
-
ψ
j
h
)
2
+
κ
h
2
d
d
t
∑
j
=
0
J
‍
[
ψ
j
+
1
2
+
2
ψ
j
+
1
ψ
j
+
ψ
j
2
4
W
W
W
r
r
r
.
W
W
+
4
ψ
j
+
1
2
-
8
ψ
j
+
1
ψ
j
4
]
=
-
κ
h
3
2
d
d
t
∑
j
=
0
J
‍
(
ψ
j
+
1
-
ψ
j
h
)
2
+
κ
h
2
d
d
t
∑
j
=
0
J
‍
[
5
ψ
j
+
1
2
-
6
ψ
j
+
1
ψ
j
+
ψ
j
2
4
]
.
On the other hand, having in mind the boundary conditions (41), we have that
(50)
∑
j
=
0
J
‍
|
ψ
j
+
1
|
2
=
∑
j
=
0
J
‍
|
ψ
j
|
2
,
and then
(51)
𝒥
h
(
t
)
=
-
κ
h
3
2
d
d
t
∑
j
=
0
J
‍
(
ψ
j
+
1
-
ψ
j
h
)
2
+
κ
h
2
d
d
t
∑
j
=
0
J
‍
[
3
ψ
j
+
1
2
-
6
ψ
j
+
1
ψ
j
+
3
ψ
j
2
4
]
=
-
κ
h
3
2
d
d
t
∑
j
=
0
J
‍
(
ψ
j
+
1
-
ψ
j
h
)
2
+
κ
3
h
3
8
d
d
t
∑
j
=
0
J
‍
(
ψ
j
+
1
-
ψ
j
h
)
2
=
-
κ
h
3
8
d
d
t
∑
j
=
0
J
‍
(
ψ
j
+
1
-
ψ
j
h
)
2
.
Therefore, we obtain that(52)
h
2
d
d
t
∑
j
=
0
J
‍
[
ρ
1
|
φ
j
′
|
2
+
ρ
2
|
ψ
j
′
|
2
]
=
-
h
2
κ
d
d
t
∑
j
=
0
J
‍
[
(
φ
j
+
1
-
φ
j
h
)
+
(
ψ
j
+
1
+
ψ
j
2
)
]
2
-
[
b
+
κ
h
2
(
1
4
-
θ
)
]
h
2
d
d
t
∑
j
=
0
J
‍
|
ψ
j
+
1
-
ψ
j
h
|
2
,
and then we get
(53)
d
d
t
E
θ
h
(
t
)
=
0
⟹
E
θ
h
(
t
)
=
E
θ
h
(
0
)
,
∀
t
∈
[
0
,
T
]
,
where E
θ
h
(
t
) is given by (42).Remark 3.
Deriving the energyE
θ
h
(
t
), collecting the common terms, and using the conservation law (53) we obtain the semidiscrete equations given by
(54)
ρ
1
φ
j
′′
-
κ
Δ
h
φ
j
-
κ
ψ
j
+
1
-
ψ
j
-
1
2
h
=
0
,
ρ
2
ψ
j
′′
-
b
Δ
h
ψ
j
-
κ
h
2
(
1
4
-
θ
)
Δ
h
ψ
j
+
κ
φ
j
+
1
-
φ
j
-
1
2
h
+
κ
ψ
j
+
1
+
2
ψ
j
+
ψ
j
-
1
4
=
0
.
Naturally, we can notice that (54) is inspired in a dispersive approximation of the Timoshenko equations given by
(55)
ρ
1
φ
t
t
-
κ
(
φ
x
+
ψ
)
x
=
0
,
(56)
ρ
2
ψ
t
t
-
b
ψ
x
x
-
κ
h
2
(
1
4
-
θ
)
ψ
x
x
+
κ
(
φ
x
+
ψ
)
=
0
.
For anyθ
∈
[
0,1
/
4
[, (56) is a typical example of a modified PDE. In [18], LeVeque offers an excellent approach of numerical methods applied to PDE that modify the original PDE.


## 3.1. The Locking Number fromθ-Scheme

In general, discrete models may not share the same qualitative behavior as their continuous counterparts. This is the case of locking phenomenon on shear force in Timoshenko beams characterized by overestimation (6).Here we show forθ-scheme how locking number appears. We assume a convex linear combination to angle rotation ψ in the following manner: we replace ψ
j
(
t
) by
(39)
θ
ψ
j
+
1
(
t
)
+
(
1
-
2
θ
)
ψ
j
(
t
)
+
θ
ψ
j
-
1
(
t
)
,
∀
t
∈
[
0
,
T
]
.Then we establish a difference scheme for (3)-(4) as follows:
(40)
ρ
1
φ
j
′′
-
κ
Δ
h
φ
j
-
κ
ψ
j
+
1
-
ψ
j
-
1
2
h
=
0
,
ρ
2
ψ
j
′′
-
b
Δ
h
ψ
j
+
κ
φ
j
+
1
-
φ
j
-
1
2
h
+
κ
(
θ
ψ
j
+
1
+
(
1
-
2
θ
)
ψ
j
+
θ
ψ
j
-
1
)
=
0
,
where θ
∈ is a positive parameter. For simplicity we consider homogeneous Dirichlet boundary conditions
(41)
φ
0
(
t
)
=
φ
J
+
1
(
t
)
=
ψ
0
(
t
)
=
ψ
J
+
1
(
t
)
=
0
,
∀
t
≥
0
.The energy of (40)-(41) is given by
(42)
E
θ
h
(
t
)
:
=
h
2
∑
j
=
0
J
‍
[
ρ
1
|
φ
j
′
|
2
+
ρ
2
|
ψ
j
′
|
2
+
b
θ
h
|
ψ
j
+
1
-
ψ
j
h
|
2
n
n
n
l
l
l
+
κ
|
φ
j
+
1
-
φ
j
h
+
ψ
j
+
1
+
ψ
j
2
|
2
]
,
where we have that
(43)
b
θ
h
=
b
[
1
+
κ
b
h
2
(
1
4
-
θ
)
]
.It is clear that the positivity ofE
θ
h
(
t
) holds for θ
∈
[
0,1
/
4
] and only for θ
=
1
/
4 we have that b
(
1
/
4
)
h
=
b. To see how the locking number b
θ
h for any θ
∈
[
0,1
/
4
[ affects the numerical solution, again we consider a rectangular beam, where A
=
a
ϵ
,
I
=
a
ϵ
3
/
12 with width a and thickness ϵ and then
(44)
b
θ
h
=
b
[
1
+
κ
b
h
2
(
1
4
-
θ
)
]
=
E
I
[
1
+
12
k
′
G
E
(
h
ϵ
)
2
(
1
4
-
θ
)
]
.Clearly for0
≤
θ
<
1
/
4 we do not have accurate numerical approximations when the beam thickness is small. This is the case of θ
=
1
/
6 (equivalent to standard finite element method) such that ψ
(
x
j
,
t
)
≈
(
ψ
j
+
1
(
t
)
+
4
ψ
j
(
t
)
+
ψ
j
-
1
(
t
)
)
/
6 and the same problem occurs with standard finite difference given by ψ
(
x
j
,
t
)
≈
ψ
j
(
case
θ
=
0
).Proposition 2.
For anyh
>
0 and (
φ
j
(
t
)
,
φ
j
′
(
t
)
,
ψ
j
(
t
)
,
ψ
j
′
(
t
)
) solutions of (40)-(41) we have
(45)
E
θ
h
(
t
)
=
E
θ
h
(
0
)
,
∀
t
∈
[
0
,
T
]
.Proof.
In what follows we use the results from Proposition1. Then we have
(46)
∑
m
=
0
2
‍
𝒮
m
h
(
t
)
=
-
κ
h
2
d
d
t
∑
j
=
0
J
‍
[
(
φ
j
+
1
-
φ
j
h
)
+
(
ψ
j
+
1
+
ψ
j
2
)
]
2
+
κ
h
2
d
d
t
∑
j
=
0
J
‍
(
ψ
j
+
1
+
ψ
j
2
)
2
.
On the other hand,(47)
κ
h
∑
j
=
1
J
‍
[
θ
ψ
j
+
1
+
(
1
-
2
θ
)
ψ
j
+
θ
ψ
j
-
1
]
ψ
j
′
=
κ
θ
h
∑
j
=
1
J
‍
[
ψ
j
+
1
-
2
ψ
j
+
ψ
j
-
1
]
ψ
j
′
+
κ
h
∑
j
=
1
J
‍
ψ
j
ψ
j
′
,
=
-
κ
θ
h
3
2
d
d
t
∑
j
=
0
J
‍
|
ψ
j
+
1
-
ψ
j
h
|
2
+
κ
h
∑
j
=
1
J
‍
ψ
j
ψ
j
′
.
Therefore,(48)
h
2
d
d
t
∑
j
=
0
J
‍
[
ρ
1
|
φ
j
′
|
2
+
ρ
2
|
ψ
j
′
|
2
]
=
-
h
2
κ
d
d
t
∑
j
=
0
J
‍
[
(
φ
j
+
1
-
φ
j
h
)
+
(
ψ
j
+
1
+
ψ
j
2
)
]
2
-
(
b
-
κ
θ
h
2
)
h
2
d
d
t
∑
j
=
0
J
‍
|
ψ
j
+
1
-
ψ
j
h
|
2
+
κ
h
2
d
d
t
∑
j
=
0
J
‍
|
ψ
j
+
1
+
ψ
j
2
|
2
-
κ
h
2
d
d
t
∑
j
=
0
J
‍
ψ
j
2
︸
:
=
𝒥
h
(
t
)
Now𝒥
h
(
t
) is rewritten as
(49)
𝒥
h
(
t
)
=
κ
h
2
d
d
t
∑
j
=
0
J
‍
|
ψ
j
+
1
+
ψ
j
2
|
2
-
κ
h
2
d
d
t
∑
j
=
0
J
‍
(
ψ
j
+
1
-
ψ
j
)
2
+
κ
h
2
d
d
t
∑
j
=
0
J
‍
[
ψ
j
+
1
2
-
2
ψ
j
+
1
ψ
j
]
=
-
κ
h
3
2
d
d
t
∑
j
=
0
J
‍
(
ψ
j
+
1
-
ψ
j
h
)
2
+
κ
h
2
d
d
t
∑
j
=
0
J
‍
[
ψ
j
+
1
2
+
2
ψ
j
+
1
ψ
j
+
ψ
j
2
4
W
W
W
r
r
r
.
W
W
+
4
ψ
j
+
1
2
-
8
ψ
j
+
1
ψ
j
4
]
=
-
κ
h
3
2
d
d
t
∑
j
=
0
J
‍
(
ψ
j
+
1
-
ψ
j
h
)
2
+
κ
h
2
d
d
t
∑
j
=
0
J
‍
[
5
ψ
j
+
1
2
-
6
ψ
j
+
1
ψ
j
+
ψ
j
2
4
]
.
On the other hand, having in mind the boundary conditions (41), we have that
(50)
∑
j
=
0
J
‍
|
ψ
j
+
1
|
2
=
∑
j
=
0
J
‍
|
ψ
j
|
2
,
and then
(51)
𝒥
h
(
t
)
=
-
κ
h
3
2
d
d
t
∑
j
=
0
J
‍
(
ψ
j
+
1
-
ψ
j
h
)
2
+
κ
h
2
d
d
t
∑
j
=
0
J
‍
[
3
ψ
j
+
1
2
-
6
ψ
j
+
1
ψ
j
+
3
ψ
j
2
4
]
=
-
κ
h
3
2
d
d
t
∑
j
=
0
J
‍
(
ψ
j
+
1
-
ψ
j
h
)
2
+
κ
3
h
3
8
d
d
t
∑
j
=
0
J
‍
(
ψ
j
+
1
-
ψ
j
h
)
2
=
-
κ
h
3
8
d
d
t
∑
j
=
0
J
‍
(
ψ
j
+
1
-
ψ
j
h
)
2
.
Therefore, we obtain that(52)
h
2
d
d
t
∑
j
=
0
J
‍
[
ρ
1
|
φ
j
′
|
2
+
ρ
2
|
ψ
j
′
|
2
]
=
-
h
2
κ
d
d
t
∑
j
=
0
J
‍
[
(
φ
j
+
1
-
φ
j
h
)
+
(
ψ
j
+
1
+
ψ
j
2
)
]
2
-
[
b
+
κ
h
2
(
1
4
-
θ
)
]
h
2
d
d
t
∑
j
=
0
J
‍
|
ψ
j
+
1
-
ψ
j
h
|
2
,
and then we get
(53)
d
d
t
E
θ
h
(
t
)
=
0
⟹
E
θ
h
(
t
)
=
E
θ
h
(
0
)
,
∀
t
∈
[
0
,
T
]
,
where E
θ
h
(
t
) is given by (42).Remark 3.
Deriving the energyE
θ
h
(
t
), collecting the common terms, and using the conservation law (53) we obtain the semidiscrete equations given by
(54)
ρ
1
φ
j
′′
-
κ
Δ
h
φ
j
-
κ
ψ
j
+
1
-
ψ
j
-
1
2
h
=
0
,
ρ
2
ψ
j
′′
-
b
Δ
h
ψ
j
-
κ
h
2
(
1
4
-
θ
)
Δ
h
ψ
j
+
κ
φ
j
+
1
-
φ
j
-
1
2
h
+
κ
ψ
j
+
1
+
2
ψ
j
+
ψ
j
-
1
4
=
0
.
Naturally, we can notice that (54) is inspired in a dispersive approximation of the Timoshenko equations given by
(55)
ρ
1
φ
t
t
-
κ
(
φ
x
+
ψ
)
x
=
0
,
(56)
ρ
2
ψ
t
t
-
b
ψ
x
x
-
κ
h
2
(
1
4
-
θ
)
ψ
x
x
+
κ
(
φ
x
+
ψ
)
=
0
.
For anyθ
∈
[
0,1
/
4
[, (56) is a typical example of a modified PDE. In [18], LeVeque offers an excellent approach of numerical methods applied to PDE that modify the original PDE.


## 4. Conclusion

In this work, we have developed nonstandard finite-difference semidiscretizations applied to the Timoshenko equations (3)-(4) by using the θ-method. We showed by the discrete energy method that the numerical equations preserve the positivity and the conservation to numerical energy and we identified a numerical anomaly known as locking phenomenon for θ
∈
[
0,1
/
4
[.The results presented here are important from point of view of the theoretical numerical analysis. Actually, an important problem from stabilization theory of the Timoshenko equations (3)-(4) says respect to exponential decay by taking into account some dissipative mechanism. In particular, Muñoz Rivera and Racke [19, 20] showed that (3)-(4) on effect of an only damping of the type ψ
t on angle rotation equations are exponentially stable if and only if κ
/
ρ
1
=
b
/
ρ
2. Therefore, it is interesting to know if the locking free scheme (
θ
=
1
/
4
) preserves the same qualitative behavior.Important issues to be investigated concern observability and control of numerical solutions of Timoshenko equations. To the best of our knowledge, we did not find in the literature numerical problems concerning observability uniform to semidiscrete schemes applied to the Timoshenko equations.



---

*Source: 686421-2014-05-07.xml*