ES 128: Homework 2

Solutions
Problem 1
Show that the weak form of

d du
( AE ) + 2x = 0 on 1 < x < 3 ,
dx dx
 du 
σ (1) =  E  = 0.1 ,
 dx  x =1
u(3) = 0.001

is given by

3 dw du 3
∫
1 dx
AE dx = −0.1( wA ) x =1 + ∫ 2 xwdx
dx 1
∀w with w(3) = 0 .

Solution
We multiply the governing equation and the natural boundary condition over the
domain [1, 3] by an arbitrary weight function:

  d 
3 du  
∫ 1


w  AE
 dx  dx


+ 2 x dx = 0
 
∀w( x ) , (1.1)

  du 
 wA E − 0.1   = 0 ∀w(1) . (1.2)
  dx   x =1

We integrate (1.1) by parts as

x =3
3   d  du    du  3  dw du 
∫1
w  AE
  dx 
 dx =  wAE
dx   
 − ∫1 
dx  x =1  dx
AE
dx 
dx . (1.3)

Substituting (1.3) into (1.1) gives

3  dw du  3  du   du 
−∫  AE dx + ∫ 2wxdx +  wAE  −  wAE  =0 ∀w( x ) . (1.4)
1
 dx dx  1
 dx  x =3  dx  x =1

With w(3) = 0 and σ (1) = 0.1 , we obtain

 dw du 
∫1 2wxdx − 0.1(wA ) x=1
3 3
∫1  dx

AE
dx 
dx = ∀w(x ) with w(3) = 0 . (1.5)

1
Problem 2
Consider the (steel) bar in Figure 1. The bar has
a uniform thickness t=1cm, Young’s modulus
E=200 × 109 Pa, and weight density
ρ =7 × 10 kg / m . In addition to its self-weight,
3 3

the bar is subjected to a point load P=100N at

its midpoint.
(a) Model the bar with two finite elements.
(b) Write down expressions for the element
stiffness matrices and element body force
vectors.
(c) Assemble the structural stiffness matrix
K and global load vector F .
(d) Solve for the global displacement vector
d.
(e) Evaluate the stresses in each element.
(f) Determine the reaction force at the
support.

Solution
(a) Using two elements, each of 0.3m in length, we
obtain the finite element model in Figure 1a. In this 1
model, x 1(1 ) = 0 , x 2(1 ) = 0.3 , x 1( 2 ) = 0.3 , x 2( 2 ) = 0.6 ,
and A(x)=0.0012-0.001x .
El (1)
1
(b) For element 1, N (1) = [0.3 − x x ] ,
0.3
1
B (1) = [− 1 1] , the element stiffness matrix is 2
0.3 El (2)
0.3
K (1) = ∫ B (1)T AEB (1 ) dx
0

200 × 10 9 0.3  1 − 1
=
0.09 ∫
(0.0012 − 0.001 x )
0 dx
− 1 1  3
 0.7 − 0.7  x
= 10 9  ,
− 0.7 0.7  Figure 1a
the element body force vector is
f (1) = ∫
0
0.3
(
N (1)T ρAdx + N (1)T P )
x = 0.3

7 × 10 3
(0.3 − x )(0.0012 − 0.001 x )
0.3  0 
=
0.3 0 

∫
x (0.0012 − 0.001 x ) 

dx +  
100
 1.155   1 0 0
=  , and the scatter matrix is L(1) =  .
101.05 0 1 0

2
 1 1
For element 2, N (2) = [0.6 − x x − 0.3] , B (2) = [− 1 1] , the element
0.3 0.3
stiffness matrix is
0.6 200 × 10 9 0.6  1 − 1
K (2) = ∫ B (2)T AEB ( 2 ) dx = ∫0.3 (0. 0012 − 0.001 x ) − 1 1 dx
0.3 0.09  
 0.5 − 0.5
= 10 9  , the element body force vector
− 0.5 0.5 
0.6 7 × 10 3 0.6(0.6 − x )(0.0012 − 0.001 x )  0.84 
is f (2) = ∫ N (2)T ρAdx =
0.3 ∫0.3 ( x − 0.3)(0.0012 − 0.001 x )
 dx = 0.735 , and
0.3
 
0 1 0
the scatter matrix is L(2) =  .
0 0 1 
(c) The global stiffness matrix is
2
 0.7 − 0.7 0 
K = ∑L K L = L
eT e e (1)T
K L +L
(1) (1) (2)T (2)
K L (2) 9  
= 10 − 0.7 1.2 − 0.5 .
e =1
 0 − 0.5 0.5 
The global load vector is
2
 1.155 
f = ∑L f = LeT e (1)T (1)
f +L (2)T (2)
f

= 101.89 .

