11.2.2. Drums getting excited by a guitar
A band has a rest during rehearsal. As usual, the guitar guy cannot resist playing. By the incident sound wave, the drums are put into motion. We solve the drum excitation height \(h\) by a driven damped wave equation on a circular membrane, where we assume axisymmetric modes only, i.e.
where \(\omega=2\pi f\) is the frequency of the guitar sound incident as a planar wave with a forcing amplitude of \(F\). \(c\) depends on the drum and the surrounding gas and \(\delta\) is a damping coefficient. We demand that \(h(r=R)=0\) at the radius \(R\) of the drum.
As usual, we start with the drum equation in a coordinate-system-independent weak formulation:
from pyoomph import *
from pyoomph.expressions import *
from pyoomph.expressions.units import *
# Load tools for periodic driving response and text file output
from pyoomph.utils.periodic_driving_response import *
from pyoomph.utils.num_text_out import *
# Driven damped wave equation
class DrumEquation(Equations):
def __init__(self,c,damping,driving):
super().__init__()
self.c,self.damping,self.driving=c,damping,driving
def define_fields(self):
# Must be formulated first order here
self.define_scalar_field("h","C2",testscale=scale_factor("temporal")**2/scale_factor("h"))
self.define_scalar_field("hdot","C2",scale=scale_factor("h")/scale_factor("temporal"),testscale=scale_factor("temporal")/scale_factor("h"))
def define_residuals(self):
h,h_test=var_and_test("h")
hdot,hdot_test=var_and_test("hdot")
self.add_weak(partial_t(h)-hdot,hdot_test)
self.add_weak(partial_t(hdot)+self.damping*hdot-self.driving,h_test)
self.add_weak(self.c**2*grad(h),grad(h_test))
For the problem, we just define it on a line mesh in axisymmetric coordinates, which are indeed axisymmetric polar coordinates. Thereby, only axisymmetric modes are allowed. Alternatively, one could use a 2d circular Cartesian domain to also allow for nonaxisymmetric modes.
class DrumProblem(Problem):
def __init__(self):
super().__init__()
self.c=60*meter/second
self.damping=50/second
self.R_drum=30*centi*meter
self.driving=meter/second**2 *cos(2*pi*440*hertz*var("time"))
def define_problem(self):
self.set_coordinate_system("axisymmetric")
self+=LineMesh(N=100,size=self.R_drum,name="drum",left_name="center",right_name="rim")
self.set_scaling(h=1*centi*meter,temporal=0.001*second,spatial=self.R_drum)
eqs=DrumEquation(self.c,self.damping,self.driving)
eqs+=DirichletBC(h=0)@"rim"+AxisymmetryBC()@"center"
eqs+=TextFileOutput()
self+=eqs@"drum"
The general procedure is the same as before for a simple harmonic oscillator. However, here, we project the response on different Bessel modes:
with DrumProblem() as problem:
# Create the PeriodicDrivingResponse before the problem is initialized
pdr=PeriodicDrivingResponse(problem)
F=10*meter/second**2 # Driving amplitude, does not really matter, will cancel out in the normalized response
problem.driving=F*pdr.get_driving_mode() # means F*cos(omega*t)
# solve for a stationary state
problem.solve()
# Scan the frequency and write output
numbessel=10 # Number of bessel modes to project
bessel_roots=scipy.special.jn_zeros(0,numbessel)
outfile=NumericalTextOutputFile(problem.get_output_directory("response.txt"))
reson_freqs=[float(root*problem.c/problem.R_drum/hertz/(2*pi)) for root in bessel_roots]
print("RESONANT UNDAMPED FREQS",reson_freqs)
outfile.header("freq[Hz]",*["mode_"+str(i)+"[mm/(m/s^2)](f={:.4f})".format(reson_freqs[i]) for i in range(numbessel)])
# Output the resonant undamped frequencies
freqs=numpy.linspace(1,1000,1000) # Use frequencies f instead of omega
# Add interpolators to perform the Bessel projection. Mesh data comes in the mesh's own
# node order, which is not the order along the radius - and a spline with s=0 insists on a
# strictly increasing x - so sort by the coordinate first rather than trusting the order.
def radial_spline(data):
r=numpy.asarray(data.get_data("coordinate_x"))
order=numpy.argsort(r)
return scipy.interpolate.UnivariateSpline(r[order],numpy.asarray(data.get_data("h"))[order],k=3,s=0)
def project_and_write():
# Get the response as nondimensional data. The response is stored as eigenvector, so we split it in real and imaginary part.
# global_mesh=True because the projection integrates over the WHOLE radius: with --distribute
# each rank holds only its own arc of it, and a spline through a partial radius would give a
# different (and wrong) amplitude on every rank. These two calls are collective, i.e. every
# rank must reach them with the same arguments, and the merged result arrives on rank 0 only.
nd_resp_real=problem.get_cached_mesh_data("drum",eigenmode="real",eigenvector=0,nondimensional=True,global_mesh=True)
nd_resp_imag=problem.get_cached_mesh_data("drum",eigenmode="imag",eigenvector=0,nondimensional=True,global_mesh=True)
# From here on, everything is rank-0 work only: the other ranks got None above and
# NumericalTextOutputFile writes from rank 0 anyway. No collective call may follow.
if get_mpi_rank()!=0:
return
interr=radial_spline(nd_resp_real)
interi=radial_spline(nd_resp_imag)
# Calculate the Bessel decomposition of the response
bessel_data=[]
for i in range(numbessel):
bess_proj=lambda r : r*scipy.special.j0(bessel_roots[i]*r)*(interr(r)+1j*interi(r))
numer=scipy.integrate.quad(bess_proj,0,1,complex_func=True)[0] # Integrate Bessel projection
denom=1/2*scipy.special.j1(bessel_roots[i])**2 # Denominator of the Fourier-Bessel transform
nd_response_ampl=numpy.absolute(numer/denom)
dim_response_ampl_by_F=(nd_response_ampl*problem.get_scaling("h")/(milli*meter))/(F/(meter/second**2))
bessel_data.append(dim_response_ampl_by_F)
# add a row to the output
outfile.add_row(pdr.get_driving_frequency()/hertz,*bessel_data)
To obtain the full response data, we can access the eigenvector of the problem. Here, the response is stored in the eigenvector=0. The nondimensional eigenfunction is stored as real and imaginary contributions, which can be accessed by get_cached_mesh_data() for each mesh. Then, cubic interpolators are generated and in the loop, the individual Bessel modes are projected. The result is plotted in Fig. 11.2.
Fig. 11.2 Response of the drum separated into different Bessel modes that fulfill \(h(r=R)=0\). The dashed lines indicate the analytical eigenfrequencies of the undamped wave equation.