Coupled acoustic-structural medium analysis

Acoustic equations

The equilibrium equation for small motions of a compressible, adiabatic fluid with velocity-dependent momentum losses is taken to be

\begin{matrix} (1) & \frac{\partial p}{\partial x} + γ (x, θ_{i}) {\dot{u}}^{f} + ρ_{f} (x, θ_{i}) {\ddot{u}}^{f} = 0, \end{matrix}

where p is the excess pressure in the fluid (the pressure in excess of any static pressure); $x$ is the spatial position of the fluid particle; ${\dot{u}}^{f}$ is the fluid particle velocity; ${\ddot{u}}^{f}$ is the fluid particle acceleration; $ρ_{f}$ is the density of the fluid; $γ$ is the “volumetric drag” (force per unit volume per velocity); and $θ_{i}$ are i independent field variables such as temperature, humidity of air, or salinity of water on which $ρ_{f}$ and $γ$ may depend (see Acoustic medium). The d'Alembert term has been written without convection on the assumption that there is no steady flow of the fluid. This is usually considered sufficiently accurate for steady fluid velocities up to Mach 0.1.

The constitutive behavior of the fluid is assumed to be inviscid, linear, and compressible, so

\begin{matrix} (2) & p = - K_{f} (x, θ_{i}) \frac{\partial}{\partial x} \cdot u^{f}, \end{matrix}

where $K_{f}$ is the bulk modulus of the fluid.

For an acoustic medium capable of undergoing cavitation, the absolute pressure (sum of the static pressure and the excess dynamic pressure) cannot drop below the specified cavitation limit. When the absolute pressure drops to this limit value, the fluid is assumed to undergo free expansion without a corresponding drop in the dynamic pressure. The pressure would rebuild in the acoustic medium once the free expansion that took place during the cavitation is reversed sufficiently to reduce the volumetric strain to the level at the cavitation limit. The constitutive behavior for an acoustic medium capable of undergoing cavitation can be stated as

p = \max {p_{v}, p_{c} - p_{0}},

where a pseudopressure $p_{v}$ , a measure of the volumetric strain, is defined as

p_{v} = - K_{f} (x, θ_{i}) \frac{\partial}{\partial x} \cdot u^{f},

where $p_{c}$ is the fluid cavitation limit and $p_{0}$ is the initial acoustic static pressure. A total wave formulation is used for a nonlinear acoustic medium undergoing cavitation. This formulation is very similar to the scattered wave formulation presented below except that the pseudopressure, defined as the product of the bulk modulus and the compressive volumetric strain, plays the role of the material state variable instead of the acoustic excess pressure. The acoustic excess pressure is readily available from this pseudopressure subject to the cavitation condition.

Physical boundary conditions in acoustic analysis

Acoustic fields are strongly dependent on the conditions at the boundary of the acoustic medium. The boundary of a region of acoustic medium that obeys Equation 1 and Equation 2 can be divided into subregions S on which the following conditions are imposed:

$S_{fp}$ ,: where the value of the acoustic pressure p is prescribed.
$S_{ft}$ ,: where we prescribe the normal derivative of the acoustic medium. This condition also prescribes the motion of the fluid particles and can be used to model acoustic sources, rigid walls (baffles), incident wave fields, and symmetry planes.
$S_{fr}$ ,: the “reactive” acoustic boundary, where there is a prescribed linear relationship between the fluid acoustic pressure and its normal derivative. Quite a few physical effects can be modeled in this manner: in particular, the effect of thin layers of material, whose own motions are unimportant, placed between acoustic media and rigid baffles. An example is the carpet glued to the floor of a room or car interior that absorbs and reflects acoustic waves. This thin layer of material provides a “reactive surface,” or impedance boundary condition, to the acoustic medium. This type of boundary condition is also referred to as an imposed impedance, admittance, or a “Dirichlet to Neumann map.”
$S_{fi}$ ,: the “radiating” acoustic boundary. Often, acoustic media extend sufficiently far from the region of interest that they can be modeled as infinite in extent. In such cases it is convenient to truncate the computational region and apply a boundary condition to simulate waves passing exclusively outward from the computational region.
$S_{fs}$ ,: where the motion of an acoustic medium is directly coupled to the motion of a solid. On such an acoustic-structural boundary the acoustic and structural media have the same displacement normal to the boundary, but the tangential motions are uncoupled.
$S_{frs}$ ,: an acoustic-structural boundary, where the displacements are linearly coupled but not necessarily identically equal due to the presence of a compliant or reactive intervening layer. This layer induces an impedance condition between the relative normal velocity between acoustic fluid and solid structure and the acoustic pressure. It is analogous to a spring and dashpot interposed between the fluid and solid particles. As implemented in Abaqus, an impedance boundary condition surface does not model any mass associated with the reactive lining; if such a mass exists, it should be incorporated into the boundary of the structure.
$S_{ff}$ ,: a boundary between acoustic fluids of possibly differing material properties. On such an interface, displacement continuity requires that the normal forces per unit mass on the fluid particles be equal. This quantity is the natural boundary traction in Abaqus, so this condition is enforced automatically during element assembly. This is also true in one-dimensional analysis (i.e., piping or ducts), where the relevant acoustic properties include the cross-sectional areas of the elements. Consequently, fluid-fluid boundaries do not require special treatment in Abaqus.

Formulation for direct integration transient dynamics

In Abaqus the finite element formulations are slightly different in direct integration transient and steady-state or modal analyses, primarily with regard to the treatment of the volumetric drag loss parameter and spatial variations of the constitutive parameters. To derive a symmetric system of ordinary differential equations for implicit integration, some approximations are made in the transient case that are not needed in steady state. For linear transient dynamic analysis, the modal procedure can be used and is much more efficient.

To derive the partial differential equation used in direct integration transient analysis, we divide Equation 1 by $ρ_{f}$ , take its gradient with respect to $x$ , neglect the gradient of $γ / ρ_{f}$ , and combine the result with the time derivatives of Equation 2 to obtain the equation of motion for the fluid in terms of the fluid pressure:

\begin{matrix} (3) & \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. \end{matrix}

The assumption that the gradient of $γ / ρ_{f}$ is small is violated where there are discontinuities in the quantity $γ / ρ_{f}$ (for example, on the boundary between two elements that have a different $γ / ρ_{f}$ value).

Variational statement

An equivalent weak form for the equation of motion, Equation 3, is obtained by introducing an arbitrary variational field, $δ p$ , and integrating over the fluid:

\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.

Green's theorem allows this to be rewritten as

\begin{matrix} (4) & \int_{V_{f}} [δ p (\frac{1}{K_{f}} \ddot{p} + \frac{γ}{ρ_{f} K_{f}} \dot{p}) + \frac{1}{ρ_{f}} \frac{\partial δ p}{\partial x} \cdot \frac{\partial p}{\partial x}] d V + \int_{S} δ p (\frac{1}{ρ_{f}} n^{-} \cdot \frac{\partial p}{\partial x}) d S = 0. \end{matrix}

Assuming that p is prescribed on $S_{fp}$ , the equilibrium equation, Equation 1, is used on the remainder of the boundary to relate the pressure gradient to the motion of the boundary:

\begin{matrix} (5) & n^{-} \cdot (\frac{1}{ρ_{f}} \frac{\partial p}{\partial x} + \frac{γ}{ρ_{f}} {\dot{u}}^{f} + {\ddot{u}}^{f}) = 0 on S - S_{fp} . \end{matrix}

Using this equation, the term $n^{-} \cdot \frac{\partial p}{\partial x}$ is eliminated from Equation 4 to produce

\begin{matrix} (6) & \int_{V_{f}} [δ p (\frac{1}{K_{f}} \ddot{p} + \frac{γ}{ρ_{f} K_{f}} \dot{p}) + \frac{1}{ρ_{f}} \frac{\partial δ p}{\partial x} \cdot \frac{\partial p}{\partial x}] d V - \int_{S - S_{fp}} δ p (T (x)) d S = 0, \end{matrix}

where, for convenience, the boundary “traction” term

\begin{matrix} (7) & T (x) = n^{-} \cdot ({\ddot{u}}^{f} + \frac{γ}{ρ_{f}} {\dot{u}}^{f}) = - n^{-} \cdot (\frac{1}{ρ_{f}} \frac{\partial p}{\partial x}) on S - S_{fp} \end{matrix}

has been introduced.

Except for the imposed pressure on $S_{fp}$ , all the other boundary conditions described above can be formulated in terms of $T (x)$ . This term has dimensions of acceleration; in the absence of volumetric drag this boundary traction is equal to the inward acceleration of the particles of the acoustic medium:

\begin{matrix} (8) & T (x) = n^{-} \cdot {\ddot{u}}^{f} on S - S_{fp} . \end{matrix}

When volumetric drag is present, the boundary traction is the normal derivative of the pressure field, divided by the true mass density: a force per unit mass of fluid. Consequently, when volumetric drag exists in a transient acoustic model, a unit of $T (x)$ yields a lower local volumetric acceleration, due to drag losses.

In direct integration transient dynamics we enforce the acoustic boundary conditions as follows:

On $S_{fp}$ ,

p is prescribed and $δ p = 0$ .

On $S_{ft}$ ,

where we prescribe the normal derivative of the acoustic pressure per unit density:

T_{ft} (x) \equiv T_{0} .

