MAE 4700/5700 (Fall 2011) Exam 1 Solutions

October 18th, 7:30-9:30 pm

Problem 1 (35 points) – Consider a one-dimensional elasticity problem with axial

deformation as shown in the figure below. The bar is constrained at both ends (A and C).
Its cross sectional area is constant (A = 0.1 m2 ) on segment AB and varies linearly A =
0.05(x - 1) on BC. The Young’s Modulus is E = 2⋅107 Pa. A distributed load b = 10 N/m
is applied along the left portion of the bar AB and a point force P = 150 N acts at point B.
Use a three-node element on AB and a two-node element on BC as shown below.
Consider the global node numbering as shown. The dimensions in the figure are in
meters.

a) Provide a precise weak statement for the problem and based on this give the final
expressions for the element stiffness matrix and load vector.
b) For each element, compute the element force vector Fe and element stiffness
e
matrix K and then provide the assembled F and K.
c) Apply the essential boundary conditions and compute the finite element nodal
displacements.
d) Provide expressions for the finite element stress in each element.

Possibly useful formulas: Using standard notation, the basis functions in x-coordinate
1 e
Ne =  x2 − x x − x1e  .
Le 
for a two-node linear element are:

Similarly for a 3-node quadratic element

2
Ne= 2
( x − x2e )( x − x3e ) −2( x − x1e )( x − x3e ) ( x − x1e )( x − x2e ) 
Le

Finite Element Analysis for Mechanical & Aerospace Design Page 1 of 10

MAE 4700/5700 (Fall 2011) Exam 1 Solutions
Solution:
(a) The weak statement is:
Find a solution u ( x) with u= ( xc ) 0 such that for any w( x) ∈ H 1 with w=
( x A ) u= ( x A ) w=
( xc ) 0
the following equation holds:
xc du dw xc
∫ xA
EA
dx dx
dx = ∫ bw dx
xA

Based on this weak form, the element stiffness matrix and load vector are
=K e EA∫ BT Bdx
= , fe ∫ bN
T
dx + N T P |concentrated
load
location

where B and N follow expressions derived in class.

(b) For element 1:

1 1 
N 1 =  ( x − 2)( x − 3) −( x − 1)( x − 3) ( x − 1)( x − 2) 
2 2 
dN 1  5 3
B1 = =  x− 4 − 2x x − 
dx  2 2
Therefore, the element stiffness matrix is
 2
5  5  5  3 
 x −   x −  ( 4 − 2 x )  x −  x −  
 2  2  2  2 
3  3 
( 4 − 2x) ( 4 − 2 x )  x −  dx
3
E1 A1 ∫ B1 B1dx =
K1 = 2 × 106 ∫  •
T 2
1 1
  2 
 2 
  3 
• • x− 
  2 
 2.3333 −2.6667 0.333 
= 10 ×  −2.6667 5.333 −2.6667 
6

 0.333 −2.6667 2.3333 

The force vector is

1  1 
 2 ( x − 2 ) ( x − 3)   2 ( x − 2 ) ( x − 3)   3.333 
3 3     
f 1 ∫ N 1 bdx + N 1 P |x=
= ∫1 − − − + − − − 0 |=
T T

1 x=B
10  ( x 1)( x 3)  dx  ( x 1)( x 3)  x 3 13.333
1  1   3.333 
 ( x − 1) ( x − 2)   ( x − 1) ( x − 2) 
2  2 

For element 2:

Finite Element Analysis for Mechanical & Aerospace Design Page 2 of 10

MAE 4700/5700 (Fall 2011) Exam 1 Solutions
1
N 2 = [5 − x x − 3]
2
 1 1
B2 = −
 2 2 
Then the element stiffness matrix is

1  −1 7 1
51  1 −1
  0.05( x − 1)(210 ) ( −1=
1) dx 106 ∫ 
5 5
= ∫ B 2 B 2 E 2 A= ∫ ( x − 1)dx
T
K2
1 
( x)dx
2 1  3 4 −1
3 3 2 
 1.5 −1.5
= 106 ×  
 −1.5 1.5 

The load vector is

1 5 − x 
( ) 150 
3
= ∫ dx + N 2 P
= = =
× 150  
T T T
f2 N 2 b = |x x= N2 P |x 3  
1
=0
B
2  x − 3 x = 3  0 

Thus, the global K is

 2.3333 0 −2.6667 0.3333 
 0 1.5 0 −1.5 
