CH4 StressAnalsis ME160
CH4 StressAnalsis ME160
CH4 StressAnalsis ME160
Chapter 4
Finite Element Analysis in Stress
Analysis of Elastic Solid Structures
Instructor
Tai-Ran Hsu, Professor
Part 1
Introduction to
Fundamentals of Theory of Linear Elasticity
What defines elastic solids?
A solid deforms in response to external actions (e.g., forces, heat, etc.). The deformation is completely reversible – meaning the
solid returns to its original shape after the removal of the external actions.
This type of elasticity occurs to solids undergoing small deformations, such as springs that exhibit linear relationships
between the applied force (F) and the induced elongation (x) that can be represented by a mathematical formula as: F = kx
where k is a constant known as the rate or spring constant. Many metallic materials fall into the category of linear elastic solids
It can also be stated as a linear relationship between stress (σ) and strain (ε) in stretching or compressing a thin rod by
The expression: σ = Eε where E Is known as the elastic modulus or Young's modulus
This type of solids behaves as elastic materials as described above, but can exhibit large deformations, such as rubbers and polymers
The FE formulation presented in this course will be based on linear elasticity theory
Linear Elastic Behavior of Solids
Fundamental assumptions
(1) The material is treated as a continuous medium (or a continuum). In the words, the material is homogeneous with no
internal defects or voids of significant sizes
(2) The material is isotropic – meaning its properties are uniform in all directions
(3) The material has no memory
(4) The material exhibits the same properties in tension and compression
L
ε =
0 1 n
∫
L0 L
= n
L0
(3) Relationship between ε and e: ε = ℓn(1+e)
True
Designations:
A’ = proportional limit
A = elastic limit
B = yield point
Stresses
m =necking point
Engineering
f = rupture point
S0 = Yield strength of material
Su = ultimate tensile strength of material
0 Strains
RESPONDES
To Allied
Foreces:
2 things will
happen:
(Element) Displacements:
{U}T: {Ux,Uy,Uz}
in a deformed solid
A small element
located at x,y,z:
Result in:
Because the forces applied to a general 3-D solid, the induced stresses are MULTI-directional designated by σαβ:
σ = magnitude, subscript α = the axis normal to the plane of action, subscript β = the direction the stress component points to
So, stress component σyy = stress component acting on the plane normal to the y-axis and pointing to the y-direction, whereas
stress component σxy = stress components acting on the face normal to the x-axis but points to the y-direction.
σ 1
In theory, there are nine (9) But in reality there σ
stress components are only six (6)
For FE formulation: 2
σ xx σ xy σ xz σ xx σ xy σ xz σ1=σxx,σ2=σyy,, σ 3 (4.1)
everywhere
independent stress
{σ } =
[ ]
inside the solid: σ = σ yx σ yy σ yz components with: [σ ] = σ yy σ yz σ3=σzz, σ4=σxy,
σ5=σyz, σ6=σxz σ 4
σxy =σyx , σxz =σzx , SYM σ zz σ 5
σ zx
σ zy σ zz
σyz = σzx
σ 6
Induced Stresses in Deformable Solids subjected to External Forces – Cont’d
z Stress components with same subscripts such as:
σzz σxx, σyy, σzz are called NORMAL stress components,
D C with unit: psi, or Pascal (Pa) = N/m2 .
σzy
σzx σyz Stress components with different subscripts such
σyy
A B as: σxy, σxz, σyx, σyz, σzx, σzy are SHEARING stress
σxz σyx
σxy components. Shear stress has the unit of change of
G
σxx y angle from the original right angle to the angle with
the stress, i.e (π/2) – θ. The unit is thus “angle
change” in radians (rad)
E
F
x
Effects of Normal and Shear Stresses to Solid Deformation
Normal stresses change size: Shear stresses change shape:
y y
contraction Exist NORMAL
Exist on the SURFACES
to surfaces
-σxx σxx
-σxy σxy
elongation
x x
Relationships between Primary unknown and Secondary Unknowns in FE analyses
Primary Unknown:
Element displacements: {U(x,y,z)} and
Nodal displacements: [N(x,y,z)]{u}T
in FE formulation: {U} = [N]{u}
where {N(x,y,z)} = Interpolation function
Strain-displacement relations
F 0 F
x
Poisson’s ratio defines as:
Lateral contraction strain ε yy
Longitudinal stretching: ν = =
+ εxx Longitudinal stretching strain ε xx
● Poisson ratio exists in most multi-axially loaded elastic solids, and this effect needs to be included in all
FE formulations. This effect also leads to the use of Generalize Hooke’s Law for the formulations
Generalized Hooke’s Law for Solids with Multi-axial Stresses:
Elongation in one-direction causes contractions
Stress vs. strain relationship for deformed solids with 3-D deformation in other directions, and vise versa.
in x-direction
in y-direction
in z-direction
The following expression express the stresses in terms of strains – the Hooke’s law:
One may derive the uniaxial stress situation as a special case from the above expression to obtain:
σxx = Eεx
by substituting: σyy = σzz = 0, and εy = -γεx and εz = -γεx in the generalized Hooke’s Law
Element Strain-Displacement Relations for FE Formulation
There are six (6) strain components corresponding to the stress components in interior of the deformed solid.
These strain components are related to the displacements of the solid induced by the external forces.
Because deformation of the solid CONTINUOUSLY varying throughout the solid, the following relationship exits:
By theory of elasticity: ∂ Element
∂x 0 0 or: {ε(x,y,z)} = [D]{U(x,y,z)} (4.3)
Displacements
∂
0 U(x,y,z) ∂
ε xx 0 ∂x 0 0
ε ∂y
yy ∂ U (x, y, z ) ∂
ε zz 0 0 x 0 0
∂z U (x, y, z ) (4.2)
∂y
=∂ ∂ y
ε xy
0 U (x, y, z ) ∂
ε yz ∂y ∂x z 0 0
∂ ∂ [D] = ∂ ∂z
ε xz 0 ∂ (4.4)
0
∂z ∂y
∂ ∂ ∂y ∂x
Element Strains
∂z
0
∂x
∂ ∂
Corresponding to 0
Element stresses ∂z ∂y
∂ 0
∂
Example: For uni-axial elongation or contraction of a rod:
∂U x ( x ) dU x ( x )
∂z ∂x
Element displacement {U} = {Ux(x)}, The corresponding strain in element is: ε xx = =
∂x dx
Element Stress-Strain Relations – the Generalized Hooke’s Law
According to generalized Hooke’s law for MULTI-Axial stress state,
the following relationship between the element stress and strain exists:
Displacements:
{Φ(x,y,z)}T: {Φx(x,y,z,Φy(x,y,z,Φz(x,y,z}
We will select tetrahedron elements as the basis for FE formulation for general 3-D solid structures
[K ] = ∑ [K em ]
m
where m=total number of elements in the FE model
FE Model for s Structure made of 1
Tetrahedron Elements
Derivation of Element Equations-Cont’d Principle of deriving element equation using Rayleigh-Ritz variational principle
From Chapter 2: Let us determine a suitable “functional” to derive the element equation.
∂{φ } ∂{φ }
A general form of functional: χ (φ ) = f {
∫v ∂rφ }, , .......... dv + ∫ g {φ }, ,......... . ds v = volume, s = surface
s
∂r
and then apply the Variational principle on:
∂χ
∂φ
1 from which equations of each element are derived:
∂χ
∂χ (φ ) ∂{χ } ∂{χ } ∂{χ }
= ∂φ2 = 0 = 0, = 0, = 0,..............
∂{φ } • ∂φ1 ∂φ2 ∂φ3
•
•
The “functional” for a deformed solid subjected to external forces is “POTENTIAL ENERGY” (P) in the situation.
The potential energy associated with a deformed solid can be defined as:
P=U-W (4.8)
where U = the strain energy in a deformed solid, and W = the work done to the deformed solid by external forces acting on
the solid body in the volume and surface of the solid
Derivation of Element Equations-Cont’d Strain energy in a deformed solid:
As we mentioned early in this Chapter that a solid in (a) deforms
into a new shape in (c) – but not indefinitely.
● It stopped further deformation after deformed by certain
(a) (C) amount.
● It reaches a new state of equilibrium. WHY???
Imagine the following phenomenon:
Stretch a free-hung spring by a weight W.
The spring will elongate, but only by a finite amount.
(b) Ask yourself: WHY?
Answer: the elongation of the spring develops a “resistance,”
which increases as the spring elongates. The spring ceases
further elongation when the “resistance” in the spring balances
the applied weight (force). We say the spring – and the applied
weight reaches a new state of equilibrium, which stop the spring
(d) From further elongation.
Next: what will happen to the spring after the weight is removed?
You will say that the spring returns to its original length, but WHY??
Answer: because a form of ENERGY was stored in the stretched
spring when it is elongated. This ENBERGY is released to restore
Now, you know why the solid in (a) ceases to deform further after the spring to its original shape after the external force (the weight)
the application of the system of external forces {p} has been applied was removed.
to the solid in (b). And you would know that there is such ENERGY associated with the solid deformation developed in the solid called
“STRAIN ENERGY.” which is responsible for restoring the solid to its original shape after the applied forces are removed in “ELASTIC” solids.
U = ∫ {ε } {σ }dv
Mathematical expression of strain energy in State (c) is: 1 T
(3.9) textbook
2 v
Derivation of Element Equations-Cont’d Potential energy in a deformed solid subjected to external forces:
Strain energy:
σ xx
σ
yy
σ zz
U = ∫ {ε } {σ }dv = ∫ {ε xx ε yy ε zz ε xy ε yz ε xz } dv
1 T 1
2 v 2 v σ xy
σ yz
σ xz
= ∫ (ε xxσ xx + ε yyσ yy + ε zzσ zz + ε xyσ xy + ε yzσ yz + ε xzσ xz )dv
1
(4.9)
2 v
Both the strain and stress components are function of (x,y,z), and dv = (dx)(dy)(dz) = the volume of given points
in the deformed solid.
Work done to deform the solid: Definition of “work”: Work (W) = Force x Displacement (deformation)
Two kinds of forces: (1) body forces (uniformly distributed throughout the volume of the solid (v)), e.g., the weight,
(2) surface tractions, e.g., the pressure or concentrated forces acting on the boundary surface (s)
fx
= ∫ {φ x (x, y , z ) φ y (x, y , z ) φ z (x, y, z )} f y dv
f
v
z
t x
+ ∫ {φ x (x, y, z ) φ y (x, y, z ) φ z (x, y, z )}t y ds
s
t (4.10)
z
where {Φ(x,y,z)} = the displacement of the solid at (x,y,z), {f} = body forces, and {t} = the surface tractions , and
ds = the part of the surface boundary on which the surface tractions apply
Derivation of Element Equations – cont’d Potential energy in a deformed solid subjected to external forces:
P = U − W = ∫ {ε ( x, y, z )} {σ ( x, y, z )}dv −
1
2 v
T
(∫ {φ (x, y, z )} { f }dv + ∫ {φ (x, y, z )} {t}ds )
v
T
s
T
(4.11)
Following the Rayleigh-Ritz Variational principle, the equilibrium condition for the deformed solid should
satisfy the following conditions:
∂P(φ )
= 0
∂{φ }
What we had formulated was for continuum solids. We will now derive the ELEMENT EQUATION for discretized solids in FE mesh:
We need to make distinction between the ELEMENT quantities and the NODAL quantities.
The primary quantity in FE analysis is DISPLACEMENTS. We need teo make distinction between the Element displacements and
the Nodal displacements.
The element displacement is: {Φ(x,y,z)} with three components:
Φx(x,y,z) = the element displacement component along the x-direction
Φy (x,y,z) = the element displacement component along the y-direction, and
Φz(x,y,z) = the element displacement component along the z-direction
Derivation of Element Equations – cont’d for FE mesh of discretized 3-D solids with tetrahedron elements
We realize that TETRAHEDRON and HEXAHEDRN elements are used in FE models for general 3-D solid structures.
The tetrahedron elements are the “basic elements” for this type of structures, because hexahedron elements
are made up by 4 or more tetrahedron elements.
Our FE formulation for general 3-D solid structures will thus be based on TETRAHEDRON elements
We notice that tetrahedron elements has four (4) associate nodes: Φ1, Φ2, Φ3, and Φ4 with FIXED (specified) COORDINATES.
Each node has three (3) displacement components too. These nodal displacement components are:
{φ }T = {φ1x φ1 y φ1z φ2 x φ2 y φ2 z φ3 x φ3 y φ3 z φ4 x φ4 y φ4 z }T
where Φ1x, Φ1y, Φ1z = displacements in Node 1 in 3 directions; Φ2x, Φ2y, Φ2z = displacements in Node 2 in 3 directions;
Φ3x, Φ3y, Φ3z = displacements in Node 3 in 3 directions; Φ4x, Φ4y, Φ4z = displacements in Node 4 in 3 directions
Derivation of Element Equations – cont’d for FE mesh of discretized 3-D solids with tetrahedron elements
We mentioned previously that the functional that we will use to derive the element equations in FE formulation of solid structures
Is the potential function in the solid as show below:
P = U −W =
1
2 ∫v
{ε ( x , y , z )}
T
{σ (x, y, z )}dv − (∫ {φ (x, y, z )} { f }dv + ∫ {φ (x, y, z )} {t}ds )
v
T
s
T
(4.11)
Now because the ELEMENTS in discretized solid) Also, these elements are interconnected by the NODES to simulate the original
solid structures. This “link” requires the FE formulation involves the functional with both the ELEMENT and NODAL quantities in
the formulation.
This “link” is established using the INTERPOLATION FUNCTION that relates the
element quantities with the corresponding nodal quantities such s:
{φ }T = {φ1x φ1 y φ1z φ2 x φ2 y φ2 z φ3 x φ3 y φ3 z φ4 x φ4 y φ4 z }T
Derivation of Element Equations – cont’d Key equations to construct the functional - the potential energy
v 2
Typical tetrahedron element
for 3-D FE models
Hence U=
1
∫ ([B(x, y, z )]{φ })T [C ]([B(x, y, z ){φ }])dv (4.15)
2 v
1
{ }
T
[ ] [C ][B( x, y, z )]{φ }dv (4.16)
2 ∫v
or U= φ B ( x , y , z )
T
Derivation of Element Equations – cont’d The functional for Variational process
We mentioned that the functional for deriving the element equations for
discretized solid structure is the POTENTIAL ENERGY (P) as shown below:
P = U −W =
1
2 ∫v
{ε ( x , y , z )}
T
{σ (x, y, z )}dv − (∫ {φ (x, y, z )} { f }dv + ∫ {φ (x, y, z )} {t}ds )
v
T
s
T
(4.11)
By substituting the Strain energy expressed in Equation (4.16) into the above equation, we get:
1
{ } [ ] [C ][B( x, y, z )]{φ }dv
2 ∫v
P (φ ) = φ T T
B ( x , y , z )
Due to the fact that nodal displacement {Φ} have “fixed value” but not a function og (x,y,z), so they can be factore
Out of the integration with respect to (x,y,z). We thus have the following:
T1
P (φ ) = {φ } ∫ {φ } [B ( x, y, z )] [C ][B ( x, y, z )]dv {φ }
T T
2 v
( v
T
)
− {φ } ∫ [N ( x, y, z )] { f }dv − {φ } ∫ [N ( x, y, z )] {t}ds
T T
s
T
(4.18)
Derivation of Element Equations – cont’d Element equation by Variational process
By substituting the potential energy P in Equation (4.18) into the above equation:
T1
{φ } ∫ {φ } [B ( x, y, z )] [C ][B ( x, y, z )]dv {φ }
T T
∂P(φ ) ∂ 2 v
=0
=
∂φ ∂φ
T
( v
T
) s
T
− {φ } ∫ [N ( x, y, z )] { f }dv − {φ } ∫ [N ( x, y, z )] {t}ds
T
Derivation of Element Equations – cont’d Element equation by Variational process
NOTE: In FE stress analysis of solid structures, it is customer to represent the element displacements by:
U x ( x, y, z )
{U ( x, y, z )} = U y (x, y, z )
U ( x, y, z )
and nodal displacements by: {u}. z
The relationship between the element displacements and the nodal displacements are:
{U} = [N(x,y,z)] {u}
where [N(x,y,z)] = the Interpolation function. It is a row matrix for 1-D bar elements,
rectangular matrices for 2- or 3-D elements.
● Interpolation function enables the determination of the primary quantities in the element with specified coordinate (x,y,z)
with the same primary quantities of the associate nodes
Part 3
Finite Element Formulation for
One-dimensional Bar elements
1. Bar elements subjected to unidirectional load
[A] = 1 and
1 x 1 x1
x1 = 0 x2 = L A = = x2 − x1
1 x2 1 x2
Figure 1.6 Linear Interpolation Function
of a Bar Element 1 x2 − x1
The inverse matrix of [A] is: [A]−1 = = [h]
x2 − x1 − 1 1
By following Equation (1.15), we have the interpolation function of a bar element to be:
1 x2 x1 1
N ( x ) = {R} [h] = {1 x} − 1 1 = L {( x2 − x ) ( x1 + x )}
T
x −
2 1x
For the present case with x1 = 0, and x2 = L, we have the interpolation function to be:
x x
N ( x ) = 1 − (1.7)
L L
Derive the element equation for a typical 1-D bar element: a bar element made of one material with Young’s modulus E.
The bar is subjected two forces F1 and F2 as shown in the figure below. Establish the element equation
u1 Node 1 Node 2 u2 Because the bar is made of one material, and the applied forces are along
F1 ● ● F2
x the length of the bar. So the bar will only deform along the length of the bar
L e.g., the element displacement U is a function of coordinate x only:
x=x1 x=x2
We thus have the element displacement to be U = U(x). The corresponding displacements at the two nodes are:
{u}T = {u1 u2}.
The relationship between the element displacement and the corresponding nodal displacements are:
{U(x)} = {N(x)} {u} (4.24)
in which {N(x)} is the interpolation function for one-D bar element. It is a row matrix with two columns
We had derived this interpolation function before with assumption that the element displacement {U(x)} follows
a linear polynomial function: U(x)= α1 + α2x, and found the interpolation function using this liner function to be
In the same form as in Equation (1.7) we derived using general formulation method:
{N (x )} = 1 − x x
L L (4.25)
with L = x2 – x1
Because the bar element is subjected to force along the x-coordinate only, and both the induced stress and strain are
in the direction of x-coordinate only. We thus obtain from Equation (4.4) and obtain:
∂
[D] = =
d
∂x dx
From Equation (4.13), we thus have the [B] matrix to be:
For the same reason of being uniaxial stress distribution, i.e., σxx = Eεxx by uniaxial Hooke’s law, we have the
elasticity matrix [C] = E – the Young’s modulus from Equation (4.7) with Poisson’s ratio γ = 0. We thus has the
element equation of the deformed bar expressed in the form:
[Ke] {Φ} = {q} (4.26)
x2 1 T 1 EA x2 − 1
= ∫ {− 1 1} (E ) {− 1 1} Adx = 2 ∫ {− 1 1}dx
x1
L L L x1 1
EA x2 1 − 1 EA ( x2 − x1 ) − ( x2 − x1 ) EA 1 − 1
= 2 ∫ dx = 2 =
L x1 − 1 1 L − ( x2 − x1 ) (x2 − x1 ) L − 1 1 (4.27)
Example 4-1 Show the element equation for the bar element in the figure, and determine displacement at Node 2
The nodal forces {q} for the bar element are expressed as: {q}T = {-F1 F2}
u1 Node 1 Node 2 u2
F1 ● ● F2
x
We may thus express the element equation for the bar element to be:
L
EA 1 − 1 u1 − F1
= (a)
x=x1 x=x2 L − 1 1 u2 F2
Because there is only one bar element in the structure, the overall stiffness equation in Equation (1.28)
is identical to that of element equation in Equation (4.24).
However, adjustments of the now overall stiffness equation with the three matrices in Equation (4.24) would be required,
as expressed:
K aa K ab qa Ra
K =
ba K bb qb Rb
where {qa} = specified (known) nodal quantities; {Rb} = specified (known) applied resulting actions, from which we may obtain:
The specified nodal unknown {qa} in this case is u1 = 0 and F1=0, the required unknown {qb} is u2, We thus need no adjustment of the
Matrices in the overall stiffness equation in (4.24). The required unknown quantity u2 is obtained from Equation (6.12, REF) to be:
u2 =
EA 1
(F2 − (− 1)(0)) = EA F2
L (1) L
Example 4-2 Deformation and stress in a compound bar made of two different materials
Use the FEM to determine the displacements at the joint of a compound road made of copper and aluminum induced by a
uniaxial force P = 30000 N at of the end of the rod as shown in the Figure A below:
Our solution begins with developing the “element equations” for both elements in the FE model:
Due to the fact that the present case involves TWO elements with Node 2 being common to both these element,
we need to assemble the coefficient matrices by following the established rule by summing up the values of Node 2
from both element coefficient matrices.
Example 4-2 – Cont’d
[ ]
We thus assemble the overall stiffness matrix by adding K e1 In Equation (a) and [K e2 ] in Equation (b) in the following way:
The numbers in boldface in Equation (c) are those associated with Node 2.
The overall stiffness equation for the bar structure with specified loading/boundary conditions is thus expressed in the
following partitioned matrices as:
7.317 − 7.317 0 u1 = 0 p1 = 0
(d)
10 6 0 154.367 − 147.05 u2 = p2 = 0
0 0 147.05 u3 p3 = 30000
The two unknown nodal displacements u2 and u3 may be obtained by the following equations using the partitions
in Equation (d):
154.367 − 147.05 u2 0
106 =
0 147 . 05 3
u 30000
resulting in the displacement of the compound bar at the joint (Node 2) to be u2 = 1.94 mm and the displacement
at the free end u3 = 2.04 mm. The total elongation of the rod is 3.98 mm.
Example 4-2 – Cont’d
Strain in Elements
Now that we have solved the displacements at the 3 nodes, we may use the train-displacement relations to
determine the induced strains in both these elements:
{ε }T {
= ε 1xx ε xx2 }
The strain-displacement relationship is available from the expression: {ε} = [B]{u} in Equation (4.3), in which the matrix {B} is:
x − x1 1 1
[B(x )] = d
{N (x )} = d x − x2 − = −
dx dx x1 − x2 x1 − x2 x1 − x2 x1 − x2
We have: Node 1 at x1 = 0; Node 2 at x2 = 915 mm, and Node 3 at x3 = 1220 mm, leading to: the length of
Element 1 = L1 = 915 mm; the length of Element 2 = L2 = 305 mm. We may thus express the [B(x)] for both
elements to be:
[B1 ] = −
1 1 1
, and [B2 ] = −
1
915 915 305 305
Example 4-2 – Cont’d
From the [B] matrices for both elements, we may compute the strains in Element 1 and 2 as follows:
1 1 u1 = 0 1 1.94
ε 1xx = − = u2 = = 0.21% for element 1, and
915 915 u2 915 915
Stresses in elements
We may use the Hooke’s Law to determine the stresses in each of these two elements from their corresponding strains
For element 1 with Ecu = 10300 MPa:
Common Trusses
FE Formulation of Truss Elements Using 1-D Bar Elements - Cont’d
The characteristics of a truss element can be summarized as follows (a quote from Dr. Agarwal’s
lecture notes):
2. It is a two-force member i.e. it can only support an axial load and cannot support a
bending load. Members are joined by pins (no translation at the constrained node,
but free to rotate in any direction). Being a two-force member, the force acting in the member is in the
length-direction only.
3. The cross-sectional dimensions and elastic properties of each member are constant along its length.
5. The bar elements for truss structures is mechanically equivalent to a spring, since it has
no stiffness against applied loads except those acting along the axis (the length direction) of the member.
6. However, unlike a spring element, discussed in previous chapters, a truss element can be oriented in any
direction in a plane, and the same element is capable of tension as well as compression along its longitudinal
directions.
FE Formulation of Truss Elements Using 1-D Bar Elements – Cont’d
Simple Bridge Structure:
FE model with 1-D bar elements with designated element and node numbers. All elements are
interconnected by “frictionless hinges”:
3 5
ο (6)
ο
Element (1) (5) (11)
(7)
(3) (9)
Node 1
ο (2)
ο2 (4)
4
ο (8)
ο6 (10) ο7
ο ο οο
P1 P2 P3
FE Formulation of Truss Elements Using 1-D Bar Elements – Cont’d
y
Truss members along x-direction of y-direction only:
The or these element equations for these elements will be in similar forms
x as for the 1-D bar elements in Examples 4-1 and 4-2.
The element equation for those truss members that are inclined with the x-coordinate with angle θ,
However, will be in different forms, as will be derived in the following procedures.
FE Formulation of Inclined Truss Elements Using 1-D Bar Elements
Inclined truss bar element in x-y plane: Element (1), (5), (7) and (11):
y
vj Fj
j
uj
(xj,yj)
vi
The inclined truss bar elements such as shown in the figure at right
may have TWO displacement components, but these elements can only i ui
elongate or contract in the LONGITUDINAL direction only.
(xi,yi)
Fi
The same applies to the induced strain and stress.
x
0
So, These are regarded as “a special bar elements.”
These bar elements remain having two nodes: Node i and Node j
Located at specified coordinates (xi, yi) and (xj, yj) respectively.
The displacement components in an inclined bar element in the plane defined by the (x,y) coordinate system in the figure below do not
contribute in producing the displacement, strain and stress in the bar element for this reason . However their “equivalence” to the ones
along the length direction of the bar element do.
Consequently, we need to convert the current situation in the inclines bar element in the left to the equivalent situation in the right
of the figure through a coordinate transformation process of the nodal displacement components and the applied nodal forces
Converting
“truss bar element”
to equivalent “typical
bar element”
By referring to the above figure, we have the following relationship in transforming the nodal displacement components from
The global coordinates (x,y) in the left figure to the local coordinate (x) in the right of the figure:
Transformation of nodal forces can be done in a similar way as we did with nodal displacements. However, it would be easier to do it using
the work done to the element by these forces, as work done is a scalar quantity, which is easier to transform in space than the vector
quantities such as displacements.
Since the work done to the element by nodal forces may be expressed as: W = {δ }T { f } = {u}T {F } for the same element in both local and global
coordinate systems, or in a long-hand form:
Fix
F
f1
W = {δ1 δ 2 } = {ui vi u j v j } or {δ } { f } = {u} {F }
iy T T
f2 F jx
F jy
Substituting {δ} = [T]{u} in Equation (26b) into the above expression, we get: ([T]{u}])T{f} = {u}T{F}, leading to:{u}T[T]T{f} = {u}T{F}
We will thus have the nodal force transformation by the following relationship:
[T]T {f} = {F} (4.30)
where {f} = nodal forces in “typical bar element”, {F} = nodal forces in “inclined truss element”
FE Formulation of Truss Elements Using 1-D Bar Elements – Cont’d
Derivation of element equation – Cont’d
Element equation in Local Coordinate System
Now that we have made the real situation in the “Global” coordinates (x,y) to be equal to the situation in 1-D “Local” coordinate (x) situation
through a “transformation” process:
=
Matrix:
[T ] =
s 0 0
c
0 0 c s
(4.26a) Previously derived interpolation function:
x x
{N (x )} = 1 −
L L
Previously derived “element equation”
Real situation in Coordinate Transform: with modified notation of nodal quantities:
“Global” coordinates
EA 1 − 1 δ i − f i
= (4.24)
L − 1 1 δ
j j
f
AE − AE
or [K e ]{δ } = { f } with [K e ] = L L (4.27)
−
L
AE AE
L
FE Formulation of Truss Elements Using 1-D Bar Elements – Cont’d
Derivation of element equation – Cont’d
Element equation in Local Coordinate System
EA 1 − 1 δ i − f i
= (4.24)
L − 1 1 δ j f j
AE − AE
[K e ]{δ } = { f } with [K e ] = − AE L L (4.27)
L
AE
L
Relationship of nodal forces between the two coordinate systems: {F} in global coordinate system = [T]T{f} in Local coordinate system .
The element equation in global coordinate system thus have the form:
[K e ]{δ } = { f } (4.31)
The nodal forces in the Local coordinate system {f} can be viewed as a transformed forces {F} rom the global coordinate systems with
the relationship: [T]T {f} = {F} as shown in Equation (4.27), or {f} = ([T]T)-1, leading to:
[Ke]{δ} = ([T]T)-1 {F}
But we also have already derived the following relationships: {δ} = [T]{u} in Equation (4/26b), and [T]T = {F} in Equation (4.27).
By substituting these relationships into the above expression, we will get the following “reversed transformation” of element equation
from the Local coordinate system to Global coordinate systems:
[K ][T]{u} = ([T]T)-1 {F} (4.32)
e
FE Formulation of Truss Elements Using 1-D Bar Elements – Cont’d
Derivation of element equation – Cont’d
Element equation in Local Coordinate System
The Element Equation of inclined truss elements is: [Ke][T]{u} = ([T]T)-1 {F} (4.33)
We realize the fact that the transformation matrix: [T]T = [T]-1 (Reference No. 2). Thus Equation (4.33) has a new form of:
[Ke][T]{u} = ([T]-1)-1 {F} leads to: [Ke][T]{u} = [T]{F} with ([T]-1)-1 = [T] in the last expression.
Now, if we multiply both sides by [T]-1, and use the relationship of [T]-1 = [T]T, we have the element equation of
the inclined truss element (i.e., the bar element in the global coordinate system) to be:
[T]T[Ke][T]{u} = {F} (4.34)
ui Fix
v F
where {u} = i and
iy
{F } = (4.36a)
u j F jx (4.36b)
v j F jy
c2 cs − c 2 − cs
AE cs s 2 − cs − s 2
The stiffness matrix: [K e ] = 2 (4.37)
L − c − cs c
2
cs
2
− cs − s s
2
cs
1) This stiffness matrix represents the stiffness of a single element with incline angle θ
2) It is symmetric about the diagonal line of the square matrix
3) Since there 4 unknown nodal displacements - meaning 4 degree-of-freedom (dof), the matrix is of the size
of (4x4)
4) The terms c and s represent the cosine and sine values of the orientation of the element with the
horizontal plane, rotated in a counter clockwise direction (+ve direction)
Example 4-3 FE Stress Analysis of a Truss Structure:
Solution:
We realize the fact that there are 3 members, each has its own dimension, and material with properties listed in the above table.
So, we may conveniently construct the FE model of the truss structure by using 3 elements 1 with node pair 1-3, element 2
with node pair 2-1, and element 3 with node pair 2-3, as shown in the figure.
We will derive the element equations for all these 3 elements by using the formulations for truss elements with: Equation (4-24) for
the horizontal Element 1, and Equations (4.35) and (4.37) for the inclined elements 2 and 3.
Example 4-3 FE Stress Analysis of a Truss Structure - cont’d
Derivation of element equations
v1y v3y Being horizontal, we have the incline angle θ = 0, which leads to:
c = cosθ=1, c2 = 1, s=sinθ = 0, s2 = 0 and cs = 0
u1x (1)
u3x We also calculate the coefficient of the stiffness matrix of element 1 as follows:
Node 1 Node 3
L1 = 0.026 m ( )(
E1 A1 70000 ×106 200 ×10 −6
=
)
= 53.84 ×106 N m
L1 0.26
c2 cs − c 2 − cs 1 0 −1 0
2 0 0
[ ] − −
2
cs s cs s 0 0
K e1 = 1 1 × 2 = 53.84 ×106
E A
(a)
L1 − c − cs c 2 cs − 1 0 1 0
2
− cs − s 2
cs s 0 0 0 0
Example 4-3 FE Stress Analysis of a Truss Structure - cont’d
Derivation of element equations
v1y
In this Element 2, we have the inclined angle θ=90, which lead to:
Node 1 u1x
c = cos 90o = 0, c2 = 0; s = sin 90o = cos0o = 1, s2 = 1 and cs = 0 ,
L2=150 mm
and
(2) ( )(
E2 A2 70000 × 106 200 × 10 −6
=
)
= 31.06 × 106 N m
L2 0.15
0 − 1 0 1
Example 4-3 FE Stress Analysis of a Truss Structure - cont’d
Derivation of element equations
u2x X =
( )(
E3 A3 200000 ×106 100 ×10 −6 )
= 66.67 ×106 N m
Node 2
L3 0.3
From Equation (4.21): [Ke] {Φ} = {q} with a general expression for element equations, we are now in the position to express the
same equations for the 3 elements in the current truss structure as follows.
The 3 Element equations for the truss structure:
Element (1) with Node 1 and 3 Element (2) with Node 1 and 2
By following the description on “assembly of element stiffness matrices” in Step 5 with diagram of Chapter 3 Steps in Finite
Element Analysis, we may assemble the three (3) truss element stiffness matrices shown above in the following form:
(d)
With this overall stiffness matrix, we may establish the overall stiffness equation for the truss structure as shown below:
(e)
Example 4-3 FE Stress Analysis of a Truss Structure - cont’d
Boundary and Loading conditions of the Truss Structure
We recognize the following specified conditions for the discretized truss structures:
The overall stiffness equation in Equation (e) after the substitution of the above specified boundary and loading conditions has the form:
We will get the 3 unknown nodal displacements from the portioned over
stiffness equation by the expression: {qb} = [Kbb]-1 {Rb}:
−1
v2 y 48.27 − 28.87 − 16.67 0 3.1646 867.36 ×10 −19 3.1646 0
−6 −8 −19
u3 x = 10 − 28.87 103.84 28.87 0 = 10 867.36 ×10 1.86 − 3.2166 0
v − 16.67 28.87 16.67 − 400 3.1646 − 3,2166 14.734 − 400
3y
u2y = -12.65x10-6 m, u3x = 12.86x10-6 m and v3y = -58.94x10-6 m with negative signs meaning downward direction
Example 4-3 FE Stress Analysis of a Truss Structure - cont’d
Solution of secondary unknowns of induced strains and stresses in all three elements in the truss structure
We will determine the strain components and then stress components in each of the 3 elements
in the truss structure.
We should bear in mind that truss members are “two force members,” in which only the in-line
displacement components will produce strains and stresses.
We should bear in mind that truss members are “two force members,” in which only the
in-line displacement components will produce strains and stresses.
with θ=0o
X=0 X=L1
1 0 0 0
The transformation matrix: [T ] =
0 0 1 0 ui
δ i 1 0 0 0 vi
The relationship of nodal displacements in Local coordinate system {δ } = =
and the corresponding global coordinate systems is: {δ} = [T]{u}, or: δ j 0 0 1 0 u j
But we obtained the nodal displacements from the previous calculations to be: v j
ui = u1x =0, vi = v1x =0, uj = u3x = 12.86x10—6m, vj = v3y = -58.94x10-6m
Example 4-3 FE Stress Analysis of a Truss Structure - cont’d
Solution of secondary unknowns of induced strains and stresses in all three elements in the truss structure
0
δ i 1 0 0 0 0
0
{δ } = = −6
= −6
δ j 0 0 1 0 12.86 ×10 12.86 ×10
− 58.94 ×10 −6
We have the strain in the element is obtained by using Equation (4.12): {ϵ} = [B(x)]{Φ}, and [B(x)] = [D]{N(x)] (3.13)
We thus have:
[B] = [D]{N (x )} = d 1 − x x = − 1 1
dx L1 L1 L1 L1
The strain in Element 1 is thus:
1 δ i 1 1 0
{ε xx } = − 1 = − −6
= 49.46 ×10 −6 m / m
L1 L1 δ j 0.26 0.26 12.86 ×10
The corresponding stress in Element 1 can be obtained from Equation (4.6): {σ}=[C]{ϵ}=(70000x106)(49.46x10-6)=3.43 MPa
Example 4-3 FE Stress Analysis of a Truss Structure - cont’d
Solution of secondary unknowns of induced strains and stresses in all three elements in the truss structure - Cont’d
Element 2: with Node 1 and Node 2
We have the strain in the element is obtained by using Equation (4.12): {ϵ} = [B(x)]{Φ}, and [B(x)] = [D]{N(x)] (3.13)
v2 d x 1 1
We thus have: [B ] = [D ]{N (x )} = 1 −
x
= −
dx L2 L2 L2 L2
1 δ i 1 1 − 12 ×10 −6
The strain in Element 2 is thus: {ε xx } = − 1 = −
−6
= 80 ×10 m / m
L2 L2 δ j 0.15 0.15 0
The corresponding stress in Element 2 can be obtained from Equation (4.6): {σ}=[C]{ϵ}=(70000x106)(80x10-6)=5.6 MPa
Example 4-3 FE Stress Analysis of a Truss Structure - cont’d
Solution of secondary unknowns of induced strains and stresses in all three elements in the truss structure - Cont’d
Element 3: with Node 2 and Node 3
v3j = -58.94x10-6
We obtained the nodal y
Node 3 u3j=12.86x10-6
Displacements from
The previous calculation v2i = -12.65x10-6
to be as shown
Node 2 x
u2i = 0
δi=u2x δj=u3x
Transformation matrix ο ο x
L3=0.3 m
[T]
The transformation matrix [T] has the following form with θ = 30o:
The nodal displacements in the equivalent Element 3 in the Local coordinates can be related to the real Element 3 in the global
coordinates by the transformation matrix to be:
ui = u 2 i = 0
−6
δ i 0.866 0.5 0 0 vi = v2i = −12.65 ×10 6.325 ×10 −6 6.325 −6
{δ } = = u = u = 12.86 ×10 −6 = −6
= x10
δ
j 0 0 0 . 866 0 .5 j 3j (11 . 1368 − 29 . 47 )10 − 18 . 3332
v j = v3 j = −58.94 ×10
−6
x
But since the interpolation function of the equivalent Element 3 in the local coordinate is: {N (x )} = 1 − x
(1.7)
L3 L3
d x 1 1
And the matrix [B(x)] = [D][N(x,y,z)] from Equation (4.13) in the present case to be: [B(x )] = 1 −
x
= −
dx L3 L3 L3 L3
Example 4-3 FE Stress Analysis of a Truss Structure - cont’d
Solution of secondary unknowns of induced strains and stresses in all three elements in the truss structure - Cont’d
Element 3: with Node 2 and Node 3
The strain {ε} in terms of nodal displacements in Element 3 may be obtained by using Equation (4.12): {ε}=[D]{N(x)}:
1 1 δ i 1 1 6.325
{ε } = ε xx = − = − × 10 −6
{σ} = σxx = E3εxx = (200000x106)(-82.194x10-6) = - 16438800.33 N/m2 or -16.44 MPa (a compressive stress)
Finite Element Formulation of Elastic Solid Structures
Induced
y Original straight beam: deflection by y P – concentrated force
applied forces W(x) – distributed load
per unit length
x x
x Deflected state
Induced deflection: y(x)
3) Beams can also resist the applied moments that cause the beams to bend, i.e. beams can have deflection
induced by applied moments (e.g., torque in the x-y plane)
4) Beams may be straight or curved. We will focus on straight beams only .
5) The application of unilateral forces or moments will not only cause the beam to bend in shapes, but also
having deflections, y(x) perpendicular to the beam axes, as shown in the above figure.
6) So, the actions to beam can be: concentrated forces, distributed load and moments
The induced reactions are: Deflection of the beam y(x), bending moments M(x) and shear forces V(x),
and the bending stresses (normal stress σn(x) – normal to the cross- section of the beam, and the shear
stress σs(x) - acting on the cross-section of the beam)
Quick Review of Simple Beam Theory – Cont’d
A B x B
A x
Pb a b Pa x
Ra = Rb = Ra =
wL wL
L L Rb =
L 2 L 2
There are two types of stresses induced to the beam subjected to external forces:
(1) The normal stress along the x-coordinate. It exists in perpendicular to the
cross-section of the beam (σxx), and
(2) The shear stress (σxz) that exists on the surfaces of the beam cross-section.
It also exists on the face along the x-coordinate (σyx) with the same magnitude
Expression for normal stress (σxx): σ xx (x ) = − M (x ) y (4.38)
I
where M(x) = bending moment at location x, y = distance from center of the beam
cross-section to the depth of the cross-section in y direction, and I = section
moment of inertia y
3
b
For beams with rectangular cross-section: bh
I= z H
12
πd 4 y
For beams with circular cross-section: I=
64
These stresses are induced by the shear force V(x) at the various cross-sections along the beam
y
dA
c1
y y1 H
z
V (x ) c1 V (x )
σ yx =
Ib ∫
y1
ydA =
Ib
Q (4.39)
V (x ) h 2
We have σ yx ( x ) = − y12 with (σ ) yx max =
3V
at the center line of the beam where y1=0
2I 4 bh
Quick Review of Simple Beam Theory – Cont’d
This theory relates the induced deflection (deformation of beams) and the applied forces
Induced
y Original straight beam: deflection by y P – concentrated force
applied forces w(x) – distributed load Applied
Forces
per unit length
x x
x Deflected state Induced
Deflection and Bending Moment M(x) Induced deflection: v(x) Deflection v(x)
and Shear force V(x) relations: & M(x), and V(x)
d 2 v( x ) M ( x ) d 2 v( x ) d 3v ( x )
2
= M ( x ) = EI and V ( x ) = EI (4.40a,b)
dx EI dx 2
dx 3
Induced Deflection v(x) may e obtained by solving the 4th order differential equation:
d 4 v( x )
EI 4
=0
dx
where E = Young’s modulus of the beam material
I = Section moment of inertia
FE Formulation of Beam Elements
Principal reference: “A First Course in the FEM” 6th Edition, Cengage Learning by Daryl Logan, 2017
ρ = Radius of curvature
y Uniform distribution load: w N/m deflection curve at x
y Θ = slope of deflection
Beam curve at x
Element V(x) = Deflection at x
ρ
A B x Node i
A B Node j
v(x)
x
x
L x
v(x) Datum
L
line
x
Primary Quantities L
vi x=0 x=L
In the element: The deflection v(x) θ
i
At the nodes: The deflections v and slope θ, or as expressed as: {d } = at Node i and Node j
v j
θ j
FE Formulation of Beam Elements
Principal reference: “A First Course in the FEM” 6th Edition, Cengage Learning by Daryl Logan, 2017
Beam elements
Contrary to the “truss elements,” beam elements deforms from its original straight shape into curved bent shape due to lateral
(or transverse) forces or applied couples (moments).
Deformed shape
Mi Node j Induced
Actions Reactions
Node i X
Original shape θ(x) in the element at Node i at Node j
vi vj Mj
x Lateral forces Vi Linear Linear Linear
and Vj displacement displacement displacement
Vi L Vj v(x) vi vj
v( x ) = a1 x 3 + a2 x 2 + a3 x + a4 (4.41)
dv(x )
At Node i: vi = v( x ) x =0 = a4 θi = = a3
dx x =0
x=0 x=L
v j = v( x ) x = L = a1 L3 + a2 L2 + a3 L + a4
2
At Node j:
dv( x )
θj = = 3a1 L2 + 2a2 L + a3
dx x = L
FE Formulation of Beam Elements – Cont’d
Principal reference: “A First Course in the FEM” 6th Edition, Cengage Learning by Daryl Logan, 2017
Derive Interpolation Function
We may thus express the element deflection v(x) in terms of the nodal deflections by
substituting the 4 constant coefficients into Equation (3.38) and after re-arranged terms
to yield:
v( x ) = 3 (vi − v j ) + 2 (θ i + θ j ) x 3
2 1
L L
+ − 2 (vi − v j ) − (2θ i + θ j ) x 2 + θ i x + vi
3 1
L L vi
θ
The above expression can also be expressed as: v( x ) = N iv { N iθ N jv N jθ }
i
(4.42)
v j
θ j
where N iv =
1
2 x(3
− 3 x 2
L + L3
) N iθ =
1 3
L3
x(L − 2 x 2 2
L + xL3
)
L3
(4.43)
1
(
N jv = 3 − 2 x 3 + 3 x 2 L ) N jθ
1
(
= 3 x 3 L − x 2 L2
L
)
L
We thus have: v(x) = [N(x)]{d} with Interpolation function: [N ( x )] = {N iv N jθ }
N iθ N jv (4.44)
We realize the fact that the rotation (or the slope) of the deflected beam is: θ = dv(x ) for 0 ≤ x ≤ L
dx
An expanded view of a deformed beam in x-direction:
Contracted top edge
-u
Undeformed neutral axis
Rotation of the beam section at x
X dx
dv( x)
We may find from the above diagram of expanded beam that the stretch u can be obtained by: u ( x) = − yθ = − y
dx
d dv( x ) d 2 v( x )
or ε xx ( x, y ) = − y
du
But from theory of elasticity: ε xx = = −y (4.45)
dx dx dx dx 2
FE Formulation of Beam Elements – Cont’d
Principal reference: “A First Course in the FEM” 6th Edition, Cengage Learning by Daryl Logan, 2017
M (x ) y
We will have the normal stress σxx= Eεxx = − (4.38)
I
V (x )
and the shear stress: σ yx (x ) = Q (4.39) with Q to be the shear moment
Ib
The relationships between the bending moment M(x) and Shear force V(x) are expressed in Equation (4.40a) and (4.40b) respectively:
d 2 v( x ) d 3v ( x )
M ( x ) = EI and V ( x ) = EI
dx 2 dx 3
FE Formulation of Beam Elements – Cont’d
Principal reference: “A First Course in the FEM” 6th Edition, Cengage Learning by Daryl Logan, 2017
We may express the applied actions to the beam element in applied forces fiy, fjy and moments mi, mj in terms of the
element deflections V(x) represented by the assumed polynomial function in Equation (4. 41) to be:
The beam element: Deformed shape
mi Node j
Node i X
Original shape θ(x)
vi
x v j mj
fiy L -fjy
X=0 X=L
Nodal forces:
d 3v ( x )
f iy = V = EI 3
=
EI
3
(12vi + 6Lθi − 12v j + 6Lθ j )
dx x =0 L
d 3v ( x )
= 3 (− 12vi − 6 Lθ i + 12v j − 6 Lθ j )
EI
f jy = −V = − EI
d 2 v(x )
3
dx x = L L
Nodal moments: mi = − M = − EI =
EI
dx 2 x =0 L3
6 Lvi(+ 4 L2
θ i − 6 Lv j + 2 L2
θj )
d 2 v(x )
m j = M (x ) = EI 2
=
dx x = L L
EI
3
(
6 Lvi + 2 L2θ i − 6 Lv j + 4 L2θ j )
FE Formulation of Beam Elements – Cont’d
Principal reference: “A First Course in the FEM” 6th Edition, Cengage Learning by Daryl Logan, 2017
The above expressions may e express in the following matrix form for the element stiffness equations:
vi 12 6 L − 12 6 L f iy
θ 6 L 4 L2 − 6 L 2 L2 m
i EI i (4.46)
=
v
j L3
− 12 − 6 L 12 − 6 L f jy
θ j 2
− m j
2
6 L 2 L 6 L 4 L
Nodal unknown Applied nodal transverse forces
quantities(transverse = Stiffness matrix [Ke] x
& moments
displacements & rotations
We thus have the stiffness equation of a beam element to be:
12 6 L − 12 6 L
2
EI 6 L 4 L − 6 L 2 L
2
[ke ] = 3 (4.47)
L − 12 − 6 L 12 − 6 L
2
6 L 2 L2
− 6 L 4 L
FE Formulation of Beam Elements – Cont’d
Principal reference: “A First Course in the FEM” 6th Edition, Cengage Learning by Daryl Logan, 2017
Example 4 Determine the deflection of a cantilever beam at the half span and at the point under the load. Dimensions and
applied lateral force are shown in the figure. The beam is made of a material with Young’s modulus E = 10000 MPa
ℓ/2=0.5 m
Node 2 P = 20 N
Node 1
• x
• • Node 3 Solution:
We realize the fact that we are seeking solutions at the mid-span and at the point
ℓ=1m under the applied force. It id reasonable to discretize the beam into to two (2) elements,
as shown below:
Beam Cross-section:
Element 1: Element 2:
b = 1 cm
Node 1 Node 2 Node 2 Node 3
h = 2 cm
• L = 0.5 m • • L = 0.5 m •
v1,, θ1 v2,θ2 v2,, θ2 v3,θ3
= 0.667 ×10 −8 m 4 matrices shown in Equation (4.47) for both Element 1 and 2.
12 12
EI = (1010)(0.667x10-8) = 66.7
FE Formulation of Beam Elements – Cont’d
Example 4
12 3 − 12 3 v1 f1
3 − θ m
533.6
1 3 0. 5 1 = 1 (b)
− 12 − 3 12 − 3 v2 f 2
3 0. 5 − 3 1 θ 2 m2
in which v1,θ1, v2, θ2 are the respective deflections and rotations in Node 1 and Node 2 respectively, whereas
f1, m1 and f2, m2 are the applied lateral forces and moments at Node 1 and Node 2 respectively
FE Formulation of Beam Elements – Cont’d
Example 4
Due to the fact that Element 2 has the same length and is made of the same material as Element 1, with the only difference
of the associated nodes, we can express the same element equation for Element 1 for element 2 as shown below:
12 3 − 12 3 v2 f 2
3 − θ m
533.6
1 3 0. 5 2 = 2
− 12 − 3 12 − 3 v3 f 3 (c)
3 0 .5 − 3 1 θ 3 m3
Assembly of element equations in (b) and (c) for Overall stiffness equations:
Because Node 2 happens to be common node shared by Element 1 and 2, we need to
assemble the element equation following the “MAP” stipulated in Step 5 in Chapter 3:
[K ] = [K e1 ]+ [K e2 ] = Sum
Node 2
Applied loads at Node 2 in Element 1 and 2 should be summed up too in the load matrix
FE Formulation of Beam Elements – Cont’d
12 3 − 12 3 0 0 v1 f1
3 1 −3 0.5 0 0 θ1 m1
− 12 − 3 (12 + 12) (− 3 + 3) − 12 3 v2 f 2
533.6 =
3 0.5 (− 3 + 3) (1 + 1) − 3 0.5 θ 2 m2
0 0 − 12 −3 12 − 3 v3 f 3
0 0 3 0.5 − 3 1 θ 3 m3
The assembled overall stiffness equation for the beam structure is:
12 3 − 12 3 0 0 v1 f1
3 1 − 3 0. 5 0 0 θ m
1 1
− 12 − 3 24 0 − 12 3 v2 f 2
533.6 = (d)
3 0.5 0 2 − 3 0.5 θ 2 m2
0 0 − 12 − 3 12 − 3 v3 f 3
0 0 3 0.5 − 3 1 θ 3 m3
We may solve the six unknown responses at the three nodes in the beam structure from Equation (d)
FE Formulation of Beam Elements – Cont’d
We realize that the following boundary and applied loading conditions apply:
Example 4
v1 and θ1 = 0 for having Node 1 being fixed at the built-in end, and
f1, m1, f2,m2 and m3 = 0. the only non-zero load is the applied force f3 = -20 N at Node 3
We thus have the overall stiffness equation of the beam structure with s[ecified boundary
and loading conditions take the form:
12 3 − 12 3 0 0 0 0
3 1 − 3 0. 5 0 0 0 0
− 12 − 3 24 0 − 12 3 v2 0
533.6 = (e)
3 0. 5 0 2 − 3 0.5 θ
2 0
0 0 − 12 − 3 12 − 3 v3 − 20
0 0 3 0.5 − 3 1 θ 3 0
The four unknown at Node 2 and 3 can be solved from Equation (e) by partitioning the above matrix equations following
Step 6 in Chapter 3.
FE Formulation of Beam Elements – Cont’d
Example 4 3 − 12 3 0 0 0
12 0
3 1 − 3 0.5 0 0 0 0
− 12 − 3 24 0 − 12 3 v2 0
533.6 =
3 0. 5 0 2 − 3 0.5 θ
2 0
0 0 − 12 − 3 12 − 3 v3 − 20
0 0 3 0.5 − 3 1 θ 3 0
K aa K ab qa Ra
In the form of: K =
ba K bb qb Rb
where {qa} = specified (known) nodal quantities; {Rb} = specified (known) applied resulting actions, from which
we may obtain:
{qb} = [Kbb]-1 ({Rb} – [Kba]{qa}) (6.12, REF)
In the present case, we have {qa} = {0 0}T We thus have the 4 unknowns obtained by: {qb} = [Kbb]-1 {Rb}
FE Formulation of Beam Elements – Cont’d
Example 4 3 − 12 3 0 0 0
12 0
3 1 − 3 0.5 0 0 0 0
− 12 − 3 24 0 − 12 3 v2 0
533.6 =
3 0. 5 0 2 − 3 0.5 θ
2 0
0 0 − 12 − 3 12 − 3 v3 − 20
0 0 3 0.5 − 3 1 θ 3 0
v2 24 0 − 12 3 0 0.3333 1 0.8333 1
θ 0 − 3 3 0 1 4 3 4
{qb } = 2 {Rb } = [K bb ]−1 =
2
[K bb ] =
One will find that:
− 12 − 3 12 − 3
v3 − 20 0.8333 3 2.6667 4
θ 3 3 0.5 − 3 1 0
1 4 4 8
FE Formulation of Beam Elements – Cont’d
Example 4 We may thus solve for the 4 primary unknown quantities in {qb} from the following equations:
Conventional wisdom shows that serious “stress concentration” can occur in the vicinity of a structure with drastic geometry changes.
The parts near the holes and the tapered areas and the notches in both plates are vulnerable for stress concentrations.
We have mentioned that the problem of plate with a small hole may be solved by the theory of “advanced strength of materials.”
However, plates with tapered and notched edges as shown above can not be handled by any existing method from classical
theories of elasticity, and the FEM appears the only available method of the solution.
Examples of FE Stress Analysis of Solid Structures
ο x
The strain – displacement relation in Equation (4.3) is modified to the new form:
∂u (x, y ) ∂
0
∂x ∂x
∂v( x, y ) ∂ u ( x, y )
{ε (x, y )} = =0
(4.48)
∂y ∂y v( x, y )
( ) ( )
∂u x, y + ∂v x, y ∂ ∂
∂y
∂x ∂y ∂x
1 ν 0 ε xx
{σ ( x, y )} = E 2 ν 1 0 ε yy (4.49)
1 −ν 1 −ν
0 0 ε xy
2
Formulation FE analysis for Two-Dimensional Stress Analysis of Solids with
Plate Elements
Let us now formulating the FE for plate structures such as with the discretization in a tapered plate as illustrated
In the figure below:
We begin our FE formulation of the plate structure with the expression of the “element displacement” components
{U(x,y)} in terms of the corresponding “nodal displacements” using an interpolation function {N(x,y)} as follows:
u (x, y )
{U (x, y )} = = {N (x, y )}{u} (4.)
v ( x , y ) ui
Element displacements
u
Interpolation Nodal displacements j
function u m
where the nodal displacement components, {u } = (4.)
vi
v j
vm
Formulation FE analysis for Two-Dimensional Stress Analysis of Solids with Plate Elements – Cont’d
We assume the elements used in this FE analysis are the “simplex” elements, meaning that the
Element displacements follow linear polynomial functions in relating their nodal displacements:
u ( x, y ) = α1 + α 2 x + α 3 y along the x-coordinate
u ( x, y )
We will thus have: {U (x, y )} =
v( x, y )
1 x y 0 0 0
= {α
1 α 2 α 3 α 4 α 5 α 6 }T (4.28)
0 0 0 1 x y
[R(x, y )] =
1 x y 0 0 0
The matrix [R(x,y)] in Equation (4.28) has the form: (4.29)
0 0 0 1 x y
Formulation FE analysis for Two-Dimensional Stress Analysis of Solids with Plate Elements – Cont’d
Derivation of interpolation function {N(x,y)} - cont’d
We are able to expand Equation (4.28) into the following form for the nodal displacements:
ui 1 xi yi 0 0 0 α1
u 1 x yj 0 0 0 α 2
j j
um 1 xm ym 0 0 0 α 3
{u} = = (4.30)
vi 0 0 0 1 xi yi α 4
v j 0 0 0 1 xj y j α 5
vm 0 0 0 1 xm ym α 6
Displacement components
of 3 nodes in the element Constant coefficients
Specified coordinates of
the 3 nodes in the element
Formulation FE analysis for Two- Dimensional Stress Analysis of Solids with Plate Elements – Cont’d
Derivation of interpolation function {N(x,y)} - cont’d
with [A]-1 = the inverse of the nodal coordinate matrix [A] in Equation (4.32)
Formulation FE analysis for Two- Dimensional Stress Analysis of Solids with Plate Elements – Cont’d
Derivation of interpolation function {N(x,y)} - cont’d
[R(x, y )] =
1 x y 0 0 0 (4.29)
Also with Equation (4.29) with:
0 0 0 1 x y
We will obtain the following expression after substituting the matric {α} in Equation (4.33) into Equation (4.28), yielding:
By comparing Equation (4.35) with Equation (4.26), we have the interpolation function of this simplex element to be:
[N(x,y)] = [R(x,y)][h] (4.36)
with Matrices [R(x,y)] in Equation (4.29) and [h] in Equation (4.34)
Formulation FE analysis for Two- Dimensional Stress Analysis of Solids with Plate Elements – Cont’d
The interpolation function {N(x,y)} - cont’d
We may refer to the interpolation function {N(x,y)} that we derived for the 3-node
ui plate elements in ui
Chapter 3 for a comparable case with the current analysis with the form of: v
vi i
u (x, y ) u N i 0 u j
{U (x, y )} = = {N (x, y )}{u} = {N i (x, y ) N j (x, y ) N m (x, y )} j =
0 Nj 0 Nm
v ( x , y ) v j 0 Ni 0 Nj 0 N m v j (4. a)
um um
vm vm
ui
u
An alternative j
u ( x, y ) Ni 0 um
{U (x, y )} =
Nj Nm 0 0
expression: = { ( )}{ } = 0 (4. b)
N m vi
N x , y u
v ( x , y ) 0 0 Ni Nj
v j
where N i ( x, y ) =
1
A
[
(x j ym − xm y j )+ (y j − ym )x + (xm − x j )y ]
vm
N j ( x, y ) =
1
[(xm yi − xi ym )+ ( ym − yi )x + (xi − xm )y ]
A
N m ( x, y ) =
1
A
[
j (xi y j − x j yi )+ ( yi − y j )x + (x j − xi )y ] (3.2)
Because the element equation is derived by minimizing the Potential energy in gthe deformed solid, we need
to derive the expression of “strain energy” in terms of nodal displacements (the primary quantities in the analysis).
We will first express the element strain vs. nodal displacements by Equation (4.12) as:
By using Equation (4.4) for [D] and Equation (4.36b) for [N(x,y)}, we get the matrix [B(x,y)] in the following:
( y j − y m ) ( ym − yi ) (yi − y j ) 0 0 0
[B] = 1 0 0 0 ( xm − x j ) ( xi − xm ) (x j − xi ) (4.)
(xm − x j ) (xi − xm ) (x j − xi ) (y j − ym ) ( ym − yi ) (yi − y j )
2A
The potential energy in the deformed element , according to Equation (4.18) is:
P({u}) =
1
{ } [ ] [C ][B]{u}dv
2 ∫v
T T
u B
(4.)
− ∫ {u} [N ( x, y )] { f }dv − ∫ {u} [N (x, y )] {t}ds
T T T T
v s
Formulation FE analysis for Two- Dimensional Stress Analysis of Solids with Plate Elements – Cont’d
Derivation of Element equation
The element equation is obtain by minimizing the potential energy in Equation (4.38) as:
∂P({u})
=0
∂{u}
Leading to the following element equation:
[K e ]{u} = {p} (4.)
{p} = [ ( )] { } [ ( )] {t}ds
T
Nodal forcwe matrix = ∫v + ∫s (4.)
T
N x , y f dv N x , y
The integration in Equation (4.40) with respect to the volume of the element may turn out to be tedious. However, if the size of
the element is not too large, this integration may be approximated by the following expression without significant error:
in which w is the thickness of the plane element, and A is the plane area that can be computed by Equation (3.2a).
A = (xi y j − x j yi ) + (x j ym − xm y j ) + ( xm yi − xi ym ) = the area of the element made of triangle (ijm)
Numerical example on FE analysis for Two-Dimensional Stress Analysis of Solids with a Plate Element
A structure made of a triangular plate defined by three corners at A,B and C. A force P is applied at corner C as shownin
the figure.
Find the following:
Solution:
And then using Equation (4.36) to obtain the interpolation function in the form of Equation (4.36b) as:
We are now ready to derive the element matrix for the structure with the newly derived interpolation
function. We will obtain first the [B] matrix from using Equation (4.37) as:
− 4 4 0 0 0 0
[B ] = 1 0 0 0 − 3 − 3 6
24
− 3 − 3 6 − 4 4 0
and the [C] matrix FROM Equation (4.7) and the coefficient matrix of Equation (4.25):
1 0.3 0
[C ] = 10.99 ×106 0.3 1 0
0 0 0.35
We are now ready to determine the element stiffness matrix in Equation (4.39) from Equations (4.40) and (4.42) as:
Because there is only one element in the structure, the above element equation is also the “overall stiffness
equation” of the structure, from which we may solve the “displacement components of ALL nodes after
Making necessary interchange of rows and columns in the above equation (according to the rule stipulatedin
Step 5 in Chapter 3.
Numerical example on FE analysis for Two-Dimensional Stress Analysis of Solids with a Plate Element – Cont’d
Thus, after the necessary interchanges of rows and columns, we reached the partition of the element equation
In the form of Equation (6.12, Ref) as:
The above simultaneous equations may be solved by either matrix inversion method or Gaussian elimination method, with:
uj = 0.16x10-3 inch, um = 0.38x10-3 inch and vm = -0.093x10-3 inch
Numerical example on FE analysis for Two-Dimensional Stress Analysis of Solids with a Plate Element – Cont’d
(b) The displacements in the element:
We may use the interpolation function {N(x,y)} to determine the displacement components everywhere in the element:
For the displacement in the x-direction:
●P u(x,y) = (1 –- 0.167x – 0.125y)ui + (0.167x -0.125y)uj + 0.25 um = 0 +(0.167x-0.125y)x0.16x10-3+0.25x0.38x10-3y
V(x,y) = (1 –- 0.167x – 0.125y)vi + (0.167x -0.125y)vj + 0.25 vm = -0.09264x0.25x10-3y
For instance the displacements at Point P(3,2) have the values of:
u(3,2) = 0 +(0.167x3-0.125x2)x0.16x10-3+0.25x0.38x10-3 x2 = 0.23x10-3 inch
v(3,2) = -0.09264x0.25x10-3x2 = -.04632x10-3 inch
(c) The strain components in the element by using Equation (4.12) with {ε} = [B]{u}:
ui = 0
u = 0.16
ε xx (x, y ) − 4 4 0 0 0 0 j 26.65
u = 0.38
{ε (x, y )} = ε yy (x, y ) = 1 0 0 0 − 3 − 3 6 m −3
×10 = − 23.16 ×10
−6
ε (x, y ) 24 − 3 − 3 6 − 4 4 0 vi = 0 75
xy v = 0
j
vm = −0.09
Numerical example on FE analysis for Two-Dimensional Stress Analysis of Solids with a Plate Element – Cont’d
The stresses in the element may be obtained by using the generalized Hooke’s Law in Equation (4.25):
1 ν 0 ε xx 1 0.3 0 26.65 216.5
{σ ( x, y)} = E 2 ν 1 0 ε yy = 10.99 ×106 0.3 1 0 − 23.16 ×10−6 = − 166.67 psi
1 −ν 1 −ν
0 0 ε xy 0 0 0.35 75 288.66
2
Rix − 4 0 − 3 − 866
R 4 0
jx
0 − 3 216.5
Rmx 1 0 0 6 866
{R} = = − 166.67 ×12 ×1 =
R
iy 24 0 − 3 − 4 288.66 − 327 . 3
R jy 0 −3 4 827.3
Rmy 0 6 0 − 500
The reactions at Node i therefore have numerical values at: Rix = 866 lb towards left, and Riy=327.3 lb in the downward direction
Numerical example on FE analysis for Two-Dimensional Stress Analysis of Solids with a Plate Element – Cont’d
This is good lesson for any intelligent FE user to learn and exercise
SUMMARY
1. This chapter presents FE formulation for solids deformed by applied forces, either in the form of body force
or as surface tractions.
2. The deformation of the solids is supposed to be small so that the induced stresses are within the elastic
limit of the material
4. General FE formulation for “simplex element” (with interpolation functions follow linear polynomial functions)
were derived.
5. FE formulations for special cases for: (a) one-dimensional bar elements, (b) 2-dimensional plate elements,
and 3-dimensional torus elements are derived.
6. FE formulations of stress analysis of deformable solids uses minimization of potential energy in the deformed
state to derive the element equations.
7. Because simplex elements are used for all these cases for simplicity in deriving element equations. However,
these simplex element result in “constant” stresses and strains in each individual elements. Intelligent users
will thus need to assign many more smaller elements in the regions where hgi gradients of primary
quantities are expected to achieve better result.
FE Formulation of Truss Elements Using 1-D Bar Elements – Cont’d
We assume the element displacement components follow linear polynomial functions as shown below:
where α1, α2, α3, α4, α5, α6 in the matrix {α}T are constants to be determined with specified nodal coordinates later.
[R(x, y )] =
The matrix [R(x,y)] in Equation (4.26) has the form: 1 x y 0 0 0 (4.27)
0 0 0 1 x y