# A Finite Difference Solution of a Simply Supported Beam of Orthotropic Composite Materials Using Displacement Potential Formulation

Authors: S. K. Deb Nath
Journal: Chinese Journal of Engineering
 (2014)

---

## Abstract

Here an efficient displacement potential formulation based finite difference technique is used to solve the elastic field of a simply supported beam of orthotropic composite materials. A simply supported beam made of orthotropic composite material under uniformly distributed loading is considered and its elastic behaviors under such loading conditions are analyzed considering plane stress condition. The solutions of the problem satisfy the force equilibrium conditions as well as boundary conditions. For understanding the elastic behavior of a simply supported beam, the displacement and stress components of some important sections of the beam are shown graphically. Effects of different orthotropic composite materials on the solutions are also analyzed. Besides, at a particular section of the beam, the comparative analysis of the elastic field is carried out by using the FDM and FEM methods.


---

## Body

## 1. Introduction

Strength-to-weight ratio of a fiber reinforced composite material is usually higher than that of the corresponding isotropic material. The use of composite materials gradually increases with time, especially to satisfy the demand of dynamic vehicles such as aircrafts, rockets and so forth, to reduce energy consumptions and to absorb vibrational energy during moving. The use of composite simply supported beams in the construction of engineering structures is quite extensive. It is known that the mechanical properties like strength and toughness of a fiber reinforced composite differ significantly from those of the isotropic materials, which eventually play an important role in defining the state of stress and displacements of the corresponding structure under loading. It is noted that stress analysis of composite structures is usually carried out by the numerical method like the finite element method (FEM) but computations are extremely higher rather than the finite difference based displacement potential numerical approach [1–6]. It is worth mentioning that if the computations of numerical solutions are decreased, the accuracy of the solution will increase and thus the present finite difference solution using the present displacement potential formulation can be used as a standard guideline for checking the reliability and accuracy of approximate solutions.In the solution of structures, the physical conditions of bounding surfaces are mathematically modeled in terms of a mixed mode of boundary conditions, where one of the components of stress as well as displacement on the boundary is assumed to be known. However, the earlier mathematical models of elasticity were inadequate in handling the practical stress problems, as most of them are of mixed-boundary-value type. Since the numerical solution of mixed-boundary-value problems considering low computations, especially with nonisotropic materials, is beyond the scope of the existing mathematical models of elasticity, the use of a new mathematical formulation is investigated to analyze the elastic field of a simply supported beam under uniformly distributed loading.Although the stress analysis has now become a classical subject in the field of solid mechanics, somehow these stress analysis problems are still suffering from many shortcomings and thus are being constantly looked into [7–14]. Elasticity problems are usually formulated in terms of either deformation parameters or stress parameters. Among the existing mathematical models of plane boundary-value stress problems, the stress function approach [15] and the displacement potential formulation [16] are noticeable. The application of the stress function formulation in conjunction with finite difference technique has been reported for the solution of plane elastic problems where all of the boundary conditions are prescribed in terms of stresses only [8, 17]. Further, Conway and Ithaca [18] extended the stress function formulation in the form of Fourier integrals to the case where the material is orthotropic and obtained analytical solutions for a number of ideal problems. The shortcoming of stress approach is that it accepts boundary conditions only in terms of loadings. Boundary restraints specified in terms of displacement components cannot be satisfactorily imposed on the stress function. As most of the practical problems of elasticity are of mixed-boundary type, the stress function approach fails to provide any explicit understanding of the state of stresses at the critical regions of supports and stiffeners. The displacement formulation, on the other hand, involves finding two displacement functions simultaneously from the two second-order elliptic partial differential equations of equilibrium, which is extremely difficult, and this problem becomes more serious when the boundary conditions are mixed [16]. The elastic field of some important structural elements made of isotropic and anisotropic materials is studied analytically and numerically using displacement potential formulation [1–6, 19–38]. It is noted that the recent research and developments in using the displacement potential boundary modeling approach [1–6, 19–38] have generated renewed interest in the field of numerical solutions of stress problems. Recently, Deb Nath [3] proposed the displacement potential formulations for the solution of general anisotropic composite structures and solved cantilever beams with uniformly distributed and point loadings [4], columns with variable compressive loads at its top end [5], a one fixed panel with shear load at its ends [6], and both end fixed beams with distributed loading on its span [19] using finite difference technique. As simply supported composite beams are widely used in different structures, so detailed investigations of the stress field for its perfect design is necessary. Here an alternative numerical solution of a simply supporting beam with uniformly distributed loading using finite difference technique is obtained. Here this paper presents numerical solution of composite simply supported beam subjected to a uniformly distributed loading. Effect of different composites on the solution is also analyzed for simply supported beam. The comparisons of the elastic field obtained by the FDM and FEM of different particular sections of a simply supported beam are also reported.


## 2. Numerical Model of the Problem

