In the Acoustic Modes of a Rectangular Room episode we explored the analytical model of a rigid walled room with some Julia code. We focused on finding the resonance frequencies (or eigenfrequencies) of the room and calculating the related modal patterns (eigenfunctions). Now that, thanks to The Pulsating Sphere episode, we know how to setup Helmholtz problems with Elmer we can approach the problem with the FEM method. In this episode we will solve for the modal superposition in a rectangular rigid walled room and use the results from the Acoustic Modes of a Rectangular Room episode to check the accuracy.

Project Files

All the files used for this project are available at the repositories below:

• Source in the corner.
• Source in the middle.

Approach

As we were able to calculate the exact modal frequencies and patterns of a rectangular room ideally we would use a Helmholtz Eigensolver to find the resonance frequencies and their associated modal patterns through FEM. Unfortunately Elmer does not have a Helmholtz Eigensolver (at the time of writing), so we will have to approach the problem in a different way. We will:

1. Calculate the expected resonance frequencies with Julia.
2. Stimulate our room model at the various exact resonance frequencies by using a sound source.
3. Compare the results with the modal patterns.

Point 2 implies that we will need an active flux boundary condition on our domain, so we will need to place a source in our room. We will use a small pulsating sphere with uniform displacement in the corner of the room (spoiler: the placement of the source is important).

Point 3 assumes that the response of the room, when driven exactly at a modal frequency, is the related modal shape. In reality, as we seen in the Acoustic Modes of a Rectangular Room episode, this is not the case, as the response is instead a modal superposition. We will see where this assumption breaks down.

Whilst this approach is useful in showing part of the nature of acoustic fields in enclosed spaces, that is, the nature of modal superposition, it is a sideways and, as we will see, not very effective way of benchmark a solver accuracy as the benchmark solution (modal patterns) and solver solution (modal superposition) are inherently different. It is true that the Helmholtz solver does not support eigen analysis, but there are Elmer acoustic solvers that do. For a more direct approach, which makes use of eigen analysis, see the episodes below:

• Rigid Walled Room Revisited - Part 1
• Rigid Walled Room Revisited - Part 2
• Rigid Walled Room Revisited - Part 3

The Room

The parameters of the room are reported below:

Room Geometry

Medium Properties

Sphere Source Parameters

Building the FEM Model

Geometry

As always, let’s open FreeCAD and follow the procedure below:

• Select File from the top menu and then New.
• Switch to the Part Workbench.
• In the toolbar you should see a button with a yellow cube. Its tooltip says Create a cube solid. Click on it. This will create a cube.
• On the Combo View on the left, select the newly created cube. On the bottom you will see the cube properties. Input the following properties:
  • Under Box:
    • Length: 5000 mm
    • Width: 4000 mm
    • Height: 3000 mm
  • Under Base > Placement > Position (we do this to centre the room in the origin of the coordinate system):
    • x: -2500 mm
    • y: -2000 mm
    • z: -1500 mm
• In the toolbar you should see a button with a yellow sphere. Its tooltip says Create a sphere solid. Click on it. This will create a sphere.
• On the Combo View on the left, select the newly created sphere. On the bottom you will see the sphere properties.
  • This sphere will represent our source, so under Sphere > Radius we will keep the radius of 5 mm.
  • Under Base > Placement > Position (we do this to centre the sphere on one corner of the room):
    • x: 2500 mm
    • y: 2000 mm
    • z: 1500 mm
• Now, click in any empty space in the Combo View to un-select all objects. Then, while holding Ctrl, select the room first and the small sphere second. On the toolbar you will see a button showing two overlapping balls, one blue and the other one white. Its tooltip says Make a cut of two shapes. Click this button. The resulting object will be our room with the volume of our smaller sphere cut out of it.
• We are done. You can save the FreeCAD project, but do not forget to export our final object as BREP. To do so, select the final object in the Combo View on the left. Then, select File from the top men and then Export. On the Files of type combo, select BREP format and save your file to disk.

