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

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

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

Post by CrocoDuck »

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: 79
Joined: 12 May 2016, 13:15
Antispam: Yes

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

Post by CrocoDuck »

I almost forgot my sif file!

Code: Select all

Header
  CHECK KEYWORDS Warn
  Mesh DB "." "."
  Include Path ""
  Results Directory ""
! Simulation Frequency [Hz]
  $ f = 300
  ! Walls normalized admittance:
  $ 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: 79
Joined: 12 May 2016, 13:15
Antispam: Yes

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

Post by CrocoDuck »

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: 79
Joined: 12 May 2016, 13:15
Antispam: Yes

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

Post by CrocoDuck »

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: 79
Joined: 12 May 2016, 13:15
Antispam: Yes

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

Post by CrocoDuck »

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
  ! Walls normalized admittance:
  $ 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
(21.58 KiB) Not downloaded yet
Cfirminy
Posts: 14
Joined: 16 Sep 2015, 10:54
Antispam: Yes

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

Post by Cfirminy »

Hi Crocoduck,

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

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

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

Post by CrocoDuck »

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

Post by Cfirminy »

Ok that would be great if you can! Thanks
CrocoDuck
Posts: 79
Joined: 12 May 2016, 13:15
Antispam: Yes

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

Post by CrocoDuck »

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

Post by pilafa »

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
Post Reply