PROCEEDINGS of the
23rd International Congress on Acoustics
9 to 13 September 2019 in Aachen, Germany

Integration of FEM shell elements as a “boundary condition”
in BEM calculations using different solution methods
Ralf BURGSCHWEIGER1; Ingo SCHÄFER2; Martin OCHMANN1
1

2

Beuth Hochschule für Technik Berlin, University of Applied Sciences, Faculty II,
Research Group Computational Acoustics, 13353 Berlin,

Bundeswehr Technical Center for Ships and Naval Weapons, Maritime Technology and Research(WTD71), Acoustic modelling (GF640), 24148 Kiel

ABSTRACT
Within the framework of research projects in the field of "Computational Acoustics", a powerful BEM-based
code for the determination of the backscattered sound pressure level in the far field has been developed inrecent years, using different solution methods (direct and iterative solver, fast multipole method, ray tracing
method etc.).
This code has been extended so that in addition to the "classic" boundary conditions (pressure, normal
velocity, impedance, inertial coupling, etc.), additional FEM shell conditions and thus elastic material
properties can be taken into account. For this purpose, an additional FEM equation system is set up for theelements concerned which is integrated directly into the BEM equation system via appropriate transformation
matrices.
The solution of this overall system can be done by means of direct equation solvers (based on a complete
sparse matrix), by eigenvalue calculations or by the so-called Schur method (using iterative solvers).
The paper deals with the mathematical and physical fundamentals, shows the differences between the solution
methods of the entire system and presents the results of test structures. As far as available until then, acomparison is also made with results from "pure" FEM calculations.
Keywords: FEM shell elements, BEM

1. INTRODUCTION
The calculation of the sound pressure in the far field for more complex structures of thin elastic
shells, which are in liquid media and interact with them, is done for the low and middle frequency
range mostly using the Finite Element Method (FEM). However, this is often no longer feasible withjustifiable effort, especially for objects with larger dimensions.
Therefore, an already developed powerful BEM application has been extended to the possibility ofdefining thin shells as a "boundary condition" in the form of FEM shell elements, without the needfor changes to the surface structures used.
This takes into account the elastic material properties of the thin shells and binds them directly
into the BEM, without the need for additional applications or calculation steps.

2. BASICS AND FORMULATIONS
2.1 Used FEM shell elements and matrices
The presented method uses simple triangular planar FEM shell elements, which are assigned anelastic material and a shell thickness. The evaluation of the associated physical quantities(displacement ܹ, rotation angles ߮, forces ‫ )ܨ‬takes place at the "corners" (nodes) of the elements.
For the calculation of the element-specific mass matrices ‫ۻ‬௩ and stiffness matrices ۹ ௩ , theKirchhoff method, as described in the formula work in [1], Chap. 7 and 8, is currently used.
These submatrices are inserted into the global sparse FEM matrices ‫ۻ‬ௗ and ۹ ௗ , ordered
1
2

7462

burgi@beuth-hochschule.de, ochmann@beuth-hochschule.deingo5schaefer@bundeswehr.org

according to the above physical quantities, and combined within a frequency-dependent global FEM
matrix ۹௙ ሺ߱ሻ of order ߍ൫۹௙ ൯ ൌ ͸ܰ௡௢ǡ௙௘௠  ൈ ͸ܰ௡௢ǡ௙௘௠ (Eq. (1)).
ଷே೙೚ǡ೑೐೘

ฑ۹௙ ሺ߱ሻ  ൌ ۹ ௗ െ ߱ ‫ۻ‬ௗ ൌ ቎ ۹ఝఝ۹ ௪ఝଷே೙೚ǡ೑೐೘

ฑ
۹
ఝ௪ ቏
۹ ௪௪

ଶ

(1)

This results in the frequency-dependent FEM equation system used for all solution variants:
۹ ఝఝ
൤
۹ ௪ఝ
2.2

۹ ఝ௪ ߮ሬԦͲ൨ ൤ ൨ ൌ ቂ ቃ۹ ௪௪ ܹ
ሬሬሬԦ
െ‫ܨ‬Ԧ

(2)

Used signs and symbols

۳
‫܁‬௘௟

