PRIORITY CLAIM

[0001] This application claims the benefit of co-pending U.S. Provisional Patent Application 61/732,508, filed 3 December, 2012 and entitled "Solution of the Lossy Nonlinear Tricomi Equation with Application to Sonic Boom Focusing", which is hereby incorporated herein by reference in its entirety.

TECHNICAL FIELD

[0002] The present invention generally relates to supersonic flight and more particularly relates to a method for determining a pressure time history of a pressure wave originating from a vehicle traveling at or above supersonic speeds when the pressure wave undergoes a focusing due to convergence and intersection with a neighboring pressure wave originating from the vehicle.

BACKGROUND

[0003] The acceleration and turning maneuvers performed by vehicles travelling at supersonic speeds generate caustic surfaces that cause amplification of the propagating waves. These caustic surfaces often reach the ground and can lead to much greater acoustic pressure amplitudes than the sonic booms generated from a vehicle cruising at a steady supersonic trajectory. An additional source of increased acoustic pressure amplitude can occur because of how the rays in the wavefronts converge as a result of the maneuver, causing the ray tube area to decrease compared to the ray tube area of an aircraft travelling in a steady supersonic trajectory. That is, even though the rays have not reached the focus, waves undergo "focusing" as they approach the caustic due to the larger Blohkintzev invariant scaling factor. This increase in acoustic pressure is amplified even further once propagated into the caustic. Thus, it is important to quantify the effects of sonic boom focusing at a caustic to understand the overall amplification caused by unsteady aircraft maneuvers at supersonic speeds.

[0004] Focusing at the ground is unavoidable for any supersonic vehicle that accelerates to a speed greater than the cutoff Mach number. During an acceleration maneuver, a caustic surface is created at the vehicle the instant its speed becomes supersonic. The Mach cone of the vehicle becomes steeper as it continues to accelerate. A continuum of rays produced by the vehicle is launched normal to the Mach cone. Thus, portions of the wavefront generated later continually overtake portions of the wavefront generated at an earlier instant causing neighboring rays to converge and intersect. The caustic surface is formed at locations where the continuous convergence of rays occurs as the wavefront propagates away from the accelerating supersonic vehicle.

[0005] FIG. 1. Illustrates two-dimensional wavefronts propagated to the ground in a standard atmosphere due to supersonic acceleration/steady Mach/deceleration at a constant flight altitude. A line 100 indicates the portion of the wavefront that has not propagated through a caustic. A dashed line 102 indicates the portion of the wavefront that has already passed through the caustic. The dashed line 104 indicates the caustic location during wavefront propagation to the ground.

[0006] FIG. 1 illustrates a two-dimensional snapshot of propagating wavefronts from numerical simulations generated by an accelerating supersonic vehicle traveling through the standard atmosphere. The plot shows a planar cut from the caustic and wavefront surfaces in the plane vertical to the ground. In this example, the aircraft had accelerated at a constant

Mach rate from Mach 0.9 to 1.3, held steady at Mach 1.3 for 20 seconds, and decelerated at a constant Mach rate to Mach 0.9 while at a constant altitude of 45,000 ft. The wavefront represented by line 100 indicates the continuum of rays that have not passed through a caustic. Dashed line 102 indicates the continuum of rays that have already passed through the caustic. When the rays converge as the wavefront represented by line 100 propagates to the ground, the caustic surface is continuously formed as indicated by line 104. The shape and curvature of the wavefronts are due to the Mach rate of the aircraft and the speed of sound gradient in the atmosphere below the tropopause. The formation of the caustic is directly related to the aircraft trajectory and the atmospheric conditions. The amplitude at the caustic involves coupling between the caustic geometry parameters, the atmospheric conditions and the vehicle characteristics such as weight and shape. According to geometric acoustics, however, the amplitude would be infinite at the caustic. The caustic is located at the edge of the wavefield, violating the smooth gradient assumption of geometrical acoustics and making diffraction effects important. It is also important to note that the caustic itself is a byproduct of the mathematical approach when using geometric acoustics in describing the focusing physics. In the physical realm, the amplitude along the wavefront remains finite and also continuously and smoothly transitions from the illuminated region and into to the shadow zone due to diffraction.

[0007] In the past, attempts have been made to model sonic boom focusing (i.e., to predict a time pressure history of a pressure wave as it undergoes focusing) using the conventional nonlinear Tricomi equation which was derived to predict the acoustic field in the vicinity of a caustic. The theory and numerical implementation has been developed by Guiraud, Gill and Seebass, Plotkin, Auger and Coulouvrat, Kandil, and Sescu. Guiraud developed his well-known scaling law for the focusing amplitude of a step shock. Plotkin used this method and coded it into PCBoom for focus boom prediction. Auger and Coulouvrat developed the pseduospectral method for solving the nonlinear Tricomi equation. Their algorithm obtains a solution by adding an unsteady term and iterating in pseudotime with a split-step method. Nonlinear effects are solved in the time domain using a shock-capturing scheme and the diffraction effects are solved in the frequency domain. This method was later modified by Marchiano et al. by formulating the nonlinear Tricomi equation in terms of the acoustic potential and using the shock fitting method of Hayes et al. for solving the nonlinear term in their splitting method. Using the shock fitting method resulted in a significant improvement in reducing the number of iterations required to achieve a solution since their method was not constrained by a Courant-Friedrichs-Lewy (CFL) condition. Kandil closely replicated the results of Marchiano and Coulouvrat using a pseudospectral method and time domain finite differencing methods for his numerical implementation. More recently, Sescu solved the nonlinear Tricomi equation in conservative form using a weighted essentially non-oscillatory (WENO) scheme and a Galerkin method with a nonlinear limiting scheme to control the numerically induced oscillations that occur at sharp pressure gradients.

