A grid model generation method for fluid displacement process in cementing engineering

Abstract

A method for generating a grid model for a fluid displacement process in a cementing project, comprising: step one, determining the node parameter values ​​of each grid node along the flow direction of the well depth according to the measured data obtained during the cementing process; step two, generating a grid model for the fluid displacement process in the cementing project according to the node parameter values. This method greatly improves the efficiency of numerical simulation by optimizing the number of grids along the flow direction of the well depth. At the same time, this method makes the grid orthogonality better by exploring the correlation between the centering degree and the circumferential encryption method, and by adjusting the circumferential grid distribution law, thus ensuring the robustness of the simulation. In addition, this method also adjusts the well inclination angle of the grid coordinates, so that the grid model is more in line with the actual situation, and thus the simulation results are more in line with objective reality.

Description

Technical Field

[0001] The invention relates to the technical field of oil and gas exploration and development, and in particular to a grid model generation method for a fluid displacement process in a cementing engineering. Background Art

[0002] Cementing is the last step in the drilling process of oil and gas wells and is the key process connecting drilling and oil production. Improving displacement efficiency is the basic prerequisite for ensuring cementing quality. Among the methods for studying displacement efficiency, cementing displacement numerical simulation has the advantages of not considering experimental conditions and low cost.

[0003] With the development of hardware and software such as computer memory, the numerical simulation capabilities of computers have been enhanced, and numerical simulation has gradually become an important means to study the cementing displacement process. In the numerical simulation of cementing displacement, the flow law of the fluid in the annulus of actual size requires a huge number of grids to describe physical parameters such as wellbore trajectory, centering, and wellbore diameter. Based on this, a physical model can be established to perform numerical simulation of the slurry displacement process in cementing. In this process, the quality of the grid plays a decisive role in the accuracy and efficiency of the simulation. Therefore, based on the actual flow law of the slurry in cementing displacement, establishing a grid model that meets accuracy and efficiency is an important guarantee for simulating the slurry flow process in the entire well section. Summary of the invention

[0004] In order to solve the above problems, the present invention provides a method for generating a grid model for a fluid displacement process in a cementing engineering, the method comprising:

[0005] Step 1: Determine the node parameter values ​​of each grid node along the well depth flow direction according to the measured data obtained during the cementing process;

[0006] Step 2: Generate a grid model for the fluid displacement process of the cementing engineering according to the node parameter values.

[0007] According to an embodiment of the present invention, in step 1, the step of determining the node parameter value of each grid node includes:

[0008] determining the horizontal displacement of the wellbore trajectory according to the measured data;

[0009] According to the measured data and the horizontal displacement of the wellbore, the node parameter values ​​of each grid node along the well depth flow direction are determined by using the difference method.

[0010] According to one embodiment of the present invention, the horizontal displacement of the wellbore trajectory is determined according to the following expression:

[0011] hd(i)=hd(i-1)+[vd(i)-vd(i-1)]×tan[θ(i)+θ(i-1) / 2]

[0012] Among them, hd represents the horizontal displacement of the wellbore trajectory, vd represents the vertical depth, θ represents the well inclination angle, and i and i-1 both represent the measured node numbers.

[0013] According to an embodiment of the present invention, the node parameter includes at least one of the following items:

[0014] Well depth, well inclination, vertical depth, well diameter, centering and horizontal displacement.

[0015] According to one embodiment of the present invention, the step of determining the grid nodes along the well depth flow direction includes:

[0016] According to the length of the starting section of cementing annulus displacement and the grid step length, the number of grids in the starting section is determined;

[0017] Determine the number of grids in the sealing section according to the length of the sealing section to be sealed and the grid step length;

[0018] Determine the number of grids in the development section according to the remaining length of the development section to the wellhead and the grid step length;

[0019] Determine the number of grids in the first transition section according to the grid step length of the starting section, the step length of the sealing section, and the length of the first transition section between the starting section and the sealing section;

[0020] Determining the number of grids of the second transition section according to the step length of the sealing section, the step length of the development section, and the length of the second transition section between the sealing section and the development section;

[0021] Among them, the grid step size gradually increases from the initial stage to the development stage, and the grid density gradually becomes sparse.

[0022] According to one embodiment of the present invention, the step size and the number of grids of the first transition section are determined according to the following expression:

[0023] x 12i =x1+(x2-x1)×i / (m 12 +1)i=1,2,...(m 12 +1)

[0024]

[0025] Among them, x 12i represents the step length between the i-th node and its previous node in the first transition segment, x1 represents the grid step length of the starting segment, x2 represents the step length of the sealing segment, m 12 Indicates the number of grids in the first transition section, L 12 Indicates the length of the first transition segment.

[0026] According to one embodiment of the present invention, the total number of grids along the well depth flow direction is equal to the sum of the number of grids in each well section plus 1.

[0027] According to one embodiment of the present invention, the step 2 includes:

[0028] Generate the initial annulus inner boundary grid and the initial annulus outer boundary grid of each grid node respectively according to the node parameter value of each grid node;

