CN104331579A - Simulation method of low-permeability reservoir crude oil boundary layer - Google Patents

Abstract

The invention discloses a simulation method of a low-permeability reservoir crude oil boundary layer. The simulation method comprises the following steps: (1) building a digital rock core by a CT (Computed Tomography) scanning method; (2) dividing the built digital rock core by a Delauny triangular dividing algorithm to obtain a three-dimensional geometric model; and (3) carrying out numerical simulation of the divided three-dimensional description geometric model of the fluid by using a finite element method under the condition that the boundary layer influence is considered, namely solving an incompressible equation which is N-S equation of the fluid with micro flowing in holes to obtain a flowing rate field and a pressure distribution field of the fluid. According to the simulation method of the low-permeability reservoir crude oil boundary layer, comsol simulation is used; the manpower, material resources and time can be saved; a regular hole mesh model is used, the pore throat condition of the real rock core can be accurately simulated; the condition is more accordance with the actual condition; data results are processed, the thickness of the boundary layer can be accurately calculated and is closer to the thickness of the crude oil boundary layer of the real rock core.

Description

A kind of analogy method of low-permeability oil reservoirs boundary layer

Technical field

The present invention relates to a kind of analogy method of low-permeability oil reservoirs boundary layer, belong to the technical field of the simulated experiment of natural petroleum gas field.

Background technology

Current research shows: existing recovery mechanism research majority is all based on macroscopical seepage theory, macroscopic view seepage theory is all based on darcy flow equation, can only the seepage flow characteristics of macroscopic description reservoir fluid so apply the percolation equationk set up based on Darcy's equation in field of oil development, the description in more accurate hole rank cannot be provided.Particularly in the research of low permeability reservoir, many scholars are shown by flow in low permeability core displacement test and Development Practice: in low permeability reservoir, and oil and gas flow does not meet Darcy's law.For characterizing the seepage flow characteristics of low permeability reservoir, many scholars propose various mathematical model, wherein mainly contain quasi-threshold pressure gradient model, segmented model, continuous model.Above-mentioned mathematical model is all the Function Fitting based on experimental data, is short of all to some extent the explanation of seepage flow characteristics.Most macroscopic properties of reservoir, as in fact capillary pressure, permeability, relative permeability etc. depend on the micromechanism of rock.Therefore, just rest on things surface only by the research launched the analysis of macroscopic property, cannot be familiar with to some extent the essence of its inside, the technological means of the raising recovery ratio proposed thus also has its limitation in essence.

If make raising oil recovery factor technology increase substantially, only on macro-level, carry out theoretical research and technological development is far from being enough, must break through from microcosmic point, the rational microscopic seepage of Erecting and improving is theoretical.The essential problem of macroscopical percolation phenomenon can be determined by the space distribution, the characteristic distributions of fluid in space, duct, interaction mechanism etc. studying porous medium internal void space, fundamentally could being familiar with the contact of Macrocosm and microcosm, finding better technological means for improving recovery ratio.Digital cores technology carries out microscopic seepage to study necessary platform, is the key point solved the problem.The pore texture of porous medium can be fully realized, the essential problems such as the rule that fluid flows in porous medium and interaction mechanism by digital cores technology.The simultaneously utilization of CT technology in digital cores builds is improving and strengthening digital cores technology, and the digital cores set up by CT technology is had closest to the feature such as true.

For middle low permeability reservoir, the principal element controlling rock sample percolation ability is throat character instead of pore character, therefore, is deeply familiar with the micropore structure of low-permeability oil deposit, is the key improving low permeability reservoir use rate.

The feature of tight sand medium: Compacted rock, pore-size is tiny, and permeability is low; Pore constriction complex geometry, capillarity is serious; Fluid and oil reservoir medium mechanism of action complexity, non newtonian phenomenon is relatively remarkable.

At present, tight sandstone reservoir has become the important object of oil-gas mining, but does not meet darcy rule due to the seepage flow characteristics (apparent permeability etc. of free-boundary problem, change) of wherein fluid, and a lot of microscopic mechanism problem not yet solves.Therefore, by digital cores technology, reservoir micropore structure to be evaluated, and simulating boundary effect layer, and then analyze the impact on non-darcy flow curve and oil-water two-phase flow curve of micropore structure and effect of boundary layer.

At present, domestic most employing experimental technique Measured Boundary layer thickness, mainly comprises and measures the effective radius of glass capillary and try to achieve mean boundary-layer thickness and measure rock capillary model and ask for mean boundary-layer thickness.The weak point of these two kinds of models is respectively, and glass capillary model method is based on Darcy's law, and for hyposmosis, extra-low permeability oil reservoirs, the flowing of fluid in rock core pore throat does not meet Darcy's law; Represent rock core pore throat with glass capillary simultaneously, can not simulate when fluid flows in pore throat and the interaction of core surface.Rock capillary model is poor according to the crude quality before and after rock core displacement, calculates the average thickness in boundary layer, and due to the impact of experiment condition, error is comparatively large, and model is too idealized.

Summary of the invention

For the deficiencies in the prior art, the invention provides a kind of analogy method of low-permeability oil reservoirs boundary layer.This analogy method can simulate the crude oil of different viscosities, under different pressure gradients, flows through the situation of different pore throat radius, can calculate the net thickness of Crude Oil in Boundary Layer by data processing.

Technical scheme of the present invention is as follows:

A kind of low-permeability oil reservoirs boundary layer simulating method, comprises step as follows:

1) CT scan method structure digital cores comprises step (a)-(c):

A () obtains the multi-faceted CT projecting image data of rock sample by CT scan; Wherein the shape of rock sample does not affect CT projected image, and the CT that the present invention uses is MicroXCT-400 high-resolution 3DX ray microscope, can carry out 360 degree of imagings to rock sample, so the angle of shooting rock sample does not affect for rock sample projected image.In theory, the MicroXCT-400CT scanning device that the present invention uses supports that sample size reaches 300mm, and weight reaches 15Kg; Utilize CT to scan, individual CT projecting image data after scanning is automatically caught by Image Acquisition software and stores; By described rock sample after its axial-rotation fixed angle again to its CT scan, individual CT projecting image data after scanning is automatically caught by Image Acquisition software and stores, until described rock sample terminates to gather CT projecting image data after adding up rotating 360 degrees vertically;

The image-forming principle of CT: when X ray is through any material, its can with the atomic interaction of material and cause energy attenuation, and the different constituent of material has different absorption coefficients (i.e. attenuation coefficient) to X ray, say it on the contrary, the constituent of material can be judged by measurement of species to the absorption coefficient of X ray.When beam of x-rays is through object, it penetrate all substances on path to the summation of X-ray absorption coefficient all by the measurement result that is reflected in X-ray intensity, as shown in (1.1):