[0008] While these past attempts to model sonic boom focusing have been adequate, none of these past attempts have taken into consideration actual atmospheric dissipative effects on the pressure wave as it focuses. Accordingly, it would be desirable to predict/calculate the finite amplitudes at the caustic and to model the effects of diffraction in the acoustic pressure field at the caustic in a manner that takes actual atmospheric dissipative effects into consideration. It would further be desirable to provide a computer code that permits the computation of focus boom predictions that take atmospheric dissipative effects into consideration. Furthermore, other desirable features and characteristics will become apparent from the subsequent detailed description and the appended claims, taken in conjunction with the accompanying drawings and the foregoing technical field and background.

BRIEF SUMMARY

[0009] Various methods for determining a pressure time history of a pressure wave originating from a vehicle traveling at or above supersonic speeds when the pressure wave undergoes a focusing due to convergence and intersection with a neighboring pressure wave originating from the vehicle are disclosed herein.

[0010] In a first non-limiting embodiment, the method includes, but is not limited to, solving a lossy nonlinear Tricomi equation with a processor. The method further includes, but is not limited to, reporting an output from the processor, the output containing a solution to the lossy nonlinear Tricomi equation. The lossy nonlinear Tricomi equation includes a variable relating to an actual atmospheric dissipative effect in a vicinity of the focusing. The solution includes, but is not limited to, a prediction of the pressure time history of the pressure wave as the pressure wave undergoes focusing. The prediction reflects the actual atmospheric dissipative effect.

[0011] In a second non-limiting embodiment, the method includes, but is not limited to receiving an input at a processor, the processor configured to solve a lossy nonlinear Tricomi equation. The method further includes, but is not limited to, solving the lossy nonlinear Tricomi equation with the processor. The method still further includes, but is not limited to, reporting an output from the processor, the output containing a solution to the lossy nonlinear Tricomi equation. The lossy nonlinear Tricomi equation includes a variable relating to an actual atmospheric dissipative effect in a vicinity of the focusing. The solution includes, but is not limited to, a prediction of the pressure time history of the pressure wave as the pressure wave undergoes focusing. The prediction reflects the actual atmospheric dissipative effect.

BRIEF DESCRIPTION OF THE DRAWINGS

[0012] The present invention will hereinafter be described in conjunction with the following drawing figures, wherein like numerals denote like elements, and [0013] FIG. 1 is a graph illustrating propagation of two-dimensional wavefronts;

[0014] FIG. 2 is a graph illustrating a continuum of rays in a vicinity of a caustic;

[0015] FIG. 3 is a graph illustrating convergence monitoring as a function of pseudotime for four convergence parameters;

[0016] FIG. 4 contains multiple graphs illustrating the results of the Lossy Nonlinear Tricomi Equation code;

[0017] FIG. 5 contains multiple graphs illustrating the results of lossy nonlinear propagation of a sinusoid through a homogeneous medium with two relaxation processes;

[0018] FIG. 6 is a graph illustrating a comparison of the front shock from N-wave focusing;

[0019] FIG. 7 is a graph illustrating a comparison of the rear shock from N-wave focusing;

[0020] FIG. 8 is a graph illustrating a comparison of Lossy Nonlinear Tricomi Equation code solutions with losses versus without losses at a dimensionless Z location of 0.5;

[0021] FIG. 9 is a graph illustrating a comparison of Lossy Nonlinear Tricomi Equation code solutions with losses versus without losses at a dimensionless Z location of 0.15;

[0022] FIG. 10 is a graph illustrating sources of the acoustic input and caustic geometry for the Lossy Nonlinear Tricomi Equation code;

[0023] FIG. 11 is a graph illustrating an incoming N-wave form to the Lossy Nonlinear Tricomi Equation code;

[0024] FIG. 12 is a graph illustrating an exemplary pressure field solution from the Lossy Nonlinear Tricomi Equation;

[0025] FIG. 13 is a graph illustrating an exemplary Lossy Nonlinear Tricomi Equation prediction at a dimensionless Z location of 0.97 compared to a microphone measurement from the SCAMP flight test; [0026] FIG. 14 is a graph illustrating an exemplary Lossy Nonlinear Tricomi Equation prediction at a dimensionless Z location of 0.67 compared to a microphone measurement from the SCAMP flight test;

[0027] FIG. 15 is a graph illustrating an exemplary Lossy Nonlinear Tricomi Equation prediction at a dimensionless Z location of 0.15 compared to a microphone measurement from the SCAMP flight test;

[0028] FIG. 16 is a graph illustrating an exemplary Lossy Nonlinear Tricomi Equation prediction at a dimensionless Z location of -0.22 compared to a microphone measurement from the SCAMP flight test; and

[0029] FIG. 17 is a flow diagram illustrating a non-limiting embodiment of a method for determining a pressure time history of a pressure wave as it focuses.

DETAILED DESCRIPTION

[0030] The following detailed description is merely exemplary in nature and is not intended to limit the invention or the application and uses of the invention. Furthermore, there is no intention to be bound by any theory presented in the preceding background or the following detailed description.

[0031] Various symbols will be used herein. A summary explanation of those symbols are provided here:

[0032] p = acoustic pressure, Pa

[0033] fac = characteristic acoustic frequency, Hz

[0034] l _ac = characteristic acoustic wavelength, m

[0035] p _ac = characteristic acoustic pressure, Pa

[0036] ε = dimensionless diffraction parameter

[0037] R _C au = radius of curvature of the line defined as the intersection of the caustic surface and the plane directed by both the tangent ray and the caustic normal, m

[0038] Ryay = radius of curvature of the ray tangent to the caustic, m [0039] R _tot = relative radius of curvature between the radius of curvature for the ray tangent to the caustic and the radius of curvature of the line defined as the intersection of the caustic surface and the plane directed by both the tangent ray and the caustic normal, m

[0040] d= diffraction boundary layer thickness, m

[0041] J = dimensionless distance from the caustic in the direction normal to the caustic

