Carbon dioxide phase change blasting numerical simulation method based on FDEM

Abstract

The invention provides a carbon dioxide phase change blasting numerical simulation method based on an FDEM. The carbon dioxide phase change blasting numerical simulation method comprises the following steps that a carbon dioxide adiabatic expansion p-V relation is determined according to a Span-Wagner fluid model; establishing a finite-discrete element (FDEM) numerical model, and setting a non-reflection boundary; determining mesoscopic parameters in the FDEM model according to the physical and mechanical parameters of the rock mass; and applying a carbon dioxide phase change blasting load, and solving a finite-discrete element (FDEM) model. According to the method, the Span-Wagner fluid model recommended by the American national standard is adopted to replace a traditional ideal gas model to calculate the adiabatic p-V relation of the carbon dioxide, the description of the carbon dioxide pressure is closer to the reality, the value of the FDEM mesoscopic parameter is associated with the physical and mechanical indexes of the rock mass, a theoretical basis is provided, the reasonability of the parameter value is ensured, and the method has the advantages of high reliability and high reliability. Therefore, a carbon dioxide blasting numerical simulation result is close to an actual application result, and the accuracy of the numerical simulation result is improved.

Description

Technical Field

[0001] The invention relates to the technical field of numerical simulation, in particular to a numerical simulation method of carbon dioxide phase change explosion based on FDEM. Background Art

[0002] Carbon dioxide phase change blasting is to use liquid carbon dioxide at a certain pressure and temperature (when the temperature exceeds 31.1℃ and the pressure exceeds 7.38MPa, carbon dioxide becomes a supercritical state) to quickly heat it through a heating tube to make it expand rapidly, generate strong physical pressure, and thus achieve the effect of breaking rocks. Compared with explosive blasting, carbon dioxide blasting has low power and long action time. It does not produce short waves during the blasting process and has little impact on the surrounding environment. The application of this technology is becoming more and more extensive. Carbon dioxide blasting technology has been increasingly used in tunnel (lane) excavation, coal seam permeability enhancement, subway station foundation pit excavation and other fields. However, there is currently insufficient understanding of the rock mass fracture mechanism under the action of carbon dioxide phase change blasting, and it is difficult to capture the rock mass fracture process through experimental means. The finite-discrete element method (FDEM) created by Munjiza can effectively simulate the rock mass fracture process, but how to apply the load of carbon dioxide phase change blasting in the algorithm, how to match the rock mass mechanical parameters with the mesoscopic parameters in the algorithm, so as to ensure the reliability of the numerical simulation results, there are relatively few reports on this aspect.

[0003] The Chinese invention patent "A general numerical simulation method for blasting" (application number 202111534256.7) proposes a numerical simulation method for simulating rock fracture under blasting, but the selection of the peak value and duration of the pressure load is somewhat arbitrary, and does not consider the role of high-pressure gas flowing into the crack in promoting crack expansion (often referred to as "gas wedge effect"). The invention patent "A simulation method for the entire process of rock blasting based on FEM-DEM joint simulation" proposes a simulation method for simulating explosive blasting, but its load form is quite different from carbon dioxide phase change blasting, and does not consider the gas wedge effect. The paper "Finite element-discrete element simulation of rock fracture driven by explosive gas" (DOI: 10.16285 / j.rsm.2015.08.039) realizes the simulation of gas wedge effect based on the FDEM algorithm, but it uses an ideal gas model to describe the explosion gas, which is not suitable for high-pressure and high-temperature carbon dioxide. The paper also does not explain how to match the microscopic parameters in FDEM with the rock parameters.

[0004] How to accurately simulate the rock-breaking process of carbon dioxide phase change blasting, so as to enhance the understanding of the rock-breaking mechanism of carbon dioxide phase change blasting, guide the design of construction parameters, and promote the application of this technology is an important issue. Summary of the invention

[0005] The purpose of the present invention is to provide a numerical simulation method for carbon dioxide phase change explosion based on FDEM to solve the problems raised in the above technical background.