Here, a rectangular simply supported beam of orthotropic composite material under uniformly distributed loading is shown in Figure1. The detailed discretization and application of the boundary conditions and governing equation in finite difference form are shown in Figure 2. The aspect ratio of the simply supported beam is 3.0. The orientation of fibers along the length of the beam is θ
=
0
°. The boundary conditions, which are used to solve the problem, are shown in Tables 1
(
a
) and 1
(
b
). The used material properties are shown in Table 2. Here, it is worthy to mention that although the formulations can be applied to any composites, boron/epoxy composite is chosen merely as an example. Furthermore, all the results presented in the study correspond to the value of applied uniformly distributed loading; σ
y
y
0 is 41.4 MPa.Table 1
(a) Boundary conditions of different boundary segments of the simply supported beam as observed in Figure1. (b) Boundary conditions of corner points of the simply supported beam problem as observed in Figure 1.
(a)
Boundary segment
Given boundary conditions
Correspondence between mesh points and given boundary conditions
Mesh point on the physical boundary conditions
Mesh point on the imaginary boundary conditions
A
B
σ
n
,
σ
t
σ
n
=
0
σ
t
=
0
B
C
u
n, σ
n
,
σ
t
At pointE and F, u
n
=
0 Except point E and F, σ
n
=
0
σ
t
=
0
C
D
σ
n
,
σ
t
σ
n
=
0
σ
t
=
0
D
A
σ
n
,
σ
t
σ
n
=
-
σ
y
y
0
σ
t
=
0
(b)
Corner points
Given boundary conditions
Used boundary conditions
Correspondence between mesh points and given boundary conditions
Mesh point on the physical boundary conditions
Mesh point on the imaginary boundary conditions
A
[
(
σ
n
,
σ
t
)
,
(
σ
n
,
σ
t
)
]
σ
n
,
σ
t
σ
n
=
-
σ
y
y
0
σ
t
=
0
B
[
(
σ
n
,
σ
t
)
,
(
σ
n
,
σ
t
)
]
σ
n
,
σ
t
σ
n
=
0
σ
t
=
0
C
[
(
σ
n
,
σ
t
)
,
(
σ
n
,
σ
t
)
]
σ
n
,
σ
t
σ
n
=
0
σ
t
=
0
D
[
(
σ
n
,
σ
t
)
,
(
σ
n
,
σ
t
)
]
σ
n
,
σ
t
σ
n
=
-
σ
y
y
0
σ
t
=
0Table 2
Properties of composites used to obtain numerical results.
Material
Property
Boron/epoxy
Fiber
E
f (103 MPa)
414
ν
f
0.20
Resin
E
r (103 MPa)
3.45
ν
r
0.35
Composite
E
1 (103 MPa)
282.9
E
2 (103 MPa)
24.2
G
12 (103 MPa)
10.4
ν
12
0.27
ν
21
0.023Figure 1
A simply supported beam of orthotropic composite material under uniformly distributed loading.Figure 2
Application of the normal and tangential displacement, stress components, and the governing equation as finite difference forms.


## 3. Displacement Potential Function Formulation

With reference to a rectangular coordinate system(
x
,
y
), the differential equations of equilibrium for the plane stress problems of orthotropic composite materials having the fiber orientation, θ
=
0
0, are as follows [3, 4, 39]:
(1)
E
1
2
E
1
-
ν
12
2
E
2
∂
2
u
x
∂
x
2
+
(
ν
12
E
1
E
2
E
1
-
ν
12
2
E
2
+
G
12
)
×
∂
2
u
y
∂
x
∂
y
+
G
12
∂
2
u
x
∂
y
2
=
0
,
(2)
E
1
E
2
E
1
-
ν
12
2
E
2
∂
2
u
y
∂
y
2
+
(
ν
12
E
1
E
2
E
1
-
ν
12
2
E
2
+
G
12
)
×
∂
2
u
x
∂
x
∂
y
+
G
12
∂
2
u
y
∂
x
2
=
0
.
The displacement components are expressed in terms of a displacement potential function, ψ
(
x
,
y
), as follows:
(3)
u
x
(
x
,
y
)
=
α
1
∂
2
ψ
∂
x
2
+
α
2
∂
2
ψ
∂
x
∂
y
+
α
3
∂
2
ψ
∂
y
2
,
u
y
(
x
,
y
)
=
α
4
∂
2
ψ
∂
x
2
+
α
5
∂
2
ψ
∂
x
∂
y
+
α
6
∂
2
ψ
∂
y
2
,
where α
i are unknown material constants.The values ofα are obtained by substituting the displacement components given by (3) into the equilibrium equation (1). The values of the constants thus obtained for the orthotropic materials with θ
=
0
0 are as follows [3, 4]:
(4)
α
1
=
α
3
=
α
5
=
0
,
α
2
=
1
,
α
4
=
-
E
1
2
Z
11
,
α
6
=
-
G
12
(
E
1
-
ν
12
2
E
2
)
Z
11
,
where Z
11
=
ν
12
E
1
E
2
+
G
12
(
E
1
-
ν
12
2
E
2
).When the values of the constantα are substituted in the equilibrium equation (2), the equation turns into a fourth-order partial differential equation that governs the displacement potential function, ψ
(
x
,
y
). The governing equation in terms of the function ψ is as follows:
(5)
E
1
G
12
∂
4
ψ
∂
x
4
+
E
2
(
E
1
-
2
ν
12
G
12
)
∂
4
ψ
∂
x
2
∂
y
2
+
E
2
G
12
∂
4
ψ
∂
y
4
=
0
.
Similarly, substituting the value of α in (3) gives the two displacement components in terms of the displacement potential function as follows:
(6)
u
x
(
x
,
y
)
=
∂
2
ψ
∂
x
∂
y
,
u
y
(
x
,
y
)
=
-
1
Z
11
[
E
1
2
∂
2
ψ
∂
x
2
+
G
12
(
E
1
-
ν
12
2
E
2
)
∂
2
ψ
∂
y
2
]
.
Using (6), the stress-strain relations, and reciprocal relationships [2] having fiber orientation θ
=
0
0, the stress components in terms of the potential function ψ are obtained as follows [3, 4]:
(7)
σ
x
x
(
x
,
y
)
=
E
1
G
12
Z
11
[
E
1
∂
3
ψ
∂
x
2
∂
y
-
ν
12
E
2
∂
3
ψ
∂
y
3
]
,
σ
y
y
(
x
,
y
)
=
E
1
E
2
Z
11
[
(
ν
12
G
12
-
E
1
)
∂
3
ψ
∂
x
2
∂
y
-
G
12
∂
3
ψ
∂
y
3
]
,
σ
x
y
(
x
,
y
)
=
-
E
1
G
12
Z
11
[
E
1
∂
3
ψ
∂
x
3
-
ν
12
E
2
∂
3
ψ
∂
x
∂
y
2
]
.
From these expressions it is found that, as far as boundary conditions are concerned, either known restraints or known stresses or any combination of these can readily be converted to finite difference expressions in terms of ψ on the boundary. To analyze the state of displacements and stresses for two-dimensional regular/irregular-shaped bodies with mixed boundary condition, the fourth-order-partial differential equation (5) together with the displacement and stress components given by (6) to (7) is considered. The expressions for the displacement and stress components as given by (6) to (7) are valid for points within the body as well as on the boundary. But, for the case of practical problems, the boundary conditions are known in terms of their normal and tangential components on the boundary. Thus, the boundary equations for the known conditions on the boundary having arbitrary shape are
(8)
u
n
(
x
,
y
)
=
u
x
·
l
+
u
y
·
m
=
-
E
1
2
m
Z
11
∂
2
ψ
∂
x
2
-
G
12
(
E
1
-
ν
12
2
E
2
)
m
Z
11
∂
2
ψ
∂
y
2
+
l
∂
2
ψ
∂
x
∂
y
,
u
t
(
x
,
y
)
=
u
y
·
l
-
u
x
·
m
=
-
E
1
2
l
Z
11
∂
2
ψ
∂
x
2
-
l
G
12
(
E
1
-
ν
12
2
E
2
)
Z
11
∂
2
ψ
∂
y
2
-
m
∂
2
ψ
∂
x
∂
y
,
σ
n
(
x
,
y
)
=
l
2
σ
x
x
(
x
,
y
)
+
2
l
m
σ
x
y
(
x
,
y
)
+
m
2
σ
y
y
(
x
,
y
)
=
E
1
Z
11
[
-
2
l
m
E
1
G
12
∂
3
ψ
∂
x
3
+
{
E
1
G
12
l
2
+
(
ν
12
G
12
-
E
1
)
m
2
E
2
}
∂
3
ψ
∂
x
2
∂
y
+
2
l
m
E
2
G
2
ν
12
∂
3
ψ
∂
x
∂
y
2
-
E
2
(
G
12
ν
12
l
2
+
G
12
m
2
)
∂
3
ψ
∂
y
3
]
,
σ
t
(
x
,
y
)
=
(
l
2
-
m
2
)
σ
x
y
(
x
,
y
)
+
l
m
{
σ
y
y
(
x
,
y
)
-
σ
x
x
(
x
,
y
)
}
=
E
1
Z
11
[
-
(
l
2
-
m
2
)
E
1
G
12
∂
3
ψ
∂
x
3
+
{
E
2
(
ν
12
G
12
-
E
1
)
-
E
1
G
12
}
l
m
∂
3
ψ
∂
x
2
∂
y
+
(
l
2
-
m
2
)
E
2
G
12
ν
12
∂
3
ψ
∂
x
∂
y
2
+
l
m
E
2
G
12
(
ν
12
-
1
)
∂
3
ψ
∂
y
3
]
.


