Analysis and Comparison of Three Bending Tests on Phosphogypsum-Based Material According to Peridynamic Theory

Abstract

Phosphogypsum-based materials have gained much attention in the field of road infrastructure from the economic and sustainable perspectives. The Three-point bending test, the Four-point bending test and the Semi-circular bending test are three typical test methods applied for fracture energy measurement. However, the optimal test method for fracture energy evaluation has not been determined for phosphogypsum-based materials. To contribute to the gap, this study aims to analyze and compare the three test methods for fracture energy evaluation of phosphogypsum materials based on the peridynamic theory. For this purpose, the load–displacement, vertical displacement–Crack Mouth Opening Displacement (CMOD) and fracture energy of the phosphogypsum-based materials were measured and calculated from the three test methods. The simulated load–displacement and vertical displacement–CMOD by PD numerical models, with different fracture energy as inputs, were compared to the corresponding tested values according to simulation error results. The results showed that the Four-point bending test led to minimized errors lower than 0.189 and indicators lower than 0.124, demonstrating the most optimal test method for the fracture energy measurement of phosphogypsum-based material. The results of this study can provide new methodological references for the selection of material fracture energy measurement tests.

1. Introduction

Phosphogypsum, as a raw material of base layers, is one of the primary concerns for road technology from the economic and sustainable perspectives [1,2]. Most studies have focused on the mix design and mechanical characterization of phosphogypsum. Some scholars verified that phosphogypsum mixed with crushed rocks, clay, lime, fly ash and curing agents in certain proportions could be a substitute for traditional semi-rigid bases [3,4,5]. However, limited analysis has been carried out on the crack resistance of phos-phogypsum-based base materials, while extensive investigations have been conducted on the performance of traditional semi-rigid base materials [6,7]. To promote the widespread application of phosphogypsum as a semi-rigid base for road construction, it is imperative to concisely characterize the crack resistance of phosphogypsum-based material.

The fracture energy is a critical parameter of materials in evaluating the ability to resist damage-induced cracking. Moreover, Mode I cracks are the most common form of base materials. Thus, Mode I fracture energy represents the focus of fracture research [8]. Currently, three experimental methods are conducted for obtaining Mode I fracture energy, including direct tension testing [9], indirect tension testing [10], and bending testing [11,12]. Among the three test methods, an appropriate cracking test method for accurately measuring the fracture energy of phosphogypsum has not yet been determined. In this study, identifying the most suitable method for acquiring fracture energy will not only provide a reliable experimental methodology, but also propose a guideline for further research on the fracture failure mechanisms of phosphogypsum and the rational design of structures.

Theoretically, fracture energy derived from crack initiation tests on standardized cross-sections should exhibit uniformity. However, the results obtained from different tests vary in practice [13]. Such discrepancies may stem from dissimilarities in energy dissipation, deviations in crack propagation from the principal tensile stress direction, and other factors [14,15]. To address these root causes, Zhang et al. [16] investigated kinetic energy dissipation, crack propagation trajectories, and deformation energy beyond the fracture surface during loading. These critical sources of error were thoroughly scrutinized, serving as the foundation for selecting experimental methods. Nonetheless, this research paradigm, grounded in extensive theoretical groundwork and experimental methodologies, presented formidable challenges and substantial cost.

To tackle these deficiencies, another approach to evaluate optimal fracture energy was proposed in this study, grounded in numerical simulation. Accurate fracture energy values, as inputs of damaged numerical models, can significantly enhance the accuracy of numerical simulation outcomes, rendering them more faithful to practical damage scenarios [17]. If fracture energy data measured from a specific test type can effectively accommodate all fracture test methods within numerical models, it means the test method is the most optimal method for fracture energy evaluation. This approach, underpinned by insights gleaned from numerical model results fueled by fracture energy, explores the veracity of fracture energy data with avoiding the test errors.

Regarding traditional damaged numerical models, numerical outcomes are influenced not only by fracture energy but also by additional factors, including tensile strength, criteria for crack propagation deflection, and handling of singular regions [18,19]. This fact indicates that fracture energy is insufficient as the sole input parameter for traditional numerical models. In light of this challenge, this study adopts peridynamic (PD) as the analytical framework, as the fracture energy is the only input characterizing the entire material damage process [20]. This theory has developed from its initial applicability to models for linear elastic–brittle materials to diverse theoretical models suitable for non-linear elastic–brittle, elastic–plastic, viscoelastic, and viscoplastic materials [21,22,23]. Moreover, it has demonstrated successful application in the fracture testing of engineering materials such as asphalt mixtures and cement concrete with excellent simulation outcomes [24,25].

