CN104657581A - Water-sand two-phase permeability parameter calculation method in water driving sand test - Google Patents

Abstract

The invention relates to a water-sand two-phase permeability parameter calculation method in a water driving sand test. The fluidity of a water phase and a sand phase, a non-Darcy flow beta factor and an accelerated speed coefficient are calculated according to a pressure gradient time sequence, a water flow time sequence and a sand flow time sequence. The water-sand two-phase permeability parameter calculation method is realized by the following steps: calculating a water permeability speed time sequence, a sand permeability speed time sequence and the pressure gradient time sequence; establishing time sequences of water-phase and sand-phase permeability speed local derivatives; calculating volume fractions of a water phase and a sand phase of each sampling moment; establishing an algebraic equation of a water flow speed and a sand flow speed of each sampling moment; calculating time sequences of the non-Darcy flow beta factors and the accelerated speed coefficients of the water phase and the sand phase; calculating the time sequences of water-phase and sand-phase permeability speeds; and comparing errors of calculation values of the water-phase and sand-phase permeability speeds and actually-measured values to obtain a testing value of a permeability parameter of water-phase and sand-phase flows of crushed rocks.

Description

Water husky two-phase perviousness Parameters Calculation method in the test of a kind of water drive sand

Technical field

The present invention relates to collery's surface movement and control, be specifically related to water husky two-phase flow perviousness Parameters Calculation method in the test of a kind of water drive sand.

Background technology

At Mu us dese, coal seam has the advantages that buried depth is little, earth's surface desilting is thick.In coal mining process, the layer of sand on earth's surface is migrated to goaf along coal seam overlying strata with surface water.In order to study in the recovery process of coal seam because Overburden Rock Failure causes burst husky mechanism, control gushing water of gushing water to burst husky generation or alleviate gushing water and to burst husky harm, carry out fractured rock water drive sand test (namely utilizing current to expel fine sand test in fractured rock hole).Due in the process of water drive sand, in fractured rock hole, the volume fraction of sand reduces gradually, and perviousness parameter that is mutually husky and aqueous phase changes in time.We calculate perviousness parameter that is mutually husky in the husky process of water drive and aqueous phase, and this patent proposes a kind of method of mobility based on pressure gradient time series, discharge time series and sand drift amount time series calculating aqueous phase and husky phase, non-Darcy flow β-factor and acceleration factor.

At present, in the husky process of the test of water drive, the computing method of perviousness parameter are carried out with reference to the method for water displacing oil test.Because sand belongs to Discontinuous transmission in essence, there is significant difference with the seepage flow of oil in the migration of sand in fractured rock hole in movement mechanism, mode of motion and drag characteristic etc.Therefore, it is not proper that the method with reference to water displacing oil test calculates this way of perviousness parameter in the husky process of the test of water drive, and it is necessary for building a kind of computing method being different from perviousness parameter in the husky process of the test of water drive of water displacing oil test, is also urgent.

Summary of the invention

The object of this invention is to provide water husky two-phase flow perviousness Parameters Calculation method in the test of a kind of water drive sand, calculating the mobility of aqueous phase and husky phase, non-Darcy flow β-factor and acceleration factor according to pressure gradient time series, discharge time series and sand drift amount time series, providing condition for analyzing water husky two-phase perviousness parameter Changing Pattern in the husky process of water drive; Use computing method of the present invention, can set up the quantitative relationship between sand volumes mark in the mobility of aqueous phase and husky phase, non-Darcy flow β-factor and acceleration factor and fractured rock, in order to go together, similar test provides technical support.

The present invention is achieved by the following technical solutions: water husky two-phase flow perviousness Parameters Calculation method in the test of a kind of water drive sand, the test of water drive sand is the fine sand mixing certain mass in fractured rock, put into the permeameter of special purpose, utilize vane pump, pressurized strut or other equipment to inject the current of constant pressure or constant rate continuously to the upper end of permeameter or lower end, fine sand is with the process of the migration of current in fractured rock hole.The most significant feature of the husky process of water drive is that the volume fraction of sand in fractured rock hole changes in time, and thus, aqueous phase and husky phase perviousness parameter (mobility, non-Darcy flow β-factor and acceleration factor) also change in time.

Water is Newton fluid, and momentum conservation equation is

$m_1{c}_{1a}\frac{\∂V_1}{\∂t}=-\frac{\∂p}{\∂x}-\frac{{\mathrm{\μ}}_1}{k_1}{V}_1-m_1{\mathrm{\β}}_1\left|V_1\right|V_1---\left(1\right)$

Wherein, t is the time; X is volume coordinate; P is water pressure; m ₁for the mass density of water; V ₁for the percolation flow velocity of water; for percolation flow velocity V ₁local derivative; μ ₁for the kinetic viscosity of water; k ₁, β ₁and c _1abe respectively water phase permeability, non-Darcy flow β-factor and acceleration factor. for pressure gradient, use symbol G _prepresent.

Sand is non-Newtonian fluid, and momentum conservation equation is

$m_2{c}_{2a}\frac{\∂V_2}{\∂t}=-\frac{\∂p}{\∂x}-\frac{{\mathrm{\μ}}_{2e}}{k_{2e}}{V}_2^n-m_2{\mathrm{\β}}_2\left|V_2\right|V_2---\left(2\right)$

Wherein, m ₂for the mass density of wet sand; V ₂for the percolation flow velocity of sand; for percolation flow velocity V ₂local derivative; μ _2efor the virtual viscosity of sand; k _2e, β ₂and c _2abe respectively husky phase effective permeability, non-Darcy flow β-factor and acceleration factor.

Use aqueous phase mobility

$I_1=\frac{k_1}{{\mathrm{\μ}}_1}---\left(3\right)$

Replace the kinetic viscosity of water phase permeability and water, and use G _prepresent pressure gradient, formula (1) can be deformed into

$m_1{c}_{1a}\frac{\∂V_1}{\∂t}=-\frac{\∂p}{\∂x}-\frac{1}{I_1}{V}_1-m_1{\mathrm{\β}}_1\left|V_1\right|V_1---\left(4\right)$

With sand effective mobility mutually

$I_{2e}=\frac{k_{2e}}{{\mathrm{\μ}}_{2e}}---\left(5\right)$

Replace husky phase effective permeability and wet husky virtual viscosity, and use G _prepresent pressure gradient, formula (2) can be deformed into

$m_2{c}_{2a}\frac{\∂V_2}{\∂t}=-G_p-\frac{1}{I_{2e}}{V}_2^n+m_2{\mathrm{\β}}_2\left|V_2\right|V_2---\left(6\right)$

In permeability test, be the sampling period gather pressure and flow with τ, through simple operation, easily obtain pressure gradient time series

$G_p^{t_i}(t_i=\mathrm{i\τ},i=\mathrm{1,2},...,N)$

Aqueous phase percolation flow velocity time series

$V_1^{t_i}(t_i=\mathrm{i\τ},i=\mathrm{1,2},...,N)$

Husky phase percolation flow velocity time series

