Acoustics - Large-amplitude Wave Motion in Air

Numerical methods and mathematical models of Elmer
mzenker
Posts: 1999
Joined: 07 Dec 2009, 11:49
Location: Germany

Re: Acoustics - Large-amplitude Wave Motion in Air

Post by mzenker »

Hi Zoltan,

I did a quick search and found this: https://arxiv.org/abs/1707.05876. Just quickly browsed it and saw that apparently they do not expect so large differences in bulk viscosity when humidity changes. But not being an expert in the field I may have overlooked things. Also I didn't take the time to study it in detail. But maybe it helps a bit.

Matthias
Z.Losonc
Posts: 40
Joined: 21 Dec 2017, 16:05
Antispam: Yes

Re: Acoustics - Large-amplitude Wave Motion in Air

Post by Z.Losonc »

Hi Matthias,

Thanks for the link. Unfortunately it does not cure the problem. The most reasonable average value for the bulk viscosity at low frequencies is between 1e-3 and 1e-2. Let's see the results when using Bulk Viscosity = Real 5.0e-3.

The geometry is a 2D wave guide 1m long, 0.1m high, filled with air. It uses a 40x4 second order mesh. The simulation is done in time domain using 200x4.54545e-4 s time steps. The input velocity is a sinusoidal time function of 1 m/s amplitude and 34.3 Hz frequency. Here is a snapshot of the acoustic velocity at step 199:
velo_x.png
velo_x.png (21.9 KiB) Viewed 5764 times
And here is how the wave propagates in time and space:
step-50.png
step-50.png (5.43 KiB) Viewed 5764 times
time step 50


step-100.png
step-100.png (6.15 KiB) Viewed 5764 times
time step 100

Only 3 attachments are allowed per post, therefore this text will continue in the next post below…

Zoltan
Last edited by Z.Losonc on 13 Feb 2018, 16:34, edited 1 time in total.
Z.Losonc
Posts: 40
Joined: 21 Dec 2017, 16:05
Antispam: Yes

Re: Acoustics - Large-amplitude Wave Motion in Air

Post by Z.Losonc »

...continued from the previous post:
step-150.png
step-150.png (6.73 KiB) Viewed 5764 times
time step 150

step-185.png
step-185.png (7.11 KiB) Viewed 5764 times
time step 185

The last image shows the moment when the signal front reaches the 1m distance from the input. The time it took to get there is 185*4.54545e-4=0.084090825 s, which corresponds to a sound speed of 11.89 m/s. The expected speed of sound is about 343 m/s, so this is much below the expected value.

The pdf of Allan J. Zuckerwar was used as a possible source of higher values for the bulk viscosity, because these realistic low values lead to very low sound speeds in the simulation. Even when a value of 0.635 was used (which is 127 times greater than the value used in this simulation) the speed of sound was still only about 18 m/s. Something is not right either with the sif file, or with the solver. The question is only, where is the problem?

Here is the (well guarded 'secret') sif file for those who want to experiment with this 'advanced' solver:

Code: Select all

Header
  CHECK KEYWORDS Warn
  Mesh DB "." "."
  Include Path ""
  Results Directory ""
End

Simulation
  Max Output Level = 10
  Coordinate System = "Cartesian 2D"
  Coordinate Mapping(3) = 1 2 3
  Steady State Max Iterations = 1
  BDF Order = 2

  Simulation Type = Transient
  Timestepping Method = BDF

  Timestep Intervals = 200
  Timestep Size = 4.54545e-4
  Output Intervals = 1

  Adaptive Timestepping = False
  Adaptive Time Error = 0.05
  Adaptive Keep Smallest = 3
  Adaptive Min Timestep = 1.0e-6
  Adaptive Max Timestep = 0.01

$w = 215.5

  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) = 1
  Name = "Body Property 1"
  Equation = 1
  Material = 1
  Initial Condition = 1
!  Body Force = 1
End

Equation 1
  Name = "CompressibleNS"
  Active Solvers(1) = 1
