Mechanical Systems and Signal Processing 148 (2021) 107135
Contents lists available at ScienceDirect
Mechanical Systems and Signal Processing
journal homepage: www.elsevier.com/locate/ymssp
Application of a Krylov subspace method for an efficient
solution of acoustic transfer functions
Christopher Sittl a, Steffen Marburg b, Marcus Wagner a,⇑
a
b
Laboratory Finite-Element-Method, Faculty of Mechanical Engineering, OTH Regensburg, 93053 Regensburg, Germany
Chair of Vibroacoustics of Vehicles and Machines, Department of Mechanical Engineering, Technical University of Munich, 85748 Garching, Germany
a r t i c l e
i n f o
Article history:
Received 27 January 2020
Received in revised form 4 June 2020
Accepted 13 July 2020
2010 MSC:
00-01
99-00
Keywords:
Lanczos algorithm
Krylov-subspace projection
Dirichlet-to-Neumann map
Padé approximation
Fluid–structure interaction
Acoustics
a b s t r a c t
Solving acoustic radiation problems, arising from systems including fluid–structure interaction, is of interest in many engineering applications. Computing frequency response
functions over a large frequency range is a concern in such applications. A method which
solves the Helmholtz equation for multiple frequencies in one step is the matrix-Padé-viaLanczos connection for unsymmetric systems, as presented by Wagner et al. [1]. The present work is based on Ref. [1] and presents a method for efficiently computing frequency
responses over a frequency range for coupled structural-acoustic problems, where the
structure and the acoustic near field are discretized with finite elements and an analytical
Dirichlet-to-Neumann map approximates the far field. The method is based on a Krylovsubspace projection technique which derives a matrix-valued Padé approximation for a
restricted area in the near field and the pressure field on a spherical boundary. On the
spherical boundary, where the finite domain is truncated, the non-local modified
Dirichlet-to-Neumann operator is applied as a low-rank update matrix. The present contribution extends this method and incorporates new techniques for a more stable model
reduction through the Lanczos algorithm and a novel weighted adaptive windowing technique. Further, structural damping is incorporated, for computing the acoustic radiation of
a harmonically excited plate. These computed results are compared with acoustic measurements in an anechoic chamber and verified with computational results obtained with
a commercial code that uses the perfectly matched layer method.
Ó 2020 Elsevier Ltd. All rights reserved.
1. Introduction
The dynamic response of linear systems is of interest in many engineering applications, such as acoustic design optimization for interior and exterior problems or transfer path analysis. A transfer function provides insight in the dynamic behavior
of a dynamical system, but generally requires the evaluation of the dynamic system equations over a frequency range.
Depending on the size of the desired frequency range and the system size, the computational effort for solving these systems
can be prohibitively expensive, especially if the size of the frequency range exceeds a particular limit [2]. To circumvent this
problem, an approach which provides solutions for multiple frequencies, can be taken into consideration. Such a method for
solving linear systems arising from a finite element discretization for multiple frequencies is presented.
⇑ Corresponding author.
E-mail address: marcus.wagner@oth-regensburg.de (M. Wagner).
https://doi.org/10.1016/j.ymssp.2020.107135
0888-3270/Ó 2020 Elsevier Ltd. All rights reserved.
2
C. Sittl et al. / Mechanical Systems and Signal Processing 148 (2021) 107135
A number of approaches are known for efficient frequency sweeps in structural acoustic radiation problems. Modal analysis and modal superposition, which are very popular for pure structural analysis and interior acoustic problems, have been
applied rather seldom probably due to the difficulty of formulating the acoustic eigenvalue problem for exterior problems. A
solution applying infinite elements for discretization was provided in [3]. It could be shown, that only a few modes are relevant for a model-order reduction, but it has been impossible to identify these modes a priori. This approach has been revisited in recent years [4], where further ideas have been given to identify modes which are relevant for modal superposition.
The technique had been applied to an optimization problem in which the radiated sound power was efficiently evaluated by
superposition of all modes and over a large frequency range [5]. However, practicability of the modal approach using infinite
elements remains limited so far. Peters et al. [6] proposed a different modal approach, which formulates a fully coupled problem using finite elements for the structure and boundary elements for the fluid. A superposition of the modes of the coupled
system has shown substantial speed-ups compared to the solution of the entire system solved at many individual frequencies successively. Wixom and McDaniels [7] have presented another technique for efficient frequency sweeps even with a
large number of load cases. Very recently, Baydoun et al. [8] developed a greedy algorithm using a low rank approximation
for an efficient frequency sweep.
A method which is well studied for solutions over a frequency window is the Padé-via-Lanczos connection (PVL). This
connection between the Lanczos process and the Padé approximation, has been observed by Gragg [9] and has also been
extended to a block case by Freund [10] to the so-called matrix-valued PVL (MPVL). The MPVL method is well known and
has been applied to several acoustic and fluid–structure interaction problems [1,11–14]. There is also a connection to the
block Arnoldi method for multi-input multi-output approximations [15,16].
In several applications only partial field solutions of the computational domain are of interest. These partial fields include
single points on a vibrating structure or in the near field of a vibrating structure. Furthermore, points on an enclosing surface
such as a sphere for far-field computations can be included in these partial fields. Hence, a method like the MPVL method is
suitable in these applications, since it provides such a partial field solution, instead of the solution of the complete domain.
This contribution extends the formulations from Wagner et al. [1] and Tuck-Lee et al. [12] for fluid–structure interaction
computations to three dimensions and introduces an adaptive windowing algorithm with weighting. It is based on an error
estimation relative to previously constructed Krylov subspaces. Van de Walle proposed a similar method in [16], for detecting the window width and stagnation of the algorithm. He uses this for model order reduction with an Arnoldi-type algorithm and also uses a time integration method with the reduced models for solutions in the time domain [16]. In the
present paper, a new scheme with overlapping frequency windows is introduced, in order to increase the robustness of
the used method.
In terms of far-field simulations, the Sommerfeld radiation condition needs to be incorporated on the boundary of the
finite elements (FE) model. There are various types of boundary conditions like the Perfectly-Matched-Layer (PML) technique
[17] or the non-local Dirichlet-to-Neumann (DtN) boundary condition [18]. Also, infinite elements can be considered for a FE
model as a non-reflecting boundary formulation. These methods have been used also in conjunction with model reduction
techniques [1,14–16]. A drawback of the PML method, in the context of the MPVL, is that the system size grows, due to the
elements needed in the PML region. The approach with infinite elements would also lead to an increased system size (see
[1,14]). For the non-local DtN map, the system size remains the same although, due to its non-local character, the mapping
matrix becomes fully populated, similar to coupling the finite element method (FEM) with the boundary element method
(BEM) [19]. To take advantage of the preservation of the system size, the DtN map is chosen in this work. With the FEMBEM-coupling, the vibrating structural FEM part is coupled with the BEM for far field computations. The PVL technique is
also applicable to a BEM approach, which is shown in Ref. [20]. However, the coefficients of the Padé approximation are computed directly, which makes the computation of several derivatives of the dense system matrices necessary.
The second new contribution in this work is the introduction of damping effects in the MPVL approach, to improve the
estimated window width of the reduced system. An option for incorporating damping in the Lanczos algorithm is the reformulation of the second-order differential equation in two first-order differential equations [21]. Another possibility is the
second order Arnoldi reduction, which was presented by Bai and Su [22]. Incorporating Rayleigh damping in the Lanczos
algorithm for estimating a Padé approximation has been studied by Meerbergen [23]. The present work investigates the
effect of damping on the projection onto the Krylov subspace and the resulting Padé approximation, by incorporating structural hysteretic damping, as this approach does not necessitate the reformulation of the resulting matrix system, see also
[24].
In the following, the matrix representation of the Helmholtz equation for exterior acoustics and the coupling of the fluid
and structural domain, are formulated by incorporating a DtN map on the spherical boundary, suitable for an FE discretization. Afterwards, the modified version of the DtN map is introduced, which is then reformulated as a transfer function representation, suitable for the MPVL algorithm.
After this, the new weighted adaptive windowing technique with overlapping window ranges is introduced. Next, the
hysteretic damping model for the structural domain is established. Finally, acoustic measurements of the radiation of a plate,
carried out in the frame of this work, are presented. These measurements validate the presented approach. Moreover, computational results obtained with the PML method demonstrate its robustness and efficiency.
C. Sittl et al. / Mechanical Systems and Signal Processing 148 (2021) 107135
3
2. Dirichlet-to-Neumann map for exterior acoustics
Considering the Helmholtz problem, the infinite domain has to be truncated in the FE approach by an artificial boundary.
In this section, the problem for the finite domain is stated and the DtN map is introduced as a non-reflecting boundary condition for the truncated finite domain. Further, the modified version for the DtN radiation condition is introduced to ensure
uniqueness and regularity of the solution. For further details on the modified version of the DtN map, the reader is referred to
Refs. [25,26]. Spherical coordinates are defined as x ¼ r sin h cos / ; y ¼ r sin h sin / and z ¼ r cos h. In Fig. 1, the domain for
exterior structural acoustic computations is illustrated. The infinite domain X1 is truncated by the surface CDtN . The surface
encloses the radiation sources completely. The domain between the elastic body Xs and the truncation surface CDtN is
denoted as Xf and represents the fluid domain. The interface for fluid–structure interaction is denoted as CI .
2.1. Problem statement
Considering harmonic time dependence eixt for the compressible fluid in an infinite domain X1 , the boundary value problem for the fluid domain is stated as:
r2 p þ j2 p ¼ 0 in Xf ;
p ¼ g f on Cg;f and
rp n ¼ hf on Ch;f ;
ð1Þ
ð2Þ
ð3Þ
The above equations find the solution of the acoustic pressure p, satisfying the Helmholtz equation with Dirichlet and Neumann boundary conditions, denoted as g and h, respectively, and j ¼ x=c is the wave number for the fluid domain, with x
denoting the angular frequency, and c the wave propagation velocity. Further, n is defining the outward normal of the coupling interface. In addition, the Sommerfeld radiation condition is enforced by
lim r
r!1
@p
@r
ijp
¼ 0;
ð4Þ
which ensures uniqueness of the solution by excluding reflections from infinity. In the acoustic domain, an elastic structural
body with the volume Xs and the boundary CI is submerged. The structural part is described by standard linear elasticity
equations in terms of displacements u
De u þ qs x2 u ¼ f in Xs ;
u ¼ g s on Cg;s ;
r n ¼ hs on Ch;s ;
ð5Þ
ð6Þ
ð7Þ
where qs is the structural density and De is the elasticity operator, cf. [19]. The Cauchy stress tensor is denoted as r ¼ C :
and is related to the strain via the elasticity tensor C. The coupling interface between the structure and the fluid is constituted by CI ¼ Cf and is governed by the two following equations:
Fig. 1. Computational domain for exterior structural acoustics.
4
C. Sittl et al. / Mechanical Systems and Signal Processing 148 (2021) 107135
r n on CI and
ð8Þ
rp n ¼ qf x2 n u on CI :
pn¼
ð9Þ
This results in a coupling between the fluid and the structural domain. The above Eqs. (1)–(7), (9) state the mathematical
problem for finding p and u in the case of exterior vibro-acoustic computations in an unbounded domain.
2.2. Dirichlet-to-Neumann operator
For the solution of the exterior problem with finite elements, a spherical boundary CDtN is introduced, which truncates the
unbounded fluid domain. By imposing a so-called Dirichlet-to-Neumann (DtN) map on the spherical boundary, as proposed
by Keller and Givoli [26], the mathematical statements from Eqs. (1)–(4) can be stated as an equivalent bounded problem by
adding
rp n ¼ BDtN ðpÞ on CDtN ;
ð10Þ
to Eqs. (1)–(3), where the operator BDtN is the so-called DtN map. The DtN operator maps the solution, which is radiated by a
sphere onto the surface of this sphere. The solution in the far-field for the sound pressure in X1 , thus the outside of a sphere,
is
rffiffiffi 1 n 0 ð1Þ
Z
RXX hn ðkrÞ ð2n þ 1Þðn jÞ!
pðr; h; /Þ ¼
Pjn ðcos hÞP jn ðcos h0 Þ cos jð/
r n¼0 j¼0 hðn1Þ ðkRÞ 4pR2 ðn þ jÞ!
CDtN
/0 ÞpðR; h0 ; /0 ÞdC0DtN ;
ð11Þ
based on the formulations of Keller and Givoli [26], where pðR; h0 ; /0 Þ is the solution on the spherical surface with radius R and
pðr; h; /Þ the solution in X1 , hence r > R. The DtN boundary condition is an exact non-reflecting boundary condition for time
harmonic problems. Based on the definition from Keller and Givoli [26], the operator BDtN ð pÞ in three dimensions for the DtN
map is formulated as:
BDtN ð pÞ ¼
1 X
n0
X
ð2n þ 1Þðn
n¼0 j¼0
ð1Þ0
jÞ! jhn ðjRÞ
4pR2 ðn þ jÞ!
ð1 Þ
hn ðjRÞ
Z
CDtN
Pjn ðcos hÞPjn ðcos h0 Þ cos jð/
/0 ÞpðR; h0 ; /0 ÞdC0DtN ;
ð12Þ
where dC0DtN ¼ R2 sin h0 dh0 d/0 . Here P jn ðcos hÞ are the associated Legendre polynomials of degree n and order j. Furthermore,
ð1Þ0
ð1Þ
hn ðjRÞ are the spherical Hankel functions of the first kind and hn ðjRÞ the derivatives with respect to the argument. The
prime on the sum in (12) indicates that for j > 0 a factor of two is pre-multiplied. p contains the solutions on the spherical
surface of radius R. Note, that the DtN-operator maps the Dirichlet datum outside of the sphere to a Neumann datum on the
sphere.
2.3. Modified Dirichlet-to-Neumann operator
The infinite series in (12) is truncated in practical applications and in the following, it is truncated at the number of series
terms N DtN . Due to the truncation of the series, the boundary condition is no longer exact. Furthermore, free vibration modes
are introduced for n > N DtN , because the boundary condition then reduces to rp n ¼ 0 on the surface. Grote and Keller present in Ref. [25] a so-called modified version of the DtN-operator which utilizes a local boundary condition in spherical coordinates on the boundary. The local boundary condition in the present approach is, as presented by Grote and Keller, the first
order Bayliss-Gunzburger-Turkel boundary condition with its operator
B1 ðpÞ ¼ ij
1
p:
R
ð13Þ
By modifying the DtN-operator with the local boundary condition B1 , the boundary condition for n > N DtN is now rp n ¼ B1 ,
which eliminates the not physically motivated real eigenvalues for n > N DtN . For n 6 N DtN the DtN-operator remains exact,
because the local boundary condition cancels out:
Bmod ¼ ðBDtN B1 ÞðpÞ þ B1 ðpÞ:
|fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflffl}
ð14Þ
BDtN;mod
The modified version of the DtN-operator BDtN;mod is then derived by inserting (13) and (12) in (14):
BDtN;mod ð pÞ ¼
N
n 0
DtN X
X
n¼0 j¼0
zn;j ðjRÞ
Z
CDtN
Pjn ðcos hÞPjn ðcos h0 Þ cos jð/
/0 ÞpðR; h0 ; /0 ÞdC0DtN :
ð15Þ
Here, zn;j ðjRÞ is introduced, which can be interpreted as modified impedance coefficients on the spherical surface, defined
by:
C. Sittl et al. / Mechanical Systems and Signal Processing 148 (2021) 107135
zn;j ðjRÞ ¼
ð2n þ 1Þðn
jÞ!
4pR2 ðn þ jÞ!
!
0
jhðn1Þ ðjRÞ
1
ij þ
:
ð1 Þ
R
hn ðjRÞ
5
ð16Þ
2.4. Finite element discretization
For the derivation of a transfer function, the above stated problem is discretized with a standard Galerkin scheme. Introducing N as a vector of finite element shape functions for the exterior acoustic problem and N s as the matrix of shape functions for the structural domain, one arrives at the matrix form:
1
"
#
"
#
B
C
B Ks
Ms
0
0
0
0 0 C u
fs
RT
C
B
x2
þ
þ
C p ¼ f ;
B 0
q
R
M
0
K
0
B
f
K
DtN
1
f
f
A
@
f
|fflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflffl}
|fflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflffl} |fflfflfflfflfflffl{zfflfflfflfflfflffl} |ffl{zffl} |fflffl{zfflffl}
0
^
K
^
M
^
K
DtN
^
B
1
x
ð17Þ
f
^ 2 RNN contains the stiffness matrix for the structural domain K s 2 RNs Ns and the acoustic domain K f 2 RNf Nf , as
where K
well as the coupling matrix R 2 RNf Ns , defined in its integral form as:
R¼
Z
CI
N f nT N Ts dC:
ð18Þ
Further, x represents the discrete displacements u of the structural domain and sound pressure values p for the fluid domain.
f contains loads on the structural domain and the fluid domain. Note, that the system size is not increasing, due to the presence of the DtN boundary condition. Thus, K DtN represents the boundary operator for the DtN map in the fluid domain Cf and
can be sorted in a block form
K DtN ¼
0
0
0 K CDtN
ð19Þ
;
if the degrees of freedom on the boundary are sorted last. K CDtN is a square matrix of size N b N b , with N b , denoting the
degrees of freedom on CDtN . Its integral form is
K CDtN ¼
Z
NBDtN;mod N T dCDtN :
ð20Þ
CDtN
^ and M
^ are unsymmetric, but remain, under proper ordering, sparse and banded matrices, whereas this is not the
Note that K
^
^ and M
^ are frequency independent, which is not
case for K DtN , due to the non-local character of the DtN map. Furthermore, K
^ DtN and B1 , whereas B1 is the matrix representation of the local boundary condition and can be
the case for the matrices K
written as:
B1 ¼
Z
CDtN
ij
1
NN T dCDtN :
R
ð21Þ
To derive a transfer function formulation for N E near field points, the restriction operator E is introduced as
E ¼ e1 . . . eNE 2 RNNE ;
ð22Þ
where ej ; j ¼ 1; . . . ; N E are basis vectors for each near-field point of interest. This leads to a transfer function for the coupled
exterior structural domain for a restricted area, inside the bounded domain:
h
^
H ðxi Þ ¼ ET x ¼ ET K
^ þK
^ DtN þ B
^1
x2i M
i
1
f 2 CNE 1 ;
ð23Þ
where xi defines the angular frequency for every frequency point i of interest. Due to the restriction operator, H now contains pressure values or displacement values for restricted points inside the computational domain.
2.5. Reformulation of the Dirichlet-to-Neumann operator
Recall, that K DtN is a square dense matrix of size N f N f . Further, the DtN-operator is applied on the discretized boundary
CDtN . In Ref. [27] it is shown that rankðK DtN Þ ¼ Nmod ; Nmod Nb . In the following, it is introduced how to reformulate K DtN , to
^ DtN as:
arrive at a formulation for K
^ DtN ¼ V
^V
^T;
K
^ as a N N mod matrix. Starting with inserting the modified DtN operator (14) in (20), we obtain
with V
ð24Þ
6
C. Sittl et al. / Mechanical Systems and Signal Processing 148 (2021) 107135
K CDtN ¼
Z
N
CDtN
"
N
n 0
DtN X
X
znj
n¼0 j¼0
Z
CDtN
Pjn ðcos hÞPjn ðcos h0 Þ cos jð/
0
T
0
DtN
/ ÞN dC
#
dCDtN :
ð25Þ
By further applying a trigonometric summation rule and interchanging integration and summation, the matrix representation of the DtN boundary condition is derived as
K CDtN ¼
N
n 0
DtN X
X
znj c n;j c Tn;j þ sn;j sTn;j ;
ð26Þ
n¼0 j¼0
where the vectors cn;j 2 RNb 1 and sn;j 2 RNb 1 are defined by
c n;j ¼ R2
Z p Z 2p
N ðh; /ÞPjn ðcos hÞ cos j/ sin h d/ dh
ð27Þ
Z p Z 2p
N ðh; /ÞP jn ðcos hÞ sin j/ sin h d/ dh:
ð28Þ
0
0
and
sn;j ¼ R2
0
0
These vectors contain the discrete surface harmonics for the discrete points on the spherical boundary CDtN and are zero
otherwise. Therefore, each vector contains N C nonzero elements, where N C is the number of points on CDtN , due to discretization. Furthermore, the matrices F and K are introduced:
F ¼ c 0;0 ; c 1;0 ; c 1;1 ; . . . ; cNDtN ;NDtN ; s1;1 ; s2;1 ; s2;2 ; . . . ; sNDtN ;NDtN 2 RNb Nmod
K ¼ diag z0;0 ; z1;0 ; z1;1 ; . . . ; zNDtN ;NDtN ; z1;1 ; z2;1 ; z2;2 ; . . . ; zNDtN ;NDtN 2 CNmod Nmod ;
ð29Þ
ð30Þ
sorting the coefficients from (16), (27) and (28) in matrices to arrive at a matrix representation of K CDtN . The rank of F is
N mod ¼ ðN DtN þ 1Þ2 and K is a diagonal matrix. With F and K it follows for (26)
K CDtN ¼ FKF T ¼ VV T with V ¼ FK1=2 :
ð31Þ
Note, that the columns of F build an orthogonal set, due to the spherical harmonics and they are therefore linearly independent [28]. For the coupled problem the matrix can be reformulated to derive at Eq. (24) with
^ ¼ 0ðNs þNf
V
N b ÞNmod
V
2 CðNs þNf ÞNmod ;
ð32Þ
with the same rank N mod as V.
2.6. Numerical integration
In (27) and (28), the shape functions N ðh; /Þ are defined in global spherical coordinates, as well as the Legendre polynomials in conjunction with sin j/ and cos j/, respectively. In the present approach, the shape functions are defined and integrated in their local coordinate system with a Gauss–Legendre quadrature formula
c n;j ¼
nw X
nw
X
N ðni ; gk Þwi wk detðJ ÞPjn ðcos hi Þ cos j/k sin hi ;
ð33Þ
i¼1 k¼1
where, hi and /i are the corresponding global coordinates with respect to the integration points ni and gi . Note, that the
radius R cancels out with the orthonormalization coefficient. Furthermore, the associated Legendre functions have to be calculated for each degree n, order j and location h. The integrals are evaluated via Gaussian quadrature with the integration
points ni and gk and their corresponding weights wi and wk (see Eq. (33)). The shape functions are integrated in their local
coordinate system, utilizing the isoparametric concept with the Jacobian J. Due to the products in Eq. (33) the number of
integration points has to be sufficient. The Gaussian quadrature formula is exact for polynomials up to a degree of
p ¼ 2nw 1, where nw is the number of points. From Eq. (33), it is clear, that the number of points nw has to be chosen higher,
than the polynomial degree of the shape functions. For a validation which nw is sufficient for the integration of cn;j and sn;j
respectively, the error for different nw are computed for c n;j and sn;j compared to the analytic solution:
¼
kc n;j;gauss c n;j;analytic k2
:
kc n;j;analytic k2
ð34Þ
The results are computed up to a degree n ¼ 10 and order j ¼ 10 of the associated Legendre polynomials, which are the highest degree and order in the present work. The sufficient choice for nw depends on the location and the size of the global element. Numerical tests have shown, that for nw P 9 the relative error is well below 1 10 12, for a linear quadrilateral
element, whereas the largest error appears near the poles of the sphere at h 0 and h p. In Fig. 2 the relative error is shown
C. Sittl et al. / Mechanical Systems and Signal Processing 148 (2021) 107135
Fig. 2. relative error
7
of the numeric solution in comparison to the analytic solution, over the number of integration points nw .
over the number of points, for a linear quadrilateral element with edge lengths of Dh ¼ D/ ¼ 0:05 and the midpoint at
h ¼ / ¼ 0:075.
3. Spherical harmonic expansion
Since, our motivation is not only to solve for a restricted area in the finite computational domain X, but also for points
outside the spherical boundary CDtN , we follow the approach from Refs. [1,11], where the pressure field on the enclosing surface CDtN is expressed in a series representation. Recall, that the pressure field outside the spherical boundary can be
expressed with (11), if the pressure-field on the boundary is known. Instead of solving for the discrete pressure-field, the
solution is obtained by finding the modal coefficients of the spherical harmonic expansion:
pðR; h; /Þ ¼
N DtN X
n
X
an;j Pjn ðcos hÞ cos j/ þ bn;j Pjn ðcos hÞ sin j/;
ð35Þ
n¼0 j¼0
with the modal coefficients an;j
an;j ¼
ð2n þ 1Þðn
jÞ!
2
4pR ðn þ jÞ!
Z p Z 2p
N ðh; /ÞPjn ðcos hÞ cos j/ sin h d/ dh p
ð36Þ
Z p Z 2p
N ðh; /ÞPjn ðcos hÞ sin j/ sin h d/ dh p;
ð37Þ
0
0
and bn;j :
bn;j ¼
ð2n þ 1Þðn
jÞ!
4pR2 ðn þ jÞ!
0
0
with p ¼ Np. By choosing the number of terms in the spherical harmonic expansion equal to the order of the DtN series
expression, the transfer function for points on the enclosing surface can be rewritten, by recalling Eqs. (27)–(29):
H ðxÞ ¼ a0;0 ; a1;0 ; a1;1 ; . . . ; aNDtN ;NDtN ; b1;1 ; b2;1 ; b2;2 ; . . . ; bNDtN ;NDtN
T
¼ DF T x:
ð38Þ
Nmod Nmod
Here, D 2 R
is a diagonal matrix, which represents the orthonormalization coefficients of the spherical harmonic
coefficients. The transfer function from Eq. (23) then yields:
H ðxi Þ ¼
"
#
DF T h ^
K
ET
^V
^T
^ þB
^1 þ V
x2i M
i
1
f:
ð39Þ
Now, we have arrived at a transfer function formulation for exterior structural acoustic problem, which solves for restricted
points of interest in the computational domain, due to the restriction operator E and solutions on the spherical boundary.
With the above transfer function, the modal coefficients on the spherical boundary are computed, instead of discrete solutions on the boundary. Therefore, the sound pressure on the boundary is computed in a post-processing step by Eq. (35). For
solving the matrix system for multiple frequencies xi simultaneously, the frequency dependence is simplified in the
following.
4. Multi-frequency analysis
For the application of the multi-frequency analysis, the complex frequency dependency of the transfer function in Eq. (39)
is simplified. For further investigations a frequency-shift parameter ri is introduced
ri ¼ x2i
x20 ;
ð40Þ
8
C. Sittl et al. / Mechanical Systems and Signal Processing 148 (2021) 107135
which shifts the matrix system by the square of the angular frequency x0 , which is in the following referred to as the reference frequency. Now, the transfer function from Eq. (39) is stated in a shifted form as
H ðri Þ ¼
"
#
DF T h
E
T
^ T ðri Þ
^ þV
^ ðri ÞV
ri M
A0
i
1
f;
ð41Þ
where A0 contains frequency independent parts
^
A0 ¼ K
^ þB
^ 1 ðx0 Þ:
x20 M
ð42Þ
Note, the modification with the local boundary condition in B1 is set on the reference frequency x0 . Recall, B1 ensures that A0
is non-singular for x0 > 0, as it was introduced in Section 2.4. This also holds by setting B1 constant at the reference fre^V
^ T for the matrix
quency x0 (see Grote [25]). What remains is the complex frequency dependence of the low rank update V
representation of the DtN-operator. Due to the low rank update structure, the Sherman–Morrison–Woodbury formula is
applicable for the inverse. With the application of the Sherman–Morrison–Woodbury formula a 2 2 block partitioned
matrix W ðri Þ can be obtained:
W ð ri Þ ¼
wF
WF
wE
WE
¼
"
#
FT h
E
T
I
ri A0 1 M
i
1
A0 1 ½ f
F ¼ LT ½I
ri A 1 R:
ð43Þ
Here, I is a square identity matrix of order N. The block partitioned matrix is reformulated with L ¼ ½F E; R ¼ A0 1 ½f F and
A ¼ A0 1 M. For a detailed description of the Sherman–Morrison–Woodbury formula application for the proposed method
and the derivation of the above relation, the reader is referred to Ref. [1]. To derive a transfer function representation
H ðri Þ, two separate steps are performed. First, solving the small dense system of size N mod N mod
K
1
þ W F q ¼ wF
ð44Þ
for q and second, inserting q in
H ðri Þ ¼
DðwF
wE
W F qÞ
WEq
2 CNmod þNE 1 :
ð45Þ
Recall, A0 is a complex non-symmetric and M is a non-symmetric matrix, therefore A 2 CNN is complex non-symmetric.
L 2 CNp and R 2 CNm are called left and right starting blocks in the following. Hence, Eq. (43) is a matrix-valued transfer
function of a time-invariant linear dynamical system with m ¼ N f þ N mod inputs and p ¼ N E þ N mod outputs, including Eqs.
(44) and (45). In practical applications A can be very large, therefore the factorization for each ri of interest can be prohibitively expensive, if the frequency range and the number of frequencies are large. Furthermore, the solution of systems
arising from fluid–structure-interaction are often ill-conditioned, which makes the application of iterative methods difficult.
Also, finding a sufficient preconditioner for A can be expensive, due to the inversion of A0 . A method which provides a solution for multiple frequencies is the matrix-valued Padé-via-Lanczos (MPVL) algorithm, thus the solution of Eq. (43) for multiple ri simultaneously, see [29].
4.1. Matrix-valued Padé-via-Lanczos connection
Considering W ðrÞ as a matrix valued rational function, its approximation can be obtained by exploiting its series representation. Since the rational function has poles, which are related to the eigenvalues of A, a matrix valued Padé approximation is an auspicious approach. In the following, the connection between the Lanczos algorithm and Padé approximation is
utilized. This connection is well known and has been observed by Gragg [9] and has also been extended to a block version by
Freund [10]. The motivation for the MPVL algorithm is to obtain an approximation of W which is referred to as the nth
matrix-Padé approximant W n . By further utilizing the connection to the Lanczos process, a representation for W n can be
derived
W n ðri Þ ¼ gTn Dn 1
ri T n Dn 1
1
qn ;
ð46Þ
which is expressed by the matrices gn ; Dn ; T n and qn , which are obtained after n iterations of the Lanczos process. We obtain
W n as a reduced-order transfer function (since n N), with respect to the reference frequency x0 . The order of approximation is then denoted as
qðnÞ 6
jnk
m
þ
n
:
p
ð47Þ
Therefore, in the scalar case, setting m ¼ p ¼ 1, the maximal order of approximation qðnÞ ¼ 2n. For a proof of the connection
between the matrix Padé approximant and the Lanczos process see [10]. Since, the projection of A onto T n results in a
reduced dimension system of size n n, the Padé approximant belongs to the group of reduced-order modeling approaches.
C. Sittl et al. / Mechanical Systems and Signal Processing 148 (2021) 107135
9
4.2. Block Lanczos algorithm
In this section, a Krylov subspace method is introduced, to obtain the coefficient matrices, which are introduced for the
reduced-order transfer function W n . This exploits the well known connection between the Lanczos process and Padé approximation. For the following discussion, we recall that
R ¼ ½r 1 r 2 ; . . . ; rm 2 CNm
and L ¼ l1 l2 ; . . . ; r p 2 CNp
ð48Þ
are block matrices. In the following the m vectors in R are referred to right starting vectors and R as the right starting block
and the p vectors in L as left starting vectors and the left starting block, respectively. Aliaga et al. [30] developed a Lanczostype algorithm for multiple starting vectors, which also resolves the non-symmetric structure of A. Furthermore, this algorithm performs so-called deflation steps, that deflate linear dependent or nearly linear dependent vectors from the Lanczos
blocks. Due to finite arithmetic, the Lanczos-type process can breakdown, due to divisions close to zero. To circumvent this
problem, proper look-ahead steps are performed by the algorithm. In the following, the properties of the introduced matrices
in the reduced-order transfer function are introduced in absence of deflation and look-ahead steps, as their relation to the
Lanczos-type process.
Starting with the introduction of the block Krylov matrices
h
i
K ðA; RÞ ¼ R; AR; A2 R; . . . ; Aj 1 R and K AT ; L ¼ L; AT L; AT
2
L; . . . ; AT
k 1
L :
ð49Þ
The Lanczos-type algorithm generates two sequences of basis vectors via three-term recurrences. These vectors are referred
to left and right Lanczos vectors and span the nth block Krylov subspaces Kn after n iterations in a vector wise manner, such
that:
spanfw1 ; w2 ; . . . ; wn g ¼ Kn AT ; L ;
ð50Þ
spanfv 1 ; v 2 ; . . . ; v n g ¼ Kn ðA; RÞ:
ð51Þ
These Lanczos vectors satisfy the biorthogonality condition
W Tn V n ¼ ½w1 ; . . . ; wn T ½v 1 ; . . . ; v n ¼ Dn
ð52Þ
where the matrices W n and V n are introduced and Dn denotes a n n matrix. Furthermore, the recurrences after n iterations
of the Lanczos-type algorithm, are summed up as:
V nT n
W n T~ n
m
p
¼
¼
½r 1 r2 . . . r n if 1 6 n 6 m;
(
AV n
m
if n > m;
½l1 l2 . . . ln if 1 6 n 6 p;
AT W n
m
if n > p:
ð53Þ
ð54Þ
~ n p are rectangular matrices, containing the recurrence coefficients of the biorthogonalizaIn the above equations T n m and T
tion process. In the following, just the square parts of both matrices are used, as it is recommended by Aliaga et al. for MPVL
[30]. Furthermore, gn 2 Cnm and qn 2 Cnp are introduced as upper triangular matrices, which are related to the generated
Lanczos vectors, the starting blocks and the diagonal matrix Dn with the following two equations:
V Tn L ¼ DTn gn and W Tn R ¼ Dn qn :
ð55Þ
For further information about this Lanczos-type algorithm, the reader is referred to [30].
Note, the above equations are introduced in absence of look-ahead and deflation steps. By incorporating look-ahead steps,
the biorthogonalization of the generated Lanczos vectors, via the modified Gram-Schmidt process is relaxed to a cluster-wise
biorthogonalization of the previous generated Lanczos vectors. If look-ahead steps are performed, Dn is no longer a diagonal
matrix, but a block-diagonal matrix. Furthermore, the bandwidth of T n is increased, if look-ahead steps occur. In absence of
look-ahead steps the bandwidth is m þ p þ 1.
In each iteration of the Lanczos-type algorithm two matrix vector products with A are performed. Since A ¼ A0 1 M, a LU
factorization of A0 is performed before the Lanczos-type algorithm and the two matrix vector products with two backward
and forward substitutions
Avn ¼ A0 1 Mvn ¼ U 0 1 L0 1 Mvn
ð56Þ
AT wn ¼ M T A0 T wn ¼ M T L0 T U 0 T wn
ð57Þ
and
10
C. Sittl et al. / Mechanical Systems and Signal Processing 148 (2021) 107135
are performed, which advance the left and right Krylov subspace. Therefore, the main costs are one LU factorization before
the Lanczos-type algorithm and then two simple forward and backward substitutions in each iteration, whereby the matrix
A0 1 M is never computed explicitly. Liew and Pinsky presented a complexity estimation of the MPVL algorithm for single
starting vectors, by means of flop count, in comparison to conventional LU factorization [13]. They stated an improvement
when i > n=k, where i is the number of frequencies, n is the number of Lanczos iterations and k the bandwidth of A0 .
5. Weighted adaptive windowing
As explained before, the rational function of the reduced-dimension system is expanded via the MPVL algorithm at a reference frequency x0 . Therefore, the range around the reference frequency, in which the algorithm converges, is expected to
increase in conjunction with increasing order of the reduced-dimension system. Thus, by further advancing the Krylov subspace, via the Lanczos-type algorithm, a wider range of convergence should be expected. Since the range of convergence,
which is denoted as frequency window in the following, is not known beforehand, an appropriate method which detects
the edges of this window is necessary.
Tuck-Lee and Pinsky presented an adaptive windowing algorithm in Ref. [12], which incorporates an error estimation,
based on the approximation convergence during the generation of the Krylov subspace. This algorithm estimates the relative
error between an actual reduced-order model (ROM), constructed at iteration n of the Lanczos-type algorithm, and a previous constructed ROM of size q < n. This method is applied to the same MPVL expansion which is used in the present work.
For rational Krylov subspace techniques, hence for multiple expansion points, other error estimates have been published
which are based on a relative error of different constructed Krylov subspaces or on a construction of a residual [16]. These
methods are not applicable in a straight forward manner to our present approach, which generates the Krylov subspace just
for a single expansion point at a time and are not further considered here.
Not only the estimation of the frequency window, where the constructed ROM results in an appropriate approximation of
the original system is intended of the adaptive windowing algorithm, but also detecting stagnation during the Lanczos-type
process. Stagnation occurs, if the frequency window is not increasing, by further advancing the Krylov subspace [12,16].
We use a similar error estimation as presented in Ref. [12], which is based on the relative error between two subsequent
ROM
en ðri Þ ¼
H n ðri Þ
H q ðri Þ2
;
jjH n ðri Þjj2
ð58Þ
where H n is the representation of the transfer function for the actual Lanczos iteration n and H q the transfer function of the
ROM at a previous iteration q < n. If this relative error exceeds a user-defined limit, the algorithm detects one side of the
window. For this, the small dense system from Eq. (46) has to be solved at the end of certain iterations in the Lanczos process, for each ri in the frequency window. Further, the transfer function for the restricted area and the modal coefficients are
derived via Eqs. (44) and (45). Here is where the approach differs from Ref. [12], where the relative error is estimated via the
solution of (46) for n and q, as it is introduced in 4.1. Therefore, we estimate the relative error between actual transfer functions of the ROM, instead of the block matrix representation W n . In the following we choose q ¼ n m p, which results in
an error estimation at every second order of the matrix-valued Padé-approximation (see Eq. (47)).
The adaptive windowing method is based on the convergence of the MPVL algorithm by further advancing Lanczos iterations. In finite arithmetic this convergence can not be guaranteed. For this reason, look-ahead and deflation methods, based
on Aliaga et al. [30] are implemented, as mentioned in Section 4.2. Also, the method is similar to the approach in Ref. [16],
where the frequency response function is compared for two subsequent ROM.
In Fig. 3 the estimated errors en are shown over the frequency, for different expansion points f p . The expansion points are
selected in a heuristic manner, following the approach in Ref. [12]. In the illustrated example, the algorithm starts with an
expansion frequency at f 1 ¼ 350 Hz and generates the Krylov subspaces via the Lanczos-type algorithm. At each m þ p iteration, the relative error en is computed. This is continued until one of the following cases is determined:
1. In the Lanczos-type process a breakdown occurs, due to division by zero or close to zero.
2. The Krylov subspace is exhausted, due to deflation of linear dependent or nearly linear dependent Lanczos vectors.
3. The adaptive windowing method detects stagnation, hence the size of the frequency window is not increasing any longer,
by generating new Lanczos vectors.
If one of the above cases occurs, a ROM is obtained, which approximates the original system in the frequency range
f min;p ; f max;p around the expansion point f p . Unless, the frequency window is smaller than the frequency range of interest,
a new expansion point is selected and a new Krylov subspace generated. Since, the choice for the expansion point is heuristic,
an overlap of frequency windows occurs. This can be seen in Fig. 3. The transfer function of the ROM in the overlapping area,
should estimate a similar value, at least in the range of the user-defined limit for en . To ensure this, we now introduce a
weighted adaptive windowing algorithm, which applies a weighting function to the results of the MPVL algorithm for each
estimated frequency window. As already shown, the accuracy of the Padé approximation around an expansion frequency is
decreasing with increasing distance to the expansion point. Therefore, a weighting function which decreases towards its
C. Sittl et al. / Mechanical Systems and Signal Processing 148 (2021) 107135
11
Fig. 3. Relative error en plotted over frequency for each expansion point f p in the frequency range. The error bound for the relative error is 1e-3. Further, the
sum of all weighting functions is illustrated.
upper and lower bound is appreciated. Hence, the solutions near the expansion point contribute more significantly than the
solution near the window edges. This is due to the increasing error distant from the expansion point. The weighting function
wp , which is used in the present approach is defined by:
8
a
a
< sin2 pf 1þ
p N
for f min;p 6 f p 6 f max;p ;
N 1
wp ¼
:
0
otherwise;
ð59Þ
where N is the estimated window width and a a parameter, which shifts the maximum of the window function to the expansion point f p and is defined as:
!
N 1 1
a ¼ ln
:
f
2f p
ln Np
ð60Þ
The weighting function has a maximum value of 1 at the expansion point and decays towards 0 at the window edges. The
weighting function is similar to a von Hann window, which is used in signal processing applications, but its maximum is
shifted with a, to hold the above properties. Hence, the function is not symmetrical to the maximum, because in most cases
the left and right edges of the estimated windows are at different distances to the expansion point. In Fig. 3 the sum of all
window functions wp is illustrated. Depending on the overlap of the frequency windows, the sum of the weighting function
can exceed 1. By assessing the sum of this weighting function, further expansion points can be chosen, where the sum is
beneath a lower bound. Thus, generating a specific amount of overlap, which leads to a more robust method. In the illustrated example in Fig. 3, an additional expansion point at 1400 Hz can be selected and a supplementary ROM generated. Furthermore, the weighting function ensures smoothness between the different windows.
6. Structural damping
Incorporating damping is limited to constant non-proportional damping models, due to the MPVL method, used in the
present work. If more general damping models are of interest, a doubling of the systems dimension is provoked in the present approach, due to the reformulation of the second-order problem in a system of first order equations. In Ref. [23] it is
shown, how to incorporate Rayleigh damping for the Lanczos method. For more general damping models second-order
Arnoldi-based methods, can be considered, which circumvent the drawback of doubling the systems dimensions [22]. Furthermore, the effort, for generating the orthogonal basis, via an Arnoldi-based method, is increasing, with accumulation of
the vectors, compared to a Lanczos-based method. Hence, the Lanczos method is used in the present approach. In the present
work we consider constant structural damping, where cs is used as the constant damping ratio on the structural part of the
system as:
^s ¼
K
"
ð1 þ ics ÞK s
0
RT
Kf
#
;
ð61Þ
12
C. Sittl et al. / Mechanical Systems and Signal Processing 148 (2021) 107135
^ s is now complex valued, which has no impact on the number of operations
where i is the imaginary unit. Note, the matrix K
and storage requirements in the present work, since the system is already complex valued, due to the modified DtN boundary condition. The constant damping factor is chosen to be frequency independent and is therefore added to the frequency
independent part in Eq. (42) as:
^s
A0 ¼ K
^ þB
^ 1 ðx0 Þ:
x20 M
ð62Þ
and the MPVL algorithm can be applied, without any further modifications. In addition, damping can cluster the eigenvalues
of the matrix A0 , which can be favorable for the Lanczos-type method. Therefore, damping can lead to broader estimated
window width and thus improves the computation time.
7. Numerical examples
In the following, we show numerical results for the sound wave radiation from a plate. The presented MPVL algorithm,
including the weighted adaptive windowing technique and structural damping, is applied to evaluate the example shown in
Fig. 4 consisting of a plate in air mounted in the center of the plate on a vibration exciter, see also Fig. 8. Due to a harmonic
structural load, applied at the center of the coordinate system, the plate radiates acoustic waves. The amplitude for the harmonic load is constant throughout the frequency range. The material of the plate is aluminum with the following properties:
Young modulus is E ¼ 69 GPa, the density is qs ¼ 2770 kg/m3, the Poisson ratio is m ¼ 0:33 and the constant structural damping factor is set to cs ¼ 0:01. For the acoustic medium the following material properties for air are used: The density is
qf ¼ 1:225 kg/m3 and the speed of sound is c ¼ 343 m/s.
The dimensions of the plate are 376 376 2 mm and the radius of the sphere is R ¼ 270:9 mm. Note, that the distance
between the structural sound source and the truncating DtN surface can be chosen to be very small. For the finite element
discretization a maximum element edge length h ¼ 25 mm is employed, which provides a minimum of 9 elements per wavelength in the fluid domain. The fluid domain is discretized with quadratic tetrahedral elements with 10 nodes each and the
structural domain is discretized with hexahedral elements. This results in a system with 139,511 degrees of freedom. The
frequency range of interest is set to f 2 ½100; 1500 Hz in steps of Df ¼ 5 Hz. For the solution, point 1 on the structure, point
2 on the spherical surface and point 3 in the exterior are selected, see Fig. 4. In the following, the relative error tolerance
en ¼ 1 10 3 and the tolerances for look-ahead and deflation are set to 1 10 10, which is a magnitude of 100 times smaller than the spectral norm jjAjj2 , which can be estimated in the first few runs during the Lanczos-type algorithm. The presented method is implemented in MATLABÒ. For the MPVL method, the main costs are the LU-factorization for each ROM
and the forward- and backward substitutions in each iteration of the Lanczos process. The LU-factorization is carried out
with row and column permutations to keep the fill-in small. Moreover, a diagonal scaling for A0 is incorporated to improve
the convergence of the method [13]. Furthermore, depending on the expansion point, A0 has a condition number of around
1 1016.
Fig. 4. Numerical example, which illustrates the aluminum plate surrounded by air The points 1, 2 and 3 represent the locations where the solution is
computed.
C. Sittl et al. / Mechanical Systems and Signal Processing 148 (2021) 107135
13
7.1. Near-field solutions
In the following we compare the solutions for point 1 and 2, obtained via the presented method, with a solution computed with ANSYSÒ. Since, as far as we know, the DtN boundary condition is not implemented in ANSYSÒ, we selected
the PML method for the boundary condition on the spherical surface, because it is a well known and common method
and therefore implemented in commercially used simulation software. The thickness of the PML regions are chosen in
two variants as 200mm with 9 layers through the thickness and 400 mm with 17 layers. Further, to keep the results comparable, the system matrices are exported from ANSYSÒ and used in our MATLABÒ code. Thus, the discretization inside the
spherical boundary is identical, and so are the system matrices. The results for the PML boundary condition are computed by
a solution at every frequency of interest, with a direct solver. For the MPVL method three different numbers of series terms
are computed, setting N DtN ¼ 4; 8 and 10. During the three computations, the number of ROMs, constructed via the MPVL, is
4, 2 and 5 respectively for the number of series terms. With the weighted adaptive windowing scheme, these ROMs are combined for the solution over the frequency range.
In Fig. 5 the displacement amplitude and phase for point 1 are illustrated over the frequency for the MPVL method and the
PML solution. An excellent agreement of the results can be observed. Due to the elements in the PML region the systems
dimension is increasing, which is not the case for the DtN boundary condition. Due to the increased number of degrees of
freedom in the PML solution, a computational advantage of the proposed method is to be expected. This will be investigated
in an upcoming publication.
The solution for point 2, which is located on the spherical surface at h ¼ p=4 and / ¼ p=2 is illustrated in Fig. 6, for both
methods. The PML region has to be sufficiently large, in comparison to the wavelength. Therefore, the number of elements in
the PML region can be very large, if one discretized model should be used over a wide frequency range, which is not the case
for the DtN boundary condition. The results for the transfer function amplitudes fit well over the frequency range. Further,
the number of series terms does not show significant variations for the amplitudes. In terms of the phase, the results show
significant variations, not only for different numbers of series terms, but also compared to the PML solutions. A reason for
this deviation is the local boundary condition B1 ð pÞ, which is introduced in Section 2.3, since it is kept constant for each computed expansion point (see Section 4). For the given problem, the derivative of the absolute value of the local boundary condition is nearly constant for frequencies above 500 Hz, but decreasing in the range from 500 Hz to 100 Hz. Therefore, an
approximation error occurs in this frequency range.
In addition the results are obtained directly on the enclosing surface, where either PML or the DtN-map are applied.
Hence, the discrepancies are related more to the differences in these boundary conditions, instead of the solution process
(either the direct solution or the MPVL). For this reason, the results in Fig. 5, obtained in the interior of the sphere, are matching well.
As a result of the weighted adaptive windowing scheme, the results preserve smoothness over the evaluated frequencies.
Computational times are exemplified in Fig. 7 in relation to the number of computed frequencies in the given frequency
range. The time is scaled to the maximum computation time for the maximum number of frequencies. The computations are
made with the same CPU (Intel Xeon E-2144G). Recognizable is that a speed up with the presented method is achieved, if the
number of computed frequencies exceeds a certain limit, which is dependent on the number of series terms for the DtN
Fig. 5. Comparison of the PML solution with 200mm thickness and the MPVL method with N DtN ¼ 8 for the displacement amplitude and phase at point 1 on
the structure.
14
C. Sittl et al. / Mechanical Systems and Signal Processing 148 (2021) 107135
Fig. 6. Comparison of PML solutions and the MPVL method with N DtN ¼ 4; 8; 10 for the sound pressure level and phase at point 2 on the spherical surface.
Fig. 7. Comparison in computation time over the number of computed frequencies.
boundary condition. The offset at the origin for the MPVL method is due to the construction of the ROMs. Therefore, the presented method is effective, if an appropriate number of frequencies is of interest.
7.2. Acoustic measurement
Further, we compare the exterior solution, obtained with the MPVL method, with an acoustic measurement. In this setup,
the aluminium plate is excited via a vibration exciter, mounted to the center of the plate, which is illustrated in Fig. 8. The
excitation force is estimated by a measured acceleration signal and the mass of the vibration exciter iron core. The solution is
Fig. 8. Measurement setup in the anechoic chamber. The aluminium plate is mounted on the vibration exciter and the acoustic pressure signal is recorded
via a calibrated condenser microphone.
C. Sittl et al. / Mechanical Systems and Signal Processing 148 (2021) 107135
15
Fig. 9. Computed and measured sound pressure levels at point 3.
evaluated at point 3, which is located at h ¼ p=4; / ¼ p=2 and radius r ¼ 300 mm. Since the solution on the boundary is
obtained via the modal coefficients, the solution at any point outside the sphere can be computed with Eq. (11), in a post
processing step. In Fig. 9 the computed sound pressure level via the MPVL method and the measured sound pressure level
are illustrated. The method predicts the measured response with very good accuracy except before the first illustrated eigenfrequency. This discrepancy occurs, due to the non-linear frequency response of the vibration exciter below 200 Hz.
8. Conclusion
In this article, we introduced a weighted adaptive windowing method, which is applied to the matrix-valued Padé-viaLanczos algorithm for multi-frequency solutions, including exterior structural acoustic domains. The method is suitable
where only a reduced number of solutions are of interest, which we denoted as partial field solution, whether in the
near-field or in the exterior domain. Applications include virtual sensing techniques or problems, where the partial field
location is known beforehand, as the drivers ear position in an automobile.
With the proposed adaptive windowing technique, multiple Krylov subspaces can be constructed which makes it well
suited for parallelization, further increasing the numerical efficiency of the approach. The weighted adaptive windowing
does not only estimate the frequency range for each ROM, but also allows for a simple summation of these frequency windows, preserving smoothness between the frequency windows. Further, the sum of the weighting functions can be assessed
for estimating the reliability of the generated ROM, within the complete frequency range of interest. This is achieved by
assuming a specific amount of overlap, thus the sum of the weighting functions can be set to a lower bound.
In the present approach, hysteretic structural damping is incorporated. Due to the property of the structural damping, it
can be added without further modification to the MPVL method. In conjunction with a diagonal scaling of the frequency
independent matrix A0 , this contributes to the stability of the method by clustering the eigenvalue distribution. Therefore,
a better conditioned matrix is obtained, as also pointed out in [13]. Balancing the system matrix A directly, is computationally expensive and numerically unstable, due to the inversion of the poorly conditioned matrix A0 . In a forthcoming publication more general damping methods in conjunction with the DtN boundary condition will be investigated, incorporating
these in the present approach, as well as applying a second-order Arnoldi scheme [22].
By comparing the proposed method to the PML method, which solves the system with a standard frequency sweep technique, good agreement for the displacement amplitudes on the structure is obtained. But there is a significant speedup in
solution times with the proposed method, since the system matrices need not to be factorized for every frequency in the
desired range. The speedup is achieved, if the number of frequencies which are computed, is sufficiently large.
The obtained results of the MPVL algorithm correlate well with measured pressure values, radiated from a plate.
CRediT authorship contribution statement
Christopher Sittl: Software, Data curation, Writing - original draft, Writing - review & editing. Steffen Marburg: Supervision. Marcus Wagner: Funding acquisition, Methodology, Writing - review & editing.
Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have
appeared to influence the work reported in this paper.
16
C. Sittl et al. / Mechanical Systems and Signal Processing 148 (2021) 107135
Acknowledgements
This work was supported by the Federal Ministry for Economic Affairs and Energy on the basis of a decision of the German
Bundestag. This support is gratefully acknowledged.
References
[1] M.M. Wagner, P.M. Pinsky, A.A. Oberai, M. Malhotra, A Krylov subspace projection method for simultaneous solution of Helmholtz problems at
multiple frequencies, Comput. Methods Appl. Mech. Eng. 192 (41–42) (2003) 4609–4640, https://doi.org/10.1016/S0045-7825(03)00429-8.
[2] S. Marburg, Developments in structural-acoustic optimization for passive noise control, Arch. Comput. Methods Eng. 9 (4) (2002) 291–370, https://doi.
org/10.1007/BF03041465.
[3] S. Marburg, Normal modes in external acoustic. Part III: Sound power evaluation based on superposition of frequency-independent modes, Acta Acust.
united Ac 92 (2006) 296–311.
[4] L. Moheit, S. Marburg, Normal modes and modal reduction in exterior acoustics, J. Theor. Comput. Acout. (2018) 1850029, https://doi.org/10.1142/
S2591728518500299.
[5] S. Marburg, F. Dienerowitz, D. Fritze, H.-J. Hardtke, Case studies on structural-acoustic optimization of a finite beam, Acta Acust united Ac 92 (3) (2006)
427–439.
[6] Herwig Peters, Nicole Kessissoglou, Steffen Marburg, Modal decomposition of exterior acoustic-structure interaction problems with model order
reduction, J. Acoust. Soc. Am 135 (2014) 2706–2717.
[7] A.S. Wixom, J.G. McDaniel, Fast frequency sweeps with many forcing vectors through adaptive interpolatory model order reduction, Int. J. Numer.
Methods Eng. 100 (6) (2014) 442–457, https://doi.org/10.1002/nme.4742.
[8] S.K. Baydoun, M. Voigt, C. Jelich, S. Marburg, A greedy reduced basis scheme for multifrequency solution of structural acoustic systems, Int. J. Numer.
Methods Eng. 121 (2) (2020) 187–200, https://doi.org/10.1002/nme.6205.
[9] W.B. Gragg, Matrix interpretations and applications of the continued fraction algorithm, Rocky Mt. J. Math. 4 (2) (1974) 213–226, https://doi.org/
10.1216/RMJ-1974-4-2-213.
[10] R.W. Freund, Computation of matrix Padé approximations of transfer functions via a Lanczos-type process, in: C.K. Chui (Ed.), Approximation Theory
VIII, Series in Approximations and Decompositions, World Scientific, Singapore, 1995, pp. 215–222.
[11] M. Malhotra, P.M. Pinsky, Efficient computation of multi-frequency far-field solutions of the Helmholtz euquation using Padé approximation, J.
Comput. Acoust. 08 (01) (2000) 223–240, https://doi.org/10.1142/S0218396X00000145.
[12] J.P. Tuck-Lee, P.M. Pinsky, Adaptive frequency windowing for multifrequency solutions in structural acoustics based on the matrix Padé-via-Lanczos
algorithm, Int. J. Numer. Methods Eng. 73 (5) (2008) 728–746, https://doi.org/10.1002/nme.2102.
[13] H.-L. Liew, P.M. Pinsky, Matrix-Padé via Lanczos solutions for vibrations of fluid-structure interaction, Int. J. Numer. Methods Eng. 84 (10) (2010) 1183–
1204, https://doi.org/10.1002/nme.2936.
[14] J. Baumgart, S. Marburg, S. Schneider, Efficient sound power computation of open structures with infinite/finite elements and by means of the Padévia-Lanczos algorithm, J. Comput. Acoust. 15 (4) (2007) 557–577, https://doi.org/10.1142/S0218396X07003494.
[15] S. van Ophem, O. Atak, E. Deckers, W. Desmet, Stable model order reduction for time-domain exterior vibro-acoustic finite element simulations,
Comput. Methods Appl. Mech. Eng. 325 (2017) 240–264, https://doi.org/10.1016/j.cma.2017.06.022.
[16] A. van de Walle, The power of model order reduction in vibroacoustics, PhD Thesis, KU Leuven, Leuven (2018)..
[17] J.-P. Berenger, A perfectly matched layer for the absorption of electromagnetic waves, J. Comput. Phys. 114 (114) (1994) 185–200.
[18] D. Givoli, I. Patlashenko, J.B. Keller, Discrete Dirichlet-to-Neumann maps for unbounded domains, Comput. Methods Appl. Mech. Eng. 164 (1–2) (1998)
173–185, https://doi.org/10.1016/S0045-7825(98)00053-X.
[19] F. Ihlenburg, Finite Element Analysis of Acoustic Scattering, vol. 132 of Applied Mathematical Sciences, Springer-Verlag, New York Inc, New York, NY,
1998, doi: 10.1007/b98828..
[20] J.-P. Coyette, C. Lecomte, J.-L. Migeot, J. Blanche, M. Rochette, G. Mirkovic, Calculation of vibro-acoustic frequency response functions using a single
frequency boundary element solution and a Padé expansion, Acoustica 85 (3) (1999) 371–377.
[21] B. Salimbahrami, B. Lohmann, Order reduction of large scale second-order systems using Krylov subspace methods, Linear Algebra Appl. 415 (2–3)
(2006) 385–405, https://doi.org/10.1016/j.laa.2004.12.013.
[22] Z. Bai, Y. Su, Dimension reduction of large-scale second-order dynamical systems via a second-order Arnoldi method, SIAM J. Sci. Comput. 26 (5) (2005)
1692–1709, https://doi.org/10.1137/040605552.
[23] K. Meerbergen, Fast frequency response computation for Rayleigh damping, Int. J. Numer. Methods Eng. 73 (1) (2008) 96–106, https://doi.org/10.1002/
nme.2058.
[24] R.S. Puri, D. Morrey, A Krylov-Arnoldi reduced order modelling framework for efficient, fully coupled, structural–acoustic optimization, Struct.
Multidiscip. Optim. 43 (4) (2011) 495–517, https://doi.org/10.1007/s00158-010-0588-5.
[25] M.J. Grote, J.B. Keller, On nonreflecting boundary conditions, J. Comput. Phys. 122 (2) (1995) 231–243, https://doi.org/10.1006/jcph.1995.1210.
[26] J.B. Keller, D. Givoli, Exact non-reflecting boundary conditions, J. Comput. Phys. 82 (1989) 172–192.
[27] M. Malhotra, Iterative solution of large-scale acoustics problems, PhD Thesis, Stanford University, 1996..
[28] G.H. Golub, C.F. van Loan, Matrix Computations, third ed., Johns Hopkins studies in the mathematical sciences, Johns Hopkins Univ. Press, Baltimore,
1996..
[29] M.M. Wagner, P.M. Pinsky, M. Malhotra, An efficient algorithm for the simultaneous solution of exterior acoustic problems over a frequency band, J.
Acoust. Soc. Am. 110 (5) (2001) 2735, https://doi.org/10.1121/1.4777490.
[30] J.I. Aliaga, D.L. Boley, R.W. Freund, V. Hernández, A Lanczos-type method for multiple starting vectors, Math. Comput. 69 (232) (2000) 1577–1602,
https://doi.org/10.1090/S0025-5718-99-01163-1.