identity matrix
Areas of the surface elements۰࢏ࢊο‫݌‬Ԧ௜ௗߩܿߙ஻ெ‫ݒ‬Ԧ௜௡௖
ܰ௡௢ǡ௙௘௠
ܰ௘௟ǡ௜ௗ۹௙ ሺ߱ሻ
‫ݒ‬Ԧெ௏௉ǡ‫כ‬

BEM matrix (indirect BEM, frequency-dependent, dense)
Pressure difference (indirect BEM)
Density of the surrounding materialSound velocity in the surrounding material
Burton-Miller-factor
incident normal velocity
Number of nodes at the FEM elements
Number of indirect BEM elements
FEM matrix (frequency dependent), also ۹௙
Matrix vector product part (depending on context)

2.3 Structure of the complete system of equations (in combination with the indirect BEM)
Within the BEM calculations, constant plane elements are used (collocation) and the physicalquantities (ȟ‫݌‬௜ௗ , ‫ݒ‬௡ ) are evaluated in the element center. The coupling with the node-based FEM
components is made with two sparse transformation matrices (‫܂‬ଵ and ‫܂‬ଶ).
The method couples the FEM with the indirect BEM (identical fluid on both sides of the elements,
as used within the inertial coupling method, see [2]) and leads to the following system of equations:
߮
ሬԦ
Ͳ
Ͳ
‫ې‬
‫ۍې‬
ሬሬሬԦ ‫ې Ͳ ۍ ۑ‬Ͳ ‫ܹ ێ‬‫ۑ‬‫ێ‬‫ۑ‬‫܁‬௘௟ ‫܂‬ଶ ‫ܨ ێ ۑ‬Ԧ ‫ ۑ‬ൌ ‫ۑ Ͳ ێ‬‫ێ‬‫ۑ‬Ͳ ‫ݒ ێ ۑ‬Ԧ௜ௗ ‫ۑ Ͳ ێ ۑ‬۰௜ௗ ‫ۏ ے‬ο‫݌‬Ԧ ‫ߙۏ ے‬஻ெ ‫ݒ‬Ԧ௜௡௖ ‫ے‬(3)This can be simplified by eliminating ‫ܨ‬Ԧ to۹ ఝఝ ۹ ఝ௪
Ͳ
Ͳ
߮
ሬԦ
Ͳ
‫ې‬
‫ۍ‬
‫ۍې‬
ሬሬሬԦ
۹
۹
Ͳ
െ‫܁‬
‫܂‬
Ͳ
࢝࢝
௘௟ ଶ ‫ۑ ܹ ێ ۑ‬‫ ێ‬௪ఝ
ൌ൦
൪
‫ێ‬
‫ۑ‬
Ͳ
ܑ߱‫܂‬ଵ۳Ͳ ‫ݒ ێ ۑ‬Ԧ௜ௗ ‫ۑ‬‫Ͳ ێ‬ߙ஻ெ ‫ݒ‬Ԧ௜௡௖
‫Ͳ ۏ‬
Ͳ
െ۳ߩܿ
۰௜ௗ ‫ۏ ے‬ο‫݌‬Ԧ௜ௗ ‫ے‬and by eliminating ‫ݒ‬Ԧ௜ௗ to

(4)

۹ ఝఝ
‫ۍ‬
۹ ௪ఝ
‫ێ‬
‫Ͳ ێ‬
‫Ͳ ێ‬
‫Ͳ ۏ‬

۹ ఝ௪
۹ ௪௪
Ͳ
ܑ߱‫܂‬ଵ
Ͳଷൈே೙೚ǡ೑೐೘

‫ ۍ‬ฑ
۹ ఝఝ
‫ێ‬
‫ ێ‬۹ ௪ఝ‫ۏ‬
Ͳor in compressed form:

Ͳ
Ͳ
۳
Ͳ
۳
Ͳ
Ͳ
۳
Ͳ െ۳ߩܿ

ଷൈே೙೚ǡ೑೐೘

ฑ
۹
ఝ௪
۹ ௪௪
ܑ߱ߩܿ‫܂‬ଵ۹௙൤ܑ߱ߩܿ‫܂‬૚௜ௗே೐೗ǡ೔೏