## 4. Finite Difference Modeling of the Problem

### 4.1. Finite Difference Discretization of the Governing Equation

The present computational scheme involves evaluation of a single function,ψ, at the mesh points of a uniform two-dimensional rectangular mesh network used to discretize the domain. The governing equation (5), used to evaluate the function ψ only at the internal mesh points, can be expressed in its corresponding difference form when all the derivatives are replaced by corresponding central difference expressions. The finite difference form of (5) is as follows:
(9)
z
k
1
{
ψ
(
i
-
2
,
j
)
+
ψ
(
i
+
2
,
j
)
}
-
z
k
2
{
ψ
(
i
-
1
,
j
)
+
ψ
(
i
+
1
,
j
)
}
-
z
k
3
{
ψ
(
i
,
j
+
1
)
+
ψ
(
i
,
j
-
1
)
}
+
z
k
4
ψ
(
i
,
j
)
+
z
k
5
{
ψ
(
i
+
1
,
j
+
1
)
+
ψ
(
i
+
1
,
j
-
1
)
+
ψ
(
i
-
1
,
j
+
1
)
+
ψ
(
i
-
1
,
j
-
1
)
}
+
Q
{
ψ
(
i
,
j
-
2
)
+
ψ
(
i
,
j
+
2
)
}
=
0
,
where z
k
1
=
R
4, z
k
2
=
4
R
4
+
2
P
R
2, z
k
3
=
2
P
R
2
+
4
Q, z
k
4
=
6
R
4
+
4
P
R
2
+
6
Q
(10)
z
k
5
=
P
R
2
,
P
=
(
E
2
G
12
-
2
μ
12
E
2
E
1
)
,
Q
=
E
2
E
1
,
R
=
k
h
.


### 4.2. Application of the Governing Equation

In finite difference method, the region occupied by the structure under consideration is divided into fine meshes to get nodal points as shown in Figure2. The solid line represents the physical boundary conditions of the structure. The discretized form of governing differential equation is applied at each of the nodal points (
i
,
j
) inside the physical boundary. This gives a set of algebraic equations for the determination of an unknown function, ψ. However, the number of unknown is greater than the number of equations, which makes the problem intractable. This is because the application of the governing equation to the interior points also involves the points on and outside the physical boundary. The line formed by connecting the points outside the physical boundary is called imaginary boundary as shown in Figure 2. The number of unknowns is equal to the number of nodal points including boundary and exterior points. To make the problem tractable, it is necessary to generate more equations. This requirement is met up by applying the discretized boundary conditions at the nodal points on the boundary and this gives a complete set of simultaneous equations which is solved by LU decomposition.


### 4.3. Boundary Conditions in Finite Difference Form

The boundary conditions are discretized in such a way that their applications does not involve additional points outside the imaginary boundary. The whole body of the beam is divided into four segments. To satisfy this requirement, the same boundary conditions must be required to discretize in different forms such asi-forward and j-forward; i-backward and j-forward; i-backward and j-backward; i-forward and j-backward for bottom left, bottom right, top right, and top left segments, respectively.


#### 4.3.1. Displacement and Stress Components at the Bottom Left Segment of the Boundary in Finite Difference Forms

Consider(11)
u
x
(
x
,
y
)
=
c
1
[
9
ψ
(
i
,
j
)
-
12
{
ψ
(
i
,
j
+
1
)
+
ψ
(
i
+
1
,
j
)
}
+
16
ψ
(
i
+
1
,
j
+
1
)
+
3
{
ψ
(
i
,
j
+
2
)
+
ψ
(
i
+
2
,
j
)
}
-
4
{
ψ
(
i
+
1
,
j
+
2
)
+
ψ
(
i
+
2
,
j
+
1
)
}
+
ψ
(
i
+
2
,
j
+
2
)
]
,
u
y
(
x
,
y
)
=
c
2
{
ψ
(
i
-
1
,
j
)
+
ψ
(
i
+
1
,
j
)
}
+
c
3
{
ψ
(
i
,
j
+
1
)
+
ψ
(
i
,
j
-
1
)
}
-
2
(
c
2
+
c
3
)
ψ
(
i
,
j
)
,
σ
x
x
(
x
,
y
)
=
-
3
c
4
{
ψ
(
i
-
1
,
j
)
+
ψ
(
i
+
1
,
j
)
}
+
4
c
4
{
ψ
(
i
-
1
,
j
+
1
)
+
ψ
(
i
+
1
,
j
+
1
)
}
+
(
6
c
4

+

10
c
5
)
ψ
(
i
,
j
)

