Large-Scale Structures in Acoustics and Electromagnetics: Proceedings of a Symposium (1996)

J. Tinsley Oden

Andrzej Safjan

Po Geng

Leszek Demkowicz

University of Texas at Austin

Much of structural acoustics is concerned with the classical problem of interaction of the motion of an elastic structure with that of a perfect inviscid fluid into which the structure is immersed. Small perturbations in the velocity or pressure fields of the fluid create waves that impinge upon the structure and scatter into the fluid-structure domain while small motions of the structure may radiate energy into

the surrounding fluid media. The radiation and scattering of waves in such fluid-structure interaction has been the subject of research for many decades.

Mathematically, the problem can be adequately modeled by a system of linear conservation laws, which, with further manipulations, can be reduced to the wave equation for the acoustical pressure coupled to the equations of elastodynamics for the structure. Traditionally, approaches for the computer simulations of structural acoustics phenomena have been primarily based on frequency-domain formulations of such problems, in which a Fourier transformation of the dynamic equations is made, and this results in a system of partial-differential equations in complex-valued field variables with frequency-dependent coefficients. The fluid may be modeled with the exterior Helmholtz equation or, as is commonly done, by an equivalent boundary-integral equation. These classical approaches have been very successful for a large and important class of problems, particularly those involving scales which are not too small in relation to characteristic dimensions of the structure. In recent times, attempts at large-scale simulations involving high ka values (k being the wave number and a the characteristic dimension) have encountered formidable computational difficulties and have exposed serious and fundamental limitations of some classical frequency-domain approaches.

In the present work, some new high-order time-domain methods are presented that have proven to be effective for simulating certain classes of problems in structural acoustics and that may have potential for overcoming many of the shortcomings of the frequency-domain methods for high ka. These techniques are built around several special ideas and algorithms:

1. the structural acoustics problem is formulated as a hyperbolic system of conservation laws and this leads to an abstract Cauchy problem in common Hilbert space settings in which the spectral theory of linear operators is applicable;

2. new, high-order, multistage Taylor-Galerkin (TG) approximations of these equations in the time-domain are constructed which deliver high-order temporal accuracy with unconditional stability;

3. spatial approximations are constructed using hp-finite element methods in which local mesh-size h and local spectral order p are varied in a way to adaptively control error, which provides multi-level approximations of the unknowns;

4. a posteriori error estimates are constructed over each element and time step to provide a running account of solution quality and to provide a basis for the adaptive control strategy;

5. applications to representative model problems reveal that the methods are robust and provide remarkably accurate results for a class of two-dimensional cases.

In addition to the time-domain methods, a brief account of a coupled boundary element/finite element method is also given. This approach also employs hp-methods and is based on a new version of the Burton-Miller integral equations for the acoustical fluid. Results obtained by using large-scale parallel computers for such problems are presented.

We also make comparisons of time-domain and frequency-domain approaches for a model class of problems. In these computations, some advantages of the time-domain approach over traditional frequency-domain approaches are observed.

Our investigation focuses on the simulation of structural acoustics phenomena using classical models of acoustics and elastic bodies characterized by systems of hyperbolic equations. The elastic structure is modeled by the classical Navier-Lamé equations of linear elasticity, and the fluid is a perfect acoustical medium modeled as a quiescent, inviscid fluid. The mathematical setting is described as follows, and is depicted in Figure 1.1.

Figure 1.1 The solid-structure interaction problem.

Let be a region occupied by an elastic body and and let its boundary consist of three distinct parts, where Г_u and Г_t are portions of with prescribed displacements and tractions, respectively, and Г_sf is the the portion of subjected to the action of the fluid (the contact surface). Similarly, let be a region occupied by acoustic fluid and let the boundary consist of four distinct parts, , where Г_v and Г_p are portions of with prescribed velocities and pressure, respectively, denotes the truncated exterior boundary, and Г_fs denotes the part of the boundary which is in contact with the elastic body. We further denote Г_fs = Г_sf = Г_c' and we use superscripts s and f to distinguish between variables having Ω^(s) and Ω^(f) as domains.

The fluid-solid interaction problem consists of solving:

• Equations of linear acoustics:

• Equations of linear elastodynamics:

• Contact conditions:

Here:

= u^(s) (x,t) is the displacement vector at particle at time t

s(u^(s)) = (s_ij) = CDu^(s) = the stress tensor evaluated at the displacement u^(s)

ρ_s = mass density of the elastic body

C = 6 × 6 symmetric positive definite matrix of elastic constants

n = (n_i) = the unit outward normal (at the interface, it is the unit outward normal to the elastic body)

p = p(x,t) is the acoustic pressure at at time t

 = v^(f) (x,t) is the velocity in the acoustic fluid