End

Solver 1
  Exec Solver = String "always"
  Equation = "CompressibleNS"
  Procedure = "CompressibleNS" "CompressibleNS"
  Variable = -dofs 4 Velocity[Velo:2 Dens:1 Temp:1]
!  Element = p:1 b:1
   Element = p:2

  Stabilize = True
  Bubbles = False

  Steady State Convergence Tolerance = 1.0e-5

  Nonlinear System Convergence Tolerance = 1.0e-5
  Nonlinear System Max Iterations = 30
  Nonlinear System Newton After Iterations = 3
  Nonlinear System Newton After Tolerance = 1.0e-5
  Nonlinear System Relaxation Factor = 1
  Nonlinear System Convergence Measure = solution   

  Linear System Solver = Iterative
  Linear System Iterative Method = BiCGStab
  Linear System Max Iterations = 500
  Linear System Convergence Tolerance = 1.0e-5
  Linear System Preconditioning = ILU1
  Linear System ILUT Tolerance = 1.0e-4
  Linear System Abort Not Converged = False
  Linear System Residual Output = 10
  Linear System Precondition Recompute = 1
End

Material 1
  Name = "Air (room temperature)"
  Heat expansion Coefficient = Real 3.43e-3
  Heat Conductivity = Real 0.0257
  Equilibrium Density = Real 1.205
  Equilibrium Temperature = Real 293.0
  Specific Heat = Real 0.718
  Specific Heat Ratio = Real 1.401
  Viscosity = Real 1.81e-5
  Bulk Viscosity = Real 5.0e-3
End

Boundary Condition 1
  Target Boundaries(1) = 2
  Name = "Out"
  Velo 1 = Real 0.0
  Velo 2 = Real 0.0
End

Boundary Condition 2
  Target Boundaries(2) = 1 3
  Name = "Walls"
!  Velo 1 = Real 0.0
  Velo 2 = Real 0.0
End

Boundary Condition 3
  Target Boundaries(1) = 4
  Name = "In"
  Velo 1 = Variable Time
	Real MATC "1.0*sin(w*tx)"
!  Velo 2 = Real 0.0
End

Initial Condition 1
  Name = "Initial Condition 1"
  Velo 1 = Real 0.0
  Velo 2 = Real 0.0
!  Dens = Real 0.0
!  Temp = Real 0.0
End
The project file with the mesh is also attached to this post.

If anybody knows the solution to this riddle, please don't hesitate to enlighten us. ;-)

Zoltan
Attachments
CompressibleNS.7z
(18.72 KiB) Downloaded 336 times
Z.Losonc
Posts: 40
Joined: 21 Dec 2017, 16:05
Antispam: Yes

Re: Acoustics - Large-amplitude Wave Motion in Air

Post by Z.Losonc »

Hi,

Found the cause of wrong results; the “Specific Heat” parameter was the culprit, not the “Bulk Viscosity”. The specific heat capacity is usually given in tables as kJ/(kg*K), but the sif file requires pure SI units in J/(kg*K). Therefore instead of
Specific Heat = Real 0.718 we should write
Specific Heat = Real 718

If (for the sake of validation) we want to approximate small signal linear acoustics, the “Bulk Viscosity” can be even neglected (as it is done in Helmholtz Solver) by giving it an extremely small value like 1e-10 (perhaps even 0 would work). But a more realistic value of Bulk Viscosity = Real 5.0e-3 works fine as well.

In order to clear the confusion let me present here a simple validation setup, which produces correct results.

The geometry is a 2D wave guide 1m long, 0.1m high, filled with air. It uses a 40x4 second order mesh. The simulation is done in time domain using 200x1.4577e-5 s time steps. The input velocity is a sinusoidal time function of 1e-3 m/s amplitude and 343 Hz frequency. The input signal amplitude was chosen to be small in order to approximate a linear wave propagation, where the dispersion can be neglected at short distances.