$V_2^{t_i}(t_i=\mathrm{i\τ},i=\mathrm{1,2},...,N)$

Utilize difference coefficient to replace difference quotient, can percolation flow velocity V be obtained ₁and V ₂local derivative with time series

$a_1^{t_i}=\frac{\∂V_1}{\∂t}=\frac{V_1^{t_{i+1}}-V_1^{t_i}}{\mathrm{\τ}},i=\mathrm{1,2},...,N-1$

With

$a_2^{t_i}=\frac{\∂V_2}{\∂t}=\frac{V_2^{t_{i+1}}-V_2^{t_i}}{\mathrm{\τ}},i=\mathrm{1,2},...,N-1$

Because current continuous driving is for the sand in fractured rock hole, therefore the volume fraction Y of sand in hole ₂change in time, and the mobility of aqueous phase and husky phase (effective mobility), non-Darcy flow β-factor and acceleration factor are with volume fraction Y ₂change.Therefore, need to build sampling instant volume fraction Y ₂, the mobility (effective mobility) of aqueous phase and husky phase, non-Darcy flow β-factor and acceleration factor

$Y_2^{t_i},I_1^{t_i},{\mathrm{\β}}_1^{t_i},c_{1a}^{t_i},I_{2e}^{t_i},{\mathrm{\β}}_2^{t_i},c_{2a}^{t_i}(t_i=\mathrm{i\τ},i=\mathrm{1,2},...,N-1)$

Computing method.

At each sampling instant t _i=i τ, by momentum conservation equation, obtains

$m_1{c}_{1a}{a}^{t_i}=-G^{t_i}-\frac{1}{I_1^{t_i}}{V}_1^{t_i}-m_1{\mathrm{\β}}_1^{t_i}\left|V_1^{t_i}\right|V_1^{t_i},i=\mathrm{1,2},...,N-1---\left(7\right)$

$m_2{c}_{2a}^{t_i}{a}_2^{t_i}=-G_p^{t_i}-\frac{1}{I_{2e}}{\left(V_2^{t_i}\right)}^n-m_2{\mathrm{\β}}_2^{t_i}\left|V_2^{t_i}\right|V_2^{t_i},i=\mathrm{1,2},...,N-1---\left(8\right)$

In system of equations (7), unknown quantity number is 3 (N-1), and the number of equation is N-1; Therefore need supplementary equation, just can obtain this solution of equations.The situation of system of equations (8) is the same with system of equations (7).

Suppose mobility, meet power exponent relation between non-Darcy flow β-factor and acceleration factor, namely

$\left\{\begin{array}{c}\frac{{\mathrm{\β}}_1}{{\mathrm{\β}}_{1r}}={\left(\frac{I_1}{I_{1r}}\right)}^{-n_{1\mathrm{\β}}}\\ \frac{c_{1a}}{c_{1a}^r}={\left(\frac{I_1}{I_{1r}}\right)}^{-n_{1c}}\end{array}\right.---\left(9\right)$

With

$\left\{\begin{array}{c}\frac{{\mathrm{\β}}_2}{{\mathrm{\β}}_{2r}}={\left(\frac{I_2}{I_{2r}}\right)}^{-n_{2\mathrm{\β}}}\\ \frac{c_{2a}}{c_{2a}^r}={\left(\frac{I_1}{I_{2r}}\right)}^{-n_{2c}}\end{array}\right.---\left(10\right)$

Then all have in each sampling instant

$\left\{\begin{array}{c}\frac{{\mathrm{\β}}_1^{t_i}}{{\mathrm{\β}}_{1r}}={\left(\frac{I_1^{t_i}}{I_{1r}}\right)}^{-n_{1\mathrm{\β}}}\\ \frac{c_{1a}^{t_i}}{c_{1a}^r}={\left(\frac{I_1^{t_i}}{I_{1r}}\right)}^{-n_{1c}}\end{array}\right.,i=\mathrm{0,1,2},...,N-1---\left(11\right)$

With

$\left\{\begin{array}{c}\frac{{\mathrm{\β}}_2^{t_i}}{{\mathrm{\β}}_{2r}}={\left(\frac{I_{2e}^{t_i}}{I_{2r}}\right)}^{-n_{2\mathrm{\β}}}\\ \frac{c_{2a}^{t_i}}{c_{2a}^r}={\left(\frac{I_{2e}^{t_i}}{I_{2r}}\right)}^{-n_{2c}}\end{array}\right.,i=\mathrm{0,1,2},...,N-1---\left(12\right)$

Wherein, the reference value I of aqueous phase mobility, non-Darcy flow β-factor and acceleration factor _1r, β _1rwith power exponent n _{1 β}and n _1c, husky phase mobility, non-Darcy flow β-factor and acceleration factor reference value I _2r, β _2rwith power exponent n _{2 β}and n _2cthe optimization method such as genetic algorithm, ant group algorithm can be utilized to determine.

Respectively formula (11) and formula (12) are substituted into formula (7) and formula (8), obtain

$m_1{c}_{1a}^r{\left(\frac{I_1^{t_i}}{I_r}\right)}^{-n_{1c}}{a}_1^{\left(i\right)}=-G_p^{t_i}-\frac{1}{I_1^{t_i}}{V}_1^{t_i}-m_1{\mathrm{\β}}_{1r}^{t_i}{\left(\frac{I_1^{t_i}}{I_{1r}}\right)}^{-n_{1\mathrm{\β}}}\left|V_1^{t_i}\right|V_1^{t_i},i=\mathrm{0,1,2},...,N-1---\left(13\right)$

With

$m_2{c}_{2a}^r{\left(\frac{I_{2e}^{t_i}}{I_{2r}}\right)}^{-n_{2c}}{a}_2^{t_i}=-G_p^{t_i}-\frac{1}{I_{2e}^{t_i}}{\left(V_2^{t_i}\right)}^n-m_2{\mathrm{\β}}_{2r}^{t_i}{\left(\frac{I_{2e}^{t_i}}{I_{2r}}\right)}^{-n_{2\mathrm{\β}}}\left|V_2^{t_i}\right|V_2^{t_i},i=\mathrm{0,1,2},...,N-1---\left(14\right)$

Utilize Newton tangent method or other solve the root of unitary nonlinear equation, obtain each sampling instant aqueous phase mobility with effective mobility of husky phase i=0,1,2 ..., N-1.According to formula (11) and formula (12), calculate non-Darcy flow β-factor and the acceleration factor of each sampling instant.

Water husky two-phase flow perviousness Parameters Calculation method in the test of a kind of water drive sand, namely the mobility of aqueous phase and husky phase, non-Darcy flow β-factor and acceleration factor is calculated according to pressure gradient time series, discharge time series and sand drift amount time series, it is characterized in that, it is realized by following steps:

Step 1) to gather discharge, sand drift amount and pressure differential at equal intervals, calculate the water percolation flow velocity of each sampling instant, husky percolation flow velocity and pressure gradient, obtain water percolation flow velocity time series, husky percolation flow velocity time series and pressure gradient time series; According to kinematic relation, and replace difference quotient by difference coefficient, set up the time series of aqueous phase and husky phase percolation flow velocity local derivative;

Step 1.1) gather discharge, sand drift amount and pressure differential with τ at equal intervals, pressure gradient G can be obtained respectively by simple algebraic operation _p, water percolation flow velocity V ₁with husky percolation flow velocity V ₂time series, namely

$G_p^{\left(i\right)}=G_p{|}_{t=\mathrm{i\τ}},(i=\mathrm{0,1,2},...N)---\left(18\right)$

$V_1^{\left(i\right)}=V_1{|}_{t=\mathrm{i\τ}},(i=\mathrm{0,1,2},...N)---\left(19\right)$

$V_2^{\left(i\right)}=V_2{|}_{t=\mathrm{i\τ}},(i=\mathrm{0,1,2},...N)---\left(20\right)$

Step 1.2) replace difference quotient by difference coefficient, calculate percolation flow velocity V ₁and V ₂local derivative with time series, namely

$a_1^{t_i}=\frac{\∂V_1}{\∂t}=\frac{V_1^{t_{i+1}}-V_1^{t_i}}{\mathrm{\τ}},i=\mathrm{1,2},...,N-1---\left(21\right)$

With

$a_2^{t_i}=\frac{\∂V_2}{\∂t}=\frac{V_2^{t_{i+1}}-V_2^{t_i}}{\mathrm{\τ}},i=\mathrm{1,2},...,N-1---\left(22\right)$

Need in the above process to measure wet husky viscosity index n, namely

The first step, utilizes rotational viscosimeter to test fine sand in different strain rate under moment of torsion, shear stress τ _sandwith apparent viscosity μ _a, draw shear stress-Strain rate curve.

Second step, carries out power exponential function matching with Excel software to shear stress-rate of strain scatter diagram,

${\mathrm{\τ}}_{\mathrm{sand}}=C_{\mathrm{sand}}{\stackrel{\·}{\mathrm{\γ}}}^n$

Obtain consistency index C _sandand viscosity index.

Step 2) volume fraction of each sampling instant aqueous phase and husky phase is calculated according to permeameter size, rock sample quality, sandy amount etc.;

