Loading due to an incident dilatational wave field

Scattering formulation

The equations for coupled fluid-solid interaction in Abaqus are developed in Coupled acoustic-structural medium analysis. Here, we proceed from the coupled system:

\begin{matrix} (1) & \begin{matrix} (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}] T^{M} + P_{f}^{P} \\ M^{N M} {\ddot{u}}^{M} + C_{(m)}^{N M} {\dot{u}}^{M} + I^{N} & ​ = - {[S_{fs}^{Q N}]}^{T} p^{Q} + P^{N} \end{matrix} \end{matrix}

These equations couple the total pressure in the fluid to the displacements in the solid. They are valid for any combination of fluid and solid domains in a particular model: the matrix $S_{fs}$ is defined on all the interacting fluid and solid surfaces. The fluid traction

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

is a quantity (with dimensions of acceleration) that describes the mechanism by which the solid motion drives the fluid; the fluid drives the solid by the pressure on the solid surface. In Equation 2 $n^{-}$ is the inward normal on the boundary.

To proceed, we decompose the total pressure into the known incident wave component and the unknown scattered component:

p \equiv p_{I} + p_{S} .

Then

\begin{matrix} (3) & \begin{matrix} (M_{f}^{P Q} + M_{fr}^{P Q}) {\ddot{p}}_{S}^{Q} + (C_{f}^{P Q} + C_{fr}^{P Q}) {\dot{p}}_{S}^{Q} + (K_{f}^{P Q} + K_{fr}^{P Q} + K_{fi}^{P Q}) p_{S}^{Q} = + [S_{fs}^{P M}] T_{S}^{M} + P_{f}^{P} - \\ ((M_{f}^{P Q} + M_{fr}^{P Q}) {\ddot{p}}_{I}^{Q} + (C_{f}^{P Q} + C_{fr}^{P Q}) {\dot{p}}_{I}^{Q} + (K_{f}^{P Q} + K_{fr}^{P Q} + K_{fi}^{P Q}) p_{I}^{Q}) + [S_{fs}^{P M}] T_{I}^{M}, \end{matrix} \end{matrix}

and

\begin{matrix} (4) & M^{N M} {\ddot{u}}^{M} + C_{(m)}^{N M} {\dot{u}}^{M} + I^{N} = - {[S_{fs}^{Q N}]}^{T} [p_{S}^{Q} + p_{I}^{Q}] + P^{N} . \end{matrix}

We see that the displacements in the solid are driven by the sum of the incident pressure, which forms an applied boundary traction, and the scattered pressure, the unknown field in the fluid.

The incident field is independent of the scattered field by convention. Therefore, it can be shown that it is a solution to the equation

(M_{f}^{P Q} + M_{fr}^{P Q}) {\ddot{p}}_{I}^{Q} + (C_{f}^{P Q} + C_{fr}^{P Q}) {\dot{p}}_{I}^{Q} + (K_{f}^{P Q} + K_{fr}^{P Q} + K_{fi}^{P Q}) p_{I}^{Q} = [S_{fs}^{P M}] T_{I}^{M},

so the fluid domain equation reduces to

\begin{matrix} (5) & \begin{matrix} (M_{f}^{P Q} + & M_{fr}^{P Q}) \ddot{p} ​_{S} ​^{Q} + (C_{f}^{P Q} + C_{fr}^{P Q}) \dot{p} ​_{S} ​^{Q} + (K_{f}^{P Q} + K_{fr}^{P Q} + K_{fi}^{P Q}) p ​_{S} ​^{Q} \\ ​ = [S_{fs}^{P M}] T_{S}^{M} + P_{f}^{P} . \end{matrix} \end{matrix}

This equation, using the scattered pressure as the unknown field variable, is solved together with Equation 4 using explicit or implicit integration or using the steady-state formulation (Coupled acoustic-structural medium analysis).

The scattered fluid traction, $T_{S}$ , depends on the incident pressure through the decomposition above and the solid motion at the boundary. In the absence of an impedance condition at this boundary, this results in the relation

\begin{matrix} (6) & \begin{matrix} T_{S} & = \frac{- 1}{ρ_{f}} \frac{\partial p_{S}}{\partial x} \cdot n^{-} = (\ddot{u} + \frac{γ}{ρ_{f}} \dot{u}) \cdot n^{-} - T_{I} \\ ​ = (\ddot{u} + \frac{γ}{ρ_{f}} \dot{u}) \cdot n^{-} - ({\ddot{u}}_{I} + \frac{γ}{ρ_{f}} {\dot{u}}_{I}) \cdot n^{-} \\ ​ = (\ddot{u} + \frac{γ}{ρ_{f}} \dot{u}) \cdot n^{-} + \frac{1}{ρ_{f}} \frac{\partial p_{I}}{\partial x} \cdot n^{-} . \end{matrix} \end{matrix}

If an impedance condition is applied at the fluid-solid boundary, application of the total pressure decomposition and Equation 13 results in the relation

