Shale oil reservoir multi-scale mechanical property prediction method

Abstract

The invention discloses a shale oil reservoir multi-scale mechanical property prediction method. Comprising the following steps: S1, carrying out multi-scale structure characterization on a shale oil reservoir; s2, carrying out a microscopic mechanical test; s3, carrying out primary homogenization analysis to obtain an equivalent elastic modulus of the three-dimensional grain model; s4, carrying out a macroscopic physical experiment; s5, determining a mean value and a variance of lamina distribution; and S6, carrying out secondary homogenization analysis, and determining the equivalent elastic modulus of the lamina distribution model. The method has the advantages that the equivalent elastic modulus of the strata distribution model can be predicted in a cross-scale mode, and a more accurate mechanical basis is provided for shale oil reservoir development and engineering practice.

Description

Technical Field

[0001] This invention relates to the field of rock mechanics and numerical simulation, and in particular to a method for predicting the multi-scale mechanical properties of shale oil reservoirs. Background Technology

[0002] In rock mechanics research and engineering applications, accurately predicting the mechanical properties of rock samples is crucial for ensuring engineering safety and optimizing design. Traditional methods for predicting rock mechanics parameters, such as those based on classical homogeneous isotropic elastic wave theory, suffer from low accuracy due to the significant differences in complex rock mechanics parameters and the intricate nature of rock physical response laws. Therefore, a multi-scale mechanical property prediction method for shale oil reservoirs is proposed. Summary of the Invention

[0003] The purpose of this invention is to provide a multi-scale mechanical property prediction method for shale oil reservoirs, which can predict the equivalent elastic modulus of the lamination distribution model across scales, providing a more accurate mechanical basis for the development and engineering practice of shale oil reservoirs.

[0004] To achieve the above objectives, the present invention adopts the following technical solution, comprising the following steps:

[0005] S1. Multi-scale structural characterization of shale oil reservoirs;

[0006] High-resolution electron microscopy was used to acquire microscopic images of grain morphology at the micro-nano scale of shale oil reservoir samples, and the grain morphology and arrangement patterns were observed and statistically analyzed. X-ray CT scanning was used to acquire three-dimensional structural images of full-diameter shale oil reservoir cores at the millimeter scale, and image segmentation algorithms were used to extract lamellar thickness, lamellar lateral continuity length, and lamellar distribution patterns. Lamellar characteristics were obtained from the grain microscopic morphology images and the three-dimensional core structural images; these lamellar characteristics include lamellar thickness, lamellar lateral continuity length, lamellar distribution patterns, grain morphology, and arrangement patterns.

[0007] S2. Conduct micromechanical experiments;

[0008] Based on the laminar features obtained in S1, typical laminar types were determined; full-diameter core samples from shale oil reservoirs were analyzed to identify typical regions; rock fragments were obtained from typical regions, and mineral dissociation analysis was used to separate and identify mineral components for each rock fragment to determine their types; microscopic morphology images of the grains obtained in S1 for this typical region were selected, and combined with the identified mineral types, a two-dimensional image that can distinguish mineral types was generated; a nanoindenter was used to test a single mineral region to obtain the microscopic elastic modulus of different mineral grains;

[0009] S3. Conduct a homogenization analysis to obtain the equivalent elastic modulus of the three-dimensional grain model;

[0010] Based on the two-dimensional image obtained by mineral dissociation analysis, pixels are converted into mesh units to construct a three-dimensional grain model; then, based on the value of the micro-elastic modulus of the mineral grains of each mineral in step S2, a micro-representative volume element model reflecting the periodic microstructure of the grain model is constructed, and periodic boundary conditions are applied by the finite element method to ensure the continuity of boundary displacement; the equivalent elastic modulus of the three-dimensional grain model is obtained.

[0011] S4. Conduct macroscopic physics experiments;

[0012] Shale oil reservoir cores corresponding to multi-scale structural characterization were selected, and cylindrical samples conforming to industry standards were prepared, covering different bedding angles of 0°, 30°, 45°, 60°, and 90°. A servo-controlled triaxial rock testing system was used, with confining pressure gradients set to 5MPa, 10MPa, 15MPa, and 20MPa, and axial strain rate of 0.001mm / s was applied until the sample failed. The failure mode was recorded, and the macroscopic elastic modulus and Poisson's ratio were calculated based on the linear segment of the stress-strain curve.

[0013] S5. Determine the mean μ of the striation distribution. X and variance

[0014] The three-dimensional core structure image from step S1 is analyzed to determine the relevant lengths of the lamellar layers. The core three-dimensional structure image is placed horizontally, with the average height of the lamellar layers as Lz and the horizontal extension lengths as Lx and Ly. The dip direction (DD) and dip angle (Dip) are further determined based on the lamellar direction. The mean value (μ) of the lamellar distribution is determined based on the analysis results of the core three-dimensional structure image. X and variance A three-dimensional random field modeling method is used to construct a striation distribution model;

[0015] S6. Conduct a secondary homogenization analysis to determine the equivalent elastic modulus of the lamellar distribution model;

[0016] Based on the equivalent elastic modulus of the laminar distribution model constructed in step S5 and the three-dimensional grain model corresponding to different laminar layers obtained in S3, a microscopic representative volume element model reflecting the periodic microstructure of the laminar distribution model is constructed, periodic boundary conditions are applied to ensure the continuity of boundary displacement, and the equivalent elastic modulus of the laminar distribution model is determined.

[0017] Preferably, step S1 includes:

[0018] S101. Select shale oil reservoir rock samples and prepare ultrathin sections with a thickness ≤5μm after argon ion polishing. Observe the samples using a field emission scanning electron microscope or transmission electron microscope at a magnification of 5000-50000 times to obtain images of the grain micromorphology. Analyze the grain micromorphology images and statistically record the following parameters: grain morphology includes the grain size distribution, aspect ratio, and angularity coefficient of different mineral grains; arrangement pattern includes grain distribution and aggregation degree.

[0019] S102. X-ray CT scanning was performed on the full-diameter core of the shale oil reservoir with a scanning resolution ≥50μm to obtain a three-dimensional structural image of the core. A deep learning-based image segmentation algorithm was used to identify the lamellar interface and extract and quantify the following features: geometric parameters include lamellar thickness, lamellar lateral continuous length, and lamellar distribution pattern; among which, lamellar thickness is divided into millimeter-level lamellar (<1cm), centimeter-level lamellar (1-10cm), and decimeter-level blocky (>10cm); lamellar distribution pattern includes the number of lamellars per unit length, i.e., lamellar density, and the distribution uniformity index in different directions.

[0020] Preferably, step S2 includes:

[0021] S201. Based on the described laminar characteristics, determine the typical laminar type; analyze the full-diameter core of the shale oil reservoir to determine the typical region; obtain rock cuttings from the typical region;