Here is a snapshot of the acoustic velocity at time step 199:
velo_x.png
velo_x.png (13.75 KiB) Viewed 5738 times

And here is how the wave propagates in time and space:
step_50.png
step_50.png (6.92 KiB) Viewed 5738 times
time step 50


step_100.png
step_100.png (7.29 KiB) Viewed 5738 times
time step 100

Only 3 attachments are allowed per post, therefore this text will continue in the next post below…

Zoltan
Z.Losonc
Posts: 40
Joined: 21 Dec 2017, 16:05
Antispam: Yes

Re: Acoustics - Large-amplitude Wave Motion in Air

Post by Z.Losonc »

...continued from the previous post:
step_150.png
step_150.png (7.7 KiB) Viewed 5738 times
time step 150

step_199.png
step_199.png (8.16 KiB) Viewed 5738 times
time step 199

The last image shows the moment when the signal front reaches the 1m distance from the input. The time it took to get there is 199*1.4577e-5=0.002900823 s, which corresponds to a sound speed of 344.7 m/s. This is very close to the expected 343 m/s, therefore the solver passed this simple preliminary test of checking the sound speed.

Here is the sif file with the correct “Specific Heat” parameter value, a negligibly small “Bulk Viscosity”, and a small input signal amplitude to simulate linear wave propagation:

Code: Select all

Header
  CHECK KEYWORDS Warn
  Mesh DB "." "."
  Include Path ""
  Results Directory ""
End

Simulation
  Max Output Level = 10
  Coordinate System = "Cartesian 2D"
  Coordinate Mapping(3) = 1 2 3
  Steady State Max Iterations = 1
  BDF Order = 2

  Simulation Type = Transient
  Timestepping Method = BDF

  Timestep Intervals = 200
  Timestep Size = 1.4577e-5
  Output Intervals = 1

  Adaptive Timestepping = False
  Adaptive Time Error = 0.05
  Adaptive Keep Smallest = 3
  Adaptive Min Timestep = 1.0e-6
  Adaptive Max Timestep = 0.01

$w = 2155

  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) = 1
  Name = "Body Property 1"
  Equation = 1
  Material = 1
  Initial Condition = 1
End

Equation 1
  Name = "CompressibleNS"
  Active Solvers(1) = 1
End

Solver 1
  Exec Solver = String "always"
  Equation = "CompressibleNS"
  Procedure = "CompressibleNS" "CompressibleNS"
  Variable = -dofs 4 Velocity[Velo:2 Dens:1 Temp:1]
!  Element = p:1 b:1
   Element = p:2

  Stabilize = True
  Bubbles = False

  Steady State Convergence Tolerance = 1.0e-5

  Nonlinear System Convergence Tolerance = 1.0e-5
  Nonlinear System Max Iterations = 30
  Nonlinear System Newton After Iterations = 3
  Nonlinear System Newton After Tolerance = 1.0e-5
  Nonlinear System Relaxation Factor = 1
  Nonlinear System Convergence Measure = solution   

  Linear System Solver = Iterative
  Linear System Iterative Method = BiCGStab
  Linear System Max Iterations = 500
  Linear System Convergence Tolerance = 1.0e-5
  Linear System Preconditioning = ILU1
  Linear System ILUT Tolerance = 1.0e-4
  Linear System Abort Not Converged = False
  Linear System Residual Output = 10
  Linear System Precondition Recompute = 1
End

Material 1
  Name = "Air (room temperature)"
  Heat expansion Coefficient = Real 3.43e-3
  Heat Conductivity = Real 0.0257
  Equilibrium Density = Real 1.205
  Equilibrium Temperature = Real 293.0
  Specific Heat = Real 718 		!Don't mix J/(kg*K) with kJ/(kg*K)!
  Specific Heat Ratio = Real 1.404
  Viscosity = Real 1.81e-5
  Bulk Viscosity = Real 1.0e-10
!  Bulk Viscosity = Real 5.0e-3	!use this for nonlinear cases at low frequencies
End