Based on the above background, few studies focused on selecting a suitable fracture energy testing method for phosphogypsum-based materials among different test modes. Thus, this study addresses this gap, as well as delving into the difference in fracture energy between these tests, based on PD simulation. In this study, suitable test methodologies for phosphogypsum materials across various fracture energy testing scenarios are delineated, employing numerical simulations and aiming to provide an optimal fracture energy testing method. The fracture energy of phosphogypsum-based specimens will be measured by the Three-point bending test [26], the Four-point bending test [27], and the Semi-circular bending test [28] and will be the input for peridynamic numerical simulations. The optimal fracture energy testing method will be determined by analyzing the simulation error between tested and simulated results.

2. PD Theory

Peridynamic is a method based on the concept of non-local interactions, and describes the mechanical behavior of materials by solving spatial integral equations. To provide theoretical and technical support for determining the optimal test method in the following sections, this section briefly introduces the theory of PD and its numerical implementation process. Firstly, the form of the motion equations is introduced. Secondly, the specific form of the force density function in the motion equations is provided, which are the fundamental equations of the PD theory. Finally, the methods for solving the equations in the numerical computation process are explained.

2.1. Kinematic Equation

The two-parameter bond-based PD model was selected due to its capacity to address the limitations associated with fixed Poisson’s ratios in traditional bond-based models. Additionally, it offers notable computational advantages compared to state-based models [29].

As depicted in Figure 1a, at any time instance t, each particle x_k within a spatial domain R encompasses a family domain with a radius of δ. Within this domain, other material points are interconnected to the central point in a bond-like fashion, facilitating the transmission of forces. The specific formulation of motion equations for the two-parameter model is derived from the principle of virtual work, shown in Equation (1) [30]:

$$\begin{gathered} \rho \frac{{\partial }^{2}u\left(x_k,t\right)}{\partial t^2}={\displaystyle {\int }_H{f}_c\left(u\left(x_j,t\right)-u\left(x_k,t\right),x_j-x_k,{\theta }_k,{\theta }_j\right)}d{V}_{x_j}+b\left(x_k,t\right), \\ I\frac{{\partial }^2\theta \left(x_k,t\right)}{\partial t^2}={\displaystyle {\int }_H{f}_{\kappa }\left(u\left(x_j,t\right)-u\left(x_k,t\right),x_j-x_k,{\theta }_k,{\theta }_j\right)}d{V}_{x_j}+M\left(x_k,t\right) \end{gathered}$$

(1)

where ρ represents material density, I denotes moment of inertia; u signifies the displacement of material points, and θ represents the rotation angle of material points; b stands for body force density, and M represents external torque; f_c and f_κ, respectively, represent normal bond force and tangential bond force.

2.2. Force Function

In this study, the phosphogypsum-based materials used for road base layer were assumed as ideal elastic–brittle materials [31]. Subsequently, the ensuing discussion introduces the theoretical framework of the elastic–brittle model within the context of the two-parameter bond-based PD.

In bond-based PD models, the constitutive force function solely delineates the radial displacement and radial force interactions among material points, thereby inadequately capturing the comprehensive motion state of these points. In contrast, the two-parameter model discerns between shape alterations and volume variations by concurrently incorporating axial forces resulting from axial deformation and bending moments arising from rotational deformation between material points. This approach effectively addresses the constraint of fixed Poisson’s ratio [29]. The definitions of the normal bond force and tangential bond force are as follows:

$$\begin{gathered} {\overrightarrow{f}}_c〈u^{\prime }-u,X^{\prime }-X〉=c\overrightarrow{s}〈u^{\prime }-u,X^{\prime }-X〉, \\ {\overrightarrow{f}}_{\kappa }〈u^{\prime }-u,X^{\prime }-X〉=\kappa \overrightarrow{\beta }〈u^{\prime }-u,X^{\prime }-X〉 \end{gathered}$$

(2)