-

(
8
c
4

+

12
c
5
)
ψ
(
i
,
j

+
1
)
+
(
2
c
4
+
6
c
5
)
ψ
(
i
,
j
+
2
)
-
c
4
ψ
(
i
-
1
,
j
+
2
)
-
c
4
ψ
(
i
+
1
,
j
+
2
)
-
3
c
5
ψ
(
i
,
j
-
1
)
-
c
5
ψ
(
i
,
j
+
3
)
,
σ
y
y
(
x
,
y
)
=
-
3
c
6
{
ψ
(
i
-
1
,
j
)
+
3
c
6
ψ
(
i
+
1
,
j
)
}
+
4
c
6
{
ψ
(
i
-
1
,
j
+
1
)
+
ψ
(
i
+
1
,
j
+
1
)
}
+
(
6
c
6

+

10
c
7
)
ψ
(
i
,
j
)

-

(
8
c
6

+

12
c
7
)
ψ
(
i
,
j

+
1
)
+
(
2
c
6
+
6
c
7
)
ψ
(
i
,
j
+
2
)
-
c
6
ψ
(
i
-
1
,
j
+
2
)
-
c
6
ψ
(
i
+
1
,
j
+
2
)
-
3
c
7
ψ
(
i
,
j
-
1
)
-
c
7
ψ
(
i
,
j
+
3
)
,
σ
x
y
(
x
,
y
)
=
-
3
c
9
ψ
(
i
-
1
,
j
)
+
(
10
c
9
+
6
c
8
)
ψ
(
i
,
j
)
-
(
12
c
9
+
8
c
8
)
ψ
(
i
+
1
,
j
)
+
(
2
c
8
+
6
c
9
)
×
ψ
(
i
+
2
,
j
)
-
c
9
ψ
(
i
+
3
,
j
)
-
3
c
8
{
ψ
(
i
,
j
-
1
)
+
ψ
(
i
,
j
+
1
)
}
+
4
c
8
{
ψ
(
i
+
1
,
j
-
1
)
+
ψ
(
i
+
1
,
j
+
1
)
}
-
c
8
{
ψ
(
i
+
2
,
j
-
1
)
+
ψ
(
i
+
2
,
j
+
1
)
}
,
where
(12)
c
1
=
1
4
k
h
,
c
2
=
-
E
1
2
Z
11
h
2
,
c
3
=
-
G
12
(
E
1
-
ν
12
2
E
2
)
Z
11
k
2
,
c
4
=
E
1
2
G
12
2
Z
11
h
2
k
,
c
5
=
-
E
1
E
2
ν
12
G
12
2
Z
11
k
3
,
c
6
=
-
E
1
E
2
(
E
1
-
ν
12
G
12
)
2
Z
11
h
2
k
,
c
7
=
-
E
1
E
2
G
12
2
Z
11
k
3
,
c
8
=
E
1
E
2
ν
12
G
12
2
Z
11
h
k
2
,
c
9
=
-
E
1
2
G
12
2
Z
11
h
3
.
The discretized forms of boundary conditions derived above cannot be conveniently used for boundaries which are not parallel to the x- and y-axis. In general, the components of stress and displacements which are normal and tangential to the boundaries are used as boundary conditions. The normal and tangential form of displacement and stress components are below:
(13)
u
n
(
i
,
j
)
=
u
x
(
i
,
j
)
·
l
+
u
y
(
i
,
j
)
·
m
,
u
t
(
i
,
j
)
=
u
y
(
i
,
j
)
·
l
-
u
x
(
i
,
j
)
·
m
,
σ
n
(
i
,
j
)
=
l
2
σ
x
x
(
i
,
j
)
+
2
l
m
σ
x
y
(
i
,
j
)
+
m
2
σ
y
y
(
i
,
j
)
,
σ
t
(
i
,
j
)
=
(
l
2
-
m
2
)
σ
x
y
(
i
,
j
)
+
l
m
{
σ
y
y
(
i
,
j
)
-
σ
x
x
(
i
,
j
)
}
.


## 4.1. Finite Difference Discretization of the Governing Equation

The present computational scheme involves evaluation of a single function,ψ, at the mesh points of a uniform two-dimensional rectangular mesh network used to discretize the domain. The governing equation (5), used to evaluate the function ψ only at the internal mesh points, can be expressed in its corresponding difference form when all the derivatives are replaced by corresponding central difference expressions. The finite difference form of (5) is as follows:
(9)
z
k
1
{
ψ
(
i
-
2
,
j
)
+
ψ
(
i
+
2
,
j
)
}
-
z
k
2
{
ψ
(
i
-
1
,
j
)
+
ψ
(
i
+
1
,
j
)
}
-
z
k
3
{
ψ
(
i
,
j
+
1
)
+
ψ
(
i
,
j
-
1
)
}
+
z
k
4
ψ
(
i
,
j
)
+
z
k
5
{
ψ
(
i
+
1
,
j
+
1
)
+
ψ
(
i
+
1
,
j
-
1
)
+
ψ
(
i
-
1
,
j
+
1
)
+
ψ
(
i
-
1
,
j
-
1
)
}
+
Q
{
ψ
(
i
,
j
-
2
)
+
ψ
(
i
,
j
+
2
)
}
=
0
,
where z
k
1
=
R
4, z
k
2
=
4
R
4
+
2
P
R
2, z
k
3
=
2
P
R
2
+
4
Q, z
k
4
=
6
R
4
+
4
P
R
2
+
6
Q
(10)
z
k
5
=
P
R
2
,
P
=
(
E
2
G
12
-
2
μ
12
E
2
E
1
)
,
Q
=
E
2
E
1
,
R
=
k
h
.


## 4.2. Application of the Governing Equation

