Abstract
The presence of nonlinear discontinuities in a rockmass makes the stress wave propagation rules in a continuous medium not applicable. In an attempt to reveal the transmission laws of a one-dimensional P-wave across a single nonlinear joint in a rockmass, a recurrence equation is deduced using a semianalytical and seminumerical method from the nonlinear wave equation by introducing the static Bandis−Barton (BB) model for a single rock joint. Parametric studies are conducted to analyze the effects of the joint position and incident wave frequency. Results demonstrated that incident one-cycle sinusoidal pulse shifted into a wave that had two lobes after propagation in nonlinear rock within a certain distance. The wave then became a wave with two lobes in a different shape after normal transmission across a nonlinear single joint. The amplitude of the recorded waveform before the joint had an obvious inverse correlation with the frequency of the incident pulse. In addition, the amplitude of the transmitted waveform had a positive correlation with the incident pulse within a fixed distance from the joint. The magnitude of the transmission coefficient increased with the incident wave frequency. The conclusions drawn from this study provide a reference for the assessment of the stability of rock structures when they are subjected to dynamic disturbance.
1. Introduction
Rockmass is a material that consists of many defects and discontinuities, such as pores, cracks, fractures, joints, faults, and structural planes. These defects, which usually exhibit nonlinearity, directly affect the static and dynamic characteristic of rockmass [1, 2]. Due to these nonlinear discontinuities, the propagation rules of stress waves, such as seismic waves, blast waves, microseisms, and others, in a continuous medium are not applicable. Therefore, the study of wave transmissions across a single nonlinear joint is of great engineering significance for the assessment of the stability of rock structures after they undergo dynamic disturbance induced by excavations, vibrations, or explosions.
Over the last several decades, a great deal of studies have investigated wave propagation rules in rockmass. These studies have generally centered on two aspects: theoretical research and experimental investigations. Mccall and Guyer [3] performed a theoretical study that used a hysteric nonlinear dynamic modulus as input for a model of elastic wave propagation in rock. They also employed a Green function method to obtain its approximate analytical solution. Cheng [4] presented a numerical investigation of nonlinear wave propagation in sandstone by taking third-order elastic coefficients into consideration. Hokstad [5] extended the constitutive relationship to include third- and fourth-order nonlinear terms and a linear dispersive term, and he obtained numerical results using finite-difference modeling. Zheng et al. [6, 7] calculated similar nonlinear wave equations using high-order finite-difference operators with a simple flux-limiter algorithm. These studies deepened the understanding of wave propagation in nonlinear rock, but they did not consider the rock joint. To remedy this limitation, Zhao et al. [8, 9] reported the transmission process of a normally incident elastic P-wave across single-dry fractures using a nonlinear normal deformational behavior. The results showed that the transmission coefficient and the reflection coefficient of the linear model were special solutions of the nonlinear model, and these two parameters were directly related to fracture features. Li et al. [10] conducted an analytical study of underground explosion-induced ground motion according to the principle of conservation of momentum at wave fronts. They also subsequently derived the particle velocity and particle acceleration on the ground surface. Zou et al. [11] analyzed the propagation paths of waves using obliquely incident P- or S-wave propagations across a single joint. The superposition states in different areas adjacent to the fracture were obtained based on the incident wave propagation that was used.
In addition to the theoretical research, other scholars have examined the issue of transmission of waves across nonlinear fractures or joints using experimental investigations. Pyrak-Nolte et al. [12–14] performed a series of laboratory experiments on P- and S-wave transmissions across different natural rock fractures. They found that the effect of fractures on the spectral amplitudes of P- and S-wave pulses transmitted across fractures could be described well using a displacement discontinuity for compressional pulses under dry and saturated conditions and by using a combined displacement and velocity discontinuity for shear wave pulses under dry and saturated conditions. Ju et al. [15] performed an experimental and theoretical investigation of the properties of stress wave propagation in jointed rocks by means of the split Hopkinson pressure bar (SHPB) technique and fractal geometry method to explore the influence of a rough joint surface configuration on stress wave propagation. They concluded that a rough joint surface distinctly affected stress wave propagation and energy dissipation in jointed rocks. Huang et al. [16] performed experiments to investigate the relationship between the transmission ratio and normal stress, joint roughness, joint number, and frequency of incident waves when ultrasonic waves were transmitted across a rockmass with one joint and multiple parallel joints oriented normally. The experimental results supported similar conclusions based on analytical results drawn in [17–19] and [20–22]. The authors of [23, 24] conducted an experimental study on wave propagation across an artificial rock joint, which revealed the relationship between the transmission coefficient and the contact situation of the joint surface. The results showed that the contact area ratio and thickness both had an effect on P- and S-wave transmission.
However, the established numerical models and the experimental results from these studies failed to meet the criteria for actual situations. Therefore, in this study, a recurrence equation is deduced using a semianalytical and seminumerical method from the nonlinear wave equation by introducing a static Bandis−Barton (BB) model for a single rock joint. The aim of this study was to reveal the transmission laws of a one-dimensional P-wave across a single nonlinear joint in a rockmass. This work provides a reference for assessing the stability of rock structures when they undergo dynamic disturbance.
The study is organized as follows: Section 2 describes the nonlinear wave equation and static BB model. Section 3 presents theoretical formulas and a numerical solution to the nonlinear equation. Section 4 conducts parametric studies to explore the effects of joint position and incident wave frequency. Section 5 summarizes the study.
2. Nonlinear Wave Equation and Static BB Model
2.1. Nonlinear Wave Equation
In the past, there have been numerous scholars who have presented theories of nonlinear wave propagation in elastic mediums [25, 26]. Assuming the direction of the P-wave propagation is in the -direction, the one-dimensional stress-strain relationship in the classical theory of nonlinear acoustics is often written as follows:where is the stress; is the linear elastic constant in which is the density and is the constant speed of the P-wave traveling in the direction of the axis; is the small displacement of a particle parallel to the coordinate axis ; and the dimensionless parameters and are the second-order and third-order coefficients of nonlinearities, respectively. The nonlinear phenomena of a small amplitude wave can be well described using this model, but it is not suitable to reflect the frequency dispersion effect.
Extensive theoretical work has been performed to describe nonlinear propagation in a rockmass [27, 28]. Based on the nonlinear features of rock, Hokstad [5] proposed a nonlinear elastic model that contains a term related to the frequency dispersion effect and two other nonlinear terms related to the elastic constants. The model can be written as follows:where is the strain; the dimensionless parameters and are coefficients of nonlinearity related to the elastic constants; and is the dispersive term related to the frequency dispersion effect, in which is the coefficient of frequency dispersion.
Traditionally, the one-dimensional equation of motion for propagation in an infinite elastic medium in the absence of dissipation is expressed as follows:where is time.
Substituting Equation (2) in Equation (3) and assuming that the rock density remains unchanged before and after deformation, the following is obtained:wherewhere , and , , are the Lamé and Murnaghan constants, respectively.
In this study, the higher-order nonlinearity and frequency dispersion effect are not taken into consideration. Thus, the coefficient of frequency dispersion, , and the coefficients of third-order nonlinearity, , are taken to be zero. Then, Equation (4) degenerates to [4]:where . Interestingly, Johnson [29] derived the same equation from the Lagrangian function.
By assuming that the wave function at the wave source is , in which is the amplitude; is the wavenumber and ; and is the angular frequency, and the approximate analytical solution to Equation (6) can be obtained using the approximate analytical method and parametric variation method as follows:
In wave propagation experiments in rock [30], the nonlinear coefficients were calculated to be . By implanting Equation (7) with , , , and , the approximate analytical solution was obtained and plotted, as shown in Figure 1.