[0029] According to the initial annular inner boundary grid and the initial annular outer boundary grid of each grid node, an initial intermediate layer grid between the initial annular inner boundary grid and the initial annular outer boundary grid is generated, and then an initial grid model of each grid node is obtained;

[0030] The initial annular grid model is adjusted according to the node parameters of each grid node to obtain the required grid model.

[0031] According to one embodiment of the present invention, the initial annulus inner boundary grid is generated according to the following expression:

[0032] α1(j)=180×(j-1) / (j0-1) / (180×π)j=1,2...,j0

[0033] x(i,j,k)=0

[0034] y(i,j,k)=sin(α1)×D in

[0035] z(i,j,k)={cos(α1)×D in -[1-c(i)]×(wd(i)-D in )}×sin(θ(i))

[0036] Where α1 represents the step size of the initial annular inner boundary grid along the circumferential direction, j0 represents the total number of grid nodes in the circumferential direction, i represents the flow direction along the well depth, j represents the circumferential direction, k represents the radial direction, (x, y, z) represents the grid coordinates, and D in represents the outer diameter of the casing, c represents the centering degree of the cross section where the grid node is located, wd represents the well diameter, and θ represents the well inclination angle.

[0037] According to one embodiment of the present invention, the initial annulus outer boundary grid is generated according to the following expression:

[0038] α2=180×((j-1) / (j0-1)) f(c(i)) / (180×π)

[0039] x(i,j,k)=0

[0040] y(i,j,k)=sin(α2)×wd(i)

[0041] z(i,j,k)={cos(α1)×D in -[1-c(i)]×(wd(i)-D in )}×sin(θ(i))

[0042] f(c(i))=0.873+0.127×c(i)

[0043] Where α2 represents the step size of the initial annulus outer boundary grid along the circumferential direction, j0 represents the total number of grid nodes in the circumferential direction, j represents the circumferential direction, k represents the radial direction, (x, y, z) represents the grid coordinates, and D in represents the outer diameter of the casing, c represents the centering of the section where the grid node is located, f(c) represents the centering density factor, and wd represents the well diameter.

[0044] According to one embodiment of the present invention, the initial intermediate layer grid is generated according to the following expression:

[0045]

[0046]

[0047] Where k (k∈[1,N]) represents the radial grid number, β represents the compression factor, (x, y, z) represents the grid coordinates, represents the unit normal vector of the casing outer diameter, represents the unit vector of the line connecting the point on the outer diameter of the casing and the point on the wellbore wall, f represents the parameter that determines the degree of grid curvature, N represents the total number of radial grids, (x1, y1, z1) represents the grid coordinates of the point on the inner diameter of the annulus, (x N ,y N ,z N ) represents the grid coordinates of the points on the outer diameter of the annulus, (x k ,y k ,z k ) represents the grid coordinates of the point numbered k on the radial grid.

[0048] According to one embodiment of the present invention, the initial annulus grid model is adjusted using the well inclination angle, the horizontal displacement and the vertical depth corresponding to each grid node.

[0049] According to one embodiment of the present invention, the initial annulus grid model is adjusted according to the following expression:

[0050] x″(i,j,k)=x′(i,j,k)+hd(i)

[0051] y″(i,j,k)=y′(i,j,k)

[0052] z″(i,j,k)=z′(i,j,k)+vd(i)

[0053] x′(i,j,k)=z(i,j,k)×cos(θ(i))

[0054] y′(i,j,k)=y(i,j,k)

[0055] z′(i,j,k)=z(i,j,k)×sin(θ(i))

[0056] Among them, (x, y, z) and (x″, y″, z″) represent the grid coordinates before and after adjustment, respectively, θ represents the well inclination angle, hd represents the horizontal displacement, vd represents the vertical depth, i represents the flow direction along the well depth, j represents the circumferential direction, and k represents the radial direction.

[0057] The grid model generation method for the fluid displacement process of cementing engineering provided by the present invention greatly improves the efficiency of numerical simulation by optimizing the number of grids along the flow direction of the well depth. At the same time, the method makes the grid orthogonality better by exploring the correlation between the centering degree and the circumferential encryption method, and by adjusting the circumferential grid distribution law, thus ensuring the robustness of the simulation. In addition, the method also adjusts the well inclination angle of the grid coordinates, so that the grid model is more in line with the actual situation, and thus the simulation results are more in line with objective reality.

[0058] Other features and advantages of the present invention will be described in the following description, and partly become apparent from the description, or understood by practicing the present invention. The purpose and other advantages of the present invention can be realized and obtained by the structures particularly pointed out in the description, claims, and drawings. BRIEF DESCRIPTION OF THE DRAWINGS

[0059] In order to more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following briefly introduces the drawings required in the embodiments or the prior art description:

[0060] Figure 1 It is a schematic diagram of the implementation process of a grid model generation method for a fluid displacement process in a cementing engineering according to an embodiment of the present invention;

[0061] Figure 2 is a schematic diagram of an implementation flow of determining node parameter values ​​of each grid node according to an embodiment of the present invention;

[0062] Figure 3 is a schematic diagram of a related coordinate system in a grid model generation process according to an embodiment of the present invention;