In finite difference method, the region occupied by the structure under consideration is divided into fine meshes to get nodal points as shown in Figure2. The solid line represents the physical boundary conditions of the structure. The discretized form of governing differential equation is applied at each of the nodal points (
i
,
j
) inside the physical boundary. This gives a set of algebraic equations for the determination of an unknown function, ψ. However, the number of unknown is greater than the number of equations, which makes the problem intractable. This is because the application of the governing equation to the interior points also involves the points on and outside the physical boundary. The line formed by connecting the points outside the physical boundary is called imaginary boundary as shown in Figure 2. The number of unknowns is equal to the number of nodal points including boundary and exterior points. To make the problem tractable, it is necessary to generate more equations. This requirement is met up by applying the discretized boundary conditions at the nodal points on the boundary and this gives a complete set of simultaneous equations which is solved by LU decomposition.


## 4.3. Boundary Conditions in Finite Difference Form

The boundary conditions are discretized in such a way that their applications does not involve additional points outside the imaginary boundary. The whole body of the beam is divided into four segments. To satisfy this requirement, the same boundary conditions must be required to discretize in different forms such asi-forward and j-forward; i-backward and j-forward; i-backward and j-backward; i-forward and j-backward for bottom left, bottom right, top right, and top left segments, respectively.


### 4.3.1. Displacement and Stress Components at the Bottom Left Segment of the Boundary in Finite Difference Forms

Consider(11)
u
x
(
x
,
y
)
=
c
1
[
9
ψ
(
i
,
j
)
-
12
{
ψ
(
i
,
j
+
1
)
+
ψ
(
i
+
1
,
j
)
}
+
16
ψ
(
i
+
1
,
j
+
1
)
+
3
{
ψ
(
i
,
j
+
2
)
+
ψ
(
i
+
2
,
j
)
}
-
4
{
ψ
(
i
+
1
,
j
+
2
)
+
ψ
(
i
+
2
,
j
+
1
)
}
+
ψ
(
i
+
2
,
j
+
2
)
]
,
u
y
(
x
,
y
)
=
c
2
{
ψ
(
i
-
1
,
j
)
+
ψ
(
i
+
1
,
j
)
}
+
c
3
{
ψ
(
i
,
j
+
1
)
+
ψ
(
i
,
j
-
1
)
}
-
2
(
c
2
+
c
3
)
ψ
(
i
,
j
)
,
σ
x
x
(
x
,
y
)
=
-
3
c
4
{
ψ
(
i
-
1
,
j
)
+
ψ
(
i
+
1
,
j
)
}
+
4
c
4
{
ψ
(
i
-
1
,
j
+
1
)
+
ψ
(
i
+
1
,
j
+
1
)
}
+
(
6
c
4

+

10
c
5
)
ψ
(
i
,
j
)

-

(
8
c
4

+

12
c
5
)
ψ
(
i
,
j

+
1
)
+
(
2
c
4
+
6
c
5
)
ψ
(
i
,
j
+
2
)
-
c
4
ψ
(
i
-
1
,
j
+
2
)
-
c
4
ψ
(
i
+
1
,
j
+
2
)
-
3
c
5
ψ
(
i
,
j
-
1
)
-
c
5
ψ
(
i
,
j
+
3
)
,
σ
y
y
(
x
,
y
)
=
-
3
c
6
{
ψ
(
i
-
1
,
j
)
+
3
c
6
ψ
(
i
+
1
,
j
)
}
+
4
c
6
{
ψ
(
i
-
1
,
j
+
1
)
+
ψ
(
i
+
1
,
j
+
1
)
}
+
(
6
c
6

+

10
c
7
)
ψ
(
i
,
j
)

-

(
8
c
6

+

12
c
7
)
ψ
(
i
,
j

+
1
)
+
(
2
c
6
+
6
c
7
)
ψ
(
i
,
j
+
2
)
-
c
6
ψ
(
i
-
1
,
j
+
2
)
-
c
6
ψ
(
i
+
1
,
j
+
2
)
-
3
c
7
ψ
(
i
,
j
-
1
)
-
c
7
ψ
(
i
,
j
+
3
)
,
σ
x
y
(
x
,
y
)
=
-
3
c
9
ψ
(
i
-
1
,
j
)
+
(
10
c
9
+
6
c
8
)
ψ
(
i
,
j
)
-
(
12
c
9
+
8
c
8
)
ψ
(
i
+
1
,
j
)
+
(
2
c
8
+
6
c
9
)
×
ψ
(
i
+
2
,
j
)
-
c
9
ψ
(
i
+
3
,
j
)
-
3
c
8
{
ψ
(
i
,
j
-
1
)
+
ψ
(
i
,
j
+
1
)
}
+
4
c
8
{
ψ
(
i
+
1
,
j
-
1
)
+
ψ
(
i
+
1
,
j
+
1
)
}
-
c
8
{
ψ
(
i
+
2
,
j
-
1
)
+
ψ
(
i
+
2
,
j
+
1
)
}
,
where
(12)
c
1
=
1
4
k
h
,
c
2
=
-
E
1
2
Z
11
h
2
,
c
3
=
-
G
12
(
E
1
-
ν
12
2
E
2
)
Z
11
k
2
,
c
4
=
E
1
2
G
12
2
Z
11
h
2
k
,
c
5
=
-
E
1
E
2
ν
12
G
12
2
Z
11
k
3
,
c
6
=
-
E
1
E
2
(
E
1
-
ν
12
G
12
)
2
Z
11
h
2
k
,
c
7
=
-
E
1
E
2
G
12
2
Z
11
k
3
,
c
8
=
E
1
E
2
ν
12
G
12
2
Z
11
h
k
2
,
c
9
=
-
E
1
2
G
12
2
Z
11
h
3
.
The discretized forms of boundary conditions derived above cannot be conveniently used for boundaries which are not parallel to the x- and y-axis. In general, the components of stress and displacements which are normal and tangential to the boundaries are used as boundary conditions. The normal and tangential form of displacement and stress components are below:
(13)
u
n
(
i
,
j
)
=
u
x
(
i
,
j
)
·
l
+
u
y
(
i
,
j
)
·
m
,
u
t
(
i
,
j
)
=
u
y
(
i
,
j
)
·
l
-
u
x
(
i
,
j
)
·
m
,
σ
n
(
i
,
j
)
=
l
2
σ
x
x
(
i
,
j
)
+
2
l
m
σ
x
y
(
i
,
j
)
+
m
2
σ
y
y
(
i
,
j
)
,
σ
t
(
i
,
j
)
=
(
l
2
-
m
2
)
σ
x
y
(
i
,
j
)
+
l
m
{
σ
y
y
(
i
,
j
)
-
σ
x
x
(
i
,
j
)
}
.


