Part 1a: Plate Clamped at One EndVibrational analysis was performed on a plate of length L = 1m, width W = 1m, and height H = 0.001m. In this first part, the plate is clamped at one end, and an initial force of F = 100N is applied to the other end over the area 0.98L x L, 0.48W y 0.51W. To display the vibration of the plate, the displacements at several points on the plate are plotted over time. The respective natural frequencies at each point are also plotted. See Figures 1-4.Please note that the code finds the natural frequencies by determining the indices of each relative maximum in the list of displacements. Then, it calls times using the first two indices and subtracts them to obtain the period. Then the natural frequency n = 2pi/T. Because the plate has so many degrees of freedom, there are many modes of vibration. In early cycles, there are multiple dominating modes of vibration, and the method by which the code used in this project determines natural frequency can find multiple frequencies. Thus, plots of natural frequency presented in this document may show unexpected jumps in frequency between points on the plate.Part 1b: Plate Clamped at One End (but bigger)The same analysis was performed for a plate of length L = 2m and width W = 2m. See Figures 5-8. The amplitudes appear to have doubled, while the natural frequencies are the same.
Figure 1: 1A: Displacements along width of 1m plate at L = L
Figure 2: 1A: Natural Frequencies along width of 1m plate at L = L
Figure 3: 1A: Displacements along width of 1m plate at L = 0.5L
Figure 4: 1A: Natural Frequencies along width of 1m plate at L = 0.5L
Figure 5: 1B: Displacements along width of 2m plate at L = L
Figure 6: 1B: Natural Frequencies along width of 2m plate at L = L
Figure 7: 1B: Displacements along width of 2m plate at L = 0.5L
Figure 8: 1B: Natural Frequencies along width of 2m plate at L = 0.5LPart 1c: Plate Clamped at Two EndsThe 1m plate was then clamped at 2 ends (x = 0 and x = L). See Figures 9-12. The amplitudes can be seen on the displacement plots.
Figure 9: 1C: Displacements along width of 1m plate at L = LPart 1d: Plate Clamped on All EndsVibrational analysis was then performed for plates of length/width 1m and 2m. See Figures 13-20. Again, the amplitudes appear to double while the frequencies do not change from the 1m to the 2m simulations.
Figure 10: 1C: Natural Frequencies along width of 1m plate at L = L
Figure 11: 1C: Displacements along width of 1m plate at L = 0.5L
Figure 12: 1C: Natural Frequencies along width of 1m plate at L = 0.5L
Figure 13: 1D: Displacements along width of 1m plate at L = L
Figure 14: 1D: Natural Frequencies along width of 1m plate at L = L
Figure 15: 1D: Displacements along width of 1m plate at L = 0.5L
Figure 16: 1D: Natural Frequencies along width of 1m plate at L = 0.5L
Figure 17: 1D: Displacements along width of 2m plate at L = L
Figure 18: 1D: Natural Frequencies along width of 2m plate at L = L
Figure 19: 1D: Displacements along width of 2m plate at L = 0.5L
Figure 20: 1D: Natural Frequencies along width of 2m plate at L = 0.5LAppendix 1a-i: Code for Problem 1a
from fenics import import matplotlibmatplotlib . use (Agg)import matplotlib . pyplot as plt from ufl import nabla div import math import numpy as np from scipy . signal import argrelextrema from scipy . signal import savgol filter#============================================================== # Dimensional parameters#============================================================== length = 1.0W = 1.0 H = 0.001nx = 20ny = 20nz = 6youngs = 200e9 # Youngs nu = 0.3 # Poisson rho = 7800 # Density# Lame parameters mu = (youngs)/(2(1+nu) )lambda = (nuyoungs)/((1+nu)(12nu) ) F = 100; traction applied = F/((1 .98) (.51 .49) )#============================================================== # Dimensionless parameters#============================================================== youngs = (mu(3.0 lambda +2.0mu) ) /(lambda +mu) bar speed = math. sqrt (youngs/rho)lnd = length/lengthwnd = W/length hnd = H/lengtht char = length/bar speed t = 0 ti = 0.5dt = 0.1 num steps = 2000mu nd = mu/youngs lambda nd = lambda /youngs traction nd = traction applied /youngs#============================================================ mesh = BoxMesh( Point (0 ,0 ,0) , Point ( l nd , w nd , h nd) ,nx ,ny , nz) V = VectorFunctionSpace (mesh , P ,1)boundary left = near (x [0] ,0) bc left = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundary left ) tol = 1E14#============================================================ def epsilon (u) : return 0.5( nabla grad (u) + nabla grad (u) .T)def sigma(u) :return lambda ndnabla 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 , 0.0 , (x [0] >= 0.98 l && x [0] <= l )&& (x [1] >= 0.49w && x [1] <= 0.51w) && near (x [2] , h)? A : 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)dsa init , L init = lhs ( F init ) , rhs ( F init )print( First solving the i n i t i a l cantelever problem)u init = Function (V) solve ( a init==L init , u init , bc left )#============================================================# 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 = Expression (( near (x [0] , l ) ? ( t <= t i ? A : 0.0) :0.0 , 0.0 , 0.0 ) , degree=1, l=l nd , A=traction nd , t=t , t i= t i )T n = Constant ((0.0 ,0.0 ,0.0) )u = TrialFunction (V) d = u. geometric dimension () v = TestFunction (V)F = (dtdt)inner (sigma(u) , epsilon (v) )dx + dot(u, v)dx (dtdt) dot( f , v)dx (dtdt)dot (T n , v)ds 2.0 dot(u n , v)dx + dot( u n 1 , v)dxa ,L = lhs (F) , rhs (F)# xdmffileu = XDMFFile( results / solution . xdmf ) # xdmffile s = XDMFFile( results / stress . xdmf )u = Function (V)Q = TensorFunctionSpace (mesh , Lagrange , 1) stress proj = Function (Q)index = 0 displacements0 = [0] numsteps ; displacements1 = [0] numsteps ; displacements2 = [0] numsteps ; displacements3 = [0] numsteps ; displacements4 = [0] numsteps ; widths = [0 ,0.25 ,0.5 ,0.75 ,1] u store = [0] num steps time = [0] num stepsfor n in range( num steps ) : print(time = %.2f % t ) T n . t = tu = Function (V)solve (a == L, u, bc left )ugrab = u( length , widths [0]W,H) u store [n] = u grab [2] displacements0 [n] = u store [n]ugrab = u( length , widths [1]W,H) u store [n] = u grab [2] displacements1 [n] = u store [n]ugrab = u( length , widths [2]W,H) u store [n] = u grab [2] displacements2 [n] = u store [n]ugrab = u( length , widths [3]W,H) u store [n] = u grab [2] displacements3 [n] = u store [n]ugrab = u( length , widths [4]W,H) u store [n] = u grab [2] displacements4 [n] = u store [n]# i f ( abs ( tindex ) <0.01) :# print ( Calculating stress . . . )# xdmffile u . write (ulength , t )# stress = lambda nabla div (u) Identity (d) + mu( epsilon (u) + epsilon (u) .T)# stress proj . vector () [ : ] = project ( stress ,Q) . vector ()# xdmffile s . write ( stress proj , t )# index += 1t+=dtu n 1 . assign (u n) u n . assign (u) time [n] = tplt . figure (1) plt . plot (time , displacements0 , label= (L,0 ,W) ) plt . plot (time , displacements1 , label= (L,0.25H,W) ) plt . plot (time , displacements2 , label= (L,0.5H,W) ) plt . plot (time , displacements3 , label= (L,0.75H,W) ) plt . plot (time , displacements4 , label= (L,H,W) ) plt . xlabel ( Time [ s ] ) plt . ylabel ( Vertical Deflection [m] ) plt . legend ( loc= best ) plt . savefig ( results a /1 a i disps . png , bbox inches= tight )nat freq = [0] len( widths )displacements0 = savgol filter ( displacements0 , 51 , 4) u np = np. array ( displacements0 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [0] = 2math. pi /perioddisplacements1 = savgol filter ( displacements1 , 51 , 4) u np = np. array ( displacements1 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [1] = 2math. pi /perioddisplacements2 = savgol filter ( displacements2 , 51 , 4) u np = np. array ( displacements2 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [2] = 2math. pi /perioddisplacements3 = savgol filter ( displacements3 , 51 , 4) u np = np. array ( displacements3 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [3] = 2math. pi /perioddisplacements4 = savgol filter ( displacements4 , 51 , 4) u np = np. array ( displacements4 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [4] = 2math. pi /periodplt . figure (2) plt . plot (widths , nat freq , bx ) plt . xlabel ( Plate Width [m] ) plt . ylabel ( Natural Frequency [ rad/s ] ) plt . savefig ( results a /1 a i natfreq . png )
Appendix 1a-ii: Code for Problem 1a
from fenics import import matplotlibmatplotlib . use (Agg)import matplotlib . pyplot as plt from ufl import nabla div import math import numpy as np from scipy . signal import argrelextrema from scipy . signal import savgol filter#============================================================== # Dimensional parameters#============================================================== length = 1.0W = 1.0 H = 0.001nx = 20ny = 20nz = 6youngs = 200e9 # Youngs nu = 0.3 # Poisson rho = 7800 # Density# Lame parameters mu = (youngs)/(2(1+nu) )lambda = (nuyoungs)/((1+nu)(12nu) ) F = 100; traction applied = F/((1 .98) (.51 .49) )#============================================================== # Dimensionless parameters#============================================================== youngs = (mu(3.0 lambda +2.0mu) ) /(lambda +mu) bar speed = math. sqrt (youngs/rho)lnd = length/lengthwnd = W/length hnd = H/lengtht char = length/bar speed t = 0 ti = 0.5dt = 0.1 num steps = 2000mu nd = mu/youngs lambda nd = lambda /youngs traction nd = traction applied /youngs#============================================================ mesh = BoxMesh( Point (0 ,0 ,0) , Point ( l nd , w nd , h nd) ,nx ,ny , nz) V = VectorFunctionSpace (mesh , P ,1)boundary left = near (x [0] ,0) bc left = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundary left ) tol = 1E14#============================================================ def epsilon (u) : return 0.5( nabla grad (u) + nabla grad (u) .T)def sigma(u) :return lambda ndnabla 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 , 0.0 , (x [0] >= 0.98 l && x [0] <= l )&& (x [1] >= 0.49w && x [1] <= 0.51w) && near (x [2] , h)? A : 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)dsa init , L init = lhs ( F init ) , rhs ( F init )print( First solving the i n i t i a l cantelever problem)u init = Function (V) solve ( a init==L init , u init , bc left )#============================================================# 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 = Expression (( near (x [0] , l ) ? ( t <= t i ? A : 0.0) :0.0 , 0.0 , 0.0 ) , degree=1, l=l nd , A=traction nd , t=t , t i= t i )T n = Constant ((0.0 ,0.0 ,0.0) )u = TrialFunction (V) d = u. geometric dimension () v = TestFunction (V)F = (dtdt)inner (sigma(u) , epsilon (v) )dx + dot(u, v)dx (dtdt) dot( f , v)dx (dtdt)dot (T n , v)ds 2.0 dot(u n , v)dx + dot( u n 1 , v)dxa ,L = lhs (F) , rhs (F)# xdmffileu = XDMFFile( results / solution . xdmf ) # xdmffile s = XDMFFile( results / stress . xdmf )u = Function (V)Q = TensorFunctionSpace (mesh , Lagrange , 1) stress proj = Function (Q)index = 0 displacements0 = [0] numsteps ; displacements1 = [0] numsteps ; displacements2 = [0] numsteps ; displacements3 = [0] numsteps ; displacements4 = [0] numsteps ; widths = [0 ,0.25 ,0.5 ,0.75 ,1] u store = [0] num steps time = [0] num stepsfor n in range( num steps ) : print(time = %.2f % t ) T n . t = tu = Function (V)solve (a == L, u, bc left )ugrab = u(0.5 length , widths [0]W,H) u store [n] = u grab [2] displacements0 [n] = u store [n]ugrab = u(0.5 length , widths [1]W,H) u store [n] = u grab [2] displacements1 [n] = u store [n]ugrab = u(0.5 length , widths [2]W,H) u store [n] = u grab [2] displacements2 [n] = u store [n]ugrab = u(0.5 length , widths [3]W,H) u store [n] = u grab [2] displacements3 [n] = u store [n]ugrab = u(0.5 length , widths [4]W,H) u store [n] = u grab [2] displacements4 [n] = u store [n]# i f ( abs ( tindex ) <0.01) :# print ( Calculating stress . . . )# xdmffile u . write (ulength , t )# stress = lambda nabla div (u) Identity (d) + mu( epsilon (u) + epsilon (u) .T)# stress proj . vector () [ : ] = project ( stress ,Q) . vector ()# xdmffile s . write ( stress proj , t )# index += 1t+=dtu n 1 . assign (u n) u n . assign (u) time [n] = tplt . figure (1) plt . plot (time , displacements0 , label= (0.5L,0 ,W) ) plt . plot (time , displacements1 , label= (0.5L,0.25H,W) ) plt . plot (time , displacements2 , label= (0.5L,0.5H,W) ) plt . plot (time , displacements3 , label= (0.5L,0.75H,W) ) plt . plot (time , displacements4 , label= (0.5L,H,W) ) plt . xlabel ( Time [ s ] ) plt . ylabel ( Vertical Deflection [m] ) plt . legend ( loc= best ) plt . savefig ( results a /1 a ii disps . png , bbox inches= tight )nat freq = [0] len( widths )displacements0 = savgol filter ( displacements0 , 51 , 4) u np = np. array ( displacements0 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [0] = 2math. pi /perioddisplacements1 = savgol filter ( displacements1 , 51 , 4) u np = np. array ( displacements1 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [1] = 2math. pi /perioddisplacements2 = savgol filter ( displacements2 , 51 , 4) u np = np. array ( displacements2 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [2] = 2math. pi /perioddisplacements3 = savgol filter ( displacements3 , 51 , 4) u np = np. array ( displacements3 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [3] = 2math. pi /perioddisplacements4 = savgol filter ( displacements4 , 51 , 4) u np = np. array ( displacements4 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [4] = 2math. pi /periodplt . figure (2) plt . plot (widths , nat freq , bx ) plt . xlabel ( Plate Width [m] ) plt . ylabel ( Natural Frequency [ rad/s ] ) plt . savefig ( results a /1 a ii natfreq . png )
Appendix 1b-i: Code for Problem 1b
from fenics import import matplotlibmatplotlib . use (Agg)import matplotlib . pyplot as plt from ufl import nabla div import math import numpy as np from scipy . signal import argrelextrema from scipy . signal import savgol filter#============================================================== # Dimensional parameters#============================================================== length = 2.0W = 2.0 H = 0.001nx = 20 ny = 20nz = 6youngs = 200e9 # Youngs nu = 0.3 # Poisson rho = 7800 # Density# Lame parameters mu = (youngs)/(2(1+nu) )lambda = (nuyoungs)/((1+nu)(12nu) ) F = 100; traction applied = F/((1 .98) (.51 .49) )#============================================================== # Dimensionless parameters#============================================================== youngs = (mu(3.0 lambda +2.0mu) ) /(lambda +mu) bar speed = math. sqrt (youngs/rho)lnd = length/lengthwnd = W/length hnd = H/lengtht char = length/bar speed t = 0 ti = 0.5dt = 0.1 num steps = 2000mu nd = mu/youngs lambda nd = lambda /youngs traction nd = traction applied /youngs#============================================================ mesh = BoxMesh( Point (0 ,0 ,0) , Point ( l nd , w nd , h nd) ,nx ,ny , nz) V = VectorFunctionSpace (mesh , P ,1)boundary left = near (x [0] ,0) bc left = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundary left ) tol = 1E14#============================================================ def epsilon (u) : return 0.5( nabla grad (u) + nabla grad (u) .T)def sigma(u) :return lambda ndnabla 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 , 0.0 , (x [0] >= 0.98 l && x [0] <= l )&& (x [1] >= 0.49w && x [1] <= 0.51w) && near (x [2] , h)? A : 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)dsa init , L init = lhs ( F init ) , rhs ( F init )print( First solving the i n i t i a l cantelever problem)u init = Function (V) solve ( a init==L init , u init , bc left )#============================================================# 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 = Expression (( near (x [0] , l ) ? ( t <= t i ? A : 0.0) :0.0 , 0.0 , 0.0 ) , degree=1, l=l nd , A=traction nd , t=t , t i= t i )T n = Constant ((0.0 ,0.0 ,0.0) )u = TrialFunction (V) d = u. geometric dimension () v = TestFunction (V)F = (dtdt)inner (sigma(u) , epsilon (v) )dx + dot(u, v)dx (dtdt) dot( f , v)dx (dtdt)dot (T n , v)ds 2.0 dot(u n , v)dx + dot( u n 1 , v)dxa ,L = lhs (F) , rhs (F)# xdmffileu = XDMFFile( results / solution . xdmf ) # xdmffile s = XDMFFile( results / stress . xdmf )u = Function (V)Q = TensorFunctionSpace (mesh , Lagrange , 1) stress proj = Function (Q)index = 0 displacements0 = [0] numsteps ; displacements1 = [0] numsteps ; displacements2 = [0] numsteps ; displacements3 = [0] numsteps ; displacements4 = [0] numsteps ; widths = [0 ,0.25 ,0.5 ,0.75 ,1] u store = [0] num steps time = [0] num stepsfor n in range( num steps ) : print(time = %.2f % t ) T n . t = tu = Function (V)solve (a == L, u, bc left )ugrab = u( l nd , widths [0] w nd , h nd) u store [n] = u grab [2] displacements0 [n] = u store [n]ugrab = u( l nd , widths [1] w nd , h nd) u store [n] = u grab [2] displacements1 [n] = u store [n]ugrab = u( l nd , widths [2] w nd , h nd) u store [n] = u grab [2] displacements2 [n] = u store [n]ugrab = u( l nd , widths [3] w nd , h nd) u store [n] = u grab [2] displacements3 [n] = u store [n]ugrab = u( l nd , widths [4] w nd , h nd) u store [n] = u grab [2] displacements4 [n] = u store [n]# i f ( abs ( tindex ) <0.01) :# print ( Calculating stress . . . )# xdmffile u . write (ulength , t )# stress = lambda nabla div (u) Identity (d) + mu( epsilon (u) + epsilon (u) .T)# stress proj . vector () [ : ] = project ( stress ,Q) . vector ()# xdmffile s . write ( stress proj , t )# index += 1t+=dtu n 1 . assign (u n) u n . assign (u) time [n] = tplt . figure (1) plt . plot (time , displacements0 , label= (L,0 ,W) ) plt . plot (time , displacements1 , label= (L,0.25H,W) ) plt . plot (time , displacements2 , label= (L,0.5H,W) ) plt . plot (time , displacements3 , label= (L,0.75H,W) ) plt . plot (time , displacements4 , label= (L,H,W) ) plt . xlabel ( Time [ s ] ) plt . ylabel ( Vertical Deflection [m] ) plt . legend ( loc= best ) plt . savefig ( results b /1 a i disps . png , bbox inches= tight )nat freq = [0] len( widths )displacements0 = savgol filter ( displacements0 , 51 , 4) u np = np. array ( displacements0 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [0] = 2math. pi /perioddisplacements1 = savgol filter ( displacements1 , 51 , 4) u np = np. array ( displacements1 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [1] = 2math. pi /perioddisplacements2 = savgol filter ( displacements2 , 51 , 4) u np = np. array ( displacements2 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [2] = 2math. pi /perioddisplacements3 = savgol filter ( displacements3 , 51 , 4) u np = np. array ( displacements3 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [3] = 2math. pi /perioddisplacements4 = savgol filter ( displacements4 , 51 , 4) u np = np. array ( displacements4 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [4] = 2math. pi /periodplt . figure (2) plt . plot (widths , nat freq , bx ) plt . xlabel ( Plate Width [m] ) plt . ylabel ( Natural Frequency [ rad/s ] ) plt . savefig ( results b /1 a i natfreq . png )
Appendix 1b-ii: Code for Problem 1b
from fenics import import matplotlibmatplotlib . use (Agg)import matplotlib . pyplot as plt from ufl import nabla div import math import numpy as np from scipy . signal import argrelextrema from scipy . signal import savgol filter#============================================================== # Dimensional parameters#============================================================== length = 2.0W = 2.0 H = 0.001nx = 20 ny = 20nz = 6youngs = 200e9 # Youngs nu = 0.3 # Poisson rho = 7800 # Density# Lame parameters mu = (youngs)/(2(1+nu) )lambda = (nuyoungs)/((1+nu)(12nu) ) F = 100; traction applied = F/((1 .98) (.51 .49) )#============================================================== # Dimensionless parameters#============================================================== youngs = (mu(3.0 lambda +2.0mu) ) /(lambda +mu) bar speed = math. sqrt (youngs/rho)lnd = length/lengthwnd = W/length hnd = H/lengtht char = length/bar speed t = 0 ti = 0.5dt = 0.1 num steps = 2000mu nd = mu/youngs lambda nd = lambda /youngs traction nd = traction applied /youngs#============================================================ mesh = BoxMesh( Point (0 ,0 ,0) , Point ( l nd , w nd , h nd) ,nx ,ny , nz) V = VectorFunctionSpace (mesh , P ,1)boundary left = near (x [0] ,0) bc left = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundary left ) tol = 1E14#============================================================ def epsilon (u) : return 0.5( nabla grad (u) + nabla grad (u) .T)def sigma(u) :return lambda ndnabla 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 , 0.0 , (x [0] >= 0.98 l && x [0] <= l )&& (x [1] >= 0.49w && x [1] <= 0.51w) && near (x [2] , h)? A : 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)dsa init , L init = lhs ( F init ) , rhs ( F init )print( First solving the i n i t i a l cantelever problem)u init = Function (V) solve ( a init==L init , u init , bc left )#============================================================# 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 = Expression (( near (x [0] , l ) ? ( t <= t i ? A : 0.0) :0.0 , 0.0 , 0.0 ) , degree=1, l=l nd , A=traction nd , t=t , t i= t i )T n = Constant ((0.0 ,0.0 ,0.0) )u = TrialFunction (V) d = u. geometric dimension () v = TestFunction (V)F = (dtdt)inner (sigma(u) , epsilon (v) )dx + dot(u, v)dx (dtdt) dot( f , v)dx (dtdt)dot (T n , v)ds 2.0 dot(u n , v)dx + dot( u n 1 , v)dxa ,L = lhs (F) , rhs (F)# xdmffileu = XDMFFile( results / solution . xdmf ) # xdmffile s = XDMFFile( results / stress . xdmf )u = Function (V)Q = TensorFunctionSpace (mesh , Lagrange , 1) stress proj = Function (Q)index = 0 displacements0 = [0] numsteps ; displacements1 = [0] numsteps ; displacements2 = [0] numsteps ; displacements3 = [0] numsteps ; displacements4 = [0] numsteps ; widths = [0 ,0.25 ,0.5 ,0.75 ,1] u store = [0] num steps time = [0] num stepsfor n in range( num steps ) : print(time = %.2f % t ) T n . t = tu = Function (V)solve (a == L, u, bc left )ugrab = u(0.5 l nd , widths [0] w nd , h nd) u store [n] = u grab [2] displacements0 [n] = u store [n]ugrab = u(0.5 l nd , widths [1] w nd , h nd) u store [n] = u grab [2] displacements1 [n] = u store [n]ugrab = u(0.5 l nd , widths [2] w nd , h nd) u store [n] = u grab [2] displacements2 [n] = u store [n]ugrab = u(0.5 l nd , widths [3] w nd , h nd) u store [n] = u grab [2] displacements3 [n] = u store [n]ugrab = u(0.5 l nd , widths [4] w nd , h nd) u store [n] = u grab [2] displacements4 [n] = u store [n]# i f ( abs ( tindex ) <0.01) :# print ( Calculating stress . . . )# xdmffile u . write (ulength , t )# stress = lambda nabla div (u) Identity (d) + mu( epsilon (u) + epsilon (u) .T)# stress proj . vector () [ : ] = project ( stress ,Q) . vector ()# xdmffile s . write ( stress proj , t )# index += 1t+=dtu n 1 . assign (u n) u n . assign (u) time [n] = tplt . figure (1) plt . plot (time , displacements0 , label= (0.5L,0 ,W) ) plt . plot (time , displacements1 , label= (0.5L,0.25H,W) ) plt . plot (time , displacements2 , label= (0.5L,0.5H,W) ) plt . plot (time , displacements3 , label= (0.5L,0.75H,W) ) plt . plot (time , displacements4 , label= (0.5L,H,W) ) plt . xlabel ( Time [ s ] ) plt . ylabel ( Vertical Deflection [m] ) plt . legend ( loc= best ) plt . savefig ( results b /1 a ii disps . png , bbox inches= tight )nat freq = [0] len( widths )displacements0 = savgol filter ( displacements0 , 51 , 4) u np = np. array ( displacements0 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [0] = 2math. pi /perioddisplacements1 = savgol filter ( displacements1 , 51 , 4) u np = np. array ( displacements1 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [1] = 2math. pi /perioddisplacements2 = savgol filter ( displacements2 , 51 , 4) u np = np. array ( displacements2 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [2] = 2math. pi /perioddisplacements3 = savgol filter ( displacements3 , 51 , 4) u np = np. array ( displacements3 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [3] = 2math. pi /perioddisplacements4 = savgol filter ( displacements4 , 51 , 4) u np = np. array ( displacements4 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [4] = 2math. pi /periodplt . figure (2) plt . plot (widths , nat freq , bx ) plt . xlabel ( Plate Width [m] ) plt . ylabel ( Natural Frequency [ rad/s ] ) plt . savefig ( results b /1 a ii natfreq . png )
Appendix 1c-i: Code for Problem 1c
from fenics import import matplotlibmatplotlib . use (Agg)import matplotlib . pyplot as plt from ufl import nabla div import math import numpy as np from scipy . signal import argrelextrema from scipy . signal import savgol filter#============================================================== # Dimensional parameters#============================================================== length = 1.0W = 1.0 H = 0.001nx = 20 ny = 20nz = 6youngs = 200e9 # Youngs nu = 0.3 # Poisson rho = 7800 # Density# Lame parameters mu = (youngs)/(2(1+nu) )lambda = (nuyoungs)/((1+nu)(12nu) ) F = 100; traction applied = F/((1 .98) (.51 .49) )#============================================================== # Dimensionless parameters#============================================================== youngs = (mu(3.0 lambda +2.0mu) ) /(lambda +mu) bar speed = math. sqrt (youngs/rho)lnd = length/lengthwnd = W/length hnd = H/lengtht char = length/bar speed t = 0 ti = 0.5dt = 0.1 num steps = 2000mu nd = mu/youngs lambda nd = lambda /youngs traction nd = traction applied /youngs#============================================================ mesh = BoxMesh( Point (0 ,0 ,0) , Point ( l nd , w nd , h nd) ,nx ,ny , nz) V = VectorFunctionSpace (mesh , P ,1)# boundary left = near (x [0] ,0) # b c l e f t = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundary left )# boundary right = near (x [0] , l nd ) # bc right = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundary right ) 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 ) tol = 1E14#============================================================ def epsilon (u) : return 0.5( nabla grad (u) + nabla grad (u) .T)def sigma(u) :return lambda ndnabla 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 , 0.0 , (x [0] >= 0.98 l && x [0] <= l )&& (x [1] >= 0.49w && x [1] <= 0.51w) && near (x [2] , h)? A : 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( First 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 ])#============================================================# 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 = Expression (( near (x [0] , l ) ? ( t <= t i ? A : 0.0) :0.0 , 0.0 , 0.0 ) , degree=1, l=l nd , A=traction nd , t=t , t i= t i )T n = Constant ((0.0 ,0.0 ,0.0) )u = TrialFunction (V) d = u. geometric dimension () v = TestFunction (V)F = (dtdt)inner (sigma(u) , epsilon (v) )dx + dot(u, v)dx (dtdt) dot( f , v)dx (dtdt)dot (T n , v)ds 2.0 dot(u n , v)dx + dot( u n 1 , v)dxa ,L = lhs (F) , rhs (F)# xdmffileu = XDMFFile( results / solution . xdmf ) # xdmffile s = XDMFFile( results / stress . xdmf )u = Function (V)Q = TensorFunctionSpace (mesh , Lagrange , 1)stress proj = Function (Q)index = 0 displacements0 = [0] numsteps ; displacements1 = [0] numsteps ; displacements2 = [0] numsteps ; displacements3 = [0] numsteps ; displacements4 = [0] numsteps ; widths = [0 ,0.25 ,0.5 ,0.75 ,1] u store = [0] num steps time = [0] num stepsfor n in range( num steps ) : print(time = %.2f % t ) T n . t = tu = Function (V)solve (a == L, u, [ bc left , bc right ])ugrab = u( length , widths [0]W,H) u store [n] = u grab [2] displacements0 [n] = u store [n]ugrab = u( length , widths [1]W,H) u store [n] = u grab [2] displacements1 [n] = u store [n]ugrab = u( length , widths [2]W,H) u store [n] = u grab [2] displacements2 [n] = u store [n]ugrab = u( length , widths [3]W,H) u store [n] = u grab [2] displacements3 [n] = u store [n]ugrab = u( length , widths [4]W,H) u store [n] = u grab [2] displacements4 [n] = u store [n]# i f ( abs ( tindex ) <0.01) :# print ( Calculating stress . . . )# xdmffile u . write (ulength , t )# stress = lambda nabla div (u) Identity (d) + mu( epsilon (u) + epsilon (u) .T)# stress proj . vector () [ : ] = project ( stress ,Q) . vector ()# xdmffile s . write ( stress proj , t )# index += 1t+=dtu n 1 . assign (u n) u n . assign (u) time [n] = tplt . figure (1) plt . plot (time , displacements0 , label= (L,0 ,W) ) plt . plot (time , displacements1 , label= (L,0.25H,W) ) plt . plot (time , displacements2 , label= (L,0.5H,W) ) plt . plot (time , displacements3 , label= (L,0.75H,W) ) plt . plot (time , displacements4 , label= (L,H,W) ) plt . xlabel ( Time [ s ] ) plt . ylabel ( Vertical Deflection [m] ) plt . legend ( loc= best ) plt . savefig ( results c /1 c i disps . png , bbox inches= tight )nat freq = [0] len( widths )displacements0 = savgol filter ( displacements0 , 51 , 4) u np = np. array ( displacements0 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [0] = 2math. pi /perioddisplacements1 = savgol filter ( displacements1 , 51 , 4) u np = np. array ( displacements1 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [1] = 2math. pi /perioddisplacements2 = savgol filter ( displacements2 , 51 , 4) u np = np. array ( displacements2 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [2] = 2math. pi /perioddisplacements3 = savgol filter ( displacements3 , 51 , 4) u np = np. array ( displacements3 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [3] = 2math. pi /perioddisplacements4 = savgol filter ( displacements4 , 51 , 4) u np = np. array ( displacements4 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [4] = 2math. pi /periodplt . figure (2) plt . plot (widths , nat freq , bx ) plt . xlabel ( Plate Width [m] ) plt . ylabel ( Natural Frequency [ rad/s ] ) plt . savefig ( results c /1 c i natfreq . png )Appendix 1c-ii: Code for Problem 1c
from fenics import import matplotlibmatplotlib . use (Agg)import matplotlib . pyplot as plt from ufl import nabla div import math import numpy as np from scipy . signal import argrelextrema from scipy . signal import savgol filter#============================================================== # Dimensional parameters#============================================================== length = 1.0W = 1.0 H = 0.001nx = 20 ny = 20nz = 6youngs = 200e9 # Youngs nu = 0.3 # Poisson rho = 7800 # Density# Lame parameters mu = (youngs)/(2(1+nu) )lambda = (nuyoungs)/((1+nu)(12nu) ) F = 100; traction applied = F/((1 .98) (.51 .49) )#============================================================== # Dimensionless parameters#============================================================== youngs = (mu(3.0 lambda +2.0mu) ) /(lambda +mu) bar speed = math. sqrt (youngs/rho)lnd = length/lengthwnd = W/length hnd = H/lengtht char = length/bar speed t = 0 ti = 0.5dt = 0.1 num steps = 2000mu nd = mu/youngs lambda nd = lambda /youngs traction nd = traction applied /youngs#============================================================ mesh = BoxMesh( Point (0 ,0 ,0) , Point ( l nd , w nd , h nd) ,nx ,ny , nz) V = VectorFunctionSpace (mesh , P ,1)# boundary left = near (x [0] ,0) # b c l e f t = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundary left )# boundary right = near (x [0] , l nd ) # bc right = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundary right ) 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 ) tol = 1E14#============================================================ def epsilon (u) : return 0.5( nabla grad (u) + nabla grad (u) .T)def sigma(u) :return lambda ndnabla 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 , 0.0 , (x [0] >= 0.98 l && x [0] <= l )&& (x [1] >= 0.49w && x [1] <= 0.51w) && near (x [2] , h)? A : 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( First 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 ])#============================================================# 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 = Expression (( near (x [0] , l ) ? ( t <= t i ? A : 0.0) :0.0 , 0.0 , 0.0 ) , degree=1, l=l nd , A=traction nd , t=t , t i= t i )T n = Constant ((0.0 ,0.0 ,0.0) )u = TrialFunction (V) d = u. geometric dimension () v = TestFunction (V)F = (dtdt)inner (sigma(u) , epsilon (v) )dx + dot(u, v)dx (dtdt) dot( f , v)dx (dtdt)dot (T n , v)ds 2.0 dot(u n , v)dx + dot( u n 1 , v)dxa ,L = lhs (F) , rhs (F)# xdmffileu = XDMFFile( results / solution . xdmf ) # xdmffile s = XDMFFile( results / stress . xdmf )u = Function (V)Q = TensorFunctionSpace (mesh , Lagrange , 1)stress proj = Function (Q)index = 0 displacements0 = [0] numsteps ; displacements1 = [0] numsteps ; displacements2 = [0] numsteps ; displacements3 = [0] numsteps ; displacements4 = [0] numsteps ; widths = [0 ,0.25 ,0.5 ,0.75 ,1] u store = [0] num steps time = [0] num stepsfor n in range( num steps ) : print(time = %.2f % t ) T n . t = tu = Function (V)solve (a == L, u, [ bc left , bc right ])ugrab = u(0.5 length , widths [0]W,H) u store [n] = u grab [2] displacements0 [n] = u store [n]ugrab = u(0.5 length , widths [1]W,H) u store [n] = u grab [2] displacements1 [n] = u store [n]ugrab = u(0.5 length , widths [2]W,H) u store [n] = u grab [2] displacements2 [n] = u store [n]ugrab = u(0.5 length , widths [3]W,H) u store [n] = u grab [2] displacements3 [n] = u store [n]ugrab = u(0.5 length , widths [4]W,H) u store [n] = u grab [2] displacements4 [n] = u store [n]# i f ( abs ( tindex ) <0.01) :# print ( Calculating stress . . . )# xdmffile u . write (ulength , t )# stress = lambda nabla div (u) Identity (d) + mu( epsilon (u) + epsilon (u) .T)# stress proj . vector () [ : ] = project ( stress ,Q) . vector ()# xdmffile s . write ( stress proj , t )# index += 1t+=dtu n 1 . assign (u n) u n . assign (u) time [n] = tplt . figure (1) plt . plot (time , displacements0 , label= (0.5L,0 ,W) ) plt . plot (time , displacements1 , label= (0.5L,0.25H,W) ) plt . plot (time , displacements2 , label= (0.5L,0.5H,W) ) plt . plot (time , displacements3 , label= (0.5L,0.75H,W) ) plt . plot (time , displacements4 , label= (0.5L,H,W) ) plt . xlabel ( Time [ s ] ) plt . ylabel ( Vertical Deflection [m] ) plt . legend ( loc= best ) plt . savefig ( results c /1 c ii disps . png , bbox inches= tight )nat freq = [0] len( widths )displacements0 = savgol filter ( displacements0 , 51 , 4) u np = np. array ( displacements0 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [0] = 2math. pi /perioddisplacements1 = savgol filter ( displacements1 , 51 , 4) u np = np. array ( displacements1 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [1] = 2math. pi /perioddisplacements2 = savgol filter ( displacements2 , 51 , 4) u np = np. array ( displacements2 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [2] = 2math. pi /perioddisplacements3 = savgol filter ( displacements3 , 51 , 4) u np = np. array ( displacements3 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [3] = 2math. pi /perioddisplacements4 = savgol filter ( displacements4 , 51 , 4) u np = np. array ( displacements4 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [4] = 2math. pi /periodplt . figure (2) plt . plot (widths , nat freq , bx ) plt . xlabel ( Plate Width [m] ) plt . ylabel ( Natural Frequency [ rad/s ] ) plt . savefig ( results c /1 c ii natfreq . png )Appendix 1d-i: Code for Problem 1d
from fenics import import matplotlibmatplotlib . use (Agg)import matplotlib . pyplot as plt from ufl import nabla div import math import numpy as np from scipy . signal import argrelextrema from scipy . signal import savgol filter#============================================================== # Dimensional parameters#============================================================== length = 1.0W = 1.0 H = 0.001nx = 20 ny = 20nz = 6youngs = 200e9 # Youngs nu = 0.3 # Poisson rho = 7800 # Density# Lame parameters mu = (youngs)/(2(1+nu) )lambda = (nuyoungs)/((1+nu)(12nu) ) F = 100; traction applied = F/((1 .98) (.51 .49) )#============================================================== # Dimensionless parameters#============================================================== youngs = (mu(3.0 lambda +2.0mu) ) /(lambda +mu) bar speed = math. sqrt (youngs/rho)lnd = length/lengthwnd = W/length hnd = H/lengtht char = length/bar speed t = 0 ti = 0.5dt = 0.1 num steps = 2000mu nd = mu/youngs lambda nd = lambda /youngs traction nd = traction applied /youngs#============================================================ mesh = BoxMesh( Point (0 ,0 ,0) , Point ( l nd , w nd , h nd) ,nx ,ny , nz) V = VectorFunctionSpace (mesh , P ,1)# boundary left = near (x [0] ,0) # b c l e f t = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundary left )# boundary right = near (x [0] , l nd ) # bc right = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundary right ) 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 )def boundary front (x , on boundary) : return on boundary and near (x [1] ,0 , tol )def boundary back (x , on boundary) : return on boundary and near (x [1] , w nd , tol )bcleft = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundaryleft ) bcright = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundary right ) bcfront = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundary front ) bcback = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundary back ) tol = 1E14#============================================================ def epsilon (u) :return 0.5( nabla grad (u) + nabla grad (u) .T)def sigma(u) :return lambda ndnabla 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 , 0.0 , (x [0] >= 0.98 l && x [0] <= l )&& (x [1] >= 0.49w && x [1] <= 0.51w) && near (x [2] , h)? A : 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( First 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 , bc front , bc back ])#============================================================# 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 = Expression (( near (x [0] , l ) ? ( t <= t i ? A : 0.0) :0.0 , 0.0 , 0.0 ) , degree=1, l=l nd , A=traction nd , t=t , t i= t i )T n = Constant ((0.0 ,0.0 ,0.0) )u = TrialFunction (V) d = u. geometric dimension () v = TestFunction (V)F = (dtdt)inner (sigma(u) , epsilon (v) )dx + dot(u, v)dx (dtdt) dot( f , v)dx (dtdt)dot (T n , v)ds 2.0 dot(u n , v)dx + dot(u n 1 , v)dxa ,L = lhs (F) , rhs (F)# xdmffileu = XDMFFile( results / solution . xdmf ) # xdmffile s = XDMFFile( results / stress . xdmf )u = Function (V)Q = TensorFunctionSpace (mesh , Lagrange , 1) stress proj = Function (Q)index = 0 displacements0 = [0] numsteps ; displacements1 = [0] numsteps ; displacements2 = [0] numsteps ; displacements3 = [0] numsteps ; displacements4 = [0] numsteps ; widths = [0 ,0.25 ,0.5 ,0.75 ,1] u store = [0] num steps time = [0] num stepsfor n in range( num steps ) : print(time = %.2f % t ) T n . t = tu = Function (V)solve (a == L, u, [ bc left , bc right , bc front , bc back ])ugrab = u( length , widths [0]W,H) u store [n] = u grab [2] displacements0 [n] = u store [n]ugrab = u( length , widths [1]W,H) u store [n] = u grab [2] displacements1 [n] = u store [n]ugrab = u( length , widths [2]W,H) u store [n] = u grab [2] displacements2 [n] = u store [n]ugrab = u( length , widths [3]W,H) u store [n] = u grab [2] displacements3 [n] = u store [n]ugrab = u( length , widths [4]W,H) u store [n] = u grab [2] displacements4 [n] = u store [n]# i f ( abs ( tindex ) <0.01) :# print ( Calculating stress . . . )# xdmffile u . write (ulength , t )# stress = lambda nabla div (u) Identity (d) + mu( epsilon (u) + epsilon (u) .T)# stress proj . vector () [ : ] = project ( stress ,Q) . vector ()# xdmffile s . write ( stress proj , t )# index += 1t+=dtu n 1 . assign (u n) u n . assign (u) time [n] = tplt . figure (1) plt . plot (time , displacements0 , label= (L,0 ,W) ) plt . plot (time , displacements1 , label= (L,0.25H,W) ) plt . plot (time , displacements2 , label= (L,0.5H,W) ) plt . plot (time , displacements3 , label= (L,0.75H,W) ) plt . plot (time , displacements4 , label= (L,H,W) ) plt . xlabel ( Time [ s ] ) plt . ylabel ( Vertical Deflection [m] ) plt . legend ( loc= best ) plt . savefig ( results d /1 d i disps . png , bbox inches= tight )widths = [0.25 ,0.5 ,0.75 ,1] nat freq = [0] len( widths )# displacements0 = s a v g o l f i l t e r ( displacements0 , 51 , 4)# u np = np . array ( displacements0 )# min args = argrelextrema (u np , np . greater )# period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] )# nat freq [0] = 2math . pi /perioddisplacements1 = savgol filter ( displacements1 , 51 , 4) u np = np. array ( displacements1 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [0] = 2math. pi /perioddisplacements2 = savgol filter ( displacements2 , 51 , 4) u np = np. array ( displacements2 ) min args = argrelextrema (u np ,np. greater )period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [1] = 2math. pi /perioddisplacements3 = savgol filter ( displacements3 , 51 , 4) u np = np. array ( displacements3 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [2] = 2math. pi /perioddisplacements4 = savgol filter ( displacements4 , 51 , 4) u np = np. array ( displacements4 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [3] = 2math. pi /periodplt . figure (2) plt . plot (widths , nat freq , bx ) plt . xlabel ( Plate Width [m] ) plt . ylabel ( Natural Frequency [ rad/s ] ) plt . savefig ( results d /1 d i natfreq . png )Appendix 1d-ii: Code for Problem 1d
from fenics import import matplotlibmatplotlib . use (Agg)import matplotlib . pyplot as plt from ufl import nabla div import math import numpy as np from scipy . signal import argrelextrema from scipy . signal import savgol filter#============================================================== # Dimensional parameters#============================================================== length = 1.0W = 1.0 H = 0.001nx = 20 ny = 20nz = 6youngs = 200e9 # Youngs nu = 0.3 # Poisson rho = 7800 # Density# Lame parameters mu = (youngs)/(2(1+nu) )lambda = (nuyoungs)/((1+nu)(12nu) ) F = 100; traction applied = F/((1 .98) (.51 .49) )#============================================================== # Dimensionless parameters#============================================================== youngs = (mu(3.0 lambda +2.0mu) ) /(lambda +mu) bar speed = math. sqrt (youngs/rho)lnd = length/lengthwnd = W/length hnd = H/lengtht char = length/bar speed t = 0 ti = 0.5dt = 0.1 num steps = 2000mu nd = mu/youngs lambda nd = lambda /youngs traction nd = traction applied /youngs#============================================================ mesh = BoxMesh( Point (0 ,0 ,0) , Point ( l nd , w nd , h nd) ,nx ,ny , nz) V = VectorFunctionSpace (mesh , P ,1)# boundary left = near (x [0] ,0) # b c l e f t = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundary left )# boundary right = near (x [0] , l nd ) # bc right = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundary right ) 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 )def boundary front (x , on boundary) : return on boundary and near (x [1] ,0 , tol )def boundary back (x , on boundary) : return on boundary and near (x [1] , w nd , tol )bcleft = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundaryleft ) bcright = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundary right ) bcfront = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundary front ) bcback = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundary back ) tol = 1E14#============================================================ def epsilon (u) :return 0.5( nabla grad (u) + nabla grad (u) .T)def sigma(u) :return lambda ndnabla 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 , 0.0 , (x [0] >= 0.98 l && x [0] <= l )&& (x [1] >= 0.49w && x [1] <= 0.51w) && near (x [2] , h)? A : 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( First 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 , bc front , bc back ])#============================================================# 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 = Expression (( near (x [0] , l ) ? ( t <= t i ? A : 0.0) :0.0 , 0.0 , 0.0 ) , degree=1, l=l nd , A=traction nd , t=t , t i= t i )T n = Constant ((0.0 ,0.0 ,0.0) )u = TrialFunction (V) d = u. geometric dimension () v = TestFunction (V)F = (dtdt)inner (sigma(u) , epsilon (v) )dx + dot(u, v)dx (dtdt) dot( f , v)dx (dtdt)dot (T n , v)ds 2.0 dot(u n , v)dx + dot(u n 1 , v)dxa ,L = lhs (F) , rhs (F)# xdmffileu = XDMFFile( results / solution . xdmf ) # xdmffile s = XDMFFile( results / stress . xdmf )u = Function (V)Q = TensorFunctionSpace (mesh , Lagrange , 1) stress proj = Function (Q)index = 0 displacements0 = [0] numsteps ; displacements1 = [0] numsteps ; displacements2 = [0] numsteps ; displacements3 = [0] numsteps ; displacements4 = [0] numsteps ; widths = [0 ,0.25 ,0.5 ,0.75 ,1] u store = [0] num steps time = [0] num stepsfor n in range( num steps ) : print(time = %.2f % t ) T n . t = tu = Function (V)solve (a == L, u, [ bc left , bc right , bc front , bc back ])ugrab = u(0.5 length , widths [0]W,H) u store [n] = u grab [2] displacements0 [n] = u store [n]ugrab = u(0.5 length , widths [1]W,H) u store [n] = u grab [2] displacements1 [n] = u store [n]ugrab = u(0.5 length , widths [2]W,H) u store [n] = u grab [2] displacements2 [n] = u store [n]ugrab = u(0.5 length , widths [3]W,H) u store [n] = u grab [2] displacements3 [n] = u store [n]ugrab = u(0.5 length , widths [4]W,H) u store [n] = u grab [2] displacements4 [n] = u store [n]# i f ( abs ( tindex ) <0.01) :# print ( Calculating stress . . . )# xdmffile u . write (ulength , t )# stress = lambda nabla div (u) Identity (d) + mu( epsilon (u) + epsilon (u) .T)# stress proj . vector () [ : ] = project ( stress ,Q) . vector ()# xdmffile s . write ( stress proj , t )# index += 1t+=dtu n 1 . assign (u n) u n . assign (u) time [n] = tplt . figure (1) plt . plot (time , displacements0 , label= (0.5L,0 ,W) ) plt . plot (time , displacements1 , label= (0.5L,0.25H,W) ) plt . plot (time , displacements2 , label= (0.5L,0.5H,W) ) plt . plot (time , displacements3 , label= (0.5L,0.75H,W) ) plt . plot (time , displacements4 , label= (0.5L,H,W) ) plt . xlabel ( Time [ s ] ) plt . ylabel ( Vertical Deflection [m] ) plt . legend ( loc= best ) plt . savefig ( results d /1 d ii disps . png , bbox inches= tight )nat freq = [0] len( widths )displacements0 = savgol filter ( displacements0 , 51 , 4) u np = np. array ( displacements0 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [0] = 2math. pi /perioddisplacements1 = savgol filter ( displacements1 , 51 , 4) u np = np. array ( displacements1 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [1] = 2math. pi /perioddisplacements2 = savgol filter ( displacements2 , 51 , 4) u np = np. array ( displacements2 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [2] = 2math. pi /perioddisplacements3 = savgol filter ( displacements3 , 51 , 4) u np = np. array ( displacements3 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [3] = 2math. pi /perioddisplacements4 = savgol filter ( displacements4 , 51 , 4) u np = np. array ( displacements4 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [4] = 2math. pi /periodplt . figure (2) plt . plot (widths , nat freq , bx ) plt . xlabel ( Plate Width [m] ) plt . ylabel ( Natural Frequency [ rad/s ] ) plt . savefig ( results d /1 d ii natfreq . png )Appendix 1d2-i: Code for Problem 1d
from fenics import import matplotlibmatplotlib . use (Agg)import matplotlib . pyplot as plt from ufl import nabla div import math import numpy as np from scipy . signal import argrelextrema from scipy . signal import savgol filter#============================================================== # Dimensional parameters#============================================================== length = 2.0W = 2.0 H = 0.001nx = 20 ny = 20nz = 6youngs = 200e9 # Youngs nu = 0.3 # Poisson rho = 7800 # Density# Lame parameters mu = (youngs)/(2(1+nu) )lambda = (nuyoungs)/((1+nu)(12nu) ) F = 100; traction applied = F/((1 .98) (.51 .49) )#============================================================== # Dimensionless parameters#============================================================== youngs = (mu(3.0 lambda +2.0mu) ) /(lambda +mu) bar speed = math. sqrt (youngs/rho)lnd = length/lengthwnd = W/length hnd = H/lengtht char = length/bar speed t = 0 ti = 0.5dt = 0.1 num steps = 2000mu nd = mu/youngs lambda nd = lambda /youngs traction nd = traction applied /youngs#============================================================ mesh = BoxMesh( Point (0 ,0 ,0) , Point ( l nd , w nd , h nd) ,nx ,ny , nz) V = VectorFunctionSpace (mesh , P ,1)# boundary left = near (x [0] ,0) # b c l e f t = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundary left )# boundary right = near (x [0] , l nd ) # bc right = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundary right ) 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 )def boundary front (x , on boundary) : return on boundary and near (x [1] ,0 , tol )def boundary back (x , on boundary) : return on boundary and near (x [1] , w nd , tol )bcleft = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundaryleft ) bcright = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundary right ) bcfront = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundary front ) bcback = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundary back ) tol = 1E14#============================================================ def epsilon (u) :return 0.5( nabla grad (u) + nabla grad (u) .T)def sigma(u) :return lambda ndnabla 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 , 0.0 , (x [0] >= 0.98 l && x [0] <= l )&& (x [1] >= 0.49w && x [1] <= 0.51w) && near (x [2] , h)? A : 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( First 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 , bc front , bc back ])#============================================================# 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 = Expression (( near (x [0] , l ) ? ( t <= t i ? A : 0.0) :0.0 , 0.0 , 0.0 ) , degree=1, l=l nd , A=traction nd , t=t , t i= t i )T n = Constant ((0.0 ,0.0 ,0.0) )u = TrialFunction (V) d = u. geometric dimension () v = TestFunction (V)F = (dtdt)inner (sigma(u) , epsilon (v) )dx + dot(u, v)dx (dtdt) dot( f , v)dx (dtdt)dot (T n , v)ds 2.0 dot(u n , v)dx + dot(u n 1 , v)dxa ,L = lhs (F) , rhs (F)# xdmffileu = XDMFFile( results / solution . xdmf ) # xdmffile s = XDMFFile( results / stress . xdmf )u = Function (V)Q = TensorFunctionSpace (mesh , Lagrange , 1) stress proj = Function (Q)index = 0 displacements0 = [0] numsteps ; displacements1 = [0] numsteps ; displacements2 = [0] numsteps ; displacements3 = [0] numsteps ; displacements4 = [0] numsteps ; widths = [0 ,0.25 ,0.5 ,0.75 ,1] u store = [0] num steps time = [0] num stepsfor n in range( num steps ) : print(time = %.2f % t ) T n . t = tu = Function (V)solve (a == L, u, [ bc left , bc right , bc front , bc back ])ugrab = u( l nd , widths [0] w nd , h nd) u store [n] = u grab [2] displacements0 [n] = u store [n]ugrab = u( l nd , widths [1] w nd , h nd) u store [n] = u grab [2] displacements1 [n] = u store [n]ugrab = u( l nd , widths [2] w nd , h nd) u store [n] = u grab [2] displacements2 [n] = u store [n]ugrab = u( l nd , widths [3] w nd , h nd) u store [n] = u grab [2] displacements3 [n] = u store [n]ugrab = u( l nd , widths [4] w nd , h nd) u store [n] = u grab [2] displacements4 [n] = u store [n]# i f ( abs ( tindex ) <0.01) :# print ( Calculating stress . . . )# xdmffile u . write (ulength , t )# stress = lambda nabla div (u) Identity (d) + mu( epsilon (u) + epsilon (u) .T)# stress proj . vector () [ : ] = project ( stress ,Q) . vector ()# xdmffile s . write ( stress proj , t )# index += 1t+=dtu n 1 . assign (u n) u n . assign (u) time [n] = tplt . figure (1) plt . plot (time , displacements0 , label= (L,0 ,W) ) plt . plot (time , displacements1 , label= (L,0.25H,W) ) plt . plot (time , displacements2 , label= (L,0.5H,W) ) plt . plot (time , displacements3 , label= (L,0.75H,W) ) plt . plot (time , displacements4 , label= (L,H,W) ) plt . xlabel ( Time [ s ] ) plt . ylabel ( Vertical Deflection [m] ) plt . legend ( loc= best ) plt . savefig ( results d2 /1 d2 i disps . png , bbox inches= tight )widths = [0.25 ,0.5 ,0.75 ,1] nat freq = [0] len( widths )# displacements0 = s a v g o l f i l t e r ( displacements0 , 51 , 4)# u np = np . array ( displacements0 )# min args = argrelextrema (u np , np . greater )# period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] )# nat freq [0] = 2math . pi /perioddisplacements1 = savgol filter ( displacements1 , 51 , 4) u np = np. array ( displacements1 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [0] = 2math. pi /perioddisplacements2 = savgol filter ( displacements2 , 51 , 4) u np = np. array ( displacements2 ) min args = argrelextrema (u np ,np. greater )period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [1] = 2math. pi /perioddisplacements3 = savgol filter ( displacements3 , 51 , 4) u np = np. array ( displacements3 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [2] = 2math. pi /perioddisplacements4 = savgol filter ( displacements4 , 51 , 4) u np = np. array ( displacements4 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [3] = 2math. pi /periodplt . figure (2) plt . plot (widths , nat freq , bx ) plt . xlabel ( Plate Width [m] ) plt . ylabel ( Natural Frequency [ rad/s ] ) plt . savefig ( results d2 /1 d2 i natfreq . png )Appendix 1d2-ii: Code for Problem 1d
from fenics import import matplotlibmatplotlib . use (Agg)import matplotlib . pyplot as plt from ufl import nabla div import math import numpy as np from scipy . signal import argrelextrema from scipy . signal import savgol filter#============================================================== # Dimensional parameters#============================================================== length = 2.0W = 2.0 H = 0.001nx = 20 ny = 20nz = 6youngs = 200e9 # Youngs nu = 0.3 # Poisson rho = 7800 # Density# Lame parameters mu = (youngs)/(2(1+nu) )lambda = (nuyoungs)/((1+nu)(12nu) ) F = 100; traction applied = F/((1 .98) (.51 .49) )#============================================================== # Dimensionless parameters#============================================================== youngs = (mu(3.0 lambda +2.0mu) ) /(lambda +mu) bar speed = math. sqrt (youngs/rho)lnd = length/lengthwnd = W/length hnd = H/lengtht char = length/bar speed t = 0 ti = 0.5dt = 0.1 num steps = 2000mu nd = mu/youngs lambda nd = lambda /youngs traction nd = traction applied /youngs#============================================================ mesh = BoxMesh( Point (0 ,0 ,0) , Point ( l nd , w nd , h nd) ,nx ,ny , nz) V = VectorFunctionSpace (mesh , P ,1)# boundary left = near (x [0] ,0) # b c l e f t = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundary left )# boundary right = near (x [0] , l nd ) # bc right = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundary right ) 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 )def boundary front (x , on boundary) : return on boundary and near (x [1] ,0 , tol )def boundary back (x , on boundary) : return on boundary and near (x [1] , w nd , tol )bcleft = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundaryleft ) bcright = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundary right ) bcfront = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundary front ) bcback = DirichletBC (V, Constant ((0 ,0 ,0) ) , boundary back ) tol = 1E14#============================================================ def epsilon (u) :return 0.5( nabla grad (u) + nabla grad (u) .T)def sigma(u) :return lambda ndnabla 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 , 0.0 , (x [0] >= 0.98 l && x [0] <= l )&& (x [1] >= 0.49w && x [1] <= 0.51w) && near (x [2] , h)? A : 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( First 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 , bc front , bc back ])#============================================================# 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 = Expression (( near (x [0] , l ) ? ( t <= t i ? A : 0.0) :0.0 , 0.0 , 0.0 ) , degree=1, l=l nd , A=traction nd , t=t , t i= t i )T n = Constant ((0.0 ,0.0 ,0.0) )u = TrialFunction (V) d = u. geometric dimension () v = TestFunction (V)F = (dtdt)inner (sigma(u) , epsilon (v) )dx + dot(u, v)dx (dtdt) dot( f , v)dx (dtdt)dot (T n , v)ds 2.0 dot(u n , v)dx + dot(u n 1 , v)dxa ,L = lhs (F) , rhs (F)# xdmffileu = XDMFFile( results / solution . xdmf ) # xdmffile s = XDMFFile( results / stress . xdmf )u = Function (V)Q = TensorFunctionSpace (mesh , Lagrange , 1) stress proj = Function (Q)index = 0 displacements0 = [0] numsteps ; displacements1 = [0] numsteps ; displacements2 = [0] numsteps ; displacements3 = [0] numsteps ; displacements4 = [0] numsteps ; widths = [0 ,0.25 ,0.5 ,0.75 ,1] u store = [0] num steps time = [0] num stepsfor n in range( num steps ) : print(time = %.2f % t ) T n . t = tu = Function (V)solve (a == L, u, [ bc left , bc right , bc front , bc back ])ugrab = u(0.5 l nd , widths [0] w nd , h nd) u store [n] = u grab [2] displacements0 [n] = u store [n]ugrab = u(0.5 l nd , widths [1] w nd , h nd) u store [n] = u grab [2] displacements1 [n] = u store [n]ugrab = u(0.5 l nd , widths [2] w nd , h nd) u store [n] = u grab [2] displacements2 [n] = u store [n]ugrab = u(0.5 l nd , widths [3] w nd , h nd) u store [n] = u grab [2] displacements3 [n] = u store [n]ugrab = u(0.5 l nd , widths [4] w nd , h nd) u store [n] = u grab [2] displacements4 [n] = u store [n]# i f ( abs ( tindex ) <0.01) :# print ( Calculating stress . . . )# xdmffile u . write (ulength , t )# stress = lambda nabla div (u) Identity (d) + mu( epsilon (u) + epsilon (u) .T)# stress proj . vector () [ : ] = project ( stress ,Q) . vector ()# xdmffile s . write ( stress proj , t )# index += 1t+=dtu n 1 . assign (u n) u n . assign (u) time [n] = tplt . figure (1) plt . plot (time , displacements0 , label= (0.5L,0 ,W) ) plt . plot (time , displacements1 , label= (0.5L,0.25H,W) ) plt . plot (time , displacements2 , label= (0.5L,0.5H,W) ) plt . plot (time , displacements3 , label= (0.5L,0.75H,W) ) plt . plot (time , displacements4 , label= (0.5L,H,W) ) plt . xlabel ( Time [ s ] ) plt . ylabel ( Vertical Deflection [m] ) plt . legend ( loc= best ) plt . savefig ( results d2 /1 d2 ii disps . png , bbox inches= tight )nat freq = [0] len( widths )displacements0 = savgol filter ( displacements0 , 51 , 4) u np = np. array ( displacements0 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [0] = 2math. pi /perioddisplacements1 = savgol filter ( displacements1 , 51 , 4) u np = np. array ( displacements1 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [1] = 2math. pi /perioddisplacements2 = savgol filter ( displacements2 , 51 , 4) u np = np. array ( displacements2 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [2] = 2math. pi /perioddisplacements3 = savgol filter ( displacements3 , 51 , 4) u np = np. array ( displacements3 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [3] = 2math. pi /perioddisplacements4 = savgol filter ( displacements4 , 51 , 4) u np = np. array ( displacements4 ) min args = argrelextrema (u np ,np. greater ) period = ( time [ min args [ 0 ] [ 1 ] ] time [ min args [ 0 ] [ 0 ] ] ) nat freq [4] = 2math. pi /periodplt . figure (2) plt . plot (widths , nat freq , bx ) plt . xlabel ( Plate Width [m] ) plt . ylabel ( Natural Frequency [ rad/s ] ) plt . savefig ( results d2 /1 d2 ii natfreq . png )
Reviews
There are no reviews yet.