Acoustic Boundary Element Solver.
Dr Dan J. O’Boy
http://www.djoboy.co.uk
R1
All rights reserved.
This manual is provided for reference and may not be copied and / or altered without the permission of the author. References to the manual should be made whenever material is used.
The manual provides guidance and information in order to use the Sonitus BE program and may be altered without notice being provided. No responsibility is assumed or given for any damage, commitments, liabilities or inaccuracies incurred as a result of the information provided.
All material is the intellectual property of this company through the owner Dr Dan J. O'Boy, who applies the property rights. The terms and conditions are applied through and in accordance with English law and any disputes will be subject to the jurisdiction of the courts of England and Wales of the United Kingdom.

Revision 1
Copyright ©2008
Dr Dan J O’Boy
4 The Green
Exton
Oakham
Rutland
LE158AP
United Kingdom
enquires@djoboy.co.uk
The input files can be generated automatically using the software Sonitus PreBE, which formats the geometry and composes an appropriate structured mesh. As an alternative, the Matlab code for defined geometry can be utilised for the mesh generation.
Figure 12 Generation of a surface mesh using computer aided design tools
In order that users may determine their own software for input file generation, sample input files are provided showing the order of the file structure. The input data file (.dat) has the parameters described in Table 6.
Description of file information. 
File heading, including creation time and date. 
A list of fixed frequencies the solver will excite the domain with. 
Information lines [4]. 
Solution type (indirect or direct). 
Tolerances for the Sonitus PreBE geometry creator. 
Reference values for the pressure, intensity and work. 
Nodal numbers and positions (x,y,z). 
Number and coordinates of the data recovery points (x,y,z). 
Details of how the nodes join together to create elements (triangular or quadrilateral). 
Details of how the nodes of the data recovery mesh join together to create elements. 
Material description: what is the domain made of, density, speed of sound, real and imaginary terms are applicable here in the case of damping being studied. 
Boundary conditions, velocities or velocity slopes applied to each structural element. 
Monopole source location and strength. 
Display information on the elements. 
Sonitus PreBE display standards. 
Input file end line. 

Table 6 Input file description
The results file is generated automatically using the software Sonitus PostBE, which loads the results file (*.rslt), and displays the data on the screen. A Matlab code is also available for this purpose, in order to allow integration with other CAE programs.
In order that users may determine their own results analysis software, sample results files are provided showing the order of the file structure. The output data file (.rslt) has the information provided in Table 7.
Description of file information. 
File heading including generation date and time. 
Number of frequencies which have been analysed. 
A list of the frequencies together with the solution type (direct or indirect). 
For each frequency, a list of the following parameters are provided for every data recovery point. 
Sound Pressure Level (dB). 
Acoustic Pressure, Magnitude. 
Acoustic Pressure, Phase Angle. 
XComponent of Velocity, Magnitude. 
XComponent of Velocity, Phase Angle. 
YComponent of Velocity, Magnitude. 
YComponent of Velocity, Phase Angle. 
ZComponent of Velocity, Magnitude. 
ZComponent of Velocity, Phase Angle. 
Magnitude of Velocity. 
XComponent of Intensity. 
YComponent of Intensity. 
ZComponent of Intensity. 
Output file end line. 

