Time-harmonic solution for acousto-elastic interaction with controllability and spectral elements

The classical way of solving the time-harmonic linear acousto-elastic wave problem is to discretize the equations with ﬁnite elements or ﬁnite diﬀerences. This approach leads to large-scale indeﬁnite complex-valued linear systems. For this kind of systems, it is diﬃcult to construct eﬃcient iterative solution methods. That is why we use an alternative approach and solve the time-harmonic problem by controlling the solution of the corresponding time dependent wave equation. In this paper, we use an unsymmetric formulation, where ﬂuid-structure interaction is modeled as a coupling between pressure and displacement. The coupled problem is discretized in space domain with spectral elements and in time domain with central ﬁnite diﬀerences. After discretization, exact controllability problem is reformulated as a least squares problem, which is solved by the conjugate gradient method.


Introduction
Acoustic waves are small oscillations of pressure, which are associated with local motions of particles in fluid domain Ω f .The linear theory of elasticity models mechanical properties in structure Ω s assuming small deformations.
Acousto-elastic interaction between these two media constitutes a coupled problem.Several phenomena, such as seismic waves in the earth and ultrasonic waves used to detect flaws in materials, can be described by an acoustoelastic model.Two approaches, in which the displacement is solved in the elastic structure, predominate in modeling the interaction between acoustic and elastic waves.Expressing the acoustic wave equation by the velocity potential results in a symmetric system of equations (see, e.g., [1,2,3,4]), while using the pressure in the fluid domain leads to an unsymmetric formulation (see, e.g., [5,6,7,8]).
In this paper, we present the acousto-elastic interaction between pressure and displacement, and thereby consentrate on the unsymmetric approach.We formulate the time-harmonic acousto-elastic interaction as an exact controllability problem [9] via the corresponding time dependent system.The time dependent problem is discretized in space domain with a higher-order spectral element method (SEM) and in time domain with second-order central finite differences.The combination of these discretization methods is well known with wave equations (see, e.g., [10]).The methods related to spectral elements are studied in context of the time dependent acousto-elastic problem with second-order time-stepping schemes for instance in references [11,12,6].
After discretization, we solve the control problem by a conjugate gradient (CG) algorithm which is related to the one developed in [13] for the acoustic wave equation.If an unpreconditioned CG algorithm is used, the number of iterations grows rapidly with the order of spectral element [14].That is why we use a modification of Kickinger's [15] algebraic multigrid (AMG), introduced in [16], for preconditioning the conjugate gradient algorithm.
The rest of this paper is organized as follows.First, the mathematical model is presented in Section 2.Then, we discretize the coupled problem in space domain with spectral elements in Section 3.For time discretization we use central finite differences in Section 4. In Section 5, we present the control problem and the preconditioned conjugate gradient algorithm.Finally, we show some numerical experiments in Section 6.

Mathematical model
We consider the use of a control algorithm to solve the time-harmonic acoustoelastic problem in the domain Ω ⊂ R 2 , which is divided into the solid part Ω s and the fluid part Ω f by the interface Γ i (see Figure 1).Instead of solving directly the time-harmonic equation, we return to the corresponding time dependent equation (see, e.g., [17,10]) and look for time-periodic solution.The convergence is accelerated with a control technique by representing the origi-nal time-harmonic equation as an exact controllability problem [18,19] for the time dependent wave equation where f , y ext , f , and g ext are the source terms.Length of the time interval is marked as T , p f denotes the pressure, and u s = (u s1 , u s2 ) T is the displacement field depending on the spatial variable x = (x 1 , x 2 ) T ∈ R 2 .Coefficients ρ f (x) and ρ s (x) represent the densities of media in domains Ω f and Ω s , respectively, and c(x) is the speed of sound in fluid domain.The stress tensor is expressed as σ(u s ) = ρ s (x) (c p (x) 2 − 2c s (x) 2 ) (∇ • u s )I + 2ρ s (x)c s (x) 2 (u s ) with the speed of the pressure wave c p (x), the speed of the shear wave c s (x), the identity matrix I, and the linearized strain tensor = 1 2 ∇u s + (∇u s ) T .The outward normal vectors to domains Ω f and Ω s are marked as n f = (n f 1 , n f 2 ) T and n s = (n s1 , n s2 ) T .
The fluid domain is bounded by Γ f = Γ 0f Γ ef Γ i , and Γ s = Γ 0s Γ es Γ i constitutes the boundary for the solid domain Ω s .The boundaries Γ 0f and Γ 0s are assumed to be rigid, whereas on the artificial boundaries Γ ef and Γ es we impose the conventional first-order absorbing boundary conditions [20,21], where B is a symmetric positive definite 2 × 2-matrix defined by In addition to the system (1)-( 8), we take into account the initial conditions e = (e 0 , e 1 ) T such that e 0 = (e 0f , e 0s ) T and e 1 = (e 1f , e 1s ) T , and For existence and uniqueness of the solution for the problem ( 1)-( 10), we refer to [22], and for the corresponding time-harmonic problem to [23].

