Parametric model order reduction for acoustic metamaterials based on local thickness variations

The use of lightweight materials is obligatory in the design of economic structures, but with decreasing mass, the vibro-acoustic properties of the structures become unfavorable. Acoustic metamaterials based on the concept of acoustic black holes (ABH) can improve the vibrational behavior by adding local thickness variations to the host structure. To make computational models of metamaterials also valid for high frequencies, a fine discretization has to be used. This leads to high computation times. A parametric model order reduction (PMOR) approach is presented, which is able to create an efficient reduced model of ABHs. In the first step, one structure preserving reduced model for each parameter is generated using the iterative rational Krylov algorithm (IRKA). The transfer functions of all reduced models are used in the second step to create a single parametric model using the Loewner framework. The resulting model can be used to efficiently evaluate the frequency response of a structure equipped with ABHs for different parameter sets.


Introduction
Lightweight materials are demanded in many fields of engineering, ranging from automotive and aerospace to civil engineering applications, but decreasing the mass of a system comes at the price of introducing resonance or noise problems. By applying arrays of locally resonant substructures of relatively low mass, compared to the host structure, the overall vibration response of the structure can be modified and resonances can be avoided. A structure equipped with such arrays is called acoustic metamaterial. Different concepts exist and show interesting effects, for example negative effective density or stiffness [1], which can generate stop bands at specific frequencies [2]. A different approach to influence the wave propagation in lightweight materials is based on local thickness variations. This contribution focuses on structures equipped with acoustic black holes (ABH), which are areas where the structure's thickness is gradually reduced, theoretically until it vanishes. With reducing thickness, the vibration amplitudes increase while the wave speed decreases. Because of the high amplitudes in the ABH region, damping material applied here can effectively dissipate much of the structure's vibration energy [3]. Arrays of ABHs can be considered as acoustic metamaterials [4] and have shown to have a great effect on the radiated sound power of plates [5]. In this work, we want to investigate the possibility to generate efficient models of plates equipped with ABHs which can be used in the design process of an acoustic metamaterial.
To be able to computationally model the frequency response of structures validly also for high frequencies, a very fine discretization has to be used. This results in large and computationally expensive models. An efficient handling of such models is often not possible, so methods to RASD IOP Conf. Series: Journal of Physics: Conf. Series 1264 (2019) 012014 IOP Publishing doi: 10.1088/1742-6596/1264/1/012014 2 reduce the order of the model while preserving the desired frequency response are sought. Various methods for model order reduction (MOR) exist and are used in different fields [6]. The methods approximate the full system's transfer function using either approaches based on singular value decompositions (SVD) or moment-matching methods based on Krylov subspaces; the latter will be considered in this contribution. There exist generalizations of the classic Krylov methods for differential equations of second order, so that the second order structure of the dynamic system is preserved [7]. This is very desirable, as the reduced models can be used, for example, to compute the energies in the system, which would not be possible if the general structure is lost in the reduction process. Additionally, the conversion of the second order system to an equivalent first order system would double the initial system size, counteracting the goal of reducing the system size. Contrary to classic reduction methods like the modal superposition, Krylov based MOR techniques can be used with arbitrarily damped systems or systems with complex material moduli [8,9]. Having a reduced model depending on one or more parameters from the full model is a key ingredient for the efficient optimization of systems. Many parametric model order reduction (PMOR) methods exist and some of them have also been extended to second order systems [10,11]. In this contribution, a PMOR approach based on the Loewner framework, which generates parametrized reduced models from transfer function measurements [12], is combined with the iterative rational Krylov algorithm (IRKA) [13]. Here, the transfer function measurements are obtained from the system reduced by IRKA.
The paper is structured as follows: The considered MOR methods, the iterative rational Krylov algorithm (IRKA) for second order systems and the parametric Loewner framework, are presented in section 2. Based on this, the proposed combination of IRKA and the parametric Loewner framework is presented in section 3. Section 4 describes the concept of acoustic black holes in detail. Numerical examples for all MOR methods and the application to design acoustic black holes are presented in section 5.