K 106 × 
=
 −2.6667 0 5.3333 −2.6667 
 
 0.333 −1.5 −2.6667 3.8333 
The global f is
 3.333 
 0 
f = 
 13.333 
 
153.333

(c) Apply the essential boundary conditions, we have

 2.3333 0 −2.6667 0.3333   u1 = 0   r1 + 3.333
 0
 1.5 0 −1.5  u2 = 0   r2 

10 ×
6
=
 −2.6667 0 5.3333 −2.6667   u3   13.333 
    
 0.333 −1.5 −2.6667 3.8333   u4   153.333 

where r1 and r2 are reaction forces.

 5.333 −2.6667  u3   13.333  u3  3.45
⇒ 106     =   ⇒   =   × 10−5 m
 −2.6667 3.8333  u4  153.333 u4   6.4 
The reaction forces can be computed as
 −2.6667 0.333 u3   r1 + 3.333  r1   −74.02 
106  =
 u    ⇒=   −96  N
 0 − 1.5  4  r2   2 
r 
(d) The finite element stresses in each element are

Finite Element Analysis for Mechanical & Aerospace Design Page 3 of 10

MAE 4700/5700 (Fall 2011) Exam 1 Solutions
 u1 = 0 
 3 −5 
EB {d } =
5
σ ( x) =
1 1
2 × 10  x −
1 7
4 − 2x x −   u3 =3.45 × 10 
 2 2
 u=4 6.4 × 10−5 
= 840 − 100 x Pa

 1 1   u= 6.4 × 10−5 
EB 2 {d 2 } =
σ 2 ( x) = 2 × 107  − 
4
=−640 Pa
 2 2   u2 = 0 

Problem 2 (10 points) – Examine if the following mapping is invertable.

Solution:
Using the four basis functions,
1 1
N1 (ξ ,η ) = (1 − ξ )(1 − η ), N 2 (ξ ,η ) = (1 + ξ )(1 − η )
4 4
1 1
N 3 (ξ ,η ) = (1 + ξ )(1 + η ), N 4 (ξ ,η ) = (1 − ξ )(1 + η )
4 4
the mapping is
4
3
x =∑ x j N j =3 N 3 + 3 N 4 = (1 + η )
j =1 2

4
1
y = ∑ y j N j = N 2 + N3 = (1 + ξ )
j =1 2

The Jacobian for this mapping is

3
0
2 = 3
| J |= − <0
1 4
0
2

Therefore, this mapping is not invertable.

Finite Element Analysis for Mechanical & Aerospace Design Page 4 of 10

MAE 4700/5700 (Fall 2011) Exam 1 Solutions
Problem 3 (5 points) – Verify that the shape functions of four-node elements derived in
class satisfy the following condition
Nnodes = 4
 i (ξ ,η ) = 1
∑i =1
N

Explain why the above condition is satisfied (and not only for this element!).

Solution:
According to the basis function defined in problem 2, we have
4
1 1 1 1
∑ N (ξ ,η )=
i =1
i
4
(1 − η − ξ + ξη ) + (1 − η + ξ − ξη ) + (1 + η + ξ + ξη ) + (1 + η − ξ − ξη )
4 4 4
1 1 1 1
= + + + =1
4 4 4 4

Generally, the above condition is always satisfied. Let’s consider a constant function
f (ξ ,η ) = a . By FEM interpolation, this function value at any point inside the element is
Nnodes Nnodes
=f i
=i 1 =i 1
∑ a N (ξ ,η ) a
= ∑ N i (ξ ,η )

Since f is a constant field, and we use at least linear interpolation, this interpolation
should be exact, i.e. f= f= a .
Therefore, we have
Nnodes Nnodes
a= a
=i 1 =i 1
∑ N i (ξ ,η ) ⇒ ∑ N i (ξ ,η ) = 1