In the absence of volumetric drag in the medium, this enforces a value of fluid particle acceleration, $n^{-} \cdot {\ddot{u}}^{f} = T_{0} = a_{i n}$ . An imposed $T_{0} = a_{i n}$ can be used to model the oscillations of a rigid plate or body exciting a fluid, for example. A special case of this boundary condition is $a_{i n} = 0$ , which represents a rigid immobile boundary. As mentioned above, if the medium has nonzero volumetric drag, a unit of $T_{0}$ imposed at the boundary will result in a relatively lower imposed particle acceleration. Incident wave fields on a boundary of a fluid are modeled with a $T_{0}$ that varies in space and time, corresponding to the effect of the arrival of the wave on the boundary. See Loading due to an incident dilatational wave field.

On $S_{fr}$ ,

the reactive boundary between the acoustic medium and a rigid baffle, we apply a condition that relates the velocity of the acoustic medium to the pressure and rate of change of pressure:

\begin{matrix} (9) & - n^{-} \cdot {\dot{u}}^{f} = (\frac{1}{k_{1}} \dot{p} + \frac{1}{c_{1}} p) on S_{fr}, \end{matrix}

where $1 / k_{1}$ and $1 / c_{1}$ are user-prescribed parameters at the boundary. This equation is in the form of an admittance relation; the impedance expression is simply the inverse. The layer of material, in admittance form, acts as a spring and dashpot in series distributed between the acoustic medium and the rigid wall. The spring and dashpot parameters are $k_{1}$ and $c_{1}$ , respectively; they are per unit area of the acoustic boundary. Using this definition for the fluid velocity, the boundary tractions in the variational statement become

\begin{matrix} (10) & T_{fr} (x) \equiv - (\frac{γ}{ρ_{f}} \frac{1}{c_{1}} p + (\frac{γ}{ρ_{f}} \frac{1}{k_{1}} + \frac{1}{c_{1}}) \dot{p} + \frac{1}{k_{1}} \ddot{p}) . \end{matrix}

On $S_{fi}$ ,

the radiating boundary, we apply the radiation boundary condition by specifying the corresponding impedance:

T_{fi} (x) \equiv - (\frac{1}{c_{1}} \dot{p} + \frac{1}{a_{1}} p),

using the admittance parameters of Equation 47 and Equation 48, defined below.

On $S_{fs}$ ,

the acoustic-structural interface, we apply the acoustic-structural interface condition by equating displacement of the fluid and solid, which enforces the condition

n^{-} \cdot u^{f} = n^{-} \cdot u^{m},

where $u^{m}$ is the displacement of the structure. In the presence of volumetric drag it follows that the acoustic boundary traction coupling fluid to solid is

T (x) = n^{-} \cdot ({\ddot{u}}^{m} + \frac{γ}{ρ_{f}} {\dot{u}}^{m}) .

In Abaqus/Standard the formulation of the transient coupled problem would be made nonsymmetric by the presence of the term $n^{-} \cdot (\frac{γ}{ρ_{f}} {\dot{u}}^{m})$ . In the great majority of practical applications the acoustic tractions associated with volumetric drag are small compared to those associated with fluid inertia,

{\ddot{u}}^{m} ≫ \frac{γ}{ρ_{f}} {\dot{u}}^{m} \forall u^{m} (t),

so this term is ignored in transient analysis:

T_{fs} (x) \equiv n^{-} \cdot {\ddot{u}}^{m} .

On $S_{frs}$ ,

the mixed impedance boundary and acoustic-structural boundary, we apply a condition that relates the relative outward velocity between the acoustic medium and the structure to the pressure and rate of change of pressure:

\begin{matrix} (11) & n^{-} \cdot ({\dot{u}}^{m} - {\dot{u}}^{f}) = \frac{1}{k_{1}} \dot{p} + \frac{1}{c_{1}} p on S_{frs} . \end{matrix}

This relative normal velocity represents a rate of compression (or extension) of the intervening layer. Applying this equation to the definition of $T (x)$ , we obtain for the transient case

\begin{matrix} (12) & T_{frs} (x) = n^{-} \cdot ({\ddot{u}}^{m}) - \frac{γ}{ρ_{f}} \frac{1}{c_{1}} p - (\frac{γ}{ρ_{f}} \frac{1}{k_{1}} + \frac{1}{c_{1}}) \dot{p} - \frac{1}{k_{1}} \ddot{p} . \end{matrix}

This expression for $T (x)$ is the sum of its definitions for $S_{fs}$ and $S_{fr}$ . In the steady-state case the effect of volumetric drag on the structural displacement term in the acoustic traction is included:

\begin{matrix} (13) & T_{frs} (x) = n^{-} \cdot ({\ddot{u}}^{m} + \frac{γ}{ρ_{f}} {\dot{u}}^{m}) - \frac{γ}{ρ_{f}} \frac{1}{c_{1}} p - (\frac{γ}{ρ_{f}} \frac{1}{k_{1}} + \frac{1}{c_{1}}) \dot{p} - \frac{1}{k_{1}} \ddot{p} . \end{matrix}

These definitions for the boundary term, $T (x)$ , are introduced into Equation 6 to give the final variational statement for the acoustic medium (this is the equivalent of the virtual work statement for the structure):

\int_{V_{f}} [δ p (\frac{1}{K_{f}} \ddot{p} + \frac{γ}{ρ_{f} K_{f}} \dot{p}) + \frac{1}{ρ_{f}} \frac{\partial δ p}{\partial x} \cdot \frac{\partial p}{\partial x}] d V - \int_{S_{ft}} δ p T_{0} d S

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

+ \int_{S_{fi}} δ p (\frac{1}{c_{1}} \dot{p} + \frac{1}{a_{1}} p) d S - \int_{S_{fs}} δ p n^{-} \cdot {\ddot{u}}^{m} d S

\begin{matrix} (14) & + \int_{S_{frs}} δ p (\frac{γ}{ρ_{f} c_{1}} p + (\frac{γ}{ρ_{f} k_{1}} + \frac{1}{c_{1}}) \dot{p} + \frac{1}{k_{1}} \ddot{p} - n^{-} \cdot {\ddot{u}}^{m}) d S = 0. \end{matrix}

The structural behavior is defined by the virtual work equation,

\int_{V} δ ε : σ d V + \int_{V} α_{c} ρ δ u^{m} \cdot {\dot{u}}^{m} d V + \int_{V} ρ δ u^{m} \cdot {\ddot{u}}^{m} d V

\begin{matrix} (15) & + \int_{S_{fs}} p δ u^{m} \cdot n d S - \int_{S_{t}} δ u^{m} \cdot t d S = 0, \end{matrix}

where $σ$ is the stress at a point in the structure, p is the pressure acting on the fluid-structural interface, $n$ is the outward normal to the structure, $ρ$ is the density of the material, $α_{c}$ is the mass proportional damping factor (part of the Rayleigh damping assumption for the structure), ${\ddot{u}}^{m}$ is the acceleration of a point in the structure, $t$ is the surface traction applied to the structure, $δ u^{m}$ is a variational displacement field, and $δ ε$ is the strain variation that is compatible with $δ u^{m}$ . For simplicity in this equation all other loading terms except the fluid pressure and surface traction $t$ have been neglected: they are imposed in the usual way.

The discretized finite element equations

Equation 14 and Equation 15 define the variational problem for the coupled fields $u^{m}$ and p. The problem is discretized by introducing interpolation functions: in the fluid $p = H^{P} p^{P}$ , $P = 1, 2 \dots$ up to the number of pressure nodes and in the structure $u^{m} = N^{N} u^{N}$ , $N = 1, 2 \dots$ up to the number of displacement degrees of freedom. In these and the following equations we assume summation over the superscripts that refer to the degrees of freedom of the discretized model. We also use the superscripts $P$ , $Q$ to refer to pressure degrees of freedom in the fluid and $N$ , $M$ to refer to displacement degrees of freedom in the structure. We use a Galerkin method for the structural system; the variational field has the same form as the displacement: $δ u^{m} = N^{N} δ u^{N}$ . For the fluid we use $δ p = H^{P} δ p^{P}$ but with the subsequent Petrov-Galerkin substitution

\begin{matrix} (16) & δ p^{P} = \frac{d^{2}}{d t^{2}} (δ {\hat{p}}^{P}) . \end{matrix}

The new function $δ {\hat{p}}^{P}$ makes the single variational equation obtained from summing Equation 14 and Equation 15 dimensionally consistent:

\begin{matrix} (17) & \begin{matrix} - δ {\hat{p}}^{P} & {(M_{f}^{P Q} + M_{fr}^{P Q}) {\ddot{p}}^{Q} + (C_{f}^{P Q} + C_{fr}^{P Q}) {\dot{p}}^{Q} + (K_{f}^{P Q} + K_{fr}^{P Q} + K_{fi}^{P Q}) p^{Q} - S_{fs}^{P M} {\ddot{u}}^{M} - P_{f}^{P}} \\ + δ u^{N} & {I^{N} + M^{N M} {\ddot{u}}^{M} + C_{(m)}^{N M} {\dot{u}}^{M} + {[S_{fs}^{Q N}]}^{T} p^{Q} - P^{N}} = 0, \end{matrix} \end{matrix}

where, for simplicity, we have introduced the following definitions:

\begin{matrix} M_{f}^{P Q} & = \int_{V_{f}} \frac{1}{K_{f}} H^{P} H^{Q} d V, \\ M_{fr}^{P Q} & = \int_{S_{fr} \cup S_{frs}} \frac{1}{k_{1}} H^{P} H^{Q} d S, \\ C_{f}^{P Q} & = \int_{V_{f}} \frac{γ}{ρ_{f}} \frac{1}{K_{f}} H^{P} H^{Q} d V, \\ C_{fr}^{P Q} & = \int_{S_{fr} \cup S_{frs}} (\frac{γ}{ρ_{f}} \frac{1}{k_{1}} + \frac{1}{c_{1}}) H^{P} H^{Q} d S + \int_{S_{fi}} \frac{1}{c_{1}} H^{P} H^{Q} d S, \\ K_{f}^{P Q} & = \int_{V_{f}} \frac{1}{ρ_{f}} \frac{\partial H^{P}}{\partial x} \cdot \frac{\partial H^{Q}}{\partial x} d V, \\ K_{fr}^{P Q} & = \int_{S_{fr} \cup S_{frs}} \frac{γ}{ρ_{f}} \frac{1}{c_{1}} H^{P} H^{Q} d S, \end{matrix}

\begin{matrix} K_{fi}^{P Q} & = \int_{S_{fi}} \frac{1}{a_{1}} H^{P} H^{Q} d S, \\ S_{fs}^{P M} & = \int_{S_{fs} \cup S_{frs}} H^{P} n^{-} \cdot N^{M} d S, \\ P_{f}^{P} & = \int_{S_{ft}} H^{P} T_{0} d S, \\ M^{N M} & = \int_{V} ρ N^{N} \cdot N^{M} d V, \\ C_{(m)}^{N M} & = \int_{V} α_{c} ρ N^{N} \cdot N^{M} d V, \\ I^{N} & = \int_{V} β^{N} : σ d V, \\ P^{N} & = \int_{S_{t}} N^{N} \cdot t d S, \end{matrix}

where $β^{N}$ is the strain interpolator. This equation defines the discretized model. We see that the volumetric drag-related terms are “mass-like”; i.e., proportional to the fluid element mass matrix.

The term $P_{f}^{P}$ is the nodal right-hand-side term for the acoustical degree of freedom $p^{P}$ , or the applied “force” on this degree of freedom. This term is obtained by integration of the normal derivative of pressure per unit density of the acoustic medium over the surface area tributary to a boundary node.

In the case of coupled systems where the forces on the structure due to the fluid— ${[S_{fs}^{Q N}]}^{T} p^{Q}$ are very small compared to the rest of the structural forces—the system can be solved in a “sequentially coupled” manner. The structural equations can be solved with the ${[S_{fs}^{Q N}]}^{T} p^{Q}$ term omitted; i.e., in an analysis without fluid coupling. Subsequently, the fluid equations can be solved, with $[S_{fs}^{P M}] {\ddot{u}}^{M}$ imposed as a boundary condition. This two-step analysis is less expensive and advantageous for systems such as metal structures in air.

Time integration

The equations are integrated through time using the standard implicit (Abaqus/Standard) and explicit (Abaqus/Explicit) dynamic integration options. From the implicit integration operator we obtain relations between the variations of the solution variables (here represented by $f$ ) and their time derivatives:

D_{a} \overset{def}{=} \frac{δ \ddot{f}}{δ f} = \frac{δ p^{P}}{δ {\hat{p}}^{P}}; D_{v} \overset{def}{=} \frac{δ \dot{f}}{δ f} .

The equations of evolution of the degrees of freedom can be written for the implicit case as

\begin{matrix} - δ p^{P} & \frac{1}{D_{a}} {(M_{f}^{P Q} + M_{fr}^{P Q}) {\ddot{p}}^{Q} + (C_{f}^{P Q} + C_{fr}^{P Q}) {\dot{p}}^{Q} \\ + (K_{f}^{P Q} + K_{fr}^{P Q} + K_{fi}^{P Q}) p^{Q} - S_{fs}^{P M} {\ddot{u}}^{M} - P_{f}^{P}} \\ + δ u^{N} {I^{N} + M^{N M} {\ddot{u}}^{M} + C_{(m)}^{N M} {\dot{u}}^{M} + {[S_{fs}^{Q N}]}^{T} p^{Q} - P^{N}} = 0. \end{matrix}

The linearization of this equation is

\begin{matrix} - δ p^{P} & {(M_{f}^{P Q} + M_{fr}^{P Q}) + \frac{D_{v}}{D_{a}} (C_{f}^{P Q} + C_{fr}^{P Q}) + \frac{1}{D_{a}} (K_{f}^{P Q} + K_{fr}^{P Q} + K_{fi}^{P Q})} d p^{Q} \\ + δ p^{P} S_{fs}^{P M} d u^{M} + δ u^{N} {[S_{fs}^{Q N}]}^{T} d p^{Q} \\ + δ u^{N} {K^{N M} + D_{a} M^{N M} + D_{v} (C_{(m)}^{N M} + C_{(k)}^{N M})} d u^{M} = 0, \end{matrix}

where $d p$ and $d u$ are the corrections to the solution obtained from the Newton iteration, $K^{N M}$ is the structural stiffness matrix, and $C_{(k)}^{N M}$ is the structural damping matrix. These equations are symmetric if the constituent stiffness, damping, and mass matrices are symmetric.

For explicit integration the fluid mass matrix is diagonalized in a manner similar to the treatment of structural mass. The explicit central difference procedure described in Explicit dynamic analysis is applied to the coupled system of equations.

Summary of additional approximations of the direct integration transient formulation

As mentioned above, derivation of symmetric ordinary differential equations in the presence of volumetric drag requires some approximations in addition to those inherent in any finite element method. First, the spatial gradients of the ratio of volumetric drag to mass density in the fluid are neglected. This may be important in lossy, inhomogeneous acoustic media. Second, to maintain symmetry, the effect of volumetric drag on the fluid-solid boundary terms is neglected. Finally, the effect of volumetric drag on the radiation boundary conditions is approximate. If any of these effects is expected to be significant in an analysis, the user should realize that the results obtained are approximate.

Formulation for steady-state response using nodal degrees of freedom

The direct-solution steady-state dynamic analysis procedure is to be preferred over the transient formulation if volumetric drag is significant. This formulation uses the nodal degrees of freedom in the solid and acoustic regions directly to form a large linear system of equations defining the coupled structural-acoustic mechanics at a single frequency. If volumetric drag effects are not significant, the mode-based procedures (see below) are preferred because of their efficiency.

All model degrees of freedom and loads are assumed to be varying harmonically at an angular frequency $Ω$ , so we can write

f = \tilde{f} \exp i Ω t,

where $\tilde{f}$ is the constant complex amplitude of the variable $f$ . Thus,

\dot{f} = i Ω f, \ddot{f} = - Ω^{2} f .

We begin with the equilibrium equation

\frac{\partial p}{\partial x} + γ {\dot{u}}^{f} + ρ_{f} {\ddot{u}}^{f} = 0

and use the harmonic time-derivative relations to obtain

\frac{\partial \tilde{p}}{\partial x} - Ω^{2} (ρ_{f} + \frac{γ}{i Ω}) {\tilde{u}}^{f} = 0.

We define the complex density, $\tilde{ρ}$ , as

\begin{matrix} (18) & \tilde{ρ} \equiv ρ_{f} + \frac{γ}{i Ω} \end{matrix}

and, thus, write

\begin{matrix} (19) & \frac{\partial \tilde{p}}{\partial x} - Ω^{2} (\tilde{ρ}) {\tilde{u}}^{f} = 0. \end{matrix}

The equilibrium equation is now in a form where the density is complex and the acoustic medium velocity does not enter. We divide this equation by $\tilde{ρ}$ and combine it with the second time derivative of the constitutive law, Equation 2, to obtain

\begin{matrix} (20) & - Ω^{2} \frac{1}{K_{f}} \tilde{p} - \frac{\partial}{\partial x} \cdot (\frac{1}{\tilde{ρ}} \frac{\partial \tilde{p}}{\partial x}) = 0. \end{matrix}

We have not used the assumption that the spatial gradient of $r / ρ_{f}$ is small, as was done in the transient dynamics formulation.

Variational statement

The development of the variational statement parallels that for the case of transient dynamics, as though the volumetric drag were absent and the density complex. The variational statement is

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

Integrating by parts, we have

- \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.

In steady state the boundary traction is defined as

\tilde{T} (x) \equiv - \frac{1}{\tilde{ρ}} \frac{\partial \tilde{p}}{\partial x} \cdot n^{-} = - Ω^{2} n^{-} \cdot {\tilde{u}}^{f} = n^{-} \cdot {\ddot{u}}^{f} .

This expression is not the Fourier transform of the boundary traction defined above for the transient case. The steady-state definition is based on the complex density and includes the volumetric drag effect in such a way that it is always equal to the acceleration of the fluid particles. The application of boundary conditions may be slightly different for some cases in steady state due to this definition of the traction.

On $S_{fp}$ ,

$\tilde{p}$ is prescribed, analogous to transient analysis.

On $S_{ft}$ ,

we prescribe

{\tilde{T}}_{ft} (x) \equiv T_{0} .

The condition $n^{-} \cdot {\ddot{u}}^{f} = T_{0} = a_{i n}$ is enforced, even in the presence of volumetric drag.

On $S_{fr}$ ,

the reactive boundary between the acoustic medium and a rigid baffle, we apply

\begin{matrix} (21) & {\tilde{T}}_{fr} (x) \equiv - (\frac{i Ω}{c_{1}} - \frac{Ω^{2}}{k_{1}}) \tilde{p} . \end{matrix}

On $S_{fi}$ ,

the radiating boundary, we apply the radiation boundary condition impedance in the same form as for the reactive boundary but with the parameters as defined in Equation 42 and Equation 43.

On $S_{fs}$ ,

the acoustic-structural interface, we equate the displacement of the fluid and solid as in the transient case. However, the acoustic boundary traction coupling fluid to solid,