\begin{matrix} (7) & \begin{matrix} T_{S} & ​ = (\ddot{u} + \frac{γ}{ρ_{f}} \dot{u}) \cdot n^{-} + \frac{1}{ρ_{f}} \frac{\partial p_{I}}{\partial x} \cdot n^{-} \\ - \frac{γ}{ρ_{f}} \frac{1}{c_{1}} p_{S} - (\frac{γ}{ρ_{f}} \frac{1}{k_{1}} + \frac{1}{c_{1}}) {\dot{p}}_{S} - \frac{1}{k_{1}} {\ddot{p}}_{S} \\ - \frac{γ}{ρ_{f}} \frac{1}{c_{1}} p_{I} - (\frac{γ}{ρ_{f}} \frac{1}{k_{1}} + \frac{1}{c_{1}}) {\dot{p}}_{I} - \frac{1}{k_{1}} {\ddot{p}}_{I} . \end{matrix} \end{matrix}

Virtual mass formulation

The “virtual mass approximation” is a simplification of the general scattered-field form of the coupled system Equation 1. Physically, it corresponds to considering the fluid wave speed to be very large compared to the characteristic structural wave speeds of interest; that is, fluid disturbances propagate and distribute in the field at an extremely high rate compared to the disturbances in the structure. This can be modeled by imposing an incompressibility condition on the fluid. Formally, applying such a condition results in the suppression of the terms related to time derivatives, volumetric drag, and impedance in the fluid equation for the coupled system:

\begin{matrix} (8) & \begin{matrix} K_{f}^{P Q} p^{Q} & ​ = + [S_{fs}^{P M}] {\ddot{u}}^{M} + \frac{1}{ρ_{f}} [S_{fs}^{P M}] n^{-} \cdot \frac{\partial p_{I}^{M}}{\partial x} \\ M^{N M} {\ddot{u}}^{M} + C_{(m)}^{N M} {\dot{u}}^{M} + I^{N} & ​ = - {[S_{fs}^{Q N}]}^{T} p^{Q} - {[S_{fs}^{Q N}]}^{T} p_{I} + P^{N} . \end{matrix} \end{matrix}

The fluid will also be assumed to have homogeneous mass density, $ρ_{f}$ . The fluid pressures, $p^{Q}$ , can be eliminated from the structural equation by inverting the fluid “stiffness,” $K_{f}^{P Q}$ to yield

\begin{matrix} (9) & \begin{matrix} M^{N M} {\ddot{u}}^{M} + C_{(m)}^{N M} {\dot{u}}^{M} + I^{N} = ​ & - {[S_{fs}^{Q N}]}^{T} {(K_{f}^{P Q})}^{- 1} [S_{fs}^{P M}] {\ddot{u}}^{M} \\ - {[S_{fs}^{Q N}]}^{T} p_{I} - {[S_{fs}^{Q N}]}^{T} {(K_{f}^{P Q})}^{- 1} [S_{fs}^{P M}] \frac{1}{ρ_{f}} n^{-} \cdot \frac{\partial p_{I}^{M}}{\partial x} + P^{N} . \end{matrix} \end{matrix}

Defining

\begin{matrix} (10) & M_{VMA}^{N M} \equiv {[S_{fs}^{Q N}]}^{T} {(K_{f}^{P Q})}^{- 1} [S_{fs}^{P M}], \end{matrix}

the equation for the structure can be written

\begin{matrix} (11) & (M^{N M} + M_{VMA}^{N M}) {\ddot{u}}^{M} + C_{(m)}^{N M} {\dot{u}}^{M} + I^{N} = - {[S_{fs}^{Q N}]}^{T} p_{I} - \frac{1}{ρ_{f}} M_{VMA}^{N M} n^{-} \cdot \frac{\partial p_{I}^{M}}{\partial x} + P^{N} . \end{matrix}

This equation makes the term “virtual mass” clear: the effect of an incompressible fluid surrounding a structure is to add inertia to the structural motion, $M_{VMA}^{N M}$ . For some shapes—in particular, cylinders, disks, and spheres—the virtual mass is a priori known from classical solutions.

The matrix operator $S_{fs}^{Q N}$ is the discrete form of the integral over the wetted or loaded surface (see Coupled acoustic-structural medium analysis). If the structure is long and thin, it is often modeled with beam elements. Then it is appropriate to replace the surface integral that generates the ${[S_{fs}^{Q N}]}^{T} p_{I}$ term with a more convenient volume integral term. Since

\begin{matrix} (12) & {[S_{fs}]}^{T} p_{I} = \int_{S_{fs}} H^{N} p_{I} n^{-} d S, \end{matrix}

where $H^{N}$ is a suitable interpolation function, Green's theorem can be applied to yield

\begin{matrix} (13) & {[S_{fs}]}^{T} p_{I} = - \int_{V_{wetted}} H^{N} \frac{\partial p_{I}}{\partial x} d S, \end{matrix}

where $V_{wetted}$ is the volume enclosed by the wetted surface. If $H^{N}$ is chosen to be the same interpolation function used for the beam and if the loading is similarly interpolated, this term becomes

\begin{matrix} (14) & \begin{matrix} {[S_{fs}]}^{T} p_{I} & ​ = - \int_{V_{solid}} H^{N} H^{M} d S \frac{\partial p_{I}^{M}}{\partial x} \\ ​ = - \frac{1}{ρ_{s}} M^{N M} \frac{\partial p_{I}^{M}}{\partial x} . \end{matrix} \end{matrix}

