Abstract
A robust algorithm for solving the Bancroft version of the Pochhammer–Chree (PC) equation is developed based on the iterative root-finding process. The formulated solver not only obtains the conventional n-series solutions but also derives a new series of solutions, named m-series solutions. The n-series solutions are located on the PC function surface that relatively gradually varies in the vicinity of the roots, whereas the m-series solutions are located between two PC function surfaces with (nearly) positive and negative infinity values. The proposed solver obtains a series of sound speeds at exactly the frequencies necessary for dispersion correction, and the derived solutions are accurate to the ninth decimal place. The solver is capable of solving the PC equation up to n = 20 and m = 20 in the ranges of Poisson’s ratio (ν) of 0.02
Keywords
Introduction
The Hopkinson bar (HB)1–4 has been traditionally used to measure dynamic transient loads at the near field of high-energy sources. Currently, the HB is also used as a stress transducer for the non-destructive assessment of materials and structures. 5 In contrast, the split Hopkinson bar (SHB)6–16 has been widely used to measure the uniaxial stress–strain and strain rate–strain curves of versatile materials at strain rates of approximately 102–104 s−1.
In both the HB and SHB, a transient elastic wave travels the bar. The detailed profiles of the wave, i.e., the rise time and mode of fluctuation in the plateau region, vary with travel. Therefore, the elastic wave profiles in the HB and SHB must be predicted at an axial position of interest from the measured profile (at a different point); this is called dispersion correction.9,17–24 The distortion of the elastic wave shape with travel mainly results from the fact that a high-frequency wave component that constitutes the overall wave is sluggish compared to a low-frequency component. Therefore, in standard dispersion correction, the phases of a series of wave components (with different frequencies) are corrected at a given axial position. In the literature,25–28 the measured magnitude of the travelling wave is also corrected because the wave magnitude depends on the axial position (i.e., phase) and radial position in the bar.
To perform standard dispersion correction and magnitude correction of the measured elastic wave in the HB and SHB, the sound speeds of a series of elastic waves with different frequencies must be known a priori for a given bar. The necessary sound speed values can be obtained for a range of frequencies by solving the Pochhammer–Chree (PC) equation.29–32
According to Love, 29 Pochhammer (1876) 30 and Chree (1886, 1889)31,32 independently derived the wave equation of motion (vibration equation) for an infinite, isotropic, and elastic bar. The PC equation is valid for a cylinder with finite length if the radius is small compared with the wavelength. 29 Its solution is the relation between c and Λ (see the nomenclature in the Supplemental Material (Section 1) for the notations used herein). This c–Λ relation can be used to obtain the frequency via the equation c = fΛ. While the solutions to the PC equation were sought via various methods, most studies that solved the PC equation for dispersion correction used the root-finding process based on iterative techniques.4,23,24,40–43 This study also employs iterative techniques for root-finding. The challenges in solving the PC equation and the resultant state-of-the-art PC equation solutions are described in the Supplementary Material (Sections 2.1 and 2.2, respectively).
Despite the difficulty of solving the PC equation, obtaining a detailed algorithm to solve the PC equation is challenging. Moreover, the studies reporting the solutions did not reveal the details of their numerical schemes.4,23,24 Because of the cumbersome nature of handling the PC equation to obtain the necessary c–f relation for a specific bar, only a few studies involving dispersion correction4,24,25 have recently solved the PC equation. In the foregoing references,4,23,24 solutions for some considered Poisson's ratios (0.2, 0.287, 0.3, and 0.35) were obtained based on in-house schemes (the details are undisclosed) and were provided in graphical form within limited ranges of frequencies and sound speeds.
The objective of this study is to provide an open-source solver of the PC equation based on a robust algorithm for the following reasons. First, there were no open-source solvers widely used; since the introduction of the PC equation in 1876, the availability of an open-source solver with proven reliability is long overdue. Such an open-source solver will afford researchers with a verifiable tool for solving the PC equation, which will facilitate dispersion correction of measured elastic waves in a (split) Hopkinson bar. It will also encourage program development and formulation of candidate algorithms.
Second, most studies that provided graphical or tabular solutions4,23,24,33–45 (without the detailed solution process) reported sound speeds at a constant a/L interval or at an arbitrary f interval. Therefore, the solution frequencies were unnecessary frequencies for the dispersion correction of the user that necessitated the interpolation of the solutions at the required frequencies. Some examples of the required constant frequency intervals for dispersion correction are summarised in Table 1. Accordingly, this study aims to formulate a solver that yields solutions exactly at the frequencies necessary for dispersion correction.
Fdc_max and dFdc for a range of a values when To = 200 µs, dt = 0.2 µs, and co = 5000 m/s (Ny = fNy/dfdc = Fdc_max/dFdc = 500).
Third, most studies that afforded graphical or tabular solutions4,23,24,33–45 (without the solution process) derived the solutions in limited ranges of approximately 0 ≤ c/co ≤ 2 and 0 ≤ a/Λ ≤ 2 up to n = 4 (corresponding to the range 0 ≤ af/co ≤ 4). The studies that only afforded graphical or tabular solutions4,23,24,33–45 (without the solution process) were also performed in similar ranges. As summarised in Table 1, wider ranges of frequencies (and sound speeds) are necessary for dispersion correction. The solutions should also be derived up to higher orders for a wider range of Poisson’s ratios than those derived in existing studies.4,23,24,33–45 However, studies that determined solutions across wide ranges of frequency, sound speed, order, and Poisson’s ratio are rare. Therefore, this study aims to provide a solver that yields solutions across wider ranges of the foregoing quantities simultaneously (0.02
Finally, as described in the Supplemental Material (Section 3), the PC equation derived by Bancroft35 is advantageous over other versions in developing a solver in that the normalized solutions (af/co vs. c/co) can be derived by inputting only one variable (Poisson's ratio) that describes medium characteristics. However, it is difficult to find the solvers that exploited such an advantage explicitly. Accordingly, this research aims to provide a solver that exploits such an advantage of the Bancroft version.
Pochhammer–Chree equation
Among the many versions of the PC equation, those introduced by Love,
29
Redwood,
44
and Bancroft
35
are introduced in the Supplemental Material (Sections 3.1–3.3) with an emphasis on the benefit of the Bancroft version in the solution process and post-processing stage. To aid in comprehending the dependent nature of the PC function on frequency, sound speed, and wavelength (which eventually aids in developing the scheme for solving the equation), this study expresses the Bancroft version using physics-friendly non-dimensional variables, instead of x,
The physics-friendly non-dimensional variables employed herein are C (=c/co), L (=a/Λ), and F (=af/co), regarded as normalised sound speed, wavelength, and frequency, respectively. According to c = fΛ, the relation between F, C, and L is given by F = CL. This relation is illustrated in the Supplemental Material (Section 3.4). In terms of F and C, the Bancroft version of the PC equation is expressed as (see Section 3.4 of the Supplemental Material)
If an F value is considered as the independent variable, the resultant C value that satisfies equation (1) can be obtained under a fixed value of ν. Once the C (=c/co) vs. F (=af/co) relation is obtained, the c–f relation necessary for dispersion correction can be obtained in the post-processing stage by using the information on the medium, i.e., co and a values. This study solves equation (1).
Solution process
Backgrounds: Fmax and dFdc
As a background for the PC equation solution process, reviewing the necessary values of maximum F (Fdc_max) and F interval (dFdc) in a typical dispersion correction of elastic waves in the (split) Hopkinson bar may be informative. Consider a linearly elastic wave travelling along the axial direction of a circular bar with a radius (a) of 10 mm and a thin-bar sound speed (co) of 5000 m/s (co =
The Fdc_max and dFdc values calculated using the above approach are listed in Table 1 for a range of a values in the case where To = 200 µs, dt = 0.2 µs, and co = 5000 m/s. The list indicates that if the bar radius is increased to 50 mm, an Fdc_max value as high as 25 is required under the considered condition. This study aims to solve an Fdc_max value of up to 30. The table further indicates that the frequency interval necessary for dispersion correction, dFdc, increases in proportion to the bar radius. According to the procedure described in the preceding paragraph, an increase in dt by x times results in 1/x times the value of Fdc_max; increasing the time window (To) by x times yields 1/x times the value of dFdc. A review paper on the dispersion correction of elastic waves in the (split) Hopkinson bar will be published elsewhere.
Scheme
The shapes of the PC function curves in two dimension and PC function surfaces in three dimension are illustrated in the Supplemental Material (Sections 4.1 and 4.2). This study derives the conventionally known n-series solutions which are located on the PC function surface that relatively gradually vary in the vicinity of the roots. It further obtains new-type solutions, named m-series solutions, located between two PC function surfaces with (nearly) positive and negative infinity values (see Section 4.1 of the Supplemental Material).
Based on the illustrated function shapes in the Supplemental Material, a schematic of the dispersion curve (with n

Algorithm (bisection method combined with linear extrapolation) employed in routines X (Xa and Xb) and Y for solving the PC equation for n
In this research, the orders of dispersion curves are counted at intermediate F values: Fstart(1) and Fstart(2). The pre-processing details for determining Fstart(1) and Fstart(2) and the routine for counting the order are described in the Supplemental Material (Sections 5.1 and 5.2). Once the orders are determined, routines
The algorithm of routine
In the low F regime (Figure 1), where the dispersion curve increases extremely rapidly towards infinity with the decrease in F (routine
After routines
The Cg (normalised group velocity) vs. F plots are prepared in the post-processing stage using the obtained (C, F) dataset. The definition of Cg in the C–F domain and the method of obtaining it are available in the Supplemental Material (Section 5.6).
If the solution itself is required for the purpose of plotting the C vs. F diagram, the (F, C) dataset is then sufficient. Furthermore, the C vs. F diagram can be prepared via a graphical method. 45 As mentioned, for the purpose of dispersion correction, the exact Cdc values at exactly necessary Fdc values must be obtained, while the (F, C) dataset generally differs from the (Fdc, Cdc) dataset. The (Fdc, Cdc) dataset is obtained herein via linear interpolation of the (F, C) set, followed by the application of the bisection method (see Section 5.7 of the Supplemental Material).
Results and discussion
Solver verification
The proposed solver was verified first using tabular solutions available in the literature. Figure 2(a) plots the n = 1 solutions available in Bancroft’s study (Table I in 35 ) at varying Poisson’s ratios. In Ref. 35, the solutions in the C–L domain were listed at L intervals of 0.025 or 0.1. Accordingly, this study transformed the tabulated L values to F values in order to obtain the C–F domain solutions; the results are presented in Figure 2(a) as open circles. In this figure, the intervals in Bancroft’s data are irregular in the F axis owing to the transformation from L to F. The figure also shows the solutions obtained using the proposed solver at a constant F interval (dF) value of 0.05 (marked as “×”); the C solutions at a regular F interval are necessary for dispersion correction.

Figure 2(b) plots the n = 2 and n = 3 solutions available in Davies’s study (Table 11
As shown in Figure 2, reasonable consistency is observed between the results reported in the literature (open circles) and those herein (marked as “×”), indicating the reliability of the solutions obtained in this research. In the studies of Bancroft and Davies, solutions (the C value at a given L) are available down to fifth and third (or fourth) decimal places, respectively. In the proposed solver, the tolerance limit (error) in the bisection algorithm is set as 10−10; the solutions obtained herein are accurate down to the ninth decimal place.
Solutions of PC equation
The proposed solver yielded n-series solutions over a range of Poisson’s ratios up to n = 20. The results presented in Figure 3 are for an F value of up to 10 and a C value of 2 or 3; the insets plot the curves in wider ranges. Figure 3 shows points around which the given dispersion curve abruptly changes; herein, these points are referred to as abrupt-change points (ACPs). At a given Poisson’s ratio (Figure 3(a)), the distance between two ACPs in adjacent curves decreases as the order (n) increases. At higher orders, this distance becomes extremely small such that the adjacent curves seem to be connected, but a magnified view (not shown) verifies that they are not. An overly close ACP distance causes a deviation of the curve-tracking process; nearby curves are falsely tracked around the ACPs. Accordingly, more refined dC and dF values are necessary around the ACPs to keep track of higher-order dispersion curves, especially at a small Poisson’s ratio; however, this increases the computational burden. The proposed solver considers the positions of the ACPs in the C–F domain while setting up the dC and dF values for tracking the dispersion curves.

‘n-series’ solutions at Poisson’s ratios of (a) 0.05, (b) 0.15, (c) 0.30 and (d) 0.45; Fcr is explained in Supplemental Material (Section 6.1).
The proposed solver also yielded m-series solutions up to an order of 20 for a Poisson’s ratio of 0.3 (Figure 4). The curves for other Poisson’s ratios (0.05

‘m-series’ solutions at a Poisson’s ratio of 0.30.
The proposed solver obtains the Cg vs. F plot using the matrices of C and F (Figures 3 and 4). The results (Cg vs. F plots) for n-series and m-series solutions at a Poisson’s ratio of 0.3 are available in the Supplemental Material (Sections 6.4 and 6.5) up to order 20.
The solver finally calculates the Fdc and Cdc values by using the (C, F) dataset. The results (Cdc vs. Fdc plots) for n-series and m-series solutions at a Poisson’s ratio of 0.3 (dFdc = 0.02) are presented in the Supplemental Material (Section 6.6) up to order 20. The proposed program writes the (Cdc, Fdc) dataset of each order to an Excel® file (Cdc-Fdc.xlsx).
Overall discussion
The proposed solver does not request the user to input any bar information other than the Poisson’s ratio. Once the Cdc vs. Fdc dataset (non-dimensional quantities) at a given ν is derived using the solver, in the second step, the cdc vs. fdc set (dimensional quantities) necessary for dispersion correction is obtained by the user using further specific information on considered bar (co and a). Therefore, only the influence of Poisson's ratio on dispersion curve shapes needs to be considered in developing a solver when the two-step approach is employed. On the other hand, in the one-step approach[40–43], the c vs. f set (or c vs. Λ set) is directly solved under fixed medium variables (ν, a, and co in the Bancroft version; a, ρ, λ, and µ (or a, ν, ρ, and Ε) in Love and other versions); the solver requests to input the foregoing variable values. Therefore, the influences of these three or four medium variables on dispersion curve shapes must be considered in developing a solver if the one-step approach is to be employed. In this regard, the two-step approach based on the Bancroft version has an advantage in developing a robust PC equation solver.
The values of user variables such as dC and dF are pre-set in the proposed program. These values reflect the characteristics of the complicated shape of the PC function. The validity of the pre-set user variable values in tracking the dispersion curves was checked up to the order of 20 (for both n- and m-series) for a wide range of Poisson’s ratios (0.02
Conclusion
A robust algorithm for solving the Bancroft version of the PC equation was developed based on an iterative root-finding process that refers to the illustrated function shape and combines the bisection, linear extrapolation, and linear interpolation methods. The developed algorithm was implemented in MATLAB®; it is available in the Supplemental Material accessible online.
The proposed solver not only obtains the conventionally known n-series solutions but also a new series of solutions, named m-series. Whereas the n-series solutions are located on the PC function surfaces that vary relatively gradually in the proximity of the roots, the m-series solutions are located between two PC function surfaces with nearly positive and negative infinity values. The solver derives a series of sound speeds exactly at the necessary frequencies for dispersion correction and yields solutions that are accurate to the ninth decimal place. Furthermore, it is capable of solving the PC equation up to n = 20 and m = 20 in the Poisson’s ratio range of 0.02
The validity of the pre-set user variable (e.g., dC and dF) values in the proposed solver in tracking the dispersion curves was verified up to the order of 20 (for both n- and m-series) in a wide range of Poisson’s ratios (0.02
Supplemental Material
sj-zip-2-pic-10.1177_0954406220980509 - Supplemental material for Pochhammer–Chree equation solver for dispersion correction of elastic waves in a (split) Hopkinson bar
Supplemental material, sj-zip-2-pic-10.1177_0954406220980509 for Pochhammer–Chree equation solver for dispersion correction of elastic waves in a (split) Hopkinson bar by Hyunho Shin in Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science
Supplemental Material
sj-pdf-1-pic-10.1177_0954406220980509 - Supplemental material for Pochhammer–Chree equation solver for dispersion correction of elastic waves in a (split) Hopkinson bar
Supplemental material, sj-pdf-1-pic-10.1177_0954406220980509 for Pochhammer–Chree equation solver for dispersion correction of elastic waves in a (split) Hopkinson bar by Hyunho Shin in Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science
Footnotes
Declaration of Conflicting Interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This study was financially supported by two grants from the National Research Foundation of Korea under contract Nos. 2019R1H1A2039686 and 2020R1A2C2009083, funded by the Ministry of Science and Technology (Korea). The author expresses gratitude to Hyekyung Park for technical support.
Supplemental material
Supplemental material for this article is available online.
References
Supplementary Material
Please find the following supplemental material available below.
For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.
For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.