Compared to a linear wave, the nonlinear wave is shown in Figure 2. This figure shows the approximate analytical solution at , and it exhibits different characteristics. A high-order harmonic wave is generated in the nonlinear wave, and the amplitude of the particle displacement becomes increasingly stronger with an increase in the propagation distance. This is because, in the nonlinear wave, the phase velocities are variable due to the nonlinearity of rock and strains.

2.2. Static BB Model
To describe the deformational behavior of a rock joint under normal stress, a number of models [31–36], have been proposed based on experimental investigations. These include the power function model, the modified exponential function model, and the hyperbolic model. These models chiefly emphasize and describe the nonlinear characteristics of dry joints in rock materials well. The static BB model, proposed by Bandis et al. [31, 32], was established using a hyperbolic function that is written as follows:where is the joint closure; is the allowable joint closure that corresponds to the maximum of ; is the normal stress; and is the initial normal joint stiffness. The parameters and are directly affected by the joint roughness coefficient (JRC), mean aperture thickness, and joint compressive strength [31, 32].
Figure 3 shows a sketch of the static BB model for a rock joint. Compared to the linear model, the static BB model considers the joint stiffness, which results in a more accurate relationship between the normal stress and the joint closure of the rock. In addition, each parameter in this model has a definite physical meaning, making it more extensively theoretically applicable to actual projects in the field of rock mechanics and rock engineering.

