Bild von Institut mit Institutslogo
homeicon uni kontakticon kontakt
unilogo Universität Stuttgart
Boundary Element Methods 

Examples and Demonstrations


Examples and demonstrations in MATLAB



Besides the theoretical foundations of boundary element method and some of its numerical aspects, the BEM serves as an engineering tool for physical problems overcoming the restriction of an analytical solution. In this BEM class, Laplace’s equation is an intensively studied governing equation solved by using a BEM formulation. The physical background of this equation is in this case the stationary heat transfer equation. Adequate Dirichlet and Neumann boundary conditions have to be applied to state a well posed problem – if this is the case, we can implement the BEM to solve for the solution of the actual problem. If source terms are added to Laplace’s equation, the so called Poisson equation evolves. For example, Poissson equation models heat transfer taking into account heat sources or sinks respectively  inside the domain.

In the course homework project, a BEM program code in MATLAB is implemented to solve Laplace’s equation in 2D. The details are announced in class and not part of this sections. However, having formulated this BEM program code, it is relatively easy to extend it in such a way to obtain results for the solution of Poisson equation.

A further topic in this class is called “BEM in Acoustics”, and some modifications of the program code  - as for example another choice of the fundamental solution and the introduction of complex numbers – enable us to use it for Helmholtz equation.

On this page our course homepage, we want to demonstrate some of the results for these three equation, which in case of Laplace’s and Poisson equation both strongly refer to problems covered in class hours.

Concerning Helmholtz equation, there are results for the Helmholtz resonator as an example for an interior problem.

Sound radiation from both rectangular domain give insight into exterior problems which are surely the stronghold of BEM as compared to FEM.

For all of the problems, linear interpolation functions have been used along with a double nodes concept for the representation of the Neumann variable on the boundary.


Laplace’s equation


In all of the following examples, the non dimensional heat transfer equation is considered as an example for Laplace’s equation


We are dealing with a rectangular domain as depicted in Figure 1. There are eight equal elements used implying eight nodes for the temperature and sixteen for the interpolation of the flux according to the double nodes concept. The nodes are represented as blue circles, whereas the dots represent Gauss point coordinates for the numerical integration.

Figure 1: Meshing of the rectangular domain. Geometrical Nodes (o) and Gauss points.


First, the boundary conditions are

They correspond to the example from class.

As result, we get the temperature field from Figure 2.


Figure 2: Temperature field for given BC.


Due to symmetry of both domain and BC, the (analytical) solution is not dependent on spatial coordinate x  and is just a linear function of y .  Figure 3 shows the same information as a contour plot. Besides, the flux vectors indicate the direction and the quantity of the heat transfer. The calculation of the flux is done by finite differences, therefore the results get inaccurate on the boundary.


Figure 3: Temperature and flux field


The next plot in Figures 4 shows the solution of the temperature field for the boundary conditions from the course project as a contour plot.

The boundary conditions

are applied, leading again to a linear solution as a function of  and . In the interior of the domain, the temperature is evaluated using the integral equation whereas values on the boundary are interpolated from the nodal solution. The flux is determined numerically by the gradient of the temperature field, depicted by blue arrows. Again, the boundary element error gets obvious.



Figure 4: Temperature and flux field.


This boundary element error is reduced by introducing smaller element sizes. The boundary of the rectangle is now divided in sixteen equal elements (Figure 5).



Figure 5: Finer mesh of the domain.


Again, the same boundary conditions are imposed, and if Figures 4 and 6 are compared, it gets obvious that the boundary element error has been reduced a lot.


Figure 6: Temperature field with a finer mesh size, compare with Figure 4.


The two problems that have been presented so far both possess linear solution. Finally, a problem is presented where we cannot expect the solution to be a linear function of spatial coordinates x and y .

The following BC are considered:

This means that on the right boundary on the bottom part, some strong heat flux is entering the domain. Besides, flux can only enter at the top and bottom boundary, since the rest of the boundary is fully isolated. The evaluation is performed for the interior of the domain, giving the temperature profile in Figure 7 and the contour plot in Figure 8.


Figure 7: Temperature field for given BC


Figure 8: temperature and fluxes for given BC.


Poisson Equation


If Laplace’s equation is extended by a inhomogeneous part , one obtains Poissons equation


In the following consideration, is considered to be a constant. If is positive, it is representing a heat sink, whereas a negative value stands for a heat source. If the boundary conditions as in the first examples

are applied again, and  ,

one obtains the temperature field in Figure 9. The solution is not linear as a function of y , but behaves in a parabolic way. It gets clear from this solution that a positive takes out energy of the system.