Efficient modeling of dynamic systems
The physical phenomena of interest are modeled as second order dynamic systems of the form with M, C, K ∈ R n×n , u, f, y, b ∈ R n , n is the dimension of the system, ω the exciting frequency, and i the imaginary unit. M, C, and K describe the system's mass, damping, and stiffness characteristics, f are the external forces, and u the resulting displacements. y is the output of the system and b the output measurement vector. The transfer function of such systems is where M r , C r , K r ∈ R r×r , u r , f r , y r , b r ∈ R r , which approximates the original system. The reduced system's transfer function is defined as H r (ω) = b T r −ω 2 M r + iωC r + K r −1 f r and the reduced model is suitable, if H r (ω) = H(ω) in the desired frequency range. Moment matching methods seek a reduced transfer function H r (ω) rationally interpolating the original transfer function H(ω) by matching a number of poles (moments) and their derivatives. A realization of the reduced system can be constructed using Petrov-Galerkin projections, where the full order system is projected onto r-dimensional subspaces V, W with bases V r , W r ∈ R n×r [14]. The reduced system is then given by The subspaces V, W, respectively their bases V r , W r can be found using various methods; [6] provides an overview. In this work, we use rational interpolation based on Krylov subspace methods. The general procedure for first order systems is discussed in [15]. For the proposed second order system, the approach has to be adopted to the structure of the problem [7,16,17,18].
Large models depending on several parameters further increase the computational effort in the design process. The second order system of (1) depending on a set of parameters p i , i = 1, . . . , m can be expressed as It is desirable to reduce parametrized models retaining their dependence on p. Various MOR methods have been extended for the use of parametrized models, for an overview see [10]. Parametrized reduced order models are used, where models have to be evaluated for many values of certain parameters. The construction of such a reduced model is more expensive than using a standard MOR method and likely also more expensive than solving the full system for one parameter. The benefit of such a method is therefore depending on the number of different parameter configurations, which have to be evaluated. In the proposed approach, we want to generate parametric reduced order models of acoustic metamaterials using a two-step approach: First, a reduced model for each parameter is created using a rational Krylov method. The frequency response of these models is calculated and in the second step used as input data for the parametric Loewner framework presented in [12]. The result is a single reduced model depending on the desired set of parameters.

The iterative rational Krylov algorithm (IRKA)
The iterative rational Krylov algorithm (IRKA) was originally presented by [13] as a method to generate H 2 optimal reduced order models of first order dynamic systems. The algorithm starts with an initial selection of expansion points s i , i = 1, . . . , r. Based on this selection, a reduced order model is generated. The eigenvalues λ i , i = 1, . . . , r of this reduced system are mirrored across the imaginary axis and then used as expansion points for the next iteration. Upon convergence, the resulting reduced system is H 2 optimal. The algorithm has been extended for second order systems for example by [18,19]. We used a slightly modified approach of the SO-IRKA algorithm presented in [18]. Our implementation is shown in algorithm 2.1.

Algorithm 2.1
The iterative rational Krylov algorithm for second order systems, SO-IRKA Require: M, C, K, f , initial expansion points s i , i = 1, . . . , r ∈ C, and tol Ensure: Solve the eigenvalue problem of second order M r λ 2 + C r λ + K r x = 0 6: Choose r eigenvalues from λ j , j = 1, . . . , 2r 7: Update expansion points s i = −λ r 9: end while 10: We want to preserve the symmetric positive definiteness of the input matrices, so we set W r = V r as proposed in [18]. The computationally expensive part of the algorithm is the solution of r n × n linear equation systems per iteration; the eigenvalue problem is cheap, as it is performed on the reduced system. Note, that the quadratic eigenvalue problem produces 2r potential expansion points, so a subset of them has to be chosen. If the eigenvalues closest to the imaginary axis are selected, the first r poles of the original transfer function can be approximated. It is, however, also possible to choose eigenvalues in a specific range of interest, for example to match poles in a certain frequency range. Upon convergence, the transfer function H r (ω i ) of the reduced system interpolates the transfer function of the full system H(ω) for r poles.