ρ_f = mass density of the acoustic fluid

c₀ = small signal sound speed in the acoustic fluid

= boundary operator at the truncation boundary

D = generalized gradient operator defined below

Problem (1.1)-(1.2) is completed by specifying initial conditions:

Here are given boundary and initial data.

Equations (1.1) and (1.2) are well known, so we comment only on the interface conditions (1.3). On Г_c, we have continuity of the normal component of velocity in both media, and continuity of the normal tractions. In addition, since the acoustic fluid is inviscid, we must have zero tangential traction on the elastic body on Г_c.

We notice that the pressure p in (1.1) is a primitive variable which is coupled with stresses s of (1.2) through contact conditions, and the stress is not a primitive variable in (1.2). These incompatibilities prompt us to seek alternative formulations of the problem.

The formulation is similar to that reported in Hubert and Palencia (1989). The principal steps in the proposed methodology are as follows.

1. Introducing an auxiliary variable u^(f),

and assuming the compatibility condition

equations (1.1) are transformed to the following form:

In addition, we make the usual assumption that u^(f) is a potential vector field:

for some scalar field Φ.

Finally, the transition conditions (1.3) take the form

and the fluid-solid interaction problems amount to solving equations (1.2), (1.8), (1.9), (1.10).

2. Weak formulation. We now restrict attention to the case in which and , as indicated in Figure 1.2. Then Ω^(f) and Ω^(s) are assumed to be disjoint open connected sets with smooth boundaries which intersect in a C² manifold Г_c of dimension N - 1. Let u denote the displacement vector in the whole region , and u^(f) and u^(s) denote, respectively, its restrictions to Ω^(f) and Ω^(s).

We start by introducing the spaces of kinematically admissible displacements

Figure 1.2 A submerged elastic structure.

Next, we define the following Hermitian forms

The space X, when endowed with the following inner product,

is a Hilbert space. In addition, we define space H as the completion of X with inner product (u, w)_H, and, finally, we define space D(Ã) as follows

Now, we consider the following weak form of the fluid-solid interaction problem (1.2), (1.8), (1.9), (1.10):

Find

such that

and

It is easily seen that (1.19) is a classical virtual work formulation of (1.2), (1.8), (1.9), (1.10). In fact, (1.19) is valid for any test functions of appropriate regularity, not only for those of the form of gradient.

Let , where denotes the space of test functions. Then (1.18) implies that (1.2) and (1.8) are satisfied in a distributional sense. Next, define a as

Green's formula applied to Ω(^l) and Ω(^s) gives

which gives transition conditions (1.10). Here n = (n_i) is the unit outward normal to the elastic body and τ = (τ_i) is the corresponding unit outward tangent.

3. Within the above formalism, the initial boundary value problem of structural acoustics can be reinterpreted as an abstract second order Cauchy problem:

where à satisfies

u = u(t) is an H- valued function of time, and u₀ and v₀ are initial conditions.

To proceed further, we reinterpret (1.23) as a first-order system in time. To this end, we introduce:

• the group variable

• the Hilbert space

• the operator A

Then, (1.22) becomes:

Here U₀ specifies initial conditions, .

To show the well-posedness of (1.28) we consider an operator

and we show that

(i) is accretive, i.e.,

and,

(ii) is surjective, i.e., the range of is the whole space

Indeed,

and, hence,

To prove that is surjective, we need to show that for each there exists a solution of (A + 2I) = F, or equivalently, of the following system

System (1.34) can be reduced to the following equation

which, by virtue of the Lax-Milgram Theorem, possesses a unique solution . This, together with (1.35) and the definition of D(Ã), implies that . Then (1.34)₂ implies that .

Summarizing, is accretive, is surjective, and D(A) is dense in , which implies that the Cauchy problem for operator is well posed (cf., the Lummer-Philips Theorem in Yosida, 1973). Finally, by using the ansatz

we get the well posedness of (1.28). Thus, given , there exists a unique such that (1.28) is satisfied. Moreover,

Problem (1.28) is to be solved numerically by using high-order Taylor-Galerkin methods and hp-adaptive finite element methods.

The second approach consists of converting the system of equations of elastodynamics into an equivalent first-order system, and coupling it with (1.1). Thus, we solve the following first-order system. (We label it the ''velocity-stress'' formulation.)

• equations of linear acoustics:

• equations of elastodynamics:

• contact conditions:

together with appropriate boundary and initial conditions. Here: , and s_ij are the components of the stress tensors.

As before, (1.38), (1.39), (1.40) can be cast in the form of a Cauchy problem

where U is the group variable U = (v^T,w^T)^T, i is the imaginary unit, and A is an appropriately defined operator. Again, (1.41) is solved numerically by using high-order Taylor-Galerkin methods and hp-adaptive finite element methods.

