<!-- source-page: 201 -->

# Procedures

$$
\begin{array}{l} - \delta \hat {p} ^ {P} \Bigg \{\left[ - \Omega^ {2} (M _ {\mathrm{f}} ^ {P Q} + M _ {\mathrm{fr}} ^ {P Q}) + i \Omega (C _ {\mathrm{f}} ^ {P Q} + C _ {\mathrm{fr}} ^ {P Q}) + K _ {\mathrm{f}} ^ {P Q} \right] \Delta \tilde {p} ^ {Q} + \Omega^ {2} S _ {\mathrm{fs}} ^ {P M} \Delta \tilde {u} ^ {M} - \Delta \tilde {P} _ {f} ^ {P} \Bigg \} \\ + \delta u ^ {N} \Bigg \{\left[ - \Omega^ {2} M ^ {N M} + i \Omega (C _ {(m)} ^ {N M} + C _ {(k)} ^ {N M}) + K ^ {N M} \right] \Delta \tilde {u} ^ {M} + \left[ S _ {\mathrm{fs}} ^ {Q N} \right] ^ {T} \Delta \tilde {p} ^ {Q} - \Delta \tilde {P} ^ {N} \Bigg \} = 0, \\ \end{array}
$$

with

$$
K ^ {N M} = \int_ {V} \left[ \frac {\partial \pmb {\beta} ^ {N}}{\partial u ^ {M}}: \pmb {\sigma} _ {0} + \pmb {\beta} ^ {N}: \mathbf {D} ^ {e l}: \pmb {\beta} ^ {M} \right] 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} \left[ \beta_ {c} \pmb {\beta} ^ {N}: \mathbf {D} ^ {e l}: \pmb {\beta} ^ {M} \right] d V.
$$

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

$$
\Delta \tilde {p} ^ {Q} = \left(\Re \left(\tilde {p} ^ {Q}\right) + i \Im \left(\tilde {p} ^ {Q}\right)\right) \exp i \Omega t
$$

$$
\Delta \tilde {u} ^ {M} = \left(\Re \left(\tilde {u} ^ {M}\right) + i \Im \left(\tilde {u} ^ {M}\right)\right) \exp i \Omega t
$$

and

$$
\Delta \tilde {P} ^ {N} = \left(\Re \left(\tilde {P} ^ {N}\right) + i \Im \left(\tilde {P} ^ {N}\right)\right) \exp i \Omega t
$$

$$
\Delta \tilde {P} _ {f} ^ {P} = \left(\Re \left(\tilde {P} _ {f} ^ {P}\right) + i \Im \left(\tilde {P} _ {f} ^ {P}\right)\right) \exp i \Omega t,
$$