$I=I_o{e}^{-\underset{i}{\mathrm{\Σ}}{\mathrm{\μ}}_i{x}_i}---\left(1.1\right)$

In formula, I _ofor the initial strength of X ray, I is the intensity after X ray passes object, namely X ray be attenuated after intensity, i represents a certain constituent on the path passed at this X ray in material, μ _i, x _ibe respectively the i-th component to the attenuation coefficient of X ray and this component in the current length on path of X ray;

The image-forming principle of CT is set up on the basis of this just, by measuring the X ray through object section, afterwards, adopt certain method for reconstructing to calculate and object tomography locus absorption coefficient one to one, thus recover the structural information of object section;

During analog computation of the present invention, the concrete size of rock sample will be determined according to the attribute of rock sample, to reach best technique effect.For ensureing that rock core CT image can clear reflection rock core internal microstructure feature, experiment rock core size can not be excessive, reason mainly contains two aspects: on the one hand, because rock core density is larger, ray penetrate in rock core process be attenuated to relatively large extent cause its arrive detector system time signal very faint, if rock core is oversize, the complete conductively-closed of ray may be caused; On the other hand, the maximum lattice point quantity adopted when solving due to software systems is certain, and rock core size increases the physical size increase meaning that each lattice point is corresponding, and namely the resolution of image reduces.For existing industrial CT equipment, obtain CT image more clearly, rock core cylinder largest diameter is centimetre-sized;

B () utilizes image rebuilding method that described CT projecting image data is converted to the gray level image of described rock sample;

The core of CT imaging is how by reconstructs projection data image.The essence of image reconstruction is the absorption coefficient being solved each pixel in image array by the CT data for projection after gathering, and then re-constructs the process of image; The method for solving of image reconstruction problem can be divided into two large classes according to its feature: the first kind is converter technique (being also analytical method), is characterized in that last discretize utilizes computer calculate first in continuous domain dissection process; Equations of The Second Kind is process of iteration (being also algebraic reconstruction technique, Series Expansion Method, optimisation technique etc.), is characterized in that directly carrying out discretize analysis obtains numerical solution.Comparatively speaking, converter technique advantages: realize simple, speed is fast, can obtain good reconstruction quality to enough accurate data for projection.Therefore, extensively converter technique is adopted in current practical CT system.

Described image rebuilding method is converter technique, comprise not containing the Fourier method for reconstructing FRM of back projection's step with based on the convolution of back projection or filtering method: the former with direct Fourier transform method DFM for representative, the latter with convolution or filtered back projection C/FBP method for representative;

The basic thought of DFM is according to the mechanism of action between radiation source and imageable target, search out objective function f (x, y) relation between two-dimensional fourier transform (FT) and the one dimension FT of projection function, utilize this relation, although objective function f is (x, y) be unknown, but its two-dimentional FT obtains by the one dimension FT of projection function in different azimuth, once obtain f (x, y) all projections, also two-dimentional FT, the f (x, y) that just obtain it then can be brought by the contravariant of two-dimentional FT to be obtained;