Therefore, the virtual mass equations for a beam immersed in an incompressible fluid are

\begin{matrix} (15) & (M^{N M} + M_{VMA}^{N M}) {\ddot{u}}^{M} + C_{(m)}^{N M} {\dot{u}}^{M} + I^{N} = - (\frac{1}{ρ_{f}} M_{VMA}^{N M} + \frac{1}{ρ_{s}} M^{N M}) \frac{\partial p_{I}^{M}}{\partial x} + P^{N} . \end{matrix}

This equation states that an immersed beam is loaded by a term due to the pressure gradient in the fluid due to the incident wave and that it experiences an additional amount of inertia due to the presence of the fluid.

Scattering load amplitude

The incident field in the fluid is assumed known; moreover, it is assumed to be associated with a separable solution to the scalar wave equation of the form

p_{I} (x, τ) \equiv p_{t} (τ) p_{x} (x),

where $τ$ is the “retarded time,” which is a function of the true time, t, and the position, $x$ . In applications of interest here, the incident wave is produced by an excitation outside the computational domain. However, its temporal part, $p_{t} (t)$ , is specified through an amplitude time history at a single point inside the computational domain, referred to as the “standoff point,” $x_{0}$ . The spatial term, $p_{x} (x)$ , is restricted to forms modeling plane and spherical wavefronts only, emanating from a specified “source” point, $x_{s}$ .

An incident wave produces a time-varying pressure at a spatial point of interest $x_{j}$ of the form

\begin{matrix} (16) & p_{I} (x_{j}, t) \equiv p_{t} (t) p_{x} (x_{j}) . \end{matrix}

The spatial variation of the field at point $x_{j}$ is expressed by either

p_{x} (x_{j}) = {[\frac{∥ x_{s} - x_{0} ∥}{∥ x_{s} - x_{j} ∥}]}^{α}

for spherically symmetric waves or simply

p_{x} (x_{j}) = 1

for plane waves. The amplitude observed at a point $x_{j}$ is related to that at the standoff point through this spatial variation, delayed by the time required for the wave to travel from the standoff point to the point $x_{j}$ :

\begin{matrix} (17) & \begin{matrix} p_{I} (x_{j}, t) & ​ = p_{t} (t - \frac{R_{j} - R_{0}}{c_{0}}) p_{x} (x_{j}) \\ ​ \equiv p_{t} (τ_{j}) p_{x} (x_{j}), \end{matrix} \end{matrix}

where $c_{0}$ is the wave speed in the fluid. We have defined

R_{0} \equiv ∥ x_{s} - x_{0} ∥,

and

R_{j} \equiv ∥ x_{s} - x_{j} ∥

for spherical waves and

R_{j} \equiv \frac{| (x_{j} - x_{s}) \cdot (x_{0} - x_{s}) |}{∥ x_{s} - x_{0} ∥}

for plane waves. $τ_{j}$ is referred to as “retarded time” because it includes a shift corresponding to the time required for the wave to move from the standoff point to $x_{j}$ . Of course, the delay may be positive or negative, depending on the relative distances between the point of interest, $x_{j}$ , and the source and between the standoff point and the source.

The incident wave boundary tractions in the solid equation Equation 4 result from direct substitution of Equation 17. The fluid boundary tractions in Equation 4 and the tractions in Equation 15 require the evaluation of the gradient of the incident wave pressure, assuming a source point with negligible speed:

\frac{\partial p_{I} (x_{j}, t)}{\partial x} = [- p_{x} {\dot{p}}_{t} \frac{1}{c_{0}} + p_{t} \frac{\partial p_{x}}{\partial R_{j}}] \frac{\partial R_{j}}{\partial x} .

First consider the case of generalized spatial decay, spherically symmetric waves with an exponent varying with the distance:

p_{x} = {(\frac{R_{0}}{R_{j}})}^{α (R_{j}, R_{0})},

where the exponent is chosen to be

α (R_{j}, R_{0}) \equiv \frac{(A + 1) R_{j}}{C R_{0} + (B + 1) R_{j}} .

This function reverts to acoustic decay for spherical waves if the constants $A$ , $B$ , and $C$ are zero and approaches plane wave propagation as $A \to - 1$ ; however, it allows some degree of flexibility in adapting the decay to experimental data. The gradient of this term is

\frac{\partial p_{x}}{\partial R_{j}} = - {(\frac{R_{0}}{R_{j}})}^{α} (\frac{α}{R_{j}}) [1 - \log (\frac{R_{0}}{R_{j}}) \frac{C R_{0}}{C R_{0} + (B + 1) R_{j}}] .

The gradient of the time-varying incident pressure simplifies for spherical acoustic waves to

\frac{\partial p_{I} (x_{j}, t)}{\partial x} = [\frac{- 1}{c_{0}} {\dot{p}}_{t} - \frac{1}{R_{j}} p_{t}] \frac{R_{0}}{R_{j}^{2}} (x_{j} - x_{s})

and to