e =1
 0.735 
(d) Note that only the reaction force at node 1 is not zero, thus
r1 + 1.155 
f + r =  101.89  .
 0.735 
The resulting global system of equations is
 0.7 − 0.7 0   0  r1 + 1.155
 
10 − 0.7 1.2 − 0.5 u2  =  101.89  .
9

 0 − 0.5 0.5  u3   0.735 

Solving the above equation,
u2 = 1.46607 × 10 _ 7 (m), u3 = 1.48077 × 10 _ 7 (m), and r1 =-103.78(N).
(e) The stress field in element 1 is given by
1  0 
σ (1 ) ( x ) = EB (1) d (1) = 200 × 10 9 × [ − 1 1]  −7 
=9.7738 × 10 4 (Pa).
0.3 1.46607 × 10 
The stress field in element 2 is given by
1 1.46607 × 10 −7 
σ ( 2 ) ( x ) = EB (2) d (2) = 200 × 10 9 × [ − 1 1] −7 
=980 (Pa).
0.3 1.48077 × 10 
(f) The reaction force at the support Node 1 is -103.78N.

3
Problem 3
Consider the mesh shown in Figure 2. The model consists of two linear
displacement constant strain elements. The cross-sectional area is A=1, Young’s
modulus is E; both are constant. A body force b(x)=cx is applied.

(a) Solve and plot u(x) and ε ( x ) for

the FEM solution.
(b) Compare (by plotting) the finite
element solution against the exact
solution for the equation
d2u
Figure 2 E 2 = −b( x ) = −cx .
dx
(c) Solve the above problem using a single quadratic displacement element.
(d) Compare the accuracy of stress and displacement at the right end with that
of two linear displacement elements.
(e) Check whether the equilibrium equation and traction boundary condition
are satisfied for the two meshes.

Solution
(a) Using two linear displacement constant strain elements, each of l in length,
we obtain the finite element model with x 1(1 ) = 0 , x 2(1 ) = l , x 1( 2 ) = l , and x 2( 2 ) = 2l .
1 1
For element 1, N (1) = [l − x x ] , B (1) = [− 1 1] , the element stiffness matrix is
l l
EA  1 − 1  l cl 2 / 6
K (1) = = ∫0 =
(1) (1)T
, the element body force vector is f N cxdx  2 ,
l − 1 1  cl / 3
 1 0 0
and the scatter matrix is L(1) =  .
0 1 0
1 1
For element 2, N (2) = [2l − x x − l ] , B (2) = [− 1 1] , the element stiffness matrix
l l
EA  1 − 1 
is K (2) = , the element body force vector
l − 1 1 
2l 2cl 2 / 3 0 1 0
is f = ∫ N cxdx =  2  , and the scatter matrix is L(2) = 
(2) (2)T
.
l
5cl / 6 0 0 1 
The global stiffness matrix is
2
 1 −1 0 
EA 
K = ∑L K L = L K L + L K L =
eT e e (1)T (1) (1) (2)T (2) (2)
 − 1 2 − 1 .
e =1 l
 0 − 1 1 
The global load vector is

4
 2
 0.1667 
eT e (1)T (1) 
f = ∑ L f = L f + L f = cl  1.0  .
(2)T (2) 2 
e =1
0.8333
The resulting global system of equations is

 1 − 1 0   0  r1 + 0.1667 cl 2 
EA     
 − 1 2 − 1 u2  =  cl 2 .
l
 0 − 1 1  u3   0.8333cl  2

Solving the above equation,

cl 3 cl 3
u2 = 1.8333 , u3 = 2.6666 , and r1 =-2 cl 2 .
EA EA
When 0 ≤ x ≤ l
1  0  cl 2 x
u( x ) = N d = l − x x
(1) (1)
[ ] cl 3
 = 1.8333 ,
l 1.8333  EA
 EA 
2
du cl
and ε ( x ) = = 1.8333 .
dx EA
When l ≤ x ≤ 2l
 cl 3 
1  1.8333  3 2
u( x ) = N (2) d (2) = [2l − x x − l ] EA  = cl + 0.8333 cl x ,
3
l 2.6666 cl  EA EA
 EA 
du cl 2
and ε ( x ) = = 0.8333 .
dx EA

(b). The governing equation is

d 2u
EA 2 = −b( x ) = −cx ,
dx
where A=1. The boundary condition is u(0)=0, and σ (2l ) = 0 .
c
Solving this linear ODE, we obtain the exact solution u( x ) = − x 3 + c1 x + c2 .
6 EA
2cl 2
Since u(0)=0, c 2 = 0 . Since σ (2l ) = 0 , c 1 = . Thus
EA
c 2cl 2
u( x ) = − x +
3
x.
6 EA EA
cx 2 2cl 2
ε (x) = − +
2 EA EA
The comparisons between the approximation results and the exact solutions are
shown in Figures 2a (displacement), and 2b (strain).

