In this episode we will build a model of a pulsating sphere source. The pulsating sphere source is an ideal source which forms the base for the development of point sources. In essence, a point source is a pulsating sphere in the limit of $a$, the radius of the sphere, approaching $0$. For this reason, although abstract, the pulsating sphere is a very powerful theoretical tool that enables the study of point sources which in turn, through integration and wave propagation principles, enable to study of any arbitrary acoustic field source.

We will study the pulsating sphere source radiating in infinite space in the domain of FEM analysis. Again, the study of simplified problems allows us to understand, through comparison with theory, how to tune our FEM problem to maximise accuracy. We will be then able to carry over this knowledge to more advanced problems.

As discussed in The FEM Pipeline episode, we will build a development pipeline to search for the best way to solve single physics steady state acoustic problems with Elmer. Hence we will:

• Dig out a theoretical model for our system from literature.
• Use the Julia programming language to create an exact model.
• Use Elmer to create an approximate model.

We will then study the effect on accuracy of few parameters, mainly mesh size.

Project Files

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

• Coarse Mesh.
• Fine Mesh.

Theoretical Model

Many equivalent formulations of the field radiated by a pulsating sphere exist. In the context of this post we will make use of the formulation presented in the first chapter of Active Control of Sound by P.A. Nelson and S.J. Elliot, to which you can refer for more information, although similar derivations will be available for sure in many other books and websites. We will not derive the result, but just take the bits we need.

When we think about acoustic wave propagation, two quantities are of interest. One is the pressure disturbance, or acoustic pressure, $p$, which is the pressure disturbance with respect the equilibrium pressure $p_{0}$. The other is the particle velocity $u$, which is the velocity of an element of fluid deriving from the propagation of the acoustic field (remember that we are talking about continuum fluid dynamics, so our fluid is modelled as a continuum and particles are not molecules, but infinitesimal fluid elements). Both $p$ and $u$ will change in space and time, and the amplitude and phase relationship between them will also change in space and time. The most convenient mathematical way to represent them is then complex numbers, as they “encode” amplitude and phase information simultaneously in a single complex value. Hence:

$$p=p\left(x,y,z,t\right)\in\mathbb{C}$$

$$u=u\left(x,y,z,t\right)\in\mathbb{C}$$

Both the relationships above are valid for every $x$, $y$, $z$ and $t$. The amplitude and phase relationship between acoustic pressure and particle velocity is nontrivial, and depends on the nature of the source and domain and, in essence, to the field itself. To capture this relationship, one can define the specific acoustic impedance as follows:

$$z=\frac{p}{u}$$

Since $p$ and $u$ are both complex functions of space and time so is the specific acoustic impedance. It is possible to show that, for spherical waves propagating in an infinite uniform medium governed by the linear wave equation, the specific acoustic impedance is:

$$z\left(r\right)=\rho_{0}c_{0}\frac{jkr}{1+jkr}$$

Where $r=\sqrt{x^{2}+y^{2}+z^{2}}$ is the distance from the centre of the source emanating the spherical wave, such is our pulsating sphere. $j$ is the imaginary unit and $k$ is the wavenumber:

$$k=\frac{2\pi f}{c_{0}}$$

With $f$ is frequency at which the source is pulsating, and $c_{0}$ the phase speed of sound in the medium. Finally, $\rho_{0}$ is the equilibrium density of the medium. Note that, for a spherical wave propagating in free field (infinite space without any other obstacle) the specific acoustic impedance does not depend on time.

The specific acoustic impedance is an ingredient which defines the sound radiated from a spherical source, as we will see just below, but it is also a property that we can use to make a finite domain to behave as an infinite one when we do FEM, as we will see later.

As we said, a pulsating sphere will emit a spherical wave, which then will be characterised by the specific acoustic impedance defined above. Let’s assume that we have a sphere of radius $a$ pulsating at the frequency $f$ with a surface velocity $U$. Note that the surface velocity is, by definition, orthogonal to the surface of the source. Then the pressure disturbance is given by:

$$p=p\left(r,t\right)=z\left(a\right)\frac{aU}{r}\exp\left(j \left[ \omega t-k\left( r-a\right) \right] \right)$$

where $\omega=2\pi f$ as the angular frequency. To be noted that, due to the symmetry of the problem, the spatial dependency of the acoustic pressure field is radial, that is, points at the same distance $r$ from the source centre have the same pressure disturbance.