[0042] z = distance from the caustic in the direction normal to the caustic, m

[0043] x = distance from the caustic in the direction tangent to the caustic, m

[0044] t = dimensionless retarded time

[0045] t = time variable, s

[0046] = pseudotime variable

[0047] Co = ambient speed of sound, m/s

[0048] po = ambient density, kg/m ³

[0049] β= coefficient of nonlinearity

[0050] M _x = ratio of the wind speed, uox, in the direction tangent to the caustic to the ambient speed of sound

[0051] M _z = ratio of the wind speed, uoz, in the direction normal to the caustic to the ambient speed of sound

[0052] δ = diffusivity of sound, m ² /s

[0053] m _v = dispersion parameter for the v-th relaxation component

[0054] T _V = relaxation time for the v-th relaxation component, s

[0055] c _∞iV = frozen speed of sound for the v-th relaxation component, m/s

[0056] The caustic is formed by rays whose ray tube area goes to zero and infinitesimally adjacent rays converge. In three dimensions, the caustic is a surface and in two dimensions the caustic is a line.

[0057] FIG. 2 illustrates the continuum of rays that propagate concavely upward and successively generate a caustic line 106 that is concave downward. The z-axis is normal to the caustic at point O, the x-axis is tangential to the caustic at point O and the y-axis is into the plane of the figure. Above caustic line 106, in a region 108 referred to as an illuminated zone, there are two rays that pass through each location - an incoming ray and an outgoing ray. The incoming ray has not yet passed through the caustic and the outgoing ray has already passed through the caustic. As an acoustic disturbance propagates along each ray and approaches the caustic, amplification and phase changes of the pressure waveform occur. The phase change contributes to transforming the classical "N" shape of a conventional sonic boom into a waveform shaped like a "U". For the linear case, the "U" shape of the outgoing waveform is the result of a Hilbert transform of the incoming N-wave. Below caustic line 106 is region 110, referred to as a shadow zone. Geometric acoustics predicts that no rays pass into the shadow zone. In the physical world, however, sound does actually penetrate into the shadow zone due to diffraction. The result is an evanescent wave whose amplitude decays exponentially away from the caustic. Conventional ray theory is insufficient to predict the acoustic field in the diffraction region in the vicinity of the caustic in both regions 108 and 1 10, the illuminated zone and the shadow zone, respectively.

[0058] As shown in FIG. 2, a continuum of rays pass in the vicinity of caustic line 106. The rays refract upward due to a local speed of sound gradient. The illuminated zone (region 108) lies above the caustic. The shadow zone (region 110) lies below the caustic. The outgoing ray has already passed through the caustic. The incoming ray has not yet passed through the caustic.

[0059] The nonlinear Tricomi equation from the previous work of others has been augmented to include atmospheric loss mechanisms and an additional term to account for wind in the direction tangent to the caustic. The acoustic pressure perturbations and shocks in a sonic boom are typically modified due to the effects of absorption and dispersion as it propagates in the far-field from the vehicle in flight down to the caustic. Absorption due to the thermoviscous effects and molecular relaxation of oxygen and nitrogen in the atmosphere contributes to the finite rise time and thickness of shocks. Molecular relaxation effects cause dispersion, producing asymmetry in the shocks. The absorption and dispersion processes must be included in the incoming waveform to capture how they influence the acoustic pressure field near the caustic.

[0060] Below is presented an equation that may be used to determine the time pressure history of a pressure wave as it undergoes focusing: d ² p _ d ² p M _z d ² p ( 2 1/, M 2 Λ d ² p

≡ ^{2 z} dP ε dtdz 2 2

2

[0061] The summation over v is from 1 to 2 for oxygen and nitrogen relaxation, p = p/p _ac is the dimensionless pressure, M _ac = p _ac /p ₀ c% is the characteristic acoustic Mach number and p _ac is the characteristic acoustic pressure which is usually specified as the peak pressure of the incoming waveform corresponding to a distance of one diffraction boundary layer thickness. The first two terms in Equation 1 account for diffraction and correspond to the linear Tricomi equation, the third term accounts for the z-component of wind, the fourth term accounts for the x-component of wind, the fifth term accounts for nonlinearity and the sixth term accounts for thermoviscous absorption and molecular relaxation due to oxygen and nitrogen. This equation shall be referred to herein as the "lossy nonlinear Tricomi equation" or "LNTE".

[0062] Auger's derivation of his lossless nonlinear Tricomi equation provides an explanation to support the assumption that the diffraction in the y-direction is negligible, which is typically valid when the 3c - z plane coincides in the region of the zero-degree azimuth from underneath the supersonic vehicle. It also assumes the propagation through the caustic occurs in a preferred direction along the x -axis and, consequently, the 3c - direction terms have been put in terms of the dimensionless time variable. The dimensionless time variable that includes the eikonal function is defined as:

[0063] The eikonal function used here was obtained by analyzing the local caustic geometry in terms of curvilinear coordinates. There are similar derivations for determining a suitable eikonal function. The eikonal function used here was version derived by Auger in his PhD thesis. The diffraction parameter is defined as

[0064] ε is a ratio of the characteristic acoustic wavelength, X _ac , to the diffraction

1 boundary layer thickness, d. The total relative radius of curvature is— !— ¹ (see

^tot ^cau ^ray

FIG. 2). The dimensionless distance to the caustic is defined by

[0065] This is the distance away from the caustic in the direction normal to the caustic divided by the diffraction boundary layer thickness. The normalized component of wind in the x- and z-directions are defined as M _x = u _ox l c ₀ and M _z = u _oz / c ₀ , respectively. The dimensionless thermoviscous absorption and dispersion coefficients are given by a = δ f _ac l c and θ _ν = f _v m _v = f _ac r _v (c _v - cl ) lcl , respectively, and the dimensionless relaxation time is F _v = f _ac r _v [24]. The absorption and dispersion terms are newly added terms in the nonlinear Tricomi equation to model the acoustic field near a caustic. Comparison of the relative orders of magnitudes of the terms as presented in Equation 1 reveals that the loss terms are several orders of magnitude lower than the diffraction terms, nonlinear term and wind terms when ranked in magnitude according to the diffraction parameter, ε. However, the terms that include absorption and dispersion also contain a third-derivative with respect to time and, as a consequence, in the frequency domain would grow in proportion to frequency cubed. The frequency range of interest here (up to 10 kHz) is 10 to 20 times higher than the frequency range of interest from the work of previous authors. Thus, when high frequency signal content is present in certain parts of the waveform (such as the shocks in N-waves) the loss terms can reach within one to two orders of magnitude of the diffraction and nonlinear terms. A strong motivating factor for including the loss terms is it is well understood that high frequency content, manifested in the shocks, can have an important impact on perceptual metrics of sonic booms used for subjective evaluation. Modeling the shocks as accurately as possible is important to permit accurate assessments of how "loud" a focus boom will sound to human observers. [0066] The time-wise boundary conditions assume that the waveform is undisturbed for both large negative and large positive time, which is interpreted as:

[0067] The pressure field decays exponentially and goes to zero as z approaches negative infinity. For the domain that extends infinitely away from the caustic, the boundary condition is: p(t, z→-∞) = 0 (6)

[0068] With continuing reference to FIG. 2, as one moves away from the caustic along the positive z -axis the solution must match the incoming and outgoing rays from geometric acoustics. That is, the Tricomi formulation in the vicinity of the caustic must also match the pressure field consisting of an incoming ray and an outgoing ray when there is sufficient distance away from the caustic. The boundary condition for the positive z -axis becomes

^-, z _→∞ ) _¾ z- ^1/4 [^ ₊ f z ^3/2 ) _{+ G} ( -f z ^3/2 )J (7) where F is the incoming waveform and G is the outgoing waveform, both of which have been normalized by the characteristic acoustic pressure, p _ac . The complication with Equation 7 being that only the incoming wave, F, is known and G is unknown. In the literature, the boundary condition for large, positive z can be written as a radiation condition only in terms of F.