where < $\left( \tilde { p } ^ { Q } \right) , \Re \left( \tilde { u } ^ { M } \right) , \Im \left( \tilde { p } ^ { Q } \right)$ , and $\mathfrak { F } \left( \tilde { u } ^ { M } \right)$ are the real and imaginary parts of the amplitudes of the response; $\Re \left( \tilde { P } ^ { N } \right)$ and $\Im \left( \tilde { P } ^ { N } \right)$ are the real and imaginary parts of the amplitude of the force applied to the structure; $\Re ( \tilde { P } _ { f } ^ { P } )$ and $\Im ( \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 2.9.1-22 and use the time-harmonic form of Equation $2 . 9 . 1 \substack { - 1 - 1 4 , \delta \hat { p } ^ { P } } = - \Omega ^ { - 2 } \delta p ^ { P }$ , which yields the coupled complex linear equation system

Equation 2.9.1-23

$$
\left[ \begin{array}{c c c c} \Re \Big [ A _ {\mathrm{f}} ^ {P Q} \Big ] & \Im \Big [ A _ {\mathrm{f}} ^ {P Q} \Big ] & \Big [ S _ {\mathrm{fs}} ^ {P M} \Big ] & 0 \\ \Im \Big [ A _ {\mathrm{f}} ^ {P Q} \Big ] & - \Re \Big [ A _ {\mathrm{f}} ^ {P Q} \Big ] & 0 & - \Big [ S _ {\mathrm{fs}} ^ {P M} \Big ] \\ \Big [ S _ {\mathrm{fs}} ^ {Q N} \Big ] ^ {T} & 0 & \Re \Big [ A _ {\mathrm{s}} ^ {N M} \Big ] & \Im \Big [ A _ {\mathrm{s}} ^ {N M} \Big ] \\ 0 & - \Big [ S _ {\mathrm{fs}} ^ {Q N} \Big ] ^ {T} & \Im \Big [ A _ {\mathrm{s}} ^ {N M} \Big ] & - \Re \Big [ A _ {\mathrm{s}} ^ {N M} \Big ] \end{array} \right] \left\{ \begin{array}{l} \Re \left(\tilde {p} ^ {Q}\right) \\ \Im \left(\tilde {p} ^ {Q}\right) \\ \Re \left(\tilde {u} ^ {M}\right) \\ \Im \left(\tilde {u} ^ {M}\right) \end{array} \right\} = \left\{ \begin{array}{c} \Omega^ {- 2} \Re \left(\tilde {P} _ {f} ^ {P}\right) \\ - \Omega^ {- 2} \Im \left(\tilde {P} _ {f} ^ {P}\right) \\ \Re \left(\tilde {P} ^ {N}\right) \\ - \Im \left(\tilde {P} ^ {N}\right) \end{array} \right\},
$$

<!-- source-page: 202 -->

where

$$
\Re \Big [ A _ {\mathrm{f}} ^ {P Q} \Big ] = \Omega^ {- 2} (K _ {\mathrm{f}} ^ {P Q} + K _ {\mathrm{fr}} ^ {P Q}) - (M _ {\mathrm{f}} ^ {P Q} + M _ {\mathrm{fr}} ^ {P Q})
$$

$$
\Im \Big [ A _ {\mathrm{f}} ^ {P Q} \Big ] = - \Omega^ {- 1} (C _ {\mathrm{f}} ^ {P Q} + C _ {\mathrm{fr}} ^ {P Q})
$$

and

$$
\Re \Big [ A _ {\mathrm{s}} ^ {N M} \Big ] = K ^ {N M} - \Omega^ {2} M ^ {N M}
$$

$$
\Im \left[ A _ {\mathrm{s}} ^ {N M} \right] = - \Omega \left(C _ {(m)} ^ {N M} + C _ {(k)} ^ {N M}\right).
$$

$\mathbf { I f } K ^ { N M }$ is symmetric, Equation 2.9.1-23 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 \*STEADY STATE DYNAMICS, DIRECT procedure. If damping is absent, the \*STEADY STATE DYNAMICS, DIRECT=REAL ONLY procedure can be used; in this case a smaller, real matrix equation is solved. Nonzero r 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.

# 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 r 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 { \mathrm { ~ \bf ~ u ~ } } ^ { f }$ loses energy at a rate

Equation 2.9.1-24

$$
\dot {E} = - r \left| \dot {\mathbf {u}} ^ {f} \right| ^ {2}.
$$

Fluids also exhibit momentum losses without a porous matrix resistive medium, through coefficients of shear viscosity $\mu$ 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 $\mu$ is usually more important that 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):

$$
\frac {\partial}{\partial \mathbf {x}} \cdot \frac {\partial p}{\partial \mathbf {x}} = \frac {\rho_ {f}}{K _ {f}} \ddot {p} - \frac {\eta + \frac {4}{3} \mu}{K _ {f}} \frac {\partial}{\partial \mathbf {x}} \cdot \frac {\partial \dot {p}}{\partial \mathbf {x}}.
$$

Equation 2.9.1-25

<!-- source-page: 203 -->

If the combined viscosity effects are small,

Equation 2.9.1-26