[0022] S202. The rock fragments are crushed to a particle size ≤0.1mm and automatically identified using a mineral dissociation analyzer. The elemental composition of each mineral is determined by scanning electron microscopy combined with energy dispersive spectroscopy, thereby identifying the mineral types. The microscopic morphology image of the grains obtained in S1 of the typical area is selected and combined with the mineral types identified above to generate a two-dimensional image that can distinguish the mineral types.

[0023] S203. A nanoindenter was used to test a single mineral region to obtain the micro-elastic modulus of different mineral grains. The micro-elastic modulus Er of the mineral grains was calculated as follows:

[0024] The unloading curve is an exponential function, as shown in formula (1):

[0025] P = B(hh) f ) m (1)

[0026] In equation (1), P represents the load applied during the test, in mN; h represents the indentation depth, in nm; h f The indentation depth at the end of the test is represented in nm; B and m are empirical constants obtained by fitting the test data.

[0027] The contact stiffness S is calculated as shown in formula (2):

[0028]

[0029] In equation (2), h max The maximum indentation depth of the indenter is expressed in nm; P represents the load applied during the test, expressed in mN; h represents the test indentation depth, expressed in nm. f The indentation depth at the end of the test is represented in nm; B and m are empirical constants obtained by fitting the test data.

[0030] Microscopic elastic modulus Er of mineral grains:

[0031]

[0032] In equation (3), β is generally a constant value, which is related to the pressure head selected in the experiment; A is the contact area, in nm. 2 S represents the contact stiffness.

[0033] S204. Perform statistical analysis on the test data for each type of lamination, calculate the average value and standard deviation of the elastic modulus of each mineral grain, and finally establish a correlation table of microscopic elastic modulus between lamination type, mineral composition and mineral grain.

[0034] Preferably, step S3 includes:

[0035] S301. The two-dimensional image with a resolution ≥1024×1024 pixels is preprocessed, and the boundaries of different mineral phases are identified by threshold segmentation method;

[0036] S302. Map each pixel to a square grid cell with a side length of 0.1-1μm to construct a three-dimensional grain model;

[0037] S303. Based on the value of the micro-elastic modulus of the mineral grains in step S2, construct a micro-representative volume element model that reflects the periodic microstructure of the grain model. Apply periodic boundary conditions through the finite element method to ensure the continuity of boundary displacement. The applied periodic boundary conditions satisfy equation (4).

[0038] u i (x+Y j )=u i (x) (4)

[0039] In equation (4), u i Let i = x, y, z, representing displacement components, corresponding to displacements in the x, y, and z directions respectively; x is the spatial coordinate, representing the position of any point within the RVE; y... j Let be the period length of RVE in the j direction, where j = x, y, z;

[0040] S304. Based on the assumption of asymptotic expansion at two scales, the displacement field is decomposed into macroscopic average displacement and microscopic perturbation displacement; based on two-scale expansion, the displacement field is expressed as:

[0041] u ξ (x)=u 0 (x,y)+ξu 1 (x,y)+…… (5)

[0042] In equation (5), y is the microscopic coordinate, y = x / ξ, ξ is the scale parameter, which is the ratio of the macroscopic feature size to the microscopic feature size, ξ << 1; ξ (x) represents the displacement field including scale effects; u 0 (x, y) represents the macroscopic average displacement field; u 1 (x, y) represents the microscopic perturbation displacement field;

[0043] And combined with linear elastic constitutive model:

[0044] σ ij =D ijkl ε kl (6)

[0045] In equation (6), σ ij For the stress tensor components; i,j=x,y,z, representing normal stress or shear stress in different directions; D ijkl ε is the fourth-order elastic stiffness tensor; kl For the strain tensor components, k,l=x,y,z, representing normal strain or shear strain in different directions;

[0046] And the equilibrium equations:

[0047] σ ij,j +f i =0 (7)

[0048] In equation (7), σ ij,j Let f be the divergence of the stress tensor, representing the rate of change of stress in space, and "," denote the partial derivative with respect to coordinates; i Let i be the volume force components, i = x, y, z;

[0049] Establish homogenized control equations and solve for the microstress σ under a given macroscopic strain using the finite element method. ij and strain ε ij Distribution, then by volume average:

[0050]

[0051] In equations (8) and (9), This is the macroscopic average stress tensor; ∫ is the macroscopic average strain tensor; V is the volume of the RVE;V σ ij dV, ∫ V ε ij dV represents the integrals of microscopic stress and strain within the volume of RVE, respectively;

[0052] According to the generalized Hooke's law, equation (10) is:

[0053]

[0054] In equation (12), C is the macroscopic average stress tensor; ijkl This is the equivalent elastic modulus of the three-dimensional grain model; For macroscopic average strain tensor;

[0055] The mechanical response of the RVE under the above boundary conditions was solved using the finite element method, and the results were obtained for each point within the RVE at each strain tensor ε. ij The stress tensor σ under (x) ij (x), and then the equivalent elastic modulus C of the three-dimensional grain model is derived. ijkl .

[0056] Preferably, in step S4, the macroscopic elastic modulus E is calculated as shown in formula (11):

[0057]

[0058] In equation (11), E is the macroscopic elastic modulus, with units of GPa; σ c (50) represents 50% of the uniaxial compressive strength of the specimen, in MPa; ε c (50) is σ c (50) corresponds to the axial strain;

[0059] The commonly used formula for calculating Poisson's ratio μ is shown in formula (12):

[0060]

[0061] In equation (12), μ is the Poisson's ratio of the specimen, and εc(50) and εt(50) are the axial strain and transverse strain corresponding to σc(50), respectively.

[0062] Preferably, step S5 includes:

[0063] S501. Based on the three-dimensional core structure image obtained in step S1, Gaussian filtering is used for noise reduction. The processed three-dimensional core structure image is placed horizontally, and the thickness of the laminae is counted along the direction perpendicular to the laminae extension, and the average laminae height Lz is calculated. Along the horizontal directions x and y, the continuous extension length of the laminae is counted, and the horizontal extension lengths Lx and Ly are obtained respectively. Combined with the core attitude data, the dip direction DD and dip angle Dip of the laminae are determined. The dip direction DD ranges from 0° to 360°, and the dip angle Dip ranges from 0° to 90°.

[0064] S502. Perform grayscale processing on the three-dimensional structure image of the core, and calculate the correlation length of the laminae in the x, y, and z directions using the autocorrelation function; for the x direction, the autocorrelation function is used.

[0065]

[0066] In equation (13), ρ(Δx) is the value of the autocorrelation function, which characterizes the strength of the correlation between two points spaced Δx apart in space; Δx is the distance between the two points in space along the x-direction, in mm; c is the lower limit of the autocorrelation function, which is a constant, 0≤c<1, representing the limit of the autocorrelation function when the distance between the two points Δx→∞; Lx is the correlation length in the x-direction. This is the exponentially decaying term, describing how the autocorrelation function changes with the distance Δx between the two points;

[0067] S503. Treat the thickness or spacing of the texture as a discrete random variable X, and its probability distribution law is as follows:

[0068] P{X=x t}=p t (t=1,2,...n) (14)

