Acoustic infinite elements

Transient formulation

The solution in the unbounded acoustic medium is assumed to be linear and governed by the same equations as the finite acoustic region:

\frac{1}{K_{f}} \ddot{p} + \frac{γ}{ρ_{f} K_{f}} \dot{p} - \frac{\partial}{\partial x} \cdot (\frac{1}{ρ_{f}} \frac{\partial p}{\partial x}) = 0.

$γ$ is used here to denote the volumetric drag parameter. Consider the infinite exterior of a region of acoustic fluid bounded by a convex surface and a conventional finite element mesh defined on this surface. Each facet of this surface mesh, together with the normal vectors at the nodes, defines a subdivision of the infinite exterior that will be referred to as the “infinite element.” Application of the method of weighted residuals results in a weak form of this equation over the infinite element volume:

\int_{V_{f}} δ p (\frac{1}{K_{f}} \ddot{p} + \frac{γ}{ρ_{f} K_{f}} \dot{p} - \frac{\partial}{\partial x} \cdot (\frac{1}{ρ_{f}} \frac{\partial p}{\partial x})) d V = 0.

This equation is formally identical to that used in the finite element region (see Coupled acoustic-structural medium analysis); however, the choice of basis functions for the weight, $δ p$ , and for the solution field, $p$ , will be different. Selection of these functions has been the subject of considerable experimentation and analysis (for example, see Allik, 1991, and Burnett, 1994). In Abaqus considerations of accuracy, numerical stability, time-domain well-posedness, and economy have led to the selection of a form proposed by Astley (Astley, 1994). Here, the Fourier transform is applied to the equilibrium equation, and the functions for the solution field are chosen as tensor products of conventional finite element shape functions in the directions tangential to the terminating surface and polynomials in $1 / r$ , where r is a coordinate in the infinite direction, and a spatially oscillatory factor. The weights are chosen as complex conjugates of the solution field functions times a $1 / r^{2}$ factor. This combination has shown to result in an element with several desirable properties. First, if the material properties are constant with respect to frequency in steady-state analysis, the mass, damping, and stiffness element matrices are constant as well; equivalently, in transient analysis the element results in a second-order differential equation for the pressure. Moreover, the damping matrix has positive eigenvalues, a necessary condition for well-posedness in transient analysis. Analytical investigations of the formulation demonstrate that the element captures the exact solutions for radiation impedances for modes of a sphere; the order of the spherical mode modeled exactly is equal to the order of the polynomial used in the infinite direction. The element integrals do not contain oscillatory kernels and can be evaluated using standard Gauss quadrature methods. Finally, the element can be formulated using the usual coordinate map due to Bettess (1984) so that arbitrary convex terminating surfaces can be used.

To continue with the derivation, we transform the weighted residual statement into the frequency domain:

\begin{matrix} (1) & \int_{V_{f}} δ \tilde{p} (- Ω^{2} \frac{1}{K_{f}} \tilde{p} + i Ω \frac{γ}{ρ_{f} K_{f}} \tilde{p} - \frac{\partial}{\partial x} \cdot (\frac{1}{ρ_{f}} \frac{\partial \tilde{p}}{\partial x})) d V = 0. \end{matrix}

Now the weight and solution interpolation functions are defined as

\begin{matrix} δ \tilde{p} & ​ = δ p_{i} F ϕ_{i} (g, h, ξ) e^{i k (r - r_{0})}, \\ \tilde{p} & ​ = p_{j} ϕ_{j} (g, h, ξ) e^{- i k (r - r_{0})}, \end{matrix}

where $g (x)$ , $h (x)$ , and $ξ (x)$ are the parent element coordinates; $F \equiv {(r_{0} / r)}^{2}$ ; $k \equiv Ω / c_{0}$ ; $c_{0}^{2} \equiv K_{f} / ρ_{f}$ ; and $ϕ_{i}$ are element shape functions with indices varying over the number of degrees of freedom of the element. The mapping and shape functions are described below.

