## [SOLVED] Comparing with Comsol: many questions for you :D

Numerical methods and mathematical models of Elmer
CrocoDuck
Posts: 50
Joined: 12 May 2016, 13:15
Antispam: Yes

### [SOLVED] Comparing with Comsol: many questions for you :D

Hey cool people!

I suddenly remembered that I had data sampled from a Comsol simulation I run during my Master Course in Acoustics. So I decided to run the same simulation on Elmer and have a look at how the two solutions compare.

It is a simple model of a rectangular room (4 m X 5 m X 3 m) with a single window. An acoustic plane wave of given SPL is entering the window and the walls have a given normalized admittance a. The simulation was run on Comsol for values of frequency of 55 Hz, 100 Hz and 300 Hz. The solver to use is the Helmholtz Solver.

On Comsol I set the normal acceleration of the window boundary in accordance with the wave amplitude and frequency. On Elmer I used equation (9.7), where the normal velocity was chosen as the velocity relative to the acceleration as used in Comsol. I feel like this should be OK. I am not sure about the boundary conditions on the walls. I can calculate the specific acoustic impedance of the walls using the normalized admittance:

z = (rho * c) / a [rho = air density, c = phase speed of sound in air, at room temperature]

Which has units of metric Rayls: [kg s^-1 m^-2]. In Comsol I could use this value directly. On ElmerGUI I put this value in the "Real Part of the Impedance" in the Boundary Condition window. Here my first question:
• Is the "Real Part of the Impedance" in that window relative to the specific acoustic impedance of the walls? I suspect it could actually be instead the quantity Z in equation 9.6...
I meshed the volume with Salome, run Elmer and used ParaView to sample the data on the same grid I used to sample Comsol data. Then, I am using Matlab to compare the SPL in the room (I have only these data from Comsol and I have not access to Comsol anymore).

I found that I can find agreement within a couple of dB when the frequency is 55 Hz and 100 Hz. However, the meshing needs to be finer in Elmer. I still have to find a mesh fine enough to give satisfying similar results at 300 Hz, but seems to be getting there the finer the mesh. Still, differences are pretty big also with a mesh as twice as fine (maximum size 0.075 m VS 0.15 m). This is made harder by the fact that the problem hardly converges on the meshes I try. I found that running the solver in parallel helps convergence. This spans few questions:
• Imaging that I set up the Elmer model exactly as the Comsol model, is the need for higher meshing on Elmer expected?
• Do you have any rule of thumb or documentation to point me at in order to figure out how to make good meshing and help convergence? It seems to get harder the finer the mesh...
• I reckon that Comsol Pressure Wave solutions appeared to have small imaginary part, while Elmer ones seem to have it higher... Should I be afraid? I guess not, but it would be interesting for me to see if you have something to say about this.
• About the parallel solver, I clearly see that it changes convergence. Also, the solution has slight differences with respect the not parallel run. I imagined that it is unwise to slice the volume into parts that are too small to guarantee accurate solution of the governing equation... Do you have some suggestion to operate proper "slicing"? At the moment, I am using the ElmerGUI defaults, which seemed very appropriate to me after reading this.
This should be about it. I couldn't find an explanation on the manuals (although I feel it is hiding somewhere) so I hope you guys can point me in the right direction. Elmer appears to be much faster than Comsol on this problem, by the way. Still, the need of higher mesh density is now making the time to solve the problem similar... Not that I am disappointed (actually it is the other way around: Elmer is showing to be cooler than what I expected!).
Last edited by CrocoDuck on 22 Jul 2016, 11:13, edited 1 time in total.

CrocoDuck
Posts: 50
Joined: 12 May 2016, 13:15
Antispam: Yes

### Re: Comparing with Comsol: many questions for you :D

I almost forgot my sif file!

Code: Select all

``````Header
CHECK KEYWORDS Warn
Mesh DB "." "."
Include Path ""
Results Directory ""
! Simulation Frequency [Hz]
\$ f = 300
\$ a = 0.02
! SPL of incoming plane wave [dB re 20uPa]
\$ inspl = 80
End

Simulation
Max Output Level = 5
Coordinate System = Cartesian
Coordinate Mapping(3) = 1 2 3
Simulation Type = Scanning
Steady State Max Iterations = 1
Output Intervals = 1
Timestepping Method = BDF
BDF Order = 1
Timestep intervals = 1
Solver Input File = case.sif
Post File = case.vtu
End

Constants
Gravity(4) = 0 -1 0 9.82
Stefan Boltzmann = 5.67e-08
Permittivity of Vacuum = 8.8542e-12
Boltzmann Constant = 1.3807e-23
Unit Charge = 1.602e-19
End