Come from Radon based on the convolution of back projection or filtering method C/FBP method to negate the direct realization of formula; C/FBP formula utilizes FT and projection theorem to derive; If the most high spatial frequency of p (x', θ) is B;

The mathematical formulae of C/FBP is write as

$f(x,y)={\∫}_0^{\mathrm{\π}}q(x^{\′},\mathrm{\θ})\mathrm{d\θ}---\left(1.2\right)$

Wherein,

$q(x^{\′},\mathrm{\θ})={\∫}_{-\∞}^{\∞}\left|R\right|P(R,\mathrm{\θ})\exp \left(j2\mathrm{\π}{x}^{\′}R\right)\mathrm{dR}---\left(1.3\right)$

Or

q(x',θ)＝c(x')*p(x',θ) (1.4)

In formula (1.2), (1.3), (1.4), the objective function that f (x, y) is image reconstruction, namely rock sample is to the distribution function of gamma ray absorption coefficient, and x, y are respectively the rectangular coordinate of certain point on rock sample scanning section; P (R, θ), p (x', θ) be f (x, y) along the projection function on θ direction, i.e. the transmitted intensity signal that detects of CT detector, x' is certain point (x on rock sample scanning section, y) horizontal ordinate in the rotated coordinate system, R for this point is to rotation center, i.e. true origin, distance; C (x') is convolution function, and it is by wave filter | the Fourier inverse transformation of R| (| R|≤B) is tried to achieve;

Q (x', θ) in formula (1.3), (1.4) is called projection P (R, θ), the filter projection of p (x', θ) and convolution projection; The image rebuilding method realized by formula (1.2), (1.3) is called filtered back-projection method FBP, and the image rebuilding method realized by formula (1.2), (1.4) is called the convolution back projection method CBP; Therefore, as obtained the projection of image, then oppositely solved by convolution or filtered back-projection method and obtain the distribution function f (x, y) of rock sample section to gamma ray absorption coefficient, the scan image of rock sample section is just obtained after adopting the mode of gray level image to express f (x, y);

C () adopts image binary segmentation method to be separated space, gray level image mesoporosity and the rock sample skeleton of described rock sample, set up digital cores;

In gray level image, pore space and rock skeleton are divided it by image partition method: Iamge Segmentation is choosing of segmentation threshold, if threshold value is chosen unreasonable, be easy to destination object to understand become background image, or background is interpreted as destination object, so just cannot effective evaluating objects feature when carrying out signature analysis to target image;

Threshold definitions is as follows: suppose that F (i, j) represents image pixel gray level value in gray level image, and the pixel that gray-scale value is greater than threshold value T after being separated represents destination object, i.e. hole, represents with 1; Other pixel represents image background, i.e. skeleton, represents with 0, then can obtain destination object by following formula:

$g(i,j)=\left\{\begin{array}{cc}1;& F(i,j)\≥T\\ 0;& \mathrm{else}\end{array}\right.---\left(1.5\right)$

In formula, T is threshold value, also known as threshold value;

Based on the image binary segmentation method of core porosity

The fundamental physical quantity that core porosity must measure when being and carrying out reservoir physics experiment; Under existing experiment condition, measurement accurately can be facilitated to obtain by gas expansion method, hold-up method etc.; Because factor of porosity gives the share shared by space, rock core mesoporosity, therefore, in known cases, factor of porosity is attached to the segmentation quality that can improve image in the cutting procedure of rock core CT image undoubtedly as constraint condition;

If the factor of porosity of rock sample is maximum, the minimum gradation value of image are respectively I _max, I _mIN, gray-scale value is the pixel count of i is p (i), then meet the gray-scale value k of formula (1.6) ^*i.e. required segmentation threshold:

Image binary segmentation method based on rock sample factor of porosity is actually selects suitable segmentation threshold, and the ratio making the core image mesoporosity pixel count obtained by this Threshold segmentation account for total pixel number equals core porosity.Because grey scale pixel value is by the integer representation in 0 to 255 intervals, therefore, splitting the image factor of porosity obtained may exist nuance with core experiment factor of porosity, but this does not affect Iamge Segmentation quality.

Find by carrying out segmentation test to experimental image, image binary segmentation method based on core porosity can obtain high-quality segmentation result, therefore, when core experiment factor of porosity is known, binary segmentation can be carried out by the method to rock core CT image.But, and the rock core that not all is used for CT experimental study has experimental port porosity, when not having factor of porosity data, need split rock core CT image by other image binary segmentation method.In image binary segmentation method, maximum kind spacing method is because principle is simple, program design facilitates and can obtain good segmentation effect and be widely used.

When not having factor of porosity, maximum kind spacing method is utilized to split gray level image; What formed by the gray-scale value of pixels all in gray level image to be split is totally considered as a set, assuming that the set threshold value of gray-scale value is divided into two groups, then when the variance (between-group variance) of mean value of two groups and the ratio of the variance (interclass variance) of each group reach maximum, the segmentation threshold corresponding to it is reasonable value; The pore space of described rock sample inside respectively becomes one group with skeleton, and the difference of each group inside is very little and nature difference between two groups is very large; After two components cut open by the rational threshold value of employing, nature difference between two groups will reach maximum value with the ratio of difference in group, visible, feature and the segmentation thinking of maximum kind spacing method to image of rock sample CT scan projected image itself are very identical, and this is also that maximum kind spacing method segmentation rock core CT image can obtain good segmentation result and the major reason that is widely adopted;

The principle of maximum kind spacing method selected threshold is as follows.If maximum, the minimum gradation value of image are respectively I _max, I _mIN, gray-scale value is the pixel count of i is p (i); Gray threshold is k, and pixel is divided into two groups by it: group 1 (being assumed to be pore space) and group 2 (being assumed to be rock skeleton).Add up respectively and calculate the pixel count N of two groups ₁(k), N ₂(k), average gray value μ ₁(k), μ ₂(k) and variance and the average gray value μ of whole pixel _t, then interclass variance and between-group variance be respectively:

${\mathrm{\σ}}_W^2=N_1{\mathrm{\σ}}_1^2+N_2{\mathrm{\σ}}_2^2---\left(1.7\right)$

${\mathrm{\σ}}_B^2=N_1{({\mathrm{\μ}}_1-{\mathrm{\μ}}_T)}^2+N_2{({\mathrm{\μ}}_2-{\mathrm{\μ}}_T)}^2---\left(1.8\right)$

In formula:

${\mathrm{\μ}}_1\left(k\right)=\frac{\underset{i=I_{\mathrm{MIN}}}{\overset{k}{\mathrm{\Σ}}}i\·p\left(i\right)}{N_1\left(k\right)}$ ${\mathrm{\μ}}_2\left(k\right)=\frac{\underset{i=k+1}{\overset{I_{\mathrm{MAX}}}{\mathrm{\Σ}}}i\·p\left(i\right)}{N_2\left(k\right)}$

${\mathrm{\μ}}_T\left(k\right)=\frac{\underset{i=I_{\mathrm{MIN}}}{\overset{I_{\mathrm{MAX}}}{\mathrm{\Σ}}}i\·p\left(i\right)}{N_1\left(k\right)+N_2\left(k\right)}$

${\mathrm{\σ}}_1^2\left(k\right)=\frac{\underset{i=I_{\mathrm{MIN}}}{\overset{k}{\mathrm{\Σ}}}{[i-{\mathrm{\μ}}_1\left(k\right)]}^2\·p\left(i\right)}{N_1\left(k\right)}$ ${\mathrm{\σ}}_2^2\left(k\right)=\frac{\underset{i=k+1}{\overset{I_{\mathrm{MAX}}}{\mathrm{\Σ}}}{[i-{\mathrm{\μ}}_2\left(k\right)]}^2\·p\left(i\right)}{N_2\left(k\right)}$

Select optimum gradation value k ^*, make it meet (1.9) formula, then k ^*i.e. required threshold value.

$f\left(k^{*}\right)=\max \{f\left(k\right)={\mathrm{\σ}}_B^2\left(k\right)/{\mathrm{\σ}}_W^2\left(k\right)\}---\left(1.9\right);$

Adopt the algorithm flow of maximum kind spacing method segmentation image, as shown in Figure 1;

2) use Delaunay triangulation to carry out subdivision to the above-mentioned digital cores built, obtain 3-D geometric model;

Wherein, described Delaunay triangulation is also a kind of iso-surface patch method, its advantage be in the process building geometric model by controlling the parameter of geometric units, make the geometric model formed be more conducive to after analysis and calculation.Due to the plurality of advantages of Delaunay triangulation, this method has been widely used in multiple field such as analysis, medical science of finite element ^[1].[1] Ding Yongxiang, the summer huge Zan. the Delaunay triangulation [J] of arbitrary polygon. Chinese journal of computers, 1994,17 (4): 270-275.

Delaunay triangulation is the triangulation optimized, and has a lot of good characteristic ^{[ 2]} .[2] Zhang Huayan. the curve reestablishing [D] based on Delaunay triangulation and field represent: Zhengzhou University, 2011.

I geometric model that () utilizes Delaunay triangulation to set up is unique;

(ii) edge of what geometric model outermost layer border after Delaunay triangulation was formed an is convex polygon;

(iii) only can have an impact to closing on delta-shaped region when some summits being processed to the formation triangulation network, and other triangle is not affected;

(iv) the leg-of-mutton inside of Delaunay exists without any summit, and namely leg-of-mutton inside is hollow;

V leg-of-mutton minimum angle that () Delaunay triangulation obtains is maximum compared with other triangulation;

(vi) if mutually can exchange the diagonal line of any two adjacent convex quadrangles, therefore minimum interior angle angle can't become large;

(vii) Delaunay triangulation is formed triangle with nearest three.

Compared with the Delaunay triangulation of two-dimensional space, no matter three-dimensional Delaunay triangulation far exceedes two-dimensional space in difficulty or complicacy.Calling the CGAL algorithms library Delaunay subdivision carried out in three dimensions is a kind of very efficient method.Computational Geometry Algorithms Library, is called for short CGAL, is a C++ algorithm storehouse of being developed jointly by 8 research institutions of several countries and regions of Israel and European Union ^{[ 3]} .[3]Fabri A，Pion S.CGAL:the computational geometry algorithms library:proceedings of the Proceedings of the 17th ACM SIGSPATIAL international conference on advances in geographic information systems，2009[C].ACM.

CGAL provides the 3D surface grids subdivision bag of the geometric model reconstruct that can be directly used in 3-D view, and this subdivision bag is based on the Dlaunay triangulation of band restriction.

The key step that 3D surface grids subdivision bag subdivision algorithm carries out 3-D view geometry reconstruction is as follows:

1) select a part to do sample spot on the surface of object to calculate;

2) curved surface is extracted in the space body utilizing the 1st step to calculate;

3) insertion point is in the surface mesh extracted, and is optimized process to the new space body formed;

4) judging whether new space body network meets and build criterion, continuing repeated optimization if do not met;