where ${\overrightarrow{f}}_c$ and ${\overrightarrow{f}}_{\kappa }$ represent normal bond force and tangential bond force, respectively, and $c$ and $\kappa$ are the PD parameters in axial direction and tangential direction, respectively; $\overrightarrow{s}$ represents the bond stretch, and $\overrightarrow{\beta }$ represents the change in the bond rotation.

Following the principle of equal strain energy density, the relationship expression between PD parameters and Young’s modulus E, as well as Poisson’s ratio ν, can be derived. The specific expressions are detailed in

Table 1. Parameters of the PD Model.

Model Parameters	Three-Dimensional	Plane Stress	Plane Strain
c	$\frac{6E}{\pi {\delta }^3\left(1-2v\right)\left(1+v\right)}$	$\frac{6E}{\pi {\delta }^{3}t\left(1-v\right)}$	$\frac{6E}{\pi {\delta }^{3}t\left(1+v\right)\left(1-2v\right)}$
k	$\frac{E\left(1-4v\right)}{4\pi {\delta }^2\left(1-2v\right)\left(1+v\right)}$	$\frac{E\left(1-3v\right)}{6\pi \delta t\left(1-v^2\right)}$	$\frac{E\left(1-4v\right)}{6\pi \delta t\left(1-2v\right)\left(1+v\right)}$

, with a comprehensive derivation process provided in ref. [32].

2.3. Failure Criterion

Combining Silling and Askari’s derivation of fracture criteria [33], the fracture energy of the material is expressed as follows:

$$G^{PD}=2h{\displaystyle {\int }_0^{\delta }{\displaystyle {\int }_z^{\delta }{\displaystyle {\int }_0^{\arccos \frac{z}{\xi }}\left(\frac{1}{2}c{s}_0^2\xi +\frac{1}{2}\kappa {\beta }_0^2\xi \right)}}}\xi d\phi d\xi dz$$

(3)

where S₀ represents the critical bond stretch, and β₀ represents the critical bond rotation, as illustrated in Figure 1b.

For a Mode I crack, the material undergoes tensile failure, and the energy at this juncture is solely attributed to the energy derived from normal bond deformation. Consequently, the critical bond stretch for tensile failure is determined as follows:

$$s_0=\sqrt{\frac{4G_I}{hc{\delta }^4}}$$

(4)

where $G_{\mathrm{I}}$ represents the experimentally measured type I fracture energy of the material, and h represents the material thickness.

In summary, three mandatory mechanical parameters, the Young’s modulus, Poisson’s ratio, and fracture energy, are required for PD numerical model simulations. The acquisition of these parameters will be elaborated upon in Section 3.

2.4. Numerical Method

To obtain the solution of the PD motion equations, it is necessary to discretize the spatial domain occupied by the material into a finite number of material points during the numerical computation process. Thus, the spatial integral form of the motion equations is transformed into a discrete sum form. For the two-parameter model, according to Equation (1), its motion equations can be discretized as [34]:

$$\begin{gathered} \rho \frac{{\partial }^{2}u\left(x_i^n,t\right)}{\partial t^2}={\displaystyle \sum _{j=1}^p{f}_c\left(u\left(x_j^n,t\right)-u\left(x_i^n,t\right),x_j^n-x_i^n,{\theta }_i^n,{\theta }_j^n\right)}{V}_{x_j}+b\left(x_i^n,t\right), \\ I\frac{{\partial }^2\theta \left(x_i^n,t\right)}{\partial t^2}={\displaystyle \sum _{j=1}^p{f}_{\kappa }\left(u\left(x_j^n,t\right)-u\left(x_i^n,t\right),x_j^n-x_i^n,{\theta }_i^n,{\theta }_j^n\right)}{V}_{x_j}+M\left(x_i^n,t\right) \end{gathered}$$

(5)

The time iteration in this study adopts an explicit finite difference method, where the differentiation of displacement with respect to time is converted into a finite difference formula. The velocity and acceleration at any time step are represented by the following equations:

$$\begin{gathered} \frac{{\partial }^2{u}_i^n}{\partial t^2}=\frac{u_i^{n+1}-2u_i^n+u_i^{n-1}}{\Delta t^2}, \\ \frac{\partial u_i^n}{\partial t}=\frac{u_i^{n+1}-u_i^{n-1}}{2\Delta t} \end{gathered}$$

(6)

Based on Equation (6), after obtaining the acceleration of a material point at the nth time step, the velocity and displacement at the next time step can be computed.

