Part 1a: Non-dimensionalization and Weak Form
The governing equation
divσ + f = ρu¨ (1)
can be put in its weak form. By redefining the second derivative and converting to index notation, the governing equation becomes
(2)
which can be rearranged to
(3)
By integrating over the body Ω, the integral becomes
(4)
We can multiply by a perturbation vi and apply the chain rule to the first term.
(5)
Applying the divergence theorem yields
(6)
The governing equation in its weak form becomes
Length terms, traction and stress terms, and time terms in the weak form equation can be nondimensionalized by the beam length L, Young’s modulus E, and characteristic time t˜.
Part 1b: Just Displacement
Deflection was determined at the center of a beam with length L = 1m, width W = 0.2m, height H = 0.2m, Youngs Modulus E = 200GPa, Poisson ratio ν = 0.3, and density 7800kg m−3. The beam is clamped at both ends, and a downward force of F = 100N is applied over the area on the top surface near the center, within 0.2 m along the width, and within 0.02 m along the length of the beam.
The deflection was determined to be 3.75 × 10−8 m downward, and the Paraview output is shown in Fig 1
Figure 1: Paraview output for problem 1b
The code used for part 1b can be seen in Appendix 1b.
Part 1c: Free Vibration
The force was removed, and the calculated deflection was set as an initial condition of the beam. A free vibration was modeled, and a time-plot can be seen in Fig 2.
The period of oscillation was determined to be 0.0011 seconds, and the natural frequency was determined to be 5681 rad/s or 904.2 Hz.
Figure 2: Vertical Deflection vs Time
The code used for Problem 1c is shown in Appendix 1c.
Part 1d: Natural Frequency, Changing Dimensions
Changing Width
Next, the beam width was varied to assess how natural frequency changes with beam width. Widths of 0.2, 0.4, 0.6, 0.8, and 1.0 m were assessed. A plot of natural frequencies over beam width is shown in Fig 3.
The Python code used for part 1d-i is shown in Appendix 1d-i.
Figure 3: Natural Frequency over Beam Width for Problem 1d
Changing Height
Next, the beam height was varied to assess how natural frequency changes with beam height. Heights of 0.2, 0.4, 0.6, 0.8, and 1.0 m were assessed. A plot of natural frequencies over beam Height is shown in Fig 4.
The Python code used for part 1d-ii is shown in Appendix 1d-ii.
Figure 4: Natural Frequency over Beam Height for Problem 1d
Appendix 1b: Code for Problem 1b
”””
Python script for Part 1b of Project 2a
Original Author : Vinamra Agrawal
Date : January 25 , 2019
Edited By: Omkar Mulekar
Date : February 10 , 2019
”””
from future import print function from fenics import ∗ import matplotlib
matplotlib . use (”Agg”)
import matplotlib . pyplot as plt from ufl import nabla div import math
#============================================================== # Define System Properties
#============================================================== length = 1;
W = 0.2;
H = 0.2;
a = 0.04/ length ; b = 0.4∗H/length ;
area = a∗b; F = −100
youngs = 200e9 # Youngs nu = 0.3 # Poisson rho = 7800 # Density
# Lame parameters mu = (youngs)/(2∗(1+nu) )
lambda = (nu∗youngs)/((1+nu)∗(1−2∗nu) )
g = 10
traction applied = F/area
#============================================================== # Dimensionless parameters
#============================================================== l nd = length/length
w nd = W/length h nd = H/length
mu nd = mu/youngs lambda nd = lambda /youngs traction nd = traction applied /youngs
#============================================================ # Boundaries and Geometry
#============================================================ mesh = BoxMesh( Point (0 ,0 ,0) , Point ( l nd , w nd , h nd) ,20 ,6 ,6) V = VectorFunctionSpace (mesh , ’P’ ,1) tol = 1E−14
def boundary left (x , on boundary) : return (onboundary and near (x [0] ,0 , tol ) )
def boundary right (x , on boundary) : return on boundary and near (x [0] , l nd , tol )
bcleft = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundaryleft ) bc right = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundary right )
#============================================================ def epsilon (u) : return 0.5∗( nabla grad (u) + nabla grad (u) .T)
def sigma(u) :
return lambda nd∗nabla div (u)∗Identity (d) + mu nd∗( epsilon (u) + epsilon (u) .T)
#============================================================ # First we solve the problem of a cantelever beam under fixed # load .
#============================================================
u init = TrialFunction (V)
d = u init . geometric dimension () v = TestFunction (V)
f = Constant ((0.0 ,0.0 ,0.0) )
T init = Expression (( ’ 0.0 ’ , ’x [0] >= 0.48∗ l && x [0] <= .52∗ l && near (x [1] ,w) && x [2] >= 0.3∗h && x [2] <= 0.7∗h? A : 0.0 ’ , ’ 0.0 ’
) , degree=1, l=l nd , w=w nd ,h=h nd , A=traction nd )
F init = inner (sigma( u init ) , epsilon (v) )∗dx − dot( f , v)∗dx − dot(
T init , v)∗ds a init , L init = lhs ( F init ) , rhs ( F init )
print(”Solving the i n i t i a l cantelever problem”) u init = Function (V)
solve ( a init==L init , u init , [ bc left , bc right ])
w nd = u init ( lnd /2.0 ,w nd/2.0 , h nd /2.0)
w = w nd ∗ length print(w[1])
vtkfile u = File ( ’ deflection . pvd ’ ) vtkfile u << u init
Appendix 1c: Code for Problem 1c
”””
Python script for Part 1c of Project 2a
Original Author : Vinamra Agrawal
Date : January 25 , 2019
Edited By: Omkar Mulekar
Date : February 28 , 2019
”””
from future import print function from fenics import ∗ import matplotlib
matplotlib . use (”Agg”)
import matplotlib . pyplot as plt from ufl import nabla div import math import numpy as np from scipy . signal import argrelextrema
#============================================================== # Define System Properties
#============================================================== length = 1;
W = 0.2;
H = 0.2;
a = 0.04∗ length ; b = 0.4∗H; area = a∗b; F = −100
youngs = 200e9 # Youngs nu = 0.3 # Poisson rho = 7800 # Density
# Lame parameters mu = (youngs)/(2∗(1+nu) )
lambda = (nu∗youngs)/((1+nu)∗(1−2∗nu) )
g = 10
traction applied = F/area
#============================================================== # Dimensionless parameters
#============================================================== lnd = length/length
wnd = W/length hnd = H/length
bar speed = math. sqrt (youngs/rho) t char = length/bar speed t = 0 t i = 0.5
dt = 0.1 num steps = 150
mu nd = mu/youngs lambda nd = lambda /youngs traction nd = traction applied /youngs
#============================================================ # Boundaries and Geometry
#============================================================ mesh = BoxMesh( Point (0 ,0 ,0) , Point ( l nd , w nd , h nd) ,20 ,6 ,6) V = VectorFunctionSpace (mesh , ’P’ ,1) tol = 1E−14
def boundary left (x , on boundary) : return (onboundary and near (x [0] ,0 , tol ) )
def boundary right (x , on boundary) : return on boundary and near (x [0] , l nd , tol )
bcleft = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundaryleft ) bc right = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundary right )
#============================================================ def epsilon (u) :
return 0.5∗( nabla grad (u) + nabla grad (u) .T)
def sigma(u) :
return lambda nd∗nabla div (u)∗Identity (d) + mu nd∗( epsilon (u) + epsilon (u) .T)
#============================================================ # First we solve the problem of a cantelever beam under fixed # load .
#============================================================ u init = TrialFunction (V)
d = u init . geometric dimension () v = TestFunction (V)
f = Constant ((0.0 ,0.0 ,0.0) )
T init = Expression (( ’ 0.0 ’ , ’x [0] >= 0.48∗ l && x [0] <= .52∗ l && near (x [1] ,w) && x [2] >= 0.3∗h && x [2] <= 0.7∗h? A : 0.0 ’ , ’ 0.0 ’
) , degree=1, l=l nd , w=w nd ,h=h nd , A=traction nd )
F init = inner (sigma( u init ) , epsilon (v) )∗dx − dot( f , v)∗dx − dot(
T init , v)∗ds a init , L init = lhs ( F init ) , rhs ( F init )
print(”Solving the i n i t i a l cantelever problem”) u init = Function (V) solve ( a init==L init , u init , [ bc left , bc right ])
w nd = u init ( lnd /2.0 ,w nd/2.0 , h nd /2.0)
w = w nd ∗ length print(w[1])
#============================================================
# Next we use this as i n i t i a l condition , l e t the force go and
# study the vertical vibrations of the beam
#============================================================ u n = interpolate (Constant ((0.0 ,0.0 ,0.0) ) ,V) u n 1 = interpolate (Constant ((0.0 ,0.0 ,0.0) ) ,V) u n . assign ( u init ) u n 1 . assign ( u init ) T n = Constant ((0.0 ,0.0 ,0.0) )
u = TrialFunction (V) d = u. geometric dimension () v = TestFunction (V)
F = (dt∗dt)∗inner (sigma(u) , epsilon (v) )∗dx + dot(u, v)∗dx − (dt∗dt)∗ dot( f , v)∗dx − (dt∗dt)∗dot (T n , v)∗ds − 2.0∗ dot(u n , v)∗dx + dot( u n 1 , v)∗dx
a ,L = lhs (F) , rhs (F)
xdmffileu = XDMFFile( ’ results / solution . xdmf ’ ) xdmffiles = XDMFFile( ’ results / stress . xdmf ’ ) u = Function (V)
u store = [0] ∗ num steps time = [0] ∗ num steps
index = 0 for n in range( num steps ) : print(”time = %.2f ” % t ) Tn . t = t
solve (a == L, u, [ bc left , bc right ]) ugrab = u(0.5 ,0.1 ,0.1) u store [n] = u grab [1]
if (abs(t−index ) <0.01) :
print(”Writing output f i l e s . . . ”)
xdmffile u . write (u∗length , t )
W = TensorFunctionSpace (mesh , ”Lagrange” , 1) stress = lambda ∗nabla div (u)∗Identity (d) + mu∗( epsilon (u)
+ epsilon (u) .T) xdmffile s . write ( project ( stress ,W) , t ) index += 1
time [n] = t t+=dt
u n 1 . assign (u n) u n . assign (u)
# Get period of o s c i l l a t i o n u np = np. array ( u store ) min args = argrelextrema (u np ,np. less ) period = ( time [ min args [ 0 ] [ 1 ] ] − time [ min args [ 0 ] [ 0 ] ] ) ∗t char nat freq = 2∗math. pi /period
print(”Period of Oscillation ” , period , ” seconds”) print(”Natural Frequency : ” , nat freq , ” rad/s”)
plt . figure (1) plt . plot (time , u store ) plt . xlabel ( ’ time [ s ] ’ ) plt . ylabel ( ’ Vertical Deflection [m] ’ ) plt . savefig ( ’1 cfig . png ’ )
Appendix 1d-i: Code for Problem 1d-i
”””
Python script for Part 1d . i of Project 2a
Original Author : Vinamra Agrawal
Date : January 25 , 2019
Edited By: Omkar Mulekar
Date : February 28 , 2019
”””
from future import print function from fenics import ∗ import matplotlib
matplotlib . use (”Agg”)
import matplotlib . pyplot as plt from ufl import nabla div import math import numpy as np from scipy . signal import argrelextrema
#============================================================== # Define System Properties
#============================================================== length = 1;
W1 = 0.2; H = 0.2;
alpha = np. linspace (0.2 ,1 ,num=5) print(”alpha = ” , alpha )
W = alpha ∗ W1
a = 0.04∗ length ; b = 0.4∗H; area = a∗b; F = −100
youngs = 200e9 # Youngs nu = 0.3 # Poisson rho = 7800 # Density
# Lame parameters mu = (youngs)/(2∗(1+nu) )
lambda = (nu∗youngs)/((1+nu)∗(1−2∗nu) )
g = 10
traction applied = F/area
nat freq = [0] ∗ len( alpha )
print(”Beginning Loop . . . ”) for i in range(len( alpha ) ) : print(”alpha = ” , alpha [ i ])
#============================================================== # Dimensionless parameters
#============================================================== lnd = length/length
wnd = W[ i ]/ length hnd = H/length
bar speed = math. sqrt (youngs/rho) t char = length/bar speed t = 0 t i = 0.5
dt = 0.1 num steps = 275
mu nd = mu/youngs lambda nd = lambda /youngs traction nd = traction applied /youngs
#============================================================ # Boundaries and Geometry
#============================================================ mesh = BoxMesh( Point (0 ,0 ,0) , Point ( l nd , w nd , h nd) ,20 ,6 ,6) V = VectorFunctionSpace (mesh , ’P’ ,1) tol = 1E−14
def boundaryleft (x , on boundary) :
return (onboundary and near (x [0] ,0 , tol ) )
def boundary right (x , on boundary) : return on boundary and near (x [0] , l nd , tol )
bcleft = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundaryleft ) bc right = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundary right )
#============================================================ def epsilon (u) : return 0.5∗( nabla grad (u) + nabla grad (u) .T)
def sigma(u) :
return lambda nd∗nabla div (u)∗Identity (d) + mu nd∗( epsilon (u
) + epsilon (u) .T)
#============================================================ # First we solve the problem of a cantelever beam under fixed # load .
#============================================================ u init = TrialFunction (V)
d = u init . geometric dimension () v = TestFunction (V)
f = Constant ((0.0 ,0.0 ,0.0) )
T init = Expression (( ’ 0.0 ’ , ’x [0] >= 0.48∗ l && x [0] <= .52∗ l && near (x [1] ,w) && x [2] >= 0.3∗h && x [2] <= 0.7∗h? A : 0.0 ’ , ’
0.0 ’ ) , degree=1, l=l nd , w=w nd ,h=h nd , A=traction nd )
F init = inner (sigma( uinit ) , epsilon (v) )∗dx − dot( f , v)∗dx − dot
( T init , v)∗ds a init , L init = lhs ( F init ) , rhs ( F init )
print(”Solving the i n i t i a l cantelever problem”) u init = Function (V) solve ( a init==L init , u init , [ bc left , bc right ])
down nd = u init ( l nd /2.0 ,w nd/2.0 , h nd /2.0)
w = down nd ∗ length
print(” I n i t i a l Displacement : ” ,w[1])
#============================================================
# Next we use this as i n i t i a l condition , l e t the force go and
# study the vertical vibrations of the beam
#============================================================ u n = interpolate (Constant ((0.0 ,0.0 ,0.0) ) ,V) u n 1 = interpolate (Constant ((0.0 ,0.0 ,0.0) ) ,V) u n . assign ( u init ) u n 1 . assign ( u init )
T n = Constant ((0.0 ,0.0 ,0.0) )
u = TrialFunction (V) d = u. geometric dimension () v = TestFunction (V)
F = (dt∗dt)∗inner (sigma(u) , epsilon (v) )∗dx + dot(u, v)∗dx − (dt∗ dt)∗dot( f , v)∗dx − (dt∗dt)∗dot (T n , v)∗ds − 2.0∗ dot(u n , v)∗dx
+ dot( u n 1 , v)∗dx a ,L = lhs (F) , rhs (F)
u store = [0] ∗ num steps time = [0] ∗ num steps
index = 0 for n in range( num steps ) : print(”time = %.2f ” % t ) T n . t = t
u = Function (V)
solve (a == L, u, [ bcleft , bc right ]) ugrab = u( l nd /2.0 ,wnd/2.0 , h nd /2.0) ustore [n] = u grab [1]
if (abs(t−index ) <0.01) :
# print (” Writing output f i l e s . . . ” )
W = TensorFunctionSpace (mesh , ”Lagrange” , 1) stress = lambda ∗nabla div (u)∗Identity (d) + mu∗( epsilon (
- u) + epsilon (u) .T)
index += 1
time [n] = t t+=dt
u n 1 . assign (u n) u n . assign (u)
plt . figure (1) plt . plot (time , u store ) plt . xlabel ( ’ time [ s ] ’ ) plt . ylabel ( ’ Vertical Deflection [m] ’ ) plt . savefig ( ’1 dfig test . png ’ )
# Get period of o s c i l l a t i o n u np = np. array ( u store )
min args = argrelextrema (u np ,np. greater ) print(”min args” , min args [0]) period = ( time [ min args [ 0 ] [ 1 ] ] − time [ min args [ 0 ] [ 0 ] ] ) ∗t char nat freq [ i ] = 2∗math. pi /period
print(”Period of Oscillation ” , period , ” seconds”) print(”Natural Frequency : ” , nat freq , ” rad/s”)
plt . figure (2) plt . plot (W, nat freq , ’b−x ’ ) plt . xlabel ( ’Beam Width [m] ’ ) plt . ylabel ( ’ Natural Frequency [ rad/s ] ’ ) plt . savefig ( ’1 dfig alpha . png ’ )
Appendix 1d-ii: Code for Problem 1d-ii
”””
Python script for Part 1d . i i of Project 2a
Original Author : Vinamra Agrawal
Date : January 25 , 2019
Edited By: Omkar Mulekar
Date : February 28 , 2019
”””
from future import print function from fenics import ∗ import matplotlib
matplotlib . use (”Agg”)
import matplotlib . pyplot as plt from ufl import nabla div import math import numpy as np from scipy . signal import argrelextrema
#============================================================== # Define System Properties
#============================================================== length = 1;
W = 0.2; H1 = 0.2; beta = np. linspace (0.2 ,1 ,num=5) print(”beta = ” , beta )
H = beta ∗ H1
print(”H = ” , H) youngs = 200e9 # Youngs nu = 0.3 # Poisson rho = 7800 # Density
# Lame parameters mu = (youngs)/(2∗(1+nu) ) lambda = (nu∗youngs)/((1+nu)∗(1−2∗nu) )
g = 10
nat freq = [0] ∗ len( beta )
print(”Beginning Loop . . . ”) for i in range(len( beta ) ) : print(”beta = ” , beta [ i ])
#============================================================== # Dimensionless parameters
#============================================================== lnd = length/length
wnd = W/length hnd = H[ i ]/ length
bar speed = math. sqrt (youngs/rho) t char = length/bar speed t = 0 t i = 0.5
dt = 0.1 num steps = 275
mu nd = mu/youngs lambda nd = lambda /youngs
F = −100 a = 0.04∗ length ; b = 0.4∗H[ i ] ;
tractionapplied = F/(a∗b) traction nd = traction applied /youngs
#============================================================ # Boundaries and Geometry
#============================================================ mesh = BoxMesh( Point (0 ,0 ,0) , Point ( l nd , w nd , h nd) ,20 ,6 ,6) V = VectorFunctionSpace (mesh , ’P’ ,1) tol = 1E−14
def boundaryleft (x , on boundary) :
return (onboundary and near (x [0] ,0 , tol ) )
def boundary right (x , on boundary) : return on boundary and near (x [0] , l nd , tol )
bcleft = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundaryleft ) bc right = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundary right )
#============================================================ def epsilon (u) : return 0.5∗( nabla grad (u) + nabla grad (u) .T)
def sigma(u) :
return lambda nd∗nabla div (u)∗Identity (d) + mu nd∗( epsilon (u
) + epsilon (u) .T)
#============================================================ # First we solve the problem of a cantelever beam under fixed # load .
#============================================================ u init = TrialFunction (V)
d = u init . geometric dimension () v = TestFunction (V)
f = Constant ((0.0 ,0.0 ,0.0) )
T init = Expression (( ’ 0.0 ’ , ’x [0] >= 0.48∗ l && x [0] <= .52∗ l && near (x [1] ,w) && x [2] >= 0.3∗h && x [2] <= 0.7∗h? A : 0.0 ’ , ’
0.0 ’ ) , degree=1, l=l nd , w=w nd ,h=h nd , A=traction nd )
F init = inner (sigma( uinit ) , epsilon (v) )∗dx − dot( f , v)∗dx − dot
( T init , v)∗ds a init , L init = lhs ( F init ) , rhs ( F init )
print(”Solving the i n i t i a l cantelever problem”) u init = Function (V) solve ( a init==L init , u init , [ bc left , bc right ])
down nd = u init ( l nd /2.0 ,w nd/2.0 , h nd /2.0)
w = down nd ∗ length
print(” I n i t i a l Displacement : ” ,w[1])
#============================================================
# Next we use this as i n i t i a l condition , l e t the force go and
# study the vertical vibrations of the beam
#============================================================ u n = interpolate (Constant ((0.0 ,0.0 ,0.0) ) ,V) u n 1 = interpolate (Constant ((0.0 ,0.0 ,0.0) ) ,V) u n . assign ( u init )
un 1 . assign ( u init ) Tn = Constant ((0.0 ,0.0 ,0.0) )
u = TrialFunction (V) d = u. geometric dimension () v = TestFunction (V)
F = (dt∗dt)∗inner (sigma(u) , epsilon (v) )∗dx + dot(u, v)∗dx − (dt∗ dt)∗dot( f , v)∗dx − (dt∗dt)∗dot (T n , v)∗ds − 2.0∗ dot(u n , v)∗dx
+ dot( u n 1 , v)∗dx a ,L = lhs (F) , rhs (F)
u store = [0] ∗ num steps time = [0] ∗ num steps
index = 0 for n in range( num steps ) : print(”time = %.2f ” % t ) T n . t = t
u = Function (V)
solve (a == L, u, [ bcleft , bc right ]) ugrab = u( l nd /2.0 ,wnd/2.0 , h nd /2.0) ustore [n] = u grab [1]
if (abs(t−index ) <0.01) :
W = TensorFunctionSpace (mesh , ”Lagrange” , 1) stress = lambda ∗nabla div (u)∗Identity (d) + mu∗( epsilon (
- u) + epsilon (u) .T)
index += 1
time [n] = t t+=dt
u n 1 . assign (u n) u n . assign (u)
plt . figure (1) plt . plot (time , u store ) plt . xlabel ( ’ time [ s ] ’ ) plt . ylabel ( ’ Vertical Deflection [m] ’ ) plt . savefig ( ’1 dfig test2 . png ’ )
# Get period of o s c i l l a t i o n u np = np. array ( u store )
min args = argrelextrema (u np ,np. greater ) print(”min args” , min args [0]) period = ( time [ min args [ 0 ] [ 1 ] ] − time [ min args [ 0 ] [ 0 ] ] ) ∗t char nat freq [ i ] = 2∗math. pi /period
print(”Period of Oscillation ” , period , ” seconds”) print(”Natural Frequency : ” , nat freq , ” rad/s”)
plt . figure (2) plt . plot (H, nat freq , ’b−x ’ ) plt . xlabel ( ’Beam Height [m] ’ ) plt . ylabel ( ’ Natural Frequency [ rad/s ] ’ ) plt . savefig ( ’1 dfig beta . png ’ )
Reviews
There are no reviews yet.