[0063] Figure 4 is a schematic diagram of segmentation of a grid along a flow direction of a well depth according to an embodiment of the present invention;

[0064] Figure 5 is a schematic diagram of a process for determining the number of grids along the flow direction of the well depth according to an embodiment of the present invention;

[0065] Figure 6 It is a specific flow chart of generating a grid model for a fluid displacement process of a cementing project according to an embodiment of the present invention;

[0066] Figure 7 is a schematic diagram of grid coordinate transformation when the well inclination angle is θ according to one embodiment of the present invention;

[0067] Figure 8 is a schematic diagram of a grid node of a well according to an embodiment of the present invention;

[0068] Fig. 9 is a schematic diagram of an annulus inner boundary grid of a well according to one embodiment of the present invention;

[0069] Fig.10 A schematic diagram of an annulus outer boundary grid of a well according to an embodiment of the present invention;

[0070] Fig.11 is a schematic diagram of an intermediate layer grid of a well according to an embodiment of the present invention;

[0071] Figures 12 to 15 is a schematic diagram of adjusting grid coordinates according to an embodiment of the present invention. DETAILED DESCRIPTION

[0072] The following will describe the implementation methods of the present invention in detail with reference to the accompanying drawings and embodiments, so that the implementation process of how the present invention applies technical means to solve technical problems and achieve technical effects can be fully understood and implemented accordingly. It should be noted that as long as there is no conflict, the various embodiments of the present invention and the various features in the embodiments can be combined with each other, and the technical solutions formed are all within the protection scope of the present invention.

[0073] Meanwhile, in the following description, many specific details are set forth for the purpose of explanation to provide a thorough understanding of the embodiments of the present invention. However, it is obvious to those skilled in the art that the present invention can be implemented without the specific details or the specific manner described herein.

[0074] In addition, the steps shown in the flowcharts of the accompanying drawings can be executed in a computer system such as a set of computer-executable instructions, and although a logical order is shown in the flowcharts, in some cases, the steps shown or described can be performed in an order different from that shown here.

[0075] In the existing technology, cementing displacement numerical simulation grid model is usually established through commercial software to study the cementing displacement process. However, there are many disadvantages in using commercial software for grid generation. For example:

[0076] 1) Commercial software cannot use the measured data of the well, and the generated grid model is too simplified

[0077] For the cementing data on site, the parameters such as well diameter, well inclination, and centering are not monotonically changing. Because commercial software cannot use the measured data of the well, the grid it generates can only simplify the model to the situation where the parameters change monotonically, so it is impossible to describe the annular state of the actual parameters.

[0078] 2) Commercial software cannot reasonably encrypt according to the flow rules

[0079] Although commercial software can improve calculation accuracy by encrypting the grid, it is unable to carry out targeted encryption along the flow direction of the well depth and adjust the circumferential orthogonality results according to the specific flow law of cementing displacement.

[0080] In view of the above problems existing in the prior art, the present invention provides a new grid model generation method for the fluid displacement process of cementing engineering. The method can generate a grid model describing the measured data, and can carry out targeted partition encryption along the well depth flow direction and circumferential orthogonal node adjustment according to the flow law.

[0081] Figure 1 The schematic diagram of the implementation flow of the grid model generation method for the fluid displacement process of cementing engineering provided by this embodiment is shown.

[0082] like Figure 1 As shown, the grid model generation method for the cementing engineering fluid displacement process provided in this embodiment first determines the node parameter values ​​of each grid node along the well depth flow direction according to the measured data obtained during the cementing process in step S101.

[0083] Specifically, in this embodiment, the measured data obtained by the method in step S101 preferably include data such as well depth, well inclination angle, vertical depth, well diameter, centering degree and casing outer diameter. Of course, in other embodiments of the present invention, according to actual needs, the measured data obtained by the method in step S101 may include only one or several of the listed items, or may include other reasonable parameters, and the present invention is not limited to this.

[0084] Figure 2 A schematic diagram of the implementation flow of determining the node parameter values ​​of each grid node in this embodiment is shown.

[0085] like Figure 2 As shown, in this embodiment, when determining the node parameter values ​​of each grid node, it is preferred to first determine the horizontal displacement of the wellbore trajectory in step S201 based on the measured data obtained in step S101, and then determine the node parameter values ​​of each grid node along the well depth flow direction by using the difference method based on the above-mentioned measured data and the horizontal displacement of the wellbore in step S202.

[0086] In this embodiment, the method preferably determines the horizontal displacement of the wellbore trajectory according to the vertical depth and the well inclination angle in the measured data in step S201. For example, the method can determine the horizontal displacement of the wellbore trajectory according to the following expression:

[0087] hd(i)=hd(i-1)+[vd(i)-vd(i-1)]×tan[θ(i)+θ(i-1) / 2] (1)

