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 : p
n¼ ð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.