- Email: [email protected]
Numerical solution of the wave equation describing acoustic scattering and propagation through complex dispersive moving media
Nonlinear Analysis 71 (2009) e849–e854
Contents lists available at ScienceDirect
Nonlinear Analysis journal homepage: www.elsevier.com/locate/na
Numerical solution of the wave equation describing acoustic scattering and propagation through complex dispersive moving media Guy V. Norton ∗ Naval Research Laboratory, SSC, MS 39529-5004, USA
article
info
PACS: 43.20 43.20.B 43.20.H 43.20.M 43.30 Keywords: Time domain Dispersion Finite difference
abstract The numerical solution of acoustic pulse propagation through dispersive moving media requires the inclusion of attenuation and its causal companion, phase velocity. For acoustic propagation in a linear medium, Szabo [T.L. Szabo, J. Acoust. Soc. Am. 96 (1994) 491] introduced the concept of a convolutional propagation operator that plays the role of a casual propagation factor in the time domain. This operator was originally proposed to replace the term responsible for including losses from thermal conduction and viscosity of the fluid in the Westervelt equation. The operator has been successfully used in the linear wave equation for quiescent media. Development of the linear wave equation for sound in dispersive fluids with inhomogeneous flow will be described, along with a one-dimensional example. The resulting modified linear wave equation is solved via the method of finite differences. Published by Elsevier Ltd
1. Introduction Traditionally underwater acoustic propagation assumes a quiescent environment. That is, since the acoustic signal in sea-water travels nominally at 1500 m/s, the travel time between the source and some arbitrary receiver is less than, or is assumed to be less than, the time it takes the environment to change appreciably enough such that the acoustic field experiences this change. There are however situations where this assumption does not hold. For example when a benthic nepheloid layer (BNL) is present. A BNL is formed when sediments lying on the ocean floor are resuspended due to water motion. The water could be in motion due to currents, tides or via the shoaling and breaking of internal waves on the continental shelf. The BNL is dynamic, characterized by changing vertical thickness, concentration and speed of resuspended material. This increase in concentration of material suspended in the water column creates a dispersive environment. When the environment/medium is dispersive, causality dictates that there must be an attenuation of the propagating field. That is, the sound attenuation is frequency dependent. In addition to this, frequency dependent attenuation is the causal phase velocity. When the acoustic source is time limited, i.e. an acoustic pulse, it is composed of many frequencies. Thus to accurately model numerically such a situation requires that the solution be able to take into account the appropriate attenuation and phase velocity of all frequencies making up the source signature. Additionally direct modeling in the time domain is highly desirable, since it allows for more direct and efficient numerical solutions, and causality is always fulfilled. Based on an idea set forth by Blackstock [1] in the realm of nonlinear acoustics, Szabo proposed a way to include attenuation and dispersion effects directly in the time domain for both nonlinear [2] and linear propagation [3–5] in linear media through the inclusion of the so-called causal convolutional propagation operator into the corresponding threedimensional wave equation. By deriving a time domain version [3] of the Kramer–Krönig relationships (K–K), [6] he arrived at a general form for the operator. Szabo’s operator was originally defined in the context of lossy media obeying a frequency
∗
Tel.: +1 228 688 4686; fax: +1 228 688 5049. E-mail address: [email protected].
0362-546X/$ – see front matter. Published by Elsevier Ltd doi:10.1016/j.na.2008.12.009
e850
G.V. Norton / Nonlinear Analysis 71 (2009) e849–e854
power law attenuation. Waters et al., [7] showed that Szabo’s operator could be used for a broader class of media (as originally postulated by Blackstock), provided the attenuation possess a Fourier transform in a distributional sense. Norton and Novarini [8] clarified the use of the operator and its relationship to the time domain propagation factor. The time domain propagation factor is the parameter that is directly connected with the dispersive properties (attenuation and dispersion) of the medium. The modified wave equation that contains the causal convolutional propagation operator is rewritten incorporating the time domain propagation factor explicitly. The capability of the time domain propagation factor to correctly incorporate the dispersive traits of the medium into the linear wave equation was verified by solving the one-dimensional scalar, inhomogeneous wave equation in a dispersive medium via a finite-difference-time-domain (FDTD) scheme [8]. It was shown that for propagation in an isotropic dispersive (homogeneous) medium, attenuation and dispersion could correctly be taken into account while staying in the time domain. Therefore, if attenuation versus frequency is known (either from measurements or theoretically) and posses a generalized Fourier transform, then the time domain propagation factor (and so the corresponding causal convolutional propagation operator) can be obtained. And hence causal propagation can be modeled directly in the time domain. Norton and Novarini [9] have extended the technique to non-isotropic media (i.e. media in which the dispersive effects vary spatially). Excellent agreement between the FDTD with dispersion and the analytic result were observed for both the forward and backscattered field. Enhancements to the numerical code included making all difference approximations of the derivatives (both space and time) fourth order accurate as well as making the absorbing boundary conditions (utilizing the Complementary Operator Method) fourth order accurate [10]. Norton and Novarini [11] have utilized this code to investigate the scattering from and propagation though bubble clouds in the ocean. The results show that the presence of the bubble clouds introduces a great deal of ringing that stays near the surface. And that there can be a scattering effect from the bubble cloud before the scattering event at the air/water interface, resulting in returns from the bubble cloud arriving at a receiver prior to the interface scattering return. Norton and Novarini [12] have also applied this code to scattering and propagation through biological tissues, showing the effect that dispersion has on the backscattered signal. The target composition varied from brain, heart, kidney, liver and tendon. And finally Norton [13,14] extended the technique to include heterogeneous dispersive media, recasting the modified wave equation so that it had a spatially varying bulk modulus and density. The paper is divided into the following sections. First, the linear wave equation for sound in fluids with inhomogeneous flow will be developed. Next, the introduction of dispersive effects into the previous developed wave equation will be accomplished through the development and inclusion of the causal time domain propagation factor. Followed by a onedimensional example, showing proof of concept. The resulting modified linear wave equation is solved via the method of finite differences. 2. Wave equation for sound in fluids with unsteady inhomogeneous flow The following development follows from Refs. [15,16]. The full nonlinear fluid-dynamic equations for compressible fluid of uniform composition in the absence of dissipation can be written as Dv Dt Dρ
+
1
ρ
∇ p − g = 0,
+ ρ∇ • v = 0,
Dt Ds
(1)
= 0, Dt p = p (ρ, s) . The field variables are expressed as sums of ambient quantities (Subscript 0) and acoustic perturbations (primed quantities) so that the linearized acoustic equations can be formed. That is p = p0 + p0 (x, t )
ρ = ρ0 + ρ 0 (x, t )
(2)
v = v0 + v0 (x, t ) s = s0 + s0 (x, t ) . Thus one obtains,
ρ0 Dt v + v • ∇ v0 + ∇p − ∇ p0 = 0, ρ0 ρ02 Dt ρ 0 + v0 • ∇ρ0 + ρ 0 ∇ • v0 + ρ0 ∇ • v0 = 0, 0
0
1
Dt s0 + v0 • ∇ s0 = 0, p0 = c 2 ρ 0 +
∂p ∂s
s0 . 0
0
(3)
G.V. Norton / Nonlinear Analysis 71 (2009) e849–e854
e851
Where, Dt = ∂∂t + v0 • ∇ is the time derivative following the ambient flow. Simplification of the above linearized equations results when neglecting terms of 2nd order and higher in 1/L and 1/T , where L is the length scale over which ambient fields quantities have appreciable spatial variation and T is the corresponding time scale. This along with the introduction of a potential resulting in p = ADt Φ
(4)
where A is a function of ρ (density) and/or c (ambient sound speed), gives the following,
Dt
1 A
∇•
A2
A A ∇ Φ − Dt D Φ − D ∇ 2 Φ = 0. t t 2 ρ ρc ρ
(5)
The acoustic part of the flow velocity perturbation is given by
A ∇ Φ + O L−1 + O T −1 . ρ Setting A = −ρ yields for the acoustic pressure and the acoustic perturbation to the flow velocity, 0
v =
(6)
p = −ρ Dt Φ , v = −∇ Φ + O L 0
(7) −1
+O T
−1
.
(8)
One finally obtains the following expression describing acoustic propagation in unsteady inhomogeneous flow, 1
ρ
∇ • (ρ∇ Φ ) − Dt
1 c2
Dt Φ
=0
(9)
where c is a reference velocity (usually the thermodynamic sound speed in the medium assumed lossless [17], or the phase velocity at a limiting frequency) and ρ is the ambient density. 3. The causal time domain propagation factor A brief review for developing the time domain propagation factor is now given (see Refs. [4,8] for details). Assuming that propagation occurs through a lossy linear medium, the propagation can be shown to begoverned by a modified Eq. (9). The modification consists of including a causal convolutional propagation operator Lγ (t ) which controls attenuation and dispersion. It plays the role of a generalized dissipative term in the time domain. 1
∇ • (ρ (x) ∇ Φ (x, t )) − Dt
1
Dt Φ (x, t )
+
1
Lγ (t ) ∗ Φ (x, t ) = δ (x − xs ) s (t )
(10)
ρ (x) c (x) (x) where δ (x − xs ) is the Dirac delta function, s(t ) is the source signature at location (xs ), and the ∗ represents a convolution. c2
In the framework of generalized functions, assuming that the pressure field is a distribution and recalling that is equivalent to a convolution with a derivative of the Dirac delta function (1) distributional(1)differentiation δ (t ) ∗ f (t ) = f (t ) , the operator can be defined as, Lγ (t ) ≡ Γ (t ) ∗ δ (1) (t ) .
(11)
The function Γ (t ) represents the causal time domain propagation factor that accounts for causal attenuation, that is, it also governs dispersion in the system to insure causality. Eq. (10) can now be written in terms of this function 1
ρ (x)
∇ • (ρ (x) ∇ Φ (x, t )) − Dt
1 c2
(x)
Dt Φ (x, t )
+
∂ [Γ (t ) ∗ Φ (x, t )] = δ (x − xs ) s (t ) c (x) ∂t 1
(12)
Note, that the time domain propagation factor is the parameter that needs to be evaluated and not Lγ (t ) . The time domain propagation factor is directly connected with the dispersive properties of the medium. Szabo [3] defined the causal time domain propagation factor Γ (τ ) as the Fourier transform of the dispersive component of the complex propagation factor and showed it is given by
Γ (τ ) = −2 1+ (τ ) FT −1 {α (ω)},
(13)
where 1+ (τ ) represents the step function defined as [18]
1+ =
0
1
2 1
τ <0 τ =0 τ > 0.
(14)
It should be noted that Szabo’s original use of the operator was limited to a power law attenuation form, excluding some values of the exponent. Waters et al. [7], using distributional analysis, showed that Eq. (13) is more general and can be applied to any attenuation form for which an associated causal phase velocity exists.
e852
G.V. Norton / Nonlinear Analysis 71 (2009) e849–e854
4. Numerical experiment In order to validate the use of the time domain propagation factor in a media where fluid flow is present, a onedimensional numerical experiment will be performed. The input signal is a doublet of the form s (t ) = te−µt
2
(15)
where µ is a constant governing the time interval between the negative and positive peaks of the doublet. The amplitude spectrum of the signal is given by the Fourier transform of s(t ),
π 1/2 ω ω2 exp − . 2µ3/2 4µ
A (ω) =
(16)
The peak frequency (in Hz) is given by
√ fpk =
2µ
2π
.
(17)
The constant is chosen equal to 3.15 × 106 and therefore, fpk = 400 Hz. The maximum frequency of the source was approximately 1500 Hz. Eq. (12) was solved using the Finite-Difference-Time-Domain method. Numerical solutions at two receiver locations were Fourier transformed into the frequency domain. Both the attenuation and phase velocity versus frequency were determine, utilizing the following expressions,
α (f ) = c (f ) =
8.68
P (x2 ,f ) P (x1 ,f )
d 2π fd ARG (P (x2 , f ) /P (x1 , f ))
(18) (19)
where d is the distance in meters separating the two receivers. The modeling is carried out in one dimension (1D). The numerical grid consist of 1500 points with ∆X = 0.1 m. The time increment ∆t equaled 1.7539 × 10−5 s. The source is located at X = 5 m. The first receiver is located at X = 100.0 m and the second is located at X = 140.0 m. The non-dispersive sound speed is 1370.0 m/s and the density of water is taken as 1000 kg/m3 . The flow velocity is +50 m/s and the attenuation was assumed to have a linear power law dependence,
α (ω) = α0 ωy = 2.35 × 10−5 Np/m/rad/s ω1 .
(20)
A maximum of 1000 points were used for the convolution. Fig. 1 depicts the time series (signal) at the second receiver for various conditions. The solid line is for the case when the medium is non-dispersive and the flow velocity is 50 m/s. The dashed line represents the case when the medium is dispersive and the flow velocity is 50 m/s. Note that this signal arrives later in time, is highly attenuated and is spread out over a longer time. These effects are the consequence of propagating in a dispersive media. The dotted line is the case when the media is dispersive but the flow velocity is 0 m/s or a quiescent dispersive medium. Note that this signal arrives later in time than the signal whose flow velocity is 50 m/s. The time delay is exactly that which is due to the absence of the water flow. The dispersive effects are similar to the case when the medium is dispersive and the flow velocity is 50 m/s. Fig. 2 compares the retrieved attenuation versus frequency profiles for two cases. The first case is when the medium is non-dispersive or lossless. The retrieved attenuation should be 0 Np/m across the frequency content of the source and as shown by the dotted line it is clearly so. The analytic attenuation profile is shown by the solid line. The retrieved attenuation profile from the case when the medium is dispersive and having a flow velocity of 50 m/s is shown by the dashed line. Note the excellent comparison demonstrating that Eq. (12) accurately introduces dispersion effects due to the medium. Fig. 3 compares the retrieved phase velocity versus frequency for two cases. The first case is for a quiescent dispersive medium. The retrieved phase velocity (dotted line) should match that which is causally connected to the analytic attenuation profile (solid line). The comparison is again excellent. The dashed line is the retrieved phase velocity for the case when the flow velocity is 50 m/s. The resulting phase velocity is approximately 50 m/s faster for this case than when the flow velocity is 0 m/s as would be expected using Eq. (19) to determine the phase velocity. Fig. 4 compares the signal spectrum versus frequency at the two receiver locations (dashed and dotted lines) to the nondispersive signal spectrum (solid line), showing the effect that the frequency dependent attenuation has on the signal. The non-dispersive signal spectrum has a peak at 400 Hz and at the first receiver the peak is at 292 Hz and the peak at the second receiver is 257 Hz. This is a consequence of the medium having an attenuation based on a linear power law. The higher frequencies suffer more attenuated than the lower frequencies.
G.V. Norton / Nonlinear Analysis 71 (2009) e849–e854
Fig. 1. Received signals at second receiver.
Fig. 2. Comparison of the analytic (solid line) and retrieved (dotted and dashed line) attenuation.
Fig. 3. Comparison of the analytic (solid line) and retrieved (dotted and dashed line) phase velocity.
e853
e854
G.V. Norton / Nonlinear Analysis 71 (2009) e849–e854
Fig. 4. Signal spectrum for the non-dispersive (solid line) and dispersive (dotted and dashed lines) media versus frequency.
5. Concluding Remarks Direct time domain modeling of a one-dimensional pulse propagating in a dispersive inhomogeneous moving fluid has been performed via numerical (FDTD) solution of a modified linear wave equation. It was shown that the inclusion of a causal time domain propagation factor successfully accounts for attenuation and dispersion in moving medium. The fact that a causal time domain propagation factor can be defined for wave propagation in any weakly dissipative media allows for inclusion of attenuation and dispersive effects directly in the time domain. Acknowledgements This work has been supported by the Office of Naval Research, (Program Element No. 61153N) and by a grant of computer time at the DoD High Performance Computing Shared Resource Center (NRL-DC). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18]
D.T. Blackstock, J. Acoust. Soc. Am. 77 (1985) 2050. L. Szabo, Time domain nonlinear wave equations for lossy media, in: S. Hobaek (Ed.), in: Proc. 13th ISNA, vol. 89, Bergen, Norway, 1993. T.L. Szabo, J. Acoust. Soc. Am. 96 (1994) 491. T.L. Szabo, J. Acoust. Soc. Am. 97 (1995) 14. T.L. Szabo, Junuru Wu, J. Acoust. Soc. Am. 107 (2000) 2437. R.D.L. Krönig, J. Opt. Soc. Am. 12 (1926) 547. K.R. Waters, M.S. Hughes, G.H. Brandenburger, J.G. Miller, J. Acoust. Soc. Am. 108 (2000) 2114. G.V. Norton, J.C. Novarini, J. Acoust. Soc. Am. 113 (2003) 3024. G.V. Norton, J.C. Novarini, J. Comp. Acoust. 12 (2004) 501. G.V. Norton, J.C. Novarini, Math. Comp. Sim. 69 (2005) 467. G.V. Norton, J.C. Novarini, Compt. Phys. Comm. 174 (2006) 961. G.V. Norton, J.C. Novarini, MCB 4 (2007) 75. G.V. Norton, Numer. Methods Partial Differential Equations 23 (2007) 1420. G.V. Norton, Math. Comput. Simul. (in press). A.D. Pierce, J. Acoust. Soc. Am. 87 (1990) 2292. O.A. Godin, J. Acoust. Soc. Am. 112 (2002) 1269. L.E. Kinsler, A.R. Frey, A.B. Coppens, J.B. Sanders, Fundamental of Acoustics, Wiley, New York, 1982 (Chapter 7). R.N. Bracewell, The Fourier Transform and its Applications, 2nd ed., McGraw-Hill, New York, 1986, 57–61.
Contents lists available at ScienceDirect
Nonlinear Analysis journal homepage: www.elsevier.com/locate/na
Numerical solution of the wave equation describing acoustic scattering and propagation through complex dispersive moving media Guy V. Norton ∗ Naval Research Laboratory, SSC, MS 39529-5004, USA
article
info
PACS: 43.20 43.20.B 43.20.H 43.20.M 43.30 Keywords: Time domain Dispersion Finite difference
abstract The numerical solution of acoustic pulse propagation through dispersive moving media requires the inclusion of attenuation and its causal companion, phase velocity. For acoustic propagation in a linear medium, Szabo [T.L. Szabo, J. Acoust. Soc. Am. 96 (1994) 491] introduced the concept of a convolutional propagation operator that plays the role of a casual propagation factor in the time domain. This operator was originally proposed to replace the term responsible for including losses from thermal conduction and viscosity of the fluid in the Westervelt equation. The operator has been successfully used in the linear wave equation for quiescent media. Development of the linear wave equation for sound in dispersive fluids with inhomogeneous flow will be described, along with a one-dimensional example. The resulting modified linear wave equation is solved via the method of finite differences. Published by Elsevier Ltd
1. Introduction Traditionally underwater acoustic propagation assumes a quiescent environment. That is, since the acoustic signal in sea-water travels nominally at 1500 m/s, the travel time between the source and some arbitrary receiver is less than, or is assumed to be less than, the time it takes the environment to change appreciably enough such that the acoustic field experiences this change. There are however situations where this assumption does not hold. For example when a benthic nepheloid layer (BNL) is present. A BNL is formed when sediments lying on the ocean floor are resuspended due to water motion. The water could be in motion due to currents, tides or via the shoaling and breaking of internal waves on the continental shelf. The BNL is dynamic, characterized by changing vertical thickness, concentration and speed of resuspended material. This increase in concentration of material suspended in the water column creates a dispersive environment. When the environment/medium is dispersive, causality dictates that there must be an attenuation of the propagating field. That is, the sound attenuation is frequency dependent. In addition to this, frequency dependent attenuation is the causal phase velocity. When the acoustic source is time limited, i.e. an acoustic pulse, it is composed of many frequencies. Thus to accurately model numerically such a situation requires that the solution be able to take into account the appropriate attenuation and phase velocity of all frequencies making up the source signature. Additionally direct modeling in the time domain is highly desirable, since it allows for more direct and efficient numerical solutions, and causality is always fulfilled. Based on an idea set forth by Blackstock [1] in the realm of nonlinear acoustics, Szabo proposed a way to include attenuation and dispersion effects directly in the time domain for both nonlinear [2] and linear propagation [3–5] in linear media through the inclusion of the so-called causal convolutional propagation operator into the corresponding threedimensional wave equation. By deriving a time domain version [3] of the Kramer–Krönig relationships (K–K), [6] he arrived at a general form for the operator. Szabo’s operator was originally defined in the context of lossy media obeying a frequency
∗
Tel.: +1 228 688 4686; fax: +1 228 688 5049. E-mail address: [email protected].
0362-546X/$ – see front matter. Published by Elsevier Ltd doi:10.1016/j.na.2008.12.009
e850
G.V. Norton / Nonlinear Analysis 71 (2009) e849–e854
power law attenuation. Waters et al., [7] showed that Szabo’s operator could be used for a broader class of media (as originally postulated by Blackstock), provided the attenuation possess a Fourier transform in a distributional sense. Norton and Novarini [8] clarified the use of the operator and its relationship to the time domain propagation factor. The time domain propagation factor is the parameter that is directly connected with the dispersive properties (attenuation and dispersion) of the medium. The modified wave equation that contains the causal convolutional propagation operator is rewritten incorporating the time domain propagation factor explicitly. The capability of the time domain propagation factor to correctly incorporate the dispersive traits of the medium into the linear wave equation was verified by solving the one-dimensional scalar, inhomogeneous wave equation in a dispersive medium via a finite-difference-time-domain (FDTD) scheme [8]. It was shown that for propagation in an isotropic dispersive (homogeneous) medium, attenuation and dispersion could correctly be taken into account while staying in the time domain. Therefore, if attenuation versus frequency is known (either from measurements or theoretically) and posses a generalized Fourier transform, then the time domain propagation factor (and so the corresponding causal convolutional propagation operator) can be obtained. And hence causal propagation can be modeled directly in the time domain. Norton and Novarini [9] have extended the technique to non-isotropic media (i.e. media in which the dispersive effects vary spatially). Excellent agreement between the FDTD with dispersion and the analytic result were observed for both the forward and backscattered field. Enhancements to the numerical code included making all difference approximations of the derivatives (both space and time) fourth order accurate as well as making the absorbing boundary conditions (utilizing the Complementary Operator Method) fourth order accurate [10]. Norton and Novarini [11] have utilized this code to investigate the scattering from and propagation though bubble clouds in the ocean. The results show that the presence of the bubble clouds introduces a great deal of ringing that stays near the surface. And that there can be a scattering effect from the bubble cloud before the scattering event at the air/water interface, resulting in returns from the bubble cloud arriving at a receiver prior to the interface scattering return. Norton and Novarini [12] have also applied this code to scattering and propagation through biological tissues, showing the effect that dispersion has on the backscattered signal. The target composition varied from brain, heart, kidney, liver and tendon. And finally Norton [13,14] extended the technique to include heterogeneous dispersive media, recasting the modified wave equation so that it had a spatially varying bulk modulus and density. The paper is divided into the following sections. First, the linear wave equation for sound in fluids with inhomogeneous flow will be developed. Next, the introduction of dispersive effects into the previous developed wave equation will be accomplished through the development and inclusion of the causal time domain propagation factor. Followed by a onedimensional example, showing proof of concept. The resulting modified linear wave equation is solved via the method of finite differences. 2. Wave equation for sound in fluids with unsteady inhomogeneous flow The following development follows from Refs. [15,16]. The full nonlinear fluid-dynamic equations for compressible fluid of uniform composition in the absence of dissipation can be written as Dv Dt Dρ
+
1
ρ
∇ p − g = 0,
+ ρ∇ • v = 0,
Dt Ds
(1)
= 0, Dt p = p (ρ, s) . The field variables are expressed as sums of ambient quantities (Subscript 0) and acoustic perturbations (primed quantities) so that the linearized acoustic equations can be formed. That is p = p0 + p0 (x, t )
ρ = ρ0 + ρ 0 (x, t )
(2)
v = v0 + v0 (x, t ) s = s0 + s0 (x, t ) . Thus one obtains,
ρ0 Dt v + v • ∇ v0 + ∇p − ∇ p0 = 0, ρ0 ρ02 Dt ρ 0 + v0 • ∇ρ0 + ρ 0 ∇ • v0 + ρ0 ∇ • v0 = 0, 0
0
1
Dt s0 + v0 • ∇ s0 = 0, p0 = c 2 ρ 0 +
∂p ∂s
s0 . 0
0
(3)
G.V. Norton / Nonlinear Analysis 71 (2009) e849–e854
e851
Where, Dt = ∂∂t + v0 • ∇ is the time derivative following the ambient flow. Simplification of the above linearized equations results when neglecting terms of 2nd order and higher in 1/L and 1/T , where L is the length scale over which ambient fields quantities have appreciable spatial variation and T is the corresponding time scale. This along with the introduction of a potential resulting in p = ADt Φ
(4)
where A is a function of ρ (density) and/or c (ambient sound speed), gives the following,
Dt
1 A
∇•
A2
A A ∇ Φ − Dt D Φ − D ∇ 2 Φ = 0. t t 2 ρ ρc ρ
(5)
The acoustic part of the flow velocity perturbation is given by
A ∇ Φ + O L−1 + O T −1 . ρ Setting A = −ρ yields for the acoustic pressure and the acoustic perturbation to the flow velocity, 0
v =
(6)
p = −ρ Dt Φ , v = −∇ Φ + O L 0
(7) −1
+O T
−1
.
(8)
One finally obtains the following expression describing acoustic propagation in unsteady inhomogeneous flow, 1
ρ
∇ • (ρ∇ Φ ) − Dt
1 c2
Dt Φ
=0
(9)
where c is a reference velocity (usually the thermodynamic sound speed in the medium assumed lossless [17], or the phase velocity at a limiting frequency) and ρ is the ambient density. 3. The causal time domain propagation factor A brief review for developing the time domain propagation factor is now given (see Refs. [4,8] for details). Assuming that propagation occurs through a lossy linear medium, the propagation can be shown to begoverned by a modified Eq. (9). The modification consists of including a causal convolutional propagation operator Lγ (t ) which controls attenuation and dispersion. It plays the role of a generalized dissipative term in the time domain. 1
∇ • (ρ (x) ∇ Φ (x, t )) − Dt
1
Dt Φ (x, t )
+
1
Lγ (t ) ∗ Φ (x, t ) = δ (x − xs ) s (t )
(10)
ρ (x) c (x) (x) where δ (x − xs ) is the Dirac delta function, s(t ) is the source signature at location (xs ), and the ∗ represents a convolution. c2
In the framework of generalized functions, assuming that the pressure field is a distribution and recalling that is equivalent to a convolution with a derivative of the Dirac delta function (1) distributional(1)differentiation δ (t ) ∗ f (t ) = f (t ) , the operator can be defined as, Lγ (t ) ≡ Γ (t ) ∗ δ (1) (t ) .
(11)
The function Γ (t ) represents the causal time domain propagation factor that accounts for causal attenuation, that is, it also governs dispersion in the system to insure causality. Eq. (10) can now be written in terms of this function 1
ρ (x)
∇ • (ρ (x) ∇ Φ (x, t )) − Dt
1 c2
(x)
Dt Φ (x, t )
+
∂ [Γ (t ) ∗ Φ (x, t )] = δ (x − xs ) s (t ) c (x) ∂t 1
(12)
Note, that the time domain propagation factor is the parameter that needs to be evaluated and not Lγ (t ) . The time domain propagation factor is directly connected with the dispersive properties of the medium. Szabo [3] defined the causal time domain propagation factor Γ (τ ) as the Fourier transform of the dispersive component of the complex propagation factor and showed it is given by
Γ (τ ) = −2 1+ (τ ) FT −1 {α (ω)},
(13)
where 1+ (τ ) represents the step function defined as [18]
1+ =
0
1
2 1
τ <0 τ =0 τ > 0.
(14)
It should be noted that Szabo’s original use of the operator was limited to a power law attenuation form, excluding some values of the exponent. Waters et al. [7], using distributional analysis, showed that Eq. (13) is more general and can be applied to any attenuation form for which an associated causal phase velocity exists.
e852
G.V. Norton / Nonlinear Analysis 71 (2009) e849–e854
4. Numerical experiment In order to validate the use of the time domain propagation factor in a media where fluid flow is present, a onedimensional numerical experiment will be performed. The input signal is a doublet of the form s (t ) = te−µt
2
(15)
where µ is a constant governing the time interval between the negative and positive peaks of the doublet. The amplitude spectrum of the signal is given by the Fourier transform of s(t ),
π 1/2 ω ω2 exp − . 2µ3/2 4µ
A (ω) =
(16)
The peak frequency (in Hz) is given by
√ fpk =
2µ
2π
.
(17)
The constant is chosen equal to 3.15 × 106 and therefore, fpk = 400 Hz. The maximum frequency of the source was approximately 1500 Hz. Eq. (12) was solved using the Finite-Difference-Time-Domain method. Numerical solutions at two receiver locations were Fourier transformed into the frequency domain. Both the attenuation and phase velocity versus frequency were determine, utilizing the following expressions,
α (f ) = c (f ) =
8.68
P (x2 ,f ) P (x1 ,f )
d 2π fd ARG (P (x2 , f ) /P (x1 , f ))
(18) (19)
where d is the distance in meters separating the two receivers. The modeling is carried out in one dimension (1D). The numerical grid consist of 1500 points with ∆X = 0.1 m. The time increment ∆t equaled 1.7539 × 10−5 s. The source is located at X = 5 m. The first receiver is located at X = 100.0 m and the second is located at X = 140.0 m. The non-dispersive sound speed is 1370.0 m/s and the density of water is taken as 1000 kg/m3 . The flow velocity is +50 m/s and the attenuation was assumed to have a linear power law dependence,
α (ω) = α0 ωy = 2.35 × 10−5 Np/m/rad/s ω1 .
(20)
A maximum of 1000 points were used for the convolution. Fig. 1 depicts the time series (signal) at the second receiver for various conditions. The solid line is for the case when the medium is non-dispersive and the flow velocity is 50 m/s. The dashed line represents the case when the medium is dispersive and the flow velocity is 50 m/s. Note that this signal arrives later in time, is highly attenuated and is spread out over a longer time. These effects are the consequence of propagating in a dispersive media. The dotted line is the case when the media is dispersive but the flow velocity is 0 m/s or a quiescent dispersive medium. Note that this signal arrives later in time than the signal whose flow velocity is 50 m/s. The time delay is exactly that which is due to the absence of the water flow. The dispersive effects are similar to the case when the medium is dispersive and the flow velocity is 50 m/s. Fig. 2 compares the retrieved attenuation versus frequency profiles for two cases. The first case is when the medium is non-dispersive or lossless. The retrieved attenuation should be 0 Np/m across the frequency content of the source and as shown by the dotted line it is clearly so. The analytic attenuation profile is shown by the solid line. The retrieved attenuation profile from the case when the medium is dispersive and having a flow velocity of 50 m/s is shown by the dashed line. Note the excellent comparison demonstrating that Eq. (12) accurately introduces dispersion effects due to the medium. Fig. 3 compares the retrieved phase velocity versus frequency for two cases. The first case is for a quiescent dispersive medium. The retrieved phase velocity (dotted line) should match that which is causally connected to the analytic attenuation profile (solid line). The comparison is again excellent. The dashed line is the retrieved phase velocity for the case when the flow velocity is 50 m/s. The resulting phase velocity is approximately 50 m/s faster for this case than when the flow velocity is 0 m/s as would be expected using Eq. (19) to determine the phase velocity. Fig. 4 compares the signal spectrum versus frequency at the two receiver locations (dashed and dotted lines) to the nondispersive signal spectrum (solid line), showing the effect that the frequency dependent attenuation has on the signal. The non-dispersive signal spectrum has a peak at 400 Hz and at the first receiver the peak is at 292 Hz and the peak at the second receiver is 257 Hz. This is a consequence of the medium having an attenuation based on a linear power law. The higher frequencies suffer more attenuated than the lower frequencies.
G.V. Norton / Nonlinear Analysis 71 (2009) e849–e854
Fig. 1. Received signals at second receiver.
Fig. 2. Comparison of the analytic (solid line) and retrieved (dotted and dashed line) attenuation.
Fig. 3. Comparison of the analytic (solid line) and retrieved (dotted and dashed line) phase velocity.
e853
e854
G.V. Norton / Nonlinear Analysis 71 (2009) e849–e854
Fig. 4. Signal spectrum for the non-dispersive (solid line) and dispersive (dotted and dashed lines) media versus frequency.
5. Concluding Remarks Direct time domain modeling of a one-dimensional pulse propagating in a dispersive inhomogeneous moving fluid has been performed via numerical (FDTD) solution of a modified linear wave equation. It was shown that the inclusion of a causal time domain propagation factor successfully accounts for attenuation and dispersion in moving medium. The fact that a causal time domain propagation factor can be defined for wave propagation in any weakly dissipative media allows for inclusion of attenuation and dispersive effects directly in the time domain. Acknowledgements This work has been supported by the Office of Naval Research, (Program Element No. 61153N) and by a grant of computer time at the DoD High Performance Computing Shared Resource Center (NRL-DC). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18]
D.T. Blackstock, J. Acoust. Soc. Am. 77 (1985) 2050. L. Szabo, Time domain nonlinear wave equations for lossy media, in: S. Hobaek (Ed.), in: Proc. 13th ISNA, vol. 89, Bergen, Norway, 1993. T.L. Szabo, J. Acoust. Soc. Am. 96 (1994) 491. T.L. Szabo, J. Acoust. Soc. Am. 97 (1995) 14. T.L. Szabo, Junuru Wu, J. Acoust. Soc. Am. 107 (2000) 2437. R.D.L. Krönig, J. Opt. Soc. Am. 12 (1926) 547. K.R. Waters, M.S. Hughes, G.H. Brandenburger, J.G. Miller, J. Acoust. Soc. Am. 108 (2000) 2114. G.V. Norton, J.C. Novarini, J. Acoust. Soc. Am. 113 (2003) 3024. G.V. Norton, J.C. Novarini, J. Comp. Acoust. 12 (2004) 501. G.V. Norton, J.C. Novarini, Math. Comp. Sim. 69 (2005) 467. G.V. Norton, J.C. Novarini, Compt. Phys. Comm. 174 (2006) 961. G.V. Norton, J.C. Novarini, MCB 4 (2007) 75. G.V. Norton, Numer. Methods Partial Differential Equations 23 (2007) 1420. G.V. Norton, Math. Comput. Simul. (in press). A.D. Pierce, J. Acoust. Soc. Am. 87 (1990) 2292. O.A. Godin, J. Acoust. Soc. Am. 112 (2002) 1269. L.E. Kinsler, A.R. Frey, A.B. Coppens, J.B. Sanders, Fundamental of Acoustics, Wiley, New York, 1982 (Chapter 7). R.N. Bracewell, The Fourier Transform and its Applications, 2nd ed., McGraw-Hill, New York, 1986, 57–61.