Source and Medium Parameters

For the rest of the post, we will assume the following source and medium parameters. You can repeat the study with any choice of parameters.

The chosen medium parameters above are that for air in equilibrium at ordinary room temperature ($20$ $\text{°C}$).

Julia Model

The code implementing the equations above is implemented in the AcousticModels.jl package. The following paragraphs illustrate how to use the functions. These functions are used by the validate.jl scripts in the repositories listed at the beginning of the article.

Open Julia by issuing the command julia. To install the package, hit the ] key to enter package mode and issue:

add https://gitlab.com/computational-acoustics/acousticmodels.jl.git

After installation you can use an using statement to make the functions available in your environment:

using AcousticModels

For example, given the parameters in the table above, to compute the impedance at $1$ $\text{m}$ from the source:

spherical_impedance(1000, 1, 1.205, 343)
412.08694629151245 + 22.49588634867694im

Or, to compute the impedance at a set of equally spaced distances between $0.005$ $\text{m}$ and $1$ $\text{m}$, with a resolution of $0.001$ $\text{m}$:

spherical_impedance.(1000, 0.005:0.001:1, 1.205, 343)
996-element Array{Complex{Float64},1}:
  3.438464634081654 + 37.54125692082816im
  4.933330796152837 + 44.88520764425224im
 6.6859932275747145 + 52.14133471072512im
  8.689694549377409 + 59.296461203974424im
 10.936816885631023 + 66.33804720893404im
 13.418952281912064 + 73.25425572657367im
 16.126978587079495 + 80.03400974817295im
 19.051139570106578 + 86.66704005848428im
                    ⋮
  412.0721206647585 + 22.63086218868773im
 412.07461016560563 + 22.60825418871258im
 412.07709220166583 + 22.585691177187382im
 412.07956680272395 + 22.563173020371515im
  412.0820339984162 + 22.540699585052277im
 412.08449381823146 + 22.51827073854228im
 412.08694629151245 + 22.49588634867694im
 

The function spherical_wave can be similarly used. For example, this is the complex pressure disturbance at $1$ $\text{m}$ and $0$ $\text{s}$:

spherical_wave(0.75, 0.005, 1000, 0, 1, 1.205, 343) 
-0.07164793924384269 + 0.12186780546373502im

And at the previous set of equally spaced distances:

spherical_wave.(0.75, 0.005, 1000, 0, 0.005:0.001:1, 1.205, 343)
996-element Array{Complex{Float64},1}:
   2.5788484755612404 + 28.15594269062112im
   2.5784638811553324 + 23.419984385921428im
    2.577447564216482 + 20.030421070628968im
    2.575928568123367 + 17.48236749878261im
    2.573978809358829 + 15.495333621210861im
   2.5716416698416644 + 13.901026983151116im
   2.5689449933704016 + 12.592353813596183im
    2.565907583772156 + 11.497919569528525im
                      ⋮
 -0.08509370871416262 + 0.11395722009711984im
 -0.08290863436648038 + 0.11538071065217743im
 -0.08070015552278063 + 0.11676266547865974im
 -0.07846905501937268 + 0.11810270747393418im
 -0.07621612162316302 + 0.11940047416733444im
   -0.073942149760358 + 0.12065561781066901im
 -0.07164793924384269 + 0.12186780546373502im

Building the FEM Model

We will now build the FEM model. As always we start with the Geometry.

Geometry

This is really very simple to build. 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 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 we will keep the radius as 5 mm.
• Click again on the Create a sphere solid button. This will create a new sphere.
• On the Combo View on the left, select the newly created sphere. On the bottom you will see the sphere properties. Click on the Radius value, type 100 and press Enter. This bigger sphere will model the extent of air that we include in the model.
• Now, click in any empty space in the Combo View to un-select all objects. Then, while holding Ctrl, select the bigger sphere first and the smaller 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 bigger sphere 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.

So, it will be clear to you that essentially our geometry is a big sphere with a spherical hole in the centre. So, in essence, we are only modelling the body of air around the source, up to $0.1$ $\text{m}$, omitting the source itself. This because the source is, in reality, nothing but a boundary condition for the air in contact with it. This is how we will model the source with Elmer.

Meshing

We will now mesh our geometry. Open your Salome installation 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 two new face objects under the solid object tree. One is the inner face, the other is the outer face.