Note that, similarly to what we did in The Pulsating Sphere episode, we are only modelling the air enclosed within the room, so the source itself is omitted from the model. We will model the influence of the source on the air as a boundary condition.

Also note that our model is essentially a rectangular room with a small corner cut off that will act as the source. Why it is important to put the source in one corner? We will see below.

Meshing

We need to make some choices here. What we want to model? Let’s first have a look at our modal frequencies. We can generate a list of many modal frequencies as follows:

First, open Julia and install this package the AcousticModels.jl package (if you didn’t already). Then, type the following lines (copy and paste should work):

using AcousticModels

A, B, C = index_grid(100, 100, 100)

Lx = 5
Ly = 4
Lz = 3

F = mode_frequency.(A, B, C, Lx, Ly, Lz)

The 3D array F contains the modal frequencies found in a cube in the mode-numbers space. If we sort this array we get the various modal frequencies of our room in increasing order.

sort(F[:])
1030301-element Array{Float64,1}:
    0.0
   34.3
   42.87500000000001
   54.90679033598667
   57.16666666666666
   66.66721666439793
   68.6
   71.45833333333333
    ⋮
 7882.4189458425635
 7885.2708577448375
 7888.466115690857
 7896.957894517217
 7903.291666666668
 7911.618830415036
 7926.401076641138

We can try to understand what story these frequencies are telling us by plotting them as follows:

On the top panel we just reported the first $1000$ elements of the sorted F array. As we can see, the modal frequencies rise steeply, telling us that there are more and more modes per unity frequency the higher the mode number (in here by mode number we mean a number by which we order the frequencies in increasing order).

On the bottom panel instead we reported the distance, in cents, between adjacent modal frequencies, which in Julia we can compute with the following lines:

# First 1000 modal frequencies:
f = sort(F[:])[1:1000] 

# Distance in cents between adjacent frequencies:
c = 1200 .* log2.(f[2:end] ./ f[1:(end - 1)])

We can see that the modal frequencies get very close together as the mode number (and hence the frequency) rises, the difference being around $100$ $\text{cent}$ ($1$ semitone distance) already at the $10$-th mode. This is why modal frequencies are normally an issue for room acoustics only at low frequency. In fact, as we already noted in the Acoustic Modes of a Rectangular Room episode, at low frequency (low mode number) the resonance frequencies are far away from each other, so they produce strong filtering (in human perception) while the various resonances sort of merge together at high frequency, producing a more uniform perception.

Hence, we will focus only on the first $10$ nonzero modes, which we can get like this:

sort(F[:])[2:11]
10-element Array{Float64,1}:
 34.3             
 42.87500000000001
 54.90679033598667
 57.16666666666666
 66.66721666439793
 68.6             
 71.45833333333333
 79.26401076641136
 80.89638820738537
 85.75000000000001

As we learned in The Pulsating Sphere episode, our maximum mesh size $s$ should be lower than one tenth of the wavelength:

$$s<\frac{\lambda}{10}=\frac{c_{0}}{10f}$$

Here we have many frequencies, so we will use the highest frequency as a reference for our mesh calculation ($85.75$ $\text{Hz}$). If we choose a mesh size of $0.3$ $\text{m}$ we will be able to model up to $114.34$ $\text{Hz}$, which is safely higher than $85.75$ $\text{Hz}$. So, we will go with that.

Now that we know what to do, let’s open Salome and follow the procedure below:

• On the toolbar, you will see the Modules combo box. Select the Geometry module.
• Select File from the top menu, then Import and then BREP. This will open the Import BREP window. Select the BREP file you exported when following the Geometry section.
• In the Object Browser on the left you will see your newly imported object. Select it.
• Select New Entity from the top menu, and then Explode. This will open the Sub-shapes Selection window. On the Sub-shapes Type combo, select Solid. Then choose Apply and Close. This will create a new solid object under our original object tree.
• Select the newly created solid object from the Object Browser (expand the original object tree if necessary).
• Select New Entity from the top menu, and then Explode. This will open the Sub-shapes Selection window. On the Sub-shapes Type combo, select Face. Then choose Apply and Close. This will create new face objects under the solid object tree.
• On the toolbar, you will see the Modules combo box. Select the Mesh module.
• From the object browser, make sure to select the solid object that we created during our first explosion, and not the top root imported BREP object. Expand the tree if necessary.
• Once the correct object is selected, select Mesh from the top menu and then Create Mesh. This will open the Create mesh window. In the Algorithm combo, select NETGEN 1D-2D-3D Then click on the button with the gear wheel next to the Hypothesis combo. From the popup menu, select Netgen 3D Parameters.
  • Input the following parameters. Type 300 in the Max. size edit box. Then type 1 in the Min. size edit box. You should interpret these sizes as millimetres. Select Fine for the Fineness. Finally, tick the Second order tickbox and click OK.
• Once you done with the options above, click on Apply and Close in the Create mesh window.
• Select your newly created mesh on the Object Browser and right click with your mouse. A menu will popup. Select Create Groups from Geometry. This will open a new window. With this window open, select the solid and face entities that were exploded in the previous steps. Then click on Apply and Close.
• At this point we are ready to compute the mesh. Select and right click on the mesh on the Object Browser. A menu will popup. Select Compute and wait until it is done.
• After the computation is complete, select and right click on the mesh on the Object Browser. A menu will popup. Select Export and then UNV File. Save the UNV file wherever you want on your hard drive.

Solving

We are now ready to solve our problem. We will make use of MATC expressions again, as in The Pulsating Sphere episode, but we will also do a scanning simulation to solve for all of our frequencies one after the other. A scanning simulation works by creating a dummy time variable that takes discrete values ($1$, $2$, …) up until a maximum that we specify. Then, simulations will be repeated for each value of the dummy time. If we do something with this dummy time at each iteration, then we will have a nice parameter sweep. Here we will put our $10$ resonance frequencies into a vector and, at each iteration, use the dummy time to index into the vector and use the frequency value as input to the solver, so that at the end we have a collection of solutions, each one for each of our target resonance frequencies.

As a note, we mentioned above that we will be setting our source to be an uniform displacement source. This is in contrast with The Pulsating Sphere episode, where we used uniform velocity. The Elmer Model Manual tells us how to setup such a boundary conditions. The flux will be given by:

$$g=-\left(2\pi f\right)^{2}\rho_{0}d$$

Where $g$ is the flux and $d$ is the normal surface displacement (we do not have any other component of the displacement in our simulation). We will set $d$ to $1$ $\text{mm}$. Note that this flux is purely real.

Finally, to model a perfectly rigid wall, we will impose a null flux on the boundaries.

Follow the procedure below to setup our model.

Open up ElmerGUI by typing ElmerGUI in the terminal. Then, proceed as follows.

• Select File from the top menu and then Open. Navigate to the location of your UNV mesh to open it. This will load the mesh into the 3D view.
• Select Model from the top menu and then Setup.
  • In the Coordinate Scaling edit box, type 0.001, so that Elmer understands that the units of our UNV mesh are millimetres.
  • In the Simulation type combo, select Scanning.
  • In the Timestep intervals editbox, type 10.
  • In the Free text edit box type the code in Snippet 1 below. This code contains the simulation variables. To change anything about the simulation, it is sufficient to change any value here.
  • Click Apply.
• Select Model from the top menu, then Equation then Add. This will open the Equation window.
  • Select the Helmholtz Equation tab.
  • Click the tick beside Active.
  • At the bottom, under Apply to bodies, tick the body to which apply our equation (there should be only one).
  • In Free text input type the code in Snippet 2 below. This code sets the simulation frequency at each scanning step.
  • Click the Edit Solver Settings button.
    • In the Linear system tab, select the Iterative method and, from the combo, select BiCGStabl. Below, from the Preconditioning combo, select ILUT. In my experience, these parameters work well with Helmholtz simulations.
    • I suggest you also tick the box Abort if the solution did not converge, so that ugly solutions will not be saved if the solver cannot converge.
    • In the Nonlinear system tab, type 1 in the Max. iterations editbox (our problem is linear).
    • Click on Apply.
  • Click on OK.
