{'fraction_non_alphanumeric': 0.030961791831357048, 'fraction_numerical': 0.010910737812911726, 'mean_word_length': 4.615953757225434, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
How the three-dimensional geometry of computational domain(s) affects the accuracy of non-reflective boundary conditions in acoustic simulation

Describing and simulating acoustic wave propagation can be difficult and time consuming; especially when modeling three-dimensional (3D) problems. As the propagating waves exit the computational domain, the amplitude needs to be sufficiently small otherwise reflections can occur from the boundary influencing the numerical solution. This paper will attempt to quantify what is meant by `sufficiently small' and investigate whether the geometry of the computational boundary can be manipulated to reduce reflections at the outer walls. The 3D compressible Euler equations were solved using the discontinuous Galerkin method on a graphical processing unit. A pressure pulse with an amplitude equivalent to 10% of atmospheric pressure was simulated through a modified trumpet within seven different geometries. The numerical results indicate that if the amplitude of the pulse is less than 0.5% of atmospheric pressure, reflections are minimal and do not significantly influence the solution in the domain. Furthermore, the computational region behind the bell can be neglected greatly reducing the required memory and run time.


Introduction and Background
Systems of partial differential equations (PDEs) which model acoustic wave propagation can typically be written as a hyperbolic system of conservative laws. Numerically solving these equations of motion can be expensive because discontinuities can arise in the solution, even if the initial wave profile is smooth [1]. A common approach when numerically simulating small-amplitude wave propagation (i.e., acoustic waves) is to use a finite volume (FV) method; whereas for finite-amplitude waves (i.e., nonlinear acoustic waves), finite element (FE) methods are more standard [2]. A high-order accurate method which combines aspects from both FV and FE methods is the discontinous Galerkin method (DGM) which was developed by Cockburn and Shu [1] [3]. This method is particularly useful for modeling hyperbolic systems because numerical fluxes defined along the cell interface guarantees stability and local solvability, i.e., at each time step, only information about the neighbouring cells are required. DGM can also accommodate unstructured meshes, adaptive refinement, [4] and is open to parallel implementations [5].
To ensure acoustic simulations can be completed in a reasonable amount of time, the computational domain needs to be truncated in an efficient manner. A computational boundary, denoted by ∂Ω, is defined to construct a computational domain, denoted by Ω. In order to have a well-posed problem, boundary conditions (BCs) must be defined on ∂Ω [6]. In aeroacoustic simulations, this is one of the most difficult endeavors [8]. The outflow BCs must be imposed such that the acoustic waves can propagate directly out of the domain with minimal reflections. Such BCs are typically referred to as pass-through or non-reflecting [7].
If the BCs defined on ∂Ω are poorly prescribed, spurious waves may be generated and reflected from the boundary [7]. This can specifically be an issue for nonlinear problems or higher-order schemes due to the low dispersion and dissipation properties of the DGM [8]. Spurious waves are typically classified as "unphysical and supposedly make elements of high-order useless for accurate computations" [2]. Physically, they represent energy from the far-field solution propagating into Ω [10]. Though in [2], it has been suggested by Mulder that for linear systems, the spurious modes contribute "to the numerical error that behaves in a reasonable manner, and that higher-order elements can be more accurate than lower-order elements." Nonetheless, over the years several non-reflecting or pass-through BCs have been implemented to better model far-field solutions for various problems [9].
One of the main approaches to prescribe far-field BCs is by using characteristics or Riemann invariants [7]. This approach however is intrinsically best suited for one-dimensional (1D) problems. Thomas [11] was one of the first to use characteristics to examine multi-dimensional non-reflecting BCs for time dependent hyperbolic systems. Such an approach is useful for disturbances that propagate through the boundary at a 90 degree angle. However, using characteristics is not very efficient for problems where slanted perturbations (i.e., oblique incidences) hit the boundary [7].
Another technique is to use asymptotic expressions of the equations of motions for the far-field solution. Considering the asymptotic expansion allows for a more continuous transition as waves propagate through the boundary. Bayliss and Turkel [10] considered this approach for non-reflecting BCs for far-field radiation. Others have further extended such work. For instance, Bogey and Bailly formulated non-reflecting BCs for three-dimensional (3D) geometries for the linearized Euler equations. For the radiation condition however, linearizing the BC usually is not accurate enough when modeling physical wave phenomena [12]. To deal with this, a buffer zone or absorbing layer around Ω can be considered as Colonius presents in [12]. In this buffer region, damping is applied to the far-field solution in which the ghost cells or points calculate the pressure gradient normal to the wall [12]. Regardless of the chosen approach to impose pass-through BCs, if the acoustic waves which reach ∂Ω have a sufficiently small amplitude, the disturbance should leave Ω with little to no reflections.
Another factor that influences the accuracy of pass-through BCs is the geometric shape of the boundary itself. In fact, one of the most challenging aspects when imposing pass-through BCs is how to deal with corners effectively [7]. In 3D this can be especially difficult since reflections can occur from multiple directions [13]. Therefore, depending on the problem being studied, using a sphere for instance as the ∂Ω shape rather than a cube or rectangular prism may be more beneficial. Although, unless curved elements are considered at the outer boundary, it is likely that the mesh may need to be refined potentially slowing down run time. Furthermore, deciding on the appropriate geometry can be challenging for acoustic simulations through ducts (for instance, in musical instruments) because the "position and geometry of the domain boundary is governed by the shape of the duct" [14].

Overall Objectives
It is of interest to the author to study the production and propagation of finite-amplitude sound waves in brass musical instruments, specifically the trumpet. To simulate such wave propagation in 3D within a timely manner, it would be optimal to construct the smallest possible Ω. In addition, manipulating the geometric shape of ∂Ω could further reduce run time while producing a more accurate solution. However, we would need to quantify how small the parameters of the outgoing waves must be upon interacting with ∂Ω to ensure that the solution inside Ω is not influenced by any reflections from the boundary as the waves exit the domain. This is precisely the the primary objective of this paper: to quantify how 'small' the amplitude of the outgoing waves must be. The second purpose of this paper is to determine whether there is an optimal geometric shape for the boundary of the computational domain, i.e., if reflections from ∂Ω can be reduced by manipulating the geometry so that the angle of the outgoing waves is small relative to the boundary.
The numerical method that will be used for this investigation will be the DGM due to its usefulness and gaining popularity in aeroacoustic simulations. For instance, the DGM has been successfully used for acoustic simulations on the trumpet and bassoon [15] [16] [17] as well as acoustic pulse studies [15] [18] [19]. Since long simulations need to be run in 3D, serial computation is too slow. Instead, a graphical processing unit (GPU) implementation of the method is used. The GPU implemention used is described in [20] and takes advantage of the numerical method's parallelization features. It was done using NVIDIA's Compute Unified Device Architecture (CUDA). All simulations presented in this paper were run on the NVIDIA GeForce 750 Ti 2 GB graphics card. A brief summary of the numerical scheme is provided in Appendix A and a more detailed description can be found in [20] [23]. Numerical experiments presented in Section 3 were computed using a second order accurate linear approximation in space and verified by running the same tests with a quadratic basis. Due to the very fine meshes used, the results were similar. Runge-Kutta (RK2) discretization in time was used with the local Lax-Friedrichs Riemann solver.

Geometric Shapes Considered
Since the trumpet is one of the main instruments being studied by the author, a modified geometry shown in Fig. 1 will be considered in various Ω. The geometry of the shortened instrument was made to accurately model the slowly increasing diameter of the bell and corresponding cylindrical tube leading up to the flare. It was important to have a realistic representation of the trumpet bell because as mentioned in [14], the optimal geometry of the domain boundary will be governed by the shape of the wavefront which exists the duct. The portion of the waveform that radiates from the flare mostly contains higher frequency components and is roughly spherically symmetric with a curved wavefront [24]. However, there will be a bias as the energy leaves the flare: more will be focused along the central axis rather than along the edge of the flare. The computational trumpet will be modeled with various domains. In particular, two sets of geometries for Ω will be considered: Type A: Domains that completely encapsulate the modified trumpet.
Type B : Domains that only have an area directly in front of the modified trumpet. The Type A geometries include a cube, a sphere, a cone, and a spherical cone. Whereas, the Type B geometries consist of a rectangular prism, a rotated hemisphere and a cone all to the right of the trumpet flare. The goal of examining these different volumes for Ω is to address the following inquiries: 1. Is it advantageous to simulate a curved acoustic wavefront in a domain with spherical symmetry? 2. Is it necessary to consider Ω to the left of the flare, or can it be neglected since only a small portion of energy propagates behind the bell?
With the same modified trumpet placed in the same position for all the Ω considered, as shown in Fig. 2, seven different meshes were generated using the GMSH software. All meshes consist of tetrahedral elements with adaptive element sizes to accurately resolve the geometric features of the modified trumpet. The points that were used to trace out the trumpet boundary are depicted in Fig. 1. The trumpet shape in each domain was created as follows: the cylindrical bore region connects two points with a line, while the remaining points outlining the flare were interpolated using cubic splines. A rotational extrusion about x = 0 was applied to the curve outlining the entire trumpet shape (above the central axis) to create the full 3D instrument. The meshes were scaled so that one unit on the axes represents 1 cm. The meshes constructed from the Type A and Type B geometries can be seen in Fig. 3 and Fig. 4, respectively. The names of each computational domain and corresponding number of elements can be found in Table 1.   2 Numerical Set-Up

Euler Equations
The mathematical model used in this paper was based from the 3D compressible Euler equations. The nonlinear equations are given by where ρ is the density of air; (ρu, ρv, ρw) are the momenta in the x, y and z direction, respectively; p is the internal pressure; and E is the total energy. The equation of state for the ideal gas connects E to the other variables and closes the system For air, the specific heat ratio γ ≈ 1.4. The system (1) describes the motion of an inviscid, compressible gas. Hence, the model neglects viscosity of the air and wall losses in the trumpet. To solve equations (1 -2), initial and boundary conditions need to be specified.

Initial Conditions and Pass-Through Boundary Condition
It was assumed initially that the flow was at rest. The remaining variables are scaled from physical values to ones more convenient for computation. In particular, the speed of sound c 0 , which is approximately 343 m/s in air, and atmospheric pressure p 0 , which is 101,325 Pa (i.e., one atmosphere), are scaled to be equal to 1.
Assuming that the flow is isentropic, i.e., c 0 = γp0 ρ0 , the initial density should then be taken to be 1.4. In summary, Initial State: Along the computational domain ∂Ω, pass-through boundary conditions are used in which the ghost state is prescribed to be the initial state given in (3). If the computational domain is large enough, the amplitude of the wave should be small enough to completely pass through ∂Ω without being reflected back into the domain.

Solid Wall Boundary Conditions
On the inner and outer walls of the trumpet, excluding the left vertical boundary, reflective boundary conditions were prescribed. A ghost state is specified so that the normal velocity is reflected with respect to the wall, i.e., taken with a change of sign. The density, pressure and tangential velocity are unchanged from the corresponding value inside the cell. In summary, Reflective Condition:

Inflow Conditions
At the left vertical boundary of the bore, the ghost state is specified to the be the inflow condition which will now be discussed. For all numerical test cases considered, a simple pressure pulse is generated at the left boundary of the modified instrument and is defined as  As this pulse propagates through the modified trumpet, a portion of the incident wave will be reflected before and within the bell and the other portion will be transmitted from the flare into Ω. As the transmitted wave propagates toward ∂Ω, its amplitude will drop off inversely proportional to the distance traveled. The purpose of considering such a large pulse is to ensure it does not completely diffuse before reaching the computational boundary.
The wave entering the bore is assumed to be planar. To relate pressure and velocity at the inlet boundary, the 1D expression derived from linear acoustic theory is used, i.e., Note that this linearization is used only for the mouthpiece boundary condition, not inside the domain where the velocity is governed by (1 -2). The density is prescribed by assuming the adiabatic relation between pressure and density from compressible flow theory, i.e., where C = γ is the proportionality constant [25]. Therefore, in summary, the dimensionless boundary condition at the mouthpiece of the computational trumpet is given by Inflow Condition:

Numerical Results and Discussion
The numerical results of solving (1) with the initial and boundary conditions on all the computational domains mentioned above will now be examined. For each simulation, the generated pressure was sampled at various positions throughout the computational domain. Although all the sensors' data will not be shown here, Fig. 6 (left) illustrates the placement of the sensors on the 3D coordinate grid. The junction of the flare and room is positioned at x = 1.48 m and the central axis of the flare is at (x, y, z) = (1.48 m, 0 m, 0 m). The flare reaches the height of y = 0.06738 m. For convenience, we will name the sensors that will be examined here. The position of these specific sensors are depicted in Fig. 6 (right) and the corresponding positions relative to the bell and boundary are listed in Table 2

Name
Location from bell & boundary (x-axis) Location from bell & boundary (y-axis) sensor 1 1cm from bell, 101cm from boundary 99cm from bell, 1cm from boundary sensor 2 100cm from bell, 2cm from boundary 0cm from bell, 100cm from boundary sensor 3 101cm from bell, 1cm from boundary 0cm from bell, 100cm from boundary sensor 4 37cm from bell, 65cm from boundary 57cm from bell, 43cm from boundary For both geometry types, the portion of the pulse that will propagate out of the bell along the central axis has a maximum and minimum pressure peak of approximately 400 Pa and -400 Pa, respectively. This corresponds to a deviation from atmospheric pressure by approximately 0.3947%. However, only half that amount of energy travels vertically from the bell (i.e., from very top or bottom of the bell) compared to that along the central axis. Plots illustrating this in the time domain can be seen in Appendix B in Figs. 16  maximum and minimum pressure peak of approximately 10 Pa and -10 Pa, respectively. This is equivalent to 0.009869% variation from atmospheric pressure. With these values in mind, we will now examine the numerical results near ∂Ω. In Figs. 7 and 8, the frequency spectra of the simulation results for the Type A and Type B geometries, respectively, at sensors 1 -3 are shown. With the exception of Ω B-half-cone , we observe that the spectrum for each Ω at sensor 1 and sensor 3 are almost identical, especially for components less than 2200 Hz. The entire spectra however at sensor 1 and sensor 3 for Ω A-sphere and Ω B-half-sphere match for all the frequency components (which is not surprising because of the spherical symmetry). Ω A-cube and Ω B-half-cube also give extremely similar results at these sensor locations with slight deviations for a few of the components between 3000 Hz -4000 Hz. We see the most discrepancy in the spectra for Ω A-scone , Ω B-half-cone and Ω A-scone (though Ω A-scone is more accurate). Figures 7 and 8 also depict the spectrum at sensor 2, which is 2 cm from the rightmost boundary along the central axis. Comparing sensors 2 and 3 can give some indication whether or not reflections occurred from ∂Ω. In particular, for Ω A-sphere and Ω A-scone , we observe that the spectra sampled 1 cm from the boundary is roughly 1-3 dB larger for components above 300 Hz relative to the spectra sampled 2 cm from the boundary. This implies that a portion of the pressure pulse was reflected from the boundary in the mentioned domains. Potentially this observation could be due to refinement since curved elements along the boundary were not considered. Reflections do not seem to be evident from any of the Type B geometries. However, the harmonic components corresponding to the Ω B-half-cone solution are almost 20 dB higher than all other solutions. We will discuss some possible reasons for this shortly. In Fig. 9, a comparison of the solution curves obtained from sensors 2 -4 are given for Ω A-sphere and Ω B-half-sphere (top left); Ω A-cone and Ω B-half-cone (top right); Ω A-scone and Ω B-half-cone (bottom left); and Ω A-cube and Ω B-half-cube (bottom right) is given. The Type A geometries are denoted with a solid line, whereas Type B is illustrated with a dotted line. We observe that for the most part, the data for both geometries match well for the sphere/hemisphere and cube/rectangular prism. However, Ω B-half-cone does not align with Ω A-cone or Ω A-scone .
Reflections seem to occur from the boundary of the domains with spherical symmetry, though this does not appear to greatly influence the solution inside Ω. This observation can be made by examining the data sampled inside the modified instrument. If the simulated pressure for all Ω is consistent, this implies the numerical solutions were not influenced by any reflections that may have taken place. Thus, the frequency spectrum of each pressure pulse was sampled at the left vertical wall inside the bore after the disturbances were reflected from the bell. These results are shown in Fig. 10. In particular, solutions of Ω A-sphere and Ω B-half-sphere (top left); Ω A-cone and Ω B-half-cone (top right); Ω A-scone and Ω B-half-cone (bottom left); and Ω A-cube and Ω B-half-cube (bottom right) where juxtaposed. Although the Ω's with flat boundary walls (i.e., plane surfaces) match slightly better than the Ω's with curved boundary walls (i.e., ruled surfaces), there is little variation overall.
Next, we will examine the numerical solutions at the bell for each domain. If reflections from the boundary influence the numerical solution in Ω, it will be more evident at the bell rather than inside the smaller cylindricaltube portion. Figure 11 plots all the simulation results for both geometry types at the midpoint of the bell (top left), 1 cm outside the midpoint of the bell (top right), at the top of the bell (bottom left) and 1 cm outside the top of the bell (bottom right). We notice that all the simulation results seem to align with each other rather well. This should not be surprising for stable simulations because one would expect that the simulated pulse propagating out of the flare should be the same for each Ω. We do however notice that the solution for Ω B-half-cone deviates rather drastically from the other solutions for frequency components less than 500 Hz. This indicates that the solution in Ω B-half-cone was influenced by reflections from the domain boundary. And we will now discuss some possible reasons why such reflections took place in Ω B-half-cone . The angle of the cone in Ω B-half-cone which extends from the flare (relative to the top of the bell) is 45 degrees. In the top right and bottom left plot of Fig. 9, we saw that the numerical solution of Ω B-half-cone at sensor 4 (which is only 37 cm outside the bell) was approximately 40 dB higher than any of the other solutions. However, 100 cm outside the bell, the harmonic distribution was only 20 dB higher. This indicates that the amplitude of the pressure pulse in Ω B-half-cone was increasing before it reached ∂Ω. The only way for this to happen is if a portion of the disturbance was being reflected along the side of the cone. If this was the case, reflections would occur in all directions due to spherical symmetry. It would then seem reasonable that the amplitude increases so quickly. If the cone instead followed the angle of the flare, such a computational domain may potentially give the best results because the propagating pressure would interact with the side of the cone at (roughly) a 90 degree angle. By the time the pressure pulse would reach the rightmost boundary, it would be equivalent to a plane wave propagating into a plane surface, i.e., it could propagate through the wall without reflections occurring.

Conclusion
A comparison of all the numerical results obtained from simulating a pressure pulse in the Type A and Type B geometries is shown in Fig. 12 (in addition to Fig. 11). In particular, the top left and top right plots in Fig.  12 shows all the Type A results at sensor 3 and at the left vertical boundary inside the bore, respectively. The bottom left and bottom right graphs illustrate the corresponding results for the Type B geometries. With the exception of Ω B-half-cone , we see little variation between the two different computational domain types. Therefore, it does not seem necessary to construct Ω to encapsulate the entire instrument. This is very encouraging because this significantly reduces the required memory and run time. In fact, from examining Table  3, we can see that the number of cells and memory required is almost cut in half if the domain to the left of the flare is neglected. Furthermore, if the amplitude of the pressure pulse is less than 0.5% of atmospheric pressure, and the mesh is sufficiently refined, there is little to no reflections from the boundaries of the computational domain. More specifically, the reflections that do occur do not seem to greatly influence the solution inside of Ω.
However, it is inconclusive from the results presented here that considering a computational domain with spherical symmetry is advantageous. But, the spherically symmetric domains that were considered were not well refined (i.e., the refinement was the same for all geometries for easier comparison) and curved elements along ∂Ω were not used. If the meshes for Ω A-sphere , Ω B-half-sphere and Ω A-scone were not so coarse, better results may have been produced. Further simulations will have to be done using a GPU with more memory.