Boundary Condition 1
  Target Boundaries(1) = 2
  Name = "Out"
  Velo 1 = Real 0.0
  Velo 2 = Real 0.0
End

Boundary Condition 2
  Target Boundaries(2) = 1 3
  Name = "Walls"
!  Velo 1 = Real 0.0
  Velo 2 = Real 0.0
End

Boundary Condition 3
  Target Boundaries(1) = 4
  Name = "In"
  Velo 1 = Variable Time
	Real MATC "1.0e-3*sin(w*tx)"
!  Velo 2 = Real 0.0
End

Initial Condition 1
  Name = "Initial Condition 1"
  Velo 1 = Real 0.0
  Velo 2 = Real 0.0
  Dens = Real 0.0
  Temp = Real 0.0
End

The project file with the mesh is also attached to this post, which can be used as a starting point for practical applications. Some more in-depth validation is required for strongly nonlinear cases, before the solver can be announced reliable. Hopefully this case will be of use to people in future.

Zoltan
Attachments
CompressibleNS1.7z
(30.17 KiB) Downloaded 288 times
mzenker
Posts: 1999
Joined: 07 Dec 2009, 11:49
Location: Germany

Re: Acoustics - Large-amplitude Wave Motion in Air

Post by mzenker »

Hi Zoltan,

glad to hear that you got reasonable results now. So there is still hope... ;)

Matthias
mika
Posts: 230
Joined: 15 Sep 2009, 07:44

Re: Acoustics - Large-amplitude Wave Motion in Air

Post by mika »

Hi,

It seems that the BCs specified in the sif file might be too slippery to make the temperature fluctuation unique (do-nothing BC implies that grad T.n = dT/dn = 0 is given throughout the boundary, so a solution to the heat equation can only be unique up to an additive constant). In addition to making the temperature unique, I would also define (change the order of "Dens" and "Temp")

Variable = -dofs 4 Velocity[Velo:2 Temp:1 Dens:1]

so that the temperature fluctuation is found as the variable "Temp".

Elmer sources come with a test case

elmerfem/fem/tests/DNS_WaveSimulation

where this solver is used to resolve also the viscous and thermal boundary layers.

-- Mika
cohor
Posts: 4
Joined: 08 Apr 2018, 17:28
Antispam: Yes

Re: Acoustics - Large-amplitude Wave Motion in Air

Post by cohor »

Hello,
Z.Losonc wrote:Some more in-depth validation is required for strongly nonlinear cases, before the solver can be announced reliable. Hopefully this case will be of use to people in future.
Thank you for starting this thread, it is useful. Did you test it for strongly non-linear cases?
mika wrote: Elmer sources come with a test case
elmerfem/fem/tests/DNS_WaveSimulation
Mika, probably there is an error in the models manual where it says:
p = RρT
R = (γ − 1)Cv
in the first equation R supposed to be the specific gas constant. But if that is so, then the second equation is wrong.
Can you explain that?

-cohor
raback
Site Admin
Posts: 4801
Joined: 22 Aug 2009, 11:57
Antispam: Yes
Location: Espoo, Finland
Contact:

Re: Acoustics - Large-amplitude Wave Motion in Air

Post by raback »

Hi Zoltan

If you don't need the thermal effects you could also try to use the standard FlowSolve with density (and its sensitivity) given as function of pressure. This recovers quite nicely the compressibility waves and should be a few times faster. There is an example on the proper material definitions in test case CompressIdealGas2.

-Peter
cohor
Posts: 4
Joined: 08 Apr 2018, 17:28
Antispam: Yes

Re: Acoustics - Large-amplitude Wave Motion in Air

Post by cohor »

Hello Raback,

Mika has written the Model 12 in the models manual, but you might also know the answer to these questions.
Do you agree that the equation R = (γ − 1)Cv is wrong?
Is this same equation used in the program (or is this an error only in the manual)?
If the program doesn't use this, then what equation it uses instead?

-cohor
Post Reply