Spatial discretization
For space discretization, we use the spectral elements method (see, e.g., [24,17,6,10]), which is based on the weak formulation of the system (1)- (8).That is why we introduce the function spaces V and V by Multiplication of the equation ( 1) with any test function v in the space V , and (5) with any test function v in the space V, and use of the Green's formula result in the following weak formulation: for any (v, v) ∈ (V × V) and t ∈ [0, T ] with The computational domain Ω is divided into N e quadrilateral elements Ω i , i = 1, . . ., N e such that Ω = Ne i=1 Ω i .For the discrete formulation, we define the reference element Ω ref = [0, 1] 2 and affine mappings Then, the spectral element discretization is obtained from the weak formulation by restricting the problem defined in the infinite dimensional spaces V and V into finite dimensional subspaces V r h and V r h .These discrete subspaces are given by V a pq ∈ R} is the set of Lagrange interpolation polynomials of order r in R 2 .In one space dimension, the nodes of the basis functions are placed at the rth order Gauss-Lobatto (GL) discretization points that are the zeroes of x 1 (1 − x 1 )L r (2x 1 − 1), x 1 ∈ [0, 1], where L r is the derivative of the rth degree Legendre polynomial L r .The GL points in R 2 are given by the tensor product of the one-dimensional GL points.The integrals in the weak form of the equation are evaluated with the corresponding Gauss-Lobatto quadrature formulas.This leads to the semi-discretized coupled problem where U is the global vector containing the values of the displacement u s (x, t) and the pressure p f (x, t) at time t at the Gauss-Lobatto points of the spectral element mesh.The entries of the matrices M, S, and K, and the right hand side vector F, are given by the formulas where matrix and vector blocks M s , S s , K s , and F s represent the elastic waves, M f , S f , K f , and F f correspond to the fluid domain, and matrices A fs and A sf arise from the coupling between acoustic and elastic wave equations.Since M s and M f are diagonal matrices, the inverse of the lower diagonal block matrix M is easily computed, and explicit time stepping with central finite differences requires only matrix-vector multiplications.

Time discretization
The time discretization of the semi-discrete equation ( 19) is performed with the central finite differences (see, e.g., [10]).This method is second-order accurate with respect to the time step ∆t and leads to an explicit time-stepping scheme.Both properties are essential for computational efficiency.
The time interval [0, T ] is divided into N time steps, each of size ∆t = T /N .After replacing the time derivatives in the semidiscretized form (19) by the approximations and taking into account the initial conditions ( 9)- (10), we obtain the fully discrete state equation where contains the vectors U i , F i is the vector F at t = i∆t, and e 0 = (e 0s , e 0f ) T and e 1 = (e 1s , e 1f ) T are the initial conditions.The matrix blocks B, C, and D are given by the formulas while I is the identity matrix.