• Select Model from the top menu, then Material then Add.
  • Click the Material library button.
    • Select Air (room temperature).
    • Click OK.
  • In the General tab, replace the value of Density with the MATC expression in Snippet 3.
  • At the bottom, under Apply to bodies, tick the body to which apply our equation (there should be only one).
  • Click OK.
• Select Model from the top menu, then Boundary condition then Add.
  • This will be our radiation condition, so I suggest you edit the Name edit box at the bottom and type Radiator.
  • Then, select the Helmholtz Equation tab.
  • Under Real part of the flux, insert the code in Snippet 4.
  • Click OK.
• Select Model from the top menu, then Boundary condition then Add.
  • This will be our rigid wall condition, so I suggest you edit the Name edit box at the bottom and type Rigid.
  • Then, select the Helmholtz Equation tab.
  • Under Real part of the flux, type 0.
  • Under Imag part of the flux, type 0.
  • Click OK.
• Select Model from the top menu, then Set boundary properties.
• Rotate the domain until you find the corner where you cut out the spherical source. Double click on it (see below for an example). A window will popup. From the Boundary condition combo, select Radiator.
• Double click on any wall. A window will popup. From the Boundary condition combo, select Rigid. Repeat this for every remaining wall.
• Select File from the top menu and then Save project. A window will popup.
  • In the new window, click on the Create Folder button. This folder is actually the project folder, and will contain all of our files. Give it a catchy name.
  • Once in the new folder click on the Open button. This will actually save all the project files in the new folder.

The last points about saving the project are a bit strange, so ensure you follow them.

To check that everything is alright, select Sif from the top menu and then Edit. This will open the Solver Input File generated by the GUI. If everything is OK, it should look like this.

If everything is OK, you are ready to press the Start solver button on the toolbar of the ElmerGUI program. This will start the solver and open the converge plot viewer and the log file viewer. Sit back while the simulation runs!

You will probably see a convergence plot like the one below, but that is not really a convergence plot as it reports results from unrelated solver runs, each one with different frequency. You can dismiss the plot.

Snippet 1

MATC definitions for the simulation. The code below defines a vector of $10$ elements f, which contains the first $10$ resonance frequencies of our room. For this reason, we set the Timestep intervals to 10 in the Setup window. p is the medium density, in kilograms per cubic meter, and d is the uniform source surface displacement, in meters.

$ f = 34.3 42.87500000000001 54.90679033598667 57.16666666666666 66.66721666439793 68.6 71.45833333333333 79.26401076641136 80.89638820738537 85.75000000000001
$ p  = 1.205
$ d = 0.001

Snippet 2

MATC expression to set the frequency for the Helmholtz solver. The scanning simulation works by advancing a variable, called time, by $1$ in each step. This variable starts at $1$, but MATC arrays are $0$ indexed. For this reason, we index into f by using tx - 1. So, in order to select the correct value of frequency for each step, we do the following:

Frequency = Variable time; Real MATC "f(tx - 1)"

Snippet 3

This code is used to assign the Real value resulting from the MATC expression between quotes to the required parameter. In this case, the medium density.

Real MATC "p"

Snippet 4

This code is used to assign the Real value resulting from the MATC expression between quotes to the required parameter. In this case, the real part of the Flux boundary condition. This is the only part of the Flux boundary condition for this problem.

Variable time; Real MATC "-((2 * pi * f(tx - 1))^2) * p * d""

Results