\tilde{T} (x) = - Ω^{2} n^{-} \cdot {\tilde{u}}^{m},

can be applied without affecting the symmetry of the overall formulation. Consequently, the acoustic tractions in the steady-state case make no assumptions about volumetric drag.

On $S_{frs}$ ,

the mixed impedance boundary and acoustic-structural boundary, the condition

n^{-} \cdot ({\dot{u}}^{m} - {\dot{u}}^{f}) = \frac{1}{k_{1}} \dot{p} + \frac{1}{c_{1}} p on S_{frs}

results in the definition

{\tilde{T}}_{frs} (x) = - Ω^{2} n^{-} \cdot ({\tilde{u}}^{m}) - \frac{i Ω}{c_{1}} \tilde{p} + \frac{Ω^{2}}{k_{1}} \tilde{p} .

In this case the effect of volumetric drag is included without approximation.

The final variational statement becomes

\begin{matrix} \int_{V_{f}} [- Ω^{2} δ p (\frac{1}{K_{f}} \tilde{p}) + \frac{1}{\tilde{ρ}} \frac{\partial δ p}{\partial x} \cdot \frac{\partial \tilde{p}}{\partial x}] d V \\ - \int_{S_{ft}} δ p a_{i n} d S + \int_{S_{fr} \cup S_{fi}} δ p (\frac{i Ω}{c_{1}} - \frac{Ω^{2}}{k_{1}}) \tilde{p} d S \\ + \int_{S_{fs}} δ p Ω^{2} n^{-} \cdot {\tilde{u}}^{m} d S \\ + \int_{S_{frs}} δ p (\frac{i Ω}{c_{1}} \tilde{p} - \frac{Ω^{2}}{k_{1}} \tilde{p} + Ω^{2} n^{-} \cdot {\tilde{u}}^{m}) d S = 0. \end{matrix}

This equation is formally identical to Equation 4, except for the pressure “stiffness” term, the radiation boundary conditions, and the imposed boundary traction term. Because the volumetric drag effect is contained in the complex density, the acoustic-structural boundary term in this formulation does not have the limitation that the volumetric drag must be small compared to other effects in the acoustic medium. In addition, in this formulation the applied flux on an acoustic boundary represents the inward acceleration of the acoustic medium, whether or not the volumetric drag is large. Finally, the radiation boundary conditions do not make any approximations with regard to the volumetric drag parameter.

The above equation uses the complex density, $1 / \tilde{ρ}$ . We manipulate it into a form that has real coefficients and an additional time derivative through the relations

\frac{1}{\tilde{ρ}} = \frac{ρ_{f}}{ρ_{f}^{2} + γ^{2} / Ω^{2}} + i \frac{γ / Ω}{ρ_{f}^{2} + γ^{2} / Ω^{2}}, i \frac{\partial p}{\partial x} = \frac{1}{Ω} \frac{\partial \dot{p}}{\partial x},

to obtain

\begin{matrix} (22) & \begin{matrix} \int_{V_{f}} [- δ p (\frac{Ω^{2}}{K_{f}} \tilde{p}) & + \frac{ρ_{f}}{ρ_{f}^{2} + γ^{2} / Ω^{2}} \frac{\partial δ p}{\partial x} \cdot \frac{\partial \tilde{p}}{\partial x} + (i Ω) \frac{γ / Ω^{2}}{ρ_{f}^{2} + γ^{2} / Ω^{2}} \frac{\partial δ p}{\partial x} \cdot \frac{\partial \tilde{p}}{\partial x}] d V \\ - \int_{S_{ft}} δ p a_{i n} d S + \int_{S_{fs} \cup S_{frs}} δ p Ω^{2} n^{-} \cdot {\tilde{u}}^{m} d S \\ + \int_{S_{fr} \cup S_{fi} \cup S_{frs}} δ p (\frac{i Ω}{c_{1}} \tilde{p} - \frac{Ω^{2}}{k_{1}} \tilde{p}) d S = 0. \end{matrix} \end{matrix}

The discretized finite element equations

Applying Galerkin's principle, the finite element equations are derived as before. We arrive again at Equation 17 with the same matrices except for the damping and stiffness matrices of the acoustic elements and the surfaces that have imposed impedance conditions, which now appear as

\begin{matrix} C_{f}^{P Q} & = \int_{V_{f}} \frac{γ / Ω^{2}}{ρ_{f}^{2} + γ^{2} / Ω^{2}} \frac{\partial H^{P}}{\partial x} \cdot \frac{\partial H^{Q}}{\partial x} d V, \\ K_{f}^{P Q} & = \int_{V_{f}} \frac{ρ_{f}}{ρ_{f}^{2} + γ^{2} / Ω^{2}} \frac{\partial H^{P}}{\partial x} \cdot \frac{\partial H^{Q}}{\partial x} d V, \\ C_{fr}^{P Q} & = \int_{S_{fr}} \frac{1}{c_{1}} H^{P} H^{Q} d S, \\ K_{fr}^{P Q} & = 0, \\ K_{fi}^{P Q} & = 0. \end{matrix}

The matrix modeling loss to volumetric drag is proportional to the fluid stiffness matrix in this formulation.

For steady-state harmonic response we assume that the structure undergoes small harmonic vibrations, identified by the prefix $Δ$ , about a deformed, stressed base state, which is identified by the subscript $0$ . Hence, the total stress can be written in the form

σ = σ_{0} + Δ σ = σ_{0} + D^{e l} : (Δ ε + i s Δ ε + β_{c} Δ \dot{ε}),

where $σ_{0}$ is the stress in the base state; $D^{e l}$ is the elasticity matrix for the material; $β_{c}$ is the stiffness proportional damping factor chosen for the material (to give the stiffness proportional contribution to the Rayleigh damping, thus introducing the viscous part of the material behavior); and, from the discretization assumption,

Δ ε = β^{M} Δ u^{M} .

To solve the steady-state problem, we assume that the governing equations are satisfied in the base state, and we linearize these equations in terms of the harmonic oscillations. For the internal force vector this yields

Δ I^{N} = (i C_{(s)}^{N M} + K^{N M}) Δ u^{M} + C_{(k)}^{N M} Δ {\dot{u}}^{M};

and Equation 17 can be rewritten, using the time-harmonic relations, as

\begin{matrix} (23) & \begin{matrix} - δ {\hat{p}}^{P} & {[- Ω^{2} (M_{f}^{P Q} + M_{fr}^{P Q}) + i Ω (C_{f}^{P Q} + C_{fr}^{P Q}) + i C_{f (s)}^{P Q} + K_{f}^{P Q}] Δ {\tilde{p}}^{Q} + Ω^{2} S_{fs}^{P M} Δ {\tilde{u}}^{M} - Δ {\tilde{P}}_{f}^{P}} \\ + δ u^{N} & {[- Ω^{2} M^{N M} + i Ω (C_{(m)}^{N M} + C_{(k)}^{N M}) + i C_{(s)}^{N M} + K^{N M}] Δ {\tilde{u}}^{M} + {[S_{fs}^{Q N}]}^{T} Δ {\tilde{p}}^{Q} - Δ {\tilde{P}}^{N}} = 0, \end{matrix} \end{matrix}

with

K^{N M} = \int_{V} [\frac{\partial β^{N}}{\partial u^{M}} : σ_{0} + β^{N} : D^{e l} : β^{M}] d V

(this stiffness includes the initial stress matrix, so “stress stiffening” and “load stiffness” effects associated with the base state stress and loads are included) and

C_{(k)}^{N M} = \int_{V} [β_{c} β^{N} : D^{e l} : β^{M}] d V .

We have also added the “global” parts of the “structural damping” terms

C_{(s)}^{N M} = \int_{V} s β^{N} : D^{e l} : β^{M} d V + s_{g} K^{N M}

and

C_{f (s)}^{P Q} \equiv s_{g} K_{f}^{P Q}

to the equation. These damping terms model finite energy loss in the low-frequency limit in steady-state dynamics—see Direct steady-state dynamic analysis and Subspace-based steady-state dynamic analysis. It should be noted that the acoustic “structural damping” operator inherits the frequency dependence of the acoustic stiffness matrix caused by volumetric drag.

We assume that the loads and (because of linearity) the response are harmonic; hence, we can write

\begin{matrix} Δ {\tilde{p}}^{Q} & = (ℜ ({\tilde{p}}^{Q}) + i ℑ ({\tilde{p}}^{Q})) \exp i Ω t \\ Δ {\tilde{u}}^{M} & = (ℜ ({\tilde{u}}^{M}) + i ℑ ({\tilde{u}}^{M})) \exp i Ω t \end{matrix}

and

\begin{matrix} Δ {\tilde{P}}^{N} & = (ℜ ({\tilde{P}}^{N}) + i ℑ ({\tilde{P}}^{N})) \exp i Ω t \\ Δ {\tilde{P}}_{f}^{P} & = (ℜ ({\tilde{P}}_{f}^{P}) + i ℑ ({\tilde{P}}_{f}^{P})) \exp i Ω t, \end{matrix}

where $ℜ ({\tilde{p}}^{Q})$ , $ℜ ({\tilde{u}}^{M})$ , $ℑ ({\tilde{p}}^{Q})$ , and $ℑ ({\tilde{u}}^{M})$ are the real and imaginary parts of the amplitudes of the response; $ℜ ({\tilde{P}}^{N})$ and $ℑ ({\tilde{P}}^{N})$ are the real and imaginary parts of the amplitude of the force applied to the structure; $ℜ ({\tilde{P}}_{f}^{P})$ and $ℑ ({\tilde{P}}_{f}^{P})$ are the real and imaginary parts of the amplitude of the acoustic traction (dimensions of volumetric acceleration) applied to the fluid; and $Ω$ is the circular frequency. We substitute these equations into Equation 23 and use the time-harmonic form of Equation 16, $δ {\hat{p}}^{P} = - Ω^{- 2} δ p^{P}$ , which yields the coupled complex linear equation system

