MKC525E Midterm Exam #1
1. Obtaining Strong Form Of The Problem
k=50 N/mm
Definition of stress resultant;
𝑀 = ∫ 𝑧. 𝜎𝑥𝑥 𝑑𝐴 𝑉 = ∫ 𝜎𝑥𝑧 𝑑𝐴 (1)
𝐴 𝐴
Equilibrium equation;
𝑑𝑉 𝑑𝑀
+ 𝑞(𝑥) − 𝑘𝑢 = 0 , =𝑉
𝑑𝑥 𝑑𝑥
(2)
𝑑2𝑀
+ 𝑞(𝑥) − 𝑘𝑢 = 0
𝑑𝑥 2
Kinematic and constitutive relations;
𝑑2𝑢 (3)
𝜀𝑥𝑥 = −𝑧 2
𝑑𝑥
𝑑2𝑢
𝜎𝑥𝑥 = 𝐸. 𝜀𝑥𝑥 = −𝐸𝑧 (4)
𝑑𝑥 2
Substitute equation (4) in Equation (1) and get moment deflection equation;
𝑑2𝑢 𝑑 2 𝑢(5)
𝑀 = ∫ 𝑧. 𝜎𝑥𝑥 𝑑𝐴 = ∫−𝑧 2 𝐸 𝑑𝐴 = −𝐸𝐼(𝑥)
𝐴 𝑑𝑥 2 𝑑𝑥 2
𝐴
In order to obtain strong form put equation (5) in equilibrium equation (2) and we will get governing
equilibrium equation;
STRONG FORM
𝑑2 𝑑2𝑢 (6)
(𝐸𝐼(𝑥) ) + 𝑘𝑢 − 𝑞(𝑥) = 0
𝑑𝑥 2 𝑑𝑥 2
2. Obtaining Weak Form Of The Problem
θ1 θ2 Q2 Q4
xe xe+1 xe xe+1
u1 u2 Q1 Q3
Generalized displacements Generalized forces
Consider notation given above and, multiply the governing differential equation with a weight function w
and integrate over a typical element given above;
𝑥𝑒+1 𝑥𝑒+1 𝑥𝑒+1
𝑑2 𝑑2𝑢 (7)
∫ 𝑤. (𝐸𝐼(𝑥) ) 𝑑𝑥 + ∫ 𝑤. 𝑘𝑢. 𝑑𝑥 − ∫ 𝑤. 𝑞(𝑥)𝑑𝑥 = 0
𝑥𝑒 𝑑𝑥 2 𝑑𝑥 2 𝑥𝑒 𝑥𝑒
With help of product rule ;
′
𝑑 𝑑2𝑢 𝑑𝑤 𝑑 𝑑2𝑢 𝑑2 𝑑2𝑢
[𝑤. (𝐸𝐼(𝑥) 2 )] = (𝐸𝐼(𝑥) 2 ) + 𝑤 2 (𝐸𝐼(𝑥) 2 ) Integrate this expression and ;
𝑑𝑥 𝑑𝑥 𝑑𝑥 𝑑𝑥 𝑑𝑥 𝑑𝑥 𝑑𝑥
𝑥
𝑥𝑒+1
𝑑2 𝑑2𝑢 𝑑 𝑑 2 𝑢 𝑒+1 𝑥𝑒+1
𝑑𝑤 𝑑 𝑑2 𝑢
∫ 𝑤. 2 (𝐸𝐼(𝑥) 2 ) 𝑑𝑥 = 𝑤. (𝐸𝐼(𝑥) 2 ) | − ∫ (𝐸𝐼(𝑥) 2 ) 𝑑𝑥 (8)
𝑥𝑒 𝑑𝑥 𝑑𝑥 𝑑𝑥 𝑑𝑥 𝑥𝑒 𝑥𝑒 𝑑𝑥 𝑑𝑥 𝑑𝑥
Substitute equation (8) in equation (7);
𝑥
𝑑 𝑑2 𝑢 𝑒+1 𝑥𝑒+1
𝑑𝑤 𝑑 𝑑2𝑢 𝑥𝑒+1 𝑥𝑒+1
(9)= 0
𝑤. (𝐸𝐼(𝑥) 2 ) | − ∫ (𝐸𝐼(𝑥) 2 ) 𝑑𝑥 + ∫ 𝑤. 𝑘𝑢. 𝑑𝑥 − ∫ 𝑤. 𝑞(𝑥)𝑑𝑥
𝑑𝑥 𝑑𝑥 𝑥𝑒 𝑥𝑒 𝑑𝑥 𝑑𝑥 𝑑𝑥 𝑥𝑒 𝑥𝑒
Again with help of product rule;
′
𝑑𝑤 𝑑2𝑢 𝑑2 𝑤 𝑑2 𝑢 𝑑𝑤 𝑑 𝑑2𝑢 Integrate this expression and ;
[ (𝐸𝐼(𝑥) 2 )] = 𝐸𝐼(𝑥) + (𝐸𝐼(𝑥) )
𝑑𝑥 𝑑𝑥 𝑑𝑥 2 𝑑𝑥 2 𝑑𝑥 𝑑𝑥 𝑑𝑥 2
𝑥
𝑥𝑒+1
𝑑𝑤 𝑑 𝑑2𝑢 𝑑𝑤 𝑑 2 𝑢 𝑒+1 𝑥𝑒+1
𝑑2𝑤 𝑑2𝑢
∫ (𝐸𝐼(𝑥) 2 ) 𝑑𝑥 = (𝐸𝐼(𝑥) 2 ) | − ∫ 𝐸𝐼(𝑥) 2 2 𝑑𝑥 (10)
𝑥𝑒 𝑑𝑥 𝑑𝑥 𝑑𝑥 𝑑𝑥 𝑑𝑥 𝑥𝑒 𝑥𝑒 𝑑𝑥 𝑑𝑥
Substitute equation (10) in equation (9);
𝑥 𝑥
𝑥𝑒+1
𝑑2 𝑤 𝑑2 𝑢 𝑥𝑒+1 𝑙
𝑑𝑤 𝑑2 𝑢 𝑒+1 𝑑 𝑑2 𝑢 𝑒+1
∫ (𝐸𝐼(𝑥) 2 2 + 𝑤. 𝑘𝑢) 𝑑𝑥 − ∫ 𝑤. 𝑞(𝑥)𝑑𝑥 − (𝐸𝐼(𝑥) 2 ) | + 𝑤. (𝐸𝐼(𝑥) 2 ) | = 0 (11)
𝑥𝑒 𝑑𝑥 𝑑𝑥 𝑥𝑒 𝑑𝑥 𝑑𝑥 𝑥 𝑑𝑥 𝑑𝑥 𝑥 𝑒 𝑒
Where ;
𝑑 𝑑2 𝑢 𝑑 𝑑2𝑢
𝑄1 𝑒 = [ (𝐸𝐼(𝑥) 2 )] 𝑄3 𝑒 = − [ (𝐸𝐼(𝑥) 2 )]
𝑑𝑥 𝑑𝑥 𝑥 𝑑𝑥 𝑑𝑥 𝑥
𝑒 𝑒+1
𝑒 𝑑2 𝑢 𝑒 𝑑2𝑢
𝑄2 = [𝐸𝐼(𝑥) 2 ] 𝑄4 = − [𝐸𝐼(𝑥) 2 ]
𝑑𝑥 𝑥 𝑑𝑥 𝑥
𝑒 𝑒+1
Rearrange equation (11) and we will get corresponding weak form of the Bernoulli-Euler Beam on elastic
foundation problem;
WEAK FORM
𝑥𝑒+1 𝑥𝑒+1
𝑑2𝑤 𝑑2 𝑢 𝑑𝑤
∫ (𝐸𝐼(𝑥) 2 2
+ 𝑤. 𝑘𝑢) 𝑑𝑥 − ∫ 𝑤. 𝑞(𝑥)𝑑𝑥 − 𝑤(𝑥𝑒 )𝑄1 𝑒 − (− ( ) 𝑄2 𝑒 )
𝑥𝑒 𝑑𝑥 𝑑𝑥 𝑥𝑒 𝑑𝑥 𝑥𝑒 (12)
𝑑𝑤
−𝑤(𝑥𝑒+1 )𝑄3 𝑒 − (− ( ) 𝑄4 𝑒 ) = 0
𝑑𝑥 𝑥𝑒+1
Where Q2 and Q4 corresponds Moment (M) , Q1 and Q3 corresponds Shear forces (Q).
In terms of bilinear and linear operators ;
𝑎(𝑤, 𝑢) = (𝑤, 𝑙) − 𝑤,𝑥 (𝑥𝑒 )𝑀 + 𝑤(𝑥𝑒 )𝑄 − 𝑤,𝑥 (𝑥𝑒+1 )𝑀 + 𝑤(𝑥𝑒+1 )𝑄
𝑥𝑒+1
𝑑2𝑤 𝑑2𝑢
𝑎(𝑤, 𝑢) = ∫ (𝐸𝐼(𝑥) + 𝑤. 𝑘𝑢) 𝑑𝑥
𝑥𝑒 𝑑𝑥 2 𝑑𝑥 2
𝑥𝑒+1
(𝑤, 𝑙) = ∫ 𝑤. 𝑞(𝑥)𝑑𝑥
𝑥𝑒
Apply following boundary conditions on equation (12) then we will obtain corresponding weak form of
the bernoulli-euler beam on elastic foundation problem with 1 element;
𝑑𝑤
𝑤(𝑥𝑒+1 ) = ( ) = 0 (𝑓𝑖𝑥𝑒𝑑)
𝑑𝑥 𝑥𝑒+1
𝑥𝑒+1 𝑥𝑒+1
𝑑2𝑤 𝑑2 𝑢 𝑑𝑤
∫ (𝐸𝐼(𝑥) 2 2
+ 𝑤. 𝑘𝑢) 𝑑𝑥 − ∫ 𝑤. 𝑞(𝑥)𝑑𝑥 − 𝑤(𝑥𝑒 )𝑄1 𝑒 − (− ( ) 𝑄2 𝑒 ) = 0
𝑥𝑒 𝑑𝑥 𝑑𝑥 𝑥𝑒 𝑑𝑥 𝑥𝑒
𝑎(𝑤, 𝑢) = (𝑤, 𝑙) − 𝑤,𝑥 (𝑥𝑒 )𝑀 + 𝑤(𝑥𝑒 )𝑄
3. Finite Element Analysis
Interpolation Function
Cubic interpolation functions are derived by interpolation weight function and its derivative at the nodes.
Such polynomials are known as the Hermite family of interpolation functions. (Reddy,Introduction to Finite
Element)
𝑥 − 𝑥𝑒 2 𝑥 − 𝑥𝑒 3 𝑥 − 𝑥𝑒 2
𝜑1 = 1 − 3 ( ) + 2( ) 𝜑2 = −(𝑥 − 𝑥𝑒 ) (1 − )
ℎ𝑒 ℎ𝑒 ℎ𝑒 (13)
𝑥 − 𝑥𝑒 2 𝑥 − 𝑥𝑒 3 𝑥 − 𝑥𝑒 2 𝑥 − 𝑥𝑒
𝜑3 = 3 ( ) − 2( ) 𝜑4 = −(𝑥 − 𝑥𝑒 ) [( ) − ]
ℎ𝑒 ℎ𝑒 ℎ𝑒 ℎ𝑒
Finite Element Model
The finite element model of the Euler-Bernoulli beam is obtained by substituting the finite element
interpolation functions (𝜑𝑖 ) for w and u. Stiffness matrix will be obtained by using equations given in (13)
into a(w,u) operator. Remember weak form;
𝑥𝑒+1
𝑑2𝑤 𝑑2𝑢 𝑥𝑒+1
𝑑𝑤 𝑑𝑤
∫ (𝐸𝐼(𝑥) + 𝑤. 𝑘𝑢) 𝑑𝑥 − ∫ 𝑤. 𝑞(𝑥)𝑑𝑥 − 𝑤(𝑥𝑒 )𝑄1 𝑒 − (− ( ) 𝑄2 𝑒 ) − 𝑤(𝑥𝑒+1 )𝑄3 𝑒 − (− ( ) 𝑄 𝑒)
𝑥𝑒
2
𝑑𝑥 𝑑𝑥 2
𝑥𝑒 𝑑𝑥 𝑥𝑒 𝑑𝑥 𝑥𝑒+1 4
4 𝑥𝑒+1
𝑑2 𝜑𝑖 𝑑2 𝜑𝑗
𝐾𝑖𝑗 = ∑ ∫ (𝐸𝐼(𝑥) + 𝜑𝑖 . 𝑘𝜑𝑗 ) 𝑑𝑥 (14)
𝑥𝑒 𝑑𝑥2 𝑑𝑥2
𝑗=1
4 𝑥𝑒+1
𝑑 𝜑𝑖 𝑑𝜑𝑖
𝐹𝑖 = ∑ ∫ 𝜑𝑖 . 𝑞(𝑥)𝑑𝑥 + 𝜑𝑖 (𝑥𝑒 )𝑄1 𝑒 − ( ) 𝑄2 𝑒 + 𝜑𝑖 (𝑥𝑒+1 )𝑄3 𝑒 − ( ) 𝑄4 𝑒 Use equation 13 and rearrange Fi
𝑑𝑥 𝑥𝑒 𝑑𝑥 𝑥𝑒+1
𝑗=1 𝑥𝑒
4 𝑥𝑒+1
𝐹𝑖 = ∑ ∫ 𝜑𝑖 . 𝑞(𝑥)𝑑𝑥 + 𝑄𝑖 (15)
𝑥𝑒
𝑗=1
∑ 𝐾𝑖𝑗 𝑢𝑗 − 𝐹𝑖 = 0 (16)
𝑗=1
According to equation (14) the element stiffness matrix is obtained as ;
6 −3ℎ𝑒 −6 −3ℎ𝑒 156 −22ℎ𝑒 54 13ℎ𝑒
2𝐸𝐼 −3ℎ𝑒 2ℎ𝑒2 3ℎ𝑒 2
ℎ𝑒 𝑘ℎ𝑒 −22ℎ𝑒 4ℎ𝑒2 −13ℎ𝑒 −3ℎ𝑒2
[𝐾 𝑒 ] = 3 +
3ℎ𝑒 −6 3ℎ𝑒2 6 3ℎ𝑒 420 54 −13ℎ𝑒 156 22ℎ𝑒
[−3ℎ𝑒 ℎ𝑒2 3ℎ𝑒 2
2ℎ𝑒 ] [ 13ℎ𝑒 −3ℎ𝑒2 22ℎ𝑒 4ℎ𝑒2 ]
Comes from beam Comes from elastic foundation
According to equation (15) the force vector is obtained as ;
6 𝑄1
𝑞 ℎ
𝑒 𝑒 −ℎ𝑒 𝑄
{𝐹 𝑒 } = { } + { 2}
12 6 𝑄3
ℎ𝑒 𝑄4
Consistent load vector for cantilever beam
Question 3.1
Divide beam into 1 element as follows ;
θ1 θ2 Q2 Q4
xe xe+1 xe xe+1
u1 u2 Q1 Q3
Generalized displacements Generalized forces
Where he=2m
I=9.10-9 m4
E= 200.109 (N/m2)
q=1000 N/m
37149 −1.0483𝑒 + 007 12851 6.1841𝑒 + 006
−1.0483𝑒 + 007 3.8181𝑒 + 009 −6.1841𝑒 + 006 −2.8529𝑒 + 009
[𝐾 𝑔𝑙𝑜𝑏𝑎𝑙 ] = [ ] 𝑁/𝑚𝑚
12851 −6.1841𝑒 + 006 37149 1.0483𝑒 + 007
6.1841𝑒 + 006 −2.8529𝑒 + 009 1.0483𝑒 + 007 3.8181𝑒 + 009
𝑢1 𝑢1
𝜃1 𝜃
𝑢 = {𝑢 } = { 1 }
2 0
𝜃2 0
𝑄1 1000
𝑄 333.33
{𝐹 𝑒 } = { 2 } = { }
𝑄3 𝑄3
𝑄4 𝑄4
37149 −1.0483𝑒 + 007 12851 6.1841𝑒 + 006 𝑢1 1000
−1.0483𝑒 + 007 3.8181𝑒 + 009 −6.1841𝑒 + 006 −2.8529𝑒 + 009 𝜃1 333.33
[ ]{ } = { }
12851 −6.1841𝑒 + 006 37149 1.0483𝑒 + 007 0 𝑄3
6.1841𝑒 + 006 −2.8529𝑒 + 009 1.0483𝑒 + 007 3.8181𝑒 + 009 0 𝑄4
𝑢1 = 0.22864 𝑚𝑚 𝜃1 = 0.0007156 𝑟𝑎𝑑
Question 3.2
Divide beam into 2 element
θ1 θ2 θ2 θ3
1 2 2 3
u1 u2 u2 u3
18623 −2.6446𝑒 + 006 6377.4 1.522𝑒 + 006 0 0
−2.6446𝑒 + 006 4.9326𝑒 + 008 −1.522𝑒 + 006 −3.4861𝑒 + 008 0 0
[𝐾 𝑔𝑙𝑜𝑏𝑎𝑙 ]= 6377.4 −1.522𝑒 + 006 55772 0 6377.4 1.522𝑒 + 006
𝑁/𝑚𝑚
1.522𝑒 + 006 −3.4861𝑒 + 008 0 4.3113𝑒 + 009 −1.522𝑒 + 006 −3.4861𝑒 + 008
0 0 6377.4 −1.522𝑒 + 006 37149 2.6446𝑒 + 006
[ 0 0 1.522𝑒 + 006 −3.4861𝑒 + 008 2.6446𝑒 + 006 3.8181𝑒 + 009 ]
𝑢1 𝑢1
𝜃1 𝜃1
𝑢2 𝑢2
𝑢= 𝜃 =
2 𝜃2
𝑢3 0
{𝜃3 } { 0 }
𝑄1 1000
𝑄2 333.33
𝑄3 0
{𝐹 𝑒 } = =
𝑄4 0
𝑄5 𝑄 5
{𝑄6 } { 𝑄6 }
18623 −2.6446𝑒 + 006 6377.4 1.522𝑒 + 006 0 0 𝑢1 1000
−2.6446𝑒 + 006 4.9326𝑒 + 008 −1.522𝑒 + 006 −3.4861𝑒 + 008 0 0 𝜃1 333.33
6377.4 −1.522𝑒 + 006 55772 0 6377.4 1.522𝑒 + 006 𝑢2 0
=
1.522𝑒 + 006 −3.4861𝑒 + 008 0 4.3113𝑒 + 009 −1.522𝑒 + 006 −3.4861𝑒 + 008 𝜃2 0
0 0 6377.4 −1.522𝑒 + 006 37149 2.6446𝑒 + 006 0 𝑄5
[ 0 0 1.522𝑒 + 006 −3.4861𝑒 + 008 2.6446𝑒 + 006 3.8181𝑒 + 009 ] { 0 } { 𝑄6 }
𝑢1 = 0.7294 𝑚𝑚 𝜃1 = 0.005467 𝑟𝑎𝑑
𝑢2 = 0.0909 𝑚𝑚 𝜃2 = 0.000802 𝑟𝑎𝑑
Question 3.3
Divide beam into 10 element and calculate the nodal deflections and rotations by following Matlab Code.
Results;
u1 1.0096 u6 -0.00013175
Rot1 0.012759 Rot6 4.5482e-006
u2 -0.14208 u7 -0.00017162
Rot2 0.0007538 Rot7 -1.3775e-006
u3 -0.058032 u8 -7.9194e-008
Rot3 -0.00063094 Rot8 -2.9097e-007
u4 0.005031 u9 8.8756e-006
Rot4 -6.4418e-005 Rot9 5.949e-008
u5 0.0032096 u10 5.4387e-007
Rot5 3.0109e-005 Rot10 1.5496e-008