Body 1
Target Bodies(1) = 2
Name = "Body 1"
Equation = 1
Material = 1
End

Solver 1
Equation = Helmholtz Equation
Variable = -dofs 2 Pressure Wave
Procedure = "HelmholtzSolve" "HelmholtzSolver"
Exec Solver = Always
Stabilize = True
Bubbles = False
Lumped Mass Matrix = False
Optimize Bandwidth = True
Steady State Convergence Tolerance = 1.0e-5
Nonlinear System Convergence Tolerance = 1.0e-7
Nonlinear System Max Iterations = 100
Nonlinear System Newton After Iterations = 3
Nonlinear System Newton After Tolerance = 1.0e-3
Nonlinear System Relaxation Factor = 1
Linear System Solver = Iterative
Linear System Iterative Method = BiCGStab
Linear System Max Iterations = 500
Linear System Convergence Tolerance = 1.0e-10
BiCGstabl polynomial degree = 2
Linear System Preconditioning = Diagonal
Linear System ILUT Tolerance = 1.0e-3
Linear System Abort Not Converged = False
Linear System Residual Output = 1
Linear System Precondition Recompute = 1
End

Equation 1
Name = "Helmholtz EQ"
Angular Frequency = Variable time; Real MATC "2 * pi * f(tx - 1)"
Active Solvers(1) = 1
End

Material 1
Name = "Air (room temperature)"
Relative Permeability = 1.00000037
Heat Conductivity = 0.0257
!Turbulent Prandtl Number = 0.713 Need to comment this or simulation won't work
Heat Capacity = 1005.0
Relative Permittivity = 1.00059
Relative Permeability = 1.00000037
Viscosity = 1.983e-5
Viscosity = 1.983e-5
Sound speed = 343.0
Heat expansion Coefficient = 3.43e-3
Relative Permittivity = 1.00059
Porosity Model = Always saturated
Relative Permittivity = 1.00059
Density = 1.21
Relative Permeability = 1.00000037
End

Boundary Condition 1
Target Boundaries(6) = 1 2 3 4 5 6
Name = "Walls"
Wave impedance 1 = Real MATC "1.21 * 343.0 / a"
End

Boundary Condition 2
Target Boundaries(1) = 7
Name = "Window"
Wave Flux 2 = Variable time; Real MATC "2 * sqrt(2) * pi * f(tx - 1) * 0.00002 * pow(10, inspl / 20) / 343"
End
``````
Let me know if you want to have a look at the mesh as well...

CrocoDuck
Posts: 50
Joined: 12 May 2016, 13:15
Antispam: Yes

### Re: Comparing with Comsol: many questions for you :D

Update: after reading this I tried BiCGStabl with l = 2,3. It helped convergence a great deal: it us much faster and monotonic. I also raised the mesh density to the point that the maximum size is 5 cm. Results are approaching the Comsol results. Still a bit of a mystery why I need such a bigger mesh density... I wonder whether there is something I can do in the solver settings.

CrocoDuck
Posts: 50
Joined: 12 May 2016, 13:15
Antispam: Yes

### Re: Comparing with Comsol: many questions for you :D

Hey there! New update:

Using ILUT as a preconditioner, with tolerance 1e-3, improves convergence greatly. It is still a pretty fast solution. I manage to get a good agreement with Comsol using a low mesh density (0.15 m maximum element size) but second order. Now the agreement is back to 2 dB also at 300 Hz. It appears that Comsol sneaked in an higher order meshing without telling it to me... the cons of user friendliness I would say. Elmer took 10 minutes for giving me the solution. Not bad at all!

I will play with this more and mark as Solved if everything works as expected. Feel free to reply here if you think you have something to suggest/comment!

CrocoDuck
Posts: 50
Joined: 12 May 2016, 13:15
Antispam: Yes

### Re: Comparing with Comsol: many questions for you :D

Solved indeed! I also made a little more complex model adding a plate of glass in front of the window (modelled just by its impedance) to realize a little plenum window. Again, very good agreement with Comsol data overall. Noticeable Disagreement is where SPL is very low or very high, as one would expect. The final .sif file:

Code: Select all

``````Header
CHECK KEYWORDS Warn
Mesh DB "." "."
Include Path ""
Results Directory ""
! Simulation Frequency [Hz]
\$ f = 300
\$ a = 0.02
! SPL of incoming plane wave [dB re 20uPa]
\$ inspl = 80
End

Simulation
Max Output Level = 5
Coordinate System = Cartesian
Coordinate Mapping(3) = 1 2 3
Simulation Type = Scanning
Steady State Max Iterations = 1
Output Intervals = 1
Timestepping Method = BDF
BDF Order = 1
Timestep intervals = 1
Solver Input File = case.sif
Post File = case.vtu
End