ሬԦ
‫߮ ې‬
Ͳ
ฎ
Ͳ
‫ۑ‬቎ ܹ
ሬሬሬԦ ቏ ൌ ൥ Ͳ ൩െ‫܁‬௘௟ ‫܂‬ଶ ‫ۑ‬ߙ஻ெ ‫ݒ‬Ԧ௜௡௖ο‫݌‬Ԧ௜ௗ۰௜ௗ ‫ے‬െ‫܁‬௘௟ ‫܂‬૛
Ͳ
ሬԦ
‫ݑ‬
൨൤
൨ൌ൤
൨Ǥ
ߙ஻ெ ‫ݒ‬Ԧ௜௡௖ο‫݌‬Ԧ௜ௗ
۰௜ௗ

(5)

(6)

7463

3. USED SOLUTION METHODS
Three different solution methods were implemented to take into account the FEM -specific material
properties:
3.1 Combined sparse matrix (FULL SPARSE mode)
In this mode the total matrix of FEM and BEM coefficients according to Eq. (5) is completely
build and solved using a direct solver for sparse matrices (Intel DSS direct sparse solver or PARDISO
sparse solver).
Advantages:
x For so-called monostatic calculations with many right-hand sides according to the incident soundsource position, only a one-time factorization of the combined matrix is required, and thesubsequent solution process for each right-hand side is very fast.
Disadvantages:
x The memory requirement for the coefficients in a sparse matrix is higher (i.g. twice the size)
since the BEM fraction ۰௜ௗ is fully populated and an additional column index is required per
coefficient entry.
x During the factorization phase, sparse solvers require about 4 to 5 times the memory requirementof the output matrix. Since the solvers generally do not check this in advance, they may crash at
this point.
3.2 Using the Schur method (SCHUR mode)
The basis of the Schur method is the simplified FEM equation system ۹௙ according to Eq. (6).
From this, ‫ݑ‬ሬԦ can be extracted‫ݑ‬ሬԦ ൌ െ۹௙ିଵ ‫܁‬௘௟ ‫܂‬૛ ο‫݌‬Ԧ௜ௗ (7)

and the system of equations in Eq. (6) can be converted intoൣെܑɘɏ ‫܂‬૚ ۹ ି૚ሬԦ୧ୢ ൌ Ƚ୆୑ ሬԦ୧୬ୡǤ୤ ‫܁‬ୣ୪ ‫܂‬૛ ൅ ۰୧ୢ ൧ο

(8)

An FEM coefficient matrix ۸௦௖௛ǡிாெ with۸ୱୡ୦ǡ୊୉୑ ൌ ൣെܑɘɏ ‫܂‬૚ ۹ ି૚୤ ‫܁‬ୣ୪ ‫܂‬૛ ൧Ǥ(9)is build, which is added to the BEM matrix ۰௜ௗ , and onlyሾ۰௜ௗ ൅ ۸௦௖௛ǡிாெ ሿሾο‫݌‬Ԧ௜ௗ ሿ ൌ ߙ஻ெ ‫ݒ‬Ԧ௜௡௖ (10)

has to be solved with a direct equation solver.
However, the formation of the inverse ۹௙ି૚ is very time-consuming and hardly usable in practice.
Alternatively, if an iterative solver is used for which matrix vector products are required per
iteration vector ‫ݔ‬Ԧ௜ according to(11)ܾሬԦ௜ ሺ‫ݔ‬Ԧ௜ ሻ ൌ  ‫ݒ‬Ԧெ௏௉ǡ஻ாெ ሺ‫ݔ‬Ԧ௜ ሻ ൅ ‫ݒ‬Ԧெ௏௉ǡ௦௖௛ ሺ‫ݔ‬Ԧ௜ǡிாெ ሻǡthe FEM matrix vector product part ‫ݒ‬Ԧெ௏௉ǡ௦௖௛ can be calculated much more simply by means of‫ݒ‬Ԧெ௏௉ǡ௦௖௛ ൫‫ݔ‬Ԧ௜ǡிாெ ൯ ൌ െܑ߱ߩܿ‫܂‬ଵ ۹ ିଵԦ௜ǡிாெ Ǥࢌ ‫܁‬ୣ୪ ‫܂‬૛ ‫ݔ‬(12)In order to avoid the formation of ۹௙ିଵ at this point, an auxiliary vector‫ݒ‬Ԧ்ଶ௫ ൌ ‫܁‬௘௟ ‫܂‬૛ ‫ݔ‬Ԧ௜ǡிாெ (13)