$$
\frac {\partial}{\partial \mathbf {x}} \cdot \frac {\partial p}{\partial \mathbf {x}} \approx \frac {\rho_ {f}}{K _ {f}} \ddot {p},
$$

so that we can write

Equation 2.9.1-27

$$
K _ {f} \frac {\partial}{\partial \mathbf {x}} \cdot \frac {\partial p}{\partial \mathbf {x}} - \rho_ {f} \frac {\partial^ {2}}{\partial t ^ {2}} p + (\eta + \frac {4}{3} \mu) \left(\frac {\rho_ {f}}{K _ {f}}\right) \frac {\partial^ {3}}{\partial t ^ {3}} p = 0.
$$

This equation involves third-order time derivatives, which we do not solve in transient analyses. However, in steady state we see that

Equation 2.9.1-28

$$
\frac {1}{\rho_ {f}} \frac {\partial}{\partial \mathbf {x}} \cdot \frac {\partial p}{\partial \mathbf {x}} + \Omega^ {2} \left(\frac {1}{K _ {f}} - i \Omega \frac {\eta + \frac {4}{3} \mu}{K _ {f} ^ {2}}\right) p = 0,
$$

where − is the forcing frequency, which leads to the following analogy between viscous fluid losses and volumetric drag or flow resistance:

Equation 2.9.1-29

$$
r = \frac {\Omega^ {2} \rho_ {f}}{K _ {f}} \left(\eta + \frac {4}{3} \mu\right).
$$

The energy loss rate for a propagating plane wave in this linearized, adiabatic, small-viscosity case is

Equation 2.9.1-30

$$
\dot {E} = - \left(\eta + \frac {4}{3} \mu\right) \frac {\Omega^ {2} \rho_ {f}}{K _ {f}} \left| \dot {\bf {u}} ^ {f} \right| ^ {2}.
$$

# Impedance and admittance at fluid boundaries

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

Equation 2.9.1-31

$$
\mathbf {n} ^ {-} \cdot (\dot {\mathbf {u}} ^ {m} - \dot {\mathbf {u}} ^ {f}) = (\frac {1}{c _ {1}} + \frac {i \Omega}{k _ {1}}) p = \frac {1}{Z (\Omega)} p = - T (\mathbf {x}) (i \Omega) ^ {- 1},
$$

where we define

Equation 2.9.1-32

$$
\frac {1}{Z (\Omega)} \equiv \frac {1}{c _ {1}} + \frac {i \Omega}{k _ {1}}.
$$

<!-- source-page: 204 -->

# Procedures

The term $1 / Z ( \Omega )$ is the complex admittance of the boundary, and $Z ( \Omega )$ is the corresponding complex impedance. Thus, a required complex impedance or admittance value can be entered for a given frequency by fitting to the parameters $1 / c _ { 1 }$ and $1 / k _ { 1 }$ using Equation 2.9.1-32.

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

$$
Z (\Omega) = \sqrt {K _ {f} \tilde {\rho}} = \sqrt {K _ {f} (\rho_ {f} + \frac {r}{i \Omega})}.
$$