\begin{matrix} (24) & [\begin{matrix} ℜ [A_{f}^{P Q}] & ℑ [A_{f}^{P Q}] & [S_{fs}^{P M}] & 0 \\ ℑ [A_{f}^{P Q}] & - ℜ [A_{f}^{P Q}] & 0 & - [S_{fs}^{P M}] \\ {[S_{fs}^{Q N}]}^{T} & 0 & ℜ [A_{s}^{N M}] & ℑ [A_{s}^{N M}] \\ 0 & - {[S_{fs}^{Q N}]}^{T} & ℑ [A_{s}^{N M}] & - ℜ [A_{s}^{N M}] \end{matrix}] {\begin{matrix} ℜ ({\tilde{p}}^{Q}) \\ ℑ ({\tilde{p}}^{Q}) \\ ℜ ({\tilde{u}}^{M}) \\ ℑ ({\tilde{u}}^{M}) \end{matrix}} = {\begin{matrix} Ω^{- 2} ℜ ({\tilde{P}}_{f}^{P}) \\ - Ω^{- 2} ℑ ({\tilde{P}}_{f}^{P}) \\ ℜ ({\tilde{P}}^{N}) \\ - ℑ ({\tilde{P}}^{N}) \end{matrix}}, \end{matrix}

where

\begin{matrix} ℜ [A_{f}^{P Q}] & = Ω^{- 2} (K_{f}^{P Q} + K_{fr}^{P Q}) - (M_{f}^{P Q} + M_{fr}^{P Q}) \\ ℑ [A_{f}^{P Q}] & = - Ω^{- 1} (C_{f}^{P Q} + C_{fr}^{P Q}) - Ω^{- 2} C_{f (s)}^{P Q} \end{matrix}

and

\begin{matrix} ℜ [A_{s}^{N M}] & = K^{N M} - Ω^{2} M^{N M} \\ ℑ [A_{s}^{N M}] & = - Ω (C_{(m)}^{N M} + C_{(k)}^{N M}) - C_{(s)}^{N M} . \end{matrix}

If $K^{N M}$ is symmetric, Equation 24 is symmetric. The system may be quite large, because the real and imaginary parts of the structural degrees of freedom and of the pressure in the fluid all appear in the system. This set of equations is solved for each frequency requested in the direct-solution steady-state dynamics procedure. If damping is absent, the user can specify that only the real matrix equation should be factored in the analysis. Nonzero volumetric drag values ( $γ$ ) for the acoustic medium and nonzero $1 / c_{1}$ values for the impedances represent damping. As mentioned above for the transient case, the coupled system can be split into an uncoupled structural analysis and an acoustic analysis driven by the structural response, provided the fluid forces on the structure are small.

Formulation for eigenvalue extraction and mode-based procedures

From the discretized equation, Equation 17, we can write the frequency domain problem as

[\begin{matrix} K_{s} & S_{fs}^{T} \\ 0 & K_{f} \end{matrix}] {\begin{matrix} u \\ p \end{matrix}} + i ω [\begin{matrix} C_{(m)} & 0 \\ 0 & C_{f} \end{matrix}] {\begin{matrix} u \\ p \end{matrix}} - ω^{2} [\begin{matrix} M_{s} & 0 \\ - S_{fs} & M_{f} \end{matrix}] {\begin{matrix} u \\ p \end{matrix}} = {\begin{matrix} P_{s} \\ P_{s} \end{matrix}},

where $ω$ is a natural (as opposed to forced response) frequency. The indices have been suppressed for brevity. This system is due to Zienkiewicz and Newton (1969) and is used in Abaqus as the starting point for mode-based procedures. Suppressing any damping terms, forcing, and any terms associated with a reactive surface,

\begin{matrix} (25) & [\begin{matrix} K_{s} & S_{fs}^{T} \\ 0 & K_{f} \end{matrix}] {\begin{matrix} u \\ p \end{matrix}} - ω^{2} [\begin{matrix} M_{s} & 0 \\ - S_{fs} & M_{f} \end{matrix}] {\begin{matrix} u \\ p \end{matrix}} = 0. \end{matrix}

Interpreted as a linear eigenvalue problem (where $ω^{2}$ is the eigenvalue), this equation cannot be solved directly in Abaqus due to the unsymmetric stiffness and mass matrices. However, it can be shown that these equations do yield real-valued natural frequencies and modes, suggesting that they can be rewritten in symmetric forms.

Application of the modes of Equation 25 to form a reduced system (see below) must be done with some caution, since this unsymmetric system has distinct left and right eigenvector sets. In particular, the “singular modes” associated with zero frequency are of interest because they describe the low-frequency limiting behavior of the system (or the “rigid-body motion” in a kinematic sense) and are, therefore, essential for the construction of an accurate projected frequency domain operator. The right singular modes of the coupled system are

ϕ_{0}^{R, S} \equiv (K_{s} u_{0}^{R, S} = 0, p_{0}^{R, S} = 0),

ϕ_{0}^{R, A} \equiv (K_{s} u_{0}^{R, A} = - S_{fs}^{T} p_{0}^{R, A}, K_{f} p_{0}^{R, A} = 0) .

In other words, there is a “structural” singular right mode $ϕ_{0}^{R, S}$ associated with the kernel of $K_{s}$ and an “acoustic” singular right mode $ϕ_{0}^{R, A}$ associated with the kernel of $K_{f}$ . The left singular modes are solutions to

ϕ_{0}^{T} [\begin{matrix} K_{s} & S_{fs}^{T} \\ 0 & K_{f} \end{matrix}] = 0

and are

ϕ_{0}^{L, S} \equiv (K_{s} u_{0}^{L, S} = 0, K_{f} p_{0}^{L, S} = - S_{fs} u_{0}^{L, S}),

ϕ_{0}^{L, A} \equiv (u_{0}^{L, A} = 0, K_{f} p_{0}^{L, A} = 0) .

The right acoustic and left structural singular modes are coupled, with nontrivial fields on the structural and acoustic domains. These coupled singular modes are a consequence of the stiffness term in Equation 25, and they must be computed if this system is to be projected.

An alternative frequency domain formulation, due to Everstine (1981), involves the substitution $p \mapsto i ω q$ and results in a formally symmetric system:

[\begin{matrix} K_{s} & 0 \\ 0 & K_{f} \end{matrix}] {\begin{matrix} u \\ q \end{matrix}} + i ω [\begin{matrix} C_{(m)} & S_{fs}^{T} \\ S_{fs} & C_{f} \end{matrix}] {\begin{matrix} u \\ q \end{matrix}} - ω^{2} [\begin{matrix} M_{s} & 0 \\ 0 & M_{f} \end{matrix}] {\begin{matrix} u \\ q \end{matrix}} = {\begin{matrix} P_{s} \\ \frac{P_{s}}{i ω} \end{matrix}} .

The corresponding coupled eigenproblem is quadratic, but the singular mode structure of this system is much simpler—the left and right pairs are identical due to symmetry, and the singular modes are uncoupled due to the diagonal structure of the stiffness matrix. The modes are simply

ϕ_{0}^{S} \equiv (K_{s} u_{0}^{S} = 0, p_{0}^{S} = 0),

ϕ_{0}^{A} \equiv (u_{0}^{A} = 0, K_{f} p_{0}^{A} = 0) .

Lanczos formulation

Introducing an auxiliary variable, $ψ = p / ω^{2}$ , augmenting the system of equations with $K_{f} p = ω^{2} K_{f} ψ$ , and manipulating the equations yields

\begin{matrix} (26) & ([\begin{matrix} K_{s} & S_{fs}^{T} & 0 \\ S_{fs} & - M_{f} & K_{f} \\ 0 & K_{f}^{T} & 0 \end{matrix}] - ω^{2} [\begin{matrix} M_{s} & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & K_{f} \end{matrix}]) {\begin{matrix} u \\ p \\ ψ \end{matrix}} = 0. \end{matrix}

This augmented system of equations is due to Ohayon and is used only for Lanczos eigenvalue extraction. The auxiliary variable $ψ$ is internal to Abaqus/Standard and is not available for output. If $K_{f}$ is singular, orthogonalization against the singular acoustic modes is done in the Lanczos eigensolver.

The left and right eigenvectors for the original system of equations, Equation 25, can be constructed from the Lanczos solution. As mentioned above, the singular modes are essential for construction of an accurate projected operator. It is easy to verify that the Lanczos system has the following structural singular mode:

ϕ_{0}^{S} \equiv (K_{s} u_{0}^{S} = 0, p_{0}^{S} = 0, K_{f} ψ_{0}^{S} = - S_{fs} u_{0}^{S}) .

However, if we seek nontrivial acoustic singular modes (i.e., $p_{0}^{A} \neq 0$ , such that $K_{f} p_{0}^{A} = 0$ ), we easily find that $K_{s} u_{0}^{A} = - S_{fs}^{T} p_{0}^{A}$ but also that

K_{f} ψ_{0}^{A} = M_{f} p_{0}^{A} - S_{fs} u_{0}^{A} .