Inserting the shape functions into Equation 1 and integrating by parts, we obtain the following element equation:

\begin{matrix} (2) & \begin{matrix} \int_{V_{f}} (ϕ_{i} \frac{\partial ϕ_{j}}{\partial x} \cdot (\frac{1}{ρ_{f}} \frac{\partial F}{\partial x}) + F \frac{\partial ϕ_{i}}{\partial x} \cdot (\frac{1}{ρ_{f}} \frac{\partial ϕ_{j}}{\partial x})) d V p_{j} \\ - i k \int_{V_{f}} (\frac{\partial (r - r_{0})}{\partial x} \cdot (\frac{1}{ρ_{f}} \frac{\partial F}{\partial x}) ϕ_{i} ϕ_{j}) d V p_{j} \\ + i k \int_{V_{f}} \frac{F}{ρ_{f}} \frac{γ}{\sqrt{K_{f} ρ_{f}}} ϕ_{i} ϕ_{j} d V p_{j} \\ - i k \int_{V_{f}} \frac{F}{ρ_{f}} (\frac{\partial (r - r_{0})}{\partial x} \cdot \frac{\partial ϕ_{i}}{\partial x} ϕ_{j} - ϕ_{i} \frac{\partial ϕ_{j}}{\partial x} \cdot \frac{\partial (r - r_{0})}{\partial x}) d V p_{j} \\ - k^{2} \int_{V_{f}} \frac{F}{K_{f}} (1 - \frac{\partial (r - r_{0})}{\partial x} \cdot \frac{\partial (r - r_{0})}{\partial x}) ϕ_{i} ϕ_{j} d V p_{j} = 0. \end{matrix} \end{matrix}

This equation is clearly nonsymmetric, due to the fact that the weight and trial functions are not identical. Nevertheless, if the material properties are constant as a function of frequency, each element matrix corresponding to the terms above is constant as well. Gradients of density inside the element volume have been ignored in the formulation.

The element shape functions are defined as follows:

ϕ_{i} (g, h, r) \equiv f b_{α} B_{β},

where

f \equiv \sqrt{\frac{r_{0} (g, h)}{r}}

in two spatial dimensions and

f \equiv \frac{r_{0} (g, h)}{r}

in three. This function serves to specify the minimum rate of decay of the acoustic field inside the infinite element domain. The subindex $α$ ranges over the n nodes of the infinite element on the terminating surface, while the subindex $β$ ranges over the number of functions used in the infinite direction. The function index i is equal to the subindex $α$ for the first n functions, $i \equiv n + α$ for the second n functions, and so on.

The functions $b_{α} (g, h)$ are conventional two-dimensional shape functions (in three dimensions) or, with $h \equiv 0$ , one-dimensional shape functions for axisymmetric or two-dimensional elements. The role of these functions is to specify the variation of the acoustic field in directions tangent to the terminating surface. The variation of the acoustic field in the infinite direction is given by the functions $B_{β}$ , which are members of a set of ten ninth-order polynomials in $ξ \equiv 1 - \frac{2 r_{0} (g, h)}{r}$ . The members of this set are constructed to correspond to the Legendre modes of a sphere. That is, if infinite elements are placed on a sphere and if tangential refinement is adequate, an ith order acoustic infinite element will absorb waves associated with the ( $i - 1$ )th Legendre mode. The first member of this set corresponds to the value of the acoustic pressures on the terminating surface; the other functions are generalized degrees of freedom. Specifically,

B_{1} \equiv 1,

and

B_{β} \equiv (1 + ξ) {\hat{B}}_{β}, for β > 1,

where

\int_{- 1}^{1} {\hat{B}}_{β} {\hat{B}}_{γ} {(1 + ξ)}^{2} f^{2} F d ξ = 0, for β \neq γ .