[0006] To achieve the above object, the present invention is implemented through the following technical solutions:

[0007] A numerical simulation method for carbon dioxide phase change explosion based on FDEM, characterized by comprising the following steps:

[0008] The first step is to determine the pV relationship of adiabatic expansion of carbon dioxide based on the Span-Wagner fluid model;

[0009] The second step is to establish a finite-discrete element (FDEM) numerical model and set a non-reflecting boundary;

[0010] The third step is to determine the microscopic parameters in the FDEM model based on the physical and mechanical parameters of the rock mass;

[0011] The fourth step is to apply the carbon dioxide phase change explosion load and solve the finite-discrete element (FDEM) model.

[0012] After high-temperature and high-pressure carbon dioxide is released, its pressure decays as its volume increases, and this process can be approximately regarded as an adiabatic process (isentropic process). The mathematical model that describes the state of the medium and its physical parameters (such as pressure, temperature, density, viscosity, specific heat, and throttling effect coefficient, etc.) is called the state equation. 2 Near supercriticality or at high temperature and high pressure, the physical property parameters show highly nonlinear characteristics. Conventional ideal gas state equations and multi-parameter third-order state equations, such as van der Waals, Redlich-Kwong (RK), Soave (SRK) and Peng-Robinson (PR) equations, cannot meet the requirements of CO2 phase change problems in terms of calculation accuracy and scope of application. In contrast, Span-Wagner starts with the energy equation of the medium and calculates the state characteristics through the Helmholtz free energy. Its physical property parameters can be derived from the energy equation (rather than the explicit function of the state parameter). In the temperature range of 216.6 to 1100K and the pressure range of 0 to 800MPa, the error of the theoretical calculation results is within the engineering allowable range, so it is used as the recommended method of the American National Standard. For this reason, the Span-Wagner fluid model is used to determine the PV relationship of the adiabatic expansion of carbon dioxide.

[0013] The Helmholtz free energy can be represented by two relatively independent variables, namely density ρ and temperature T. Its dimensionless form can be divided into the energy corresponding to the ideal state φ 0 and the residual energy φ r :

[0014] φ(δ, τ)=φ 0(δ, τ)-φ r (δ, τ) (1)

[0015] In the formula, δ = ρ / ρ c is the specific density, τ = T c / T is the reciprocal of the specific temperature; ρ c =468kg / m 3 and T c =304.41K are the critical density and critical temperature of carbon dioxide respectively.

[0016] According to CO 2 The thermodynamic properties of the two parts of energy φ can be obtained o and φ r They are as follows

[0017]

[0018] In the formula, n i d i ,t i 、c i , α i , β i , ε i , γ i 、a i 、b i , A i , B i , C i and D i is a constant, and its value can be found in manuals or literature (Span R, Wagner W. A new equation of state for carbon dioxide covering the fluid region from the triple-point temperature to 1100K at pressures up to 800MPa [J]. Journal of physical and chemical reference data, 1996, 25 (6): 1509-1596.)

[0019] According to the above Span-Wagner equation, CO 2 Expressions for pressure p and specific entropy s:

[0020]

[0021] Where, R = 0.1889241 kJ / (kg·K), Its expanded expression can still be found in the above literature.

[0022] According to the density of carbon dioxide near release ρ 0 and pressure p 0 (Determined by the thickness of the energy dissipation plate or the wall thickness of the fracture tube) Substitute into equations (4) and (5) and solve them together to obtain the specific entropy s 0 Since the release of carbon dioxide is a rapid process, it can generally be regarded as an adiabatic process, and the specific entropy of carbon dioxide remains unchanged, that is, it satisfies the following formula:

[0023]