If a nontrivial $p_{0}^{A}$ exists, $K_{f}$ is singular; therefore, for a solution $ψ_{0}^{A}$ to exist, the right-hand-side must be orthogonal to the null space of $K_{f}$ . But we quickly observe that

p_{0}^{A^{T}} (M_{f} p_{0}^{A} - S_{fs} u_{0}^{A}) = p_{0}^{A^{T}} M_{f} p_{0}^{A} + p_{0}^{A^{T}} S_{fs}^{T} K_{s}^{- 1} S_{fs} p_{0}^{A} > 0 \forall p_{0}^{A} .

Consequently, to find an acoustic singular mode using the Lanczos formulation, we construct a perturbation “force” $f_{null}$ such that $f_{null}^{T} p_{0}^{A} \neq 0.$ The Lanczos formulation will yield the nontrivial singular acoustic mode

ϕ_{0}^{A} \equiv (K_{s}^{- 1} u_{0}^{A} = - S_{fs}^{T} p_{0}^{A}, K_{f} p_{0}^{A} = 0, K_{f} ψ_{0}^{A} = M_{f} p_{0}^{A} - S_{fs} u_{0}^{A} + f_{null}) .

The left and right eigenvectors of the original, unsymmetric system Equation 25, including the singular modes, can be constructed from the Lanczos solutions $(u, p, ψ)$ :

ϕ^{L} = {\begin{matrix} v \\ ψ \end{matrix}}; ϕ^{R} = {\begin{matrix} u \\ p \end{matrix}},

where

v = {\begin{matrix} 0, if
  a singular acoustic
  mode \\ u, otherwise \end{matrix}} .

For any nonsingular acoustic mode $α$ , $p_{α} = ω_{α}^{2} ψ_{α}$ , where $ω_{α}$ is the circular eigenfrequency. The left and right eigenvector subspaces are then used to compute modal quantities (generalized mass, participation factors, and effective mass) and to project the mass, stiffness, and damping matrices in mode-based procedures (such as subspace-based steady-state dynamic analysis or transient modal dynamic analysis) to obtain a reduced system of equations. Most of these computations are conducted in a very similar fashion to the way they are carried out in the pure structural problem and will not be discussed here. In addition, for each mode an acoustic fraction of the generalized mass is computed as the ratio between acoustic contributions to the generalized mass and to the total generalized mass.

The only exception worth a brief discussion is the choice for the calculation of the acoustic participation factors and effective masses, as follows. First, a “rigid body” acoustic mode, $T_{a}$ , analogous to the rigid body modes for the structural problem outlined in Variables associated with the natural modes of a model, is chosen to be a constant pressure field of unity. A total “acoustic mass” is then defined as $m_{ac} = T_{a}^{T} M_{f} T_{a}$ . Left and right acoustic participation factors are defined as

Γ_{α}^{L} = \frac{1}{m_{α}} {[\begin{matrix} v^{M} \\ ψ^{Q} \end{matrix}]}^{T} [\begin{matrix} M^{N M} & 0 \\ - S_{fs} & M_{f} \end{matrix}] [\begin{matrix} 0 \\ T_{a} \end{matrix}] = \frac{1}{m_{α}} ψ^{Q} M_{f} T_{a}

and

Γ_{α}^{R} = \frac{1}{m_{α}} {[\begin{matrix} 0 \\ T_{a} \end{matrix}]}^{T} [\begin{matrix} M^{N M} & 0 \\ - S_{fs} & M_{f} \end{matrix}] [\begin{matrix} u^{M} \\ p^{Q} \end{matrix}] = \frac{1}{m_{α}} (- T_{a} S_{fs} u^{M} + T_{a} M_{f} p^{Q}) .

Abaqus/Standard will then report the acoustic participation factor computed as

Γ_{α}^{ac} = \sqrt{Γ_{α}^{L} Γ_{α}^{R} / m_{ac}}

and an acoustic effective mass computed as

m_{α ac}^{eff} = {(Γ_{α}^{ac})}^{2} m_{α} .

The scaling by $m_{ac}$ in the equation for $Γ_{α}^{ac}$ is arbitrary. However, this scaling ensures that when all eigenmodes are extracted, the sum of all acoustic effective masses is 1.0 (minus the contributions from nodes constrained in the acoustic degree of freedom).

Frequency-domain solution using projections onto modal spaces

Distinct modal space projection methods for coupled forced structural-acoustic response exist in Abaqus for the following cases: using coupled modes from Lanczos, using uncoupled modes from Lanczos, and using uncoupled modes from Abaqus/AMS. In the Lanczos mode cases the forced response is computed using the Zienkiewicz-Newton equation, with separate right and left projection operators. In the Abaqus/AMS uncoupled mode case the Everstine equation is used for the forced response and the right and left projection operators are identical. This case is described in more detail below.

Using uncoupled Abaqus/AMS modes

In this case the Everstine equation is used for the coupled forced response problem and modes are computed from decoupled structural and acoustic Abaqus/AMS runs. In nodal degrees of freedom the forced response is governed by

[\begin{matrix} K_{s} & 0 \\ 0 & K_{f} \end{matrix}] {\begin{matrix} u \\ q \end{matrix}} + i [\begin{matrix} C_{s} & Ω S_{fs}^{T} \\ Ω S_{fs} & C_{f} \end{matrix}] {\begin{matrix} u \\ q \end{matrix}} - Ω^{2} [\begin{matrix} M_{s} & 0 \\ 0 & M_{f} \end{matrix}] {\begin{matrix} u \\ q \end{matrix}} = {\begin{matrix} P_{s} \\ \frac{P_{f}}{i Ω} \end{matrix}},

where $C_{s}$ and $C_{f}$ here are the complete assembled damping matrices for the structure and fluid, including viscous and structural damping, as well as boundary impedance effects. Using transformations constructed from the acoustic and structural modes,

{ϕ_{i}^{S}} \mapsto Φ_{S},

{ϕ_{j}^{A}} \mapsto Φ_{A},

and representations of the structural and acoustic fields in the spaces spanned by these modes,

u \mapsto Φ_{S} (\tilde{u} + i \tilde{v}),

q \mapsto Φ_{A} (\tilde{q} + i \tilde{r}),

we obtain

[\begin{matrix} k_{s} - Ω^{2} m_{s} & Ω c_{s} & 0 & - Ω s^{T} \\ Ω c_{s} & - k_{s} + Ω^{2} m_{s} & - Ω s^{T} & 0 \\ 0 & - Ω s_{fs} & k_{f} - Ω^{2} m_{f} & Ω c_{f} \\ - Ω s & 0 & Ω c_{f} & - k_{f} + Ω^{2} m_{f} \end{matrix}] {\begin{matrix} \tilde{u} \\ \tilde{v} \\ \tilde{q} \\ \tilde{r} \end{matrix}} = {\begin{matrix} ℜ ({\tilde{f}}_{s}) \\ - ℑ ({\tilde{f}}_{s}) \\ ℜ ({\tilde{f}}_{f}) \\ - ℑ ({\tilde{f}}_{f}) \end{matrix}} .

The terms in this matrix correspond to the nodal degree-of-freedom operators, projected onto the modal spaces. The damping and coupling matrices in modal coordinates are full and unsymmetric.

Volumetric drag and fluid viscosity

The medium supporting acoustic waves may be flowing through a porous matrix, such as fiberglass used for sound deadening. In this case the parameter $γ$ is the flow resistance, the pressure drop required to force a unit flow through the porous matrix. A propagating plane wave with nominal particle velocity ${\dot{u}}^{f}$ loses energy at a rate

\begin{matrix} (27) & \dot{E} = - γ {[{\dot{u}}^{f}]}^{2} . \end{matrix}

Fluids also exhibit momentum losses without a porous matrix resistive medium through coefficients of shear viscosity $μ$ and bulk viscosity $η$ . These are proportionality constants between components of the stress and spatial derivatives of the shear strain rate and volumetric strain rate, respectively. In fluid mechanics the shear viscosity term $μ$ is usually more important than the bulk term $η$ ; however, acoustics is the study of volumetrically straining flows, so both constants can be important. The linearized Navier-Stokes equations for adiabatic perturbations about a base state can be expressed in terms of the pressure field alone (Morse and Ingard, 1968):

\begin{matrix} (28) & \frac{\partial}{\partial x} \cdot \frac{\partial p}{\partial x} = \frac{ρ_{f}}{K_{f}} \ddot{p} - \frac{η + \frac{4}{3} μ}{K_{f}} \frac{\partial}{\partial x} \cdot \frac{\partial \dot{p}}{\partial x} . \end{matrix}

In steady state this linearized equation can be written in the form of Equation 19, with

ℜ ({\tilde{ρ}}_{f}) = [\frac{1}{1 + {(\frac{Ω}{K_{f}})}^{2} {(η + \frac{4}{3} μ)}^{2}}] ρ_{f},

ℑ ({\tilde{ρ}}_{f}) = - [\frac{(\frac{Ω}{K_{f}}) (η + \frac{4}{3} μ)}{1 + {(\frac{Ω}{K_{f}})}^{2} {(η + \frac{4}{3} μ)}^{2}}] ρ_{f},

so that the viscosity effects can be modeled as a volumetric drag parameter with the value

γ = [\frac{{(\frac{Ω}{K_{f}})}^{2} (η + \frac{4}{3} μ)}{1 + {(\frac{Ω}{K_{f}})}^{2} {(η + \frac{4}{3} μ)}^{2}}] ρ_{f} K_{f} .

If the combined viscosity effects are small,

\begin{matrix} (29) & \frac{\partial}{\partial x} \cdot \frac{\partial p}{\partial x} \approx \frac{ρ_{f}}{K_{f}} \ddot{p}, \end{matrix}