is generated and from this a Schur solution vector ‫ݒ‬Ԧ௦௖௛ is calculated by solving the sparse systemof equations۹௙ ‫ݒ‬Ԧ௦௖௛ ൌ ‫ݒ‬Ԧ்ଶ௫ (14)with the aid of a sparse solver. This simplifies Eq. (12) to‫ݒ‬Ԧெ௏௉ǡ௦௖௛ ൫‫ݔ‬Ԧ௜ǡிாெ ൯ ൌ െ݅߱ߩܿ‫܂‬ଵ ‫ݒ‬Ԧ௦௖௛ Ǥ

(15)

Since ۹௙ from Eq. (14) has to be factorized only once and then only has to be solved per iteration
vector ‫ݔ‬Ԧ௜ , the computing time per iteration is much lower.

7464

Advantages:
x Only the memory required for the BEM matrix ۰௜ௗ is needed, the additional ۸௦௖௛ǡிாெ values areadded up. Also, the amount of memory required for factorization usually does not increase for an
equation solver for dense matrices.
x When using iterative solvers, the required matrix vector product parts of Eq. (11) can be preparedby rapid approximate methods, such as the fast multipole method (MLFMM). This will makesignificantly larger numbers of elements possible.
Disadvantages:
x For frequency sweeps, ۹௙ and ‫ݒ‬Ԧ௦௖௛ from Eq. (14) must be reassembled in accordance with Eq.
(1) and factorized for each frequency.

3.3 Using eigenvalues (EIGEN mode)
In this mode, a predefined number ܰ௘௩ of eigenvalues or vectors based on the frequencyindependent FEM matrices ‫ۻ‬ௗ and ۹ ௗ from Eq. (1) is determined by a so-called eigenvalue solver
and stored in an eigenvalue matrix શ௘௩ .
Using these matrices, two "small" square frequency-independent diagonal matrices ۹ ௘௩ and ‫ۻ‬௘௩
are formed by means of்۹ ௘௩ ൌ  શ௘௩ή ۹ ௗ ή શ௘௩ and