## 4.3.1. Displacement and Stress Components at the Bottom Left Segment of the Boundary in Finite Difference Forms

Consider(11)
u
x
(
x
,
y
)
=
c
1
[
9
ψ
(
i
,
j
)
-
12
{
ψ
(
i
,
j
+
1
)
+
ψ
(
i
+
1
,
j
)
}
+
16
ψ
(
i
+
1
,
j
+
1
)
+
3
{
ψ
(
i
,
j
+
2
)
+
ψ
(
i
+
2
,
j
)
}
-
4
{
ψ
(
i
+
1
,
j
+
2
)
+
ψ
(
i
+
2
,
j
+
1
)
}
+
ψ
(
i
+
2
,
j
+
2
)
]
,
u
y
(
x
,
y
)
=
c
2
{
ψ
(
i
-
1
,
j
)
+
ψ
(
i
+
1
,
j
)
}
+
c
3
{
ψ
(
i
,
j
+
1
)
+
ψ
(
i
,
j
-
1
)
}
-
2
(
c
2
+
c
3
)
ψ
(
i
,
j
)
,
σ
x
x
(
x
,
y
)
=
-
3
c
4
{
ψ
(
i
-
1
,
j
)
+
ψ
(
i
+
1
,
j
)
}
+
4
c
4
{
ψ
(
i
-
1
,
j
+
1
)
+
ψ
(
i
+
1
,
j
+
1
)
}
+
(
6
c
4

+

10
c
5
)
ψ
(
i
,
j
)

-

(
8
c
4

+

12
c
5
)
ψ
(
i
,
j

+
1
)
+
(
2
c
4
+
6
c
5
)
ψ
(
i
,
j
+
2
)
-
c
4
ψ
(
i
-
1
,
j
+
2
)
-
c
4
ψ
(
i
+
1
,
j
+
2
)
-
3
c
5
ψ
(
i
,
j
-
1
)
-
c
5
ψ
(
i
,
j
+
3
)
,
σ
y
y
(
x
,
y
)
=
-
3
c
6
{
ψ
(
i
-
1
,
j
)
+
3
c
6
ψ
(
i
+
1
,
j
)
}
+
4
c
6
{
ψ
(
i
-
1
,
j
+
1
)
+
ψ
(
i
+
1
,
j
+
1
)
}
+
(
6
c
6

+

10
c
7
)
ψ
(
i
,
j
)

-

(
8
c
6

+

12
c
7
)
ψ
(
i
,
j

+
1
)
+
(
2
c
6
+
6
c
7
)
ψ
(
i
,
j
+
2
)
-
c
6
ψ
(
i
-
1
,
j
+
2
)
-
c
6
ψ
(
i
+
1
,
j
+
2
)
-
3
c
7
ψ
(
i
,
j
-
1
)
-
c
7
ψ
(
i
,
j
+
3
)
,
σ
x
y
(
x
,
y
)
=
-
3
c
9
ψ
(
i
-
1
,
j
)
+
(
10
c
9
+
6
c
8
)
ψ
(
i
,
j
)
-
(
12
c
9
+
8
c
8
)
ψ
(
i
+
1
,
j
)
+
(
2
c
8
+
6
c
9
)
×
ψ
(
i
+
2
,
j
)
-
c
9
ψ
(
i
+
3
,
j
)
-
3
c
8
{
ψ
(
i
,
j
-
1
)
+
ψ
(
i
,
j
+
1
)
}
+
4
c
8
{
ψ
(
i
+
1
,
j
-
1
)
+
ψ
(
i
+
1
,
j
+
1
)
}
-
c
8
{
ψ
(
i
+
2
,
j
-
1
)
+
ψ
(
i
+
2
,
j
+
1
)
}
,
where
(12)
c
1
=
1
4
k
h
,
c
2
=
-
E
1
2
Z
11
h
2
,
c
3
=
-
G
12
(
E
1
-
ν
12
2
E
2
)
Z
11
k
2
,
c
4
=
E
1
2
G
12
2
Z
11
h
2
k
,
c
5
=
-
E
1
E
2
ν
12
G
12
2
Z
11
k
3
,
c
6
=
-
E
1
E
2
(
E
1
-
ν
12
G
12
)
2
Z
11
h
2
k
,
c
7
=
-
E
1
E
2
G
12
2
Z
11
k
3
,
c
8
=
E
1
E
2
ν
12
G
12
2
Z
11
h
k
2
,
c
9
=
-
E
1
2
G
12
2
Z
11
h
3
.
The discretized forms of boundary conditions derived above cannot be conveniently used for boundaries which are not parallel to the x- and y-axis. In general, the components of stress and displacements which are normal and tangential to the boundaries are used as boundary conditions. The normal and tangential form of displacement and stress components are below:
(13)
u
n
(
i
,
j
)
=
u
x
(
i
,
j
)
·
l
+
u
y
(
i
,
j
)
·
m
,
u
t
(
i
,
j
)
=
u
y
(
i
,
j
)
·
l
-
u
x
(
i
,
j
)
·
m
,
σ
n
(
i
,
j
)
=
l
2
σ
x
x
(
i
,
j
)
+
2
l
m
σ
x
y
(
i
,
j
)
+
m
2
σ
y
y
(
i
,
j
)
,
σ
t
(
i
,
j
)
=
(
l
2
-
m
2
)
σ
x
y
(
i
,
j
)
+
l
m
{
σ
y
y
(
i
,
j
)
-
σ
x
x
(
i
,
j
)
}
.


## 5. Results and Discussion

### 5.1. Analysis of Elastic Behavior of the Simply Supported Beam