\frac{\partial p_{I} (x_{j}, t)}{\partial x} = [\frac{- 1}{c_{0}} {\dot{p}}_{t}] \frac{1}{R_{0}} (x_{0} - x_{s})

for plane waves. The direction of the wave travel is

\frac{1}{R_{j}} (x_{j} - x_{s})

for spherical waves and equal to the constant vector

\frac{1}{R_{0}} (x_{0} - x_{s})

for plane waves. The time derivatives of the incident pressure field can be written

{\dot{p}}_{I} (x_{j}, t) = {\dot{p}}_{t} p_{x}

and

{\ddot{p}}_{I} (x_{j}, t) = {\ddot{p}}_{t} p_{x},

again neglecting the small terms associated with source point motion.

Steady-state dynamics

Development of incident wave loading in steady-state dynamics is as for the transient case, with a few exceptions. First, no source point motion is included in the steady-state formulation. Second, only spherical and planar sources are supported. The Fourier transform is applied to the excitation: with the usual assumptions on $p_{t} (τ_{j})$ , this is formally equivalent to substitution of a complex exponential in the steady-state formulation:

p_{t} (τ_{j}) \mapsto \tilde{p_{t}} e^{i Ω τ_{j}},

where $\tilde{p_{t}}$ is the complex magnitude of the incident pressure at the standoff point, as a function of frequency. Then

p_{I} (x_{j}, t) \mapsto \tilde{p_{t}} e^{i Ω τ_{j}} e^{- i Ω [\frac{R_{j} - R_{0}}{c_{0}}]} p_{x} (x_{j});

and the gradient becomes

\frac{\partial p_{I} (x_{j}, t)}{\partial x} \mapsto [\frac{- i Ω}{c_{0}} - \frac{1}{R_{j}}] \tilde{p_{t}} e^{- i Ω [\frac{R_{j} - R_{0}}{c_{0}}]} \frac{R_{0}}{R_{j}^{2}} (x_{s} - x_{j})

for spherical waves and

\frac{\partial p_{I} (x_{j}, t)}{\partial x} \mapsto [\frac{- i Ω}{c_{0}}] \tilde{p_{t}} \frac{1}{R_{0}} (x_{0} - x_{s})

for plane waves. The time derivatives of incident pressure, needed for the computation of the fluid traction on impedance boundaries, are

{\dot{p}}_{t} (t) \mapsto i Ω \tilde{p_{t}} e^{i Ω τ_{j}} e^{- i Ω [\frac{R_{j} - R_{0}}{c_{0}}]} p_{x} (x_{j}),

and

{\ddot{p}}_{t} (t) \mapsto - Ω^{2} \tilde{p_{t}} e^{i Ω τ_{j}} e^{- i Ω [\frac{R_{j} - R_{0}}{c_{0}}]} p_{x} (x_{j}) .

The value of the load applied at a point $x_{j}$ depends on the specified value $\tilde{p_{t}}$ and the spatial variation but includes a phase shift relative to the standoff point, $e^{- i Ω [\frac{R_{j} - R_{0}}{c_{0}}]}$ , corresponding to the time shifts in transient analysis.

Effect of structural motion

Generally, the effects of the motion of the source point, $x_{s} (t)$ , are included in the defined amplitude $p_{t} (t)$ of the incident wave field (specified at the standoff point $x_{0}$ ). That is, it is assumed that the stated amplitude reflects the time history at a specific fixed standoff point, including the effects of the source point motion, if any. For example, in underwater shock loading the incident wave field is produced by a pulsating gas bubble, which migrates toward the free surface of the water. Therefore, the corresponding incident wave amplitude model for this effect includes a moving source point.

The structural model affected by the incident wave may be in motion, as well. For example, a ship can move with respect to an explosive charge source during the loading. Because the wave motion is usually extremely fast compared to the ship motion and because incident wave pulse durations are typically short, some approximations apply. First, it can be assumed that the motion of the standoff point, $x_{0}$ , can be described by its position at time zero and by the average velocity during the incident wave loading, $v_{0}$ . In addition, the magnitude of this velocity can be assumed small compared to the speed of wave propagation, so the wave equation itself is unaffected by the relative motion. Furthermore, it can be assumed that the specified amplitude history, $p_{t} (t)$ , corresponds to observations made at the position of the standoff point at time zero, $x_{0} (0)$ .

Under these assumptions the effect of structural motion during incident wave loading can be modeled in a reference frame fixed with the standoff point, so the source point appears to move with $v_{0}$ :

x_{s} (t) \mapsto x_{s} (0) - v_{0} t .

This motion also affects the perceived arrival time of the wave at the standoff point. Defining the time of arrival at $x_{0} (0)$ as t and the time of arrival at $x_{0} (t)$ as $t_{I}$ , we can express the difference between these times in terms of $v_{0}$ :

c_{0} | t_{I} - t | = ∥ x_{s} (t) - v_{0} t - x_{0} (0) ∥ - ∥ x_{s} (t) - x_{0} (0) ∥ .

This time shift is maximum when the movement $v_{0}$ is aligned with $x_{s} (t) - x_{0} (0)$ . Therefore,