Quality is M by the first step _rockbroken rock sample and quality be M _sandfine sand blending evenly as the sample permeating rock sample;

Second step, puts into the permeameter of cylinder barrel diameter 2a by sample, and measures the height of the stacking naturally h of broken rock sample fine sand admixture _uncomp;

3rd step, utilizes Material Testing Machine to load admixture, when sample height reaches preset value h ₀in time, stops loading, and keeps rock sample constant height;

4th step, is calculated as follows fractured rock factor of porosity

$\mathrm{\φ}=1-\frac{M_{\mathrm{rock}}}{m_{\mathrm{rock}}\mathrm{\π}{a}^2{h}_0}---\left(15\right)$

Wherein, m _rockfor broken rock sample parent, the i.e. mass density of rock block;

5th step, calculates the volume fraction before infiltration.Before infiltration in fractured rock hole, the volume fraction of air is calculated as follows

$Y_1^0=1-\frac{\frac{M_{\mathrm{sand}}}{m_2}}{\mathrm{\π}{a}^2{h}_0-\frac{M_{\mathrm{rock}}}{m_{\mathrm{rock}}}}---\left(16\right)$

Husky volume fraction is

$Y_2^0=\frac{\frac{M_{\mathrm{sand}}}{m_2}}{\mathrm{\π}{a}^2{h}_0-\frac{M_{\mathrm{rock}}}{m_{\mathrm{rock}}}}---\left(17\right)$

6th step, calculates initial volume mark

Before seepage flow starts, carry out water saturation to rock sample, water drive has replaced air; If water injection pressure is very little, tiny loss is little, can think the initial volume mark of water equal the volume fraction of fine sand equal

Step 3) utilize the momentum conservation equation of aqueous phase, husky phase, and utilize mobility, power exponent relation between non-Darcy flow β-factor and acceleration factor, set up the algebraic equation of each sampling instant current degree and sand drift degree respectively; Utilize Newton tangent method Solving Algebraic Equation, obtain the time series of current degree and sand drift degree;

Step 3.1) relational expression between setting perviousness parameter

The perviousness parameter of aqueous phase and husky phase medium is set as power exponent relation, namely all has in each sampling instant

$\left\{\begin{array}{c}\frac{{\mathrm{\β}}_1^{t_i}}{{\mathrm{\β}}_{1r}}={\left(\frac{I_1^{t_i}}{I_{1r}}\right)}^{-n_{1\mathrm{\β}}}\\ \frac{c_{1a}^{t_i}}{c_{1a}^r}={\left(\frac{I_1^{t_i}}{I_{1r}}\right)}^{-n_{1c}}\end{array}\right.,i=\mathrm{0,1,2},...,N-1---\left(23\right)$

With

$\left\{\begin{array}{c}\frac{{\mathrm{\β}}_2^{t_i}}{{\mathrm{\β}}_{2r}}={\left(\frac{I_{2e}^{t_i}}{I_{2r}}\right)}^{-n_{2\mathrm{\β}}}\\ \frac{c_{2a}^{t_i}}{c_{2a}^r}={\left(\frac{I_{2e}^{t_i}}{I_{2r}}\right)}^{-n_{2c}}\end{array}\right.,i=\mathrm{0,1,2},...,N-1---\left(24\right)$

Wherein, the reference value of aqueous phase mobility, non-Darcy flow β-factor and acceleration factor is I _1r, β _1rwith , power exponent n _{1 β}and n _1c, husky phase mobility, non-Darcy flow β-factor and acceleration factor reference value I _2r, β _2rwith , power exponent n _{2 β}and n _2c;

Step 3.2) set up the closed system of equations of mobility

Power exponent relation between the perviousness parameter of aqueous phase and husky phase medium and momentum conservation equation simultaneous, by formula (23), formula (24) and respectively formula (13), formula (14) simultaneous, obtain the system of equations that aqueous phase and husky phase mobility meet respectively

$m_1{c}_{1a}^r{\left(\frac{I_1^{t_i}}{I_r}\right)}^{-n_{1c}}{a}_1^{\left(i\right)}=-G_p^{t_i}-\frac{1}{I_1^{t_i}}{V}_1^{t_i}-m_1{\mathrm{\β}}_{1r}^{t_i}{\left(\frac{I_1^{t_i}}{I_{1r}}\right)}^{-n_{1\mathrm{\β}}}\left|V_1^{t_i}\right|V_1^{t_i},i=\mathrm{0,1,2},...,N-1---\left(25\right)$

With

$m_2{c}_{2a}^r{\left(\frac{I_{2e}^{t_i}}{I_{2r}}\right)}^{-n_{2c}}{a}_2^{t_i}=-G_p^{t_i}-\frac{1}{I_{2e}^{t_i}}{\left(V_2^{t_i}\right)}^n-m_2{\mathrm{\β}}_{2r}^{t_i}{\left(\frac{I_{2e}^{t_i}}{I_{2r}}\right)}^{-n_{2\mathrm{\β}}}\left|V_2^{t_i}\right|V_2^{t_i},i=\mathrm{0,1,2},...,N-1---\left(26\right)$

Step 3.3) primary Calculation sampling instant aqueous phase and husky phase mobility