[0088] Among them, hd represents the horizontal displacement of the wellbore trajectory, vd represents the vertical depth, θ represents the well inclination, and i and i-1 represent the measured node numbers. In this way, hd(i) and hd(i-1) represent the horizontal displacement of the i-th and i-1-th nodes, vd(i) and vd(i-1) represent the vertical depth of the i-th and i-1-th nodes, and θ(i) and θ(i-1) represent the well inclination of the i-th and i-1-th nodes.

[0089] Of course, in other embodiments of the present invention, the method may also adopt other reasonable methods to determine the horizontal displacement of the wellbore trajectory, and the present invention does not specifically limit this.

[0090] In step S202, in this embodiment, the method preferably determines the node distribution of the grid along the well depth flow direction according to the lengths of the starting section, the sealing section and the development section.

[0091] like Figure 3As shown in the figure, in the process of grid generation, there are two sets of coordinate systems, one is the Cartesian coordinate system, and the other is the grid coordinate system. In the Cartesian coordinate system, the x direction is the horizontal direction, the z direction is the vertical direction, the y direction is the direction perpendicular to the xz plane, and the coordinates of the center of the wellhead are (0, 0, 0).

[0092] The grid coordinate system is a cylindrical coordinate system. In the cylindrical coordinate system, the i direction is the flow direction along the well depth, and its zero point is the well mouth position; the j direction is the circumferential direction, and its zero point is the zero point position of the circumference as the dial, and its value range is 0° to 360°; the k direction is the radial direction, and its zero point is the axis position. The starting position of the grid in this direction is the innermost position of the model, that is, the outer diameter position of the casing. i0 is the number of grids along the well depth flow direction, j0 is the number of circumferential grids, and k0 is the number of radial grids.

[0093] like Figure 4 In the grid, the section at the beginning of the displacement flow (i.e. the starting section) is encrypted to capture the interface at the beginning; followed by the sealing section, whose grid is sparser than that of the starting section, and this section is used to capture the interface in the middle of the displacement; then the development section is located near the wellhead, and the grid in this section is particularly sparse, which improves the calculation efficiency while ensuring the calculation accuracy.

[0094] Among them, the lengths of the starting section, the sealing section, and the development section are preferably calculated from a large number of numerical simulation experiments, which can not only ensure the calculation accuracy and accurately describe the displacement interface morphology, but also ensure the calculation efficiency.

[0095] Figure 5 A schematic diagram of the implementation process of determining the number of grids in the starting section of cementing annulus displacement in this embodiment is shown.

[0096] like Figure 5 As shown, in this embodiment, the method preferably determines the number of grids m1 of the starting segment according to the length L1 of the starting segment and the grid step x1 in step S501, and determines the number of grids m2 of the sealing segment according to the length L2 of the sealing segment and the grid step x2 in step S502, and determines the number of grids m3 of the development segment according to the length L3 of the development segment and the grid step x3 in step S503.

[0097] For example, the number of grids in the starting segment m1 can be calculated according to the following expression:

[0098]

[0099] The grid step length x2 of the sealing section can be calculated according to the following expression:

[0100]

[0101] The development segment grid step x3 can be calculated according to the following expression:

[0102]

[0103] In step S504, the method preferably determines the first transition section step length and its grid number m according to the obtained starting section grid step length x1, the sealing section step length x2 and the length of the first transition section between the starting section and the sealing section. 12 .

[0104] For example, in this embodiment, the method may determine the first transition step length and its length according to the following expression:

[0105] x 12i =x1+(x2-x1)×i / (m 12 +1)i=1,2,...(m 12 +1) (5)

[0106]

[0107] Among them, x 12i represents the step length between the i-th node and its previous node in the first transition segment, m 12 Indicates the number of grids in the first transition section, L 12 Indicates the length of the first transition segment.

[0108] Similarly, the method will also determine the second transition section step length and its grid number m in step S505 according to the sealing section step length x2, the development section step length x3 and the length of the second transition section between the sealing section and the development section. 23 .

[0109] For example, in this embodiment, the method can determine the second transition segment step length and its length according to the following expression:

[0110] x 23i =x2+(x3-x2)×i / (m 23 +1)i=1,2,...(m 23 +1) (7)

[0111]

[0112] Among them, x 23i represents the step length between the ith node and its previous node in the second transition segment, m 23 Indicates the number of grids in the second transition section, L 23 Indicates the length of the second transition segment.

[0113] It should be pointed out that due to the possible existence of remainders and other problems, in order to ensure that the grids along the well depth flow direction can cover the entire length of the well, in the design of the number of grids along the well depth flow direction, the total number of grids needs to be added by 1 to deal with the possible problem of remainders that cannot be divided evenly. Therefore, the total number of grids along the well depth flow direction is equal to the sum of the number of grids in each section plus 1.

[0114] Since the actually measured data do not necessarily correspond to each grid node, that is, the number of grid nodes is greater than the number of nodes of the measured data, this method preferably determines the node parameter values ​​of each grid node along the well depth flow direction using the difference method based on the measured data and the horizontal displacement of the wellbore.

[0115] For example, the node parameter values ​​of a grid node located between node i1 and node i2 along the well depth flow direction can be calculated according to the following expressions:

[0116] d(i)=(i-i1) / (i2-i1)×(d(i2)-d(i1)) (9)