This last condition is enforced to improve conditioning of the element matrices at higher order and results in different functions in two and three spatial dimensions. All of the $B_{β}$ are equal to zero at the terminating surface, except for $B_{1}$ .

The coordinate map is described in part by the element shape functions, in the usual isoparametric manner, and in part by a singular function. Together they map the true, semi-infinite domain onto the parent element square or cube. To specify the map for a given element, we first define distances $r_{α}$ between each infinite element node $α$ on the terminating surface and the element's reference node, located at $x_{0}$ :

r_{α} \equiv | x_{α} - x_{0} | .

Intermediate points in the infinite direction are defined as offset replicas of the nodes on the terminating surface:

y_{α} \equiv x_{α} + n_{α} r_{α},

where $n_{α}$ is the “node ray” at the node. The “node ray” at a node is the unit vector in the direction of the line between the reference point and the node.

We use the interpolated reference distance on the terminating surface,

r_{0} (g, h) \equiv \sum_{α} r_{α} b_{α} (g, h),

in the definition of the parent coordinate, $ξ$ , corresponding to the infinite direction:

ξ \equiv 1 - \frac{2 r_{0} (g, h)}{r},

where r is the distance between an arbitrary point in the infinite element volume and the reference point, $x_{0}$ . With these definitions, the geometric map can be specified as

x (g, h, ξ) \equiv \sum_{α} b_{α} (g, h) [x_{α} + (\frac{ξ + 1}{ξ - 1}) r_{α} n_{α}] .

The infinite elements are not isoparametric, since the map uses a lower-order function of the parent coordinates than the interpolation scheme does. However, this singular mapping is convenient and invertible.

The phase factor $(r - r_{0})$ appears in the element formulation. This factor models the oscillatory nature of the solution inside the infinite element volume but must also satisfy some additional properties. First, it must be zero at the face of the infinite element that connects to the finite element mesh. Second, it must be $C^{0}$ continuous across infinite element lateral boundaries. Third, it should be such that the mass-like term in the infinite element equation,

\int_{V_{f}} \frac{F}{K_{f}} (1 - \frac{\partial (r - r_{0})}{\partial x} \cdot \frac{\partial (r - r_{0})}{\partial x}) ϕ_{i} ϕ_{j} d V,

be positive semi-definite. This last criterion is not essential for execution of a steady-state analysis, but it is essential for the well-posedness of a transient analysis. Defining

q_{α} \equiv \min_{E} (n_{α} \cdot n_{E}),

where the subscript E refers to all elements at node $α$ , the definition

r - r_{0} (g, h) \equiv \sum_{α} r_{α} q_{α} b_{α} (g, h) (\frac{ξ}{ξ - 1})

satisfies these requirements. The inclusion of the factor $q_{α}$ is the only departure from the definition used in Astley (1994).

In a transient analysis the element matrix equation is transformed back to the time domain, using constant material properties taken from the value at zero frequency. The resulting second-order ordinary differential equation for the pressure degrees of freedom is added to the overall system in the model and integrated in the usual manner.

Steady-state analysis

In steady-state analysis the formulation for the acoustic infinite elements is consistent with that used for acoustic finite elements (see Coupled acoustic-structural medium analysis and Acoustic, Shock, and Coupled Acoustic-Structural Analysis). First, the complex density

\tilde{ρ} \equiv ρ_{f} + \frac{γ}{i Ω}

is defined. The resulting weighted residual statement is

- \int_{V_{f}} δ p \frac{Ω^{2}}{K_{f}} \tilde{p} d V + \int_{V_{f}} \frac{1}{\tilde{ρ}} \frac{\partial δ p}{\partial x} \cdot \frac{\partial \tilde{p}}{\partial x} d V + \int_{S} δ p \frac{1}{\tilde{ρ}} \frac{\partial \tilde{p}}{\partial x} \cdot n^{-} d S = 0

after integration by parts.