[0069] In equation (14), n is the sample size, p t xt represents the frequency; xt is the t-th possible value of the discrete random variable X, where t = 1, 2, ..., n.

[0070] The mean is calculated as follows:

[0071]

[0072] In equation (15), μ X (or E(X)) represents the mean of the random variable X; X is a discrete random variable; n is the total number of elements contained in the discrete random variable X within the region D; t is the summation index, used to iterate through the counts from 1 to n; xt is the t-th possible value of the discrete random variable X, t = 1, 2, ..., n; p t Let P{X=xt} be the probability that the random variable X takes the t-th value xt, and satisfy that the sum of all such values ​​is 1.

[0073] Calculate the variance;

[0074] Var(x) = E{[XE(X)]} 2} (16)

[0075] In equation (16), Var(x) represents the variance of the random variable X; X is a random variable; E(X) is the mean of the random variable X;

[0076] Similarly, when the random variable X is discrete, the variance formula can be simplified to:

[0077]

[0078] In equation (17), Let X be the variance of the discrete random variable X; X is a discrete random variable; n is the number of all possible values ​​of the discrete random variable X; t is the summation index, used to iterate through each possible value of X, from 1 to n; xt is the t-th possible value of the discrete random variable X (t = 1, 2, ..., n); pt is the probability that the discrete random variable X takes the t-th value xt, i.e., P{X = xt}; E(X) is the mean of the discrete random variable X.

[0079] S504. A three-dimensional random field modeling method is adopted to construct a striatal distribution model based on exponential or quadratic exponential autocorrelation functions.

[0080] S504a, the relevant length Lx obtained in step S502, and the mean μ obtained in S503. X and variance Using the input parameters, a random field is generated through Fast Fourier Transform;

[0081] S504b discretizes the random field into a standard grid model, adjusts the grid spatial orientation according to the inclination DD and dip angle Dip of the striation, and finally obtains the three-dimensional spatial distribution of striation parameters, realizing the quantitative characterization of striation structure.

[0082] Preferably, step S6 includes: Step S6 includes:

[0083] S601. Based on the laminar distribution model obtained by three-dimensional random field modeling in step S5, and combined with the equivalent elastic modulus of the three-dimensional grain model obtained in step S3, construct a microscopic representative volume element model that reflects the periodic microstructure of the laminar distribution model. Apply periodic boundary conditions through the finite element method to ensure the continuity of boundary displacement. The applied periodic boundary conditions satisfy the following conditions.

[0084] u′ i (x′+Y′)=u′ i (x′) (18)

[0085] In equation (18), u′ i Let i = x, y, z, corresponding to displacements in the x, y, and z directions, respectively; x′ is the spatial coordinate, representing the position of any point within the RVE; Y′... j Let be the period length of RVE in the j direction, where j = x, y, z;

[0086] S602. Based on the dual-scale asymptotic expansion assumption, the displacement field is decomposed into macroscopic average displacement and microscopic perturbation displacement; based on dual-scale expansion, the displacement field is expressed as:

[0087] u ξ′ (x′)=u 0 (x′,y′)+ξ′u 1 (x′,y′)+…… (19) In equation (19), y' is the microscopic coordinate, y' = x' / ξ', ξ' is the scale parameter, which is the ratio of the macroscopic feature size to the microscopic feature size, and ξ' << 1; ξ′ (x') represents the displacement field with scale effects; u 0′ (x',y′) represents the macroscopic average displacement field; u 1′ (x',y') represents the microscopic perturbation displacement field;

[0088] And combined with linear elastic constitutive model:

[0089] σ′ ij =D′ ijkl ε′ kl (20)

[0090] In equation (20), σ′ ij For the stress tensor components; i, j = x, y, z, representing normal or shear stresses in different directions; D′ ijkl For the fourth-order elastic stiffness tensor; ε' kl For the strain tensor components, k,l=x,y,z, representing normal strain or shear strain in different directions;

[0091] And the equilibrium equations:

[0092] σ' ij，j +f' i =0 (21)

[0093] In equation (21), σ' ij，j f' is the divergence of the stress tensor, representing the rate of change of stress in space; , "," denotes the partial derivative with respect to coordinates; i Let i be the volume force components, i = x, y, z;

[0094] Establish homogenized control equations and solve for the microstress σ' under a given macroscopic strain using the finite element method. ijand strain ε' ij Distribution, then by volume average:

[0095]

[0096] In equations (22) and (23), This is the macroscopic average stress tensor; V is the macroscopic average strain tensor; V' is the volume of the RVE; ∫ V′ σ′ ij dV, ∫ V′ ε′ ij dV represents the integrals of microscopic stress and strain within the volume of RVE, respectively;

[0097] According to the generalized Hooke's law, equation (12) is:

[0098]

[0099] In equation (24), C′ is the macroscopic average stress tensor. ijkl This is the equivalent elastic modulus of the lamellar distribution model; For macroscopic average strain tensor;

[0100] The mechanical response of the RVE under the above boundary conditions was solved using the finite element method, and the results were obtained for each point within the RVE at each strain tensor ε′. ij The stress tensor σ′ under (x) ij (x), from which the equivalent elastic modulus C' of the lamellar distribution model is derived. ijkl The macroscopic elastic modulus E obtained in step S4 is compared with the macroscopic elastic modulus E. If the error is within 20%, the elastic modulus of shale is considered to be accurately predicted.

[0101] The beneficial effects of this invention are: it can predict the equivalent elastic modulus of the laminar distribution model across scales, providing a more accurate mechanical basis for the development and engineering practice of shale oil reservoirs. Attached Figure Description

[0102] Figure 1 This is a flowchart of the multi-scale mechanical property prediction method for shale oil reservoirs according to the present invention;

[0103] Figure 2 This is the three-dimensional grain model in this invention;

[0104] Figure 3 This is a schematic diagram of the homogenization analysis results of the three-dimensional grain model in this invention;

[0105] Figure 4 This is the laminar distribution model in this invention;

[0106] Figure 5This is a schematic diagram of the homogenization analysis results of the laminar distribution model in this invention. Detailed Implementation

[0107] The invention will now be described in further detail with reference to the accompanying drawings, so that those skilled in the art can implement it based on the description.

[0108] It should be understood that terms such as “having,” “comprising,” and “including” as used herein do not exclude the presence or addition of one or more other elements or combinations thereof.

[0109] like Figure 1 As shown, to achieve the above objectives, the present invention adopts the following technical solution, including the following steps:

[0110] S1. Multi-scale structural characterization of shale oil reservoirs;

[0111] High-resolution electron microscopy was used to acquire microscopic images of grain morphology at the micro-nano scale of shale oil reservoir samples, and the grain morphology and arrangement patterns were observed and statistically analyzed. X-ray CT scanning was used to acquire three-dimensional structural images of full-diameter shale oil reservoir cores at the millimeter scale, and image segmentation algorithms were used to extract lamellar thickness, lamellar lateral continuity length, and lamellar distribution patterns. Lamellar characteristics were obtained from the grain microscopic morphology images and the three-dimensional core structural images; these lamellar characteristics include lamellar thickness, lamellar lateral continuity length, lamellar distribution patterns, grain morphology, and arrangement patterns.

