[SOLVED] AERO4630 Project3- Non-dimensionalization and Weak Form

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 v_i 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⁻³. 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⁻⁸ 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 (

1. 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 (

1. 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 ’ )