Abstract
The intensive pulse sound wave can be generated by the underwater plasma sound source (UPSS) based on the discharge of the underwater high voltage. The distribution of the sound field is prominently nonlinear. In this paper, the sound field of the intensive UPSS is described by the integral two-dimensional axisymmetric unsteady Euler equations firstly. In order to solve the Euler equations numerically, an optimized fifth-order symmetric WENO (weighted essentially nonoscillatory) method based on the three templates is proposed which is called WENO-SYM3. Without increasing the number of candidate templates, a new symmetric template structure can be obtained by expanding the second template and shifting the third one backwards for one space. The method is validated through numerical examples and experiments, and the results show that WENO-SYM3 has a high distinguished accuracy; meanwhile, its nonphysical oscillations are not obvious. The experimental results are basically the same as the calculation results, and the maximum error is around 3%.
1. Introduction
After the development in recent half a century, the technology of high intense pulse sound wave, generated by the discharge of underwater plasma based on the liquid-electric effect [1], has been applied in many fields, such as extracorporeal shock wave lithotripsy on health [2], nontraditional machining technology in industry [3], and sewage disposal in environmental [4]. Nowadays, the underwater plasma sound source (UPSS) also serves as the intense sound source, with the advantages of high emission acoustic power, high peak energy, wide frequency band, and long distance propagation, which can be used in ocean geological survey, underwater target detection [5], wideband acoustic interference, and long-range underwater acoustic communication.
The sound field distribution of a large amplitude pulse sound wave generated by UPSS is prominently nonlinear. Usually, the Schlieren method of photography is used to observe the propagation of intense sound field by the means of high-speed camera or pressure sensors which are used to measure the sound pressure of the sound field [6]. Unfortunately, due to the more expensive equipment and many different kinds of suitable equipment needed, as well as the high sensitivity of equipment required, it results in a big cost of the experiments. Therefore, the modeling and calculation of the intense sound field distribution of underwater plasma sound source have become especially important [7].
In this paper, an optimized fifth-order symmetric WENO (weighted essentially nonoscillatory) method (WENO-SYM3) based on the three templates is proposed [8–11]. The basic idea of WENO is to obtain a higher approximation by a linear combination of low order numerical fluxes and nonoscillation across discontinuities. The sound fields of the underwater plasma intense sound source are simulated by WENO-SYM3 scheme, and some related experiments to verify the accuracy of the simulation results are conducted.
2. Brief Introduction of WENO Scheme and WENO-SYMOO Scheme
Generally, the sound field modeling of the underwater plasma intense sound adopts Euler equations. The integral Euler equations can be expressed as in which is the flow field parameters and is the flux vectors. and can be expressed in 2D form as where is the flow density, is the flow pressure, and are the velocity components, is the energy, and is the change-resisting velocity which can be expressed as In the above equation, is the unit vector along the x-axes and is the unit vector along the y-axes.
Since there is no analytic solution of the Euler equations, numerical method can be used to solve it. As a high-resolution numerical method, WENO has been more and more widely used to solve the Euler equations. The alternative template of the normal WENO scheme [11] developed by Jiang and Shu is shown in Figure 1(b). The template of the symmetric WENO scheme optimized by the bandwidth is shown in Figure 1(a).
To improve the performance of the WENO scheme, Weirs [12] chooses two very different ways to optimize the WENO scheme proposed by Jiang and Shu. Considering the slight upwind characteristics taken by Jiang and Shu while they selected a candidate template, the first optimization method adds a candidate template S3 on the existing templates to get the zero dissipation effect theoretically. Another optimization method with bandwidth might choose a right weight to maximize Lele’s bandwidth-resolution efficiency index [13]. By comparing the two optimization methods in [14], it was found that the WENO-JS method often suffers from nonphysical oscillations in discontinuous points. Therefore, the first optimization method can be chosen as an optimization method of bandwidth, such as WENO-SYMOO.
Flow field parameters can be stored in the center of the control volume if the cell-centered method is chosen. According to the sum of flux density on the control volume boundary, the derivation of the shift can be calculated. By considering the following examples focusing on the right side of the control volume boundary, on which the flow field parameter is , the construction process of the WENO scheme could be analyzed.
Assuming that the initial template is , after four times expansion and optimization, the candidate templates of the symmetric WENO scheme are shown in Figure 1(a). can be expressed as in which is equal to 2 in WENO-JS scheme, is equal to 3 in WENO-SYMOO scheme, is the candidate template number, and and are the interpolation function and the weight coefficient, respectively. The expression of is can be expressed as in which is where is a smoothness measurement coefficient of the flux function, which can be expressed as
The definitions of these coefficients , , and in (5), (7), and (8) can be found in [8–11, 14].
3. A Novel Symmetry Optimization Method WENO-SYM3
Optimized symmetry of WENO-SYMOO scheme is realized by adding a candidate template. In this section, a novel optimized fifth-order symmetric WENO method based on the three templates (WENO-SYM3) is proposed, which easily realizes optimized symmetry of the algorithm (as shown in Figure 2) and improves shock capture capability by adding the number of elements in the template S1 without changing the number of candidate template.
is calculated by (4) based on the template of Figure 2. In (4), and the interpolation function can be expressed as
The smoothness measurement coefficient can be written as
The coefficients are given by
Finally, can be expressed as in which the weight coefficient is calculated by (6). As for the flow field parameter, , of the left side of the control volume boundary, a similar method could be used. and are symmetry at . So we can simply replace in (9) and (10) with in order. Thus, the expression for are finally obtained.
After the flow field parameters of the control volume boundary are calculated, this problem could be translated into solving the Riemann problem, and this method was firstly proposed by Godunov [13]. The flux on both sides of some control volume boundary can be obtained by accurately solving the Riemann problem. Similarly, the flux of other control volume boundary also might be obtained. Thus, the variations of the flow field parameters of the control volume stored versus time can be deduced. In this way, the pressure distribution, density distribution, and velocity distribution in the whole flow field could be finally calculated. Because of the calculation and storage complexity in solving the analytic solution of the Riemann problem, some approximate solutions of the Riemann problem are usually used in practical application. Commonly used methods include Roe scheme [15] and HLL/HLLC approximate Riemann solver. Roe scheme has been widely adopted for high precision and good shock capturing skills. In this paper, we have chosen the Roe scheme as the method to solve the approximate Riemann solution, that is, in which is called the Roe matrix and expressed as where
Flow field parameters, such as the flow density and the velocity components and , can be replaced by average variables of Roe when calculating the Roe matrix as in which the subscripts “ and ” represent the left side flow field parameters of the control volume boundary and the right side flow field parameters, respectively.
As for the time discretization, the TVD Runge-Kutta method with three order accuracy is more popular [16]. The expressions are as follows:
4. Validation of Numerical Method’s Effectiveness
In order to validate whether the numerical method is correct or not, we need to compare the numerical result with the exact solution and to find some difficult numerical examples (close to real conditions) to test the effectiveness of the performance of the numerical method. It is more helpful to show the superiority of this algorithm. The correctness and effectiveness of numerical method are validated by Lax shock tube problem and Shu-Osher problem. For the first one, the exact solution existed. Shu-Osher problem contains the fine structure of shock wave and smooth flow field, which is a simple model of shock wave and turbulence effect. The effectiveness and computational accuracy can be tested by this example.
4.1. Lax Shock Tube Problem
The mathematical model of Lax shock tube problem is as follows:
The initial state is
WENO-JS method, WENO-SYMOO method, WENO-SYM3 method, and MUSCL method are used to solve this problem in simulation with grid number 400, CFL = 0.1, and total time . The density and partially enlarged figure in flow field are shown in Figure 3. WENO-SYM3 method is compared with WENO-SYMOO method, and the computational complexity of different numbers of grids is analyzed. The comparison of results before and after optimization is given by Table 1.
(a)
(b)
In Figure 3, whether for WENO method or for MUSCL method, the results obtained are in good agreement with those obtained by the theoretical solution, and there is no any oscillatory phenomenon. There is a phenomenon of smoothing and intermittent to some degree at the point where density exists intermittent. However, as for the explicit solution, the calculating error of WENO method is less than that of MUSCL method. Compared to calculation results in Table 1, WENO-SYM3 method is better than other methods and improves the computation speed.
4.2. Shu-Osher Problem
The initial conditions of one-dimensional Euler (18) are as follows:
Computational domain is [−5,10]. This example is a model that contains the fine structure of shock wave and smooth flow field, which is a simple model of shock wave and turbulence effect. So shock wave and small flow structure appear together as it contains dextral positive shock wave entering flow field with density fluctuate. There are high requirements for resolution of format and diffusion of value in numerical simulation. In this paper, speeds in different positions are given when grid number is 300, CFL = 0.1, and , as shown in Figure 4. Figure 4(a) shows the global figure of simulation results. Figure 4(b) gives the corresponding partial enlargement of global figure. The “explicit solution” is identified as the result when grid number is 1800. Among them, WENO-SYMOO is the symmetric WENO algorithm described in [14]. WENO-JS is shown in [11], and WENO-SYM3 is a new method introduced in this paper.
(a)
(b)
For one-dimensional shock tube problem, there is a precision difference in the results of the two methods, but it is not too big. As shown in Figure 4, the calculating result of WENO method for the Osher problem is obviously better than that of MUSCL method. Meanwhile, there are also differences within WENO methods mentioned before. The result of WENO-SYMOO method is better than that of WENO-JS based on the number of discrete points in this paper. However, the differences are not obvious. Compared to the other two methods, the result of WENO-SYM3 proposed is much closer to theoretical solution. The computation speed of the WENO-SYM3 method equals to that by the normal fifth-order symmetric WENO method. But it is better than the symmetric WENO method that is formed by increasing the number of the templates.
5. Propagation Calculation and Experimental Validation of UPSS
In consideration of the features of UPSS studied in this paper, namely, high pressure peak and narrow pulse width, we choose WENO-SYM3 method as numerical calculation method to calculate propagation characteristics of UPSS field by comparison of the two examples above, and certain related experiments were designed to validate the accuracy of the results.
5.1. Experimental Setup
A schematic diagram of the experimental setup is shown in Figure 5. The experimental system mainly consists of high-voltage charging system, discharge circuit, trigger circuit, charge-discharge control system, discharge electrode, and pressure sensors. The discharge electrode is located in the center of a water tank, which is filled up with tap water.
The function of power-storage capacitor C is to store and discharge the voltage power. The charge-discharge control system would control the trigger circuit to generate trigger pulse, which would break down and turn on the spark gap switch G. Once the spark gap switch G is triggered, the discharge electrode S will be broken down. Then, a mass of electrical power is dissipated into the discharge channel. As a consequence, the sound wave is generated and measured at the main axis by pressure sensors. In the same time, the data is recorded by using a multichannel oscilloscope.
5.2. Simulation Calculation
Due to the symmetry of UPSS propagates in open water, we choose a quarter of physical region as the calculating region, in which , , and the number of grids was 500 × 500 with equidistant division. The origin pulse sound source was located in [0,0] as shown in Figure 6.
The nephograms of the sound field calculated at different times are shown in Figure 7.
(a)
(b)
The attenuation rule of pulse sound propagation in x-axis is shown in Figure 8.
5.3. Experimental Validation
In order to validate the correctness of simulation results, related experiment was designed. Experiment was carried out in a water tank of 3 × 1 × 1 m. The capacitance of capacitor was 5 μF. Discharge voltage was recorded as 20 kV, and gap distance of discharge electrode was 4 mm. The pressure of direct wave was measured by PCB pressure sensors in acoustic axis at a distance of 0.877 m from the discharge electrode. Both experimental wave and simulation wave are plotted in Figure 9, from which we can find a good agreement between the results. The main reasons for the difference of the experimental and simulation waves are the simulation results which can be affected by accuracy rating of the scheme, discreteness of experiment measurement, and the sound field grid division.
The capacitance of capacitor was chosen as 5 μF. Discharge voltages were 12 kV, 16 kV, and 20 kV, respectively. The gap distance of discharge electrode was 4 mm. The pressure of direct wave was measured by three PCB pressure sensors in acoustic axis at the distances of 0.24 m, 0.48 m, and 0.877 m, respectively, from the discharge electrode. The result obtained by multimeasurement average in the same discharging voltages is shown in Figure 10.
In the experiment, the pressure at the sound source varied with the discharge voltage. Thus, the simulation of different discharge voltages could be realized by setting different initial pressures in simulation. The simulation results and experimental data of direct wave amplitude at the distances of 0.24 m, 0.48 m, and 0.877 m were listed in Figure 11. It shows that the simulation results were consistent with experimental result in the condition that rectangular grid was used, and the maximum error was about 3%.
6. Conclusions
An optimized fifth-order symmetric WENO method (WENO-SYM3) is proposed based on the two-dimensional axisymmetric unsteady Euler equations in order to calculate UPSS sound field. Compared with the common WENO method and symmetry method in [14] by using typical numerical example, the WENO-SYM3 has a high distinguished resolution and accuracy; also, its nonphysical oscillations are not obvious.
In the last of the paper, UPSS sound field was calculated by WENO-SYM3 method, and related experiment was conducted to validate the results of calculation. The results show that calculation results are consistent with experimental data with the maximum error about 3%.
Conflicts of Interest
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (61701529, 61471298) and the Science and Technology on Electronic Information Control Laboratory Foundation of China (N2016KC0029).