- _{1/4 + rl/4 + r5M} _ _{P = 2f} ,(- , . _3/2 |

dt 3Έ 4 ³ '

where the prime symbol applied to F refers to a derivative with respect to its argument.

[0069] The numerical solution to Equation (1) includes newly added terms to account for the loss mechanisms and the x-wind term. A pseudotime variable, σ , and unsteady term were added to Equation (1). The equation implemented in the numerical algorithm to solve the LNTE was dtda cE ² dt ² ε dtdz _ε ² ε ² dt ²

[0070] As the computations converge towards the numerical solution, the unsteady term on the left hand side of Equation (9) tends towards zero and the desired solution to the LNTE is achieved. In some embodiments, an approach to solving Equation (9) is a splitting method that solves for the diffraction terms, the nonlinear term and the absorption/dispersion terms in separate steps. The nonlinear term is solved in the time domain and the linear terms are solved in the frequency domain. The frequency domain steps are solved for each harmonic component. The time domain step is solved for each l coordinate. The discretized computational domain in both the z - and t -directions have uniform spacing.

[0071] The upper and lower boundary conditions must be modified accordingly to be consistent in both the frequency domain and time domain and also be suitable for a computational domain of finite extent. The boundary condition in Equation (8) becomes

[0072] The boundary condition in Equation (10) is applied at the maximum z location, z _max , in the computational domain, typically at z _max = 1 which corresponds to one diffraction boundary layer thickness. p _n is the w-th harmonic of the pressure field at a particular z - coordinate, which in the case of the above boundary condition is z _max . F _n is the w-th harmonic of the FFT of the incoming waveform, F, and ω„ is the w-th harmonic of the dimensionless angular frequency. The z derivative in Equation (10) is discretized to second order.

[0073] For the lower boundary condition, the solution to the lossless, windless, linear Tricomi equation includes an asymptotic expression for the Airy function. The spectral solution of the linear Tricomi equation is: p _n (z) =

where sgn is the signum function, Ai is the Airy function and ξ is the argument of the Airy function. Equations 11a and l ib were combined and manipulated to obtain the boundary condition for the lower part of the computational domain:

[0074] As in Equation (10), Equation (12) is applied for each harmonic component, n. The J derivative in Equation (12) is discretized to first order. The lower edge of the boundary condition typically extends to a minimum z location of z _min = -1 . The boundary condition in Equation (12) is different than what has been specified previously where pit , z _min ) = 0 for the minimum z boundary location. The disadvantage to using that boundary condition is the z _min = -1 pressure amplitudes are typically not close enough to zero for that approximation to be valid in a finite computational domain. Applying a zero pressure-amplitude at that location will influence the pressure amplitudes in the vicinity of the caustic leading to underestimated values for the pressure. The z _min value should be at least 2 in order to avoid underestimating the pressure amplitudes near the caustic, leading to a considerably larger computational domain when using a p(t, z _min ) = 0 boundary condition. Thus, as also mentioned by Auger in his PhD thesis, it is more computationally efficient and more accurate to use the linear Tricomi solution instead of a zero-pressure value for the boundary condition at z _min .

[0075] The first step in the computational process transforms the pressure field to the frequency domain via Fast Fourier Transform (FFT) and then solves for the diffraction and z-wind component terms.

ιω„— VPζn — 1/ dp

ζ- = ω _η 2ζρ _η + , 2

- ιω, (13)

ΪΓ2 ε dz [0076] Equation (13) is solved using a fully implicit numerical discretization scheme that is second order in space and first order discretization in pseudotime. The corresponding boundary conditions are Equations (10) and (12). The solution for Equation (13) is determined for each harmonic component, n, of the pressure field. Next, the x-wind, absorption and dispersion terms are solved in the second step.