so we can write

\begin{matrix} (30) & K_{f} \frac{\partial}{\partial x} \cdot \frac{\partial p}{\partial x} - ρ_{f} \frac{\partial^{2}}{\partial t^{2}} p + (η + \frac{4}{3} μ) (\frac{ρ_{f}}{K_{f}}) \frac{\partial^{3}}{\partial t^{3}} p = 0. \end{matrix}

In steady-state form

\begin{matrix} (31) & \frac{1}{ρ_{f}} \frac{\partial}{\partial x} \cdot \frac{\partial p}{\partial x} + Ω^{2} (\frac{1}{K_{f}} - i Ω \frac{η + \frac{4}{3} μ}{K_{f}^{2}}) p = 0, \end{matrix}

where $Ω$ is the forcing frequency. This leads to the following analogy between viscous fluid losses and volumetric drag or flow resistance:

\begin{matrix} (32) & γ = \frac{Ω^{2} ρ_{f}}{K_{f}} (η + \frac{4}{3} μ), \end{matrix}

with density constant with respect to frequency. The energy loss rate for a propagating plane wave in this linearized, adiabatic, small-viscosity case is

\begin{matrix} (33) & \dot{E} = - (η + \frac{4}{3} μ) \frac{Ω^{2} ρ_{f}}{K_{f}} {[{\dot{u}}^{f}]}^{2} . \end{matrix}

Acoustic output quantities

Several secondary quantities are useful in acoustic analysis, derived from the fundamental acoustic pressure field variable. In steady-state dynamics the acoustic particle velocity at any field point is

\tilde{v} = \frac{i}{\tilde{ρ} Ω} \frac{\partial \tilde{p}}{\partial x} .

The acoustic intensity vector, a measure of the rate of flow of energy at a point, is

\tilde{I} = - \frac{1}{2} σ \cdot \hat{\tilde{v}} .

In an acoustic medium the stress tensor is simply the acoustic pressure times the identity tensor, $σ = - p I$ , so this expression simplifies to

\tilde{I} = \frac{- 1}{2 i \hat{\tilde{ρ}} Ω} [\tilde{p} \frac{\partial \hat{\tilde{p}}}{\partial x}] .

The hats denote complex conjugation. The real part of the intensity is referred to as the “active intensity,” and the imaginary part is the “reactive intensity.”

Acoustic contribution factors

Acoustic contribution factors help the user interpret the behavior of a coupled structural-acoustic system by showing the relationship between the acoustic pressure and either specific structural surfaces or specific structural modes. In the literature they are sometimes referred to as acoustic “participation factors,” but since this term is used in Abaqus to describe characteristics of modes (see Variables associated with the natural modes of a model), a different nomenclature is chosen here.

First, consider an acoustic medium in contact with a structure undergoing time-harmonic vibration. The structure exerts a traction on the fluid at each point on the wetted surface, causing harmonic pressure in the acoustic medium. In a given solution to a coupled forced response problem, it is sometimes useful to separate the pressure into constituent parts, each due to the vibration of a portion of the wetted surface. For example, in an automotive acoustic problem it may be useful to determine the parts of the acoustic pressure field due to the windows, floor, and other panels separately. The pressure field $π_{J}$ generated by some given structural vibration acting only on structural surface $J$ , with all other parts of the wetted surface held fixed, is defined as the acoustic contribution factor of that surface:

\begin{matrix} (34) & A_{f} π_{J} \equiv - Ω^{2} S_{J} u, \end{matrix}

where $A_{f} \equiv K_{f} + i Ω C_{f} - Ω^{2} M_{f}$ and $S_{J}$ is the coupling matrix associated with surface partition $J$ . $S_{J}$ can correspond to a group of disjoint surfaces (for example, all the window glass in an automobile) or to a single node. Because the natural boundary condition in Abaqus for acoustic elements is a rigid wall, Equation 34 corresponds physically to an acoustic field coupled to the structure only at surface $J$ , with all other bounding surfaces rigid.

For example, if a single panel's acoustic contribution is separated from the total acoustic pressure,

p = p_{0} + π_{1},

the coupled system of equations for the structural acoustic problem can be written

\begin{matrix} (35) & [\begin{matrix} A_{s} & S_{fs}^{T} & S_{1}^{T} \\ Ω^{2} S_{fs} & A_{f} & A_{f} \\ Ω^{2} S_{1} & 0 & A_{f} \end{matrix}] (\begin{matrix} u \\ p_{0} \\ π_{1} \end{matrix}) = (\begin{matrix} P \\ P_{f} \\ 0 \end{matrix}), \end{matrix}

where $A_{s} \equiv K_{s} + i Ω C_{s} - Ω^{2} M_{s}$ . This equation makes it clear that the panel's acoustic contribution factor depends on the solution to the specific coupled harmonic forced response problem. However, it is more efficient to solve for $p$ and $u$ instead and then solve for $π_{1}$ using Equation 34.

When subspace-based steady-state dynamics or mode-based steady-state dynamics is used, $A_{f}$ and $S_{J}$ are projected; in turn, these projected matrices depend on whether the preceding eigenanalysis step was coupled or uncoupled. For the uncoupled case separate modal transformations $Φ_{f}$ and $Φ_{s}$ correspond to the acoustic and structural modes, and

π_{J} \equiv Φ_{f} {\tilde{π}}_{J},

u \equiv Φ_{s} \tilde{u} .

The transformed equation defining $π_{J}$ becomes

Φ_{f}^{T} A_{f} Φ_{f} {\tilde{π}}_{J} = - Ω^{2} Φ_{f}^{T} S_{J} Φ_{s} \tilde{u} .

The contribution of a specific mode to the acoustic pressure of a forced harmonic coupled system may be of interest as well. Physically, a modal acoustic contribution factor is the part of the acoustic field in a forced response problem due to the action of one structural (or coupled) mode on the acoustic fluid. The calculation of a modal acoustic contribution factor depends on whether the modes in question are uncoupled or coupled structural-acoustic modes. However, its definition is analogous to the surface or panel acoustic contribution factor: it is the acoustic response due to forcing on the wetted surface due only to a single mode of interest, with all other modes held fixed. Starting from Equation 34, but using the entire wetted surface coupling operator $S$ ,

\begin{matrix} (36) & A_{f} π_{J} \equiv - Ω^{2} S u_{J}, \end{matrix}

where $u_{J}$ is the structural response of the coupled problem, restricted to mode $J$ . If coupled mode transformations are used, this equation becomes

\begin{matrix} (37) & {\hat{Φ}}_{f}^{T} A_{f} {\hat{Φ}}_{f} {\hat{π}}_{J} = - Ω^{2} {\hat{Φ}}_{f}^{T} S {\hat{Φ}}_{s} {\hat{u}}_{J} . \end{matrix}

If there is no acoustic force $P_{f}$ in the coupled system of interest and no damping or boundary impedances in the acoustic fluid, this equation is simply the Jth row of the acoustic part of the projected coupled harmonic forced response problem. Consequently, the modal acoustic contribution ${\hat{π}}_{J}$ due to mode J is simply equal to the Jth modal coefficient of the solution to the coupled problem, ${\hat{u}}_{J}$ , times the Jth column of the pressure partition of the modal transformation, ${\hat{Φ}}_{f}$ . Thus, no additional solution is required to obtain modal acoustic contribution factors when using coupled mode projections if acoustic forcing is absent. If acoustic forcing or damping is present in the coupled response problem defining ${\hat{u}}_{J}$ , Equation 37 must be solved after the solution ${\hat{u}}_{J}$ is obtained.

When uncoupled modes are used in the projection for the solution of a coupled system, there is no direct relationship between acoustic and structural mode shapes. Therefore, application of the uncoupled modal transformations to the harmonic forced response problem does not produce the same trivial result as in the coupled mode case. The system resulting from the application of the separate uncoupled mode transformations $Φ_{f}$ and $Φ_{s}$ to Equation 36 must be solved for the modal coefficients ${\tilde{π}}_{J}$ corresponding to forcing via structural mode $J$ :

Φ_{f}^{T} A_{f} Φ_{f} {\tilde{π}}_{J} = - Ω^{2} Φ_{f}^{T} S Φ_{s} {\tilde{u}}_{J} .

Impedance and admittance at fluid boundaries

Equation 11 (or alternatively Equation 9) can be written in a complex admittance form for steady-state analysis:

\begin{matrix} (38) & n^{-} \cdot ({\dot{u}}^{m} - {\dot{u}}^{f}) = (\frac{1}{c_{1}} + \frac{i Ω}{k_{1}}) p = \frac{1}{Z (Ω)} p = - T (x) {(i Ω)}^{- 1}, \end{matrix}

where we define

\begin{matrix} (39) & \frac{1}{Z (Ω)} \equiv \frac{1}{c_{1}} + \frac{i Ω}{k_{1}} . \end{matrix}

The term $1 / Z (Ω)$ is the complex admittance of the boundary, and $Z (Ω)$ is the corresponding complex impedance. Thus, a required complex impedance or admittance value can be entered for a given frequency by fitting data to the parameters $1 / c_{1}$ and $1 / k_{1}$ using Equation 39.

For absorption of plane waves in an infinite medium with volumetric drag, the complex impedance can be shown to be

\begin{matrix} (40) & Z (Ω) = \sqrt{K_{f} \tilde{ρ}} = \sqrt{K_{f} (ρ_{f} + \frac{γ}{i Ω})} . \end{matrix}