Parametric model order reduction using the Loewner framework
In the following section, the Loewner framework introduced in [20] and its extension to parametric models introduced in [12] is presented. IRKA and the non-parametric Loewner framework both have a second order realization [18,21], while the parametric Loewner framework does not have such a realization. The Loewner framework is a data driven approach based on rational interpolation, which generates a reduced model based on N transfer function measurements. Contrary to the aforementioned Krylov methods, where the quality of the reduced model depends on the choice of a suitable set of expansion points s i , sets of data (i.e. transfer function measurements) exhibiting the important features of the full order model are required for the Loewner framework to generate proper reduced models. The framework finds a rational function that interpolates a number of transfer function measurements H (s), i.e. r (s i ) = H (s i ) [22]. With the Loewner framework, it is possible to generate reduced systems without knowledge of the full system (data driven model order reduction). A rational function (5) of order k can be constructed, given some expansion points λ i , using the rational barycentric formula The Loewner framework is now used to find the coefficients α i , β i . First, transfer function measurements are concatenated in the Loewner matrix L. Given N transfer function evaluations φ i at frequencies s i , they are partitioned into two sets, with n + n = N and n = N/2 , which populate the Loewner matrix [22] The partitioning in (7) can be arbitrary. However, [23] shows, that an alternating partition leads to better conditioned interpolation matrices. A vector c satisfying Lc = 0 contains the coefficients for (6) with α i = c i , β i = c i w i . The rank of matrix L holds information about the minimal order of the interpolation problem. A reduced model of order k = rank L + 1 fully interpolates all information contained in the transfer function measurements, while a reduced order model of smaller order approximates the measurements [12]. The Loewner framework has been extended for parametric systems in [12]. The barycentric interpolation formula (6) is extended to the two variables frequency s with expansion points λ i and parameter p with expansion points π j :

RASD
The parameters c ij are obtained as follows: First, transfer function measurements for a set of parameters p j , j = 1, . . . , m at frequencies s i , i = 1, . . . , n have to be generated and partitioned into a matrix φ ij = H (s i , p j ): [s 1 , · · · , s N ] = [λ 1 , · · · , λ n ] ∪ µ 1 , · · · , µ n , The transfer function samples φ ij can have various origins: results from numerical simulations can be used as well as measurement data or a combination of both. We now determine k = max j rank L p j from the Loewner matrices L p j associated with the columns of φ ij and q = max i rank L s i from the Loewner matrices L s i associated with the rows of φ ij . A function of order k +1, q +1 interpolates all data given in φ ij . Now, the two-variable Loewner matrix L 2 can be computed as in [12, (4.10)] and its null space contains the coefficients for (9), i.e. L 2 c = 0. If no null space is found, the coefficients can be obtained from a singular value decomposition of L 2 . In this case, the reduced model approximates the measurements and the approximation error is proportional to the smallest singular value of L 2 . Note, that of course only data contained in the measurements can be interpolated and a too coarse sampling in either frequency or parameter domain results in reduced models of bad quality compared to the full system response. The number of transfer function measurements N and M can not be directly related to the order of the full model n, but can be adjusted using the ranks k and q of the Loewner matrices. If k + 1 = n (respective q + 1 = m), not all features of the transfer function can be kept in the reduced model and a finer grind in frequency (respective parameter) space should be considered.

Parametric model order reduction using IRKA and the Loewner framework
A main issue in parametric model order reduction using the Loewner framework is the need for many transfer function evaluations. In cases they cannot be obtained from experimental measurements, they are computed numerically. This can take a serious amount of time, so it is advisable to use already reduced models (which are parameter independent) in this step. We therefore combine IRKA and the parametric Loewner framework as suggested in [24]. The idea is to generate reduced models for a set of parameters using IRKA, and then use the reduced models to compute the frequency response for each parameter, which are then used as input for the parametric Loewner framework. The convergence for the different IRKA runs can be accelerated, if the expansion points used at convergence for one parameter are reused as starting point for the following parameters. Often, the overall system response does not change drastically, if a parameter is slightly changed. Figure 1 shows transfer function evaluations of cantilevered beams of length l = 0.8 m with quadratic cross sections of different heights a i = [0.005, . . . , 0.05] m. The system is exposed to a harmonic loading at the free end. It can be observed, that the poles of two parameters next to each other lie in a similar region. So even while the best IRKA expansion points for one parameter are not necessarily the best for a second parameter, the reused expansion points are a good start for the following iterations and the algorithm converges faster. The used version of IRKA is structure preserving, i.e. results in reduced mass, damping, and stiffness matrices M r , C r , K r . These matrices are used to compute the transfer function measurements, which are the basis for the Loewner framework. It is also possible, for example, to compute the frequency dependent potential and kinetic energy in the model and use this data as input to the Loewner framework. The parametric model resulting from the Loewner framework does not have second order structure, so post-processing methods requiring the system matrices have to be computed prior to applying the Loewner framework. The procedure is shown in algorithm 4.1.