[0077] Equation (14) is solved using an exact solution. This is a computational step not performed by previous authors of sonic boom focusing work because the terms are newly derived for this invention. In the third step, an inverse FFT is applied to the pressure field and the nonlinear step is solved using: dp _ fi M _ac dp ²

σ ε ² dt

[0078] Equation (15) is solved using an exact solution at each Έ location and is different from the methods used by other authors when solving the nonlinear Tricomi equation. The physical dimensionless time base is adjusted according to the Poisson solution and is linearly interpolated back to the uniform time base.

[0079] The pseudotime step size is chosen using the criteria by Cleveland for the splitting method in his thesis. The pseudotime increment value is applied equally for each of the above three computational steps describing Equations (13) to (15). The numerical shock "distance," which in this case the "distance" is the pseudotime, is calculated for each z location. The pseudotime step increment is chosen to be 1/10 ^th the minimum shock "distance" calculated from the entire solution domain. As a result, the pseudotime step increment can vary from one iteration to the next depending on the magnitude of pressure gradients in the solution. The approach here for determining the pseudotime step increment to account for the nonlinear solution step is conservative. This method, in conjunction with the step that accounts for absorption and dispersion, is an attempt to best preserve the resolution of any shocks that might develop as a result of focusing. [0080] The above three steps from Equations (13) to (15) are repeated iteratively until convergence is achieved. The starting pressure field is zero. The code iterates to determine the linear Tricomi solution, then uses the linear solution as the initial guess to iterate towards the solution of the lossy nonlinear Tricomi equation. The iterative convergence of the solution is monitored in several ways as the pressure field evolves in pseudotime. The code performs a global search over all of the z locations and determines the z -coordinate that shows the largest change in pressure value during the nonlinear solution step for all of the time-values at that z -coordinate. The convergence parameter for that z location is the maximum difference in the solution between the current and the next iterations divided by the pseudotime step size. The normalization by the pseudotime step size was chosen because the pseudotime increment can vary and the solution difference per unit pseudotime places the amount of change per iteration in the same relative value from one iterative step to the next. Additionally, the code also performs a global search over all of the ω frequencies and finds the frequency with the largest change in spectral amplitude during the diffraction solution step for all of the z locations at that frequency. The convergence parameter for that ω frequency is the maximum difference in the solution from the current iteration to the next iteration, again normalized by the pseudotime step increment, for the spectral amplitudes at that ω for all of the corresponding z locations. An additional convergence parameter corresponds to the top edge of the domain to monitor the evolution of the outgoing waveform. The convergence parameter for the z _max location is the maximum value of the difference of the solution at z _max between the current and previous iterations divided by the pseudotime step size. Lastly, the pseudotime step value itself is a convergence parameter that correlates to the largest slope in the solution domain. For the nonlinear solution step, the pseudotime step increment is determined based on a search for the maximum pressure gradient, with respect to dimensionless time, in the solution. Thus, the stabilization of the pseudotime step increment as the solution is iterated is also an indicator for convergence since it represents how the maximum slope in the solution domain evolves in pseudotime.

[0081] An example of convergence monitoring using these parameters is given in FIG. 3. In FIG. 3, convergence monitoring as a function of pseudotime for the four convergence parameters is illustrated. Observe that the upper boundary difference 1 12, the nonlinear steepening difference 1 14 and the diffraction spectral difference 1 16 all reduce by four orders of magnitude and level off after a pseudotime of 6. The pseudotime step increment 118 remains fairly constant after a pseudotime of 5.

[0082] The values for these convergence parameters reduce in value by four orders of magnitude after a pseudotime value of 7 over the course of approximately 23,000 iteration steps. The convergence parameter values shown in FIG. 3 correspond to the experimental validation case described below.

[0083] Two numerical validation tests were performed for the coded implementation of the LNTE. The first check examines the ability of the code to model the diffraction effects. The second check examines the code's ability to model the nonlinear, absorption and dispersion effects. The first check compared the computational output from the code to the analytical solution of the linear Tricomi equation by only enabling the diffraction effects in the computer code. The second check compared the computational output of the code to a spectral solution of lossy nonlinear propagation of a sinusoid in a homogenous medium. This check was accomplished numerically by disabling the diffraction effects in the code and retaining the nonlinear, absorption and dispersion effects.

[0084] The first numerical validation compares the analytical solution for the linear Tricomi equation to the LNTE computation with only the diffraction effects enabled in the code. The time-domain version of the analytical solution was obtained from the expression in:

where 'IFFT is the inverse Fourier transform, sgn is the signum function and Ai is the Airy function.

[0085] The comparison plots illustrated in FIG. 4 show the results of the diffraction numerical validation at four different j locations. FIG. 4 illustrates the results of the LNTE code with only the diffraction terms included in the computations (dotted line 120) compared to the analytical solution of the lossless linear Tricomi equation (line 122) at four different Έ locations: a) z = 1.0 , b) z = 0.5 , c) z = 0. and d) z = -0.5 .

[0086] An N-wave with finite rise time at the shocks was used as the incoming waveform. The incoming N-wave had a characteristic pressure of 50 Pa and a waveform duration of 0.3 seconds. The agreement is very good in the illuminated zone (the two upper graphs illustrated in FIG. 4), the shadow zone (the lower right hand graph illustrated in FIG. 4) and on the caustic line (the lower left hand graph illustrated in FIG. 4).

[0087] The second validation check examines the lossy nonlinear propagation of a sinusoid through a homogeneous medium. This situation can be adequately modeled by a lossy Burgers equation. A spectral solution for this appears in the nonlinear acoustics book edited by Hamilton and Blackstock. They solve the Burgers equation in the frequency domain for a finite number of harmonics using a Runge-Kutta scheme. This is the same numerical check Cleveland used to validate his code that implements the augmented Burgers equation; however, the Hamilton-Blackstock solution is stated for only one relaxation process. Applicant modified this to include two relaxation processes for numerical validation of his augmented KZK code.