Table 7 Output file description
Syntax is SonitusBE INPUTFILE OUTPUTFILE
The Sonitus BE program code is comprised of a main file (SonitusBE.cpp) and a series of header files supplying generic subroutines. The pseudocode program listing is as follows:
Display welcome message to the screen.
Ensure correct number of arguments are passed at the command line.
Argument one is the input filename, two is the output filename, three is an unused protocol identifier.
Obtain input and output filenames from the command line arguments.
Call the routine to read in all input data, positions, boundary conditions and source locations.
Call the routine to check the status of the input mesh for consistency.
Call the routine to check that the required number of boundary conditions have been supplied.
Generate an empty list of output result variables.
Generate the output file (.rslt format).
Call the routine holding the numerical solver (direct or indirect solver).
Repeat for each excitation frequency.
Print output messages to the screen.
Release memory and close program safely.
The primary namespaces for the boundary element solver are SonitusBem and BemRoutine
The subroutine header files contain a series of both generic and specific functions which manipulate the data in different ways. In this section, the individual routines are discussed and function calls detailed.
This file contains a routine which manages the error messages sent to it, and uses a built in function to print them to the screen. 
BemRoutine::GeneralException(const string &mess){message=mess;}
This file contains the routines which interpolate data from numerical recipes type functions. It exists to contain all of the routines which allow interpolation on triangular and quadrilateral elements. 
The primary namespace for the routines are Recipes, as the functions are derived both from the Numerical Recipes formulations, and formulae provided in Abramowitz and Stegun eds. (1972), Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, New York: Dover Publications, ISBN 9780486612720.
The following routines are provided based on an interpolation of double precision numbers.
· Interpolation routine for a triangular element with nodes _{} to _{} and nondimensionalised surface coordinates _{} and _{}, which vary from 0 to 1.
double interpTRI( const double &N1,
const double &N2,
const double &N3,
const double &zeta1,
const double &zeta2);
Define a set of weighting functions, _{}, _{}, _{}, equated to _{}, _{}, _{}.
Interpolated return value = _{}.
· Interpolation routine for a triangular element with nodes _{} to _{} and nondimensionalised surface coordinates _{} and _{} which vary from 0 to 1. This provides the derivative _{} at the values of _{} and _{}.
double DZETA1interpTRI(
const double &N1,
const double &N2,
const double &N3,
const double &zeta1,
const double &zeta2);
Return _{}= N1N3.
· Interpolation routine for a triangular element with nodes _{} to _{} and nondimensionalised surface coordinates _{} and _{} which vary from 0 to 1. This provides the derivative _{} at the values of _{} and _{}.
double DZETA2interpTRI(
const double &N1,
const double &N2,
const double &N3,
const double &zeta1,
const double &zeta2);
Return _{}= N2N3.
· Interpolation routine for a quadrilateral element with nodes _{} to _{} and nondimensionalised surface coordinates _{} and _{} which vary from 1 to 1.
double interpQUAD(const double &N1,
const double &N2,
const double &N3,
const double &N4,
const double &zeta1,
const double &zeta2);
Define a set of weighting functions, _{}, _{}, _{}, _{} equated to _{}, _{}, _{}, _{}.
Interpolated return value = _{}.
· Interpolation routine for a quadrilateral element with nodes _{} to _{} and nondimensionalised surface coordinates _{} and _{} which vary from 1 to 1. This provides the derivative _{} at the values of _{} and _{}.
double DZETA1interpQUAD(
const double &N1,
const double &N2,
const double &N3,
const double &N4,
const double &zeta1,
const double &zeta2);
Define a set of weighting functions, _{}, _{}, _{}, _{} equated to _{}, _{}, _{}, _{}.
Interpolated return value = _{}.
· Interpolation routine for a quadrilateral element with nodes _{} to _{} and nondimensionalised surface coordinates _{} and _{} which vary from 1 to 1. This provides the derivative _{} at the values of _{} and _{}.
double DZETA2interpQUAD(
const double &N1,
const double &N2,
const double &N3,
const double &N4,
const double &zeta1,
const double &zeta2);
Define a set of weighting functions, _{}, _{}, _{}, _{} equated to _{}, _{}, _{}, _{}.
Interpolated return value = _{}.
A similar set of routines are also provided for the case where the interpolated values are complex and of double precision.
ComplexMath_std::complex <double> CinterpTRI(
const ComplexMath_std::complex <double> &N1,
const ComplexMath_std::complex <double> &N2,
const ComplexMath_std::complex <double> &N3,
const double &zeta1,
const double &zeta2);
ComplexMath_std::complex <double> CinterpQUAD(
const ComplexMath_std::complex < double > &N1,
const ComplexMath_std::complex < double > &N2,
const ComplexMath_std::complex < double > &N3,
const ComplexMath_std::complex < double > &N4,
const double &zeta1,
const double &zeta2);
The primary namespace for this set of subroutine functions is SontitusBem.
This file contains the routines which print messages of welcome to the screen and to an output file. 
There are three functions provided which print information to the control console screen, showing the requirements of the program, the welcome message with copyright information and a final screen which lets the user know that the program has finished running for all frequencies.
void requirements(void);
void welcome(void);
void goodbye(void);
The primary namespace for this routine is BemRoutine.
This file returns the method used in the boundary element calculation. Only the indirect method is implemented in this code. 
The only function call is to determine what solver is used, either direct formulation or indirect. At this time, the only function only returns the indirect solver choice, regardless of the input conditions. If returned integer is 1, the direct solver would be chosen, if 2 the solver is indirect formulation.
int SolutionMethod(void);
The primary namespace for this routine is SonitusBemFormat.
This file checks the meshes for consistency. 
The function checkmesh provides basic error checking for the input structural mesh, in order to ensure that the correct number of nodes, elements, boundary conditions are provided, as without these the solver would be unable to reliably complete the calculation and would perform in an unpredictable fashion.
void checkmesh(const int & NumberofNodes,
const vector <int> & NodeNum,
const vector <MeshTypes::Node <double> > &NodeCoord,
const vector <MeshTypes::Element <int> > &ElemStruct,
const vector <int> &ElemNum,
const int &NumofElems)
This routine checks for repeated node numbers in the same element. It also compares sequential node numbers in different elements to ascertain whether the nodes are all pointing in consistent directions. Essentially, the element direction is determined by the order in which the nodes are arranges, either clockwise or counterclockwise. If two different elements have two of the same nodes in the same direction, then the two elements are facing in opposite directions.
This routine contains a series of if then else statements which test the above condition, whilst taking into account the possibility of triangular or quadrilateral element types.
void checkBCnumbers(const int &NumofElems,
const int &BCnum)
This routine compares the number of elements and the number of boundary conditions available. If these are the same, then we have a fully constrained problem and can send the problem to the boundary element solver for completion. If not, these we have a badly posed input problem and must correct the input files before proceeding.
The primary namespace for these subroutines is MeshTypes.
Routine allows vectors to be manipulated. 
This routine contains the definitions of different classes used in the main program.
Definition of a node, with three double elements representing _{}.
Definition of an element, defined by five fields (integers). The first field is a flag which says what element type is described.
If the flag is 1, this means a triangular element, and the following three fields are integers containing the node numbers defining that triangle.
If the flag is 2, this means a quadrilateral element, and the following four fields are integers containing the node numbers defining the corners of the quadrilateral.
Definition of a boundary condition number class. This has two entries, the first is the boundary element number (this should coincide with the element number) and the second is the type (1 is velocity and 2 is impedance).
The primary namespace for this routine is SonitusBemFormat.
This file contains the routines which read in the data from an input file. 
The routine call is as follows. Input the filename as a string and it will find the entries for the input data file.
void FilereaderSonitus::getinputdata(const string filename)
The routine by checking whether the input file can actually be opened or whether the routine needs to terminate an report a failure to complete. If no errors are reported, it means we can then generate and initialise variables.
The first seven lines are assumed to be organised header information describing the kind of input geometry, so that the user can keep track of what the file is representing.
Note that the maximum length of any input data line is 200 characters. Any lines exceeding this will cause a termination error.
The next lines should start with the word ACOUSTIC, then listing the acoustic frequency to investigate using the boundary element solver. Then read in and interpret control flags and data control cards.
Next the reference pressure (PREF), reference intensity (IREF) and reference work (WREF) are read in. The next four lines are then ignored.
The nodal coordinates are then read in, with each line starting with the word NODE. The node number is found, then the three coordinates interpreted. The next four lines are ignored.
The number and location of the data recovery nodes are then read in, with each line starting with DRNODE. The next four lines are ignored.
Once the nodes are known and organised, the elements can be defined. A triangular element is defined through three nodes, while a quadrilateral requires four nodes.
The data entry for a triangular element starts with the word TRIA3. If the triangular elements are not found, check for the word QUAD4, and read in the nodes which create these. The next four lines are ignored.
Although the data recovery nodes do not need to be generated as a mesh, it is found that this is more useful for processing results over an area. Given this is a possibility, the program now looks for the elements which make up the data recovery surface. Again, for each element, the control word TRIA3 or QUAD4 defines the type of element. The next two lines are ignored. At this point, the geometric information has been completely read in. The next information which is required is the type of material that the acoustic waves are propagating through, whether it is air or liquid. The primary difference is the density, with the control line starting with the word MATERIAL. The density has real and imaginary components, as does the subsequent speed of sound. The following six lines are then ignored.
The boundary condition information is then provided, such as whether for each element the boundary condition is a velocity (control word EVELO) or an impedance (EIMPE) and what value it takes.
Finally, in order to obtain a nontrivial solution, where no sound field is present, a source must be introduced into the calculation. At present, only one sound source can be introduced, although this may be modified should the requirement exist. The control word for this is SOURCE, followed by the source strength, including real and imaginary components, and the source location.
The routine ends by closing the input file and returning to the main program with the input data.
The primary namespace is BemRoutine.
Outputs the file results including sound pressure over the data recovery mesh. 
Once the boundary element solver has completed all of the calculations, it needs to output the results to a file which can be postprocessed at a later date without the worry of losing information held in memory. Rather than using a proprietary format for the generation and storage of results, the output files are held in a plain text ASCII format.
There are three routines in this header file, these are:
void ResultWritingInit(const string &filename,
const int &Nfreqs,
const vector <double> &Freqs,
const int &NumDRNodes,
const vector <int> &DRNodeNums);
void ResultWritingFreq(const string &filename,
const int &ifreq,
const int &NumDRNodes,
const vector <int> &DRNodeNums);
void ResultWritingFinish(const string &filename);
The first function opens a results file for outputting the results into, including checking whether the file is valid and successfully opened. At this point the “offline” information is written, including the date and time of the run, the number and value of the frequencies, whether the solver is an indirect or direct boundary element type. The file is then closed.
The second function opens the file again for each frequency (when the solver has generated a set of results at a specific frequency, this function is called). Then for each data recovery node, the sound pressure level in dB, the acoustic pressure (magnitude and phase angle), the _{} magnitude and phase of the velocities are written, the magnitude of the velocity, the _{} values of the intensity are output.
Finally, the third function opens the results file, places a termination character and closes the file again.
The default spacing for each entry is 12 characters using a 5 digit precision, should alterations be required, this is a simple modification to enact.
The default namespace for this function is ComplexMath_std.
Provides complex number handling for a range of functions. 
This provides handling for complex numbers, including returning real and imaginary components, performing addition, subtraction, multiplication, division and power calculations.
The primary namespace for this routine is BemRoutine.
This is the file which solves the boundary element matrix. It assembles all terms and requests the solution to the matrix problem.