The Taylor-Galerkin (TG) schemes represent Galerkin approximations of temporal Taylor expansions of the solution. The first TG scheme was presented by Oden (1974) in connection with a Lax-Wendroff approximation of nonlinear waves by finite element methods. The TG schemes were made popular by a series of important papers by Donea and collaborators (Donea, 1984; Selmin et al., 1985), but were not applied to schemes of higher "order" than three. By a scheme of "order s" we shall understand that the one-step truncation error measured in L² -norm is bounded by CΔt^s+1, where At is a time-step size and C is a constant independent of Δt. The idea of developing multi-stage high-order TG schemes was presented in the dissertation by Safjan (1993) and in subsequent papers by Safjan and Oden (1993, 1994, 1995a,b). A brief summary of the algorithm is presented below.

Let H be a Hilbert space with inner product (·,·) and corresponding norm , and be a self-adjoint operator in H. Being self-adjoint, A admits the spectral decomposition of the form (Yosida, 1973)

where E_λ is a uniquely defined spectral family. In this section we introduce TG schemes for solving abstract Cauchy problems of the form

Any approximation to (1.44) must involve discretization both in space and time variables. We will adopt the assumption here that the final approximation is obtained by using finite differences in time and

finite elements in space variables. Then, two different approaches are possible. In the classical method of lines, an approximation in space variables converts the original initial-boundary-value problem into a system of ordinary differential equations (ODEs), which next is discretized in time using one of many time integration schemes for ODE.

The TG schemes belong to a different class of methods, known as themethod of discretization in time, which consists of the same two steps but done in the reverse order. By discretizing in time first, the initial-boundary-value problem (1.4) is converted into a sequence of boundary-value problems (1.45)

which, in turn, give a basis for a spatial approximation and, consequently, a fully discretized scheme,

Here U_τ and U_τh are the solutions of the semi-discrete problem (1.45) and the fully discrete problem (1.46), respectively, T is an appropriate transient operator, T_h denotes its finite-dimensional realization, and U_0h is a suitable approximation of U₀.

Let N be a positive integer. We define a partition in time as Δt= t */N, t_n = nΔt, n = 0, 1,..., N and consider a typical time-step . Given the solution U_τ (t_n), we seek the next time step solution U_τ(t_n + Δt) through an s-stage scheme of the following form

where

Coefficients µ_ij, v_ij, µ_i0, v_io, c_i, are to be chosen so as to obtain the highest possible order of accuracy, subject to stability or other constraints. A free parameter η is to be chosen from stability considerations.

It is convenient to rewrite (1.47) in the following compact form:

where

Matrices N,M, and K together with vector c completely characterize the difference scheme (1.49).

A distinct feature of (1.49) is that coefficient matrices N and M are lower triangular which makes the resulting scheme semi-implicit (i.e., to compute Z_i it is necessary to know Z_i-1,Z_i-2,...,Z₁, but it is not necessary to know Z_i+1 ). Moreover, all diagonal elements of N are equal (v_ii = η, i = 1, 2,..., s) and so are those of M (μ_ii = 0, i = 1, 2,..., s). This makes the operator defining the left-hand side of each stage of (1.46) identical for linear (or linearized) problems, and, hence, significantly reduces the cost of applying the method. A particular choice of zero diagonal elements of M is made with an eye on well-posedness of a typical one-stage problem and a possible splitting of the operator defining the left-hand side of each stage.

To make the i-th stage solution Z_im-th order accurate, it is necessary to satisfy the order conditions for Z_i and to make the previous stage solutions Z_i-1,Z_i-2,...,Z₁, at least of the order m-1, (otherwise, some coefficients have to be set to 0; e.g., if Z_i-1 is of the order m-2, then, necessarily μ_i,_ι-1 = 0 ). The order conditions for Z_i are obtained by expanding it in Taylor series about Z₀:

Introducing this expression into the left-hand side of the i-th stage equation,

and equating coefficients of like powers of At to zero, leads to the following system of nonlinear algebraic equations:

Equations (1.53) are referred to herein as the order conditions.

It is important to realize that the particular forms of the coefficient matrices M and N limit the attainable order of (direct) TG schemes. For example, it is not possible to attain order 6 in 3 stages. To remedy this situation, a more general class of transformed Taylor-Galerkin schemes has been introduced in Safjan and Oden (1994, 1995a).

First, using the original equations (1.44), we calculate the time derivatives in terms of spatial derivatives as follows:

Next, replacing the time derivatives in (1.49) by formulas (1.54) and (1.55), we arrive at the following system of equations,

where the indices i and j are used to denote a particular stage. Finally, multiplying (1.56) by a vector-valued test function W, integrating over Ω and integrating by parts, we arrive at a variational formulation of the form:

where X = D(A), and the bilinear (sesquilinear) forms A, B, and C are defined by

Weak form (1.57) is the basis for finite element approximations.

The properties of TG schemes are summarized as follows:

1. the TG methods lead to a series of well-posed problems that can be solved using hp-adaptive finite element methods with very high accuracy and high rates of convergence,

2. the TG methods have excellent accuracy; some of them possess a built-in temporal error estimate,

3. with a proper choice of parameters, TG methods can be unconditionally stable for certain classes of problems,

4. the TG methods can be classified as semi-implicit schemes, so that they are generally more efficient than fully implicit schemes, and provide a natural splitting of the operator that preserves accuracy, and

5. in most cases, the resulting system of linear algebraic equations is symmetric and positive definite, which is a comfortable setting for most iterative solvers.

For more details on TG schemes, we refer to Safjan and Oden (1993, 1994, 1995a,b).

In many situations, we wish to solve a transient problem in an infinite domain. However, any computational domain is necessarily bounded, so it is essential to introduce appropriate boundary conditions at the boundaries. In particular, we are interested in radiating boundary conditions that have the following features:

1. the boundary conditions substantially reduce the (unphysical) reflections from the artificial boundaries,

2. the boundary conditions are local,

3. the boundary conditions together with governing differential equation result in a well-posedproblem,

4. the boundary conditions fit into the framework of our Taylor-Galerkin schemes.

We first develop radiating boundary conditions for the scalar wave equation and then deduce boundary conditions suitable for the equivalent 1st-order hyperbolic system. Our development follows Enquist and Majda (1977, 1979).

Let p denote the acoustical pressure and let us consider the following initial boundary-value problem for the wave equation

Here g(y, t) is an arbitrary square-integrable source function and H(t) is the Heaviside step function. Let δ be an acceptable error and suppose that we are interested in the solution of (1.59) in the strip , rather than in the whole right half-plane. Our objective is to find a boundary operator, , on the wall x = a with , but as close to a* as possible, so that if is the solution of

then

where time T is so large that M reflections from the artificial boundary at x = a might have occurred.

We first derive a theoretically exact nonreflecting boundary condition. By taking the Fourier transform with respect to y and the Laplace transform with respect to t, we arrive at the following solution of (1.59):

where

ξ and ω to are dual variables to y and t, respectively, and ^ denotes the Fourier transform. We notice the correspondences and .

Let us introduce the operator which is non-local in space and time, defined by

Then, a simple computation reveals that the solution p satisfies the following non-local theoretical radiation condition at x = a:

By taking the Fourier transform on (1.65) we arrive at

Introducing the approximation,

(1.66) reads

By taking the inverse Fourier transform on (1.68), we get the first radiating boundary condition :

Using the next approximation to (1 - ξ² / ω²)^1/2:

which, after taking the inverse Fourier transform, yields the boundary operator :

To get stable boundary operators of arbitrary order, we consider the Padé approximation to (1 - ξ² / ω² )^1/2 of the form

Here P_N (z) are recursively defined by

and T_N+1 and S_N+1 are even polynomials in z, satisfying the following relations:

Introducing (1.72) into (1.66) we get

By clearing the denominator we obtain

and, observing from (1.74) that deg and deg , we get

Using identities (1.74), it may be verified that satisfies

Finally, defining via the correspondence and , we get

or

We now record the following fundamental result (Enquist and Majda, 1977).

1. Every is a local strongly well-posed boundary condition for the wave equation.

2. Given δ, an acceptable error, a* > 0, an integer M, and an arbitrary source function satisfying

there is an and N so that if solves the boundary value problem in (1.60) with then

for any time T with .

1. Operator is exact for a single space dimension and can easily be derived from physical considerations.

2. The distance a dictating the size of computational domain, and the order of approximation N are critical parameters in practical computations.

3. Using a Taylor expansion of (1 - ξ² / ω² )^1/2 rather than Padé approximations may result in strongly unstable boundary conditions (Enquist and Majda, 1979).

We first notice that wave equation (1.59) can be derived from the 1st-order system

and vice versa. Here u and v denote the x- and the y-component of velocity vector, respectively, and we notice that (1.83) is nothing but the system of equations of linear acoustics in appropriately defined nondimensional variables. Then the boundary operator ,

can be written as

Next, assuming that the acoustic field is quiescent in the vicinity of the artificial boundary at time t = 0,

we integrate (1.85) with respect to time and obtain the following boundary condition:

Similarly, boundary operator for system (1.83) reads

Finally, we notice that boundary operator is not easily implementable. Indeed, using recurrence relation (1.80) we compue as

Next, introducing (1.83) into (1.89) and integrating with respect to time we get