[0088] In this case, the dimensionless angular frequency is ω„ = 2rf ₀ n / f _ac , f ₀ is the fundamental frequency of the initial sinusoid, A„ is the absorption coefficient, D„ is the dispersion coefficient and p _n is the complex pressure at the w-th harmonic, for n=\..M harmonics. The ambient pressure, temperature and relative humidity were 101.325 kPa, 284.5 K, and 70%, respectively, for this numerical validation. The initial sinusoidal amplitude, p _ac , was 12.09 Pa and had a total of M=100 harmonics. The total relative radius of curvature was 40,000 meters.

[0089] FIG. 5 shows comparisons between the spectral solution 124 (the dashed lines) and the numerical propagation 126 from the LNTE code (solid line) at four different pseudotimes. FIG. 5 illustrates the results of lossy nonlinear propagation of a sinusoid through a homogeneous medium with two relaxation processes. The LNTE code output (numerical propagation 126) is compared to a spectral solution (spectral solution 124) from equation 16 after propagation for a pseudotime of four different shock formation "distances". The upper left plot is a shock formation "distance" of 0.5, pseudotime of 0.1131. The upper right plot is a shock formation "distance" of 1.0, pseudotime of 0.2266. The lower left plot is a shock formation "distance" of 2.0, pseudotime of 0.4525. The lower right plot is a shock formation "distance" of 3.0, pseudotime of 0.6788.

[0090] For this numerical validation the pseudotime variable represents a propagation "distance" variable. The pseudotimes for each plot in FIG. 5 corresponded to a series of lossless shock formation "distances", where one shock formation "distance" in this formulation corresponds to the pseudotime value it would take for a shock to form in the absence of losses. From the top row to the bottom row starting at the upper left corner, the plots correspond to shock formation "distances" of 0.5, 1.0, 2.0, and 3.0, respectively, equal to pseudotime values of 0.1131, 0.2266, 0.4525 and 0.6788. The agreement is very good between the spectral solution and the LNTE output for all four cases.

[0091] A check was also performed to examine the influence of grid discretization on the numerical solution for N-wave focusing. The characteristic acoustic frequency, z _min and z _max were held constant at 8.174 Hz, -1.0, and 1.0, respectively. The total relative radius of curvature, ambient pressure, temperature and relative humidity were also held constant at 48,769 meters, 91.3 kPa, 284 K and 40%, respectively. Table 1 summarizes the grid parameters that were varied to verify the solution converges towards a solution as the grid resolution is increased. The first column of Table 1 was the sample rate in dimensional units of the t -axis and the incoming N-wave which had a characteristic acoustic pressure of 57 Pa. The second column was the number of points in the Έ -axis. The third column was the number of points in the t -axis, which was always a power of two due to the FFT component of the numerical implementation. The next two columns are the <E and dt that pertain to the parameters in the first three columns. Each of the four solutions was determined by iterating to a pseudotime of 5.

[0092] FIG. 6 illustrates a comparison of the front shock from N-wave focusing as the grid resolution is increased to finer values for the location corresponding to the maximum peak pressure of the front shock. The pressure peak and rise profile approach similar values as the grid density is increased. Specifically, FIG. 6 shows a comparison of the solutions from the four different grid resolutions that correspond to the z = 0.154 location and closely zoomed in to the front shock. The z -coordinate was chosen because it was the highest pressure value found for the front shock in the solution at the highest grid resolution. The solutions approached a consistent peak pressure amplitude and rise profile as the grid resolution was refined to finer values.

[0093] FIG. 7 Illustrates a comparison of the rear shock from N-wave focusing as the grid resolution is increased to finer values for the location corresponding to the maximum peak pressure of the rear shock. The pressure peak and rise profile reach similar values as the grid density is increased. Specifically, FIG. 7 shows a comparison of the same solutions for the same acoustic and caustic parameters and the same four grid resolutions as in Table 1 and Figure 6, but is a close-up of the rear shock for the solution at the z = 0.062 location. This z -coordinate corresponds to the location for the highest pressure peak for the rear shock at the highest grid resolution. Again, the solutions converge to a consistent peak amplitude and pressure rise as the grid resolution is improved. Additionally, the "Gibbs- like" oscillations just prior to the onset of the observed shock reduced in both FIGS. 6 and 7 as the grid resolution becomes finer.

Table 1. Parameters for demonstrating grid discretization convergence.

[0094] Lastly, comparisons were made between LNTE solutions with versus without the loss mechanisms included in the numerical computations. The motivation behind the comparisons was to illustrate the influence of the absorption and dispersion terms in equation 1 for an N-wave focusing case. FIG. 8 illustrates a comparison at z = 0.5 of the LNTE code solutions with versus without losses included in the computational solution. Line 128 is the LNTE code with losses included and line 130 is the LNTE solution with the absorption and dispersion terms excluded in the numerical computations. Line 132 is the incoming N-wave scaled and phased according to Equation 7 for the function, F. FIG. 8 is a comparison of the LNTE solution with losses (line 128) versus without losses (130) at z = 0.5 for the same highest resolution solution case from FIGS. 6 and 7. For reference, the incoming N-wave (line 132) is included in the comparison and was scaled and phased according to the incoming function, F, given in Equation 7. The distance of J = 0.5 was chosen for the comparison because it is far enough from the caustic for the front part of the incoming N-wave to not show significant signs of amplification due to the focusing (matching with geometric acoustics), and also far enough from the upper boundary for lossy propagation effects to be observable. In FIG. 8, the LNTE solution with losses included in the numerical computations shows a front shock profile very similar to the incoming N- wave, which also included absorption and dispersion in its propagation. The "lossless" N- wave shock is noticeably steeper than the LNTE solution with losses. Thus, the loss terms in the LNTE maintain a shock profile for the N-wave in its solution that is consistent with the shock profile of the incoming N-wave.