[0024] Next, enter a series of density values ​​ρ = ρ 1 , ρ 2 ,…,ρ i To the above equation (6), solve the equation to obtain a series of corresponding temperature values ​​T = T 1 , T 2 , …, T i , and then substitute the density and temperature data points into equation (4) to obtain a series of corresponding pressure values ​​p = p 1 , p 2 ,…,p i Since the expansion process of carbon dioxide also follows the law of conservation of mass, when the density is ρ = ρ 1 , ρ 2 ,…,ρ i When , the corresponding volume is V = V 1 , V 2 , …, V i =m / ρ 1 , m / ρ 2 ,…m / ρ i , m is the mass of carbon dioxide. i , p i ) are plotted on the coordinate graph and connected to obtain the adiabatic pV curve of carbon dioxide.

[0025] For the second step, the geometric model of the rock mass is established, the solid triangular unit finite element mesh is divided, and the four-node joint unit is inserted between all adjacent triangular units. This is a common step in FDEM numerical simulation. Then, the non-reflecting boundary is set as needed, and the setting method is referenced (LisjakA, Grasselli G. Combined finite-discrete element analysis of rock slope stability under dynamic loading [C] / / Pan-Am Canadian Geotechnical Society Geotechnical Conference. 2011.).

[0026] For the third step, the microscopic physical and mechanical parameters involved in the FDEM model are solid density ρ, elastic modulus E, Poisson's ratio ν, joint element normal penalty p fn , joint unit tangential penalty p fs , tensile strength of joint element f t , joint unit cohesion c, joint unit internal friction angle φ, joint unit mode I fracture energy density G fI , joint unit type II fracture energy density G fII , contact penalty p n and p s , the key to improving the accuracy of simulation results is to reasonably select the values ​​of these microscopic parameters. In the present invention, the values ​​of the above microscopic parameters are as follows: ρ = ρ m , E=E m , ν=ν m , p fn and p fs The value is 30E or above, f t =σ t,m , c=c m ,φ=φ m , G fII The value of has little effect on the simulation results. fII =3G fI , p n =p s =0.1448p fn . Among them G fI,r is the mode I fracture energy density of the rock mass, K fI,r is the mode I fracture toughness of the rock, E 0，r is the initial modulus of the rock mass, K v is the rock mass integrity index, ρ m 、E m , ν m , σ t,m 、cm ,φ m They are the density, deformation modulus, Poisson's ratio, tensile strength, cohesion, and internal friction angle of the rock mass. When direct measurement data is lacking, the physical and mechanical parameters of the rock mass and the rock mass integrity index K can be used to calculate the rock mass integrity index. v To get the value: ρ m =ρ r , E m =K v E r , ν m =v r , σ t,m =K v σ t,r , c m =c r ,φ m =φ r , ρ r 、E r , ν r , σ t,r 、c r ,φ r They are the density, deformation modulus, Poisson's ratio, tensile strength, cohesion, and internal friction angle of the rock.

[0027] After the carbon dioxide is released from the fracturing tube, it impacts the borehole wall and generates initial radial cracks. Subsequently, the carbon dioxide flows into the cracks, which will cause the cracks to expand further, and the crack expansion interacts with the flow of carbon dioxide. The present invention adopts a simplified method to describe the expansion flow of carbon dioxide, that is, it is assumed that the gas seepage rate is constant in the radial direction of the borehole. According to the experimental results of explosive blasting, the average radial flow velocity of the explosive gas in the cracks is in the range of 160 to 240 m / s. Since there is no measurement result on the flow rate for carbon dioxide phase change blasting as a reference, considering that the carbon dioxide pressure is lower than the explosive gas pressure, its radial flow velocity is taken as 100 m / s.

[0028] Corresponding to the above process, in the FDEM numerical model, the load application method of carbon dioxide phase change explosion is shown in Figure 1 and Figure 2 . Figure 1 is the load time history curve. The front section of the peak is approximated by a straight line to represent the rising process of the hole wall pressure. 0 The value is determined based on the measured results. If there is no measured data, it is 0.2-0.3ms. m The value is taken according to the measured results. When the measured results are missing, the volume of carbon dioxide filling the borehole is substituted into the pV relationship obtained in the first step for calculation; the dotted line after the peak represents the attenuation process of carbon dioxide pressure during the calculation process. Its specific pt relationship cannot be determined before the program starts running. The pressure p is calculated according to the pV relationship obtained in the first step based on the volume V of carbon dioxide at each calculation step. Figure 1 The post-peak section in the process is applied to the crack surface through which the drilling and carbon dioxide flow. The schematic diagram is shown in Figure 2 For details of the fracture network search process, please refer to the literature “Finite element-discrete element simulation of rock mass fracture driven by explosive gas” (DOI:10.16285 / j.rsm.2015.08.039).

