# [Solved] AERO4630 Project3- Non-dimensionalization and Weak Form

Click And Check all our Assignments

10 USD \$

Category: Tag:

# 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) :

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

## Reviews

There are no reviews yet.

Only logged in customers who have purchased this product may leave a review.

Shopping Cart
[Solved] AERO4630 Project3- Non-dimensionalization and Weak Form
10 USD \$