Equation 2.9.1-33

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 2.9.1-39. For more information, see \`\`Acoustic and coupled acoustic-structural analysis,'' Section 6.9.1 of the ABAQUS/Standard User's Manual.

# 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 / 4$ 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

Equation 2.9.1-34

$$
\mathbf {n} ^ {-} \cdot \frac {\partial p}{\partial \mathbf {x}} = M p = f (i \tilde {k} + \beta) p,
$$

where $\tilde { k } = \Omega \sqrt { \tilde { \rho } / K _ { f } }$ is the wave number, and $\tilde { \rho }$ is the complex density (see Equation 2.9.1-16). f is a geometric factor related to the metric factors of the curvilinear coordinate system used on the boundary, and $\beta$ is a spreading loss term (see Table 2.9.1-1). Comparison of Equation 2.9.1-34 and Equation 2.9.1-9 reveals that, for steady-state analysis, there exists a direct analogy to the reactive boundary equation, Equation 2.9.1-19, with

$$
\frac {1}{k _ {1}} = \Im \bigl (\frac {f}{\Omega \sqrt {\tilde {\rho} K _ {f}}} \bigr) - \frac {f \beta}{\Omega^ {2} \rho_ {f} \left(1 + (r / \Omega \rho_ {f} \right. ^ {2})},
$$

Equation 2.9.1-35

and

$$
\frac {1}{c _ {1}} = \Re \bigl (\frac {f}{\sqrt {\tilde {\rho} K _ {f}}} \bigr) + \frac {f \beta r / \rho_ {f}}{\Omega^ {2} \rho_ {f} (1 + (r / \Omega \rho_ {f}) ^ {2})}.
$$

Equation 2.9.1-36

<!-- source-page: 205 -->

# Procedures

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

Equation 2.9.1-37

$$
- \mathbf {n} ^ {-} \cdot \frac {\partial p}{\partial \mathbf {x}} \frac {1}{\rho_ {f}} = \mathbf {n} ^ {-} \cdot \left(\ddot {\mathbf {u}} ^ {f} + \frac {r}{\rho_ {f}} \dot {\mathbf {u}} ^ {f}\right).
$$

Combining Equation 2.9.1-34 with Equation 2.9.1-37, expanding about r = 0, and retaining only first-order terms leads to

Equation 2.9.1-38

$$
\left(\mathbf {n} ^ {-} \cdot \frac {\partial p}{\partial \mathbf {x}}\right) / \rho_ {f} = f \left[ \frac {i \Omega}{\sqrt {\rho_ {f} K _ {f}}} \right] p + f \left[ \frac {\beta}{\rho_ {f}} - \left(\frac {r}{2 \rho_ {f} \sqrt {\rho_ {f} K _ {f}}}\right) - \frac {\beta r}{i \Omega \rho_ {f} ^ {2}} \right] p.
$$

The Fourier inverse of the steady-state form contains a time convolution term, which is not implemented. Dropping this term, retaining only differential terms, is equivalent to making the physical assumption that the volumetric drag is small compared to $\Omega ^ { 2 } \rho _ { f } \left( \sqrt { \rho _ { f } / K _ { f } } \right) / \beta$ . Since this is a common case, we have implemented the transient boundary condition

Equation 2.9.1-39

$$
\left(\mathbf {n} ^ {-} \cdot \frac {\partial p}{\partial \mathbf {x}}\right) / \rho_ {f} = f \left[ \frac {1}{\sqrt {\rho_ {f} K _ {f}}} \right] \dot {p} + f \left[ \frac {\beta}{\rho_ {f}} - \frac {r}{2 \rho_ {f} \sqrt {\rho_ {f} K _ {f}}} \right] p.
$$

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

Equation 2.9.1-40

$$
\frac {1}{c _ {1}} = \left[ \frac {f}{\sqrt {\rho_ {f} K _ {f}}} \right]
$$

and

Equation 2.9.1-41

$$
\frac {1}{a _ {1}} = f \left[ \frac {\beta}{\rho_ {f}} - \frac {r}{2 \rho_ {f} \sqrt {\rho_ {f} K _ {f}}} \right],
$$

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

$$
- \mathbf {n} ^ {-} \cdot \frac {\partial p}{\partial \mathbf {x}} \frac {1}{\rho_ {f}} = \frac {1}{c _ {1}} \dot {p} + \frac {1}{a _ {1}} p.
$$

Equation -

<!-- source-page: 206 -->

The values of the parameters $f$ and $\beta$ vary with the geometry of the boundary of the radiating surface of the acoustic medium. The geometries supported in ABAQUS/Standard are summarized in Table 2.9.1-1.

Table 2.9.1-1. Boundary condition parameters. 

<table><tr><td>Geometry</td><td>f</td><td>β</td></tr><tr><td>Plane</td><td>1</td><td>0</td></tr><tr><td>Circle or circular</td><td>1</td><td> $\frac{1}{2r_1}$ </td></tr><tr><td>cylinder</td><td></td><td></td></tr><tr><td>Ellipse or elliptical cylinder</td><td><img src="images/59ce17c3b367db7f7716f71bfad72fb42d2ada00132b0986addf55e54e98fa82.jpg"/></td><td> $\frac{1}{2r_1}$ </td></tr><tr><td>Sphere</td><td>1</td><td> $\frac{1}{r_1}$ </td></tr><tr><td>Prolate spheroid</td><td><img src="images/b44b5f7a298099cb1e8c5569d2e9828e3956a643aea894b30ddad580ee087b5e.jpg"/></td><td> $\frac{1}{r_1}$ </td></tr></table>

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; $\mathbf { x } _ { g }$ is the vector locating the integration point on the ellipse or spheroid; $\mathbf { x } _ { 0 }$ is the vector locating the center of the ellipse or spheroid, and $\mathbf { 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.

# 2.9.2 Underwater shock analysis

The underwater shock analysis capability is provided by the coupling of ABAQUS/Standard with Unique Software Applications' USA code (DeRuntz et al., 1980). This coupling is accomplished by a staggered solution method in which ABAQUS/Standard is used to calculate the structural response and USA calculates the fluid pressure response at the interaction surface.

This section provides a summary of equations involved in the staggered solution method. The notation used for the various matrices, vectors, and scalars is somewhat different from the notation used in the other ABAQUS manuals but corresponds closely to that used in the two references of DeRuntz.

# Fluid-surface equations--DAA

The USA solution of fluid-surface interaction is based on the Doubly Asymptotic Approximations (DAA): "Doubly asymptotic approximations are differential equations for the simplified analysis of the transient interaction between a flexible structure and a surrounding infinite medium." ( Geers, 1978) More simply, the DAA are surface interaction approximations. As their name suggests, these approximations approach the exact relationships in the limit of both low- and high-frequency response.

<!-- source-page: 207 -->

# Procedures

The first-order approximation, called $\mathbf { D A A _ { 1 } }$ , is given in classical matrix notation by (DeRuntz, 1989):

Equation 2.9.2-1

$$
[ M _ {f} ] \{\dot {p} _ {S} \} + \rho c [ A _ {f} ] \{p _ {S} \} = \rho c [ M _ {f} ] \{\dot {u} _ {S} \},
$$

where $[ M _ { f } ]$ is the symmetric fluid mass matrix; $\{ p _ { S } \}$ is the scattered-wave pressure vector; $\rho$ and c are the density and sound velocity of the fluid; $\left[ A _ { f } \right]$ is the (diagonal) matrix of fluid areas; and $\{ u _ { S } \}$ is the vector of scattered-wave fluid particle velocities normal to the structure's surface and it is this term that is coupled to the structural response by the following:

Equation 2.9.2-2

$$
\{u _ {S} \} = [ G ] ^ {T} \{\dot {x} \} - \{u _ {I} \}.
$$

$[ G ] ^ { T }$ transforms nodal structural velocities, $\{ \dot { x } \}$ , into normal surface velocities at the centroid of each fluid element, and $\{ u _ { I } \}$ is the fluid incident velocity. The second-order approximation, called $\mathbf { D A A } _ { 2 }$ , is given by (DeRuntz, 1989):

Equation 2.9.2-3

$$
[ M _ {f} ] \{\ddot {p} _ {S} \} + \rho c [ A _ {f} ] \{\dot {p} _ {S} \} + \rho c [ \Omega_ {f} ] [ A _ {f} ] \{p _ {S} \} = \rho c \big ([ M _ {f} ] \{\ddot {u} _ {S} \} + [ \Omega_ {f} ] [ M _ {f} ] \{\dot {u} _ {S} \} \big),
$$

where

Equation 2.9.2-4

$$
[ \Omega_ {f} ] = \eta \rho c [ A _ {f} ] [ M _ {f} ] ^ {- 1}.
$$

The scalar parameter $\eta$ is bounded as $0 \leq \eta \leq 1 . \mathrm { D A A _ { 1 } }$ is recovered by setting $\scriptstyle \eta = 0$ . Typically $\scriptstyle \eta = 1 / 2$ for an infinite cylinder and $\eta { = } 1$ for a sphere. Equation 2.9.2-3 is usually transformed by integrating with respect to time to yield

Equation 2.9.2-5

$$
[ M _ {f} ] \{\ddot {q} _ {S} \} + \rho c [ A _ {f} ] \{\dot {q} _ {S} \} + \rho c [ \Omega_ {f} ] [ A _ {f} ] \{q _ {S} \} = \rho c \big ([ M _ {f} ] \{\dot {u} _ {S} \} + [ \Omega_ {f} ] [ M _ {f} ] \{u _ {S} \} \big),
$$

where

Equation 2.9.2-6

$$
\{q \} = \int_ {0} ^ {t} \{p (\tau) \} d \tau .
$$

# Coupled fluid-structure interaction equations

The DAA equation is used in conjunction with the standard structural response equation of motion,

Equation 2.9.2-7

$$
[ M _ {s} ] \{\ddot {x} \} + [ C _ {s} ] \{\dot {x} \} + [ K _ {s} ] \{x \} = \{f \}.
$$

<!-- source-page: 208 -->

The coupling to the fluid response is through the right-hand-side force. The structural force vector ffg due to fluid loading can be written as

Equation 2.9.2-8

$$
\{f \} = - [ G ] [ A _ {f} ] \big (\{p _ {I} \} + \{p _ {S} \} \big),
$$

where $\{ p _ { I } \}$ is the incident pressure. In the absence of nonfluid loading, the coupled fluid-structure interaction equations for $\mathbf { D A A _ { 1 } }$ are

Equation 2.9.2-9

$$
[ M _ {s} ] \{\ddot {x} \} + [ C _ {s} ] \{\dot {x} \} + [ K _ {s} ] \{x \} = - [ G ] [ A _ {f} ] \big (\{p _ {I} \} + \{p _ {S} \} \big)
$$

Equation 2.9.2-10

$$
[ M _ {f} ] \{\dot {p} _ {S} \} + \rho c [ A _ {f} ] \{p _ {S} \} = \rho c [ M _ {f} ] \big ([ G ] ^ {T} \{\ddot {x} \} - \{\dot {u} _ {I} \} \big).
$$

In the absence of nonfluid loading, the coupled fluid-structure interaction equations for $\mathrm { \ D A A _ { 2 } }$ are

Equation 2.9.2-11

$$
[ M _ {s} ] \{\ddot {x} \} + [ C _ {s} ] \{\dot {x} \} + [ K _ {s} ] \{x \} = - [ G ] [ A _ {f} ] \big (\{p _ {I} \} + \{p _ {S} \} \big)
$$

Equation 2.9.2-12

$$
[ M _ {f} ] \{\ddot {q} _ {S} \} + \rho c [ A _ {f} ] \{\dot {q} _ {S} \} + \rho c [ \Omega_ {f} ] [ A _ {f} ] \{q _ {S} \} =
$$

$$
\rho c \bigl [ [ M _ {f} ] \bigl ([ G ] ^ {T} \{\ddot {x} \} - \{\dot {u} _ {I} \} \bigr) + [ \Omega_ {f} ] [ M _ {f} ] \bigl ([ G ] ^ {T} \{\dot {x} \} - \{u _ {I} \} \bigr) \bigr ].
$$

# The staggered solution scheme

A simple staggered solution scheme would be as follows. Solve Equation 2.9.2-9 and Equation 2.9.2-11 for the structural unknowns, assuming that all right-hand side terms are known (explicitly or by extrapolation). Then all right-hand terms in Equation 2.9.2-10 and Equation 2.9.2-12 are known, and the scattered pressures can be found. Unfortunately this scheme is conditionally stable. Park et al. (1977) stabilized the staggered scheme by a process called augmentation, which involved solving algebraically for the structural accelerations in Equation 2.9.2-9 and Equation 2.9.2-11 and substituting them into Equation 2.9.2-10 and Equation 2.9.2-12. The augmented fluid equation is then solved first, using extrapolated values of the structural vector $[ C _ { s } ] \{ \dot { x } \} + [ K _ { s } ] \{ x \}$ . Then Equation 2.9.2-9 and Equation 2.9.2-11 are solved for the structural unknowns. The same process is used for each time increment. For the DAA2 equation there is an additional structural velocity term that must also be extrapolated.

Additional modifications are usually made to the fluid equation to remove singularities in the incident fluid velocities. The details can be found in DeRuntz (1989). The augmented, coupled fluid-structure

<!-- source-page: 209 -->

interaction equations are summarized below:

Equation 2.9.2-13

$$
\begin{array}{l} \left[ \begin{array}{c c} M _ {s} & 0 \\ 0 & A _ {f} \end{array} \right] \left\{ \begin{array}{c} \ddot {x} \\ \ddot {q} _ {M} \end{array} \right\} \\ + \left[ \begin{array}{c c} C _ {s} & G A _ {f} \\ \rho c A _ {f} G ^ {T} M _ {s} ^ {- 1} C _ {s} - \eta \rho c D _ {f _ {1}} G ^ {T} & D _ {f _ {1}} + D _ {s} \end{array} \right] \left\{ \begin{array}{c} \dot {x} \\ \dot {q} _ {M} \end{array} \right\} \\ + \left[ \begin{array}{c c} K _ {s} & 0 \\ \rho c A _ {f} G ^ {T} M _ {s} ^ {- 1} K _ {s} & \eta D _ {f _ {2}} \end{array} \right] \left\{ \begin{array}{l} x \\ q _ {M} \end{array} \right\} = \left\{ \begin{array}{l} - G A _ {f} (I - \Gamma) p _ {I} \\ - H _ {1} p _ {I} + H _ {2} \stackrel {\star} {p} _ {I} \end{array} \right\}, \\ \end{array}
$$

where

$$
[ D _ {s} ] = \rho c [ A _ {f} ] [ G ] ^ {T} [ M _ {s} ] ^ {- 1} [ G ] [ A _ {f} ]
$$

$$
[ D _ {f _ {1}} ] = \rho c [ A _ {f} ] [ M _ {f} ] ^ {- 1} [ A _ {f} ]
$$

$$
[ D _ {f _ {2}} ] = \rho^ {2} c ^ {2} [ A _ {f} ] [ M _ {f} ] ^ {- 1} [ A _ {f} ] [ M _ {f} ] ^ {- 1} [ A _ {f} ] = \rho c [ D _ {f _ {1}} ] [ M _ {f} ] ^ {- 1} [ A _ {f} ]
$$

$$
[ H _ {1} ] = [ D _ {s} ] - [ D _ {s} + (1 - \eta) D _ {f _ {1}} ] [ \Gamma ]
$$

$$
[ H _ {2} ] = \eta [ D _ {f _ {2}} ] [ \Gamma ],
$$

and where

$$
\{\stackrel {\star} {p} (t) \} = \int_ {0} ^ {t} \{p (\tau) \} d \tau
$$

$$
\{q \} = \{\stackrel {\star} {p} \}
$$

$$
\{\dot {q} \} = \{p \}
$$

$$
\{p \} = \{p _ {S} \} + \{p _ {I} \}
$$

$$
\{p _ {M} \} = \{p _ {S} \} + [ \Gamma ] \{p _ {I} \}
$$

$$
\{p \} = [ I - \Gamma ] \{p _ {I} \} + \{p _ {M} \},
$$

where [¡] is a diagonal matrix of cosines $\gamma _ { i } ,$ , where $\gamma _ { i }$ is the cosine of the angle between the vector from the charge source to the position of the ith fluid element and the normal to the ith fluid element. ´ is the $\mathrm { \ D A A _ { 2 } }$ parameter $( 0 \leq \eta \leq 1 )$ .

The initial conditions are $\{ x ( 0 ) \} = 0 , \{ \dot { x } ( 0 ) \} = 0 , \{ q _ { M } ( 0 ) \} = 0$ , and $\{ { \dot { q } } _ { M } ( 0 ) \} = 0$ . In addition, $\{ { \ddot { q } } _ { M } ( 0 ) \} = 0$ .

# 2.10 Piezoelectric analysis

# 2.10.1 Piezoelectric analysis

The piezoelectrical effect is the coupling of stress and electrical field in a material: an electrical field causes the material to strain, and vice versa. ABAQUS/Standard has the capability to perform fully coupled piezoelectric analysis. The elements that are used in this case contain both displacement degrees of freedom and the electric potential as nodal variables.

<!-- source-page: 210 -->

# Equilibrium and flux conservation

The piezoelectric effect is governed by coupled mechanical equilibrium and electric flux conservation equations.

The mechanical equilibrium equation is

Equation 2.10.1-1

$$
\int_ {V} \boldsymbol {\sigma}: \delta \boldsymbol {\varepsilon} d V = \int_ {S} \mathbf {t} \cdot \delta \mathbf {u} d S + \int_ {V} \mathbf {f} \cdot \delta \mathbf {u} d V,
$$

where $\pmb { \sigma }$ is the "true" (Cauchy) stress at a point currently at $\mathbf { x } ;$ t is the traction across a point of the surface of the body; f is the body force per unit volume in the body (such as the d'Alembert force $\mathbf { f } = - \rho \ddot { \mathbf { u } }$ , in which $\rho$ is the density of the body); and $\delta \pmb { \varepsilon } \overset { \mathrm { d e f } } { = } \mathrm { s y m } ( \partial \delta \mathbf { u } / \partial \mathbf { x } )$ , where ±u is an arbitrary, continuous vector field (the virtual velocity field).

The electrical flux conservation equation is

Equation 2.10.1-2

$$
\int_ {V} \mathbf {q} \cdot \delta \pmb {E} d V + \int_ {S} q _ {S} \delta \varphi d S + \int_ {V} q _ {V} \delta \varphi d V = 0,
$$

where $\mathbf { q }$ is the electric flux vector; $q _ { S }$ is the electric flux per unit area entering the body at a point on its surface; $q _ { V }$ is the electric flux entering the body per unit volume; and $\delta \boldsymbol { \mathbf { { z } } } \overset { \mathrm { d e f } } { = } - \partial \delta \varphi / \partial \mathbf { x }$ , where $\delta \varphi$ is an arbitrary, continuous, scalar field (the virtual potential). These quantities are also known by other terms that are frequently used within electrical engineering references. The electric flux vector $q$ is known as the electrical displacement, and the potential gradient $E$ is known as the electrical field.

# Constitutive behavior: material coupling

Currently the assumption of linear materials is utilized. The basic equations for a piezoelectric linear medium are defined in the following. The mechanical behavior is

$$
\sigma_ {i j} = D _ {i j k l} \varepsilon_ {k l} - D _ {i j k l} d _ {m k l} ^ {\varphi} E _ {m},
$$

which is given in terms of the piezoelectric strain coefficient matrix $d _ { m k l } ^ { \varphi }$ . In terms of the piezoelectric stress coefficient matrix $e _ { m k l } ^ { \varphi }$ , the mechanical behavior is

Equation 2.10.1-3

$$
\sigma_ {i j} = D _ {i j k l} \varepsilon_ {k l} - e _ {m i j} ^ {\varphi} E _ {m}.
$$

The electrical behavior is

$$
q _ {i} = d _ {i j k} ^ {\varphi} D _ {j k l m} \varepsilon_ {l m} + D _ {i j} ^ {\varphi} E _ {j},
$$

which can be given as