Constants
Gravity(4) = 0 -1 0 9.82
Stefan Boltzmann = 5.67e-08
Permittivity of Vacuum = 8.8542e-12
Boltzmann Constant = 1.3807e-23
Unit Charge = 1.602e-19
End

Body 1
Target Bodies(1) = 2
Name = "Body 1"
Equation = 1
Material = 1
End

Solver 1
Equation = Helmholtz Equation
Variable = -dofs 2 Pressure Wave
Procedure = "HelmholtzSolve" "HelmholtzSolver"
Exec Solver = Always
Stabilize = True
Bubbles = False
Lumped Mass Matrix = False
Optimize Bandwidth = True
Steady State Convergence Tolerance = 1.0e-5
Nonlinear System Convergence Tolerance = 1.0e-7
Nonlinear System Max Iterations = 100
Nonlinear System Newton After Iterations = 3
Nonlinear System Newton After Tolerance = 1.0e-3
Nonlinear System Relaxation Factor = 1
Linear System Solver = Iterative
Linear System Iterative Method = BiCGStabl
Linear System Max Iterations = 500
Linear System Convergence Tolerance = 1.0e-10
BiCGstabl polynomial degree = 3
Linear System Preconditioning = Diagonal
Linear System ILUT Tolerance = 1.0e-3
Linear System Abort Not Converged = False
Linear System Residual Output = 1
Linear System Precondition Recompute = 1
End

Equation 1
Name = "Helmholtz EQ"
Angular Frequency = Variable time; Real MATC "2 * pi * f(tx - 1)"
Active Solvers(1) = 1
End

Material 1
Name = "Air (room temperature)"
Relative Permeability = 1.00000037
Heat Conductivity = 0.0257
!Turbulent Prandtl Number = 0.713
Heat Capacity = 1005.0
Relative Permittivity = 1.00059
Relative Permeability = 1.00000037
Viscosity = 1.983e-5
Viscosity = 1.983e-5
Sound speed = 343.0
Heat expansion Coefficient = 3.43e-3
Relative Permittivity = 1.00059
Porosity Model = Always saturated
Relative Permittivity = 1.00059
Density = 1.21
Relative Permeability = 1.00000037
End

Boundary Condition 1
Target Boundaries(6) = 1 2 3 4 5 10
Name = "Walls"
Wave impedance 1 = Real MATC "1.21 * 343.0 / a"
End

Boundary Condition 2
Target Boundaries(1) = 11
Name = "Window"
Wave Flux 2 = Variable time; Real MATC "2 * sqrt(2) * pi * f(tx - 1) * 0.00002 * pow(10, inspl / 20) / 343"
End

Boundary Condition 3
Target Boundaries(4) = 6 7 8 9
Name = "Glass"
Wave impedance 1 = 13000000
End
``````
There is still something I am curious about tho. Like this: by sorting the SPL values on the sample points I found the sort of behaviour you can see in the attached picture. What is producing the drop in the Elmer Curve? Is it a numerical thing?
Attachments
Sorted300HzPLE.png

Cfirminy
Posts: 14
Joined: 16 Sep 2015, 10:54
Antispam: Yes

### Re: [SOLVED] Comparing with Comsol: many questions for you :D

Hi Crocoduck,

Your study is very interesting, could you send your mesh file ?

CFirminy

CrocoDuck
Posts: 50
Joined: 12 May 2016, 13:15
Antispam: Yes

### Re: [SOLVED] Comparing with Comsol: many questions for you :D

Cfirminy wrote:Hi Crocoduck,

Your study is very interesting, could you send your mesh file ?

CFirminy
Unfortunately the file is to big to be sent over the forums. I will try to provide a link after work!

Cfirminy
Posts: 14
Joined: 16 Sep 2015, 10:54
Antispam: Yes

### Re: [SOLVED] Comparing with Comsol: many questions for you :D

Ok that would be great if you can! Thanks

CrocoDuck
Posts: 50
Joined: 12 May 2016, 13:15
Antispam: Yes

### Re: [SOLVED] Comparing with Comsol: many questions for you :D

Here a link. Hope it helps.

pilafa
Posts: 56
Joined: 17 Apr 2012, 00:33
Antispam: Yes
Location: Barcelona

### Re: [SOLVED] Comparing with Comsol: many questions for you :D

Still a bit of a mystery why I need such a bigger mesh density.
I observe the same here.

To fit Comsol results, Elmer needs a much (much!) finer mesh (even if quadratic mesh is used). Any explanation?

Otherwise, the solution is very "noisy". Comsol must do something to improve that because it can find a smooth solution is few seconds.

Any idea of what could improve that in Elmer?

Thanks,
PL