The expression for the basis and weight functions is formally identical to the preceding transient development, but with the complex wave number

\tilde{k} \equiv Ω \sqrt{\frac{\tilde{ρ}}{K_{f}}}

used in the definition of the spatially oscillatory term in the weight and basis functions:

e^{\pm \tilde{k} (r - r_{0})} .

After substitution and some simplification, the acoustic infinite element expression becomes

\begin{matrix} (3) & \begin{matrix} \int_{V_{f}} (ϕ_{i} \frac{\partial ϕ_{j}}{\partial x} \cdot (\frac{1}{\tilde{ρ}} \frac{\partial F}{\partial x}) + F \frac{\partial ϕ_{i}}{\partial x} \cdot (\frac{1}{\tilde{ρ}} \frac{\partial ϕ_{j}}{\partial x})) d V p_{j} \\ - i \int_{V_{f}} \tilde{k} (\frac{\partial (r - r_{0})}{\partial x} \cdot (\frac{1}{\tilde{ρ}} \frac{\partial F}{\partial x}) ϕ_{i} ϕ_{j}) d V p_{j} \\ - i \int_{V_{f}} \tilde{k} \frac{F}{\tilde{ρ}} (\frac{\partial (r - r_{0})}{\partial x} \cdot \frac{\partial ϕ_{i}}{\partial x} ϕ_{j} - ϕ_{i} \frac{\partial ϕ_{j}}{\partial x} \cdot \frac{\partial (r - r_{0})}{\partial x}) d V p_{j} \\ - \int_{V_{f}} \frac{{\tilde{k}}^{2}}{\tilde{ρ}} (1 - \frac{\partial (r - r_{0})}{\partial x} \cdot \frac{\partial (r - r_{0})}{\partial x}) F ϕ_{i} ϕ_{j} d V p_{j} = 0. \end{matrix} \end{matrix}

Now

\frac{{\tilde{k}}^{2}}{\tilde{ρ}} = \frac{Ω^{2}}{K_{f}},

so the mass matrix $M_{f}^{i j}$ is the same as for the transient case. It is real-valued.

The remaining terms are complex; and, as in the case for acoustic finite elements, they are manipulated into real and imaginary parts. We use

\frac{1}{\tilde{ρ}} = \frac{ρ_{f}}{ρ_{f}^{2} + γ^{2} / Ω^{2}} + i \frac{γ / Ω}{ρ_{f}^{2} + γ^{2} / Ω^{2}} \equiv R + i I

to obtain

\frac{\tilde{k}}{\tilde{ρ}} = [\frac{Ω}{{\sqrt{K}}_{f}}] {[\frac{1}{ρ_{f}^{2} + γ^{2} / Ω^{2}}]}^{(1 / 4)} [\cos θ / 2 + i \sin θ / 2] \equiv [\frac{Ω}{{\sqrt{K}}_{f}}] [R_{2} + i I_{2}],

where

\tan θ = (\frac{I}{R}) .

Using also

i \frac{\partial p}{\partial x} = \frac{1}{Ω} \frac{\partial \dot{p}}{\partial x},

we obtain real-valued element matrices as follows:

\begin{matrix} C_{f}^{i j} & ​ = \int_{V_{f}} \frac{I}{Ω} \frac{\partial [F ϕ_{i}]}{\partial x} \cdot \frac{\partial ϕ_{j}}{\partial x} d V, \\ + \int_{V_{f}} \frac{R_{2}}{{\sqrt{K}}_{f}} [- \frac{\partial F ϕ_{i}}{\partial x} \cdot \frac{\partial (r - r_{0})}{\partial x} ϕ_{j} + \frac{\partial ϕ_{j}}{\partial x} \cdot \frac{\partial (r - r_{0})}{\partial x} F ϕ_{i}] d V, \\ K_{f}^{i j} & ​ = \int_{V_{f}} R \frac{\partial [F ϕ_{i}]}{\partial x} \cdot \frac{\partial ϕ_{j}}{\partial x} d V, \\ - \int_{V_{f}} \frac{Ω I_{2}}{{\sqrt{K}}_{f}} [- \frac{\partial F ϕ_{i}}{\partial x} \cdot \frac{\partial (r - r_{0})}{\partial x} ϕ_{j} + \frac{\partial ϕ_{j}}{\partial x} \cdot \frac{\partial (r - r_{0})}{\partial x} F ϕ_{i}] d V . \end{matrix}