This is the main routine which solves the boundary element acoustic problem. Given a surface region Gamma, with boundary conditions defined as velocity terms for each element [u1 u2 u3...]^T, we require the amplitude of source terms sigma on each element which cause the boundary conditions to be satisfied.
So ui = Gij sigma_j.
So the boundary condition at node i has to be satisfied with a distribution of sources at elements j. The weightings are given by fundamental solutions.
This header file contains a number of subroutines, which will be described in turn. The most significant is the indirect solver routine, as shown below.
void IndirectSolution(const double &Cspeed,
const double &rhodensity,
const int &Nfreqs,
const int & ifreq,
const double &Freq,
const double &Pref,
const double &Iref,
const int &NumNodes,
const vector <int> &NodeNum,
const vector <MeshTypes::Node <double> > &NodeCoord,
const int &NumDRNodes,
const vector <int> &DRNodeNums,
const vector <MeshTypes::Node <double> > &DRNodeCoord,
const int &NumofElems,
const vector <int> &ElemNum,
const vector <MeshTypes::Element <int> > &ElemStruct,
const vector <MeshTypes::BCdetails <int> > &BCnum,
const vector < ComplexMath_std::complex <double> > &BCval,
const ComplexMath_std::complex <double> &ST,
const vector <MeshTypes::Node <double> > &SL,
vector <double> &OutputdB,
vector <vector <double> > &OutputVar);
The matrices are initialised, for F = A Sigma, where F are the array of boundary conditions, Sigma is the vectors of unknown potentials and A the matrix of boundary integrals. For every element, we perform a boundary integral for every other element. As the elements are all at different angles to each other, we perform a Jacobian transformation based on nondimensionalised coordinates on the surface of one element to obtain the integral.
For the nontrivial solution, we must also add the boundary integral contribution of a discrete point source, which is calculated using the IntBHelm function.
We then solve the system of equations and based on this solution, generate a series of physical parameters at every single data recovery point.
In order to calculate the integral term for each element, we call the routine IntSourceHelm. This can be described as finding the contribution to the node i from a source placed on element j. Therefore we are observing the contribution of element j. This is determined by examining the amplitudes of the fundamental solutions at each of the element j's node locations.
ComplexMath_std::complex <double> IntSourceHelm(
const MeshTypes::Node <double> &iNode,
const MeshTypes::Element <int> &ObsElemStruct,
const MeshTypes::Node <double> &ObsNode1,
const MeshTypes::Node <double> &ObsNode2,
const MeshTypes::Node <double> &ObsNode3,
const MeshTypes::Node <double> &ObsNode4,
const double &Freq, const double &C,const double &rho);
The code first calculates the location of the node i. It then immediately branches into a loop depending on whether the element j is triangular or quadrilateral. This then follows the same concept described in section VI.
The value of u* for each node making up element j is obtained by examination of the distance between this element’s node and the original i node.
For each node i, the contribution to the boundary condition value must be included. This routine calculates the required contribution of the pressure at node i from the source.
ComplexMath_std::complex <double> IntBHelm(
const MeshTypes::Node <double> &iNode,
const ComplexMath_std::complex <double> &Bsourcestr,
const MeshTypes::Node <double> &Bsourceloc,
const double &Freq, const double &speed,const double &rho);
The location of the source is provided by the variables (Bsourceloc.x, Bsourceloc.y, Bsourceloc.z). The distance between the source and the node at location i is then,
r=( (Bsourceloc.xiNode.x)^2
+ (Bsourceloc.yiNode.y)^2
+ (Bsourceloc.ziNode.z)^2)^0.5;
r = abs(r);
_{}
The fundamental solution is provided through the HelmFunSol function, as follows.
ComplexMath_std::complex <double> HelmFunSol( const double &R,
const double &speed,
const double &f);
void IntV(
const MeshTypes::Node <double> &iNode,
const ComplexMath_std::complex <double> &Bsourcestr,
const MeshTypes::Node <double> &Bsourceloc,
const double &Freq, const double &speed,const double &rho,
ComplexMath_std::complex <double > &Velx,
ComplexMath_std::complex <double > &Vely,
ComplexMath_std::complex <double > &Velz,
ComplexMath_std::complex < double > &V);
This provides
_{}
This file contains the make instructions for the overall program executable. Currently includes options for operating system specific choices, i.e. Windows or Linux directives. 
This routine is designed based on the Numerical Recipes subroutines to solve the system of equations A x = ( L U ) x = L ( U x ) = B.
L y = B
U x = y
Ax = b
We assume we know the square matrix A[0..n1][0..n1] and also the vector b[0..n1]. We return the inverse matrix A^1 and the solution vectors x.
Require assessment. 
This routine performs an LU decomposition on a matrix. This is a slightly more numerically stable method of doing an analysis rather than inverting the main matrix. Given a typically square matrix a[1..n][1..n] this routine replaces it by the LU decomposition of a rowwise permutation of itself. The matrix A and and size n are input, a is output. indx[1..n] is an output vector that records the row permutation effected by the partial pivoting; d is output as + 1 depending on whether the number of row interchanges was even or odd, respectively. This routine is used in combination with lubksb to solve linear equations or invert a matrix.
void DJOLudcmp( vector <vector <double> > &A,
const int n, vector <int> &indx,double d)
void DJOLubksb( vector < vector <double> > &A,
const int n, vector <int> &indx, vector <double> &B)
The versions of these for complex variables and matrices are provides also with the same input format.
void DJOLudcmpc( vector <vector <ComplexMath_std::complex < double> > > &A,
const int n, vector <int> &indx, double d)
void DJOLubksbc( vector < vector <ComplexMath_std::complex < double> > > &A,
const int n, vector <int> &indx,
vector <ComplexMath_std::complex < double > > &B)
This is a simple Gauss Jordan solver based on Numerical Recipes method for solving systems of equations of the type Ax = b.
We assume we know the square matrix A[0..n1][0..n1] and the vector b[0..n1]. We return the inverse matrix A^1 and the solution vectors x.
Require assessment. 
void Gaussjc( vector <vector <ComplexMath_std::complex < double > > > &A,
const int n,
vector <ComplexMath_std::complex < double> > &B)
Three sets of results are provided for illustration of the capability of the Sonitus BE software. These are all based on a solid tyre geometrical shape, resting on an infinite plane. The data recovery mesh is also based on the automotive case. The first case is where a source is placed in the farfield and the surface pressure is calculated when there is no solid geometry anywhere in the domain. This relatively trivial solution, replicates the standard Green’s function solution for a monopole source propagating outwards, with a mirror image source correcting the result for the ground plane.
The second result is the case where the solid tyre geometry is input. At this point a sound amplification is possible due to the diverging surfaces near to the contact patch where amplifications of up to 20dB are possible for certain frequencies.
The final case is the sound pressure inside a tyre surface with a sound source and a damping layer (modelled as a surface with a prescribed complex impedance). This not only demonstrates the use of the solver for interior against exterior domains, but also demonstrates the ability to mix boundary conditions, with different surfaces containing different damping materials.
Figure 13 Examples of solid geometry defined by an unstructured surface mesh composed of triangular elements and one with a mixture of quadrilateral structured elements (on the wheel arch) and a mixture of triangular and quadrilateral elements on the wheel surface.
The solver has the option (through the input data file) to allow the incorporation of an infinite plane to provide a ground plane for any solid geometry. A finite area section of this is shown in Figure 14 where the normal impedance relation between pressure and velocity can be prescribed as a complex value.
The use of a complex impedance allows the efficient inclusion of a ground surface which has variable levels of acoustic damping.
Figure 14 Ground plane (similar to an image source)
The data recovery array is a series of nodes in the xyz domain, see Figure 15.
Figure 15 Data recovery mesh for automobile tyre noise
The freefield solution is a relatively trivial case, where the solid geometry in the exterior domain is extremely small, (a cube surface comprised of fifty nodes is used as a basis for the solution), with quadrilateral regular mesh spacing for integration of surface integrals. The tyre data recovery mesh shown in Figure 15 is used to generate the sound pressure for a range of excitation frequencies between 0.1 and 2kHz for a source located at (X,Y,Z)=(1.0,0.0,0.1)m..
The strength of the source is unity, providing a relatively large sound pressure at the data recovery positions, which are fairly close to the source position. The sound pressure at an excitation frequency of 100, 400 and 1300 Hz are provided in Figure 16.
Figure 16 Sound pressure for a source in the freefield (note the scale of the amplitude is different at each frequency).
The solution for a sound source in a free field with no solid geometry leads to a spherically propagating pressure distribution. As soon as a solid geometry is introduced into the environment, it becomes possible for the sound field to scatter, leading to regions where for some frequencies, an amplification effect is possible.
In this example, the amplification due to the geometry of a tyre is investigated, as the tyre surface diverges from the road in a pattern similar to a gramophone horn. When the tyre vibrates close to the road surface, this leads to sound close to the contact patch at the throat of the horn. Any sound source is therefore amplified, see Figure 17.
Figure 17 Diverging surfaces leading to a horn geometry where sound sources may be amplified
In this case, rather than placing the sound source in the horn area, it is placed in the farfield and the tyre surface becomes the observer point. This reciprocal technique is used to simplify the number of calculations required.
The sound pressure at an excitation frequency of 100, 400 and 1300Hz are shown in
Figure 18 Sound pressure for a source with a solid geometry in the domain (note the scale of the amplitude is different at each frequency).
The solver code can calculate the sound pressure in an enclosed domain or an infinite domain (where we assume Sommerfeld’s radiation condition holds). However, it is essential that in the creation of the solid geometry mesh and data recovery meshes, that the direction of element normals are consistent (either all pointing into the domain or out of the domain). The code provides a check to ensure that the elements are all consistent with each other and will generate an error message should this not be the case.
There are two main types of geometrical element which generate a surface area. These are triangular and quadrilateral. The solver can integrate both of these and the choice is largely dictated by the mesh generation software, the mesh density and shape of the solid geometry.
A number of monopoles can be prescribed at different locations, with different strengths, to generate an excitation sound field. The amplitude of each source can be complex, allowing a different phase relation to be produced between different sources.
In order to simulate a ground surface, an infinite impedance ground plane can be introduced into the calculation, similar to the effect of adding image sources into the model.
FreeField.dat
21 Frequencies
50 Nodes
4160 Data recovery nodes
Quad elements for structure
Quad data recovery elements
Density 1.12kgm3
Speed of sound 343m/s
Zero velocity boundary conditions
SOURCE 1 1.00000E+000 0.00000E+000 1.00000E+000 0.00000E+000 1.00000E001
Debug problems:
i) Windows version of the LU decomposition requires checking with new maths header files.
ii) Sample files need validating for inclusion in final release. Freefield.dat needs two more elements to ensure an enclosed domain.
iii) Preand Post processing scripts.
iv) The velocity term in the boundary condition specifies the magnitude of the velocity on the element. It could be an improvement to only specify the normal velocity on the element.
v) The contribution of the mesh to the velocity in the data recovery section is not coded at this time. Therefore, the only contribution is from the monopole source.
This manual has provided the details of a numerical solver for a boundary element solver code designed specifically for acoustics applications. Although the method is similar to many other engineering applications, for example, heat transfer, the variables are all related to acoustic fields.
The Sonitus BE acoustic solver provides a comprehensive tool for the investigation of acoustic problems within a wide range of industrial sectors. This manual provides a selection of the numerical techniques used to produce this program. The manual sets out the physical processes which occur in the code, and provides guidance on the physical parameters being calculated.
Sample files have been provided so that a user may write their own input geometry and mesh files for processing. For further information, see the contact details at the beginning of this document.
C. A. Brebbia, “The Boundary Element Method for Engineers”. Pentech Press, 1980.
C. A. Brebbia and S. Walker, “Boundary Element Techniques in Engineering”, NewnesButterworths, 1980.
Abramowitz and Stegun eds. (1972), Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, New York: Dover Publications, ISBN 9780486612720.
F. Fahy and P. Gardonio, “Sound and structural vibration”. Academic Press, 1985.
K. F. Graff, “Wave Motion in Elastic Solids”. Clarendon Press, Oxford, 1975.
D. G. Crighton, A. P. Dowling, J. E. Ffowcs Williams, M. Heckl, and F. G. Leppington,
“Modern methods in Analytical Acoustics”. SpringerVerlag, 1992.
M. Heckl, L. Cremer and E. E. Ungar, “Structure Borne Sound”, 2nd Edition.
SpringerVerlag, Berlin, 1988.
D. E. Bourne and P. C. Kendall, “Vector Analysis and Cartesian Tensors”, 3rd Edition. Chapman and Hall, 1992
A. P. Dowling and J. E. Ffowcs Williams, “Sound and Sources of Sound”. Ellis Horwood, 1989.
W. H. Press, W. T. Vetterling and B. P. Flannery, “Numerical recipes in Fortran: the art of scientific computing”. Cambridge University Press, 1993.
M. R. Spiegel, “Mathematical Handbook of Formulas and Tables”. Schaum, 1998.
E. Kreyszig, “Advanced Engineering Mathematics”, 8^{th} Edition. Wiley, 1998.
Flowers, “An Introduction to Numerical Methods in C++”, Clarendon Press, Oxford, 1995.
A
Abramowitz and Stegun, 37, 54
acoustic, 11, 13, 16, 17, 18, 20, 22, 41, 42, 43, 47
Acoustic Pressure, 36
amplitude, 17, 18, 19, 20, 21, 43, 50, 51, 52
ASCII, 42
B
BEM, 11
boundary, 11, 13, 16, 17, 18, 19, 20, 21, 22, 26, 28, 36, 39, 40, 41, 42, 43, 44, 47, 52
C
CAD, 11, 13
complex, 13, 17, 39, 43, 44, 45, 46, 47, 52
D
Data recovery mesh, 18, 48
discretisation, 13, 22
domain, 11, 13, 16, 17, 18, 19, 20, 21, 22, 24, 35, 46, 48, 49, 51
E
Element, 1, 11, 16, 20, 40, 43, 44, 52, 54
error, 13, 21, 37, 40, 41, 51
Excitation frequency, 17
F
finite element, 11, 13, 19
Fundamental solution, 17
G
Gauss Jordan solver, 46
Green's function, 14, 17, 21, 46
H
Helmholtz, 17, 18, 21, 22
I
impedance, 13, 41, 42, 47, 52
Input, 33, 35, 41, 52
Interpolation, 37, 38
J
Jacobian transformation, 44
L
Linux, 11, 12
LU decomposition, 46
M
MATERIAL, 41
Matlab, 11, 33, 35
monopole, 17, 18, 19, 46
N
Node, 16, 40, 43, 44, 45
Numerical Recipes, 37, 45, 46
O
Output, 35, 36
P
phase, 18, 42, 52
potential, 14
precision, 37, 39, 43
pressure, 11, 17, 18, 19, 20, 33, 35, 42, 44, 46, 47, 49, 50, 51
Q
quadrilateral, 16, 35, 37, 38, 41, 47, 49, 52
R
reference intensity, 41
reference pressure, 41
reference work, 41
Requirements, 12, 14, 15, 33, 35, 36
S
Sample files, 46
SOURCE, 42, 52
Subroutine, 37
surface, 11, 13, 14, 16, 17, 18, 19, 20, 21, 22, 24, 34, 37, 38, 41, 43, 44, 46, 47, 49, 51, 52
Surface, 16
T
triangular, 16, 35, 37, 38, 40, 41, 47, 52
V
velocity, 11, 17, 19, 35, 41, 42, 43, 47, 52
Velocity, 36
W
weighting, 14, 17, 21, 37, 38, 39
Windows, 11