## 4 Study of the Bearing Characteristics of NCPDRE
### Numerical Simulation Modelling  
#### 4.1.1 Determination of the Parameters of the Goaf Collapse Zone  
The process of goaf gangue compaction needs to be considered when using numerical simulation, for which the double-yield model is widely used [32; 33; 34; 35]. If the model is applied to the numerical simulation of the 2109 working face of Qingwa Coal Mine, the parameters of the model need to be determined according to the actual engineering situation. The goaf material should conform to the stress-strain relationship of the broken rock mass derived by [PERSON] [36; 37; 38; 39]:  
\[\sigma_{z}=\frac{E_{0}\varepsilon_{z}}{1-\frac{\varepsilon_{z}}{\varepsilon_{ max}}} \tag{10}\]  
where \(\sigma_{z}\) is the vertical stress of the rock body in the extraction zone, MPa; \(E_{0}\) is the initial modulus of elasticity of the rock body, MPa; \(\varepsilon_{z}\) is the body strain under the action of vertical stress; and \(\varepsilon_{max}\) is the maximum body strain. \(E_{0}\) and \(\varepsilon_{max}\) can be calculated using the following equation:  
\[E_{0}=\frac{10.39\sigma_{c}^{1.042}}{b^{7.7}} \tag{11}\]  
\[\varepsilon_{max}=\frac{b-1}{b} \tag{12}\]  
where \(\sigma_{c}\) is the uniaxial compressive strength of the rock mass, MPa; \(b\) is the coefficient of rock mass fragmentation, calculated as follows:  
\[b=1+\frac{c_{1}m+c_{2}}{100} \tag{13}\]  
where \(m\) is the mining height of the coal seam, m; \(c_{1}\) and \(c_{2}\) are the lithology coefficients of the roof [40], which can be calculated from Table 1. The calculated value of b was 1.37. Figure 4: Coal pillar width and its bearing stress curve. Inputting the parameters into Equation (10), the calculation yields the following:  
\[\sigma_{z}=\frac{26.63\varepsilon_{z}}{1-3.7\varepsilon_{z}} \tag{14}\]  
According to the above equation, the stress and strain correspondence in the process of rock mass collapse and compaction can be calculated, as shown in Table 2.  
In the FLAC3D6.0 numerical simulation software, the cap pressure \(p_{c}\) and plastic volume strain \(\varepsilon^{p}\) of the double-yield model need to be defined by the user, so that the stress-strain curves during the collapse and compaction of the rock mass can be simulated in accordance with the Salamon formula. In this formula, \(p_{c}\), \(\sigma_{Z}\), \(\varepsilon^{p}\), and \(\varepsilon_{Z}\) are determined using the correlation below, and for any set of \(\sigma_{Z}\) and \(\varepsilon_{Z}\), the cap pressure \(p_{c}\) and plastic volume strain \(\varepsilon^{P}\) can be calculated using the following formula [35]:  
\[p_{c}=\frac{1+2\eta}{3}\sigma_{z} \tag{15}\]  
\[\varepsilon^{p}=\frac{R}{R+1}\varepsilon_{z} \tag{16}\]  
where \(p_{c}\) is the cap pressure, MPa; \(\eta\) is the lateral pressure coefficient; \(\varepsilon^{p}\) is a plastic volumetric strain; \(R\) is the proportionality factor. According to the inversion data in Table 3, a FLAC3D numerical simulation model of a cube with an edge length of 1 m was established, the velocity boundary was fixed around and at the bottom, and then a fixed velocity of 1 \(\times\) 10\({}^{-5}\) m/step was applied to the top for loading. The vertical stress of the Z = 0 plane was obtained by writing the FISH function, and then the model stress and its displacement at each time step were recorded using the history command of FLAC3D. The vertical stress-strain curve of the numerical model was then drawn. Comparing the numerically inverted curve with the theoretically calculated curve in Figure 5, it can be seen that the numerical simulation curve obtained using the cap pressure method is close to the theoretical curve. The parameters of the goat double-yield model are shown in \[E_{0}=\frac{10.39\sigma_{c}^{1.042}}{b^{7.7}} \tag{11}\]0