These are used directly in Equation 24.

Lateral face terms due to impedance conditions

An infinite element volume may have an impedance boundary condition defined on a face extending to infinity. Such a definition can be used to model, for example, an infinite half-space of water above a lossy seabed or under a layer of ice. The derivation of the contribution of the impedance condition to the infinite element matrix is straightforward and follows the derivation for finite elements (see Coupled acoustic-structural medium analysis). The contributing term in the transient case is

\int_{S_{lateral}} δ p (\frac{γ}{ρ_{f} c_{1}} p + (\frac{γ}{ρ_{f} k_{1}} + \frac{1}{c_{1}}) \dot{p} + \frac{1}{k_{1}} \ddot{p}) d S .

Applying the infinite element basis and weight functions, the lateral face contributes to the mass, damping, and stiffness matrices of the acoustic infinite element:

\begin{matrix} \int_{S_{lateral}} F ϕ_{i} ϕ_{j} (\frac{γ}{ρ_{f} c_{1}}) d S p_{j} + \int_{S_{lateral}} F ϕ_{i} ϕ_{j} (\frac{γ}{ρ_{f} k_{1}} + \frac{1}{c_{1}}) d S {\dot{p}}_{j} \\ + \int_{S_{lateral}} F ϕ_{i} ϕ_{j} (\frac{1}{k_{1}}) d S {\ddot{p}}_{j} . \end{matrix}

In steady-state and eigenvalue analysis the corresponding frequency-domain representation is used.

Coupling between acoustic and solid infinite elements

Lateral faces of acoustic and solid infinite elements (that is, those that extend to infinity) can be coupled to each other. For example, two infinite half-spaces may interact at their interface. As in the case of lateral-face impedance terms, the derivation of the contribution of the coupling condition to the infinite element matrix follows the derivation for finite elements to the greatest degree possible. We restrict the consideration to solid medium infinite elements that share an edge with the acoustic medium infinite elements and where the acoustic infinite element is the secondary in a tied contact condition. The contributing terms in the transient case are

- \int_{S_{lateral}} δ p (n^{-} \cdot {\ddot{u}}^{m}) d S

on the acoustic infinite element and

+ \int_{S_{lateral}} (n^{-} \cdot δ u^{m}) p d S

on the solid infinite element.

In tied contact between fluid and solid medium elements, a tributary area of the secondary node is defined, and the conservation of momentum and continuity of displacement is enforced between the secondary node and the main surface. On a lateral face of an infinite element, the “area” is unbounded, so an area measure based on the acoustic infinite element formulation,

A_{lateral} \equiv \int_{S_{lateral}} F f^{2} d S,

is defined. This area measure is used to compute a volumetric acceleration load on the acoustic infinite element:

- \int_{S_{lateral}} δ p (n^{-} \cdot {\ddot{u}}^{m}) d S \mapsto A_{lateral} (n^{-} \cdot {\ddot{u}}^{m}) .

The structural displacement, ${\ddot{u}}^{m}$ , is evaluated at the shared edge to compute this coupling condition. For consistency the boundary traction on the solid medium infinite element due to the acoustic pressure is defined as

+ \int_{S_{lateral}} (n^{-} \cdot δ u^{m}) p d S \mapsto A_{lateral} (n^{-} \cdot u^{m}) .

In steady-state and eigenvalue analysis the corresponding frequency-domain representation is used.