The key idea here is that the FEM basis functions used form a complete set of functions –
can represent both rigid body motion and `constant strain’. It is the completeness of the
FE basis that makes this equation true.

Problem 4 (10 points) – For the 4-node element Ω shown in the figure below, compute
2
∂N
at the point .
(ξ ,η ) = (0.5,0.5) in Ω
∂y

Finite Element Analysis for Mechanical & Aerospace Design Page 5 of 10

MAE 4700/5700 (Fall 2011) Exam 1 Solutions
Solution:

The mapping here is given as:

1
x = N 2 + 2 N3 = ( 3 + 3ξ + η + ξη )
4
1
y = N 4 + 2 N3 = ( 3 + ξ + 3η + ξη )
4
Therefore, the Jacobian matrix is
 ∂x ∂x  7 3
 ∂ξ ∂η  1 3 + η 1 + ξ  8 8
=J =   =  
 ∂y ∂y  4 1 + η 3 + ξ=ξ 0.5,
= η 0.5 3 7
 ∂ξ ∂η   8 8 
 
5
The determinate is | J |= .
8
From the formulae derived in class, we have
∂ξ 1 ∂x 3 ∂η 1 ∂x 7
=− =
− , = =
∂y | J | ∂η 5 ∂y | J | ∂ξ 5
∂N 2 1 1 ∂N 2 1 3
= (1 − η ) = , − (1 + ξ )
= =
−
=∂ξ 4 η =
0.5 8 ∂η 4 ξ 0.5 8

Finally, we can get

∂N 2 ∂N 2 ∂ξ ∂N 2 ∂η 1  3   3  7 3
= + = ×  −  +  −  × =−
∂y ∂ξ ∂y ∂η ∂y 8  5   8  5 5

Problem 5 (40 points) - Consider a heat conduction problem on a rectangular (2 x 1) m

domain as shown in the Figure below (the governing equations are summarized for your
convenience). The conductivity is k=4 WoC-1, T = 10 oC is prescribed along edge BC.
Edges AB and AD are insulated, i.e. q = 0 W m-1 ; along edge DC the boundary flux is
q = 30 W m-1 . A constant heat source is given f = 50 W m-2 . We will analyze this
problem using one 4-node finite element. Use as (local and global) node numbering the
one shown on the right of the Figure.

Finite Element Analysis for Mechanical & Aerospace Design Page 6 of 10

MAE 4700/5700 (Fall 2011) Exam 1 Solutions

a) Provide a precise weak statement for the problem and based on this give the
final expressions for the element stiffness matrix and load vectors.
b) Compute the element stiffness matrix. This will require a number of tasks
summarized below
1. Using the basis functions, give an expression in terms of (ξ ,η ) for
 ∂x e ∂y e 
the Jacobian matrix  
 ∂ξ ∂ξ 
J =
e
 ∂x e ∂y e 
 
 ∂η ∂η 
2. Compute expressions for the determinant of the Jacobian and the
inverse of this Jacobian matrix.
3. Compute the matrix of the derivates of the shape functions wrt natural
coordinates, i.e.
∂N1e ∂ N e2 ∂ N  3e ∂ N e4 
 
 ∂ξ ∂ξ ∂ξ ∂ξ 
 e 
∂N1 ∂ N  e2 ∂ N  3e ∂ N e4 
 ∂η ∂η ∂η ∂η 

 ∂T e 
 
∂x
Give an expression for the B matrix relating   and the
e
4.
 ∂T e 
 
 ∂y 
nodal temperatures.
5. Use Gauss quadrature to compute the stiffness matrix and load vector.
These will also be your global stiffness and load vector.

Finite Element Analysis for Mechanical & Aerospace Design Page 7 of 10

MAE 4700/5700 (Fall 2011) Exam 1 Solutions
Note: The location of the 4 Gauss points for 2D integration is
1 1
(ξ ,η ) =
(± , ± ) and each weight is equal to 1.
3 3
c) Compute the element load vector f e . There are three contributions – one
from the heat source f Ωe , another f qe from the natural boundary condition
and the third (unknown) from the normal heat flux at the nodes with
prescribed temperature (`reaction fluxes’).

1. Compute the first term f Ωe either analytically or using 2D Gauss

integration.

2. Compute the load contribution f qe from the flux boundary

condition. This requires boundary (one-dimensional) integration that
for this problem you can perform analytically.

d) Apply essential BCs in the final system of equations and solve for the
unknown nodal temperatures.
1 1
e) Compute the heat flux components at the Gauss point ( ,− ).
3 3

Solution:

(a) The weak statement for this problem takes the form:
Find a function T ( x, y ) ∈ H 1 ( Ω ) , with =
T ( s ) T ( s ), s ∈ ΓT such that for any w ∈ H 1 ( Ω )
with w = 0 on ΓT the following equations holds:

∫ Ω
k ∇T ⋅ ∇=
wd Ω ∫Ω
fwd Ω − ∫ qwd Γ
Γq

From the weak form, we can easily identify:

=Ke ∫Ω e
k BT Bd=
Ω, f e ∫Ω e
NT f dΩ − ∫ e NT q dΓ
Γq