\begin{matrix} c_{0} | t_{I} - t | & ​ \leq ∥ x_{s} (t) - \frac{t ∥ v_{0} ∥ (x_{s} (t) - x_{0} (0))}{∥ x_{s} (t) - x_{0} (0) ∥} - x_{0} (0) ∥ - ∥ x_{s} (t) - x_{0} (0) ∥ \\ ​ \leq | t | ∥ v_{0} ∥ . \end{matrix}

Since we have assumed that the standoff motion is slow compared to the wave speed, we neglect this small time shift.

Underwater explosion—bubble load amplitude

An underwater explosion can lead to the formation of a highly compressed bubble that propels the surrounding water radially outward, generating an outward-propagating shockwave. As the bubble expands, the pressure inside decreases until it is considerably below the ambient pressure. After reaching a maximum radius with a minimum pressure, the bubble begins to contract. The contraction proceeds until the bubble collapses to a minimum radius. Because of a large pressure generated inside the bubble during this stage, the bubble begins to expand again, generating a second outward-propagating wave. Once the bubble expands to a second maximum radius, it contracts again. The same expansion-contraction sequence can repeat many times. At the same time during this pulsation process, the bubble migrates upward under the force of buoyancy. As long as the bubble does not reach the free water surface, the expansion-contraction sequence continues, each time with reduced amplitude.

Geers-Hunter model

Geers and Hunter (2002) proposed a phenomenological model that treats an underwater explosion as a single bubble event comprised of a shockwave phase and an oscillation phase, with the first phase providing initial conditions to the second. This model was slightly refined in a later paper, Geers and Park (2005).

Based on this model, the volume acceleration during the shockwave phase is given by

\begin{matrix} (18) & \ddot{V} (t) = \frac{4 π a_{c}}{ρ_{f}} P_{c} [0.8251 \exp (- 1.338 t / T_{c}) + 0.1749 \exp (- 0.1805 t / T_{c})], \end{matrix}

in which $P_{c} = K {(m_{c}^{1 / 3} / a_{c})}^{1 + A}$ and $T_{c} = k m_{c}^{1 / 3} {(a_{c} / m_{c}^{1 / 3})}^{B}$ , where $m_{c}$ and $a_{c}$ are the mass and radius of the explosive charge, respectively, and K, k, A, and B are constants for the charge material; $ρ_{f}$ is the mass density of the fluid.

Integration with $\dot{V} (0) = 0$ and further integration with $V (0) = (4 / 3) π a_{c}^{3}$ then yields

\begin{matrix} \dot{V} (t) & ​ = \frac{4 π a_{c}}{ρ_{f}} P_{c} [1.5857 - 0.6167 \exp (- 1.338 t / T_{c}) - 0.9690 \exp (- 0.1805 t / T_{c})] T_{c}, \\ V (t) & = \frac{4}{3} π a_{c}^{3} + \frac{4 π a_{c}}{ρ_{f}} P_{c} [1.5857 t / T_{c} - 5.8293 + 0.4609 \exp (- 1.338 t / T_{c}) + \\ + 5.3684 \exp (- 0.1805 t / T_{c})] T ​_{c} ​^{2}, \end{matrix}

and radial displacement and velocity follow as

a = {(\frac{3}{4 π} V)}^{1 / 3} and \dot{a} = \frac{1}{4 π} \dot{V} / a^{2} .

These expressions are evaluated at $t_{I} = 7 T_{c}$ to determine the initial conditions for the subsequent bubble response calculations during the oscillation phase. This choice was validated since, for a single set of charge constants, the initial condition values for values of $t_{I}$ between $3 T_{c}$ and $7 T_{c}$ produce essentially the same response during the oscillation phase, as demonstrated by Geers and Hunter (2002).

The following are the equations of motion for the doubly asymptotic approximation model to describe the evolution of the bubble radius, a, and migration, u, during the oscillation phase:

\begin{matrix} (19) & \begin{matrix} \dot{a} = - \frac{ϕ_{l 0}}{a} - \frac{1}{c_{f}} ({\dot{ϕ}}_{l 0} - {\dot{a}}^{2} - \frac{1}{3} {\dot{u}}^{2} - \frac{2}{3} \dot{u} \frac{ϕ_{l 1}}{a}), \\ \dot{u} = - 2 \frac{ϕ_{l 1}}{a} - \frac{1}{c_{f}} ({\dot{ϕ}}_{l 1} - 2 \dot{a} \dot{u}), \\ {\dot{ϕ}}_{l 0} = \frac{1}{(1 + ζ)} [(\frac{1}{2} + \frac{1}{2} \frac{ρ_{g}}{ρ_{f}} + ζ) ({\dot{a}}^{2} + \frac{1}{3} {\dot{u}}^{2}) - \frac{ρ_{g}}{ρ_{f}} c_{g} \frac{ϕ_{l 0}}{a} + \frac{2}{3} (1 + ζ) \dot{u} \frac{ϕ_{l 1}}{a} - Z], \\ {\dot{ϕ}}_{l 1} = \frac{1}{(1 + ζ)} [(1 + \frac{ρ_{g}}{ρ_{f}} + 2 ζ) \dot{a} \dot{u} - (1 - \frac{ρ_{g}}{ρ_{f}}) g a - \frac{ρ_{g}}{ρ_{f}} c_{g} (2 \frac{ϕ_{l 1}}{a} + \frac{ϕ_{g 1}}{a}) + \frac{3}{8} C_{D} {| \dot{u} |}^{E_{D}}], \\ {\dot{ϕ}}_{g 1} = \frac{1}{(1 + ζ)} [(2 + \frac{c_{g}}{c_{f}} + ζ) \dot{a} \dot{u} + \frac{c_{g}}{c_{f}} (1 - \frac{ρ_{g}}{ρ_{f}}) g a - c_{g} (2 \frac{ϕ_{l 1}}{a} + \frac{ϕ_{g 1}}{a}) - \frac{c_{g}}{c_{f}} \frac{3}{8} C_{D} {| \dot{u} |}^{E_{D}}], \end{matrix} \end{matrix}