3. Theoretical Calculation
3.1. Theoretical Formulas
Figure 4 gives the mechanical calculated model of transmission of a one-dimensional P-wave across a single nonlinear joint in a rock. It is assumed that the joint conforms to the static BB model at in the elastic half-space, and it has not been affected by the initial disturbance. The rock is assumed to be elastic, isotropic, and homogeneous. At an initial time (), a pulse in the normal direction is produced on the boundary of this elastic half-space (). After moving out from the wave source, the pulse, which is in the form of a nonlinear stress wave, travels in the nonlinear rock, and the propagation is governed by Equation (6). As the stress wave transmits across the joint, the seismic stress, , and particle displacement, , are assumed to be continuous and discontinuous, respectively. Specifically, the assumption can be given by the following two equations [8]:where the superscripts (+) and (−) denote the wave fields in the right- and left-hand of the joint, respectively.

By deducing the partial derivatives of Equation (7) with respect to and , the following two derivative expressions are obtained:where is the particle velocity. Thus,
Substituting the expression for , the relationship between the strain, , and the particle velocity, , becomes
According to [8, 17]; the method of characteristics has always been used to solve the problem of a one-dimensional wave transmission in a continuous medium. This method has been confirmed to be suitable for use with a one-dimensional P-wave normally propagating across a single joint in a rockmass. Therefore, the characteristic line equation of Equation (7) can be obtained. The method of characteristics derives , in which is a function of and . Therefore, the total differential of is as follows:
Comparing the coefficients of each differential term in Equation (14) with Equation (6), we find that
Substituting Equation (15) in Equation (14), the characteristic compatibility relation can be obtained as follows:
Therefore, the left-running characteristic line that is in the wave field after the joint can be expressed as follows:
In addition, the right- and left-running characteristic lines that are in the wave field before the joint can be expressed as follows:where the superscripts (+) and (−) are the wave fields after and before the joint, respectively, and is the particle velocity input to the boundary of the elastic half-space at time .
Substituting Equations (17) and (19) in Equation (18), the following is obtained:
This equation represents the relationship of particle velocities before and after the joint. Integrating both sides of this equation with respect to time, , the following is obtained:where and are the particle displacement before and after the joint, respectively. Substituting Equation (8) in Equation (10) and combining the result with Equation (21), the following is obtained:
Since the coefficients, and , are taken to be zero in this study, Equation (2) can be simplified as follows:
Substituting Equation (23) in Equation (22), the following is obtained:
By deducing the partial derivative of Equation (24) with respect to , the following is obtained:
Substituting Equation (13) in Equation (25), a complex differential equation with respect to can be derived. However, it is difficult to solve this equation directly, so the finite-difference method is employed to find a numerical solution for .
3.2. Numerical Solution
The period of the pulse that moves out from the wave source is denoted as , and it can be divided into equal time steps. The time after each time step is denoted as (), and the time increment, , for each time step is denoted as . Within each time interval, the left side of the equality sign in Equation (25), , can be written as follows:
Therefore, Equation (25) becomes
Equation (27) is a recurrence equation, and can be numerically calculated using it.
4. Parametric Studies
4.1. Effect of the Joint Position,
Parametric studies were conducted to gain a better understanding of the effects of a nonlinear joint in a rockmass on the characteristics of a P-wave normal transmission. The characteristics investigated included the distance of the joint to the wave source, the frequency of the incident wave, as well as the waveform of the incident wave. In the numerical calculation, the incident wave adopted was a one-cycle sinusoidal wave in the form of . The incident P-wave velocity () and the amplitude () of the incident wave were set to 5,000 and 0.005 mm, respectively. The coefficient of nonlinearity () was taken to be 7,000 based on the experimental result [30]. The density of the rock was 2,650 . The allowable joint closure () and the initial normal joint stiffness () were assumed to be 0.57 mm and 2 × 109 , respectively. The number of time increments () in a period () of an incident wave was selected to be 5 × 106, which is large enough to ensure a steady differential calculation. All the parameters related to the incident wave and the rock joint used for numerical calculation are summarized in Table 1.
Figure 5 illustrates the waveforms recorded in the two sides of the joint after propagation at the different distances of , , , , , and . It was deduced from the plot that the incident one-cycle sinusoidal pulse shifted into a wave that had two lobes after propagation in the nonlinear rock with a certain distance. Then, this wave transmitted across the nonlinear joint and became a wave with two lobes with a different shape. At a relatively low distance, the major lobe of the waveform before the joint occurred at the first trough and the minor lobe occurred at the second trough. With an increase in , the amplitude of the minor lobe became larger. After increased up to a certain value, the second trough gradually exceeded the first one and became the major lobe of the waveform. However, there was an inverse pattern of the waveform recorded after transmission across the joint. In addition, in all these cases, the absolute values of the maximum and the minimum particle displacement for the two waves recorded both before and after the joint increased with the distance. References [6, 7] studied one-dimensional nonlinear P-wave propagation in a solid without any joints or fractures, and they found that the particle displacement increased with the distance due to the formation of a second-order harmonic wave. In this study, the finding before the joint was in agreement with their conclusion. This phenomenon can be explained from the view of the energy conservation. Since the frequency dispersion effect was not considered in this study, the energy would not have dissipated during the entire process of wave propagation. All the energy of nonlinear wave propagation comes from the incident wave. When the amplitude of the harmonic wave gradually increased with distance, this led to an increase in potential energy, and the fundamental wave decreased to keep the total energy constant. Because the frequency of the fundamental wave did not change in this study, the particle velocity, which reflects the magnitude of the kinetic energy, then decreased with distance. Particle velocity attenuation has been proven to be a common phenomenon in the propagation of a nonlinear wave in a solid [6, 37, 38].