[0112] Step S1 includes:

[0113] S101. Select shale oil reservoir rock samples and prepare ultrathin sections with a thickness ≤5μm after argon ion polishing. Observe the samples using a field emission scanning electron microscope or transmission electron microscope at a magnification of 5000-50000 times to obtain images of the grain micromorphology. Analyze the grain micromorphology images and statistically record the following parameters: grain morphology includes the grain size distribution, aspect ratio, and angularity coefficient of different mineral grains; arrangement pattern includes grain distribution and aggregation degree.

[0114] S102. X-ray CT scanning was performed on the full-diameter core of the shale oil reservoir with a scanning resolution ≥50μm to obtain a three-dimensional structural image of the core. A deep learning-based image segmentation algorithm was used to identify the lamellar interface and extract and quantify the following features: geometric parameters include lamellar thickness, lamellar lateral continuous length, and lamellar distribution pattern; among which, lamellar thickness is divided into millimeter-level lamellar (<1cm), centimeter-level lamellar (1-10cm), and decimeter-level blocky (>10cm); lamellar distribution pattern includes the number of lamellars per unit length, i.e., lamellar density, and the distribution uniformity index in different directions.

[0115] S2. Conduct micromechanical experiments;

[0116] Based on the laminar features obtained in S1, typical laminar types were determined; full-diameter core samples from shale oil reservoirs were analyzed to identify typical regions; rock fragments were obtained from typical regions, and mineral dissociation analysis was used to separate and identify mineral components for each rock fragment to determine their types; microscopic morphology images of the grains obtained in S1 for this typical region were selected, and combined with the identified mineral types, a two-dimensional image that can distinguish mineral types was generated; a nanoindenter was used to test a single mineral region to obtain the microscopic elastic modulus of different mineral grains;

[0117] Step S2 includes:

[0118] S201. Based on the described laminar characteristics, determine the typical laminar type; analyze the full-diameter core of the shale oil reservoir to determine the typical region; obtain rock cuttings from the typical region;

[0119] S202. The rock fragments are crushed to a particle size ≤0.1mm and automatically identified using a mineral dissociation analyzer. The elemental composition of each mineral is determined by scanning electron microscopy combined with energy dispersive spectroscopy, thereby identifying the mineral types. The microscopic morphology image of the grains obtained in S1 of the typical area is selected and combined with the mineral types identified above to generate a two-dimensional image that can distinguish the mineral types.

[0120] S203. A nanoindenter was used to test a single mineral region to obtain the micro-elastic modulus of different mineral grains. The micro-elastic modulus Er of the mineral grains was calculated as follows:

[0121] The unloading curve is an exponential function, as shown in formula (1):

[0122] P = B(hh) f ) m (1)

[0123] In equation (1), P represents the load applied during the test, in mN; h represents the indentation depth, in nm; h f The indentation depth at the end of the test is represented in nm; B and m are empirical constants obtained by fitting the test data.

[0124] The contact stiffness S is calculated as shown in formula (2):

[0125]

[0126] In equation (2), h max The maximum indentation depth of the indenter is given by P, measured in nm; the applied load during the test is given by h, measured in mN; h is the test indentation depth, measured in nm. f The indentation depth at the end of the test is represented in nm; B and m are empirical constants obtained by fitting the test data.

[0127] Microscopic elastic modulus Er of mineral grains:

[0128]

[0129] In equation (3), β is generally a constant value, which is related to the pressure head selected in the experiment; A is the contact area, in nm. 2 S represents the contact stiffness.

[0130] S204. Perform statistical analysis on the test data for each type of lamination, calculate the average value and standard deviation of the elastic modulus of each mineral grain, and finally establish a correlation table of microscopic elastic modulus between lamination type, mineral composition and mineral grain.

[0131] S3. Conduct a homogenization analysis to obtain the equivalent elastic modulus of the three-dimensional grain model;

[0132] Based on the two-dimensional image obtained using mineral dissociation analysis technology, pixels are converted into mesh cells to construct a three-dimensional grain model, such as... Figure 2 Then, based on the values ​​of the micro-elastic modulus of the mineral grains of each mineral in step S2, a micro-representative volume element model reflecting the periodic microstructure of the grain model is constructed. Periodic boundary conditions are applied using the finite element method to ensure continuous boundary displacement, thereby obtaining the equivalent elastic modulus of the three-dimensional grain model.

[0133] Step S3 includes:

[0134] S301. The two-dimensional image with a resolution ≥1024×1024 pixels is preprocessed, and the boundaries of different mineral phases are identified by threshold segmentation method;

[0135] S302. Map each pixel to a square grid cell with a side length of 0.1-1μm to construct a three-dimensional grain model;

[0136] S303. Based on the value of the micro-elastic modulus of the mineral grains in step S2, construct a micro-representative volume element model that reflects the periodic microstructure of the grain model. Apply periodic boundary conditions through the finite element method to ensure the continuity of boundary displacement. The applied periodic boundary conditions satisfy equation (4).

[0137] u i (x+Y j )=u i (x) (4)

[0138] In equation (4), u i Let i = x, y, z, representing displacement components, corresponding to displacements in the x, y, and z directions respectively; x is the spatial coordinate, representing the position of any point within the RVE; y... j Let RVE be the period length in the j-direction, where j = x, y, z, as shown below. Figure 3 ;

[0139] S304. Based on the assumption of asymptotic expansion at two scales, the displacement field is decomposed into macroscopic average displacement and microscopic perturbation displacement; based on two-scale expansion, the displacement field is expressed as:

[0140] u ξ (x)=u 0 (x,y)+ξu 1 (x,y)+…… (5)

[0141] In equation (5), y is the microscopic coordinate, y = x / ξ, ξ is the scale parameter, which is the ratio of the macroscopic feature size to the microscopic feature size, ξ << 1; ξ (x) represents the displacement field including scale effects; u 0 (x, y) represents the macroscopic average displacement field; u 1 (x, y) represents the microscopic perturbation displacement field;

[0142] And combined with linear elastic constitutive model:

[0143] σ ij =D ijkl ε kl (6)

[0144] In equation (6), σ ij For the stress tensor components; i,j=x,y,z, representing normal stress or shear stress in different directions; D ijkl ε is the fourth-order elastic stiffness tensor; kl For the strain tensor components, k,l=x,y,z, representing normal strain or shear strain in different directions;

[0145] And the equilibrium equations:

[0146] σ ij,j +f i =0 (7)

[0147] In equation (7), σ ij,j Let f be the divergence of the stress tensor, representing the rate of change of stress in space, and "," denote the partial derivative with respect to coordinates; i Let i be the volume force components, i = x, y, z;

[0148] Establish homogenized control equations and solve for the microstress σ under a given macroscopic strain using the finite element method. ij and strain ε ij Distribution, then by volume average:

[0149]

