A Krylov Space-Based Finite Element Time Domain Method for Broadband Frequency Domain Solutions

A Krylov space based time domain method for wave propagation problems is introduced. The proposed method uses the Arnoldi algorithm to obtain broadband frequency domain solutions. This method is especially advantageous in cases where slow convergence is observed when using traditional time domain methods. The efficiency of the method is examined in several test cases to show its fast convergence in such problems.


Introduction
Wave equations can be generally solved using two categories of methods: time-domain and frequency-domain methods.Finite-difference time-domain (FDTD) methods have been widely adopted for solving different kinds of wave propagation problems.In these methods, the field is discretized into a series of uniform hexahedral volumes.The spatial derivatives are approximated with second-order central differences and the leapfrog method is used for temporal discretization [1].The interests in the time domain methods have been expanded over the past few decades for several reasons, such as the natural treatments of impulsive behaviors, nonlinear behaviors, and so on.
The solution to linear wave equations can be approximated by linear combinations of sinusoidal waves at various frequencies, which leads to a frequency-domain formulation for different time-harmonic sources.The time-domain solution is of interest in this research for its capability of obtaining frequency W. Y. Lin DOI: 10.4236/oja.2017.7400996 Open Journal of Acoustics domain solutions across a frequency range.Besides, in the time domain, the residual of the linear system is reduced several orders of magnitude, and requires much less search direction in a Krylov subspace [2].
The time domain methods, however, have been found to suffer from convergence problems for physical models that include locally resonant structures.
Such structures may be the result of large material mismatches, or complex geometries.The convergence problems have also been a concern in structural optimizations such as the designs of metamaterials [3] [4], where the geometry constantly changes during the design iterations and is not predictable.
The purpose of this paper is to introduce a Krylov space based time domain approach for frequency domain solutions.The proposed algorithm is suitable for physical models where slow convergence is found with traditional time domain methods.Related research can be found in [5] [6] [7], and so on.Different from existing work, the proposed algorithm is derived directly from discrete Fourier transform of time domain data, by which the corresponding frequency domain solutions of a wide range can be obtained with negligible computational costs once the projection of the time domain solution onto the Krylov subspace is obtained.No matrix operation with complex numbers is required in the algorithm.
Several numerical cases are examined to demonstrate the efficiency of the method, in which cases the spatial discretization is done with finite element methods.

Numerical Model
Wave propagations can be modeled by first order partial differential equations (PDEs) in general, which may be written as In acoustics, the wave equations can be written in a non-conservative form, known as the linearized Euler equations 1 0 where p is the acoustic pressure, u is the velocity vector, K e is the bulk modulus of compressibility of the material, and ρ e is the density.
The wave propagation problems in this paper are modeled using a stabilized finite element formulation: Streamlined Upwind/Petrov Galerkin (SUPG) method [8]  Finite element methods typically start with formulating the equations using a weighted residual method, which can be cast in the form of [ ] where φ is a weighting function.The temporal discretization of the governing equation is advanced by implementing a backward differentiation formula (BDF), and the fully discretized equation is solved with an implicit time marching approach (Newton's method) so that the time-step size can be determined by accuracy considerations instead of stability limitations as in explicit methods.

A Broadband Frequency Domain Method Based on Krylov Space Projections of Time Domain Solutions
While the time domain formulation can be used to obtain the frequency domain solution across a frequency range, the convergence of the frequency domain solution is not always as rapid as expected.Especially in physical problems where locally resonant structures exist, the solution procedure can easily take numerous time-marching steps before the errors are reduced to a tolerance.
To circumvent the problem of slow convergence in certain situations, the time domain method is reformulated with a Krylov space projection.Consider the fully discretized equations with BDF-1 scheme, the residual at time step j (t = jΔt) can be written as [ ] ( ) Since the equations are linear, the residual can be expressed as where and are the linearization matrices.Substitute this into the Newton's step, and use deduction, given an initial condition Q 0 , for any given time j, the time dependent solution can be described as which indicates that the time domain solution forms a Krylov basis with the iter- where nt is the number of time steps.It also indicates that, the damping of the numerical error is strongly related to the eigenvalues of the matrix [M]; the closer the eigenvalues are to 1.0, the slower the convergence would be.
Since the frequency domain solutions are of interest, the Fourier transform can be applied to the time domain solutions.At each frequency point f k , the discrete Fourier transform may be expressed as where ι 2 = −1 and Denote and substitute Equation (11) into Equation ( 13), the frequency domain solutions may now be written as The initial condition Q 0 can be projected into the eigen-space of matrix [M].
That is, it can be expressed as where nn is the number of nodes (degree of freedom), v i is the i th eigenvector of [M], and a i is the corresponding coefficient.Therefore the solution at time step j can be re-written as then the frequency domain solution becomes Switch the sum operators, and this can be expressed as W. Y. Lin DOI: 10.4236/oja.2017.7400999 Open Journal of Acoustics The set of the eigenvalues and eigenvectors can be approximated by the Arnoldi algorithm [9].The algorithm can be summarized with the following: Algorithm.A Krylov space-based time domain method for broadband frequency domain solutions 1) Use Q 0 as the initial condition for time-domain method to get matrices [A] and [B] as in Equations (7-8).
2) Set i = 1, and set and start the Arnoldi iteration.Store the q i vectors in the [Q] matrix, and h j,i items in the [H] matrix.