(a)

(b)

(c)

(d)

(e)

(f)
4.2. Effect of Frequency
By adopting the same parameters found in Table 1 and using frequencies of the incident wave ranging from 0.2 kHz to 4.0 kHz, the waveforms before and after the joint at for a one-cycle sinusoidal incident pulse were calculated, as shown in Figure 6. It can be observed from the figure that when the frequency of the incident pulse increases, the amplitude of the recorded waveform before the joint decreases, but the amplitude of the transmitted wave increases. Because of the nonlinearity of the rock and the joint, the peak and valley of the waves both before and after the joint at a certain frequency are unequal. In addition, it can be seen from this plot that an obvious distortion of the waveform begins to occur when the frequency is greater than 4.0 kHz. Also, the waveform before the joint has very little distortion and is approximate to a sinusoidal pulse when the frequency is at a lower value, for example, 0.2 kHz. From previous studies, the frequency of blasting vibration signals is commonly no more than 0.2 kHz [39, 40], and the microseismic precursor signals for rockbursts collected in the field [41–43] do not always exceed 1.0 kHz. Thus, the signals collected from safety monitoring for blasting excavation and rockburst in deep-buried tunnels and mines can reflect the damage and fracture of the inner structure more accurately in good integrity rock than it can with rocks with discontinuities, such as joints and fractures. In addition, the amplitude of the transmitted wave recorded after the joint at is much less than that of the wave recorded before the joint. Hence, sensors should be installed at more appropriate positions by considering the integrity of the rock when safety monitoring is conducted in situ.

