CN113919126A - 油气层损害类型和程度时空演化4d定量与智能诊断方法及其*** - Google Patents
油气层损害类型和程度时空演化4d定量与智能诊断方法及其*** Download PDFInfo
- Publication number
- CN113919126A CN113919126A CN202110989424.5A CN202110989424A CN113919126A CN 113919126 A CN113919126 A CN 113919126A CN 202110989424 A CN202110989424 A CN 202110989424A CN 113919126 A CN113919126 A CN 113919126A
- Authority
- CN
- China
- Prior art keywords
- reservoir
- determining
- damage
- equation
- fluid
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 319
- 238000003745 diagnosis Methods 0.000 title abstract description 14
- 239000004215 Carbon black (E152) Substances 0.000 title description 8
- 229930195733 hydrocarbon Natural products 0.000 title description 8
- 150000002430 hydrocarbons Chemical class 0.000 title description 8
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims abstract description 275
- 238000004088 simulation Methods 0.000 claims abstract description 247
- 238000002347 injection Methods 0.000 claims abstract description 71
- 239000007924 injection Substances 0.000 claims abstract description 71
- 239000012530 fluid Substances 0.000 claims description 423
- 239000011148 porous material Substances 0.000 claims description 313
- 230000035699 permeability Effects 0.000 claims description 311
- 239000002245 particle Substances 0.000 claims description 288
- 238000009826 distribution Methods 0.000 claims description 211
- 239000003921 oil Substances 0.000 claims description 181
- 229920000642 polymer Polymers 0.000 claims description 179
- 239000012071 phase Substances 0.000 claims description 166
- 230000008569 process Effects 0.000 claims description 122
- 239000004576 sand Substances 0.000 claims description 117
- 238000009792 diffusion process Methods 0.000 claims description 116
- 241000894006 Bacteria Species 0.000 claims description 114
- 239000011435 rock Substances 0.000 claims description 113
- 238000001179 sorption measurement Methods 0.000 claims description 108
- 238000004519 manufacturing process Methods 0.000 claims description 96
- 150000002500 ions Chemical class 0.000 claims description 92
- 239000010779 crude oil Substances 0.000 claims description 74
- 230000000694 effects Effects 0.000 claims description 70
- 238000013508 migration Methods 0.000 claims description 67
- 230000005012 migration Effects 0.000 claims description 65
- 230000006870 function Effects 0.000 claims description 62
- 230000008859 change Effects 0.000 claims description 52
- 239000004927 clay Substances 0.000 claims description 48
- 239000007790 solid phase Substances 0.000 claims description 38
- 238000001556 precipitation Methods 0.000 claims description 37
- 230000008961 swelling Effects 0.000 claims description 36
- 238000009736 wetting Methods 0.000 claims description 34
- 230000006735 deficit Effects 0.000 claims description 32
- 238000004945 emulsification Methods 0.000 claims description 31
- 235000015097 nutrients Nutrition 0.000 claims description 29
- 230000033001 locomotion Effects 0.000 claims description 27
- 238000006243 chemical reaction Methods 0.000 claims description 26
- 239000007788 liquid Substances 0.000 claims description 26
- 238000005315 distribution function Methods 0.000 claims description 25
- 230000035945 sensitivity Effects 0.000 claims description 25
- 239000008346 aqueous phase Substances 0.000 claims description 23
- 239000013049 sediment Substances 0.000 claims description 20
- 239000007787 solid Substances 0.000 claims description 19
- 238000003795 desorption Methods 0.000 claims description 13
- 239000000839 emulsion Substances 0.000 claims description 12
- 206010013496 Disturbance in attention Diseases 0.000 claims description 11
- 238000003860 storage Methods 0.000 claims description 10
- 239000000758 substrate Substances 0.000 claims description 6
- 238000005553 drilling Methods 0.000 claims description 5
- 230000001010 compromised effect Effects 0.000 claims 1
- 238000011161 development Methods 0.000 abstract description 50
- 210000003491 skin Anatomy 0.000 description 70
- 239000002244 precipitate Substances 0.000 description 29
- 238000010586 diagram Methods 0.000 description 27
- 239000003129 oil well Substances 0.000 description 27
- 238000000151 deposition Methods 0.000 description 26
- 230000014509 gene expression Effects 0.000 description 25
- 210000002615 epidermis Anatomy 0.000 description 18
- 230000000875 corresponding effect Effects 0.000 description 17
- 230000008021 deposition Effects 0.000 description 17
- 239000000126 substance Substances 0.000 description 15
- 230000000903 blocking effect Effects 0.000 description 13
- 238000004364 calculation method Methods 0.000 description 13
- 239000011159 matrix material Substances 0.000 description 13
- 235000020681 well water Nutrition 0.000 description 12
- 239000002349 well water Substances 0.000 description 12
- 238000013461 design Methods 0.000 description 11
- 238000005516 engineering process Methods 0.000 description 11
- 230000001580 bacterial effect Effects 0.000 description 10
- 230000002829 reductive effect Effects 0.000 description 10
- 239000000243 solution Substances 0.000 description 10
- 238000007906 compression Methods 0.000 description 8
- 238000006073 displacement reaction Methods 0.000 description 8
- 230000006835 compression Effects 0.000 description 7
- 230000001186 cumulative effect Effects 0.000 description 7
- 230000009545 invasion Effects 0.000 description 7
- VTYYLEPIZMXCLO-UHFFFAOYSA-L Calcium carbonate Chemical compound [Ca+2].[O-]C([O-])=O VTYYLEPIZMXCLO-UHFFFAOYSA-L 0.000 description 6
- TZCXTZWJZNENPQ-UHFFFAOYSA-L barium sulfate Chemical compound [Ba+2].[O-]S([O-])(=O)=O TZCXTZWJZNENPQ-UHFFFAOYSA-L 0.000 description 6
- 230000015572 biosynthetic process Effects 0.000 description 6
- OSGAYBCDTDRGGQ-UHFFFAOYSA-L calcium sulfate Chemical compound [Ca+2].[O-]S([O-])(=O)=O OSGAYBCDTDRGGQ-UHFFFAOYSA-L 0.000 description 6
- 238000004090 dissolution Methods 0.000 description 6
- 230000005484 gravity Effects 0.000 description 6
- 238000010521 absorption reaction Methods 0.000 description 5
- XDFCIPNJCBUZJN-UHFFFAOYSA-N barium(2+) Chemical compound [Ba+2] XDFCIPNJCBUZJN-UHFFFAOYSA-N 0.000 description 5
- 238000004134 energy conservation Methods 0.000 description 5
- 230000003902 lesion Effects 0.000 description 5
- 239000000203 mixture Substances 0.000 description 5
- 230000009467 reduction Effects 0.000 description 5
- 238000011160 research Methods 0.000 description 5
- PWYYWQHXAPXYMF-UHFFFAOYSA-N strontium(2+) Chemical compound [Sr+2] PWYYWQHXAPXYMF-UHFFFAOYSA-N 0.000 description 5
- 238000009825 accumulation Methods 0.000 description 4
- 230000008901 benefit Effects 0.000 description 4
- 238000002474 experimental method Methods 0.000 description 4
- 210000000540 fraction c Anatomy 0.000 description 4
- 230000006872 improvement Effects 0.000 description 4
- 230000007246 mechanism Effects 0.000 description 4
- 230000000704 physical effect Effects 0.000 description 4
- 238000011084 recovery Methods 0.000 description 4
- UBXAKNTVXQMEAG-UHFFFAOYSA-L strontium sulfate Chemical compound [Sr+2].[O-]S([O-])(=O)=O UBXAKNTVXQMEAG-UHFFFAOYSA-L 0.000 description 4
- 102000004190 Enzymes Human genes 0.000 description 3
- 108090000790 Enzymes Proteins 0.000 description 3
- 238000005411 Van der Waals force Methods 0.000 description 3
- 230000009471 action Effects 0.000 description 3
- 229910052925 anhydrite Inorganic materials 0.000 description 3
- 239000011575 calcium Substances 0.000 description 3
- 229910000019 calcium carbonate Inorganic materials 0.000 description 3
- 238000012512 characterization method Methods 0.000 description 3
- 238000011156 evaluation Methods 0.000 description 3
- 230000004907 flux Effects 0.000 description 3
- 230000000977 initiatory effect Effects 0.000 description 3
- 238000005259 measurement Methods 0.000 description 3
- 238000005457 optimization Methods 0.000 description 3
- 230000036961 partial effect Effects 0.000 description 3
- 230000003068 static effect Effects 0.000 description 3
- 238000012360 testing method Methods 0.000 description 3
- OKTJSMMVPCPJKN-UHFFFAOYSA-N Carbon Chemical compound [C] OKTJSMMVPCPJKN-UHFFFAOYSA-N 0.000 description 2
- 230000001133 acceleration Effects 0.000 description 2
- 150000001450 anions Chemical class 0.000 description 2
- 238000013459 approach Methods 0.000 description 2
- 230000006399 behavior Effects 0.000 description 2
- 238000009835 boiling Methods 0.000 description 2
- 229910052799 carbon Inorganic materials 0.000 description 2
- 150000001768 cations Chemical class 0.000 description 2
- 210000004027 cell Anatomy 0.000 description 2
- 239000000084 colloidal system Substances 0.000 description 2
- 238000000354 decomposition reaction Methods 0.000 description 2
- 238000005421 electrostatic potential Methods 0.000 description 2
- 210000003495 flagella Anatomy 0.000 description 2
- 230000010354 integration Effects 0.000 description 2
- 230000003993 interaction Effects 0.000 description 2
- 230000001788 irregular Effects 0.000 description 2
- 230000000670 limiting effect Effects 0.000 description 2
- 239000007791 liquid phase Substances 0.000 description 2
- 239000012528 membrane Substances 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 230000037361 pathway Effects 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 229920006395 saturated elastomer Polymers 0.000 description 2
- 238000009991 scouring Methods 0.000 description 2
- 238000004062 sedimentation Methods 0.000 description 2
- 239000008247 solid mixture Substances 0.000 description 2
- 230000007480 spreading Effects 0.000 description 2
- 238000003892 spreading Methods 0.000 description 2
- 230000009466 transformation Effects 0.000 description 2
- 230000005514 two-phase flow Effects 0.000 description 2
- 238000012795 verification Methods 0.000 description 2
- BVKZGUZCCUSVTD-UHFFFAOYSA-M Bicarbonate Chemical compound OC([O-])=O BVKZGUZCCUSVTD-UHFFFAOYSA-M 0.000 description 1
- 230000005653 Brownian motion process Effects 0.000 description 1
- OYPRJOBELJOOCE-UHFFFAOYSA-N Calcium Chemical compound [Ca] OYPRJOBELJOOCE-UHFFFAOYSA-N 0.000 description 1
- 241000084490 Esenbeckia delta Species 0.000 description 1
- 241000224466 Giardia Species 0.000 description 1
- 238000004510 Lennard-Jones potential Methods 0.000 description 1
- QAOWNCQODCNURD-UHFFFAOYSA-L Sulfate Chemical compound [O-]S([O-])(=O)=O QAOWNCQODCNURD-UHFFFAOYSA-L 0.000 description 1
- 238000005076 Van der Waals potential Methods 0.000 description 1
- 230000004913 activation Effects 0.000 description 1
- 239000008186 active pharmaceutical agent Substances 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000037358 bacterial metabolism Effects 0.000 description 1
- 229910052788 barium Inorganic materials 0.000 description 1
- DSAJWYNOEDNPEQ-UHFFFAOYSA-N barium atom Chemical compound [Ba] DSAJWYNOEDNPEQ-UHFFFAOYSA-N 0.000 description 1
- 238000005452 bending Methods 0.000 description 1
- 238000005537 brownian motion Methods 0.000 description 1
- 229910052791 calcium Inorganic materials 0.000 description 1
- 230000003197 catalytic effect Effects 0.000 description 1
- 229910052923 celestite Inorganic materials 0.000 description 1
- 230000035605 chemotaxis Effects 0.000 description 1
- 230000002301 combined effect Effects 0.000 description 1
- 239000002131 composite material Substances 0.000 description 1
- 230000001276 controlling effect Effects 0.000 description 1
- 230000002596 correlated effect Effects 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 235000011850 desserts Nutrition 0.000 description 1
- 230000006866 deterioration Effects 0.000 description 1
- 238000002405 diagnostic procedure Methods 0.000 description 1
- 230000004069 differentiation Effects 0.000 description 1
- 239000006185 dispersion Substances 0.000 description 1
- 229920001971 elastomer Polymers 0.000 description 1
- 239000000806 elastomer Substances 0.000 description 1
- 238000000556 factor analysis Methods 0.000 description 1
- 239000010419 fine particle Substances 0.000 description 1
- 229910052602 gypsum Inorganic materials 0.000 description 1
- XLYOFNOQVPJJNP-ZSJDYOACSA-N heavy water Substances [2H]O[2H] XLYOFNOQVPJJNP-ZSJDYOACSA-N 0.000 description 1
- 238000011065 in-situ storage Methods 0.000 description 1
- 230000002427 irreversible effect Effects 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 238000013178 mathematical model Methods 0.000 description 1
- 230000002503 metabolic effect Effects 0.000 description 1
- 230000004060 metabolic process Effects 0.000 description 1
- 229910052751 metal Inorganic materials 0.000 description 1
- 239000002184 metal Substances 0.000 description 1
- 239000011859 microparticle Substances 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000005325 percolation Methods 0.000 description 1
- 238000005381 potential energy Methods 0.000 description 1
- 239000000047 product Substances 0.000 description 1
- 230000001681 protective effect Effects 0.000 description 1
- 230000002441 reversible effect Effects 0.000 description 1
- 239000002904 solvent Substances 0.000 description 1
- 241000894007 species Species 0.000 description 1
- 229910052712 strontium Inorganic materials 0.000 description 1
- CIOAGBVUUVVLOB-UHFFFAOYSA-N strontium atom Chemical compound [Sr] CIOAGBVUUVVLOB-UHFFFAOYSA-N 0.000 description 1
- 230000008093 supporting effect Effects 0.000 description 1
- 239000000725 suspension Substances 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
Images
Classifications
-
- E—FIXED CONSTRUCTIONS
- E21—EARTH OR ROCK DRILLING; MINING
- E21B—EARTH OR ROCK DRILLING; OBTAINING OIL, GAS, WATER, SOLUBLE OR MELTABLE MATERIALS OR A SLURRY OF MINERALS FROM WELLS
- E21B47/00—Survey of boreholes or wells
- E21B47/006—Detection of corrosion or deposition of substances
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N15/00—Investigating characteristics of particles; Investigating permeability, pore-volume or surface-area of porous materials
- G01N15/08—Investigating permeability, pore-volume, or surface area of porous materials
-
- E—FIXED CONSTRUCTIONS
- E21—EARTH OR ROCK DRILLING; MINING
- E21B—EARTH OR ROCK DRILLING; OBTAINING OIL, GAS, WATER, SOLUBLE OR MELTABLE MATERIALS OR A SLURRY OF MINERALS FROM WELLS
- E21B43/00—Methods or apparatus for obtaining oil, gas, water, soluble or meltable materials or a slurry of minerals from wells
- E21B43/16—Enhanced recovery methods for obtaining hydrocarbons
- E21B43/20—Displacing by water
-
- E—FIXED CONSTRUCTIONS
- E21—EARTH OR ROCK DRILLING; MINING
- E21B—EARTH OR ROCK DRILLING; OBTAINING OIL, GAS, WATER, SOLUBLE OR MELTABLE MATERIALS OR A SLURRY OF MINERALS FROM WELLS
- E21B47/00—Survey of boreholes or wells
- E21B47/06—Measuring temperature or pressure
-
- E—FIXED CONSTRUCTIONS
- E21—EARTH OR ROCK DRILLING; MINING
- E21B—EARTH OR ROCK DRILLING; OBTAINING OIL, GAS, WATER, SOLUBLE OR MELTABLE MATERIALS OR A SLURRY OF MINERALS FROM WELLS
- E21B49/00—Testing the nature of borehole walls; Formation testing; Methods or apparatus for obtaining samples of soil or well fluids, specially adapted to earth drilling or wells
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N11/00—Investigating flow properties of materials, e.g. viscosity, plasticity; Analysing materials by determining flow properties
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N13/00—Investigating surface or boundary effects, e.g. wetting power; Investigating diffusion effects; Analysing materials by determining surface, boundary, or diffusion effects
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N15/00—Investigating characteristics of particles; Investigating permeability, pore-volume or surface-area of porous materials
- G01N15/08—Investigating permeability, pore-volume, or surface area of porous materials
- G01N15/082—Investigating permeability by forcing a fluid through a sample
- G01N15/0826—Investigating permeability by forcing a fluid through a sample and measuring fluid flow rate, i.e. permeation rate or pressure change
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
- G16C10/00—Computational theoretical chemistry, i.e. ICT specially adapted for theoretical aspects of quantum chemistry, molecular mechanics, molecular dynamics or the like
-
- E—FIXED CONSTRUCTIONS
- E21—EARTH OR ROCK DRILLING; MINING
- E21B—EARTH OR ROCK DRILLING; OBTAINING OIL, GAS, WATER, SOLUBLE OR MELTABLE MATERIALS OR A SLURRY OF MINERALS FROM WELLS
- E21B2200/00—Special features related to earth drilling for obtaining oil, gas or water
- E21B2200/20—Computer models or simulations, e.g. for reservoirs under production, drill bits
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N11/00—Investigating flow properties of materials, e.g. viscosity, plasticity; Analysing materials by determining flow properties
- G01N2011/006—Determining flow properties indirectly by measuring other parameters of the system
- G01N2011/0093—Determining flow properties indirectly by measuring other parameters of the system thermal properties
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N13/00—Investigating surface or boundary effects, e.g. wetting power; Investigating diffusion effects; Analysing materials by determining surface, boundary, or diffusion effects
- G01N2013/003—Diffusion; diffusivity between liquids
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2113/00—Details relating to the application field
- G06F2113/08—Fluids
Landscapes
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Life Sciences & Earth Sciences (AREA)
- Mining & Mineral Resources (AREA)
- Geology (AREA)
- Chemical & Material Sciences (AREA)
- Health & Medical Sciences (AREA)
- General Health & Medical Sciences (AREA)
- General Physics & Mathematics (AREA)
- Fluid Mechanics (AREA)
- Theoretical Computer Science (AREA)
- Pathology (AREA)
- Immunology (AREA)
- Analytical Chemistry (AREA)
- Biochemistry (AREA)
- General Life Sciences & Earth Sciences (AREA)
- Geochemistry & Mineralogy (AREA)
- Environmental & Geological Engineering (AREA)
- Computing Systems (AREA)
- Dispersion Chemistry (AREA)
- Geophysics (AREA)
- Computer Hardware Design (AREA)
- Bioinformatics & Computational Biology (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
- Testing Resistance To Weather, Investigating Materials By Mechanical Methods (AREA)
Abstract
本发明涉及油田勘探技术领域，公开一种油气层损害类型和程度时空演化4D定量与智能诊断方法及其***。所述方法包括：基于多个因素损害储层的时空演化模拟方程，确定表征多个因素各自损害所述储层的特征参数；及基于表征所述多个因素各自损害所述储层的特征参数，确定表征储层的损害程度的有效特征参数。本发明可定量模拟由各个因素引起的储层损害的特征参数及对储层损害的总特征参数，从而对未发生储层损害的井进行储层损害定量预测和损害规律时空推演，为预防储层损害、制定油藏的开发方案及之后增产措施具有科学指导意义，以及对已损害井也进行储层损害定量诊断和损害规律时空推演，进而优化设计解堵措施、提高或恢复油气井产量和水井注水能力。
Description
技术领域
本发明涉及油田勘探技术领域，具体地涉及一种油气层损害类型和程度时空演化4D定量与智能诊断方法及其***。
背景技术
油田勘探开发的各个时期，由于受到多种内外因素影响，会导致储层原有的物理、化学、热力学和水动力学平衡状态变化，不可避免的使储层近井壁区乃至远井壁区的储层内部渗透率降低，堵塞流体流动，造成储层损害和油井产量下降，甚至“枪毙”储层，使油气井不出一滴油和一口气。造成储层损害的原因是多样的、复杂的，特别是在生产过程中，储层岩石储渗空间、表面润湿性、水动力学场、温度场、岩石种类等不断发生变化，使损害机理随时间而变，且损害周期长、范围宽，损害更具复杂性和叠加性。储层损害一旦发生，必须根据储层损害情况采取相应的解堵措施恢复流体流动通道，以便提高油气井产量和水井注入能力。
通常情况下，导致油气层损害的原因和类型是多种多样的、是多种类型并存的、各类型导致的油气层损害程度是互不相同的，弄清究竟是哪几种或者哪十几种类型导致的油气层损害、每种损害类型导致的损害程度或者每种损害类型对总损害程度的贡献率、每种损害类型的损害带半径和总损害带半径、每种损害类型对油气层损害的时空演化规律等，难度非常大、但非常重要，它不仅对正常生产井避免油气层损害的发生而提高产量(或注入能力)、对已停产井采取解堵措施恢复产量、对面临停产的井采取保护措施延长生产寿命，而且对提高油藏数值模拟精度，准确预测“甜点”等都是不可或缺的，不仅具有较强的理论意义，而且具有很强的实用价值，意义重大。
目前，诊断储层损害的方法可分为矿场诊断法和室内评价法。其中，所述矿场诊断法包括试井法。虽然所述试井法可以定量给出表征待诊断井的预设区域内的储层的损害程度的表皮系数、堵塞比、附加压降等重要参数，但由其表征的表皮系数与其它参数相互联系。也就是说，通过所述试井法得出的表皮系数并不仅仅反映真实储层损害特征，还是各个环节、多因素的综合表现(即其是真实损害表皮系数和由井斜表皮系数、储层形状表皮系数、打开储层不完善表皮系数、非达西流表皮系数、射孔表皮系数等组成的拟表皮系数之和)，必须进行表皮系数分解才能得到真实损害表皮系数。其中，所述室内评价法包括岩心流动实验法。所述岩心流动实验法是通过岩心驱替前后的渗透率变化来了解损害程度大小，虽然比较适合研究单因素储层损害，但难以反映较大尺度上储层损害规律。并且，因室内岩心实验条件比较理想化、评价用岩心都是原始状态岩心、无法考虑储层特性动态变化，使实验结果与井下储层真实损害出入较大。也就是说，迄今为止，国内外历经几十年研究都未对储层损害原因、类型和程度实现准确的定量、快速诊断，更未实现时空演化诊断，成为国际上“增储上产”的关键重大核心难题。
发明内容
本发明的目的是提供一种油气层损害类型和程度时空演化4D定量与智能诊断方法及其***，其可定量模拟由各个相关因素引起的储层损害的特征参数及对所述储层损害的总特征参数，使油气生产井或注入井在勘探开发各个时期或各个环节的油气层损害情况在空间域和时间域上“透明化、数字化、动态化、可视化、智能化”，并具有身临其境的体验，从而对未发生储层损害的井进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义，以及对已损害井也进行储层损害定量诊断和损害规律时空推演，进而优化设计解堵措施、提高或恢复油气井产量和水井注水能力，以及提高油藏数值模拟精度都具有十分重大意义，并使本发明成为油气藏勘探开发过程中不可或缺的技术。
为了实现上述目的，本发明第一方面提供一种确定储层损害程度的方法，所述方法包括：基于多个因素损害储层的时空演化模拟方程，确定表征所述多个因素各自损害所述储层的特征参数，其中，所述储层位于待诊断井的预设区域内；以及基于所述表征所述多个因素各自损害所述储层的特征参数，确定表征所述储层的损害程度的有效特征参数。
通过上述技术方案，本发明创造性地基于多个因素损害储层的时空演化模拟方程，确定表征所述多个因素各自损害所述储层的特征参数；基于所述表征所述多个因素各自损害所述储层的特征参数，确定表征所述储层的损害程度的有效特征参数。由此，通过各个相关因素损害储层的时空演化模拟方程可定量模拟由其各自引起的储层损害的特征参数(例如渗透率)及多个相关因素对所述储层损害的总特征参数(例如总渗透率或有效渗透率)，从而对未发生储层损害的井进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义，以及对已损害井优化设计解堵措施、提高或恢复油井产量和水井注水能力，以及提高油藏数值模拟精度都具有十分重大意义。
本发明第二方面提供一种确定储层损害程度的***，所述***包括：第一参数确定装置，用于基于多个因素损害储层的时空演化模拟方程，确定表征所述多个因素各自损害所述储层的特征参数，其中，所述储层位于待诊断井的预设区域内；以及第二参数确定装置，用于基于所述表征所述多个因素各自损害所述储层的特征参数，确定表征所述储层的损害程度的有效特征参数。
所述确定储层损害程度的***与上述确定储层损害程度的方法相对于现有技术所具有的优势相同，在此不再赘述。
相应地，本发明第三方面还提供一种机器可读存储介质，所述机器可读存储介质上存储有指令，该指令用于使得机器执行所述的确定储层损害程度的方法。
本发明实施例的其它特征和优点将在随后的具体实施方式部分予以详细说明。
附图说明
附图是用来提供对本发明实施例的进一步理解，并且构成说明书的一部分，与下面的具体实施方式一起用于解释本发明实施例，但并不构成对本发明实施例的限制。在附图中：
图1-0是本发明一实施例提供的确定储层损害程度的方法的流程图；
图1-1是本发明一实施例提供的外来固相颗粒损害储层的建模方法的流程图；
图1-2是本发明一实施例提供的确定储层损害程度的方法的流程图；
图2-1是本发明一实施例提供的黏土膨胀损害储层的建模方法的流程图；
图2-2是本发明一实施例提供的储层中的孔隙内的水分子向岩石内部扩散的示意图；
图2-3是本发明一实施例提供的确定黏土膨胀损害储层的时空演化模拟方程的流程图；
图2-4是本发明一实施例提供的确定储层损害程度的方法的流程图；
图2-5是本发明一实施例提供的吸水速率随时间演化的示意图；
图2-6是本发明一实施例提供的表皮系数随时间演化的示意图；
图2-7是本发明一实施例提供的由储层渗透率损害率表征的在第500天黏土膨胀损害储层的半径的示意图；
图3-1是本发明一实施例提供的无机沉淀损害储层的建模方法的流程图；
图3-2是本发明一实施例提供的确定储层损害程度的方法的流程图；
图4-1是本发明一实施例提供的储层内部微粒损害储层的建模方法的流程图；
图4-2是岩石孔隙内表面微粒起动时的受力示意图；
图4-3是本发明一实施例提供的确定储层损害程度的方法的流程图；
图4-4是本发明一实施例提供的表皮系数随时间演化的示意图；
图4-5是本发明一实施例提供的由储层渗透率损害率表征的在第40天储层内部微粒运移损害储层的半径的示意图；
图5-1是本发明一实施例提供的水锁效应损害储层的建模方法的流程图；
图5-2是本发明一实施例提供的建立所述储层的渗透率的分布方程的流程图；
图5-3是本发明一实施例提供的确定孔径为λ的孔隙的体积密度函数的流程图；
图5-4是本发明一实施例提供的确定所述储层的水相饱和度的孔径分布方程的流程图；
图5-5是本发明一实施例提供的渗透率随时间演化的示意图；
图5-6是本发明一实施例提供的表皮系数随时间演化的示意图；
图6-1A是本发明一实施例提供的由毛管束组成的多孔介质的模型示意图；
图6-1B是本发明一实施例提供的在有效应力下的图6-1A所示的毛管束的横截面的示意图；
图6-1C是本发明一实施例提供的在有效应力下的图6-1A所示的某固体基质的纵向示意图；
图6-2是本发明一实施例提供的应力敏感性损害储层的建模方法的流程图；
图6-3是本发明一实施例提供的确定所述储层中的流体的流量的流程图；
图6-4是本发明一实施例提供的表皮系数随时间演化的示意图；
图6-5是本发明一实施例提供的由储层渗透率损害率表征的在第40天应力敏感性损害储层的半径的示意图；
图7-1是本发明一实施例提供的出砂损害储层的建模方法的流程图
图7-2是本发明一实施例提供的确定储层损害程度的方法的流程图；
图7-3是本发明一实施例提供的表皮系数随时间演化的示意图；
图7-4是本发明一实施例提供的由储层渗透率损害率表征的在第40天储层出砂损害储层的半径的示意图；
图8-1是本发明一实施例提供的润湿反转损害储层的建模方法的流程图；
图8-2是本发明一实施例提供的确定储层损害程度的方法的流程图；
图8-3是本发明一实施例提供的表皮系数随时间演化的示意图；
图8-4是本发明一实施例提供的由储层渗透率损害率表征的在第365天润湿反转损害储层的半径的示意图；
图9-1是本发明一实施例提供的乳化堵塞损害储层的建模方法的流程图；
图9-2是本发明一实施例提供的确定储层损害程度的方法的流程图；
图9-3是本发明一实施例提供的表皮系数随时间演化的示意图；
图9-4是本发明一实施例提供的由储层渗透率损害率表征的在第40天乳化堵塞损害储层的半径的示意图；
图10-1是本发明一实施例提供的有机垢损害储层的建模方法的流程图；
图10-2是本发明一实施例提供的储层模型的示意图；
图10-3是本发明一实施例提供的确定有机垢最大溶解量随压力变化的关系式的流程图；
图10-4是本发明一实施例提供的确定储层损害程度的方法的流程图；
图10-5是本发明一实施例提供的表皮系数随时间演化的示意图；
图10-6是本发明一实施例提供的由储层渗透率损害率表征的在第40天有机垢损害储层的半径的示意图；
图11-1是本发明一实施例提供的贾敏效应损害储层的建模方法的流程图；
图11-2是本发明一实施例提供的建立所述储层的渗透率的分布方程的流程图；
图11-3是本发明一实施例提供的确定孔径为λ的孔隙的体积密度函数的流程图；
图11-4是本发明一实施例提供的确定所述储层的水相饱和度的孔径分布方程的流程图；
图11-5是本发明一实施例提供的渗透率随含水饱和度演化的示意图；
图11-6是本发明一实施例提供的含水饱和度随空间演化的示意图；
图11-7是本发明一实施例提供的渗透率损害率随含水饱和度的示意图；
图12-1是本发明一实施例提供的细菌损害储层的建模方法的流程图；
图12-2是本发明一实施例提供的确定所述细菌的生长速率的流程图；
图12-3是本发明一实施例提供的确定储层损害程度的方法的流程图；
图13-1A至图13-1C是本发明一实施例提供的聚合物的层吸附模式的示意图；
图13-1D是本发明一实施例提供的聚合物的桥接吸附模式的示意图；
图13-2是本发明一实施例提供的聚合物损害储层的建模方法的流程图；
图13-3是本发明一实施例提供的建立吸附聚合物的分子链的占比分布方程；
图13-4是本发明一实施例提供的确定储层损害程度的方法的流程图；
图13-5是本发明一实施例提供的表皮系数随时间演化的示意图；
图13-6是本发明一实施例提供的由储层渗透率损害率表征的在第40天聚合物损害储层的半径的示意图；
图14-1至图14-4是本发明一实施例提供的由储层渗透率损害率表征的在第0.6天、第2.4天、第8.4天及第30天储层的损害半径与总表皮系数的示意图；
图14-5是本发明一实施例提供的由储层渗透率损害率表征的在第40天储层的总损害半径的示意图；
图14-6是本发明一实施例提供的总表皮系数随时间演化的示意图；
图14-7是本发明一实施例提供的各类型因素造成的损害在总损害中所占比例随时间演化规律。
具体实施方式
以下结合附图对本发明的具体实施方式进行详细说明。应当理解的是，此处所描述的具体实施方式仅用于说明和解释本发明，并不用于限制本发明。
油田开发过程中会因各种原因造成不同类型的储层损害，使油井产量降低甚至停产、水井注水能力下降，给油田造成巨大经济损失。对已损害储层实施解堵是增加产量、提高采收率的重要措施，但解堵措施的优化决策必须以储层损害类型和程度的定量诊断为前提。因油田开发过程中的储层损害周期长、范围宽，损害更具复杂性和叠加性，给储层损害诊断带来了更大难度，至今国内外缺少准确度高和实用性强的诊断方法和技术，无法定量诊断各损害类型在待诊断井总损害程度中的比例，更无法实现时空演化定量模拟。
针对影响待诊断井周围区域内的储层的渗透率等特征参数的复杂机理及相关参数，本发明各个实施例(例如实施例1-13)创造性地提供了定量模拟各个相关因素(因素1-外来固相颗粒、因素2-黏土膨胀、因素3-无机沉淀、因素4-储层内部微粒运移、因素5-水锁效应、因素6-应力敏感性、因素7- 出砂、因素8-润湿反转、因素9-乳化、因素10-有机垢、因素11-贾敏效应、因素12-细菌、因素13-聚合物)损害储层的时空演化模拟方程，然后通过各个时空演化模拟方程可定量确定对应因素损害所述储层的渗透率等参数，最后结合叠加性原理确定所述带诊断井在不同时空的总损害程度和各个相关因素对总损害程度的贡献。
图1-0是本发明一实施例提供的确定储层损害程度的方法的流程图。如图1-0所示，所述方法可包括：步骤S101，基于多个因素损害储层的时空演化模拟方程，确定表征所述多个因素各自损害所述储层的特征参数，其中，所述储层位于待诊断井的预设区域内；以及步骤S102，基于所述表征所述多个因素各自损害所述储层的特征参数，确定表征所述储层的损害程度的有效特征参数。
具体而言，对于不同类型的待诊断井而言，影响其周围区域内的储层的渗透率等特征参数的相关因素明显不同，如表1所示；并且对于同一类型的待诊断井而言，不同的开发阶段(例如，钻井阶段、注水阶段/注聚阶段、采油阶段)影响其周围区域内的粗曾的渗透率等特征参数的相关因素也明显不同，如表1所示。也就是说，所述多个因素对于不同类型的待诊断井而言是不同的；并且所述多个因素对于处在不同开发阶段的同一类型的待诊断井而言也是不同的。
表1影响待诊断井周围区域内的储层的渗透率的相关因素
对于步骤S102，所述确定表征所述储层的损害程度的有效特征参数可包括：基于所述表征所述多个因素中的第i个因素损害所述储层的特征参数
及下式，确定表征所述储层的损害程度的有效特征参数
其中，Li为
的权重(根据实际情况对其进行合理设置，例如可将其均设置为1)；
为所述储层中的任意位置；t为任意时刻；以及n为所述多个因素的个数。
并且，在获得所述表征所述多个因素中的第i个因素损害所述储层的特征参数
的情况下，基于所述表征所述多个因素中的第i个因素损害所述储层的特征参数
及下式，还可确定
占表征所述储层的损害程度的有效特征参数的比例
其中，所述特征参数可为渗透率、表皮系数和/或渗透率损害率等。在实际应用中，可根据实际需要进行合理设置，可计算上述三个特征参数中的部分或者全部。
在所述待诊断井为注水井、注聚井或采油井并且处于钻井阶段的情况下，所述多个因素包括以下各项中的至少两者：外来固相颗粒、黏土膨胀、储层内部微粒运移、无机沉淀与水锁效应。在实际应用中，可根据实际需要进行合理设置。
在一实施例中，在所述待诊断井为注水井、注聚井或采油井并且处于钻井阶段的情况下，所述多个因素可包括外来固相颗粒、黏土膨胀与储层内部微粒运移。
对于步骤S102，基于所述表征外来固相颗粒、黏土膨胀与储层内部微粒运移各自损害所述储层的特征参数
(其下角标对应于表1中的钻井阶段中的第1个因素-外来固相颗粒、第2个因素-黏土膨胀和第4个因素-储层内部微粒运移，并该3个特征参数可分别通过以下实施例1、2、4所示的时空演化模拟方程)及下式，确定表征所述储层的损害程度的有效特征参数
在另一实施例中，所述多个因素还可包括以下各项中的至少一者：无机沉淀与水锁效应。例如，所述多个因素还可包括无机沉淀与水锁效应两者。
对于步骤S102，基于所述表征外来固相颗粒、黏土膨胀、储层内部微粒运移、无机沉淀与水锁效应各自损害所述储层的特征参数
(其可分别通过以下实施例1、2、4、3、5所示的时空演化模拟方程)及下式，确定表征所述储层的损害程度的有效特征参数
在所述待诊断井为采油井并且处于采油阶段的情况下，所述多个因素包括以下各项中的至少两者：储层内部微粒运移、出砂、乳化、贾敏效应、应力敏感性、润湿反转与有机垢。在实际应用中，可根据实际需要进行合理设置。
在一实施例中，在所述待诊断井为采油井并且处于采油阶段的情况下，所述多个因素可包括储层内部微粒运移、出砂、乳化与有机垢。
对于步骤S102，基于所述表征储层内部微粒运移、出砂、乳化与有机垢各自损害所述储层的特征参数
(其下角标对应于表1中的采油阶段中的第1个因素-储层内部微粒运移、第3个因素-出砂、第5个因素-乳化和第6个因素-贾敏效应，并该4个特征参数可分别通过以下实施例4、7、9、10所示的时空演化模拟方程)及下式，确定表征所述储层的损害程度的有效特征参数
在另一实施例中，所述多个因素还可包括以下各项中的至少一者：应力敏感性、润湿反转与有机垢。例如，所述多个因素还可包括应力敏感性、润湿反转与贾敏效应三者。
对于步骤S102，基于所述表征储层内部微粒运移、出砂、乳化、贾敏效应、应力敏感性、润湿反转与有机垢各自损害所述储层的特征参数
(其可分别通过以下实施例4、7、9、10、6、8、11所示的时空演化模拟方程)及下式，确定表征所述储层的损害程度的有效特征参数
在所述待诊断井为注水井并且处于注水阶段的情况下，所述多个因素包括以下各项中的至少两者：黏土膨胀、细菌、水锁效应、外来固相颗粒、微粒运移与无机沉淀。在实际应用中，可根据实际需要进行合理设置。
在一实施例中，在所述待诊断井为注水井并且处于注水阶段的情况下，所述多个因素可包括黏土膨胀、细菌与水锁效应。
对于步骤S102，基于所述表征黏土膨胀、细菌与水锁效应各自损害所述储层的特征参数
(其下角标对应于表1中的注水阶段中的第2个因素-黏土膨胀、第5个因素-水锁效应、第6个因素-细菌，并该3个特征参数可分别通过以下实施例2、12、5所示的时空演化模拟方程)及下式，确定表征所述储层的损害程度的有效特征参数
在另一实施例中，所述多个因素还可包括以下各项中的至少一者：外来固相颗粒、储层内部微粒运微粒运移与无机沉淀。例如，所述多个因素还可包括外来固相颗粒、储层内部微粒运微粒运移与无机沉淀三者。
对于步骤S102，基于所述表征黏土膨胀、细菌、水锁效应、外来固相颗粒、储层内部微粒运微粒运移与无机沉淀各自损害所述储层的特征参数
(其可分别通过以下实施例2、12、5、1、4、3所示的时空演化模拟方程)及下式，确定表征所述储层的损害程度的有效特征参数
在所述待诊断井为注聚井并且处于注聚阶段的情况下，所述多个因素包括以下各项中的至少两者：聚合物、黏土膨胀、外来固相颗粒、微粒运移与无机沉淀。在实际应用中，可根据实际需要进行合理设置。
在一实施例中，在所述待诊断井为注聚井并且处于注聚阶段的情况下，所述多个因素可包括聚合物与黏土膨胀。
对于步骤S102，基于所述表征聚合物与黏土膨胀各自损害所述储层的特征参数
(其下角标对应于表1中的注聚阶段中的第2个因素-黏土膨胀、第5个因素-聚合物，并该2个特征参数可分别通过以下实施例13、2所示的时空演化模拟方程)及下式，确定表征所述储层的损害程度的有效特征参数
在另一实施例中，所述多个因素还可包括以下各项中的至少一者：外来固相颗粒、储层内部微粒运微粒运移与无机沉淀。例如，所述多个因素还可包括外来固相颗粒、储层内部微粒运微粒运移与无机沉淀三者。
对于步骤S102，基于所述表征聚合物、黏土膨胀、外来固相颗粒、储层内部微粒运微粒运移与无机沉淀各自损害所述储层的特征参数
(其可分别通过以下实施例2、 13、1、4、3所示的时空演化模拟方程)及下式，确定表征所述储层的损害程度的有效特征参数
在实际应用中，可根据实际需要对所述多个因素进行合理设置，即并不限于上述各个实施例的组合。
下面分别介绍建立上述表1中的13种相关因素中的每个因素损害储层的时空间演化模拟方程及通过各个时空间演化模拟方程确定相关特征参数的过程(详见实施例1-实施例13，该13个实施例分别对应于表1中的13 个因素)。针对不同阶段中所涉及的同一因素，可通过相应实施例确定该同一因素损害储层的渗透率等参数；只不过针对不同阶段所述相应实施例中的某些参数的具体数值可能会不同。
需要说明的是，为了简单描述起见，在本发明的各个实施例中的随时空演化的物理量、化学量可省略变量
例如
可简写为K。
实施例1-外来固相颗粒
外来固相颗粒侵入储层并产生堵塞的本质是外源性固相颗粒侵入介质之后的运移和沉积。因此，本发明各实施例的核心是要建立固相颗粒运移和沉积的动力学模型。具体地，基于质量守恒、扩散关系等建立外来固相颗粒在待诊断井周围的储层中的浓度分布的时空演化控制唯象模型(该模型包含流动颗粒的浓度C和沉积颗粒的浓度Cd)，再结合沉积浓度和渗透率等储层损害特征参数间的关系，就可诊断渗透率等储层损害特征参数的时空场分布。
图1-1是本发明一实施例提供的外来固相颗粒损害储层的建模方法的流程图。所述建模方法可包括步骤S1101-S1104。
步骤S1101，确定储层中的含有流动颗粒的流体的速度。
其中，所述储层位于待诊断井(例如，注水井、采油井等)的预设区域内。
对于步骤S1101，所述确定储层中的含有流动颗粒的流体的速度可包括：建立所述流体进入所述储层的压力传导方程；以及根据所述压力传导方程及达西公式，确定所述流体的速度。
具体地，压力是驱动固-液混合液(即含有流动颗粒的流体)从注水井的井筒持续侵入待诊断井周围的储层的动力，由此可建立如公式(1-1)的所述流体进入储层的压力传导方程：
再根据公式(1-1)及达西公式(如下式(1-2))可确定所述流体的速度，
其中，
为所述流体的压力；φ为所述储层的孔隙度；μ为流体粘度；ct为流体-岩石综合压缩系数；
为所述储层的渗透率；以及τ为所述储层的迂曲度。
步骤S1102，基于所述流体的对流参数与扩散参数，建立所述流体与所述储层上的沉积颗粒之间的质量平衡方程。
根据质量守恒，所述流体的质量的变化等于沉积颗粒的质量的负变化。对于步骤S1102，所述建立所述流体与所述储层上的沉积颗粒之间的质量平衡方程可包括：基于所述流体的对流参数与扩散参数，建立下式(1-3)表示的所述质量平衡方程，
其中，ρ为所述流体的密度；φ为所述储层的孔隙度；
为所述流动颗粒的质量分数(也可称为质量浓度)；u为达西表观速度；
为扩散流量，
其中ρL为所述流体的密度，
为扩散系数，
α为垂直扩散率，
为所述流体的速度；
为单位时间内的所述沉积颗粒的累积质量；t为时间。
步骤S1103，基于所述流体的对流参数与扩散参数，建立所述沉积颗粒的体积浓度与所述流体的体积浓度之间的连接条件方程。
对于步骤S1103，所述建立所述沉积颗粒的体积浓度与所述流体的体积浓度之间的连接条件方程可包括：基于所述流体的对流参数与扩散参数，建立下式(1-4)表示的所述连接条件方程，
其中，ρp为颗粒密度；
为所述沉积颗粒的体积浓度；
其中k0为原始滤失系数，
以及
由于F1(T)和温度的相关性是由exp(1/T)来度量的，在常见的温度范围(比如300K～400K)内，这个函数的变化其实是很缓慢的，实际上接近等温过程，故
其中
为所述沉积颗粒的体积浓度，Cdmax为所述沉积颗粒的最大体积浓度，以及mk为第一经验值。上述所有参数都既可以是常数，也可以是随空间变化的参数，即非均质的情况。
步骤S1104，根据所述流动颗粒的质量分数与所述流动颗粒的体积浓度之间的关系、所述流体的速度、所述质量平衡方程及所述连接条件方程，确定所述颗粒损害储层的时空演化模拟方程。
其中，所述时空演化模拟方程用于模拟由外来固相颗粒引起的储层损害特征的四维时空演化过程。
其中，所述流动颗粒的质量分数与所述流动颗粒的体积浓度之间的关系可为
其中ρp为所述沉积颗粒的密度；ρL为所述流体的密度；
为所述流动颗粒的体积浓度。所述颗粒损害储层的时空演化模拟方程可包括：公式(1-5)所示的颗粒运移损害储层的时空演化模拟方程，以及公式(1-6)所示的颗粒沉积损害储层的时空演化模拟方程。
对于步骤S1104，所述确定所述颗粒损害储层的时空演化模拟方程可包括：根据所述流动颗粒的质量分数与所述流动颗粒的体积浓度之间的关系、所述流体的速度及公式(1-3)所示的质量平衡方程，确定由下式(1-5)所示的颗粒运移损害储层的时空演化模拟方程：
以及根据所述流动颗粒的质量分数与所述流动颗粒的体积浓度之间的关系、所述流体的速度及公式(1-4)所示的连接条件方程，确定由公式(1-6)所示的颗粒沉积损害储层的时空演化模拟方程：
其中，
为所述流动颗粒的体积浓度；
为所述流体的速度；τ为所述储层的迂曲度；ρp为所述沉积颗粒的密度；ρL为所述流体的密度；
所述储层的滤失系数的初始值；
为所述沉积颗粒的体积浓度；Cdmax为所述沉积颗粒的最大体积浓度；mk为第一经验值；α为垂直扩散率；以及φ为所述储层的孔隙度。其中
其中 NR、NPe、NA、NDL、NE1、NE2、NG、NLo、NvdW、ζp(g)分别位半径数、皮克雷数、吸引数、双电层数、第一电位力数、第二电位力数、重力数、伦敦力数、范德瓦尔斯力数及流动颗粒和基质颗粒的电势(各个参数的相关表达式详见表2)；
表2包含固相沉积驱动因素的无量纲参数表及其表达式
注：D∞为颗粒自由扩散率。H为Hamaker数。Dp、Dg分别为流动颗粒直径与基质颗粒直径。μ为流体粘度。kB为Boltzmann常数。ζp、ζg分别为流动颗粒与基质颗粒的电势。
综上所述，本发明创造性地建立含有流动颗粒的流体与储层上的沉积颗粒之间的质量平衡方程；建立所述沉积颗粒的体积浓度与所述流体的体积浓度之间的连接条件方程；根据所述流动颗粒的质量分数与所述流动颗粒的体积浓度之间的关系、所述流体的速度、所述质量平衡方程及所述连接条件方程，确定所述颗粒损害储层的时空演化模拟方程。由此，通过所确定的时空演化模拟方程可定量模拟由外来固相颗粒引起的储层损害特征的四维时空演化过程，从而对已损害井优化设计解堵措施、提高或恢复油井产量和水井注水能力，以及提高油藏数值模拟精度都具有十分重大意义。
图1-2是本发明一实施例提供的确定储层损害程度的方法的流程图。如图1-2所示，所述确定储层损害程度的方法可包括步骤S1201-S1202。
步骤S1201，基于所述的外来固相颗粒损害储层的建模方法建立的时空演化模拟方程，确定所述沉积颗粒的体积浓度。
对于上述公式(1-5)所示的颗粒运移损害储层的时空演化模拟方程，在一维情形下，该类方程可以整理为以下的一般形式：
其中，aa,bb,cc可为常数(如扩散系数)，也可为函数(如所述流体的速度)；f可为压力、物质浓度、应力等。对时间采用向后差分，空间采用中心差分。则上述方程可以有如下差分格式：
其中，i＝1,2,3...Ni，
n＝1,2,3...，t＝nΔt，Ni为离散空间点个数。
求解区间为x∈(0,xmax)，Δx、Δt为空间、时间步长。同时，考虑初始条件
和边界条件(fi n|i＝1＝f0,n＝1,2,3...(井壁处)以及
)(构造了一个虚拟网格i+1，预设范围的边界处或距井壁数米处)。
首先，对于i＝2,3,...,Ni-1整理上述差分格式有：
其中，A1i，A2i，A3i分别为，
同时，据公式(1-5)可知：
将公式(1-11)代入公式(1-10)可得到迭代关系式(1-9)的具体表现形式，由于该迭代关系式(1-9)的具体表现形式复杂，故在此不对其进行列出。然后，利用初始条件和边界条件进行迭代计算就可得到场f的值。
接着，对说明边界条件的差分求解过程进行说明。
上述迭代关系式(1-9)适用于非边界网格。而对于i＝1(井壁处)而言，因为采用的是点中心网格，且其为狄利克雷(Dirichlet)边界条件，故直接可得到以下关系式：
f1 n＝f0(常数)，i＝1 (1-12)
对于i＝N(预设范围的边界处距井壁数米处)而言，其为诺伊曼或第二类(Neumann)边界条件，增加一个虚拟网格i＝Ni+1，由
知
将其代入式(1-9)可知：
根据上述过程可求解出场函数f的时空变化情况。由于上述数值模型是针对待诊断井的井筒附近储层建立的，在求解某物理量f在井周的分布时，需要采用柱坐标系。由此，式
需要变换为
这种形式不利于等距差分，可以引入坐标变换： r＝rwex′，其中，rw为井筒半径，x′为一个无量纲的空间坐标。将这个变换代入一般方程中，可以得到关于x′的方程：
如果将
和
作为新的方程系数，则上式和
相比，本质上是一样的。因此，便可以在x′坐标进行等距差分并沿用前述的迭代格式。计算完f的值后，再将空间坐标从x′映射回r即可得到f(r,t)。
通过上述外来固相颗粒损害储层的建模方法建立的时空演化模拟方程综合考虑了固相颗粒侵入储层时多种物理化学因素对储层损害的影响，由此通过该步骤S1201求解得到的沉积颗粒的体积浓度非常精确。
步骤S1202，基于所述沉积颗粒的体积浓度，确定表征待诊断井的预设区域内的储层的损害程度的特征参数。
其中，所述特征参数可为所述储层的渗透率和/或所述储层的滤失系数。
在一实施例中，所述特征参数可为所述储层的渗透率。
对于步骤S1202，所述确定表征待诊断井的预设区域内的储层的损害程度的特征参数可包括：基于所述沉积颗粒的体积浓度
及公式(1-15)，确定所述储层的渗透率
在一实施例中，所述特征参数可为所述储层的滤失系数。
对于步骤S1202，所述确定表征待诊断井的预设区域内的储层的损害程度的特征参数可包括：基于所述沉积颗粒的体积浓度
及公式(1-16)，确定所述储层的滤失系数
其中，φ0为孔隙度的初始值；Cdmax为所述沉积颗粒的最大体积浓度；mk与mK分别为第一经验值与第二经验值；
为所述储层的渗透率的初始值；以及
所述储层的滤失系数的初始值。
其中，所述特征参数可为所述储层的表皮系数。
对于步骤S1202，所述确定表征待诊断井的预设区域内的储层的损害程度的特征参数可包括：基于所述沉积颗粒的体积浓度
及公式
确定所述储层的渗透率
以及基于所述储层的渗透率
及公式(1-17)，确定所述储层的表皮系数
其中，
为所述储层的渗透率的初始值；以及
通过该步骤S1202得到的特征参数(例如所述储层的渗透率
与表皮系数
)是时空演化4D定量模拟的结果，可根据其演化特点进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义。
综上所述，通过所确定的时空演化模拟方程可确定所述沉积颗粒的体积浓度，再根据所述沉积颗粒的体积浓度可确定表征待诊断井的预设区域内的储层的损害程度的特征参数(例如所述储层的渗透率和/或表皮系数)，由此可定量模拟由外来固相颗粒引起的储层损害特征的四维时空演化过程，从而对已发生储层损害的井，利用历史参数实现储层损害定量模拟和时空演化，为解堵措施优化设计和提高油藏数值模拟精度具有重要意义；对未发生储层损害的井，利用物性参数和将实施的工程参数，可进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义。
实施例2-黏土膨胀
外来流体(例如注入水)的水分子通过固-液(储层中的岩石与流体之间的)界面向固相介质(储层中的岩石)内部的扩散过程可认为是固-液界面的任意局部区域的情况，对于该界面上的任意足够小尺度的局域区域而言，可认为水分子的扩散方向垂直于该区域的某点(例如O点)的切线方向(即扩散方向与该区域所在的平面相垂直)，如图2-2(其中阴影部分表示岩石，其他空白部分表示储层中的孔隙)所示。所述储层中的岩石包括黏土，在水分子向岩石扩散的过程中，所述黏土会发生膨胀，进而可能会导致储层的渗透率降低(甚至堵塞)。因此，本发明各实施例的核心是要建立水分子向岩石内部扩散和储层中的孔隙内的水含量变化的动力学模型(即水分子通过固 -液界面由孔隙中的液相向固相内部扩散的扩散方程与孔隙中的流体的对流扩散方程)。具体地，基于菲克扩散定律以及储层中的孔隙内的流体的对流扩散关系等建立黏土膨胀影响待诊断井周围的储层中的孔隙度分布的时空演化控制唯象模型(该模型包含储层内的孔隙的含水体积分数c1与储层中的岩石的含水体积分数的初始值c0)，再结合储层的孔隙度和渗透率等储层损害特征参数间的关系，就可诊断渗透率等储层损害特征参数的时空场分布。
图2-1是本发明一实施例提供的黏土膨胀损害储层的建模方法的流程图。所述建模方法可包括步骤S2101-S2104。
步骤S2101，确定待诊断井的预设区域内的储层中的流体的达西表观速度。
其中，所述待诊断井可例如为注水井或采油井。
对于步骤S2101，所述确定储层中的流体的速度可包括：建立所述流体进入所述储层的压力传导方程；以及根据所述压力传导方程及达西公式，确定所述流体的达西表观速度。
具体地，压力是驱动固-液混合液从注水井的井筒持续侵入周围的储层的动力，由此可建立如公式(2-1)的所述流体进入储层的压力传导方程：
再根据公式(2-1)及达西公式(如下式(2-2))可确定所述流体的达西表观速度，
其中，
为所述流体的压力；φ为所述储层的孔隙度；μ为流体粘度；ct为流体-岩石综合压缩系数以及
为所述储层的渗透率。
步骤S2102，根据所述流体的达西表观速度与所述流体中的水分子的扩散系数，建立所述流体中的水分子的质量平衡方程。
在储层条件下，所述储层中的孔隙内的不同位置的含水量满足质量守恒方程。其中，所述储层内的含水量变化主要由对流与扩散两个过程决定。具体地，对于步骤S2102，所述建立所述流体中的水分子的质量平衡方程可包括：根据所述流体的达西表观速度u与所述水分子的扩散系数Dw，建立下式 (2-3)表示的所述质量平衡方程：
其中，φ0为所述储层的孔隙度的初始值；
为所述储层内的孔隙的含水体积分数；以及
为所述储层内的任意点的空间位置。
所述流体中的水分子的质量平衡方程的初始条件为
以及所述流体中的水分子的质量平衡方程的边界条件为
(也就是说，在注水井的井壁处的储层孔隙完全被水充满，即孔隙中的水饱和度为1)。其中，φ0为所述储层的孔隙度的初始值；rw为所述待诊断井的井筒半径；以及Swc为所述储层中的束缚水饱和度。
步骤S2103，根据菲克扩散定律，建立所述流体中的水分子向储层中的岩石内部扩散的扩散方程。
需要说明的是，c(n，t)是储层中的岩石在时刻t、在以位置
(例如图2-2 中的O点)为原点建立的一维坐标系(坐标轴的方向为固-液界面在该位置
处的法线方向)中的坐标n处的含水体积分数；相应地，储层中的孔隙在位置
处的含水体积分数为
可建立垂直固-液界面指向固相内部的一维坐标系n，界面处有n＝0，固相内部n>0，如图2-2所示。对于步骤S2103，所述建立所述流体中的水分子向储层中的岩石内部扩散的扩散方程可包括：根据所述菲克扩散定律，建立下式(2-4)表示的所述流体中的水分子向储层内部扩散的扩散方程：
其中，n为以储层中的位置
为原点且以所述流体与所述岩石的界面在位置
处的法线方向为坐标轴的方向建立的一维坐标系中的坐标；t为时间；Dw为所述水分子的扩散系数；以及c(n,t)为所述储层中的岩石的含水体积分数。
其中，所述扩散方程的初始条件为c(n,t＝0)＝c0；以及所述扩散方程的边界条件为
Dw为所述水分子的扩散系数；c(n,t)为所述储层中的岩石的含水体积分数；以及kf为膜交换系数。
步骤S2104，根据所述扩散方程及所述质量平衡方程，确定黏土膨胀损害储层的时空演化模拟方程。
对于步骤S2104，如图2-3所示，所述确定黏土膨胀损害储层的时空演化模拟方程可包括步骤S2301-S2303。
步骤S2301，根据所述质量平衡方程、该质量平衡方程的边界条件与初始条件，确定所述储层内的孔隙的含水体积分数。
根据上述公式(2-3)及所述流体中的水分子的质量平衡方程的初始条件与边界条件，可求解得到所述储层内的孔隙的含水体积分数
步骤S2302，根据所述扩散方程、所述扩散方程的边界条件、所述扩散方程的初始条件及所述储层内的孔隙的含水体积分数，确定所述储层中的岩石的吸水速率。
首先，根据上述公式(2-4)及所述扩散方程的边界条件、所述扩散方程的初始条件，可求解得到所述储层中的岩石的含水体积分数c(n，t)，
其中，
即余误差函数，L-1{·}表示逆拉普拉斯变换。
然后，根据
c(n，t)及所述储层中的岩石的吸水速率的定义
可求解得到
具体地，c1、c0的相对大小决定了
的正负。若c1＞c0，则
说明孔隙中水分含量大于固相中水分含量，将向固相内扩散；反之，若c1＜c0，两式将给出
这意味着岩石固相失水。最后，在储层条件下，孔隙中的含水量总是要大于或等于固相中的含水量，需对
加限制条件，使其在c1＜c0时等于0。因此，
表达式为：
步骤S303，根据所述储层中的岩石的吸水速率，确定所述黏土膨胀损害储层的时空演化模拟方程。
对于步骤S103，所述确定所述黏土膨胀损害储层的时空演化模拟方程可包括：根据所述储层的吸水速率
确定由下式(2-6)表示的所述黏土膨胀损害储层的时空演化模拟方程：
其中，
为所述储层的孔隙度；以及λ为黏土膨胀系数。
具体地，黏土膨胀系数
其中Cc为岩石中的黏土质量百分数；PI为岩石的塑性系数(无量纲)，若PI＜1～2，该岩石为脆性岩石，若2＜PI＜6，该岩石为塑脆性岩石，若PI＞6，该岩石为塑性岩石；k′为经验参数。
根据上述公式(2-6)可求解得到
若
(即吸水速率为正)，黏土就会膨胀，使得
即孔隙度下降。
综上所述，本发明创造性地确定待诊断井的预设区域内的储层中的流体的达西表观速度；根据所述流体的达西表观速度与所述流体中的水分子的扩散系数，建立所述流体中的水分子的质量平衡方程；建立所述流体中的水分子向储层中的岩石内部扩散的扩散方程；以及根据所述扩散方程及所述质量平衡方程，确定黏土膨胀损害储层的时空演化模拟方程。由此，通过所确定的时空演化模拟方程可定量模拟由黏土膨胀引起的储层损害特征的四维时空演化过程，从而对未发生储层损害的井进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义，以及对已损害井优化设计解堵措施、提高或恢复油井产量和水井注水能力，以及提高油藏数值模拟精度都具有十分重大意义。
图2-4是本发明一实施例提供的确定储层损害程度的方法的流程图。如图2-4所示，所述确定储层损害程度的方法可包括步骤S2401-S2402。
步骤S2401，基于所述的黏土膨胀损害储层的建模方法建立的时空演化模拟方程，确定所述储层的孔隙度。
对于上述公式(2-6)所示的黏土膨胀损害储层的时空演化模拟方程的求取，需要根据公式(2-3)计算得到
具体求解过程可参见上述实施例1求解沉积颗粒的体积浓度的过程，于此不再进行赘述。
在通过上述方法计算得到所述储层内的孔隙的含水体积分数c1(r，t)之后，再根据上述公式(2-5)可计算得到所述储层的吸水速率
(图2-5示出了储层中的特定位置r处的吸水速率
随时间变化情况)，由此通过上述黏土膨胀损害储层的建模方法建立的时空演化模拟方程综合考虑了黏土膨胀时多种物理化学因素对储层损害的影响，由此通过该步骤S2401求解得到的储层的孔隙度非常精确。
步骤S2402，基于所确定的所述储层的孔隙度，确定表征待诊断井的预设区域内的储层的损害程度的特征参数。
其中，所述特征参数可为所述储层的渗透率。
对于步骤S2402，所述确定表征待诊断井的预设区域内的储层的损害程度的特征参数可包括：基于所述储层的孔隙度
及公式(2-7)，确定所述储层的渗透率
其中，φ0为孔隙度的初始值；mK为第二经验值；以及
为所述储层的渗透率的初始值。
其中，所述特征参数可为所述储层的滤失系数。
对于步骤S2402，所述确定表征待诊断井的预设区域内的储层的损害程度的特征参数可包括：基于所述储层的孔隙度
及公式(2-8)，确定所述储层的滤失系数
其中，φ0为孔隙度的初始值；φdmax为所述储层的最大孔隙度；mk为第一经验值；以及
所述储层的滤失系数的初始值。
在所述特征参数为所述储层的渗透率与所述储层的滤失系数的情况下，可通过公式(2-7)确定所述储层的渗透率，并通过公式(2-8)所述储层的渗透率所述储层的渗透率。
其中，所述特征参数可为所述储层的表皮系数。
对于步骤S2402，所述确定表征待诊断井的预设区域内的储层的损害程度的特征参数可包括：基于所述储层的孔隙度φ(r,t及公式
确定所述储层的渗透率
以及基于所述储层的渗透率
及公式(2-9)，确定所述储层的表皮系数
其中，
为所述储层的渗透率的初始值；以及
rw为所述待诊断井的井筒半径，以及rsw为所述储层的损害半径。
通过该步骤S2402得到的特征参数(例如所述储层的渗透率
与表皮系数
)是时空演化4D定量模拟的结果(如图2-6所示)。更具体地，图2-7示出了由储层渗透率损害率(基于所述储层的渗透率
及公式
确定所述储层的渗透率损害率I(ri,t)，其中
为
的最大值)表征的在第40天黏土膨胀损害储层的半径的示意图(如箭头所指示的半径)，相关工作人员可通过该图2-7直观地确认储层被损害的程度。因此，可根据渗透率或表皮系数的演化特点进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义。
综上所述，本发明创造性地通过所确定的时空演化模拟方程可确定所述储层的孔隙度，再根据所述储层的孔隙度可确定表征待诊断井的预设区域内的储层的损害程度的特征参数(例如所述储层的渗透率和/或表皮系数)，由此可定量模拟由黏土膨胀引起的储层损害特征的四维时空演化过程，从而对未发生储层损害的井进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义，以及对已损害井优化设计解堵措施、提高或恢复油井产量和水井注水能力，以及提高油藏数值模拟精度都具有十分重大意义。
实施例3-无机沉淀
当外来流体与储层内的流体不配伍时，可能会生成无机沉淀。该无机沉淀会堵塞流体流动通道，从而造成储层损害。通常的无机沉淀可包括：碳酸钙(CaCO3)、硫酸钙(CaSO4)、硫酸锶(SrSO4)、硫酸钡(BaSO4)等无机沉沉淀物。
图3-1是本发明一实施例提供的无机沉淀损害储层的建模方法的流程图。如图3-1所示，所述建模方法可包括以下步骤S3101-S3104。
步骤S3101，确定待诊断井(例如，注水井或采油井)的预设区域内的储层中的流体的达西表观速度。
对于步骤S3101，所述确定储层中的流体的速度可包括：建立所述流体进入所述储层的压力传导方程；以及根据所述压力传导方程及达西公式，确定所述流体的达西表观速度。
具体确定过程可参见上述实施例2中的确定达西表观速度的过程(即上式(2-1)与(2-2)及其相关描述)。
步骤S3102，确定与外来流体中的多个离子中的各个离子相对应的离子浓度损失。
其中，所述离子浓度损失是由所述各个离子与所述储层中的流体中的相应离子发生沉淀反应引起的。
对于步骤S3102，所述确定与外来流体中的多个离子中的各个离子相对应的离子浓度损失可包括：根据所述各个离子i的浓度
与所述各个离子i发生沉淀反应的至少一个离子中的离子j的浓度
及下式(3-1)，确定与所述各个离子i相对应的离子浓度损失
其中，kij为所述各个离子i与所述离子j的反应速率；Ni为与所述各个离子 i发生沉淀反应的至少一个离子的个数；以及βi、βj分别为离子i与离子j 的反应系数。
以外来流体中的
为例，该
可与储层内的流体中的Ca2+、Ba2+与 Sr2+发生反应而产生沉淀物质。t时刻，在储层内的任一空间位置
处，因外来流体中的
与所述储层内的流体中的Ca2+、Ba2+与Sr2+发生沉淀反应而产生的离子的浓度损失
其中，
分别为
与Ca2+、Ba2+与Sr2+发生沉淀反应的反应速率；
分别为t时刻在储层空间
处的离子
Ca2+、Ba2+与Sr2+的浓度；
分别为离子
与离子Ca2+的反应系数；离子
与离子Ba2+的反应系数；以及离子
与离子Sr2+的反应系数。类似地，可确定由外来流体中的其他各个离子与储层内的相应离子发生沉淀反应而导致的储层内的离子的浓度损失。
具体地，所述各个离子i与所述离子j的反应速率kij由所述各个离子i 与所述离子j产生的相应沉淀物的结垢指数决定。例如，所述各个离子i与所述离子j的反应速率kij可满足以下关系式(3-2)：
其中，kij0为常数；以及ISij为所述各个离子i与所述离子j产生的相应沉淀物的结垢指数。更为具体地，所述结垢指数ISij可由所述各个离子i的浓度、所述流体的离子强度、所述离子j的浓度、所述流体的温度及所述流体的压力决定。
具体而言，在一定的温度T、压力P及离子浓度([Me]为(自由)阳离子的浓度，其中，阳离子可为钙离子、锶离子或钡离子等；[An]为(自由) 阴离子的浓度，其中，阴离子可为碳酸氢根离子或硫酸根离子等)下，储层溶液中某一沉淀反应能否发生通常采用结垢指数来判断，其表达式为：
其中，Si为所述流体的离子强度；以及KC为沉淀反应的溶度积系数。
若IS≤0，溶液处于欠饱和或者饱和状态，无机沉淀物生成；若IS＞0，溶液处于过饱和状态，有生成无机沉淀物的趋势。Is随着储层位置、离子浓度、温度、压力的变化而变化，是一个与时间、空间有关的函数。
根据Tomson-Oddo计算法，在储层条件下，四种无机沉淀物的结垢指数分别为：
①碳酸钙CaCO3：
②硫酸钙CaSO4：
a、T＜80℃时，形成的沉淀主要是CaSO4·2H2O，其结垢指数表达式为：
b、80℃＜T＜121℃时，形成的沉淀主要是CaSO4·1/2H2O，其结垢指数表达式为：
c、T＞121℃时，形成的沉淀主要是CaSO4，其结垢指数表达式为：
③硫酸钡BaSO4：
④硫酸锶SrSO4：
上述各式中的PH为储层内的液体(储层内的原液体与外来液体混合而成的液体)的PH值。
步骤S3103，根据所述流体的达西表观速度、与所述各个离子相对应的离子浓度损失及所述各个离子的扩散系数，建立所述各个离子的运移方程。
在离子浓度控制方程中，离子运移的流量J包含对流和扩散两部分，以离子i为例：
其中，Ji为离子i的运移流量；Jid为离子i的扩散流量；Jic为离子i的对流流量；u为所述流体的达西表观速度；Ci为所述各个离子i的浓度；以及Di为所述各个离子i的扩散系数。
对于步骤S3103，沉淀反应会导致在时刻t的某储层空间
处的离子浓度损失
根据质量守恒定律可建立所述各个离子i的运移方程可包括：
其中，u为所述流体的达西表观速度；
为所述储层的孔隙度；以及Di为所述各个离子i的扩散系数。
步骤S3104，根据所述各个离子的运移方程及所述各个离子所产生的沉淀物的反应系数，确定所述无机沉淀损害储层的时空演化模拟方程。
其中，所述时空演化模拟方程用于模拟由所述多个离子产生的相应沉淀引起的储层损害特征的四维时空演化过程。
对于步骤S104，所述确定所述无机沉淀损害储层的时空演化模拟方程可包括：根据所述各个离子的质量平衡方(3-4)，确定所述各个离子i的浓度
以及根据所述各个离子i的浓度
及所述各个离子所产生的沉淀物的反应系数coij，确定下式(3-5)表示的所述无机沉淀损害储层的时空演化模拟方程，
其中，
为在时刻t且在储层空间
处，由离子i与离子j发生沉淀产生的沉淀物的浓度累积量；以及
为在时刻t+dt且在储层空间
处，由离子i与离子j发生沉淀产生的沉淀物的浓度累积量。
综上所述，本发明创造性地根据待诊断井的预设区域内的储层中的流体的达西表观速度、与外来流体中的多个离子中的各个离子相对应的离子浓度损失及所述各个离子的扩散系数，建立所述各个离子的运移方程；以及根据所述各个离子的运移方程及所述各个离子所产生的沉淀物的反应系数，确定所述无机沉淀损害储层的时空演化模拟方程。由此，通过所确定的时空演化模拟方程可定量模拟由无机沉淀引起的储层损害特征的四维时空演化过程，从而对未发生储层损害的井进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义，以及对已损害井优化设计解堵措施、提高或恢复油井产量和水井注水能力，以及提高油藏数值模拟精度都具有十分重大意义。
图3-2是本发明一实施例提供的确定储层损害程度的方法的流程图。如图3-2所示，所述方法可包括步骤S3201-S3202。
步骤S3201，基于所述的无机沉淀损害储层的建模方法建立的时空演化模拟方程及储层内产生的各个沉淀物的摩尔质量与密度，确定无机沉淀的体积浓度。
其中，所述各个沉淀物由外来流体中的多个离子中的各个离子与储层中的流体中的相应离子发生沉淀反应而产生，以及所述无机沉淀的体积浓度为所述各个沉淀物的总体积浓度。
对于上述公式(3-4)所示的无机沉淀损害储层的各个离子的运移方程，可参见上述实施例1求解沉积颗粒的体积浓度的过程来求解得到所述储层内的孔隙的含水体积分数
于此不再进行赘述。
在通过上述方法计算得到所述各个离子i的浓度
之后，再根据上述公式(3-5)可计算得到各个沉淀物的浓度累积量
然后根据所述各个沉淀物的浓度累积量
摩尔质量及密度，可确定所述各个沉淀物的体积浓度
最后可确定所有沉淀物的体积浓度
其中，Ni为离子i与其发生沉淀反应的位于储层中的离子的个数；以及M为所述外来流体中的多个离子的个数。通过上述无机沉淀损害储层的建模方法建立的时空演化模拟方程综合考虑了储层内发生无机沉淀时多种物理化学因素对储层损害的影响，由此通过该步骤S3201求解得到的无机沉淀的体积浓度非常精确。
步骤S3202，基于所确定的所述无机沉淀的体积浓度，确定表征待诊断井的预设区域内的储层的损害程度的特征参数。
其中，所述特征参数可为所述储层的渗透率和/或所述储层的滤失系数。
在一实施例中，所述特征参数可为所述储层的渗透率。
对于步骤S3202，所述确定表征待诊断井的预设区域内的储层的损害程度的特征参数包括：基于所述无机沉淀的体积浓度
及公式(3-6)，确定所述储层的渗透率
在一实施例中，所述特征参数可为所述储层的渗透率。
对于步骤S202，所述确定表征待诊断井的预设区域内的储层的损害程度的特征参数包括：基于所述无机沉淀的体积浓度
及公式(3-7)，确定所述储层的滤失系数
其中，φ0为所述储层的孔隙度的初始值；Cdmax为所述无机沉淀的最大体积浓度；mk与mK分别为第一经验值与第二经验值；
为所述储层的渗透率的初始值；以及
所述储层的滤失系数的初始值。
其中，所述特征参数为所述储层的表皮系数。
对于步骤S202，所述确定表征待诊断井的预设区域内的储层的损害程度的特征参数可包括：基于所述无机沉淀的体积浓度
及公式
确定所述储层的渗透率
以及基于所述储层的渗透率
及公式(3-8)，确定所述储层的表皮系数
其中，
为所述储层的渗透率的初始值，
rw为所述待诊断井的井筒半径，以及rsw为所述储层的损害半径。
通过该步骤S3202得到的特征参数(例如所述储层的渗透率
与表皮系数
)是时空演化4D定量模拟的结果(未示出)。因此，可根据渗透率或表皮系数的演化特点进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义。
综上所述，通过所确定的时空演化模拟方程可确定无机沉淀的体积浓度，再基于所确定的所述无机沉淀的体积浓度，确定表征待诊断井的预设区域内的储层的损害程度的特征参数(例如所述储层的渗透率和/或表皮系数)，由此可定量模拟由无机沉淀引起的储层损害特征的四维时空演化过程，从而对未发生储层损害的井进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义，以及对已损害井优化设计解堵措施、提高或恢复油井产量和水井注水能力，以及提高油藏数值模拟精度都具有十分重大意义。
实施例4-储层内部微粒
储层内部微粒(即粒径小于预设尺寸(例如37微米)的固相颗粒)产生堵塞的本质是所述储层内部微粒的运移和沉积。因此，本发明各实施例的核心是要建立储层内部微粒的运移和沉积的动力学模型。具体地，基于质量守恒、扩散关系等建立储层内部微粒在待诊断井周围的储层中的浓度分布的时空演化控制唯象模型(该模型包含运移微粒的浓度C和沉积微粒的浓度 Cd)，再结合沉积浓度和渗透率等储层损害特征参数间的关系，就可诊断渗透率等储层损害特征参数的时空场分布。
需要说明的是，为了简单描述起见，在本发明的各个实施例中的随时空演化的物理量、化学量可省略变量
例如
可简写为φw；以及
可简写为K。
图4-1是本发明一实施例提供的储层内部微粒损害储层的建模方法的流程图。所述建模方法可包括步骤S4101-S4104。
步骤S4101，确定储层中的流体的速度。
其中，所述储层位于待诊断井(例如注水井、采油井)的预设区域内。
对于步骤S4101，所述确定储层中的流体的速度可包括：建立所述流体进入所述储层的压力传导方程；以及根据所述压力传导方程及达西公式，确定所述流体的速度。
具体地，压力是驱动固-液混合液(即含有所述运移微粒的流体)从注水井的井筒持续侵入周围的储层的动力，由此可建立如公式(4-1)的所述流体进入储层的压力传导方程：
再根据公式(4-1)及达西公式(如下式(4-2))可确定所述流体的速度，
其中，
为所述流体的压力；φ为所述储层的孔隙度；μ为流体粘度；ct为流体-岩石综合压缩系数；
为所述储层的渗透率；以及τ为所述储层的迂曲度。
步骤S4102，基于所述流体的对流参数与扩散参数及所述流体内的运移微粒的质量变化率，建立所述流体与所述储层中的岩石上的沉积微粒之间的质量平衡方程。
其中，所述运移微粒的质量变化率与所述流体的速度之间具有关联关系。下文会详细描述获取所述运移微粒的质量变化率的过程。
在建立微粒运移损害模型时，首先考虑沉积微粒起动运移时流体的临界速度，然后考虑运移微粒如何改变固液流动沉积方程。
根据微粒起动模型，半径为rs的弹性固体微粒在岩石孔隙粗糙内表面的受力情况，如图4-2所示。沉积微粒在流体的冲刷及与岩石表面的相互作用下，受到力和力矩的作用，临界速度即为这些力和力矩刚好平衡时所对应的流体速度。所述临界速度可通过以下方式获取：根据所述沉积微粒的受力情况，建立所述沉积微粒的力矩平衡方程，其中所述沉积微粒的受力情况与所述流体的速度相关；以及根据所述沉积微粒的力矩平衡方程，确定所述临界速度。
具体地，所述沉积微粒受到与流速v同向的拖拽力Fd、重力Fg、静电力Fe、举升力Fl。由于拖拽力Fd与举升力Fl的大小都是流速v的函数，因此，在它们力矩平衡时对应的流速即为临界速度vcr。根据所述沉积微粒的受力情况，建立下式(4-3)所示的所述沉积微粒的力矩平衡方程：
Fd·ld＝(Fe-Fl+Fg)·ln， (4-3)
岩石孔隙粗糙内表面附近的流场对附着在该表面上的沉积微粒产生的拖拽力Fd可通过Navier-Stokes方程的渐近解获得，其表示如下：
Fd＝ωπμrsvcr， (4-4)
其中，μ为所述流体的粘度，rs为所述沉积微粒的半径，vcr为距离表面rs处的流速，ω为拖拽系数(例如，ω＝6×1.7，特别的，如果ω取6，这时的拖拽力对应自由无边界流动中的固相微粒受到的拖拽力)。
剪切流场对于该沉积微粒的举升力Fl可表示如下：
Fl＝χ[ρμ(rsvcr)3]1/2， (4-5)
其中，ρ、μ分别为流体的密度与粘度，χ为举升系数。
所述沉积微粒的重力Fg可表示如下：
其中，g为重力加速度，ρs为所述沉积微粒的密度，以及ρ为流体的密度。
通常，所述沉积微粒的静电力Fe的大小是由静电势对空间的导数求得：
其中，V(h)是总的静电势能，其包含三个部分：VLVA(伦敦-范德瓦尔斯势)、 VDLR(双电层势)及VBR(波恩势)。即V(h)可表示为：
V＝VLVA+VDLR+VBR
VLVA和VDLR可以通过著名的DLVO理论得出：
VBR可表示为：
其中，
H为Hamaker常数，h是沉积微粒表面与介质(例如岩石)表面分开的距离，ε0为沉积微粒的介电常数，De为沉积微粒的双电层常数，ψ01、ψ02分别为沉积微粒和岩石骨架微粒的表面势能，σLJ为原子分子作用 Lennard-Jones势常数，κ为倒德拜长度(即其量纲为长度分之一)。
vcr和所述流体的流速v(真实达西速度，微粒运移模型中的v)有如下关系：
式中，rc为储层孔喉平均半径(即孔隙的平均半径)。联立式(4-4)、(4-8) 即可确定所述流体的流速v和其对表面微粒拖拽力的关系。
(法向力Fn的)力臂ln约为在法向力Fn(Fn＝Fe-Fl+Fg)作用下沉积微粒和基质(例如岩石)接触变形面的半径：
rs为起动微粒(沉积微粒中的部分微粒在流体的冲涮作用下开始运移并合并加入流体中成为运移微粒，其中所述部分微粒成为起动微粒)的半径，K为复合杨氏模量，
其中，E1、E2分别为起动微粒和基质的杨氏模量，ν1、ν2分别为起动微粒和基质的泊松比。
一旦求得ln，ld可由简单的几何关系确定：
通过联立上述公式(4-3)-(4-10)可得到临界速度vcr，其表示如下：
根据公式(4-11)可知，临界速度vcr与沉积微粒和介质(例如岩石)的力学、物理、化学性质有关。只有当储层中的流体的实际速度超过临界速度 vcr，在该流体的作用下沉积微粒才会运移而成为运移微粒(或运移物质源)。一般情况下，越接近待诊断井的井筒中心的流体流速越大，故微粒运移区域应该为井筒附近的一个圆环带。
根据质量方程，设运移微粒的质量变化率(即释放微粒量)为qs，则qs具有如下性质：
也就是说，只有当流体的速度(也可称为流体流速)超过临界速度时，微粒才会起动并加入到流体中参与运移，从而增加流-固混合液的质量。因此，对于步骤S4102，所述建立所述流体与所述储层中的岩石上的沉积微粒之间的质量平衡方程可包括：基于所述流体的对流参数与扩散参数，建立下式表示的所述质量平衡方程，
其中，ρ为所述流体的密度；φ为所述储层的孔隙度；
为所述沉积微粒 的质量分数(也可称为质量浓度)；u为达西表观速度；
为扩散流量，
其中ρL为所述流体的密度，
为运移微粒的扩散系 数，
α为垂直扩散率，
为所述流体的速度；
为单位时间内的所述沉积微粒的累积 质量；t为时间；以及qs为所述运移微粒的质量变化率。
其中，所述运移微粒的质量变化率qs通过以下方式获取：确定所述沉积微粒的释放场的强度Q(r)；确定所述释放场的强度的衰减函数Y(t)；以及根据所述释放场的强度Q(r)与所述释放场的强度的衰减函数Y(t)，确定所述运移微粒的质量变化率qs＝Q(r)Y(t)。具体地，所述释放场的强度Q(r)可为一个常数(q0)，衰减函数Y(t)可为随时间变化的指数衰减函数(例如e-λt，其中λ为衰减常数)。
步骤S4103，基于所述流体的对流参数与扩散参数，建立所述沉积微粒的体积浓度与所述流体的体积浓度之间的连接条件方程。
对于步骤S4103，所述建立所述沉积微粒的体积浓度与所述流体的体积浓度之间的连接条件方程可包括：基于所述流体的对流参数与扩散参数，建立下式(4-13)表示的所述连接条件方程，
其中，ρp为沉积微粒的密度；
为所述沉积微粒的体积浓度；
其中k0为原始滤失系数，
以及
由于F1(T)(F1(T)为与温度相关的指数函数) 和温度的相关性是由exp(1/T)来度量的，在常见的温度范围(比如300K～400K) 内，这个函数的变化其实是很缓慢的，实际上接近等温过程，故
其中
为所述沉积微粒的体积浓度， Cdmax为所述沉积微粒的最大体积浓度，以及mk为第一经验值。上述所有参数都既可以是常数，也可以是随空间变化的参数，即非均质的情况。
步骤S4104，根据所述运移微粒的质量分数与所述运移微粒的体积浓度之间的关系、所述流体的速度、所述质量平衡方程及所述连接条件方程，确定所述储层内部微粒损害储层的时空演化模拟方程。
其中，所述时空演化模拟方程用于模拟由微粒引起的储层损害特征的四维时空演化过程。
其中，所述运移微粒的质量分数与所述运移微粒的体积浓度之间的关系可为
其中ρp为所述沉积微粒的密度；ρL为所述流体的密度；
为所述运移微粒的质量分数；以及
为所述运移微粒的体积浓度。所述微粒损害储层的时空演化模拟方程可包括：公式(4-14)所示的微粒运移损害储层的时空演化模拟方程，以及公式(4-15)所示的微粒沉积损害储层的时空演化模拟方程。
对于步骤S4104，所述确定所述微粒损害储层的时空演化模拟方程可包括：根据所述运移微粒的质量分数与所述运移微粒的体积浓度之间的关系、所述流体的速度及公式(4-12)所示的质量平衡方程，确定由下式(4-14) 所示的微粒运移损害储层的时空演化模拟方程：
以及根据所述运移微粒的质量分数与所述运移微粒的体积浓度之间的关系、所述流体的速度及公式(4-13)所示的连接条件方程，确定由公式(4-15) 所示的微粒沉积损害储层的时空演化模拟方程：
其中，
为所述运移微粒的体积浓度；
为所述流体的速度；τ为所述储层的迂曲度；ρp为所述沉积微粒的密度；ρL为所述流体的密度；
所述储层的滤失系数的初始值；
为所述沉积微粒的体积浓度；Cdmax为所述沉积微粒的最大体积浓度；mk为第一经验值；α为垂直扩散率；φ为所述储层的孔隙度；以及qs为所述运移微粒的质量变化率。其中
其中NR、VPe、NA、NDL、NE1、NE2、NG、NLo、NvdW、ζp(g)分别位半径数、皮克雷数、吸引数、双电层数、第一电位力数、第二电位力数、重力数、伦敦力数、范德瓦尔斯力数及运移微粒和基质颗粒(即沉积在所述岩石上的颗粒)的电势(各个参数的相关表达式详见表2)；
综上所述，本发明创造性地根据储层内的流体的对流参数与扩散参数及运移微粒的质量变化率，建立所述流体与储层中的岩石上的沉积微粒之间的质量平衡方程；根据所述流体的对流参数与扩散参数，建立所述沉积微粒的体积浓度与所述流体的体积浓度之间的连接条件方程；根据所述运移微粒的质量分数与所述运移微粒的体积浓度之间的关系、所述流体的速度、所述质量平衡方程及所述连接条件方程，确定所述储层内部微粒损害储层的时空演化模拟方程。由此，通过所确定的时空演化模拟方程可定量模拟由储层内部微粒引起的储层损害特征的四维时空演化过程，从而对未发生储层损害的井进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义，以及对已损害井优化设计解堵措施、提高或恢复油井产量和水井注水能力，以及提高油藏数值模拟精度都具有十分重大意义。
图4-3是本发明一实施例提供的确定储层损害程度的方法的流程图。如图4-3所示，所述确定储层损害程度的方法可包括步骤S4301-S4302。
步骤S4301，基于所述的储层内部微粒损害储层的建模方法建立的时空演化模拟方程，确定所述沉积微粒的体积浓度。
对于上述公式(4-14)所示的微粒运移损害储层的时空演化模拟方程，可参见上述实施例1求解沉积颗粒的体积浓度的过程来确定所述沉积微粒的体积浓度
于此不再进行赘述。
在通过上述方法计算得到所述运移微粒的体积浓度
之后，再根据上述公式(4-15)可计算得到所述沉积微粒的体积浓度
由此通过上述微粒损害储层的建模方法建立的时空演化模拟方程综合考虑了储层内部微粒运移时多种物理化学因素对储层损害的影响，由此通过该步骤S4301求解得到的沉积微粒的体积浓度非常精确。
步骤S4302，基于所述沉积微粒的体积浓度，确定表征待诊断井的预设区域内的储层的损害程度的特征参数。
其中，所述特征参数可为所述储层的渗透率和/或所述储层的滤失系数。
在一实施例中，所述特征参数可为所述储层的渗透率。
对于步骤S4302，所述确定表征待诊断井的预设区域内的储层的损害程度的特征参数可包括：基于所述沉积微粒的体积浓度
及公式(1-15)，确定所述储层的渗透率
在一实施例中，所述特征参数可为所述储层的渗透率。
对于步骤S302，所述确定表征待诊断井的预设区域内的储层的损害程度的特征参数可包括：基于所述沉积微粒的体积浓度
及公式(1-16)，确定所述储层的滤失系数
其中，所述特征参数可为所述储层的表皮系数。
对于步骤S302，所述确定表征待诊断井的预设区域内的储层的损害程度的特征参数可包括：基于所述沉积微粒的体积浓度
及公式
确定所述储层的渗透率
以及基于所述储层的渗透率
及公式(1-17)，确定所述储层的表皮系数
通过该步骤S4302得到的特征参数(例如所述储层的渗透率
与表皮系数
)是时空演化4D定量模拟的结果(如图4-4所示)。更具体地，图4-5示出了由储层渗透率损害率(基于所述储层的渗透率
及公式
确定所述储层的渗透率损害率I(ri,t)，其中
为
的最大值)表征的在第40天储层内部微粒运移损害储层的半径的示意图(如箭头所指示的半径)，相关工作人员可通过该图4-5直观地确认储层被损害的程度。因此，可根据渗透率或表皮系数的演化特点进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义。
综上所述，通过所确定的时空演化模拟方程可确定所述沉积微粒的体积浓度，再根据所述沉积微粒的体积浓度可确定表征待诊断井的预设区域内的储层的损害程度的特征参数(例如所述储层的渗透率和/或表皮系数)，由此可定量模拟由储层内部微粒引起的储层损害特征的四维时空演化过程，从而对未发生储层损害的井进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义，以及对已损害井优化设计解堵措施、提高或恢复油井产量和水井注水能力，以及提高油藏数值模拟精度都具有十分重大意义。
实施例5-水锁效应
水锁效应主要发生在水湿性储层岩石中，由于水为润湿相，水总是率先占据小孔隙直至中大孔隙，由此会将油气分割为分散相，从而显著降低油气在储层中的渗透情况(例如降低储层的渗透率)。水锁效应受到储层岩性、物性、孔隙结构、侵入流体等多种因素的控制，尤其是储层介质的几何特征对于水锁效应损害储层的影响很大。不同的孔喉结构分布模式和复杂性会导致水润湿相在岩石中的展布模式显著变化，从而影响储层的渗透率。
因此，本发明各实施例的核心是要建立水分子向岩石内部扩散和储层中的孔隙内的水含量变化的动力学模型(即水分子通过固-液界面由孔隙中的液相向固相内部扩散的扩散方程与孔隙中的流体的对流扩散方程)。具体地，基于菲克扩散定律以及储层中的孔隙内的流体的对流扩散关系等建立水锁效应影响待诊断井周围的储层中的孔隙度分布的时空演化控制唯象模型(该模型包含储层内的孔隙的含水体积分数c1与储层中的岩石的含水体积分数的初始值c0)，再结合储层的孔隙度和渗透率等储层损害特征参数间的关系，就可诊断渗透率等储层损害特征参数的时空场分布。
图5-1是本发明一实施例提供的水锁效应损害储层的建模方法的流程图。所述建模方法可包括步骤S5101-S5104。
步骤S5101，确定待诊断井的预设区域内的储层中的流体的达西表观速度。
其中，所述待诊断井可例如为注水井。
对于步骤S5101，所述确定储层中的流体的速度可包括：建立所述流体进入所述储层的压力传导方程；以及根据所述压力传导方程及达西公式，确定所述流体的达西表观速度。
具体确定过程可参见上述实施例2中的确定达西表观速度的过程(即上式(2-1)与(2-2)及其相关描述)。
步骤S5102，根据所述流体的达西表观速度与所述流体中的水分子的扩散系数，建立所述储层的水相运动方程。
在储层条件下，所述储层中的孔隙内的不同位置的含水量满足质量守恒方程。其中，外来水相在所述储层内的运动主要由对流与扩散两个过程决定。具体地，对于步骤S5102，所述建立所述储层的水相运动方程可包括：根据所述流体的达西表观速度u与所述水分子的扩散系数Dw，建立下式表示的所述质量平衡方程：
其中，φ0为所述储层的孔隙度的初始值；
为所述储层内的孔隙被水相占据的绝对孔隙度；以及
为所述储层内的任意点的空间位置(例如，以所述待诊断井的中心为原点)。
根据所述质量平衡方程及所述储层的水相饱和度的时空分布函数
建立下式(5-1)表示的所述水相运动方程：
所述水相运动方程的初始条件为
以及所述水相运动方程的边界条件为
(也就是说，在注水井的井壁处的储层孔隙完全被水充满，即孔隙中的水相饱和度为1)。其中，φ0为所述储层的孔隙度的初始值；rw为所述待诊断井的井筒半径；以及Swc为所述储层中的束缚水饱和度。
步骤S5103，根据所述储层的孔隙的孔径分布特征及所述储层的预设渗透率模型，建立所述储层的渗透率的分布方程。
对于步骤S5103，如图5-2所示，所述建立所述储层的渗透率的分布方程可包括步骤S5201-S5202。
步骤S5201，根据所述储层的孔隙的孔径分布特征，确定孔径为λ的孔隙的体积密度函数与所述储层的水相饱和度的孔径分布方程。
为了定量描述储层的孔隙结构，采用分形理论来研究孔隙结构中的水锁效应。根据分形的几何原理，若储层的孔径分布具有分形特征，则储层中孔径大于λ的孔隙的数量N(>λ)与λ有如下幂函数关系：
其中，D为孔隙的分形维数(0<D<3)。
在所述储层的孔隙的孔径分布特征为所述储层中的孔径大于λ的孔隙的数量N(＞λ)满足上述公式(5-2)的情况下，如图5-3所示，所述步骤S5201 中的确定孔径为λ的孔隙的体积密度函数可包括步骤S5301-S5302。
步骤S5301，根据所述储层中的孔径大于λ的孔隙的数量N(＞λ)，确定所述储层中的孔隙的总体积为
具体地，根据上式公式(5-2)进而可得到储层中的孔隙的总数目N(＞λmin)：
根据公式(5-2)与(5-3)可得到下式(5-4)
则孔径大于λ的孔隙的数量N(＞λ)与λ之间的关系为下式(5-2)所表示的幂函数关系：
公式(5-2)中λ、λmin与λmax分别为所述孔隙的孔径、最小孔径与最大孔径(λmin与λmax可以由平均孔径与孔径分布标准差得出，一般而言有
)；a 为比例常数。
接着，由公式(5-5)可得到储层的孔径分布密度函数f(λ)满足下式(5-6)：
式中a'＝-Da为比例常数。
由上式(5-6)表示的所述孔径分布密度函数可以得到储层中的孔隙的总体积分形表达式：
其中α为与孔隙的形状有关的常数(孔隙的形状为立方体时α＝1，或者孔隙的形状为球体时α＝π/6)，通过积分可得：
类似地，所述储层中的孔径小于λ的孔隙的体积
步骤S5302，根据所述储层中的孔隙的总体积Φmax及所述储层中的孔径小于λ的孔隙的体积
确定所述孔径为λ的孔隙的体积密度函数为
其中，D为所述孔隙的分形维数；以及λ、λmin与λmax分别为所述孔隙的孔径、最小孔径与最大孔径；以及A＝αa'/(3-D)(常数)。
步骤S5202，根据所述预设渗透率模型、所述孔径为λ的孔隙的体积密度函数及所述水相饱和度的孔径分布方程，建立所述储层的渗透率的分布方程。
如图5-4所示，步骤S5202中的确定所述储层的水相饱和度的孔径分布方程可包括步骤S5401-S5402。
步骤S5401，根据所述储层中的孔径大于λ的孔隙的数量N(＞λ)，确定被非水相占据的孔隙的体积为
假设孔径小于λ的孔隙完全被水相占据，而孔径大于λ的孔隙完全被非水相占据(亦即储层中的岩石为水湿(即亲水性))，则结合上述公式(5-8) 可得到非水相占据的孔隙体积Φnw(λ)，
步骤S5402，根据所述储层中的孔隙的总体积Φmax与所述被非水相占据的孔隙的体积Φnw(λ)，确定下式表示的所述水相饱和度的孔径分布方程
其中，D为所述孔隙的分形维数；λ、λmin与λmax分别为所述孔隙的孔径、最小孔径与最大孔径；以及A＝αa'/(3-D)。
具体地，根据公式(5-8)与公式(5-10)可确定非水相饱和度的孔径分布方程为
再根据所述非水相饱和度的孔径分布方程Snw(λ)可确定所述水相饱和度的孔径分布方程，
由线性Hagen-Poiseuille粘性流可知，毛管束模型的渗透率可表示为
在本发明实施例中，可将储层的渗透通道看做是多个毛管束累积起来的。由于孔隙的尺寸分布具有连续性，因此可以将所述毛管束模型的渗透率的表达式用积分形式写出：
在所述储层的预设渗透率模型满足
的情况下，所述建立所述储层的渗透率的分布方程可包括：根据所述储层的预设渗透率模型
所述孔径为λ的孔隙的体积密度函数dξ及所述水相饱和度的孔径分布方程，建立渗透率的分布方程(未示出)。
根据所建立的渗透率的分布方程，可进一步建立下式表示的所述渗透率损害率建立下式表示的所述渗透率损害率的分布方程，
具体地，首先，将公式(5-9)代入公式(5-12)可得出所述储层的渗透率的孔径分布函数：
，然后，根据式(5-11)中的λ与含水饱和度Sw(λ)之间的关系式将上述渗透率的孔径分布函数K(λ)中的变量λ替换为Sw，可得到公式(5-13)。
步骤S5104，根据所述渗透率的分布方程及所述水相运动方程，确定水锁效应损害储层的时空演化模拟方程，其中所述时空演化模拟方程用于模拟由水锁效应引起的储层损害特征的四维时空演化过程。
具体地，根据公式(5-1)可求解到所述储层的水相饱和度的时空分布函数
并将
代入所述储层的渗透率在四维时空上的分布形式，即水锁效应损害储层的时空演化模拟方程。
综上所述，本发明创造性地确定待诊断井的预设区域内的储层中的流体的达西表观速度；根据所述流体的达西表观速度与所述流体中的水分子的扩散系数，建立所述储层的水相运动方程；建立所述储层的渗透率的分布方程；以及根据所述渗透率的分布方程及所述水相运动方程，确定水锁效应损害储层的时空演化模拟方程。由此，通过所确定的时空演化模拟方程可定量模拟由水锁效应引起的储层损害特征的四维时空演化过程，从而对未发生储层损害的井进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义，以及对已损害井优化设计解堵措施、提高或恢复油井产量和水井注水能力，以及提高油藏数值模拟精度都具有十分重大意义。
相应地，本发明另一实施例还提供一种确定储层损害程度的方法，所述方法可包括：基于所述的水锁效应损害储层的建模方法建立的时空演化模拟方程，确定表征待诊断井的预设区域内的储层的损害程度的特征参数。
对于上述所示的水锁效应损害储层的时空演化模拟方程的求取，需要根据公式(5-1)计算得到
具体求解过程可参见上述实施例1求解沉积颗粒的体积浓度的过程，于此不再进行赘述。
在通过上述方法计算得到所述储层的水相饱和度
之后，再根据上述公式(5-13)可计算得到所述储层的渗透率
(当然，也可得到储层中的特定位置r处的渗透率
随时间变化情况，如图5-5所示)，由此通过上述水锁效应损害储层的建模方法建立的时空演化模拟方程综合考虑了水锁损害时多种物理化学因素对储层损害的影响，由此通过本实施例求解得到的储层的渗透率非常精确。
在得到所述储层的渗透率的基础上，可计算表征待诊断井的预设区域内的储层的损害程度的特征参数。
在一实施例中，所述特征参数可为所述储层的渗透率损害率。
相应地，所述确定表征待诊断井的预设区域内的储层的损害程度的特征参数可包括：基于所述时空演化模拟方程，确定所述储层的渗透率
以及基于所述储层的渗透率
及公式(5-14)，确定所述储层的渗透率损害率
其中，
为
的最大值。
在另一实施例中，所述特征参数可为所述储层的表皮系数。
所述确定表征待诊断井的预设区域内的储层的损害程度的特征参数可包括：基于所述时空演化模拟方程，确定所述储层的渗透率
以及基于所述储层的渗透率
及公式(5-15)，确定所述储层的表皮系数
其中，
为所述储层的渗透率的初始值；以及
rw为所述待诊断井的井筒半径，以及rsw为所述储层的损害半径。
通过上述各个实施例得到的特征参数(例如所述储层的渗透率
与表皮系数
)是时空演化4D定量模拟的结果(图5-6示出了位置
处的表皮系数随时间变化的情形)。因此，可根据渗透率或表皮系数的演化特点进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义。
实施例6-应力敏感性
储层内部的孔隙结构和固体基质(即岩石基质或骨架颗粒)的分布及其复杂。微观的一个固体基质(即团簇)形成一条渗流通道，多个固体基质(即团簇)之间夹杂着多个孔隙(相当于水力管束)，而所述多个固体基质与所述多个孔隙形成如图6-1A所示的毛管束(也称为毛细管束)。这些孔隙具有特定的迂曲度，使固体基质与孔隙的实际长度Lc0比由毛管束形成的多孔介质的长度L0(即毛管束的表观长度)更长，如图6-1A所示。
并且，一方面因储层的上覆岩层的作用而使储层受到向下的压力，另一方面因储层内部的流体的支撑作用而使储层受到向上的压力，由此，储层在这两个压力(即有效应力)作用下，毛管束收缩(如图6-1B所示)、毛管束延展(如图6-1C所示)，进而会影响液体在储层中的渗透率(即简称为所述储层的渗透率)。应力敏感性受储层的孔隙结构、岩性、物性、侵入流体等多种因素的影响。因此，本发明各实施例的核心是要建立有效应力作用下的储层中的流体的流量与渗透率的分布方程。具体地，根据储层的孔隙的孔径分布特征、有效应力下的储层的各个毛管束的直径与长度及流体流量公式，确定所述储层中的流体的流量，再结合渗透率模型，就可诊断渗透率等储层损害特征参数的时空场分布。
图6-2是本发明一实施例提供的应力敏感性损害储层的建模方法的流程图。如图6-2所示，所述建模方法包括步骤S6201-S6203。
步骤S6201，确定待诊断井的预设区域内的储层受到的有效应力。
对于步骤S6201，所述确定待诊断井的预设区域内的储层受到的有效应力可包括：建立所述流体进入所述储层的压力传导方程；以及根据所述储层的压力传导方程与所述储层的上覆岩层的压力，确定所述储层受到的有效应力。
具体地，压力是驱动固-液混合液从注水井的井筒持续侵入待诊断井周围的储层的动力，由此可建立如公式(6-1)的所述流体进入储层的压力传导方程：
其中，所述压力传导方程的初始条件为
以及边界条件为
(也就是说，在注水井的井壁处的压力为Pw)。
然后根据公式(6-1)可得到所述流体对所述储层产生的压力
(即孔隙压力)，再根据
与上覆岩层的压力
确定所述储层受到的有效应力
步骤S6202，根据所述储层的孔隙的孔径分布特征、所述有效应力下的所述储层的各个毛管束的直径与长度、及流体流量公式，确定所述储层中的流体的流量。
其中，所述毛管束由多个固体基质及所述多个固体基质之间的孔隙组成。
对于步骤S6202，如图6-3所示，所述确定所述储层中的流体的流量可包括步骤S6301-步骤S6302。
在执行步骤S6301之前，所述确定所述储层中的流体的流量还可包括：根据所述储层的弹性模量，确定所述有效应力下的所述储层的各个毛管束的直径；以及根据所述储层的弹性模量与泊松比，确定所述有效应力下的所述储层的各个毛管束的长度。
下面，分别介绍关于所述有效应力下的所述储层的各个毛管束的直径与长度的确定过程。
待诊断井的横波速度是确定储层的泊松比、弹性模量等重要数据的关键特征。而在现场测井时，一般不会针对待诊断井进行横波测井，而是通过含泥量和关于含泥量的线性公式确定泊松比。其中，所述线性公式通过特定岩石样本拟合而成，故针对不同的储层，通过所述线性公式确定的泊松比的精确性非常差。在本实施例中，根据已获取的待诊断井周围的多个邻井的横波速度、纵波速度及对待诊断井进行现场测井所获取的纵波速度，可间接计算得到所述诊断井的横波速度。
在执行所述确定所述有效应力下的所述储层的各个毛管束的长度的步骤之前，所述确定所述储层中的流体的流量还可包括：根据位于所述预设区域内的多个特定相邻井中的每一者的横波速度与纵波速度，确定所述储层的横波速度；以及根据所述储层的横波速度，确定所述储层的所述弹性模量与所述泊松比。
其中，所述特定相邻井可为已通过测井方式获取横、纵波速度的井。
具体地，所述待诊断井的横波速度υt可由纵波速度υl及关系式υt＝aυl-b 确定。其中，所述关系式中的系数a、b可基于统计的邻井的横波速度与纵波速度计算得到。例如，通过最小二乘法来线性回归出所述系数a、b：
其中，n为通过测井得到其横、波速度的所述预设区域内的邻井的数目；υti为第i个邻井的横波速度；υli为第i个邻井的纵波速度；
为n个邻井的横波速度的平均值；
为n个邻井的纵波速度的平均值。
与上述现有的方式获取的泊松比相比，由于测井数据非常精确，故本实施例通过邻井的横、纵波速度等测井数据获取的待诊断井的横波速度的结果非常精确，从而后续获取的泊松比等数据也非常准确。接着下面将介绍如何通过待诊断井的横波速度获取储层的泊松比与弹性模量。
具体地，可通过下式(6-3)与所述待诊断井的横波速度，得到储层的泊松比ν与弹性模量E：
其中，ρ为所述储层的体积密度；Δtt为横波时差(可由单发射双接收声速测量装置测量得到，即垂直方向上相距L的两个接收器接收到声波所形成的滑行波横波的时间差)；以及Δt为声波时差(可由单发射双接收声速测量装置测量得到，即垂直方向上相距L的两个接收器接收到声波所形成的滑行波的时间差)。
下面基于所获取的弹性模量与泊松比，确定一定应力作用下的储层的任一毛管束的直径与长度。
第一，所述确定所述有效应力下的所述储层的各个毛管束的直径可包括：根据所述储层的弹性模量E及下式，确定所述有效应力
下的所述储层的特定毛管束的直径λp，
其中，F为所述特定毛管束的形状因子；以及λc0为无有效应力下的所述储层的特定毛管束的初始直径。
第二，所述确定所述有效应力下的所述储层的各个毛管束的长度可包括： 根据无有效应力下的所述储层的特定毛管束的表观长度、无有效应力下的所 述储层的特定毛管束的初始直径及分形尺度法则，确定无有效应力下的所述 储层的特定毛管束的初始长度Lc0；以及根据所述储层的弹性模量E、所述储 层的泊松比ν、无有效应力下的所述储层的特定毛管束的初始长度Lc0及下式， 确定所述有效应力
下的所述储层的特定毛管束的长度Lp， 
其中，所述确定无有效应力下的所述储层的特定毛管束的初始长度Lc0可包括：根据无有效应力下的所述储层的特定毛管束的表观长度L0、无有效应力下的所述储层的特定毛管束的初始直径λc0及分形尺度法则
确定所述特定毛管束的初始长度Lc0，其中，DcT为所述特定毛管束的迂回分形维度。
步骤S6301，根据所述有效应力下的所述储层的各个毛管束的直径与长度及流体流量比率公式，确定所述流体的流量比率。
根据上述得到的所述有效应力
下的所述储层的特定毛管束的直径λp与长度Lp、及流体流量比率公式
确定所述流体的流量比率
其中，G为所述特定毛管束的横截面的几何因子；Δp为所述储层的横截面两端的驱替压力；以及F为所述特定毛管束的形状因子。
步骤S6302，根据所述储层的孔隙的孔径分布特征、所述储层中的流体的流量比率及所述流体流量公式，确定所述流体的流量。
对于步骤S6302，所述确定所述流体的流量包括：根据所述储层的孔隙的孔径分布特征
所述储层中的流体的流量比率
及所述流体流量公式
确定所述流体的流量
其中，λc0为所述储层的特定毛管束的初始直径；λc0max为所述储层的最大毛管束与最小毛管束的初始直径；Dcf为所述孔隙的分形维度；以及
为所述有效应力。
具体地，基于岩石分形几何理论，直径为λc0的固体基质的孔隙的数密度 Nc(即单位体积内某一直径固体基质的数量累积分布，或孔隙的孔径分布特征)满足关系：
对于二维(储层)模型，0＜Dcf＜2；对于三维(储层)模型，0＜Dcf＜3。
根据上述储层的孔隙的孔径分布特征，可得到
然后将dNc及公式(6-4)代入所述流体流量公式，可得到下式(6-5)表示的所述流体的流量
步骤S6203，根据所述储层的渗透率模型及所述储层中的流体的流量，确定应力敏感性损害储层的时空演化模拟方程。
其中，所述时空演化模拟方程用于模拟由应力敏感性引起的储层损害特征的四维时空演化过程。
对于步骤S6203，所述确定应力敏感性损害储层的时空演化模拟方程可包括：根据所述储层的渗透率模型
及所述储层中的流体的流量
确定应力敏感性损害储层的时空演化模拟方程。
其中，μ为所述流体的粘度；L0为无有效应力下的所述储层的特定毛管束的表观长度；A为所述储层的横截面的面积；以及Δp为所述储层的横截面两端的驱替压力。具体地，
然后代入上述dNc可得到A。
具体地，根据所述储层的渗透率模型
及上式(6-5) 表示的流量
可得到下式(6-6)表示的应力敏感性损害储层的时空演化模拟方程，
综上所述，本发明创造性地确定待诊断井的预设区域内的储层受到的有效应力；根据所述储层的孔隙的孔径分布特征、所述有效应力下的所述储层的各个毛管束的直径与长度、及流体流量公式，确定所述储层中的流体的流量；以及根据所述储层的渗透率模型及所述储层中的流体的流量，确定应力敏感性损害储层的时空演化模拟方程。由此，通过所确定的时空演化模拟方程可定量模拟由应力敏感性引起的储层损害特征的四维时空演化过程，从而对未发生储层损害的井进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义，以及对已损害井优化设计解堵措施、提高或恢复油井产量和水井注水能力，以及提高油藏数值模拟精度都具有十分重大意义。
相应地，本发明另一实施例还提供一种确定储层损害程度的方法，所述方法包括：基于所述的应力敏感性损害储层的建模方法建立的时空演化模拟方程，确定表征待诊断井的预设区域内的储层的损害程度的特征参数。
对于上述公式(6-6)所示的压力敏感性损害储层的时空演化模拟方程的求取，需要根据公式(6-1)计算得到
具体求解过程可参见上述实施例1求解沉积颗粒的体积浓度的过程，于此不再进行赘述。然后再根据计算得到的
确定所述储层的水相饱和度
在通过上述方法计算得到
与所述储层的水相饱和度
之后，再根据上述公式(6-6)可计算得到所述储层的渗透率
由此通过上述压力敏感性损害储层的建模方法建立的时空演化模拟方程综合考虑了压力敏感性损害时多种物理化学因素对储层损害的影响，由此通过本实施例求解得到的储层的渗透率非常精确。
在得到所述储层的渗透率的基础上，可计算表征待诊断井的预设区域内的储层的损害程度的特征参数。
在一实施例中，所述特征参数可为所述储层的滤失系数。
相应地，所述确定表征待诊断井的预设区域内的储层的损害程度的特征参数可包括：基于所述时空演化模拟方程，确定所述储层的渗透率
以及基于所述储层的渗透率
及公式(5-14)，确定所述储层的渗透率损害率
在另一实施例中，所述特征参数可为所述储层的表皮系数。
所述确定表征待诊断井的预设区域内的储层的损害程度的特征参数可包括：基于所述时空演化模拟方程，确定所述储层的渗透率
以及基于所述储层的渗透率
及公式(5-15)，确定所述储层的表皮系数S。
通过上述各个实施例得到的特征参数(例如所述储层的渗透率
表皮系数
)与渗透率损害率是时空演化4D定量模拟的结果(如图6-4 所示)。更具体地，图6-5示出了由储层渗透率损害率表征的在第40天应力敏感性损害储层的半径的示意图(如箭头所指示的半径)，相关工作人员可通过该图5直观地确认储层被损害的程度。因此，可根据渗透率或表皮系数的演化特点进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义。
实施例7-出砂
根据出砂机理可知，出砂导致储层损害的关键参数为门限流速和临界流速，那么出砂损害储层时空演化定量模拟的核心就是求解门限流速和临界流速。门限流速为砂粒起动的流速，当流体流速大于门限流速，部分砂粒则开始运动，此时出的砂是附着砂；若流体流速超过临界流速，则岩石骨架剪切破坏，开始大量出砂，此时出的砂包含附着砂和骨架砂。因为在流体流速超过门限流速且低于临界流速时，储层才开始部分出砂，其对储层(例如渗透率)几乎无影响，且适当出砂反而会对渗透率有帮助；而在流体流速大于临界流速时，储层大量出砂时才会对储层(例如渗透率)造成较大的损害，故在本发明各个实施例中仅考虑流体流速大于临界流速对储层产生影响的情况。临界流速与临界生产压差(或临界产量)相关，因此，本文涉及如何确定临界生产压差(或临界产量)的具体过程。在此基础上，基于质量守恒、扩散关系等建立砂粒在待诊断井周围的储层中的浓度分布的时空演化控制唯象模型(该模型包含砂粒的浓度C和沉积砂粒的浓度Cd)，再结合沉积浓度和渗透率等储层损害特征参数间的关系，就可诊断渗透率等储层损害特征参数的时空场分布。
图7-1是本发明一实施例提供的出砂损害储层的建模方法的流程图。所述建模方法可包括步骤S7101-S7104。
步骤S7101，确定储层中的流体的速度。
其中，所述储层位于待诊断井(例如注水井)的预设区域内。
对于步骤S7101，所述确定储层中的流体的速度可包括：建立所述流体进入所述储层的压力传导方程；以及根据所述压力传导方程及达西公式，确定所述流体的速度。
具体确定过程可参见上述实施例1中的确定流体速度的过程(即上式 (1-1)与(1-2)及其相关描述)。
步骤S7102，基于所述流体的对流参数与扩散参数及所述流体内的砂粒的质量变化率，建立所述流体与所述储层中的岩石上的沉积砂粒之间的质量平衡方程。
其中，所述砂粒的质量变化率与所述储层的原油产量之间具有关联关系。具体地，所述砂粒的质量变化率与所述储层的原油产量之间的关联关系包括：在所述储层的原油产量大于临界产量的情况下，所述砂粒的质量变化率大于0。
在一实施例中，所述临界产量通过以下方式获取：根据所述压力传导方程确定所述流体的压力；根据所述储层的有效径向应力与有效周向应力及摩尔库伦准则，确定所述储层开始出砂时所述流体的临界井底流压；以及根据所述流体的临界井底流压、所述流体的压力及裘比公式，确定所述临界产量。
在建立出砂损害模型时，首先考虑砂粒起动运移时流体的临界速度，然后考虑砂粒如何改变固液流动沉积方程。
地层被钻开后，井周的应力分布会发生变化，且在井壁处产生应力集中。
在井内液柱压力(下式右侧第一项)、地应力(下式右侧第二项与第三) 和流体渗流(下式第四项)的联合作用下，所述储层(例如储层岩石是各项同性的均匀弹性体)的径向应力
与周向应力
分别为：
其中，rw为所述待诊断井的井筒半径；
为所述流体的压力；σH、σh分别为最大水平地应力与最小水平地应力；
为所述流体的流速；φ为所述储层的孔隙度；β为孔隙Biot系数(
其中，Cr为岩石基岩压缩系数；Cb为岩石体积压缩系数)；θ为储层中
处的径向与水平最大地应力方向的夹角；以及当井壁有渗透时δ为1，当井壁无渗透时δ为0。这里考虑井壁有渗透的情况(即δ为1)，
对于多孔渗透性储层，所述储层的径向应力与有效径向应力σ′r之间的关系满足下式(7-3)，以及所述储层的周向应力与有效周向应力σ′θ之间的关系满足下式(7-4)：
忽略构造应力，由Anderson单轴应变模型(σH＝σh)知，油田地应力为：
其中，σv为所述储层的上覆岩层的压力。由密度测井资料可知，上覆岩层压力σv可通过以下等式获取：
若采用上覆岩层梯度理论值，σv＝22.7H，其中，H为深度；若假设上覆岩层压力随深度均匀增加，则σv＝[ρS(1-φ)+ρLφ]gH，其中，ρS为岩石骨架的平均密度；ρL为所述流体的密度；以及H为深度。
根据上述公式(8-3)-(7-5)可确定所述储层的有效径向应力与有效轴向应力分别满足以下等式：
当r＝rw且深度为井底所在的深度时，
由于井壁处的压差最大，故只有井壁处的储层出砂之后，其他位置才会出砂，在本实施例中，考虑井壁表面的有效径向应力与有效周向应力为：
地层岩石在地应力、地层中孔隙压力和流体渗流时的摩擦力等综合作用下，当有效径向应力增大σ′θ超出非固结砂岩储层中屈服条件会造成储层岩石的不稳定性和塑形流动，从而引发出砂，储层岩石破坏遵循摩尔库伦准则。具体地，当用最大主应力σmax与最小主应力σmin表示时，考虑储层孔隙中的流体压力，摩尔库伦准则可表示为：
其中，C为砂泥岩内聚力；地层内摩擦角
νpo为岩石泊松比；vp为纵波速度；vmcl为泥质含量(％)。
对于砂泥岩内聚力C，所述砂泥岩内聚力C与声波传播速度vp满足以下关系：
其中，ρb(H)为深度H处的所述储层中的岩石体积密度。对于泥质含量vmcl，可根据自然电位测井资料并利用经验公式计算得到：
其中，PSP为含泥质砂岩的自然电位；以及SSP为所述预设区域内的含水纯砂岩的静自然电位。对于地层内摩擦角
其可通过以下等式计算得到：
σmax＝σ′θ且σmin＝σ′r，则将等式(7-8)-(7-9)代入上式(7-10)可以得到所述待诊断井开始出砂时的临界井底流压：
由此，可根据上式(7-11)得到临界生产压差ΔPcr为：
在得出临界生产压差ΔPcr之后，再结合裘比公式(7-13)可确定临界产量Qcr为：
其中，Bo为油相体积系数；re为油藏半径；rw为所述待诊断井的井筒半径；μo为原油粘度；以及K为所述储层的渗透率。
根据公式(7-13)可知，临界产量Qcr与流体的速度密切有关。只有当储层中的流体的实际速度超过临界速度(或者原油的实际产量超过临界产量)，在该流体的作用下所述储层才会出砂。
根据质量方程，设砂粒的质量变化率(即释放砂粒量)为qs，则qs具有如下性质：
也就是说，只有当原油的产量超过临界产量时，所述储层才会出砂并且因出砂而形成的砂粒加入到流体中参与运移，从而增加流-固混合液的质量。因此，对于步骤S102，所述建立所述流体与所述储层中的岩石上的沉积砂粒之间的质量平衡方程可包括：基于所述流体的对流参数与扩散参数，建立下式表示的所述质量平衡方程，
其中，ρ为所述流体的密度；φ为所述储层的孔隙度；
为所述沉积砂粒的质量分数(也可称为质量浓度)；u为达西表观速度；
为扩散流量，
其中ρL为所述流体的密度，
为砂粒的扩散系数，
α为垂直扩散率，
为所述流体的速度；
为单位时间内的所述沉积砂粒的累积质量；t为时间；以及qs为所述砂粒的质量变化率。
其中，所述砂粒的质量变化率qs通过以下方式获取：确定所述沉积砂粒的释放场的强度Q(r)；确定所述释放场的强度的衰减函数Y(t)；以及根据所述释放场的强度Q(r)与所述释放场的强度的衰减函数Y(t)，确定所述砂粒的质量变化率qs＝Q(r)Y(t)。具体地，所述释放场的强度Q(r)可为一个常数(q0)，衰减函数Y(t)可为随时间变化的指数衰减函数(例如e-λt，其中λ为衰减常数)。
步骤S7103，基于所述流体的对流参数与扩散参数，建立所述沉积砂粒的体积浓度与所述流体的体积浓度之间的连接条件方程。
对于步骤S7103，所述建立所述沉积砂粒的体积浓度与所述流体的体积浓度之间的连接条件方程可包括：基于所述流体的对流参数与扩散参数，建立下式(7-15)表示的所述连接条件方程，
其中，ρp为沉积砂粒的密度；Cd为所述沉积砂粒的体积浓度；
其中k0为原始滤失系数，
以及
由于F1(T)和温度的相关性是由exp(1/T)来度量的，在常见的温度范围(比如300K～400K)内，这个函数的变化其实是很缓慢的，实际上接近等温过程，故
其中
为所述沉积砂粒的体积浓度，Cdmax为所述沉积砂粒的最大体积浓度，以及mk为第一经验值。上述所有参数都既可以是常数，也可以是随空间变化的参数，即非均质的情况。
步骤S7104，根据所述砂粒的质量分数与所述砂粒的体积浓度之间的关系、所述流体的速度、所述质量平衡方程及所述连接条件方程，确定所述出砂损害储层的时空演化模拟方程。
其中，所述运移砂粒的质量分数与所述运移砂粒的体积浓度之间的关系可为
其中ρp为所述沉积砂粒的密度；ρL为所述流体的密度；
为所述砂粒的质量分数；以及
为所述砂粒的体积浓度。所述出砂损害储层的时空演化模拟方程可包括：公式(7-16)所示的出砂损害储层的时空演化模拟方程，以及公式(7-17)所示的砂粒沉积损害储层的时空演化模拟方程。
对于步骤S7104，所述确定所述出砂损害储层的时空演化模拟方程可包括：根据所述砂粒的质量分数与所述砂粒的体积浓度之间的关系、所述流体的速度及公式(7-14)所示的所述质量平衡方程，确定由下式(7-16)表示的所述出砂损害储层的时空演化模拟方程：
以及根据所述砂粒的质量分数与所述砂粒的体积浓度之间的关系、所述流体的速度及公式(7-15)所示的连接条件方程，确定由公式(7-17)所示的砂粒沉积损害储层的时空演化模拟方程：
其中，
为所述砂粒的体积浓度；
为所述流体的速度；τ为所述储层的迂曲度；ρp为所述沉积砂粒的密度；ρL为所述流体的密度；
所述储层的滤失系数的初始值，
为所述沉积砂粒的体积浓度，Cdmax为所述沉积砂粒的最大体积浓度，mk为第一经验值；α为垂直扩散率；φ为所述储层的孔隙度；以及qs为所述砂粒的质量变化率。
其中
其中 NR、NPe、NA、NDL、NE1、NE2、NG、NLo、NvdW、ζp(g)分别位半径数、皮克雷数、吸引数、双电层数、第一电位力数、第二电位力数、重力数、伦敦力数、范德瓦尔斯力数及砂粒和基质颗粒(即沉积在所述岩石上的颗粒)的电势(各个参数的相关表达式详见表2)。
综上所述，本发明创造性地基于所述流体的对流参数与扩散参数及所述流体内的砂粒的质量变化率，建立所述流体与所述储层中的岩石上的沉积砂粒之间的质量平衡方程，其中所述砂粒的质量变化率与所述储层的原油产量之间具有关联关系；基于所述流体的对流参数与扩散参数，建立所述沉积砂粒的体积浓度与所述流体的体积浓度之间的连接条件方程；以及根据所述砂粒的质量分数与所述砂粒的体积浓度之间的关系、所述流体的速度、所述质量平衡方程及所述连接条件方程，确定所述出砂损害储层的时空演化模拟方程。由此，通过所确定的时空演化模拟方程可定量模拟由出砂引起的储层损害特征的四维时空演化过程，从而对未发生储层损害的井进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义，以及对已损害井优化设计解堵措施、提高或恢复油井产量和水井注水能力，以及提高油藏数值模拟精度都具有十分重大意义。
图7-2是本发明一实施例提供的确定储层损害程度的方法的流程图。如图7-2所示，所述确定储层损害程度的方法可包括步骤S7201-S7202。
步骤S7201，基于所述的出砂损害储层的建模方法建立的时空演化模拟方程，确定所述沉积砂粒的体积浓度。
对于上述公式(7-16)所示的出砂损害储层的时空演化模拟方程，可参见上述实施例1求解沉积颗粒的体积浓度的过程来计算得到所述砂粒的体积浓度
在通过上述方法计算得到所述砂粒的体积浓度
之后，再根据上述公式(7-17)可计算得到所述沉积砂粒的体积浓度
由此通过上述砂粒损害储层的建模方法建立的时空演化模拟方程综合考虑了出砂时多种物理化学因素对储层损害的影响，由此通过该步骤S7201求解得到的沉积砂粒的体积浓度非常精确。
步骤S7202，基于所述沉积砂粒的体积浓度，确定表征待诊断井的预设区域内的储层的损害程度的特征参数。
在一实施例中，所述特征参数可为所述储层的渗透率。
对于步骤S7202，所述确定表征待诊断井的预设区域内的储层的损害程度的特征参数可包括：基于所述沉积砂粒的体积浓度
及(1-15)，确定所述储层的渗透率
在一实施例中，所述特征参数可为所述储层的滤失系数。
对于步骤S7202，所述确定表征待诊断井的预设区域内的储层的损害程度的特征参数可包括：基于所述沉积砂粒的体积浓度
及(1-15)，确定所述储层的渗透率
以及基于所述沉积砂粒的体积浓度
及公式 (1-16)，确定所述储层的滤失系数
其中，所述特征参数可为所述储层的表皮系数。
对于步骤S7202，所述确定表征待诊断井的预设区域内的储层的损害程度的特征参数可包括：基于所述沉积砂粒的体积浓度
及公式
确定所述储层的渗透率
以及基于所述储层的渗透率
及公式(1-17)，确定所述储层的表皮系数
通过该步骤S7202得到的特征参数(例如所述储层的渗透率
与表皮系数
)是时空演化4D定量模拟的结果(如图7-3所示)。更具体地，图7-4示出了由储层渗透率损害率(基于所述储层的渗透率
及公式
确定所述储层的渗透率损害率I(ri,t)，其中
为
的最大值)表征的在第40天出砂损害储层的半径的示意图(如箭头所指示的半径)，相关工作人员可通过该图7-4直观地确认储层被损害的程度。因此，可根据渗透率或表皮系数的演化特点进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义。
综上所述，通过所确定的时空演化模拟方程可确定所述沉积砂粒的体积浓度，再根据所述沉积砂粒的体积浓度可确定表征待诊断井的预设区域内的储层的损害程度的特征参数(例如所述储层的渗透率和/或表皮系数)，由此可定量模拟由出砂引起的储层损害特征的四维时空演化过程，从而对未发生储层损害的井进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义，以及对已损害井优化设计解堵措施、提高或恢复油井产量和水井注水能力，以及提高油藏数值模拟精度都具有十分重大意义。
实施例8-润湿反转
润湿反转是指储层的孔隙表面由亲水性变为亲油性的现象，由此油相在孔隙中的流通性变弱，从而导致储层的渗透性变差。油相饱和度高时，油相占据大孔隙并且相连通性较好，其流动具有接近毛细管流的特征，渗透率和饱和度的二次方呈线性关系；油相饱和度低时，油相主要呈分散相附着在较小孔隙壁面，相连通性相对较差，随着饱和度的降低，油相渗透率的下降速率更快(例如呈4次方律变化)。
一方面，储层润湿性损害程度由油相相对渗透率和含油(或含水)饱和度之间的关系决定；另一方面，储层的压力分布也影响着孔隙内部的流体流速与渗透率。因此，本发明各实施例的核心是要建立油相的压力分布场(分别考虑流体中的油相与水相的压力分布情况)与油相的对流扩散定律。具体地，根据待诊断井的预设区域内的储层的压力分布方程与毛细管的受力平衡条件，确定所述油相的压力分布场；再结合达西公式与油相的对流扩散定律，就可诊断渗透率等储层损害特征参数的时空场分布。
图8-1是本发明一实施例提供的润湿反转损害储层的建模方法的流程图。如图8-1所示，所述建模方法包括步骤S8101-S8104。
步骤S8101，根据待诊断井的预设区域内的储层的压力分布方程，确定所述储层中的水相的压力分布场与毛细管的压力分布场之间的关系。
其中，所述毛细管由所述水相与所述储层中的油相的接触界面的润湿反转而形成。
通常，采用液体(油相与水相的混合物)的压力场来整体描述储层的压力分布场，而在本实施例中将油相与水相进行分离，分别通过油相的压力分布与水相的压力分布来考虑储层的压力分布情况，由此，可更贴近储层中的实际情况来模拟储层的压力分布，从而通过其中油相的压力分布场可精确地模拟润湿反转损害储层的时空演化模拟方程(即获得非常精确的储层渗透率结果)。
对于步骤S8101，根据下式(8-1)表示的待诊断井的预设区域内的储层的压力分布方程，确定所述储层中的水相的压力分布场
与毛细管的压力分布场
之间的关系，
其中，φ为储层的孔隙度(常数)；ct为储层的综合压缩系数(常数)；k、 krw、kro为所述储层的渗透率、储层油相相对渗透率、水相相对渗透率；μw、μo为储层水相粘度、油相粘度。
步骤S8102，根据所述水相的压力分布场与所述毛细管的压力分布场之间的关系及所述毛细管的受力平衡条件，确定所述油相的压力分布场。
所述毛细管的受力平衡条件可为下式(8-2)表示的三力平衡条件，
其中，
为所述毛细管的压力且
由所述毛细管内的有效水饱和度决定；
为所述油相的压力分布场；以及
为所述水相的压力分布场。
其中，所述
由所述毛细管内的有效水饱和度决定可包括：根据所述有效水饱和度及下式(8-3)，确定
其中，
为所述毛细管内的有效水饱和度；m为科里(corey)常数；以及Pce为所述毛细管的压力阈值。
具体地，所述有效水饱和度
可通过以下方式确定：根据所述油相的饱和度与所述水相的饱和度，确定下式(8-4)表示的所述有效水饱和度，
其中，
为所述油相的饱和度；
为所述水相的饱和度且
以及Swir为所述毛细管内的束缚水饱和度。
也就是说，根据上式(8-1)-(8-4)可得到油相的压力分布场
步骤S8103，根据所述油相的压力分布场及达西公式，确定所述油相的速度分布场。
具体地，将油相的压力分布场
代入下式(8-5)表示的达西公式，可确定油相的速度分布场，
其中，μo为油相的粘度；以及
为所述储层的渗透率。
步骤S8104，根据所述油相的对流扩散定律、所述速度分布场及所述油相的弥散系数，确定润湿反转损害储层的时空演化模拟方程。
其中，所述时空演化模拟方程用于模拟由润湿反转引起的储层损害特征的四维时空演化过程。
对于步骤S8104，所述确定润湿反转损害储层的时空演化模拟方程可包括：根据所述油相的对流扩散定律、所述速度分布场
及所述油相的弥散系数Do，确定下式表示的润湿反转损害储层的时空演化模拟方程，
其中，
为所述储层的孔隙度；
为所述油相的饱和度。
也就是说，根据公式(8-1)-(8-6)可确定所述油相的饱和度
(即润湿反转损害储层的时空演化模拟方程)。
综上所述，本发明创造性地根据待诊断井的预设区域内的储层的压力分布方程，确定所述储层中的水相的压力分布场与毛细管的压力分布场之间的关系；然后根据所述水相的压力分布场与所述毛细管的压力分布场之间的关系及所述毛细管的受力平衡条件，确定所述油相的压力分布场；接着根据所述油相的压力分布场及达西公式，确定所述油相的速度分布场，最后根据所述油相的对流扩散定律、所述速度分布场及所述油相的弥散系数，确定润湿反转损害储层的时空演化模拟方程。由此，通过所确定的时空演化模拟方程可定量模拟由润湿反转引起的储层损害特征的四维时空演化过程，从而对未发生储层损害的井进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义，以及对已损害井优化设计解堵措施、提高或恢复油井产量和水井注水能力，以及提高油藏数值模拟精度都具有十分重大意义。
图8-2是本发明一实施例提供的确定储层损害程度的方法的流程图。如图8-2所示，所述方法包括可步骤S8201-S8202。
步骤S8201，基于根据所述的润湿反转损害储层的建模方法建立的时空演化模拟方程，确定待诊断井的预设区域内的储层中的油相的饱和度。
对于上述公式(8-6)所示的润湿反转损害储层的时空演化模拟方程的求取，可参见上述实施例1求解沉积颗粒的体积浓度的过程，由此可计算得到油相的饱和度
通过上述方法可计算得到油相的饱和度
由于通过上述润湿反转损害储层的建模方法建立的时空演化模拟方程综合考虑了润湿反转时多种物理化学因素对储层损害的影响，由此通过该步骤S8201求解得到的油相的饱和度非常精确。
步骤S8202，基于所确定的所述油相的饱和度，确定表征所述储层的损害程度的特征参数。
对于步骤S8202，所述特征参数可为所述储层的相对渗透率。相应地，所述确定所述储层的损害程度的特征参数可包括：基于所述油相的饱和度
及下式(8-7)表示的所述油相的相对渗透率与饱和度之间的关系，确定所述油相的相对渗透率
其中，α1、α2、α3、α4、α5为常数，
进一步，可根据相对渗透率
确定所述油相的渗透率
在一实施例中，所述特征参数可为所述储层的渗透率损害率。
相应地，所述确定表征待诊断井的预设区域内的储层的损害程度的特征参数可包括：基于所述储层的渗透率
及公式(8-8)，计算所述储层的渗透率损害率
其中，
为
的最大值。
在另一实施例中，所述特征参数可为所述储层的表皮系数。相应地，所述确定表征待诊断井的预设区域内的储层的损害程度的特征参数可包括：基于所述储层的渗透率
及公式(8-9)，计算所述储层的表皮系数
其中，
为所述储层的渗透率的初始值；以及
rw为所述待诊断井的井筒半径，以及rsw为所述储层的损害半径。
通过该步骤S8202得到的特征参数(例如所述储层的渗透率
与表皮系数
)是时空演化4D定量模拟的结果(如图8-3所示)。更具体地，图8-4示出了由储层渗透率损害率表征的在第365天润湿反转损害储层的半径的示意图(如箭头所指示的半径)，相关工作人员可通过该图8-4直观地确认储层被损害的程度。因此，可根据渗透率或表皮系数的演化特点进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义。
综上所述，本发明创造性地通过所确定的时空演化模拟方程可计算油相的饱和度，再基于所确定的所述油相的饱和度，确定表征待诊断井的预设区域内的储层的损害程度的特征参数(例如所述储层的渗透率和/或表皮系数)，由此可定量模拟由润湿反转引起的储层损害特征的四维时空演化过程，从而对未发生储层损害的井进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义，以及对已损害井优化设计解堵措施、提高或恢复油井产量和水井注水能力，以及提高油藏数值模拟精度都具有十分重大意义。
实施例9-乳化
在储层中的多孔介质中，油、水两相间的低界面张力和高机械剪切力是引起乳状液形成的主要因素。因此，基于储层温度场对油相粘度的影响和流体的乳化条件，确定形成的乳化液滴的半径；然后基于储层的孔径分布函数与所述乳化液滴的半径，确定堵塞概率的时空演化控制唯象模型；接着，再结合堵塞概率和渗透率等储层损害特征参数间的关系，就可诊断渗透率等储层损害特征参数的时空场分布。
图9-1是本发明一实施例提供的乳化堵塞损害储层的建模方法的流程图。所述建模方法可包括步骤S9101-S9104。
步骤S9101，确定待诊断井的预设区域内的储层中的流体的达西表观速度。
其中，所述待诊断井可例如为采油井。
对于步骤S9101，所述确定储层中的流体的速度可包括：建立所述流体进入所述储层的压力传导方程；以及根据所述压力传导方程及达西公式，确定所述流体的达西表观速度。
具体确定过程可参见上述实施例2中的确定达西表观速度的过程(即上式(2-1)与(2-2)及其相关描述)。
步骤S9102，根据所述储层的温度场与所述储层中的油相的粘度与温度之间的函数关系，确定所述油相的粘度。
在执行步骤S9102之前，所述建模方法还可包括：根据所述流体的导热系数、所述储层的热扩散系数、外来流体的平均流速及所述储层的热平衡方程，确定所述储层的温度场。
具体地，根据所述流体的导热系数Df、所述储层的热扩散系数Dl、外来流体的平均流速um及下式(9-1)表示的所述储层的热平衡方程，确定所述储层的温度场
结合初始条件以及边界条件，通过该公式(9-1)可以计算所述储层的不同时间-空间下的温度分布(即温度场)，具体计算过程可详见下文描述。
对于步骤S9102，根据所述储层的温度场
及下式(9-2)表示的所述油相的粘度与温度之间的函数关系，确定所述油相的粘度
其中，γAPI为所述油相的重力参数，a、b均为常数；以及
为所述储层内的任意点的空间位置。可根据公式(9-1)与(9-2)计算得到所述储层的不同时间-空间条件下的油相的粘度分布。由于油相的粘度分布受储层的温度分布的影响较大，故考虑了储层的温度场的油相粘度的计算结果更为准确，根据该油相粘度可更精确地模拟乳化堵塞对储层造成的损害。
步骤S9103，根据所述流体的达西表观速度、所述油相的粘度及所述流体的乳化条件，确定由所述流体乳化形成的乳化液滴的半径。
其中，所述流体的乳化条件可为下式(9-3)表示的临界条件：
其中，μw为所述流体中的水相的粘度；μo为所述油相的粘度；Kw为事所述水相的渗透率；Ko为所述油相的渗透率；ρw为所述水相的密度；ρo为所述油相的密度；g为重力加速度；
为所述油相的粘度；
为所述流体的达西表观速度；以及
为所述储层内的任意点的空间位置。
也就是说，如果不等式(9-3)成立，则说明储层中发生了油水乳化。
对于步骤S9103，所述确定由所述流体乳化形成的乳化液滴的半径可包括：根据所述流体的达西表观速度
所述油相的粘度
及所述流体的乳化条件，确定下式(9-4)表示的所述乳化液滴的半径：
其中，σ为油水界面张力；
为所述储层的孔径的均值；φ为所述储层的孔隙度；K为所述储层的渗透率；ρ为所述油相的密度；以及iΔt为第i个时间增量，i为非负整数。
步骤S9104，根据所述储层的孔隙的孔径分布函数及所述乳化液滴的半径，确定乳化堵塞损害储层的时空演化模拟方程。
其中，所述时空演化模拟方程可用于模拟由乳化堵塞引起的储层损害特征的四维时空演化过程。
对于步骤S9104，所述确定乳化堵塞损害储层的时空演化模拟方程可包括：根据所述储层的孔隙的孔径分布函数N(λ,μs,σs)及所述乳化液滴的半径
确定下式(9-5)表示的所述乳化堵塞损害储层的时空演化模拟方程：
其中，
为所述储层的堵塞概率；μs、σs分别为第一孔径分布特征参数与第二孔径分布特征参数；iΔt为第i个时间增量，i为非负整数；以及
为所述储层内的任意点的空间位置。
在一实施例中，可将所述储层中的孔隙的孔径分布函数N(λ)近似为下式 (9-6)表示的log-normal函数：
关于上式中的μs、σs可具体通过以下过程进行计算：根据现场数据可计算得到孔径均值和标准差分别为
(即
)与
进而可根据孔径均值和标准差得到下式(9-7) 表示的μs、σs的表达式：
由于孔径均值和标准差是已知量，故可得到相应的μs、σs，然后将μs、σs代入上式(9-6)可得出孔径分布函数的具体形式。
对于某一半径为λo的乳化液滴，只有分布函数中小于λo的那部分孔径才会被乳化液滴堵塞，因此，上式(9-5)中的堵塞概率β(即0到λo的累积分布)可具体表示为下式(9-8)：
其中，erfc()即为余误差函数：
综上所述，本发明创造性地确定待诊断井的预设区域内的储层中的流体的达西表观速度；根据所述储层的温度场与所述储层中的油相的粘度与温度之间的函数关系，确定所述油相的粘度；根据所述流体的乳化条件，确定由所述流体乳化形成的乳化液滴的半径；根据所述储层的孔隙的孔径分布函数及所述乳化液滴的半径，确定乳化堵塞损害储层的时空演化模拟方程。由此，通过所确定的时空演化模拟方程可定量模拟由乳化堵塞引起的储层损害特征的四维时空演化过程，从而对未发生储层损害的井进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义，以及对已损害井优化设计解堵措施、提高或恢复油井产量和水井注水能力，以及提高油藏数值模拟精度都具有十分重大意义。
图9-2是本发明一实施例提供的确定储层损害程度的方法的流程图。如图9-2所示，所述确定储层损害程度的方法可包括步骤S9201-S9202。
步骤S9201，基于所述的乳化堵塞损害储层的建模方法建立的时空演化模拟方程，确定所述储层的堵塞概率。
对于上述公式(9-8)所示的乳化堵塞损害储层的时空演化模拟方程的求取，需要根据公式(9-1)计算得到
具体求解过程可参见上述实施例1求解沉积颗粒的体积浓度的过程，于此不再进行赘述。
在通过上述方法计算得到所述储层的温度场
之后，再根据上述公式(9-2)、(9-4)、(9-8)可计算得到所述储层的堵塞概率，由此通过上述乳化堵塞损害储层的建模方法建立的时空演化模拟方程综合考虑了乳化堵塞时多种物理化学因素对储层损害的影响，由此通过该步骤S9201求解得到的储层的堵塞概率非常精确。
步骤S9202，基于所确定的所述储层的堵塞概率，确定表征待诊断井的预设区域内的储层的损害程度的特征参数。
其中，所述特征参数可为所述储层的渗透率。
对于步骤S9202，所述确定表征待诊断井的预设区域内的储层的损害程度的特征参数可包括：基于所确定的所述储层的堵塞概率
及下式 (9-9)，计算所述储层的渗透率
其中，μs为第一孔径分布特征参数；
为所述乳化液滴的半径；mK为第一经验值；以及n为时间增量Δt的总个数。
其中，所述特征参数可为所述储层的表皮系数。
对于步骤S9202，所述确定表征待诊断井的预设区域内的储层的损害程度的特征参数可包括：基于所确定的所述储层的堵塞概率
及下式，计算所述储层的渗透率
以及基于所述储层的渗透率
及公式(9-10)，计算所述储层的表皮系数
其中，μs为第一孔径分布特征参数；
为所述乳化液滴的半径；mK为第一经验值；n为时间增量Δt的总个数；rw为所述待诊断井的井筒半径，以及rsw为所述储层的损害半径。
通过该步骤S9202得到的特征参数(例如所述储层的渗透率
与表皮系数
)是时空演化4D定量模拟的结果(如图9-3所示)。更具体地，图9-4示出了由储层渗透率损害率(基于所述储层的渗透率
及公式
确定所述储层的渗透率损害率I(ri,t)，其中
为
的最大值)表征的在第40天乳化损害储层的半径的示意图(如箭头所指示的半径)，相关工作人员可通过该图9-4直观地确认储层被损害的程度。因此，可根据渗透率或表皮系数的演化特点进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义。
综上所述，本发明创造性地通过所确定的时空演化模拟方程可计算储层的堵塞概率，然后根据所确定的堵塞概率确定表征待诊断井的预设区域内的储层的损害程度的特征参数(例如所述储层的渗透率和/或表皮系数)，由此可定量模拟由乳化堵塞引起的储层损害特征的四维时空演化过程，从而对未发生储层损害的井进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义，以及对已损害井优化设计解堵措施、提高或恢复油井产量和水井注水能力，以及提高油藏数值模拟精度都具有十分重大意义。
实施例10-有机垢
在油田开发过程中，油层中的压力平衡状态被破坏，导致原油中的部分轻质组分溢出、部分重质组分沉淀，所述沉淀形成有机垢，进而堵塞油气流动通道，对储层导致严重损害。通常而言，沥青质等有机垢在原油中的赋存形式为部分溶解与部分胶体悬浮，沥青质等有机垢的胶体析出与其溶解能力有关。因此，本发明各实施例的核心是要建立有机垢的溶解能力随储层压力变化的动力学模型。具体地，基于压力对有机垢溶解能力的影响等建立有机垢影响待诊断井周围的储层中的有机垢颗粒分布的时空演化控制唯象模型，再结合有机垢颗粒分布与储层的孔隙度和渗透率等储层损害特征参数间的关系，就可诊断渗透率等储层损害特征参数的时空场分布。
图10-1是本发明一实施例提供的有机垢损害储层的建模方法的流程图。如图10-1所示，所述建模方法包括以下步骤S10101-S10104。
步骤S10101，确定位于待诊断井的预设区域内的储层的压力。
首先，可建立如图10-2所示的储层(或油藏)圆柱模型，即以所述待诊断井(例如采油井)的中心轴O为参照中心轴，将所述储层划分为多个圆柱壳，且每个圆柱壳的厚度为预设厚度dr(例如，该厚度的值非常小)。
由流体的流动对称性可知，每个圆柱壳内的压力大小相等(例如，第i 个圆柱壳内的压力均为P(ri))。根据下式(10-1)表示的所述储层的压力的连续性方程、内边界条件P(rw)×Pw(即所述待诊断井的井壁到所述参照中心轴的水平距离为rw，在所述带诊断井的井壁处的压力为Pw)及外边界条件 P(re)＝Pe(即所述储层的最***所在的圆柱壳到所述参照中心轴的水平距离为re，所述储层的最***所在的圆柱壳处的压力为Pe)，可得到P(ri)的表达式(10-2)，
结合上式(10-2)与下式(10-3)表示的裘比公式，可得到下式(10-4) 表示的第i个圆柱壳处的储层压力P(ri，t)，
其中，Qo(t)为所述待诊断井的原油产量；K为所述储层的渗透率；H为所述储层的厚度；以及μ为所述储层内的流体的粘度。
对于步骤S10101，所述确定位于待诊断井的预设区域内的储层的压力可包括：在所述储层被划分为以所述待诊断井的中心轴为参照中心轴且具有预设厚度的多个圆柱壳的情况下，根据所述储层的压力的连续性方程及裘比公式，确定下式(10-4)表示的所述储层的压力P(ri)，
其中，ri为所述多个圆柱壳中的第i个圆柱壳距所述参照中心轴的平均水平距离；re为所述储层的油藏半径；Qo(t)为所述待诊断井的原油产量；K为所述储层的渗透率；H为所述储层的厚度；以及μ为所述储层内的流体的粘度。
步骤S10102，根据所述储层的泡点压力、所述储层产出的原油在所述泡点压力下的摩尔体积、所述原油的溶解度参数、所述有机垢在所述原油中的溶解度参数及所述有机垢的摩尔体积，确定所述有机垢在所述原油中的最大溶解量随所述储层的压力变化的第一关系式。
其中，所述原油的溶解度参数包括以下方式获取：根据所述多个预设组分的沸点温度、临界温度与摩尔体积及所述储层的温度，确定所述多个预设组分的溶解度参数，其中，所述多个预设组分为具有多个预设碳含量的沥青质(例如C7、C8等C7+沥青质)；以及根据所述原油中的多个预设组分的溶解度参数与体积分数，确定所述原油的溶解度参数。
具体地，单位体积内原油分子内聚能量(即原油溶解度参数δL)的行为最复杂，需要先单独对各组分(例如，具有多个预设碳含量的沥青质(例如 C7、C8等C7+沥青质))的溶解度参数δi(P)进行求解。对于组分i，有：
其中：
其中，Tci、Tbi分别为所述组分i的临界温度与沸点温度；Vi为所述组分i的摩尔体积(Vi(P)＝xiV(P)，其中，xi为所述组分i的摩尔分数；V(P)可通过
计算得到，其中a、b分别为第一经验系数与第二经验系数)；T为所述储层的温度；R为气体常量。
然后，结合上式(13-3)并利用
计算得到δL，其中φi为所述组分i的体积分数(其可从油品物性分析数据得到)；以及n为所述组分的数目。
对于步骤S10102，所述确定所述有机垢在所述原油中的最大溶解量随所述储层的压力变化的第一关系式包括以下步骤S10301-S10302，如图10-3 所示。
步骤S10301，根据所述储层的泡点压力及所述储层产出的原油在所述泡点压力下的摩尔体积，确定所述原油在所述储层的压力下的摩尔体积。
对于步骤S10301，所述确定所述原油在所述储层的压力下的摩尔体积可包括：根据所述储层的泡点压力Pb及所述储层产出的原油在所述泡点压力下的摩尔体积VLb，确定下式(10-6)表示的所述原油在所述储层的压力P(ri，t) 下的摩尔体积VL(P(ri，t))，
其中，Cf为所述原油的压缩系数(例如，(10～140)×10-4MPa-1)；以及 Bo为所述原油的压缩系数(通常为1.0～1.2)。
当压力P高于泡点压力时，是纯压缩过程，压力增大则体积减小；当压力P低于泡点压力时，一方面压力减小、体积增大，另一方面原油解析气相、体积缩小，在某个压力点形成一个极值。
步骤S10302，根据所述原油在所述储层的压力下的摩尔体积、所述原油的溶解度参数、所述有机垢在所述原油中的溶解度参数及所述有机垢的摩尔体积，确定所述第一关系式。
对于步骤S10302，所述确定所述第一关系式可包括：根据所述原油在所述储层的压力下的摩尔体积VL(P(ri，t))、所述原油的溶解度参数δL(P(ri，t))、所述有机垢在所述原油中的溶解度参数δa及所述有机垢的摩尔体积Va，确定下式(10-7)表示的所述第一关系式，
其中，exp{x}是以自然常数e为底的指数函数，T为所述储层的温度； R为气体常量；以及(φa)max的单位可为％。在一实施例中，所述可按照以下方式进行取值：δa＝9.99×(1-5.94×10-4T)。
步骤S10103，根据所述有机垢中的有机垢颗粒的分布函数及所述原油的摩尔数，确定所述原油中的有机垢颗粒的摩尔数随所述储层的压力与所述有机垢在所述原油中的最大溶解量两者变化的第二关系式。
其中，所述分布函数为粒径大于预设粒径的有机垢颗粒的摩尔数占有机垢颗粒的总摩尔数的比例函数。
对于步骤S10103，所述确定所述原油中的有机垢颗粒的摩尔数随所述储层的压力变化的第二关系式可包括：根据所述有机垢中的有机垢颗粒的分布函数ftrap(Ro)及所述原油的摩尔数ηo，确定下式(10-8)表示的所述原油中的有机垢颗粒的摩尔数η(P(ri，t)，(φa)max(P(ri，t)))随所述储层的压力P(ri，t)与所述有机垢在所述原油中的最大溶解量(φa)max(P(ri，t))两者变化的第二关系式，
其中，ηo为所述原油的摩尔数；φa为所述原油中的有机垢的总含量；以及Rp为所述预设粒径。
下面对上式中的各个参数进行解释和说明。所述原油的摩尔数
其中，所述ρo为所述原油的密度；Qo(t)为所述待诊断井的原油产量；以及MWo为原油的平均摩尔质量。ηo[φa-(φa)max(P(ri，t))]表示在t时刻第i个圆柱壳内的有机垢颗粒的摩尔数变化速率；对于分布函数ftrap(Rp)，
其中f(r)为有机垢(例如沥青质)中的有机垢颗粒(例如沥青质颗粒)的密度分布函数(其可为正态分布函数)。由于ftrap(Rp)表示粒径大于所述预设粒径(例如，所述储层的孔隙的平均孔径)的有机垢颗粒的摩尔数占比，故上式(10-8)表示第i个圆柱壳内的粒径大于所述预设粒径的有机垢颗粒的摩尔数。
步骤S10104，根据所述第一关系式、所述第二关系式及所述储层的压力，确定有机垢损害储层的时空演化模拟方程。
其中，所述时空演化模拟方程用于模拟由有机垢引起的储层损害特征的四维时空演化过程。
具体地，将公式(10-4)与(10-8)代入公式(10-7)，可确定所述有机垢损害储层的时空演化模拟方程。由此，该时空演化模拟方程的具体形式非常复杂，于此不对其进行展示。也就是说，所述有机垢损害储层的时空演化模拟方程相当于由公式(10-4)、(10-7)-(10-8)组成的方程组。
在本申请的各个实施例中主要讨论的是关于沥青质的具体情况，也就是说，有机垢与沥青质两个概念可进行互换。
综上所述，本发明创造性地根据所述储层的泡点压力、所述储层产出的原油在所述泡点压力下的摩尔体积、所述原油的溶解度参数、所述有机垢在所述原油中的溶解度参数及所述有机垢的摩尔体积，确定所述有机垢在所述原油中的最大溶解量随所述储层的压力变化的第一关系式；根据所述有机垢中的有机垢颗粒的分布函数及所述原油的摩尔数，确定所述原油中的有机垢颗粒的摩尔数随所述储层的压力与所述有机垢在所述原油中的最大溶解量两者变化的第二关系式；根据所述第一关系式、所述第二关系式及所述储层的压力，确定有机垢损害储层的时空演化模拟方程。由此，通过所确定的时空演化模拟方程可定量模拟由有机垢引起的储层损害特征的四维时空演化过程，从而对未发生储层损害的井进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义，以及对已损害井优化设计解堵措施、提高或恢复油井产量和水井注水能力，以及提高油藏数值模拟精度都具有十分重大意义。
图10-4是本发明一实施例提供的确定储层损害程度的方法的流程图。如图10-4所示，所述方法可包括步骤S10401-S10402。
步骤S10401，基于所述的有机垢损害储层的建模方法建立的时空演化模拟方程，确定所述原油中的有机垢颗粒的摩尔数。
对于上述公式(10-1)所示的压力方程，可参见上述实施例1求解沉积颗粒的体积浓度的过程来求解所述储层的压力P(ri，t)及所述原油中的有机垢颗粒的摩尔数η(P(ri，t)，(φa)max(P(ri，t)))。
通过上述方法计算得到所述储层的压力P(ri，t)及所述原油中的有机垢颗粒的摩尔数η(P(ri，t)，(φa)max(P(ri，t)))，由于通过上述有机垢损害储层的建模方法建立的时空演化模拟方程综合考虑了有机垢颗粒堵塞储层时多种物理化学因素对储层损害的影响，由此通过该步骤S10401求解得到的所述原油中的有机垢颗粒的摩尔数非常精确。
步骤S10402，基于所述原油中的有机垢颗粒的摩尔数，确定表征待诊断井的预设区域内的储层的损害程度的特征参数。
在一实施例中，所述特征参数为所述储层的渗透率。
对于步骤S10402，所述确定表征待诊断井的预设区域内的储层的损害程度的特征参数可包括：基于所述原油中的有机垢颗粒的摩尔数η(P(ri，t)，(φa)max(P(ri，t)))及公式
确定所述储层的孔隙度φ(ri，t)，其中，dr为所述圆柱壳的预设厚度；基于所述储层的孔隙度φ(ri，t)及公式
确定所述储层的渗透率K(ri,t)。
其中，φ0为孔隙度的初始值；mK为第二经验值；以及K0(ri)为所述储层的渗透率的初始值。
在一实施例中，所述特征参数为所述储层的滤失系数。
对于步骤S10402，所述确定表征待诊断井的预设区域内的储层的损害程度的特征参数可包括：基于所述原油中的有机垢颗粒的摩尔数η(P(ri，t)，(φa)max(P(ri，t)))及公式
确定所述储层的孔隙度φ(ri，t)，其中，dr为所述圆柱壳的预设厚度；以及基于所述储层的孔隙度φ(ri，t)及公式
确定所述储层的滤失系数 k(ri,t)。
其中，φ0为孔隙度的初始值；mk为第一经验值；以及k0(ri)所述储层的滤失系数的初始值。
在一实施例中，所述特征参数为所述储层的表皮系数。
对于步骤S10402，所述确定表征待诊断井的预设区域内的储层的损害程度的特征参数包括：基于所述原油中的有机垢颗粒的摩尔数η(P(ri，t)，(φa)max(P(ri，t)))及公式
确定所述储层的孔隙度φ(ri，t)；基于所述储层的孔隙度φ(ri，t)及公式
确定所述储层的渗透率K(ri,t)；以及基于所述储层的渗透率K(ri,t)及公式
确定所述储层的表皮系数S(ri,t)。
其中，φ0为孔隙度的初始值；mK为第二经验值；K0(ri)为所述储层的渗透率的初始值；
rw为所述待诊断井的井筒半径；以及rsw为所述储层的损害半径。
通过该步骤S10402得到的特征参数(例如所述储层的渗透率K(ri,t)与表皮系数S(ri,t))是时空演化4D定量模拟的结果(如图10-5所示)。更具体地，图10-6示出了由储层渗透率损害率(基于所述储层的渗透率K(ri,t)及公式
确定所述储层的渗透率损害率I(ri,t)，其中Kmax(ri,t)为 K(ri,t)的最大值)表征的在第40天有机垢损害储层的半径的示意图(如箭头所指示的半径)，相关工作人员可通过该图10-6直观地确认储层被损害的程度。因此，可根据渗透率或表皮系数的演化特点进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义。
综上所述，本发明创造性地根据所述的有机垢损害储层的建模方法建立的时空演化模拟方程，确定所述原油中的有机垢颗粒的摩尔数；以及基于所述原油中的有机垢颗粒的摩尔数，确定表征待诊断井的预设区域内的储层的损害程度的特征参数。由此可定量模拟由黏土膨胀引起的储层损害特征的四维时空演化过程，从而对未发生储层损害的井进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义，以及对已损害井优化设计解堵措施、提高或恢复油井产量和水井注水能力，以及提高油藏数值模拟精度都具有十分重大意义。
实施例11-贾敏效应
在低渗透储层内的油、水两相流动过程中，通常含有大量分散的油滴和气泡。当这些个体运动到储层中的狭小孔喉(例如孔径为2.5μm的孔隙)时，由于这些个体的直径大于孔喉直径，其流动受到阻力，迫使油滴或气泡等发生变形，故将油滴或气泡等在直径分布不均匀的储层中运动时受到的附加的毛细管阻力的现象称为贾敏效应，如图1所示。与水锁效应相比，微小孔喉中受到阻力被迫变形的是油滴或者气泡，此时润湿相是水，非润湿相是油滴或气泡。
通常情况下，产生贾敏效应需同时满足以下条件：首先，原油应该处于分散状态，不能以连续的油流形式运动；其次，地质储层的孔道直径应该比较细小。当分散的油滴运动到油层狭小孔道时，由于其直径大于孔道直径，就会受到阻力而致使其无法继续运动，或者油滴在变直径孔道中流动时会产生毛细管效应，这样就会导致贾敏损害(因贾敏效应导致储层的渗透率下降的情形)。
具体地，如图1所示，在储层左右两侧的驱替压力P的作用下，液滴从左向右运动。当油滴运动到孔隙1处时，因孔隙1的阻挡，对于油滴前端 (2R1＝λ，其中R1为前端液面的曲率半径，λ为孔隙1的孔径)，其受到泡前压力P1的作用；对于油滴后端(其直径
其中R2为后端液面的曲率半径，
为孔隙的平均孔径)，其受到P3的作用；以及对于油滴内部，其泡内压力为P2，并且，P1、P2、P3满足以下两个等式：
由此，因贾敏效应油滴受到的附加阻力为：
考虑液滴前端
液滴后端
其中σ为液滴的表面张力。当驱替压力(即孔喉外压力差)至少达到Pc时，气泡或液滴才能通过孔喉，否则就会被堵住。
贾敏效应受储层的孔隙结构、岩性、物性、侵入流体等多种因素的影响。贾敏损害同岩石的孔隙介质的几何特性关系密切，不同的孔喉结构分布模式和复杂性会导致水润湿相在岩石中的展布模式显著变化，从而影响储层的渗透率。因此，本发明各实施例的核心是要建立储层中的孔隙内的水相饱和度变化的动力学模型(即储层中的水相运动方程及渗透率的分布方程)。具体地，基于储层中的孔隙内的流体的对流扩散关系(即质量平衡方程)等建立储层中的水相运动方程；基于孔隙的孔径分布特征及预设渗透率模型，建立渗透率的分布方程；然后根据所述水相运动方程及所述渗透率的分布方程，确定贾敏效应影响待诊断井周围的储层中的渗透率分布的时空演化控制唯象模型，由此就可诊断渗透率等储层损害特征参数的时空场分布。
图11-1是本发明一实施例提供的贾敏效应损害储层的建模方法的流程图。所述建模方法可包括步骤S11101-S11104。
步骤S11101，确定待诊断井的预设区域内的储层中的流体的达西表观速度。
其中，所述储层的渗透率低于预设渗透率；以及所述待诊断井可例如为采油井。
对于步骤S11101，所述确定储层中的流体的速度可包括：建立所述流体进入所述储层的压力传导方程；以及根据所述压力传导方程及达西公式，确定所述流体的达西表观速度。
具体确定过程可参见上述实施例2中的确定达西表观速度的过程(即上式(2-1)与(2-2)及其相关描述)。
步骤S11102，根据所述流体的达西表观速度与所述流体中的水分子的扩散系数，建立所述储层的水相运动方程。
在储层条件下，所述储层中的孔隙内的不同位置的含水量满足质量守恒方程。其中，外来水相在所述储层内的运动主要由对流与扩散两个过程决定。具体地，对于步骤S11102，所述建立所述储层的水相运动方程可包括：根据所述流体的达西表观速度u与所述水分子的扩散系数Dw，建立下式表示的所述质量平衡方程：
其中，φ0为所述储层的孔隙度的初始值；
为所述储层内的孔隙被水相占据的绝对孔隙度；以及
为所述储层内的任意点的空间位置(例如，以所述待诊断井的中心为原点)。
根据所述质量平衡方程及所述储层的水相饱和度的时空分布函数
建立下式(11-1)表示的所述水相运动方程：
所述水相运动方程的初始条件为
以及所述水相运动方程的边界条件为
(也就是说，在注水井的井壁处的储层孔隙完全被水充满，即孔隙中的水相饱和度为1)。其中，φ0为所述储层的孔隙度的初始值；rw为所述待诊断井的井筒半径；以及Swc为所述储层中的束缚水饱和度。
步骤S11103，根据所述储层的孔隙的孔径分布特征及所述储层的预设渗透率模型，建立所述储层的渗透率的分布方程。
对于步骤S11103，如图2所示，所述建立所述储层的渗透率的分布方程可包括步骤S11201-S11202。
步骤S11201，根据所述储层的孔隙的孔径分布特征，确定孔径为λ的孔隙的体积密度函数与所述储层的水相饱和度的孔径分布方程。
储层岩石的孔隙结构具有强烈的不规则性，受到多种因素的控制，使得经典几何学为基础的孔隙结构模型与真实的孔隙相差太远。为了定量描述储层的孔隙结构，采用分形理论来研究孔隙结构中两相流的贾敏效应，并忽略水相在孔隙通道内表面的润湿薄膜影响。根据分形的几何原理，若储层的孔径分布具有分形特征，则储层中孔径大于λ的孔隙的数量N(>λ)与λ有如下幂函数关系：
其中，D为孔隙的分形维数(2<D<3)。
在所述储层的孔隙的孔径分布特征为所述储层中的孔径大于λ的孔隙的数量N(＞λ)满足上式
的情况下，如图11-3所示，所述步骤 S11201中的确定孔径为λ的孔隙的体积密度函数可包括步骤S11301-S11302。
步骤S11301，根据所述储层中的孔径大于λ的孔隙的数量N(＞λ)，确定所述储层中的孔隙的总体积为
具体地，根据上式
可得到储层中的孔隙的总数目N(＞λmin)及孔径小于λ的孔隙的数目：
根据公式(11-2)与(11-3)可得到下式(11-4)
上式(11-4)中的f(λ)为储层的孔径分布密度函数，则孔径小于λpc的孔隙的数量N(＜λpc)与λ之间的关系为下式(11-5)所表示的幂函数关系：
上述相关公式中的λ、λmin、λmax与λpc分别为所述孔隙的孔径、最小孔径、最大孔径(λmin与λmax可以由平均孔径与孔径分布标准差得出，一般而言有
)与发生贾敏效应的最大孔径(即非水相受到所述最小毛细管阻力处的孔隙的特定孔径)；a为比例常数。
接着，由公式(11-5)可得到储层的孔径分布密度函数f(λ)满足下式 (11-6)：
式中a'＝-Da为比例常数。
由上式(11-6)表示的所述孔径分布密度函数可以得到储层中的孔隙的总体积分形表达式：
其中α为与孔隙的形状有关的常数(孔隙的形状为立方体时α＝1，或者孔隙的形状为球体时α＝π/6)，通过积分可得：
类似地，所述储层中的孔径小于λ的孔隙的体积
步骤S11302，根据所述储层中的孔隙的总体积Φmax及所述储层中的孔径小于λ的孔隙的体积
确定所述孔径为λ的孔隙的体积密度函数为
其中，D为所述孔隙的分形维数；以及λ、λmin与λmax分别为所述孔隙的孔径、最小孔径与最大孔径；以及A＝αa'/(3-D)(常数)。
步骤S11202，根据所述预设渗透率模型、所述孔径为λ的孔隙的体积密度函数及所述水相饱和度的孔径分布方程，建立所述储层的渗透率的分布方程。
如图11-4所示，步骤S11202中的确定所述储层的水相饱和度的孔径分布方程可包括步骤S11401-S11402。
步骤S11401，根据所述储层中的孔径大于λ的孔隙的数量N(＞λ)，确定被非水相占据的孔隙的体积为
其中，λpc为非水相受到所述最小毛细管阻力处的孔隙的特定孔径(即发生贾敏效应时的最大孔喉直径)。
在储层为水湿(即亲水性)的情况下，向所述储层的两侧施加一个驱替压力P(从零逐渐增加到P)，根据贾敏效应产生条件，油滴或气泡可顺利通过孔径处于λmax到λpc范围之间的孔隙(即连通孔径由最大λmax到λpc逐次打开)。当压力P未到Pc时，λ大于λpc的孔喉完全被水占据，而λ小于λpc的孔喉还被油滴或气泡堵塞。
由于根据所述储层中的孔径大于λ的孔隙的数量N(＞λ)，可确定上式 (13-6)确定的孔径分布密度函数f(λ)，则进而可根据f(λ)得到压力P未到 Pc时的非水相占据的孔隙体积Φnw(λ)，
步骤S11402，根据所述储层中的孔隙的总体积Φmax与所述被非水相占据的孔隙的体积Φnw(λ)，确定下式表示的所述水相饱和度的孔径分布方程
其中，D为所述孔隙的分形维数；λ、λmin与λmax分别为所述孔隙的孔径、最小孔径与最大孔径；以及A＝αa'/(3-D)。
具体地，根据公式(11-8)与公式(11-10)可确定水相占据的孔隙体积为：
然后再结合公式(11-8) 可得到下式：
由于λmin<<λmax，
故上式可写为下式(11-11)，即可确定所述水相饱和度的孔径分布方程(11-11)，
根据流体力学理论，在压差ΔP作用下，粘度为μ的流体流过毛细管束总流量可以由Hagen-poiseuille方程描述，
上式中，λi为孔隙通道i的直径(即孔隙i的孔径)，Vi为孔隙通道i的体积， Vi＝Φ0MLmξi，Φ0为储层的孔隙度，ξi为孔隙通道i所占孔隙体积的分数，Lc为孔隙通道i的弯曲长度，Lm为孔隙通道i的直线长度，M为平均孔喉截面积
为孔隙通道i的迂曲度(同时可由经验公式
计算)，将上述各量代入上式得到，
结合达西定律，可得到毛管束模型的渗透率可表示为
在本发明实施例中，可将储层的渗透通道看做是多个毛管束累积起来的。由于孔隙的尺寸分布具有连续性，因此可以将所述毛管束模型的渗透率的表达式用积分形式写出：
在所述储层的预设渗透率模型满足
的情况下，所述建立所述储层的渗透率的分布方程可包括：根据所述储层的预设渗透率模型
所述孔径为λ的孔隙的体积密度函数dξ及所述水相饱和度的孔径分布方程，建立下式表示的所述渗透率的分布方程，
其中，所述压力传导方程的初始条件为K(t＝0)＝1。K(Sw)为渗透率关于含水饱和度的函数(如图11-5所示，该函数可简称为无因次渗透率)；φ0为储层的初始孔隙度，Sw为含水饱和度；D为所述孔隙的分形维数，其可由
计算(D越大说明小孔喉所占比例越大，常规储层λmax＞ 20μm，D<2.6，低渗储层λmax＜0.2μm，D>2.8)，τ为迂曲度，λmax为最大孔喉直径，λmin为最小孔喉直径)。
具体地，首先，逐渐增加驱替压力直至P，此时直径大于λ的孔喉被水水占据，被水相占据的绝对孔隙度为φw，则将公式(11-9)代入公式(11-12) 可得出所述储层的渗透率的孔径分布函数：
然后，根据式(11-11)中的λ与含水饱和度Sw(λ)之间的关系式将上述渗透率的孔径分布函数K(λ)中的变量λ替换为Sw，并结合φw＝φ0Sw，可得到公式 (11-13)(即得到无因次渗透率的表达式，该无因次渗透率为含水饱和度的函数)。
步骤S11104，根据所述渗透率的分布方程及所述水相运动方程，确定贾敏效应损害储层的时空演化模拟方程，其中所述时空演化模拟方程用于模拟由贾敏效应引起的储层损害特征的四维时空演化过程。
具体地，根据公式(11-1)可求解到所述储层的水相饱和度的时空分布函数
(如图11-6所示)，并将
代入公式(11-13)可得到所述储层的渗透率在四维时空上的分布形式，即贾敏效应损害储层的时空演化模拟方程。
综上所述，本发明创造性地确定待诊断井的预设区域内的储层中的流体的达西表观速度，其中所述储层的渗透率低于预设渗透率；根据所述流体的达西表观速度与所述流体中的水分子的扩散系数，建立所述储层的水相运动方程；建立所述储层的渗透率的分布方程；以及根据所述渗透率的分布方程及所述水相运动方程，确定贾敏效应损害储层的时空演化模拟方程。由此，通过所确定的时空演化模拟方程可定量模拟由贾敏效应引起的储层损害特征的四维时空演化过程，从而对未发生储层损害的井进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义，以及对已损害井优化设计解堵措施、提高或恢复油井产量和水井注水能力，以及提高油藏数值模拟精度都具有十分重大意义。
相应地，本发明另一实施例还提供一种确定储层损害程度的方法，所述方法可包括：基于所述的贾敏效应损害储层的建模方法建立的时空演化模拟方程，确定表征待诊断井的预设区域内的储层的损害程度的特征参数。
对于上述公式(11-13)所示的贾敏效应损害储层的时空演化模拟方程的求取，需要根据公式(11-1)计算得到
具体求解过程可参见上述实施例1求解沉积颗粒的体积浓度的过程，于此不再进行赘述。
在通过上述方法计算得到所述储层的水相饱和度
之后，再根据上述公式(11-13)可计算得到所述储层的渗透率
由此通过上述贾敏效应损害储层的建模方法建立的时空演化模拟方程综合考虑了贾敏损害时多种物理化学因素对储层损害的影响，由此通过本实施例求解得到的储层的渗透率非常精确。
在得到所述储层的渗透率的基础上，表征待诊断井的预设区域内的储层的损害程度的特征参数。
在一实施例中，所述特征参数可为所述储层的渗透率损害率。
相应地，所述确定表征待诊断井的预设区域内的储层的损害程度的特征参数可包括：
基于所述时空演化模拟方程，确定所述储层的渗透率
以及基于所述储层的渗透率
及公式(5-14)，确定所述储层的渗透率损害率
在另一实施例中，所述特征参数可为所述储层的表皮系数。
所述确定表征待诊断井的预设区域内的储层的损害程度的特征参数可包括：基于所述时空演化模拟方程，确定所述储层的渗透率
以及基于所述储层的渗透率
及公式(5-15)，确定所述储层的表皮系数
通过上述各个实施例得到的特征参数(例如所述储层的渗透率损害率 I(Sw)是时空演化4D定量模拟的结果(分别如图11-7所示)。因此，可根据渗透率或表皮系数的演化特点进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义。
实施例12-细菌
细菌的代谢由其细胞内各种酶的催化作用来控制，而酶的活性对温度极为敏感，例如过高或过低的温度都会使细胞内的酶失活。在注水过程中，会使井内与储层间形成温度变化区间，其严重影响细菌代谢速率，细菌附着在岩石表面会形成生物膜，从而降低储层的孔隙度。因此，本发明各实施例的核心是要建立储层内部流体中的营养物的表观浓度分布方程与温度之间的变化关系以及所述流体中的细菌的表观浓度分布方程与温度之间的变化关系。具体地，基于能量守恒、质量守恒、扩散关系等建立待诊断井周围的储层内流体中的营养物与细菌的浓度分布的时空演化控制唯象模型，再结合浓度分布和渗透率等储层损害特征参数间的关系，就可诊断渗透率等储层损害特征参数的时空场分布。
图12-1是本发明一实施例提供的细菌损害储层的建模方法的流程图。所述建模方法可包括步骤S12101-S12105。
在执行步骤S12101之前，所述建模方法还包括：确定待诊断井的预设区域内的储层中的流体的达西表观速度。
所述确定储层中的流体的速度可包括：建立所述流体进入所述储层的压力传导方程；以及根据所述压力传导方程及达西公式，确定所述流体的达西表观速度。
具体确定过程可参见上述实施例2中的确定达西表观速度的过程(即上式(2-1)与(2-2)及其相关描述)。
在执行步骤S12101之前，所述建模方法还包括：根据所述流体的达西表观速度、所述流体的导热系数与热扩散率、及能量守恒定理，确定所述储层的温度分布场。
在注水过程中，由于注入水温度与储层温度存在温度差，能量在储层岩石与流体之间以温度变化的形式传递。在这种情况下，可根据待诊断井的预设区域内的储层中的流体的达西表观速度
所述流体的导热系数Dcon与热扩散率Ddis、及能量守恒定理，建立矿场尺度储层温度分布数学模型，得到储层温度分布控制方程表达式：
其中，
为储层温度分布；所述热扩散率可采用导热率进行表示。
步骤S12101，根据所述储层的温度分布场及所述流体中的营养物的实际浓度，确定所述细菌的生长速率。
对于步骤S12101，所述确定所述细菌(包括附着在岩石表面的细菌与流体中的细菌)的生长速率可包括以下步骤S12201-S12202，如图12-2所示。
步骤S12201，根据所述储层的温度分布场及细菌最大生长速率公式，确定所述细菌的最大生长速率。
目前，通常采用含有活化能、频率因数等参数的平方根模型来描述细菌的最大生长速率，但该模型的现场适用性差(主要表现为预测结果不准确)。然而，发明人经过不断研究发现储层的温度是影响细菌最大生长速率和细菌在储层中的分布情况的主要因素。因此，在注水水源和细菌营养源充足的条件下，本实施例采用温度为主变量的细菌(包括附着在岩石表面的细菌与流体中的细菌)的最大生长速率公式，来模拟细菌在储层中最大生长速率的空间分布情况。
对于步骤S12201，所述确定所述细菌的最大生长速率包括：根据所述储层的温度分布场
及下式(12-2)表示的细菌最大生长速率公式，确定所述细菌的最大生长速率
其中，b1、c1分别为第一细菌生长经验参数与第二细菌生长经验参数； Tmax、Tmin分别为细菌生长的最大温度与最小温度；
步骤S12202，根据所述细菌的最大生长速率及所述流体中的营养物的实际浓度，确定所述细菌的生长速率。
考虑到细菌生长为不可逆一级反应或者细菌代谢会消耗营养物，对于步骤S12202，所述确定所述细菌的生长速率包括：根据所述细菌的最大生长速率
所述流体中的营养物的实际浓度Cnu、莫诺半生长系数kS及下式(12-3)，确定所述细菌的生长速率gactual，
步骤S12102，根据与所述流体中的细菌的表观浓度及所述细菌在所述岩石表面上的总量两者相关联的所述细菌在所述储层中的岩石表面上的附着量、所述细菌的生长速率与衰亡速率，确定所述细菌在所述岩石表面上的总量。
具体地，影响细菌在岩石表面附着量的主要因素为细菌吸附与解吸附速率，且细菌附着于储层岩石表面会形成生物膜，则可根据堵塞速率kclogging(其为常数，单位可为1/day)、解堵速率kdeclogging(其为常数，单位可为1/day)、所述流体中的细菌的表观浓度Cbacteriatran、所述细菌在所述岩石表面上的总量 Vbacteriatran及下式(12-4)，确定所述细菌在所述储层中的岩石表面上的附着量Cdeposition，
Cdeposition＝kcloggingCbacteriatran-kdecloggingVbacteriatran。(12-4)
附着在岩石表面的细菌会形成生物膜，从而降低储层的孔隙度。在本实施例中，岩石表面的细菌量主要由细菌生长、衰亡净增值与附着量两大因素决定。
对于步骤S12102，所述确定所述细菌在所述岩石表面上的总量可包括：根据所述细菌在所述储层中的岩石表面上的附着量Cdeposition、所述细菌的生长速率gactual与衰亡速率kdecay及下式，确定所述细菌在所述岩石表面上的总量Vbacteriatran，
上述(12-5)的左侧描述细菌在岩石表面的总量随时间的变化，右侧分别是细菌在岩石表面上代谢衰亡总生长值和细菌的沉积量。
步骤S12103，根据所述流体的达西表观速度、所述细菌的弥散系数、所述细菌的生长速率与衰亡速率及所述细菌在所述储层中的岩石表面上的附着量，建立所述流体中的细菌的表观浓度分布方程。
在本实施例中，主要考虑因浓度梯度而使得细菌发生扩散、因注入水而使得细菌发生对流运移、以及细菌的生长与衰亡，前四项因素导致细菌浓度的变化(其中包括因细菌附着于储层岩石上形成生物膜而导致浓度的减少)。由于细菌鞭毛振动的无规则运动在储层中的宏观作用效果不明显，模型中将鞭毛振动产生的无规则运动效果合并进入对流项；而细菌布朗运动作用被细菌扩散作用覆盖。此外，由于细菌趋化性作用对储层中细菌分布规律影响较小，该影响被对流效果覆盖，将其与对流项合并为一项。
对于步骤S12103，所述建立所述流体中的细菌的表观浓度分布方程可包括：根据所述流体的达西表观速度u、所述细菌的弥散系数Dsum、所述细菌的生长速率gactual与衰亡速率kdecay及所述细菌在所述储层中的岩石表面上的附着量Cdeposition，建立下式(12-6)表示的所述流体中的细菌的表观浓度分布方程，
其中，Cnutran为所述流体中的细菌的表观浓度分布。
步骤S12104，根据所述流体的达西表观速度、所述营养物的弥散系数、所述细菌在所述岩石表面上的总量及所述细菌的表观浓度，建立所述营养物的表观浓度分布方程。
在本实施例中，主要考虑因浓度梯度而使得营养物发生扩散与因注入水而使得营养物发生对流运移，这两项因素导致营养物浓度的变化。
对于步骤S12104，所述建立所述营养物的表观浓度分布方程可包括：根据所述流体的达西表观速度
所述营养物的弥散系数Dnusum、所述细菌在所述岩石表面上的总量Vbacteriatran及所述细菌的表观浓度Cbacteriatran，建立下式(12-7)表示的所述营养物的表观浓度分布方程，
其中，gactual为所述细菌的生长速率；Y为所述细菌的产率系数；以及 Cnutran为所述营养物的表观浓度分布。上式(12-7)的一维形式可写为
步骤S12105，根据所述营养物的表观浓度分布方程与所述流体中的细菌的表观浓度分布方程，确定细菌损害储层的时空演化模拟方程。
其中，所述时空演化模拟方程用于模拟由细菌引起的储层损害特征的四维时空演化过程。
具体地，根据上式(12-7)表示的所述营养物的表观浓度分布方程与上式(12-6)表示的所述流体中的细菌的表观浓度分布方程，再结合其他公式 (12-1)-(12-5)，可得到所述细菌损害储层的时空演化模拟方程。也就是说，所述细菌损害储层的时空演化模拟方程相当于由公式(12-1)-(12-7) 组成的方程组。
综上所述，本发明创造性地根据所述储层的温度分布场及所述流体中的营养物的实际浓度，确定所述细菌的生长速率；根据所述细菌在所述储层中的岩石表面上的附着量、所述细菌的生长速率与衰亡速率，确定所述细菌在所述岩石表面上的总量；根据所述流体的达西表观速度、所述细菌的弥散系数、所述细菌的生长速率与衰亡速率及所述细菌在所述储层中的岩石表面上的附着量，建立所述流体中的细菌的表观浓度分布方程；根据所述流体的达西表观速度、所述营养物的弥散系数、所述细菌在所述岩石表面上的总量及所述细菌的表观浓度，建立所述营养物的表观浓度分布方程；以及根据所述营养物的表观浓度分布方程与所述流体中的细菌的表观浓度分布方程，确定细菌损害储层的时空演化模拟方程。由此，通过所确定的时空演化模拟方程可定量模拟由细菌引起的储层损害特征的四维时空演化过程，从而对未发生储层损害的井进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义，以及对已损害井优化设计解堵措施、提高或恢复油井产量和水井注水能力，以及提高油藏数值模拟精度都具有十分重大意义。
图12-3是本发明一实施例提供的确定储层损害程度的方法的流程图。如图12-3所示，所述方法包括步骤S12301-S12302。
步骤S12301，基于所述的细菌损害储层的建模方法建立的时空演化模拟方程，确定所述细菌在所述岩石表面上的附着量。
对于上述公式(12-6)-(12-7)所示的细菌损害储层的时空演化模拟方程，可参见上述实施例1求解沉积颗粒的体积浓度的过程来求解所述细菌的附着量Cdeposition。
在通过上述方法计算得到所述细菌的附着量Cdeposition，由此通过上述细菌损害储层的建模方法建立的时空演化模拟方程综合考虑了细菌与营养物在流体中运移时多种物理化学因素对储层损害的影响，由此通过该步骤 S12301求解得到的细菌的附着(即沉积)量非常精确。
步骤S12302，基于所述细菌在所述岩石表面上的附着量，确定表征待诊断井的预设区域内的储层的损害程度的特征参数。
其中，所述特征参数为所述储层的渗透率。
对于步骤S12302，所述确定表征待诊断井的预设区域内的储层的损害程度的特征参数包括：基于所述细菌在所述储层中的岩石表面上的附着量 Cdeposition与所述细菌的密度ρ下式(12-8)，确定所述储层的渗透率
其中，φ0为所述储层的孔隙度的初始值；以及
为所述储层的渗透率的初始值。
其中，所述特征参数为所述储层的表皮系数，
对于步骤S12302，所述确定表征待诊断井的预设区域内的储层的损害程度的特征参数包括：基于述流体中的营养物的表观浓度
与实际浓度
及公式
确定所述储层的渗透率
以及基于所述储层的渗透率
及下式(12-9)，确定所述储层的表皮系数
其中，
为所述储层的渗透率的初始值，
rw为所述待诊断井的井筒半径，以及rsw为所述储层的损害半径。
通过该步骤S12302得到的特征参数(例如所述储层的渗透率
与表皮系数
)是时空演化4D定量模拟的结果(未示出)。因此，可根据渗透率或表皮系数的演化特点进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义。
综上所述，通过所确定的时空演化模拟方程可确定所述流体中的营养物的表观浓度，再基于所述流体中的营养物的表观浓度与实际浓度，可确定表征待诊断井的预设区域内的储层的损害程度的特征参数(例如所述储层的渗透率和/或表皮系数)，由此可定量模拟由细菌引起的储层损害特征的四维时空演化过程，从而对未发生储层损害的井进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义，以及对已损害井优化设计解堵措施、提高或恢复油井产量和水井注水能力，以及提高油藏数值模拟精度都具有十分重大意义。
实施例13-聚合物
根据聚合物在多孔介质表面的吸附机理可知，聚合物在储层上的吸附包括层吸附(如图(13-1A)-图(13-1C)为例)和桥接吸附(如图13-1D所示)。聚合物和介质壁面有着非常丰富的物理化学相互作用，吸附行为也是非常复杂。在聚合物吸附起初，聚合物发生层吸附作用而使得越来越多的分子链吸附在一起，形成吸附聚合物(即沉淀聚合物)；当越来越多的聚合物的分子链吸附在储层的同一位置时，此处的吸附聚合物的质量越来越大，当此位置上的吸附聚合物的质量大于临界质量(m＞mc)时，聚合物发生桥接吸附而使得更多的分子链吸附在一起，形成更大质量或更大分子链数目的吸附聚合物(即沉淀聚合物)。在实际应用中，吸附聚合物的质量很容易满足所述临界质量，故在下述各个实施例中，以发生桥接吸附的场景进行描述与说明。
聚合物吸附产生堵塞的本质是储层内的流体中的聚合物的运移和吸附。因此，本发明各实施例的核心是要根据质量守恒定律建立聚合物的运移和吸附的动力学模型。具体地，基于质量守恒、扩散关系等建立聚合物在待诊断井周围的储层中的浓度分布的时空演化控制唯象模型(该模型包含流体中的聚合物的浓度C)，再结合该浓度C和渗透率等储层损害特征参数间的关系，就可诊断渗透率等储层损害特征参数的时空场分布。
图13-2是本发明一实施例提供的聚合物损害储层的建模方法的流程图。如图13-2所示，所述建模方法可包括以下步骤S13201-S13204。
步骤S13201，确定储层中的流体的速度。
其中，所述储层位于待诊断井(例如注聚井)的预设区域内。
对于步骤S13201，所述确定储层中的含有所述聚合物的流体的速度包括：建立所述流体进入所述储层的压力传导方程；以及根据所述压力传导方程及达西公式，确定所述流体的速度。
具体确定过程可参见上述实施例1中的确定流体速度的过程(即上式 (1-1)与(1-2)及其相关描述)。
在t时刻，位于储层多孔介质中
位置处的单位控制体积的流体包含聚合物的分子链共
条，该单位控制体积以速度v在孔隙内流动，经过了时间Δt，经历位移
因层吸附与桥接吸附作用使得吸附在储层上的聚合物的分子链的条数变化，故该单位控制体积内包含聚合物的分子链变为
条。在上述吸附过程中，整个储层内的聚合物分子链的总条数不变。
在执行步骤S13202之前，所述建模方法还可包括：根据层吸附在所述储层上的聚合物的分子链的数目与桥接吸附在所述储层上的聚合物的分子链的数目，建立下式(13-1)表示的所述流体中的聚合物的分子链的数目分布方程，
其中，χ为在t时刻层吸附在所述储层的任意位置
处的聚合物的质量占比(即层吸附在所述储层的任意位置
处的聚合物的质量占在t时刻在该位置处的所有吸附聚合物的总质量的比例)；
m为吸附聚合物的分子链的平均质量；ρs为吸附聚合物的密度；
分别为发生层吸附的储层和桥接吸附的储层的比表面积；φ0为所述储层的孔隙度(在本发明各个实施例中可认为所述储层的孔隙度保持不变)；Γl、Γb分别为所述层吸附密度与所述桥接吸附密度；以及
为在t 时刻在所述储层的任意位置
处的所述流体中的聚合物的分子链的数目，
步骤S202，根据所述流体的速度、所述流体中的聚合物的分子链的数目分布方程及所述聚合物的扩散系数，建立所述流体中的聚合物的质量浓度的平衡方程及所述聚合物中吸附在所述储层上的吸附聚合物的分子链的占比分布方程。
其中，所述吸附聚合物的分子链的占比为所述吸附聚合物的分子链的数目占所述流体中的聚合物的分子链的初始数目的比例。
对于步骤S13202，所述建立所述流体中的聚合物的质量浓度的平衡方程及所述聚合物中吸附在所述储层上的吸附聚合物的分子链的占比分布方程可包括以下步骤S13301-S13303，如图13-3所示。
步骤S13301，根据所述流体中的聚合物的分子链的数目分布方程，确定所述吸附聚合物的质量浓度随所述层吸附密度与所述桥接吸附密度变化的函数关系。
上式(13-1)的左右两端均表示单位控制体积内发生聚合物吸附的分子链数，将(13-1)式两端乘以m(发生吸附的聚合物的质量)，并取极限Δt→0，则两端可表示单位控制体积内发生吸附的聚合物的分子质量(即质量浓度 Cd(Γl，Γb))，
其中，
表示发生层吸附的聚合物的质量浓度；
表示发生桥接吸附的聚合物质量浓度。
步骤S13302，根据所述流体的速度、所述聚合物的扩散系数及所述吸附聚合物的质量浓度随所述层吸附密度与所述桥接吸附密度变化的函数关系，建立所述流体中的聚合物的质量浓度的平衡方程。
对于聚合物的吸附过程，聚合物分子链遵循质量守恒定律，即流动的聚合物分子质量的变化等于吸附在储层的多孔介质表面的聚合物分子质量的负变化：
其中，
为聚合物的扩散流量，
D表示聚合物的扩散系数；
表示发生吸附的聚合物的质量；
为所述流体的速度；ρ为所述流体的密度；
为流动的聚合物分子在流体中的质量分数；以及φ0为所述储层的孔隙度的初始值。若考虑单位控制体积的质量变化(即质量浓度的变化)，并且由于孔隙内的微观变化(无需考虑孔隙度φ0)，可将上式换为下式(13-3)所示的质量浓度。
对于步骤S13302，所述建立所述流体中的聚合物的质量浓度的平衡方程可包括：根据所述流体的速度
所述聚合物的扩散系数D及所述吸附聚合物的质量浓度随所述层吸附密度与所述桥接吸附密度变化的函数关系Cd(Γl，Γb)，建立下式(13-3)表示的所述流体中的聚合物的质量浓度的平衡方程，
其中，
为所述流体中的聚合物的质量浓度
以及Γl与Γb分别为所述层吸附密度与所述桥接吸附密度。
步骤S13303，根据所述流体中的聚合物的分子链的数目分布方程，建立所述吸附聚合物的分子链的占比分布方程。
对于步骤S13303，所述建立所述吸附聚合物的分子链的占比分布方程可包括：根据所述流体中的聚合物的分子链的数目分布方程及
建立下式(13-4)表示的所述吸附聚合物的分子链的占比分布方程，
其中，C0为所述流体中的聚合物的初始质量浓度；N0为所述流体中的聚合物的分子链的初始数目；Ppoly为吸附聚合物的分子链数占所述流体中的聚合物的分子链的初始数目的比例；以及
为所述流体的速度。
具体地，将等式(13-1)表示的所述流体中的聚合物的分子链的数目分布方程的左端第一项做泰勒一阶展开(忽略二阶项)，然后方程两边同时除以Δt，其中
当Δt→0时，上式(13-1)变为：
其中，
分别表示发生层吸附的吸附率与发生桥接吸附的吸附率(即单位时间内吸附的聚合物分子链数量)：
其中，ρs表示吸附聚合物的密度；
分别表示发生层吸附和桥接吸附的介质比表面积(可视为
)；φ0表示所述储层的孔隙度；Γl表示聚合物的层吸附密度；Γb表示聚合物的桥接吸附密度。将上两式代入(7)式得到下式：
然后，将上式(13-6)两端同时除以N0(即所述流体中的聚合物的分子链的初始数目)，得到公式(13-4)表示的所述吸附聚合物的分子链的占比分布方程。
步骤S13203，根据层吸附速率与层解吸附速率、桥接吸附速率与桥接解吸附速率及动能守恒定律，确定在所述聚合物的吸附过程未达到饱和状态时所述聚合物的层吸附密度与桥接吸附密度。
聚合物的吸附过程同时遵循动能守恒定律，即单位时间聚合物分子发生吸附与解吸附满足动态平衡，并最终达到吸附饱和状态。对于层吸附和桥接吸附的解、吸附动态平衡可表示为式(13-7)和式(13-8)。
对于步骤S13203中的部分内容，所述确定在所述聚合物的吸附过程未达到饱和状态时所述聚合物的层吸附密度可包括：根据所述层吸附速率与所述层解吸附速率及下式表示的动能守恒定律，确定所述聚合物的所述层吸附密度，
其中，
为所述层吸附速率；
为所述层解吸附速率；
为在所述聚合物的吸附过程达到饱和状态时所述聚合物的饱和层吸附密度；以及Γl为所述层吸附密度，
对于步骤S13203中的部分内容，所述确定在所述聚合物的吸附过程未达到饱和状态时所述聚合物的桥接吸附密度可包括：根据所述桥接吸附速率与所述桥接解吸附速率及下式表示的动能守恒定律，确定在所述聚合物的吸附过程未达到饱和状态时所述聚合物的桥接吸附密度，
其中，
为所述桥接吸附速率；
为所述桥接解吸附速率；
为在所述聚合物的吸附过程达到饱和状态时所述聚合物的饱和桥接吸附密度；Γb为所述桥接吸附密度；以及Ppoly为吸附聚合物的分子链数占所述流体中的聚合物的分子链的初始数目的比例。
上两式(13-7)-(13-8)左侧项
表示聚合物吸附过程未达到饱和状态时单位时间内的聚合物的层吸附密度，右侧
表示单位时间内聚合物的动态吸附密度，
表示单位时间内的动态解吸附密度，其中i＝l(层吸附)、b(桥接吸附)。
步骤S13204，根据所述流体中的聚合物的质量浓度的平衡方程、所述吸附聚合物的分子链的占比分布方程、所述层吸附密度及所述桥接吸附密度，确定所述聚合物损害储层的时空演化模拟方程。
将公式(13-2)代入公式(13-3)的右侧，得到下式(13-9)，
然后将上式(13-7)-(13-8)代入公式(13-4)、(13-9)，可得到所述聚合物损害储层的时空演化模拟方程。也就是说，所述聚合物损害储层的时空演化模拟方程相当于由公式(13-4)、(13-7)-(13-9)组成的方程组。
综上所述，本发明创造性地本发明创造性地根据所述流体的速度、所述流体中的聚合物的分子链的数目分布方程及所述聚合物的扩散系数，建立所述流体中的聚合物的质量浓度的平衡方程及所述聚合物中吸附在所述储层上的吸附聚合物的分子链的占比分布方程；根据层吸附速率与层解吸附速率、桥接吸附速率与桥接解吸附速率及动能守恒定律，确定在所述聚合物的吸附过程未达到饱和状态时所述聚合物的层吸附密度与桥接吸附密度；以及根据所述流体中的聚合物的质量浓度的平衡方程、所述吸附聚合物的分子链的占比分布方程、所述层吸附密度及所述桥接吸附密度，确定所述聚合物损害储层的时空演化模拟方程。由此，通过所确定的时空演化模拟方程可定量模拟由聚合物吸附引起的储层损害特征的四维时空演化过程，从而对未发生储层损害的井进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义，以及对已损害井优化设计解堵措施、提高或恢复油井产量和水井注水能力，以及提高油藏数值模拟精度都具有十分重大意义。
图13-4是本发明一实施例提供的确定储层损害程度的方法的流程图。如图13-4所示，所述方法可包括步骤S13401-S13402。
步骤S13401，基于所述的聚合物损害储层的建模方法建立的时空演化模拟方程，确定所述流体中的聚合物的质量浓度。
对于上述公式(13-7)所示的聚合物吸附损害储层的聚合物的运移方程，可参见上述实施例1求解沉积颗粒的体积浓度的过程来计算得到所述流体中的聚合物的质量浓度
通过上述方法计算得到所述流体中的聚合物的质量浓度
由于通过上述聚合物吸附损害储层的建模方法建立的时空演化模拟方程综合考虑了储层内发生聚合物吸附时多种物理化学因素对储层损害的影响，由此通过该步骤S13401求解得到的所述流体中的聚合物的质量浓度非常精确。
步骤S13402，基于所述流体中的聚合物的质量浓度，确定表征待诊断井的预设区域内的储层的损害程度的特征参数。
其中，所述特征参数为所述储层的渗透率。
对于步骤S13402，所述确定表征待诊断井的预设区域内的储层的损害程度的特征参数可包括：基于所述流体中的聚合物的质量浓度
及公式 (13-10)，确定所述储层的渗透率
其中，φ0为孔隙度的初始值；
为所述流体中的聚合物的初始质量浓度；C∞为所述聚合物的最大吸附质量浓度；mK为第二经验值；以及
为所述储层的渗透率的初始值。
其中，所述特征参数为所述储层的滤失系数。
对于步骤S13402，所述确定表征待诊断井的预设区域内的储层的损害程度的特征参数可包括：基于所述流体中的聚合物的质量浓度
及公式 (13-11)，确定所述储层的滤失系数
其中，φ0为孔隙度的初始值；
为所述流体中的聚合物的初始质量浓度；C∞为所述聚合物的最大吸附质量浓度；mk为第一经验值；以及
所述储层的滤失系数的初始值。
其中，所述特征参数为所述储层的表皮系数。
对于步骤S13402，所述确定表征待诊断井的预设区域内的储层的损害程度的特征参数包括：基于所述流体中的聚合物的质量浓度
及公式
确定所述储层的渗透率
以及基于所述储层的渗透率
及公式(13-12)，确定所述储层的表皮系数
其中，
为所述储层的渗透率的初始值；φ0为孔隙度的初始值；
为所述流体中的聚合物的初始质量浓度；mk为第一经验值；
rw为所述待诊断井的井筒半径；以及rsw为所述储层的损害半径。
通过该步骤S13402得到的特征参数(例如所述储层的渗透率
与表皮系数
)是时空演化4D定量模拟的结果(如图13-5所示)。更具体地，图13-6示出了由储层渗透率损害率(基于所述储层的渗透率
及公式
确定所述储层的渗透率损害率I(ri,t)，其中
为
的最大值)表征的在第365天聚合物吸附损害储层的半径的示意图(如箭头所指示的半径)，相关工作人员可通过该图13-6直观地确认储层被损害的程度。因此，可根据渗透率或表皮系数的演化特点进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义。
在通过上述各个实施例得到各个因素损害储层时储层的渗透率的情况下，可根据渗透率与表皮系数(或渗透率损害率)之间的关系式(例如公式 (3-8)或公式(5-14))得到各个因素损害储层时储层的表皮系数(或渗透率损害率)。然后，针对同一诊断井的各个阶段，对多个因素各自损害储层时储层的渗透率(或表皮系数，或渗透率损害率)进行加权求和，以确定多个因素同时损害储层时储层的渗透率(或表皮系数，或渗透率损害率)。进一步，还可确定各个因素对总渗透率(或总表皮系数，或总渗透率损害率) 的贡献比例。
综上所述，本发明创造性地通过所确定的时空演化模拟方程可确定所述流体中的聚合物的质量浓度，再基于所确定的所述聚合物的质量浓度，确定表征待诊断井的预设区域内的储层的损害程度的特征参数(例如所述储层的渗透率和/或表皮系数)，由此可定量模拟由聚合物吸附引起的储层损害特征的四维时空演化过程，从而对未发生储层损害的井进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义，以及对已损害井优化设计解堵措施、提高或恢复油井产量和水井注水能力，以及提高油藏数值模拟精度都具有十分重大意义。
现场验证与应用
使用现场提供的相关参数，计算待诊断井若干时间后的总表皮系数(如图14-1至图14-4所示)，并与实测表皮系数对比，来验证本文建立的时空演化4D定量模拟技术的正确性和准确性。因实测表皮系数是由真实表皮系数和拟表皮系数之和，因此需要对实测表皮系数进行分解，利用分解后得到的真实表皮系数与计算的表皮系数进行对比，判断4D定量模拟技术的正确性和准确性。
下面以绥中36-1油田G01s1井为例。该井相关参数和实测表皮系数如表3所示。
表3绥中36-1油田G01s1井基本参数及表皮系数实测结果
| 参数 | 取值 | 参数 | 取值 |
| 实际生产时长(d) | 40 | 储层原始渗透率(10<sup>-3</sup>μm<sup>2</sup>) | 1660 |
| 实际生产压差(MPa) | 5.42 | 储层温度(℃) | 65 |
| 井筒半径(m) | 0.09 | 储层深度(m) | 1800 |
| 井斜角(°) | 40 | 储层原始孔隙度(％) | 32 |
| 储层厚度(m) | 38.3 | 储层孔喉平均半径(μm) | 14 |
| 实际完井厚度(m) | 33.7 | 储层油水界面张力(mN/m) | 13.5m |
| 采液速率(m<sup>3</sup>/d) | 45m<sup>3</sup>/d | 出砂砂粒密度(kg/m<sup>3</sup>) | 2.65 |
| 储层油相粘度(mPa.s) | 120 | 储层泥质含量(％) | 15.3 |
| 储层压力(MPa) | 13.2 | 岩石动态泊松比 | 0.13 |
| 上覆岩层压力(MPa) | 1800m | 泡点压力MPa | 11 |
| 岩石体积密度(kg/m<sup>3</sup>) | 2.5 | 原油密度g/cm<sup>3</sup> | 0.96 |
| 岩石杨氏模量(MPa) | 500 | 体系系数 | 1.1 |
| 岩石(静态)泊松比 | 0.13 | 胶质沥青质含量％ | 40.9 |
| 储层岩石压缩系数(MPa<sup>-1</sup>) | 5e-3 | 表皮系数 | 36.5 |
将表3中各参数输入加载有上述方法的软件，得到“出砂、应力敏感性、乳化、储层内部微粒运移、润湿反转、有机垢”6种损害类型导致损害带半径和损害程度(例如渗透率损害率)随时间的变化情况(分别如图7-4、图6-5、图9-4、图4-5、图8-4、图10-6)；以及40天后的总损害带半径(如图14-5 所示，其为各个损害因素对应的损害半径的最大值)和总损害程度(例如无因次渗透率，其等于1-渗透率损害率)随时间分布规律(如图14-6所示)；各类型因素造成的损害在总损害中所占比例随时间演化规律如图14-7所示。
从模拟计算结果可知，第40天时的总表皮系数35.65、实测表皮系数36.5、井斜和部分射孔导致的拟表皮系数0.49，相对误差仅0.99％，吻合度极高。从损害因素分析，乳化、储层内部微粒运移、有机垢、出砂导致的储层损害为主要损害因素，而润湿反转和应力敏感性损害很小。当然，根据需要还可模拟计算不同深度点、不同方向等4D空间任一储层损害参数分布和演化情况。
此外，本文建立的时空演化4D定量模拟技术在中海油等我国不同类型油田进行了广泛验证与推广应用，平均准确率达95％以上，充分证明本模拟技术的可靠性、准确性和实用性。部分井的模拟结果对比如表4所示。
表4部分井表皮系数模拟计算值与实测值对比
因此上述各个实施例可得到以下几个方面的优势：
(1)通过现场验证表明，采用连续介质力学和概率过程相结合的方法，建立的13种常见储层损害类型导致储层损害程度的时空演化4D定量数值模拟技术可用于油气田勘探开发全过程储层损害原因和程度时空演化4D定量模拟，精确度和准确度高、实用性和可操作性强；同时，证明采用连续介质力学和概率过程相融合的方法进行储层损害定量模拟研究，具有更高的科学性、可靠性和可行性，将是未来的发展方向。
(2)本文建立的储层损害时空演化4D定量模拟技术不仅可实现各损害类型和总损害程度的时空演化模拟，而且可提供每种损害类型对总损害的敏感程度，定量提供每种损害类型在总损害程度中的占比，为精准控制与解除储层损害、恢复油井产量和水井注水能力提供了核心技术。
(3)建立的储层损害时空演化4D定量模拟技术不仅可用于储层损害定量模拟也可用于储层损害定量预测。对已发生储层损害的井，利用历史参数实现储层损害定量模拟和时空演化，为解堵措施优化设计和提高油藏数值模拟精度具有重要意义；对未发生储层损害的井，利用物性参数和将实施的工程参数，可进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义。
(4)储层损害原因不仅具有多样性和动态变化性，而且各损害因素之间相互联系、相互制约，从而影响油水井总损害程度和时空演化规律，后期还应加强各损害类型相互影响和制约的研究工作，以进一步完善和提高储层损害定量模拟准确性。
综上所述，本发明创造性地基于多个因素损害储层的时空演化模拟方程，确定表征所述多个因素各自损害所述储层的特征参数；基于所述表征所述多个因素各自损害所述储层的特征参数，确定表征所述储层的损害程度的有效特征参数。由此，通过各个相关因素损害储层的时空演化模拟方程可定量模拟由其各自引起的储层损害的特征参数(例如渗透率)及多个相关因素对所述储层损害的总特征参数(例如总渗透率或有效渗透率)，从而对未发生储层损害的井进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义，以及对已损害井优化设计解堵措施、提高或恢复油井产量和水井注水能力，以及提高油藏数值模拟精度都具有十分重大意义。
相应地，本发明一实施例还提供一种确定储层损害程度的***，所述***包括：第一参数确定装置(未示出)，用于基于多个因素损害储层的时空演化模拟方程，确定表征所述多个因素各自损害所述储层的特征参数，其中，所述储层位于待诊断井的预设区域内；以及第二参数确定装置(未示出)，用于基于所述表征所述多个因素各自损害所述储层的特征参数，确定表征所述储层的损害程度的有效特征参数。
所述确定储层损害程度的***与上述确定储层损害程度的方法相对于现有技术所具有的优势相同，在此不再赘述。
相应地，本发明一实施例还提供一种机器可读存储介质，该机器可读存储介质上存储有指令，该指令用于使得机器执行所述的确定储层损害程度的方法。
所述机器可读存储介质包括但不限于相变内存(相变随机存取存储器的简称，Phase Change Random Access Memory，PRAM，亦称为 RCM/PCRAM)、静态随机存取存储器(SRAM)、动态随机存取存储器 (DRAM)、其他类型的随机存取存储器(RAM)、只读存储器(ROM)、电可擦除可编程只读存储器(EEPROM)、快闪记忆体(Flash Memory)或其他内存技术、只读光盘只读存储器(CD-ROM)、数字多功能光盘(DVD)或其他光学存储、磁盒式磁带，磁带磁盘存储或其他磁性存储设备等各种可以存储程序代码的介质。
上述各个实施例中的各个步骤均可通过计算机来执行，且某些步骤(例如步骤S1101-S1104等)所涉及的各种物理化学量的处理过程实现了对各个因素损害储层的时空演化场的模拟，以及某些步骤(例如步骤S1201-S1202) 所涉及的各种物理化学量的处理过程实现了对各个因素损害储层的时空演化的预测。上述各个实施例中的各个步骤均可由处理器执行完成。
以上结合附图详细描述了本发明的优选实施方式，但是，本发明并不限于上述实施方式中的具体细节，在本发明的技术构思范围内，可以对本发明的技术方案进行多种简单变型，这些简单变型均属于本发明的保护范围。
另外需要说明的是，在上述具体实施方式中所描述的各个具体技术特征，在不矛盾的情况下，可以通过任何合适的方式进行组合。为了避免不必要的重复，本发明对各种可能的组合方式不再另行说明。
此外，本发明的各种不同的实施方式之间也可以进行任意组合，只要其不违背本发明的思想，其同样应当视为本发明所公开的内容。
Claims (21)
1.一种确定储层损害程度的方法，其特征在于，所述方法包括：
基于多个因素损害储层的时空演化模拟方程，确定表征所述多个因素各自损害所述储层的特征参数，其中，所述储层位于待诊断井的预设区域内；以及
基于所述表征所述多个因素各自损害所述储层的特征参数，确定表征所述储层的损害程度的有效特征参数。
2.根据权利要求1所述的确定储层损害程度的方法，其特征在于，在所述待诊断井为注水井、注聚井或采油井并且处于钻井阶段的情况下，所述多个因素包括以下各项中的至少两者：外来固相颗粒、黏土膨胀、储层内部微粒运移、无机沉淀与水锁效应。
3.根据权利要求1所述的确定储层损害程度的方法，其特征在于，在所述待诊断井为采油井并且处于采油阶段的情况下，所述多个因素包括以下各项中的至少两者：储层内部微粒运移、出砂、乳化、贾敏效应、应力敏感性、润湿反转与有机垢。
4.根据权利要求1所述的确定储层损害程度的方法，其特征在于，在所述待诊断井为注水井并且处于注水阶段的情况下，所述多个因素包括以下各项中的至少两者：黏土膨胀、细菌、水锁效应、外来固相颗粒、微粒运移与无机沉淀。
5.根据权利要求1所述的确定储层损害程度的方法，其特征在于，在所述待诊断井为注聚井并且处于注聚阶段的情况下，所述多个因素包括以下各项中的至少两者：聚合物、黏土膨胀、外来固相颗粒、微粒运移与无机沉淀。
6.根据权利要求2、4-5中任一项所述的确定储层损害程度的方法，其特征在于，所述方法还包括：通过以下建模过程确定所述外来固相颗粒损害储层的时空演化模拟方程：
确定所述储层中的含有流动颗粒的流体的速度；
基于所述流体的对流参数与扩散参数，建立所述流体与所述储层中的岩石上的沉积颗粒之间的质量平衡方程；
基于所述流体的对流参数与扩散参数，建立所述沉积颗粒的体积浓度与所述流体的体积浓度之间的连接条件方程；以及
根据所述流动颗粒的质量分数与所述流动颗粒的体积浓度之间的关系、所述流体的速度、所述质量平衡方程及所述连接条件方程，确定所述颗粒损害储层的时空演化模拟方程，其中所述时空演化模拟方程用于模拟由外来固相颗粒引起的储层损害特征的四维时空演化过程。
7.根据权利要求2-5中任一项所述的确定储层损害程度的方法，其特征在于，所述方法还包括：通过以下建模过程确定所述储层内部微粒运移损害储层的时空演化模拟方程：
确定所述储层中的流体的速度；
基于所述流体的对流参数与扩散参数及所述流体内的运移微粒的质量变化率，建立所述流体与所述储层中的岩石上的沉积微粒之间的质量平衡方程，其中所述运移微粒的质量变化率与所述流体的速度之间具有关联关系；
基于所述流体的对流参数与扩散参数，建立所述沉积微粒的体积浓度与所述流体的体积浓度之间的连接条件方程；以及
根据所述运移微粒的质量分数与所述运移微粒的体积浓度之间的关系、所述流体的速度、所述质量平衡方程及所述连接条件方程，确定所述储层内部微粒损害储层的时空演化模拟方程，其中所述时空演化模拟方程用于模拟由微粒引起的储层损害特征的四维时空演化过程。
8.根据权利要求2或4所述的确定储层损害程度的方法，其特征在于，所述方法还包括：通过以下建模过程确定所述黏土膨胀损害储层的时空演化模拟方程：
确定所述储层中的流体的达西表观速度；
根据所述流体的达西表观速度与所述流体中的水分子的扩散系数，建立所述流体中的水分子的质量平衡方程；
根据菲克扩散定律，建立所述流体中的水分子向储层中的岩石内部扩散的扩散方程；以及
根据所述扩散方程及所述质量平衡方程，确定黏土膨胀损害储层的时空演化模拟方程，其中所述时空演化模拟方程用于模拟由黏土膨胀引起的储层损害特征的四维时空演化过程，以及所述黏土为所述岩石的组成成分。
9.根据权利要求2或4中所述的确定储层损害程度的方法，其特征在于，所述方法还包括：通过以下建模过程确定所述水锁效应损害储层的时空演化模拟方程：
确定所述储层中的流体的达西表观速度；
根据所述流体的达西表观速度与所述流体中的水分子的扩散系数，建立所述储层的水相运动方程；
根据所述储层的孔隙的孔径分布特征及所述储层的预设渗透率模型，建立所述储层的渗透率的分布方程；以及
根据所述渗透率的分布方程及所述水相运动方程，确定水锁效应损害储层的时空演化模拟方程，其中所述时空演化模拟方程用于模拟由水锁效应引起的储层损害特征的四维时空演化过程。
10.根据权利要求2、4-5中任一项所述的确定储层损害程度的方法，其特征在于，所述方法还包括：通过以下建模过程确定所述无机沉淀损害储层的时空演化模拟方程：
确定所述储层中的流体的达西表观速度；
确定与外来流体中的多个离子中的各个离子相对应的离子浓度损失，其中所述离子浓度损失是由所述各个离子与所述储层中的流体中的相应离子发生沉淀反应引起的；
根据所述流体的达西表观速度、与所述各个离子相对应的离子浓度损失、所述各个离子的扩散系数，建立所述各个离子的运移方程；以及
根据所述各个离子的运移方程及所述各个离子所产生的沉淀物的反应系数，确定所述无机沉淀损害储层的时空演化模拟方程，其中所述时空演化模拟方程用于模拟由所述多个离子产生的相应沉淀引起的储层损害特征的四维时空演化过程。
11.根据权利要求3所述的确定储层损害程度的方法，其特征在于，所述方法还包括：通过以下建模过程确定所述出砂损害储层的时空演化模拟方程：
确定所述储层中的流体的速度；
基于所述流体的对流参数与扩散参数及所述流体内的砂粒的质量变化率，建立所述流体与所述储层中的岩石上的沉积砂粒之间的质量平衡方程，其中所述砂粒的质量变化率与所述储层的原油产量之间具有关联关系；
基于所述流体的对流参数与扩散参数，建立所述沉积砂粒的体积浓度与所述流体的体积浓度之间的连接条件方程；以及
根据所述砂粒的质量分数与所述砂粒的体积浓度之间的关系、所述流体的速度、所述质量平衡方程及所述连接条件方程，确定所述出砂损害储层的时空演化模拟方程。
12.根据权利要求3所述的确定储层损害程度的方法，其特征在于，所述方法还包括：通过以下建模过程确定所述乳化损害储层的时空演化模拟方程：
确定所述储层中的流体的达西表观速度；
根据所述储层的温度场与所述储层中的油相的粘度与温度之间的函数关系，确定所述油相的粘度；
根据所述流体的达西表观速度、所述油相的粘度及所述流体的乳化条件，确定由所述流体乳化形成的乳化液滴的半径；以及
根据所述储层的孔隙的孔径分布函数及所述乳化液滴的半径，确定乳化堵塞损害储层的时空演化模拟方程，其中所述时空演化模拟方程用于模拟由乳化堵塞引起的储层损害特征的四维时空演化过程。
13.根据权利要求3所述的确定储层损害程度的方法，其特征在于，所述方法还包括：通过以下建模过程确定所述贾敏损害储层的时空演化模拟方程：
确定所述储层中的流体的达西表观速度，其中所述储层的渗透率低于预设渗透率；
根据所述流体的达西表观速度与所述流体中的水分子的扩散系数，建立所述储层的水相运动方程；
根据所述储层的孔隙的孔径分布特征及所述储层的预设渗透率模型，建立所述储层的渗透率的分布方程；以及
根据所述渗透率的分布方程及所述水相运动方程，确定贾敏效应损害储层的时空演化模拟方程，其中所述时空演化模拟方程用于模拟由贾敏效应引起的储层损害特征的四维时空演化过程。
14.根据权利要求3所述的确定储层损害程度的方法，其特征在于，所述方法还包括：通过以下建模过程确定所述应力敏感性损害储层的时空演化模拟方程：
确定所述储层受到的有效应力；
根据所述储层的孔隙的孔径分布特征、所述有效应力下的所述储层的各个毛管束的直径与长度、及流体流量公式，确定所述储层中的流体的流量，其中，所述毛管束由多个固体基质及所述多个固体基质之间的孔隙组成；以及
根据所述储层的渗透率模型及所述储层中的流体的流量，确定应力敏感性损害储层的时空演化模拟方程，其中，所述时空演化模拟方程用于模拟由应力敏感性引起的储层损害特征的四维时空演化过程。
15.根据权利要求3所述的确定储层损害程度的方法，其特征在于，所述方法还包括：通过以下建模过程确定所述有机垢损害储层的时空演化模拟方程：
确定所述储层的压力；
根据所述储层的泡点压力、所述储层产出的原油在所述泡点压力下的摩尔体积、所述原油的溶解度参数、所述有机垢在所述原油中的溶解度参数及所述有机垢的摩尔体积，确定所述有机垢在所述原油中的最大溶解量随所述储层的压力变化的第一关系式；
根据所述有机垢中的有机垢颗粒的分布函数及所述原油的摩尔数，确定所述原油中的有机垢颗粒的摩尔数随所述储层的压力与所述有机垢在所述原油中的最大溶解量两者变化的第二关系式，其中，所述分布函数为粒径大于预设粒径的有机垢颗粒的摩尔数占有机垢颗粒的总摩尔数的比例函数；以
根据所述第二关系式、所述第一关系式及所述储层的压力，确定有机垢损害储层的时空演化模拟方程，其中所述时空演化模拟方程用于模拟由有机垢引起的储层损害特征的四维时空演化过程。
16.根据权利要求3所述的确定储层损害程度的方法，其特征在于，所述方法还包括：通过以下建模过程确定所述润湿反转损害储层的时空演化模拟方程：
根据所述储层的压力分布方程，确定所述储层中的水相的压力分布场与毛细管的压力分布场之间的关系，其中，所述毛细管由所述水相与所述储层中的油相的接触界面的润湿反转而形成；
根据所述水相的压力分布场与所述毛细管的压力分布场之间的关系及所述毛细管的受力平衡条件，确定所述油相的压力分布场；
根据所述油相的压力分布场及达西公式，确定所述油相的速度分布场；以及
根据所述油相的对流扩散定律、所述速度分布场及所述油相的弥散系数，确定润湿反转损害储层的时空演化模拟方程，其中所述时空演化模拟方程用于模拟由润湿反转引起的储层损害特征的四维时空演化过程。
17.根据权利要求4所述的确定储层损害程度的方法，其特征在于，所述方法还包括：通过以下建模过程确定所述细菌损害储层的时空演化模拟方程：
根据所述储层的温度分布场及所述储层中的流体中的营养物的实际浓度，确定所述细菌的生长速率；
根据与所述流体中的细菌的表观浓度及所述细菌在所述岩石表面上的总量两者相关联的所述细菌在所述储层中的岩石表面上的附着量、所述细菌的生长速率与衰亡速率，确定所述细菌在所述岩石表面上的总量；
根据所述流体的达西速度、所述细菌的弥散系数、所述细菌的生长速率与衰亡速率及所述细菌在所述储层中的岩石表面上的附着量，建立所述流体中的细菌的表观浓度分布方程；
根据所述流体的达西速度、所述营养物的弥散系数、所述细菌在所述岩石表面上的总量及所述细菌的表观浓度，建立所述营养物的表观浓度分布方程；以及
根据所述营养物的表观浓度分布方程与所述流体中的细菌的表观浓度分布方程，确定细菌损害储层的时空演化模拟方程，其中所述时空演化模拟方程用于模拟由细菌引起的储层损害特征的四维时空演化过程。
18.根据权利要求5所述的确定储层损害程度的方法，其特征在于，所述方法还包括：通过以下建模过程确定所述聚合物损害储层的时空演化模拟方程：
确定所述储层中的流体的速度；
根据所述流体的速度、所述流体中的聚合物的分子链的数目分布方程及所述聚合物的扩散系数，建立所述流体中的聚合物的质量浓度的平衡方程及所述聚合物中吸附在所述储层上的吸附聚合物的分子链的占比分布方程，其中所述吸附聚合物的分子链的占比为所述吸附聚合物的分子链的数目占所述流体中的聚合物的分子链的初始数目的比例；
根据层吸附速率与层解吸附速率、桥接吸附速率与桥接解吸附速率及动能守恒定律，确定在所述聚合物的吸附过程未达到饱和状态时所述聚合物的层吸附密度与桥接吸附密度；以及
根据所述流体中的聚合物的质量浓度的平衡方程、所述吸附聚合物的分子链的占比分布方程、所述层吸附密度及所述桥接吸附密度，确定所述聚合物损害储层的时空演化模拟方程。
19.根据权利要求1所述的确定储层损害程度的方法，其特征在于，所述确定表征所述储层的损害程度的有效特征参数包括：
基于所述表征所述多个因素中的第i个因素损害所述储层的特征参数
及下式，确定表征所述储层的损害程度的有效特征参数
其中，Li为
的权重；以及n为所述多个因素的个数。
20.一种确定储层损害程度的***，其特征在于，所述***包括：
第一参数确定装置，用于基于多个因素损害储层的时空演化模拟方程，确定表征所述多个因素各自损害所述储层的特征参数，其中，所述储层位于待诊断井的预设区域内；以及
第二参数确定装置，用于基于所述表征所述多个因素各自损害所述储层的特征参数，确定表征所述储层的损害程度的有效特征参数。
21.一种机器可读存储介质，其特征在于，所述机器可读存储介质上存储有指令，该指令用于使得机器执行上述权利要求1-19中任一项所述的确定储层损害程度的方法。
Applications Claiming Priority (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010873143 | 2020-08-26 | ||
| CN2020108731439 | 2020-08-26 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113919126A true CN113919126A (zh) | 2022-01-11 |
| CN113919126B CN113919126B (zh) | 2022-10-04 |
Family
ID=79233262
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110989424.5A Active CN113919126B (zh) | 2020-08-26 | 2021-08-26 | 油气层损害类型和程度时空演化4d定量与智能诊断方法及其*** |
Country Status (3)
| Country | Link |
|---|---|
| US (1) | US11686190B2 (zh) |
| CN (1) | CN113919126B (zh) |
| WO (1) | WO2022042663A1 (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116307280A (zh) * | 2023-05-18 | 2023-06-23 | 成都理工大学 | 一种酸性气藏气井储层硫堵损伤量化评估方法 |
Families Citing this family (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2019045682A1 (en) * | 2017-08-29 | 2019-03-07 | Intel Corporation | SELECTOR DEVICES |
| CN112983368B (zh) * | 2021-03-25 | 2022-10-11 | 中国石油大学(华东) | 一种注采优化协同化学驱实现原油均衡驱替的方法 |
| CN115898375B (zh) * | 2022-12-20 | 2024-06-18 | 西南石油大学 | 一种模拟裂缝流固耦合变形的颗粒运移可视化实验装置及方法 |
| CN116127251B (zh) * | 2023-04-17 | 2023-09-29 | 西南石油大学 | 一种注二氧化碳的有水气藏的地层压力监测方法 |
| CN117130068A (zh) * | 2023-10-26 | 2023-11-28 | 东北石油大学三亚海洋油气研究院 | 断-砂配置侧向分流运移油气的判别方法及装置 |
| CN118095576B (zh) * | 2024-04-26 | 2024-07-09 | 中国石油大学(华东) | 一种疏松砂岩油藏储层微粒运移预测方法及评价方法 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103806904A (zh) * | 2012-11-12 | 2014-05-21 | 中国石油天然气集团公司 | 砂岩储层非均匀损害深度的计算方法 |
| CN105181927A (zh) * | 2015-08-05 | 2015-12-23 | 河南能源化工集团研究院有限公司 | 多场耦合低渗煤层水力压裂模拟试验方法 |
| CN107526892A (zh) * | 2017-08-30 | 2017-12-29 | 广州海洋地质调查局 | 一种海洋天然气水合物试采储层的稳定性评估方法 |
| CN108804779A (zh) * | 2018-05-24 | 2018-11-13 | 西安石油大学 | 一种基于数字岩心模型的稠油油藏储层解堵过程模拟方法 |
| CN109946436A (zh) * | 2019-04-15 | 2019-06-28 | 西南石油大学 | 一种兼顾基块与裂缝的裂缝性致密气层工作液损害评价方法 |
Family Cites Families (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108240955A (zh) * | 2016-12-27 | 2018-07-03 | 中国石油天然气股份有限公司 | 水泥浆对储层损害的评价***和方法 |
| CN109594982A (zh) * | 2018-12-19 | 2019-04-09 | 中国科学院广州能源研究所 | 一种含水合物储层损害的评价装置及评价方法 |
| CN110376112A (zh) * | 2019-07-29 | 2019-10-25 | 中国地质大学（北京） | 一种评价致密气储层水锁伤害程度的方法 |
-
2021
- 2021-08-26 CN CN202110989424.5A patent/CN113919126B/zh active Active
- 2021-08-26 WO PCT/CN2021/114851 patent/WO2022042663A1/zh active Application Filing
-
2022
- 2022-07-13 US US17/864,224 patent/US11686190B2/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103806904A (zh) * | 2012-11-12 | 2014-05-21 | 中国石油天然气集团公司 | 砂岩储层非均匀损害深度的计算方法 |
| CN105181927A (zh) * | 2015-08-05 | 2015-12-23 | 河南能源化工集团研究院有限公司 | 多场耦合低渗煤层水力压裂模拟试验方法 |
| CN107526892A (zh) * | 2017-08-30 | 2017-12-29 | 广州海洋地质调查局 | 一种海洋天然气水合物试采储层的稳定性评估方法 |
| CN108804779A (zh) * | 2018-05-24 | 2018-11-13 | 西安石油大学 | 一种基于数字岩心模型的稠油油藏储层解堵过程模拟方法 |
| CN109946436A (zh) * | 2019-04-15 | 2019-06-28 | 西南石油大学 | 一种兼顾基块与裂缝的裂缝性致密气层工作液损害评价方法 |
Non-Patent Citations (5)
| Title |
|---|
| LEILI等: "A review of the current status of induced seismicity monitoring for hydraulic fracturing in unconventional tight oil and gas reservoirs", 《FUEL》 * |
| 宁伏龙等: "水合物开采储层出砂数值模拟研究进展", 《石油科学通报》 * |
| 易飞等: "基于水动力不稳定性的油水乳化储层伤害数值模拟", 《钻井液与完井液》 * |
| 王松等: "预测钻井完井液对油气层损害的数学模型", 《钻采工艺》 * |
| 蒋官澄等: "低渗特低渗储层水锁损害定量预测方法", 《石油钻探技术》 * |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116307280A (zh) * | 2023-05-18 | 2023-06-23 | 成都理工大学 | 一种酸性气藏气井储层硫堵损伤量化评估方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| US11686190B2 (en) | 2023-06-27 |
| CN113919126B (zh) | 2022-10-04 |
| US20230026538A1 (en) | 2023-01-26 |
| WO2022042663A1 (zh) | 2022-03-03 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN113919126B (zh) | 油气层损害类型和程度时空演化4d定量与智能诊断方法及其*** | |
| US9898559B2 (en) | Method for predicting changes in properties of a formation in a near-wellbore area exposed to a drilling mud | |
| Izadi et al. | A new approach in permeability and hydraulic-flow-unit determination | |
| CN109815516A (zh) | 对页岩气井产能进行预测的方法及装置 | |
| Mahmoudi et al. | An experimental evaluation of pore plugging and permeability reduction near SAGD sand control liners | |
| GB2515411A (en) | Well simulation | |
| CN110263362B (zh) | 基于孔隙密度流的岩土体离散元流固耦合数值模拟方法 | |
| CN105629327B (zh) | 一种针对弱胶结、深层致密砂岩储层成岩相定量表征方法 | |
| Yang et al. | Influence of reservoirs/interlayers thickness on hydraulic fracture propagation laws in low-permeability layered rocks | |
| Hou et al. | Numerical simulation and evaluation of the fracturing and tight gas production with a new dimensionless fracture conductivity (f cd) model | |
| Mojarad et al. | Coupled numerical simulation of reservoir flow with formation plugging | |
| Ma et al. | Experimental investigation of the plugging mechanisms of non-consolidated prepacked gravel screens | |
| Liu et al. | Formation damage by fines migration including effects of filter cake, pore compressibility, and non-Darcy flow—a modeling approach to scaling from core to field | |
| CN113655546B (zh) | 出砂损害油气层建模方法、损害程度时空演化4d定量与智能诊断方法及其*** | |
| Mitchel et al. | xzsEffect of Drilling Fluid (Water-Based vs Oil-Based) on Phase Trap Damage in Tight Sand Gas Reservoirs | |
| Yao et al. | Estimation of shale intrinsic permeability with process-based pore network modeling approach | |
| CN113705123B (zh) | 外来颗粒损害油气层建模方法、损害程度时空演化4d定量与智能诊断方法及其*** | |
| CN115828440A (zh) | 页岩油水平井水力裂缝延伸模型的构建方法 | |
| Haugland | Spontaneous Imbibition in Sand-Viscosity Effects on Oil Recovery and Flow Using Polymer and Glycerol | |
| CN113919242B (zh) | 微粒运移损害油气层建模方法、损害程度时空演化4d定量与智能诊断方法及其*** | |
| Zhou et al. | Quantitative Prediction of Sand Production in Weakly Consolidated Sandstone Through Numerical Simulation | |
| Won et al. | Investigation of mud-filtrate invasion using computational fluid dynamics | |
| Wennberg | Particle Retention in Porous Media: Applications to Water Injectivity Decline | |
| Dai et al. | Recognition of the Lower Limit of Reservoir Physical Properties | |
| Chan et al. | From Wellbore Instability and Grain Mixing to Injectivity Reduction |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| CB03 | Change of inventor or designer information | ||
| CB03 | Change of inventor or designer information |
Inventor after: Jiang Guancheng Inventor after: Wang Yong Inventor after: Fu Daqi Inventor after: Geng Tie Inventor after: Ran Qihua Inventor after: Liu Xiaobo Inventor after: Cheng Rongchao Inventor after: Wang Zenglin Inventor after: Chen Gang Inventor after: Zhu Honghao Inventor after: Li Yizheng Inventor after: He Kenbo Inventor after: Zhao Li Inventor after: Cui Kaixiao Inventor after: Dong Tengfei Inventor after: Peng Chunyao Inventor after: Yang Lili Inventor after: Luo Xiaohu Inventor after: Luo Xuwu Inventor after: Liang Xing Inventor after: Tan Bin Inventor before: Jiang Guancheng Inventor before: Wang Yong Inventor before: Fu Daqi Inventor before: Geng Tie Inventor before: Ran Qifa Inventor before: Liu Xiaobo Inventor before: Cheng Rongchao Inventor before: Wang Zenglin Inventor before: Chen Gang Inventor before: Zhu Honghao Inventor before: Li Yizheng Inventor before: He Kenbo Inventor before: Zhao Li Inventor before: Cui Kaixiao Inventor before: Dong Tengfei Inventor before: Peng Chunyao Inventor before: Yang Lili Inventor before: Luo Xiaohu Inventor before: Luo Xuwu Inventor before: Liang Xing Inventor before: Tan Bin |
|
| GR01 | Patent grant | ||
| GR01 | Patent grant |