Acoustic metamaterials based on local thickness variations
Thickness variations in thin-walled structures have a strong influence on the propagation of bending waves. Theoretically, a wave traveling through a structure with decaying wall thickness comes to a halt when the structure thickness smoothly reduces to zero. These local thickness variations are known as acoustic black holes (ABH) and an optimal thickness profile has been proposed by [25]. Practically, the thickness of a structure can not vanish completely, but due to larger amplitudes in the ABH region, more energy can be dissipated here compared to a uniform plate. The ABH can also be used to focalize the bending waves to certain regions of the structure, where they can be damped efficiently [26]. Applying damping material inside the ABH region has been proven to be an efficient way of reducing the vibration of the structure as well as the structure-borne noise [4,27]. The damping effect can be increased, if the damping material is constrained by a more rigid material. Figure 2 shows the total kinetic and potential energies, normalized with respect to the input power of a titanium plate (E = 1.04e11 N/m 2 , ν = 0.31, ρ = 4430 kg/m 3 ) under harmonic loading compared to the same plate equipped with an ABH with constrained layer damping (CLD). The damping material properties are E = 1.0e7 N/m 2 , ν = 0.45, ρ = 1000 kg/m 3 . A material loss factor η = 0.001 is assumed for the titanium parts in both models, the damping layer has η = 0.1. The damping Call SO-IRKA (Algorithm 2.1) with M, C, K, f, s i to obtain reduced system matrices M r , C r , K r , f r and updated expansion points s i

4:
Compute the frequency response function in the frequency range of interest of the reduced model using ω n and (2) and save it in the measurement matrix φ nm 5: end for 6: Partition ω and p: [ω 1 , · · · , ω N ] = [λ 1 , · · · , λ n ] ∪ µ 1 , · · · , µ n , [p 1 , · · · , p M ] = [π 1 , · · · , π m ] ∪ ν 1 , · · · , ν m 7: Construct the parametric Loewner matrix L 2 as shown in [12, (4.10)] 8: Compute c so that L 2 c = 0 9: Use the two-variable barycentric formula (9)   It can be observed, that the normalized total energy in the plate is lower in the complete frequency range. Thus, a considerable reduction of energy in the plate can be achieved using only a relatively low amount of damping material and without significantly increasing the total mass of the system. It is still an open question, where ABHs should be placed on a structure to optimally reduce the energy contained in the system. Also the properties of the damping material have an influence on the resulting energies. Due to small wavelengths corresponding to higher frequencies, the numerical models need a fine discretization and the computations are expensive. Reduced models depending on some parameters, for example location or damping layer thickness, can help reducing the computation time and therefore a more efficient design process is possible.

Numerical examples
In the following, we will present numerical experiments using IRKA, the parametric Loewner framework, and the combination thereof. IRKA and the parametric Loewner framework will be tested on a model of a clamped beam, the combined framework will also be used to reduce models of acoustic black holes to show the performance of the method in the design process. All numerical experiments are computed on a Linux system using an Intel R Xeon R W-2135 CPU at 3.70 GHz and 32 GB RAM. All codes are implemented with MATLAB R version R2018a, the FEM matrices for the ABH models were obtained from KRATOS Multiphysics [28].

Beam model
To demonstrate the capabilities of the MOR methods, a cantilevered beam is discretized using finite element modeling. We are using proportional damping C = αM + βK; the symmetric positive definite matrices have a dimension of 300 × 300. The beam has a length of l = 0.8 m and a quadratic cross section with height a = 0.01 m; material parameters used for discretization are E = 2.1e11 N/m 2 , ν = 0.3, ρ = 7850 kg/m 3 . The free end is loaded with 1 N, the frequency response functions H(ω) and H r (ω) show the displacement at the tip. Using IRKA, the model is reduced to a dimension of 14 × 14 with error norm ||H(ω)−Hr(ω)|| ||H(ω)|| = 1.153e−6 in the frequency range ω = 10 1 , . . . , 10 5 rad/s. Figure 4 shows the frequency response function and the error. The algorithm converges within 3 iterations for initial expansion points distributed logarithmically in the frequency range of interest. The computation times are t f = 0.243 s for the full model and t r = 0.056 s for the IRKA reduced model (reduction and evaluation time). Given the small errors, IRKA provides a reliable method to reduce the beam model.
The beam model described above is now parametrized with respect to the cross section height a, all other parameters remain constant. Frequency responses for a set of 12 linearly distributed cross-section sizes a i = [0.005, . . . , 0.05] m evaluated at 400 frequencies logarithmically distributed between 10 and 10 5 rad/s are given as input to the parametric Loewner framework. The framework returns the coefficients for (9), which can be evaluated for any frequency and parameter value in range of the snapshot data. The reduced model has dimensions of k = 66 in the frequency space and q = 5 in the parameter space. The frequency response of the reduced model is now evaluated for a cross section height of a = 0.0425 m, which is exactly in the middle between two parameters in the set the reduced model is built upon. A good approximation of the full order model can be achieved in the frequency range of interest for any given parameter a in the defined parameter space. The frequency response and model error evaluated for a = 0.0425 m are given in figure 5. Due to the fact, that the transfer function for each parameter a i has to be computed to populate the Loewner matrix, the computation time for the reduction phase t r = 3.966 s is significantly higher than the evaluation of one specific transfer function t f = 0.244 s. After computing the reduced order model, the evaluation time of the transfer function for any parameter in the parameter space a i reduces to t r = 0.065 s. The method presented in section 3.1 aims at reducing the initial computation time of the reduced Loewner model by obtaining the transfer functions from already reduced models.