[0029] After the load conditions are set, the model is solved according to the finite-discrete element FDEM solution format.

[0030] Compared with the prior art, the present invention has the following beneficial effects:

[0031] The present invention adopts the Span-Wagner fluid model recommended by the American national standard to replace the traditional ideal gas model to calculate the adiabatic pV relationship of carbon dioxide. The description of carbon dioxide pressure is closer to reality. The values ​​of FDEM microscopic parameters are linked to the physical and mechanical indicators of the rock mass, which has a theoretical basis and ensures the rationality of the parameter values, so that the numerical simulation results of carbon dioxide blasting are close to the actual application results, thereby improving the accuracy of the numerical simulation results. BRIEF DESCRIPTION OF THE DRAWINGS

[0032] Figure 1 is the load time history curve of the load application mode of carbon dioxide phase change explosion;

[0033] Figure 2 It is a schematic diagram of gas pressure application;

[0034] Figure 3 is the adiabatic pV curve of carbon dioxide according to an embodiment of the present invention;

[0035] Figure 4 is a numerical model of an embodiment of the present invention;

[0036] Figure 5 is the appearance of the rock mass corresponding to the embodiment of the present invention;

[0037] Figure 6 is the numerical simulation result of the embodiment of the present invention;

[0038] Figure 7 It is the field test result compared with the embodiment of the present invention. DETAILED DESCRIPTION

[0039] In order to better understand the present invention, the present invention is further described below in conjunction with the accompanying drawings and specific embodiments.

[0040] The present invention provides a numerical simulation method for carbon dioxide phase change explosion based on FDEM, and the specific steps are as follows:

[0041] In the first step, the pV relationship of the adiabatic expansion of carbon dioxide is determined based on the Span-Wagner fluid model.

[0042] Taking the Φ89×4×1200 (outer diameter×wall thickness×length, mm) disposable carbon dioxide fracturing tube as an example, its bursting pressure (i.e., the pressure when carbon dioxide is about to be released) is 42.47MPa, the volume is 6.18L, and the density of carbon dioxide when it fills the fracturing tube is 835.3kg / m 3 , then the mass of carbon dioxide is m = 5.162kg. 0 =42.47MPa,ρ=ρ 0 =835.3kg / m 3 Substituting into equations (4) and (5), we can obtain the specific entropy s = s 0 =1.3143 kJ / (kg·K). Next, enter a series of density values ​​ρ=ρ 1 , ρ 2 ,…,ρ i To equation (6), solve the equation to obtain a series of corresponding temperature values ​​T = T 1 , T 2 ,…,T i , and then substitute the density and temperature data points into equation (4) to obtain a series of corresponding pressure values ​​p = p 1 , p 2 ,…,p i . With density ρ = ρ 1 , ρ 2 ,…,ρ i The corresponding volume is V = V 1 , V 2 , …, V i =m / ρ 1 , m / ρ 2 ,…m / ρ i , m=5.126kg. i , p i ) are plotted on the coordinate graph and connected to obtain Figure 3 Adiabatic pV curve of carbon dioxide shown.

[0043] The second step is to establish a finite-discrete element (FDEM) numerical model and set a non-reflecting boundary.

[0044] A 12m×7m numerical model with 21550 triangular elements (the smallest size is about 9mm) and 32167 joint elements is established. Figure 4 The borehole diameter is 0.1 m, the upper side of the model is a free surface, 1 m away from the center of the borehole, and non-reflecting boundaries are set on the other three sides to simulate the propagation of stress waves in infinite rock.