Our simulation is solved. But is it good? Let’s figure it out. Let’s head over again to the Elmer Models Manual to have a better look at what we calculated. The solver gave us a complex field defined over the whole of the room, called $P$. But what does a complex field mean? First of all, let’s remember that we actually solved for $10$ different fields, at $10$ different frequencies, so let’s add a suffix $n$ to $P$, so that we can identify them all: $P_{n}$, for $n$ going from $1$ to $10$. Then, remember that that the Helmholtz solver gives us the spatial part of the steady state solution, and the spatial part only depends on space variables: $P_{n} = P_{n}\left(x,y,z\right)$. So, If we want to include time, we need to multiply this spatial part by this complex sine wave:

$$\exp\left(j2\pi f_{n}t\right)$$

Where $j$ is the imaginary unit, $f_{n}$ are the resonance frequencies we used as input for the system, and $t$ is time, measured in seconds. Now we have all the ingredients to derive our steady state pressure waves $p_{n}\left(x,y,z,t\right)$. According to the Elmer Models Manual the physical solution is the real part of the complex field after multiplication with the complex sine wave:

$$p_{n}\left(x,y,z,t\right)=\Re\left[P_{n}\left(x,y,z\right)\exp\left(j2\pi f_{n}t\right)\right]$$

This is a very common convention, and people often say “we use complex numbers as they make algebra simpler, but the physical solution is the real part, so at the end we discard the imaginary part”. Well, this is not completely accurate. Let’s expand the equation above:

$$p_{n}\left(x,y,z,t\right)=\Re\left[P_{n}\left(x,y,z\right)\right]\cos\left(2\pi f_{n}t\right)-\Im\left[P_{n}\left(x,y,z\right)\right]\sin\left(2\pi f_{n}t\right)$$

Where $\Re$ and $\Im$ denote the real and imaginary parts respectively. As you can see, to build our final solution we use both real and imaginary parts of $P_{n}$. But also note this: if $t=\frac{k}{2f}$, with $k$ any integer number, then the sine term will be $0$, while the cosine term will be either $1$ or $-1$. This means the real part of $P_{n}$ is controlling the vibration at those times. The opposite happens instead at $t=\frac{2k+1}{4f}$, where instead the cosine is $0$ and the sine is either $1$ or $-1$. For all the other times $p_{n}$ is a sum of the real and imaginary parts of $P_{n}$, with weights decided by the trigonometric functions. So, the imaginary parts of things matter, if we discarded the imaginary part of $P_{n}$ prematurely we would have ended up with a solution which is good only for certain discrete times (for example, $t=0$). So, don’t think that complex description of wave is just mathematical convenience: for most of your quantities both the real and imaginary part carry physical meaning.

Cool, now let’s build the benchmark. In the Acoustic Modes of a Rectangular Room episode we seen that the modal shape of the room is given by:

$$ p\left(x,y,z,t\right)=S\left(x,y,z\right)T\left(t\right)$$

Where $S$ is the spatial part of the mode, and $T$ is a cosine wave. Remembering that we have $10$ of them, and to avoid confusion the FEM solution above, let’s denote the mode shape with an overline:

$$\overline{p}_{n}\left(x,y,z,t\right)=S_{n}\left(x,y,z\right)T_{n}\left(t\right)$$

Where $S$ and $T$ are given as follows (note that $S$ is different with respect the Acoustic Modes of a Rectangular Room episode to account for the fact that the coordinate system is, in this episode, centred in the middle of the room rather than on a corner):

$$S\left(x,y,z\right)=\cos\left(n_{x}\pi\left[\frac{x}{L_{x}}+\frac{1}{2}\right]\right)\cos\left(n_{y}\pi\left[\frac{y}{L_{y}}+\frac{1}{2}\right]\right)\cos\left(n_{z}\pi\left[\frac{z}{L_{z}}+\frac{1}{2}\right]\right)$$

$$T_{n}\left(t\right)=\cos\left(2\pi f_{n}t\right)$$