[0117] θ(i)=(i-i1) / (i2-i1)×(θ(i2)-θ(i1)) (10)

[0118] vd(i)=(i-i1) / (i2-i1)×(vd(i2)-vd(i1)) (11)

[0119] wd(i)=(i-i1) / (i2-i1)×(wd(i2)-wd(i1)) (12)

[0120] c(i)=(i-i1) / (i2-i1)×(c(i2)-c(i1)) (13)

[0121] hd(i)=(i-i1) / (i2-i1)×(hd(i2)-hd(i1)) (14)

[0122] Of course, in other embodiments of the present invention, the method also adopts other reasonable methods to determine the node parameter values ​​of each grid node along the well depth flow direction, and the present invention does not specifically limit this.

[0123] Again Figure 1 As shown, in this embodiment, after obtaining the node parameter values ​​of each grid node along the flow direction of the well depth, the method preferably generates a grid model for the cementing engineering fluid displacement process according to the above node parameter values ​​in step S102.

[0124] Figure 6A specific schematic diagram of the process of generating a grid model for the fluid displacement process of cementing engineering in this embodiment is shown.

[0125] like Figure 6 As shown, in this embodiment, the method preferably first generates an initial annulus inner boundary grid and an initial annulus outer boundary grid of each grid node according to the node parameter values ​​of each grid node in step S601.

[0126] The inner boundary of the annulus is the outer diameter of the casing. Its initial position is set at the wellhead, that is, the position where its horizontal displacement and vertical depth are 0; the initial well inclination of the preset grid is 90°, and it is adjusted to the actual angle in the subsequent well inclination adjustment; the grid section is perpendicular to the flow direction along the well depth and is a circle. In order to improve the calculation efficiency, the 180-degree grid is used instead of the circular 360° grid by using its axisymmetric properties. When generating the inner boundary grid, the coordinates are adjusted considering the influence of the casing centering.

[0127] Based on the above principles, in this embodiment, the method can preferably generate the initial annulus inner boundary grid according to the following expression in step S601:

[0128] α1(j)=180×(j-1) / (j0-1) / (180×π)j=1,2...,j0 (15)

[0129] x(i,j,k)=0 (16)

[0130] y(i,j,k)=sin(α1)×D in (17)

[0131] z(i,j,k)={cos(α1)×D in -[1-c(i)]×(wd(i)-D in )}×sin(θ(i)) (18)

[0132] Where α represents the step size of the initial annular inner boundary grid along the circumferential direction, j0 represents the total number of grid nodes in the circumferential direction, i represents the flow direction along the well depth, j represents the circumferential direction, k represents the radial direction, (xy,z) represents the grid coordinates, and D in represents the outer diameter of the casing, c represents the centering degree of the cross section where the grid node is located, wd represents the well diameter, and θ represents the well inclination angle.

[0133] Subsequently, the method preferably generates an initial intermediate layer grid between the initial annular inner boundary grid and the initial annular outer boundary grid according to the initial annular inner boundary grid and the initial annular outer boundary grid of each grid node in step S602, and then obtains the initial grid model of each grid node.

[0134] The outer boundary of the annulus is the wellbore, and its initial position is set at the wellhead, that is, the position where its horizontal displacement and vertical depth are 0; the initial well inclination of the preset grid is 90°, and it is adjusted to the actual angle in the subsequent well inclination adjustment; the grid section is perpendicular to the flow direction along the well depth and is a circle. In order to improve the calculation efficiency, the axial symmetry property is used to use a 180-degree circular grid instead of a 360° circular grid; according to the nodes and eccentricity on the inner boundary, the grid distribution on the outer boundary is adjusted so that the unit vector connecting the inner boundary points and the outer boundary points is as close as possible to the direction of the inner boundary normal vector, thereby improving the grid quality.

[0135] Based on the above principles, in this embodiment, the method can preferably generate the initial annulus outer boundary grid according to the following expression in step S602:

[0136] α2=180×((j-1) / (j0-1)) f(c(i)) / (180×π) (19)

[0137] x(i,j,k)=0 (20)

[0138] y(i,j,k)=sin(α2)×wd(i) (21)

[0139] z(i,j,k)={cos(α1)×D in -[1-c(i)]×(wd(i)-D in )}×sin(θ(i)) (22)

[0140] f(c(i))=0.873+0.127×c(i) (23)

[0141] Among them, α2 represents the step size of the initial annulus outer boundary grid along the circumferential direction, j0 represents the total number of grid nodes in the circumferential direction, j represents the circumferential direction, k represents the radial direction, (x, y, z) represents the grid coordinates, c represents the centering of the section where the grid node is located, f(c) represents the centering density factor, and wd represents the wellbore diameter.

[0142] The method preferably generates an initial intermediate layer grid between the initial annular inner boundary grid and the initial annular outer boundary grid according to the initial annular inner boundary grid and the initial annular outer boundary grid of each grid node, and then obtains the initial grid model of each grid node;