5) as satisfied condition, then the process after continuing, once completes all grid reconstructions.

As can be seen here, building criterion is the key determining subdivision quality and result, and affects the final shape and size of grid.

3) under considering the condition of boundary layer influence, use in the three-dimensional description geometric model of Finite Element Method convection cell after above-mentioned subdivision and carry out numerical simulation, namely convection cell incompressible equation N-S equation of micro flow in duct solves, and obtains velocity field and the pressure field distribution of fluid flowing;

A. the basic thought of finite element method: by the continuous solving discrete region of research object be first a group limited and the unit combination body be connected with each other by certain way.Because unit can combine by different bind modes, and unit itself can have difformity again, and what therefore can be modeled to different geometries solves zonule; Then mechanical analysis is carried out to unit (zonule), finally holistic approach again.Thisly to break the whole up into parts, the basic ideas of collection zero to be whole method be exactly finite element ^[4].[4] Chen Xidong, Yang Jie. the current situation of finite element method and application [J]. Chinese manufacturing is information-based, and 2010,39 (11): 6-8

Finite element method is one of conventional computational fluid dynamics simulation method, there is strict mathematical theory, and do not need to simplify pore space, but the calculating of finite element carries out on the basis of geometric model, the digital cores by we being obtained someway is therefore needed to carry out geometry reconstruction.Simultaneously in order to reduce the calculated amount in three-dimensional geometry restructuring procedure, we also will take certain rock core hole treatment measures usually, and such as feature unit body is analyzed, the removal etc. of isolated pore.

Wherein said N-S equation:

Macroscopic view continuous model is based on continuous hypothesis, and fluid follows the mass conservation, momentum conservation and energy conservation, and its descriptive equation is Navier-Stokes equation (N-S equation):

$\frac{\∂\mathrm{\ρ}}{\∂t}+\▿\·\left(\mathrm{\ρu}\right)=0---\left(1.10\right)$

$\frac{\∂\left(\mathrm{\ρu}\right)}{\∂t}+\▿\·\left(\mathrm{\ρuu}\right)=\▿\·\mathrm{\σ}---\left(1.11\right)$

$\frac{\∂\left(\mathrm{\ρe}\right)}{\∂t}+\▿\·\left(\mathrm{\ρue}\right)=\mathrm{\σ}:\▿u-\▿\·q---\left(1.12\right)$

In formula, ρ is fluid macroscopic density; U is fluid velocity; T is temperature; E is interior energy; σ is stress tensor; Q is thermoflux; Above-mentioned system of equations is not closed, and requires that the parameter such as decompression force, speed must supplement the relation such as constitutive relation, state equation.Macroscopic view continuous model develops comparative maturity at present, the method be most widely used.

The concrete form of the N-S equation solved is as follows:

$\begin{array}{c}\mathrm{\ρ}\frac{\∂u}{\∂t}+\mathrm{\ρ}(u\·\▿)u=\▿\·[-\mathrm{\ρI}+\mathrm{\μ}(\▿u+{(\▿u)}^T)]\\ \mathrm{\ρ}\▿\·u=0\end{array}---\left(1.13\right)$

Wherein, ρ is the density of fluid; U is the speed of fluid; μ is the viscosity of fluid; Ι is unit matrix;

Boundary condition is:

$\begin{array}{c}p=p_0\\ [\mathrm{\μ}(\▿u)+{(\▿u)}^T]\·n=0\end{array}---\left(1.14\right)$

Wherein, p is pressure, p ₀for original pressure;

Described simulation steps comprise step 1.-5.:

1. the good 3-D geometric model of subdivision is imported;

2. univers parameter is arranged, and comprises fluid density, viscosity of crude and pressure;

3. select corresponding rock sample outlet, entrance and wall, calculate, obtain the speed of rock sample outlet and entrance;

4. integration is carried out to the speed of rock sample outlet and entrance, obtain the flow of rock sample outlet and entrance;

5. the result of simulation is drawn, obtain rock sample velocity field and pressure field distribution.

According to data obtained above, velocity field, pressure field distribution and streamline etc. can be drawn, the mobility status of fluid in digital cores can be found out by process.

By said method, what can contrast boundary layer in dissimilar rock affects situation.

The thickness in boundary layer and viscosity have a significant impact the seepage flow of fluid in microscopic void.The drag effect in boundary layer is more obvious in the rock that characteristic dimension is smaller, and main manifestations is that speed declines sooner, and pressure reduces faster; Boundary layer is more obvious on the impact of low permeability media, studies boundary layer even more important in least permeable medium.

Relative to prior art, the present invention has following beneficial effect:

1, the present invention utilizes Comsol to simulate, and can use manpower and material resources sparingly and the time.

2, the present invention utilizes digital cores technology, can the pore throat situation of accurate simulation true core and pressure, fluid distrbution, more tallies with the actual situation.

3, the present invention processes data result, can accurate computation bound layer thickness, closer to the Crude Oil in Boundary Layer thickness of true core.

Accompanying drawing explanation

Fig. 1 is that the present invention adopts maximum kind spacing method to split the algorithm flow chart of image.

Fig. 2 is the process flow diagram of simulation method step of the present invention.

Embodiment

Below in conjunction with embodiment, the present invention is described in detail, but is not limited thereto.

Embodiment,

A kind of low-permeability oil reservoirs boundary layer simulating method, comprises step as follows:

1) CT scan method structure digital cores comprises step (a)-(c):

A () obtains the multi-faceted CT projecting image data of rock sample by CT scan; Wherein the shape of rock sample does not affect CT projected image, and the CT that the present invention uses is MicroXCT-400 high-resolution 3DX ray microscope, can carry out 360 degree of imagings to rock sample, so the angle of shooting rock sample does not affect for rock sample projected image.In theory, the MicroXCT-400CT scanning device that the present invention uses supports that sample size reaches 300mm, and weight reaches 15Kg; Utilize CT to scan, individual CT projecting image data after scanning is automatically caught by Image Acquisition software and stores; By described rock sample after its axial-rotation fixed angle again to its CT scan, individual CT projecting image data after scanning is automatically caught by Image Acquisition software and stores, until described rock sample terminates to gather CT projecting image data after adding up rotating 360 degrees vertically;