After this geometry manipulation, we are ready to mesh. As a rule of thumb, our maximum mesh size $s$ should always be smaller, along every direction, than one tenth of the wavelength $\lambda$. Hence:

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

At $1$ $\text{kHz}$ this “critical size” is $3.43$ $\text{cm}$. We will see that this size is already sufficient to provide accurate results, but we can also run a higher density mesh study to see what we gain from it. Let’s then move on and create our mesh.

• 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.
  • For a coarse mesh, proceed as follows. Type 29 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.
  • For a fine mesh, proceed as follows. Type 5 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.

Both the coarse and fine mesh parameters outlined above produce meshes with maximum sizes smaller than a tenth of a wavelength.

Pictures of the coarse and fine mesh are shown below. They are obtained with Salome. If you want to get similar views, just right click on your mesh in the 3D view and choose Clipping. A window will popup allowing to introduce any clipping plane you might like. Note that this clipping plane is for visualisation only, and does not affect the mesh. You can see that NETGEN meshes adapt their size where needed, and are finer where the curvature is higher (like in the centre).

Solving

It is finally the time to solve our study. For this study we will use a little more advanced feature of Elmer, which is MATC expressions. MATC is a library for the evaluation of numerical expressions that is included with Elmer. Its documentation is available here. With MATC it is possible to add flexibility to Elmer and do all sorts of cool things, like parametric sweep studies or post-processing. You can run MATC expressions on the fields, at each solver iteration, but beware that MATC is pretty slow. It is best used to set up the simulation parameters, as we will do here.

Also, be aware that recent Elmer versions can employ Lua for the same scripting tasks, but faster. For more information, see here

Understanding the Problem

But first, what equation we have to solve? Look again at our theoretical solution:

$$p\left(r,t\right)=z\left(a\right)\frac{aU}{r}\exp\left(j\left[\omega t-k\left( r-a\right) \right]\right)$$

This solution holds for a spherical source at steady state, and it is found by using the linear wave equation which we already met in the What is Acoustic Modelling episode. Turns out that Elmer can solve for that equation. But also, notice this property:

$$p\left(r,t\right)=z\left(a\right)\frac{aU}{r}\exp\left(-jk\left[r-a\right]\right)\exp\left(j\omega t\right)$$

Thanks to the properties of exponential functions, we can factor the function in one part, up to the first exponential, that depends on space only ($r$) and a second part (the second exponential) that depends on time only which is just a complex sine wave at our source frequency. When this happens, we can actually use a variant of the linearised wave equation that allows to solve for the spatial part only of our disturbance. This equation is the Helmholtz equation that we already encountered in the Acoustic Modes of a Rectangular Room episode, and we will use the Elmer Helmholtz model.

Note that this factorisation is typical of steady state propagation, and it is possible to demonstrate that the time dependant part is always the same complex exponential function for every problem (for this type of PDE). So you should use the Helmholtz solver every time you are after a steady state solution, so that time will not complicate your analysis and potentially introduce more errors. If you need to model steady state time behaviour all you have to do is multiply your final space solution by $\exp\left(j\omega t\right)$. Time dependent solvers are best suited to study transients, that is, the evolution between steady states.

Since we will be using the Helmholtz model, we should always refer to its documentation here.

Now, in order to set up the model properly we will use a FEM trick. If you were reading carefully, you must have observed that our theoretical model is for a spherical source radiating in free field, which is an infinite space. However, our geometry is a ball of air of $0.1$ $\text{m}$ radius! How can we simulate our spherical source realistically with such a small volume?

Well, remember the spherical wave specific acoustic impedance defined above? That holds for a spherical wave propagating in an infinite medium. So, at a certain distance $r$ from the source the impedance will be $z\left(r\right)$. Now, imagine there was a wall in our infinite medium, at some distance from the sphere source. A wall is some object that offer its own impedance to the wave. This will cause part of the wave to be transmitted into the wall, and part to be reflected. So, impedance controls what happens to the field. Let’s think about the outer boundary of our domain again. If we assign to it the very same impedance that the infinite spherical wave would have at that radius, that is $z\left(r\right)$, with $r$ equal to $0.1$ $\text{m}$, than the outer boundary behaves exactly as if it was perfectly transparent to the wave and there was infinite space behind it. This because its impedance matches exactly that of the wave, meaning that no reflection happens, and transmission through the boundary is the same as that in free field. Hence, we will do just that: we will set the impedance of the outer layer of our domain to that of the wave, so that our finite domain behaves exactly as an infinite one.