Figure 9: Interior temperature field for  (sink).


If , heat is introduced into the domain of the system, and therefore the temperature level increases. (Figure 10)


Figure 10: Interior temperature field for  (source).


Helmholtz equation


The BEM finds one of its major application in fluid acoustics, where time harmonic behavior is described by Helmholtz equation


where is the acoustic (excess) pressure and the acoustic wave number. In general, there are two classes of acoustic problems referring to the domain. One of them has a domain which is enclosed completely by its boundary, called an Interior Problem. In the other class, the domain is unbounded, thus going to infinity. In this presentation, examples of both classes are shown in 2-D. A 2-D application stands also for a domain with a constant thickness and corresponding boundary conditions, such that the 3-D physics can be reduced to two dimensions.


Interior Problem


The so called Helmholtz resonator just represents a tube with rigid walls and certain boundary conditions at the two ends. Even though it can be threatened as a one dimensional problem, it is also possible to calculate it in two dimensions. The domain is shown in Figure 1 as well as nodal points. Between nodes 1-57, a rigid wall is assumed, leading to Neumann boundary conditions which are equal to zero. Between nodes 57-1, the system is excited by applying Dirichlet BC with the amplitude equal to one. Note that this amplitude is normalized. The excitation frequency is .



Figure 1: Domain of the Helmholtz resonator


The medium in the domain is considered to be air. As a result, a sinusoidal   standing wave builds up in the domain with maxima and minima (Figure 2). The gradient at the rigid wall ( ) is equal to zero. The peaks have an amplitude of about 1. In the contour plot (Figure 3), the pressure gradient is shown as well. Maxima (red and blue) and nodes (yellow) are a result of interference of incident wave and its reflection at the rigid wall. When the frequency is lowered to 100Hz, one obtains the solution from Figures 4 and 5. One of the resonance frequencies is at 99Hz, thus we are operating in the close vicinity of that, and therefore the amplitude increases from about 1 to nearly 15. All of these numerical results in BEM have been verified by comparing its solution to the analytical one. However, a discussion is beyond the purpose of this presentation.


Figure 2: Acoustic pressure in the Helmholtz resonator at 200Hz


Figure 3: Acoustic pressure in the Helmholtz resonator at 200Hz


100Hz (Resonance at 99Hz)


Figure 4: Acoustic pressure in the Helmholtz resonator at 100Hz


Figure 5: Acoustic pressure in the Helmholtz resonator at 100Hz


Exterior Problem


Now a quadratic domain is considered (Figure 6) as well as element nodes (o).

Figure 6: Domain for the exterior acoustic field and BEM discretisation.


Neumann BC are applied on the whole boundary as follows: On the upper, lower and right side, a rigid wall is assumed thus leading to zero boundary conditions. The left wall vibrates with a normalized amplitude of one with different frequencies. Therefore, an outgoing wave is generated and propagating mainly to the left. However, the wave starts diffracting around the corners and continues finally in all directions, even though only a very low wave is propagating to the right. The blank rectangle is at the position of the excluded domain.

In the following pictures, the acoustic pressure field and the intensity field is depicted in the figures for different frequencies. It can be seen easily that the near field distance is increased by increasing the frequency. The near field characteristics are underlined by the numerous side lobes occurring at high frequencies ( ). Above 300Hz, the number of elements has to be increased in such a way that sufficient results are obtained.


Figure 7: Acoustic pressure field at 100Hz




Figure 8: Intensity field at 100Hz




Figure 10: Acoustic pressure field at 200Hz


Figure11: Intensity field at 200Hz


300 Hz


Figure 12: Acoustic pressure field at 300Hz


Figure 13: Intensity field at 300Hz


520 Hz

Figure 14: Acoustic pressure field at 520Hz with 16 elements per side length.



Figure 15: Intensity field at 520Hz with 16 elements per side length












Figure 16: Intensity field at 900Hz with 32 elements per side length.


Figure 17: Intensity field at 1800Hz with a 0.05 spatial field point resolution and 64 elements per side length.


Time animations - Shallow water waves


Solutions of Helmholtz equation render a picture for one moment of time. The solution can be time animated by multiplying the complex solution with the formerly separated harmonic time behavior. In those animations, a blue color map has been chosen because Helmholtz equation

represents a model for shallow water waves as well, where is the actual water depth, the mean water depth and the gravity acceleration.

As in the plots above, one of the four side boundaries is vibrating harmonically.




Click on the links to start animations

                  Wave at 100 Hz        

                  Wave at 200 Hz        

                   Wave at 300 Hz