Step 3.3.1), design iteration form, introduces function

$f_{1i}\left(I_1^{t_i}\right)=m_1{\mathrm{\β}}_{1r}^{t_i}{\left(\frac{I_1^{t_i}}{I_{1r}}\right)}^{-n_{1\mathrm{\β}}}\left|V_1^{t_i}\right|V_1^{t_i}+\frac{1}{I_1^{t_i}}{V}_1^{t_i}+m_1{c}_{1a}^r{\left(\frac{I_1^{t_i}}{I_r}\right)}^{-n_{1c}}{a}_1^{\left(i\right)}+G_p^{t_i},i=\mathrm{0,1,2},...,N-1---\left(27\right)$

$f_{2i}\left(I_{2e}^{t_i}\right)=m_2{\mathrm{\β}}_{2r}^{t_i}{\left(\frac{I_{2e}^{t_i}}{I_{2r}}\right)}^{-n_{2\mathrm{\β}}}\left|V_2^{t_i}\right|V_2^{t_i}+\frac{1}{I_{2e}^{t_i}}{\left(V_2^{t_i}\right)}^n+m_2{c}_{2a}^r{\left(\frac{I_{2e}^{t_i}}{I_{2r}}\right)}^{-n_{2c}}{a}_2^{t_i}+G_p^{t_i},i=\mathrm{0,1,2},...,N-1---\left(28\right)$

System of equations (27) and (28) are deformed into respectively

$f_{1i}\left(I_1^{t_i}\right)=0,i=\mathrm{0,1,2},...,N-1---\left(29\right)$

$f_{2i}\left(I_{2e}^{t_i}\right)=0,i=\mathrm{0,1,2},...,N-1---\left(30\right)$

Utilize the Iteration of Newton tangent method equationof structure group (29) and (30);

Step 3.3.2) determine iteration initial value with

Get the stable percolation equation root of Newton fluid as iteration initial value, namely

$-G^{t_i}-\frac{1}{I_1^{*}}{V}_1^{t_i}-m_1{\mathrm{\β}}_1^{t_i}\left|V_1^{t_i}\right|V_1^{t_i}=0,i=\mathrm{1,2},...,N-1---\left(31\right)$

$-G_p^{t_i}-\frac{1}{I_{2e}^{*}}{V}_2^{t_i}-m_2{\mathrm{\β}}_2^{t_i}\left|V_2^{t_i}\right|V_2^{t_i}=0,i=\mathrm{1,2},...,N-1---\left(32\right)$

Step 3.3.3) according to the Iteration of first two steps and iteration initial value, calculating sampling moment aqueous phase and husky phase mobility with

Step 4) utilize mobility, power exponent relation between non-Darcy flow β-factor and acceleration factor, calculate the non-Darcy flow β-factor of aqueous phase and husky phase and the time series of acceleration factor;

Step 4.1) according to formula (23) and formula (24), calculate aqueous phase and husky phase non-Darcy flow β-factor and acceleration factor time series, (i=1,2 ..., N-1).

Step 5) utilize four step Runge-Kutta respectively to the numerical solution of the momentum conservation equation of aqueous phase, husky phase, obtain the time series of aqueous phase and husky phase percolation flow velocity;

Step 5.1) calculate aqueous phase and husky phase percolation flow velocity

Step 5.1.1) build the external function of momentum conservation equation, by momentum conservation equation (13) and (14) distortion, obtain

$\frac{{\∂V}_1}{\∂t}=-\frac{1}{m_1{c}_{1a}}(G_p+\frac{1}{I_1}{V}_1+m_1{\mathrm{\β}}_1|V_1\left|V_1\right)---\left(33\right)$

$\frac{{\∂V}_2}{\∂t}=-\frac{1}{m_2{c}_{2a}}(G_p+\frac{1}{I_{2e}}{V}_2^n+m_2{\mathrm{\β}}_2|V_2\left|V_2\right)---\left(34\right)$

Respectively to pressure gradient G _p, mobility I ₁and I _2e, non-Darcy flow β-factor β ₁and β ₂, acceleration factor c _1aand c _2acarry out three knot interpolations between the whole district, the external function of equationof structure (33) and (34), for any time t ∈ (t _k, t _k+2), have

$G_p=G_p^{t_k}\frac{(t-t_{k+1})(t-t_{k+2})}{(t_k-t_{k+1})(t_k-t_{k+2})}+G_p^{t_{k+1}}\frac{(t-t_k)(t-t_{k+2})}{(t_{k+1}-t_k)(t_{k+1}-t_{k+2})}+G_p^{t_{k+2}}\frac{(t-t_k)(t-t_{k+1})}{(t_{k+2}-t_k)(t_{k+2}-t_{k+1})}---\left(35\right)$

$I_1=I_1^{t_k}\frac{(t-t_{k+1})(t-t_{k+2})}{(t_k-t_{k+1})(t_k-t_{k+2})}+I_1^{t_{k+1}}\frac{(t-t_k)(t-t_{k+2})}{(t_{k+1}-t_k)(t_{k+1}-t_{k+2})}+I_1^{t_{k+2}}\frac{(t-t_k)(t-t_{k+1})}{(t_{k+2}-t_k)(t_{k+2}-t_{k+1})}---\left(36\right)$

$I_{2e}=I_{2e}^{t_k}\frac{(t-t_{k+1})(t-t_{k+2})}{(t_k-t_{k+1})(t_k-t_{k+2})}+I_{2e}^{t_{k+1}}\frac{(t-t_k)(t-t_{k+2})}{(t_{k+1}-t_k)(t_{k+1}-t_{k+2})}+I_{2e}^{t_{k+2}}\frac{(t-t_k)(t-t_{k+1})}{(t_{k+2}-t_k)(t_{k+2}-t_{k+1})}---\left(37\right)$

${\mathrm{\β}}_1={\mathrm{\β}}_1^{t_k}\frac{(t-t_{k+1})(t-t_{k+2})}{(t_k-t_{k+1})(t_k-t_{k+2})}+{\mathrm{\β}}_1^{t_{k+1}}\frac{(t-t_k)(t-t_{k+2})}{(t_{k+1}-t_k)(t_{k+1}-t_{k+2})}+{\mathrm{\β}}_1^{t_{k+2}}\frac{(t-t_k)(t-t_{k+1})}{(t_{k+2}-t_k)(t_{k+2}-t_{k+1})}---\left(38\right)$