For the impedance-based nonreflective boundary condition in Abaqus/Standard, the equations above are used to determine the required constants $1 / c_{1}$ and $1 / k_{1}$ . They are a function of frequency if the volumetric drag is nonzero. The small-drag versions of these equations are used in the direct time integration procedures, as in Equation 46.

Radiation boundary conditions

Many acoustic studies involve a vibrating structure in an infinite domain. In these cases we model a layer of the acoustic medium using finite elements, to a thickness of $1 / 3$ to a full wavelength, out to a “radiating” boundary surface. We then impose a condition on this surface to allow the acoustic waves to pass through and not reflect back into the computational domain. For radiation boundaries of simple shapes—such as planes, spheres, and the like—simple impedance boundary conditions can represent good approximations to the exact radiation conditions. In particular, we include local algebraic radiation conditions of the form

\begin{matrix} (41) & n^{-} \cdot \frac{\partial p}{\partial x} = M p = f (i \tilde{k} + β) p, \end{matrix}

where $\tilde{k} = Ω \sqrt{\tilde{ρ} / K_{f}}$ is the wave number and $\tilde{ρ}$ is the complex density (see Equation 18). f is a geometric factor related to the metric factors of the curvilinear coordinate system used on the boundary, and $β$ is a spreading loss term (see Table 1).

Table 1. Boundary condition parameters.
Geometry	f	$β$
Plane	1	0
Circle or circular cylinder	1	$\frac{1}{2 r_{1}}$
Ellipse or elliptical cylinder	$\sqrt{\frac{1 - ϵ^{2}}{1 - ϵ^{2} {[(x_{g} - x_{0}) \cdot e_{m} / r_{1}]}^{2}}}$	$\frac{1}{2 r_{1}}$
Sphere	1	$\frac{1}{r_{1}}$
Prolate spheroid	$\sqrt{\frac{1 - ϵ^{2}}{1 - ϵ^{2} {[(x_{g} - x_{0}) \cdot e_{m} / r_{1}]}^{2}}}$	$\frac{1}{r_{1}}$

Comparison of Equation 41 and Equation 9 reveals that, for steady-state analysis, there exists a direct analogy to the reactive boundary equation, Equation 21, with

\begin{matrix} (42) & \frac{1}{k_{1}} = ℑ (\frac{f}{Ω \sqrt{\tilde{ρ} K_{f}}}) - \frac{f β}{Ω^{2} ρ_{f} (1 + {(γ / Ω ρ_{f})}^{2})}, \end{matrix}

and

\begin{matrix} (43) & \frac{1}{c_{1}} = ℜ (\frac{f}{\sqrt{\tilde{ρ} K_{f}}}) + \frac{f β r / ρ_{f}}{Ω^{2} ρ_{f} (1 + {(γ / Ω ρ_{f})}^{2})} . \end{matrix}

For transient procedures the treatment of volumetric drag in the acoustic equations and the radiation conditions necessitates an approximation. In the acoustics equation we use the boundary term

\begin{matrix} (44) & - n^{-} \cdot \frac{\partial p}{\partial x} \frac{1}{ρ_{f}} = n^{-} \cdot ({\ddot{u}}^{f} + \frac{γ}{ρ_{f}} {\dot{u}}^{f}) . \end{matrix}

Combining Equation 41 with Equation 44, expanding about $γ = 0$ , and retaining only first-order terms leads to

\begin{matrix} (45) & n^{-} \cdot \frac{\partial p}{\partial x} \frac{1}{ρ_{f}} = f [\frac{i Ω}{\sqrt{ρ_{f} K_{f}}}] p + f [\frac{β}{ρ_{f}} + (\frac{γ}{2 ρ_{f} \sqrt{ρ_{f} K_{f}}})] p . \end{matrix}

The Fourier inverse of the steady-state form results in the transient boundary condition

\begin{matrix} (46) & n^{-} \cdot \frac{\partial p}{\partial x} \frac{1}{ρ_{f}} = f [\frac{1}{\sqrt{ρ_{f} K_{f}}}] \dot{p} + f [\frac{β}{ρ_{f}} + \frac{γ}{2 ρ_{f} \sqrt{ρ_{f} K_{f}}}] p . \end{matrix}

This expression involves independent coefficients for pressure and its first derivative in time, unlike the transient reactive boundary expression (Equation 10), which includes independent coefficients for the first and second derivatives of pressure only. Consequently, to implement this expression, we define the admittance parameters

\begin{matrix} (47) & \frac{1}{c_{1}} = [\frac{f}{\sqrt{ρ_{f} K_{f}}}] \end{matrix}

and

\begin{matrix} (48) & \frac{1}{a_{1}} = f [\frac{β}{ρ_{f}} + \frac{γ}{2 ρ_{f} \sqrt{ρ_{f} K_{f}}}], \end{matrix}

so the boundary traction for the transient radiation boundary condition can be written

n^{-} \cdot \frac{\partial p}{\partial x} \frac{1}{ρ_{f}} = \frac{1}{c_{1}} \dot{p} + \frac{1}{a_{1}} p .

The values of the parameters f and $β$ vary with the geometry of the boundary of the radiating surface of the acoustic medium. The geometries supported in Abaqus are summarized in Table 1. In the table $ϵ$ refers to the eccentricity of the ellipse or spheroid; $r_{1}$ refers to the radius of the circle, sphere, or the semimajor axis of the ellipse or spheroid; $x_{g}$ is the vector locating the integration point on the ellipse or spheroid; $x_{0}$ is the vector locating the center of the ellipse or spheroid; and $e_{m}$ is the vector that orients the major axis.

These algebraic boundary conditions are approximations to the exact impedance of a boundary radiating into an infinite exterior. The plane wave condition is the exact impedance for plane waves normally incident to a planar boundary. The spherical condition exactly annihilates the first Legendre mode of a radiating spherical surface; the circular condition is asymptotically correct for the first mode (Bayliss et al., 1982). The elliptical and prolate spheroidal conditions are based on expansions of elliptical and prolate spheroidal wave functions in the low-eccentricity limit (Grote and Keller, 1995); the prolate spheroidal condition exactly annihilates the first term of its expansion, while the elliptical condition is asymptotic.

An improvement on radiation boundary conditions for plane waves

As already pointed out, the radiation boundary conditions derived in the previous section for plane waves are actually based on the presumption that the sound wave impinges on the boundary from an orthogonal direction. But this is not always the case. Figure 1 shows a general example for plane waves in which the sound wave direction differs from the boundary normal by an angle of $θ$ .

Figure 1. A plane wave not normally incident to the boundary.

To consider this situation accurately, we adopt the plane-wave radiation equation used in Sandler (1998); i.e.,

\begin{matrix} (49) & n^{-} \cdot \frac{\partial p}{\partial x} = \frac{1}{C_{n}} \dot{p} = \frac{\cos θ}{C} \dot{p}, \end{matrix}

where $C = \sqrt{K_{f} / \tilde{ρ}}$ is the sound speed with $\tilde{ρ} = ρ_{f} + γ p / \dot{p}$ and $C_{n}$ is the corresponding speed normal to the boundary. This exact description of the geometry is the starting point for many developments of approximate absorbing boundary conditions (see, for example, Engquist and Majda, 1977). Thus, we have

\begin{matrix} (50) & n^{-} \cdot \frac{\partial p}{\partial x} \frac{1}{ρ_{f}} = \frac{\cos θ}{\sqrt{ρ_{f} K_{f}}} \sqrt{1 + \frac{γ}{ρ_{f}} \frac{p}{\dot{p}}} \cdot \dot{p} . \end{matrix}

Using the first-order expanding approximation to the second term in the square root in the above equation (similar to what we did to reach Equation 45), we can obtain an improved radiation boundary condition

\begin{matrix} (51) & n^{-} \cdot \frac{\partial p}{\partial x} \frac{1}{ρ_{f}} = \frac{\cos θ}{\sqrt{ρ_{f} K_{f}}} (\dot{p} + \frac{γ}{2 ρ_{f}} p) . \end{matrix}

It can be found from comparison that this equation differs from Equation 46 only by a factor of $\cos θ$ for plane waves. In two dimensions the $\cos θ$ can be calculated as

\begin{matrix} (52) & \cos θ = | \frac{\partial p}{\partial n} | / \sqrt{{(\frac{\partial p}{\partial n})}^{2} + {(\frac{\partial p}{\partial s})}^{2}} . \end{matrix}

The normal and tangential derivatives $\partial p / \partial n$ and $\partial p / \partial s$ at the integration points can be evaluated using the corresponding element along the radiation boundary surface (see Figure 2); i.e.,

\begin{matrix} (53) & \frac{\partial p}{\partial n} = n \cdot \frac{\partial H^{P}}{\partial x} p^{P} and \frac{\partial p}{\partial s} = s \cdot \frac{\partial H^{P}}{\partial x} p^{P}, \end{matrix}

where $p^{P}$ are the nodal pressure values of the element.

Figure 2. An element along the boundary.

The method described in this section can be used only for direct integration transient dynamics; it cannot be used with steady-state or modal response. In addition, it is available for planar, axisymmetric, and three-dimensional geometries.

Finally, the method makes the equilibrium equations nonlinear, as shown in Equation 52. Although in theory the iteration process in Abaqus/Standard can solve the nonlinear equilibrium equations accurately, the use of a small half-increment residual tolerance is strongly suggested since in many cases the pressure and its related residual along the radiation boundaries are very weak relative to the other places in the modeled domain. The computation of $\cos θ$ at the integration point is based on the nodal pressures. The nodal pressures are updated using the explicit central difference procedure described in Explicit dynamic analysis.