Conjugate gradient algorithm
Essentially, the solution procedure of the exact controllability problem is similar to those presented for the Helmholtz equation in [19,25] and for the Navier equation in [26].After discretization, the exact controllability problem is reformulated as a least-squares optimization problem 1 2 where we use a short notation L for the block-diagonal matrix containing the non-coupling terms of the matrices K and M such that The minimization problem ( 26) is solved with a preconditioned conjugate gradient algorithm.Each conjugate gradient iteration requires computation of the gradient of the discretized least-squares functional, solution of a linear system with the block-diagonal preconditioner L, and some matrix-vector operations.Computation of the gradient is an essential point of the method.By following the adjoint state technique (see, e.g., [25]), we obtain the gradient where Z 0 and Z 1 are solutions of the adjoint state equation at time t = 0 and t = ∆t, respectively.
For preconditioning the algorithm, we use the algebraic multigrid (AMG) method [16] (see also [27,19]).As a smoother for the AMG we use the successive over relaxation (SOR) method with relaxation factor equal to 1.2.One iteration of the SOR is used as pre-and post-smoothing.Additionally, in the beginning of every multigrid iteration, four iterations of the SOR are used to smooth the solution initially.So called W-cycle [28] is utilized as a multigrid iteration until the residual norm of the solution is smaller than 10 −6 .

Numerical experiments
In this section, we show some numerical results in order to validate the method discussed in previous sections.The material parameters in fluid domain are ρ f (x) = 1.0 and c(x) = 1.0.In solid domain, we use the values c p (x) = 6.20, c s (x) = 3.12, and ρ s (x) = 2.7.Angular frequency ω is the same for both media, and throughout the tests we set the propagation direction (-1,0) by the vector ω = (ω 1 , ω 2 ) = (−1, 0)ω.For test problems, we use the right hand side functions Algorithm 1 Preconditioned CG algorithm Initialize e = 0. Solve the state equation (22) with initial condition e. Solve the corresponding adjoint state equation.Compute the gradient g(e, Û(e)).
Solve linear system with the preconditioner Lw = −g.Set ς 0 = −(w, g) and ς = ς 0 .repeat Solve the state equation ( 22) with initial condition w and F = 0. Solve the corresponding adjoint state equation.

Table 1
The ratio between time step and mesh stepsize for different element orders.
Element order r 1 2 3 ∆t/h 0.1250 0.0227 0.0101 sin(ω • x) cos(ωt), , where According to the results reported in [25,26], we choose the length of time step, for a certain element order r, to reduce the temporal error to a lower level than the spatial error.These values are presented in Table 1.For each element order r, we construct square element meshes, which are matching on the interface Γ i .Numerical experiments are carried out on an AMD Opteron 885 processor at 2.6 GHz, and iterations are continued until the stopping criterion ε = 10 −4 is reached.

Accuracy
The analytical solution of the first test problem is known to be    2 shows how the accuracy improves when element order grows with ω = 8π and h = 1/20.The number of iterations required to attain the stopping criterion was sufficiently large (see Figure 3).

Scattering
This scattering test is solved in domain Ω, where the absorbing boundary coincides with the perimeter of the rectangle [0, 5] × [0, 3.75].In the center of the computational domain Ω, we have two rigid non-convex semi-open cavities.These reflectors are located at perpendicular distance of 1 from the absorbing boundary.Thus, the lower left corner of the left side cavity is at the point (1, 1) and the lower right corner of the right side cavity is at the point (4, 1).Internal width and height of each cavity are 0.75 and 1.25, respectively.Thickness of  4-5, where displacements are presented as vector fields and pressures are presented as contours.

Conclusions
We considered the spectral element solution of time-harmonic acousto-elastic interaction problems by the exact controllability method.Simulation results show that the number of iterations required to attain the stopping criterion is large but independent of the element order.In future work, the central finite difference time discretization will be substituted by higher-order time discretization methods to reduce the number of iterations and improve the accuracy.Also symmetric fluid-structure interaction formulations, with coupling between velocity potential and displacement, will be discussed.

Figure 1 .
Figure 1.The domain Ω is divided into the solid part Ω s and the fluid part Ω f .