$$\begin{gathered} \frac{\partial u_{x_i}^{n+1}}{\partial t}=\frac{{\partial }^2{u}_{x_i}^n}{\partial t^2}\Delta t+\frac{\partial u_{x_i}^n}{\partial t}, \\ {u}_{x_i}^{n+1}=\frac{\partial u_{x_i}^n}{\partial t}\Delta t+u_{x_i}^n \end{gathered}$$

(7)

where $\Delta t$ represents the time step length during the numerical computation process.

The selection of time step length is crucial in explicit time differencing. To obtain convergent computational results, the time step must satisfy certain stability conditions. The choice in this study is based on the work of Silling and Askari [30].

2.5. Simulation Result Processing Method

In order to simulate the bending experiment of phosphogypsum-based material in PD theory, virtual material layers are set at the top and bottom of the sample. During the loading process, the virtual material layers will bear the reaction forces applied by the actual materials, and the relationship between the axial force P and the vertical component $F_y^i$ of the PD force density is expressed as follows:

$$P=\left|{\displaystyle {\sum }_{i=1}^N{F}_y^i{V}_i}\right|$$

(8)

where N represents the total number of material points in the virtual layer.

3. Methodology

In order to ascertain the fracture energy of phosphogypsum and furnish the material parameters for the numerical model delineated in Section 4, this section introduces the preparation of phosphogypsum specimens, the mechanical test process, and establishment of numerical models.

3.1. Materials and Sample Preparation

In this study, the phosphogypsum-based material applied for the road base layer is the pre-treated material with optimized water content, as determined and studied within the framework of the project titled “Research on Environmentally Friendly and Durable Phosphogypsum-Based Assembled Road Pavement Material and Structure”.

As presented in Figure 2, the preparation process of this material was divided into five steps: First, the phosphogypsum powder was completely mixed with water for about 1.5 min until the lumps disappear. Second, the mixed phosphogypsum material was poured into beam and cylindrical molds and pressed for 15 min into a mold. Third, using a demolding machine, beam specimens with dimensions of 50 mm × 50 mm × 200 mm and cylindrical specimens with dimensions of Φ150 mm × 100 mm were formed. Fourth, the demolded specimens were cut into beam-shaped and semi-circular specimens with dimensions of 50 mm × 50 mm × 100 mm and Φ150 mm × 40 mm, respectively. Then, pre-cut cracks of 10 mm and 15 mm, respectively, were made at the midspan of the beam-shaped and semi-circular specimens for better control of the crack propagation path. Lastly, all specimens were cured at 20 °C for seven days before testing.

3.2. Tests

The bending test is a commonly utilized and effective method for evaluating the fracture performance of brittle materials [35]. In this study, Three-point bending, Four-point bending, and Semi-circular bending tests were performed for characterizing the fracture energy of phosphogypsum-based materials. Additionally, the uniaxial compression tests were conducted to measure Young’s modulus and Poisson’s ratio of the phosphogypsum-based materials. Each type of experiment is conducted with six parallel tests, and the three groups with the least variability are selected as the experimental results. The instruments, test parameters, and procedures employed in the compression and bending tests are outlined as follows.

3.2.1. Bending Test

The bending tests were performed by a servo-pneumatic Mechanical Testing System (MTS) with a capacity of 50 KN, produced by Cooper Technology (Borken, Germany). The MTS exerted a controlled linear axial load with a steady rate of 0.2 mm/min to ensure the quasi-static loading.

The illustration of the bending tests is shown in Figure 3. The beam specimen dimensions refer to the recommendations for fracture energy testing in the “Test Methods of Materials Stabilized with Inorganic Binders for Highway Engineering” (JTG 3441-2024) [36]. The semi-circular specimen dimensions follow the International Society for Rock Mechanics (ISRM) guidelines for Semi-Circular Bend (SCB) testing [37]. During the loading process, vertical displacement and vertical load were recorded in the MTS, CMOD was measured by a YYJ-12/10 clip gauge produced by the Central Iron & Steel Research Institute, Beijing, which was used to form the load–displacement curve and the displacement–CMOD curve. The fracture energy is the area below the load–displacement curve until the specimen is broken, and describes the energy released during cracking.

$$G_f=\frac{W}{A_{lig}}$$

(9)

where W is the fracture’s work and Alig is the ligament area, defined as:

$$A_{lig}=\left(r-a\right)xh$$

(10)

where r is the specimen radius or height, in m, a is the notch length, in m, and h specimen thickness, in m.

3.2.2. Uniaxial Compression Test

The uniaxial compression test was performed using the Wdw-300 Series Computerized Electronic Universal Testing Machine manufactured by Kexin (Shenzhen, China). The machine is capable of exerting a linear axial load. Four strain gauges were affixed vertically and horizontally on opposite sides of the specimen. According to the “Test Methods of Materials Stabilized with Inorganic Binders for Highway Engineering” (JTG 3441-2024), the UC modulus test was conducted in a displacement-controlled loading mode at a constant rate of 1 mm/min. The decision to use prismatic specimens is due to challenges in measuring Poisson’s ratio with cylindrical specimens.

The illustration of the UC test is shown in Figure 4. The UC test determines the Young’s modulus and Poisson’s ratio uniaxially, as it measures the stress σy(t) and the strain εx(t) and εy(t), which is expressed as follows:

$$E=\frac{{\sigma }_y\left(t\right)}{{\epsilon }_y\left(t\right)},\nu =-\frac{{\epsilon }_x\left(t\right)}{{\epsilon }_y\left(t\right)}$$

(11)

In the UC test, the value of axial stress, axial strain and transverse strain were recorded, which constitute the stress–strain curves and bi-directional strain curves. The value of stiffness from the linear segment of the stress–strain curve determined the elastic modulus of the UC test specimen at a certain test condition. The peak value of axial strain and transverse strain were chosen for calculating the Poisson’s ratio.

3.3. Numerical Model Establishment

To fulfill the aim of identifying the optimal fracture energy testing approach in this study, this section established distinct PD numerical models for the Three-point bending test, the Four-point bending test, and the Semi-circular bending test. The fracture energy data garnered from varied bending tests served as the input parameters in their respective numerical models. Subsequently, the simulated outcomes were compared to the actual test results for analyses.

The virtual experiments conducted in this study were addressed as plane stress problems. Two rectangular plates measuring 50 mm × 200 mm and one semi-circular plate with a diameter of Φ150 mm were configured. They were discretized into particles, forming a simple cubic lattice with ∆x = 0.5 mm and δ = 3.015∆x. Cracks of 10 mm and 15 mm were, respectively, set at the midspan of the rectangular plates and the semi-circular plate by truncating the PD bonds. The schematic diagram of the virtual model is shown in Figure 5.

To achieve a more precise simulation of the phosphogypsum-based material in this study, all material property parameters utilized in the model were derived from the experimental tests conducted in Section 3.2. The loading and support areas consisted of 6 × 6 particles. For boundary conditions, a constant loading rate of 0.2 mm/min was applied downward at the upper boundary. Through a combined optimization test ensuring integral stability and computational efficiency, the time step was established as ∆t = 1.33×10⁻⁶ s in this study. Considering the fracture energy results from the three types of bending tests, the critical bond stretch was computed by Equation (4). Subsequently, three critical bond stretch values were sequentially input into every numerical model to explore the simulation outcomes of each type of virtual test under varied fracture energy inputs. In this study, the load–displacement curve and the displacement–CMOD curve were used for test scheme comparison.

4. Results and Discussion

This section constitutes the core of determining the optimal testing method. Firstly, the experimental results were presented and analyzed. Subsequently, by inputting the material parameters obtained in Section 3 into the PD numerical model, the simulation results for various virtual bending tests under different fracture energy inputs were obtained. Finally, three error metrics were selected to compare the test results and the simulation results. The optimal test was selected based on a comprehensive analysis of simulation error analyses.

4.1. Test Results

4.1.1. Load–Displacement Curves and Fracture Energy

Figure 6a–c show the three sets of load–displacement curves for the phosphogypsum-based specimens under three types of bending tests. The pre-peak curves of the Four-point bending test demonstrated a linear ascending trend. In contrast, the pre-peak curves of the Three-point bending test and the Semi-circular bending test showed an exponential uptrend. Specifically, the slope gradually increased as the curve progresses and approaches the peak point, indicating a gradual acceleration of the upward trend. For the three tests, the curve experienced a sudden decline after the peak point, which signifies the occurrence of fracture. This phenomenon verifies the brittle nature of the phosphogypsum materials and the assumption of the ideal elastic–brittle materials. With the increase in load, the stress concentration at the specimen crack leads to the appearance of small cracks, which rapidly propagate into through cracks. The sudden decline after the peak point indicated that this process was brief. The rapid transition from microcrack generation to through crack formation in brittle materials emphasizes the criticality of crack-resistant design for engineering brittle materials. The averaged results of the calculated and analyzed test outcomes, encompassing fracture energy, failure load, and failure opening displacement, are presented in