${\mathrm{\β}}_2={\mathrm{\β}}_2^{t_k}\frac{(t-t_{k+1})(t-t_{k+2})}{(t_k-t_{k+1})(t_k-t_{k+2})}+{\mathrm{\β}}_2^{t_{k+1}}\frac{(t-t_k)(t-t_{k+2})}{(t_{k+1}-t_k)(t_{k+1}-t_{k+2})}+{\mathrm{\β}}_2^{t_{k+2}}\frac{(t-t_k)(t-t_{k+1})}{(t_{k+2}-t_k)(t_{k+2}-t_{k+1})}---\left(39\right)$

$c_{a1}=c_{a1}^{t_k}\frac{(t-t_{k+1})(t-t_{k+2})}{(t_k-t_{k+1})(t_k-t_{k+2})}+c_{a1}^{t_{k+1}}\frac{(t-t_k)(t-t_{k+2})}{(t_{k+1}-t_k)(t_{k+1}-t_{k+2})}+c_{a1}^{t_{k+2}}\frac{(t-t_k)(t-t_{k+1})}{(t_{k+2}-t_k)(t_{k+2}-t_{k+1})}---\left(40\right)$

$c_{a2}=c_{a2}^{t_k}\frac{(t-t_{k+1})(t-t_{k+2})}{(t_k-t_{k+1})(t_k-t_{k+2})}+c_{a2}^{t_{k+1}}\frac{(t-t_k)(t-t_{k+2})}{(t_{k+1}-t_k)(t_{k+1}-t_{k+2})}+c_{a2}^{t_{k+2}}\frac{(t-t_k)(t-t_{k+1})}{(t_{k+2}-t_k)(t_{k+2}-t_{k+1})}---\left(41\right)$

Step 5.1.2) be step-length with τ, utilize the four step Runge-Kutta of variable step to ask the numerical solution of equation (33) and (34), obtain percolation flow velocity calculated value time series with (i=0,1,2, L, N-1).

Step 6) compare aqueous phase and the husky calculated value of phase percolation flow velocity and the error of measured value, if each sampling instant percolation flow velocity error sum of squares is greater than certain preset value, parameter in adjustment mobility, non-Darcy flow β-factor and acceleration factor relational expression, recalculates the time series of the mobility of aqueous phase and husky phase, non-Darcy flow β-factor and acceleration factor; If each sampling instant percolation flow velocity error sum of squares is less than or equal to certain preset value, then using the test value of the result of calculation of mobility, non-Darcy flow β-factor and acceleration factor as the perviousness parameter of the husky two-phase flow of fractured rock water.

Calculate the numerical solution of percolation flow velocity (i=0,1,2, L, N-1) and test figure difference between (i=0,1,2, L, N-1)

$E_{\mathrm{rr}}=\max \{\frac{1}{N-1}\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{(1-\frac{{\tilde{V}}_1^{t_i}}{V_1^{t_i}})}^2,\frac{1}{N-1}\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{(1-\frac{{\tilde{V}}_2^{t_i}}{V_2^{t_i}})}^2\}---\left(42\right)$

If E _rrbe greater than the permissible value of a certain prior setting then adjust the reference value I of aqueous phase mobility, non-Darcy flow β-factor and acceleration factor _1r, β _1rwith power exponent n _{1 β}and n _1c, husky phase mobility, non-Darcy flow β-factor and acceleration factor reference value I _2r, β _2rwith power exponent n _{2 β}and n _2cvalue, recalculate mobility non-Darcy flow β-factor, acceleration factor and percolation flow velocity until E _rrbe not more than

The invention has the beneficial effects as follows: (1) integrated application Newton tangent method, Runge-Kutta method and unitary 3 Lagrange interpolation methods calculate the husky process water of water drive husky two-phase perviousness parameter, i.e. mobility, non-Darcy flow β-factor and acceleration factor; (2) pressure in process of osmosis (poor), discharge and husky turnover rate only need be carried out periodic sampling by experimenter, produce one to comprise four column datas (first is classified as the time, second is classified as pressure/pressure differential, 3rd is classified as discharge, 4th is classified as sand drift vector) text or Data file, just can revise some parameter in a program and carry out parameter optimization and computing permeability.These parameters comprise the reference value I of the mass density of liquid, kinetic viscosity, rock sample diameter, thickness, quality, sampling period, data length, aqueous phase mobility, non-Darcy flow β-factor and acceleration factor _1r, β _1rwith power exponent n _{1 β}and n _1c, husky phase mobility, non-Darcy flow β-factor and acceleration factor reference value I _2r, β _2rwith power exponent n _{2 β}and n _2cvalue etc.

Accompanying drawing explanation

Below in conjunction with drawings and Examples, the invention will be further described.

Fig. 1 is the time-serial position figure of pressure gradient, water percolation flow velocity and husky percolation flow velocity.

Fig. 2 is aqueous phase mobility, non-Darcy flow β-factor and acceleration factor time series.

Fig. 3 is husky phase mobility, non-Darcy flow β-factor and acceleration factor time series.

Fig. 4 is water percolation flow velocity time series.

Embodiment

1 experimental data processing

Test figure is changed into four row, first is classified as the time, and unit is s; Second is classified as pressure (poor), and unit is MPa; 3rd is classified as discharge, and unit is l/h, and the 4th is classified as husky turnover rate kg/h.According to this three columns group, the time-serial position of pressure gradient, water percolation flow velocity and husky percolation flow velocity can be drawn, see Fig. 1.

2 input relevant parameters

Input relevant parameters, comprises following a few class:

(1) the mass density m of water ₁, kinetic viscosity μ ₁;

(2) husky mass density m ₂, consistency index C and power exponent n;

(2) rock sample radius r _rock, naturally stack height h ₀, height h, fractured rock mass M after compression _rock, husky mass M _sand;

(3) sampling period τ, data length N;

(4) limits of error ε of Newton tangent method ₁, Runge-Kutta method limits of error ε, the numerical solution of percolation flow velocity (i=0,1,2, L, N-1) and test figure the permissible value of the error between (i=0,1,2, L, N-1)

(5) the reference value I of aqueous phase mobility, non-Darcy flow β-factor and acceleration factor _1r, β _1rwith power exponent n _{1 β}and n _1c; The reference value I of husky phase mobility, non-Darcy flow β-factor and acceleration factor _2r, β _2rwith power exponent n _{2 β}and n _2c.

In this example, these parameters are respectively

(1) the mass density ρ=1000 (kg/m of liquid ³), kinetic viscosity μ=1.01 × 10 ^-3(Pas)

(2) rock sample radius r _s=0.50 × 10 ^-2(m), naturally stacking height h ₀=1.25 × 10 ^-1height h=1.00 × 10 after (m), compression ^-1(m), fractured rock mass M _rock=1500 (kg), fine sand quality is M _sand=430 (kg);

(3) sampling period τ=1.0 (s), data length N=301;

(4) limits of error ε of Newton tangent method ₁=1.0 × 10 ^-9, Runge-Kutta method limits of error ε=1.0 × 10 ^-4, the numerical solution of percolation flow velocity (i=0,1,2, L, N-1) and test figure the permissible value of the error between (i=0,1,2, L, N-1) .

