## Acoustics – Basic Postprocessing With ElmerSolver & Paraview

Numerical methods and mathematical models of Elmer
Z.Losonc
Posts: 40
Joined: 21 Dec 2017, 16:05
Antispam: Yes

### Re: Acoustics – Basic Postprocessing With ElmerSolver & Paraview

Hi,

Thank you Peter for revealing this “special” syntax of MATC expression. There is nothing written about this in the MATCManual.pdf, and (unless analysing the source code) I would have perhaps never used this form. More accurately, I would have never used this form again, because I have attempted to use it already on the top of the sif file, outside of sections a while ago, and it did not work. In our example this form on the top of the sif file:

Density_M = \$ 1.205

stops the solver with the error message:

“ERROR:: Model Input: Unknown input field in header section: density_m”

The alternate form of:

Density_M = Real \$ 1.205

Doesn't work either. This experience lead me to conclude that this syntax that you have presented was incorrect, and the \$ sign must be always at the beginning of the instruction (which was a wrong conclusion). Without detailed explanations I am having hard time to grasp the logic of Elmer. I hope that you continue to reveal the occult (hidden) rules and techniques of practical Elmer implementation.

To summarize my present understanding of MATC syntax:

When using the short \$ form of expression:
- the \$ sign must be at the beginning of the expression when it is outside of sections like this:
\$ Density_M = 1.205
and the variable type (Real, Int, etc.) may not be used (“\$ Density_M = Real 1.205” won't work).
- within the sections the \$ sign must come after the = mark like this:
Density = \$ Density_M
in this case the variable type may (or in some cases must) also be used:
Density = Real \$ Density_M
- this short \$ form of expression will get evaluated only once, when the sif file is being read by Elmer. Therefore, if an expression needs to be evaluated only once, we can save processor time by using this economical short form instead of the MATC form.

When using the long MATC form of expression:
- this form may not be used outside of sections. Therefore this form at the top of the sif file would not work:
Density_M = Real MATC "1.205" neither this:
\$ Density_M = Real MATC "1.205"
- inside the sections the keyword MATC must come always after the = sign, and may be (or sometimes must be) preceded with the variable type specifier.
- the mathematical expression must be between two “ marks, and after the keyword MATC:
Cv = Real MATC "1/(Density_M * Angular_Frequency_M)"
- this expression form will be evaluated (only if and) every time when the parameter it was assigned to (Cv in this case) is called/used somewhere during the computation. If it is never called, then the expression will never get evaluated; and if it contains errors, these errors will never get reported.
- this feature of being re-evaluated every time the parameter is used consumes more processor time, but it enables us to change its value during the simulation, which can be useful f. ex. during parametric sweeps.
- here is an example (from FEM-CurvedFlowinaPipe.pdf) how to specify a variable, and assign a value to it with a MATC expression:

Code: Select all

``````Material 1
Name = "Water (room temperature)"
Viscosity = Variable Temperature
Real MATC "mu0 * relativevisc(tx)"``````
In this form the tx represents the parameter Temperature as an argument of the function relativevisc. The function can be defined as:

Code: Select all

``````\$ mu0 = 1.788e-3
\$ function relativevisc(T){\
a = -1.704;\
b = -5.306;\
c = 7.003;\
z = 273.15/T;\
_relativevisc = exp(a + b * z + c *(z^2));\
}``````
If there is something wrong with these explanations (or need to be extended), please let me know. Based on the simpler \$ form of MATC expressions, here is an alternative sif for our sample project that works as well:

Code: Select all

``````\$ Density_M = 1.205
\$ Angular_Frequency_M = 2155

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

Simulation
Max Output Level = 10
Coordinate System = Cartesian
Coordinate Mapping(3) = 1 2 3
Steady State Max Iterations = 1
Output Intervals = 1
Timestepping Method = BDF
BDF Order = 1
Solver Input File = tube.sif
Post File = tube.ep
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

Solver 1
Equation = Helmholtz Equation
Procedure = "HelmholtzSolve" "HelmholtzSolver"
Variable = -dofs 2 Pressure Wave
Stabilize = True
Bubbles = False
Lumped Mass Matrix = False
Optimize Bandwidth = True
Steady State Convergence Tolerance = 1.0e-5
Nonlinear System Convergence Tolerance = 1.0e-8
Nonlinear System Max Iterations = 20
Nonlinear System Newton After Iterations = 3
Nonlinear System Newton After Tolerance = 1.0e-3
Nonlinear System Relaxation Factor = 1
Linear System Solver = Direct
Linear System Direct Method = Banded
End

Solver 2
Equation = "flux compute 1"
Procedure = "FluxSolver" "FluxSolver"
Calculate Flux = Logical True
Target Variable = String "Pressure Wave 2"
Flux Coefficient = String "Cv"
Linear System Solver = Direct
Linear System Direct Method = Banded
End

Solver 3
Equation = "flux compute 2"
Procedure = "FluxSolver" "FluxSolver"
Calculate Flux = Logical True
Target Variable = String "Pressure Wave 1"
Flux Coefficient = String "Cv"
Linear System Solver = Direct
Linear System Direct Method = Banded
End

Solver 4
Equation = Result Output
Procedure = "ResultOutputSolve" "ResultOutputSolver"
Save Geometry Ids = False
Output File Name = tube
Output Format = Vtu
Scalar Field 2 = Pressure wave 2
Scalar Field 1 = Pressure wave 1
Vector Field 1 = Pressure wave 2 flux
Vector Field 2 = Pressure wave 1 flux
End

Equation 1
Name = "Helmholtz"
Angular Frequency =  \$ Angular_Frequency_M
Active Solvers(4) = 1 2 3 4
End

Material 1
Name = "Air (room temperature)"
Viscosity = 1.983e-5
Heat expansion Coefficient = 3.43e-3
Heat Conductivity = 0.0257
Relative Permittivity = 1.00059
Sound speed = 343.0
Heat Capacity = 1005.0
Density = \$ Density_M
Cv = Real \$ 1/(Density_M * Angular_Frequency_M)
End

Initial Condition 1
Name = "InitialCondition 1"
Pressure Wave 1 = 0
Pressure Wave 2 = 0
End

Boundary Condition 1
Target Boundaries(1) = 6
Name = "In"
Pressure Wave 2 = 0
Pressure Wave 1 = 100
End

Boundary Condition 2
Target Boundaries(4) = 2 3 4 5
Name = "Walls"
Wave Flux 1 = 0
Wave Flux 2 = 0
End

Boundary Condition 3
Target Boundaries(1) = 1
Name = "Out"
Wave impedance 1 = 343
Wave impedance 2 = 0
End``````
Good to know that Elmer can do parametric sweeps when the dimensions of the geometry are swept. I will study the recommended RigidMeshMapper and MeshUpdate later, when the more basic features are understood.

In my earlier post dated 06 Jan 2018, 03:37 I have discussed that the default iterative solver of the FluxSolver, which can be found in the ElmerModelsManual.pdf did not work for me:

Code: Select all

``````Solver 3
Exec Solver = after all
Equation = "flux compute"
Procedure = "FluxSolver" "FluxSolver"
Calculate Flux = Logical True
Flux Variable = String Temperature
Flux Coefficient = String "Heat Conductivity"
Linear System Solver = "Iterative"
Linear System Iterative Method = "cg"
Linear System Preconditioning = ILU0
Linear System Residual Output = 10
Linear System Max Iterations = Integer 500
Linear System Convergence Tolerance = 1.0e-10
End
``````
Is there a fix for this problem?

The alternative BiCGStab did produce correct results, but with warnings in the log file about convergence not being reached:

Code: Select all

`````` 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 = 20
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
``````
This problem looks very similar to the one reported by adamo in this thread viewtopic.php?f=3&t=4930&start=0. It looks like the solver residua are too small, and the solver can not handle it. Is there a way to overcome these problems?

Which iterative solver would be the best to use in FluxSolver for acoustic simulations?

Zoltan

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

### Re: Acoustics – Basic Postprocessing With ElmerSolver & Paraview

Hi,

The questions in this post are primarily to the Elmer team, but if anybody else also knows the correct answers, then please comment.

The next step in post-processing acoustic simulations with ElmerSolver would be to calculate the acoustic (sound) Intensity (RMS), which can be calculated from the variables we have derived so far using the equation:

I = (pressure wave 1 * pressure wave 2 flux - pressure wave 2 * pressure wave 1 flux)/2

The MATC does not seem to be able to use the parameters and variables (used by and) computed by other solvers, therefore the intensity can't be calculated using MATC. Is this correct?

Is there any other way to compute sound intensity in ElmerSolver?

By the way, the post processor can't be launched from the ElmerGUI. Was this deliberately disabled in Elmer 8.3, or just an installation error?

Zoltan

mzenker
Posts: 1954
Joined: 07 Dec 2009, 11:49
Location: Germany

### Re: Acoustics – Basic Postprocessing With ElmerSolver & Paraview

Hi,

it is recommended to use ParaView as postprocessor. It has many features allowing you to calculate whatever you want provided that it can be done using the results calculated by Elmer.
The Elmer postprocessors are not under active development any more AFAIK.
To use ParaView, simply replace the .ep extension of the output file by .vtu in the simulation section (Model->Setup in ElmerGUI).
You can also calculate things in Elmer using the SaveScalars solver, see models manual.

HTH,

Matthias

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

### Re: Acoustics – Basic Postprocessing With ElmerSolver & Paraview

Hi,

I am using Paraview, as it has been already discussed in this thread, and the vtu file is written in every sif I posted. I don't mind if the Elmer Postprocessor is not included in the installation; it is just useful to know whether it was deliberately left out from the installation binary, or it can't be started due to an error. If it was deliberately removed, then it would be more user friendly to remove the starting button and menu option from ElmerGUI as well.

The Elmer VTK is very useful though for having a quick look at the results when experimenting with different settings in the sif file. Please don't remove that from ElmerGUI, not even if Paraview is better for post processing.

As described earlier, the computing of acoustic velocity in Elmer (instead of doing it in Paraview) is necessary at this time (due to Paraview filter bug), and we have done that successfully. The computation of acoustic intensity on the other hand, can be done in Paraview. The only reason I am trying to do it in Elmer is to learn whether it can be done at all, or not; and what are the exact capabilities of Elmer. We were discussing the practical usage of MATC, and it would make sense to see if MATC is able to multiply a vector field computed by the FluxSolver with a scalar filed, which is the output of the main Helmholtz solver. I was reading the manuals, and read all acoustic related threads in this forum, but could not find anything about this.

There is a way to calculate a scalar field (f. ex. the absolute pressure in dB) from the variables of the main solver like this:

Code: Select all

``````Body 1
…
Body Force = 1
End

Body Force 1
! This is p_abs in dB
Press_dB = Variable Pressure Wave 1, Pressure Wave 2
Real MATC "20*log(((sqrt(tx(0)^2+tx(1)^2))/sqrt(2)))"
End

Solver 1
Equation = Helmholtz Equation
Procedure = "HelmholtzSolve" "HelmholtzSolver"
...
Nonlinear Update Exported Variables = Logical True
Update Exported Variables = Logical True
Exported Variable 1 = Press_dB
End

Solver 4
Equation = Result Output
Procedure = "ResultOutputSolve" "ResultOutputSolver"
...
Scalar Field 3 = "Press_dB"
End``````
See: viewtopic.php?f=4&t=2120&p=6288&hilit=a ... 8c15#p6288

But this worked for me only as long as the result is a scalar field. If the result supposed to be a vector field computed from another vector field, I could not get it to work. But there might be some hidden syntax, or trick that we are not aware of, and perhaps the Elmer developers would be kind enough to share such knowledge with us… If it is not possible do such computation, no problem. Just please, let us know explicitly what we can do, and what we can't do in Elmer.

The SaveScalars solver may be used to compute derived quantities and save scalar values to an external file. But in our case, we have to compute and save vector values (not scalar values) to the output file, therefore this won't help, as far as I can tell.

Zoltan

CWeng
Posts: 4
Joined: 04 Sep 2019, 18:48
Antispam: Yes

### Re: Acoustics – Basic Postprocessing With ElmerSolver & Paraview

Hello,
I am new to Elmer. In the beginning I found it hard to set up acoustic simulation for multiple frequencies and to export boundary data for post processing. Thanks to your great posts, I finally got what I need! Here I share one of my post processing code for the calculation of transmission loss (TL) of acoustic two ports by using only the inlet and outlet data; hope it is helpful for the community.
Summary of the .sif file (enclosed):
• A test case for computing the transmission loss of a simple area expansion (See Bild1.png) is studied.
• Simulation type is "scanning" so we can run the simulation for multiple frequencies.
• The grad. of pressure is computed, based on what I learnt from the previous posts in this topic.
• Pressure and its gradient on the inlet and outlet surface are saved for post processing. For the computation of TL we don't need to save all the data from the computational domain. If the mesh size is large, saving all volume data for multiple frequencies may be storage demanding.
Bild1.png
The data are post processed using numpy. The python file is enclosed. Lines 38-44 may be modified to use the code for other simulation cases.

The computed TL is compared with the low-frequency theory ( https://en.wikipedia.org/wiki/Transmiss ... acoustics) ); see validation.png.
validation.png (56.9 KiB) Viewed 75 times
I have to say that the post-processing code is not well tested yet. It is simply an example showing how to deal with acoustic data. Please tell me if there is any mistake in the code.
Limitations of the current post-processing code:
1, The inlet and outlet are with same cross-section area. However, the code can be easily modified to deal with different inlet/outlet pipe sizes, by simply including the area into the calculation of the TL.
2, The "planeWaveOutlet" boundary condition is assumed to be "non-reflecting". However, this boundary condition may fail in preventing reflection in some cases (e.g., when higher order modes exist or when duct wall impedance is not symmetric). In these cases the "Two source method" (https://www.researchgate.net/publicatio ... _Acoustics) may be used to obtain the scattering matrix.
3, Now the sound wave at the inlet is decomposed into upstream and downstream propagating wave using the pressure and velocity by assuming plane-wave propagation. If higher order acoustic modes exist, we may need to use the multiple-microphone method to decompose the sound. I will investigate this soon (I think it is possible for Elmer to position some virtual microphone in the computational domain).

C.Weng

Code: Select all

``````# -*- coding: utf-8 -*-
# Computing transmission loss from 2-port simulation.
# 1st version: 11/09/2019 at Tuve, CWeng
import numpy as np
import matplotlib.pyplot as plt

#************ Parse *.dat.names  #************#************
# Should return res as a Str list which is used
# later for mapping the nemerical data in the .dat files...
def parseDatNames(filename):
res = [] # Initialize
if filename[-6:] != ".names":
filename = filename + ".names"
with open(filename, 'r') as f:
for line in f:
if 'Variables in columns' in line:
continue
continue
# finally we can do some work here:
# strParse1: get rid of 1: boundary mean, and split over "over"
strParse1 = line.strip().split(':')[-1].strip().split('over')
if len(strParse1) == 1:
res.append(strParse1[0])
else:
if len(strParse1[0]) > 0:
itemName = strParse1[0]
else:
strParse1[0] = itemName
res.append(strParse1[0].strip() + '.' + strParse1[1].strip())
return res
#************#************#************#************

## User input:
filename = 'BC0D.dat'
inletBCid = 1
outletBCid = 2
rho0 = 1.205
c0 = 343
flowDir = np.array([1, 0, 0]);

# Parse input names:
resNames = parseDatNames(filename)

with open(filename, 'r') as f:
res = np.array([np.fromstring(line.strip(), dtype=float, sep=' ') for line in f])

#############################     Assign values #########################
## Get freq.
freqVec = res[:, resNames.index('frequency')]
omegaVec = 2*np.pi*freqVec
## Get pressure at inlet and outlet
pIn = res[:, resNames.index('pressure wave 1.bc ' + str(inletBCid))] \
+ 1j*res[:, resNames.index('pressure wave 2.bc ' + str(inletBCid))]
pOut = res[:, resNames.index('pressure wave 1.bc ' + str(outletBCid))] \
+ 1j*res[:, resNames.index('pressure wave 2.bc ' + str(outletBCid))]

## Get velocity at inlet. Note that uIn = 1j/rho0/omegaVec*dp/dn
# To decompose the p to p+ and p-, we notice that
#     p+ + p- = p
#     p+ - p- = uIn*rho0*c0
# so, p+ = (p + uIn*rho0*c0)/2

#  Check if it is 3D
dim = 3
if "pressure wave 1 grad 3.bc 1" not in resNames:
dim = 2
flowDir = flowDir[:-1]

fmt = 'pressure wave {0:d} grad {1:d}.bc {2:s}'
for ind in range(dim):
# Inlet
strReal = fmt.format(1, ind+1, str(inletBCid))
strImag = fmt.format(2, ind+1, str(inletBCid))
pGradIn[:, ind] = res[:, resNames.index(strReal)] \
+ 1j*res[:, resNames.index(strImag)]
# outlet
strReal = fmt.format(1, ind+1, str(outletBCid))
strImag = fmt.format(2, ind+1, str(outletBCid))
pGradOut[:, ind] = res[:, resNames.index(strReal)] \
+ 1j*res[:, resNames.index(strImag)]

# Decomposition at inlet
uIn = 1j/rho0/omegaVec*dpdnIn
pInPlus = (pIn + uIn*rho0*c0)/2
pInMinus = pIn - pInPlus

# Decomposition at outlet
uOut = 1j/rho0/omegaVec*dpdnOut
# %% TL, assuming same cross section area in inlet and outlet

TL = 20*np.log10(np.abs(pInPlus)/np.abs(pOut))

plt.close("all")
fig, ax = plt.subplots()
ax.minorticks_on()
#ax.grid(which='both')
ln = plt.plot(freqVec, TL, '-o', label='Simulation')
plt.xlabel('Frequncy, Hz'), plt.ylabel('TL, dB')
plt.title(filename)
plt.show()

``````
Attachments
areaExpansion.sif
Last edited by CWeng on 05 Oct 2019, 22:40, edited 1 time in total.
Best regards,
C. Weng

raback
Posts: 3461
Joined: 22 Aug 2009, 11:57
Antispam: Yes
Location: Espoo, Finland
Contact:

### Re: Acoustics – Basic Postprocessing With ElmerSolver & Paraview

Hi CWeng

Nice work!

Would you have some development ideas for the Helmholtz solver that would make the code more usable in this application area? Often there are some specific details in each field that may be pretty easy to implement in the code.

-Peter

CWeng
Posts: 4
Joined: 04 Sep 2019, 18:48
Antispam: Yes

### Re: Acoustics – Basic Postprocessing With ElmerSolver & Paraview

Hello Peter,
Thank you for the feedback, and I do appreciate you and your team members for the efforts you have put on Elmer. Regarding your question on the development ideas, I will provide my thoughts based on my limited knowledge and experiences.
If we only talk about extending the current Helmholtz solver, instead of including new wave equations (such as the linearized Euler's equation), then I currently can come up with three ideas. I try to order them based on their ease of implementation (from easy to difficult):

1, Including temperature gradient in the wave propagation.
Based on Eq. (11.2) in the Elmer Model Manual, the current Helmholtz equation is based on the assumption of constant background temperature, and hence a constant speed of sound. In some applications, however, the temperature gradient is large and the resulting spatial variation of speed of sound is not negligible. For example, in automobile exhaust systems the variation of the temperature from the hot-end to the tail pipe exist may reach hundreds of degree. To estimate acoustic resonance of the exhaust systems accurately, one may have to consider the variation of the temperature, and hence the speed of sound, in his acoustic model.
To include the variation of speed of sound in Elmer's model, Eq.(2.50) in the following book may be used:
https://www.win.tue.nl/~sjoerdr/papers/boek.pdf
This equation can be easily transformed into the frequency domain. The difficulty might be that, Elmer should allow the user to provide the non-constant speed of sound c0 as a background quantities, either as numerical values of as analytical expressions. Usually c0 is obtained numerically by solving the steady heat equation and using the ideal gas law c0= sqrt(gamma*R*T). If the user choose to use the same mesh in both the Helmholtz model and the heat equation model, then Elmer may possibly read the values of c0 directly (similar to the "Convection Computed" in the Heat Equation model, maybe?); otherwise Elmer may implement certain interpolation scheme while reading c0 to the Helmholtz solver from a different mesh.

2, Including background potential flow in the wave equation.
The currently Helmholtz model in Elmer assumes quiescent background flow. However, in many applications the acoustic wave is propagating in moving medium and the propagation may therefore be influenced by the convection and refraction effects of the mean flow.
The Helmholtz equation can be extended to include background mean flow, in order to take into account the convection effect (but perhaps not the refraction effect, unfortunately). This may be done by deriving the model from the last equation on page 3 in the following document,
https://arxiv.org/pdf/1009.3766.pdf

ie.,
ConvectiveWaveEquation.PNG (5.96 KiB) Viewed 53 times
where phi is the acoustic potential. The background mean flow in this equation is potential flow, which can be obtained by Elmer potential flow solver. It should be noted that, if this equation is implemented, then the variation of the speed of sound (denoted as a0 in this equation) can also be considered.
One application of this extended Helmholtz equation is to reversely compute the impedance of aero-engine acoustic liners from experimental data, i.e., the so-called "impedance eduction" which has been investigated exhaustively by NASA, DLR, ONERA and many other institutes and universities. As an example, an early paper from NASA,
https://ntrs.nasa.gov/archive/nasa/casi ... 090459.pdf
uses the convective Helmholtz equation with uniform background flow, Eq. (4), and the Myers boundary condition Eq. (10) to educe the impedance of the linear. One may see the effect of mean flow in Fig. 7. More recent research on the impedance eduction can be found from the google citation page of the paper shown here:
The difficulty of implementing this model may be similar to that in the previous idea. In addition, changing the variable from pressure to potential may require changes of the boundary-condition formulation that has been already implemented in Elmer.

3, Computational aeroacoustics (CAA)
As I know, the computation of aerodynamically generated sound remains a hot topic in the acoustic community. One may read the following Wikipedia pages for a brief introduction of the topic.
https://en.wikipedia.org/wiki/Aeroacoustics
https://en.wikipedia.org/wiki/Computati ... oacoustics
In summary, the CAA is about determining the sound source from the turbulent (sometimes laminar) flow, and propagate the sound generated by the source. The Helmholtz equation in Elmer is a model that can propagate the sound, but it is not a model that can estimate the aeroacoustic source (i.e., the right hand side).
I think it would be a good news for the CAA community if Elmer can extend its model to include aeroacoustic sources. The Lighthill analogy that is introduced in the Wiki page might be a good model to start with. Usually the integral solution of the Lighthill analogy is limited to free-space radiation, e.g., no solid boundaries are allowed while using the model though the boundaries usually "enhance" sound generation at low Mach numbers. However, if Elmer implements the FEM version of the analogy, then solid boundaries may be taken into account. Details can be found in this paper:
https://www.sciencedirect.com/science/a ... 2500002061
where the variational form of the analogy is derived. In the implementation Elmer may allow the user to input the flow quantities (E.g., the Lighthill stress tensor Tij in Eq. 19 in the paper) obtained from CFD.
Last edited by CWeng on 05 Oct 2019, 23:05, edited 2 times in total.
Best regards,
C. Weng

CWeng
Posts: 4
Joined: 04 Sep 2019, 18:48
Antispam: Yes

### Re: Acoustics – Basic Postprocessing With ElmerSolver & Paraview

Hello Peter,
Another useful feature is the non-reflection condition (maybe it already has been implemented but I haven't read it yet; I am still learning Elmer). If I am correct, in the current version of Elmer the non-reflection boundary condition is realised by setting the boundary impedance to the characteristic impedance (density*speed of sound). This may work well for normal incident but not for oblique incident, and the latter is encountered in many applications. For example, in duct acoustics if the frequency is high enough or if the duct wall impedance are not symmetric, etc, non-plane waves are able to propagate. Then oblique incident may occur at the boundaries and the wave is reflected. Another example is sound radiation in free space. If the user want to compute the radiated acoustic power he has to make sure there is no reflection from the outer boundaries which is generally hard to achieve by using a non reflection boundary condition (Actran uses the infinite element and the perfectly matched layer to realise the radiation condition, see https://www.fft.be/product/actran-acoustics)

To create more general non-reflecting condition, I think the buffer-zone approach might be adopted. This may be realised using the perfectly matched layer (PML). I find the following introduction to PML very helpful
https://math.mit.edu/~stevenj/18.369/pml.pdf
and figure 1 in the above tutorial explains the idea behind the buffer zone clearly
pml_fromStevenJ.PNG (117.95 KiB) Viewed 50 times
Since the Helmholtz model is a frequency solver so the implementation of PML might be straightforward.

Another buffer-zone approach might be the implementation the artificial viscosity. In the Helmholtz solver the artificial viscosity may be given as the imaginary part of the speed of sound. I also saw that in the WaveSolver the viscosity has already been implemented and can be given as parameters:
artificialViscosity.PNG (8.4 KiB) Viewed 50 times
In the region of interest (ROI) the viscosity can be set to zero. In the buffer zone the viscosity can be set to an nonphysically large value to damp the wave quickly within the zone. To avoid reflection due to sudden change in the viscosity at the interface between the ROI and buffer zone, the viscosity may be turned on gradually from zero to the large value at the outer boundary (linearly or quadratically, etc..). I think this function would be easier to implement if Elmer simply allows the users to define the viscosity as a function of space, i.e., eta = eta(x,y,z); then, the users have the freedom to create the buffer zone and tune parameters such as the width of the zone and the turn-on function of the viscosity.
Best regards,
C. Weng

raback
Posts: 3461
Joined: 22 Aug 2009, 11:57
Antispam: Yes
Location: Espoo, Finland
Contact:

### Re: Acoustics – Basic Postprocessing With ElmerSolver & Paraview

Hi C. Weng,

Great development ideas! Some of them would be quick to implement, others might take a little longer.

Conserning idea 1) Note that even the current implementation allows for density to be a function of something. You could test it by simply setting

Code: Select all

``````  Density = Variable "Temperature"
Real MATC "1.29/(1+tx/273.15)"
``````
where Temperature would be the precomputed temperature of air using the ideal gas law. MATC is notoriously slow, so for real use better use piece of code.

Would there be interest to test and demonstrate some new features the coding work might not be a bottle-neck.

-Peter