[0150] In equations (8) and (9), This is the macroscopic average stress tensor; ∫ is the macroscopic average strain tensor; V is the volume of the RVE; V σ ij dV, ∫ V εd ij dV represents the integrals of microscopic stress and strain within the volume of RVE, respectively;

[0151] According to the generalized Hooke's law, equation (10) is:

[0152]

[0153] In equation (12), C is the macroscopic average stress tensor; ijkl This is the equivalent elastic modulus of the three-dimensional grain model; For macroscopic average strain tensor;

[0154] The mechanical response of the RVE under the above boundary conditions was solved using the finite element method, and the results were obtained for each point within the RVE at each strain tensor ε. ij The stress tensor σ under (x) ij (x), and then the equivalent elastic modulus C of the three-dimensional grain model is derived. ijkl .

[0155] S4. Conduct macroscopic physics experiments;

[0156] Shale oil reservoir cores corresponding to multi-scale structural characterization were selected, and cylindrical samples conforming to industry standards were prepared, covering different bedding angles of 0°, 30°, 45°, 60°, and 90°. A servo-controlled triaxial rock testing system was used, with confining pressure gradients set to 5MPa, 10MPa, 15MPa, and 20MPa, and axial strain rate of 0.001mm / s was applied until the sample failed. The failure mode was recorded, and the macroscopic elastic modulus and Poisson's ratio were calculated based on the linear segment of the stress-strain curve.

[0157] In step S4, the macroscopic elastic modulus E is calculated as shown in formula (11):

[0158]

[0159] In equation (11), E is the macroscopic elastic modulus, with units of GPa; σ c (50) represents 50% of the uniaxial compressive strength of the specimen, in MPa; ε c (50) is σ c (50) corresponds to the axial strain;

[0160] The commonly used formula for calculating Poisson's ratio μ is shown in formula (12):

[0161]

[0162] In equation (12), μ is the Poisson's ratio of the specimen, and εc(50) and εt(50) are the axial strain and transverse strain corresponding to σc(50), respectively.

[0163] S5. Determine the mean μ of the striation distribution. X and variance

[0164] The three-dimensional core structure image from step S1 is analyzed to determine the relevant lengths of the lamellar layers. The core three-dimensional structure image is placed horizontally, with the average height of the lamellar layers as Lz and the horizontal extension lengths as Lx and Ly. The dip direction (DD) and dip angle (Dip) are further determined based on the lamellar direction. The mean value (μ) of the lamellar distribution is determined based on the analysis results of the core three-dimensional structure image. X and variance A three-dimensional random field modeling method is used to construct a striation distribution model, such as... Figure 4 ;

[0165] Step S5 includes:

[0166] S501. Based on the three-dimensional core structure image obtained in step S1, Gaussian filtering is used for noise reduction. The processed three-dimensional core structure image is placed horizontally, and the thickness of the laminae is counted along the direction perpendicular to the laminae extension, and the average laminae height Lz is calculated. Along the horizontal directions x and y, the continuous extension length of the laminae is counted, and the horizontal extension lengths Lx and Ly are obtained respectively. Combined with the core attitude data, the dip direction DD and dip angle Dip of the laminae are determined. The dip direction DD ranges from 0° to 360°, and the dip angle Dip ranges from 0° to 90°.

[0167] S502. Perform grayscale processing on the three-dimensional structure image of the core, and calculate the correlation length of the laminae in the x, y, and z directions using the autocorrelation function; for the x direction, the autocorrelation function is used.

[0168]

[0169] In equation (13), ρ(Δx) is the value of the autocorrelation function, which characterizes the strength of the correlation between two points spaced Δx apart in space; Δx is the distance between the two points in space along the x-direction, in mm; c is the lower limit of the autocorrelation function, which is a constant, 0≤c<1, representing the limit of the autocorrelation function when the distance between the two points Δx→∞; Lx is the correlation length in the x-direction. This is the exponentially decaying term, describing how the autocorrelation function changes with the distance Δx between the two points;

[0170] S503. Treat the thickness or spacing of the texture as a discrete random variable X, and its probability distribution law is as follows:

[0171] P{X=xt}=p t (t=1,2,...n) (14)

[0172] In equation (14), n is the sample size, p t xt represents the frequency; xt is the t-th possible value of the discrete random variable X, where t = 1, 2, ..., n.

[0173] The mean is calculated as follows:

[0174]

[0175] In equation (15), μ X (or E(X)) represents the mean of the random variable X; X is a discrete random variable; n is the total number of elements contained in the discrete random variable X within the region D; t is the summation index, used to iterate through the counts from 1 to n; xt is the t-th possible value of the discrete random variable X, t = 1, 2, ..., n; p t Let P{X=xt} be the probability that the random variable X takes the t-th value xt, and satisfy that the sum of all such values ​​is 1.

[0176] Calculate the variance;

[0177] Var(x) = E{[XE(X)]} 2} (16)

[0178] In equation (16), Var(x) represents the variance of the random variable X; X is a random variable; E(X) is the mean of the random variable X;

[0179] Similarly, when the random variable X is discrete, the variance formula can be simplified to:

[0180]

[0181] In equation (17), Let X be the variance of the discrete random variable X; X is a discrete random variable; n is the number of all possible values ​​of the discrete random variable X; t is the summation index, used to iterate through each possible value of X, from 1 to n; xt is the t-th possible value of the discrete random variable X (t = 1, 2, ..., n); pt is the probability that the discrete random variable X takes the t-th value xt, i.e., P{X = xt}; E(X) is the mean of the discrete random variable X.

[0182] S504. A three-dimensional random field modeling method is adopted to construct a striatal distribution model based on exponential or quadratic exponential autocorrelation functions.

[0183] S504a, the relevant length Lx obtained in step S502, and the mean μ obtained in S503. X and variance Using the input parameters, a random field is generated through Fast Fourier Transform;

[0184] S504b discretizes the random field into a standard grid model, adjusts the grid spatial orientation according to the inclination DD and dip angle Dip of the striation, and finally obtains the three-dimensional spatial distribution of striation parameters, realizing the quantitative characterization of striation structure.

[0185] S6. Conduct a secondary homogenization analysis to determine the equivalent elastic modulus of the lamellar distribution model;

[0186] Based on the equivalent elastic modulus of the laminar distribution model constructed in step S5 and the three-dimensional grain model corresponding to different laminar layers obtained in S3, a micro-representative volume element model reflecting the periodic microstructure of the laminar distribution model is constructed, periodic boundary conditions are applied to ensure the continuity of boundary displacement, and the equivalent elastic modulus of the laminar distribution model is determined.

[0187] Step S6 includes:

[0188] S601. Based on the laminar distribution model obtained by three-dimensional random field modeling in step S5, and combined with the equivalent elastic modulus of the three-dimensional grain model obtained in step S3, construct a microscopic representative volume element model that reflects the periodic microstructure of the laminar distribution model. Apply periodic boundary conditions through the finite element method to ensure the continuity of boundary displacement. The applied periodic boundary conditions satisfy the following conditions.