The image-forming principle of CT: when X ray is through any material, its can with the atomic interaction of material and cause energy attenuation, and the different constituent of material has different absorption coefficients (i.e. attenuation coefficient) to X ray, say it on the contrary, the constituent of material can be judged by measurement of species to the absorption coefficient of X ray.When beam of x-rays is through object, it penetrate all substances on path to the summation of X-ray absorption coefficient all by the measurement result that is reflected in X-ray intensity, as shown in (1.1):

$I=I_o{e}^{-\underset{i}{\mathrm{\Σ}}{\mathrm{\μ}}_i{x}_i}---\left(1.1\right)$

In formula, I _ofor the initial strength of X ray, I is the intensity after X ray passes object, namely X ray be attenuated after intensity, i represents a certain constituent on the path passed at this X ray in material, μ _i, x _ibe respectively the i-th component to the attenuation coefficient of X ray and this component in the current length on path of X ray;

The image-forming principle of CT is set up on the basis of this just, by measuring the X ray through object section, afterwards, adopt certain method for reconstructing to calculate and object tomography locus absorption coefficient one to one, thus recover the structural information of object section;

During analog computation of the present invention, the concrete size of rock sample will be determined according to the attribute of rock sample, to reach best technique effect.For ensureing that rock core CT image can clear reflection rock core internal microstructure feature, experiment rock core size can not be excessive, reason mainly contains two aspects: on the one hand, because rock core density is larger, ray penetrate in rock core process be attenuated to relatively large extent cause its arrive detector system time signal very faint, if rock core is oversize, the complete conductively-closed of ray may be caused; On the other hand, the maximum lattice point quantity adopted when solving due to software systems is certain, and rock core size increases the physical size increase meaning that each lattice point is corresponding, and namely the resolution of image reduces.For existing industrial CT equipment, obtain CT image more clearly, rock core cylinder largest diameter is centimetre-sized;

B () utilizes image rebuilding method that described CT projecting image data is converted to the gray level image of described rock sample;

The core of CT imaging is how by reconstructs projection data image.The essence of image reconstruction is the absorption coefficient being solved each pixel in image array by the CT data for projection after gathering, and then re-constructs the process of image; The method for solving of image reconstruction problem can be divided into two large classes according to its feature: the first kind is converter technique (being also analytical method), is characterized in that last discretize utilizes computer calculate first in continuous domain dissection process; Equations of The Second Kind is process of iteration (being also algebraic reconstruction technique, Series Expansion Method, optimisation technique etc.), is characterized in that directly carrying out discretize analysis obtains numerical solution.Comparatively speaking, converter technique advantages: realize simple, speed is fast, can obtain good reconstruction quality to enough accurate data for projection.Therefore, extensively converter technique is adopted in current practical CT system.

Described image rebuilding method is converter technique, comprise not containing the Fourier method for reconstructing FRM of back projection's step with based on the convolution of back projection or filtering method: the former with direct Fourier transform method DFM for representative, the latter with convolution or filtered back projection C/FBP method for representative;

The basic thought of DFM is according to the mechanism of action between radiation source and imageable target, search out objective function f (x, y) relation between two-dimensional fourier transform (FT) and the one dimension FT of projection function, utilize this relation, although objective function f is (x, y) be unknown, but its two-dimentional FT obtains by the one dimension FT of projection function in different azimuth, once obtain f (x, y) all projections, also two-dimentional FT, the f (x, y) that just obtain it then can be brought by the contravariant of two-dimentional FT to be obtained;