[0095] FIG. 9 illustrates a comparison at z = 0.15 of the LNTE code solutions with versus without losses included in the computational solution. Line 134 is the LNTE code with losses included and line 136 is the LNTE solution with the absorption and dispersion terms excluded in the numerical computations. Specifically, FIG. 9 shows a comparison between the LNTE solution (line 134) and the "lossless" LNTE solution (line 136) at 1 = 0.15 in the vicinity of the maximum peak pressure of the LNTE solution for the same case computed in FIG. 8. Observe that the front part of the waveform has a lower peak pressure but comparable rise times when the absorption and dispersion terms are included in the numerical solution. The comparison at the rear part of the waveform shows the "lossless" LNTE solution has shorter rise times and slightly higher pressure rises at the shocks than the LNTE solution. In examining both FIGS. 8 and 9, the "lossless" LNTE solutions show signs of numerical oscillations, so it is likely that including the loss mechanisms also provides some advantageous suppression of non-physical numerical oscillations. The comparison in FIG. 8 does show that the influence of the absorption and dispersion terms is consistent with the physics demonstrated in the shock profile of the incoming N-wave.

[0096] The SCAMP flight test successfully measured dozens of focus booms from a maneuvering NASA Dryden F-18. The main purpose of the SCAMP flight test was to conduct a flight test that provided measured data for the large-scale experimental validation of focus boom prediction codes. The ground array that captured the focus booms consisted of eighty one microphones linearly spaced over two-miles at the Cuddeback Gunnery Range, California. Trajectory information from the F-18 was simultaneously recorded to log the time-space data of the aircraft for synchronization with the acoustic measurement systems at the ground. Upper-air meteorological data was also captured by sounding the atmosphere with GPS-sonde weather balloons. Ground weather conditions were logged with weather stations. The combination of the aircraft position data, the upper-air meteorological data and ground weather data was used to facilitate focusing predictions at the ground.

[0097] Focus signatures in the SCAMP flight test were very consistent from pass to pass, tending to vary only for adverse atmospheric conditions. One particular test point, "maneuver C", with the aircraft accelerating at a Mach rate of 0.0035 Mach/second and pushing downward at a pitch rate of -0.25 degrees/second, was in the middle of the overall test matrix. This maneuver (flight 1264, pass number 4) was selected for detailed comparisons with the LNTE predictions. Measured acoustic data from the other "maneuver C" flight passes under favorable atmospheric conditions would generally overlay with the selected pass.

[0098] The predicted incoming waveform for the selected flight pass was obtained from PCBoom propagation down to a distance within one diffraction boundary layer thickness normal from the caustic intercept at the ground. The propagation from the near-field to the far-field included the same atmospheric loss mechanisms as Equation (1) and accounts for the ray tube area change resulting from a maneuvering supersonic aircraft relative to an aircraft at a steady supersonic trajectory. The near- field pressure waveform was obtained from CFD at the approximate altitude and Mach number that corresponded to the onset of the ray that was tangent to the caustic at the ground. This corresponds to flight conditions of a Mach number of 1.23, an altitude of 42,000 ft and an angle of attack of 2.3 degrees. The CFD was computed out to three body lengths and served as input to PCBoom. PCBoom computes focus points from a sequence of rays, and computes the radius of curvature by fitting a circle to those points. The rays are chosen such that they lie in a plane normal to the 3-D caustic surface. R _cau for use in LNTE is that at the caustic-ground intercept. This is shown in FIG. 10.

[0099] FIG. 10 illustrates sources of the acoustic input and caustic geometry for the LNTE code. A near- field CFD pressure signature 138 is propagated along a "delta tangent ray" 140 down to a distance of one diffraction boundary layer thickness from the caustic intercept at the ground. The caustic geometry is determined from a caustic tangent ray 142. [00100] The ray curvature, R _ray , necessary to obtain R _tot , is computed by fitting a circle to points along caustic tangent ray 142. The incoming waveform, caustic geometry at the ground intercept and the ground meteorological data were provided to the LNTE code by PCBoom as inputs to compute the focused pressure field at the ground for one of the SCAMP flight test passes.

[00101] FIG. 1 1 is a plot of the incoming waveform that was the predicted output from PCBoom for flight number 1264, pass number 4 from the SCAMP test. The focusing parameters and conditions were a characteristic acoustic frequency of 8.17 Hz, a total relative radius of curvature of 48,769 meters, ambient temperature of 284 K, ambient pressure of 91.3 kPa, and relative humidity of 40%. Since there was no wind at the ground, the wind terms were omitted from the calculations. The computational domain extended from z _max = 1 to z _min = -1 with 12,000 points in the z -direction and 16,384 points in the t direction at a dimensional sample rate of 20,480 Hz. The solution was iterated up to a pseudotime value of 7 (over 23,000 iteration steps). The convergence monitoring for the solution is the example shown in FIG. 3.

[00102] FIG. 12 illustrates the pressure field solution to the LNTE prediction for the NASA SCAMP test, flight 1264, pass 4. More generally, FIG. 12 shows the pressure contour of the LNTE predicted acoustic field. The upper left part of the solution shows the incoming N-wave propagating through to the caustic, the outgoing wave reflecting off of the caustic and propagating to the upper right part of the solution. Observe that the predicted peak pressure does not coincide with the z = 0 location that geometric acoustics predicts, but instead is predicted to be located slightly above it for both the front and rear shocks.

