3D Helmholtz wave equation by fictitious domain method

A fictitious domain method for solving acoustic scattering problems is considered. A nonsymmetric complex-valued augmentation of the stiffness matrix is proposed. The augmentation block corresponds to the same wave operator inside the scatterer with the absorbing boundary conditions posed on its boundary. An iterative method in a subspace of /i-harmonic functions is used. The arithmetic complexity of each step of the method (except an initialization step) depends on the number of mesh nodes on the scatterer boundary. The numerous experiments justify a good convergence rate of the method. The goal of this article is to present an effective numerical method for solving the exterior boundary value problem of acoustic scattering. The solution of the problem satisfies the Helmholtz wave equation outside the scatterer Ο as well as an outgoing radiation condition at infinity. The unboundedness of the exterior problem is treated by introducing an artificial sphere of sufficiently large radius R and an approximate local radiation boundary condition on this sphere. The local radiation conditions lead to wellposed problems [3,7,11] and have a simple form. Unfortunately, the radiation conditions of only the first and second orders are suitable for a finite element approximation. The first-order radiation boundary condition results in an error O(R~) as R increases [11,12]. Thus, the benefits of approximating of local radiation conditions lead to difficulties in solving large algebraic problems. One of the possibilities to overcome these difficulties is to use iterative methods in subspaces of small dimension [13,14]. We show in Section 3.3 that the subspace in our method is associated with the scatterer boundary dO, its dimension is proportional to the number NQO of mesh nodes on dO and consequently independent of R. Moreover, the arithmetic complexity of each step of the iterative procedure is also independent of R (except for an initialization step in the subspace entering) and is of order N^ In Nd0. The methodology which we follow in this paper was originally proposed within algebraic framework in [1,18] for solving the Poisson equation on rectangular grids. It has many common features with other methodologies based on the fictitious domain or domain embedding procedures in particular with matrix capacitance methods [8,20] and fictitious domain methods with Lagrange multipliers [2]. The first attempts to apply fictitious domain methods for solving the Helmholtz equation were done in [4,15]. The basic idea of our algorithm is to modify the standard non-symmetric version of the fictitious domain method by involving absorbing boundary conditions in constructing an augmented block in the stiffness matrix. In other words, to overcome difficulties related to indefiniteness of the Helmholtz operator, we add to sesquilinear form a complexvalued term associated with the scatterer boundary dO. In Section 3.5 we demonstrate the influence of absorbing boundary conditions in the case of a model spherically symmetrical problem. In the asymptotical case, when the incident wavelength tends to infinity, the additional term disappears reducing our method to that used for solving problems with the Laplace operator (see [19] and the references therein). 'Institute of Numerical Mathematics, Russian Academy of Sciences, Moscow 117951, GSP-1, Russia 372 Yu.A. Kiiznetsov and K.N. Lipnikov From the physical viewpoint our method can be treated as follows. It is well-known that non-symmetric augmentation of the stiffness matrix would be optimal if we could construct on dO an absorbing boundary condition transparent for all waves from the sources placed inside the scatterer. Therefore the idea of our algorithm is to approximate an unknown perfectly absorbing condition by a simple one. In numerical experiments we used the first-order local radiation condition. The transparency properties of this condition are aggravated when we are moving from convex to concave boundaries. The results of the numerical experiments presented in Section 4 support this fact. The convergence rate of the iterative method is changed significantly when we are moving from ellipsoids (scatterers with the concave boundary) through star-shaped obstacles to scatterers with convex boundaries. Nevertheless, in the last case we also observe good convergence properties of our method. 1. 3D ACOUSTIC SCATTERING PROBLEM Let us consider a bounded simply-connected domain Ο £ Re with a piecewise smooth boundary dO. In the unbounded domain Ω = Re \ 0 we formulate the following problem: Au, + kus = 0 in Ω u3 = g on dO /^ ^ du3 . . — iku3 = o(r~) uniformly as r -» oo dr where k is a positive constant, g is a smooth complex-valued function, i is the imaginary unit, and r = \x\. The problem formulated describes the scattering of time-harmonic acoustic waves propagating in a homogeneous medium Ω. The function us is a complex amplitude of the scattered wave. The wave number k and the wavelength L satisfy the relation k = 2π/1. The Dirichlet boundary condition on dO is physically defined by an incident wave Ui, i.e. g = -ui on dO. In order to treat the unboundedness of the domain Ω, we approximate the problem (1.1) by introducing a sphere TR of sufficiently large radius R as well as a radiation boundary condition on this sphere. Let us denote by ΩΗ the domain bounded by T R and dO and replace the problem (1.1) by Δι/f + ku* = 0 in ΩΛ uf = g on dO (1.2) = 0 on TR where the operator Mk corresponds to the chosen radiation boundary condition. In the scattering problems such boundary conditions are often spoken of as non-reflecting or absorbing ones [10]. On a spherical boundary we use the first-order radiation condition given by The approximate problem is well-posed [3,11,12], and