where

ζ = \frac{ρ_{g} c_{g}}{ρ_{f} c_{f}}, and Z = \frac{1}{ρ_{f}} (P_{g} - p_{I} + ρ_{f} g u) + \frac{1}{3} [{(\frac{ϕ_{l 1}}{a})}^{2} - \frac{ρ_{g}}{ρ_{f}} {(\frac{ϕ_{g 1}}{a})}^{2}];

P_{g} = K_{c} {(\frac{V_{c}}{V})}^{γ}, ρ_{g} = ρ_{c} (\frac{V_{c}}{V}), c_{g} = c_{c} {(\frac{V_{c}}{V})}^{\frac{1}{2} (γ - 1)}, and c_{c} = \sqrt{\frac{γ K_{c}}{ρ_{c}}};

and $ρ_{c}$ is the charge mass density, $c_{f}$ is the sound speed in fluid, $V = (4 / 3) π a^{3}$ is the current volume of the bubble, $K_{c}$ is the adiabatic charge constant, $V_{c}$ is the volume of charge, $γ$ is the ratio of specific heats for gas, g is the acceleration due to gravity, and $p_{I} = p_{a t m} + ρ_{f} g d_{I}$ (where $p_{a t m}$ as the atmospheric pressure and $d_{I}$ as the depth of the charge's center). $C_{D}$ is an empirical flow drag parameter, which impedes the bubble's migration, and $E_{D}$ is an exponent that can be tuned to match experimental migration data [Geers and Park (2005)].

Seven initial conditions are needed. The first two are $a (t_{I}) = a_{I}$ , $\dot{a} (t_{I}) = {\dot{a}}_{I}$ , the second two are $u (t_{I}) = 0$ , $\dot{u} (t_{I}) = 0$ , the fifth one is $ϕ_{l 1} (t_{I}) = 0$ , and the remaining two are determined as

ϕ_{l 0} (t_{I}) = - a_{I} {\dot{a}}_{I} [1 - \frac{1}{2} (1 - \frac{ρ_{g I}}{ρ_{f}}) \frac{{\dot{a}}_{I}}{c_{f}} + ζ_{I}] + \frac{a_{I}}{c_{f}} Z_{I}

and

ϕ_{g 1} (t_{I}) = - ζ_{I}^{- 1} (1 - \frac{ρ_{g I}}{ρ_{f}}) \frac{g}{c_{f}} a_{I}^{2}

with

Z_{I} = ρ_{f}^{- 1} (P_{g I} - p_{I}) - \frac{1}{3} {(1 - \frac{ρ_{g I}}{ρ_{f}})}^{2} \frac{ρ_{f}}{ρ_{g I}} {(\frac{g a_{I}}{c_{g I}})}^{2} .

Then Equation 19 can be solved by using any suitable method for nonlinear ordinary differential equations; in Abaqus, a fourth-order Runge-Kutta integrator with variable time steps is used.

The incident pressure induced during the bubble response can be expressed as

\begin{matrix} (20) & p_{I} (x_{j}, t) \equiv p_{t} (t) p_{x} (x_{j}), \end{matrix}

where

p_{x} (x_{j}) = \frac{1}{R_{j}}

and

\begin{matrix} (21) & p_{t} (t) = \frac{ρ_{f}}{4 π} {(\frac{a_{c}}{R_{j}})}^{A} \ddot{V} ({(a_{c} / R_{j})}^{B} t), \end{matrix}

with $\ddot{V}$ given by Equation 18 for the shockwave phase ( $t_{I} < 7 T_{c}$ ), and

p_{t} (t) = \frac{ρ_{f}}{4 π} \ddot{V} (t) = ρ_{f} (a^{2} \ddot{a} + 2 a {\dot{a}}^{2})

for the oscillation phase ( $t_{I} \geq 7 T_{c}$ ). Here, $R_{j} \equiv ∥ x_{s} - x_{j} ∥$ .

“Waveless” model

The model of Geers and Hunter (2002) can be simplified to ignore the energy losses due to waves in the fluid and the gases. This “waveless” form of the equations more closely reproduces the results of earlier research on this phenomenon.

In the waveless model the shockwave phase of the loading is unchanged from the previous case. The equations of motion for the evolution of the bubble radius and migration during the oscillation phase are simplified, however, using the assumptions

\frac{1}{c_{f}} \mapsto 0,

\frac{ρ_{g}}{ρ_{f}} \mapsto 0,

and

c_{g} \mapsto 0.

Under these assumptions, the waveless equations of motion for the bubble dynamics are

\begin{matrix} (22) & \begin{matrix} \dot{a} & ​ = - \frac{ϕ_{l 0}}{a}, \\ \dot{u} & ​ = - 2 \frac{ϕ_{l 1}}{a}, \\ {\dot{ϕ}}_{l 0} & ​ = \frac{1}{2} ({\dot{a}}^{2} + \frac{1}{3} {\dot{u}}^{2}) + \frac{2}{3} \dot{u} \frac{ϕ_{l 1}}{a} - Z, \\ {\dot{ϕ}}_{l 1} & ​ = \dot{a} \dot{u} - g a + \frac{3}{8} C_{D} {| \dot{u} |}^{E_{D}}, \\ {\dot{ϕ}}_{g 1} & ​ = 2 \dot{a} \dot{u}, \end{matrix} \end{matrix}

where

Z = \frac{1}{ρ_{f}} (P_{g} - p_{I} + ρ_{f} g u) + \frac{1}{3} {(\frac{ϕ_{l 1}}{a})}^{2};

P_{g} = K_{c} {(\frac{V_{c}}{V})}^{γ}, ρ_{g} = 0, c_{g} = 0.

The initial conditions are determined as for the previously discussed case. Since the equation for $ϕ_{g 1}$ is now decoupled from the others, it does not need to be solved.

Reflections outside the computational domain

Abaqus allows the user to specify planes outside the computational domain that reflect the incident wave. This reflected wave is superposed, with a suitable time delay or phase shift, onto the wave arriving at the standoff point via the direct path, forming the total incident wave field for that source. This functionality is implemented for spherical or planar waves and reflection planes at any orientation, which may have arbitrary impedance properties. Only a single reflection from each plane is considered. However, multiple reflections inside the finite element computational domain are considered automatically.

Spherical waves

Consider the schematic arrangement of source point, standoff point, and reflection plane A as shown in Figure 1. The reflection plane can be specified by a normal vector, $n_{A}$ , and some point on the plane, $x_{A}$ . The distance from the source to the plane is $d_{A} \equiv ∥ (x_{A} - x_{S}) \cdot n_{A} ∥$ , and the plane is assumed to have a specific characteristic impedance, $Z_{A}$ . This information, together with the behavior of the source, is sufficient to characterize the reflected field in terms of an image source, located at the point $x_{S A} \equiv x_{S} + 2 d_{A} n_{A}$ . The source will be derived in steady state and then extended to the transient case.

The actual source is given by

P_{S} \equiv \tilde{p_{t}} e^{- i Ω [\frac{∥ x_{S} - x ∥}{c_{0}}]} \frac{∥ x_{S} - x_{0} ∥}{∥ x_{S} - x ∥},

and the image source is given by

P_{S A} \equiv {\tilde{p}}_{A} e^{- i Ω [\frac{∥ x_{S A} - x ∥}{c_{0}}]} \frac{∥ x_{S} - x_{0} ∥}{∥ x_{S A} - x ∥},

where the normalization length $∥ x_{S} - x_{0} ∥$ is chosen for convenience. Generally, a nonzero imaginary part of the plane's impedance $Z_{A}$ will imply that ${\tilde{p}}_{A}$ will also have a nonzero imaginary part. The unknown ${\tilde{p}}_{A}$ is found by enforcing the impedance condition at the plane:

P_{S A} + P_{S} = i Ω Z_{A} ({\tilde{u}}_{f} \cdot n_{A}),

and the conservation of linear momentum:

- Ω^{2} \tilde{ρ} {\tilde{u}}_{f} = \frac{\partial}{\partial x} [P_{S A} + P_{S}],

where $\tilde{ρ}$ is the complex density of the medium, including volumetric drag effects. Combining these two equations and simplifying, we have

{\tilde{p}}_{A} = \tilde{p_{t}} [\frac{1 - Z_{A} (\frac{1}{c_{0} \tilde{ρ}} + \frac{1}{i Ω \tilde{ρ} d_{A}})}{- 1 - Z_{A} (\frac{1}{c_{0} \tilde{ρ}} + \frac{1}{i Ω \tilde{ρ} d_{A}})}] .

In Abaqus the impedance property is specified using admittance constants $1 / c_{1}$ and $1 / k_{1}$ ,

\frac{1}{Z_{A}} = \frac{1}{c_{1}} + \frac{i Ω}{k_{1}},

so it is convenient to use

{\tilde{p}}_{A} \equiv [\frac{Q_{+}}{Q_{-}}] \tilde{p_{t}} = \tilde{p_{t}} [\frac{\frac{1}{c_{1}} + \frac{i Ω}{k_{1}} - (\frac{1}{c_{0} \tilde{ρ}} + \frac{1}{i Ω \tilde{ρ} d_{A}})}{- \frac{1}{c_{1}} - \frac{i Ω}{k_{1}} - (\frac{1}{c_{0} \tilde{ρ}} + \frac{1}{i Ω \tilde{ρ} d_{A}})}] .

This describes the source strength for a spherical image source in steady state, but in the transient case the complex density needs to be eliminated using

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

The expression for the image source can be expanded to

[- \frac{γ}{c_{1}} - \frac{1}{d_{A}} - i Ω \frac{ρ_{f}}{c_{1}} - i Ω \frac{γ}{k_{1}} - i Ω \frac{1}{c_{0}} + Ω^{2} \frac{ρ_{f}}{k_{1}}] {\tilde{p}}_{A}

​ = [\frac{γ}{c_{1}} - \frac{1}{d_{A}} + i Ω \frac{ρ_{f}}{c_{1}} + i Ω \frac{γ}{k_{1}} - i Ω \frac{1}{c_{0}} - Ω^{2} \frac{ρ_{f}}{k_{1}}] \tilde{p_{t}} .

In Abaqus the effects of fluid volumetric drag, the spreading loss due to the distance to the reflection plane, and the complex part of the admittance $1 / k_{1}$ are ignored, so the amplitude of the reflected wave is related to the amplitude of the incident wave by

{\tilde{p}}_{A} \equiv \tilde{p_{t}} [\frac{\frac{1}{c_{1}} - \frac{1}{c_{0} ρ_{f}}}{- \frac{1}{c_{1}} - \frac{1}{c_{0} ρ_{f}}}] \equiv Q \tilde{p_{t}} .

Therefore, the reflected spherical load is similar to the direct load, with magnitude reduced by the reflection impedance effect and by the greater distance travelled. In the time domain inversion of the steady-state result, retaining a phase shift corresponding to the arrival time at the standoff point, yields

p_{A} (t, x_{j}) \equiv Q p_{t} (τ_{j}) \frac{∥ x_{S} - x_{0} ∥}{∥ x_{S A} - x ∥},

where

τ_{j} \equiv t - \frac{∥ x_{S A} - x_{j} ∥}{c_{0}} + \frac{∥ x_{S} - x_{0} ∥}{c_{0}} .

The direction of this wave varies at each field point $x_{j}$ :

n_{s a} \equiv \frac{x_{S A} - x_{j}}{∥ x_{S A} - x_{j} ∥} .

“Soft” and “hard” limits

If a reflecting plane is considered “soft,” the pressure is zero there. Then $Z_{A} \equiv 0$ and

p_{A} = - p_{t} .

A “hard” reflecting plane refers to a zero fluid particle motion condition, implying $1 / Z_{A} \equiv 0$ . Then

p_{A} = p_{t} .

Planar waves

The schematic for reflected plane waves is similar to that for spherical waves, except for the geometry of the reflection. In contrast to the spherical case, when planar waves reflect from a boundary (see Figure 2), the direction of the reflected wave is common for all points on the wavefront. The location of the perceived source of the reflected wave is different, however. To calculate the phase shift for the reflected wave incident at an arbitrary point $x_{j}$ correctly, the locations of these perceived sources need to be calculated.

The reflected plane wave load can be calculated on the basis of the reflected wave direction and the time delay associated with the additional path length of the reflected wave. Referring to Figure 2, we see that the reflected wave direction is given by

n_{v} = n_{0} - 2 (n_{0} \cdot n_{s a}) n_{s a} .

The direct path wave travels the distance

∥ {\hat{x}}_{S} - x_{S} ∥ + ∥ {\hat{x}}_{S} - x_{0} ∥ + | (x_{j} - x_{0}) \cdot n_{0} |,

and the reflected path travels the distance

∥ {\hat{x}}_{S} - x_{S} ∥ + ∥ {\hat{x}}_{j} - x_{v} ∥ .

The time delay of interest for point $x_{j}$ , $τ_{j}$ , is related to the difference of these distances via the acoustic wave speed:

c_{0} τ_{j} = | (x_{j} - x_{0}) \cdot n_{0} | + ∥ {\hat{x}}_{S} - x_{0} ∥ - ∥ {\hat{x}}_{j} - x_{v} ∥ .

The first term on the right is the same as the delay due to the direct path; the rest is the increment due to the additional path length of the reflection. Some geometric manipulations allow this expression to be written in terms of the locations of $x_{j}$ , $x_{S}$ , $x_{0}$ , and the location of the spherical-wave virtual source, $x_{s a}$ :

c_{0} τ_{j} = | (x_{j} - x_{0}) \cdot n_{0} | + 2 r_{v} {(\frac{d_{s a} \cdot d_{0}}{∥ d_{s a} ∥ r_{0}})}^{2},

where

r_{0} = ∥ x_{S} - x_{0} ∥,

and

r_{v} = \frac{∥ d_{s a} ∥ r_{0}}{| n_{0} \cdot n_{s a} |} (\frac{1}{2} - \frac{(x_{j} - {\hat{x}}_{S}) \cdot d_{s a}}{{∥ d_{s a} ∥}^{2}}) .

Reflected planar waves are subject to reduction in amplitude due to impedance conditions at the reflection plane in the same manner as spherical waves, but there is no spatial attenuation. Unlike spherical waves, planar waves may yield no reflections for some orientations of the specified reflection plane: if the direction of the wave is parallel to the reflection plane, or strikes at an obtuse angle, no reflections are generated.