[0045] The third step is to determine the microscopic parameters in the FDEM model based on the physical and mechanical parameters of the rock mass.

[0046] Taking the rock mass in a mining area as an example, the rock type is tuff, the average saturated uniaxial compressive strength of the rock block is Rc=22.6MPa, and the density, deformation modulus, Poisson's ratio, tensile strength, cohesion, and internal friction angle of the rock block are ρ r =2500kg / m 3 , E r =3.91GPa, v = 0.15, σ t，r =1.94MPa,c r =5.6MPa,φ r =38°, the I-mode fracture energy density of the rock is G fI，r =105J / m 2 . Qualitative evaluation of the rock mass from the side appearance of the step (e.g. Figure 5 As shown in the figure, the rock mass integrity is determined to be “relatively complete”, and the rock mass integrity index is estimated to be Kv=0.65.

[0047] According to the above rock mass physical and mechanical parameters, the corresponding FDEM microscopic physical and mechanical parameters are as follows: ρ = ρ m =ρ r =2500kg / m 3 , E=E m =K v E r =0.65×3.91=2.5415GPa, v=v m =v r =0.15, p fn =p fs =30E=30×3.91=117.3GPa, f t =σ t，m =K v σ t，r =0.65×1.94=1.261MPa, c=c m =K v c r =0.65×5.6=3.64MPa, φ=φ m =φ r =38°, p n =p s =0.1448p fn =0.1448×117.3=16.985GPa,G fI =K v ·G fI，r =0.65×105=68.25J / m 2 , G fII =3G fI=3×105=315J / m 2 .

[0048] The fourth step is to apply the carbon dioxide phase change explosion load and solve the finite-discrete element (FDEM) model.

[0049] Carbon dioxide phase change explosion load Figure 1 As shown, the pore wall pressure rise time is taken as t m =0.24ms. The peak value of the borehole pressure is calculated based on the volume of carbon dioxide when the borehole is filled with pm. The volume of carbon dioxide when the borehole is filled with V = 6.833 × 10 -3 m 3 Substitute the volume of the borehole into the pV relationship obtained in the first step to obtain p m =21.0MPa. The radial velocity of carbon dioxide is taken as 100m / s. Figure 1 The dotted line after the peak represents the attenuation process of the carbon dioxide pressure during the calculation process. Its specific pt relationship cannot be determined before the program starts running. The pressure p is calculated based on the pV relationship obtained in the first step by substituting the carbon dioxide volume V at each calculation step into the calculation. Figure 1 The post-peak section in the process is applied to the crack surface through which the drilling and carbon dioxide flow. The schematic diagram is shown in Figure 2 For details of the fracture network search process, please refer to the literature “Finite element-discrete element simulation of rock mass fracture driven by explosive gas” (DOI:10.16285 / j.rsm.2015.08.039).

[0050] After the load conditions are set, the model is solved according to the finite element discrete element FDEM solution format. When the crack no longer expands with the increase of the calculation time step, the program is terminated. The numerical simulation results are as follows: Figure 6 As shown in the figure, the loosening area of ​​carbon dioxide phase change blasting on site is as follows Figure 7 .from Figure 6 It can be seen that the crack extension depth obtained by numerical simulation (i.e., the maximum distance that the crack extends perpendicular to the free surface) is about 2m, which is close to the on-site carbon dioxide blasting result, and the crack width obtained by numerical simulation (i.e., the maximum distance that the crack extends parallel to the free surface) is about 2.5m, which is close to the on-site carbon dioxide blasting drilling spacing. It can be seen that the numerical simulation method provided by the present invention has high accuracy.

[0051] It will be easily understood by those skilled in the art that the above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions and improvements made within the spirit and principles of the present invention should be included in the protection scope of the present invention.

Claims

