Home

Sonitus BE R1 User Manual

1. 16 FIGURE 2 A TRIANGULAR ELEMENT CHARACTERISED BY THREE INDEPENDENT NODES AND A QUADRILATERAL ELEMENT CHARACTERIZED BY FOUR NOTE THAT THE DIRECTION NORMAL OF ALL ELEMENTS MUST BE CONSISTENT oaa 16 FIGURE 3 OMNI DIRECTIONAL SOUND PRESSURE PROPAGATING FROM FUNDAMENTAL SOLUTION m 19 FIGURE 4 DIRECTIONAL SOUND FIELD PROPAGATING FROM A DISTRIBUTION OF MONOPOLES AT LOW AND HIGH FREQUENCIBS dra este a vet EUNT ae ee aee e CENE Toad Ree Ere ushu 20 FIGURE 5 DEFINITION OF A DOMAIN FOR ACOUSTIC CALCULATIONS eene enne 21 FIGURE 6 DISCRETISATION OF THE SURFACE BOUNDING THE SOLUTION DOMAIN THE SURFACE IS CREATED BY JOINING NODES TOGETHER WITH LINEAR SEGMENTS WHEN COMPLETED OVER A THREE DIMENSIONAL DOMAIN THE LINEAR SEGMENTS CREATE SURFACE AREAS 5 23 FIGURE 7 BOUNDARY INTEGRAL SURFACE ILLUSTRATING HOW EACH ELEMENT TAKES INTO ACCOUNT THE EFFECT OF A SOURCE OF UNKNOWN AMPLITUDE LOCATED AT EACH ELEMENT THE PROBLEM IS THEN POSED AS WHAT AMPLITUDE OF ALL OF THE SOURCES YIELDS THE CORRECT KNOWN BOUNDARY CONDITION AT THIS ELEMENT IN MATRIX 000 00 n 26 FIGURE 8 QUADRILATERAL BOUNDARY INTEGRAL SURFACE SHOWING NON DIMENSIONAL COORDINATE FORM OF THE SUREACE DIR
2. SonitusBE P Now we introduce the discretisation of the source term These can be used to include either monopole sources or body forces These are written as B bu dO Our final discretised indirect boundary integral equation is left as N _ j l N q gt B j l In this code the matrices are hard coded for N elements as F Ao u G Gi Uy _ G Gi C Guy KNOWN KNOWN UNKNOWN Element J Element I Figure 7 Boundary integral surface illustrating how each element takes into account the effect of a source of unknown amplitude located at each element The problem is then posed as what amplitude of all of the sources yields the correct known boundary condition at this element in matrix form The effect of each fundamental solution located at the nodes are shown in Figure 7 where it may be seen that the vector r is never zero as two nodes are never co located Sonitus BE Page 26 of 58 c SonitusBE iN INTEGRATING A BOUNDARY ELEMENT The integration of the boundary element is usually the most complex part of the boundary element code due to the fact that the surface is generally not aligned with the principle axes of the problem see for example the example in any of the previous figures but rather with the normals pointing in an arbitrary direction yet consistent with each other In thi
3. _ _ SonitusBE Acoustic Boundary Element Solver Dr Dan J O Boy http www djoboy co uk R1 Revision 1 Copyright 2008 Dr Dan J 4 The Green Exton Oakham Rutland LE158AP United Kingdom enquires djoboy co uk 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 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 TABLE OF CONTENTS TABLE OF CONTENTS LIST OF FIGURES LIST OF TABLES 1 INTRODUCTION REQUIREMENTS BACKGROUND OF THE BOUNDARY ELEMENT METHOD IV NOMENCLATURE V KEY TERMINOLOGY VI MATHEMATICAL FOUNDATION OF THE BOUNDARY ELEMENT METHOD PRINCIPLE OF THE INDIRECT BOUNDARY ELEMENT METHOD DISCRETISATION OF THE SURFACE DOMAIN INTEGRATING A BOUNDARY ELEMENT
4. 2010210 Outputvar i 1 Sound pressure magnitude N m 2 y ia 2 Outputvar i 2 Sound pressure phase in degrees rsd acun Sonitus BE Page 32 of 58 The velocity is given by V Qliolc 1 where is distance from source element to the observer data recovery node The velocity is different from the pressure in that the former is a scalar whereas the velocity can be defined in terms of components in the three Cartesian axes SonitusBE V V MIS x RC Outputvar 1 3 X component of the velocity magnitude m s JRelV Im V Outputvar 1i 5 Y component of the velocity magnitude m s JRe v Outputvar 1 7 Z component of the velocity magnitude m s JRe V Im v Outputvar i 9 Magnitude of the velocity magnitude m s V Outputvar i 4 X component of velocity phase degrees 180 arctan Re V 1 6 component of velocity phase degrees Im 180 arctan Re V Outputvar 1 8 Z component of velocity phase degrees oe The acoustic intensity can be determined for each of these coordinate axes through the power dissipation analogy in electrical components The intensity may be written in terms of the pressure and velocity using 1 I pV 5 Rena V where the overbar denotes the complex con
5. to and non dimensionalised surface coordinates 4 and 6 which vary from 0 to 1 This provides the derivative at 1 the values of and double DZETAlinterpTRI const doubl const doubl const doubl const doubl const doubl Return N1 N3 0 amp N1 amp N2 amp N3 amp zetal amp zeta2 Sonitus BE Page 38 of 58 SonitusBE iN b e Interpolation routine for a triangular element with nodes N to N and non dimensionalised surface coordinates and which vary from 0 to 1 This provides the derivative at 2 the values of and double DZETA2interpTRI const double amp N1 const double amp N2 const double amp N3 const double amp zetal const double amp zeta2 Return 0 N2 N3 2 e Interpolation routine for a quadrilateral element with nodes to N and non dimensionalised surface coordinates 6 and which vary from 1 to 1 double interpQUAD const double amp N1 const double amp N2 const double amp N3 const double amp 4 const double amp zetal const double amp zeta2 Define a set of weighting functions Wi W3 W3 W equated to Vla 20a 25 y Va eoa e Kapag a oa Interpolated return value v4 N w N e Interpolation routine for a quadrilateral element with nodes N to N and non dimensionalised surface coordinates and 4 which vary fr
6. QUADRILATERAL ELEMENTS TRIANGULAR ELEMENTS SOURCE DEFINITION AND LOCATION DATA RECOVERY PARAMETERS THE SOLUTION METHOD FOR A MATRIX PROBLEM VII INPUT FILE DEFINITION III OUTPUT FILE DEFINITION IX PROGRAM OVERVIEW SONITUSBE CPP X SUBROUTINE FUNCTIONS ERRORHANDLE H COMPLETE INTERP H COMPLETE MESSAGES H COMPLETE SOLUTIONCHOICE H COMPLETE VERIFICATION H COMPLETE STRUCTURES H COMPLETE FILEREADERSONITUS H COMPLETE SonitusBE s i 14 15 20 23 25 27 27 29 30 30 34 34 36 37 37 38 38 38 40 40 41 41 42 Sonitus Page 5 of 58 SonitusBE s RESULTOUTPUT H COMPLETE 43 COMPLEX H COMPLETE 44 INDIRECTSOLVER H COMPLETE 44 MAKEFILE COMPLETE 46 DJO_LU H INGONIPIER 46 GAUSSJC H COMPLETE 47 XI SAMPLE FILES 47 SOLID GEOMETRY 48 INFINITE PLANE SURFACE 48 DATA RECOVERY MESH EXTERIOR 49 DATA RECOVERY MESH INTERIOR 49 FREE FIELD SOLUTION 50 SOUND PRESSURE IN A DOMAIN WITH A SOLID GEOMETRY 52 SOUND PRESSURE INSIDE AN ENCLOSURE WITH DAMPING MATERIAL INCLUDED ON ONE SURFACE 52 NOTES 52 EXTERIOR AND INTERIOR DOMAINS 52 ELEMENT TYPES 53 SOUND SOURCES 53 INFINITE GROUND PLANE 53 INPUT FILE TYPES DAT 53 XIV DEBUG REQUIREMENTS 53 XV SUMMARY 53 XVI REFERENCES 55 Sonitus BE Page 6 of 58 SonitusBE iN a LIST OF FIGURES FIGURE 1 NODE CHARACTERISED BY THREE INDEPENDENT COORDINATES
7. 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 ul 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 amp Cspeed const double amp rhodensity const int amp Nfreqs const int amp ifreq const double amp Freq const double amp Pref const double amp Iref const int amp NumNodes const vector int amp NodeNum const vector lt MeshTypes Node double gt amp NodeCoord const int amp NumDRNodes const vector int amp DRNodeNums const vector lt MeshTypes Node double gt amp DRNodeCoord const int amp NumofElems const vector int amp ElemNum const vector lt MeshTypes Element int gt amp ElemStruct Sonitus BE P
8. distance between the source and the node at location 1 is then Bsourceloc x iNode x 2 Bsourceloc y iNode y 2 Bsourceloc z iNode z 2 0 5 abs r Bsourcestr HelmFunSol r speed Freq The fundamental solution is provided through the HelmFunSol function as follows ComplexMath std complex double HelmFunSol const double amp R const double amp speed const double amp f void IntV const MeshTypes Node double amp iNode const ComplexMath std complex double amp Bsourcestr const MeshTypes Node double amp Bsourceloc const double const double amp speed const double amp rho ComplexMath std complex double gt amp Velx ComplexMath std complex double gt amp Vely ComplexMath std complex double amp Velz ComplexMath std complex lt double gt amp V This provides Makefile COMPLETE 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 DJO LU h This routine is designed based on the Numerical Recipes subroutines to solve the system of equations Ax LU x L Ux B Ly B Ux y Ax b We assume we know the square matrix A 0 n 1 0 n 1 and also the vector b 0 n 1 We return the inverse matrix 1 and the solution vectors x Requ
9. each of the element j s node locations ComplexMath std complex double IntSourceHelm const MeshTypes Node double amp iNode const MeshTypes Element int amp ObsElemStruct const MeshTypes Node double amp ObsNodel const MeshTypes Node double amp ObsNode2 const MeshTypes Node double amp ObsNode3 const MeshTypes Node double amp ObsNode4 const double amp Freq const double amp C const double amp 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 amp iNode const ComplexMath std complex double amp Bsourcestr const MeshTypes Node double amp Bsourceloc Sonitus BE Page 45 of 58 SonitusBE iN bam const double amp Freg const double amp speed const double amp rho The location of the source is provided by the variables Bsourceloc x Bsourceloc y Bsourceloc z
10. 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 y z magnitude and phase of the velocities are written the magnitude of the velocity the x 5 values of the intensity output Finally the third function opens the results file places a termination character and closes the file again Sonitus Page 43 of 58 SonitusBE The default spacing for each entry is 12 characters using 5 digit precision should alterations be required this is a simple modification to enact complex h COMPLETE 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 IndirectSolver h COMPLETE
11. hypothetical predictions The main prediction tools for the prediction of the acoustic field include theoretical models finite elements methods FE ray methods geometrical acoustics approximations statistical energy analysis variational energy methods or the boundary element method The main features leading to advantages when the methods are compared are provided in this section Finite element FE methods are the most widely used due to the efficient and simple integration of computer aided design tools in the engineering workplace As computer aided engineering design usually involves the creation of CAD drawings of the object it is easy to take those drawings and incorporate them into the acoustic analysis Finite element tools require the discretisation of a full volume and subsequent meshing of that volume Although it produces extremely good solutions it can tend to be complex requiring high computing power and for higher frequency prediction accurate geometry The computing requirement increases with the upper analysis frequency and size of the geometry The statistical energy method was designed to produce approximate acoustic prediction for segmented rooms using diffuse field assumption Due to the averaging process the prediction can be error prone however the method is fast and efficient As diffuse theory is assumed leading to average energy flow between compartments the method is not suitable for low frequencies yet this then
12. means that the design process is tolerant to variations in geometry Ray tracing ray theory or geometrical acoustics approximations are extremely useful for enclosed domains bounded exterior spaces or rooms where diffuse theory could be applied Each acoustic wave is represented as a ray tracing a predictable path through space with the incident and reflection angles from a surface taken into account Damping can be applied at the reflective surface through a variable impedance boundary condition As ray theory has a lower frequency limit it is usually used for high frequency studies of room acoustics and reverberation times The computational time varies depending on the level of complexity of the room geometry as small changes in room angles can lead to large changes in the ray transmission paths hence the importance of overall diffuse field assumption Sonitus BE Page 13 of 58 bem Theoretical models are used in studies to gain an understanding of the physical sound transmission processes that can occur By idealising a geometry using simple geometrical primitives or through using significant assumptions it may be possible to reduce the complexity of the problem down to a level where a simple theoretical model can provide insight into the problem order of magnitude calculations and trend analysis in a parametric study As these models are often much simpler than reality they often are far less computationally expensive or time consumi
13. of challenges with this method most notably that few geometries exist in reality where there is simple analytical solution for the boundary integral It is more likely that the geometry is a component in an automotive aeronautical aerospace or industrial environment For this reason the integral has to be completed using a numerical method which requires that the surface be discretised into a series of nominally evenly spaced nodes see Figure 1 for a description of the node and Figure 6 to see how the surface is broken up into elements separated by these nodes There N elements from i 1 to where we may describe two different elements by their indices i and j The value of on element i is denoted u while the value attached to element j is We may then begin to write the discretised formulation in matrix notation using a summation over N elements At this point we do not consider the implications of j in the mathematical derivation In addition we temporarily leave the source term B out of the equations for brevity N c 1 j l 1 N 4 7 20 jr The boundary integrals are provided with their own specific notation for the matrices such that G var and H 1 1 q gt c HR j l Finally for the boundary integral for the fluxes we may write that if i j then H j gt Whereas for i j Hii Hii 1 2 Sonitus BE Page 25 of 58
14. 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 Horn geometry 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 far field 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 Sound pressure inside an enclosure with damping material included on one surface XII NOTES Exterior and interior domains 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 soli
15. surface with the first section denoted and the second as shown in Figure 5 On each surface are located surface potentials of unknown amplitude which give rise to certain physical conditions on the surface On the first I for the acoustic calculation is a term while in on the surface we define the property of interest to be the normal derivative the flux q Ou On It is possible to write the physical boundary conditions of pressure or velocity in terms of the amplitudes of the fundamental solutions By converting physical conditions such as pressure or velocity into just the amplitudes of the fundamental solutions the problem can be decomposed into a generic solution to a matrix The exact solutions for the physical properties are denoted by the functions and q on the sections of the surface respectively This allows an error on the surfaces to be defined as g 6 4 4 I Sonitus Page 21 of 58 c SonitusBE iN The aim of the boundary element method is to determine the amplitude of a distribution of potentials so as to minimise the error This error is distributed over the surface according to a weighting function w which may be a fundamental solution to the governing equation in the domain The distribution of the error in the domain is written For the Helmholtz problem the governing equation is V u 0 wh
16. what value it takes Finally in order to obtain a non trivial 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 Resultoutput h COMPLETE 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 post processed 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 amp filename const int amp const vector double amp Fregs const int amp NumDRNodes const vector int amp DRNodeNums void ResultWritingFreq const string amp filename const int amp ifreq const int amp NumDRNodes const vector int amp DRNodeNums void ResultWritingFinish const string amp filename The first
17. ECTION 28 FIGURE 9 QUADRILATERAL LINEAR INTEGRATION METHOD 28 FIGURE 10 TRIANGULAR LINEAR INTEGRATION METHOD 29 FIGURE 11 BOUNDARY INTEGRALS FROM THE SURFACE THE DATA RECOVERY MESH 32 FIGURE 12 GENERATION OF A SURFACE MESH USING COMPUTER AIDED DESIGN TOOLIS 35 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 48 FIGURE 14 GROUND PLANE SIMILAR AN IMAGE 5 49 FIGURE 15 DATA RECOVERY MESH FOR AUTOMOBILE NOISE 0 0 eene eene nte enne titan nis 49 FIGURE 16 SOUND PRESSURE FOR A SOURCE IN THE FREE FIELD NOTE THE SCALE OF THE AMPLITUDE IS DIFFERENT AT EACH seen ene eti GIS 51 FIGURE 17 DIVERGING SURFACES LEADING TO A HORN GEOMETRY WHERE SOUND SOURCES MAY BE 52 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 52 Sonitus BE Page 7 of 58
18. EMATICAL FOUNDATION OF THE BOUNDARY ELEMENT METHOD A typical boundary element code has several advantages over traditional finite differencing calculation methods As with finite element codes a geometry is described and a mesh generated however only the bounded surfaces require meshing which allows not just a speed advantage but also the possibility to investigate domains which are relatively large in comparison to the geometry being examined Broadly speaking boundary conditions are applied to the surface under investigation and solutions found for the amplitude of potentials applied to the surface which yield those boundary conditions Once these amplitudes are known they may be used to obtain the displacement velocity or other parameter inside the domain The aim of this section is to provide an illustration of the boundary element method through the work of Brebbia and Walker C A Brebbia The Boundary Element Method for Engineers Pentech Press 1980 C A Brebbia and S Walker Boundary Element Techniques in Engineering Newnes Butterworths 1980 Sonitus BE Page 20 of 58 SonitusBE P e Surface on Problem domain exterior Ou on q 5 2 Figure 5 Definition of a domain for acoustic calculations A domain is defined by the symbol and a surface which is closed inside it This surface is separated into two sections joined together to create a closed
19. ETE 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 eac
20. For further information see the contact details at the beginning of this document Sonitus BE Page 54 of 58 SonitusBE XVI REFERENCES C Brebbia The Boundary Element Method for Engineers Pentech Press 1980 C A Brebbia and S Walker Boundary Element Techniques in Engineering Newnes Butterworths 1980 Abramowitz and Stegun eds 1972 Handbook of Mathematical Functions with Formulas Graphs and Mathematical Tables New York Dover Publications ISBN 978 0 486 61272 0 F Fahy and P Gardonio Sound and structural vibration Academic Press 1985 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 Springer Verlag 1992 M Heckl L Cremer and E E Ungar Structure Borne Sound 2nd Edition Springer Verlag 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 Edition Wi
21. LUTION 2 FREQUENCY RAD SEC Q DOMAIN Table 3 Notation for Greek notation Sonitus BE Page 14 of 58 SonitusBE ae A MATRIX OF BOUNDARY INTEGRALS FOR ____ 7 BOUNDARYELEMENTMETHOD COMPUTER AIDED ENGINEERING COMPILER LANGUAGE F VECTOR OF PRESCRIBED BOUNDARY CONDITIONS FOR EACH ELEMENT FE FINTEELEMENT O IREF REFERENCEINTENSTY O o PREF REFERENCEPRESSURE 6 SPEEDOFSOUND 00000 E ELEMENTNUMBR 00000004 1 H 2 HANKEL FUNCTION IMAGINARY NUMBER JACOBIAN WAVENUMBER NODE NUMBER Ou On FLUX OF PRESSURE DISTANCE FROM SOURCE ELEMENT TO OBSERVER ELEMENT VARIABLE PRESSURE FUNDAMENTAL SOLUTION PRESSURE FUNDAMENTAL SOLUTION NORMAL PRESSURE GRADIENT WEIGHTING FACTOR x COORDINATES IN CARTESIAN SPACE Table 4 Notation for English notation V KEY TERMINOLOGY Acoustic The medium in which a sound wave may propagate is an acoustic medium A problem relating to sound pressure perturbations in the audible frequency range is considered an acoustic problem Typically this medium cannot sustain shear waves only longitudinal waves compression Sonitus BE Page 15 of 58 Speed of sound It is assumed that in the elastic acoustic medium distance travelled sound wave in a given time does not vary this is the speed of sound Fo
22. SonitusBE LIST OF TABLES TABLE 1 MINIMUM RECOMMENDED COMPUTER SPECIFICATIONS 00000ssesseseseeeseseseseseseseseseseseseseseseseseeeeens 12 TABLE 2 COMPILER RECOMMENDATIONS FOR SOURCE CODE 0 0 20 13 TABLE 3 NOTATION FOR GREEK enn n nennen enne penes enne ene n nenne nnne 14 TABLE 4 NOTATION FOR ENGLISH NOTATION 15 TABLE 5 FUNDAMENTAL SOLUTIONS FOR ACOUSTIC WAVE PROPAGATION IN FREE SPACE 18 TABLE 6 INPUT FILE iie custards ente oe eine uy dya 36 TABLE 7 OUTPUT FILE DESCRIPTION 37 Sonitus BE Page 9 of 58 SonitusBE P Sonitus BE _ Vl L INTRODUCTION Sonitus BE is a versatile boundary element solver for acoustic engineering problems It can provide the sound pressure fluid velocity or intensity in an enclosure or in the free field for a range of problems Boundary Element Methods BEM are a class of numerical solution method which are efficient for the solution of acoustic problems and differ when compared to finite element methods as they do not require a volume mesh to be generated only a surface mesh reducing the complexity of the problem domain This document details some of the technical code comprising the Sonitus BE program including subroutin
23. age 44 of 58 SonitusBE ao const vector lt MeshTypes BCdetails lt int gt gt amp BCnum const vector lt ComplexMath std complex lt double gt gt amp BCval const ComplexMath std complex lt double gt amp ST const vector lt MeshTypes Node double gt amp SL vector double amp OutputdB vector vector double amp 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 non dimensionalised coordinates on the surface of one element to obtain the integral For the non trivial 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
24. al element types void checkBCnumbers const int amp NumofElems const int amp 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 Structures h COMPLETE 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 x z 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 Sonitus BE Page 41 of 58 SonitusBE iN bm sg 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 FileReaderSonitus h COMPL
25. ample 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 triangular 16 35 37 38 40 41 47 52 velocity 11 17 19 35 41 42 43 47 52 Velocity 36 Sonitus Page 56 of 58 SonitusBE weighting 14 17 21 37 38 39 Windows 11 Sonitus BE Page 57 of 58 SonitusBE es Sonitus BE Page 58 of 58
26. ate 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 Sonitus Pre BE display standards Input file end line Table 6 Input file description VIII OUTPUT FILE DEFINITION The results file is generated automatically using the software Sonitus Post BE 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 Sonitus BE Page 36 of 58 0 SonitusBE iN dm 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 X Component of Velocity Magnitude X Component of Velocity Phase Angle z Component of Intensity Component of Intensity Output file en
27. cal Recipes method for solving systems of equations of the type Ax b We assume we know the square matrix 0 1 0 1 and the vector b 0 n 1 We return the inverse matrix 1 and the solution vectors x Require assessment void Gaussjc vector vector ComplexMath std complex lt double amp A const int n vector ComplexMath std complex double amp B XI SAMPLE FILES 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 far field 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 Sonitus BE Page 47 of 58 SonitusBE t e sr 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 prescribe
28. can be directly linked to the parameters of interest for acousticians such as fluid velocity or pressure To take an example a sphere of radius a pulsating with a frequency c 2zF and a radial velocity v a at the surface generates a spherical wave The oscillatory volume flow at the surface has an amplitude Q that is given by the monopole strength 4za v a In the limiting case of a monopole point source the wavelength dimension ka lt lt 1 we obtain for the acoustic pressure p r V plico exp ikR In the code listed below the source strength is user provided as a complex constant B Therefore both the velocity V and the pressure p can be written in terms of this fundamental solution p ipoQu 9 intensity also be related to acoustic pressure velocity through analogy with power dissipation of electrical circuits 1 I where the overbar denotes the complex conjugate is taken Data recovery mesh The boundary element method starts with a known boundary condition on a geometrical surface and the unknown distribution of fundamental solutions across the domain which yield that boundary condition is formulated in terms of a matrix problem Once this distribution is known the sound pressure can be obtained everywhere in the domain Data recovery points are the locations where the user is interested in knowing these resulting sound pressures intensiti
29. cial packages their minimum specifications should be consulted For post processing in Matlab in Linux it is recommended that a good graphics card be included with OpenGL support ATTRIBUTE MINIMUM SPECIFICATION MICROSOFT WINDOWS WINDOWS SERVICE PACK3 OPENSUSE 10 0 MEMORY DEPENDENT ON 2GB RAM MINIMUM GEOMETRY SIZE AND MESH DENSITY FILES CPU2GHZ GRAPHICS FOR MATLAB V7 POSTPROCESSING Table 1 Minimum recommended computer specifications To compile the source code in Linux type make at the console To compile the source code in Windows using MinGW type g SonitusBE cpp o SonitusBE exe Alternatively type at the console To compile using operating specific directives used for pausing the screen to show an error message include the flag dWINDOWS or dLINUX in the compile command string See makefile for example Under Windows XP and using the MinGW gcc compiler it is possible to create a structure which is large enough to cause memory issues The solution is to increase the available stack size by supplying the linker with a sufficiently large maximum memory allocation so that a stack overflow never occurs This is currently available for MinGW by including the flag W1 stack 1073741824 at compile time This option leads to a maximum Windows stack size of 1Gb Note that due to dynamic allocation Linux does not require this maximum stack size to be increased The
30. computer code has been written using dynamically allocated memory rather than static matrix sizes and through passing a memory reference pointer to subroutines rather than copying variables however engineering problems typically require the solution of large matrices While there are methods to decompose large matrices into smaller sizes for more appropriate handling this requires application specific analysis which can often lead to longer development time for diminishing returns Sonitus BE Page 12 of 58 dm SonitusBE For compiler recommendations using modified source code see Table 2 OPERATING SYSTEM COMPILER MICROSOFT WINDOWS MINGW PORTING LINUX GCC C IMPLEMENTATION Table 2 Compiler recommendations for source code modification The command line is SonitusBE INPUTFILE OUTPUTFILE For Microsoft Windows running in a console the typical output would be SonitusBE exe SampleFiles FreeField dat output rslt or SonitusBE out SampleFiles FreeField dat output rslt The output is directed to the screen by default To direct this to an output file SonitusBE exe SampleFiles FreeField dat output rslt gt screenoutput txt BACKGROUND OF THE BOUNDARY ELEMENT METHOD There are several different methods available for the prediction of the acoustic field in either an enclosed interior space or an exterior domain These methods are usually designed to support experimental measurements or
31. d 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 Sonitus BE Page 52 of 58 SonitusBE Element types 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 Sound sources 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 Infinite ground plane 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 XIII INPUT FILE TYPES DAT FreeField dat 2 Frequencies 50 Nodes 4160 Data recovery nodes Quad elements for structure Quad data recovery elements Density 1 12kgm 3 Speed of sound 343m s Zero velocity boundary conditions SOURCE 1 1 00000E 000 0 00000 000 1 00000 000 0 00000E 000 1 00000 001 XIV DEBUG REQUIREMENTS Debu
32. d 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 Solid geometry 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 Infinite plane 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 Sonitus BE Page 48 of 58 SonitusBE s 008 Figure 14 Ground plane similar to an image source Data recovery mesh exterior The data recovery array is a series of nodes in the xyz domain see Figure 15 Data recovery mesh for extenor tyre hom amplification Z distance from ground level m 04 Y distance 24 8 X distance m Figure 15 Data recovery mesh for automobile tyre noise Data recovery mesh interior Sonit
33. d line Table 7 Output file description IX PROGRAM OVERVIEW Sonitus BE cpp 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 pseudo code 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 Sonitus BE Page 37 of 58 X SUBROUTINE FUNCTIONS SonitusBE aw cu 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 ro
34. e calls The program is available in Windows or Linux specific versions Sonitus BE is the numerical solver that takes the input data files describing the acoustic problem and constructs the matrix problem It then automatically solves this matrix problem and uses the data recovery points specified in the input file to generate a results file The program runs as a standalone application without graphics interface This deliberately allows the user to choose the graphical user interface and results processing software from a range of sources e Commercial CAD packages Commercial meshing packages e Matlab standalone code for input geometry Recommended e Sonitus Pre BE input file processor Recommended e Commercial result processing software e Matlab standalone code for results processing Recommended e Sonitus Post BE result processing software Recommended In this manual the requirements of the boundary element program are provided in terms of the minimum recommended specifications in section II A background of the different methods available for the study of acoustic problems is provided in section which shows the advantages and disadvantages of the boundary element method As the boundary element method is described in this manual is mathematical terms a nomenclature is provided in section IV and a terminology page summarising important concepts in section V The mathematical formulation of a boundary integral method is
35. e particular numerical problems when attempting to calculate the steady state acoustic pressure deflection the frequency equal to zero Hertz 2 In order for the calculation to converge to a satisfactory solution it is required that the zero frequency instead be replaced by a relatively small value for example 0 01Hz Boundary condition Each element has an associated boundary condition which is an initial condition required to be able to solve the problem The types of physical condition which can be applied to the element include a prescribed velocity or pressure including complex impedance These are proportional to the fundamental solution In the input file each element has an associated boundary condition In the code these are denoted boundary condition types 1 or 2 or 3 The pressure on an observer element can be written in terms of a fundamental solution on a source element a distance r away This is a type 1 boundary condition pressure As will be described in the following sections the velocity applied to the element type 2 can be written in terms of a fundamental solution as can the pressure from this observer element to a source element a distance r away The normal velocity can be applied to the element to represent a solid surface or prescribed flow condition The amplitude of the fundamental solution provided by the source element yielding the appropriate boundary condition on the observer element is denoted Q p i
36. ed without this the only sound field that exists is silence The amplitude of this source is b with a location inside Sonitus BE Page 23 of 58 SonitusBE gudr bu aQ shown in Figure 6 surface either integrated as continuous surface we will use the continuous form for the present section or discretised by nodes Joined with linear elements In order to enforce the boundary condition on the surface it becomes necessary to also consider a domain exterior to the surface where no sources exist and the properties u q are governed by the equation 4 0 as the surface is pointing the opposite direction the sign is reversed r As part of the derivation we assume that the solution generated for the exterior field is identical to that generated for the interior field on the surface Therefore we are able to write gudr bu ao qu ar qwar On the surface it is possible to write that as exact terms Rearranging the above equation these terms cancel with each other fou u q d ou do 0 4 0 I If we examine the form of this equation we can define 4 4 as the unknown density of the fundamental solutions or amplitudes of the fundamental solutions across the surface of the domain
37. equencies as seen in Figure 4 the monopoles are spaced within an acoustic wavelength of each other such that the overall sound field is omni directional However at higher frequencies as also shown in Figure 4 the spacing between each monopole becomes more significant compared to the acoustic wavelength and the phase change between each monopole results in a more complicated sound field with lobe patterns This is the primary reason why loudspeaker packages comprise many different speaker sizes each providing sound in a specified frequency range In the boundary element method we assume that we know the boundary conditions applied to the surface of a geometry and require knowledge of the amplitude and phase of a distribution of fundamental solutions on each node that satisfies this boundary condition Sound pressure from monopole dB Frequency 500 Hz Rigid wall m ik v 8 d Distance from monopole m 0 1 2 Figure 3 Omni directional sound pressure propagating from a monopole fundamental solution Sonitus BE Page 19 of 58 m SonitusBE es Sound pressure from baffled piston dB Frequencyz1000 Hz 1 istance troi PAR Distance from piston Sound pressure from baffled piston dB Frequency 8000 Hz 5 1 0 Rigid wall m Rigid wall m 4 Figure 4 Directional sound field propagating from a distribution of monopoles at low and high frequencies VI MATH
38. ere the wavelength is denoted 2 Substitution of the error residuals and weightings fundamental solutions provides the starting weighting residual statement V2u a zwar f u ar A where u describes the fundamental solution to the Helmholtz equation V u also termed the Green s function The first term may be integrated by parts Zu aQ f q q u ar a ar It is noted that we have specified that the boundary conditions on the two parts of the surface must exactly meet the exact solution Therefore this simplifies the equation to Zu han ug ar r integration over domain of the fundamental solution yields only delta term located at point of interest 1 i f qu ar uq ar 66 25 1 This final equation relates value of at the location of point 1 with values of and q over the domain and it applies when the point is inside the interior domain The advantage of boundary integrals appears when we only consider points on the surface of the geometry therefore the point needs to be moved to this boundary Due to the singularity which occurs at this point we analytically represent this through a constant c For a point inside the boundary c 1 for a point the boundary it is 0 5 and outside the boundary it is zero T
39. es or acoustic fluid velocities Although it is described as a mesh for graphical interpretive purposes it can equally be defined as a set of single unconnected points In designing the data recovery mesh or locations of the recovery points it may be of interest to obtain the acoustic pressure close to the solid geometry In order to avoid singularities in the matrix solution it is necessary to locate each data recovery point at least 1 4 of an acoustic wavelength away from the solid geometry element or node The sound field from multiple monopoles A monopole is defined as a solution to the Helmholtz equation which yields a sound field which propagates omni directionally For a three dimensional field this results in an omni directional field where the amplitude falls away as a function of the inverse of the radius This is illustrated in Figure 3 for a monopole originating on the left hand side of the screen and shown propagating to the right hand side of the diagram Sonitus BE Page 18 of 58 SonitusBE p 7 When several monopoles situated domain it is possible to obtain very complicated sound fields through linear superposition leading to regions of high sound pressures from constructive interference and regions of low sound pressure from waves destructively interfering This is most easily illustrated using a large number of monopoles situated close to each other in a baffled piston arrangement At lower fr
40. g problems 1 Windows version of the LU decomposition requires checking with new maths header files 11 Sample files need validating for inclusion in final release Freefield dat needs two more elements to ensure an enclosed domain iii Pre and 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 XV SUMMARY Sonitus BE Page 53 of 58 r 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 SonitusBE iN 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
41. h 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 QUADA 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 QUADA 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 Sonitus BE Page 42 of 58 m density has real and imaginary components as does the subsequent speed of sound The following six lines are then ignored SonitusBE as 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
42. his is the governing boundary element equation for use with a fundamental Helmholtz solution Such solutions for the fundamental equation are shown in Table 5 Sonitus BE Page 22 of 58 SonitusBE s bs s gi The boundary element equation is a continuous integral however the numerical implementation is a discrete equation with truncation errors through the continuous domain being represented by a series of linear joins as shown in Figure 6 The requirements for this discretisation the type of elements which can be created and the number of elements for a given acoustic frequency are discussed later in this manual Surface xu I Fundamental solution domain Figure 6 Discretisation of the surface bounding the solution domain The surface is created by joining nodes together with linear segments When completed over a three dimensional domain the linear segments create surface areas elements PRINCIPLE OF THE INDIRECT BOUNDARY ELEMENT METHOD The formulation of the indirect boundary element method starts with the boundary integral equation qu applied to an interior domain Q The indirect method starts with a solution which satisfies the governing equation in the domain but has unknown coefficients at a number of points At this point we also introduce to the formulation a monopole source located inside the domain which will allow the non trivial solution to be obtain
43. ire assessment 2 III T Sonitus BE Page 46 of 58 This routine performs an LU decomposition on a matrix This is 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 row wise permutation of itself The matrix 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 SonitusBE AN void DJOLudcmp vector vector double gt const int n vector int amp indx double d void DJOLubksb vector vector double gt amp const int n vector int amp indx vector double amp 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 lt double amp A const int n vector int amp indx double d void DJOLubksbc vector vector ComplexMath std complex double amp A const int n vector int amp vector ComplexMath std complex lt double gt gt amp B Gaussjc h COMPLETE This is a simple Gauss Jordan solver based on Numeri
44. jugate Outputvar 1 10 X component of the intensity 1 T Pu V Outputvar 1 11 Y component of the intensity 1 Outputvar 1 12 Z component of the intensity Sonitus BE Page 33 of 58 SonitusBE e THE SOLUTION METHOD FOR A MATRIX PROBLEM Two numerical methods are available for the solution of an matrix problem however the computer code has been written so that a different solution method may be substituted with no significant alteration required The two methods are Gaussian Jordan and LU decomposition both validated through the use of complex value matrices The method is based on the classical equation Ao where the vector of boundary conditions and sound source contributions are given by known the amplitudes of the fundamental solutions o unknown and the square matrix of boundary integrals A known In the code these variables are copied into duplicated variables for processing The solution function is then called which replaces the square matrix with A and the vector with o Details of solution methods are available in Flowers Introduction to Numerical Methods in Clarendon Press Oxford 1995 The Gauss Jordan elimination is not the most efficient solver however is well known and many varieties of the code exist It differs from a standard Gauss solver as back substitution to reduce a matrix to a triangular form i
45. ley 1998 Flowers Introduction to Numerical Methods Clarendon Press Oxford 1995 Sonitus BE Page 55 of 58 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 BEM 11 boundary 11 13 16 17 18 19 20 21 22 26 28 36 39 40 41 42 43 44 47 52 CAD 11 13 complex 13 17 39 43 44 45 46 47 52 Data recovery mesh 18 48 discretisation 13 22 domain 11 13 16 17 18 19 20 21 22 24 35 46 48 49 51 Element 1 11 16 20 40 43 44 52 54 error 13 21 37 40 41 51 Excitation frequency 17 finite element 11 13 19 Fundamental solution 17 Gauss Jordan solver 46 Green s function 14 17 21 46 Helmholtz 17 18 21 22 impedance 13 41 42 47 52 Input 33 35 41 52 Interpolation 37 38 Jacobian transformation 44 SonitusBE Linux 11 12 LU decomposition 46 MATERIAL 41 Matlab 11 33 35 monopole 17 18 19 46 Node 16 40 43 44 45 Numerical Recipes 37 45 46 Output 35 36 phase 18 42 52 potential 14 precision 37 39 43 pressure 11 17 18 19 20 33 35 42 44 46 47 49 50 quadrilateral 16 35 37 38 41 47 49 52 reference intensity 41 reference pressure 41 reference work 41 Requirements 12 14 15 33 35 36 S
46. mp zeta2 ComplexMath std complex double CinterpQUAD const ComplexMath std complex lt double gt 851 const ComplexMath std complex double amp N2 const ComplexMath std complex lt double gt amp N3 const ComplexMath std complex lt double gt amp NA4 const double amp zetal const double amp zeta2 Messages h COMPLETE 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 SolutionChoice h COMPLETE 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 Sonitus BE Page 40 of 58 SonitusBE aa conditions If returned integer is 1 the direct solver
47. necessary to generate u The constants can be taken inside this variable and summarised as u woar pw aQ An alternative formulation is clearly gained by writing the sum of the derivatives of the fluxes as equal to zero If q q 0 then the boundary integral equation becomes a uar w 4 pw aQ 0 4 pu aQ 0 Writing the dipole term u u relates the value of in the domain to the amplitude of a series of dipoles distributed over the whole surface u q udr bi aQ Q The former method has greater advantages as lower orders of derivatives of the fundamental solution are utilised If we choose this former method then we can write and the derivative the flux as Sonitus BE Page 24 of 58 SonitusBE AN u w oar bu U Zoar bi aQ r r This applies point inside of domain however in order to find value of a point boundary the location i needs to be moved to the surface A singularity exists therefore the integration is carried out analytically to obtain u woar bu aq 4 2 7 q odr bu aQ r 2 DISCRETISATION OF THE SURFACE DOMAIN Up until this point the boundary integral equation has been considered as a continuous integral over the surface in the domain There are a number
48. ng making them a popular first choice for initial investigations SonitusBE N The boundary element is closely related to the finite element method and exploits the principle that in acoustic predictions it is often the case that the governing transmission medium is non dispersive and isotropic The method shares some similarities to the finite element method in that it can utilise information from a computer aided engineering process however the key advantage is that rather than discretising and meshing the volume only the surfaces need treating This leads to a reduction in computational effort and memory for a given problem As the surfaces only need meshing it is relatively easy to join dissimilar domains through common boundary conditions The boundary element method requires the use of a fundamental solution governing the propagation of sound in free space applied as a weighting function to every node on the surface mesh This fundamental solution takes the form of a Green s function found in theoretical analysis This usage of ideal propagation solutions provides the potential for extremely accurate solutions IV NOMENCLATURE GREEK NOTATION EXPLANATION VECTOR OF SOURCE AMPLITUDES OF THE ERE FUNDAMENTAL SOLUTIONS TD SRE 0 __ 6 BRROR _ _ 59 250 6 NONDINENIONALCOORDINNIE A WAVELENGTH UNKNOWN AMPLITUDE OF THE FLUXES OF THE FUNDAMENTAL SOLUTIONS DENSITY UNKNOWN AMPLITUDE OF THE FUNDAMENTAL SO
49. non parabolic surface can be defined by three distinct nodes for a triangular element or four for a quadrilateral element Through a regular distribution of node points the elements are evenly distributed and spaced over the entire surface of the geometry Elements may share common nodes and it is important that in the definition of the nodes creating all elements that the normal directions are consistent This is shown on the figure as the node connections all pointing clockwise or anticlockwise Ni Vice N x yas Z4 N x y 2 N x Z2 N3 Z4 23 N UX L 5525 Figure 2 triangular element characterised by three independent nodes and a quadrilateral element characterized by four Note that the direction normal of all elements must be consistent Sonitus BE Page 16 of 58 cu SonitusBE e The density of the mesh is equivalent to the number of elements for a given surface area The mesh should be smooth and the density approximately constant Excitation frequency In order that a non zero sound pressure field is calculated the non trivial solution the geometric domain must contain at least one sound source otherwise the only sound field that can exist is silence This monopole source is located at an arbitrary position The boundary element calculation routine calculates a predicted sound field for each discrete excitation frequency up to a finite upper limit There ar
50. nts first the value of the fundamental solution is obtained between each node i and nodes which make up the three corners of the triangular element 1 2 3 Figure 10 Triangular linear integration method The position of a point on the triangular element is r xi 16 L6 L 8 S pp 1 L L These are rearranged to obtain any point on the triangular element as a function of non dimensional coordinates 69x 0 6 y 6 6y tk L 0 D x 0 2 05 65 6 x A S 1562 z a xe Sulle 3 u G T US 9274 Iz Sonitus BE Page 29 of 58 bam Once the values at an arbitrary location on the triangular element are known the integration over the surface area can be found using a simple Gaussian integration method for triangular elements where the non dimensional coordinate varies from 0 to 1 SonitusBE as u J G fwar 1 2 16 Je 92492 2 Oz 06 06 Oz Ox Ox ag Ox Ox BU ag 22 1 1 D 26 ma 11 Lm Wim 1 2 3 1 3 2 1 6 1 3 This is Gauss Integration method for linear elements in two dimensions non dimensionalised See Abramowitz and Stegun P887 P916 for weightings So i C dede HIE 1 6 f 1 6 1 6 71 6413 where f 6
51. om 1 to 1 This provides the derivative at 1 the values of and double DZETAlinterpQUAD const double amp 1 const double amp N2 const double amp N3 const double amp N4 const double amp zetal const double amp zeta2 Define a set of weighting functions y V W equated to y V 1 4 001 6 W3 1 4 V Y 4 Interpolated return value y N 0 e Interpolation routine for a quadrilateral element with nodes N to and non dimensionalised surface coordinates and which vary from 1 to 1 This provides the derivative 4 at 2 the values of double DZETA2interpQUAD const double amp 1 Sonitus BE Page 39 of 58 SonitusBE e const double amp N2 const double amp N3 const double amp 4 const double amp zetal const double amp zeta2 Define a set of weighting functions W3 W equated to y 1 4 V 1 4 V 1 4 V 1 41 Interpolated return value v4 N v N V N 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 amp N1 const ComplexMath std complex double amp N2 const ComplexMath std complex double amp N3 const double amp zetal const double a
52. poQu V oft ze r C Note Determine if it is normal velocity only or magnitude of velocity which is zero Difference is rx Irl type term but in the unit normal direction of the element The alternative boundary condition is a type 3 impedance condition where the complex value of normal pressure to normal velocity is specified Although this can also be written in terms of a fundamental solution it may be shown that these terms cancel Fundamental solution boundary element method relies use of fundamental solutions to applied to each node as a weighting factor In the case of acoustic predictions the fundamental solutions take the form of the Green s function solution to the Helmholtz equation This represents the propagation of a wave through free space and can be applied in either two dimensional or three dimensional forms the two dimensional form utilises Hankel functions rather than exponential functions Further details are provided in the section on monopole sound fields Sonitus BE Page 17 of 58 SonitusBE HELMHOLTZ FUNDAMENTAL SOLUTIONS MONOPOLES E exp kr r Table 5 Fundamental Solutions for acoustic wave propagation in free space These Green s functions represent how a wave would propagate through the free space between the source position and observer location a distance r away It is important to understand that the fundamental solution
53. r dry air at a temperature of twenty degrees Centigrade this equates to approximately 343m s SonitusBE e Surface The surface is a boundary between interior and exterior solution domains The surface must be a continuous function in both displacement and slope the method can fail should sharp corners exist in the geometry Node The surface is discretised into a sertes of evenly spaced points Each point on the surface is defined by three coordinates x y z where the join between nodes are assumed to be a linear segment It is important that when discretising the surface the nodes are spaced evenly across the whole domain with the spacing dictated by the maximum excitation frequency of interest Thus for calculations there should be at least four nodes per acoustic wavelength to avoid aliasing while for interior domains the recommended minimum number rises to eight nodes per wavelength The spacing between nodes does not need to be the same only that the macro spacing is similar N X Yiz Figure 1 A node characterised by three independent coordinates Wavelength The acoustic frequency f is related to the acoustic wavelength A by the speed of sound c fA It is assumed that the acoustic medium is an isotropic material where shear forces cannot be sustained Element The nodes are just points in Cartesian space however by joining three or more unique nodes together a linear flat surface is created This
54. s avoided Instead additional operations reduce the matrix to diagonal formation It is assumed that the matrix A is non singular with a well defined inverse VII INPUT FILE DEFINITION The input files can be generated automatically using the software Sonitus Pre BE 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 Sonitus BE Page 34 of 58 SonitusBE Figure 12 Generation of 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 Sonitus BE Page 35 of 58 SonitusBE ees DESCRIPTION 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 Pre BE geometry creator Reference values for the pressure intensity and work Nodal numbers and positions 2 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 cre
55. s document two integral methods are shown as examples relating to the integration of a quadrilateral and also triangular element respectively Both element types are calculated in a similar manner The components of the matrix Ali G dV where the fundamental solution is assumed to be Dj the form u i uL exp icr c where the frequency in radians second is related to the frequency in Hz Ar by 2 and is the distance between the source element and the observer element The element is defined using non dimensional coordinates and an interpolation routine is used to obtain the values of the fundamental solution across the whole element A linear integration method is then applied to determine the value of the integration over the area of the element QUADRILATERAL ELEMENTS The principle axes are x y z however the linear element may be sloping in directions which are different to this Therefore there needs to be an integration along non dimensional coordinates This involves a change of variable in the double integral using a Jacobian method The change of variable in a double integral using a Jacobian can be summarised from Kreyszig as b b f G dx du 212135 1562 1262 d d 2 e AS Ox 206 Ela 24 90 G wictacac 62 16 Jg 92 g2 Sonitus BE Page 27 of 58 SonitusBE scs We now look a
56. t the element j with four labels denoting nodes k 1 2 3 and 4 Element j u k 2008 2 3 G k 3 i ro 4 0 4 7 Node i e Figure 8 Quadrilateral boundary integral surface showing non dimensional coordinate form of the surface direction For each coordinate node on the element the product u 1 obtained where the surface directions are obtained through an interpolation routine the details are provided in the coding explanation _ Oz amp P 2 Ox 4 7710505 0500 _ Ox Ox 5 2222 24 22 Once the product of is found at each of the corners the numerical integration can take place For a linear surface the following is appropriate f G 1 1 y 21e 1 Ls J 9 12 1 Jib Le e Figure 9 Quadrilateral linear integration method 1 1 2 yz gt Y ww f Cari ES 1 Lm d w i m 1 1d 2 Sonitus Page 28 of 58 This is Gauss integration method for linear elements in two dimensions non dimensionalised See Abramowitz and Stegun P887 P916 for weightings SonitusBE So FCU 8 8 1 3 f 45 45 f 5 1 3 where 1 i 2 i5 u G ou This is the final boundary integral for the linear quadrilateral element TRIANGULAR ELEMENTS The process is similar for triangular eleme
57. their associated amplitudes F Ao At this point we assume that we have a suitable solution for the amplitudes o which will allow all physical parameters to be determined inside the domain If the amplitudes of the fundamental solutions are known across the surface as shown in Figure 11 then the corresponding boundary integrals can be calculated between each data recovery node and each boundary element on the surface This is written as Sonitus BE Page 31 of 58 SonitusBE Surface nodes fe A Surface 9 N j ks K E B 3 YN x e 9 ET j MER Data recovery nodes and mesh j 16 Figure 11 Boundary integrals from the surface to the data recovery mesh In order to obtain the pressure and velocity at each data recovery point the relationship between the fundamental solution and these physical relations are re stated _ 1 Azr p iQpou Qu V2Q 9 or The pressure may be calculated at each of the data recovery points thus iP Fp as the amplitude of the source usually Q is given by In addition to this the contribution from the source location is also added The output variables relating to the sound pressure are denoted for the r th node Outputvar 1 0 to Outputvar i 2 Outputvar 1 0 Sound pressure in dB re Pref T Pane
58. then described in section VI specifically the indirect boundary element method the method of discretising the surface of the geometry in the acoustic domain the mathematical routines used to integrate a surface element both quadrilateral and triangular The source location and amplitude are then described as are the data recovery parameters A numerical routine for the solution to the matrix problem is briefly summarised As the Sonitus BE program is a stand alone console based program the input and output file definitions are described in sections VII and VIII respectively An overview of the program code is provided so that a user may understand the organisational hierarchy of the program section IX functions and subroutines section X Sonitus BE Page 11 of 58 S A number of sample files are provided for reference and are detailed in section XI Notes on the program are found in section XII A summary of this manual is provided in section XV SonitusBE e REQUIREMENTS The requirements are highly dependent on the geometry size mesh density and frequency range which is required The program runs as a standalone application without graphics interface and can therefore run either as a background activity or offline monitored through a network connection on a server The recommended starting points are shown in Table 1 For pre processing requirements Table 1 provides useful guidance however for other commer
59. tries into the correct surface mesh format through the use of bespoke program One such example is provided here in the subdirectory DRMGen The Patran file format contains information about the location of the nodes making up the data recovery mesh and also how they are connected to create a structured orderly design for graphical processing A file may also include data on material properties and node loading however these are not of interest for this application The scripted mesh generator provides a structured grid of quadrilateral elements with a fixed number in the x and y directions It models a mesh surrounding an automotive tyre see Figure 15 for an illustration The first few lines of the file describe how many nodes are to be defined in the program and a date and time of creation The nodal data number is defined followed by it s position in three dimensional Cartesian coordinates x y z Following this a list of the elements are defined These quadrilateral elements connect four nodes together As previously described it is vital that the normal directions of conjoining elements be consistent The C data recovery mesh program be compiled using the command g DRMGen Tyre cpp a out Once the location of the data recovery nodes are known it remains to determine the physical parameters of interest We start with the equation for the known boundary condition relating to the matrix of boundary integrals and
60. u G This l is the final boundary integral for the linear A SOURCE DEFINITION AND LOCATION For a non trivial solution the domain needs to include at least one source with defined amplitude Q and location If the distance between the source and the observer element is given by Rs then the pressure and velocity may be written in terms of this distance and the fundamental solution u exp ioR c 4 p iQoou V Q0 9 or DATA RECOVERY PARAMETERS The data recovery points are incorporated into the input file as a list of points where the user is interested in knowing the acoustic pressure intensity or velocity They must be at least one quarter of an acoustic Sonitus BE Page 30 of 58 wavelength away from the solid geometry boundary in order to avoid singularities occurring in the boundary integral formulation SonitusBE e The data recovery points can either be designed as single points or through the use of a meshing program which allows a fuller graphical representation of the end data The preferred format for the data recovery mesh is through the use of a Patran Neutral file with extension out 1 Although Patran is a commercial program many converters exist for moving between one neutral format to another As a result it is relatively easy to obtain the data recovery mesh in this appropriate format however it is also possible to directly program typical geome
61. us BE Page 49 of 58 SonitusBE oss ss Free field solution The free field 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 Sonitus BE Page 50 of 58 SonitusBE Qn P dB F Hz 100 1 128 05 N 126 0 1 2 0 us 124 Y mF 0 1 0 5 X m P dB F Hz 400 at E 140 0 5 N 138 0 1 0 5 0 136 Y m 0 1_ 0 5 X m P dB F Hz 1300 152 1 2 150 05 N 148 22 05 Y mF 0 1 0 5 X m Figure 16 Sound pressure for a source in the free field note the scale of the amplitude is different at each frequency Sonitus BE Page 51 of 58 SonitusBE Sound pressure domain with solid geometry The solution for a sound source in a free field with no solid geometry leads to a spherically propagating
62. utines are discussed and function calls detailed ErrorHandle h COMPLETE 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 amp mess message mess Interp h COMPLETE type functions and quadrilateral elements This file contains the routines which interpolate data from numerical recipes It exists to contain all of the routines which allow interpolation on triangular 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 978 0 486 61272 0 The following routines are provided based on an interpolation of double precision numbers e Interpolation routine for a triangular element with nodes N to and non dimensionalised surface coordinates and which vary from 0 to 1 double interpTRI const double 651 const doubl const doubl const doubl const doubl e e e e amp N2 amp N3 amp zetal amp zeta2 Define a set of weighting functions y v Y3 equated to y V c V4 1 6 6 Interpolated return value y N y N y N e Interpolation routine for a triangular element with nodes
63. would be chosen if 2 the solver is indirect formulation int SolutionMethod void Verification h COMPLETE 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 amp NumberofNodes const vector int amp NodeNum const vector MeshTypes Node double gt amp NodeCoord 4 const vector MeshTypes Element int gt amp ElemStruct Em const vector int amp ElemNum const int amp 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 15 determined by the order in which the nodes are arranges either clockwise or counter clockwise 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 quadrilater

Sonitus BE R1 User Manual

Contents

Download Pdf Manuals

Related Search

Related Contents