[0143] In order to accurately simulate the flow of the boundary layer near the outer diameter of the inner boundary casing, a dense grid is required near the wall. The present invention uses an exponential compression function to encrypt the wall near the outer diameter of the casing. Specifically, the method can generate an initial intermediate layer grid according to the following expression in step S602:

[0144]

[0145] Where k (k∈[1,N]) represents the radial grid number. When k is 1, it corresponds to the outer diameter of the casing (i.e., the inner diameter of the annulus). When k is N, it corresponds to the wellbore diameter (i.e., the outer diameter of the annulus). β represents the compression factor, (x, y, z) represents the grid coordinates, represents the unit normal vector of the casing outer diameter, represents the unit vector of the line connecting the point on the outer diameter of the casing and the point on the wellbore wall, f represents the parameter that determines the degree of grid curvature, N represents the total number of radial grids, (x1, y1, z1) represents the grid coordinates of the point on the inner diameter of the annulus, (x N ,y N ,z N ) represents the grid coordinates of the points on the outer diameter of the annulus, (x k ,y k ,z k ) represents the grid coordinates of the point numbered k on the radial grid.

[0146] In this embodiment, the casing outer diameter unit normal vector It can be determined according to the following expression:

[0147]

[0148] The unit vector of the line connecting the point on the outer diameter of the casing and the corresponding point on the well wall is It can be determined according to the following expression:

[0149]

[0150] Among them, R bN It represents the length of the line connecting a point on the outer diameter of the casing and the corresponding point on the well wall.

[0151] In this embodiment, the parameter f that determines the degree of mesh curvature preferably satisfies the following conditions:

[0152] 0≤f(k)≤1,k=1,2,…,N (29)

[0153] f(k)=1,k=1 (30)

[0154] Obviously, if f(k) = 1-a(k), then when k = 1, f(1) = 1, and the grid is always perpendicular to the object surface; and when k = N, f(N) = 1, and the grid is always perpendicular to the object surface. The direction is consistent.

[0155] Of course, in other embodiments of the present invention, the method may also adopt other reasonable methods to determine the initial grid model according to actual needs, and the present invention does not specifically limit this.

[0156] like Figure 6 As shown, after the initial network model is obtained, the method preferably adjusts the obtained initial annular grid model according to the node parameters of each grid node in step S603, so as to obtain the final required grid model.

[0157] Specifically, in this embodiment, the method preferably uses the well inclination angle, the horizontal displacement and the vertical depth corresponding to each grid node to adjust the initial annulus grid model in step S603.

[0158] The grid coordinate pair generated in step S602 corresponds to the coordinates of the initial position of the grid. In this case, the well inclination angle is 90°, so at this time there is:

[0159] x(i,j,k)=0 (31)

[0160] like Figure 7 As shown in the figure, when the well inclination angle is θ, the corresponding grid coordinates are transformed into:

[0161] x′(i,j,k)=z(i,j,k)×cos(θ(i)) (32)

[0162] y′(i,j,k)=y(i,j,k) (33)

[0163] z′(i,j,k)=z(i,j,k)×sin(θ(i)) (34)

[0164] Based on the grid coordinates calculated by expressions (32) to (34), the method will also superimpose the horizontal displacement and vertical depth corresponding to the grid origin on the above grid coordinates in step S603, reflect the influence of the spatial position of each node grid in the grid three-dimensional coordinates, and then generate an annular space grid model in a unified coordinate system.

[0165] Specifically, in this embodiment, the method can preferably obtain the final grid model according to the following expression:

[0166] x″(i,j,k)=x′(i,j,k)+hd(i) (35)

[0167] y″(i,j,k)=y′(i,j,k) (36)

[0168] z″(i,j,k)=z′(i,j,k)+vd(i) (37)

[0169] Among them, (x, y, z) and (x″, y″, z″) represent the grid coordinates before and after adjustment, respectively, θ represents the well inclination angle, hd represents the horizontal displacement, and vd represents the vertical depth.

[0170] Of course, in other embodiments of the present invention, the method may also adopt other reasonable methods according to actual needs to adjust the initial annular grid model according to the node parameters of each grid node in step S603 to obtain the required grid model.

[0171] In this embodiment, the flow area of ​​the cementing displacement process is an annular space composed of the outer diameter of the casing as the inner layer and the wellbore diameter as the outer layer. The well inclination angle refers to the angle between the central axis of a point in the oil well and the plumb line of the earth, ranging from 0° to 180°. The well inclination angle is used to indicate the inclination of the wellbore trajectory. The centering degree describes the degree to which the casing deviates downward in the wellbore due to the effect of gravity.

[0172] In order to demonstrate the reliability, applicability and advantages of the grid model generation method provided in this embodiment, a certain well is taken as an example for further explanation below.

[0173] Take the logging data of a well as an example. The following table shows the well depth, well inclination, vertical depth, well diameter, and centering data obtained during the cementing process.

[0174] Table 1

[0175]

[0176] Using the well depth, vertical depth and well inclination angle, this method can calculate the horizontal displacement of the wellbore trajectory according to expression (1), thereby obtaining the data shown in Table 2:

[0177] Table 2