1. A numerical simulation method for carbon dioxide phase change explosion based on FDEM, characterized in that: The following steps are involved: The first step is to determine the pV relationship of adiabatic expansion of carbon dioxide based on the Span-Wagner fluid model; The second step is to establish a finite-discrete element (FDEM) numerical model and set a non-reflecting boundary; The third step is to determine the microscopic parameters in the FDEM model based on the physical and mechanical parameters of the rock mass; The fourth step is to apply the carbon dioxide phase change explosion load and solve the finite-discrete element (FDEM) model.

2. A numerical simulation method for carbon dioxide phase change explosion based on FDEM as claimed in claim 1, characterized in that: The expressions of CO2 pressure p and specific entropy s are: Where, R = 0.1889241 kJ / (kg·K), According to the density ρ0 and pressure p0 of carbon dioxide near release, substitute them into equations (4) and (5) and solve them together to obtain the specific entropy s0. Since the release of carbon dioxide is a rapid process, it is approximately regarded as an adiabatic process, and the specific entropy of carbon dioxide remains unchanged, that is, it satisfies the following equation: Next, enter a series of density values ​​ρ = ρ1, ρ2, ..., ρ i To equation (6), where δ = ρ / ρ c is the specific density, τ = T c / T is the reciprocal of the specific temperature; ρ c =468kg / m 3 and T c =304.41K are the critical density and critical temperature of carbon dioxide respectively; solving the equations we can get a series of corresponding temperature values ​​T=T1, T2, ...T i , and then substitute the density and temperature data points into equation (4) to obtain a series of corresponding pressure values ​​p = p1, p2, ..., p i ; Since the expansion process of carbon dioxide also follows the law of conservation of mass, when the density is ρ = ρ1, ρ2, ..., ρ i When , the corresponding volume is V = V1, V2, ..., V i =m / ρ1, m / ρ2, …m / ρ i , m is the mass of carbon dioxide; the data point (V i , p i ) are plotted on the coordinate graph and connected to obtain the adiabatic pV curve of carbon dioxide.

3. The FDEM-based numerical simulation method for carbon dioxide phase change explosion according to claim 1, characterized in that: The microscopic physical and mechanical parameters involved in the FDEM model include solid density ρ, elastic modulus E, Poisson's ratio ν, joint element normal penalty p fn , joint unit tangential penalty p fs , tensile strength of joint element f t , joint unit cohesion c, joint unit internal friction angle φ, joint unit mode I fracture energy density G fI , joint unit type II fracture energy density G fII , contact penalty p n and p s , Where ρ = ρ m , E=E m , ν=ν m , p fn and p fs The value is 30E or above, f t =σ t,m , c=c m ,φ=φ m , G fII =3G fI , p n =p s =0.1448p fn ; where G fI,r is the mode I fracture energy density of the rock mass, K fI,r is the mode I fracture toughness of the rock, E 0，r is the initial modulus of the rock mass, K v is the rock mass integrity index, ρ m 、E m , ν m , σ t,m 、c m ,φ m They are the density, deformation modulus, Poisson's ratio, tensile strength, cohesion and internal friction angle of the rock mass.

4. A numerical simulation method for carbon dioxide phase change explosion based on FDEM as claimed in claim 3, characterized in that: When direct measurement data are lacking, the physical and mechanical parameters of the rock mass and the rock mass integrity index K are used to estimate the rock mass integrity. v To get the value: ρ m =ρ r , E m =K v E r , ν m =v r , σ t,m =K v σ t,r , c m =c r ,φ m =φ r , ρ r 、E r , ν r , σ t,r 、c r ,φ r They are the density, deformation modulus, Poisson's ratio, tensile strength, cohesion, and internal friction angle of the rock.

5. The FDEM-based numerical simulation method for carbon dioxide phase change explosion as claimed in claim 3, characterized in that: For the fourth step, the CO2 phase change explosion load is applied based on the CO2 expansion flow described by the following method: Assuming that the gas seepage rate is constant in the radial direction of the borehole, the average radial flow velocity of the explosion gas in the fracture is in the range of 160 to 240 m / s, and considering that the CO2 pressure is lower than the explosion gas pressure, its radial flow velocity is taken as 100 m / s.

Images