This chapter introduced basic one-dimensional modeling techniques for applications in linear acoustics. From the governing equations based on the principles of continuum mechanics, a unidimensional modeling framework was deduced. Intrinsic and wall losses as well as radiation effects from different openings were incorporated into the model. As an example application, the equivalent acoustic circuit of the resonator of an open cylindrical flue pipe was constructed. It was demonstrated that the one-dimensional description can provide a lot of useful information on the acoustical system at hand by means of a few straightforward calculation steps.
The limitations of the one-dimensional model are that (1) it is only capable of handling sim-ple geometries, and (2) it is only applicable under the cutoff frequency of the system. In order to overcome the limitations of one-dimensional models, numerical techniques can be applied for the simulation of two- or three-dimensional systems. Some of these techniques are discussed in the next chapter.
Chapter 4
An introduction to finite element methods in acoustics
This chapter introduces the numerical techniques used in the thesis for three-dimensional acous-tic simulations. Using numerical models have various advantages, some of them are listed below.
1. By means of numerical simulations approximative solutions of problems whose analytical solution is not known—or too cumbersome to evaluate—can be attained. The error of these approximations only depend on the accuracy of the underlying physical model and the computational capacity available at hand.
2. By hybrid
models decomposing a problem into sub-problems of different complexity, results of numerical
simulations can also be incorporated into analytical models. The resulting hybrid models lead to an advantageous combination of accuracy and computational performance.
3. Numerical techniques can also be used in order to assess the validity of analytical approxi-mations of the solutions of certain problems.
The aim of this chapteris to summarize the most important features of the techniques and objective their implementation in order to facilitate a better understanding of their applications. The
dis-cussion, however, is presented as a general framework and is not limited by the specific field of application utilized later on in the thesis. Explaining the background of numerical analysis methods and discussing each technique in detail are certainly not the objective of this chapter.
For more explanation and in-depth analyses the reader is referred to the cited sources.
As far as three-dimensional acoustic simulation is concerned, there are quite a few options to choose from, the most common being theboundary element method (BEM), the finite element
method(FEM), and thefinite difference time domain method(FDTD). Among choice of the FEM these the FEM seems
to provide the most flexibility for the applications discussed in the dissertation because of the following reasons. (1) It can easily handle elements of arbitrary shape and arbitrarily varying material properties, unlike the BEM or FDTD. (2) The same formulation can be used both in the time and the frequency domain, which is not the case neither for the BEM, nor for FDTD.
(3) By means of the extensions introduced in Sections 4.3 and 4.4 the FEM can be applied for unbounded and semi-infinite problems. Since this versatility is exploited in the studies presented in Chapters 7 and 8, the FEM seems a reasonable choice for solving the corresponding problems.
Thischapter first presents the finite element formulation of the Helmholtz equation and the structure corresponding boundary value problem in Section 4.1. The problem of open boundaries is
ad-dressed in Section 4.2. The following two sections discuss two different approaches for finite element simulations with open boundaries. The infinite element method is introduced in Sec-tion 4.3, whereas the perfectly matched layer technique is presented in SecSec-tion 4.4.
37
4.1 Finite element formulation for acoustic problems
This section presents the general formulation of the acousticfinite element method(FEM) using the Galerkin variational procedure. The finite element method is a general technique for obtaining the approximative solution of boundary value problems of partial differential equations. Although the techniques presented in the sequel can be applied for various problems, the methodology is discussed using the acoustic Helmholtz equation.
4.1.1 The boundary value problem for the Helmholtz equation
Our discussion starts with the homogeneous Helmholtz equation (3.21) for the sound pressure p(x, ω)ind-dimensional space, that is
Helmholtz
equation ∇2p(x, ω) +k2p(x, ω) = 0 x∈Ω⊆Rd. (4.1)
It is assumed that the domainΩis bounded by the boundary Γ. In order to get a well-posed problem the boundary conditions for the partial differential equation must be known in all points ofΓ. For the Helmholtz equation (4.1) three types of boundary conditions can be defined.
Dirichlet boundary conditions also known as first-type boundary conditions for the Helmholtz equation define the pressure on one part of the boundaryΓp⊆Γas
Dirichlet
condition p(x, ω) = ¯p(x, ω) x∈Γp. (4.2)
Neumann boundary conditions also called second-type boundary conditions for the Helmholtz equation define the normal derivative of the pressure, or equivalently the normal particle velocityvnon a part of the boundaryΓv⊆Γas
Neumann condition
∂p(x, ω)
∂n = jωρ0v¯n(x, ω) x∈Γv. (4.3) Robin boundary conditions are obtained as the linear combination of Dirichlet and Neumann boundary conditions. The Robin boundary condition prescribes the specific impedanceZ relating the pressure and its normal derivative on a subdomain of the boundaryΓz⊆Γ.
Boundary conditions must be specified on the whole boundary, such thatΓp∪Γv∪Γz = Γ. With the boundary conditions defined, the uniqueness of the solution to the Helmholtz equa-tion (4.1) is ensured.
boundary value problem
The differential equation together with its boundary conditions form a boundary value problem.
4.1.2 The weak form of the boundary value problem
Theweak formof a differential equation is its transformation into an integral equation. The weak form is obtained by multiplying the PDE by an arbitrary test functionφ(x)and integrating the result over the whole domainΩ. The integral form is called “weak” since—after integration by parts—it defines less strict criteria for the solution than the original PDE, see e.g. [151, chapter 3].
For the homogeneous Helmholtz equation (4.1) this gives Z
Ω
φ(x)∇2p(x) dx+ Z
Ω
φ(x)k2p(x) dx= 0. (4.4) For the sake of convenience the notation of the dependence onωis omitted hereafter.
4.1. FINITE ELEMENT FORMULATION FOR ACOUSTIC PROBLEMS 39 Making use of the identity that ∇ ·[φ(x)∇p(x)] = ∇φ(x)· ∇p(x) +φ(x)∇2p(x), the first integral can be extracted into two parts, which gives
Z The right hand side can be transformed into a surface integral using the divergence theorem:
Z
Expressing the normal derivative of the pressure from the Neumann boundary condition (4.3) in the right hand side, and multiplying byρ0c2gives the final equation of the weak form:
Helmholtz
To arrive at a system of equations with a finite number of unknowns, the weak form must be discretized. The discretization is done by representing the solution and test functions—p(x)and φ(x)for the Helmholtz equation—by a finite space ofshape functions. When the approximation is substituted into the original integral statement—(4.4) in our case—the result is a weighted
integral of the residual error over the test function space. Thus, weighted residual this approximation can also
be called themethod of weighted residuals. There are different approaches for choosing the shape function spaces, see e.g. [151, chapter 3] for details. Here only the most commonly used Galerkin procedure is introduced.
In the Galerkin variational method both the solution function and the test function are approx-imated by the linear combination of a finite set ofelementary shape functionsof the same function space. The approximation using the shape functionsNj(x)is given as
p(x)≈
In the above notation the vector of shape functionsN(x)is defined as a row vector, whereas the
vector of test and solution weightsΦand p, respectively, are column vectors. The number degrees of freedom of
shape functions (and corresponding weights)nis the number ofdegrees of freedom(DOF) of the system. Similarly, the gradient of the above variables is attained as
∇p(x)≈
with∇Ndenoting ad×nmatrix,(∇N)ij representing the derivative of thejth shape function with respect to theith coordinate∂Nj/∂xi.
In the following the Galerkin formulation of the weak form of the Helmholtz equation is dis-cussed. Substituting equations (4.8) and (4.9) into the weak form of the Helmholtz equation (4.7) and approximatingv¯n(x)also by the same shape functionsNj(x)asp(x)andφ(x)we get
Finally, taking the spatially independent terms out of the integrals, and noticing that the common multiplierΦTcan be omitted, the Galerkin weak form reads as
Galerkin Equation (4.11) is usually given in its matrix form as
matrix form Kp−ω2Mp= jωAv. (4.12)
In the above equation the matrices
system Then×ntype matricesKandMare referred to as the acousticstiffnessandmass system matrices, respectively. From equation (4.13) it can be seen that the resulting matrices are hermitian, which is an advantageous property of the Galerkin method, since such matrices require less storage area in memory and enable the usage of faster matrix solver routines [129, p. 153].
4.1.4 Spatial discretization
Spatial discretizationrefers to the process in which the solution domain is decomposed into a finite number of non-overlapping subdomainsΩkas
domain
It should be noted that the approximation in the above equation means that the numerical so-lution domain is not necessarily congruent with the original domain; however, the error of the approximation is usually negligibly small.
The subdomainsΩkare calledelementsin the finite element framework, and they are defined bymapping transformationsfrom finite parent domainsOe. The coordinate transformation from localξ∈ Oeto globalx∈Ωcoordinates is given as withLi(ξ)denoting the mapping functions, also referred to as geometry shape functions, andxi
being theithnodeof the element, which is defined by a total number ofnenodes and the corre-sponding mapping functions. As a shorter notation the row vectorL(ξ)contains the geometry shape functions and thene×dmatrixXis composed of the spatial coordinates of the nodes.
It is usually useful to define the shape functionsNj in the local coordinate systemξ. Then, by
matrix assembly
applying the domain decomposition defined in equation (4.14) and using the element map-ping (4.15) the system matrices in equation 4.13 can be assembled in an element-by-element man-ner with evaluating integrals only over the parent domainsOe. As for one element type the parent domains are the same, the integration can be performed using simple numerical quadra-ture rules. In some special cases even analytical integration is possible. Furthermore, with the proper—and usual—choice ofNj(ξ), the system matrices become sparse, which significantly re-duces the computational effort required for the solution.