Come from Radon based on the convolution of back projection or filtering method C/FBP method to negate the direct realization of formula; C/FBP formula utilizes FT and projection theorem to derive; If the most high spatial frequency of p (x', θ) is B;

The mathematical formulae of C/FBP is write as

$f(x,y)={\∫}_0^{\mathrm{\π}}q(x^{\′},\mathrm{\θ})\mathrm{d\θ}---\left(1.2\right)$

Wherein,

$q(x^{\′},\mathrm{\θ})={\∫}_{-\∞}^{\∞}\left|R\right|P(R,\mathrm{\θ})\exp \left(j2\mathrm{\π}{x}^{\′}R\right)\mathrm{dR}---\left(1.3\right)$

Or

q(x',θ)＝c(x')*p(x',θ) (1.4)

In formula (1.2), (1.3), (1.4), the objective function that f (x, y) is image reconstruction, namely rock sample is to the distribution function of gamma ray absorption coefficient, and x, y are respectively the rectangular coordinate of certain point on rock sample scanning section; P (R, θ), p (x', θ) be f (x, y) along the projection function on θ direction, i.e. the transmitted intensity signal that detects of CT detector, x' is certain point (x on rock sample scanning section, y) horizontal ordinate in the rotated coordinate system, R for this point is to rotation center, i.e. true origin, distance; C (x') is convolution function, and it is by wave filter | the Fourier inverse transformation of R| (| R|≤B) is tried to achieve;

Q (x', θ) in formula (1.3), (1.4) is called projection P (R, θ), the filter projection of p (x', θ) and convolution projection; The image rebuilding method realized by formula (1.2), (1.3) is called filtered back-projection method FBP, and the image rebuilding method realized by formula (1.2), (1.4) is called the convolution back projection method CBP; Therefore, as obtained the projection of image, then oppositely solved by convolution or filtered back-projection method and obtain the distribution function f (x, y) of rock sample section to gamma ray absorption coefficient, the scan image of rock sample section is just obtained after adopting the mode of gray level image to express f (x, y);

C () adopts image binary segmentation method to be separated space, gray level image mesoporosity and the rock sample skeleton of described rock sample, set up digital cores;

In gray level image, pore space and rock skeleton are divided it by image partition method: Iamge Segmentation is choosing of segmentation threshold, if threshold value is chosen unreasonable, be easy to destination object to understand become background image, or background is interpreted as destination object, so just cannot effective evaluating objects feature when carrying out signature analysis to target image;

Threshold definitions is as follows: suppose that F (i, j) represents image pixel gray level value in gray level image, and the pixel that gray-scale value is greater than threshold value T after being separated represents destination object, i.e. hole, represents with 1; Other pixel represents image background, i.e. skeleton, represents with 0, then can obtain destination object by following formula:

$g(i,j)=\left\{\begin{array}{cc}1;& F(i,j)\≥T\\ 0;& \mathrm{else}\end{array}\right.---\left(1.5\right)$

In formula, T is threshold value, also known as threshold value;

Based on the image binary segmentation method of core porosity

The fundamental physical quantity that core porosity must measure when being and carrying out reservoir physics experiment; Under existing experiment condition, measurement accurately can be facilitated to obtain by gas expansion method, hold-up method etc.; Because factor of porosity gives the share shared by space, rock core mesoporosity, therefore, in known cases, factor of porosity is attached to the segmentation quality that can improve image in the cutting procedure of rock core CT image undoubtedly as constraint condition;

If the factor of porosity of rock sample is maximum, the minimum gradation value of image are respectively I _max, I _mIN, gray-scale value is the pixel count of i is p (i), then meet the gray-scale value k of formula (1.6) ^*i.e. required segmentation threshold:

Image binary segmentation method based on rock sample factor of porosity is actually selects suitable segmentation threshold, and the ratio making the core image mesoporosity pixel count obtained by this Threshold segmentation account for total pixel number equals core porosity.Because grey scale pixel value is by the integer representation in 0 to 255 intervals, therefore, splitting the image factor of porosity obtained may exist nuance with core experiment factor of porosity, but this does not affect Iamge Segmentation quality.

Find by carrying out segmentation test to experimental image, image binary segmentation method based on core porosity can obtain high-quality segmentation result, therefore, when core experiment factor of porosity is known, binary segmentation can be carried out by the method to rock core CT image.But, and the rock core that not all is used for CT experimental study has experimental port porosity, when not having factor of porosity data, need split rock core CT image by other image binary segmentation method.In image binary segmentation method, maximum kind spacing method is because principle is simple, program design facilitates and can obtain good segmentation effect and be widely used.

When not having factor of porosity, maximum kind spacing method is utilized to split gray level image; What formed by the gray-scale value of pixels all in gray level image to be split is totally considered as a set, assuming that the set threshold value of gray-scale value is divided into two groups, then when the variance (between-group variance) of mean value of two groups and the ratio of the variance (interclass variance) of each group reach maximum, the segmentation threshold corresponding to it is reasonable value; The pore space of described rock sample inside respectively becomes one group with skeleton, and the difference of each group inside is very little and nature difference between two groups is very large; After two components cut open by the rational threshold value of employing, nature difference between two groups will reach maximum value with the ratio of difference in group, visible, feature and the segmentation thinking of maximum kind spacing method to image of rock sample CT scan projected image itself are very identical, and this is also that maximum kind spacing method segmentation rock core CT image can obtain good segmentation result and the major reason that is widely adopted;

The principle of maximum kind spacing method selected threshold is as follows.If maximum, the minimum gradation value of image are respectively I _max, I _mIN, gray-scale value is the pixel count of i is p (i); Gray threshold is k, and pixel is divided into two groups by it: group 1 (being assumed to be pore space) and group 2 (being assumed to be rock skeleton).Add up respectively and calculate the pixel count N of two groups ₁(k), N ₂(k), average gray value μ ₁(k), μ ₂(k) and variance and the average gray value μ of whole pixel _t, then interclass variance and between-group variance be respectively:

${\mathrm{\σ}}_W^2=N_1{\mathrm{\σ}}_1^2+N_2{\mathrm{\σ}}_2^2---\left(1.7\right)$

${\mathrm{\σ}}_B^2=N_1{({\mathrm{\μ}}_1-{\mathrm{\μ}}_T)}^2+N_2{({\mathrm{\μ}}_2-{\mathrm{\μ}}_T)}^2---\left(1.8\right)$

In formula:

${\mathrm{\μ}}_1\left(k\right)=\frac{\underset{i=I_{\mathrm{MIN}}}{\overset{k}{\mathrm{\Σ}}}i\·p\left(i\right)}{N_1\left(k\right)}$ ${\mathrm{\μ}}_2\left(k\right)=\frac{\underset{i=k+1}{\overset{I_{\mathrm{MAX}}}{\mathrm{\Σ}}}i\·p\left(i\right)}{N_2\left(k\right)}$

${\mathrm{\μ}}_T\left(k\right)=\frac{\underset{i=I_{\mathrm{MIN}}}{\overset{I_{\mathrm{MAX}}}{\mathrm{\Σ}}}i\·p\left(i\right)}{N_1\left(k\right)+N_2\left(k\right)}$

${\mathrm{\σ}}_1^2\left(k\right)=\frac{\underset{i=I_{\mathrm{MIN}}}{\overset{k}{\mathrm{\Σ}}}{[i-{\mathrm{\μ}}_1\left(k\right)]}^2\·p\left(i\right)}{N_1\left(k\right)}$ ${\mathrm{\σ}}_2^2\left(k\right)=\frac{\underset{i=k+1}{\overset{I_{\mathrm{MAX}}}{\mathrm{\Σ}}}{[i-{\mathrm{\μ}}_2\left(k\right)]}^2\·p\left(i\right)}{N_2\left(k\right)}$

Select optimum gradation value k ^*, make it meet (1.9) formula, then k ^*i.e. required threshold value.

$f\left(k^{*}\right)=\max \{f\left(k\right)={\mathrm{\σ}}_B^2\left(k\right)/{\mathrm{\σ}}_W^2\left(k\right)\}---\left(1.9\right);$

Adopt the algorithm flow of maximum kind spacing method segmentation image, as shown in Figure 1;

2) use Delaunay triangulation to carry out subdivision to the above-mentioned digital cores built, obtain 3-D geometric model;

Wherein, described Delaunay triangulation is also a kind of iso-surface patch method, its advantage be in the process building geometric model by controlling the parameter of geometric units, make the geometric model formed be more conducive to after analysis and calculation.Due to the plurality of advantages of Delaunay triangulation, this method has been widely used in multiple field such as analysis, medical science of finite element ^[1].[1] Ding Yongxiang, the summer huge Zan. the Delaunay triangulation [J] of arbitrary polygon. Chinese journal of computers, 1994,17 (4): 270-275.

Delaunay triangulation is the triangulation optimized, and has a lot of good characteristic ^{[ 2]} .[2] Zhang Huayan. the curve reestablishing [D] based on Delaunay triangulation and field represent: Zhengzhou University, 2011.

I geometric model that () utilizes Delaunay triangulation to set up is unique;

(ii) edge of what geometric model outermost layer border after Delaunay triangulation was formed an is convex polygon;

(iii) only can have an impact to closing on delta-shaped region when some summits being processed to the formation triangulation network, and other triangle is not affected;

(iv) the leg-of-mutton inside of Delaunay exists without any summit, and namely leg-of-mutton inside is hollow;

V leg-of-mutton minimum angle that () Delaunay triangulation obtains is maximum compared with other triangulation;

(vi) if mutually can exchange the diagonal line of any two adjacent convex quadrangles, therefore minimum interior angle angle can't become large;

(vii) Delaunay triangulation is formed triangle with nearest three.

Compared with the Delaunay triangulation of two-dimensional space, no matter three-dimensional Delaunay triangulation far exceedes two-dimensional space in difficulty or complicacy.Calling the CGAL algorithms library Delaunay subdivision carried out in three dimensions is a kind of very efficient method.Computational Geometry Algorithms Library, is called for short CGAL, is a C++ algorithm storehouse of being developed jointly by 8 research institutions of several countries and regions of Israel and European Union ^{[ 3]} .[3]Fabri A，Pion S.CGAL:the computational geometry algorithms library:proceedings of the Proceedings of the 17th ACM SIGSPATIAL international conference on advances in geographic information systems，2009[C].ACM.

CGAL provides the 3D surface grids subdivision bag of the geometric model reconstruct that can be directly used in 3-D view, and this subdivision bag is based on the Dlaunay triangulation of band restriction.

The key step that 3D surface grids subdivision bag subdivision algorithm carries out 3-D view geometry reconstruction is as follows:

1) select a part to do sample spot on the surface of object to calculate;

2) curved surface is extracted in the space body utilizing the 1st step to calculate;

3) insertion point is in the surface mesh extracted, and is optimized process to the new space body formed;

4) judging whether new space body network meets and build criterion, continuing repeated optimization if do not met;