Since (1.90) contains spatial derivatives of order 2, it is not acceptable for a second order operator. A new class of split boundary operators circumvents this difficulty (Safjan and Oden, 1995b).

In this section we give a brief outline of the hp-finite element schemes used in construction of approximation spaces .

First, we introduce a partition of into N_ELEM subdomains (elements) with disjoint interiors

and each element K is the image of a master element under an invertible map F_K,

On we define the following three categories of shape functions:

• Node functions. These are bilinear:

• Edge functions. Defining nodes at the midpoints of each side of we define the edge functions

Here is associated with side of , Δ = 1, 2, 3, 4and Φ_k are integrated Legendre polynomials; e.g.,

P_k-1 (s) = the Legendre polynomial of degree k - 1. Thus, polynomials of differing degree k_Δ can be assigned to each edge .

• Bubble functions. On the interior of , at its centroid, we assign the shape functions,

Again, different orders can be assigned to different directions at (0, 0).

The master element shape functions define a space,

which, together with (1.92), define P_K, a space of local shape functions for element K:

We now define X_h,p, a global finite element approximation space:

The space X_h,p possesses standard interpolation properties of hp-finite element methods (see, e.g., Babuška and Suri, 1987). In particular, for element K which is affine equivalent to a master element (i.e., T_K is an affine map), if h_K = dia (K) and p_K is the highest degree of the complete polynomials contained in the span of shape functions associated with this element, then, for any function , an interpolant exists, such that

where P_pK(K)= space of polynomials of degree defined on K, constant C is independent of u, p_k or h_k, and denotes the usual Sobolev norm. Then, globally, for uniform h and p, and for , there exists an interpolant , such that

where, again, C is independent of u, p or h.

Finally, we define as

As mentioned above, (1.94) and (1.95) allow us to use shape functions of different orders for different edges (and different directions), and, hence, to construct a global mesh with a non-uniform distribution of spectral orders p. In addition, we treat mesh size h as a parameter by using h-adaptive techniques for bisecting elements and imposing hp-constraints to maintain continuity across element edges. In particular, we employ unstructured, anisotropic, 1-irregular meshes. For the mesh is composed of quadrilateral elements which are refined by dividing appropriate elements into two siblings, as opposed to four, in such a way that only one constrained node exists between any two active vertices. This provides for directional refinement whenever refinement in one direction is needed but none is called for in an orthogonal direction. An h-refined mesh of this type is depicted in Figure 1.3.

Continuity of global basis functions is maintained by imposing constraints at the interface of elements supporting edge functions of different order. Thus, for example, one encounters cases of the type illustrated in Figure 1.4 where an isotropic h-refinement has resulted in two small elements interfacing a large element with p-shape functions of differing order assigned to the edge (and interior) nodes. In the figure, open circles denote active vertex nodes which support bilinear shape functions, open rectangles denote edge and interior nodes which support various high-order shape functions, and darkened circles and rectangles denote constrained nodes. Thus, along the interface, values of the order of p at the p-nodes of the smaller elements (in this example) are constrained to exactly fit the shape function of the larger element to produce full continuity across this interface. For details on hp-constraints, see Demkowicz et al. (1989, 1991 a,b).

Figure 1.3 Unstructured, 1-irregular mesh of quadrilateral elements with anisotropic h-refinements. Different shadings correspond to different element spectral orders p.

Figure 1.4 Mesh topology in two dimensions showing constrained h- and p-nodes.

For completeness, we briefly summarize parallel work we have done on the analysis of structural acoustic problems using new coupled boundary-element and finite-element approaches in which the problem is formulated in the frequency domain.

Frequency-domain formulations are obtained by Fourier transformations of the time-domain equations and variables. For a function , its Fourier transform in time is a function û(x,ω) defined by

where ω is the frequency. We recover u(x, t) from û(x, t) by the inverse transform

Denoting the transformed displacement and pressure fields (amplitudes at the frequency co) in (1.38) and (1.39) by u(x) and p(x), respectively, we observe that a frequency-domain formulation of dynamical problem (1.27) is of the form

and

Also, the following boundary conditions must be enforced:

where k = ω / c is the wave number, Г_c is the interface boundary between the interior domain Ω(^s) and the exterior domain Г(^f), and σ_n,στ, and p^inc represent amplitudes of normal stress, tangential stress, and incident wave pressure at the frequency ω, respectively.

A major advantage of solving problems in the frequency-domain is that the Helmholtz equation on the exterior domain Ω(^f) can be converted into an integral equation on the boundary Г_c, i.e.,

In this way, the dimensionality of the problem on Ω(^f) can be reduced by 1 and the infinite boundary condition (1.108) is satisfied automatically. In (1.107), Φ represents the fundamental solution of the Helmholtz equation and is of the form