The limits of error ε of Newton tangent method ₁, Runge-Kutta method limits of error ε,

(5) utilize genetic algorithm, through several times breeding, be met the reference value I of the aqueous phase mobility of error requirements, non-Darcy flow β-factor and acceleration factor _1r, β _1rwith , power exponent n _{1 β}and n _1c, the reference value I of husky effective mobility, non-Darcy flow β-factor and acceleration factor mutually _2r, β _2rwith , power exponent n _{2 β}and n _2c, be respectively

I _1r＝0.1×10 ^-13，β _1r＝0.22×10 ⁺¹⁰， ，n _1β＝-1.6，n _1c＝-1.4；

I _2r＝0.869×10 ^-10，β _2r＝0.653×10 ⁺⁰⁷， ，n _b＝-1.86，n _c＝-1.31。

3 perviousness Parameters Calculation

Utilize the present invention, calculate aqueous phase and husky phase perviousness parameter, see Fig. 2 and Fig. 3 respectively.

4 percolation flow velocities and calculated value compare with test value

Fig. 4 gives the reference value I of aqueous phase mobility, non-Darcy flow β-factor and acceleration factor _1r, β _1rwith , power exponent n _{1 β}and n _1cget the time series of the water percolation flow velocity of two groups of different numerical value.I in Fig. 4 (a) _1r=0.821 × 10 ^-11, β _1r=0.241 × 10 ^{+ 11}, , n _β=-1.3, n _c=-1.7, there is appreciable error in the computable value with test value of water percolation flow velocity; In Fig. 4 (a), I _1r=0.1 × 10 ^-13, β _1r=0.22 × 10 ^{+ 10}, , n _β=-1.6, n _c=-1.4, the computable value with test value error of water percolation flow velocity is minimum.

Claims (7)

1. water husky two-phase flow perviousness Parameters Calculation method in water drive sand test, namely the mobility of aqueous phase and husky phase, non-Darcy flow β-factor and acceleration factor is calculated according to pressure gradient time series, discharge time series and sand drift amount time series, it is characterized in that, it is realized by following steps:

Step 1) to gather discharge, sand drift amount and pressure differential at equal intervals, calculate the water percolation flow velocity of each sampling instant, husky percolation flow velocity and pressure gradient, obtain water percolation flow velocity time series, husky percolation flow velocity time series and pressure gradient time series; According to kinematic relation, and replace difference quotient by difference coefficient, set up the time series of aqueous phase and husky phase percolation flow velocity local derivative;

Step 2) volume fraction of each sampling instant aqueous phase and husky phase is calculated according to permeameter size, rock sample quality, sandy amount etc.;

Step 3) utilize the momentum conservation equation of aqueous phase, husky phase, and utilize mobility, power exponent relation between non-Darcy flow β-factor and acceleration factor, set up the algebraic equation of each sampling instant current degree and sand drift degree respectively; Utilize Newton tangent method Solving Algebraic Equation, obtain the time series of current degree and sand drift degree;

Step 4) utilize mobility, power exponent relation between non-Darcy flow β-factor and acceleration factor, calculate the non-Darcy flow β-factor of aqueous phase and husky phase and the time series of acceleration factor;

Step 5) utilize four step Runge-Kutta respectively to the numerical solution of the momentum conservation equation of aqueous phase, husky phase, obtain the time series of aqueous phase and husky phase percolation flow velocity;

Step 6) compare aqueous phase and the husky calculated value of phase percolation flow velocity and the error of measured value, if each sampling instant percolation flow velocity error sum of squares is greater than certain preset value, parameter in adjustment mobility, non-Darcy flow β-factor and acceleration factor relational expression, recalculates the time series of the mobility of aqueous phase and husky phase, non-Darcy flow β-factor and acceleration factor; If each sampling instant percolation flow velocity error sum of squares is less than or equal to certain preset value, then using the test value of the result of calculation of mobility, non-Darcy flow β-factor and acceleration factor as the perviousness parameter of the husky two-phase flow of fractured rock water.

2. water husky two-phase flow perviousness Parameters Calculation method in a kind of water drive sand according to claim 1 test, it is characterized in that, the concrete steps of step 1 are:

Step 1.1) gather discharge, sand drift amount and pressure differential with τ at equal intervals, pressure gradient G can be obtained respectively by simple algebraic operation _p, water percolation flow velocity V ₁with husky percolation flow velocity V ₂time series, namely

$G_p^{\left(i\right)}=G_p{|}_{t=\mathrm{i\τ}},(i=\mathrm{0,1,2},\·\·\·N)---\left(18\right)$

$V_1^{\left(i\right)}=V_1{|}_{t=\mathrm{i\τ}},(i=\mathrm{0,1,2},\·\·\·N)---\left(19\right)$

$V_2^{\left(i\right)}=V_2{|}_{t=\mathrm{i\τ}},(i=\mathrm{0,1,2},\·\·\·N)---\left(20\right)$

Step 1.2) replace difference quotient by difference coefficient, calculate percolation flow velocity V ₁and V ₂local derivative with time series, namely

${\mathrm{\α}}_1^{t_i}=\frac{{\∂V}_1}{\∂t}=\frac{V_1^{t_{i+1}}-V_1^{t_i}}{\mathrm{\τ}},i=\mathrm{1,2},\·\·\·,N-1---\left(21\right)$

With

${\mathrm{\α}}_2^{t_i}=\frac{{\∂V}_2}{\∂t}=\frac{V_2^{t_{i+1}}-V_2^{t_i}}{\mathrm{\τ}},i=\mathrm{1,2},\·\·\·,N-1---\left(22\right)$

3. water husky two-phase flow perviousness Parameters Calculation method in a kind of water drive sand according to claim 1 test, it is characterized in that, the computing method of the husky volume fraction initial value described in step 2 are:

Quality is M by the first step _rockbroken rock sample and quality be M _sandfine sand blending evenly as the sample permeating rock sample;

Second step, puts into the permeameter of cylinder barrel diameter 2a by sample, and measures the height of the stacking naturally h of broken rock sample fine sand admixture _uncomp;

3rd step, utilizes Material Testing Machine to load admixture, when sample height reaches preset value h ₀in time, stops loading, and keeps rock sample constant height;

4th step, is calculated as follows fractured rock factor of porosity

$\mathrm{\φ}=1-\frac{M_{\mathrm{rock}}}{m_{\mathrm{rock}}{\mathrm{\π\α}}^2{h}_0}---\left(15\right)$

Wherein, m _rockfor broken rock sample parent, the i.e. mass density of rock block;

5th step, calculates the volume fraction before infiltration.Before infiltration in fractured rock hole, the volume fraction of air is calculated as follows

$Y_1^0=1-\frac{\frac{M_{\mathrm{sand}}}{m_2}}{{\mathrm{\π\α}}^2{h}_0-\frac{M_{\mathrm{rock}}}{m_{\mathrm{rock}}}}---\left(16\right)$

Husky volume fraction is