Table 2. Bending tests results.

Types	Peak Load/KN	Failure Displacement/mm	Failure CMOD/mm	Fracture Energy (J/m2)
Three-point bending	1.009	0.4512	0.014	91.58
Four-point bending	1.286	0.3377	0.017	109.33
Semi-circular bending	1.691	0.5027	0.021	130.25

. The fracture energy data revealed significant disparities of approximately 10–12% in the results across different tests, surpassing the variability within each test type. Among the three tests, the Semi-circular bending test resulted in highest value in fracture energy, followed by Four-point bending and Three-point bending. This fact is connected with the stress distribution and crack propagation paths, which are influenced by the shape and size of the specimen, as well as the loading method.

4.1.2. Displacement–CMOD Curves

Figure 7 summaries the displacement–CMOD curves obtained from the three types of bending tests. The displacement–CMOD curve of the Four-point bending test showed a linear rise at the beginning and kept stable after 0.017 mm. However, the pre-peak segments of the curves for the Three-point bending test and the Semi-circular bending test grew with a gradually decreasing slope. Moreover, the Semi-circular bending test led to deepest vertical displacement compared to the Four-point bending test and the Three-point bending test. These results are consistent with the results of load–displacement curves.

4.1.3. Young’s Modulus and Poisson’s Ratio

Figure 8 depicts the load–displacement curve and bi-directional strain curve derived from the uniaxial compression test. Young’s modulus and Poisson’s ratio are individually calculated by analyzing the linear segments with curves. As a result, 3.3 GPa and 0.29 were calculated and used for the Young’s modulus and Possion’s ratio of phosphogypsum-based materials, respectively. The stress–strain curve exhibited a brief yielding phase, followed by a rapid decrease upon surpassing the elastic limit. This is consistent with the characteristics of elastic–brittle materials.

4.2. Comparison between the Three Tests Based on PD Simulation

4.2.1. Numerical Results

Figure 9a–c illustrate the comparison between the simulated load–displacement curves under different fracture energy index inputs and the test results for various bending tests. The selected interval for the simulation results correspond to the range where the load demonstrates a stable increasing trend. Upon analyzing the simulated curves, it is shown that variations in fracture energy influenced the endpoint of the pre-peak curve, while the curve’s trend was hardly affected by the changes in fracture energy. This result suggests that the characteristics of the pre-peak curve, particularly the linear segment, are primarily dependent on Young’s modulus and Poisson’s ratio for the same type of fracture test. Additionally, the Four-point bending test resulted in similar load–displacement curves between the simulated and test results among the three tests, followed by the Three-point bending test. However, the simulated results of the Semi-circular bending test were significantly different from the test results compared to that of the other tests, which indicates that the evolving traits of the Semi-circular bending test curve is partially characterized by the numerical model.

Figure 9d–f show the comparison between the simulated results and the actual experimental results of the vertical displacement–CMOD curve. The curves under different fracture energy input values exhibited a similar changing trend to the load–displacement curves, entering the failure stage successively after undergoing approximately the same trend. Among them, the simulated results of the Semi-circular bending model captured the trend of curvature change during the experimental process, corresponding to the simulated results of the load–displacement curve.

By comparing experimental and simulation results, it is observed that there are discrepancies in the rising stage between the experimental and simulated results for both the Three-point bending and Semi-circular bending tests according to the load–displacement curves or the displacement–CMOD curves. The main reason for the error is believed to be the difference between the selected constitutive model and the actual material characteristics. The ductility characteristics of the phosphogypsum-based material used in this study can be reflected in the descending section of Figure 8a. However, the bond-based brittle peridynamic constitutive model neglects the ductility of cementitious materials. Additionally, the consistency between the test results and the simulation results of the Four-point bending test suggests that the loading method of symmetrical double pressure heads can, to some extent, take the ductility characteristics of the material into consideration.

4.2.2. Quantitative Comparison between the Three Test Methods Based on Simulation Error