5
 Exact Approximation
solution results
Exact
solution
u( x ) EA EA
ε (x)
l cl 2 Approximation cl 2
results

Figure 2a Figure 2b

(c) Using a single quadratic displacement element with x 1(1 ) = 0 , x 2(1 ) = l , and
x 3(1 ) = 2l , the element shape functions are

N (1)
=
(x − x )(x − x ) = ( x − l )( x − 2l ) = ( x − l )( x − 2l ) ,
(1)
2
(1 )
3
1
(x − x )(x − x ) (−l )(−2l )
(1 )
1
(1)
2
(1 )
1
(1)
3 2l 2

N (1)
=
(x − x )(x − x ) = x( x − 2l ) = x( x − 2l ) ,
(1)
1
(1 )
3
2
(x − x )(x − x ) l(−l )
(1 )
2
(1)
1
(1 )
2
(1)
3 −l 2

N (1) =
(x − x )(x − x ) = x( x − l ) = x( x − l ) .
(1)
1
(1 )
2
3
(x − x )(x − x ) 2l
(1 )
3
(1)
1
(1 )
3
(1)
2
2
2l 2

The corresponding B-matrix is

dN (1) x − 3l / 2
B (1)
1 = 1
= ,
dx l2
dN (1) 2 x − 2l
B (1)
2 = 2
= ,
dx − l2
dN (1)
3 x − l /2
B3 =
(1)
= .
dx l2
The element stiffness matrix is
 x − 3l / 2 
 l2 
2l 2l  2 x − 2l   x − 3l / 2 2 x − 2l x − l /2
K (1) = ∫ B(1)T EAB(1)dx = EA ∫   dx
0 2
 −l
2
 l
2
− l2 l 2 
 x − l /2 
 l 2 

6
   x − 3l / 2 
2
 x − 3l / 2  2 x − 2l   x − 3l / 2  x − l / 2 
        
l2 l2  − l  
2
    l2  l
2

2 l   x − 3l / 2  2 x − 2l   2 x − 2l  x − l / 2  
2
 2 x − 2l 
= EA ∫          dx
0
  l 2
 − l 2
  − l 2
  − l 2
  l 2

 x − 3 l / 2  x − l / 2   2 x − 2l  x − l / 2  x − l / 2 
      
  − l  l
2 2 2 2 2
 l  l  l 
 7 /6 − 4/ 3 1/6 
EA 
= − 4 / 3 8 / 3 − 4 / 3.
l 
 1 / 6 − 4 / 3 7 / 6 
 0 
cxdx = 4cl 2 / 3 .
2l
The element body force vector is f (1)
=∫ N (1)T
0
2cl 2 / 3
The resulting global system of equations is

 7 / 6 − 4 / 3 1 / 6   0   r1 
EA   
 − 4 / 3 8 / 3 − 4 / 3 u2  = 4cl 2 / 3 .
l
 1 / 6 − 4 / 3 7 / 6  u3  2cl 2 / 3
Solving the above equation, we obtain
cl 3 cl 3
u2 = 1.8333 , u3 = 2.6666 , and r1 =-2 cl 2 .
EA EA
 
 0 
 ( x − l )( x − 2l ) x( x − 2l ) x( x − l )   cl 3 
u( x ) = N (1) d (1) =    1.8333 ,
 2l 2 − l2 2l 2
 EA 
 cl 3 
2 .6666
EA 
0.5cl 2 cl 2
=− x + 2.3333 x.
EA EA

cl cl 2
ε (x) = − x + 2.3333 .
EA EA
(d) At the right end with x=2l, for both of two linear displacement elements and a
single quadratic displacement element, the displacements are same to the exact
cl 3
result ( u3 = 2.6666 ). As for the stress and strain, with two linear
EA
cl 2 cl 2
displacement elements, we obtain ε (2l ) = 0.8333 and σ (2l ) = 0.8333 .
EA A

7
 cl 2
With a single quadratic displacement element, ε (2l ) = 0.3333 , and
EA
cl 2
σ (2l ) = 0.3333 . It is found that the approximation result of the stress and
A
strain with a single quadratic displacement element is closer to the exact result
( σ = 0 and ε = 0 ).

(e). For both of two linear displacement elements and a single quadratic
displacement element, the reaction forces acting on node 1 are r1 =-2 cl 2 , which
2l
satisfy the equilibrium equation ( r1 + ∫ cxdx = 0 ). However, as shown in (d), the
0
stress boundary conditions are not satisfied for the two meshes.