[0189] u′ i (x′+Y′)=u′ i (x′) (18)

[0190] In equation (18), u′ i Let i = x, y, z, corresponding to displacements in the x, y, and z directions, respectively; x′ is the spatial coordinate, representing the position of any point within the RVE; Y′... j Let RVE be the period length in the j-direction, where j = x, y, z; for example... Figure 5 .

[0191] S602. Based on the dual-scale asymptotic expansion assumption, the displacement field is decomposed into macroscopic average displacement and microscopic perturbation displacement; based on dual-scale expansion, the displacement field is expressed as:

[0192] u ξ′ (x')=u 0 (x',y')+ξ'u 1 (x',y')+…… (19) In equation (19), y' is the microscopic coordinate, y' = x' / ξ', ξ' is the scale parameter, which is the ratio of the macroscopic feature size to the microscopic feature size, and ξ' << 1; ξ′ (x') represents the displacement field with scale effects; u0′ (x′,y′) represents the macroscopic average displacement field; u 1′ (x′,y′) represents the microscopic perturbation displacement field;

[0193] And combined with linear elastic constitutive model:

[0194] σ′ ij =D′ ijkl ε′ kl (20)

[0195] In equation (20), σ′ ij For the stress tensor components; i,j=x,y,z, representing normal stress or shear stress in different directions; D′ ijkl ε′ is the fourth-order elastic stiffness tensor; kl For the strain tensor components, k,l=x,y,z, representing normal strain or shear strain in different directions;

[0196] And the equilibrium equations:

[0197] σ′ ij,j +f′ i =0 (21)

[0198] In equation (21), σ′ ij,j f' is the divergence of the stress tensor, representing the rate of change of stress in space; , "," denotes the partial derivative with respect to coordinates; i Let i be the volume force components, i = x, y, z;

[0199] Establish homogenized control equations and solve for the microstress σ′ under a given macroscopic strain using the finite element method. ij and strain ε′ ij Distribution, then by volume average:

[0200]

[0201] In equations (22) and (23), This is the macroscopic average stress tensor; V is the macroscopic average strain tensor; V′ is the volume of the RVE; ∫ V′ σ′ ij dV, ∫ V′ ε′ ij dV represents the integrals of microscopic stress and strain within the volume of RVE, respectively;

[0202] According to the generalized Hooke's law, equation (12) is:

[0203]

[0204] In equation (24), C′ is the macroscopic average stress tensor. ijklThis is the equivalent elastic modulus of the lamellar distribution model; For macroscopic average strain tensor;

[0205] The mechanical response of the RVE under the above boundary conditions was solved using the finite element method, and the results were obtained for each point within the RVE at each strain tensor ε′. ij The stress tensor σ′ under (x) ij (x), from which the equivalent elastic modulus C′ of the lamellar distribution model is derived. ijkl The macroscopic elastic modulus E obtained in step S4 is compared with the macroscopic elastic modulus E. If the error is within 20%, the elastic modulus of shale is considered to be accurately predicted.

[0206] Data Analysis

[0207] Multiple shale oil reservoir samples were selected for multi-scale mechanical property prediction of shale oil reservoirs, and the data are shown in Table 1.

[0208] Table 1

[0209]

[0210]

[0211] As shown in Table 1, the multi-scale mechanical property prediction method for shale oil reservoirs of this invention derives the equivalent elastic modulus C′ of the lamination distribution model. ijkl When compared with the macroscopic elastic modulus E, the error is less than 20%, indicating good accuracy.

[0212] In summary, this invention provides a multi-scale mechanical property prediction method for shale oil reservoirs, which can predict the equivalent elastic modulus of the lamination distribution model across scales, providing a more accurate mechanical basis for the development and engineering practice of shale oil reservoirs.

[0213] Although embodiments of the present invention have been disclosed above, they are not limited to the applications listed in the specification and embodiments. They can be applied to various fields suitable for the present invention. For those skilled in the art, other modifications can be easily made. Therefore, without departing from the general concept defined by the claims and their equivalents, the present invention is not limited to the specific details and illustrations shown and described herein.

Claims

1. A method for predicting the multi-scale mechanical properties of shale oil reservoirs, characterized in that, Includes the following steps: S1. Multi-scale structural characterization of shale oil reservoirs; High-resolution electron microscopy was used to acquire microscopic images of grain morphology at the micro-nano scale of shale oil reservoir samples, and the grain morphology and arrangement patterns were observed and statistically analyzed. X-ray CT scanning was used to acquire three-dimensional structural images of full-diameter shale oil reservoir cores at the millimeter scale, and image segmentation algorithms were used to extract lamellar thickness, lamellar lateral continuity length, and lamellar distribution patterns. Lamellar characteristics were obtained from the grain microscopic morphology images and the three-dimensional core structural images; these lamellar characteristics include lamellar thickness, lamellar lateral continuity length, lamellar distribution patterns, grain morphology, and arrangement patterns. S2. Conduct micromechanical experiments; Based on the laminar features obtained from S1, typical laminar types were determined; full-diameter core samples from shale oil reservoirs were analyzed to identify typical regions; rock fragments were obtained from typical regions, and mineral dissociation analysis was used to separate and identify mineral components for each rock fragment to determine its type; The microstructure images of the grains in the typical region obtained in S1 are selected, and combined with the mineral types identified above, a two-dimensional image that can distinguish mineral types is generated; a nanoindenter is used to test a single mineral region to obtain the microelastic modulus of different mineral grains. S3. Conduct a homogenization analysis to obtain the equivalent elastic modulus of the three-dimensional grain model; Based on the two-dimensional image obtained by mineral dissociation analysis, pixels are converted into mesh units to construct a three-dimensional grain model; then, based on the value of the micro-elastic modulus of the mineral grains of each mineral in step S2, a micro-representative volume element model reflecting the periodic microstructure of the grain model is constructed, and periodic boundary conditions are applied by the finite element method to ensure the continuity of boundary displacement; the equivalent elastic modulus of the three-dimensional grain model is obtained. S4. Conduct macroscopic physics experiments; Shale oil reservoir cores corresponding to multi-scale structural characterization were selected, and cylindrical samples conforming to industry standards were prepared, covering different bedding angles of 0°, 30°, 45°, 60°, and 90°. A servo-controlled triaxial rock testing system was used, with confining pressure gradients set to 5MPa, 10MPa, 15MPa, and 20MPa, and axial strain rate of 0.001mm / s was applied until the sample failed. The failure mode was recorded, and the macroscopic elastic modulus and Poisson's ratio were calculated based on the linear segment of the stress-strain curve. S5. Determine the mean μ of the striation distribution. X and variance The three-dimensional core structure image from step S1 is analyzed to determine the relevant lengths of the lamellar layers. The core three-dimensional structure image is placed horizontally, with the average height of the lamellar layers as Lz and the horizontal extension lengths as Lx and Ly. The dip direction (DD) and dip angle (Dip) are further determined based on the lamellar direction. The mean value (μ) of the lamellar distribution is determined based on the analysis results of the core three-dimensional structure image. X and variance A three-dimensional random field modeling method is used to construct a striation distribution model; S6. Conduct a secondary homogenization analysis to determine the equivalent elastic modulus of the lamellar distribution model; Based on the equivalent elastic modulus of the laminar distribution model constructed in step S5 and the three-dimensional grain model corresponding to different laminar layers obtained in S3, a microscopic representative volume element model reflecting the periodic microstructure of the laminar distribution model is constructed, periodic boundary conditions are applied to ensure the continuity of boundary displacement, and the equivalent elastic modulus of the laminar distribution model is determined.

2. The method for predicting the multi-scale mechanical properties of shale oil reservoirs according to claim 1, characterized in that, Step S1 includes: S101. Select shale oil reservoir rock samples and prepare ultrathin sections with a thickness ≤5μm after argon ion polishing. Observe the samples using a field emission scanning electron microscope or transmission electron microscope at a magnification of 5000-50000 times to obtain images of the grain micromorphology. Analyze the grain micromorphology images and statistically record the following parameters: grain morphology includes the grain size distribution, aspect ratio, and angularity coefficient of different mineral grains; arrangement pattern includes grain distribution and aggregation degree. S102. X-ray CT scanning was performed on the full-diameter core of the shale oil reservoir with a scanning resolution ≥50μm to obtain a three-dimensional structural image of the core. A deep learning-based image segmentation algorithm was used to identify the lamellar interface and extract and quantify the following features: geometric parameters include lamellar thickness, lamellar lateral continuous length, and lamellar distribution pattern; among which, lamellar thickness is divided into millimeter-level lamellar (<1cm), centimeter-level lamellar (1-10cm), and decimeter-level blocky (>10cm); lamellar distribution pattern includes the number of lamellars per unit length, i.e., lamellar density, and the distribution uniformity index in different directions.

3. The method for predicting the multi-scale mechanical properties of shale oil reservoirs according to claim 2, characterized in that, Step S2 includes: S201. Based on the described laminar characteristics, determine the typical laminar type; analyze the full-diameter core of the shale oil reservoir to determine the typical region; obtain rock cuttings from the typical region; S202. The rock fragments are crushed to a particle size ≤0.1mm and automatically identified using a mineral dissociation analyzer. The elemental composition of each mineral is determined by scanning electron microscopy combined with energy dispersive spectroscopy, thereby identifying the mineral types. The microscopic morphology image of the grains obtained in S1 of the typical area is selected and combined with the mineral types identified above to generate a two-dimensional image that can distinguish the mineral types. S203. A nanoindenter was used to test a single mineral region to obtain the micro-elastic modulus of different mineral grains. The micro-elastic modulus Er of the mineral grains was calculated as follows: The unloading curve is an exponential function, as shown in formula (1): P=B(h-h f ) m (1) In equation (1), P represents the load applied during the test, in mN; h represents the indentation depth, in nm; h f The indentation depth at the end of the test is represented in nm; B and m are empirical constants obtained by fitting the test data. The contact stiffness S is calculated as shown in formula (2): In equation (2), h max The maximum indentation depth of the indenter is given by P, measured in nm; the applied load during the test is given by h, measured in mN; h is the test indentation depth, measured in nm. f The indentation depth at the end of the test is represented in nm; B and m are empirical constants obtained by fitting the test data. Microscopic elastic modulus Er of mineral grains: In equation (3), β is generally a constant value, which is related to the pressure head selected in the experiment; A is the contact area, in nm. 2 S represents the contact stiffness. S204. Perform statistical analysis on the test data for each type of lamination, calculate the average value and standard deviation of the elastic modulus of each mineral grain, and finally establish a correlation table of microscopic elastic modulus between lamination type, mineral composition and mineral grain.

4. The method for predicting the multi-scale mechanical properties of shale oil reservoirs according to claim 3, characterized in that, Step S3 includes: S301. The two-dimensional image with a resolution ≥1024×1024 pixels is preprocessed, and the boundaries of different mineral phases are identified by threshold segmentation method; S302. Map each pixel to a square grid cell with a side length of 0.1-1μm to construct a three-dimensional grain model; S303. Based on the value of the micro-elastic modulus of the mineral grains in step S2, construct a micro-representative volume element model that reflects the periodic microstructure of the grain model. Apply periodic boundary conditions through the finite element method to ensure the continuity of boundary displacement. The applied periodic boundary conditions satisfy equation (4). u i (x+Y j )=u i (x) (4) In equation (4), u i Let i = x, y, z, representing displacement components, corresponding to displacements in the x, y, and z directions respectively; x is the spatial coordinate, representing the position of any point within the RVE; y... j Let be the period length of RVE in the j direction, where j = x, y, z; S304. Based on the assumption of asymptotic expansion at two scales, the displacement field is decomposed into macroscopic average displacement and microscopic perturbation displacement; based on two-scale expansion, the displacement field is expressed as: u ξ (x)=u 0 (x,y)+ξu 1 (x,y)+…… (5) In equation (5), y is the microscopic coordinate, y = x / ξ, ξ is the scale parameter, which is the ratio of the macroscopic feature size to the microscopic feature size, ξ << 1; ξ (x) represents the displacement field with scale effects; u 0 (x,y) represents the macroscopic average displacement field; u 1 (x,y) represents the microscopic perturbation displacement field; And combined with linear elastic constitutive model: s ij =D ijkl e kl (6) In equation (6), σ ij For the stress tensor components; i,j=x,y,z, representing normal stress or shear stress in different directions; D ijkl ε is the fourth-order elastic stiffness tensor; kl For the strain tensor components, k,l=x,y,z, representing normal strain or shear strain in different directions; And the equilibrium equations: s ij,j +f i =0 (7) In equation (7), σ ij,j Let f be the divergence of the stress tensor, representing the rate of change of stress in space, and "," denote the partial derivative with respect to coordinates; i Let i be the volume force components, i = x, y, z; Establish homogenized control equations and solve for the microstress σ under a given macroscopic strain using the finite element method. ij and strain ε ij Distribution, then by volume average: In equations (8) and (9), This is the macroscopic average stress tensor; ∫ is the macroscopic average strain tensor; V is the volume of the RVE; V σ ij dV, ∫ V ε ij dV represents the integrals of microscopic stress and strain within the volume of RVE, respectively; According to the generalized Hooke's law, equation (10) is: In equation (12), C is the macroscopic average stress tensor; ijkl This is the equivalent elastic modulus of the three-dimensional grain model; For macroscopic average strain tensor; The mechanical response of the RVE under the above boundary conditions was solved using the finite element method, and the results were obtained for each point within the RVE at each strain tensor ε. ij The stress tensor σ under (x) ij (x), and then the equivalent elastic modulus C of the three-dimensional grain model is derived. ijkl .

5. The method for predicting the multi-scale mechanical properties of shale oil reservoirs according to claim 4, characterized in that, In step S4, the macroscopic elastic modulus E is calculated as shown in formula (11): In equation (11), E is the macroscopic elastic modulus, with units of GPa; σ c (50) represents 50% of the uniaxial compressive strength of the specimen, in MPa; ε c (50) is σ c (50) corresponds to the axial strain; The commonly used formula for calculating Poisson's ratio μ is shown in formula (12): In equation (12), μ is the Poisson's ratio of the specimen, and εc(50) and εt(50) are the axial strain and transverse strain corresponding to σc(50), respectively.

6. The method for predicting the multi-scale mechanical properties of shale oil reservoirs according to claim 5, characterized in that, Step S5 includes: S501. Based on the three-dimensional core structure image obtained in step S1, Gaussian filtering is used for noise reduction. The processed three-dimensional core structure image is placed horizontally, and the thickness of the laminae is counted along the direction perpendicular to the laminae extension, and the average laminae height Lz is calculated. Along the horizontal directions x and y, the continuous extension length of the laminae is counted, and the horizontal extension lengths Lx and Ly are obtained respectively. Combined with the core attitude data, the dip direction DD and dip angle Dip of the laminae are determined. The dip direction DD ranges from 0° to 360°, and the dip angle Dip ranges from 0° to 90°. S502. Perform grayscale processing on the three-dimensional structure image of the core, and calculate the correlation length of the laminae in the x, y, and z directions using the autocorrelation function; for the x direction, the autocorrelation function is used. In equation (13), ρ(Δx) is the value of the autocorrelation function, which characterizes the strength of the correlation between two points spaced Δx apart in space; Δx is the distance between the two points in space along the x-direction, in mm; c is the lower limit of the autocorrelation function, which is a constant, 0≤c<1, representing the limit of the autocorrelation function when the distance between the two points Δx→∞; Lx is the correlation length in the x-direction. This is the exponentially decaying term, describing how the autocorrelation function changes with the distance Δx between the two points; S503. Treat the thickness or spacing of the texture as a discrete random variable X, and its probability distribution law is as follows: P{X=x t }=p t (t=1,2,...n)(14) In equation (14), n is the sample size, p t xt represents the frequency; xt is the t-th possible value of the discrete random variable X, where t = 1, 2, ..., n. The mean is calculated as follows: In equation (15), μ X (or E(X)) represents the mean of the random variable X; X is a discrete random variable; n is the total number of elements contained in the discrete random variable X within the region D; t is the summation index, used to iterate through the counts from 1 to n; xt is the t-th possible value of the discrete random variable X, t = 1, 2, ..., n; p t Let P{X=xt} be the probability that the random variable X takes the t-th value xt, and satisfy that the sum of all such values ​​is 1. Calculate the variance; Var(x) = E{[XE(X)] 2 } (16) In equation (16), Var(x) represents the variance of the random variable X; X is a random variable; E(X) is the mean of the random variable X; Similarly, when the random variable X is discrete, the variance formula can be simplified to: In equation (17), Let X be the variance of the discrete random variable X; X is a discrete random variable; n is the number of all possible values ​​of the discrete random variable X; t is the summation index, used to iterate through each possible value of X, from 1 to n; xt is the t-th possible value of the discrete random variable X (t = 1, 2, ..., n); pt is the probability that the discrete random variable X takes the t-th value xt, i.e., P{X = xt}; E(X) is the mean of the discrete random variable X. S504. A three-dimensional random field modeling method is adopted to construct a striatal distribution model based on exponential or quadratic exponential autocorrelation functions. S504a, the relevant length Lx obtained in step S502, and the mean μ obtained in S503. X and variance Using the input parameters, a random field is generated through Fast Fourier Transform; S504b discretizes the random field into a standard grid model, adjusts the grid spatial orientation according to the inclination DD and dip angle Dip of the striation, and finally obtains the three-dimensional spatial distribution of striation parameters, realizing the quantitative characterization of striation structure.

7. The method for predicting the multi-scale mechanical properties of shale oil reservoirs according to claim 6, characterized in that, Step S6 includes: Step S6 includes: S601. Based on the laminar distribution model obtained by three-dimensional random field modeling in step S5, and combined with the equivalent elastic modulus of the three-dimensional grain model obtained in step S3, construct a microscopic representative volume element model that reflects the periodic microstructure of the laminar distribution model. Apply periodic boundary conditions through the finite element method to ensure the continuity of boundary displacement. The applied periodic boundary conditions satisfy the following conditions. u′ i (x′+Y j ′)=u′ i (x′) (18) In equation (18), u′ i Let i = x, y, z, corresponding to displacements in the x, y, and z directions, respectively; x′ is the spatial coordinate, representing the position of any point within the RVE; Y′... j Let be the period length of RVE in the j direction, where j = x, y, z; S602. Based on the dual-scale asymptotic expansion assumption, the displacement field is decomposed into macroscopic average displacement and microscopic perturbation displacement; based on dual-scale expansion, the displacement field is expressed as: u ξ′ (x′)=u 0 (x′,y')+ξ'u 1 (x',y')+…… (19) In equation (19), y' is the microscopic coordinate, y' = x' / ξ', ξ' is the scale parameter, which is the ratio of the macroscopic feature size to the microscopic feature size, and ξ' << 1; ξ′ (x') represents the displacement field with scale effects; u 0′ (x',y′) represents the macroscopic average displacement field; u 1′ (x',y') represents the microscopic perturbation displacement field; And combined with linear elastic constitutive model: in ij =D' ijkl e' kl (20) In equation (20), σ' ij For the stress tensor components; i,j=x,y,z, representing normal stress or shear stress in different directions; D′ ijkl ε′ is the fourth-order elastic stiffness tensor. kl For the strain tensor components, k,l=x,y,z, representing normal strain or shear strain in different directions; And the equilibrium equations: in ij,j +f′ i =0 (21) In equation (21), σ' ij,j Let f' be the divergence of the stress tensor, representing the rate of change of stress in space, and "," denote the partial derivative with respect to coordinates; i Let i be the volume force components, i = x, y, z; Establish homogenized control equations and solve for the microstress σ' under a given macroscopic strain using the finite element method. ij and strain ε' ij Distribution, then by volume average: In equations (22) and (23), This is the macroscopic average stress tensor; V is the macroscopic average strain tensor; V' is the volume of the RVE; ∫ V' σ' ij dV, ∫ V' ε' ij dV represents the integrals of microscopic stress and strain within the volume of RVE, respectively; According to the generalized Hooke's law, equation (12) is: In equation (24), C' is the macroscopic average stress tensor; ijkl This is the equivalent elastic modulus of the lamellar distribution model; For macroscopic average strain tensor; The mechanical response of the RVE under the above boundary conditions was solved using the finite element method, and the results were obtained for each point within the RVE at each strain tensor ε'. ij The stress tensor σ' under (x) ij (x), from which the equivalent elastic modulus C' of the lamellar distribution model is derived. ijkl The macroscopic elastic modulus E obtained in step S4 is compared with the macroscopic elastic modulus E. If the error is within 20%, the elastic modulus of shale is considered to be accurately predicted.

Images