Where $n_x$, $n_y$ and $n_z$ are the mode numbers associated with $f_n$ (for more info refer back to the Acoustic Modes of a Rectangular Room episode.

Now, we should really not forget that the benchmark field we put together above is a single mode (for each $n$). But the solution from FEM is different. What we solved for with Elmer is not a single modal shape, but a field driven by a source instead. Turns out that for modal systems the steady state vibration, when driven by a source, is a weighted sum of all the possible modal shapes:

$$p_{\text{Expected}}\left(x,y,z,t\right)=\sum_{n=0}^{\infty}a_{n}\overline{p}_{n}\left(x,y,z,t\right)$$

Where $a_{n}$ are the coefficients of the superposition, and depend on where the source is located, the properties of the source and the boundary conditions. They are very hard to calculate analytically.

Still, two things are going on in this simulation:

• The walls are rigid.
• We are modelling low frequency.

The two things above, in combination, tell us that the modal superposition will be either pure (i.e. the only mode in the superposition is that related to the resonance frequency in input) or only “contaminated” by neighbouring modal shapes. This because rigid walls make resonances sharp, so the influence of each mode is narrow in frequency, making its ability to contaminate other modes reduced. But also, and more importantly, because at low frequency the distance between modes is larger, further a part than the resonance width, so modes do not really interact much in the superposition.

So, in other words: the response of the room is a modal superposition, but for rigid rooms at low frequency the superposition is going to be pure, or just a little bit dirty. But at what frequency this will stop? It is hard to tell, but we can find it by comparing the FEM result against the pure modal shapes as computed by the exact model. The lowest frequency modes should match the pure modes closely. If at some frequency we see mismatch, then we know that the superposition is not pure anymore.

Wait, but could other things cause mismatch, like solver accuracy? As we seen in The Pulsating Sphere episode we know that our mesh size should grant us low errors. Moreover, if the error patterns that we observe have a modal taste to them, and are not random or uniform, we know that the source of mismatch must be contamination from neighbouring modes.

So, to sum it up, Elmer cannot solve for the exact modal shapes (with the Helmholtz solver), but can solve for the steady state pressure as driven by a source. This pressure is given by a modal superposition, and by comparing it with the pure modes, given that our room is rigid and our frequencies are low, we will:

• Understand how Elmer is accurate by checking how close the lowest frequency modes are to the pure modal shapes.
• Understand at which frequency the modal superposition stops being pure, and more than one modes contribute to the steady state pressure.

Even though modal shapes and modal superposition are not the same thing, we will still refer to the comparison between the superposition delivered by Elmer and the pure modes delivered by Julia as error.

Only one thing is left to figure out. Since in our FEM model we have a driving source, the amplitude of the resulting field is dictated by the source. But in our benchmark solution (the pure mode) the amplitude is $1$. In order to compare the two, we will normalise the FEM solution as follows:

• Find the peak of the exact mode shape and the position at which this is reached.
• Multiply the FEM solution by a factor that brings the value at the same position to that of the exact solution.

So, let’s go ahead and let’s have a look at our FEM solution first with ParaView. Open ParaView, select File from the top menu and then Open. Navigate to the Elmer project folder. You should see the entry case..vtu. If you expand it you will see that it is just a collection of vtu files (case0001.vtu, case0002.vtu, …). These files are generated by Elmer for each time step of our scanning simulation. If we select case..vtu we will import them all into ParaView. Remember to click on Apply in the Properties browser on the left.

On the toolbar, you will see a combo that, by default, says Solid Color. Click on it and select pressure wave 1. This is the real part of the pressure field, while pressure wave 2 is the imaginary part. If you use the play/stop buttons at the top you can cycle through all 10 solutions. This will show you all the different pressure fields. The animation below shows all the fields, but normalised to take only values between $1$ and $-1$ so that they can be plotted on the same scale.

The animation was performed with ParaView by selecting only the real part of the field. As we seen above, we can think of this as being the result field at $t=0$. The various surfaces are surfaces of constant normalised pressure. They help us understanding the shape of the field. As you probably remember from the Acoustic Modes of a Rectangular Room episode, the shapes do really look similar to modes. You will probably also see some strange surfaces close to the corner (the source location) or the walls. Those exists due to the negative peak value of the field not being exactly $-1$.

However, there is something important to note. Depending on the Elmer version, and details of the meshing and boundary conditions, the results might out of phase. For example, in the video above, done with Elmer 8.4 (if my memory serves me well) the real part bears a high resemblance to the modal patterns. But with with Elmer 9.0 it will be the imaginary part to do so instead. In other words, there is an unknown phase factor $\theta_{?}$ which depends on the solver and boundary conditions, so our solution is more like this:

$$p_{n}\left(x,y,z,t\right)=\Re\left[P_{n}\left(x,y,z\right)\exp\left(j2\pi f_{n}t\right)\exp\left(\theta_{?}\right)\right]$$

This can make validation impervious, due to the value of $\theta_{?}$ that maximise the similarity of the field to a given modal shape (for a certain time $t$) being unknown. We will not try to determine this $\theta_{?}$ value here. Instead, the following will treat the lucky case (as it was in the solution used to write this article) in which the maximal similarity between solver solution and modal shapes is found at $\theta_{?}=0$ and $t=0$.

Now to validation. To validate, we should export the solution into CSV files, so that we can read it with Julia. In the Pipeline Browser, select the data you imported. Then, from the top menu choose File and Save Data. If you want to use the validate.jl script, save them in the export sub-directory of the repository, type data in the File name: editbox and choose csv from the File of type: combo. This will open the Configure Writer (CSVWriter) window. Tick Write Time Steps. Tick Choose Arrays to Write and untick GeometryIds. In Precision type 30. Click OK. The exported CSV files will contain the value of the fields at the nodal points of the mesh.

The validation is carried on by the Julia script validate.jl by using the exported CSV files. The script will use the equations above to calculate two metrics of error. One is the norm. This metric is calculated by putting the FEM solution and the exact one in two arrays, taking the difference, summing the square of the difference and taking the squared root of the result and diving over the room volume:

$$\text{Norm}_{n}=\frac{1}{V}\sqrt{\sum{|P_{n} - S_{n}|^2}}$$

Where $V$ is the room volume, $P_{n}$ was normalised as discussed above and the sum has to be though as performed over the fields $P_{n}$ and $S_{n}$ sampled at the mesh nodes. The norm has the units of $\frac{\text{Pa}}{\text{m}^3}$ and can be though of as a metric of average error. Remember that in this case by error we are meaning the degree of dissimilarity from the modal superposition $P_{n}$ and the pure modal shape $S_{n}$. The results are reported below for each of the study frequencies $f_{n}$:

It is possible to see that the error is overall rising with frequency. This is to be expected because the higher the frequency the higher the contamination from neighbouring modes that we should expect. On a similar note, error is also locally higher for those frequencies closest to other modal frequencies, such as $68.60$ $\text{Hz}$ and $54.91$ $\text{Hz}$.

The validation script also calculates the ratio in $\text{dB}$ between the FEM solution and the exact one is exported in output CSV files. This allows us to see where the FEM solution differs, in $\text{dB}$, most from the exact solution. The error field used to compute the norm are also saved in the CSV files.

We can inspect the error and the $\text{dB}$ error of the fields by importing it into Paraview. Select File, then Open. Then, navigate to the error subdirectory in the repository, where validate.jl has exported all the error fields as CSV files. Select error_..csv to import them all. This will import a data table into ParaView available in the Pipeline Browser as error*. Select it, then go Filters > Alphabetical > Table To Points. In the Properties tab, for the X Column, Y Column and Z Column combos select x, y and z respectively. You will now see all the nodal points of the mesh displayed in space. In the colouring combo at the top, you can select either error_field or dB_error_field to colour the points according to how big the selected error field is. My results are shown in the animations below.

It is possible to see that the error in $\text{Pa}$ (left panel) is very low except at round $66$ $\text{Hz}$ and $80$ $\text{Hz}$. When we look at the pattern of the error we clearly see something that is not random, and that has a striking resemblance to the modal patterns themselves. As we look at the frequencies at which the highest errors are recorded, we clearly see that there are other modal frequencies in close proximity. This suggests that the resulting field is deviating from the target not because model or computational errors, but because neighbouring modes are activated by the source. If we look at the $\text{dB}$ error field on the right we will see that the error is always very small a part for the highest frequency solution field. This tells us that the contribution from other modes, even though we can detect it at around $66$ $\text{Hz}$, starts becoming actually significant only at around $80$ $\text{Hz}$. So, for this rigid walled rectangular room, we can say that at $80$ $\text{Hz}$ the modal superposition cannot be considered pure.

By the way, note that if we had to drive the room at frequencies different than exact modal frequencies we will never get a pure superposition, not even at low frequency, as we will activate the modes closest to the driving frequency.

Now, you could do the same analysis as we did for any other source location. Will anything change? The animation below compares two identical simulations, the only difference being that in the one on the left the source is at the centre of the room.

This animation, still prepared with ParaView, shows the sound pressure level associated to the various solution fields. To give some context to the numbers, $40$ $\text{dB}$ SPL is the level you would record in a very quiet room, while $130$ $\text{dB}$ SPL is the threshold of pain. You can see that if we put the source in the middle the room is extremely quiet, only reaching high level SPL for the last frequency. Why?

Well, from the Acoustic Modes of a Rectangular Room episode you will probably remember that modes have nodal surfaces, surfaces in which the mode shape is $0$. If we put a source there, we can hardly activate the mode at all, as we are trying to drive it when it can only be $0$, so no wave can start. Turns out that the middle of the room is in a nodal surface for most of the low frequency modes we are considering. We can still see a bit of sound developing because our source is a sphere, not an ideal point, so it is driving air also just outside the nodal surface. The higher the frequency the sharper the variation of pressure close to the nodal surface, so high pressure can still be found close to the nodal surface. Hence, at higher frequency the sphere will become more able to drive pressure even by being located in a nodal surface due to its non infinitesimal extent. By contrast, the modes always have a peak on the corners of the room, so it is the ideal place to drive all modes or, if you want to control them, to put a bass trap in. This has the con that, in our simulation, the modal superposition gets dirty with many modes pretty soon in frequency.

Conclusion

This was another long one, it took more time than what I anticipated to finish. Still, we learned something valuable:

• Elmer does not have, at the moment, an Helmholtz Eigensolver. This complicates things for us when we want to study modal behaviour of enclosures, but we can work around it. So, even when your solver does not support what you want to do, you can be creative. It should be mentioned though that other Elmer solvers for acoustics do have an eigensolver, and these will be covered in future episodes.
• When we drive an enclosure with a sound source the steady state vibration of the enclosed air is a modal superposition. Even for rigid walled rooms, the superposition will be pure (one single mode) only at when driving the room at the lowest modal frequencies. For our room, we can expect one single mode only when driving at resonance frequencies lower than $80$ $\text{Hz}$.
• The Elmer Helmholtz solver solution might have an initial phase applied to it which can be hard to predict. In the solution presented in this article, they luckily didn’t.
• The response of a room depends on its shape (which dictates the shape of the modes) but also on the position of the source. Also, the pressure is not uniform in the room, meaning that the pressure recorded by different microphones at different points will be most likely different. In other words, the transmission line between a source and a receiver, when placed inside a room, depends on the locations of the source and receiver, even when they are omnidirectional!

In this episode we used ParaView a lot to make cool animations without mentioning how to do it. In the next episodes we will then have a look at ParaView and the basics of how to use it for postprocessing. Then, we will make more realistic rooms and investigate their response.

License Information

This work is licensed under a Creative Commons Attribution 4.0 International License.

• Physics
• Elmer
• FEM
• Linear System
• Single Physics
• Steady State
• Modal
• Series: Rigid Walled Room