The magnitude of the transmission coefficient in acoustics [44] is commonly defined as the ratio of the transmitted energy flux to the incident energy flux. But this definition is inapplicable when the frequency dispersion effect is not considered. Accordingly, [8, 9] defined the magnitude of the transmission coefficient as the ratio of the transmitted wave amplitude ( = peak value of to the incident wave amplitude (). According to this definition, the magnitude of the transmission coefficient () is defined in this study aswhere , , , and denote the maximum and minimum values of particle displacement before and after the joint.
The results of the magnitude of the transmission coefficient as a function of the incident wave frequency are illustrated in Figure 7. It can be seen from the plot that the magnitude of the transmission coefficient increases with the incident wave frequency. Zhao et al. [9] showed the relationship between the transmission coefficient and the amplitude of the incident wave. The finding herein gives the relationship between the transmission coefficient and the incident wave frequency, which is a supplement of previous studies.

5. Conclusion
To investigate the transmission laws of a one-dimensional P-wave across a single nonlinear joint in a rockmass, a recurrence equation was deduced using a semianalytical and seminumerical method obtained from the nonlinear wave equation by introducing the static BB model for a single rock joint. On the basis of the parametric studies conducted in this study, the following conclusions can be drawn:(1)The incident one-cycle sinusoidal pulse shifts into a wave that has two lobes after propagation in the nonlinear rock within a certain distance. The wave then becomes a wave with two lobes in a different shape after normal transmission across the nonlinear single joint. The amplitudes of the major lobe and the minor lobe of the waveforms recorded before and after the joint change with the distance of the nonlinear single joint from the wave source.(2)The remaining distance of the joint and the amplitude of the recorded waveform before the joint have an obvious inverse correlation with the frequency of the incident pulse. In addition, the amplitude of the transmitted waveform has a positive correlation with the incident pulse.(3)The magnitude of the transmission coefficient grows with the incident wave frequency.
The conclusions of this study, which revealed the propagation characteristics of a one-dimensional P-wave across a single nonlinear joint in a rockmass, are beneficial for applications in seismic modeling. However, the normal incidence of a P-wave is just a special case that cannot fully reflect the transmission features of a one-dimensional P-wave across joints. Future studies are needed to explore the results for a nonnormal transmission of a P-wave across a single nonlinear joint or a group of joints. It is also necessary to conduct experiments in the future to verify the conclusions of this study with SHPB.
Abbreviation
α: | Coefficient of frequency dispersion |
, , : | Coefficients of nonlinearity |
: | Strain |
, : | Second-order and third-order coefficient of nonlinearity |
, : | Lamé constants |
: | A function of particle velocity and strain |
: | Stress |
: | Normal stress |
: | Density |
: | Wave angular frequency |
: | Time increment |
: | Amplitude of incident wave |
: | Speed of P-wave traveling in the direction of the axis |
: | Joint closure |
: | Allowable joint closure |
: | Frequency of incident wave |
: | Wavenumber |
: | Initial normal joint stiffness |
, , : | Murnaghan constants |
: | Number of time increment in a period of incident wave |
: | Particle velocity input to the boundary of the elastic half-space at time |
: | Time |
: | Time after each time step |
: | Coordinate |
: | Distance of joint to wave source |
: | Displacement of particle |
, : | Particle displacement before and after the joint |
: | Particle velocity |
, : | Particle velocity before and after the joint |
: | Linear elastic constant |
: | Transmission coefficient |
: | Period of incident wave |
, : | Peak values of particle velocity before and after the joint |
, , , : | Maximum and minimum value of particle displacement before and after the joint |
BB: | Bandis–Barton |
SHPB: | Split Hopkinson pressure bar. |
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (grant no. 51674061), the State Key Research Development Program of China (grant no. 2018YFC0604601), and the program for innovative talents in Liaoning University (grant no. LR2016024).