FINITE ELEMENT SOLUTION OF SCATTERING
IN COUPLED FLUIDSOLID SYSTEMS
by
Mirela O. Popa
B.S., Harvey Mudd College, Claremont CA, 1995
M.S., University of Colorado at Denver, 1997
A thesis submitted to the
University of Colorado at Denver
in partial fulfillment
of the requirements for the degree of
Doctor of Philosophy
Applied Mathematics
2002
This thesis for the Doctor of Philosophy
degree by
Mirela 0. Popa
has been approved
by
Jan Mandel v
12, 2/&02
Date
Popa, Mirela O. (Ph.D., Applied Mathematics)
Finite element solution of scattering in coupled fluidsolid systems
Thesis directed by Professor Jan Mandel
ABSTRACT
In this thesis we investigate the mathematical theory of wave scatter
ing by an obstacle. The obstacle considered is a bounded elastic body in a fluid
domain. We analyse finite element methods and show existence and uniqueness
of the solution for the coupled fluidsolid interaction problem in more than one
dimension. We study the stability and regularity properties of acoustic wave
scattering by introducing interpolation of spaces and scaled norms. This type
of analysis, to the authors knowledge, has not been investigated. Then we con
sider a multigrid method for the coupled fluidsolid interaction model problem
in higher dimension.
We use the Garding Inequality to obtain coercivity, then we use the
Fredholm Alternative to analyze spectral properties and show uniqueness and
existence of solutions. We need only weak regularity assumptions, and give
a rigorous treatment of scale of spaces with constants independent of wave
number k. A new approach is used to show stability of the coupled problem
by using intermediate spaces and norms.
m
Finally, we present a multigrid algorithm to solve a coupled solid
fluid interface problem. To the authors knowledge, multigrid methods for the
coupled acousticelastic problem have not been investigated. In this thesis, we
formulate such a method and present numerical experiments from a prototype
implementation in MATLAB.
This abstract accurately represents the content of the candidates thesis. I
recommend its publication.
Signed
Jan Mandel
IV
DEDICATION
To my family.
ACKNOWLEDGMENTS
I would like to thank my advisor, Prof. Jan Mandel, for his support,
guidance, and generosity.
This research was supported by the National Science Foundation un
der grants ECS9725504 and DMS0074278, and by the Office of Naval Research
under grant N000149510663.
CONTENTS
Figures ............................................................ ix
Tables.............................................................. xi
1. Introduction......................................................... 1
2. Existing Methods..................................................... 5
3. Theoretical Preliminaries........................................... 15
3.1 Sobolev Spaces.................................................... 15
3.2 Finite Element Approximation...................................... 16
3.3 Garding Inequality................................................ 19
3.4 Generalized Korn Inequality....................................... 20
3.5 Fredholm Alternative.............................................. 20
3.6 Trace Theorem..................................................... 21
3.7 Riesz Representation Theorem...................................... 22
3.8 LaxMilgram Theorem............................................... 22
3.9 Hilbert Interpolation Spaces...................................... 22
4. Statement of Coupled Problem........................................ 24
4.1 Derivation of the Coupled Problem................................. 24
4.1.1 Elastic Waves................................................... 24
4.1.2 Acoustic Waves.................................................. 28
4.1.3 Boundary Conditions............................................. 31
4.1.4 SolidFluid Interface Conditions............................... 32
vii
4.2 Variational Form of Coupled Problem............................... 34
4.3 Weak Formulation ................................................. 36
4.3.1 Hilbert scale .................................... 40
5. Analysis............................................................ 42
5.1 Garding Inequality for the Coupled Problem....................... 42
5.2 Existence of Solution............................................. 46
5.3 Intermediate Spaces .............................................. 56
5.3.1 Intermediate norms................................. 57
5.3.2 Regularity of Solution for Coupled Problem...................... 63
5.4 Discretization and Error Bound.................................... 65
6. Multigrid Method.................................................... 71
6.1 Multigrid for the Coupled Problem................................. 72
6.1.1 GMRES.......................................,................ 76
6.1.2 BICGSTAB .................................... 76
6.1.3 GaussSeidel .................................... 77
7. Numerical Results.................................................. 78
7.1 Numerical Verification of the Discretization ..................... 78
7.2 Computational Results with Multigrid.............................. 89
7.3 Conclusion and Future Work....................................... 114
References............................................................ 118
vm
FIGURES
Figure
4.1 Problem setup....................................................... 34
7.1 Configuration of sample points for Table 7.1. Sample points
1,2,3,4 are fluid pressure values, and sample points 5 and 6 are
ux and uy displacements, respectively............................. 80
7.2 Log log plot of difference of solution values xh at sample points
for mesh size 10 x 10 to 320 x 320 and extrapolated exact solution
x*, k=5............................................................. 82
7.3 Solution for a 40 x 40 mesh, k = 5, righthand side modified for
the Dirichlet boundary condition p = p0 on ..................... 83
7.4 Exact solution for a 64 x 64 mesh, k = 15, righthand side is
modified for the Dirichlet boundary condition p = p0 on Fd . . 84
7.5 Three multigrid iterations, 20 smoothing steps, smoother GM
RES, 2 levels to solve a 64 x 64 mesh, k = 15, righthand side
is modified for the Dirichlet boundary condition p = p0 on . 85
7.6 Contour of an obstacle (0.4 m) and (0.2 m) in the x and the
y direction, respectively. The gap is on the y axis. The size of
the gap is 0.5 or 50% in the x direction and 0.4 or 40% in the y
direction....................................................... 87
IX
7.7 Contour of an obstacle (0.2 m) and (0.4 m) in the x and the y
direction, respectively. The gap is on x axis. The size of the gap
is 0.4 or 40% in the x and 0.5 or 50% in the y direction........ 87
7.8 Solution for a 64 x 64 mesh with ushaped obstacle (0.2 m) by
(0.4 m). The obstacle has a gap of size 0.4 by 0.5 on the y axis.
The wave number is k = 15, and the righthand side is modified
for the Dirichlet boundary condition p = p0 on Td............... 88
7.9 Solution for a 64 x 64 mesh ushaped obstacle (0.2 m) by (0.4
m). The obstacle has a gap of 0.4 by 0.5 on the x axis. The
wave number is k = 15, and the righthand side is modified for
the Dirichlet boundary condition p = p0 on T^................... 90
7.10 Residual reduction as a function of adding coarse meshes, de
creasing mesh size h while keeping k3h2 constant................... 98
7.11 Residual reduction as a function of adding coarse meshes, de
creasing mesh size h while keeping k3h2 constant................... 99
7.12 Residual reduction as a function of adding smoothing steps, de
creasing mesh size h while keeping k3h2 constant.................. 100
7.13 Residual reduction as a function of adding smoothing steps, de
creasing mesh size h while keeping k3h2 constant.................. 101
7.14 Relative residual varies as mesh size h is decreased for the case
in which k3h2 is constant and in the case of resonance.......... 112
7.15 Solution for a 64 x 64 mesh, k = Anr, and righthand side mod
ified for the Dirichlet boundary condition p = po on T^. A gap
of half wavelength in x and y direction......................... 116
x
TABLES
Table
7.1 Values of the solution at the 6 sample points as explained in
(Fig.7.1), mesh size h is halved at each run and wave number is
kept constant k 5............................................. 79
7.2 Exact solution is Â£* and solution at mesh size his Xh........... 81
7.3 h = step size, k = wave number, sm number of smoothing
steps, lv = number of levels, mth = iterative method used for the
smoothing algorithm, it = number of multigrid iterations, res red
= residual reduction, rel rs = relative residual, res fl = residual
in the fluid part, and res el = residual in the elastic part. BC
is BICGSTAB, GM is GMRES, GD is GMRES preconditioned
by D_1 (as explained in (7.5)), GA is GaussSeidel .......... 92
7.4 h step size, k = wave number, sm = number of smoothing
steps, lv = number of levels, mth = iterative method used for the
smoothing algorithm, it = number of multigrid iterations, res red
= residual reduction, rel rs = relative residual, res fl = residual
in the fluid part, and res el = residual in the elastic part. BC
is BICGSTAB, GM is GMRES, GD is GMRES preconditioned
by D~l (as explained in (7.5)), GA is GaussSeidel .......... 93
XI
7.5 Iteration counts for multigrid Vcycle to achieve a relative resid
ual of order le6, for smoothers GMS=GMRES; BGS=BICG
STAB; GMD=GMRES preconditioned by the inverse of the
lower triangular part of A.................................... 95
7.6 h = step size, k = wave number, sm = number of smoothing
steps, lv = number of levels, mth = iterative method used for the
smoothing algorithm, it = number of multigrid iterations, res red
= residual reduction, rel rs = relative residual, res fl = residual
in the fluid part, and res el = residual in the elastic part. BC
is BICGSTAB, GM is GMRES, GD is GMRES preconditioned
by D~x (as explained in (7.5))................................ 96
7.7 h = step size, k = wave number, sm = number of smoothing
steps, lv = number of levels, it = number of multigrid iterations,
res red = residual reduction, rel rs = relative residual, res fl =
residual in the fluid part, and res el = residual in the elastic part. 103
7.8 h step size, k wave number, sm = number of smoothing
steps, lv = number of levels, mth = iterative method used for the
smoothing algorithm, it = number of multigrid iterations, res red
= residual reduction, rel rs = relative residual, res fl = residual
in the fluid part, and res el = residual in the elastic part. BC is
BICGSTAB, GM is GMRES, GT is GMRES preconditioned by
the inverse of the lower triangular part of A, GA=GaussSeidel. 104
Xll
7.9 h = the step size, k = the wave number, sm = the number
of smoothing steps, lv = the number of levels, sm =smoother,
mth = the iterative method used for the smoothing algorithm,
it = the number of multigrid iterations, res red = the residual
reduction, rl res = the relative residual, rs fl = the residual in
the fluid part, and rs el = the residual in the elastic part. BC
is BICGSTAB, GT is GMRES preconditioned by the inverse of
the lower triangular part of A................................. 106
7.10 h = the step size, k = the wave number, sm = the number of
smoothing steps, lv = the number of levels, mth = the iterative
method used for the smoothing algorithm, it = the number of
multigrid iterations, res red = the residual reduction, rl res =
the relative residual, rs fl = the residual in the fluid part, and
rs el = the residual in the elastic part. BC is BICGSTAB, GT
is GMRES preconditioned by the inverse of the lower triangular
part of A ..................................................... 107
7.11 h = the step size, x is size of obstacle in x direction, y size of
obstacle in y direction, gap x is size of the gap on the x axis as
a percentage, gap y is size of gap on the y axis, k = the wave
number, rs red = the residual reduction, rl rs = the relative
residual, rs fl = the residual in the fluid part, and rs el = the
residual in the elastic part........................................ 108
xm
7.12 h = the step size, x is size of obstacle in x direction, y size of
obstacle in y direction, gp x is size of the gap on the x axis as
a percentage, gp y is size of gap on the y axis, k = the wave
number, rs red = the residual reduction, rl rs = the relative
residual, rs fl = the residual in the fluid part, and rs el = the
residual in the elastic part..................................... 109
7.13 h = the step size, k = the wave number, rel res = the relative
residual, rel res fl = the relative residual in the fluid part, and
rel res el = the relative residual in the elastic part........... 110
7.14 GMRES preconditioned by ILU h = the step size, k = the wave
number, rel res = the relative residual, res fl = the residual in
the fluid part, and res el = the residual in the elastic part. ... Ill
7.15 h = the step size, k = the wave number, conv fctr = convergence
factor or the residual reduction, rel res = the relative residual,
res fl = the residual in the fluid part, and res el = the residual
in the elastic part............................................... 112
7.16 h = the step size, k = the wave number, conv fctr = convergence
factor or the residual reduction, rel res = the relative residual,
res fl = the residual in the fluid part, and res el = the residual
in the elastic part.............................................. 113
xiv
7.17 Residual Reduction per unit of work for one Multigrid cycle, h
= the step size, x is size of obstacle in x direction, y size of
obstacle in y direction, gp x is size of the gap on the x axis as
a percentage, gp y is size of gap on the y axis, k = the wave
number, rs red = the residual reduction, rl rs = the relative
residual, rs fl = the residual in the fluid part, and rs el = the
residual in the elastic part......................................... 115
xv
1. Introduction
Time harmonic acoustics of coupled fluidsolid systems in fluid pres
sure and solid displacement formulation has many wide ranging applications.
For this reason it has been extensively studied in recent years in an effort
to predict the dynamic response of elastic objects immersed in fluids. Wave
propagation in fluid, assuming timeharmonic behavior, is described by the
Helmholtz equation,
Ap + k2p = 0, (1.1)
where the the critical parameter is the wave number k. In elastic media, waves
propagate in the form of oscillations of the stress field, and are described by
the elastodynamic Helmholtz equation,
V r + uj2pu = 0, (1.2)
where u is displacement, p is density of the elastic medium, and r is the stress
tensor. Wave propagation in a composite medium, consisting of a fluid part
and an elastic part, is described by (1.1) and (1.2), in the acoustic and in
the solid part, respectively, together with boundary conditions imposed on the
fluid solid interface.
Important applications of elastic wave scattering are found in geo
physical exploration, seismology, the response of structures to seismic waves,
1
acoustic response of elastic objects immersed in fluids, response of aircraft parts
to dynamic loads, and applications of ultrasound to biological systems for di
agnostic of therapeutic purposes [64, 69]. A significant growth in the literature
employing scattering problems can be seen in recent years; [14, 48, 50, 34, 60,
59, 21]. The question of existence and uniqueness of solutions of Helmholtz
problems was addressed by the end of the 1950s; [46, 13, 58].
In this thesis, we analyze a timeharmonic solution of coupled fluid
solid interaction model problem in more than one dimension. We analyse finite
element methods, show existence and uniqueness of the solution, and study its
stability and regularity properties by introducing interpolation of spaces and
scaled norms. We give a rigorous treatment of scale of spaces with constants
independent of the wave number k. We present a multigrid method for the
solution of linear systems arising from the finite element discretization of the
coupled fluidsolid system. The obstacle considered is a bounded elastic body
embedded in fluid.
It is known, that numerical solutions to the Helmholtz equation de
teriorate for increasing wave numbers k, [37, 36], which is called the pollution
effect. The effect of pollution is that the wave number of the finite element
solution is different from the wave number of the exact solution, and pollution
means that the local error has a global effect, [4, 20]. In one dimension the
pollution effect has been extensively studied, and different approaches have
been proposed that lead to solutions that do not suffer the pollution effect,
2
[3, 4]. In two and three dimensions it has been shown that pollution cannot be
avoided, [4]. In this thesis we will not be concerned with the pollution effect
or dispersion, and we will keep in mind that for large wave numbers, numer
ical pollution in the error dominates the error of interpolation. Completely
different methods are required to address the pollution effect.
Iterative methods consisting of alternating solution in the fluid and
the solid region are known [14, 50]. Numerical solution of elliptic partial differ
ential equations, in two or three dimensions, is a typical application for iterative
solvers based on the multilevel paradigm. In contrast to other methods, multi
grid does not depend on the separability of the equations; thus using multigrid
methods for solving the Helmholtz equation in modeling acoustic scattering is
classical [25, 31, 63, 62]; for more recent developments, see [22, 44, 45, 62].
The remainder of the thesis is organized as follows. In Chapter 2,
we give an overview of some existing methods and approaches to solving the
problem of scattering by an elastic obstacle in an acoustic fluid. In particular,
we describe the analysis of a fluidsolid interaction problem in one dimension
studied by Makridakis et. al. [48]. In Chapter 3, we present some theoretical
preliminaries, and describe notation we use later in this thesis. In Chapter 4,
we present the relations of linear wave physics, derive the differential equations
with an emphasis on the boundary conditions that describe the solidfluid inter
action, and then we define the spaces necessary to derive the weak formulation
of the coupled problem. Chapter 5 is concerned with analyzing the existence
3
and uniqueness of the solution of scattering by an elastic obstacle in an acoustic
fluid in a bounded region. We use the Garding Inequality to obtain coercivity,
then we use the Fredholm Alternative to analyze the spectral properties and
show existence of solution. We show that if a solution exists, then the solution
is unique. By using intermediate spaces, we show stability of the solution.
In Chapter 6, we present a multigrid algorithm to solve a coupled solidfluid
interface problem. In Chapter 7, we present numerical experiments from a
prototype implementation in MATLAB.
4
2. Existing Methods
The scattering and propagation of waves is a classical area, which
has engaged the interest of several generations of mathematicians, physicists,
and engineers. Extensive research in the area of acoustics, electromagnetics,
and elastodynamics has resulted in many mathematical and computational
techniques [5, 8, 47, 64, 68, 69, 70].
Bielak, MacCamy, and Zeng [8] studied acoustic scattering in an un
bounded domain, representing the solution in the unbounded domain by a
combination of potentials. Demkowicz [15, 16] considered elastic scattering in
unbounded domains for spherical shells submerged in fluid, and showed that
an infsup condition holds for the reduced problem on the boundary of the
scatterer. In [21], Djellouli, Farhat, Tezaur, and Macedo apply a finite element
method coupled with the BaylissTurkellike non reflecting boundary condition
to solve direct acoustic problems. Domain decomposition approaches to solv
ing the Helmholtz equation are in papers by Hetmaniuk and Farhat [34], and
Tezaur, Macedo and Farhat [59]. Analysis and domain decomposition methods
for the time dependent coupled fluidsolid interaction problem are found in pa
pers by Cummings and Feng [14, 24]. Mandel in [50] presents a domain decom
position method for solving the time harmonic coupled fluidsolid interaction
5
problem. Ohayon and Valid [55] present symmetric variational formulations for
the problem of transient and modal analysis of bounded coupled fluidstructure
linear systems, taking into account gravity and compressibility effects. Gener
alized added mass operators, independent of time and circular frequency, are
introduced. Numerical results are presented for an incompressible hydroelas
tic modal analysis of a liquid propelled launch vehicle, elastoacoustic modal
analysis of an incompressible structure containing a compressible gas, and an
elastic cylinder partially filled with liquid under gravity effects. In [22], Elman
et. al. present a multigrid method enhanced by Krylov subspace iterations
for the discrete Helmholtz equation. In this thesis, we extend the approach
of [22] to the coupled problem. Makridakis et. al. [48] proved existence and
regularity of the solution of elastic scattering in bounded domains in one di
mension, using the Garding inequality and an infsup condition, which leads
to asymptotically optimal estimates. Our results extend [48] to more than one
dimension; however, some tools which work in one dimension are not available
here, so we do not obtain asymptotically optimal estimates. More details on
the approaches outlined above follow.
Bielak, MacCamy, and Zeng [8] discuss a scattering problem and its
solution by a coupling method. Existence, uniqueness, convergence, as well
as accuracy of the numerical approximations is also presented. This is ac
complished by using potential theory and three different representations of the
pressure p to solve the coupled fluidsolid problem. The domain is unbounded,
6
and it is separated into a bounded region Cl in J23, with boundary T, and exte
rior Cl+. The domain Cl represents an inhomogeneous elastic obstacle and Cl+
is a compressible, non viscous, homogeneous fluid. Three different coupling
problems are presented, of which, the first two both fail for different sets of
frequencies. The third coupling method is based on potential theory, and is
shown to be stable.
Another approach to the problem of scattering by an elastic obstacle
in an acoustic fluid was presented by Demkowicz [15, 16]. The main difference
between what Demkowicz did and this thesis is that Demkowicz investigates
rigid scattering and vibrations of an elastic submerged shell. From the spectral
decomposition of the operator for an elastic spherical shell in fluid, he com
putes the LBB constant as a function of the wave number k, and shows its
effects on convergence. This work is motivated by the earlier paper [17, 18],
where asymptotic convergence was studied for both finite and boundary el
ement methods on acoustic problems, and a 1D model acoustic interaction
problem is presented. Elastic scattering problems investigated by Demkowicz
in [15, 16] assume an elastic spherical shell freely floating in fluid. The spectral
decomposition for this operator is computed, and finding the LBB constant
reduces to solving a saddle point problem. Pointwise infimnm of the spectral
decomposition show dependence of the LBB constant on the wave number k.
Demkowicz is able to prove that the magnitude of the LBB constant depends
upon the distance from the nearest resonant frequency, and that without a
7
strict control of the discrete LBB constant during the solution process, the
results may be unreliable.
Cummings and Feng [14] present a domain decomposition method for
the time dependent system of coupled acoustic and elastic interaction prob
lems. Two classes of iterative methods are proposed for decoupling the domain
problem into fluid and solid subdomains, and replacing the physical interface
conditions with equivalent relaxation conditions as the transmission conditions.
The nonoverlapping domain decomposition methods developed are regarded as
Jacobi and GaussSeidel type algorithms, and they use convex combinations
of the original physical interface conditions to transmit information between
the subdomains. Strong convergence in the energy norm for the fluidsolid
interaction problem is shown for the iterative methods, and their findings are
supported by numerical test results.
Feng [24] analyzes some finite Galerkin approximations for the time
dependent fluidsolid interaction model, and presents a domain decomposition
method for the time dependent system of coupled elastic acoustic problem. An
optimal order apriori error estimates in L(H1)norm and in L(L2)norm for
the semidiscrete and fully discrete Galerkin approximation to the solution of
the model is established. Higher order time derivatives of the errors are used
for the error estimates of the interface conditions describing the interaction
between the fluid and the solid. To handle the terms involving the interface
conditions, the boundary duality argument due to Douglas and Dupont [38]
8
is used. The error estimates for the fully discrete methods are obtained by
averaging error equations at different time steps and choosing some nonstan
dard test functions. Feng defines a second order discretetime Galerkin method,
and presents a generalized parallelizable nonoverlapping domain decomposition
method for solving the coupled fluidsolid interaction problem.
Another domain decomposition method for time harmonic coupled
fluidsolid systems that decomposes the fluid and solid domains into nonover
lapping subdomains is presented by Mandel in [49]. The continuity of the
solution uses Lagrange multipliers to ensure that the values of the degrees of
freedom coincide on the interface between the subdomains. This method is
known as the FETIH domain decomposition method originally proposed by
Farhat and Roux [23]. Mandel finds that the division into subdomains does
not need to match across the wet interface. The system is augmented by dupli
cating the degrees of freedom on the wet interface, and the original degrees of
freedom are eliminated. The intersub domain Lagrange multipliers and the du
plicates of the degrees of freedom on the wet interface are retained and form the
reduced problem. The resulting system is solved by iterations preconditioned
by a coarse space correction. Numerical results are presented for a bounded
2D region.
Elman et. al. [22] study the exterior Helmholtz problem in an un
bounded domain. The unbounded domain is truncated to a finite domain
9
by introducing an artificial boundary on which the radiation boundary condi
tion approximates the outgoing Sommerfeld radiation condition. In this pa
per multigrid methods are used to solve the discretized Helmholtz equation.
The authors identify difficulties arising in a standard multigrid iteration for
the Helmholtz equation, and analyze and test techniques designed to address
these difficulties. Some difficulties encountered are with the smoothing and
coarse grid corrections. Standard smoothers such as Jacobi and GaussSeidel
relaxation become unstable for indefinite problems, since there are error com
ponents that are amplified by these smoothers. The difficulties with the coarse
grid corrections are due to the poor approximation of the Helmholtz operator
on very coarse meshes. The approach used in [22] for smoothing is to use stan
dard damped Jacobi relaxation when it works reasonably well, on fine enough
grids, and then, to replace it with a Krylov subspace iteration when standard
damped Jacobi fails as smoother. For the coarse grid correction, the number of
eigenvalues that are handled poorly during the correction is identified, and an
acceleration for multigrid is introduced by using multigrid as a preconditioner
for an outer Krylov subspace iteration. The authors observe that multigrid
does a poor job of eliminating some modes from the error, so it converges
slowly or even diverges in some cases, and an outer Krylov subspace iteration
is needed for the method to be robust. GMRES is used as the Krylov sub
space method. In this thesis we extended the approach of [22] to the coupled
problem.
10
The analysis for a fluidsolid, interaction problem in one dimension
done by Makridakis, Ihlenburg, and Babuska [48], is the closest to the analysis
presented in this thesis. The approach taken in [48] focuses on the stability
of the continuous problem (2.1) and on the stability and convergence of the
discrete problem (2.9) with respect to the wave number k. The problem is a ID
layered fluidsolidfluid medium with configuration = fii U Â£l2 U = [0, L],
L > 0,
Pxx + k2p = <71, in fti
(aux)x + k2gu = /, in 0,2 (2.1)
Pxx + k2p = <72, in ft3
with radiation boundary conditions at the boundary of fl,
Px(0)+ikp(0) = 0,
px{L) ikp(L) = 0,
transmission conditions on the interface boundary at x\,
px{x 1) k2u(xi) = 0,
p{x 1) + aux(x 1) = 0,
and at x2,
Px{x2) k2u(x2) = 0,
p(x2) + aux(x2) 0.
11
The space considered is H = x i71(^2) x i?1^). With the weak
formulation, we find U = (p, u,p) Â£ H such that
B(U,V) = (F,V)0,h WÂ£H (2.2)
where V (q, v, q) Â£ H,
B(U, V) PxQxdx ~k2 pqdx k2u(x1)q(xi) ikp(0)q(0)
+k2 / auxvxdx k4 guvdx k2p(xi)v{xi) + k2p(x2)v{x2) (2.3)
n2 n2
+ f Mxdx k2 [ p\dx + k2u{x2)n{x2) ikp(L)ftL)
n3 q3
and
(F, V'io.H = (gi, q) + k2(f, v) + (32, q).
(2.4)
The standard definition for the inner product is used, (u, v) = f uv dx. With
a suitable choice of norms,
(U,V)0,h = (p,q)+k2(u,v) + (p,q), (2.5)
(U,V)hH = (Ux,Vx)0,H + k2(U,V)o,H, (2.6)
it was shown using variational techniques, that if the righthand side F Â£ L2,
then the solution satisfies a regularity estimate of the form
Â£/i,*<^11*110,*, (2.7)
where C\ is a constant independent of k. The estimate (2.7) establishes the
uniqueness of the solution of (2.2), and is needed in the analysis of the discrete
12
problem. Uniqueness combined with the fact that the variational form satisfies
a Garding type inequality, yields the existence of the solution U. Using an
estimate of the form (2.7) for a properly chosen auxiliary problem, the Babuska
Brezzi condition for the bilinear form B(, ) is proved:
sup Re^f'V) >7^11.,H, w H. (2.8)
veH K
This approach follows earlier Babuska work relying on the LBB constant. LBB
condition (2.8) is equivalent to a bound on the norm of the solution operator,
which depends on k linearly. Using (2.7) and (2.8), other similar regularity
estimates are derived. These estimates are useful in the convergence analysis
of the numerical methods.
If Sh is a suitable finite element subspace of H consisting of piecewise
polynomial functions, approximation Uh E Sh to U is defined as the solution
of
B(Uh,(j)) (F,E Sh. (2.9)
In [48], the authors showed the existence of a unique solution of (2.9),
provided the quantity ^k2 is small enough, where h is the maximum mesh size
and d is the polynomial degree of the functions of Sh It is also shown that the
discrete analog of the LBB condition (2.8) is satisfied on Sh The authors used
(2.7) to show that an optimal estimate of the form
IIu Uh\\^H < (7* inf IIU 0i(2.10)
13
holds, where C* is a positive constant independent of h and k. The approxi
mation properties of Sh are known, thus (2.10) implies convergence of optimal
order of the finite element approximations in the  i,rr norm. In ID, the final
estimate has no dependence on wave number k. In more than ID, the equiv
alent of (2.8) does not hold; instead we use estimates and FEM error bounds
that rely on the Garding inequality. We have dependence on wave number
k, because the estimate relies on the boundedness of (/ (3k2G)~l, where (3
is a constant independent of k and G is a compact operator with finite but
unknown norm.
14
3. Theoretical Preliminaries
In this chapter, we describe notation and summarize some concepts
and algorithms used throughout this thesis. Let Q, be a bounded open con
nected domain in Rn, n E 1,2,3, with Lipschitzcontinuous boundary T. Let
dfi = rdurn with rd, rn disjoint, and meas(Trf) > 0. We will denote by rd, Fn
the parts of the boundary with Dirichlet and Neumann boundary conditions,
respectively.
3.1 Sobolev Spaces
Let Wq(Â£l) be the usual Sobolev space of complex valued functions
IV*(Q) = {/ e IL : ll/ll?(n) < oo} ,
where q is a nonnegative integer, and \\f\\w';[C,i = (E*a, " ^'
We use the multiindex notation a for denoting partial derivatives by an n
tuple with nonnegative integer components, a = (ax, ...,an), with the length
of a given by \a\ = Y2= i ai For p = 2, the space Wq(fl) is a Hilbert space,
and we denote W% (Q) by Hq(Â£l). For q = 1, we obtain the Hilbert space Hl
equipped with the norm \\f\\Hi = (/ia(n) + lV/2(n))1/2. In addition to
integerorder Sobolev spaces, there are fractionalorder Sobolev spaces. For
15
seR, 0 0, we recall
Mlwp9+S(fi)
where Cl C Rn. For more details, see [2, 10, 54].
3.2 Finite Element Approximation
The finite element method (FEM) is a general technique for the nu
merical solution of partial differential equations in structural engineering. From
the engineering point of view, the method was thought of as a generalization
of earlier methods in structural engineering for beams, frames, and plates,
where the structure was subdivided into small parts, called finite elements,
with known simple behavior. The presentation here follows [10], where more
details can be found. To fix ideas, as a simple example, consider the second
order elliptic equation with Dirichlet boundary condition,
Au = f in Q
(3.1)
u = 0 on dCl,
where Cl is a Lipschitz domain in R2 or R3, and
16
with the matrix [a(i,j)\ symmetric, positive definite, and bounded on fh Mul
tiplying by a test function v G V(fi), where the Sobolev space V(Q) = Hl(Q),
is a given set of admissible functions, and integrating by parts, the model
problem (3.1) can be written in the variational form,
Find u G V(Q) such that a(u,v) = F(v) for all v G 17(f)), (3.3)
where F : V > R is a functional with F(v) = (/, v), and the bilinear form
a(, ) : V x V > R given by
, \ f du dv du
, a(u,v):=J
^ ij=1 ux3 k=l UJ'k
du
v + bouv dx
(3.4)
is defined for all u and v in the Sobolev space V(fi) = i?1(fl).
The functions v G V(Q) usually represent a continuously varying
quantity, such as a displacement in an elastic body, and F(v) is the total en
ergy associated with v. In general, the functions in V(Q) cannot be described
by a finite number of parameters, and so the problem cannot be solved directly.
The standard Galerkin approximation approach is to look for an approximate
solution of (3.3) in a finite dimensional subspace Vh(Q) of the space V(fl), in
which the weak form is posed. The space Vh (fi) consists of simple functions
only depending on finitely many parameters, usually chosen to be piecewise
polynomials. This leads to a finitedimensional problem. The Galerkin ap
proximation is the solution to the following problem:
Find Uh G 14(fl) such that a(uh, Vf,) = F(vh) for all Vh G 14(Q). (3.5)
17
When a basis is chosen for 14, the Galerkin approximation leads to a system
of equations. Let {4>i\f=i be a basis for 14(H). Assuming that
N
i1
equation (3.5) becomes
N
a(<^> fa)** = (/> fa)i j = !> I N
2=1
The lefthand side matrix K = [a(0i, j)\ is called the stiffness matrix, the
righthand side vector [(/, (pj)\ is called the load vector, and the vector of
unknowns [if*] is referred to as the vector of degrees of freedom.
For a wide class of approximation spaces 14(H), Uh is a good ap
proximation for u. The choice of the finite dimensional subspace 14(^) is
influenced by the variational formulation, accuracy requirements, and regular
ity properties of the exact solution. The construction of suitable spaces 14 uses
triangulation, which splits the domain into small disjoint simple geometries.
In R2, triangular and rectangular shapes are considered, and for R3, tetra
hedrons and hexahedrons are used. By imposing certain assumptions on the
triangulation, the finer the triangulation, the closer a finite element Galerkin
solution is to the exact solution.
As in the example above, functions in 14 arise from a polynomial
interpolation on the elements of the triangulation, and are generated by basis
functions that are usually polynomial on each element of the triangulation of
the domain. The supports of basis functions have only small overlaps. Every
18
polynomial defined on a given region is uniquely determined by its values and
perhaps the values of its derivatives at some nodal points. Thus, each function
in Vh is determined by a set of values at nodal points, so called degrees of
freedom. Simple examples of finite element spaces include spaces formed by
continuous linear or triangle regions in R2, or tetrahedrons in R?. These spaces
are referred to as Pi and Qi respectively. In this thesis we use standard Qx
finite element spaces.
For finite element spaces, the standard basis functions chosen are
the continuous piecewise linear functions that take the value 1 at one node
point and the value 0 at other node points. Thus, the unknowns in the linear
system arising from the discretization are the degrees of freedom of the Galerkin
approximation. Since the overlaps of the supports of the basis functions are
small, the stiffness matrix is sparse.
3.3 Garding Inequality
A bilinear form a(, ) on a normed linear space H, is said to be
bounded (or continuous), if 3cx < oo such that
\a(u,v)\ < CilltiHfflMIfl
and coercive on V C H if 3 C2 > 0 such that
\a(v,v)\ > C2\\v\\2h WveV. (3.6)
There can be wellposed elliptic problems (3.1) for which the corresponding
variational problem (3.3) is not coercive, although a suitably large additive
19
constant can always make it coercive as follows. By Garding inequality there
is a constant, k < oo, such that
a(v,v) + kM2(q) > a\\v\\2HHn), (3.7)
for some a > 0, where a(, ) is defined in equation (3.4).
3.4 Generalized Korn Inequality
Korns inequality plays an important role in establishing the existence
and uniqueness of a solution in linearized elasticity, and it is used to establish
coercivity of the operator. Let Q G Hn be a domain with Lipschitz boundary.
There exists a constant C* > 0, such that for all u G [40, 53]
lle(M)ll(L2(Q))x + IMI(L2(0))" ^ ^k\\u\\fH!^n, (3.8)
where e(u) = (Vu + (Vu)T) is the strain tensor (see section 4.1).
3.5 Fredholm Alternative
Compactness of a linear operator is essential in Fredholms theory. For
X and Y normed spaces, an operator T : X > Y is called a compact linear
operator if T is linear, and if for every bounded subset M of X, the closure
T(M) is compact, i.e. every sequence in T(M) has a convergent subsequence
whose limit is an element of T(M).
A bounded linear operator A : X > X on a normed space X is said
to satisfy the Fredholm alternative [42] if A is such that either (I) or (II) holds:
(I) The nonhomogeneous equations Ax = y, Ax f = g, have solutions x and
20
/, respectively, for every given y E X and g E X' the dual space of X, the
solutions being unique.
(II) The homogeneous equations Ax = 0, Axf = 0, have the same number of
linearly independent solutions x\,...,xn and respectively. The non
homogeneous equations Ax = y, Ax f = g, have a solution if and only if y and
g are such that fk(y) = 0, g(xk) = 0 (k = 1,..., n), respectively. The particular
important case is that in which both Ax 0 and Af = 0 admit only the trivial
solution; then there is a solution of Ax = y for any y E Y.
Compact operators and the Fredholm alternative theory are related
by the following result. Let T : X Y be a compact linear operator on a
normed space X, and let A ^ 0. Then T\ = T XI satisfies the Fredholm
alternative.
3.6 Trace Theorem
Let be a bounded open set of Rn with a Lipshitz boundary dll and
T a subset of the boundary. Consider the Sobolev space Hm(fl); the trace yu =
Ur of a function u 6 Hm(Il) on T is a bounded linear operator, and is defined
as a restriction of the function u to the boundary 7 : Hm(Q) > Hm~1/2(T),
m> \ [27]. There is a constant C'r(m), such that
(3.9)
21
3.7 Riesz Representation Theorem
Any continuous linear functional F on a Hilbert space H, can be
represented uniquely in terms of the inner product, F{x) {x,z), where z
depends on F, is uniquely determined by /, and has norm \\z\\h = \F\\#/.
In general, let Hi, be Hilbert spaces and, h : Hi x H2 K a bounded
sesquihnear form. Then h has a representation (x,y) = (Sx, y), where S :
Hi > H2 is bounded linear operator, uniquely determined by h, and has norm
m = mi
3.8 LaxMilgram Theorem
Given a Hilbert space (V, (, )), a continuous, coercive bilinear form
a(, ) and a continuous linear functional F Â£ V', there exists a unique u Â£ V
such that a(u,v) = F(v), \/v Â£ V. Since the bilinear form a(, ) is coercive,
we have the estimate wy < A F(y/, where c2 is the coercivity constant (3.6).
3.9 Hilbert Interpolation Spaces
There are many ways in which one can define an interpolation space,
with the most common methods being the real and the complex interpolation
methods. We will describe the complex interpolation method, and we will
follow [41], for the real interpolation method [10, 61].
Here and in section 5.2, we use the notation Xa X& to denote the
continuous embedding of Xa in Xb (i.e. convergence in Xa strongly implies
convergence in Xb strongly), [58]. Similarly, we let Xa Xb denote compact
22
embedding of Xa in Xb (i.e. the embedding operator is compact). We are
interested in the interpolation theory of Hilbert spaces. Given two Hilbert
spaces Xq and X\ ^ X0 with inner products (, )x0, (, )xi and norms  11*0 >
 xi; respectively, we will define Hilbert spaces that interpolate between
them. The spaces Xg = [X0, X\\g, 0 < 6 < 1 are called interpolation spaces
between X0 and X\.
Let an unbounded positive definite selfadjoint operator A exist in
X0, with dense domain D{A) = X\, such that H^Hx! = >1xoj u Â£ D{A).
The existence of operator A follows from the Riesz representation theorem.
For 0 < 9 < 1, the set
[X0,Xi]e = Xg := ju G X0 : M[x0,Xi]fl < (3.10)
forms a Hilbert scale of spaces (Xg,  xj given by Xe = D(Ae), with norm
IMI[Xb,*i] = ll^wIUo
The following theorem is known as the Convexity Theorem, and will
be used later.
Theorem 3.1 [10] Suppose that Xi and Yi (i 0,1) are two pairs of Banach
spaces, and that A is a linear operator that maps Xi to Y^. Then A maps Xg
to Yg for 0 < 6 < 1. Moreover,
IMky. < M&YMWkYi (3n)
23
4. Statement of Coupled Problem
4.1 Derivation of the Coupled Problem
In this section, we summarize some of the basic relations of linear
wave physics, starting with elastic waves, proceeding to acoustic waves, and
then fluidsolid interaction. Derivation of the coupled problem is standard, and
this section is intended for completeness only. We will derive the differential
equations and the boundary conditions that describe the fluidsolid interaction.
We are interested in the timeharmonic case, and we assume that all waves are
steadystate with circular frequency u. An acoustic wave that is incident onto
an elastic obstacle is not totally reflected, and part of the incident energy is
transmitted in the form of elastic vibrations. Basic ideas regarding the nature
of the elastic field were established long before the application of elasticity
theory [64]. We will use the word field to refer to any physical quantities
that vary in space and time. For time harmonic fields, the scalar and vector
potentials satisfy the scalar and vector Helmholtz equations.
4.1.1 Elastic Waves In an elastic medium, waves propagate in
the form of small oscillations of the stress field. In a solid elastic material,
the physical quantities that are of interest are the displacement, the stress, the
24
strain tensors and, if any are present, the body forces. The speed of propaga
tion is usually denoted by c. The governing equations for the elastic medium
are obtained from the basic relations of continuum mechanics. Consider the
momentum equation
+ V (pvv) V t = f.
Assuming no external forces, the momentum equation becomes
+ V (pvv) V r = 0.
Assume that the density, p(x,t), and the velocity, v(x,t), undergo small fluc
tuations about equilibrium values, po and v0, respectively. The domain as a
whole is assumed to be at rest (i.e. v0 = 0). We linearize the momentum
equation by writing
V = V + v0
p p + po
where v, p represent small perturbations about the equilibrium values. The
term pvv can be linearized as
pvv = p(v + v0)(v + v o)
= p(vv + VVq I V0V + VqVq) 0,
25
since vv is negligible and u0 = 0.
Similarly,
pv = (p + po)(v + v0)
= pv + pv o + Pqv + p0v0
~ Pov,
since pv is also assumed negligible. The momentum equation then becomes
d(p0v)
dt
V r = 0,
(4.1)
where we have dropped the tilde notation for simplicity. Since po is the constant
equilibrium value, equation (4.1) becomes
dv
^wv'T=a
(4.2)
To linearize equation (4.2) and write it in terms of displacement we assume
small oscillations. The partial derivative, J^, is related to the total derivative
^ by the nonlinear expression [43]
d d
H = dt+V'v'
If u(r,Â£) is the vector field of particle displacement at position r, then in order
to formulate equation (4.2) in terms of displacement, we write
/ \ du. ,
(r,i) = (r ,t)
du
~dt
du
~dt
+ v Vm
26
where Vn is assumed to be of the same order as v so that v Vm is negligible.
So,
dv d2u
~dt = W'
(4.3)
and we obtain the linearized form of the momentum equation in terms of
displacement
d2u
PQW
 V r = 0.
(4.4)
For wave propagation, it is generally assumed that a timedependent scalar
field F(x,t) can be separated as
F(x,t) = f(x)erM,
where / is a stationary amplitude function. A similar convention holds for
vector fields. Assuming a time harmonic solution to equation (4.4) and letting
u(x,t) = u{x)e~'UlJt, we obtain
where r is the stress tensor. The symmetric stress tensor r(x) is also known
as the Cauchy stress tensor at the point x 6 fie.
With the assumption of small deformations, the strains are related to
the displacements by the linearized equations, and are defined by
e(u) = hvu + (Vu)T). (4.6)
z
27
By Generalized Hookes Law, stress is a linear function of strain,
where the strain assumes small displacements. We then have
(4.7)
where Cijki is a fourth order elastic stiffness tensor which, for an isotropic
medium, is invariant under rotations and reflections, so that it takes on the
form
where A, p are the Lame coefficients. Substituting equation (4.8) into (4.7) we
get
7" ij X5ij Cfck 4
XSijdkUk 4 /7(iij 4 djUi).
Thus we obtain the equations governing the elastic medium
u2pu + V t = 0
r = AI(V u) + 2pe(u),
where r is the stress tensor, e(u) is the strain tensor, u is displacement, p is
density of the elastic medium, and A, p, are the Lame coefficients of the elastic
4.1.2 Acoustic Waves Acoustic waves (sound) are small oscil
lations of pressure in a compressible ideal fluid (acoustic medium) [35], and are
Cijki ASijSki 4~ fJtSjiSjk
(4.8)
medium.
28
associated with local motions of the particles of the fluid and not with bodily
motion of the fluid itself. The difference between fluids and elastic solids is
that fluids cannot support shear stresses in the absence of internal friction.
In viscous fluids, frictional forces are generated when gradients in velocity are
present. Here we consider inviscid media, so the fluid cannot support any
shear forces. We assume that the fluid is compressible, i.e. the density of the
fluid changes as a consequence of flow processes. The equations governing the
acoustic medium are obtained from fundamental laws for compressible fluids.
In this section, the differential equations governing acoustic wave propagation
in a liquid or gaseous medium are described following the derivation in [65].
As with the derivation of the differential equations governing the elas
tic medium, we assume that the density, velocity, and pressure, undergo small
fluctuations about their respective mean variables: Po,Vq, and po We also
assume that the fluid as a whole is at rest, i.e. vQ = 0. The simplest con
stitutive equations encountered in continuum mechanics are those for an ideal
fluid, where the pressure field is isotropic and depends only on density and
temperature; then the stress field r is represented by
T = p(p,T) I. (4.9)
where p is the density of the fluid and T is the absolute temperature.
Substituting equation (4.9) into the linearized momentum equation (4.2) we
29
obtain
Pof' V (pi) = 0.
Taking the curl of equation (4.10), we see [65]
(4.10)
Vxv = 0 => v = V4>,
(4,11)
where $ is a scalar field called the velocity potential. Equation (4.11) holds
because of the assumption of a nonviscous fluid. Using equation (4.11) in
equation (4.10), we obtain
dÂ§
P = Po
dt
(4.12)
From thermodynamics, we can write the pressure as a function of entropy and
density. We assume that the acoustic wave propagation is an adiabatic process
at constant entropy, and that the changes in density are small. Then using the
Taylor series representation, we have
dp'
p=l^sdp
(4.13)
where p' p0 + p. We define the adiabatic compressibility Ks via
'dp'
 c
(4.14)
\dpjs Kspo
where c is the speed of the acoustic waves and depends on material properties.
Using equation (4.13), (4.14), and (4.12) in equation (4.10), we obtain the time
dependent wave equation for the velocity potential
V2$
1 d2$
c2 dt2
= 0.
(4.15)
30
With the assumption of timeharmonic waves of frequency u, the wave equation
(4.15) can be transformed to the scalar Helmholtz equation
V2fiVt = 0, or
c2
V2V + k2V = 0 (4.16)
where ^ is amplitude, i = >/T is the imaginary unit, and
k = ~. (4.17)
c
The physical parameter k is called the wave number. The physical interpre
tation of the parameter k is the number of waves per 2tt units. Hence, k
characterizes the oscillatory behavior of the exact solution. The larger the
value of k, the greater the spacial frequency of the waves.
4.1.3 Boundary Conditions In order to set up a wellposed
problem and to model the physics, the equations need to be complemented by
boundary conditions. The most common boundary conditions imposed on a
model problem are the Dirichlet and Neumann boundary conditions. Holding
the pressure constant on the boundary is an essential, or Dirichlet, boundary
condition of the form p = p0. For a rigid surface, Neumann boundary conditions
are imposed, ^ = 0, where n is the outer unit normal to dCl. This condition
models a free boundary, i.e., no external forces are acting at the boundary [58].
The physical requirement that all radiated waves are outgoing leads to
the Sommerfeld radiation condition, which can be interpreted as a boundary
condition at infinity. The Sommerfeld condition is usually replaced by its
31
approximation on an artificial boundary of a finite computational domain Q,
is is given by
dp ., n
t;1" 'IKP = 0.
on
This condition is also called a radiating boundary condition [58].
4.1.4 SolidFluid Interface Conditions When an elastody
namic wave, which propagates through a homogeneous medium, encounters
a region that is characterized by different physical properties, it will experi
ence changes in amplitude, wave numbers, and direction of propagation. We
want to impose physically motivated conditions for the description of solid
fluid interaction. Most often, boundary conditions in elastodynamics require
the displacements and surface traction to be continuous across the boundary
region. The traction vector for an arbitrary region with unit outer normal n
may be written as
t = r n, (418)
where r is the stress tensor as in equation (4.5). Let ue,Uf,te,tf denote the
displacement and traction in the elastic and fluid medium, respectively. We are
dealing with an inviscid fluid, so we must allow slip in the direction tangential
to the interface. Then the only condition on u is the continuity across the
interface,
n ue = n Uf. (419)
32
The displacement of the fluid is obtained by eliminating Uf between the lin
earized momentum equation (4.10) and equation (4.3). Assuming a time har
monic solution, we obtain
Uf = ^Vp. (4.20)
PfLO2
Combining equations (4.19) and (4.20) we obtain the first interface condition
n ue =  Vp n. (4.21)
Pfu2
The last two interface conditions can be obtained from the the balance of forces,
namely the pressure is in static equilibrium with the traction normal to the
solid boundary,
n te = n tf.
From (4.9) and (4.18), we obtain
n r = pn. (4.22)
Crossing both sides of equation (4.22) with n, gives
n x t n = 0. (4.23)
Hence, the fluid cannot support any shearing force. By dotting both sides of
equation (4.22) with n, we obtain
nr n = p (4.24)
which means that at the wet interface, the fluid pressure is in equilibrium with
the normal traction of the solid.
33
Tn
Me . n
Q f
rn
Figure 4.1. Problem setup
4.2 Variational Form of Coupled Problem
We consider the propagation of waves in a composite medium Q C
Rn ,n = 2 or 3, which consists of a fluid part Clf and a solid part Qe, so that
= Â£lf U f2e. Denote the interface between the two media by V = ddj fi dVte,
where the boundary V is assumed piecewise smooth, and the outward unit
normal is denoted by n.
Our model problem is a channel with rigid walls and an elastic obsta
cle in the middle, Dirichlet boundary condition for the incoming wave on Td,
Neumann boundary conditions at rigid surfaces on the sides on Tn, an absorb
ing boundary at the outgoing side Ta, and interface condition on the boundary
Ti of the scattering obstacle, (see Fig. 4.2).
Mathematically, the model is described by the coupled system of the
elastic wave equations and the acoustic wave equation, derived in Section 4.1.2.
The acoustic field in the fluid is governed by the Helmholtz equation for the
34
pressure p, (4.16),
Ap + k2p = 0 in Â£lf, (4.25)
with the boundary conditions derived in section 4.1.3,
= p0 on Yd (excitation),
= 0 on Tn (sound hard surface),
= 0 on Ta (radiation boundary condition). (4.26)
Here, A is the Laplace operator, and k E R is the wave number. The Som
merfeld radiation condition prevents reflection of outgoing waves from infinity,
and the radiation boundary condition (4.26) is the approximation to the Som
merfeld radiation condition on an artifical boundary that delimits the finite
computational domain Q.
Using (4.5) and (4.6) for the isotropic elastic medium we have
V r + u)2peu = 0 in Qe, (4.27)
dp
dn
P
dp
dn
+ ikp
where
r = A I(Vu) + 2pe(u) is the stress tensor,
6ij{u) =  (fy1 + is the strain tensor,
u is the displacement, pe is the density, to is the wave frequency, I is the n x n
identity matrix, and A and p are the Lame coefficients of the elastic medium.
35
The solidfluid interface conditions given by (4.21), (4.23), and (4.24),
are [65]
n U = 12"
pfUJJ oft
n r n = p
n x r n = 0
(continuity)
(balance of normal forces)
(zero tangential tension)
on rv
(4.28)
4.3 Weak Formulation
To derive a weak formulation of the acoustic elastodynamic model
problem, we define the spaces
v, = {(p.pr) I, pe p = 0 on D,}, (4.29)
Ve = {(.., ) (/?'(SI.))3, r (Hi(r))3}. (4.30)
where the restrictions pr = Trrp and r = Trrtt are the traces of p and u on
r.
Multiplying equation (4.25) by a test function (q, qr) G Vf, equation
(4.27) by a test function (v, n vr) G Ve, and integrating by parts over Q we
obtain the following variational form of (4.25), (4.27), and (4.28):
Find (p,pr), (p Po,Pr) G Vf and (u, n ur) G Ve such that
J VpVq + k2 Jpq ik Jprqr J PfUJ2{n ur)qr = 0 (4.31)
i.f $if r
 J (A(V it)(V v) + 2//e(u) : e(v)) + J2 J peu v
fig fie
JPr(n vr) = 0 (4.32)
Vq,qr G Vf and \/v,nvr G K,
36
where po on Td is understood to be extended to a function in //1(0/).
For the weak formulation of the elastodynamic equation, we use Greens for
mula [12] to obtain
/
fie
V
r v
Jr: Vu + J n t v
fie 9fie
Jt : + VvT) + J
n t v.
By multiplying equation (4.32) by pjuj2, the coupled fluidelastic problem be
comes
(v,n vr,pp0,pr) eVexVf :
a(u,nuT,p,pr;v,n vr,q,qr) = ((f,n fr,r,rr), (v,n vr,q,qr))
V(n, n ur,p,pr) e Ve x Vf, (4.33)
where
pq
a(u,nur,p,pr;v,nvr,q,qr) = J VpVqk2 J'
+ik Jprqr + P/^2 J{n ur)qr
ra r
+PfU>2 j (A(V tt)(V v) + 2pe(u) : e(v))
fie
co4pfpe J U V + pfUJ2 JPr(n vT),
r
and the righthand side (f, n fr, r, rr) ho is a functional defined by
((f,nfr,r,rr),(v,nvT,q,qr)}v0 = Jrq + k2jiv
(4.34)
37
+
Jrrqr + j n frn vrU435)
r r
with f G (L2(Qe))3, r G L2(ftf), fr G (L2(r))3 and rr G L2(T) given. Since A
and p, may be very large, we use a scaling of the form u = su' and v = sv',
where s is a scalar, so that the leading terms of equation (4.34) differ by a
factor of k2. For the purpose of obtaining error estimates, this scaling is used
in analysis only and not in computations. For computation, we use (4.43)
below, for better covergence of multigrid. Then
Using (4.17) and by dropping the in equation (4.36), we obtain the variational
form
+PfUJ2s2 /(A (V u')(y v') + 2p,e(u') : e(v'))
(4.36)
+c2pfk2s2 /(A(V w)(V v) + 2p,e(u) : e(v)) (4.39)
38
We choose s so that c2pjs2 max{A, 2p,} = 1. Hence,
s = . 1 . (4.41)
cJpf max{A, 2p}
Replacing V) and Ve with conforming finite element spaces, we obtain the
algebraic system
( Sf + k2Mf ikGf
v pf^sT1
pfuPsT
pfu2s2Se + pepfU>As2Me J
\ ( ~\
p
\u )
R (4.42)
In computations we scale to unit diagonal, and s in (4.42) is of the form
s =. 1 =. (4.43)
ckJpf max{A, 2p,}
In (4.42), p and u are the algebraic representations of p and u, i.e., p and u are
the finite element interpolations of p and u, respectively. The matrix blocks in
(4.42) are defined by
S/ = j VphVqh, Mf = Jphqh, Gf = Jphqh,
cif cif ra
Se = /,A(V"Uft)(VVfc) + 2jue(ttfc):e(wfc),
oe
Me = J uhvh, T = /ph(n vh).
cie r
where we followed the notation of section 3.2. The righthand side R of the
system (4.42) is defined from the zero righthand side of the variational form
with the modification for the Dirichlet boundary condition p = po on T^ as
follows. Listing the degrees of freedom on L first and the remaining variables
39
as second, (4.42) becomes AU = R, where the solution, the righthand side,
and the matrix are, respectively,
U =
uA
\U2j
R =
( Ri A
\R2 /
A =
^ An A12 ^
^21 A22
Here, R\ are the degrees of freedom for pQ, and R2 = 0. We impose the
constraint U\ = R\ and solve for U2 from the equations in the second block,
which gives the system of equations that is actually solved,
\ (
(
\
CO J1 OS ^
V2 J ^ A2iRi f
I 0
0 A22 )
This formulation is standard, [55].
4.3.1 Hilbert scale In the analysis, we will need a scale of
spaces analogous to the intermediate spaces between L2(Q) and H'1(0) for
a scalar problem. Let
V0 = (L2mf x L2(T) x L2(flf) x L2(r),
be a Hilbert space equipped with the inner product
((u,nur,p,Pr),(.v,nvr,q,qr))Vo = J pq + k2 j u v
+ J prqr + / n Mrn ur[,4.44)
r r
where the restrictions pr and are the traces of p and u on T, respectively.
Let
V! = VexVf
40
be a Hilbert space equipped with the inner product
{(u,n uT,p,pr),(v,n vr,q,qr))Vl = j VpVq + k2 j(Vu)(Vv)
Q.f
+k2((u, n ur,p,pr), (v, n vr, q, ^r))v0,(445)
with the associated norms denoted by \\(u,nur,p,pr)\\v0 and (ii, nur,p,pr)\\v,.
In general, define for
Vj = {(u,n uT,p,pr) e (tf^e))3 x H^(T) x x H^(T)
p = 0 on Td, ur e (H^T))3, pr *(r)}, (4.46)
where j > \ so traces exist. Define the scaled norm on Vj,j > \ by
(w,nwr,p,pr)y. = k2\u\2(Hj{Qe))3 + p^(n/) + k2j\\(u, n ur,p,pr)y0,
= k2\u\2fjj + \p\ffj(cif) + k2j+2\M2(L^ne)r
+k2j\\p\\h(cif) + k2j J \pr\2 + k2j J n urf4.47)
r r
41
5. Analysis
A boundary value problem is well posed if, for a given class of data,
the solution exists, is unique, and depends continuously on the data (i.e. it is
stable). In this chapter we will show that the bilinear form (4.40) is coercive
via the Garding inequality, and that it is continuous. We will also show that
if a solution exists, then the solution is unique. By using intermediate spaces
we will show stability of the solution.
5.1 Garding Inequality for the Coupled Problem
coercive, but we will show that adding a sufficiently large constant multiplied
by the Vo inner product makes it coercive over all of V\. Prom (4.44) and (4.45),
the norms on V0 and Vi are
The bilinear form (4.40) associated with the coupled problem is not
(u,nt*r,P,Pr)vi
(u,n ur,p,pr)\\v0
+k2\\(u,n ur,p,pr)\\l0
lVp2(%) + ^Ibllis^)
+&2 Vu(L2(fie))3 + &4MI(L2(nc))3
+k2 J \p\2 + k2 J nw2.
r
r
42
An analogous inequality to the Garding inequality (3.7) is
Lemma 5.1 (Garding Inequality for the coupled problem) Let 0 < k0 <
oo; then there exists (3 < oo and a > 0 such that
Re[a(u,n uT,p,pr\u, n ur,P,Pr)] + /3k2\\(u,n ur,p,pr)\\v0
>a(,riur,pJpr)vr1 (51)
for all (u, n Ur,p,pr) U and k > k0. We can choose (3 independently of k,
and
_ L a , 2(iCk (?pe
a mm l,/3 1, max{A, 2^} ^ max{A,2//}
H___2Mgfc ~ 1) n_c I Pf _] _
k2 max{A, 2p} y max{A, 2p,} j
Proof: We have
Re [a(u,n ur,p,pr;u,n ur,P,Pr)] + (3k2\\(u, n ur,P,Pr)v0
= / \Vp\2k2J p2 + Re
ck
Pf
max{A, 2/x}
J (n u)p
+^max{A,2rf/(AV u? + 2#e() : e(u))
c2pek4
1
max{A, 2/x}
[ \u\2 + Re
ck
Pf
max{A, 2p}
j p(nu)
\(3k2 J \p\2 4 (3kA j \u\2 + (3k2 J \p\2 + (3k2 J \n u\2. (5.3)
n f Qe r r
Using the inequality 2ab > a2 b2, we have
Re
J (n u)p >^J\nu\2^J\p\2.
.r J r r
(54)
43
Using (5.4) and the generalized Korns inequality (3.8) for the integrals over
f2e, it follows from (5.3) that
Re [a(u,n ur,p,pr\u,n ur,P,Pr)] + (3k2\\(u,p)\\
2fiCk
Vo
 IIvpIIÂ£2(%) + (P ~ WMhw + max{A,2/i}^Vw^L2(ne
2 iik2
C2pe
max
+ Pc
+ U max{X,2fi}) k ''w
/W + /'
Pf
max{A, 2/j,}
TO TO
2 pGk
IIVp,(%, + w i)*2lbll!=(%, + ma^{A,;M}fc2llv^ll?^(0.
C2pe
+ 1 P ~ max{\2/i} + F
+ \P ~ C\
Pf
max{A, 2//}/
k2
/w2+/
TO TO
> all (to, to TOr,p,pr)y1
The theorem follows by choosing (3 so that a given by (5.2) with k = ko satisfies
a > 0.
The Garding inequality provides a lower bound on the form (4.40). We also
need an upper bound, which is the subject of the next lemma.
Lemma 5.2 The bilinear form a(u, n ur,p,pr', v, to Ur, q, qr) is continuous
on Vi, and
u(to, to TOr,p,pr; v,n wr,g,gr)
< Cb\\{u, to TOr,p,pr)v1('y, to vr, q, 9r)vu (55)
V(u,nuT,p,pr),(v,nvT,q,qr) 6 Vi,
44
where
C = f 1 I Pf + c2pe
B \ 2 y max{A, 2/i} max{ A, 2/x} max{A, 2/a}
Proof: For any n x n square matrix A [a^], from the Cauchy inequality it
follows that
M! = W = <4
(5.6)
where repeating indices imply summation, and is the Kronecker symbol.
Since V u = tr(e(u)), it follows from (5.6) that
V u\2 = tr(e(u))2 = (Sueu)2 < < tfeye# = ne(u) : e{u). (5.7)
Using (4.40), (5.7), and Schwarz inequality, we obtain:
Ia(u,n ur,p,pr,v,n vr,q,qr)\ =
+<*2 I ^
[ VpVq k2 I
nf (if
pq
max{A, 2/i}
J(n u)q + J p(n v)
+k2m^{\2p} /(A(V ")(V v) + 2/e() : e(v))
2 7,4
 c pe&
1
max{A, 2//}
juv
< lVp]L2(n/)VgL2(Q/) + A:2pL2(n/)g'L2(fj/)
J(n u)g + J p(n v)
+ck
Pf
max{A, 2/x}
+k2 2// + An
max{A, 2/i}
r r
e() II (L2((ie))3\\e(v) II (L2(O0)3
+k
,4 C Pe
max{A, 2/x}
u
(L2(ne))3ll,yll(L2(ne))3
(5.8)
45
Prom (4.6), (5.8), and Cauchys inequality, we obtain
\a(u,n ur,p,pr)v,n vT,q,qr)\ <  VpL2(f2/) VgL2(n/)
+^2blU2(fi/)IMU2(n/)
2V'==fhrf(/'n'2+/lrf
Ifcfh^ (/2+/h2
+fc2
+&5
+*2 , VWIIt^n.)). 11'V(f)ll(J.>(n.))
max{A, 2/x}
+k4
C*pe
max{A, 2/i}
M!(L2(fie))3IMI(L2(ne))3
< CB\\{u,n uT,p,pT)\\W\(v,n vr,?,gr)
Vi
Lemma 5.3 The bilinear form
a{u, n uT,p,pr;v,n vr,q,qr) + pk2((u,n ur,p,pr), (v, n vr, q, qr))v0
is continuous on V\.
Proof: The proof is similar to the proof of Lemma 5.2.
5.2 Existence of Solution
To establish that the variational problem (4.33) has a solution, we
will show that the coupled fluidelastic variational problem has at most one
solution. Then the existence of a solution will be proven by converting the
46
boundary value problem into an equivalent integral equation and invoking the
Fredholm alternative. Uniqueness of the solution will prove existence.
Consider the variational problem
a{u, n wr,p,pr; v,n vT,q,qr)
+/3k2((u, n ur,p,pr), (v, n vr, q, qr))Vo (59)
= ((f, n fr,UUr), (v, n vr, q, qr))v0,
where (f, n fr, r, rr) G V0.
Lemma 5.4 Variational problem (5.9) has a unique solution (w, n ur,p,pr) G
Vi
Proof: Define F by
F : (v,nvr,q,qr) i ((f, n fr, r, rr), (v, n vT, q, qr))v0 (5.10)
Since Vi is Hilbert space, by Lemma 5.1 and Lemma 5.3 the bilinear form
a(u,n ur,p,pr,v,n vr,q,qr) + Pk2((u,n ur,p,pr),{v,n r,g,9r))v0 is
coercive and continuous. Because F is a continuous linear functional, from
the Lax Milgram theorem, (Section 3.8), it follows that there exists a unique
(it, n ur,p,pr) G Vi such that
a(u,n ur,p,pr,v,n vr, q, qr) + (3k2(((u, n ur,p,pr), (v,n vr,q,qr))v0
= F(v, n vr, q, qr) V(v, n vv, q, qr) G V0.
m
Let the linear operator G : Vo > Vi be defined by the solution of the variational
47
problem (5.9)
G : (f,n fr,r,rr) (u,n ur,p,pr) (5.11)
Lemma 5.5 The linear operator G defined in (5.11) is bounded from V0 to
K.
Proof: By Lemma 5.4, the variational form (5.9) has a unique solution
(u,n itr,p,pr) Â£ V\. Operator G maps the righthand side to the solution
(u,n ur,p,pr) Thus proving that G is bounded from Vo to Vi is equivalent
to proving \\(u,n ur,p,pr\\vi < G\\i,n fr,r,rry0. From the LaxMilgram
theorem, definition of the dual norm, and the Riesz representation theorem for
Hilbert spaces, we obtain
\\(u,n uT,p,pr)\\Vl < C\\F\\V> < C\\F\\V> < C(f, n fr,r,rr)Vo.
We use the notation A B as described in section 3.9.
Lemma 5.6 It holds that
14 A V0.
Proof: From the definition of the norms on V0 and Vi,
(w,nwr,p,pr)^ = &2VtifL2(a0)3 + VpÂ£2(fi/)
+k2\\(u,n ur,p,Pr)\\vQ
> C\\(u,nur,p,pr)\\^0,
48
where C = k2. Hence, V\ Vo. Now we show that the embedding is compact.
We need to show that any sequence in Vi has a convergent subsequence in Vo
Let {(wn, (n ur)n,Pn, (pr)n)} be a bounded sequence in V\. Then {un} is a
bounded sequence in (i?1(f2e))3, and has a convergent subsequence {unk} in
(L2(Qe))3. Similarly, since {pn} is a bounded sequence in i71(ri/), there is a
subsequence {pmk} of {pn} convergent in L2(flf). From the trace theorem and
the compact embedding of Sobolev spaces [27],
Tr : H\Q) > H1/2(T) ^ L2(T).
Since {pm,:} is a bounded sequence in i71(0/), there is a subsequence {pik}
of {Pmk} such that the traces are convergent in L2(r), and Tr pik > Tr p in
By the same argument, there is a subsequence {uik} of {w;fc} such that
Tr uik > Tr u in ^L2(r)j .
From CauchySchwartz inequality,
n Tr uik n Tr wL2(r) < jnL2(r)Tr uik Tr ti(i2(r))3.
Hence, nTr Uik nTr u in L2(r). We can conclude that (uik,pik) > (u,p)
in V0. m
Lemma 5.7 The operator G : V0 > Vo is compact.
49
Proof: By Lemma 5.6, the injection / : Vj. > Vo is compact. Prom Lemma
5.5, G : Vq Vi is bounded. Since / : V\ > Vo is compact and G : Vq > V\ is
bounded, the operator I o G : Vq > Vo is compact.
Lemma 5.8 The operator G : Vj V\ is compact.
Proof: Consider a bounded sequence {(un, (n uT)n,pn, (Pr)n)} in V\. Prom
Lemma 5.6, V Vo, so there is a subsequence of {(un, (n ur)n,Pn, (Pr)n)}>
such that {(unfc,(n ur)nk,Pnk, (pr)n*)}  {(o,( ur)o,Po,(pr)o)} in V0.
Since G is bounded from V0 to Vi, it holds that G{unk, (nur)nk,Pnk, (.Pr)nk) >
G(u0, (n ur)o,Po, (pr)o) in V\. Hence, G : V > V\ is compact.
Lemma 5.9 The variational problem (4.33) is equivalent to
{w'tyw0'
where U = (u,p),V = (v, q), and F satisfies (5.10).
Proof: The original variational formulation was
a{u, n ur,p,Pr',n vr, q, qr) = ((f, n fr, r, rr), (v, nvr,q, qr))Vo
or, a(U, V) = (F, V)y0. By adding to both sides the term f3k2{U, V)y0, we find
that
a(U, V) + (3k2(U, V)Vo = (F + (3k2U, V)Vo
50
Then by the definition of G, U = G {F + (3k2U), hence
[wIG)u = WGR
Finally, since G : V0 > V\ Vo, we have F0.
A relation between the eigenvalues of the bilinear form a(, ) and the operator
G follows.
Lemma 5.10 Let Aa Pk2 ^ 0. Then Aa (3k2 is an eigenvalue of the bilinear
form a(u, n up,p,pr; v, n vp, q, ?r); that is
a(u,nur,p,pr,v,nvT,q,qr) =
(Aa f3k2)((u,n ur,p,pr), (v,n vT,q,qr))Vo,
V(v,nvv,q,qr) G Vlt
for some (u, n itp,p,pr) ^ (0,0,0,0), if and only if is an eigenvalue of G.
Proof: From Lemma 5.9 with F = (Aa (3k2)U, we obtain
1 U GU = Xa~?k2GU,
/3k2
/3k2
which is equivalent to jU = GU.
The goal of the following lemma is to show that the variational form a(, ) does
not have real eigenvalues.
Lemma 5.11 Let A G R be an eigenvalue of G with associated eigenvector
(u, n up,p,pr) Â£ Vq. Then p = 0 and ^ = 0 on Ta.
Proof: Assume j is a real eigenvalue of G. Then from Lemma 5.10 the form
a(, ) has eigenvalue Aa (3k2 with eigenvector (u, n up,p,pr) Â£ Vo, and
a(u, nur,p,pr; v, n vr, q, qr) = (Aa(3k2)({u, nur,p,pr), (v, nvT, q, qr))v0
51
Choosing the test functions q p and v = it, we obtain
J I Vp2 + (pk2 Aa k2) J \p\2 + ikj \p\2
:if nf rtt
+ PmaJI{1A,2
+ y3k2 \a k4
+ ck2
c2pe \ r
x{A,2 p})J
u
max{A,
^ l J(nu)p + jp{nu)
max{A, 2p}
+
C0k2 A0) \p\2 + J\n w2j = 0, V(u,n ur,p,pr) e V0.
The terms /(n u)p and / p(n u) are complex conjugates of each other; there
r r
fore, their sum is a real number. All other terms are real numbers with the
exception of ik f \p\2, so ik f \p\2 = 0. Hence, p = 0 on Ta. From the radiation
ra ra
boundary condition (4.26), p ikijr = 0, we also have that = 0 on Ta.
Lemma 5.12 Let p G and Ap + k2p = 0 in Qf. If p = 0 and = 0
on ra, then p = 0 in Qf.
Proof: Since p is a solution of the elliptic equation with constant coefficients
Ap+k2p 0 in Clf, it holds that p is analytic in Clf, [39]. That is, at any point
x of fIf, the Taylor series of p about x converges to p in some neighborhood
Nq(x). Consider x G ra C dVLf. By the CauchyKowalevski theorem [39, 67],
p = 0 is a unique solution in the class of analytic functions of Ap + k2p = 0 in
some neighborhood N0(x) nTa. This solution satisfies the boundary conditions
52
p 0 and ^ = 0 on C dttf. Since the neighborhood N0(x), x G Qf
is compact, there exists a finite subset Ni, of the neighborhoods that
form an open covering of Qf. We will use proof by contradiction to show
that p = 0 on all neighborhoods Assume that p 0 on exactly
n < m neighborhoods. Since the domain Qf is connected, at least one of the
remaining neighborhoods, Nk, has a nonempty intersection with at least one
of the n neighborhoods where p = 0. Since p is analytic and by the Cauchy
Kowalevski theorem, p = 0 on Nk. Hence p = 0 on n +1 of the neighborhoods,
which is a contradiction.
Consider (u, n u^) G Ve, u the displacement vector and itr the trace of u.
Then from the solidfluid interface conditions (4.28), it follows that on Tj,
nu = 0 (5.12)
nTn = 0 (5.13)
n x t n = 0 (514)
Equations (5.13) and (5.14) are equivalent to
t n = 0. (5.15)
Assumption 5.13 If u is solution to the elastodynamic equation (4.27) with
the interface conditions (5.12) and (5.15) on dfle, then u = 0 in fle.
The elastodynamic equation (4.27) along with boundary condition (5.15),
VT + u)2peu = 0 in f^e,
t n = 0 on T, (5.16)
53
does not necessarily imply that u vanishes in Qe. It is known [47] that, for
certain geometries and frequencies, there are nontrivial solutions to this prob
lem. This is equivalent to the assumption that there are no nonradiating
modes [33, 47]. The eigenvalue problem (5.16) has by the Fredholm Alter
native countably many eigenvalues peu2, and its only accumulation point is
infinity. If u2pe is not an eigenvalue of (5.16), then we can conclude that u = 0
on dÂ£le. However, if u2pe is an eigenvalue of (5.16), then the question is if
there exist eigenfunctions u/0 such that u n = 0. These eigenfunctions are
called nonradiating modes. It has been shown [47] that there are bodies such
that nonradiating modes exist. These bodies must have certain symmetries.
An example of such a body would be a sphere. Another example is a shear
stationary wave in half space for which the amplitude is orthogonal to both
the direction of propagation and to the normal, n. It has been shown [33] that
arbitrary elastic bodies have no nonradiating modes.
Lemma 5.14 The operator G has no real eigenvalues.
Proof: If \a /3k2 Â£ jR is an eigenvalue of a(, ), by Lemma 5.12 and Lemma
5.13 it follows that the eigenvector associated with Xa (5k2 must be the zero
vector, which is not possible. Therefore, by Lemma 5.10 j is not an eigenvalue
of G. m
We are now ready to prove the main result of this section.
Theorem 5.15 The variational problem (4.33) has a unique solution for every
given righthand side F Â£ Vq.
54
Proof: From Lemma 5.14, j 7^ 0 is not an eigenvalue of G. The operator
G is compact. Then, (yI G) satisfies the Fredholm alternative. If the
homogeneous problem
Ap + k2p
P =
dp
(in
V r + co2peu =
t n
0 in flf
0 on Td
0 on rOJ
0 in Qe,
0 on r,
has only the trivial solution, then from section 3.5, part (I) holds, and (j^I
G)U F has a unique solution JJ for every given righthand side F G Vq.
We have a standard stability estimate.
Theorem 5.16 Let k be fixed, (u,n ur,p,pr) Vj., and
a(u, n uT,p,pr, v, n vr, q, qT) = ((f, n fr, r, rr), (v, n vT, q, qr))v0 (5.17)
V(v,nvr,q,qr) e Vi.
Then
(uJnur,p,pr)v1 < C'jl(f,nfr,r,rr)vo, (5.18)
holds for all (f, n fr, r, rr)
Cr = (/^2G)1y1^y1Gy0^y1 (5.19)
depending on k and other problem data.
55
Proof: By definition of G, we have
(u,nur,p,pr) = G((f, n fr, r, rr) + /3k2(u, n ur,p,pr))
Hence, by linearity of G
(u,nur,p,pr) = (/ pk2G) 1 G(f, n fr, r, rr).
Consequently,
\\(u,nuT,p,pr)\\Vl
(5.20)
< (I /?A:2G)_1vi_+y1Gvb^v1(f,n fr,r,rr)y0.
From Lemma 5.8, G is compact from Vi to V\. By Lemma 5.14, operator G
has no real eigenvalues, so is not an eigenvalue of G. From the Fredholm
Alternative, it follows that IG is bounded from V\ onto Vi. From the open
mapping theorem [42], G) is bounded from Vi to Vi. Therefore,
(//?A;2G)1y1_iy1
From Lemma 5.5 operator G is bounded from Vo to Vi, and using equation
(5.21) in (5.20) we get the desired result.
5.3 Intermediate Spaces
In this section we introduce intermediate spaces and norms, which are necessary
for the stability arguments in the next section.
56
5.3.1 Intermediate norms For 0 < j < 1, define the seminorm
\{u,n wr,P,Pr)y3. = k2\u\2(Hj{Qe))3 + \p\2Hj{nf). (5.22)
Recall the intermediate norms from (4.47),
\\(u,n ur,p,pr)\\l. = \(u,n uT,p,pr)\l. + k2j\\{u, n uT,p,pr)\\lo,
= k2\u\fHj^e^3 + \p\2Hj(af) + ^2j+2''(i,2(ne))3
+k23\\p\\h(nf)+k23 f \p\2 + k2j J \n u\2. (5.23)
r r
We show that, at least in a certain range of j, the spaces Vj form a Hilbert
scale. The constants we obtain are of great importance, since we want uniform
equivalence for k > oo. We need the following inequality
Lemma 5.17 For O<0l and A > 0,
21(k26 + \2e) < (k2 + \2)6 < {k2e + \29). (5.24)
Proof: Function f(x) = is a differentiable function. Analyzing its first
derivative, fx = ^1+x^i+xe^xwith the constraint 0 < 6 < 1 and x > 0, we
find that x = 1 is an absolute minimum, so f(x) > /(1) for all x, and x = 0 is
a maximum, so f(x) < /(0) for all x. This means that 2e~l < < 1 for
all 0 < 9 < 1 and all x > 0, which is equivalent to inequality (5.24) for x
Lemma 5.18 Let Xg, 0 < 8 < 1, be a Hilbert scale. Define the scaled norms
on Xe by
= \\u\\2xe + k26\\u\\2Xo. (5.25)
57
Then
[(*0, II \\x0,k) pfl, II xi,fe)]0 (Xg, II Wxe,k) ,
(5.26)
with equivalent norms and the constant of equivalence independent of k > 1:
u
Xe,k
> IIu
^o,x0.fc).(^i>IHIx1,fc)]9
> 2
161
u
Xe,k
Proof: Let A be the positive definite, selfadjoint, unbounded operator from
X0 to X(h with domain D{A) = Xx such that uxi = II^IUoj as presented
in section (3.9), [41]. From the spectral theorem [42], there exists a spectral
decomposition of operator
CO
A = J XdE\.
o
Then by definition, we have
IMIxo = / d(Exu,u),
0
oo
IMHi = (A2u,u)Xo = J \2d(Exu,u), (5.27)
0
oo
\\u\\2Xg = \\A0u\\xo = (A20u,u)xo = J \26d(Exu,u). (5.28)
o
From 5.25, we obtain
\u\
2
Xe,k
\\A0u\\2Xo+k20\\u\\2Xo>k
(A20u, u)x o + k20(u, u)xo
= ((A20 + k20I)u, u)xQ
OO
= j(X20+ k20)d(Exu,u).
o
(5.29)
58
The interpolated norm is
IMI
2
((A2 + k2I)eu, u)Xo
(5.30)
= J (A2 + k2)dd(E\u, u).
o
From (5.30) and (5.29) we only need to prove that the norms are equivalent as
presented in section 3.9, i.e.
((.A2 + k2I)9u, uj ((A20 + k29)u, uj . (5.31)
Using the inequality (5.24) in (5.29), it follows that
oo oo
121e(k29 + A 26)d(Exu, u) < J (A2 + k2)ed{E\u, u)
0 0
oo
< J(k29 + X2d)d(E\u, u),
o
which is equivalent with (5.31).
Lemma 5.19 If Xg,Yg, 0 < d < 1 are two Hilbert scales of spaces, then
[Xo X To, x1 x = X0xY0 (5.32)
and
\\{u,v)\\2XexYe = IMIx, + IMI ye (533)
Proof: wx9 = Aflit^0 and Myfl = Bewy0, where Xg = D(A9) and
Yg = D(B9). Now define the operator C such that C : (u,v) i (Au, Bv).
59
Then C9 : (u, v) e> Bev) and
II Mile. =
((Aeu, Bfy;), (u, w)j = (Aeu, u) + (Bev, u)
IMIx, + IMIy, = IKMIlw,
We need the following property of interpolation of a subspace, which is known
from the theory of categories [61]. In the theory of categories, an interpolation
method is called a functor, if it maps the pair of spaces (X0, Xi) to an inter
polated space Xg. We modify and extend category results to get norms with
explicit equivalence constants.
Lemma 5.20 Consider an interpolation functor (Xq,Xi) i Xg such that
(3.11) holds. Let (.Xo,Xl) be an interpolation pair, Xo and X\ Banach spaces,
Y c X0nXi, and suppose there is a linear map T : X0 * X0, T : Xi > X\,
such that Range(T) = Y and Ty = /. Denote by YDXg the space Y equipped
with the Xg norm. Then for 0 < 9 < 1 and x G Y,
fyynxe < IfylliynXo.ynXiie < \\T\\]{o0,Xo\\T\\x1^x1\\x\\YnXg (5.34)
In addition, Y n Xg is a closed subspace of Xg.
Proof: We proceed as in [6] but use norm arguments instead of the theory of
categories. Interpolating the operator I : Y D X0 X0 and / : Y fl X\ Xi,
the convexity property (3.11) gives the lefthand side inequality in (5.34). The
righthand side inequality also follows from (3.11) by interpolating the map
T : Xo * YnX0, T : Xx  Y nXv Finally, ynl=:{xG Xg\(I T)x 0}
60
and T : Xe > Xg is continuous by (3.11); hence, Y n is a closed subspace
of X0.
Lemma 5.21 Let  < a < b < +oo. Then the spaces Vj, a < j < b, form
a Hilbert scale. Define the scaled norms on Vj by \\(u, n wr,p,pr)y, as in
(4.47). Then [(K,  vk)> (H,  yj] = (Vj,  y) with equivalent norms and
the constant of equivalence independent of k > 1, and
\\(u,nur,p,pr)\\Vj < (M,nur,p,jjr)[v.,H]i
< Ci\\(u,n ur,p,pr)\\yj. (5.35)
Proof: It is known that Sobolev spaces LP form a Hilbert scale, [61]. Define
the spaces
Wj = (Hj(ne)f x l2(t) x Hj(nf) x l2(t)
equipped with the norm
\\(u,ur,p,pr)\\2w. = k2\u\2{Hj{ne))3 + \p\2HHQf)
+k2j+2\\u\\2L2(Qe^3 + k2j\\p\\2L2inf).
+^lnr2L2(r) + ^br!2(r)
The space Vj equals ImT, where T : Wj > Wj, T(u,n Wr,p,pr) = (u,Tru
n,p, Tip). Define scales of norms
IMI(ffJ(n.)) = k2\u\2iHj(cle))3 + k2j+2\\u\\2m(le))a, and
\\p\\Hi(nf) = \p\2Hi(nf) + k2j\\p\\2L2{Qf].
61
j/m
Let m> 1. Prom Lemma 5.18 it follows that
[((L2(Qe))3, [ (L2(fie))3,fc) ((Hm(Q,e))3,  (H"*(fie))3,A:) Jj
= ({Hj(fle))3, II (tfj(Qe))3jfe) and
[(L2(fi/),  11^(0/),ft) (Hm(Qf),  Hm(n/)ifc)]^m = (i/*(fy),  llfrjpi,),*) ,
with equivalent norms and the constant of equivalence independent of k > 0.
Since the spaces Wj are a cross product of the above spaces, it follows from
Lemma 5.19, that the spaces Wj form a Hilbert scale. Let a and b be fixed,
\ < a < b. Prom the trace and Cauchy inequalities
lln urf,2(r) < n nr^a_i (r) < Cx(a, r)IM^(a:))3,
lbrli,2(r) < pr^0_i(r) < C'2(a,r)bi/a(n/),
and from the definition of the norm in Wj, it follows that
ll^llwaWa < Ci(a, T), and HTHw^Wt ^(^T).
Hence, by equality of Vj and ImT we have
(5.36)
ll(*.p)llv> = \\(u,Tiunp,Tr:p)\\lmTnWj. (5.37)
Prom Lemma 5.20 for 0 < 6 < 1, it follows for j = a + 6(b a) that
ll(,Tru n,p,Trp)\\imTnWj < (^, Tr w n,p, Trp)[mTniya,imTnw^ <
ll^llwa^W0ll^llwi)VV(1ll('U?Ti:u ,P,^P)limTnwr
The space Wa is isometric with ImT and the space Wb is isometric with ImT,
and by (5.36) and (5.37) we obtain
(u,nur,p,j7r)lbJ < (ti1nr,P,pr)[v,vi]J <
62
lie'll nil1 e\\C2(a,r)\f\\(u,n Ur^.p^llvj,
where the equivalence constants depend on a and T but not on k.
5.3.2 Regularity of Solution for Coupled Problem. Regu
larity results for the Helmholtz equation and elasticity are well known. We get
regularity results for the coupled problem by using known regularity for the
fluid domain Qj and the elastic obstacle domain He.
Lemma 5.22 Let domains Â£lf and f2e be Lipschitz. The solution of the cou
pled problem
Ap = r k2p in Qf,
7P = pd on rd
dpr dn = 0 on r
dpr dn = ikp on ra
Vr = f uj2peu in Â£le,
n ur 1 dpr PfU>2 dn on Tj (5.38)
n Trn = pr on r* (5.39)
n x Tp n = 0 on Tj
exists, where p on is given in H2+Q,  < a < 1. Then (u, n ur,p,pr) is in
Vi+a space with a >
Proof: Let u E (f?1(Qe))3 be solution to the elasticity equation with f E
(L2(f2e))3. Prom the trace theorem, (section 3.6), Ur E (/fs(r))3, so n r G
63
H^(T). From the wet interface condition (5.38), we observe that pn = G
H^{r). Since in general, H% c H~5+ai for \ < aii < 1, then pn G 1T_5+Qi(r).
It is known [27, 54, 30] that the solution of Ap+k2p r with nonhomogeneous
boundary conditions pd G H^+ai(T), pn G H~^+ai(T), and r G (L2(f2e))3, ex
ists and belongs to H1+ai(Qf), for constant  < ai < 1. The following bound
on the norm of the solution holds [27],
lbll/r1+i(S2/) < C1MU1+<*i(n/) + ^Ibdlltf5+i(rd) + (5.40)
where pd and pn are the Dirichlet and Neumann boundary conditions. For
the purpose of this presentation constant C is used to denote a general con
stant that is not the same for all inequalities. Let p G be the so
lution to the Helmholtz equation with righthand side r G L2(fiy); then by
the trace theorem pr G H^(T). From the wet interface condition (5.39),
we then observe that n Tp n G H^(T). There exists  < ai < 1, so that
n Tp n G H~^+a2(T) C We know [28, 29] that the solution of the
nonhomogeneous elasticity operator with constant coefficients with righthand
side f G (L2(fle))3 and normal traction boundary condition n Ty n G
tff+2(r), is regular, and belongs to (H1+c*2 (p,e))3 [30, 54]. There exists
a constant C such that
ll'wll(ff1+2(ne))3 < C]f (ffi+a2(f2e))3 I Cwr(ff_i+Q2(r))3. (5.41)
Because p G by the trace theorem, pr G H%+otl(T). Since u G
(i/1+2(Oe))3, then wr G (H%+a2(T))3 and Mr n G H%+<*2{T). Then we
64
conclude that (u, n ur,p,pr) E 14+a> where a min{o!i, a2}.
We need to assume a stronger property, namely the following bound on the
14+ norm of the solution of the coupled problem.
Assumption 5.23 Assume that there exists a constant Cwr such that
\u,nur,p,pr\Vl+a < CWfif,nfr,r,rrv& (5.42)
5.4 Discretization and Error Bound
In Theorem 5.15, we proved that the variational problem (4.33) has a
unique solution (it, nur,p,pr) for every given righthand side. In this section,
we show that the discrete solution {u^n urh,PinPrh) is also unique. Let 14
be a finite dimensional subspace of 14 as introduced in section (3.2), and define
the Galerkin approximation (uh, n urhiPhiPvh) E 14 to (w, n u?,p,pr) as
the solution of
a{uh, n uTh,Ph,Prh; v, n vr, q, qr) =
((f, n fr, r, rr), (w, n vr, 9, qr)) (5.43)
V(v,nvr,q,qr) e 14 C 14
Let Ie
interpolant of u, and p, respectively [41, 10]. The next lemma gives a bound
on the interpolation error.
65
Lemma 5.24 Let j > If the linear interpolation satisfies
\\U Ie^uW^H1^))3 ^ C'l^W(/irJ+1(^e))3
\\pIf,hP\\Hi(nf) < C2hj\p\Hj+i(Q.f),
(5.45)
(5.44)
then,
\\(u,n ur,p,pr) Ih(u,n wr,P,Pr)lk < CAh(u,n ur,p,pr)\vj+1
where CA = max{C2(max{l, &2} + CrÂ£4(r)), Cf (max{l, Â£;2} + &:2Cr)}. Here Cr
is the constant from the trace inequality (3.9), and p,(L) is the measure of T.
Proof: Let /^(w, n ur,p,pr) be the linear interpolant of (u, n ur,p,pr) on
a mesh with mesh size h. From the trace inequality (3.9) and the assumptions
(5.44) and (5.45), we have
(u,nur,p,pr) Ih(u,nuv,p,pT)\W =
V(u,n ur,p,pr) e V),(5.46)
= k2\\P ~ Ifrfp\\h(nf) + lV(p If,hp)\\h{Qf)
+ k4\\U Ie,hU\\lL2(ne))Z + k2\\V(u Ie,hU)\\2L^ne])3
r
r
< max{l,A:2}p/Mpk1(Q/)
+ k2 max{l, k2} (u Ie,hu) (//i(ne))3
+ k2Cr\pr If,hPr\m(nf) + k2Cr^(T)\u
66
< A;2(max{l, k2} + C'rMr))llw Ie,hu\\\Hi(ne))s
+ (max{l, A:2} + A;2Cr)p Ifthp\\2Hi^f)
< C2h2:>k2(max{l, A:2} + Cr^(T))\u\2Hj+i{ne))3
+ Clh2]{ max{l, A:2} + k2Cr)\p\2Hj+1(nf)
< Cah2j (A:2tt(/fi+i(ne))3 + blfls+qn,))
< CAh2j\(u,n ur,p,pr)\\+1.
To simplify notation we denote by U = (u, n ur,p,pr), V (u, n vr, Q, Qr),
Uh = {uh,n uVh,ph,Prh), Vh = {vhi n vFh, Qh, qrh), W = (w,n wr,y,yr),
and T = (f, n fy, r, ry). We need to assume a bound for the solution of the
adjoint problem, similar to (5.42).
Assumption 5.25 Assume that the adjoint problem
o(V;W) = {{Ul4),V)Vo VV e Vlt (5.47)
has a unique solution W G Vj+i for every V G Vi, and that there exists a
constant Cwr and j > \ such that
Wv5+1 < Cwr\\(U Â£4)v0 (5.48)
We have proved that (5.1), (5.5), (5.18), (5.42), and (5.46) hold. We are
now ready to prove the main result of this chapter. In the following, a is
the constant from the Garding inequality (5.1), Cb is the constant from the
continuity inequality (5.5), CA is the constant from the linear interpolation
67
lemma 5.24, and Cwr is the constant from the regularity estimate (5.42).
They depend on k as follows:
CB
a
CA
fl Â£ / pf + \n c2pe \ =
maX[ 2y max{A, 2/j,} max{A, 2/x} max{A, 2/i} J
j I a 1 2/xCfe ~_____________^ Pe
mm \ max{A, 2p} max{A, 2/x}
2/^(Cfc 1)
k2 max{A, 2/x}
Pf
max{A, 2p}
= 0(1) as k > oo,
max{C2(max{l, A;2} + Cr/x(T)), (^(maxjT, A:2} + k2Cr)} ~ k2,
where p is from Lemma 5.1, independent of A; as A: becomes unbounded, Cr is
the constant from trace inequality (3.9), and p(T) is the measure of T.
Theorem 5.26 Let (5.1), (5.5), (5.18), (5.42), (5.46), and Assumption 5.25
hold. If
khj
then there is a unique solution to
a(Uh; Vfc) = (F, Vh)y0, VVh e Vh, (5.50)
and the solution satisfies
\\{U Uh)\\Vl < Cqinfyheyh (W Vfe)v^, (5.51)
where we may take Cq = 2CB/a, and
\\{UUh)\\Vo < CBCACWRh?\\{UUh)\\vVj > (552)
68
Proof: The proof follows a standard argument [10]. We proved that the
variational problem (4.33) has a unique solution U for every given righthand
side T. First, assume that the discrete solution Uh exists. Since 14 C hi, we
have the standard orthogonality condition obtained by subtracting (5.17) from
(5.43),
a(UUh]Vh) = 0 VVhEVh. (5.53)
From the Garding inequality (Lemma 5.1) and from (5.53), it follows that for
any Vh E 14,
a\\{UUh)fVl < a(UUh,UUh)+(3k2\\(UUh)\\2Vo
= a{UUh]U Uh) + a{U Uh\Uh Vh) + (3k2\\{U Uh)\\2Vo
= a{UUh\U Vh) + (3k2\\(U Uh)\\2Va. (5.54)
Using Lemma 5.2 in inequality (5.54), we obtain
a\\{UUh)\\\ < CB\\(UUh)\\vA(U H)v, + Pk2\\(U Uh)\\\. (5.55)
Then, for any V14 E 14, from (5.53) and Lemma 5.2, applied to the adjoint
problem (5.47), we obtain
\muh)\\2Vo = ((UUh),(UUh))VlB
= a(UUh,W) = a(UUh;WWh)
< CbUUUMKWW^Wv,. (5.56)
From Lemma 5.24, applied to the adjoint problem
IICW MWIk < CAW\W\v, Vj > i (5.57)
69
From Assumption 5.25, substituting (5.48) into (5.57) and using this result in
(5.56), we obtain
(ZYZ4)lk < CbCaCwrVWUh)Ik. (5.58)
Hence substituting (5.58) in (5.55) we obtain,
ot\\(U 74)vi < CB\\{U Vfc)k + (kh?)2[5(CBCACwR)2\\(y Â£4)vi
From (5.49), we find
a\\{UUh)\\v1<2CB\\{UVh)\\v1 Wh Â£ Vh,
and we obtain the desired result. So far, we have assumed the existence of a
solution Uh. By Theorem 5.15, existence and uniqueness are equivalent. If the
solution is not unique, there is a nontrivial solution, Uh, for the righthand side
T = 0. By (5.18), we have U = 0. However, inequality (5.51) then implies
that Uh = 0, provided that h is sufficiently small. Therefore, we conclude that
(5.50) has a unique solution for h sufficiently small.
Theorem 5.27 From Theorem 5.26, there is a unique solution to (5.50) and
the solution satisfies (5.51). If in addition U Â£ Vi+j, j > \, then
\\UUh\\v0 < CBCACWRhj\\U Uh\\Vl < Ch2j\\U\\v1+r
Proof: The result follows immediately by substituting (5.46) into (5.52).
70
6. Multigrid Method
Multigrid or multilevel algorithms are very fast solvers used in numer
ical analysis, physics, dynamics, and computing. In this chapter, we provide
a brief description of multigrid, defining terms, and providing comments on
the structure of multigrid with the goal of setting the stage for multigrid for
coupled fluidsolid scattering. Multigrid has been an active area of research for
almost 30 years, and much literature can be found on the subject.
For the multigrid method, a number of different grids are used on
the domain ranging from coarse to fine. Iterations on different grids reduce
different components of the error: lowfrequency errors are damped on coarse
grids, and highfrequency errors are damped on fine grids. An approximate
numerical solution is computed on a coarse grid, which is used as a starting
point for an iterative method on a finer grid. The method alternates between
smoothers and coarse grid corrections, and the recursive combination of the
two results in the multigrid method. A more detailed description of multigrid
methods can be found in [11, 51].
Using multigrid methods for solving the Helmholtz equation modeling
acoustic scattering is classical [9, 32, 44, 25, 26], but it appears that multigrid
methods for the coupled acousticelastic problem have not been investigated.
71
In this chapter we formulate such method.
6.1 Multigrid for the Coupled Problem
In this section, we present the multigrid algorithm developed to solve
the coupled problem (4.42). This section provides an overview of the imple
mented multigrid algorithm. The setting follows mostly Briggs [11] and Mandel
et. al. [51].
Discretization by finite elements of (4.42), leads to a linear system of
equations Au = /, where the coefficient matrix A is complexsymmetric, i.e.,
not Hermitian. For large values of the wave number k, the system of equations
becomes highly indefinite. This indefiniteness has prevented multigrid methods
from being applied to the discrete Helmholtz equation with the same success
as for symmetricpositive definite problems.
Denote by u the exact solution of the algebraic system, and by v the
approximation to the exact solution; then the algebraic error is e = u v. By
rewriting the original problem, we obtain the residual equation Ae = r, where
r represents residual given by r f Av. Let AiVi = fi denote the system of
equations on grid Ql. Smoothers are simple and inexpensive iterative methods
that reduce the high energy components of the error. The low frequency error
is reduced by projecting the solution onto a smaller space. The grid transfer
operators, the restrictions R and prolongations P, are standard and based on
the finite element space. We will present a variational multigrid method, i.e.
72
the restriction operator is the transpose of the prolongation, and the size of the
prolongation operator P[+1 from coarse level l+l to fine level l is computed from
the size of xe, ye, the number of elements in x and y directions, respectively. All
nodes have three degrees of freedom, pressure p, displacement in x direction,
ux, and displacement in y direction, uy. A node is active, if it has nonzero
row or column contributions. For the fluid nodes only p is active, for the
displacement nodes, ux is active in the displacement in the x direction block,
and uy is active in the displacement in the y direction block. Prolongation P
for the coupled problem for level l has size (3 (1 + xe/2l~l) (1 + ye/2l~1) by
3 (1 + xe/2l) (1 + ye/21)). A temporary matrix Ptemp is assembled by using
the following bilinear interpolation by averaging neighboring values according
to the stencil,
111
4 2 4
 1  . (6.1)
111
4 2 4
which transports the correction obtained on the coarser grid l + 1 to the
fine grid l. We compute the vectors of all odd indices, vx = (1,3, ...,xn),
vy = (1,3, ...,yn), and then use the Kronecker tensor product to create a ma
trix by taking all possible products between the elements of vx and those of vy,
to obtain a large matrix that reproduces all coarse nodes. We use piecewise
bilinear interpolation inside each of the three fields for pressure and displace
ment separately. The elastic obstacle is rectangular in shape in the middle of
73
fluid, and for a coarse node at a straight part of the wet interface we use the
stencil
1
4
1
2
1
4
to average neighboring values. For
we use
1
2
1
1
2
a coarse
0
0 ,
0
(6.2)
corner node of the wet interface,
5 1 0
0 0 0
(6.3)
and for an inside corner of the wet interface we use the stencil
ill
4 2 4
1 1 i
2 x 2
I I o
To obtain the matrix Pfiuid, we zero out the entries in Ptemp that correspond
to the obstacle, while keeping contributions from the wet interface. To obtain
matrix Psoud we zero out the entries in Ptemp that correspond to the fluid, and
leave the wet interface contributions unchanged. Final prolongation matix P
has the form
(6.4)
Pfiuid 0 0
0 Psolid 0
0 0 Psolid
(6.5)
74
where 0 represents the zero matrix of size ((1 + xe/2l~l) (1 + ye/2l~1) by
(1 +xe/2l) (1 + ye/21)). The interface is always on the mesh boundary of
the coarse mesh. Given the finest level matrix Aq from the finite element
application, the coarse matrices are created variationally from the recurrence
A+i < (Pl+i)TA,Pl+v One iteration of the standard Vcycle multigrid algo
rithm x < MG(x, 6), solving A0x = b, can be described in abstract terms as
follows. Set MG = MG;, where l = 1,..., L 1
Algorithm 6.1
Presmoothing: Perform mi smoothing steps on AiUi = fi
Coarse grid correction
compute residual ri = f\ Atui
restrict ft = (P/+1)Tn
if 1+1 = L then solve coarse grid equation by a direct method, else perform
mc iterations of ui+1 < MGi+1(0, fi+1)
interpolate Vi = PiVi+1
correct the solution on level l by ui < ui + Vi
Postsmoothing: Perform m2 smoothing steps on AiUi = fi
Here, mi and m2 represent the number of presmoothing and post smoothing
steps per multigrid level, and mc is the multigrid cycle parameter.
We observe that the key to satisfactory multigrid performance is the
design and choice of suitable smoothing, which is the central component of
a multigrid algorithm. We test various smoothers on a squared diagonally
75
preconditioned system. A brief description of the iterative methods used as
smoothers follows.
6.1.1 GMRES. Generalized Minimum Residual method is a pop
ular Krylov method for symmetric and unsymmetric problems [1]. GMRES is
a projection method that approximates the solution of a system of linear equa
tions A u = / for it, by minimizing the residual norm / Au\\ over the
rnth Krylov subspace {xo +Span{r,Ar,..., Am_1r}} [7, 56, 57]. The approxi
mation can be obtained as xm = x0 + Vmym, where ym minimizes the function
J(y) = 6 A{xo + Vmy)\\2. The minimizer is inexpensive to compute, since it
requires the solution of an (m + 1) x m least squares problem, where m is typ
ically small [56]. GMRES requires more memory than other Krylov methods
because it stores all search directions.
6.1.2 BICGSTAB BICGSTAB was developed to solve non
symmetric linear systems while avoiding the often irregular convergence pat
terns of the Conjugate Gradient Squared method, CGS [19, 56, 57]. The CGS
method is based on squaring the residual polynomial; however, sometimes this
is slow, and it may lead to substantial buildup of rounding errors or even over
flow. CGS guarantees that the norm of the residual decreases. BICGSTAB is
a variation of CGS, and was developed as a stabilization of the Arnoldi process.
It uses a 3term recursion, but a decrease of the norm of the residual is not
guaranteed [56]. BICGSTAB produces iterates whose residual vectors are of
76
the form
rj = ^(A)^(A)r0,
where is the residual polynomial associated with the Biconjugate Gradient
algorithm, and ipj(t) is a new polynomial that is defined recursively, so that it
stabilizes or smoothes the convergence behavior of the original algorithm.
6.1.3 GaussSeidel GaussSeidel is a popular multigrid
smoother. We have used GaussSeidel in the form Xk+i = Xk + Axk),
where the matrix Q is the lower triangular part of A that was found by using
the MATLAB function tril(A).
We present results and convergence factors for GaussSeidel, GMRES,
and BICGSTAB used as smoothers. Multigrid is often used as a preconditioner
for an iterative method. In our numerical experiments, we use GMRES pre
conditioned with one full multigrid iteration. Tests are run for different mesh
sizes h; when the mesh is refined the frequency is increased, so that the value
of k3h2 is kept constant. This is needed to avoid pollution by the phase error
and to keep the error decreasing with h for the Helmholtz equation [37]. The
essential ingredients of the multigrid algorithm are as follows. In all the prob
lems below the standard Vcycle algorithm described in Algorithm 6.1 is used
for the multigrid method. At each multigrid level the number of presmoothing
steps always equals the number of postsmoothing steps. The problem on the
coarsest level is solved directly.
77
7. Numerical Results
7.1 Numerical Verification of the Discretization
We now present preliminary numerical experiments from a prototype
implementation in MATLAB. We consider a 2D model problem as in Fig. 4.2.
The domain is the square (0,1) x (0,1). The obstacle in the channel is set up
in the center of the fluid domain as a square of size 0.2 m, unless otherwise
specified. In all computational examples, the size of the domain in the x and
ydirection is chosen to be 1 m. The boundary condition on is 'Po(x, y) = 1.
The origin of the coordinate system is assumed to be in the lower left corner
of Q. The fluid medium is water with density pj = 1000 kgrn~3 and speed
of sound Cf = 1500 ms1. The elastic medium is aluminum with density
pe = 2700 kgm~3 and Lame elasticity coefficients axe A = 5.5263 1010 Nm~2,
and p = 2.595 1010 Nm~2. The speed of sound is Cp = 6300 m s1 for pressure
waves and cs = 3100 ms1 for shear waves. The fluid domain and the solid
domain are discretized by standard bilinear square finite elements Ql on a
uniform mesh with mesh size h.
Problem (4.42) is solved directly using sparse LU decomposition in
MATLAB with SYMRCM reordering to reduce the profile bandwidth of the
matrix [52]. Table 7.1 demonstrates the convergence of the solution. The wave
78
h=To h
1.02134046416 0.85252801718i 1.02134046416 0.85252801718i 0.25394900255 + 1.07327043738i 0.25394900255 + 1.07327043738i 1.51245645214 + 2.40851542272i 0.00000000000 + O.OOOOOOOOOOOi 1.04992768430 0.83392861 lOli 1.04992768430 0.83392861 lOli 0.30531713828 + 1.06650378005i 0.30531713828 + 1.06650378005i 1.48018801699 + 2.54931071269i 0.00000000000 O.OOOOOOOOOOOi h = 11 80
1.05757604357 0.82882031388i 1.05757604357 0.82882031388i 0.31891339863 + 1.06408783048i 0.31891339863 + 1.06408783048i 1.47408877517 + 2.59422086989i 0.00000000000 + O.OOOOOOOOOOOi n 160 1.05965748637187 0.82752600255i 1.05965748637186 0.82752600255i 0.32248585398500 + 1.06339927283i 0.32248585398500 + 1.06339927283i 1.47404124589866 + 2.60934086845i 0.00000000000001 + O.OOOOOOOOOOOi h = 11 320
1.06024737503 0.82720711572i 1.06024737503 0.82720711572i 0.32354182387 + 1.06320783223i 0.32345282387 + 1.06320783223i 1.47467790973 + 2.61466232513i 0.00000000000 O.OOOOOOOOOOOi 1.06024377503023 0.82720711572i 1.06024377503010 0.82720711572i 0.32343982387092 + 1.06320783223i 0.32343982387133 + 1.06320783223i 1.47466790973365 + 2.61466232513i 0.00000000000005 O.OOOOOOOOOOOi
Table 7.1. Values of the solution at the 6 sample points as explained in
(Fig.7.1), mesh size h is halved at each run and wave number is kept constant
k = 5.
79
w 2 1 f 4
rd r. a*66 n
f 1 A 3 d k
wv r Figure 7.1. Configuration of sample points for Table 7.1
1,2,3,4 are fluid pressure values, and sample points 5 and
displacements, respectively.
Sample points
are ux and uy
80
Error decreases as 0(h2)
h Xh Eft logh logxft 
1/10 0.25394900255948 6.9457e02 2.30258 2.66704
1/20 0.30531713828081 1.8089e02 2.99573 4.01246
1/40 0.31891339863201 4.4924e03 3.68887 5.40536
1/80 0.32248585398500 9.1997e04 4.38202 6.99116
1/160 0.32354182387092 1.6600e04 5.07517 8.70352
1/320 0.32343982387092 3.4000e05 5.76832 10.28915
Table 7.2. Exact solution is z* and solution at mesh size h is x^.
number is chosen to be k = 5. The results are values at sample points for
different mesh sizes, from 10 x 10 to 320 x 320. Configuration of sample points
for Table 7.1 are displayed in Fig. 7.1. Sample points 1,2,3,4 are fluid pressure
values, and sample points 5,6 are ux and uy displacements, respectively. Denote
the unknown exact solution by Â£* and the solution at mesh size h by Xh Then,
Xh x* Ch2 => loglz/j x*\ log/i+const
To find rr* by extrapolation, we solved the system of equations
x* + Ch? xh
x* + C(2h)2 = x2 h
We use as an example solution points 3 in the lower left corner from Fig. 7.1
to show that the solution converges as 0 (h2).
We display the solution for some cases considered. We display the
real part of the pressure amplitude as a function of elevation above the x y
plane, and the displacement in the x and the y directions as equally spaced
arrows based at the original configuration. The solution obtained by solving
81
Figure 7.2. Log log plot of difference of solution values Xh at sample points
for mesh size 10 x 10 to 320 x 320 and extrapolated exact solution a:*, k=5.
82
fluid pressure
solid displacement
Figure 7.3. Solution for a 40 x 40 mesh, k = 5, righthand side modified for
the Dirichlet boundary condition p = p0 on
83
fluid pressure
Figure 7.4. Exact solution for a 64 x 64 mesh, k = 15, righthand side is
modified for the Dirichlet boundary condition p = p0 on Td
84
fluid pressure
Figure 7.5. Three multigrid iterations, 20 smoothing steps, smoother GM
RES, 2 levels to solve a 64 x 64 mesh, k = 15, righthand side is modified for
the Dirichlet boundary condition p = po on Td.
85