(b) 1. The basis functions in natural coordinates are:

1 1
N1 (ξ ,η ) = (1 − ξ )(1 − η ), N 2 (ξ ,η ) = (1 + ξ )(1 − η )
4 4
1 1
N 3 (ξ ,η ) = (1 + ξ )(1 + η ), N 4 (ξ ,η ) = (1 − ξ )(1 + η )
4 4
Then the mapping is

Finite Element Analysis for Mechanical & Aerospace Design Page 8 of 10

MAE 4700/5700 (Fall 2011) Exam 1 Solutions
4

∑ xi Ni =
x= 1 −η
2 N1 + 2 N 2 =
i =1
4
1
y = ∑ yi N i = N 2 + N 3 = (1 + ξ )
i =1 2
The Jacobian matrix can be obtained as
 ∂x ∂y 
 ∂ξ ∂ξ   1
=J =   0 2
 ∂x ∂y   
 ∂η ∂η   −1 0 
 
1
2. The determinant of the Jacobian is | J |= .
2
The inverse of this Jacobian matrix is
 1
−1  0 −   0 −1
= 2=
  2 0 
J 2

 1 0 

 ∂N1 ∂N 2 ∂N 3 ∂N 4 
 ∂ξ ∂ξ ∂ξ ∂ξ  1 η − 1 1 − η 1 + η −η − 1
 =
3.  ∂N
1 ∂N 2 ∂N 3 ∂N 4  4 ξ − 1 −1 − ξ 1 + ξ 1 − ξ 
 ∂η ∂η ∂η ∂η 


 ∂N1 ∂N 2 ∂N 3 ∂N 4 
 ∂ξ  1  1 − ξ 1+ ξ −1 − ξ ξ −1 
−1  ∂ξ ∂ξ ∂ξ

4. B J=  2 η − 1 2 1 − η 2 1 + η −2 η + 1 
 ∂N1 ∂N 2 ∂N 3 ∂N 4  4  ( ) ( ) ( ) ( )
 ∂η ∂η ∂η ∂η  


5. The stiffness matrix is

∑∑ k ( B B | J |)
1 1 2 2
= ∫ kB=
Bd Ω ∫ ∫ | d ξ dη
kB B | J=
T T T
K WiW j
Ω −1 −1 ξ ξi=
= ,η ηi
=i 1 =i 1

 3.3333 −2.3333 −1.6667 0.6667 

 −2.3333 3.3333 0.6667 −1.6667 
= 
 −1.6667 0.6667 3.3333 −2.3333
 
 0.6667 −1.6667 −2.3333 3.3333 

(c) 1.

Finite Element Analysis for Mechanical & Aerospace Design Page 9 of 10

MAE 4700/5700 (Fall 2011) Exam 1 Solutions
 25
 25
f N | J | d ξ dη  
1 1
f Ω ∫=
−1 ∫−1
= e T
 25
 
 25
2. We now need to compute the contribution to element load from the natural boundary
conditions at η = 1 .
 0 
 0   0 
   0 
1   2−0
− ∫ q ( N | J |) d ξ =  
1
fq =
e T
−30 ∫  (1 + ξ ) 
1 dξ =
−1 η =1 −1 2  −30 
2   
1   −30 
 2 (1 − ξ ) 
(d) Assemble all the node contributions give:

 3.3333 −2.3333 −1.6667 0.6667   T1 = 10   r1 + 25

   = 10   
 −2.3333 3.3333 0.6667 −1.6667  T=2   r2 − 5  ⇒
 −1.6667 0.6667 3.3333 −2.3333  T3   −5 
    
 0.6667 −1.6667 −2.3333 3.3333   T4   25 

 3.3333 −2.3333 T3   −5  −2.3333 −1.6667  10   5 

 −2.3333 3.3333  T  = −   =  ⇒
   4   25   −1.6667 −2.3333 10  35

T3  17.3529  
T  =  22.6471 C
 4  

 1 1 
(e) The heat flux components at the Gauss point  ,−  are
 3 3
 ∂T   T1   10 
 − k ∂x  T  
 0.1057 0.3943 −0.3943 −0.1057   10 

q=  = − kB 1  2 = − 4  −0.7887 0.7887 0.2113 −0.2113 17.3529 
 − k ∂T  T3 
1
ξ = ,η = −
 
 ∂y 
3 3
   
T4   22.6471
16.9434 
= W /m
 4.4752 

Finite Element Analysis for Mechanical & Aerospace Design Page 10 of 10