Analytical Beam Example¶
In this example, we’ll go through the process of setting up a model of an Euler-Bernoulli beam using the pywfe package.
Introduction¶
An Euler-Bernoulli beam can be described with a finite element approximation giving the mass and stiffness matrices:
For a beam segment of length \(l\), cross-sectional area \(A\) made from a material with Young’s modulus and density \(E, \rho\), and second moment of area \(I\). These matrices relate the displacement/rotation vector \([w_1, \theta_1, w_2, \theta_2]^T\) with the force/moment vector \([F_1, M_1, f_2, F_2]^T\) by
The FE model only has two nodes with two degrees of freedom each. The analytical formulation of an infinite beam has well known solutions. The dispersion relation for transverse waves is
The transfer mobility is subject to a transverse point force at \(x=0\) is
Creating pywfe Model of Beam¶
To begin with we define the system parameters
import numpy as np
import pywfe
import matplotlib.pyplot as plt
E = 2.1e11 # young mod
rho = 7850 # density
h = 0.1 # bean cross section side length length
A = h**2 # beam cross sectional area
I = h**4 / 12 # second moment of area
a = np.sqrt(E*I/(rho*A)) # factor in dispersion relation
and define the known solutions for the analytical dispersion relation and transfer mobility
def euler_wavenumber(f):
# wavenumber of euler bernoulli beam
return np.sqrt(2*np.pi*f/a)
def transfer_velocity(f, x):
# transfer velocity for beam x > 0
k = euler_wavenumber(f)
omega = 2*np.pi*f
return -omega/(4*E*I*k**3) * (1j*np.exp(-k*x) - np.exp(-1j*k*x))
For the FE discretisation, the beam length must be significantly shorter than the minimum wavelength. We define maximum frequency and find the maximum wavenumber analytically to set the beam length for WFE modelling.
f_max = 1e3 # maximum frequency
lambda_min = 2*np.pi/euler_wavenumber(f_max) # mimimum wavelength
l_max = lambda_min / 10 # unit cell length max - 10 unit cells per wavelength
l = np.round(l_max, decimals=1) # rounded unit cell length chosen
Now the mass and stiffness matrices can be defined
# stiffness matrix
K = E*I/(l**3) * np.array([
[12, 6*l, -12, 6*l],
[6*l, 4*l**2, -6*l, 2*l**2],
[-12, -6*l, 12, -6*l],
[6*l, 2*l**2, -6*l, 4*l**2]
])
# mass matrix
M = rho*A*l/420 * np.array([
[156, 22*l, 54, -13*l],
[22*l, 4*l**2, 13*l, -3*l**2],
[54, 13*l, 156, -22*l],
[-13*l, -3*l**2, -22*l, 4*l**2]
])
These, along with the ‘mesh’ information are all that are needed to create the pywfe.Model object. The mesh information is given with a dictionary with three keys node, fieldvar and coord.
These specify the node number, field variable, and coordinates in 1-3D of each degree of freedom in the model. The beam has 4 degrees of freedom, ordered as in the displacement vectors. Thus we define the dof dictionary
dof = {'node': [0, 0, 1, 1],
'fieldvar': ['w', 'phi']*2,
'coord': [
[0, 0, l, l],
[0, 0, 0, 0]
]
}
which describes the two nodes, the field quantities w, phi (repeated on each node), and the coordinates of each degree of freedom.
The coordinates are given in x and y with two lists for demonstrative purposes. Only the first is required for this 1D model.
NOTE: It is important that the waveguide axis is correctly chosen. By default the first axis in the coordinate set is used as the waveguide axis.
If your mesh information contains coordinates [x_arr, y_arr, z_arr] for example, and z_arr is the waveguide axis, then set axis = 2 or axis = -1 when initialising the model.
See ‘Initialisation’ in pywfe.Model.
The pywfe.Model object can now be created
beam_model = pywfe.Model(K, M, dof)
At this point, you might want to check the model with pywfe.Model.see(), which creates an interactive matplotlib view of the nodes in the mesh.
In this case however there is only one node to look at. The model can now be saved with pywfe.save().
Usage¶
Free Waves¶
Firstly let’s check the dispersion relation with the analytical solution
#create frequency array
f_arr = np.linspace(10, f_max, 100)
# calculate the wfe wavenumbers
k_wfe = beam_model.dispersion_relation(f_arr)
plt.plot(f_arr, euler_wavenumber(f_arr), '.', color='red', label='analytical')
plt.plot(f_arr, k_wfe, color='black')
plt.legend(loc='best')
plt.xlabel("frequency (Hz)")
plt.ylabel("wavenumber (1/m)")
Forcing¶
Forces can be added to degrees of freedom by changing elements of the pywfe.Model.force array. We compare the mobility in the WFE model with the known solution
beam_model.force[0] = 1
x_r = 0
w = beam_model.transfer_function(f_arr, x_r=x_r, dofs=[0], derivative=1)
plt.semilogy(f_arr, abs(transfer_velocity(f_arr, x_r)), '.', color='red', label='analytical')
plt.semilogy(f_arr, abs(w), color='black', label='WFE')
plt.legend(loc='best')
plt.xlabel("frequency (Hz)")
plt.ylabel("abs(mobility) (m/(Ns)")
The pywfe.Model.transfer_function() method calculates the response over all frequencies at the response distance x_r. The response distance can also be a list or array, in which case a higher dimensional array will be returned.
The dofs keyword argument specifies for which degrees of freedom the output should be returned. In this case we want the same dof as the one we’re forcing. The derivative keyword argument applies n derivatives in the
frequency domain, i.e a multiplication of the displacement by \(i \omega\). So the output of the method call is the transverse velocity at x=0 for a transverse unit point force. This is the mobility of the beam and is compared
with the analytical solution.
See pywfe.Model.transfer_function() for more information
To save the model for later, use pywfe.utils.io_utils.save(), where it may be saved in either in the current working directory or the database folder which is located in the user folder by default.
Before saving, you can add a description to the model to store metadata or other information with beam_model.description = '[description string]'.
More Functionality¶
For more functionality see pywfe.Model