5) as satisfied condition, then the process after continuing, once completes all grid reconstructions.

As can be seen here, building criterion is the key determining subdivision quality and result, and affects the final shape and size of grid.

3) under considering the condition of boundary layer influence, use in the three-dimensional description geometric model of Finite Element Method convection cell after above-mentioned subdivision and carry out numerical simulation, namely convection cell incompressible equation N-S equation of micro flow in duct solves, and obtains velocity field and the pressure field distribution of fluid flowing;

A. the basic thought of finite element method: by the continuous solving discrete region of research object be first a group limited and the unit combination body be connected with each other by certain way.Because unit can combine by different bind modes, and unit itself can have difformity again, and what therefore can be modeled to different geometries solves zonule; Then mechanical analysis is carried out to unit (zonule), finally holistic approach again.Thisly to break the whole up into parts, the basic ideas of collection zero to be whole method be exactly finite element ^[4].[4] Chen Xidong, Yang Jie. the current situation of finite element method and application [J]. Chinese manufacturing is information-based, and 2010,39 (11): 6-8

Finite element method is one of conventional computational fluid dynamics simulation method, there is strict mathematical theory, and do not need to simplify pore space, but the calculating of finite element carries out on the basis of geometric model, the digital cores by we being obtained someway is therefore needed to carry out geometry reconstruction.Simultaneously in order to reduce the calculated amount in three-dimensional geometry restructuring procedure, we also will take certain rock core hole treatment measures usually, and such as feature unit body is analyzed, the removal etc. of isolated pore.

Wherein said N-S equation:

Macroscopic view continuous model is based on continuous hypothesis, and fluid follows the mass conservation, momentum conservation and energy conservation, and its descriptive equation is Navier-Stokes equation (N-S equation):

$\frac{\∂\mathrm{\ρ}}{\∂t}+\▿\·\left(\mathrm{\ρu}\right)=0---\left(1.10\right)$

$\frac{\∂\left(\mathrm{\ρu}\right)}{\∂t}+\▿\·\left(\mathrm{\ρuu}\right)=\▿\·\mathrm{\σ}---\left(1.11\right)$

$\frac{\∂\left(\mathrm{\ρe}\right)}{\∂t}+\▿\·\left(\mathrm{\ρue}\right)=\mathrm{\σ}:\▿u-\▿\·q---\left(1.12\right)$

In formula, ρ is fluid macroscopic density; U is fluid velocity; T is temperature; E is interior energy; σ is stress tensor; Q is thermoflux; Above-mentioned system of equations is not closed, and requires that the parameter such as decompression force, speed must supplement the relation such as constitutive relation, state equation.Macroscopic view continuous model develops comparative maturity at present, the method be most widely used.

The concrete form of the N-S equation solved is as follows:

$\begin{array}{c}\mathrm{\ρ}\frac{\∂u}{\∂t}+\mathrm{\ρ}(u\·\▿)u=\▿\·[-\mathrm{\ρI}+\mathrm{\μ}(\▿u+{(\▿u)}^T)]\\ \mathrm{\ρ}\▿\·u=0\end{array}---\left(1.13\right)$

Wherein, ρ is the density of fluid; U is the speed of fluid; μ is the viscosity of fluid; Ι is unit matrix;

Boundary condition is:

$\begin{array}{c}p=p_0\\ [\mathrm{\μ}(\▿u)+{(\▿u)}^T]\·n=0\end{array}---\left(1.14\right)$

Wherein, p is pressure, p ₀for original pressure;

Described simulation steps comprise step 1.-5.:

1. the good 3-D geometric model of subdivision is imported;

2. univers parameter is arranged, and comprises fluid density, viscosity of crude and pressure;

3. select corresponding rock sample outlet, entrance and wall, calculate, obtain the speed of rock sample outlet and entrance;

4. integration is carried out to the speed of rock sample outlet and entrance, obtain the flow of rock sample outlet and entrance;

5. the result of simulation is drawn, obtain rock sample velocity field and pressure field distribution.

According to data obtained above, velocity field, pressure field distribution and streamline etc. can be drawn, the mobility status of fluid in digital cores can be found out by process.

By said method, what can contrast boundary layer in dissimilar rock affects situation.

The thickness in boundary layer and viscosity have a significant impact the seepage flow of fluid in microscopic void.The drag effect in boundary layer is more obvious in the rock that characteristic dimension is smaller, and main manifestations is that speed declines sooner, and pressure reduces faster; Boundary layer is more obvious on the impact of low permeability media, studies boundary layer even more important in least permeable medium.

Claims (1)

1. a low-permeability oil reservoirs boundary layer simulating method, is characterized in that, it is as follows that described method comprises step:

1) CT scan method structure digital cores comprises step (a)-(c):

A () obtains the multi-faceted CT projecting image data of rock sample by CT scan; Utilize CT to scan, individual CT projecting image data after scanning is automatically caught by Image Acquisition software and stores; By described rock sample after its axial-rotation fixed angle again to its CT scan, individual CT projecting image data after scanning is automatically caught by Image Acquisition software and stores, until described rock sample terminates to gather CT projecting image data after adding up rotating 360 degrees vertically;