$Y_2^0=\frac{\frac{M_{\mathrm{sand}}}{m_2}}{{\mathrm{\π\α}}^2{h}_0-\frac{M_{\mathrm{rock}}}{m_{\mathrm{rock}}}}---\left(17\right)$

6th step, calculates initial volume mark

Before seepage flow starts, carry out water saturation to rock sample, water drive has replaced air; If water injection pressure is very little, tiny loss is little, can think the initial volume mark of water equal the volume fraction of fine sand equal

4. water husky two-phase flow perviousness Parameters Calculation method in a kind of water drive sand according to claim 1 test, it is characterized in that, the concrete steps of step 3 are:

Step 3.1) relational expression between setting perviousness parameter

The perviousness parameter of aqueous phase and husky phase medium is set as power exponent relation, namely all has in each sampling instant

$\left\{\begin{array}{c}\frac{{\mathrm{\β}}_1^{t_i}}{{\mathrm{\β}}_{1r}}={\left(\frac{I_1^{t_i}}{I_{1r}}\right)}^{-n_{1\mathrm{\β}}}\\ \frac{c_{1\mathrm{\α}}^{t_i}}{c_{1\mathrm{\α}}^r}={\left(\frac{I_1^{t_i}}{I_{1r}}\right)}^{-n_{1c}}\end{array}\right.,i=\mathrm{0,1,2},\·\·\·,N-1---\left(23\right)$

With

$\left\{\begin{array}{c}\frac{{\mathrm{\β}}_2^{t_i}}{{\mathrm{\β}}_{2r}}={\left(\frac{I_{2e}^{t_i}}{I_{2r}}\right)}^{-n_{2\mathrm{\β}}}\\ \frac{c_{2\mathrm{\α}}^{t_i}}{c_{2\mathrm{\α}}^r}={\left(\frac{I_{2e}^{t_i}}{I_{2r}}\right)}^{-n_{2c}}\end{array}\right.,i=\mathrm{0,1,2},\·\·\·,N-1---\left(24\right)$

Wherein, the reference value of aqueous phase mobility, non-Darcy flow β-factor and acceleration factor is I _1r, β _1rwith power exponent n _{1 β}and n _1c, husky phase mobility, non-Darcy flow β-factor and acceleration factor reference value I _2r, β _2rwith power exponent n _{2 β}and n _2c;

Step 3.2) set up the closed system of equations of mobility

Power exponent relation between the perviousness parameter of aqueous phase and husky phase medium and momentum conservation equation simultaneous, by formula (23), formula (24) and respectively formula (13), formula (14) simultaneous, obtain the system of equations that aqueous phase and husky phase mobility meet respectively

$m_1{c}_{1\mathrm{\α}}^r{\left(\frac{I_1^{t_i}}{I_r}\right)}^{-n_{1c}}{\mathrm{\α}}_1^{\left(i\right)}={-G}_p^{t_i}-\frac{1}{I_1^{t_i}}{V}_1^{t_i}-m_1{\mathrm{\β}}_{1r}^{t_i}{\left(\frac{I_1^{t_i}}{I_{1r}}\right)}^{-n_{1\mathrm{\β}}}\left|V_1^{t_i}\right|V_1^{t_i},i=\mathrm{0,1,2},\·\·\·,N-1---\left(25\right)$

With

$m_2{c}_{2\mathrm{\α}}^r{\left(\frac{I_{2e}^{t_i}}{I_{2r}}\right)}^{-n_{2c}}{\mathrm{\α}}_2^{t_i}={-G}_p^{t_i}-\frac{1}{I_{2e}^{t_i}}{\left(V_2^{t_i}\right)}^n-m_2{\mathrm{\β}}_{2r}^{t_i}{\left(\frac{I_{2e}^{t_i}}{I_{2r}}\right)}^{-n_{2\mathrm{\β}}}\left|V_2^{t_i}\right|V_2^{t_i},i=\mathrm{0,1,2},\·\·\·,N-1---\left(26\right)$

Step 3.3) primary Calculation sampling instant aqueous phase and husky phase mobility

Step 3.3.1), design iteration form, introduces function

$f_{1i}\left(I_1^{t_i}\right)=m_1{\mathrm{\β}}_{1r}^{t_i}{\left(\frac{I_1^{t_i}}{I_{1r}}\right)}^{-n_{1\mathrm{\β}}}\left|V_1^{t_i}\right|V_1^{t_i}+\frac{1}{I_1^{t_i}}{V}_1^{t_i}+m_1{c}_{1\mathrm{\α}}^r{\left(\frac{I_1^{t_i}}{I_r}\right)}^{-n_{1c}}{\mathrm{\α}}_1^{\left(i\right)}+G_p^{t_i},i=\mathrm{0,1,2},\·\·\·N-1---\left(27\right)$

$f_{2i}\left(I_{2e}^{t_i}\right)=m_2{\mathrm{\β}}_{2r}^{t_i}{\left(\frac{I_{2e}^{t_i}}{I_{2r}}\right)}^{-n_{2\mathrm{\β}}}\left|V_2^{t_i}\right|V_2^{t_i}+\frac{1}{I_{2e}^{t_i}}{V}_2^{t_i}+m_2{c}_{2\mathrm{\α}}^r{\left(\frac{I_{2e}^{t_i}}{I_{2r}}\right)}^{-n_{2c}}{\mathrm{\α}}_2^{t_i}+G_p^{t_i},i=\mathrm{0,1,2},\·\·\·N-1---\left(28\right)$

System of equations (27) and (28) are deformed into respectively

$f_{1i}\left(I_1^{t_i}\right)=0,i=\mathrm{0,1,2},\·\·\·,N-1---\left(29\right)$

$f_{2i}\left(I_{2e}^{t_i}\right)=0,i=\mathrm{0,1,2},\·\·\·,N-1---\left(30\right)$

Utilize the Iteration of Newton tangent method equationof structure group (29) and (30);

Step 3.3.2) determine iteration initial value with

Get the stable percolation equation root of Newton fluid as iteration initial value, namely

$-G^{t_i}-\frac{1}{I_1^{*}}{V}_1^{t_i}-m_1{\mathrm{\β}}_1^{t_i}\left|V_1^{t_i}\right|V_1^{t_i}=0,i=\mathrm{1,2},\·\·\·,N-1---\left(31\right)$

$-G_p^{t_i}-\frac{1}{I_{2e}^{*}}{V}_2^{t_i}-m_2{\mathrm{\β}}_2^{t_i}\left|V_2^{t_i}\right|V_2^{t_i}=0,i=\mathrm{1,2},\·\·\·,N-1---\left(32\right)$

Step 3.3.3) according to the Iteration of first two steps and iteration initial value, calculating sampling moment aqueous phase and husky phase mobility with

5. water husky two-phase flow perviousness Parameters Calculation method in a kind of water drive sand according to claim 1 test, it is characterized in that, the concrete steps of step 4 are:

Step 4.1) according to formula (23) and formula (24), calculate aqueous phase and husky phase non-Darcy flow β-factor and acceleration factor time series, ${\mathrm{\β}}_1^{t_i},c_{\mathrm{\α}1}^{t_i},{\mathrm{\β}}_2^{t_i},c_{\mathrm{\α}2}^{t_i}(i=\mathrm{1,2},...,N-1).$

6. water husky two-phase flow perviousness Parameters Calculation method in a kind of water drive sand according to claim 1 test, it is characterized in that, the concrete steps of step 5 are:

Step 5.1) calculate aqueous phase and husky phase percolation flow velocity

Step 5.1.1) build the external function of momentum conservation equation, by momentum conservation equation (13) and (14) distortion, obtain

$\frac{{\∂V}_1}{\∂t}=-\frac{1}{m_1{c}_{1\mathrm{\α}}}(G_p+\frac{1}{I_1}{V}_1+m_1{\mathrm{\β}}_1|V_1\left|V_1\right)---\left(33\right)$

$\frac{{\∂V}_2}{\∂t}=-\frac{1}{m_2{c}_{2\mathrm{\α}}}(G_p+\frac{1}{I_{2e}}{V}_2^n+m_2{\mathrm{\β}}_2|V_2\left|V_2\right)---\left(34\right)$

Respectively to pressure gradient G _p, mobility I ₁and I _2e, non-Darcy flow β-factor β ₁and β ₂, acceleration factor c _1aand c _2acarry out three knot interpolations between the whole district, the external function of equationof structure (33) and (34), for any time t ∈ (t _k, t _k+2), have

$G_p=G_p^{t_k}\frac{(t-t_{k+1})(t-t_{k+2})}{(t_k-t_{k+1})(t_k-t_{k+2})}+G_p^{t_{k+1}}\frac{(t-t_k)(t-t_{k+2})}{(t_{k+1}-t_k)(t_{k+1}-t_{k+2})}+G_p^{t_{k+2}}\frac{(t-t_k)(t-t_{k+1})}{(t_{k+2}-t_k)(t_{k+2}-t_{k+1})}---\left(35\right)$

$I_1=I_1^{t_k}\frac{(t-t_{k+1})(t-t_{k+2})}{(t_k-t_{k+1})(t_k-t_{k+2})}+I_1^{t_{k+1}}\frac{(t-t_k)(t-t_{k+2})}{(t_{k+1}-t_k)(t_{k+1}-t_{k+2})}+I_1^{t_{k+2}}\frac{(t-t_k)(t-t_{k+1})}{(t_{k+2}-t_k)(t_{k+2}-t_{k+1})}---\left(36\right)$

$I_{2e}=I_{2e}^{t_k}\frac{(t-t_{k+1})(t-t_{k+2})}{(t_k-t_{k+1})(t_k-t_{k+2})}+I_{2e}^{t_{k+1}}\frac{(t-t_k)(t-t_{k+2})}{(t_{k+1}-t_k)(t_{k+1}-t_{k+2})}+I_{2e}^{t_{k+2}}\frac{(t-t_k)(t-t_{k+1})}{(t_{k+2}-t_k)(t_{k+2}-t_{k+1})}---\left(37\right)$

${\mathrm{\β}}_1={\mathrm{\β}}_1^{t_k}\frac{(t-t_{k+1})(t-t_{k+2})}{(t_k-t_{k+1})(t_k-t_{k+2})}+{\mathrm{\β}}_1^{t_{k+1}}\frac{(t-t_k)(t-t_{k+2})}{(t_{k+1}-t_k)(t_{k+1}-t_{k+2})}+{\mathrm{\β}}_1^{t_{k+2}}\frac{(t-t_k)(t-t_{k+1})}{(t_{k+2}-t_k)(t_{k+2}-t_{k+1})}---\left(38\right)$

${\mathrm{\β}}_2={\mathrm{\β}}_2^{t_k}\frac{(t-t_{k+1})(t-t_{k+2})}{(t_k-t_{k+1})(t_k-t_{k+2})}+{\mathrm{\β}}_2^{t_{k+1}}\frac{(t-t_k)(t-t_{k+2})}{(t_{k+1}-t_k)(t_{k+1}-t_{k+2})}+{\mathrm{\β}}_2^{t_{k+2}}\frac{(t-t_k)(t-t_{k+1})}{(t_{k+2}-t_k)(t_{k+2}-t_{k+1})}---\left(39\right)$

$c_{\mathrm{\α}1}=c_{\mathrm{\α}1}^{t_k}\frac{(t-t_{k+1})(t-t_{k+2})}{(t_k-t_{k+1})(t_k-t_{k+2})}+c_{\mathrm{\α}1}^{t_{k+1}}\frac{(t-t_k)(t-t_{k+2})}{(t_{k+1}-t_k)(t_{k+1}-t_{k+2})}+c_{\mathrm{\α}1}^{t_{k+2}}\frac{(t-t_k)(t-t_{k+1})}{(t_{k+2}-t_k)(t_{k+2}-t_{k+1})}---\left(40\right)$

$c_{\mathrm{\α}2}=c_{\mathrm{\α}2}^{t_k}\frac{(t-t_{k+1})(t-t_{k+2})}{(t_k-t_{k+1})(t_k-t_{k+2})}+c_{\mathrm{\α}2}^{t_{k+1}}\frac{(t-t_k)(t-t_{k+2})}{(t_{k+1}-t_k)(t_{k+1}-t_{k+2})}+c_{\mathrm{\α}2}^{t_{k+2}}\frac{(t-t_k)(t-t_{k+1})}{(t_{k+2}-t_k)(t_{k+2}-t_{k+1})}---\left(41\right)$

Step 5.1.2) be step-length with τ, utilize the four step Runge-Kutta of variable step to ask the numerical solution of equation (33) and (34), obtain percolation flow velocity calculated value time series with

7. water husky two-phase flow perviousness Parameters Calculation method in a kind of water drive sand according to claim 1 test, it is characterized in that, the concrete steps of step 6 are:

Calculate the numerical solution of percolation flow velocity ${\tilde{V}}_1^{t_i},{\tilde{V}}_2^{t_i}(i=\mathrm{0,1,2},L,N-1)$ With test figure

Difference between (i=0,1,2, L, N-1)

$E_{\mathrm{rr}}=\max \{\frac{1}{N-1}\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{(1-\frac{{{\tilde{V}}_1}^{t_i}}{{V_1}^{t_i}})}^2,\frac{1}{N-1}\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{(1-\frac{{\tilde{V}}_2^{t_i}}{V_2^{t_i}})}^2\}---\left(42\right)$

If E _rrbe greater than the permissible value of a certain prior setting then adjust the reference value I of aqueous phase mobility, non-Darcy flow β-factor and acceleration factor _1r, β _1rwith power exponent n _{1 β}and n _1c, husky phase mobility, non-Darcy flow β-factor and acceleration factor reference value I _2r, β _2rwith power exponent n _{2 β}and n _2cvalue, recalculate mobility non-Darcy flow β-factor, acceleration factor and percolation flow velocity until E _rrbe not more than