By referring to the Helmholtz model manual, we can see that the model allows two types of boundary conditions, Flux and impedance. Clearly, we want the second one for our outer surface, but note that the quantity $Z$ evidenced in the manual is not our specific acoustic impedance $z$. However, thanks to some dimentional considerations, it is easy to calculate $Z$: it just the specific acoustic impedance divided by the density of the medium:

$$Z\left(r\right)=\frac{z\left(r\right)}{\rho_{0}}=c_{0}\frac{jkr}{1+jkr}$$

So, all we need to do is calculate this for $r$ equal to $0.1$ $\text{m}$ (the radius at which our outer spherical boundary is placed) and use it as boundary condition. Since Elmer requires the quantity to be input as a real-imaginary parts pair, here the same quantity but with the two parts separated:

$$Z\left(r\right)=c_{0}\frac{\left(kr\right)^{2}}{1+\left(kr\right)^{2}}+jc_{0}\frac{kr}{1+\left(kr\right)^{2}}$$

This is an expression we will evaluate with MATC.

Finally, we need to assign a boundary condition at the inner surface of our medium, that in contact with the spherical source. This boundary condition must be that a spherical source would produce. The Elmer model manual already gives us the Flux boundary condition that we need:

$$g=j\omega\rho_{0}U$$

We will implement this too as a MATC expression.

Model Setup

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 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.
  • 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 no ugly solution will be saved if the solver cannot converge.
    • 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 below.
  • In the Helmholtz Equation replace the value of Sound speed with the MATC expression in Snippet 4 below.
  • 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 field radiation condition, so I suggest you edit the Name edit box at the bottom and type Radiator.
  • Then, select the Helmholtz Equation tab.
  • Under Imag part of the flux, insert the code in Snippet 5 below.
  • Click OK.
• Select Model from the top menu, then Boundary condition then Add.
  • This will be our far field condition, so I suggest you edit the Name edit box at the bottom and type Outlet.
  • Then, select the Helmholtz Equation tab.
  • Under Real part of the impedance, insert the code in Snippet 6 below.
  • Under Imag part of the impedance, insert the code in Snippet 7 below.
  • Click OK.
• Select Model from the top menu, then Set boundary properties.
• Zoom into the domain with the mouse wheel until the view penetrates the outer surface. The inner surface should now be visible. Double click on it. A window will popup. From the Boundary condition combo, select Radiator.
• Zoom out until the outer surface of the domain is visible. Double click on it. A window will popup. From the Boundary condition combo, select Outlet.
• 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 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!

Snippet 1

Variables of the simulation. From top to bottom: frequency ($\text{Hz}$), source surface velocity ($\frac{\text{m}}{\text{s}}$), phase speed of sound in the medium ($\frac{\text{m}}{\text{s}}$), density of the medium ($\frac{\text{kg}}{\text{m}^3}$) and radius at which the domain terminates ($\text{m}$).

$ f = 1000.0
$ U = 0.75
$ c = 343.0
$ p = 1.205
$ r = 0.1

Snippet 2

This code instructs to assign the Real value resulting from the MATC expression between quotes to the required parameter. In this case, the frequency.

Frequency = Real MATC "f"

Snippet 3

This code instructs 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 instructs to assign the Real value resulting from the MATC expression between quotes to the required parameter. In this case, the medium phase speed of sound.

Real MATC "c"

Snippet 5

This code instructs to assign the Real value resulting from the MATC expression between quotes to the required parameter. In this case, the imaginary part of the Flux boundary condition. This is the only part of the Flux boundary condition for this problem. This expression is the same as that outlined in the section Understanding the Problem above.

Real MATC "2 * pi * f * p * U"

Snippet 6

This code instructs to assign the Real value resulting from the MATC expression between quotes to the required parameter. In this case, the real part of the impedance boundary condition. This expression is the same as that outlined in the section Understanding the Problem above.

Real MATC "((2 * pi * f * r)^2) / (c * (1 + (((2 * pi * f * r)^2) / (c^2))))"

Snippet 7

This code instructs to assign the Real value resulting from the MATC expression between quotes to the required parameter. In this case, the imaginary part of the impedance boundary condition. This expression is the same as that outlined in the section Understanding the Problem above.

Real MATC "(2 * pi * f * r) / (1 + (((2 * pi * f * r)^2) / (c^2)))"

Results