and

where is the Hankel function of the first kind of 0 order and .

As is well known, a major problem in using integral forms of the Helmholtz equation on exterior domains is that its solutions are nonunique at certain frequencies (forbidden frequencies). The CHIEF (Combined Helmholtz Integral Equation Formulation) method, proposed by Schenck (1967), is one of the earliest and simplest methods to overcome the nonuniqueness problem. In this method, the system of linear equations is combined with a few additional equations generated from the Helmholtz integral equation with the source points in the interior of the closed boundary. However, the CHIEF method is not reliable or effective in problems with complicated boundary geometry and frequencies in the higher range.

Burton and Miller (1971) proved that a complex linear combination of (1.109) and a so-called hypersingular integral equation,

is equivalent to the original differential equation. Replacing (1.106) and (1.108) by the Burton-Miller formulation leads to the following problem

Find amplitudeu = u(x)and p = p(x) such that

and

where α is a complex constant whose imaginary part cannot be 0.

Equations (1.111) through (1.113) define the mathematical model which we use to solve the acoustic scattering problem at a given frequency ω. It can be proven that the problem defined in (1.111-1.113) is strongly elliptic on the space of , and at non-resonant frequency the problem is well defined.

A special problem that arises in the implementation of the Burton-Miller formulation is that the integral

in (1.112) cannot be evaluated easily in a meaningful way. Symbolically, it is customary to write (1.114) in the form

The integral kernel contains a hypersingular term and it easily verifies that as in 2 dimensions

and in 3 dimensions

Although the hypersingular integrals can be computed numerically in the sense of the finite part integral (see Kutt, 1975 and Demkowicz et al., 1991a for details), the implementation for such procedures is complicated and the accuracy of its result is usually very poor.

The method used in this study to avoid hypersingular integral kernels associated with (1.115) was proposed in Karafiat et al. (1993) and Demkowicz et al. (1992). The method requires that the problem be formulated in a weak (Galerkin) approach and is based on the transformation

and, in 2 dimensions, this reduces to

where q(x) is the test function and t(·) is the positive tangential direction at the corresponding position of x and y.

With (1.116) or (1.117), the variational formulation of (1.116)-(1.118) can be written as:

Findsuch that for any ,

and

or

where .

Equations (1.119) and (1.120) contain only weak singular and regular items, which can be treated easily in numerical implementations; and the whole problem can be solved by the coupled finite element and boundary element methods.

The coupled boundary element and finite element codes in this research are developed in a special hp-hierarchical data structure proposed in Demkowicz et al. (1991b), and, as mentioned, this data structure allows the mesh size and order of polynomial interpolation to be varied over the mesh and provides a good basis for the development of an adaptive method. Comparisons with our time-domain approach are given in the next section for the elastic and rigid-scattering problems, respectively.

We consider the problem of diffraction of a plane wave consisting of a short train of pulses, by an elastic shell, see Figure 1.5. The problem was solved for the following data:

• Small signal sound speed, c₀ = 1500 m / s

• Ambient mass density, P_f = 7700 kg / m³

Figure 1.5 Plane progressive wave impinging on a cylindrical elastic shell. Problem definition.

• Young modulus, E = 19.5 × 10¹⁰Pa,

• Poisson's ratio, v = 0.28,

• Mass density, Ps = 7700 kg / m³

The initial condition function:

where , and A = 0.005.

The problem was treated as a 1st-order system, and was solved by a 2-stage 4th-order Taylor-Galerkin method with stability parameter η = 0.471 and constant time step Δt = 0.01. The solution algorithm described by Safjan and Oden (1993) was used. As the solution of the problem is fairly smooth, we chose to use p-enrichments/unenrichments only. The only exception from this rule was made in elements with contribution to the global error exceeding 50 percent; those elements were first h-refined. In addition, only p- or anisotropic h-refinements were allowed at the interface.

Figures 1.6 through 1.16 show the computed spatial waveforms of various acoustical and elastic variables in nondimensional form. Figure 1.6 shows the initial finite element mesh and Figure 1.7 shows the initial condition function in the form of contour maps. Figures 1.8, 1.9, and 1.10 show pressure

distribution at time t = 0.5, 1.0, and 1.5, respectively. A forward radiation from the shell can be observed, as the wave propagation speed in the shell is approximately 4 times greater than that in acoustic fluid. Figures 1.11-1.13 show the propagation of θ - θ (bending) stress waves in the elastic shell. Finally, Figures 1.14-1.16 show the evolution of finite element meshes.

The problem is defined as follows (see Figure 1.17):

1. Governing equations of linear acoustics are to be solved in the domain

   
   
   
   

2. Boundary condition at r = a,

   

   where A = -0.005, ω = 2π, T = 1 and H(•) is the Heaviside step function.

