University of Wollongong MECH419/928 Assignment 1: Simulation Project

This assignment presents a comprehensive solution to a simulation project in MECH419/928, focusing on finite element analysis (FEA) of a plate. The solution details the application of FEA principles, including stress and strain calculations, stiffness matrix determination, and the use of MATLAB for simulation. The assignment covers theoretical calculations, a MATLAB program, and an Excel spreadsheet with detailed calculations for plate force constants, stresses, and strains. It also includes an analysis of the constraints, governing equations, and a discussion of the Finite Element Analysis method and its application. The solution uses triangular elements for modeling and references relevant literature. Finally, the assignment provides an analysis of the results and conclusions regarding the simulation and analysis.

MECH419/928 1
Name
Institution

MECH419/928 2
Simulation assignment
Part a
The stresses and strains on common edges, including the ‘cut’ edges are continuous and uniform
across the boundaries (L18 and L34). The plate is transversely homogeneous.Forces acting
perpendicular to the plane are zero or may be neglected, as the ratio of thickness to a
fundamental x or y of any section is < 0.1. The plate is, assumed transversely homogeneous.The
plate may be modeled using a quarter sections because of reflection symmetry on the x- and y-
axes.
o Each edge and shared node experiences equal displacement (e.g., L28, L36, L18, or
L34).
o The shape functions, N, must sum to 1[Nh + Ni + Nj = 1].
The plate is free floating. The distributed load on the annulus affects it uniformly.Each edge and
shared node experiences equal displacement.A06 = A01, A02 = A04, and A03=A05 by symmetry. If L
denotes the line between two nodes: L18 = L27, L12 = L34 = L78, L26 = L35, and L23 = L56 by
symmetry.
Displacement Matrix, {d}, and Displacement function, 
For each node h, i, or j on the triangular element, the displacement
vector, {d} is { d }=
{dh
di
d j
}, where {dh }= {uh
vh }, {di }= {ui
vi },∧ {d j }= {u j
v j }.
The linear displacement function for triangular element n, n, is
Ψ n = {u (x , y )
v (x , y ) }=
{
a1
a2
a1 xh
a2 xi
a1 yh
a2 y j
a3 a3 x j a3 y j
a4
a5
a6
a4 xh
a5 xi
a6 x j
a4 yh
a5 yi
a6 y j
}. This gives {u} and {v}:, where ( xh , y h ), ( xi , yi ), and
( x j . y j ) are the coordinates of node h, i, or j. Thus {u}= {uh
ui
u j
}=
{1 xh yh
1 xi yi
1 x j y j
}{a1
a2
a3
} and