To quantitatively evaluate the fitting degree between the simulated curve and the actual experimental curve, this study introduces three parameters: Peak error, Difference error, and Absolute error [38]. Peak error reflects the accuracy of the value at the peak, that is, the accuracy of the value at the failure point. Difference error represents the error of the peak of the curve on the abscissa. Absolute error reflects the accuracy of the simulation results throughout the process. The results of the above three types of errors can be seen in Figure 10. The lower values of Peak error, Difference error and Absolute error indicate a better fitting degree between the simulated curve and the actual experimental curve, and vice versa.

Figure 11 displays the heatmap of the three parameters between the simulated results and the test results. The depth of color corresponds to the magnitude of the error value. In the coordinates, T, F, and S represent the Three-point, Four-point, and Semi-circular bending tests, respectively. For example, the coordinate (G(T), Accuracy(F)) means the error value between the simulated and experimental results of the Four-point bending test when using the fracture energy obtained from the Three-point bending test as the input into the numerical model. Most of the values on the main diagonal are the smallest in their respective rows, and the value at the center point is generally the smallest on the main diagonal. These results indicate that the fracture energy measured in each test ensures the accuracy when simulating the corresponding experiment, and the Four-point bending simulation has the highest precision. The two ends of the sub-diagonal have the largest values, implying that the fracture energy results obtained from the Three-point bending test and the Semi-circular bending test are not suitable for each other’s numerical models for the materials. From the perspective of row distribution, the second row in each graph has relatively smaller values compared to the other rows, suggesting that the fracture energy measured from the Four-point bending test can achieve relatively high simulation accuracy for all test types. These facts preliminarily indicate that using the fracture energy obtained from the Four-point bending test as the input for the numerical models of different types of tests generally leads to relatively higher simulation accuracy. To make a more informed selection of the optimal experimental plan, this study supplemented three additional judgment criteria.

The rationale for selecting the optimal experimental plan for obtaining fracture energy is as follows: The fracture energy data should not only serve the numerical model of the bending test, from which the data are obtained effectively, but also be applicable to other bending test numerical models. Additionally, when this fracture energy data are used as input, the difference in simulation accuracy between the numerical models of different types of tests should not be too large. Otherwise, it can be inferred that the fracture energy data are not capable of serving numerical models of different types of bending tests. Based on this premise, this study relies on three indicators which were used for characterizing the simulation error to select the best experimental plan.

Indicator 1 represents the simulation error between the test results and the numerical model results, where the input fracture energy is obtained from the corresponding test. Indicator 2 indicates the average of two simulation errors between the test results and the numerical model results, where the input fracture energy is obtained from the other two tests. Indicator 3 represents the discrepancy of three simulation errors between the test results and the numerical model results, where the input fracture energy is obtained from the three tests.

Figure 12 summarizes the ranking sequences for all indicators. It is obviously observed that the fracture energy of the Four-point bending test generally resulted in the smallest values for indicators 1, 2, 3 of the three errors among the three tests. Additionally, the fracture energy of the Three-point bending test resulted in smallest value for Indicator 3 of Peak error. The above facts demonstrate that the Four-point bending test is the most suitable test for fracture energy characterization.

5. Conclusions

This study aims to analyze and compare the three bending tests based on the peridynamic theory and the numerical simulation method. According to the attained results, the conclusions are drawn as follows:

• Regarding the phosphogypsum-based materials, there were significant differences in the load–displacement curves, displacement–CMOD curves and fracture energy between different bending tests, resulting in a variation of approximately 10–12% in the obtained material fracture energy.
• The simulated load–displacement and displacement–CMOD values were consistent with the test results based on the elastic–brittle PD model. And the results of load–displacement and fracture energy of the three tests were in good agreement with that of the displacement–CMOD of the three tests.
• The three error and indicator analyses indicate that the Four-point bending test is the best test method for fracture energy evaluation, followed by the Three-point bending test and the Semi-circular bending test.

Given that road base materials are commonly assumed to exhibit elastic–brittle behavior, this study adopts a peridynamic elastic–brittle constitutive model. It is worth noting that different conclusions may arise when employing alternative constitutive models. Moreover, peridynamic theory has also developed a criterion that controls the entire process of material damage based on tensile strength. Future research based on this criterion is recommended to investigate the optimal methods for tensile strength evaluation.