5) Calculate
, ) If j = i, stop the inner loop; otherwise set j = j + 1 and go to step 5.

11) Calculate the eigenvalues [Λ] and eigenvectors [V H ] of the Hessenberg matrix [H] obtained by the Arnoldi iteration, and the eigenvectors [V] of the [M]
matrix is given by  13) and add it to the primitive frequency domain solution.
For physical problems which require a lot of iterations before the frequency domain solutions reach convergence, this algorithm is recommended.In fact, the eigenvalues in these problems are found to be close to 1.0, which can be effectively approximated by the Arnoldi algorithm since these eigenvalues are typically extremes in the eigensystems.The advantage of using this algorithm is that, the number of time steps nt can be chosen to be an arbitrary number, so that the number of iterations is reduced from nt to niter as is needed for the Arnoldi iteration.

Numerical Examples
An acoustic wave propagation case with artificial material properties is designed to examine the efficiency of the proposed algorithm.The geometry is shown in Figure 1, where the background material (material 1) is shown in gray color, and the inclusion (material 2) is shown in black color.This geometry is a simple locally resonant structure.The unit cell is a square with edge length a = 1 cm.
The frequency range of interest is from 0.5 kHz to 3 kHz as reported in a design of acoustic metamaterials [4].Three cases with different material properties are considered, and the relative material properties of the inclusion are shown in Table 1.
Figure 1.The geometry of the acoustic wave propagation case.
Table 1.The relative material properties of the inclusion.
Cases and material properties  The simulation is configured as initiating a Gaussian pulse in the left end, and sensors are used to collect information for use in the frequency domain solutions.
In addition, periodic boundary conditions are applied to the upper and lower parts of the domain, and absorbing boundary conditions [10] are applied to the left and right ends.
The frequency domain solutions for a sensor located in the center of the computational domain are used as the reference.Since the convergence of the solutions is of primary concern, only the frequency domain solutions after the completion of the Gaussian pulse are considered.It is shown in Figure 2 that the time domain solutions for cases 1 and 2 take a long time before settling down.It is thus to be expected that the convergences of the frequency solution will be slow, and the errors to the converged numerical solutions are plotted in Figure 3.
In practice, the numbers of time steps required for cases 1 and 2 to reach the error level of 10 −10 are 441,769 and 26,280, respectively.In comparison, in both cases the solutions converge to the level of 10 −10 within 100 Arnoldi iterations with the proposed time domain-based broadband frequency domain algorithm, as shown in Figure 4.However, it is also noted that the solution may diverge due to the fact that the eigenvalues approximated by the Arnoldi algorithm can be bigger than 1, and their exponentials may be exceedingly large.It is thus necessary to monitor the approximated eigenvalues in some cases to ensure their values are all below 1, and if one or more of them are larger than 1, more Arnoldi iterations need to be taken.promptly with time marching method, and the resonant frequency is not obvious, due to the fact that most of the vibrations are damped after the excitation is completed.In addition, the Arnoldi iteration is also able to reduce the error level for this case.

Conclusion
A Krylov space-based time domain method is introduced for broadband frequency domain solutions.This method is useful for dealing with physical models for which slow convergences are observed with traditional time marching methods.The efficiency of the method is examined in several test cases to show its fast convergence in such problems.An acoustic wave propagation problem is modeled with a stabilized finite element method.In the test cases, the proposed method uses less than 100 iterations before an error level is reached.This decrease is orders of magnitude less than the time-marching approach.
a. K e is the bulk modulus of compressibility of the material, and ρ e is the density.b.The superscripts represent the acoustic materials 1 (background) and 2 (inclusion).

Figure 5 Figure 2 .
Figure 2. The time domain solutions at a sensor located in the center of the geometry.

Figure 3 .
Figure 3.The convergence of frequency domain solutions with the time-marching approach.

Figure 4 .
Figure 4.The convergence of frequency domain solutions with the Krylov space-based time domain algorithm.

Figure 5 .
Figure 5.The Frequency domain solutions at a sensor located in the center of the geometry.