MECH419/928 3
{v} = {
vh
vi
v j
}=
{
1 xh yh
1 xi yi
1 x j y j
}{
a4
a5
a6
} , where { Χ } =
{
1 xh yh
1 xi yi
1 x j y j
}= {1 X Y }. If we define
{ a }=
{a1
a2
a3
}∧ {a ' }=
{a4
a5
a6
}, then { u } = { X } { a } and { v }= { X } { a' }. Transposing, { a }= { Χ }−1 { u } and
{ a ' }= { Χ }−1 { v }.
The inverse of {X} is { Χ }−1= 1
2 A {α h α i α j
βh βi β j
γ h γi γj
}, where
{ α }=
{αh
α i
α j
}=
{ xi y j − yi x j
xh y j − yh x j
xh yi − yh xi
}, { β }=
{βh
βi
β j
}=
{ yi− y j
y j− yh
yh− yi
},∧ { γ }=
{γ h
γ i
γ j
}=
{ x j−xi
xh−x j
xi−xh
}, and
2 A=xh ( yi− y j ) + xi ( y j− yh ) + x j ( yh− yi ), where A is area of the triangular element
and the (xm, ym) represent the Cartesian coordinates of node “m”.
Substituting, the expressions for u(x,y) and v(x,y) become:
u ( x , y ) = { Nh uh+ N i ui + N j u j } and v ( x , y ) = { Nh vh +Ni vi +N j v j }, where Nh, Ni, and Nj are
the shape functions.
Nh= 1
2 A ( αh + βh x + γh y ), Ni= 1
2 A ( αi + βi x +γ i y ), and N j= 1
2 A ( α j +β j x+ γ j y ). This gives i
in terms of N and d: {i}= [N]{d}, [ N ] = [ N h 0 Ni 0 N j 0
0 Nh 0 Ni 0 N j ]. Determining
equations for {σ x
σ y
τ xy
}, using the definition of the components of [N}, We now
write:

MECH419/928 4
ε x= ∂ u
∂ x = ∂ Nh
∂ x uh+ ∂ Ni
∂ x ui + ∂ N j
∂ x u j,
ε y= ∂ v
∂ y = ∂ Nh
∂ y vh+ ∂ Ni
∂ y vi+ ∂ N j
∂ y v j, and
γ xy= ∂u
∂ y + ∂ v
∂ x = ( ∂ N h
∂ y uh + ∂ Ni
∂ y ui+ ∂ N j
∂ y uj )+( ∂ Nh
∂ x vh + ∂ Ni
∂ x vi + ∂ N j
∂ x v j ). The
derivatives of the Ns for node
m are given by:
∂ N m
∂ x =
βm
2 A ∧∂ Nm
∂ y = γm
2 A
which can be used to construct the matrix,
[B] ¿ ¿ where { Bm }= 1
2 a {βm 0
0 γ m
γm βm
} for each node m. This gives
{ εm }= { ε x
ε y
γ xy
}=
{ ∂u
∂ x
∂ v
∂ y
∂ u
∂ y + ∂ v
∂ x
}= { Bm } { dm }, for each node
m, and
{ σ }=
{σ x
σ y
τ xy
}= [ D ] [ B ] {d }.
Theoretical claculations
I10= 1
12∗b h3= 1
12∗0.75∗0.0613=0.000014 ∈4 ¿
I 12= 1
12∗b h3= 1
12∗0.752∗0.0623=0.000015∈4 ¿
I 15= 1
12∗b h3= 1
12∗0.75∗0.0623 =0.000015∈4 ¿
I 18= 1
12∗b h3= 1
12∗0.75∗0.0623=0.000015∈4 ¿

MECH419/928 5
Pcr. Exp= 4 π2 ( 29∗106 ) ( 0.0000144 )
102 =160.282=712.96 N
Pcr. Exp= 4 π2 ( 29∗106 ) ( 0.0000154 )
122 =119.258 PSI =530.486 N
Pcr. Exp= 4 π2 ( 29∗106 ) ( 0.0000154 )
152 =76.32 PSI =339.48 N
Pcr. Exp= 4 π2 ( 29∗106 ) ( 0.0000154 )
182 =53.00 PSI =235.75 N
σ cr .10=119.282
0.00575 =3503.43 PSI
σ cr .12= 119.282
0.046624 =2557.87 PSI
σ cr.15 = 76.32
0.0465 =1641.29 PSI
σ cr .18= 53
0.0465 =1139 PSI
Part b: MATLAB program

MECH419/928 6
Part c

MECH419/928 7
The stiffness matrix is defined as {F} = [K]{d} where [K] describes the
global stiffness matrix: [ K ] = ∑
n=1
N
kn. The kn are determined from [ k ]=tA [ B ] T [ D ] [ B ],
where t = object thickness, A = element area, and [B], [D], and [B]t can be
calculated.
The potential energy function, comprised of strain energy and
external forces energy, is: Π=U −W , -W = + b + p + s, where b, p, and s
refer to body forces, concentrated forces, and distributed forces,
respectively. The strain energy is the integral of the normal, the shear strain
matrix, {}, and stress matrix,
{}: U =1
2 ∫
V
❑
{ ε }
T { σ } dV =1
2 ∫
V
❑
{ ε }
T [ D ] { ε } dV .
the potential energy becomes Ωs=−∫
Ω
❑
{ Ψ }T {T } dS, where T represents the
surface tractions or forces that act on the plate boundary (in our case, the
distributed load) and  describes the general displacement functions.
Substituting, into the total potential energy:
π p= 1
2 ∫
V
❑
{d }T [ B ] T [ D ] [ B ] { d } dV −∫
V
❑
{ d }T [ N ]T { X } dV − { d }T { P } − { d }T
∫
S
❑
[ N ]T {T } dS.
Note that {d} is not dependent on the node (x,y) coordinates and can be
pulled out of the integrals, giving
π p= 1
2 {d }T
∫
V
❑
[ B ] T
[ D ] [ B ] dV { d }− { d }T
(∫
V
❑
[ N ] T { X } dV + { P } +∫
S
❑
[ N ] T { T } dS ).

MECH419/928 8
The last three terms in the previous expression define the nodal force
vector, f,.
{ f }=∫
V
❑
[ N ]T { X } dV ± { P }+∫
S
❑
[ N ]T { T } dS.
Writing potential energy as π p= 1
2 {d }T
(∫
V
❑
[ B ]T
[ D ] [ B ] dV ) { d }− { d }T
{f }, it can be
minimized to determine the force vector ( ∂ π p
∂ { d } =0 ¿. The expression for { f } is
{ f }=∫
V
❑
[ B ] t [ D ] [ B ] dV { d }, where the stiffness matrix, k, is defined: [ k ]=∫
V
❑
[ B ] t [ D ] [ B ] dV
.
For a body with constant thickness, t can be pulled out as a constant.
The integral over the volume, V, becomes the integral over area, A, and the
expression for [k] becomes [ k ]=t∫
A
❑
[ B ]t [ D ] [ B ] dXdY . The two dimensional
integral has no x or y dependence except dXdY, so the integral just adds
area term to the integrand. The new expression becomes [ k ]=t [ B ]t [ D ] [ B ] A,
where A is the element area.
The [k] matrix for all three nodes is:
[ k ]=
[ [ khh ] [ khi ] [ khj ]
[ kih ] [ kii ] [ kij ]
[ k jh ] [ k ji ] [ k jj ] ], where the expression for any element in the matrix
follows this convention for solution: [ k hi ] =t [ Bh ] t [ D ] [ Bi ] A.

MECH419/928 9
Once the [k]s are calculated for each element, the global [K] may be
determined as the sum of the [k] for each node: [ K ] =∑
m=1
n
[k ]m. Finally, the
expression for {f} becomes simplified: { f }=[k ] { d } and { F }=∑
m=1
n
{ f }m.
When a rectangular plate is subjected to a uniformly distributed load p (x, y) =Po
Mid span deflection (UDL) = 5qL4/ (384EI)
The load at point D is at the center hence deflection calculated as under;
Since this is a uniformly-distributed load we will have the equations:
The dry unit of the specimen is determined using the following formulae
γd =( M ¿¿ 1−M2)
3.808 V ¿
The modulus of elasticity is taken to be equal to the stiffness of the material and it is expressed
using the following equation.
E= ∆ σ
∆ ε ……………………………….Equation 1
Where σ is the stress on the cantilever and ε is the displacement. The stress of the cantilever at
any point is mathematically expressed asσ = MC
I , σ = 6 PL
b t2
……………………………….Equation 2
The dry unit for the of the specimen for the three tests is as shown below
Test 1

MECH419/928 10
γd = ( 498.16−336.36 )
3.308 × 2.36× 2.36 ×0.998 =110.846
Test 2
γd = ( 494−319.6 )
3.308 × 2.36× 2.36 ×1.05 =113.5612372
Test 3
γd = ( 377.06−204.31 )
3.308 × 2.36× 2.36 ×1.059 =111.5308545
Average γd= 111.98
The maximum and minimum stresses are found by the following:
σ max= σ x +σ y
2 + √ ( σ x−σ y
2 )2
+τ xy
2 andσ min= σ x +σ y
2 − √ ( σx −σ y
2 )2
+ τxy
2 ,
and may be derived from the original Cauchy tensor. These are the
eigenvalues or invariants and may be used to calculate the von Mises stress
(leads to an estimate of safety factor). After the method discussed in
Rathkjen (1994), as well as the treatment in Decamp (2013), the von Mises
stress will be used to draw contour lines to capture regions of constant
stress. The von Mises stress can be evaluated using the principal stresses,
in the plane stress case, as: σ v=σ max
2 −σmax σ min+ ρmin
2 (Rathkjen 1994, Decamp
2014, Nasim 2014).
Part d
Excel spreadsheet attached with calculations for the plate force
constants, stresses and strains.

MECH419/928 11
Relevant Equations and Constants, Initial quantities
Statement for all Elements 01-06: defining base of triangle as bbq and
height of triangle as hhq, where q is the Element number. Area of a right
triangle: A = 0.5 (bbq)( hhq).
Node 01 - description
 Element 01 is a right triangle bounded by nodes 1, 2, and 8
(moving CCW) .
o Node 1: (0,35.35534), Node 2: (50,35.35534), and Node 8:
(0,85.35534)
 Area by two methods:
o A01( = A06):
 base, bb01 = 50mm, height hh01 = 50mm
 A01 = 0.5(50)(50) = 1250 mm2.
o
2 A=xh ( yi− y j ) +xi ( y j− yh ) + x j ( yh− yi ) =0+2500+ 0=2500 mm2 → A=1250 mm2
- ˇ!
 On Element 1, the nodes are, in CCW order, 1, 2, and 8, as
previously stated.
 [ D ] =2.2 x 105 MPa { 1 0.3 0
0.3 1 0
0 0 0.35 }= {2.2 x 105 6.6 x 104 0
6.6 x 104 2.2 x 105 0
0 0 7.7 x 104 }

MECH419/928 12
 [ B ] =
{
{ Bh }
{ Bi }
{ Bj } } where { Bq } = 1
2 a { βq 0
0 γ q
γ q βq
} (see spreadsheet)
o { α01 }=
{ xi y j− yi x j
xh y j− yh x j
xh yi− yh xi
}= { 4267.767
0
−1767.767 }, see spreadsheet for more detail
o { β01 }=
{ yi− y j
y j− yh
yh− yi
}= {−50
50
0 }
o { γ01 }=
{x j −xi
xh−x j
xi−xh
}= {−50
0
50 }
o { Bh }= 1
2 a {βh 0
0 γh
γh βh
}= {−0.02 0
0 −0.02
−0.02 −0.02 }
o { Bi }= 1
2 a {βi 0
0 γi
γi βi
}= { 0.02 0
0 −0.02
−0.02 0.02 }
o { B j }= 1
2 a { β j 0
0 γ j
γ j β j
}= { 0 0
0 0.02
0.02 0 }, see spreadsheet summary
matrices
 { Χ }=
{1 xh yh
1 xi yi
1 x j y j
}= {1 0 35.35534
1 50 35.35534
1 0 85.35534 }
 { X }01
−1= 1
2 A {αh αi α j
βh βi β j
γh γi γ j
}=¿