Acoustic black hole
The previous section showed, that IRKA and the parametric Loewner framework can produce good quality reduced models of dynamic structures. Algorithm 4.1 will now be used to generate reduced models of a plate with acoustic black holes equipped with different damping materials. Contrary to the examples above, the ABH model is of larger scale, so the applicability of the reduction method in a design context is also evaluated.
To show the method's capabilities to generate models of ABHs, we equip an aluminum plate (E = 6.9e10 N/m 2 , ν = 0.22, ρ = 2650 kg/m 3 ), dimensions 0.6 × 0.5 × 0.003 m, with an ABH with diameter d = 0.2 m and remaining plate thickness in the middle of the ABH h = 0.001 m. In the ABH region, a damping material is applied on the opposite side of the plate and clamped by a layer of aluminum, as shown in figure 3. The model uses separate damping coefficients for plate and damping material, which is problematic for classical approaches like the modal superposition, where proportional damping is required [29]  The resulting parametric model has the dimensions k = 47 in the frequency space and q = 26 in the parameter space, which indicates, that snapshots for more parameters are required to interpolate the measurement data, rather than approximate it. Nevertheless, the reduced model is able to approximate the full model up to a frequency of 5000 rad/s. The resulting error norm ||H(ω)−Hr(ω)|| ||H(ω)|| = 6.377e−4 shows, that the reduced model is valid in the frequency range ω = 10 1 , . . . , 10 5 rad/s. The reduced parametric model is generated within t r = 14.7 s, the evaluation of the reduced model for a parameter takes t e = 0.23 s. It is clear, that the generation of the transfer function snapshots for 50 parameters results in a longer computation time than solving the system for one distinct parameter. But if the reduced model is to be evaluated for many parameters, for example in an optimization study, the initial overhead of generating the reduced model is easily compensated.

Concluding remarks
We presented an algorithm for the reduction of parametric second order dynamic systems. The individual MOR methods in the algorithm have been evaluated and show a good robustness and accuracy. The combined algorithm was tested on a large-scale model and also yielded acceptable results. A main advantage of the algorithm is the ability to reduce non-proportionally damped systems. Faster post-processing is also possible, as the reduced matrices can be used here.
Algorithm 4.1 already exploits the fact, that poles of parametric models evaluated at some parameters lie in a similar region and reuses the expansion points at convergence of one IRKA run as starting point for the next. A-priori optimal pole selection has the potential to speed up convergence also in the first IRKA run. The estimated modal density presented, for example, in [30] can be an indicator for the proper choice of the initial expansion points. Using appropriately preconditioned iterative solvers in the IRKA part of the algorithm for finding Krylov subspaces as suggested in [31] or the efficient generation of orthogonal subspaces using an Arnoldi procedure as shown in [32] would also help reducing computation times and reaching even faster convergence for the reduced models.
Regarding the Loewner framework part in algorithm 4.1, it will be crucial to find a rule on the spacing between parameters. The rank for the parameter space encoded in the Loewner matrix only has a valid meaning, if all features of the original models have been captured in the frequency response function snapshots. Without prior knowledge of the system response, valid parametric models can only be created using a rather dense parameter space. An adaptive greedy procedure to detect optimal sampling points in the parameter space could, for example, be used.