(16)்‫ۻ‬௘௩ ൌ  શ௘௩ή ‫ۻ‬ௗ ή શ௘௩ (17)with the order of ࣩሺ۹ ௘௩ ሻ ൌ ࣩሺ‫ۻ‬௘௩ ሻ ൌ ܰ௘௩ .
ሬሬሬԦ )The eigenvalue matrix શ௘௩ is then reduced to its "lower" half (belonging to the displacement ܹ(શ௘௩ǡ௪ ), since the rotation angles are now no longer needed, and two auxiliary matrices ‫ۿ‬ଵ and ‫ۿ‬ଶare generated by‫ۿ‬ଵ ൌ  ‫܂‬ଵ ή શ௘௩ǡ௪ and(18)்‫ۿ‬ଶ ൌ  શ௘௩ǡ௪ή ‫܁‬௘௟ ή ‫܂‬ଶ.

(19)

Now, a frequency-dependent matrix ۲௘௩ is generated by means of۲௘௩ ൌ ሺ۹ ௘௩ െ ߱ଶ ‫ۻ‬௘௩ ሻିଵ,

(20)

whereby the formation of the inverse here is not so time-consuming due to the "small" order ofࣩሺ۲௘௩ ሻ ൌ ܰ௘௩ ൈ ܰ௘௩ .
From this, an auxiliary matrix ‫ۿ‬஽ଶ is generated by means of‫ۿ‬஽ଶ ൌ ܑ߱ሺ۲௘௩ ή ‫ۿ‬ଶ ሻ.

(21)

The matrix-matrix product of these auxiliary matrices is given by۸௘௩ǡிாெ ൌ ߩܿ‫ۿ‬ଵ ‫ۿ‬஽ଶ .

(22)

As in the Schur method, it corresponds to an FEM coefficient matrix of the same order as that of
the BEM matrix and is added to it, so that finally only the system of equationsሾ۰௜ௗ ൅ ۸௘௩ǡிாெ ሿሾο‫݌‬Ԧ௜ௗ ሿ ൌ ߙ஻ெ ‫ݒ‬Ԧ௜௡௖ (23)

must be solved.
Advantages:
x As with the Schur method, only the memory needed for the fully populated BEM matrix ۰௜ௗ isrequired.
x For frequency sweeps, the eigenvalues need to be determined only once (based on the highest
desired frequency), accordingly, the time required decreases after the first for all following
frequencies. Furthermore, only the frequency-independent auxiliary matrices ‫ۿ‬ଵ and ‫ۿ‬ଶ mustbe "saved", all other variables can be discarded to save memory space.
x Using an iterative solver (e.g., GMRES), the time-consuming formation of the ۸௘௩ǡிாெ coefficients can be replaced by a FEM MVP part ‫ݒ‬Ԧெ௏௉ǡ௘௜௚ :ܾሬԦ௜ ሺ‫ݔ‬Ԧ௜ ሻ ൌ  ‫ݒ‬Ԧெ௏௉ǡ஻ாெ ሺ‫ݔ‬Ԧ௜ ሻ ൅ ‫ݒ‬Ԧெ௏௉ǡ௘௜௚ ሺ‫ݔ‬Ԧ௜ǡிாெ ሻ

(24)

7465

This can be calculated very quickly by means of‫ݒ‬Ԧெ௏௉ǡ௘௜௚ ൫‫ݔ‬Ԧ௜ǡிாெ ൯ ൌ ߩܿ‫ۿ‬ଵ ‫ۿ‬஽ଶ ‫ݔ‬Ԧ௜ǡிாெ ,

(25)

since the order of the matrices ‫ۿ‬ଵ and ‫ۿ‬஽ଶ is significantly smaller than that of ۰௜ௗ , while theBEM component‫ݒ‬Ԧெ௏௉ǡ஻ாெ ሺ‫ݔ‬Ԧ௜ ሻ ൌ ۰௜ௗ ‫ݔ‬Ԧ௜ (26)

is conventionally calculated.
x Here too, in the case of iterative solvers, the required matrix vector products of Eq. (26) can beformed by rapid approximate methods.
x It is possible to determine the eigenvalues for individual FEM areas separately, this can lead toa lower number of required eigenvalues and thus reduce the time required.
Disadvantages:
x The number of eigenvalues ܰ௘௩ required for a good quality of the solution can be estimatedpoorly in advance, since this strongly depends on the model and the boundary conditions.
Accordingly, this specification must be checked manually and, if necessary, be increased
adaptively, in which case the time required to determine the eigenvalues also increases
disproportionately.

4. RESULTS FOR A FREQUENCY SWEEP
4.1 Used model
The model used (Figure 1) consists of a conical shell made of 6,058 triangular elements with aradius of 1 m, a thickness of 1 cm and steel as material, surrounded by water.

Figure 1 - used model
A plane incident wave in the negative X direction is used as a sound source to compare the resultswith the FEM COMSOL application, which gives a shorter calculation time due to the rotational
symmetry.

7466

4.2 Results
For a frequency range of 100 Hz to 1 kHz, the normalized sound pressure level in the far field
was determined in 0.5 Hz increments ( ܰ௙ = 1,801 frequencies) at a monostatic evaluation point
at position [10, 0, 0 km] in order to be able to clearly recognize occurring resonances ( Figure 2).
The maximum allowable relative error for the iterative solvers (b , d, g) has been set to 10 -4.

10norm. sound pressure level [dB]

0
-10
-20
-30
-40
-50
-60

0

0,1

0,2

0,3

0,4
0,5
0,6
Frequenz [kHz]

0,7

0,8

0,9

1

a) FEM, full sparse matrix (direct sparse solver, IMKL/DSS)b) FEM, Schur mode (iterative solver, GMRES-HH)c) FEM, eigenvalues,Nev = 500 (direct solver, IMKL/LA+ARP/DSS)d) FEM, eigenvalues, Nev = 500 (iterative solver, GMRES-HH+ARP/UMF)e) rotationally symmetrical 2D-FEM (, COMSOL)f) 3D-FEM (COMSOL)g) only indirect BEM (inertial coupling, iterative solver, GMRES-HH)
Figure 2 - Sound pressure level at a monostatic evaluation point at position [10, 0, 0] km(݂ = 100 … 1,000 Hz, ȟ݂ = 0,5 Hz, ܰ௙ = 1,801)
All three FEM calculation variants (full sparse matrix, eigenvalues and Schur method, a ... d) agree
very well, the differences are, except for some isolated spikes at some frequencies, below one dB.
The difference to the indirect BEM (g) using the thin-shell inertial coupling method according to[2] is clearly seen and shows the great influence of the elastic material properties.
In comparison to the FEM results of COMSOL (e, f), a frequency shift appears with increasing
frequency (with similar maxima and minima values). This is due to the different methods forcalculating the FEM matrices (Kirchhoff method in our inhouse code, Reissner-Mindlin method in
COMSOL) and requires further investigation. For better comparability, therefore, an additional
implementation of the Reissner-Mindlin method for the own code is currently being prepared.
For information, Table 1 gives the calculation times for the different results.
The solution time of the FULL SPARSE variant is not significantly greater here than in the othermethods due to the "small" number of elements.