The geometry of the problem is shown in Figure1. Figure 3 presents the normalized displacement components u
x as a function of x. At x
/
b
=
0.5, the axial displacement of each section is zero. The nature of the axial displacement above and below the neutral axis is opposite and antisymmetric. Towards the center (x
/
b
=
0.5) of each section, the magnitude of u
x gradually decreases to zero. On the other hand, the axial displacement component of each section with respect to x is the highest at the left and right ends of the beam, which gradually decreases to zero towards the center x
/
b
=
0.5. Figure 4 describes the normalized lateral displacement component u
y with respect to x
/
b at different sections of the beam. Middle region (x
/
b
=
0.5) of each section shows the highest lateral displacement. Figure 5 shows the deformed shape of the beam which shows that deformation in y direction is more significant.Figure 3
Distribution of the normalized displacement component,
u
x, as a function of x
/
b at different sections of the beam.Figure 4
Distribution of the normalized displacement component,u
y, as a function of x
/
b at different sections of the beam.Figure 5
Deformed shape of the cantilever beam under uniformly distributed loading.Figure6 describes the normalized axial stress component, σ
x
x, as a function of x
/
b. The stress σ
x
x is zero at the neutral axis y
/
a
=
0.5, above which σ
x
x is negative and below which σ
x
x is positive. The negative and positive stresses correspond that the regions are in compression and tension, respectively. This agrees with the physical characteristics of the problem. Figure 7 shows that σ
y
y
/
σ
y
y
0 is unity at y
/
a
=
1.0 for any value of x
/
b. This tells that the computed stress is equal to the applied stress at the top surface which verifies the boundary condition. Further, σ
y
y
/
σ
y
y
0 is maximum at the two supports of the bottom surface (y
/
a
=
0). This is due to the fact that the two supports equalize the total load applied at the top surface. Figure 8 exhibits the distribution of normalized shear stress σ
x
y as a function of x
/
b. It is observed that the shear stress at the entire boundary is zero which conforms to the physical boundary conditions of the problem. The shear stress distribution as a function of x
/
a is antisymmetric as shown in Figure 8.Figure 6
Distribution of the normalized normal stress component,σ
y
y
/
σ
y
y
0, as a function of x
/
b at different sections of the beam.Figure 7
Distribution of the normalized normal stress component,σ
y
y
/
σ
y
y
0, as a function of y
/
a at different sections of the beam.Figure 8
Distribution of the normalized shear stress component,σ
x
y
/
σ
y
y
0, as a function of x
/
b at different sections of the beam.Figure9 shows the distribution of normalized displacement component, u
x, at the section x
/
b
=
1.0 of different composites. The axial displacement component, u
x, of glass/epoxy is the highest among those of boron/epoxy and graphite/epoxy. From Figure 10, it is observed that the displacement component, u
y, at the section x
/
b
=
1.0 is the highest for glass/epoxy and the lowest for boron epoxy, and that of graphite/epoxy remains in between them. Figures 11 and 12 describe the effect of different composites on the normalized stress components, σ
x
x and σ
y
y, at the section x
/
b
=
0.5. Effect of different composites is slightly observed on the above stress distributions.Figure 9
Distribution of the normalized displacement component,u
x, as a function of y
/
a at the section x
/
b
=
1.0 of the beam under different composites.Figure 10
Distribution of normalized displacement component,u
y, as a function of y
/
a at the section x
/
b
=
1.0 of the beam under different composites.Figure 11
Distribution of the normalized axial stress component,σ
x
x
/
σ
y
y
0, as a function of y
/
a at the section x
/
b
=
0.5 of the beam under different composites.Figure 12
Distribution of the normalized lateral stress component,σ
y
y
/
σ
y
y
0, as a function of y
/
a at the section x
/
b
=
0.5 of the beam under different composites.


### 5.2. Comparison of the Present Solution with Available Results

Three-dimensional beam is considered two-dimensional problem considering plane stress condition. The number of elements used to mesh the geometry of the simply supported beam as observed in Figure1 is 45
×
15. Figure 13 illustrates the distribution of normalized displacement component, u
x, at the section y
/
a
=
0 by the FDM and FEM methods. With increase of x
/
b ratio, the variation between the FDM and FEM solutions gradually increases. Figure 14 illustrates the distribution of normalized lateral displacement component, u
y, at the section y
/
a
=
1.0 by the FDM and FEM methods. At the region, x
/
b
=
0.5, the variation of the displacement component obtained by both of the methods is the highest and, away from the region x
/
b
=
0.5, this variation decreases. From Figure 15, it is observed that the stress component, σ
x
x, at the sectiony
/
a
=
0.5 obtained by the FDM is significantly higher than that of the FEM. But the stress component inydirection, at the sectiony
/
a
=
0.5obtained by the FD and FE methods show negligible variation as observed in Figure 16. Figure 17 illustrates the shear stress component at the section,y
/
a
=
0.5by the FD and FE methods. Here the shear stress obtained by the FD method is significantly higher than that of FEM method. The present FD method results come directly from the theoretical expressions of elasticity but the FE solutions somewhat depend on the accurateness of the assumption of the shape functions. The disadvantage of FDM is that sometimes it is difficult to handle complex boundary conditions of loaded structures.Figure 13
Distribution of the normalized axial displacement component,u
x, as a function of x
/
b at the section y
/
a
=
0 of the beam under different methods.Figure 14
Distribution of the normalized lateral displacement component,u
y, as a function of x
/
b at the section y
/
a
=
1.0 of the beam under different methods.Figure 15
Distribution of the normalized axial stress component,σ
x
x, as a function of x
/
b at the section y
/
a
=
0.5 of the beam under different methods.Figure 16
Distribution of the normalized lateral stress component,σ
y
y
/
σ
y
y
0, as a function of x
/
b at the section y
/
a
=
0.5 of the beam under different methods.Figure 17
Distribution of the normalized shear stress component,σ
x
y
/
σ
y
y
0, as a function of x
/
b at the section y
/
a
=
0.5 of the beam under different methods.


## 5.1. Analysis of Elastic Behavior of the Simply Supported Beam