[0178]

[0179]

[0180] The length of the sealing section corresponding to the well is 600m. According to the length of the sealing section, this method can determine the number of grid nodes along the flow direction of the well depth.

[0181] Figure 8 Indicates that one dot is displayed every 50 i dots. Figure 8 It can be seen that the number of grids in the starting section exceeds 500, which is used to simulate the initial formation process of the displacement interface, and its length is less than 20m; the number of grids in the sealing section exceeds 1500, which is used to quickly simulate the development process of the displacement interface in the sealing section, and its length is 600m; the development section has the least number of grid nodes, less than 300, which is used to quickly simulate the process of the displacement interface moving to the wellhead, and the simulation length is approximately 700m.

[0182] Based on the well depth, well inclination, vertical depth, well diameter, centering, casing outer diameter, and horizontal displacement data, this method can use the interpolation method to calculate the above data at the nodes along the well depth flow direction. Specifically, the data shown in Table 3 can be obtained:

[0183] Table 3

[0184]

[0185] This method can then construct the cross-sectional mesh of each node according to the parameters of each node, thereby further generating an annulus mesh model. Fig. 9 The generated inner boundary grid of each node is shown (i.e., inner layer grid), and the grid coordinates are adjusted according to the centering degree. Fig.10 As shown, this method can define the outer boundary grid of the annulus (i.e., the outer grid) according to the wellbore data, and adjust the circumferential distribution of the outer boundary grid of the annulus according to the inner boundary grid of the annulus, thereby improving the orthogonality of the grid.

[0186] Then, if Fig.11 As shown, this method generates an intermediate layer grid between the annulus inner boundary grid and the annulus outer convenient grid corresponding to each node according to the annulus inner boundary grid and the annulus outer boundary grid of each node.

[0187] This method adjusts the grid coordinates according to the well inclination angle of each node. The specific results are as follows Fig.12 and Fig.13 As shown. Fig.12 and Fig.13 It can be seen that in the yz function of the well inclination adjustment, the z coordinate is -6 to 6, and the y coordinate is -250 to 250. This is because the influence of the well inclination causes the coordinates to have projections in both the yz plane and the xy plane.

[0188] Finally, if Fig.14 and Fig.15 As shown, this method will further adjust the grid coordinates according to the spatial coordinates of each node to generate the final required grid model.

[0189] From the above description, it can be seen that the grid model generation method for the fluid displacement process of cementing engineering provided by the present invention greatly improves the efficiency of numerical simulation by optimizing the number of grids along the flow direction of the well depth. At the same time, this method makes the grid orthogonality better by exploring the correlation between the centering degree and the circumferential encryption method, and by adjusting the circumferential grid distribution law, thus ensuring the robustness of the simulation. In addition, the method also adjusts the well inclination angle of the grid coordinates, so that the grid model is more in line with the actual situation, and thus the simulation results are more in line with objective reality.

[0190] It should be understood that the embodiments disclosed in the present invention are not limited to the specific structures or processing steps disclosed herein, but should be extended to equivalent substitutions of these features understood by ordinary technicians in the relevant field. It should also be understood that the terms used herein are only used for the purpose of describing specific embodiments and are not meant to be limiting.

[0191] The "one embodiment" or "embodiment" mentioned in the specification means that the specific features, structures or characteristics described in conjunction with the embodiment are included in at least one embodiment of the present invention. Therefore, the phrases "one embodiment" or "embodiment" appearing in various places throughout the specification do not necessarily refer to the same embodiment.

[0192] Although the above examples are used to illustrate the principles of the present invention in one or more applications, it is obvious to those skilled in the art that various modifications can be made in form, usage and implementation details without inventive effort without departing from the principles and ideas of the present invention. Therefore, the present invention is defined by the appended claims.

Claims

1. A method for generating a grid model for a fluid displacement process in a cementing engineering, characterized in that: The method comprises: Step 1: Determine the node parameter values ​​of each grid node along the well depth flow direction according to the measured data obtained during the cementing process, wherein the horizontal displacement of the wellbore trajectory is first determined according to the measured data, and then the node parameter values ​​of each grid node along the well depth flow direction are determined by using the difference method according to the measured data and the wellbore horizontal displacement; Step 2: Generate a grid model for the fluid displacement process of the cementing engineering according to the node parameter values, wherein the initial annulus inner boundary grid and the initial annulus outer boundary grid of each grid node are first generated according to the node parameter values ​​of each grid node, and then the initial intermediate layer grid between the initial annulus inner boundary grid and the initial annulus outer boundary grid is generated according to the initial annulus inner boundary grid and the initial annulus outer boundary grid of each grid node, and then the initial grid model of each grid node is obtained, and the initial annulus outer boundary grid model is further adjusted according to the node parameters of each grid node to obtain the required grid model.

2. The method according to claim 1, characterized in that The horizontal displacement of the wellbore trajectory is determined according to the following expression: hd(i)=hd(i-1)+[vd(i)-vd(i-1)]×tan[θ(i)+θ(i-1) / 2] Among them, hd represents the horizontal displacement of the wellbore trajectory, vd represents the vertical depth, θ represents the well inclination angle, and i and i-1 both represent the measured node numbers.