Once the simulation is finished, we can use ParaView to look at the results. Open ParaView.

Select File from the top menu and then Open. Navigate to your ElemerFEM project directory. In the directory, select the vtu file and then click OK.

In the properties browser on the left, you will see that the vtu file contains the pressure wave 1 and pressure wave 2 fields. They are the real and imaginary part of the complex acoustic pressure field respectively. Click on the Apply button. Now the data is loaded into ParaView.

There are tons of good things you can do with ParaView. Here I will just show you how we can export the field values along a line. If you look at the exact solution for the field, you will see that it depends only from $r$. So, to do a simple comparison with theory, we can simply export the field along a radius of the spherical domain.

• From the top menu, select Filters.
• Select Alphabetical and then Plot Over Line.
• In the Line Parameters, under the Properties tab on the left, input the values in the line parameters table below. This will make a radial line along the $x$ axis, extending from the source to the edge of the domain.
• Click Apply. A new view will open showing a plot of all the fields in the vtu file as a function of the position along the line.

Line Parameters Table

This data can be exported. To do so, select the result of the Plot Over Line filter from the Pipeline Browser on the left. Then, click on File from the top menu and choose Save Data. If you follow the steps below, you will be able to use the Julia code validate.jl to calculate the error between the fields sampled along the line with the theoretical expected values.

• In the File name edit box type a name with csv extension, for example linedata.csv.
• Click OK.
• In the Precision editbox, type 30. This is actually overkill, but it will ensure that all digits of our numbers are written to file.
• Then tick the Choose Arrays To Write tickbox.
• Un-select everything but pressure wave 1 and pressure wave 2.
• Click OK.

You could now edit validate.jl so that it points to your csv file. I suggest you put the csv file in the same directory of validate.jl and change line 14 with the name of your file. For example, if you called your file linedata.csv:

linedata = readdlm("linedata.csv", ',')

You can run this file from the Julia REPL by issuing this command:

include("validate.jl")

This will work only if Julia was open in the same directory as validate.jl.

The script will produce a plot and hang after it displays it until Enter is pressed. It might take some time to run the code for the first time, as Julia will have to compile it. The plots show various metrics of error. I will present my results below.

Coarse Mesh

Fine Mesh

Discussion

The best way to compute the error is by taking the ratio of the FEM complex field and the exact complex field. This will give us an error field which, in the case the FEM solution was the same as the exact one, would be $1$ on all over. This means that its magnitude would be $1$ and its phase would be $0$. Bigger magnitude means bigger amplitude for the FEM solution (with respect the exact one) while bigger phase means a higher phase for the FEM solution (with respect the exact one).

The plots above show the magnitude of this ratio (expressed in $\text{dB}$) in blue and the phase in orange.

It is possible to see that the error is in all cases vanishing small, only a tiny fraction of a $\text{dB}$ or a $\text{rad}$. These errors are, in all cases, way smaller than the highest accuracies one can typically hope to obtain in real life measurements ($0.1$ $\text{dB}$ for magnitude, whilst phase typically is extremely sensitive to temperature). Hence, if we had a real spherical source in an anechoic space the values we would measure would agree with all $3$ of our models: the exact one, the coarse FEM one and the fine FEM one, telling us that they are all equivalent.

We can then conclude that, despite the gain in accuracy by using the fine mesh is definitely measurable (the error for the finer mesh is indeed lower, especially for phase), the coarse solution is already very accurate, suggesting that using mesh sizes much lower than one tenth of a wavelength does not bring particular advantages.

Conclusion

This was quite a long one. However, we learned a few things:

• How to use MATC expression to configure Elmer.
• How to setup impedance boundary conditions with Elmer the right way.
• How to setup source boundary conditions with Elmer the right way.
• That the solution will most likely be good already with a mesh size just smaller than one tenth of a wavelength.
• That we can simulate infinite spaces with finite geometries if we are clever with the boundary conditions.

In reality, the symmetry of the problem would have allowed for many more FEM hacks, like solving the problem only on one eight of the domain. Also, we checked the solution only along one direction and it would be more interesting to check it in the entire domain. But we will do that in the next episodes, when we get back to the modes of a room.

But perhaps the next episode should focus on more solver details. Why we set the solver settings the way we did? What is the convergence plot that Elmer shows to us during the solution process? How do we use it?

License Information

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

• Physics
• Elmer
• FEM
• Linear System
• Single Physics
• Steady State