The geometry of the problem is shown in Figure1. Figure 3 presents the normalized displacement components u
x as a function of x. At x
/
b
=
0.5, the axial displacement of each section is zero. The nature of the axial displacement above and below the neutral axis is opposite and antisymmetric. Towards the center (x
/
b
=
0.5) of each section, the magnitude of u
x gradually decreases to zero. On the other hand, the axial displacement component of each section with respect to x is the highest at the left and right ends of the beam, which gradually decreases to zero towards the center x
/
b
=
0.5. Figure 4 describes the normalized lateral displacement component u
y with respect to x
/
b at different sections of the beam. Middle region (x
/
b
=
0.5) of each section shows the highest lateral displacement. Figure 5 shows the deformed shape of the beam which shows that deformation in y direction is more significant.Figure 3
Distribution of the normalized displacement component,
u
x, as a function of x
/
b at different sections of the beam.Figure 4
Distribution of the normalized displacement component,u
y, as a function of x
/
b at different sections of the beam.Figure 5
Deformed shape of the cantilever beam under uniformly distributed loading.Figure6 describes the normalized axial stress component, σ
x
x, as a function of x
/
b. The stress σ
x
x is zero at the neutral axis y
/
a
=
0.5, above which σ
x
x is negative and below which σ
x
x is positive. The negative and positive stresses correspond that the regions are in compression and tension, respectively. This agrees with the physical characteristics of the problem. Figure 7 shows that σ
y
y
/
σ
y
y
0 is unity at y
/
a
=
1.0 for any value of x
/
b. This tells that the computed stress is equal to the applied stress at the top surface which verifies the boundary condition. Further, σ
y
y
/
σ
y
y
0 is maximum at the two supports of the bottom surface (y
/
a
=
0). This is due to the fact that the two supports equalize the total load applied at the top surface. Figure 8 exhibits the distribution of normalized shear stress σ
x
y as a function of x
/
b. It is observed that the shear stress at the entire boundary is zero which conforms to the physical boundary conditions of the problem. The shear stress distribution as a function of x
/
a is antisymmetric as shown in Figure 8.Figure 6
Distribution of the normalized normal stress component,σ
y
y
/
σ
y
y
0, as a function of x
/
b at different sections of the beam.Figure 7
Distribution of the normalized normal stress component,σ
y
y
/
σ
y
y
0, as a function of y
/
a at different sections of the beam.Figure 8
Distribution of the normalized shear stress component,σ
x
y
/
σ
y
y
0, as a function of x
/
b at different sections of the beam.Figure9 shows the distribution of normalized displacement component, u
x, at the section x
/
b
=
1.0 of different composites. The axial displacement component, u
x, of glass/epoxy is the highest among those of boron/epoxy and graphite/epoxy. From Figure 10, it is observed that the displacement component, u
y, at the section x
/
b
=
1.0 is the highest for glass/epoxy and the lowest for boron epoxy, and that of graphite/epoxy remains in between them. Figures 11 and 12 describe the effect of different composites on the normalized stress components, σ
x
x and σ
y
y, at the section x
/
b
=
0.5. Effect of different composites is slightly observed on the above stress distributions.Figure 9
Distribution of the normalized displacement component,u
x, as a function of y
/
a at the section x
/
b
=
1.0 of the beam under different composites.Figure 10
Distribution of normalized displacement component,u
y, as a function of y
/
a at the section x
/
b
=
1.0 of the beam under different composites.Figure 11
Distribution of the normalized axial stress component,σ
x
x
/
σ
y
y
0, as a function of y
/
a at the section x
/
b
=
0.5 of the beam under different composites.Figure 12
Distribution of the normalized lateral stress component,σ
y
y
/
σ
y
y
0, as a function of y
/
a at the section x
/
b
=
0.5 of the beam under different composites.


## 5.2. Comparison of the Present Solution with Available Results

Three-dimensional beam is considered two-dimensional problem considering plane stress condition. The number of elements used to mesh the geometry of the simply supported beam as observed in Figure1 is 45
×
15. Figure 13 illustrates the distribution of normalized displacement component, u
x, at the section y
/
a
=
0 by the FDM and FEM methods. With increase of x
/
b ratio, the variation between the FDM and FEM solutions gradually increases. Figure 14 illustrates the distribution of normalized lateral displacement component, u
y, at the section y
/
a
=
1.0 by the FDM and FEM methods. At the region, x
/
b
=
0.5, the variation of the displacement component obtained by both of the methods is the highest and, away from the region x
/
b
=
0.5, this variation decreases. From Figure 15, it is observed that the stress component, σ
x
x, at the sectiony
/
a
=
0.5 obtained by the FDM is significantly higher than that of the FEM. But the stress component inydirection, at the sectiony
/
a
=
0.5obtained by the FD and FE methods show negligible variation as observed in Figure 16. Figure 17 illustrates the shear stress component at the section,y
/
a
=
0.5by the FD and FE methods. Here the shear stress obtained by the FD method is significantly higher than that of FEM method. The present FD method results come directly from the theoretical expressions of elasticity but the FE solutions somewhat depend on the accurateness of the assumption of the shape functions. The disadvantage of FDM is that sometimes it is difficult to handle complex boundary conditions of loaded structures.Figure 13
Distribution of the normalized axial displacement component,u
x, as a function of x
/
b at the section y
/
a
=
0 of the beam under different methods.Figure 14
Distribution of the normalized lateral displacement component,u
y, as a function of x
/
b at the section y
/
a
=
1.0 of the beam under different methods.Figure 15
Distribution of the normalized axial stress component,σ
x
x, as a function of x
/
b at the section y
/
a
=
0.5 of the beam under different methods.Figure 16
Distribution of the normalized lateral stress component,σ
y
y
/
σ
y
y
0, as a function of x
/
b at the section y
/
a
=
0.5 of the beam under different methods.Figure 17
Distribution of the normalized shear stress component,σ
x
y
/
σ
y
y
0, as a function of x
/
b at the section y
/
a
=
0.5 of the beam under different methods.


## 6. Conclusions

Applying the finite difference technique on the displacement potential function formulation, a simply supported orthotropic composite beam is solved to obtain the elastic field. Simply supported composite beams are widely used in many structures. The solutions of a simply supported orthotropic composite beam using FDM technique identify the critical sections of the beam under uniformly distributed loading. Effect of different composites on the solutions is also analyzed by the present finite difference technique. The displacement and stress components at a particular section of a simply supported beam are analyzed by using FDM and FEM methods in a comparative fashion. From the comparison between the present and finite element solutions, it is observed that there are some differences in some solutions of the present problem obtained by the FDM and FEM although the trend of the curves of the solutions obtained by both of the methods is similar. The present studies will be helpful for the designing of simply supported beams made of orthotropic composite materials.



---

*Source: 961503-2014-02-13.xml*