3. Sommerfeld radiation condition at

4. Initial condition

   

5. Unit material constants: c₀ = 1 and P_f = 1.

For computations we accept as the computational domain and we use the second highly absorbing boundary condition , on truncation boundary . The problem was solved using a 2-stage 4th-order Taylor-Galerkin method with stability parameter η = 0.471, constant time step Δt = 0.01, and a fixed finite element mesh. Figure 1.18 shows the finite element mesh, and the consecutive Figures 1.19, 1.20, 1.21, and 1.22 present the computed pressure distributions at time stages t = 1, 1.5, 2.0, and 2.5, respectively. As can be seen, the wave is leaving cleanly the computational domain. For comparison, the same experiment is performed except that the first highly absorbing boundary condition is used on truncation boundary . Figure 1.23 presents the computed pressure distribution at time t = 2.5. From a comparison of the performance of with that of , it may be inferred that represents a significant improvement over .

The Galerkin-Burton-Miller formulation developed here is intended to be used to solve problems at very high frequencies. To verify this, we consider an example of a harmonic plane incident wave

impinging on a cylinder with a boundary condition on the surface of the cylinder (as shown in Figure 1.24). The analytic solution is available in the case of β being a constant (Demkowicz et al., 1991b). Clearly, for β = 0, we have the problem of acoustical wave scattering on a rigid body. Figure 1.25 compares the numerical solution and corresponding analytic solution for the pressure distribution around the surface of the cylinder with ka = 400 and . With a uniform mesh of 1,024 quadratic elements and 2,048 degrees of freedom, the result illustrated in Figure 1.25 shows a good correlation between the numerical solutions and analytic solution and the L² -norm of relative error, computed by

is around 2 percent in this example where P_e represents the exact solution and p^h is the numerical solution.

Figure 1.24 The test problem: a plane harmonic wave impinging on a cylinder.

Equations (1.105)-(1.108) define a problem of an elastic structure surrounded by an inviscid fluid field. Because of the absence of damping, the problem is not uniquely solvable at certain frequencies (resonant frequencies) and the stiffness matrix becomes singular or ill-conditioned when the frequency ω is at or close to one of these resonant frequencies.

Figure 1.25 The comparison between numerical solution and analytic solution around the cylindrical boundary with ka = 400 and . Because of symmetry, only the solution from 0° to 180° around the cylindrical boundary is plotted.

To verify this, the coupled elastic cylinder and fluid model given in Example 1 is solved for the same incident wave given in (1.125) and the wave number here varies from k = 0.02 to 14 with Δk = 0.02. A uniform mesh of 256 quadratic finite elements and 128 quadratic boundary elements is used in the computation. Figure 1.26 plots the relationship between the amplitude of total pressure at a far field point (1,0) and wave number k. The experimental results show that the test problem possesses several resonant frequencies in the range of . As a comparison, the analytic solution based on elastic thin shell theory (see Junger and Feit, 1972) is also given in Figure 1.26, and for the small wave number, the elastic solution and thin shell solution are close to each other.

To understand the influence of resonant frequencies on the numerical solution, the test problem given in Example 1 is solved again by the frequency-domain approach. The following incident wave is used

Figure 1.26 The resonant frequencies of the coupled elastic cylinder and fluid model. The plain harmonic incident wave is used and the problem is solved from k = 0.02 to k = 14 with Δk = 0.02.

where τ = x - ct and the time unit is 1/c. As shown in (1.126), the incident wave pulse must be supposed as a periodical function and T here must be large enough so that the solution at the time period in question is not affected by the solution wave from other periods. In this example, we are interested in the solution in the time period and T = 5.2 is selected. Then, the incident wave in (1.126) is decomposed into the sum of a series of sine functions, i.e.,

and N = 2,048. At most frequencies, p_i is in fact almost zero; here, we solve the problem only at a frequency when its corresponding amplitude satisfies

Then, instead of 2,048 different frequencies, we need to solve the problem only in 124 different frequencies.

Figures 1.27 and 1.28 present the results at t = 0 and t = 2/c for the rigid scattering (, upper half part) and elastic scattering (lower half part), respectively. For the rigid scattering, because the

Helmholtz equation with Neumann boundary condition is well defined on the exterior domain (no resonant frequencies), the frequency-domain method produces essentially the same result obtained by the time-domain method. However, the frequency-domain method does not provide a satisfactory result for elastic scattering, because the method falls at the frequencies close to the resonant frequencies of the coupled solid-fluid system. Figure 1.27 shows spurious scattering waves produced by aliasing in the frequency-domain approach. These nonexistent waves are produced ahead of the plane wave that actually excites the structure. It is likely that these types of erroneous excitations could be controlled, if viscous and material damping effects were added to the model.