3. The method according to claim 1 or 2, characterized in that The node parameters include at least one of the following items: Well depth, well inclination, vertical depth, well diameter, centering and horizontal displacement.

4. The method according to claim 3, characterized in that The steps of determining the grid nodes along the flow direction of the well depth include: According to the length of the starting section of cementing annulus displacement and the grid step length, the number of grids in the starting section is determined; Determine the number of grids in the sealing section according to the length of the sealing section to be sealed and the grid step length; Determine the number of grids in the development section according to the remaining length of the development section to the wellhead and the grid step length; Determine the number of grids in the first transition section according to the grid step length of the starting section, the step length of the sealing section, and the length of the first transition section between the starting section and the sealing section; Determining the number of grids of the second transition section according to the step length of the sealing section, the step length of the development section, and the length of the second transition section between the sealing section and the development section; Among them, the grid step size gradually increases from the initial stage to the development stage, and the grid density gradually becomes sparse.

5. The method according to claim 4, characterized in that The step size and number of grids of the first transition segment are determined according to the following expressions: x 12i =x1+(x2-x1)×i / (m 12 +1)i=1,2,...(m 12 +1) Among them, x 12i represents the step length between the i-th node and its previous node in the first transition segment, x1 represents the grid step length of the starting segment, x2 represents the step length of the sealing segment, m 12 Indicates the number of grids in the first transition section, L 12 Indicates the length of the first transition segment.

6. The method according to claim 4 or 5, characterized in that The total number of grids along the well depth flow direction is equal to the sum of the number of grids in each well section plus 1.

7. The method according to claim 6, characterized in that The initial annulus inner boundary grid is generated according to the following expression: α1(j)=180×(j-1) / (j0-1) / (180×π)j=1,2...,j0 x(i,j,k)=0 y(i,j,k)=sin(α1)×D in z(i,j,k)={cos(α1)×D in -[1-c(i)]×(wd(i)-D in )}×sin(θ(i)) Where α1 represents the step size of the initial annular inner boundary grid along the circumferential direction, j0 represents the total number of grid nodes in the circumferential direction, i represents the flow direction along the well depth, j represents the circumferential direction, k represents the radial direction, (x, y, z) represents the grid coordinates, and D in represents the outer diameter of the casing, c represents the centering degree of the cross section where the grid node is located, wd represents the well diameter, and θ represents the well inclination angle.

8. The method according to claim 7, characterized in that The initial annulus outer boundary grid is generated according to the following expression: α2=180×((j-1) / (j0-1)) f(c(i)) / (180×π) x(i,j,k)=0 y(i,j,k)=sin(α2)×wd(i) z(i,j,k)={cos(α1)×D in -[1-c(i)]×(wd(i)-D in )}×sin(θ(i)) f(c(i))=0.873+0.127×c(i) Where α2 represents the step size of the initial annulus outer boundary grid along the circumferential direction, j0 represents the total number of grid nodes in the circumferential direction, i represents the flow direction along the well depth, j represents the circumferential direction, k represents the radial direction, (x, y, z) represents the grid coordinates, and D in represents the outer diameter of the casing, c represents the centering of the section where the grid node is located, f(c) represents the centering density factor, wd represents the well diameter, θ represents the well inclination angle, and α1 represents the step size of the initial annulus inner boundary grid along the circumferential direction.

9. The method according to claim 7 or 8, characterized in that The initial intermediate layer grid is generated according to the following expression: Where k (k∈[1,N]) represents the radial grid number, β represents the compression factor, (x, y, z) represents the grid coordinates, represents the unit normal vector of the casing outer diameter, represents the unit vector of the line connecting the point on the outer diameter of the casing and the point on the wellbore wall, f represents the parameter that determines the degree of grid curvature, N represents the total number of radial grids, (x1, y1, z1) represents the grid coordinates of the point on the inner diameter of the annulus, (x N ,y N ,z N ) represents the grid coordinates of the points on the outer diameter of the annulus, (x k ,y k ,z k ) represents the grid coordinates of the point numbered k on the radial grid.

10. The method according to claim 9, characterized in that The initial annulus grid model is adjusted using the well inclination angle, the horizontal displacement and the vertical depth corresponding to each grid node.

11. The method according to claim 10, characterized in that The initial annulus grid model is adjusted according to the following expression: x″(i,j,k)=x′(i,j,k)+hd(i) y″(i,j,k)=y′(i,j,k) z″(i,j,k)=z′(i,j,k)+vd(i) x′(i,j,k)=z(i,j,k)×cos(θ(i)) y′(i,j,k)=y(i,j,k) z′(i,j,k)=z(i,j,k)×sin(θ(i)) Among them, (x, y, z) and (x″, y″, z″) represent the grid coordinates before and after adjustment, respectively, θ represents the well inclination angle, hd represents the horizontal displacement, vd represents the vertical depth, i represents the flow direction along the well depth, j represents the circumferential direction, and k represents the radial direction.

Images