[00103] FIGS. 13 - 16 show a comparison of four of the microphone measurements to their respective LNTE predictions. FIG. 13 illustrates the LNTE prediction at z = 0.97 compared to SCAMP flight test measurement at microphone 71. FIG. 14 illustrates the LNTE prediction at z = 0.67 compared to SCAMP flight test measurement at microphone 67. FIG. 15 illustrates the LNTE prediction at i ^" = 0.15 compared to SCAMP flight test measurement at microphone 60. The microphone at this location measured the highest peak pressure from any of the microphones in the array for this flight pass. FIG. 16 illustrates the LNTE prediction at z = -0.22 compared to SCAMP flight test measurement at microphone 55. [00104] With respect to FIGS. 13-16, the focus predictions have been multiplied by a scale factor of 1.9 to account for ground reflection. In principle, ground impedance effects should be accounted for when comparing predictions with ground measurements. However, for the SCAMP test, the microphones were not placed directly on the ground but instead were placed on ¾"-inch plywood boards that were wrapped by a thin plastic moisture barrier. The measurement site and ground boards presented a hard surface, for which the detailed effects of impedance on sonic boom amplitude tend to be small, so for the current analysis reflection was addressed by the traditional simple factor of 1.9. In FIG. 13, the LNTE prediction of the waveform duration for the incoming N-wave shows good agreement with the microphone measurement. The predicted amplitude is less than the measured data, but the mismatch in pressure levels can be attributed to the influence of atmospheric turbulence, which is not accounted for in the model equation. The predicted outgoing U- wave is in good agreement with the amplitude of the measured U-wave, but the predicted phase is slightly shorter than the measurement. In FIG. 14, the comparison between the LNTE prediction and the measured data shows good agreement for the incoming N-wave. The outgoing U-wave waveform period is slightly under-predicted by the LNTE code, but the amplitude of the U-wave is in good agreement. In FIG. 15, the comparison location corresponds to the microphone that measured the highest peak pressure for any of the microphones in the array for that flight pass. The prediction for the front shock slightly overestimates the peak amplitude. The LNTE predicts the two closely spaced rear shocks similar to what is observed in the measurement, except that the shock spacing predicted by the LNTE code is slightly less than what was measured. FIG. 16 is a comparison of the LNTE prediction to a microphone in the shadow zone. The LNTE prediction shows the amplitudes of positive and negative pressure peaks are less than the measured data, however, the predicted LNTE waveform period closely matches the measured data. It is likely that atmospheric turbulence in the microphone measurement caused some of the observed amplitude perturbations. Overall, there is strong agreement between the LNTE predictions and the microphone measurements at the four locations in the comparisons.

[00105] FIG. 17 illustrates a non-limiting embodiment of a method for determining a pressure time history of a pressure wave originating from a vehicle traveling at or above supersonic speeds when the pressure wave undergoes a focusing due to convergence and intersection with a neighboring pressure wave originating from the vehicle. At step 144, an input is received at a processor. The processor is configured to solve the Lossy Nonlinear Tricomi Equation (described in detail above).

[00106] In some embodiments, the input may comprise a time pressure history of the pressure wave at a predetermined location upstream of the focusing. In some examples, the location is one diffraction boundary layer thickness above the caustic. In some embodiments, the input may comprise a relative radius of curvature between the pressure wave and a caustic. In some embodiments, the input may comprise an atmospheric condition proximate the focusing. In some examples, the atmospheric condition may be measured at a ground surface, may be an ambient pressure proximate the focusing, may be an ambient temperature proximate the focusing; and/or may be a relative percent humidity proximate the focusing. In some embodiments, the input may comprise a distance between the pressure wave prior to focusing and a caustic. In some embodiments, the input may comprise a distance below the caustic. In an example, the distance below the caustic may be selected so as to account for substantially all of the diffraction effects arising out of the focusing. In some embodiments, the input may comprise a desired amount of discretization between a designated upper limit of a region proximate the focusing and a designated lower limit of the region proximate the focusing. In some embodiments, the input may comprise a desired sample rate along a time axis of the pressure time history of the pressure wave in a region of the focusing. In some embodiments, the input may comprise a wind speed. In an example, the wind speed may comprise a wind speed in an X direction with respect to the direction of the propagation of the pressure wave. In an example, the wind speed may comprise a wind in a Z direction with respect to a direction of the propagation of the pressure wave. In some embodiments, the input may comprise one, several, or all of the above mentioned inputs.

[00107] At step 146, the processor is used to solve the Lossy Nonlinear Tricomi

Equation. The Lossy Nonlinear Tricomi Equation includes a variable that relates to an actual atmospheric dissipative effect in a vicinity of the focusing. In some embodiments, the variable relates to heat conduction. In some embodiments, the variable relates to viscosity. In some embodiments, the variable relates to molecular relaxation. In some examples, the variable relates to absorption caused by the molecular relaxation. In some examples, the molecular relaxation comprises relaxation of nitrogen molecules. In some examples, the molecular relaxation comprises relaxation of oxygen molecules. In some examples, the variable relates to dispersion caused by molecular relaxation. In some examples, he dispersion arises, at least in part, out of relaxation of nitrogen molecules. In some examples, the dispersion arises, at least in part, out of relaxation of oxygen molecules. In some embodiments, multiple variables are input corresponding to heat conduction, viscosity, and molecular relaxation, respectively.

[00108] At step 148, an output from the processor is reported. The output contains the solution to the Lossy Nonlinear Tricomi Equation determined by the processor. The solution is reflective of the actual atmospheric dissipative effect. In some embodiments, the output may take the form of a graph. In some embodiments, the graph may depict time along one axis, distance from the caustic along a second axis, and may depict the pressure of the pressure wave in a color coded manner.

[00109] While at least one exemplary embodiment has been presented in the foregoing detailed description of the invention, it should be appreciated that a vast number of variations exist. It should also be appreciated that the exemplary embodiment or exemplary embodiments are only examples, and are not intended to limit the scope, applicability, or configuration of the invention in any way. Rather, the foregoing detailed description will provide those skilled in the art with a convenient road map for implementing an exemplary embodiment of the invention. It being understood that various changes may be made in the function and arrangement of elements described in an exemplary embodiment without departing from the scope of the invention as set forth in the appended claims.