7467

Table 1 - Calculation times for the frequency sweep
Method / modea) FULL SPARSE (direct)b) SCHUR-Modus (iterative)c) EIGEN-Modus (direct)d) EIGEN-Modus (iterative)e) 2D-FEM (COMSOL, rotationally symmetrical!)f) 3D-FEM (COMSOL)g) only indirect BEM (iterative)

Calculation time [h:m:s]
3:48:55
2:57:57
2:17:00
1:37:16
0:07:12ca. 13 h
0:21:52

4.3 Pressure differences at resonance points
To illustrate the changes in the difference of the pressure on both sides of the elements at theresonance points, the following figures show the real part of the pressure difference Ը൫ȟ‫݌‬௘௟ǡிாெ ൯ on
the surface in the region around the first resonance point between 294 and 299 Hz.

݂ = 294 Hz

݂ = 295 Hz

݂ = 296 Hz

݂ = 298 Hz
݂ = 299 Hz
݂ = 297 Hz
Figure 3 - Ը൫ȟ‫݌‬௘௟ǡிாெ ൯, real part of the pressure difference on the surface
The same value and color range (-15 ... +30 N/m²) is used in all 6 subfigures, so that it is easy tosee how the pressure difference in the area around the resonance point increases significantly at approx.
݂ = 297.15 Hz and then decreases again. Also the sign of the real part changes.

7468

5. SUMMARY AND OUTLOOKThe implemented FEM coupling methods show good agreement with each other and in comparison
with commercial FEM applications.
The speed advantage when using the BEM instead of the FEM in 3D is clearly visible.
The use of iterative equation solvers in conjunction with the implemented SCHUR or EIGEN
solution variants enables the combination with the multi-level fast multipole method.
The code is currently being expanded for FEM shell boundary conditions with respect to closed
structures, i.e., there is a vacuum inside the corresponding structural part.
Also in the future, the combination of arbitrary boundary conditions within the same model should
be possible, as well as the "tight clamping" of existing edges.

REFERENCES
1. Steinke, P., Finite-Elemente-Methode; Rechnergestützte Einführung, 5. Auflage, Springer Vieweg Verlag,
2015, ISBN 978-3-642-53936-7.
2. Burgschweiger, R., Schäfer, I. and Ochmann, M., Implementation and results of a mass inertia couplingas an extension of the BEM for thin shells, Proceedings of the 22nd International Congress on Sound andVibration (ICSV22), Florence, Italy (07/2015).
3. Burgschweiger, R., Schäfer, I. and Ochmann, M., Integration von FEM-Schalenelementen als Randbedingungen in BEM-Berechnungen unter Verwendung verschiedener Lösungsmethoden, Proceedings ofthe 45nd Deutsche Jahrestagung für Akustik (DAGA2019), Rostock, Germany (03/2019).

7469