Figure 1.27 The acoustic pressure contours for the rigid scattering (upper half) and elastic scattering (lower half) at t = 0. The frequency-domain approach is used and the problem is solved in 124 different frequencies. The elastic cylinder produces spurious waves before the incident plane wave excites it.

Figure 1.28 The acoustic pressure contours for the rigid scattering (upper half) and elastic scattering (lower half) at t = 2. The frequency-domain approach is used and the problem is solved in 124 different frequencies.

Support of this work by the Office of Naval Research (ONR) under #N00014-89-J-1451 is gratefully acknowledged.

Babuška, I., and M. Suri, 1987, ''The hp version of the finite element method with quasiuniform meshes,''RAIRO Math. Model. Numer. Anal.21, 199-238.

Burton, A.J., and G.F. Miller, 1971, "The application of integral equation methods to the numerical solution of some exterior boundary-value problems,"Proc. R. Soc. London Ser.A 323,201-210.

Demkowicz, L., J.T. Oden, W. Rachowicz, and O. Hardy, 1989, "Toward a universal hp adaptive finite element strategy. part 1: constrained approximation and data structure,"Comput. Meth. Appl. Mech. Eng.77, 79-112.

Demkowicz, L., J.T. Oden, M. Ainsworth, and P. Geng,1991 a, "Solution of elastic scattering problems in linear acoustics using hp boundary element methods,"J. Comput. Appl. Math.36, 29-63.

Demkowicz, L., J.T. Oden, W. Rachowicz, and O. Hardy,1991b, "An hp Taylor-Galerkin finite method for compressible Euler equations,"Comput. Meth. Appl. Mech. Eng.88, 363-396.

Demkowicz, L., A. Karafiat, and J.T. Oden,1992, "Solution of elastic scattering problems in linear acoustics using hp boundary element method,"Comput. Meth. Appl. Mech. Eng.101, 251-282.

Donea, J., 1984, "A Taylor-Galerkin method for convective transport problems,"Int. J. Numer. Meth. Eng.88, 101-120.

Enquist, B., and A. Majda,1977, "Radiation boundary conditions for the numerical simulation of waves,"Math. Comp.31, 629-651.

Enquist, B., and A. Majda, 1979, "Radiation boundary conditions for acoustic and elastic wave calculations,"Commun. Pure Appl. Math.32, 313-357.

Hubert, J.H., and E.S. Palencia, 1989, Vibration and Coupling ofContinuous Systems, Berlin, Heidelberg, New York: Springer-Verlag.

Junger, M.C., and D. Feit,1972, Sound, Structures, and Their Interaction , Cambridge, Mass.: MIT Press.

Karafiat, A., J.T. Oden, and P. Geng, 1993, "Variational formulations and hp-boundary element approximation of hypersingular integral equations for Helmholtz-exterior boundary-value problems in two dimensions,"Int. J. Eng. Sci.31, 649-672.

Kutt, H.R., 1975, "On the numerical evaluation of finite part integrals involving an algebraic singularity,"Report WISK 179, Pretoria: The National Research Institute for Mathematical Sciences.

Oden, J.T., 1974, "Formulation and application of certain primal and mixed finite element models of finite deformations of elastic bodies." In: Computing Methods in Applied Sciences and Engineering , R. Glowinski and J.L. Lions (eds.), Berlin, Heidelberg, New York: Springer-Verlag.

Safjan, A., 1993, "High-order Taylor-Galerkin and adaptive hp methods for hyperbolic systems with application to structural acoustics,"Ph.D. Dissertation, University of Texas at Austin.

Safjan, A., and J.T. Oden, 1993, "High-order Taylor-Galerkin and adaptive hp methods for second-order hyperbolic systems: application to elastodynamics,"Comput. Meth. Appl. Mech. Eng.103, 187-230.

Safjan, A., and J.T. Oden, 1994, "High-order Taylor-Galerkin and adaptive hp methods for first-order hyperbolic systems,"TICAM Report.

Safjan, A., and J.T. Oden, 1995a, "High-order Taylor-Galerkin and adaptive hp methods for hyperbolic systems," submitted to J. Comput.Phys.

Safjan, A., and J.T. Oden, 1995b, "Highly accurate nonreflecting boundary conditions for wave propagation problems," to be submitted.

Schenck, A.H., 1967, "Improved integral formation for acoustic radiation problems,"J. Acoust. Soc. Am.44(1), 41-58.

Selmin, V., J. Donea, and L. Quartapelle, 1985, "Finite element methods for nonlinear advection,"Comput. Meth. Appl. Mech. Eng.52, 817-845.

Yosida, K., 1973, Functional Analysis, Fourth edition, Berlin, Heidelberg, New York: Springer-Verlag.