B () utilizes image rebuilding method that described CT projecting image data is converted to the gray level image of described rock sample;

Described image rebuilding method is converter technique, comprise not containing the Fourier method for reconstructing FRM of back projection's step with based on the convolution of back projection or filtering method: the former with direct Fourier transform method DFM for representative, the latter with convolution or filtered back projection C/FBP method for representative;

The basic thought of DFM is according to the mechanism of action between radiation source and imageable target, search out objective function f (x, y) relation between two-dimensional fourier transform (FT) and the one dimension FT of projection function, utilize this relation, although objective function f is (x, y) be unknown, but its two-dimentional FT obtains by the one dimension FT of projection function in different azimuth, once obtain f (x, y) all projections, also two-dimentional FT, the f (x, y) that just obtain it then can be brought by the contravariant of two-dimentional FT to be obtained;

Come from Radon based on the convolution of back projection or filtering method C/FBP method to negate the direct realization of formula; C/FBP formula utilizes FT and projection theorem to derive; If the most high spatial frequency of p (x', θ) is B;

The mathematical formulae of C/FBP is write as

$f(x,y)={\∫}_0^{\mathrm{\π}}q(x^{\′},\mathrm{\θ})\mathrm{d\θ}---\left(1.2\right)$

Wherein,

$q(x^{\′},\mathrm{\θ})={\∫}_{-\∞}^{\∞}\left|R\right|P(R,\mathrm{\θ})\exp \left(j2\mathrm{\π}{x}^{\′}R\right)\mathrm{dR}---(1.3)$

Or

q(x',θ)＝c(x')*p(x',θ) (1.4)

In formula (1.2), (1.3), (1.4), the objective function that f (x, y) is image reconstruction, namely rock sample is to the distribution function of gamma ray absorption coefficient, and x, y are respectively the rectangular coordinate of certain point on rock sample scanning section; P (R, θ), p (x', θ) be f (x, y) along the projection function on θ direction, i.e. the transmitted intensity signal that detects of CT detector, x' is certain point (x on rock sample scanning section, y) horizontal ordinate in the rotated coordinate system, R for this point is to rotation center, i.e. true origin, distance; C (x') is convolution function, and it is by wave filter | the Fourier inverse transformation of R| (| R|≤B) is tried to achieve;

Q (x', θ) in formula (1.3), (1.4) is called projection P (R, θ), the filter projection of p (x', θ) and convolution projection; The image rebuilding method realized by formula (1.2), (1.3) is called filtered back-projection method FBP, and the image rebuilding method realized by formula (1.2), (1.4) is called the convolution back projection method CBP; Therefore, as obtained the projection of image, then oppositely solved by convolution or filtered back-projection method and obtain the distribution function f (x, y) of rock sample section to gamma ray absorption coefficient, the scan image of rock sample section is just obtained after adopting the mode of gray level image to express f (x, y);

C () adopts image binary segmentation method to be separated space, gray level image mesoporosity and the rock sample skeleton of described rock sample, set up digital cores;

In gray level image, pore space and rock skeleton are divided it by image partition method: Iamge Segmentation is choosing of segmentation threshold;

Threshold definitions is as follows: suppose that F (i, j) represents image pixel gray level value in gray level image, and the pixel that gray-scale value is greater than threshold value T after being separated represents destination object, i.e. hole, represents with 1; Other pixel represents image background, i.e. skeleton, represents with 0, then can obtain destination object by following formula:

$g(i,j)=\left\{\begin{array}{cc}1;& F(i,j)\≥T\\ 0;& \mathrm{else}\end{array}\right.---\left(1.5\right)$

In formula, T is threshold value, also known as threshold value;

Based on the image binary segmentation method of core porosity

If the factor of porosity of rock sample is maximum, the minimum gradation value of image are respectively I _max, I _mIN, gray-scale value is the pixel count of i is p (i), then meet the gray-scale value k of formula (1.6) ^*i.e. required segmentation threshold:

When not having factor of porosity, maximum kind spacing method is utilized to split gray level image;

2) use Delaunay triangulation to carry out subdivision to the above-mentioned digital cores built, obtain 3-D geometric model;

3) under considering the condition of boundary layer influence, use in the three-dimensional description geometric model of Finite Element Method convection cell after above-mentioned subdivision and carry out numerical simulation, namely convection cell incompressible equation N-S equation of micro flow in duct solves, and obtains velocity field and the pressure field distribution of fluid flowing;

Wherein said N-S equation:

Macroscopic view continuous model is based on continuous hypothesis, and fluid follows the mass conservation, momentum conservation and energy conservation, and its descriptive equation is Navier-Stokes equation (N-S equation):

$\frac{\∂\mathrm{\ρ}}{\∂t}+\▿\·\left(\mathrm{\ρu}\right)=0---\left(1.10\right)$

$\frac{\∂\left(\mathrm{\ρu}\right)}{\∂t}+\▿\·\left(\mathrm{\ρuu}\right)=\▿\·\mathrm{\σ}---\left(1.11\right)$

$\frac{\∂\left(\mathrm{\ρe}\right)}{\∂t}+\▿\·\left(\mathrm{\ρue}\right)=\mathrm{\σ}:\▿u-\▿\·q---\left(1.12\right)$

In formula, ρ is fluid macroscopic density; U is fluid velocity; T is temperature; E is interior energy; σ is stress tensor; Q is thermoflux;

The concrete form of the N-S equation solved is as follows:

$\begin{array}{c}\mathrm{\ρ}\frac{\∂u}{\∂t}+\mathrm{\ρ}(u\·\▿)u=\▿\·[-\mathrm{\ρI}+\mathrm{\μ}(\▿u+{(\▿u)}^T)]\\ \mathrm{\ρ}\▿\·u=0\end{array}---\left(1.13\right)$

Wherein, ρ is the density of fluid; U is the speed of fluid; μ is the viscosity of fluid; Ι is unit matrix;

Boundary condition is:

$\begin{array}{c}p=p_0\\ [\mathrm{\μ}(\▿u)+{(\▿u)}^T]\·n=0\end{array}---\left(1.14\right)$

Wherein, p is pressure, p ₀for original pressure;

Described simulation steps comprise step 1.-5.:

1. the good 3-D geometric model of subdivision is imported;

2. univers parameter is arranged, and comprises fluid density, viscosity of crude and pressure;

3. select corresponding rock sample outlet, entrance and wall, calculate, obtain the speed of rock sample outlet and entrance;

4. integration is carried out to the speed of rock sample outlet and entrance, obtain the flow of rock sample outlet and entrance;

5. the result of simulation is drawn, obtain rock sample velocity field and pressure field distribution.